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