Merge pull request #74 from cshanahan1/fix_jump

Jump step: fix median, add tests, add logging

CHANGES.rst

@@ -1,11 +1,18 @@
-0.6.0 (unreleased)
-==================
+0.6.0 (22-01-14)
+================
 
 ramp_fitting
 ------------
 
 - Adding GLS code back to ramp fitting. [#64]
 
+jump
+----
+
+- Fix issue in jump detection that occured when there were only 2 usable
+  differences with no other groups flagged. This PR also added tests and
+  fixed some of the logging statements in twopoint difference. [#74]
+
 0.5.1 (2022-01-07)
 ==================
 

src/stcal/jump/twopoint_difference.py

@@ -1,5 +1,9 @@
+import logging
 import numpy as np
 
+log = logging.getLogger(__name__)
+log.setLevel(logging.DEBUG)
+
 
 def find_crs(dataa, group_dq, read_noise, normal_rej_thresh,
              two_diff_rej_thresh, three_diff_rej_thresh, nframes,
@@ -80,6 +84,8 @@ def find_crs(dataa, group_dq, read_noise, normal_rej_thresh,
 
     for integ in range(nints):
 
+        log.info(f'Working on integration {integ + 1}:')
+
         # get data, gdq for this integration
         dat = dataa[integ]
         gdq_integ = gdq[integ]
@@ -119,6 +125,11 @@ def find_crs(dataa, group_dq, read_noise, normal_rej_thresh,
         row2cr, col2cr = np.where(np.logical_and(ndiffs - num_unusable_groups == 2,
                                   max_ratio > two_diff_rej_thresh))
 
+        log_str = 'From highest outlier, two-point found {} pixels with at least one CR from {} groups.'
+        log.info(log_str.format(len(row4cr), 'five'))
+        log.info(log_str.format(len(row3cr), 'four'))
+        log.info(log_str.format(len(row2cr), 'three'))
+
         # get the rows, col pairs for all pixels with at least one CR
         all_crs_row = np.concatenate((row4cr, row3cr, row2cr))
         all_crs_col = np.concatenate((col4cr, col3cr, col2cr))
@@ -263,8 +274,8 @@ def calc_med_first_diffs(first_diffs):
             return np.nanmedian(first_diffs[mask])
         elif num_usable_groups == 3:  # if 3, no clipping just return median
             return(np.nanmedian(first_diffs))
-        elif num_usable_groups == 2:  # if 2, return minimum
-            return(np.nanmin(np.abs(first_diffs)))
+        elif num_usable_groups == 2:  # if 2, return diff with minimum abs
+            return(first_diffs[np.nanargmin(np.abs(first_diffs))])
         else:
             return(np.nan)
 
@@ -276,22 +287,25 @@ def calc_med_first_diffs(first_diffs):
 
     # process groups with >=4 usable groups
     row4, col4 = np.where(num_usable_groups >= 4)  # locations of >= 4 usable group pixels
-
-    four_slice = first_diffs[:, row4, col4]
-    four_slice[np.nanargmax(np.abs(four_slice), axis=0),
-               np.arange(four_slice.shape[1])] = np.nan  # mask largest group in slice
-    median_diffs[row4, col4] = np.nanmedian(four_slice, axis=0)  # add median to return arr for these pix
+    if len(row4) > 0:
+        four_slice = first_diffs[:, row4, col4]
+        four_slice[np.nanargmax(np.abs(four_slice), axis=0),
+                   np.arange(four_slice.shape[1])] = np.nan  # mask largest group in slice
+        median_diffs[row4, col4] = np.nanmedian(four_slice, axis=0)  # add median to return arr for these pix
 
     # process groups with 3 usable groups
     row3, col3 = np.where(num_usable_groups == 3)  # locations of >= 4 usable group pixels
-    three_slice = first_diffs[:, row3, col3]
-    median_diffs[row3, col3] = np.nanmedian(three_slice, axis=0)  # add median to return arr for these pix
+    if len(row3) > 0:
+        three_slice = first_diffs[:, row3, col3]
+        median_diffs[row3, col3] = np.nanmedian(three_slice, axis=0)  # add median to return arr for these pix
 
     # process groups with 2 usable groups
     row2, col2 = np.where(num_usable_groups == 2)  # locations of >= 4 usable group pixels
     if len(row2) > 0:
         two_slice = first_diffs[:, row2, col2]
-        median_diffs[row2, col2] = np.nanmin(two_slice)  # add median to return arr for these pix
+        two_slice[np.nanargmax(np.abs(two_slice), axis=0),
+                  np.arange(two_slice.shape[1])] = np.nan  # mask larger abs. val
+        median_diffs[row2, col2] = np.nanmin(two_slice, axis=0)  # add med. to return arr
 
     # set the medians all groups with less than 2 usable groups to nan to skip further
     # calculations for these pixels

tests/test_twopoint_difference.py

@@ -1,7 +1,7 @@
 import pytest
 import numpy as np
 
-from stcal.jump.twopoint_difference import find_crs
+from stcal.jump.twopoint_difference import find_crs, calc_med_first_diffs
 
 
 DQFLAGS = {'JUMP_DET': 4, 'SATURATED': 2, 'DO_NOT_USE': 1}
@@ -833,3 +833,75 @@ def test_first_last_3group(setup_cube):
     assert outgdq[0, 0, 0, 0] == 0
     assert outgdq[0, 2, 0, 0] == DQFLAGS['DO_NOT_USE']
     assert outgdq[0, 1, 0, 0] == 0
+
+
+def test_median_func():
+
+    """
+      Test the function `calc_med_first_diffs` that computes median of pixels.
+      Ensure that the correct treatment based on number of non-nan diffs
+      is being done, and that it works for individual pixels as well as
+      pixels embedded in 3d arrays, and that it works for arrays with or
+      without nans (which represent masked pixels)."""
+
+    # single pix with 5 good diffs, should clip 1 pix and return median
+    # 1d, no nans
+    arr = np.array([1., 2., 3., 4., 5])
+    assert calc_med_first_diffs(arr) == 2.5
+    # 3d array, no nans
+    arr = np.zeros(5 * 2 * 2).reshape(5, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3., 4., 5])
+    assert calc_med_first_diffs(arr)[0, 0] == 2.5
+    # 1d, with nans
+    arr = np.array([1., 2., 3., np.nan, 4., 5, np.nan])
+    assert calc_med_first_diffs(arr) == 2.5
+    # 3d, with nans
+    arr = np.zeros(7 * 2 * 2).reshape(7, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3., np.nan, 4., 5, np.nan])
+    assert calc_med_first_diffs(arr)[0, 0] == 2.5
+
+    # single pix with exactly 4 good diffs, should also clip 1 pix and return median
+    # 1d, no nans
+    arr = np.array([1., 2., 3., 4.])
+    assert calc_med_first_diffs(arr) == 2
+    # 3d array, no nans
+    arr = np.zeros(4 * 2 * 2).reshape(4, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3., 4.])
+    assert calc_med_first_diffs(arr)[0, 0] == 2
+    # 1d, with nans
+    arr = np.array([1., 2., 3., np.nan, 4., np.nan])
+    assert calc_med_first_diffs(arr) == 2
+    # 3d, with nans
+    arr = np.zeros(6 * 2 * 2).reshape(6, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3., np.nan, 4., np.nan])
+    assert calc_med_first_diffs(arr)[0, 0] == 2
+
+    # single pix with exactly 3 good diffs, should compute median without clipping
+    arr = np.array([1., 2., 3.])
+    assert calc_med_first_diffs(arr) == 2
+    # 3d array, no nans
+    arr = np.zeros(3 * 2 * 2).reshape(3, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3.])
+    assert calc_med_first_diffs(arr)[0, 0] == 2
+    # 1d, with nans
+    arr = np.array([1., 2., 3., np.nan, np.nan])
+    assert calc_med_first_diffs(arr) == 2
+    # 3d, with nans
+    arr = np.zeros(5 * 2 * 2).reshape(5, 2, 2)
+    arr[:, 0, 0] = np.array([1., 2., 3., np.nan, np.nan])
+    assert calc_med_first_diffs(arr)[0, 0] == 2
+
+    # # single pix with exactly 2 good diffs, should return the element with the minimum abs val
+    arr = np.array([-1., -2.])
+    assert calc_med_first_diffs(arr) == -1
+    # 3d array, no nans
+    arr = np.zeros(2 * 2 * 2).reshape(2, 2, 2)
+    arr[:, 0, 0] = np.array([-1., -2.])
+    assert calc_med_first_diffs(arr)[0, 0] == -1
+    # 1d, with nans
+    arr = np.array([-1., -2., np.nan, np.nan])
+    assert calc_med_first_diffs(arr) == -1
+    # 3d, with nans
+    arr = np.zeros(4 * 2 * 2).reshape(4, 2, 2)
+    arr[:, 0, 0] = np.array([-1., -2., np.nan, np.nan])
+    assert calc_med_first_diffs(arr)[0, 0] == -1