File tree Expand file tree Collapse file tree
Expand file tree Collapse file tree Original file line number Diff line number Diff 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 ,
Original file line number Diff line number Diff 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 )
You can’t perform that action at this time.
0 commit comments