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

[seanlaw]

See issue #1112

This is a second attempt that follows failed PR #1131

Pull Request Checklist

Below is a simple checklist but please do not hesitate to ask for assistance!

• Fork, clone, and checkout the newest version of the code
• Create a new branch
• Make necessary code changes
• Install black (i.e., python -m pip install black or conda install -c conda-forge black)
• Install flake8 (i.e., python -m pip install flake8 or conda install -c conda-forge flake8)
• Install pytest-cov (i.e., python -m pip install pytest-cov or conda install -c conda-forge pytest-cov)
• Run black --exclude=".*\.ipynb" --extend-exclude=".venv" --diff ./ in the root stumpy directory
• Run flake8 --extend-exclude=.venv ./ in the root stumpy directory
• Run ./setup.sh dev && ./test.sh in the root stumpy directory
• Reference a Github issue (and create one if one doesn't already exist)

To Do List

1. Search for "np.random" in all stumpy/*.py and tests/test_*py files
2. Replace np.random.rand(x) with rng.RNG.random(x)
3. Replace np.random.rand(x*y).reshape(x,y) with rng.RNG.random(size=(x, y))
4. Replace np.random.uniform(x, y, [z]) with rng.RNG.uniform(x, y, size=z)
5. Replace np.random.permutation([x, y, z]) with rng.RNG.permutation([x, y, z])
6. Replace np.random.randint(x, y, z) with rng.RNG.integers(x, y, z)
7. Replace np.random.choice([x, y], l, replace=True) with rng.RNG.choice([x, y], l, replace=True)
8. Replace np.random.normal(size=x) with rng.RNG.normal(size=x)
9. Replace all np.random from docstrings, README, tutorials, etc
10. Fix tests/test_non_normalized.py
11. Search for"seed" in all stumpy/*.py and tests/test_*py files
12. Run unit tests on Google Colab

[gitnotebooks]

Review these changes at https://app.gitnotebooks.com/stumpy-dev/stumpy/pull/1147

[seanlaw]

@NimaSarajpoor I've created a newstumpy/rng.py module that will manage our random number usage and tested the numpy.random migration on test_core.py as an initial trial run. Would you mind taking a look and providing any feedback before I do the same thing on the other tests?

I will be adding a tests/test_rng.py next

[NimaSarajpoor]

@seanlaw Please allow me some time to take a look. Will provide an update in a couple of days.

[seanlaw]

This function in tests/test_scrump.py is going to need more attention in order to migrate:

stumpy/tests/test_scrump.py

Lines 963 to 1065 in 86145c2

	def test_prescrump_self_join_KNN_no_overlap():
	# This test is particularly designed to raise error in a rare case described
	# as follows: Let's denote `I[i]` as the array with length `k` that contains
	# the start indices of the best-so-far top-k nearest neighbors of `subseq i`,
	# (`S_i`). Also, we denote `P[i]` as their corresponding distances sorted in
	# ascending order. Let's denote `d` as the distance between `S_i` and `S_j`. P[i]
	# and I[i] must be updated if (1) `j` is not in I[i] and (2) `d` < P[i,-1].
	# Regarding the former condition, one needs to check the whole array I[i]. Checking
	# the array I[i, :idx], where `idx = np.searchsorted(P[i], 'd', side='right')` is
	# not completely correct and that is due to imprecision in numerical calculation.
	# It may happen that `j` is not in `I[i, :idx]`, but it is in fact at `I[i, idx]`
	# (or any other position in array I[i]). And, its corresponding distance, i.e
	# P[i, idx], is d + 1e-5, for instance. In theory, this should be exactly `d`.
	# However, due to imprecision, we may calculated a slightly different value
	# for such distance in one of previous iterations in function prescrump. This
	# test results in error if someone tries to change the performant code of prescrump
	# function and check `I[i, :idx]` rather than the full array `I[i]`.
	T = np.array(
	[
	-916.64703784,
	-327.42056679,
	379.19386284,
	-281.80427628,
	-189.85401773,
	-38.69610569,
	187.89889345,
	578.65862523,
	528.09687811,
	-667.42973795,
	-285.27749324,
	-211.28930925,
	-703.93802657,
	-820.53780562,
	-955.91174663,
	383.65471851,
	932.08809422,
	-563.57569746,
	784.0546579,
	-343.14886064,
	-612.72329848,
	-270.09273091,
	-448.39346549,
	578.03202014,
	867.15436674,
	-783.55167049,
	-494.78062922,
	-311.18567747,
	522.70052256,
	933.45474094,
	192.34822368,
	-162.11374908,
	-612.95359279,
	-449.62297051,
	-351.79138459,
	-77.70189101,
	-439.46519487,
	-660.48431174,
	548.69362177,
	485.36004744,
	-535.3566627,
	-568.0955257,
	755.26647273,
	736.1079588,
	-597.65672557,
	379.3299783,
	731.38211912,
	247.34827447,
	545.41888454,
	644.94300763,
	20.99042666,
	788.19859515,
	-898.24325898,
	-929.47841134,
	-738.45875181,
	66.01030291,
	512.945841,
	-44.07720164,
	302.97141464,
	-696.95271302,
	662.98385163,
	-712.3807531,
	-43.62688539,
	74.16927482,
	]
	)
	# test_cases: dict() with `key: value` pair, where key is `(m, k)`, and value
	# is a list of random `seeds`
	test_cases = {
	(3, 2): [4279, 9133, 8190],
	(3, 5): [1267, 4016, 4046],
	(5, 2): [6327, 4926, 3712],
	(5, 5): [3032, 3032, 8117],
	}
	for (m, k), specified_seeds in test_cases.items():
	zone = int(np.ceil(m / 4))
	for seed in specified_seeds:
	np.random.seed(seed)
	ref_P, ref_I = naive.prescrump(T, m, T, s=1, exclusion_zone=zone, k=k)
	comp_P, comp_I = prescrump(T, m, s=1, k=k)
	npt.assert_almost_equal(ref_P, comp_P)
	npt.assert_almost_equal(ref_I, comp_I)

[seanlaw]

The function get_central_motif (and also get_aamp_central_motif) in naive.py is missing coverage in some occasions. However, when the state is set to:

RNG.bit_generator.state = {
    "bit_generator": "PCG64",
    "state": {
        "state": 62474265484091438068208639126596717761,
        "inc": 18021712563375961045953296450069431849
    },
    "has_uint32": 0,
    "uinteger": 0
}

This allows it to be covered (likely by tests/test_ostinato.py)

Update

This was resolved by setting an RNG state that passes through the missing coverage code

[seanlaw]

tests/test_precision.py has a lot of seeds set and causing missing coverage

Update
Resolved by hardcoding input array T

[seanlaw]

The _iac function in stumpy/floss.py sets a random seed AND the scipy.stats.rv_histogram uses a random number generator:

stumpy/stumpy/floss.py

Lines 87 to 106 in 86145c2

	np.random.seed(seed)
	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] = np.random.randint(i + 1, width, dtype=np.int64)
	target_AC = _nnmark(I)
	params = np.empty((n_iter, 2), dtype=np.float64)
	for i in range(n_iter):
	hist_dist = scipy.stats.rv_histogram(
	(target_AC, np.append(np.arange(width), width))
	)
	data = hist_dist.rvs(size=n_samples)
	a, b, c, d = scipy.stats.beta.fit(data, floc=0, fscale=width)
	params[i, 0] = a
	params[i, 1] = b

Update

Luckily, the rvs method inside of scipy.stats.rv_histogram supports a random_state parameter and you can simply pass in your own random_state=rng and this resolves the issue of reproducibility