Metadata-Version: 2.4
Name: cosmofit
Version: 2.0.0
Summary: A modular Python library for cosmological parameter estimation.
Author: Salih Yeşil
License-Expression: MIT
Project-URL: Homepage, https://github.com/salihyesil59/CosmoFit
Project-URL: Repository, https://github.com/salihyesil59/CosmoFit
Project-URL: Issues, https://github.com/salihyesil59/CosmoFit/issues
Project-URL: Documentation, https://salihyesil59.github.io/CosmoFit/
Project-URL: Changelog, https://github.com/salihyesil59/CosmoFit/blob/main/CHANGELOG.md
Project-URL: Examples, https://github.com/salihyesil59/CosmoFit/tree/main/examples
Keywords: cosmology,astronomy,dark-energy,bao,supernova,cc,cmb,planck,camb,boltzmann,power-spectrum,cmb-lensing,reionization,act,wcdm,cpl,lcdm,jbp,chaplygin-gas,des-sn5yr,union3,sdss,eboss,boss,6dfgs,bbn,sh0es,pede,dgp,cardassian,running-vacuum,interacting-dark-energy,sound-horizon,drag-epoch,neutrinos,likelihood,mcmc,hdf5
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: Intended Audience :: Education
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Topic :: Scientific/Engineering :: Astronomy
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=2.0
Requires-Dist: scipy>=1.14
Requires-Dist: pandas>=2.2
Requires-Dist: matplotlib>=3.9
Requires-Dist: emcee>=3.1
Requires-Dist: h5py>=3.10
Requires-Dist: corner>=2.2
Requires-Dist: tqdm>=4.66
Requires-Dist: threadpoolctl>=3.0
Requires-Dist: pyyaml>=6.0
Provides-Extra: dev
Requires-Dist: pytest; extra == "dev"
Requires-Dist: ruff; extra == "dev"
Requires-Dist: mypy; extra == "dev"
Requires-Dist: coverage; extra == "dev"
Requires-Dist: ipykernel; extra == "dev"
Requires-Dist: jupyterlab; extra == "dev"
Provides-Extra: speed
Requires-Dist: numba>=0.60; extra == "speed"
Provides-Extra: evidence
Requires-Dist: dynesty>=2.1; extra == "evidence"
Provides-Extra: theory
Requires-Dist: sympy>=1.12; extra == "theory"
Provides-Extra: cmb
Requires-Dist: camb>=1.5; extra == "cmb"
Provides-Extra: gui
Requires-Dist: streamlit>=1.38; extra == "gui"
Provides-Extra: docs
Requires-Dist: sphinx>=8; extra == "docs"
Requires-Dist: furo; extra == "docs"
Requires-Dist: myst-parser>=4; extra == "docs"
Dynamic: license-file

# CosmoFit

> **Modern Cosmological Parameter Estimation in Python**

[![PyPI](https://img.shields.io/pypi/v/cosmofit.svg)](https://pypi.org/project/cosmofit/)
[![Python](https://img.shields.io/pypi/pyversions/cosmofit.svg)](https://pypi.org/project/cosmofit/)
[![Tests](https://github.com/salihyesil59/CosmoFit/actions/workflows/tests.yml/badge.svg?branch=main)](https://github.com/salihyesil59/CosmoFit/actions/workflows/tests.yml)
[![Docs](https://img.shields.io/badge/docs-API%20reference-blue.svg)](https://salihyesil59.github.io/CosmoFit/)
[![License](https://img.shields.io/pypi/l/cosmofit.svg)](LICENSE)

**CosmoFit** is an open-source Python library for cosmological parameter estimation and Bayesian inference. It provides a modular framework for fitting cosmological models to observational data using Markov Chain Monte Carlo (MCMC) techniques.

The project is designed to make cosmological analyses simple, reproducible, and extensible while remaining flexible for research applications.

```bash
pip install cosmofit
```

---

## Features

* Modular cosmological model framework, with curvature (Omega_k) supported end-to-end (flat/open/closed E(z) and D_M(z))

  * **LCDM** -- flat/curved ΛCDM
  * **wCDM** -- constant dark-energy equation of state w0
  * **CPL** (Chevallier-Polarski-Linder) -- w(z) = w0 + wa z/(1+z)
  * **JBP** (Jassal-Bagla-Padmanabhan) -- w(z) = w0 + wa z/(1+z)^2
  * **BA** (Barboza-Alcaniz) -- w(z) = w0 + wa z(1+z)/(1+z^2)
  * **LogarithmicDE** (Efstathiou) -- w(z) = w0 + wa ln(1+z); the one w0-wa form here
    that does *not* saturate at high z, so it is the control case for asking whether a
    measured `wa` is the data or the assumed shape
  * **PEDE** (Phenomenologically Emergent Dark Energy) -- Omega_de(z) = Omega_de0
    [1 - tanh(log10(1+z))], with **no free dark-energy parameter at all**: the same
    parameter count as LCDM and a completely different expansion history, so an AIC/BIC
    comparison against LCDM is a pure comparison of fit
  * **GEDE** (Generalized EDE) -- contains both LCDM (Delta -> 0) and PEDE
    (Delta = 1, z_t = 0) as exact limits, so `Delta` measures the distance from a
    cosmological constant on a continuous scale
  * **LsCDM** (sign-switching Lambda, Akarsu et al.) -- Lambda flips sign at
    z_dagger ~ 2 (AdS below, dS above), shrinking r_d and so raising the BAO-inferred H0;
    a route to the H0 tension that late-time-only dark-energy models cannot take
  * **GCG** (Generalized Chaplygin Gas) -- unified dark matter/dark energy fluid,
    p = -A/rho^alpha
  * **IDE** (Interacting Dark Energy) -- Q = 3 xi H rho_de in closed form; changes how
    *matter* dilutes, which no w(z) parametrization does, so it has its own
    growth-of-structure signature
  * **RunningVacuum** -- Lambda(H) = c0 + 3 nu H^2, from renormalization-group running of
    the vacuum energy; one of the few extensions whose extra parameter has a *predicted*
    magnitude (|nu| ~ 10^-3) rather than an arbitrary one

* Models where acceleration comes from a modified Friedmann equation rather than any dark
  energy at all

  * **Cardassian** (modified polytropic) -- H^2 = A rho + B rho^n, a flat, matter-dominated,
    accelerating universe
  * **DGP** (Dvali-Gabadadze-Porrati braneworld, self-accelerating branch) -- gravity leaks
    into a fifth dimension; LCDM's parameter count, and a *suppressed* growth history
    (mu ~ 0.72 today) that is the real observational handle on it

* Modified-gravity models, at both background *and* growth-of-structure level (see
  [CHANGELOG.md](CHANGELOG.md) at v0.16.0 and v0.17.0 for exactly what's real vs.
  simplified in each)

  * **FQExponential** -- f(Q) gravity (symmetric teleparallel), f(Q) = Q exp(lambda Q0/Q)
  * **FTPowerLaw** -- f(T) gravity (metric teleparallel), f(T) = T + alpha T^n; n = 0 is LCDM exactly
  * **FRTLinear** -- f(R,T) gravity (linear), f(R,T) = R + 2 lambda T
  * **FRHuSawicki** -- f(R) gravity (Hu-Sawicki); background is identical to LCDM's by
    construction (stated explicitly, see [CHANGELOG.md](CHANGELOG.md) at v0.16.0), but
    growth of structure -- where this model's `f_R0`/`n` parameters actually show up --
    is now real: a scale- and time-dependent `mu(a,k)` -- the linear result, with no chameleon
    screening (which is non-linear)

* Growth of structure: every model (not just the three above) gets a linear growth factor
  D(z), growth rate f(z), and fsigma8(z)/S8 via a generic `mu(a,k)` hook on top of its own
  E(z)/dE(z)dz (`cosmology.calculators.growth.GrowthCalculator`) -- `mu = 1` (standard GR
  growth) for LCDM/wCDM/CPL/JBP/BA/GCG, a real derived `mu(a,k)` for the three
  modified-gravity models above
* **The BAO sound horizon `r_d` computed from scratch** (`Fitter(..., compute_rd=True)`)
  rather than fitted as a free nuisance parameter -- the integral
  `r_d = int c_s/H dz` with photons, massless neutrinos and massive neutrinos carrying
  their exact Fermi-Dirac energy density, validated against CAMB's `rdrag` to 5e-5. This
  is what turns BAO from a *relative* distance measurement (which constrains only
  `H0 * r_d`) into an absolute one, and it is how "BAO + BBN gives H0" works
* **σ₈ derived from the CMB rather than fitted** (`Fitter(..., derive_sigma8=True)`), when the fit
  contains a from-scratch CMB likelihood. Without it, `sigma8` is a free parameter *and* the CMB
  fixes an amplitude of its own through `ln1e10As` -- two unrelated numbers for one quantity, with
  the sampler reporting a posterior for each. It also changes what a fit can ask: with `sigma8`
  free, the S₈ measurement is absorbed (the parameter slides onto it and χ² goes to zero, which
  looks like agreement); derived, the CMB's own prediction meets the measurement and the S₈
  tension is visible
* Flexible parameter management
* Built-in observational datasets

  * Cosmic Chronometers (CC)
  * BAO (**DESI DR2 2025** or DESI DR1 2024; SDSS BOSS DR12 + eBOSS DR16 LRG/QSO;
    **low-z 6dFGS + SDSS DR7 MGS**) -- the low-z pair is independent of the other two and
    can join either; DESI and SDSS cannot be combined with each other, see the note below
  * The **SDSS BAO + full-shape consensus**: D_M/r_d, D_H/r_d *and* f*sigma8 at
    z = 0.38, 0.51, 0.698, 1.48 with the covariance between them. The same galaxies
    as the BAO-only entry above, so use one or the other -- and prefer this one over
    pairing `sdss_bao` with `fsigma8`, which covers the same galaxies while treating
    growth and geometry as uncorrelated when they are not
  * BAO as a **tabulated likelihood surface** rather than a mean and a covariance:
    **eBOSS DR16 ELG** (`D_V/r_d` at z=0.845, a 399-point curve) and **eBOSS DR16
    Lyman-alpha** (`(D_M/r_d, D_H/r_d)` at z=2.334, a 50x50 surface, the
    highest-redshift BAO here outside the CMB). eBOSS released both as grids because
    a Gaussian misrepresents them -- the ELG BAO is a 1.4-sigma detection whose
    likelihood is still rising at the low edge of the table
  * The **eBOSS DR16 ELG full-shape** analysis of the same galaxies: a 100x100x100
    grid in `(D_M/r_d, D_H/r_d, f*sigma8)`, so it constrains the **growth rate**
    alongside the geometry and carries the fsigma8/Alcock-Paczynski degeneracy
    exactly. Mutually exclusive with the BAO-only ELG entry
  * **Holographic dark energy** (Li 2004) -- the one background model here with
    no closed-form `E(z)`: `Omega_DE` comes from an ODE solved on every parameter
    change, with the future event horizon as the infrared cutoff. Its single
    parameter `c` decides whether `w` crosses below -1, as a prediction rather
    than a parametrization choice
  * **Agegraphic** and **Ricci** dark energy, the other two members of the holographic
    family -- same `rho_DE = 3c^2 M_p^2 / L^2`, different cutoff (conformal age, Ricci
    scalar). ADE has **one fewer free parameter than LCDM**: its early-time condition
    fixes `Omega_m` from `n`, so n = 2.8 *predicts* `Omega_m = 0.280`
  * Supernova (Pantheon+, DES-SN5YR, **Union3**) -- one per fit, see the note below
  * CMB, four ways: the compressed distance priors (Planck 2018 R, l_A, omega_b_h2 -- fast,
    works for every model); the **full Planck 2018 TT/TE/EE spectra** (615 `plik_lite`
    bandpowers including the two Commander low-l temperature bins, computed from scratch with
    CAMB -- slow, LCDM and w(z) models only); **Planck 2018 CMB lensing** (9 bandpowers of the
    reconstructed lensing potential, `8 <= L <= 400`), the CMB's own measurement of how much
    structure grew; and **Planck 2018 low-l EE** as the tabulated, *non-Gaussian* likelihood it
    actually is rather than a Gaussian shorthand (which is still available as the `"tau"`
    dataset: `0.0506 +- 0.0086` from lowE alone by default, or the familiar `0.0544 +- 0.0073`,
    which is the TT,TE,EE+lowE posterior and must not be combined with `planck_lite`)
  * **ACT DR6 CMB lensing** -- a second, independent lensing reconstruction, tighter than
    Planck's (2.3% on the amplitude), built on the lensing convergence rather than the potential
  * Growth rate fsigma8(z) (Gold-2018 RSD compilation, 22 points)
  * S8 weak-lensing prior (KiDS-1000 or DES Y3, Gaussian) -- don't combine the two versions,
    see the note below
  * **External single-number measurements**, entering as datasets rather than as priors so
    they show up in the chi2 breakdown and the degrees-of-freedom count: local **H0**
    (SH0ES 2022/2024, or TDCOSMO 2025 time-delay lensing -- independent of the Cepheid
    ladder), a **BBN** constraint on omega_b h^2 (Schoeneberg 2024 or Cooke 2018), and the
    Planck lowE **tau** prior

* Modular likelihood architecture
* Bayesian parameter estimation on the 2.0 core: adaptive Metropolis-Hastings stopped by the
  Gelman-Rubin `R - 1`, or emcee's ensemble stopped by the autocorrelation time, with
  getdist-format output, resuming, fast/slow parameter blocking and chains in several
  processes or under `mpirun` (`cosmofit run`, or `cosmofit.compat.sample_on_core(fitter, ...)`)
* **Bayesian evidence by nested sampling** (the core's `nested` sampler, or
  `cosmofit.compat.evidence_on_core(fitter)`; needs `pip install
  "cosmofit[evidence]"`): `ln Z` and a proper Bayes factor, for the comparisons a
  likelihood-ratio test cannot make. `LsCDM` reduces to `LCDM` only as `z_dagger -> infinity`
  and `DGP` is not nested at all, so Wilks' theorem does not apply to either -- an evidence
  ratio is defined regardless, and integrates rather than maximizes, so it charges a model
  for prior volume it does not use. Validated against an analytically integrable Gaussian.
  Read `stats.evidence`'s note on prior sensitivity before quoting one
* **Profile likelihood** (the core's `profile` sampler, or
  `cosmofit.compat.profile_on_core(fitter, "z_dagger", values)`): `chi2` minimized over every
  other parameter at each fixed value. The honest tool where Wilks fails, and where a marginal
  posterior would smooth over structure -- it is how `lscdm_mcmc.ipynb` found a 28-unit cliff
* **Tension statistics** (`stats.tension`): the `np.hypot` this repo used to write by hand,
  named and with its assumptions stated -- plus the alternatives for when they fail.
  `sample_tension` works from posterior samples with no Gaussian assumption;
  `gaussian_tension_nd` because a tension in a plane is not the larger of its projections
  (two posteriors can sit 0.35 sigma apart in each parameter and 1.74 sigma apart jointly);
  and `suspiciousness`, which divides out the prior dependence a Bayes factor carries.
  Checked against the published Hubble (4.85 sigma) and S8 (2.67 sigma) tensions, and
  against the analytic parameter-difference chi2
* **Fisher matrix** (the core's `fisher` sampler, or `cosmofit.compat.fisher_on_core(fitter)`):
  parameter errors from the curvature at the best fit, in `~2n^2` evaluations rather than a
  chain, with each step sized by the likelihood and checked by halving. What `s8_tension_cmb.ipynb` needed when every
  likelihood call is a CAMB call and a converged chain is thirteen hours
* Consolidated result object (`fitter.result`): best-fit + MCMC posterior in one printable,
  JSON-serializable snapshot
* Saved chains (an `output` prefix): written as they are sampled, in getdist's format, and
  read back instead of re-sampled next time -- a run that had met its stopping rule costs not
  a single new step, an interrupted one continues -- see [Saved Chains](#saved-chains)
* Custom models (`define_model`): fit a brand-new, not-in-the-library `E(z)` -- with its own
  extra parameters -- against every built-in dataset/likelihood/MCMC, no library changes needed
* **Models from an action** (`cosmofit.theory`, `pip install -e ".[theory]"`): give a
  gravitational action on an FLRW metric and the library does the variational calculus --
  reduces it to a point-like Lagrangian, varies the lapse for the Friedmann constraint,
  solves or integrates it, and hands back an ordinary model every dataset already works on

  * **Undeformed gravity in three sectors** -- `R`, `T` (torsion), `Q` (non-metricity);
    the sign conventions are fixed by requiring each to reproduce General Relativity
    exactly, and the tests assert it rather than trusting the convention
  * **`f(T)` / `f(Q)`** -- solved pointwise, transcendental constraints included: the
    exponential `f(Q)` model's Lambert-`W` Friedmann equation is rederived from the action
  * **Quintessence and k-essence**, any `L(X, phi)` and any number of fields -- integrated
    forwards from `z_init`, with `E(0) = 1` as a shooting condition. Validated against
    *both* Copeland-Liddle-Wands attractors to five decimals
  * **Scalar-tensor `F(phi) R`** -- the field sets the strength of gravity, bringing the
    `3 H dF/dt` term into the Friedmann equation and a moving `G_eff/G_N` into the growth
  * **A general `f(R)`** -- fourth-order, reduced by promoting `R` to an independent
    variable held to its geometric value by a Lagrange multiplier
* Graphical interface (`cosmofit gui`, `pip install "cosmofit[gui]"`): tick datasets, configure
  one or more models (built-in or your own), edit free parameters, and run the MCMC fit(s) + plots
  with one click, no code -- compare models side by side statistically (AIC/BIC/a likelihood-ratio
  test) and on the same figures, all from the browser, with every figure downloadable as SVG, PNG,
  or PDF. It is also **self-explaining**: every dataset and model carries a note saying what it
  measures, over what redshift range, what it constrains and where it comes from; six presets
  configure a whole analysis in one click; the parameter table shows only the parameters *this*
  fit actually uses; and the app warns about combinations that will not work (or, worse, will
  quietly produce a posterior for something the data cannot constrain) before you run them
* Covariance matrix support, including precision (inverse-covariance) matrices as shipped
  directly by some data releases
* Dedicated plotting module (`fitter.plots`): MCMC chain/corner plots, Hubble diagram, H(z)
  diagram, BAO distance plot, Planck pull plot, w(z) evolution, deceleration parameter, and the
  w0-wa dark-energy plane (2D credible contours over the phantom/quintessence/quintom regions,
  with ΛCDM marked) -- the headline figure of the DESI evolving-dark-energy results
* Publication-ready figure text: every axis, title, legend and tick renders as LaTeX
  (`$\Omega_b$`, `$r_d$`, `$\Lambda$CDM`, `$f(R)$ Hu-Sawicki`) rather than as Python
  identifiers, corner-plot titles quote each parameter to its own precision, and each label
  states the quantity actually plotted (`$m_B$` with `$M_B$` marginalized, `$D_V/r_d$` with the
  `r_d` the code really divides by)
* Model comparison plots (`fitter.plots.compare_*`): any of the figures above, overlaying this
  fit's curve with one or more other models' curves on the same data/axes -- defaults to this
  model vs. a quick LCDM reference, or an arbitrary N-model comparison via `other_fits=[...]`
* Derived-quantity posteriors (`cosmofit.stats.derived`): the acceleration transition
  redshift `z_t` (where q(z) changes sign) and today's deceleration parameter `q0`, with
  proper error bars -- every posterior sample is pushed back through the model's own
  E(z)/dE(z)dz rather than the quantity being evaluated once at the best fit, and it works
  for every model, not just CPL (whose w0/wa-specific diagnostics live in
  `stats.cpl_diagnostics`)

---

## Installation

```bash
pip install cosmofit
```

Everything in the core works with no extras at all -- twenty-one
bundled datasets, twenty models, the samplers and every figure. The
four optional extras are each there for one thing:

```bash
pip install "cosmofit[cmb]"       # CAMB: the CMB computed rather than compressed
pip install "cosmofit[theory]"    # sympy: derive a model from its action
pip install "cosmofit[evidence]"  # dynesty: nested sampling, Bayes factors
pip install "cosmofit[speed]"     # numba: ~1.7x on growth-heavy fits, nothing else
```

Everything is importable from the top-level package:

```python
from cosmofit import CPL, Fitter
```

To work on the library itself, clone it instead:

```bash
git clone https://github.com/salihyesil59/cosmofit.git
cd CosmoFit
pip install -e ".[dev,cmb,theory,evidence,speed,docs]"
```

See [CONTRIBUTING.md](CONTRIBUTING.md).

---

## Project Structure

CosmoFit uses a standard `src` layout. The importable package
(`cosmofit`, the same name as on PyPI) is separate from the repository
root, which keeps the package import path unambiguous and matches how
the library is meant to be used -- `pip install`ed, then imported by
name, not by adding the repo root to `sys.path`.

Before 2.0 the package was imported as `CosmoFit`. That name still
works for one release -- `import CosmoFit` and every
`from CosmoFit.<module> import ...` resolve to the same `cosmofit`
modules -- with a `DeprecationWarning` saying to switch.

```text
CosmoFit/
├── src/cosmofit/
│   ├── __init__.py    # unified public API: from cosmofit import ...
│   ├── core/          # 2.0 core: theories, likelihoods, parameters, priors, YAML inputs
│   ├── samplers/      # 2.0 samplers: evaluate, minimize, profile, fisher, mcmc, emcee, nested
│   ├── theories/      # 2.0 theories: background (radiation, neutrinos), early universe
│   ├── cosmology/     # models (LCDM, wCDM, CPL), parameters, distances, background
│   ├── data/          # dataset loaders + bundled CC/DESI/Pantheon+/Planck data files
│   ├── likelihoods/   # per-dataset chi2/likelihood classes + joint likelihood
│   ├── stats/         # Fitter (MCMC/best-fit), priors, posterior, saved chains, model comparison
│   ├── theory/        # build a model from an action: reduce, vary, solve (optional: sympy)
│   ├── gui/           # the Streamlit GUI: `cosmofit gui`
│   └── plots/         # FitPlotter -- every figure, attached as fitter.plots
├── examples/          # example notebooks, in five sections
├── tests/             # the test suite
├── app/               # the GUI's entry point in a checkout (the page is cosmofit.gui)
├── CHANGELOG.md       # the release history
└── pyproject.toml
```

The subpackages remain directly importable too
(`cosmofit.stats.model_comparison`, `cosmofit.data.loader`, ...) for
anything not re-exported at the top level -- `from cosmofit import ...`
is a convenience layer over them, not a replacement.

---

## Example Notebooks

Seventeen notebooks under [`examples/`](examples/), organised by what you are
trying to do — see [`examples/README.md`](examples/README.md) for the full
index. All are Colab-ready: click a badge and *Runtime → Run all*, no local
setup.

| | |
|---|---|
| **[01 · Getting started](examples/01-getting-started)** | [`quickstart`](examples/01-getting-started/quickstart.ipynb) — nothing to a real MCMC posterior in a couple of minutes; [`dataset_zoo`](examples/01-getting-started/dataset_zoo.ipynb) — all 21 datasets, the three non-Gaussian ones, and the sixteen pairs not to combine; [`cmb_from_scratch`](examples/01-getting-started/cmb_from_scratch.ipynb) — the CMB computed rather than compressed, and what that costs. |
| **[02 · Models](examples/02-models)** | [`model_zoo_comparison`](examples/02-models/model_zoo_comparison.ipynb) — six dark-energy parametrizations head to head; [`modified_gravity_growth`](examples/02-models/modified_gravity_growth.ipynb) — f(Q)/f(R,T)/f(R) with `mu(a,k)` and growth data; [`holographic_family`](examples/02-models/holographic_family.ipynb) — HDE/ADE/RDE against published constraints. |
| **[03 · Building models](examples/03-building-models)** | [`custom_models`](examples/03-building-models/custom_models.ipynb) — three routes to your own `E(z)`; [`models_from_an_action`](examples/03-building-models/models_from_an_action.ipynb) — give an *action* and let the library derive `E(z)`; [`scalar_field_models`](examples/03-building-models/scalar_field_models.ipynb) — quintessence, k-essence and scalar-tensor gravity, integrated. |
| **[04 · Inference](examples/04-inference)** | [`evidence_and_model_selection`](examples/04-inference/evidence_and_model_selection.ipynb) — AIC/BIC, likelihood-ratio, and Bayesian evidence; [`profile_likelihood_and_fisher`](examples/04-inference/profile_likelihood_and_fisher.ipynb) — three ways to an error bar; [`tension_statistics`](examples/04-inference/tension_statistics.ipynb) — four definitions of "they disagree at 4σ". |
| **[05 · Case studies](examples/05-case-studies)** | [`cpl_mcmc_analysis`](examples/05-case-studies/cpl_mcmc_analysis.ipynb) and [`cpl_mcmc_tfd42`](examples/05-case-studies/cpl_mcmc_tfd42.ipynb) — the CPL deep dive and its publication-scale variant; [`lscdm_mcmc`](examples/05-case-studies/lscdm_mcmc.ipynb) — does ΛsCDM still relieve the H₀ tension?; [`s8_tension_cmb`](examples/05-case-studies/s8_tension_cmb.ipynb) — the S₈ tension from a from-scratch Planck spectrum; [`dark_energy_evidence_audit`](examples/05-case-studies/dark_energy_evidence_audit.ipynb) — how much of the dark-energy evidence is a choice, across 20 combinations. |

---

## Quick Example

```python
from cosmofit import CPL, Fitter

fitter = Fitter(
    model=CPL,
    datasets=["cc", "desi", "pantheon", "planck"],
    free_params=[
        "H0",
        "Omega_m",
        "w0",
        "wa",
        "rd",
    ],
    initial={
        "H0": 67.4,
        "Omega_m": 0.315,
        "w0": -1.0,
        "wa": 0.0,
        "rd": 147.1,
        "Omega_b": 0.0493,
    },
)

fitter.best_fit()

# The 2.0 core samples it: adaptive Metropolis-Hastings, stopped once
# its chains agree (R - 1 < 0.01), written to chains/cpl.* in getdist's
# format. The chains become the fitter's own, so everything below reads
# them.
from cosmofit.compat import sample_on_core

sample_on_core(fitter, "mcmc", {"Rminus1_stop": 0.01}, output="chains/cpl")

fitter.summary()
fitter.convergence()

print(fitter.result)

fitter.plots.corner()
fitter.plots.hubble_diagram()
fitter.plots.w_of_z()
```

The same fit with no Python at all -- the input the 2.0 core runs, which
`fitter.to_info()` writes for you:

```bash
cosmofit run examples/yaml/cpl_bao_sn_cmb.yaml
```

---

## Saved Chains

The sampling is the expensive part of a fit; everything built on top of it -- summaries,
convergence diagnostics, corner plots, derived quantities, model comparison -- is seconds.
Give the 2.0 core an output prefix and the chains go to disk as they are sampled, so the
expensive part happens once:

```python
from cosmofit.compat import sample_on_core

sample_on_core(fitter, "mcmc", {"Rminus1_stop": 0.01}, output="chains/cpl")
```

Run that same code again -- next week, in a new session, after adding three more plots below
it -- and nothing is re-sampled: the saved chains are read back, checked against their
stopping rule, and `summary()`, `convergence()`, `best_fit()` and every `plots.*` figure work
off them exactly as before.

| what you do | what happens |
|---|---|
| run it again unchanged | the chains had met their rule: read back, not a single new step |
| tighten `Rminus1_stop`, or raise `max_samples` | the saved chains continue from where they stopped |
| Ctrl-C halfway through | every check's rows are already on disk; run again to carry on |
| change the model, datasets, free parameters, or priors | refused, naming what differs -- samples from two different posteriors must never be merged |

The files are getdist's format -- `chains/cpl.1.txt` to `.4.txt` (one per chain), with
`.paramnames`, `.ranges`, the learned `.covmat`, the input as given and as run
(`.input.yaml`, `.updated.yaml`) and a `.progress` log -- so getdist, or anything else that
reads CosmoMC/cobaya chains, opens them directly:

```python
from getdist import loadMCSamples

samples = loadMCSamples("chains/cpl", settings={"ignore_rows": 0.3})
```

The same holds for a run started from the command line (`cosmofit run input.yaml` with
`output: chains/cpl`), and `resume: true` / `force: true` in an input say what to do with
output already there. A finished run can be reweighted by another likelihood without
sampling again -- a `post:` block, see `cosmofit.core.post`.

### Chains saved by 1.x

`Fitter.run_mcmc(save="chains/cpl.h5")` -- emcee, written to HDF5 -- still works for one
release, with a `DeprecationWarning`, and the files it wrote stay readable:
`Fitter.from_chain("chains/cpl.h5")` rebuilds the whole fit from what the file records, and
`cosmofit.stats.chains.open_chain` reads the samples without evaluating a likelihood.

---

## The w0-wa Plane

`fitter.plots.w0_wa_plane()` draws the figure the recent
evolving-dark-energy literature is built around: the 2D posterior of the CPL parameters over the
four dark-energy regions of the (w0, wa) plane, with ΛCDM marked at (-1, 0).

```python
from cosmofit import CPL, Fitter

fit = Fitter(
    model=CPL,
    datasets=["cc", "desi", "pantheon"],
    free_params=["H0", "Omega_m", "w0", "wa"],
    initial={"H0": 67.4, "Omega_m": 0.315, "w0": -1.0, "wa": 0.0, "rd": 147.1},
)
from cosmofit.compat import sample_on_core

sample_on_core(fit, "mcmc", {"Rminus1_stop": 0.005}, output="chains/cpl")

fit.plots.w0_wa_plane()
```

The plane is cut by two lines -- `w0 = -1` (the equation of state today) and `wa = -1 - w0`
(its high-z limit `w0 + wa`) -- into the four cases the classification is made of:

| region | w(z) | meaning |
|---|---|---|
| **Phantom** | below -1 at all times | ρ grows with expansion; a Big Rip future |
| **Quintessence** | above -1 at all times | reachable by an ordinary canonical scalar field |
| **Quintom-A** | quintessence in the past, phantom today | crosses w = -1 |
| **Quintom-B** | phantom in the past, quintessence today | crosses w = -1 the other way |

The crossing is the physically loaded part: no single canonical (or single phantom) scalar field
can cross `w = -1` at all, so a posterior sitting in either quintom region calls for something
more -- two fields, a non-canonical kinetic term, or modified gravity.

Useful options:

```python
fit.plots.w0_wa_plane(show_fractions=True)          # posterior probability of each region
fit.plots.w0_wa_plane(levels=(0.68, 0.95, 0.997))   # add a third contour
fit.plots.w0_wa_plane(w0_range=(-1.3, 0.0), wa_range=(-3, 1))
fit.plots.w0_wa_plane(bins=120, smooth=1.0)         # finer density estimate
```

`levels` are 2D credible **probabilities** (the smallest area holding that fraction of the
samples), not sigmas -- in two dimensions the familiar "1σ"/"2σ" contours enclose only 39% and
86.5%, and conflating the two is the usual way this figure gets over-read.

The same numbers without the picture:

```python
from cosmofit.stats import cpl_diagnostics

cpl_diagnostics.region_fractions(*[fit.samples_dict()[k] for k in ("w0", "wa")])
# {'phantom': 0.004, 'quintessence': 0.221, 'quintom-a': 0.012, 'quintom-b': 0.763}

cpl_diagnostics.classify_region(-0.85, -0.6)   # 'quintom-b'
```

To overlay several dataset combinations -- the "what does adding this data do to w0-wa?"
figure -- every fit needs its own chain, and then:

```python
fit_desi.plots.compare_w0_wa_plane(other_fits=[fit_desi_sn, fit_desi_sn_cmb])
```

Each contour set labels itself with its own dataset combination (`CC + DESI + Pantheon+`);
pass `labels=[...]` to override.

Both methods need `w0` and `wa` to be free parameters of a completed MCMC, and a model whose
w(z) really does tend to `w0 + wa` at high z (CPL, BA). For JBP -- whose w(z) returns to `w0`
instead -- the diagonal boundary would not be that model's own limit, so they refuse rather than
mislabel the regions; use `plots.w_of_z()` there.

---

## Model Comparison Plots

Every `fitter.plots.*` figure above has a `compare_*` counterpart that overlays this fit's curve
with one or more other models' curves on the same data/axes -- the "model A vs model B" figures
cosmology papers use (`compare_hz`, `compare_deceleration`, `compare_w_of_z`,
`compare_hubble_diagram`, `compare_des_hubble_diagram`, `compare_bao_distances`,
`compare_sdss_bao_distances`):

```python
fitter.plots.compare_hz()             # vs. a quick LCDM reference, built automatically
fitter.plots.compare_deceleration()   # ditto, each curve's own transition redshift z_t marked
```

Called with no arguments, `other_fits` defaults to a quick best-fit-only LCDM reference (no MCMC
-- just enough for a comparison curve) built from this fit's own datasets, unless this fit already
*is* LCDM. Pass an already-fit `Fitter` (or a list of them, for more than two models at once) to
compare against instead:

```python
fitter.plots.compare_hz(other_fits=fit_lcdm, labels=["CPL", "LCDM"])
fitter.plots.compare_hz(other_fits=[fit_lcdm, fit_wcdm], labels=["CPL", "LCDM", "WCDM"])
```

This works for any model CosmoFit knows about, built-in or [custom](#custom-models).

---

## Custom Models

Testing a model that isn't in the literature -- and isn't one of
CosmoFit's twenty built-ins -- doesn't require touching the library's
internals. `define_model` builds a usable model from a single
`E(z)` function (that alone is enough to fit against every dataset
and produce every plot except `w_of_z()`/`deceleration()`); any new
parameter it needs beyond the standard set (`H0`, `Omega_m`,
`Omega_k`, `w0`, `wa`, ...) is declared right there, with a default
and prior bounds:

```python
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", "pantheon"],
    free_params=["H0", "Omega_m", "w0", "beta"],
    initial={"H0": 67.4, "Omega_m": 0.315, "w0": -1.0, "beta": 0.0},
)
fit.best_fit()
sample_on_core(fit, "mcmc", {"Rminus1_stop": 0.01})   # from cosmofit.compat
print(fit.result)
fit.plots.corner()
```

For more control, the same mechanism is available by subclassing
`Cosmology` directly and declaring `EXTRA_PARAMS` on the class --
`define_model` is a thin convenience wrapper around exactly this:

```python
from cosmofit import Cosmology

class MyModel(Cosmology):
    EXTRA_PARAMS = {"beta": {"default": 0.0, "bounds": (-2.0, 2.0)}}

    def _dark(self, z):
        return (1 - self.Omega_m) * (1 + z) ** (3 * (1 + self.w0)) * (1 + self.beta * z)

    def E(self, z):
        z = np.asarray(z, dtype=float)
        return np.sqrt(self.Omega_m * (1 + z) ** 3 + self._dark(z))

    def dEdz(self, z, h=1e-5):
        z = np.asarray(z, dtype=float)
        return (self.E(z + h) - self.E(z - h)) / (2 * h)      # or by hand
```

A direct subclass **must** define `dEdz`, and from the moment it is
constructed: the distance integrator interpolates `1/E(z)` with a Hermite
spline built from that derivative, which is what makes it exact to fourth
order instead of second. `define_model` installs a central-difference
fallback for you; subclassing does not, so that a model can supply the exact
derivative -- as every built-in one does -- rather than silently give up
about four orders of magnitude of accuracy in every distance.

---

## Models From an Action

`define_model` still asks for `E(z)`, which means somebody has
already done the variational calculus by hand. `cosmofit.theory`
takes the other end: give it a gravitational action on an FLRW
metric and it derives `E(z)` itself.

The method is the standard minisuperspace one. FLRW is written
with an explicit lapse `N(t)`, the action is reduced to a
point-like Lagrangian in `a`, `adot` and `N`, and varying the
lapse produces the Friedmann *constraint* -- the lapse is a
non-dynamical gauge degree of freedom, which is exactly why its
variation gives a constraint rather than an evolution equation.
Setting `N = 1` afterwards recovers the familiar form. Writing
`N = 1` from the start would lose the equation altogether.

```python
from cosmofit import Fitter
from cosmofit.theory import Action

# Power-law f(T) gravity (Bengochea & Ferraro 2009).
model = Action(
    "T + A0*(-T)**b",
    geometry="teleparallel",
    params={
        "A0": {"default": -4.2, "bounds": (-30.0, 0.0)},
        "b": {"default": 0.0, "bounds": (-2.0, 0.9), "label": r"$b$"},
    },
    closure="A0",           # fixed by E(0) = 1, not fit
    growth="quasi_static",  # mu = 1/f', for fsigma8/s8
).build("PowerLawFT")

fit = Fitter(
    model=model,
    datasets=["cc", "desi"],
    free_params=["H0", "Omega_m", "b", "rd"],
    initial={"H0": 70.0, "Omega_m": 0.3, "b": 0.0, "rd": 147.0},
)
fit.best_fit()
```

What comes back is an ordinary `Cosmology` subclass. Every
dataset, likelihood, sampler and plot in the library works on it
unchanged.

`Action.constraint()` returns the derived Friedmann equation
symbolically, if the derivation is what you wanted rather than the
fit.

### What it checks, and what it refuses

The three geometry scalars carry different signs across the
literature, and picking one wrongly inverts every modification
built on it while leaving nothing downstream to complain. They are
fixed here by requiring that an undeformed `f` reproduce General
Relativity exactly, and the test suite asserts that in all three
sectors rather than trusting the convention.

* An action that does not satisfy `E(0) = 1` is **refused**. It
  would predict every distance wrong by a constant factor without
  looking broken. `closure=` names the parameter that condition
  fixes -- in Lambda-CDM, `"R - 2*Lam"` with `closure="Lam"` is
  what makes `Lam = 3 (1 - Omega_m - Omega_k)`.
* `growth="quasi_static"` is **opt-in**. `mu = 1/f'` is a
  statement about perturbations, which a background action does
  not by itself determine.

### A general `f(R)`

Fourth-order, and it gets its own reduction — applied
automatically, with nothing to ask for:

```python
model = Action(
    "R - 2*Lam + alpha_fr*R**2",
    params={"Lam":      {"default": 2.1, "bounds": (0.0, 6.0)},
            "alpha_fr": {"default": 1e-3, "bounds": (1e-6, 1.0)}},
).build("Starobinsky")
```

The ordinary reduction removes the `addot` in the Einstein–Hilbert
term by integrating by parts, which is legitimate only while the
Lagrangian is *linear* in it. For a general `f(R)` it is not — the
term that would be dropped is not a total derivative. Promoting `R`
to an **independent variable**, held to its geometric value by a
Lagrange multiplier, gives

```
L = (1/2) N a^3 [ f(R) - f'(R) R + f'(R) R_geom ]
```

which is linear in `addot` again, at the cost of one extra dynamical
variable. That variable is the theory's fourth order made visible,
and it appears as a parameter: **`R_0`**, the Ricci scalar today —
an initial condition General Relativity does not have. There is no
`closure` here, because `E(0) = 1` holds by construction.

**Checked three ways.** The derived Friedmann constraint equals the
textbook `3 f_R H² = (f_R R − f)/2 − 3H d(f_R)/dt + rho`
symbolically. `alpha -> 0` reproduces ΛCDM, and does so smoothly:
over `z <= 5` the departure falls 6.4e-01 → 1.8e-01 → 4.2e-03 as
`alpha` goes 1e-1 → 1e-3 → 1e-5. And the accuracy measure is the
**third** equation of
motion, from varying `a` — the one the integration never uses,
which holds by the Bianchi identity and comes out at 1e-15.

Unlike a scalar field, this integrates **backwards** from today, and
the difference is measured rather than assumed: a `1e-8` kick to
`R_0` moves `E(z)` by less than that out to `z = 1100`, where a
scalar field would have run away.

#### Which direction the integration runs

Backwards from today is the default, and it is well conditioned only
while the scalaron's oscillating mode does not grow into the past.
That holds for `R + alpha R^2`. It fails for the whole
"disappearing cosmological constant" family — Hu–Sawicki,
Starobinsky 2007, Tsujikawa, the arctan models — where `f_RR` falls
steeply with `R`: measured on the arctan model, no achievable
precision in `R_0` gets past `z ≈ 1.2`.

Those integrate the other way:

```python
model = Action(
    "R - (4*Lam/pi)*atan(R/Rw)",          # arXiv:1601.07928
    params={"Lam": {"default": 2.25, "bounds": (1.0, 4.0)},
            "Rw":  {"default": 1.0,  "bounds": (0.05, 20.0)}},
    closure="Lam",
    background="forward",
    growth="quasi_static",
).build("ArctanFR")
```

Forwards the oscillating mode **decays**, so a slightly wrong initial
condition relaxes onto the attractor instead of leaving it. Two
things swap over. `R_0` stops being a parameter and becomes derived
— on the attractor the curvature today is not free — and a `closure`
becomes *required*, because `E(0) = 1` is now a condition to satisfy
rather than the point one starts from.

Neither direction is universally better, and neither is guessed:
which one a theory needs is a property of the theory.

Below the starting redshift the history is continued analytically,
because forward integration cannot start deep enough for the growth
ODE and this family is General Relativity there anyway. That
assumption is checked rather than trusted — moving the junction from
`z = 20` to `z = 60` does not move `E(z)` or `fσ8`, and a model that
has *not* returned to matter domination where it starts is refused
rather than spliced.

#### Growth, which is where `f(R)` actually differs

`growth="quasi_static"` works here too, but it cannot be the
`mu = 1/f'` of the teleparallel sectors: `f(R)`'s scalaron has a
Compton wavelength, so the coupling is **scale-dependent**,

```
mu = (1/f_R) (1 + 4m)/(1 + 3m),    m = (k/a)^2 f_RR/f_R
```

with `R` read off the model's own integrated background. Far
outside the Compton wavelength this is `1/f_R`; far inside it is
`4/(3 f_R)`, the extra third being the scalar fifth force. `mu(a, k)`
takes `k` in h/Mpc.

The algebra is not the risky part — the GR-05 notebook in the
[companion toolkit](https://github.com/salihyesil59/wljs-gr-toolkit)
derives `G_eff` from the perturbed field equations and its result is
*identically* the expression above, checked by subtracting the two
symbolically and getting zero. The **units** are the risky part:
`m` is dimensionless only if `(k/a)^2` and `f_RR` are in matching
ones, and getting that wrong moves the Compton wavelength without
making anything raise. So it is held against `FRHuSawicki`, which
performs the same conversion through independently written code, and
the two agree to 2e-16 across three scale factors and three
wavenumbers.

This is the **linear** result. Chameleon screening is non-linear and
is not in it, so where screening matters this overstates the
departure from GR.

#### Is the theory worth fitting at all?

A free function can be tuned to fit almost anything, so what
separates a proposal from a curve is whether its perturbations
behave. Two conditions decide it, and a compiled `f(R)` checks both
along its own background:

```python
model.viability()
# {"ok": False, "failed": ["f_RR"], "reasons": [...]}
```

* **`f_R > 0`** — the effective gravitational coupling. Negative
  makes the graviton a ghost; zero makes the coupling diverge.
* **`f_RR >= 0`** — the scalaron's mass squared goes as `1/f_RR`, so
  negative makes it tachyonic. That is the Dolgov–Kawasaki
  instability, and it is fast enough that such a background would
  not survive to be observed. Zero is *not* flagged: that is the
  general-relativity limit, where the scalaron is absent rather than
  sick.

`model.scalaron(z)` returns the pair directly.

This matters because a non-viable `f(R)` does not announce itself.
Flip the sign of `alpha` in the example above and `E(z)` stays
perfectly smooth and finite — nothing raises, nothing looks wrong,
and the model will happily be fitted. What does break is `mu`: with
`f_RR < 0` the denominator `1 + 3m` passes through zero, and beyond
that pole the formula still evaluates to a finite, plausible number.
It is not returned. `mu` refuses and names the condition, on the
same principle as the Boltzmann guard — a number that means nothing
is worse than an error.

### Scalar fields

An action can carry dynamical fields, in which case the expansion
history is integrated rather than solved pointwise:

```python
model = Action(
    "R",
    fields={"phi": "X - V0*exp(-lam*phi)"},   # any L(X, phi)
    params={"V0":  {"default": 2.1, "bounds": (0.05, 50.0)},
            "lam": {"default": 0.5, "bounds": (0.0, 1.7)}},
    closure="V0",
).build("ExponentialQuintessence")
```

More than one field is fine, and a field's name is also in scope in
`gravity` -- which is how **scalar-tensor gravity** is written:

```python
Action(
    "(1 + xi*phi**2)*R",                       # F(phi) R, not R + phi
    fields={"phi": "X - V0"},
    params={"xi": {"default": 0.02, "bounds": (-0.5, 0.5)},
            "V0": {"default": 2.1, "bounds": (0.05, 20.0)}},
    closure="V0",
)
```

That is a different theory from a field sitting on top of General
Relativity: the field sets the strength of gravity, and the
Friedmann equation gains the `3 H dF/dt` term. The derivation
reproduces `3 F H^2 + 3 H dF/dt = rho` symbolically, and `xi = 0`
returns ΛCDM to 1e-9.

It also changes how structure grows, and that is easy to get
silently wrong. `growth="quasi_static"` gives the scalar-tensor
`mu = G_eff/G_N` of Boisseau, Esposito-Farèse, Polarski &
Starobinsky (2000), evaluated on the field's own solution so that
it moves as the field rolls. Left at the default `"gr"`, `mu` is 1
-- correct for a *minimally* coupled field, and wrong here -- so
`Fitter` warns if such a model meets growth data.

One more thing to know: a field's history starts at `z_init`
(default 3000), and the growth ODE starts at `z = 9999`. Fitting
`fsigma8` therefore needs `Action(z_init=20000)` or similar, and
saying so is now an error rather than an infinite chi-squared.

Each field adds two parameters, `phi_i` and `dphi_i` -- its value
and `dphi/dN` at `z_init` (default 3000, and the earliest
redshift the model can be evaluated at).

**They are set early, not today, and that is the whole design.**
Setting them at `a = 1` would make `E(0) = 1` algebraic and skip
the shooting entirely. It is also wrong: integrating backwards
from a field at rest today, the Hubble friction that damps the
field forwards becomes anti-friction, the generic past solution
is kinetic-dominated, and `rho_phi` grows as `a^-6`. For
exponential quintessence normalized to today's dark-energy
density that gives `E(2.5) = 9.3` where ΛCDM gives 3.7. Nothing
fails -- it is the correct past of those initial conditions, and
those initial conditions are not a universe.

So the state is given where a quintessence model actually gives
it, early and typically frozen, and `E(0) = 1` becomes a shooting
condition on `closure`. Forwards is also the numerically stable
direction.

Validated against Copeland, Liddle & Wands (1998), whose
late-time attractors are exact functions of the slope alone. On a
matter background there are two, selected by `lambda^2` against 3:

| | `w` | predicted | `Omega_phi` | predicted |
|---|---|---|---|---|
| λ = 0.5 | −0.916667 | −0.916667 | 1.00000 | 1 |
| λ = 1.0 | −0.666667 | −0.666667 | 1.00000 | 1 |
| λ = 1.9 | −0.000035 | 0 | 0.83097 | 0.83102 |
| λ = 2.0 | +0.000035 | 0 | 0.75006 | 0.75000 |

The first two are the field-dominated attractor
(`w = -1 + lambda^2/3`); the last two are the scaling attractor,
where the field tracks matter at a fixed fraction `3/lambda^2` --
which is why a steep potential cannot accelerate. Ask such a
model for more dark energy than its attractor allows and `H(0)`
saturates below 1 with no potential scale reaching it; that
arrives as a statement about the model rather than as a stalled
integrator.

A constant potential pins the other end: it reproduces `LCDM` to
1e-10 in `E(z)`, with `w = -1` and `Omega_de` constant to machine
precision at every redshift.

Cost is about 37 ms per parameter set -- a shooting solve plus an
integration, against ~150 µs for a closed-form model. `w` and
`Omega_de` are read off the field's own Lagrangian
(`p = L`, `rho = 2 X L_X - L`) rather than by subtracting the
fluids from `E(z)^2`, which at z = 2000 would be 2.4e9 minus
2.4e9 to get 0.7.

### Validation

Because this derives what the rest of the library has typed in by
hand, it can be checked against it:

| action | reproduces | agreement |
|---|---|---|
| `R - 2*Lam` | `LCDM`, curvature included | `E(z)`, `dE/dz` to 1e-16 |
| `Q*exp(lam*Q0/Q)` | `FQExponential` | constraint identical; `lam`, `E(z)`, `mu` to 1e-13 |
| `T + A0*(-T)**b`, `b = 0` | `LCDM`, via a different sector | `E(z)` to 1e-13 |
| `T + A0*(-T)**b`, any `b` | `FTPowerLaw` | `E(z)` and `mu` to 1e-15, across two opposite sign conventions for `T` |

The `f(Q)` case is the demanding one. Its constraint is
transcendental -- the hand-written model inverts a Lambert `W` to
solve it -- so the derivation has to produce
`(E^2 - 2 lam) exp(lam/E^2) = Omega_m (1+z)^3` from the action
alone, and the solver has to land on the same branch `W_0` picks
rather than on any root of that equation.

Installed as an extra, since nothing else in the library needs
sympy:

```bash
pip install "cosmofit[theory]"
```

---

## Graphical Interface

For datasets/models/parameters by clicking rather than coding, a
[Streamlit](https://streamlit.io) app (the `cosmofit.gui` package) sits
on top of the same library. It builds a `Fitter` for its checks and
figures, and runs the calculations on the 2.0 core
(`cosmofit.compat`): the adaptive Metropolis-Hastings sampler, stopped
by R − 1, or emcee's ensemble, and the core's profile, Fisher and
nested samplers.

```bash
pip install "cosmofit[gui]"
cosmofit gui
```

Anything after `cosmofit gui` goes to `streamlit run`
(`cosmofit gui --server.port 8600`). From a checkout,
`streamlit run app/streamlit_app.py` runs the same page.

Or, without touching a terminal: double-click `run_gui.sh` (Linux/macOS)
or `run_gui.bat` (Windows) in the repository root. Either installs
whatever's missing on first run and then opens the app in your
browser; safe to double-click again any time to relaunch it.

Chains are saved and reused (the "Saved chains" box in the sidebar,
on by default), in getdist's format, one folder per distinct
configuration: add a second model to compare against and the first
one comes back instantly instead of being re-sampled, tighten R − 1 or
raise Max steps and the saved chains continue rather than restart,
close the app and it's all still there. Changing the model, its
datasets, free parameters, priors, the sampler or the seed makes a
different folder, so nothing is reused when it shouldn't be.

It lets you: tick which built-in datasets to fit; pick one of the
seventeen built-in models, or write a custom one directly as an
`E(z)` expression (e.g. `sqrt(Omega_m*(1+z)**3 + ...)`, using the
same mechanism as `model_from_expression()` below) with your own
extra parameters; tick which parameters are free and edit their
initial values/bounds in a table; and, with one click, run the MCMC
fit and render whichever result plots apply to your model/dataset
choice. Results (best fit, posterior summary, convergence) can be
downloaded as JSON via `FitResult.save_json()`.

**It also explains itself.** Fourteen datasets and seventeen models
is a lot to face cold, and a wrong combination does not announce
itself -- it produces a perfectly ordinary-looking posterior. So:

* **Every dataset carries a note** -- what it measures, over what
  redshift range, how many points, what it actually constrains, and
  the paper it comes from. They are grouped by probe (expansion rate,
  BAO, supernovae, CMB, growth, external measurements), because a fit
  is normally built by taking one from each family rather than by
  ticking everything.
* **Every model carries a note** -- what it is, what its extra
  parameters mean, and *which parameter values reduce it to ΛCDM*,
  which is the number an AIC/BIC comparison is measuring the distance
  from. Capability badges say whether it has its own `w(z)`, whether
  it modifies growth, and whether the full CMB spectra can be
  computed for it.
* **Six presets** configure a whole analysis in one click --
  including "DESI DR2 + BBN → H₀ without the CMB", which sets the
  datasets, picks DR2, and turns on the computed sound horizon
  together, since none of those three is useful without the others.
* **The parameter table shows only what this fit uses.** The shared
  container carries every parameter any model needs; for an LCDM fit
  against CC and BAO all but four are inert. Relevance depends on the
  datasets too -- `rd` appears only with BAO, `sigma8` only with
  growth data, `n_s`/`tau` only with the CMB spectra.
* **It warns before you run, not after.** A model the CMB spectra
  cannot be computed for is an error you would otherwise meet as a
  stack trace minutes in. The quieter ones matter more: a
  modified-gravity model with no growth data ticked, a model whose
  own parameters are left fixed, a free `sigma8` that nothing in the
  fit constrains. Each of those runs, converges, and returns a
  posterior that looks like a measurement.
* **A χ² breakdown per dataset** on the results tab, which is what
  turns a total χ² into "the local H₀ measurement is contributing 24
  of it, on one data point" -- the entire content of a tension.
* **Derived quantities** (`q₀`, `z_t`, and `r_d` computed from the
  densities) with real error bars, pushed sample by sample back
  through the model's own `E(z)`.
* A **📖 Guide** panel listing every dataset and model in one
  browsable table, with the conflict rules and a short "how to use
  this" walkthrough.

The custom-model expression box evaluates with `eval()` but with
Python's builtins removed and only whitelisted `numpy` math plus the
model's own parameter names reachable -- appropriate for running
locally on your own machine, not for a publicly hosted, multi-tenant
deployment.

---

## Project Status

**v2.0.0** turns the library into a set of components: theories,
likelihoods and samplers assembled from a YAML input, run with
`cosmofit run`, written as getdist chains -- see
[Upgrading from 1.x](#upgrading-from-1x). It is on
[PyPI](https://pypi.org/project/cosmofit/) with an
[API reference](https://salihyesil59.github.io/CosmoFit/). The full
release history -- what each version changed, and for many of them
how the bug was found rather than only that it was fixed -- lives in
**[CHANGELOG.md](CHANGELOG.md)**.

Where it stands today:

| | |
|---|---|
| datasets | **21**, from cosmic chronometers to the from-scratch CMB, each a native likelihood |
| models | **20** written out by hand, plus three routes to one that is not here |
| samplers | adaptive Metropolis-Hastings, emcee, dynesty, minimize, profile, Fisher, post-processing |
| tests | **1360**, on Python 3.11-3.13 on Linux, Windows and macOS, with and without every optional extra |
| notebooks | **17**, in five sections under [`examples/`](examples/), and YAML inputs beside them |

---

## Upgrading from 1.x

Most 1.x code keeps working, with warnings that say what to change.

* **The import name is `cosmofit`**, the same as on PyPI. `import CosmoFit`
  still works for this release, with a `DeprecationWarning`, and resolves
  to the same modules; pickles and saved chains that name `CosmoFit.*`
  load.
* **`Fitter` stays; its samplers move to the core.** `run_mcmc`,
  `run_nested`, `profile` and `fisher` warn and keep working. Their
  replacements take the same fitter: `cosmofit.compat.sample_on_core(fit,
  "mcmc", ...)`, `evidence_on_core`, `profile_on_core`, `fisher_on_core`.
  Everything after them -- `summary()`, `convergence()`, the derived
  posteriors, every `fit.plots.*` figure -- reads the result as before.
* **Saved chains** are now the core's getdist output under an `output`
  prefix. Chains 1.x saved as HDF5 still open with `Fitter.from_chain`;
  the GUI does not reuse them and samples those fits again once.
* **A `Fitter` becomes a 2.0 input** with `fit.to_info()` -- write it to
  YAML, add a `sampler` block, and `cosmofit run` takes over.
  `to_info(exact=True)` runs the fitter's own classes and gives its
  `chi2` exactly; the default runs the native theories, which put
  radiation into the late-time expansion and so move `chi2` by up to
  about one. Four parameters are spelled differently there: GCG's `A_s`
  and `alpha` are `A_gcg` and `alpha_gcg`, the power-law `f(T)`'s `n` is
  `n_ft`, Hu-Sawicki's `n` is `n_hs`.
* **Some numbers change because they were wrong.** Before the rewrite,
  every model and likelihood was audited, and the fixes move results:
  Pantheon+ is the published 1590-supernova sample (and its SH0ES
  Cepheids now calibrate), DES-SN5YR is the Dovekie recalibration and
  cited as one, the CC covariance keeps its systematic correlations,
  the default `tau` prior no longer double-counts `plik_lite`, DESI
  defaults to DR2, `f(R,T)`'s `H0` is the Hubble rate today, and a
  computed `r_d` sees each model's own early universe. Each is a
  CHANGELOG entry saying how it was found.
* **CAMB**: the CMB likelihoods need a CAMB with cmbant/CAMB#211, which
  fixed lensed spectra that turned NaN once a figure had been drawn. It
  is merged but not yet released -- install CAMB from its `main` branch
  until the next release.

---

## Roadmap

What each release meant, and what it took to get there, is its
entry in [CHANGELOG.md](CHANGELOG.md) rather than a section here --
a roadmap should say where the road goes next.

Most of what this section used to list has been built. What is left:

### Datasets

* **Strong-lensing time delays** (TDCOSMO / H0LiCOW) -- the one genuinely
  independent probe of `H0` still missing, and the third leg of the H0
  tension alongside the CMB and the distance ladder. Its per-lens
  likelihoods are skewed log-normals, which the tabulated machinery
  already here would handle.
* **DESI DR1/DR2 full-shape** -- deliberately *not* planned. DESI publish
  MCMC chains and the inputs to their modelling pipeline, but no
  compressed Gaussian summary, so using it means implementing an EFT
  model with its own nuisance parameters. That is outside a
  background-plus-linear-growth library, and it is not cheaply
  validatable. See the `2d59fe7` commit message.

### Physics

* **Massive neutrinos and `N_eff` in each model's own `E(z)`.** They
  reach the Boltzmann backend and the sound horizon and are treated
  exactly there, but the models' `E(z)` counts massive neutrinos inside
  `Omega_m` as pure matter -- right below z ~ 100, and increasingly wrong
  above it.

  Worth saying what this is *not*, because the three places that consume
  `E(z)` above z ~ 100 are each already handled. The sound horizon
  integrates the exact Fermi-Dirac density. The compressed Planck priors
  deliberately keep massive neutrinos inside `Omega_m` at every redshift,
  because a prediction has to share the compression's definitions rather
  than improve on them -- see `likelihoods/planck.py` on what mixing
  conventions costs. And the growth ODE solves a matter-plus-dark-energy
  equation, so starting it deep in *that* matter era is self-consistent.
  So this is narrower than it looks, and doing it carelessly would break
  agreements the library has validated.
* **`Sum m_nu` as a fitted parameter.** Everything needed is in place;
  what is missing is enough growth-suppression signal to constrain it.
  The compressed Planck priors cannot supply it -- they are blind to
  `m_nu` by construction, and `Fitter` says so.

### Performance

* **Batched log-posterior** (`emcee`'s `vectorize=True`): evaluating every
  walker in one call rather than one at a time. The one structural lever
  left -- worth perhaps 2-4x, and the only one that helps inside a
  Jupyter kernel, where `n_processes` does not. It needs a batch axis
  through every calculator, so it is a real refactor rather than a flag.

## Contributing

**[CONTRIBUTING.md](CONTRIBUTING.md)** -- how to get set up, why there
is no code formatter, and what a finished change looks like here.
The short version of the last one: a test that fails before the fix,
and a validation against something that does not share the machinery
being validated.

---

## References

CosmoFit bundles real observational data and implements published
cosmological models and statistical methods -- see
**[REFERENCES.md](REFERENCES.md)** for every dataset, model, and
methodology paper used, with links, and where each is used in the
code. If you use CosmoFit in a publication, cite the underlying
data/method papers relevant to what you used.

[CITATION.cff](CITATION.cff) describes the software itself, if you
also want to cite the library -- but the papers behind whatever you
actually used matter more.

---

## License

This project is licensed under the MIT License.
