_SIMULATED ANNEALING_
by Michael P. McLaughlin


[LISTING ONE]

/* ANNEAL.C -- Author: Dr. Michael P. McLaughlin
*/
#include <stdio.h>
#include <math.h>

#define ABS(x) ((x<0)?(-(x)):(x))
#define BOOLEAN int
#define TRUE 1
#define FALSE 0
#define FORWARD 1
#define BACK 0

struct cardtype {
   char card[3];
   int code;
   } cards[25];
float ratios[10][2];
int tableau[25],best_tableau[25];
float temperature,best_temperature,t_low,t_min,t_ratio;
long seed,score,best_score,tot_successes,tot_failures,scor,exg_cutoff,
     report_time,report_interval,modulus = 2147483399,step = 1;
int next;
BOOLEAN equilibrium,frozen = FALSE,quench=FALSE,final=FALSE;

/* variables needed to assess equilibrium */
float totscore,totscore2,avg_score,sigma,half_sigma,score_limit;
long bscore,wscore,max_change,nsuccesses,nfailures,scale,
     incount,outcount;
int incount_limit,outcount_limit,succ_min,first_limit,ultimate_limit;

/* poker hands in tableau */
int hand[12][5] = {0,1,2,3,4,      /* rows */
                   5,6,7,8,9,
                   10,11,12,13,14,
                   15,16,17,18,19,
                   20,21,22,23,24,
                   0,5,10,15,20,   /* columns */
                   1,6,11,16,21,
                   2,7,12,17,22,
                   3,8,13,18,23,
                   4,9,14,19,24,
                   0,6,12,18,24,   /* diagonals */
                   4,8,12,16,20};
                /* 0,4,12,20,24     = corners */

long rand1()
/* Get uniform 31-bit integer -- Ref. CACM, June 1988, pg. 742 */
{
   register long k;
   k = seed/52774;
   seed = 40692*(seed-k*52774)-k*3791;
   if (seed < 0)
      seed += modulus;
   return(seed);
}
double uni(z)
int z;
{
   return((double) z*rand1()/modulus);
}
void initialize()
/* Set up entire simulated annealing run. */
{
   char label[20];
   double exgc;
   register i;
   FILE *fp;
   fp = fopen("a.dat","r");
   fscanf(fp,"%f %s\n",&temperature,label);
   fscanf(fp,"%f %s\n",&t_low,label);
   fscanf(fp,"%f %s\n",&t_min,label);
   fscanf(fp,"%f %s\n",&avg_score,label);
   fscanf(fp,"%ld %s\n",&scale,label);
   fscanf(fp,"%f %s\n",&t_ratio,label);
   fscanf(fp,"%f %s\n",&sigma,label);
   fscanf(fp,"%lf %s\n",&exgc,label);
   fscanf(fp,"%d %s\n",&succ_min,label);
   fscanf(fp,"%d %s\n",&incount_limit,label);
   fscanf(fp,"%d %s\n",&outcount_limit,label);
   fscanf(fp,"%d %s\n",&first_limit,label);
   fscanf(fp,"%d %s\n",&ultimate_limit,label);
   fscanf(fp,"%ld %s\n",&report_interval,label);
   fscanf(fp,"%ld %s\n",&seed,label);
   half_sigma = 0.5*sigma; report_time = report_interval;
   score_limit = 20; exg_cutoff = modulus*exgc;
   for (i=0;i<10;i++)
      fscanf(fp,"%f %f\n",&ratios[i][0],&ratios[i][1]);
   for (i=0;i<25;i++) {
      fscanf(fp,"%d %s\n",&cards[i].code,cards[i].card);
      tableau[i] = cards[i].code;
   }
   fclose(fp);
   printf("                   TOTAL     TOTAL  CURRENT AVERAGE\n");
   printf("STEP TEMPERATURE SUCCESSES FAILURES  SCORE   SCORE ");
   printf(" SIGMA\n\n");
}
void init()
/* set up for new temperature */
{
   nsuccesses = 0; nfailures = 0; equilibrium = FALSE; incount = 0;
   outcount = 0; bscore = 0; wscore = 100000; max_change = 0;
   totscore = 0; totscore2 = 0;
   if (temperature < t_min)
      frozen = TRUE;
}
void get_new_temperature()
{
   if (temperature < ratios[next][0]) {
      t_ratio = ratios[next][1];
      next++;
   }
   temperature = t_ratio*temperature;
}
void check_equilibrium()
/* Determine whether equilibrium has been reached. */
{
   if ((nsuccesses+nfailures) > ultimate_limit)
      equilibrium = TRUE;
   else if (nsuccesses >= succ_min) {
      if (incount > incount_limit)
         equilibrium = TRUE;
      else {
         if (outcount > outcount_limit) {
            if (nsuccesses > first_limit)
               equilibrium = TRUE;
            else {
               incount = 0;
               outcount = 0;
            }
         }
      }
   }
}
void update()
/* Compute statistics, etc. */
{
   float ascore,s;
   if (nsuccesses > 0) {
      ascore = totscore/nsuccesses;
      avg_score = ascore+scale;
      s = totscore2/nsuccesses-ascore*ascore;
      if (s > 0.0) {
         sigma = sqrt(s); half_sigma = 0.5*sigma;
         score_limit = avg_score-half_sigma;
      }
   }
   if ((temperature < t_low)&&((bscore-wscore) == max_change))
      frozen = TRUE;
}
void selectwr(array,mode,nchoices)
/* Select from array without replacement. */
int *array,mode,nchoices;
{
   int i,temp,size,index;
   size = mode==1 ? 12 : 25;
   for (i=0;i<nchoices;i++) {
      index = uni(size);
      temp = array[--size];
      array[size] = array[index];
      array[index] = temp;
   }
}
void reconfigure(direction)
/* If direction is FORWARD, make a new move; if it is BACK, restore
   the tableau to its previous configuration. */
int direction;
{
   static long partition[6] = {0,0,715827799,1296031795,1766307855,
                               2147483400};
   static int hindex[12] = {0,1,2,3,4,5,6,7,8,9,10,11};
   static int cindex[25] = {0,1,2,3,4,5,6,7,8,9,10,11,12,13,
                            14,15,16,17,18,19,20,21,22,23,24};
   static int overlap[66] = {-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,0,5,10,15,
                             20,1,6,11,16,21,-1,2,7,12,17,22,-1,-1,3,
                             8,13,18,23,-1,-1,-1,4,9,14,19,24,-1,-1,-1,
                             -1,0,6,12,18,24,0,6,12,18,24,4,8,12,16,20,
                             20,16,12,8,4,12};
   static int mode,nchoices,common,last;   /* remember changes */
   int temp,hi,lo,i,j,k;
   if (direction == FORWARD) {
      if (rand1() < exg_cutoff) {   /* giant step */
         mode = 1; nchoices = 2;
         selectwr(hindex,mode,nchoices);
         hi = hindex[11]; lo = hindex[10];   /* order hand indices */
         if (hi < lo) {
            hi = lo;
            lo = hindex[11];
         }
         common = overlap[hi*(hi-1)/2+lo];   /* triangular matrix */
      }
      else {   /* baby step */
         mode = 0; nchoices = 2;
         rand1();   /* How many cards to rotate? */
         while (seed > partition[nchoices])
            nchoices++;
         selectwr(cindex,mode,nchoices);
         temp = tableau[cindex[24]];   /* rotate forward */
         for (i=1,j=24;i<nchoices;i++,j--)
            tableau[cindex[j]] = tableau[cindex[j-1]];
         tableau[cindex[j]] = temp;
         last = j;
      }
   }
   if (mode == 1) {   /* swapping hands is symmetrical */
      register first,second,c;
      first = hindex[11]; second = hindex[10];
      k = (common<0) ? 5 : 4;
      for (c=0,i=0,j=0;c<k;c++,i++,j++) {
         if (hand[first][i] == common)
            i++;
         if (hand[second][j] == common)
            j++;
         temp = tableau[hand[first][i]];
         tableau[hand[first][i]] = tableau[hand[second][j]];
         tableau[hand[second][j]] = temp;
      }
   }
   else if (direction == BACK) {   /* rotation is not */
      temp = tableau[cindex[last]];
      for (i=1,j=last;i<nchoices;i++,j++)
         tableau[cindex[j]] = tableau[cindex[j+1]];
      tableau[cindex[24]] = temp;
   }
}
long get_score(crds)
/* Return score of one hand. */
int *crds;
{
   static long scores[18] = {0,1,1,20,1,9,20,1760,1,9,9,
                             293,20,293,1760,108,215,27456};
   int flush,i,index = 0;
   flush = crds[0]&crds[1]&crds[2]&crds[3]&crds[4]&0xff00;
   for (i=0;i<5;i++)
      crds[i] = crds[i]&0xff;   /* ignore suits */
   {
      register b,k,t;   /* sort cards from high to low */
      for (b=4;b>0;b--)
         for (k=0;k<b;k++)
            if (crds[k] < crds[k+1]) {
               t = crds[k];
               crds[k] = crds[k+1];
               crds[k+1] = t;
            }
   }
   if (!flush) {   /* multiple cards ? */
      if (crds[0] == crds[1]) index += 8;
      if (crds[1] == crds[2]) index += 4;
      if (crds[2] == crds[3]) index += 2;
      if (crds[3] == crds[4]) index += 1;
   }
   if (!index) {   /* straight ? */
      if (((crds[0]-crds[4]) == 4)||((crds[0] == 14)&&(crds[1] == 5)))
         index = 15;
      if (flush)
         index = (index == 15) ? 17 : 16;
   }
   return(scores[index]);
}
long evaluate1()
/* Return total score of tableau (no corners). */
{
   int temp[5],h,c;
   long scr = 0;
   for (h=0;h<12;h++) {   /* h<13 for corners */
      for (c=0;c<5;c++)
         temp[c] = tableau[hand[h][c]];
      scr += get_score(temp);
   }
   return(scr);
}
void decide()
/* Either accept new configuration or restore old configuration */
{
   float pdt1,pdt2,sscor;
   long s;
   BOOLEAN acceptable = FALSE;
   scor = evaluate1();
   if (scor >= score) {
      if (scor > best_score) {
         register i;
         best_score = scor;
         best_temperature = temperature;
         for (i=0;i<25;i++)
            best_tableau[i] = tableau[i];
      }
      acceptable = TRUE;
   }
   else if (temperature > 0.0) {   /* uphill movement? */
      if (uni(1) < exp((scor-score)/temperature))   /* Equation 2 */
         acceptable = TRUE;
   }
   if (acceptable) {   /* statistics, etc. */
      s = ABS(score-scor);
      if (s > max_change)
         max_change = s;
      score = scor;
      sscor = scor-scale;   /* to aid precision of sigma */
      totscore += sscor;
      totscore2 += sscor*sscor;
      nsuccesses++;
      if (ABS(avg_score-score) < half_sigma)
         incount++;
      else
         outcount++;
      if (score > bscore)   /* maximization */
         bscore = score;
      else if (score < wscore)
         wscore = score;
   }
   else {   /* unacceptable */
      reconfigure(BACK);
      nfailures++;
   }
}
void report()
{
   tot_successes += nsuccesses; tot_failures += nfailures;
   printf("%3ld %10.1f %9ld %9ld %7ld %7.0f  %5.1f\n",step,temperature,
          tot_successes,tot_failures,score,avg_score,sigma);
   report_time += report_interval;
   if (frozen) {
      int temp[25],k,kk;
      register i,j;
      BOOLEAN ok;
      if (final)
         printf("\nFINAL -- ");
      else if (quench)
         printf("\nQUENCH -- ");
      else
         printf("\nFROZEN -- ");
      printf("BEST SCORE = %5ld  BEST TEMPERATURE = %5.1f\n\n",
              best_score,best_temperature);
      for (i=0;i<25;i++) {
         k = best_tableau[i];
         j = 0;
         ok = FALSE;
         while (!ok) {
            if (cards[j].code == k) {
               temp[i] = j;
               ok = TRUE;
            }
            else j++;
         }
      }
      for (i=0;i<25;i+=5) {
         for (j=i;j<(i+5);j++) {
            k = 4-strlen(cards[temp[j]].card);
            for (kk=0;kk<k;kk++) printf(" ");
            printf("%s",cards[temp[j]].card);
         }
         printf("\n");
      }
      printf("\n");
   }
}
void try_new_temperature()
{
   init();
   while (!equilibrium) {
      reconfigure(FORWARD);
      decide();
      check_equilibrium();
   }
   update();
   if (!quench) {   /* ratchet */
      while ((float) score < score_limit) {   /* maximization */
         reconfigure(FORWARD);
         decide();
      }
   }
}
main()
{
   initialize();
   while (!frozen) {
      try_new_temperature();
      update();
      if ((step == report_time)||(frozen))
         report();
      step++;
      if (temperature > 0.0)
         get_new_temperature();
   }
/* Quench frozen configuration. */
   temperature = 0;
   quench = TRUE;
   try_new_temperature();
   update();
   report();
   step++;
/* Quench best configuration (rarely useful). */
   {
      register i;
      final = TRUE;
      for (i=0;i<25;i++)
         tableau[i] = best_tableau[i];
   }
   try_new_temperature();
   update();
   report();
   printf("SEED = %ld\n",seed);   /* seed of next run */
}



[LISTING TWO]

     if ((nsuccesses+nfailures) > ULTIMATE_LIMIT)
       equilibrium = true;
     else if (nsuccesses >= SUCC_MIN) {
       if (incount > INCOUNT_LIMIT)
        equilibrium = true;
       else {
        if (outcount > OUTCOUNT_LIMIT) {
         if (nsuccesses > FIRST_LIMIT)
          equilibrium = true;
         else {
          incount = 0;
          outcount = 0;
         }
        }
       }
     }


[LISTING THREE]

2150 temperature
1.5 t_low
0.1 t_min
40 avg_score
4400 scale
0.7 t_ratio
150 sigma
0.2 exg_cutoff
200 succ_min
250 incount_limit
400 outcount_limit
2700 first_limit
10000 ultimate_limit
1 report_interval
1234567890 seed
360 0.8
215 0.7
85 0.8
60 0.9
30 0.95
15 0.9
7 0.8
3 0.7
0 0.7
0 0.7
1032 8H
2061 KS
270 AC
526 AD
522 10D
1029 5H
265 9C
1035 JH
1037 KH
525 KD
515 3D
263 7C
2058 10S
2062 AS
518 6D
1036 QH
267 JC
259 3C
2053 5S
269 KC
520 8D
1028 4H
1038 AH
519 7D
1026 2H



[EXAMPLE 1] 

Equation 1:  N /su/hi//N /su/lo/ = exp((E lo-E hi)/kT)     

Equation 2:  Prob(accept worse score) = exp((S /su/lo/ -S /su/hi/)/T)   


Example 1: Equations that describe the Boltzman distribution