Return ragged arrays in Environment Matching Code

ChangeLog.md

@@ -15,6 +15,7 @@ and this project adheres to
 
 ### Removed
 * The `global_search` flag in `freud.environment.EnvironmentCluster`.
+* Zero padding from properties of `freud.environment.EnvironmentCluster` and `freud.environment.EnvironmentMotifMatch`.
 
 ## v2.12.1 -- 2022-12-05
 

cpp/environment/MatchEnv.cc

@@ -14,8 +14,7 @@ namespace freud { namespace environment {
 /*****************
  * EnvDisjoinSet *
  *****************/
-EnvDisjointSet::EnvDisjointSet(unsigned int Np) : rank(std::vector<unsigned int>(Np, 0)), m_max_num_neigh(0)
-{}
+EnvDisjointSet::EnvDisjointSet(unsigned int Np) : rank(std::vector<unsigned int>(Np, 0)) {}
 
 void EnvDisjointSet::merge(const unsigned int a, const unsigned int b,
                            BiMap<unsigned int, unsigned int> vec_map, rotmat3<float>& rotation)
@@ -192,7 +191,7 @@ std::vector<vec3<float>> EnvDisjointSet::getAvgEnv(const unsigned int m)
 {
     bool invalid_ind = true;
 
-    std::vector<vec3<float>> env(m_max_num_neigh, vec3<float>(0.0, 0.0, 0.0));
+    std::vector<vec3<float>> env;
     unsigned int N = 0;
 
     // loop over all the environments in the set
@@ -211,7 +210,15 @@ std::vector<vec3<float>> EnvDisjointSet::getAvgEnv(const unsigned int m)
                 for (unsigned int proper_ind = 0; proper_ind < i.vecs.size(); proper_ind++)
                 {
                     unsigned int relative_ind = i.vec_ind[proper_ind];
-                    env[proper_ind] += i.proper_rot * i.vecs[relative_ind];
+                    vec3<float> proper_vec = i.proper_rot * i.vecs[relative_ind];
+                    if (proper_ind < env.size())
+                    {
+                        env[proper_ind] += proper_vec;
+                    }
+                    else
+                    {
+                        env.push_back(proper_vec);
+                    }
                 }
                 ++N;
                 invalid_ind = false;
@@ -228,10 +235,10 @@ std::vector<vec3<float>> EnvDisjointSet::getAvgEnv(const unsigned int m)
 
     // loop through the vectors in env now, dividing by the total number
     // of contributing particle environments to make an average
-    for (unsigned int n = 0; n < m_max_num_neigh; n++)
+    for (auto& vec : env)
     {
-        vec3<float> normed = env[n] / static_cast<float>(N);
-        env[n] = normed;
+        vec3<float> normed = vec / static_cast<float>(N);
+        vec = normed;
     }
     return env;
 }
@@ -246,17 +253,13 @@ std::vector<vec3<float>> EnvDisjointSet::getIndividualEnv(const unsigned int m)
     }
 
     std::vector<vec3<float>> env;
-    for (unsigned int n = 0; n < m_max_num_neigh; n++)
-    {
-        env.emplace_back(0.0, 0.0, 0.0);
-    }
 
     // loop through the vectors, getting them properly indexed
     // add them to env
     for (unsigned int proper_ind = 0; proper_ind < s[m].vecs.size(); proper_ind++)
     {
         unsigned int relative_ind = s[m].vec_ind[proper_ind];
-        env[proper_ind] += s[m].proper_rot * s[m].vecs[relative_ind];
+        env.push_back(s[m].proper_rot * s[m].vecs[relative_ind]);
     }
 
     return env;
@@ -539,13 +542,8 @@ void EnvironmentCluster::compute(const freud::locality::NeighborQuery* nq,
     {
         Environment ei = buildEnv(&env_nlist, env_num_bonds, env_bond, i, i);
         dj.s.push_back(ei);
-        dj.m_max_num_neigh = std::max(dj.m_max_num_neigh, ei.num_vecs);
-        ;
     }
 
-    // reallocate the m_point_environments array
-    m_point_environments.prepare({Np, dj.m_max_num_neigh});
-
     size_t bond(0);
     // loop through points
     for (unsigned int i = 0; i < Np; i++)
@@ -614,9 +612,11 @@ unsigned int EnvironmentCluster::populateEnv(EnvDisjointSet dj)
 
             // label this particle in m_env_index
             m_env_index[particle_ind] = label_ind;
-            for (unsigned int m = 0; m < part_vecs.size(); m++)
+
+            m_point_environments.emplace_back();
+            for (auto& part_vec : part_vecs)
             {
-                m_point_environments(particle_ind, m) = part_vecs[m];
+                m_point_environments[particle_ind].push_back(part_vec);
             }
             particle_ind++;
         }
@@ -654,18 +654,6 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq,
     // because we're inserting the motif into it.
     EnvDisjointSet dj(Np + 1);
 
-    // The NeighborList may contain different numbers of neighbors for each particle, so
-    // we must determine the maximum programmatically to ensure that the disjoint set
-    // operations always allocate enough memory for the largest possible local environment.
-    auto counts = nlist.getCounts();
-    auto* begin = counts.get();
-    auto* end = begin + counts.size();
-    auto max_num_neigh = *std::max_element(begin, end);
-    dj.m_max_num_neigh = max_num_neigh;
-
-    // reallocate the m_point_environments array
-    m_point_environments.prepare({Np, max_num_neigh});
-
     // create the environment characterized by motif. Index it as 0.
     // set the IGNORE flag to true, since this is not an environment we have
     // actually encountered in the simulation.
@@ -714,9 +702,10 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq,
         // grab the set of vectors that define this individual environment
         std::vector<vec3<float>> part_vecs = dj.getIndividualEnv(dummy);
 
-        for (unsigned int m = 0; m < part_vecs.size(); m++)
+        m_point_environments.emplace_back();
+        for (auto& part_vec : part_vecs)
         {
-            m_point_environments(i, m) = part_vecs[m];
+            m_point_environments[i].push_back(part_vec);
         }
     }
 }
@@ -738,10 +727,6 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq,
     // this has to have ONE MORE environment than there are actual particles,
     // because we're inserting the motif into it.
     EnvDisjointSet dj(Np + 1);
-    dj.m_max_num_neigh = motif_size;
-
-    // reallocate the m_point_environments array
-    m_point_environments.prepare({Np, motif_size});
 
     // create the environment characterized by motif. Index it as 0.
     // set the IGNORE flag to true, since this is not an environment we
@@ -797,9 +782,10 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq,
         // grab the set of vectors that define this individual environment
         std::vector<vec3<float>> part_vecs = dj.getIndividualEnv(dummy);
 
-        for (unsigned int m = 0; m < part_vecs.size(); m++)
+        m_point_environments.emplace_back();
+        for (auto& part_vec : part_vecs)
         {
-            m_point_environments(i, m) = part_vecs[m];
+            m_point_environments[i].push_back(part_vec);
         }
     }
 }

cpp/environment/MatchEnv.h

@@ -88,7 +88,6 @@ struct EnvDisjointSet
 
     std::vector<Environment> s;     //!< The disjoint set data
     std::vector<unsigned int> rank; //!< The rank of each tree in the set
-    unsigned int m_max_num_neigh;   //!< The maximum number of neighbors in any environment in the set
 };
 
 /*****************************************************************************
@@ -209,14 +208,14 @@ class MatchEnv
                                 unsigned int i, unsigned int env_ind);
 
     //! Returns the entire Np by m_num_neighbors by 3 matrix of all environments for all particles
-    const util::ManagedArray<vec3<float>>& getPointEnvironments()
+    const std::vector<std::vector<vec3<float>>>& getPointEnvironments()
     {
         return m_point_environments;
     }
 
 protected:
-    util::ManagedArray<vec3<float>> m_point_environments; //!< m_NP by m_max_num_neighbors by 3 matrix of all
-                                                          //!< environments for all particles
+    //!< number of particles x env size x 3 matrix to hold all particle environments
+    std::vector<std::vector<vec3<float>>> m_point_environments;
 };
 
 //! Cluster particles with similar environments.

doc/source/reference/credits.rst

@@ -338,6 +338,7 @@ Tommy Waltmann
 * Extending plotting options for the ``Voronoi`` module.
 * Work with Bradley Dice on adding neighbor vectors to ``freud.locality.NeighborList``.
 * Remove `global_search` flag in ``freud.environment.EnvironmentCluster``.
+* Remove zero-pading from arrays in ``freud.environment.EnvironmentCluster`` and ``freud.environment.EnvironmentMotifMatch``.
 
 Maya Martirossyan
 

freud/_environment.pxd

@@ -71,7 +71,7 @@ cdef extern from "MatchEnv.h" namespace "freud::environment":
 
     cdef cppclass MatchEnv:
         MatchEnv() except +
-        const freud.util.ManagedArray[vec3[float]] &getPointEnvironments()
+        vector[vector[vec3[float]]] &getPointEnvironments()
 
     cdef cppclass EnvironmentMotifMatch(MatchEnv):
         EnvironmentMotifMatch() except +

freud/environment.pyx

@@ -500,10 +500,10 @@ cdef class _MatchEnv(_PairCompute):
     @_Compute._computed_property
     def point_environments(self):
         """:math:`\\left(N_{points}, N_{neighbors}, 3\\right)`
-        :class:`numpy.ndarray`: All environments for all points."""
-        return freud.util.make_managed_numpy_array(
-            &self.matchptr.getPointEnvironments(),
-            freud.util.arr_type_t.FLOAT, 3)
+        list[:class:`numpy.ndarray`]: All environments for all points."""
+        envs = self.matchptr.getPointEnvironments()
+        return [np.asarray([[p.x, p.y, p.z] for p in env])
+                for env in envs]
 
     def __repr__(self):
         return ("freud.environment.{cls}()").format(
@@ -672,7 +672,7 @@ cdef class EnvironmentCluster(_MatchEnv):
     @_Compute._computed_property
     def cluster_environments(self):
         """:math:`\\left(N_{clusters}, N_{neighbors}, 3\\right)`
-        :class:`numpy.ndarray`): The environments for all clusters."""
+        list[:class:`numpy.ndarray`]: The environments for all clusters."""
         envs = self.thisptr.getClusterEnvironments()
         return [np.asarray([[p.x, p.y, p.z] for p in env])
                 for env in envs]
@@ -774,14 +774,8 @@ cdef class EnvironmentMotifMatch(_MatchEnv):
         motif = freud.util._convert_array(motif, shape=(None, 3))
         if (motif == 0).all(axis=1).any():
             warnings.warn(
-                "You are attempting to match a motif containing the zero "
-                "vector. This is likely because you are providing a motif "
-                "created using EnvironmentCluster, which ensures that all "
-                "computed environments are the same size by padding smaller "
-                "environments with the 0 vector. You should remove these "
-                "vectors from the motif before running using "
-                "EnvironmentMotifMatch.compute, since these zero vectors are "
-                "likely to prevent you from finding a match.",
+                "Attempting to match a motif containing the zero "
+                "vector is likely to result in zero matches.",
                 RuntimeWarning
             )
 

tests/test_environment_MatchEnv.py

@@ -11,6 +11,13 @@
 matplotlib.use("agg")
 
 
+def assert_ragged_array(arr):
+    """Assert the given array is a list of numpy arrays."""
+    assert type(arr) == list
+    for a in arr:
+        assert type(a) == np.ndarray
+
+
 class TestCluster:
     # by Chrisy
     test_folder = os.path.join(os.path.dirname(__file__), "numpy_test_files")
@@ -204,6 +211,17 @@ def test_cluster_hardr(self):
             err_msg="BCC Cluster Environment fail",
         )
 
+    def test_ragged_properties(self):
+        """Assert that some properties are returned as ragged arrays."""
+        N = 100
+        L = 10
+        sys = freud.data.make_random_system(L, N)
+        env_cluster = freud.environment.EnvironmentCluster()
+        qargs = dict(r_max=2.0)  # Using r_max ensures different env sizes
+        env_cluster.compute(sys, threshold=0.8, neighbors=qargs)
+        assert_ragged_array(env_cluster.point_environments)
+        assert_ragged_array(env_cluster.cluster_environments)
+
     def _make_global_neighborlist(self, box, points):
         """Get neighborlist where all particles are neighbors."""
         # pairwise distances after wrapping
@@ -474,6 +492,22 @@ def test_warning_motif_zeros(self):
         with pytest.warns(RuntimeWarning):
             match.compute((box, motif), motif, 0.1, neighbors=query_args)
 
+    def test_ragged_properties(self):
+        """Assert that some properties are returned as ragged arrays."""
+        N = 100
+        L = 10
+
+        # sc motif
+        motif = np.array(
+            [[0, 1, 0], [1, 0, 0], [0, 0, 1], [0, 0, -1], [-1, 0, 0], [0, -1, 0]]
+        )
+
+        sys = freud.data.make_random_system(L, N)
+        env_mm = freud.environment.EnvironmentMotifMatch()
+        qargs = dict(r_max=2.0)  # Using r_max ensures different env sizes
+        env_mm.compute(sys, motif, threshold=0.8, neighbors=qargs)
+        assert_ragged_array(env_mm.point_environments)
+
 
 class TestEnvironmentRMSDMinimizer:
     def test_api(self):