Gromacs  2022.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Functions | Variables
expfit.cpp File Reference
#include "gmxpre.h"
#include "expfit.h"
#include <cstring>
#include <algorithm>
#include "gromacs/correlationfunctions/integrate.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/math/vec.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
#include "gmx_lmcurve.h"
+ Include dependency graph for expfit.cpp:

Description

Implements routine for fitting a data set to a curve.

Author
David van der Spoel david.nosp@m..van.nosp@m.dersp.nosp@m.oel@.nosp@m.icm.u.nosp@m.u.se

Functions

int effnNparams (int effn)
 Returns number of function parameters associated with a fitting function. More...
 
const char * effnDescription (int effn)
 Returns description corresponding to the enum above, or NULL if out of range. More...
 
int sffn2effn (const char **sffn)
 Returns corresponding to the selected enum option in sffn. More...
 
static double myexp (double x, double A, double tau)
 Compute exponential function A exp(-x/tau)
 
static double lmc_erffit (double x, const double *a)
 Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2)
 
static double safe_exp (double x)
 Exponent function that prevents overflow.
 
static double safe_expm1 (double x)
 Exponent minus 1 function that prevents overflow.
 
static double lmc_exp_one_parm (double x, const double *a)
 Compute y = exp(-x/|a0|)
 
static double lmc_exp_two_parm (double x, const double *a)
 Compute y = a1 exp(-x/|a0|)
 
static double lmc_exp_exp (double x, const double *a)
 Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|)
 
static double lmc_exp_5_parm (double x, const double *a)
 Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4.
 
static double lmc_exp_7_parm (double x, const double *a)
 Compute 7 parameter exponential function value. More...
 
static double lmc_exp_9_parm (double x, const double *a)
 Compute 9 parameter exponential function value. More...
 
static double lmc_pres_6_parm (double x, const double *a)
 Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|)
 
static double lmc_vac_2_parm (double x, const double *a)
 Compute vac function.
 
static double lmc_errest_3_parm (double x, const double *a)
 Compute error estimate.
 
double fit_function (const int eFitFn, const double parm[], const double x)
 Returns the value of fit function eFitFn at x. More...
 
static void initiate_fit_params (int eFitFn, double params[])
 Ensure the fitting parameters are well-behaved. More...
 
static void extract_fit_params (int eFitFn, double params[])
 Process the fitting parameters to get output parameters. More...
 
static void print_chi2_params (FILE *fp, const int eFitFn, const double fitparms[], const char *label, const int nfitpnts, const double x[], const double y[])
 Print chi-squared value and the parameters.
 
real do_lmfit (int ndata, const real c1[], real sig[], real dt, const real *x0, real begintimefit, real endtimefit, const gmx_output_env_t *oenv, gmx_bool bVerbose, int eFitFn, double fitparms[], int fix, const char *fn_fitted)
 Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential or to a transverse current autocorrelation function. More...
 
real fit_acf (int ncorr, int fitfn, const gmx_output_env_t *oenv, gmx_bool bVerbose, real tbeginfit, real tendfit, real dt, real c1[], real *fit)
 Fit an autocorrelation function to a pre-defined functional form. More...
 

Variables

static const int nfp_ffn [effnNR] = { 0, 1, 2, 3, 5, 7, 9, 2, 4, 3, 6 }
 Number of parameters for each fitting function.
 
const char * s_ffn [effnNR+2]
 Short name of each function type. This is exported for now in order to use when calling parse_common_args. More...
 
static const char * longs_ffn [effnNR]
 Long description for each fitting function type. More...
 
t_lmcurve lmcurves [effnNR+1]
 array of fitting functions corresponding to the pre-defined types More...
 

Function Documentation

real do_lmfit ( int  ndata,
const real  c1[],
real  sig[],
real  dt,
const real x0,
real  begintimefit,
real  endtimefit,
const gmx_output_env_t *  oenv,
gmx_bool  bVerbose,
int  eFitFn,
double  fitparms[],
int  fix,
const char *  fn_fitted 
)

Use Levenberg-Marquardt method to fit to a nfitparm parameter exponential or to a transverse current autocorrelation function.

If x == NULL, the timestep dt will be used to create a time axis.

Parameters
[in]ndataNumber of data points
[in]c1The data points
[in]sigThe standard deviation in the points (can be NULL)
[in]dtThe time step
[in]x0The X-axis (may be NULL, see above)
[in]begintimefitStarting time for fitting
[in]endtimefitEnding time for fitting
[in]oenvOutput formatting information
[in]bVerboseShould the routine write to console?
[in]eFitFnFitting function (0 .. effnNR)
[in,out]fitparms[]Fitting parameters, see printed manual for a detailed description. Note that in all implemented cases the parameters corresponding to time constants will be generated with increasing values. Such input parameters should therefore be provided in increasing order. If this is not the case or if subsequent time parameters differ by less than a factor of 2, they will be modified to ensure tau_i+1 >= 2 tau_i.
[in]fixConstrains fit parameter i at it's starting value, when the i'th bit of fix is set. This works only when the N last parameters are fixed but not when a parameter somewhere in the middle needs to be fixed.
[in]fn_fittedIf not NULL file to print the data and fitted curve to
Returns
integral.
const char* effnDescription ( int  effn)

Returns description corresponding to the enum above, or NULL if out of range.

Parameters
[in]effnIndex
Returns
Description or NULL
int effnNparams ( int  effn)

Returns number of function parameters associated with a fitting function.

Parameters
[in]effnIndex
Returns
number or -1 if index out of range
static void extract_fit_params ( int  eFitFn,
double  params[] 
)
static

Process the fitting parameters to get output parameters.

See comment at the previous function.

real fit_acf ( int  ncorr,
int  fitfn,
const gmx_output_env_t *  oenv,
gmx_bool  bVerbose,
real  tbeginfit,
real  tendfit,
real  dt,
real  c1[],
real fit 
)

Fit an autocorrelation function to a pre-defined functional form.

Todo:
check parameters
Parameters
[in]ncorr
[in]fitfnFitting function (0 .. effnNR)
[in]oenvOutput formatting information
[in]bVerboseShould the routine write to console?
[in]tbeginfitStarting time for fitting
[in]tendfitEnding time for fitting
[in]dtThe time step
[in]c1The data points
[in,out]fitThe fitting parameters
Returns
the integral over the autocorrelation function?
double fit_function ( int  eFitFn,
const double  parm[],
double  x 
)

Returns the value of fit function eFitFn at x.

Parameters
[in]eFitFnthe index to the fitting function (0 .. effnNR)
[in]parmArray of parameters, the length of which depends on eFitFn
[in]xThe value of x
Returns
the value of the fit
static void initiate_fit_params ( int  eFitFn,
double  params[] 
)
static

Ensure the fitting parameters are well-behaved.

In order to make sure that fitting behaves according to the constraint that time constants are positive and increasing we enforce this here by setting all time constants to their absolute value and by adding e.g. |a_0| to |a_2|. This is done in conjunction with the fitting functions themselves. When there are multiple time constants we make sure that the successive time constants are at least double the previous ones and during fitting we enforce the they remain larger. It may very well help the convergence of the fitting routine.

static double lmc_exp_7_parm ( double  x,
const double *  a 
)
static

Compute 7 parameter exponential function value.

Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 exp(-x/(|a1|+|a3|+|a5|)) + a6

static double lmc_exp_9_parm ( double  x,
const double *  a 
)
static

Compute 9 parameter exponential function value.

Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8

int sffn2effn ( const char **  sffn)

Returns corresponding to the selected enum option in sffn.

Parameters
[in]sffnTwo dimensional string array coming from parse_common_args
Returns
the ffn enum

Variable Documentation

t_lmcurve lmcurves[effnNR+1]
Initial value:
static double lmc_vac_2_parm(double x, const double *a)
Compute vac function.
Definition: expfit.cpp:303
static double lmc_exp_one_parm(double x, const double *a)
Compute y = exp(-x/|a0|)
Definition: expfit.cpp:206
static double lmc_erffit(double x, const double *a)
Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2)
Definition: expfit.cpp:142
static double lmc_exp_5_parm(double x, const double *a)
Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4.
Definition: expfit.cpp:228
static double lmc_errest_3_parm(double x, const double *a)
Compute error estimate.
Definition: expfit.cpp:351
static double lmc_exp_two_parm(double x, const double *a)
Compute y = a1 exp(-x/|a0|)
Definition: expfit.cpp:212
static double lmc_exp_9_parm(double x, const double *a)
Compute 9 parameter exponential function value.
Definition: expfit.cpp:261
static double lmc_pres_6_parm(double x, const double *a)
Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|)
Definition: expfit.cpp:279
static double lmc_exp_exp(double x, const double *a)
Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|)
Definition: expfit.cpp:218
static double lmc_exp_7_parm(double x, const double *a)
Compute 7 parameter exponential function value.
Definition: expfit.cpp:242

array of fitting functions corresponding to the pre-defined types

const char* longs_ffn[effnNR]
static
Initial value:
= { "no fit",
"y = exp(-x/|a0|)",
"y = a1 exp(-x/|a0|)",
"y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0",
"y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4, a3 >= a1",
"y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1",
"y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1",
"y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)",
"y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)",
"y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)",
"y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)" }

Long description for each fitting function type.

const char* s_ffn[effnNR+2]
Initial value:
= { nullptr, "none", "exp", "aexp", "exp_exp", "exp5", "exp7",
"exp9", nullptr, nullptr, nullptr, nullptr, nullptr }

Short name of each function type. This is exported for now in order to use when calling parse_common_args.