src/gromacs/correlationfunctions/expfit.cpp

   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,2014,2015,2016,2017,2018, 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 <string.h>
  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] = {
  73     nullptr, "none", "exp", "aexp", "exp_exp",
  74     "exp5", "exp7", "exp9",
  75     nullptr, nullptr, nullptr, nullptr, nullptr
  76 };
  77
  78 /*! \brief Long description for each fitting function type */
  79 static const char *longs_ffn[effnNR] = {
  80     "no fit",
  81     "y = exp(-x/|a0|)",
  82     "y = a1 exp(-x/|a0|)",
  83     "y = a1 exp(-x/|a0|) + (1-a1) exp(-x/(|a2|)), a2 > a0 > 0",
  84     "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4, a3 >= a1",
  85     "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6, a5 >= a3 >= a1",
  86     "y = a0 exp(-x/|a1|) +  a2 exp(-x/|a3|) + a4 exp(-x/|a5|) + a6 exp(-x/|a7|) + a8, a7 >= a5 >= a3 >= a1",
  87     "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a0), w = sqrt(1 - a1)",
  88     "y = 1/2*(a0+a1) - 1/2*(a0-a1)*erf( (x-a2) /a3^2)",
  89     "y = a1 *2*a0*((exp(-x/a0)-1)*(a0/x)+1)+(1-a1)*2*a2*((exp(-x/a2)-1)*(a2/x)+1)",
  90     "y = (1-a0)*cos(x*a1)*exp(-(x/a2)^a3) + a0*exp(-(x/a4)^a5)"
  91 };
  92
  93 int effnNparams(int effn)
  94 {
  95     if ((0 <= effn) && (effn < effnNR))
  96     {
  97         return nfp_ffn[effn];
  98     }
  99     else
 100     {
 101         return -1;
 102     }
 103 }
 104
 105 const char *effnDescription(int effn)
 106 {
 107     if ((0 <= effn) && (effn < effnNR))
 108     {
 109         return longs_ffn[effn];
 110     }
 111     else
 112     {
 113         return nullptr;
 114     }
 115 }
 116
 117 int sffn2effn(const char **sffn)
 118 {
 119     int eFitFn, i;
 120
 121     eFitFn = 0;
 122     for (i = 0; i < effnNR; i++)
 123     {
 124         if (sffn[i+1] && strcmp(sffn[0], sffn[i+1]) == 0)
 125         {
 126             eFitFn = i;
 127         }
 128     }
 129
 130     return eFitFn;
 131 }
 132
 133 /*! \brief Compute exponential function A exp(-x/tau) */
 134 static double myexp(double x, double A, double tau)
 135 {
 136     if ((A == 0) || (tau == 0))
 137     {
 138         return 0;
 139     }
 140     return A*exp(-x/tau);
 141 }
 142
 143 /*! \brief Compute y=(a0+a1)/2-(a0-a1)/2*erf((x-a2)/a3^2) */
 144 static double lmc_erffit (double x, const double *a)
 145 {
 146     double erfarg;
 147     double myerf;
 148
 149     if (a[3] != 0)
 150     {
 151         erfarg = (x-a[2])/(a[3]*a[3]);
 152         myerf  = std::erf(erfarg);
 153     }
 154     else
 155     {
 156         /* If a[3] == 0, a[3]^2 = 0 and the erfarg becomes +/- infinity */
 157         if (x < a[2])
 158         {
 159             myerf = -1;
 160         }
 161         else
 162         {
 163             myerf = 1;
 164         }
 165     }
 166     return 0.5*((a[0]+a[1]) - (a[0]-a[1])*myerf);
 167 }
 168
 169 /*! \brief Exponent function that prevents overflow */
 170 static double safe_exp(double x)
 171 {
 172     double exp_max = 200;
 173     double exp_min = -exp_max;
 174     if (x <= exp_min)
 175     {
 176         return exp(exp_min);
 177     }
 178     else if (x >= exp_max)
 179     {
 180         return exp(exp_max);
 181     }
 182     else
 183     {
 184         return exp(x);
 185     }
 186 }
 187
 188 /*! \brief Exponent minus 1 function that prevents overflow */
 189 static double safe_expm1(double x)
 190 {
 191     double exp_max = 200;
 192     double exp_min = -exp_max;
 193     if (x <= exp_min)
 194     {
 195         return -1;
 196     }
 197     else if (x >= exp_max)
 198     {
 199         return exp(exp_max);
 200     }
 201     else
 202     {
 203         return std::expm1(x);
 204     }
 205 }
 206
 207 /*! \brief Compute y = exp(-x/|a0|) */
 208 static double lmc_exp_one_parm(double x, const double *a)
 209 {
 210     return safe_exp(-x/fabs(a[0]));
 211 }
 212
 213 /*! \brief Compute y = a1 exp(-x/|a0|) */
 214 static double lmc_exp_two_parm(double x, const double *a)
 215 {
 216     return a[1]*safe_exp(-x/fabs(a[0]));
 217 }
 218
 219 /*! \brief Compute y = a1 exp(-x/|a0|) + (1-a1) exp(-x/|a2|) */
 220 static double lmc_exp_exp(double x, const double *a)
 221 {
 222     double e1, e2;
 223
 224     e1      = safe_exp(-x/fabs(a[0]));
 225     e2      = safe_exp(-x/(fabs(a[0])+fabs(a[2])));
 226     return a[1]*e1 + (1-a[1])*e2;
 227 }
 228
 229 /*! \brief Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) + a4 */
 230 static double lmc_exp_5_parm(double x, const double *a)
 231 {
 232     double e1, e2;
 233
 234     e1      = safe_exp(-x/fabs(a[1]));
 235     e2      = safe_exp(-x/(fabs(a[1])+fabs(a[3])));
 236     return a[0]*e1 + a[2]*e2 + a[4];
 237 }
 238
 239 /*! \brief Compute 7 parameter exponential function value.
 240  *
 241  * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
 242  * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6
 243  */
 244 static double lmc_exp_7_parm(double x, const double *a)
 245 {
 246     double e1, e2, e3;
 247     double fa1, fa3, fa5;
 248
 249     fa1 = fabs(a[1]);
 250     fa3 = fa1 + fabs(a[3]);
 251     fa5 = fa3 + fabs(a[5]);
 252     e1  = safe_exp(-x/fa1);
 253     e2  = safe_exp(-x/fa3);
 254     e3  = safe_exp(-x/fa5);
 255     return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6];
 256 }
 257
 258 /*! \brief Compute 9 parameter exponential function value.
 259  *
 260  * Compute y = a0 exp(-x/|a1|) + a2 exp(-x/(|a1|+|a3|)) +
 261  * a4 exp(-x/(|a1|+|a3|+|a5|)) + a6 exp(-x/(|a1|+|a3|+|a5|+|a7|)) + a8
 262  */
 263 static double lmc_exp_9_parm(double x, const double *a)
 264 {
 265     double e1, e2, e3, e4;
 266     double fa1, fa3, fa5, fa7;
 267
 268     fa1 = fabs(a[1]);
 269     fa3 = fa1 + fabs(a[3]);
 270     fa5 = fa3 + fabs(a[5]);
 271     fa7 = fa5 + fabs(a[7]);
 272
 273     e1      = safe_exp(-x/fa1);
 274     e2      = safe_exp(-x/fa3);
 275     e3      = safe_exp(-x/fa5);
 276     e4      = safe_exp(-x/fa7);
 277     return a[0]*e1 + a[2]*e2 + a[4]*e3 + a[6]*e4 + a[8];
 278 }
 279
 280 /*! \brief Compute y = (1-a0)*exp(-(x/|a2|)^|a3|)*cos(x*|a1|) + a0*exp(-(x/|a4|)^|a5|) */
 281 static double lmc_pres_6_parm(double x, const double *a)
 282 {
 283     double term1, term2, term3;
 284     double pow_max = 10;
 285
 286     term3  = 0;
 287     if ((a[4] != 0) && (a[0] != 0))
 288     {
 289         double power = std::min(fabs(a[5]), pow_max);
 290         term3        = a[0] * safe_exp(-pow((x/fabs(a[4])), power));
 291     }
 292
 293     term1  = 1-a[0];
 294     term2  = 0;
 295     if ((term1 != 0) && (a[2] != 0))
 296     {
 297         double power = std::min(fabs(a[3]), pow_max);
 298         term2        = safe_exp(-pow((x/fabs(a[2])), power)) * cos(x*fabs(a[1]));
 299     }
 300
 301     return term1*term2 + term3;
 302 }
 303
 304 /*! \brief Compute vac function */
 305 static double lmc_vac_2_parm(double x, const double *a)
 306 {
 307     /* Fit to function
 308      *
 309      * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 310      *
 311      *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 312      *
 313      *    v = x/(2 a0)
 314      *    w = sqrt(1 - a1)
 315      *
 316      *    For tranverse current autocorrelation functions:
 317      *       a0 = tau
 318      *       a1 = 4 tau (eta/rho) k^2
 319      *
 320      */
 321
 322     double y, v, det, omega, wv, em, ec, es;
 323     double wv_max = 100;
 324
 325     v   = x/(2*fabs(a[0]));
 326     det = 1 - a[1];
 327     em  = safe_exp(-v);
 328     if (det != 0)
 329     {
 330         omega = sqrt(fabs(det));
 331         wv    = std::min(omega*v, wv_max);
 332
 333         if (det > 0)
 334         {
 335             ec = em*0.5*(safe_exp(wv)+safe_exp(-wv));
 336             es = em*0.5*(safe_exp(wv)-safe_exp(-wv))/omega;
 337         }
 338         else
 339         {
 340             ec = em*cos(wv);
 341             es = em*sin(wv)/omega;
 342         }
 343         y      = ec + es;
 344     }
 345     else
 346     {
 347         y      = (1+v)*em;
 348     }
 349     return y;
 350 }
 351
 352 /*! \brief Compute error estimate */
 353 static double lmc_errest_3_parm(double x, const double *a)
 354 {
 355     double e1, e2, v1, v2;
 356     double fa0 = fabs(a[0]);
 357     double fa1;
 358     double fa2 = fa0+fabs(a[2]);
 359
 360     if (a[0] != 0)
 361     {
 362         e1 = safe_expm1(-x/fa0);
 363     }
 364     else
 365     {
 366         e1 = 0;
 367     }
 368     if (a[2] != 0)
 369     {
 370         e2 = safe_expm1(-x/fa2);
 371     }
 372     else
 373     {
 374         e2 = 0;
 375     }
 376
 377     if (x > 0)
 378     {
 379         v1      = 2*fa0*(e1*fa0/x + 1);
 380         v2      = 2*fa2*(e2*fa2/x + 1);
 381         /* We need 0 <= a1 <= 1 */
 382         fa1     = std::min(1.0, std::max(0.0, a[1]));
 383
 384         return fa1*v1 + (1-fa1)*v2;
 385     }
 386     else
 387     {
 388         return 0;
 389     }
 390 }
 391
 392 /*! \brief array of fitting functions corresponding to the pre-defined types */
 393 t_lmcurve lmcurves[effnNR+1] = {
 394     lmc_exp_one_parm, lmc_exp_one_parm, lmc_exp_two_parm,
 395     lmc_exp_exp,      lmc_exp_5_parm,   lmc_exp_7_parm,
 396     lmc_exp_9_parm,
 397     lmc_vac_2_parm,   lmc_erffit,  lmc_errest_3_parm, lmc_pres_6_parm
 398 };
 399
 400 double fit_function(const int eFitFn, const double parm[], const double x)
 401 {
 402     if ((eFitFn < 0) || (eFitFn >= effnNR))
 403     {
 404         fprintf(stderr, "fitfn = %d, should be in the range 0..%d\n",
 405                 eFitFn, effnNR-1);
 406         return 0.0;
 407     }
 408     return lmcurves[eFitFn](x, parm);
 409 }
 410
 411 /*! \brief Ensure the fitting parameters are well-behaved.
 412  *
 413  * In order to make sure that fitting behaves according to the
 414  * constraint that time constants are positive and increasing
 415  * we enforce this here by setting all time constants to their
 416  * absolute value and by adding e.g. |a_0| to |a_2|. This is
 417  * done in conjunction with the fitting functions themselves.
 418  * When there are multiple time constants we make sure that
 419  * the successive time constants are at least double the
 420  * previous ones and during fitting we enforce the they remain larger.
 421  * It may very well help the convergence of the fitting routine.
 422  */
 423 static void initiate_fit_params(int    eFitFn,
 424                                 double params[])
 425 {
 426     int i, nparm;
 427
 428     nparm = effnNparams(eFitFn);
 429
 430     switch (eFitFn)
 431     {
 432         case effnVAC:
 433             GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
 434             break;
 435         case effnEXP1:
 436         case effnEXP2:
 437         case effnEXPEXP:
 438             GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
 439             if (nparm > 2)
 440             {
 441                 GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
 442                 /* In order to maintain params[2] >= params[0] in the final
 443                  * result, we fit the difference between the two, that is
 444                  * params[2]-params[0] and in the add add in params[0]
 445                  * again.
 446                  */
 447                 params[2] = std::max(fabs(params[2])-params[0], params[0]);
 448             }
 449             break;
 450         case effnEXP5:
 451         case effnEXP7:
 452         case effnEXP9:
 453             GMX_ASSERT(params[1] >= 0, "parameters should be >= 0");
 454             params[1] = fabs(params[1]);
 455             if (nparm > 3)
 456             {
 457                 GMX_ASSERT(params[3] >= 0, "parameters should be >= 0");
 458                 /* See comment under effnEXPEXP */
 459                 params[3] = std::max(fabs(params[3])-params[1], params[1]);
 460                 if (nparm > 5)
 461                 {
 462                     GMX_ASSERT(params[5] >= 0, "parameters should be >= 0");
 463                     /* See comment under effnEXPEXP */
 464                     params[5] = std::max(fabs(params[5])-params[3], params[3]);
 465                     if (nparm > 7)
 466                     {
 467                         GMX_ASSERT(params[7] >= 0, "parameters should be >= 0");
 468                         /* See comment under effnEXPEXP */
 469                         params[7] = std::max(fabs(params[7])-params[5], params[5]);
 470                     }
 471                 }
 472             }
 473             break;
 474         case effnERREST:
 475             GMX_ASSERT(params[0] >= 0, "parameters should be >= 0");
 476             GMX_ASSERT(params[1] >= 0 && params[1] <= 1, "parameter 1 should in 0 .. 1");
 477             GMX_ASSERT(params[2] >= 0, "parameters should be >= 0");
 478             /* See comment under effnEXPEXP */
 479             params[2] = fabs(params[2])-params[0];
 480             break;
 481         case effnPRES:
 482             for (i = 1; (i < nparm); i++)
 483             {
 484                 GMX_ASSERT(params[i] >= 0, "parameters should be >= 0");
 485             }
 486             break;
 487         default:
 488             break;
 489     }
 490 }
 491
 492 /*! \brief Process the fitting parameters to get output parameters.
 493  *
 494  * See comment at the previous function.
 495  */
 496 static void extract_fit_params(int    eFitFn,
 497                                double params[])
 498 {
 499     int i, nparm;
 500
 501     nparm = effnNparams(eFitFn);
 502
 503     switch (eFitFn)
 504     {
 505         case effnVAC:
 506             params[0] = fabs(params[0]);
 507             break;
 508         case effnEXP1:
 509         case effnEXP2:
 510         case effnEXPEXP:
 511             params[0] = fabs(params[0]);
 512             if (nparm > 2)
 513             {
 514                 /* Back conversion of parameters from the fitted difference
 515                  * to the absolute value.
 516                  */
 517                 params[2] = fabs(params[2])+params[0];
 518             }
 519             break;
 520         case effnEXP5:
 521         case effnEXP7:
 522         case effnEXP9:
 523             params[1] = fabs(params[1]);
 524             if (nparm > 3)
 525             {
 526                 /* See comment under effnEXPEXP */
 527                 params[3] = fabs(params[3])+params[1];
 528                 if (nparm > 5)
 529                 {
 530                     /* See comment under effnEXPEXP */
 531                     params[5] = fabs(params[5])+params[3];
 532                     if (nparm > 7)
 533                     {
 534                         /* See comment under effnEXPEXP */
 535                         params[7] = fabs(params[7])+params[5];
 536                     }
 537                 }
 538             }
 539             break;
 540         case effnERREST:
 541             params[0] = fabs(params[0]);
 542             if (params[1] < 0)
 543             {
 544                 params[1] = 0;
 545             }
 546             else if (params[1] > 1)
 547             {
 548                 params[1] = 1;
 549             }
 550             /* See comment under effnEXPEXP */
 551             params[2] = params[0]+fabs(params[2]);
 552             break;
 553         case effnPRES:
 554             for (i = 1; (i < nparm); i++)
 555             {
 556                 params[i] = fabs(params[i]);
 557             }
 558             break;
 559         default:
 560             break;
 561     }
 562 }
 563
 564 /*! \brief Print chi-squared value and the parameters */
 565 static void print_chi2_params(FILE        *fp,
 566                               const int    eFitFn,
 567                               const double fitparms[],
 568                               const char  *label,
 569                               const int    nfitpnts,
 570                               const double x[],
 571                               const double y[])
 572 {
 573     int    i;
 574     double chi2 = 0;
 575
 576     for (i = 0; (i < nfitpnts); i++)
 577     {
 578         double yfit = lmcurves[eFitFn](x[i], fitparms);
 579         chi2 += gmx::square(y[i] - yfit);
 580     }
 581     fprintf(fp, "There are %d data points, %d parameters, %s chi2 = %g\nparams:",
 582             nfitpnts, effnNparams(eFitFn), label, chi2);
 583     for (i = 0; (i < effnNparams(eFitFn)); i++)
 584     {
 585         fprintf(fp, "  %10g", fitparms[i]);
 586     }
 587     fprintf(fp, "\n");
 588 }
 589
 590 real do_lmfit(int ndata, const real c1[], real sig[], real dt, const real *x0,
 591               real begintimefit, real endtimefit, const gmx_output_env_t *oenv,
 592               gmx_bool bVerbose, int eFitFn, double fitparms[], int fix,
 593               const char *fn_fitted)
 594 {
 595     FILE    *fp;
 596     int      i, j, nfitpnts;
 597     double   integral, ttt;
 598     double  *x, *y, *dy;
 599
 600     if (0 != fix)
 601     {
 602         fprintf(stderr, "Using fixed parameters in curve fitting is temporarily not working.\n");
 603     }
 604     if (debug)
 605     {
 606         fprintf(debug, "There are %d points to fit %d vars!\n", ndata, effnNparams(eFitFn));
 607         fprintf(debug, "Fit to function %d from %g through %g, dt=%g\n",
 608                 eFitFn, begintimefit, endtimefit, dt);
 609     }
 610
 611     snew(x, ndata);
 612     snew(y, ndata);
 613     snew(dy, ndata);
 614     j    = 0;
 615     for (i = 0; i < ndata; i++)
 616     {
 617         ttt = x0 ? x0[i] : dt*i;
 618         if (ttt >= begintimefit && ttt <= endtimefit)
 619         {
 620             x[j]  = ttt;
 621             y[j]  = c1[i];
 622             if (nullptr == sig)
 623             {
 624                 // No weighting if all values are divided by 1.
 625                 dy[j] = 1;
 626             }
 627             else
 628             {
 629                 dy[j] = std::max(1.0e-7, (double)sig[i]);
 630             }
 631             if (debug)
 632             {
 633                 fprintf(debug, "j= %d, i= %d, x= %g, y= %g, dy=%g, ttt=%g\n",
 634                         j, i, x[j], y[j], dy[j], ttt);
 635             }
 636             j++;
 637         }
 638     }
 639     nfitpnts = j;
 640     integral = 0;
 641     if (nfitpnts < effnNparams(eFitFn))
 642     {
 643         fprintf(stderr, "Not enough (%d) data points for fitting, dt = %g!\n",
 644                 nfitpnts, dt);
 645     }
 646     else
 647     {
 648         gmx_bool bSuccess;
 649
 650         if (bVerbose)
 651         {
 652             print_chi2_params(stdout, eFitFn, fitparms, "initial", nfitpnts, x, y);
 653         }
 654         initiate_fit_params(eFitFn, fitparms);
 655
 656         bSuccess = lmfit_exp(nfitpnts, x, y, dy, fitparms, bVerbose, eFitFn, fix);
 657         extract_fit_params(eFitFn, fitparms);
 658
 659         if (!bSuccess)
 660         {
 661             fprintf(stderr, "Fit failed!\n");
 662         }
 663         else
 664         {
 665             if (bVerbose)
 666             {
 667                 print_chi2_params(stdout, eFitFn, fitparms, "final", nfitpnts, x, y);
 668             }
 669             switch (eFitFn)
 670             {
 671                 case effnEXP1:
 672                     integral = fitparms[0]*myexp(begintimefit, 1,  fitparms[0]);
 673                     break;
 674                 case effnEXP2:
 675                     integral = fitparms[0]*myexp(begintimefit, fitparms[1],  fitparms[0]);
 676                     break;
 677                 case effnEXPEXP:
 678                     integral = (fitparms[0]*myexp(begintimefit, fitparms[1],  fitparms[0]) +
 679                                 fitparms[2]*myexp(begintimefit, 1-fitparms[1], fitparms[2]));
 680                     break;
 681                 case effnEXP5:
 682                 case effnEXP7:
 683                 case effnEXP9:
 684                     integral = 0;
 685                     for (i = 0; (i < (effnNparams(eFitFn)-1)/2); i++)
 686                     {
 687                         integral += fitparms[2*i]*myexp(begintimefit, fitparms[2*i+1], fitparms[2*i]);
 688                     }
 689                     break;
 690                 default:
 691                     /* Do numerical integration */
 692                     integral = 0;
 693                     for (i = 0; (i < nfitpnts-1); i++)
 694                     {
 695                         double y0 = lmcurves[eFitFn](x[i], fitparms);
 696                         double y1 = lmcurves[eFitFn](x[i+1], fitparms);
 697                         integral += (x[i+1]-x[i])*(y1+y0)*0.5;
 698                     }
 699                     break;
 700             }
 701
 702             if (bVerbose)
 703             {
 704                 printf("FIT: Integral of fitted function: %g\n", integral);
 705                 if ((effnEXP5 == eFitFn) || (effnEXP7 == eFitFn) || (effnEXP9 == eFitFn))
 706                 {
 707                     printf("FIT: Note that the constant term is not taken into account when computing integral.\n");
 708                 }
 709             }
 710             /* Generate debug output */
 711             if (nullptr != fn_fitted)
 712             {
 713                 fp = xvgropen(fn_fitted, "Data + Fit", "Time (ps)",
 714                               "Data (t)", oenv);
 715                 for (i = 0; (i < effnNparams(eFitFn)); i++)
 716                 {
 717                     fprintf(fp, "# fitparms[%d] = %g\n", i, fitparms[i]);
 718                 }
 719                 for (j = 0; (j < nfitpnts); j++)
 720                 {
 721                     real ttt = x0 ? x0[j] : dt*j;
 722                     fprintf(fp, "%10.5e  %10.5e  %10.5e\n",
 723                             x[j], y[j], (lmcurves[eFitFn])(ttt, fitparms));
 724                 }
 725                 xvgrclose(fp);
 726             }
 727         }
 728     }
 729
 730     sfree(x);
 731     sfree(y);
 732     sfree(dy);
 733
 734     return integral;
 735 }
 736
 737 real fit_acf(int ncorr, int fitfn, const gmx_output_env_t *oenv, gmx_bool bVerbose,
 738              real tbeginfit, real tendfit, real dt, real c1[], 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, (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",
 765                0.0, 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", tbeginfit, std::min(ncorr*dt, tendfit));
 769     }
 770
 771     tStart = 0;
 772     if (bPrint)
 773     {
 774         printf("COR:%11s%11s%11s%11s%11s%11s%11s\n",
 775                "Fit from", "Integral", "Tail Value", "Sum (ps)", " a1 (ps)",
 776                (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, (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,
 847                              bDebugMode(), 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 }