Merge pull request #871 from AlexN1234/UpdateRNGCalls839

Update Numpy RNG calls.  Fixes #839

AUTHORS.rst

@@ -47,6 +47,7 @@ Alphabetical list of code contributors
 * James McCormac (@jmccormac01)
 * Abigale Moen (@AbigaleMoen)
 * Stefan Nelson (@stefannelson)
+* Alex Niemi (@AlexN1234)
 * Joe Philip Ninan (@indiajoe)
 * Punyaslok Pattnaik (@Punyaslok)
 * Adrian Price-Whelan (@adrn)

ccdproc/tests/make_mef.py

@@ -1,6 +1,5 @@
 import numpy as np
 from astropy.io import fits
-from astropy.utils.misc import NumpyRNGContext
 
 
 def make_sample_mef(science_name, flat_name, size=10, dtype="float32"):
@@ -23,32 +22,31 @@ def make_sample_mef(science_name, flat_name, size=10, dtype="float32"):
     dtype : str or numpy dtype, optional
         dtype of the generated images.
     """
-    with NumpyRNGContext(1234):
-        number_of_image_extensions = 3
-        science_image = [fits.PrimaryHDU()]
-        flat_image = [fits.PrimaryHDU()]
-        for _ in range(number_of_image_extensions):
-            # Simulate a cloudy night, average pixel
-            # value of 100 with a read_noise of 1 electron.
-            data = np.random.normal(100.0, 1.0, [size, size]).astype(dtype)
-            hdu = fits.ImageHDU(data=data)
-            # Make a header that is at least somewhat realistic
-            hdu.header["unit"] = "electron"
-            hdu.header["object"] = "clouds"
-            hdu.header["exptime"] = 30.0
-            hdu.header["date-obs"] = "1928-07-23T21:03:27"
-            hdu.header["filter"] = "B"
-            hdu.header["imagetyp"] = "LIGHT"
-            science_image.append(hdu)
-
-            # Make a perfect flat
-            flat = np.ones_like(data, dtype=dtype)
-            flat_hdu = fits.ImageHDU(data=flat)
-            flat_hdu.header["unit"] = "electron"
-            flat_hdu.header["filter"] = "B"
-            flat_hdu.header["imagetyp"] = "FLAT"
-            flat_hdu.header["date-obs"] = "1928-07-23T21:03:27"
-            flat_image.append(flat_hdu)
+    number_of_image_extensions = 3
+    science_image = [fits.PrimaryHDU()]
+    flat_image = [fits.PrimaryHDU()]
+    for _ in range(number_of_image_extensions):
+        # Simulate a cloudy night, average pixel
+        # value of 100 with a read_noise of 1 electron.
+        data = np.random.default_rng().normal(100.0, 1.0, [size, size]).astype(dtype)
+        hdu = fits.ImageHDU(data=data)
+        # Make a header that is at least somewhat realistic
+        hdu.header["unit"] = "electron"
+        hdu.header["object"] = "clouds"
+        hdu.header["exptime"] = 30.0
+        hdu.header["date-obs"] = "1928-07-23T21:03:27"
+        hdu.header["filter"] = "B"
+        hdu.header["imagetyp"] = "LIGHT"
+        science_image.append(hdu)
+
+        # Make a perfect flat
+        flat = np.ones_like(data, dtype=dtype)
+        flat_hdu = fits.ImageHDU(data=flat)
+        flat_hdu.header["unit"] = "electron"
+        flat_hdu.header["filter"] = "B"
+        flat_hdu.header["imagetyp"] = "FLAT"
+        flat_hdu.header["date-obs"] = "1928-07-23T21:03:27"
+        flat_image.append(flat_hdu)
 
     science_image = fits.HDUList(science_image)
     science_image.writeto(science_name)

ccdproc/tests/pytest_fixtures.py

@@ -6,7 +6,6 @@
 import pytest
 from astropy import units as u
 from astropy.nddata import CCDData
-from astropy.utils import NumpyRNGContext
 
 from ..utils.sample_directory import directory_for_testing
 
@@ -54,8 +53,10 @@ def ccd_data(
     scale = data_scale
     mean = data_mean
 
-    with NumpyRNGContext(rng_seed):
-        data = np.random.normal(loc=mean, size=[size, size], scale=scale)
+    ##Create random number generator with a specified state
+    rng = np.random.default_rng(seed=rng_seed)
+
+    data = rng.normal(loc=mean, size=[size, size], scale=scale)
 
     fake_meta = {"my_key": 42, "your_key": "not 42"}
     ccd = CCDData(data, unit=u.adu)

ccdproc/tests/run_for_memory_profile.py

@@ -34,19 +34,17 @@ def generate_fits_files(n_images, size=None, seed=1523):
     else:
         use_size = (size, size)
 
-    np.random.seed(seed)
+    rng = np.random.default_rng(seed=seed)
 
     base_name = "test-combine-{num:03d}.fits"
 
     for num in range(n_images):
-        data = np.random.normal(size=use_size)
+        data = rng.normal(size=use_size)
         # Now add some outlying pixels so there is something to clip
         n_bad = 50000
-        bad_x = np.random.randint(0, high=use_size[0] - 1, size=n_bad)
-        bad_y = np.random.randint(0, high=use_size[1] - 1, size=n_bad)
-        data[bad_x, bad_y] = np.random.choice([-1, 1], size=n_bad) * (
-            10 + np.random.rand(n_bad)
-        )
+        bad_x = rng.integers(0, high=use_size[0] - 1, size=n_bad)
+        bad_y = rng.integers(0, high=use_size[1] - 1, size=n_bad)
+        data[bad_x, bad_y] = rng.choice([-1, 1], size=n_bad) * (10 + rng.random(n_bad))
         hdu = fits.PrimaryHDU(data=np.asarray(data, dtype="float32"))
         hdu.header["for_prof"] = "yes"
         hdu.header["bunit"] = "adu"

ccdproc/tests/run_with_file_number_limit.py

@@ -42,7 +42,7 @@ def generate_fits_files(number, size=None):
 def generate_plain_files(number):
     for i in range(number):
         file = TMPPATH / (f"{i:03d}." + ALLOWED_EXTENSIONS["plain"])
-        file.write_bytes(np.random.random(100))
+        file.write_bytes(np.random.default_rng().random(100))
 
 
 def open_files_with_open(kind):

ccdproc/tests/test_bitfield.py

@@ -8,71 +8,71 @@
 
 def test_bitfield_not_integer():
     with pytest.raises(TypeError):
-        bitfield_to_boolean_mask(np.random.random((10, 10)))
+        bitfield_to_boolean_mask(np.random.default_rng().random((10, 10)))
 
 
 def test_bitfield_negative_flags():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, [-1])
 
 
 def test_bitfield_non_poweroftwo_flags():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, [3])
 
 
 def test_bitfield_flipbits_when_no_bits():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(TypeError):
         bitfield_to_boolean_mask(bm, None, flip_bits=1)
 
 
 def test_bitfield_flipbits_when_stringbits():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(TypeError):
         bitfield_to_boolean_mask(bm, "3", flip_bits=1)
 
 
 def test_bitfield_string_flag_flip_not_start_of_string():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, "1, ~4")
 
 
 def test_bitfield_string_flag_unbalanced_parens():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, "(1, 4))")
 
 
 def test_bitfield_string_flag_wrong_positioned_parens():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, "((1, )4)")
 
 
 def test_bitfield_string_flag_empty():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(ValueError):
         bitfield_to_boolean_mask(bm, "~")
 
 
 def test_bitfield_flag_non_integer():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.raises(TypeError):
         bitfield_to_boolean_mask(bm, [1.3])
 
 
 def test_bitfield_duplicate_flag_throws_warning():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     with pytest.warns(UserWarning):
         bitfield_to_boolean_mask(bm, [1, 1])
 
 
 def test_bitfield_none_identical_to_strNone():
-    bm = np.random.randint(0, 10, (10, 10))
+    bm = np.random.default_rng().integers(0, 10, (10, 10))
     m1 = bitfield_to_boolean_mask(bm, None)
     m2 = bitfield_to_boolean_mask(bm, "None")
     np.testing.assert_array_equal(m1, m2)

ccdproc/tests/test_ccdproc.py

@@ -321,7 +321,7 @@ def test_trim_image_fits_section(mask_data, uncertainty):
     if mask_data:
         ccd_data.mask = np.zeros_like(ccd_data)
     if uncertainty:
-        err = np.random.normal(size=ccd_data.shape)
+        err = np.random.default_rng().normal(size=ccd_data.shape)
         ccd_data.uncertainty = StdDevUncertainty(err)
 
     trimmed = trim_image(ccd_data, fits_section="[20:40,:]")
@@ -531,7 +531,7 @@ def test_flat_correct():
     ccd_data.header["my_key"] = 42
     size = ccd_data.shape[0]
     # create the flat, with some scatter
-    data = 2 * np.random.normal(loc=1.0, scale=0.05, size=(size, size))
+    data = 2 * np.random.default_rng().normal(loc=1.0, scale=0.05, size=(size, size))
     flat = CCDData(data, meta=fits.header.Header(), unit=ccd_data.unit)
     flat_data = flat_correct(ccd_data, flat, add_keyword=None)
 
@@ -555,7 +555,7 @@ def test_flat_correct_min_value():
     size = ccd_data.shape[0]
 
     # Create the flat
-    data = 2 * np.random.normal(loc=1.0, scale=0.05, size=(size, size))
+    data = 2 * np.random.default_rng().normal(loc=1.0, scale=0.05, size=(size, size))
     flat = CCDData(data, meta=fits.header.Header(), unit=ccd_data.unit)
     flat_orig_data = flat.data.copy()
     min_value = 2.1  # Should replace some, but not all, values
@@ -589,7 +589,7 @@ def test_flat_correct_norm_value():
     # Note that mean value of flat is set below and is different than
     # the mean of the flat data.
     flat_mean = 5.0
-    data = np.random.normal(loc=1.0, scale=0.05, size=ccd_data.shape)
+    data = np.random.default_rng().normal(loc=1.0, scale=0.05, size=ccd_data.shape)
     flat = CCDData(data, meta=fits.Header(), unit=ccd_data.unit)
     flat_data = flat_correct(ccd_data, flat, add_keyword=None, norm_value=flat_mean)
 
@@ -608,7 +608,7 @@ def test_flat_correct_norm_value_bad_value():
     # it is given a bad norm_value. Bad means <=0.
 
     # Create the flat, with some scatter
-    data = np.random.normal(loc=1.0, scale=0.05, size=ccd_data.shape)
+    data = np.random.default_rng().normal(loc=1.0, scale=0.05, size=ccd_data.shape)
     flat = CCDData(data, meta=fits.Header(), unit=ccd_data.unit)
     with pytest.raises(ValueError) as e:
         flat_correct(ccd_data, flat, add_keyword=None, norm_value=-7)
@@ -694,7 +694,7 @@ def test_transform_image(mask_data, uncertainty):
         ccd_data.mask = np.zeros_like(ccd_data)
         ccd_data.mask[10, 10] = 1
     if uncertainty:
-        err = np.random.normal(size=ccd_data.shape)
+        err = np.random.default_rng().normal(size=ccd_data.shape)
         ccd_data.uncertainty = StdDevUncertainty(err)
 
     def tran(arr):
@@ -948,7 +948,7 @@ def test_wcs_project_onto_shifted_wcs():
     target_wcs = wcs_for_testing(ccd_data.shape)
     target_wcs.wcs.crpix += [1, 1]
 
-    ccd_data.mask = np.random.choice([0, 1], size=ccd_data.shape)
+    ccd_data.mask = np.random.default_rng().choice([0, 1], size=ccd_data.shape)
 
     new_ccd = wcs_project(ccd_data, target_wcs)
 

ccdproc/tests/test_cosmicray.py

@@ -3,7 +3,6 @@
 import numpy as np
 import pytest
 from astropy import units as u
-from astropy.utils import NumpyRNGContext
 from astropy.utils.exceptions import AstropyDeprecationWarning
 
 from ccdproc.core import (
@@ -22,14 +21,17 @@
 
 def add_cosmicrays(data, scale, threshold, ncrays=NCRAYS):
     size = data.shape[0]
-    with NumpyRNGContext(125):
-        crrays = np.random.randint(0, size, size=(ncrays, 2))
-        # use (threshold + 15) below to make sure cosmic ray is well above the
-        # threshold no matter what the random number generator returns
-        crflux = 10 * scale * np.random.random(NCRAYS) + (threshold + 15) * scale
-        for i in range(ncrays):
-            y, x = crrays[i]
-            data.data[y, x] = crflux[i]
+    rng = np.random.default_rng(99)
+    crrays = rng.integers(0, size, size=(ncrays, 2))
+    # use (threshold + 15) below to make sure cosmic ray is well above the
+    # threshold no matter what the random number generator returns
+    # add_cosmicrays is highly sensitive to the seed
+    # ideally threshold should be set so it is not sensitive to seed, but
+    # this is not working right now
+    crflux = 10 * scale * rng.random(NCRAYS) + (threshold + 15) * scale
+    for i in range(ncrays):
+        y, x = crrays[i]
+        data.data[y, x] = crflux[i]
 
 
 def test_cosmicray_lacosmic():
@@ -266,33 +268,29 @@ def test_cosmicray_median_background_deviation():
 
 
 def test_background_deviation_box():
-    with NumpyRNGContext(123):
-        scale = 5.3
-        cd = np.random.normal(loc=0, size=(100, 100), scale=scale)
+    scale = 5.3
+    cd = np.random.default_rng(seed=123).normal(loc=0, size=(100, 100), scale=scale)
     bd = background_deviation_box(cd, 25)
     assert abs(bd.mean() - scale) < 0.10
 
 
 def test_background_deviation_box_fail():
-    with NumpyRNGContext(123):
-        scale = 5.3
-        cd = np.random.normal(loc=0, size=(100, 100), scale=scale)
+    scale = 5.3
+    cd = np.random.default_rng(seed=123).normal(loc=0, size=(100, 100), scale=scale)
     with pytest.raises(ValueError):
         background_deviation_box(cd, 0.5)
 
 
 def test_background_deviation_filter():
-    with NumpyRNGContext(123):
-        scale = 5.3
-        cd = np.random.normal(loc=0, size=(100, 100), scale=scale)
+    scale = 5.3
+    cd = np.random.default_rng(seed=123).normal(loc=0, size=(100, 100), scale=scale)
     bd = background_deviation_filter(cd, 25)
     assert abs(bd.mean() - scale) < 0.10
 
 
 def test_background_deviation_filter_fail():
-    with NumpyRNGContext(123):
-        scale = 5.3
-        cd = np.random.normal(loc=0, size=(100, 100), scale=scale)
+    scale = 5.3
+    cd = np.random.default_rng(seed=123).normal(loc=0, size=(100, 100), scale=scale)
     with pytest.raises(ValueError):
         background_deviation_filter(cd, 0.5)
 

ccdproc/tests/test_rebin.py

@@ -59,7 +59,7 @@ def test_rebin_ccddata(mask_data, uncertainty):
     if mask_data:
         ccd_data.mask = np.zeros_like(ccd_data)
     if uncertainty:
-        err = np.random.normal(size=ccd_data.shape)
+        err = np.random.default_rng().normal(size=ccd_data.shape)
         ccd_data.uncertainty = StdDevUncertainty(err)
 
     with pytest.warns(AstropyDeprecationWarning):

ccdproc/tests/test_wrapped_external_funcs.py

@@ -69,7 +69,7 @@ def test_medianfilter_unusued():
 
 
 def test_medianfilter_ndarray():
-    arr = np.random.random((5, 5))
+    arr = np.random.default_rng().random((5, 5))
     result = core.median_filter(arr, 3)
     reference = ndimage.median_filter(arr, 3)
     # It's a wrapped function so we can use the equal comparison.

docs/ccddata.rst

@@ -139,7 +139,7 @@ Uncertainty
 
 Pixel-by-pixel uncertainty can be calculated for you:
 
-    >>> data = np.random.normal(size=(10, 10), loc=1.0, scale=0.1)
+    >>> data = np.random.default_rng().normal(size=(10, 10), loc=1.0, scale=0.1)
     >>> ccd = CCDData(data, unit="electron")
     >>> ccd_new = ccdproc.create_deviation(ccd, readnoise=5 * u.electron)
 

docs/getting_started.rst

@@ -26,7 +26,7 @@ See the documentation for `~astropy.nddata.CCDData` for a complete list of attri
 
 Most operations are performed by functions in `ccdproc`:
 
-    >>> dark = CCDData(np.random.normal(size=(10, 10)), unit="adu")
+    >>> dark = CCDData(np.random.default_rng().normal(size=(10, 10)), unit="adu")
     >>> dark_sub = ccdproc.subtract_dark(image_1, dark,
     ...                                  dark_exposure=30*u.second,
     ...                                  data_exposure=15*u.second,
@@ -63,7 +63,7 @@ with `~ccdproc.gain_correct` that you do not get with multiplication:
 The same advantages apply to operations that are more complex, like flat
 correction, in which one image is divided by another:
 
-    >>> flat = CCDData(np.random.normal(1.0, scale=0.1, size=(10, 10)),
+    >>> flat = CCDData(np.random.default_rng().normal(1.0, scale=0.1, size=(10, 10)),
     ...                        unit='adu')
     >>> image_1_flat = ccdproc.flat_correct(image_1, flat)