diff --git a/stumpy/core.py b/stumpy/core.py index 09fd80d8f..0f85effce 100644 --- a/stumpy/core.py +++ b/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( diff --git a/tests/test_precision.py b/tests/test_precision.py index 043406bba..8338f2548 100644 --- a/tests/test_precision.py +++ b/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.