Improve MASS precision for near-perfect correlations

stumpy/core.py

@@ -1576,16 +1576,43 @@ def _mass(Q, T, QT, μ_Q, σ_Q, M_T, Σ_T, Q_subseq_isconstant, T_subseq_isconst
 
     See Table II
 
-    Note that Q, T are not directly required to calculate D
-
     Note: Unlike the Matrix Profile I paper, here, M_T, Σ_T can be calculated
     once for all subsequences of T and passed in so the redundancy is removed
     """
     m = Q.shape[0]
+    k = M_T.shape[0]
+    distance_profile = np.empty(k, dtype=np.float64)
+    # When `1 - ρ` approaches machine precision, relative error from cancellation
+    # dominates. Recompute those rare distances from pointwise z-normalized
+    # differences once `1 - ρ` falls below `sqrt(eps)`.
+    threshold = 2.0 * m * np.sqrt(np.finfo(np.float64).eps)
 
-    return calculate_distance_profile(
-        m, QT, μ_Q, σ_Q, M_T, Σ_T, Q_subseq_isconstant, T_subseq_isconstant
-    )
+    for i in range(k):
+        D_squared = _calculate_squared_distance(
+            m,
+            QT[i],
+            μ_Q,
+            σ_Q,
+            M_T[i],
+            Σ_T[i],
+            Q_subseq_isconstant,
+            T_subseq_isconstant[i],
+        )
+
+        if (
+            D_squared < threshold
+            and not Q_subseq_isconstant
+            and not T_subseq_isconstant[i]
+        ):
+            D_squared = 0.0
+            for j in range(m):
+                Q_norm = (Q[j] - μ_Q) / σ_Q
+                T_norm = (T[i + j] - M_T[i]) / Σ_T[i]
+                D_squared += (Q_norm - T_norm) ** 2
+
+        distance_profile[i] = np.sqrt(D_squared)
+
+    return distance_profile
 
 
 @non_normalized(

tests/test_precision.py

@@ -74,6 +74,41 @@ def test_distace_profile():
         npt.assert_almost_equal(D_ref, D_comp)
 
 
+def test_mass_near_perfect_correlation():
+    # Regression data from #1171 contains two nearly affine-equivalent subsequences.
+    Q = np.array(
+        [
+            0.93406220661832629,
+            0.94732544958491172,
+            0.94288413821602846,
+        ]
+    )
+    T = np.array(
+        [
+            0.25089093970715337,
+            0.44046885357649379,
+            0.37698714499372254,
+        ]
+    )
+    m = Q.shape[0]
+
+    ref = naive.distance_profile(Q, T, m)
+    comp = core.mass(Q, T)
+
+    npt.assert_allclose(ref, comp, rtol=1e-12, atol=1e-12)
+
+
+def test_mass_near_perfect_correlation_rounded_to_zero():
+    Q = np.array([1.0, 2.0, 4.0])
+    T = np.array([-4.0, -0.99999999, 5.0])
+    m = Q.shape[0]
+
+    ref = naive.distance_profile(Q, T, m)
+    comp = core.mass(Q, T)
+
+    npt.assert_allclose(ref, comp, rtol=1e-12, atol=1e-12)
+
+
 def test_calculate_squared_distance():
     # This test function raises an error if the distance between a subsequence
     # and another does not satisfy the symmetry property.