src/tools/expfit.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Green Red Orange Magenta Azure Cyan Skyblue
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39
  40 #include <sysstuff.h>
  41 #include <string.h>
  42 #include <math.h>
  43 #include "typedefs.h"
  44 #include "smalloc.h"
  45 #include "xvgr.h"
  46 #include "futil.h"
  47 #include "gstat.h"
  48 #include "vec.h"
  49 #include "statutil.h"
  50 #include "index.h"
  51
  52 int  nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 3 };
  53
  54 const char *s_ffn[effnNR+2] = { NULL, "none", "exp", "aexp", "exp_exp", "vac",
  55                                 "exp5", "exp7", "exp9", NULL, NULL };
  56 /* We don't allow errest as a choice on the command line */
  57
  58 const char *longs_ffn[effnNR] = {
  59   "no fit",
  60   "y = exp(-x/a1)",
  61   "y = a2 exp(-x/a1)",
  62   "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
  63   "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
  64   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5",
  65   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
  66   "y = a1 exp(-x/a2) +  a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
  67   "y = a2*ee(a1,x) + (1-a2)*ee(a2,x)"
  68 };
  69
  70 extern bool mrqmin(real x[],real y[],real sig[],int ndata,real a[],
  71                    int ma,int lista[],int mfit,real **covar,real **alpha,
  72                    real *chisq,
  73                    void (*funcs)(real x,real a[],real *y,real dyda[]),
  74                    real *alamda);
  75
  76 extern bool mrqmin_new(real x[],real y[],real sig[],int ndata,real a[],
  77                        int ia[],int ma,real **covar,real **alpha,real *chisq,
  78                        void (*funcs)(real, real [], real *, real []),
  79                        real *alamda);
  80
  81 static real myexp(real x,real A,real tau)
  82 {
  83   if ((A == 0) || (tau == 0))
  84     return 0;
  85   return A*exp(-x/tau);
  86 }
  87
  88 static void exp_one_parm(real x,real a[],real *y,real dyda[])
  89 {
  90   /* Fit to function
  91    *
  92    * y = exp(-x/a1)
  93    *
  94    */
  95
  96   real e1;
  97
  98   e1      = exp(-x/a[1]);
  99   *y      = e1;
 100   dyda[1] = x*e1/sqr(a[1]);
 101 }
 102
 103 static void exp_two_parm(real x,real a[],real *y,real dyda[])
 104 {
 105   /* Fit to function
 106    *
 107    * y = a2 exp(-x/a1)
 108    *
 109    */
 110
 111   real e1;
 112
 113   e1      = exp(-x/a[1]);
 114   *y      = a[2]*e1;
 115   dyda[1] = x*a[2]*e1/sqr(a[1]);
 116   dyda[2] = e1;
 117 }
 118
 119 static void exp_3_parm(real x,real a[],real *y,real dyda[])
 120 {
 121   /* Fit to function
 122    *
 123    * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
 124    *
 125    */
 126
 127   real e1,e2;
 128
 129   e1      = exp(-x/a[1]);
 130   e2      = exp(-x/a[3]);
 131   *y      = a[2]*e1 + (1-a[2])*e2;
 132   dyda[1] = x*a[2]*e1/sqr(a[1]);
 133   dyda[2] = e1-e2;
 134   dyda[3] = x*(1-a[2])*e2/sqr(a[3]);
 135 }
 136
 137 static void exp_5_parm(real x,real a[],real *y,real dyda[])
 138 {
 139   /* Fit to function
 140    *
 141    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5
 142    *
 143    */
 144
 145   real e1,e2;
 146
 147   e1      = exp(-x/a[2]);
 148   e2      = exp(-x/a[4]);
 149   *y      = a[1]*e1 + a[3]*e2 + a[5];
 150
 151   if (debug)
 152     fprintf(debug,"exp_5_parm called: x = %10.3f  y = %10.3f\n"
 153             "a = ( %8.3f  %8.3f  %8.3f  %8.3f  %8.3f)\n",
 154             x,*y,a[1],a[2],a[3],a[4],a[5]);
 155   dyda[1] = e1;
 156   dyda[2] = x*e1/sqr(a[2]);
 157   dyda[3] = e2;
 158   dyda[4] = x*e2/sqr(a[4]);
 159   dyda[5] = 0;
 160 }
 161
 162 static void exp_7_parm(real x,real a[],real *y,real dyda[])
 163 {
 164   /* Fit to function
 165    *
 166    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 167    *
 168    */
 169
 170   real e1,e2,e3;
 171
 172   e1      = exp(-x/a[2]);
 173   e2      = exp(-x/a[4]);
 174   e3      = exp(-x/a[6]);
 175   *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7];
 176
 177   dyda[1] = e1;
 178   dyda[2] = x*e1/sqr(a[2]);
 179   dyda[3] = e2;
 180   dyda[4] = x*e2/sqr(a[4]);
 181   dyda[5] = e3;
 182   dyda[6] = x*e3/sqr(a[6]);
 183   dyda[7] = 0;
 184 }
 185
 186 static void exp_9_parm(real x,real a[],real *y,real dyda[])
 187 {
 188   /* Fit to function
 189    *
 190    * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
 191    *
 192    */
 193
 194   real e1,e2,e3,e4;
 195
 196   e1      = exp(-x/a[2]);
 197   e2      = exp(-x/a[4]);
 198   e3      = exp(-x/a[6]);
 199   e4      = exp(-x/a[8]);
 200   *y      = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7]*e4 + a[9];
 201
 202   dyda[1] = e1;
 203   dyda[2] = x*e1/sqr(a[2]);
 204   dyda[3] = e2;
 205   dyda[4] = x*e2/sqr(a[4]);
 206   dyda[5] = e3;
 207   dyda[6] = x*e3/sqr(a[6]);
 208   dyda[7] = e4;
 209   dyda[8] = x*e4/sqr(a[8]);
 210   dyda[9] = 0;
 211 }
 212
 213 static void vac_2_parm(real x,real a[],real *y,real dyda[])
 214 {
 215   /* Fit to function
 216    *
 217    * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
 218    *
 219    *   = exp(-v) (cosh(wv) + 1/w sinh(wv))
 220    *
 221    *    v = x/(2 a1)
 222    *    w = sqrt(1 - a2)
 223    *
 224    *    For tranverse current autocorrelation functions:
 225    *       a1 = tau
 226    *       a2 = 4 tau (eta/rho) k^2
 227    *
 228    */
 229
 230   double v,det,omega,omega2,em,ec,es;
 231
 232   v   = x/(2*a[1]);
 233   det = 1 - a[2];
 234   em  = exp(-v);
 235   if (det != 0) {
 236     omega2  = fabs(det);
 237     omega   = sqrt(omega2);
 238     if (det > 0) {
 239       ec = em*0.5*(exp(omega*v)+exp(-omega*v));
 240       es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
 241     } else {
 242       ec = em*cos(omega*v);
 243       es = em*sin(omega*v)/omega;
 244     }
 245     *y      = ec + es;
 246     dyda[2] = (v/det*ec+(v-1/det)*es)/(-2.0);
 247     dyda[1] = (1-det)*v/a[1]*es;
 248   } else {
 249     *y      = (1+v)*em;
 250     dyda[2] = -v*v*em*(0.5+v/6);
 251     dyda[1] = v*v/a[1]*em;
 252   }
 253 }
 254
 255 static void errest_3_parm(real x,real a[],real *y,real dyda[])
 256 {
 257   real e1,e2,v1,v2;
 258
 259   if (a[1])
 260     e1 = exp(-x/a[1]) - 1;
 261   else
 262     e1 = 0;
 263   if (a[3])
 264     e2 = exp(-x/a[3]) - 1;
 265   else
 266     e2 = 0;
 267
 268   if (x > 0) {
 269     v1 = 2*a[1]*(e1*a[1]/x + 1);
 270     v2 = 2*a[3]*(e2*a[3]/x + 1);
 271     *y      = a[2]*v1 + (1-a[2])*v2;
 272     dyda[1] = 2*     a[2] *(v1/a[1] + e1);
 273     dyda[3] = 2*(1 - a[2])*(v2/a[3] + e2);
 274     dyda[2] = (v1 - v2);
 275   } else {
 276     *y      = 0;
 277     dyda[1] = 0;
 278     dyda[3] = 0;
 279     dyda[2] = 0;
 280   }
 281 }
 282
 283 typedef void (*myfitfn)(real x,real a[],real *y,real dyda[]);
 284 myfitfn mfitfn[effnNR] = {
 285   exp_one_parm, exp_one_parm, exp_two_parm, exp_3_parm, vac_2_parm,
 286   exp_5_parm,   exp_7_parm,   exp_9_parm,   errest_3_parm
 287 };
 288
 289 real fit_function(int eFitFn,real *parm,real x)
 290 {
 291   static real y,dum[8];
 292
 293   mfitfn[eFitFn](x,parm-1,&y,dum);
 294
 295   return y;
 296 }
 297
 298 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
 299 static bool lmfit_exp(int nfit,real x[],real y[],real dy[],real ftol,
 300                       real parm[],real dparm[],bool bVerbose,
 301                       int eFitFn,int fix)
 302 {
 303   real chisq,ochisq,alamda;
 304   real *a,**covar,**alpha,*dum;
 305   bool bCont;
 306   int  i,j,ma,mfit,*lista,*ia;
 307
 308   if ((eFitFn < 0) || (eFitFn >= effnNR))
 309     gmx_fatal(FARGS,"fitfn = %d, should be in 0..%d (%s,%d)",
 310                 effnNR-1,eFitFn,__FILE__,__LINE__);
 311
 312   ma=mfit=nfp_ffn[eFitFn];         /* number of parameters to fit */
 313   snew(a,ma+1);
 314   snew(covar,ma+1);
 315   snew(alpha,ma+1);
 316   snew(lista,ma+1);
 317   snew(ia,ma+1);
 318   snew(dum,ma+1);
 319   for(i=1; (i<ma+1); i++) {
 320     lista[i] = i;
 321     ia[i] = i;
 322     snew(covar[i],ma+1);
 323     snew(alpha[i],ma+1);
 324   }
 325   if (fix) {
 326     if (bVerbose)
 327       fprintf(stderr,"Will keep parameters fixed during fit procedure: %d\n",
 328               fix);
 329     for(i=0; i<ma; i++)
 330       if (fix & 1<<i)
 331         ia[i+1] = 0;
 332   }
 333   if (debug)
 334     fprintf(debug,"%d parameter fit\n",mfit);
 335
 336   /* Initial params */
 337   alamda = -1;    /* Starting value   */
 338   chisq  = 1e12;
 339   for(i=0; (i<mfit); i++)
 340     a[i+1] = parm[i];
 341
 342   j = 0;
 343   if (bVerbose)
 344     fprintf(stderr,"%4s  %10s  %10s  %10s  %10s  %10s\n",
 345             "Step","chi^2","Lambda","A1","A2","A3");
 346   do {
 347     ochisq = chisq;
 348     /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 349      *   &chisq,expfn[mfit-1],&alamda)
 350      */
 351     if (!mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
 352                     mfitfn[eFitFn],&alamda))
 353       return FALSE;
 354
 355     if (bVerbose) {
 356       fprintf(stderr,"%4d  %10.5e  %10.5e  %10.5e",
 357               j,chisq,alamda,a[1]);
 358       if (mfit > 1)
 359         fprintf(stderr,"  %10.5e",a[2]);
 360       if (mfit > 2)
 361         fprintf(stderr,"  %10.5e",a[3]);
 362       fprintf(stderr,"\n");
 363     }
 364     j++;
 365     bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
 366              ((ochisq == chisq)));
 367   } while (bCont && (alamda != 0.0) && (j < 50));
 368   if (bVerbose)
 369     fprintf(stderr,"\n");
 370
 371   /* Now get the covariance matrix out */
 372   alamda = 0;
 373
 374   /*  mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
 375    * &chisq,expfn[mfit-1],&alamda)
 376    */
 377   if ( !mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
 378                    mfitfn[eFitFn],&alamda))
 379     return FALSE;
 380
 381   for(j=0; (j<mfit); j++) {
 382     parm[j]  = a[j+1];
 383     dparm[j] = covar[j+1][j+1];
 384   }
 385
 386   for(i=0; (i<ma+1); i++) {
 387     sfree(covar[i]);
 388     sfree(alpha[i]);
 389   }
 390   sfree(a);
 391   sfree(covar);
 392   sfree(alpha);
 393   sfree(lista);
 394   sfree(dum);
 395
 396   return TRUE;
 397 }
 398
 399 real do_lmfit(int ndata,real c1[],real sig[],real dt,real x0[],
 400               real begintimefit,real endtimefit,const output_env_t oenv,
 401               bool bVerbose, int eFitFn,real fitparms[],int fix)
 402 {
 403   FILE *fp;
 404   char buf[32];
 405
 406   int  i,j,nparm,nfitpnts;
 407   real integral,ttt;
 408   real *parm,*dparm;
 409   real *x,*y,*dy;
 410   real ftol = 1e-4;
 411
 412   nparm = nfp_ffn[eFitFn];
 413   if (debug) {
 414     fprintf(debug,"There are %d points to fit %d vars!\n",ndata,nparm);
 415     fprintf(debug,"Fit to function %d from %g thru %g, dt=%g\n",
 416             eFitFn,begintimefit,endtimefit,dt);
 417   }
 418
 419   snew(x,ndata);
 420   snew(y,ndata);
 421   snew(dy,ndata);
 422
 423   j=0;
 424   for(i=0; i<ndata; i++) {
 425     ttt = x0 ? x0[i] : dt*i;
 426     if (ttt>=begintimefit && ttt<=endtimefit) {
 427       x[j] = ttt;
 428       y[j] = c1[i];
 429
 430       /* mrqmin does not like sig to be zero */
 431       if (sig[i]<1.0e-7)
 432         dy[j]=1.0e-7;
 433       else
 434         dy[j]=sig[i];
 435       if (debug)
 436         fprintf(debug,"j= %d, i= %d, x= %g, y= %g, dy= %g\n",
 437                 j,i,x[j],y[j],dy[j]);
 438       j++;
 439     }
 440   }
 441   nfitpnts = j;
 442   integral = 0;
 443   if (nfitpnts < nparm)
 444     fprintf(stderr,"Not enough data points for fitting!\n");
 445   else {
 446     snew(parm,nparm);
 447     snew(dparm,nparm);
 448     if (fitparms)
 449       for(i=0; (i < nparm); i++)
 450         parm[i]=fitparms[i];
 451
 452     if (!lmfit_exp(nfitpnts,x,y,dy,ftol,parm,dparm,bVerbose,eFitFn,fix))
 453       fprintf(stderr,"Fit failed!\n");
 454     else if (nparm <= 3) {
 455       /* Compute the integral from begintimefit */
 456       if (nparm == 3)
 457         integral=(parm[0]*myexp(begintimefit,parm[1],  parm[0]) +
 458                   parm[2]*myexp(begintimefit,1-parm[1],parm[2]));
 459       else if (nparm == 2)
 460         integral=parm[0]*myexp(begintimefit,parm[1],  parm[0]);
 461       else if (nparm == 1)
 462         integral=parm[0]*myexp(begintimefit,1,  parm[0]);
 463       else
 464         gmx_fatal(FARGS,"nparm = %d in file %s, line %d",
 465                     nparm,__FILE__,__LINE__);
 466
 467       /* Generate THE output */
 468       if (bVerbose) {
 469         fprintf(stderr,"FIT: # points used in fit is: %d\n",nfitpnts);
 470         fprintf(stderr,"FIT: %21s%21s%21s\n",
 471                 "parm0     ","parm1 (ps)   ","parm2 (ps)    ");
 472         fprintf(stderr,"FIT: ------------------------------------------------------------\n");
 473         fprintf(stderr,"FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
 474                 parm[0],dparm[0],parm[1],dparm[1],parm[2],dparm[2]);
 475         fprintf(stderr,"FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
 476                 begintimefit,integral);
 477
 478         sprintf(buf,"test%d.xvg",nfitpnts);
 479         fp = xvgropen(buf,"C(t) + Fit to C(t)","Time (ps)","C(t)",oenv);
 480         fprintf(fp,"# parm0 = %g, parm1 = %g, parm2 = %g\n",
 481                 parm[0],parm[1],parm[2]);
 482         for(j=0; (j<nfitpnts); j++) {
 483           ttt = x0 ? x0[j] : dt*j;
 484           fprintf(fp,"%10.5e  %10.5e  %10.5e\n",
 485                   ttt,c1[j],fit_function(eFitFn,parm,ttt));
 486         }
 487         ffclose(fp);
 488       }
 489     }
 490     for(i=0;(i<nparm);i++)
 491       fitparms[i] = parm[i];
 492     sfree(parm);
 493     sfree(dparm);
 494   }
 495
 496   sfree(x);
 497   sfree(y);
 498   sfree(dy);
 499
 500   return integral;
 501 }
 502
 503 void do_expfit(int ndata,real c1[],real dt,real begintimefit,real endtimefit)
 504 {
 505   int i,n;
 506   real *x,*y,*Dy;
 507   real aa,bb,saa,sbb,A,tau,dA,dtau;
 508
 509   fprintf(stderr,"Will fit data from %g (ps) to %g (ps).\n",
 510           begintimefit,endtimefit);
 511
 512   snew(x,ndata);   /* allocate the maximum necessary space */
 513   snew(y,ndata);
 514   snew(Dy,ndata);
 515   n=0;
 516
 517   for(i=0; (i<ndata); i++) {
 518     if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) ) {
 519       x[n]=dt*i;
 520       y[n]=c1[i];
 521       Dy[n]=0.5;
 522       fprintf(stderr,"n= %d, i= %d, x= %g, y= %g\n",n,i,x[n],y[n]);
 523       n++;
 524     }
 525   }
 526   fprintf(stderr,"# of data points used in the fit is : %d\n\n",n);
 527   expfit(n,x,y,Dy,&aa,&saa,&bb,&sbb);
 528
 529   A=exp(aa);
 530   dA=exp(aa)*saa;
 531   tau=-1.0/bb;
 532   dtau=sbb/sqr(bb);
 533   fprintf(stderr,"Fitted to y=exp(a+bx):\n");
 534   fprintf(stderr,"a = %10.5f\t b = %10.5f",aa,bb);
 535   fprintf(stderr,"\n");
 536   fprintf(stderr,"Fitted to y=Aexp(-x/tau):\n");
 537   fprintf(stderr,"A  = %10.5f\t tau  = %10.5f\n",A,tau);
 538   fprintf(stderr,"dA = %10.5f\t dtau = %10.5f\n",dA,dtau);
 539 }
 540
 541
 542 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
 543             real *b, real *sb)
 544 {
 545   real *w,*ly,A,SA,B,SB;
 546   int  i;
 547   real sum,xbar,ybar,Sxx,Sxy,wr2,chi2,gamma,Db;
 548
 549 #define ZERO 0.0
 550 #define ONE 1.0
 551 #define ONEP5 1.5
 552 #define TWO 2.0
 553
 554 #define sqr(x) ((x)*(x))
 555
 556   /*allocate memory */
 557   snew(w,n);
 558   snew(ly,n);
 559
 560   /* Calculate weights and values of ln(y). */
 561   for(i=0;(i<n); i++){
 562     w[i]=sqr(y[i]/Dy[i]);
 563     ly[i]=log(y[i]);
 564   }
 565
 566   /* The weighted averages of x and y: xbar and ybar. */
 567   sum=ZERO;
 568   xbar=ZERO;
 569   ybar=ZERO;
 570   for(i=0;(i<n);i++){
 571     sum+=w[i];
 572     xbar+=w[i]*x[i];
 573     ybar+=w[i]*ly[i];
 574   }
 575   xbar/=sum;
 576   ybar/=sum;
 577
 578   /* The centered product sums Sxx and Sxy, and hence A and B. */
 579   Sxx=ZERO;
 580   Sxy=ZERO;
 581   for(i=0;(i<n);i++){
 582     Sxx+=w[i]*sqr(x[i]-xbar);
 583     Sxy+=w[i]*(x[i]-xbar)*(ly[i]-ybar);
 584   }
 585   B=Sxy/Sxx;
 586   A=ybar-B*xbar;
 587
 588   /* Chi-squared (chi2) and gamma. */
 589   chi2=ZERO;
 590   gamma=ZERO;
 591   for(i=0;(i<n);i++){
 592     wr2=w[i]*sqr(ly[i]-A-B*x[i]);
 593     chi2+=wr2;
 594     gamma+=wr2*(x[i]-xbar);
 595   }
 596
 597   /* Refined values of A and B. Also SA and SB. */
 598   Db=-ONEP5*gamma/Sxx;
 599   B+=Db;
 600   A-=ONEP5*chi2/sum-xbar*Db;
 601   SB=sqrt((chi2/(n-2))/Sxx);
 602   SA=SB*sqrt(Sxx/sum+sqr(xbar));
 603   *a=A;
 604   *b=B;
 605   *sa=SA;
 606   *sb=SB;
 607 }
 608
 609