diff --git a/stumpy/core.py b/stumpy/core.py index 57d05f989..320885999 100644 --- a/stumpy/core.py +++ b/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 - ρ)) diff --git a/stumpy/gpu_stump.py b/stumpy/gpu_stump.py index d66283d7f..7f103b7ee 100644 --- a/stumpy/gpu_stump.py +++ b/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 - ρ) diff --git a/tests/test_precision.py b/tests/test_precision.py index c819c7e22..c87985092 100644 --- a/tests/test_precision.py +++ b/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)