Fixed #880 Add param isconstant to gpu_ostinato and gpu_stimp (#883)

* add param isconstant to gpu_ostinato

* add isconstant to gpu_stimp and add test functions

* minor changes

* minor change

* Removed bad test function

Author: NimaSarajpoor

stumpy/gpu_ostinato.py

@@ -8,8 +8,11 @@
 from .ostinato import _get_central_motif, _ostinato
 
 
-@core.non_normalized(gpu_aamp_ostinato)
-def gpu_ostinato(Ts, m, device_id=0, normalize=True, p=2.0):
+@core.non_normalized(
+    gpu_aamp_ostinato,
+    exclude=["normalize", "p", "Ts_subseq_isconstant"],
+)
+def gpu_ostinato(Ts, m, device_id=0, normalize=True, p=2.0, Ts_subseq_isconstant=None):
     """
     Find the z-normalized consensus motif of multiple time series with one or more GPU
     devices
@@ -43,6 +46,9 @@ def gpu_ostinato(Ts, m, device_id=0, normalize=True, p=2.0):
         and the Euclidean distance, respectively. This parameter is ignored when
         `normalize == True`.
 
+    Ts_subseq_isconstant : list, default None
+        A list of rolling window isconstant for each time series in `Ts`.
+
     Returns
     -------
     central_radius : float
@@ -99,11 +105,15 @@ def gpu_ostinato(Ts, m, device_id=0, normalize=True, p=2.0):
     if not isinstance(Ts, list):  # pragma: no cover
         raise ValueError(f"`Ts` is of type `{type(Ts)}` but a `list` is expected")
 
-    Ts_subseq_isconstant = [None] * len(Ts)
+    if Ts_subseq_isconstant is None:
+        Ts_subseq_isconstant = [None] * len(Ts)
+
     M_Ts = [None] * len(Ts)
     Σ_Ts = [None] * len(Ts)
     for i, T in enumerate(Ts):
-        Ts[i], M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(T, m)
+        Ts[i], M_Ts[i], Σ_Ts[i], Ts_subseq_isconstant[i] = core.preprocess(
+            T, m, T_subseq_isconstant=Ts_subseq_isconstant[i]
+        )
 
     bsf_radius, bsf_Ts_idx, bsf_subseq_idx = _ostinato(
         Ts, m, M_Ts, Σ_Ts, Ts_subseq_isconstant, device_id=device_id, mp_func=gpu_stump

stumpy/gpu_stimp.py

@@ -10,7 +10,7 @@
 
 @core.non_normalized(
     gpu_aamp_stimp,
-    exclude=["pre_scrump", "normalize", "p", "pre_scraamp"],
+    exclude=["pre_scrump", "normalize", "p", "pre_scraamp", "T_subseq_isconstant_func"],
     replace={"pre_scrump": "pre_scraamp"},
 )
 class gpu_stimp(_stimp):
@@ -24,16 +24,16 @@ class gpu_stimp(_stimp):
     T : numpy.ndarray
         The time series or sequence for which to compute the pan matrix profile
 
-    m_start : int, default 3
+    min_m : int, default 3
         The starting (or minimum) subsequence window size for which a matrix profile
         may be computed
 
-    m_stop : int, default None
+    max_m : int, default None
         The stopping (or maximum) subsequence window size for which a matrix profile
         may be computed. When `m_stop = Non`, this is set to the maximum allowable
         subsequence window size
 
-    m_step : int, default 1
+    step : int, default 1
         The step between subsequence window sizes
 
     device_id : int or list, default 0
@@ -53,6 +53,14 @@ class gpu_stimp(_stimp):
         and the Euclidean distance, respectively. This parameter is ignored when
         `normalize == True`.
 
+    T_subseq_isconstant_func : function, default None
+        A custom, user-defined function that returns a boolean array that indicates
+        whether a subsequence in `T` is constant (True). The function must only take
+        two arguments, `a`, a 1-D array, and `w`, the window size, while additional
+        arguments may be specified by currying the user-defined function using
+        `functools.partial`. Any subsequence with at least one np.nan/np.inf will
+        automatically have its corresponding value set to False in this boolean array.
+
     Attributes
     ----------
     PAN_ : numpy.ndarray
@@ -106,6 +114,7 @@ def __init__(
         device_id=0,
         normalize=True,
         p=2.0,
+        T_subseq_isconstant_func=None,
     ):
         """
         Initialize the `stimp` object and compute the Pan Matrix Profile
@@ -144,6 +153,15 @@ def __init__(
             is typically used with `p` being 1 or 2, which correspond to the Manhattan
             distance and the Euclidean distance, respectively. This parameter is ignored
             when `normalize == True`.
+
+        T_subseq_isconstant_func : function, default None
+            A custom, user-defined function that returns a boolean array that indicates
+            whether a subsequence in `T` is constant (True). The function must only take
+            two arguments, `a`, a 1-D array, and `w`, the window size, while additional
+            arguments may be specified by currying the user-defined function using
+            `functools.partial`. Any subsequence with at least one np.nan/np.inf will
+            automatically have its corresponding value set to False in this boolean
+            array.
         """
         super().__init__(
             T,
@@ -153,5 +171,6 @@ def __init__(
             percentage=1.0,
             pre_scrump=False,
             device_id=device_id,
+            T_subseq_isconstant_func=T_subseq_isconstant_func,
             mp_func=gpu_stump,
         )

tests/test_gpu_ostinato.py

@@ -1,3 +1,5 @@
+import functools
+
 import numpy as np
 import numpy.testing as npt
 from numba import cuda
@@ -52,3 +54,63 @@ def test_deterministic_gpu_ostinato(seed):
     npt.assert_almost_equal(ref_radius, comp_radius)
     npt.assert_almost_equal(ref_Ts_idx, comp_Ts_idx)
     npt.assert_almost_equal(ref_subseq_idx, comp_subseq_idx)
+
+
+@pytest.mark.filterwarnings("ignore", category=NumbaPerformanceWarning)
+@pytest.mark.parametrize(
+    "seed", np.random.choice(np.arange(10000), size=25, replace=False)
+)
+@patch("stumpy.config.STUMPY_THREADS_PER_BLOCK", TEST_THREADS_PER_BLOCK)
+def test_random_gpu_ostinato_with_isconstant(seed):
+    isconstant_custom_func = functools.partial(
+        naive.isconstant_func_stddev_threshold, quantile_threshold=0.05
+    )
+
+    m = 50
+    np.random.seed(seed)
+    Ts = [np.random.rand(n) for n in [64, 128, 256]]
+    Ts_subseq_isconstant = [isconstant_custom_func for _ in range(len(Ts))]
+
+    ref_radius, ref_Ts_idx, ref_subseq_idx = naive.ostinato(
+        Ts, m, Ts_subseq_isconstant=Ts_subseq_isconstant
+    )
+    comp_radius, comp_Ts_idx, comp_subseq_idx = gpu_ostinato(
+        Ts, m, Ts_subseq_isconstant=Ts_subseq_isconstant
+    )
+
+    npt.assert_almost_equal(ref_radius, comp_radius)
+    npt.assert_almost_equal(ref_Ts_idx, comp_Ts_idx)
+    npt.assert_almost_equal(ref_subseq_idx, comp_subseq_idx)
+
+
+@pytest.mark.filterwarnings("ignore", category=NumbaPerformanceWarning)
+@pytest.mark.parametrize("seed", [79, 109, 112, 133, 151, 161, 251, 275, 309, 355])
+@patch("stumpy.config.STUMPY_THREADS_PER_BLOCK", TEST_THREADS_PER_BLOCK)
+def test_deterministic_gpu_ostinato_with_isconstant(seed):
+    isconstant_custom_func = functools.partial(
+        naive.isconstant_func_stddev_threshold, quantile_threshold=0.05
+    )
+
+    m = 50
+    np.random.seed(seed)
+    Ts = [np.random.rand(n) for n in [64, 128, 256]]
+
+    l = 64 - m + 1
+    subseq_isconsant = np.full(l, 0, dtype=bool)
+    subseq_isconsant[np.random.randint(0, l)] = True
+    Ts_subseq_isconstant = [
+        subseq_isconsant,
+        None,
+        isconstant_custom_func,
+    ]
+
+    ref_radius, ref_Ts_idx, ref_subseq_idx = naive.ostinato(
+        Ts, m, Ts_subseq_isconstant=Ts_subseq_isconstant
+    )
+    comp_radius, comp_Ts_idx, comp_subseq_idx = gpu_ostinato(
+        Ts, m, Ts_subseq_isconstant=Ts_subseq_isconstant
+    )
+
+    npt.assert_almost_equal(ref_radius, comp_radius)
+    npt.assert_almost_equal(ref_Ts_idx, comp_Ts_idx)
+    npt.assert_almost_equal(ref_subseq_idx, comp_subseq_idx)

tests/test_gpu_stimp.py

@@ -1,3 +1,4 @@
+import functools
 from unittest.mock import patch
 
 import numpy as np
@@ -71,3 +72,60 @@ def test_gpu_stimp(T):
     naive.replace_inf(cmp_pan)
 
     npt.assert_almost_equal(ref_pan, cmp_pan)
+
+
+@pytest.mark.filterwarnings("ignore", category=NumbaPerformanceWarning)
+@patch("stumpy.config.STUMPY_THREADS_PER_BLOCK", TEST_THREADS_PER_BLOCK)
+def test_gpu_stimp_with_isconstant():
+    T = np.random.uniform(-1, 1, [64])
+    isconstant_func = functools.partial(
+        naive.isconstant_func_stddev_threshold, stddev_threshold=0.5
+    )
+
+    threshold = 0.2
+    min_m = 3
+    n = T.shape[0] - min_m + 1
+
+    pan = gpu_stimp(
+        T,
+        min_m=min_m,
+        max_m=None,
+        step=1,
+        # normalize=True,
+        T_subseq_isconstant_func=isconstant_func,
+    )
+
+    for i in range(n):
+        pan.update()
+
+    ref_PAN = np.full((pan.M_.shape[0], T.shape[0]), fill_value=np.inf)
+
+    for idx, m in enumerate(pan.M_[:n]):
+        zone = int(np.ceil(m / 4))
+        ref_mp = naive.stump(
+            T,
+            m,
+            T_B=None,
+            exclusion_zone=zone,
+            T_A_subseq_isconstant=isconstant_func,
+        )
+        ref_PAN[pan._bfs_indices[idx], : ref_mp.shape[0]] = ref_mp[:, 0]
+
+    # Compare raw pan
+    cmp_PAN = pan._PAN
+
+    naive.replace_inf(ref_PAN)
+    naive.replace_inf(cmp_PAN)
+
+    npt.assert_almost_equal(ref_PAN, cmp_PAN)
+
+    # Compare transformed pan
+    cmp_pan = pan.PAN_
+    ref_pan = naive.transform_pan(
+        pan._PAN, pan._M, threshold, pan._bfs_indices, pan._n_processed
+    )
+
+    naive.replace_inf(ref_pan)
+    naive.replace_inf(cmp_pan)
+
+    npt.assert_almost_equal(ref_pan, cmp_pan)