Merge pull request #1708 from cta-observatory/toymodel_rng

Use new random API

ctapipe/calib/camera/tests/test_calibrator.py

@@ -154,7 +154,7 @@ def test_dl1_charge_calib(example_subarray):
     n_samples = 96
     mid = n_samples // 2
     pulse_sigma = 6
-    random = np.random.RandomState(1)
+    random = np.random.default_rng(1)
     x = np.arange(n_samples)
 
     # Randomize times and create pulses

ctapipe/image/muon/tests/test_muon_features.py

@@ -13,31 +13,32 @@ def test_ring_containment():
     ring_center_x = 0 * u.deg
     ring_center_y = 0 * u.deg
     containment = ring_containment(
-        ring_radius, ring_center_x, ring_center_y, cam_radius,
+        ring_radius, ring_center_x, ring_center_y, cam_radius
     )
     assert containment == 1.0
 
     ring_center_x = 0 * u.deg
     ring_center_y = cam_radius
     containment = ring_containment(
-        ring_radius, ring_center_x, ring_center_y, cam_radius,
+        ring_radius, ring_center_x, ring_center_y, cam_radius
     )
     assert 0.4 <= containment <= 0.5
 
     ring_center_x = 0 * u.deg
     ring_center_y = cam_radius + 1.1 * ring_radius
     containment = ring_containment(
-        ring_radius, ring_center_x, ring_center_y, cam_radius,
+        ring_radius, ring_center_x, ring_center_y, cam_radius
     )
     assert containment == 0.0
 
 
 def test_ring_completeness():
+    rng = np.random.default_rng(0)
 
     angle_ring = np.linspace(0, 2 * math.pi, 360)
     x = np.cos(angle_ring) * u.m
     y = np.sin(angle_ring) * u.m
-    pe = np.random.uniform(0, 100, len(x))
+    pe = rng.uniform(0, 100, len(x))
     ring_radius = 1.0 * u.m
 
     ring_center_x = 0.0 * u.m

ctapipe/image/muon/tests/test_muon_fitting.py

@@ -3,19 +3,18 @@
 
 from ctapipe.image.muon import kundu_chaudhuri_circle_fit
 
-np.random.seed(0)
-
 
 def test_kundu_chaudhuri():
 
     num_tests = 10
-    center_xs = np.random.uniform(-1000, 1000, num_tests)
-    center_ys = np.random.uniform(-1000, 1000, num_tests)
-    radii = np.random.uniform(10, 1000, num_tests)
+    rng = np.random.default_rng(0)
+    center_xs = rng.uniform(-1000, 1000, num_tests)
+    center_ys = rng.uniform(-1000, 1000, num_tests)
+    radii = rng.uniform(10, 1000, num_tests)
 
     for center_x, center_y, radius in zip(center_xs, center_ys, radii):
 
-        phi = np.random.uniform(0, 2 * np.pi, 100)
+        phi = rng.uniform(0, 2 * np.pi, 100)
         x = center_x + radius * np.cos(phi)
         y = center_y + radius * np.sin(phi)
 
@@ -34,7 +33,8 @@ def test_kundu_chaudhuri_with_units():
     center_y = 0.5 * u.meter
     radius = 1 * u.meter
 
-    phi = np.random.uniform(0, 2 * np.pi, 100)
+    rng = np.random.default_rng(0)
+    phi = rng.uniform(0, 2 * np.pi, 100)
     x = center_x + radius * np.cos(phi)
     y = center_y + radius * np.sin(phi)
 

ctapipe/image/tests/test_hillas.py

@@ -6,7 +6,6 @@
 from astropy import units as u
 import numpy as np
 from numpy import isclose, zeros_like
-from numpy.random import seed
 from pytest import approx
 import itertools
 import pytest
@@ -20,19 +19,15 @@ def create_sample_image(
     length=0.15 * u.m,
     intensity=1500,
 ):
-    seed(10)
 
     geom = CameraGeometry.from_name("LSTCam")
 
     # make a toymodel shower model
     model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)
 
     # generate toymodel image in camera for this shower model.
-    image, _, _ = model.generate_image(
-        geom,
-        intensity=1500,
-        nsb_level_pe=3,
-    )
+    rng = np.random.default_rng(0)
+    image, _, _ = model.generate_image(geom, intensity=1500, nsb_level_pe=3, rng=rng)
 
     # calculate pixels likely containing signal
     clean_mask = tailcuts_clean(geom, image, 10, 5)
@@ -113,7 +108,7 @@ def test_hillas_container():
 
 
 def test_with_toy():
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
 
     geom = CameraGeometry.from_name("LSTCam")
 
@@ -131,18 +126,10 @@ def test_with_toy():
         for psi in psis:
 
             # make a toymodel shower model
-            model = toymodel.Gaussian(
-                x=x,
-                y=y,
-                width=width,
-                length=length,
-                psi=psi,
-            )
+            model = toymodel.Gaussian(x=x, y=y, width=width, length=length, psi=psi)
 
             image, signal, noise = model.generate_image(
-                geom,
-                intensity=intensity,
-                nsb_level_pe=5,
+                geom, intensity=intensity, nsb_level_pe=5, rng=rng
             )
 
             result = hillas_parameters(geom, signal)
@@ -162,7 +149,7 @@ def test_with_toy():
 
 
 def test_skewness():
-    np.random.seed(42)
+    rng = np.random.default_rng(42)
 
     geom = CameraGeometry.from_name("LSTCam")
 
@@ -178,18 +165,11 @@ def test_skewness():
     for x, y, psi, skew in itertools.product(xs, ys, psis, skews):
         # make a toymodel shower model
         model = toymodel.SkewedGaussian(
-            x=x,
-            y=y,
-            width=width,
-            length=length,
-            psi=psi,
-            skewness=skew,
+            x=x, y=y, width=width, length=length, psi=psi, skewness=skew
         )
 
         _, signal, _ = model.generate_image(
-            geom,
-            intensity=intensity,
-            nsb_level_pe=5,
+            geom, intensity=intensity, nsb_level_pe=5, rng=rng
         )
 
         result = hillas_parameters(geom, signal)
@@ -222,7 +202,7 @@ def test_straight_line_width_0():
     trans = np.zeros(len(long))
     pix_id = np.arange(len(long))
 
-    np.random.seed(0)
+    rng = np.random.default_rng(0)
 
     for dx in (-1, 0, 1):
         for dy in (-1, 0, 1):
@@ -239,7 +219,7 @@ def test_straight_line_width_0():
                     pix_area=1 * u.m ** 2,
                 )
 
-                img = np.random.poisson(5, size=len(long))
+                img = rng.poisson(5, size=len(long))
                 result = hillas_parameters(geom, img)
                 assert result.width.value == 0
                 assert np.isnan(result.width_uncertainty.value)

ctapipe/image/tests/test_reducer.py

@@ -28,7 +28,8 @@ def subarray_lst():
 
 def test_null_data_volume_reducer(subarray_lst):
     subarray, _, _, _, _ = subarray_lst
-    waveforms = np.random.uniform(0, 1, (2048, 96))
+    rng = np.random.default_rng(0)
+    waveforms = rng.uniform(0, 1, (2048, 96))
     reducer = NullDataVolumeReducer(subarray=subarray)
     reduced_waveforms_mask = reducer(waveforms)
     reduced_waveforms = waveforms.copy()

ctapipe/image/tests/test_statistics.py

@@ -5,8 +5,8 @@
 def test_statistics():
     from ctapipe.image import descriptive_statistics
 
-    np.random.seed(0)
-    data = np.random.normal(5, 2, 1000)
+    rng = np.random.default_rng(0)
+    data = rng.normal(5, 2, 1000)
 
     stats = descriptive_statistics(data)
 
@@ -19,8 +19,8 @@ def test_statistics():
 def test_skewness():
     from ctapipe.image.statistics import skewness
 
-    np.random.seed(0)
-    data = np.random.normal(5, 2, 1000)
+    rng = np.random.default_rng(0)
+    data = rng.normal(5, 2, 1000)
     mean = np.mean(data)
     std = np.std(data)
 
@@ -33,8 +33,8 @@ def test_skewness():
 def test_kurtosis():
     from ctapipe.image.statistics import kurtosis
 
-    np.random.seed(0)
-    data = np.random.normal(5, 2, 1000)
+    rng = np.random.default_rng(0)
+    data = rng.normal(5, 2, 1000)
 
     mean = np.mean(data)
     std = np.std(data)
@@ -52,8 +52,8 @@ def test_return_type():
     from ctapipe.containers import PeakTimeStatisticsContainer, StatisticsContainer
     from ctapipe.image import descriptive_statistics
 
-    np.random.seed(0)
-    data = np.random.normal(5, 2, 1000)
+    rng = np.random.default_rng(0)
+    data = rng.normal(5, 2, 1000)
 
     stats = descriptive_statistics(data)
     assert isinstance(stats, StatisticsContainer)

ctapipe/image/tests/test_timing_parameters.py

@@ -19,7 +19,7 @@ def test_psi_0():
     geom = CameraGeometry.from_name("LSTCam")
     hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg)
 
-    random = np.random.RandomState(1)
+    random = np.random.default_rng(0)
     peak_time = intercept + grad * geom.pix_x.value
     peak_time += random.normal(0, deviation, geom.n_pixels)
 
@@ -49,7 +49,7 @@ def test_psi_20():
     psi = 20 * u.deg
     hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=psi)
 
-    random = np.random.RandomState(1)
+    random = np.random.default_rng(0)
     peak_time = intercept + grad * (
         np.cos(psi) * geom.pix_x.value + np.sin(psi) * geom.pix_y.value
     )
@@ -79,7 +79,7 @@ def test_ignore_negative():
     geom = CameraGeometry.from_name("LSTCam")
     hillas = HillasParametersContainer(x=0 * u.m, y=0 * u.m, psi=0 * u.deg)
 
-    random = np.random.RandomState(1)
+    random = np.random.default_rng(0)
     peak_time = intercept + grad * geom.pix_x.value
     peak_time += random.normal(0, deviation, geom.n_pixels)
 

ctapipe/image/tests/test_toy.py

@@ -4,12 +4,13 @@
 from pytest import approx
 from scipy.stats import poisson, skewnorm, norm
 import astropy.units as u
+import pytest
 
 
-def test_intensity():
+@pytest.mark.parametrize("seed", [None, 0])
+def test_intensity(seed):
     from ctapipe.image.toymodel import Gaussian
 
-    np.random.seed(0)
     geom = CameraGeometry.from_name("LSTCam")
 
     x, y = u.Quantity([0.2, 0.3], u.m)
@@ -21,7 +22,13 @@ def test_intensity():
     # make a toymodel shower model
     model = Gaussian(x=x, y=y, width=width, length=length, psi=psi)
 
-    _, signal, _ = model.generate_image(geom, intensity=intensity, nsb_level_pe=5)
+    if seed is None:
+        _, signal, _ = model.generate_image(geom, intensity=intensity, nsb_level_pe=5)
+    else:
+        rng = np.random.default_rng(seed)
+        _, signal, _ = model.generate_image(
+            geom, intensity=intensity, nsb_level_pe=5, rng=rng
+        )
 
     # test if signal reproduces given cog values
     assert np.average(geom.pix_x.to_value(u.m), weights=signal) == approx(0.2, rel=0.15)
@@ -43,7 +50,7 @@ def test_skewed():
 
     # test if the parameters we calculated for the skew normal
     # distribution produce the correct moments
-    np.random.seed(0)
+    rng = np.random.default_rng(0)
     geom = CameraGeometry.from_name("LSTCam")
 
     x, y = u.Quantity([0.2, 0.3], u.m)
@@ -56,7 +63,7 @@ def test_skewed():
     model = SkewedGaussian(
         x=x, y=y, width=width, length=length, psi=psi, skewness=skewness
     )
-    model.generate_image(geom, intensity=intensity, nsb_level_pe=5)
+    model.generate_image(geom, intensity=intensity, nsb_level_pe=5, rng=rng)
 
     a, loc, scale = model._moments_to_parameters()
     mean, var, skew = skewnorm(a=a, loc=loc, scale=scale).stats(moments="mvs")

ctapipe/image/toymodel.py

@@ -23,6 +23,7 @@
 from scipy.stats import multivariate_normal, skewnorm, norm
 from scipy.ndimage import convolve1d
 from abc import ABCMeta, abstractmethod
+from numpy.random import default_rng
 
 __all__ = [
     "WaveformModel",
@@ -33,6 +34,9 @@
 ]
 
 
+TOYMODEL_RNG = default_rng(0)
+
+
 @u.quantity_input(
     x=u.m,
     y=u.m,
@@ -182,7 +186,7 @@ def pdf(self, x, y):
         """Probability density function.
         """
 
-    def generate_image(self, camera, intensity=50, nsb_level_pe=20):
+    def generate_image(self, camera, intensity=50, nsb_level_pe=20, rng=None):
         """Generate a randomized DL1 shower image.
         For the signal, poisson random numbers are drawn from
         the expected signal distribution for each pixel.
@@ -205,10 +209,13 @@ def generate_image(self, camera, intensity=50, nsb_level_pe=20):
         noise: only the noise part of image
 
         """
+        if rng is None:
+            rng = TOYMODEL_RNG
+
         expected_signal = self.expected_signal(camera, intensity)
 
-        signal = np.random.poisson(expected_signal)
-        noise = np.random.poisson(nsb_level_pe, size=signal.shape)
+        signal = rng.poisson(expected_signal)
+        noise = rng.poisson(nsb_level_pe, size=signal.shape)
         image = (signal + noise) - np.mean(noise)
 
         return image, signal, noise

ctapipe/instrument/tests/test_subarray.py

@@ -15,14 +15,16 @@
 
 def example_subarray(n_tels=10):
     """ generate a simple subarray for testing purposes """
+    rng = np.random.default_rng(0)
+
     pos = {}
     tel = {}
 
     for tel_id in range(1, n_tels + 1):
         tel[tel_id] = TelescopeDescription.from_name(
             optics_name="MST", camera_name="NectarCam"
         )
-        pos[tel_id] = np.random.uniform(-100, 100, size=3) * u.m
+        pos[tel_id] = rng.uniform(-100, 100, size=3) * u.m
 
     return SubarrayDescription("test array", tel_positions=pos, tel_descriptions=tel)