2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
40 * Implements routine for fitting a data set to a curve
42 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
43 * \ingroup module_correlationfunctions
52 #include "external/lmfit/lmcurve.h"
54 #include "gromacs/correlationfunctions/integrate.h"
55 #include "gromacs/fileio/xvgr.h"
56 #include "gromacs/legacyheaders/macros.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/real.h"
61 #include "gromacs/utility/smalloc.h"
63 /*! \brief Number of parameters for each fitting function */
64 static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3, 6};
66 const char *s_ffn[effnNR+2] = {
67 NULL, "none", "exp", "aexp", "exp_exp", "vac",
68 "exp5", "exp7", "exp9", "erffit", "effnPRES", NULL, NULL
70 /* We don't allow errest as a choice on the command line */
72 /*! \brief Long description for each fitting function type */
73 static const char *longs_ffn[effnNR] = {
77 "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
78 "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
79 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5",
80 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
81 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
82 "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
83 "y = a2 *2*a1*((exp(-x/a1)-1)*(a1/x)+1)+(1-a2)*2*a3*((exp(-x/a3)-1)*(a3/x)+1)",
84 "y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6)"
87 int effnNparams(int effn)
89 if ((0 <= effn) && (effn < effnNR))
99 const char *effnDescription(int effn)
101 if ((0 <= effn) && (effn < effnNR))
103 return longs_ffn[effn];
111 int sffn2effn(const char **sffn)
116 for (i = 0; i < effnNR; i++)
118 if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
127 /*! \brief Compute exponential function A exp(-x/tau) */
128 static double myexp(double x, double A, double tau)
130 if ((A == 0) || (tau == 0))
134 return A*exp(-x/tau);
137 /*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
138 static double lmc_erffit (double x, const double *a)
142 erfarg = (x-a[2])/(a[3]*a[3]);
143 return (a[0]+a[1])/2-(a[0]-a[1])*gmx_erf(erfarg);
146 /*! \brief Compute y = exp(-x/a0) */
147 static double lmc_exp_one_parm(double x, const double *a)
152 /*! \brief Compute y = a1 exp(-x/a0) */
153 static double lmc_exp_two_parm(double x, const double *a)
155 return a[1]*exp(-x/a[0]);
158 /*! \brief Compute y = a1 exp(-x/a0) + (1-a1) exp(-x/a2) */
159 static double lmc_exp_3_parm(double x, const double *a)
165 return a[1]*e1 + (1-a[1])*e2;
168 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 */
169 static double lmc_exp_5_parm(double x, const double *a)
175 return a[0]*e1 + a[2]*e2 + a[4];
178 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 */
179 static double lmc_exp_7_parm(double x, const double *a)
186 return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6];
189 /*! \brief Compute y = a0 exp(-x/a1) + a2 exp(-x/a3) + a4 exp(-x/a5) + a6 exp(-x/a7) + a8 */
190 static double lmc_exp_9_parm(double x, const double *a)
192 double e1, e2, e3, e4;
198 return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8];
201 /*! \brief Compute y = (1-a1) * exp(-(x/a3)^a4)*cos(x*a2) + a1*exp(-(x/a5)^a6) */
202 static double effnPRES_fun(double x, const double *a)
204 double term1, term2, term3;
208 term3 = a[0] * exp(-pow((x/a[4]), a[5]));
218 term2 = exp(-pow((x/a[2]), a[3])) * cos(x*a[1]);
225 return term1*term2 + term3;
228 /*! \brief Compute vac function */
229 static double lmc_vac_2_parm(double x, const double *a)
233 * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
235 * = exp(-v) (cosh(wv) + 1/w sinh(wv))
240 * For tranverse current autocorrelation functions:
242 * a2 = 4 tau (eta/rho) k^2
246 double y, v, det, omega, omega2, em, ec, es;
254 omega = sqrt(omega2);
257 ec = em*0.5*(exp(omega*v)+exp(-omega*v));
258 es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
262 ec = em*cos(omega*v);
263 es = em*sin(omega*v)/omega;
274 /*! \brief Compute error estimate */
275 static double lmc_errest_3_parm(double x, const double *a)
278 e1 = (exp(-x/a1) - 1)
279 e2 = (exp(-x/a3) - 1)
280 v1= 2*a1 * (e1*a1/x + 1)
281 v2 = 2*a3 * (e2*a3/x + 1)
282 fun = a2*v1 + (1 - a2) * v2
284 double e1, e2, v1, v2;
288 e1 = gmx_expm1(-x/a[0]);
296 e2 = gmx_expm1(-x/a[2]);
305 v1 = 2*a[0]*(e1*a[0]/x + 1);
306 v2 = 2*a[2]*(e2*a[2]/x + 1);
307 return a[1]*v1 + (1-a[1])*v2;
315 /*! \brief function type for passing to fitting routine */
316 typedef double (*t_lmcurve)(double x, const double *a);
318 /*! \brief array of fitting functions corresponding to the pre-defined types */
319 t_lmcurve lmcurves[effnNR] = {
320 lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
321 lmc_exp_3_parm, lmc_vac_2_parm,
322 lmc_exp_5_parm, lmc_exp_7_parm,
323 lmc_exp_9_parm, lmc_erffit, lmc_errest_3_parm, effnPRES_fun
326 double fit_function(int eFitFn, double *parm, double x)
328 // TODO: check that eFitFn is within bounds
329 return lmcurves[eFitFn](x, parm);
332 /*! \brief lmfit_exp supports up to 3 parameter fitting of exponential functions */
333 static gmx_bool lmfit_exp(int nfit, double x[], double y[],
334 real ftol, int maxiter,
335 double parm[], gmx_bool bVerbose,
338 double chisq, ochisq;
341 lm_control_struct *control;
342 lm_status_struct *status;
344 if ((eFitFn < 0) || (eFitFn >= effnNR))
346 gmx_fatal(FARGS, "fitfn = %d, should be in 0..%d (%s,%d)",
347 effnNR-1, eFitFn, __FILE__, __LINE__);
351 control->ftol = ftol;
352 control->xtol = ftol;
353 control->gtol = ftol;
354 control->epsilon = 0.1;
355 control->stepbound = 100;
356 control->patience = maxiter;
357 control->scale_diag = 1;
358 control->msgfile = stdout;
359 control->verbosity = (bVerbose ? 1 : 0);
360 control->n_maxpri = nfp_ffn[eFitFn];
361 control->m_maxpri = 0;
368 fprintf(stderr, "%4s %10s Parameters\n",
375 lmcurve(nfp_ffn[eFitFn], parm, nfit, x, y,
376 lmcurves[eFitFn], control, status);
377 chisq = sqr(status->fnorm);
381 fprintf(stderr, "%4d %10.5e", j, chisq);
382 for (mmm = 0; (mmm < nfp_ffn[eFitFn]); mmm++)
384 fprintf(stderr, " %10.5e", parm[mmm]);
386 fprintf(stderr, "\n");
389 bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
390 ((ochisq == chisq)));
392 while (bCont && (j < maxiter));
395 fprintf(stderr, "\n");
404 /*! \brief See description in header file. */
405 real do_lmfit(int ndata, real c1[], real sig[], real dt, real x0[],
406 real begintimefit, real endtimefit, const output_env_t oenv,
407 gmx_bool bVerbose, int eFitFn, double fitparms[], int fix)
412 int i, j, nparm, nfitpnts;
413 double integral, ttt;
414 double *parm, *dparm;
425 fprintf(stderr, "Using fixed parameters in curve fitting is not working anymore\n");
427 nparm = nfp_ffn[eFitFn];
430 fprintf(debug, "There are %d points to fit %d vars!\n", ndata, nparm);
431 fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
432 eFitFn, begintimefit, endtimefit, dt);
440 for (i = 0; i < ndata; i++)
442 ttt = x0 ? x0[i] : dt*i;
443 if (ttt >= begintimefit && ttt <= endtimefit)
448 /* mrqmin does not like sig to be zero */
459 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy= %g\n",
460 j, i, x[j], y[j], dy[j]);
467 if (nfitpnts < nparm)
469 fprintf(stderr, "Not enough data points for fitting!\n");
477 for (i = 0; (i < nparm); i++)
479 parm[i] = fitparms[i];
483 if (!lmfit_exp(nfitpnts, x, y, ftol, maxiter, parm, bVerbose, eFitFn))
485 fprintf(stderr, "Fit failed!\n");
489 /* Compute the integral from begintimefit */
492 integral = (parm[0]*myexp(begintimefit, parm[1], parm[0]) +
493 parm[2]*myexp(begintimefit, 1-parm[1], parm[2]));
497 integral = parm[0]*myexp(begintimefit, parm[1], parm[0]);
501 integral = parm[0]*myexp(begintimefit, 1, parm[0]);
505 gmx_fatal(FARGS, "nparm = %d in file %s, line %d",
506 nparm, __FILE__, __LINE__);
509 /* Generate THE output */
512 fprintf(stderr, "FIT: # points used in fit is: %d\n", nfitpnts);
513 fprintf(stderr, "FIT: %21s%21s%21s\n",
514 "parm0 ", "parm1 (ps) ", "parm2 (ps) ");
515 fprintf(stderr, "FIT: ------------------------------------------------------------\n");
516 fprintf(stderr, "FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
517 parm[0], dparm[0], parm[1], dparm[1], parm[2], dparm[2]);
518 fprintf(stderr, "FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
519 begintimefit, integral);
521 sprintf(buf, "test%d.xvg", nfitpnts);
522 fp = xvgropen(buf, "C(t) + Fit to C(t)", "Time (ps)", "C(t)", oenv);
523 fprintf(fp, "# parm0 = %g, parm1 = %g, parm2 = %g\n",
524 parm[0], parm[1], parm[2]);
525 for (j = 0; (j < nfitpnts); j++)
527 ttt = x0 ? x0[j] : dt*j;
528 fprintf(fp, "%10.5e %10.5e %10.5e\n",
530 lmcurves[eFitFn](ttt, parm));
537 for (i = 0; (i < nparm); i++)
539 fitparms[i] = parm[i];
553 /*! See description in header file. */
554 real fit_acf(int ncorr, int fitfn, const output_env_t oenv, gmx_bool bVerbose,
555 real tbeginfit, real tendfit, real dt, real c1[], real *fit)
558 double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
560 int i, j, jmax, nf_int;
563 bPrint = bVerbose || bDebugMode();
574 nf_int = min(ncorr, (int)(tendfit/dt));
575 sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
577 /* Estimate the correlation time for better fitting */
578 ct_estimate = 0.5*c1[0];
579 for (i = 1; (i < ncorr) && (c1[i] > 0); i++)
581 ct_estimate += c1[i];
583 ct_estimate *= dt/c1[0];
587 printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n",
588 0.0, dt*nf_int, sum);
589 printf("COR: Relaxation times are computed as fit to an exponential:\n");
590 printf("COR: %s\n", longs_ffn[fitfn]);
591 printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n", tbeginfit, min(ncorr*dt, tendfit));
597 printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
598 "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
599 (nfp_ffn[fitfn] >= 2) ? " a2 ()" : "",
600 (nfp_ffn[fitfn] >= 3) ? " a3 (ps)" : "");
613 for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr*dt)); j++)
615 /* Estimate the correlation time for better fitting */
618 for (i = 0; (i < ncorr) && (dt*i < tStart || c1[i] > 0); i++)
625 ct_estimate = 0.5*c1[i];
630 ct_estimate += c1[i];
635 ct_estimate *= dt/c_start;
639 /* The data is strange, so we need to choose somehting */
640 ct_estimate = tendfit;
644 fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
647 if (fitfn == effnEXP3)
649 fitparm[0] = 0.002*ncorr*dt;
651 fitparm[2] = 0.2*ncorr*dt;
655 /* Good initial guess, this increases the probability of convergence */
656 fitparm[0] = ct_estimate;
661 /* Generate more or less appropriate sigma's */
662 for (i = 0; i < ncorr; i++)
664 sig[i] = sqrt(ct_estimate+dt*i);
667 nf_int = min(ncorr, (int)((tStart+1e-4)/dt));
668 sum = print_and_integrate(debug, nf_int, dt, c1, NULL, 1);
669 tail_corr = do_lmfit(ncorr, c1, sig, dt, NULL, tStart, tendfit, oenv,
670 bDebugMode(), fitfn, fitparm, 0);
671 sumtot = sum+tail_corr;
672 if (fit && ((jmax == 1) || (j == 1)))
675 for (i = 0; (i < 3); i++)
679 for (i = 0; (i < ncorr); i++)
681 fit[i] = lmcurves[fitfn](i*dt, mfp);
686 printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
687 for (i = 0; (i < nfp_ffn[fitfn]); i++)
689 printf(" %11.4e", fitparm[i]);