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

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}}}}
 
 
@@ -297,24 +297,23 @@ in the "rate" product they contain values for the exposure as a whole.
 
 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
+For a given pixel, if all groups in an integration are flagged as DO_NOT_USE,
+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,20 +337,20 @@ 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
 differences.  Because the covariance matrix is tridiagonal, the  computational
 complexity reduces from :math:`O(n^3)` to :math:`O(n)`.  To see the detailed
 derivation and computations implemented, refer to the links above.
 The Poisson and read noise  computations are based on equations (27) and (28), in
-`Brandt (2024) <https://iopscience.iop.org/article/10.1088/1538-3873/ad38d9>`__
+`Brandt (2024) <https://iopscience.iop.org/article/10.1088/1538-3873/ad38d9>`__.
 
 For more details, especially for the jump detection portion in the liklihood
 algorithm, see

src/stcal/ramp_fitting/ramp_fit.py

@@ -27,6 +27,8 @@
 
 BUFSIZE = 1024 * 300000  # 300Mb cache size for data section
 
+__all__ = ["create_ramp_fit_class", "ramp_fit", "ramp_fit_data", "suppress_one_good_group_ramps"]
+
 
 def create_ramp_fit_class(model, algorithm, dqflags=None, suppress_one_group=False):
     """
@@ -199,7 +201,8 @@ def ramp_fit_data(
     depending on the choice of ramp fitting algorithms (which is only ordinary
     least squares right now) and the choice of single or muliprocessing.
 
-
+    Parameters
+    ----------
     ramp_data : RampData
         Input data necessary for computing ramp fitting.
 
@@ -247,7 +250,7 @@ def ramp_fit_data(
                   " but ngroups = {ngroups}.  Due to this, the ramp fitting algorithm"
                   " is being changed to OLS_C")
         algorithm = "OLS_C"
-        
+
     if algorithm.upper() == "LIKELY" and ngroups >= likely_fit.LIKELY_MIN_NGROUPS:
         image_info, integ_info, opt_info = likely_fit.likely_ramp_fit(
             ramp_data, readnoise_2d, gain_2d

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

@@ -12,6 +12,8 @@
 LARGE_VARIANCE = 1.0e8
 LARGE_VARIANCE_THRESHOLD = 0.01 * LARGE_VARIANCE
 
+__all__ = ["set_if_total_ramp", "set_if_total_integ", "dq_compress_sect", "dq_compress_final"]
+
 
 def set_if_total_ramp(pixeldq_sect, gdq_sect, flag, set_flag):
     """
@@ -72,11 +74,10 @@ def set_if_total_integ(final_dq, integ_dq, flag, set_flag):
 
 def dq_compress_sect(ramp_data, num_int, gdq_sect, pixeldq_sect):
     """
-    This sets the integration level flags for DO_NOT_USE, JUMP_DET and
-    SATURATED.  If any ramp has a jump, this flag will be set for the
+    This sets the integration level flags for DO_NOT_USE, JUMP_DET, and
+    SATURATED.  If any ramp has a jump or saturated, the respective flag will be set for the
     integration.  If all groups in a ramp are flagged as DO_NOT_USE, then the
-    integration level DO_NOT_USE flag will be set.  If a ramp is saturated in
-    group 0, then the integration level flag is marked as SATURATED.  Also, if
+    integration level DO_NOT_USE flag will be set.  Also, if
     all groups are marked as DO_NOT_USE or SATURATED (as in suppressed one
     groups), then the DO_NOT_USE flag is set.
 
@@ -107,9 +108,10 @@ def dq_compress_sect(ramp_data, num_int, gdq_sect, pixeldq_sect):
     # Check total SATURATED or DO_NOT_USE
     set_if_total_ramp(pixeldq_sect, gdq_sect, sat | dnu, dnu)
 
-    # Assume total saturation if group 0 is SATURATED.
-    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)
@@ -142,18 +144,14 @@ def dq_compress_final(dq_int, ramp_data):
         final_dq = np.bitwise_or(final_dq, dq_int[integ, :, :])
 
     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)
-
-    # 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)
+    not_dnu = np.uint32(~dnu)
+    final_dq = np.bitwise_and(final_dq, not_dnu)
 
-    # If all integrations have SATURATED, then set DO_NOT_USE and SATURATED.
-    set_if_total_integ(final_dq, dq_int, sat, dnu | sat)
+    # If all integrations are DO_NOT_USE, then set DO_NOT_USE.
+    set_if_total_integ(final_dq, dq_int, dnu, dnu)
 
     return final_dq