Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

JP-1860: Fix ramp_fit counting of segments #5653

Merged
merged 6 commits into from
Jan 22, 2021
Merged

JP-1860: Fix ramp_fit counting of segments #5653

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
136c9d2
JP-1860: Fix ramp_fit counting of segments
hbushouse Jan 21, 2021
39cea44
Update change log
hbushouse Jan 21, 2021
cdbd200
Fix logic for counting flagged groups
hbushouse Jan 21, 2021
2e84039
add unit test for calc_num_seg function
hbushouse Jan 22, 2021
e79bc4e
a little cleanup of tests
hbushouse Jan 22, 2021
31ac154
Merge branch 'master' into jp1860
hbushouse Jan 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
7 changes: 6 additions & 1 deletion CHANGES.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,14 +1,19 @@
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
-----------

- Fixed error causing multi-segment data to reject int_times
for MJDs [#5566]



0.18.2 (2021-01-19)
===================

Expand Down
15 changes: 9 additions & 6 deletions jwst/ramp_fitting/ramp_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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()
Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

note that max_cr is bumped up by 1 later in line 2956, to account for the fact that n flags implies n+1 segments.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor comment: Is there a reason not to bump now? Or is the later use a conceptually different value, segments, as opposed to "maximum rejects".

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No really good reason I suppose, although they are a bit different conceptually. This line simply computes the max number "interrupts" in any ramp, and then that number of interrupts turns into n+1 segments for fitting. Besides, the line is so long already, I didn't want to add yet another operation on the end of it!


# 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
Expand Down
120 changes: 82 additions & 38 deletions jwst/ramp_fitting/tests/test_ramp_fit.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -322,25 +364,25 @@ 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)
model1.data[0, 0, 50, 50] = 10.0
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])

Expand All @@ -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)
Expand All @@ -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)

Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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.,
Expand All @@ -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')
Expand Down Expand Up @@ -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')
Expand Down Expand Up @@ -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'
Expand Down Expand Up @@ -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'
Expand Down
2 changes: 1 addition & 1 deletion jwst/ramp_fitting/tests/test_ramp_fit_cases.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)

Expand Down