JP-3583: Long integrations truncation in ramp fitting (#251)

* Update utils.py

* Update test_ramp_fitting.py

* Update test_ramp_fitting.py

* update test

* Update CHANGES.rst

* Update CHANGES.rst

---------

Co-authored-by: Howard Bushouse <bushouse@stsci.edu>

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
 -----
 

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

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
+
 
 # -----------------------------------------------------------------------------