Re: [Help-gsl] Simulated Annealing

[John Gehman]

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.
> >
> > Regards,
> > john
> >
> > > 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
> > >  
> > http://www.gnu.org/software/gsl/manual/html_node/Trivial- 
> > 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
> = {N_TRIES, ITERS_FIXED_T, STEP_SIZE,
>    K, T_INITIAL, MU_T, T_MIN};
>     
> /* 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));
> }
>     
> int
> main(int argc, char *argv[])
> {
>   const gsl_rng_type * T;
>   gsl_rng *  r;
>   double x_initial = 15.5;
>   gsl_rng_env_setup();
>   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
>           M1,        
>           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,
john

int SimanObject::Unknown_siman( Paar &P ) {

        const gsl_rng_type * Tr;
        gsl_rng * r;

        gsl_rng_env_setup();

        Tr = gsl_rng_taus2;
        r = gsl_rng_alloc(Tr);

        ofstream saout; 
        char saname[179];
        strcpy(saname,P.path);
        strcat(saname,P.outbase);
        strcat(saname,".sa"); 
        saout.open(saname, ios::out);

        sprintf(P.str,"#%4s %7s %12s %10s %4s %4s %4s\n",
                "iter","evals","temp","chi2","#lss","#acc","$rej");
        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",
                        iterT,ntry,boltzmann,maxstep,T,mintemp,attentemp);
                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) {
                                new_x.copy_min(x);
                                // take a step
                                for (int lm=0; lm<P.nc; 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;
                                        ++n_rejects;
                                        continue;
                                }
                                if(new_E <= best_E){
                                        best_x.copy_min(new_x);
                        best_E=new_E;
                                        magicTemp = T;
                                }
                                // keep track of Ef() evaluations
                                ++n_evals;
                                //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
                                        x.copy_min(new_x);
                                        E = new_E;
                                        ++n_eless;
                                        
                                } else if
                                        (gsl_rng_uniform(r) < exp 
(-(new_E-E)/(boltzmann*T)) ) {
                                        // also take a step */
                                        x.copy_min(new_x);
                                        E = new_E;
                                        ++n_accepts;
                                } else {
                                        ++n_rejects;
                                }
                }
        
                        //print report
                        char format[79];
                        strcpy(format,"%5d %7d %12g %10.8f %4.2f %4.2f %4.2f ");
                        sprintf(P.str,format,n_iter,n_evals,T,x.chi2,
                                100.0*double(n_eless)/double(iterT),
                                100.0*double(n_accepts)/double(iterT),
                                100.0*double(n_rejects)/double(iterT) );
                        saout << P.str;
                        for (int i=0; i<P.nc; 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 
%14.6g",T,E,best_E);
                        }
        
                // apply the cooling schedule to the temperature
                        T /= attentemp;
                        ++n_iter;
                        if (T < mintemp) done = 1;
                }
        }
        //at the end, copy the result onto the initial point
        P.copy_min(best_x);
        P.calc_state_dataonly_nosim();

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

        saout.close();
        return 0;
}

---------------------------------------------------------

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