Version 2.0 (#11)

* Add support random_state argument

* require python3

* require numpy >= 1.17 for default_rng

Co-authored-by: akiyuki ishikawa <[email protected]>

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).

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

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

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',

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)