siman.cpp

#include "qscr.h"
#include "siman.h"
#include "linrdist.h"
#include "util.h"
#include <iostream>
#include <cmath>
#include <stdlib.h>
#include <cassert>



SimAnneal::SimAnneal(Function *ff, matrix simplex, 
                     double temptr, unsigned it, double tol, long s,
                     unsigned steps) {
  unsigned n = (*ff).getdim();
  assert(simplex.row()==int(n+1) && simplex.column()==int(n));
  p = simplex;
  funk = ff;
  ndim = n;
  tt = -temptr;
  seed = s;
  ftol = tol;
  iter = it;
  static_iter = it;
  nosteps = steps;

  psum.newsize(ndim);
  calc_psum();

  pb.newsize(ndim);
  pb = p.getrow(0).trans();
  yb = (*funk)(pb);

  y.newsize(ndim);
  y = pval(p);  
}


matrix SimAnneal::pval(matrix const& pp) {
  assert(pp.row()==int(ndim+1) && pp.column()==int(ndim));
  matrix res(ndim+1);
  for (unsigned i=0; i<ndim+1; i++) {
    res.set((*funk)(pp.getrow(i).trans()),i);
  }
  return res;
}


void SimAnneal::calc_psum()
{
 
  // ndim = p.rows
  // ndim+1 = p.cols
  for (unsigned j=0;j<ndim;j++) {
    double sum = 0.0;
    for (unsigned i=0;i<ndim+1;i++)
      sum += p.get(i,j);
    psum.set(sum,j);
  }
}


double SimAnneal::amotsa(int ihi, double &yhi, 
                         double const fac) {
  // Extrapolates by a factor fac through the face of the simplex across from the high point,
  // tries it, and replaces the high point if the new point is better.
  
  unsigned j;
  double fac1, fac2, yflu, ytry;
  matrix ptry(ndim);

  fac1 = (1.0-fac)/ndim;
  fac2 = fac1-fac;
  for (j=0; j<ndim; j++) {
    double val = psum.get(j)*fac1 - p.get(ihi,j)*fac2;
    ptry.set(val,j);
  }
  //  cerr << "ptry="; ptry.print();
  ytry = (*funk)(ptry);
  //  cerr << "f(ptry)=ytry=" << ytry << endl;
  if (ytry<=yb) {
    for (j=0; j<ndim; j++) { // Save the best ever
      pb.set(ptry.get(j),j);
    }
    yb = ytry; 
  }
  // We added a thermal fluctuation to all the current vertices, but we subtract
  // it here, so as to give the simplex a thermal Brownian motion: It "likes" to
  // accept any suggested change.
  yflu = ytry-tt*log(ran1(&seed)); 
  if (yflu<yhi) {
    y.set(ytry,ihi);
    yhi=yflu;
    for (j=0; j<ndim; j++) {
      double val = ptry.get(j)-p.get(ihi,j);
      psum.set(psum.get(j) + val,j);
      p.set(ptry.get(j),ihi,j);
    }
  }
  return yflu;
}


/*
  Multidimensional minimization of the function funk(x) where x[1..ndim] is a 
  vector in ndim dimensions, by simulated annealing combined with the downhill
  simplex method of Nelder and Mead. The input matrix p[1..ndim+1][1..ndim] has
  ndim+1 rows, each an ndim-dimensional vector which is a vertex of the starting
  simplex. Also input are the following: the vector y[1..ndim+1] whose components
  must be pre-initialized to the values of funk evaluated at the ndim+1 vertices
  (rows) of p; ftol, the fractional convergence tolerance to be archived in the
  function value for an early return; iter, and temptr. The routine makes iter
  function evaulations at an annealing temperature temptr, the returns. You
  should then decrease temptr according to your annealing schedule, reset iter
  and call the routine again (leaving other arguments unaltered between calls).
  If iter is returned with a postive value, then early convergence and return
  occurred . If you initialize yb to a very large value on the first call, 
  then yb and pb[1..ndim] will subsequently return the best function value
  and point ever encountered (even if it is no longer a point in the simplex).
*/

void SimAnneal::amebsa() {

  unsigned ihi,ilo,mpts=ndim+1;
  double rtol, swap, ylo, ynhi, ysave, yt, ytry, yhi;
  matrix psum(ndim);
  unsigned cnt = 1;
  std::string st;

   // Determine which point is the highest (worst), next-highest and lowest (best).
  calc_psum();
  for (;;) {

    st = "Min-value so far=" + numStr(yb) + 
      " found at ( ";
    for (unsigned i=0; i<ndim; i++) {
      st += numStr(pb.get(i)) + " ";
    }
    st += ")";
    qmvaddstr(22,5, st.c_str()); qclrtoeol();

    switch(cnt) {
    case 1: 
      qmvaddstr(24,5, "---");
      break;
    case 2:
      qmvaddstr(24,5, "///");
      break;
    case 3: 
      qmvaddstr(24,5, "|||");
      break;
    case 4: 
      qmvaddstr(24,5, "\\\\\\");
      break;
    case 5: 
      qmvaddstr(24,5, "---");
      break;
    case 6: 
      qmvaddstr(24,5, "\\\\\\");
      break;
    case 7: 
      qmvaddstr(24,5, "|||");
      break;
    case 8: 
      qmvaddstr(24,5, "///");
      cnt = 0;
      break;
    }

    cnt++;
    
//     cerr << "Min-value so far=" << yb 
//       << " found at (" << pb.get(0) << ", " << pb.get(1) << ", " 
//       << pb.get(2) << ")" 
//       << endl;


    
    
    qmvaddstr(6,55, "Current simplex:");
    for (int i=0; i<p.row(); i++) {
      qmvaddstr(7+i,56, p.row2string(i)); qclrtoeol();
    }
    st = "Temperature: " + numStr(-tt); 
    qmvaddstr(12,55, st.c_str()); qclrtoeol();
    st = "Iteration: " + numStr(iter);
    qmvaddstr(13,55, st.c_str()); qclrtoeol();

//     cerr << "Current simplex="; p.print();
//     cerr << "values here="; y.print();

    qmove(24,8);

    ilo=0;
    ihi=1;
    // Whenever we "look at" a vertex, it gets a random thermal fluctuation.
    ynhi=ylo=y.get(0)+tt*log(ran1(&seed));
    yhi=y.get(1)+tt*log(ran1(&seed));
    if (ylo>yhi) {
      ihi=0;
      ilo=1;
      swap = ynhi; // ??
      ynhi=yhi;
      yhi=ylo;
      ylo=ynhi;  // OBS: skulle dette være den gamle værdi af ynhi?!?!?
      //      ylo = swap;
    }
    for (unsigned i=2; i<mpts; i++) {
      yt=y.get(i)+tt*log(ran1(&seed));      
      if (yt<=ylo) {
        ilo=i;
        ylo=yt;
      }
      if (yt>yhi) {
        ynhi=yhi;
        ihi=i;
        yhi=yt;
      } else if (yt>ynhi) {
        ynhi = yt;
      } 
    }

    // Compute the fractional range from highest to lowest and return if satisfactory.
    rtol = 2.0*fabs(yhi-ylo)/(fabs(yhi)+fabs(ylo));
    if (rtol<ftol || iter<0) {
      // If returning put best point and value in first slot
      swap=y.get(0);
      y.set(y.get(ilo),0);
      y.set(swap,ilo);
      for (unsigned n=0;n<ndim;n++) {
        swap = p.get(0,n);
        p.set(p.get(ilo,n),0,n);
        p.set(swap,ilo,n);
      }
      break;
    }
  
    iter -= 2;
    // Begin a new iteration. First extrapolate by a factor -1 through the face of the simplex
    // across from the high point, i.e. reflect the simplex from the high point.
    ytry = amotsa(ihi,yhi,-1.0);
    if (ytry<=ylo) {
      // Gives a result better than the best point, so try an additional extrapolation by a
      // factor 2
      ytry = amotsa(ihi,yhi,2.0);
      
    } else if (ytry>=ynhi) {
      // The reflected point is worse than the second-highest, so look for an intermediate
      // lower point, i.e., do a one-dimensional contraction.
      ysave=yhi;
      ytry = amotsa(ihi,yhi,0.5);
      
      if (ytry>=ysave) {
        // Can't seem to get rid of that high point. Better contract around the lowest (best) point.
        for (unsigned i=0;i<mpts;i++) {
          if (i != ilo) {
            for (unsigned j=0;j<ndim;j++) {
              double val = 0.5*(p.get(i,j)+p.get(ilo,j));
              psum.set(val,j);
              p.set(val,i,j);
            }
            y.set((*funk)(psum),i);
          }
        }
        iter -= ndim;
        calc_psum(); // Recompute psum    
      }
    } else ++(iter);  // Correct the evaluation count.
  }
  
}
 
void SimAnneal::optimize() {
  amebsa();
  iter = static_iter;
}



void SimAnneal::anneal() {
  double step = tt/10;
  double origtt = tt;
  double curtt = tt;
  for (unsigned i=0; i<nosteps; i++) {
    tt = curtt;
    p = newsimplex(pb);
    optimize();
    curtt -= step;
  }
  tt = origtt;
}




/////////////// NON-member function follows:


matrix newsimplex(matrix const& x, double epsilon) {
  // P = (pb0,pb1,pb2) + lambda*e_i

  unsigned nn = x.row();

  matrix p(nn+1,nn);
  for (unsigned i=0; i<nn; i++) {
     p.set(x.get(i),0,i);
  }

  for (unsigned j=1; j<=nn; j++) {
    for (unsigned i=0; i<nn; i++) {
      p.set(x.get(i),j,i);
    }
    p.set(p.get(j,j-1)+epsilon,j,j-1);
  }

  return p;
}



#ifdef TEST

// Functions to minimize:

// Assumes input is a 2x1-matrix. f1(x,y) = x^2+y^2
double f2(matrix const& x) {
  double xx = (x.get(0)-1);
  double xy = (x.get(1)-1);
  return xx*xx+xy*xy;
}

// Assumes input is a singleton matrix. Function with many local minima. Global minimum somewhere in [1.2, 1.5]
double f1(matrix const& x) {
  double xx = x.get(0);
  return exp(-pow((xx-1.0),2.0))*sin(8*xx);
}


int main(int argc, char *argv[]) {
  try {

    qinit();

//     matrix p(3,2);
//     p.set(6,0,0);  p.set(1,0,1);
//     p.set(2,1,0);  p.set(2,1,1);
//     p.set(3,2,0);  p.set(3,2,1);

//     SimAnneal sa(f2,p);


    Function FF(f1,1);
    matrix simplex(2);  // Sets simplex to the interval [6,7]
    simplex.set(6,0); simplex.set(7,1);

    /*
    SimAnneal(double (*f)(matrix), matrix simplex, 
            double temptr=0.01, int it=100, double tol=0.001, long s=12345);
    */
    std::cerr << std::endl;

    matrix minP_anneal, minP_simplex;
    double minV_anneal, minV_simplex;

    { // Demonstrates simulated annealing (SA)
      SimAnneal sa(&FF, simplex, 0.01, 500);
      sa.anneal();
      minP_anneal = sa.get_pb();
      minV_anneal = sa.get_yb();
    } // Destruct sa

    std::cerr << "balblablalblaba" << std::endl;
    { // Demonstrates simplex algorithm (temperature=0.0 in SA-algorithm)
      SimAnneal sa(&FF, simplex, 0.0, 2000);
      sa.optimize();
      minP_simplex = sa.get_pb();
      minV_simplex = sa.get_yb();
    } // Destruct sa
    

    std::string line = "";
    for (unsigned i=0; i<50; i++) 
      line += "-";
    line += "\n";

    qend();


    std::cerr << line;
    std::cerr << "Minimum found with SA-method " <<
      (minP_anneal.trans()).row2string(0) << std::endl; 
    std::cerr << "With function value " << minV_anneal << std::endl << std::endl;

    std::cerr << line;
    std::cerr << "Minimum found with Simplex-method " << 
      (minP_simplex.trans()).row2string(0) << std::endl; 
    std::cerr << "With function value " << minV_simplex << std::endl << std::endl;

    std::cerr << line;

    } catch(std::string const& err) {
    std::cerr << err << std::endl;
    return 1;
  }
  return 0;
}



#endif//TEST

---