Simulated Annealing

Detailed Description

Simulated annealing is a generic probabilistic algorithm for global optimization of a function in a large search space. In this case, we are optimizing the Chi squared measure of the fit for a given set of parameters. The search space is the space of possible parameter values. The simulated annealing routine here is implemented using the GNU Scientific Library. Please refer to the Wikipedia article on simulated annealing for more information. (http://en.wikipedia.org/wiki/Simulated_annealing)

When using this technique, it can be useful to direct the printout to be saved in a file. It will contain the system energy and parameter values over the course of the simulation, which can be analyzed later.

Functions

int siman_main (Parratt *d)
The main function for simulated annealing.
double siman_energy (void *data)
Calculates the energy of the configuration.
void siman_step (const gsl_rng *rng, void *data, double step_size)
Steps the parameters randomly.
double siman_metric (void *data1, void *data2)
Calculates the metric for the configuration.
void siman_print (void *data)
Prints the current configuration.

Function Documentation

 int siman_main ( Parratt * d )

There are several parameters here that really should be configurable, such as the number of trial steps, and most importantly, the step size. Perhaps a later version will allow user configuration of these.

The initial temperature is set to be four times initial energy, which is just evaluated with the initial parameters, where is given by

The final temperature is times smaller. This works well when the Boltzmann constant is set to 1.

As in the other fitting methods, the system configuration, (void *)d, is set to be the Parratt struct, so that the data is totally accessible. The parameters are overwritten with this technique, however, since the Parratt struct is used as workspace.

The difficulty is choosing the best step size, since there is only one that must be shared among all the variables. When the SLD and background is given in the scale of , then all variables hopefully are of the similar order of magnitude, and can be varied by a similar step size.

The code requires one external function, calcReflectivity(), which takes the form

double calcReflectivity(double q, double *a, Parratt *p)

Parameters:
 q The value to evaluate . a The set of all adjustable parameters in one array. p The Parratt struct containing all the other information necessary for calculation of .
Returns:
The reflectivity .
00023 {
00024     size_t              size;
00025     double              st, ft;
00026     const gsl_rng      *rng;
00027     gsl_siman_params_t  params;
00028
00029     rng = gsl_rng_alloc(gsl_rng_rand);
00030     gsl_rng_set(rng, time(NULL));
00031
00032
00033     st = 4.0*siman_energy((void *)d);
00034     ft = st / 1e6;
00035
00036     params.n_tries       = 100;
00037     params.iters_fixed_T =  25;
00038     params.step_size     =   0.001;
00039     params.k             =   1.0;
00040     params.t_initial     =   st;
00041     params.mu_t          =   1.02;
00042     params.t_min         =   ft;
00043
00044     size = sizeof(d[0]);
00045     fprintf(stdout,"size: %d\n", size);
00046
00047     gsl_siman_solve(rng, (void *)d, \
00048                     siman_energy, \
00049                     siman_step,   \
00050                     siman_metric, \
00051                     siman_print,  \
00052                     NULL, NULL, NULL,
00053                     size, params);
00054
00055     return GSL_SUCCESS;
00056 }


 double siman_energy ( void * data )

Returns:
The energy
The energy is defined as , which is given above.
00059 {
00060     Parratt *d;
00061     int i;
00062     double tmp, tmp1;
00063     double sum1 = 0.0;
00064     double sum2 = 0.0;
00065
00066     d = (Parratt *)data;
00067
00068
00069     d->g.plz = POLARIZED_UP;
00070     for (i=0; i<d->d1.num; i++)
00071     {
00072         tmp = calcReflectivity(d->d1.xpt[i], d->m.par, (Parratt *)data);
00073         tmp1 = (tmp - d->d1.ypt[i]) / (d->d1.sig[i]);
00074         sum1 += (tmp1 * tmp1);
00075     }
00076     sum2 = (double)(d->d1.num - d->m.num);
00077
00078     if (d->s.pntp == POLARIZED)
00079     {
00080         d->g.plz = POLARIZED_DOWN;
00081         for (i=0; i<d->d2.num; i++)
00082         {
00083             tmp = calcReflectivity(d->d2.xpt[i], d->m.par, (Parratt *)data);
00084             tmp1 = (tmp - d->d2.ypt[i]) / (d->d2.sig[i]);
00085             sum1 += (tmp1 * tmp1);
00086         }
00087         sum2 += (double)d->d2.num;
00088     }
00089
00090     tmp =  sum1 / sum2;
00091     return tmp;
00092 }


 void siman_step ( const gsl_rng * rng, void * data, double step_size )

The parameters are changed by a maximum of step_size . The variable s is randomly chosen from the uniform interval .

00095 {
00096     Parratt *d;
00097     int i;
00098     double s;
00099
00100     d = (Parratt *)data;
00101
00102     for (i=0; i<d->m.num; i++)
00103     {
00104         if (d->m.fix[i] == 0)
00105         {
00106             s = gsl_rng_uniform(rng);
00107             d->m.par[i] = step_size * (2.0*s-1.0) + d->m.par[i];
00108         }
00109     }
00110     return;
00111 }


 double siman_metric ( void * data1, void * data2 )

This returns the metric , defined as the distance between parameters between this step () and the last step (),

00114 {
00115     Parratt *d1;
00116     Parratt *d2;
00117     double sum;
00118     int i;
00119
00120     d1 = (Parratt *)data1;
00121     d2 = (Parratt *)data2;
00122     sum = 0.0;
00123     for (i=0; i<d1->m.num; i++)
00124     {
00125         sum += fabs(d1->m.par[i] - d2->m.par[i]);
00126     }
00127     return sum;
00128 }


 void siman_print ( void * data )

The GSL requests this function to print the output to the stdout. In this case we only print those variables not fixed.

00131 {
00132     Parratt *d;
00133     int i;
00134
00135     d = (Parratt *)data;
00136
00137     for (i=0; i<d->m.num; i++)
00138     {
00139         if (d->m.fix[i] == 0)
00140         {
00141             printf("%8.4g", d->m.par[i]);
00142         }
00143     }
00144     return;
00145 }


Generated on Wed Mar 14 13:24:56 2007 for Parratt by  1.4.7