diff --git a/CHANGES.rst b/CHANGES.rst index f45add44..b4faab90 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -17,6 +17,13 @@ ramp_fitting Bug Fixes --------- +ramp_fitting +~~~~~~~~~~~~ + +- Changed the data type of three variables that are used in measuring + the jump free segments of integrations. The variables were uint8 and + they would yield wrong results for integrations with more than 256 + groups. [#251] Other ----- diff --git a/src/stcal/ramp_fitting/utils.py b/src/stcal/ramp_fitting/utils.py index 73f4b5e4..89515aad 100644 --- a/src/stcal/ramp_fitting/utils.py +++ b/src/stcal/ramp_fitting/utils.py @@ -520,10 +520,10 @@ def calc_slope_vars(ramp_data, rn_sect, gain_sect, gdq_sect, group_time, max_seg gdq_2d_nan[np.bitwise_and(gdq_2d, ramp_data.flags_saturated).astype(bool)] = np.nan # Get lengths of semiramps for all pix [number_of_semiramps, number_of_pix] - segs = np.zeros_like(gdq_2d) + segs = np.zeros_like(gdq_2d).astype(np.uint16) # Counter of semiramp for each pixel - sr_index = np.zeros(npix, dtype=np.uint8) + sr_index = np.zeros(npix, dtype=np.uint16) pix_not_done = np.ones(npix, dtype=bool) # initialize to True i_read = 0 @@ -558,7 +558,7 @@ def calc_slope_vars(ramp_data, rn_sect, gain_sect, gdq_sect, group_time, max_seg i_read += 1 - segs = segs.astype(np.uint8) + segs = segs.astype(np.uint16) segs_beg = segs[:max_seg, :] # the leading nonzero lengths # Create reshaped version [ segs, y, x ] to simplify computation diff --git a/tests/test_ramp_fitting.py b/tests/test_ramp_fitting.py index a6edbe33..93674e09 100644 --- a/tests/test_ramp_fitting.py +++ b/tests/test_ramp_fitting.py @@ -28,6 +28,23 @@ # ----------------------------------------------------------------------------- # Test Suite +def test_long_integration(): + nints, nrows, ncols = 1, 1, 1 + rnoise_val, gain_val = 0.1, 40.0 + nframes, gtime, ftime = 1, 3, 3 + tm = (nframes, gtime, ftime) + num_grps1 = 301 + num_grps2 = 20 + ramp_data, rnoise_array, gain_array = create_test_2seg_obs(rnoise_val, nints, num_grps1, num_grps2, ncols, + nrows, tm, rate=0, Poisson=True, grptime=gtime, + gain=gain_val, bias=0) + ramp_data.data[0, 291:, 0, 0] = 320 * 3 + # Run ramp fit on RampData + buffsize, save_opt, algo, wt, ncores = 512, True, "OLS", "optimal", "none" + slopes, cube, optional, gls_dummy = ramp_fit_data( + ramp_data, buffsize, save_opt, rnoise_array, gain_array, algo, wt, ncores, dqflags) + np.testing.assert_almost_equal(slopes[0], .65, 2) + def base_neg_med_rates_single_integration(): """ @@ -1506,6 +1523,49 @@ def setup_inputs(dims, var, tm): return ramp_data, rnoise, gain +def create_test_2seg_obs(readnoise, num_ints, num_grps1, num_grps2, ncols, + nrows, tm, rate=0, Poisson=True, grptime=2.77, + gain=4.0, bias=3000, sat_group=0, sat_value=100000): + nframes, gtime, dtime = tm + rng = np.random.default_rng() + outcube1a = np.zeros(shape=(num_ints, num_grps1 + num_grps2, ncols, nrows), dtype=np.float32) + outcube1 = np.random.normal(loc=0.0, scale=readnoise / np.sqrt(2), + size=(num_ints, num_grps1 + num_grps2 + 1, ncols, ncols)) + if rate > 0: + pvalues = grptime * rate + (rng.poisson(lam=gain * rate * grptime, + size=(num_ints, num_grps1 + num_grps2, ncols, + nrows)) - gain * rate * grptime) / gain + for intg in range(num_ints): + outcube1a[intg, 0, :, :] = outcube1[intg, 0, :, :] + for grp in range(1, num_grps1 + num_grps2): + outcube1a[intg, grp, :, :] = outcube1[intg, grp, :, :] + np.sum(pvalues[intg, 0:grp, :, :], axis=0) + outcube1f = outcube1a + else: + outcube1f = outcube1 + outdata = outcube1f + bias + # print("cube mean values", np.mean(outdata[0,:,:,:], axis=(2, 3))) + outgdq = np.zeros_like(outdata, dtype=np.uint8) + outgdq[:, 0, :, :] = 1 + outgdq[:, -1, :, :] = 1 + if num_grps2 > 0: + outgdq[:, num_grps1, :, :] = 4 + if sat_group > 0: + outgdq[:, sat_group:, :, :] = 2 + outdata[:, sat_group:, :, :] = sat_value + pixdq = np.zeros(shape=(ncols, nrows), dtype=np.int32) + err = np.ones(shape=(num_ints, num_grps1 + num_grps2 + 1, nrows, ncols), dtype=np.float32) + ramp_data = RampData() + dark_current = np.zeros((nrows, ncols)) + ramp_data.set_arrays( + data=outdata, err=err, groupdq=outgdq, pixeldq=pixdq, average_dark_current=dark_current) + ramp_data.set_meta( + name="MIRI", frame_time=dtime, group_time=gtime, groupgap=0, + nframes=nframes, drop_frames1=None) + ramp_data.set_dqflags(dqflags) + readnoise_array = np.ones_like(pixdq) * readnoise + gain_array = np.ones_like(pixdq) * gain + return ramp_data, readnoise_array, gain_array + # -----------------------------------------------------------------------------