From f1314ae0cdcb9a2826cdcc53f0fd7d7653eb5563 Mon Sep 17 00:00:00 2001 From: David Law Date: Fri, 7 Feb 2025 10:26:35 -0500 Subject: [PATCH 1/5] Fix flagged groups counting bug --- src/stcal/jump/twopoint_difference.py | 39 ++++++++++++++++++--------- 1 file changed, 27 insertions(+), 12 deletions(-) diff --git a/src/stcal/jump/twopoint_difference.py b/src/stcal/jump/twopoint_difference.py index 511761bb..a38e87cf 100644 --- a/src/stcal/jump/twopoint_difference.py +++ b/src/stcal/jump/twopoint_difference.py @@ -225,24 +225,39 @@ def find_crs_old( num_flagged_grps = 0 # determine the number of groups with all pixels set to DO_NOT_USE ngrps = dat.shape[1] + max_flagged_grps = 0 + total_flagged_grps = 0 for integ in range(nints): + num_flagged_grps = 0 for grp in range(dat.shape[1]): if np.all(np.bitwise_and(gdq[integ, grp, :, :], dnu_flag)): num_flagged_grps += 1 + if num_flagged_grps > max_flagged_grps: + max_flagged_grps = num_flagged_grps + total_flagged_grps += num_flagged_grps if only_use_ints: total_sigclip_groups = nints else: - total_sigclip_groups = nints * (ngrps - num_flagged_grps) - total_groups = nints * (ngrps - num_flagged_grps) - total_diffs = nints * (ngrps - 1 - num_flagged_grps) - total_usable_diffs = total_diffs - num_flagged_grps - if ((ngrps < minimum_groups and only_use_ints and nints < minimum_sigclip_groups) or - (not only_use_ints and nints * ngrps < minimum_sigclip_groups and - total_groups < minimum_groups)): + total_sigclip_groups = nints * ngrps - num_flagged_grps + + min_usable_groups = ngrps - max_flagged_grps + total_groups = nints * ngrps - total_flagged_grps + min_usable_diffs = min_usable_groups - 1 + sig_clip_grps_fails = False + total_noise_min_grps_fails = False + + # Determine whether there are enough usable groups the two sigma clip options + if (only_use_ints and nints < minimum_sigclip_groups) \ + or \ + (not only_use_ints and total_sigclip_groups < minimum_sigclip_groups): + sig_clip_grps_fails = True + if min_usable_groups < minimum_groups: + total_noise_min_grps_fails = True + + if total_noise_min_grps_fails and sig_clip_grps_fails: log.info("Jump Step was skipped because exposure has less than the minimum number of usable groups") - dummy = np.zeros((dataa.shape[1] - 1, dataa.shape[2], dataa.shape[3]), - dtype=np.float32) - return gdq, row_below_gdq, row_above_gdq, 0, dummy + dummy = np.zeros((ngroups - 1, nrows, ncols), dtype=np.float32) + return gdq, row_below_gdq, row_above_gdq, -99, dummy else: # set 'saturated' or 'do not use' pixels to nan in data dat[gdq & (dnu_flag | sat_flag) != 0] = np.nan @@ -301,8 +316,8 @@ def find_crs_old( warnings.resetwarnings() else: # There are not enough groups for sigma clipping - - if total_usable_diffs >= min_diffs_single_pass: + if min_usable_diffs >= min_diffs_single_pass: + # There are enough diffs in all ints to look for more than one jump # compute 'ratio' for each group. this is the value that will be # compared to 'threshold' to classify jumps. subtract the median of From b39c96d4cca6de6acbcf5bdbe7f9fa0e24bd723a Mon Sep 17 00:00:00 2001 From: David Law Date: Fri, 7 Feb 2025 10:38:06 -0500 Subject: [PATCH 2/5] Add unit tests --- tests/test_jump_twopoint_difference.py | 92 ++++++++++++++++++++++++++ 1 file changed, 92 insertions(+) diff --git a/tests/test_jump_twopoint_difference.py b/tests/test_jump_twopoint_difference.py index ae95fa78..91e3c56d 100644 --- a/tests/test_jump_twopoint_difference.py +++ b/tests/test_jump_twopoint_difference.py @@ -22,6 +22,98 @@ def _cube(ngroups, nints=1, nrows=204, ncols=204, readnoise=10): return _cube +def test_sigclip_not_enough_groups(setup_cube): + ngroups = 10 + nints = 9 + nrows = 200 + ncols = 200 + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=10, minimum_groups=5 + ) + assert total_crs == -99 + +def test_sigclip_not_enough_groups(setup_cube): + ngroups = 5 + nints = 9 + nrows = 200 + ncols = 200 + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=11, minimum_groups=6 + ) + assert total_crs == -99 + +def test_sigclip_not_enough_groups(setup_cube): + ngroups = 12 + nints = 9 + nrows = 200 + ncols = 200 + + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) + data[0, 0:2, 0:, 0] = 1 + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=11, minimum_groups=16 + ) + assert total_crs == -99 + +def test_sigclip_with_num_small_groups(setup_cube): + ngroups = 12 + nints = 190 + + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=nints, nrows=2, ncols=2, readnoise=8) + data[0, 0:2, 0:, 0] = 1 + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=11, minimum_groups=16 + ) + assert total_crs != -99 + +def test_nosigclip_with_enough_groups(setup_cube): + ngroups = 12 + nints = 1 + + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=nints, nrows=2, ncols=2, readnoise=8) + data[0, 0:2, 0:, 0] = 1 + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=11, minimum_groups=5 + ) + assert total_crs != -99 + +def test_nosigclip_with_num_small_groups(setup_cube): + ngroups = 4 + nints = 1 + + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=nints, nrows=2, ncols=2, readnoise=8) + data[0, 0:2, 0:, 0] = 1 + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=11, minimum_groups=16 + ) + assert total_crs == -99 + +def test_nosigclip_with_num_small_groups(setup_cube): + ngroups = 4 + nints = 100 + + data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=nints, nrows=2, ncols=2, readnoise=8) + data[0, 0:2, 0:, 0] = 1 + out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( + data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, + nframes, False, 200, 10, DQFLAGS, + minimum_sigclip_groups=100, minimum_groups=16 + ) + assert total_crs != -99 def test_varying_groups(setup_cube): ngroups = 5 From 723134d2aa694d6a35ee039a894a74c9fcfb25db Mon Sep 17 00:00:00 2001 From: David Law Date: Fri, 7 Feb 2025 10:57:42 -0500 Subject: [PATCH 3/5] Update tests --- tests/test_jump_twopoint_difference.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/tests/test_jump_twopoint_difference.py b/tests/test_jump_twopoint_difference.py index 91e3c56d..491e2550 100644 --- a/tests/test_jump_twopoint_difference.py +++ b/tests/test_jump_twopoint_difference.py @@ -24,9 +24,6 @@ def _cube(ngroups, nints=1, nrows=204, ncols=204, readnoise=10): def test_sigclip_not_enough_groups(setup_cube): ngroups = 10 - nints = 9 - nrows = 200 - ncols = 200 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, @@ -37,9 +34,6 @@ def test_sigclip_not_enough_groups(setup_cube): def test_sigclip_not_enough_groups(setup_cube): ngroups = 5 - nints = 9 - nrows = 200 - ncols = 200 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( data, gdq, read_noise, rej_threshold, rej_threshold, rej_threshold, @@ -50,9 +44,6 @@ def test_sigclip_not_enough_groups(setup_cube): def test_sigclip_not_enough_groups(setup_cube): ngroups = 12 - nints = 9 - nrows = 200 - ncols = 200 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) data[0, 0:2, 0:, 0] = 1 From be438abeee1a15475fbb5de77be651c91f72f23c Mon Sep 17 00:00:00 2001 From: David Law Date: Fri, 7 Feb 2025 11:11:55 -0500 Subject: [PATCH 4/5] Change log entry --- changes/338.bugfix.rst | 1 + 1 file changed, 1 insertion(+) create mode 100644 changes/338.bugfix.rst diff --git a/changes/338.bugfix.rst b/changes/338.bugfix.rst new file mode 100644 index 00000000..48e5bfac --- /dev/null +++ b/changes/338.bugfix.rst @@ -0,0 +1 @@ +Fix a bug in the jump step with counting the number of usable groups. From ec9b74cbe71e59bec467faf94fcf878da6900078 Mon Sep 17 00:00:00 2001 From: David Law Date: Fri, 7 Feb 2025 11:46:01 -0500 Subject: [PATCH 5/5] Unique test names --- tests/test_jump_twopoint_difference.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/test_jump_twopoint_difference.py b/tests/test_jump_twopoint_difference.py index 491e2550..8bd62b74 100644 --- a/tests/test_jump_twopoint_difference.py +++ b/tests/test_jump_twopoint_difference.py @@ -22,7 +22,7 @@ def _cube(ngroups, nints=1, nrows=204, ncols=204, readnoise=10): return _cube -def test_sigclip_not_enough_groups(setup_cube): +def test_sigclip_not_enough_groups_ng10(setup_cube): ngroups = 10 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( @@ -32,7 +32,7 @@ def test_sigclip_not_enough_groups(setup_cube): ) assert total_crs == -99 -def test_sigclip_not_enough_groups(setup_cube): +def test_sigclip_not_enough_groups_ng5(setup_cube): ngroups = 5 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) out_gdq, row_below_gdq, rows_above_gdq, total_crs, stddev = find_crs( @@ -42,7 +42,7 @@ def test_sigclip_not_enough_groups(setup_cube): ) assert total_crs == -99 -def test_sigclip_not_enough_groups(setup_cube): +def test_sigclip_not_enough_groups_ng12(setup_cube): ngroups = 12 data, gdq, nframes, read_noise, rej_threshold = setup_cube(ngroups, nints=8, nrows=2, ncols=2, readnoise=8) @@ -54,7 +54,7 @@ def test_sigclip_not_enough_groups(setup_cube): ) assert total_crs == -99 -def test_sigclip_with_num_small_groups(setup_cube): +def test_sigclip_with_num_small_groups_ni190(setup_cube): ngroups = 12 nints = 190 @@ -67,7 +67,7 @@ def test_sigclip_with_num_small_groups(setup_cube): ) assert total_crs != -99 -def test_nosigclip_with_enough_groups(setup_cube): +def test_nosigclip_with_enough_groups_ni1(setup_cube): ngroups = 12 nints = 1 @@ -80,7 +80,7 @@ def test_nosigclip_with_enough_groups(setup_cube): ) assert total_crs != -99 -def test_nosigclip_with_num_small_groups(setup_cube): +def test_nosigclip_with_num_small_groups_ni1(setup_cube): ngroups = 4 nints = 1 @@ -93,7 +93,7 @@ def test_nosigclip_with_num_small_groups(setup_cube): ) assert total_crs == -99 -def test_nosigclip_with_num_small_groups(setup_cube): +def test_nosigclip_with_num_small_groups_ni100(setup_cube): ngroups = 4 nints = 100