Skip to content

Commit fc00636

Browse files
committed
fix: coerce fnnls_cholesky inputs to numpy — JAX ZTx broke the numba buffer kernels
The sparse-operator inversion path hands fnnls_cholesky JAX arrays even when the fit itself runs the numba CPU path. Indexing a JAX array yields another JAX array, which numba maps to a *readonly* buffer — and the in-place buffer kernels introduced by #453 (_cho_solve_buffer via _solve_upper_transposed_buffer) overwrite their vector argument, so kernel compilation fails with "Cannot modify readonly array" the first time the solver runs on such input. The scipy solvers #453 replaced tolerated JAX input by copying internally, which is why this never bit before. Visible failure: HowToLens smoke on main, red since 2026-08-20 22:30 UTC — tutorial_8_adaptive_pixelization and tutorial_11_brightness_adaption (the two tutorials whose fits route through the sparse operator) died with "During: Pass nopython_type_inference" on both CI Pythons, 21 minutes after the #453 merge. This RED is what blocked the 2026-08-21 nightly release. Coerce ZTZ / ZTx / P_initial once at the function boundary: fancy indexing a numpy parent hands every downstream kernel a fresh writeable array. The new regression test calls fnnls_cholesky with jnp arrays (importorskip'd) and is verified to fail without the coercion. test_autoarray: 1063 passed; the 3 test_transformer pynufft failures reproduce identically on unmodified main in this environment (missing [optional] extras) and are unrelated. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_015GsUfbCwPd4XC8kpsiUJp7
1 parent d29ad97 commit fc00636

2 files changed

Lines changed: 40 additions & 0 deletions

File tree

autoarray/util/fnnls.py

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -32,6 +32,20 @@ def fnnls_cholesky(
3232
"""
3333
from scipy import linalg as slg
3434

35+
# The buffer kernels below (_cho_solve_buffer / cholinsertlast_inplace)
36+
# overwrite their vector argument in place, so every slice handed to them
37+
# must be a writeable numpy array. A JAX ZTZ / ZTx — the sparse-operator
38+
# inversion path hands one over even when the fit itself runs the numba
39+
# CPU path — breaks that contract: indexing a JAX array yields another
40+
# JAX array, which numba maps to a *readonly* buffer and rejects at
41+
# compile time ("Cannot modify readonly array"). Coerce once at the
42+
# boundary — fancy indexing a numpy parent then hands the kernels fresh
43+
# writeable arrays, exactly as the scipy solvers this replaced tolerated
44+
# by copying internally.
45+
ZTZ = np.asarray(ZTZ)
46+
ZTx = np.asarray(ZTx)
47+
P_initial = np.asarray(P_initial)
48+
3549
lstsq = lambda A, x: slg.solve(
3650
A,
3751
x,

test_autoarray/util/test_cholesky_inplace.py

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -167,3 +167,29 @@ def test__fnnls_cholesky__warm_start_matches_cold_start(seed):
167167
d_warm = fnnls_cholesky(ZTZ, ZTx, P_initial=P_initial)
168168

169169
assert d_warm == pytest.approx(d_cold, rel=1e-8, abs=1e-10)
170+
171+
172+
@pytest.mark.parametrize("seed", [0, 1])
173+
def test__fnnls_cholesky__accepts_jax_arrays(seed):
174+
"""
175+
The sparse-operator inversion path hands fnnls_cholesky JAX arrays even
176+
when the fit runs the numba CPU path. Indexing a JAX array yields another
177+
JAX array, which numba maps to a readonly buffer — before the boundary
178+
coercion in fnnls_cholesky this failed kernel compilation with
179+
"Cannot modify readonly array" (HowToLens smoke, 2026-08-20).
180+
"""
181+
jnp = pytest.importorskip("jax.numpy")
182+
183+
rng = np.random.default_rng(seed)
184+
n = 30
185+
Z = rng.normal(size=(50, n))
186+
x = Z @ rng.normal(size=n) + rng.normal(size=50)
187+
188+
ZTZ = Z.T @ Z
189+
ZTx = Z.T @ x
190+
191+
d_np = fnnls_cholesky(ZTZ, ZTx)
192+
d_jax = fnnls_cholesky(jnp.asarray(ZTZ), jnp.asarray(ZTx))
193+
194+
assert np.all(np.asarray(d_jax) >= 0.0)
195+
assert np.asarray(d_jax) == pytest.approx(d_np, rel=1e-6, abs=1e-8)

0 commit comments

Comments
 (0)