Remove max_seen_particle global from python interface (espressomd#3505)

Description of changes:
 - Added particle_data interface functions to get particle ids
   and maximal id
 - Removed global usage

The rest is pending on the collision detection refactoring, and other stuff...

src/core/particle_data.cpp

@@ -51,9 +51,11 @@
 #include <boost/serialization/vector.hpp>
 #include <boost/variant.hpp>
 
+#include <boost/range/numeric.hpp>
 #include <cmath>
 #include <unordered_map>
 #include <unordered_set>
+#include <utils/keys.hpp>
 /************************************************
  * defines
  ************************************************/
@@ -1465,3 +1467,23 @@ bool particle_exists(int part_id) {
     build_particle_node();
   return particle_node.count(part_id);
 }
+
+std::vector<int> get_particle_ids() {
+  if (particle_node.empty())
+    build_particle_node();
+
+  auto ids = Utils::keys(particle_node);
+  boost::sort(ids);
+
+  return ids;
+}
+
+int get_maximal_particle_id() {
+  if (particle_node.empty())
+    build_particle_node();
+
+  return boost::accumulate(particle_node, -1,
+                           [](int max, const std::pair<int, int> &kv) {
+                             return std::max(max, kv.first);
+                           });
+}

src/core/particle_data.hpp

@@ -625,4 +625,16 @@ bool particle_exists(int part);
  */
 int get_particle_node(int id);
 
+/**
+ * @brief Get all particle ids.
+ *
+ * @return Sorted ids of all existing particles.
+ */
+std::vector<int> get_particle_ids();
+
+/**
+ * @brief Get maximal particle id.
+ */
+int get_maximal_particle_id();
+
 #endif

src/core/reaction_ensemble.cpp

@@ -595,7 +595,7 @@ void ReactionAlgorithm::hide_particle(int p_id, int previous_type) {
  * avoid the id range becoming excessively huge.
  */
 int ReactionAlgorithm::delete_particle(int p_id) {
-  int old_max_seen_id = max_seen_particle;
+  int old_max_seen_id = get_maximal_particle_id();
   if (p_id == old_max_seen_id) {
     // last particle, just delete
     remove_particle(p_id);
@@ -694,7 +694,7 @@ int ReactionAlgorithm::create_particle(int desired_type) {
     p_id = *p_id_iter;
     m_empty_p_ids_smaller_than_max_seen_particle.erase(p_id_iter);
   } else {
-    p_id = max_seen_particle + 1;
+    p_id = get_maximal_particle_id() + 1;
   }
   Utils::Vector3d pos_vec;
 
@@ -1648,10 +1648,10 @@ int WangLandauReactionEnsemble::load_wang_landau_checkpoint(
 }
 
 int ConstantpHEnsemble::get_random_valid_p_id() {
-  int random_p_id = i_random(max_seen_particle + 1);
+  int random_p_id = i_random(get_maximal_particle_id() + 1);
   // draw random p_ids till we draw a pid which exists
   while (not particle_exists(random_p_id))
-    random_p_id = i_random(max_seen_particle + 1);
+    random_p_id = i_random(get_maximal_particle_id() + 1);
   return random_p_id;
 }
 

src/python/espressomd/globals.pxd

@@ -74,7 +74,6 @@ cdef extern from "particle_data.hpp":
 
 cdef extern from "nonbonded_interactions/nonbonded_interaction_data.hpp":
     extern double max_cut
-    extern int max_seen_particle
     extern int max_seen_particle_type
     extern double min_global_cut
     double recalc_maximal_cutoff_bonded()

src/python/espressomd/particle_data.pxd

@@ -194,6 +194,10 @@ cdef extern from "particle_data.hpp":
 
     const particle & get_particle_data(int id) except +
 
+    vector[int] get_particle_ids() except +
+
+    int get_maximal_particle_id()
+
 cdef extern from "virtual_sites.hpp":
     IF VIRTUAL_SITES_RELATIVE == 1:
         void vs_relate_to(int part_num, int relate_to)

src/python/espressomd/particle_data.pyx

@@ -25,7 +25,7 @@ from .interactions import BondedInteraction
 from .interactions import BondedInteractions
 from .interactions cimport bonded_ia_params
 from copy import copy
-from .globals cimport max_seen_particle, max_seen_particle_type, n_part, \
+from .globals cimport max_seen_particle_type, n_part, \
     n_rigidbonds
 import collections
 import functools
@@ -1562,7 +1562,7 @@ cdef class _ParticleSliceImpl:
 
         # We start with a full list of possible particle ids and then
         # remove ids of non-existing particles
-        id_list = np.arange(max_seen_particle + 1, dtype=int)
+        id_list = np.arange(get_maximal_particle_id() + 1, dtype=int)
         id_list = id_list[slice_]
 
         # Generate a mask which will remove ids of non-existing particles
@@ -1818,7 +1818,7 @@ cdef class ParticleList:
         # Handling of particle id
         if "id" not in P:
             # Generate particle id
-            P["id"] = max_seen_particle + 1
+            P["id"] = get_maximal_particle_id() + 1
         else:
             if particle_exists(P["id"]):
                 raise Exception("Particle %d already exists." % P["id"])
@@ -1874,9 +1874,10 @@ Set quat and scalar dipole moment (dipm) instead.")
 
     # Iteration over all existing particles
     def __iter__(self):
-        for i in range(max_seen_particle + 1):
-            if particle_exists(i):
-                yield self[i]
+        ids = get_particle_ids()
+
+        for i in ids:
+            yield self[i]
 
     def exists(self, idx):
         if is_valid_type(idx, int):
@@ -1901,12 +1902,7 @@ Set quat and scalar dipole moment (dipm) instead.")
         remove_all_particles()
 
     def __str__(self):
-        res = ""
-        for i in range(max_seen_particle + 1):
-            if self.exists(i):
-                res += str(self[i]) + ", "
-        # Remove final comma
-        return "ParticleList([" + res[:-2] + "])"
+        return "ParticleList([" + ",".join(get_particle_ids()) + "])"
 
     def writevtk(self, fname, types='all'):
         """
@@ -1975,7 +1971,7 @@ Set quat and scalar dipole moment (dipm) instead.")
         """
 
         def __get__(self):
-            return max_seen_particle
+            return get_maximal_particle_id()
 
     property n_part_types:
         """
@@ -1995,27 +1991,16 @@ Set quat and scalar dipole moment (dipm) instead.")
         def __get__(self):
             return n_rigidbonds
 
-    # property max_part:
-    #     def __get__(self):
-    #         return max_seen_particle
-
-    # # property n_part:
-    #     def __get__(self):
-    #         return n_part
-
     def pairs(self):
         """
         Generator returns all pairs of particles.
 
         """
 
-        cdef int i
-        cdef int j
-        for i in range(max_seen_particle + 1):
-            for j in range(i + 1, max_seen_particle + 1, 1):
-                if not (particle_exists(i) and particle_exists(j)):
-                    continue
+        ids = get_particle_ids()
 
+        for i in ids:
+            for j in ids[i + 1:]:
                 yield (self[i], self[j])
 
     def select(self, *args, **kwargs):