Skip to content

Commit

Permalink
Fix all negative indexing.
Browse files Browse the repository at this point in the history
This enables removing bounds/wrap checks which are slow
  • Loading branch information
WilliamJamieson committed Oct 4, 2023
1 parent a7dc16b commit 1f72ae2
Show file tree
Hide file tree
Showing 3 changed files with 32 additions and 11 deletions.
8 changes: 8 additions & 0 deletions src/stcal/ramp_fitting/ols_cas22/_core.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,10 @@ Functions
from libcpp.stack cimport stack
from libcpp.deque cimport deque
from libc.math cimport log10

import numpy as np
cimport numpy as np
cimport cython

from stcal.ramp_fitting.ols_cas22._core cimport Thresh, DerivedData

Expand All @@ -57,6 +59,8 @@ cdef float[2][6] PTABLE = [
[0, 0.4, 1, 3, 6, 10]]


@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline float get_power(float s):
"""
Return the power from Casertano+22, Table 2
Expand Down Expand Up @@ -96,6 +100,8 @@ cdef inline float threshold(Thresh thresh, float slope):
return thresh.intercept - thresh.constant * log10(slope)


@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline deque[stack[RampIndex]] init_ramps(int[:, :] dq):
"""
Create the initial ramp stack for each pixel
Expand Down Expand Up @@ -166,6 +172,8 @@ cdef inline deque[stack[RampIndex]] init_ramps(int[:, :] dq):
return pixel_ramps


@cython.boundscheck(False)
@cython.wraparound(False)
cdef DerivedData read_data(list[list[int]] read_pattern, float read_time):
"""
Derive the input data from the the read pattern
Expand Down
30 changes: 20 additions & 10 deletions src/stcal/ramp_fitting/ols_cas22/_fixed.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@ make_fixed : function
"""
import numpy as np
cimport numpy as np
cimport cython

from stcal.ramp_fitting.ols_cas22._core cimport Thresh, DerivedData, Diff
from stcal.ramp_fitting.ols_cas22._fixed cimport Fixed
Expand Down Expand Up @@ -63,6 +64,8 @@ cdef class Fixed:
from pre-computing the values and reusing them.
"""

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline float[:, :] t_bar_diff_val(Fixed self):
"""
Compute the difference offset of t_bar
Expand All @@ -81,15 +84,18 @@ cdef class Fixed:
# 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 = <float [:self.data.t_bar.size()]> self.data.t_bar.data()
cdef int end = len(t_bar)

cdef np.ndarray[float, ndim=2] t_bar_diff = np.zeros((2, self.data.t_bar.size() - 1), dtype=np.float32)

t_bar_diff[Diff.single, :] = np.subtract(t_bar[1:], t_bar[:-1])
t_bar_diff[Diff.double, :-1] = np.subtract(t_bar[2:], t_bar[:-2])
t_bar_diff[Diff.double, -1] = np.nan # last double difference is undefined
t_bar_diff[Diff.single, :] = np.subtract(t_bar[1:], t_bar[:end - 1])
t_bar_diff[Diff.double, :end - 2] = np.subtract(t_bar[2:], t_bar[:end - 2])
t_bar_diff[Diff.double, end - 2] = np.nan # last double difference is undefined

return t_bar_diff

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline float[:, :] recip_val(Fixed self):
"""
Compute the reciprical sum values
Expand All @@ -109,18 +115,21 @@ cdef class Fixed:
# 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 = <int [:self.data.n_reads.size()]> self.data.n_reads.data()
cdef int end = len(n_reads)

cdef np.ndarray[float, ndim=2] recip = np.zeros((2, self.data.n_reads.size() - 1), dtype=np.float32)

recip[Diff.single, :] = (np.divide(1.0, n_reads[1:], dtype=np.float32) +
np.divide(1.0, n_reads[:-1], dtype=np.float32))
recip[Diff.double, :-1] = (np.divide(1.0, n_reads[2:], dtype=np.float32) +
np.divide(1.0, n_reads[:-2], dtype=np.float32))
recip[Diff.double, -1] = np.nan # last double difference is undefined
np.divide(1.0, n_reads[:end - 1], dtype=np.float32))
recip[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))
recip[Diff.double, end - 2] = np.nan # last double difference is undefined

return recip


@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline float[:, :] slope_var_val(Fixed self):
"""
Compute slope part of the variance
Expand All @@ -140,12 +149,13 @@ cdef class Fixed:
# stays local to the function (numpy operations create brand new objects)
cdef float[:] t_bar = <float [:self.data.t_bar.size()]> self.data.t_bar.data()
cdef float[:] tau = <float [:self.data.tau.size()]> self.data.tau.data()
cdef int end = len(t_bar)

cdef np.ndarray[float, ndim=2] slope_var = np.zeros((2, self.data.t_bar.size() - 1), dtype=np.float32)

slope_var[Diff.single, :] = (np.add(tau[1:], tau[:-1]) - np.minimum(t_bar[1:], t_bar[:-1]))
slope_var[Diff.double, :-1] = (np.add(tau[2:], tau[:-2]) - np.minimum(t_bar[2:], t_bar[:-2]))
slope_var[Diff.double, -1] = np.nan # last double difference is undefined
slope_var[Diff.single, :] = (np.add(tau[1:], tau[:end - 1]) - np.minimum(t_bar[1:], t_bar[:end - 1]))
slope_var[Diff.double, :end - 2] = (np.add(tau[2:], tau[:end - 2]) - np.minimum(t_bar[2:], t_bar[:end - 2]))
slope_var[Diff.double, end - 2] = np.nan # last double difference is undefined

return slope_var

Expand Down
5 changes: 4 additions & 1 deletion src/stcal/ramp_fitting/ols_cas22/_pixel.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -209,6 +209,8 @@ cdef class Pixel:

return ramp_fit

@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline float correction(Pixel self, RampIndex ramp, int index, int diff):
"""
Compute the correction factor for the variance used by a statistic
Expand Down Expand Up @@ -274,7 +276,6 @@ cdef class Pixel:

@cython.boundscheck(False)
@cython.wraparound(False)
@cython.cdivision(True)
cdef inline float[:] stats(Pixel self, float slope, RampIndex ramp):
"""
Compute fit statistics for jump detection on a single ramp
Expand Down Expand Up @@ -469,6 +470,8 @@ cdef class Pixel:
return ramp_fits


@cython.boundscheck(False)
@cython.wraparound(False)
cdef inline Pixel make_pixel(Fixed fixed, float read_noise, float [:] resultants):
"""
Fast constructor for the Pixel C class.
Expand Down

0 comments on commit 1f72ae2

Please sign in to comment.