Preserve correlation for nonconstant small-angle samples

pycircstat2/correlation.py

@@ -107,6 +107,14 @@ def _coerce_angles(
     return array, None
 
 
+def _has_zero_pairwise_variation(angles: np.ndarray) -> bool:
+    """Detect constant or antipodal directions without an absolute variance floor."""
+    differences = angles - angles[0]
+    # sin(k*pi) has rounding residuals proportional to the argument's magnitude.
+    tolerance = np.finfo(float).eps * np.abs(differences)
+    return bool(np.all(np.abs(np.sin(differences)) <= tolerance))
+
+
 def _circ_corrcc_fl(
     a: Union[Circular, np.ndarray, Sequence[float]],
     b: Union[Circular, np.ndarray, Sequence[float]],
@@ -146,7 +154,7 @@ def _corr(a, b):
         diff_b = np.sin(np.subtract.outer(b, b)[np.triu_indices(len(b), k=1)])
         num = np.sum(diff_a * diff_b)
         den = np.sqrt(np.sum(diff_a**2) * np.sum(diff_b**2))
-        if np.isclose(den, 0.0):
+        if den == 0.0 or _has_zero_pairwise_variation(a) or _has_zero_pairwise_variation(b):
             raise ValueError("Degenerate data produced zero variance in pairwise differences.")
         return num / den
 
@@ -223,7 +231,7 @@ def _circ_corrcc_js(
     Sb = np.sin(bbar)
     num = np.sum(Sa * Sb)
     den = np.sqrt(np.sum(Sa**2) * np.sum(Sb**2))
-    if np.isclose(den, 0.0):
+    if den == 0.0 or _has_zero_pairwise_variation(a_alpha) or _has_zero_pairwise_variation(b_alpha):
         raise ValueError("Degenerate data produced zero variance around the mean direction.")
     r = num / den
 
@@ -232,7 +240,7 @@ def _circ_corrcc_js(
         l20 = np.mean(Sa**2)
         l02 = np.mean(Sb**2)
         l22 = np.mean(Sa**2 * Sb**2)
-        if np.isclose(l22, 0.0):
+        if l22 == 0.0:
             raise ValueError("Degenerate data caused division by zero in variance term.")
         test_stat = np.sqrt(n) * np.sqrt(l20 * l02 / l22) * r
         p_val = 2 * (1 - norm.cdf(np.abs(test_stat)))

tests/test_correlation_small_angles.py

@@ -0,0 +1,30 @@
+import numpy as np
+import pytest
+
+from pycircstat2.correlation import circ_corrcc
+
+
+@pytest.mark.parametrize("method", ["fl", "js"])
+@pytest.mark.parametrize("test", [False, True])
+@pytest.mark.parametrize("scale", [1e-3, 1e-5, 1e-20])
+def test_nonconstant_small_angles_have_unit_self_correlation(method, test, scale):
+    angles = scale * np.array([0.0, 1.0, 4.0, 2.0, 3.0])
+    result = circ_corrcc(angles, angles, method=method, test=test)
+    r = result.r if test else result
+    assert r == pytest.approx(1.0)
+    if method == "js" and test:
+        assert np.isfinite(result.test_stat)
+        assert 0.0 < result.p_value < 1.0
+
+
+@pytest.mark.parametrize("method", ["fl", "js"])
+@pytest.mark.parametrize("test", [False, True])
+@pytest.mark.parametrize("angles", [
+    np.zeros(5),
+    np.full(5, 0.7),
+    np.array([0.0, np.pi, 0.0, np.pi, 0.0]),
+    np.array([0.7, 0.7 + 2 * np.pi, 0.7, 0.7 + 2 * np.pi, 0.7]),
+])
+def test_constant_or_antipodal_angles_remain_degenerate(method, test, angles):
+    with pytest.raises(ValueError, match="zero variance"):
+        circ_corrcc(angles, angles, method=method, test=test)