Gromacs
2022.2
|
#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"
Implements routine for fitting a data set to a curve.
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... | |
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.
[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 |
const char* effnDescription | ( | int | effn | ) |
Returns description corresponding to the enum above, or NULL if out of range.
[in] | effn | Index |
int effnNparams | ( | int | effn | ) |
Returns number of function parameters associated with a fitting function.
[in] | effn | Index |
|
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.
[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 |
double fit_function | ( | int | eFitFn, |
const double | parm[], | ||
double | x | ||
) |
Returns the value of fit function eFitFn at x.
[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 |
|
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 |
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 |
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.
[in] | sffn | Two dimensional string array coming from parse_common_args |
t_lmcurve lmcurves[effnNR+1] |
array of fitting functions corresponding to the pre-defined types
|
static |
Long description for each fitting function type.
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.