expfit.h

用gsl实现Levenberg-Marquardt算法

/* expfit.h -- model functions for exponential + background */
#ifndef _EXPFIT
#define _EXPFIT

struct data {
	size_t n;
	double * y;
	double * sigma;
};


/*设置需要优化参数的函数*/
/* 这里函数设置为 Yi = A * exp(-lambda * i) + b */
/*需要优化A、lambda、b这三个参数*/
int
expb_f (const gsl_vector * x, void *data, 
		gsl_vector * f)
{
	size_t n = ((struct data *)data)->n;
	double *y = ((struct data *)data)->y;
	double *sigma = ((struct data *) data)->sigma;
	
	double A = gsl_vector_get (x, 0);
	double lambda = gsl_vector_get (x, 1);
	double b = gsl_vector_get (x, 2);
	
	size_t i;
	
	for (i = 0; i < n; i++)
	{
		/* Model Yi = A * exp(-lambda * i) + b */
		double t = i;
		double Yi = A * exp (-lambda * t) + b;
		gsl_vector_set (f, i, (Yi - y[i])/sigma[i]);
		/*This function sets the value of the i-th element of a vector v to x. 
		If i lies outside the allowed range of 0 to n-1 then the error handler is invoked. */
	}
	
	return GSL_SUCCESS;
}


/*计算函数的导数*/
int
expb_df (const gsl_vector * x, void *data, 
		 gsl_matrix * J)
{
	size_t n = ((struct data *)data)->n;
	double *sigma = ((struct data *) data)->sigma;
	
	double A = gsl_vector_get (x, 0);
	double lambda = gsl_vector_get (x, 1);
	
	size_t i;
	
	for (i = 0; i < n; i++)
	{
		/* Jacobian matrix J(i,j) = dfi / dxj, */
		/* where fi = (Yi - yi)/sigma[i],      */
		/*       Yi = A * exp(-lambda * i) + b  */
		/* and the xj are the parameters (A,lambda,b) */
		double t = i;
		double s = sigma[i];
		double e = exp(-lambda * t);
		gsl_matrix_set (J, i, 0, e/s); 
		gsl_matrix_set (J, i, 1, -t * A * e/s);
		gsl_matrix_set (J, i, 2, 1/s);
	}
	return GSL_SUCCESS;
}

int
expb_fdf (const gsl_vector * x, void *data,
		  gsl_vector * f, gsl_matrix * J)
{
	expb_f (x, data, f);
	expb_df (x, data, J);
	
	return GSL_SUCCESS;
}
#endif