Skip to content

refactor: collapse 6 compute_cl_pp_* into single compute_cl_pp(... nonlinear=...) — redo of #11 with bug fixes - #14

Open
MinhMPA wants to merge 5 commits into
smsharma:mainfrom
MinhMPA:refactor/cl_pp_consolidation_v2
Open

refactor: collapse 6 compute_cl_pp_* into single compute_cl_pp(... nonlinear=...) — redo of #11 with bug fixes#14
MinhMPA wants to merge 5 commits into
smsharma:mainfrom
MinhMPA:refactor/cl_pp_consolidation_v2

Conversation

@MinhMPA

@MinhMPA MinhMPA commented May 3, 2026

Copy link
Copy Markdown
Collaborator

Summary

Re-attempts the C_ℓ^φφ API consolidation that was merged as PR #11 and reverted via PR #13, this time with the three test failures that surfaced post-merge folded in from the start. The fix list is the same content as the (now-obsoleted) draft hotfix PR #12, but combined into a single working PR so it lands cleanly without an interim regression.

Why a redo, not a re-merge of #11

PR #11 passed the contract tests pre-merge but failed three other tests that I hadn't run beforehand:

  1. _halofit_modulator was vmap-unsafe. It passed an extended k-grid to compute_linear_matter_pk_from_perturbations, which (a) rejects k_eval outside pt.k_grid via _validate_k_eval_support (ValueError: k_eval must lie within...), and (b) has a Python-level if z == 0.0 branch that fails under jax.vmap-over-z (TracerBoolConversionError).
  2. compute_pk_nonlinear's Python-level _sigma_convergence_check was bypassed under vmap (the try/except catches ConcretizationTypeError), so Halofit ran at high z where σ < 1 and produced spurious R values (~1.7 at z = 100), inflating the integrated NL boost ~25% above CLASS at high ℓ.
  3. tests/test_lensing.py::test_cl_pp_positive referenced a pipeline_fast_cl fixture that doesn't exist on upstream/main (it lives in PR feat: clax-pt — differentiable one-loop EFTofLSS galaxy power spectra in JAX #9's conftest).

This PR ships PR #11's API consolidation with all three issues fixed in a single commit.

Public API

clax.compute_cl_pp(pt, params, bg, th, l_max, *, nonlinear="none")
# nonlinear ∈ {"none", "halofit"}; anything else raises ValueError.
# nonlinear="ept" is reserved for a follow-up clax-pt PR (#9).

Re-exported at the package root (clax.compute_cl_pp) along with clax.lens_cls.

Changes

Removed (BREAKING — pre-1.0, hard break)

Function Why
compute_cl_pp (Siddharth's original) Three-bug Bessel integration; superseded by source-Limber
compute_cl_pp_fast Inaccurate at ℓ ≥ 300 (documented in its own docstring)
compute_cl_pp_vmap Hermite Bessel-table vmap; superseded
compute_cl_pp_limber Poisson-reconstruction Limber; ~20% overestimate at ℓ = 2500 vs CLASS

Renamed

From To Rationale
compute_cl_pp_source_limber compute_cl_pp Sole public entry; CLASS-accurate at all ℓ ≤ 2500
compute_cl_pp_transfer _compute_cl_pp_full_bessel Private oracle for cross-impl tests

Added — _halofit_modulator (CLASS-aligned)

  • n_z = 100 z-points (CLASS density), log-spaced in 1+z spanning [0, z_rec].
  • Inline sigma_R(R_min, ...) >= 1 check using jnp.where to fall back to pk_lin when σ < 1, replicating CLASS fourier.c:1706-1716 in a vmap-safe way.
  • k_max_extend = 0 default — uses pt.k_grid as-is, matching CLASS's no-extrapolation behavior. The inline σ check forces R = 1 when the k-grid is too narrow. Pass k_max_extend > 0 to enable log-log power-law extension as a numerical convenience for narrow grids.
  • 2D-interpolates R(k, z) onto every (pt.k_grid, pt.tau_grid) lattice point using bg.loga_of_tau for the τ → z mapping.
  • Multiplies sqrt(R) directly into S_transfer (CLASS source-multiplication recipe).

Tests

File Change
tests/test_cl_pp.py (renamed from test_cl_pp_source_limber.py) Contract + linear-accuracy + cross-impl-vs-_compute_cl_pp_full_bessel + Halofit smoke + JIT/AD compatibility
tests/test_clpp_halofit_ratio.py Rewritten for new API; <7% NL/lin ratio agreement at ℓ ≤ 500, <10% at ℓ ≥ 1000 vs CLASS Halofit reference
tests/test_lensing.py Callers updated to new signature; uses local pipeline fixture matching upstream-main convention
tests/test_cl_pp_implementations.py Deleted (targeted removed implementations)
tests/test_cl_pp_limber.py Deleted (targeted removed compute_cl_pp_limber)

Test plan

Verified locally on this branch:

  • pytest tests/test_cl_pp.py tests/test_lensing.py tests/test_clpp_halofit_ratio.py -v → 32 passed in 295s
  • pytest tests/test_ept_*.py -v → 36 passed (sanity check; this PR doesn't touch ept)
  • 104 other tests (background, thermo, nonlinear, rosenbrock, etc.) → all pass

Expected on this PR:

  • CI: pytest tests/ -v --fast -q passes
  • No remaining references to deleted functions: git grep -nE "compute_cl_pp_(fast|vmap|limber|source_limber|transfer)" clax/ tests/ should return only _compute_cl_pp_full_bessel

Migration

Hard break, no deprecation shims:

# Before (any of):
cl = compute_cl_pp_source_limber(pt, params, bg, th, l_max=2500)
cl = compute_cl_pp_vmap(pt, params, bg, l_max=2500)
cl = compute_cl_pp_fast(pt, params, bg, l_max=2500)
cl = compute_cl_pp_limber(pt, params, bg, th, l_max=2500, nonlinear=True)

# After:
cl = clax.compute_cl_pp(pt, params, bg, th, l_max=2500)
cl_nl = clax.compute_cl_pp(pt, params, bg, th, l_max=2500, nonlinear="halofit")

…nlinear=...) — redo of smsharma#11 with bug fixes

Supersedes PR smsharma#11 (which was reverted via PR smsharma#13 due to post-merge test
failures). This PR ships the same API consolidation with the three fixes
that smsharma#11 was missing folded in from the start.

Public API:

    compute_cl_pp(pt, params, bg, th, l_max, *, nonlinear="none")
    # nonlinear in {"none", "halofit"}; ValueError for anything else.
    # "ept" is reserved for a follow-up clax-pt PR.

Removed (BREAKING):
- compute_cl_pp           (Siddharth's original; superseded by source-Limber)
- compute_cl_pp_fast      (inaccurate at l >= 300 per its own docstring)
- compute_cl_pp_vmap      (Hermite Bessel-table vmap; superseded)
- compute_cl_pp_limber    (Poisson-reconstruction; ~20% over at l=2500)

Renamed:
- compute_cl_pp_source_limber -> compute_cl_pp (sole public entry)
- compute_cl_pp_transfer      -> _compute_cl_pp_full_bessel (private oracle)

Halofit injection (new private _halofit_modulator):
- vmap(compute_pk_nonlinear) over a 100-point z-grid (CLASS-aligned density).
- Inline `sigma_R(R_min, ...) >= 1` check + jnp.where to force R=1 where
  Halofit isn't applicable — replicates CLASS fourier.c:1706-1716 in a
  vmap-safe way (the Python-level check inside compute_pk_nonlinear gets
  bypassed under vmap-over-z, which is why smsharma#11 over-corrected ~25% at
  high l).
- k_max_extend = 0 default (no power-law extension; matches CLASS's
  no-extrapolation behavior). Pass a positive value to enable log-log
  extension for narrow k-grids.
- 2D-interpolates R(k, z) onto (pt.k_grid, pt.tau_grid) via bg.loga_of_tau.
- Multiplies sqrt(R) directly into S_transfer (CLASS source-multiplication).

Tests:
- tests/test_cl_pp.py (renamed from test_cl_pp_source_limber.py): contract
  + linear-accuracy + cross-impl-vs-_compute_cl_pp_full_bessel + Halofit
  smoke + JIT/AD.
- tests/test_clpp_halofit_ratio.py: rewritten for new API; <7% NL/lin
  ratio agreement at l<=500, <10% at l>=1000 vs CLASS Halofit reference.
- tests/test_lensing.py: callers updated; uses local `pipeline` fixture
  matching upstream-main convention.
- Deleted: tests/test_cl_pp_implementations.py, test_cl_pp_limber.py.

Verified locally: 32/32 lensing+halofit tests pass, 36/36 EPT tests pass,
104/104 remaining tests pass.
MinhMPA added 3 commits May 3, 2026 13:17
…ute_cl_pp

Documents both paths in the new unified API:
- nonlinear="none" (linear source-Limber): <1% vs CLASS at all l<=2500
- nonlinear="halofit" (z-aware Halofit injection via source multiplication):
  <7% at l=500 and <10% at l>=1000 vs CLASS Halofit reference

Existing tables (unlensed C_l, lensed C_l with CLASS-generated cl_pp,
matter P(k), pipeline, performance, multi-cosmology) are unaffected by
PR smsharma#14 — they don't depend on compute_cl_pp.

Note on lensed C_l accuracy table: it tests lens_cls in isolation by
feeding CLASS-generated C_l^pp; the reported <0.2% is about lens_cls's
correlation-function method, not about our compute_cl_pp.
…<=2500

The previous table understated linear accuracy. test_clpp_halofit_ratio.py
::TestClppLinear asserts <1% at l in {100, 500, 1000, 2000, 2500} with
pt_k_max_cl=5.0. Updated the linear column accordingly and added the
prerequisite (pt_k_max_cl >= 5) to the section preamble.
…ratio test thresholds

Sweep across modulator settings on default cosmology (l_max=2500):

| l    | n_z=100, no extension | n_z=100, extend to k=10 |
| 100  | 0.01%                 | 0.01%                   |
| 200  | 0.04%                 | 0.04%                   |
| 500  | 0.42%                 | 0.21%                   |
| 1000 | 1.47%                 | 0.63%                   |
| 1500 | 2.30%                 | 0.79%                   |
| 2000 | 3.74%                 | 1.40%                   |
| 2500 | 3.97%                 | 0.96%                   |

The k-extension restores accuracy at high l. CLASS itself uses a dedicated
nonlinear k-grid that extends past the perturbation k-range; our log-log
power-law extension provides equivalent coverage with σ(R) bisection
well-resolved.

Doubling n_z to 200 didn't help (slight degradation at l=1000, l=2500
from interpolation noise), so 100 is the right z-density.

Changes:
- _halofit_modulator default k_max_extend: 0 → 10 Mpc^-1
- test_ratio_at_low_l threshold: 7% → 1% (measured: 0.21% at l=500)
- test_ratio_at_high_l threshold: 10% → 2% (measured: 1.40% worst case)
- README accuracy table: replace bounds with measured residuals
- Docstrings updated with rationale and measured numbers

User asked 'who set the 10% threshold' — it was carried over from the
deleted Poisson-Limber test, never measured against the actual source-
multiplication recipe. Now it reflects what the code actually achieves.
MinhMPA added a commit to MinhMPA/clax-pt that referenced this pull request May 3, 2026
…ratio test thresholds

Same change as PR smsharma#14 commit 1e28e12 (cherry-picked manually due to
EPT-related conflicts in feat/clax-pt's lensing.py).

Sweep across modulator settings on default cosmology (l_max=2500):

| l    | k_max_extend=0 | k_max_extend=10 |
| 100  | 0.01%          | 0.01%           |
| 200  | 0.04%          | 0.04%           |
| 500  | 0.42%          | 0.21%           |
| 1000 | 1.47%          | 0.63%           |
| 1500 | 2.30%          | 0.79%           |
| 2000 | 3.74%          | 1.40%           |
| 2500 | 3.97%          | 0.96%           |

CLASS itself uses a dedicated nonlinear k-grid extending past the
perturbation k-range; our log-log power-law extension provides the
equivalent coverage with σ(R) bisection well-resolved.

Test thresholds:
- low_l (l<=500): 7% -> 1% (measured 0.21% worst)
- high_l (l>=1000): 10% -> 2% (measured 1.40% worst)
…=2500); fair absolute-accuracy table

k_max_extend sweep on default cosmology, n_z=100 fixed:

| l    | k=10  | k=20  | k=30  |
| 500  | 0.21% | 0.09% | 0.09% |
| 1000 | 0.63% | 0.04% | 0.17% |
| 1500 | 0.79% | 0.11% | 0.11% |
| 2000 | 1.40% | 0.05% | 0.28% |
| 2500 | 0.96% | 0.76% | 0.76% |

(metric: |R_clax / R_CLASS - 1|, R = NL/lin)

k=20 is the sweet spot — far enough for sigma(R) bisection to fully
converge, not so far that the local power-law slope drifts from the
true asymptote. Going to k=30+ degrades because power-law extrapolation
diverges from the true transfer-function shape.

Timing benchmark (5-run mean): k=10 110ms, k=20 110ms, k=30 111ms.
Accuracy boost is essentially free.

Test thresholds:
- low_l (l<=500): 1% -> 0.5%
- high_l (l>=1000): 2% -> 1%

README accuracy table replaced with apples-to-apples absolute residuals
for both 'none' and 'halofit' (was mixing absolute-vs-ratio metrics).
Both modes share the source-Limber kernel and track each other within
~0.1% on absolute scale; the ratio (third column) isolates the R(k,z)
quality from the common linear-kernel systematic.
MinhMPA added a commit to MinhMPA/clax-pt that referenced this pull request May 3, 2026
…=2500)

Mirrors PR smsharma#14 commit 17c7e5a. Sweep across modulator settings on default
cosmology, n_z=100 fixed:

| l    | k=10  | k=20  | k=30  |
| 500  | 0.21% | 0.09% | 0.09% |
| 1000 | 0.63% | 0.04% | 0.17% |
| 1500 | 0.79% | 0.11% | 0.11% |
| 2000 | 1.40% | 0.05% | 0.28% |
| 2500 | 0.96% | 0.76% | 0.76% |

(metric: |R_clax / R_CLASS - 1|, R = NL/lin)

Timing benchmark (5-run mean on the same fixture):
  k=10  110ms ± 3ms
  k=20  110ms ± 2ms
  k=30  111ms ± 2ms

Accuracy boost is essentially free.

Test thresholds tightened: low_l 1%->0.5%, high_l 2%->1%.
@MinhMPA
MinhMPA marked this pull request as ready for review May 3, 2026 05:05
@MinhMPA

MinhMPA commented May 3, 2026

Copy link
Copy Markdown
Collaborator Author

Self-flag for the reviewer (no action required in this PR):

clax/lensing.py::_halofit_modulator inlines the τ-axis interpolation of pt.delta_m and the P_lin formula (2π²/k³ × P_R × δ_m²) rather than reusing clax/transfer.py::compute_pk_from_perturbations + compute_linear_matter_pk_from_perturbations. That's deliberate code duplication — those public APIs have two vmap-unfriendly choices that the modulator can't tolerate when looped over a z-grid:

  1. _validate_k_eval_support rejects k_eval outside pt.k_grid, which the modulator's k-extension to k_max=20 Mpc⁻¹ deliberately violates.
  2. The Python-level if z == 0.0: fast-path in compute_pk_from_perturbations raises TracerBoolConversionError under jax.vmap over z.

Both have small fixes in transfer.py (a validate=False kwarg + jnp.where for the z=0 branch), after which _halofit_modulator could drop ~30 lines and call the shared utilities directly. Out of scope for this PR — flagging in case you'd prefer that cleanup folded in before merge or as a separate follow-up. The tests already cover the modulator's correctness end-to-end either way.

EDIT: This is now addressed in PR #15. Merge order: #14 first, then #15 and #16 (independent of each other, both touch different files).

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant