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 <[email protected]>

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)
 ==================

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:
 

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, :, :, :]
 

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

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):
         """

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
     )

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,

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(