#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
-static int gmx_dnint(double x)
-{
- return gmx::roundToInt(x);
-}
typedef struct gmx_stats
{
return static_cast<gmx_stats_t>(stats);
}
-StatisticsStatus gmx_stats_get_npoints(gmx_stats_t gstats, int* N)
-{
- gmx_stats* stats = static_cast<gmx_stats*>(gstats);
-
- *N = stats->np;
-
- return StatisticsStatus::Ok;
-}
-
void gmx_stats_free(gmx_stats_t gstats)
{
gmx_stats* stats = static_cast<gmx_stats*>(gstats);
return StatisticsStatus::Ok;
}
-StatisticsStatus gmx_stats_get_point(gmx_stats_t gstats, real* x, real* y, real* dx, real* dy, real level)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
- int outlier;
- real rmsd, r;
-
- if ((ok = gmx_stats_get_rmsd(gstats, &rmsd)) != StatisticsStatus::Ok)
- {
- return ok;
- }
- outlier = 0;
- while ((outlier == 0) && (stats->np_c < stats->np))
- {
- r = std::abs(stats->x[stats->np_c] - stats->y[stats->np_c]);
- outlier = static_cast<int>(r > rmsd * level);
- if (outlier)
- {
- if (nullptr != x)
- {
- *x = stats->x[stats->np_c];
- }
- if (nullptr != y)
- {
- *y = stats->y[stats->np_c];
- }
- if (nullptr != dx)
- {
- *dx = stats->dx[stats->np_c];
- }
- if (nullptr != dy)
- {
- *dy = stats->dy[stats->np_c];
- }
- }
- stats->np_c++;
-
- if (outlier)
- {
- return StatisticsStatus::Ok;
- }
- }
-
- stats->np_c = 0;
-
- return StatisticsStatus::NoPoints;
-}
-
-StatisticsStatus gmx_stats_add_points(gmx_stats_t gstats, int n, real* x, real* y, real* dx, real* dy)
-{
- for (int i = 0; (i < n); i++)
- {
- StatisticsStatus ok;
- if ((ok = gmx_stats_add_point(
- gstats, x[i], y[i], (nullptr != dx) ? dx[i] : 0, (nullptr != dy) ? dy[i] : 0))
- != StatisticsStatus::Ok)
- {
- return ok;
- }
- }
- return StatisticsStatus::Ok;
-}
-
static StatisticsStatus gmx_stats_compute(gmx_stats* stats, int weight)
{
double yy, yx, xx, sx, sy, dy, chi2, chi2aa, d2;
return StatisticsStatus::Ok;
}
-StatisticsStatus gmx_stats_get_a(gmx_stats_t gstats, int weight, real* a, real* da, real* chi2, real* Rfit)
-{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, weight)) != StatisticsStatus::Ok)
- {
- return ok;
- }
- if (nullptr != a)
- {
- *a = stats->aa;
- }
- if (nullptr != da)
- {
- *da = stats->sigma_aa;
- }
- if (nullptr != chi2)
- {
- *chi2 = stats->chi2aa;
- }
- if (nullptr != Rfit)
- {
- *Rfit = stats->Rfitaa;
- }
-
- return StatisticsStatus::Ok;
-}
-
StatisticsStatus gmx_stats_get_average(gmx_stats_t gstats, real* aver)
{
gmx_stats* stats = gstats;
return StatisticsStatus::Ok;
}
-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;
-}
-
static const char* enumValueToString(StatisticsStatus enumValue)
{
constexpr gmx::EnumerationArray<StatisticsStatus, const char*> statisticsStatusNames = {
GMX_ASSERT(estats == StatisticsStatus::Ok, enumValueToString(estats));
}
-/* Old convenience functions, should be merged with the core
- statistics above. */
-StatisticsStatus lsq_y_ax(int n, real x[], real y[], real* a)
-{
- gmx_stats_t lsq = gmx_stats_init();
- real da, chi2, Rfit;
-
- gmx_stats_add_points(lsq, n, x, y, nullptr, nullptr);
- StatisticsStatus ok = gmx_stats_get_a(lsq, elsqWEIGHT_NONE, a, &da, &chi2, &Rfit);
- gmx_stats_free(lsq);
-
- return ok;
-}
-
static StatisticsStatus
low_lsq_y_ax_b(int n, const real* xr, const double* xd, real yr[], real* a, real* b, real* r, real* chi2)
{
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff