if (fnadip)
{
- real aver;
- gmx_stats_get_average(muframelsq, &aver);
- fprintf(adip, "%10g %f \n", t, aver);
+ fprintf(adip, "%10g %f \n", t, gmx_stats_get_average(muframelsq));
}
/*if (dipole)
printf("%f %f\n", norm(dipole[0]), norm(dipole[1]));
}
if (!bMU)
{
- real aver, sigma, error;
-
- gmx_stats_get_ase(mulsq, &aver, &sigma, &error);
+ auto [aver, sigma, error] = gmx_stats_get_ase(mulsq);
printf("\nDipole moment (Debye)\n");
printf("---------------------\n");
printf("Average = %8.4f Std. Dev. = %8.4f Error = %8.4f\n", aver, sigma, error);
if (bQuad)
{
- rvec a, s, e;
- for (m = 0; (m < DIM); m++)
- {
- gmx_stats_get_ase(mulsq, &(a[m]), &(s[m]), &(e[m]));
- }
+ auto [averageXX, sigmaXX, errorXX] = gmx_stats_get_ase(Qlsq[XX]);
+ auto [averageYY, sigmaYY, errorYY] = gmx_stats_get_ase(Qlsq[YY]);
+ auto [averageZZ, sigmaZZ, errorZZ] = gmx_stats_get_ase(Qlsq[ZZ]);
printf("\nQuadrupole moment (Debye-Ang)\n");
printf("-----------------------------\n");
- printf("Averages = %8.4f %8.4f %8.4f\n", a[XX], a[YY], a[ZZ]);
- printf("Std. Dev. = %8.4f %8.4f %8.4f\n", s[XX], s[YY], s[ZZ]);
- printf("Error = %8.4f %8.4f %8.4f\n", e[XX], e[YY], e[ZZ]);
+ printf("Averages = %8.4f %8.4f %8.4f\n", averageXX, averageYY, averageZZ);
+ printf("Std. Dev. = %8.4f %8.4f %8.4f\n", sigmaXX, sigmaYY, sigmaZZ);
+ printf("Error = %8.4f %8.4f %8.4f\n", errorXX, errorYY, errorZZ);
}
printf("\n");
}
*/
#include "gmxpre.h"
-#include "gromacs/utility/enumerationhelpers.h"
#include "statistics.h"
#include <cmath>
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
-#include "gromacs/utility/enumerationhelpers.h"
-#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
sfree(stats);
}
-StatisticsStatus gmx_stats_add_point(gmx_stats_t gstats, double x, double y, double dx, double dy)
+void gmx_stats_add_point(gmx_stats_t gstats, double x, double y, double dx, double dy)
{
gmx_stats* stats = gstats;
stats->dy[stats->np] = dy;
stats->np++;
stats->computed = 0;
-
- return StatisticsStatus::Ok;
}
-static StatisticsStatus gmx_stats_compute(gmx_stats* stats, int weight)
+static void gmx_stats_compute(gmx_stats* stats, int weight)
{
double yy, yx, xx, sx, sy, dy, chi2, chi2aa, d2;
double ssxx, ssyy, ssxy;
if (stats->computed == 0)
{
- if (N < 1)
- {
- return StatisticsStatus::NoPoints;
- }
+ GMX_RELEASE_ASSERT(N >= 1, "Must have points to work on");
xx = xx_nw = 0;
yy = yy_nw = 0;
stats->computed = 1;
}
-
- return StatisticsStatus::Ok;
}
-StatisticsStatus
-gmx_stats_get_ab(gmx_stats_t gstats, int weight, real* a, real* b, real* da, real* db, real* chi2, real* Rfit)
+void gmx_stats_get_ab(gmx_stats_t gstats, int weight, real* a, real* b, real* da, real* db, real* chi2, real* Rfit)
{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
+ gmx_stats* stats = gstats;
- if ((ok = gmx_stats_compute(stats, weight)) != StatisticsStatus::Ok)
- {
- return ok;
- }
+ gmx_stats_compute(stats, weight);
if (nullptr != a)
{
*a = stats->a;
{
*Rfit = stats->Rfit;
}
-
- return StatisticsStatus::Ok;
}
-StatisticsStatus gmx_stats_get_average(gmx_stats_t gstats, real* aver)
+real gmx_stats_get_average(gmx_stats_t gstats)
{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
+ gmx_stats* stats = gstats;
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
+ if (gstats->np < 1)
{
- return ok;
+ GMX_THROW(gmx::InconsistentInputError("No points to average"));
}
+ gmx_stats_compute(stats, elsqWEIGHT_NONE);
- *aver = stats->aver;
-
- return StatisticsStatus::Ok;
+ return stats->aver;
}
-StatisticsStatus gmx_stats_get_ase(gmx_stats_t gstats, real* aver, real* sigma, real* error)
+std::tuple<real, real, real> gmx_stats_get_ase(gmx_stats_t gstats)
{
- gmx_stats* stats = gstats;
- StatisticsStatus ok;
-
- if ((ok = gmx_stats_compute(stats, elsqWEIGHT_NONE)) != StatisticsStatus::Ok)
- {
- return ok;
- }
+ gmx_stats* stats = gstats;
- if (nullptr != aver)
- {
- *aver = stats->aver;
- }
- if (nullptr != sigma)
+ if (gstats->np < 1)
{
- *sigma = stats->sigma_aver;
+ GMX_THROW(gmx::InconsistentInputError("No points to average"));
}
- if (nullptr != error)
- {
- *error = stats->error;
- }
-
- return StatisticsStatus::Ok;
-}
-
-static const char* enumValueToString(StatisticsStatus enumValue)
-{
- constexpr gmx::EnumerationArray<StatisticsStatus, const char*> statisticsStatusNames = {
- "All well in STATS land", "No points"
- };
- return statisticsStatusNames[enumValue];
-}
+ gmx_stats_compute(stats, elsqWEIGHT_NONE);
-void gmx_stats_message([[maybe_unused]] StatisticsStatus estats)
-{
- GMX_ASSERT(estats == StatisticsStatus::Ok, enumValueToString(estats));
+ return std::make_tuple(stats->aver, stats->sigma_aver, stats->error);
}
-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)
+// When using GMX_DOUBLE=OFF, some callers want to analyse x values
+// that are already computed in double precision. So we need to
+// compile two versions, so that the promotion to double is used when
+// needed.
+template<typename RealT>
+void low_lsq_y_ax_b(int n, const RealT* xr, real yr[], real* a, real* b, real* r, real* chi2)
{
- gmx_stats_t lsq = gmx_stats_init();
- StatisticsStatus ok;
-
+ gmx_stats_t lsq = gmx_stats_init();
for (int i = 0; (i < n); i++)
{
- double pt;
-
- if (xd != nullptr)
- {
- pt = xd[i];
- }
- else if (xr != nullptr)
- {
- pt = xr[i];
- }
- else
- {
- gmx_incons("Either xd or xr has to be non-NULL in low_lsq_y_ax_b()");
- }
-
- if ((ok = gmx_stats_add_point(lsq, pt, yr[i], 0, 0)) != StatisticsStatus::Ok)
- {
- gmx_stats_free(lsq);
- return ok;
- }
+ gmx_stats_add_point(lsq, double(xr[i]), yr[i], 0, 0);
}
- ok = gmx_stats_get_ab(lsq, elsqWEIGHT_NONE, a, b, nullptr, nullptr, chi2, r);
+ gmx_stats_get_ab(lsq, elsqWEIGHT_NONE, a, b, nullptr, nullptr, chi2, r);
gmx_stats_free(lsq);
-
- return ok;
}
-StatisticsStatus lsq_y_ax_b(int n, real x[], real y[], real* a, real* b, real* r, real* chi2)
+void lsq_y_ax_b(int n, real x[], real y[], real* a, real* b, real* r, real* chi2)
{
- return low_lsq_y_ax_b(n, x, nullptr, y, a, b, r, chi2);
+ low_lsq_y_ax_b(n, x, y, a, b, r, chi2);
}
-StatisticsStatus lsq_y_ax_b_xdouble(int n, double x[], real y[], real* a, real* b, real* r, real* chi2)
+void lsq_y_ax_b_xdouble(int n, double x[], real y[], real* a, real* b, real* r, real* chi2)
{
- return low_lsq_y_ax_b(n, nullptr, x, y, a, b, r, chi2);
+ low_lsq_y_ax_b(n, x, y, a, b, r, chi2);
}
-StatisticsStatus
-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)
+void 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)
{
- gmx_stats_t lsq = gmx_stats_init();
- StatisticsStatus ok;
+ if (n < 1)
+ {
+ GMX_THROW(gmx::InconsistentInputError("No points to fit"));
+ }
+ gmx_stats_t lsq = gmx_stats_init();
for (int i = 0; (i < n); i++)
{
- ok = gmx_stats_add_point(lsq, x[i], y[i], 0, dy[i]);
- if (ok != StatisticsStatus::Ok)
- {
- gmx_stats_free(lsq);
- return ok;
- }
+ gmx_stats_add_point(lsq, x[i], y[i], 0, dy[i]);
}
- ok = gmx_stats_get_ab(lsq, elsqWEIGHT_Y, a, b, da, db, chi2, r);
+ gmx_stats_get_ab(lsq, elsqWEIGHT_Y, a, b, da, db, chi2, r);
gmx_stats_free(lsq);
-
- return ok;
}
#include <cstdio>
+#include <tuple>
+
#include "gromacs/utility/real.h"
//! Abstract container type
typedef struct gmx_stats* gmx_stats_t;
-//! Error codes returned by the routines
-enum class StatisticsStatus : int
-{
- Ok,
- NoPoints,
- Count
-};
-
//! Enum for statistical weights
enum
{
* \param[in] y The y value
* \param[in] dx The error in the x value
* \param[in] dy The error in the y value
- * \return error code
*/
-StatisticsStatus gmx_stats_add_point(gmx_stats_t stats, double x, double y, double dx, double dy);
+void gmx_stats_add_point(gmx_stats_t stats, double x, double y, double dx, double dy);
/*! \brief
* Fit the data to y = ax + b, possibly weighted, if uncertainties
* \param[out] db sigma in b
* \param[out] chi2 normalized quality of fit
* \param[out] Rfit correlation coefficient
- * \return error code
*/
-StatisticsStatus
-gmx_stats_get_ab(gmx_stats_t stats, int weight, real* a, real* b, real* da, real* db, real* chi2, real* Rfit);
+void gmx_stats_get_ab(gmx_stats_t stats, int weight, real* a, real* b, real* da, real* db, real* chi2, real* Rfit);
/*! \brief
* Computes and returns the average value.
* \param[in] stats The data structure
* \param[out] aver Average value
- * \return error code
+ * \return Average value
+ * \throws InconsistentInputError if given no points to average
*/
-StatisticsStatus gmx_stats_get_average(gmx_stats_t stats, real* aver);
+real gmx_stats_get_average(gmx_stats_t stats);
/*! \brief
* Pointers may be null, in which case no assignment will be done.
* \param[in] stats The data structure
- * \param[out] aver Average value
- * \param[out] sigma Standard deviation
- * \param[out] error Standard error
- * \return error code
+ * \return Tuple of (average value, its standard deviation, its standard error)
+ * \throws InconsistentInputError if given no points to analyze
*/
-StatisticsStatus gmx_stats_get_ase(gmx_stats_t stats, real* aver, real* sigma, real* error);
-
-/*! \brief Assert that statistics are OK
- *
- * \param[in] estats error code
- */
-void gmx_stats_message([[maybe_unused]] StatisticsStatus estats);
+std::tuple<real, real, real> gmx_stats_get_ase(gmx_stats_t gstats);
/****************************************************
* Some statistics utilities for convenience: useful when a complete data
* \param[out] b intercept
* \param[out] r correlation coefficient
* \param[out] chi2 quality of fit
- * \return error code
+ *
+ * \throws InconsistentInputError if given no points to fit
*/
-StatisticsStatus lsq_y_ax_b(int n, real x[], real y[], real* a, real* b, real* r, real* chi2);
+void lsq_y_ax_b(int n, real x[], real y[], real* a, real* b, real* r, real* chi2);
/*! \copydoc lsq_y_ax_b
+ * Suits cases where x is already always computed in double precision
+ * even in a mixed-precision build configuration.
*/
-StatisticsStatus lsq_y_ax_b_xdouble(int n, double x[], real y[], real* a, real* b, real* r, real* chi2);
+void lsq_y_ax_b_xdouble(int n, double x[], real y[], real* a, real* b, real* r, real* chi2);
/*! \brief
* Fit a straight line y=ax+b thru the n data points x, y.
* \param[out] db error in intercept
* \param[out] r correlation coefficient
* \param[out] chi2 quality of fit
- * \return error code
+ *
+ * \throws InconsistentInputError if given no points to fit
*/
-StatisticsStatus
-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);
+void 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);
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff