Fixed #1112 Migrate to New NumPy Random Number Generator (RNG) API

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,15 @@
 # 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 json
+
+from stumpy import rng
+
+
+def pytest_configure(config):
+    """
+    Called after command line options have been parsed
+    and all plugins and initial conftest files been loaded.
+    """
+    state = json.dumps(rng.RNG.bit_generator.state, indent=4)
+    print(f"stumpy/rng.py: RNG.bit_generator.state = {state}")

docs/Tutorial_Fast_Approximate_Matrix_Profiles.ipynb

@@ -334,11 +334,14 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Please keep in mind that the approximate matrix profile is computed by randomly computing distances along a subset of diagonals. So, each time you initialize a new `scrump` object by calling `stumpy.scrump`, this will randomly shuffle the order that the distances are computed, which inevitably results in different approximate matrix profiles (except when `percentage=1.0`). Depending on your use case, to ensure reproducible results, you may consider setting random seed prior to calling `stumpy.scrump`:\n",
+    "Please keep in mind that the approximate matrix profile is computed by randomly computing distances along a subset of diagonals. So, each time you initialize a new `scrump` object by calling `stumpy.scrump`, this will randomly shuffle the order that the distances are computed, which inevitably results in different approximate matrix profiles (except when `percentage=1.0`). Depending on your use case, to ensure reproducible results, you may consider fixing the random state prior to calling `stumpy.scrump`:\n",
     "\n",
     "```\n",
-    "seed = np.random.randint(100000)\n",
-    "np.random.seed(seed)\n",
+    "# Get and set random state\n",
+    "from stumpy import rng\n",
+    "state = rng.get_state()\n",
+    "rng.set_state(state)\n",
+    "\n",
     "approx = stumpy.scrump(steam_df['steam flow'], m, percentage=0.01, pre_scrump=False)\n",
     "```\n",
     "\n",
@@ -503,7 +506,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.3"
+   "version": "3.14.6"
   }
  },
  "nbformat": 4,

docs/Tutorial_Matrix_Profiles_For_Streaming_Data.ipynb

@@ -47,7 +47,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "T = np.random.rand(336)"
+    "rng = np.random.default_rng()\n",
+    "T = rng.random(336)"
    ]
   },
   {
@@ -122,7 +123,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "t = np.random.rand()"
+    "t = rng.random()"
    ]
   },
   {
@@ -155,7 +156,7 @@
    "outputs": [],
    "source": [
     "for i in range(1000):\n",
-    "    t = np.random.rand()\n",
+    "    t = rng.random()\n",
     "    stream.update(t)"
    ]
   },
@@ -183,7 +184,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "T_full = np.random.rand(64)\n",
+    "T_full = rng.random(64)\n",
     "m = 8\n",
     "\n",
     "mp = stumpy.stump(T_full, m)\n",
@@ -280,7 +281,7 @@
     }
    ],
    "source": [
-    "T_full = np.random.rand(1_000)\n",
+    "T_full = rng.random(1_000)\n",
     "T_stream = T_full[:20].copy()\n",
     "m = 10\n",
     "\n",
@@ -266241,7 +266242,7 @@
     }
    ],
    "source": [
-    "T_full = np.random.rand(64)\n",
+    "T_full = rng.random(64)\n",
     "m = 8\n",
     "\n",
     "T_stream = T_full[:10].copy()\n",
@@ -266302,7 +266303,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.4"
+   "version": "3.14.6"
   }
  },
  "nbformat": 4,

docs/Tutorial_Multidimensional_Motif_Discovery.ipynb

@@ -571,7 +571,7 @@
    ],
    "source": [
     "for i in range(4, 11):\n",
-    "    df[f'T{i}'] = np.random.uniform(0.1, -0.1, size=df.shape[0]).cumsum()\n",
+    "    df[f'T{i}'] = np.random.default_rng().uniform(0.1, -0.1, size=df.shape[0]).cumsum()\n",
     "    \n",
     "df = df.sample(frac=1, axis=\"columns\")  # Randomly shuffle the columns\n",
     "\n",
@@ -804,7 +804,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.4"
+   "version": "3.14.6"
   },
   "widgets": {
    "application/vnd.jupyter.widget-state+json": {

docs/Tutorial_Pattern_Matching.ipynb

@@ -645,8 +645,9 @@
     }
    ],
    "source": [
-    "Q_random = np.random.rand(100)\n",
-    "T_random = np.random.rand(1_000_000)\n",
+    "rng = np.random.default_rng()\n",
+    "Q_random = rng.random(100)\n",
+    "T_random = rng.random(1_000_000)\n",
     "\n",
     "naive_distance_profile = compute_naive_distance_profile(Q_random, T_random)"
    ]
@@ -736,7 +737,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.4"
+   "version": "3.14.6"
   }
  },
  "nbformat": 4,

docs/Tutorial_Time_Series_Chains.ipynb

@@ -70,11 +70,12 @@
     }
    ],
    "source": [
-    "x = np.random.rand(20)\n",
-    "y = np.random.rand(20)\n",
+    "rng = np.random.default_rng()\n",
+    "x = rng.random(20)\n",
+    "y = rng.random(20)\n",
     "n = 10\n",
-    "motifs_x = 0.5 * np.ones(n) + np.random.uniform(-0.05, 0.05, n)\n",
-    "motifs_y = 0.5 * np.ones(n) + np.random.uniform(-0.05, 0.05, n)\n",
+    "motifs_x = 0.5 * np.ones(n) + rng.uniform(-0.05, 0.05, n)\n",
+    "motifs_y = 0.5 * np.ones(n) + rng.uniform(-0.05, 0.05, n)\n",
     "sin_x = np.linspace(0, np.pi/2, n+1)\n",
     "sin_y = np.sin(sin_x)/4\n",
     "chains_x = 0.5 * np.ones(n+1) + 0.02 * np.arange(n+1)\n",
@@ -605,7 +606,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.4"
+   "version": "3.14.6"
   }
  },
  "nbformat": 4,

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)
+    local_rng = np.random.default_rng(seed)
 
-    I = np.random.randint(0, width, size=width, dtype=np.int64)
+    I = local_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] = local_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=local_rng)
         a, b, c, d = scipy.stats.beta.fit(data, floc=0, fscale=width)
 
         params[i, 0] = a