Remove pure python implementation of cas22 ramp fitting

src/stcal/ramp_fitting/ols_cas22_fit.py

@@ -110,99 +110,3 @@ def fit_ramps_casertano(
         parameters = parameters * resultants_unit
 
     return parameters, variances
-
-
-def fit_ramps_casertano_no_dq(resultants, read_noise, ma_table):
-    """Fit ramps following Casertano+2022, only using full ramps.
-
-    This is a simpler implementation of fit_ramps_casertano, which doesn't
-    address the case of partial ramps broken by CRs.  This case is easier
-    and can be done reasonably efficiently in pure python; results can be
-    compared with fit_ramps_casertano in for the case of unbroken ramps.
-
-    Parameters
-    ----------
-    resultants : np.ndarry[nresultants, npixel]
-        the resultants in electrons
-    read noise: float
-        the read noise in electrons
-    ma_table : list[list[int]]
-        the ma table prescription
-
-    Returns
-    -------
-    par : np.ndarray[nx, ny, 2] (float)
-        the best fit pedestal and slope for each pixel
-    var : np.ndarray[nx, ny, 3, 2, 2] (float)
-        the covariance matrix of par, for each of three noise terms:
-        the read noise, Poisson source noise, and total noise.
-    """
-    nadd = len(resultants.shape) - 1
-    if np.ndim(read_noise) <= 1:
-        read_noise = np.array(read_noise).reshape((1,) * nadd)
-    smax = resultants[-1]
-    s = smax / np.sqrt(read_noise**2 + smax)  # Casertano+2022 Eq. 44
-    ptable = np.array([  # Casertano+2022, Table 2
-        [-np.inf, 0], [5, 0.4], [10, 1], [20, 3], [50, 6], [100, 10]])
-    pp = ptable[np.searchsorted(ptable[:, 0], s) - 1, 1]
-    nn = np.array([x[1] for x in ma_table])  # number of reads in each resultant
-    tbar = ma_table_to_tbar(ma_table)
-    tau = ma_table_to_tau(ma_table)
-    tbarmid = (tbar[0] + tbar[-1]) / 2
-    if nadd > 0:
-        newshape = ((-1,) + (1,) * nadd)
-        nn = nn.reshape(*newshape)
-        tbar = tbar.reshape(*newshape)
-        tau = tau.reshape(*newshape)
-        tbarmid = tbarmid.reshape(*newshape)
-    ww = (  # Casertano+22, Eq. 45
-        (1 + pp)[None, ...] * nn
-        / (1 + pp[None, ...] * nn)
-        * np.abs(tbar - tbarmid) ** pp[None, ...])
-
-    # Casertano+22 Eq. 35
-    f0 = np.sum(ww, axis=0)
-    f1 = np.sum(ww * tbar, axis=0)
-    f2 = np.sum(ww * tbar**2, axis=0)
-    # Casertano+22 Eq. 36
-    dd = f2 * f0 - f1 ** 2
-    bad = dd == 0
-    dd[bad] = 1
-    # Casertano+22 Eq. 37
-    kk = (f0[None, ...] * tbar - f1[None, ...]) * ww / (
-        dd[None, ...])
-    # shape: [n_resultant, ny, nx]
-    ff = np.sum(kk * resultants, axis=0)  # Casertano+22 Eq. 38
-    # Casertano+22 Eq. 39
-    vr = np.sum(kk**2 / nn, axis=0) * read_noise**2
-    # Casertano+22 Eq. 40
-    vs1 = np.sum(kk**2 * tau, axis=0)
-    vs2inner = np.cumsum(kk * tbar, axis=0)
-    vs2inner = np.concatenate([0 * vs2inner[0][None, ...], vs2inner[:-1, ...]], axis=0)
-    vs2 = 2 * np.sum(vs2inner * kk, axis=0)
-    # sum_{i=1}^{j-1} K_i \bar{t}_i
-    # this is the inner of the two sums in the 2nd term of Eq. 40
-    # Casertano+22 has some discussion of whether it's more efficient to do
-    # this as an explicit double sum or to construct the inner sum separately.
-    # We've made a lot of other quantities that are [nr, ny, nx] in size,
-    # so I don't feel bad about making another.  Clearly a memory optimized
-    # code would work a lot harder to reuse a lot of variables above!
-
-    vs = (vs1 + vs2) * ff
-    vs = np.clip(vs, 0, np.inf)
-    # we can estimate negative flux, but we really shouldn't add variance for
-    # that case!
-
-    # match return values from RampFitInterpolator.fit_ramps
-    # we haven't explicitly calculated here the pedestal, its
-    # uncertainty, or covariance terms.  We just fill
-    # with zeros.
-
-    par = np.zeros(ff.shape + (2,), dtype='f4')
-    var = np.zeros(ff.shape + (3, 2, 2), dtype='f4')
-    par[..., 1] = ff
-    var[..., 0, 1, 1] = vr
-    var[..., 1, 1, 1] = vs
-    var[..., 2, 1, 1] = vr + vs
-
-    return par, var

tests/test_jump_cas22.py

@@ -118,7 +118,7 @@ def ramp_data(base_ramp_data):
             The read pattern used for testing
         t_bar:
             The t_bar values for the read pattern
-        tau: 
+        tau:
             The tau values for the read pattern
         n_reads:
             The number of reads for the read pattern
@@ -243,7 +243,7 @@ def pixel_data(ramp_data):
             Resultants for a single pixel
         t_bar:
             The t_bar values for the read pattern used for the resultants
-        tau: 
+        tau:
             The tau values for the read pattern used for the resultants
         n_reads:
             The number of reads for the read pattern used for the resultants
@@ -277,7 +277,8 @@ def test_make_pixel(pixel_data, use_jump):
         double_gen = zip(pixel['local_slopes'][Diff.double], pixel['var_read_noise'][Diff.double])
 
         for index, (local_slope_1, var_read_noise_1) in enumerate(single_gen):
-            assert local_slope_1 == (resultants[index + 1] - resultants[index]) / (t_bar[index + 1] - t_bar[index])
+            assert local_slope_1 == (
+                (resultants[index + 1] - resultants[index]) / (t_bar[index + 1] - t_bar[index]))
             assert var_read_noise_1 == READ_NOISE * (
                 np.float32(1 / n_reads[index + 1]) + np.float32(1 / n_reads[index])
             )
@@ -314,7 +315,6 @@ def detector_data(ramp_data):
             The read noise vector for those pixels
         read_pattern:
             The read pattern used for the resultants
-
     """
     read_pattern, *_ = ramp_data
     read_noise = np.ones(N_PIXELS, dtype=np.float32) * READ_NOISE
@@ -367,7 +367,7 @@ def test_fit_ramps_no_dq(detector_data, use_jump):
         total_var = fit['average']['read_var'] + fit['average']['poisson_var']
         chi2 += (fit['average']['slope'] - FLUX)**2 / total_var
 
-    chi2 /= N_PIXELS 
+    chi2 /= N_PIXELS
 
     assert np.abs(chi2 - 1) < CHI2_TOL
 
@@ -413,6 +413,18 @@ def jump_data(detector_data):
         Note this specifically checks that we can detect jumps in any read, meaning
         it has an insurance check that a jump has been placed in every single
         read position.
+
+    Returns:
+        resultants
+            The resultants for a large number of pixels
+        read_noise:
+            The read noise vector for those pixels
+        read_pattern:
+            The read pattern used for the resultants
+        jump_reads:
+            Index of read where a jump occurs for each pixel
+        jump_resultants:
+            Index of resultant where a jump occurs for each pixel
     """
     resultants, read_noise, read_pattern = detector_data