Simulated annealing cleanup

Simulated annealing cleanup

[Heikki Orsila]

[-- Attachment #1: Type: text/plain, Size: 1578 bytes --]

Hello. I propose the attached patch to cleanup and optimize the 
simulated annealing algorithm.

Diffstat:

siman.c     |  112 ++++++++++++++++++++++++++++++------------------------------
siman_tsp.c |   11 ++++-
2 files changed, 66 insertions(+), 57 deletions(-)

Here's the description of changes:

siman/siman.c:
- remove useless "done" variable
- optimization: new_E <= best_E is only necessary when new_E < E
- new_E < best_E is the proper comparison, not new_E <= best_E. This
  reduces useless copying of state.
- fix whitespace use
- replace repeated occurence of:
    if (copyfunc) {
      copyfunc(src, dst);
    } else {
      memcpy(dst, src, size);
    }
  with copy_state().
- spelling: boltzman -> boltzmann
- make boltzmann function an inline function -> easier to modularize 
  the transition probability function later
- Use T *= T_factor instead of T /= params.mu_t. Avoids division -> faster.
- malloc() returns void * always, remove useless cast
- Print the best energy so far on progress. In the TSP example it reveals
  that optimization continues long after the optimum point has been found.
  On one test the optimum was reached with just 28% of total iterations.

siman/siman_tsp.c:
- initialize parameters using C99 initialization style. This way params
  structure may change inside and siman_tsp.c need not be updated. It
  is also more clear.

You should attribute the patch to "Heikki Orsila <heikki.orsila@iki.fi>"
if you apply it.

-- 
Heikki Orsila			Barbie's law:
heikki.orsila@iki.fi		"Math is hard, let's go shopping!"
http://www.iki.fi/shd

[-- Attachment #2: simulated_annealing_cleanup.diff --]
[-- Type: text/plain, Size: 7778 bytes --]

diff -urp gsl-1.9-org/siman/siman.c gsl-1.9/siman/siman.c
--- gsl-1.9-org/siman/siman.c	2005-11-15 18:22:47.000000000 +0200
+++ gsl-1.9/siman/siman.c	2007-05-25 22:25:02.000000000 +0300
@@ -33,6 +33,22 @@ double safe_exp (double x) /* avoid unde
 { 
   return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x);
 }
+
+static inline double
+boltzmann(double E, double new_E, double T, gsl_siman_params_t *params)
+{
+  return safe_exp (-(new_E - E) / (params->k * T));
+}
+
+static inline void
+copy_state(void *src, void *dst, size_t size, gsl_siman_copy_t copyfunc)
+{
+  if (copyfunc) {
+    copyfunc(src, dst);
+  } else {
+    memcpy(dst, src, size);
+  }
+}
       
 /* implementation of a basic simulated annealing algorithm */
 
@@ -49,8 +65,8 @@ gsl_siman_solve (const gsl_rng * r, void
 {
   void *x, *new_x, *best_x;
   double E, new_E, best_E;
-  int i, done;
-  double T;
+  int i;
+  double T, T_factor;
   int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless;
 
   /* this function requires that either the dynamic functions (copy,
@@ -67,67 +83,57 @@ gsl_siman_solve (const gsl_rng * r, void
     new_x = copy_constructor(x0_p);
     best_x = copy_constructor(x0_p);
   } else {
-    x = (void *) malloc (element_size);
+    x = malloc (element_size);
     memcpy (x, x0_p, element_size);
-    new_x = (void *) malloc (element_size);
-    best_x =  (void *) malloc (element_size);
+    new_x = malloc (element_size);
+    best_x =  malloc (element_size);
     memcpy (best_x, x0_p, element_size);
   }
 
   best_E = E;
 
   T = params.t_initial;
-  done = 0;
+  T_factor = 1.0 / params.mu_t;
 
   if (print_position) {
-    printf ("#-iter  #-evals   temperature     position   energy\n");
+    printf ("#-iter  #-evals   temperature     position   energy   bestenergy\n");
   }
 
-  while (!done) {
+  while (1) {
 
     n_accepts = 0;
     n_rejects = 0;
     n_eless = 0;
+
     for (i = 0; i < params.iters_fixed_T; ++i) {
-      if (copyfunc) {
-        copyfunc(x, new_x);
-      } else {
-        memcpy (new_x, x, element_size);
-      }
+
+      copy_state(x, new_x, element_size, copyfunc);
 
       take_step (r, new_x, params.step_size);
       new_E = Ef (new_x);
 
-      if(new_E <= best_E){
-        if (copyfunc) {
-          copyfunc(new_x,best_x);
-        } else {
-          memcpy (best_x, new_x, element_size);
-        }
-        best_E=new_E;
-      }
-
       ++n_evals;                /* keep track of Ef() evaluations */
+
       /* now take the crucial step: see if the new point is accepted
-         or not, as determined by the boltzman probability */
+         or not, as determined by the boltzmann probability */
       if (new_E < E) {
+
+	if (new_E < best_E) {
+	  copy_state(new_x, best_x, element_size, copyfunc);
+	  best_E = new_E;
+	}
+
         /* yay! take a step */
-        if (copyfunc) {
-          copyfunc(new_x, x);
-        } else {
-          memcpy (x, new_x, element_size);
-        }
+	copy_state(new_x, x, element_size, copyfunc);
         E = new_E;
         ++n_eless;
-      } else if (gsl_rng_uniform(r) < safe_exp (-(new_E - E)/(params.k * T)) ) {
+
+      } else if (gsl_rng_uniform(r) < boltzmann(E, new_E, T, &params)) {
         /* yay! take a step */
-        if (copyfunc) {
-          copyfunc(new_x, x);
-        } else {
-          memcpy(x, new_x, element_size);
-        }
+	copy_state(new_x, x, element_size, copyfunc);
         E = new_E;
         ++n_accepts;
+
       } else {
         ++n_rejects;
       }
@@ -140,25 +146,21 @@ gsl_siman_solve (const gsl_rng * r, void
       /*           100*n_rejects/n_steps); */
       printf ("%5d   %7d  %12g", n_iter, n_evals, T);
       print_position (x);
-      printf ("  %12g\n", E);
+      printf ("  %12g  %12g\n", E, best_E);
     }
 
     /* apply the cooling schedule to the temperature */
     /* FIXME: I should also introduce a cooling schedule for the iters */
-    T /= params.mu_t;
+    T *= T_factor;
     ++n_iter;
     if (T < params.t_min) {
-      done = 1;
+      break;
     }
   }
 
   /* at the end, copy the result onto the initial point, so we pass it
      back to the caller */
-  if (copyfunc) {
-    copyfunc(best_x, x0_p);
-  } else {
-    memcpy (x0_p, best_x, element_size);
-  }
+  copy_state(best_x, x0_p, element_size, copyfunc);
 
   if (copyfunc) {
     destructor(x);
@@ -185,8 +187,8 @@ gsl_siman_solve_many (const gsl_rng * r,
   void *x, *new_x;
   double *energies, *probs, *sum_probs;
   double Ex;                    /* energy of the chosen point */
-  double T;                     /* the temperature */
-  int i, done;
+  double T, T_factor;           /* the temperature and a step multiplier */
+  int i;
   double u;                     /* throw the die to choose a new "x" */
   int n_iter;
 
@@ -195,19 +197,19 @@ gsl_siman_solve_many (const gsl_rng * r,
     printf ("         delta_pos        energy\n");
   }
 
-  x = (void *) malloc (params.n_tries * element_size);
-  new_x = (void *) malloc (params.n_tries * element_size);
+  x = malloc (params.n_tries * element_size);
+  new_x = malloc (params.n_tries * element_size);
   energies = (double *) malloc (params.n_tries * sizeof (double));
   probs = (double *) malloc (params.n_tries * sizeof (double));
   sum_probs = (double *) malloc (params.n_tries * sizeof (double));
 
   T = params.t_initial;
-/*    memcpy (x, x0_p, element_size); */
+  T_factor = 1.0 / params.mu_t;
+
   memcpy (x, x0_p, element_size);
-  done = 0;
 
   n_iter = 0;
-  while (!done)
+  while (1)
     {
       Ex = Ef (x);
       for (i = 0; i < params.n_tries - 1; ++i)
@@ -217,12 +219,12 @@ gsl_siman_solve_many (const gsl_rng * r,
           memcpy ((char *)new_x + i * element_size, x, element_size);
           take_step (r, (char *)new_x + i * element_size, params.step_size);
           energies[i] = Ef ((char *)new_x + i * element_size);
-          probs[i] = safe_exp (-(energies[i] - Ex) / (params.k * T));
+          probs[i] = boltzmann(Ex, energies[i], T, &params);
         }
       /* now add in the old value of "x", so it is a contendor */
       memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size);
       energies[params.n_tries - 1] = Ex;
-      probs[params.n_tries - 1] = safe_exp (-(energies[i] - Ex) / (params.k * T));
+      probs[params.n_tries - 1] = boltzmann(Ex, energies[i], T, &params);
 
       /* now throw biased die to see which new_x[i] we choose */
       sum_probs[0] = probs[0];
@@ -235,7 +237,7 @@ gsl_siman_solve_many (const gsl_rng * r,
         {
           if (u < sum_probs[i])
             {
-              memcpy (x, (char *)new_x + i * element_size, element_size);
+              memcpy (x, (char *) new_x + i * element_size, element_size);
               break;
             }
         }
@@ -245,11 +247,11 @@ gsl_siman_solve_many (const gsl_rng * r,
           print_position (x);
           printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex);
         }
-      T /= params.mu_t;
+      T *= T_factor;
       ++n_iter;
       if (T < params.t_min)
-        {
-          done = 1;
+	{
+	  break;
         }
     }
 
diff -urp gsl-1.9-org/siman/siman_tsp.c gsl-1.9/siman/siman_tsp.c
--- gsl-1.9-org/siman/siman_tsp.c	2005-06-26 16:25:35.000000000 +0300
+++ gsl-1.9/siman/siman_tsp.c	2007-05-25 21:56:30.000000000 +0300
@@ -36,8 +36,15 @@
 #define MU_T 1.002              /* damping factor for temperature */
 #define T_MIN 5.0e-1
 
-gsl_siman_params_t params = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
-                             K, T_INITIAL, MU_T, T_MIN};
+gsl_siman_params_t params = {
+  .n_tries = N_TRIES,
+  .iters_fixed_T = ITERS_FIXED_T,
+  .step_size = STEP_SIZE,
+  .k = K,
+  .t_initial = T_INITIAL,
+  .mu_t = MU_T,
+  .t_min = T_MIN
+};
 
 struct s_tsp_city {
   const char * name;

Re: Simulated annealing cleanup

[Brian Gough]

At Fri, 25 May 2007 22:43:34 +0300,
Heikki Orsila wrote:
> Hello. I propose the attached patch to cleanup and optimize the 
> simulated annealing algorithm.

Thanks for the patch.  I will apply those changes.  Regarding C99, we
have to stick to C89 for compatibility for this project.

-- 
best regards,

Brian Gough

Network Theory Ltd,
Publishing Free Software Manuals --- http://www.network-theory.co.uk/