#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\|) + a0exp(-(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]	ndata	Number of data points
[in]	c1	The data points
[in]	sig	The standard deviation in the points (can be NULL)
[in]	dt	The time step
[in]	x0	The X-axis (may be NULL, see above)
[in]	begintimefit	Starting time for fitting
[in]	endtimefit	Ending time for fitting
[in]	oenv	Output formatting information
[in]	bVerbose	Should the routine write to console?
[in]	eFitFn	Fitting 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]	fix	Constrains 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_fitted	If 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] effn Index

Returns: Description or NULL

int effnNparams ( int effn )

Returns number of function parameters associated with a fitting function.

Parameters

[in] effn Index

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]	fitfn	Fitting function (0 .. effnNR)
[in]	oenv	Output formatting information
[in]	bVerbose	Should the routine write to console?
[in]	tbeginfit	Starting time for fitting
[in]	tendfit	Ending time for fitting
[in]	dt	The time step
[in]	c1	The data points
[in,out]	fit	The 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]	eFitFn	the index to the fitting function (0 .. effnNR)
[in]	parm	Array of parameters, the length of which depends on eFitFn
[in]	x	The 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] sffn Two dimensional string array coming from parse_common_args

Returns: the ffn enum

Variable Documentation

t_lmcurve lmcurves[effnNR+1]

Initial value:

= { lmc_exp_one_parm,  lmc_exp_one_parm, lmc_exp_two_parm,
                                   lmc_exp_exp,       lmc_exp_5_parm,   lmc_exp_7_parm,
                                   lmc_exp_9_parm,    lmc_vac_2_parm,   lmc_erffit,
                                   lmc_errest_3_parm, lmc_pres_6_parm }

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.

Description

Functions

Variables

Function Documentation

Variable Documentation