Observables and correlators/accumulators maintainance

src/core/accumulators/Correlator.cpp

@@ -65,7 +65,7 @@ std::vector<double> compress_discard2(std::vector<double> const &A1,
 
 std::vector<double> scalar_product(std::vector<double> const &A,
                                    std::vector<double> const &B,
-                                   Utils::Vector3d) {
+                                   Utils::Vector3d const &) {
   if (A.size() != B.size()) {
     throw std::runtime_error(
         "Error in scalar product: The vector sizes do not match");
@@ -77,7 +77,7 @@ std::vector<double> scalar_product(std::vector<double> const &A,
 
 std::vector<double> componentwise_product(std::vector<double> const &A,
                                           std::vector<double> const &B,
-                                          Utils::Vector3d) {
+                                          Utils::Vector3d const &) {
   std::vector<double> C(A.size());
   if (A.size() != B.size()) {
     throw std::runtime_error(
@@ -91,13 +91,13 @@ std::vector<double> componentwise_product(std::vector<double> const &A,
 
 std::vector<double> tensor_product(std::vector<double> const &A,
                                    std::vector<double> const &B,
-                                   Utils::Vector3d) {
+                                   Utils::Vector3d const &) {
   std::vector<double> C(A.size() * B.size());
   auto C_it = C.begin();
 
-  for (auto A_it = A.begin(); A_it != A.begin(); ++A_it) {
-    for (double B_it : B) {
-      *(C_it++) = *A_it * B_it;
+  for (double a : A) {
+    for (double b : B) {
+      *(C_it++) = a * b;
     }
   }
   assert(C_it == C.end());
@@ -107,7 +107,7 @@ std::vector<double> tensor_product(std::vector<double> const &A,
 
 std::vector<double> square_distance_componentwise(std::vector<double> const &A,
                                                   std::vector<double> const &B,
-                                                  Utils::Vector3d) {
+                                                  Utils::Vector3d const &) {
   if (A.size() != B.size()) {
     throw std::runtime_error(
         "Error in square distance componentwise: The vector sizes do not "
@@ -127,7 +127,7 @@ std::vector<double> square_distance_componentwise(std::vector<double> const &A,
 // sets w
 std::vector<double> fcs_acf(std::vector<double> const &A,
                             std::vector<double> const &B,
-                            Utils::Vector3d wsquare) {
+                            Utils::Vector3d const &wsquare) {
   if (A.size() != B.size()) {
     throw std::runtime_error(
         "Error in fcs acf: The vector sizes do not match.");

src/core/accumulators/Correlator.hpp

@@ -259,9 +259,9 @@ class Correlator : public AccumulatorBase {
   unsigned int dim_A; ///< dimensionality of A
   unsigned int dim_B; ///< dimensionality of B
 
-  using correlation_operation_type =
-      std::vector<double> (*)(std::vector<double> const &,
-                              std::vector<double> const &, Utils::Vector3d);
+  using correlation_operation_type = std::vector<double> (*)(
+      std::vector<double> const &, std::vector<double> const &,
+      Utils::Vector3d const &);
 
   correlation_operation_type corr_operation;
 

testsuite/python/correlation.py

@@ -27,38 +27,144 @@
 
 class CorrelatorTest(ut.TestCase):
     # Handle for espresso system
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
+    system = espressomd.System(box_l=[10, 10, 10])
+    system.cell_system.skin = 0.4
+    system.time_step = 0.01
 
-    def test(self):
+    def tearDown(self):
+        self.system.part.clear()
+
+    def calc_tau(self, time_step, tau_lin, length):
+        tau = []
+        for i in range(tau_lin):
+            tau.append(i)
+        factor = 1
+        while len(tau) < length:
+            p = tau[-1] + factor * 1
+            for i in range(0, tau_lin, 2):
+                tau.append(p + factor * i)
+            factor *= 2
+        return time_step * np.array(tau[:length])
+
+    def test_square_distance_componentwise(self):
         s = self.system
-        s.box_l = [10, 10, 10]
-        s.cell_system.skin = 0.4
-        # s.periodicity=0,0,0
-        s.time_step = 0.01
-        s.thermostat.turn_off()
         s.part.add(id=0, pos=(0, 0, 0), v=(1, 2, 3))
 
-        O = espressomd.observables.ParticlePositions(ids=(0,))
-        C2 = espressomd.accumulators.Correlator(
-            obs1=O, tau_lin=10, tau_max=10.0, delta_N=1,
+        obs = espressomd.observables.ParticlePositions(ids=(0,))
+        acc = espressomd.accumulators.Correlator(
+            obs1=obs, tau_lin=10, tau_max=10.0, delta_N=1,
             corr_operation="square_distance_componentwise")
 
         s.integrator.run(1000)
-        s.auto_update_accumulators.add(C2)
+        s.auto_update_accumulators.add(acc)
         s.integrator.run(20000)
 
-        corr = C2.result()
+        corr = acc.result()
 
         # Check pickling
-        C_unpickeled = pickle.loads(pickle.dumps(C2))
-        np.testing.assert_array_equal(corr, C_unpickeled.result())
+        acc_unpickeled = pickle.loads(pickle.dumps(acc))
+        np.testing.assert_array_equal(corr, acc_unpickeled.result())
 
+        tau = self.calc_tau(s.time_step, acc.tau_lin, corr.shape[0])
+        np.testing.assert_array_almost_equal(corr[:, 0], tau)
         for i in range(corr.shape[0]):
             t = corr[i, 0]
             self.assertAlmostEqual(corr[i, 2], t * t, places=3)
             self.assertAlmostEqual(corr[i, 3], 4 * t * t, places=3)
             self.assertAlmostEqual(corr[i, 4], 9 * t * t, places=3)
 
+    def test_tensor_product(self):
+        s = self.system
+        v = np.array([1, 2, 3])
+        s.part.add(id=0, pos=(0, 0, 0), v=v)
+
+        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        acc = espressomd.accumulators.Correlator(
+            obs1=obs, tau_lin=10, tau_max=10.0, delta_N=1,
+            corr_operation="tensor_product")
+
+        s.integrator.run(1000)
+        s.auto_update_accumulators.add(acc)
+        s.integrator.run(20000)
+
+        corr = acc.result()
+
+        # Check pickling
+        acc_unpickeled = pickle.loads(pickle.dumps(acc))
+        np.testing.assert_array_equal(corr, acc_unpickeled.result())
+
+        tau = self.calc_tau(s.time_step, acc.tau_lin, corr.shape[0])
+        np.testing.assert_array_almost_equal(corr[:, 0], tau)
+        for i in range(corr.shape[0]):
+            np.testing.assert_array_almost_equal(corr[i, 2:], np.kron(v, v))
+
+    def test_componentwise_product(self):
+        s = self.system
+        v = np.array([1, 2, 3])
+        s.part.add(id=0, pos=(0, 0, 0), v=v)
+
+        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        acc = espressomd.accumulators.Correlator(
+            obs1=obs, tau_lin=10, tau_max=10.0, delta_N=1,
+            corr_operation="componentwise_product")
+
+        s.integrator.run(1000)
+        s.auto_update_accumulators.add(acc)
+        s.integrator.run(20000)
+
+        corr = acc.result()
+
+        # Check pickling
+        acc_unpickeled = pickle.loads(pickle.dumps(acc))
+        np.testing.assert_array_equal(corr, acc_unpickeled.result())
+
+        tau = self.calc_tau(s.time_step, acc.tau_lin, corr.shape[0])
+        np.testing.assert_array_almost_equal(corr[:, 0], tau)
+        for i in range(corr.shape[0]):
+            np.testing.assert_array_almost_equal(corr[i, 2:], v**2)
+
+    def test_scalar_product(self):
+        s = self.system
+        v = np.array([1, 2, 3])
+        s.part.add(id=0, pos=(0, 0, 0), v=v)
+
+        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        acc = espressomd.accumulators.Correlator(
+            obs1=obs, tau_lin=10, tau_max=10.0, delta_N=1,
+            corr_operation="scalar_product")
+
+        s.integrator.run(1000)
+        s.auto_update_accumulators.add(acc)
+        s.integrator.run(20000)
+
+        corr = acc.result()
+
+        # Check pickling
+        acc_unpickeled = pickle.loads(pickle.dumps(acc))
+        np.testing.assert_array_equal(corr, acc_unpickeled.result())
+
+        tau = self.calc_tau(s.time_step, acc.tau_lin, corr.shape[0])
+        np.testing.assert_array_almost_equal(corr[:, 0], tau)
+        for i in range(corr.shape[0]):
+            np.testing.assert_array_almost_equal(corr[i, 2:], np.sum(v**2))
+
+    def test_correlator_interface(self):
+        # test setters and getters
+        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        acc = espressomd.accumulators.Correlator(
+            obs1=obs, tau_lin=10, tau_max=12.0, delta_N=1,
+            corr_operation="scalar_product")
+        # check tau_lin
+        self.assertEqual(acc.tau_lin, 10)
+        # check tau_max
+        self.assertEqual(acc.tau_max, 12.)
+        # check delta_N
+        self.assertEqual(acc.delta_N, 1)
+        acc.delta_N = 2
+        self.assertEqual(acc.delta_N, 2)
+        # check corr_operation
+        self.assertEqual(acc.corr_operation, "scalar_product")
+
 
 if __name__ == "__main__":
     ut.main()

testsuite/python/observable_cylindrical.py

@@ -46,6 +46,9 @@ class TestCylindricalObservable(ut.TestCase):
         'max_z': 5.0,
     }
 
+    def tearDown(self):
+        self.system.part.clear()
+
     def swap_axis(self, arr, axis):
         if axis == 'x':
             arr = np.dot(
@@ -208,6 +211,70 @@ def test_hist_z(self):
         self.flux_density_profile_test()
         self.density_profile_test()
 
+    def test_cylindrical_pid_profile_interface(self):
+        # test setters and getters
+        params = self.params.copy()
+        params['n_r_bins'] = 4
+        params['n_phi_bins'] = 6
+        params['n_z_bins'] = 8
+        params['ids'] = [0, 1]
+        params['axis'] = [0.0, 1.0, 0.0]
+        self.system.part.add(id=0, pos=[0, 0, 0], type=0)
+        self.system.part.add(id=1, pos=[0, 0, 0], type=1)
+        observable = espressomd.observables.CylindricalDensityProfile(**params)
+        # check pids
+        self.assertEqual(observable.ids, params['ids'])
+        new_pids = [params['ids'][0]]
+        observable.ids = new_pids
+        self.assertEqual(observable.ids, new_pids)
+        # check bins
+        self.assertEqual(observable.n_r_bins, params['n_r_bins'])
+        self.assertEqual(observable.n_phi_bins, params['n_phi_bins'])
+        self.assertEqual(observable.n_z_bins, params['n_z_bins'])
+        obs_data = observable.calculate()
+        np.testing.assert_array_equal(obs_data.shape, [4, 6, 8])
+        observable.n_r_bins = 1
+        observable.n_phi_bins = 2
+        observable.n_z_bins = 3
+        self.assertEqual(observable.n_r_bins, 1)
+        self.assertEqual(observable.n_phi_bins, 2)
+        self.assertEqual(observable.n_z_bins, 3)
+        obs_data = observable.calculate()
+        np.testing.assert_array_equal(obs_data.shape, [1, 2, 3])
+        # check edges lower corner
+        self.assertEqual(observable.min_r, params['min_r'])
+        self.assertEqual(observable.min_phi, params['min_phi'])
+        self.assertEqual(observable.min_z, params['min_z'])
+        observable.min_r = 4
+        observable.min_phi = 5
+        observable.min_z = 6
+        self.assertEqual(observable.min_r, 4)
+        self.assertEqual(observable.min_phi, 5)
+        self.assertEqual(observable.min_z, 6)
+        obs_bin_edges = observable.bin_edges()
+        np.testing.assert_array_equal(obs_bin_edges[0, 0, 0], [4, 5, 6])
+        # check edges upper corner
+        self.assertEqual(observable.max_r, params['max_r'])
+        self.assertEqual(observable.max_phi, params['max_phi'])
+        self.assertEqual(observable.max_z, params['max_z'])
+        observable.max_r = 7
+        observable.max_phi = 8
+        observable.max_z = 9
+        self.assertEqual(observable.max_r, 7)
+        self.assertEqual(observable.max_phi, 8)
+        self.assertEqual(observable.max_z, 9)
+        obs_bin_edges = observable.bin_edges()
+        np.testing.assert_array_equal(obs_bin_edges[-1, -1, -1], [7, 8, 9])
+        # check center
+        np.testing.assert_array_equal(
+            np.copy(observable.center), params['center'])
+        observable.center = [3, 2, 1]
+        np.testing.assert_array_equal(np.copy(observable.center), [3, 2, 1])
+        # check axis
+        np.testing.assert_array_equal(np.copy(observable.axis), params['axis'])
+        observable.axis = [6, 5, 4]
+        np.testing.assert_array_equal(np.copy(observable.axis), [6, 5, 4])
+
 
 if __name__ == "__main__":
     ut.main()

testsuite/python/observable_cylindricalLB.py

@@ -61,6 +61,9 @@ class CylindricalLBObservableCommon:
         'max_z': 5.0,
     }
 
+    def tearDown(self):
+        self.system.part.clear()
+
     def swap_axis(self, arr, axis):
         if axis == 'x':
             arr = np.dot(tests_common.rotation_matrix(
@@ -214,6 +217,69 @@ def test_z_axis(self):
         self.params['axis'] = 'z'
         self.perform_tests()
 
+    def test_cylindrical_lb_profile_interface(self):
+        # test setters and getters
+        params = self.params.copy()
+        params['n_r_bins'] = 4
+        params['n_phi_bins'] = 6
+        params['n_z_bins'] = 8
+        params['axis'] = [0.0, 1.0, 0.0]
+        params['sampling_density'] = 2
+        del params['ids']
+        observable = espressomd.observables.CylindricalLBVelocityProfile(
+            **params)
+        # check bins
+        self.assertEqual(observable.n_r_bins, params['n_r_bins'])
+        self.assertEqual(observable.n_phi_bins, params['n_phi_bins'])
+        self.assertEqual(observable.n_z_bins, params['n_z_bins'])
+        obs_data = observable.calculate()
+        np.testing.assert_array_equal(obs_data.shape, [4, 6, 8, 3])
+        observable.n_r_bins = 1
+        observable.n_phi_bins = 2
+        observable.n_z_bins = 3
+        self.assertEqual(observable.n_r_bins, 1)
+        self.assertEqual(observable.n_phi_bins, 2)
+        self.assertEqual(observable.n_z_bins, 3)
+        obs_data = observable.calculate()
+        np.testing.assert_array_equal(obs_data.shape, [1, 2, 3, 3])
+        # check edges lower corner
+        self.assertEqual(observable.min_r, params['min_r'])
+        self.assertEqual(observable.min_phi, params['min_phi'])
+        self.assertEqual(observable.min_z, params['min_z'])
+        observable.min_r = 4
+        observable.min_phi = 5
+        observable.min_z = 6
+        self.assertEqual(observable.min_r, 4)
+        self.assertEqual(observable.min_phi, 5)
+        self.assertEqual(observable.min_z, 6)
+        obs_bin_edges = observable.bin_edges()
+        np.testing.assert_array_equal(obs_bin_edges[0, 0, 0], [4, 5, 6])
+        # check edges upper corner
+        self.assertEqual(observable.max_r, params['max_r'])
+        self.assertEqual(observable.max_phi, params['max_phi'])
+        self.assertEqual(observable.max_z, params['max_z'])
+        observable.max_r = 7
+        observable.max_phi = 8
+        observable.max_z = 9
+        self.assertEqual(observable.max_r, 7)
+        self.assertEqual(observable.max_phi, 8)
+        self.assertEqual(observable.max_z, 9)
+        obs_bin_edges = observable.bin_edges()
+        np.testing.assert_array_equal(obs_bin_edges[-1, -1, -1], [7, 8, 9])
+        # check center
+        np.testing.assert_array_equal(
+            np.copy(observable.center), params['center'])
+        observable.center = [3, 2, 1]
+        np.testing.assert_array_equal(np.copy(observable.center), [3, 2, 1])
+        # check axis
+        np.testing.assert_array_equal(np.copy(observable.axis), params['axis'])
+        observable.axis = [6, 5, 4]
+        np.testing.assert_array_equal(np.copy(observable.axis), [6, 5, 4])
+        # check sampling_density
+        self.assertEqual(observable.sampling_density, 2)
+        observable.sampling_density = 3
+        self.assertEqual(observable.sampling_density, 3)
+
 
 class CylindricalLBObservableCPU(ut.TestCase, CylindricalLBObservableCommon):