htswanalysis/MACS/lib/gsl/gsl-1.11/siman/siman.c

   1 /* siman/siman.c
   2  *
   3  * Copyright (C) 2007 Brian Gough
   4  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Mark Galassi
   5  *
   6  * This program is free software; you can redistribute it and/or modify
   7  * it under the terms of the GNU General Public License as published by
   8  * the Free Software Foundation; either version 3 of the License, or (at
   9  * your option) any later version.
  10  *
  11  * This program is distributed in the hope that it will be useful, but
  12  * WITHOUT ANY WARRANTY; without even the implied warranty of
  13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14  * General Public License for more details.
  15  *
  16  * You should have received a copy of the GNU General Public License
  17  * along with this program; if not, write to the Free Software
  18  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  19  */
  20
  21 #include <config.h>
  22 #include <stdio.h>
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include <string.h>
  26 #include <assert.h>
  27
  28 #include <gsl/gsl_machine.h>
  29 #include <gsl/gsl_rng.h>
  30 #include <gsl/gsl_siman.h>
  31
  32 static inline double
  33 boltzmann(double E, double new_E, double T, gsl_siman_params_t *params)
  34 {
  35   double x = -(new_E - E) / (params->k * T);
  36   /* avoid underflow errors for large uphill steps */
  37   return (x < GSL_LOG_DBL_MIN) ? 0.0 : exp(x);
  38 }
  39
  40 static inline void
  41 copy_state(void *src, void *dst, size_t size, gsl_siman_copy_t copyfunc)
  42 {
  43   if (copyfunc) {
  44     copyfunc(src, dst);
  45   } else {
  46     memcpy(dst, src, size);
  47   }
  48 }
  49
  50 /* implementation of a basic simulated annealing algorithm */
  51
  52 void
  53 gsl_siman_solve (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
  54                  gsl_siman_step_t take_step,
  55                  gsl_siman_metric_t distance,
  56                  gsl_siman_print_t print_position,
  57                  gsl_siman_copy_t copyfunc,
  58                  gsl_siman_copy_construct_t copy_constructor,
  59                  gsl_siman_destroy_t destructor,
  60                  size_t element_size,
  61                  gsl_siman_params_t params)
  62 {
  63   void *x, *new_x, *best_x;
  64   double E, new_E, best_E;
  65   int i;
  66   double T, T_factor;
  67   int n_evals = 1, n_iter = 0, n_accepts, n_rejects, n_eless;
  68
  69   /* this function requires that either the dynamic functions (copy,
  70      copy_constructor and destrcutor) are passed, or that an element
  71      size is given */
  72   assert((copyfunc != NULL && copy_constructor != NULL && destructor != NULL)
  73          || (element_size != 0));
  74
  75   distance = 0 ; /* This parameter is not currently used */
  76   E = Ef(x0_p);
  77
  78   if (copyfunc) {
  79     x = copy_constructor(x0_p);
  80     new_x = copy_constructor(x0_p);
  81     best_x = copy_constructor(x0_p);
  82   } else {
  83     x = (void *) malloc (element_size);
  84     memcpy (x, x0_p, element_size);
  85     new_x = (void *) malloc (element_size);
  86     best_x =  (void *) malloc (element_size);
  87     memcpy (best_x, x0_p, element_size);
  88   }
  89
  90   best_E = E;
  91
  92   T = params.t_initial;
  93   T_factor = 1.0 / params.mu_t;
  94
  95   if (print_position) {
  96     printf ("#-iter  #-evals   temperature     position   energy\n");
  97   }
  98
  99   while (1) {
 100
 101     n_accepts = 0;
 102     n_rejects = 0;
 103     n_eless = 0;
 104
 105     for (i = 0; i < params.iters_fixed_T; ++i) {
 106
 107       copy_state(x, new_x, element_size, copyfunc);
 108
 109       take_step (r, new_x, params.step_size);
 110       new_E = Ef (new_x);
 111
 112       if(new_E <= best_E){
 113         if (copyfunc) {
 114           copyfunc(new_x,best_x);
 115         } else {
 116           memcpy (best_x, new_x, element_size);
 117         }
 118         best_E=new_E;
 119       }
 120
 121       ++n_evals;                /* keep track of Ef() evaluations */
 122       /* now take the crucial step: see if the new point is accepted
 123          or not, as determined by the boltzmann probability */
 124       if (new_E < E) {
 125
 126         if (new_E < best_E) {
 127           copy_state(new_x, best_x, element_size, copyfunc);
 128           best_E = new_E;
 129         }
 130
 131         /* yay! take a step */
 132         copy_state(new_x, x, element_size, copyfunc);
 133         E = new_E;
 134         ++n_eless;
 135
 136       } else if (gsl_rng_uniform(r) < boltzmann(E, new_E, T, &params)) {
 137         /* yay! take a step */
 138         copy_state(new_x, x, element_size, copyfunc);
 139         E = new_E;
 140         ++n_accepts;
 141
 142       } else {
 143         ++n_rejects;
 144       }
 145     }
 146
 147     if (print_position) {
 148       /* see if we need to print stuff as we go */
 149       /*       printf("%5d %12g %5d %3d %3d %3d", n_iter, T, n_evals, */
 150       /*           100*n_eless/n_steps, 100*n_accepts/n_steps, */
 151       /*           100*n_rejects/n_steps); */
 152       printf ("%5d   %7d  %12g", n_iter, n_evals, T);
 153       print_position (x);
 154       printf ("  %12g  %12g\n", E, best_E);
 155     }
 156
 157     /* apply the cooling schedule to the temperature */
 158     /* FIXME: I should also introduce a cooling schedule for the iters */
 159     T *= T_factor;
 160     ++n_iter;
 161     if (T < params.t_min) {
 162       break;
 163     }
 164   }
 165
 166   /* at the end, copy the result onto the initial point, so we pass it
 167      back to the caller */
 168   copy_state(best_x, x0_p, element_size, copyfunc);
 169
 170   if (copyfunc) {
 171     destructor(x);
 172     destructor(new_x);
 173     destructor(best_x);
 174   } else {
 175     free (x);
 176     free (new_x);
 177     free (best_x);
 178   }
 179 }
 180
 181 /* implementation of a simulated annealing algorithm with many tries */
 182
 183 void
 184 gsl_siman_solve_many (const gsl_rng * r, void *x0_p, gsl_siman_Efunc_t Ef,
 185                       gsl_siman_step_t take_step,
 186                       gsl_siman_metric_t distance,
 187                       gsl_siman_print_t print_position,
 188                       size_t element_size,
 189                       gsl_siman_params_t params)
 190 {
 191   /* the new set of trial points, and their energies and probabilities */
 192   void *x, *new_x;
 193   double *energies, *probs, *sum_probs;
 194   double Ex;                    /* energy of the chosen point */
 195   double T, T_factor;           /* the temperature and a step multiplier */
 196   int i;
 197   double u;                     /* throw the die to choose a new "x" */
 198   int n_iter;
 199
 200   if (print_position) {
 201     printf ("#-iter    temperature       position");
 202     printf ("         delta_pos        energy\n");
 203   }
 204
 205   x = (void *) malloc (params.n_tries * element_size);
 206   new_x = (void *) malloc (params.n_tries * element_size);
 207   energies = (double *) malloc (params.n_tries * sizeof (double));
 208   probs = (double *) malloc (params.n_tries * sizeof (double));
 209   sum_probs = (double *) malloc (params.n_tries * sizeof (double));
 210
 211   T = params.t_initial;
 212   T_factor = 1.0 / params.mu_t;
 213
 214   memcpy (x, x0_p, element_size);
 215
 216   n_iter = 0;
 217   while (1)
 218     {
 219       Ex = Ef (x);
 220       for (i = 0; i < params.n_tries - 1; ++i)
 221         {                       /* only go to N_TRIES-2 */
 222           /* center the new_x[] around x, then pass it to take_step() */
 223           sum_probs[i] = 0;
 224           memcpy ((char *)new_x + i * element_size, x, element_size);
 225           take_step (r, (char *)new_x + i * element_size, params.step_size);
 226           energies[i] = Ef ((char *)new_x + i * element_size);
 227           probs[i] = boltzmann(Ex, energies[i], T, &params);
 228         }
 229       /* now add in the old value of "x", so it is a contendor */
 230       memcpy ((char *)new_x + (params.n_tries - 1) * element_size, x, element_size);
 231       energies[params.n_tries - 1] = Ex;
 232       probs[params.n_tries - 1] = boltzmann(Ex, energies[i], T, &params);
 233
 234       /* now throw biased die to see which new_x[i] we choose */
 235       sum_probs[0] = probs[0];
 236       for (i = 1; i < params.n_tries; ++i)
 237         {
 238           sum_probs[i] = sum_probs[i - 1] + probs[i];
 239         }
 240       u = gsl_rng_uniform (r) * sum_probs[params.n_tries - 1];
 241       for (i = 0; i < params.n_tries; ++i)
 242         {
 243           if (u < sum_probs[i])
 244             {
 245               memcpy (x, (char *) new_x + i * element_size, element_size);
 246               break;
 247             }
 248         }
 249       if (print_position)
 250         {
 251           printf ("%5d\t%12g\t", n_iter, T);
 252           print_position (x);
 253           printf ("\t%12g\t%12g\n", distance (x, x0_p), Ex);
 254         }
 255       T *= T_factor;
 256       ++n_iter;
 257       if (T < params.t_min)
 258         {
 259           break;
 260         }
 261     }
 262
 263   /* now return the value via x0_p */
 264   memcpy (x0_p, x, element_size);
 265
 266   /*  printf("the result is: %g (E=%g)\n", x, Ex); */
 267
 268   free (x);
 269   free (new_x);
 270   free (energies);
 271   free (probs);
 272   free (sum_probs);
 273 }