src/gromacs/statistics/statistics.cpp

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