src/tools/expfit.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_expfit_c = "$Id$";
  30
  31 #include <sysstuff.h>
  32 #include <string.h>
  33 #include <math.h>
  34 #include "typedefs.h"
  35 #include "smalloc.h"
  36 #include "xvgr.h"
  37 #include "futil.h"
  38 #include "gstat.h"
  39 #include "vec.h"
  40 #include "statutil.h"
  41 #include "rdgroup.h"
  42
  43 extern void mrqmin(real x[],real y[],real sig[],int ndata,real a[],
  44                    int ma,int lista[],int mfit,real **covar,real **alpha,
  45                    real *chisq,
  46                    void (*funcs)(real x,real a[],real *y,real dyda[]),
  47                    real *alamda);
  48
  49 extern void mrqmin_new(real x[],real y[],real sig[],int ndata,real a[],
  50                        int ia[],int ma,real **covar,real **alpha,real *chisq,
  51                        void (*funcs)(real, real [], real *, real []),
  52                        real *alamda);
  53
  54 static real myexp(real x,real A,real tau)
  55 {
  56   if ((A == 0) || (tau == 0))
  57     return 0;
  58   return A*exp(-x/tau);
  59 }
  60
  61 static void exp_one_parm(real x,real a[],real *y,real dyda[])
  62 {
  63   /* Fit to function
  64    *
  65    * y = exp(-a1 x)
  66    *
  67    */
  68
  69   real e1;
  70
  71   e1      = exp(-x/a[1]);
  72   *y      = e1;
  73   dyda[1] = x*e1/(a[1]*a[1]);
  74 }
  75
  76 static void exp_two_parm(real x,real a[],real *y,real dyda[])
  77 {
  78   /* Fit to function
  79    *
  80    * y = a2 exp(-x/a1)
  81    *
  82    */
  83
  84   real e1;
  85
  86   e1      = exp(-x/a[1]);
  87   *y      = a[2]*e1;
  88   dyda[1] = x*a[2]*e1/(a[1]*a[1]);
  89   dyda[2] = e1;
  90 }
  91
  92 static void exp_3_parm(real x,real a[],real *y,real dyda[])
  93 {
  94   /* Fit to function
  95    *
  96    * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
  97    *
  98    */
  99
 100   real e1,e2;
 101
 102   e1      = exp(-x/a[1]);
 103   e2      = exp(-x/a[3]);
 104   *y      = a[2]*e1 + (1-a[2])*e2;
 105   dyda[1] = x*a[2]*e1/(a[1]*a[1]);
 106   dyda[2] = e1-e2;
 107   dyda[3] = x*(1-a[2])*e2/(a[3]*a[3]);
 108   /* fprintf(stderr,"exp3: x=%10.3e *y=%10.3e dyda=%10.3e %10.3e %10.3e\n",
 109     x,*y,dyda[1],dyda[2],dyda[3]);  */
 110 }
 111
 112 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
 113 static void lmfit_exp(int nframes,real x[],real y[],real dy[],real ftol,
 114                       real parm[],real dparm[],bool bVerbose,int nfitparm,
 115                       char *fix)
 116 {
 117   typedef void (*myfitfn)(real x,real a[],real *y,real dyda[]);
 118   myfitfn fitfn[3] = { exp_one_parm, exp_two_parm, exp_3_parm };
 119   real chisq,ochisq,alamda;
 120   real *a,**covar,**alpha;
 121   bool bCont;
 122   int  i,j,ma,mfit,*lista,*ia;
 123
 124   if ((nfitparm < 1) || (nfitparm > 3))
 125     fatal_error(0,"nfitparm = %d, should be in 1..3 (%s,%d)",
 126                 nfitparm,__FILE__,__LINE__);
 127
 128   ma=mfit=nfitparm;         /* number of parameters to fit */
 129   snew(a,ma+1);
 130   snew(covar,ma+1);
 131   snew(alpha,ma+1);
 132   snew(lista,ma+1);
 133   snew(ia,ma+1);
 134   for(i=1; (i<ma+1); i++) {
 135     lista[i] = i;
 136     ia[i] = i;
 137     snew(covar[i],ma+1);
 138     snew(alpha[i],ma+1);
 139   }
 140   if (fix != NULL) {
 141     fprintf(stderr,"Will keep %s fixed during fit procedure\n",fix);
 142     if (strcmp(fix,"tau1") == 0)
 143       ia[1]=0;
 144     else if (strcmp(fix,"A") == 0)
 145       ia[2]=0;
 146     else if (strcmp(fix,"tau2") == 0)
 147       ia[3]=0;
 148   }
 149   fprintf(stderr,"%d parameter fit\n",mfit);
 150
 151   /* Initial params */
 152   alamda = -1;    /* Starting value   */
 153   chisq  = 1e12;
 154   a[1]   = parm[0];    /* tau1 */
 155   if (nfitparm > 1)
 156     a[2]   = parm[1];  /* AA   */
 157   if (nfitparm > 2)
 158     a[3]   = parm[2];  /* tau2 */
 159
 160   j = 0;
 161   if (bVerbose)
 162     fprintf(stderr,"%4s  %10s  %10s  %10s  %10s  %10s\n",
 163             "Step","chi^2","Lambda","A1","A2","A3");
 164   do {
 165     ochisq = chisq;
 166     /* mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
 167      *   &chisq,fitfn[mfit-1],&alamda);
 168      */
 169     mrqmin_new(x-1,y-1,dy-1,nframes,a,ia,ma,covar,alpha,&chisq,
 170                fitfn[mfit-1],&alamda);
 171
 172     if (bVerbose) {
 173       fprintf(stderr,"%4d  %10.5e  %10.5e  %10.5e",
 174               j,chisq,alamda,a[1]);
 175       if (mfit > 1)
 176         fprintf(stderr,"  %10.5e",a[2]);
 177       if (mfit > 2)
 178         fprintf(stderr,"  %10.5e",a[3]);
 179       fprintf(stderr,"\n");
 180     }
 181     j++;
 182     bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
 183              ((ochisq == chisq)));
 184   } while (bCont && (alamda != 0.0) && (j < 50));
 185   if (bVerbose)
 186     fprintf(stderr,"\n");
 187
 188   /* Now get the covariance matrix out */
 189   alamda = 0;
 190
 191   /*  mrqmin(x-1,y-1,dy-1,nframes,a,ma,lista,mfit,covar,alpha,
 192    * &chisq,fitfn[mfit-1],&alamda);
 193    */
 194   mrqmin_new(x-1,y-1,dy-1,nframes,a,ia,ma,covar,alpha,&chisq,
 195              fitfn[mfit-1],&alamda);
 196
 197   for(j=0; (j<mfit); j++) {
 198     parm[j]  = a[j+1];
 199     dparm[j] = covar[j+1][j+1];
 200   }
 201
 202   for(i=0; (i<ma+1); i++) {
 203     sfree(covar[i]);
 204     sfree(alpha[i]);
 205   }
 206   sfree(a);
 207   sfree(covar);
 208   sfree(alpha);
 209   sfree(lista);
 210 }
 211
 212 real do_lmfit(int ndata,real c1[],real sig[],real dt,
 213               real begintimefit,real endtimefit,bool bVerbose,int nfitparm,
 214               real fit[],real fitparms[],char *fix)
 215 {
 216   FILE *fp;
 217   char buf[32];
 218
 219   int  i,j,nfitpnts;
 220   real integral,fittedfunc;
 221   real *parm,*dparm;
 222   real AA=0,tau1=0,tau2=0,srAA=0,srtau1,srtau2=0;
 223   real *x,*y,*dy;
 224   real ftol = 1e-4;
 225
 226   fprintf(stderr,"There are %d points to fit %d vars!\n",
 227           ndata,nfitparm);
 228   fprintf(stderr,"Fit from %g thru %g, dt=%g\n",
 229           begintimefit,endtimefit,dt);
 230
 231   snew(x,ndata);
 232   snew(y,ndata);
 233   snew(dy,ndata);
 234
 235   j=0;
 236   for(i=0; (i<ndata); i++) {
 237     if ( ((dt*i) >= begintimefit) && ((dt*i) <= endtimefit) ) {
 238       x[j]=dt*i;
 239       y[j]=c1[i];
 240
 241       /* mrqmin does not like sig to be zero */
 242       if (sig[i]<1.0e-7)
 243         sig[i]=1.0e-7;
 244       dy[j]=sig[i];
 245 #ifdef DEBUG
 246       fprintf(stderr,"j= %d, i= %d, x= %g, y= %g, dy= %g\n",
 247               j,i,x[j],y[j],dy[j]);
 248 #endif
 249       j++;
 250     }
 251   }
 252   nfitpnts=j;
 253   if (j < nfitparm) {
 254     fprintf(stderr,"Not enough data points for fitting!\n");
 255     integral = 0;
 256   }
 257   else {
 258     snew(parm,4);
 259     snew(dparm,4);
 260
 261     parm[0]=parm[1]=parm[2] = 1.0;
 262     if (fitparms) {
 263       parm[0]=fitparms[0];
 264       if (nfitparm == 2)
 265         parm[1]=fitparms[1];
 266       else {
 267         parm[1]=fitparms[1];
 268         parm[2]=fitparms[2];
 269       }
 270     }
 271
 272     lmfit_exp(nfitpnts,x,y,dy,ftol,parm,dparm,bVerbose,nfitparm,fix);
 273
 274     tau1 = parm[0];
 275     srtau1 = dparm[0];
 276     if (nfitparm > 1) {
 277       AA = parm[1];
 278       srAA = dparm[1];
 279     }
 280     else
 281       AA = 1.0;
 282     if (nfitparm > 2) {
 283       tau2 = parm[2];
 284       srtau2 = dparm[2];
 285     }
 286     else
 287       tau2 = 0.0;
 288
 289     /* Compute the integral from begintimefit to endtimefit
 290      */
 291     integral=(tau1*(myexp(begintimefit,AA,  tau1)   -
 292                     myexp(endtimefit,  AA,  tau1)) +
 293               tau2*(myexp(begintimefit,1-AA,tau2) -
 294                     myexp(endtimefit,  1-AA,tau2)));
 295
 296     /* Generate THE output */
 297     if (bVerbose) {
 298       fprintf(stderr,"FIT: # points used in fit is: %d\n",nfitpnts);
 299       fprintf(stderr,"FIT: %21s%21s%21s\n",
 300               "   A      ","tau1 (ps)    ","tau2 (ps)     ");
 301       fprintf(stderr,"FIT: ------------------------------------------------------------\n");
 302       fprintf(stderr,"FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
 303               AA,srAA,tau1,srtau1,tau2,srtau2);
 304       fprintf(stderr,"FIT: Integral (calc with fitted function) from %g to %g ps is: %g\n",
 305               begintimefit,endtimefit,integral);
 306
 307       sprintf(buf,"test%d.xvg",nfitpnts);
 308       fp = xvgropen(buf,"C(t) + Fit to C(t)","Time (ps)","C(t)");
 309       fprintf(fp,"# AA = %g, tau1 = %g, tau2 = %g\n",AA,tau1,tau2);
 310       for(j=0; (j<ndata); j++) {
 311         fittedfunc = myexp(dt*j,AA,tau1) + myexp(dt*j,1-AA,tau2);
 312         fprintf(fp,"%10.5e  %10.5e  %10.5e\n",dt*j,c1[j],fittedfunc);
 313       }
 314       fclose(fp);
 315     }
 316     /* Now copy the fitted func into the fit array, if not null */
 317     if (fit)
 318       for(i=0; (i<ndata); i++)
 319         fit[i] = myexp(i*dt,AA,tau1) + myexp(i*dt,1-AA,tau2);
 320   }
 321
 322   fitparms[0]=tau1;
 323   fitparms[1]=AA;
 324   fitparms[2]=tau2;
 325
 326   sfree(x);
 327   sfree(y);
 328   sfree(dy);
 329
 330   return integral;
 331 }
 332
 333 void do_expfit(int ndata,real c1[],real dt,real begintimefit,real endtimefit)
 334 {
 335   int i,n;
 336   real *x,*y,*Dy;
 337   real aa,bb,saa,sbb,A,tau,dA,dtau;
 338
 339   fprintf(stderr,"Will fit data from %g (ps) to %g (ps).\n",
 340           begintimefit,endtimefit);
 341
 342   snew(x,ndata);   /* allocate the maximum necessary space */
 343   snew(y,ndata);
 344   snew(Dy,ndata);
 345   n=0;
 346
 347   for(i=0; (i<ndata); i++) {
 348     if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) ) {
 349       x[n]=dt*i;
 350       y[n]=c1[i];
 351       Dy[n]=0.5;
 352       fprintf(stderr,"n= %d, i= %d, x= %g, y= %g\n",n,i,x[n],y[n]);
 353       n++;
 354     }
 355   }
 356   fprintf(stderr,"# of data points used in the fit is : %d\n\n",n);
 357   expfit(n,x,y,Dy,&aa,&saa,&bb,&sbb);
 358
 359   A=exp(aa);
 360   dA=exp(aa)*saa;
 361   tau=-1.0/bb;
 362   dtau=sbb/sqr(bb);
 363   fprintf(stderr,"Fitted to y=exp(a+bx):\n");
 364   fprintf(stderr,"a = %10.5f\t b = %10.5f",aa,bb);
 365   fprintf(stderr,"\n");
 366   fprintf(stderr,"Fitted to y=Aexp(-x/tau):\n");
 367   fprintf(stderr,"A  = %10.5f\t tau  = %10.5f\n",A,tau);
 368   fprintf(stderr,"dA = %10.5f\t dtau = %10.5f\n",dA,dtau);
 369 }
 370
 371
 372 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
 373             real *b, real *sb)
 374 {
 375   real *w,*ly,A,SA,B,SB;
 376   int  i;
 377   real sum,xbar,ybar,Sxx,Sxy,wr2,chi2,gamma,Db;
 378
 379 #define ZERO 0.0
 380 #define ONE 1.0
 381 #define ONEP5 1.5
 382 #define TWO 2.0
 383
 384 #define sqr(x) ((x)*(x))
 385
 386   /*allocate memory */
 387   snew(w,n);
 388   snew(ly,n);
 389
 390   /* Calculate weights and values of ln(y). */
 391   for(i=0;(i<n); i++){
 392     w[i]=sqr(y[i]/Dy[i]);
 393     ly[i]=log(y[i]);
 394   }
 395
 396   /* The weighted averages of x and y: xbar and ybar. */
 397   sum=ZERO;
 398   xbar=ZERO;
 399   ybar=ZERO;
 400   for(i=0;(i<n);i++){
 401     sum+=w[i];
 402     xbar+=w[i]*x[i];
 403     ybar+=w[i]*ly[i];
 404   }
 405   xbar/=sum;
 406   ybar/=sum;
 407
 408   /* The centered product sums Sxx and Sxy, and hence A and B. */
 409   Sxx=ZERO;
 410   Sxy=ZERO;
 411   for(i=0;(i<n);i++){
 412     Sxx+=w[i]*sqr(x[i]-xbar);
 413     Sxy+=w[i]*(x[i]-xbar)*(ly[i]-ybar);
 414   }
 415   B=Sxy/Sxx;
 416   A=ybar-B*xbar;
 417
 418   /* Chi-squared (chi2) and gamma. */
 419   chi2=ZERO;
 420   gamma=ZERO;
 421   for(i=0;(i<n);i++){
 422     wr2=w[i]*sqr(ly[i]-A-B*x[i]);
 423     chi2+=wr2;
 424     gamma+=wr2*(x[i]-xbar);
 425   }
 426
 427   /* Refined values of A and B. Also SA and SB. */
 428   Db=-ONEP5*gamma/Sxx;
 429   B+=Db;
 430   A-=ONEP5*chi2/sum-xbar*Db;
 431   SB=sqrt((chi2/(n-2))/Sxx);
 432   SA=SB*sqrt(Sxx/sum+sqr(xbar));
 433   *a=A;
 434   *b=B;
 435   *sa=SA;
 436   *sb=SB;
 437 }
 438
 439