From 06387ddf634253c6a089dd33b7545547bc5a9e6c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 13 Dec 2023 17:38:45 -0500 Subject: [PATCH 01/40] test --- src/stcal/jump/jump.py | 38 +++++++++++++++++++++++++++++-- tests/test_jump.py | 18 +++++++++++++++ tests/test_twopoint_difference.py | 5 ++++ 3 files changed, 59 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 6139d103..5a2ec5ac 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1,11 +1,20 @@ import logging import multiprocessing import time +<<<<<<< Updated upstream import cv2 as cv import numpy as np from astropy import stats from astropy.convolution import Ring2DKernel, convolve +======= +import numpy as np +import cv2 as cv +import astropy.stats as stats +from astropy.io import fits +from astropy.convolution import Ring2DKernel +from astropy.convolution import convolve +>>>>>>> Stashed changes from . import constants from . import twopoint_difference as twopt @@ -471,7 +480,7 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - + print("jumpgdq.fits", gdq, overwrite=True) # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev @@ -842,6 +851,7 @@ def find_faint_extended( all_ellipses = [] first_diffs = np.diff(data, axis=1) first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) + fits.writeto("first_diffs_masked.fits", first_diffs_masked.filled(np.nan), overwrite=True) nints = data.shape[0] if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) @@ -849,16 +859,19 @@ def find_faint_extended( # calculate sigma for each pixel if nints <= minimum_sigclip_groups: median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) + fits.writeto("median_diffs.fits", median_diffs, overwrite=True) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) # The difference from the median difference for each group e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - + fits.writeto("e_jump.fits", e_jump.filled(np.nan), overwrite=True) + fits.writeto("ratio.fits", ratio.filled(np.nan), overwrite=True) # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] for grp in range(1, ngrps): + print("starting group ", grp) if nints > minimum_sigclip_groups: median_diffs = median[grp - 1] sigma = stddev[grp - 1] @@ -876,11 +889,15 @@ def find_faint_extended( # mask pix. that are already flagged as sat. masked_ratio[saty, satx] = np.nan masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) + if grp == 2: + fits.writeto("masked_ratio.fits", masked_ratio.filled(np.nan), overwrite=True) + fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True) nrows = ratio.shape[1] ncols = ratio.shape[2] extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 + # find the contours of the extended emission contours, hierarchy = cv.findContours(extended_emission, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE) # get the contours that are above the minimum size @@ -888,7 +905,13 @@ def find_faint_extended( # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] +<<<<<<< Updated upstream +======= + if grp == 2: + fits.writeto("extended_emmission.fits", extended_emission, overwrite=True) + print("number of ellises",len(ellipses)) +>>>>>>> Stashed changes expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] @@ -942,6 +965,7 @@ def find_faint_extended( intg = showers[0] grp = showers[1] ellipses = showers[2] +<<<<<<< Updated upstream gdq, num = extend_ellipses( gdq, intg, @@ -954,6 +978,16 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) +======= + gdq, num = extend_ellipses(gdq, intg, grp, ellipses, sat_flag, + jump_flag, expansion=ellipse_expand, + expand_by_ratio=True, + num_grps_masked=num_grps_masked, + max_extended_radius=max_extended_radius) + if grp == 10: + fits.writeto("gdq10.fits", gdq, overwrite=True) + fits.writeto("gdqall.fits", gdq, overwrite=True) +>>>>>>> Stashed changes return gdq, len(all_ellipses) diff --git a/tests/test_jump.py b/tests/test_jump.py index 0ddbefb1..c2fc54da 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -1,5 +1,6 @@ import numpy as np import pytest +<<<<<<< Updated upstream from stcal.jump.jump import ( calc_num_slices, @@ -9,6 +10,11 @@ flag_large_events, point_inside_ellipse, ) +======= +from astropy.io import fits +from stcal.jump.jump import flag_large_events, find_ellipses, extend_saturation, \ + point_inside_ellipse, find_faint_extended, calc_num_slices +>>>>>>> Stashed changes DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} @@ -29,6 +35,18 @@ def _cube(ngroups, readnoise=10): return _cube +def test_slowmode(): + hdul = fits.open("miri_1264_00_jump.fits") + data = hdul["sci"].data + gdq = hdul["groupdq"].data + readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") + gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, + snr_threshold=1.3, + min_shower_area=20, inner=1, + outer=2, sat_flag=2, jump_flag=4, + ellipse_expand=1.1, num_grps_masked=3) + fits.writeto("outgdq.fits", gdq, overwrite=True) + print("number of showers", num_showers) def test_find_simple_ellipse(): plane = np.zeros(shape=(5, 5), dtype=np.uint8) diff --git a/tests/test_twopoint_difference.py b/tests/test_twopoint_difference.py index c6443bc7..87bb79d3 100644 --- a/tests/test_twopoint_difference.py +++ b/tests/test_twopoint_difference.py @@ -1,9 +1,14 @@ import numpy as np +<<<<<<< Updated upstream import pytest from stcal.jump.twopoint_difference import calc_med_first_diffs, find_crs DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1} +======= +from astropy.io import fits +from stcal.jump.twopoint_difference import find_crs, calc_med_first_diffs +>>>>>>> Stashed changes @pytest.fixture() From 88adb7364ad287958f386a63ebd930b667bc23df Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 10:41:00 -0500 Subject: [PATCH 02/40] fix conflicts --- src/stcal/jump/jump.py | 21 +-------------------- 1 file changed, 1 insertion(+), 20 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 5a2ec5ac..7d584773 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -1,20 +1,13 @@ import logging import multiprocessing import time -<<<<<<< Updated upstream - -import cv2 as cv -import numpy as np -from astropy import stats -from astropy.convolution import Ring2DKernel, convolve -======= import numpy as np import cv2 as cv import astropy.stats as stats from astropy.io import fits from astropy.convolution import Ring2DKernel from astropy.convolution import convolve ->>>>>>> Stashed changes + from . import constants from . import twopoint_difference as twopt @@ -905,13 +898,9 @@ def find_faint_extended( # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] -<<<<<<< Updated upstream - -======= if grp == 2: fits.writeto("extended_emmission.fits", extended_emission, overwrite=True) print("number of ellises",len(ellipses)) ->>>>>>> Stashed changes expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] @@ -965,7 +954,6 @@ def find_faint_extended( intg = showers[0] grp = showers[1] ellipses = showers[2] -<<<<<<< Updated upstream gdq, num = extend_ellipses( gdq, intg, @@ -978,16 +966,9 @@ def find_faint_extended( num_grps_masked=num_grps_masked, max_extended_radius=max_extended_radius, ) -======= - gdq, num = extend_ellipses(gdq, intg, grp, ellipses, sat_flag, - jump_flag, expansion=ellipse_expand, - expand_by_ratio=True, - num_grps_masked=num_grps_masked, - max_extended_radius=max_extended_radius) if grp == 10: fits.writeto("gdq10.fits", gdq, overwrite=True) fits.writeto("gdqall.fits", gdq, overwrite=True) ->>>>>>> Stashed changes return gdq, len(all_ellipses) From 7a9ceb0df143ee1befa071ff740eb25b430e3880 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 11:42:06 -0500 Subject: [PATCH 03/40] Update jump.py --- src/stcal/jump/jump.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 7d584773..07521b5c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -473,8 +473,7 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - print("jumpgdq.fits", gdq, overwrite=True) - # Return the updated data quality arrays + # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev From 4bcd4461f3c6c178aa61c92bcaa876cc7da240c4 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 11:49:25 -0500 Subject: [PATCH 04/40] Update jump.py --- src/stcal/jump/jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 07521b5c..7c6c2811 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -473,6 +473,7 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d + print("good version") # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev From 183305b18b493f2879d27fe44b7aa4e308c3fccf Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 11:54:46 -0500 Subject: [PATCH 05/40] Update jump.py --- src/stcal/jump/jump.py | 1 - 1 file changed, 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 7c6c2811..07521b5c 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -473,7 +473,6 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - print("good version") # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev From 5d8d6d8102dd5e4926882029242825e1e2fca0b2 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 12:46:43 -0500 Subject: [PATCH 06/40] tests2 --- src/stcal/jump/jump.py | 7 ++++++- tests/test_jump.py | 5 ----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 07521b5c..c936cc8e 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -463,6 +463,7 @@ def detect_jumps( num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, ) + print("after find showers multiproc") log.info("Total showers= %i", num_showers) number_extended_events = num_showers elapsed = time.time() - start @@ -473,6 +474,8 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d + print("end of jump.py") + fits.writeto("finalgdq.fits", gdq, overwrite=True) # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev @@ -949,10 +952,12 @@ def find_faint_extended( # This is deferred until all showers are detected. because the showers # can flag future groups and would confuse the detection algorithm if # we worked on groups that already had some flagged showers. + total_showers = 0 for showers in all_ellipses: intg = showers[0] grp = showers[1] ellipses = showers[2] + total_showers += len(ellipses) gdq, num = extend_ellipses( gdq, intg, @@ -968,7 +973,7 @@ def find_faint_extended( if grp == 10: fits.writeto("gdq10.fits", gdq, overwrite=True) fits.writeto("gdqall.fits", gdq, overwrite=True) - return gdq, len(all_ellipses) + return gdq, total_showers def calc_num_slices(n_rows, max_cores, max_available): diff --git a/tests/test_jump.py b/tests/test_jump.py index c2fc54da..17debc31 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -1,6 +1,5 @@ import numpy as np import pytest -<<<<<<< Updated upstream from stcal.jump.jump import ( calc_num_slices, @@ -10,11 +9,7 @@ flag_large_events, point_inside_ellipse, ) -======= from astropy.io import fits -from stcal.jump.jump import flag_large_events, find_ellipses, extend_saturation, \ - point_inside_ellipse, find_faint_extended, calc_num_slices ->>>>>>> Stashed changes DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} From c099e12976d9f3b5448b5e5c6468a988dacbb810 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 13:05:19 -0500 Subject: [PATCH 07/40] Update jump.py --- src/stcal/jump/jump.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index c936cc8e..3f7694d1 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -889,6 +889,7 @@ def find_faint_extended( fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True) nrows = ratio.shape[1] ncols = ratio.shape[2] + print("snr_threshold", snr_threshold) extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 @@ -902,7 +903,7 @@ def find_faint_extended( ellipses = [cv.minAreaRect(con) for con in bigcontours] if grp == 2: fits.writeto("extended_emmission.fits", extended_emission, overwrite=True) - print("number of ellises",len(ellipses)) + print("number of ellipses", len(ellipses)) expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] From a7f70a4236102c54d757fa90edcaceb39da360ce Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 13:08:59 -0500 Subject: [PATCH 08/40] changes --- src/stcal/jump/jump.py | 2 +- tests/test_jump.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 3f7694d1..e6e9a7d9 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -903,7 +903,7 @@ def find_faint_extended( ellipses = [cv.minAreaRect(con) for con in bigcontours] if grp == 2: fits.writeto("extended_emmission.fits", extended_emission, overwrite=True) - print("number of ellipses", len(ellipses)) + print("number of ellipses in grp", len(ellipses)) expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] diff --git a/tests/test_jump.py b/tests/test_jump.py index 17debc31..4ad80fbb 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -41,7 +41,7 @@ def test_slowmode(): outer=2, sat_flag=2, jump_flag=4, ellipse_expand=1.1, num_grps_masked=3) fits.writeto("outgdq.fits", gdq, overwrite=True) - print("number of showers", num_showers) + print("total number of showers", num_showers) def test_find_simple_ellipse(): plane = np.zeros(shape=(5, 5), dtype=np.uint8) From 5a5db9494c06246371a03685c6e45e25bc60cace Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 13:14:49 -0500 Subject: [PATCH 09/40] Update jump.py --- src/stcal/jump/jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index e6e9a7d9..78db145f 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -867,7 +867,9 @@ def find_faint_extended( ngrps = data.shape[1] for grp in range(1, ngrps): print("starting group ", grp) + print("snr_threshold 1", snr_threshold) if nints > minimum_sigclip_groups: + print("inside sigclipped groups") median_diffs = median[grp - 1] sigma = stddev[grp - 1] # The difference from the median difference for each group From 5ddb6c867ff389c3391d738d79ffd3c39cc17fa6 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 13:48:19 -0500 Subject: [PATCH 10/40] tests --- src/stcal/jump/jump.py | 2 +- tests/test_jump.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 78db145f..5be32302 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -838,7 +838,7 @@ def find_faint_extended( Total number of showers detected. """ - read_noise_2 = readnoise_2d**2 + read_noise_2 = readnoise_2d data = indata.copy() data[gdq == sat_flag] = np.nan data[gdq == 1] = np.nan diff --git a/tests/test_jump.py b/tests/test_jump.py index 4ad80fbb..7617027f 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -34,7 +34,7 @@ def test_slowmode(): hdul = fits.open("miri_1264_00_jump.fits") data = hdul["sci"].data gdq = hdul["groupdq"].data - readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") + readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") * 3.9 gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, snr_threshold=1.3, min_shower_area=20, inner=1, From 4569b111120758139088637508d636f3054f9a3d Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 14:10:52 -0500 Subject: [PATCH 11/40] fixes --- src/stcal/jump/jump.py | 6 ++++-- tests/test_jump.py | 2 +- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 5be32302..b18983cd 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -838,8 +838,10 @@ def find_faint_extended( Total number of showers detected. """ - read_noise_2 = readnoise_2d + read_noise_2d_sqr = readnoise_2d**2 + print("median readnoise sqr", np.nanmedian(read_noise_2d_sqr)) data = indata.copy() + print("median readnoise data", np.nanmedian(data)) data[gdq == sat_flag] = np.nan data[gdq == 1] = np.nan data[gdq == jump_flag] = np.nan @@ -855,7 +857,7 @@ def find_faint_extended( if nints <= minimum_sigclip_groups: median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) fits.writeto("median_diffs.fits", median_diffs, overwrite=True) - sigma = np.sqrt(np.abs(median_diffs) + read_noise_2 / nframes) + sigma = np.sqrt(np.abs(median_diffs) + read_noise_2d_sqr / nframes) # The difference from the median difference for each group e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] # SNR ratio of each diff. diff --git a/tests/test_jump.py b/tests/test_jump.py index 7617027f..d5b991c0 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -32,7 +32,7 @@ def _cube(ngroups, readnoise=10): def test_slowmode(): hdul = fits.open("miri_1264_00_jump.fits") - data = hdul["sci"].data + data = hdul["sci"].data * 3.9 gdq = hdul["groupdq"].data readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") * 3.9 gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, From f5fed788da40f07a1745097e8a18a221fb6e7642 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:10:46 -0500 Subject: [PATCH 12/40] changes --- src/stcal/jump/jump.py | 8 ++++---- tests/test_jump.py | 9 +++++---- 2 files changed, 9 insertions(+), 8 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index b18983cd..6cd87a99 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -838,10 +838,11 @@ def find_faint_extended( Total number of showers detected. """ + fits.writeto("gdqin.fits", gdq, overwrite=True) read_noise_2d_sqr = readnoise_2d**2 - print("median readnoise sqr", np.nanmedian(read_noise_2d_sqr)) + print("median readnoise sqr in e-", np.nanmedian(read_noise_2d_sqr)) data = indata.copy() - print("median readnoise data", np.nanmedian(data)) + print("median data in e=", np.nanmedian(data)) data[gdq == sat_flag] = np.nan data[gdq == 1] = np.nan data[gdq == jump_flag] = np.nan @@ -869,9 +870,8 @@ def find_faint_extended( ngrps = data.shape[1] for grp in range(1, ngrps): print("starting group ", grp) - print("snr_threshold 1", snr_threshold) if nints > minimum_sigclip_groups: - print("inside sigclipped groups") + print("inside sig clipped groups") median_diffs = median[grp - 1] sigma = stddev[grp - 1] # The difference from the median difference for each group diff --git a/tests/test_jump.py b/tests/test_jump.py index d5b991c0..b18a69c0 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -32,13 +32,14 @@ def _cube(ngroups, readnoise=10): def test_slowmode(): hdul = fits.open("miri_1264_00_jump.fits") - data = hdul["sci"].data * 3.9 + gain = fits.getdata("jwst_miri_gain_0026.fits") + data = hdul["sci"].data * gain gdq = hdul["groupdq"].data - readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") * 3.9 + readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") * gain gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, snr_threshold=1.3, - min_shower_area=20, inner=1, - outer=2, sat_flag=2, jump_flag=4, + min_shower_area=90, inner=1, + outer=2.6, sat_flag=2, jump_flag=4, ellipse_expand=1.1, num_grps_masked=3) fits.writeto("outgdq.fits", gdq, overwrite=True) print("total number of showers", num_showers) From 1251f170bd3d367cfc431ca3d22121acfd43edab Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:26:24 -0500 Subject: [PATCH 13/40] Update jump.py --- src/stcal/jump/jump.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 6cd87a99..aa834dc6 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -635,7 +635,8 @@ def extend_ellipses( # For a given DQ plane it will use the list of ellipses to create # expanded ellipses of pixels with # the jump flag set. - plane = gdq_cube[intg, grp, :, :] + out_gdq_cube = gdq_cube.copy() + plane = gdq_cube[intg, grp, :, :].copy() ncols = plane.shape[1] nrows = plane.shape[0] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) @@ -681,8 +682,8 @@ def extend_ellipses( sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 - gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) - return gdq_cube, num_ellipses + out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) + return out_gdq_cube, num_ellipses def find_circles(dqplane, bitmask, min_area): From affff199a14d54e1e80c9cf754d497d572a5101f Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:30:57 -0500 Subject: [PATCH 14/40] Update jump.py --- src/stcal/jump/jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index aa834dc6..46c1944f 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -683,6 +683,8 @@ def extend_ellipses( saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) + diff_cube = out_gdq_cube - gdq_cube + fits.writeto("diff_cube.fits", diff_cube, overwrite=True) return out_gdq_cube, num_ellipses From 0fd250a42ba02fc8bb568e29bd05aa2ed1469f2b Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:44:07 -0500 Subject: [PATCH 15/40] Update jump.py --- src/stcal/jump/jump.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 46c1944f..a6b64e85 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -641,6 +641,7 @@ def extend_ellipses( nrows = plane.shape[0] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) num_ellipses = len(ellipses) + # print("num ellipses", num_ellipses) for ellipse in ellipses: ceny = ellipse[0][0] cenx = ellipse[0][1] @@ -676,12 +677,13 @@ def extend_ellipses( ) jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] - last_grp = min(grp + num_grps_masked, ngrps) + last_grp = min(grp + 1, min(grp + num_grps_masked, ngrps)) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 + print("jump ellipse max", np.max(jump_ellipse)) out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) diff_cube = out_gdq_cube - gdq_cube fits.writeto("diff_cube.fits", diff_cube, overwrite=True) From 22dab7e1343f4f4f5de5e85b6b584bfc12088f35 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:49:17 -0500 Subject: [PATCH 16/40] Update jump.py --- src/stcal/jump/jump.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index a6b64e85..fbcff5cb 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -678,12 +678,13 @@ def extend_ellipses( jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] last_grp = min(grp + 1, min(grp + num_grps_masked, ngrps)) + print("range", grp, last_grp) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 - print("jump ellipse max", np.max(jump_ellipse)) + print("jump ellipse sum", np.sum(jump_ellipse)) out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) diff_cube = out_gdq_cube - gdq_cube fits.writeto("diff_cube.fits", diff_cube, overwrite=True) From 0617d35d0940151fddd62f7034cdcfcc4a5cebf0 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:53:40 -0500 Subject: [PATCH 17/40] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index fbcff5cb..16bb3ecd 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -641,7 +641,7 @@ def extend_ellipses( nrows = plane.shape[0] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) num_ellipses = len(ellipses) - # print("num ellipses", num_ellipses) + print("num ellipses", num_ellipses) for ellipse in ellipses: ceny = ellipse[0][0] cenx = ellipse[0][1] From 6806a7c67eaf083ea3305f685c94edb83e42aa57 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 14 Dec 2023 15:57:06 -0500 Subject: [PATCH 18/40] Update jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 16bb3ecd..773fac62 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -677,7 +677,7 @@ def extend_ellipses( ) jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] - last_grp = min(grp + 1, min(grp + num_grps_masked, ngrps)) + last_grp = max(grp + 1, min(grp + num_grps_masked, ngrps)) print("range", grp, last_grp) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): From 5008db3056d35cb7ecafce5e1d8a3bb460ecd285 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 13:41:30 -0500 Subject: [PATCH 19/40] add tests --- src/stcal/jump/jump.py | 7 ++++++- tests/test_jump.py | 9 +++++++++ 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 773fac62..c9e9c6dd 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -677,7 +677,7 @@ def extend_ellipses( ) jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] - last_grp = max(grp + 1, min(grp + num_grps_masked, ngrps)) + last_grp = find_last_grp(grp, ngrps, num_grps_masked) print("range", grp, last_grp) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): @@ -690,6 +690,11 @@ def extend_ellipses( fits.writeto("diff_cube.fits", diff_cube, overwrite=True) return out_gdq_cube, num_ellipses +def find_last_grp(grp, ngrps, num_grps_masked): + if num_grps_masked == 0: + num_grps_masked = 1 + last_grp = min(grp + num_grps_masked, ngrps) + return last_grp def find_circles(dqplane, bitmask, min_area): # Using an input DQ plane this routine will find the groups of pixels with at least the minimum diff --git a/tests/test_jump.py b/tests/test_jump.py index b18a69c0..4b8fc373 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -8,6 +8,7 @@ find_faint_extended, flag_large_events, point_inside_ellipse, + find_last_grp ) from astropy.io import fits @@ -366,3 +367,11 @@ def test_calc_num_slices(): assert calc_num_slices(n_rows, "3/4", max_available_cores) == 1 n_rows = 9 assert calc_num_slices(n_rows, "21", max_available_cores) == 9 + +def test_find_last_grp(): + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=2) == 7) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=3) == 7) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=1) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=2) == 6) From d7f1513c520548ab0a4dcce64c299efbcd4abc2e Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 13:48:15 -0500 Subject: [PATCH 20/40] Update test_jump.py --- tests/test_jump.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_jump.py b/tests/test_jump.py index 4b8fc373..ae9276e3 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -368,6 +368,7 @@ def test_calc_num_slices(): n_rows = 9 assert calc_num_slices(n_rows, "21", max_available_cores) == 9 + def test_find_last_grp(): assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=0) == 6) assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=2) == 7) From 1d8105bf26a0507b3eb704343b0b3b6ef7bcd90a Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 14:06:01 -0500 Subject: [PATCH 21/40] clean up --- src/stcal/jump/jump.py | 33 +++------------------------------ tests/test_jump.py | 3 +++ 2 files changed, 6 insertions(+), 30 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index c9e9c6dd..43bfa66e 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -4,7 +4,7 @@ import numpy as np import cv2 as cv import astropy.stats as stats -from astropy.io import fits + from astropy.convolution import Ring2DKernel from astropy.convolution import convolve @@ -463,7 +463,6 @@ def detect_jumps( num_grps_masked=grps_masked_after_shower, max_extended_radius=max_extended_radius, ) - print("after find showers multiproc") log.info("Total showers= %i", num_showers) number_extended_events = num_showers elapsed = time.time() - start @@ -474,8 +473,6 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - print("end of jump.py") - fits.writeto("finalgdq.fits", gdq, overwrite=True) # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev @@ -641,7 +638,6 @@ def extend_ellipses( nrows = plane.shape[0] image = np.zeros(shape=(nrows, ncols, 3), dtype=np.uint8) num_ellipses = len(ellipses) - print("num ellipses", num_ellipses) for ellipse in ellipses: ceny = ellipse[0][0] cenx = ellipse[0][1] @@ -678,21 +674,17 @@ def extend_ellipses( jump_ellipse = image[:, :, 2] ngrps = gdq_cube.shape[1] last_grp = find_last_grp(grp, ngrps, num_grps_masked) - print("range", grp, last_grp) # This loop will flag the number of groups for flg_grp in range(grp, last_grp): sat_pix = np.bitwise_and(gdq_cube[intg, flg_grp, :, :], sat_flag) saty, satx = np.where(sat_pix == sat_flag) jump_ellipse[saty, satx] = 0 - print("jump ellipse sum", np.sum(jump_ellipse)) out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) diff_cube = out_gdq_cube - gdq_cube - fits.writeto("diff_cube.fits", diff_cube, overwrite=True) return out_gdq_cube, num_ellipses def find_last_grp(grp, ngrps, num_grps_masked): - if num_grps_masked == 0: - num_grps_masked = 1 + num_grps_masked += 1 last_grp = min(grp + num_grps_masked, ngrps) return last_grp @@ -849,18 +841,14 @@ def find_faint_extended( Total number of showers detected. """ - fits.writeto("gdqin.fits", gdq, overwrite=True) read_noise_2d_sqr = readnoise_2d**2 - print("median readnoise sqr in e-", np.nanmedian(read_noise_2d_sqr)) data = indata.copy() - print("median data in e=", np.nanmedian(data)) data[gdq == sat_flag] = np.nan data[gdq == 1] = np.nan data[gdq == jump_flag] = np.nan all_ellipses = [] first_diffs = np.diff(data, axis=1) first_diffs_masked = np.ma.masked_array(first_diffs, mask=np.isnan(first_diffs)) - fits.writeto("first_diffs_masked.fits", first_diffs_masked.filled(np.nan), overwrite=True) nints = data.shape[0] if nints > minimum_sigclip_groups: mean, median, stddev = stats.sigma_clipped_stats(first_diffs_masked, sigma=5, axis=0) @@ -868,21 +856,16 @@ def find_faint_extended( # calculate sigma for each pixel if nints <= minimum_sigclip_groups: median_diffs = np.nanmedian(first_diffs_masked[intg], axis=0) - fits.writeto("median_diffs.fits", median_diffs, overwrite=True) sigma = np.sqrt(np.abs(median_diffs) + read_noise_2d_sqr / nframes) # The difference from the median difference for each group e_jump = first_diffs_masked[intg] - median_diffs[np.newaxis, :, :] # SNR ratio of each diff. ratio = np.abs(e_jump) / sigma[np.newaxis, :, :] - fits.writeto("e_jump.fits", e_jump.filled(np.nan), overwrite=True) - fits.writeto("ratio.fits", ratio.filled(np.nan), overwrite=True) # The convolution kernel creation ring_2D_kernel = Ring2DKernel(inner, outer) ngrps = data.shape[1] for grp in range(1, ngrps): - print("starting group ", grp) if nints > minimum_sigclip_groups: - print("inside sig clipped groups") median_diffs = median[grp - 1] sigma = stddev[grp - 1] # The difference from the median difference for each group @@ -899,12 +882,8 @@ def find_faint_extended( # mask pix. that are already flagged as sat. masked_ratio[saty, satx] = np.nan masked_smoothed_ratio = convolve(masked_ratio, ring_2D_kernel) - if grp == 2: - fits.writeto("masked_ratio.fits", masked_ratio.filled(np.nan), overwrite=True) - fits.writeto("masked_smoothed_ratio.fits", masked_smoothed_ratio, overwrite=True) nrows = ratio.shape[1] ncols = ratio.shape[2] - print("snr_threshold", snr_threshold) extended_emission = np.zeros(shape=(nrows, ncols), dtype=np.uint8) exty, extx = np.where(masked_smoothed_ratio > snr_threshold) extended_emission[exty, extx] = 1 @@ -916,9 +895,6 @@ def find_faint_extended( # get the minimum enclosing rectangle which is the same as the # minimum enclosing ellipse ellipses = [cv.minAreaRect(con) for con in bigcontours] - if grp == 2: - fits.writeto("extended_emmission.fits", extended_emission, overwrite=True) - print("number of ellipses in grp", len(ellipses)) expand_by_ratio = True expansion = 1.0 plane = gdq[intg, grp, :, :] @@ -984,11 +960,8 @@ def find_faint_extended( expansion=ellipse_expand, expand_by_ratio=True, num_grps_masked=num_grps_masked, - max_extended_radius=max_extended_radius, + max_extended_radius=max_extended_radius ) - if grp == 10: - fits.writeto("gdq10.fits", gdq, overwrite=True) - fits.writeto("gdqall.fits", gdq, overwrite=True) return gdq, total_showers diff --git a/tests/test_jump.py b/tests/test_jump.py index ae9276e3..417f9c97 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -376,3 +376,6 @@ def test_find_last_grp(): assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=1) == 6) assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=0) == 6) assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=2) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=1) == 7) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=2) == 8) From 541ceed416deaa4cc75256ec550f3b68ebebd2aa Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 14:18:02 -0500 Subject: [PATCH 22/40] Update test_jump.py --- tests/test_jump.py | 14 -------------- 1 file changed, 14 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 417f9c97..5068ee95 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -10,7 +10,6 @@ point_inside_ellipse, find_last_grp ) -from astropy.io import fits DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1, "GOOD": 0, "NO_GAIN_VALUE": 8} @@ -31,19 +30,6 @@ def _cube(ngroups, readnoise=10): return _cube -def test_slowmode(): - hdul = fits.open("miri_1264_00_jump.fits") - gain = fits.getdata("jwst_miri_gain_0026.fits") - data = hdul["sci"].data * gain - gdq = hdul["groupdq"].data - readnoise = fits.getdata("jwst_miri_readnoise_0050.fits") * gain - gdq, num_showers = find_faint_extended(data, gdq, readnoise, 1, 100, - snr_threshold=1.3, - min_shower_area=90, inner=1, - outer=2.6, sat_flag=2, jump_flag=4, - ellipse_expand=1.1, num_grps_masked=3) - fits.writeto("outgdq.fits", gdq, overwrite=True) - print("total number of showers", num_showers) def test_find_simple_ellipse(): plane = np.zeros(shape=(5, 5), dtype=np.uint8) From 93d3b56f69b45d1ddacbaba9b7f27dce8f472747 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 14:20:51 -0500 Subject: [PATCH 23/40] Update test_twopoint_difference.py --- tests/test_twopoint_difference.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/tests/test_twopoint_difference.py b/tests/test_twopoint_difference.py index 87bb79d3..c6443bc7 100644 --- a/tests/test_twopoint_difference.py +++ b/tests/test_twopoint_difference.py @@ -1,14 +1,9 @@ import numpy as np -<<<<<<< Updated upstream import pytest from stcal.jump.twopoint_difference import calc_med_first_diffs, find_crs DQFLAGS = {"JUMP_DET": 4, "SATURATED": 2, "DO_NOT_USE": 1} -======= -from astropy.io import fits -from stcal.jump.twopoint_difference import find_crs, calc_med_first_diffs ->>>>>>> Stashed changes @pytest.fixture() From 07badf3e0c659340ad959e531c6f384e306ca98b Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Thu, 28 Dec 2023 14:29:19 -0500 Subject: [PATCH 24/40] Update CHANGES.rst --- CHANGES.rst | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index b22ce3e2..32e0e6d9 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,7 +1,10 @@ 1.5.3 (unreleased) ================== -- +jump +------------ +Fix the code to at least always flag the group with the shower and the requested +groups after the primary shower. [#237] 1.5.2 (2023-12-13) ================== From 28bff55770c78dde9f0df2ff19213bfe6ca1bf8a Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 5 Jan 2024 12:33:09 -0500 Subject: [PATCH 25/40] Update jump.py --- src/stcal/jump/jump.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 43bfa66e..66c72fa6 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -683,9 +683,9 @@ def extend_ellipses( diff_cube = out_gdq_cube - gdq_cube return out_gdq_cube, num_ellipses -def find_last_grp(grp, ngrps, num_grps_masked): - num_grps_masked += 1 - last_grp = min(grp + num_grps_masked, ngrps) +def find_last_grp(grp, ngrps, num_grps_masked_after): + num_grps_masked_after += 1 + last_grp = min(grp + num_grps_masked_after, ngrps) return last_grp def find_circles(dqplane, bitmask, min_area): From b0f06b80c7ffa8ac372232a5f447762ca925c13c Mon Sep 17 00:00:00 2001 From: mwregan2 Date: Fri, 26 Jan 2024 15:30:51 -0500 Subject: [PATCH 26/40] Update CHANGES.rst based on suggestion Co-authored-by: Howard Bushouse --- CHANGES.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGES.rst b/CHANGES.rst index 32e0e6d9..2ea09d4f 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -2,7 +2,7 @@ ================== jump ------------- +---- Fix the code to at least always flag the group with the shower and the requested groups after the primary shower. [#237] From 4f1baf1cd96e66dc937f565849605971893c4608 Mon Sep 17 00:00:00 2001 From: mwregan2 Date: Fri, 26 Jan 2024 15:31:18 -0500 Subject: [PATCH 27/40] Update CHANGES.rst uses suggestion Co-authored-by: Howard Bushouse --- CHANGES.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 2ea09d4f..8281b2e8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -3,8 +3,8 @@ jump ---- -Fix the code to at least always flag the group with the shower and the requested -groups after the primary shower. [#237] +- Fix the code to at least always flag the group with the shower and the requested + groups after the primary shower. [#237] 1.5.2 (2023-12-13) ================== From 7205fdba4b795a13ca69954625a689b1c265c70b Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Fri, 26 Jan 2024 15:31:53 -0500 Subject: [PATCH 28/40] Update jump.py addressed comments in PR --- src/stcal/jump/jump.py | 18 +++++++++++++++++- 1 file changed, 17 insertions(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 66c72fa6..73e5bbe1 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -473,7 +473,7 @@ def detect_jumps( data /= gain_2d err /= gain_2d readnoise_2d /= gain_2d - # Return the updated data quality arrays + # Return the updated data quality arrays return gdq, pdq, total_primary_crs, number_extended_events, stddev @@ -684,6 +684,22 @@ def extend_ellipses( return out_gdq_cube, num_ellipses def find_last_grp(grp, ngrps, num_grps_masked_after): + """ + Parameters + _________ + grp : int + The location of the shower + ngrps : + The number of groups in the integration + num_grps_masked_after : + The requested number of groups to be flagged after the shower + + Returns + _______ + last_grp : int + The index of the last group to flag for the shower + + """ num_grps_masked_after += 1 last_grp = min(grp + num_grps_masked_after, ngrps) return last_grp From 5b15a2420926edf8d7009ab93f74a12fd7e8c53d Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:40:04 -0500 Subject: [PATCH 29/40] Update src/stcal/jump/jump.py --- src/stcal/jump/jump.py | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 73e5bbe1..5cde6522 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -686,18 +686,20 @@ def extend_ellipses( def find_last_grp(grp, ngrps, num_grps_masked_after): """ Parameters - _________ + ---------- grp : int - The location of the shower - ngrps : - The number of groups in the integration - num_grps_masked_after : - The requested number of groups to be flagged after the shower + The location of the shower + + ngrps : int + The number of groups in the integration + + num_grps_masked_after : int + The requested number of groups to be flagged after the shower Returns - _______ + ------- last_grp : int - The index of the last group to flag for the shower + The index of the last group to flag for the shower """ num_grps_masked_after += 1 From e4017c36da212b0b9731f22baf90166bf489cfe5 Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:45:05 -0500 Subject: [PATCH 30/40] Update tests/test_jump.py --- tests/test_jump.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 5068ee95..700807a8 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -356,12 +356,12 @@ def test_calc_num_slices(): def test_find_last_grp(): - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=0) == 6) - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=2) == 7) - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=3) == 7) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=1) == 6) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=0) == 6) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=2) == 6) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=0) == 6) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=1) == 7) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=2) == 8) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=0) == 6) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=2) == 7) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=3) == 7) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=1) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=0) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=2) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=0) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=1) == 7) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=2) == 8) From a3cbe27a4eabcdc9b1f198aa157ed4fa99f20fd9 Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:50:56 -0500 Subject: [PATCH 31/40] Update tests/test_jump.py --- tests/test_jump.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 700807a8..5068ee95 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -356,12 +356,12 @@ def test_calc_num_slices(): def test_find_last_grp(): - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=0) == 6) - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=2) == 7) - assert (find_last_grp(grp=5, ngrps=7, num_grps_masked_after=3) == 7) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=1) == 6) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=0) == 6) - assert (find_last_grp(grp=5, ngrps=6, num_grps_masked_after=2) == 6) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=0) == 6) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=1) == 7) - assert (find_last_grp(grp=5, ngrps=8, num_grps_masked_after=2) == 8) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=2) == 7) + assert (find_last_grp(grp=5, ngrps=7, num_grps_masked=3) == 7) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=1) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=6, num_grps_masked=2) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=0) == 6) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=1) == 7) + assert (find_last_grp(grp=5, ngrps=8, num_grps_masked=2) == 8) From e69926abc22452a3f573deee0e18be45996099c1 Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:51:07 -0500 Subject: [PATCH 32/40] Update src/stcal/jump/jump.py --- src/stcal/jump/jump.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 5cde6522..647d0867 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -702,8 +702,8 @@ def find_last_grp(grp, ngrps, num_grps_masked_after): The index of the last group to flag for the shower """ - num_grps_masked_after += 1 - last_grp = min(grp + num_grps_masked_after, ngrps) + num_grps_masked += 1 + last_grp = min(grp + num_grps_masked, ngrps) return last_grp def find_circles(dqplane, bitmask, min_area): From 07186c4a3e9e41b1b4db47ef080ab622c980c025 Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:51:17 -0500 Subject: [PATCH 33/40] Update src/stcal/jump/jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 647d0867..06c60f93 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -693,7 +693,7 @@ def find_last_grp(grp, ngrps, num_grps_masked_after): ngrps : int The number of groups in the integration - num_grps_masked_after : int + num_grps_masked : int The requested number of groups to be flagged after the shower Returns From 304a4b0b0ff7040c27ebe40faeb29659dcea5475 Mon Sep 17 00:00:00 2001 From: Howard Bushouse Date: Wed, 31 Jan 2024 11:51:27 -0500 Subject: [PATCH 34/40] Update src/stcal/jump/jump.py --- src/stcal/jump/jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 06c60f93..945cf69a 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -683,7 +683,7 @@ def extend_ellipses( diff_cube = out_gdq_cube - gdq_cube return out_gdq_cube, num_ellipses -def find_last_grp(grp, ngrps, num_grps_masked_after): +def find_last_grp(grp, ngrps, num_grps_masked): """ Parameters ---------- From d4b7e0f2905799b023dcbc7d9850c938f65b7786 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 15:13:23 -0500 Subject: [PATCH 35/40] Update test_jump.py --- tests/test_jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/test_jump.py b/tests/test_jump.py index 5068ee95..67a55820 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -127,6 +127,7 @@ def test_flag_large_events_withsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] + print(cube[0, 2, :, :]) flag_large_events( cube, DQFLAGS["JUMP_DET"], @@ -142,6 +143,7 @@ def test_flag_large_events_withsnowball(): assert cube[0, 1, 2, 2] == 0 assert cube[0, 1, 3, 5] == 0 assert cube[0, 2, 0, 0] == 0 + print(cube[0, 2, :, :]) assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended assert cube[0, 2, 3, 6] == DQFLAGS["JUMP_DET"] From 89d7ef56d2a42c6ab8e447be19238e45b55ca45c Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 16:16:17 -0500 Subject: [PATCH 36/40] fixes --- src/stcal/jump/jump.py | 8 ++++---- tests/test_jump.py | 7 +++---- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 945cf69a..d03dfa66 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -281,7 +281,7 @@ def detect_jumps( ) # This is the flag that controls the flagging of snowballs. if expand_large_events: - total_snowballs = flag_large_events( + gdq, total_snowballs = flag_large_events( gdq, jump_flag, sat_flag, @@ -431,7 +431,7 @@ def detect_jumps( k += 1 # This is the flag that controls the flagging of snowballs. if expand_large_events: - total_snowballs = flag_large_events( + gdq, total_snowballs = flag_large_events( gdq, jump_flag, sat_flag, @@ -581,7 +581,7 @@ def flag_large_events( expansion=expand_factor, max_extended_radius=max_extended_radius, ) - return total_snowballs + return gdq, total_snowballs def extend_saturation( @@ -957,12 +957,12 @@ def find_faint_extended( if len(ellipses) > 0: # add all the showers for this integration to the list all_ellipses.append([intg, grp, ellipses]) + total_showers = 0 if all_ellipses: # Now we actually do the flagging of the pixels inside showers. # This is deferred until all showers are detected. because the showers # can flag future groups and would confuse the detection algorithm if # we worked on groups that already had some flagged showers. - total_showers = 0 for showers in all_ellipses: intg = showers[0] grp = showers[1] diff --git a/tests/test_jump.py b/tests/test_jump.py index 67a55820..b92479ee 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -127,8 +127,7 @@ def test_flag_large_events_withsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - print(cube[0, 2, :, :]) - flag_large_events( + cube, total_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], @@ -143,7 +142,6 @@ def test_flag_large_events_withsnowball(): assert cube[0, 1, 2, 2] == 0 assert cube[0, 1, 3, 5] == 0 assert cube[0, 2, 0, 0] == 0 - print(cube[0, 2, :, :]) assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended assert cube[0, 2, 3, 6] == DQFLAGS["JUMP_DET"] @@ -162,7 +160,7 @@ def test_flag_large_events_groupedsnowball(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - flag_large_events( + cube, total_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], @@ -174,6 +172,7 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) + print(cube[0, 2, :, :]) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 assert cube[0, 2, 0, 0] == 0 From f50aaadd4eb9f083b991921c06fe178cbfdafd65 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 16:31:09 -0500 Subject: [PATCH 37/40] Update test_jump.py --- tests/test_jump.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index b92479ee..7281ab04 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -172,7 +172,8 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) - print(cube[0, 2, :, :]) + print("") + print("final corner", cube[0, 2, 0, 0]) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 assert cube[0, 2, 0, 0] == 0 @@ -192,7 +193,7 @@ def test_flag_large_events_withsnowball_noextension(): cube[0, 2, 5, 1:6] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 1] = DQFLAGS["JUMP_DET"] cube[0, 2, 1:6, 5] = DQFLAGS["JUMP_DET"] - flag_large_events( + cube, num_snowballs = flag_large_events( cube, DQFLAGS["JUMP_DET"], DQFLAGS["SATURATED"], From 26274dfd373b657a1937cb00cf98f0c5c64ed352 Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 16:37:13 -0500 Subject: [PATCH 38/40] Update jump.py --- src/stcal/jump/jump.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index d03dfa66..280ddb42 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -681,6 +681,8 @@ def extend_ellipses( jump_ellipse[saty, satx] = 0 out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) diff_cube = out_gdq_cube - gdq_cube + print("") + print("extend ellipse corner", out_gdq_cube[0, 2, 0, 0]) return out_gdq_cube, num_ellipses def find_last_grp(grp, ngrps, num_grps_masked): From 9ec7676fed991b51248f295b7ca38750b57aa92a Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 16:46:10 -0500 Subject: [PATCH 39/40] Update test_jump.py --- tests/test_jump.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_jump.py b/tests/test_jump.py index 7281ab04..7de81085 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -176,7 +176,7 @@ def test_flag_large_events_groupedsnowball(): print("final corner", cube[0, 2, 0, 0]) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 - assert cube[0, 2, 0, 0] == 0 + assert cube[0, 2, 0, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended From 037966e8197f08768744bdd2e03aa7c8fa0641ae Mon Sep 17 00:00:00 2001 From: Mike Regan Date: Wed, 31 Jan 2024 17:05:59 -0500 Subject: [PATCH 40/40] fix num_grps_masked for snowballs --- src/stcal/jump/jump.py | 3 +-- tests/test_jump.py | 4 +--- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/src/stcal/jump/jump.py b/src/stcal/jump/jump.py index 280ddb42..d197ce87 100644 --- a/src/stcal/jump/jump.py +++ b/src/stcal/jump/jump.py @@ -579,6 +579,7 @@ def flag_large_events( sat_flag, jump_flag, expansion=expand_factor, + num_grps_masked=0, max_extended_radius=max_extended_radius, ) return gdq, total_snowballs @@ -681,8 +682,6 @@ def extend_ellipses( jump_ellipse[saty, satx] = 0 out_gdq_cube[intg, flg_grp, :, :] = np.bitwise_or(gdq_cube[intg, flg_grp, :, :], jump_ellipse) diff_cube = out_gdq_cube - gdq_cube - print("") - print("extend ellipse corner", out_gdq_cube[0, 2, 0, 0]) return out_gdq_cube, num_ellipses def find_last_grp(grp, ngrps, num_grps_masked): diff --git a/tests/test_jump.py b/tests/test_jump.py index 7de81085..995a5367 100644 --- a/tests/test_jump.py +++ b/tests/test_jump.py @@ -172,11 +172,9 @@ def test_flag_large_events_groupedsnowball(): min_sat_radius_extend=0.5, sat_expand=1.1, ) - print("") - print("final corner", cube[0, 2, 0, 0]) # assert cube[0, 1, 2, 2] == 0 # assert cube[0, 1, 3, 5] == 0 - assert cube[0, 2, 0, 0] == DQFLAGS["JUMP_DET"] # Jump was extended + assert cube[0, 2, 0, 0] == 0 assert cube[0, 2, 1, 0] == DQFLAGS["JUMP_DET"] # Jump was extended assert cube[0, 2, 2, 2] == DQFLAGS["SATURATED"] # Saturation was extended