From e499f2ad0c8601bd5b5385a9b2dc11b98f5534d2 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 11 Sep 2024 11:45:17 -0400 Subject: [PATCH 01/28] Begin porting environment to nanobind --- freud/CMakeLists.txt | 15 +++- freud/environment/AngularSeparation.cc | 21 ++--- freud/environment/AngularSeparation.h | 18 ++--- .../export-AngularSeparationGlobal.cc | 81 +++++++++++++++++++ .../export-AngularSeparationNeighbor.cc | 31 +++++++ 5 files changed, 145 insertions(+), 21 deletions(-) create mode 100644 freud/environment/export-AngularSeparationGlobal.cc create mode 100644 freud/environment/export-AngularSeparationNeighbor.cc diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index b84bbb65c..54839ba8f 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -43,8 +43,8 @@ add_library( # diffraction/StaticStructureFactorDirect.h # diffraction/StaticStructureFactorDirect.cc # environment - # environment/AngularSeparation.h - # environment/AngularSeparation.cc + environment/AngularSeparation.h + environment/AngularSeparation.cc # environment/BondOrder.h # environment/BondOrder.cc # environment/LocalBondProjection.h @@ -154,6 +154,17 @@ 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 +) +target_link_libraries(_environment PUBLIC freud TBB:tbb) +target_set_install_rpath(_environment) + + # locality nanobind_add_module( _locality diff --git a/freud/environment/AngularSeparation.cc b/freud/environment/AngularSeparation.cc index 6646491be..13a2ab0c0 100644 --- a/freud/environment/AngularSeparation.cc +++ b/freud/environment/AngularSeparation.cc @@ -49,32 +49,33 @@ float computeMinSeparationAngle(const quat& ref_q, const quat& q, return min_angle; } -void AngularSeparationNeighbor::compute(const locality::NeighborQuery* nq, const quat* orientations, +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, 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]; - 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)); + const size_t j((*m_nlist->getNeighbors())(bond, 1)); quat query_q = query_orientations[j]; float theta = computeMinSeparationAngle(q, query_q, equiv_orientations, n_equiv_orientations); - m_angles[bond] = theta; + (*m_angles)[bond] = theta; } } }); @@ -85,7 +86,7 @@ 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) @@ -96,7 +97,7 @@ void AngularSeparationGlobal::compute(const quat* global_orientations, un quat global_q = global_orientations[j]; 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..b7909d0bb 100644 --- a/freud/environment/AngularSeparation.h +++ b/freud/environment/AngularSeparation.h @@ -36,13 +36,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 +60,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, + unsigned int n_equiv_orientations, const std::shared_ptr &nlist, 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/export-AngularSeparationGlobal.cc b/freud/environment/export-AngularSeparationGlobal.cc new file mode 100644 index 000000000..a39c08580 --- /dev/null +++ b/freud/environment/export-AngularSeparationGlobal.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 // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include + +#include "BondHistogramCompute.h" + +namespace nb = nanobind; + +namespace freud { namespace locality { + +namespace wrap { + +/** + * Here we convert a vector of vectors into a list of lists for returning to python. + * */ +template +inline nb::object vectorVectorsToListLists(const std::vector>& vectorOfVectors) +{ + nb::list outer_python_list; + for (const auto& vector : vectorOfVectors) + { + nb::list inner_python_list; + for (const auto& element : vector) + { + inner_python_list.append(element); + } + outer_python_list.append(inner_python_list); + } + return outer_python_list; +} + +nb::object getBinCenters(const std::shared_ptr& bondHist) +{ + auto bin_centers_cpp = bondHist->getBinCenters(); + return vectorVectorsToListLists(bin_centers_cpp); +} + +nb::object getBinEdges(const std::shared_ptr& bondHist) +{ + auto bin_edges_cpp = bondHist->getBinEdges(); + return vectorVectorsToListLists(bin_edges_cpp); +} + +nb::object getBounds(const std::shared_ptr& bondHist) +{ + auto bounds_cpp = bondHist->getBounds(); + + // convert the vector of pairs to a list of tuples + nb::list python_list; + for (const auto& pair : bounds_cpp) + { + nb::tuple python_tuple = nb::make_tuple(pair.first, pair.second); + python_list.append(python_tuple); + } + return python_list; +} + +}; // namespace wrap + +namespace detail { + +void export_BondHistogramCompute(nb::module_& module) +{ + nb::class_(module, "BondHistogramCompute") + .def("getBox", &BondHistogramCompute::getBox) + .def("reset", &BondHistogramCompute::reset) + .def("getBinCounts", &BondHistogramCompute::getBinCounts) + .def("getAxisSizes", &BondHistogramCompute::getAxisSizes) + .def("getBinCenters", &wrap::getBinCenters) + .def("getBinEdges", &wrap::getBinEdges) + .def("getBounds", &wrap::getBounds); +} + +}; // namespace detail + +}; }; // namespace freud::locality diff --git a/freud/environment/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc new file mode 100644 index 000000000..028cd8962 --- /dev/null +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -0,0 +1,31 @@ +// 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 // NOLINT(misc-include-cleaner): used implicitly +#include // NOLINT(misc-include-cleaner): used implicitly +#include + +#include "AngularSeparation.h" + +namespace nb = nanobind { + +namespace freud { namespace environment { + +namespace wrap { + + +namespace detail { + +void export_AngularSeparationNeighbor(nb::module_& module) +{ + nb::class_(module, "AngularSeparationNeighbor") + .def("getNList", &AngularSeparationNeighbor::getNList) + .def("getAngles", &AngularSeparationNeighbor::getAngles) + .def("compute", &AngularSeparationNeighbor::compute); +} + +}; // namespace detail + +}; }; }; }; // namespace freud::locality From 61827eed620d3b2f7a239db671cfe558f07898ca Mon Sep 17 00:00:00 2001 From: Alain Kadar Date: Wed, 11 Sep 2024 13:12:31 -0400 Subject: [PATCH 02/28] Export classes --- freud/CMakeLists.txt | 11 +-- .../export-AngularSeparationGlobal.cc | 73 +++---------------- .../export-AngularSeparationNeighbor.cc | 5 +- freud/environment/module-environment.cc | 18 +++++ 4 files changed, 35 insertions(+), 72 deletions(-) create mode 100644 freud/environment/module-environment.cc diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 54839ba8f..38c3fea2d 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -145,10 +145,6 @@ target_include_directories(freud SYSTEM PUBLIC ${PROJECT_SOURCE_DIR}/extern/) # 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) @@ -159,9 +155,8 @@ nanobind_add_module( _environment environment/module-environment.cc environment/export-AngularSeparationNeighbor.cc - #environment/export-AngularSeparationGlobal.cc -) -target_link_libraries(_environment PUBLIC freud TBB:tbb) + environment/export-AngularSeparationGlobal.cc) +target_link_libraries(_environment PUBLIC freud TBB::tbb) target_set_install_rpath(_environment) @@ -222,7 +217,7 @@ if(SKBUILD) install(TARGETS _box DESTINATION freud) # install(TARGETS _cluster DESTINATION freud) install(TARGETS _density # DESTINATION freud) install(TARGETS _diffraction DESTINATION freud) - # install(TARGETS _environment 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/environment/export-AngularSeparationGlobal.cc b/freud/environment/export-AngularSeparationGlobal.cc index a39c08580..5e2dde451 100644 --- a/freud/environment/export-AngularSeparationGlobal.cc +++ b/freud/environment/export-AngularSeparationGlobal.cc @@ -7,75 +7,24 @@ #include // NOLINT(misc-include-cleaner): used implicitly #include -#include "BondHistogramCompute.h" +#include "AngularSeparation.h" -namespace nb = nanobind; - -namespace freud { namespace locality { - -namespace wrap { - -/** - * Here we convert a vector of vectors into a list of lists for returning to python. - * */ -template -inline nb::object vectorVectorsToListLists(const std::vector>& vectorOfVectors) -{ - nb::list outer_python_list; - for (const auto& vector : vectorOfVectors) - { - nb::list inner_python_list; - for (const auto& element : vector) - { - inner_python_list.append(element); - } - outer_python_list.append(inner_python_list); - } - return outer_python_list; -} - -nb::object getBinCenters(const std::shared_ptr& bondHist) -{ - auto bin_centers_cpp = bondHist->getBinCenters(); - return vectorVectorsToListLists(bin_centers_cpp); -} - -nb::object getBinEdges(const std::shared_ptr& bondHist) -{ - auto bin_edges_cpp = bondHist->getBinEdges(); - return vectorVectorsToListLists(bin_edges_cpp); -} -nb::object getBounds(const std::shared_ptr& bondHist) -{ - auto bounds_cpp = bondHist->getBounds(); +namespace nb = nanobind; - // convert the vector of pairs to a list of tuples - nb::list python_list; - for (const auto& pair : bounds_cpp) - { - nb::tuple python_tuple = nb::make_tuple(pair.first, pair.second); - python_list.append(python_tuple); - } - return python_list; -} +namespace freud { namespace environment { -}; // namespace wrap +namespace wrap {}; namespace detail { -void export_BondHistogramCompute(nb::module_& module) +void export_AngularSeparationGlobal(nb::module_& module) { - nb::class_(module, "BondHistogramCompute") - .def("getBox", &BondHistogramCompute::getBox) - .def("reset", &BondHistogramCompute::reset) - .def("getBinCounts", &BondHistogramCompute::getBinCounts) - .def("getAxisSizes", &BondHistogramCompute::getAxisSizes) - .def("getBinCenters", &wrap::getBinCenters) - .def("getBinEdges", &wrap::getBinEdges) - .def("getBounds", &wrap::getBounds); + nb::class_(module, "AngularSeparationGlobal") + .def(nb::init<>()) + .def("getAngles", &AngularSeparationGlobal::getAngles) + .def("compute", &AngularSeparationGlobal::compute); } -}; // namespace detail - -}; }; // namespace freud::locality +}; }; // namespace detail +}; // namespace freud::locality diff --git a/freud/environment/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc index 028cd8962..f2cb4f25f 100644 --- a/freud/environment/export-AngularSeparationNeighbor.cc +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -9,12 +9,13 @@ #include "AngularSeparation.h" -namespace nb = nanobind { +namespace nb = nanobind; namespace freud { namespace environment { namespace wrap { +} namespace detail { @@ -28,4 +29,4 @@ void export_AngularSeparationNeighbor(nb::module_& module) }; // namespace detail -}; }; }; }; // namespace freud::locality +}; }; // namespace freud::locality diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc new file mode 100644 index 000000000..2e3e11b2c --- /dev/null +++ b/freud/environment/module-environment.cc @@ -0,0 +1,18 @@ +// 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); +}; // 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); +} From babe6581f22de88a6234a1eef2c8270c0ee92239 Mon Sep 17 00:00:00 2001 From: Tim Moore Date: Wed, 11 Sep 2024 15:29:59 -0400 Subject: [PATCH 03/28] Tests are passing for angular separation classes --- freud/__init__.py | 4 +- freud/environment.pyx | 80 ++++--------------- .../export-AngularSeparationGlobal.cc | 22 ++++- .../export-AngularSeparationNeighbor.cc | 23 ++++++ 4 files changed, 62 insertions(+), 67 deletions(-) diff --git a/freud/__init__.py b/freud/__init__.py index 4e84eb160..3aeab0065 100644 --- a/freud/__init__.py +++ b/freud/__init__.py @@ -3,7 +3,7 @@ # cluster,; density,; diffraction,; environment,; interface,; msd,; order, -from . import box, data, locality, parallel, pmft +from . import box, data, environment, locality, parallel, pmft from .box import Box from .locality import AABBQuery, LinkCell, NeighborList from .parallel import NumThreads, get_num_threads, set_num_threads @@ -21,7 +21,7 @@ "data", # "density", # "diffraction", - # "environment", + "environment", # "interface", "locality", # "msd", diff --git a/freud/environment.pyx b/freud/environment.pyx index d7a67fb10..31883bd65 100644 --- a/freud/environment.pyx +++ b/freud/environment.pyx @@ -882,19 +882,11 @@ cdef class _EnvironmentRMSDMinimizer(_MatchEnv): freud.util.arr_type_t.FLOAT) -cdef class AngularSeparationNeighbor(_PairCompute): +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 + self._cpp_obj = freud._environment.AngularSeparationNeighbor() def compute(self, system, orientations, query_points=None, query_orientations=None, @@ -935,14 +927,8 @@ cdef class AngularSeparationNeighbor(_PairCompute): `_ (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 = \ + nq, nlist, qargs, query_points, num_query_points = \ self._preprocess_arguments(system, query_points, neighbors) orientations = freud.util._convert_array( @@ -956,22 +942,16 @@ cdef class AngularSeparationNeighbor(_PairCompute): 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 - 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)) + self._cpp_obj.compute( + nq._cpp_obj, + l_orientations, + query_points, + query_orientations, + equiv_orientations, + nlist._cpp_obj, + qargs._cpp_obj) return self @_Compute._computed_property @@ -979,9 +959,7 @@ cdef class AngularSeparationNeighbor(_PairCompute): """: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) + return self._cpp_obj.getAngles().toNumpyArray() def __repr__(self): return "freud.environment.{cls}()".format( @@ -991,24 +969,16 @@ cdef class AngularSeparationNeighbor(_PairCompute): 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 -cdef class AngularSeparationGlobal(_Compute): +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 + self._cpp_obj = freud._environment.AngularSeparationGlobal() def compute(self, global_orientations, orientations, equiv_orientations=np.array([[1, 0, 0, 0]])): @@ -1039,30 +1009,14 @@ cdef class AngularSeparationGlobal(_Compute): 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) + 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 freud.util.make_managed_numpy_array( - &self.thisptr.getAngles(), - freud.util.arr_type_t.FLOAT) + return self._cpp_obj.getAngles().toNumpyArray() def __repr__(self): return "freud.environment.{cls}()".format( diff --git a/freud/environment/export-AngularSeparationGlobal.cc b/freud/environment/export-AngularSeparationGlobal.cc index 5e2dde451..43a049a49 100644 --- a/freud/environment/export-AngularSeparationGlobal.cc +++ b/freud/environment/export-AngularSeparationGlobal.cc @@ -3,9 +3,9 @@ #include #include +#include #include // NOLINT(misc-include-cleaner): used implicitly #include // NOLINT(misc-include-cleaner): used implicitly -#include #include "AngularSeparation.h" @@ -14,7 +14,25 @@ namespace nb = nanobind; namespace freud { namespace environment { -namespace wrap {}; +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute(const std::shared_ptr &self, + 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()); + self->compute(global_orientations_data, n_global, orientations_data, n_points, equiv_orientations_data, n_equiv_orientations); + } + +}; namespace detail { diff --git a/freud/environment/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc index f2cb4f25f..e911f6108 100644 --- a/freud/environment/export-AngularSeparationNeighbor.cc +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -3,6 +3,7 @@ #include #include +#include #include // NOLINT(misc-include-cleaner): used implicitly #include // NOLINT(misc-include-cleaner): used implicitly #include @@ -13,8 +14,30 @@ namespace nb = nanobind; namespace freud { namespace environment { +template +using nb_array = nanobind::ndarray; + namespace wrap { +void compute( + std::shared_ptr &self, + std::shared_ptr &nq, + nb_array> & orientations, + nb_array> & query_points, + nb_array> & query_orientations, + nb_array> & equiv_orientations, + 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()); + unsigned int n_equiv_orientations = equiv_orientations.shape(0); + unsigned int n_query_points = query_orientations.shape(0); + self->compute(nq, orientations_data, query_points_data, query_orientations_data, n_query_points, equiv_orientations_data, n_equiv_orientations, nlist, qargs); + } } namespace detail { From 8a47eae20cf19dec1bcde60385b4d1beaf208274 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 16:32:11 -0400 Subject: [PATCH 04/28] move environment to pure python --- freud/environment.py | 1155 +++++++++++++++++++++++++++++++++++++++++ freud/environment.pyx | 1144 ---------------------------------------- 2 files changed, 1155 insertions(+), 1144 deletions(-) create mode 100644 freud/environment.py delete mode 100644 freud/environment.pyx diff --git a/freud/environment.py b/freud/environment.py new file mode 100644 index 000000000..fc2c74bfb --- /dev/null +++ b/freud/environment.py @@ -0,0 +1,1155 @@ +# 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. + +r""" +The :class:`freud.environment` module contains functions which characterize the +local environments of particles in the system. These methods use the positions +and orientations of particles in the local neighborhood of a given particle to +characterize the particle environment. +""" + +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 + + +class AngularSeparationNeighbor(_PairCompute): + r"""Calculates the minimum angles of separation between orientations and + query orientations.""" + + def __init__(self): + self._cpp_obj = freud._environment.AngularSeparationNeighbor() + + 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 + + nq, nlist, qargs, 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) + ) + + self._cpp_obj.compute( + nq._cpp_obj, + orientations, + query_points, + query_orientations, + equiv_orientations, + nlist._cpp_obj, + qargs._cpp_obj, + ) + return self + + @_PairCompute._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__}()" + + @_PairCompute._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 + + +class AngularSeparationGlobal(_Compute): + r"""Calculates the minimum angles of separation between orientations and + global orientations.""" + + def __init__(self): + self._cpp_obj = freud._environment.AngularSeparationGlobal() + + 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) + ) + + 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. +# +# The bond orientational order diagram (BOOD) is a way of studying the +# average local environments experienced by particles. In a BOOD, a particle +# and its nearest neighbors (determined by either a prespecified number of +# neighbors or simply a cutoff distance) are treated as connected by a bond +# joining their centers. All of the bonds in the system are then binned by +# their azimuthal (:math:`\theta`) and polar (:math:`\phi`) angles to +# indicate the location of a particle's neighbors relative to itself. The +# distance between the particle and its neighbor is only important when +# determining whether it is counted as a neighbor, but is not part of the +# BOOD; as such, the BOOD can be viewed as a projection of all bonds onto the +# unit sphere. The resulting 2D histogram provides insight into how particles +# are situated relative to one-another in a system. +# +# This class provides access to the classical BOOD as well as a few useful +# variants. These variants can be accessed via the :code:`mode` arguments to +# the :meth:`~BondOrder.compute` method. Available modes of calculation are: +# +# * :code:`'bod'` (Bond Order Diagram, *default*): +# This mode constructs the default BOOD, which is the 2D histogram +# containing the number of bonds formed through each azimuthal +# :math:`\left( \theta \right)` and polar :math:`\left( \phi \right)` +# angle. +# +# * :code:`'lbod'` (Local Bond Order Diagram): +# In this mode, a particle's neighbors are rotated into the local frame of +# the particle before the BOOD is calculated, *i.e.* the directions of +# bonds are determined relative to the orientation of the particle rather +# than relative to the global reference frame. An example of when this mode +# would be useful is when a system is composed of multiple grains of the +# same crystal; the normal BOOD would show twice as many peaks as expected, +# but using this mode, the bonds would be superimposed. +# +# * :code:`'obcd'` (Orientation Bond Correlation Diagram): +# This mode aims to quantify the degree of orientational as well as +# translational ordering. As a first step, the rotation that would align a +# particle's neighbor with the particle is calculated. Then, the neighbor +# is rotated **around the central particle** by that amount, which actually +# changes the direction of the bond. One example of how this mode could be +# useful is in identifying plastic crystals, which exhibit translational +# but not orientational ordering. Normally, the BOOD for a plastic crystal +# would exhibit clear structure since there is translational order, but +# with this mode, the neighbor positions would actually be modified, +# resulting in an isotropic (disordered) BOOD. +# +# * :code:`'oocd'` (Orientation Orientation Correlation Diagram): +# This mode is substantially different from the other modes. Rather than +# compute the histogram of neighbor bonds, this mode instead computes a +# histogram of the directors of neighboring particles, where the director +# is defined as the basis vector :math:`\hat{z}` rotated by the neighbor's +# quaternion. The directors are then rotated into the central particle's +# reference frame. This mode provides insight into the local orientational +# environment of particles, indicating, on average, how a particle's +# neighbors are oriented. +# +# Args: +# bins (unsigned int or sequence of length 2): +# If an unsigned int, the number of bins in :math:`\theta` and +# :math:`\phi`. If a sequence of two integers, interpreted as +# :code:`(num_bins_theta, num_bins_phi)`. +# mode (str, optional): +# 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} +# +# def __cinit__(self, bins, str mode="bod"): +# try: +# n_bins_theta, n_bins_phi = bins +# except TypeError: +# 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) +# +# def __dealloc__(self): +# del self.thisptr +# +# @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): +# r"""Calculates the correlation function and adds to the current +# histogram. +# +# 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. Uses identity quaternions if :code:`None` +# (Default value = :code:`None`). +# 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`). +# 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). +# reset (bool): +# 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) +# if orientations is None: +# orientations = np.array([[1, 0, 0, 0]] * nq.points.shape[0]) +# if query_orientations is None: +# query_orientations = orientations +# +# orientations = freud.util._convert_array( +# 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)) +# 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) +# +# @_Compute._computed_property +# def box(self): +# """:class:`freud.box.Box`: Box used in the calculation.""" +# return freud.box.BoxFromCPP(self.thisptr.getBox()) +# +# def __repr__(self): +# 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)) +# +# @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 +# +# +# cdef class LocalDescriptors(_PairCompute): +# r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's +# local environment. +# +# The resulting spherical harmonic array will be a complex-valued +# array of shape :code:`(num_bonds, num_sphs)`. Spherical harmonic +# calculation can be restricted to some number of nearest neighbors +# through the :code:`max_num_neighbors` argument; if a particle has more +# bonds than this number, the last one or more rows of bond spherical +# harmonics for each particle will not be set. This feature is useful for +# computing descriptors on the same system but with different subsets of +# neighbors; a :class:`freud.locality.NeighborList` with the correct +# ordering can then be reused in multiple calls to :meth:`~.compute` +# with different values of :code:`max_num_neighbors` to compute descriptors +# for different local neighborhoods with maximum efficiency. +# +# Args: +# l_max (unsigned int): +# Maximum spherical harmonic :math:`l` to consider. +# negative_m (bool, optional): +# True if we should also calculate :math:`Y_{lm}` for negative +# :math:`m`. (Default value = :code:`True`) +# mode (str, optional): +# Orientation mode to use for environments, either +# :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 (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} +# +# def __cinit__(self, l_max, negative_m=True, mode='neighborhood'): +# cdef freud._environment.LocalDescriptorOrientation l_mode +# 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): +# r"""Calculates the local descriptors of bonds from a set of source +# points to a set of destination points. +# +# Args: +# system: +# Any object that is a valid argument to +# :class:`freud.locality.NeighborQuery.from_system`. +# 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`). +# orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): +# Orientations associated with system points that are used to +# calculate bonds. +# 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). +# max_num_neighbors (unsigned int, optional): +# Hard limit on the maximum number of neighbors to use for each +# 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) +# +# # 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 orientations is None: +# raise RuntimeError( +# ('Orientations must be given to orient LocalDescriptors ' +# 'with particles\' orientations')) +# +# orientations = freud.util._convert_array( +# orientations, shape=(nq.points.shape[0], 4)) +# +# l_orientations = orientations +# l_orientations_ptr = &l_orientations[0, 0] +# +# 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) +# 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._compute = self +# return nlist +# +# @_Compute._computed_property +# 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) +# +# @_Compute._computed_property +# def num_sphs(self): +# """unsigned int: The last number of spherical harmonics computed. This +# is equal to the number of bonds in the last computation, which is at +# most the number of :code:`points` multiplied by the lower of the +# :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() +# +# @property +# def l_max(self): +# """unsigned int: The maximum spherical harmonic :math:`l` calculated +# for.""" +# return self.thisptr.getLMax() +# +# @property +# def negative_m(self): +# """bool: True if we also calculated :math:`Y_{lm}` for negative +# :math:`m`.""" +# return self.thisptr.getNegativeM() +# +# @property +# def mode(self): +# """str: Orientation mode to use for environments, either +# :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() +# for key, value in self.known_modes.items(): +# if value == mode: +# return key +# +# 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) +# +# +# def _minimize_RMSD(box, ref_points, points, registration=False): +# r"""Get the somewhat-optimal RMSD between the set of vectors ref_points +# and the set of vectors points. +# +# Args: +# ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): +# Vectors that make up motif 1. +# points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): +# Vectors that make up motif 2. +# registration (bool, optional): +# If true, first use brute force registration to orient one set +# of environment vectors with respect to the other set such that +# it minimizes the RMSD between the two sets +# (Default value = :code:`False`). +# +# Returns: +# tuple (float, (:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): +# A triplet that gives the associated min_rmsd, rotated (or not) +# 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) +# +# 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] +# +# if nRef1 != nRef2: +# raise ValueError( +# ("The number of vectors in ref_points must MATCH" +# "the number of vectors in points")) +# +# 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] +# +# +# def _is_similar_motif(box, ref_points, points, threshold, registration=False): +# r"""Test if the motif provided by ref_points is similar to the motif +# provided by points. +# +# Args: +# ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): +# Vectors that make up motif 1. +# points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): +# Vectors that make up motif 2. +# threshold (float): +# Maximum magnitude of the vector difference between two vectors, +# below which they are "matching". Typically, a good choice is +# between 10% and 30% of the first well in the radial +# distribution function (this has distance units). +# registration (bool, optional): +# If True, first use brute force registration to orient one set +# of environment vectors with respect to the other set such that +# it minimizes the RMSD between the two sets +# (Default value = :code:`False`). +# +# Returns: +# tuple ((:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): +# A doublet that gives the rotated (or not) set of +# :code:`points`, and the mapping between the vectors of +# :code:`ref_points` and :code:`points` that will make them +# 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) +# +# 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 +# +# if nRef1 != nRef2: +# raise ValueError( +# ("The number of vectors in ref_points must match" +# "the number of vectors in points")) +# +# 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] +# +# +# cdef class _MatchEnv(_PairCompute): +# r"""Parent for environment matching methods. """ +# cdef freud._environment.MatchEnv * matchptr +# +# def __cinit__(self, *args, **kwargs): +# # Abstract class +# pass +# +# @_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] +# +# def __repr__(self): +# return ("freud.environment.{cls}()").format( +# cls=type(self).__name__) +# +# +# cdef 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): +# r"""Determine clusters of particles with matching environments. +# +# An environment is defined by the bond vectors between a particle and its +# neighbors as defined by :code:`env_neighbors`. +# For example, :code:`env_neighbors= {'num_neighbors': 8}` means that every +# particle's local environment is defined by its 8 nearest neighbors. +# Then, each particle's environment is compared to the environments of +# particles that satisfy a different cutoff parameter :code:`cluster_neighbors`. +# For example, :code:`cluster_neighbors={'r_max': 3.0}` +# means that the environment of each particle will be compared to the +# environment of every particle within a distance of 3.0. +# +# Two environments are compared using `point-set registration +# `_. +# +# The thresholding criterion we apply in order to determine if the two point sets +# match is quite conservative: the two point sets match if and only if, +# for every pair of matched points between the sets, the distance between +# the matched pair is less than :code:`threshold`. +# +# When :code:`registration=False`, environments are not rotated prior to comparison +# between them. Which pairs of vectors are compared depends on the order in +# which the vectors of the environment are traversed. +# +# When :code:`registration=True`, we use the `Kabsch algorithm +# `_ to find the optimal +# rotation matrix mapping one environment onto another. The Kabsch +# algorithm assumes a one-to-one correspondence between points in the two +# environments. However, we typically do not know which point in one environment +# corresponds to a particular point in the other environment. The brute force +# solution is to try all possible correspondences, but the resulting +# combinatorial explosion makes this problem infeasible, so we use the thresholding +# criterion above to terminate the search when we have found a permutation +# of points that results in a sufficiently low RMSD. +# +# .. note:: +# +# Using a distance cutoff for :code:`env_neighbors` could +# lead to situations where the :code:`cluster_environments` +# contain different numbers of neighbors. In this case, the +# environments which have a number of neighbors less than +# the environment with the maximum number of neighbors +# :math:`k_{max}` will have their entry in :code:`cluster_environments` +# padded with zero vectors. For example, a cluster environment +# with :math:`m < k` neighbors, will have :math:`k - m` zero +# vectors at the end of its entry in :code:`cluster_environments`. +# +# +# .. warning:: +# +# All vectors of :code:`cluster_environments` and :code:`point_environments` +# are defined with respect to the query particle. +# Zero vectors are only used to pad the cluster vectors so that they +# have the same shape. +# In a future version of freud, zero-padding will be removed. +# +# .. warning:: +# +# Comparisons between two environments are only made when both +# environments contain the same number of neighbors. +# However, no warning will be given at runtime if mismatched +# environments are provided for comparison. +# +# Example:: +# +# >>> import freud +# >>> # Compute clusters of particles with matching environments +# >>> box, points = freud.data.make_random_system(10, 100, seed=0) +# >>> env_cluster = freud.environment.EnvironmentCluster() +# >>> env_cluster.compute( +# ... system = (box, points), +# ... threshold=0.2, +# ... cluster_neighbors={'num_neighbors': 6}, +# ... registration=False) +# freud.environment.EnvironmentCluster() +# +# Args: +# system: +# Any object that is a valid argument to +# :class:`freud.locality.NeighborQuery.from_system`. +# threshold (float): +# Maximum magnitude of the vector difference between two vectors, +# below which they are "matching". Typically, a good choice is +# between 10% and 30% of the first well in the radial +# distribution function (this has distance units). +# cluster_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 +# `_. +# Defines the neighbors used for comparing different particle +# environments (Default value: None). +# env_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 +# `_. +# Defines the neighbors used as the environment of each particle. +# If ``None``, the value provided for ``neighbors`` will be used +# (Default value: None). +# registration (bool, optional): +# If True, first use brute force registration to orient one set +# of environment vectors with respect to the other set such that +# it minimizes the RMSD between the two sets. Enabling this +# 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) +# +# 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) +# 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) +# +# @_Compute._computed_property +# def num_clusters(self): +# """unsigned int: The number of clusters.""" +# return self.thisptr.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] +# +# def plot(self, ax=None): +# """Plot cluster distribution. +# +# Args: +# ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If +# :code:`None`, make a new figure and axis. +# (Default value = :code:`None`) +# +# Returns: +# (: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) +# +# 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): +# r"""Find matches between local arrangements of a set of points and +# a provided motif, as done in :cite:`Teich2019`. +# +# A particle's environment can only match the motif if it contains the +# same number of neighbors as the motif. Any environment with a +# 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 +# +# def compute(self, system, motif, threshold, env_neighbors=None, +# registration=False): +# r"""Determine which particles have local environments +# matching the given environment motif. +# +# .. warning:: +# +# Comparisons between two environments are only made when both +# environments contain the same number of neighbors. +# However, no warning will be given at runtime if mismatched +# environments are provided for comparison. +# +# Args: +# system: +# Any object that is a valid argument to +# :class:`freud.locality.NeighborQuery.from_system`. +# motif ((:math:`N_{points}`, 3) :class:`numpy.ndarray`): +# Vectors that make up the motif against which we are matching. +# threshold (float): +# Maximum magnitude of the vector difference between two vectors, +# below which they are "matching". Typically, a good choice is +# between 10% and 30% of the first well in the radial +# distribution function (this has distance units). +# env_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 +# `_. +# Defines the environment of the query particles +# (Default value: None). +# registration (bool, optional): +# If True, first use brute force registration to orient one set +# of environment vectors with respect to the other set such that +# 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) +# +# 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 +# ) +# +# 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) +# 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) +# +# +# cdef class _EnvironmentRMSDMinimizer(_MatchEnv): +# r"""Find linear transformations that map the environments of points onto a +# motif. +# +# In general, it is recommended to specify a number of neighbors rather than +# just a distance cutoff as part of your neighbor querying when performing +# this computation since it can otherwise be very sensitive. Specifically, it +# is highly recommended that you choose a number of neighbors that you +# specify a number of neighbors query that requests at least as many +# neighbors as the size of the motif you intend to test against. Otherwise, +# you will struggle to match the motif. However, this is not currently +# 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 +# +# 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. +# +# Args: +# system: +# Any object that is a valid argument to +# :class:`freud.locality.NeighborQuery.from_system`. +# motif ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): +# Vectors that make up the motif against which we are matching. +# 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). +# registration (bool, optional): +# If True, first use brute force registration to orient one set +# of environment vectors with respect to the other set such that +# it minimizes the RMSD between the two sets +# (Default value = :code:`False`). +# Returns: +# :math:`\left(N_{particles}\right)` :class:`numpy.ndarray`: +# 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) +# +# 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) +# +# 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 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): +# 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 +# frames of the :code:`points` as defined by the orientations +# :code:`orientations`. This computation accounts for the underlying +# symmetries of the reference frame as encoded in :code:`equiv_orientations`. +# +# 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. +# proj_vecs ((:math:`N_{vectors}`, 3) :class:`numpy.ndarray`): +# The set of projection vectors, defined in the query +# particles' reference frame, to calculate maximal local bond +# projections onto. +# 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`). +# (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=(None, 4)) +# +# equiv_orientations = freud.util._convert_array( +# 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)) +# 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._compute = self +# return nlist +# +# @_Compute._computed_property +# 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) +# +# @_Compute._computed_property +# def normed_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, 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) +# +# def __repr__(self): +# return ("freud.environment.{cls}()").format(cls=type(self).__name__) +# diff --git a/freud/environment.pyx b/freud/environment.pyx deleted file mode 100644 index 31883bd65..000000000 --- a/freud/environment.pyx +++ /dev/null @@ -1,1144 +0,0 @@ -# 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. - -r""" -The :class:`freud.environment` module contains functions which characterize the -local environments of particles in the system. These methods use the positions -and orientations of particles in the local neighborhood of a given particle to -characterize the particle environment. -""" - -from freud.errors import NO_DEFAULT_QUERY_ARGS_MESSAGE - -from cython.operator cimport dereference -from libcpp.map cimport map - -from freud.locality cimport _PairCompute, _SpatialHistogram -from freud.util cimport _Compute, quat, vec3 - -import warnings - -import numpy as np - -import freud.locality - -cimport numpy as np - -cimport freud._environment -cimport freud.box -cimport freud.locality -cimport freud.util - -# numpy must be initialized. When using numpy from C or Cython you must -# _always_ do that, or you will have segfaults -np.import_array() - - -cdef class BondOrder(_SpatialHistogram): - r"""Compute the bond orientational order diagram for the system of - particles. - - The bond orientational order diagram (BOOD) is a way of studying the - average local environments experienced by particles. In a BOOD, a particle - and its nearest neighbors (determined by either a prespecified number of - neighbors or simply a cutoff distance) are treated as connected by a bond - joining their centers. All of the bonds in the system are then binned by - their azimuthal (:math:`\theta`) and polar (:math:`\phi`) angles to - indicate the location of a particle's neighbors relative to itself. The - distance between the particle and its neighbor is only important when - determining whether it is counted as a neighbor, but is not part of the - BOOD; as such, the BOOD can be viewed as a projection of all bonds onto the - unit sphere. The resulting 2D histogram provides insight into how particles - are situated relative to one-another in a system. - - This class provides access to the classical BOOD as well as a few useful - variants. These variants can be accessed via the :code:`mode` arguments to - the :meth:`~BondOrder.compute` method. Available modes of calculation are: - - * :code:`'bod'` (Bond Order Diagram, *default*): - This mode constructs the default BOOD, which is the 2D histogram - containing the number of bonds formed through each azimuthal - :math:`\left( \theta \right)` and polar :math:`\left( \phi \right)` - angle. - - * :code:`'lbod'` (Local Bond Order Diagram): - In this mode, a particle's neighbors are rotated into the local frame of - the particle before the BOOD is calculated, *i.e.* the directions of - bonds are determined relative to the orientation of the particle rather - than relative to the global reference frame. An example of when this mode - would be useful is when a system is composed of multiple grains of the - same crystal; the normal BOOD would show twice as many peaks as expected, - but using this mode, the bonds would be superimposed. - - * :code:`'obcd'` (Orientation Bond Correlation Diagram): - This mode aims to quantify the degree of orientational as well as - translational ordering. As a first step, the rotation that would align a - particle's neighbor with the particle is calculated. Then, the neighbor - is rotated **around the central particle** by that amount, which actually - changes the direction of the bond. One example of how this mode could be - useful is in identifying plastic crystals, which exhibit translational - but not orientational ordering. Normally, the BOOD for a plastic crystal - would exhibit clear structure since there is translational order, but - with this mode, the neighbor positions would actually be modified, - resulting in an isotropic (disordered) BOOD. - - * :code:`'oocd'` (Orientation Orientation Correlation Diagram): - This mode is substantially different from the other modes. Rather than - compute the histogram of neighbor bonds, this mode instead computes a - histogram of the directors of neighboring particles, where the director - is defined as the basis vector :math:`\hat{z}` rotated by the neighbor's - quaternion. The directors are then rotated into the central particle's - reference frame. This mode provides insight into the local orientational - environment of particles, indicating, on average, how a particle's - neighbors are oriented. - - Args: - bins (unsigned int or sequence of length 2): - If an unsigned int, the number of bins in :math:`\theta` and - :math:`\phi`. If a sequence of two integers, interpreted as - :code:`(num_bins_theta, num_bins_phi)`. - mode (str, optional): - 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} - - def __cinit__(self, bins, str mode="bod"): - try: - n_bins_theta, n_bins_phi = bins - except TypeError: - 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) - - def __dealloc__(self): - del self.thisptr - - @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): - r"""Calculates the correlation function and adds to the current - histogram. - - 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. Uses identity quaternions if :code:`None` - (Default value = :code:`None`). - 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`). - 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). - reset (bool): - 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) - if orientations is None: - orientations = np.array([[1, 0, 0, 0]] * nq.points.shape[0]) - if query_orientations is None: - query_orientations = orientations - - orientations = freud.util._convert_array( - 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)) - 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) - - @_Compute._computed_property - def box(self): - """:class:`freud.box.Box`: Box used in the calculation.""" - return freud.box.BoxFromCPP(self.thisptr.getBox()) - - def __repr__(self): - 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)) - - @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 - - -cdef class LocalDescriptors(_PairCompute): - r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's - local environment. - - The resulting spherical harmonic array will be a complex-valued - array of shape :code:`(num_bonds, num_sphs)`. Spherical harmonic - calculation can be restricted to some number of nearest neighbors - through the :code:`max_num_neighbors` argument; if a particle has more - bonds than this number, the last one or more rows of bond spherical - harmonics for each particle will not be set. This feature is useful for - computing descriptors on the same system but with different subsets of - neighbors; a :class:`freud.locality.NeighborList` with the correct - ordering can then be reused in multiple calls to :meth:`~.compute` - with different values of :code:`max_num_neighbors` to compute descriptors - for different local neighborhoods with maximum efficiency. - - Args: - l_max (unsigned int): - Maximum spherical harmonic :math:`l` to consider. - negative_m (bool, optional): - True if we should also calculate :math:`Y_{lm}` for negative - :math:`m`. (Default value = :code:`True`) - mode (str, optional): - Orientation mode to use for environments, either - :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 (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} - - def __cinit__(self, l_max, negative_m=True, mode='neighborhood'): - cdef freud._environment.LocalDescriptorOrientation l_mode - 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): - r"""Calculates the local descriptors of bonds from a set of source - points to a set of destination points. - - Args: - system: - Any object that is a valid argument to - :class:`freud.locality.NeighborQuery.from_system`. - 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`). - orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): - Orientations associated with system points that are used to - calculate bonds. - 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). - max_num_neighbors (unsigned int, optional): - Hard limit on the maximum number of neighbors to use for each - 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) - - # 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 orientations is None: - raise RuntimeError( - ('Orientations must be given to orient LocalDescriptors ' - 'with particles\' orientations')) - - orientations = freud.util._convert_array( - orientations, shape=(nq.points.shape[0], 4)) - - l_orientations = orientations - l_orientations_ptr = &l_orientations[0, 0] - - 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) - 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._compute = self - return nlist - - @_Compute._computed_property - 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) - - @_Compute._computed_property - def num_sphs(self): - """unsigned int: The last number of spherical harmonics computed. This - is equal to the number of bonds in the last computation, which is at - most the number of :code:`points` multiplied by the lower of the - :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() - - @property - def l_max(self): - """unsigned int: The maximum spherical harmonic :math:`l` calculated - for.""" - return self.thisptr.getLMax() - - @property - def negative_m(self): - """bool: True if we also calculated :math:`Y_{lm}` for negative - :math:`m`.""" - return self.thisptr.getNegativeM() - - @property - def mode(self): - """str: Orientation mode to use for environments, either - :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() - for key, value in self.known_modes.items(): - if value == mode: - return key - - 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) - - -def _minimize_RMSD(box, ref_points, points, registration=False): - r"""Get the somewhat-optimal RMSD between the set of vectors ref_points - and the set of vectors points. - - Args: - ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): - Vectors that make up motif 1. - points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): - Vectors that make up motif 2. - registration (bool, optional): - If true, first use brute force registration to orient one set - of environment vectors with respect to the other set such that - it minimizes the RMSD between the two sets - (Default value = :code:`False`). - - Returns: - tuple (float, (:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): - A triplet that gives the associated min_rmsd, rotated (or not) - 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) - - 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] - - if nRef1 != nRef2: - raise ValueError( - ("The number of vectors in ref_points must MATCH" - "the number of vectors in points")) - - 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] - - -def _is_similar_motif(box, ref_points, points, threshold, registration=False): - r"""Test if the motif provided by ref_points is similar to the motif - provided by points. - - Args: - ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): - Vectors that make up motif 1. - points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): - Vectors that make up motif 2. - threshold (float): - Maximum magnitude of the vector difference between two vectors, - below which they are "matching". Typically, a good choice is - between 10% and 30% of the first well in the radial - distribution function (this has distance units). - registration (bool, optional): - If True, first use brute force registration to orient one set - of environment vectors with respect to the other set such that - it minimizes the RMSD between the two sets - (Default value = :code:`False`). - - Returns: - tuple ((:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): - A doublet that gives the rotated (or not) set of - :code:`points`, and the mapping between the vectors of - :code:`ref_points` and :code:`points` that will make them - 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) - - 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 - - if nRef1 != nRef2: - raise ValueError( - ("The number of vectors in ref_points must match" - "the number of vectors in points")) - - 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] - - -cdef class _MatchEnv(_PairCompute): - r"""Parent for environment matching methods. """ - cdef freud._environment.MatchEnv * matchptr - - def __cinit__(self, *args, **kwargs): - # Abstract class - pass - - @_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] - - def __repr__(self): - return ("freud.environment.{cls}()").format( - cls=type(self).__name__) - - -cdef 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): - r"""Determine clusters of particles with matching environments. - - An environment is defined by the bond vectors between a particle and its - neighbors as defined by :code:`env_neighbors`. - For example, :code:`env_neighbors= {'num_neighbors': 8}` means that every - particle's local environment is defined by its 8 nearest neighbors. - Then, each particle's environment is compared to the environments of - particles that satisfy a different cutoff parameter :code:`cluster_neighbors`. - For example, :code:`cluster_neighbors={'r_max': 3.0}` - means that the environment of each particle will be compared to the - environment of every particle within a distance of 3.0. - - Two environments are compared using `point-set registration - `_. - - The thresholding criterion we apply in order to determine if the two point sets - match is quite conservative: the two point sets match if and only if, - for every pair of matched points between the sets, the distance between - the matched pair is less than :code:`threshold`. - - When :code:`registration=False`, environments are not rotated prior to comparison - between them. Which pairs of vectors are compared depends on the order in - which the vectors of the environment are traversed. - - When :code:`registration=True`, we use the `Kabsch algorithm - `_ to find the optimal - rotation matrix mapping one environment onto another. The Kabsch - algorithm assumes a one-to-one correspondence between points in the two - environments. However, we typically do not know which point in one environment - corresponds to a particular point in the other environment. The brute force - solution is to try all possible correspondences, but the resulting - combinatorial explosion makes this problem infeasible, so we use the thresholding - criterion above to terminate the search when we have found a permutation - of points that results in a sufficiently low RMSD. - - .. note:: - - Using a distance cutoff for :code:`env_neighbors` could - lead to situations where the :code:`cluster_environments` - contain different numbers of neighbors. In this case, the - environments which have a number of neighbors less than - the environment with the maximum number of neighbors - :math:`k_{max}` will have their entry in :code:`cluster_environments` - padded with zero vectors. For example, a cluster environment - with :math:`m < k` neighbors, will have :math:`k - m` zero - vectors at the end of its entry in :code:`cluster_environments`. - - - .. warning:: - - All vectors of :code:`cluster_environments` and :code:`point_environments` - are defined with respect to the query particle. - Zero vectors are only used to pad the cluster vectors so that they - have the same shape. - In a future version of freud, zero-padding will be removed. - - .. warning:: - - Comparisons between two environments are only made when both - environments contain the same number of neighbors. - However, no warning will be given at runtime if mismatched - environments are provided for comparison. - - Example:: - - >>> import freud - >>> # Compute clusters of particles with matching environments - >>> box, points = freud.data.make_random_system(10, 100, seed=0) - >>> env_cluster = freud.environment.EnvironmentCluster() - >>> env_cluster.compute( - ... system = (box, points), - ... threshold=0.2, - ... cluster_neighbors={'num_neighbors': 6}, - ... registration=False) - freud.environment.EnvironmentCluster() - - Args: - system: - Any object that is a valid argument to - :class:`freud.locality.NeighborQuery.from_system`. - threshold (float): - Maximum magnitude of the vector difference between two vectors, - below which they are "matching". Typically, a good choice is - between 10% and 30% of the first well in the radial - distribution function (this has distance units). - cluster_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 - `_. - Defines the neighbors used for comparing different particle - environments (Default value: None). - env_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 - `_. - Defines the neighbors used as the environment of each particle. - If ``None``, the value provided for ``neighbors`` will be used - (Default value: None). - registration (bool, optional): - If True, first use brute force registration to orient one set - of environment vectors with respect to the other set such that - it minimizes the RMSD between the two sets. Enabling this - 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) - - 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) - 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) - - @_Compute._computed_property - def num_clusters(self): - """unsigned int: The number of clusters.""" - return self.thisptr.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] - - def plot(self, ax=None): - """Plot cluster distribution. - - Args: - ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If - :code:`None`, make a new figure and axis. - (Default value = :code:`None`) - - Returns: - (: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) - - 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): - r"""Find matches between local arrangements of a set of points and - a provided motif, as done in :cite:`Teich2019`. - - A particle's environment can only match the motif if it contains the - same number of neighbors as the motif. Any environment with a - 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 - - def compute(self, system, motif, threshold, env_neighbors=None, - registration=False): - r"""Determine which particles have local environments - matching the given environment motif. - - .. warning:: - - Comparisons between two environments are only made when both - environments contain the same number of neighbors. - However, no warning will be given at runtime if mismatched - environments are provided for comparison. - - Args: - system: - Any object that is a valid argument to - :class:`freud.locality.NeighborQuery.from_system`. - motif ((:math:`N_{points}`, 3) :class:`numpy.ndarray`): - Vectors that make up the motif against which we are matching. - threshold (float): - Maximum magnitude of the vector difference between two vectors, - below which they are "matching". Typically, a good choice is - between 10% and 30% of the first well in the radial - distribution function (this has distance units). - env_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 - `_. - Defines the environment of the query particles - (Default value: None). - registration (bool, optional): - If True, first use brute force registration to orient one set - of environment vectors with respect to the other set such that - 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) - - 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 - ) - - 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) - 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) - - -cdef class _EnvironmentRMSDMinimizer(_MatchEnv): - r"""Find linear transformations that map the environments of points onto a - motif. - - In general, it is recommended to specify a number of neighbors rather than - just a distance cutoff as part of your neighbor querying when performing - this computation since it can otherwise be very sensitive. Specifically, it - is highly recommended that you choose a number of neighbors that you - specify a number of neighbors query that requests at least as many - neighbors as the size of the motif you intend to test against. Otherwise, - you will struggle to match the motif. However, this is not currently - 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 - - 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. - - Args: - system: - Any object that is a valid argument to - :class:`freud.locality.NeighborQuery.from_system`. - motif ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): - Vectors that make up the motif against which we are matching. - 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). - registration (bool, optional): - If True, first use brute force registration to orient one set - of environment vectors with respect to the other set such that - it minimizes the RMSD between the two sets - (Default value = :code:`False`). - Returns: - :math:`\left(N_{particles}\right)` :class:`numpy.ndarray`: - 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) - - 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) - - 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) - - -class AngularSeparationNeighbor(_PairCompute): - r"""Calculates the minimum angles of separation between orientations and - query orientations.""" - def __init__(self): - self._cpp_obj = freud._environment.AngularSeparationNeighbor() - - 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 - - nq, nlist, qargs, 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 unsigned int n_equiv_orientations = l_equiv_orientations.shape[0] - - self._cpp_obj.compute( - nq._cpp_obj, - l_orientations, - query_points, - query_orientations, - equiv_orientations, - nlist._cpp_obj, - qargs._cpp_obj) - 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 self._cpp_obj.getAngles().toNumpyArray() - - 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._cpp_obj.getNList()) - nlist._compute = self - return nlist - - -class AngularSeparationGlobal(_Compute): - r"""Calculates the minimum angles of separation between orientations and - global orientations.""" - def __init__(self): - self._cpp_obj = freud._environment.AngularSeparationGlobal() - - 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)) - - 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 "freud.environment.{cls}()".format( - cls=type(self).__name__) - - -cdef 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): - 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 - frames of the :code:`points` as defined by the orientations - :code:`orientations`. This computation accounts for the underlying - symmetries of the reference frame as encoded in :code:`equiv_orientations`. - - 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. - proj_vecs ((:math:`N_{vectors}`, 3) :class:`numpy.ndarray`): - The set of projection vectors, defined in the query - particles' reference frame, to calculate maximal local bond - projections onto. - 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`). - (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=(None, 4)) - - equiv_orientations = freud.util._convert_array( - 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)) - 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._compute = self - return nlist - - @_Compute._computed_property - 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) - - @_Compute._computed_property - def normed_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, 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) - - def __repr__(self): - return ("freud.environment.{cls}()").format(cls=type(self).__name__) From 32f7bbf80523a2cb00412ae145bee83e74e47bbe Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 16:51:50 -0400 Subject: [PATCH 05/28] fix cmake config for environment --- freud/CMakeLists.txt | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index a9e47d3c4..133e8140c 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -204,6 +204,7 @@ set(python_files box.py cluster.py data.py + environment.py errors.py locality.py msd.py @@ -212,7 +213,7 @@ set(python_files interface.py plot.py util.py) -# density.py diffraction.py environment.py order.py) +# density.py diffraction.py order.py) copy_files_to_build("${python_files}" "freud" "*.py") From 128c40603c8c31a973725a8ba1e9e1ed91490187 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 17:13:34 -0400 Subject: [PATCH 06/28] tests on angular separations pass --- freud/environment/AngularSeparation.cc | 4 ++-- freud/environment/AngularSeparation.h | 6 ++--- .../export-AngularSeparationGlobal.cc | 13 +++++----- .../export-AngularSeparationNeighbor.cc | 24 +++++++++++-------- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/freud/environment/AngularSeparation.cc b/freud/environment/AngularSeparation.cc index 13a2ab0c0..1d36fa7fe 100644 --- a/freud/environment/AngularSeparation.cc +++ b/freud/environment/AngularSeparation.cc @@ -49,12 +49,12 @@ float computeMinSeparationAngle(const quat& ref_q, const quat& q, return min_angle; } -void AngularSeparationNeighbor::compute(const std::shared_ptr &nq, const quat* orientations, +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 std::shared_ptr &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. diff --git a/freud/environment/AngularSeparation.h b/freud/environment/AngularSeparation.h index b7909d0bb..ff921debd 100644 --- a/freud/environment/AngularSeparation.h +++ b/freud/environment/AngularSeparation.h @@ -60,11 +60,11 @@ class AngularSeparationNeighbor ~AngularSeparationNeighbor() = default; //! Compute the angular separation between neighbors - void compute(const std::shared_ptr &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 std::shared_ptr &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 std::shared_ptr> getAngles() const diff --git a/freud/environment/export-AngularSeparationGlobal.cc b/freud/environment/export-AngularSeparationGlobal.cc index 43a049a49..ab984cc33 100644 --- a/freud/environment/export-AngularSeparationGlobal.cc +++ b/freud/environment/export-AngularSeparationGlobal.cc @@ -18,10 +18,10 @@ template using nb_array = nanobind::ndarray; namespace wrap { -void compute(const std::shared_ptr &self, - const nb_array> &global_orientations, - const nb_array> & orientations, - const nb_array> & equiv_orientations) +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); @@ -29,7 +29,7 @@ void compute(const std::shared_ptr &self, 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()); - self->compute(global_orientations_data, n_global, orientations_data, n_points, equiv_orientations_data, n_equiv_orientations); + angular_separation->compute(global_orientations_data, n_global, orientations_data, n_points, equiv_orientations_data, n_equiv_orientations); } }; @@ -41,7 +41,8 @@ void export_AngularSeparationGlobal(nb::module_& module) nb::class_(module, "AngularSeparationGlobal") .def(nb::init<>()) .def("getAngles", &AngularSeparationGlobal::getAngles) - .def("compute", &AngularSeparationGlobal::compute); + .def("compute", &wrap::compute, nb::arg("global_orientations"), nb::arg("orientations"), + nb::arg("equiv_orientations")); } }; }; // namespace detail diff --git a/freud/environment/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc index e911f6108..31707ce0b 100644 --- a/freud/environment/export-AngularSeparationNeighbor.cc +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -20,14 +20,14 @@ using nb_array = nanobind::ndarray &self, - std::shared_ptr &nq, - nb_array> & orientations, - nb_array> & query_points, - nb_array> & query_orientations, - nb_array> & equiv_orientations, - std::shared_ptr &nlist, - const locality::QueryArgs qargs + std::shared_ptr& angular_separation, + std::shared_ptr nq, + nb_array>& orientations, + nb_array>& query_points, + nb_array>& query_orientations, + nb_array>& equiv_orientations, + std::shared_ptr nlist, + const locality::QueryArgs& qargs ) { auto* orientations_data = reinterpret_cast*>(orientations.data()); @@ -36,7 +36,7 @@ void compute( auto* equiv_orientations_data = reinterpret_cast*>(equiv_orientations.data()); unsigned int n_equiv_orientations = equiv_orientations.shape(0); unsigned int n_query_points = query_orientations.shape(0); - self->compute(nq, orientations_data, query_points_data, query_orientations_data, n_query_points, equiv_orientations_data, n_equiv_orientations, nlist, qargs); + angular_separation->compute(nq, orientations_data, query_points_data, query_orientations_data, n_query_points, equiv_orientations_data, n_equiv_orientations, nlist, qargs); } } @@ -45,9 +45,13 @@ 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", &AngularSeparationNeighbor::compute); + .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 From 19265bc3d64ed6b5af889bee23aeabc8f4521472 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 18:10:42 -0400 Subject: [PATCH 07/28] port local bond projection, tests pass --- freud/CMakeLists.txt | 7 +- freud/environment.py | 221 ++++++++---------- freud/environment/LocalBondProjection.cc | 25 +- freud/environment/LocalBondProjection.h | 20 +- .../export-AngularSeparationNeighbor.cc | 2 +- .../environment/export-LocalBondProjection.cc | 61 +++++ freud/environment/module-environment.cc | 2 + 7 files changed, 190 insertions(+), 148 deletions(-) create mode 100644 freud/environment/export-LocalBondProjection.cc diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 133e8140c..4e83cc062 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -47,8 +47,8 @@ add_library( environment/AngularSeparation.cc # environment/BondOrder.h # environment/BondOrder.cc - # environment/LocalBondProjection.h - # environment/LocalBondProjection.cc + environment/LocalBondProjection.h + environment/LocalBondProjection.cc # environment/LocalDescriptors.h # environment/LocalDescriptors.cc # environment/MatchEnv.h @@ -158,7 +158,8 @@ nanobind_add_module( _environment environment/module-environment.cc environment/export-AngularSeparationNeighbor.cc - environment/export-AngularSeparationGlobal.cc) + environment/export-AngularSeparationGlobal.cc + environment/export-LocalBondProjection.cc) target_link_libraries(_environment PUBLIC freud TBB::tbb) target_set_install_rpath(_environment) diff --git a/freud/environment.py b/freud/environment.py index fc2c74bfb..410bb429d 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -111,7 +111,7 @@ class attribute. ) return self - @_PairCompute._computed_property + @_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 @@ -121,7 +121,7 @@ def angles(self): def __repr__(self): return f"freud.environment.{type(self).__name__}()" - @_PairCompute._computed_property + @_Compute._computed_property def nlist(self): """:class:`freud.locality.NeighborList`: The neighbor list from the last compute.""" @@ -1033,123 +1033,100 @@ def __repr__(self): # # # -# cdef 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): -# 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 -# frames of the :code:`points` as defined by the orientations -# :code:`orientations`. This computation accounts for the underlying -# symmetries of the reference frame as encoded in :code:`equiv_orientations`. -# -# 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. -# proj_vecs ((:math:`N_{vectors}`, 3) :class:`numpy.ndarray`): -# The set of projection vectors, defined in the query -# particles' reference frame, to calculate maximal local bond -# projections onto. -# 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`). -# (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=(None, 4)) -# -# equiv_orientations = freud.util._convert_array( -# 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)) -# 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._compute = self -# return nlist -# -# @_Compute._computed_property -# 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) -# -# @_Compute._computed_property -# def normed_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, 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) -# -# def __repr__(self): -# return ("freud.environment.{cls}()").format(cls=type(self).__name__) -# +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. + """ + + def __init__(self): + self._cpp_obj = freud._environment.LocalBondProjection() + + def compute(self, system, orientations, proj_vecs, + query_points=None, equiv_orientations=np.array([[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 + frames of the :code:`points` as defined by the orientations + :code:`orientations`. This computation accounts for the underlying + symmetries of the reference frame as encoded in :code:`equiv_orientations`. + + 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. + proj_vecs ((:math:`N_{vectors}`, 3) :class:`numpy.ndarray`): + The set of projection vectors, defined in the query + particles' reference frame, to calculate maximal local bond + projections onto. + 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`). + (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 + nq, nlist, qargs, query_points, num_query_points = \ + self._preprocess_arguments(system, query_points, neighbors) + + orientations = freud.util._convert_array( + orientations, shape=(None, 4)) + + equiv_orientations = freud.util._convert_array( + equiv_orientations, shape=(None, 4)) + proj_vecs = freud.util._convert_array(proj_vecs, shape=(None, 3)) + + 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._cpp_obj.getNList()) + nlist._compute = self + return nlist + + @_Compute._computed_property + 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 self._cpp_obj.getProjections().toNumpyArray() + + @_Compute._computed_property + def normed_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, normalized by the length of the + bond.""" # noqa: E501 + return self._cpp_obj.getNormedProjections().toNumpyArray() + + def __repr__(self): + return (f"freud.environment.{type(self).__name__}()") + diff --git a/freud/environment/LocalBondProjection.cc b/freud/environment/LocalBondProjection.cc index 33d8e3f12..b6bec3871 100644 --- a/freud/environment/LocalBondProjection.cc +++ b/freud/environment/LocalBondProjection.cc @@ -46,45 +46,46 @@ float computeMaxProjection(const vec3& proj_vec, const vec3& local return max_proj; } -void LocalBondProjection::compute(const locality::NeighborQuery* nq, const quat* orientations, +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 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); // 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]; + 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; + (*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..01e4ab198 100644 --- a/freud/environment/LocalBondProjection.h +++ b/freud/environment/LocalBondProjection.h @@ -32,35 +32,35 @@ 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/export-AngularSeparationNeighbor.cc b/freud/environment/export-AngularSeparationNeighbor.cc index 31707ce0b..1a1697ae8 100644 --- a/freud/environment/export-AngularSeparationNeighbor.cc +++ b/freud/environment/export-AngularSeparationNeighbor.cc @@ -38,7 +38,7 @@ void compute( 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 detail { diff --git a/freud/environment/export-LocalBondProjection.cc b/freud/environment/export-LocalBondProjection.cc new file mode 100644 index 000000000..993928a21 --- /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" + + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { + +void compute( + std::shared_ptr& local_bond_projection, + std::shared_ptr nq, + nb_array>& orientations, + nb_array>& query_points, + nb_array>& projected_vectors, + nb_array>& equiv_orientations, + 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()); + unsigned int n_proj_vec = projected_vectors.shape(0); + unsigned int n_query_points = query_points.shape(0); + 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 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::locality diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc index 2e3e11b2c..68d498923 100644 --- a/freud/environment/module-environment.cc +++ b/freud/environment/module-environment.cc @@ -8,6 +8,7 @@ namespace freud::environment::detail { void export_AngularSeparationNeighbor(nanobind::module_& m); void export_AngularSeparationGlobal(nanobind::module_& m); +void export_LocalBondProjection(nanobind::module_& m); }; // namespace freud::environment::detail using namespace freud::environment::detail; @@ -15,4 +16,5 @@ NB_MODULE(_environment, module) // NOLINT(misc-use-anonymous-namespace): caused { export_AngularSeparationNeighbor(module); export_AngularSeparationGlobal(module); + export_LocalBondProjection(module); } From 6d22fd07e21c73ccd58955f91c13059eac60a7d8 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 18:14:47 -0400 Subject: [PATCH 08/28] update testing framework --- .github/workflows/test.yaml | 1 + 1 file changed, 1 insertion(+) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index e2b0c5535..f7a4cd0bb 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_interface.py -v pytest tests/test_msd_msd.py -v pytest tests/test_cluster.py -v From 8ef7d689136e926cd814f1c205d030147f5100f0 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Fri, 1 Nov 2024 18:22:30 -0400 Subject: [PATCH 09/28] begin porting local descriptors, python class is ported --- freud/environment.py | 314 ++++++++++++++++++++----------------------- 1 file changed, 149 insertions(+), 165 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index 410bb429d..d626da4d1 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -378,171 +378,155 @@ def __repr__(self): # return key # # -# cdef class LocalDescriptors(_PairCompute): -# r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's -# local environment. -# -# The resulting spherical harmonic array will be a complex-valued -# array of shape :code:`(num_bonds, num_sphs)`. Spherical harmonic -# calculation can be restricted to some number of nearest neighbors -# through the :code:`max_num_neighbors` argument; if a particle has more -# bonds than this number, the last one or more rows of bond spherical -# harmonics for each particle will not be set. This feature is useful for -# computing descriptors on the same system but with different subsets of -# neighbors; a :class:`freud.locality.NeighborList` with the correct -# ordering can then be reused in multiple calls to :meth:`~.compute` -# with different values of :code:`max_num_neighbors` to compute descriptors -# for different local neighborhoods with maximum efficiency. -# -# Args: -# l_max (unsigned int): -# Maximum spherical harmonic :math:`l` to consider. -# negative_m (bool, optional): -# True if we should also calculate :math:`Y_{lm}` for negative -# :math:`m`. (Default value = :code:`True`) -# mode (str, optional): -# Orientation mode to use for environments, either -# :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 (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} -# -# def __cinit__(self, l_max, negative_m=True, mode='neighborhood'): -# cdef freud._environment.LocalDescriptorOrientation l_mode -# 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): -# r"""Calculates the local descriptors of bonds from a set of source -# points to a set of destination points. -# -# Args: -# system: -# Any object that is a valid argument to -# :class:`freud.locality.NeighborQuery.from_system`. -# 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`). -# orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): -# Orientations associated with system points that are used to -# calculate bonds. -# 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). -# max_num_neighbors (unsigned int, optional): -# Hard limit on the maximum number of neighbors to use for each -# 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) -# -# # 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 orientations is None: -# raise RuntimeError( -# ('Orientations must be given to orient LocalDescriptors ' -# 'with particles\' orientations')) -# -# orientations = freud.util._convert_array( -# orientations, shape=(nq.points.shape[0], 4)) -# -# l_orientations = orientations -# l_orientations_ptr = &l_orientations[0, 0] -# -# 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) -# 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._compute = self -# return nlist -# -# @_Compute._computed_property -# 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) -# -# @_Compute._computed_property -# def num_sphs(self): -# """unsigned int: The last number of spherical harmonics computed. This -# is equal to the number of bonds in the last computation, which is at -# most the number of :code:`points` multiplied by the lower of the -# :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() -# -# @property -# def l_max(self): -# """unsigned int: The maximum spherical harmonic :math:`l` calculated -# for.""" -# return self.thisptr.getLMax() -# -# @property -# def negative_m(self): -# """bool: True if we also calculated :math:`Y_{lm}` for negative -# :math:`m`.""" -# return self.thisptr.getNegativeM() -# -# @property -# def mode(self): -# """str: Orientation mode to use for environments, either -# :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() -# for key, value in self.known_modes.items(): -# if value == mode: -# return key -# -# 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) -# -# +class LocalDescriptors(_PairCompute): + r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's + local environment. + + The resulting spherical harmonic array will be a complex-valued + array of shape :code:`(num_bonds, num_sphs)`. Spherical harmonic + calculation can be restricted to some number of nearest neighbors + through the :code:`max_num_neighbors` argument; if a particle has more + bonds than this number, the last one or more rows of bond spherical + harmonics for each particle will not be set. This feature is useful for + computing descriptors on the same system but with different subsets of + neighbors; a :class:`freud.locality.NeighborList` with the correct + ordering can then be reused in multiple calls to :meth:`~.compute` + with different values of :code:`max_num_neighbors` to compute descriptors + for different local neighborhoods with maximum efficiency. + + Args: + l_max (unsigned int): + Maximum spherical harmonic :math:`l` to consider. + negative_m (bool, optional): + True if we should also calculate :math:`Y_{lm}` for negative + :math:`m`. (Default value = :code:`True`) + mode (str, optional): + Orientation mode to use for environments, either + :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 (Default value = :code:`'neighborhood'`). + """ # noqa: E501 + + known_modes = {'neighborhood': freud._environment.LocalNeighborhood, + 'global': freud._environment.Global, + 'particle_local': freud._environment.ParticleLocal} + + def __init__(self, l_max, negative_m=True, mode='neighborhood'): + try: + l_mode = self.known_modes[mode] + except KeyError: + 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. + + Args: + system: + Any object that is a valid argument to + :class:`freud.locality.NeighborQuery.from_system`. + 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`). + orientations ((:math:`N_{points}`, 4) :class:`numpy.ndarray`): + Orientations associated with system points that are used to + calculate bonds. + 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). + max_num_neighbors (unsigned int, optional): + Hard limit on the maximum number of neighbors to use for each + particle for the given neighbor-finding algorithm. Uses + all neighbors if set to 0 (Default value = 0). + """ # noqa: E501 + 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. + if self.mode == 'particle_local': + if orientations is None: + 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)) + + 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._cpp_obj.getNList()) + nlist._compute = self + return nlist + + @_Compute._computed_property + def sph(self): + """:math:`\\left(N_{bonds}, \\text{SphWidth} \\right)` + :class:`numpy.ndarray`: The last computed spherical harmonic array.""" + return self._cpp_obj.getSph().toNumpyArray() + + @_Compute._computed_property + def num_sphs(self): + """unsigned int: The last number of spherical harmonics computed. This + is equal to the number of bonds in the last computation, which is at + most the number of :code:`points` multiplied by the lower of the + :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._cpp_obj.getNSphs() + + @property + def l_max(self): + """unsigned int: The maximum spherical harmonic :math:`l` calculated + for.""" + 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._cpp_obj.getNegativeM() + + @property + def mode(self): + """str: Orientation mode to use for environments, either + :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._cpp_obj.getMode() + for key, value in self.known_modes.items(): + if value == mode: + return key + return None + + def __repr__(self): + 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): # r"""Get the somewhat-optimal RMSD between the set of vectors ref_points # and the set of vectors points. From 4484dd566be362b09c37e20c9d6cb81059173bd8 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Mon, 4 Nov 2024 13:29:21 -0500 Subject: [PATCH 10/28] update build for LocalDescriptors --- freud/CMakeLists.txt | 8 +++++--- freud/environment/module-environment.cc | 2 ++ 2 files changed, 7 insertions(+), 3 deletions(-) diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 4e83cc062..471f25264 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -49,8 +49,8 @@ add_library( # environment/BondOrder.cc environment/LocalBondProjection.h environment/LocalBondProjection.cc - # environment/LocalDescriptors.h - # environment/LocalDescriptors.cc + environment/LocalDescriptors.h + environment/LocalDescriptors.cc # environment/MatchEnv.h # environment/MatchEnv.cc # environment/Registration.h @@ -159,7 +159,9 @@ nanobind_add_module( environment/module-environment.cc environment/export-AngularSeparationNeighbor.cc environment/export-AngularSeparationGlobal.cc - environment/export-LocalBondProjection.cc) + environment/export-LocalBondProjection.cc + environment/export-LocalDescriptors.cc +) target_link_libraries(_environment PUBLIC freud TBB::tbb) target_set_install_rpath(_environment) diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc index 68d498923..9a26428b5 100644 --- a/freud/environment/module-environment.cc +++ b/freud/environment/module-environment.cc @@ -9,6 +9,7 @@ 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); }; // namespace freud::environment::detail using namespace freud::environment::detail; @@ -17,4 +18,5 @@ NB_MODULE(_environment, module) // NOLINT(misc-use-anonymous-namespace): caused export_AngularSeparationNeighbor(module); export_AngularSeparationGlobal(module); export_LocalBondProjection(module); + export_LocalDescriptors(module); } From a7286ae84fff653899e856cb9d99caf24a44f0db Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Mon, 4 Nov 2024 15:16:06 -0500 Subject: [PATCH 11/28] port LocalDescriptors, enum in constructor is still TODO --- freud/environment/LocalDescriptors.cc | 31 +++++----- freud/environment/LocalDescriptors.h | 21 ++++--- freud/environment/export-LocalDescriptors.cc | 64 ++++++++++++++++++++ 3 files changed, 93 insertions(+), 23 deletions(-) create mode 100644 freud/environment/export-LocalDescriptors.cc diff --git a/freud/environment/LocalDescriptors.cc b/freud/environment/LocalDescriptors.cc index 1ddfee993..f26ac3f9a 100644 --- a/freud/environment/LocalDescriptors.cc +++ b/freud/environment/LocalDescriptors.cc @@ -18,9 +18,12 @@ 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, +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. @@ -30,14 +33,15 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3::max(); } - m_sphArray.prepare({m_nlist.getNumBonds(), getSphWidth()}); + m_sphArray = std::make_shared>>( + 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 +51,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 +102,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)); @@ -125,13 +128,13 @@ void LocalDescriptors::compute(const locality::NeighborQuery* nq, const vec3getNumBonds(); } }; }; // end namespace freud::environment diff --git a/freud/environment/LocalDescriptors.h b/freud/environment/LocalDescriptors.h index 153a6425d..742d46bf8 100644 --- a/freud/environment/LocalDescriptors.h +++ b/freud/environment/LocalDescriptors.h @@ -52,13 +52,16 @@ 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, - unsigned int n_query_points, const quat* orientations, - const freud::locality::NeighborList* nlist, locality::QueryArgs qargs, - unsigned int max_num_neighbors = 0); + void 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_neighbor); //! Get a reference to the last computed spherical harmonic array - const util::ManagedArray>& getSph() const + std::shared_ptr>> getSph() { return m_sphArray; } @@ -70,9 +73,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 @@ -89,11 +92,11 @@ class LocalDescriptors 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. + 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/export-LocalDescriptors.cc b/freud/environment/export-LocalDescriptors.cc new file mode 100644 index 000000000..495cb8370 --- /dev/null +++ b/freud/environment/export-LocalDescriptors.cc @@ -0,0 +1,64 @@ +// 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" + + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute(const std::shared_ptr& local_descriptors, + std::shared_ptr nq, + const nb_array>& query_points, + const unsigned int n_query_points, + const nb_array>& orientations, + 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 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("getSphArray", &LocalDescriptors::getSphArray) + .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"), + nb::arg("nlist").none(), + nb::arg("qargs"), nb::arg("max_num_neighbors")); +} + +}; }; // namespace detail +}; // namespace freud::locality From d01949cb1febeaf9458163fc38d05280c61bd92b Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 5 Nov 2024 13:45:06 -0500 Subject: [PATCH 12/28] export complex managed array --- freud/util/export-ManagedArray.cc | 1 + freud/util/export-ManagedArray.h | 1 + 2 files changed, 2 insertions(+) diff --git a/freud/util/export-ManagedArray.cc b/freud/util/export-ManagedArray.cc index 7e2a03e35..2a9a6e2b1 100644 --- a/freud/util/export-ManagedArray.cc +++ b/freud/util/export-ManagedArray.cc @@ -12,6 +12,7 @@ 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"); }; diff --git a/freud/util/export-ManagedArray.h b/freud/util/export-ManagedArray.h index 87d932c96..09a8cd9b7 100644 --- a/freud/util/export-ManagedArray.h +++ b/freud/util/export-ManagedArray.h @@ -7,6 +7,7 @@ #include "ManagedArray.h" #include "VectorMath.h" +#include #include #include #include From 0208c7418b302c05394bb6be30de6a4f14d250e4 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 5 Nov 2024 13:45:29 -0500 Subject: [PATCH 13/28] correctly export constructor and fix other exports --- freud/environment/export-LocalDescriptors.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/freud/environment/export-LocalDescriptors.cc b/freud/environment/export-LocalDescriptors.cc index 495cb8370..15b56d2ff 100644 --- a/freud/environment/export-LocalDescriptors.cc +++ b/freud/environment/export-LocalDescriptors.cc @@ -47,15 +47,15 @@ void export_LocalDescriptors(nb::module_& module) .export_values(); nb::class_(module, "LocalDescriptors") - .def(nb::init<>()) + .def(nb::init()) .def("getNList", &LocalDescriptors::getNList) - .def("getSphArray", &LocalDescriptors::getSphArray) + .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"), + 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")); } From 95499882687c8b4782e540a9f66e4e6fdda17305 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 5 Nov 2024 13:47:53 -0500 Subject: [PATCH 14/28] fix copy with new managed array --- freud/environment/LocalDescriptors.cc | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/freud/environment/LocalDescriptors.cc b/freud/environment/LocalDescriptors.cc index f26ac3f9a..c83e1db49 100644 --- a/freud/environment/LocalDescriptors.cc +++ b/freud/environment/LocalDescriptors.cc @@ -127,8 +127,7 @@ void LocalDescriptors::compute(const std::shared_ptr nq } sph_eval.compute(phi, theta); - - std::copy(sph_eval.begin(m_negative_m), sph_eval.end(), (*m_sphArray)[sphCount]); + std::copy(sph_eval.begin(m_negative_m), sph_eval.end(), m_sphArray->data() + sphCount); } } }); From f5f7b07d0829f45f158ea1ee0b84017473fef1f4 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 5 Nov 2024 13:49:36 -0500 Subject: [PATCH 15/28] add back missing const --- freud/environment/LocalDescriptors.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freud/environment/LocalDescriptors.h b/freud/environment/LocalDescriptors.h index 742d46bf8..5bf19b050 100644 --- a/freud/environment/LocalDescriptors.h +++ b/freud/environment/LocalDescriptors.h @@ -61,7 +61,7 @@ class LocalDescriptors unsigned int max_num_neighbor); //! Get a reference to the last computed spherical harmonic array - std::shared_ptr>> getSph() + std::shared_ptr>> getSph() const { return m_sphArray; } From 5eff9b0f6790694386b5ec6eb4975184621e7594 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 8 Nov 2024 14:14:29 -0500 Subject: [PATCH 16/28] Translate BondOrder --- freud/environment/BondOrder.cc | 28 ++++++++++++++++------------ freud/environment/BondOrder.h | 12 ++++++------ 2 files changed, 22 insertions(+), 18 deletions(-) diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index 376182c3a..340c35374 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -33,8 +33,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,14 +46,15 @@ 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), @@ -65,22 +66,25 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd void BondOrder::reduce() { - m_histogram.prepare(m_histogram.shape()); - m_bo_array.prepare(m_histogram.shape()); + // m_histogram.prepare(m_histogram.shape()); + // m_histogram = std::make_shared>(std::vector {m_histogram.shape()}); + // m_bo_array.prepare(m_histogram.shape()); + m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); + 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] = m_histogram[i] / (*m_sa_array)[i] / static_cast(m_frame_counter); }); } -const util::ManagedArray& BondOrder::getBondOrder() +const std::shared_ptr> BondOrder::getBondOrder() { return reduceAndReturn(m_bo_array); } -void BondOrder::accumulate(const locality::NeighborQuery* neighbor_query, quat* orientations, +void BondOrder::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, + unsigned int n_query_points, const std::shared_ptr& nlist, freud::locality::QueryArgs qargs) { accumulateGeneral(neighbor_query, query_points, n_query_points, nlist, qargs, diff --git a/freud/environment/BondOrder.h b/freud/environment/BondOrder.h index 45ca18644..060d4cd43 100644 --- a/freud/environment/BondOrder.h +++ b/freud/environment/BondOrder.h @@ -40,14 +40,14 @@ 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; - //! Get a reference to the last computed bond order - const util::ManagedArray& getBondOrder(); + //! Get a shared_ptr to the last computed bond order + const std::shared_ptr> getBondOrder(); BondOrderMode getMode() const { @@ -55,8 +55,8 @@ class BondOrder : public locality::BondHistogramCompute } private: - util::ManagedArray m_bo_array; //!< bond order array computed - util::ManagedArray m_sa_array; //!< surface area array computed + 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. }; From e72d8fce67b1d894fbae4cdd950a74eacf8fab84 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 8 Nov 2024 14:14:35 -0500 Subject: [PATCH 17/28] Progress on BondOrder export --- freud/environment/export-BondOrder.cc | 76 +++++++++++++++++++++++++++ 1 file changed, 76 insertions(+) create mode 100644 freud/environment/export-BondOrder.cc diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc new file mode 100644 index 000000000..3d9be677e --- /dev/null +++ b/freud/environment/export-BondOrder.cc @@ -0,0 +1,76 @@ +// 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 + +#include "BondOrder.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, + std::shared_ptr nlist, + const locality::QueryArgs& qargs + ) +{ + unsigned int const n_query_points = query_points.shape(0); + + 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("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::locality From 408b470bf20bca5372a8ae475fc6c976c8432f79 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 8 Nov 2024 14:14:47 -0500 Subject: [PATCH 18/28] Compile BondOrder --- freud/CMakeLists.txt | 5 +++-- freud/environment/module-environment.cc | 2 ++ 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 471f25264..f6ea06431 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -45,8 +45,8 @@ add_library( # environment environment/AngularSeparation.h environment/AngularSeparation.cc - # environment/BondOrder.h - # environment/BondOrder.cc + environment/BondOrder.h + environment/BondOrder.cc environment/LocalBondProjection.h environment/LocalBondProjection.cc environment/LocalDescriptors.h @@ -161,6 +161,7 @@ nanobind_add_module( environment/export-AngularSeparationGlobal.cc environment/export-LocalBondProjection.cc environment/export-LocalDescriptors.cc + environment/export-BondOrder.cc ) target_link_libraries(_environment PUBLIC freud TBB::tbb) target_set_install_rpath(_environment) diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc index 9a26428b5..a2672dbc1 100644 --- a/freud/environment/module-environment.cc +++ b/freud/environment/module-environment.cc @@ -10,6 +10,7 @@ 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); }; // namespace freud::environment::detail using namespace freud::environment::detail; @@ -19,4 +20,5 @@ NB_MODULE(_environment, module) // NOLINT(misc-use-anonymous-namespace): caused export_AngularSeparationGlobal(module); export_LocalBondProjection(module); export_LocalDescriptors(module); + export_BondOrder(module); } From 295a523dd8cebd3761172432d00eae03c219406a Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 8 Nov 2024 14:14:54 -0500 Subject: [PATCH 19/28] Python for BondOrder --- freud/environment.py | 385 +++++++++++++++++++++---------------------- 1 file changed, 191 insertions(+), 194 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index d626da4d1..0afae1949 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -8,6 +8,7 @@ characterize the particle environment. """ +import collections.abc import warnings import numpy as np @@ -184,200 +185,196 @@ 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. -# -# The bond orientational order diagram (BOOD) is a way of studying the -# average local environments experienced by particles. In a BOOD, a particle -# and its nearest neighbors (determined by either a prespecified number of -# neighbors or simply a cutoff distance) are treated as connected by a bond -# joining their centers. All of the bonds in the system are then binned by -# their azimuthal (:math:`\theta`) and polar (:math:`\phi`) angles to -# indicate the location of a particle's neighbors relative to itself. The -# distance between the particle and its neighbor is only important when -# determining whether it is counted as a neighbor, but is not part of the -# BOOD; as such, the BOOD can be viewed as a projection of all bonds onto the -# unit sphere. The resulting 2D histogram provides insight into how particles -# are situated relative to one-another in a system. -# -# This class provides access to the classical BOOD as well as a few useful -# variants. These variants can be accessed via the :code:`mode` arguments to -# the :meth:`~BondOrder.compute` method. Available modes of calculation are: -# -# * :code:`'bod'` (Bond Order Diagram, *default*): -# This mode constructs the default BOOD, which is the 2D histogram -# containing the number of bonds formed through each azimuthal -# :math:`\left( \theta \right)` and polar :math:`\left( \phi \right)` -# angle. -# -# * :code:`'lbod'` (Local Bond Order Diagram): -# In this mode, a particle's neighbors are rotated into the local frame of -# the particle before the BOOD is calculated, *i.e.* the directions of -# bonds are determined relative to the orientation of the particle rather -# than relative to the global reference frame. An example of when this mode -# would be useful is when a system is composed of multiple grains of the -# same crystal; the normal BOOD would show twice as many peaks as expected, -# but using this mode, the bonds would be superimposed. -# -# * :code:`'obcd'` (Orientation Bond Correlation Diagram): -# This mode aims to quantify the degree of orientational as well as -# translational ordering. As a first step, the rotation that would align a -# particle's neighbor with the particle is calculated. Then, the neighbor -# is rotated **around the central particle** by that amount, which actually -# changes the direction of the bond. One example of how this mode could be -# useful is in identifying plastic crystals, which exhibit translational -# but not orientational ordering. Normally, the BOOD for a plastic crystal -# would exhibit clear structure since there is translational order, but -# with this mode, the neighbor positions would actually be modified, -# resulting in an isotropic (disordered) BOOD. -# -# * :code:`'oocd'` (Orientation Orientation Correlation Diagram): -# This mode is substantially different from the other modes. Rather than -# compute the histogram of neighbor bonds, this mode instead computes a -# histogram of the directors of neighboring particles, where the director -# is defined as the basis vector :math:`\hat{z}` rotated by the neighbor's -# quaternion. The directors are then rotated into the central particle's -# reference frame. This mode provides insight into the local orientational -# environment of particles, indicating, on average, how a particle's -# neighbors are oriented. -# -# Args: -# bins (unsigned int or sequence of length 2): -# If an unsigned int, the number of bins in :math:`\theta` and -# :math:`\phi`. If a sequence of two integers, interpreted as -# :code:`(num_bins_theta, num_bins_phi)`. -# mode (str, optional): -# 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} -# -# def __cinit__(self, bins, str mode="bod"): -# try: -# n_bins_theta, n_bins_phi = bins -# except TypeError: -# 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) -# -# def __dealloc__(self): -# del self.thisptr -# -# @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): -# r"""Calculates the correlation function and adds to the current -# histogram. -# -# 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. Uses identity quaternions if :code:`None` -# (Default value = :code:`None`). -# 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`). -# 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). -# reset (bool): -# 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) -# if orientations is None: -# orientations = np.array([[1, 0, 0, 0]] * nq.points.shape[0]) -# if query_orientations is None: -# query_orientations = orientations -# -# orientations = freud.util._convert_array( -# 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)) -# 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) -# -# @_Compute._computed_property -# def box(self): -# """:class:`freud.box.Box`: Box used in the calculation.""" -# return freud.box.BoxFromCPP(self.thisptr.getBox()) -# -# def __repr__(self): -# 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)) -# -# @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 -# -# +class BondOrder(_SpatialHistogram): + r"""Compute the bond orientational order diagram for the system of + particles. + + The bond orientational order diagram (BOOD) is a way of studying the + average local environments experienced by particles. In a BOOD, a particle + and its nearest neighbors (determined by either a prespecified number of + neighbors or simply a cutoff distance) are treated as connected by a bond + joining their centers. All of the bonds in the system are then binned by + their azimuthal (:math:`\theta`) and polar (:math:`\phi`) angles to + indicate the location of a particle's neighbors relative to itself. The + distance between the particle and its neighbor is only important when + determining whether it is counted as a neighbor, but is not part of the + BOOD; as such, the BOOD can be viewed as a projection of all bonds onto the + unit sphere. The resulting 2D histogram provides insight into how particles + are situated relative to one-another in a system. + + This class provides access to the classical BOOD as well as a few useful + variants. These variants can be accessed via the :code:`mode` arguments to + the :meth:`~BondOrder.compute` method. Available modes of calculation are: + + * :code:`'bod'` (Bond Order Diagram, *default*): + This mode constructs the default BOOD, which is the 2D histogram + containing the number of bonds formed through each azimuthal + :math:`\left( \theta \right)` and polar :math:`\left( \phi \right)` + angle. + + * :code:`'lbod'` (Local Bond Order Diagram): + In this mode, a particle's neighbors are rotated into the local frame of + the particle before the BOOD is calculated, *i.e.* the directions of + bonds are determined relative to the orientation of the particle rather + than relative to the global reference frame. An example of when this mode + would be useful is when a system is composed of multiple grains of the + same crystal; the normal BOOD would show twice as many peaks as expected, + but using this mode, the bonds would be superimposed. + + * :code:`'obcd'` (Orientation Bond Correlation Diagram): + This mode aims to quantify the degree of orientational as well as + translational ordering. As a first step, the rotation that would align a + particle's neighbor with the particle is calculated. Then, the neighbor + is rotated **around the central particle** by that amount, which actually + changes the direction of the bond. One example of how this mode could be + useful is in identifying plastic crystals, which exhibit translational + but not orientational ordering. Normally, the BOOD for a plastic crystal + would exhibit clear structure since there is translational order, but + with this mode, the neighbor positions would actually be modified, + resulting in an isotropic (disordered) BOOD. + + * :code:`'oocd'` (Orientation Orientation Correlation Diagram): + This mode is substantially different from the other modes. Rather than + compute the histogram of neighbor bonds, this mode instead computes a + histogram of the directors of neighboring particles, where the director + is defined as the basis vector :math:`\hat{z}` rotated by the neighbor's + quaternion. The directors are then rotated into the central particle's + reference frame. This mode provides insight into the local orientational + environment of particles, indicating, on average, how a particle's + neighbors are oriented. + + Args: + bins (unsigned int or sequence of length 2): + If an unsigned int, the number of bins in :math:`\theta` and + :math:`\phi`. If a sequence of two integers, interpreted as + :code:`(num_bins_theta, num_bins_phi)`. + mode (str, optional): + Mode to calculate bond order. Options are :code:`'bod'`, + :code:`'lbod'`, :code:`'obcd'`, or :code:`'oocd'` + (Default value = :code:`'bod'`). + """ + + known_modes = {'bod': freud._environment.bod, + 'lbod': freud._environment.lbod, + 'obcd': freud._environment.obcd, + 'oocd': freud._environment.oocd} + + def __init__(self, bins, mode="bod"): + if isinstance(bins, collections.abc.Sequence): + n_bins_theta, n_bins_phi = bins + else: + n_bins_theta = n_bins_phi = bins + + try: + l_mode = self.known_modes[mode] + except KeyError as err: + msg = f"Unknown BondOrder mode: {mode}" + raise ValueError(msg) from err + + 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): + r"""Calculates the correlation function and adds to the current + histogram. + + 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. Uses identity quaternions if :code:`None` + (Default value = :code:`None`). + 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`). + 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). + reset (bool): + Whether to erase the previously computed values before adding + the new computation; if False, will accumulate data (Default + value: True). + """ + if reset: + self._reset() + + 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 + + orientations = freud.util._convert_array(orientations, shape=(nq.points.shape[0], 4)) + query_orientations = freud.util._convert_array(query_orientations, shape=(num_query_points, 4)) + + # print(f"orientations {orientations}") + print(f"qp {query_points}") + # print(f"qo {query_orientations}") + + + self._cpp_obj.accumulate( + # &l_orientations[0, 0], + # &l_query_points[0, 0], + # &l_query_orientations[0, 0], + nq._cpp_obj, + orientations, + query_points, + query_orientations, + # num_query_points, + 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._cpp_obj.getBox()) + + def __repr__(self): + 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)) + + @property + def mode(self): + """str: Bond order mode.""" + mode = self._cpp_obj.getMode() + # for key, value in self.known_modes.items(): + # if value == mode: + # return key + return mode + + class LocalDescriptors(_PairCompute): r"""Compute a set of descriptors (a numerical "fingerprint") of a particle's local environment. From 494fa15fd51c56bb75b71852b997380864fc7390 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 8 Nov 2024 15:10:16 -0500 Subject: [PATCH 20/28] Progress on BondOrder --- freud/environment.py | 19 +++++++------------ freud/environment/BondOrder.cc | 13 ++++++++----- freud/environment/export-BondOrder.cc | 18 +++++++++++++++--- 3 files changed, 30 insertions(+), 20 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index 0afae1949..6f709d579 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -323,24 +323,19 @@ def compute(self, system, orientations=None, query_points=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)) query_orientations = freud.util._convert_array(query_orientations, shape=(num_query_points, 4)) + query_points = freud.util._convert_array(query_points, shape=(num_query_points, 3)) - # print(f"orientations {orientations}") - print(f"qp {query_points}") - # print(f"qo {query_orientations}") - self._cpp_obj.accumulate( - # &l_orientations[0, 0], - # &l_query_points[0, 0], - # &l_query_orientations[0, 0], nq._cpp_obj, orientations, query_points, query_orientations, - # num_query_points, nlist._cpp_obj, qargs._cpp_obj ) @@ -352,7 +347,7 @@ def bond_order(self): # return freud.util.make_managed_numpy_array( # &self.thisptr.getBondOrder(), # freud.util.arr_type_t.FLOAT) - return self._cpp_obj.getBondOrder().toNumpyArray() + return self._cpp_obj.getBondOrder()#.toNumpyArray() @_Compute._computed_property def box(self): @@ -369,9 +364,9 @@ def __repr__(self): def mode(self): """str: Bond order mode.""" mode = self._cpp_obj.getMode() - # for key, value in self.known_modes.items(): - # if value == mode: - # return key + for key, value in self.known_modes.items(): + if value == mode: + return key return mode diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index 340c35374..11bc50712 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -57,18 +57,21 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd (*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); + 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_local_histograms = BondHistogram::ThreadLocalHistogram(m_histogram); } void BondOrder::reduce() { + // TODO: previously, we could prepare the histogram: but not anymore? // m_histogram.prepare(m_histogram.shape()); - // m_histogram = std::make_shared>(std::vector {m_histogram.shape()}); - // m_bo_array.prepare(m_histogram.shape()); + // m_histogram = std::make_shared>(std::vector {m_histogram.shape()}); + m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc index 3d9be677e..31fba4b50 100644 --- a/freud/environment/export-BondOrder.cc +++ b/freud/environment/export-BondOrder.cc @@ -1,12 +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. +#include #include #include #include #include // NOLINT(misc-include-cleaner): used implicitly -// #include // NOLINT(misc-include-cleaner): used implicitly -// #include +#include // NOLINT(misc-include-cleaner): used implicitly +#include #include "BondOrder.h" @@ -30,11 +31,20 @@ void accumulateBondOrder(const std::shared_ptr& self, ) { 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 ); @@ -57,6 +67,8 @@ void export_BondOrder(nb::module_& module) nb::class_(module, "BondOrder") .def(nb::init()) .def("getBondOrder", &BondOrder::getBondOrder) + .def("getBinCounts", &BondOrder::getBinCounts) + .def("getAxisSizes", &BondOrder::getAxisSizes) .def("getMode", &BondOrder::getMode) .def("accumulate", &wrap::accumulateBondOrder, nanobind::arg("nq").none(), From 42b24e40b4ccb493c41e5befdd44ed34218021e5 Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 12 Nov 2024 15:06:12 -0500 Subject: [PATCH 21/28] fix managed array issue for LocalDescriptors --- freud/environment/LocalDescriptors.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/freud/environment/LocalDescriptors.cc b/freud/environment/LocalDescriptors.cc index c83e1db49..df9d11010 100644 --- a/freud/environment/LocalDescriptors.cc +++ b/freud/environment/LocalDescriptors.cc @@ -34,7 +34,7 @@ void LocalDescriptors::compute(const std::shared_ptr nq max_num_neighbors = std::numeric_limits::max(); } m_sphArray = std::make_shared>>( - m_nlist->getNumBonds() * getSphWidth()); + std::vector{m_nlist->getNumBonds(), getSphWidth()}); util::forLoopWrapper(0, nq->getNPoints(), [&](size_t begin, size_t end) { fsph::PointSPHEvaluator sph_eval(m_l_max); @@ -127,6 +127,7 @@ void LocalDescriptors::compute(const std::shared_ptr nq } sph_eval.compute(phi, theta); + // copy results from sph_eval to the sphArray std::copy(sph_eval.begin(m_negative_m), sph_eval.end(), m_sphArray->data() + sphCount); } } From e6fa0ba909aaa7b414d81b6668fe30f991234a8f Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 12 Nov 2024 17:07:17 -0500 Subject: [PATCH 22/28] start porting environment matching, python file ported --- freud/CMakeLists.txt | 7 +- freud/environment.py | 938 +++++++++++------------- freud/environment/export-MatchEnv.cc | 65 ++ freud/environment/module-environment.cc | 2 + 4 files changed, 518 insertions(+), 494 deletions(-) create mode 100644 freud/environment/export-MatchEnv.cc diff --git a/freud/CMakeLists.txt b/freud/CMakeLists.txt index 9ea831571..86a97011f 100644 --- a/freud/CMakeLists.txt +++ b/freud/CMakeLists.txt @@ -50,9 +50,9 @@ add_library( environment/LocalBondProjection.cc environment/LocalDescriptors.h environment/LocalDescriptors.cc - # environment/MatchEnv.h - # environment/MatchEnv.cc - # environment/Registration.h + environment/MatchEnv.h + environment/MatchEnv.cc + environment/Registration.h # locality locality/AABB.h locality/AABBQuery.cc @@ -166,6 +166,7 @@ nanobind_add_module( 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) diff --git a/freud/environment.py b/freud/environment.py index 6f709d579..b00a90a5a 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -398,7 +398,7 @@ class LocalDescriptors(_PairCompute): neighborhood, :code:`'particle_local'` to use the given particle orientations, or :code:`'global'` to not rotate environments (Default value = :code:`'neighborhood'`). - """ # noqa: E501 + """ known_modes = {'neighborhood': freud._environment.LocalNeighborhood, 'global': freud._environment.Global, @@ -519,496 +519,452 @@ def __repr__(self): f"negative_m={self.negative_m}, mode='{self.mode}')") -# def _minimize_RMSD(box, ref_points, points, registration=False): -# r"""Get the somewhat-optimal RMSD between the set of vectors ref_points -# and the set of vectors points. -# -# Args: -# ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): -# Vectors that make up motif 1. -# points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): -# Vectors that make up motif 2. -# registration (bool, optional): -# If true, first use brute force registration to orient one set -# of environment vectors with respect to the other set such that -# it minimizes the RMSD between the two sets -# (Default value = :code:`False`). -# -# Returns: -# tuple (float, (:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): -# A triplet that gives the associated min_rmsd, rotated (or not) -# 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) -# -# 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] -# -# if nRef1 != nRef2: -# raise ValueError( -# ("The number of vectors in ref_points must MATCH" -# "the number of vectors in points")) -# -# 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] -# -# -# def _is_similar_motif(box, ref_points, points, threshold, registration=False): -# r"""Test if the motif provided by ref_points is similar to the motif -# provided by points. -# -# Args: -# ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): -# Vectors that make up motif 1. -# points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): -# Vectors that make up motif 2. -# threshold (float): -# Maximum magnitude of the vector difference between two vectors, -# below which they are "matching". Typically, a good choice is -# between 10% and 30% of the first well in the radial -# distribution function (this has distance units). -# registration (bool, optional): -# If True, first use brute force registration to orient one set -# of environment vectors with respect to the other set such that -# it minimizes the RMSD between the two sets -# (Default value = :code:`False`). -# -# Returns: -# tuple ((:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): -# A doublet that gives the rotated (or not) set of -# :code:`points`, and the mapping between the vectors of -# :code:`ref_points` and :code:`points` that will make them -# 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) -# -# 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 -# -# if nRef1 != nRef2: -# raise ValueError( -# ("The number of vectors in ref_points must match" -# "the number of vectors in points")) -# -# 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] -# -# -# cdef class _MatchEnv(_PairCompute): -# r"""Parent for environment matching methods. """ -# cdef freud._environment.MatchEnv * matchptr -# -# def __cinit__(self, *args, **kwargs): -# # Abstract class -# pass -# -# @_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] -# -# def __repr__(self): -# return ("freud.environment.{cls}()").format( -# cls=type(self).__name__) -# -# -# cdef 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): -# r"""Determine clusters of particles with matching environments. -# -# An environment is defined by the bond vectors between a particle and its -# neighbors as defined by :code:`env_neighbors`. -# For example, :code:`env_neighbors= {'num_neighbors': 8}` means that every -# particle's local environment is defined by its 8 nearest neighbors. -# Then, each particle's environment is compared to the environments of -# particles that satisfy a different cutoff parameter :code:`cluster_neighbors`. -# For example, :code:`cluster_neighbors={'r_max': 3.0}` -# means that the environment of each particle will be compared to the -# environment of every particle within a distance of 3.0. -# -# Two environments are compared using `point-set registration -# `_. -# -# The thresholding criterion we apply in order to determine if the two point sets -# match is quite conservative: the two point sets match if and only if, -# for every pair of matched points between the sets, the distance between -# the matched pair is less than :code:`threshold`. -# -# When :code:`registration=False`, environments are not rotated prior to comparison -# between them. Which pairs of vectors are compared depends on the order in -# which the vectors of the environment are traversed. -# -# When :code:`registration=True`, we use the `Kabsch algorithm -# `_ to find the optimal -# rotation matrix mapping one environment onto another. The Kabsch -# algorithm assumes a one-to-one correspondence between points in the two -# environments. However, we typically do not know which point in one environment -# corresponds to a particular point in the other environment. The brute force -# solution is to try all possible correspondences, but the resulting -# combinatorial explosion makes this problem infeasible, so we use the thresholding -# criterion above to terminate the search when we have found a permutation -# of points that results in a sufficiently low RMSD. -# -# .. note:: -# -# Using a distance cutoff for :code:`env_neighbors` could -# lead to situations where the :code:`cluster_environments` -# contain different numbers of neighbors. In this case, the -# environments which have a number of neighbors less than -# the environment with the maximum number of neighbors -# :math:`k_{max}` will have their entry in :code:`cluster_environments` -# padded with zero vectors. For example, a cluster environment -# with :math:`m < k` neighbors, will have :math:`k - m` zero -# vectors at the end of its entry in :code:`cluster_environments`. -# -# -# .. warning:: -# -# All vectors of :code:`cluster_environments` and :code:`point_environments` -# are defined with respect to the query particle. -# Zero vectors are only used to pad the cluster vectors so that they -# have the same shape. -# In a future version of freud, zero-padding will be removed. -# -# .. warning:: -# -# Comparisons between two environments are only made when both -# environments contain the same number of neighbors. -# However, no warning will be given at runtime if mismatched -# environments are provided for comparison. -# -# Example:: -# -# >>> import freud -# >>> # Compute clusters of particles with matching environments -# >>> box, points = freud.data.make_random_system(10, 100, seed=0) -# >>> env_cluster = freud.environment.EnvironmentCluster() -# >>> env_cluster.compute( -# ... system = (box, points), -# ... threshold=0.2, -# ... cluster_neighbors={'num_neighbors': 6}, -# ... registration=False) -# freud.environment.EnvironmentCluster() -# -# Args: -# system: -# Any object that is a valid argument to -# :class:`freud.locality.NeighborQuery.from_system`. -# threshold (float): -# Maximum magnitude of the vector difference between two vectors, -# below which they are "matching". Typically, a good choice is -# between 10% and 30% of the first well in the radial -# distribution function (this has distance units). -# cluster_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 -# `_. -# Defines the neighbors used for comparing different particle -# environments (Default value: None). -# env_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 -# `_. -# Defines the neighbors used as the environment of each particle. -# If ``None``, the value provided for ``neighbors`` will be used -# (Default value: None). -# registration (bool, optional): -# If True, first use brute force registration to orient one set -# of environment vectors with respect to the other set such that -# it minimizes the RMSD between the two sets. Enabling this -# 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) -# -# 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) -# 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) -# -# @_Compute._computed_property -# def num_clusters(self): -# """unsigned int: The number of clusters.""" -# return self.thisptr.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] -# -# def plot(self, ax=None): -# """Plot cluster distribution. -# -# Args: -# ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If -# :code:`None`, make a new figure and axis. -# (Default value = :code:`None`) -# -# Returns: -# (: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) -# -# 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): -# r"""Find matches between local arrangements of a set of points and -# a provided motif, as done in :cite:`Teich2019`. -# -# A particle's environment can only match the motif if it contains the -# same number of neighbors as the motif. Any environment with a -# 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 -# -# def compute(self, system, motif, threshold, env_neighbors=None, -# registration=False): -# r"""Determine which particles have local environments -# matching the given environment motif. -# -# .. warning:: -# -# Comparisons between two environments are only made when both -# environments contain the same number of neighbors. -# However, no warning will be given at runtime if mismatched -# environments are provided for comparison. -# -# Args: -# system: -# Any object that is a valid argument to -# :class:`freud.locality.NeighborQuery.from_system`. -# motif ((:math:`N_{points}`, 3) :class:`numpy.ndarray`): -# Vectors that make up the motif against which we are matching. -# threshold (float): -# Maximum magnitude of the vector difference between two vectors, -# below which they are "matching". Typically, a good choice is -# between 10% and 30% of the first well in the radial -# distribution function (this has distance units). -# env_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 -# `_. -# Defines the environment of the query particles -# (Default value: None). -# registration (bool, optional): -# If True, first use brute force registration to orient one set -# of environment vectors with respect to the other set such that -# 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) -# -# 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 -# ) -# -# 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) -# 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) -# -# -# cdef class _EnvironmentRMSDMinimizer(_MatchEnv): -# r"""Find linear transformations that map the environments of points onto a -# motif. -# -# In general, it is recommended to specify a number of neighbors rather than -# just a distance cutoff as part of your neighbor querying when performing -# this computation since it can otherwise be very sensitive. Specifically, it -# is highly recommended that you choose a number of neighbors that you -# specify a number of neighbors query that requests at least as many -# neighbors as the size of the motif you intend to test against. Otherwise, -# you will struggle to match the motif. However, this is not currently -# 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 -# -# 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. -# -# Args: -# system: -# Any object that is a valid argument to -# :class:`freud.locality.NeighborQuery.from_system`. -# motif ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): -# Vectors that make up the motif against which we are matching. -# 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). -# registration (bool, optional): -# If True, first use brute force registration to orient one set -# of environment vectors with respect to the other set such that -# it minimizes the RMSD between the two sets -# (Default value = :code:`False`). -# Returns: -# :math:`\left(N_{particles}\right)` :class:`numpy.ndarray`: -# 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) -# -# 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) -# -# 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) -# -# -# +def _minimize_RMSD(box, ref_points, points, registration=False): + r"""Get the somewhat-optimal RMSD between the set of vectors ref_points + and the set of vectors points. + + Args: + ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): + Vectors that make up motif 1. + points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): + Vectors that make up motif 2. + registration (bool, optional): + If true, first use brute force registration to orient one set + of environment vectors with respect to the other set such that + it minimizes the RMSD between the two sets + (Default value = :code:`False`). + + Returns: + tuple (float, (:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): + A triplet that gives the associated min_rmsd, rotated (or not) + set of points, and the mapping between the vectors of + ref_points and points that somewhat minimizes the RMSD. + """ # noqa: E501 + 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)) + + nRef1 = ref_points.shape[0] + nRef2 = points.shape[0] + + if nRef1 != nRef2: + msg = ( + "The number of vectors in ref_points must MATCH" + "the number of vectors in points" + ) + raise ValueError(msg) + + min_rmsd = -1 + 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): + r"""Test if the motif provided by ref_points is similar to the motif + provided by points. + + Args: + ref_points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): + Vectors that make up motif 1. + points ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): + Vectors that make up motif 2. + threshold (float): + Maximum magnitude of the vector difference between two vectors, + below which they are "matching". Typically, a good choice is + between 10% and 30% of the first well in the radial + distribution function (this has distance units). + registration (bool, optional): + If True, first use brute force registration to orient one set + of environment vectors with respect to the other set such that + it minimizes the RMSD between the two sets + (Default value = :code:`False`). + + Returns: + tuple ((:math:`\left(N_{particles}, 3\right)` :class:`numpy.ndarray`), map[int, int]): + A doublet that gives the rotated (or not) set of + :code:`points`, and the mapping between the vectors of + :code:`ref_points` and :code:`points` that will make them + correspond to each other. Empty if they do not correspond to + each other. + """ # noqa: E501 + 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)) + + nRef1 = ref_points.shape[0] + nRef2 = points.shape[0] + threshold_sq = threshold*threshold + + if nRef1 != nRef2: + msg = ( + "The number of vectors in ref_points must match" + "the number of vectors in points" + ) + raise ValueError(msg) + + vec_map = freud._environment.isSimilar( + b._cpp_obj, + ref_points, + points, + nRef1, + threshold_sq, + registration) + return [np.asarray(points), vec_map] + + +class _MatchEnv(_PairCompute): + r"""Parent for environment matching methods. """ + + def __init__(self, *args, **kwargs): + # Abstract class + 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._cpp_obj.getPointEnvironments() + return [np.asarray([[p.x, p.y, p.z] for p in env]) + for env in envs] + + def __repr__(self): + return (f"freud.environment.{type(self).__name__}()") + + +class EnvironmentCluster(_MatchEnv): + r"""Clusters particles according to whether their local environments match + or not, using various shape matching metrics defined in :cite:`Teich2019`. + """ + + def __init__(self): + 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 + neighbors as defined by :code:`env_neighbors`. + For example, :code:`env_neighbors= {'num_neighbors': 8}` means that every + particle's local environment is defined by its 8 nearest neighbors. + Then, each particle's environment is compared to the environments of + particles that satisfy a different cutoff parameter :code:`cluster_neighbors`. + For example, :code:`cluster_neighbors={'r_max': 3.0}` + means that the environment of each particle will be compared to the + environment of every particle within a distance of 3.0. + + Two environments are compared using `point-set registration + `_. + + The thresholding criterion we apply in order to determine if the two point sets + match is quite conservative: the two point sets match if and only if, + for every pair of matched points between the sets, the distance between + the matched pair is less than :code:`threshold`. + + When :code:`registration=False`, environments are not rotated prior to comparison + between them. Which pairs of vectors are compared depends on the order in + which the vectors of the environment are traversed. + + When :code:`registration=True`, we use the `Kabsch algorithm + `_ to find the optimal + rotation matrix mapping one environment onto another. The Kabsch + algorithm assumes a one-to-one correspondence between points in the two + environments. However, we typically do not know which point in one environment + corresponds to a particular point in the other environment. The brute force + solution is to try all possible correspondences, but the resulting + combinatorial explosion makes this problem infeasible, so we use the thresholding + criterion above to terminate the search when we have found a permutation + of points that results in a sufficiently low RMSD. + + .. note:: + + Using a distance cutoff for :code:`env_neighbors` could + lead to situations where the :code:`cluster_environments` + contain different numbers of neighbors. In this case, the + environments which have a number of neighbors less than + the environment with the maximum number of neighbors + :math:`k_{max}` will have their entry in :code:`cluster_environments` + padded with zero vectors. For example, a cluster environment + with :math:`m < k` neighbors, will have :math:`k - m` zero + vectors at the end of its entry in :code:`cluster_environments`. + + + .. warning:: + + All vectors of :code:`cluster_environments` and :code:`point_environments` + are defined with respect to the query particle. + Zero vectors are only used to pad the cluster vectors so that they + have the same shape. + In a future version of freud, zero-padding will be removed. + + .. warning:: + + Comparisons between two environments are only made when both + environments contain the same number of neighbors. + However, no warning will be given at runtime if mismatched + environments are provided for comparison. + + Example:: + + >>> import freud + >>> # Compute clusters of particles with matching environments + >>> box, points = freud.data.make_random_system(10, 100, seed=0) + >>> env_cluster = freud.environment.EnvironmentCluster() + >>> env_cluster.compute( + ... system = (box, points), + ... threshold=0.2, + ... cluster_neighbors={'num_neighbors': 6}, + ... registration=False) + freud.environment.EnvironmentCluster() + + Args: + system: + Any object that is a valid argument to + :class:`freud.locality.NeighborQuery.from_system`. + threshold (float): + Maximum magnitude of the vector difference between two vectors, + below which they are "matching". Typically, a good choice is + between 10% and 30% of the first well in the radial + distribution function (this has distance units). + cluster_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 + `_. + Defines the neighbors used for comparing different particle + environments (Default value: None). + env_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 + `_. + Defines the neighbors used as the environment of each particle. + If ``None``, the value provided for ``neighbors`` will be used + (Default value: None). + registration (bool, optional): + If True, first use brute force registration to orient one set + of environment vectors with respect to the other set such that + it minimizes the RMSD between the two sets. Enabling this + option incurs a significant performance penalty. + (Default value = :code:`False`) + """ # noqa: E501 + 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._cpp_obj.compute( + nq._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 self._cpp_obj.getClusterIdx().toNumpyArray() + + @_Compute._computed_property + def num_clusters(self): + """unsigned int: The number of clusters.""" + 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._cpp_obj.getClusterEnvironments() + return [np.asarray([[p.x, p.y, p.z] for p in env]) + for env in envs] + + def plot(self, ax=None): + """Plot cluster distribution. + + Args: + ax (:class:`matplotlib.axes.Axes`, optional): Axis to plot on. If + :code:`None`, make a new figure and axis. + (Default value = :code:`None`) + + Returns: + (: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) + + def _repr_png_(self): + try: + import freud.plot + return freud.plot._ax_to_bytes(self.plot()) + except (AttributeError, ImportError): + return None + + +class EnvironmentMotifMatch(_MatchEnv): + r"""Find matches between local arrangements of a set of points and + a provided motif, as done in :cite:`Teich2019`. + + A particle's environment can only match the motif if it contains the + same number of neighbors as the motif. Any environment with a + different number of neighbors than the motif will always fail to match + the motif. See :class:`freud.environment.EnvironmentCluster.compute` for + the matching criterion. + """ + + def __init__(self): + self._cpp_obj = freud._environment.EnvironmentMotifMatch() + + def compute(self, system, motif, threshold, env_neighbors=None, + registration=False): + r"""Determine which particles have local environments + matching the given environment motif. + + .. warning:: + + Comparisons between two environments are only made when both + environments contain the same number of neighbors. + However, no warning will be given at runtime if mismatched + environments are provided for comparison. + + Args: + system: + Any object that is a valid argument to + :class:`freud.locality.NeighborQuery.from_system`. + motif ((:math:`N_{points}`, 3) :class:`numpy.ndarray`): + Vectors that make up the motif against which we are matching. + threshold (float): + Maximum magnitude of the vector difference between two vectors, + below which they are "matching". Typically, a good choice is + between 10% and 30% of the first well in the radial + distribution function (this has distance units). + env_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 + `_. + Defines the environment of the query particles + (Default value: None). + registration (bool, optional): + If True, first use brute force registration to orient one set + of environment vectors with respect to the other set such that + it minimizes the RMSD between the two sets + (Default value = False). + """ + 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, + stacklevel=2 + ) + 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 self._cpp_obj.getMatches().toNumpyArray() + + +class _EnvironmentRMSDMinimizer(_MatchEnv): + r"""Find linear transformations that map the environments of points onto a + motif. + + In general, it is recommended to specify a number of neighbors rather than + just a distance cutoff as part of your neighbor querying when performing + this computation since it can otherwise be very sensitive. Specifically, it + is highly recommended that you choose a number of neighbors that you + specify a number of neighbors query that requests at least as many + neighbors as the size of the motif you intend to test against. Otherwise, + you will struggle to match the motif. However, this is not currently + enforced (but we could add a warning to the compute...). + """ + + def __init__(self): + self._cpp_obj = freud._environment.EnvironmentRMSDMinimizer() + + 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. + + Args: + system: + Any object that is a valid argument to + :class:`freud.locality.NeighborQuery.from_system`. + motif ((:math:`N_{particles}`, 3) :class:`numpy.ndarray`): + Vectors that make up the motif against which we are matching. + 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). + registration (bool, optional): + If True, first use brute force registration to orient one set + of environment vectors with respect to the other set such that + it minimizes the RMSD between the two sets + (Default value = :code:`False`). + Returns: + :math:`\left(N_{particles}\right)` :class:`numpy.ndarray`: + Vector of minimal RMSD values, one value per particle. + + """ + nq, nlist, qargs, l_query_points, num_query_points = \ + self._preprocess_arguments(system, neighbors=neighbors) + + motif = freud.util._convert_array(motif, shape=(None, 3)) + 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 self._cpp_obj.getRMSDs().toNumpyArray() + + 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' diff --git a/freud/environment/export-MatchEnv.cc b/freud/environment/export-MatchEnv.cc new file mode 100644 index 000000000..fb1927e54 --- /dev/null +++ b/freud/environment/export-MatchEnv.cc @@ -0,0 +1,65 @@ +// 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 "MatchEnv.h" +#include "Registration.h" + + +namespace nb = nanobind; + +namespace freud { namespace environment { + +template +using nb_array = nanobind::ndarray; + +namespace wrap { +void compute(const std::shared_ptr& match_env, + std::shared_ptr nq, + const nb_array>& query_points, + const unsigned int n_query_points, + const nb_array>& orientations, + 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()); + match_env->compute(nq, query_points_data, n_query_points, orientations_data, nlist, qargs, max_num_neighbors); + } + +}; + +namespace detail { + +void export_MatchEnv(nb::module_& module) +{ + // export minimizeRMSD function + // export isSimilar function + // export MatchEnv class + // export getPointEnvironments fn + nb::class_(module, "MatchEnv") + .def(nb::init<>) + .def("getPointEnvironments", &MatchEnv::getPointEnvironments) + // export EnvironmentCluster class + // export compute fn + // export getClusterIdx fn + // export getClusterEnvironments fn + // export getNumClusters fn + // export EnvironmentMotifMatch class + // export compute fn + // export getMatches fn + // export EnvironmentRMSDMinimizer class + // export compute fn + // export getRMSDs fn + +} + +}; }; // namespace detail +}; // namespace freud::locality diff --git a/freud/environment/module-environment.cc b/freud/environment/module-environment.cc index a2672dbc1..725cbf534 100644 --- a/freud/environment/module-environment.cc +++ b/freud/environment/module-environment.cc @@ -11,6 +11,7 @@ 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; @@ -21,4 +22,5 @@ NB_MODULE(_environment, module) // NOLINT(misc-use-anonymous-namespace): caused export_LocalBondProjection(module); export_LocalDescriptors(module); export_BondOrder(module); + export_MatchEnv(module); } From ccb80055b9002e5eb188f7c5022589f4e61a11cb Mon Sep 17 00:00:00 2001 From: Domagoj Fijan Date: Tue, 12 Nov 2024 17:54:57 -0500 Subject: [PATCH 23/28] update computes for env motif match and rmsd min with a wrapper --- freud/environment.py | 2 +- freud/environment/export-MatchEnv.cc | 70 ++++++++++++++++++---------- 2 files changed, 46 insertions(+), 26 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index b00a90a5a..af87f3bb4 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -775,7 +775,7 @@ def compute(self, system, threshold, cluster_neighbors=None, def cluster_idx(self): """:math:`\\left(N_{particles}\\right)` :class:`numpy.ndarray`: The per-particle index indicating cluster membership.""" - return self._cpp_obj.getClusterIdx().toNumpyArray() + return self._cpp_obj.getClusters().toNumpyArray() @_Compute._computed_property def num_clusters(self): diff --git a/freud/environment/export-MatchEnv.cc b/freud/environment/export-MatchEnv.cc index fb1927e54..a0507f22a 100644 --- a/freud/environment/export-MatchEnv.cc +++ b/freud/environment/export-MatchEnv.cc @@ -19,20 +19,33 @@ template using nb_array = nanobind::ndarray; namespace wrap { -void compute(const std::shared_ptr& match_env, +void compute_env_motif_match(const std::shared_ptr& env_motif_match, std::shared_ptr nq, - const nb_array>& query_points, - const unsigned int n_query_points, - const nb_array>& orientations, std::shared_ptr nlist, const locality::QueryArgs& qargs, - const unsigned int max_num_neighbors + const nb_array>& motif, + const unsigned int motif_size, + const float threshold, + const bool registration ) - { - auto* query_points_data = reinterpret_cast*>(query_points.data()); - auto* orientations_data = reinterpret_cast*>(orientations.data()); - match_env->compute(nq, query_points_data, n_query_points, orientations_data, nlist, qargs, max_num_neighbors); - } +{ + 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, + std::shared_ptr nq, + 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_rmsd_min->compute(nq, nlist, qargs, motif_data, motif_size, threshold, registration); +} }; @@ -40,24 +53,31 @@ namespace detail { void export_MatchEnv(nb::module_& module) { - // export minimizeRMSD function - // export isSimilar function - // export MatchEnv class - // export getPointEnvironments fn + // export minimizeRMSD function, move convenience to wrap? TODO + nb::function("minimizeRMSD", &minimizeRMSD);//carefull ths fn is overloaded for easier python interactivity. You should use the one that takes box etc in. + // export isSimilar function, move convenience to wrap? TODO + nb::function("isSimilar", &isSimilar); //carefull ths fn is overloaded for easier python interactivity. You should use the one that takes box etc in. + nb::class_(module, "MatchEnv") .def(nb::init<>) .def("getPointEnvironments", &MatchEnv::getPointEnvironments) - // export EnvironmentCluster class - // export compute fn - // export getClusterIdx fn - // export getClusterEnvironments fn - // export getNumClusters fn - // export EnvironmentMotifMatch class - // export compute fn - // export getMatches fn - // export EnvironmentRMSDMinimizer class - // export compute fn - // export getRMSDs fn + + nb::class_(module, "EnvironmentCluster") + .def(nb::init<>()) + .def("compute", &EnvironmentCluster::compute) + .def("getClusters", &EnvironmentCluster::getClusterIdx) + .def("getClusterEnvironments", &EnvironmentCluster::getClusterEnvironments) + .def("getNumClusters", &EnvironmentCluster::getNumClusters) + + nb::class_(module, "EnvironmentMotifMatch") + .def(nb::init<>()) + .def("compute", &wrap::compute_env_motif_match, nb::arg("nq"), nb::arg("nlist"), nb::arg("qargs"), nb::arg("motif"), nb::arg("motif_size"), nb::arg("threshold"), nb::arg("registration")) + .def("getMatches", &EnvironmentMotifMatch::getMatches) + + nb::class_(module, "EnvironmentRMSDMinimizer") + .def(nb::init<>()) + .def("compute", &wrap::compute_env_rmsd_min, nb::arg("nq"), nb::arg("nlist"), nb::arg("qargs"), nb::arg("motif"), nb::arg("motif_size"), nb::arg("threshold"), nb::arg("registration")) + .def("getRMSDs", &EnvironmentRMSDMinimizer::getRMSDs) } From c15a6349738a901bbc19a2a1d2f8d6d2b10f4a3f Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 15 Nov 2024 14:23:05 -0500 Subject: [PATCH 24/28] Minor BondOrder fixes --- freud/environment.py | 2 +- freud/environment/BondOrder.cc | 18 +++++++++--------- freud/environment/BondOrder.h | 1 + freud/environment/export-BondOrder.cc | 2 ++ 4 files changed, 13 insertions(+), 10 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index 6f709d579..263663420 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -347,7 +347,7 @@ def bond_order(self): # return freud.util.make_managed_numpy_array( # &self.thisptr.getBondOrder(), # freud.util.arr_type_t.FLOAT) - return self._cpp_obj.getBondOrder()#.toNumpyArray() + return self._cpp_obj.getBondOrder().toNumpyArray() @_Compute._computed_property def box(self): diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index 11bc50712..002a2d283 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -1,6 +1,6 @@ // 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 #ifdef __SSE2__ @@ -61,20 +61,20 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd std::make_shared(n_bins_theta, 0, constants::TWO_PI), std::make_shared(n_bins_phi, 0, M_PI) }; - // m_histogram = BondHistogram(axes); + m_histogram = BondHistogram(axes); m_local_histograms = BondHistogram::ThreadLocalHistogram(m_histogram); + } -void BondOrder::reduce() -{ - // TODO: previously, we could prepare the histogram: but not anymore? - // m_histogram.prepare(m_histogram.shape()); - // m_histogram = std::make_shared>(std::vector {m_histogram.shape()}); - - m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); +void BondOrder::reset() { + 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); }); diff --git a/freud/environment/BondOrder.h b/freud/environment/BondOrder.h index 060d4cd43..d5cb97c4c 100644 --- a/freud/environment/BondOrder.h +++ b/freud/environment/BondOrder.h @@ -45,6 +45,7 @@ class BondOrder : public locality::BondHistogramCompute const std::shared_ptr& nlist, freud::locality::QueryArgs qargs); void reduce() override; + void reset() override; //! Get a shared_ptr to the last computed bond order const std::shared_ptr> getBondOrder(); diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc index 31fba4b50..5f789488d 100644 --- a/freud/environment/export-BondOrder.cc +++ b/freud/environment/export-BondOrder.cc @@ -68,6 +68,8 @@ void export_BondOrder(nb::module_& module) .def(nb::init()) .def("getBondOrder", &BondOrder::getBondOrder) .def("getBinCounts", &BondOrder::getBinCounts) + // .def("getBinCenters", &BondOrder::getBinCenters) + .def("getBox", &BondOrder::getBox) .def("getAxisSizes", &BondOrder::getAxisSizes) .def("getMode", &BondOrder::getMode) .def("accumulate", &wrap::accumulateBondOrder, From a1ed072721775c54024ff0f97fdeef8ba81961db Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 15 Nov 2024 14:32:02 -0500 Subject: [PATCH 25/28] Clean up working functions --- freud/environment.py | 9 +-------- 1 file changed, 1 insertion(+), 8 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index 263663420..08c03330e 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -344,9 +344,6 @@ def compute(self, system, orientations=None, query_points=None, @_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 @@ -363,11 +360,7 @@ def __repr__(self): @property def mode(self): """str: Bond order mode.""" - mode = self._cpp_obj.getMode() - for key, value in self.known_modes.items(): - if value == mode: - return key - return mode + return self._cpp_obj.getMode().name class LocalDescriptors(_PairCompute): From 373139aa92f72361277ec6db37409ad0e6335611 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 15 Nov 2024 16:45:34 -0500 Subject: [PATCH 26/28] Fix BondOrder --- freud/environment/BondOrder.cc | 7 ++++++- freud/environment/BondOrder.h | 1 + freud/environment/export-BondOrder.cc | 5 +++-- freud/locality.py | 7 ++++--- tests/test_environment_bond_order.py | 6 ++++-- 5 files changed, 18 insertions(+), 8 deletions(-) diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index 002a2d283..c5abc942e 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -1,6 +1,5 @@ // 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 #ifdef __SSE2__ @@ -62,6 +61,7 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd 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); @@ -80,6 +80,11 @@ void BondOrder::reduce() }); } +std::vector> BondOrder::getBinCenters() +{ + return m_histogram.getBinCenters(); +} + const std::shared_ptr> BondOrder::getBondOrder() { return reduceAndReturn(m_bo_array); diff --git a/freud/environment/BondOrder.h b/freud/environment/BondOrder.h index d5cb97c4c..7b3628e44 100644 --- a/freud/environment/BondOrder.h +++ b/freud/environment/BondOrder.h @@ -46,6 +46,7 @@ class BondOrder : public locality::BondHistogramCompute void reduce() override; void reset() override; + std::vector> getBinCenters(); //! Get a shared_ptr to the last computed bond order const std::shared_ptr> getBondOrder(); diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc index 5f789488d..d4ed0f1f8 100644 --- a/freud/environment/export-BondOrder.cc +++ b/freud/environment/export-BondOrder.cc @@ -1,7 +1,6 @@ // 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 @@ -50,6 +49,7 @@ void accumulateBondOrder(const std::shared_ptr& self, ); } + }; // namespace wrap namespace detail { @@ -68,7 +68,8 @@ void export_BondOrder(nb::module_& module) .def(nb::init()) .def("getBondOrder", &BondOrder::getBondOrder) .def("getBinCounts", &BondOrder::getBinCounts) - // .def("getBinCenters", &BondOrder::getBinCenters) + .def("getBinCenters", &BondOrder::getBinCenters) + .def("getBinEdges", &BondOrder::getBinEdges) .def("getBox", &BondOrder::getBox) .def("getAxisSizes", &BondOrder::getAxisSizes) .def("getMode", &BondOrder::getMode) diff --git a/freud/locality.py b/freud/locality.py index 06114ade4..fa6008bdb 100644 --- a/freud/locality.py +++ b/freud/locality.py @@ -992,17 +992,18 @@ 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 + @property # TODO: Actual shape is 2, (d1=1, d2+1) and type is list def bin_centers(self): """:class:`numpy.ndarray`: The centers of each bin in the histogram (has the same shape as the histogram itself).""" centers = self._cpp_obj.getBinCenters() + # TODO: update docs to list the correct dimensions and type return [np.array(c) for c in centers] - @property + @property # TODO: Actual shape is 2, (d1=1, d2+1) and type is list 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 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") From e472950c0bbfe0b7073d8cde1b2591285948ed6f Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 15 Nov 2024 16:47:03 -0500 Subject: [PATCH 27/28] pre-commit working --- freud/environment/BondOrder.cc | 19 +++++------ freud/environment/BondOrder.h | 5 +-- freud/environment/export-BondOrder.cc | 49 ++++++++++----------------- freud/locality.py | 4 +-- 4 files changed, 32 insertions(+), 45 deletions(-) diff --git a/freud/environment/BondOrder.cc b/freud/environment/BondOrder.cc index c5abc942e..3c838ff66 100644 --- a/freud/environment/BondOrder.cc +++ b/freud/environment/BondOrder.cc @@ -1,5 +1,6 @@ // 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 #ifdef __SSE2__ @@ -56,19 +57,16 @@ BondOrder::BondOrder(unsigned int n_bins_theta, unsigned int n_bins_phi, BondOrd (*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) - }; + 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::reset() { +void BondOrder::reset() +{ BondHistogramCompute::reset(); m_bo_array = std::make_shared>(std::vector {m_histogram.shape()}); } @@ -90,9 +88,10 @@ const std::shared_ptr> BondOrder::getBondOrder() return reduceAndReturn(m_bo_array); } -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, +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, diff --git a/freud/environment/BondOrder.h b/freud/environment/BondOrder.h index 7b3628e44..6db775a9a 100644 --- a/freud/environment/BondOrder.h +++ b/freud/environment/BondOrder.h @@ -42,7 +42,8 @@ class BondOrder : public locality::BondHistogramCompute //! Accumulate the bond order void 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); + const std::shared_ptr& nlist, + freud::locality::QueryArgs qargs); void reduce() override; void reset() override; @@ -59,7 +60,7 @@ class BondOrder : public locality::BondHistogramCompute private: 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. + BondOrderMode m_mode; //!< The mode to calculate with. }; }; }; // end namespace freud::environment diff --git a/freud/environment/export-BondOrder.cc b/freud/environment/export-BondOrder.cc index d4ed0f1f8..b566b19b2 100644 --- a/freud/environment/export-BondOrder.cc +++ b/freud/environment/export-BondOrder.cc @@ -14,20 +14,16 @@ namespace nb = nanobind; namespace freud { namespace environment { -template -using nb_array = nb::ndarray; +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, - std::shared_ptr nlist, - const locality::QueryArgs& qargs - ) + const std::shared_ptr& nq, + const nb_array>& orientations, + const nb_array>& query_points, + const nb_array>& query_orientations, + 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; @@ -38,32 +34,28 @@ void accumulateBondOrder(const std::shared_ptr& self, // 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_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 - ); + 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) @@ -73,19 +65,14 @@ void export_BondOrder(nb::module_& module) .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) - ; + .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::locality +}; }; // namespace freud::environment diff --git a/freud/locality.py b/freud/locality.py index fa6008bdb..a66b61ef7 100644 --- a/freud/locality.py +++ b/freud/locality.py @@ -995,7 +995,7 @@ def bin_counts(self): """: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 # TODO: Actual shape is 2, (d1=1, d2+1) and type is list + @property # TODO: Actual shape is 2, (d1=1, d2+1) and type is list def bin_centers(self): """:class:`numpy.ndarray`: The centers of each bin in the histogram (has the same shape as the histogram itself).""" @@ -1003,7 +1003,7 @@ def bin_centers(self): # TODO: update docs to list the correct dimensions and type return [np.array(c) for c in centers] - @property # TODO: Actual shape is 2, (d1=1, d2+1) and type is list + @property # TODO: Actual shape is 2, (d1=1, d2+1) and type is list 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 From b7558970c6f33195a2ebb5edd2bebc2d637eddf8 Mon Sep 17 00:00:00 2001 From: janbridley Date: Fri, 15 Nov 2024 16:49:34 -0500 Subject: [PATCH 28/28] Lint BondOrder part of environment.py --- freud/environment.py | 49 +++++++++++++++++++++++++++++--------------- 1 file changed, 32 insertions(+), 17 deletions(-) diff --git a/freud/environment.py b/freud/environment.py index 08c03330e..ff38ba115 100644 --- a/freud/environment.py +++ b/freud/environment.py @@ -254,10 +254,12 @@ class BondOrder(_SpatialHistogram): (Default value = :code:`'bod'`). """ - 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 __init__(self, bins, mode="bod"): if isinstance(bins, collections.abc.Sequence): @@ -273,7 +275,6 @@ def __init__(self, bins, mode="bod"): self._cpp_obj = freud._environment.BondOrder(n_bins_theta, n_bins_phi, l_mode) - @property def default_query_args(self): """No default query arguments.""" @@ -282,8 +283,15 @@ def default_query_args(self): 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): + 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. @@ -317,8 +325,9 @@ def compute(self, system, orientations=None, query_points=None, if reset: self._reset() - 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: @@ -326,10 +335,15 @@ def compute(self, system, orientations=None, query_points=None, if query_points is None: query_points = nq.points - orientations = freud.util._convert_array(orientations, shape=(nq.points.shape[0], 4)) - query_orientations = freud.util._convert_array(query_orientations, shape=(num_query_points, 4)) - query_points = freud.util._convert_array(query_points, shape=(num_query_points, 3)) - + orientations = freud.util._convert_array( + orientations, shape=(nq.points.shape[0], 4) + ) + query_orientations = freud.util._convert_array( + 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, @@ -337,7 +351,7 @@ def compute(self, system, orientations=None, query_points=None, query_points, query_orientations, nlist._cpp_obj, - qargs._cpp_obj + qargs._cpp_obj, ) return self @@ -352,10 +366,11 @@ def box(self): 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):