From 0ea174a167911c2137d96fd49ad7dca2afe118f4 Mon Sep 17 00:00:00 2001 From: Tyler Pauly Date: Thu, 29 Feb 2024 13:08:12 -0500 Subject: [PATCH 1/2] JP-3463: Utilize average_dark_current in variance calculations (#243) * utilize dark current in variance calcs * further docs and changelog * fix test calls to ramp data array settting * add array to docstring * bump minor version over patch in changelog * Update CHANGES.rst Get the stupid rst parser to not complain about heading level inconsistencies! * propose fix for roman side * fix slicing for starmap --------- Co-authored-by: Howard Bushouse --- CHANGES.rst | 14 +++++++++++++- docs/stcal/ramp_fitting/description.rst | 6 ++++-- src/stcal/ramp_fitting/ols_fit.py | 11 +++++++---- src/stcal/ramp_fitting/ramp_fit.py | 10 ++++++++-- src/stcal/ramp_fitting/ramp_fit_class.py | 8 +++++++- tests/test_ramp_fitting.py | 20 ++++++++++++++------ tests/test_ramp_fitting_cases.py | 3 ++- tests/test_ramp_fitting_gls_fit.py | 4 +++- 8 files changed, 58 insertions(+), 18 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 664e60d0..90f96964 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,7 +1,19 @@ -1.6.1 (unreleased) +1.7.0 (unreleased) ================== +Changes to API +-------------- + +ramp_fitting +~~~~~~~~~~~~ + +- Add ``average_dark_current`` to calculations of poisson variance. [#243] +Bug Fixes +--------- + +Other +----- 1.6.0 (2024-02-15) ================== diff --git a/docs/stcal/ramp_fitting/description.rst b/docs/stcal/ramp_fitting/description.rst index 37fc625c..9ea6830d 100644 --- a/docs/stcal/ramp_fitting/description.rst +++ b/docs/stcal/ramp_fitting/description.rst @@ -195,14 +195,16 @@ time in seconds (from the keyword TGROUP). The variance of the slope in a segment due to Poisson noise is: .. math:: - var^P_{s} = \frac{ slope_{est} }{ tgroup \times gain\ (ngroups_{s} -1)} \,, + var^P_{s} = \frac{ slope_{est} + darkcurrent}{ tgroup \times gain\ (ngroups_{s} -1)} \,, where :math:`gain` is the gain for the pixel (from the GAIN reference file), in e/DN. The :math:`slope_{est}` is an overall estimated slope of the pixel, calculated by taking the median of the first differences of the groups that are unaffected by saturation and cosmic rays, in all integrations. This is a more robust estimate of the slope than the segment-specific slope, which may be noisy -for short segments. +for short segments. The contributions from the dark current are added when present; +the value can be provided by the user during the `jwst.dark_current.DarkCurrentStep`, +or it can be specified in scalar or 2D array form by the dark reference file. The combined variance of the slope of a segment is the sum of the variances: diff --git a/src/stcal/ramp_fitting/ols_fit.py b/src/stcal/ramp_fitting/ols_fit.py index f2c6d5a2..a52e2fba 100644 --- a/src/stcal/ramp_fitting/ols_fit.py +++ b/src/stcal/ramp_fitting/ols_fit.py @@ -542,8 +542,9 @@ def slice_ramp_data(ramp_data, start_row, nrows): err = ramp_data.err[:, :, start_row : start_row + nrows, :].copy() groupdq = ramp_data.groupdq[:, :, start_row : start_row + nrows, :].copy() pixeldq = ramp_data.pixeldq[start_row : start_row + nrows, :].copy() + average_dark_current = ramp_data.average_dark_current[start_row : start_row + nrows, :].copy() - ramp_data_slice.set_arrays(data, err, groupdq, pixeldq) + ramp_data_slice.set_arrays(data, err, groupdq, pixeldq, average_dark_current) if ramp_data.zeroframe is not None: ramp_data_slice.zeroframe = ramp_data.zeroframe[:, start_row : start_row + nrows, :].copy() @@ -1138,7 +1139,8 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans) # Suppress harmless arithmetic warnings for now warnings.filterwarnings("ignore", ".*invalid value.*", RuntimeWarning) warnings.filterwarnings("ignore", ".*divide by zero.*", RuntimeWarning) - var_p4[num_int, :, rlo:rhi, :] = den_p3 * med_rates[rlo:rhi, :] + var_p4[num_int, :, rlo:rhi, :] = (den_p3 * med_rates[rlo:rhi, :]) + \ + ramp_data.average_dark_current[rlo:rhi, :] # Find the segment variance due to read noise and convert back to DN var_r4[num_int, :, rlo:rhi, :] = num_r3 * den_r3 / gain_sect**2 @@ -1175,10 +1177,11 @@ def ramp_fit_compute_variances(ramp_data, gain_2d, readnoise_2d, fit_slopes_ans) # Huge variances correspond to non-existing segments, so are reset to 0 # to nullify their contribution. var_p3[var_p3 > utils.LARGE_VARIANCE_THRESHOLD] = 0.0 - var_p3[:, med_rates <= 0.0] = 0.0 + med_rate_mask = med_rates <= 0.0 + var_p3[:, med_rate_mask] = 0.0 warnings.resetwarnings() - var_p4[num_int, :, med_rates <= 0.0] = 0.0 + var_p4[num_int, :, med_rate_mask] = ramp_data.average_dark_current[med_rate_mask][..., np.newaxis] var_both4[num_int, :, :, :] = var_r4[num_int, :, :, :] + var_p4[num_int, :, :, :] inv_var_both4[num_int, :, :, :] = 1.0 / var_both4[num_int, :, :, :] diff --git a/src/stcal/ramp_fitting/ramp_fit.py b/src/stcal/ramp_fitting/ramp_fit.py index 514429f1..11ba5b9e 100755 --- a/src/stcal/ramp_fitting/ramp_fit.py +++ b/src/stcal/ramp_fitting/ramp_fit.py @@ -53,10 +53,16 @@ def create_ramp_fit_class(model, dqflags=None, suppress_one_group=False): """ ramp_data = ramp_fit_class.RampData() + if not hasattr(model, 'average_dark_current'): + dark_current_array = np.zeros_like(model.pixeldq) + else: + dark_current_array = model.average_dark_current + if isinstance(model.data, u.Quantity): - ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, model.pixeldq) + ramp_data.set_arrays(model.data.value, model.err.value, model.groupdq, + model.pixeldq, dark_current_array) else: - ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq) + ramp_data.set_arrays(model.data, model.err, model.groupdq, model.pixeldq, dark_current_array) # Attribute may not be supported by all pipelines. Default is NoneType. drop_frames1 = model.meta.exposure.drop_frames1 if hasattr(model, "drop_frames1") else None diff --git a/src/stcal/ramp_fitting/ramp_fit_class.py b/src/stcal/ramp_fitting/ramp_fit_class.py index f8a78efd..ef2df6bc 100644 --- a/src/stcal/ramp_fitting/ramp_fit_class.py +++ b/src/stcal/ramp_fitting/ramp_fit_class.py @@ -6,6 +6,7 @@ def __init__(self): self.err = None self.groupdq = None self.pixeldq = None + self.average_dark_current = None # Meta information self.instrument_name = None @@ -41,7 +42,7 @@ def __init__(self): self.current_integ = -1 - def set_arrays(self, data, err, groupdq, pixeldq): + def set_arrays(self, data, err, groupdq, pixeldq, average_dark_current): """ Set the arrays needed for ramp fitting. @@ -62,12 +63,17 @@ def set_arrays(self, data, err, groupdq, pixeldq): pixeldq : ndarray (uint32) 4-D array containing the pixel data quality information. It has dimensions (nintegrations, ngroups, nrows, ncols) + + average_dark_current : ndarray (float32) + 2-D array containing the average dark current. It has + dimensions (nrows, ncols) """ # Get arrays from the data model self.data = data self.err = err self.groupdq = groupdq self.pixeldq = pixeldq + self.average_dark_current = average_dark_current def set_meta(self, name, frame_time, group_time, groupgap, nframes, drop_frames1=None): """ diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index d8e90610..aab91a12 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -322,12 +322,13 @@ def jp_2326_test_setup(): gdq = np.zeros((nints, ngroups, nrows, ncols), dtype=np.uint8) err = np.zeros((nints, ngroups, nrows, ncols)) pdq = np.zeros((nrows, ncols), dtype=np.uint32) + dark_current = np.zeros((nrows, ncols)) data[0, :, 0, 0] = ramp.copy() gdq[0, :, 0, 0] = dq.copy() ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pdq, average_dark_current=dark_current) ramp_data.set_meta( name="MIRI", frame_time=2.77504, group_time=2.77504, groupgap=0, nframes=1, drop_frames1=None ) @@ -417,6 +418,7 @@ def test_2_group_cases(): rnoise = np.ones((1, npix)) * rnoise_val gain = np.ones((1, npix)) * gain_val pixeldq = np.zeros((1, npix), dtype=np.uint32) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) data = np.zeros(dims, dtype=np.float32) # Science data for k in range(npix): @@ -433,7 +435,7 @@ def test_2_group_cases(): # Setup the RampData class to run ramp fitting on. ramp_data = RampData() - ramp_data.set_arrays(data, err, groupdq, pixeldq) + ramp_data.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRSPEC", frame_time=14.58889, group_time=14.58889, groupgap=0, nframes=1, drop_frames1=None @@ -730,6 +732,8 @@ def create_zero_frame_data(): pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) zframe = np.ones(shape=(nints, nrows, ncols), dtype=np.float32) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) + # Create base ramps for each pixel in each integration. base_slope = 2000.0 @@ -756,7 +760,7 @@ def create_zero_frame_data(): # Create RampData for testing. ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRCam", frame_time=frame_time, @@ -855,6 +859,8 @@ def create_only_good_0th_group_data(): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros((nrows, ncols), dtype=np.float32) + # Create base ramps for each pixel in each integration. base_slope = 2000.0 @@ -877,7 +883,7 @@ def create_only_good_0th_group_data(): # Create RampData for testing. ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRCam", frame_time=frame_time, @@ -1422,9 +1428,10 @@ def create_blank_ramp_data(dims, var, tm): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype = np.float32) ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="NIRSpec", frame_time=frame_time, @@ -1476,6 +1483,7 @@ def setup_inputs(dims, var, tm): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) base_array = np.array([k + 1 for k in range(ngroups)]) base, inc = 1.5, 1.5 @@ -1488,7 +1496,7 @@ def setup_inputs(dims, var, tm): data[c_int, :, :, :] = data[0, :, :, :].copy() ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0, nframes=nframes, drop_frames1=None ) diff --git a/tests/test_ramp_fitting_cases.py b/tests/test_ramp_fitting_cases.py index 675e6a75..9e32813b 100644 --- a/tests/test_ramp_fitting_cases.py +++ b/tests/test_ramp_fitting_cases.py @@ -786,9 +786,10 @@ def create_blank_ramp_data(dims, var, timing, ts_name="NIRSpec"): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) ramp_data = RampData() - ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq) + ramp_data.set_arrays(data=data, err=err, groupdq=gdq, pixeldq=pixdq, average_dark_current=dark_current) ramp_data.set_meta( name=ts_name, frame_time=frame_time, diff --git a/tests/test_ramp_fitting_gls_fit.py b/tests/test_ramp_fitting_gls_fit.py index e940bd94..fef7c0a1 100644 --- a/tests/test_ramp_fitting_gls_fit.py +++ b/tests/test_ramp_fitting_gls_fit.py @@ -63,9 +63,11 @@ def setup_inputs(dims, gain, rnoise, group_time, frame_time): err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) groupdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) pixeldq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) + dark_current = np.zeros(shape=(nrows, ncols), dtype=np.float32) + # Set clas arrays - ramp_class.set_arrays(data, err, groupdq, pixeldq) + ramp_class.set_arrays(data, err, groupdq, pixeldq, average_dark_current=dark_current) # Set class meta ramp_class.set_meta( From 4ddddf8ed0e035b1524c4255abea286971c730eb Mon Sep 17 00:00:00 2001 From: Zach Burnett Date: Fri, 1 Mar 2024 06:41:40 -0500 Subject: [PATCH 2/2] change log for release 1.6.1 (#244) * freeze change log entry * patch release --- CHANGES.rst | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 90f96964..cf7bc1f8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,20 +1,26 @@ -1.7.0 (unreleased) +1.6.2 (unreleased) ================== Changes to API -------------- -ramp_fitting -~~~~~~~~~~~~ - -- Add ``average_dark_current`` to calculations of poisson variance. [#243] - Bug Fixes --------- Other ----- +1.6.1 (2024-02-29) +================== + +Changes to API +-------------- + +ramp_fitting +~~~~~~~~~~~~ + +- Add ``average_dark_current`` to calculations of poisson variance. [#243] + 1.6.0 (2024-02-15) ==================