Observables sanity checks (#4255)

Fixes #4212

Description of changes:
- catch memory errors when passing an invalid observable to the correlator framework
- raise python exceptions when instantiating observables with invalid parameters

src/core/observables/BondAngles.hpp

@@ -25,8 +25,11 @@
 
 #include <utils/Vector.hpp>
 
+#include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <stdexcept>
+#include <utility>
 #include <vector>
 
 namespace Observables {
@@ -38,6 +41,10 @@ namespace Observables {
 class BondAngles : public PidObservable {
 public:
   using PidObservable::PidObservable;
+  explicit BondAngles(std::vector<int> ids) : PidObservable(std::move(ids)) {
+    if (this->ids().size() < 3)
+      throw std::runtime_error("At least 3 particles are required");
+  }
 
   std::vector<double>
   evaluate(Utils::Span<std::reference_wrapper<const Particle>> particles,
@@ -67,7 +74,10 @@ class BondAngles : public PidObservable {
     }
     return res;
   }
-  std::vector<size_t> shape() const override { return {ids().size() - 2}; }
+  std::vector<size_t> shape() const override {
+    assert(ids().size() >= 2);
+    return {ids().size() - 2};
+  }
 };
 
 } // Namespace Observables

src/core/observables/BondDihedrals.hpp

@@ -26,8 +26,11 @@
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
 
+#include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <stdexcept>
+#include <utility>
 #include <vector>
 
 namespace Observables {
@@ -44,6 +47,10 @@ namespace Observables {
 class BondDihedrals : public PidObservable {
 public:
   using PidObservable::PidObservable;
+  explicit BondDihedrals(std::vector<int> ids) : PidObservable(std::move(ids)) {
+    if (this->ids().size() < 4)
+      throw std::runtime_error("At least 4 particles are required");
+  }
 
   std::vector<double>
   evaluate(Utils::Span<std::reference_wrapper<const Particle>> particles,
@@ -67,7 +74,10 @@ class BondDihedrals : public PidObservable {
     }
     return res;
   }
-  std::vector<size_t> shape() const override { return {ids().size() - 3}; }
+  std::vector<size_t> shape() const override {
+    assert(ids().size() >= 3);
+    return {ids().size() - 3};
+  }
 };
 
 } // Namespace Observables

src/core/observables/CosPersistenceAngles.hpp

@@ -26,8 +26,11 @@
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
 
+#include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <stdexcept>
+#include <utility>
 #include <vector>
 
 namespace Observables {
@@ -39,6 +42,12 @@ namespace Observables {
 class CosPersistenceAngles : public PidObservable {
 public:
   using PidObservable::PidObservable;
+  explicit CosPersistenceAngles(std::vector<int> ids)
+      : PidObservable(std::move(ids)) {
+    if (this->ids().size() < 3)
+      throw std::runtime_error("At least 3 particles are required");
+  }
+
   std::vector<double>
   evaluate(Utils::Span<std::reference_wrapper<const Particle>> particles,
            const ParticleObservables::traits<Particle> &traits) const override {
@@ -65,7 +74,10 @@ class CosPersistenceAngles : public PidObservable {
 
     return angles;
   }
-  std::vector<size_t> shape() const override { return {ids().size() - 2}; }
+  std::vector<size_t> shape() const override {
+    assert(ids().size() >= 2);
+    return {ids().size() - 2};
+  }
 };
 
 } // Namespace Observables

src/core/observables/LBProfileObservable.hpp

@@ -25,8 +25,10 @@
 #include <utils/Vector.hpp>
 
 #include <array>
+#include <cassert>
 #include <cmath>
 #include <cstddef>
+#include <stdexcept>
 
 namespace Observables {
 
@@ -44,6 +46,18 @@ class LBProfileObservable : public ProfileObservable {
         sampling_offset{sampling_offset_x, sampling_offset_y,
                         sampling_offset_z},
         allow_empty_bins(allow_empty_bins) {
+    if (sampling_delta[0] <= 0.)
+      throw std::domain_error("sampling_delta_x has to be > 0");
+    if (sampling_delta[1] <= 0.)
+      throw std::domain_error("sampling_delta_y has to be > 0");
+    if (sampling_delta[2] <= 0.)
+      throw std::domain_error("sampling_delta_z has to be > 0");
+    if (sampling_offset[0] < 0.)
+      throw std::domain_error("sampling_offset_x has to be >= 0");
+    if (sampling_offset[1] < 0.)
+      throw std::domain_error("sampling_offset_y has to be >= 0");
+    if (sampling_offset[2] < 0.)
+      throw std::domain_error("sampling_offset_z has to be >= 0");
     calculate_sampling_positions();
   }
   std::array<double, 3> sampling_delta;
@@ -53,9 +67,10 @@ class LBProfileObservable : public ProfileObservable {
   void calculate_sampling_positions() {
     auto const lim = limits();
     sampling_positions.clear();
-    if (sampling_delta[0] == 0 or sampling_delta[1] == 0 or
-        sampling_delta[2] == 0)
-      throw std::runtime_error("Parameter delta_x/y/z must not be zero!");
+    assert(sampling_delta[0] > 0. and sampling_delta[1] > 0. and
+           sampling_delta[2] > 0.);
+    assert(sampling_offset[0] >= 0. and sampling_offset[1] >= 0. and
+           sampling_offset[2] >= 0.);
     const auto n_samples_x = static_cast<size_t>(
         std::rint((lim[0].second - lim[0].first) / sampling_delta[0]));
     const auto n_samples_y = static_cast<size_t>(

src/core/observables/ParticleDistances.hpp

@@ -26,7 +26,10 @@
 #include <utils/Span.hpp>
 #include <utils/Vector.hpp>
 
+#include <cassert>
 #include <cstddef>
+#include <stdexcept>
+#include <utility>
 #include <vector>
 
 namespace Observables {
@@ -38,6 +41,11 @@ namespace Observables {
 class ParticleDistances : public PidObservable {
 public:
   using PidObservable::PidObservable;
+  explicit ParticleDistances(std::vector<int> ids)
+      : PidObservable(std::move(ids)) {
+    if (this->ids().size() < 2)
+      throw std::runtime_error("At least 2 particles are required");
+  }
 
   std::vector<double>
   evaluate(Utils::Span<std::reference_wrapper<const Particle>> particles,
@@ -51,7 +59,10 @@ class ParticleDistances : public PidObservable {
     }
     return res;
   }
-  std::vector<size_t> shape() const override { return {ids().size() - 1}; }
+  std::vector<size_t> shape() const override {
+    assert(!ids().empty());
+    return {ids().size() - 1};
+  }
 };
 
 } // Namespace Observables

src/core/observables/ProfileObservable.hpp

@@ -27,6 +27,7 @@
 
 #include <array>
 #include <cstddef>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -54,6 +55,12 @@ class ProfileObservable : virtual public Observable {
       throw std::runtime_error("max_y has to be > min_y");
     if (max_z <= min_z)
       throw std::runtime_error("max_z has to be > min_z");
+    if (n_x_bins <= 0)
+      throw std::domain_error("n_x_bins has to be >= 1");
+    if (n_y_bins <= 0)
+      throw std::domain_error("n_y_bins has to be >= 1");
+    if (n_z_bins <= 0)
+      throw std::domain_error("n_z_bins has to be >= 1");
   }
 
   std::vector<size_t> shape() const override {

src/core/observables/RDF.hpp

@@ -26,6 +26,7 @@
 #include <utils/Span.hpp>
 
 #include <cstddef>
+#include <stdexcept>
 #include <utility>
 #include <vector>
 
@@ -54,7 +55,12 @@ class RDF : public Observable {
   explicit RDF(std::vector<int> ids1, std::vector<int> ids2, int n_r_bins,
                double min_r, double max_r)
       : m_ids1(std::move(ids1)), m_ids2(std::move(ids2)), min_r(min_r),
-        max_r(max_r), n_r_bins(n_r_bins) {}
+        max_r(max_r), n_r_bins(n_r_bins) {
+    if (max_r <= min_r)
+      throw std::runtime_error("max_r has to be > min_r");
+    if (n_r_bins <= 0)
+      throw std::domain_error("n_r_bins has to be >= 1");
+  }
   std::vector<double> operator()() const final;
 
   std::vector<int> &ids1() { return m_ids1; }

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -302,7 +302,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
 
     // check integrated trajectory; the time step is chosen
     // small enough so that particles don't travel too far
+#ifndef NDEBUG
     auto const pos_com = Utils::Vector3d{box_center, box_center, 1.0};
+#endif
     auto const pids = std::vector<int>{pid1, pid2, pid3};
     for (int i = 0; i < 10; ++i) {
       std::unordered_map<int, Utils::Vector3d> expected;

testsuite/python/observable_chain.py

@@ -89,6 +89,11 @@ def test_ParticleDistances(self):
                 res_obs_chain, distances, decimal=9,
                 err_msg="Data did not agree for observable ParticleDistances")
 
+        # check exceptions
+        for i in range(2):
+            with self.assertRaises(RuntimeError):
+                espressomd.observables.ParticleDistances(ids=np.arange(i))
+
     def test_BondAngles(self):
         """
         Check BondAngles, for a particle triple and for a chain.
@@ -127,6 +132,11 @@ def test_BondAngles(self):
                 res_obs_chain, angles, decimal=9,
                 err_msg="Data did not agree for observable BondAngles")
 
+        # check exceptions
+        for i in range(3):
+            with self.assertRaises(RuntimeError):
+                espressomd.observables.BondAngles(ids=np.arange(i))
+
     def test_BondDihedrals(self):
         """
         Check BondDihedrals, for a particle quadruple and for a chain.
@@ -201,6 +211,11 @@ def place_particles(bl, offset):
                         res_obs_chain, [dih1, dih2], decimal=9,
                         err_msg="Data did not agree for observable BondDihedrals")
 
+        # check exceptions
+        for i in range(4):
+            with self.assertRaises(RuntimeError):
+                espressomd.observables.BondDihedrals(ids=np.arange(i))
+
     def test_CosPersistenceAngles(self):
         # First test: compare with python implementation
         self.system.part.clear()
@@ -222,6 +237,11 @@ def test_CosPersistenceAngles(self):
         expected = np.arange(1, 9) * delta_phi
         np.testing.assert_allclose(obs.calculate(), np.cos(expected))
 
+        # check exceptions
+        for i in range(3):
+            with self.assertRaises(RuntimeError):
+                espressomd.observables.CosPersistenceAngles(ids=np.arange(i))
+
 
 if __name__ == "__main__":
     ut.main()

testsuite/python/observable_profile.py

@@ -143,6 +143,16 @@ def test_pid_profile_interface(self):
         obs_bin_edges = observable.bin_edges()
         np.testing.assert_array_equal(obs_bin_edges[-1, -1, -1], [7, 8, 9])
 
+    def test_exceptions(self):
+        params = self.kwargs.copy()
+        for axis in 'xyz':
+            with self.assertRaises(RuntimeError):
+                espressomd.observables.DensityProfile(
+                    **{**params, f'min_{axis}': 100.})
+            with self.assertRaises(ValueError):
+                espressomd.observables.DensityProfile(
+                    **{**params, f'n_{axis}_bins': 0})
+
 
 if __name__ == '__main__':
     ut.main()

testsuite/python/observable_profileLB.py

@@ -95,15 +95,6 @@ def test_velocity_profile(self):
                          LB_VELOCITY_PROFILE_PARAMS['n_y_bins'] *
                          LB_VELOCITY_PROFILE_PARAMS['n_z_bins'] * 3)
 
-    def test_error_sampling_delta_of_0(self):
-        lb_velocity_params_local = copy.copy(LB_VELOCITY_PROFILE_PARAMS)
-        lb_velocity_params_local['sampling_delta_x'] = 0.0
-        lb_velocity_params_local['sampling_delta_y'] = 0.0
-        lb_velocity_params_local['sampling_delta_z'] = 0.0
-        with self.assertRaises(RuntimeError):
-            _ = espressomd.observables.LBVelocityProfile(
-                **lb_velocity_params_local)
-
     def test_error_if_no_LB(self):
         self.system.actors.clear()
         obs = espressomd.observables.LBVelocityProfile(
@@ -119,6 +110,22 @@ def test_error_if_empty_bin(self):
         with self.assertRaises(RuntimeError):
             obs.calculate()
 
+    def test_exceptions(self):
+        params = LB_VELOCITY_PROFILE_PARAMS.copy()
+        for axis in 'xyz':
+            with self.assertRaisesRegex(RuntimeError, f'max_{axis} has to be > min_{axis}'):
+                espressomd.observables.LBVelocityProfile(
+                    **{**params, f'min_{axis}': 100.})
+            with self.assertRaisesRegex(ValueError, f'n_{axis}_bins has to be >= 1'):
+                espressomd.observables.LBVelocityProfile(
+                    **{**params, f'n_{axis}_bins': 0})
+            with self.assertRaisesRegex(ValueError, f'sampling_delta_{axis} has to be > 0'):
+                espressomd.observables.LBVelocityProfile(
+                    **{**params, f'sampling_delta_{axis}': 0})
+            with self.assertRaisesRegex(ValueError, f'sampling_offset_{axis} has to be >= 0'):
+                espressomd.observables.LBVelocityProfile(
+                    **{**params, f'sampling_offset_{axis}': -1e-8})
+
     def test_lb_profile_interface(self):
         # test setters and getters
         params = LB_VELOCITY_PROFILE_PARAMS.copy()

testsuite/python/rdf.py

@@ -139,6 +139,11 @@ def test_rdf_interface(self):
         # check bin centers
         obs_bin_centers = observable.bin_centers()
         np.testing.assert_array_almost_equal(obs_bin_centers, [1.25, 1.75])
+        # check exceptions
+        with self.assertRaises(RuntimeError):
+            espressomd.observables.RDF(**{**params, 'min_r': 100.})
+        with self.assertRaises(ValueError):
+            espressomd.observables.RDF(**{**params, 'n_r_bins': 0})
 
 
 if __name__ == "__main__":