add param isconstant to ostinato (#864)

* add param isconstant to naive function ostinato

* add test func, expecting error

* fix performant func and update test func

* minor changes

* add param isconstant to ostinatoed, also fix converage

* fix decorator

* fix coverage

Author: NimaSarajpoor

stumpy/ostinato.py

@@ -236,7 +236,14 @@ def _ostinato(
         else:
             h = 0
 
-        mp = partial_mp_func(Ts[j], m, Ts[h], ignore_trivial=False)
+        mp = partial_mp_func(
+            Ts[j],
+            m,
+            Ts[h],
+            ignore_trivial=False,
+            T_A_subseq_isconstant=Ts_subseq_isconstant[j],
+            T_B_subseq_isconstant=Ts_subseq_isconstant[h],
+        )
         si = np.argsort(mp[:, 0])
         for q in si:
             radius = mp[q, 0]
@@ -271,8 +278,11 @@ def _ostinato(
     return bsf_radius, bsf_Ts_idx, bsf_subseq_idx
 
 
-@core.non_normalized(aamp_ostinato)
-def ostinato(Ts, m, normalize=True, p=2.0):
+@core.non_normalized(
+    aamp_ostinato,
+    exclude=["normalize", "p", "Ts_subseq_isconstant"],
+)
+def ostinato(Ts, m, normalize=True, p=2.0, Ts_subseq_isconstant=None):
     """
     Find the z-normalized consensus motif of multiple time series
 
@@ -297,6 +307,9 @@ def ostinato(Ts, m, normalize=True, p=2.0):
         The p-norm to apply for computing the Minkowski distance. 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
@@ -352,9 +365,12 @@ def ostinato(Ts, m, normalize=True, p=2.0):
 
     M_Ts = [None] * len(Ts)
     Σ_Ts = [None] * len(Ts)
-    Ts_subseq_isconstant = [None] * len(Ts)
+    if Ts_subseq_isconstant is None:
+        Ts_subseq_isconstant = [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
@@ -371,8 +387,11 @@ def ostinato(Ts, m, normalize=True, p=2.0):
     return central_radius, central_Ts_idx, central_subseq_idx
 
 
-@core.non_normalized(aamp_ostinatoed)
-def ostinatoed(client, Ts, m, normalize=True, p=2.0):
+@core.non_normalized(
+    aamp_ostinatoed,
+    exclude=["normalize", "p", "Ts_subseq_isconstant"],
+)
+def ostinatoed(client, Ts, m, normalize=True, p=2.0, Ts_subseq_isconstant=None):
     """
     Find the z-normalized consensus motif of multiple time series with a distributed
     cluster
@@ -403,6 +422,9 @@ def ostinatoed(client, Ts, m, normalize=True, p=2.0):
         The p-norm to apply for computing the Minkowski distance. 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
@@ -461,9 +483,13 @@ def ostinatoed(client, Ts, m, normalize=True, p=2.0):
 
     M_Ts = [None] * len(Ts)
     Σ_Ts = [None] * len(Ts)
-    Ts_subseq_isconstant = [None] * len(Ts)
+
+    if Ts_subseq_isconstant is None:
+        Ts_subseq_isconstant = [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,

tests/naive.py

@@ -1063,7 +1063,7 @@ def left_I_(self):
         return self._left_I.astype(np.int64)
 
 
-def across_series_nearest_neighbors(Ts, Ts_idx, subseq_idx, m):
+def across_series_nearest_neighbors(Ts, Ts_idx, subseq_idx, m, Ts_subseq_isconstant):
     """
     For multiple time series find, per individual time series, the subsequences closest
     to a query.
@@ -1085,6 +1085,9 @@ def across_series_nearest_neighbors(Ts, Ts_idx, subseq_idx, m):
     m : int
         Subsequence window size
 
+    Ts_subseq_isconstant : list
+        A list of `T_subseq_isconstant`, where the i-th item corresponds to `Ts[i]`
+
     Returns
     -------
     nns_radii : ndarray
@@ -1096,19 +1099,33 @@ def across_series_nearest_neighbors(Ts, Ts_idx, subseq_idx, m):
         `Ts[Ts_idx][subseq_idx : subseq_idx + m]`
     """
     k = len(Ts)
+
     Q = Ts[Ts_idx][subseq_idx : subseq_idx + m]
+    Q_subseq_isconstant = Ts_subseq_isconstant[Ts_idx][subseq_idx]
+
     nns_radii = np.zeros(k, dtype=np.float64)
     nns_subseq_idx = np.zeros(k, dtype=np.int64)
 
     for i in range(k):
         dist_profile = distance_profile(Q, Ts[i], len(Q))
+        for j in range(len(dist_profile)):
+            if np.isfinite(dist_profile[j]):
+                if Q_subseq_isconstant and Ts_subseq_isconstant[i][j]:
+                    dist_profile[j] = 0
+                elif Q_subseq_isconstant or Ts_subseq_isconstant[i][j]:
+                    dist_profile[j] = np.sqrt(m)
+                else:  # pragma: no cover
+                    pass
+
         nns_subseq_idx[i] = np.argmin(dist_profile)
         nns_radii[i] = dist_profile[nns_subseq_idx[i]]
 
     return nns_radii, nns_subseq_idx
 
 
-def get_central_motif(Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m):
+def get_central_motif(
+    Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m, Ts_subseq_isconstant
+):
     """
     Compare subsequences with the same radius and return the most central motif
 
@@ -1129,6 +1146,9 @@ def get_central_motif(Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m):
     m : int
         Window size
 
+    Ts_subseq_isconstant : list
+        A list of boolean arrays, each corresponds to a time series in `Ts`
+
     Returns
     -------
     bsf_radius : float
@@ -1142,7 +1162,7 @@ def get_central_motif(Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m):
         series `bsf_Ts_idx` that contains it
     """
     bsf_nns_radii, bsf_nns_subseq_idx = across_series_nearest_neighbors(
-        Ts, bsf_Ts_idx, bsf_subseq_idx, m
+        Ts, bsf_Ts_idx, bsf_subseq_idx, m, Ts_subseq_isconstant
     )
     bsf_nns_mean_radii = bsf_nns_radii.mean()
 
@@ -1151,7 +1171,7 @@ def get_central_motif(Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m):
 
     for Ts_idx, subseq_idx in zip(candidate_nns_Ts_idx, candidate_nns_subseq_idx):
         candidate_nns_radii, _ = across_series_nearest_neighbors(
-            Ts, Ts_idx, subseq_idx, m
+            Ts, Ts_idx, subseq_idx, m, Ts_subseq_isconstant
         )
         if (
             np.isclose(candidate_nns_radii.max(), bsf_radius)
@@ -1164,7 +1184,7 @@ def get_central_motif(Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m):
     return bsf_radius, bsf_Ts_idx, bsf_subseq_idx
 
 
-def consensus_search(Ts, m):
+def consensus_search(Ts, m, Ts_subseq_isconstant):
     """
     Brute force consensus motif from
     <https://www.cs.ucr.edu/~eamonn/consensus_Motif_ICDM_Long_version.pdf>
@@ -1184,7 +1204,13 @@ def consensus_search(Ts, m):
         radii = np.zeros(len(Ts[j]) - m + 1)
         for i in range(k):
             if i != j:
-                mp = stump(Ts[j], m, Ts[i])
+                mp = stump(
+                    Ts[j],
+                    m,
+                    Ts[i],
+                    T_A_subseq_isconstant=Ts_subseq_isconstant[j],
+                    T_B_subseq_isconstant=Ts_subseq_isconstant[i],
+                )
                 radii = np.maximum(radii, mp[:, 0])
         min_radius_idx = np.argmin(radii)
         min_radius = radii[min_radius_idx]
@@ -1196,10 +1222,24 @@ def consensus_search(Ts, m):
     return bsf_radius, bsf_Ts_idx, bsf_subseq_idx
 
 
-def ostinato(Ts, m):
-    bsf_radius, bsf_Ts_idx, bsf_subseq_idx = consensus_search(Ts, m)
+def ostinato(Ts, m, Ts_subseq_isconstant=None):
+    if Ts_subseq_isconstant is None:
+        Ts_subseq_isconstant = [None] * len(Ts)
+
+    Ts_subseq_isconstant = [
+        rolling_isconstant(Ts[i], m, Ts_subseq_isconstant[i]) for i in range(len(Ts))
+    ]
+
+    bsf_radius, bsf_Ts_idx, bsf_subseq_idx = consensus_search(
+        Ts, m, Ts_subseq_isconstant
+    )
     radius, Ts_idx, subseq_idx = get_central_motif(
-        Ts, bsf_radius, bsf_Ts_idx, bsf_subseq_idx, m
+        Ts,
+        bsf_radius,
+        bsf_Ts_idx,
+        bsf_subseq_idx,
+        m,
+        Ts_subseq_isconstant,
     )
     return radius, Ts_idx, subseq_idx
 

tests/test_ostinato.py

@@ -1,3 +1,5 @@
+import functools
+
 import naive
 import numpy as np
 import numpy.testing as npt
@@ -78,3 +80,60 @@ def test_deterministic_ostinatoed(seed, dask_cluster):
         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.parametrize(
+    "seed", np.random.choice(np.arange(10000), size=25, replace=False)
+)
+def test_random_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 = stumpy.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.parametrize("seed", [79, 109, 112, 133, 151, 161, 251, 275, 309, 355])
+def test_deterministic_ostinatoed_with_isconstant(seed, dask_cluster):
+    isconstant_custom_func = functools.partial(
+        naive.isconstant_func_stddev_threshold, quantile_threshold=0.05
+    )
+
+    with Client(dask_cluster) as dask_client:
+        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 = stumpy.ostinatoed(
+            dask_client, 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)