Merge pull request #494 from glotzerlab/feature/wont-you-be-my-neighb…

…orquery

Use NeighborQuery in more classes

ChangeLog.md

@@ -19,6 +19,7 @@ and this project adheres to
 * C++ BondHistogramCompute class encapsulates logic of histogram-based methods.
 * NeighborLists and query arguments are now accepted on equal footing by compute methods that involve neighbor finding.
 * 2D PMFTs accept quaternions as well as angles for their orientations.
+* Added NeighborQuery support to ClusterProperties, GaussianDensity, Voronoi.
 
 ### Changed
 * All compute objects that perform neighbor computations now use NeighborQuery internally.

cpp/cluster/Cluster.h

@@ -6,7 +6,6 @@
 
 #include <vector>
 
-#include "Box.h"
 #include "ManagedArray.h"
 #include "NeighborList.h"
 #include "NeighborQuery.h"

cpp/cluster/ClusterProperties.cc

@@ -4,28 +4,28 @@
 #include <vector>
 
 #include "ClusterProperties.h"
+#include "NeighborComputeFunctional.h"
 
 /*! \file ClusterProperties.cc
     \brief Routines for computing properties of point clusters.
 */
 
 namespace freud { namespace cluster {
 
-/*! \param box Box containing the points
-    \param points Positions of the points making up the clusters
+/*! \param nq NeighborQuery containing the points making up the clusters
     \param cluster_idx Index of which cluster each point belongs to
-    \param Np Number of points (length of \a points and \a cluster_idx)
 
-    computeClusterProperties loops over all points in the given array and
-    determines the center of mass of the cluster as well as the gyration
-    tensor. These can be accessed after the call to compute with
-    getClusterCenters() and getClusterGyrations().
+    compute loops over all points in the given array and determines the center
+    of mass of the cluster as well as the gyration tensor. These can be
+    accessed after the call to compute with getClusterCenters() and
+    getClusterGyrations().
 */
-void ClusterProperties::compute(const box::Box& box, const vec3<float>* points,
-                                const unsigned int* cluster_idx, unsigned int Np)
+
+void ClusterProperties::compute(const freud::locality::NeighborQuery* nq,
+                                const unsigned int* cluster_idx)
 {
     // determine the number of clusters
-    const unsigned int* max_cluster_id = std::max_element(cluster_idx, cluster_idx + Np);
+    const unsigned int* max_cluster_id = std::max_element(cluster_idx, cluster_idx + nq->getNPoints());
     const unsigned int num_clusters = *max_cluster_id + 1;
 
     // allocate memory for the cluster properties and temporary arrays
@@ -37,32 +37,30 @@ void ClusterProperties::compute(const box::Box& box, const vec3<float>* points,
     // ref_pos is the first point found in a cluster, it is used as a reference
     // to compute the center in relation to, for handling of the periodic
     // boundary conditions
-    std::vector<vec3<float>> ref_pos(num_clusters, vec3<float>(0.0f, 0.0f, 0.0f));
+    std::vector<vec3<float>> ref_pos(num_clusters);
     // determine if we have seen this cluster before or not (used to initialize ref_pos)
     std::vector<bool> cluster_seen(num_clusters, false);
 
     // Start by determining the center of mass of each cluster. Since we are
     // given an array of points, the easiest way to do this is to loop over
     // all points and add the appropriate information to m_cluster_centers as
     // we go.
-    for (unsigned int i = 0; i < Np; i++)
+    for (unsigned int i = 0; i < nq->getNPoints(); i++)
     {
-        unsigned int c = cluster_idx[i];
-        vec3<float> pos = points[i];
+        const unsigned int c = cluster_idx[i];
 
         // The first time we see the cluster, mark a reference position
         if (!cluster_seen[c])
         {
-            ref_pos[c] = pos;
+            ref_pos[c] = (*nq)[i];
             cluster_seen[c] = true;
         }
 
         // To compute the center in periodic boundary conditions, compute all
         // reference vectors as wrapped vectors relative to ref_pos. When we
         // are done, we can add the computed center to ref_pos to get the
         // center in the space frame.
-        vec3<float> delta = pos - ref_pos[c];
-        delta = box.wrap(delta);
+        const vec3<float> delta(bondVector(locality::NeighborBond(c, i), nq, ref_pos.data()));
 
         // Add the vector into the center tally so far
         m_cluster_centers[c] += delta;
@@ -76,16 +74,16 @@ void ClusterProperties::compute(const box::Box& box, const vec3<float>* points,
     {
         float s = float(m_cluster_sizes[c]);
         vec3<float> v = m_cluster_centers[c] / s + ref_pos[c];
-        m_cluster_centers[c] = box.wrap(v);
+        m_cluster_centers[c] = nq->getBox().wrap(v);
     }
 
     // Now that we have determined the centers of mass for each cluster, tally
     // up the gyration tensor. This has to be done in a loop over the points.
-    for (unsigned int i = 0; i < Np; i++)
+    for (unsigned int i = 0; i < nq->getNPoints(); i++)
     {
         unsigned int c = cluster_idx[i];
-        vec3<float> pos = points[i];
-        vec3<float> delta = box.wrap(pos - m_cluster_centers[c]);
+        vec3<float> pos = (*nq)[i];
+        vec3<float> delta = nq->getBox().wrap(pos - m_cluster_centers[c]);
 
         // get the start pointer for our 3x3 matrix
         m_cluster_gyrations(c, 0, 0) += delta.x * delta.x;

cpp/cluster/ClusterProperties.h

@@ -5,7 +5,7 @@
 #define CLUSTER_PROPERTIES_H
 
 #include "ManagedArray.h"
-#include "Box.h"
+#include "NeighborQuery.h"
 
 /*! \file ClusterProperties.h
     \brief Routines for computing properties of point clusters.
@@ -32,8 +32,8 @@ class ClusterProperties
     ClusterProperties() {}
 
     //! Compute properties of the point clusters
-    void compute(const box::Box& box, const vec3<float>* points,
-            const unsigned int* cluster_idx, unsigned int Np);
+    void compute(const freud::locality::NeighborQuery* nq,
+            const unsigned int* cluster_idx);
 
     //! Get a reference to the last computed cluster centers
     const util::ManagedArray<vec3<float>> &getClusterCenters() const

cpp/density/GaussianDensity.cc

@@ -36,8 +36,10 @@ vec3<unsigned int> GaussianDensity::getWidth()
 //! internal
 /*! \brief Function to compute the density array
  */
-void GaussianDensity::compute(const box::Box& box, const vec3<float>* points, unsigned int n_points)
+void GaussianDensity::compute(const freud::locality::NeighborQuery* nq)
 {
+    auto box = nq->getBox();
+    auto n_points = nq->getNPoints();
     m_box = box;
 
     vec3<unsigned int> width(m_width);
@@ -69,10 +71,11 @@ void GaussianDensity::compute(const box::Box& box, const vec3<float>* points, un
         // for each reference point
         for (size_t idx = begin; idx < end; ++idx)
         {
+            const vec3<float> point = (*nq)[idx];
             // Find which bin the particle is in
-            int bin_x = int((points[idx].x + lx / 2.0f) / grid_size_x);
-            int bin_y = int((points[idx].y + ly / 2.0f) / grid_size_y);
-            int bin_z = int((points[idx].z + lz / 2.0f) / grid_size_z);
+            int bin_x = int((point.x + lx / 2.0f) / grid_size_x);
+            int bin_y = int((point.y + ly / 2.0f) / grid_size_y);
+            int bin_z = int((point.z + lz / 2.0f) / grid_size_z);
 
             // In 2D, only loop over the z=0 plane
             if (m_box.is2D())
@@ -83,16 +86,16 @@ void GaussianDensity::compute(const box::Box& box, const vec3<float>* points, un
             // Only evaluate over bins that are within the cutoff
             for (int k = bin_z - bin_cut_z; k <= bin_z + bin_cut_z; k++)
             {
-                const float dz = float((grid_size_z * k + grid_size_z / 2.0f) - points[idx].z - lz / 2.0f);
+                const float dz = float((grid_size_z * k + grid_size_z / 2.0f) - point.z - lz / 2.0f);
 
                 for (int j = bin_y - bin_cut_y; j <= bin_y + bin_cut_y; j++)
                 {
-                    const float dy = float((grid_size_y * j + grid_size_y / 2.0f) - points[idx].y - ly / 2.0f);
+                    const float dy = float((grid_size_y * j + grid_size_y / 2.0f) - point.y - ly / 2.0f);
 
                     for (int i = bin_x - bin_cut_x; i <= bin_x + bin_cut_x; i++)
                     {
                         // Calculate the distance from the particle to the grid cell
-                        const float dx = float((grid_size_x * i + grid_size_x / 2.0f) - points[idx].x - lx / 2.0f);
+                        const float dx = float((grid_size_x * i + grid_size_x / 2.0f) - point.x - lx / 2.0f);
                         vec3<float> delta = m_box.wrap(vec3<float>(dx, dy, dz));
 
                         const float r_sq = dot(delta, delta);

cpp/density/GaussianDensity.h

@@ -5,9 +5,10 @@
 #define GAUSSIAN_DENSITY_H
 
 #include "Box.h"
+#include "ManagedArray.h"
+#include "NeighborQuery.h"
 #include "ThreadStorage.h"
 #include "VectorMath.h"
-#include "ManagedArray.h"
 
 /*! \file GaussianDensity.h
     \brief Routines for computing Gaussian smeared densities from points.
@@ -48,7 +49,7 @@ class GaussianDensity
     }
 
     //! Compute the Density
-    void compute(const box::Box& box, const vec3<float>* points, unsigned int n_points);
+    void compute(const freud::locality::NeighborQuery* nq);
 
     //! Get a reference to the last computed Density
     const util::ManagedArray<float> &getDensity() const;

cpp/locality/Voronoi.cc

@@ -15,8 +15,11 @@
 namespace freud { namespace locality {
 
 // Voronoi calculations should be kept in double precision.
-void Voronoi::compute(const box::Box &box, const vec3<double>* points, unsigned int n_points)
+void Voronoi::compute(const freud::locality::NeighborQuery* nq)
     {
+        auto box = nq->getBox();
+        auto n_points = nq->getNPoints();
+
         m_polytopes.resize(n_points);
         m_volumes.prepare(n_points);
 
@@ -43,7 +46,8 @@ void Voronoi::compute(const box::Box &box, const vec3<double>* points, unsigned
         );
 
         for (size_t query_point_id = 0; query_point_id < n_points; query_point_id++) {
-            container.put(query_point_id, points[query_point_id].x, points[query_point_id].y, points[query_point_id].z);
+            auto point = (*nq)[query_point_id];
+            container.put(query_point_id, double(point.x), double(point.y), double(point.z));
         }
 
         voro::voronoicell_neighbor cell;

cpp/locality/Voronoi.h

@@ -6,8 +6,9 @@
 
 #include "Box.h"
 #include "ManagedArray.h"
-#include "VectorMath.h"
 #include "NeighborList.h"
+#include "NeighborQuery.h"
+#include "VectorMath.h"
 #include <voro++/src/voro++.hh>
 
 namespace freud { namespace locality {
@@ -20,7 +21,7 @@ class Voronoi
     {
     }
 
-    void compute(const box::Box &box, const vec3<double>* points, unsigned int n_points);
+    void compute(const freud::locality::NeighborQuery* nq);
 
     std::shared_ptr<NeighborList> getNeighborList() const
     {

doc/source/credits.rst

@@ -124,6 +124,7 @@ Bradley Dice - **Lead developer**
 * Added checks to C++ for 2D boxes where required.
 * Refactored cluster module.
 * Standardized vector directionality in computes.
+* NeighborQuery support to ClusterProperties, GaussianDensity, Voronoi.
 * Standardized APIs for order parameters.
 
 Richmond Newman

freud/_cluster.pxd

@@ -21,11 +21,8 @@ cdef extern from "Cluster.h" namespace "freud::cluster":
 cdef extern from "ClusterProperties.h" namespace "freud::cluster":
     cdef cppclass ClusterProperties:
         ClusterProperties()
-        void compute(
-            const freud._box.Box &,
-            const vec3[float]*,
-            const unsigned int*,
-            unsigned int) except +
+        void compute(const freud._locality.NeighborQuery*,
+                     const unsigned int*) except +
         const freud.util.ManagedArray[vec3[float]] &getClusterCenters() const
         const freud.util.ManagedArray[float] &getClusterGyrations() const
         const freud.util.ManagedArray[unsigned int] &getClusterSizes() const

freud/_density.pxd

@@ -28,10 +28,7 @@ cdef extern from "GaussianDensity.h" namespace "freud::density":
         GaussianDensity(vec3[unsigned int], float, float) except +
         const freud._box.Box & getBox() const
         void reset()
-        void compute(
-            const freud._box.Box &,
-            const vec3[float]*,
-            unsigned int) except +
+        void compute(const freud._locality.NeighborQuery*) except +
         const freud.util.ManagedArray[float] &getDensity() const
         vec3[unsigned int] getWidth() const
         float getSigma() const

freud/_locality.pxd

@@ -111,10 +111,7 @@ cdef extern from "AABBQuery.h" namespace "freud::locality":
 cdef extern from "Voronoi.h" namespace "freud::locality":
     cdef cppclass Voronoi:
         Voronoi()
-        void compute(
-            const freud._box.Box &,
-            const vec3[double]*,
-            const unsigned int) nogil except +
+        void compute(const NeighborQuery*) nogil except +
         vector[vector[vec3[double]]] getPolytopes() const
         const freud.util.ManagedArray[double] &getVolumes() const
         shared_ptr[NeighborList] getNeighborList() const

freud/cluster.pyx

@@ -179,7 +179,7 @@ cdef class ClusterProperties(Compute):
     def __dealloc__(self):
         del self.thisptr
 
-    def compute(self, box, points, cluster_idx):
+    def compute(self, neighbor_query, cluster_idx):
         R"""Compute properties of the point clusters.
         Loops over all points in the given array and determines the center of
         mass of the cluster as well as the gyration tensor. These can be
@@ -194,19 +194,14 @@ cdef class ClusterProperties(Compute):
             cluster_idx ((:math:`N_{points}`,) :class:`np.ndarray`):
                 Cluster indexes for each point.
         """
-        cdef freud.box.Box b = freud.util._convert_box(box)
-
-        points = freud.util._convert_array(points, shape=(None, 3))
+        cdef freud.locality.NeighborQuery nq = \
+            freud.locality._make_default_nq(neighbor_query)
         cluster_idx = freud.util._convert_array(
-            cluster_idx, shape=(points.shape[0], ), dtype=np.uint32)
-        cdef const float[:, ::1] l_points = points
+            cluster_idx, shape=(nq.points.shape[0], ), dtype=np.uint32)
         cdef const unsigned int[::1] l_cluster_idx = cluster_idx
-        cdef unsigned int Np = l_points.shape[0]
         self.thisptr.compute(
-            dereference(b.thisptr),
-            <vec3[float]*> &l_points[0, 0],
-            <unsigned int*> &l_cluster_idx[0],
-            Np)
+            nq.get_ptr(),
+            <unsigned int*> &l_cluster_idx[0])
         return self
 
     @Compute._computed_property

freud/density.pyx

@@ -257,7 +257,7 @@ cdef class GaussianDensity(Compute):
     def box(self):
         return freud.box.BoxFromCPP(self.thisptr.getBox())
 
-    def compute(self, box, points):
+    def compute(self, neighbor_query):
         R"""Calculates the Gaussian blur for the specified points. Does not
         accumulate (will overwrite current image).
 
@@ -267,12 +267,9 @@ cdef class GaussianDensity(Compute):
             points ((:math:`N_{points}`, 3) :class:`numpy.ndarray`):
                 Points to calculate the local density.
         """
-        cdef freud.box.Box b = freud.util._convert_box(box)
-        points = freud.util._convert_array(points, shape=(None, 3))
-        cdef const float[:, ::1] l_points = points
-        cdef unsigned int n_p = points.shape[0]
-        self.thisptr.compute(dereference(b.thisptr),
-                             <vec3[float]*> &l_points[0, 0], n_p)
+        cdef freud.locality.NeighborQuery nq = \
+            freud.locality._make_default_nq(neighbor_query)
+        self.thisptr.compute(nq.get_ptr())
         return self
 
     @Compute._computed_property