Simulated annealing: Difference between revisions

Content added Content deleted

Inline

Latest revision as of 05:59, 1 May 2024

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 them. 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.

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

----------------------------------------------------------------------

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.

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

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)

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

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

(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))

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)

Fortran

Translation of: Ada

Works with: gfortran version 11.3.0

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

Output:

$ 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

FreeBASIC

Uses 'LCS' function from Longest common subsequence#FreeBASIC:

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

Output:

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

Go

Translation of: zkl

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)
}

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

Icon

Translation of: Fortran

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

Output:

An example run:

$ 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

J

Implementation:

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=. (<s i.u,v) C. s
    dk=. pathcost sk
    dE=. dk-d
    if. (^-dE%T) >?0 do.
      s=.sk
      d=.dk
    end.
    if. 0=step|k do.
      echo k,T,d
    end.
  end.
  0,s,0
)

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:

   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

jq

Adapted from 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:

< /dev/urandom tr -cd '0-9' | fold -w 1 | jq -Rcnr -f sa.jq

Since gojq does not include jq's `_nwise/1`, here is a suitable def:

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

## 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)

Output:

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

Julia

Translation of: EchoLisp

Module:

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

Main:

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)

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

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()

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

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

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
(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))

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

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

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

RATFOR

Translation of: Fortran

Works with: ratfor77 version public domain 1.0

Works with: gfortran version 11.3.0

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

Output:

An example run:

$ 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

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.

(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)

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

(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)

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

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)

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

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)

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

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(*));
}

sa(0d1_000_000,1);

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)