From 0266cb165cc034cc8f9fa09197f5388e598ae7d8 Mon Sep 17 00:00:00 2001 From: William Jamieson Date: Fri, 6 Oct 2023 15:03:05 -0400 Subject: [PATCH] Add further documentation --- src/stcal/ramp_fitting/ols_cas22/_fixed.pyx | 18 +++++---- src/stcal/ramp_fitting/ols_cas22/_pixel.pyx | 10 +++-- tests/test_jump_cas22.py | 45 +++++++++++++++------ 3 files changed, 49 insertions(+), 24 deletions(-) diff --git a/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx b/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx index fe0dbdd58..6296bbaf6 100644 --- a/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx +++ b/src/stcal/ramp_fitting/ols_cas22/_fixed.pyx @@ -55,24 +55,26 @@ cdef class FixedValues: the jump detection statistics. These are formed from the reciprocal sum of the number of reads. single sum of reciprocal n_reads: - recip[Diff.single, :] = ((1/n_reads[i+1]) + (1/n_reads[i])) + read_recip_coeffs[Diff.single, :] = ((1/n_reads[i+1]) + (1/n_reads[i])) double sum of reciprocal n_reads: - recip[Diff.double, :] = ((1/n_reads[i+2]) + (1/n_reads[i])) + 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: - slope_var[Diff.single, :] = ([tau[i] + tau[i+1] - min(t_bar[i], t_bar[i+1])) + var_slope_coeffs[Diff.single, :] = (tau[i] + tau[i+1] + - min(t_bar[i], t_bar[i+1])) double of slope variance term: - slope_var[Diff.double, :] = ([tau[i] + tau[i+2] - min(t_bar[i], t_bar[i+2])) + var_slope_coeffs[Diff.double, :] = (tau[i] + tau[i+2] + - min(t_bar[i], t_bar[i+2])) Notes ----- - - t_bar_diffs, 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. + - 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. diff --git a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx index e376123e6..6991811b9 100644 --- a/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx +++ b/src/stcal/ramp_fitting/ols_cas22/_pixel.pyx @@ -46,15 +46,17 @@ cdef class Pixel: local_slopes : float [:, :] These are the local slopes between the resultants for the pixel. single difference local slope: - delta[Diff.single, :] = (resultants[i+1] - resultants[i]) / (t_bar[i+1] - t_bar[i]) + local_slopes[Diff.single, :] = (resultants[i+1] - resultants[i]) + / (t_bar[i+1] - t_bar[i]) double difference local slope: - delta[Diff.double, :] = (resultants[i+2] - resultants[i]) / (t_bar[i+2] - t_bar[i]) + 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: - sigma[Diff.single, :] = read_noise * ((1/n_reads[i+1]) + (1/n_reads[i])) + var_read_noise[Diff.single, :] = read_noise * ((1/n_reads[i+1]) + (1/n_reads[i])) double difference read_noise variance: - sigma[Diff.doule, :] = read_noise * ((1/n_reads[i+2]) + (1/n_reads[i])) + var_read_noise[Diff.doule, :] = read_noise * ((1/n_reads[i+2]) + (1/n_reads[i])) Notes ----- diff --git a/tests/test_jump_cas22.py b/tests/test_jump_cas22.py index e8c566020..347e253f2 100644 --- a/tests/test_jump_cas22.py +++ b/tests/test_jump_cas22.py @@ -9,12 +9,27 @@ from stcal.ramp_fitting.ols_cas22 import fit_ramps, Parameter, Variance, Diff, RampJumpDQ +# Purposefully set a fixed seed so that the tests in this module are deterministic RNG = np.random.default_rng(619) -ROMAN_READ_TIME = 3.04 -READ_NOISE = np.float32(5) -N_PIXELS = 100_000 + +# The read time is constant for the given telescope/instrument so we set it here +# to be the one for Roman as it is known to be a reasonable value +READ_TIME = 3.04 + +# Choose small read noise relative to the flux to make it extremely unlikely +# that the random process will "accidentally" generate a set of data, which +# can trigger jump detection. This makes it easier to cleanly test jump +# detection is doing what we expect. FLUX = 100 +READ_NOISE = np.float32(5) + +# Set a value for jumps which makes them obvious relative to the normal flux JUMP_VALUE = 10_000 + +# Choose reasonable values for arbitrary test parameters, these are kept the same +# across all tests to make it easier to isolate the effects of something using +# multiple tests. +N_PIXELS = 100_000 CHI2_TOL = 0.03 GOOD_PROB = 0.7 @@ -40,7 +55,7 @@ def base_ramp_data(): [22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36] ] - yield read_pattern, metadata_from_read_pattern(read_pattern, ROMAN_READ_TIME) + yield read_pattern, metadata_from_read_pattern(read_pattern, READ_TIME) def test_metadata_from_read_pattern(base_ramp_data): @@ -215,7 +230,7 @@ def _generate_resultants(read_pattern, n_pixels=1): for _ in reads: # Compute the next value of the ramp # Using a Poisson process for the flux - ramp_value += RNG.poisson(FLUX * ROMAN_READ_TIME, size=n_pixels).astype(np.float32) + ramp_value += RNG.poisson(FLUX * READ_TIME, size=n_pixels).astype(np.float32) # Add to running total for the resultant resultant_total += ramp_value @@ -350,13 +365,19 @@ def test_fit_ramps(detector_data, use_jump, use_dq): if not use_dq: assert okay.all() - output = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=use_jump) + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=use_jump) assert len(output.fits) == N_PIXELS # sanity check that a fit is output for each pixel chi2 = 0 for fit, use in zip(output.fits, okay): if not use_dq: - assert len(fit['fits']) == 1 # only one fit per pixel since no dq/jump in this case + # Check that the data generated does not generate any false positives + # for jumps as this data is reused for `test_find_jumps` below. + # This guarantees that all jumps detected in that test are the + # purposefully placed ones which we know about. So the `test_find_jumps` + # can focus on checking that the jumps found are the correct ones, + # and that all jumps introduced are detected properly. + assert len(fit['fits']) == 1 if use: # Add okay ramps to chi2 @@ -382,7 +403,7 @@ def test_fit_ramps_array_outputs(detector_data, use_jump): resultants, read_noise, read_pattern = detector_data dq = np.zeros(resultants.shape, dtype=np.int32) - output = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=use_jump) + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=use_jump) for fit, par, var in zip(output.fits, output.parameters, output.variances): assert par[Parameter.intercept] == 0 @@ -454,7 +475,7 @@ def test_find_jumps(jump_data): resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data dq = np.zeros(resultants.shape, dtype=np.int32) - output = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=True) + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) assert len(output.fits) == len(jump_reads) # sanity check that a fit/jump is set for every pixel chi2 = 0 @@ -513,8 +534,8 @@ def test_override_default_threshold(jump_data): resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data dq = np.zeros(resultants.shape, dtype=np.int32) - standard = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=True) - override = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=True, + standard = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) + override = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True, intercept=0, constant=0) # All this is intended to do is show that with all other things being equal passing non-default @@ -529,7 +550,7 @@ def test_jump_dq_set(jump_data): resultants, read_noise, read_pattern, jump_reads, jump_resultants = jump_data dq = np.zeros(resultants.shape, dtype=np.int32) - output = fit_ramps(resultants, dq, read_noise, ROMAN_READ_TIME, read_pattern, use_jump=True) + output = fit_ramps(resultants, dq, read_noise, READ_TIME, read_pattern, use_jump=True) for fit, pixel_dq in zip(output.fits, output.dq.transpose()): # Check that all jumps found get marked