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 .

{
    size_t              size;
    double              st, ft;
    const gsl_rng      *rng;
    gsl_siman_params_t  params;

    rng = gsl_rng_alloc(gsl_rng_rand);
    gsl_rng_set(rng, time(NULL));

    
    st = 4.0*siman_energy((void *)d);
    ft = st / 1e6;

    params.n_tries       = 100;
    params.iters_fixed_T =  25;
    params.step_size     =   0.001;
    params.k             =   1.0;
    params.t_initial     =   st;
    params.mu_t          =   1.02;
    params.t_min         =   ft;
    
    size = sizeof(d[0]);
    fprintf(stdout,"size: %d\n", size);
    
    gsl_siman_solve(rng, (void *)d, \
                    siman_energy, \
                    siman_step,   \
                    siman_metric, \
                    siman_print,  \
                    NULL, NULL, NULL,
                    size, params);
    
    return GSL_SUCCESS;
}

double siman_energy

(

void *

data

)

Returns:

The energy

The energy is defined as , which is given above.

{
    Parratt *d;
    int i;
    double tmp, tmp1;
    double sum1 = 0.0;
    double sum2 = 0.0;

    d = (Parratt *)data;
    
    
    d->g.plz = POLARIZED_UP;
    for (i=0; i<d->d1.num; i++)
    {
        tmp = calcReflectivity(d->d1.xpt[i], d->m.par, (Parratt *)data);
        tmp1 = (tmp - d->d1.ypt[i]) / (d->d1.sig[i]);
        sum1 += (tmp1 * tmp1);
    }
    sum2 = (double)(d->d1.num - d->m.num);
    
    if (d->s.pntp == POLARIZED)
    {
        d->g.plz = POLARIZED_DOWN;
        for (i=0; i<d->d2.num; i++)
        {
            tmp = calcReflectivity(d->d2.xpt[i], d->m.par, (Parratt *)data);
            tmp1 = (tmp - d->d2.ypt[i]) / (d->d2.sig[i]);
            sum1 += (tmp1 * tmp1);
        }
        sum2 += (double)d->d2.num;
    }
    
    tmp =  sum1 / sum2;
    return tmp;
}

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 .

{
    Parratt *d;
    int i;
    double s;
    
    d = (Parratt *)data;
    
    for (i=0; i<d->m.num; i++)
    {
        if (d->m.fix[i] == 0) 
        {
            s = gsl_rng_uniform(rng);
            d->m.par[i] = step_size * (2.0*s-1.0) + d->m.par[i];
        }
    }
    return;
}

double siman_metric	(	void *	data1,
		void *	data2
	)

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

{
    Parratt *d1;
    Parratt *d2;
    double sum;
    int i;

    d1 = (Parratt *)data1;
    d2 = (Parratt *)data2;
    sum = 0.0;
    for (i=0; i<d1->m.num; i++)
    {
        sum += fabs(d1->m.par[i] - d2->m.par[i]);
    }
    return sum;
}

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.

{
    Parratt *d;
    int i;
    
    d = (Parratt *)data;
    
    for (i=0; i<d->m.num; i++)
    {
        if (d->m.fix[i] == 0) 
        {
            printf("%8.4g", d->m.par[i]);
        }
    }
    return;
}

---