python: Cast away array_locked on matrix operation

With the new universal functions framework introduced in numpy 1.13,
in-place operations such as `obj1 |= obj2` (with obj1 an np.ndarray
and obj2 an espressomd.utils.array_locked) are now possible by
casting the array_locked to np.ndarray on the fly. For more details,
see https://numpy.org/devdocs/reference/ufuncs.html

doc/sphinx/particles.rst

@@ -68,12 +68,13 @@ Similarly, the position can be set::
 Vectorial properties
 ~~~~~~~~~~~~~~~~~~~~
 
-For vectorial particle properties, component-wise manipulation like ``system.part[0].pos[0]
-= 1`` or in-place operators like ``+=`` or ``*=`` are not allowed and result in an error.
-This behavior is inherited, so the same applies to ``a`` after ``a =
-system.part[0].pos``. If you want to use a vectorial property for further
-calculations, you should explicitly make a copy e.g. via
-``a = numpy.copy(system.part[0].pos)``.
+For vectorial particle properties, component-wise manipulation such as
+``system.part[0].pos[0] = 1`` or in-place operators like ``+=`` or ``*=``
+are not allowed and raise an error. This behavior is inherited, so
+the same applies to ``a`` after ``a = system.part[0].pos`` but not to
+``b`` after ``b = np.array([1. ,1. , 1.]); b *= system.part[0].pos``.
+If you want to use a vectorial property for further calculations, you
+should explicitly make a copy e.g. via ``a = numpy.copy(system.part[0].pos)``.
 
 .. _Interacting with groups of particles:
 

doc/sphinx/system_setup.rst

@@ -12,11 +12,10 @@ The global variables in Python are controlled via the
 :class:`espressomd.system.System` class.
 Global system variables can be read and set in Python simply by accessing the
 attribute of the corresponding Python object. Those variables that are already
-available in the Python interface are listed in the following. Note that for the
-vectorial properties ``box_l`` and ``periodicity``, component-wise manipulation
-like ``system.box_l[0] = 1`` or in-place operators like ``+=`` or ``*=`` are not
-allowed and result in an error. This behavior is inherited, so the same applies
-to ``a`` after ``a = system.box_l``. If you want to use a vectorial property
+available in the Python interface are listed in the following. Note that for
+vectorial properties such as ``box_l`` and ``periodicity``, component-wise
+manipulation and in-place operations are disabled
+(see :ref:`Vectorial properties`). If you want to use a vectorial property
 for further calculations, you should explicitly make a copy e.g. via
 ``a = numpy.copy(system.box_l)``.
 

src/python/espressomd/utils.pyx

@@ -152,17 +152,15 @@ Use numpy.copy(<ESPResSo array property>) to get a writable copy."
         obj.flags.writeable = False
         return obj
 
-    def __add__(self, other):
-        return np.copy(self) + other
-
-    def __radd__(self, other):
-        return other + np.copy(self)
-
-    def __sub__(self, other):
-        return np.copy(self) - other
-
-    def __rsub__(self, other):
-        return other - np.copy(self)
+    def __array_ufunc__(self, ufunc, method, *args, **kwargs):
+        inputs = []
+        for i, input_ in enumerate(args):
+            if isinstance(input_, array_locked):
+                # cast away the array_locked type
+                inputs.append(input_.view(np.ndarray))
+            else:
+                inputs.append(input_)
+        return super().__array_ufunc__(ufunc, method, *inputs, **kwargs)
 
     def __repr__(self):
         return repr(np.array(self))

testsuite/python/array_properties.py

@@ -87,11 +87,27 @@ def test_setter(self):
         array = [4, 5, 6]
         np.testing.assert_array_equal(array, [4, 5, 6])
 
+    def test_in_place_operations(self):
+        array1 = espressomd.utils.array_locked([1, 2, 3])
+        array2 = espressomd.utils.array_locked([1, 2, 3])
+        # in-place operations on a regular numpy array should not raise
+        lhs = np.array([1, 2, 3])
+        lhs += array1
+        lhs -= array1
+        lhs *= array1
+        lhs //= array1
+        np.testing.assert_array_equal(lhs, [1, 2, 3])
+        # numpy testing should not fail (involves in-place logical operations)
+        lhs = np.array([True, True, False])
+        lhs &= espressomd.utils.array_locked([True, False, False])
+        np.testing.assert_array_almost_equal(lhs, [True, False, False])
+        np.testing.assert_array_almost_equal(array1, array2)
+
 
 def check_array_writable(array):
     value = np.random.random(array.shape[0]).astype(type(array[0]))
     array = value
-    np.testing.assert_array_almost_equal(np.copy(array), value)
+    np.testing.assert_array_almost_equal(array, value)
 
 
 class ArrayPropertyTest(ArrayCommon):
@@ -111,6 +127,21 @@ def assert_copy_is_writable(self, array):
         cpy = np.copy(array)
         self.assertTrue(cpy.flags.writeable)
 
+    def assert_immutable(self, array):
+        original_value = array[0]
+        # check a simple operation
+        mutable_copy = 1 * array
+        self.assertTrue(mutable_copy.flags.writeable)
+        mutable_copy[0] = original_value + 1
+        self.assertEqual(mutable_copy[0], original_value + 1)
+        self.assertEqual(array[0], original_value)
+        # check an in-place operation
+        mutable_array = np.core.float_(0.)
+        mutable_array += array
+        mutable_array[0] = original_value + 1
+        self.assertEqual(mutable_array[0], original_value + 1)
+        self.assertEqual(array[0], original_value)
+
     def test_common(self):
         self.assert_operator_usage_raises(self.system.part[0].pos)
         self.assert_operator_usage_raises(self.system.part[0].v)
@@ -119,6 +150,8 @@ def test_common(self):
 
         self.assert_operator_usage_raises(self.system.box_l)
 
+        self.assert_immutable(self.system.box_l)
+
         check_array_writable(self.system.part[0].pos)
         check_array_writable(self.system.part[0].v)
         check_array_writable(self.system.part[0].f)
@@ -131,6 +164,10 @@ def test_common(self):
 
         self.assert_copy_is_writable(self.system.box_l)
 
+        original_box_l = self.system.box_l
+        self.system.box_l = original_box_l + 1
+        np.testing.assert_array_equal(self.system.box_l, original_box_l + 1)
+
     @utx.skipIfMissingFeatures(["ROTATION"])
     def test_rotation(self):
         self.assert_operator_usage_raises(self.system.part[0].omega_lab)