Skip to content

fix: JAX-trace CSE mass-profile decomposition (Sersic/NFW deflections) #499

Description

@Jammy2211

Overview

group/features/advanced/mass_stellar_dark/chaining.py fails fast (~26 s) on the
release stack with a TracerArrayConversionError, not a timeout. Diagnosed as a
real library bug, not a jax-0.10.2 regression: the CSE decomposition path used
by Sérsic stellar-mass and NFW dark-mass deflections was only partially ported
to JAX (commit 35e6812), so it cannot be jax.jit-traced when the mass params are
free. It surfaced now because this is the first release script to jit a CSE mass
profile with free params (isothermal deflects analytically and never hits CSE).

Plan

  • Thread xp into the deflection entry so the immediate crash clears.
  • Make cse_settings_from traceable: xp ops + jnp.where/jnp.minimum instead
    of np.log10/np.min, and freeze total_cses/sample_points to a static
    conservative max (~50/80)
    so decomposition array shapes are constant under jit
    (no Python branch on the sersic_index tracer, no value-dependent shapes).
  • Port _decompose_convergence_via_cse_from to xp.logspace, a vectorized
    coefficient-matrix build, and jnp.linalg.lstsq (replacing scipy.linalg.lstsq).
  • Port both decompose_convergence_via_cse impls (Sérsic + NFW) to xp.
  • Keep the NumPy path (xp=np) behaviourally identical; validate the JAX path via
    the autogalaxy_workspace_test parity scripts (jit round-trip + fitness._vmap),
    never test_autogalaxy/; confirm mass_stellar_dark/chaining.py runs green.

Phased: P1 — threading + traceable cse_settings_from/decompose to unblock the
release tail. P2 — parity hardening across the Sérsic-index range and NFW.

Detailed implementation plan

Affected Repositories

  • PyAutoGalaxy (primary)
  • PyAutoArray, PyAutoLens (reproduction only — no changes expected)

Branch Survey

Repository Current Branch Dirty?
./PyAutoGalaxy main clean
./PyAutoArray main clean
./PyAutoLens main clean

Suggested branch: feature/mass-cse-jax-decompose

Root cause (reproduced on jax 0.10.2, JAX_TRACEBACK_FILTERING=off)

search[1] (lens light) cached; search[2] (lens mass = Sérsic + NFW) crashes on
its first jitted likelihood eval. The decorator AbstractMaker._xp
(PyAutoArray/.../decorators/abstract.py) resolves the backend from the xp passed
into the decorated call, so a dropped xp silently becomes numpy on a tracer.

  1. cse.py:197 radial_grid_from(grid=grid, **kwargs) drops xpnp.sqrt on a
    float64[225] tracer.
  2. After threading that, sersic.py:58 cse_settings_fromnp.log10((23.0/ sersic_constant)**sersic_index) on a scalar tracer; whole decompose uses
    np.logspace on tracer bounds + scipy.linalg.lstsq (untraceable). NFW decompose
    (np.max(grid_radii), complex coord_func_g) has the same shape.

Implementation Steps (P1)

  1. autogalaxy/profiles/mass/abstract/cse.py:197 — pass xp=xp to radial_grid_from.
  2. autogalaxy/profiles/mass/stellar/sersic.py cse_settings_from (lines ~1–82) —
    remove Python branching on sersic_index; compute upper_dex/lower_dex with
    xp/jnp.where/jnp.minimum; return constant total_cses/sample_points
    (conservative max, ~50/80) so shapes are static.
  3. autogalaxy/profiles/mass/abstract/cse.py _decompose_convergence_via_cse_from
    (lines ~85–147) — xp.logspace; build coefficient_matrix vectorized (no
    np.zeros+loop); jnp.linalg.lstsq when xp is jnp, else keep scipy for numpy
    parity/speed.
  4. sersic.py:245 + dark/nfw.py:177 decompose_convergence_via_csexp the
    np.sqrt/np.max/np.real ops and radii bounds.

Validation

  • autogalaxy_workspace_test JAX parity (jit round-trip + fitness._vmap) on a
    Sérsic mass + NFW model; NumPy path unchanged (test_autogalaxy/ stays green,
    NumPy-only).
  • End-to-end: mass_stellar_dark/chaining.py under the release profile
    (PYAUTO_TEST_MODE=1) completes without a traceback.

Key Files

  • autogalaxy/profiles/mass/abstract/cse.py — decomposition + deflection driver.
  • autogalaxy/profiles/mass/stellar/sersic.pycse_settings_from, Sérsic decompose.
  • autogalaxy/profiles/mass/dark/nfw.py — NFW decompose.
  • autogalaxy/profiles/geometry_profiles.py:113radial_grid_from (crash site).

Classification

Real library bug (incomplete JAX port), not a jax-0.10.2 regression. Surfaced
on release-validation run 29279095224 / TestPyPI 2026.7.13.1.dev65601
(PyAutoHeart#72).

Original Prompt

Click to expand starting prompt

group/mass_stellar_dark chaining fails with a JAX exception on release wheels

Type: bug
Target: autolens_workspace
Repos:

  • autolens_workspace
  • PyAutoGalaxy
  • PyAutoArray
    Difficulty: hard
    Autonomy: supervised
    Priority: normal
    Status: diagnosed
    Fix-locus: PyAutoGalaxy (autogalaxy/profiles/mass — CSE decomposition path)
    Workflow: library
    Strategy: full jnp trace of the CSE decomposition; freeze total_cses/sample_points
    to a conservative static max (~50/80) on both paths (static shapes under jit);
    phased — P1 unblock chaining, P2 Sersic-range + NFW parity harden.

Surfaced burning down the 2026-07-13 release-validation tail (PyAutoHeart#72). On the
release wheels under the release profile:

scripts/group/features/advanced/mass_stellar_dark/chaining.py ... FAIL (26.6s)
  For simplicity, JAX has removed its internal frames from the traceback of the following
  exception. Set JAX_TRACEBACK_FILTERING=off to include these.

The run-runner captured only the filtered banner, not the actual exception — it is a real
JAX error
during a group-scale mass_stellar_dark chaining pipeline (real sampler,
PYAUTO_TEST_MODE=1 for user workspaces), and fails fast (26.6 s), so it is NOT a timeout
(distinct from the other group//chaining scripts that TIMEOUT and are handled by the
release timeout bump PyAutoHeart#c19c948).

Reproduce from the autolens_workspace checkout on current main under the release profile,
with JAX_TRACEBACK_FILTERING=off to get the real traceback:

cd autolens_workspace
JAX_TRACEBACK_FILTERING=off PYAUTO_TEST_MODE=1 PYAUTO_SKIP_WORKSPACE_VERSION_CHECK=1 \
  python scripts/group/features/advanced/mass_stellar_dark/chaining.py

The dev venv is now jax 0.10.2 (matches the release stack — see
[[reference_laptop_gpu_jax_setup]]), so this should reproduce locally. Candidate causes given
the release-validation context: a jax-0.10.2 numerical/behaviour change in the chained
group-scale fit (cf. the jax-drift and LinAlg-not-PD siblings), a NaN/inf in a chained-prior
pass, or a real group/mass_stellar_dark modelling issue. Classify (jax regression vs real
bug) and fix at the right layer. Release run: PyAutoHeart workspace-validation 29279095224,
TestPyPI 2026.7.13.1.dev65601. See [[project_release_2026_07_13_blocked_3bugs]] and
PyAutoHeart#72.

Diagnosis (2026-07-14, reproduced locally on jax 0.10.2)

Classification: real library bug (incomplete JAX port of the CSE mass-profile
decomposition), NOT a jax-0.10.2 regression.
jax 0.10.2 is incidental — the
release stack simply jits the likelihood; the same code fails on any jax version.

Reproduced on autolens_workspace main with JAX_TRACEBACK_FILTERING=off.
search[1] (lens light) is cached; search[2]__mass_stellar_dark crashes on the
first jitted likelihood eval. The lens mass is a Sersic (stellar) + NFW (dark)
model; the Sersic mass profile deflections go through the CSE decomposition
path, which is only partially JAX-threaded:

  1. PyAutoArray .../decorators/abstract.py (AbstractMaker._xp) resolves the
    backend from the xp passed into the decorated call. So a dropped xp at a
    call site silently becomes xp=np.
  2. First crash — PyAutoGalaxy/autogalaxy/profiles/mass/abstract/cse.py:197:
    self.radial_grid_from(grid=grid, **kwargs) drops xp (it is a named param
    of _deflections_2d_via_cse_from, not in **kwargs), so radial_grid_from
    (geometry_profiles.py:113 xp.sqrt(...)) runs with xp=np on a JAX tracer →
    TracerArrayConversionError (float64[225]). Threading xp=xp here clears it.
  3. Second crash (after fixing Coordinates #2) — the decomposition itself is not traceable.
    With sersic_index free in search[2], decompose_convergence_via_cse runs on
    tracer params via numpy/scipy:
    • sersic.py:58 cse_settings_fromnp.log10((23.0/sersic_constant)**sersic_index)
      on a scalar tracer → TracerArrayConversionError (float64[]).
    • sersic.py:280 np.sqrt(self.axis_ratio()); param-dependent radii_min/max.
    • cse.py:122/126 np.logspace(np.log10(radii_min), ...) on tracer bounds.
    • cse.py:143 scipy.linalg.lstsq — cannot trace at all.
      The dark NFW decompose (dark/nfw.py) has the same shape (np.max(grid_radii),
      scipy.lstsq) and would fail identically once reached.

Isothermal traces fine because it deflects analytically (xp.arctan/arctanh),
never touching CSE — so the CSE JAX port (commit 35e6812a, 2026-05-26 "port CSE
module to support JAX via xp parameter") was never actually exercised under jit
until this chaining script. This is the first release script to jit a CSE mass
profile with free params.

Fix locus: PyAutoGalaxy (profiles/mass/abstract/cse.py, stellar/sersic.py
incl. cse_settings_from, dark/nfw.py). Real fix ≈ thread xp through
cse_settings_from + both decompose_convergence_via_cse + _decompose_convergence_via_cse_from,
and replace scipy.linalg.lstsq/np.logspace with jnp.linalg.lstsq/xp.logspace
so the decomposition is jit-traceable. Validate via the JAX parity scripts in
autogalaxy_workspace_test (jit round-trip + fitness._vmap), never test_autogalaxy/.
Difficulty raised medium → hard (numerical path + lstsq-in-trace, affects all
CSE mass profiles under JAX).

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions