Commit 136826b2 authored by DOE CODE's avatar DOE CODE Committed by Tim Sowers
Browse files

Initial commit.

parents
In the these C code files, we implemented a version
of the Levenberg-Marquardt algorithm for fitting a
sum of exponentials convolved with an instrument
response function (IRF).
The file “convolvedexponential.c” contains functions
for calculating the model as well as the entry point
for calling fitting routines from external software.
Note that we use squared parameters to enforce
positivity. The main functions are:
fit_convolved_exponential entry point for fitting data
convolved_exp_f calculates function
convolved_exp_df calculates Jacobian
These routines must be called by external software since
there is no main() function. We implemented our fitting
routine as a DLL on a PC, and called the routine from
National Instrument LabVIEW, although there is no reason
the routines cannot be used with other compilers.
The function "fit_convolved_exponential" allows the user
to select which parameters in the function will be fitted.
Each element in the array "ip" is the index of the elements
of the parameter array "a" that will be fitted.
The parameters in the array "a" are:
a[0] = background constant
a[1] = amplitude of first lifetime
a[2] = first lifetime
a[3] = amplitude of second lifetime
a[4] = second lifetime
...
a[2*n-1] = amplitude of nth lifetime
a[2*n] = nth lifetime
The parameters for the function "fit_convolved_exponential"
are the following:
int n, // number of measurement bins
double deltaT, // time resolution for each bin
double *x, // Data Measured (size n array)
double *f, // Space for function (size n array)
double *irf, // Instrument Response (size n array)
int m_total, // Total number of parameters (including those not fitted)
double *a, // Parameter array (length m_total)
int m, // Number of fitted parameters
int *ip, // Array of length m with indices of fitted parameters in a
double *chisq, // Result for chisq
int *nIterations, // Number of iterations performed
int fitType, // 0 = MLE, 1 = Neyman weighting least squares, 2 = Equal Weighting least squares(sigma=1)
double deltaChisqLimit, // Stop criterion for change in chisq
int maxIterations // maximum number of iterations
The actual minimization routine dlevmar_mle_der is contained
in file “lm_mle.c”, and requires the two headers
“lm_mle_compiler.h” and “levmar_mle.h”. These files are
significantly modified forms of the code in the levmar package
(http://www.ics.forth.gr/~lourakis/levmar/). In the levmar
routines, there are many options for single or double
precision floating point arithmetic and constraints.
For our purposes, those were unnecessary complications.
We therefore stripped the code of all of those options,
and simply have a double precision implementation without
fitting constraints. The main differences in the code
for use with the MLE vs. least squares can be found in
the switch statements in the “lm_mle.c” file.
We use the GNU Scientific Library (GSL) routines for linear
algebra in the L-M routine, and for the FFT routine
for convolution calculations in accounting for
instrument response (http://www.gnu.org/software/gsl/).
For this reason, these routines must be linked with the
GNU scientific library before use.
We used Microsoft Visual Studio, and linked these routines
to the port of GSL 1.13 provided by David Geldreich at
http://david.geldreich.free.fr/dev.html
/* convolvedexponential.c -- model functions for multiple exponentials + background, convolved with an instrument response */
/////////////////////////////////////////////////////////////////////////////////
//
// Calculation of multiple exponentials convolved with an instrument response.
// using GNU Scientific Library. Also, calls Levenberg-Marquardt minimization routine
// for fitting event counting histograms. Options include using the Maximum Likelihood Estimator
// for Poisson deviates and least squares fitting.
// Copyright (C) 2010 Ted Laurence
// Lawrence Livermore National Laboratory
// Livermore, California
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_fft_halfcomplex.h>
#include <gsl/gsl_vector.h>
#include "levmar_mle.h"
struct data {
double deltaT; /*Spacing between bins*/
size_t n; /*Number of data elements*/
size_t nRadix2; /*Next higher power of 2 for n, giving size of irf and fWork arrays*/
double *irfWork;/*Instrument Response Function*/
double *fWork; /*Workspace for function calculation*/
size_t m_total; /*Number of parameters (number of exponentials = (m-1)/2 ) */
double *a; /*Parameter Array; fitted parameters x are subsituted into this array */
size_t m; /*Number of fitted parameters*/
int *ip; /*Mapping of fitted parameter array x to parameter array a*/
};
void complex_convolution_multiply(double *a, const double *b, int n)
{
int i,j;
double re,im,factor;
factor = 1.0/n; /* to normalize inverse FFT */
a[0]*=b[0]*factor;
a[n/2]*=b[n/2]*factor;
for (i=1;i<n/2;i++)
{
j=n-i;
re=a[i];
im=a[j];
a[i]=(re*b[i]-im*b[j])*factor;
a[j]=(im*b[i]+re*b[j])*factor;
}
}
void convolved_exp_f(double *p, double *f, int m, int n, void *data)
{
double deltaT=((struct data *)data)->deltaT;
size_t nRadix2 = ((struct data *)data)->nRadix2;
double *irfWork = ((struct data *) data)->irfWork;
double *fWork = ((struct data *) data)->fWork;
size_t m_total = ((struct data *)data)->m_total;
double *a = ((struct data *) data)->a;
int *ip = ((struct data *) data)->ip;
int i,j;
int numExponentials = (m_total-1)/2;
double tau, amp, sqrtTau,sqrtAmp;
double constantValue;
double totalT=deltaT*n;
double factor1;
double factor2;
if ( (n != ((struct data *)data)->n ) || (m != ((struct data *)data)->m) ) return;
for (i=0; i<m; i++)
a[ ip[i] ] = p[i];
constantValue=a[0]*a[0]/n;
for (i=0; i<n; i++) fWork[i] = constantValue;
for (i=n; i<nRadix2; i++) fWork[i] = 0.0;
for (j=0;j<numExponentials;j++)
{
sqrtAmp=a[2*j+1];
sqrtTau=a[2*j+2];
amp=sqrtAmp*sqrtAmp;
tau=sqrtTau*sqrtTau;
factor1=1-exp(-deltaT/tau);
factor2=1-exp(-totalT/tau);
for (i=0;i<n;i++)
fWork[i]+=amp*factor1/factor2*exp(-deltaT*i/tau);
}
gsl_fft_real_radix2_transform(fWork,1,nRadix2);
complex_convolution_multiply(fWork,irfWork,nRadix2);
gsl_fft_halfcomplex_radix2_backward(fWork,1,nRadix2);
// This loop uses 0 offset indexing
for (i=n;i<nRadix2;i++) fWork[i%n]+=fWork[i];
for (i=0;i<n;i++) f[i]=fWork[i];
return;
}
void convolved_exp_df(double *p, double *jac, int m, int n, void *data)
{
double deltaT=((struct data *)data)->deltaT;
size_t nRadix2 = ((struct data *)data)->nRadix2;
double *irfWork = ((struct data *) data)->irfWork;
double *fWork = ((struct data *) data)->fWork;
size_t m_total = ((struct data *)data)->m_total;
double *a = ((struct data *) data)->a;
int *ip = ((struct data *) data)->ip;
double *jac_row;
int i,k;
if ( (n != ((struct data *)data)->n ) || (m != ((struct data *)data)->m) ) return;
for (k=0; k<m; k++)
{
double tau, amp, sqrtTau,sqrtAmp;
double temp, totalT=deltaT*n;
double factor1;
double factor2;
if (ip[k]==0)
{
double d_constant=2*a[0]/n;
jac_row=jac;
for (i=0;i<n;i++,jac_row+=m)
jac_row[k]=d_constant;
continue;
}
for (i=0; i<nRadix2; i++) fWork[i] = 0.0;
if ((ip[k]-1)%2==0) /* Amplitude */
{
sqrtAmp=a[ip[k]];
sqrtTau=a[ip[k]+1];
amp=sqrtAmp*sqrtAmp;
tau=sqrtTau*sqrtTau;
factor1=1-exp(-deltaT/tau);
factor2=1-exp(-totalT/tau);
for (i=0;i<n;i++)
fWork[i]=2*sqrtAmp*factor1/factor2*exp(-deltaT*i/tau);
}
else /* lifetime */
{
sqrtAmp=a[ip[k]-1];
sqrtTau=a[ip[k]];
amp=sqrtAmp*sqrtAmp;
tau=sqrtTau*sqrtTau;
factor1=1-exp(-deltaT/tau);
factor2=1-exp(-totalT/tau);
for (i=0;i<n;i++)
{
temp=factor1/factor2*exp(-deltaT*i/tau);
fWork[i]= -2*amp*deltaT*exp(-(deltaT*i+deltaT)/tau)/(factor2*tau*sqrtTau)
+2*temp*amp*deltaT*i/(tau*sqrtTau)
+2*temp*amp*totalT*exp(-totalT/tau)/(factor2*tau*sqrtTau);
}
}
gsl_fft_real_radix2_transform(fWork,1,nRadix2);
complex_convolution_multiply(fWork,irfWork,nRadix2);
gsl_fft_halfcomplex_radix2_backward(fWork,1,nRadix2);
// This loop uses 0 offset indexing
for (i=n;i<nRadix2;i++) fWork[i%n]+=fWork[i];
jac_row=jac;
for (i=0;i<n;i++,jac_row+=m)
jac_row[k]=fWork[i];
}
return;
}
int fit_convolved_exponential(int n, // number of measurement bins
double deltaT, // time resolution for each bin
double *x, // Data Measured
double *f, // Space for function
double *irf, // Instrument Response
int m_total, // Total number of parameters (including those not fitted)
double *a, // Parameter array (length m_total)
int m, // Number of fitted parameters
int *ip, // Array of length m with indices of fitted parameters in a
double *chisq, // Result for chisq
int *nIterations, // Number of iterations performed
int fitType, // 0 = MLE, 1 = Neyman weighting, 2 = Equal Weighting (sigma=1)
double deltaChisqLimit, // Stop criterion for change in chisq
int maxIterations) // maximum number of iterations
{
struct data importedData;
double tempN,nRadix2;
double *irfWork, *fWork;
double opts[LM_OPTS_SZ],info[LM_INFO_SZ];
int i;
int status;
unsigned iter = 0;
double *p;
opts[0]=LM_INIT_MU;
opts[1]=LM_STOP_THRESH;
opts[2]=LM_STOP_THRESH;
opts[3]=LM_STOP_THRESH;
opts[4]=deltaChisqLimit;
importedData.n=n;
importedData.deltaT=deltaT;
importedData.m_total=m_total;
importedData.a=a;
importedData.m=m;
importedData.ip=ip;
for (i=0; i<m_total; i++) a[i]=sqrt(a[i]);
/* Determine power of 2 higher than n required for performing convolutions via FFT */
tempN = pow(2, 1.0+ceil( log((double)n)/log(2.0) ) );
nRadix2 = ( (int) tempN );
irfWork=malloc(nRadix2*sizeof(double));
if (irfWork==NULL) return GSL_ENOMEM;
fWork=malloc(nRadix2*sizeof(double));
if (fWork==NULL)
{
free(irfWork);
return GSL_ENOMEM;
}
p=malloc(m*sizeof(double));
if (p==NULL)
{
free(irfWork);
free(fWork);
return GSL_ENOMEM;
}
for (i=0; i<n; i++)
irfWork[i]=irf[i];
for (i=n; i<nRadix2; i++)
irfWork[i]=0.0;
gsl_fft_real_radix2_transform(irfWork,1,nRadix2);
importedData.nRadix2=nRadix2;
importedData.irfWork=irfWork;
importedData.fWork=fWork;
for (i=0; i<m; i++)
p[i]=a[ip[i]];
status = dlevmar_mle_der( &convolved_exp_f, &convolved_exp_df, p, x, m, n, maxIterations,opts,info,NULL,NULL,&importedData,fitType);
*chisq = info[1]/(n-m);
*nIterations=info[5];
convolved_exp_f(p, f, m, n, &importedData); /* Evaluate function for external use in array f */
for (i=0; i<m; i++)
a[ip[i]]=p[i];
free(p);
free(fWork);
free(irfWork);
for (i=0; i<m_total; i++) a[i]=a[i]*a[i];
return info[6];
}
This diff is collapsed.
/*
////////////////////////////////////////////////////////////////////////////////////
//
// Prototypes and definitions for the Levenberg - Marquardt minimization algorithm
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
////////////////////////////////////////////////////////////////////////////////////
*/
#include <float.h>
#ifndef _LEVMAR_H_
#define _LEVMAR_H_
#ifdef __cplusplus
extern "C" {
#endif
#define FABS(x) (((x)>=0.0)? (x) : -(x))
/* work arrays size for ?levmar_der and ?levmar_dif functions.
* should be multiplied by sizeof(double) or sizeof(float) to be converted to bytes
*/
#define LM_MLE_WORKSZ(npar, nmeas) (3*(nmeas) + 4*(npar) + (nmeas)*(npar) + 2*(npar)*(npar))
#define LM_OPTS_SZ 5
#define LM_INFO_SZ 10
#define LM_ERROR -1
#define LM_INIT_MU 1E-03
#define LM_STOP_THRESH 1E-20
#define LM_DIFF_DELTA 1E-20
#define LM_VERSION "2.5 (December 2009)"
#define LM_CHISQ_MLE 0
#define LM_CHISQ_NEYMAN 1
#define LM_CHISQ_EQUAL_WT 2
/* double precision LM, with Jacobian */
/* unconstrained minimization */
extern int dlevmar_mle_der(
void (*func)(double *p, double *hx, int m, int n, void *adata),
void (*jacf)(double *p, double *j, int m, int n, void *adata),
double *p, double *x, int m, int n, int itmax, double *opts,
double *info, double *work, double *covar, void *adata, int fitType);
extern double compute_chisq_measure(double *e, double *x, double *hx, int n, int fitType);
#ifdef __cplusplus
}
#endif
#endif /* _LEVMAR_H_ */
/////////////////////////////////////////////////////////////////////////////////
//
// Levenberg - Marquardt non-linear minimization algorithm
// Modified and simplified by Ted Laurence to use for MLE of Poisson-distributed data; Used only for
// double precision, without constraints
// Copyright (C) 2004 Manolis Lourakis (lourakis at ics forth gr)
// Institute of Computer Science, Foundation for Research & Technology - Hellas
// Heraklion, Crete, Greece.
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include "levmar_mle.h"
#include "lm_mle_compiler.h"
#define EPSILON 1E-12
#define ONE_THIRD 0.3333333334 /* 1.0/3.0 */
#define LM_REAL_MAX FLT_MAX
#define LM_REAL_MIN -FLT_MAX
double compute_chisq_measure(double *e, double *x, double *hx, int n, int fitType)
{
int i;
double chisq=0.0;
switch (fitType)
{
case LM_CHISQ_MLE:
for (i=0; i<n; i++)
{
if ( hx[i] > 0)
{
if (x[i]==0)
chisq += 2 * hx[i];
else
chisq += 2 * (hx[i]-x[i]-x[i]*log(hx[i]/x[i]));
e[i]=x[i]/hx[i]-1.0;
}
else
e[i]=0.0;
}
break;
case LM_CHISQ_NEYMAN:
for (i=0; i<n; i++)
{
if (x[i]==0)
{
chisq += (x[i]-hx[i])*(x[i]-hx[i]);
e[i]=x[i]-hx[i];
}
else
{
chisq += (x[i]-hx[i])*(x[i]-hx[i])/x[i];
e[i]=(x[i]-hx[i])/x[i];
}
}
break;
case LM_CHISQ_EQUAL_WT:
for (i=0; i<n; i++)
{
e[i]=x[i]-hx[i];
chisq += e[i]*e[i];
}
break;
}
return chisq;
}
/*
* This function seeks the parameter vector p that best describes the measurements vector x.
* More precisely, given a vector function func : R^m --> R^n with n>=m,
* it finds p s.t. func(p) ~= x, i.e. the squared second order (i.e. L2) norm of
* e=x-func(p) is minimized.
*
* This function requires an analytic Jacobian.
*
* Returns the number of iterations (>=0) if successful, LM_ERROR if failed
*
* For more details, see K. Madsen, H.B. Nielsen and O. Tingleff's lecture notes on
* non-linear least squares at http://www.imm.dtu.dk/pubdb/views/edoc_download.php/3215/pdf/imm3215.pdf
*/
int dlevmar_mle_der(
void (*func)(double *p, double *hx, int m, int n, void *adata), /* functional relation describing measurements. A p \in R^m yields a \hat{x} \in R^n */
void (*jacf)(double *p, double *j, int m, int n, void *adata), /* function to evaluate the Jacobian \part x / \part p */
double *p, /* I/O: initial parameter estimates. On output has the estimated solution */
double *x, /* I: measurement vector. NULL implies a zero vector */
int m, /* I: parameter vector dimension (i.e. #unknowns) */
int n, /* I: measurement vector dimension */
int itmax, /* I: maximum number of iterations */
double opts[5], /* I: minim. options [\mu, \epsilon1, \epsilon2, \epsilon3,\epsilon4]. Respectively the scale factor for initial \mu,
* stopping thresholds for ||J^T e||_inf, ||Dp||_2, chisq, and delta_chisq. Set to NULL for defaults to be used
*/
double info[LM_INFO_SZ],
/* O: information regarding the minimization. Set to NULL if don't care
* info[0]= ||e||_2 at initial p.
* info[1-4]=[ ||e||_2, ||J^T e||_inf, ||Dp||_2, mu/max[J^T J]_ii ], all computed at estimated p.
* info[5]= # iterations,
* info[6]=reason for terminating: 1 - stopped by small gradient J^T e
* 2 - stopped by small Dp
* 3 - stopped by itmax
* 4 - singular matrix. Restart from current p with increased mu
* 5 - no further error reduction is possible. Restart with increased mu
* 6 - stopped by small ||e||_2
* 7 - stopped by invalid (i.e. NaN or Inf) "func" values. This is a user error
* 8 - stopped by small change in chisq on successful iteration (dF<eps4)
* info[7]= # function evaluations
* info[8]= # Jacobian evaluations
* info[9]= # linear systems solved, i.e. # attempts for reducing error
*/
double *work, /* working memory at least LM_DER_WORKSZ() reals large, allocated if NULL */
double *covar, /* O: Covariance matrix corresponding to LS solution; mxm. Set to NULL if not needed. */
void *adata, /* pointer to possibly additional data, passed uninterpreted to func & jacf.
* Set to NULL if not needed
*/
int fitType) /* Which type of fit to perform
* LM_CHISQ_MLE is MLE for Poisson distribution
* LM_CHISQ_NEYMAN is least squares where sigma^2 is set to max(x,1)
* LM_CHISQ_EQUAL_WT is least squares where sigma^2 is set to 1 */
{
register int i, j, k, l;
int worksz, freework=0, issolved;
/* temp work arrays */
double *e, /* nx1 */
*e_test, /* nx1 */
*hx, /* \hat{x}_i, nx1 */
*jacTe, /* J^T e_i mx1 */
*jac, /* nxm */
*jacTjac, /* mxm */
*LUdecomp, /* mxm */
*Dp, /* mx1 */
*diag_jacTjac, /* diagonal of J^T J, mx1 */
*pDp; /* p + Dp, mx1 */
register double mu, /* damping constant */
tmp; /* mainly used in matrix & vector multiplications */
double p_chisq, jacTe_inf, pDp_chisq; /* ||e(p)||_2, ||J^T e||_inf, ||e(p+Dp)||_2 */
double p_L2, Dp_L2=LM_REAL_MAX, dF, dL;