Skip to content

Commit

Permalink
Add further documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
WilliamJamieson committed Oct 6, 2023
1 parent 9d1e54a commit 0266cb1
Show file tree
Hide file tree
Showing 3 changed files with 49 additions and 24 deletions.
18 changes: 10 additions & 8 deletions src/stcal/ramp_fitting/ols_cas22/_fixed.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
10 changes: 6 additions & 4 deletions src/stcal/ramp_fitting/ols_cas22/_pixel.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -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
-----
Expand Down
45 changes: 33 additions & 12 deletions tests/test_jump_cas22.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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):
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 0266cb1

Please sign in to comment.