-StatisticsStatus gmx_stats_get_sigma(gmx_stats_t gstats, real* sigma)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
- {
- return ok;
- }
-
- *sigma = stats->sigma_aver;
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_get_error(gmx_stats_t gstats, real* error)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
- {
- return ok;
- }
-
- *error = stats->error;
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_get_corr_coeff(gmx_stats_t gstats, real* R)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
- {
- return ok;
- }
-
- *R = stats->Rdata;
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_get_rmsd(gmx_stats_t gstats, real* rmsd)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
- {
- return ok;
- }
-
- *rmsd = stats->rmsd;
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_dump_xy(gmx_stats_t gstats, FILE* fp)
-{
- gmx_stats* stats = gstats;
-
- for (int i = 0; (i < stats->np); i++)
- {
- fprintf(fp, "%12g %12g %12g %12g\n", stats->x[i], stats->y[i], stats->dx[i], stats->dy[i]);
- }
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_remove_outliers(gmx_stats_t gstats, double level)
-{
- gmx_stats* stats = gstats;
- int iter = 1, done = 0;
- StatisticsStatus ok;
- real rmsd, r;
-
- while ((stats->np >= 10) && !done)
- {
- if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != StatisticsStatus::Ok)
- {
- return ok;
- }
- done = 1;
- for (int i = 0; (i < stats->np);)
- {
- r = std::abs(stats->x[i] - stats->y[i]);
- if (r > level * rmsd)
- {
- fprintf(stderr,
- "Removing outlier, iter = %d, rmsd = %g, x = %g, y = %g\n",
- iter,
- rmsd,
- stats->x[i],
- stats->y[i]);
- if (i < stats->np - 1)
- {
- stats->x[i] = stats->x[stats->np - 1];
- stats->y[i] = stats->y[stats->np - 1];
- stats->dx[i] = stats->dx[stats->np - 1];
- stats->dy[i] = stats->dy[stats->np - 1];
- }
- stats->np--;
- done = 0;
- }
- else
- {
- i++;
- }
- }
- iter++;
- }
-
- return StatisticsStatus::Ok;
-}
-
-StatisticsStatus gmx_stats_make_histogram(gmx_stats_t gstats,
- real binwidth,
- int* nb,
- int ehisto,
- int normalized,
- real** x,
- real** y)
-{
- gmx_stats* stats = gstats;
- int index = 0, nbins = *nb, *nindex;
- double minx, maxx, maxy, miny, delta, dd, minh;
-
- if (((binwidth <= 0) && (nbins <= 0)) || ((binwidth > 0) && (nbins > 0)))
- {
- return StatisticsStatus::InvalidInput;
- }
- if (stats->np <= 2)
- {
- return StatisticsStatus::NoPoints;
- }
- minx = maxx = stats->x[0];
- miny = maxy = stats->y[0];
- for (int i = 1; (i < stats->np); i++)
- {
- miny = (stats->y[i] < miny) ? stats->y[i] : miny;
- maxy = (stats->y[i] > maxy) ? stats->y[i] : maxy;
- minx = (stats->x[i] < minx) ? stats->x[i] : minx;
- maxx = (stats->x[i] > maxx) ? stats->x[i] : maxx;
- }
- if (ehisto == ehistoX)
- {
- delta = maxx - minx;
- minh = minx;
- }
- else if (ehisto == ehistoY)
- {
- delta = maxy - miny;
- minh = miny;
- }
- else
- {
- return StatisticsStatus::InvalidInput;
- }
-
- if (binwidth == 0)
- {
- binwidth = (delta) / nbins;
- }
- else
- {
- nbins = gmx_dnint((delta) / binwidth + 0.5);
- }
- snew(*x, nbins);
- snew(nindex, nbins);
- for (int i = 0; (i < nbins); i++)
- {
- (*x)[i] = minh + binwidth * (i + 0.5);
- }
- if (normalized == 0)
- {
- dd = 1;
- }
- else
- {
- dd = 1.0 / (binwidth * stats->np);
- }
-
- snew(*y, nbins);
- for (int i = 0; (i < stats->np); i++)
- {
- if (ehisto == ehistoY)
- {
- index = static_cast<int>((stats->y[i] - miny) / binwidth);
- }
- else if (ehisto == ehistoX)
- {
- index = static_cast<int>((stats->x[i] - minx) / binwidth);
- }
- if (index < 0)
- {
- index = 0;
- }
- if (index > nbins - 1)
- {
- index = nbins - 1;
- }
- (*y)[index] += dd;
- nindex[index]++;
- }
- if (*nb == 0)
- {
- *nb = nbins;
- }
- for (int i = 0; (i < nbins); i++)
- {
- if (nindex[i] > 0)
- {
- (*y)[i] /= nindex[i];
- }
- }
-
- sfree(nindex);
-
- return StatisticsStatus::Ok;
-}
-