CosmoFit.cosmology

Calculators

The pieces every model shares: the background evolution, the distance integrals, recombination, the sound horizon and linear growth. A model supplies E(z); these supply everything derived from it.

BackgroundCalculator
DistanceCalculator	Cosmological distance calculator.
SoundHorizon	The sound horizon at the drag epoch for a given cosmology.
RecombinationCalculator	Recombination-epoch quantities for a given cosmology, following the CMB distance-prior conventions (see the module docstring).
GrowthCalculator	Fast evaluator for the linear growth factor D(z), growth rate f(z) = dlnD/dlna, and fsigma8(z), given a Cosmology.
DistanceIntegrator	Fast evaluator for the dimensionless comoving distance.

Custom models

Minimal-code custom cosmological models.

CosmoFit’s built-in models (LCDM, CPL, …) are hand-written Cosmology subclasses. Testing a genuinely new model – one not in the literature, invented to check an idea – against CosmoFit’s existing datasets/likelihoods/MCMC machinery shouldn’t require writing one. define_model() builds a Cosmology subclass from a single E(z) function (everything downstream – distances, every likelihood, every plot except deceleration()/w_of_z() – only ever needs E(z), see Cosmology.__init_subclass__).

Example

>>> from CosmoFit import define_model, Fitter
>>> import numpy as np
>>>
>>> MyModel = define_model(
...     "MyModel",
...     E=lambda p, z: np.sqrt(
...         p["Omega_m"] * (1 + z) ** 3
...         + (1 - p["Omega_m"]) * (1 + z) ** (3 * (1 + p["w0"])) * (1 + p["beta"] * z)
...     ),
...     extra_params={"beta": {"default": 0.0, "bounds": (-2.0, 2.0), "label": r"$\beta$"}},
... )
>>>
>>> fit = Fitter(
...     model=MyModel,
...     datasets=["cc", "desi"],
...     free_params=["H0", "Omega_m", "w0", "beta"],
...     initial={"H0": 67.4, "Omega_m": 0.315, "w0": -1.0, "beta": 0.0},
... )
>>> fit.run_mcmc(nwalkers=48, nsteps=3000, burnin=500)
>>> fit.best_fit()
>>> fit.plots.corner()

E (and, if given, w/dEdz/Omega_de) receives (p, z): p is a plain dict of every current parameter value (standard ones – H0, Omega_m, Omega_k, w0, wa, … – plus any extra_params), and z is already an ndarray. It must be written with numpy operations (it is called on a whole redshift grid at once, not per-point) and return an ndarray the same shape as z.

If dEdz isn’t supplied, a central-finite-difference fallback is installed so fitter.plots.deceleration() works without deriving one by hand – direct Cosmology subclasses don’t get this fallback (they keep the base class’s strict NotImplementedError, same as every built-in model).

CosmoFit.cosmology.custom.define_model(name, E, *, extra_params=None, label=None, w=None, dEdz=None, Omega_de=None, mu=None)

Build a new Cosmology subclass from a plain E(z) function.

Parameters:

• name (str) – Model name (model.MODEL_NAME and the generated class’s __name__).

• E (callable(params: dict, z: ndarray) -> ndarray) – Dimensionless Hubble parameter. Required – this alone is enough to fit the model against every CosmoFit dataset and produce every plot except w_of_z() / deceleration().

• extra_params (dict[str, dict], optional) – New parameters this model needs beyond the standard set (H0, Omega_m, Omega_k, w0, wa, rd, MB, Omega_b, A_s, alpha). Maps name -> spec dict with keys "default" (float, default 0.0), "bounds" ((lower, upper), needed only if the parameter will be fit rather than fixed, unless bounds are instead passed to Fitter(bounds=...)), and "label" (str, optional, LaTeX label for corner plots). Names must not collide with the standard set.

• label (str, optional) – How the model’s name should appear in figure legends and titles (model.MODEL_LABEL), when name is really a set of symbols spelled out as an identifier – e.g. name="MyQuintessence", label=r"$\phi$CDM". Matplotlib renders the $...$ spans with mathtext. Defaults to name used verbatim, which is right for an acronym.

• w (callable(params, z) -> ndarray, optional) – Dark-energy equation of state, for w_of_z(). Not required for fitting.

• dEdz (callable(params, z) -> ndarray, optional) – Derivative of E(z), for deceleration() / background.q(). Defaults to a numerical (central finite-difference) approximation of E if omitted.

• Omega_de (callable(params, z) -> ndarray, optional) – Dark-energy density parameter, for background.Omega_de().

• mu (callable(params: dict, a: ndarray, k: float or None) -> ndarray, optional) – Effective-to-Newtonian gravitational coupling G_eff(a,k)/G_N, for growth-of-structure predictions (background.{growth_rate,sigma8,fsigma8}, the "fsigma8"/"s8" datasets) – see mu(). a is the scale factor (an ndarray, evaluated on a whole grid at once, same convention as E); k is a single wavenumber [h/Mpc] (or None, for a scale-independent mu that never reads it). Defaults to 1 everywhere (standard GR growth) if omitted – correct for any model that reparametrizes dark energy without touching gravity itself; only give this for a genuinely modified-gravity model, exactly as FQExponential/FRTLinear/ FRHuSawicki do internally.

Returns:

type – A new Cosmology subclass, usable directly as Fitter(model=..., ...).

Return type:

type

CosmoFit.cosmology.custom.model_from_expression(name, E, *, extra_params=None, label=None, w=None, dEdz=None, Omega_de=None, mu=None)

Same as define_model(), but E/w/dEdz/ Omega_de/mu are given as Python expression strings (e.g. "sqrt(Omega_m*(1+z)**3 + (1-Omega_m)*(1+z)**(3*(1+w0)))") instead of callables – convenient for text-entry UIs (see the Streamlit app under app/), where asking for a Python function isn’t practical. E/w/dEdz/Omega_de expressions see z and every current parameter value (standard ones plus any extra_params) as plain names; mu instead sees a (scale factor) and k (wavenumber [h/Mpc], or None if not supplied by the caller) – see _compile_expression() for exactly what else is available.

Example

>>> MyModel = model_from_expression(
...     "MyModel",
...     E="sqrt(Omega_m*(1+z)**3 + (1-Omega_m)*(1+z)**(3*(1+w0))*(1+beta*z))",
...     extra_params={"beta": {"default": 0.0, "bounds": (-2.0, 2.0)}},
... )

A modified-gravity example (custom growth on top of an otherwise LCDM background, in the spirit of FRTLinear):

>>> MyMG = model_from_expression(
...     "MyMG",
...     E="sqrt(Omega_m*(1+z)**3 + (1-Omega_m))",
...     mu="1 + 3*beta",
...     extra_params={"beta": {"default": 0.0, "bounds": (-0.2, 0.2)}},
... )

Parameters:

• name (str)

• E (str)

• extra_params (dict | None)

• label (str | None)

• w (str | None)

• dEdz (str | None)

• Omega_de (str | None)

• mu (str | None)

Return type:

type

The Boltzmann backend

Boltzmann-code backend: CMB angular power spectra from scratch.

Everything else in CosmoFit is a background calculation – distances and expansion rates from E(z), plus a one-equation linear growth ODE. That is enough for CC, BAO, supernovae and the compressed CMB distance priors, and it is why the library runs an MCMC in minutes without a Fortran dependency.

It is not enough to predict C_l. The CMB anisotropy spectrum comes out of the coupled Boltzmann hierarchy for photons, neutrinos, baryons and cold dark matter, integrated through recombination with a full ionization history – thousands of coupled ODEs per wavenumber, over hundreds of wavenumbers. Reimplementing that here would be a Boltzmann code, not a feature of one, and a pure-Python one would be far too slow to put inside an MCMC.

So this module does the honest thing and calls one. CAMBBackend translates a Cosmology into CAMB’s parameter conventions, runs it, and hands back binned-ready C_l arrays. CAMB is an optional dependency (pip install "cosmofit[cmb]"); nothing else in the library imports this module, and a fit that does not include the "planck_lite" dataset never touches it.

What can and cannot be pushed through it

CAMB solves the perturbation equations for a specific set of physical models. A CosmoFit model that is only an E(z) cannot be handed to it – the same E(z) is consistent with many different perturbation histories, and picking one silently would be inventing physics the model never specified. Three cases:

• LCDM maps exactly onto CAMB’s default.

• Any model exposing a dark-energy equation of state ``w(z)`` (wCDM, CPL, JBP, BA, GCG) is passed through CAMB’s PPF dark-energy module as a tabulated w(a). This is exact at the background level and is the standard treatment of the perturbations for a smooth dark-energy fluid, including across w = -1, where a quintessence-fluid treatment breaks down.

• Modified-gravity models (f(Q), f(R,T), f(R)) are refused. Their whole content is that the field equations differ from GR, which is exactly what CAMB’s perturbation solver assumes. A model like FRHuSawicki would run – its background is LCDM’s by construction – and would return LCDM’s C_l while its f_R0 did nothing, which is worse than an error. So this raises instead.

Where a model is refused, the compressed distance priors ("planck") still work: they only need E(z).

CosmoFit.cosmology.boltzmann.NEUTRINO_MASS_DENOM = 93.14

Neutrino mass-to-density conversion, omega_nu h^2 = m_nu / NEUTRINO_MASS_DENOM [eV]. The standard value for the temperature and degeneracy CAMB assumes.

CosmoFit.cosmology.boltzmann.supports_cmb_spectra(model)

Whether CMB power spectra can be computed for a model, and why not when they cannot – without importing CAMB.

Answers the question CAMBBackend’s constructor answers by raising, but as a value, for callers that need to decide before offering the choice: a GUI greying out a dataset, a script picking between the compressed priors and the full spectra, a table of what each model supports.

Parameters:

model (type or Cosmology) – A model class or an instance of one.

Returns:

(bool, str) – (True, "") if supported, else (False, reason).

Return type:

tuple[bool, str]

exception CosmoFit.cosmology.boltzmann.BoltzmannError

Bases: RuntimeError

Raised when CMB spectra cannot be computed for this cosmology – because CAMB is not installed, because the model’s perturbations are outside what CAMB can represent, or because CAMB itself rejected the parameters.

class CosmoFit.cosmology.boltzmann.CAMBBackend(cosmology, lmax=2508, lens_potential_accuracy=1)

Bases: object

Computes CMB angular power spectra for a CosmoFit cosmology by calling CAMB.

Parameters:

• cosmology (Cosmology) – The model to compute spectra for. Read live on every call, so an in-place params.update(theta) during an MCMC step is picked up without rebuilding the backend.

• lmax (int, optional) – Maximum multipole to compute. CAMB is asked for somewhat more than this internally, since its high-l accuracy degrades near the requested limit.

• lens_potential_accuracy (int, optional) – CAMB’s lensing accuracy setting. 1 (the default) is what Planck analyses use for parameter estimation; 0 disables lensing entirely and is wrong at the 10-sigma level for Planck’s error bars, so it is not offered as a shortcut.

Notes

Parameter translation is where the errors hide, so it is spelled out:

• H0, Omega_b h^2, Omega_k, tau, n_s map across directly.

• A_s is exp(ln1e10As) * 1e-10. The library’s A_s field is the Chaplygin gas parameter and is not this; see CosmologyParameters.

• omch2 is Omega_m h^2 - Omega_b h^2 - Omega_nu h^2. CosmoFit counts massive neutrinos inside Omega_m (they are non-relativistic across the whole redshift range its background calculations cover), while CAMB counts them separately – so the neutrino density is subtracted here rather than added. Getting this backwards shifts Omega_c h^2 by ~0.0006, which is ~0.5 sigma of Planck’s constraint on it.

classmethod attached(cosmology)

The backend already attached to cosmology, or None.

Unlike shared() this never creates one. It exists for callers that want the Boltzmann code’s answer if it is already being computed, and must not trigger a CAMB run of their own if it is not – deriving sigma8 for the growth machinery, in particular, which would otherwise make a growth-only fit pay for a CMB calculation nothing asked for.

Return type:

CAMBBackend | None

classmethod shared(cosmology, lmax=2508, lens_potential_accuracy=1)

The backend attached to cosmology, creating it on first use and widening it if this caller needs more than the last one did.

Two likelihoods can want CMB spectra from the same cosmology – plik_lite and the lensing reconstruction do, in any fit that uses the full Planck data. Giving each its own backend runs CAMB twice per MCMC step for two views of one calculation, which doubles the cost of the single most expensive thing in the library. Measured on a 615-bandpower + lensing fit: 1.93 s per evaluation with separate backends, 1.36 s sharing one. Not the clean halving it looks like it should be, because sharing also raises the accuracy the bandpower half is computed at – see the widening rule below.

Widening rather than asserting equality, because the two callers legitimately differ: the bandpower likelihood needs lmax = 2508 at accuracy 1, the lensing one lmax = 2500 at accuracy 4. The union – the larger of each – is correct for both, since more multipoles and more accuracy can only help. Widening invalidates the cache, so a already-computed result is never served at the lower setting it was computed with.

Parameters:

• lmax (int)

• lens_potential_accuracy (int)

Return type:

CAMBBackend

property omega_nu_h2: float

Physical density of the massive neutrinos, from m_nu.

sigma8()

sigma_8 as the Boltzmann code derives it, from the primordial amplitude and the transfer function.

This is not cosmology.sigma8. That one is a free parameter the growth machinery (GrowthCalculator, and the "fsigma8"/"s8" likelihoods) uses to normalize its own scale-independent growth factor. This one is a derived quantity, fixed by ln1e10As, n_s, tau_reio and the densities.

A fit that varies the free sigma8 while also using a CAMB-based CMB likelihood is therefore carrying two different amplitudes that nothing forces to agree. Fitter warns about that combination; this method is how to check it.

Return type:

float

lensing_spectra(lmax)

The inputs Planck’s lensing likelihood is defined on, indexed from l = 0 so array index and multipole coincide.

Returns:

dict – {"TT", "EE", "TE"} as D_l = l(l+1) C_l / 2 pi in muK^2, and "PP" as [L(L+1)]^2 C_L^{phiphi} / 2 pi.

Parameters:

lmax (int)

Return type:

dict[str, ndarray]

Notes

These scalings are not stylistic – they are what the bundled window functions were built against. The PP one is the easy one to get wrong: Planck’s bandpowers are ~1.5e-7 at L ~ 30, where C_L^{phiphi} itself is ~1e-8 and L(L+1)C_L/2pi is ~1.3e-6, so those two would land orders of magnitude out and be caught immediately – but a stray 2 pi would not, which is why the convention is written down rather than left to be inferred.

cls(lmin=2)

Lensed CMB angular power spectra.

Parameters:

lmin (int, optional) – First multipole to return.

Returns:

dict – {"ell", "TT", "TE", "EE"}, with the spectra given as C_l in muK^2 (not D_l), which is what the Planck bandpower windows are defined on.

Return type:

dict[str, ndarray]