simann.cxx

// simanneal.c++   Implementation of a General Purpose Simulated Annealing Class
//                      (c) Copyright 1994, Everett F. Carter Jr.
//                      Permission is granted by the author to use
//			this software for any application provided this
//			copyright notice is preserved.

// Updated to 2003 C++ standard by Shawn Waldon in 2014

/* Uses Cauchy training		*/

static const char rcsid[] = "@(#)simann.c++	1.4 15:02:21 7/25/94   EFC";

#include <cstdlib>
#include <cmath>

#include "simann.hpp"

#ifndef PI2
#define PI2 (PI / 2.0)
#endif

SimAnneal::SimAnneal(const CostFunction& f, const int d)
    : func(f),
      dimension(d),
      maxit(400),
      ddwell(20),
      dt(0.1),
      c_jump(100.0),
      rrange(PI2),
      K(1.0),
      rho(0.5),
      t0(0.0),
      tscale(0.1)

{

  x = new double[dimension];
  xnew = new double[dimension];
  xbest = new double[dimension];
  y = ybest = HUGE;

  if (x == NULL || xnew == NULL || xbest == NULL)
    err = -1;
  else
    err = 0;
}

int SimAnneal::set_up(CostFunction f, const int d)
{

  dimension = d;

  func = f;

  x = new double[dimension];
  xnew = new double[dimension];
  xbest = new double[dimension];
  y = ybest = HUGE;

  if (x == NULL || xnew == NULL || xbest == NULL)
    err = -1;
  else
    err = 0;

  return err;
}

/* increase the temperature until the system "melts" */
double SimAnneal::melt(const int iters)
{
  int i, j, ok = 0;
  double xc, ynew, t, cold, c = 0.0;

  int n = iters;
  if (n < 1) n = maxit;

  t = t0;

  for (i = 0; i < n; i++) {
    if (i > 0 && c > 0.0) {
      cold = c;
      ok = 1;
    }

    t += dt;

    for (j = 0; j < dimension; j++) {
      xc = rho * t * tan(uniform.number(-rrange, rrange));
      x[j] += xc;
    }

    equilibrate(t, ddwell);

    /* "energy" */
    ynew = func(x);
    c = ynew - y;

    if (c < 0.0 && ynew < ybest) {
      for (j = 0; j < dimension; j++) xbest[j] = x[j];
      ybest = ynew;
    }

    y = ynew;

    if (ok && c > (c_jump * cold)) /* phase transition */
      break;
  }

  return t0 = t;
}

/* iterate a few times at the present temperature */
/* to get to thermal equilibrium */
int SimAnneal::equilibrate(const double t, const int n)
{
  int i, j, equil = 0;
  double xc, ynew, c, delta;
  double* xtmp;
  //	double p;

  delta = 1.0;
  for (i = 0; i < n; i++) {
    for (j = 0; j < dimension; j++) {
      xc = rho * t * tan(uniform.number(-rrange, rrange));
      xnew[j] = x[j] + xc;
    }
    /* "energy" */
    ynew = func(xnew);
    c = ynew - y;

    if (c < 0.0) /* keep xnew if energy is reduced */
    {
      xtmp = x;
      x = xnew;
      xnew = xtmp;

      y = ynew;
      if (y < ybest) {
        for (j = 0; j < dimension; j++) xbest[j] = x[j];
        ybest = y;
      }

      delta = fabs(c);
      delta = (y != 0.0) ? delta / y : (ynew != 0.0) ? delta / ynew : delta;

      /* equilibrium is defined as a 10% or smaller change
         in 10 iterations */
      if (delta < 0.10)
        equil++;
      else
        equil = 0;

    } else {
      /* keep xnew with probability, p, if ynew is increased */
      /*
                              p = exp( - (ynew - y) / (K * t) );

                              if ( p > uniform.number(0.0, 1.0) )
                              {
                                      xtmp = x;
                                      x = xnew;
                                      xnew = xtmp;
                                      y = ynew;
                                      equil = 0;
                              }
                              else
      */

      equil++;
    }

    if (equil > 9) break;
  }

  return i + 1;
}

/* cool the system with annealing */
double SimAnneal::anneal(const int iters)
{
  int i, j;
  double xc, p, ynew, t, c, dt, told;
  double* xtmp;

  int n = iters;
  if (n < 1) n = maxit;

  equilibrate(t0, 10 * ddwell);

  told = t0;
  for (i = 0; i < n; i++) {
    t = t0 / (1.0 + i * tscale);
    dt = t - told;
    told = t;

    equilibrate(t, ddwell);

    for (j = 0; j < dimension; j++) {
      xc = rho * t * tan(uniform.number(-rrange, rrange));
      xnew[j] = x[j] + xc;
    }
    /* "energy" */
    ynew = func(xnew);
    c = ynew - y;

    if (ynew <= y) /* keep xnew if energy is reduced */
    {
      xtmp = x;
      x = xnew;
      xnew = xtmp;
      y = ynew;

      if (y < ybest) {
        for (j = 0; j < dimension; j++) xbest[j] = x[j];
        ybest = y;
      }
      continue;
    } else {
      /* keep xnew with probability, p, if ynew is increased */
      p = exp(-(ynew - y) / (K * t));

      if (p > uniform.number(0.0, 1.0)) {
        xtmp = x;
        x = xnew;
        xnew = xtmp;
        y = ynew;
      }
    }
  }

  equilibrate(t, 10 * ddwell);

  t0 = t;

  return t;
}

void SimAnneal::initial(const double* xi)
{
  for (int k = 0; k < dimension; k++) x[k] = xi[k];
}

void SimAnneal::current(double* xc)
{
  for (int k = 0; k < dimension; k++) xc[k] = x[k];
}

void SimAnneal::optimum(double* xb)
{
  for (int k = 0; k < dimension; k++) xb[k] = xbest[k];
}