Handle symmetry ops that produce fractional HKLs in apply_to_hkl()

reciprocalspaceship/utils/symop.py

@@ -14,10 +14,29 @@ def apply_to_hkl(H, op):
 
     Returns
     -------
-    result : array
+    result : np.ndarray(int32)
         n x 3 array of Miller indices after operator application
+
+    Raises
+    ------
+    RuntimeError
+        If `op` generates fractional Miller indices when applied to `H`
     """
-    return np.floor_divide(np.matmul(H, op.rot), op.DEN)
+    # Case 1: No risk of fractional Miller indices
+    if ((np.array(op.rot) / op.DEN) % 1 == 0).all():
+        return np.floor_divide(np.matmul(H, op.rot), op.DEN).astype("int32")
+
+    # Case 2: Depends on input Miller indices
+    else:
+        Hnew = np.divide(np.matmul(H, op.rot), op.DEN)
+        # Check for issues
+        if np.any(np.mod(Hnew, 1)):
+            raise RuntimeError(
+                f"Applying {op} to Miller indices produced non-integer results. "
+                f"Fractional Miller indices are not currently supported."
+            )
+        else:
+            return Hnew.astype("int32")
 
 
 def phase_shift(H, op):

tests/conftest.py

@@ -4,12 +4,21 @@
 
 import gemmi
 import numpy as np
-import pandas as pd
 import pytest
 
 import reciprocalspaceship as rs
 
 
+@pytest.fixture
+def hkls():
+    """
+    Return all Miller indices with H, K, L values between [-5, 5]
+    """
+    hmin, hmax = -5, 5
+    H = np.mgrid[hmin : hmax + 1, hmin : hmax + 1, hmin : hmax + 1].reshape((3, -1)).T
+    return H
+
+
 @pytest.fixture
 def dataset_hkl():
     """

tests/utils/test_symop.py

@@ -0,0 +1,32 @@
+import gemmi
+import numpy as np
+import pytest
+
+import reciprocalspaceship as rs
+
+
+@pytest.mark.parametrize("H_even", [True, False])
+@pytest.mark.parametrize(
+    "op_str", ["x,y,z", "2*x,2*y,2*z", "x,z,y", "1/2*x,y,z", "1/2*x,1/2*y,1/2*z"]
+)
+def test_apply_to_hkl(hkls, H_even, op_str):
+    """
+    Test rs.utils.apply_to_hkl() detects symops that yield fractional Miller indices.
+
+    apply_to_hkl() should raise a RuntimeError if the combination of `H` and `op`
+    yield fractional Miller indices, and should return new Miller indices all other
+    cases.
+    """
+    if H_even:
+        hkls = hkls[~np.any(hkls % 2, axis=1)]
+
+    op = gemmi.Op(op_str)
+
+    if ((np.array(op.rot) / op.DEN) % 1 == 0).all() or H_even:
+        H_result = rs.utils.apply_to_hkl(hkls, op)
+        H_expected = np.array([op.apply_to_hkl(hkl) for hkl in hkls])
+        assert np.array_equal(H_expected, H_result)
+        assert H_result.dtype is np.dtype(np.int32)
+    else:
+        with pytest.raises(RuntimeError):
+            H_result = rs.utils.apply_to_hkl(hkls, op)