[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.
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)
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- Re: [Help-gsl] Simulated Annealing,
John Gehman <=