diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd index 19b201896ed..ce75327f4a9 100644 --- a/src/python/espressomd/particle_data.pxd +++ b/src/python/espressomd/particle_data.pxd @@ -227,7 +227,7 @@ cdef class ParticleHandle(object): cdef int update_particle_data(self) except -1 -cdef class ParticleSlice: +cdef class _ParticleSliceImpl: cdef particle particle_data cdef int update_particle_data(self, id) except -1 diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx index daf01216466..fac24edcd38 100644 --- a/src/python/espressomd/particle_data.pyx +++ b/src/python/espressomd/particle_data.pyx @@ -29,6 +29,7 @@ from .interactions import BondedInteractions from copy import copy from globals cimport max_seen_particle, time_step, smaller_time_step, box_l, n_part, n_rigidbonds, n_particle_types import collections +import functools PARTICLE_EXT_FORCE = 1 @@ -1148,10 +1149,11 @@ cdef class ParticleHandle: setattr(self, k, P[k]) -cdef class ParticleSlice: +cdef class _ParticleSliceImpl: """ Handles slice inputs e.g. part[0:2]. Sets values for selected slices or returns values as a single list. + This is the base class that should not be used directly as it lacks some properties. Use :class:`espressomd.ParticleSlice` instead. """ def __cinit__(self, slice_): @@ -1178,59 +1180,12 @@ cdef class ParticleSlice: else: return 0 - # Particle Type - property type: - """ - Particle type. - - """ - def __get__(self): - type_list = [] - for id in self.id_selection: - type_list.append(ParticleHandle(id).type) - return type_list - - def __set__(self, _type_list): - if isinstance(_type_list, int): - for id in self.id_selection: - ParticleHandle(id).type = _type_list - return - if len(self.id_selection) != len(_type_list): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_type_list), len(self.id_selection))) - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[i]).type = _type_list[i] - - # Position - property pos: - """ - Particle position (not folded into central image). - - """ - def __set__(self, _pos_array): - if len(self.id_selection) != len(_pos_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_pos_array), len(self.id_selection))) - - cdef double mypos[3] - for i in range(len(_pos_array)): - ParticleHandle(self.id_selection[i]).pos = _pos_array[i] - - def __get__(self): - pos_array = np.zeros((len(self.id_selection), 3)) - for i in range(len(self.id_selection)): - pos_array[i, :] = ParticleHandle(self.id_selection[i]).pos - return pos_array - property pos_folded: """ Particle position (folded into central image). """ - def __set__(self, d): - raise Exception("setting a folded position is not implemented.") - def __get__(self): pos_array = np.zeros((len(self.id_selection), 3)) for i in range(len(self.id_selection)): @@ -1238,172 +1193,7 @@ cdef class ParticleSlice: self.id_selection[i]).pos_folded return pos_array - # Velocity - property v: - """ - Particle velocity. - - """ - - def __set__(self, _v_array): - if len(np.array(_v_array).shape) == 1: - for id in self.id_selection: - ParticleHandle(id).v = _v_array - return - - if len(self.id_selection) != len(_v_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_v_array), len(self.id_selection))) - - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[i]).v = _v_array[i] - - def __get__(self): - v_array = np.zeros((len(self.id_selection), 3)) - for i in range(len(self.id_selection)): - v_array[i, :] = ParticleHandle(self.id_selection[i]).v - return v_array - - # Force - property f: - """ - Particle force. - - """ - - def __set__(self, _f_array): - if len(np.array(_f_array).shape) == 1: - for id in self.id_selection: - ParticleHandle(id).f = _f_array - return - - if len(self.id_selection) != len(_f_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_f_array), len(self.id_selection))) - for i in range(len(_f_array)): - ParticleHandle(self.id_selection[i]).f = _f_array[i] - - def __get__(self): - f_array = np.zeros((len(self.id_selection), 3)) - for i in range(len(self.id_selection)): - f_array[i, :] = ParticleHandle(self.id_selection[i]).f - return f_array - - property mass: - """ - Particle mass. - - .. note:: - - If not set the particle mass is ``1`` in reduced units. - - """ - - def __set__(self, _mass_array): - IF MASS: - if isinstance(_mass_array, int) or isinstance(_mass_array, float): - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[i]).mass = _mass_array - return - if len(self.id_selection) != len(_mass_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_mass_array), len(self.id_selection))) - for i in range(len(_mass_array)): - ParticleHandle(self.id_selection[i]).mass = _mass_array[i] - ELSE: - raise Exception("You are trying to set the particle mass \ - but the mass feature is not compiled in.") - - def __get__(self): - mass_array = np.zeros_like(self.id_selection) - for i in range(len(self.id_selection)): - mass_array[i] = ParticleHandle(self.id_selection[i]).mass - return mass_array - - IF ELECTROSTATICS == 1: - property q: - """ - Particle charge. - - """ - - def __set__(self, _q_array): - if isinstance(_q_array, int) or isinstance(_q_array, float): - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[i]).q = _q_array - return - - if len(self.id_selection) != len(_q_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_q_array), len(self.id_selection))) - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[i]).q = _q_array[i] - - def __get__(self): - q_array = np.zeros_like(self.id_selection) - for i in range(len(self.id_selection)): - q_array[i] = ParticleHandle(self.id_selection[i]).q - return q_array - - IF EXTERNAL_FORCES: - property ext_force: - """ - External force on a particle defined by a vector. - - """ - - def __set__(self, _ext_f_array): - if len(np.array(_ext_f_array).shape) == 1: - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[ - i]).ext_force = _ext_f_array - return - - if len(self.id_selection) != len(_ext_f_array): - raise Exception("Input list size (%i) does not match slice size (%i)." % ( - len(_ext_f_array), len(self.id_selection))) - - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[ - i]).ext_force = _ext_f_array[i] - - def __get__(self): - ext_f_array = np.zeros((len(self.id_selection), 3)) - for i in range(len(self.id_selection)): - ext_f_array[i, :] = ParticleHandle( - self.id_selection[i]).ext_force - - return ext_f_array - IF EXCLUSIONS: - property exclude: - """ - Exclude particle from interaction. - - """ - - def __set__(self, _partners): - if not isinstance(_partners, list): - raise Exception( - "List object expected for exclusion partners.") - if isinstance(_partners[0], list): - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[ - i]).exclude = _partners[i] - elif isinstance(_partners[0], int): - for i in range(len(self.id_selection)): - ParticleHandle(self.id_selection[ - i]).exclude = _partners - else: - raise TypeError("Unexpected exclusion partner type.") - - def __get__(self): - _exclude_array = [] - for i in range(len(self.id_selection)): - _exclude_array.append(ParticleHandle( - self.id_selection[i]).exclude) - return _exclude_array - def add_exclusion(self, _partners): self.exclude = _partners @@ -1430,9 +1220,9 @@ cdef class ParticleSlice: pl = ParticleList() for i in self.id_selection: if pl.exists(i): - res += str(pl[i]) + "\n" - # Remove final newline - return res[:-1] + res += str(pl[i]) + ", " + # Remove final comma + return "ParticleSlice([" + res[:-2] + "])" def update(self, P): if "id" in P: @@ -1479,6 +1269,13 @@ cdef class ParticleSlice: ParticleHandle(id).remove() +class ParticleSlice(_ParticleSliceImpl): + """ + Handles slice inputs e.g. part[0:2]. Sets values for selected slices or returns values as a single list. + + """ + pass + cdef class ParticleList: """ Provides access to the particles via [i], where i is the particle id. Returns a ParticleHandle object. @@ -1625,7 +1422,7 @@ cdef class ParticleList: yield self[i] def exists(self, idx): - if isinstance(idx, int): + if isinstance(idx, int) or issubclass(type(idx), np.integer): return particle_exists(idx) if isinstance(idx, slice) or isinstance(idx, tuple) or isinstance(idx, list) or isinstance(idx, np.ndarray): tf_array = np.zeros(len(idx), dtype=np.bool) @@ -1640,9 +1437,9 @@ cdef class ParticleList: res = "" for i in range(max_seen_particle + 1): if self.exists(i): - res += str(self[i]) + "\n" - # Remove final newline - return res[:-1] + res += str(self[i]) + ", " + # Remove final comma + return "ParticleList([" + res[:-2] + "])" def writevtk(self, fname, types='all'): """ @@ -1721,3 +1518,61 @@ cdef class ParticleList: if not (self.exists(i) and self.exists(j)): continue yield (self[i], self[j]) + +# auto-add missing particle properties to ParticleSlice +for attribute_name in particle_attributes: + if attribute_name in dir(ParticleSlice): + continue + + def seta(particle_slice, values, attribute): + """Setter function that sets attribute on every member of particle_slice. + If values contains only one element, all members are set to it. If it + contains as many elements as there are members, each of them gets set + to the corresponding one.""" + target = getattr(ParticleHandle(particle_slice.id_selection[0]), attribute) + target_shape = np.shape(target) + N = len(particle_slice.id_selection) + + if not target_shape: # scalar quantity + if not np.shape(values): # one value provided + for i in range(N): + setattr(ParticleHandle(particle_slice.id_selection[i]), attribute, values) + elif np.shape(values)[0] == N: # one value for each particle provided + for i in range(N): + setattr(ParticleHandle(particle_slice.id_selection[i]), attribute, values[i]) + else: + raise Exception("Shape of value (%s) does not broadcast to shape of attribute (%s)." % ( + np.shape(values), target_shape)) + return + + if target_shape == np.shape(values): # one value provided + for i in range(N): + setattr(ParticleHandle(particle_slice.id_selection[i]), attribute, values) + elif target_shape == tuple(np.shape(values)[1:]) and np.shape(values)[0] == N: # one value for each particle provided + for i in range(N): + setattr(ParticleHandle(particle_slice.id_selection[i]), attribute, values[i]) + else: + raise Exception("Shape of value (%s) does not broadcast to shape of attribute (%s)." % ( + np.shape(values), target_shape)) + + def geta(particle_slice, attribute): + """Getter function that copies attribute from every member of particle_slice into an array.""" + N = len(particle_slice.id_selection) + if N == 0: + return np.empty(0, dtype=type(None)) + + target = getattr(ParticleHandle(particle_slice.id_selection[0]), attribute) # get first slice member to determine its type + if type(target) is np.ndarray: # vectorial quantity + target_type = target.dtype + else: # scalar quantity + target_type = type(target) + + values = np.empty((N,) + np.shape(target), dtype=target_type) + for i in range(N): + values[i] = getattr(ParticleHandle(particle_slice.id_selection[i]), attribute) + return values + + # synthesize a new property + new_property = property(functools.partial(geta, attribute=attribute_name), functools.partial(seta, attribute=attribute_name), doc=getattr(ParticleHandle, attribute_name).__doc__) + # attach the property to ParticleSlice + setattr(ParticleSlice, attribute_name, new_property) diff --git a/src/python/espressomd/utils.pyx b/src/python/espressomd/utils.pyx index 34987079793..bf2a8ca2047 100644 --- a/src/python/espressomd/utils.pyx +++ b/src/python/espressomd/utils.pyx @@ -61,7 +61,7 @@ cdef check_type_or_throw_except(x, n, t, msg): if hasattr(x, "__getitem__"): for i in range(len(x)): if not isinstance(x[i], t): - if not (t == float and isinstance(x[i], int)): + if not (t == float and isinstance(x[i], int)) and not (t == int and issubclass(type(x[i]), np.integer)): raise ValueError( msg + " -- Item " + str(i) + " was of type " + type(x[i]).__name__) else: @@ -71,7 +71,7 @@ cdef check_type_or_throw_except(x, n, t, msg): else: # N=1 and a single value if not isinstance(x, t): - if not (t == float and isinstance(x, int)): + if not (t == float and isinstance(x, int)) and not (t == int and issubclass(type(x), np.integer)): raise ValueError(msg + " -- Got an " + type(x).__name__) diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index c6f2bec9d4b..83a9a6f8237 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -15,6 +15,7 @@ set(py_tests bondedInteractions.py observables.py p3m_gpu.py particle.py + particle_slice.py rotational_inertia.py lbgpu_remove_total_momentum.py tabulated.py diff --git a/testsuite/python/observables.py b/testsuite/python/observables.py index 172c259ec69..1cbe92e576d 100644 --- a/testsuite/python/observables.py +++ b/testsuite/python/observables.py @@ -104,16 +104,15 @@ def func(self): com_force =generate_test_for_pid_observable(ComForce,"f","sum") - # Disabled untile ParticleSlice.dip is implemented - #if espressomd.has_features(["DIPOLES"]): - #test_mag_dip =generate_test_for_pid_observable(MagneticDipoleMoment,"dip","sum") + if espressomd.has_features(["DIPOLES"]): + test_mag_dip = generate_test_for_pid_observable(MagneticDipoleMoment,"dip","sum") - # This is disabled until ParticleSlice.omega_body is implemented + # This is disabled as it does not currently work #if espressomd.has_features(["ROTATION"]): - #test_omega_body= generate_test_for_pid_observable(ParticleBodyVelocities,"omega_body") + # test_omega_body = generate_test_for_pid_observable(ParticleBodyVelocities,"omega_body") def test_stress_tensor(self): s=self.es.analysis.stress_tensor()["total"].reshape(9) diff --git a/testsuite/python/particle_slice.py b/testsuite/python/particle_slice.py new file mode 100644 index 00000000000..25a97cdaacb --- /dev/null +++ b/testsuite/python/particle_slice.py @@ -0,0 +1,76 @@ +from __future__ import print_function +import unittest as ut +import espressomd +from espressomd import has_features +import numpy as np + +class ParticleSliceTest(ut.TestCase): + + state = [[0,0,0], [0,0,1]] + system = espressomd.System() + + def __init__(self, *args, **kwargs): + super(ParticleSliceTest, self).__init__(*args, **kwargs) + self.system.part.clear() + self.system.part.add(pos=[0,0,0]) + if has_features(["EXTERNAL_FORCES"]): + self.system.part.add(pos=[0,0,1], fix=self.state[1]) + self.assertTrue( np.array_equal(self.system.part[0].fix, self.state[0]) ) + self.assertTrue( np.array_equal(self.system.part[1].fix, self.state[1]) ) + self.assertTrue( np.array_equal(self.system.part[:].fix, self.state) ) + xs = self.system.part[:].pos + for i in range(len(xs)): + self.assertTrue( np.array_equal(xs[i], self.system.part[i].pos) ) + + @ut.skipIf(not has_features(["EXTERNAL_FORCES"]),"Features not available, skipping test!") + def test_1_set_different_values(self): + self.state[0] = [1,0,0] + self.state[1] = [1,0,0] + self.system.part[:].fix = self.state + self.assertTrue( np.array_equal(self.system.part[:].fix, self.state) ) + + @ut.skipIf(not has_features(["EXTERNAL_FORCES"]),"Features not available, skipping test!") + def test_2_set_same_value(self): + self.state[0] = [0,1,0] + self.state[1] = [0,1,0] + self.system.part[:].fix = self.state[1] + self.assertTrue( np.array_equal(self.system.part[:].fix, self.state) ) + + @ut.skipIf(not has_features(["EXTERNAL_FORCES"]),"Features not available, skipping test!") + def test_3_set_one_value(self): + self.state[1] = [0,0,1] + self.system.part[1:].fix = self.state[1] + self.assertTrue( np.array_equal(self.system.part[:].fix, self.state) ) + + @ut.skipIf(not has_features(["EXTERNAL_FORCES"]),"Features not available, skipping test!") + def test_4_str(self): + self.assertEqual( repr(self.system.part[0].fix), repr(np.array([0,1,0])) ) + self.assertEqual( repr(self.system.part[:].fix), repr(np.array([[0,1,0], [0,0,1]])) ) + + def test_pos_str(self): + self.assertEqual( repr(self.system.part[0].pos), repr(np.array([0.0,0.0,0.0])) ) + if (len(self.system.part)) > 1: + self.assertEqual( repr(self.system.part[:].pos), repr(np.array([[0.0,0.0,0.0], [0.0,0.0,1.0]])) ) + + @ut.skipIf(not has_features(["ELECTROSTATICS"]),"Features not available, skipping test!") + def test_scalar(self): + self.system.part[:1].q = 1.3 + self.assertEqual(self.system.part[0].q, 1.3) + self.system.part[:].q = 2.0 + self.assertEqual(self.system.part[0].q, 2) + self.assertEqual(self.system.part[1].q, 2) + self.system.part[:].q = 3 + self.assertEqual(self.system.part[0].q, 3) + self.assertEqual(self.system.part[1].q, 3) + self.system.part[:].q = [-1,1.0] + self.assertEqual(self.system.part[0].q, -1) + self.assertEqual(self.system.part[1].q, 1) + qs = self.system.part[:].q + self.assertEqual(qs[0], -1) + self.assertEqual(qs[1], 1) + + def test_empty(self): + self.assertTrue( np.array_equal(self.system.part[0:0].pos, np.empty(0)) ) + +if __name__ == "__main__": + ut.main()