Parratt: Non-linear least squares

These robust routines are implemented using the GNU Scientific Library. By using the fit_method struct, the user can choose which of the two Levenberg-Marquardt algorithms available in the GSL to use.


Functions
int	lstsq_main (Parratt *d)
	The main function for Non-linear least squares.
int	lstsq_fdf (const gsl_vector params, void data, gsl_vector f, gsl_matrix J)
	Required by the GSL library.
int	lstsq_f (const gsl_vector params, void data, gsl_vector *f)
	Calculates the residuals.
int	lstsq_df (const gsl_vector params, void data, gsl_matrix *J)
	Calculates the Jacobian matrix.
double	lstsq_derivatives1 (double p, void *data)
	Numerically calculate the residual gradients.
double	lstsq_derivatives2 (double p, void *data)
	Numerically calculate the residual gradients.
void	lstsq_status (int i, gsl_multifit_fdfsolver s, gsl_matrix c, Parratt *d)
	Updates the fit to `stdout`.

Function Documentation

int lstsq_main ( Parratt * d )

The method minimizes the sum of residuals $r_i = r(q_i,R_i,\sigma_i;\mathbf{a})$ , with respect to a set of parameters $\mathbf{a}=[a_1, a_2, \ldots]$ , which define the layered structure. The function to be minimized is given by,

$\begin{eqnarray*} \Phi &=& \frac{1}{2} \left( \sum_{i}^{\mathrm{up}} r_i^2(q_i,R_i,\sigma_i;\mathbf{a}) \right) \\ r_i &=& \frac{R(q_i;\mathbf{a}) - R_i}{\sigma_i} \end{eqnarray*}$

where the sum is over the data points $(q_i,R_i,\sigma_i)$ .

The iteration can be stopped by an error, or by successful completion. There are two stopping criterion, and the iteration stops on the success of either one. The first is when the gradient falls below a tolerance level. The test can be summarized as,

$\left|\frac{\partial\Phi}{\partial a_i}\right| \leq \epsilon.$

The second stop test compares the residuals, $r$ , of the last step with the residuals of the current step according to,

$|r^{i-1}| \leq \epsilon_{\mathrm{absolute}} + \epsilon_{\mathrm{relative}}|r^i|$

Perhaps a later version will allow user configuration of the tolerances.

As in the other fitting methods, the GSL library takes a (void *) struct as the fitting parameters, so this is set to be the Parratt struct. In this way the data is totally accessible. The parameters are overwritten with this technique, however, since the Parratt struct is used as workspace.

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     unsigned int  iter;
00025     int           status, status2, ndat;
00026     gsl_vector                         *guess;
00027     gsl_vector                         *gradt;
00028     gsl_matrix                         *covar;
00029     gsl_multifit_function_fdf           my_func;
00030     gsl_multifit_fdfsolver             *solv;
00031     
00032     /* Initial guess */
00033     guess = gsl_vector_alloc(d->m.num);
00034     for (iter=0; iter<d->m.num; iter++) 
00035     {
00036         gsl_vector_set(guess, iter, d->m.par[iter]);
00037     }
00038 
00039     /* Need for stop test #2 */
00040     gradt = gsl_vector_alloc(d->m.num);
00041     
00042     /* Co-variant matrix      */
00043     covar = gsl_matrix_alloc(d->m.num,d->m.num);
00044 
00045     my_func.f      = &lstsq_f;
00046     my_func.df     = &lstsq_df;
00047     my_func.fdf    = &lstsq_fdf;
00048     /* Number of data points */
00049     if (d->s.pntp == POLARIZED) ndat = d->d1.num + d->d2.num;
00050     else                        ndat = d->d1.num;
00051     my_func.n      = ndat;
00052     /* Number of parameters */
00053     my_func.p      = d->m.num;
00054     my_func.params = (void *)d;
00055     
00056     solv = gsl_multifit_fdfsolver_alloc(fit_method.lstsq_type, \
00057                                         ndat, d->m.num);
00058 
00059     /*==================================================*/
00060     iter = 0;
00061     gsl_multifit_fdfsolver_set(solv, &my_func, guess);
00062     gsl_multifit_covar(solv->J, 0.0, covar);
00063     lstsq_status(iter, solv, covar, d);
00064 
00065     do
00066     {
00067         iter++;
00068         status = gsl_multifit_fdfsolver_iterate(solv);      
00069         gsl_multifit_covar(solv->J, 0.0, covar);
00070         lstsq_status(iter, solv, covar, d);
00071         
00072         if (status) 
00073         {
00074             fprintf(stderr,"STATUS: %s\n", gsl_strerror(status));
00075             break;
00076         }
00077         
00078         /* Stop test #1 */
00079         status = gsl_multifit_test_delta(solv->dx, solv->x, 1e-7, 1e-7);
00080 
00081         /* Stop test #2 */      
00082         gsl_multifit_gradient(solv->J, solv->f, gradt);
00083         status2 = gsl_multifit_test_gradient(gradt, 1e-7);
00084         
00085         if ((status==GSL_SUCCESS) || (status2==GSL_SUCCESS))
00086         {
00087             fprintf(stderr,"STATUS1: %s\n", gsl_strerror(status));
00088             fprintf(stderr,"STATUS2: %s\n", gsl_strerror(status2));
00089         }
00090 
00091     } while ((status == GSL_CONTINUE)  && \
00092              (status2 == GSL_CONTINUE) && \
00093              (iter < MAX_ITER));
00094 
00095     /*==================================================*/
00096     gsl_multifit_covar(solv->J, 0.0, covar);
00097     lstsq_status(-1, solv, covar, d);
00098     
00099     /*==================================================*/
00100 
00101     gsl_multifit_fdfsolver_free(solv);
00102     gsl_matrix_free(covar);
00103     gsl_vector_free(guess);
00104     gsl_vector_free(gradt);
00105 
00106     return GSL_SUCCESS;
00107 }

int lstsq_fdf	(	const gsl_vector *	params,
		void *	data,
		gsl_vector *	f,
		gsl_matrix *	J
	)

This function is required by the GSL, and merely calls the functions that fills the vector f with the residuals and the matrix J with the derivatives.

00110 {
00111     lstsq_f(params, data, f);
00112     lstsq_df(params, data, J);
00113     return GSL_SUCCESS;
00114 }

int lstsq_f	(	const gsl_vector *	params,
		void *	data,
		gsl_vector *	f
	)

The array f is the residuals $r_i$ .

00117 {
00118     int     i, j;
00119     double  tmp;
00120     Parratt    *d;
00121     
00122     d = (Parratt *)data;
00123     
00124     for (i=0; i<d->m.num; i++) 
00125     {
00126         d->m.par[i] = gsl_vector_get(params,i);
00127     }
00128     
00129     d->g.plz = POLARIZED_UP;
00130     for(i=0; i<d->d1.num; i++)
00131     {
00132         tmp = (calcReflectivity(d->d1.xpt[i], d->m.par, (Parratt *)data) \
00133                               - d->d1.ypt[i]) \
00134                               / d->d1.sig[i];    
00135         gsl_vector_set(f, i, tmp);
00136     }
00137 
00138     if (d->s.pntp == POLARIZED)
00139     {
00140         d->g.plz = POLARIZED_DOWN;
00141         for(i=0, j=d->d1.num; i<d->d2.num; i++, j++)
00142         {
00143             tmp = (calcReflectivity(d->d2.xpt[i], d->m.par, (Parratt *)data) \
00144                                   - d->d2.ypt[i]) \
00145                                   / d->d2.sig[i];    
00146             gsl_vector_set(f, j, tmp);
00147         }
00148     }
00149     
00150     return GSL_SUCCESS;
00151 }

int lstsq_df	(	const gsl_vector *	params,
		void *	data,
		gsl_matrix *	J
	)

The gsl_matrix J is the Jacobian matrix, whose elements are

$J_{ij} = \partial r_i/\partial a_j.$

This is calculated using numerical derivatives. If the parameter is to be held fixed, then the element is set to zero.

00154 {
00155     int i, j, k;
00156     gsl_function F;
00157     double res, err;
00158     Parratt   *d;
00159     
00160     d = (Parratt *)data;
00161     F.params   = data;
00162 
00163     for (i=0; i<d->m.num; i++) 
00164     {
00165         d->m.par[i] = gsl_vector_get(params,i);
00166     }
00167 
00168     
00169     F.function = &lstsq_derivatives1;
00170     d->g.plz = POLARIZED_UP;
00171     for(i=0; i<d->d1.num; i++)
00172     {
00173         d->d1.idx = i;
00174         
00175         for (j=0; j<d->m.num; j++) 
00176         {
00177             d->m.idx = j;
00178             if (d->m.fix[j]) res=0.0;
00179             else gsl_deriv_central(&F, d->m.par[j], 1e-4, &res, &err);
00180             gsl_matrix_set(J, i, j, res);
00181         }
00182     }
00183         
00184     if (d->s.pntp == POLARIZED)
00185     {
00186         F.function = &lstsq_derivatives2;
00187         d->g.plz = POLARIZED_DOWN;
00188         for(i=0, k=d->d1.num; i<d->d2.num; i++, k++)
00189         {
00190             d->d2.idx = i;
00191         
00192             for (j=0; j<d->m.num; j++) 
00193             {
00194                 d->m.idx = j;
00195                 if (d->m.fix[j]) res=0.0;
00196                 else gsl_deriv_central(&F, d->m.par[j], 1e-4, &res, &err);
00197                 gsl_matrix_set(J, k, j, res);
00198             }
00199         }
00200     }
00201     
00202     return GSL_SUCCESS;
00203 }

double lstsq_derivatives1	(	double	p,
		void *	data
	)

This function is needed by gsl_deriv_central to return the value the residual $r_i$ at slightly changing values of the parameter $a_j$ .

In order for this trick to work, the parameter that the derivative is with respect to is flagged by d->m.idx and the data point the derivative is taken at is d->d.idx . The value of the other parameters are unchanged in this derivative, since they were put into d->m.par , back in the function mdmin_df() , they are available to fill the vector a.

00206 {
00207     int     i;
00208     double  tmp, a[MAX_PARAM];
00209     Parratt *d;
00210     
00211     d = (Parratt *)data;    
00212     
00213     for (i=0; i<d->m.num; i++) 
00214     {
00215         if (d->m.idx==i) a[i] = p;
00216         else a[i] = d->m.par[i];
00217     }
00218             
00219     tmp = (calcReflectivity(d->d1.xpt[d->d1.idx], a, (Parratt *)data) \
00220                           - d->d1.ypt[d->d1.idx]) \
00221                           / d->d1.sig[d->d1.idx];    
00222     return tmp;
00223 }

void lstsq_status	(	int	i,
		gsl_multifit_fdfsolver *	s,
		gsl_matrix *	c,
		Parratt *	d
	)

Parameters:

	i	Current iteration
	s	Current state of the solution

This is a report of the current state of the solution.

00246 {
00247     int      j;
00248     double   chisq;
00249     
00250     if (i >= 0)
00251     {
00252         fprintf(stderr,"ITER:%3u ", i);
00253         fprintf(stderr,"===================\n");
00254     }
00255     else
00256     {
00257         fprintf(stderr,"RESULT ");
00258         fprintf(stderr,"===================\n");
00259     }
00260     
00261     /*
00262      * Chisq
00263      */
00264     if (d->s.pntp == POLARIZED)
00265     {
00266         chisq = gsl_blas_dnrm2(s->f) / \
00267                (double)((d->d1.num)+(d->d2.num)-(d->m.num));
00268     }
00269     else
00270     {
00271         chisq = gsl_blas_dnrm2(s->f) / \
00272                (double)((d->d1.num)-(d->m.num));
00273     }
00274     
00275     /*
00276      *  Report results
00277      */
00278     fprintf(stderr,"ChiSq: %12.6e\n", chisq);
00279     for (j=0; j<d->m.num; j++)
00280     {
00281         if (d->m.fix[j] == 0) 
00282         {
00283             fprintf(stderr,"%3d: %12.6f +/- %8.6f\n",
00284                     j,
00285                     gsl_vector_get(s->x, j),
00286                     sqrt(gsl_matrix_get(c, j, j)));
00287             d->m.par[j] = gsl_vector_get(s->x, j);
00288         }
00289     }
00290     fprintf(stderr,"\n");
00291     
00292     
00293     return;
00294 }

Non-linear least squares

Detailed Description

Functions

Function Documentation