diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index 17ba73ba7..835fd477c 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -65,6 +65,7 @@ jobs: pytest tests/test_data.py -v pytest tests/test_pmft.py -v pytest tests/test_util.py -v + pytest tests/test_environment*.py -v pytest tests/test_diffraction*.py -v pytest tests/test_interface.py -v pytest tests/test_msd_msd.py -v diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 787166f26..3b8777978 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -84,7 +84,6 @@ repos: ^extern/| ^freud/density/| ^freud/diffraction/| - ^freud/environment/| ^freud/order/ ) args: diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 08711d4bd..e78725232 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -42,17 +42,17 @@ add_library( diffraction/StaticStructureFactorDirect.h diffraction/StaticStructureFactorDirect.cc # environment - # environment/AngularSeparation.h - # environment/AngularSeparation.cc - # environment/BondOrder.h - # environment/BondOrder.cc - # environment/LocalBondProjection.h - # environment/LocalBondProjection.cc - # environment/LocalDescriptors.h - # environment/LocalDescriptors.cc - # environment/MatchEnv.h - # environment/MatchEnv.cc - # environment/Registration.h + environment/AngularSeparation.h + environment/AngularSeparation.cc + environment/BondOrder.h + environment/BondOrder.cc + environment/LocalBondProjection.h + environment/LocalBondProjection.cc + environment/LocalDescriptors.h + environment/LocalDescriptors.cc + environment/MatchEnv.h + environment/MatchEnv.cc + environment/Registration.h # locality locality/AABB.h locality/AABBQuery.cc @@ -152,15 +152,24 @@ nanobind_add_module( target_link_libraries(_diffraction PUBLIC freud TBB::tbb) target_set_install_rpath(_diffraction) -# environment nanobind_add_module(_environment environment/...) -# target_link_libraries(_environment PUBLIC freud TBB::tbb) -# target_set_install_rpath(_environment) - # box nanobind_add_module(_box box/module-box.cc box/export-Box.cc box/export-Box.h) target_link_libraries(_box PUBLIC TBB::tbb) target_set_install_rpath(_box) +# environment +nanobind_add_module( + _environment + environment/module-environment.cc + environment/export-AngularSeparationNeighbor.cc + environment/export-AngularSeparationGlobal.cc + environment/export-LocalBondProjection.cc + environment/export-LocalDescriptors.cc + environment/export-BondOrder.cc + environment/export-MatchEnv.cc) +target_link_libraries(_environment PUBLIC freud TBB::tbb) +target_set_install_rpath(_environment) + # locality nanobind_add_module( _locality @@ -213,6 +222,7 @@ set(python_files box.py cluster.py data.py + environment.py diffraction.py errors.py locality.py @@ -223,7 +233,7 @@ set(python_files interface.py plot.py util.py) -# cluster.py density.py diffraction.py environment.py interface.py msd.py) +# density.py) copy_files_to_build("${python_files}" "freud" "*.py") @@ -232,9 +242,9 @@ if(SKBUILD) install(FILES ${python_files} DESTINATION freud) install(TARGETS _box DESTINATION freud) install(TARGETS _cluster DESTINATION freud) + # install(TARGETS _density DESTINATION freud) + install(TARGETS _environment DESTINATION freud) install(TARGETS _diffraction DESTINATION freud) - # install(TARGETS _density DESTINATION freud) install(TARGETS _environment - # DESTINATION freud) install(TARGETS _locality DESTINATION freud) install(TARGETS _order DESTINATION freud) install(TARGETS _parallel DESTINATION freud) diff --git a/freud/__init__.py b/freud/__init__.py index 618852f3c..9a8b67dda 100644 --- a/freud/__init__.py +++ b/freud/__init__.py @@ -1,11 +1,13 @@ # Copyright (c) 2010-2024 The Regents of the University of Michigan # This file is from the freud project, released under the BSD 3-Clause License. +# density, from . import ( box, cluster, data, diffraction, + environment, interface, locality, msd, @@ -30,7 +32,7 @@ "data", # "density", "diffraction", - # "environment", + "environment", "interface", "locality", "msd", diff --git a/freud/box.py b/freud/box.py index 4edbb9db3..08c711327 100644 --- a/freud/box.py +++ b/freud/box.py @@ -63,7 +63,9 @@ def __init__(self, Lx, Ly, Lz=0, xy=0, xz=0, yz=0, is2D=None): elif not (Lx and Ly and Lz): msg = "Lx, Ly, and Lz must be nonzero for 3D boxes." raise ValueError(msg) - self._cpp_obj = freud._box.Box(Lx, Ly, Lz, xy, xz, yz, is2D) + self._cpp_obj = freud._box.Box( + *[float(x) for x in [Lx, Ly, Lz, xy, xz, yz]], bool(is2D) + ) @property def L(self): diff --git a/freud/environment.pyx b/freud/environment.py similarity index 68% rename from freud/environment.pyx rename to freud/environment.py index d7a67fb10..a5e2e3071 100644 --- a/freud/environment.pyx +++ b/freud/environment.py @@ -8,33 +8,185 @@ characterize the particle environment. """ +import collections.abc +import warnings + +import numpy as np + +import freud._environment +import freud.box +import freud.locality +import freud.util +from freud._util import ( # noqa F401 + ManagedArray_double, + ManagedArray_float, + ManagedArray_unsignedint, + ManagedArrayVec3_float, + Vector_double, + Vector_float, + Vector_unsignedint, + VectorVec3_float, +) from freud.errors import NO_DEFAULT_QUERY_ARGS_MESSAGE +from freud.locality import _Compute, _PairCompute, _SpatialHistogram -from cython.operator cimport dereference -from libcpp.map cimport map -from freud.locality cimport _PairCompute, _SpatialHistogram -from freud.util cimport _Compute, quat, vec3 +class AngularSeparationNeighbor(_PairCompute): + r"""Calculates the minimum angles of separation between orientations and + query orientations.""" -import warnings + def __init__(self): + self._cpp_obj = freud._environment.AngularSeparationNeighbor() + + def compute( + self, + system, + orientations, + query_points=None, + query_orientations=None, + equiv_orientations=((1, 0, 0, 0),), + neighbors=None, + ): + r"""Calculates the minimum angles of separation between :code:`orientations` + and :code:`query_orientations`, checking for underlying symmetry as encoded + in :code:`equiv_orientations`. The result is stored in the :code:`neighbor_angles` + class attribute. -import numpy as np + Args: + system: + Any object that is a valid argument to + :class:`freud.locality.NeighborQuery.from_system`. + orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): + Orientations associated with system points that are used to + calculate bonds. + query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, + optional): + Query points used to calculate the correlation function. Uses + the system's points if :code:`None` (Default + value = :code:`None`). + query_orientations ((:math:`N_{query\_points}`, 4) :class:`numpy.ndarray`, + optional): + Query orientations used to calculate bonds. Uses + :code:`orientations` if :code:`None`. (Default + value = :code:`None`). + equiv_orientations ((:math:`N_{equiv}`, 4) :class:`numpy.ndarray`, + optional): + The set of all equivalent quaternions that map the particle + to itself (the elements of its rotational symmetry group). + Important: :code:`equiv_orientations` must include both + :math:`q` and :math:`-q`, for all included quaternions. Note + that this calculation assumes that all points in the system + share the same set of equivalent orientations. + (Default value = :code:`((1, 0, 0, 0),)`) + neighbors (:class:`freud.locality.NeighborList` or dict, optional): + Either a :class:`NeighborList ` of + neighbor pairs to use in the calculation, or a dictionary of + `query arguments + `_ + (Default value: None). + """ # noqa: E501 + equiv_orientations = np.asarray(equiv_orientations) + nq, nlist, qargs, query_points, num_query_points = self._preprocess_arguments( + system, query_points, neighbors + ) -import freud.locality + orientations = freud.util._convert_array( + orientations, shape=(nq.points.shape[0], 4) + ) + if query_orientations is None: + query_orientations = orientations + else: + query_orientations = freud.util._convert_array( + query_orientations, shape=(query_points.shape[0], 4) + ) -cimport numpy as np + equiv_orientations = freud.util._convert_array( + equiv_orientations, shape=(None, 4) + ) + + self._cpp_obj.compute( + nq._cpp_obj, + orientations, + query_points, + query_orientations, + equiv_orientations, + nlist._cpp_obj, + qargs._cpp_obj, + ) + return self -cimport freud._environment -cimport freud.box -cimport freud.locality -cimport freud.util + @_Compute._computed_property + def angles(self): + """:math:`\\left(N_{bonds}\\right)` :class:`numpy.ndarray`: The + neighbor angles in radians. The angles are stored in the order of the + neighborlist object.""" + return self._cpp_obj.getAngles().toNumpyArray() + + def __repr__(self): + return f"freud.environment.{type(self).__name__}()" + + @_Compute._computed_property + def nlist(self): + """:class:`freud.locality.NeighborList`: The neighbor list from the + last compute.""" + nlist = freud.locality._nlist_from_cnlist(self._cpp_obj.getNList()) + nlist._compute = self + return nlist -# numpy must be initialized. When using numpy from C or Cython you must -# _always_ do that, or you will have segfaults -np.import_array() +class AngularSeparationGlobal(_Compute): + r"""Calculates the minimum angles of separation between orientations and + global orientations.""" -cdef class BondOrder(_SpatialHistogram): + def __init__(self): + self._cpp_obj = freud._environment.AngularSeparationGlobal() + + def compute( + self, global_orientations, orientations, equiv_orientations=((1, 0, 0, 0),) + ): + r"""Calculates the minimum angles of separation between + :code:`global_orientations` and :code:`orientations`, checking for + underlying symmetry as encoded in :code:`equiv_orientations`. The + result is stored in the :code:`global_angles` class attribute. + + Args: + global_orientations ((:math:`N_{global}`, 4) :class:`numpy.ndarray`): + Set of global orientations to calculate the order + parameter. + orientations ((:math:`N_{particles}`, 4) :class:`numpy.ndarray`): + Orientations to calculate the order parameter. + equiv_orientations ((:math:`N_{equiv}`, 4) :class:`numpy.ndarray`, optional): + The set of all equivalent quaternions that map the particle + to itself (the elements of its rotational symmetry group). + Important: :code:`equiv_orientations` must include both + :math:`q` and :math:`-q`, for all included quaternions. Note + that this calculation assumes that all points in the system + share the same set of equivalent orientations. + (Default value = :code:`((1, 0, 0, 0),)`) + """ # noqa: E501 + equiv_orientations = np.asarray(equiv_orientations) + global_orientations = freud.util._convert_array( + global_orientations, shape=(None, 4) + ) + orientations = freud.util._convert_array(orientations, shape=(None, 4)) + equiv_orientations = freud.util._convert_array( + equiv_orientations, shape=(None, 4) + ) + + self._cpp_obj.compute(global_orientations, orientations, equiv_orientations) + return self + + @_Compute._computed_property + def angles(self): + """:math:`\\left(N_{orientations}, N_{global\\_orientations}\\right)` :class:`numpy.ndarray`: + The global angles in radians.""" # noqa: E501 + return self._cpp_obj.getAngles().toNumpyArray() + + def __repr__(self): + return f"freud.environment.{type(self).__name__}()" + + +class BondOrder(_SpatialHistogram): r"""Compute the bond orientational order diagram for the system of particles. @@ -101,42 +253,46 @@ Mode to calculate bond order. Options are :code:`'bod'`, :code:`'lbod'`, :code:`'obcd'`, or :code:`'oocd'` (Default value = :code:`'bod'`). - """ # noqa: E501 - cdef freud._environment.BondOrder * thisptr + """ - known_modes = {'bod': freud._environment.bod, - 'lbod': freud._environment.lbod, - 'obcd': freud._environment.obcd, - 'oocd': freud._environment.oocd} + known_modes = { + "bod": freud._environment.bod, + "lbod": freud._environment.lbod, + "obcd": freud._environment.obcd, + "oocd": freud._environment.oocd, + } - def __cinit__(self, bins, str mode="bod"): - try: + def __init__(self, bins, mode="bod"): + if isinstance(bins, collections.abc.Sequence): n_bins_theta, n_bins_phi = bins - except TypeError: + else: n_bins_theta = n_bins_phi = bins - cdef freud._environment.BondOrderMode l_mode try: l_mode = self.known_modes[mode] - except KeyError: - raise ValueError( - 'Unknown BondOrder mode: {}'.format(mode)) - - self.thisptr = self.histptr = new freud._environment.BondOrder( - n_bins_theta, n_bins_phi, l_mode) + except KeyError as err: + msg = f"Unknown BondOrder mode: {mode}" + raise ValueError(msg) from err - def __dealloc__(self): - del self.thisptr + self._cpp_obj = freud._environment.BondOrder(n_bins_theta, n_bins_phi, l_mode) @property def default_query_args(self): """No default query arguments.""" # Must override the generic histogram's defaults. raise NotImplementedError( - NO_DEFAULT_QUERY_ARGS_MESSAGE.format(type(self).__name__)) - - def compute(self, system, orientations=None, query_points=None, - query_orientations=None, neighbors=None, reset=True): + NO_DEFAULT_QUERY_ARGS_MESSAGE.format(type(self).__name__) + ) + + def compute( + self, + system, + orientations=None, + query_points=None, + query_orientations=None, + neighbors=None, + reset=True, + ): r"""Calculates the correlation function and adds to the current histogram. @@ -148,11 +304,13 @@ def compute(self, system, orientations=None, query_points=None, Orientations associated with system points that are used to calculate bonds. Uses identity quaternions if :code:`None` (Default value = :code:`None`). - query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, optional): + query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, + optional): Query points used to calculate the correlation function. Uses the system's points if :code:`None` (Default value = :code:`None`). - query_orientations ((:math:`N_{query\_points}`, 4) :class:`numpy.ndarray`, optional): + query_orientations ((:math:`N_{query\_points}`, 4) :class:`numpy.ndarray`, + optional): Query orientations used to calculate bonds. Uses :code:`orientations` if :code:`None`. (Default value = :code:`None`). @@ -166,69 +324,64 @@ def compute(self, system, orientations=None, query_points=None, Whether to erase the previously computed values before adding the new computation; if False, will accumulate data (Default value: True). - """ # noqa: E501 + """ if reset: self._reset() - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, query_points, neighbors) + nq, nlist, qargs, l_query_points, num_query_points = self._preprocess_arguments( + system, query_points, neighbors + ) if orientations is None: orientations = np.array([[1, 0, 0, 0]] * nq.points.shape[0]) if query_orientations is None: query_orientations = orientations + if query_points is None: + query_points = nq.points orientations = freud.util._convert_array( - orientations, shape=(nq.points.shape[0], 4)) + orientations, shape=(nq.points.shape[0], 4) + ) query_orientations = freud.util._convert_array( - query_orientations, shape=(num_query_points, 4)) - - cdef const float[:, ::1] l_orientations = orientations - cdef const float[:, ::1] l_query_orientations = query_orientations - - self.thisptr.accumulate( - nq.get_ptr(), - &l_orientations[0, 0], - &l_query_points[0, 0], - &l_query_orientations[0, 0], - num_query_points, - nlist.get_ptr(), dereference(qargs.thisptr)) + query_orientations, shape=(num_query_points, 4) + ) + query_points = freud.util._convert_array( + query_points, shape=(num_query_points, 3) + ) + + self._cpp_obj.accumulate( + nq._cpp_obj, + orientations, + query_points, + query_orientations, + nlist._cpp_obj, + qargs._cpp_obj, + ) return self @_Compute._computed_property def bond_order(self): """:math:`\\left(N_{\\phi}, N_{\\theta} \\right)` :class:`numpy.ndarray`: Bond order.""" # noqa: E501 - return freud.util.make_managed_numpy_array( - &self.thisptr.getBondOrder(), - freud.util.arr_type_t.FLOAT) + return self._cpp_obj.getBondOrder().toNumpyArray() @_Compute._computed_property def box(self): """:class:`freud.box.Box`: Box used in the calculation.""" - return freud.box.BoxFromCPP(self.thisptr.getBox()) + return freud.box.BoxFromCPP(self._cpp_obj.getBox()) def __repr__(self): - return ("freud.environment.{cls}(bins=({bins}), mode='{mode}')".format( + return "freud.environment.{cls}(bins=({bins}), mode='{mode}')".format( cls=type(self).__name__, - bins=', '.join([str(b) for b in self.nbins]), - mode=self.mode)) + bins=", ".join([str(b) for b in self.nbins]), + mode=self.mode, + ) @property def mode(self): """str: Bond order mode.""" - mode = self.thisptr.getMode() - for key, value in self.known_modes.items(): - if value == mode: - return key + return self._cpp_obj.getMode().name -cdef class LocalDescriptors(_PairCompute): +class LocalDescriptors(_PairCompute): r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's local environment. @@ -256,29 +409,31 @@ def mode(self): neighborhood, :code:`'particle_local'` to use the given particle orientations, or :code:`'global'` to not rotate environments (Default value = :code:`'neighborhood'`). - """ # noqa: E501 - cdef freud._environment.LocalDescriptors * thisptr + """ - known_modes = {'neighborhood': freud._environment.LocalNeighborhood, - 'global': freud._environment.Global, - 'particle_local': freud._environment.ParticleLocal} + known_modes = { + "neighborhood": freud._environment.LocalNeighborhood, + "global": freud._environment.Global, + "particle_local": freud._environment.ParticleLocal, + } - def __cinit__(self, l_max, negative_m=True, mode='neighborhood'): - cdef freud._environment.LocalDescriptorOrientation l_mode + def __init__(self, l_max, negative_m=True, mode="neighborhood"): try: l_mode = self.known_modes[mode] except KeyError: - raise ValueError( - 'Unknown LocalDescriptors orientation mode: {}'.format(mode)) - - self.thisptr = new freud._environment.LocalDescriptors( - l_max, negative_m, l_mode) - - def __dealloc__(self): - del self.thisptr - - def compute(self, system, query_points=None, orientations=None, - neighbors=None, max_num_neighbors=0): + msg = f"Unknown LocalDescriptors orientation mode: {mode}" + raise ValueError(msg) from None + + self._cpp_obj = freud._environment.LocalDescriptors(l_max, negative_m, l_mode) + + def compute( + self, + system, + query_points=None, + orientations=None, + neighbors=None, + max_num_neighbors=0, + ): r"""Calculates the local descriptors of bonds from a set of source points to a set of destination points. @@ -304,44 +459,39 @@ def compute(self, system, query_points=None, orientations=None, particle for the given neighbor-finding algorithm. Uses all neighbors if set to 0 (Default value = 0). """ # noqa: E501 - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, query_points, neighbors) + nq, nlist, qargs, query_points, num_query_points = self._preprocess_arguments( + system, query_points, neighbors + ) # The l_orientations_ptr is only used for 'particle_local' mode. - cdef const float[:, ::1] l_orientations - cdef quat[float] *l_orientations_ptr = NULL - if self.mode == 'particle_local': + if self.mode == "particle_local": if orientations is None: - raise RuntimeError( - ('Orientations must be given to orient LocalDescriptors ' - 'with particles\' orientations')) + msg = ( + "Orientations must be given to orient LocalDescriptors " + "with particles' orientations" + ) + raise RuntimeError(msg) orientations = freud.util._convert_array( - orientations, shape=(nq.points.shape[0], 4)) - - l_orientations = orientations - l_orientations_ptr = &l_orientations[0, 0] + orientations, shape=(nq.points.shape[0], 4) + ) - self.thisptr.compute( - nq.get_ptr(), - &l_query_points[0, 0], num_query_points, - l_orientations_ptr, - nlist.get_ptr(), dereference(qargs.thisptr), - max_num_neighbors) + self._cpp_obj.compute( + nq._cpp_obj, + query_points, + num_query_points, + orientations, + nlist._cpp_obj, + qargs._cpp_obj, + max_num_neighbors, + ) return self @_Compute._computed_property def nlist(self): """:class:`freud.locality.NeighborList`: The neighbor list from the last compute.""" - nlist = freud.locality._nlist_from_cnlist(self.thisptr.getNList()) + nlist = freud.locality._nlist_from_cnlist(self._cpp_obj.getNList()) nlist._compute = self return nlist @@ -349,9 +499,7 @@ def nlist(self): def sph(self): """:math:`\\left(N_{bonds}, \\text{SphWidth} \\right)` :class:`numpy.ndarray`: The last computed spherical harmonic array.""" - return freud.util.make_managed_numpy_array( - &self.thisptr.getSph(), - freud.util.arr_type_t.COMPLEX_FLOAT) + return self._cpp_obj.getSph().toNumpyArray() @_Compute._computed_property def num_sphs(self): @@ -361,19 +509,19 @@ def num_sphs(self): :code:`num_neighbors` arguments passed to the last compute call or the constructor (it may be less if there are not enough neighbors for every particle).""" - return self.thisptr.getNSphs() + return self._cpp_obj.getNSphs() @property def l_max(self): """unsigned int: The maximum spherical harmonic :math:`l` calculated for.""" - return self.thisptr.getLMax() + return self._cpp_obj.getLMax() @property def negative_m(self): """bool: True if we also calculated :math:`Y_{lm}` for negative :math:`m`.""" - return self.thisptr.getNegativeM() + return self._cpp_obj.getNegativeM() @property def mode(self): @@ -381,16 +529,17 @@ def mode(self): :code:`'neighborhood'` to use the orientation of the local neighborhood, :code:`'particle_local'` to use the given particle orientations, or :code:`'global'` to not rotate environments.""" - mode = self.thisptr.getMode() + mode = self._cpp_obj.getMode() for key, value in self.known_modes.items(): if value == mode: return key + return None def __repr__(self): - return ("freud.environment.{cls}(l_max={l_max}, " - "negative_m={negative_m}, mode='{mode}')").format( - cls=type(self).__name__, l_max=self.l_max, - negative_m=self.negative_m, mode=self.mode) + return ( + f"freud.environment.{type(self).__name__}(l_max={self.l_max}, " + f"negative_m={self.negative_m}, mode='{self.mode}')" + ) def _minimize_RMSD(box, ref_points, points, registration=False): @@ -414,29 +563,26 @@ def _minimize_RMSD(box, ref_points, points, registration=False): set of points, and the mapping between the vectors of ref_points and points that somewhat minimizes the RMSD. """ # noqa: E501 - cdef freud.box.Box b = freud.util._convert_box(box) + b = freud.util._convert_box(box) ref_points = freud.util._convert_array(ref_points, shape=(None, 3)) points = freud.util._convert_array(points, shape=(None, 3)) - cdef const float[:, ::1] l_ref_points = ref_points - cdef const float[:, ::1] l_points = points - cdef unsigned int nRef1 = l_ref_points.shape[0] - cdef unsigned int nRef2 = l_points.shape[0] + nRef1 = ref_points.shape[0] + nRef2 = points.shape[0] if nRef1 != nRef2: - raise ValueError( - ("The number of vectors in ref_points must MATCH" - "the number of vectors in points")) + msg = ( + "The number of vectors in ref_points must MATCH" + "the number of vectors in points" + ) + raise ValueError(msg) - cdef float min_rmsd = -1 - cdef map[unsigned int, unsigned int] results_map = \ - freud._environment.minimizeRMSD( - dereference(b.thisptr), - &l_ref_points[0, 0], - &l_points[0, 0], - nRef1, min_rmsd, registration) - return [min_rmsd, np.asarray(l_points), results_map] + min_rmsd = -1 + min_rmsd, results_map = freud._environment.minimizeRMSD( + b._cpp_obj, ref_points, points, nRef1, min_rmsd, registration + ) + return [min_rmsd, np.asarray(points), results_map] def _is_similar_motif(box, ref_points, points, threshold, registration=False): @@ -467,70 +613,62 @@ def _is_similar_motif(box, ref_points, points, threshold, registration=False): correspond to each other. Empty if they do not correspond to each other. """ # noqa: E501 - cdef freud.box.Box b = freud.util._convert_box(box) + b = freud.util._convert_box(box) ref_points = freud.util._convert_array(ref_points, shape=(None, 3)) points = freud.util._convert_array(points, shape=(None, 3)) - cdef const float[:, ::1] l_ref_points = ref_points - cdef const float[:, ::1] l_points = points - cdef unsigned int nRef1 = l_ref_points.shape[0] - cdef unsigned int nRef2 = l_points.shape[0] - cdef float threshold_sq = threshold*threshold + nRef1 = ref_points.shape[0] + nRef2 = points.shape[0] + threshold_sq = threshold * threshold if nRef1 != nRef2: - raise ValueError( - ("The number of vectors in ref_points must match" - "the number of vectors in points")) + msg = ( + "The number of vectors in ref_points must match" + "the number of vectors in points" + ) + raise ValueError(msg) - cdef map[unsigned int, unsigned int] vec_map = \ - freud._environment.isSimilar( - dereference(b.thisptr), &l_ref_points[0, 0], - &l_points[0, 0], nRef1, threshold_sq, - registration) - return [np.asarray(l_points), vec_map] + vec_map = freud._environment.isSimilar( + b._cpp_obj, ref_points, points, nRef1, threshold_sq, registration + ) + return [np.asarray(points), vec_map] -cdef class _MatchEnv(_PairCompute): - r"""Parent for environment matching methods. """ - cdef freud._environment.MatchEnv * matchptr +class _MatchEnv(_PairCompute): + r"""Parent for environment matching methods.""" - def __cinit__(self, *args, **kwargs): + def __init__(self, *args, **kwargs): # Abstract class - pass + self._cpp_obj = freud._environment.MatchEnv() @_Compute._computed_property def point_environments(self): """:math:`\\left(N_{points}, N_{neighbors}, 3\\right)` 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] + envs = self._cpp_obj.getPointEnvironments() + return [np.asarray(env) for env in envs] def __repr__(self): - return ("freud.environment.{cls}()").format( - cls=type(self).__name__) + return f"freud.environment.{type(self).__name__}()" -cdef class EnvironmentCluster(_MatchEnv): +class EnvironmentCluster(_MatchEnv): r"""Clusters particles according to whether their local environments match or not, using various shape matching metrics defined in :cite:`Teich2019`. """ - cdef freud._environment.EnvironmentCluster * thisptr - - def __cinit__(self): - self.thisptr = self.matchptr = \ - new freud._environment.EnvironmentCluster() - def __init__(self): - pass - - def __dealloc__(self): - del self.thisptr - - def compute(self, system, threshold, cluster_neighbors=None, - env_neighbors=None, registration=False): + self._cpp_obj = freud._environment.EnvironmentCluster() + + def compute( + self, + system, + threshold, + cluster_neighbors=None, + env_neighbors=None, + registration=False, + ): r"""Determine clusters of particles with matching environments. An environment is defined by the bond vectors between a particle and its @@ -638,46 +776,42 @@ def compute(self, system, threshold, cluster_neighbors=None, option incurs a significant performance penalty. (Default value = :code:`False`) """ # noqa: E501 - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist, env_nlist - freud.locality._QueryArgs qargs, env_qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, neighbors=cluster_neighbors) + nq, nlist, qargs, l_query_points, num_query_points = self._preprocess_arguments( + system, neighbors=cluster_neighbors + ) if env_neighbors is None: env_neighbors = cluster_neighbors env_nlist, env_qargs = self._resolve_neighbors(env_neighbors) - self.thisptr.compute( - nq.get_ptr(), nlist.get_ptr(), dereference(qargs.thisptr), - env_nlist.get_ptr(), dereference(env_qargs.thisptr), threshold, - registration) + self._cpp_obj.compute( + nq._cpp_obj, + nlist._cpp_obj, + qargs._cpp_obj, + env_nlist._cpp_obj, + env_qargs._cpp_obj, + threshold, + registration, + ) return self @_Compute._computed_property def cluster_idx(self): """:math:`\\left(N_{particles}\\right)` :class:`numpy.ndarray`: The per-particle index indicating cluster membership.""" - return freud.util.make_managed_numpy_array( - &self.thisptr.getClusters(), - freud.util.arr_type_t.UNSIGNED_INT) + return self._cpp_obj.getClusters().toNumpyArray() @_Compute._computed_property def num_clusters(self): """unsigned int: The number of clusters.""" - return self.thisptr.getNumClusters() + return self._cpp_obj.getNumClusters() @_Compute._computed_property def cluster_environments(self): """:math:`\\left(N_{clusters}, N_{neighbors}, 3\\right)` 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] + envs = self._cpp_obj.getClusterEnvironments() + return [np.asarray(env) for env in envs] def plot(self, ax=None): """Plot cluster distribution. @@ -691,23 +825,26 @@ def plot(self, ax=None): (:class:`matplotlib.axes.Axes`): Axis with the plot. """ import freud.plot + try: values, counts = np.unique(self.cluster_idx, return_counts=True) except ValueError: return None else: return freud.plot.clusters_plot( - values, counts, num_clusters_to_plot=10, ax=ax) + values, counts, num_clusters_to_plot=10, ax=ax + ) def _repr_png_(self): try: import freud.plot + return freud.plot._ax_to_bytes(self.plot()) except (AttributeError, ImportError): return None -cdef class EnvironmentMotifMatch(_MatchEnv): +class EnvironmentMotifMatch(_MatchEnv): r"""Find matches between local arrangements of a set of points and a provided motif, as done in :cite:`Teich2019`. @@ -716,19 +853,12 @@ def _repr_png_(self): different number of neighbors than the motif will always fail to match the motif. See :class:`freud.environment.EnvironmentCluster.compute` for the matching criterion. - """ # noqa: E501 - - cdef freud._environment.EnvironmentMotifMatch * thisptr - - def __cinit__(self): - self.thisptr = self.matchptr = \ - new freud._environment.EnvironmentMotifMatch() + """ def __init__(self): - pass + self._cpp_obj = freud._environment.EnvironmentMotifMatch() - def compute(self, system, motif, threshold, env_neighbors=None, - registration=False): + def compute(self, system, motif, threshold, env_neighbors=None, registration=False): r"""Determine which particles have local environments matching the given environment motif. @@ -763,44 +893,40 @@ def compute(self, system, motif, threshold, env_neighbors=None, it minimizes the RMSD between the two sets (Default value = False). """ - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, neighbors=env_neighbors) + nq, nlist, qargs, l_query_points, num_query_points = self._preprocess_arguments( + system, neighbors=env_neighbors + ) motif = freud.util._convert_array(motif, shape=(None, 3)) if (motif == 0).all(axis=1).any(): warnings.warn( "Attempting to match a motif containing the zero " "vector is likely to result in zero matches.", - RuntimeWarning + RuntimeWarning, + stacklevel=2, ) - - cdef const float[:, ::1] l_motif = motif - cdef unsigned int nRef = l_motif.shape[0] - - self.thisptr.compute( - nq.get_ptr(), nlist.get_ptr(), dereference(qargs.thisptr), - - &l_motif[0, 0], nRef, - threshold, registration) + nRef = motif.shape[0] + + self._cpp_obj.compute( + nq._cpp_obj, + nlist._cpp_obj, + qargs._cpp_obj, + motif, + nRef, + threshold, + registration, + ) return self @_Compute._computed_property def matches(self): """(:math:`N_{points}`) :class:`numpy.ndarray`: A boolean array indicating whether each point matches the motif.""" - return freud.util.make_managed_numpy_array( - &self.thisptr.getMatches(), - freud.util.arr_type_t.BOOL) + # NOTE: Numpy stores bools as a byte each, so this cast should be free + return self._cpp_obj.getMatches().toNumpyArray().astype(bool) -cdef class _EnvironmentRMSDMinimizer(_MatchEnv): +class _EnvironmentRMSDMinimizer(_MatchEnv): r"""Find linear transformations that map the environments of points onto a motif. @@ -814,17 +940,10 @@ def matches(self): enforced (but we could add a warning to the compute...). """ - cdef freud._environment.EnvironmentRMSDMinimizer * thisptr - - def __cinit__(self): - self.thisptr = self.matchptr = \ - new freud._environment.EnvironmentRMSDMinimizer() - def __init__(self): - pass + self._cpp_obj = freud._environment.EnvironmentRMSDMinimizer() - def compute(self, system, motif, neighbors=None, - registration=False): + def compute(self, system, motif, neighbors=None, registration=False): r"""Rotate (if registration=True) and permute the environments of all particles to minimize their RMSD with respect to the motif provided by motif. @@ -851,243 +970,43 @@ def compute(self, system, motif, neighbors=None, Vector of minimal RMSD values, one value per particle. """ - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, neighbors=neighbors) + nq, nlist, qargs, l_query_points, num_query_points = self._preprocess_arguments( + system, neighbors=neighbors + ) motif = freud.util._convert_array(motif, shape=(None, 3)) - cdef const float[:, ::1] l_motif = motif - cdef unsigned int nRef = l_motif.shape[0] - - self.thisptr.compute( - nq.get_ptr(), nlist.get_ptr(), dereference(qargs.thisptr), - - &l_motif[0, 0], nRef, - registration) + nRef = motif.shape[0] + self._cpp_obj.compute( + nq._cpp_obj, nlist._cpp_obj, qargs._cpp_obj, motif, nRef, registration + ) return self @_Compute._computed_property def rmsds(self): """:math:`(N_p, )` :class:`numpy.ndarray`: A boolean array of the RMSDs found for each point's environment.""" - return freud.util.make_managed_numpy_array( - &self.thisptr.getRMSDs(), - freud.util.arr_type_t.FLOAT) - - -cdef class AngularSeparationNeighbor(_PairCompute): - r"""Calculates the minimum angles of separation between orientations and - query orientations.""" - cdef freud._environment.AngularSeparationNeighbor * thisptr - - def __cinit__(self): - self.thisptr = new freud._environment.AngularSeparationNeighbor() - - def __init__(self): - pass - - def __dealloc__(self): - del self.thisptr - - def compute(self, system, orientations, query_points=None, - query_orientations=None, - equiv_orientations=np.array([[1, 0, 0, 0]]), - neighbors=None): - r"""Calculates the minimum angles of separation between :code:`orientations` - and :code:`query_orientations`, checking for underlying symmetry as encoded - in :code:`equiv_orientations`. The result is stored in the :code:`neighbor_angles` - class attribute. - - Args: - system: - Any object that is a valid argument to - :class:`freud.locality.NeighborQuery.from_system`. - orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): - Orientations associated with system points that are used to - calculate bonds. - query_points ((:math:`N_{query\_points}`, 3) :class:`numpy.ndarray`, optional): - Query points used to calculate the correlation function. Uses - the system's points if :code:`None` (Default - value = :code:`None`). - query_orientations ((:math:`N_{query\_points}`, 4) :class:`numpy.ndarray`, optional): - Query orientations used to calculate bonds. Uses - :code:`orientations` if :code:`None`. (Default - value = :code:`None`). - equiv_orientations ((:math:`N_{equiv}`, 4) :class:`numpy.ndarray`, optional): - The set of all equivalent quaternions that map the particle - to itself (the elements of its rotational symmetry group). - Important: :code:`equiv_orientations` must include both - :math:`q` and :math:`-q`, for all included quaternions. Note - that this calculation assumes that all points in the system - share the same set of equivalent orientations. - (Default value = :code:`[[1, 0, 0, 0]]`) - neighbors (:class:`freud.locality.NeighborList` or dict, optional): - Either a :class:`NeighborList ` of - neighbor pairs to use in the calculation, or a dictionary of - `query arguments - `_ - (Default value: None). - """ # noqa: E501 - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points - - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, query_points, neighbors) - - orientations = freud.util._convert_array( - orientations, shape=(nq.points.shape[0], 4)) - if query_orientations is None: - query_orientations = orientations - else: - query_orientations = freud.util._convert_array( - query_orientations, shape=(query_points.shape[0], 4)) - - equiv_orientations = freud.util._convert_array( - equiv_orientations, shape=(None, 4)) - - cdef const float[:, ::1] l_orientations = orientations - cdef const float[:, ::1] l_query_orientations = query_orientations - cdef const float[:, ::1] l_equiv_orientations = equiv_orientations + return self._cpp_obj.getRMSDs().toNumpyArray() - cdef unsigned int n_equiv_orientations = l_equiv_orientations.shape[0] - self.thisptr.compute( - nq.get_ptr(), - &l_orientations[0, 0], - &l_query_points[0, 0], - &l_query_orientations[0, 0], - num_query_points, - &l_equiv_orientations[0, 0], - n_equiv_orientations, - nlist.get_ptr(), - dereference(qargs.thisptr)) - return self - - @_Compute._computed_property - def angles(self): - """:math:`\\left(N_{bonds}\\right)` :class:`numpy.ndarray`: The - neighbor angles in radians. The angles are stored in the order of the - neighborlist object.""" - return freud.util.make_managed_numpy_array( - &self.thisptr.getAngles(), - freud.util.arr_type_t.FLOAT) - - def __repr__(self): - return "freud.environment.{cls}()".format( - cls=type(self).__name__) - - @_Compute._computed_property - def nlist(self): - """:class:`freud.locality.NeighborList`: The neighbor list from the - last compute.""" - nlist = freud.locality._nlist_from_cnlist(self.thisptr.getNList()) - nlist._compute = self - return nlist - - -cdef class AngularSeparationGlobal(_Compute): - r"""Calculates the minimum angles of separation between orientations and - global orientations.""" - cdef freud._environment.AngularSeparationGlobal * thisptr - - def __cinit__(self): - self.thisptr = new freud._environment.AngularSeparationGlobal() - - def __init__(self): - pass - - def __dealloc__(self): - del self.thisptr - - def compute(self, global_orientations, orientations, - equiv_orientations=np.array([[1, 0, 0, 0]])): - r"""Calculates the minimum angles of separation between - :code:`global_orientations` and :code:`orientations`, checking for - underlying symmetry as encoded in :code:`equiv_orientations`. The - result is stored in the :code:`global_angles` class attribute. - - Args: - global_orientations ((:math:`N_{global}`, 4) :class:`numpy.ndarray`): - Set of global orientations to calculate the order - parameter. - orientations ((:math:`N_{particles}`, 4) :class:`numpy.ndarray`): - Orientations to calculate the order parameter. - equiv_orientations ((:math:`N_{equiv}`, 4) :class:`numpy.ndarray`, optional): - The set of all equivalent quaternions that map the particle - to itself (the elements of its rotational symmetry group). - Important: :code:`equiv_orientations` must include both - :math:`q` and :math:`-q`, for all included quaternions. Note - that this calculation assumes that all points in the system - share the same set of equivalent orientations. - (Default value = :code:`[[1, 0, 0, 0]]`) - """ # noqa: E501 - global_orientations = freud.util._convert_array( - global_orientations, shape=(None, 4)) - orientations = freud.util._convert_array( - orientations, shape=(None, 4)) - equiv_orientations = freud.util._convert_array( - equiv_orientations, shape=(None, 4)) - - cdef const float[:, ::1] l_global_orientations = global_orientations - cdef const float[:, ::1] l_orientations = orientations - cdef const float[:, ::1] l_equiv_orientations = equiv_orientations - - cdef unsigned int n_global = l_global_orientations.shape[0] - cdef unsigned int n_points = l_orientations.shape[0] - cdef unsigned int n_equiv_orientations = l_equiv_orientations.shape[0] - - self.thisptr.compute( - &l_global_orientations[0, 0], - n_global, - &l_orientations[0, 0], - n_points, - &l_equiv_orientations[0, 0], - n_equiv_orientations) - return self - - @_Compute._computed_property - def angles(self): - """:math:`\\left(N_{orientations}, N_{global\\_orientations}\\right)` :class:`numpy.ndarray`: - The global angles in radians.""" # noqa: E501 - return freud.util.make_managed_numpy_array( - &self.thisptr.getAngles(), - freud.util.arr_type_t.FLOAT) - - def __repr__(self): - return "freud.environment.{cls}()".format( - cls=type(self).__name__) - - -cdef class LocalBondProjection(_PairCompute): +class LocalBondProjection(_PairCompute): r"""Calculates the maximal projection of nearest neighbor bonds for each particle onto some set of reference vectors, defined in the particles' local reference frame. """ - cdef freud._environment.LocalBondProjection * thisptr - - def __cinit__(self): - self.thisptr = new freud._environment.LocalBondProjection() def __init__(self): - pass - - def __dealloc__(self): - del self.thisptr - - def compute(self, system, orientations, proj_vecs, - query_points=None, equiv_orientations=np.array([[1, 0, 0, 0]]), - neighbors=None): + self._cpp_obj = freud._environment.LocalBondProjection() + + def compute( + self, + system, + orientations, + proj_vecs, + query_points=None, + equiv_orientations=((1, 0, 0, 0),), + neighbors=None, + ): r"""Calculates the maximal projections of nearest neighbor bonds (between :code:`points` and :code:`query_points`) onto the set of reference vectors :code:`proj_vecs`, defined in the local reference @@ -1118,7 +1037,7 @@ def compute(self, system, orientations, proj_vecs, :math:`q` and :math:`-q`, for all included quaternions. Note that this calculation assumes that all points in the system share the same set of equivalent orientations. - (Default value = :code:`[[1, 0, 0, 0]]`) + (Default value = :code:`((1, 0, 0, 0),)`) neighbors (:class:`freud.locality.NeighborList` or dict, optional): Either a :class:`NeighborList ` of neighbor pairs to use in the calculation, or a dictionary of @@ -1126,44 +1045,34 @@ def compute(self, system, orientations, proj_vecs, `_ (Default value: None). """ # noqa: E501 - cdef: - freud.locality.NeighborQuery nq - freud.locality.NeighborList nlist - freud.locality._QueryArgs qargs - const float[:, ::1] l_query_points - unsigned int num_query_points + equiv_orientations = np.asarray(equiv_orientations) + nq, nlist, qargs, query_points, num_query_points = self._preprocess_arguments( + system, query_points, neighbors + ) - nq, nlist, qargs, l_query_points, num_query_points = \ - self._preprocess_arguments(system, query_points, neighbors) - - orientations = freud.util._convert_array( - orientations, shape=(None, 4)) + orientations = freud.util._convert_array(orientations, shape=(None, 4)) equiv_orientations = freud.util._convert_array( - equiv_orientations, shape=(None, 4)) + equiv_orientations, shape=(None, 4) + ) proj_vecs = freud.util._convert_array(proj_vecs, shape=(None, 3)) - cdef const float[:, ::1] l_orientations = orientations - cdef const float[:, ::1] l_equiv_orientations = equiv_orientations - cdef const float[:, ::1] l_proj_vecs = proj_vecs - - cdef unsigned int n_equiv = l_equiv_orientations.shape[0] - cdef unsigned int n_proj = l_proj_vecs.shape[0] - - self.thisptr.compute( - nq.get_ptr(), - &l_orientations[0, 0], - &l_query_points[0, 0], num_query_points, - &l_proj_vecs[0, 0], n_proj, - &l_equiv_orientations[0, 0], n_equiv, - nlist.get_ptr(), dereference(qargs.thisptr)) + self._cpp_obj.compute( + nq._cpp_obj, + orientations, + query_points, + proj_vecs, + equiv_orientations, + nlist._cpp_obj, + qargs._cpp_obj, + ) return self @_Compute._computed_property def nlist(self): """:class:`freud.locality.NeighborList`: The neighbor list from the last compute.""" - nlist = freud.locality._nlist_from_cnlist(self.thisptr.getNList()) + nlist = freud.locality._nlist_from_cnlist(self._cpp_obj.getNList()) nlist._compute = self return nlist @@ -1172,9 +1081,7 @@ def projections(self): """:math:`\\left(N_{bonds}, N_{projection\\_vecs} \\right)` :class:`numpy.ndarray`: The projection of each bond between query particles and their neighbors onto each of the projection vectors.""" # noqa: E501 - return freud.util.make_managed_numpy_array( - &self.thisptr.getProjections(), - freud.util.arr_type_t.FLOAT) + return self._cpp_obj.getProjections().toNumpyArray() @_Compute._computed_property def normed_projections(self): @@ -1182,9 +1089,7 @@ def normed_projections(self): The projection of each bond between query particles and their neighbors onto each of the projection vectors, normalized by the length of the bond.""" # noqa: E501 - return freud.util.make_managed_numpy_array( - &self.thisptr.getNormedProjections(), - freud.util.arr_type_t.FLOAT) + return self._cpp_obj.getNormedProjections().toNumpyArray() def __repr__(self): - return ("freud.environment.{cls}()").format(cls=type(self).__name__) + return f"freud.environment.{type(self).__name__}()" diff --git a/freud/environment/AngularSeparation.cc b/freud/environment/AngularSeparation.cc index 6646491be..2693b7c6a 100644 --- a/freud/environment/AngularSeparation.cc +++ b/freud/environment/AngularSeparation.cc @@ -1,8 +1,17 @@ // Copyright (c) 2010-2024 The Regents of the University of Michigan // This file is from the freud project, released under the BSD 3-Clause License. +#include +#include +#include +#include + #include "AngularSeparation.h" +#include "ManagedArray.h" #include "NeighborComputeFunctional.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" #include "utils.h" /*! \file AngularSeparation.cc @@ -13,7 +22,7 @@ namespace freud { namespace environment { float computeSeparationAngle(const quat& ref_q, const quat& q) { - quat R = q * conj(ref_q); + const quat R = q * conj(ref_q); auto theta = float(2.0 * std::acos(util::clamp(R.s, -1, 1))); return theta; } @@ -25,9 +34,9 @@ float computeSeparationAngle(const quat& ref_q, const quat& q) float computeMinSeparationAngle(const quat& ref_q, const quat& q, const quat* equiv_qs, unsigned int n_equiv_quats) { - quat qconst = equiv_qs[0]; + const quat qconst = equiv_qs[0]; // here we undo a rotation represented by one of the equivalent orientations - quat qtemp = q * conj(qconst); + const quat qtemp = q * conj(qconst); // start with the quaternion before it has been rotated by equivalent rotations float min_angle = computeSeparationAngle(ref_q, q); @@ -35,10 +44,10 @@ float computeMinSeparationAngle(const quat& ref_q, const quat& q, // loop through all equivalent rotations and see if they have smaller angles with ref_q for (unsigned int i = 0; i < n_equiv_quats; ++i) { - quat qe = equiv_qs[i]; - quat qtest = qtemp * qe; + const quat qe = equiv_qs[i]; + const quat qtest = qtemp * qe; - float angle_test = computeSeparationAngle(ref_q, qtest); + const float angle_test = computeSeparationAngle(ref_q, qtest); if (angle_test < min_angle) { @@ -49,32 +58,35 @@ float computeMinSeparationAngle(const quat& ref_q, const quat& q, return min_angle; } -void AngularSeparationNeighbor::compute(const locality::NeighborQuery* nq, const quat* orientations, - const vec3* query_points, +void AngularSeparationNeighbor::compute(const std::shared_ptr& nq, + const quat* orientations, const vec3* query_points, const quat* query_orientations, unsigned int n_query_points, const quat* equiv_orientations, unsigned int n_equiv_orientations, - const freud::locality::NeighborList* nlist, locality::QueryArgs qargs) + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs) { // This function requires a NeighborList object, so we always make one and store it locally. + m_nlist = locality::makeDefaultNlist(nq, nlist, query_points, n_query_points, qargs); - const size_t tot_num_neigh = m_nlist.getNumBonds(); - m_angles.prepare(tot_num_neigh); + const size_t tot_num_neigh = m_nlist->getNumBonds(); + m_angles = std::make_shared>(std::vector {tot_num_neigh}); util::forLoopWrapper(0, nq->getNPoints(), [&](size_t begin, size_t end) { - size_t bond(m_nlist.find_first_index(begin)); + size_t bond(m_nlist->find_first_index(begin)); for (size_t i = begin; i < end; ++i) { - quat q = orientations[i]; + const quat q = orientations[i]; - for (; bond < tot_num_neigh && m_nlist.getNeighbors()(bond, 0) == i; ++bond) + for (; bond < tot_num_neigh && (*m_nlist->getNeighbors())(bond, 0) == i; ++bond) { - const size_t j(m_nlist.getNeighbors()(bond, 1)); - quat query_q = query_orientations[j]; + const size_t j((*m_nlist->getNeighbors())(bond, 1)); + const quat query_q = query_orientations[j]; - float theta = computeMinSeparationAngle(q, query_q, equiv_orientations, n_equiv_orientations); - m_angles[bond] = theta; + const float theta + = computeMinSeparationAngle(q, query_q, equiv_orientations, n_equiv_orientations); + (*m_angles)[bond] = theta; } } }); @@ -85,18 +97,18 @@ void AngularSeparationGlobal::compute(const quat* global_orientations, un const quat* equiv_orientations, unsigned int n_equiv_orientations) { - m_angles.prepare({n_points, n_global}); + m_angles = std::make_shared>(std::vector {n_points, n_global}); util::forLoopWrapper(0, n_points, [&](size_t begin, size_t end) { for (size_t i = begin; i < end; ++i) { - quat q = orientations[i]; + const quat q = orientations[i]; for (unsigned int j = 0; j < n_global; j++) { - quat global_q = global_orientations[j]; - float theta + const quat global_q = global_orientations[j]; + const float theta = computeMinSeparationAngle(q, global_q, equiv_orientations, n_equiv_orientations); - m_angles(i, j) = theta; + (*m_angles)(i, j) = theta; } } }); diff --git a/freud/environment/AngularSeparation.h b/freud/environment/AngularSeparation.h index a45168959..0a20e144c 100644 --- a/freud/environment/AngularSeparation.h +++ b/freud/environment/AngularSeparation.h @@ -4,6 +4,9 @@ #ifndef ANGULAR_SEPARATION_H #define ANGULAR_SEPARATION_H +#include + +#include "ManagedArray.h" #include "NeighborList.h" #include "NeighborQuery.h" #include "VectorMath.h" @@ -36,13 +39,13 @@ class AngularSeparationGlobal const quat* equiv_orientations, unsigned int n_equiv_orientations); //! Returns the last computed global angle array - const util::ManagedArray& getAngles() const + std::shared_ptr> getAngles() const { return m_angles; } private: - util::ManagedArray m_angles; //!< Global angle array computed + std::shared_ptr> m_angles; //!< Global angle array computed }; //! Compute the difference in orientation between pairs of points. @@ -60,27 +63,27 @@ class AngularSeparationNeighbor ~AngularSeparationNeighbor() = default; //! Compute the angular separation between neighbors - void compute(const locality::NeighborQuery* nq, const quat* orientations, + void compute(const std::shared_ptr& nq, const quat* orientations, const vec3* query_points, const quat* query_orientations, unsigned int n_query_points, const quat* equiv_orientations, - unsigned int n_equiv_orientations, const freud::locality::NeighborList* nlist, - locality::QueryArgs qargs); + unsigned int n_equiv_orientations, const std::shared_ptr& nlist, + const locality::QueryArgs& qargs); //! Returns the last computed neighbor angle array - const util::ManagedArray& getAngles() const + std::shared_ptr> getAngles() const { return m_angles; } //! Return a pointer to the NeighborList used in the last call to compute. - locality::NeighborList* getNList() + std::shared_ptr getNList() { - return &m_nlist; + return m_nlist; } private: - util::ManagedArray m_angles; //!< neighbor angle array computed - locality::NeighborList m_nlist; //!< The NeighborList used in the last call to compute. + std::shared_ptr> m_angles; //!< neighbor angle array computed + std::shared_ptr m_nlist; //!< The NeighborList used in the last call to compute. }; }; }; // end namespace freud::environment diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index 376182c3a..686287e28 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -6,9 +6,20 @@ #ifdef __SSE2__ #include #endif +#include +#include // NOLINT(modernize-deprecated-headers): Use std::numbers when c++20 is default. +#include +#include +#include "BondHistogramCompute.h" #include "BondOrder.h" -#include "NeighborComputeFunctional.h" +#include "Box.h" +#include "Histogram.h" +#include "ManagedArray.h" +#include "NeighborBond.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" #include "utils.h" /*! \file BondOrder.h @@ -33,8 +44,8 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd /* 0 < \theta < 2PI; 0 < \phi < PI */ - float dt = constants::TWO_PI / float(n_bins_theta); - float dp = M_PI / float(n_bins_phi); + const float dt = constants::TWO_PI / float(n_bins_theta); + const float dp = M_PI / float(n_bins_phi); // this shouldn't be able to happen, but it's always better to check if (dt > constants::TWO_PI) { @@ -46,41 +57,52 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd } // precompute the surface area array - m_sa_array.prepare({n_bins_theta, n_bins_phi}); + // m_sa_array.prepare({n_bins_theta, n_bins_phi}); + m_sa_array = std::make_shared>(std::vector {n_bins_theta, n_bins_phi}); for (unsigned int i = 0; i < n_bins_theta; i++) { for (unsigned int j = 0; j < n_bins_phi; j++) { - float phi = (float) j * dp; - float sa = dt * (std::cos(phi) - std::cos(phi + dp)); - m_sa_array(i, j) = sa; + const float phi = (float) j * dp; + const float sa = dt * (std::cos(phi) - std::cos(phi + dp)); + (*m_sa_array)(i, j) = sa; } } const auto axes = util::Axes {std::make_shared(n_bins_theta, 0, constants::TWO_PI), std::make_shared(n_bins_phi, 0, M_PI)}; m_histogram = BondHistogram(axes); + m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); m_local_histograms = BondHistogram::ThreadLocalHistogram(m_histogram); } -void BondOrder::reduce() +void BondOrder::reset() { - m_histogram.prepare(m_histogram.shape()); - m_bo_array.prepare(m_histogram.shape()); + BondHistogramCompute::reset(); + m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); +} +void BondOrder::reduce() +{ m_histogram.reduceOverThreadsPerBin(m_local_histograms, [&](size_t i) { - m_bo_array[i] = m_histogram[i] / m_sa_array[i] / static_cast(m_frame_counter); + (*m_bo_array)[i] = float(m_histogram[i]) / (*m_sa_array)[i] / static_cast(m_frame_counter); }); } -const util::ManagedArray& BondOrder::getBondOrder() +std::vector> BondOrder::getBinCenters() +{ + return m_histogram.getBinCenters(); +} + +std::shared_ptr> BondOrder::getBondOrder() { return reduceAndReturn(m_bo_array); } -void BondOrder::accumulate(const locality::NeighborQuery* neighbor_query, quat* orientations, - vec3* query_points, quat* query_orientations, - unsigned int n_query_points, const freud::locality::NeighborList* nlist, +void BondOrder::accumulate(const std::shared_ptr& neighbor_query, + quat* orientations, vec3* query_points, + quat* query_orientations, unsigned int n_query_points, + const std::shared_ptr& nlist, freud::locality::QueryArgs qargs) { accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs, @@ -121,7 +143,7 @@ void BondOrder::accumulate(const locality::NeighborQuery* neighbor_query, quat +#include + #include "BondHistogramCompute.h" -#include "Box.h" -#include "Histogram.h" #include "ManagedArray.h" #include "NeighborList.h" #include "NeighborQuery.h" @@ -40,14 +41,17 @@ class BondOrder : public locality::BondHistogramCompute ~BondOrder() override = default; //! Accumulate the bond order - void accumulate(const locality::NeighborQuery* neighbor_query, quat* orientations, + void accumulate(const std::shared_ptr& neighbor_query, quat* orientations, vec3* query_points, quat* query_orientations, unsigned int n_query_points, - const freud::locality::NeighborList* nlist, freud::locality::QueryArgs qargs); + const std::shared_ptr& nlist, + freud::locality::QueryArgs qargs); void reduce() override; + void reset() override; + std::vector> getBinCenters(); - //! Get a reference to the last computed bond order - const util::ManagedArray& getBondOrder(); + //! Get a shared_ptr to the last computed bond order + std::shared_ptr> getBondOrder(); BondOrderMode getMode() const { @@ -55,9 +59,9 @@ class BondOrder : public locality::BondHistogramCompute } private: - util::ManagedArray m_bo_array; //!< bond order array computed - util::ManagedArray m_sa_array; //!< surface area array computed - BondOrderMode m_mode; //!< The mode to calculate with. + std::shared_ptr> m_bo_array; //!< bond order array computed + std::shared_ptr> m_sa_array; //!< surface area array computed + BondOrderMode m_mode; //!< The mode to calculate with. }; }; }; // end namespace freud::environment diff --git a/freud/environment/LocalBondProjection.cc b/freud/environment/LocalBondProjection.cc index 33d8e3f12..c98277ab1 100644 --- a/freud/environment/LocalBondProjection.cc +++ b/freud/environment/LocalBondProjection.cc @@ -1,8 +1,17 @@ // Copyright (c) 2010-2024 The Regents of the University of Michigan // This file is from the freud project, released under the BSD 3-Clause License. +#include +#include +#include + #include "LocalBondProjection.h" +#include "ManagedArray.h" #include "NeighborComputeFunctional.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" +#include "utils.h" /*! \file LocalBondProjection.h \brief Compute the projection of nearest neighbor bonds for each particle onto some @@ -22,7 +31,7 @@ namespace freud { namespace environment { float computeMaxProjection(const vec3& proj_vec, const vec3& local_bond, const quat* equiv_qs, unsigned int n_equiv_qs) { - quat qconst = equiv_qs[0]; + const quat qconst = equiv_qs[0]; // start with the reference vector before it has been rotated by equivalent quaternions float max_proj = dot(proj_vec, local_bond); @@ -30,12 +39,12 @@ float computeMaxProjection(const vec3& proj_vec, const vec3& local // loop through all equivalent rotations and see if they have a larger projection onto local_bond for (unsigned int i = 0; i < n_equiv_qs; i++) { - quat qe = equiv_qs[i]; + const quat qe = equiv_qs[i]; // here we undo a rotation represented by one of the equivalent orientations - quat qtest = conj(qconst) * qe; - vec3 equiv_proj_vec = rotate(qtest, proj_vec); + const quat qtest = conj(qconst) * qe; + const vec3 equiv_proj_vec = rotate(qtest, proj_vec); - float proj_test = dot(equiv_proj_vec, local_bond); + const float proj_test = dot(equiv_proj_vec, local_bond); if (proj_test > max_proj) { @@ -46,45 +55,50 @@ float computeMaxProjection(const vec3& proj_vec, const vec3& local return max_proj; } -void LocalBondProjection::compute(const locality::NeighborQuery* nq, const quat* orientations, - const vec3* query_points, unsigned int n_query_points, - const vec3* proj_vecs, unsigned int n_proj, - const quat* equiv_orientations, unsigned int n_equiv_orientations, - const freud::locality::NeighborList* nlist, locality::QueryArgs qargs) +void LocalBondProjection::compute(const std::shared_ptr& nq, + const quat* orientations, const vec3* query_points, + unsigned int n_query_points, const vec3* proj_vecs, + unsigned int n_proj, const quat* equiv_orientations, + unsigned int n_equiv_orientations, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs) { // This function requires a NeighborList object, so we always make one and store it locally. m_nlist = locality::makeDefaultNlist(nq, nlist, query_points, n_query_points, qargs); // Get the maximum total number of bonds in the neighbor list - const unsigned int tot_num_neigh = m_nlist.getNumBonds(); + const unsigned int tot_num_neigh = m_nlist->getNumBonds(); + const auto& neighbors = m_nlist->getNeighbors(); - m_local_bond_proj.prepare({tot_num_neigh, n_proj}); - m_local_bond_proj_norm.prepare({tot_num_neigh, n_proj}); + m_local_bond_proj + = std::make_shared>(std::vector {tot_num_neigh, n_proj}); + m_local_bond_proj_norm + = std::make_shared>(std::vector {tot_num_neigh, n_proj}); // compute the order parameter util::forLoopWrapper(0, n_query_points, [&](size_t begin, size_t end) { - size_t bond(m_nlist.find_first_index(begin)); + size_t bond(m_nlist->find_first_index(begin)); for (size_t i = begin; i < end; ++i) { - for (; bond < tot_num_neigh && m_nlist.getNeighbors()(bond, 0) == i; ++bond) + for (; bond < tot_num_neigh && (*neighbors)(bond, 0) == i; ++bond) { - const size_t j(m_nlist.getNeighbors()(bond, 1)); + const size_t j((*neighbors)(bond, 1)); // compute bond vector between the two particles - vec3 local_bond(m_nlist.getVectors()[bond]); + vec3 local_bond((*m_nlist->getVectors())(bond)); // rotate bond vector into the local frame of particle p local_bond = rotate(conj(orientations[j]), local_bond); // store the length of this local bond - float local_bond_len = m_nlist.getDistances()[bond]; + const float local_bond_len = (*m_nlist->getDistances())(bond); for (unsigned int k = 0; k < n_proj; k++) { - vec3 proj_vec = proj_vecs[k]; - float max_proj = computeMaxProjection(proj_vec, local_bond, equiv_orientations, - n_equiv_orientations); - m_local_bond_proj(bond, k) = max_proj; - m_local_bond_proj_norm(bond, k) = max_proj / local_bond_len; + const vec3 proj_vec = proj_vecs[k]; + const float max_proj = computeMaxProjection(proj_vec, local_bond, equiv_orientations, + n_equiv_orientations); + (*m_local_bond_proj)(bond, k) = max_proj; + (*m_local_bond_proj_norm)(bond, k) = max_proj / local_bond_len; } } } diff --git a/freud/environment/LocalBondProjection.h b/freud/environment/LocalBondProjection.h index a016e1b10..000e79c32 100644 --- a/freud/environment/LocalBondProjection.h +++ b/freud/environment/LocalBondProjection.h @@ -4,7 +4,8 @@ #ifndef LOCAL_BOND_PROJECTION_H #define LOCAL_BOND_PROJECTION_H -#include "Box.h" +#include + #include "ManagedArray.h" #include "NeighborList.h" #include "NeighborQuery.h" @@ -32,35 +33,37 @@ class LocalBondProjection ~LocalBondProjection() = default; //! Compute the maximal local bond projection - void compute(const locality::NeighborQuery* nq, const quat* orientations, + void compute(const std::shared_ptr& nq, const quat* orientations, const vec3* query_points, unsigned int n_query_points, const vec3* proj_vecs, unsigned int n_proj, const quat* equiv_orientations, - unsigned int n_equiv_orientations, const freud::locality::NeighborList* nlist, - locality::QueryArgs qargs); + unsigned int n_equiv_orientations, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs); //! Get a reference to the last computed maximal local bond projection array - const util::ManagedArray& getProjections() const + std::shared_ptr> getProjections() const { return m_local_bond_proj; } //! Get a reference to the last computed normalized maximal local bond projection array - const util::ManagedArray& getNormedProjections() const + std::shared_ptr> getNormedProjections() const { return m_local_bond_proj_norm; } //! Return a pointer to the NeighborList used in the last call to compute. - locality::NeighborList* getNList() + std::shared_ptr getNList() { - return &m_nlist; + return m_nlist; } private: - locality::NeighborList m_nlist; //!< The NeighborList used in the last call to compute. + std::shared_ptr m_nlist; //!< The NeighborList used in the last call to compute. - util::ManagedArray m_local_bond_proj; //!< Local bond projection array computed - util::ManagedArray m_local_bond_proj_norm; //!< Normalized local bond projection array computed + std::shared_ptr> m_local_bond_proj; //!< Local bond projection array computed + std::shared_ptr> + m_local_bond_proj_norm; //!< Normalized local bond projection array computed }; }; }; // end namespace freud::environment diff --git a/freud/environment/LocalDescriptors.cc b/freud/environment/LocalDescriptors.cc index 1ddfee993..9e5cab83c 100644 --- a/freud/environment/LocalDescriptors.cc +++ b/freud/environment/LocalDescriptors.cc @@ -1,11 +1,25 @@ // Copyright (c) 2010-2024 The Regents of the University of Michigan // This file is from the freud project, released under the BSD 3-Clause License. +#include +#include +#include +#include +#include +#include // NOLINT(modernize-deprecated-headers): Use std::numbers when c++20 is default. +#include +#include #include +#include "Box.h" #include "LocalDescriptors.h" +#include "ManagedArray.h" #include "NeighborComputeFunctional.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" #include "diagonalize.h" +#include "utils.h" /*! \file LocalDescriptors.cc \brief Computes local descriptors. @@ -18,10 +32,11 @@ LocalDescriptors::LocalDescriptors(unsigned int l_max, bool negative_m, : m_l_max(l_max), m_negative_m(negative_m), m_nSphs(0), m_orientation(orientation) {} -void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3* query_points, - unsigned int n_query_points, const quat* orientations, - const freud::locality::NeighborList* nlist, locality::QueryArgs qargs, - unsigned int max_num_neighbors) +void LocalDescriptors::compute(const std::shared_ptr& nq, + const vec3* query_points, unsigned int n_query_points, + const quat* orientations, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs, unsigned int max_num_neighbors) { // This function requires a NeighborList object, so we always make one and store it locally. m_nlist = locality::makeDefaultNlist(nq, nlist, query_points, n_query_points, qargs); @@ -30,14 +45,15 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3::max(); } - m_sphArray.prepare({m_nlist.getNumBonds(), getSphWidth()}); + m_sphArray = std::make_shared>>( + std::vector {m_nlist->getNumBonds(), getSphWidth()}); util::forLoopWrapper(0, nq->getNPoints(), [&](size_t begin, size_t end) { fsph::PointSPHEvaluator sph_eval(m_l_max); for (size_t i = begin; i < end; ++i) { - size_t bond(m_nlist.find_first_index(i)); + size_t bond(m_nlist->find_first_index(i)); unsigned int neighbor_count(0); vec3 rotation_0; @@ -47,12 +63,11 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3 inertiaTensor = util::ManagedArray({3, 3}); - - for (size_t bond_copy(bond); bond_copy < m_nlist.getNumBonds() - && m_nlist.getNeighbors()(bond_copy, 0) == i && neighbor_count < max_num_neighbors; + for (size_t bond_copy(bond); bond_copy < m_nlist->getNumBonds() + && (*m_nlist->getNeighbors())(bond_copy, 0) == i && neighbor_count < max_num_neighbors; ++bond_copy, ++neighbor_count) { - const vec3 r_ij(m_nlist.getVectors()[bond_copy]); + const vec3 r_ij((*m_nlist->getVectors())(bond_copy)); const float r_sq(dot(r_ij, r_ij)); for (size_t ii(0); ii < 3; ++ii) @@ -99,13 +114,13 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3getNumBonds() && (*m_nlist->getNeighbors())(bond, 0) == i && neighbor_count < max_num_neighbors; ++bond, ++neighbor_count) { const unsigned int sphCount(bond * getSphWidth()); - const vec3 r_ij(m_nlist.getVectors()[bond]); - const float magR(m_nlist.getDistances()[bond]); + const vec3 r_ij((*m_nlist->getVectors())(bond)); + const float magR((*m_nlist->getDistances())(bond)); const vec3 bond_ij(dot(rotation_0, r_ij), dot(rotation_1, r_ij), dot(rotation_2, r_ij)); @@ -124,14 +139,15 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3data() + sphCount); } } }); // save the last computed number of particles - m_nSphs = m_nlist.getNumBonds(); + m_nSphs = m_nlist->getNumBonds(); } }; }; // end namespace freud::environment diff --git a/freud/environment/LocalDescriptors.h b/freud/environment/LocalDescriptors.h index 153a6425d..918988fa4 100644 --- a/freud/environment/LocalDescriptors.h +++ b/freud/environment/LocalDescriptors.h @@ -5,8 +5,8 @@ #define LOCAL_DESCRIPTORS_H #include +#include -#include "Box.h" #include "ManagedArray.h" #include "NeighborList.h" #include "NeighborQuery.h" @@ -52,13 +52,13 @@ class LocalDescriptors //! Compute the local neighborhood descriptors given some //! positions and the number of particles - void compute(const locality::NeighborQuery* nq, const vec3* query_points, + void compute(const std::shared_ptr& nq, const vec3* query_points, unsigned int n_query_points, const quat* orientations, - const freud::locality::NeighborList* nlist, locality::QueryArgs qargs, - unsigned int max_num_neighbors = 0); + const std::shared_ptr& nlist, const locality::QueryArgs& qargs, + unsigned int max_num_neighbor); //! Get a reference to the last computed spherical harmonic array - const util::ManagedArray>& getSph() const + std::shared_ptr>> getSph() const { return m_sphArray; } @@ -70,9 +70,9 @@ class LocalDescriptors } //! Return a pointer to the NeighborList used in the last call to compute. - locality::NeighborList* getNList() + std::shared_ptr getNList() { - return &m_nlist; + return m_nlist; } bool getNegativeM() const @@ -86,14 +86,14 @@ class LocalDescriptors } private: - unsigned int m_l_max; //!< Maximum spherical harmonic l to calculate - bool m_negative_m; //!< true if we should compute Ylm for negative m - unsigned int m_nSphs; //!< Last number of bond spherical harmonics computed - locality::NeighborList m_nlist; //!< The NeighborList used in the last call to compute. - LocalDescriptorOrientation m_orientation; //!< The orientation mode to compute with. + unsigned int m_l_max; //!< Maximum spherical harmonic l to calculate + bool m_negative_m; //!< true if we should compute Ylm for negative m + unsigned int m_nSphs; //!< Last number of bond spherical harmonics computed + std::shared_ptr m_nlist; //!< The NeighborList used in the last call to compute. + LocalDescriptorOrientation m_orientation; //!< The orientation mode to compute with. //! Spherical harmonics for each neighbor - util::ManagedArray> m_sphArray; + std::shared_ptr>> m_sphArray; }; }; }; // end namespace freud::environment diff --git a/freud/environment/MatchEnv.cc b/freud/environment/MatchEnv.cc index 0a88b51b5..1da138250 100644 --- a/freud/environment/MatchEnv.cc +++ b/freud/environment/MatchEnv.cc @@ -2,12 +2,24 @@ // This file is from the freud project, released under the BSD 3-Clause License. #include +#include +#include +#include #include #include +#include +#include +#include +#include "BiMap.h" +#include "Box.h" +#include "ManagedArray.h" #include "MatchEnv.h" - #include "NeighborComputeFunctional.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "Registration.h" +#include "VectorMath.h" namespace freud { namespace environment { @@ -23,9 +35,9 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, if (rank[s[a].env_ind] == rank[s[b].env_ind]) { // Get the ENTIRE set that corresponds to head_b. - unsigned int head_b = find(b); - std::vector m_set = findSet(head_b); - for (unsigned int node : m_set) + const unsigned int head_b = find(b); + const std::vector m_set = findSet(head_b); + for (const unsigned int node : m_set) { // Go through the entire tree/set. // Make a copy of the old set of vector indices for this @@ -39,7 +51,7 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, // and set it properly. for (unsigned int proper_a_ind = 0; proper_a_ind < vec_map.size(); proper_a_ind++) { - unsigned int proper_b_ind = vec_map.left[proper_a_ind]; + const unsigned int proper_b_ind = vec_map.left[proper_a_ind]; // old_node_vec_ind[proper_b_ind] is "relative_b_ind" s[node].vec_ind[proper_a_ind] = old_node_vec_ind[proper_b_ind]; @@ -62,9 +74,9 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, if (rank[s[a].env_ind] > rank[s[b].env_ind]) { // Get the ENTIRE set that corresponds to head_b. - unsigned int head_b = find(b); - std::vector m_set = findSet(head_b); - for (unsigned int node : m_set) + const unsigned int head_b = find(b); + const std::vector m_set = findSet(head_b); + for (const unsigned int node : m_set) { // Go through the entire tree/set. // Make a copy of the old set of vector indices for this @@ -78,7 +90,7 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, // and set it properly. for (unsigned int proper_a_ind = 0; proper_a_ind < vec_map.size(); proper_a_ind++) { - unsigned int proper_b_ind = vec_map.left[proper_a_ind]; + const unsigned int proper_b_ind = vec_map.left[proper_a_ind]; // old_node_vec_ind[proper_b_ind] is "relative_b_ind" s[node].vec_ind[proper_a_ind] = old_node_vec_ind[proper_b_ind]; @@ -96,11 +108,11 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, } else { - rotmat3 rotationT = transpose(rotation); + const rotmat3 rotationT = transpose(rotation); // Get the ENTIRE set that corresponds to head_a. - unsigned int head_a = find(a); - std::vector m_set = findSet(head_a); - for (unsigned int node : m_set) + const unsigned int head_a = find(a); + const std::vector m_set = findSet(head_a); + for (const unsigned int node : m_set) { // Go through the entire tree/set. // Make a copy of the old set of vector indices for this @@ -114,7 +126,7 @@ void EnvDisjointSet::merge(const unsigned int a, const unsigned int b, // and set it properly. for (unsigned int proper_b_ind = 0; proper_b_ind < vec_map.size(); proper_b_ind++) { - unsigned int proper_a_ind = vec_map.right[proper_b_ind]; + const unsigned int proper_a_ind = vec_map.right[proper_b_ind]; // old_node_vec_ind[proper_a_ind] is "relative_a_ind" s[node].vec_ind[proper_b_ind] = old_node_vec_ind[proper_a_ind]; @@ -151,7 +163,7 @@ unsigned int EnvDisjointSet::find(const unsigned int c) unsigned int i = c; while (i != r) { - unsigned int j = s[i].env_ind; + const unsigned int j = s[i].env_ind; s[i].env_ind = r; i = j; } @@ -168,7 +180,7 @@ std::vector EnvDisjointSet::findSet(const unsigned int m) for (unsigned int i = 0; i < s.size(); i++) { // get the head environment index - unsigned int head_env = find(s[i].env_ind); + const unsigned int head_env = find(s[i].env_ind); // if we are part of the environment m, add the vectors to m_set if (head_env == m) { @@ -180,7 +192,7 @@ std::vector EnvDisjointSet::findSet(const unsigned int m) if (invalid_ind) { std::ostringstream msg; - msg << "Index " << m << " must be a head index in the environment set!" << std::endl; + msg << "Index " << m << " must be a head index in the environment set!\n"; throw std::invalid_argument(msg.str()); } @@ -201,7 +213,7 @@ std::vector> EnvDisjointSet::getAvgEnv(const unsigned int m) if (!i.ghost) { // get the head environment index - unsigned int head_env = find(i.env_ind); + const unsigned int head_env = find(i.env_ind); // if we are part of the environment m, add the vectors to env if (head_env == m) { @@ -209,8 +221,8 @@ std::vector> EnvDisjointSet::getAvgEnv(const unsigned int m) // add them to env for (unsigned int proper_ind = 0; proper_ind < i.vecs.size(); proper_ind++) { - unsigned int relative_ind = i.vec_ind[proper_ind]; - vec3 proper_vec = i.proper_rot * i.vecs[relative_ind]; + const unsigned int relative_ind = i.vec_ind[proper_ind]; + const vec3 proper_vec = i.proper_rot * i.vecs[relative_ind]; if (proper_ind < env.size()) { env[proper_ind] += proper_vec; @@ -229,7 +241,7 @@ std::vector> EnvDisjointSet::getAvgEnv(const unsigned int m) if (invalid_ind) { std::ostringstream msg; - msg << "Index " << m << " must be a head index in the environment set!" << std::endl; + msg << "Index " << m << " must be a head index in the environment set!\n"; throw std::invalid_argument(msg.str()); } @@ -246,7 +258,7 @@ std::vector> EnvDisjointSet::getIndividualEnv(const unsigned int m) if (m >= s.size()) { std::ostringstream msg; - msg << "Index " << m << " must be less than the size of the environment set!" << std::endl; + msg << "Index " << m << " must be less than the size of the environment set!\n"; throw std::invalid_argument(msg.str()); } @@ -256,7 +268,7 @@ std::vector> EnvDisjointSet::getIndividualEnv(const unsigned int m) // 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]; + const unsigned int relative_ind = s[m].vec_ind[proper_ind]; env.push_back(s[m].proper_rot * s[m].vecs[relative_ind]); } @@ -308,8 +320,8 @@ std::pair, BiMap> isSimilar(Environme // minimal RMSD, as best as it can figure out. // Does this vector mapping pass the more stringent criterion // imposed by the threshold? - vec3 delta = v1[registered_pair->first] - v2[registered_pair->second]; - float r_sq = dot(delta, delta); + const vec3 delta = v1[registered_pair->first] - v2[registered_pair->second]; + const float r_sq = dot(delta, delta); if (r_sq < threshold_sq) { vec_map.emplace(registered_pair->first, registered_pair->second); @@ -324,8 +336,8 @@ std::pair, BiMap> isSimilar(Environme { for (unsigned int j = 0; j < e2.vecs.size(); j++) { - vec3 delta = v1[i] - v2[j]; - float r_sq = dot(delta, delta); + const vec3 delta = v1[i] - v2[j]; + const float r_sq = dot(delta, delta); if (r_sq < threshold_sq) { // these vectors are deemed "matching" @@ -344,7 +356,7 @@ std::pair, BiMap> isSimilar(Environme return std::pair, BiMap>(rotation, vec_map); } // otherwise, return an empty bimap - BiMap empty_map; + const BiMap empty_map; return std::pair, BiMap>(rotation, empty_map); } @@ -357,9 +369,9 @@ std::map isSimilar(const box::Box& box, const vec3, BiMap> mapping + const std::pair, BiMap> mapping = isSimilar(e0, e1, threshold_sq, registration); - rotmat3 rotation = mapping.first; + const rotmat3 rotation = mapping.first; BiMap vec_map = mapping.second; // update refPoints2 in case registration has taken place @@ -392,8 +404,8 @@ std::pair makeEnvironments(const box::Box& box, const // be wrapped into the box as well. for (unsigned int i = 0; i < numRef; i++) { - vec3 p0 = box.wrap(refPoints1[i]); - vec3 p1 = box.wrap(refPoints2[i]); + const vec3 p0 = box.wrap(refPoints1[i]); + const vec3 p1 = box.wrap(refPoints2[i]); e0.addVec(p0); e1.addVec(p1); } @@ -459,9 +471,9 @@ std::map minimizeRMSD(const box::Box& box, const vec std::tie(e0, e1) = makeEnvironments(box, refPoints1, refPoints2, numRef); float tmp_min_rmsd = -1.0; - std::pair, BiMap> mapping + const std::pair, BiMap> mapping = minimizeRMSD(e0, e1, tmp_min_rmsd, registration); - rotmat3 rotation = mapping.first; + const rotmat3 rotation = mapping.first; BiMap vec_map = mapping.second; min_rmsd = tmp_min_rmsd; @@ -488,20 +500,20 @@ MatchEnv::~MatchEnv() = default; EnvironmentCluster::~EnvironmentCluster() = default; -Environment MatchEnv::buildEnv(const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond, - unsigned int i, unsigned int env_ind) +Environment MatchEnv::buildEnv(const std::shared_ptr& nlist, size_t num_bonds, + size_t& bond, unsigned int i, unsigned int env_ind) { Environment ei = Environment(); // set the environment index equal to the particle index ei.env_ind = env_ind; - for (; bond < num_bonds && nlist->getNeighbors()(bond, 0) == i; ++bond) + for (; bond < num_bonds && (*nlist->getNeighbors())(bond, 0) == i; ++bond) { // compute vec{r} between the two particles - const size_t j(nlist->getNeighbors()(bond, 1)); + const size_t j((*nlist->getNeighbors())(bond, 1)); if (i != j) { - vec3 delta(nlist->getVectors()[bond]); + const vec3 delta((*nlist->getVectors())[bond]); ei.addVec(delta); } } @@ -509,25 +521,27 @@ Environment MatchEnv::buildEnv(const freud::locality::NeighborList* nlist, size_ return ei; } -void EnvironmentCluster::compute(const freud::locality::NeighborQuery* nq, - const freud::locality::NeighborList* nlist_arg, locality::QueryArgs qargs, - const freud::locality::NeighborList* env_nlist_arg, +void EnvironmentCluster::compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, + locality::QueryArgs qargs, + const std::shared_ptr& env_nlist_arg, locality::QueryArgs env_qargs, float threshold, bool registration) { - const locality::NeighborList nlist + const std::shared_ptr nlist = locality::makeDefaultNlist(nq, nlist_arg, nq->getPoints(), nq->getNPoints(), qargs); - const locality::NeighborList env_nlist + const std::shared_ptr env_nlist = locality::makeDefaultNlist(nq, env_nlist_arg, nq->getPoints(), nq->getNPoints(), env_qargs); - unsigned int Np = nq->getNPoints(); - m_env_index.prepare(Np); + const unsigned int Np = nq->getNPoints(); + m_env_index = std::make_shared>(std::vector {Np}); + // m_env_index.prepare(Np); - float m_threshold_sq = threshold * threshold; + const float m_threshold_sq = threshold * threshold; - nlist.validate(Np, Np); - env_nlist.validate(Np, Np); + nlist->validate(Np, Np); + env_nlist->validate(Np, Np); size_t env_bond(0); - const size_t env_num_bonds(env_nlist.getNumBonds()); + const size_t env_num_bonds(env_nlist->getNumBonds()); // create a disjoint set where all particles belong in their own cluster EnvDisjointSet dj(Np); @@ -538,7 +552,7 @@ void EnvironmentCluster::compute(const freud::locality::NeighborQuery* nq, // if you don't do this, things will get screwy. for (unsigned int i = 0; i < Np; i++) { - Environment ei = buildEnv(&env_nlist, env_num_bonds, env_bond, i, i); + const Environment ei = buildEnv(env_nlist, env_num_bonds, env_bond, i, i); dj.s.push_back(ei); } @@ -547,21 +561,21 @@ void EnvironmentCluster::compute(const freud::locality::NeighborQuery* nq, for (unsigned int i = 0; i < Np; i++) { // loop over the neighbors - for (; bond < nlist.getNumBonds() && nlist.getNeighbors()(bond, 0) == i; ++bond) + for (; bond < nlist->getNumBonds() && (*nlist->getNeighbors())(bond, 0) == i; ++bond) { - const size_t j(nlist.getNeighbors()(bond, 1)); - std::pair, BiMap> mapping + const size_t j((*nlist->getNeighbors())(bond, 1)); + const std::pair, BiMap> mapping = isSimilar(dj.s[i], dj.s[j], m_threshold_sq, registration); rotmat3 rotation = mapping.first; - BiMap vec_map = mapping.second; + const BiMap vec_map = mapping.second; // if the mapping between the vectors of the environments // is NOT empty, then the environments are similar, so // merge them. if (!vec_map.empty()) { // merge the two sets using the disjoint set - unsigned int a = dj.find(i); - unsigned int b = dj.find(j); + const unsigned int a = dj.find(i); + const unsigned int b = dj.find(j); if (a != b) { dj.merge(i, j, vec_map, rotation); @@ -590,15 +604,15 @@ unsigned int EnvironmentCluster::populateEnv(EnvDisjointSet dj) if (!dj.s[i].ghost) { // grab the set of vectors that define this individual environment - std::vector> part_vecs = dj.getIndividualEnv(i); + const std::vector> part_vecs = dj.getIndividualEnv(i); - unsigned int c = dj.find(i); + const unsigned int c = dj.find(i); // insert the set into the mapping if we haven't seen it before. // also grab the vectors that define the set and insert them into cluster_env if (label_map.count(c) == 0) { label_map[c] = cur_set; - std::vector> vecs = dj.getAvgEnv(c); + const std::vector> vecs = dj.getAvgEnv(c); label_ind = label_map[c]; cluster_env[label_ind] = vecs; cur_set++; @@ -609,7 +623,7 @@ unsigned int EnvironmentCluster::populateEnv(EnvDisjointSet dj) } // label this particle in m_env_index - m_env_index[particle_ind] = label_ind; + (*m_env_index)[particle_ind] = label_ind; m_point_environments.emplace_back(); for (const auto& part_vec : part_vecs) @@ -634,18 +648,18 @@ unsigned int EnvironmentCluster::populateEnv(EnvDisjointSet dj) /************************* * EnvironmentMotifMatch * *************************/ -void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq, - const freud::locality::NeighborList* nlist_arg, locality::QueryArgs qargs, - const vec3* motif, unsigned int motif_size, float threshold, - bool registration) +void EnvironmentMotifMatch::compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, + locality::QueryArgs qargs, const vec3* motif, + unsigned int motif_size, float threshold, bool registration) { - const locality::NeighborList nlist + const std::shared_ptr nlist = locality::makeDefaultNlist(nq, nlist_arg, nq->getPoints(), nq->getNPoints(), qargs); - unsigned int Np = nq->getNPoints(); - float m_threshold_sq = threshold * threshold; + const unsigned int Np = nq->getNPoints(); + const float m_threshold_sq = threshold * threshold; - nlist.validate(Np, Np); + nlist->validate(Np, Np); // create a disjoint set where all particles belong in their own cluster. // this has to have ONE MORE environment than there are actual particles, @@ -663,7 +677,7 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq, // be wrapped into the box as well. for (unsigned int i = 0; i < motif_size; i++) { - vec3 p = nq->getBox().wrap(motif[i]); + const vec3 p = nq->getBox().wrap(motif[i]); e0.addVec(p); } @@ -671,9 +685,10 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq, dj.s.push_back(e0); size_t bond(0); - const size_t num_bonds(nlist.getNumBonds()); + const size_t num_bonds(nlist->getNumBonds()); - m_matches.prepare(Np); + m_matches = std::make_shared>(std::vector {Np}); + // m_matches.prepare(Np); // loop through the particles and add their environments to the set // take care, here: set things up s.t. the env_ind of every environment @@ -681,24 +696,24 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq, // if you don't do this, things will get screwy. for (unsigned int i = 0; i < Np; i++) { - unsigned int dummy = i + 1; - Environment ei = buildEnv(&nlist, num_bonds, bond, i, dummy); + const unsigned int dummy = i + 1; + const Environment ei = buildEnv(nlist, num_bonds, bond, i, dummy); dj.s.push_back(ei); // if the environment matches e0, merge it into the e0 environment set - std::pair, BiMap> mapping + const std::pair, BiMap> mapping = isSimilar(dj.s[0], dj.s[dummy], m_threshold_sq, registration); rotmat3 rotation = mapping.first; - BiMap vec_map = mapping.second; + const BiMap vec_map = mapping.second; // if the mapping between the vectors of the environments is NOT empty, // then the environments are similar. if (!vec_map.empty()) { dj.merge(0, dummy, vec_map, rotation); - m_matches[i] = true; + (*m_matches)[i] = 1; } // grab the set of vectors that define this individual environment - std::vector> part_vecs = dj.getIndividualEnv(dummy); + const std::vector> part_vecs = dj.getIndividualEnv(dummy); m_point_environments.emplace_back(); for (const auto& part_vec : part_vecs) @@ -711,15 +726,15 @@ void EnvironmentMotifMatch::compute(const freud::locality::NeighborQuery* nq, /**************************** * EnvironmentRMSDMinimizer * ****************************/ -void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq, - const freud::locality::NeighborList* nlist_arg, +void EnvironmentRMSDMinimizer::compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, locality::QueryArgs qargs, const vec3* motif, unsigned int motif_size, bool registration) { - const locality::NeighborList nlist + const std::shared_ptr nlist = locality::makeDefaultNlist(nq, nlist_arg, nq->getPoints(), nq->getNPoints(), qargs); - unsigned int Np = nq->getNPoints(); + const unsigned int Np = nq->getNPoints(); // create a disjoint set where all particles belong in their own cluster. // this has to have ONE MORE environment than there are actual particles, @@ -737,7 +752,7 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq, // be wrapped into the box as well. for (unsigned int i = 0; i < motif_size; i++) { - vec3 p = nq->getBox().wrap(motif[i]); + const vec3 p = nq->getBox().wrap(motif[i]); e0.addVec(p); } @@ -745,9 +760,9 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq, dj.s.push_back(e0); size_t bond(0); - const size_t num_bonds(nlist.getNumBonds()); - - m_rmsds.prepare(Np); + const size_t num_bonds(nlist->getNumBonds()); + m_rmsds = std::make_shared>(std::vector {Np}); + // m_rmsds.prepare(Np); // loop through the particles and add their environments to the set // take care, here: set things up s.t. the env_ind of every environment @@ -755,18 +770,18 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq, // if you don't do this, things will get screwy. for (unsigned int i = 0; i < Np; i++) { - unsigned int dummy = i + 1; - Environment ei = buildEnv(&nlist, num_bonds, bond, i, dummy); + const unsigned int dummy = i + 1; + const Environment ei = buildEnv(nlist, num_bonds, bond, i, dummy); dj.s.push_back(ei); // if the environment matches e0, merge it into the e0 environment set float min_rmsd = -1.0; - std::pair, BiMap> mapping + const std::pair, BiMap> mapping = minimizeRMSD(dj.s[0], dj.s[dummy], min_rmsd, registration); rotmat3 rotation = mapping.first; - BiMap vec_map = mapping.second; + const BiMap vec_map = mapping.second; // populate the min_rmsd vector - m_rmsds[i] = min_rmsd; + (*m_rmsds)[i] = min_rmsd; // if the mapping between the vectors of the environments is NOT // empty, then the environments are similar. @@ -778,10 +793,10 @@ void EnvironmentRMSDMinimizer::compute(const freud::locality::NeighborQuery* nq, } // grab the set of vectors that define this individual environment - std::vector> part_vecs = dj.getIndividualEnv(dummy); + const std::vector> part_vecs = dj.getIndividualEnv(dummy); m_point_environments.emplace_back(); - for (auto& part_vec : part_vecs) + for (const auto& part_vec : part_vecs) { m_point_environments[i].push_back(part_vec); } diff --git a/freud/environment/MatchEnv.h b/freud/environment/MatchEnv.h index 7ad35463c..4cd39a284 100644 --- a/freud/environment/MatchEnv.h +++ b/freud/environment/MatchEnv.h @@ -4,7 +4,10 @@ #ifndef MATCH_ENV_H #define MATCH_ENV_H +#include #include +#include +#include #include #include "BiMap.h" @@ -12,7 +15,6 @@ #include "ManagedArray.h" #include "NeighborList.h" #include "NeighborQuery.h" -#include "Registration.h" #include "VectorMath.h" /*! \file MatchEnv.h @@ -204,8 +206,8 @@ class MatchEnv //! Construct and return a local environment surrounding the particle indexed by i. Set the environment //! index to env_ind. - static Environment buildEnv(const freud::locality::NeighborList* nlist, size_t num_bonds, size_t& bond, - unsigned int i, unsigned int env_ind); + static Environment buildEnv(const std::shared_ptr& nlist, size_t num_bonds, + size_t& bond, 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 std::vector>>& getPointEnvironments() @@ -263,12 +265,13 @@ class EnvironmentCluster : public MatchEnv * orient the second set of vectors such that it * minimizes the RMSD between the two sets */ - void compute(const freud::locality::NeighborQuery* nq, const freud::locality::NeighborList* nlist_arg, - locality::QueryArgs qargs, const freud::locality::NeighborList* env_nlist_arg, + void compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, locality::QueryArgs qargs, + const std::shared_ptr& env_nlist_arg, locality::QueryArgs env_qargs, float threshold, bool registration = false); //! Get a reference to the particles, indexed into clusters according to their matching local environments - const util::ManagedArray& getClusters() + const std::shared_ptr>& getClusters() { return m_env_index; } @@ -303,8 +306,9 @@ class EnvironmentCluster : public MatchEnv */ unsigned int populateEnv(EnvDisjointSet dj); - unsigned int m_num_clusters {0}; //!< Last number of local environments computed - util::ManagedArray m_env_index; //!< Cluster index determined for each particle + unsigned int m_num_clusters {0}; //!< Last number of local environments computed + std::shared_ptr> + m_env_index; //!< Cluster index determined for each particle std::vector>> m_cluster_environments; //!< Dictionary of (cluster id, vectors) pairs }; @@ -350,18 +354,19 @@ class EnvironmentMotifMatch : public MatchEnv * orient the second set of vectors such that it * minimizes the RMSD between the two sets */ - void compute(const freud::locality::NeighborQuery* nq, const freud::locality::NeighborList* nlist_arg, - locality::QueryArgs qargs, const vec3* motif, unsigned int motif_size, - float threshold, bool registration = false); + void compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, locality::QueryArgs qargs, + const vec3* motif, unsigned int motif_size, float threshold, + bool registration = false); //! Return the array indicating whether each particle matched the motif or not. - const util::ManagedArray& getMatches() + const std::shared_ptr>& getMatches() { return m_matches; } private: - util::ManagedArray + std::shared_ptr> m_matches; //!< Boolean array indicating whether or not a particle's environment matches the motif. }; @@ -412,19 +417,19 @@ class EnvironmentRMSDMinimizer : public MatchEnv * orient the second set of vectors such that it * minimizes the RMSD between the two sets */ - void compute(const freud::locality::NeighborQuery* nq, const freud::locality::NeighborList* nlist_arg, - locality::QueryArgs qargs, const vec3* motif, unsigned int motif_size, - bool registration = false); + void compute(const std::shared_ptr& nq, + const std::shared_ptr& nlist_arg, locality::QueryArgs qargs, + const vec3* motif, unsigned int motif_size, bool registration = false); //! Return the array indicating whether or not a successful mapping was found between each particle and //! the provided motif. - const util::ManagedArray& getRMSDs() + const std::shared_ptr>& getRMSDs() { return m_rmsds; } private: - util::ManagedArray + std::shared_ptr> m_rmsds; //!< Boolean array indicating whether or not a particle's environment matches the motif. }; diff --git a/freud/environment/Registration.h b/freud/environment/Registration.h index fa36db401..a0996b8fd 100644 --- a/freud/environment/Registration.h +++ b/freud/environment/Registration.h @@ -4,10 +4,16 @@ #ifndef REGISTRATION_H #define REGISTRATION_H +#include #include +#include +#include #include #include +#include #include +#include +#include #include #ifdef _WIN32 @@ -18,21 +24,21 @@ #include #endif -#include "Eigen/Eigen/Dense" -#include "Eigen/Eigen/Sparse" +#include #include "BiMap.h" #include "VectorMath.h" namespace freud { namespace environment { +// NOLINTNEXTLINE(misc-include-cleaner) - Eigen's include structure confuses clang-tidy using matrix = Eigen::Matrix; inline matrix makeEigenMatrix(const std::vector>& vecs) { // build the Eigen matrix matrix mat; - unsigned int size = vecs.size(); + const unsigned int size = vecs.size(); // we know the dimension is 3 bc we're dealing with a vector of vec3's. mat.resize(size, 3); for (unsigned int i = 0; i < size; i++) @@ -54,8 +60,7 @@ inline std::vector> makeVec3Matrix(const matrix& m) if (m.cols() != 3) { std::ostringstream msg; - msg << "makeVec3Matrix requires the input matrix to have 3 columns, not " << m.cols() << "!" - << std::endl; + msg << "makeVec3Matrix requires the input matrix to have 3 columns, not " << m.cols() << "!\n"; throw std::invalid_argument(msg.str()); } std::vector> vecs; @@ -98,7 +103,7 @@ inline matrix Rotate(const matrix& R, const matrix& P) { std::ostringstream msg; msg << "Rotation matrix has " << R.cols() << " columns and point matrix has " << P.rows() - << " rows. These must be equal to perform the rotation." << std::endl; + << " rows. These must be equal to perform the rotation.\n"; throw std::invalid_argument(msg.str()); } // Apply the rotation R. @@ -111,14 +116,15 @@ inline matrix Rotate(const matrix& R, const matrix& P) inline void KabschAlgorithm(const matrix& P, const matrix& Q, matrix& Rotation) { // Preconditions: P and Q have been translated to have the same center of mass. - matrix A = P.transpose() * Q; + const matrix A = P.transpose() * Q; // singular value decomposition (~ eigen decomposition) - Eigen::JacobiSVD svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV); + // NOLINTNEXTLINE(misc-include-cleaner) - Eigen's include structure confuses clang-tidy + const Eigen::JacobiSVD svd(A, Eigen::ComputeFullU | Eigen::ComputeFullV); // A = USV^T matrix U = svd.matrixU(); matrix V = svd.matrixV(); - double det = (V * U.transpose()).determinant(); + const double det = (V * U.transpose()).determinant(); // if the rotation as we've found it, rot=VU^T, is IMPROPER, find the next best // (proper) rotation by reflecting the smallest principal axis in rot: @@ -169,13 +175,12 @@ class RegisterBruteForce points = makeEigenMatrix(pts); int num_pts; - int N = points.rows(); + const size_t N = points.rows(); if (N != m_ref_points.rows()) { std::ostringstream msg; msg << "There are " << m_ref_points.rows() << " reference points and " << N << " points. "; - msg << "Brute force matching requires the same number of reference points and points!" - << std::endl; + msg << "Brute force matching requires the same number of reference points and points!\n"; throw std::invalid_argument(msg.str()); } @@ -188,14 +193,14 @@ class RegisterBruteForce int p2 = 0; while (p0 == p1 || p0 == p2 || p1 == p2) { - p0 = rng.random_int(0, N - 1); + p0 = rng.random_int(0, int(N - 1)); if (N == 1) { p1 = -2; } else { - p1 = rng.random_int(0, N - 1); + p1 = rng.random_int(0, int(N - 1)); } if (N == 2 || N == 1) @@ -204,7 +209,7 @@ class RegisterBruteForce } else { - p2 = rng.random_int(0, N - 1); + p2 = rng.random_int(0, int(N - 1)); } } @@ -220,7 +225,7 @@ class RegisterBruteForce p.row(1) = m_ref_points.row(p1); q.resize(num_pts, m_ref_points.cols()); } - else if (N == int(1)) + else if (N == 1) { num_pts = 1; p.resize(num_pts, m_ref_points.cols()); @@ -240,20 +245,20 @@ class RegisterBruteForce { do { - if (N == int(2)) + if (N == 2) { - q.row(0) = points.row(comb[0]); - q.row(1) = points.row(comb[1]); + q.row(0) = points.row(int(comb[0])); + q.row(1) = points.row(int(comb[1])); } - else if (N == int(1)) + else if (N == 1) { - q.row(0) = points.row(comb[0]); + q.row(0) = points.row(int(comb[0])); } else { - q.row(0) = points.row(comb[0]); - q.row(1) = points.row(comb[1]); - q.row(2) = points.row(comb[2]); + q.row(0) = points.row(int(comb[0])); + q.row(1) = points.row(int(comb[1])); + q.row(2) = points.row(int(comb[2])); } // finds the optimal rotation of the FIRST input set @@ -268,7 +273,7 @@ class RegisterBruteForce // feed back in the TRANSPOSE of rot_points such that // the input matrix is (Nx3). BiMap vec_map; - float rmsd = AlignedRMSDTree(rot_points.transpose(), vec_map); + const float rmsd = AlignedRMSDTree(rot_points.transpose(), vec_map); if (rmsd < rmsd_min || rmsd_min < 0.0) { m_rmsd = rmsd; @@ -287,7 +292,7 @@ class RegisterBruteForce } } } while (std::next_permutation(comb, comb + num_pts)); - } while (NextCombination(comb, N, num_pts)); + } while (NextCombination(comb, int(N), num_pts)); } // end for loop over shuffles // The rotation that we've found from the KabschAlgorithm // actually acts on P^T. @@ -299,13 +304,13 @@ class RegisterBruteForce std::vector> getRotation() { - matrix R = m_rotation; + const matrix R = m_rotation; return makeVec3Matrix(R); } std::vector> getTranslation() { - matrix T = m_translation; + const matrix T = m_translation; return makeVec3Matrix(T); } @@ -357,14 +362,14 @@ class RegisterBruteForce for (int r = 0; r < points.rows(); r++) { // get the rotated point - vec3 pfit = make_point(points.row(r)); + const vec3 pfit = make_point(points.row(r)); // compute squared distances to all unused reference points std::vector> ref_distances; for (auto ref_index : unused_indices) { - vec3 ref_point = make_point(m_ref_points.row(ref_index)); - vec3 delta = ref_point - pfit; - float r_sq = dot(delta, delta); + const vec3 ref_point = make_point(m_ref_points.row(ref_index)); + const vec3 delta = ref_point - pfit; + const float r_sq = dot(delta, delta); ref_distances.emplace_back(ref_index, r_sq); } // sort the ref_distances from nearest to farthest @@ -386,11 +391,11 @@ class RegisterBruteForce { if (row.rows() == 2) { - return vec3(row[0], row[1], 0.0); + return vec3(float(row[0]), float(row[1]), 0.0); } if (row.rows() == 3) { - return vec3(row[0], row[1], row[2]); + return vec3(float(row[0]), float(row[1]), float(row[2])); } throw(std::runtime_error("points must 2 or 3 dimensions")); } @@ -401,7 +406,7 @@ class RegisterBruteForce return (a.second < b.second); } - static inline bool NextCombination(size_t* comb, int N, int k) + static bool NextCombination(size_t* comb, int N, int k) { // returns next combination. if (k == 0 || N == 0 || (comb == nullptr)) @@ -437,12 +442,13 @@ class RegisterBruteForce } int random_int(int a, int b) { + // NOLINTNEXTLINE(misc-const-correctness) - const causes a compile error std::uniform_int_distribution distribution(a, b); return distribution(m_generator); } private: - inline void seed_generator(const size_t& n = 100) + void seed_generator(const size_t& n = 100) { std::vector seeds; try @@ -455,7 +461,7 @@ class RegisterBruteForce } catch (...) { - std::cerr << "random_device is not available..." << std::endl; + std::cerr << "random_device is not available...\n"; seeds.push_back(size_t(std::chrono::system_clock::now().time_since_epoch().count())); seeds.push_back(size_t(getpid())); } diff --git a/freud/environment/export-AngularSeparationGlobal.cc b/freud/environment/export-AngularSeparationGlobal.cc new file mode 100644 index 000000000..56d13b323 --- /dev/null +++ b/freud/environment/export-AngularSeparationGlobal.cc @@ -0,0 +1,50 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly + +#include "AngularSeparation.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute(const std::shared_ptr& angular_separation, + const nb_array>& global_orientations, + const nb_array>& orientations, + const nb_array>& equiv_orientations) +{ + unsigned int const n_global = global_orientations.shape(0); + unsigned int const n_points = orientations.shape(0); + unsigned int const n_equiv_orientations = equiv_orientations.shape(0); + auto* global_orientations_data = reinterpret_cast*>(global_orientations.data()); + auto* orientations_data = reinterpret_cast*>(orientations.data()); + auto* equiv_orientations_data = reinterpret_cast*>(equiv_orientations.data()); + angular_separation->compute(global_orientations_data, n_global, orientations_data, n_points, + equiv_orientations_data, n_equiv_orientations); +} + +}; // namespace wrap + +namespace detail { + +void export_AngularSeparationGlobal(nb::module_& module) +{ + nb::class_(module, "AngularSeparationGlobal") + .def(nb::init<>()) + .def("getAngles", &AngularSeparationGlobal::getAngles) + .def("compute", &wrap::compute, nb::arg("global_orientations"), nb::arg("orientations"), + nb::arg("equiv_orientations")); +} + +}; // namespace detail +}; }; // namespace freud::environment diff --git a/freud/environment/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc new file mode 100644 index 000000000..b2175705a --- /dev/null +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -0,0 +1,58 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly + +#include "AngularSeparation.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { + +void compute(std::shared_ptr& angular_separation, + const std::shared_ptr& nq, + nb_array>& orientations, + nb_array>& query_points, + nb_array>& query_orientations, + nb_array>& equiv_orientations, + const std::shared_ptr& nlist, const locality::QueryArgs& qargs) +{ + auto* orientations_data = reinterpret_cast*>(orientations.data()); + auto* query_points_data = reinterpret_cast*>(query_points.data()); + auto* query_orientations_data = reinterpret_cast*>(query_orientations.data()); + auto* equiv_orientations_data = reinterpret_cast*>(equiv_orientations.data()); + const unsigned int n_equiv_orientations = equiv_orientations.shape(0); + const unsigned int n_query_points = query_orientations.shape(0); + angular_separation->compute(nq, orientations_data, query_points_data, query_orientations_data, + n_query_points, equiv_orientations_data, n_equiv_orientations, nlist, qargs); +} +}; // namespace wrap + +namespace detail { + +void export_AngularSeparationNeighbor(nb::module_& module) +{ + nb::class_(module, "AngularSeparationNeighbor") + .def(nb::init<>()) + .def("getNList", &AngularSeparationNeighbor::getNList) + .def("getAngles", &AngularSeparationNeighbor::getAngles) + .def("compute", &wrap::compute, nb::arg("nq"), nb::arg("orientations"), nb::arg("query_points"), + nb::arg("query_orientations"), nb::arg("equiv_orientations"), nb::arg("nlist").none(), + nb::arg("qargs")); +} + +}; // namespace detail + +}; }; // namespace freud::environment diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc new file mode 100644 index 000000000..c3358d26b --- /dev/null +++ b/freud/environment/export-BondOrder.cc @@ -0,0 +1,81 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly + +#include "BondOrder.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template using nb_array = nb::ndarray; + +namespace wrap { + +void accumulateBondOrder(const std::shared_ptr& self, + const std::shared_ptr& nq, + const nb_array>& orientations, + const nb_array>& query_points, + const nb_array>& query_orientations, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs) +{ + unsigned int const n_query_points = query_points.shape(0); + // std::cout << n_query_points << std::endl; + + // if (query_points.is_none()){ + // auto* query_points_data = nq->getPoints(); + // } + // else { + // auto* query_points_data= reinterpret_cast*>(query_points.data()); + // } + + auto* orientations_data = reinterpret_cast*>(orientations.data()); + auto* query_points_data = reinterpret_cast*>(query_points.data()); + auto* query_orientations_data = reinterpret_cast*>(query_orientations.data()); + + self->accumulate(nq, orientations_data, query_points_data, query_orientations_data, n_query_points, nlist, + qargs); +} + +}; // namespace wrap + +namespace detail { + +void export_BondOrder(nb::module_& module) +{ + nb::enum_(module, "BondOrderMode") + .value("bod", BondOrderMode::bod) + .value("lbod", BondOrderMode::lbod) + .value("obcd", BondOrderMode::obcd) + .value("oocd", BondOrderMode::oocd) + .export_values(); + + nb::class_(module, "BondOrder") + .def(nb::init()) + .def("getBondOrder", &BondOrder::getBondOrder) + .def("getBinCounts", &BondOrder::getBinCounts) + .def("getBinCenters", &BondOrder::getBinCenters) + .def("getBinEdges", &BondOrder::getBinEdges) + .def("getBox", &BondOrder::getBox) + .def("getAxisSizes", &BondOrder::getAxisSizes) + .def("getMode", &BondOrder::getMode) + .def("accumulate", &wrap::accumulateBondOrder, nanobind::arg("nq").none(), + nanobind::arg("orientations"), nanobind::arg("query_points"), + nanobind::arg("query_orientations"), + // nanobind::arg("n_query_points"), + nanobind::arg("nlist").none(), nanobind::arg("qargs").none()) + .def("reset", &BondOrder::reset); +} + +}; // namespace detail + +}; }; // namespace freud::environment diff --git a/freud/environment/export-LocalBondProjection.cc b/freud/environment/export-LocalBondProjection.cc new file mode 100644 index 000000000..bfb10563c --- /dev/null +++ b/freud/environment/export-LocalBondProjection.cc @@ -0,0 +1,61 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly + +#include "LocalBondProjection.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { + +void compute(std::shared_ptr& local_bond_projection, + const std::shared_ptr& nq, + nb_array>& orientations, + nb_array>& query_points, + nb_array>& projected_vectors, + nb_array>& equiv_orientations, + const std::shared_ptr& nlist, const locality::QueryArgs& qargs) +{ + auto* orientations_data = reinterpret_cast*>(orientations.data()); + auto* query_points_data = reinterpret_cast*>(query_points.data()); + auto* proj_vectors_data = reinterpret_cast*>(projected_vectors.data()); + auto* equiv_orientations_data = reinterpret_cast*>(equiv_orientations.data()); + const unsigned int n_proj_vec = projected_vectors.shape(0); + const unsigned int n_query_points = query_points.shape(0); + const unsigned int n_equiv_orientations = equiv_orientations.shape(0); + local_bond_projection->compute(nq, orientations_data, query_points_data, n_query_points, + proj_vectors_data, n_proj_vec, equiv_orientations_data, + n_equiv_orientations, nlist, qargs); +} + +}; // namespace wrap + +namespace detail { + +void export_LocalBondProjection(nb::module_& module) +{ + nb::class_(module, "LocalBondProjection") + .def(nb::init<>()) + .def("getNList", &LocalBondProjection::getNList) + .def("getProjections", &LocalBondProjection::getProjections) + .def("getNormedProjections", &LocalBondProjection::getNormedProjections) + .def("compute", &wrap::compute, nb::arg("nq"), nb::arg("orientations"), nb::arg("query_points"), + nb::arg("projected_vectors"), nb::arg("equiv_orientations"), nb::arg("nlist").none(), + nb::arg("qargs")); +} + +}; // namespace detail +}; }; // namespace freud::environment diff --git a/freud/environment/export-LocalDescriptors.cc b/freud/environment/export-LocalDescriptors.cc new file mode 100644 index 000000000..7481ad86d --- /dev/null +++ b/freud/environment/export-LocalDescriptors.cc @@ -0,0 +1,62 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly + +#include "LocalDescriptors.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute(const std::shared_ptr& local_descriptors, + const std::shared_ptr& nq, + const nb_array>& query_points, const unsigned int n_query_points, + const nb_array>& orientations, + const std::shared_ptr& nlist, const locality::QueryArgs& qargs, + const unsigned int max_num_neighbors) +{ + auto* query_points_data = reinterpret_cast*>(query_points.data()); + auto* orientations_data = reinterpret_cast*>(orientations.data()); + local_descriptors->compute(nq, query_points_data, n_query_points, orientations_data, nlist, qargs, + max_num_neighbors); +} + +}; // namespace wrap + +namespace detail { + +void export_LocalDescriptors(nb::module_& module) +{ + nb::enum_(module, "LocalDescriptorOrientation") + .value("LocalNeighborhood", LocalDescriptorOrientation::LocalNeighborhood) + .value("Global", LocalDescriptorOrientation::Global) + .value("ParticleLocal", LocalDescriptorOrientation::ParticleLocal) + .export_values(); + + nb::class_(module, "LocalDescriptors") + .def(nb::init()) + .def("getNList", &LocalDescriptors::getNList) + .def("getSph", &LocalDescriptors::getSph) + .def("getNSphs", &LocalDescriptors::getNSphs) + .def("getLMax", &LocalDescriptors::getLMax) + .def("getNegativeM", &LocalDescriptors::getNegativeM) + .def("getMode", &LocalDescriptors::getMode) + .def("compute", &wrap::compute, nb::arg("nq"), nb::arg("query_points"), nb::arg("n_query_points"), + nb::arg("orientations").none(), nb::arg("nlist").none(), nb::arg("qargs"), + nb::arg("max_num_neighbors")); +} + +}; // namespace detail +}; }; // namespace freud::environment diff --git a/freud/environment/export-MatchEnv.cc b/freud/environment/export-MatchEnv.cc new file mode 100644 index 000000000..8bd37f039 --- /dev/null +++ b/freud/environment/export-MatchEnv.cc @@ -0,0 +1,149 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include +#include +#include +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include +#include + +#include "Box.h" +#include "MatchEnv.h" +#include "NeighborList.h" +#include "NeighborQuery.h" +#include "VectorMath.h" + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute_env_motif_match(const std::shared_ptr& env_motif_match, + const std::shared_ptr& nq, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs, + const nb_array>& motif, + const unsigned int motif_size, const float threshold, const bool registration) +{ + auto* motif_data = reinterpret_cast*>(motif.data()); + env_motif_match->compute(nq, nlist, qargs, motif_data, motif_size, threshold, registration); +} + +void compute_env_rmsd_min(const std::shared_ptr& env_rmsd_min, + const std::shared_ptr& nq, + const std::shared_ptr& nlist, + const locality::QueryArgs& qargs, + const nb_array>& motif, const unsigned int motif_size, + const bool registration) +{ + auto* motif_data = reinterpret_cast*>(motif.data()); + env_rmsd_min->compute(nq, nlist, qargs, motif_data, motif_size, registration); +} + +std::pair> +compute_minimize_RMSD(const box::Box& box, const nb_array>& refPoints1, + nb_array>& refPoints2, unsigned int numRef, + float min_rmsd, bool registration) +{ + float min_rmsd_modified = min_rmsd; + auto* refPoints1_data = reinterpret_cast*>(refPoints1.data()); + auto* refPoints2_data = reinterpret_cast*>(refPoints2.data()); + return std::make_pair( + min_rmsd_modified, + minimizeRMSD(box, refPoints1_data, refPoints2_data, numRef, min_rmsd_modified, registration)); +} + +std::map +compute_is_similar(const box::Box& box, const nb_array>& refPoints1, + nb_array>& refPoints2, unsigned int numRef, + float threshold_sq, bool registration) +{ + auto* refPoints1_data = reinterpret_cast*>(refPoints1.data()); + auto* refPoints2_data = reinterpret_cast*>(refPoints2.data()); + return isSimilar(box, refPoints1_data, refPoints2_data, numRef, threshold_sq, registration); +} + +nb::list convertVectorsToLists(const std::vector>>& vecs) +{ + nb::list vecs_python; + for (const auto& vec : vecs) + { + nb::list vec_python; + for (const auto& v : vec) + { + nb::list v_python; + v_python.append(v.x); + v_python.append(v.y); + v_python.append(v.z); + vec_python.append(v_python); + } + vecs_python.append(vec_python); + } + return vecs_python; +} + +nb::object getClusterEnv(const std::shared_ptr& env_cls) +{ + return convertVectorsToLists(env_cls->getClusterEnvironments()); +} + +nb::object getPointEnv(const std::shared_ptr& env_cls) +{ + return convertVectorsToLists(env_cls->getPointEnvironments()); +} + +nb::object getPointEnvmm(const std::shared_ptr& env_cls) +{ + return convertVectorsToLists(env_cls->getPointEnvironments()); +} + +}; // namespace wrap + +namespace detail { + +void export_MatchEnv(nb::module_& module) +{ + module.def("minimizeRMSD", &wrap::compute_minimize_RMSD); + + module.def("isSimilar", &wrap::compute_is_similar); + + nb::class_(module, "MatchEnv") + .def(nb::init<>()) + .def("getPointEnvironments", &MatchEnv::getPointEnvironments); + + nb::class_(module, "EnvironmentCluster") + .def(nb::init<>()) + .def("compute", &EnvironmentCluster::compute, nb::arg("nq"), nb::arg("nlist").none(), + nb::arg("qargs"), nb::arg("env_nlist").none(), nb::arg("env_qargs"), nb::arg("threshold"), + nb::arg("registration")) + .def("getClusterEnvironments", &wrap::getClusterEnv) + .def("getPointEnvironments", &wrap::getPointEnv) + .def("getClusters", &EnvironmentCluster::getClusters) + .def("getNumClusters", &EnvironmentCluster::getNumClusters); + + nb::class_(module, "EnvironmentMotifMatch") + .def(nb::init<>()) + .def("compute", &wrap::compute_env_motif_match, nb::arg("nq"), nb::arg("nlist").none(), + nb::arg("qargs"), nb::arg("motif"), nb::arg("motif_size"), nb::arg("threshold"), + nb::arg("registration")) + .def("getPointEnvironments", &wrap::getPointEnvmm) + .def("getMatches", &EnvironmentMotifMatch::getMatches); + + nb::class_(module, "EnvironmentRMSDMinimizer") + .def(nb::init<>()) + .def("compute", &wrap::compute_env_rmsd_min, nb::arg("nq"), nb::arg("nlist").none(), nb::arg("qargs"), + nb::arg("motif"), nb::arg("motif_size"), nb::arg("registration")) + .def("getRMSDs", &EnvironmentRMSDMinimizer::getRMSDs); +} + +}; // namespace detail +}; }; // namespace freud::environment diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc new file mode 100644 index 000000000..725cbf534 --- /dev/null +++ b/freud/environment/module-environment.cc @@ -0,0 +1,26 @@ +// Copyright (c) 2010-2024 The Regents of the University of Michigan +// This file is from the freud project, released under the BSD 3-Clause License. + +#include +#include + +namespace freud::environment::detail { + +void export_AngularSeparationNeighbor(nanobind::module_& m); +void export_AngularSeparationGlobal(nanobind::module_& m); +void export_LocalBondProjection(nanobind::module_& m); +void export_LocalDescriptors(nanobind::module_& m); +void export_BondOrder(nanobind::module_& m); +void export_MatchEnv(nanobind::module_& m); +}; // namespace freud::environment::detail +using namespace freud::environment::detail; + +NB_MODULE(_environment, module) // NOLINT(misc-use-anonymous-namespace): caused by nanobind +{ + export_AngularSeparationNeighbor(module); + export_AngularSeparationGlobal(module); + export_LocalBondProjection(module); + export_LocalDescriptors(module); + export_BondOrder(module); + export_MatchEnv(module); +} diff --git a/freud/locality.py b/freud/locality.py index 06114ade4..c90c7d64f 100644 --- a/freud/locality.py +++ b/freud/locality.py @@ -992,21 +992,21 @@ def box(self): @_Compute._computed_property def bin_counts(self): - """:class:`numpy.ndarray`: The bin counts in the histogram.""" + """:math:`\\left(N_0, N_1 \\right)` :class:`numpy.ndarray`: The bin counts in the histogram.""" # noqa: E501 return self._cpp_obj.getBinCounts().toNumpyArray() @property def bin_centers(self): - """:class:`numpy.ndarray`: The centers of each bin in the histogram - (has the same shape as the histogram itself).""" + """:class:`list` (:class:`numpy.ndarray`): The centers of each bin in the + histogram (has the same shape in each dimension as the histogram itself).""" centers = self._cpp_obj.getBinCenters() return [np.array(c) for c in centers] @property def bin_edges(self): - """:class:`numpy.ndarray`: The edges of each bin in the histogram (is - one element larger in each dimension than the histogram because each - bin has a lower and upper bound).""" + """:class:`list` (:class:`numpy.ndarray`): The edges of each bin in the + histogram (is one element larger in each dimension than the histogram + because each bin has a lower and upper bound).""" edges = self._cpp_obj.getBinEdges() return [np.array(e) for e in edges] diff --git a/freud/util/BiMap.h b/freud/util/BiMap.h index 003e2fbf0..9e1373bcc 100644 --- a/freud/util/BiMap.h +++ b/freud/util/BiMap.h @@ -170,19 +170,19 @@ template class BiMap friend BiMap; private: - BiMap& b() const + BiMap& b() { BiMap t; return *reinterpret_cast(reinterpret_cast(this) - (((size_t) (&(&t)->left) - ((size_t) &t)))); } - Pair* getPairPtr(const T* Item_in) const + Pair* getPairPtr(const T* Item_in) { return b().getPairVal(Item_in); } - U& getVal(const T* Item_in) const + U& getVal(const T* Item_in) { return getPairPtr(Item_in)->second; } @@ -198,7 +198,7 @@ template class BiMap return getVal(*itr); } - U& operator[](const T& Key_in) const + U& operator[](const T& Key_in) { if (!this->has(Key_in)) { @@ -209,17 +209,17 @@ template class BiMap return getVal(*(this->b().set_A.find(&Key_in))); } - auto find(const T& Key_in) const -> decltype(this->b().set_A.find(&Key_in)) + auto find(const T& Key_in) -> decltype(this->b().set_A.find(&Key_in)) { return this->b().set_A.find(&Key_in); } - int count(const T& Key_in) const + int count(const T& Key_in) { return this->b().set_A.count(&Key_in); } - bool has(const T& Key_in) const + bool has(const T& Key_in) { return find(Key_in) != std::end(this->b().set_A); } @@ -243,7 +243,7 @@ template class BiMap friend BiMap; private: - BiMap& b() const + BiMap& b() { BiMap t; return *reinterpret_cast(reinterpret_cast(this) @@ -255,18 +255,18 @@ template class BiMap return reinterpret_cast(Item_in) - offsetof(Pair, second); } - Pair* getPairPtr(const U* Item_in) const + Pair* getPairPtr(const U* Item_in) { return b().getPairVal(getPairPtr_B(Item_in)); } - T& getVal(const U* Item_in) const + T& getVal(const U* Item_in) { return getPairPtr(Item_in)->first; } public: - const T& at(const U& Key_in) const + const T& at(const U& Key_in) { const auto& itr(this->b().set_B.find(&Key_in)); if (itr == std::end(this->b().set_B)) @@ -287,17 +287,17 @@ template class BiMap return getVal(*(this->b().set_B.find(&Key_in))); } - auto find(const U& Key_in) const -> decltype(this->b().set_B.find(&Key_in)) + auto find(const U& Key_in) -> decltype(this->b().set_B.find(&Key_in)) { return this->b().set_B.find(&Key_in); } - int count(const U& Key_in) const + int count(const U& Key_in) { return this->b().set_B.count(&Key_in); } - bool has(const U& Key_in) const + bool has(const U& Key_in) { return find(Key_in) != std::end(this->b().set_B); } diff --git a/freud/util/export-ManagedArray.cc b/freud/util/export-ManagedArray.cc index a8b8645a0..3c85991d7 100644 --- a/freud/util/export-ManagedArray.cc +++ b/freud/util/export-ManagedArray.cc @@ -13,9 +13,10 @@ void export_ManagedArray(nanobind::module_& module) { // python wrapper classes for ManagedArray export_ManagedArray(module, "ManagedArray_float"); export_ManagedArray(module, "ManagedArray_double"); + export_ManagedArray>(module, "ManagedArray_complexfloat"); export_ManagedArray(module, "ManagedArray_unsignedint"); export_ManagedArray>(module, "ManagedArrayVec3_float"); - export_ManagedArray>(module, "ManagedArray_complexfloat"); + export_ManagedArray(module, "ManagedArray_char"); }; }; // namespace freud::util::detail diff --git a/tests/test_environment_bond_order.py b/tests/test_environment_bond_order.py index 09df84b9b..c3a787c62 100644 --- a/tests/test_environment_bond_order.py +++ b/tests/test_environment_bond_order.py @@ -40,6 +40,9 @@ def test_bond_order(self): # Test access bo.box bo.bond_order + bo.bin_counts + bo.bin_edges + bo.bin_centers # Test all the basic attributes. assert bo.nbins[0] == n_bins_theta @@ -57,8 +60,7 @@ def test_bond_order(self): box, positions, positions, "nearest", r_max, num_neighbors, True ) for nq, neighbors in test_set: - # Test that lbod gives identical results when orientations are the - # same. + # Test that lbod gives identical results when orientations are the same. # TODO: Find a way to test a rotated system to ensure that lbod gives # the desired results. bo = freud.environment.BondOrder(nbins, mode="lbod") diff --git a/tests/test_environment_local_descriptors.py b/tests/test_environment_local_descriptors.py index 2ac0da5e7..d63c55eea 100644 --- a/tests/test_environment_local_descriptors.py +++ b/tests/test_environment_local_descriptors.py @@ -347,7 +347,7 @@ def test_ld(self): spherical harmonics manually and verifying them.""" sph_harm = pytest.importorskip("scipy.special").sph_harm - atol = 1e-5 + atol = 1e-4 L = 8 N = 100 box, points = freud.data.make_random_system(L, N) diff --git a/tests/test_environment_match_env.py b/tests/test_environment_match_env.py index 3ca04e968..5bbd273d6 100644 --- a/tests/test_environment_match_env.py +++ b/tests/test_environment_match_env.py @@ -483,6 +483,8 @@ def test_square(self): ) matches = match.matches + assert matches.dtype == bool + for i in range(len(motif)): assert not matches[i] assert matches[len(motif)]