JP-3967: SATURATION flag for partially saturated pixels should be propagated to the _rate DQ

C: Remove outdated comment

Possible partially saturated ramp flagging.

JP-3967 update for pllim (#1)

* Undoing test changes in prep to change the C-extension, as well as analyze the test failures in the LIKELY CI tests.

* Commenting out a part of a test that fails. The OLS SAT flagging needs to be updated before this test works.

* Updated C-extension and tests.

* Updating tests for ramp fitting cases.

* Removing preprocessing switches and adding a comment, cleaning up the code.

* Removing unnecessary comment.

* Comments on how to update slope fitter for partial saturation flagging.

* Updating the parital flagging for the rate product. Updating ramp fit tests for updating flagging.

* Updating partial saturation flagging and tests.

DOC: Reword second paragraph.

Co-authored-by: Ken MacDonald <[email protected]
m>

TST: Extend test_ramp_fitting.py to LIKELY algo

Author: pllim

changes/421.apichange.rst

@@ -0,0 +1 @@
+Set saturation flag for any saturation. This reverts back to saturation flagging behavior prior to #125

docs/stcal/ramp_fitting/description.rst

@@ -195,11 +195,11 @@ Segment-specific Computations
 +++++++++++++++++++++++++++++
 The variance of the slope of a segment due to read noise is:
 
-.. math::  
+.. math::
    var^R_{s} = \frac{12 \ R^2 }{ (ngroups_{s}^3 - ngroups_{s})(tgroup^2)(gain^2) } \,,
 
-where :math:`R` is the noise in the difference between 2 frames, 
-:math:`ngroups_{s}` is the number of groups in the segment, and :math:`tgroup` is the group 
+where :math:`R` is the noise in the difference between 2 frames,
+:math:`ngroups_{s}` is the number of groups in the segment, and :math:`tgroup` is the group
 time in seconds (from the keyword TGROUP).  The divide by gain converts to
 :math:`DN`.  For the special case where as segment has length 1, the
 :math:`ngroups_{s}` is set to :math:`2`.
@@ -272,7 +272,7 @@ The square-root of the combined variance is stored in the ERR array of the outpu
 The overall slope depends on the slope and the combined variance of the slope of each integration's
 segments, so is a sum over integration values computed from the segements:
 
-.. math::    
+.. math::
     slope_{o} = \frac{ \sum_{i}{ \frac{slope_{i}} {var^C_{i}}}} { \sum_{i}{ \frac{1} {var^C_{i}}}}
 
 
@@ -300,21 +300,20 @@ Data Quality Propagation
 For a given pixel, if all groups in an integration are flagged as DO_NOT_USE or
 SATURATED, then that pixel will be flagged as DO_NOT_USE in the corresponding
 integration in the "rateints" product.  Note this does NOT mean that all groups
-are flagged as SATURATED, nor that all groups are flagged as DO_NOT_USE.  For
+are flagged as DO_NOT_USE.  For
 example, slope calculations that are suppressed due to a ramp containing only
 one good group will be flagged as DO_NOT_USE in the
 first group, but not necessarily any other group, while only groups two and
 beyond are flagged as SATURATED.  Further, only if all integrations in the "rateints"
 product are flagged as DO_NOT_USE, then the pixel will be flagged as DO_NOT_USE
 in the "rate" product.
 
-For a given pixel, if all groups in an integration are flagged as SATURATED,
-then that pixel will be flagged as SATURATED and DO_NOT_USE in the corresponding
-integration in the "rateints" product.  This is different from the above case in
-that this is only for all groups flagged as SATURATED, not for some combination
-of DO_NOT_USE and SATURATED.  Further, only if all integrations in the "rateints"
-product are flagged as SATURATED, then the pixel will be flagged as SATURATED
-and DO_NOT_USE in the "rate" product.
+For a given pixel, if any groups in an integration are flagged as SATURATED,
+then that pixel will be flagged as SATURATED in the corresponding
+integration in the "rateints" product.
+Furthermore, if any integration in the "rateints"
+product is flagged as SATURATED, then the pixel will be flagged as SATURATED
+in the "rate" product.
 
 For a given pixel, if any group in an integration is flagged as JUMP_DET, then
 that pixel will be flagged as JUMP_DET in the corresponding integration in the
@@ -338,12 +337,12 @@ group/resultant directly, the likelihood algorithm is based on differences of
 the groups/resultants :math:`d_i = r_i - r_{i-1}`.  The model used to determine
 the slope/countrate, :math:`a`, is:
 
-.. math::    
+.. math::
     \chi^2 = ({\bf d} - a \cdot {\bf 1})^T C ({\bf d} - a \cdot {\bf 1}) \,,
 
-Differentiating, setting to zero, then solving for :math:`a` results in 
+Differentiating, setting to zero, then solving for :math:`a` results in
 
-.. math::    
+.. math::
     a = ({\bf 1}^T C {\bf d})({\bf 1}^T C {\bf 1})^T \,,
 
 The covariance matrix :math:`C` is a tridiagonal matrix, due to the nature of the

src/stcal/ramp_fitting/src/slope_fitter.c

@@ -14,17 +14,6 @@
 #include <numpy/arrayobject.h>
 #include <numpy/npy_math.h>
 
-/*
-To build C code, make sure the setup.py file is correct and
-lists all extensions, then run:
-
-python setup.py build_ext --inplace
-
-        or
-
-pip install -e .
- */
-
 /* ========================================================================= */
 /*                               TYPEDEFs                                    */
 /* ------------------------------------------------------------------------- */
@@ -2746,6 +2735,7 @@ ramp_fit_pixel(
     int ret = 0;
     npy_intp integ;
     int sat_cnt = 0, dnu_cnt = 0;
+    int set_rate_sat_flag = 0;
 
     /* Ramp fitting depends on the averaged median rate for each integration */
     if (compute_median_rate(rd, pr)) {
@@ -2777,20 +2767,28 @@ ramp_fit_pixel(
             dnu_cnt++;
             pr->rateints[integ].slope = NAN;
         }
+
+        if (rd->save_opt) {
+            get_pixel_ramp_integration_segments_and_pedestal(integ, pr, rd);
+        }
+
         if (pr->rateints[integ].dq & rd->sat) {
             sat_cnt++;
             pr->rateints[integ].slope = NAN;
         }
 
-        if (rd->save_opt) {
-            get_pixel_ramp_integration_segments_and_pedestal(integ, pr, rd);
+        // The partial saturation must go here to not mess up pedestal computations
+        if (pr->stats[integ].cnt_sat > 0) {
+            pr->rateints[integ].dq |= rd->sat;
+            set_rate_sat_flag = 1;
         }
     }
 
     if (rd->nints == dnu_cnt) {
         pr->rate.dq |= rd->dnu;
     }
-    if (rd->nints == sat_cnt) {
+
+    if (sat_cnt == rd->nints) {
         pr->rate.dq |= rd->sat;
     }
 
@@ -2815,6 +2813,11 @@ ramp_fit_pixel(
         pr->rate.var_err = 0.;
     }
 
+    // Partial saturation flagging must be done here
+    if (set_rate_sat_flag) {
+        pr->rate.dq |= rd->sat;
+    }
+
     if (!isnan(pr->rate.slope)) {
         pr->rate.slope = pr->rate.slope / pr->invvar_e_sum;
     }
@@ -2968,6 +2971,7 @@ ramp_fit_pixel_integration(
         goto END;
     }
 
+    // Whole ramp not usable
     if (rd->ngroups == pr->stats[integ].cnt_dnu_sat) {
         pr->rateints[integ].dq |= rd->dnu;
         if (rd->ngroups == pr->stats[integ].cnt_sat) {

src/stcal/ramp_fitting/utils.py

@@ -111,6 +111,11 @@ def dq_compress_sect(ramp_data, num_int, gdq_sect, pixeldq_sect):
     gdq0_sat = np.bitwise_and(gdq_sect[0], sat)
     pixeldq_sect[gdq0_sat != 0] = np.bitwise_or(pixeldq_sect[gdq0_sat != 0], sat | dnu)
 
+    # If saturation occurs mark the appropriate flag.
+    sat_loc = np.bitwise_and(gdq_sect, sat)
+    sat_check = np.where(sat_loc.sum(axis=0) > 0)
+    pixeldq_sect[sat_check] = np.bitwise_or(pixeldq_sect[sat_check], sat)
+
     # If jump occurs mark the appropriate flag.
     jump_loc = np.bitwise_and(gdq_sect, jump)
     jump_check = np.where(jump_loc.sum(axis=0) > 0)
@@ -144,11 +149,11 @@ def dq_compress_final(dq_int, ramp_data):
     dnu = np.uint32(ramp_data.flags_do_not_use)
     sat = np.uint32(ramp_data.flags_saturated)
 
-    # Remove DO_NOT_USE and SATURATED because they need special handling.
-    # These flags are not set in the final pixel DQ array by simply being set
+    # Remove DO_NOT_USE because it needs special handling.
+    # This flag is not set in the final pixel DQ array by simply being set
     # in one of the integrations.
-    not_sat_or_dnu = np.uint32(~(dnu | sat))
-    final_dq = np.bitwise_and(final_dq, not_sat_or_dnu)
+    not_dnu = np.uint32(~dnu)
+    final_dq = np.bitwise_and(final_dq, not_dnu)
 
     # If all integrations are DO_NOT_USE or SATURATED, then set DO_NOT_USE.
     set_if_total_integ(final_dq, dq_int, dnu | sat, dnu)