diff --git a/CHANGES.rst b/CHANGES.rst index bf59741484..00ffe2c961 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,12 @@ 0.18.3 (unreleased) =================== +ramp_fitting +------------ + +- Fix a bug in estimating the max number of segments that will be needed + to fit any pixel [#5653] + white_light ----------- @@ -8,7 +14,6 @@ white_light for MJDs [#5566] - 0.18.2 (2021-01-19) =================== diff --git a/jwst/ramp_fitting/ramp_fit.py b/jwst/ramp_fitting/ramp_fit.py index ce89e3f13e..8515d549d9 100755 --- a/jwst/ramp_fitting/ramp_fit.py +++ b/jwst/ramp_fitting/ramp_fit.py @@ -30,6 +30,9 @@ log = logging.getLogger(__name__) log.setLevel(logging.DEBUG) +DO_NOT_USE = dqflags.group['DO_NOT_USE'] +JUMP_DET = dqflags.group['JUMP_DET'] + BUFSIZE = 1024 * 300000 # 300Mb cache size for data section @@ -201,6 +204,7 @@ def ols_ramp_fit_multi(input_model, buffsize, save_opt, readnoise_2d, gain_2d, # Call ramp fitting for the single processor (1 data slice) case if number_slices == 1: max_segments, max_CRs = calc_num_seg(input_model.groupdq, number_of_integrations) + log.debug(f"Max segments={max_segments}") int_model, opt_model, out_model = create_output_models(input_model, number_of_integrations, save_opt, total_cols, total_rows, max_segments, max_CRs) @@ -2940,12 +2944,11 @@ def calc_num_seg(gdq, n_int): """ max_cr = 0 # max number of CRS for all integrations - # For all 2d pixels, get max number of CRs along their ramps - # to use as a surrogate for the number of semiramps along the ramps - for nint in range(n_int): # loop over integrations - gdq_cr = np.bitwise_and( gdq[nint, :, :, :], dqflags.group['JUMP_DET']) - temp_max_cr = int((gdq_cr.sum(axis=0)).max()/dqflags.group['JUMP_DET']) - max_cr = max( max_cr, temp_max_cr ) + # For all 2d pixels, get max number of CRs or DO_NOT_USE flags along their + # ramps, to use as a surrogate for the number of segments along the ramps + # Note that we only care about flags that are NOT in the first or last groups, + # because exclusion of a first or last group won't result in an additional segment. + max_cr = np.count_nonzero(np.bitwise_and(gdq[:, 1:-1], JUMP_DET | DO_NOT_USE), axis=1).max() # Do not want to return a value > the number of groups, which can occur if # this is a MIRI dataset in which the first or last group was flagged as diff --git a/jwst/ramp_fitting/tests/test_ramp_fit.py b/jwst/ramp_fitting/tests/test_ramp_fit.py index 659cbbe684..f11c5d3280 100644 --- a/jwst/ramp_fitting/tests/test_ramp_fit.py +++ b/jwst/ramp_fitting/tests/test_ramp_fit.py @@ -2,10 +2,15 @@ import numpy as np from jwst.ramp_fitting.ramp_fit import ramp_fit +from jwst.ramp_fitting.ramp_fit import calc_num_seg from jwst.datamodels import dqflags from jwst.datamodels import RampModel from jwst.datamodels import GainModel, ReadnoiseModel +DO_NOT_USE = dqflags.group['DO_NOT_USE'] +JUMP_DET = dqflags.group['JUMP_DET'] +SATURATED = dqflags.group['SATURATED'] + # single group intergrations fail in the GLS fitting # so, keep the two method test separate and mark GLS test as @@ -43,6 +48,43 @@ def test_drop_frames1_not_set(): np.testing.assert_allclose(slopes[0].data[50, 50],10.0, 1e-6) +def test_mixed_crs_and_donotuse(): + + gdq = np.zeros((3,10,3,3), dtype=np.uint32) + + # pix with only first and last group flagged DO_NOT_USE; + # results in 1 segment (flags at ends of ramp do not break the ramp) + gdq[0,0,0,0] = DO_NOT_USE + gdq[0,-1,0,0] = DO_NOT_USE + + # pix with first and last group flagged DO_NOT_USE and 1 CR in middle + # results in 2 segments + gdq[0,0,1,1] = DO_NOT_USE + gdq[0,-1,1,1] = DO_NOT_USE + gdq[0,3,1,1] = JUMP_DET + + # max segments should be 2 + max_seg, max_cr = calc_num_seg(gdq, 3) + assert(max_seg == 2) + + # pix with only 1 middle group flagged DO_NOT_USE; + # results in 2 segments + gdq[1,2,0,0] = DO_NOT_USE + + # pix with middle group flagged as CR and DO_NOT_USE; + # results in 2 segments + gdq[2,2,0,0] = DO_NOT_USE+JUMP_DET + + # pix with DO_NOT_USE and CR in different middle groups; + # results in 3 segments + gdq[2,2,1,1] = DO_NOT_USE + gdq[2,4,1,1] = JUMP_DET + + # max segments should now be 3 + max_seg, max_cr = calc_num_seg(gdq, 3) + assert(max_seg == 3) + + @pytest.mark.skip(reason="GLS code does not [yet] handle single group integrations.") def test_one_group_small_buffer_fit_gls(): model1, gdq, rnModel, pixdq, err, gain = setup_inputs(ngroups=1,gain=1,readnoise=10) @@ -281,7 +323,7 @@ def test_photon_noise_only_bad_last_frame(self, method): model1.data[0, 2, 50, 50] = 25.0 model1.data[0, 3, 50, 50] = 33.0 model1.data[0, 4, 50, 50] = 60.0 - model1.groupdq[0,4,:,:] = dqflags.group['DO_NOT_USE'] + model1.groupdq[0,4,:,:] = DO_NOT_USE slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') cds_slope = (model1.data[0,3,50,50] - model1.data[0,0,50,50])/ 3.0 np.testing.assert_allclose(slopes[0].data[50, 50], cds_slope, 1e-2) @@ -293,7 +335,7 @@ def test_photon_noise_only_bad_last_frame_two_groups(self, method): model1, gdq, rnModel, pixdq, err, gain = setup_inputs(ngroups=2,gain=1000,readnoise=1) model1.data[0, 0, 50, 50] = 10.0 model1.data[0, 1, 50, 50] = 15.0 - model1.groupdq[0,1,:,:] = dqflags.group['DO_NOT_USE'] + model1.groupdq[0,1,:,:] = DO_NOT_USE slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') cds_slope = (model1.data[0,1,50,50] - model1.data[0,0,50,50])/ 1.0 np.testing.assert_allclose(slopes[0].data[50, 50], cds_slope, 1e-6) @@ -322,17 +364,17 @@ def test_two_groups_fit(self, method): model1.data[0, 1, 50, 52] = 600.0 model1.meta.exposure.drop_frames1 = 0 #2nd group is saturated - model1.groupdq[0,1,50,51]=dqflags.group['SATURATED'] + model1.groupdq[0,1,50,51]=SATURATED #1st group is saturated - model1.groupdq[0,0,50,52]=dqflags.group['SATURATED'] - model1.groupdq[0,1,50,52]=dqflags.group['SATURATED'] #should not be set this way + model1.groupdq[0,0,50,52]=SATURATED + model1.groupdq[0,1,50,52]=SATURATED # should not be set this way slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') cds_slope = (model1.data[0,1,50,50] - model1.data[0,0,50,50]) np.testing.assert_allclose(slopes[0].data[50, 50], cds_slope, 1e-6) #expect SATURATED - assert slopes[0].dq[50, 51] == dqflags.pixel['SATURATED'] + assert slopes[0].dq[50, 51] == SATURATED #expect SATURATED and DO_NOT_USE, because 1st group is Saturated - assert slopes[0].dq[50, 52] == dqflags.pixel['SATURATED'] + dqflags.pixel['DO_NOT_USE'] + assert slopes[0].dq[50, 52] == SATURATED + DO_NOT_USE def test_four_groups_oneCR_orphangroupatend_fit(self, method): model1, gdq, rnModel, pixdq, err, gain = setup_inputs(ngroups=4,gain=1,readnoise=10) @@ -340,7 +382,7 @@ def test_four_groups_oneCR_orphangroupatend_fit(self, method): model1.data[0, 1, 50, 50] = 15.0 model1.data[0, 2, 50, 50] = 20.0 model1.data[0, 3, 50, 50] = 145.0 - model1.groupdq[0,3,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,3,50,50]=JUMP_DET slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') cds_slope = (model1.data[0,1,50,50] - model1.data[0,0,50,50]) @@ -353,8 +395,8 @@ def test_four_groups_two_CRs_at_end(self, method): model1.data[0, 1, 50, 50] = 15.0 model1.data[0, 2, 50, 50] = 25.0 model1.data[0, 3, 50, 50] = 145.0 - model1.groupdq[0,2,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,3,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,2,50,50]=JUMP_DET + model1.groupdq[0,3,50,50]=JUMP_DET slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') cds_slope = (model1.data[0,1,50,50] - model1.data[0,0,50,50]) np.testing.assert_allclose(slopes[0].data[50, 50], cds_slope, 1e-6) @@ -366,10 +408,10 @@ def test_four_groups_four_CRs(self, method): model1.data[0, 1, 50, 50] = 15.0 model1.data[0, 2, 50, 50] = 25.0 model1.data[0, 3, 50, 50] = 145.0 - model1.groupdq[0,0,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,1,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,2,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,3,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,0,50,50]=JUMP_DET + model1.groupdq[0,1,50,50]=JUMP_DET + model1.groupdq[0,2,50,50]=JUMP_DET + model1.groupdq[0,3,50,50]=JUMP_DET slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') np.testing.assert_allclose(slopes[0].data[50, 50], 0,1e-6) @@ -380,9 +422,9 @@ def test_four_groups_three_CRs_at_end(self, method): model1.data[0, 1, 50, 50] = 15.0 model1.data[0, 2, 50, 50] = 25.0 model1.data[0, 3, 50, 50] = 145.0 - model1.groupdq[0,1,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,2,50,50]=dqflags.group['JUMP_DET'] - model1.groupdq[0,3,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,1,50,50]=JUMP_DET + model1.groupdq[0,2,50,50]=JUMP_DET + model1.groupdq[0,3,50,50]=JUMP_DET slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') expected_slope=10.0 np.testing.assert_allclose(slopes[0].data[50, 50],expected_slope, 1e-6) @@ -394,7 +436,7 @@ def test_four_groups_CR_causes_orphan_1st_group(self, method): model1.data[0, 1, 50, 50] = 125.0 model1.data[0, 2, 50, 50] = 145.0 model1.data[0, 3, 50, 50] = 165.0 - model1.groupdq[0,1,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,1,50,50]=JUMP_DET slopes = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') expected_slope=20.0 np.testing.assert_allclose(slopes[0].data[50, 50],expected_slope, 1e-6) @@ -474,7 +516,7 @@ def test_oneCR_10_groups_combination(self, method): model1.data[0, 7, 50, 50] = 160.0 model1.data[0, 8, 50, 50] = 170.0 model1.data[0, 9, 50, 50] = 180.0 - model1.groupdq[0,5,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,5,50,50]=JUMP_DET slopes, int_model, opt_model, gls_opt_model = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') segment_groups = 5 @@ -507,7 +549,7 @@ def test_oneCR_10_groups_combination_noisy2ndSegment(self, method): model1.data[0, 7, 50, 50] = 160.0 model1.data[0, 8, 50, 50] = 168.0 model1.data[0, 9, 50, 50] = 180.0 - model1.groupdq[0,5,50,50]=dqflags.group['JUMP_DET'] + model1.groupdq[0,5,50,50]=JUMP_DET slopes, int_model, opt_model, gls_opt_model= ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') avg_slope = (opt_model.slope[0,0,50,50] + opt_model.slope[0,1,50,50])/2.0 @@ -528,14 +570,14 @@ def test_twenty_groups_two_segments(): # b) ramp having 1 CR at group 15; 2 segments model1.data[0,:,0,1] = np.arange(ngroups)*10. + 50. - gdq[0,15,0,1] = dqflags.group['JUMP_DET'] + gdq[0,15,0,1] = JUMP_DET model1.data[0,15:,0,1] += 1000. # c) ramp having 1 CR at group 2; SAT starting in group 15 model1.data[0,:,0,2] = np.arange(ngroups)*10. + 70. - gdq[0,2,0,2] = dqflags.group['JUMP_DET'] + gdq[0,2,0,2] = JUMP_DET model1.data[0,2:,0,2] += 2000. - gdq[0,15:,0,2] = dqflags.group['SATURATED'] + gdq[0,15:,0,2] = SATURATED model1.data[0,15:,0,2] = 25000. new_mod, int_model, opt_model, gls_opt_model = ramp_fit(model1, 1024*30000., @@ -559,7 +601,7 @@ def test_miri_all_sat(): model1, gdq, rnModel, pixdq, err, gain = setup_small_cube(ngroups, nints, nrows, ncols, deltatime) - model1.groupdq[:, :, :, :] = dqflags.group['SATURATED'] + model1.groupdq[:, :, :, :] = SATURATED new_mod, int_model, opt_model, gls_opt_model = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') @@ -612,11 +654,11 @@ def test_miri_first_last(): 190., 200., 210., -400.], dtype=np.float32) # For all pixels, set gdq for 0th and final groups to DO_NOT_USE - model1.groupdq[:,0,:,:] = dqflags.group['DO_NOT_USE'] - model1.groupdq[:,-1,:,:] = dqflags.group['DO_NOT_USE'] + model1.groupdq[:,0,:,:] = DO_NOT_USE + model1.groupdq[:,-1,:,:] = DO_NOT_USE # Put CR in 1st (0-based) group - model1.groupdq[0,1,1,1] = dqflags.group['JUMP_DET'] + model1.groupdq[0,1,1,1] = JUMP_DET new_mod, int_model, opt_model, gls_opt_model = ramp_fit(model1, 1024*30000., True, rnModel, gain, 'OLS', 'optimal', 'none') @@ -677,13 +719,14 @@ def setup_small_cube(ngroups=10, nints=1, nrows=2, ncols=2, deltatime=10., def setup_inputs(ngroups=10, readnoise=10, nints=1, nrows=103, ncols=102, nframes=1, grouptime=1.0,gain=1, deltatime=1): + data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + pixdq = np.zeros(shape=(nrows, ncols), dtype= np.uint32) + gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8) gain = np.ones(shape=(nrows, ncols), dtype=np.float64) * gain - err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float64) - data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint32) - pixdq = np.zeros(shape=(nrows, ncols), dtype= np.float64) - read_noise = np.full((nrows, ncols), readnoise, dtype=np.float64) - gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.int32) + read_noise = np.full((nrows, ncols), readnoise, dtype=np.float32) int_times = np.zeros((nints,)) + model1 = RampModel(data=data, err=err, pixeldq=pixdq, groupdq=gdq, int_times=int_times) model1.meta.instrument.name='MIRI' model1.meta.instrument.detector='MIRIMAGE' @@ -722,13 +765,14 @@ def setup_subarray_inputs(ngroups=10, readnoise=10, nints=1, subxsize=1024, subysize=1032, nframes=1, grouptime=1.0,gain=1, deltatime=1): - times = np.array(list(range(ngroups)),dtype=np.float64) * deltatime + data = np.zeros(shape=(nints, ngroups, subysize, subxsize), dtype=np.float32) + err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) + pixdq = np.zeros(shape=(subysize, subxsize), dtype= np.uint32) + gdq = np.zeros(shape=(nints, ngroups, subysize, subxsize), dtype=np.uint8) gain = np.ones(shape=(nrows, ncols), dtype=np.float64) * gain - err = np.ones(shape=(nints, ngroups, nrows, ncols), dtype=np.float64) - data = np.zeros(shape=(nints, ngroups, subysize, subxsize), dtype=np.float64) - pixdq = np.zeros(shape=(subysize, subxsize), dtype= np.float64) - read_noise = np.full((nrows, ncols), readnoise, dtype=np.float64) - gdq = np.zeros(shape=(nints, ngroups, subysize, subxsize), dtype=np.int32) + read_noise = np.full((nrows, ncols), readnoise, dtype=np.float32) + times = np.array(list(range(ngroups)),dtype=np.float64) * deltatime + model1 = RampModel(data=data, err=err, pixeldq=pixdq, groupdq=gdq, times=times) model1.meta.instrument.name='MIRI' model1.meta.instrument.detector='MIRIMAGE' diff --git a/jwst/ramp_fitting/tests/test_ramp_fit_cases.py b/jwst/ramp_fitting/tests/test_ramp_fit_cases.py index 5757523e8b..1629a15fd7 100644 --- a/jwst/ramp_fitting/tests/test_ramp_fit_cases.py +++ b/jwst/ramp_fitting/tests/test_ramp_fit_cases.py @@ -613,7 +613,7 @@ def create_mod_arrays(ngroups, nints, nrows, ncols, deltatime, gain, readnoise): gain = np.ones(shape=(nrows, ncols), dtype=np.float32) * gain err = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) data = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.float32) - pixdq = np.zeros(shape=(nrows, ncols), dtype=np.int32) + pixdq = np.zeros(shape=(nrows, ncols), dtype=np.uint32) read_noise = np.full((nrows, ncols), readnoise, dtype=np.float32) gdq = np.zeros(shape=(nints, ngroups, nrows, ncols), dtype=np.uint8)