2 * $Id: expfit.c,v 1.33 2005/08/29 19:39:11 lindahl Exp $
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * Green Red Orange Magenta Azure Cyan Skyblue
53 const int nfp_ffn[effnNR] = { 0, 1, 2, 3, 2, 5, 7, 9, 4, 3};
55 const char *s_ffn[effnNR+2] = { NULL, "none", "exp", "aexp", "exp_exp", "vac",
56 "exp5", "exp7", "exp9","erffit", NULL, NULL };
57 /* We don't allow errest as a choice on the command line */
59 const char *longs_ffn[effnNR] = {
63 "y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)",
64 "y = exp(-v) (cosh(wv) + 1/w sinh(wv)), v = x/(2 a1), w = sqrt(1 - a2)",
65 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5",
66 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7",
67 "y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7 exp(-x/a8) + a9",
68 "y = 1/2*(a1+a2) - 1/2*(a1-a2)*erf( (x-a3) /a4^2)",
69 "y = a2*ee(a1,x) + (1-a2)*ee(a2,x)"
72 extern gmx_bool mrqmin_new(real x[],real y[],real sig[],int ndata,real a[],
73 int ia[],int ma,real **covar,real **alpha,real *chisq,
74 void (*funcs)(real, real [], real *, real []),
77 static real myexp(real x,real A,real tau)
79 if ((A == 0) || (tau == 0))
84 void erffit (real x, real a[], real *y, real dyda[])
87 * y=(a1+a2)/2-(a1-a2)/2*erf((x-a3)/a4^2)
95 erfarg=(x-a[3])/(a[4]*a[4]);
96 erfarg2=erfarg*erfarg;
97 erfval=gmx_erf(erfarg)/2;
98 derf=M_2_SQRTPI*(a[1]-a[2])/2*exp(-erfarg2)/(a[4]*a[4]);
99 *y=(a[1]+a[2])/2-(a[1]-a[2])*erfval;
103 dyda[4]=2*derf*erfarg;
108 static void exp_one_parm(real x,real a[],real *y,real dyda[])
120 dyda[1] = x*e1/sqr(a[1]);
123 static void exp_two_parm(real x,real a[],real *y,real dyda[])
135 dyda[1] = x*a[2]*e1/sqr(a[1]);
139 static void exp_3_parm(real x,real a[],real *y,real dyda[])
143 * y = a2 exp(-x/a1) + (1-a2) exp(-x/a3)
151 *y = a[2]*e1 + (1-a[2])*e2;
152 dyda[1] = x*a[2]*e1/sqr(a[1]);
154 dyda[3] = x*(1-a[2])*e2/sqr(a[3]);
157 static void exp_5_parm(real x,real a[],real *y,real dyda[])
161 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5
169 *y = a[1]*e1 + a[3]*e2 + a[5];
172 fprintf(debug,"exp_5_parm called: x = %10.3f y = %10.3f\n"
173 "a = ( %8.3f %8.3f %8.3f %8.3f %8.3f)\n",
174 x,*y,a[1],a[2],a[3],a[4],a[5]);
176 dyda[2] = x*e1/sqr(a[2]);
178 dyda[4] = x*e2/sqr(a[4]);
182 static void exp_7_parm(real x,real a[],real *y,real dyda[])
186 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
195 *y = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7];
198 dyda[2] = x*e1/sqr(a[2]);
200 dyda[4] = x*e2/sqr(a[4]);
202 dyda[6] = x*e3/sqr(a[6]);
206 static void exp_9_parm(real x,real a[],real *y,real dyda[])
210 * y = a1 exp(-x/a2) + a3 exp(-x/a4) + a5 exp(-x/a6) + a7
220 *y = a[1]*e1 + a[3]*e2 + a[5]*e3 + a[7]*e4 + a[9];
223 dyda[2] = x*e1/sqr(a[2]);
225 dyda[4] = x*e2/sqr(a[4]);
227 dyda[6] = x*e3/sqr(a[6]);
229 dyda[8] = x*e4/sqr(a[8]);
233 static void vac_2_parm(real x,real a[],real *y,real dyda[])
237 * y = 1/2 (1 - 1/w) exp(-(1+w)v) + 1/2 (1 + 1/w) exp(-(1-w)v)
239 * = exp(-v) (cosh(wv) + 1/w sinh(wv))
244 * For tranverse current autocorrelation functions:
246 * a2 = 4 tau (eta/rho) k^2
250 double v,det,omega,omega2,em,ec,es;
257 omega = sqrt(omega2);
259 ec = em*0.5*(exp(omega*v)+exp(-omega*v));
260 es = em*0.5*(exp(omega*v)-exp(-omega*v))/omega;
262 ec = em*cos(omega*v);
263 es = em*sin(omega*v)/omega;
266 dyda[2] = (v/det*ec+(v-1/det)*es)/(-2.0);
267 dyda[1] = (1-det)*v/a[1]*es;
270 dyda[2] = -v*v*em*(0.5+v/6);
271 dyda[1] = v*v/a[1]*em;
275 static void errest_3_parm(real x,real a[],real *y,real dyda[])
280 e1 = exp(-x/a[1]) - 1;
284 e2 = exp(-x/a[3]) - 1;
289 v1 = 2*a[1]*(e1*a[1]/x + 1);
290 v2 = 2*a[3]*(e2*a[3]/x + 1);
291 *y = a[2]*v1 + (1-a[2])*v2;
292 dyda[1] = 2* a[2] *(v1/a[1] + e1);
293 dyda[3] = 2*(1 - a[2])*(v2/a[3] + e2);
303 typedef void (*myfitfn)(real x,real a[],real *y,real dyda[]);
304 myfitfn mfitfn[effnNR] = {
305 exp_one_parm, exp_one_parm, exp_two_parm, exp_3_parm, vac_2_parm,
306 exp_5_parm, exp_7_parm, exp_9_parm, erffit, errest_3_parm
309 real fit_function(int eFitFn,real *parm,real x)
311 static real y,dum[8];
313 mfitfn[eFitFn](x,parm-1,&y,dum);
318 /* lmfit_exp supports up to 3 parameter fitting of exponential functions */
319 static gmx_bool lmfit_exp(int nfit,real x[],real y[],real dy[],real ftol,
320 real parm[],real dparm[],gmx_bool bVerbose,
323 real chisq,ochisq,alamda;
324 real *a,**covar,**alpha,*dum;
326 int i,j,ma,mfit,*lista,*ia;
328 if ((eFitFn < 0) || (eFitFn >= effnNR))
329 gmx_fatal(FARGS,"fitfn = %d, should be in 0..%d (%s,%d)",
330 effnNR-1,eFitFn,__FILE__,__LINE__);
332 ma=mfit=nfp_ffn[eFitFn]; /* number of parameters to fit */
339 for(i=1; (i<ma+1); i++) {
341 ia[i] = 1; /* fixed bug B.S.S 19/11 */
347 fprintf(stderr,"Will keep parameters fixed during fit procedure: %d\n",
354 fprintf(debug,"%d parameter fit\n",mfit);
357 alamda = -1; /* Starting value */
359 for(i=0; (i<mfit); i++)
364 fprintf(stderr,"%4s %10s %10s %10s %10s %10s %10s\n",
365 "Step","chi^2","Lambda","A1","A2","A3","A4");
368 /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
369 * &chisq,expfn[mfit-1],&alamda)
371 if (!mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
372 mfitfn[eFitFn],&alamda))
376 fprintf(stderr,"%4d %10.5e %10.5e %10.5e",
377 j,chisq,alamda,a[1]);
379 fprintf(stderr," %10.5e",a[2]);
381 fprintf(stderr," %10.5e",a[3]);
383 fprintf(stderr," %10.5e",a[4]);
384 fprintf(stderr,"\n");
387 bCont = ((fabs(ochisq - chisq) > fabs(ftol*chisq)) ||
388 ((ochisq == chisq)));
389 } while (bCont && (alamda != 0.0) && (j < 50));
391 fprintf(stderr,"\n");
393 /* Now get the covariance matrix out */
396 /* mrqmin(x-1,y-1,dy-1,nfit,a,ma,lista,mfit,covar,alpha,
397 * &chisq,expfn[mfit-1],&alamda)
399 if ( !mrqmin_new(x-1,y-1,dy-1,nfit,a,ia,ma,covar,alpha,&chisq,
400 mfitfn[eFitFn],&alamda))
403 for(j=0; (j<mfit); j++) {
405 dparm[j] = covar[j+1][j+1];
408 for(i=0; (i<ma+1); i++) {
421 real do_lmfit(int ndata,real c1[],real sig[],real dt,real x0[],
422 real begintimefit,real endtimefit,const output_env_t oenv,
423 gmx_bool bVerbose, int eFitFn,real fitparms[],int fix)
428 int i,j,nparm,nfitpnts;
434 nparm = nfp_ffn[eFitFn];
436 fprintf(debug,"There are %d points to fit %d vars!\n",ndata,nparm);
437 fprintf(debug,"Fit to function %d from %g through %g, dt=%g\n",
438 eFitFn,begintimefit,endtimefit,dt);
446 for(i=0; i<ndata; i++) {
447 ttt = x0 ? x0[i] : dt*i;
448 if (ttt>=begintimefit && ttt<=endtimefit) {
452 /* mrqmin does not like sig to be zero */
458 fprintf(debug,"j= %d, i= %d, x= %g, y= %g, dy= %g\n",
459 j,i,x[j],y[j],dy[j]);
465 if (nfitpnts < nparm)
466 fprintf(stderr,"Not enough data points for fitting!\n");
471 for(i=0; (i < nparm); i++)
474 if (!lmfit_exp(nfitpnts,x,y,dy,ftol,parm,dparm,bVerbose,eFitFn,fix))
475 fprintf(stderr,"Fit failed!\n");
476 else if (nparm <= 3) {
477 /* Compute the integral from begintimefit */
479 integral=(parm[0]*myexp(begintimefit,parm[1], parm[0]) +
480 parm[2]*myexp(begintimefit,1-parm[1],parm[2]));
482 integral=parm[0]*myexp(begintimefit,parm[1], parm[0]);
484 integral=parm[0]*myexp(begintimefit,1, parm[0]);
486 gmx_fatal(FARGS,"nparm = %d in file %s, line %d",
487 nparm,__FILE__,__LINE__);
489 /* Generate THE output */
491 fprintf(stderr,"FIT: # points used in fit is: %d\n",nfitpnts);
492 fprintf(stderr,"FIT: %21s%21s%21s\n",
493 "parm0 ","parm1 (ps) ","parm2 (ps) ");
494 fprintf(stderr,"FIT: ------------------------------------------------------------\n");
495 fprintf(stderr,"FIT: %8.3g +/- %8.3g%9.4g +/- %8.3g%8.3g +/- %8.3g\n",
496 parm[0],dparm[0],parm[1],dparm[1],parm[2],dparm[2]);
497 fprintf(stderr,"FIT: Integral (calc with fitted function) from %g ps to inf. is: %g\n",
498 begintimefit,integral);
500 sprintf(buf,"test%d.xvg",nfitpnts);
501 fp = xvgropen(buf,"C(t) + Fit to C(t)","Time (ps)","C(t)",oenv);
502 fprintf(fp,"# parm0 = %g, parm1 = %g, parm2 = %g\n",
503 parm[0],parm[1],parm[2]);
504 for(j=0; (j<nfitpnts); j++) {
505 ttt = x0 ? x0[j] : dt*j;
506 fprintf(fp,"%10.5e %10.5e %10.5e\n",
507 ttt,c1[j],fit_function(eFitFn,parm,ttt));
514 for(i=0;(i<nparm);i++)
516 fitparms[i] = parm[i];
530 void do_expfit(int ndata,real c1[],real dt,real begintimefit,real endtimefit)
534 real aa,bb,saa,sbb,A,tau,dA,dtau;
536 fprintf(stderr,"Will fit data from %g (ps) to %g (ps).\n",
537 begintimefit,endtimefit);
539 snew(x,ndata); /* allocate the maximum necessary space */
544 for(i=0; (i<ndata); i++) {
545 if ( (dt*i >= begintimefit) && (dt*i <= endtimefit) ) {
549 fprintf(stderr,"n= %d, i= %d, x= %g, y= %g\n",n,i,x[n],y[n]);
553 fprintf(stderr,"# of data points used in the fit is : %d\n\n",n);
554 expfit(n,x,y,Dy,&aa,&saa,&bb,&sbb);
560 fprintf(stderr,"Fitted to y=exp(a+bx):\n");
561 fprintf(stderr,"a = %10.5f\t b = %10.5f",aa,bb);
562 fprintf(stderr,"\n");
563 fprintf(stderr,"Fitted to y=Aexp(-x/tau):\n");
564 fprintf(stderr,"A = %10.5f\t tau = %10.5f\n",A,tau);
565 fprintf(stderr,"dA = %10.5f\t dtau = %10.5f\n",dA,dtau);
569 void expfit(int n, real *x, real *y, real *Dy, real *a, real *sa,
572 real *w,*ly,A,SA,B,SB;
574 real sum,xbar,ybar,Sxx,Sxy,wr2,chi2,gamma,Db;
581 #define sqr(x) ((x)*(x))
587 /* Calculate weights and values of ln(y). */
589 w[i]=sqr(y[i]/Dy[i]);
593 /* The weighted averages of x and y: xbar and ybar. */
605 /* The centered product sums Sxx and Sxy, and hence A and B. */
609 Sxx+=w[i]*sqr(x[i]-xbar);
610 Sxy+=w[i]*(x[i]-xbar)*(ly[i]-ybar);
615 /* Chi-squared (chi2) and gamma. */
619 wr2=w[i]*sqr(ly[i]-A-B*x[i]);
621 gamma+=wr2*(x[i]-xbar);
624 /* Refined values of A and B. Also SA and SB. */
627 A-=ONEP5*chi2/sum-xbar*Db;
628 SB=sqrt((chi2/(n-2))/Sxx);
629 SA=SB*sqrt(Sxx/sum+sqr(xbar));