system.part does not behave like an array

[itischler]

choosing a random particle via np.random.choice will fail, if there is only one particle in the system and if this particle has not the id of 0

import espressomd
import numpy as np

sys = espressomd.System(box_l=[10,10,10])
sys.part.add(pos=[4,4,4])
for i in range(10):
  print(np.random.choice(sys.part).id) ### works fine

sys.part.add(pos=[5,5,5])
sys.part.add(pos=[6,6,6])

for i in range(10):
  print(np.random.choice(sys.part).id) ### works fine

sys.part[0].remove()

for i in range(10):
  print(np.random.choice(sys.part).id) ### still works fine

sys.part[2].remove()

print(np.random.choice(sys.part).id) ### fails but particle 1 is still in there

[KaiSzuttor]

what is "failing"?

[itischler]

print(np.random.choice(sys.part).id) ### fails 
File "particle_data.pyx", line 110, in espressomd.particle_data.ParticleHandle.id.__get__
File "particle_data.pyx", line 58, in espressomd.particle_data.ParticleHandle.update_particle_data
RuntimeError: Particle node for id 0 not found!

[KaiSzuttor]

so, your MWE can be made more "minimal":

import espressomd
import numpy as np

sys = espressomd.System(box_l=[10,10,10])
sys.part.add(pos=[4,4,4], id=1)
print(np.random.choice(sys.part))

[KaiSzuttor]

the issue is that the particle list in espresso is not "array-like" because it can have gaps and the index is interpreted as the id (which is not the best idea), right?

[itischler]

yes, and np.random.choice thinks if the length of a list is one the index must be 0

[KaiSzuttor]

okay, but is it a bug or expected behavior? looking at https://espressomd.github.io/doc/espressomd.html#espressomd.particle_data.ParticleList I'd say it's documented behavior

[jhossbach]

so, your MWE can be made more "minimal":

import espressomd
import numpy as np

sys = espressomd.System(box_l=[10,10,10])
sys.part.add(pos=[4,4,4], id=1)
print(np.random.choice(sys.part))

I am using numpy v1.20.3 and I can't reproduce this problem. However running the same snippet with numpy v1.17.4 returns the mentioned error. Seems like this has been resolved somewhere inbetween.

[RudolfWeeber]

This is duck-typing at its worst. If it has [], it must be a list. Unfortunately, dictionaries use [] too.

If we are prepared to break the interface, we can remove the [] operator and instead have
system.part.by_id(id) as a replacement for []
and
system.part.select(predicate) for more complex stuff

[KaiSzuttor]

so, your MWE can be made more "minimal":

import espressomd
import numpy as np

sys = espressomd.System(box_l=[10,10,10])
sys.part.add(pos=[4,4,4], id=1)
print(np.random.choice(sys.part))

I am using numpy v1.20.3 and I can't reproduce this problem. However running the same snippet with numpy v1.17.4 returns the mentioned error. Seems like this has been resolved somewhere inbetween.

it might be that the underlying implementation in numpy changed, but as I said it's not an issue of numpy but an issue of our ParticleList implementation that does not behave like an "array-like" which is the requirement for the numpy function. This is in principle also connected to #3992 and #3369

[KaiSzuttor]

btw: our list is also officially non-pythonic:

Sequences
These represent finite ordered sets indexed by non-negative numbers. The built-in function len() returns the number of items
of a sequence. When the length of a sequence is n, the index set contains the numbers 0, 1, …, n-1. Item i of sequence a is
selected by a[i].

... and

Mutable sequences
[...]
Lists
The items of a list are arbitrary Python objects. Lists are formed by placing a comma-separated list of expressions in square
brackets. (Note that there are no special cases needed to form lists of length 0 or 1.)

https://docs.python.org/3/reference/datamodel.html

[RudolfWeeber]

btw: our list is also officially non-pythonic:

When comparing to Py built-in types System.part should be thought of a dict, I guess.

[KaiSzuttor]

btw: our list is also officially non-pythonic:
When comparing to Py built-in types System.part should be thought of a dict, I guess.

That's why we call it ParticleList? 😃

[KaiSzuttor]

btw: our list is also officially non-pythonic:
When comparing to Py built-in types System.part should be thought of a dict, I guess.

Python dictionaries also do not support slicing.

[jngrad]

This is an issue in the OpenGL visualizer, which relies on particle ids to index arrays. Thanks @mebrito for finding this out!

Here is how the list of bonds is serialized:

espresso/src/python/espressomd/visualization_opengl.py

Lines 861 to 874 in 9f8072f

	def _update_bonds(self):
	if self.specs['draw_bonds']:
	self.bonds = []
	for i, particle in enumerate(self.system.part):
	for bond in particle.bonds:
	# b[0]: Bond, b[1:] Partners
	bond_type = bond[0].type_number()
	if len(bond) == 4:
	self.bonds.append([i, bond[1], bond_type])
	self.bonds.append([i, bond[2], bond_type])
	self.bonds.append([bond[2], bond[3], bond_type])
	else:
	for bond_partner in bond[1:]:
	self.bonds.append([i, bond_partner, bond_type])

The visualizer breaks as soon as a particle id larger or equal to the number of particles is hit (access out of bounds):

espresso/src/python/espressomd/visualization_opengl.py

Lines 1117 to 1118 in 9f8072f

	x_a = self.particles['pos'][b[0]]
	x_b = self.particles['pos'][b[1]]

[RudolfWeeber]

This is done, to my understanding