Simulated annealing

Quoted from the Wikipedia page : Simulated annealing (SA) is a probabilistic technique for approximating the global optimum of a given function. Simulated annealing interprets slow cooling as a slow decrease in the probability of temporarily accepting worse solutions as it explores the solution space.

Pseudo code from Wikipedia

Notations :
  T : temperature. Decreases to 0.
  s : a system state
  E(s) : Energy at s. The function we want to minimize
  ∆E : variation of E, from state s to state s_next
  P(∆E , T) : Probability to move from s to s_next. 
  	if  ( ∆E < 0 ) P = 1
  	      else P = exp ( - ∆E / T) . Decreases as T →  0
  
Pseudo-code:
    Let s = s0  -- initial state
    For k = 0 through kmax (exclusive):
        T ← temperature(k , kmax)
        Pick a random neighbour state , s_next ← neighbour(s)
        ∆E ← E(s) - E(s_next) 
        If P(∆E , T) ≥ random(0, 1), move to the new state:
            s ← s_next
    Output: the final state s

Problem statement

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 cities. 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).

Definition : The neighbours of a city are the closest cities at distance 1 horizontally/vertically, or √2 diagonally. A corner city (0,9,90,99) has 3 neighbours. A center city has 8 neighbours.

Distances between cities
d ( 0, 7) → 7
d ( 0, 99) → 12.7279
d ( 23, 78) → 7.0711
d ( 33, 44) → 1.4142 // sqrt(2)

Task

Apply SA to the travelling salesman problem, using the following set of parameters/functions :

kT = 1 (Multiplication by kT is a placeholder, representing computing temperature as a function of 1-k/kmax):
temperature (k, kmax) = kT * (1 - k/kmax)
neighbour (s) : Pick a random city u > 0 . Pick a random neighbour city v > 0 of u , among u's 8 (max) neighbours on the grid. Swap u and v in s . This gives the new state s_next.
kmax = 1000_000
s0 = a random permutation

For k = 0 to kmax by step kmax/10 , display k, T, E(s). Display the final state s_final, and E(s_final).

You will see that the Energy may grow to a local optimum, before decreasing to a global optimum.

Illustrated example Temperature charts

Numerical example

kT = 1
E(s0) = 529.9158

k:  0         T:  1       Es:  529.9158
k:  100000    T:  0.9     Es:  201.1726
k:  200000    T:  0.8     Es:  178.1723
k:  300000    T:  0.7     Es:  154.7069
k:  400000    T:  0.6     Es:  158.1412 <== local optimum
k:  500000    T:  0.5     Es:  133.856
k:  600000    T:  0.4     Es:  129.5684
k:  700000    T:  0.3     Es:  112.6919
k:  800000    T:  0.2     Es:  105.799
k:  900000    T:  0.1     Es:  102.8284
k:  1000000   T:  0       Es:  102.2426

E(s_final) =    102.2426    
Path  s_final =   ( 0 10 11 21 31 20 30 40 50 60 70 80 90 91 81 71 73 83 84 74 64 54 55 65 75 76 66
 67 77 78 68 58 48 47 57 56 46 36 37 27 26 16 15 5 6 7 17 18 8 9 19 29 28 38 39 49 59 69 
79 89 99 98 88 87 97 96 86 85 95 94 93 92 82 72 62 61 51 41 42 52 63 53 43 32 22 12 13 
23 33 34 44 45 35 25 24 14 4 3 2 1 0)

Extra credit

Tune the parameters kT, kmax, or use different temperature() and/or neighbour() functions to demonstrate a quicker convergence, or a better optimum.

Ada

Translation of: C

Translation of: Scheme

Works with: GNAT version 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.

<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;

</lang>

Output:

An example run. In the following, you could use gnatmake instead of gprbuild.

$ 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

C

Translation of: Scheme

For your platform you might have to modify parts, such as the call to getentropy(3).

You can 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.

<lang c>#include <math.h>

include <stdio.h>
include <stdlib.h>
include <stdint.h>
include <string.h>
include <unistd.h>

define VECSZ 100
define STATESZ 64

typedef float floating_pt;

define EXP expf
define SQRT sqrtf

static floating_pt randnum (void) {

 return (floating_pt)
   ((double) (random () & 2147483647) / 2147483648.0);

}

static void shuffle (uint8_t vec[], size_t i, size_t n) {

 /* A Fisher-Yates shuffle of n elements of vec, starting at index
    i. */
 for (size_t j = 0; j != n; j += 1)
   {
     size_t k = i + j + (random () % (n - j));
     uint8_t xi = vec[i];
     uint8_t xk = vec[k];
     vec[i] = xk;
     vec[k] = xi;
   }

}

static void init_s (uint8_t vec[VECSZ]) {

 for (uint8_t j = 0; j != VECSZ; j += 1)
   vec[j] = j;
 shuffle (vec, 1, VECSZ - 1);

}

static inline void add_neighbor (uint8_t neigh[8],

             unsigned int *neigh_size,
             uint8_t neighbor)

{

 if (neighbor != 0)
   {
     neigh[*neigh_size] = neighbor;
     *neigh_size += 1;
   }

}

static void neighborhood (uint8_t neigh[8],

             unsigned int *neigh_size,
             uint8_t city)

{

 /* Find all non-zero neighbor cities. */

 const uint8_t i = city / 10;
 const uint8_t j = city % 10;

 uint8_t c0 = 0;
 uint8_t c1 = 0;
 uint8_t c2 = 0;
 uint8_t c3 = 0;
 uint8_t c4 = 0;
 uint8_t c5 = 0;
 uint8_t c6 = 0;
 uint8_t c7 = 0;

 if (i < 9)
   {
     c0 = (10 * (i + 1)) + j;
     if (j < 9)
       c1 = (10 * (i + 1)) + (j + 1);
     if (0 < j)
       c2 = (10 * (i + 1)) + (j - 1);
   }
 if (0 < i)
   {
     c3 = (10 * (i - 1)) + j;
     if (j < 9)
       c4 = (10 * (i - 1)) + (j + 1);
     if (0 < j)
       c5 = (10 * (i - 1)) + (j - 1);
   }
 if (j < 9)
   c6 = (10 * i) + (j + 1);
 if (0 < j)
   c7 = (10 * i) + (j - 1);

 *neigh_size = 0;
 add_neighbor (neigh, neigh_size, c0);
 add_neighbor (neigh, neigh_size, c1);
 add_neighbor (neigh, neigh_size, c2);
 add_neighbor (neigh, neigh_size, c3);
 add_neighbor (neigh, neigh_size, c4);
 add_neighbor (neigh, neigh_size, c5);
 add_neighbor (neigh, neigh_size, c6);
 add_neighbor (neigh, neigh_size, c7);

}

static floating_pt distance (uint8_t m, uint8_t n) {

 const uint8_t im = m / 10;
 const uint8_t jm = m % 10;
 const uint8_t in = n / 10;
 const uint8_t jn = n % 10;
 const int di = (int) im - (int) in;
 const int dj = (int) jm - (int) jn;
 return SQRT ((di * di) + (dj * dj));

}

static floating_pt path_length (uint8_t vec[VECSZ]) {

 floating_pt len = distance (vec[0], vec[VECSZ - 1]);
 for (size_t j = 0; j != VECSZ - 1; j += 1)
   len += distance (vec[j], vec[j + 1]);
 return len;

}

static void swap_s_elements (uint8_t vec[], uint8_t u, uint8_t v) {

 size_t j = 1;
 size_t iu = 0;
 size_t iv = 0;
 while (iu == 0 || iv == 0)
   {
     if (vec[j] == u)
       iu = j;
     else if (vec[j] == v)
       iv = j;
     j += 1;
   }
 vec[iu] = v;
 vec[iv] = u;

}

static void update_s (uint8_t vec[]) {

 const uint8_t u = 1 + (random () % (VECSZ - 1));
 uint8_t neighbors[8];
 unsigned int num_neighbors;
 neighborhood (neighbors, &num_neighbors, u);
 const uint8_t v = neighbors[random () % num_neighbors];
 swap_s_elements (vec, u, v);

}

static inline void copy_s (uint8_t dst[VECSZ], uint8_t src[VECSZ]) {

 memcpy (dst, src, VECSZ * (sizeof src[0]));

}

static void trial_s (uint8_t trial[VECSZ], uint8_t vec[VECSZ]) {

 copy_s (trial, vec);
 update_s (trial);

}

static floating_pt temperature (floating_pt kT, unsigned int kmax, unsigned int k) {

 return kT * (1 - ((floating_pt) k / (floating_pt) kmax));

}

static floating_pt probability (floating_pt delta_E, floating_pt T) {

 floating_pt prob;
 if (delta_E < 0)
   prob = 1;
 else if (T == 0)
   prob = 0;
 else
   prob = EXP (-(delta_E / T));
 return prob;

}

static void show (unsigned int k, floating_pt T, floating_pt E) {

 printf (" %7u %7.1f %13.5f\n", k, (double) T, (double) E);

}

static void simulate_annealing (floating_pt kT,

                   unsigned int kmax,
                   uint8_t s[VECSZ])

{

 uint8_t trial[VECSZ];

 unsigned int kshow = kmax / 10;
 floating_pt E = path_length (s);
 for (unsigned int k = 0; k != kmax + 1; k += 1)
   {
     const floating_pt T = temperature (kT, kmax, k);
     if (k % kshow == 0)
       show (k, T, E);
     trial_s (trial, s);
     const floating_pt E_trial = path_length (trial);
     const floating_pt delta_E = E_trial - E;
     const floating_pt P = probability (delta_E, T);
     if (P == 1 || randnum () <= P)
       {
         copy_s (s, trial);
         E = E_trial;
       }
   }

}

static void display_path (uint8_t vec[VECSZ]) {

 for (size_t i = 0; i != VECSZ; i += 1)
   {
     const uint8_t x = vec[i];
     printf ("%2u ->", (unsigned int) x);
     if ((i % 8) == 7)
       printf ("\n");
     else
       printf (" ");
   }
 printf ("%2u\n", vec[0]);

}

int main (void) {

 char state[STATESZ];
 uint32_t seed[1];
 int status = getentropy (&seed[0], sizeof seed[0]);
 if (status < 0)
   seed[0] = 1;
 initstate (seed[0], state, STATESZ);

 floating_pt kT = 1.0;
 unsigned int kmax = 1000000;

 uint8_t s[VECSZ];
 init_s (s);

 printf ("\n");
 printf ("   kT: %f\n", (double) kT);
 printf ("   kmax: %u\n", kmax);
 printf ("\n");
 printf ("       k       T          E(s)\n");
 printf (" -----------------------------\n");
 simulate_annealing (kT, kmax, s);
 printf ("\n");
 display_path (s);
 printf ("\n");
 printf ("Final E(s): %.5f\n", (double) path_length (s));
 printf ("\n");

 return 0;

}</lang>

Output:

An example run:

$ cc -Ofast -march=native simanneal.c -lm && ./a.out

   kT: 1.000000
   kmax: 1000000

       k       T          E(s)
 -----------------------------
       0     1.0     383.25223
  100000     0.9     195.81190
  200000     0.8     186.58963
  300000     0.7     152.46564
  400000     0.6     143.59039
  500000     0.5     130.91815
  600000     0.4     126.53572
  700000     0.3     112.85691
  800000     0.2     103.72134
  900000     0.1     103.07108
 1000000     0.0     102.24265

 0 -> 10 -> 20 -> 21 -> 31 -> 30 -> 40 -> 50 ->
60 -> 61 -> 62 -> 72 -> 82 -> 81 -> 71 -> 70 ->
80 -> 90 -> 91 -> 92 -> 93 -> 94 -> 84 -> 83 ->
73 -> 63 -> 64 -> 65 -> 66 -> 76 -> 86 -> 87 ->
77 -> 67 -> 68 -> 58 -> 57 -> 56 -> 55 -> 45 ->
35 -> 26 -> 36 -> 46 -> 47 -> 48 -> 38 -> 37 ->
27 -> 28 -> 29 -> 39 -> 49 -> 59 -> 69 -> 79 ->
78 -> 88 -> 89 -> 99 -> 98 -> 97 -> 96 -> 95 ->
85 -> 75 -> 74 -> 54 -> 53 -> 52 -> 51 -> 41 ->
42 -> 43 -> 44 -> 34 -> 33 -> 32 -> 22 -> 23 ->
14 ->  4 ->  5 ->  6 ->  7 ->  8 ->  9 -> 19 ->
18 -> 17 -> 16 -> 15 -> 25 -> 24 -> 13 ->  3 ->
 2 -> 12 -> 11 ->  1 ->  0

Final E(s): 102.24265

C++

Compiler: MSVC (19.27.29111 for x64) <lang c++>

include<array>
include<utility>
include<cmath>
include<random>
include<iostream>

using coord = std::pair<int,int>; constexpr size_t numCities = 100;

// CityID with member functions to get position struct CityID{

   int v{-1};
   CityID() = default;
   constexpr explicit CityID(int i) noexcept : v(i){}
   constexpr explicit CityID(coord ij) : v(ij.first * 10 + ij.second) {
       if(ij.first < 0 || ij.first > 9 || ij.second < 0 || ij.second > 9){
           throw std::logic_error("Cannot construct CityID from invalid coordinates!");
       }
   }

   constexpr coord get_pos() const noexcept { return {v/10,v%10}; }

}; bool operator==(CityID const& lhs, CityID const& rhs) {return lhs.v == rhs.v;}

// Function for distance between two cities double dist(coord city1, coord city2){

   double diffx = city1.first - city2.first;
   double diffy = city1.second - city2.second;
   return std::sqrt(std::pow(diffx, 2) + std::pow(diffy,2));

}

// Function for total distance travelled template<size_t N> double dist(std::array<CityID,N> cities){

   double sum = 0;
   for(auto it = cities.begin(); it < cities.end() - 1; ++it){
       sum += dist(it->get_pos(),(it+1)->get_pos());
   }
   sum += dist((cities.end()-1)->get_pos(), cities.begin()->get_pos());
   return sum;

}

// 8 nearest cities, id cannot be at the border and has to have 8 valid neighbors constexpr std::array<CityID,8> get_nearest(CityID id){

   auto const ij = id.get_pos();
   auto const i = ij.first;
   auto const j = ij.second;
   return {
       CityID({i-1,j-1}),
       CityID({i  ,j-1}),
       CityID({i+1,j-1}),
       CityID({i-1,j  }),
       CityID({i+1,j  }),
       CityID({i-1,j+1}),
       CityID({i  ,j+1}),
       CityID({i+1,j+1}),
   };

}

// Function for formating of results constexpr int get_num_digits(int num){

   int digits = 1;
   while(num /= 10){
       ++digits;
   }
   return digits;

}

// Function for shuffeling of initial state template<typename It, typename RandomEngine> void shuffle(It first, It last, RandomEngine& rand_eng){

   for(auto i=(last-first)-1; i>0; --i){
       std::uniform_int_distribution<int> dist(0,i);
       std::swap(first[i], first[dist(rand_eng)]);
   }

}

class SA{

   int kT{1};
   int kmax{1'000'000};
   std::array<CityID,numCities> s;
   std::default_random_engine rand_engine{0};

   // Temperature 
   double temperature(int k) const { return kT * (1.0 - static_cast<double>(k) / kmax); }

   // Probabilty of acceptance between 0.0 and 1.0
   double P(double dE, double T){
       if(dE < 0){
           return 1;
       }
       else{
           return std::exp(-dE/T);
       }
   }

   // Permutation of state through swapping of cities in travel path
   std::array<CityID,numCities> next_permut(std::array<CityID,numCities> cities){
       std::uniform_int_distribution<> disx(1,8);
       std::uniform_int_distribution<> disy(1,8);
       auto randCity = CityID({disx(rand_engine),disy(rand_engine)});      // Select city which is not at the border, since all neighbors are valid under this condition and all permutations are still possible
       auto neighbors = get_nearest(randCity);                             // Get list of nearest neighbors
       std::uniform_int_distribution<> selector(0,neighbors.size()-1);     // [0,7]
       const auto [i,j] = randCity.get_pos();
       auto randNeighbor = neighbors[selector(rand_engine)];               // Since randCity is not at the border, all 8 neighbors are valid
       auto cityit1 = std::find(cities.begin(),cities.end(),randCity);     // Find selected city in state
       auto cityit2 = std::find(cities.begin(), cities.end(),randNeighbor);// Find selected neighbor in state
       std::swap(*cityit1, *cityit2);                                      // Swap city and neighbor
       return cities;
   }

   // Logging function for progress output
   void log_progress(int k, double T, double E) const {
       auto nk = get_num_digits(kmax);
       auto nt = get_num_digits(kT);
       std::printf("k: %*i | T: %*.3f | E(s): %*.4f\n", nk, k, nt, T, 3, E);
   }

public:

   // Initialize state with integers from 0 to 99
   SA() {
       int i = 0;
       for(auto it = s.begin(); it != s.end(); ++it){
           *it = CityID(i);
           ++i;
       }
       shuffle(s.begin(),s.end(),rand_engine);
   }

   // Logging function for final path
   void log_path(){
       for(size_t idx = 0; idx < s.size(); ++idx){
           std::printf("%*i -> ", 2, s[idx].v);
           if((idx + 1)%20 == 0){
               std::printf("\n");
           }
       }
       std::printf("%*i", 2, s[0].v);
   }

   // Core simulated annealing algorithm
   std::array<CityID,numCities> run(){
       std::cout << "kT == " << kT << "\n" << "kmax == " << kmax << "\n" << "E(s0) == " << dist(s) << "\n";
       for(int k = 0; k < kmax; ++k){
           auto T = temperature(k);
           auto const E1 = dist(s);
           auto s_next{next_permut(s)};
           auto const E2 = dist(s_next);
           auto const dE = E2 - E1; // lower is better
           std::uniform_real_distribution dist(0.0, 1.0);
           auto E = E1;
           if(P(dE,T) >= dist(rand_engine)){
               s = s_next;
               E = E2;
           }
           if(k%100000 == 0){
               log_progress(k,T,E1);
           }
       }
       log_progress(kmax,0.0,dist(s));
       std::cout << "\nFinal path: \n";
       log_path();
       return s;
   }

};

int main(){

   SA sa{};
   auto result = sa.run(); // Run simulated annealing and log progress and result
   std::cin.get();
   return 0;

} </lang>

Output:

kT == 1
kmax == 1000000
E(s0) == 529.423
k:       0 | T: 1.000 | E(s): 529.4231
k:  100000 | T: 0.900 | E(s): 197.5111
k:  200000 | T: 0.800 | E(s): 183.7467
k:  300000 | T: 0.700 | E(s): 165.8442
k:  400000 | T: 0.600 | E(s): 143.8588
k:  500000 | T: 0.500 | E(s): 133.9247
k:  600000 | T: 0.400 | E(s): 125.9499
k:  700000 | T: 0.300 | E(s): 115.8657
k:  800000 | T: 0.200 | E(s): 107.8635
k:  900000 | T: 0.100 | E(s): 102.4853
k: 1000000 | T: 0.000 | E(s): 102.4853

Final path:
71 -> 61 -> 51 -> 50 -> 60 -> 70 -> 80 -> 90 -> 91 -> 92 -> 82 -> 83 -> 93 -> 94 -> 84 -> 85 -> 95 -> 96 -> 86 -> 76 ->
75 -> 74 -> 64 -> 65 -> 55 -> 45 -> 44 -> 54 -> 53 -> 43 -> 33 -> 34 -> 35 -> 26 -> 16 ->  6 ->  7 -> 17 -> 27 -> 37 ->
47 -> 57 -> 58 -> 48 -> 38 -> 28 -> 18 ->  8 ->  9 -> 19 -> 29 -> 39 -> 49 -> 59 -> 69 -> 68 -> 78 -> 79 -> 88 -> 89 ->
99 -> 98 -> 97 -> 87 -> 77 -> 67 -> 66 -> 56 -> 46 -> 36 -> 25 -> 24 -> 14 -> 15 ->  5 ->  4 ->  3 ->  2 -> 11 -> 21 ->
31 -> 41 -> 40 -> 30 -> 20 -> 10 ->  0 ->  1 -> 12 -> 13 -> 23 -> 22 -> 32 -> 42 -> 52 -> 62 -> 63 -> 73 -> 72 -> 81 ->
71

EchoLisp

<lang scheme> (lib 'math)

distances

(define (d ci cj) (distance (% ci 10) (quotient ci 10) (% cj 10) (quotient cj 10))) (define _dists (build-vector 10000 (lambda (ij) (d (quotient ij 100) (% ij 100))))) (define-syntax-rule (dist ci cj) [_dists (+ ci (* 100 cj))])

E(s) = length(path)

(define (Es path) (define lpath (vector->list path)) (for/sum ((ci lpath) (cj (rest lpath))) (dist ci cj)))

temperature() function

(define (T k kmax kT) (* kT (- 1 (// k kmax))))

|

alternative temperature() must be decreasing with k increasing and → 0

(define (T k kmax kT) (* kT (- 1 (sin (* PI/2 (// k kmax)))))) |#

∆E = Es_new - Es_old > 0 probability to move if ∆E > 0, → 0 when T → 0 (frozen state)

(define (P ∆E k kmax kT) (exp (// (- ∆E ) (T k kmax kT))))

∆E from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..) ∆E before swapping (u,v) Quicker than Es(s_next) - Es(s)

(define (dE s u v)

old

(define a (dist [s (1- u)] [s u])) (define b (dist [s (1+ u)] [s u])) (define c (dist [s (1- v)] [s v])) (define d (dist [s (1+ v)] [s v]))

new

(define na (dist [s (1- u)] [s v])) (define nb (dist [s (1+ u)] [s v])) (define nc (dist [s (1- v)] [s u])) (define nd (dist [s (1+ v)] [s u]))

(cond ((= v (1+ u)) (- (+ na nd) (+ a d))) ((= u (1+ v)) (- (+ nc nb) (+ c b))) (else (- (+ na nb nc nd) (+ a b c d)))))

all 8 neighbours

(define dirs #(1 -1 10 -10 9 11 -11 -9))

(define (sa kmax (kT 10)) (define s (list->vector (cons 0 (append (shuffle (range 1 100)) 0)))) (printf "E(s0) %d" (Es s)) ;; random starter (define Emin (Es s)) ;; E0

(for ((k kmax)) (when (zero? (% k (/ kmax 10))) (printf "k: %10d T: %8.4d Es: %8.4d" k (T k kmax kT) (Es s)) )

(define u (1+ (random 99))) ;; city index 1 99 (define cv (+ [s u] [dirs (random 8)])) ;; city number #:continue (or (> cv 99) (<= cv 0)) #:continue (> (dist [s u] cv) 5) ;; check true neighbour (eg 0 9) (define v (vector-index cv s 1)) ;; city index

(define ∆e (dE s u v)) (when (or (< ∆e 0) ;; always move if negative (>= (P ∆e k kmax kT) (random))) (vector-swap! s u v) (+= Emin ∆e))

;; (assert (= (round Emin) (round (Es s)))) ) ;; for

(printf "k: %10d T: %8.4d Es: %8.4d" kmax (T (1- kmax) kmax kT) (Es s)) (s-plot s 0) (printf "E(s_final) %d" Emin) (writeln 'Path s)) </lang>

Output:

(sa 1000000 1)

E(s0) 501.0909

k:  0         T:  1       Es:  501.0909
k:  100000    T:  0.9     Es:  167.3632
k:  200000    T:  0.8     Es:  160.7791
k:  300000    T:  0.7     Es:  166.8746
k:  400000    T:  0.6     Es:  142.579
k:  500000    T:  0.5     Es:  131.0657
k:  600000    T:  0.4     Es:  116.9214
k:  700000    T:  0.3     Es:  110.8569
k:  800000    T:  0.2     Es:  103.3137
k:  900000    T:  0.1     Es:  102.4853
k:  1000000   T:  0       Es:  102.4853

E(s_final)     102.4853    
Path     #( 0 10 20 30 40 50 60 70 71 61 62 53 63 64 54 44 45 55 65
 74 84 83 73 72 82 81 80 90 91 92 93 94 95 85 75 76 86 96 97 98 99
 88 89 79 69 59 49 48 47 57 58 68 78 87 77 67 66 56 46 36 35 25 24
 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)

Go

Translation of: zkl

<lang go>package main

import (

   "fmt"
   "math"
   "math/rand"
   "time"

)

var (

   dists = calcDists()
   dirs  = [8]int{1, -1, 10, -10, 9, 11, -11, -9} // all 8 neighbors

)

// distances func calcDists() []float64 {

   dists := make([]float64, 10000)
   for i := 0; i < 10000; i++ {
       ab, cd := math.Floor(float64(i)/100), float64(i%100)
       a, b := math.Floor(ab/10), float64(int(ab)%10)
       c, d := math.Floor(cd/10), float64(int(cd)%10)
       dists[i] = math.Hypot(a-c, b-d)
   }
   return dists

}

// index into lookup table of float64s func dist(ci, cj int) float64 {

   return dists[cj*100+ci]

}

// energy at s, to be minimized func Es(path []int) float64 {

   d := 0.0
   for i := 0; i < len(path)-1; i++ {
       d += dist(path[i], path[i+1])
   }
   return d

}

// temperature function, decreases to 0 func T(k, kmax, kT int) float64 {

   return (1 - float64(k)/float64(kmax)) * float64(kT)

}

// variation of E, from state s to state s_next func dE(s []int, u, v int) float64 {

   su, sv := s[u], s[v]
   // old
   a, b, c, d := dist(s[u-1], su), dist(s[u+1], su), dist(s[v-1], sv), dist(s[v+1], sv)
   // new
   na, nb, nc, nd := dist(s[u-1], sv), dist(s[u+1], sv), dist(s[v-1], su), dist(s[v+1], su)
   if v == u+1 {
       return (na + nd) - (a + d)
   } else if u == v+1 {
       return (nc + nb) - (c + b)
   } else {
       return (na + nb + nc + nd) - (a + b + c + d)
   }

}

// probability to move from s to s_next func P(deltaE float64, k, kmax, kT int) float64 {

   return math.Exp(-deltaE / T(k, kmax, kT))

}

func sa(kmax, kT int) {

   rand.Seed(time.Now().UnixNano())
   temp := make([]int, 99)
   for i := 0; i < 99; i++ {
       temp[i] = i + 1
   }
   rand.Shuffle(len(temp), func(i, j int) {
       temp[i], temp[j] = temp[j], temp[i]
   })
   s := make([]int, 101) // all 0 by default
   copy(s[1:], temp)     // random path from 0 to 0
   fmt.Println("kT =", kT)
   fmt.Printf("E(s0) %f\n\n", Es(s)) // random starter
   Emin := Es(s)                     // E0
   for k := 0; k <= kmax; k++ {
       if k%(kmax/10) == 0 {
           fmt.Printf("k:%10d   T: %8.4f   Es: %8.4f\n", k, T(k, kmax, kT), Es(s))
       }
       u := 1 + rand.Intn(99)          // city index 1 to 99
       cv := s[u] + dirs[rand.Intn(8)] // city number
       if cv <= 0 || cv >= 100 {       // bogus city
           continue
       }
       if dist(s[u], cv) > 5 { // check true neighbor (eg 0 9)
           continue
       }
       v := s[cv] // city index
       deltae := dE(s, u, v)
       if deltae < 0 || // always move if negative
           P(deltae, k, kmax, kT) >= rand.Float64() {
           s[u], s[v] = s[v], s[u]
           Emin += deltae
       }
   }
   fmt.Printf("\nE(s_final) %f\n", Emin)
   fmt.Println("Path:")
   // output final state
   for i := 0; i < len(s); i++ {
       if i > 0 && i%10 == 0 {
           fmt.Println()
       }
       fmt.Printf("%4d", s[i])
   }
   fmt.Println()

}

func main() {

   sa(1e6, 1)

}</lang>

Output:

Sample run:

kT = 1
E(s0) 520.932463

k:         0   T:   1.0000   Es: 520.9325
k:    100000   T:   0.9000   Es: 185.1279
k:    200000   T:   0.8000   Es: 167.7657
k:    300000   T:   0.7000   Es: 158.6923
k:    400000   T:   0.6000   Es: 151.6564
k:    500000   T:   0.5000   Es: 139.9185
k:    600000   T:   0.4000   Es: 132.9964
k:    700000   T:   0.3000   Es: 121.8962
k:    800000   T:   0.2000   Es: 120.0445
k:    900000   T:   0.1000   Es: 116.8476
k:   1000000   T:   0.0000   Es: 116.5565

E(s_final) 116.556509
Path:
   0  11  21  31  41  51  52  61  62  72
  82  73  74  64  44  45  55  54  63  53
  42  32  43  33  35  34  24  23  22  13
  12   2   3   4  14  25  26   7   6  16
  15   5  17  27  36  46  56  66  65  75
  77  78  68  69  59  49  39  38  37  28
  29  19   9   8  18  47  48  58  57  67
  76  86  85  95  96  97  87  88  79  89
  99  98  84  94  83  93  92  91  90  80
  81  71  70  60  50  40  30  20  10   1
   0

J

Implementation:

<lang J>dist=: +/&.:*:@:-"1/~10 10#:i.100

satsp=:4 :0

 kT=. 1
 pathcost=. [: +/ 2 {&y@<\ 0 , ] , 0:
 neighbors=. 0 (0}"1) y e. 1 2{/:~~.,y
 s=. (?~#y)-.0
 d=. pathcost s
 step=. x%10
 for_k. i.x+1 do.
   T=. kT*1-k%x
   u=. ({~ ?@#)s
   v=. ({~ ?@#)I.u{neighbors
   sk=. (?0 do.
     s=.sk
     d=.dk
   end.
   if. 0=step|k do.
     echo k,T,d
   end.
 end.
 0,s,0

)</lang>
Notes:
E(s_final) gets displayed on the kmax progress line.
We do not do anything special for negative deltaE because the exponential will be greater than 1 for that case and that will always be greater than our random number from the range 0..1.
Also, while we leave connection distances (and, thus, number of cities) as a parameter, some other aspects of this problem made more sense when included in the implementation:
We leave city 0 out of our data structure, since it can't appear in the middle of our path. But we bring it back in when computing path distance.
Neighbors are any city which have one of the two closest non-zero distances from the current city (and specifically excluding city 0, since that is anchored as our start and end city).
Sample run:
<lang J> 1e6 satsp dist 0 1 538.409 100000 0.9 174.525 200000 0.8 165.541 300000 0.7 173.348 400000 0.6 168.188 500000 0.5 134.983 600000 0.4 121.585 700000 0.3 111.443 800000 0.2 101.657 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</lang>
~~Julia~~
Translation of: EchoLisp
Module: <lang julia>module TravelingSalesman
using Random, Printf
~~Eₛ: length(path)~~
~~Eₛ(distances, path) = sum(distances[ci, cj] for (ci, cj) in zip(path, Iterators.drop(path, 1)))~~
~~T: temperature~~
~~T(k, kmax, kT) = kT * (1 - k / kmax)~~
Alternative temperature:

T(k, kmax, kT) = kT * (1 - sin(π / 2 * k / kmax))
ΔE = Eₛ_new - Eₛ_old > 0

Prob. to move if ΔE > 0, → 0 when T → 0 (fronzen state)
~~P(ΔE, k, kmax, kT) = exp(-ΔE / T(k, kmax, kT))~~
∆E from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..)

∆E before swapping (u,v)

Quicker than Eₛ(s_next) - Eₛ(path)
~~function dE(distances, path, u, v)~~
~~a = distances[path[u - 1], path[u]] b = distances[path[u + 1], path[u]] c = distances[path[v - 1], path[v]] d = distances[path[v + 1], path[v]]~~
~~na = distances[path[u - 1], path[v]] nb = distances[path[u + 1], path[v]] nc = distances[path[v - 1], path[u]] nd = distances[path[v + 1], path[u]]~~
~~if v == u + 1 return (na + nd) - (a + d) elseif u == v + 1 return (nc + nb) - (c + b) else return (na + nb + nc + nd) - (a + b + c + d) end~~
end
const dirs = [1, -1, 10, -10, 9, 11, -11, -9]
function _prettypath(path)
~~r = IOBuffer() for g in Iterators.partition(path, 10) println(r, join(lpad.(g, 3), ", ")) end return String(take!(r))~~
end
function findpath(distances, kmax, kT)
~~n = size(distances, 1) path = vcat(1, shuffle(2:n), 1) Emin = Eₛ(distances, path) @printf("\n# Entropy(s₀) = %10.2f\n", Emin) println("# Random path: \n", _prettypath(path))~~
for k in Base.OneTo(kmax) if iszero(k % (kmax ÷ 10)) @printf("k: %10d | T: %8.4f | Eₛ: %8.4f\n", k, T(k, kmax, kT), Eₛ(distances, path)) end u = rand(2:n) v = path[u] + rand(dirs) v ∈ 2:n || continue
~~δE = dE(distances, path, u, v) if δE < 0 || P(δE, k, kmax, kT) ≥ rand() path[u], path[v] = path[v], path[u] Emin += δE end end~~
~~@printf("k: %10d | T: %8.4f | Eₛ: %8.4f\n", kmax, T(kmax, kmax, kT), Eₛ(distances, path)) println("\n# Found path:\n", _prettypath(path)) return path~~
end
end # module TravelingSalesman</lang>
Main: <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)</lang>

~~Output:~~
# Entropy(s₀) = 521.86 # Random path: 1, 2, 11, 80, 78, 73, 68, 19, 43, 69 86, 79, 66, 67, 77, 96, 26, 62, 60, 98 71, 3, 59, 37, 18, 40, 34, 92, 97, 6 84, 94, 29, 63, 36, 50, 87, 45, 83, 90 76, 28, 15, 38, 91, 58, 47, 44, 85, 17 25, 33, 31, 99, 27, 74, 53, 95, 16, 13 42, 88, 8, 4, 7, 64, 54, 9, 14, 41 5, 81, 65, 23, 75, 100, 89, 51, 20, 48 82, 12, 21, 55, 24, 70, 49, 10, 35, 72 52, 22, 61, 32, 46, 57, 30, 93, 39, 56 1 k: 100000 | T: 0.9000 | Eₛ: 184.4448 k: 200000 | T: 0.8000 | Eₛ: 175.3662 k: 300000 | T: 0.7000 | Eₛ: 169.0505 k: 400000 | T: 0.6000 | Eₛ: 160.8328 k: 500000 | T: 0.5000 | Eₛ: 147.1973 k: 600000 | T: 0.4000 | Eₛ: 132.9186 k: 700000 | T: 0.3000 | Eₛ: 126.9931 k: 800000 | T: 0.2000 | Eₛ: 122.0656 k: 900000 | T: 0.1000 | Eₛ: 119.7924 k: 1000000 | T: 0.0000 | Eₛ: 119.7924 k: 1000000 | T: 0.0000 | Eₛ: 119.7924 # Found path: 1, 2, 12, 13, 3, 4, 6, 7, 8, 9 19, 18, 17, 5, 14, 15, 16, 27, 28, 29 39, 38, 26, 25, 24, 23, 22, 10, 21, 20 30, 31, 32, 33, 34, 35, 36, 37, 49, 48 47, 46, 45, 44, 43, 42, 41, 40, 50, 51 52, 53, 54, 55, 56, 57, 58, 59, 69, 68 67, 65, 64, 63, 62, 61, 71, 60, 70, 80 81, 82, 72, 73, 74, 66, 78, 79, 89, 99 98, 97, 96, 95, 94, 85, 86, 87, 88, 77 76, 75, 84, 83, 93, 92, 91, 100, 90, 11 1
~~Nim~~
<lang Nim>import math, random, sequtils, strformat
const
~~kT = 1 kMax = 1_000_000~~
~~proc randomNeighbor(x: int): int =~~
case x of 0: sample([1, 10, 11]) of 9: sample([8, 18, 19]) of 90: sample([80, 81, 91]) of 99: sample([88, 89, 98]) elif x > 0 and x < 9: # top ceiling sample [x-1, x+1, x+9, x+10, x+11] elif x > 90 and x < 99: # bottom floor sample [x-11, x-10, x-9, x-1, x+1] elif x mod 10 == 0: # left wall sample([x-10, x-9, x+1, x+10, x+11]) elif (x+1) mod 10 == 0: # right wall sample([x-11, x-10, x-1, x+9, x+10]) else: # center sample([x-11, x-10, x-9, x-1, x+1, x+9, x+10, x+11])
~~proc neighbor(s: seq[int]): seq[int] =~~
result = s var city = sample(s) var cityNeighbor = city.randomNeighbor while cityNeighbor == 0 or city == 0: city = sample(s) cityNeighbor = city.randomNeighbor result[s.find city].swap result[s.find cityNeighbor]
~~func distNeighbor(a, b: int): float =~~
~~template divmod(a: int): (int, int) = (a div 10, a mod 10) let (diva, moda) = a.divmod (divb, modb) = b.divmod hypot((diva-divb).float, (moda-modb).float)~~
~~func temperature(k, kmax: float): float =~~
~~kT * (1 - (k / kmax))~~
~~func pdelta(eDelta, temp: float): float =~~
~~if eDelta < 0: 1.0 else: exp(-eDelta / temp)~~
~~func energy(path: seq[int]): float =~~
~~var sum = 0.distNeighbor path[0] for i in 1 ..< path.len: sum += path[i-1].distNeighbor(path[i]) sum + path[^1].distNeighbor 0~~
~~proc main =~~
randomize() var s = block: var x = toSeq(0..99) template shuffler: int = rand(1 .. x.high) for i in 1 .. x.high: x[i].swap x[shuffler()] x echo fmt"E(s0): {energy s:6.4f}" for k in 0 .. kMax: var temp = temperature(float k, float kMax) lastenergy = energy s newneighbor = s.neighbor newenergy = newneighbor.energy if k mod (kMax div 10) == 0: echo fmt"k: {k:7} T: {temp:6.2f} Es: {lastenergy:6.4f}" var deltaEnergy = newenergy - lastenergy if pDelta(deltaEnergy, temp) >= rand(1.0): s = newneighbor
~~s.add 0 echo fmt"E(sFinal): {energy s:6.4f}" echo fmt"path: {s}"~~
~~main()</lang>~~

~~Compile and run:~~

~~nim c -r -d:release --opt:speed travel_sa.nim~~

~~Output:~~
~~Sample run:~~
E(s0): 505.1591 k: 0 T: 1.00 Es: 505.1591 k: 100000 T: 0.90 Es: 196.5216 k: 200000 T: 0.80 Es: 165.6735 k: 300000 T: 0.70 Es: 159.3411 k: 400000 T: 0.60 Es: 144.8330 k: 500000 T: 0.50 Es: 131.7888 k: 600000 T: 0.40 Es: 127.6914 k: 700000 T: 0.30 Es: 113.9280 k: 800000 T: 0.20 Es: 104.7279 k: 900000 T: 0.10 Es: 103.3137 k: 1000000 T: 0.00 Es: 103.3137 E(sFinal): 103.3137 path: @[0, 10, 11, 22, 21, 20, 30, 31, 41, 40, 50, 51, 61, 60, 70, 71, 81, 80, 90, 91, 92, 93, 82, 83, 73, 72, 62, 63, 53, 52, 42, 32, 33, 23, 13, 14, 24, 34, 35, 25, 15, 16, 26, 36, 47, 48, 38, 39, 49, 59, 58, 57, 68, 69, 79, 89, 99, 98, 97, 96, 95, 94, 84, 74, 75, 85, 86, 87, 88, 78, 77, 67, 76, 66, 65, 64, 54, 43, 44, 45, 55, 56, 46, 37, 27, 28, 29, 19, 9, 8, 18, 17, 7, 6, 5, 4, 3, 2, 12, 1, 0]
~~Perl~~
Translation of: Raku
~~<lang perl>use utf8; use strict; use warnings; use List::Util qw(shuffle sum);~~
~~simulation setup~~
my $cities = 100; # number of cities my $kmax = 1e6; # iterations to run my $kT = 1; # initial 'temperature'
die 'City count must be a perfect square.' if sqrt($cities) != int sqrt($cities);
~~locations of (up to) 8 neighbors, with grid size derived from number of cities~~
~~my $gs = sqrt $cities; my @neighbors = (1, -1, $gs, -$gs, $gs-1, $gs+1, -($gs+1), -($gs-1));~~
~~matrix of distances between cities~~
~~my @D; for my $j (0 .. $cities**2 - 1) {~~
~~my ($ab, $cd) = (int($j/$cities), int($j%$cities)); my ($a, $b, $c, $d) = (int($ab/$gs), int($ab%$gs), int($cd/$gs), int($cd%$gs)); $D[$ab][$cd] = sqrt(($a - $c)**2 + ($b - $d)**2);~~
}
~~temperature function, decreases to 0~~
~~sub T {~~
~~my($k, $kmax, $kT) = @_; (1 - $k/$kmax) * $kT~~
}
~~probability to move from s to s_next~~
~~sub P {~~
~~my($ΔE, $k, $kmax, $kT) = @_; exp -$ΔE / T($k, $kmax, $kT)~~
}
~~energy at s, to be minimized~~
~~sub Es {~~
~~my(@path) = @_; sum map { $D[ $path[$_] ] [ $path[$_+1] ] } 0 .. @path-2~~
}
~~variation of E, from state s to state s_next~~
~~sub delta_E {~~
my($u, $v, @s) = @_; my ($a, $b, $c, $d) = ($D[$s[$u-1]][$s[$u]], $D[$s[$u+1]][$s[$u]], $D[$s[$v-1]][$s[$v]], $D[$s[$v+1]][$s[$v]]); my ($na, $nb, $nc, $nd) = ($D[$s[$u-1]][$s[$v]], $D[$s[$u+1]][$s[$v]], $D[$s[$v-1]][$s[$u]], $D[$s[$v+1]][$s[$u]]); if ($v == $u+1) { return ($na + $nd) - ($a + $d) } elsif ($u == $v+1) { return ($nc + $nb) - ($c + $b) } else { return ($na + $nb + $nc + $nd) - ($a + $b + $c + $d) }
}
~~E(s0), initial state~~
my @s = 0; map { push @s, $_ } shuffle 1..$cities-1; push @s, 0; my $E_min_global = my $E_min = Es(@s);
for my $k (0 .. $kmax-1) {
~~printf "k:%8u T:%4.1f Es: %3.1f\n" , $k, T($k, $kmax, $kT), Es(@s) if $k % ($kmax/10) == 0;~~
~~# valid candidate cities (exist, adjacent) my $u = 1 + int rand $cities-1; my $cv = $neighbors[int rand 8] + $s[$u]; next if $cv <= 0 or $cv >= $cities or $D[$s[$u]][$cv] > sqrt(2);~~
my $v = $s[$cv]; my $ΔE = delta_E($u, $v, @s); if ($ΔE < 0 or P($ΔE, $k, $kmax, $kT) >= rand) { # always move if negative ($s[$u], $s[$v]) = ($s[$v], $s[$u]); $E_min += $ΔE; $E_min_global = $E_min if $E_min < $E_min_global; }
}
printf "\nE(s_final): %.1f\n", $E_min_global; for my $l (0..4) {
~~printf "@{['%4d' x 20]}\n", @s[$l*20 .. ($l+1)*20 - 1];~~
~~} printf " 0\n";</lang>~~

~~Output:~~
k: 0 T: 1.0 Es: 519.2 k: 100000 T: 0.9 Es: 188.2 k: 200000 T: 0.8 Es: 178.5 k: 300000 T: 0.7 Es: 162.3 k: 400000 T: 0.6 Es: 157.0 k: 500000 T: 0.5 Es: 148.9 k: 600000 T: 0.4 Es: 128.7 k: 700000 T: 0.3 Es: 129.5 k: 800000 T: 0.2 Es: 119.8 k: 900000 T: 0.1 Es: 119.5 E(s_final): 119.1 0 12 23 24 35 36 26 27 16 7 8 9 19 29 28 18 17 6 14 13 22 32 33 34 25 15 5 4 3 2 1 11 20 21 31 30 40 51 50 60 61 62 53 43 44 54 56 57 48 49 39 38 37 46 45 55 65 64 63 74 84 83 82 81 80 90 91 92 93 94 95 85 66 47 58 59 69 89 88 87 77 67 68 78 79 99 98 97 96 86 76 75 73 72 70 71 52 42 41 10 0
~~Phix~~
Translation of: zkl
with javascript_semantics function hypot(atom a,b) return sqrt(a*a+b*b) end function function calc_dists() sequence dists = repeat(0,10000) for abcd=1 to 10000 do integer {ab,cd} = {floor(abcd/100),mod(abcd,100)}, {a,b,c,d} = {floor(ab/10),mod(ab,10), floor(cd/10),mod(cd,10)} dists[abcd] = hypot(a-c,b-d) end for return dists end function constant dists = calc_dists() function dist(integer ci,cj) return dists[cj*100+ci] end function function Es(sequence path) atom d = 0 for i=1 to length(path)-1 do d += dist(path[i],path[i+1]) end for return d end function -- temperature() function function T(integer k, kmax, kT) return (1-k/kmax)*kT end function -- deltaE = Es_new - Es_old > 0 -- probability to move if deltaE > 0, -->0 when T --> 0 (frozen state) function P(atom deltaE, integer k, kmax, kT) return exp(-deltaE/T(k,kmax,kT)) end function -- deltaE from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..) function dE(sequence s, integer u,v) -- (note that u,v are 0-based, but 1..99 here) -- integer sum1 = s[u-1], su = s[u], sup1 = s[u+1], -- svm1 = s[v-1], sv = s[v], svp1 = s[v+1] integer sum1 = s[u], su = s[u+1], sup1 = s[u+2], svm1 = s[v], sv = s[v+1], svp1 = s[v+2] -- old atom {a,b,c,d}:={dist(sum1,su), dist(su,sup1), dist(svm1,sv), dist(sv,svp1)}, -- new {na,nb,nc,nd}:={dist(sum1,sv), dist(sv,sup1), dist(svm1,su), dist(su,svp1)} return iff(v==u+1?(na+nd)-(a+d): iff(u==v+1?(nc+nb)-(c+b): (na+nb+nc+nd)-(a+b+c+d))) end function -- all 8 neighbours constant dirs = {1, -1, 10, -10, 9, 11, -11, -9} procedure sa(integer kmax, kT=10) sequence s = 0&shuffle(tagset(99))&0 atom Emin:=Es(s) -- E0 printf(1,"E(s0) %f\n",Emin) -- random starter for k=0 to kmax do if mod(k,kmax/10)=0 then printf(1,"k:%,10d T: %8.4f Es: %8.4f\n",{k,T(k,kmax,kT),Es(s)}) if k=kmax then exit end if -- avoid exp(x,-inf) end if integer u = rand(99), -- city index 1 99 cv = s[u+1]+dirs[rand(8)] -- city number if cv>0 and cv<100 -- not bogus city and dist(s[u+1],cv)<5 then -- and true neighbour integer v = s[cv+1] -- city index atom deltae := dE(s,u,v); if deltae<0 -- always move if negative or P(deltae,k,kmax,kT)>=rnd() then {s[u+1],s[v+1]} = {s[v+1],s[u+1]} Emin += deltae end if end if end for printf(1,"E(s_final) %f\n",Emin) printf(1,"Path:\n") pp(s,{pp_IntFmt,"%2d",pp_IntCh,false}) end procedure sa(1_000_000,1)

~~Output:~~
E(s0) 515.164811 k: 0 T: 1.0000 Es: 515.1648 k: 100,000 T: 0.9000 Es: 189.3123 k: 200,000 T: 0.8000 Es: 198.7498 k: 300,000 T: 0.7000 Es: 158.2189 k: 400,000 T: 0.6000 Es: 165.4813 k: 500,000 T: 0.5000 Es: 156.3467 k: 600,000 T: 0.4000 Es: 142.7928 k: 700,000 T: 0.3000 Es: 128.0352 k: 800,000 T: 0.2000 Es: 121.7794 k: 900,000 T: 0.1000 Es: 121.2328 k: 1,000,000 T: 0.0000 Es: 121.1291 E(s_final) 121.129115 Path: { 0,10,62,63,64,65,76,75,84,85,95,86,96,97,87,77,67,66,56,46,47,48,49,59,69, 79,89,99,98,88,78,68,58,57,37,38,27,26,36,35,45,55,54,53,52,43,33,23,22,32, 42,41,51,61,60,50,40,30,31,21,20,11,12, 2, 3, 4, 5, 6,17,18,28,39,29,19, 9, 8, 7,16,15,24,44,74,83,93,94,92,91,71,70,90,80,81,82,72,73,34,25,14,13, 1, 0}
~~Racket~~
~~<lang racket>~~
~~lang racket~~
(require racket/fixnum)
(define current-dim (make-parameter 10)) (define current-dim-- (make-parameter 9)) (define current-dim² (make-parameter 100)) (define current-kT (make-parameter 1)) (define current-k-max (make-parameter 1000000)) (define current-monitor-frequency (make-parameter 100000)) (define current-monitor (make-parameter
~~(λ (s k T E) (when (zero? (modulo k (current-monitor-frequency))) (printf "T:~a E:~a~%" (~r T) E)))))~~
~~Simulated Annealing Solver~~
~~(define (P ΔE T)~~
~~(if (negative? ΔE) 1 (exp (- (/ ΔE T)))))~~
~~(define (solve/SA s₀ next-s k-max temperature E monitor)~~
(for*/fold ((s s₀) (E_s (E s₀))) ((k (in-range k-max))) (define T (temperature k k-max)) (when monitor (monitor s k T E_s)) (let* ((s´ (next-s s k)) (E_s´ (E s´)) (ΔE (- E_s´ E_s))) (if (>= (P ΔE T) (random)) (values s´ E_s´) (values s E_s)))))
~~(define (temperature k k-max)~~
~~(* (current-kT) (- 1 (/ k k-max))))~~
~~TSP Problem~~
(struct tsp (path indices Es ΣE) #:transparent)
(define (y/x i d) (quotient/remainder i d))
(define (dist a b (d (current-dim)))
~~(let-values (([ay ax] (y/x a d)) ([by bx] (y/x b d))) (sqrt (+ (sqr (- ay by)) (sqr (- ax bx))))))~~
~~(define (indices->tsp indices)~~
(define path (make-fxvector (current-dim²))) (for ((i indices) (n (current-dim²))) (fxvector-set! path i n)) (define Es (for/vector #:length (fxvector-length path) ((a (in-fxvector path)) (b (in-sequences (in-fxvector path 1) (in-value (fxvector-ref path 0))))) (dist a b))) (tsp path indices Es (for/sum ((E Es)) E)))
(define (dir->delta dir (dim (current-dim))) (case dir [(l) -1] [(r) +1] [(u) (- dim)] [(d) dim]))
(define (invalid-direction? x y d (mx (current-dim--)))
~~(match* (x y d) ((0 _ 'l) #t) (((== mx) _ 'r) #t) ((_ 0 'u) #t) ((_ (== mx) 'd) #t) ((_ _ _) #f)))~~
~~extended to take k to reset numerical drift from the Δ calculation~~
~~(define (tsp:swap-one-neighbour t k)~~
(define dim (current-dim)) (define dim² (current-dim²)) (define candidate (random dim²)) (define-values [cy cx] (quotient/remainder candidate dim)) (define dir (vector-ref #(l r u d) (random 4))) (cond [(invalid-direction? cx cy dir) (tsp:swap-one-neighbour t k)] [else (define delta (dir->delta dir)) (define neighbour (+ candidate delta)) (define path´ (fxvector-copy (tsp-path t))) (define indices´ (fxvector-copy (tsp-indices t))) (define cand-idx (fxvector-ref (tsp-indices t) candidate)) (define ngbr-idx (fxvector-ref (tsp-indices t) neighbour)) (fxvector-set! path´ cand-idx neighbour) (fxvector-set! path´ ngbr-idx candidate) (fxvector-set! indices´ candidate ngbr-idx) (fxvector-set! indices´ neighbour cand-idx) (define Es (tsp-Es t)) (define Es´ (vector-copy Es))
(let* ((cand-idx++ (modulo (add1 cand-idx) dim²)) (cand-idx-- (modulo (sub1 cand-idx) dim²)) (ngbr-idx++ (modulo (add1 ngbr-idx) dim²)) (ngbr-idx-- (modulo (sub1 ngbr-idx) dim²))) (define Σold-E-around-nodes (+ (vector-ref Es cand-idx) (vector-ref Es cand-idx--) (vector-ref Es ngbr-idx) (vector-ref Es ngbr-idx--))) (define E´-at-cand (dist (fxvector-ref path´ cand-idx) (fxvector-ref path´ cand-idx++))) (define E´-pre-cand (dist (fxvector-ref path´ cand-idx) (fxvector-ref path´ cand-idx--))) (define E´-at-ngbr (dist (fxvector-ref path´ ngbr-idx) (fxvector-ref path´ ngbr-idx++))) (define E´-pre-ngbr (dist (fxvector-ref path´ ngbr-idx) (fxvector-ref path´ ngbr-idx--))) (vector-set! Es´ cand-idx E´-at-cand) (vector-set! Es´ cand-idx-- E´-pre-cand) (vector-set! Es´ ngbr-idx E´-at-ngbr) (vector-set! Es´ ngbr-idx-- E´-pre-ngbr)
(define ΔE (- (+ E´-at-cand E´-pre-cand E´-at-ngbr E´-pre-ngbr) Σold-E-around-nodes)) (tsp path´ indices´ Es´ (if (zero? (modulo k 1000)) (for/sum ((e Es´)) e) (+ (tsp-ΣE t) ΔE))))]))
~~(define (tsp:random-state)~~
~~(indices->tsp (for/fxvector ((i (shuffle (range (current-dim²))))) i)))~~
~~(define (Simulated-annealing)~~
(define-values (solution solution-E) (solve/SA (tsp:random-state) tsp:swap-one-neighbour (current-k-max) temperature tsp-ΣE (current-monitor))) (displayln (tsp-path solution)) (displayln solution-E))
~~(module+ main~~
~~(Simulated-annealing))</lang>~~

~~Output:~~
T:1 E:552.4249706051347 T:0.9 E:204.89460292101052 T:0.8 E:178.6926191428981 T:0.7 E:157.77681824512447 T:0.6 E:145.91227208091533 T:0.5 E:127.16624235784029 T:0.4 E:119.56239369288322 T:0.3 E:111.92798007771523 T:0.2 E:102.24264068711928 T:0.1 E:101.65685424949237 #fx(67 68 78 88 98 99 89 79 69 59 49 48 38 39 29 19 9 8 7 17 18 28 27 37 36 26 25 15 16 6 5 4 3 12 13 14 24 34 44 54 53 43 33 23 22 32 31 21 11 2 1 0 10 20 30 40 41 42 52 62 61 51 50 60 70 71 81 80 90 91 92 93 94 84 83 82 72 73 63 64 74 75 65 55 45 35 46 47 58 57 56 66 76 86 85 95 96 97 87 77) 101.65685424949237
~~Raku~~
~~(formerly Perl 6)~~
Translation of: Go
<lang perl6># simulation setup my \cities = 100; # number of cities my \kmax = 1e6; # iterations to run my \kT = 1; # initial 'temperature'
die 'City count must be a perfect square.' if cities.sqrt != cities.sqrt.Int;
~~locations of (up to) 8 neighbors, with grid size derived from number of cities~~
~~my \gs = cities.sqrt; my \neighbors = [1, -1, gs, -gs, gs-1, gs+1, -(gs+1), -(gs-1)];~~
~~matrix of distances between cities~~
~~my \D = [;]; for 0 ..^ cities² -> \j {~~
~~my (\ab, \cd) = (j/cities, j%cities)».Int; my (\a, \b, \c, \d) = (ab/gs, ab%gs, cd/gs, cd%gs)».Int; D[ab;cd] = sqrt (a - c)² + (b - d)²~~
}
sub T(\k, \kmax, \kT) { (1 - k/kmax) × kT } # temperature function, decreases to 0 sub P(\ΔE, \k, \kmax, \kT) { exp( -ΔE / T(k, kmax, kT)) } # probability to move from s to s_next sub Es(\path) { sum (D[ path[$_]; path[$_+1] ] for 0 ..^ +path-1) } # energy at s, to be minimized
~~variation of E, from state s to state s_next~~
~~sub delta-E(\s, \u, \v) {~~
my (\a, \b, \c, \d) = D[s[u-1];s[u]], D[s[u+1];s[u]], D[s[v-1];s[v]], D[s[v+1];s[v]]; my (\na, \nb, \nc, \nd) = D[s[u-1];s[v]], D[s[u+1];s[v]], D[s[v-1];s[u]], D[s[v+1];s[u]]; if v == u+1 { return (na + nd) - (a + d) } elsif u == v+1 { return (nc + nb) - (c + b) } else { return (na + nb + nc + nd) - (a + b + c + d) }
}
~~E(s0), initial state~~
my \s = @ = flat 0, (1 ..^ cities).pick(*), 0; my \E-min-global = my \E-min = $ = Es(s);
for 0 ..^ kmax -> \k {
~~printf "k:%8u T:%4.1f Es: %3.1f\n" , k, T(k, kmax, kT), Es(s) if k % (kmax/10) == 0;~~
~~# valid candidate cities (exist, adjacent) my \cv = neighbors[(^8).roll] + s[ my \u = 1 + (^(cities-1)).roll ]; next if cv ≤ 0 or cv ≥ cities or D[s[u];cv] > sqrt(2);~~
my \v = s[cv]; my \ΔE = delta-E(s, u, v); if ΔE < 0 or P(ΔE, k, kmax, kT) ≥ rand { # always move if negative (s[u], s[v]) = (s[v], s[u]); E-min += ΔE; E-min-global = E-min if E-min < E-min-global; }
}
say "\nE(s_final): " ~ E-min-global.fmt('%.1f'); say "Path:\n" ~ s».fmt('%2d').rotor(20,:partial).join: "\n";</lang>

~~Output:~~
k: 0 T: 1.0 Es: 522.0 k: 100000 T: 0.9 Es: 185.3 k: 200000 T: 0.8 Es: 166.1 k: 300000 T: 0.7 Es: 174.7 k: 400000 T: 0.6 Es: 146.9 k: 500000 T: 0.5 Es: 140.2 k: 600000 T: 0.4 Es: 127.5 k: 700000 T: 0.3 Es: 115.9 k: 800000 T: 0.2 Es: 111.9 k: 900000 T: 0.1 Es: 109.4 E(s_final): 109.4 Path: 0 10 20 30 40 50 60 84 85 86 96 97 87 88 98 99 89 79 78 77 67 68 69 59 58 57 56 66 76 95 94 93 92 91 90 80 70 81 82 83 73 72 71 62 63 64 74 75 65 55 54 53 52 61 51 41 31 21 22 32 42 43 44 45 46 35 34 24 23 33 25 15 16 26 36 47 37 27 17 18 28 38 48 49 39 29 19 9 8 7 6 5 4 14 13 12 11 2 3 1 0
~~Scheme~~
Works with: CHICKEN version 5.3.0
Library: r7rs
Library: srfi-1
Library: srfi-27
Library: srfi-144
Library: format

Note: Each line of the printed table gives E(s) before the next path is chosen. Thus the top line describes the initial state.

<lang scheme>(cond-expand
~~(r7rs) (chicken (import r7rs)))~~
~~(import (scheme base)) (import (scheme inexact)) (import (scheme write)) (import (only (srfi 1) delete)) (import (only (srfi 1) iota)) (import (srfi 27)) ; Random numbers.~~

~~You can do without SRFI-144 by changing fl+ to +, etc.~~

(import (srfi 144)) ; Optimizations for flonums.
(cond-expand
~~(chicken (import (format))) (else))~~
(random-source-randomize! default-random-source)
(define (n->ij n)
~~(truncate/ n 10))~~
~~(define (ij->n i j)~~
~~(+ (* 10 i) j))~~
~~(define neighbor-offsets~~
~~'((0 . 1) (1 . 0) (1 . 1) (0 . -1) (-1 . 0) (-1 . -1) (1 . -1) (-1 . 1)))~~
~~(define (neighborhood n)~~
(let-values (((i j) (n->ij n))) (let recurs ((offsets neighbor-offsets)) (if (null? offsets) '() (let* ((offs (car offsets)) (i^ (+ i (car offs))) (j^ (+ j (cdr offs)))) (if (and (not (negative? i^)) (not (negative? j^)) (< i^ 10) (< j^ 10)) (cons (ij->n i^ j^) (recurs (cdr offsets))) (recurs (cdr offsets))))))))
~~(define (distance m n)~~
~~(let-values (((im jm) (n->ij m)) ((in jn) (n->ij n))) (flsqrt (inexact (+ (square (- im in)) (square (- jm jn)))))))~~
~~(define (shuffle! vec i n)~~
;; A Fisher-Yates shuffle of n elements of vec, starting at index i. (do ((j 0 (+ j 1))) ((= j n)) (let* ((k (+ i j (random-integer (- n j)))) (xi (vector-ref vec i)) (xk (vector-ref vec k))) (vector-set! vec i xk) (vector-set! vec k xi))))
~~(define (make-s0)~~
~~(let ((vec (list->vector (iota 100)))) (shuffle! vec 1 99) vec))~~
~~(define (swap-s-elements! vec u v)~~
(let loop ((j 1) (iu 0) (iv 0)) (cond ((positive? iu) (if (= (vector-ref vec j) v) (begin (vector-set! vec iu v) (vector-set! vec j u)) (loop (+ j 1) iu iv))) ((positive? iv) (if (= (vector-ref vec j) u) (begin (vector-set! vec j v) (vector-set! vec iv u)) (loop (+ j 1) iu iv))) ((= (vector-ref vec j) u) (loop (+ j 1) j 0)) ((= (vector-ref vec j) v) (loop (+ j 1) 0 j)) (else (loop (+ j 1) 0 0)))))
~~(define (update-s! vec)~~
(let* ((u (+ 1 (random-integer 99))) (neighbors (delete 0 (neighborhood u) =)) (n (length neighbors)) (v (list-ref neighbors (random-integer n)))) (swap-s-elements! vec u v)))
~~(define (s->s vec) ; s_k -> s_(k + 1)~~
~~(let ((vec^ (vector-copy vec))) (update-s! vec^) vec^))~~
~~(define (path-length vec) ; E(s)~~
(let loop ((plen (distance (vector-ref vec 0) (vector-ref vec 99))) (x (vector-ref vec 0)) (i 1)) (if (= i 100) plen (let ((y (vector-ref vec i))) (loop (fl+ plen (distance x y)) y (+ i 1))))))
~~(define (make-temperature-procedure kT kmax)~~
~~(let ((kT (inexact kT)) (kmax (inexact kmax))) (lambda (k) (fl* kT (fl- 1.0 (fl/ (inexact k) kmax))))))~~
~~(define (probability delta-E T)~~
~~(if (flnegative? delta-E) 1.0 (if (flzero? T) 0.0 (flexp (fl- (fl/ delta-E T))))))~~
~~(define fmt10 (string-append " k T E(s)~%"~~
~~" -----------------------------~%"))~~
(define fmt20 " ~7D ~3,1F ~12,5F~%")
(define (simulate-annealing kT kmax)
(let* ((temperature (make-temperature-procedure kT kmax)) (s0 (make-s0)) (E0 (path-length s0)) (kmax/10 (truncate-quotient kmax 10)) (show (lambda (k T E) (if (zero? (truncate-remainder k kmax/10)) (cond-expand (chicken (format #t fmt20 k T E)) (else (display k) (display " ") (display T) (display " ") (display E) (newline))))))) (cond-expand (chicken (format #t fmt10)) (else)) (let loop ((k 0) (s s0) (E E0)) (if (= k (+ 1 kmax)) s (let* ((T (temperature k)) (_ (show k T E)) (s^ (s->s s)) (E^ (path-length s^)) (delta-E (fl- E^ E)) (P (probability delta-E T))) (if (or (fl=? P 1.0) (fl<=? (random-real) P)) (loop (+ k 1) s^ E^) (loop (+ k 1) s E)))))))
~~(define (display-path vec)~~
(do ((i 0 (+ i 1))) ((= i 100)) (let ((x (vector-ref vec i))) (when (< x 10) (display " ")) (display x) (display " -> ") (when (= 7 (truncate-remainder i 8)) (newline)))) (let ((x (vector-ref vec 0))) (when (< x 10) (display " ")) (display x)))
(define kT 1) (define kmax 1000000)
(newline) (display " kT: ") (display kT) (newline) (display " kmax: ") (display kmax) (newline) (newline) (define s-final (simulate-annealing kT kmax)) (newline) (display "Final path:") (newline) (display-path s-final) (newline) (newline) (cond-expand
~~(chicken (format #t "Final E(s): ~,5F~%" (path-length s-final))) (else (display "Final E(s): ") (display (path-length s-final)) (newline)))~~
~~(newline)</lang>~~

~~Output:~~
~~An example run:~~
$ csc -O5 -X r7rs -R r7rs sa.scm && ./sa kT: 1 kmax: 1000000 k T E(s) ----------------------------- 0 1.0 422.71361 100000 0.9 185.28073 200000 0.8 171.70817 300000 0.7 156.66104 400000 0.6 145.07621 500000 0.5 130.88759 600000 0.4 115.34219 700000 0.3 112.27113 800000 0.2 105.37820 900000 0.1 103.89949 1000000 0.0 103.89949 Final path: 0 -> 1 -> 2 -> 3 -> 4 -> 6 -> 7 -> 8 -> 9 -> 19 -> 29 -> 39 -> 38 -> 37 -> 47 -> 48 -> 49 -> 58 -> 59 -> 69 -> 79 -> 89 -> 99 -> 98 -> 97 -> 96 -> 95 -> 94 -> 84 -> 83 -> 93 -> 92 -> 82 -> 72 -> 62 -> 71 -> 81 -> 91 -> 90 -> 80 -> 70 -> 60 -> 61 -> 50 -> 40 -> 41 -> 51 -> 52 -> 63 -> 73 -> 74 -> 75 -> 85 -> 86 -> 76 -> 77 -> 87 -> 88 -> 78 -> 68 -> 67 -> 57 -> 56 -> 66 -> 65 -> 64 -> 55 -> 45 -> 46 -> 36 -> 35 -> 25 -> 26 -> 27 -> 28 -> 18 -> 17 -> 16 -> 15 -> 5 -> 14 -> 13 -> 23 -> 24 -> 34 -> 44 -> 54 -> 53 -> 43 -> 33 -> 42 -> 32 -> 31 -> 30 -> 20 -> 21 -> 22 -> 12 -> 11 -> 10 -> 0 Final E(s): 103.89949

~~A different E(s)~~
Library: srfi-143

Here E(s) is the sum of squares of differences, so that E(s) is an integer. Also I use kT=1.5 and twice as large a kmax.
(I also demonstrate some of SRFI 143. Note that CHICKEN has type annotations as an alternative to using SRFI 143 and SRFI 144, but the SRFI extensions are more portable.)

<lang scheme>(cond-expand
~~(r7rs) (chicken (import r7rs)))~~
~~(import (scheme base)) (import (scheme inexact)) (import (scheme write)) (import (only (srfi 1) delete)) (import (only (srfi 1) iota)) (import (srfi 27)) ; Random numbers.~~

The following import is CHICKEN-specific, but your Scheme likely

has Common Lisp formatting somewhere.

~~(import (format)) ; Common Lisp formatting.~~

You can do without SRFI-143 or SRFI-144 by changing fx+ or fl+ to

+, etc.

(import (srfi 143)) ; Optimizations for fixnums. (import (srfi 144)) ; Optimizations for flonums.
(random-source-randomize! default-random-source)
(define (n->ij n)
~~(values (fxquotient n 10) (fxremainder n 10)))~~
~~(define (ij->n i j)~~
~~(fx+ (fx* 10 i) j))~~
~~(define neighbor-offsets~~
~~'((0 . 1) (1 . 0) (1 . 1) (0 . -1) (-1 . 0) (-1 . -1) (1 . -1) (-1 . 1)))~~
~~(define (neighborhood n)~~
(let-values (((i j) (n->ij n))) (let recurs ((offsets neighbor-offsets)) (if (null? offsets) '() (let* ((offs (car offsets)) (i^ (fx+ i (car offs))) (j^ (fx+ j (cdr offs)))) (if (and (not (fxnegative? i^)) (not (fxnegative? j^)) (fx<? i^ 10) (fx<? j^ 10)) (cons (ij->n i^ j^) (recurs (cdr offsets))) (recurs (cdr offsets))))))))
~~(define (distance**2 m n)~~
~~(let-values (((im jm) (n->ij m)) ((in jn) (n->ij n))) (fx+ (fxsquare (fx- im in)) (fxsquare (fx- jm jn)))))~~
~~(define (shuffle! vec i n)~~
;; A Fisher-Yates shuffle of n elements of vec, starting at index i. (do ((j 0 (+ j 1))) ((= j n)) (let* ((k (+ i j (random-integer (- n j)))) (xi (vector-ref vec i)) (xk (vector-ref vec k))) (vector-set! vec i xk) (vector-set! vec k xi))))
~~(define (make-s0)~~
~~(let ((vec (list->vector (iota 100)))) (shuffle! vec 1 99) vec))~~
~~(define (swap-s-elements! vec u v)~~
(let loop ((j 1) (iu 0) (iv 0)) (cond ((fxpositive? iu) (if (fx=? (vector-ref vec j) v) (begin (vector-set! vec iu v) (vector-set! vec j u)) (loop (fx+ j 1) iu iv))) ((fxpositive? iv) (if (fx=? (vector-ref vec j) u) (begin (vector-set! vec j v) (vector-set! vec iv u)) (loop (fx+ j 1) iu iv))) ((fx=? (vector-ref vec j) u) (loop (fx+ j 1) j 0)) ((fx=? (vector-ref vec j) v) (loop (fx+ j 1) 0 j)) (else (loop (fx+ j 1) 0 0)))))
~~(define (update-s! vec)~~
(let* ((u (fx+ 1 (random-integer 99))) (neighbors (delete 0 (neighborhood u) fx=?)) (n (length neighbors)) (v (list-ref neighbors (random-integer n)))) (swap-s-elements! vec u v)))
~~(define (s->s vec) ; s_k -> s_(k + 1)~~
~~(let ((vec^ (vector-copy vec))) (update-s! vec^) vec^))~~
~~(define (path-length vec)~~
(let loop ((plen (flsqrt (inexact (distance**2 (vector-ref vec 0) (vector-ref vec 99))))) (x (vector-ref vec 0)) (i 1)) (if (fx=? i 100) plen (let ((y (vector-ref vec i))) (loop (fl+ plen (flsqrt (inexact (distance**2 x y)))) y (fx+ i 1))))))
~~(define (E_s vec)~~
(let loop ((E (distance**2 (vector-ref vec 0) (vector-ref vec 99))) (x (vector-ref vec 0)) (i 1)) (if (fx=? i 100) E (let ((y (vector-ref vec i))) (loop (fx+ E (distance**2 x y)) y (fx+ i 1))))))
~~(define (make-temperature-procedure kT kmax)~~
~~(let ((kT (inexact kT)) (kmax (inexact kmax))) (lambda (k) (fl* kT (fl- 1.0 (fl/ (inexact k) kmax))))))~~
~~(define (probability delta-E T)~~
~~(if (fxnegative? delta-E) 1.0 (if (flzero? T) 0.0 (flexp (fl- (fl/ (inexact delta-E) T))))))~~
~~(define fmt10 (string-append~~
~~" k T E(s) path length~%" " ---------------------------------------~%"))~~
(define fmt20 " ~7D ~7,2F ~8D ~14,5F~%")
(define (simulate-annealing kT kmax)
(let* ((temperature (make-temperature-procedure kT kmax)) (s0 (make-s0)) (E0 (E_s s0)) (kmax/10 (fxquotient kmax 10)) (show (lambda (k T E s) (when (fxzero? (fxremainder k kmax/10)) (format #t fmt20 k T E (path-length s)))))) (format #t fmt10) (let loop ((k 0) (s s0) (E E0)) (if (fx=? k (fx+ 1 kmax)) s (let* ((T (temperature k)) (_ (show k T E s)) (s^ (s->s s)) (E^ (E_s s^)) (delta-E (fx- E^ E)) (P (probability delta-E T))) (if (or (fl=? P 1.0) (fl<=? (random-real) P)) (loop (fx+ k 1) s^ E^) (loop (fx+ k 1) s E)))))))
~~(define (display-path vec)~~
(do ((i 0 (+ i 1))) ((= i 100)) (let ((x (vector-ref vec i))) (when (< x 10) (display " ")) (display x) (display " -> ") (when (= 7 (truncate-remainder i 8)) (newline)))) (let ((x (vector-ref vec 0))) (when (< x 10) (display " ")) (display x)))
(define kT 1.5) (define kmax 2000000)
(newline) (display " kT: ") (display kT) (newline) (display " kmax: ") (display kmax) (newline) (newline) (define s-final (simulate-annealing kT kmax)) (newline) (display "Final path:") (newline) (display-path s-final) (newline) (newline) (format #t "Final E(s): ~,5F~%" (E_s s-final)) (format #t "Final path length: ~,5F~%" (path-length s-final)) (newline)</lang>

~~Output:~~
$ csc -O5 -X r7rs -R r7rs sa2.scm && ./sa2 kT: 1.5 kmax: 2000000 k T E(s) path length --------------------------------------- 0 1.50 2298 384.59396 200000 1.35 180 127.94332 400000 1.20 168 124.60675 600000 1.05 148 118.07012 800000 0.90 130 111.94113 1000000 0.75 116 106.62742 1200000 0.60 114 105.55635 1400000 0.45 112 104.97056 1600000 0.30 104 101.65685 1800000 0.15 104 101.65685 2000000 0.00 104 101.65685 Final path: 0 -> 1 -> 2 -> 3 -> 4 -> 5 -> 15 -> 24 -> 34 -> 43 -> 44 -> 54 -> 53 -> 63 -> 64 -> 65 -> 66 -> 67 -> 77 -> 76 -> 75 -> 85 -> 84 -> 74 -> 73 -> 72 -> 62 -> 52 -> 42 -> 41 -> 31 -> 32 -> 33 -> 23 -> 13 -> 14 -> 25 -> 26 -> 27 -> 17 -> 16 -> 6 -> 7 -> 8 -> 9 -> 19 -> 18 -> 28 -> 29 -> 39 -> 49 -> 59 -> 58 -> 48 -> 38 -> 37 -> 47 -> 46 -> 36 -> 35 -> 45 -> 55 -> 56 -> 57 -> 68 -> 69 -> 79 -> 78 -> 88 -> 89 -> 99 -> 98 -> 97 -> 87 -> 86 -> 96 -> 95 -> 94 -> 93 -> 83 -> 82 -> 92 -> 91 -> 90 -> 80 -> 81 -> 71 -> 70 -> 60 -> 61 -> 51 -> 50 -> 40 -> 30 -> 20 -> 21 -> 22 -> 12 -> 11 -> 10 -> 0 Final E(s): 104.00000 Final path length: 101.65685

~~A second run shows E(s) temporarily increasing:~~
kT: 1.5 kmax: 2000000 k T E(s) path length --------------------------------------- 0 1.50 2132 368.24125 200000 1.35 142 115.58483 400000 1.20 146 116.75641 600000 1.05 148 117.40669 800000 0.90 124 109.27770 1000000 0.75 112 104.97056 1200000 0.60 124 109.45584 1400000 0.45 114 105.55635 1600000 0.30 108 103.31371 1800000 0.15 108 103.31371 2000000 0.00 108 103.31371 Final path: 0 -> 1 -> 11 -> 21 -> 31 -> 42 -> 52 -> 62 -> 72 -> 73 -> 63 -> 74 -> 64 -> 65 -> 55 -> 45 -> 54 -> 53 -> 43 -> 44 -> 34 -> 35 -> 36 -> 26 -> 25 -> 24 -> 23 -> 33 -> 32 -> 22 -> 12 -> 2 -> 3 -> 13 -> 14 -> 4 -> 5 -> 15 -> 16 -> 6 -> 7 -> 17 -> 27 -> 28 -> 18 -> 8 -> 9 -> 19 -> 29 -> 39 -> 38 -> 37 -> 47 -> 48 -> 49 -> 58 -> 59 -> 69 -> 68 -> 67 -> 57 -> 46 -> 56 -> 66 -> 76 -> 86 -> 87 -> 77 -> 78 -> 79 -> 89 -> 99 -> 88 -> 98 -> 97 -> 96 -> 95 -> 85 -> 75 -> 84 -> 94 -> 93 -> 83 -> 82 -> 92 -> 91 -> 90 -> 80 -> 81 -> 71 -> 70 -> 60 -> 61 -> 51 -> 50 -> 41 -> 40 -> 30 -> 20 -> 10 -> 0 Final E(s): 108.00000 Final path length: 103.31371

~~A third run shows E(s) temporarily increasing, and also achieves a path length less than 101:~~
kT: 1.5 kmax: 2000000 k T E(s) path length --------------------------------------- 0 1.50 2246 388.77550 200000 1.35 176 124.95651 400000 1.20 160 121.22854 600000 1.05 148 117.29304 800000 0.90 126 109.86348 1000000 0.75 118 107.45584 1200000 0.60 120 108.04163 1400000 0.45 108 103.31371 1600000 0.30 106 102.48528 1800000 0.15 102 100.82843 2000000 0.00 102 100.82843 Final path: 0 -> 1 -> 11 -> 12 -> 2 -> 3 -> 13 -> 14 -> 4 -> 5 -> 15 -> 16 -> 6 -> 7 -> 17 -> 27 -> 28 -> 18 -> 8 -> 9 -> 19 -> 29 -> 39 -> 38 -> 48 -> 49 -> 59 -> 69 -> 79 -> 78 -> 77 -> 67 -> 68 -> 58 -> 57 -> 47 -> 37 -> 36 -> 26 -> 25 -> 24 -> 23 -> 22 -> 21 -> 31 -> 32 -> 43 -> 44 -> 54 -> 53 -> 52 -> 62 -> 61 -> 51 -> 41 -> 42 -> 33 -> 34 -> 35 -> 45 -> 46 -> 56 -> 55 -> 65 -> 66 -> 76 -> 86 -> 87 -> 88 -> 89 -> 99 -> 98 -> 97 -> 96 -> 95 -> 94 -> 84 -> 85 -> 75 -> 74 -> 64 -> 63 -> 73 -> 83 -> 93 -> 92 -> 82 -> 72 -> 71 -> 81 -> 91 -> 90 -> 80 -> 70 -> 60 -> 50 -> 40 -> 30 -> 20 -> 10 -> 0 Final E(s): 102.00000 Final path length: 100.82843
~~Sidef~~
Translation of: Julia
~~<lang ruby>module TravelingSalesman {~~
~~# Eₛ: length(path) func Eₛ(distances, path) { var total = 0 [path, path.slice(1)].zip {|ci,cj| total += distances[ci-1][cj-1] } total }~~
~~# T: temperature func T(k, kmax, kT) { kT * (1 - k/kmax) }~~
~~# ΔE = Eₛ_new - Eₛ_old > 0 # Prob. to move if ΔE > 0, → 0 when T → 0 (fronzen state) func P(ΔE, k, kmax, kT) { exp(-ΔE / T(k, kmax, kT)) }~~
~~# ∆E from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..) # ∆E before swapping (u,v) # Quicker than Eₛ(s_next) - Eₛ(path) func dE(distances, path, u, v) {~~
~~var a = distances[path[u-1]-1][path[u]-1] var b = distances[path[u+1]-1][path[u]-1] var c = distances[path[v-1]-1][path[v]-1] var d = distances[path[v+1]-1][path[v]-1]~~
~~var na = distances[path[u-1]-1][path[v]-1] var nb = distances[path[u+1]-1][path[v]-1] var nc = distances[path[v-1]-1][path[u]-1] var nd = distances[path[v+1]-1][path[u]-1]~~
~~if (v == u+1) { return ((na+nd) - (a+d)) }~~
~~if (u == v+1) { return ((nc+nb) - (c+b)) }~~
~~return ((na+nb+nc+nd) - (a+b+c+d)) }~~
~~const dirs = [1, -1, 10, -10, 9, 11, -11, -9]~~
~~func _prettypath(path) { path.slices(10).map { .map{ "%3s" % _ }.join(', ') }.join("\n") }~~
~~func findpath(distances, kmax, kT) {~~
~~const n = distances.len const R = 2..n~~
~~var path = [1, R.shuffle..., 1] var Emin = Eₛ(distances, path)~~
~~printf("# Entropy(s₀) = s%10.2f\n", Emin) printf("# Random path:\n%s\n\n", _prettypath(path))~~
~~for k in (1 .. kmax) {~~
~~if (k % (kmax//10) == 0) { printf("k: %10d | T: %8.4f | Eₛ: %8.4f\n", k, T(k, kmax, kT), Eₛ(distances, path)) }~~
~~var u = R.rand var v = (path[u-1] + dirs.rand) v ~~ R || next~~
~~var δE = dE(distances, path, u-1, v-1) if ((δE < 0) || (P(δE, k, kmax, kT) >= 1.rand)) { path.swap(u-1, v-1) Emin += δE } }~~
~~printf("k: %10d | T: %8.4f | Eₛ: %8.4f\n", kmax, T(kmax, kmax, kT), Eₛ(distances, path)) say ("\n# Found path:\n", _prettypath(path)) return path }~~
}
var citydist = {|ci|
~~{ |cj| var v1 = Vec(ci%10, ci//10) var v2 = Vec(cj%10, cj//10) v1.dist(v2) }.map(1..100)~~
}.map(1..100)
TravelingSalesman::findpath(citydist, 1e6, 1)</lang>

~~Output:~~
# Entropy(s₀) = 520.29 # Random path: 1, 10, 79, 52, 24, 9, 58, 11, 42, 4 15, 87, 62, 88, 21, 91, 99, 84, 61, 14 5, 17, 33, 95, 74, 31, 40, 13, 37, 69 6, 22, 97, 45, 56, 63, 75, 83, 53, 41 3, 47, 89, 80, 78, 98, 46, 18, 25, 51 93, 16, 50, 30, 48, 8, 66, 68, 59, 73 49, 96, 36, 32, 100, 27, 76, 44, 64, 39 90, 82, 20, 12, 54, 86, 29, 81, 26, 72 60, 94, 35, 92, 43, 7, 85, 55, 28, 57 23, 34, 65, 71, 38, 2, 77, 70, 19, 67 1 k: 100000 | T: 0.9000 | Eₛ: 185.1809 k: 200000 | T: 0.8000 | Eₛ: 168.6262 k: 300000 | T: 0.7000 | Eₛ: 146.5948 k: 400000 | T: 0.6000 | Eₛ: 140.1441 k: 500000 | T: 0.5000 | Eₛ: 129.5132 k: 600000 | T: 0.4000 | Eₛ: 132.8942 k: 700000 | T: 0.3000 | Eₛ: 124.2865 k: 800000 | T: 0.2000 | Eₛ: 120.0859 k: 900000 | T: 0.1000 | Eₛ: 115.0771 k: 1000000 | T: 0.0000 | Eₛ: 114.9728 k: 1000000 | T: 0.0000 | Eₛ: 114.9728 # Found path: 1, 2, 13, 3, 4, 5, 6, 7, 8, 9 19, 29, 18, 28, 27, 17, 16, 26, 25, 15 14, 24, 23, 12, 11, 10, 20, 21, 30, 40 41, 31, 32, 44, 45, 46, 47, 48, 49, 39 38, 37, 36, 35, 34, 42, 51, 50, 60, 61 52, 53, 54, 55, 56, 57, 58, 59, 69, 68 77, 67, 66, 65, 64, 62, 72, 71, 70, 80 81, 82, 74, 75, 76, 87, 88, 78, 79, 89 99, 98, 97, 96, 86, 85, 83, 91, 90, 100 92, 93, 94, 95, 84, 73, 63, 43, 33, 22 1
~~Wren~~
Translation of: Go
Library: Wren-math
Library: Wren-fmt
<lang ecmascript>import "random" for Random import "/math" for Math import "/fmt" for Fmt
// distances var calcDists = Fn.new {
var dists = List.filled(10000, 0) for (i in 0..9999) { var ab = (i/100).floor var cd = i % 100 var a = (ab/10).floor var b = ab % 10 var c = (cd/10).floor var d = cd % 10 dists[i] = Math.hypot(a-c, b-d) } return dists
}
var dists = calcDists.call() var dirs = [1, -1, 10, -10, 9, 11, -11, -9] // all 8 neighbors var rand = Random.new()
// index into lookup table of Nums var dist = Fn.new { |ci, cj| dists[cj*100 + ci] }
// energy at s, to be minimized var Es = Fn.new { |path|
~~var d = 0 for (i in 0...path.count-1) d = d + dist.call(path[i], path[i+1]) return d~~
}
// temperature function, decreases to 0 var T = Fn.new { |k, kmax, kT| (1 - k / kmax) * kT }
// variation of E, from state s to state s_next var dE = Fn.new { |s, u, v|
var su = s[u] var sv = s[v] // old var a = dist.call(s[u-1], su) var b = dist.call(s[u+1], su) var c = dist.call(s[v-1], sv) var d = dist.call(s[v+1], sv) // new var na = dist.call(s[u-1], sv) var nb = dist.call(s[u+1], sv) var nc = dist.call(s[v-1], su) var nd = dist.call(s[v+1], su) if (v == u+1) return (na + nd) - (a + d) if (u == v+1) return (nc + nb) - (c + b) return (na + nb + nc + nd) - (a + b + c + d)
}
// probability to move from s to s_next var P = Fn.new { |deltaE, k, kmax, kT| (-deltaE / T.call(k, kmax, kT)).exp }
// Simulated annealing var sa = Fn.new { |kmax, kT|
var temp = List.filled(99, 0) for (i in 0..98) temp[i] = i + 1 rand.shuffle(temp) var s = List.filled(101, 0) for (i in 0..98) s[i+1] = temp[i] // random path from 0 to 0 System.print("kT = %(kT)") System.print("E(s0) %(Es.call(s))\n") // random starter var Emin = Es.call(s) // E0 for (k in 0..kmax) { if (k % (kmax/10).floor == 0) { Fmt.print("k:$10d T: $8.4f Es: $8.4f", k, T.call(k, kmax, kT), Es.call(s)) } var u = rand.int(1, 100) // city index 1 to 99 var cv = s[u] + dirs[rand.int(8)] // city number if (cv <= 0 || cv >= 100) { // bogus city continue } if (dist.call(s[u], cv) > 5) { // check true neighbor (eg 0 9) continue } var v = s[cv] // city index var deltae = dE.call(s, u, v) if (deltae < 0 || // always move if negative P.call(deltae, k, kmax, kT) >= rand.float()) { s.swap(u, v) Emin = Emin + deltae } } System.print("\nE(s_final) %(Emin)") System.print("Path:") // output final state for (i in 0...s.count) { if (i > 0 && i % 10 == 0) System.print() Fmt.write("$4d", s[i]) } System.print()
}
sa.call(1e6, 1)</lang>

~~Output:~~
~~Sample run:~~
kT = 1 E(s0) 541.82779520458 k: 0 T: 1.0000 Es: 541.8278 k: 100000 T: 0.9000 Es: 187.1429 k: 200000 T: 0.8000 Es: 191.0983 k: 300000 T: 0.7000 Es: 171.7284 k: 400000 T: 0.6000 Es: 154.0549 k: 500000 T: 0.5000 Es: 147.0249 k: 600000 T: 0.4000 Es: 123.5822 k: 700000 T: 0.3000 Es: 121.5808 k: 800000 T: 0.2000 Es: 114.0930 k: 900000 T: 0.1000 Es: 112.6788 k: 1000000 T: 0.0000 Es: 112.6788 E(s_final) 112.67876668098 Path: 0 11 10 1 2 12 13 24 25 15 14 3 4 5 6 8 9 19 18 28 29 39 38 27 17 7 16 26 36 37 47 46 45 35 34 44 43 33 23 22 32 53 52 51 41 31 21 20 30 40 50 60 61 70 80 90 91 92 93 73 63 62 72 71 81 82 83 84 94 95 85 86 96 97 98 99 89 79 69 59 49 48 58 57 68 78 88 87 77 67 56 55 54 65 66 76 75 74 64 42 0
~~zkl~~
Translation of: EchoLisp
~~<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)~~
}); fcn dist(ci,cj){ _dists[cj*100 + ci] } // index into lookup table of floats
fcn Es(path) // E(s) = length(path): E(a,b,c)--> dist(a,b) + dist(b,c)
~~{ d:=Ref(0.0); path.reduce('wrap(a,b){ d.apply('+,dist(a,b)); b }); d.value }~~
// temperature() function fcn T(k,kmax,kT){ (1.0 - k.toFloat()/kmax)*kT }
// deltaE = Es_new - Es_old > 0 // probability to move if deltaE > 0, -->0 when T --> 0 (frozen state) fcn P(deltaE,k,kmax,kT){ (-deltaE/T(k,kmax,kT)).exp() } //-->Float
// deltaE from path ( .. a u b .. c v d ..) to (.. a v b ... c u d ..) // deltaE before swapping (u,v) fcn dE(s,u,v){ su,sv:=s[u],s[v]; //-->Float
~~// old a,b,c,d:=dist(s[u-1],su), dist(s[u+1],su), dist(s[v-1],sv), dist(s[v+1],sv); // new na,nb,nc,nd:=dist(s[u-1],sv), dist(s[u+1],sv), dist(s[v-1],su), dist(s[v+1],su);~~
~~if (v==u+1) (na+nd) - (a+d); else if(u==v+1) (nc+nb) - (c+b); else (na+nb+nc+nd) - (a+b+c+d);~~
}
// all 8 neighbours var [const] dirs=ROList(1, -1, 10, -10, 9, 11, -11, -9),
~~fmt="k:%10,d T: %8.4f Es: %8.4f".fmt; // since we use it twice~~
~~fcn sa(kmax,kT=10){~~
s:=List(0, [1..99].walk().shuffle().xplode(), 0); // random path from 0 to 0 println("E(s0) %f".fmt(Es(s))); // random starter Emin:=Es(s); // E0 foreach k in (kmax){ if(0==k%(kmax/10)) println(fmt(k,T(k,kmax,kT),Es(s))); u:=(1).random(100); // city index 1 99 cv:=s[u] + dirs[(0).random(8)]; // city number if(not (0<cv<100)) continue; // bogus city if(dist(s[u],cv)>5) continue; // check true neighbour (eg 0 9) v:=s.index(cv,1); // city index deltae:=dE(s,u,v); if(deltae<0 or // always move if negative
~~P(deltae,k,kmax,kT)>=(0.0).random(1)){ s.swap(u,v); Emin+=deltae;~~
~~} // (assert (= (round Emin) (round (Es s)))) }//foreach println(fmt(kmax,T(kmax-1,kmax,kT),Es(s))); println("E(s_final) %f".fmt(Emin)); println("Path: ",s.toString(*));~~
~~}</lang> <lang zkl>sa(0d1_000_000,1);</lang>~~

~~Output:~~
E(s0) 540.897080 k: 0 T: 1.0000 Es: 540.8971 k: 100,000 T: 0.9000 Es: 181.5102 k: 200,000 T: 0.8000 Es: 167.1944 k: 300,000 T: 0.7000 Es: 159.0975 k: 400,000 T: 0.6000 Es: 170.2344 k: 500,000 T: 0.5000 Es: 130.9919 k: 600,000 T: 0.4000 Es: 115.3422 k: 700,000 T: 0.3000 Es: 113.9280 k: 800,000 T: 0.2000 Es: 106.7924 k: 900,000 T: 0.1000 Es: 103.7213 k: 1,000,000 T: 0.0000 Es: 103.7213 E(s_final) 103.721349 Path: L(0,10,11,21,20,30,40,50,60,70,80,81,71,72,73,63,52,62,61,51,41,31,32,22,12,13,14,15,25,16,17,18,28,27,26,36,35,45,34,24,23,33,42,43,44,54,53,64,74,84,83,82,90,91,92,93,94,95,85,86,96,97,87,88,98,99,89,79,69,68,78,77,67,66,76,75,65,55,56,46,37,38,48,47,57,58,59,49,39,29,19,9,8,7,6,5,4,3,2,1,0)