vppinfra/vppinfra/anneal.c

   1 /*
   2   Copyright (c) 2011 Cisco and/or its affiliates.
   3
   4   * Licensed under the Apache License, Version 2.0 (the "License");
   5  * you may not use this file except in compliance with the License.
   6  * You may obtain a copy of the License at:
   7  *
   8  *     http://www.apache.org/licenses/LICENSE-2.0
   9  *
  10  * Unless required by applicable law or agreed to in writing, software
  11  * distributed under the License is distributed on an "AS IS" BASIS,
  12  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  13  * See the License for the specific language governing permissions and
  14  * limitations under the License.
  15 */
  16
  17 #include <vppinfra/anneal.h>
  18
  19 /*
  20  * Optimize an objective function by simulated annealing
  21  *
  22  * Here are a couple of short, easily-understood
  23  * descriptions of simulated annealing:
  24  *
  25  * http://www.cs.sandia.gov/opt/survey/sa.html
  26  * Numerical Recipes in C, 2nd ed., 444ff
  27  *
  28  * The description in the Wikipedia is not helpful.
  29  *
  30  * The algorithm tries to produce a decent answer to combinatorially
  31  * explosive optimization problems by analogy to slow cooling
  32  * of hot metal, aka annealing.
  33  *
  34  * There are (at least) three problem-dependent annealing parameters
  35  * to consider:
  36  *
  37  * t0, the initial "temperature. Should be set so that the probability
  38  * of accepting a transition to a higher cost configuration is
  39  * initially about 0.8.
  40  *
  41  * ntemps, the number of temperatures to use. Each successive temperature
  42  * is some fraction of the previous temperature.
  43  *
  44  * nmoves_per_temp, the number of configurations to try at each temperature
  45  *
  46  * It is a black art to set ntemps, nmoves_per_temp, and the rate
  47  * at which the temperature drops. Go too fast with too few iterations,
  48  * and the computation falls into a local minimum instead of the
  49  * (desired) global minimum.
  50  */
  51
  52 void clib_anneal (clib_anneal_param_t * p)
  53 {
  54   f64 t;
  55   f64 cost, prev_cost, delta_cost, initial_cost, best_cost;
  56   f64 random_accept, delta_cost_over_t;
  57   f64 total_increase=0.0, average_increase;
  58   u32 i, j;
  59   u32 number_of_increases = 0;
  60   u32 accepted_this_temperature;
  61   u32 best_saves_this_temperature;
  62   int accept;
  63
  64   t = p->initial_temperature;
  65   best_cost = initial_cost = prev_cost = p->anneal_metric (p->opaque);
  66   p->anneal_save_best_configuration (p->opaque);
  67
  68   if (p->flags & CLIB_ANNEAL_VERBOSE)
  69     fformat(stdout, "Initial cost %.2f\n", initial_cost);
  70
  71   for (i = 0; i < p->number_of_temperatures; i++)
  72     {
  73       accepted_this_temperature = 0;
  74       best_saves_this_temperature = 0;
  75
  76       p->anneal_restore_best_configuration (p->opaque);
  77       cost = best_cost;
  78
  79       for (j = 0; j < p->number_of_configurations_per_temperature; j++)
  80         {
  81           p->anneal_new_configuration (p->opaque);
  82           cost = p->anneal_metric (p->opaque);
  83
  84           delta_cost = cost - prev_cost;
  85
  86           /* cost function looks better, accept this move */
  87           if (p->flags & CLIB_ANNEAL_MINIMIZE)
  88               accept = delta_cost < 0.0;
  89           else
  90               accept = delta_cost > 0.0;
  91
  92           if (accept)
  93             {
  94               if (p->flags & CLIB_ANNEAL_MINIMIZE)
  95                 if (cost < best_cost)
  96                   {
  97                     if (p->flags & CLIB_ANNEAL_VERBOSE)
  98                       fformat (stdout, "New best cost %.2f\n", cost);
  99                     best_cost = cost;
 100                     p->anneal_save_best_configuration (p->opaque);
 101                     best_saves_this_temperature++;
 102                   }
 103
 104               accepted_this_temperature++;
 105               prev_cost = cost;
 106               continue;
 107             }
 108
 109           /* cost function worse, keep stats to suggest t0 */
 110           total_increase += (p->flags & CLIB_ANNEAL_MINIMIZE) ?
 111             delta_cost : -delta_cost;
 112
 113           number_of_increases++;
 114
 115           /*
 116            * Accept a higher cost with Pr { e^(-(delta_cost / T)) },
 117            * equivalent to rnd[0,1] < e^(-(delta_cost / T))
 118            *
 119            * AKA, the Boltzmann factor.
 120            */
 121           random_accept = random_f64 (&p->random_seed);
 122
 123           delta_cost_over_t = delta_cost / t;
 124
 125           if (random_accept < exp (-delta_cost_over_t))
 126             {
 127               accepted_this_temperature++;
 128               prev_cost = cost;
 129               continue;
 130             }
 131           p->anneal_restore_previous_configuration (p->opaque);
 132         }
 133
 134       if (p->flags & CLIB_ANNEAL_VERBOSE)
 135         {
 136           fformat (stdout, "Temp %.2f, cost %.2f, accepted %d, bests %d\n", t,
 137                    prev_cost, accepted_this_temperature,
 138                    best_saves_this_temperature);
 139           fformat (stdout, "Improvement %.2f\n", initial_cost - prev_cost);
 140           fformat (stdout, "-------------\n");
 141         }
 142
 143       t = t * p->temperature_step;
 144     }
 145
 146   /*
 147    * Empirically, one wants the probability of accepting a move
 148    * at the initial temperature to be about 0.8.
 149    */
 150   average_increase = total_increase / (f64) number_of_increases;
 151   p->suggested_initial_temperature =
 152     average_increase / 0.22 ; /* 0.22 = -ln (0.8) */
 153
 154   p->final_temperature = t;
 155   p->final_metric = p->anneal_metric (p->opaque);
 156
 157   if (p->flags & CLIB_ANNEAL_VERBOSE)
 158     {
 159       fformat (stdout, "Average cost increase from a bad move: %.2f\n",
 160                average_increase);
 161       fformat (stdout, "Suggested t0 = %.2f\n",
 162                p->suggested_initial_temperature);
 163     }
 164 }