fix: derive numba PSF kernel y/x shifts from the correct kernel axes - #487
Merged
Conversation
Both numba PSF gathers in inversion_imaging_numba_util.py computed their
kernel half-widths from the transposed kernel axes -- the y shift from
kernel_native.shape[1] (the x axis) and the x shift from shape[0] (the y
axis):
kernel_shift_y = -(kernel_native.shape[1] // 2)
kernel_shift_x = -(kernel_native.shape[0] // 2)
For a square kernel the two are equal and nothing is observably wrong, which
is why every existing test passed. For a non-square odd PSF the gather is
mis-centred along both axes. This is reachable: kernels are validated as odd
per axis, never as square (exc.KernelException in operators/convolver.py:268
and structures/grids/uniform_2d.py:1153 both check parity only). With the
bounds guard from #456 in place the mis-centred reads are clipped rather than
reading uninitialized memory, so the failure mode is a silent wrong answer
rather than a crash.
In psf_precision_value_from the same variables also drive the pair-overlap
early-exit, so the transposition additionally dropped genuinely overlapping
image-pixel pairs on one axis while admitting non-overlapping ones on the
other. The swap corrects that test too; no separate change is needed.
Both sites are fixed together deliberately: they must agree on kernel
orientation, or the psf_weighted_data and psf_precision_operator paths would
disagree with each other.
Reproduced on the pre-fix source with a 3x5 kernel on a fully-unmasked 5x5
image: max|numba - numpy| = 1420.0, with the numba result being the numpy
result shifted; the square-kernel control matched to 0.0. Post-fix both
match to 0.0.
Tests:
- Parametrised the numba-vs-numpy equivalence test and the precision-operator
edge test over 3x3, 3x5, 5x3 and 5x7 kernels, so both orientations of
non-square are covered rather than only the wide one.
- Corrected the _reference_value helper in the precision-operator edge test,
which replicated the very transposition it was meant to catch and so was
not the independent reference its comment claimed.
- Added two direct orientation probes, one per gather, that read back which
pixel was actually gathered (single-tap kernel over a coordinate-encoding
map) instead of re-deriving the kernel walk.
Verified as a detector: against the unpatched source these fail 8 tests --
every non-square parametrisation plus both orientation probes -- while the
square-kernel cases still pass, which is precisely why the defect survived.
A repo-wide scan for half-width derivations confirms these were the only two
transposed sites; convolve_with_kernel_native, the numpy twin, the JAX
precision operator and array_2d_util resize centring were already correct.
Full test_autoarray suite: 1154 passed, 55 skipped.
Co-Authored-By: Claude Opus 5 <[email protected]>
Claude-Session: https://claude.ai/code/session_01KkZgzrxmQ9YQpSRVatEydp
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes #486.
The defect
Both numba PSF gathers in
autoarray/inversion/inversion/imaging_numba/inversion_imaging_numba_util.pycomputed their kernel half-widths from the transposed kernel axes — the y shift fromkernel_native.shape[1](the x axis) and the x shift fromshape[0](the y axis):at
psf_weighted_data_fromandpsf_precision_value_from. For a square kernel the two are equal and nothing is observably wrong — which is why every existing test passed. For a non-square odd PSF the gather is mis-centred along both axes.This is reachable. Kernels are validated as odd per axis, never as square:
exc.KernelException("Convolver Convolver must be odd")inoperators/convolver.py:268and the same parity-only check instructures/grids/uniform_2d.py:1153. With the bounds guard from #456 in place the mis-centred reads are clipped rather than reading uninitialized memory, so the failure mode is a silent wrong answer, not a crash.One thing the issue didn't capture. In
psf_precision_value_fromthese same variables also drive the pair-overlap early-exit:So for a 3x5 kernel the transposition also admitted y offsets up to ±4 (should be ±2) while rejecting x offsets beyond ±2 (should be ±4) — genuinely overlapping image-pixel pairs were silently dropped on one axis and non-overlapping ones admitted on the other. The swap corrects that test too; no separate change was needed.
Both sites are fixed together deliberately: they must agree on kernel orientation, or the
psf_weighted_dataandpsf_precision_operatorpaths would disagree with each other.Reproduction
On the pre-fix source, a 3x5 kernel on a fully-unmasked 5x5 image:
140, 205, 180, 146, 104, 500, 660, 540352, 500, 660, 540, 412, 723, 980, 1240max|numba - numpy| = 1420.0, with the numba result being the numpy result shifted — exactly what a transposed half-width predicts. The square-kernel control matched to0.0. Post-fix both match to0.0.Changes
autoarray/inversion/inversion/imaging_numba/inversion_imaging_numba_util.py— the two-line swap at each of the two sites, plus a comment naming the[y, x]indexing convention so the next reader does not re-transpose it. No signature, call-site or behavioural change for square kernels.test_autoarray/inversion/inversion/imaging/test_inversion_imaging_util.py_reference_valuehelper in the precision-operator edge test. Its comment called it an "independent reference", but it replicated the very transposition it was meant to catch (kernel_shift_y = -(kw // 2)), agreeing with the buggy code only because its kernel was square. Left as it was, the new non-square case would have asserted the bug.10 * (y + 1) + (x + 1)makes the returned value name the pixel actually gathered — so the test reads back the orientation instead of re-deriving the kernel walk. A wide kernel and its tall transpose must return different, individually predictable values, which is exactly what a square kernel cannot show.Validation
The tests are verified as a detector, not just as passing tests. Against the unpatched source they fail 8 tests — every non-square parametrisation of both functions, plus both orientation probes — while the 6 square-kernel cases still pass, which is precisely why this defect survived:
With the fix applied: 1154 passed, 55 skipped across the full
test_autoarraysuite.blackclean.Orientation audit
A repo-wide scan for half-width derivations (
shape[0] // 2/shape[1] // 2) confirms these were the only transposed sites in the library. Everything else was already correct and is untouched:inversion_imaging_numba_util.pyconvolve_with_kernel_native—ky, kx = psf_kernel.shape; cy, cx = ky // 2, kx // 2.inversion_imaging_numba_util.pypsf_precision_operator_sparse_from—kernel_overlap_sizeis a symmetric product of both axes, orientation-agnostic.inversion_imaging_util.py:44-45— the zero-padded numpy twin, the reference implementation.inversion_imaging_util.py:409/479/747— the JAX precision operator.structures/arrays/array_2d_util.py:311-319— resize centring.Downstream impact
None. This is an internal numba util and behaviour is bit-identical for square kernels, which is every PSF in the workspaces today. No workspace changes required.
🤖 Generated with Claude Code
https://claude.ai/code/session_01KkZgzrxmQ9YQpSRVatEydp
Generated by Claude Code