diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 821f304b..c238dd55 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -32,53 +32,16 @@ jobs: - linux: py39-oldestdeps-cov-xdist - linux: py39-xdist - linux: py310-xdist - - linux: py311-xdist - - macos: py311-xdist - - linux: py3-xdist-cov + - linux: py311-cov-xdist coverage: codecov - data: - name: retrieve WebbPSF data - runs-on: ubuntu-latest - outputs: - data_path: ${{ steps.data.outputs.path }} - webbpsf_path: ${{ steps.webbpsf_path.outputs.path }} - data_hash: ${{ steps.data_hash.outputs.hash }} - steps: - # webbpsf: - - uses: actions/checkout@v3 - with: - fetch-depth: 0 - - id: data - run: | - echo "path=/tmp/data" >> $GITHUB_OUTPUT - echo "webbpsf_url=https://stsci.box.com/shared/static/n1fealx9q0m6sdnass6wnyfikvxtc0zz.gz" >> $GITHUB_OUTPUT - - run: | - mkdir -p tmp/data/ - mkdir -p ${{ steps.data.outputs.path }} - - run: wget ${{ steps.data.outputs.webbpsf_url }} -O tmp/minimal-webbpsf-data.tar.gz - - run: tar -xzvf tmp/minimal-webbpsf-data.tar.gz -C tmp/data/ - - id: data_hash - run: echo "hash=${{ hashFiles( 'tmp/data' ) }}" >> $GITHUB_OUTPUT - - run: mv tmp/data/* ${{ steps.data.outputs.path }} - - uses: actions/cache@v3 - with: - path: ${{ steps.data.outputs.path }} - key: data-${{ steps.data_hash.outputs.hash }} - - id: webbpsf_path - run: echo "path=${{ steps.data.outputs.path }}/webbpsf-data" >> $GITHUB_OUTPUT + - macos: py311-xdist test_downstream: uses: OpenAstronomy/github-actions-workflows/.github/workflows/tox.yml@v1 - needs: [ data ] with: setenv: | - WEBBPSF_PATH: ${{ needs.data.outputs.webbpsf_path }} - CRDS_PATH: /tmp/crds_cache + CRDS_PATH: /tmp/data/crds_cache CRDS_CLIENT_RETRY_COUNT: 3 CRDS_CLIENT_RETRY_DELAY_SECONDS: 20 - WEBBPSF_PATH: ${{ needs.data.outputs.data_path }}/webbpsf-data - cache-path: ${{ needs.data.outputs.data_path }} - cache-key: data-${{ needs.data.outputs.data_hash }} envs: | - - linux: py3-jwst-cov-xdist - - linux: py3-romancal-cov-xdist - + - linux: py311-jwst-cov-xdist + - linux: py311-romancal-cov-xdist diff --git a/.github/workflows/ci_cron.yml b/.github/workflows/ci_cron.yml index f6ade4f3..1b984244 100644 --- a/.github/workflows/ci_cron.yml +++ b/.github/workflows/ci_cron.yml @@ -22,4 +22,3 @@ jobs: - windows: py310-xdist - windows: py311-xdist - windows: py3-xdist - diff --git a/.gitignore b/.gitignore index c4fc59e8..2afbcecc 100644 --- a/.gitignore +++ b/.gitignore @@ -139,6 +139,9 @@ dmypy.json # Cython debug symbols cython_debug/ +src/stcal/ramp_fitting/ols_cas22/*.c +src/stcal/ramp_fitting/ols_cas22/*.cpp +src/stcal/ramp_fitting/ols_cas22/*.html # setuptools-scm generated module src/stcal/_version.py diff --git a/CHANGES.rst b/CHANGES.rst index 66ba6c1f..dd186cf8 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -3,6 +3,12 @@ - Added ``alignment`` sub-package. [#179] +ramp_fitting +------------ + +- Refactor Casertano, et.al, 2022 uneven ramp fitting and incorporate the matching + jump detection algorithm into it. [#215] + Changes to API -------------- diff --git a/setup.py b/setup.py index 3d3d8fb6..b34e7dfa 100644 --- a/setup.py +++ b/setup.py @@ -6,9 +6,31 @@ 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()], - extra_compile_args=['-std=c99'])] +extensions = [ + Extension( + 'stcal.ramp_fitting.ols_cas22._core', + ['src/stcal/ramp_fitting/ols_cas22/_core.pyx'], + include_dirs=[np.get_include()], + language='c++' + ), + Extension( + 'stcal.ramp_fitting.ols_cas22._fixed', + ['src/stcal/ramp_fitting/ols_cas22/_fixed.pyx'], + include_dirs=[np.get_include()], + language='c++' + ), + Extension( + 'stcal.ramp_fitting.ols_cas22._pixel', + ['src/stcal/ramp_fitting/ols_cas22/_pixel.pyx'], + include_dirs=[np.get_include()], + language='c++' + ), + Extension( + 'stcal.ramp_fitting.ols_cas22._fit_ramps', + ['src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx'], + include_dirs=[np.get_include()], + language='c++' + ), +] setup(ext_modules=cythonize(extensions)) diff --git a/src/stcal/ramp_fitting/ols_cas22.pyx b/src/stcal/ramp_fitting/ols_cas22.pyx deleted file mode 100644 index a57960b1..00000000 --- a/src/stcal/ramp_fitting/ols_cas22.pyx +++ /dev/null @@ -1,220 +0,0 @@ -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/__init__.py b/src/stcal/ramp_fitting/ols_cas22/__init__.py new file mode 100644 index 00000000..d07f4ddf --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/__init__.py @@ -0,0 +1,4 @@ +from ._fit_ramps import fit_ramps, RampFitOutputs +from ._core import Parameter, Variance, Diff, RampJumpDQ + +__all__ = ['fit_ramps', 'RampFitOutputs', 'Parameter', 'Variance', 'Diff', 'RampJumpDQ'] diff --git a/src/stcal/ramp_fitting/ols_cas22/_core.pxd b/src/stcal/ramp_fitting/ols_cas22/_core.pxd new file mode 100644 index 00000000..f7fe7877 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_core.pxd @@ -0,0 +1,58 @@ +from libcpp.vector cimport vector +from libcpp.stack cimport stack +from libcpp.deque cimport deque + + +cdef struct RampIndex: + int start + int end + + +cdef struct RampFit: + float slope + float read_var + float poisson_var + + +cdef struct RampFits: + vector[RampFit] fits + vector[RampIndex] index + vector[int] jumps + RampFit average + + +cdef struct ReadPatternMetadata: + vector[float] t_bar + vector[float] tau + vector[int] n_reads + + +cdef struct Thresh: + float intercept + float constant + + +cpdef enum Diff: + single = 0 + double = 1 + + +cpdef enum Parameter: + intercept = 0 + slope = 1 + + +cpdef enum Variance: + read_var = 0 + poisson_var = 1 + total_var = 2 + + +cpdef enum RampJumpDQ: + JUMP_DET = 4 + + +cpdef float threshold(Thresh thresh, float slope) +cdef float get_power(float s) +cdef deque[stack[RampIndex]] init_ramps(int[:, :] dq) +cpdef ReadPatternMetadata metadata_from_read_pattern(list[list[int]] read_pattern, float read_time) diff --git a/src/stcal/ramp_fitting/ols_cas22/_core.pyx b/src/stcal/ramp_fitting/ols_cas22/_core.pyx new file mode 100644 index 00000000..5c10cb4b --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_core.pyx @@ -0,0 +1,274 @@ +""" +Define the basic types and functions for the CAS22 algorithm with jump detection + +Structs +------- + RampIndex + int start: starting index of the ramp in the resultants + int end: ending index of the ramp in the resultants + + Note that the Python range would be [start:end+1] for any ramp index. + RampFit + float slope: slope of a single ramp + float read_var: read noise variance of a single ramp + float poisson_var: poisson noise variance of single ramp + RampFits + vector[RampFit] fits: ramp fits (in time order) for a single pixel + vector[RampIndex] index: ramp indices (in time order) for a single pixel + RampFit average: average ramp fit for a single pixel + ReadPatternMetata + vector[float] t_bar: mean time of each resultant + vector[float] tau: variance time of each resultant + vector[int] n_reads: number of reads in each resultant + + Note that these are entirely computed from the read_pattern and + read_time (which should be constant for a given telescope) for the + given observation. + Thresh + float intercept: intercept of the threshold + float constant: constant of the threshold + +Enums +----- + Diff + This is the enum to track the index for single vs double difference related + computations. + + single: single difference + double: double difference + + Parameter + This is the enum to track the index of the computed fit parameters for + the ramp fit. + + intercept: the intercept of the ramp fit + slope: the slope of the ramp fit + + Variance + This is the enum to track the index of the computed variance values for + the ramp fit. + + read_var: read variance computed + poisson_var: poisson variance computed + total_var: total variance computed (read_var + poisson_var) + + RampJumpDQ + This enum is to specify the DQ flags for Ramp/Jump detection + + JUMP_DET: jump detected + +Functions +--------- + get_power + Return the power from Casertano+22, Table 2 + threshold + Compute jump threshold + - cpdef gives a python wrapper, but the python version of this method + is considered private, only to be used for testing + init_ramps + Find initial ramps for each pixel, accounts for DQ flags + - A python wrapper, _init_ramps_list, that adjusts types so they can + be directly inspected in python exists for testing purposes only. + metadata_from_read_pattern + Read the read pattern and derive the baseline metadata parameters needed + - cpdef gives a python wrapper, but the python version of this method + is considered private, only to be used for testing +""" +from libcpp.stack cimport stack +from libcpp.deque cimport deque +from libc.math cimport log10 + +import numpy as np +cimport numpy as np +cimport cython + +from stcal.ramp_fitting.ols_cas22._core cimport Thresh, ReadPatternMetadata + + +# Casertano+2022, Table 2 +cdef float[2][6] PTABLE = [ + [-np.inf, 5, 10, 20, 50, 100], + [0, 0.4, 1, 3, 6, 10]] + + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef inline float get_power(float signal): + """ + Return the power from Casertano+22, Table 2 + + Parameters + ---------- + signal: float + signal from the resultants + + Returns + ------- + signal power from Table 2 + """ + cdef int i + for i in range(6): + if signal < PTABLE[0][i]: + return PTABLE[1][i - 1] + + return PTABLE[1][i] + + +cpdef inline float threshold(Thresh thresh, float slope): + """ + Compute jump threshold + + Parameters + ---------- + thresh : Thresh + threshold parameters struct + slope : float + slope of the ramp in question + + Returns + ------- + intercept - constant * log10(slope) + """ + return thresh.intercept - thresh.constant * log10(slope) + + +@cython.boundscheck(False) +@cython.wraparound(False) +cdef inline deque[stack[RampIndex]] init_ramps(int[:, :] dq): + """ + Create the initial ramp stack for each pixel + if dq[index_resultant, index_pixel] == 0, then the resultant is in a ramp + otherwise, the resultant is not in a ramp + + Parameters + ---------- + dq : int[n_resultants, n_pixel] + DQ array + + Returns + ------- + deque of stacks of RampIndex objects + - deque with entry for each pixel + Chosen to be deque because need element access to loop + - stack with entry for each ramp found (top of stack is last ramp found) + - RampIndex with start and end indices of the ramp in the resultants + """ + cdef int n_pixel, n_resultants + + n_resultants, n_pixel = np.array(dq).shape + cdef deque[stack[RampIndex]] pixel_ramps + + cdef int index_resultant, index_pixel + cdef stack[RampIndex] ramps + cdef RampIndex ramp + + for index_pixel in range(n_pixel): + ramps = stack[RampIndex]() + + # Note: if start/end are -1, then no value has been assigned + # ramp.start == -1 means we have not started a ramp + # dq[index_resultant, index_pixel] == 0 means resultant is in ramp + ramp = RampIndex(-1, -1) + for index_resultant in range(n_resultants): + if ramp.start == -1: + # Looking for the start of a ramp + if dq[index_resultant, index_pixel] == 0: + # We have found the start of a ramp! + ramp.start = index_resultant + else: + # This is not the start of the ramp yet + continue + else: + # Looking for the end of a ramp + if dq[index_resultant, index_pixel] == 0: + # This pixel is in the ramp do nothing + continue + else: + # This pixel is not in the ramp + # => index_resultant - 1 is the end of the ramp + ramp.end = index_resultant - 1 + + # Add completed ramp to stack and reset ramp + ramps.push(ramp) + ramp = RampIndex(-1, -1) + + # Handle case where last resultant is in ramp (so no end has been set) + if ramp.start != -1 and ramp.end == -1: + # Last resultant is end of the ramp => set then add to stack + ramp.end = n_resultants - 1 + ramps.push(ramp) + + # Add ramp stack for pixel to list + pixel_ramps.push_back(ramps) + + return pixel_ramps + + +def _init_ramps_list(np.ndarray[int, ndim=2] dq): + """ + This is a wrapper for init_ramps so that it can be fully inspected from pure + python. A cpdef cannot be used in that case becase a stack has no direct python + analog. Instead this function turns that stack into a list ordered in the same + order as the stack; meaning that, the first element of the list is the top of + the stack. + Note this function is for testing purposes only and so is marked as private + within this private module + """ + cdef deque[stack[RampIndex]] raw = init_ramps(dq) + + # Have to turn deque and stack into python compatible objects + cdef RampIndex index + cdef stack[RampIndex] ramp + cdef list out = [] + cdef list stack_out + for ramp in raw: + stack_out = [] + while not ramp.empty(): + index = ramp.top() + ramp.pop() + # So top of stack is first item of list + stack_out = [index] + stack_out + + out.append(stack_out) + + return out + + +@cython.boundscheck(False) +@cython.wraparound(False) +cpdef ReadPatternMetadata metadata_from_read_pattern(list[list[int]] read_pattern, float read_time): + """ + Derive the input data from the the read pattern + + read pattern is a list of resultant lists, where each resultant list is + a list of the reads in that resultant. + + Parameters + ---------- + read pattern: list[list[int]] + read pattern for the image + read_time : float + Time to perform a readout. + + Returns + ------- + ReadPatternMetadata struct: + vector[float] t_bar: mean time of each resultant + vector[float] tau: variance time of each resultant + vector[int] n_reads: number of reads in each resultant + """ + cdef int n_resultants = len(read_pattern) + cdef ReadPatternMetadata data = ReadPatternMetadata(vector[float](n_resultants), + vector[float](n_resultants), + vector[int](n_resultants)) + + cdef int index, n_reads + cdef list[int] resultant + for index, resultant in enumerate(read_pattern): + n_reads = len(resultant) + + data.n_reads[index] = n_reads + data.t_bar[index] = read_time * np.mean(resultant) + data.tau[index] = np.sum((2 * (n_reads - np.arange(n_reads)) - 1) * resultant) * read_time / n_reads**2 + + return data diff --git a/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx b/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx new file mode 100644 index 00000000..ad19fb28 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx @@ -0,0 +1,135 @@ +import numpy as np +cimport numpy as np +from libcpp cimport bool +from libcpp.stack cimport stack +from libcpp.list cimport list as cpp_list +from libcpp.deque cimport deque +cimport cython + +from stcal.ramp_fitting.ols_cas22._core cimport (RampFits, RampIndex, Thresh, + metadata_from_read_pattern, init_ramps, + Parameter, Variance, RampJumpDQ) +from stcal.ramp_fitting.ols_cas22._fixed cimport fixed_values_from_metadata, FixedValues +from stcal.ramp_fitting.ols_cas22._pixel cimport make_pixel + +from typing import NamedTuple + + +# Fix the default Threshold values at compile time these values cannot be overridden +# dynamically at runtime. +DEF DefaultIntercept = 5.5 +DEF DefaultConstant = 1/3.0 + +class RampFitOutputs(NamedTuple): + """ + Simple tuple wrapper for outputs from the ramp fitting algorithm + This clarifies the meaning of the outputs via naming them something + descriptive. + + Attributes + ---------- + fits: list of RampFits + the raw ramp fit outputs, these are all structs which will get mapped to + python dictionaries. + parameters: np.ndarray[n_pixel, 2] + the slope and intercept for each pixel's ramp fit. see Parameter enum + for indexing indicating slope/intercept in the second dimension. + variances: np.ndarray[n_pixel, 3] + the read, poisson, and total variances for each pixel's ramp fit. + see Variance enum for indexing indicating read/poisson/total in the + second dimension. + dq: np.ndarray[n_resultants, n_pixel] + the dq array, with additional flags set for jumps detected by the + jump detection algorithm. + """ + fits: list + parameters: np.ndarray + variances: np.ndarray + dq: np.ndarray + + +@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, + float read_time, + list[list[int]] read_pattern, + bool use_jump=False, + float intercept=DefaultIntercept, + float constant=DefaultConstant): + """Fit ramps using the Casertano+22 algorithm. + This implementation uses the Cas22 algorithm to fit ramps, where + ramps are fit between bad resultants marked by dq flags for each pixel + which are not equal to zero. If use_jump is True, it additionally uses + jump detection to mark additional resultants for each pixel as bad if + a jump is suspected in them. + + Parameters + ---------- + resultants : np.ndarry[n_resultants, n_pixel] + the resultants in electrons + dq : np.ndarry[n_resultants, n_pixel] + the dq array. dq != 0 implies bad pixel / CR. + read_noise : np.ndarray[n_pixel] + the read noise in electrons for each pixel + read_time : float + Time to perform a readout. For Roman data, this is FRAME_TIME. + read_pattern : list[list[int]] + the read pattern for the image + use_jump : bool + If True, use the jump detection algorithm to identify CRs. + If False, use the DQ array to identify CRs. + intercept : float + The intercept value for the threshold function. Default=5.5 + constant : float + The constant value for the threshold function. Default=1/3.0 + + Returns + ------- + A RampFitOutputs tuple + """ + cdef int n_pixels, n_resultants + n_resultants = resultants.shape[0] + n_pixels = resultants.shape[1] + + if n_resultants != len(read_pattern): + raise RuntimeError(f'The read pattern length {len(read_pattern)} does not ' + f'match number of resultants {n_resultants}') + + # Pre-compute data for all pixels + cdef FixedValues fixed = fixed_values_from_metadata(metadata_from_read_pattern(read_pattern, read_time), + Thresh(intercept, constant), + use_jump) + + # Compute all the initial sets of ramps + cdef deque[stack[RampIndex]] pixel_ramps = init_ramps(dq) + + # Use list because this might grow very large which would require constant + # reallocation. We don't need random access, and this gets cast to a python + # list in the end. + cdef cpp_list[RampFits] ramp_fits + + cdef np.ndarray[float, ndim=2] parameters = np.zeros((n_pixels, 2), dtype=np.float32) + cdef np.ndarray[float, ndim=2] variances = np.zeros((n_pixels, 3), dtype=np.float32) + + # Perform all of the fits + cdef RampFits fit + cdef int index + for index in range(n_pixels): + # Fit all the ramps for the given pixel + fit = make_pixel(fixed, read_noise[index], + resultants[:, index]).fit_ramps(pixel_ramps[index]) + + parameters[index, Parameter.slope] = fit.average.slope + + variances[index, Variance.read_var] = fit.average.read_var + variances[index, Variance.poisson_var] = fit.average.poisson_var + variances[index, Variance.total_var] = fit.average.read_var + fit.average.poisson_var + + for jump in fit.jumps: + dq[jump, index] = RampJumpDQ.JUMP_DET + + ramp_fits.push_back(fit) + + return RampFitOutputs(ramp_fits, parameters, variances, dq) diff --git a/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd b/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd new file mode 100644 index 00000000..72087051 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd @@ -0,0 +1,21 @@ +from libcpp cimport bool + +from stcal.ramp_fitting.ols_cas22._core cimport Thresh, ReadPatternMetadata + + +cdef class FixedValues: + cdef bool use_jump + cdef ReadPatternMetadata data + cdef Thresh threshold + + cdef float[:, :] t_bar_diffs + cdef float[:, :] t_bar_diff_sqrs + cdef float[:, :] read_recip_coeffs + cdef float[:, :] var_slope_coeffs + + cdef float[:, :] t_bar_diff_vals(FixedValues self) + cdef float[:, :] read_recip_vals(FixedValues self) + cdef float[:, :] var_slope_vals(FixedValues self) + + +cpdef FixedValues fixed_values_from_metadata(ReadPatternMetadata data, Thresh threshold, bool use_jump) diff --git a/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx b/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx new file mode 100644 index 00000000..6bd72b07 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx @@ -0,0 +1,273 @@ +""" +Define the data which is fixed for all pixels to compute the CAS22 algorithm with + jump detection + +Objects +------- +FixedValues : class + Class to contain the data fixed for all pixels and commonly referenced + universal values for jump detection + +Functions +--------- + fixed_values_from_metadata : function + Fast constructor for FixedValues from the read pattern metadata + - cpdef gives a python wrapper, but the python version of this method + is considered private, only to be used for testing +""" +import numpy as np +cimport numpy as np +cimport cython + +from stcal.ramp_fitting.ols_cas22._core cimport Thresh, ReadPatternMetadata, Diff +from stcal.ramp_fitting.ols_cas22._fixed cimport FixedValues + +cdef class FixedValues: + """ + Class to contain all the values which are fixed for all pixels for a given + read pattern. + This class is used to pre-compute these values once so that they maybe + reused for all pixels. This is done for performance reasons. + + Parameters + ---------- + use_jump : bool + flag to indicate whether to use jump detection (user input) + + data : ReadPatternMetadata + Metadata struct created from a read pattern + + threshold : Thresh + Parameterization struct for threshold function + + t_bar_diffs : float[:, :] + These are the differences of t_bar used for jump detection. + single differences of t_bar: + t_bar_diffs[Diff.single, :] = (t_bar[i+1] - t_bar[i]) + double differences of t_bar: + t_bar_diffs[Diff.double, :] = (t_bar[i+2] - t_bar[i]) + t_bar_diff_sqrs : float[:, :] + These are the squared differnences of t_bar used for jump detection. + single differences of t_bar: + t_bar_diff_sqrs[Diff.single, :] = (t_bar[i+1] - t_bar[i])**2 + double differences of t_bar: + t_bar_diff_sqrs[Diff.double, :] = (t_bar[i+2] - t_bar[i])**2 + read_recip_coeffs : float[:, :] + Coefficients for the read noise portion of the variance used to compute + the jump detection statistics. These are formed from the reciprocal sum + of the number of reads. + single sum of reciprocal n_reads: + read_recip_coeffs[Diff.single, :] = ((1/n_reads[i+1]) + (1/n_reads[i])) + double sum of reciprocal n_reads: + read_recip_coeffs[Diff.double, :] = ((1/n_reads[i+2]) + (1/n_reads[i])) + var_slope_coeffs : float[:, :] + Coefficients for the slope portion of the variance used to compute the + jump detection statistics, which happend to be fixed for any given ramp + fit. + single of slope variance term: + var_slope_coeffs[Diff.single, :] = (tau[i] + tau[i+1] + - min(t_bar[i], t_bar[i+1])) + double of slope variance term: + var_slope_coeffs[Diff.double, :] = (tau[i] + tau[i+2] + - min(t_bar[i], t_bar[i+2])) + + Notes + ----- + - t_bar_diffs, t_bar_diff_sqrs, read_recip_coeffs, var_slope_coeffs are only + computed if use_jump is True. These values represent reused computations + for jump detection which are used by every pixel for jump detection. They + are computed once and stored in the FixedValues for reuse by all pixels. + - The computations are done using vectorized operations for some performance + increases. However, this is marginal compaired with the performance increase + from pre-computing the values and reusing them. + """ + + @cython.boundscheck(False) + @cython.wraparound(False) + cdef inline float[:, :] t_bar_diff_vals(FixedValues self): + """ + Compute the difference offset of t_bar + + Returns + ------- + [ + , + , + ] + """ + # Cast vector to memory view + # This way of doing it is potentially memory unsafe because the memory + # can outlive the vector. However, this is much faster (no copies) and + # much simpler than creating an intermediate wrapper which can pretend + # to be a memory view. In this case, I make sure that the memory view + # stays local to the function (numpy operations create brand new objects) + cdef float[:] t_bar = self.data.t_bar.data() + cdef int end = len(t_bar) + + cdef np.ndarray[float, ndim=2] t_bar_diff_vals = np.zeros((2, self.data.t_bar.size() - 1), dtype=np.float32) + + t_bar_diff_vals[Diff.single, :] = np.subtract(t_bar[1:], t_bar[:end - 1]) + t_bar_diff_vals[Diff.double, :end - 2] = np.subtract(t_bar[2:], t_bar[:end - 2]) + t_bar_diff_vals[Diff.double, end - 2] = np.nan # last double difference is undefined + + return t_bar_diff_vals + + @cython.boundscheck(False) + @cython.wraparound(False) + cdef inline float[:, :] read_recip_vals(FixedValues self): + """ + Compute the reciprical sum of the number of reads + + Returns + ------- + [ + <(1/n_reads[i+1] + 1/n_reads[i])>, + <(1/n_reads[i+2] + 1/n_reads[i])>, + ] + + """ + # Cast vector to memory view + # This way of doing it is potentially memory unsafe because the memory + # can outlive the vector. However, this is much faster (no copies) and + # much simpler than creating an intermediate wrapper which can pretend + # to be a memory view. In this case, I make sure that the memory view + # stays local to the function (numpy operations create brand new objects) + cdef int[:] n_reads = self.data.n_reads.data() + cdef int end = len(n_reads) + + cdef np.ndarray[float, ndim=2] read_recip_vals = np.zeros((2, self.data.n_reads.size() - 1), dtype=np.float32) + + read_recip_vals[Diff.single, :] = (np.divide(1.0, n_reads[1:], dtype=np.float32) + + np.divide(1.0, n_reads[:end - 1], dtype=np.float32)) + read_recip_vals[Diff.double, :end - 2] = (np.divide(1.0, n_reads[2:], dtype=np.float32) + + np.divide(1.0, n_reads[:end - 2], dtype=np.float32)) + read_recip_vals[Diff.double, end - 2] = np.nan # last double difference is undefined + + return read_recip_vals + + + @cython.boundscheck(False) + @cython.wraparound(False) + cdef inline float[:, :] var_slope_vals(FixedValues self): + """ + Compute slope part of the jump statistic variances + + Returns + ------- + [ + <(tau[i] + tau[i+1] - min(t_bar[i], t_bar[i+1])) * correction(i, i+1)>, + <(tau[i] + tau[i+2] - min(t_bar[i], t_bar[i+2])) * correction(i, i+2)>, + ] + """ + # Cast vectors to memory views + # This way of doing it is potentially memory unsafe because the memory + # can outlive the vector. However, this is much faster (no copies) and + # much simpler than creating an intermediate wrapper which can pretend + # to be a memory view. In this case, I make sure that the memory view + # stays local to the function (numpy operations create brand new objects) + cdef float[:] t_bar = self.data.t_bar.data() + cdef float[:] tau = self.data.tau.data() + cdef int end = len(t_bar) + + cdef np.ndarray[float, ndim=2] var_slope_vals = np.zeros((2, self.data.t_bar.size() - 1), dtype=np.float32) + + var_slope_vals[Diff.single, :] = (np.add(tau[1:], tau[:end - 1]) - np.minimum(t_bar[1:], t_bar[:end - 1])) + var_slope_vals[Diff.double, :end - 2] = (np.add(tau[2:], tau[:end - 2]) - np.minimum(t_bar[2:], t_bar[:end - 2])) + var_slope_vals[Diff.double, end - 2] = np.nan # last double difference is undefined + + return var_slope_vals + + def _to_dict(FixedValues self): + """ + This is a private method to convert the FixedValues object to a dictionary, + so that attributes can be directly accessed in python. Note that this + is needed because class attributes cannot be accessed on cython classes + directly in python. Instead they need to be accessed or set using a + python compatible method. This method is a pure puthon method bound + to to the cython class and should not be used by any cython code, and + only exists for testing purposes. + """ + cdef np.ndarray[float, ndim=2] t_bar_diffs + cdef np.ndarray[float, ndim=2] t_bar_diff_sqrs + cdef np.ndarray[float, ndim=2] read_recip_coeffs + cdef np.ndarray[float, ndim=2] var_slope_coeffs + + if self.use_jump: + t_bar_diffs = np.array(self.t_bar_diffs, dtype=np.float32) + t_bar_diff_sqrs = np.array(self.t_bar_diff_sqrs, dtype=np.float32) + read_recip_coeffs = np.array(self.read_recip_coeffs, dtype=np.float32) + var_slope_coeffs = np.array(self.var_slope_coeffs, dtype=np.float32) + else: + try: + self.t_bar_diffs + except AttributeError: + t_bar_diffs = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("t_bar_diffs should not exist") + + try: + self.t_bar_diff_sqrs + except AttributeError: + t_bar_diff_sqrs = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("t_bar_diff_sqrs should not exist") + + try: + self.read_recip_coeffs + except AttributeError: + read_recip_coeffs = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("read_recip_coeffs should not exist") + + try: + self.var_slope_coeffs + except AttributeError: + var_slope_coeffs = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("var_slope_coeffs should not exist") + + return dict(data=self.data, + threshold=self.threshold, + t_bar_diffs=t_bar_diffs, + t_bar_diff_sqrs=t_bar_diff_sqrs, + read_recip_coeffs=read_recip_coeffs, + var_slope_coeffs=var_slope_coeffs) + + +cpdef inline FixedValues fixed_values_from_metadata(ReadPatternMetadata data, Thresh threshold, bool use_jump): + """ + Fast constructor for FixedValues class + Use this instead of an __init__ because it does not incure the overhead + of switching back and forth to python + + Parameters + ---------- + data : ReadPatternMetadata + metadata object created from the read pattern (user input) + threshold : Thresh + threshold object (user input) + use_jump : bool + flag to indicate whether to use jump detection (user input) + + Returns + ------- + FixedValues object (with pre-computed values for jump detection if use_jump + is True) + """ + cdef FixedValues fixed = FixedValues() + + # Fill in input information for all pixels + fixed.use_jump = use_jump + fixed.threshold = threshold + + # Cast vector to a c array + fixed.data = data + + # Pre-compute jump detection computations shared by all pixels + if use_jump: + fixed.t_bar_diffs = fixed.t_bar_diff_vals() + fixed.t_bar_diff_sqrs = np.square(fixed.t_bar_diffs, dtype=np.float32) + fixed.read_recip_coeffs = fixed.read_recip_vals() + fixed.var_slope_coeffs = fixed.var_slope_vals() + + return fixed diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd b/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd new file mode 100644 index 00000000..bf390419 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd @@ -0,0 +1,23 @@ +from libcpp.stack cimport stack + +from stcal.ramp_fitting.ols_cas22._core cimport RampFit, RampFits, RampIndex +from stcal.ramp_fitting.ols_cas22._fixed cimport FixedValues + +cdef class Pixel: + cdef FixedValues fixed + cdef float read_noise + cdef float [:] resultants + + cdef float[:, :] local_slopes + cdef float[:, :] var_read_noise + + cdef float[:, :] local_slope_vals(Pixel self) + cdef RampFit fit_ramp(Pixel self, RampIndex ramp) + + 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) + + +cpdef Pixel make_pixel(FixedValues fixed, float read_noise, float [:] resultants) diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx new file mode 100644 index 00000000..88544243 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx @@ -0,0 +1,556 @@ +""" +Define the C class for the Cassertano22 algorithm for fitting ramps with jump detection + +Objects +------- +Pixel : class + Class to handle ramp fit with jump detection for a single pixel + Provides fits method which fits all the ramps for a single pixel + +Functions +--------- + make_pixel : function + Fast constructor for a Pixel class from input data. + - cpdef gives a python wrapper, but the python version of this method + is considered private, only to be used for testing +""" +from libc.math cimport sqrt, fabs +from libcpp.vector cimport vector +from libcpp.stack cimport stack + +import numpy as np +cimport numpy as np +cimport cython + + +from stcal.ramp_fitting.ols_cas22._core cimport get_power, threshold, RampFit, RampFits, RampIndex, Diff +from stcal.ramp_fitting.ols_cas22._fixed cimport FixedValues +from stcal.ramp_fitting.ols_cas22._pixel cimport Pixel + + +cdef class Pixel: + """ + Class to contain the data to fit ramps for a single pixel. + This data is drawn from for all ramps for a single pixel. + This class pre-computes jump detection values shared by all ramps + for a given pixel. + + Parameters + ---------- + fixed : FixedValues + The object containing all the values and metadata which is fixed for a + given read pattern> + read_noise : float + The read noise for the given pixel + resultants : float [:] + Resultants input for the given pixel + + local_slopes : float [:, :] + These are the local slopes between the resultants for the pixel. + single difference local slope: + local_slopes[Diff.single, :] = (resultants[i+1] - resultants[i]) + / (t_bar[i+1] - t_bar[i]) + double difference local slope: + local_slopes[Diff.double, :] = (resultants[i+2] - resultants[i]) + / (t_bar[i+2] - t_bar[i]) + var_read_noise : float [:, :] + The read noise variance term of the jump statistics + single difference read noise variance: + var_read_noise[Diff.single, :] = read_noise * ((1/n_reads[i+1]) + (1/n_reads[i])) + double difference read_noise variance: + var_read_noise[Diff.doule, :] = read_noise * ((1/n_reads[i+2]) + (1/n_reads[i])) + + Notes + ----- + - local_slopes and var_read_noise are only computed if use_jump is True. + These values represent reused computations for jump detection which are + used by every ramp for the given pixel for jump detection. They are + computed once and stored for reuse by all ramp computations for the pixel. + - The computations are done using vectorized operations for some performance + increases. However, this is marginal compaired with the performance increase + from pre-computing the values and reusing them. + + Methods + ------- + fit_ramp (ramp_index) : method + Compute the ramp fit for a single ramp defined by an inputed RampIndex + fit_ramps (ramp_stack) : method + Compute all the ramps for a single pixel using the Casertano+22 algorithm + with jump detection. + """ + + @cython.boundscheck(False) + @cython.wraparound(False) + cdef inline float[:, :] local_slope_vals(Pixel self): + """ + Compute the local slopes between resultants for the pixel + + Returns + ------- + [ + <(resultants[i+1] - resultants[i])> / <(t_bar[i+1] - t_bar[i])>, + <(resultants[i+2] - resultants[i])> / <(t_bar[i+2] - t_bar[i])>, + ] + """ + cdef float[:] resultants = self.resultants + cdef int end = len(resultants) + + # Read the t_bar_diffs into a local variable to avoid calling through Python + # multiple times + cdef np.ndarray[float, ndim=2] t_bar_diffs = np.array(self.fixed.t_bar_diffs, dtype=np.float32) + + cdef np.ndarray[float, ndim=2] local_slope_vals = np.zeros((2, end - 1), dtype=np.float32) + + local_slope_vals[Diff.single, :] = (np.subtract(resultants[1:], resultants[:end - 1]) + / t_bar_diffs[Diff.single, :]).astype(np.float32) + local_slope_vals[Diff.double, :end - 2] = (np.subtract(resultants[2:], resultants[:end - 2]) + / t_bar_diffs[Diff.double, :end-2]).astype(np.float32) + local_slope_vals[Diff.double, end - 2] = np.nan # last double difference is undefined + + return local_slope_vals + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + cdef inline RampFit fit_ramp(Pixel self, RampIndex ramp): + """ + Fit a single ramp using Casertano+22 algorithm. + + Parameters + ---------- + ramp : RampIndex + Struct for start and end of ramp to fit + + Returns + ------- + RampFit struct of slope, read_var, poisson_var + """ + cdef int n_resultants = ramp.end - ramp.start + 1 + + # Special case where there is no or one resultant, there is no fit and + # we bail out before any computations. + # Note that in this case, we cannot compute the slope or the variances + # because these computations require at least two resultants. Therefore, + # this case is degernate and we return NaNs for the values. + if n_resultants <= 1: + return RampFit(np.nan, np.nan, np.nan) + + # Start computing the fit + + # Cast vectors to memory views for faster access + # This way of doing it is potentially memory unsafe because the memory + # can outlive the vector. However, this is much faster (no copies) and + # much simpler than creating an intermediate wrapper which can pretend + # to be a memory view. In this case, I make sure that the memory view + # stays local to the function t_bar, tau, n_reads are used only for + # computations whose results are stored in new objects, so they are local + cdef float[:] t_bar_ = self.fixed.data.t_bar.data() + cdef float[:] tau_ = self.fixed.data.tau.data() + cdef int[:] n_reads_ = self.fixed.data.n_reads.data() + + # Setup data for fitting (work over subset of data) + # Recall that the RampIndex contains the index of the first and last + # index of the ramp. Therefore, the Python slice needed to get all the + # data within the ramp is: + # ramp.start:ramp.end + 1 + cdef float[:] resultants = self.resultants[ramp.start:ramp.end + 1] + cdef float[:] t_bar = t_bar_[ramp.start:ramp.end + 1] + cdef float[:] tau = tau_[ramp.start:ramp.end + 1] + cdef int[:] n_reads = n_reads_[ramp.start:ramp.end + 1] + + # Reference read_noise as a local variable to avoid calling through Python + # every time it is accessed. + cdef float read_noise = self.read_noise + + # Compute mid point time + cdef int end = len(resultants) - 1 + cdef float t_bar_mid = (t_bar[0] + t_bar[end]) / 2 + + # Casertano+2022 Eq. 44 + # Note we've departed from Casertano+22 slightly; + # there s is just resultants[ramp.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[0], 0) + s = s / sqrt(read_noise**2 + s) + cdef float power = get_power(s) + + # It's easy to use up a lot of dynamic range on something like + # (tbar - tbarmid) ** 10. Rescale these. + cdef float t_scale = (t_bar[end] - t_bar[0]) / 2 + t_scale = 1 if t_scale == 0 else t_scale + + # Initalize the fit loop + cdef int i = 0, j = 0 + cdef vector[float] weights = vector[float](n_resultants) + cdef vector[float] coeffs = vector[float](n_resultants) + cdef RampFit ramp_fit = RampFit(0, 0, 0) + cdef float f0 = 0, f1 = 0, f2 = 0 + + # Issue when tbar[] == tbarmid causes exception otherwise + with cython.cpow(True): + for i in range(n_resultants): + # Casertano+22, Eq. 45 + weights[i] = ((((1 + power) * n_reads[i]) / (1 + power * n_reads[i])) * + fabs((t_bar[i] - t_bar_mid) / t_scale) ** power) + + # Casertano+22 Eq. 35 + f0 += weights[i] + f1 += weights[i] * t_bar[i] + f2 += weights[i] * t_bar[i]**2 + + # Casertano+22 Eq. 36 + cdef float det = f2 * f0 - f1 ** 2 + if det == 0: + return ramp_fit + + for i in range(n_resultants): + # Casertano+22 Eq. 37 + coeffs[i] = (f0 * t_bar[i] - f1) * weights[i] / det + + for i in range(n_resultants): + # Casertano+22 Eq. 38 + ramp_fit.slope += coeffs[i] * resultants[i] + + # Casertano+22 Eq. 39 + ramp_fit.read_var += (coeffs[i] ** 2 * read_noise ** 2 / n_reads[i]) + + # Casertano+22 Eq 40 + ramp_fit.poisson_var += coeffs[i] ** 2 * tau[i] + for j in range(i + 1, n_resultants): + ramp_fit.poisson_var += (2 * coeffs[i] * coeffs[j] * t_bar[i]) + + return ramp_fit + + @cython.boundscheck(False) + @cython.wraparound(False) + 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 + slope : float + The computed slope for the ramp + """ + + 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)) + 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])) + + Parameters + ---------- + slope : float + The computed slope for the ramp + 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. + + Returns + ------- + Create a single instance of the stastic for the given parameters + """ + 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)) + correct + + + @cython.boundscheck(False) + @cython.wraparound(False) + cdef inline float[:] stats(Pixel self, float slope, RampIndex ramp): + """ + Compute fit statistics for jump detection on a single ramp + stats[i] = max(stat(i, 0), stat(i, 1)) + Note for i == end - 1, no stat(i, 1) exists, so its just stat(i, 0) + + Parameters + ---------- + slope : float + The computed slope for the ramp + ramp : RampIndex + Struct for start and end of ramp to fit + + Returns + ------- + list of statistics for each resultant + """ + cdef int start = ramp.start # index of first resultant for ramp + cdef int end = ramp.end # index of last resultant for ramp + + # Observe that the length of the ramp's sub array of the resultant would + # be `end - start + 1`. However, we are computing single and double + # "differences" which means we need to reference at least two points in + # this subarray at a time. For the single case, the maximum index allowed + # would be `end - 1`. Observe that `range(start, end)` will iterate over + # `start, start+1, start+1, ..., end-2, end-1` + # as the second argument to the `range` is the first index outside of the + # range + + cdef np.ndarray[float, ndim=1] stats = np.zeros(end - start, dtype=np.float32) + + cdef int index, stat + for stat, index in enumerate(range(start, end)): + if index == end - 1: + # It is not possible to compute double differences for the second + # to last resultant in the ramp. Therefore, we just compute the + # single difference for this resultant. + stats[stat] = self.stat(slope, ramp, index, Diff.single) + else: + stats[stat] = max(self.stat(slope, ramp, index, Diff.single), + self.stat(slope, ramp, index, Diff.double)) + + return stats + + + @cython.boundscheck(False) + @cython.wraparound(False) + @cython.cdivision(True) + cdef inline RampFits fit_ramps(Pixel self, stack[RampIndex] ramps): + """ + Compute all the ramps for a single pixel using the Casertano+22 algorithm + with jump detection. + + Parameters + ---------- + ramps : stack[RampIndex] + Stack of initial ramps to fit for a single pixel + multiple ramps are possible due to dq flags + + Returns + ------- + RampFits struct of all the fits for a single pixel + """ + # Setup algorithm + cdef RampFits ramp_fits + cdef RampIndex ramp + cdef RampFit ramp_fit + + ramp_fits.average.slope = 0 + ramp_fits.average.read_var = 0 + ramp_fits.average.poisson_var = 0 + + cdef float [:] stats + cdef int jump0, jump1 + cdef float weight, total_weight = 0 + + # Run while the stack is non-empty + while not ramps.empty(): + # Remove top ramp of the stack to use + ramp = ramps.top() + ramps.pop() + + # Compute fit + ramp_fit = self.fit_ramp(ramp) + + # Run jump detection if enabled + if self.fixed.use_jump: + stats = self.stats(ramp_fit.slope, ramp) + + # We have to protect against the case where the passed "ramp" is + # only a single point. In that case, stats will be empty. This + # will create an error in the max() call. + if len(stats) > 0 and max(stats) > threshold(self.fixed.threshold, ramp_fit.slope): + # Compute jump point to create two new ramps + # This jump point corresponds to the index of the largest + # statistic: + # argmax(stats) + # These statistics are indexed relative to the + # ramp's range. Therefore, we need to add the start index + # of the ramp to the result. + # + # Note that because the resultants are averages of reads, but + # jumps occur in individual reads, it is possible that the + # jump is averaged down by the resultant with the actual jump + # causing the computed jump to be off by one index. + # In the idealized case this is when the jump occurs near + # the start of the resultant with the jump. In this case, + # the statistic for the resultant will be maximized at + # index - 1 rather than index. This means that we have to + # remove argmax(stats) + 1 as it is also a possible jump. + # This case is difficult to distinguish from the case where + # argmax(stats) does correspond to the jump resultant. + # Therefore, we just remove both possible resultants from + # consideration. + jump0 = np.argmax(stats) + ramp.start + jump1 = jump0 + 1 + ramp_fits.jumps.push_back(jump0) + ramp_fits.jumps.push_back(jump1) + + # The two resultant indicies need to be skipped, therefore + # the two + # possible new ramps are: + # RampIndex(ramp.start, jump0 - 1) + # RampIndex(jump1 + 1, ramp.end) + # This is because the RampIndex contains the index of the + # first and last resulants in the sub-ramp it describes. + # Note: The algorithm works via working over the sub-ramps + # backward in time. Therefore, since we are using a stack, + # we need to add the ramps in the time order they were + # observed in. This results in the last observation ramp + # being the top of the stack; meaning that, + # it will be the next ramp handeled. + + if jump0 > ramp.start: + # Note that when jump0 == ramp.start, we have detected a + # jump in the first resultant of the ramp. This means + # there is no sub-ramp before jump0. + # Also, note that this will produce bad results as + # the ramp indexing will go out of bounds. So it is + # important that we exclude it. + # Note that jump0 < ramp.start is not possible because + # the argmax is always >= 0 + ramps.push(RampIndex(ramp.start, jump0 - 1)) + + if jump1 < ramp.end: + # Note that if jump1 == ramp.end, we have detected a + # jump in the last resultant of the ramp. This means + # there is no sub-ramp after jump1. + # Also, note that this will produce bad results as + # the ramp indexing will go out of bounds. So it is + # important that we exclude it. + # Note that jump1 > ramp.end is technically possible + # however in those potential cases it will draw on + # resultants which are not considered part of the ramp + # under consideration. Therefore, we have to exlude all + # of those values. + ramps.push(RampIndex(jump1 + 1, ramp.end)) + + continue + + # Add ramp_fit to ramp_fits if no jump detection or stats are less + # than threshold + # Note that ramps are computed backward in time meaning we need to + # reverse the order of the fits at the end + ramp_fits.fits.push_back(ramp_fit) + ramp_fits.index.push_back(ramp) + + # Start computing the averages + # Note we do not do anything in the NaN case for degenerate ramps + if not np.isnan(ramp_fit.slope): + # protect weight against the extremely unlikely case of a zero + # variance + weight = 0 if ramp_fit.read_var == 0 else 1 / ramp_fit.read_var + total_weight += weight + + ramp_fits.average.slope += weight * ramp_fit.slope + ramp_fits.average.read_var += weight**2 * ramp_fit.read_var + ramp_fits.average.poisson_var += weight**2 * ramp_fit.poisson_var + + # Reverse to order in time + ramp_fits.fits = ramp_fits.fits[::-1] + ramp_fits.index = ramp_fits.index[::-1] + + # Finish computing averages + ramp_fits.average.slope /= total_weight if total_weight != 0 else 1 + ramp_fits.average.read_var /= total_weight**2 if total_weight != 0 else 1 + ramp_fits.average.poisson_var /= total_weight**2 if total_weight != 0 else 1 + + # Multiply poisson term by flux, (no negative fluxes) + ramp_fits.average.poisson_var *= max(ramp_fits.average.slope, 0) + + return ramp_fits + + def _to_dict(Pixel self): + """ + This is a private method to convert the Pixel object to a dictionary, so + that attributes can be directly accessed in python. Note that this is + needed because class attributes cannot be accessed on cython classes + directly in python. Instead they need to be accessed or set using a + python compatible method. This method is a pure puthon method bound + to to the cython class and should not be used by any cython code, and + only exists for testing purposes. + """ + + cdef np.ndarray[float, ndim=1] resultants_ = np.array(self.resultants, dtype=np.float32) + + cdef np.ndarray[float, ndim=2] local_slopes + cdef np.ndarray[float, ndim=2] var_read_noise + + if self.fixed.use_jump: + local_slopes = np.array(self.local_slopes, dtype=np.float32) + var_read_noise = np.array(self.var_read_noise, dtype=np.float32) + else: + try: + self.local_slopes + except AttributeError: + local_slopes = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("local_slopes should not exist") + + try: + self.var_read_noise + except AttributeError: + var_read_noise = np.array([[np.nan],[np.nan]], dtype=np.float32) + else: + raise AttributeError("var_read_noise should not exist") + + return dict(fixed=self.fixed._to_dict(), + resultants=resultants_, + read_noise=self.read_noise, + local_slopes=local_slopes, + var_read_noise=var_read_noise) + + +@cython.boundscheck(False) +@cython.wraparound(False) +cpdef inline Pixel make_pixel(FixedValues fixed, float read_noise, float [:] resultants): + """ + Fast constructor for the Pixel C class. + This creates a Pixel object for a single pixel from the input data. + + This is signifantly faster than using the `__init__` or `__cinit__` + this is because this does not have to pass through the Python as part + of the construction. + + Parameters + ---------- + fixed : FixedValues + Fixed values for all pixels + read_noise : float + read noise for the single pixel + resultants : float [:] + array of resultants for the single pixel + - memoryview of a numpy array to avoid passing through Python + + Return + ------ + Pixel C-class object (with pre-computed values if use_jump is True) + """ + cdef Pixel pixel = Pixel() + + # Fill in input information for pixel + pixel.fixed = fixed + pixel.read_noise = read_noise + pixel.resultants = resultants + + # Pre-compute values for jump detection shared by all pixels for this pixel + if fixed.use_jump: + pixel.local_slopes = pixel.local_slope_vals() + pixel.var_read_noise = read_noise * np.array(fixed.read_recip_coeffs) + + return pixel diff --git a/src/stcal/ramp_fitting/ols_cas22_fit.py b/src/stcal/ramp_fitting/ols_cas22_fit.py index daadea6f..9584970e 100644 --- a/src/stcal/ramp_fitting/ols_cas22_fit.py +++ b/src/stcal/ramp_fitting/ols_cas22_fit.py @@ -33,10 +33,19 @@ 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): +def fit_ramps_casertano( + resultants, + dq, + read_noise, + read_time, + read_pattern, + use_jump=False, + *, + threshold_intercept=None, + threshold_constant=None, +): """Fit ramps following Casertano+2022, including averaging partial ramps. Ramps are broken where dq != 0, and fits are performed on each sub-ramp. @@ -54,12 +63,18 @@ def fit_ramps_casertano(resultants, dq, read_noise, read_time, ma_table=None, re 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 + read_pattern : list[list[int]] The read pattern prescription. If None, use `ma_table`. One of `ma_table` or `read_pattern` must be defined. + use_jump : bool + If True, use the jump detection algorithm to identify CRs. + If False, use the DQ array to identify CRs. + threshold_intercept : float (optional, keyword-only) + Override the intercept parameter for threshold for the jump detection + algorithm. + theshold_constant : float (optional, keyword-only) + Override the constant parameter for threshold for the jump detection + algorithm. Returns ------- @@ -70,174 +85,50 @@ def fit_ramps_casertano(resultants, dq, read_noise, read_time, ma_table=None, re 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.') + # Trickery to avoid having to specify the defaults for the threshold + # parameters outside the cython code. + kwargs = {} + if threshold_intercept is not None: + kwargs['intercept'] = threshold_intercept + if threshold_constant is not None: + kwargs['constant'] = threshold_constant 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') + resultants = np.array(resultants).astype(np.float32) + dq = np.array(dq).astype(np.int32) if np.ndim(read_noise) <= 1: read_noise = read_noise * np.ones(resultants.shape[1:]) - read_noise = np.array(read_noise).astype('f4') + read_noise = np.array(read_noise).astype(np.float32) - origshape = resultants.shape + orig_shape = 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,)) + resultants = resultants.reshape(orig_shape + (1,)) + dq = dq.reshape(orig_shape + (1,)) + read_noise = read_noise.reshape(orig_shape[1:] + (1,)) - rampfitdict = ols_cas22.fit_ramps( + output = 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 + read_pattern, + use_jump, + **kwargs) + parameters = output.parameters.reshape(orig_shape[1:] + (2,)) + variances = output.variances.reshape(orig_shape[1:] + (3,)) + dq = output.dq.reshape(orig_shape) -def fit_ramps_casertano_no_dq(resultants, read_noise, ma_table): - """Fit ramps following Casertano+2022, only using full ramps. + if resultants.shape != orig_shape: + parameters = parameters[0] + variances = variances[0] - 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 + if resultants_unit is not None: + parameters = parameters * resultants_unit - 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 + return ols_cas22.RampFitOutputs(output.fits, parameters, variances, dq) diff --git a/src/stcal/ramp_fitting/ols_cas22_util.py b/src/stcal/ramp_fitting/ols_cas22_util.py deleted file mode 100644 index 63bb9a27..00000000 --- a/src/stcal/ramp_fitting/ols_cas22_util.py +++ /dev/null @@ -1,164 +0,0 @@ -"""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_jump_cas22.py b/tests/test_jump_cas22.py new file mode 100644 index 00000000..5d49597f --- /dev/null +++ b/tests/test_jump_cas22.py @@ -0,0 +1,613 @@ +import numpy as np +import pytest +from numpy.testing import assert_allclose + +from stcal.ramp_fitting.ols_cas22._core import metadata_from_read_pattern, threshold +from stcal.ramp_fitting.ols_cas22._fixed import fixed_values_from_metadata +from stcal.ramp_fitting.ols_cas22._pixel import make_pixel + +from stcal.ramp_fitting.ols_cas22 import fit_ramps, Parameter, Variance, Diff, RampJumpDQ + + +# Purposefully set a fixed seed so that the tests in this module are deterministic +RNG = np.random.default_rng(619) + +# The read time is constant for the given telescope/instrument so we set it here +# to be the one for Roman as it is known to be a reasonable value +READ_TIME = 3.04 + +# Choose small read noise relative to the flux to make it extremely unlikely +# that the random process will "accidentally" generate a set of data, which +# can trigger jump detection. This makes it easier to cleanly test jump +# detection is doing what we expect. +FLUX = 100 +READ_NOISE = np.float32(5) + +# Set a value for jumps which makes them obvious relative to the normal flux +JUMP_VALUE = 10_000 + +# Choose reasonable values for arbitrary test parameters, these are kept the same +# across all tests to make it easier to isolate the effects of something using +# multiple tests. +N_PIXELS = 100_000 +CHI2_TOL = 0.03 +GOOD_PROB = 0.7 + + +@pytest.fixture(scope="module") +def base_ramp_data(): + """ + Basic data for simulating ramps for testing (not unpacked) + + Returns + ------- + read_pattern : list[list[int]] + The example read pattern + metadata : dict + The metadata computed from the read pattern + """ + read_pattern = [ + [1, 2, 3, 4], + [5], + [6, 7, 8], + [9, 10, 11, 12, 13, 14, 15, 16, 17, 18], + [19, 20, 21], + [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36] + ] + + yield read_pattern, metadata_from_read_pattern(read_pattern, READ_TIME) + + +def test_metadata_from_read_pattern(base_ramp_data): + """Test turning read_pattern into the time data""" + _, data = base_ramp_data + + # Basic sanity checks (structs become dicts) + assert isinstance(data, dict) + assert 't_bar' in data + assert 'tau' in data + assert 'n_reads' in data + assert len(data) == 3 + + # Check that the data is correct + assert_allclose(data['t_bar'], [7.6, 15.2, 21.279999, 41.040001, 60.799999, 88.159996]) + assert_allclose(data['tau'], [5.7, 15.2, 19.928888, 36.023998, 59.448887, 80.593781]) + assert data['n_reads'] == [4, 1, 3, 10, 3, 15] + + +def test_init_ramps(): + """ + Test turning dq flags into initial ramp splits + Note that because `init_ramps` itself returns a stack, which does not have + a direct python equivalent, we call the wrapper for `init_ramps` which + converts that stack into a list ordered in the same fashion as the stack + """ + from stcal.ramp_fitting.ols_cas22._core import _init_ramps_list + + dq = np.array([[0, 1, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1], + [0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1], + [0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1], + [0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1]], dtype=np.int32) + + ramps = _init_ramps_list(dq) + assert len(ramps) == dq.shape[1] == 16 + + # Check that the ramps are correct + + # No DQ + assert ramps[0] == [{'start': 0, 'end': 3}] + + # 1 DQ + assert ramps[1] == [{'start': 1, 'end': 3}] + assert ramps[2] == [{'start': 0, 'end': 0}, {'start': 2, 'end': 3}] + assert ramps[3] == [{'start': 0, 'end': 1}, {'start': 3, 'end': 3}] + assert ramps[4] == [{'start': 0, 'end': 2}] + + # 2 DQ + assert ramps[5] == [{'start': 2, 'end': 3}] + assert ramps[6] == [{'start': 1, 'end': 1}, {'start': 3, 'end': 3}] + assert ramps[7] == [{'start': 1, 'end': 2}] + assert ramps[8] == [{'start': 0, 'end': 0}, {'start': 3, 'end': 3}] + assert ramps[9] == [{'start': 0, 'end': 0}, {'start': 2, 'end': 2}] + assert ramps[10] == [{'start': 0, 'end': 1}] + + # 3 DQ + assert ramps[11] == [{'start': 3, 'end': 3}] + assert ramps[12] == [{'start': 2, 'end': 2}] + assert ramps[13] == [{'start': 1, 'end': 1}] + assert ramps[14] == [{'start': 0, 'end': 0}] + + # 4 DQ + assert ramps[15] == [] + + +def test_threshold(): + """ + Test the threshold object/fucnction + intercept - constant * log10(slope) = threshold + """ + + # Create the python analog of the Threshold struct + # Note that structs get mapped to/from python as dictionary objects with + # the keys being the struct members. + thresh = { + 'intercept': np.float32(5.5), + 'constant': np.float32(1/3) + } + + # Check the 'intercept' is correctly interpreted. + # Since the log of the input slope is taken, log10(1) = 0, meaning that + # we should directly recover the intercept value in that case. + assert thresh['intercept'] == threshold(thresh, 1.0) + + # Check the 'constant' is correctly interpreted. + # Since we know that the intercept is correctly identified and that `log10(10) = 1`, + # we can use that to check that the constant is correctly interpreted. + assert np.float32(thresh['intercept'] - thresh['constant']) == threshold(thresh, 10.0) + + +@pytest.fixture(scope="module") +def ramp_data(base_ramp_data): + """ + Unpacked metadata for simulating ramps for testing + + Returns + ------- + read_pattern: + The read pattern used for testing + t_bar: + The t_bar values for the read pattern + tau: + The tau values for the read pattern + n_reads: + The number of reads for the read pattern + """ + read_pattern, read_pattern_metadata = base_ramp_data + t_bar = np.array(read_pattern_metadata['t_bar'], dtype=np.float32) + tau = np.array(read_pattern_metadata['tau'], dtype=np.float32) + n_reads = np.array(read_pattern_metadata['n_reads'], dtype=np.int32) + + yield read_pattern, t_bar, tau, n_reads + + +@pytest.mark.parametrize("use_jump", [True, False]) +def test_fixed_values_from_metadata(ramp_data, use_jump): + """Test computing the fixed data for all pixels""" + _, t_bar, tau, n_reads = ramp_data + + # Create the python analog of the ReadPatternMetadata struct + # Note that structs get mapped to/from python as dictionary objects with + # the keys being the struct members. + data = { + "t_bar": t_bar, + "tau": tau, + "n_reads": n_reads, + } + + # Create the python analog of the Threshold struct + # Note that structs get mapped to/from python as dictionary objects with + # the keys being the struct members. + thresh = { + 'intercept': np.float32(5.5), + 'constant': np.float32(1/3) + } + + # Note this is converted to a dictionary so we can directly interrogate the + # variables in question + fixed = fixed_values_from_metadata(data, thresh, use_jump)._to_dict() + + # Basic sanity checks that data passed in survives + assert (fixed['data']['t_bar'] == t_bar).all() + assert (fixed['data']['tau'] == tau).all() + assert (fixed['data']['n_reads'] == n_reads).all() + assert fixed['threshold']["intercept"] == thresh['intercept'] + assert fixed['threshold']["constant"] == thresh['constant'] + + # Check the computed data + # These are computed via vectorized operations in the main code, here we + # check using item-by-item operations + if use_jump: + single_gen = zip( + fixed['t_bar_diffs'][Diff.single], + fixed['t_bar_diff_sqrs'][Diff.single], + fixed['read_recip_coeffs'][Diff.single], + fixed['var_slope_coeffs'][Diff.single] + ) + double_gen = zip( + fixed['t_bar_diffs'][Diff.double], + fixed['t_bar_diff_sqrs'][Diff.double], + fixed['read_recip_coeffs'][Diff.double], + fixed['var_slope_coeffs'][Diff.double] + ) + + for index, (t_bar_diff_1, t_bar_diff_sqr_1, read_recip_1, var_slope_1) in enumerate(single_gen): + assert t_bar_diff_1 == t_bar[index + 1] - t_bar[index] + assert t_bar_diff_sqr_1 == np.float32((t_bar[index + 1] - t_bar[index]) ** 2) + assert read_recip_1 == np.float32(1 / n_reads[index + 1]) + np.float32(1 / n_reads[index]) + assert var_slope_1 == (tau[index + 1] + tau[index] - min(t_bar[index], t_bar[index + 1])) + + for index, (t_bar_diff_2, t_bar_diff_sqr_2, read_recip_2, var_slope_2) in enumerate(double_gen): + if index == len(fixed['t_bar_diffs'][1]) - 1: + # Last value must be NaN + assert np.isnan(t_bar_diff_2) + assert np.isnan(read_recip_2) + assert np.isnan(var_slope_2) + else: + assert t_bar_diff_2 == t_bar[index + 2] - t_bar[index] + assert t_bar_diff_sqr_2 == np.float32((t_bar[index + 2] - t_bar[index])**2) + assert read_recip_2 == np.float32(1 / n_reads[index + 2]) + np.float32(1 / n_reads[index]) + assert var_slope_2 == (tau[index + 2] + tau[index] - min(t_bar[index], t_bar[index + 2])) + else: + # If not using jumps, these values should not even exist. However, for wrapping + # purposes, they are checked to be non-existent and then set to NaN + assert np.isnan(fixed['t_bar_diffs']).all() + assert np.isnan(fixed['t_bar_diff_sqrs']).all() + assert np.isnan(fixed['read_recip_coeffs']).all() + assert np.isnan(fixed['var_slope_coeffs']).all() + + +def _generate_resultants(read_pattern, n_pixels=1): + """ + Generate a set of resultants for a pixel + + Parameters: + read_pattern : list[list[int]] + The read pattern to use + n_pixels: + The number of pixels to generate resultants for. Default is 1. + + Returns: + resultants + The resultants generated + """ + resultants = np.zeros((len(read_pattern), n_pixels), dtype=np.float32) + + # Use Poisson process to simulate the accumulation of the ramp + ramp_value = np.zeros(n_pixels, dtype=np.float32) # Last value of ramp + for index, reads in enumerate(read_pattern): + resultant_total = np.zeros(n_pixels, dtype=np.float32) # Total of all reads in this resultant + for _ in reads: + # Compute the next value of the ramp + # Using a Poisson process for the flux + ramp_value += RNG.poisson(FLUX * READ_TIME, size=n_pixels).astype(np.float32) + + # Add to running total for the resultant + resultant_total += ramp_value + + # Add read noise to the resultant + resultant_total += ( + RNG.standard_normal(size=n_pixels, dtype=np.float32) * READ_NOISE / np.sqrt(len(reads)) + ) + + # Record the average value for resultant (i.e., the average of the reads) + resultants[index] = (resultant_total / len(reads)).astype(np.float32) + + if n_pixels == 1: + resultants = resultants[:, 0] + + return resultants + + +@pytest.fixture(scope="module") +def pixel_data(ramp_data): + """ + Create data for a single pixel + + Returns: + resultants + Resultants for a single pixel + t_bar: + The t_bar values for the read pattern used for the resultants + 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 + """ + + read_pattern, t_bar, tau, n_reads = ramp_data + resultants = _generate_resultants(read_pattern) + + yield resultants, t_bar, tau, n_reads + + +@pytest.mark.parametrize("use_jump", [True, False]) +def test_make_pixel(pixel_data, use_jump): + """Test computing the initial pixel data""" + resultants, t_bar, tau, n_reads = pixel_data + + # Create a fixed object to pass into the constructor + # This requires setting up some structs as dictionaries + data = { + "t_bar": t_bar, + "tau": tau, + "n_reads": n_reads, + } + thresh = { + 'intercept': np.float32(5.5), + 'constant': np.float32(1/3) + } + fixed = fixed_values_from_metadata(data, thresh, use_jump) + + # Note this is converted to a dictionary so we can directly interrogate the + # variables in question + pixel = make_pixel(fixed, READ_NOISE, resultants)._to_dict() + + # Basic sanity checks that data passed in survives + assert (pixel['resultants'] == resultants).all() + assert READ_NOISE == pixel['read_noise'] + + # the "fixed" data is not checked as this is already done above + + # Check the computed data + # These are computed via vectorized operations in the main code, here we + # check using item-by-item operations + if use_jump: + single_gen = zip(pixel['local_slopes'][Diff.single], pixel['var_read_noise'][Diff.single]) + 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 var_read_noise_1 == READ_NOISE * ( + np.float32(1 / n_reads[index + 1]) + np.float32(1 / n_reads[index]) + ) + + for index, (local_slope_2, var_read_noise_2) in enumerate(double_gen): + if index == len(pixel['local_slopes'][1]) - 1: + # Last value must be NaN + assert np.isnan(local_slope_2) + assert np.isnan(var_read_noise_2) + else: + assert local_slope_2 == ( + (resultants[index + 2] - resultants[index]) / (t_bar[index + 2] - t_bar[index]) + ) + assert var_read_noise_2 == READ_NOISE * ( + np.float32(1 / n_reads[index + 2]) + np.float32(1 / n_reads[index]) + ) + else: + # If not using jumps, these values should not even exist. However, for wrapping + # purposes, they are checked to be non-existent and then set to NaN + assert np.isnan(pixel['local_slopes']).all() + assert np.isnan(pixel['var_read_noise']).all() + + +@pytest.fixture(scope="module") +def detector_data(ramp_data): + """ + Generate a set of with no jumps data as if for a single detector as it + would be passed in by the supporting code. + + 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 + """ + read_pattern, *_ = ramp_data + read_noise = np.ones(N_PIXELS, dtype=np.float32) * READ_NOISE + + resultants = _generate_resultants(read_pattern, n_pixels=N_PIXELS) + + return resultants, read_noise, read_pattern + + +@pytest.mark.parametrize("use_jump", [True, False]) +@pytest.mark.parametrize("use_dq", [True, False]) +def test_fit_ramps(detector_data, use_jump, use_dq): + """ + Test fitting ramps + Since no jumps are simulated in the data, jump detection shouldn't pick + up any jumps. + """ + resultants, read_noise, read_pattern = detector_data + dq = ( + (RNG.uniform(size=resultants.shape) > GOOD_PROB).astype(np.int32) if use_dq else + np.zeros(resultants.shape, dtype=np.int32) + ) + + # 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. + okay = np.sum((dq[1:, :] == 0) & (dq[:-1, :] == 0), axis=0) != 0 + assert okay.dtype == bool + + # Note that for use_dq = False, okay == True for all ramps, so we perform + # a sanity check that the above criterion is correct + if not use_dq: + assert okay.all() + + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=use_jump) + assert len(output.fits) == N_PIXELS # sanity check that a fit is output for each pixel + + chi2 = 0 + for fit, use in zip(output.fits, okay): + if not use_dq: + # Check that the data generated does not generate any false positives + # for jumps as this data is reused for `test_find_jumps` below. + # This guarantees that all jumps detected in that test are the + # purposefully placed ones which we know about. So the `test_find_jumps` + # can focus on checking that the jumps found are the correct ones, + # and that all jumps introduced are detected properly. + assert len(fit['fits']) == 1 + + if use: + # Add okay ramps to chi2 + total_var = fit['average']['read_var'] + fit['average']['poisson_var'] + chi2 += (fit['average']['slope'] - FLUX)**2 / total_var + else: + # Check no slope fit for bad ramps + assert fit['average']['slope'] == 0 + assert fit['average']['read_var'] == 0 + assert fit['average']['poisson_var'] == 0 + + assert use_dq # sanity check that this branch is only encountered when use_dq = True + + chi2 /= np.sum(okay) + assert np.abs(chi2 - 1) < CHI2_TOL + + +@pytest.mark.parametrize("use_jump", [True, False]) +def test_fit_ramps_array_outputs(detector_data, use_jump): + """ + Test that the array outputs line up with the dictionary output + """ + resultants, read_noise, read_pattern = detector_data + dq = np.zeros(resultants.shape, dtype=np.int32) + + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=use_jump) + + for fit, par, var in zip(output.fits, output.parameters, output.variances): + assert par[Parameter.intercept] == 0 + assert par[Parameter.slope] == fit['average']['slope'] + + assert var[Variance.read_var] == fit['average']['read_var'] + assert var[Variance.poisson_var] == fit['average']['poisson_var'] + assert var[Variance.total_var] == np.float32( + fit['average']['read_var'] + fit['average']['poisson_var'] + ) + + +@pytest.fixture(scope="module") +def jump_data(detector_data): + """ + Generate resultants with single jumps in them for testing jump detection. + 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 + + # Choose read to place a single jump in for each pixel + num_reads = read_pattern[-1][-1] + jump_reads = RNG.integers(num_reads - 1, size=N_PIXELS) + + # This shows that a jump as been placed in every single possible + # read position. Technically, this check can fail; however, + # N_PIXELS >> num_reads so it is very unlikely in practice since + # all reads are equally likely to be chosen for a jump. + # It is a good check that we can detect a jump occurring in any read except + # the first read. + assert set(jump_reads) == set(range(num_reads - 1)) + + # Fill out jump reads with jump values + jump_flux = np.zeros((num_reads, N_PIXELS), dtype=np.float32) + for index, jump in enumerate(jump_reads): + jump_flux[jump:, index] = JUMP_VALUE + + # Average the reads into the resultants + jump_resultants = np.zeros(N_PIXELS, dtype=np.int32) + for index, reads in enumerate(read_pattern): + indices = np.array(reads) - 1 + resultants[index, :] += np.mean(jump_flux[indices, :], axis=0) + for read in reads: + jump_resultants[np.where(jump_reads == read)] = index + + return resultants, read_noise, read_pattern, jump_reads, jump_resultants + + +def test_find_jumps(jump_data): + """ + Full unit tests to demonstrate that we can detect jumps in any read (except + the first one) and that we correctly remove these reads from the fit to recover + the correct FLUX/slope. + """ + resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data + dq = np.zeros(resultants.shape, dtype=np.int32) + + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) + assert len(output.fits) == len(jump_reads) # sanity check that a fit/jump is set for every pixel + + chi2 = 0 + for fit, jump_index, resultant_index in zip(output.fits, jump_reads, jump_resultants): + + # Check that the jumps are detected correctly + if jump_index == 0: + # There is no way to detect a jump if it is in the very first read + # The very first pixel in this case has a jump in the first read + assert len(fit['jumps']) == 0 + assert resultant_index == 0 # sanity check that the jump is indeed in the first resultant + + # Test the correct ramp_index was recorded: + assert len(fit['index']) == 1 + assert fit['index'][0]['start'] == 0 + assert fit['index'][0]['end'] == len(read_pattern) - 1 + else: + # There should be a single jump detected; however, this results in + # two resultants being excluded. + assert len(fit['jumps']) == 2 + assert resultant_index in fit['jumps'] + + # The two resultants excluded should be adjacent + for jump in fit['jumps']: + assert jump == resultant_index or jump == resultant_index - 1 or jump == resultant_index + 1 + + # Test the correct ramp indexes are recorded + 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_indices = list(range(ramp_index["start"], ramp_index["end"] + 1)) + + # check that all the ramps are non-overlapping + assert set(ramp_indices).isdisjoint(new_indices) + + ramp_indices.extend(new_indices) + + # check that no ramp_index is a jump + assert set(ramp_indices).isdisjoint(fit['jumps']) + + # check that all resultant indices are either in a ramp or listed as a jump + assert set(ramp_indices).union(fit['jumps']) == set(range(len(read_pattern))) + + # Compute the chi2 for the fit and add it to a running "total chi2" + total_var = fit['average']['read_var'] + fit['average']['poisson_var'] + chi2 += (fit['average']['slope'] - FLUX)**2 / total_var + + # Check that the average chi2 is ~1. + chi2 /= N_PIXELS + assert np.abs(chi2 - 1) < CHI2_TOL + + +def test_override_default_threshold(jump_data): + """This tests that we can override the default jump detection threshold constants""" + resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data + dq = np.zeros(resultants.shape, dtype=np.int32) + + standard = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) + override = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True, + intercept=0, constant=0) + + # All this is intended to do is show that with all other things being equal passing non-default + # threshold parameters changes the results. + assert (standard.parameters != override.parameters).any() + + +def test_jump_dq_set(jump_data): + # Check the DQ flag value to start + assert RampJumpDQ.JUMP_DET == 2**2 + + resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data + dq = np.zeros(resultants.shape, dtype=np.int32) + + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) + + for fit, pixel_dq in zip(output.fits, output.dq.transpose()): + # Check that all jumps found get marked + assert (pixel_dq[fit['jumps']] == RampJumpDQ.JUMP_DET).all() + + # Check that dq flags for jumps are only set if the jump is marked + assert set(np.where(pixel_dq == RampJumpDQ.JUMP_DET)[0]) == set(fit['jumps']) diff --git a/tests/test_ramp_fitting_cas22.py b/tests/test_ramp_fitting_cas22.py index 4eb63aa7..72a7b6d1 100644 --- a/tests/test_ramp_fitting_cas22.py +++ b/tests/test_ramp_fitting_cas22.py @@ -2,10 +2,11 @@ """ Unit tests for ramp-fitting functions. """ +import astropy.units as u import numpy as np +import pytest 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 @@ -14,54 +15,75 @@ 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]] +@pytest.mark.parametrize("use_unit", [True, False]) +@pytest.mark.parametrize("use_dq", [True, False]) +def test_simulated_ramps(use_unit, use_dq): + # Perfect square like the detector, this is so we can test that the code + # reshapes the data correctly for the computation and then reshapes it back + # to the original shape. + ntrial = 320 * 320 + read_pattern, flux, read_noise, resultants = simulate_many_ramps(ntrial=ntrial) + + # So we get a detector-like input shape + resultants = resultants.reshape((len(read_pattern), 320, 320)) + + if use_unit: + resultants = resultants * u.electron + + dq = np.zeros(resultants.shape, dtype=np.int32) + read_noise = np.ones(resultants.shape[1:], dtype=np.float32) * read_noise + + # now let's mark a bunch of the ramps as compromised. When using dq flags + if use_dq: + bad = np.random.uniform(size=resultants.shape) > 0.7 + dq |= bad + + output = ramp.fit_ramps_casertano( + resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, + threshold_constant=0, threshold_intercept=0) # set the threshold parameters + # to demo the interface. This + # will raise an error if + # the interface changes, but + # does not effect the computation + # since jump detection is off in + # this case. + + # Check that the output shapes are correct + assert output.parameters.shape == (320, 320, 2) == resultants.shape[1:] + (2,) + assert output.variances.shape == (320, 320, 3) == resultants.shape[1:] + (3,) + assert output.dq.shape == dq.shape + + # check the unit + if use_unit: + assert output.parameters.unit == u.electron + parameters = output.parameters.value + else: + parameters = output.parameters + + # Turn into single dimension arrays to make the indexing for the math easier + parameters = parameters.reshape((320 * 320, 2)) + variances = output.variances.reshape((320 * 320, 3)) - 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) + okay = np.sum((dq[1:, :] == 0) & (dq[:-1, :] == 0), axis=0) != 0 + okay = okay.reshape((320 * 320)) + + # Sanity check that when no dq is used, all ramps are used + if not use_dq: + assert np.all(okay) + + chi2dof_slope = np.sum((parameters[okay, 1] - flux)**2 / variances[okay, 2]) / np.sum(okay) assert np.abs(chi2dof_slope - 1) < 0.03 - assert np.all(par[~m, 1] == 0) - assert np.all(var[~m, 1] == 0) + assert np.all(parameters[~okay, 1] == 0) + assert np.all(variances[~okay, 1] == 0) # ######### # Utilities # ######### -def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, ma_table=None): +def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, read_pattern=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 @@ -75,9 +97,8 @@ def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, ma_table=None): 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 + read_pattern : list[list] (int) + An optional read pattern Returns ------- @@ -90,18 +111,22 @@ def simulate_many_ramps(ntrial=100, flux=100, readnoise=5, ma_table=None): 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') + if read_pattern is None: + read_pattern = [[1, 2, 3, 4], + [5], + [6, 7, 8], + [9, 10, 11, 12, 13, 14, 15, 16, 17, 18], + [19, 20, 21], + [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]] + nread = np.array([len(x) for x in read_pattern]) + resultants = np.zeros((len(read_pattern), ntrial), dtype='f4') buf = np.zeros(ntrial, dtype='i4') - for i, ti in enumerate(tij): + for i, reads in enumerate(read_pattern): subbuf = np.zeros(ntrial, dtype='i4') - for t0 in ti: + for _ in reads: 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) + resultants[i] = (subbuf / len(reads)).astype('f4') + resultants += np.random.randn(len(read_pattern), ntrial) * ( + readnoise / np.sqrt(nread)).reshape(len(read_pattern), 1) + return (read_pattern, flux, readnoise, resultants)