diff --git a/CHANGES.rst b/CHANGES.rst index 14bf31e0..03572974 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -15,6 +15,8 @@ ramp_fitting - Moving some CI tests from JWST to STCAL. [#228, spacetelescope/jwst#6080] +- Significantly, improve the performance of the uneven ramp fitting algorithm. [#229] + Changes to API -------------- diff --git a/setup.py b/setup.py index b34e7dfa..ffb85902 100644 --- a/setup.py +++ b/setup.py @@ -8,26 +8,20 @@ extensions = [ Extension( - 'stcal.ramp_fitting.ols_cas22._core', - ['src/stcal/ramp_fitting/ols_cas22/_core.pyx'], + 'stcal.ramp_fitting.ols_cas22._ramp', + ['src/stcal/ramp_fitting/ols_cas22/_ramp.pyx'], include_dirs=[np.get_include()], language='c++' ), Extension( - 'stcal.ramp_fitting.ols_cas22._fixed', - ['src/stcal/ramp_fitting/ols_cas22/_fixed.pyx'], + 'stcal.ramp_fitting.ols_cas22._jump', + ['src/stcal/ramp_fitting/ols_cas22/_jump.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'], + 'stcal.ramp_fitting.ols_cas22._fit', + ['src/stcal/ramp_fitting/ols_cas22/_fit.pyx'], include_dirs=[np.get_include()], language='c++' ), diff --git a/src/stcal/ramp_fitting/ols_cas22/__init__.py b/src/stcal/ramp_fitting/ols_cas22/__init__.py index d07f4ddf..9d0e18e6 100644 --- a/src/stcal/ramp_fitting/ols_cas22/__init__.py +++ b/src/stcal/ramp_fitting/ols_cas22/__init__.py @@ -1,4 +1,4 @@ -from ._fit_ramps import fit_ramps, RampFitOutputs -from ._core import Parameter, Variance, Diff, RampJumpDQ +from ._fit import fit_ramps, RampFitOutputs, Parameter, Variance +from ._jump import JUMP_DET -__all__ = ['fit_ramps', 'RampFitOutputs', 'Parameter', 'Variance', 'Diff', 'RampJumpDQ'] +__all__ = ['fit_ramps', 'RampFitOutputs', 'Parameter', 'Variance', 'Diff', 'JUMP_DET'] diff --git a/src/stcal/ramp_fitting/ols_cas22/_core.pxd b/src/stcal/ramp_fitting/ols_cas22/_core.pxd deleted file mode 100644 index 641c206e..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_core.pxd +++ /dev/null @@ -1,56 +0,0 @@ -from libcpp.vector cimport vector - - -cdef struct RampIndex: - int start - int end - - -cdef struct RampFit: - float slope - float read_var - float poisson_var - - -cdef struct RampFits: - RampFit average - vector[int] jumps - vector[RampFit] fits - vector[RampIndex] index - - -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) -cpdef vector[RampIndex] init_ramps(int[:, :] dq, int n_resultants, int index_pixel) -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 deleted file mode 100644 index 0d312261..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_core.pyx +++ /dev/null @@ -1,235 +0,0 @@ -""" -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.vector cimport vector -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) - """ - # clip slope in 1, 1e4 - slope = slope if slope > 1 else 1 - slope = slope if slope < 1e4 else 1e4 - return thresh.intercept - thresh.constant * log10(slope) - - -@cython.boundscheck(False) -@cython.wraparound(False) -cpdef inline vector[RampIndex] init_ramps(int[:, :] dq, int n_resultants, int index_pixel): - """ - 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 - n_resultants : int - Number of resultants - index_pixel : int - The index of the pixel to find ramps for - - Returns - ------- - vector of RampIndex objects - - vector with entry for each ramp found (last entry is last ramp found) - - RampIndex with start and end indices of the ramp in the resultants - """ - cdef vector[RampIndex] ramps = vector[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 - cdef RampIndex 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_back(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_back(ramp) - - return ramps - - -@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.pyx b/src/stcal/ramp_fitting/ols_cas22/_fit.pyx new file mode 100644 index 00000000..87aaf02d --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_fit.pyx @@ -0,0 +1,235 @@ +# cython: language_level=3str + +""" +External interface module for the Casertano+22 ramp fitting algorithm with jump detection. + This module is intended to contain everything needed by external code. + +Enums +----- +Parameter : + Enumerate the index for the output parameters array. + +Variance : + Enumerate the index for the output variances array. + +Classes +------- +RampFitOutputs : NamedTuple + Simple tuple wrapper for outputs from the ramp fitting algorithm + This clarifies the meaning of the outputs via naming them something + descriptive. + +(Public) Functions +------------------ +fit_ramps : function + Fit ramps using the Castenario+22 algorithm to a set of pixels accounting + for jumps (if use_jump is True) and bad pixels (via the dq array). This + is the primary externally callable function. +""" +from __future__ import annotations + +import numpy as np +cimport numpy as cnp + +from cython cimport boundscheck, wraparound +from libcpp cimport bool +from libcpp.list cimport list as cpp_list + +from stcal.ramp_fitting.ols_cas22._jump cimport (Thresh, + JumpFits, + fill_fixed_values, + fit_jumps, + n_fixed_offsets, + n_pixel_offsets) +from stcal.ramp_fitting.ols_cas22._ramp cimport ReadPattern, from_read_pattern + +from typing import NamedTuple + + +# Initialize numpy for cython use in this module +cnp.import_array() + + +cpdef enum Parameter: + intercept + slope + n_param + + +cpdef enum Variance: + read_var + poisson_var + total_var + n_var + + +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 + ---------- + 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 of RampFits + the raw ramp fit outputs, these are all structs which will get mapped to + python dictionaries. + """ + parameters: np.ndarray + variances: np.ndarray + dq: np.ndarray + fits: list | None = None + + +@boundscheck(False) +@wraparound(False) +def fit_ramps(float[:, :] resultants, + cnp.ndarray[int, ndim=2] dq, + float[:] read_noise, + float read_time, + list[list[int]] read_pattern, + bool use_jump=False, + float intercept=5.5, + float constant=1/3, + bool include_diagnostic=False): + """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 : float[n_resultants, n_pixel] + the resultants in electrons (Note that this can be based as any sort of + array, such as a numpy array. The memmory view is just for efficiency in + cython) + dq : np.ndarry[n_resultants, n_pixel] + the dq array. dq != 0 implies bad pixel / CR. (Kept as a numpy array + so that it can be passed out without copying into new numpy array, will + be working on memmory views of this array) + read_noise : float[n_pixel] + the read noise in electrons for each pixel (same note as the resultants) + 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 + include_diagnostic : bool + If True, include the raw ramp fits in the output. Default=False + + Returns + ------- + A RampFitOutputs tuple + """ + cdef int n_pixels, n_resultants + n_resultants = resultants.shape[0] + n_pixels = resultants.shape[1] + + # Raise error if input data is inconsistent + 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}') + + # Compute the main metadata from the read pattern and cast it to memory views + cdef ReadPattern metadata = from_read_pattern(read_pattern, read_time, n_resultants) + cdef float[:] t_bar = metadata.t_bar + cdef float[:] tau = metadata.tau + cdef int[:] n_reads = metadata.n_reads + + # Setup pre-compute arrays for jump detection + cdef float[:, :] fixed + cdef float[:, :] pixel + if use_jump: + # Initialize arrays for the jump detection pre-computed values + fixed = np.empty((n_fixed_offsets, n_resultants - 1), dtype=np.float32) + pixel = np.empty((n_pixel_offsets, n_resultants - 1), dtype=np.float32) + + # Pre-compute the values from the read pattern + fixed = fill_fixed_values(fixed, t_bar, tau, n_reads, n_resultants) + else: + # "Initialize" the arrays when not using jump detection, they need to be + # initialized because they do get passed around, but they don't need + # to actually have any entries + fixed = np.empty((0, 0), dtype=np.float32) + pixel = np.empty((0, 0), dtype=np.float32) + + # Create a threshold struct + cdef Thresh thresh = Thresh(intercept, constant) + + # Create variable to old the diagnostic data + # 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[JumpFits] ramp_fits + + # Initialize the output arrays. Note that the fit intercept is currently always + # zero, where as every variance is calculated and set. This means that the + # parameters need to be filled with zeros, where as the variances can just + # be allocated + cdef float[:, :] parameters = np.zeros((n_pixels, Parameter.n_param), dtype=np.float32) + cdef float[:, :] variances = np.empty((n_pixels, Variance.n_var), dtype=np.float32) + + # Cast the enum values into integers for indexing (otherwise compiler complains) + # These will be optimized out + cdef int slope = Parameter.slope + cdef int read_var = Variance.read_var + cdef int poisson_var = Variance.poisson_var + cdef int total_var = Variance.total_var + + # Pull memory view of dq for speed of access later + # changes to this array will backpropagate to the original numpy array + cdef int[:, :] dq_ = dq + + # Run the jump fitting algorithm for each pixel + cdef JumpFits fit + cdef int index + for index in range(n_pixels): + # Fit all the ramps for the given pixel + fit = fit_jumps(resultants[:, index], + dq_[:, index], + read_noise[index], + t_bar, + tau, + n_reads, + n_resultants, + fixed, + pixel, + thresh, + use_jump, + include_diagnostic) + + # Extract the output fit's parameters + parameters[index, slope] = fit.average.slope + + # Extract the output fit's variances + variances[index, read_var] = fit.average.read_var + variances[index, poisson_var] = fit.average.poisson_var + variances[index, total_var] = fit.average.read_var + fit.average.poisson_var + + # Store diagnostic data if requested + if include_diagnostic: + ramp_fits.push_back(fit) + + # Cast memory views into numpy arrays for ease of use in python. + return RampFitOutputs(np.array(parameters, dtype=np.float32), + np.array(variances, dtype=np.float32), + dq, + ramp_fits if include_diagnostic else None) diff --git a/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx b/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx deleted file mode 100644 index 4b2e39e3..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_fit_ramps.pyx +++ /dev/null @@ -1,135 +0,0 @@ -import numpy as np -cimport numpy as np -from libcpp cimport bool -from libcpp.list cimport list as cpp_list -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, Optional - - -# 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 - ---------- - 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 of RampFits - the raw ramp fit outputs, these are all structs which will get mapped to - python dictionaries. - """ - parameters: np.ndarray - variances: np.ndarray - dq: np.ndarray - fits: Optional[list] = None - - -@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, - bool include_diagnostic=False): - """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 - include_diagnostic : bool - If True, include the raw ramp fits in the output. Default=False - - 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) - - # 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(init_ramps(dq, n_resultants, index), include_diagnostic) - - 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 - - if include_diagnostic: - ramp_fits.push_back(fit) - - # return RampFitOutputs(ramp_fits, parameters, variances, dq) - return RampFitOutputs(parameters, variances, dq, ramp_fits if include_diagnostic else None) diff --git a/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd b/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd deleted file mode 100644 index 72087051..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_fixed.pxd +++ /dev/null @@ -1,21 +0,0 @@ -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 deleted file mode 100644 index 8417507b..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx +++ /dev/null @@ -1,274 +0,0 @@ -""" -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 libcpp cimport bool - -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] - - 2 * min(t_bar[i], t_bar[i+1])) - double of slope variance term: - var_slope_coeffs[Diff.double, :] = (tau[i] + tau[i+2] - - 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]) - 2 * np.minimum(t_bar[1:], t_bar[:end - 1])) - var_slope_vals[Diff.double, :end - 2] = (np.add(tau[2:], tau[:end - 2]) - 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/_jump.pxd b/src/stcal/ramp_fitting/ols_cas22/_jump.pxd new file mode 100644 index 00000000..8693e791 --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_jump.pxd @@ -0,0 +1,62 @@ +# cython: language_level=3str + +from libcpp cimport bool +from libcpp.vector cimport vector + +from stcal.ramp_fitting.ols_cas22._ramp cimport RampFit, RampQueue + + +cpdef enum FixedOffsets: + single_t_bar_diff + double_t_bar_diff + single_t_bar_diff_sqr + double_t_bar_diff_sqr + single_read_recip + double_read_recip + single_var_slope_val + double_var_slope_val + n_fixed_offsets + + +cpdef enum PixelOffsets: + single_local_slope + double_local_slope + single_var_read_noise + double_var_read_noise + n_pixel_offsets + + +cpdef enum: + JUMP_DET = 4 + +cdef struct Thresh: + float intercept + float constant + + +cdef struct JumpFits: + RampFit average + vector[int] jumps + vector[RampFit] fits + RampQueue index + + +cpdef float[:, :] fill_fixed_values(float[:, :] fixed, + float[:] t_bar, + float[:] tau, + int[:] n_reads, + int n_resultants) + + +cdef JumpFits fit_jumps(float[:] resultants, + int[:] dq, + float read_noise, + float[:] t_bar, + float[:] tau, + int[:] n_reads, + int n_resultants, + float[:, :] fixed, + float[:, :] pixel, + Thresh thresh, + bool use_jump, + bool include_diagnostic) diff --git a/src/stcal/ramp_fitting/ols_cas22/_jump.pyx b/src/stcal/ramp_fitting/ols_cas22/_jump.pyx new file mode 100644 index 00000000..4ace423b --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_jump.pyx @@ -0,0 +1,595 @@ +# cython: language_level=3str + + +""" +This module contains all the functions needed to execute jump detection for the + Castentano+22 ramp fitting algorithm + + The _ramp module contains the actual ramp fitting algorithm, this module + contains a driver for the algoritm and detection of jumps/splitting ramps. + +Structs +------- +Thresh : struct + intercept - constant * log10(slope) + - intercept : float + The intercept of the jump threshold + - constant : float + The constant of the jump threshold + +JumpFits : struct + All the data on a given pixel's ramp fit with (or without) jump detection + - average : RampFit + The average of all the ramps fit for the pixel + - jumps : vector[int] + The indices of the resultants which were detected as jumps + - fits : vector[RampFit] + All of the fits for each ramp fit for the pixel + - index : RampQueue + The RampIndex representations correspoinding to each fit in fits + +Enums +----- +FixedOffsets : enum + Enumerate the different pieces of information computed for jump detection + which only depend on the read pattern. + +PixelOffsets : enum + Enumerate the different pieces of information computed for jump detection + which only depend on the given pixel (independent of specific ramp). + +JUMP_DET : value + A the fixed value for the jump detection dq flag. + +(Public) Functions +------------------ +fill_fixed_values : function + Pre-compute all the values needed for jump detection for a given read_pattern, + this is independent of the pixel involved. + +fit_jumps : function + Compute all the ramps for a single pixel using the Casertano+22 algorithm + with jump detection. This is a driver for the ramp fit algorithm in general + meaning it automatically handles splitting ramps across dq flags in addition + to splitting across detected jumps (if jump detection is turned on). +""" + +from cython cimport boundscheck, wraparound, cdivision + +from libcpp cimport bool +from libc.math cimport sqrt, log10, fmaxf, NAN, isnan + +from stcal.ramp_fitting.ols_cas22._jump cimport Thresh, JumpFits, JUMP_DET, FixedOffsets, PixelOffsets +from stcal.ramp_fitting.ols_cas22._ramp cimport RampIndex, RampQueue, RampFit, fit_ramp, init_ramps + + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cpdef inline float[:, :] fill_fixed_values(float[:, :] fixed, + float[:] t_bar, + float[:] tau, + int[:] n_reads, + int n_resultants): + """ + Pre-compute all the values needed for jump detection which only depend on + the read pattern. + + Parameters + ---------- + fixed : float[:, :] + A pre-allocated memoryview to store the pre-computed values in, its faster + to allocate outside this function. + t_bar : float[:] + The average time for each resultant + tau : float[:] + The time variance for each resultant + n_reads : int[:] + The number of reads for each resultant + n_resultants : int + The number of resultants for the read pattern + + Returns + ------- + [ + , + , + ** 2, + ** 2, + <(1/n_reads[i+1] + 1/n_reads[i])>, + <(1/n_reads[i+2] + 1/n_reads[i])>, + <(tau[i] + tau[i+1] - 2 * min(t_bar[i], t_bar[i+1]))>, + <(tau[i] + tau[i+2] - 2 * min(t_bar[i], t_bar[i+2]))>, + ] + """ + # Cast the enum values into integers for indexing (otherwise compiler complains) + # These will be optimized out + cdef int single_t_bar_diff = FixedOffsets.single_t_bar_diff + cdef int double_t_bar_diff = FixedOffsets.double_t_bar_diff + cdef int single_t_bar_diff_sqr = FixedOffsets.single_t_bar_diff_sqr + cdef int double_t_bar_diff_sqr = FixedOffsets.double_t_bar_diff_sqr + cdef int single_read_recip = FixedOffsets.single_read_recip + cdef int double_read_recip = FixedOffsets.double_read_recip + cdef int single_var_slope_val = FixedOffsets.single_var_slope_val + cdef int double_var_slope_val = FixedOffsets.double_var_slope_val + + # Coerce division to be using floats + cdef float num = 1 + + cdef int i + for i in range(n_resultants - 1): + fixed[single_t_bar_diff, i] = t_bar[i + 1] - t_bar[i] + fixed[single_t_bar_diff_sqr, i] = fixed[single_t_bar_diff, i] ** 2 + fixed[single_read_recip, i] = (num / n_reads[i + 1]) + (num / n_reads[i]) + fixed[single_var_slope_val, i] = tau[i + 1] + tau[i] - 2 * min(t_bar[i + 1], t_bar[i]) + + if i < n_resultants - 2: + fixed[double_t_bar_diff, i] = t_bar[i + 2] - t_bar[i] + fixed[double_t_bar_diff_sqr, i] = fixed[double_t_bar_diff, i] ** 2 + fixed[double_read_recip, i] = (num / n_reads[i + 2]) + (num / n_reads[i]) + fixed[double_var_slope_val, i] = tau[i + 2] + tau[i] - 2 * min(t_bar[i + 2], t_bar[i]) + else: + # Last double difference is undefined + fixed[double_t_bar_diff, i] = NAN + fixed[double_t_bar_diff_sqr, i] = NAN + fixed[double_read_recip, i] = NAN + fixed[double_var_slope_val, i] = NAN + + return fixed + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cpdef inline float[:, :] _fill_pixel_values(float[:, :] pixel, + float[:] resultants, + float[:, :] fixed, + float read_noise, + int n_resultants): + """ + Pre-compute all the values needed for jump detection which only depend on + the a specific pixel (independent of the given ramp for a pixel). + + Parameters + ---------- + pixel : float[:, :] + A pre-allocated memoryview to store the pre-computed values in, its faster + to allocate outside this function. + resultants : float[:] + The resultants for the pixel in question. + fixed : float[:, :] + The pre-computed fixed values for the read_pattern + read_noise : float + The read noise for the pixel + n_resultants : int + The number of resultants for the read_pattern + + 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])>, + read_noise**2 * <(1/n_reads[i+1] + 1/n_reads[i])>, + read_noise**2 * <(1/n_reads[i+2] + 1/n_reads[i])>, + ] + """ + cdef int single_t_bar_diff = FixedOffsets.single_t_bar_diff + cdef int double_t_bar_diff = FixedOffsets.double_t_bar_diff + cdef int single_read_recip = FixedOffsets.single_read_recip + cdef int double_read_recip = FixedOffsets.double_read_recip + + cdef int single_slope = PixelOffsets.single_local_slope + cdef int double_slope = PixelOffsets.double_local_slope + cdef int single_var = PixelOffsets.single_var_read_noise + cdef int double_var = PixelOffsets.double_var_read_noise + + cdef float read_noise_sqr = read_noise ** 2 + + cdef int i + for i in range(n_resultants - 1): + pixel[single_slope, i] = (resultants[i + 1] - resultants[i]) / fixed[single_t_bar_diff, i] + pixel[single_var, i] = read_noise_sqr * fixed[single_read_recip, i] + + if i < n_resultants - 2: + pixel[double_slope, i] = (resultants[i + 2] - resultants[i]) / fixed[double_t_bar_diff, i] + pixel[double_var, i] = read_noise_sqr * fixed[double_read_recip, i] + else: + # The last double difference is undefined + pixel[double_slope, i] = NAN + pixel[double_var, i] = NAN + + return pixel + + +cdef 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) + """ + slope = slope if slope > 1 else 1 + slope = slope if slope < 1e4 else 1e4 + + return thresh.intercept - thresh.constant * log10(slope) + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cdef inline float _correction(float[:] t_bar, RampIndex ramp, float slope): + """ + Compute the correction factor for the variance used by a statistic + + - slope / (t_bar[end] - t_bar[start]) + + Parameters + ---------- + t_bar : float[:] + The computed t_bar values for the ramp + ramp : RampIndex + Struct for start and end indices resultants for the ramp + slope : float + The computed slope for the ramp + """ + + cdef float diff = t_bar[ramp.end] - t_bar[ramp.start] + + return - slope / diff + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cdef inline float _statstic(float local_slope, + float var_read_noise, + float t_bar_diff_sqr, + float var_slope_coeff, + float slope, + float correct): + """ + Compute a single fit statistic + delta / sqrt(var + correct) + + where: + delta = local_slope - slope + var = (var_read_noise + slope * var_slope_coeff) / t_bar_diff_sqr + + pre-computed: + local_slope = (resultant[i + j] - resultant[i]) / (t_bar[i + j] - t_bar[i]) + var_read_noise = read_noise ** 2 * (1/n_reads[i + j] + 1/n_reads[i]) + var_slope_coeff = tau[i + j] + tau[i] - 2 * min(t_bar[i + j], t_bar[i]) + t_bar_diff_sqr = (t_bar[i + j] - t_bar[i]) ** 2 + + Parameters + ---------- + local_slope : float + The local slope the statistic is computed for + float : var_read_noise + The read noise variance for local_slope + t_bar_diff_sqr : float + The square difference for the t_bar corresponding to local_slope + var_slope_coeff : float + The slope variance coefficient for local_slope + slope : float + The computed slope for the ramp + correct : float + The correction factor needed + + Returns + ------- + Create a single instance of the stastic for the given parameters + """ + + cdef float delta = local_slope - slope + cdef float var = (var_read_noise + slope * var_slope_coeff) / t_bar_diff_sqr + + return delta / sqrt(var + correct) + + +@boundscheck(False) +@wraparound(False) +cdef inline (int, float) _fit_statistic(float[:, :] pixel, + float[:, :] fixed, + float[:] t_bar, + float slope, + RampIndex ramp): + """ + Compute the maximum index and its value over all fit statistics for a given + ramp. Each index's stat is the max of the single and double difference + statistics: + all_stats = + + Parameters + ---------- + pixel : float[:, :] + The pre-computed fixed values for a given pixel + fixed : float[:, :] + The pre-computed fixed values for a given read_pattern + t_bar : float[:, :] + The average time for each resultant + slope : float + The computed slope for the ramp + ramp : RampIndex + Struct for start and end of ramp to fit + + Returns + ------- + argmax(all_stats), max(all_stats) + """ + # Cast the enum values into integers for indexing (otherwise compiler complains) + # These will be optimized out + cdef int single_local_slope = PixelOffsets.single_local_slope + cdef int double_local_slope = PixelOffsets.double_local_slope + cdef int single_var_read_noise = PixelOffsets.single_var_read_noise + cdef int double_var_read_noise = PixelOffsets.double_var_read_noise + + cdef int single_t_bar_diff_sqr = FixedOffsets.single_t_bar_diff_sqr + cdef int double_t_bar_diff_sqr = FixedOffsets.double_t_bar_diff_sqr + cdef int single_var_slope_val = FixedOffsets.single_var_slope_val + cdef int double_var_slope_val = FixedOffsets.double_var_slope_val + + # Note that a ramp consisting of a single point is degenerate and has no + # fit statistic so we bail out here + if ramp.start == ramp.end: + return 0, NAN + + # Start computing fit statistics + cdef float correct = _correction(t_bar, ramp, slope) + + # We are computing single and double differences of using the ramp's resultants. + # Each of these computations requires two points meaning that there are + # start - end - 1 possible differences. However, we cannot compute a double + # difference for the last point as there is no point after it. Therefore, + # We use this point's single difference as our initial guess for the fit + # statistic. Note that the fit statistic can technically be negative so + # this makes it much easier to compute a "lazy" max. + cdef int index = ramp.end - 1 + cdef int argmax = ramp.end - ramp.start - 1 + cdef float max_stat = _statstic(pixel[single_local_slope, index], + pixel[single_var_read_noise, index], + fixed[single_t_bar_diff_sqr, index], + fixed[single_var_slope_val, index], + slope, + correct) + + # Compute the rest of the fit statistics + cdef float stat + cdef int stat_index + for stat_index, index in enumerate(range(ramp.start, ramp.end - 1)): + # Compute max of single and double difference statistics + stat = fmaxf(_statstic(pixel[single_local_slope, index], + pixel[single_var_read_noise, index], + fixed[single_t_bar_diff_sqr, index], + fixed[single_var_slope_val, index], + slope, + correct), + _statstic(pixel[double_local_slope, index], + pixel[double_var_read_noise, index], + fixed[double_t_bar_diff_sqr, index], + fixed[double_var_slope_val, index], + slope, + correct)) + + # If this is larger than the current max, update the max + if stat > max_stat: + max_stat = stat + argmax = stat_index + + return argmax, max_stat + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cdef inline JumpFits fit_jumps(float[:] resultants, + int[:] dq, + float read_noise, + float[:] t_bar, + float[:] tau, + int[:] n_reads, + int n_resultants, + float[:, :] fixed, + float[:, :] pixel, + Thresh thresh, + bool use_jump, + bool include_diagnostic): + """ + Compute all the ramps for a single pixel using the Casertano+22 algorithm + with jump detection. + + Parameters + ---------- + resultants : float[:] + The resultants for the pixel + dq : int[:] + The dq flags for the pixel. This is modified in place, so the external + dq flag array will be modified as a side-effect. + read_noise : float + The read noise for the pixel. + ramps : RampQueue + RampQueue for initial ramps to fit for the pixel + multiple ramps are possible due to dq flags + t_bar : float[:] + The average time for each resultant + tau : float[:] + The time variance for each resultant + n_reads : int[:] + The number of reads for each resultant + n_resultants : int + The number of resultants for the pixel + fixed : float[:, :] + The jump detection pre-computed values for a given read_pattern + pixel : float[:, :] + A pre-allocated array for the jump detection fixed values for the + given pixel. This will be modified in place, it is passed in to avoid + re-allocating it for each pixel. + thresh : Thresh + The threshold parameter struct for jump detection + use_jump : bool + Turn on or off jump detection. + include_diagnostic : bool + Turn on or off recording all the diaganostic information on the fit + + Returns + ------- + RampFits struct of all the fits for a single pixel + """ + # Find initial set of ramps + cdef RampQueue ramps = init_ramps(dq, n_resultants) + + # Initialize algorithm + cdef JumpFits 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 int argmax, jump0, jump1 + cdef float max_stat + cdef float weight, total_weight = 0 + + # Fill in the jump detection pre-compute values for a single pixel + if use_jump: + pixel = _fill_pixel_values(pixel, resultants, fixed, read_noise, n_resultants) + + # Run while the Queue is non-empty + while not ramps.empty(): + # Remove top ramp of the stack to use + ramp = ramps.back() + ramps.pop_back() + + # Compute fit using the Casertano+22 algorithm + ramp_fit = fit_ramp(resultants, + t_bar, + tau, + n_reads, + read_noise, + ramp) + + # Run jump detection if enabled + if use_jump: + argmax, max_stat = _fit_statistic(pixel, + fixed, + t_bar, + ramp_fit.slope, + ramp) + + # Note that when a "ramp" is a single point, _fit_statistic returns + # a NaN for max_stat. Note that NaN > anything is always false so the + # result drops through as desired. + if max_stat > _threshold(thresh, ramp_fit.slope): + # Compute jump point to create two new ramps + # This jump point corresponds to the index of the largest + # statistic: + # argmax = 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. + # + jump0 = argmax + ramp.start + + # 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. + jump1 = jump0 + 1 + + # Update the dq flags + dq[jump0] = JUMP_DET + dq[jump1] = JUMP_DET + + # Record jump diagnotics + if include_diagnostic: + 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_back(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_back(RampIndex(jump1 + 1, ramp.end)) + + # Skip recording the ramp as it has a detected jump + continue + + # Start recording the the fit (no jump detected) + + # Record the diagnositcs + if include_diagnostic: + ramp_fits.fits.push_back(ramp_fit) + ramp_fits.index.push_back(ramp) + + # Start computing the averages using a lazy process + # Note we do not do anything in the NaN case for degenerate ramps + if not 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 + + # Finish computing averages using the lazy proces + 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 \ No newline at end of file diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd b/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd deleted file mode 100644 index 1539724b..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_pixel.pxd +++ /dev/null @@ -1,24 +0,0 @@ -from libcpp cimport bool -from libcpp.vector cimport vector - -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, vector[RampIndex] ramps, bool include_diagnostic) - - -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 deleted file mode 100644 index ddf04971..00000000 --- a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx +++ /dev/null @@ -1,559 +0,0 @@ -""" -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 libcpp cimport bool -from libc.math cimport sqrt, fabs -from libcpp.vector cimport vector - -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, vector[RampIndex] ramps, bool include_diagnostic): - """ - Compute all the ramps for a single pixel using the Casertano+22 algorithm - with jump detection. - - Parameters - ---------- - ramps : vector[RampIndex] - Vector 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.back() - ramps.pop_back() - - # 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 - if include_diagnostic: - 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_back(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_back(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 - if include_diagnostic: - 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 - if include_diagnostic: - 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 ** 2) * np.array(fixed.read_recip_coeffs) - - return pixel diff --git a/src/stcal/ramp_fitting/ols_cas22/_ramp.pxd b/src/stcal/ramp_fitting/ols_cas22/_ramp.pxd new file mode 100644 index 00000000..de31cd6c --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_ramp.pxd @@ -0,0 +1,40 @@ +# cython: language_level=3str + +from libcpp.vector cimport vector + + +cdef struct RampIndex: + int start + int end + + +cdef struct RampFit: + float slope + float read_var + float poisson_var + + +ctypedef vector[RampIndex] RampQueue + + +cdef class ReadPattern: + cdef float[::1] t_bar + cdef float[::1] tau + cdef int[::1] n_reads + + +cpdef RampQueue init_ramps(int[:] dq, + int n_resultants) + + +cpdef ReadPattern from_read_pattern(list[list[int]] read_pattern, + float read_time, + int n_resultants) + + +cdef RampFit fit_ramp(float[:] resultants_, + float[:] t_bar_, + float[:] tau_, + int[:] n_reads, + float read_noise, + RampIndex ramp) diff --git a/src/stcal/ramp_fitting/ols_cas22/_ramp.pyx b/src/stcal/ramp_fitting/ols_cas22/_ramp.pyx new file mode 100644 index 00000000..bd60e0fc --- /dev/null +++ b/src/stcal/ramp_fitting/ols_cas22/_ramp.pyx @@ -0,0 +1,358 @@ +# cython: language_level=3str + +""" +This module contains all the functions needed to execute the Casertano+22 ramp + fitting algorithm on its own without jump detection. + + The _jump module contains a driver function which calls the `fit_ramp` function + from this module iteratively. This evvetively handles dq flags and detected + jumps simultaneously. + +Structs +------- +RampIndex : struct + - start : int + Index of the first resultant in the ramp + - end : int + Index of the last resultant in the ramp (so indexing of ramp requires end + 1) + +RampFit : struct + - slope : float + The slope fit to the ramp + - read_var : float + The read noise variance for the fit + - poisson_var : float + The poisson variance for the fit + +RampQueue : vector[RampIndex] + Vector of RampIndex objects (convienience typedef) + +Classes +------- +ReadPattern : + Container class for all the metadata derived from the read pattern, this + is just a temporary object to allow us to return multiple memory views from + a single function. + +(Public) Functions +------------------ +init_ramps : function + Create the initial ramp "queue" for each pixel in order to handle any initial + "dq" flags passed in from outside. + +from_read_pattern : function + Derive the input data from the the read pattern + This is faster than using __init__ or __cinit__ to construct the object with + these calls. + +fit_ramps : function + Implementation of running the Casertano+22 algorithm on a (sub)set of resultants + listed for a single pixel +""" +import numpy as np +cimport numpy as cnp + +from cython cimport boundscheck, wraparound, cdivision, cpow + +from libc.math cimport sqrt, fabs, INFINITY, NAN, fmaxf +from libcpp.vector cimport vector + +from stcal.ramp_fitting.ols_cas22._ramp cimport RampIndex, RampQueue, RampFit, ReadPattern + + +# Initialize numpy for cython use in this module +cnp.import_array() + + +cdef class ReadPattern: + """ + Class to contain the read pattern derived metadata + This exists only to allow us to output multiple memory views at the same time + from the same cython function. This is needed because neither structs nor unions + can contain memory views. + + In the case of this code memory views are the fastest "safe" array data structure. + This class will immediately be unpacked into raw memory views, so that we avoid + any further overhead of swithcing between python and cython. + + Attributes: + ---------- + t_bar : np.ndarray[float_t, ndim=1] + The mean time of each resultant + tau : np.ndarray[float_t, ndim=1] + The variance in time of each resultant + n_reads : np.ndarray[cnp.int32_t, ndim=1] + The number of reads in each resultant + """ + + def _to_dict(ReadPattern self): + """ + This is a private method to convert the ReadPattern 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. + """ + return dict(t_bar=np.array(self.t_bar, dtype=np.float32), + tau=np.array(self.tau, dtype=np.float32), + n_reads=np.array(self.n_reads, dtype=np.int32)) + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cpdef ReadPattern from_read_pattern(list[list[int]] read_pattern, float read_time, int n_resultants): + """ + Derive the input data from the the read pattern + This is faster than using __init__ or __cinit__ to construct the object with + these calls. + + Parameters + ---------- + read pattern: list[list[int]] + read pattern for the image + read_time : float + Time to perform a readout. + n_resultants : int + Number of resultants in the image + + Returns + ------- + ReadPattern + Contains: + - t_bar + - tau + - n_reads + """ + + cdef ReadPattern data = ReadPattern() + data.t_bar = np.empty(n_resultants, dtype=np.float32) + data.tau = np.empty(n_resultants, dtype=np.float32) + data.n_reads = np.empty(n_resultants, dtype=np.int32) + + 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 + + +@boundscheck(False) +@wraparound(False) +cpdef inline RampQueue init_ramps(int[:] dq, int n_resultants): + """ + Create the initial ramp "queue" 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] + DQ array + n_resultants : int + Number of resultants + + Returns + ------- + RampQueue + vector of RampIndex objects + - vector with entry for each ramp found (last entry is last ramp found) + - RampIndex with start and end indices of the ramp in the resultants + """ + cdef RampQueue ramps = RampQueue() + + # 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 + cdef RampIndex ramp = RampIndex(-1, -1) + cdef int index_resultant + for index_resultant in range(n_resultants): + if ramp.start == -1: + # Looking for the start of a ramp + if dq[index_resultant] == 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] == 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_back(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_back(ramp) + + return ramps + +# Keeps the static type checker/highligher happy this has no actual effect +ctypedef float[6] _row + +# Casertano+2022, Table 2 +cdef _row[2] _PTABLE = [[-INFINITY, 5, 10, 20, 50, 100], + [ 0, 0.4, 1, 3, 6, 10 ]] + + +@boundscheck(False) +@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] + + +@boundscheck(False) +@wraparound(False) +@cdivision(True) +cdef inline RampFit fit_ramp(float[:] resultants_, + float[:] t_bar_, + float[:] tau_, + int[:] n_reads_, + float read_noise, + RampIndex ramp): + """ + Fit a single ramp using Casertano+22 algorithm. + + Parameters + ---------- + resultants_ : float[:] + All of the resultants for the pixel + t_bar_ : float[:] + All the t_bar values + tau_ : float[:] + All the tau values + n_reads_ : int[:] + All the n_reads values + read_noise : float + The read noise for the pixel + ramp : RampIndex + Struct for start and end of ramp to fit + + Returns + ------- + RampFit + struct containing + - 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(NAN, NAN, NAN) + + # Compute the fit + cdef int i = 0, j = 0 + + # Setup data for fitting (work over subset of data) to make things cleaner + # 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 = 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] + + # Compute mid point time + cdef int end = n_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 power = fmaxf(resultants[end] - resultants[0], 0) + power = power / sqrt(read_noise**2 + power) + power = _get_power(power) + + # 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 + # it is faster to generate a c++ vector than a numpy array + 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 + cdef float coeff + + # Issue when tbar[] == tbarmid causes exception otherwise + with 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 + coeff = (f0 * t_bar[i] - f1) * weights[i] / det + coeffs[i] = coeff + + # Casertano+22 Eq. 38 + ramp_fit.slope += coeff * resultants[i] + + # Casertano+22 Eq. 39 + ramp_fit.read_var += (coeff ** 2 * read_noise ** 2 / n_reads[i]) + + # Casertano+22 Eq 40 + # Note that this is an inversion of the indexing from the equation; + # however, commutivity of addition results in the same answer. This + # makes it so that we don't have to loop over all the resultants twice. + ramp_fit.poisson_var += coeff ** 2 * tau[i] + for j in range(i): + ramp_fit.poisson_var += (2 * coeff * coeffs[j] * t_bar[j]) + + return ramp_fit diff --git a/tests/test_jump_cas22.py b/tests/test_jump_cas22.py index 0bf1b6ed..a6966a02 100644 --- a/tests/test_jump_cas22.py +++ b/tests/test_jump_cas22.py @@ -2,11 +2,13 @@ import pytest from numpy.testing import assert_allclose -from stcal.ramp_fitting.ols_cas22._core import init_ramps, 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._jump import (fill_fixed_values, + _fill_pixel_values, + FixedOffsets, + PixelOffsets) +from stcal.ramp_fitting.ols_cas22._ramp import from_read_pattern, init_ramps -from stcal.ramp_fitting.ols_cas22 import fit_ramps, Parameter, Variance, Diff, RampJumpDQ +from stcal.ramp_fitting.ols_cas22 import fit_ramps, Parameter, Variance, JUMP_DET # Purposefully set a fixed seed so that the tests in this module are deterministic @@ -34,47 +36,6 @@ 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 @@ -90,7 +51,7 @@ def test_init_ramps(): [0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1]], dtype=np.int32) n_resultants, n_pixels = dq.shape - ramps = [init_ramps(dq, n_resultants, index_pixel) for index_pixel in range(n_pixels)] + ramps = [init_ramps(dq[:, index], n_resultants) for index in range(n_pixels)] assert len(ramps) == dq.shape[1] == 16 @@ -123,129 +84,106 @@ def test_init_ramps(): assert ramps[15] == [] -def test_threshold(): +@pytest.fixture(scope="module") +def read_pattern(): """ - Test the threshold object/fucnction - intercept - constant * log10(slope) = threshold + 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 """ + yield [ + [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] + ] + - # 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) - } +def test_from_read_pattern(read_pattern): + """Test turning read_pattern into the time data""" + metadata = from_read_pattern(read_pattern, READ_TIME, len(read_pattern))._to_dict() - # 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) + t_bar = metadata['t_bar'] + tau = metadata['tau'] + n_reads = metadata['n_reads'] + + # Check that the data is correct + assert_allclose(t_bar, [7.6, 15.2, 21.279999, 41.040001, 60.799999, 88.159996]) + assert_allclose(tau, [5.7, 15.2, 19.928888, 36.023998, 59.448887, 80.593781]) + assert np.all(n_reads == [4, 1, 3, 10, 3, 15]) - # 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) + # Check datatypes + assert t_bar.dtype == np.float32 + assert tau.dtype == np.float32 + assert n_reads.dtype == np.int32 @pytest.fixture(scope="module") -def ramp_data(base_ramp_data): +def ramp_data(read_pattern): """ - Unpacked metadata for simulating ramps for testing + Basic data for simulating ramps for testing (not unpacked) 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 : list[list[int]] + The example read pattern + metadata : dict + The metadata computed from 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) + data = from_read_pattern(read_pattern, READ_TIME, len(read_pattern))._to_dict() - yield read_pattern, t_bar, tau, n_reads + yield data['t_bar'], data['tau'], data['n_reads'], read_pattern -@pytest.mark.parametrize("use_jump", [True, False]) -def test_fixed_values_from_metadata(ramp_data, use_jump): +def test_fill_fixed_values(ramp_data): """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'] + t_bar, tau, n_reads, _ = ramp_data - # 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] - ) + n_resultants = len(t_bar) + fixed = np.empty((FixedOffsets.n_fixed_offsets, n_resultants - 1), dtype=np.float32) + fixed = fill_fixed_values(fixed, t_bar, tau, n_reads, n_resultants) + + # Sanity check that the shape of fixed is correct + assert fixed.shape == (2 * 4, n_resultants - 1) - 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] - 2 * 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] - 2 * 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() + # Split into the different types of data + t_bar_diffs = fixed[FixedOffsets.single_t_bar_diff:FixedOffsets.double_t_bar_diff + 1, :] + t_bar_diff_sqrs = fixed[FixedOffsets.single_t_bar_diff_sqr:FixedOffsets.double_t_bar_diff_sqr + 1, :] + read_recip = fixed[FixedOffsets.single_read_recip:FixedOffsets.double_read_recip + 1, :] + var_slope_vals = fixed[FixedOffsets.single_var_slope_val:FixedOffsets.double_var_slope_val + 1, :] + + # Sanity check that these are all the right shape + assert t_bar_diffs.shape == (2, n_resultants - 1) + assert t_bar_diff_sqrs.shape == (2, n_resultants - 1) + assert read_recip.shape == (2, n_resultants - 1) + assert var_slope_vals.shape == (2, n_resultants - 1) + + # Check the computed data + # These are computed using loop in cython, here we check against numpy + # Single diffs + assert np.all(t_bar_diffs[0] == t_bar[1:] - t_bar[:-1]) + assert np.all(t_bar_diff_sqrs[0] == (t_bar[1:] - t_bar[:-1])**2) + assert np.all(read_recip[0] == np.float32(1 / n_reads[1:]) + np.float32(1 / n_reads[:-1])) + assert np.all(var_slope_vals[0] == (tau[1:] + tau[:-1] - 2 * np.minimum(t_bar[1:], t_bar[:-1]))) + + # Double diffs + assert np.all(t_bar_diffs[1, :-1] == t_bar[2:] - t_bar[:-2]) + assert np.all(t_bar_diff_sqrs[1, :-1] == (t_bar[2:] - t_bar[:-2])**2) + assert np.all(read_recip[1, :-1] == np.float32(1 / n_reads[2:]) + np.float32(1 / n_reads[:-2])) + assert np.all(var_slope_vals[1, :-1] == (tau[2:] + tau[:-2] - 2 * np.minimum(t_bar[2:], t_bar[:-2]))) + + # Last double diff should be NaN + assert np.isnan(t_bar_diffs[1, -1]) + assert np.isnan(t_bar_diff_sqrs[1, -1]) + assert np.isnan(read_recip[1, -1]) + assert np.isnan(var_slope_vals[1, -1]) def _generate_resultants(read_pattern, n_pixels=1): @@ -305,76 +243,57 @@ def pixel_data(ramp_data): n_reads: The number of reads for the read pattern used for the resultants """ + t_bar, tau, n_reads, read_pattern = ramp_data + + n_resultants = len(t_bar) + fixed = np.empty((FixedOffsets.n_fixed_offsets, n_resultants - 1), dtype=np.float32) + fixed = fill_fixed_values(fixed, t_bar, tau, n_reads, n_resultants) - read_pattern, t_bar, tau, n_reads = ramp_data resultants = _generate_resultants(read_pattern) - yield resultants, t_bar, tau, n_reads + yield resultants, t_bar, tau, n_reads, fixed -@pytest.mark.parametrize("use_jump", [True, False]) -def test_make_pixel(pixel_data, use_jump): +def test__fill_pixel_values(pixel_data): """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 + resultants, t_bar, tau, n_reads, fixed = pixel_data + + n_resultants = len(t_bar) + pixel = np.empty((PixelOffsets.n_pixel_offsets, n_resultants - 1), dtype=np.float32) + pixel = _fill_pixel_values(pixel, resultants, fixed, READ_NOISE, n_resultants) + + # Sanity check that the shape of pixel is correct + assert pixel.shape == (2 * 2, n_resultants - 1) + + # Split into the different types of data + local_slopes = pixel[PixelOffsets.single_local_slope:PixelOffsets.double_local_slope + 1, :] + var_read_noise = pixel[PixelOffsets.single_var_read_noise:PixelOffsets.double_var_read_noise + 1, :] + + # Sanity check that these are all the right shape + assert local_slopes.shape == (2, n_resultants - 1) + assert var_read_noise.shape == (2, n_resultants - 1) # 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 == np.float32(READ_NOISE ** 2)* ( - 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 == np.float32(READ_NOISE ** 2) * ( - 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() + # These are computed using loop in cython, here we check against numpy + # Single diffs + assert np.all(local_slopes[0] == (resultants[1:] - resultants[:-1]) / (t_bar[1:] - t_bar[:-1])) + assert np.all(var_read_noise[0] == np.float32(READ_NOISE ** 2) * ( + np.float32(1 / n_reads[1:]) + np.float32(1 / n_reads[:-1])) + ) + + # Double diffs + assert np.all(local_slopes[1, :-1] == (resultants[2:] - resultants[:-2]) / (t_bar[2:] - t_bar[:-2])) + assert np.all(var_read_noise[1, :-1] == np.float32(READ_NOISE ** 2) * ( + np.float32(1 / n_reads[2:]) + np.float32(1 / n_reads[:-2])) + ) + + # Last double diff should be NaN + assert np.isnan(local_slopes[1, -1]) + assert np.isnan(var_read_noise[1, -1]) @pytest.fixture(scope="module") -def detector_data(ramp_data): +def detector_data(read_pattern): """ Generate a set of with no jumps data as if for a single detector as it would be passed in by the supporting code. @@ -387,7 +306,6 @@ def detector_data(ramp_data): 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) @@ -623,7 +541,7 @@ def test_override_default_threshold(jump_data): def test_jump_dq_set(jump_data): # Check the DQ flag value to start - assert RampJumpDQ.JUMP_DET == 2**2 + assert JUMP_DET == 2**2 resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data dq = np.zeros(resultants.shape, dtype=np.int32) @@ -633,7 +551,7 @@ def test_jump_dq_set(jump_data): 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() + assert (pixel_dq[fit['jumps']] == 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']) + assert set(np.where(pixel_dq == JUMP_DET)[0]) == set(fit['jumps'])