diff --git a/src/stcal/ramp_fitting/ols_cas22_fit.py b/src/stcal/ramp_fitting/ols_cas22_fit.py index afc5b908e..bf4963567 100644 --- a/src/stcal/ramp_fitting/ols_cas22_fit.py +++ b/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 diff --git a/tests/test_jump_cas22.py b/tests/test_jump_cas22.py index f5adc45b1..a7acf0107 100644 --- a/tests/test_jump_cas22.py +++ b/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