Fixed #881 Snippets Unit Test Assertion Failure (#892)

* add test function to test_precision, expecting assertion failure

* add post processing to have symmetry property of distance without loss of precision

* retrieve code, expecting error

* reduce arithmetic operation to reduce loss of precision

* revise equation to mitigate loss of precision

* add test function to check precision for _calculate_squared_distance

* minor change

* change order of test function in test_precision

* add parentheses to make sure the ouput remains the same when values are swapped

* avoid fastmath reassoc flag to increase floating point precision

* enhance test function

* add fatmath flag reassoc and remove ninf

* minor change in test function

* increase decimal coverage in assertion

* remove nnan from fastmath to be consistent with other func in core.py

* use assert to make sure two values are identical

* re-order arithmetic operation to increase safety and stability of the code regarding precision

* add test function for distance symmetry propery in gpu

* add parentheses and change order of multiplication

* skip testing gpu-based test func if cude is not available

* perform assersion with numpy testing instead of assert to be more explicit

* reduce precision to pass tests for minimum version

* add comment to test function to explain its purpose

stumpy/core.py

@@ -1042,7 +1042,7 @@ def compute_mean_std(T, m):
 
 @njit(
     # "f8(i8, f8, f8, f8, f8, f8)",
-    fastmath=True
+    fastmath={"nsz", "arcp", "contract", "afn", "reassoc"}
 )
 def _calculate_squared_distance(
     m, QT, μ_Q, σ_Q, M_T, Σ_T, Q_subseq_isconstant, T_subseq_isconstant
@@ -1097,10 +1097,10 @@ def _calculate_squared_distance(
     elif Q_subseq_isconstant or T_subseq_isconstant:
         D_squared = m
     else:
-        denom = m * σ_Q * Σ_T
+        denom = (σ_Q * Σ_T) * m
         denom = max(denom, config.STUMPY_DENOM_THRESHOLD)
 
-        ρ = (QT - m * μ_Q * M_T) / denom
+        ρ = (QT - (μ_Q * M_T) * m) / denom
         ρ = min(ρ, 1.0)
 
         D_squared = np.abs(2 * m * (1.0 - ρ))

stumpy/gpu_stump.py

@@ -178,9 +178,9 @@ def _compute_and_update_PI_kernel(
         elif Q_subseq_isconstant[i] or T_subseq_isconstant[j]:
             p_norm = m
         else:
-            denom = m * σ_Q[i] * Σ_T[j]
+            denom = (σ_Q[i] * Σ_T[j]) * m
             denom = max(denom, config.STUMPY_DENOM_THRESHOLD)
-            ρ = (QT_out[i] - m * μ_Q[i] * M_T[j]) / denom
+            ρ = (QT_out[i] - (μ_Q[i] * M_T[j]) * m) / denom
             ρ = min(ρ, 1.0)
             p_norm = 2 * m * (1.0 - ρ)
 

tests/test_precision.py

@@ -1,10 +1,22 @@
+import functools
+from unittest.mock import patch
+
 import naive
 import numpy as np
 import numpy.testing as npt
+import pytest
+from numba import cuda
 
 import stumpy
 from stumpy import config, core
 
+try:
+    from numba.errors import NumbaPerformanceWarning
+except ModuleNotFoundError:
+    from numba.core.errors import NumbaPerformanceWarning
+
+TEST_THREADS_PER_BLOCK = 10
+
 
 def test_mpdist_snippets_s():
     # This test function raises an error if the distance between
@@ -55,3 +67,134 @@ def test_distace_profile():
         )
 
         npt.assert_almost_equal(D_ref, D_comp)
+
+
+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.
+    seed = 332
+    np.random.seed(seed)
+    T = np.random.uniform(-1000.0, 1000.0, [64])
+    m = 3
+
+    T_subseq_isconstant = core.rolling_isconstant(T, m)
+    M_T, Σ_T = core.compute_mean_std(T, m)
+
+    n = len(T)
+    k = n - m + 1
+    for i in range(k):
+        for j in range(k):
+            QT_i = core._sliding_dot_product(T[i : i + m], T)
+            dist_ij = core._calculate_squared_distance(
+                m,
+                QT_i[j],
+                M_T[i],
+                Σ_T[i],
+                M_T[j],
+                Σ_T[j],
+                T_subseq_isconstant[i],
+                T_subseq_isconstant[j],
+            )
+
+            QT_j = core._sliding_dot_product(T[j : j + m], T)
+            dist_ji = core._calculate_squared_distance(
+                m,
+                QT_j[i],
+                M_T[j],
+                Σ_T[j],
+                M_T[i],
+                Σ_T[i],
+                T_subseq_isconstant[j],
+                T_subseq_isconstant[i],
+            )
+
+            comp = dist_ij - dist_ji
+            ref = 0.0
+
+            npt.assert_almost_equal(ref, comp, decimal=14)
+
+
+def test_snippets():
+    # This test function raises an error if there is a considerable loss of precision
+    # that violates the symmetry property of a distance measure.
+    m = 10
+    k = 3
+    s = 3
+    seed = 332
+    np.random.seed(seed)
+    T = np.random.uniform(-1000.0, 1000.0, [64])
+
+    isconstant_custom_func = functools.partial(
+        naive.isconstant_func_stddev_threshold, quantile_threshold=0.05
+    )
+    (
+        ref_snippets,
+        ref_indices,
+        ref_profiles,
+        ref_fractions,
+        ref_areas,
+        ref_regimes,
+    ) = naive.mpdist_snippets(
+        T, m, k, s=s, mpdist_T_subseq_isconstant=isconstant_custom_func
+    )
+    (
+        cmp_snippets,
+        cmp_indices,
+        cmp_profiles,
+        cmp_fractions,
+        cmp_areas,
+        cmp_regimes,
+    ) = stumpy.snippets(T, m, k, s=s, mpdist_T_subseq_isconstant=isconstant_custom_func)
+
+    npt.assert_almost_equal(
+        ref_snippets, cmp_snippets, decimal=config.STUMPY_TEST_PRECISION
+    )
+    npt.assert_almost_equal(
+        ref_indices, cmp_indices, decimal=config.STUMPY_TEST_PRECISION
+    )
+    npt.assert_almost_equal(
+        ref_profiles, cmp_profiles, decimal=config.STUMPY_TEST_PRECISION
+    )
+    npt.assert_almost_equal(
+        ref_fractions, cmp_fractions, decimal=config.STUMPY_TEST_PRECISION
+    )
+    npt.assert_almost_equal(ref_areas, cmp_areas, decimal=config.STUMPY_TEST_PRECISION)
+    npt.assert_almost_equal(ref_regimes, cmp_regimes)
+
+
+@pytest.mark.filterwarnings("ignore", category=NumbaPerformanceWarning)
+@patch("stumpy.config.STUMPY_THREADS_PER_BLOCK", TEST_THREADS_PER_BLOCK)
+def test_distance_symmetry_property_in_gpu():
+    if not cuda.is_available():  # pragma: no cover
+        pytest.skip("Skipping Tests No GPUs Available")
+
+    # This test function raises an error if the distance between a subsequence
+    # and another one does not satisfy the symmetry property.
+    seed = 332
+    np.random.seed(seed)
+    T = np.random.uniform(-1000.0, 1000.0, [64])
+    m = 3
+
+    i, j = 2, 10
+    # M_T, Σ_T = core.compute_mean_std(T, m)
+    # Σ_T[i] is `650.912209452633`
+    # Σ_T[j] is `722.0717285148525`
+
+    # This test raises an error if arithmetic operation in ...
+    # ... `gpu_stump._compute_and_update_PI_kernel` does not
+    # generates the same result if values of variable for mean and std
+    # are swapped.
+
+    T_A = T[i : i + m]
+    T_B = T[j : j + m]
+
+    mp_AB = stumpy.gpu_stump(T_A, m, T_B)
+    mp_BA = stumpy.gpu_stump(T_B, m, T_A)
+
+    d_ij = mp_AB[0, 0]
+    d_ji = mp_BA[0, 0]
+
+    comp = d_ij - d_ji
+    ref = 0.0
+
+    npt.assert_almost_equal(comp, ref, decimal=15)