Refactor/sdp structure cleanup (#1126)

* added new module sdp

* Move sliding-dot-product from core to new module sdp

* move convolve_sdp and add test

* refactor naive function

Author: NimaSarajpoor

stumpy/core.py

@@ -12,10 +12,9 @@
 from numba import cuda, njit, prange
 from scipy import linalg
 from scipy.ndimage import maximum_filter1d, minimum_filter1d
-from scipy.signal import convolve
 from scipy.spatial.distance import cdist
 
-from . import config
+from . import config, sdp
 
 try:
     from numba.cuda.cudadrv.driver import _raise_driver_not_found
@@ -649,36 +648,9 @@ def check_window_size(m, max_size=None, n=None):
             warnings.warn(msg)
 
 
-@njit(fastmath=config.STUMPY_FASTMATH_TRUE)
-def _sliding_dot_product(Q, T):
-    """
-    A Numba JIT-compiled implementation of the sliding window dot product.
-
-    Parameters
-    ----------
-    Q : numpy.ndarray
-        Query array or subsequence
-
-    T : numpy.ndarray
-        Time series or sequence
-
-    Returns
-    -------
-    out : numpy.ndarray
-        Sliding dot product between `Q` and `T`.
-    """
-    m = Q.shape[0]
-    l = T.shape[0] - m + 1
-    out = np.empty(l)
-    for i in range(l):
-        out[i] = np.dot(Q, T[i : i + m])
-
-    return out
-
-
 def sliding_dot_product(Q, T):
     """
-    Use FFT convolution to calculate the sliding window dot product.
+    Calculate the sliding window dot product.
 
     Parameters
     ----------
@@ -692,27 +664,8 @@ def sliding_dot_product(Q, T):
     -------
     output : numpy.ndarray
         Sliding dot product between `Q` and `T`.
-
-    Notes
-    -----
-    Calculate the sliding dot product
-
-    `DOI: 10.1109/ICDM.2016.0179 \
-    <https://www.cs.ucr.edu/~eamonn/PID4481997_extend_Matrix%20Profile_I.pdf>`__
-
-    See Table I, Figure 4
-
-    Following the inverse FFT, Fig. 4 states that only cells [m-1:n]
-    contain valid dot products
-
-    Padding is done automatically in fftconvolve step
     """
-    n = T.shape[0]
-    m = Q.shape[0]
-    Qr = np.flipud(Q)  # Reverse/flip Q
-    QT = convolve(Qr, T)
-
-    return QT.real[m - 1 : n]
+    return sdp._sliding_dot_product(Q, T)
 
 
 @njit(
@@ -1327,7 +1280,7 @@ def _p_norm_distance_profile(Q, T, p=2.0):
             T_squared[i] = (
                 T_squared[i - 1] - T[i - 1] * T[i - 1] + T[i + m - 1] * T[i + m - 1]
             )
-        QT = _sliding_dot_product(Q, T)
+        QT = sdp._njit_sliding_dot_product(Q, T)
         for i in range(l):
             p_norm_profile[i] = Q_squared + T_squared[i] - 2.0 * QT[i]
     else:
@@ -1900,7 +1853,7 @@ def _mass_distance_matrix(
         if np.any(~np.isfinite(Q[i : i + m])):  # pragma: no cover
             distance_matrix[i, :] = np.inf
         else:
-            QT = _sliding_dot_product(Q[i : i + m], T)
+            QT = sdp._njit_sliding_dot_product(Q[i : i + m], T)
             distance_matrix[i, :] = _mass(
                 Q[i : i + m],
                 T,

stumpy/scrump.py

@@ -6,7 +6,7 @@
 import numpy as np
 from numba import njit, prange
 
-from . import config, core
+from . import config, core, sdp
 from .scraamp import prescraamp, scraamp
 from .stump import _stump
 
@@ -235,7 +235,7 @@ def _compute_PI(
     QT = np.empty(w, dtype=np.float64)
     for i in indices[start:stop]:
         Q = T_A[i : i + m]
-        QT[:] = core._sliding_dot_product(Q, T_B)
+        QT[:] = sdp._njit_sliding_dot_product(Q, T_B)
         squared_distance_profile[:] = core._calculate_squared_distance_profile(
             m,
             QT,

stumpy/sdp.py

@@ -0,0 +1,77 @@
+import numpy as np
+from numba import njit
+from scipy.signal import convolve
+
+from . import config
+
+
+@njit(fastmath=config.STUMPY_FASTMATH_TRUE)
+def _njit_sliding_dot_product(Q, T):
+    """
+    A Numba JIT-compiled implementation of the sliding window dot product.
+
+    Parameters
+    ----------
+    Q : numpy.ndarray
+        Query array or subsequence
+
+    T : numpy.ndarray
+        Time series or sequence
+
+    Returns
+    -------
+    out : numpy.ndarray
+        Sliding dot product between `Q` and `T`.
+    """
+    m = Q.shape[0]
+    l = T.shape[0] - m + 1
+    out = np.empty(l)
+    for i in range(l):
+        out[i] = np.dot(Q, T[i : i + m])
+
+    return out
+
+
+def _convolve_sliding_dot_product(Q, T):
+    """
+    Use (direct or FFT) convolution to calculate the sliding window dot product.
+
+    Parameters
+    ----------
+    Q : numpy.ndarray
+        Query array or subsequence
+
+    T : numpy.ndarray
+        Time series or sequence
+
+    Returns
+    -------
+    output : numpy.ndarray
+        Sliding dot product between `Q` and `T`.
+    """
+    n = T.shape[0]
+    m = Q.shape[0]
+    Qr = np.flipud(Q)  # Reverse/flip Q
+    QT = convolve(Qr, T)
+
+    return QT.real[m - 1 : n]
+
+
+def _sliding_dot_product(Q, T):
+    """
+    Compute the sliding dot product between `Q` and `T`
+
+    Parameters
+    ----------
+    Q : numpy.ndarray
+        Query array or subsequence
+
+    T : numpy.ndarray
+        Time series or sequence
+
+    Returns
+    -------
+    out : numpy.ndarray
+        Sliding dot product between `Q` and `T`.
+    """
+    return _convolve_sliding_dot_product(Q, T)

tests/naive.py

@@ -2379,3 +2379,11 @@ def mpdist_custom_func(P_ABBA, m, percentage, n_A, n_B):
         MPdist = P_ABBA[k]
 
     return MPdist
+
+
+def rolling_window_dot_product(Q, T):
+    window = len(Q)
+    result = np.zeros(len(T) - window + 1)
+    for i in range(len(result)):
+        result[i] = np.dot(T[i : i + window], Q)
+    return result

tests/test_core.py

@@ -35,14 +35,6 @@ def _gpu_searchsorted_kernel(a, v, bfs, nlevel, is_left, idx):
 TEST_THREADS_PER_BLOCK = 10
 
 
-def naive_rolling_window_dot_product(Q, T):
-    window = len(Q)
-    result = np.zeros(len(T) - window + 1)
-    for i in range(len(result)):
-        result[i] = np.dot(T[i : i + window], Q)
-    return result
-
-
 def naive_compute_mean_std_multidimensional(T, m):
     n = T.shape[1]
     nrows, ncols = T.shape
@@ -208,16 +200,9 @@ def test_check_window_size_excl_zone():
         core.check_window_size(m, max_size=len(T), n=len(T))
 
 
-@pytest.mark.parametrize("Q, T", test_data)
-def test_njit_sliding_dot_product(Q, T):
-    ref_mp = naive_rolling_window_dot_product(Q, T)
-    comp_mp = core._sliding_dot_product(Q, T)
-    npt.assert_almost_equal(ref_mp, comp_mp)
-
-
 @pytest.mark.parametrize("Q, T", test_data)
 def test_sliding_dot_product(Q, T):
-    ref_mp = naive_rolling_window_dot_product(Q, T)
+    ref_mp = naive.rolling_window_dot_product(Q, T)
     comp_mp = core.sliding_dot_product(Q, T)
     npt.assert_almost_equal(ref_mp, comp_mp)
 

tests/test_precision.py

@@ -8,7 +8,7 @@
 import pytest
 from numba import cuda
 
-from stumpy import cache, config, core, fastmath
+from stumpy import cache, config, core, fastmath, sdp
 
 if cuda.is_available():
     from stumpy.gpu_stump import gpu_stump
@@ -90,7 +90,7 @@ def test_calculate_squared_distance():
     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)
+            QT_i = sdp._njit_sliding_dot_product(T[i : i + m], T)
             dist_ij = core._calculate_squared_distance(
                 m,
                 QT_i[j],
@@ -102,7 +102,7 @@ def test_calculate_squared_distance():
                 T_subseq_isconstant[j],
             )
 
-            QT_j = core._sliding_dot_product(T[j : j + m], T)
+            QT_j = sdp._njit_sliding_dot_product(T[j : j + m], T)
             dist_ji = core._calculate_squared_distance(
                 m,
                 QT_j[i],

tests/test_sdp.py

@@ -0,0 +1,36 @@
+import naive
+import numpy as np
+import pytest
+from numpy import testing as npt
+
+from stumpy import sdp
+
+test_data = [
+    (np.array([-1, 1, 2], dtype=np.float64), np.array(range(5), dtype=np.float64)),
+    (
+        np.array([9, 8100, -60], dtype=np.float64),
+        np.array([584, -11, 23, 79, 1001], dtype=np.float64),
+    ),
+    (np.random.uniform(-1000, 1000, [8]), np.random.uniform(-1000, 1000, [64])),
+]
+
+
+@pytest.mark.parametrize("Q, T", test_data)
+def test_njit_sliding_dot_product(Q, T):
+    ref_mp = naive.rolling_window_dot_product(Q, T)
+    comp_mp = sdp._njit_sliding_dot_product(Q, T)
+    npt.assert_almost_equal(ref_mp, comp_mp)
+
+
+@pytest.mark.parametrize("Q, T", test_data)
+def test_convolve_sliding_dot_product(Q, T):
+    ref_mp = naive.rolling_window_dot_product(Q, T)
+    comp_mp = sdp._convolve_sliding_dot_product(Q, T)
+    npt.assert_almost_equal(ref_mp, comp_mp)
+
+
+@pytest.mark.parametrize("Q, T", test_data)
+def test_sliding_dot_product(Q, T):
+    ref_mp = naive.rolling_window_dot_product(Q, T)
+    comp_mp = sdp._sliding_dot_product(Q, T)
+    npt.assert_almost_equal(ref_mp, comp_mp)