Rosetta Code

Simulated annealing: Difference between revisions

Page Discussion

Simulated annealing (view source)

Revision as of 05:59, 1 May 2024

43,337 bytes added ,  18 days ago
→‎{{header|jq}}
 
(12 intermediate revisions by 5 users not shown)
Line 29:
We want to apply SA to the travelling salesman problem. There are 100 cities, numbered 0 to 99, located on a plane, at integer coordinates i,j : 0 <= i,j < 10 . The city at (i,j) has number 10*i + j. The cities are '''all''' connected : the graph is complete : you can go from one city to any other city in one step.
 
The salesman wants to start from city 0, visit all cities, each one time, and go back to city 0. The travel cost between two cities is the euclidian distance between there citiesthem. The total travel cost is the total path length.
 
A path '''s''' is a sequence (0 a b ...z 0) where (a b ..z) is a permutation of the numbers (1 2 .. 99). The path length = E(s) is the sum d(0,a) + d(a,b) + ... + d(z,0) , where d(u,v) is the distance between two cities. Naturally, we want to minimize E(s).
Line 88:
 
Tune the parameters kT, kmax, or use different temperature() and/or neighbour() functions to demonstrate a quicker convergence, or a better optimum.
 
=={{header|Ada}}==
{{trans|C}}
{{trans|Scheme}}
{{works with|GNAT|Community 2021}}
 
 
This implementation is adapted from the C, which was adapted from the Scheme. It uses fixed-point numbers for no better reason than to demonstrate that Ada has fixed-point numbers support built in.
 
 
<syntaxhighlight lang="ada">----------------------------------------------------------------------
--
-- The Rosetta Code simulated annealing task in Ada.
--
-- This implementation demonstrates that Ada has fixed-point numbers
-- support built in. Otherwise there is no particular reason I used
-- fixed-point instead of floating-point numbers.
--
-- (Actually, for the square root and exponential, I cheat and use the
-- floating-point functions.)
--
----------------------------------------------------------------------
 
with Ada.Numerics.Discrete_Random;
with Ada.Numerics.Elementary_Functions;
with Ada.Text_IO; use Ada.Text_IO;
 
procedure simanneal
is
 
Bigint : constant := 1_000_000_000;
Bigfpt : constant := 1_000_000_000.0;
 
-- Fixed point numbers.
type Fixed_Point is delta 0.000_01 range 0.0 .. Bigfpt;
 
-- Integers.
subtype Zero_or_One is Integer range 0 .. 1;
subtype Coordinate is Integer range 0 .. 9;
subtype City_Location is Integer range 0 .. 99;
subtype Nonzero_City_Location is City_Location range 1 .. 99;
subtype Path_Index is City_Location;
subtype Nonzero_Path_Index is Nonzero_City_Location;
 
-- Arrays.
type Path_Vector is array (Path_Index) of City_Location;
type Neighborhood_Array is array (1 .. 8) of City_Location;
 
-- Random numbers.
subtype Random_Number is Integer range 0 .. Bigint - 1;
package Random_Numbers is new Ada.Numerics.Discrete_Random
(Random_Number);
use Random_Numbers;
 
gen : Generator;
 
function Randnum
return Fixed_Point
is
begin
return (Fixed_Point (Random (gen)) / Fixed_Point (Bigfpt));
end Randnum;
 
function Random_Natural
(imin : Natural;
imax : Natural)
return Natural
is
begin
-- There may be a tiny bias in the result, due to imax-imin+1 not
-- being a divisor of Bigint. The algorithm should work, anyway.
return imin + (Random (gen) rem (imax - imin + 1));
end Random_Natural;
 
function Random_City_Location
(minloc : City_Location;
maxloc : City_Location)
return City_Location
is
begin
return City_Location (Random_Natural (minloc, maxloc));
end Random_City_Location;
 
function Random_Path_Index
(imin : Path_Index;
imax : Path_Index)
return Path_Index
is
begin
return Random_City_Location (imin, imax);
end Random_Path_Index;
 
package Natural_IO is new Ada.Text_IO.Integer_IO (Natural);
package City_Location_IO is new Ada.Text_IO.Integer_IO
(City_Location);
package Fixed_Point_IO is new Ada.Text_IO.Fixed_IO (Fixed_Point);
 
function sqrt
(x : Fixed_Point)
return Fixed_Point
is
begin
-- Cheat by using the floating-point routine. It is an exercise
-- for the reader to write a true fixed-point function.
return
Fixed_Point (Ada.Numerics.Elementary_Functions.Sqrt (Float (x)));
end sqrt;
 
function expneg
(x : Fixed_Point)
return Fixed_Point
is
begin
-- Cheat by using the floating-point routine. It is an exercise
-- for the reader to write a true fixed-point function.
return
Fixed_Point (Ada.Numerics.Elementary_Functions.Exp (-Float (x)));
end expneg;
 
function i_Coord
(loc : City_Location)
return Coordinate
is
begin
return loc / 10;
end i_Coord;
 
function j_Coord
(loc : City_Location)
return Coordinate
is
begin
return loc rem 10;
end j_Coord;
 
function Location
(i : Coordinate;
j : Coordinate)
return City_Location
is
begin
return (10 * i) + j;
end Location;
 
function distance
(loc1 : City_Location;
loc2 : City_Location)
return Fixed_Point
is
i1, j1 : Coordinate;
i2, j2 : Coordinate;
di, dj : Coordinate;
begin
i1 := i_Coord (loc1);
j1 := j_Coord (loc1);
i2 := i_Coord (loc2);
j2 := j_Coord (loc2);
di := (if i1 < i2 then i2 - i1 else i1 - i2);
dj := (if j1 < j2 then j2 - j1 else j1 - j2);
return sqrt (Fixed_Point ((di * di) + (dj * dj)));
end distance;
 
procedure Randomize_Path_Vector
(path : out Path_Vector)
is
j : Nonzero_Path_Index;
xi, xj : Nonzero_City_Location;
begin
for i in 0 .. 99 loop
path (i) := i;
end loop;
 
-- Do a Fisher-Yates shuffle of elements 1 .. 99.
for i in 1 .. 98 loop
j := Random_Path_Index (i + 1, 99);
xi := path (i);
xj := path (j);
path (i) := xj;
path (j) := xi;
end loop;
end Randomize_Path_Vector;
 
function Path_Length
(path : Path_Vector)
return Fixed_Point
is
len : Fixed_Point;
begin
len := distance (path (0), path (99));
for i in 0 .. 98 loop
len := len + distance (path (i), path (i + 1));
end loop;
return len;
end Path_Length;
 
-- Switch the index of s to switch which s is current and which is
-- the trial vector.
s : array (0 .. 1) of Path_Vector;
 
Current : Zero_or_One;
 
function Trial
return Zero_or_One
is
begin
return 1 - Current;
end Trial;
 
procedure Accept_Trial
is
begin
Current := 1 - Current;
end Accept_Trial;
 
procedure Find_Neighbors
(loc : City_Location;
neighbors : out Neighborhood_Array;
num_neighbors : out Integer)
is
i, j : Coordinate;
c1, c2, c3, c4, c5, c6, c7, c8 : City_Location := 0;
 
procedure Add_Neighbor
(neighbor : City_Location)
is
begin
if neighbor /= 0 then
num_neighbors := num_neighbors + 1;
neighbors (num_neighbors) := neighbor;
end if;
end Add_Neighbor;
 
begin
i := i_Coord (loc);
j := j_Coord (loc);
 
if i < 9 then
c1 := Location (i + 1, j);
if j < 9 then
c2 := Location (i + 1, j + 1);
end if;
if 0 < j then
c3 := Location (i + 1, j - 1);
end if;
end if;
if 0 < i then
c4 := Location (i - 1, j);
if j < 9 then
c5 := Location (i - 1, j + 1);
end if;
if 0 < j then
c6 := Location (i - 1, j - 1);
end if;
end if;
if j < 9 then
c7 := Location (i, j + 1);
end if;
if 0 < j then
c8 := Location (i, j - 1);
end if;
 
num_neighbors := 0;
Add_Neighbor (c1);
Add_Neighbor (c2);
Add_Neighbor (c3);
Add_Neighbor (c4);
Add_Neighbor (c5);
Add_Neighbor (c6);
Add_Neighbor (c7);
Add_Neighbor (c8);
end Find_Neighbors;
 
procedure Make_Neighbor_Path
is
u, v : City_Location;
neighbors : Neighborhood_Array;
num_neighbors : Integer;
j, iu, iv : Path_Index;
begin
for i in 0 .. 99 loop
s (Trial) := s (Current);
end loop;
 
u := Random_City_Location (1, 99);
Find_Neighbors (u, neighbors, num_neighbors);
v := neighbors (Random_Natural (1, num_neighbors));
 
j := 0;
iu := 0;
iv := 0;
while iu = 0 or iv = 0 loop
if s (Trial) (j + 1) = u then
iu := j + 1;
elsif s (Trial) (j + 1) = v then
iv := j + 1;
end if;
j := j + 1;
end loop;
s (Trial) (iu) := v;
s (Trial) (iv) := u;
end Make_Neighbor_Path;
 
function Temperature
(kT : Fixed_Point;
kmax : Natural;
k : Natural)
return Fixed_Point
is
begin
return
kT * (Fixed_Point (1) - (Fixed_Point (k) / Fixed_Point (kmax)));
end Temperature;
 
function Probability
(delta_E : Fixed_Point;
T : Fixed_Point)
return Fixed_Point
is
prob : Fixed_Point;
begin
if T = Fixed_Point (0.0) then
prob := Fixed_Point (0.0);
else
prob := expneg (delta_E / T);
end if;
return prob;
end Probability;
 
procedure Show
(k : Natural;
T : Fixed_Point;
E : Fixed_Point)
is
begin
Put (" ");
Natural_IO.Put (k, Width => 7);
Put (" ");
Fixed_Point_IO.Put (T, Fore => 5, Aft => 1);
Put (" ");
Fixed_Point_IO.Put (E, Fore => 7, Aft => 2);
Put_Line ("");
end Show;
 
procedure Display_Path
(path : Path_Vector)
is
begin
for i in 0 .. 99 loop
City_Location_IO.Put (path (i), Width => 2);
Put (" ->");
if i rem 8 = 7 then
Put_Line ("");
else
Put (" ");
end if;
end loop;
City_Location_IO.Put (path (0), Width => 2);
end Display_Path;
 
procedure Simulate_Annealing
(kT : Fixed_Point;
kmax : Natural)
is
kshow : Natural := kmax / 10;
E : Fixed_Point;
E_trial : Fixed_Point;
T : Fixed_Point;
P : Fixed_Point;
begin
E := Path_Length (s (Current));
for k in 0 .. kmax loop
T := Temperature (kT, kmax, k);
if k rem kshow = 0 then
Show (k, T, E);
end if;
Make_Neighbor_Path;
E_trial := Path_Length (s (Trial));
if E_trial <= E then
Accept_Trial;
E := E_trial;
else
P := Probability (E_trial - E, T);
if P = Fixed_Point (1) or else Randnum <= P then
Accept_Trial;
E := E_trial;
end if;
end if;
end loop;
end Simulate_Annealing;
 
kT : constant := Fixed_Point (1.0);
kmax : constant := 1_000_000;
 
begin
 
Reset (gen);
 
Current := 0;
Randomize_Path_Vector (s (Current));
 
Put_Line ("");
Put (" kT:");
Put_Line (Fixed_Point'Image (kT));
Put (" kmax:");
Put_Line (Natural'Image (kmax));
Put_Line ("");
Put_Line (" k T E(s)");
Simulate_Annealing (kT, kmax);
Put_Line ("");
Put_Line ("Final path:");
Display_Path (s (Current));
Put_Line ("");
Put_Line ("");
Put ("Final E(s): ");
Fixed_Point_IO.Put (Path_Length (s (Current)), Fore => 3, Aft => 2);
Put_Line ("");
Put_Line ("");
 
end simanneal;
 
----------------------------------------------------------------------</syntaxhighlight>
 
 
{{out}}
An example run. In the following, you could use gnatmake instead of gprbuild.
<pre>$ gprbuild -q simanneal.adb && ./simanneal
 
kT: 1.00000
kmax: 1000000
 
k T E(s)
0 1.0 525.23
100000 0.9 189.97
200000 0.8 180.33
300000 0.7 153.31
400000 0.6 156.18
500000 0.5 136.17
600000 0.4 119.56
700000 0.3 110.51
800000 0.2 106.21
900000 0.1 102.89
1000000 0.0 102.89
 
Final path:
0 -> 10 -> 11 -> 21 -> 20 -> 30 -> 31 -> 32 ->
22 -> 23 -> 33 -> 43 -> 42 -> 52 -> 51 -> 41 ->
40 -> 50 -> 60 -> 70 -> 80 -> 90 -> 91 -> 92 ->
93 -> 84 -> 94 -> 95 -> 85 -> 86 -> 96 -> 97 ->
98 -> 99 -> 89 -> 88 -> 87 -> 77 -> 67 -> 57 ->
58 -> 68 -> 78 -> 79 -> 69 -> 59 -> 49 -> 39 ->
29 -> 19 -> 9 -> 8 -> 7 -> 6 -> 25 -> 24 ->
34 -> 35 -> 44 -> 54 -> 53 -> 63 -> 62 -> 61 ->
71 -> 81 -> 72 -> 82 -> 83 -> 73 -> 74 -> 64 ->
65 -> 75 -> 76 -> 66 -> 56 -> 55 -> 45 -> 46 ->
47 -> 48 -> 38 -> 37 -> 36 -> 26 -> 27 -> 28 ->
18 -> 17 -> 16 -> 15 -> 5 -> 4 -> 14 -> 3 ->
13 -> 12 -> 2 -> 1 -> 0
 
Final E(s): 102.89
</pre>
 
=={{header|C}}==
Line 94 ⟶ 554:
For your platform you might have to modify parts, such as the call to getentropy(3).
 
You can also easily change the kind of floating point used. I apologize for false precision in printouts using the default single precision floating point.
 
Some might notice the calculations of random integers are done in a way that may introduce a bias, which is miniscule as long as the integer is much smaller than 2 to the 31st power. I mention this now so no one will complain about it later.
 
<langsyntaxhighlight lang="c">#include <math.h>
#include <stdio.h>
#include <stdlib.h>
Line 363 ⟶ 824:
 
return 0;
}</langsyntaxhighlight>
 
{{out}}
Line 405 ⟶ 866:
=={{header|C++}}==
'''Compiler:''' [[MSVC]] (19.27.29111 for x64)
<langsyntaxhighlight lang="c++">
#include<array>
#include<utility>
Line 578 ⟶ 1,039:
return 0;
}
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 606 ⟶ 1,067:
 
=={{header|EchoLisp}}==
<langsyntaxhighlight lang="scheme">
(lib 'math)
;; distances
Line 690 ⟶ 1,151:
(printf "E(s_final) %d" Emin)
(writeln 'Path s))
</syntaxhighlight>
</lang>
{{out}}
<pre>
Line 715 ⟶ 1,176:
34 33 32 43 42 52 51 41 31 21 11 12 22 23 13 14 15 16 17 26 27 37 38
39 29 28 18 19 9 8 7 6 5 4 3 2 1 0)
</pre>
 
=={{header|Fortran}}==
{{trans|Ada}}
{{works with|gfortran|11.3.0}}
 
 
<syntaxhighlight lang="fortran">module simanneal_support
implicit none
 
!
! The following two integer kinds are meant to be treated as
! synonyms.
!
! selected_int_kind (2) = integers in the range of at least -100 to
! +100.
!
integer, parameter :: city_location_kind = selected_int_kind (2)
integer, parameter :: path_index_kind = city_location_kind
 
!
! selected_int_kind (1) = integers in the range of at least -10 to
! +10.
!
integer, parameter :: coordinate_kind = selected_int_kind(1)
 
!
! selected_real_kind (6) = floating point with at least 6 decimal
! digits of precision.
!
integer, parameter :: float_kind = selected_real_kind (6)
 
!
! Shorthand notations.
!
integer, parameter :: clk = city_location_kind
integer, parameter :: pik = path_index_kind
integer, parameter :: cok = coordinate_kind
integer, parameter :: flk = float_kind
 
type path_vector
integer(kind = clk) :: elem(0:99)
end type path_vector
 
contains
 
function random_integer (imin, imax) result (n)
integer, intent(in) :: imin, imax
integer :: n
 
real(kind = flk) :: randnum
 
call random_number (randnum)
n = imin + floor ((imax - imin + 1) * randnum)
end function random_integer
 
function i_coord (loc) result (i)
integer(kind = clk), intent(in) :: loc
integer(kind = cok) :: i
 
i = loc / 10_clk
end function i_coord
function j_coord (loc) result (j)
integer(kind = clk), intent(in) :: loc
integer(kind = cok) :: j
 
j = mod (loc, 10_clk)
end function j_coord
 
function location (i, j) result (loc)
integer(kind = cok), intent(in) :: i, j
integer(kind = clk) :: loc
 
loc = (10_clk * i) + j
end function location
 
subroutine randomize_path_vector (path)
type(path_vector), intent(out) :: path
 
integer(kind = pik) :: i, j
integer(kind = clk) :: xi, xj
 
do i = 0_pik, 99_pik
path%elem(i) = i
end do
 
! Do a Fisher-Yates shuffle of elements 1 .. 99.
do i = 1_pik, 98_pik
j = int (random_integer (i + 1, 99), kind = pik)
xi = path%elem(i)
xj = path%elem(j)
path%elem(i) = xj
path%elem(j) = xi
end do
end subroutine randomize_path_vector
 
function distance (loc1, loc2) result (dist)
integer(kind = clk), intent(in) :: loc1, loc2
real(kind = flk) :: dist
 
integer(kind = cok) :: i1, j1
integer(kind = cok) :: i2, j2
integer :: di, dj
 
i1 = i_coord (loc1)
j1 = j_coord (loc1)
i2 = i_coord (loc2)
j2 = j_coord (loc2)
di = i1 - i2
dj = j1 - j2
dist = sqrt (real ((di * di) + (dj * dj), kind = flk))
end function distance
 
function path_length (path) result (len)
type(path_vector), intent(in) :: path
real(kind = flk) :: len
 
integer(kind = pik) :: i
 
len = distance (path%elem(0_pik), path%elem(99_pik))
do i = 0_pik, 98_pik
len = len + distance (path%elem(i), path%elem(i + 1_pik))
end do
end function path_length
 
subroutine find_neighbors (loc, neighbors, num_neighbors)
integer(kind = clk), intent(in) :: loc
integer(kind = clk), intent(out) :: neighbors(1:8)
integer, intent(out) :: num_neighbors
 
integer(kind = cok) :: i, j
integer(kind = clk) :: c1, c2, c3, c4, c5, c6, c7, c8
 
c1 = 0_clk
c2 = 0_clk
c3 = 0_clk
c4 = 0_clk
c5 = 0_clk
c6 = 0_clk
c7 = 0_clk
c8 = 0_clk
 
i = i_coord (loc)
j = j_coord (loc)
 
if (i < 9_cok) then
c1 = location (i + 1_cok, j)
if (j < 9_cok) then
c2 = location (i + 1_cok, j + 1_cok)
end if
if (0_cok < j) then
c3 = location (i + 1_cok, j - 1_cok)
end if
end if
if (0_cok < i) then
c4 = location (i - 1_cok, j)
if (j < 9_cok) then
c5 = location (i - 1_cok, j + 1_cok)
end if
if (0_cok < j) then
c6 = location (i - 1_cok, j - 1_cok)
end if
end if
if (j < 9_cok) then
c7 = location (i, j + 1_cok)
end if
if (0_cok < j) then
c8 = location (i, j - 1_cok)
end if
 
num_neighbors = 0
call add_neighbor (c1)
call add_neighbor (c2)
call add_neighbor (c3)
call add_neighbor (c4)
call add_neighbor (c5)
call add_neighbor (c6)
call add_neighbor (c7)
call add_neighbor (c8)
 
contains
 
subroutine add_neighbor (neighbor)
integer(kind = clk), intent(in) :: neighbor
 
if (neighbor /= 0_clk) then
num_neighbors = num_neighbors + 1
neighbors(num_neighbors) = neighbor
end if
end subroutine add_neighbor
 
end subroutine find_neighbors
 
function make_neighbor_path (path) result (neighbor_path)
type(path_vector), intent(in) :: path
type(path_vector) :: neighbor_path
 
integer(kind = clk) :: u, v
integer(kind = clk) :: neighbors(1:8)
integer :: num_neighbors
integer(kind = pik) :: j, iu, iv
 
neighbor_path = path
 
u = int (random_integer (1, 99), kind = clk)
call find_neighbors (u, neighbors, num_neighbors)
v = neighbors (random_integer (1, num_neighbors))
 
j = 0_pik
iu = 0_pik
iv = 0_pik
do while (iu == 0_pik .or. iv == 0_pik)
if (neighbor_path%elem(j + 1) == u) then
iu = j + 1
else if (neighbor_path%elem(j + 1) == v) then
iv = j + 1
end if
j = j + 1
end do
neighbor_path%elem(iu) = v
neighbor_path%elem(iv) = u
end function make_neighbor_path
 
function temperature (kT, kmax, k) result (temp)
real(kind = flk), intent(in) :: kT
integer, intent(in) :: kmax, k
real(kind = flk) :: temp
 
real(kind = flk) :: kf, kmaxf
 
kf = real (k, kind = flk)
kmaxf = real (kmax, kind = flk)
temp = kT * (1.0_flk - (kf / kmaxf))
end function temperature
 
function probability (delta_E, T) result (prob)
real(kind = flk), intent(in) :: delta_E, T
real(kind = flk) :: prob
 
if (T == 0.0_flk) then
prob = 0.0_flk
else
prob = exp (-(delta_E / T))
end if
end function probability
 
subroutine show (k, T, E)
integer, intent(in) :: k
real(kind = flk), intent(in) :: T, E
 
write (*, 10) k, T, E
10 format (1X, I7, 1X, F7.1, 1X, F10.2)
end subroutine show
 
subroutine display_path (path)
type(path_vector), intent(in) :: path
 
integer(kind = pik) :: i
 
999 format ()
100 format (' ->')
110 format (' ')
120 format (I2)
 
do i = 0_pik, 99_pik
write (*, 120, advance = 'no') path%elem(i)
write (*, 100, advance = 'no')
if (mod (i, 8_pik) == 7_pik) then
write (*, 999, advance = 'yes')
else
write (*, 110, advance = 'no')
end if
end do
write (*, 120, advance = 'no') path%elem(0_pik)
end subroutine display_path
 
subroutine simulate_annealing (kT, kmax, initial_path, final_path)
real(kind = flk), intent(in) :: kT
integer, intent(in) :: kmax
type(path_vector), intent(in) :: initial_path
type(path_vector), intent(inout) :: final_path
 
integer :: kshow
integer :: k
real(kind = flk) :: E, E_trial, T
type(path_vector) :: path, trial
real(kind = flk) :: randnum
 
kshow = kmax / 10
 
path = initial_path
E = path_length (path)
do k = 0, kmax
T = temperature (kT, kmax, k)
if (mod (k, kshow) == 0) call show (k, T, E)
trial = make_neighbor_path (path)
E_trial = path_length (trial)
if (E_trial <= E) then
path = trial
E = E_trial
else
call random_number (randnum)
if (randnum <= probability (E_trial - E, T)) then
path = trial
E = E_trial
end if
end if
end do
final_path = path
end subroutine simulate_annealing
 
end module simanneal_support
 
program simanneal
 
use, non_intrinsic :: simanneal_support
implicit none
 
real(kind = flk), parameter :: kT = 1.0_flk
integer, parameter :: kmax = 1000000
 
type(path_vector) :: initial_path
type(path_vector) :: final_path
 
call random_seed
 
call randomize_path_vector (initial_path)
 
10 format ()
20 format (' kT: ', F0.2)
30 format (' kmax: ', I0)
40 format (' k T E(s)')
50 format (' --------------------------')
60 format ('Final E(s): ', F0.2)
 
write (*, 10)
write (*, 20) kT
write (*, 30) kmax
write (*, 10)
write (*, 40)
write (*, 50)
call simulate_annealing (kT, kmax, initial_path, final_path)
write (*, 10)
call display_path (final_path)
write (*, 10)
write (*, 10)
write (*, 60) path_length (final_path)
write (*, 10)
 
end program simanneal</syntaxhighlight>
 
 
{{out}}
<pre>$ gfortran -std=f2018 -Ofast simanneal.f90 && ./a.out
 
kT: 1.00
kmax: 1000000
 
k T E(s)
--------------------------
0 1.0 517.11
100000 0.9 198.12
200000 0.8 169.43
300000 0.7 164.66
400000 0.6 149.10
500000 0.5 138.38
600000 0.4 119.24
700000 0.3 113.69
800000 0.2 105.80
900000 0.1 101.66
1000000 0.0 101.66
 
0 -> 10 -> 11 -> 21 -> 31 -> 20 -> 30 -> 40 ->
41 -> 51 -> 50 -> 60 -> 70 -> 71 -> 61 -> 62 ->
72 -> 82 -> 81 -> 80 -> 90 -> 91 -> 92 -> 93 ->
83 -> 73 -> 74 -> 84 -> 94 -> 95 -> 96 -> 97 ->
98 -> 99 -> 89 -> 88 -> 79 -> 69 -> 59 -> 58 ->
48 -> 49 -> 39 -> 38 -> 28 -> 29 -> 19 -> 9 ->
8 -> 18 -> 17 -> 7 -> 6 -> 16 -> 15 -> 5 ->
4 -> 14 -> 24 -> 25 -> 26 -> 27 -> 37 -> 36 ->
35 -> 45 -> 46 -> 47 -> 57 -> 67 -> 68 -> 78 ->
77 -> 87 -> 86 -> 85 -> 75 -> 76 -> 66 -> 56 ->
55 -> 65 -> 64 -> 63 -> 54 -> 53 -> 52 -> 42 ->
43 -> 44 -> 34 -> 33 -> 32 -> 22 -> 23 -> 12 ->
13 -> 3 -> 2 -> 1 -> 0
 
Final E(s): 101.66
</pre>
 
=={{header|FreeBASIC}}==
Uses 'LCS' function from [[Longest common subsequence#FreeBASIC]]:
 
<syntaxhighlight lang="vbnet">Dim Shared As Double dists(0 To 9999)
 
' index into lookup table of Nums
Function dist(ci As Integer, cj As Integer) As Double
Return dists(cj*100 + ci)
End Function
 
' energy at s, to be minimized
Function Ens(path() As Integer) As Double
Dim As Double d = 0
For i As Integer = 0 To Ubound(path) - 1
d += dist(path(i), path(i+1))
Next
Return d
End Function
 
' temperature function, decreases to 0
Function T(k As Double, kmax As Double, kT As Double) As Double
Return (1 - k / kmax) * kT
End Function
 
' variation of E, from state s to state s_next
Function dE(s() As Integer, u As Integer, v As Integer) As Double
Dim As Integer su = s(u)
Dim As Integer sv = s(v)
' old
Dim As Double a = dist(s(u-1), su)
Dim As Double b = dist(s(u+1), su)
Dim As Double c = dist(s(v-1), sv)
Dim As Double d = dist(s(v+1), sv)
' new
Dim As Double na = dist(s(u-1), sv)
Dim As Double nb = dist(s(u+1), sv)
Dim As Double nc = dist(s(v-1), su)
Dim As Double nd = dist(s(v+1), su)
If v = u+1 Then Return (na + nd) - (a + d)
If u = v+1 Then Return (nc + nb) - (c + b)
Return (na + nb + nc + nd) - (a + b + c + d)
End Function
 
' probability to move from s to s_next
Function P(deltaE As Double, k As Double, kmax As Double, kT As Double) As Double
Return Exp(-deltaE / T(k, kmax, kT))
End Function
 
' Simulated annealing
Sub sa(kmax As Double, kT As Double)
Dim As Integer s(0 To 100)
Dim As Integer temp(0 To 98)
Dim As Integer dirs(0 To 7) = {1, -1, 10, -10, 9, 11, -11, -9}
Dim As Integer i, k, u, v, cv
Dim As Double Emin
For i = 0 To 98
temp(i) = i + 1
Next
Randomize Timer
For i = 0 To 98
Swap temp(i), temp(Int(Rnd * 99))
Next
For i = 0 To 98
s(i+1) = temp(i)
Next
Print "kT = "; kT
Print "E(s0) "; Ens(s())
Print
Emin = Ens(s())
For k = 0 To kmax
If k Mod (kmax/10) = 0 Then
Print Using "k: ####### T: #.#### Es: ###.####"; k; T(k, kmax, kT); Ens(s())
End If
u = Int(Rnd * 99) + 1
cv = s(u) + dirs(Int(Rnd * 8))
If cv <= 0 Or cv >= 100 Then Continue For
If Abs(dist(s(u), cv)) > 5 Then Continue For
v = s(cv)
Dim As Double deltae = dE(s(), u, v)
If deltae < 0 Or P(deltae, k, kmax, kT) >= Rnd Then
Swap s(u), s(v)
Emin = Emin + deltae
End If
Next k
Print
Print "E(s_final) "; Emin
Print "Path:"
For i = 0 To Ubound(s)
If i > 0 And i Mod 10 = 0 Then Print
Print Using "####"; s(i);
Next
Print
End Sub
 
' distances
For i As Integer = 0 To 9999
Dim As Integer ab = (i \ 100)
Dim As Integer cd = i Mod 100
Dim As Integer a = (ab \ 10)
Dim As Integer b = ab Mod 10
Dim As Integer c = (cd \ 10)
Dim As Integer d = cd Mod 10
dists(i) = Sqr((a-c)^2 + (b-d)^2)
Next i
 
Dim As Double kT = 1, kmax = 1e6
sa(kmax, kT)
 
Sleep</syntaxhighlight>
{{out}}
<pre>kT = 1
E(s0) 510.1804163299929
 
k: 0 T: 1.0000 Es: 510.1804
k: 100000 T: 0.9000 Es: 195.1253
k: 200000 T: 0.8000 Es: 182.4579
k: 300000 T: 0.7000 Es: 153.4156
k: 400000 T: 0.6000 Es: 150.7938
k: 500000 T: 0.5000 Es: 141.6804
k: 600000 T: 0.4000 Es: 128.4290
k: 700000 T: 0.3000 Es: 123.2713
k: 800000 T: 0.2000 Es: 117.4202
k: 900000 T: 0.1000 Es: 116.0060
k: 1000000 T: 0.0000 Es: 116.0060
 
E(s_final) 116.0060090954848
Path:
0 11 10 20 21 32 22 12 2 3
13 14 34 33 23 24 35 25 16 15
4 5 6 7 9 8 18 19 29 39
49 48 38 28 27 17 26 36 47 37
45 46 57 56 55 54 44 43 42 52
51 41 31 30 40 50 60 61 83 73
63 62 72 71 70 80 90 91 81 82
92 93 94 96 97 98 99 89 79 69
59 58 68 67 77 87 88 78 76 66
65 75 86 95 85 84 74 64 53 1
0
</pre>
 
=={{header|Go}}==
{{trans|zkl}}
<langsyntaxhighlight lang="go">package main
 
import (
Line 833 ⟶ 1,823:
func main() {
sa(1e6, 1)
}</langsyntaxhighlight>
 
{{out}}
Line 866 ⟶ 1,856:
81 71 70 60 50 40 30 20 10 1
0
</pre>
 
=={{header|Icon}}==
{{trans|Fortran}}
 
 
<syntaxhighlight lang="icon">link printf
link random
 
procedure main ()
local initial_path
local final_path
local kT, kmax
 
randomize()
 
kT := 1.0
kmax := 1000000
 
write()
write(" kT: ", kT)
write(" kmax: ", kmax)
write()
write(" k T E(s)")
write(" --------------------------")
initial_path := randomize_path_vector()
final_path := simulate_annealing (kT, kmax, initial_path)
write()
display_path (final_path)
write()
write()
printf("Final E(s): %.2r\n", path_length(final_path))
write()
end
 
procedure randomize_path_vector ()
local path
local i, j
 
path := []
every put (path, 0 to 99)
 
# Shuffle elements 2 to 0.
every i := 1 to 98 do {
j := ?(99 - i) + i + 1
path[i + 1] :=: path[j + 1]
}
 
return path
end
 
procedure distance (loc1, loc2)
local i1, j1
local i2, j2
local di, dj
 
i1 := loc1 / 10
j1 := loc1 % 10
i2 := loc2 / 10
j2 := loc2 % 10
di := i1 - i2
dj := j1 - j2
return sqrt ((di * di) + (dj * dj))
end
 
procedure path_length (path)
local i
local len
 
len := distance(path[1], path[100])
every i := 1 to 99 do len +:= distance(path[i], path[i + 1])
return len
end
 
procedure find_neighbors (loc)
local c1, c2, c3, c4, c5, c6, c7, c8
local i, j
local neighbors
 
c1 := c2 := c3 := c4 := c5 := c6 := c7 := c8 := 0
 
i := loc / 10
j := loc % 10
 
if (i < 9) then {
c1 := (10 * (i + 1)) + j
if (j < 9) then c2 := (10 * (i + 1)) + (j + 1)
if (0 < j) then c3 := (10 * (i + 1)) + (j - 1)
}
if (0 < i) then {
c4 := (10 * (i - 1)) + j
if (j < 9) then c5 := (10 * (i - 1)) + (j + 1)
if (0 < j) then c6 := (10 * (i - 1)) + (j - 1)
}
if (j < 9) then c7 := (10 * i) + (j + 1)
if (0 < j) then c8 := (10 * i) + (j - 1)
 
neighbors := []
every put(neighbors, 0 ~= (c1 | c2 | c3 | c4 | c5 | c6 | c7 | c8))
return neighbors
end
 
procedure make_neighbor_path (path)
local neighbor_path
local u, v, iu, iv, j
local neighbors
 
neighbor_path := copy(path)
 
u := ?99
neighbors := find_neighbors(u)
v := neighbors[?(*neighbors)]
 
j := 2
iu := 0
iv := 0
while iu = 0 | iv = 0 do {
if neighbor_path[j] = u then {
iu := j
} else if neighbor_path[j] = v then {
iv := j
}
j +:= 1
}
neighbor_path[iu] := v
neighbor_path[iv] := u
 
return neighbor_path
end
 
procedure temperature (kT, kmax, k)
return kT * (1.0 - (real(k) / real(kmax)))
end
 
procedure my_exp (x)
# Icon's exp() might bail out with an underflow error, if we are not
# careful.
return (if x < -50 then 0.0 else exp(x))
end
 
procedure probability (delta_E, T)
return (if T = 0.0 then 0.0 else my_exp(-(delta_E / T)))
end
 
procedure show (k, T, E)
printf(" %7d %7.1r %10.2r\n", k, T, E)
return
end
 
procedure display_path (path)
local i
 
every i := 1 to 100 do {
printf("%2d ->", path[i])
if ((i - 1) % 8) = 7 then {
write()
} else {
writes(" ")
}
}
printf("%2d", path[1])
return
end
 
procedure simulate_annealing (kT, kmax, path)
local kshow
local k
local E, E_trial, T
local trial
 
kshow := kmax / 10
 
E := path_length(path)
every k := 0 to kmax do {
T := temperature(kT, kmax, k)
if (k % kshow) = 0 then show(k, T, E)
trial := make_neighbor_path(path)
E_trial := path_length(trial)
if E_trial <= E | ?0 <= probability (E_trial - E, T) then {
path := trial
E := E_trial
}
}
return path
end</syntaxhighlight>
 
{{out}}
An example run:
<pre>$ icont -s -u simanneal-in-Icon.icn && ./simanneal-in-Icon
 
kT: 1.0
kmax: 1000000
 
k T E(s)
--------------------------
0 1.0 511.67
100000 0.9 206.16
200000 0.8 186.68
300000 0.7 165.92
400000 0.6 158.49
500000 0.5 141.76
600000 0.4 122.53
700000 0.3 119.47
800000 0.2 107.56
900000 0.1 102.89
1000000 0.0 102.24
 
0 -> 10 -> 20 -> 30 -> 31 -> 41 -> 40 -> 50 ->
60 -> 70 -> 71 -> 72 -> 62 -> 61 -> 51 -> 52 ->
53 -> 63 -> 54 -> 44 -> 45 -> 35 -> 34 -> 24 ->
25 -> 26 -> 27 -> 17 -> 7 -> 8 -> 9 -> 19 ->
29 -> 39 -> 49 -> 59 -> 69 -> 79 -> 89 -> 99 ->
98 -> 97 -> 96 -> 86 -> 76 -> 75 -> 84 -> 85 ->
95 -> 94 -> 93 -> 92 -> 91 -> 90 -> 80 -> 81 ->
82 -> 83 -> 73 -> 74 -> 64 -> 55 -> 65 -> 66 ->
56 -> 46 -> 36 -> 37 -> 47 -> 57 -> 67 -> 77 ->
87 -> 88 -> 78 -> 68 -> 58 -> 48 -> 38 -> 28 ->
18 -> 16 -> 6 -> 5 -> 15 -> 14 -> 4 -> 3 ->
2 -> 12 -> 13 -> 23 -> 33 -> 43 -> 42 -> 32 ->
22 -> 21 -> 11 -> 1 -> 0
 
Final E(s): 102.24
</pre>
 
Line 872 ⟶ 2,084:
Implementation:
 
<langsyntaxhighlight Jlang="j">dist=: +/&.:*:@:-"1/~10 10#:i.100
 
satsp=:4 :0
Line 897 ⟶ 2,109:
end.
0,s,0
)</langsyntaxhighlight>
 
Notes:
Line 913 ⟶ 2,125:
Sample run:
 
<langsyntaxhighlight Jlang="j"> 1e6 satsp dist
0 1 538.409
100000 0.9 174.525
Line 925 ⟶ 2,137:
900000 0.1 101.657
1e6 0 101.657
0 1 2 3 4 13 23 24 34 44 43 33 32 31 41 42 52 51 61 62 53 54 64 65 55 45 35 25 15 14 5 6 7 17 16 26 27 37 36 46 47 48 38 28 18 8 9 19 29 39 49 59 69 79 78 68 58 57 56 66 67 77 76 75 85 86 87 88 89 99 98 97 96 95 94 84 74 73 63 72 82 83 93 92 91 90 80 81 71 70 60 50 40 30 20 21 22 12 11 10 0</langsyntaxhighlight>
 
=={{header|jq}}==
'''Adapted from [[#Wren|Wren]]'''
 
'''Works with jq, the C implementation of jq'''
 
'''Works with gojq, the Go implementation of jq'''
 
This adaptation does not cache the distances
and can be used for any square grid of cities.
 
Since jq does not include a PRN generator, we assume an
external source of randomness, such as /dev/urandom.
Specifically, the following program assumes an invocation
of jq along the lines of:
<pre>
< /dev/urandom tr -cd '0-9' | fold -w 1 | jq -Rcnr -f sa.jq
</pre>
 
Since gojq does not include jq's `_nwise/1`, here is a suitable def:
<pre>
# Require $n > 0
def _nwise($n):
def _n: if length <= $n then . else .[:$n] , (.[$n:] | _n) end;
if $n <= 0 then "_nwise: argument should be non-negative" else _n end;
</pre>
<syntaxhighlight lang="jq">
## Pseuo-random numbers and shuffling
 
# Output: a prn in range(0;$n) where $n is `.`
def prn:
if . == 1 then 0
else . as $n
| ([1, (($n-1)|tostring|length)]|max) as $w
| [limit($w; inputs)] | join("") | tonumber
| if . < $n then . else ($n | prn) end
end;
 
def randFloat:
(1000|prn) / 1000;
 
def knuthShuffle:
length as $n
| if $n <= 1 then .
else {i: $n, a: .}
| until(.i == 0;
.i += -1
| (.i + 1 | prn) as $j
| .a[.i] as $t
| .a[.i] = .a[$j]
| .a[$j] = $t)
| .a
end;
 
 
## Generic utilities
def divmod($j):
(. % $j) as $mod
| [(. - $mod) / $j, $mod] ;
 
def hypot($a;$b):
($a*$a) + ($b*$b) | sqrt;
 
def lpad($len): tostring | ($len - length) as $l | (" " * $l) + .;
 
def round($ndec): pow(10;$ndec) as $p | . * $p | round / $p;
 
def sum(s): reduce s as $x (0; . + $x);
 
def swap($i; $j):
.[$i] as $tmp
| .[$i] = .[$j]
| .[$j] = $tmp;
 
 
### The cities
 
# all 8 neighbors for an $n x $n grid
def neighbors($n): [1, -1, $n, -$n, $n-1, $n+1, -$n-1, $n+1];
 
# Distance between two cities $x and $y in an .n * .n grid
def dist($x; $y):
.n as $n
| ($x | divmod($n)) as [$xi, $xj]
| ($y | divmod($n)) as [$yi, $yj]
| hypot( $xi-$yi; $xj - $yj );
 
 
### Simulated annealing
 
# The energy of the input state (.s), to be minimized
# Input: {s, n}
def Es:
.s as $path
| sum( range(0; $path|length - 1) as $i
| dist($path[$i]; $path[$i+1]) );
 
# temperature function, decreases to 0
def T($k; $kmax; $kT):
(1 - ($k / $kmax)) * $kT;
 
# variation of E, from one state to the next state
# Input: {s, n}
def dE($u; $v):
.s as $s
| $s[$u] as $su
| $s[$v] as $sv
# old
| dist($s[$u-1]; $su) as $a
| dist($s[$u+1]; $su) as $b
| dist($s[$v-1]; $sv) as $c
| dist($s[$v+1]; $sv) as $d
# new
| dist($s[$u-1]; $sv) as $na
| dist($s[$u+1]; $sv) as $nb
| dist($s[$v-1]; $su) as $nc
| dist($s[$v+1]; $su) as $nd
| if ($v == $u+1) then ($na + $nd) - ($a + $d)
elif ($u == $v+1) then ($nc + $nb) - ($c + $b)
else ($na + $nb + $nc + $nd) - ($a + $b + $c + $d)
end;
 
# probability of moving from one state to another
def P($deltaE; $k; $kmax; $kT):
T($k; $kmax; $kT) as $T
| if $T == 0 then 0
else (-$deltaE / $T) | exp
end;
 
# Simulated annealing for $n x $n cities
def sa($kmax; $kT; $n):
def format($k; $T; $E):
[ "k:", ($k | lpad(10)),
"T:", ($T | round(2) | lpad(4)),
"Es:", $E ]
| join(" ");
 
neighbors($n) as $neighbors # potential neighbors
| ($n*$n) as $n2
# random path from 0 to 0
| {s: ([0] + ([ range(1; $n2)] | knuthShuffle) + [0]) }
| .n = $n # for dist/2
| .Emin = Es # E0
| "kT = \($kT)",
"E(s0) \(.Emin)\n",
( foreach range(0; 1+$kmax) as $k (.;
.emit = null
| if ($k % (($kmax/10)|floor)) == 0
then .emit = format($k; T($k; $kmax; $kT); Es)
else .
end
| (($n2-1)|prn + 1) as $u # a random city apart from the starting point
| (.s[$u] + $neighbors[8|prn]) as $cv # a neighboring city, perhaps
| if ($cv <= 0 or $cv >= $n2) # check the city is not bogus
then . # continue
elif dist(.s[$u]; $cv) > 5 # check true neighbor
then . # continue
else .s[$cv] as $v # city index
| dE($u; $v) as $deltae
| if ($deltae < 0 or # always move if negative
P($deltae; $k; $kmax; $kT) >= randFloat)
then .s |= swap($u; $v)
| .Emin += $deltae
end
end;
 
select(.emit).emit,
(select($k == $kmax)
| "\nE(s_final) \(.Emin)",
"Path:",
# output final state
(.s | map(lpad(3)) | _nwise(10) | join(" ")) ) ));
 
# Cities on a 10 x 10 grid
sa(1e6; 1; 10)
</syntaxhighlight>
{{output}}
<pre>
kT = 1
E(s0) 511.63434626356127
 
k: 0 T: 1 Es: 511.63434626356127
k: 100000 T: 0.9 Es: 183.44842684951274
k: 200000 T: 0.8 Es: 173.6522166458839
k: 300000 T: 0.7 Es: 191.88956498870922
k: 400000 T: 0.6 Es: 161.63509965859427
k: 500000 T: 0.5 Es: 173.6829125726551
k: 600000 T: 0.4 Es: 135.5154326151275
k: 700000 T: 0.3 Es: 174.33930236055193
k: 800000 T: 0.2 Es: 141.907500599355
k: 900000 T: 0.1 Es: 141.76740977979034
k: 1000000 T: 0 Es: 148.13930861301918
 
E(s_final) 148.13930861301935
Path:
0 1 2 11 10 20 21 22 23 24
14 5 4 3 13 12 32 42 33 25
15 16 6 7 8 9 19 29 39 28
18 17 27 26 36 46 35 34 43 52
41 30 31 44 45 55 56 38 37 49
48 47 65 64 54 53 51 72 91 81
80 71 70 61 62 40 50 60 92 82
83 73 63 57 67 66 75 74 84 93
94 95 78 77 68 58 87 76 86 99
89 79 69 59 88 98 97 96 85 90
0
</pre>
 
=={{header|Julia}}==
Line 931 ⟶ 2,350:
 
'''Module''':
<langsyntaxhighlight lang="julia">module TravelingSalesman
 
using Random, Printf
Line 1,006 ⟶ 2,425:
end
 
end # module TravelingSalesman</langsyntaxhighlight>
 
'''Main''':
<langsyntaxhighlight lang="julia">distance(a, b) = sqrt(sum((a .- b) .^ 2))
const _citydist = collect(distance((ci % 10, ci ÷ 10), (cj % 10, cj ÷ 10)) for ci in 1:100, cj in 1:100)
 
TravelingSalesman.findpath(_citydist, 1_000_000, 1)</langsyntaxhighlight>
 
{{out}}
Line 1,055 ⟶ 2,474:
 
=={{header|Nim}}==
<langsyntaxhighlight Nimlang="nim">import math, random, sequtils, strformat
 
const
Line 1,137 ⟶ 2,556:
echo fmt"path: {s}"
 
main()</langsyntaxhighlight>
 
Compile and run: <pre>nim c -r -d:release --opt:speed travel_sa.nim</pre>
Line 1,159 ⟶ 2,578:
=={{header|Perl}}==
{{trans|Raku}}
<langsyntaxhighlight lang="perl">use utf8;
use strict;
use warnings;
Line 1,237 ⟶ 2,656:
printf "@{['%4d' x 20]}\n", @s[$l*20 .. ($l+1)*20 - 1];
}
printf " 0\n";</langsyntaxhighlight>
{{out}}
<pre>k: 0 T: 1.0 Es: 519.2
Line 1,260 ⟶ 2,679:
=={{header|Phix}}==
{{trans|zkl}}
<!--<langsyntaxhighlight Phixlang="phix">(phixonline)-->
<span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span>
<span style="color: #008080;">function</span> <span style="color: #000000;">hypot</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">return</span> <span style="color: #7060A8;">sqrt</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">*</span><span style="color: #000000;">a</span><span style="color: #0000FF;">+</span><span style="color: #000000;">b</span><span style="color: #0000FF;">*</span><span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
Line 1,341 ⟶ 2,760:
<span style="color: #008080;">end</span> <span style="color: #008080;">procedure</span>
<span style="color: #000000;">sa</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1_000_000</span><span style="color: #0000FF;">,</span><span style="color: #000000;">1</span><span style="color: #0000FF;">)</span>
<!--</langsyntaxhighlight>-->
{{out}}
<pre>
Line 1,366 ⟶ 2,785:
 
=={{header|Racket}}==
<langsyntaxhighlight lang="racket">
#lang racket
(require racket/fixnum)
Line 1,479 ⟶ 2,898:
 
(module+ main
(Simulated-annealing))</langsyntaxhighlight>
{{out}}
<pre>T:1 E:552.4249706051347
Line 1,497 ⟶ 2,916:
(formerly Perl 6)
{{trans|Go}}
<syntaxhighlight lang="raku" perl6line># simulation setup
my \cities = 100; # number of cities
my \kmax = 1e6; # iterations to run
Line 1,551 ⟶ 2,970:
 
say "\nE(s_final): " ~ E-min-global.fmt('%.1f');
say "Path:\n" ~ s».fmt('%2d').rotor(20,:partial).join: "\n";</langsyntaxhighlight>
{{out}}
<pre>k: 0 T: 1.0 Es: 522.0
Line 1,572 ⟶ 2,991:
28 38 48 49 39 29 19 9 8 7 6 5 4 14 13 12 11 2 3 1
0</pre>
 
=={{header|RATFOR}}==
{{trans|Fortran}}
 
{{works with|ratfor77|[https://sourceforge.net/p/chemoelectric/ratfor77/ public domain 1.0]}}
{{works with|gfortran|11.3.0}}
 
 
<syntaxhighlight lang="ratfor">#
# The Rosetta Code simulated annealing task, in Ratfor 77.
#
# This implementation uses the RANDOM_NUMBER intrinsic and therefore
# will not work with f2c. It will work with gfortran. (One could
# substitute a random number generator from the Fullerton Function
# Library, or from elsewhere.)
#
 
function rndint (imin, imax)
implicit none
 
integer imin, imax, rndint
 
real rndnum
 
call random_number (rndnum)
rndint = imin + floor ((imax - imin + 1) * rndnum)
end
 
function icoord (loc)
implicit none
 
integer loc, icoord
 
icoord = loc / 10
end
 
function jcoord (loc)
implicit none
 
integer loc, jcoord
 
jcoord = mod (loc, 10)
end
 
function locatn (i, j) # Location.
implicit none
 
integer i, j, locatn
 
locatn = (10 * i) + j
end
 
subroutine rndpth (path) # Randomize a path.
implicit none
 
integer path(0:99)
 
integer rndint
 
integer i, j, xi, xj
 
for (i = 0; i <= 99; i = i + 1)
path(i) = i
 
# Fisher-Yates shuffle of elements 1 .. 99.
for (i = 1; i <= 98; i = i + 1)
{
j = rndint (i + 1, 99)
xi = path(i)
xj = path(j)
path(i) = xj
path(j) = xi
}
end
 
function dstnce (loc1, loc2) # Distance.
implicit none
 
integer loc1, loc2
real dstnce
 
integer icoord, jcoord
 
integer i1, j1
integer i2, j2
integer di, dj
 
i1 = icoord (loc1)
j1 = jcoord (loc1)
i2 = icoord (loc2)
j2 = jcoord (loc2)
di = i1 - i2
dj = j1 - j2
dstnce = sqrt (real ((di * di) + (dj * dj)))
end
 
function pthlen (path) # Path length.
implicit none
 
integer path(0:99)
real pthlen
 
real dstnce
 
real len
integer i
 
len = dstnce (path(0), path(99))
for (i = 0; i <= 98; i = i + 1)
len = len + dstnce (path(i), path(i + 1))
pthlen = len
end
 
subroutine addnbr (nbrs, numnbr, nbr) # Add neighbor.
implicit none
 
integer nbrs(1:8)
integer numnbr
integer nbr
 
if (nbr != 0)
{
numnbr = numnbr + 1
nbrs(numnbr) = nbr
}
end
 
subroutine fndnbr (loc, nbrs, numnbr) # Find neighbors.
implicit none
 
integer loc
integer nbrs(1:8)
integer numnbr
 
integer icoord, jcoord
integer locatn
 
integer i, j
integer c1, c2, c3, c4, c5, c6, c7, c8
 
c1 = 0
c2 = 0
c3 = 0
c4 = 0
c5 = 0
c6 = 0
c7 = 0
c8 = 0
 
i = icoord (loc)
j = jcoord (loc)
 
if (i < 9)
{
c1 = locatn (i + 1, j)
if (j < 9)
c2 = locatn (i + 1, j + 1)
if (0 < j)
c3 = locatn (i + 1, j - 1)
}
if (0 < i)
{
c4 = locatn (i - 1, j)
if (j < 9)
c5 = locatn (i - 1, j + 1)
if (0 < j)
c6 = locatn (i - 1, j - 1)
}
if (j < 9)
c7 = locatn (i, j + 1)
if (0 < j)
c8 = locatn (i, j - 1)
 
numnbr = 0
call addnbr (nbrs, numnbr, c1)
call addnbr (nbrs, numnbr, c2)
call addnbr (nbrs, numnbr, c3)
call addnbr (nbrs, numnbr, c4)
call addnbr (nbrs, numnbr, c5)
call addnbr (nbrs, numnbr, c6)
call addnbr (nbrs, numnbr, c7)
call addnbr (nbrs, numnbr, c8)
end
 
subroutine nbrpth (path, nbrp) # Make a neighbor path.
implicit none
 
integer path(0:99), nbrp(0:99)
 
integer rndint
 
integer u, v
integer nbrs(1:8)
integer numnbr
integer j, iu, iv
 
for (j = 0; j <= 99; j = j + 1)
nbrp(j) = path(j)
 
u = rndint (1, 99)
call fndnbr (u, nbrs, numnbr)
v = nbrs(rndint (1, numnbr))
 
j = 1
iu = 0
iv = 0
while (iu == 0 || iv == 0)
{
if (nbrp(j) == u)
iu = j
else if (nbrp(j) == v)
iv = j
j = j + 1
}
nbrp(iu) = v
nbrp(iv) = u
end
 
function temp (kT, kmax, k) # Temperature.
implicit none
 
real kT
integer kmax, k
real temp
 
real kf, kmaxf
 
kf = real (k)
kmaxf = real (kmax)
temp = kT * (1.0 - (kf / kmaxf))
end
 
function prob (deltaE, T) # Probability.
implicit none
 
real deltaE, T, prob
real x
 
if (T == 0.0)
prob = 0.0
else
{
x = -(deltaE / T)
if (x < -80)
prob = 0 # Avoid underflow.
else
prob = exp (-(deltaE / T))
}
end
 
subroutine show (k, T, E)
implicit none
 
integer k
real T, E
 
10 format (1X, I7, 1X, F7.1, 1X, F10.2)
 
write (*, 10) k, T, E
end
 
subroutine dsplay (path)
implicit none
 
integer path(0:99)
 
100 format (8(I2, ' -> '))
 
write (*, 100) path
end
 
subroutine sa (kT, kmax, path)
implicit none
 
real kT
integer kmax
integer path(0:99)
 
real pthlen
real temp, prob
 
integer kshow
integer k
integer j
real E, Etrial, T
integer trial(0:99)
real rndnum
 
kshow = kmax / 10
 
E = pthlen (path)
for (k = 0; k <= kmax; k = k + 1)
{
T = temp (kT, kmax, k)
if (mod (k, kshow) == 0)
call show (k, T, E)
call nbrpth (path, trial)
Etrial = pthlen (trial)
if (Etrial <= E)
{
for (j = 0; j <= 99; j = j + 1)
path(j) = trial(j)
E = Etrial
}
else
{
call random_number (rndnum)
if (rndnum <= prob (Etrial - E, T))
{
for (j = 0; j <= 99; j = j + 1)
path(j) = trial(j)
E = Etrial
}
}
}
end
 
program simanl
implicit none
 
real pthlen
 
integer path(0:99)
real kT
integer kmax
 
kT = 1.0
kmax = 1000000
 
10 format ()
20 format (' kT: ', F0.2)
30 format (' kmax: ', I0)
40 format (' k T E(s)')
50 format (' --------------------------')
60 format ('Final E(s): ', F0.2)
 
write (*, 10)
write (*, 20) kT
write (*, 30) kmax
write (*, 10)
write (*, 40)
write (*, 50)
call rndpth (path)
call sa (kT, kmax, path)
write (*, 10)
call dsplay (path)
write (*, 10)
write (*, 60) pthlen (path)
write (*, 10)
end</syntaxhighlight>
 
{{out}}
An example run:
<pre>$ ratfor77 simanneal.r > sa.f && gfortran -O3 -std=legacy sa.f && ./a.out
 
kT: 1.00
kmax: 1000000
 
k T E(s)
--------------------------
0 1.0 547.76
100000 0.9 190.62
200000 0.8 187.74
300000 0.7 171.72
400000 0.6 153.08
500000 0.5 131.15
600000 0.4 119.57
700000 0.3 111.20
800000 0.2 105.31
900000 0.1 103.07
1000000 0.0 102.89
 
0 -> 1 -> 2 -> 12 -> 11 -> 32 -> 33 -> 43 ->
42 -> 52 -> 51 -> 41 -> 31 -> 30 -> 40 -> 50 ->
60 -> 61 -> 62 -> 63 -> 53 -> 54 -> 44 -> 34 ->
24 -> 25 -> 14 -> 15 -> 16 -> 26 -> 36 -> 35 ->
45 -> 55 -> 56 -> 46 -> 47 -> 57 -> 58 -> 68 ->
67 -> 77 -> 86 -> 76 -> 66 -> 65 -> 64 -> 74 ->
75 -> 85 -> 84 -> 83 -> 73 -> 72 -> 71 -> 70 ->
80 -> 90 -> 91 -> 81 -> 82 -> 92 -> 93 -> 94 ->
95 -> 96 -> 97 -> 87 -> 98 -> 99 -> 89 -> 88 ->
78 -> 79 -> 69 -> 59 -> 49 -> 48 -> 39 -> 38 ->
37 -> 27 -> 17 -> 18 -> 28 -> 29 -> 19 -> 9 ->
8 -> 7 -> 6 -> 5 -> 4 -> 3 -> 13 -> 23 ->
22 -> 21 -> 20 -> 10 ->
 
Final E(s): 102.89
</pre>
 
 
 
=={{header|Scheme}}==
Line 1,585 ⟶ 3,394:
 
 
<langsyntaxhighlight lang="scheme">(cond-expand
(r7rs)
(chicken (import r7rs)))
Line 1,786 ⟶ 3,595:
(display (path-length s-final))
(newline)))
(newline)</langsyntaxhighlight>
 
{{out}}
Line 1,837 ⟶ 3,646:
 
 
<langsyntaxhighlight lang="scheme">(cond-expand
(r7rs)
(chicken (import r7rs)))
Line 2,044 ⟶ 3,853:
(format #t "Final E(s): ~,5F~%" (E_s s-final))
(format #t "Final path length: ~,5F~%" (path-length s-final))
(newline)</langsyntaxhighlight>
 
 
Line 2,167 ⟶ 3,976:
=={{header|Sidef}}==
{{trans|Julia}}
<langsyntaxhighlight lang="ruby">module TravelingSalesman {
 
# Eₛ: length(path)
Line 2,259 ⟶ 4,068:
}.map(1..100)
 
TravelingSalesman::findpath(citydist, 1e6, 1)</langsyntaxhighlight>
 
{{out}}
Line 2,307 ⟶ 4,116:
{{libheader|Wren-math}}
{{libheader|Wren-fmt}}
<langsyntaxhighlight ecmascriptlang="wren">import "random" for Random
import "./math" for Math
import "./fmt" for Fmt
 
// distances
Line 2,405 ⟶ 4,214:
}
 
sa.call(1e6, 1)</langsyntaxhighlight>
 
{{out}}
Line 2,442 ⟶ 4,251:
=={{header|zkl}}==
{{trans|EchoLisp}}
<langsyntaxhighlight lang="zkl">var [const] _dists=(0d10_000).pump(List,fcn(abcd){ // two points (a,b) & (c,d), calc distance
ab,cd,a,b,c,d:=abcd/100, abcd%100, ab/10,ab%10, cd/10,cd%10;
(a-c).toFloat().hypot(b-d)
Line 2,500 ⟶ 4,309:
println("E(s_final) %f".fmt(Emin));
println("Path: ",s.toString(*));
}</langsyntaxhighlight>
<syntaxhighlight lang ="zkl">sa(0d1_000_000,1);</langsyntaxhighlight>
{{out}}
<pre>
Peak
2,442

edits

Retrieved from "https://rosettacode.org/wiki/Special:MobileDiff/266339...364114"