py-SP(k) Model

[jemme07]

Overview

This PR implements the SP(k) baryon-suppression model (Salcido et al. 2023) in CAMB.

Model references:

• pyspk reference implementation: https://github.com/jemme07/pyspk
• paper: https://ui.adsabs.harvard.edu/abs/2023MNRAS.523.2247S/abstract

What is implemented

1. SP(k) non-linear wrapper

• fortran/SPkNonLinear.f90 wraps a base non-linear model (e.g. HMCode/halofit) and applies SP(k) suppression multiplicatively.
• Supports SP(k) relation modes exposed in the interface:
  • power-law
  • cosmology power-law
  • double power-law
• Calibrated-range handling consistent with model limits:
  • outside calibrated z range, SP(k) correction is skipped
  • for k above calibrated limit, suppression evaluation is clamped at calibrated k_max
• One-time runtime warnings (for FeedbackLevel > 0) are emitted for both out-of-range behaviors.

2. Python API path and integration coverage

• SP(k) usage is validated through Python-side CAMB APIs, including class-selection.
• camb/tests/camb_test.py includes coverage that set_classes(non_linear_model="SPkNonLinear") selects and runs correctly.

3. Dedicated SP(k) validation tests

• camb/tests/spk_test.py adds dedicated SP(k) tests.
• Tests compare CAMB SP(k) outputs directly against pyspk across supported relation modes.
• Includes boundary-focused regression checks for high-k behavior and warning behavior.
• Uses stable numerical tolerances and sampling choices for reproducible comparisons.

4. SP(k) demo notebook

• Adds `docs/SPk_demo.ipynb
• Notebook demonstrates:
  • baseline non-linear vs non-linear+SP(k)
  • CAMB-vs-pyspk suppression overlays
  • relative error diagnostics as a function of k and z

Supporting updates

• pyproject.toml: adds pyspk>=2.0.0 to dev extras for validation workflows.

Validation run

• Focused SP(k) suite:
  • python -m unittest camb.tests.spk_test
  • passing

Reviewer notes

• The SP(k) formulas/behavior are anchored to the published model and cross-checked against the reference pyspk implementation.

[jemme07]

Hi @cmbant ,

This PR is ready for review. Could you please assign some reviewers?

Thanks!

[jemme07]

Any chance this is reviewed @cmbant?

[cmbant]

Odd, thought I'd add least added a comment.
As I recall, that some odd incidental changes in .vscode and camb_test, and asking whether this nested structure is compatible with Cobaya?

[jemme07]

As I recall, that some odd incidental changes in .vscode and camb_test, and asking whether this nested structure is compatible with Cobaya?

Thanks @cmbant !. The changes in .vscode and camb_test came from development/testing while I was iterating on the implementation, but I’ve cleaned those up in the latest revision.

And yes, the SPkNonLinear setup works correctly with Cobaya. Cobaya passes the relevant extra_args and sampled parameters through to the CAMB non-linear model configuration. I also added a Cobaya YAML example to the SPkNonLinear.set_params() docstring showing how to configure and sample the SP(k) parameters.

[cmbant]

AI review: The SP(k) integration builds cleanly, and the in-range power-law suppression matches pyspk to roundoff in a spot check. I do think there are two issues to address before merge:

1. fortran/SPkNonLinear.f90 / fortran/spk_model.f90: SP(k) does not enforce the pyspk baryon-fraction fitting limits.

   The reference implementation checks the computed baryon fraction against the calibrated fitting limits and returns NaN with warnings when values are outside the valid range. This implementation always applies a finite correction. For example, with SPk_fb_a=100, pyspk returns all NaN, while CAMB produces finite ratios up to about 1.33.

   This could silently turn invalid sampled parameters into apparently valid scientific output. The Fortran wrapper should either reproduce the fitting-limit behavior or explicitly reject/warn and skip/mark invalid cases.

2. camb/_ini.py: write_ini does not support SPkNonLinear.

   The PR adds new .ini parameters and inifiles/params_spk.ini, but Python-side write_ini() still only serializes plain Halofit. With SPkNonLinear and nonlinear enabled, it raises:

   CAMBValueError: write_ini does not support non-linear model class SPkNonLinear
   

[jemme07]

Thanks @cmbant, these comments should now be addressed.

The SP(k) implementation now reproduces the reference pyspk fitting-limit behaviour. In particular, the baryon-fraction limits are computed from the calibrated tables and any out-of-range points are set to NaN, with a one-time warning emitted when encountered. This avoids invalid parameter combinations.

I also added full write_ini() support for SPkNonLinear, together with regression coverage in the Python test suite.

I also updated the documentation to clarify the calibrated validity domain, boundary behaviour, NaN propagation, and MCMC/Cobaya implications, including guidance on using the interpolator interface for reliable NaN detection.

[cmbant]

AI followup: P1 fortran/spk_model.f90 (line 94): SPk_GetFbLimits uses linear interpolation for the fitting-limit coefficients, but the reference pyspk implementation uses Akima1DInterpolator for each coefficient before computing min_fb/max_fb (pyspk source). Since those limits decide whether fortran/SPkNonLinear.f90 (line 182) returns finite suppression or NaN, CAMB can accept/reject the wrong parameter points at off-node redshifts. The current NaN test uses z=0.5, which is a table node, so it would not catch this.

[jemme07]

I replaced the linear interpolation with Akima for the limit coefficients. I checked this, and the discrepancy was ≤1% for z < 1 (absorbed by the 20% safety margin) but grew to ~20% at z = 2.5. Now matches pyspk exactly at all redshifts. I also added an off-node-z test to prevent regression.

I updated the documentation to clarify interpolation guidance and NaN preservation.

[cmbant]

AI review (may be intentional, but ini inconsistency ma be real?): Medium — SPkNonLinear.set_params() cannot configure the wrapped Halofit model except for halofit_version, so common CAMB/HMCode options become unusable through the standard API.
The new wrapper exposes only halofit_version for the base Halofit and then calls self.BaseModel.set_params(halofit_version=halofit_version), which resets the base model’s HMCode_A_baryon, HMCode_eta_baryon, and HMCode_logT_AGN to Halofit defaults and provides no direct way for camb.set_params()/Cobaya-style flat arguments to set them. In current CAMB, camb.set_params() only forwards keyword arguments that appear in the active NonLinearModel.set_params signature, then treats remaining flat keys as attributes of CAMBparams, so HMCode_logT_AGN/HMCode_A_baryon/HMCode_eta_baryon are no longer accepted as flat non-linear-model arguments when non_linear_model="SPkNonLinear". This is particularly surprising because _ini.py does serialize those base-model values for SPkNonLinear, so the ini path acknowledges they are part of the wrapped model state, but the Python setter path cannot set them directly. PR diff lines show those base values being written for SPk ini output. ()
Suggested fix: Add optional HMCode_A_baryon=3.13, HMCode_eta_baryon=0.603, and HMCode_logT_AGN=7.8 parameters to SPkNonLinear.set_params() and pass all four Halofit arguments through to self.BaseModel.set_params(...). Keep the existing validation that rejects double-counting baryonic corrections when SPk_feedback=True.

[jemme07]

AI review (may be intentional, but ini inconsistency ma be real?): Medium — SPkNonLinear.set_params() cannot configure the wrapped Halofit model except for halofit_version, so common CAMB/HMCode options become unusable through the standard API.
The new wrapper exposes only halofit_version for the base Halofit and then calls self.BaseModel.set_params(halofit_version=halofit_version), which resets the base model’s HMCode_A_baryon, HMCode_eta_baryon, and HMCode_logT_AGN to Halofit defaults and provides no direct way for camb.set_params()/Cobaya-style flat arguments to set them. In current CAMB, camb.set_params() only forwards keyword arguments that appear in the active NonLinearModel.set_params signature, then treats remaining flat keys as attributes of CAMBparams, so HMCode_logT_AGN/HMCode_A_baryon/HMCode_eta_baryon are no longer accepted as flat non-linear-model arguments when non_linear_model="SPkNonLinear". This is particularly surprising because _ini.py does serialize those base-model values for SPkNonLinear, so the ini path acknowledges they are part of the wrapped model state, but the Python setter path cannot set them directly. PR diff lines show those base values being written for SPk ini output. ()
Suggested fix: Add optional HMCode_A_baryon=3.13, HMCode_eta_baryon=0.603, and HMCode_logT_AGN=7.8 parameters to SPkNonLinear.set_params() and pass all four Halofit arguments through to self.BaseModel.set_params(...). Keep the existing validation that rejects double-counting baryonic corrections when SPk_feedback=True.

Good point. Originally I intended to keep this independent, but yes, the three HMCode parameters were being silently reset to defaults on every set_params() call and were completely unreachable via the Python/Cobaya API, even though _ini.py already serialised them correctly for SPk.

I added all three to the signature (matching Halofit.set_params() defaults) and forwarded them through to the base model. This also makes it straightforward to compare HMCode's baryon correction (mead2020_feedback, SPk_feedback=False) directly against an SP(k) run.

I updated the docstring and spk.rst to make this clear, and tightened the signature test to cover the new params.

[cmbant]

For the next CAMB version AI tools made some optimizations/restructuring of your changes let me know if you think any issues:
https://github.com/cmbant/CAMB/blob/a5a7e64e336f72a487e1dcb1644641fd034cba25/fortran/SPkNonLinear.f90
(mainly preventing duplicate recalculations).

[jemme07]

For the next CAMB version AI tools made some optimizations/restructuring of your changes let me know if you think any issues: https://github.com/cmbant/CAMB/blob/a5a7e64e336f72a487e1dcb1644641fd034cba25/fortran/SPkNonLinear.f90 (mainly preventing duplicate recalculations).

Hi @cmbant,

Looks good to me.

I've verified the changes against master by building both commits and comparing the SPK suppression ratio on the raw internal k-grid (using get_linear_matter_power_spectrum(nonlinear=True) to avoid interpolation artifacts). The results are identical to machine precision across all redshifts z=0–3.

The changes make sense. Happy with this.