src/gromacs/gmxlib/statistics/gmx_statistics.c

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*- */
   2 /*
   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 #include <math.h>
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "vec.h"
  44 #include "gmx_statistics.h"
  45
  46 static int gmx_dnint(double x)
  47 {
  48     return (int) (x+0.5);
  49 }
  50
  51 typedef struct gmx_stats {
  52     double aa,a,b,sigma_aa,sigma_a,sigma_b,aver,sigma_aver,error;
  53     double rmsd,Rdata,Rfit,Rfitaa,chi2,chi2aa;
  54     double *x,*y,*dx,*dy;
  55     int    computed;
  56     int    np,np_c,nalloc;
  57 } gmx_stats;
  58
  59 gmx_stats_t gmx_stats_init()
  60 {
  61     gmx_stats *stats;
  62
  63     snew(stats,1);
  64
  65     return (gmx_stats_t) stats;
  66 }
  67
  68 int gmx_stats_get_npoints(gmx_stats_t gstats, int *N)
  69 {
  70     gmx_stats *stats = (gmx_stats *) gstats;
  71
  72     *N = stats->np;
  73
  74     return estatsOK;
  75 }
  76
  77 int gmx_stats_done(gmx_stats_t gstats)
  78 {
  79     gmx_stats *stats = (gmx_stats *) gstats;
  80
  81     sfree(stats->x);
  82     stats->x = NULL;
  83     sfree(stats->y);
  84     stats->y = NULL;
  85     sfree(stats->dx);
  86     stats->dx = NULL;
  87     sfree(stats->dy);
  88     stats->dy = NULL;
  89
  90     return estatsOK;
  91 }
  92
  93 int gmx_stats_add_point(gmx_stats_t gstats,double x,double y,
  94                         double dx,double dy)
  95 {
  96     gmx_stats *stats = (gmx_stats *) gstats;
  97     int i;
  98
  99     if (stats->np+1 >= stats->nalloc)
 100     {
 101         if (stats->nalloc == 0)
 102             stats->nalloc = 1024;
 103         else
 104             stats->nalloc *= 2;
 105         srenew(stats->x,stats->nalloc);
 106         srenew(stats->y,stats->nalloc);
 107         srenew(stats->dx,stats->nalloc);
 108         srenew(stats->dy,stats->nalloc);
 109         for(i=stats->np; (i<stats->nalloc); i++)
 110         {
 111             stats->x[i]  = 0;
 112             stats->y[i]  = 0;
 113             stats->dx[i] = 0;
 114             stats->dy[i] = 0;
 115         }
 116     }
 117     stats->x[stats->np]  = x;
 118     stats->y[stats->np]  = y;
 119     stats->dx[stats->np] = dx;
 120     stats->dy[stats->np] = dy;
 121     stats->np++;
 122     stats->computed = 0;
 123
 124     return estatsOK;
 125 }
 126
 127 int gmx_stats_get_point(gmx_stats_t gstats,real *x,real *y,
 128                         real *dx,real *dy,real level)
 129 {
 130     gmx_stats *stats = (gmx_stats *) gstats;
 131     int  ok,outlier;
 132     real rmsd,r;
 133
 134     if ((ok = gmx_stats_get_rmsd(gstats,&rmsd)) != estatsOK)
 135     {
 136         return ok;
 137     }
 138     outlier = 0;
 139     while ((outlier == 0) && (stats->np_c < stats->np))
 140     {
 141         r = fabs(stats->x[stats->np_c] - stats->y[stats->np_c]);
 142         outlier = (r > rmsd*level);
 143         if (outlier)
 144         {
 145             if (NULL != x)  *x  = stats->x[stats->np_c];
 146             if (NULL != y)  *y  = stats->y[stats->np_c];
 147             if (NULL != dx) *dx = stats->dx[stats->np_c];
 148             if (NULL != dy) *dy = stats->dy[stats->np_c];
 149         }
 150         stats->np_c++;
 151
 152         if (outlier)
 153             return estatsOK;
 154     }
 155
 156     stats->np_c = 0;
 157
 158     return estatsNO_POINTS;
 159 }
 160
 161 int gmx_stats_add_points(gmx_stats_t gstats,int n,real *x,real *y,
 162                          real *dx,real *dy)
 163 {
 164     int i,ok;
 165
 166     for(i=0; (i<n); i++)
 167     {
 168         if ((ok = gmx_stats_add_point(gstats,x[i],y[i],
 169                                       (NULL != dx) ? dx[i] : 0,
 170                                       (NULL != dy) ? dy[i] : 0)) != estatsOK)
 171         {
 172             return ok;
 173         }
 174     }
 175     return estatsOK;
 176 }
 177
 178 static int gmx_stats_compute(gmx_stats *stats,int weight)
 179 {
 180     double yy,yx,xx,sx,sy,dy,chi2,chi2aa,d2;
 181     double ssxx,ssyy,ssxy;
 182     double w,wtot,yx_nw,sy_nw,sx_nw,yy_nw,xx_nw,dx2,dy2;
 183     int i,N;
 184
 185     N = stats->np;
 186     if (stats->computed == 0)
 187     {
 188         if (N < 1)
 189         {
 190             return estatsNO_POINTS;
 191         }
 192
 193         xx = xx_nw = 0;
 194         yy = yy_nw = 0;
 195         yx = yx_nw = 0;
 196         sx = sx_nw = 0;
 197         sy = sy_nw = 0;
 198         wtot = 0;
 199         d2   = 0;
 200         for(i=0; (i<N); i++)
 201         {
 202             d2 += dsqr(stats->x[i]-stats->y[i]);
 203             if ((stats->dy[i]) && (weight == elsqWEIGHT_Y))
 204             {
 205                 w = 1/dsqr(stats->dy[i]);
 206             }
 207             else
 208             {
 209                 w = 1;
 210             }
 211
 212             wtot  += w;
 213
 214             xx    += w*dsqr(stats->x[i]);
 215             xx_nw += dsqr(stats->x[i]);
 216
 217             yy    += w*dsqr(stats->y[i]);
 218             yy_nw += dsqr(stats->y[i]);
 219
 220             yx    += w*stats->y[i]*stats->x[i];
 221             yx_nw += stats->y[i]*stats->x[i];
 222
 223             sx    += w*stats->x[i];
 224             sx_nw += stats->x[i];
 225
 226             sy    += w*stats->y[i];
 227             sy_nw += stats->y[i];
 228         }
 229
 230         /* Compute average, sigma and error */
 231         stats->aver       = sy_nw/N;
 232         stats->sigma_aver = sqrt(yy_nw/N - dsqr(sy_nw/N));
 233         stats->error      = stats->sigma_aver/sqrt(N);
 234
 235         /* Compute RMSD between x and y */
 236         stats->rmsd = sqrt(d2/N);
 237
 238         /* Correlation coefficient for data */
 239         yx_nw /= N;
 240         xx_nw /= N;
 241         yy_nw /= N;
 242         sx_nw /= N;
 243         sy_nw /= N;
 244         ssxx = N*(xx_nw - dsqr(sx_nw));
 245         ssyy = N*(yy_nw - dsqr(sy_nw));
 246         ssxy = N*(yx_nw - (sx_nw*sy_nw));
 247         stats->Rdata = sqrt(dsqr(ssxy)/(ssxx*ssyy));
 248
 249         /* Compute straight line through datapoints, either with intercept
 250            zero (result in aa) or with intercept variable (results in a
 251            and b) */
 252         yx = yx/wtot;
 253         xx = xx/wtot;
 254         sx = sx/wtot;
 255         sy = sy/wtot;
 256
 257         stats->aa = (yx/xx);
 258         stats->a  = (yx-sx*sy)/(xx-sx*sx);
 259         stats->b  = (sy)-(stats->a)*(sx);
 260
 261         /* Compute chi2, deviation from a line y = ax+b. Also compute
 262            chi2aa which returns the deviation from a line y = ax. */
 263         chi2   = 0;
 264         chi2aa = 0;
 265         for(i=0; (i<N); i++)
 266         {
 267             if (stats->dy[i] > 0)
 268             {
 269                 dy = stats->dy[i];
 270             }
 271             else
 272             {
 273                 dy = 1;
 274             }
 275             chi2aa += dsqr((stats->y[i]-(stats->aa*stats->x[i]))/dy);
 276             chi2   += dsqr((stats->y[i]-(stats->a*stats->x[i]+stats->b))/dy);
 277         }
 278         if (N > 2)
 279         {
 280             stats->chi2   = sqrt(chi2/(N-2));
 281             stats->chi2aa = sqrt(chi2aa/(N-2));
 282
 283             /* Look up equations! */
 284             dx2 = (xx-sx*sx);
 285             dy2 = (yy-sy*sy);
 286             stats->sigma_a = sqrt(stats->chi2/((N-2)*dx2));
 287             stats->sigma_b = stats->sigma_a*sqrt(xx);
 288             stats->Rfit    = fabs(ssxy)/sqrt(ssxx*ssyy);
 289                 /*stats->a*sqrt(dx2/dy2);*/
 290             stats->Rfitaa  = stats->aa*sqrt(dx2/dy2);
 291         }
 292         else
 293         {
 294             stats->chi2    = 0;
 295             stats->chi2aa  = 0;
 296             stats->sigma_a = 0;
 297             stats->sigma_b = 0;
 298             stats->Rfit    = 0;
 299             stats->Rfitaa  = 0;
 300         }
 301
 302         stats->computed = 1;
 303     }
 304
 305     return estatsOK;
 306 }
 307
 308 int gmx_stats_get_ab(gmx_stats_t gstats,int weight,
 309                      real *a,real *b,real *da,real *db,
 310                      real *chi2,real *Rfit)
 311 {
 312     gmx_stats *stats = (gmx_stats *) gstats;
 313     int ok;
 314
 315     if ((ok = gmx_stats_compute(stats,weight)) != estatsOK)
 316         return ok;
 317     if (NULL != a)
 318     {
 319         *a    = stats->a;
 320     }
 321     if (NULL != b)
 322     {
 323         *b    = stats->b;
 324     }
 325     if (NULL != da)
 326     {
 327         *da   = stats->sigma_a;
 328     }
 329     if (NULL != db)
 330     {
 331         *db   = stats->sigma_b;
 332     }
 333     if (NULL != chi2)
 334     {
 335         *chi2 = stats->chi2;
 336     }
 337     if (NULL != Rfit)
 338     {
 339         *Rfit = stats->Rfit;
 340     }
 341
 342     return estatsOK;
 343 }
 344
 345 int gmx_stats_get_a(gmx_stats_t gstats,int weight,real *a,real *da,
 346                     real *chi2,real *Rfit)
 347 {
 348     gmx_stats *stats = (gmx_stats *) gstats;
 349     int ok;
 350
 351     if ((ok = gmx_stats_compute(stats,weight)) != estatsOK)
 352         return ok;
 353     if (NULL != a)    *a    = stats->aa;
 354     if (NULL != da)   *da   = stats->sigma_aa;
 355     if (NULL != chi2) *chi2 = stats->chi2aa;
 356     if (NULL != Rfit) *Rfit = stats->Rfitaa;
 357
 358     return estatsOK;
 359 }
 360
 361 int gmx_stats_get_average(gmx_stats_t gstats,real *aver)
 362 {
 363     gmx_stats *stats = (gmx_stats *) gstats;
 364     int ok;
 365
 366     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 367     {
 368         return ok;
 369     }
 370
 371     *aver = stats->aver;
 372
 373     return estatsOK;
 374 }
 375
 376 int gmx_stats_get_ase(gmx_stats_t gstats,real *aver,real *sigma,real *error)
 377 {
 378     gmx_stats *stats = (gmx_stats *) gstats;
 379     int ok;
 380
 381     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 382     {
 383         return ok;
 384     }
 385
 386     if (NULL != aver)
 387     {
 388         *aver  = stats->aver;
 389     }
 390     if (NULL != sigma)
 391     {
 392         *sigma = stats->sigma_aver;
 393     }
 394     if (NULL != error)
 395     {
 396         *error = stats->error;
 397     }
 398
 399     return estatsOK;
 400 }
 401
 402 int gmx_stats_get_sigma(gmx_stats_t gstats,real *sigma)
 403 {
 404     gmx_stats *stats = (gmx_stats *) gstats;
 405     int ok;
 406
 407     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 408         return ok;
 409
 410     *sigma = stats->sigma_aver;
 411
 412     return estatsOK;
 413 }
 414
 415 int gmx_stats_get_error(gmx_stats_t gstats,real *error)
 416 {
 417     gmx_stats *stats = (gmx_stats *) gstats;
 418     int ok;
 419
 420     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 421         return ok;
 422
 423     *error = stats->error;
 424
 425     return estatsOK;
 426 }
 427
 428 int gmx_stats_get_corr_coeff(gmx_stats_t gstats,real *R)
 429 {
 430     gmx_stats *stats = (gmx_stats *) gstats;
 431     int ok;
 432
 433     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 434         return ok;
 435
 436     *R = stats->Rdata;
 437
 438     return estatsOK;
 439 }
 440
 441 int gmx_stats_get_rmsd(gmx_stats_t gstats,real *rmsd)
 442 {
 443     gmx_stats *stats = (gmx_stats *) gstats;
 444     int ok;
 445
 446     if ((ok = gmx_stats_compute(stats,elsqWEIGHT_NONE)) != estatsOK)
 447     {
 448         return ok;
 449     }
 450
 451     *rmsd = stats->rmsd;
 452
 453     return estatsOK;
 454 }
 455
 456 int gmx_stats_dump_xy(gmx_stats_t gstats,FILE *fp)
 457 {
 458     gmx_stats *stats = (gmx_stats *) gstats;
 459     int i,ok;
 460
 461     for(i=0; (i<stats->np); i++)
 462     {
 463         fprintf(fp,"%12g  %12g  %12g  %12g\n",stats->x[i],stats->y[i],
 464                 stats->dx[i],stats->dy[i]);
 465     }
 466
 467     return estatsOK;
 468 }
 469
 470 int gmx_stats_remove_outliers(gmx_stats_t gstats,double level)
 471 {
 472     gmx_stats *stats = (gmx_stats *) gstats;
 473     int  i,iter=1,done=0,ok;
 474     real rmsd,r;
 475
 476     while ((stats->np >= 10) && !done)
 477     {
 478         if ((ok = gmx_stats_get_rmsd(gstats,&rmsd)) != estatsOK)
 479         {
 480             return ok;
 481         }
 482         done = 1;
 483         for(i=0; (i<stats->np); )
 484         {
 485             r = fabs(stats->x[i]-stats->y[i]);
 486             if (r > level*rmsd)
 487             {
 488                 fprintf(stderr,"Removing outlier, iter = %d, rmsd = %g, x = %g, y = %g\n",
 489                         iter,rmsd,stats->x[i],stats->y[i]);
 490                 if (i < stats->np-1)
 491                 {
 492                     stats->x[i]  = stats->x[stats->np-1];
 493                     stats->y[i]  = stats->y[stats->np-1];
 494                     stats->dx[i] = stats->dx[stats->np-1];
 495                     stats->dy[i] = stats->dy[stats->np-1];
 496                 }
 497                 stats->np--;
 498                 done = 0;
 499             }
 500             else
 501             {
 502                 i++;
 503             }
 504         }
 505         iter++;
 506     }
 507
 508     return estatsOK;
 509 }
 510
 511 int gmx_stats_make_histogram(gmx_stats_t gstats,real binwidth,int *nb,
 512                              int ehisto,int normalized,real **x,real **y)
 513 {
 514     gmx_stats *stats = (gmx_stats *) gstats;
 515     int    i,ok,index=0,nbins=*nb,*nindex;
 516     double minx,maxx,maxy,miny,delta,dd,minh;
 517
 518     if (((binwidth <= 0) && (nbins <= 0)) ||
 519         ((binwidth > 0) && (nbins > 0)))
 520     {
 521         return estatsINVALID_INPUT;
 522     }
 523     if (stats->np <= 2)
 524     {
 525         return estatsNO_POINTS;
 526     }
 527     minx = maxx = stats->x[0];
 528     miny = maxy = stats->y[0];
 529     for(i=1; (i<stats->np); i++)
 530     {
 531         miny = (stats->y[i] < miny) ? stats->y[i] : miny;
 532         maxy = (stats->y[i] > maxy) ? stats->y[i] : maxy;
 533         minx = (stats->x[i] < minx) ? stats->x[i] : minx;
 534         maxx = (stats->x[i] > maxx) ? stats->x[i] : maxx;
 535     }
 536     if (ehisto == ehistoX)
 537     {
 538         delta = maxx-minx;
 539         minh = minx;
 540     }
 541     else if (ehisto == ehistoY)
 542     {
 543         delta = maxy-miny;
 544         minh = miny;
 545     }
 546     else
 547         return estatsINVALID_INPUT;
 548
 549     if (binwidth == 0)
 550     {
 551         binwidth = (delta)/nbins;
 552     }
 553     else
 554     {
 555         nbins = gmx_dnint((delta)/binwidth + 0.5);
 556     }
 557     snew(*x,nbins);
 558     snew(nindex,nbins);
 559     for(i=0; (i<nbins); i++)
 560     {
 561         (*x)[i] = minh + binwidth*(i+0.5);
 562     }
 563     if (normalized == 0)
 564     {
 565         dd = 1;
 566     }
 567     else
 568     {
 569         dd = 1.0/(binwidth*stats->np);
 570     }
 571
 572     snew(*y,nbins);
 573     for(i=0; (i<stats->np); i++)
 574     {
 575         if (ehisto == ehistoY)
 576             index = (stats->y[i]-miny)/binwidth;
 577                 else if (ehisto == ehistoX)
 578             index = (stats->x[i]-minx)/binwidth;
 579                 if (index<0)
 580                 {
 581                         index = 0;
 582                 }
 583                 if (index>nbins-1)
 584                 {
 585                         index = nbins-1;
 586                 }
 587         (*y)[index] += dd;
 588         nindex[index]++;
 589     }
 590     if (*nb == 0)
 591         *nb = nbins;
 592     for(i=0; (i<nbins); i++)
 593         if (nindex[i] > 0)
 594             (*y)[i] /= nindex[i];
 595
 596     sfree(nindex);
 597
 598     return estatsOK;
 599 }
 600
 601 static const char *stats_error[estatsNR] =
 602 {
 603     "All well in STATS land",
 604     "No points",
 605     "Not enough memory",
 606     "Invalid histogram input",
 607     "Unknown error",
 608     "Not implemented yet"
 609 };
 610
 611 const char *gmx_stats_message(int estats)
 612 {
 613     if ((estats >= 0) && (estats < estatsNR))
 614     {
 615         return stats_error[estats];
 616     }
 617     else
 618     {
 619         return stats_error[estatsERROR];
 620     }
 621 }
 622
 623 /* Old convenience functions, should be merged with the core
 624    statistics above. */
 625 int lsq_y_ax(int n, real x[], real y[], real *a)
 626 {
 627     gmx_stats_t lsq = gmx_stats_init();
 628     int  ok;
 629     real da,chi2,Rfit;
 630
 631     gmx_stats_add_points(lsq,n,x,y,0,0);
 632     if ((ok = gmx_stats_get_a(lsq,elsqWEIGHT_NONE,a,&da,&chi2,&Rfit)) != estatsOK)
 633     {
 634         return ok;
 635     }
 636
 637     /*  int    i;
 638         double xx,yx;
 639
 640         yx=xx=0.0;
 641         for (i=0; i<n; i++) {
 642         yx+=y[i]*x[i];
 643         xx+=x[i]*x[i];
 644         }
 645         *a=yx/xx;
 646         */
 647     return estatsOK;
 648 }
 649
 650 static int low_lsq_y_ax_b(int n, real *xr, double *xd, real yr[],
 651                           real *a, real *b,real *r,real *chi2)
 652 {
 653     int    i,ok;
 654     gmx_stats_t lsq;
 655
 656     lsq = gmx_stats_init();
 657     for(i=0; (i<n); i++)
 658     {
 659         if ((ok = gmx_stats_add_point(lsq,(NULL != xd) ? xd[i] : xr[i],yr[i],0,0))
 660             != estatsOK)
 661         {
 662             return ok;
 663         }
 664     }
 665     if ((ok = gmx_stats_get_ab(lsq,elsqWEIGHT_NONE,a,b,NULL,NULL,chi2,r)) != estatsOK)
 666     {
 667         return ok;
 668     }
 669
 670     return estatsOK;
 671     /*
 672       double x,y,yx,xx,yy,sx,sy,chi2;
 673
 674       yx=xx=yy=sx=sy=0.0;
 675       for (i=0; i<n; i++) {
 676       if (xd != NULL) {
 677       x = xd[i];
 678       } else {
 679       x = xr[i];
 680       }
 681       y =   yr[i];
 682
 683       yx += y*x;
 684       xx += x*x;
 685       yy += y*y;
 686       sx += x;
 687       sy += y;
 688       }
 689       *a = (n*yx-sy*sx)/(n*xx-sx*sx);
 690       *b = (sy-(*a)*sx)/n;
 691       *r = sqrt((xx-sx*sx)/(yy-sy*sy));
 692
 693       chi2 = 0;
 694       if (xd != NULL) {
 695       for(i=0; i<n; i++)
 696       chi2 += dsqr(yr[i] - ((*a)*xd[i] + (*b)));
 697       } else {
 698       for(i=0; i<n; i++)
 699       chi2 += dsqr(yr[i] - ((*a)*xr[i] + (*b)));
 700       }
 701
 702       if (n > 2)
 703       return sqrt(chi2/(n-2));
 704       else
 705       return 0;
 706     */
 707 }
 708
 709 int lsq_y_ax_b(int n, real x[], real y[], real *a, real *b,real *r,real *chi2)
 710 {
 711     return low_lsq_y_ax_b(n,x,NULL,y,a,b,r,chi2);
 712 }
 713
 714 int lsq_y_ax_b_xdouble(int n, double x[], real y[], real *a, real *b,
 715                        real *r,real *chi2)
 716 {
 717     return low_lsq_y_ax_b(n,NULL,x,y,a,b,r,chi2);
 718 }
 719
 720 int lsq_y_ax_b_error(int n, real x[], real y[], real dy[],
 721                      real *a, real *b, real *da, real *db,
 722                      real *r,real *chi2)
 723 {
 724     gmx_stats_t lsq;
 725     int    i,ok;
 726
 727     lsq = gmx_stats_init();
 728     for(i=0; (i<n); i++)
 729     {
 730         if ((ok = gmx_stats_add_point(lsq,x[i],y[i],0,dy[i])) != estatsOK)
 731         {
 732             return ok;
 733         }
 734     }
 735     if ((ok = gmx_stats_get_ab(lsq,elsqWEIGHT_Y,a,b,da,db,chi2,r)) != estatsOK)
 736     {
 737         return ok;
 738     }
 739     if ((ok = gmx_stats_done(lsq)) != estatsOK)
 740     {
 741         return ok;
 742     }
 743     sfree(lsq);
 744
 745     return estatsOK;
 746     /*
 747       double sxy,sxx,syy,sx,sy,w,s_2,dx2,dy2,mins;
 748
 749       sxy=sxx=syy=sx=sy=w=0.0;
 750       mins = dy[0];
 751       for(i=1; (i<n); i++)
 752       mins = min(mins,dy[i]);
 753       if (mins <= 0)
 754       gmx_fatal(FARGS,"Zero or negative weigths in linear regression analysis");
 755
 756       for (i=0; i<n; i++) {
 757       s_2  = dsqr(1.0/dy[i]);
 758       sxx += s_2*dsqr(x[i]);
 759       sxy += s_2*y[i]*x[i];
 760       syy += s_2*dsqr(y[i]);
 761       sx  += s_2*x[i];
 762       sy  += s_2*y[i];
 763       w   += s_2;
 764       }
 765       sxx = sxx/w;
 766       sxy = sxy/w;
 767       syy = syy/w;
 768       sx  = sx/w;
 769       sy  = sy/w;
 770       dx2 = (sxx-sx*sx);
 771       dy2 = (syy-sy*sy);
 772       *a=(sxy-sy*sx)/dx2;
 773       *b=(sy-(*a)*sx);
 774
 775       *chi2=0;
 776       for(i=0; i<n; i++)
 777       *chi2+=dsqr((y[i]-((*a)*x[i]+(*b)))/dy[i]);
 778       *chi2 = *chi2/w;
 779
 780       *da = sqrt(*chi2/((n-2)*dx2));
 781       *db = *da*sqrt(sxx);
 782       *r  = *a*sqrt(dx2/dy2);
 783
 784       if (debug)
 785       fprintf(debug,"sx = %g, sy = %g, sxy = %g, sxx = %g, w = %g\n"
 786       "chi2 = %g, dx2 = %g\n",
 787       sx,sy,sxy,sxx,w,*chi2,dx2);
 788
 789       if (n > 2)
 790       *chi2 = sqrt(*chi2/(n-2));
 791       else
 792       *chi2 = 0;
 793       */
 794 }
 795
 796