[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] Simulated Annealing

From: John Gehman
Subject: Re: [Help-gsl] Simulated Annealing
Date: Thu, 02 Nov 2006 07:21:31 +1100

On 31/10/2006, at 10:42 PM, Fred J. wrote:

John Gehman <address@hidden> wrote:

I don't have the documentation handy at the moment, so just some general
guidance ...

What you've done is a systematic search of all parameter space. And even
though you've declared double type, you appear to be considering only
integer values for x and y. Maybe that's good enough, but if you consider
real double number, you might get a better extremes.

Simulated Annealing will help you if your "energy surface" (i.e. the
surface defined by z values for give (x,y) is relatively discontinuous, i.e. very very flat in most areas with narrow bands of steep descent, or strafed with lots of local minima/maxima which make it difficult to find
the global minimum/maximum.

On the other hand if your surface is reasonably smooth, and if you have an analytic expression for the derivative of your function with respect to
both x and y, (or even if not), you could look into nonlinear least
squares fitting routines with or without derivatives, (or use the
empirical derivative-finding routines). Alternatively you can use the
minimizer functions which iteratively bracket the minimum until a
satisfactory minimum is found.


> Hi
> I have a routine
> double myRou ( double x, double y ) { ... };
> the way I use it is run it with many combinations of x,y and find which
> x,y combination gives the lowest or highest retrun. so I do
> vector res;
> for(double i = 1; i < 1000; ++i) {
> for( double j = 1; j < 1000; ++j) {
> res.push_back( myRou( i, j );
> }
> }
> then I go to find the maximum or minimum values in res.
> with out the loop which ends up taking the whole day.
> is there some algorithm that can help me find the global min/max.
> with out the looping, I was reading about gsl Simulated Annealing
> example.html#Trivial-example > but realy could not understand how to apply it to my function. if someone
> could please show a code on how to use myRou.
> thanks
G'day John

thanks for you help.
I was able to compile and run the following samle code in in one dimensional cartesian space, but I am woundering how can I use it to work on 2-dim configuratrion.

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <gsl/gsl_siman.h>
/* set up parameters for this simulated annealing run */
/* how many points do we try before stepping */
#define N_TRIES 200            
/* how many iterations for each T? */
#define ITERS_FIXED_T 10       
/* max step size in random walk */
#define STEP_SIZE 10           
/* Boltzmann constant */
#define K 1.0                  
/* initial temperature */
#define T_INITIAL 0.002        
/* damping factor for temperature */
#define MU_T 1.005             
#define T_MIN 2.0e-6
gsl_siman_params_t params
/* now some functions to test in one dimension */
double E1(void *xp)
  double x = * ((double *) xp);
  return exp(-pow((x-1.0),2.0))*sin(8*x);
double M1(void *xp, void *yp)
  double x = *((double *) xp);
  double y = *((double *) yp);
  return fabs(x - y);
void S1(const gsl_rng * r, void *xp, double step_size)
  double old_x = *((double *) xp);
  double new_x;
  double u = gsl_rng_uniform(r);
  new_x = u * 2 * step_size - step_size + old_x;
  memcpy(xp, &new_x, sizeof(new_x));
void P1(void *xp)
  printf ("%12g", *((double *) xp));
main(int argc, char *argv[])
  const gsl_rng_type * T;
  gsl_rng *  r;
  double x_initial = 15.5;
  T = gsl_rng_default;
  r = gsl_rng_alloc(T);
  gsl_siman_solve(r,        // const gsl_rng *, for steps generation
          &x_initial,     // void *, points to starting configuration and best results once finished           E1,         // gsl_siman_Efunc_t, specifies the space to search
          S1,         // gsl_siman_step_t, for steps generation
          P1,        // gsl_siman_print_t, if NULL (no print out), else stdout
          NULL,     // gsl_siman_copy_t, copyfunc,
          NULL,     // gsl_siman_copy_construct_t, copy_constructor
          NULL,     // gsl_siman_destroy_t, destructor
          sizeof(double), // size_t, fixed size "updating configuration mode", zero in the variable-size mode
          params    // gsl_siman_params_t,
  return 0;

Check out the New Yahoo! Mail - Fire up a more powerful email and get things done faster.

Hi Fred,

My implementation of SImulated Annealing was rather involved; all the details of my problem are held in a C++ class object Paar, passed to the subroutine as P below, and I then have another class object (called SimanObject below) which knows (among other things) how to perform simulated annealing on P. This function is below. I wanted to do a few things that the stock gsl simulated annealing routine didn't provide (actually the authors obviously thought of it, but commented it out for distribution; I like to catch the run progress details to a file, saout), so I actually adapted the code as distributed and incorporated everything explicitly into my code below. Hopefully between this and having a look at the original gsl code (siman.c I think it was?), you'll be able to sort out what you need? P.unknown (and all of its copies) is a gsl_vector which can hold any number of unknowns you wish. If you decide to follow my code below more closely than you probably need to for just two unknowns, be very careful with your equivalent of the Paar copy constructors and destructors; be sure to copy gsl_vector elements with gsl_vector_memcpy, otherwise you're just copying the pointer, and all nominal copies then refer back to the the same data. Also be sure to gsl_vector_free when destroying objects, else you'll get crippling memory leaks.

Best of luck,

int SimanObject::Unknown_siman( Paar &P ) {

        const gsl_rng_type * Tr;
        gsl_rng * r;


        Tr = gsl_rng_taus2;
        r = gsl_rng_alloc(Tr);

        ofstream saout; 
        char saname[179];
        strcat(saname,".sa");, ios::out);

        sprintf(P.str,"#%4s %7s %12s %10s %4s %4s %4s\n",
        saout << P.str;

        Paar x(P), new_x(P), best_x(P);
        double E, new_E, best_E;
        int i, done;
        double T, u;
        int n_evals = 0, n_iter = 0, n_accepts, n_rejects, n_eless;

        E = P.calc_state_dataonly_nosim();
        best_E = E;
        double magicTemp;

        for (int ntry = 0; ntry < P.ntries; ntry++) {

                int iterT = int( gvg(P.iterT,ntry) );   
                double boltzmann = gvg(P.boltzmann,ntry);
                double maxstep = gvg(P.maxstep,ntry);
                double attentemp = gvg(P.attentemp,ntry);
                double mintemp = gvg(P.mintemp,ntry);
                double T = gvg(P.inittemp,ntry);
                if ( abs(T) < 1E-12 ) T = magicTemp;
                if ( abs(attentemp) < 12E-12 ) attentemp = gvg(P.attentemp,0);
                if ( abs(maxstep) < 12E-12 ) maxstep = gvg(P.maxstep,0);
                if ( abs(boltzmann) < 12E-12 ) boltzmann = gvg(P.boltzmann,0);
                if ( abs(mintemp) < 12E-12 ) mintemp = gvg(P.mintemp,0);
                if ( iterT == 0 ) iterT = int ( gvg(P.iterT,0) );

sprintf(P.str,"\n%5d/%d B %8.2e S %8.2e Temp: %8.2e-%8.2e by %11.5e\n",
                P.lfile << P.str;
                cout << P.str;
                done = 0;

                while (!done) {
                        n_accepts = 0;
                        n_rejects = 0;
                        n_eless = 0;
                for (i = 0; i < iterT; ++i) {
                                // take a step
                                for (int lm=0; lm<; lm++) {
                                        u = gsl_rng_uniform(r);
                                        gsl_vector_set( new_x.unknown, lm, u * 
2.0 * maxstep -
                                                maxstep + gvg( new_x.unknown, 
lm) );
                                new_E = new_x.calc_state_dataonly_nosim();
                                if ( !isfinite(new_E) ) {
                                        cout << endl;
                                        cout << "BANG. Chi2 is not finite at temp " 
<< T
                                                << endl;
                                if(new_E <= best_E){
                                        magicTemp = T;
                                // keep track of Ef() evaluations
                                //now take the crucial step: see if the new 
point is accepted
                                //or not, as determined by the boltzman 
probability. If step is
                                //is taken, write it to the .la file together 
with the
                                //Lagrange multiplier configuration.
                                if (new_E < E) {
                        // take a step
                                        E = new_E;
                                } else if
                                        (gsl_rng_uniform(r) < exp 
(-(new_E-E)/(boltzmann*T)) ) {
                                        // also take a step */
                                        E = new_E;
                                } else {
                        //print report
                        char format[79];
                        strcpy(format,"%5d %7d %12g %10.8f %4.2f %4.2f %4.2f ");
                                100.0*double(n_rejects)/double(iterT) );
                        saout << P.str;
                        for (int i=0; i<; i++) {
                                sprintf(P.str,"%10.7le ", gvg(x.unknown,i) );
                                saout << P.str;
                        saout << endl;
                        // Some updates to screen
                        if ( double(n_iter)/5.0 - (n_iter/5) < 1E-6 ) {
                                for (int bs = 0; bs<62; bs++ ) cout << "\b";
                                printf("Temp %10.8f Current E %14.6g Best E 
                // apply the cooling schedule to the temperature
                        T /= attentemp;
                        if (T < mintemp) done = 1;
        //at the end, copy the result onto the initial point

        cout << endl;
        cout << "CHI2   " <<      P.chi2 << endl;

        return 0;


Dr John Gehman (address@hidden)
Research Fellow; Separovic Lab
School of Chemistry
University of Melbourne (Australia)

reply via email to

[Prev in Thread] Current Thread [Next in Thread]