Migrate to modern NumPy random API and add STUMPY_SEED test reproducibility (#707, #1112)

Author: Vansh-Sharmaa

README.rst

@@ -82,7 +82,7 @@ Typical usage (1-dimensional time series data) with `STUMP <https://stumpy.readt
     import numpy as np
     
     if __name__ == "__main__":
-        your_time_series = np.random.rand(10000)
+        your_time_series = np.random.default_rng().random(10000)
         window_size = 50  # Approximately, how many data points might be found in a pattern 
     
         matrix_profile = stumpy.stump(your_time_series, m=window_size)
@@ -97,7 +97,7 @@ Distributed usage for 1-dimensional time series data with Dask Distributed via `
 
     if __name__ == "__main__":
         with Client() as dask_client:
-            your_time_series = np.random.rand(10000)
+            your_time_series = np.random.default_rng().random(10000)
             window_size = 50  # Approximately, how many data points might be found in a pattern 
     
             matrix_profile = stumpy.stumped(dask_client, your_time_series, m=window_size)
@@ -111,7 +111,7 @@ GPU usage for 1-dimensional time series data with `GPU-STUMP <https://stumpy.rea
     from numba import cuda
 
     if __name__ == "__main__":
-        your_time_series = np.random.rand(10000)
+        your_time_series = np.random.default_rng().random(10000)
         window_size = 50  # Approximately, how many data points might be found in a pattern
         all_gpu_devices = [device.id for device in cuda.list_devices()]  # Get a list of all available GPU devices
 
@@ -125,7 +125,7 @@ Multi-dimensional time series data with `MSTUMP <https://stumpy.readthedocs.io/e
     import numpy as np
 
     if __name__ == "__main__":
-        your_time_series = np.random.rand(3, 1000)  # Each row represents data from a different dimension while each column represents data from the same dimension
+        your_time_series = np.random.default_rng().random((3, 1000))  # Each row represents data from a different dimension while each column represents data from the same dimension
         window_size = 50  # Approximately, how many data points might be found in a pattern
 
         matrix_profile, matrix_profile_indices = stumpy.mstump(your_time_series, m=window_size)
@@ -140,7 +140,7 @@ Distributed multi-dimensional time series data analysis with Dask Distributed `M
 
     if __name__ == "__main__":
         with Client() as dask_client:
-            your_time_series = np.random.rand(3, 1000)   # Each row represents data from a different dimension while each column represents data from the same dimension
+            your_time_series = np.random.default_rng().random((3, 1000))   # Each row represents data from a different dimension while each column represents data from the same dimension
             window_size = 50  # Approximately, how many data points might be found in a pattern
 
             matrix_profile, matrix_profile_indices = stumpy.mstumped(dask_client, your_time_series, m=window_size)
@@ -153,7 +153,7 @@ Time Series Chains with `Anchored Time Series Chains (ATSC) <https://stumpy.read
     import numpy as np
     
     if __name__ == "__main__":
-        your_time_series = np.random.rand(10000)
+        your_time_series = np.random.default_rng().random(10000)
         window_size = 50  # Approximately, how many data points might be found in a pattern 
         
         matrix_profile = stumpy.stump(your_time_series, m=window_size)
@@ -174,7 +174,7 @@ Semantic Segmentation with `Fast Low-cost Unipotent Semantic Segmentation (FLUSS
     import numpy as np
 
     if __name__ == "__main__":
-        your_time_series = np.random.rand(10000)
+        your_time_series = np.random.default_rng().random(10000)
         window_size = 50  # Approximately, how many data points might be found in a pattern
 
         matrix_profile = stumpy.stump(your_time_series, m=window_size)
@@ -236,7 +236,7 @@ In order to fully understand and appreciate the underlying algorithms and applic
 Performance
 -----------
 
-We tested the performance of computing the exact matrix profile using the Numba JIT compiled version of the code on randomly generated time series data with various lengths (i.e., ``np.random.rand(n)``) along with different `CPU and GPU hardware resources <hardware_>`_. 
+We tested the performance of computing the exact matrix profile using the Numba JIT compiled version of the code on randomly generated time series data with various lengths (i.e., ``np.random.default_rng().random(n)``) along with different `CPU and GPU hardware resources <hardware_>`_. 
 
 .. image:: https://raw.githubusercontent.com/stumpy-dev/stumpy/main/docs/images/performance.png
     :alt: STUMPY Performance Plot

conftest.py

@@ -4,3 +4,18 @@
 # to fix eventual module import errors that can arise, for example when
 # running tests from inside VS code.
 # See https://stackoverflow.com/a/34520971
+
+import os
+import numpy as np
+
+# Set STUMPY_SEED for reproducible test failures.
+# If not already set (e.g., by test.sh), generate a random seed.
+# The seed is printed so it can be noted and reused to reproduce failures:
+#   STUMPY_SEED=12345 pytest tests/test_stump.py
+if "STUMPY_SEED" not in os.environ:
+    os.environ["STUMPY_SEED"] = str(np.random.default_rng().integers(2**31))
+
+
+def pytest_configure(config):
+    """Print STUMPY_SEED so failed tests can be reproduced."""
+    print(f"\nSTUMPY_SEED={os.environ['STUMPY_SEED']}")

stumpy/config.py

@@ -40,6 +40,8 @@
 STUMPY_EXCL_ZONE_DENOM = _STUMPY_DEFAULTS["STUMPY_EXCL_ZONE_DENOM"]
 STUMPY_FASTMATH_TRUE = _STUMPY_DEFAULTS["STUMPY_FASTMATH_TRUE"]
 STUMPY_FASTMATH_FLAGS = _STUMPY_DEFAULTS["STUMPY_FASTMATH_FLAGS"]
+RNG = np.random.default_rng()
+
 
 
 def _reset(var=None):

stumpy/core.py

@@ -3717,7 +3717,8 @@ def check_ignore_trivial(T_A, T_B, ignore_trivial):
     import numpy as np
     import warnings
 
-    T = np.random.rand(10_000)
+    rng = np.random.default_rng()
+    T = rng.random(10_000)
     m = 50
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="Arrays T_A, T_B are equal")
@@ -4493,7 +4494,8 @@ def check_self_join(ignore_trivial):
     import numpy as np
     import warnings
 
-    T = np.random.rand(10_000)
+    rng = np.random.default_rng()
+    T = rng.random(10_000)
     m = 50
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="`ignore_trivial` cannot be `False`")

stumpy/floss.py

@@ -84,13 +84,13 @@ def _iac(
     IAC : numpy.ndarray
         Idealized arc curve (IAC)
     """
-    np.random.seed(seed)
+    rng = np.random.default_rng(seed)
 
-    I = np.random.randint(0, width, size=width, dtype=np.int64)
+    I = rng.integers(0, width, size=width, dtype=np.int64)
     if bidirectional is False:  # Idealized 1-dimensional matrix profile index
         I[:-1] = width
         for i in range(width - 1):
-            I[i] = np.random.randint(i + 1, width, dtype=np.int64)
+            I[i] = rng.integers(i + 1, width, dtype=np.int64)
 
     target_AC = _nnmark(I)
 
@@ -99,7 +99,7 @@ def _iac(
         hist_dist = scipy.stats.rv_histogram(
             (target_AC, np.append(np.arange(width), width))
         )
-        data = hist_dist.rvs(size=n_samples)
+        data = hist_dist.rvs(size=n_samples, random_state=rng)
         a, b, c, d = scipy.stats.beta.fit(data, floc=0, fscale=width)
 
         params[i, 0] = a

stumpy/scraamp.py

@@ -78,7 +78,7 @@ def _preprocess_prescraamp(T_A, m, T_B=None, s=None):
         else:  # AB-join
             s = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
 
-    indices = np.random.permutation(range(0, l, s)).astype(np.int64)
+    indices = config.RNG.permutation(np.arange(0, l, s)).astype(np.int64)
 
     return (T_A, T_B, T_A_subseq_isfinite, T_B_subseq_isfinite, indices, s, excl_zone)
 
@@ -718,8 +718,8 @@ def __init__(
             core._merge_topk_PI(self._P, P, self._I, I)
 
         if self._ignore_trivial:
-            self._diags = np.random.permutation(
-                range(self._excl_zone + 1, self._n_A - self._m + 1)
+            self._diags = config.RNG.permutation(
+                np.arange(self._excl_zone + 1, self._n_A - self._m + 1)
             ).astype(np.int64)
             if self._diags.shape[0] == 0:  # pragma: no cover
                 max_m = core.get_max_window_size(self._T_A.shape[0])
@@ -728,8 +728,8 @@ def __init__(
                     f"Please try a value of `m <= {max_m}`"
                 )
         else:
-            self._diags = np.random.permutation(
-                range(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
+            self._diags = config.RNG.permutation(
+                np.arange(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
             ).astype(np.int64)
 
         self._n_threads = numba.config.NUMBA_NUM_THREADS

stumpy/scrump.py

@@ -116,7 +116,7 @@ def _preprocess_prescrump(
         else:  # AB-join
             s = int(np.ceil(m / config.STUMPY_EXCL_ZONE_DENOM))
 
-    indices = np.random.permutation(range(0, l, s)).astype(np.int64)
+    indices = config.RNG.permutation(np.arange(0, l, s)).astype(np.int64)
 
     return (
         T_A,
@@ -997,8 +997,8 @@ def __init__(
             core._merge_topk_PI(self._P, P, self._I, I)
 
         if self._ignore_trivial:
-            self._diags = np.random.permutation(
-                range(self._excl_zone + 1, self._n_A - self._m + 1)
+            self._diags = config.RNG.permutation(
+                np.arange(self._excl_zone + 1, self._n_A - self._m + 1)
             ).astype(np.int64)
             if self._diags.shape[0] == 0:  # pragma: no cover
                 max_m = core.get_max_window_size(self._T_A.shape[0])
@@ -1007,8 +1007,8 @@ def __init__(
                     f"Please try a value of `m <= {max_m}`"
                 )
         else:
-            self._diags = np.random.permutation(
-                range(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
+            self._diags = config.RNG.permutation(
+                np.arange(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
             ).astype(np.int64)
 
         self._n_threads = numba.config.NUMBA_NUM_THREADS

test.sh

@@ -347,6 +347,10 @@ convert_notebooks()
 #   Main  #
 ###########
 
+# Set STUMPY_SEED for reproducible test failures
+export STUMPY_SEED=${STUMPY_SEED:-$RANDOM}
+echo "STUMPY_SEED=$STUMPY_SEED"
+
 if [[ $test_mode == "show" ]]; then
     echo "Show development/test environment"
     show

tests/conftest.py

@@ -0,0 +1,13 @@
+import os
+import pytest
+import numpy as np
+from stumpy import config
+
+seed = int(os.environ.get("STUMPY_SEED", np.random.default_rng().integers(2**31)))
+RNG = np.random.default_rng(seed)
+config.RNG = RNG
+
+
+@pytest.fixture(scope="module")
+def rng_state():
+    yield RNG.bit_generator.state

tests/naive.py

@@ -1818,7 +1818,7 @@ def prescrump(
     P = np.full((l, k), np.inf, dtype=np.float64)
     I = np.full((l, k), -1, dtype=np.int64)
 
-    for i in np.random.permutation(range(0, l, s)):
+    for i in config.RNG.permutation(np.arange(0, l, s)):
         distance_profile = dist_matrix[i]
         if exclusion_zone is not None:
             apply_exclusion_zone(distance_profile, i, exclusion_zone, np.inf)
@@ -1914,13 +1914,13 @@ def scrump(
                     pass
 
     if exclusion_zone is not None:
-        diags = np.random.permutation(range(exclusion_zone + 1, n_A - m + 1)).astype(
-            np.int64
-        )
+        diags = config.RNG.permutation(
+            np.arange(exclusion_zone + 1, n_A - m + 1)
+        ).astype(np.int64)
     else:
-        diags = np.random.permutation(range(-(n_A - m + 1) + 1, n_B - m + 1)).astype(
-            np.int64
-        )
+        diags = config.RNG.permutation(
+            np.arange(-(n_A - m + 1) + 1, n_B - m + 1)
+        ).astype(np.int64)
 
     n_chunks = int(np.ceil(1.0 / percentage))
     ndist_counts = core._count_diagonal_ndist(diags, m, n_A, n_B)
@@ -1981,7 +1981,7 @@ def prescraamp(T_A, m, T_B, s, exclusion_zone=None, p=2.0, k=1):
     P = np.full((l, k), np.inf, dtype=np.float64)
     I = np.full((l, k), -1, dtype=np.int64)
 
-    for i in np.random.permutation(range(0, l, s)):
+    for i in config.RNG.permutation(np.arange(0, l, s)):
         distance_profile = distance_matrix[i]
         if exclusion_zone is not None:
             apply_exclusion_zone(distance_profile, i, exclusion_zone, np.inf)
@@ -2054,13 +2054,13 @@ def scraamp(T_A, m, T_B, percentage, exclusion_zone, pre_scraamp, s, p=2.0, k=1)
     l = n_A - m + 1
 
     if exclusion_zone is not None:
-        diags = np.random.permutation(range(exclusion_zone + 1, n_A - m + 1)).astype(
-            np.int64
-        )
+        diags = config.RNG.permutation(
+            np.arange(exclusion_zone + 1, n_A - m + 1)
+        ).astype(np.int64)
     else:
-        diags = np.random.permutation(range(-(n_A - m + 1) + 1, n_B - m + 1)).astype(
-            np.int64
-        )
+        diags = config.RNG.permutation(
+            np.arange(-(n_A - m + 1) + 1, n_B - m + 1)
+        ).astype(np.int64)
 
     n_chunks = int(np.ceil(1.0 / percentage))
     ndist_counts = core._count_diagonal_ndist(diags, m, n_A, n_B)