Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added tests and documentation for cases zero neighbor cases. #775

Draft
wants to merge 12 commits into
base: main
Choose a base branch
from
9 changes: 8 additions & 1 deletion cpp/order/HexaticTranslational.cc
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,14 @@ void HexaticTranslational<T>::computeGeneral(Func func, const freud::locality::N
}
if (normalize_by_k)
{
m_psi_array[i] /= std::complex<float>(m_k);
if (total_weight == 0.)
{
m_psi_array[i] /= std::complex<float>(total_weight);
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Force this to be nan by dividing by zero. Is there a better way to do this? I saw std::nanf but I wasn't sure what to pick as the argument.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I haven't reviewed this PR thoroughly yet, but I'm not sure this is generating the output you want. std::complex<float>(1)/std::complex<float>(0) will return (inf, nan), not just nan (here's an example)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

So tried to set m_psy_array equal to just nan like this:

m_psi_array[i] = std::complex<float>(std::nanf)

and then like this

m_psi_array[i] = std::complex<float>(std::nanf)

All tests pass on my laptop, but fail on CI, so I reverted back to a divide by zero. I double checked the values from the no neighbor test for the Translational order parameter, and in python, it registers as nan and not inf.

}
else
{
m_psi_array[i] /= std::complex<float>(m_k);
}
}
else
{
Expand Down
33 changes: 31 additions & 2 deletions freud/order.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -283,6 +283,15 @@ cdef class Hexatic(_PairCompute):
**2D:** :class:`freud.order.Hexatic` is only defined for 2D systems.
The points must be passed in as :code:`[x, y, 0]`.

.. note::
The value of per-particle order parameter will be set to NaN for
particles with no neighbors. We choose this value rather than setting
the order parameter to 0 because in more complex order parameter
calculations, it is possible to observe a value of 0 for the
per-particle order parameter even with a finite number of neighbors.
If you would like to ignore this distinction, you can mask the output
order parameter values using NumPy: :code:`numpy.nan_to_num(particle_order)`.

Args:
k (unsigned int, optional):
Symmetry of order parameter (Default value = :code:`6`).
Expand Down Expand Up @@ -420,6 +429,16 @@ cdef class Translational(_PairCompute):
.. note::
This class is slated for deprecation and will be removed in freud 3.0.

.. note::
The value of per-particle order parameter will be set to NaN for
particles with no neighbors. We choose this value rather than setting
the order parameter to 0 because in more complex order parameter
calculations, it is possible to observe a value of 0 for the
per-particle order parameter even with a finite number of neighbors.
If you would like to ignore this distinction, you can mask the output
order parameter values using NumPy: :code:`numpy.nan_to_num(particle_order)`.


Args:
k (float, optional):
Normalization of order parameter (Default value = :code:`6.0`).
Expand Down Expand Up @@ -628,7 +647,7 @@ cdef class Steinhardt(_PairCompute):
@_Compute._computed_property
def particle_harmonics(self):
""":math:`\\left(N_{particles}, 2*l+1\\right)` :class:`numpy.ndarray`:
The raw array of \\overline{q}_{lm}(i). The array is provided in the
The raw array of :math:`\\overline{q}_{lm}(i)`. The array is provided in the
order given by fsph: :math:`m = 0, 1, ..., l, -1, ..., -l`."""
return freud.util.make_managed_numpy_array(
&self.thisptr.getQlm(),
Expand Down Expand Up @@ -735,6 +754,15 @@ cdef class SolidLiquid(_PairCompute):
the particle is considered solid-like. Finally, solid-like particles are
clustered.

.. note::
The value of :math:`q_l(i, j)` will be set to NaN for
particles with no neighbors. We choose this value rather than setting
the order parameter to 0 because in more complex order parameter
calculations, it is possible to observe a value of 0 for the :math:`q_l(i, j)`
order parameter even with a finite number of neighbors.
If you would like to ignore this distinction, you can mask the output order
parameter values using NumPy: :code:`numpy.nan_to_num(particle_order)`.

Args:
l (unsigned int):
Spherical harmonic quantum number l.
Expand Down Expand Up @@ -786,6 +814,7 @@ cdef class SolidLiquid(_PairCompute):
self.thisptr.compute(nlist.get_ptr(),
nq.get_ptr(),
dereference(qargs.thisptr))
return self

@property
def l(self): # noqa: E743
Expand Down Expand Up @@ -827,7 +856,7 @@ cdef class SolidLiquid(_PairCompute):
@_Compute._computed_property
def particle_harmonics(self):
""":math:`\\left(N_{particles}, 2*l+1\\right)` :class:`numpy.ndarray`:
The raw array of \\overline{q}_{lm}(i). The array is provided in the
The raw array of :math:`\\overline{q}_{lm}(i)`. The array is provided in the
order given by fsph: :math:`m = 0, 1, ..., l, -1, ..., -l`."""
return freud.util.make_managed_numpy_array(
&self.thisptr.getQlm(),
Expand Down
10 changes: 10 additions & 0 deletions tests/test_environment_LocalDescriptors.py
Original file line number Diff line number Diff line change
Expand Up @@ -438,3 +438,13 @@ def test_query_point_ne_points(self):
"Failed for l={}, m={}, x={}, y = {}" "\ntheta={}, phi={}"
).format(l, m, scipy_val, ld_val, theta, phi)
count += 1

def test_no_neighbors(self):
l_max = 8
box = freud.box.Box.cube(10)
points = [(0, 0, 0)]

ld = freud.environment.LocalDescriptors(l_max)
ld.compute((box, points), neighbors={"r_max": 1.25})
assert ld.num_sphs == 0
assert len(ld.sph) == 0
10 changes: 10 additions & 0 deletions tests/test_order_Hexatic.py
Original file line number Diff line number Diff line change
Expand Up @@ -218,3 +218,13 @@ def test_repr_png(self):
hop._repr_png_()
hop.plot()
plt.close("all")

def test_no_neighbors(self):
"""Ensure that particles without neighbors are assigned NaN"""
box = freud.box.Box.square(10)
positions = [(0, 0, 0)]
hop = freud.order.Hexatic()
hop.compute((box, positions), neighbors={"r_max": 1.25})

assert np.all(np.isnan(hop.particle_order))
npt.assert_allclose(np.nan_to_num(hop.particle_order), 0)
15 changes: 15 additions & 0 deletions tests/test_order_SolidLiquid.py
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
import matplotlib
import numpy as np
import numpy.testing as npt
import pytest

Expand Down Expand Up @@ -95,3 +96,17 @@ def test_attribute_access(self):
def test_repr(self):
comp = freud.order.SolidLiquid(6, q_threshold=0.7, solid_threshold=6)
assert str(comp) == str(eval(repr(comp)))

def test_no_neighbors(self):
"""Ensure that particles without neighbors are assigned NaN"""
box = freud.box.Box.cube(10)
positions = [(0, 0, 0)]
comp = freud.order.SolidLiquid(6, q_threshold=0.7, solid_threshold=6)
comp.compute((box, positions), neighbors=dict(r_max=2.0))

npt.assert_equal(comp.cluster_idx, np.arange(len(comp.cluster_idx)))
npt.assert_equal(comp.cluster_sizes, np.ones_like(comp.cluster_sizes))
assert comp.largest_cluster_size == 1
npt.assert_equal(comp.num_connections, np.zeros_like(comp.cluster_sizes))
assert np.all(np.isnan(comp.particle_harmonics))
npt.assert_equal(comp.ql_ij, [])
6 changes: 6 additions & 0 deletions tests/test_order_Translational.py
Original file line number Diff line number Diff line change
Expand Up @@ -40,3 +40,9 @@ def test_repr(self):
with pytest.warns(FreudDeprecationWarning):
trans = freud.order.Translational(4)
assert str(trans) == str(eval(repr(trans)))

def test_no_neighbors(self):
box = freud.box.Box.square(10)
positions = [(0, 0, 0)]
trans = freud.order.Translational(4)
trans.compute((box, positions), neighbors={"r_max": 1.25})