From 72daefa7043f0ef1a39617ef91525543d0a2a95a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 30 Mar 2020 17:39:34 +0200 Subject: [PATCH 1/8] Refactor Profile observables Make CylinderProfile a subclass of Observable, cleanup include statements, add a Base member in the script interface classes. --- .../CylindricalLBProfileObservable.hpp | 11 ++++----- .../CylindricalPidProfileObservable.hpp | 11 ++++----- ...e.hpp => CylindricalProfileObservable.hpp} | 23 +++++++++++++------ src/core/observables/ProfileObservable.hpp | 6 ++++- .../CylindricalLBProfileObservable.hpp | 4 ++++ .../CylindricalPidProfileObservable.hpp | 4 ++++ .../observables/LBProfileObservable.hpp | 3 +++ .../observables/PidProfileObservable.hpp | 3 +++ 8 files changed, 45 insertions(+), 20 deletions(-) rename src/core/observables/{CylindricalProfile.hpp => CylindricalProfileObservable.hpp} (69%) diff --git a/src/core/observables/CylindricalLBProfileObservable.hpp b/src/core/observables/CylindricalLBProfileObservable.hpp index db64b084534..241ebc7bbde 100644 --- a/src/core/observables/CylindricalLBProfileObservable.hpp +++ b/src/core/observables/CylindricalLBProfileObservable.hpp @@ -19,8 +19,7 @@ #ifndef OBSERVABLES_CYLINDRICALLBPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALLBPROFILEOBSERVABLE_HPP -#include "CylindricalProfile.hpp" -#include "Observable.hpp" +#include "CylindricalProfileObservable.hpp" #include #include @@ -31,16 +30,16 @@ using Utils::Vector3d; namespace Observables { -class CylindricalLBProfileObservable : public Observable, - public CylindricalProfile { +class CylindricalLBProfileObservable : public CylindricalProfileObservable { public: CylindricalLBProfileObservable(Vector3d const ¢er, Vector3d const &axis, int n_r_bins, int n_phi_bins, int n_z_bins, double min_r, double min_phi, double min_z, double max_r, double max_phi, double max_z, double sampling_density) - : CylindricalProfile(center, axis, min_r, max_r, min_phi, max_phi, min_z, - max_z, n_r_bins, n_phi_bins, n_z_bins), + : CylindricalProfileObservable(center, axis, min_r, max_r, min_phi, + max_phi, min_z, max_z, n_r_bins, + n_phi_bins, n_z_bins), sampling_density(sampling_density) { calculate_sampling_positions(); } diff --git a/src/core/observables/CylindricalPidProfileObservable.hpp b/src/core/observables/CylindricalPidProfileObservable.hpp index 60932aef263..ffc576a54bb 100644 --- a/src/core/observables/CylindricalPidProfileObservable.hpp +++ b/src/core/observables/CylindricalPidProfileObservable.hpp @@ -19,15 +19,13 @@ #ifndef OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP -#include - -#include "CylindricalProfile.hpp" +#include "CylindricalProfileObservable.hpp" #include "PidObservable.hpp" namespace Observables { class CylindricalPidProfileObservable : public PidObservable, - public CylindricalProfile { + public CylindricalProfileObservable { public: CylindricalPidProfileObservable(std::vector const &ids, Utils::Vector3d const ¢er, @@ -36,8 +34,9 @@ class CylindricalPidProfileObservable : public PidObservable, double min_phi, double min_z, double max_r, double max_phi, double max_z) : PidObservable(ids), - CylindricalProfile(center, axis, min_r, max_r, min_phi, max_phi, min_z, - max_z, n_r_bins, n_phi_bins, n_z_bins) {} + CylindricalProfileObservable(center, axis, min_r, max_r, min_phi, + max_phi, min_z, max_z, n_r_bins, + n_phi_bins, n_z_bins) {} }; } // Namespace Observables diff --git a/src/core/observables/CylindricalProfile.hpp b/src/core/observables/CylindricalProfileObservable.hpp similarity index 69% rename from src/core/observables/CylindricalProfile.hpp rename to src/core/observables/CylindricalProfileObservable.hpp index a87e5e2506f..5b8e458d3dd 100644 --- a/src/core/observables/CylindricalProfile.hpp +++ b/src/core/observables/CylindricalProfileObservable.hpp @@ -19,17 +19,22 @@ #ifndef OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP -#include +#include + +#include "Observable.hpp" -#include +#include namespace Observables { -class CylindricalProfile { + +/** Cylindrical profile observable */ +class CylindricalProfileObservable : virtual public Observable { public: - CylindricalProfile(Utils::Vector3d const ¢er, Utils::Vector3d const &axis, - double min_r, double max_r, double min_phi, double max_phi, - double min_z, double max_z, int n_r_bins, int n_phi_bins, - int n_z_bins) + CylindricalProfileObservable(Utils::Vector3d const ¢er, + Utils::Vector3d const &axis, double min_r, + double max_r, double min_phi, double max_phi, + double min_z, double max_z, int n_r_bins, + int n_phi_bins, int n_z_bins) : center(center), axis(axis), min_r(min_r), max_r(max_r), min_phi(min_phi), max_phi(max_phi), min_z(min_z), max_z(max_z), n_r_bins(static_cast(n_r_bins)), @@ -37,11 +42,15 @@ class CylindricalProfile { n_z_bins(static_cast(n_z_bins)){}; Utils::Vector3d center; Utils::Vector3d axis; + // Range of the profile edges. double min_r, max_r; double min_phi, max_phi; double min_z, max_z; // Number of bins for each coordinate. size_t n_r_bins, n_phi_bins, n_z_bins; + std::vector shape() const override { + return {n_r_bins, n_phi_bins, n_z_bins}; + } }; } // Namespace Observables diff --git a/src/core/observables/ProfileObservable.hpp b/src/core/observables/ProfileObservable.hpp index f9626c836fc..26eac949aaa 100644 --- a/src/core/observables/ProfileObservable.hpp +++ b/src/core/observables/ProfileObservable.hpp @@ -19,11 +19,13 @@ #ifndef OBSERVABLES_PROFILEOBSERVABLE_HPP #define OBSERVABLES_PROFILEOBSERVABLE_HPP +#include + #include "Observable.hpp" namespace Observables { -// Observable which acts on a given list of particle ids +/** Cartesian profile observable */ class ProfileObservable : virtual public Observable { public: ProfileObservable(double min_x, double max_x, double min_y, double max_y, @@ -33,9 +35,11 @@ class ProfileObservable : virtual public Observable { max_z(max_z), n_x_bins(static_cast(n_x_bins)), n_y_bins(static_cast(n_y_bins)), n_z_bins(static_cast(n_z_bins)) {} + // Range of the profile edges. double min_x, max_x; double min_y, max_y; double min_z, max_z; + // Number of bins for each coordinate. size_t n_x_bins, n_y_bins, n_z_bins; std::vector shape() const override { return {n_x_bins, n_y_bins, n_z_bins}; diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp index 95847f79147..7dd568af562 100644 --- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp @@ -37,10 +37,14 @@ template class CylindricalLBProfileObservable : public AutoParameters, Observable> { + using Base = + AutoParameters, Observable>; + public: static_assert(std::is_base_of<::Observables::CylindricalLBProfileObservable, CoreCylLBObs>::value, ""); + using Base::Base; CylindricalLBProfileObservable() { this->add_parameters({ {"center", diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp index 4ae01fbef83..1a80d76839c 100644 --- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp @@ -37,10 +37,14 @@ template class CylindricalPidProfileObservable : public AutoParameters, Observable> { + using Base = + AutoParameters, Observable>; + public: static_assert(std::is_base_of<::Observables::CylindricalPidProfileObservable, CoreObs>::value, ""); + using Base::Base; CylindricalPidProfileObservable() { this->add_parameters({ {"ids", diff --git a/src/script_interface/observables/LBProfileObservable.hpp b/src/script_interface/observables/LBProfileObservable.hpp index 626d1deb25b..01ad5cc4198 100644 --- a/src/script_interface/observables/LBProfileObservable.hpp +++ b/src/script_interface/observables/LBProfileObservable.hpp @@ -35,10 +35,13 @@ namespace Observables { template class LBProfileObservable : public AutoParameters, Observable> { + using Base = AutoParameters, Observable>; + public: static_assert( std::is_base_of<::Observables::LBProfileObservable, CoreLBObs>::value, ""); + using Base::Base; LBProfileObservable() { this->add_parameters( {{"n_x_bins", diff --git a/src/script_interface/observables/PidProfileObservable.hpp b/src/script_interface/observables/PidProfileObservable.hpp index 499bf72c722..be9a2df7af0 100644 --- a/src/script_interface/observables/PidProfileObservable.hpp +++ b/src/script_interface/observables/PidProfileObservable.hpp @@ -36,7 +36,10 @@ namespace Observables { template class PidProfileObservable : public AutoParameters, Observable> { + using Base = AutoParameters, Observable>; + public: + using Base::Base; PidProfileObservable() { this->add_parameters( {{"ids", From 1b0653c171a3e21e2ebd61e64274e7b3552bc8f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 30 Mar 2020 19:36:05 +0200 Subject: [PATCH 2/8] Calculate histogram edges in the core --- .../CylindricalProfileObservable.hpp | 26 +++++++++++++++++++ src/core/observables/ProfileObservable.hpp | 26 +++++++++++++++++++ src/python/espressomd/observables.py | 2 +- .../CylindricalLBProfileObservable.hpp | 16 +++++++----- .../CylindricalPidProfileObservable.hpp | 14 ++++++++++ .../observables/LBProfileObservable.hpp | 16 +++++++----- .../observables/PidProfileObservable.hpp | 14 ++++++++++ testsuite/python/observable_cylindrical.py | 13 ++++++---- testsuite/python/observable_cylindricalLB.py | 5 +++- testsuite/python/observable_profile.py | 19 +++++++++++++- testsuite/python/observable_profileLB.py | 19 ++++++++++++++ 11 files changed, 148 insertions(+), 22 deletions(-) diff --git a/src/core/observables/CylindricalProfileObservable.hpp b/src/core/observables/CylindricalProfileObservable.hpp index 5b8e458d3dd..5f4459433e3 100644 --- a/src/core/observables/CylindricalProfileObservable.hpp +++ b/src/core/observables/CylindricalProfileObservable.hpp @@ -19,6 +19,7 @@ #ifndef OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP +#include #include #include "Observable.hpp" @@ -48,9 +49,34 @@ class CylindricalProfileObservable : virtual public Observable { double min_z, max_z; // Number of bins for each coordinate. size_t n_r_bins, n_phi_bins, n_z_bins; + std::vector shape() const override { return {n_r_bins, n_phi_bins, n_z_bins}; } + + /** Calculate the bin edges for each dimension */ + std::vector> edges() { + std::vector> profile_edges(3); + profile_edges[0] = std::vector(n_r_bins + 1); + profile_edges[1] = std::vector(n_phi_bins + 1); + profile_edges[2] = std::vector(n_z_bins + 1); + std::generate(profile_edges[0].begin(), profile_edges[0].end(), + [n = 0, start = min_r, + bin_size = (max_r - min_r) / n_r_bins]() mutable { + return start + bin_size * n++; + }); + std::generate(profile_edges[1].begin(), profile_edges[1].end(), + [n = 0, start = min_phi, + bin_size = (max_phi - min_phi) / n_phi_bins]() mutable { + return start + bin_size * n++; + }); + std::generate(profile_edges[2].begin(), profile_edges[2].end(), + [n = 0, start = min_z, + bin_size = (max_z - min_z) / n_z_bins]() mutable { + return start + bin_size * n++; + }); + return profile_edges; + } }; } // Namespace Observables diff --git a/src/core/observables/ProfileObservable.hpp b/src/core/observables/ProfileObservable.hpp index 26eac949aaa..9e0b0d2974a 100644 --- a/src/core/observables/ProfileObservable.hpp +++ b/src/core/observables/ProfileObservable.hpp @@ -19,6 +19,7 @@ #ifndef OBSERVABLES_PROFILEOBSERVABLE_HPP #define OBSERVABLES_PROFILEOBSERVABLE_HPP +#include #include #include "Observable.hpp" @@ -41,9 +42,34 @@ class ProfileObservable : virtual public Observable { double min_z, max_z; // Number of bins for each coordinate. size_t n_x_bins, n_y_bins, n_z_bins; + std::vector shape() const override { return {n_x_bins, n_y_bins, n_z_bins}; } + + /** Calculate the bin edges for each dimension */ + std::vector> edges() { + std::vector> profile_edges(3); + profile_edges[0] = std::vector(n_x_bins + 1); + profile_edges[1] = std::vector(n_y_bins + 1); + profile_edges[2] = std::vector(n_z_bins + 1); + std::generate(profile_edges[0].begin(), profile_edges[0].end(), + [n = 0, start = min_x, + bin_size = (max_x - min_x) / n_x_bins]() mutable { + return start + bin_size * n++; + }); + std::generate(profile_edges[1].begin(), profile_edges[1].end(), + [n = 0, start = min_y, + bin_size = (max_y - min_y) / n_y_bins]() mutable { + return start + bin_size * n++; + }); + std::generate(profile_edges[2].begin(), profile_edges[2].end(), + [n = 0, start = min_z, + bin_size = (max_z - min_z) / n_z_bins]() mutable { + return start + bin_size * n++; + }); + return profile_edges; + } }; } // Namespace Observables diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index 9dbc23b460d..f3e8534760e 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -21,7 +21,7 @@ @script_interface_register class Observable(ScriptInterfaceHelper): _so_name = "Observables::Observable" - _so_bind_methods = ("shape",) + _so_bind_methods = ("shape", "edges") _so_creation_policy = "LOCAL" def calculate(self): diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp index 7dd568af562..88bbfa3e7d4 100644 --- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp @@ -137,14 +137,16 @@ class CylindricalLBProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "calculate") { - return cylindrical_profile_observable()->operator()(); + if (method == "edges") { + auto const core_edges = cylindrical_profile_observable()->edges(); + std::vector variant_edges; + for (auto &edge : core_edges) { + variant_edges.push_back(Variant{edge}); + } + return variant_edges; + } else { + return Base::call_method(method, parameters); } - if (method == "shape") { - auto const shape = cylindrical_profile_observable()->shape(); - return std::vector{shape.begin(), shape.end()}; - } - return {}; } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp index 1a80d76839c..93510fe98a2 100644 --- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp @@ -135,6 +135,20 @@ class CylindricalPidProfileObservable "max_z"); } + Variant call_method(std::string const &method, + VariantMap const ¶meters) override { + if (method == "edges") { + auto const core_edges = cylindrical_pid_profile_observable()->edges(); + std::vector variant_edges; + for (auto &edge : core_edges) { + variant_edges.push_back(Variant{edge}); + } + return variant_edges; + } else { + return Base::call_method(method, parameters); + } + } + std::shared_ptr<::Observables::Observable> observable() const override { return m_observable; } diff --git a/src/script_interface/observables/LBProfileObservable.hpp b/src/script_interface/observables/LBProfileObservable.hpp index 01ad5cc4198..aa8105323ab 100644 --- a/src/script_interface/observables/LBProfileObservable.hpp +++ b/src/script_interface/observables/LBProfileObservable.hpp @@ -148,14 +148,16 @@ class LBProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { - if (method == "calculate") { - return profile_observable()->operator()(); + if (method == "edges") { + auto const core_edges = profile_observable()->edges(); + std::vector variant_edges; + for (auto &edge : core_edges) { + variant_edges.push_back(Variant{edge}); + } + return variant_edges; + } else { + return Base::call_method(method, parameters); } - if (method == "shape") { - auto const shape = profile_observable()->shape(); - return std::vector{shape.begin(), shape.end()}; - } - return {}; } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/PidProfileObservable.hpp b/src/script_interface/observables/PidProfileObservable.hpp index be9a2df7af0..60a4b10d81a 100644 --- a/src/script_interface/observables/PidProfileObservable.hpp +++ b/src/script_interface/observables/PidProfileObservable.hpp @@ -111,6 +111,20 @@ class PidProfileObservable "min_z", "max_x", "max_y", "max_z"); } + Variant call_method(std::string const &method, + VariantMap const ¶meters) override { + if (method == "edges") { + auto const core_edges = pid_profile_observable()->edges(); + std::vector variant_edges; + for (auto &edge : core_edges) { + variant_edges.push_back(Variant{edge}); + } + return variant_edges; + } else { + return Base::call_method(method, parameters); + } + } + std::shared_ptr<::Observables::PidProfileObservable> pid_profile_observable() const { return m_observable; diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py index 3915c07f3dc..fa27c54a0ef 100644 --- a/testsuite/python/observable_cylindrical.py +++ b/testsuite/python/observable_cylindrical.py @@ -104,7 +104,7 @@ def set_particles(self): def calculate_numpy_histogram(self): pol_positions = self.pol_coords() - np_hist, _ = np.histogramdd( + np_hist, np_edges = np.histogramdd( pol_positions, bins=(self.params['n_r_bins'], self.params['n_phi_bins'], @@ -112,7 +112,7 @@ def calculate_numpy_histogram(self): range=[(self.params['min_r'], self.params['max_r']), (self.params['min_phi'], self.params['max_phi']), (self.params['min_z'], self.params['max_z'])]) - return np_hist + return np_hist, np_edges def normalize_with_bin_volume(self, histogram): bin_volume = tests_common.get_cylindrical_bin_volume( @@ -141,9 +141,12 @@ def density_profile_test(self): local_params['axis'] = [0.0, 0.0, 1.0] obs = espressomd.observables.CylindricalDensityProfile(**local_params) core_hist = obs.calculate() - np_hist = self.calculate_numpy_histogram() + core_edges = obs.edges() + np_hist, np_edges = self.calculate_numpy_histogram() np_hist = self.normalize_with_bin_volume(np_hist) np.testing.assert_array_almost_equal(np_hist, core_hist) + for i in range(3): + np.testing.assert_array_almost_equal(np_edges[i], core_edges[i]) self.assertEqual(np.prod(obs.shape()), len(np_hist.flatten())) def velocity_profile_test(self): @@ -161,7 +164,7 @@ def velocity_profile_test(self): core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] - np_hist = self.calculate_numpy_histogram() + np_hist, _ = self.calculate_numpy_histogram() for x in np.nditer(np_hist, op_flags=['readwrite']): if x[...] > 0.0: x[...] /= x[...] @@ -187,7 +190,7 @@ def flux_density_profile_test(self): core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] - np_hist = self.calculate_numpy_histogram() + np_hist, _ = self.calculate_numpy_histogram() np_hist = self.normalize_with_bin_volume(np_hist) np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) np.testing.assert_array_almost_equal( diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index e955808f548..30aa1650d57 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -147,8 +147,9 @@ def LB_fluxdensity_profile_test(self): core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] + core_edges = p.edges() self.pol_positions = self.pol_coords() - np_hist, _ = np.histogramdd( + np_hist, np_edges = np.histogramdd( self.pol_positions, bins=(self.params['n_r_bins'], self.params['n_phi_bins'], @@ -164,6 +165,8 @@ def LB_fluxdensity_profile_test(self): np.testing.assert_array_almost_equal( np_hist * self.v_phi, core_hist_v_phi) np.testing.assert_array_almost_equal(np_hist * self.v_z, core_hist_v_z) + for i in range(3): + np.testing.assert_array_almost_equal(np_edges[i], core_edges[i]) self.assertEqual(np.prod(p.shape()), len(np_hist.flatten()) * 3) def LB_velocity_profile_at_particle_positions_test(self): diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index d28548f2a03..018fd6c4429 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -41,10 +41,27 @@ class ProfileObservablesTest(ut.TestCase): 'min_z': 0.0, 'max_z': 10.0} + def calculate_numpy_histogram(self): + np_hist, np_edges = np.histogramdd( + np.copy(self.system.part[:].pos), + bins=(self.kwargs['n_x_bins'], + self.kwargs['n_y_bins'], + self.kwargs['n_z_bins']), + range=[(self.kwargs['min_x'], self.kwargs['max_x']), + (self.kwargs['min_y'], self.kwargs['max_y']), + (self.kwargs['min_z'], self.kwargs['max_z'])], + density=True) + return np_hist, np_edges + def test_density_profile(self): density_profile = espressomd.observables.DensityProfile(**self.kwargs) obs_data = density_profile.calculate() - self.assertEqual(obs_data[0, 0, 1], 2.0 / self.bin_volume) + obs_edges = np.array(density_profile.edges()) + np_hist, np_edges = self.calculate_numpy_histogram() + np_hist *= len(self.system.part) + np.testing.assert_array_almost_equal(obs_data, np_hist) + for i in range(3): + np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i]) self.assertEqual(np.prod(obs_data.shape), self.kwargs['n_x_bins'] * self.kwargs['n_y_bins'] * self.kwargs['n_z_bins']) diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py index a058ce84f73..ccbb4836d1b 100644 --- a/testsuite/python/observable_profileLB.py +++ b/testsuite/python/observable_profileLB.py @@ -66,6 +66,21 @@ class ObservableProfileLBCommon: system.time_step = TIME_STEP system.cell_system.skin = 0.4 * AGRID + def calculate_numpy_histogram(self): + np_hist, np_edges = np.histogramdd( + np.array([3 * [0]]), + bins=(LB_VELOCITY_PROFILE_PARAMS['n_x_bins'], + LB_VELOCITY_PROFILE_PARAMS['n_y_bins'], + LB_VELOCITY_PROFILE_PARAMS['n_z_bins']), + range=[(LB_VELOCITY_PROFILE_PARAMS['min_x'], + LB_VELOCITY_PROFILE_PARAMS['max_x']), + (LB_VELOCITY_PROFILE_PARAMS['min_y'], + LB_VELOCITY_PROFILE_PARAMS['max_y']), + (LB_VELOCITY_PROFILE_PARAMS['min_z'], + LB_VELOCITY_PROFILE_PARAMS['max_z'])], + density=True) + return np_hist, np_edges + def set_fluid_velocities(self): """Set an x dependent fluid velocity.""" for x in range(int(self.system.box_l[0] / AGRID)): @@ -78,6 +93,10 @@ def test_velocity_profile(self): obs = espressomd.observables.LBVelocityProfile( **LB_VELOCITY_PROFILE_PARAMS) obs_data = obs.calculate() + obs_edges = obs.edges() + _, np_edges = self.calculate_numpy_histogram() + for i in range(3): + np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i]) for x in range(obs_data.shape[0]): for y in range(obs_data.shape[1]): for z in range(obs_data.shape[2]): From 1109b1dd09f0ef4d6f2e780bbc207dab76fbe818 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 30 Mar 2020 20:02:49 +0200 Subject: [PATCH 3/8] Simplify code and add wrapper for histogramdd --- testsuite/python/observable_cylindrical.py | 20 ++-------- testsuite/python/observable_cylindricalLB.py | 13 +------ testsuite/python/observable_profile.py | 26 +++---------- testsuite/python/observable_profileLB.py | 20 ++-------- testsuite/python/tests_common.py | 40 ++++++++++++++++++++ 5 files changed, 55 insertions(+), 64 deletions(-) diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py index fa27c54a0ef..c499ffd4a9f 100644 --- a/testsuite/python/observable_cylindrical.py +++ b/testsuite/python/observable_cylindrical.py @@ -14,7 +14,6 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import sys import numpy as np import unittest as ut import espressomd @@ -28,9 +27,8 @@ class TestCylindricalObservable(ut.TestCase): Testcase for the cylindrical observables. """ - system = espressomd.System(box_l=[1.0, 1.0, 1.0]) + system = espressomd.System(box_l=[15.0, 15.0, 15.0]) system.time_step = 0.01 - system.box_l = [15.0, 15.0, 15.0] system.cell_system.skin = 0.4 params = { @@ -104,14 +102,8 @@ def set_particles(self): def calculate_numpy_histogram(self): pol_positions = self.pol_coords() - np_hist, np_edges = np.histogramdd( - pol_positions, - bins=(self.params['n_r_bins'], - self.params['n_phi_bins'], - self.params['n_z_bins']), - range=[(self.params['min_r'], self.params['max_r']), - (self.params['min_phi'], self.params['max_phi']), - (self.params['min_z'], self.params['max_z'])]) + np_hist, np_edges = tests_common.get_histogram( + pol_positions, self.params, 'cylindrical') return np_hist, np_edges def normalize_with_bin_volume(self, histogram): @@ -218,8 +210,4 @@ def test_hist_z(self): if __name__ == "__main__": - suite = ut.TestSuite() - suite.addTests(ut.TestLoader().loadTestsFromTestCase( - TestCylindricalObservable)) - result = ut.TextTestRunner(verbosity=4).run(suite) - sys.exit(not result.wasSuccessful()) + ut.main() diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index 30aa1650d57..4ea6cd57d41 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -149,17 +149,8 @@ def LB_fluxdensity_profile_test(self): core_hist_v_z = core_hist[:, :, :, 2] core_edges = p.edges() self.pol_positions = self.pol_coords() - np_hist, np_edges = np.histogramdd( - self.pol_positions, - bins=(self.params['n_r_bins'], - self.params['n_phi_bins'], - self.params['n_z_bins']), - range=[(self.params['min_r'], - self.params['max_r']), - (self.params['min_phi'], - self.params['max_phi']), - (self.params['min_z'], - self.params['max_z'])]) + np_hist, np_edges = tests_common.get_histogram( + self.pol_positions, self.params, 'cylindrical') np_hist = self.normalize_with_bin_volume(np_hist) np.testing.assert_array_almost_equal(np_hist * self.v_r, core_hist_v_r) np.testing.assert_array_almost_equal( diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index 018fd6c4429..71189c72d54 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -14,17 +14,16 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . -import sys import unittest as ut import unittest_decorators as utx import numpy as np import espressomd import espressomd.observables +import tests_common class ProfileObservablesTest(ut.TestCase): - system = espressomd.System(box_l=[1.0, 1.0, 1.0]) - system.box_l = [10.0, 10.0, 10.0] + system = espressomd.System(box_l=[10.0, 10.0, 10.0]) system.cell_system.skin = 0.1 system.time_step = 0.01 system.part.add(id=0, pos=[4.0, 4.0, 6.0], v=[0.0, 0.0, 1.0]) @@ -41,23 +40,12 @@ class ProfileObservablesTest(ut.TestCase): 'min_z': 0.0, 'max_z': 10.0} - def calculate_numpy_histogram(self): - np_hist, np_edges = np.histogramdd( - np.copy(self.system.part[:].pos), - bins=(self.kwargs['n_x_bins'], - self.kwargs['n_y_bins'], - self.kwargs['n_z_bins']), - range=[(self.kwargs['min_x'], self.kwargs['max_x']), - (self.kwargs['min_y'], self.kwargs['max_y']), - (self.kwargs['min_z'], self.kwargs['max_z'])], - density=True) - return np_hist, np_edges - def test_density_profile(self): density_profile = espressomd.observables.DensityProfile(**self.kwargs) obs_data = density_profile.calculate() obs_edges = np.array(density_profile.edges()) - np_hist, np_edges = self.calculate_numpy_histogram() + np_hist, np_edges = tests_common.get_histogram( + np.copy(self.system.part[:].pos), self.kwargs, 'cartesian', density=True) np_hist *= len(self.system.part) np.testing.assert_array_almost_equal(obs_data, np_hist) for i in range(3): @@ -88,8 +76,4 @@ def test_flux_density_profile(self): if __name__ == '__main__': - suite = ut.TestSuite() - suite.addTests(ut.TestLoader().loadTestsFromTestCase( - ProfileObservablesTest)) - result = ut.TextTestRunner(verbosity=4).run(suite) - sys.exit(not result.wasSuccessful()) + ut.main() diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py index ccbb4836d1b..ec062fd1215 100644 --- a/testsuite/python/observable_profileLB.py +++ b/testsuite/python/observable_profileLB.py @@ -18,6 +18,7 @@ import unittest_decorators as utx import numpy as np import copy +import tests_common import espressomd import espressomd.lb @@ -66,21 +67,6 @@ class ObservableProfileLBCommon: system.time_step = TIME_STEP system.cell_system.skin = 0.4 * AGRID - def calculate_numpy_histogram(self): - np_hist, np_edges = np.histogramdd( - np.array([3 * [0]]), - bins=(LB_VELOCITY_PROFILE_PARAMS['n_x_bins'], - LB_VELOCITY_PROFILE_PARAMS['n_y_bins'], - LB_VELOCITY_PROFILE_PARAMS['n_z_bins']), - range=[(LB_VELOCITY_PROFILE_PARAMS['min_x'], - LB_VELOCITY_PROFILE_PARAMS['max_x']), - (LB_VELOCITY_PROFILE_PARAMS['min_y'], - LB_VELOCITY_PROFILE_PARAMS['max_y']), - (LB_VELOCITY_PROFILE_PARAMS['min_z'], - LB_VELOCITY_PROFILE_PARAMS['max_z'])], - density=True) - return np_hist, np_edges - def set_fluid_velocities(self): """Set an x dependent fluid velocity.""" for x in range(int(self.system.box_l[0] / AGRID)): @@ -94,7 +80,9 @@ def test_velocity_profile(self): **LB_VELOCITY_PROFILE_PARAMS) obs_data = obs.calculate() obs_edges = obs.edges() - _, np_edges = self.calculate_numpy_histogram() + _, np_edges = tests_common.get_histogram( + np.zeros([1, 3]), LB_VELOCITY_PROFILE_PARAMS, 'cartesian', + density=True) for i in range(3): np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i]) for x in range(obs_data.shape[0]): diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py index ce015a733f7..e685ec59220 100644 --- a/testsuite/python/tests_common.py +++ b/testsuite/python/tests_common.py @@ -275,6 +275,46 @@ def get_cylindrical_bin_volume( phi_bin_size / (2.0 * np.pi) * z_bin_size return bin_volume + +def get_histogram(pos, obs_params, coord_system, **kwargs): + """ + Helper function for ``np.histogramdd()`` and observables. + + Parameters + ---------- + pos : (N, 3) array_like of :obj:`float` + Particle positions. + obs_params : :obj:`dict` + Parameters of the observable. + coord_system : :obj:`str`, \{'cartesian', 'cylindrical'\} + Coordinate system. + \*\*kwargs : + Optional parameters to ``np.histogramdd()``. + + Returns + ------- + array_like + Bins and bin edges. + + """ + if coord_system == 'cartesian': + bins = (obs_params['n_x_bins'], + obs_params['n_y_bins'], + obs_params['n_z_bins']) + extent = [(obs_params['min_x'], obs_params['max_x']), + (obs_params['min_y'], obs_params['max_y']), + (obs_params['min_z'], obs_params['max_z'])] + elif coord_system == 'cylindrical': + bins = (obs_params['n_r_bins'], + obs_params['n_phi_bins'], + obs_params['n_z_bins']) + extent = [(obs_params['min_r'], obs_params['max_r']), + (obs_params['min_phi'], obs_params['max_phi']), + (obs_params['min_z'], obs_params['max_z'])] + else: + raise ValueError("Unknown coord system '{}'".format(coord_system)) + return np.histogramdd(pos, bins=bins, range=extent, **kwargs) + # # Analytical Expressions for interactions # From 9ec9025484e6756c8a2c8c9fbeae1d67cbf340c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 30 Mar 2020 23:05:42 +0200 Subject: [PATCH 4/8] Fix broken sample --- samples/lb_profile.py | 6 ++---- testsuite/scripts/samples/CMakeLists.txt | 1 + testsuite/scripts/samples/test_lb_profile.py | 2 +- 3 files changed, 4 insertions(+), 5 deletions(-) diff --git a/samples/lb_profile.py b/samples/lb_profile.py index dba07d8a31b..9d599af91fa 100644 --- a/samples/lb_profile.py +++ b/samples/lb_profile.py @@ -44,7 +44,7 @@ system.thermostat.set_lb(LB_fluid=lb_fluid, seed=23) fluid_obs = espressomd.observables.CylindricalLBVelocityProfile( center=[5.0, 5.0, 0.0], - axis='z', + axis=[0, 0, 1], n_r_bins=100, n_phi_bins=1, n_z_bins=1, @@ -54,9 +54,7 @@ max_phi=np.pi, min_z=0.0, max_z=10.0, - sampling_delta_x=0.05, - sampling_delta_y=0.05, - sampling_delta_z=1.0) + sampling_density=0.1) cylinder_shape = espressomd.shapes.Cylinder( center=[5.0, 5.0, 5.0], axis=[0, 0, 1], diff --git a/testsuite/scripts/samples/CMakeLists.txt b/testsuite/scripts/samples/CMakeLists.txt index 915f3c59fa1..b217db1fd09 100644 --- a/testsuite/scripts/samples/CMakeLists.txt +++ b/testsuite/scripts/samples/CMakeLists.txt @@ -32,6 +32,7 @@ if(HDF5_FOUND) endif() sample_test(FILE test_lbf.py SUFFIX cpu) sample_test(FILE test_lbf.py SUFFIX gpu LABELS "gpu") +sample_test(FILE test_lb_profile.py) sample_test(FILE test_lj_liquid_distribution.py) sample_test(FILE test_lj_liquid.py) sample_test(FILE test_lj_liquid_structurefactor.py) diff --git a/testsuite/scripts/samples/test_lb_profile.py b/testsuite/scripts/samples/test_lb_profile.py index 079d9339749..ecf91ade378 100644 --- a/testsuite/scripts/samples/test_lb_profile.py +++ b/testsuite/scripts/samples/test_lb_profile.py @@ -31,7 +31,7 @@ def test_fit(self): # compare the simulated data against the analytical solution sim = sample.lb_fluid_profile[:, 0, 0, 2] ana = sample.poiseuille_flow(sample.r, sample.r_max, 0.15) - self.assertLess(np.max(np.abs(sim - ana)), 0.12) + self.assertLess(np.max(np.abs(sim - ana)), 0.16) if __name__ == "__main__": From d96da5b660cbb87d84bdd942574d7255053ca73a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 31 Mar 2020 00:07:16 +0200 Subject: [PATCH 5/8] Apply code review feedback Co-authored-by: Florian Weik --- .../CylindricalProfileObservable.hpp | 33 +++++++----------- src/core/observables/ProfileObservable.hpp | 34 +++++++------------ .../CylindricalLBProfileObservable.hpp | 11 +++--- .../CylindricalPidProfileObservable.hpp | 11 +++--- .../observables/LBProfileObservable.hpp | 11 +++--- .../observables/PidProfileObservable.hpp | 12 +++---- .../observables/ProfileObservable.hpp | 14 ++++++++ 7 files changed, 60 insertions(+), 66 deletions(-) diff --git a/src/core/observables/CylindricalProfileObservable.hpp b/src/core/observables/CylindricalProfileObservable.hpp index 5f4459433e3..e443d9796ec 100644 --- a/src/core/observables/CylindricalProfileObservable.hpp +++ b/src/core/observables/CylindricalProfileObservable.hpp @@ -19,12 +19,13 @@ #ifndef OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP #define OBSERVABLES_CYLINDRICALPROFILEOBSERVABLE_HPP -#include +#include #include #include "Observable.hpp" #include +#include namespace Observables { @@ -55,26 +56,16 @@ class CylindricalProfileObservable : virtual public Observable { } /** Calculate the bin edges for each dimension */ - std::vector> edges() { - std::vector> profile_edges(3); - profile_edges[0] = std::vector(n_r_bins + 1); - profile_edges[1] = std::vector(n_phi_bins + 1); - profile_edges[2] = std::vector(n_z_bins + 1); - std::generate(profile_edges[0].begin(), profile_edges[0].end(), - [n = 0, start = min_r, - bin_size = (max_r - min_r) / n_r_bins]() mutable { - return start + bin_size * n++; - }); - std::generate(profile_edges[1].begin(), profile_edges[1].end(), - [n = 0, start = min_phi, - bin_size = (max_phi - min_phi) / n_phi_bins]() mutable { - return start + bin_size * n++; - }); - std::generate(profile_edges[2].begin(), profile_edges[2].end(), - [n = 0, start = min_z, - bin_size = (max_z - min_z) / n_z_bins]() mutable { - return start + bin_size * n++; - }); + std::array, 3> edges() { + std::array, 3> profile_edges = { + {std::vector(n_r_bins + 1), std::vector(n_phi_bins + 1), + std::vector(n_z_bins + 1)}}; + boost::copy(Utils::make_lin_space(min_r, max_r, n_r_bins + 1), + profile_edges[0].begin()); + boost::copy(Utils::make_lin_space(min_phi, max_phi, n_phi_bins + 1), + profile_edges[1].begin()); + boost::copy(Utils::make_lin_space(min_z, max_z, n_z_bins + 1), + profile_edges[2].begin()); return profile_edges; } }; diff --git a/src/core/observables/ProfileObservable.hpp b/src/core/observables/ProfileObservable.hpp index 9e0b0d2974a..81d1b459c9d 100644 --- a/src/core/observables/ProfileObservable.hpp +++ b/src/core/observables/ProfileObservable.hpp @@ -19,9 +19,11 @@ #ifndef OBSERVABLES_PROFILEOBSERVABLE_HPP #define OBSERVABLES_PROFILEOBSERVABLE_HPP -#include +#include #include +#include + #include "Observable.hpp" namespace Observables { @@ -48,26 +50,16 @@ class ProfileObservable : virtual public Observable { } /** Calculate the bin edges for each dimension */ - std::vector> edges() { - std::vector> profile_edges(3); - profile_edges[0] = std::vector(n_x_bins + 1); - profile_edges[1] = std::vector(n_y_bins + 1); - profile_edges[2] = std::vector(n_z_bins + 1); - std::generate(profile_edges[0].begin(), profile_edges[0].end(), - [n = 0, start = min_x, - bin_size = (max_x - min_x) / n_x_bins]() mutable { - return start + bin_size * n++; - }); - std::generate(profile_edges[1].begin(), profile_edges[1].end(), - [n = 0, start = min_y, - bin_size = (max_y - min_y) / n_y_bins]() mutable { - return start + bin_size * n++; - }); - std::generate(profile_edges[2].begin(), profile_edges[2].end(), - [n = 0, start = min_z, - bin_size = (max_z - min_z) / n_z_bins]() mutable { - return start + bin_size * n++; - }); + std::array, 3> edges() { + std::array, 3> profile_edges = { + {std::vector(n_x_bins + 1), std::vector(n_y_bins + 1), + std::vector(n_z_bins + 1)}}; + boost::copy(Utils::make_lin_space(min_x, max_x, n_x_bins + 1), + profile_edges[0].begin()); + boost::copy(Utils::make_lin_space(min_y, max_y, n_y_bins + 1), + profile_edges[1].begin()); + boost::copy(Utils::make_lin_space(min_z, max_z, n_z_bins + 1), + profile_edges[2].begin()); return profile_edges; } }; diff --git a/src/script_interface/observables/CylindricalLBProfileObservable.hpp b/src/script_interface/observables/CylindricalLBProfileObservable.hpp index 88bbfa3e7d4..6bceb55e14e 100644 --- a/src/script_interface/observables/CylindricalLBProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalLBProfileObservable.hpp @@ -24,6 +24,8 @@ #include "script_interface/auto_parameters/AutoParameters.hpp" +#include +#include #include #include "Observable.hpp" @@ -138,15 +140,12 @@ class CylindricalLBProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "edges") { - auto const core_edges = cylindrical_profile_observable()->edges(); std::vector variant_edges; - for (auto &edge : core_edges) { - variant_edges.push_back(Variant{edge}); - } + boost::copy(cylindrical_profile_observable()->edges(), + std::back_inserter(variant_edges)); return variant_edges; - } else { - return Base::call_method(method, parameters); } + return Base::call_method(method, parameters); } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/CylindricalPidProfileObservable.hpp b/src/script_interface/observables/CylindricalPidProfileObservable.hpp index 93510fe98a2..4ffd84aa5cb 100644 --- a/src/script_interface/observables/CylindricalPidProfileObservable.hpp +++ b/src/script_interface/observables/CylindricalPidProfileObservable.hpp @@ -22,6 +22,8 @@ #ifndef SCRIPT_INTERFACE_OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP #define SCRIPT_INTERFACE_OBSERVABLES_CYLINDRICALPIDPROFILEOBSERVABLE_HPP +#include +#include #include #include "script_interface/ScriptInterface.hpp" @@ -138,15 +140,12 @@ class CylindricalPidProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "edges") { - auto const core_edges = cylindrical_pid_profile_observable()->edges(); std::vector variant_edges; - for (auto &edge : core_edges) { - variant_edges.push_back(Variant{edge}); - } + boost::copy(cylindrical_pid_profile_observable()->edges(), + std::back_inserter(variant_edges)); return variant_edges; - } else { - return Base::call_method(method, parameters); } + return Base::call_method(method, parameters); } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/LBProfileObservable.hpp b/src/script_interface/observables/LBProfileObservable.hpp index aa8105323ab..c616f71a204 100644 --- a/src/script_interface/observables/LBProfileObservable.hpp +++ b/src/script_interface/observables/LBProfileObservable.hpp @@ -24,6 +24,8 @@ #include "script_interface/auto_parameters/AutoParameters.hpp" +#include +#include #include #include "Observable.hpp" @@ -149,15 +151,12 @@ class LBProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "edges") { - auto const core_edges = profile_observable()->edges(); std::vector variant_edges; - for (auto &edge : core_edges) { - variant_edges.push_back(Variant{edge}); - } + boost::copy(profile_observable()->edges(), + std::back_inserter(variant_edges)); return variant_edges; - } else { - return Base::call_method(method, parameters); } + return Base::call_method(method, parameters); } std::shared_ptr<::Observables::Observable> observable() const override { diff --git a/src/script_interface/observables/PidProfileObservable.hpp b/src/script_interface/observables/PidProfileObservable.hpp index 60a4b10d81a..1a817d8cb76 100644 --- a/src/script_interface/observables/PidProfileObservable.hpp +++ b/src/script_interface/observables/PidProfileObservable.hpp @@ -23,6 +23,9 @@ #define SCRIPT_INTERFACE_OBSERVABLES_PIDPROFILEOBSERVABLE_HPP #include "script_interface/auto_parameters/AutoParameters.hpp" + +#include +#include #include #include "core/observables/DensityProfile.hpp" @@ -114,15 +117,12 @@ class PidProfileObservable Variant call_method(std::string const &method, VariantMap const ¶meters) override { if (method == "edges") { - auto const core_edges = pid_profile_observable()->edges(); std::vector variant_edges; - for (auto &edge : core_edges) { - variant_edges.push_back(Variant{edge}); - } + boost::copy(pid_profile_observable()->edges(), + std::back_inserter(variant_edges)); return variant_edges; - } else { - return Base::call_method(method, parameters); } + return Base::call_method(method, parameters); } std::shared_ptr<::Observables::PidProfileObservable> diff --git a/src/script_interface/observables/ProfileObservable.hpp b/src/script_interface/observables/ProfileObservable.hpp index dc464385003..6916b6f5a43 100644 --- a/src/script_interface/observables/ProfileObservable.hpp +++ b/src/script_interface/observables/ProfileObservable.hpp @@ -34,7 +34,10 @@ namespace Observables { template class ProfileObservable : public AutoParameters, Observable> { + using Base = AutoParameters, Observable>; + public: + using Base::Base; ProfileObservable() { this->add_parameters( {{"n_x_bins", @@ -95,6 +98,17 @@ class ProfileObservable void construct(VariantMap const ¶ms) override {} + Variant call_method(std::string const &method, + VariantMap const ¶meters) override { + if (method == "edges") { + std::vector variant_edges; + boost::copy(profile_observable()->edges(), + std::back_inserter(variant_edges)); + return variant_edges; + } + return Base::call_method(method, parameters); + } + std::shared_ptr<::Observables::ProfileObservable> profile_observable() const { return m_observable; } From a304ae2a7b6b2df79867a58f69677a86daba0722 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 31 Mar 2020 00:15:37 +0200 Subject: [PATCH 6/8] Fix broken python tests Optional parameter `density=True` is not allowed in np.histogramdd() in numpy 1.11. --- testsuite/python/observable_profile.py | 3 ++- testsuite/python/observable_profileLB.py | 2 +- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index 71189c72d54..6df405d4144 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -45,7 +45,8 @@ def test_density_profile(self): obs_data = density_profile.calculate() obs_edges = np.array(density_profile.edges()) np_hist, np_edges = tests_common.get_histogram( - np.copy(self.system.part[:].pos), self.kwargs, 'cartesian', density=True) + np.copy(self.system.part[:].pos), self.kwargs, 'cartesian', + normed=True) np_hist *= len(self.system.part) np.testing.assert_array_almost_equal(obs_data, np_hist) for i in range(3): diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py index ec062fd1215..86e72ce051d 100644 --- a/testsuite/python/observable_profileLB.py +++ b/testsuite/python/observable_profileLB.py @@ -82,7 +82,7 @@ def test_velocity_profile(self): obs_edges = obs.edges() _, np_edges = tests_common.get_histogram( np.zeros([1, 3]), LB_VELOCITY_PROFILE_PARAMS, 'cartesian', - density=True) + normed=True) for i in range(3): np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i]) for x in range(obs_data.shape[0]): From 8ccfe3495a0d15f0edf752d0614dc8fbd9e34edf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 15 Apr 2020 17:56:25 +0200 Subject: [PATCH 7/8] ProfileObservable: add bin_edges(), bin_centers() Reproduce the core ProfileObservable class in the Python interface. Calculate the bin edges and centers in Python. --- src/python/espressomd/observables.py | 63 ++++++++++++++++---- testsuite/python/observable_cylindrical.py | 2 +- testsuite/python/observable_cylindricalLB.py | 2 +- testsuite/python/observable_profile.py | 12 +++- testsuite/python/observable_profileLB.py | 2 +- 5 files changed, 66 insertions(+), 15 deletions(-) diff --git a/src/python/espressomd/observables.py b/src/python/espressomd/observables.py index f3e8534760e..8ff0cc98384 100644 --- a/src/python/espressomd/observables.py +++ b/src/python/espressomd/observables.py @@ -14,20 +14,61 @@ # # You should have received a copy of the GNU General Public License # along with this program. If not, see . +import itertools import numpy as np from .script_interface import ScriptInterfaceHelper, script_interface_register @script_interface_register class Observable(ScriptInterfaceHelper): + """ + Base class for all observables. + + Methods + ------- + shape() + Return the shape of the observable. + """ _so_name = "Observables::Observable" - _so_bind_methods = ("shape", "edges") + _so_bind_methods = ("shape",) _so_creation_policy = "LOCAL" def calculate(self): return np.array(self.call_method("calculate")).reshape(self.shape()) +class ProfileObservable(Observable): + """ + Base class for histogram-based observables. + """ + + def bin_edges(self): + """ + Returns + ------- + :obj:`ndarray` of :obj:`float` + Positions between the bins. If the histogram has dimensions + ``(M,N,O)``, the bin edges have dimensions ``(M+1,N+1,O+1,3)``. + """ + edges = self.call_method("edges") + shape = list(map(len, edges)) + [len(edges)] + return np.array(list(itertools.product(*edges))).reshape(shape) + + def bin_centers(self): + """ + Returns + ------- + :obj:`ndarray` of :obj:`float` + Positions of the bins centers. If the histogram has dimensions + ``(M,N,O)``, the bin centers have dimensions ``(M,N,O,3)``. + """ + edges = self.call_method("edges") + for i, edge in enumerate(edges): + edges[i] = np.array(edge[:-1]) + (edge[1] - edge[0]) / 2 + shape = list(map(len, edges)) + [len(edges)] + return np.array(list(itertools.product(*edges))).reshape(shape) + + @script_interface_register class ComPosition(Observable): @@ -93,7 +134,7 @@ class Current(Observable): @script_interface_register -class DensityProfile(Observable): +class DensityProfile(ProfileObservable): """Calculates the particle density profile for particles with given ids. @@ -149,7 +190,7 @@ class DipoleMoment(Observable): @script_interface_register -class FluxDensityProfile(Observable): +class FluxDensityProfile(ProfileObservable): """Calculates the particle flux density for particles with given ids. @@ -187,7 +228,7 @@ class FluxDensityProfile(Observable): @script_interface_register -class ForceDensityProfile(Observable): +class ForceDensityProfile(ProfileObservable): """Calculates the force density profile for particles with given ids. @@ -225,7 +266,7 @@ class ForceDensityProfile(Observable): @script_interface_register -class LBVelocityProfile(Observable): +class LBVelocityProfile(ProfileObservable): """Calculates the LB fluid velocity profile. @@ -589,7 +630,7 @@ class DPDStress(Observable): @script_interface_register -class CylindricalDensityProfile(Observable): +class CylindricalDensityProfile(ProfileObservable): """Calculates the particle density in cylindrical coordinates. @@ -629,7 +670,7 @@ class CylindricalDensityProfile(Observable): @script_interface_register -class CylindricalFluxDensityProfile(Observable): +class CylindricalFluxDensityProfile(ProfileObservable): """Calculates the particle flux density in cylindrical coordinates. @@ -671,7 +712,7 @@ class CylindricalFluxDensityProfile(Observable): @script_interface_register -class CylindricalLBFluxDensityProfileAtParticlePositions(Observable): +class CylindricalLBFluxDensityProfileAtParticlePositions(ProfileObservable): """Calculates the LB fluid flux density at the particle positions in cylindrical coordinates. @@ -714,7 +755,7 @@ class CylindricalLBFluxDensityProfileAtParticlePositions(Observable): @script_interface_register -class CylindricalLBVelocityProfileAtParticlePositions(Observable): +class CylindricalLBVelocityProfileAtParticlePositions(ProfileObservable): """Calculates the LB fluid velocity at the particle positions in cylindrical coordinates. @@ -757,7 +798,7 @@ class CylindricalLBVelocityProfileAtParticlePositions(Observable): @script_interface_register -class CylindricalVelocityProfile(Observable): +class CylindricalVelocityProfile(ProfileObservable): """Calculates the particle velocity profile in cylindrical coordinates. @@ -799,7 +840,7 @@ class CylindricalVelocityProfile(Observable): @script_interface_register -class CylindricalLBVelocityProfile(Observable): +class CylindricalLBVelocityProfile(ProfileObservable): """Calculates the LB fluid velocity profile in cylindrical coordinates. diff --git a/testsuite/python/observable_cylindrical.py b/testsuite/python/observable_cylindrical.py index c499ffd4a9f..3be0875408e 100644 --- a/testsuite/python/observable_cylindrical.py +++ b/testsuite/python/observable_cylindrical.py @@ -133,7 +133,7 @@ def density_profile_test(self): local_params['axis'] = [0.0, 0.0, 1.0] obs = espressomd.observables.CylindricalDensityProfile(**local_params) core_hist = obs.calculate() - core_edges = obs.edges() + core_edges = obs.call_method("edges") np_hist, np_edges = self.calculate_numpy_histogram() np_hist = self.normalize_with_bin_volume(np_hist) np.testing.assert_array_almost_equal(np_hist, core_hist) diff --git a/testsuite/python/observable_cylindricalLB.py b/testsuite/python/observable_cylindricalLB.py index 4ea6cd57d41..55b9e8caed7 100644 --- a/testsuite/python/observable_cylindricalLB.py +++ b/testsuite/python/observable_cylindricalLB.py @@ -147,7 +147,7 @@ def LB_fluxdensity_profile_test(self): core_hist_v_r = core_hist[:, :, :, 0] core_hist_v_phi = core_hist[:, :, :, 1] core_hist_v_z = core_hist[:, :, :, 2] - core_edges = p.edges() + core_edges = p.call_method("edges") self.pol_positions = self.pol_coords() np_hist, np_edges = tests_common.get_histogram( self.pol_positions, self.params, 'cylindrical') diff --git a/testsuite/python/observable_profile.py b/testsuite/python/observable_profile.py index 6df405d4144..4661bf65c54 100644 --- a/testsuite/python/observable_profile.py +++ b/testsuite/python/observable_profile.py @@ -43,7 +43,9 @@ class ProfileObservablesTest(ut.TestCase): def test_density_profile(self): density_profile = espressomd.observables.DensityProfile(**self.kwargs) obs_data = density_profile.calculate() - obs_edges = np.array(density_profile.edges()) + obs_edges = density_profile.call_method("edges") + obs_bin_edges = density_profile.bin_edges() + obs_bin_centers = density_profile.bin_centers() np_hist, np_edges = tests_common.get_histogram( np.copy(self.system.part[:].pos), self.kwargs, 'cartesian', normed=True) @@ -53,6 +55,14 @@ def test_density_profile(self): np.testing.assert_array_almost_equal(obs_edges[i], np_edges[i]) self.assertEqual(np.prod(obs_data.shape), self.kwargs['n_x_bins'] * self.kwargs['n_y_bins'] * self.kwargs['n_z_bins']) + np.testing.assert_array_almost_equal( + obs_bin_edges[0, 0, 0], + [self.kwargs['min_x'], self.kwargs['min_y'], self.kwargs['min_z']]) + np.testing.assert_array_almost_equal( + obs_bin_edges[-1, -1, -1], + [self.kwargs['max_x'], self.kwargs['max_y'], self.kwargs['max_z']]) + np.testing.assert_array_almost_equal(obs_bin_centers[0, 0, 0], + [2.5, 2.5, 2.5]) @utx.skipIfMissingFeatures("EXTERNAL_FORCES") def test_force_density_profile(self): diff --git a/testsuite/python/observable_profileLB.py b/testsuite/python/observable_profileLB.py index 1cbd9905fca..fcfffb374eb 100644 --- a/testsuite/python/observable_profileLB.py +++ b/testsuite/python/observable_profileLB.py @@ -79,7 +79,7 @@ def test_velocity_profile(self): obs = espressomd.observables.LBVelocityProfile( **LB_VELOCITY_PROFILE_PARAMS) obs_data = obs.calculate() - obs_edges = obs.edges() + obs_edges = obs.call_method("edges") _, np_edges = tests_common.get_histogram( np.zeros([1, 3]), LB_VELOCITY_PROFILE_PARAMS, 'cartesian', normed=True) From 2463118fac4ea8701ae9bfb347793373800a6562 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 15 Apr 2020 18:55:05 +0200 Subject: [PATCH 8/8] Update Observables documentation Discuss ongoing efforts to convert functionality from the Analysis module to the Observable framework. Rewrite section explaining how to instantiate and use observables. Document the new bin_edges() and bin_centers() methods and provide a matplotlib script to illustrate how they differ. --- doc/sphinx/analysis.rst | 76 +++++++++++++++++++++++++++++++++++------ 1 file changed, 65 insertions(+), 11 deletions(-) diff --git a/doc/sphinx/analysis.rst b/doc/sphinx/analysis.rst index d1862881588..47a7a9cb76f 100644 --- a/doc/sphinx/analysis.rst +++ b/doc/sphinx/analysis.rst @@ -345,8 +345,10 @@ Both versions are equivalent in the :math:`N\rightarrow \infty` limit but give n Observables and correlators --------------------------- -Analysis in the core is a new concept introduced in since version 3.1. -It was motivated by the fact, that sometimes it is desirable that the +Observables extract properties of the particles and the LB fluid and +return either the raw data or a statistic derived from them. The +Observable framework is progressively replacing the Analysis framework. +This is motivated by the fact, that sometimes it is desirable that the analysis functions do more than just return a value to the scripting interface. For some observables it is desirable to be sampled every few integration steps. In addition, it should be possible to pass the @@ -380,12 +382,13 @@ not all observables carry out their calculations in parallel. Instead, the entire particle configuration is collected on the head node, and the calculations are carried out there. This is only performance-relevant if the number of processor cores is large and/or interactions are calculated very frequently. -.. _Creating an observable: +.. _Using observables: -Creating an observable -~~~~~~~~~~~~~~~~~~~~~~ +Using observables +~~~~~~~~~~~~~~~~~ -The observables are represented as Python classes derived from :class:`espressomd.observables.Observable`. They are contained in +The observables are represented as Python classes derived from +:class:`espressomd.observables.Observable`. They are contained in the ``espressomd.observables`` module. An observable is instantiated as follows @@ -397,17 +400,70 @@ follows Here, the keyword argument ``ids`` specifies the ids of the particles, which the observable should take into account. -The current value of an observable can be obtained using its calculate()-method:: +The current value of an observable can be obtained using its +:meth:`~espressomd.observables.Observable.calculate` method:: print(part_pos.calculate()) +Profile observables have additional methods +:meth:`~espressomd.observables.ProfileObservable.bin_centers` and +:meth:`~espressomd.observables.ProfileObservable.bin_edges` to facilitate +plotting of histogram slices with functions that require either bin centers +or bin edges for the axes. Example:: + + import matplotlib.pyplot as plt + import numpy as np + import espressomd + import espressomd.observables + + system = espressomd.System(box_l=[10.0, 10.0, 10.0]) + system.part.add(id=0, pos=[4.0, 3.0, 6.0]) + system.part.add(id=1, pos=[7.0, 3.0, 6.0]) + + # histogram in Cartesian coordinates + density_profile = espressomd.observables.DensityProfile( + ids=[0, 1], + n_x_bins=8, min_x=1.0, max_x=9.0, + n_y_bins=8, min_y=1.0, max_y=9.0, + n_z_bins=4, min_z=4.0, max_z=8.0) + obs_data = density_profile.calculate() + obs_bins = density_profile.bin_centers() + + # 1D slice: requires bin centers + plt.plot(obs_bins[:, 2, 2, 0], obs_data[:, 2, 2]) + plt.show() + + # 2D slice: requires extent + plt.imshow(obs_data[:, :, 2].T, origin='lower', + extent=[density_profile.min_x, density_profile.max_x, + density_profile.min_y, density_profile.max_y]) + plt.show() + + # histogram in cylindrical coordinates + density_profile = espressomd.observables.CylindricalDensityProfile( + ids=[0, 1], center=[5.0, 5.0, 0.0], axis=[0, 0, 1], + n_r_bins=8, min_r=0.0, max_r=4.0, + n_phi_bins=16, min_phi=-np.pi, max_phi=np.pi, + n_z_bins=4, min_z=4.0, max_z=8.0) + obs_data = density_profile.calculate() + obs_bins = density_profile.bin_edges() + + # 2D slice: requires bin edges + fig = plt.figure() + ax = fig.add_subplot(111, polar='True') + r = obs_bins[:, 0, 0, 0] + phi = obs_bins[0, :, 0, 1] + ax.pcolormesh(phi, r, obs_data[:, :, 2]) + plt.show() + + .. _Available observables: Available observables ~~~~~~~~~~~~~~~~~~~~~ -The following list contains some of the available observables. You can find documentation for -all available observables in :mod:`espressomd.observables`. +The following list contains some of the available observables. You can find +documentation for all available observables in :mod:`espressomd.observables`. - Observables working on a given set of particles: @@ -468,7 +524,6 @@ all available observables in :mod:`espressomd.observables`. - :class:`~espressomd.observables.CylindricalLBVelocityProfileAtParticlePositions` - - System-wide observables - :class:`~espressomd.observables.StressTensor`: Total stress tensor (see :ref:`stress tensor`) @@ -476,7 +531,6 @@ all available observables in :mod:`espressomd.observables`. - :class:`~espressomd.observables.DPDStress` - .. _Correlations: Correlations