diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index ffdbf1fe38d..2a5847c57d1 100644 --- a/doc/sphinx/analysis.rst +++ b/doc/sphinx/analysis.rst @@ -739,8 +739,8 @@ In order to calculate the running mean and variance of an observable, obs=position_observable, delta_N=2) system.auto_update_accumulators.add(accumulator) system.integrator.run(10) - print(accumulator.get_mean()) - print(accumulator.get_variance()) + print(accumulator.mean()) + print(accumulator.variance()) In the example above the automatic update of the accumulator is used. However, it's also possible to manually update the accumulator by calling diff --git a/doc/tutorials/charged_system/charged_system-1.ipynb b/doc/tutorials/charged_system/charged_system-1.ipynb index da927f3a64f..485150305b3 100644 --- a/doc/tutorials/charged_system/charged_system-1.ipynb +++ b/doc/tutorials/charged_system/charged_system-1.ipynb @@ -597,7 +597,7 @@ "**Hints:**\n", "* Don't forget to clear the system before setting up the system with a new set of parameters\n", "* Don't forget to ``tune()`` the ``p3m`` instance after each change of parameters. If we reuse the p3m that was tuned before, likely the desired accuracy will not be achieved. \n", - "* Extract the radial density profile from the accumulator via ``.get_mean()``" + "* Extract the radial density profile from the accumulator via ``.mean()``" ] }, { @@ -622,7 +622,7 @@ " print('starting simulation')\n", " integrate_system(system, N_SAMPLES * STEPS_PER_SAMPLE)\n", "\n", - " run['histogram'] = radial_profile_accs[COUNTERION_TYPE].get_mean()\n", + " run['histogram'] = radial_profile_accs[COUNTERION_TYPE].mean()\n", " print('simulation for parameters {} done\\n'.format(run['params']))\n", "```" ] @@ -862,7 +862,7 @@ "rs = bin_edges[1:, 0, 0, 0]\n", "cum_hists = {}\n", "for ion_type in all_ion_types:\n", - " hist = radial_profile_accs[ion_type].get_mean()\n", + " hist = radial_profile_accs[ion_type].mean()\n", " hist = hist[:, 0, 0] * rs\n", " cum_hist = np.cumsum(hist)\n", " cum_hist /= cum_hist[-1]\n", diff --git a/doc/tutorials/lennard_jones/lennard_jones.ipynb b/doc/tutorials/lennard_jones/lennard_jones.ipynb index c16eef66c14..73c8f7154f0 100644 --- a/doc/tutorials/lennard_jones/lennard_jones.ipynb +++ b/doc/tutorials/lennard_jones/lennard_jones.ipynb @@ -956,7 +956,7 @@ }, "source": [ "```python\n", - "rdf = rdf_acc.get_mean()\n", + "rdf = rdf_acc.mean()\n", "rs = rdf_obs.bin_centers()\n", "```" ] diff --git a/samples/drude_bmimpf6.py b/samples/drude_bmimpf6.py index a29f928610b..2b2767e07a0 100644 --- a/samples/drude_bmimpf6.py +++ b/samples/drude_bmimpf6.py @@ -402,9 +402,9 @@ def combination_rule_sigma(rule, sig1, sig2): file_traj.close() - rdf_00 = acc_00.get_mean() - rdf_11 = acc_11.get_mean() - rdf_01 = acc_01.get_mean() + rdf_00 = acc_00.mean() + rdf_11 = acc_11.mean() + rdf_01 = acc_01.mean() r = obs_01.bin_centers() np.savetxt(args.path + "rdf.dat", np.c_[r, rdf_00, rdf_11, rdf_01]) print("\n-->Done") diff --git a/samples/electrophoresis.py b/samples/electrophoresis.py index d82ac8831ca..106a7e89ec3 100644 --- a/samples/electrophoresis.py +++ b/samples/electrophoresis.py @@ -214,10 +214,10 @@ def exponential(x, lp): def persistence_length_obs( acc_bond_length, acc_persistence_angles, exponential): - bond_lengths_obs = np.array(acc_bond_length.get_mean()) + bond_lengths_obs = np.array(acc_bond_length.mean()) sampling_positions_obs = np.insert( np.cumsum(bond_lengths_obs)[:-1], 0, 0.0) - cos_thetas_obs = np.array(acc_persistence_angles.get_mean()) + cos_thetas_obs = np.array(acc_persistence_angles.mean()) cos_thetas_obs = np.insert(cos_thetas_obs, 0, 1.0) opt_obs, _ = scipy.optimize.curve_fit( diff --git a/samples/lb_profile.py b/samples/lb_profile.py index b3cb65da2f7..7f4f6235dd6 100644 --- a/samples/lb_profile.py +++ b/samples/lb_profile.py @@ -70,7 +70,7 @@ system.auto_update_accumulators.add(accumulator) system.integrator.run(n_steps) -lb_fluid_profile = accumulator.get_mean() +lb_fluid_profile = accumulator.mean() def poiseuille_flow(r, R, ext_force_density): diff --git a/src/core/accumulators/MeanVarianceCalculator.cpp b/src/core/accumulators/MeanVarianceCalculator.cpp index ddc5808940f..1b8abd3f918 100644 --- a/src/core/accumulators/MeanVarianceCalculator.cpp +++ b/src/core/accumulators/MeanVarianceCalculator.cpp @@ -33,12 +33,14 @@ namespace Accumulators { void MeanVarianceCalculator::update() { m_acc(m_obs->operator()()); } -std::vector MeanVarianceCalculator::get_mean() { - return m_acc.get_mean(); +std::vector MeanVarianceCalculator::mean() { return m_acc.mean(); } + +std::vector MeanVarianceCalculator::variance() { + return m_acc.variance(); } -std::vector MeanVarianceCalculator::get_variance() { - return m_acc.get_variance(); +std::vector MeanVarianceCalculator::std_error() { + return m_acc.std_error(); } std::string MeanVarianceCalculator::get_internal_state() const { diff --git a/src/core/accumulators/MeanVarianceCalculator.hpp b/src/core/accumulators/MeanVarianceCalculator.hpp index 8fcac7b9d81..cfc9a8ae7f6 100644 --- a/src/core/accumulators/MeanVarianceCalculator.hpp +++ b/src/core/accumulators/MeanVarianceCalculator.hpp @@ -39,8 +39,9 @@ class MeanVarianceCalculator : public AccumulatorBase { : AccumulatorBase(delta_N), m_obs(obs), m_acc(obs->n_values()) {} void update() override; - std::vector get_mean(); - std::vector get_variance(); + std::vector mean(); + std::vector variance(); + std::vector std_error(); /* Partial serialization of state that is not accessible via the interface. */ std::string get_internal_state() const; diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp index 966a4fe96e0..3e783fed4b4 100644 --- a/src/core/reaction_ensemble.cpp +++ b/src/core/reaction_ensemble.cpp @@ -1611,14 +1611,13 @@ WidomInsertion::measure_excess_chemical_potential(int reaction_id) { current_reaction.accumulator_exponentials(exponential); std::pair result = std::make_pair( - -temperature * - log(current_reaction.accumulator_exponentials.get_mean()[0]), + -temperature * log(current_reaction.accumulator_exponentials.mean()[0]), std::abs(-temperature / - current_reaction.accumulator_exponentials.get_mean()[0] * + current_reaction.accumulator_exponentials.mean()[0] * current_reaction.accumulator_exponentials - .get_std_error()[0])); // excess chemical potential; error - // excess chemical potential, - // determined via error propagation + .std_error()[0])); // excess chemical potential; error + // excess chemical potential, + // determined via error propagation return result; } diff --git a/src/python/espressomd/accumulators.py b/src/python/espressomd/accumulators.py index 3d6bdc8d697..03552551a09 100644 --- a/src/python/espressomd/accumulators.py +++ b/src/python/espressomd/accumulators.py @@ -43,18 +43,25 @@ class MeanVarianceCalculator(ScriptInterfaceHelper): ) _so_creation_policy = "LOCAL" - def get_mean(self): + def mean(self): """ Returns the samples mean values of the respective observable with which the accumulator was initialized. """ - return np.array(self.call_method("get_mean")).reshape(self.shape()) + return np.array(self.call_method("mean")).reshape(self.shape()) - def get_variance(self): + def variance(self): """ Returns the samples variance for the observable. """ - return np.array(self.call_method("get_variance")).reshape(self.shape()) + return np.array(self.call_method("variance")).reshape(self.shape()) + + def std_error(self): + """ + Returns the standard error calculated from the samples variance for the observable by + assuming uncorrelated samples. + """ + return np.array(self.call_method("std_error")).reshape(self.shape()) @script_interface_register diff --git a/src/script_interface/accumulators/MeanVarianceCalculator.hpp b/src/script_interface/accumulators/MeanVarianceCalculator.hpp index 6fbc124724b..60dd6565bb6 100644 --- a/src/script_interface/accumulators/MeanVarianceCalculator.hpp +++ b/src/script_interface/accumulators/MeanVarianceCalculator.hpp @@ -61,15 +61,14 @@ class MeanVarianceCalculator : public AccumulatorBase { Variant do_call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "update") { + if (method == "update") mean_variance_calculator()->update(); - } - if (method == "get_mean") - return mean_variance_calculator()->get_mean(); - - if (method == "get_variance") - return mean_variance_calculator()->get_variance(); - + if (method == "mean") + return mean_variance_calculator()->mean(); + if (method == "variance") + return mean_variance_calculator()->variance(); + if (method == "std_error") + return mean_variance_calculator()->std_error(); return AccumulatorBase::call_method(method, parameters); } diff --git a/src/utils/include/utils/Accumulator.hpp b/src/utils/include/utils/Accumulator.hpp index b8021081097..f805282f462 100644 --- a/src/utils/include/utils/Accumulator.hpp +++ b/src/utils/include/utils/Accumulator.hpp @@ -49,9 +49,9 @@ class Accumulator { public: explicit Accumulator(std::size_t N) : m_n(0), m_acc_data(N) {} void operator()(const std::vector &); - std::vector get_mean() const; - std::vector get_variance() const; - std::vector get_std_error() const; + std::vector mean() const; + std::vector variance() const; + std::vector std_error() const; private: std::size_t m_n; @@ -88,7 +88,7 @@ inline void Accumulator::operator()(const std::vector &data) { } } -inline std::vector Accumulator::get_mean() const { +inline std::vector Accumulator::mean() const { std::vector res; std::transform( m_acc_data.begin(), m_acc_data.end(), std::back_inserter(res), @@ -96,34 +96,30 @@ inline std::vector Accumulator::get_mean() const { return res; } -inline std::vector Accumulator::get_variance() const { +inline std::vector Accumulator::variance() const { std::vector res; if (m_n == 1) { res = std::vector(m_acc_data.size(), std::numeric_limits::max()); } else { - std::transform( - m_acc_data.begin(), m_acc_data.end(), std::back_inserter(res), - [this](const AccumulatorData &acc_data) { - return acc_data.m / - (static_cast(m_n) - - 1); // numerically stable sample variance, see - // https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance - }); + std::transform(m_acc_data.begin(), m_acc_data.end(), + std::back_inserter(res), + [this](const AccumulatorData &acc_data) { + return acc_data.m / (static_cast(m_n) - 1); + }); } return res; } /** -returns the standard error of the mean of uncorrelated data. if data are -correlated the correlation time needs to be known... -*/ -inline std::vector Accumulator::get_std_error() const { - auto const variance = get_variance(); - std::vector std_error(variance.size()); - std::transform(variance.begin(), variance.end(), std_error.begin(), + * Returns the standard error of the mean assuming uncorrelated samples. + */ +inline std::vector Accumulator::std_error() const { + auto const var = variance(); + std::vector err(var.size()); + std::transform(var.begin(), var.end(), err.begin(), [this](double d) { return std::sqrt(d / m_n); }); - return std_error; + return err; } } // namespace Utils diff --git a/src/utils/tests/accumulator.cpp b/src/utils/tests/accumulator.cpp index 550fb20000c..363746b9a0a 100644 --- a/src/utils/tests/accumulator.cpp +++ b/src/utils/tests/accumulator.cpp @@ -31,14 +31,13 @@ BOOST_AUTO_TEST_CASE(accumulator) { auto test_data1 = std::vector{{0.0, 1.0, 2.0, 3.0}}; auto test_data2 = std::vector{{1.5, 3.5, 3.5, 4.5}}; acc(test_data1); - BOOST_CHECK(acc.get_mean() == test_data1); - BOOST_CHECK(acc.get_variance() == + BOOST_CHECK(acc.mean() == test_data1); + BOOST_CHECK(acc.variance() == std::vector(4, std::numeric_limits::max())); acc(test_data2); + BOOST_CHECK((acc.mean() == std::vector{{0.75, 2.25, 2.75, 3.75}})); BOOST_CHECK( - (acc.get_mean() == std::vector{{0.75, 2.25, 2.75, 3.75}})); - BOOST_CHECK((acc.get_variance() == - std::vector{{1.125, 3.125, 1.125, 1.125}})); + (acc.variance() == std::vector{{1.125, 3.125, 1.125, 1.125}})); BOOST_CHECK( - (acc.get_std_error() == std::vector{{0.75, 1.25, 0.75, 0.75}})); + (acc.std_error() == std::vector{{0.75, 1.25, 0.75, 0.75}})); } diff --git a/testsuite/python/accumulator.py b/testsuite/python/accumulator.py index c0ac2a6bf9b..db5332c09be 100644 --- a/testsuite/python/accumulator.py +++ b/testsuite/python/accumulator.py @@ -21,9 +21,10 @@ Testmodule for the observable accumulator. """ -import sys import unittest as ut + import numpy as np + import espressomd import espressomd.observables import espressomd.accumulators @@ -62,17 +63,17 @@ def test_accumulator(self): self.system.part[:].pos = pos self.system.integrator.run(1) self.assertEqual(self.pos_obs, self.pos_obs_acc.get_params()['obs']) - acc_mean = self.pos_obs_acc.get_mean() np.testing.assert_allclose( - acc_mean, + self.pos_obs_acc.mean(), np.mean(self.positions, axis=0), atol=1e-4) + variance = np.var(self.positions, axis=0, ddof=1) + np.testing.assert_allclose( + self.pos_obs_acc.variance(), + variance, atol=1e-4) np.testing.assert_allclose( - self.pos_obs_acc.get_variance(), - np.var(self.positions, axis=0, ddof=1), atol=1e-4) + self.pos_obs_acc.std_error(), + np.sqrt(variance / len(self.positions)), atol=1e-4) if __name__ == "__main__": - suite = ut.TestSuite() - suite.addTests(ut.TestLoader().loadTestsFromTestCase(AccumulatorTest)) - result = ut.TextTestRunner(verbosity=4).run(suite) - sys.exit(not result.wasSuccessful()) + ut.main() diff --git a/testsuite/python/lb_poiseuille_cylinder.py b/testsuite/python/lb_poiseuille_cylinder.py index 88790bd992c..f0dee7984a7 100644 --- a/testsuite/python/lb_poiseuille_cylinder.py +++ b/testsuite/python/lb_poiseuille_cylinder.py @@ -159,11 +159,7 @@ def check_observable(self): self.prepare_obs() # gather some statistics for the observable accumulator self.system.integrator.run(1) - obs_result = np.array( - self.accumulator.get_mean()).reshape(OBS_PARAMS['n_r_bins'], - OBS_PARAMS['n_phi_bins'], - OBS_PARAMS['n_z_bins'], - 3) + obs_result = self.accumulator.mean() x = np.linspace( OBS_PARAMS['min_r'], OBS_PARAMS['max_r'], diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py index b90fc388415..6eb604fdda8 100644 --- a/testsuite/python/test_checkpoint.py +++ b/testsuite/python/test_checkpoint.py @@ -292,13 +292,13 @@ def test_virtual_sites(self): def test_mean_variance_calculator(self): np.testing.assert_array_equal( - acc_mean_variance.get_mean(), + acc_mean_variance.mean(), np.array([[1.0, 1.5, 2.0], [1.0, 1.0, 2.0]])) np.testing.assert_array_equal( - acc_mean_variance.get_variance(), + acc_mean_variance.variance(), np.array([[0., 0.5, 2.], [0., 0., 0.]])) np.testing.assert_array_equal( - system.auto_update_accumulators[0].get_variance(), + system.auto_update_accumulators[0].variance(), np.array([[0., 0.5, 2.], [0., 0., 0.]])) def test_time_series(self):