Add frame zero effects. (#130)

* Add frame zero effects.
* Updates for numpy 2 and ruff.

romanisim/image.py

@@ -103,6 +103,10 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None,
 
     if gain is None:
         gain = parameters.reference_data['gain']
+    try:
+        gain = gain.astype('f4')
+    except AttributeError:  # gain is not a Quantity
+        gain = np.float32(gain)
 
     if linearity is not None:
         resultants = linearity.apply(resultants)
@@ -145,7 +149,7 @@ def make_l2(resultants, read_pattern, read_noise=None, gain=None, flat=None,
     poissonvar = rampvar[..., 1, 1, 1] / gain**2 * (u.DN / u.s)**2
 
     if flat is not None:
-        flat = np.clip(flat, 1e-9, np.inf)
+        flat = np.clip(flat, 1e-9, np.inf).astype('f4')
         slopes /= flat
         readvar /= flat**2
         poissonvar /= flat**2
@@ -756,6 +760,7 @@ def simulate(metadata, objlist,
     saturation = refdata['saturation']
     reffiles = refdata['reffiles']
     flat = refdata['flat']
+    pedestal_extra_noise = parameters.pedestal_extra_noise
 
     if rng is None and seed is None:
         seed = 43
@@ -776,7 +781,9 @@ def simulate(metadata, objlist,
         im = dict(data=counts.array, meta=dict(image_mod.meta.items()))
     else:
         l1, l1dq = romanisim.l1.make_l1(
-            counts, ma_table_number, read_noise=read_noise, rng=rng, gain=gain,
+            counts, ma_table_number, read_noise=read_noise,
+            pedestal_extra_noise=pedestal_extra_noise,
+            rng=rng, gain=gain,
             crparam=crparam,
             inv_linearity=inv_linearity,
             tstart=image_mod.meta.exposure.start_time,

romanisim/l1.py

@@ -331,7 +331,7 @@ def apportion_counts_to_resultants(
 
 
 def add_read_noise_to_resultants(resultants, tij, read_noise=None, rng=None,
-                                 seed=None):
+                                 seed=None, pedestal_extra_noise=None):
     """Adds read noise to resultants.
 
     The resultants get Gaussian read noise with sigma = sigma_read/sqrt(N).
@@ -351,6 +351,9 @@ def add_read_noise_to_resultants(resultants, tij, read_noise=None, rng=None,
         Random number generator to use
     seed : int
         Seed for populating RNG.  Only used if rng is None.
+    pedestal_extra_noise : float
+        Extra read noise to add to each pixel across all groups.
+        Equivalent to noise in the reference read.
 
     Returns
     -------
@@ -366,6 +369,9 @@ def add_read_noise_to_resultants(resultants, tij, read_noise=None, rng=None,
     else:
         rng = galsim.GaussianDeviate(rng)
 
+    # separate noise generator for pedestals so we can turn it on and off.
+    pedestalrng = galsim.GaussianDeviate(rng.raw())
+
     if read_noise is None:
         read_noise = parameters.reference_data['readnoise']
     if read_noise is None:
@@ -377,6 +383,13 @@ def add_read_noise_to_resultants(resultants, tij, read_noise=None, rng=None,
     noise = noise * read_noise / np.array(
         [len(x)**0.5 for x in tij]).reshape(-1, 1, 1)
     resultants += noise
+
+    if pedestal_extra_noise is not None:
+        noise = np.zeros(resultants.shape[1:], dtype='f4')
+        amplitude = np.hypot(pedestal_extra_noise, read_noise)
+        pedestalrng.generate(noise)
+        resultants += noise[None, ...] * amplitude
+
     return resultants
 
 
@@ -476,7 +489,8 @@ def add_ipc(resultants, ipc_kernel=None):
 
 
 def make_l1(counts, read_pattern,
-            read_noise=None, rng=None, seed=None,
+            read_noise=None, pedestal_extra_noise=None,
+            rng=None, seed=None,
             gain=None, inv_linearity=None, crparam=None,
             persistence=None, tstart=None, saturation=None):
     """Make an L1 image from a counts image.
@@ -493,6 +507,8 @@ def make_l1(counts, read_pattern,
         resultant.
     read_noise : np.ndarray[nx, ny] (float) or float
         Read noise entering into each read
+    pedestal_extra_noise : np.ndarray[nx, ny] (float) or float
+        Extra noise entering into each pixel (i.e., degenerate with pedestal)
     rng : galsim.BaseDeviate
         Random number generator to use
     seed : int
@@ -549,7 +565,8 @@ def make_l1(counts, read_pattern,
     log.info('Adding read noise...')
     resultants = add_read_noise_to_resultants(
         resultants, tij, rng=rng, seed=seed,
-        read_noise=read_noise)
+        read_noise=read_noise,
+        pedestal_extra_noise=pedestal_extra_noise)
 
     # quantize
     resultants = np.round(resultants)

romanisim/parameters.py

@@ -113,6 +113,9 @@
 # arbitrary constant to add to initial L1 image so that pixels aren't clipped at zero.
 pedestal = 100 * u.DN
 
+# Addd this much extra noise as correlated extra noise in all resultants.
+pedestal_extra_noise = 4 * u.DN
+
 dqbits = dict(saturated=2, jump_det=4, nonlinear=2**16, no_lin_corr=2**20)
 dq_do_not_use = dqbits['saturated'] | dqbits['jump_det']
 

romanisim/tests/test_l1.py

@@ -70,9 +70,11 @@ def test_apportion_counts_to_resultants():
     counts_no_poisson_noise = 100
     counts = np.random.poisson(counts_no_poisson_noise, size=(100, 100))
     read_noise = 10
+    pedestal_extra_noise = 20
     res1out = []
     res2out = []
     res3out = []
+    res4out = []
     for tij in tijlist:
         resultants, dq = l1.apportion_counts_to_resultants(counts, tij)
         assert np.all(np.diff(resultants, axis=0) >= 0)
@@ -82,11 +84,15 @@ def test_apportion_counts_to_resultants():
                                                tij)
         res3 = l1.add_read_noise_to_resultants(resultants.copy(), tij,
                                                read_noise=read_noise)
+        res4 = l1.add_read_noise_to_resultants(
+            resultants.copy(), tij, read_noise=0,
+            pedestal_extra_noise=pedestal_extra_noise)
         assert np.all(res2 != resultants)
         assert np.all(res3 != resultants)
         res1out.append(resultants)
         res2out.append(res2)
         res3out.append(res3)
+        res4out.append(res4)
         for restij, plane_index in zip(tij, np.arange(res3.shape[0])):
             predcounts = (np.mean(restij) * counts_no_poisson_noise
                           / tij[-1][-1])
@@ -96,8 +102,16 @@ def test_apportion_counts_to_resultants():
             sdev = np.std(res3[plane_index] - resultants[plane_index])
             assert (np.abs(sdev - read_noise / np.sqrt(len(restij)))
                     < 20 * sdev / np.sqrt(2 * len(counts.ravel())))
+            sdev = np.std(res4[plane_index] - resultants[plane_index])
+            assert (np.abs(sdev - pedestal_extra_noise)
+                    < 20 * sdev / np.sqrt(2 * len(counts.ravel())))
+        # pedestal extra read noise should be correlated and cancel out of the
+        # first difference, and so should agree exactly with resultants
+        assert np.allclose(np.diff(res4 - resultants, axis=0), 0, atol=1e-4)
     log.info('DMS220: successfully added read noise to resultants.')
     log.info('DMS229: successfully generated ramp from counts.')
+    log.info('DMS223: successfully added correlated noise associated '
+             'with frame zero')
 
     artifactdir = os.environ.get('TEST_ARTIFACT_DIR', None)
     if artifactdir is not None:
@@ -108,6 +122,14 @@ def test_apportion_counts_to_resultants():
                    'tijlist': tijlist}
         af.write_to(os.path.join(artifactdir, 'dms229.asdf'))
 
+        af = asdf.AsdfFile()
+        af.tree = {'resultants': res4out,
+                   'resultants_nonoise': res1out,
+                   'counts': counts,
+                   'counts_no_poisson_noise': 100,
+                   'tijlist': tijlist}
+        af.write_to(os.path.join(artifactdir, 'dms223.asdf'))
+
 
 @metrics_logger("DMS222")
 @pytest.mark.soctests

romanisim/tests/test_ramp.py

@@ -87,8 +87,8 @@ def test_ramp(test_table=None):
         scale = flux_on_readvar_pts[i]
         ki, var = ramp.construct_ki_and_variances(
             atcinva, atcinv, [rcovar, fcovar * scale, acovar])
-        assert np.allclose(kigrid[i], ki, atol=1e-6)
-        assert np.allclose(vargrid[i], var, atol=1e-6)
+        assert np.allclose(kigrid[i], ki, atol=1e-5, rtol=1e-4)
+        assert np.allclose(vargrid[i], var, atol=1e-5, rtol=1e-4)
     fitter = ramp.RampFitInterpolator(test_table, flux_on_readvar_pts)
     ki = fitter.ki(flux, read_noise)
     var = fitter.variances(flux, read_noise)

tox.ini

@@ -19,7 +19,7 @@ skip_install = true
 deps =
     ruff
 commands =
-    ruff . {posargs}
+    ruff check . {posargs}
 
 [testenv:check-security]
 description = run bandit to check security compliance