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