diff --git a/romanisim/image.py b/romanisim/image.py index e5116ee7..617710bf 100644 --- a/romanisim/image.py +++ b/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, diff --git a/romanisim/l1.py b/romanisim/l1.py index cbdd660c..741cbdc6 100644 --- a/romanisim/l1.py +++ b/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) diff --git a/romanisim/parameters.py b/romanisim/parameters.py index 221a14f2..bbf79227 100644 --- a/romanisim/parameters.py +++ b/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'] diff --git a/romanisim/tests/test_l1.py b/romanisim/tests/test_l1.py index e33e79e3..d4a7aea4 100644 --- a/romanisim/tests/test_l1.py +++ b/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 diff --git a/romanisim/tests/test_ramp.py b/romanisim/tests/test_ramp.py index affbdd3e..d042c767 100644 --- a/romanisim/tests/test_ramp.py +++ b/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) diff --git a/tox.ini b/tox.ini index 24a8c13d..6fa44c55 100644 --- a/tox.ini +++ b/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