Skip to content

Commit

Permalink
Finish cubatic
Browse files Browse the repository at this point in the history
  • Loading branch information
janbridley committed Oct 26, 2024
1 parent 2cc36c2 commit 3d049a3
Show file tree
Hide file tree
Showing 4 changed files with 91 additions and 99 deletions.
139 changes: 62 additions & 77 deletions freud/order.py
Original file line number Diff line number Diff line change
Expand Up @@ -66,94 +66,79 @@ def __init__(self, t_initial, t_final, scale, n_replicates=1, seed=None):

self._cpp_obj = freud._order.Cubatic(t_initial, t_final, scale, n_replicates, seed)

# def compute(self, orientations):
# r"""Calculates the per-particle and global order parameter.

# Args:
# orientations ((:math:`N_{particles}`, 4) :class:`numpy.ndarray`):
# Orientations as angles to use in computation.
# """
# orientations = freud.util._convert_array(
# orientations, shape=(None, 4))

# cdef const float[:, ::1] l_orientations = orientations
# cdef unsigned int num_particles = l_orientations.shape[0]
def compute(self, orientations):
r"""Calculates the per-particle and global order parameter.
# self.thisptr.compute(
# <quat[float]*> &l_orientations[0, 0], num_particles)
# return self
Args:
orientations ((:math:`N_{particles}`, 4) :class:`numpy.ndarray`):
Orientations as angles to use in computation.
"""
orientations = freud.util._convert_array(orientations, shape=(None, 4))
self._cpp_obj.compute(orientations)
return self

# @property
# def t_initial(self):
# """float: The value of the initial temperature."""
# return self.thisptr.getTInitial()
@property
def t_initial(self):
"""float: The value of the initial temperature."""
return self._cpp_obj.getTInitial()

# @property
# def t_final(self):
# """float: The value of the final temperature."""
# return self.thisptr.getTFinal()
@property
def t_final(self):
"""float: The value of the final temperature."""
return self._cpp_obj.getTFinal()

# @property
# def scale(self):
# """float: The scale."""
# return self.thisptr.getScale()
@property
def scale(self):
"""float: The scale."""
return self._cpp_obj.getScale()

# @property
# def n_replicates(self):
# """unsigned int: Number of replicate simulated annealing runs."""
# return self.thisptr.getNReplicates()
@property
def n_replicates(self):
"""unsigned int: Number of replicate simulated annealing runs."""
return self._cpp_obj.getNReplicates()

# @property
# def seed(self):
# """unsigned int: Random seed to use in calculations."""
# return self.thisptr.getSeed()
@property
def seed(self):
"""unsigned int: Random seed to use in calculations."""
return self._cpp_obj.getSeed()

# @_Compute._computed_property
# def order(self):
# """float: Cubatic order parameter of the system."""
# return self.thisptr.getCubaticOrderParameter()
@_Compute._computed_property
def order(self):
"""float: Cubatic order parameter of the system."""
return self._cpp_obj.getCubaticOrderParameter()

# @_Compute._computed_property
# def orientation(self):
# """:math:`\\left(4 \\right)` :class:`numpy.ndarray`: The quaternion of
# global orientation."""
# cdef quat[float] q = self.thisptr.getCubaticOrientation()
# return np.asarray([q.s, q.v.x, q.v.y, q.v.z], dtype=np.float32)
@_Compute._computed_property
def orientation(self):
""":math:`\\left(4 \\right)` :class:`numpy.ndarray`: The quaternion of
global orientation."""
q = self._cpp_obj.getCubaticOrientation()
return np.asarray(q, dtype=np.float32)

# @_Compute._computed_property
# def particle_order(self):
# """:math:`\\left(N_{particles} \\right)` :class:`numpy.ndarray`: Order
# parameter."""
# return freud.util.make_managed_numpy_array(
# &self.thisptr.getParticleOrderParameter(),
# freud.util.arr_type_t.FLOAT)
@_Compute._computed_property
def particle_order(self):
""":math:`\\left(N_{particles} \\right)` :class:`numpy.ndarray`: Order
parameter."""
return self._cpp_obj.getParticleOrderParameter().toNumpyArray()

# @_Compute._computed_property
# def global_tensor(self):
# """:math:`\\left(3, 3, 3, 3 \\right)` :class:`numpy.ndarray`: Rank 4
# tensor corresponding to the global orientation. Computed from all
# orientations."""
# return freud.util.make_managed_numpy_array(
# &self.thisptr.getGlobalTensor(),
# freud.util.arr_type_t.FLOAT)
@_Compute._computed_property
def global_tensor(self):
""":math:`\\left(3, 3, 3, 3 \\right)` :class:`numpy.ndarray`: Rank 4
tensor corresponding to the global orientation. Computed from all
orientations."""
return self._cpp_obj.getGlobalTensor().toNumpyArray()

# @_Compute._computed_property
# def cubatic_tensor(self):
# """:math:`\\left(3, 3, 3, 3 \\right)` :class:`numpy.ndarray`: Rank 4
# homogeneous tensor representing the optimal system-wide coordinates."""
# return freud.util.make_managed_numpy_array(
# &self.thisptr.getCubaticTensor(),
# freud.util.arr_type_t.FLOAT)
@_Compute._computed_property
def cubatic_tensor(self):
""":math:`\\left(3, 3, 3, 3 \\right)` :class:`numpy.ndarray`: Rank 4
homogeneous tensor representing the optimal system-wide coordinates."""
return self._cpp_obj.getCubaticTensor().toNumpyArray()

# def __repr__(self):
# return ("freud.order.{cls}(t_initial={t_initial}, t_final={t_final}, "
# "scale={scale}, n_replicates={n_replicates}, "
# "seed={seed})").format(cls=type(self).__name__,
# t_initial=self.t_initial,
# t_final=self.t_final,
# scale=self.scale,
# n_replicates=self.n_replicates,
# seed=self.seed)
def __repr__(self):
return (
f"freud.order.{type(self).__name__}(t_initial={self.t_initial}, "
f"t_final={self.t_final}, scale={self.scale}, "
f"n_replicates={self.n_replicates}, seed={self.seed})"
)


class Nematic(_Compute):
Expand Down Expand Up @@ -565,7 +550,7 @@ def particle_order(self):
Variant of the Steinhardt order parameter for each particle (filled with
:code:`nan` for particles with no neighbors)."""
p_orders = self._cpp_obj.getParticleOrder().toNumpyArray()
return np.ravel(p_orders) if p_orders.shape[1] == 1 else p_orders
return np.ravel(p_orders) if p_orders.shape[1] == 1 else p_orders

@_Compute._computed_property
def ql(self):
Expand Down
15 changes: 8 additions & 7 deletions freud/order/Cubatic.cc
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,11 @@ tensor4 tensor4::operator*(const float& b) const

void tensor4::copyToManagedArray(util::ManagedArray<float>& ma)
{
std::copy(data.begin(), data.end(), ma.get());
// TODO: this may be possible with std::copy but this works as well.
for (unsigned int i = 0; i < 81; i++)
{
ma[i] = data[i];
}
}

//! Complete tensor contraction.
Expand Down Expand Up @@ -236,14 +240,12 @@ tensor4 Cubatic::calculateGlobalTensor(quat<float>* orientations) const
void Cubatic::compute(quat<float>* orientations, unsigned int num_orientations)
{
m_n = num_orientations;
// m_particle_order_parameter.prepare(m_n);
m_particle_order_parameter = std::make_shared<util::ManagedArray<float>>(std::vector<size_t> {m_n});

// Calculate the per-particle tensor
tensor4 global_tensor = calculateGlobalTensor(orientations);
// m_global_tensor.prepare({3, 3, 3, 3});
m_global_tensor = std::make_shared<util::ManagedArray<float>>(std::vector<size_t> {3, 3, 3, 3});
global_tensor.copyToManagedArray(m_global_tensor);
global_tensor.copyToManagedArray((*m_global_tensor));

// The paper recommends using a Newton-Raphson scheme to optimize the order
// parameter, but in practice we find that simulated annealing performs
Expand Down Expand Up @@ -331,9 +333,8 @@ void Cubatic::compute(quat<float>* orientations, unsigned int num_orientations)
}
}

// m_cubatic_tensor.prepare({3, 3, 3, 3});
m_cubatic_tensor = std::make_shared<util::ManagedArray<float>>(std::vector<size_t> {3, 3, 3, 3});
p_cubatic_tensor[max_idx].copyToManagedArray(m_cubatic_tensor);
p_cubatic_tensor[max_idx].copyToManagedArray((*m_cubatic_tensor));
m_cubatic_orientation = p_cubatic_orientation[max_idx];
m_cubatic_order_parameter = p_cubatic_order_parameter[max_idx];

Expand All @@ -344,7 +345,7 @@ void Cubatic::compute(quat<float>* orientations, unsigned int num_orientations)
// The per-particle order parameter is defined as the value of the
// cubatic order parameter if the global orientation was the
// particle orientation, so we can reuse the same machinery.
m_particle_order_parameter[i]
(*m_particle_order_parameter)[i]
= calcCubaticOrderParameter(calcCubaticTensor(orientations[i]), global_tensor);
}
});
Expand Down
6 changes: 3 additions & 3 deletions freud/order/Cubatic.h
Original file line number Diff line number Diff line change
Expand Up @@ -78,17 +78,17 @@ class Cubatic
return m_cubatic_order_parameter;
}

const util::ManagedArray<float>& getParticleOrderParameter() const
const std::shared_ptr<util::ManagedArray<float>>& getParticleOrderParameter() const
{
return m_particle_order_parameter;
}

const util::ManagedArray<float>& getGlobalTensor() const
const std::shared_ptr<util::ManagedArray<float>>& getGlobalTensor() const
{
return m_global_tensor;
}

const util::ManagedArray<float>& getCubaticTensor() const
const std::shared_ptr<util::ManagedArray<float>>& getCubaticTensor() const
{
return m_cubatic_tensor;
}
Expand Down
30 changes: 18 additions & 12 deletions freud/order/export-Cubatic.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,7 @@
#include <memory>
#include <nanobind/nanobind.h>
#include <nanobind/ndarray.h>
// #include <nanobind/stl/shared_ptr.h> // NOLINT(misc-include-cleaner): used implicitly
// #include <nanobind/stl/tuple.h> // NOLINT(misc-include-cleaner): used implicitly
#include <nanobind/stl/shared_ptr.h> // NOLINT(misc-include-cleaner): used implicitly
#include <utility>

#include "Cubatic.h"
Expand All @@ -27,24 +26,31 @@ void computeCubatic(const std::shared_ptr<Cubatic>& self,
self->compute(orientations_data, num_orientations);
}

// nanobind::tuple getNematicDirector(const std::shared_ptr<Cubatic>& self)
// {
// vec3<float> nem_d_cpp = self->getNematicDirector();
// return nanobind::make_tuple(nem_d_cpp.x, nem_d_cpp.y, nem_d_cpp.z);
// }
nanobind::tuple getCubaticOrientation(const std::shared_ptr<Cubatic>& self)
{
quat<float> q = self->getCubaticOrientation();
return nanobind::make_tuple(q.s, q.v.x, q.v.y, q.v.z);
}
}; // namespace wrap

namespace detail {

void export_Nematic(nanobind::module_& m)
void export_Cubatic(nanobind::module_& m)
{
nanobind::class_<Cubatic>(m, "Cubatic")
.def(nanobind::init<float, float, float, unsigned int, unsigned int>())
.def("compute", &wrap::computeCubatic, nanobind::arg("orientations"))
// .def("getNematicOrderParameter", &Nematic::getNematicOrderParameter)
// .def("getNematicDirector", &wrap::getNematicDirector)
// .def("getParticleTensor", &Nematic::getParticleTensor)
// .def("getNematicTensor", &Nematic::getNematicTensor);
.def("getTInitial", &Cubatic::getTInitial)
.def("getTFinal", &Cubatic::getTFinal)
.def("getScale", &Cubatic::getScale)
.def("getNReplicates", &Cubatic::getNReplicates)
.def("getSeed", &Cubatic::getSeed)
.def("getCubaticOrderParameter", &Cubatic::getCubaticOrderParameter)
.def("getCubaticOrientation", &wrap::getCubaticOrientation)
.def("getParticleOrderParameter", &Cubatic::getParticleOrderParameter)
.def("getGlobalTensor", &Cubatic::getGlobalTensor)
.def("getCubaticTensor", &Cubatic::getCubaticTensor)
;
}

} // namespace detail
Expand Down

0 comments on commit 3d049a3

Please sign in to comment.