From 0c7cb32c95c2ac4f990141838bfeabb3b4e711d6 Mon Sep 17 00:00:00 2001 From: Jonathan Eisenhamer Date: Thu, 10 Aug 2023 22:34:39 -0400 Subject: [PATCH] RCAL-511 Implement Uneven ramp fitting (#175) Changes: * initial import of the MA_TABLE version of ramp fitting from romanisim * Import ramp tests from romanisim * initial ma ramp fitting tests implemented and passing * fix circular import Changes: - fix circular import between `matable_fit` and `matable_fit_cas2022`. - Update tests from `romanisim` to work completely within `stcal`. - Allow test to use `romanisim` if present. * rename algorithm from matable to ols_cas21 * implement conversion from read pattern to multi-accum table * add conversion from multi-accum tables to read patterns * Update fitting routines to take both multi-accum and read pattern * fixed doc test format * add read pattern to RampData * refactor to use match statements * WIP: Integrating CAS21 algorithm into the ramp_fit api Initial implementation done. However, untested and tests need to be implemented. * Revert "WIP: Integrating CAS21 algorithm into the ramp_fit api" This reverts commit a34fc593ef5fff9e6bbc19b752a0e58e989bcff5. * Revert "refactor to use match statements" This reverts commit 1d310a222e407b3a82a6c66c7d24221d819d1ad2. * Revert "add read pattern to RampData" This reverts commit 63512ec7f0aa855d63f549f60ec9dca6431de60d. * fix (most) ruff issues and tox setup * remove stray import * update changelog * Rename references to the Casertano to year 2022 * update ramp fitting from romanisim memory efficiency work * Raise exception when read_pattern does not match resultant data * remove code unused by the ramp fitting code * Parameterize the read time * Remove romanisim test dependency and protect against complex results * simplify the variances return * handle special case where there are no or only one resultant --------- --- CHANGES.rst | 5 + pyproject.toml | 2 + setup.py | 13 ++ src/stcal/ramp_fitting/ols_cas22.pyx | 220 ++++++++++++++++++++ src/stcal/ramp_fitting/ols_cas22_fit.py | 243 +++++++++++++++++++++++ src/stcal/ramp_fitting/ols_cas22_util.py | 164 +++++++++++++++ tests/test_ramp_fitting_cas22.py | 107 ++++++++++ tox.ini | 1 + 8 files changed, 755 insertions(+) create mode 100644 setup.py create mode 100644 src/stcal/ramp_fitting/ols_cas22.pyx create mode 100644 src/stcal/ramp_fitting/ols_cas22_fit.py create mode 100644 src/stcal/ramp_fitting/ols_cas22_util.py create mode 100644 tests/test_ramp_fitting_cas22.py diff --git a/CHANGES.rst b/CHANGES.rst index 4a817988..e1fc4f2f 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -12,6 +12,11 @@ jump 1.4.2 (2023-07-11) ================== +ramp_fitting +~~~~~~~~~~~~ + +- Implement the Casertano, et.al, 2022 uneven ramp fitting [#175] + Bug Fixes --------- diff --git a/pyproject.toml b/pyproject.toml index 58c16764..c336af4d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -45,6 +45,8 @@ requires = [ 'setuptools >=61', 'setuptools_scm[toml] >=3.4', 'wheel', + 'Cython >=0.29.21', + 'numpy >=1.18', ] build-backend = 'setuptools.build_meta' diff --git a/setup.py b/setup.py new file mode 100644 index 00000000..a3b58d91 --- /dev/null +++ b/setup.py @@ -0,0 +1,13 @@ +from setuptools import setup, Extension +from Cython.Build import cythonize +from Cython.Compiler import Options +import numpy as np + +Options.docstrings = True +Options.annotate = False + +extensions = [Extension('stcal.ramp_fitting.ols_cas22', + ['src/stcal/ramp_fitting/ols_cas22.pyx'], + include_dirs=[np.get_include()])] + +setup(ext_modules=cythonize(extensions)) diff --git a/src/stcal/ramp_fitting/ols_cas22.pyx b/src/stcal/ramp_fitting/ols_cas22.pyx new file mode 100644 index 00000000..a57960b1 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22.pyx @@ -0,0 +1,220 @@ +from libc.math cimport sqrt, fabs +import numpy as np +cimport numpy as np +cimport cython + +from stcal.ramp_fitting.ols_cas22_util import ma_table_to_tau, ma_table_to_tbar + +# Casertano+2022, Table 2 +cdef float[2][6] PTABLE = [ + [-np.inf, 5, 10, 20, 50, 100], + [0, 0.4, 1, 3, 6, 10]] +cdef int PTABLE_LENGTH = 6 + + +cdef inline float get_weight_power(float s): + cdef int ise + for i in range(PTABLE_LENGTH): + if s < PTABLE[0][i]: + return PTABLE[1][i - 1] + return PTABLE[1][i] + + +@cython.boundscheck(False) +@cython.wraparound(False) +@cython.cdivision(True) +cdef inline (float, float, float) fit_one_ramp( + float [:] resultants, int start, int end, float read_noise, + float [:] tbar, float [:] tau, int [:] nn): + """Fit a portion of single ramp using the Casertano+22 algorithm. + + Parameters + ---------- + resultants : float [:] + array of resultants for single pixel + start : int + starting point of portion to fit within this pixel + end : int + ending point of portion to fit within this pixel + read_noise : float + read noise for this pixel + tbar : float [:] + mean times of resultants + tau : float [:] + variance weighted mean times of resultants + nn : int [:] + number of reads contributing to reach resultant + + Returns + ------- + slope : float + fit slope + slopereadvar : float + read noise induced variance in slope + slopepoissonvar : float + coefficient of Poisson-noise induced variance in slope + multiply by true flux to get actual Poisson variance. + """ + cdef int i = 0, j = 0 + cdef int nres = end - start + 1 + cdef float ww[2048] + cdef float kk[2048] + cdef float slope = 0, slopereadvar = 0, slopepoissonvar = 0 + cdef float tbarmid = (tbar[start] + tbar[end]) / 2 + + # Casertano+2022 Eq. 44 + # Note we've departed from Casertano+22 slightly; + # there s is just resultants[end]. But that doesn't seem good if, e.g., + # a CR in the first resultant has boosted the whole ramp high but there + # is no actual signal. + cdef float s = max(resultants[end] - resultants[start], 0) + s = s / sqrt(read_noise**2 + s) + cdef float weight_power = get_weight_power(s) + + # It's easy to use up a lot of dynamic range on something like + # (tbar - tbarmid) ** 10. Rescale these. + cdef float tscale = (tbar[end] - tbar[start]) / 2 + if tscale == 0: + tscale = 1 + cdef float f0 = 0, f1 = 0, f2 = 0 + + # Special case where there is no or one resultant, there is no fit. + if nres <= 1: + return 0, 0, 0 + + # Else, do the fitting. + with cython.cpow(True): # Issue when tbar[] == tbarmid causes exception otherwise + for i in range(nres): + # Casertano+22, Eq. 45 + ww[i] = ((((1 + weight_power) * nn[start + i]) / + (1 + weight_power * nn[start + i])) * + fabs((tbar[start + i] - tbarmid) / tscale) ** weight_power) + # Casertano+22 Eq. 35 + f0 += ww[i] + f1 += ww[i] * tbar[start + i] + f2 += ww[i] * tbar[start + i]**2 + # Casertano+22 Eq. 36 + cdef float dd = f2 * f0 - f1 ** 2 + if dd == 0: + return (0.0, 0.0, 0.0) + for i in range(nres): + # Casertano+22 Eq. 37 + kk[i] = (f0 * tbar[start + i] - f1) * ww[i] / dd + for i in range(nres): + # Casertano+22 Eq. 38 + slope += kk[i] * resultants[start + i] + # Casertano+22 Eq. 39 + slopereadvar += kk[i] ** 2 * read_noise ** 2 / nn[start + i] + # Casertano+22 Eq 40 + slopepoissonvar += kk[i] ** 2 * tau[start + i] + for j in range(i + 1, nres): + slopepoissonvar += 2 * kk[i] * kk[j] * tbar[start + i] + + return (slope, slopereadvar, slopepoissonvar) + + +@cython.boundscheck(False) +@cython.wraparound(False) +def fit_ramps(np.ndarray[float, ndim=2] resultants, + np.ndarray[int, ndim=2] dq, + np.ndarray[float, ndim=1] read_noise, read_time, + ma_table): + """Fit ramps using the Casertano+22 algorithm. + + This implementation fits all ramp segments between bad pixels + marked in the dq image with values not equal to zero. So the + number of fit ramps can be larger than the number of pixels. + The derived slopes, corresponding variances, and the locations of + the ramps in each pixel are given in the returned dictionary. + + Parameters + ---------- + resultants : np.ndarry[nresultants, npixel] + the resultants in electrons + dq : np.ndarry[nresultants, npixel] + the dq array. dq != 0 implies bad pixel / CR. + read noise : float + the read noise in electrons + read_time : float + Time to perform a readout. For Roman data, this is FRAME_TIME. + ma_table : list[list[int]] + the ma table prescription + + Returns + ------- + dictionary containing the following keywords: + slope : np.ndarray[nramp] + slopes fit for each ramp + slopereadvar : np.ndarray[nramp] + variance in slope due to read noise + slopepoissonvar : np.ndarray[nramp] + variance in slope due to Poisson noise, divided by the slope + i.e., the slope poisson variance is coefficient * flux; this term + is the coefficient. + pix : np.ndarray[nramp] + the pixel each ramp is in + resstart : np.ndarray[nramp] + The first resultant in this ramp + resend : np.ndarray[nramp] + The last resultant in this ramp. + """ + cdef int nresultant = len(ma_table) + if nresultant != resultants.shape[0]: + raise RuntimeError(f'MA table length {nresultant} does not match number of resultants {resultants.shape[0]}') + + cdef np.ndarray[int] nn = np.array([x[1] for x in ma_table]).astype('i4') + # number of reads in each resultant + cdef np.ndarray[float] tbar = ma_table_to_tbar(ma_table, read_time).astype('f4') + cdef np.ndarray[float] tau = ma_table_to_tau(ma_table, read_time).astype('f4') + cdef int npixel = resultants.shape[1] + cdef int nramp = (np.sum(dq[0, :] == 0) + + np.sum((dq[:-1, :] != 0) & (dq[1:, :] == 0))) + cdef np.ndarray[float] slope = np.zeros(nramp, dtype='f4') + cdef np.ndarray[float] slopereadvar = np.zeros(nramp, dtype='f4') + cdef np.ndarray[float] slopepoissonvar = np.zeros(nramp, dtype='f4') + cdef np.ndarray[int] resstart = np.zeros(nramp, dtype='i4') - 1 + cdef np.ndarray[int] resend = np.zeros(nramp, dtype='i4') - 1 + cdef np.ndarray[int] pix = np.zeros(nramp, dtype='i4') - 1 + cdef int i, j + cdef int inramp = -1 + cdef int rampnum = 0 + for i in range(npixel): + inramp = 0 + for j in range(nresultant): + if (not inramp) and (dq[j, i] == 0): + inramp = 1 + pix[rampnum] = i + resstart[rampnum] = j + elif (not inramp) and (dq[j, i] != 0): + continue + elif inramp and (dq[j, i] == 0): + continue + elif inramp and (dq[j, i] != 0): + inramp = 0 + resend[rampnum] = j - 1 + rampnum += 1 + else: + raise ValueError('unhandled case') + if inramp: + resend[rampnum] = j + rampnum += 1 + # we should have just filled out the starting and stopping locations + # of each ramp. + + cdef float slope0, slopereadvar0, slopepoissonvar0 + cdef float [:, :] resview = resultants + cdef float [:] rnview = read_noise + cdef float [:] tbarview = tbar + cdef float [:] tauview = tau + cdef int [:] nnview = nn + for i in range(nramp): + slope0, slopereadvar0, slopepoissonvar0 = fit_one_ramp( + resview[:, pix[i]], resstart[i], resend[i], rnview[pix[i]], + tbarview, tauview, nnview) + slope[i] = slope0 + slopereadvar[i] = slopereadvar0 + slopepoissonvar[i] = slopepoissonvar0 + + return dict(slope=slope, slopereadvar=slopereadvar, + slopepoissonvar=slopepoissonvar, + pix=pix, resstart=resstart, resend=resend) diff --git a/src/stcal/ramp_fitting/ols_cas22_fit.py b/src/stcal/ramp_fitting/ols_cas22_fit.py new file mode 100644 index 00000000..daadea6f --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22_fit.py @@ -0,0 +1,243 @@ +"""Ramp fitting routines. + +The simulator need not actually fit any ramps, but we would like to do a good +job simulating the noise induced by ramp fitting. That requires computing the +covariance matrix coming out of ramp fitting. But that's actually a big part +of the work of ramp fitting. + +There are a few different proposed ramp fitting algorithms, differing in their +weights. The final derived covariances are all somewhat similarly difficult +to compute, however, since we ultimately end up needing to compute + +.. math:: (A^T C^{-1} A)^{-1} + +for the "optimal" case, or + +.. math:: (A^T W^{-1} A)^{-1} A^T W^{-1} C W^{-1} A (A^T W^{-1} A)^{-1} + +for some alternative weighting. + +We start trying the "optimal" case below. + +For the "optimal" case, a challenge is that we don't want to compute +:math:`C^{-1}` for every pixel individually. Fortunately, we only +need :math:`(A^T C^{-1} A)^{-1}` (which is only a 2x2 matrix) for variances, +and only :math:`(A^T C^{-1} A)^{-1} A^T C^{-1}` for ramp fitting, which is 2xn. +Both of these matrices are effectively single parameter families, depending +after rescaling by the read noise only on the ratio of the read noise and flux. + +So the routines in these packages construct these different matrices, store +them, and interpolate between them for different different fluxes and ratios. +""" +from astropy import units as u +import numpy as np + +from . import ols_cas22 +from .ols_cas22_util import ma_table_to_tau, ma_table_to_tbar, readpattern_to_matable + + +def fit_ramps_casertano(resultants, dq, read_noise, read_time, ma_table=None, read_pattern=None): + """Fit ramps following Casertano+2022, including averaging partial ramps. + + Ramps are broken where dq != 0, and fits are performed on each sub-ramp. + Resultants containing multiple ramps have their ramp fits averaged using + inverse variance weights based on the variance in the individual slope fits + due to read noise. + + Parameters + ---------- + resultants : np.ndarry[nresultants, ...] + the resultants in electrons + dq : np.ndarry[nresultants, ...] + the dq array. dq != 0 implies bad pixel / CR. + read_noise : float + the read noise in electrons + read_time : float + Read time. For Roman data this is the FRAME_TIME keyword. + ma_table : list[list[int]] or None + The MA table prescription. If None, use `read_pattern`. + One of `ma_table` or `read_pattern` must be defined. + read_pattern : list[list[int]] or None + The read pattern prescription. If None, use `ma_table`. + One of `ma_table` or `read_pattern` must be defined. + + Returns + ------- + par : np.ndarray[..., 2] (float) + the best fit pedestal and slope for each pixel + var : np.ndarray[..., 3, 2, 2] (float) + the covariance matrix of par, for each of three noise terms: + the read noise, Poisson source noise, and total noise. + """ + + # Get the Multi-accum table, either as given or from the read pattern + if ma_table is None: + if read_pattern is not None: + ma_table = readpattern_to_matable(read_pattern) + if ma_table is None: + raise RuntimeError('One of `ma_table` or `read_pattern` must be given.') + + resultants_unit = getattr(resultants, 'unit', None) + if resultants_unit is not None: + resultants = resultants.to(u.electron).value + + resultants = np.array(resultants).astype('f4') + + dq = np.array(dq).astype('i4') + + if np.ndim(read_noise) <= 1: + read_noise = read_noise * np.ones(resultants.shape[1:]) + read_noise = np.array(read_noise).astype('f4') + + origshape = resultants.shape + if len(resultants.shape) == 1: + # single ramp. + resultants = resultants.reshape(origshape + (1,)) + dq = dq.reshape(origshape + (1,)) + read_noise = read_noise.reshape(origshape[1:] + (1,)) + + rampfitdict = ols_cas22.fit_ramps( + resultants.reshape(resultants.shape[0], -1), + dq.reshape(resultants.shape[0], -1), + read_noise.reshape(-1), + read_time, + ma_table) + + par = np.zeros(resultants.shape[1:] + (2,), dtype='f4') + var = np.zeros(resultants.shape[1:] + (3,), dtype='f4') + + npix = resultants.reshape(resultants.shape[0], -1).shape[1] + # we need to do some averaging to merge the results in each ramp. + # inverse variance weights based on slopereadvar + weight = ((rampfitdict['slopereadvar'] != 0) / ( + rampfitdict['slopereadvar'] + (rampfitdict['slopereadvar'] == 0))) + totweight = np.bincount(rampfitdict['pix'], weights=weight, minlength=npix) + totval = np.bincount(rampfitdict['pix'], + weights=weight * rampfitdict['slope'], + minlength=npix) + # fill in the averaged slopes + par.reshape(npix, 2)[:, 1] = ( + totval / (totweight + (totweight == 0))) + + # read noise variances + totval = np.bincount( + rampfitdict['pix'], weights=weight ** 2 * rampfitdict['slopereadvar'], + minlength=npix) + var.reshape(npix, 3,)[:, 0] = ( + totval / (totweight ** 2 + (totweight == 0))) + + # poisson noise variances + totval = np.bincount( + rampfitdict['pix'], + weights=weight ** 2 * rampfitdict['slopepoissonvar'], minlength=npix) + var.reshape(npix, 3)[..., 1] = ( + totval / (totweight ** 2 + (totweight == 0))) + + # multiply Poisson term by flux. Clip at zero; no negative Poisson variances. + var[..., 1] *= np.clip(par[..., 1], 0, np.inf) + var[..., 2] = var[..., 0] + var[..., 1] + + if resultants.shape != origshape: + par = par[0] + var = var[0] + + if resultants_unit is not None: + par = par * resultants_unit + + return par, var + + +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/src/stcal/ramp_fitting/ols_cas22_util.py b/src/stcal/ramp_fitting/ols_cas22_util.py new file mode 100644 index 00000000..63bb9a27 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22_util.py @@ -0,0 +1,164 @@ +"""Utility routines for Mutli-Accum Ramp Fitting +""" +import numpy as np + +__all__ = ['ma_table_to_tau', 'ma_table_to_tbar'] + + +def matable_to_readpattern(ma_table): + """Convert read patterns to multi-accum lists + + Using Roman terminology, a "read pattern" is a list of resultants. Each element of this list + is a list of reads that were combined, on-board, to create a resultant. An example read pattern is + + [[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]] + + This pattern has 6 resultants, the first consistent of the first read, the + next consisting of reads 2 and 3, the third consists of read 4, and so on. + + A "Multi-accum table" is a short-hand version of the read pattern. It is a + list of 2-tuples consisting of the following: + + (start_read, n_reads) + + For example, the above read pattern would be represented as, using lists instead of tuples: + + [[1, 1], [2, 2], [4, 1], [5, 4], [9,2], [11,1]] + + The example above, using this function, should perform as follows: + >>> matable_to_readpattern([[1, 1], [2, 2], [4, 1], [5, 4], [9,2], [11,1]]) + [[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]] + + Parameters + ---------- + ma_table : [(first_read, n_reads)[,...]] + The multi-accum table to convert. + + Returns + ------- + read_pattern : [[int[,...]][,...]] + The read pattern that represents the given multi-accum table. + + """ + read_pattern = [list(range(start, start + len)) + for start, len in ma_table] + + return read_pattern + + +def ma_table_to_tau(ma_table, read_time): + """Construct the tau for each resultant from an ma_table. + + .. math:: \\tau = \\overline{t} - (n - 1)(n + 1)\\delta t / 6n + + following Casertano (2022). + + Parameters + ---------- + ma_table : list[list] + List of lists specifying the first read and the number of reads in each + resultant. + + read_time : float + Time to perform a read out. For Roman data, this is FRAME_TIME. + + Returns + ------- + :math:`\\tau` + A time scale appropriate for computing variances. + """ + + meantimes = ma_table_to_tbar(ma_table, read_time) + nreads = np.array([x[1] for x in ma_table]) + return meantimes - (nreads - 1) * (nreads + 1) * read_time / 6 / nreads + + +def ma_table_to_tij(ma_table, read_time): + """Get the times of each read going into resultants for a MA table. + + Currently only ma_table_number = 1 is supported, corresponding to a simple + fiducial high latitude imaging MA table. + + This presently uses a hard-coded, somewhat inflexible MA table description + in the parameters file. But that seems like an okay option given that the + current 'official' file is slated for redesign when the format is relaxed. + + Parameters + ---------- + ma_table : list[list] + A list of (first_read, n_reads) tuples going into resultants. + + read_time : float + The time taken for a read-out. For Roman, this is FRAME_TIME. + + Returns + ------- + list[list[float]] + list of list of readout times for each read entering a resultant + """ + tij = [read_time * np.arange(f, f + n) for (f, n) in ma_table] + return tij + + +def ma_table_to_tbar(ma_table, read_time): + """Construct the mean times for each resultant from an ma_table. + + Parameters + ---------- + ma_table : list[list] + List of lists specifying the first read and the number of reads in each + resultant. + + Returns + ------- + tbar : np.ndarray[n_resultant] (float) + The mean time of the reads of each resultant. + """ + firstreads = np.array([x[0] for x in ma_table]) + nreads = np.array([x[1] for x in ma_table]) + meantimes = read_time * firstreads + read_time * (nreads - 1) / 2 + # at some point I need to think hard about whether the first read has + # slightly less exposure time than all other reads due to the read/reset + # time being slightly less than the read time. + return meantimes + + +def readpattern_to_matable(read_pattern): + """Convert read patterns to multi-accum lists + + Using Roman terminology, a "read pattern" is a list of resultants. Each element of this list + is a list of reads that were combined, on-board, to create a resultant. An example read pattern is + + [[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]] + + This pattern has 6 resultants, the first consistent of the first read, the + next consisting of reads 2 and 3, the third consists of read 4, and so on. + + A "Multi-accum table" is a short-hand version of the read pattern. It is a + list of 2-tuples consisting of the following: + + (start_read, n_reads) + + For example, the above read pattern would be represented as, using lists instead of tuples: + + [[1, 1], [2, 2], [4, 1], [5, 4], [9,2], [11,1]] + + The example above, using this function, should perform as follows: + >>> readpattern_to_matable([[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]]) + [[1, 1], [2, 2], [4, 1], [5, 4], [9, 2], [11, 1]] + + Parameters + ---------- + read_pattern : [[int[,...]][,...]] + The read pattern to convert. + + Returns + ------- + ma_table : [(first_read, n_reads)[,...]] + The multi-accum table that represents the given read pattern. + + """ + ma_table = [[resultant[0], len(resultant)] + for resultant in read_pattern] + + return ma_table diff --git a/tests/test_ramp_fitting_cas22.py b/tests/test_ramp_fitting_cas22.py new file mode 100644 index 00000000..4eb63aa7 --- /dev/null +++ b/tests/test_ramp_fitting_cas22.py @@ -0,0 +1,107 @@ + +""" +Unit tests for ramp-fitting functions. +""" +import numpy as np + +from stcal.ramp_fitting import ols_cas22_fit as ramp +from stcal.ramp_fitting import ols_cas22_util + +# Read Time in seconds +# For Roman, the read time of the detectors is a fixed value and is currently +# backed into code. Will need to refactor to consider the more general case. +# Used to deconstruct the MultiAccum tables into integration times. +ROMAN_READ_TIME = 3.04 + + +def test_matable_to_readpattern(): + """Test conversion from read pattern to multi-accum table""" + ma_table = [[1, 1], [2, 2], [4, 1], [5, 4], [9,2], [11,1]] + expected = [[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]] + + result = ols_cas22_util.matable_to_readpattern(ma_table) + + assert result == expected + + +def test_readpattern_to_matable(): + """Test conversion from read pattern to multi-accum table""" + pattern = [[1], [2, 3], [4], [5, 6, 7, 8], [9, 10], [11]] + expected = [[1, 1], [2, 2], [4, 1], [5, 4], [9,2], [11,1]] + + result = ols_cas22_util.readpattern_to_matable(pattern) + + assert result == expected + + +def test_simulated_ramps(): + ntrial = 100000 + ma_table, flux, read_noise, resultants = simulate_many_ramps(ntrial=ntrial) + + par, var = ramp.fit_ramps_casertano( + resultants, resultants * 0, read_noise, ROMAN_READ_TIME, ma_table=ma_table) + chi2dof_slope = np.sum((par[:, 1] - flux)**2 / var[:, 2]) / ntrial + assert np.abs(chi2dof_slope - 1) < 0.03 + + # now let's mark a bunch of the ramps as compromised. + bad = np.random.uniform(size=resultants.shape) > 0.7 + dq = resultants * 0 + bad + par, var = ramp.fit_ramps_casertano( + resultants, dq, read_noise, ROMAN_READ_TIME, ma_table=ma_table) + # only use okay ramps + # ramps passing the below criterion have at least two adjacent valid reads + # i.e., we can make a measurement from them. + m = np.sum((dq[1:, :] == 0) & (dq[:-1, :] == 0), axis=0) != 0 + chi2dof_slope = np.sum((par[m, 1] - flux)**2 / var[m, 2]) / np.sum(m) + assert np.abs(chi2dof_slope - 1) < 0.03 + assert np.all(par[~m, 1] == 0) + assert np.all(var[~m, 1] == 0) + + +# ######### +# Utilities +# ######### +def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, ma_table=None): + """Simulate many ramps with a particular flux, read noise, and ma_table. + + To test ramp fitting, it's useful to be able to simulate a large number + of ramps that are identical up to noise. This function does that. + + Parameters + ---------- + ntrial : int + number of ramps to simulate + flux : float + flux in electrons / s + read_noise : float + read noise in electrons + ma_table : list[list] (int) + list of lists indicating first read and number of reads in each + resultant + + Returns + ------- + ma_table : list[list] (int) + ma_table used + flux : float + flux used + readnoise : float + read noise used + resultants : np.ndarray[n_resultant, ntrial] (float) + simulated resultants +""" + if ma_table is None: + ma_table = [[1, 4], [5, 1], [6, 3], [9, 10], [19, 3], [22, 15]] + nread = np.array([x[1] for x in ma_table]) + tij = ols_cas22_util.ma_table_to_tij(ma_table, ROMAN_READ_TIME) + resultants = np.zeros((len(ma_table), ntrial), dtype='f4') + buf = np.zeros(ntrial, dtype='i4') + for i, ti in enumerate(tij): + subbuf = np.zeros(ntrial, dtype='i4') + for t0 in ti: + buf += np.random.poisson(ROMAN_READ_TIME * flux, ntrial) + subbuf += buf + resultants[i] = (subbuf / len(ti)).astype('f4') + resultants += np.random.randn(len(ma_table), ntrial) * ( + readnoise / np.sqrt(nread)).reshape(len(ma_table), 1) + return (ma_table, flux, readnoise, resultants) diff --git a/tox.ini b/tox.ini index 6c4a7ac1..ce419dcf 100644 --- a/tox.ini +++ b/tox.ini @@ -57,6 +57,7 @@ deps = jwst: jwst[test] @ git+https://github.com/spacetelescope/jwst.git romancal: romancal[test] @ git+https://github.com/spacetelescope/romancal.git oldestdeps: minimum_dependencies +use_develop = true pass_env = CI set_env =