From f305d704148b544fcffe5b3c0f3f53e7c614c717 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Thu, 5 Oct 2023 14:51:37 -0400 Subject: [PATCH] Use new correction factor formulation --- src/stcal/ramp_fitting/ols_cas22/_pixel.pxd | 2 +- src/stcal/ramp_fitting/ols_cas22/_pixel.pyx | 37 +++++++++------------ tests/test_jump_cas22.py | 15 +++++---- 3 files changed, 24 insertions(+), 30 deletions(-) diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd b/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd index a909b3465..326c1d1f5 100644 --- a/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd +++ b/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) diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx index b5968f4de..e376123e6 100644 --- a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx +++ b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx @@ -220,29 +220,23 @@ 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) @@ -250,13 +244,12 @@ cdef class Pixel: 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) diff --git a/tests/test_jump_cas22.py b/tests/test_jump_cas22.py index 7cfaf61dd..76eca2589 100644 --- a/tests/test_jump_cas22.py +++ b/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