diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f646c2c..87d4aa6 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -1,6 +1,17 @@ Changelog ========= +:Version: 2.0.0 of 2022-04-16 + +Allow for control over random number generator state by adding optional random_state +argument thanks to contributions from i-aki-y. +Drop Python 2.7 support to use of NumPy's recommended default_rng constructor. + +:Version: 1.2.0 of 2022-03-13 + +Improve doc strings based on user questions. +Check that fmin parameter is in the right range. + :Version: 1.1.1 of 2019-02-08 Use numpy's sum instead of python's (thank's to RuABraun). diff --git a/README.rst b/README.rst index 920d9b2..8ac90e4 100644 --- a/README.rst +++ b/README.rst @@ -23,14 +23,15 @@ Installation pip install colorednoise - + Dependencies ------------ - - Python >= 2.7 or >= 3.6 - - NumPy + - Python >= 3.6.15 + - NumPy >= 1.17.0 -Other Python versions were not tested, but are likely to work. +Older Python 3 versions were not tested, but are likely to work. +For Python 2 please use colorednoise version 1.x. Examples diff --git a/colorednoise.py b/colorednoise.py index fae2392..ee99a22 100644 --- a/colorednoise.py +++ b/colorednoise.py @@ -1,12 +1,12 @@ """Generate colored noise.""" -from numpy import sqrt, newaxis +from numpy import sqrt, newaxis, integer from numpy.fft import irfft, rfftfreq -from numpy.random import normal +from numpy.random import default_rng, Generator, RandomState from numpy import sum as npsum -def powerlaw_psd_gaussian(exponent, size, fmin=0): +def powerlaw_psd_gaussian(exponent, size, fmin=0, random_state=None): """Gaussian (1/f)**beta noise. Based on the algorithm in: @@ -46,6 +46,12 @@ def powerlaw_psd_gaussian(exponent, size, fmin=0): sample. The largest possible value is fmin = 0.5, the Nyquist frequency. The output for this value is white noise. + random_state : int, numpy.integer, numpy.random.Generator, numpy.random.RandomState, + optional + Optionally sets the state of NumPy's underlying random number generator. + Integer-compatible values or None are passed to np.random.default_rng. + np.random.RandomState or np.random.Generator are used directly. + Default: None. Returns ------- @@ -100,9 +106,12 @@ def powerlaw_psd_gaussian(exponent, size, fmin=0): dims_to_add = len(size) - 1 s_scale = s_scale[(newaxis,) * dims_to_add + (Ellipsis,)] + # prepare random number generator + normal_dist = _get_normal_distribution(random_state) + # Generate scaled random power + phase - sr = normal(scale=s_scale, size=size) - si = normal(scale=s_scale, size=size) + sr = normal_dist(scale=s_scale, size=size) + si = normal_dist(scale=s_scale, size=size) # If the signal length is even, frequencies +/- 0.5 are equal # so the coefficient must be real. @@ -118,3 +127,18 @@ def powerlaw_psd_gaussian(exponent, size, fmin=0): y = irfft(s, n=samples, axis=-1) / sigma return y + + +def _get_normal_distribution(random_state): + normal_dist = None + if isinstance(random_state, (integer, int)) or random_state is None: + random_state = default_rng(random_state) + normal_dist = random_state.normal + elif isinstance(random_state, (Generator, RandomState)): + normal_dist = random_state.normal + else: + raise ValueError( + "random_state must be one of integer, numpy.random.Generator, " + "numpy.random.Randomstate" + ) + return normal_dist diff --git a/setup.py b/setup.py index 3aaffde..0144504 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ setup( name='colorednoise', - version='1.2.0', + version='2.0.0', description='Generate Gaussian (1/f)**beta noise (e.g. pink noise)', long_description=""" Generate Gaussian distributed noise with a power law spectrum. @@ -17,10 +17,12 @@ 'Intended Audience :: Education', 'Intended Audience :: Science/Research', 'License :: OSI Approved :: MIT License', - 'Programming Language :: Python :: 2', - 'Programming Language :: Python :: 2.7', 'Programming Language :: Python :: 3', 'Programming Language :: Python :: 3.6', + 'Programming Language :: Python :: 3.7', + 'Programming Language :: Python :: 3.8', + 'Programming Language :: Python :: 3.9', + 'Programming Language :: Python :: 3.10', 'Topic :: Scientific/Engineering', 'Topic :: Software Development :: Libraries', 'Topic :: Software Development :: Libraries :: Python Modules' @@ -35,8 +37,9 @@ license='MIT', py_modules=['colorednoise'], install_requires=[ - 'numpy', + 'numpy>=1.17.0', ], + python_requires='>=3', include_package_data=True, zip_safe=False, test_suite='tests', diff --git a/tests/test_powerlaw_psd_gaussian.py b/tests/test_powerlaw_psd_gaussian.py index 9c70311..00106b4 100644 --- a/tests/test_powerlaw_psd_gaussian.py +++ b/tests/test_powerlaw_psd_gaussian.py @@ -50,3 +50,44 @@ def test_slope_distribution(self): slope_in = (exponent + fit[0] < 3 * np.sqrt(fcov[0,0])).mean() self.assertTrue(slope_in > 0.95) + + def test_random_state_type(self): + exp = 1 + n = 5 + seed = 1 + + good_random_states = [ + np.random.default_rng(seed), + np.random.RandomState(seed), + int(seed), + np.int32(1), + True, + None + ] + for random_state in good_random_states: + cn.powerlaw_psd_gaussian(exp, n, random_state=random_state) + + bad_random_states = ["1", 0.15, [1]] + for random_state in bad_random_states: + self.assertRaises( + ValueError, + cn.powerlaw_psd_gaussian, + exp, + n, + random_state=random_state + ) + + + def test_random_state_reproducibility(self): + exp = 1 + n = 5 + seed = 1 + + rs1 = np.random.default_rng(seed) + rs2 = np.random.default_rng(seed) + + y1 = cn.powerlaw_psd_gaussian(exp, n, random_state=rs1) + np.random.seed(123) + y2 = cn.powerlaw_psd_gaussian(exp, n, random_state=rs2) + + np.testing.assert_array_equal(y1, y2)