src/tools/expfit.c

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