Use new correction factor formulation

src/stcal/ramp_fitting/ols_cas22/_pixel.pxd

@@ -14,7 +14,7 @@ cdef class Pixel:
     cdef float[:, :] local_slope_vals(Pixel self)
     cdef RampFit fit_ramp(Pixel self, RampIndex ramp)
 
-    cdef float correction(Pixel self, RampIndex ramp, int index, int diff)
+    cdef float correction(Pixel self, RampIndex ramp, float slope)
     cdef float stat(Pixel self, float slope, RampIndex ramp, int index, int diff)
     cdef float[:] stats(Pixel self, float slope, RampIndex ramp)
     cdef RampFits fit_ramps(Pixel self, stack[RampIndex] ramps)

src/stcal/ramp_fitting/ols_cas22/_pixel.pyx

@@ -220,43 +220,36 @@ cdef class Pixel:
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
-    cdef inline float correction(Pixel self, RampIndex ramp, int index, int diff):
+    cdef inline float correction(Pixel self, RampIndex ramp, float slope):
         """
         Compute the correction factor for the variance used by a statistic
+
+            - slope / (t_bar[end] - t_bar[start])
         
         Parameters
         ----------
         ramp : RampIndex
             Struct for start and end indices resultants for the ramp
-        index : int
-            The main index for the resultant to compute the statistic for
-        diff : int
-            The offset to use for the delta and sigma values, this should be
-            a value from the Diff enum.
+        slope : float
+            The computed slope for the ramp
         """
-        cdef float comp = (self.fixed.t_bar_diffs[diff, index] /
-                           (self.fixed.data.t_bar[ramp.end] - self.fixed.data.t_bar[ramp.start]))
 
-        if diff == Diff.single:
-            return (1 - comp)**2
-        elif diff == Diff.double:
-            return (1 - 0.75 * comp)**2
-        else:
-            raise ValueError("diff must be Diff.single or Diff.double")
+        cdef float diff = (self.fixed.data.t_bar[ramp.end] - self.fixed.data.t_bar[ramp.start])
+
+        return - slope / diff
 
     @cython.boundscheck(False)
     @cython.wraparound(False)
     @cython.cdivision(True)
     cdef inline float stat(Pixel self, float slope, RampIndex ramp, int index, int diff):
         """
         Compute a single set of fit statistics
-            delta / sqrt(var)
+            (delta / sqrt(var)) + correction
         where
             delta = ((R[j] - R[i]) / (t_bar[j] - t_bar[i]) - slope)
                     * (t_bar[j] - t_bar[i])
             var   = sigma * (1/N[j] + 1/N[i]) 
                     + slope * (tau[j] + tau[i] - min(t_bar[j], t_bar[i]))
-                    * correction(diff)
         
         Parameters
         ----------
@@ -274,13 +267,13 @@ cdef class Pixel:
         -------
             Create a single instance of the stastic for the given parameters
         """
-        cdef float delta = ((self.local_slopes[diff, index] - slope) *
-                            fabs(self.fixed.t_bar_diffs[diff, index]))
-        cdef float var = (self.var_read_noise[diff, index] +
-                          slope * self.fixed.var_slope_coeffs[diff, index] *
-                                  self.correction(ramp, index, diff))
+        cdef float delta = (self.local_slopes[diff, index] - slope)
+        cdef float var = ((self.var_read_noise[diff, index] +
+                           slope * self.fixed.var_slope_coeffs[diff, index])
+                          / self.fixed.t_bar_diff_sqrs[diff, index]) 
+        cdef float correct = self.correction(ramp, slope)
 
-        return delta / sqrt(var)
+        return (delta / sqrt(var)) + correct
 
 
     @cython.boundscheck(False)

tests/test_jump_cas22.py

@@ -396,6 +396,7 @@ def test_find_jumps(jump_data):
 
     # Check that all the jumps have been located per the algorithm's constraints
     for index, (fit, jump) in enumerate(zip(fits, jumps)):
+        print(f"{index=}, {fit['jumps']=}, {jump=}")
         # sanity check that only one jump should have been added
         assert np.where(jump)[0].shape == (1,)
         if index == 0:
@@ -419,25 +420,25 @@ def test_find_jumps(jump_data):
             # Test that all the jumps recorded are +/- 1 of the real jump
             #    This is due to the need to exclude two resultants
             for jump_index in fit['jumps']:
-                assert jump[jump_index] or jump[jump_index + 1] or jump[jump_index - 1]
+                assert jump[jump_index] or jump[jump_index - 1] or jump[jump_index + 1]
 
             # Test that the correct indexes are recorded
-            ramp_indicies = []
+            ramp_indices = []
             for ramp_index in fit["index"]:
                 # Note start/end of a ramp_index are inclusive meaning that end
                 #    is an index included in the ramp_index so the range is to end + 1
-                new_indicies = list(range(ramp_index["start"], ramp_index["end"] + 1))
+                new_indices = list(range(ramp_index["start"], ramp_index["end"] + 1))
 
                 # check that all the ramps are non-overlapping
-                assert set(ramp_indicies).isdisjoint(new_indicies)
+                assert set(ramp_indices).isdisjoint(new_indices)
 
-                ramp_indicies.extend(new_indicies)
+                ramp_indices.extend(new_indices)
 
             # check that no ramp_index is a jump
-            assert set(ramp_indicies).isdisjoint(fit['jumps'])
+            assert set(ramp_indices).isdisjoint(fit['jumps'])
 
             # check that all resultant indicies are either in a ramp or listed as a jump
-            assert set(ramp_indicies).union(fit['jumps']) == set(range(len(read_pattern)))
+            assert set(ramp_indices).union(fit['jumps']) == set(range(len(read_pattern)))
 
     # Check that the slopes have been estimated reasonably well
     #   There are not that many pixels to test this against and many resultants