Revert modern NumPy RNG migration (wontfix #1112) and keep STUMPY_SEED reproducibility (#707)

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.default_rng().random(10000)
+        your_time_series = np.random.rand(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.default_rng().random(10000)
+            your_time_series = np.random.rand(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.default_rng().random(10000)
+        your_time_series = np.random.rand(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.default_rng().random((3, 1000))  # Each row represents data from a different dimension while each column represents data from the same dimension
+        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
         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.default_rng().random((3, 1000))   # Each row represents data from a different dimension while each column represents data from the same dimension
+            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
             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.default_rng().random(10000)
+        your_time_series = np.random.rand(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.default_rng().random(10000)
+        your_time_series = np.random.rand(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.default_rng().random(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.rand(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

@@ -16,6 +16,9 @@
 if "STUMPY_SEED" not in os.environ:
     os.environ["STUMPY_SEED"] = str(np.random.default_rng().integers(2**31))
 
+# Seed the legacy np.random globally for all tests:
+np.random.seed(int(os.environ["STUMPY_SEED"]))
+
 
 def pytest_configure(config):
     """Print STUMPY_SEED so failed tests can be reproduced."""

stumpy/config.py

@@ -40,7 +40,6 @@
 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,8 +3717,7 @@ def check_ignore_trivial(T_A, T_B, ignore_trivial):
     import numpy as np
     import warnings
 
-    rng = np.random.default_rng()
-    T = rng.random(10_000)
+    T = np.random.rand(10_000)
     m = 50
     with warnings.catch_warnings():
         warnings.filterwarnings("ignore", message="Arrays T_A, T_B are equal")
@@ -4494,8 +4493,7 @@ def check_self_join(ignore_trivial):
     import numpy as np
     import warnings
 
-    rng = np.random.default_rng()
-    T = rng.random(10_000)
+    T = np.random.rand(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)
     """
-    rng = np.random.default_rng(seed)
+    np.random.seed(seed)
 
-    I = rng.integers(0, width, size=width, dtype=np.int64)
+    I = np.random.randint(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] = rng.integers(i + 1, width, dtype=np.int64)
+            I[i] = np.random.randint(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, random_state=rng)
+        data = hist_dist.rvs(size=n_samples)
         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 = config.RNG.permutation(np.arange(0, l, s)).astype(np.int64)
+    indices = np.random.permutation(range(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 = config.RNG.permutation(
-                np.arange(self._excl_zone + 1, self._n_A - self._m + 1)
+            self._diags = np.random.permutation(
+                range(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 = config.RNG.permutation(
-                np.arange(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
+            self._diags = np.random.permutation(
+                range(-(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 = config.RNG.permutation(np.arange(0, l, s)).astype(np.int64)
+    indices = np.random.permutation(range(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 = config.RNG.permutation(
-                np.arange(self._excl_zone + 1, self._n_A - self._m + 1)
+            self._diags = np.random.permutation(
+                range(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 = config.RNG.permutation(
-                np.arange(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
+            self._diags = np.random.permutation(
+                range(-(self._n_A - self._m + 1) + 1, self._n_B - self._m + 1)
             ).astype(np.int64)
 
         self._n_threads = numba.config.NUMBA_NUM_THREADS

tests/conftest.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 config.RNG.permutation(np.arange(0, l, s)):
+    for i in np.random.permutation(range(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 = config.RNG.permutation(
-            np.arange(exclusion_zone + 1, n_A - m + 1)
-        ).astype(np.int64)
+        diags = np.random.permutation(range(exclusion_zone + 1, n_A - m + 1)).astype(
+            np.int64
+        )
     else:
-        diags = config.RNG.permutation(
-            np.arange(-(n_A - m + 1) + 1, n_B - m + 1)
-        ).astype(np.int64)
+        diags = np.random.permutation(range(-(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 config.RNG.permutation(np.arange(0, l, s)):
+    for i in np.random.permutation(range(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 = config.RNG.permutation(
-            np.arange(exclusion_zone + 1, n_A - m + 1)
-        ).astype(np.int64)
+        diags = np.random.permutation(range(exclusion_zone + 1, n_A - m + 1)).astype(
+            np.int64
+        )
     else:
-        diags = config.RNG.permutation(
-            np.arange(-(n_A - m + 1) + 1, n_B - m + 1)
-        ).astype(np.int64)
+        diags = np.random.permutation(range(-(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)

tests/naive.py

@@ -13,8 +13,8 @@
         np.array([584, -11, 23, 79, 1001, 0, -19], dtype=np.float64),
     ),
     (
-        config.RNG.uniform(-1000, 1000, [8]).astype(np.float64),
-        config.RNG.uniform(-1000, 1000, [64]).astype(np.float64),
+        np.random.uniform(-1000, 1000, [8]).astype(np.float64),
+        np.random.uniform(-1000, 1000, [64]).astype(np.float64),
     ),
 ]
 
@@ -72,7 +72,7 @@ def test_aamp_constant_subsequence_self_join():
 
 
 def test_aamp_one_constant_subsequence_A_B_join():
-    T_A = config.RNG.random(20)
+    T_A = np.random.rand(20)
     T_B = np.concatenate((np.zeros(20, dtype=np.float64), np.ones(5, dtype=np.float64)))
     m = 3
     ref_mp = naive.aamp(T_A, m, T_B=T_B)
@@ -122,8 +122,8 @@ def test_aamp_two_constant_subsequences_A_B_join():
 
 
 def test_aamp_identical_subsequence_self_join():
-    identical = config.RNG.random(8)
-    T_A = config.RNG.random(20)
+    identical = np.random.rand(8)
+    T_A = np.random.rand(20)
     T_A[1 : 1 + identical.shape[0]] = identical
     T_A[11 : 11 + identical.shape[0]] = identical
     m = 3
@@ -143,9 +143,9 @@ def test_aamp_identical_subsequence_self_join():
 
 
 def test_aamp_identical_subsequence_A_B_join():
-    identical = config.RNG.random(8)
-    T_A = config.RNG.random(20)
-    T_B = config.RNG.random(20)
+    identical = np.random.rand(8)
+    T_A = np.random.rand(20)
+    T_B = np.random.rand(20)
     T_A[1 : 1 + identical.shape[0]] = identical
     T_B[11 : 11 + identical.shape[0]] = identical
     m = 3