Repository: spacetelescope/jwst

Issue:
IRS2 intermittently bad pixels algorithm fix
_Issue [JP-3321](https://jira.stsci.edu/browse/JP-3321) was created on JIRA by [Maria Pena-Guerrero](https://jira.stsci.edu/secure/ViewProfile.jspa?name=pena):_

This ticket is related to JP-3157. The algorithm is a dynamic check on the reference pixel values. The previous overcorrection mitigation factor is turned into the number of sigmas away from the mean (more intuitive quantity). This is what the algorithm does:
 # Loop through the reference pixel sequences and calculates per column the mean, standard deviation, and differences between even and odd, and store these in their corresponding array.
 # Calculate the means and standard deviations for each of these 3 arrays
 # Flag anything  >  mean_of_array + FACTOR * standard_deviation_of_array
 # Replace these flagged (suspicious) pixels with the nearest good ref pix

Discussion:
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) Report
Patch coverage: **`80.00%`** and no project coverage change.
> Comparison is base [(`d1a60ee`)](https://app.codecov.io/gh/spacetelescope/jwst/commit/d1a60ee11df7bbeb3b7ea7075ce933b5756a3faa?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 76.60% compared to head [(`38f51f6`)](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) 76.60%.

<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##           master    #7745   +/-   ##
=======================================
  Coverage   76.60%   76.60%           
=======================================
  Files         456      456           
  Lines       36950    36950           
=======================================
  Hits        28304    28304           
  Misses       8646     8646           
```

| Flag | Coverage Δ | | *Carryforward flag |
|---|---|---|---|
| nightly | `77.43% <ø> (ø)` | | Carriedforward from [d1a60ee](https://app.codecov.io/gh/spacetelescope/jwst/commit/d1a60ee11df7bbeb3b7ea7075ce933b5756a3faa?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) |

*This pull request uses carry forward flags. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) to find out more.

| [Files Changed](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope) | Coverage Δ | |
|---|---|---|
| [jwst/refpix/refpix\_step.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9yZWZwaXgvcmVmcGl4X3N0ZXAucHk=) | `73.17% <ø> (ø)` | |
| [jwst/refpix/irs2\_subtract\_reference.py](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope#diff-andzdC9yZWZwaXgvaXJzMl9zdWJ0cmFjdF9yZWZlcmVuY2UucHk=) | `86.62% <80.00%> (ø)` | |


</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/spacetelescope/jwst/pull/7745?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=spacetelescope).

penaguerrero: ok, regression tests ran. The test shows very small differences expected due to the changes made.
hbushouse: > ok, regression tests ran. The test shows very small differences expected due to the changes made.

Please post a link to the regression test run, so that the results can be reviewed.
penaguerrero: regression test link: https://plwishmaster.stsci.edu:8081/job/RT/job/JWST-Developers-Pull-Requests/836/

Patch:
diff --git a/CHANGES.rst b/CHANGES.rst
index 4f3312c51..8a5537282 100644
--- a/CHANGES.rst
+++ b/CHANGES.rst
@@ -28,7 +28,7 @@ flat_field
 ----------
 
 - Modify the test_flatfield_step_interface unit test to prevent it from causing
-  other tests to fail [#7752]
+  other tests to fail. [#7752]
 
 general
 -------
@@ -51,6 +51,13 @@ pathloss
 
 - Fix interpolation error for point source corrections. [#7799]
 
+refpix
+------
+
+- Modified algorithm of intermittent bad pixels factor to be the number
+  of sigmas away from mean for the corresponding array (either differences,
+  means, or standard deviations arrays). [#7745]
+
 resample
 --------
 
diff --git a/docs/jwst/refpix/arguments.rst b/docs/jwst/refpix/arguments.rst
index 7f39dd962..6cad2fc8e 100644
--- a/docs/jwst/refpix/arguments.rst
+++ b/docs/jwst/refpix/arguments.rst
@@ -41,5 +41,6 @@ default value is True, and this argument applies to MIR data only.
 *  ``--ovr_corr_mitigation_ftr``
 
 This is a factor to avoid overcorrection of intermittently bad reference
-pixels in the IRS2 algorithm. The default value is 1.8, and this argument
-applies only to NIRSpec data taken with IRS2 mode.
+pixels in the IRS2 algorithm. This factor is the number of sigmas away
+from the mean. The default value is 3.0, and this argument applies
+only to NIRSpec data taken with IRS2 mode.
diff --git a/docs/jwst/refpix/description.rst b/docs/jwst/refpix/description.rst
index cf3df5ce0..635abb32f 100644
--- a/docs/jwst/refpix/description.rst
+++ b/docs/jwst/refpix/description.rst
@@ -198,11 +198,14 @@ will be set to zero if they are flagged as bad in the DQ extension.
 At this point the algorithm looks for intermittently bad (or suspicious)
 reference pixels. This is done by calculating the means and standard
 deviations per reference pixel column, as well as the difference between
-even and odd pairs; then calculates the mean of each of these arrays (the
-mean of the absolute values for the differences array), and flag all
-values greater than the corresponding mean times a factor to avoid
-overcorrection. All suspicious pixels will be replaced by their
-nearest good reference pixel.
+even and odd pairs; then calculates the mean and standard deviation of
+each of these arrays (the mean of the absolute values for the
+differences array), and flag all values greater than the corresponding
+mean plus the standard deviation times a factor to avoid overcorrection.
+All suspicious pixels will be replaced by their nearest good reference
+pixel, or set to zero if there were no good reference pixels left
+(although this is unlikely to happen as there are typically only a few
+pixels flagged as suspicious).
 
 The next step in this processing is to
 copy the science data and the reference pixel data separately to temporary
diff --git a/jwst/refpix/irs2_subtract_reference.py b/jwst/refpix/irs2_subtract_reference.py
index e49b897b1..9bdeba766 100644
--- a/jwst/refpix/irs2_subtract_reference.py
+++ b/jwst/refpix/irs2_subtract_reference.py
@@ -9,7 +9,7 @@ log.setLevel(logging.DEBUG)
 
 def correct_model(input_model, irs2_model,
                   scipix_n_default=16, refpix_r_default=4, pad=8,
-                  ovr_corr_mitigation_ftr=1.8):
+                  ovr_corr_mitigation_ftr=3.0):
     """Correct an input NIRSpec IRS2 datamodel using reference pixels.
 
     Parameters
@@ -34,7 +34,8 @@ def correct_model(input_model, irs2_model,
         the phase of temporally periodic signals.
 
     ovr_corr_mitigation_ftr: float
-        Factor to avoid overcorrection of intermittently bad reference pixels
+        Factor to avoid overcorrection of intermittently bad reference
+        pixels. This factor is the N sigmas away from the mean.
 
     Returns
     -------
@@ -146,7 +147,7 @@ def correct_model(input_model, irs2_model,
 
     # Compute and apply the correction to one integration at a time
     for integ in range(n_int):
-        log.info(f'Working on integration {integ+1}')
+        log.info(f'Working on integration {integ+1} out of {n_int}')
 
         # The input data have a length of 3200 for the last axis (X), while
         # the output data have an X axis with length 2048, the same as the
@@ -343,9 +344,9 @@ def clobber_ref(data, output, odd_even, mask, scipix_n=16, refpix_r=4):
         for k in bits:
             ref = (offset + scipix_n // 2 + k * (scipix_n + refpix_r) +
                    2 * (odd_even_row - 1))
-            log.debug("bad interleaved reference at pixels {} {}"
-                      .format(ref, ref + 1))
-            data[..., ref:ref + 2] = 0.
+            log.debug("bad interleaved reference at pixels {} through {}"
+                      .format(ref, ref + 4))
+            data[..., ref:ref + 4] = 0.
 
 
 def decode_mask(output, mask):
@@ -377,14 +378,12 @@ def decode_mask(output, mask):
 
     # The bit number corresponds to a count of groups of reads of the
     # interleaved reference pixels. The 32-bit unsigned integer encoding
-    # has increasing index, following the amplifier readout direction.
+    # has increasing index, from left to right.
 
     flags = np.array([2**n for n in range(32)], dtype=np.uint32)
     temp = np.bitwise_and(flags, mask)
     bits = np.where(temp > 0)[0]
     bits = list(bits)
-    if output // 2 * 2 == output:
-        bits = [31 - bit for bit in bits]
     bits.sort()
 
     return bits
@@ -410,7 +409,8 @@ def rm_intermittent_badpix(data, scipix_n, refpix_r, ovr_corr_mitigation_ftr):
         regular samples.
 
     ovr_corr_mitigation_ftr: float
-        Factor to avoid overcorrection of bad intermittent reference pixels
+        Factor to avoid overcorrection of intermittently bad reference
+        pixels. This factor is the N sigmas away from the mean.
 
     Returns
     -------
@@ -419,6 +419,9 @@ def rm_intermittent_badpix(data, scipix_n, refpix_r, ovr_corr_mitigation_ftr):
         science and interleaved reference pixel values.  The intermittently
         bad pixels are now set to the nearest good reference pixel value.
     """
+
+    log.info('Using overcorrection mitigation factor = {}'.format(ovr_corr_mitigation_ftr))
+
     # The intermittently bad pixels will be replaced for all integrations
     # and all groups. The last group will be used to identify them
     nints, ngroups, ny, nx = np.shape(data)
@@ -458,23 +461,26 @@ def rm_intermittent_badpix(data, scipix_n, refpix_r, ovr_corr_mitigation_ftr):
                         rp2check.append(odd_pix)
                     pair = 0
         diff_m = np.mean(np.abs(diffs))
+        std_of_diffs = np.std(diffs)
         mean_mean = np.mean(rp_means)
+        std_of_means = np.std(rp_means)
         mean_std = np.mean(rp_stds)
+        std_of_std = np.std(rp_stds)
 
         # order indexes increasing from left to right
         rp2check.sort()
 
         # find the additional intermittent bad pixels - the factor is to avoid overcorrection
-        log.info('Using overcorrection mitigation factor = {}'.format(ovr_corr_mitigation_ftr))
-        high_diffs = np.where(np.abs(diffs) > ovr_corr_mitigation_ftr * diff_m)[0]
+        high_diffs = np.where(np.abs(diffs) > ovr_corr_mitigation_ftr*std_of_diffs + diff_m)[0]
         hd_rp2replace = []
         for j in high_diffs:
             rp2r = rp2check[int(diffs.index(diffs[j]) * 2)]
             # include both even and odd
             hd_rp2replace.append(rp2r)
             hd_rp2replace.append(rp2r+1)
-        high_means_idx = np.where(np.array(rp_means) > ovr_corr_mitigation_ftr * mean_mean)[0]
-        high_std_idx = np.where(np.array(rp_stds) > ovr_corr_mitigation_ftr * mean_std)[0]
+        high_means_idx = np.where(np.array(rp_means) > ovr_corr_mitigation_ftr*std_of_means + mean_mean)[0]
+        high_std_idx = np.where(np.array(rp_stds) > ovr_corr_mitigation_ftr*std_of_std + mean_std)[0]
+        log.info('high_diffs={}  high_means={}  high_stds={}'.format(len(high_diffs), len(high_means_idx), len(high_std_idx)))
         ref_pix = np.array(ref_pix)
         rp2replace = []
         for rp in ref_pix:
@@ -499,12 +505,13 @@ def rm_intermittent_badpix(data, scipix_n, refpix_r, ovr_corr_mitigation_ftr):
                 remaining_rp = remaining_rp_even
             else:
                 remaining_rp = remaining_rp_odd
-            if len(remaining_rp) == 0:   # claims all ref pix are bad, avoid overcorrection and skip
-                continue
-            good_idx = (np.abs(remaining_rp - bad_pix)).argmin()
-            good_pix = remaining_rp[good_idx]
-            data[..., bad_pix] = data[..., good_pix]
-            log.debug('   Pixel {}'.format(bad_pix))
+            if len(remaining_rp) == 0:   # no good ref pix left, let interpolation do it's magic
+                data[..., bad_pix] = 0.
+            else:
+                good_idx = (np.abs(remaining_rp - bad_pix)).argmin()
+                good_pix = remaining_rp[good_idx]
+                data[..., bad_pix] = data[..., good_pix]
+            log.info('   Pixel {}'.format(bad_pix))
     log.info('Total intermittent bad reference pixels: {}'.format(len(total_rp2replace)))
 
 
diff --git a/jwst/refpix/refpix_step.py b/jwst/refpix/refpix_step.py
index 30f709945..a79289af0 100644
--- a/jwst/refpix/refpix_step.py
+++ b/jwst/refpix/refpix_step.py
@@ -22,7 +22,7 @@ class RefPixStep(Step):
         side_smoothing_length = integer(default=11)
         side_gain = float(default=1.0)
         odd_even_rows = boolean(default=True)
-        ovr_corr_mitigation_ftr = float(default=1.8)
+        ovr_corr_mitigation_ftr = float(default=3.0)
     """
 
     reference_file_types = ['refpix']
@@ -53,7 +53,7 @@ class RefPixStep(Step):
 
                 # Apply the IRS2 correction scheme
                 result = irs2_subtract_reference.correct_model(input_model, irs2_model,
-                                                               self.ovr_corr_mitigation_ftr)
+                                  ovr_corr_mitigation_ftr=self.ovr_corr_mitigation_ftr)
 
                 if result.meta.cal_step.refpix != 'SKIPPED':
                     result.meta.cal_step.refpix = 'COMPLETE'


Test patch:
diff --git a/jwst/refpix/tests/test_clobber_ref.py b/jwst/refpix/tests/test_clobber_ref.py
index b368a74db..88dc54bd4 100644
--- a/jwst/refpix/tests/test_clobber_ref.py
+++ b/jwst/refpix/tests/test_clobber_ref.py
@@ -4,7 +4,7 @@ from jwst.refpix.irs2_subtract_reference import decode_mask, clobber_ref
 
 
 def test_clobber_ref():
-    data = np.ones((2, 3, 5, 3200), dtype=np.float32)
+    data = np.ones((2, 3, 5, 3200))
 
     output = np.array([1, 1, 2, 2, 3, 3, 4, 4], dtype=np.int16)
     odd_even = np.array([1, 2, 1, 2, 1, 2, 1, 2], dtype=np.int16)
@@ -15,38 +15,24 @@ def test_clobber_ref():
                      2**5 + 2**7,
                      2**11 + 2**13,
                      0,
-                     2**4],
+                     0],
                     dtype=np.uint32)
 
     clobber_ref(data, output, odd_even, mask)
 
-    compare = np.ones((2, 3, 5, 3200), dtype=np.float32)
-    compare[:, :, :, 648] = 0.
-    compare[:, :, :, 649] = 0.
-    compare[:, :, :, 668] = 0.
-    compare[:, :, :, 669] = 0.
-    compare[:, :, :, 690] = 0.
-    compare[:, :, :, 691] = 0.
-    compare[:, :, :, 710] = 0.
-    compare[:, :, :, 711] = 0.
-    compare[:, :, :, 1290] = 0.
-    compare[:, :, :, 1291] = 0.
-    compare[:, :, :, 1310] = 0.
-    compare[:, :, :, 1311] = 0.
-    compare[:, :, :, 1368] = 0.
-    compare[:, :, :, 1369] = 0.
-    compare[:, :, :, 1388] = 0.
-    compare[:, :, :, 1389] = 0.
-    compare[:, :, :, 2028] = 0.
-    compare[:, :, :, 2029] = 0.
-    compare[:, :, :, 2068] = 0.
-    compare[:, :, :, 2069] = 0.
-    compare[:, :, :, 2150] = 0.
-    compare[:, :, :, 2151] = 0.
-    compare[:, :, :, 2190] = 0.
-    compare[:, :, :, 2191] = 0.
-    compare[:, :, :, 3108] = 0.
-    compare[:, :, :, 3109] = 0.
+    compare = np.ones((2, 3, 5, 3200))
+    compare[..., 648: 648+4] = 0.
+    compare[..., 668: 668+4] = 0.
+    compare[..., 690: 690+4] = 0.
+    compare[..., 710: 710+4] = 0.
+    compare[..., 1890: 1890+4] = 0.
+    compare[..., 1910: 1910+4] = 0.
+    compare[..., 1808: 1808+4] = 0.
+    compare[..., 1828: 1828+4] = 0.
+    compare[..., 2028: 2028+4] = 0.
+    compare[..., 2068: 2068+4] = 0.
+    compare[..., 2150: 2150+4] = 0.
+    compare[..., 2190: 2190+4] = 0.
 
     assert np.allclose(data, compare)
 
@@ -62,12 +48,12 @@ def test_decode_mask():
     check = np.zeros(nrows, dtype=bool)
     compare = [[5, 20],
                [],
-               [18, 23, 27],
-               [28],
+               [4, 8, 13],
+               [3],
                [6, 14],
                [9, 25],
-               [28, 31],
-               [16, 24, 31]]
+               [0, 3],
+               [0, 7, 15]]
     for row in range(nrows):
         bits = decode_mask(output[row], mask[row])
         check[row] = (bits == compare[row])
diff --git a/jwst/refpix/tests/test_rm_intermittent_badpix.py b/jwst/refpix/tests/test_rm_intermittent_badpix.py
index 29c4905d9..544cd756b 100644
--- a/jwst/refpix/tests/test_rm_intermittent_badpix.py
+++ b/jwst/refpix/tests/test_rm_intermittent_badpix.py
@@ -16,7 +16,7 @@ def test_rm_intermittent_badpix():
     data[..., 3128] = 15.
 
     scipix_n, refpix_r = 16, 4
-    ovr_corr_mitigation_ftr = 1.8
+    ovr_corr_mitigation_ftr = 3.0
     rm_intermittent_badpix(data, scipix_n, refpix_r, ovr_corr_mitigation_ftr)
 
     compare = np.ones((2, 3, 5, 3200), dtype=np.float32)