1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013-2019, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37 /*! \internal
38 * \file
39 * \brief
40 * Implements routine for fitting a data set to a curve
41 *
42 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
43 * \ingroup module_correlationfunctions
44 */
45 #include "gmxpre.h"
46
47 #include "expfit.h"
48
49 #include <cstring>
50
51 #include <algorithm>
52
53 #include "gromacs/correlationfunctions/integrate.h"
54 #include "gromacs/fileio/xvgr.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/gmxassert.h"
59 #include "gromacs/utility/real.h"
60 #include "gromacs/utility/smalloc.h"
61
62 #include "gmx_lmcurve.h"
63
64 /*! \brief Number of parameters for each fitting function */
65 static const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 5, 7, 9, 2, 4, 3, 6 };
66
67 /* +2 becuase parse_common_args wants leading and concluding NULL.
68 * We only allow exponential functions as choices on the command line,
69 * hence there are many more NULL field (which have to be at the end of
70 * the array).
71 */
72 const char* s_ffn[effnNR + 2] = { nullptr, "none", "exp", "aexp", "exp_exp", "exp5", "exp7",
73 "exp9", nullptr, nullptr, nullptr, nullptr, nullptr };
74
75 // clang-format off
76 // needed because clang-format wants to break the lines below
77 /*! \brief Long description for each fitting function type */
78 static const char* longs_ffn[effnNR]
79 = { "no fit",
80 "y = exp(-x/|a0|)",
81 "y = a1 exp(-x/|a0|)",
82 "y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0",
83 "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4, a3 >= a1",
84 "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1",
85 "y = a0 exp(-x/|a1|) + a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1",
86 "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)",
87 "y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)",
88 "y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)",
89 "y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)" };
90 // clang-format on
91
effnNparams(int effn)92 int effnNparams(int effn)
93 {
94 if ((0 <= effn) && (effn < effnNR))
95 {
96 return nfp_ffn[effn];
97 }
98 else
99 {
100 return -1;
101 }
102 }
103
effnDescription(int effn)104 const char* effnDescription(int effn)
105 {
106 if ((0 <= effn) && (effn < effnNR))
107 {
108 return longs_ffn[effn];
109 }
110 else
111 {
112 return nullptr;
113 }
114 }
115
sffn2effn(const char ** sffn)116 int sffn2effn(const char** sffn)
117 {
118 int eFitFn, i;
119
120 eFitFn = 0;
121 for (i = 0; i < effnNR; i++)
122 {
123 if (sffn[i + 1] && strcmp(sffn[0], sffn[i + 1]) == 0)
124 {
125 eFitFn = i;
126 }
127 }
128
129 return eFitFn;
130 }
131
132 /*! \brief Compute exponential function A exp(-x/tau) */
myexp(double x,double A,double tau)133 static double myexp(double x, double A, double tau)
134 {
135 if ((A == 0) || (tau == 0))
136 {
137 return 0;
138 }
139 return A * exp(-x / tau);
140 }
141
142 /*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
lmc_erffit(double x,const double * a)143 static double lmc_erffit(double x, const double* a)
144 {
145 double erfarg;
146 double myerf;
147
148 if (a[3] != 0)
149 {
150 erfarg = (x - a[2]) / (a[3] * a[3]);
151 myerf = std::erf(erfarg);
152 }
153 else
154 {
155 /* If a[3] == 0, a[3]^2 = 0 and the erfarg becomes +/- infinity */
156 if (x < a[2])
157 {
158 myerf = -1;
159 }
160 else
161 {
162 myerf = 1;
163 }
164 }
165 return 0.5 * ((a[0] + a[1]) - (a[0] - a[1]) * myerf);
166 }
167
168 /*! \brief Exponent function that prevents overflow */
safe_exp(double x)169 static double safe_exp(double x)
170 {
171 double exp_max = 200;
172 double exp_min = -exp_max;
173 if (x <= exp_min)
174 {
175 return exp(exp_min);
176 }
177 else if (x >= exp_max)
178 {
179 return exp(exp_max);
180 }
181 else
182 {
183 return exp(x);
184 }
185 }
186
187 /*! \brief Exponent minus 1 function that prevents overflow */
safe_expm1(double x)188 static double safe_expm1(double x)
189 {
190 double exp_max = 200;
191 double exp_min = -exp_max;
192 if (x <= exp_min)
193 {
194 return -1;
195 }
196 else if (x >= exp_max)
197 {
198 return exp(exp_max);
199 }
200 else
201 {
202 return std::expm1(x);
203 }
204 }
205
206 /*! \brief Compute y = exp(-x/|a0|) */
lmc_exp_one_parm(double x,const double * a)207 static double lmc_exp_one_parm(double x, const double* a)
208 {
209 return safe_exp(-x / fabs(a[0]));
210 }
211
212 /*! \brief Compute y = a1 exp(-x/|a0|) */
lmc_exp_two_parm(double x,const double * a)213 static double lmc_exp_two_parm(double x, const double* a)
214 {
215 return a[1] * safe_exp(-x / fabs(a[0]));
216 }
217
218 /*! \brief Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|) */
lmc_exp_exp(double x,const double * a)219 static double lmc_exp_exp(double x, const double* a)
220 {
221 double e1, e2;
222
223 e1 = safe_exp(-x / fabs(a[0]));
224 e2 = safe_exp(-x / (fabs(a[0]) + fabs(a[2])));
225 return a[1] * e1 + (1 - a[1]) * e2;
226 }
227
228 /*! \brief Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 */
lmc_exp_5_parm(double x,const double * a)229 static double lmc_exp_5_parm(double x, const double* a)
230 {
231 double e1, e2;
232
233 e1 = safe_exp(-x / fabs(a[1]));
234 e2 = safe_exp(-x / (fabs(a[1]) + fabs(a[3])));
235 return a[0] * e1 + a[2] * e2 + a[4];
236 }
237
238 /*! \brief Compute 7 parameter exponential function value.
239 *
240 * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
241 * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6
242 */
lmc_exp_7_parm(double x,const double * a)243 static double lmc_exp_7_parm(double x, const double* a)
244 {
245 double e1, e2, e3;
246 double fa1, fa3, fa5;
247
248 fa1 = fabs(a[1]);
249 fa3 = fa1 + fabs(a[3]);
250 fa5 = fa3 + fabs(a[5]);
251 e1 = safe_exp(-x / fa1);
252 e2 = safe_exp(-x / fa3);
253 e3 = safe_exp(-x / fa5);
254 return a[0] * e1 + a[2] * e2 + a[4] * e3 + a[6];
255 }
256
257 /*! \brief Compute 9 parameter exponential function value.
258 *
259 * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
260 * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8
261 */
lmc_exp_9_parm(double x,const double * a)262 static double lmc_exp_9_parm(double x, const double* a)
263 {
264 double e1, e2, e3, e4;
265 double fa1, fa3, fa5, fa7;
266
267 fa1 = fabs(a[1]);
268 fa3 = fa1 + fabs(a[3]);
269 fa5 = fa3 + fabs(a[5]);
270 fa7 = fa5 + fabs(a[7]);
271
272 e1 = safe_exp(-x / fa1);
273 e2 = safe_exp(-x / fa3);
274 e3 = safe_exp(-x / fa5);
275 e4 = safe_exp(-x / fa7);
276 return a[0] * e1 + a[2] * e2 + a[4] * e3 + a[6] * e4 + a[8];
277 }
278
279 /*! \brief Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|) */
lmc_pres_6_parm(double x,const double * a)280 static double lmc_pres_6_parm(double x, const double* a)
281 {
282 double term1, term2, term3;
283 double pow_max = 10;
284
285 term3 = 0;
286 if ((a[4] != 0) && (a[0] != 0))
287 {
288 double power = std::min(fabs(a[5]), pow_max);
289 term3 = a[0] * safe_exp(-pow((x / fabs(a[4])), power));
290 }
291
292 term1 = 1 - a[0];
293 term2 = 0;
294 if ((term1 != 0) && (a[2] != 0))
295 {
296 double power = std::min(fabs(a[3]), pow_max);
297 term2 = safe_exp(-pow((x / fabs(a[2])), power)) * cos(x * fabs(a[1]));
298 }
299
300 return term1 * term2 + term3;
301 }
302
303 /*! \brief Compute vac function */
lmc_vac_2_parm(double x,const double * a)304 static double lmc_vac_2_parm(double x, const double* a)
305 {
306 /* Fit to function
307 *
308 * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
309 *
310 * = exp(-v) (cosh(wv) + 1/w sinh(wv))
311 *
312 * v = x/(2 a0)
313 * w = sqrt(1 - a1)
314 *
315 * For tranverse current autocorrelation functions:
316 * a0 = tau
317 * a1 = 4 tau (eta/rho) k^2
318 *
319 */
320
321 double y, v, det, omega, wv, em, ec, es;
322 double wv_max = 100;
323
324 v = x / (2 * fabs(a[0]));
325 det = 1 - a[1];
326 em = safe_exp(-v);
327 if (det != 0)
328 {
329 omega = sqrt(fabs(det));
330 wv = std::min(omega * v, wv_max);
331
332 if (det > 0)
333 {
334 ec = em * 0.5 * (safe_exp(wv) + safe_exp(-wv));
335 es = em * 0.5 * (safe_exp(wv) - safe_exp(-wv)) / omega;
336 }
337 else
338 {
339 ec = em * cos(wv);
340 es = em * sin(wv) / omega;
341 }
342 y = ec + es;
343 }
344 else
345 {
346 y = (1 + v) * em;
347 }
348 return y;
349 }
350
351 /*! \brief Compute error estimate */
lmc_errest_3_parm(double x,const double * a)352 static double lmc_errest_3_parm(double x, const double* a)
353 {
354 double e1, e2, v1, v2;
355 double fa0 = fabs(a[0]);
356 double fa1;
357 double fa2 = fa0 + fabs(a[2]);
358
359 if (a[0] != 0)
360 {
361 e1 = safe_expm1(-x / fa0);
362 }
363 else
364 {
365 e1 = 0;
366 }
367 if (a[2] != 0)
368 {
369 e2 = safe_expm1(-x / fa2);
370 }
371 else
372 {
373 e2 = 0;
374 }
375
376 if (x > 0)
377 {
378 v1 = 2 * fa0 * (e1 * fa0 / x + 1);
379 v2 = 2 * fa2 * (e2 * fa2 / x + 1);
380 /* We need 0 <= a1 <= 1 */
381 fa1 = std::min(1.0, std::max(0.0, a[1]));
382
383 return fa1 * v1 + (1 - fa1) * v2;
384 }
385 else
386 {
387 return 0;
388 }
389 }
390
391 /*! \brief array of fitting functions corresponding to the pre-defined types */
392 t_lmcurve lmcurves[effnNR + 1] = { lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
393 lmc_exp_exp, lmc_exp_5_parm, lmc_exp_7_parm,
394 lmc_exp_9_parm, lmc_vac_2_parm, lmc_erffit,
395 lmc_errest_3_parm, lmc_pres_6_parm };
396
fit_function(const int eFitFn,const double parm[],const double x)397 double fit_function(const int eFitFn, const double parm[], const double x)
398 {
399 if ((eFitFn < 0) || (eFitFn >= effnNR))
400 {
401 fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n", eFitFn, effnNR - 1);
402 return 0.0;
403 }
404 return lmcurves[eFitFn](x, parm);
405 }
406
407 /*! \brief Ensure the fitting parameters are well-behaved.
408 *
409 * In order to make sure that fitting behaves according to the
410 * constraint that time constants are positive and increasing
411 * we enforce this here by setting all time constants to their
412 * absolute value and by adding e.g. |a_0| to |a_2|. This is
413 * done in conjunction with the fitting functions themselves.
414 * When there are multiple time constants we make sure that
415 * the successive time constants are at least double the
416 * previous ones and during fitting we enforce the they remain larger.
417 * It may very well help the convergence of the fitting routine.
418 */
initiate_fit_params(int eFitFn,double params[])419 static void initiate_fit_params(int eFitFn, double params[])
420 {
421 int i, nparm;
422
423 nparm = effnNparams(eFitFn);
424
425 switch (eFitFn)
426 {
427 case effnVAC: GMX_ASSERT(params[0] >= 0, "parameters should be >= 0"); break;
428 case effnEXP1:
429 case effnEXP2:
430 case effnEXPEXP:
431 GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
432 if (nparm > 2)
433 {
434 GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
435 /* In order to maintain params[2] >= params[0] in the final
436 * result, we fit the difference between the two, that is
437 * params[2]-params[0] and in the add add in params[0]
438 * again.
439 */
440 params[2] = std::max(fabs(params[2]) - params[0], params[0]);
441 }
442 break;
443 case effnEXP5:
444 case effnEXP7:
445 case effnEXP9:
446 GMX_ASSERT(params[1] >= 0, "parameters should be >= 0");
447 params[1] = fabs(params[1]);
448 if (nparm > 3)
449 {
450 GMX_ASSERT(params[3] >= 0, "parameters should be >= 0");
451 /* See comment under effnEXPEXP */
452 params[3] = std::max(fabs(params[3]) - params[1], params[1]);
453 if (nparm > 5)
454 {
455 GMX_ASSERT(params[5] >= 0, "parameters should be >= 0");
456 /* See comment under effnEXPEXP */
457 params[5] = std::max(fabs(params[5]) - params[3], params[3]);
458 if (nparm > 7)
459 {
460 GMX_ASSERT(params[7] >= 0, "parameters should be >= 0");
461 /* See comment under effnEXPEXP */
462 params[7] = std::max(fabs(params[7]) - params[5], params[5]);
463 }
464 }
465 }
466 break;
467 case effnERREST:
468 GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
469 GMX_ASSERT(params[1] >= 0 && params[1] <= 1, "parameter 1 should in 0 .. 1");
470 GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
471 /* See comment under effnEXPEXP */
472 params[2] = fabs(params[2]) - params[0];
473 break;
474 case effnPRES:
475 for (i = 1; (i < nparm); i++)
476 {
477 GMX_ASSERT(params[i] >= 0, "parameters should be >= 0");
478 }
479 break;
480 default: break;
481 }
482 }
483
484 /*! \brief Process the fitting parameters to get output parameters.
485 *
486 * See comment at the previous function.
487 */
extract_fit_params(int eFitFn,double params[])488 static void extract_fit_params(int eFitFn, double params[])
489 {
490 int i, nparm;
491
492 nparm = effnNparams(eFitFn);
493
494 switch (eFitFn)
495 {
496 case effnVAC: params[0] = fabs(params[0]); break;
497 case effnEXP1:
498 case effnEXP2:
499 case effnEXPEXP:
500 params[0] = fabs(params[0]);
501 if (nparm > 2)
502 {
503 /* Back conversion of parameters from the fitted difference
504 * to the absolute value.
505 */
506 params[2] = fabs(params[2]) + params[0];
507 }
508 break;
509 case effnEXP5:
510 case effnEXP7:
511 case effnEXP9:
512 params[1] = fabs(params[1]);
513 if (nparm > 3)
514 {
515 /* See comment under effnEXPEXP */
516 params[3] = fabs(params[3]) + params[1];
517 if (nparm > 5)
518 {
519 /* See comment under effnEXPEXP */
520 params[5] = fabs(params[5]) + params[3];
521 if (nparm > 7)
522 {
523 /* See comment under effnEXPEXP */
524 params[7] = fabs(params[7]) + params[5];
525 }
526 }
527 }
528 break;
529 case effnERREST:
530 params[0] = fabs(params[0]);
531 if (params[1] < 0)
532 {
533 params[1] = 0;
534 }
535 else if (params[1] > 1)
536 {
537 params[1] = 1;
538 }
539 /* See comment under effnEXPEXP */
540 params[2] = params[0] + fabs(params[2]);
541 break;
542 case effnPRES:
543 for (i = 1; (i < nparm); i++)
544 {
545 params[i] = fabs(params[i]);
546 }
547 break;
548 default: break;
549 }
550 }
551
552 /*! \brief Print chi-squared value and the parameters */
print_chi2_params(FILE * fp,const int eFitFn,const double fitparms[],const char * label,const int nfitpnts,const double x[],const double y[])553 static void print_chi2_params(FILE* fp,
554 const int eFitFn,
555 const double fitparms[],
556 const char* label,
557 const int nfitpnts,
558 const double x[],
559 const double y[])
560 {
561 int i;
562 double chi2 = 0;
563
564 for (i = 0; (i < nfitpnts); i++)
565 {
566 double yfit = lmcurves[eFitFn](x[i], fitparms);
567 chi2 += gmx::square(y[i] - yfit);
568 }
569 fprintf(fp, "There are %d data points, %d parameters, %s chi2 = %g\nparams:", nfitpnts,
570 effnNparams(eFitFn), label, chi2);
571 for (i = 0; (i < effnNparams(eFitFn)); i++)
572 {
573 fprintf(fp, " %10g", fitparms[i]);
574 }
575 fprintf(fp, "\n");
576 }
577
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)578 real do_lmfit(int ndata,
579 const real c1[],
580 real sig[],
581 real dt,
582 const real* x0,
583 real begintimefit,
584 real endtimefit,
585 const gmx_output_env_t* oenv,
586 gmx_bool bVerbose,
587 int eFitFn,
588 double fitparms[],
589 int fix,
590 const char* fn_fitted)
591 {
592 FILE* fp;
593 int i, j, nfitpnts;
594 double integral, ttt;
595 double *x, *y, *dy;
596
597 if (0 != fix)
598 {
599 fprintf(stderr, "Using fixed parameters in curve fitting is temporarily not working.\n");
600 }
601 if (debug)
602 {
603 fprintf(debug, "There are %d points to fit %d vars!\n", ndata, effnNparams(eFitFn));
604 fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n", eFitFn, begintimefit,
605 endtimefit, dt);
606 }
607
608 snew(x, ndata);
609 snew(y, ndata);
610 snew(dy, ndata);
611 j = 0;
612 for (i = 0; i < ndata; i++)
613 {
614 ttt = x0 ? x0[i] : dt * i;
615 if (ttt >= begintimefit && ttt <= endtimefit)
616 {
617 x[j] = ttt;
618 y[j] = c1[i];
619 if (nullptr == sig)
620 {
621 // No weighting if all values are divided by 1.
622 dy[j] = 1;
623 }
624 else
625 {
626 dy[j] = std::max(1.0e-7, static_cast<double>(sig[i]));
627 }
628 if (debug)
629 {
630 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy=%g, ttt=%g\n", j, i, x[j], y[j], dy[j], ttt);
631 }
632 j++;
633 }
634 }
635 nfitpnts = j;
636 integral = 0;
637 if (nfitpnts < effnNparams(eFitFn))
638 {
639 fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!\n", nfitpnts, dt);
640 }
641 else
642 {
643 gmx_bool bSuccess;
644
645 if (bVerbose)
646 {
647 print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y);
648 }
649 initiate_fit_params(eFitFn, fitparms);
650
651 bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix);
652 extract_fit_params(eFitFn, fitparms);
653
654 if (!bSuccess)
655 {
656 fprintf(stderr, "Fit failed!\n");
657 }
658 else
659 {
660 if (bVerbose)
661 {
662 print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y);
663 }
664 switch (eFitFn)
665 {
666 case effnEXP1: integral = fitparms[0] * myexp(begintimefit, 1, fitparms[0]); break;
667 case effnEXP2:
668 integral = fitparms[0] * myexp(begintimefit, fitparms[1], fitparms[0]);
669 break;
670 case effnEXPEXP:
671 integral = (fitparms[0] * myexp(begintimefit, fitparms[1], fitparms[0])
672 + fitparms[2] * myexp(begintimefit, 1 - fitparms[1], fitparms[2]));
673 break;
674 case effnEXP5:
675 case effnEXP7:
676 case effnEXP9:
677 integral = 0;
678 for (i = 0; (i < (effnNparams(eFitFn) - 1) / 2); i++)
679 {
680 integral += fitparms[2 * i]
681 * myexp(begintimefit, fitparms[2 * i + 1], fitparms[2 * i]);
682 }
683 break;
684 default:
685 /* Do numerical integration */
686 integral = 0;
687 for (i = 0; (i < nfitpnts - 1); i++)
688 {
689 double y0 = lmcurves[eFitFn](x[i], fitparms);
690 double y1 = lmcurves[eFitFn](x[i + 1], fitparms);
691 integral += (x[i + 1] - x[i]) * (y1 + y0) * 0.5;
692 }
693 break;
694 }
695
696 if (bVerbose)
697 {
698 printf("FIT: Integral of fitted function: %g\n", integral);
699 if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn))
700 {
701 printf("FIT: Note that the constant term is not taken into account when "
702 "computing integral.\n");
703 }
704 }
705 /* Generate debug output */
706 if (nullptr != fn_fitted)
707 {
708 fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)", "Data (t)", oenv);
709 for (i = 0; (i < effnNparams(eFitFn)); i++)
710 {
711 fprintf(fp, "# fitparms[%d] = %g\n", i, fitparms[i]);
712 }
713 for (j = 0; (j < nfitpnts); j++)
714 {
715 real ttt = x0 ? x0[j] : dt * j;
716 fprintf(fp, "%10.5e %10.5e %10.5e\n", x[j], y[j], (lmcurves[eFitFn])(ttt, fitparms));
717 }
718 xvgrclose(fp);
719 }
720 }
721 }
722
723 sfree(x);
724 sfree(y);
725 sfree(dy);
726
727 return integral;
728 }
729
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)730 real fit_acf(int ncorr,
731 int fitfn,
732 const gmx_output_env_t* oenv,
733 gmx_bool bVerbose,
734 real tbeginfit,
735 real tendfit,
736 real dt,
737 real c1[],
738 real* fit)
739 {
740 double fitparm[3];
741 double tStart, tail_corr, sum, sumtot = 0, c_start, ct_estimate;
742 real* sig;
743 int i, j, jmax, nf_int;
744 gmx_bool bPrint;
745
746 GMX_ASSERT(effnNparams(fitfn) < static_cast<int>(sizeof(fitparm) / sizeof(double)),
747 "Fitting function with more than 3 parameters not supported!");
748 bPrint = bVerbose || bDebugMode();
749
750 if (bPrint)
751 {
752 printf("COR:\n");
753 }
754
755 if (tendfit <= 0)
756 {
757 tendfit = ncorr * dt;
758 }
759 nf_int = std::min(ncorr, static_cast<int>(tendfit / dt));
760 sum = print_and_integrate(debug, nf_int, dt, c1, nullptr, 1);
761
762 if (bPrint)
763 {
764 printf("COR: Correlation time (plain integral from %6.3f to %6.3f ps) = %8.5f ps\n", 0.0,
765 dt * nf_int, sum);
766 printf("COR: Relaxation times are computed as fit to an exponential:\n");
767 printf("COR: %s\n", effnDescription(fitfn));
768 printf("COR: Fit to correlation function from %6.3f ps to %6.3f ps, results in a\n",
769 tbeginfit, std::min(ncorr * dt, tendfit));
770 }
771
772 tStart = 0;
773 if (bPrint)
774 {
775 printf("COR:%11s%11s%11s%11s%11s%11s%11s\n", "Fit from", "Integral", "Tail Value",
776 "Sum (ps)", " a1 (ps)", (effnNparams(fitfn) >= 2) ? " a2 ()" : "",
777 (effnNparams(fitfn) >= 3) ? " a3 (ps)" : "");
778 }
779
780 snew(sig, ncorr);
781
782 if (tbeginfit > 0)
783 {
784 jmax = 3;
785 }
786 else
787 {
788 jmax = 1;
789 }
790 for (j = 0; ((j < jmax) && (tStart < tendfit) && (tStart < ncorr * dt)); j++)
791 {
792 /* Estimate the correlation time for better fitting */
793 c_start = -1;
794 ct_estimate = 0;
795 for (i = 0; (i < ncorr) && (dt * i < tStart || c1[i] > 0); i++)
796 {
797 if (c_start < 0)
798 {
799 if (dt * i >= tStart)
800 {
801 c_start = c1[i];
802 ct_estimate = 0.5 * c1[i];
803 }
804 }
805 else
806 {
807 ct_estimate += c1[i];
808 }
809 }
810 if (c_start > 0)
811 {
812 ct_estimate *= dt / c_start;
813 }
814 else
815 {
816 /* The data is strange, so we need to choose somehting */
817 ct_estimate = tendfit;
818 }
819 if (debug)
820 {
821 fprintf(debug, "tStart %g ct_estimate: %g\n", tStart, ct_estimate);
822 }
823
824 if (fitfn == effnEXPEXP)
825 {
826 fitparm[0] = 0.002 * ncorr * dt;
827 fitparm[1] = 0.95;
828 fitparm[2] = 0.2 * ncorr * dt;
829 }
830 else
831 {
832 /* Good initial guess, this increases the probability of convergence */
833 fitparm[0] = ct_estimate;
834 fitparm[1] = 1.0;
835 fitparm[2] = 1.0;
836 }
837
838 /* Generate more or less appropriate sigma's */
839 for (i = 0; i < ncorr; i++)
840 {
841 sig[i] = sqrt(ct_estimate + dt * i);
842 }
843
844 nf_int = std::min(ncorr, static_cast<int>((tStart + 1e-4) / dt));
845 sum = print_and_integrate(debug, nf_int, dt, c1, nullptr, 1);
846 tail_corr = do_lmfit(ncorr, c1, sig, dt, nullptr, tStart, tendfit, oenv, bDebugMode(),
847 fitfn, fitparm, 0, nullptr);
848 sumtot = sum + tail_corr;
849 if (fit && ((jmax == 1) || (j == 1)))
850 {
851 double mfp[3];
852 for (i = 0; (i < 3); i++)
853 {
854 mfp[i] = fitparm[i];
855 }
856 for (i = 0; (i < ncorr); i++)
857 {
858 fit[i] = lmcurves[fitfn](i * dt, mfp);
859 }
860 }
861 if (bPrint)
862 {
863 printf("COR:%11.4e%11.4e%11.4e%11.4e", tStart, sum, tail_corr, sumtot);
864 for (i = 0; (i < effnNparams(fitfn)); i++)
865 {
866 printf(" %11.4e", fitparm[i]);
867 }
868 printf("\n");
869 }
870 tStart += tbeginfit;
871 }
872 sfree(sig);
873
874 return sumtot;
875 }
876