Metadata-Version: 2.4
Name: gaplike
Version: 0.2.0
Summary: Gapped/windowed stationary Gaussian inference: gap-pattern generation, windowed frequency-domain covariances, and exact time-domain likelihoods (closed-form or matrix-free)
Author: Federico Pozzoli, Ollie Burke
License: MIT
Project-URL: Repository, https://github.com/gaplike/gaplike
Keywords: LISA,gravitational waves,data gaps,window functions,Whittle likelihood,Toeplitz,conjugate gradients
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Astronomy
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: ipykernel>=6.31.0
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.8
Provides-Extra: pe
Requires-Dist: emcee>=3; extra == "pe"
Requires-Dist: corner; extra == "pe"
Requires-Dist: matplotlib; extra == "pe"
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Provides-Extra: docs
Requires-Dist: sphinx>=7; extra == "docs"
Requires-Dist: furo; extra == "docs"
Requires-Dist: myst-parser>=2; extra == "docs"
Requires-Dist: sphinx-copybutton; extra == "docs"
Provides-Extra: mbhb
Requires-Dist: lisabeta; extra == "mbhb"
Dynamic: license-file

<picture>
  <source media="(prefers-color-scheme: dark)" srcset="docs/_static/logo-wordmark-dark.svg">
  <img alt="gaplike" src="docs/_static/logo-wordmark.svg" width="330">
</picture>

[![tests](https://github.com/gaplike/gaplike/actions/workflows/tests.yml/badge.svg)](https://github.com/gaplike/gaplike/actions/workflows/tests.yml)
[![docs](https://github.com/gaplike/gaplike/actions/workflows/docs.yml/badge.svg)](https://github.com/gaplike/gaplike/actions/workflows/docs.yml)
[![PyPI](https://img.shields.io/pypi/v/gaplike.svg)](https://pypi.org/project/gaplike/)
[![license: MIT](https://img.shields.io/badge/license-MIT-blue.svg)](LICENSE)

**Inference on gapped / windowed stationary Gaussian data.**

`gaplike` provides everything needed to do joint signal + noise-parameter
estimation when a stationary time series is interrupted by data gaps (or
multiplied by any window): gap-pattern generation with arbitrary patterns,
windowed frequency-domain covariances, and a hierarchy of likelihoods — from
the cheap Whittle approximation to the dense windowed covariance and the
**exact time-domain likelihood**, evaluated either in closed form (two
components, simultaneous diagonalization) or **matrix-free at any scale by
preconditioned conjugate gradients** (`gaplike.cg`: the covariance is never
formed, one solve is a few hundred FFT pairs).

<img src="assets/gap_leakage_web.gif" width="600"/>

*Four acts on the same record: a raw cut, then a taper on the record edges,
then sharp-edged gaps, then tapered gap edges. Watch the spectral kernel, the
modelled PSD and the bin-to-bin correlation matrix move together. Rendered by
[`notebooks/anim_leakage.py`](notebooks/anim_leakage.py).*

Companion package to *"Zurückbleiben bitte: the impact of gaps on
noise and signal parameter inference"* (O. Burke, F. Pozzoli & M. Muratore); the `paper/`
folder reproduces every figure of that paper.

Only `numpy` and `scipy` are required. The two-component LISA TDI-2 A/E noise
model is built in; any user-defined PSD components and any waveform work the
same way.

**Documentation: <https://gaplike.github.io/gaplike/>**

## Install

`gaplike` needs only `numpy` and `scipy`. Pick **one** environment and install
everything into it: the commonest way to lose an afternoon here is to install
with one tool and run with another.

### From PyPI

```bash
pip install gaplike                 # everything except paper/ and the notebooks
pip install "gaplike[pe]"           # + emcee, corner, matplotlib
```

That is all the package itself needs. Everything below installs from a clone
instead, which you want only if you also need the paper reproduction pipeline
in `paper/`, the notebooks, or the test suite.

### Into an environment you already have (conda, system Python, ...)

```bash
python -m pip install -e ".[pe,dev,docs]"
python -m sphinx -b html docs docs/_build/html
python -m pytest
```

The `python -m` form ties the installer and the command to the same
interpreter, which is what keeps them in step.

### Into a fresh virtual environment

```bash
uv venv
source .venv/bin/activate          # <- do not skip this
uv pip install -e ".[pe,dev,docs]"
sphinx-build -b html docs docs/_build/html
pytest
```

Without the `activate` line, `uv` installs into `.venv` while `sphinx-build`
and `pytest` still come from whatever is on your `PATH`, and you get a
confusing "module not found" for something you just watched install.
`uv run <command>` does the same job without activating.

### Extras

| extra | pulls in | needed for |
|---|---|---|
| `pe` | `emcee`, `corner`, `matplotlib` | sampling, the notebooks, every paper figure |
| `dev` | `pytest` | the test suite |
| `docs` | `sphinx`, `furo`, `myst-parser`, `sphinx-copybutton` | building the documentation |
| `mbhb` | `lisabeta` | the MBHB waveform adapter only |

The MBHB waveform adapter additionally needs
[lisabeta](https://gitlab.in2p3.fr/marsat/lisabeta), which is **not** pulled
in by any of the lines above:

```bash
uv pip install -e ".[pe,mbhb]"          # + lisabeta (from PyPI)
```

`[mbhb]` is a separate extra because exactly one factory needs it. Everything
else — every likelihood, every gap builder, the CG solver, the paper figures
from the cached chains, and `notebooks/exact_inference_demo.ipynb` — runs
without lisabeta, and the import happens inside
`gaplike.waveform.lisabeta_mbhb_ae`, so a missing install fails only there.

Installing from PyPI gives wheels. Building the GitLab sources instead
requires FFTW and a CMake toolchain, and the FFTW location has to be handed
to the build explicitly:

```bash
brew install fftw                       # macOS;  apt install libfftw3-dev on Debian/Ubuntu
SKBUILD_CMAKE_DEFINE="FFTW_ROOT=$(brew --prefix fftw)" \
  uv pip install "lisabeta @ git+https://gitlab.in2p3.fr/marsat/lisabeta.git"
```

## Quickstart

```python
import numpy as np
import gaplike as gl

n, dt = 2880, 15.0                          # 12 h at 15 s cadence
f_lo, f_hi = 1e-4, 3.1e-2                   # analysis band [Hz]

# --- 1. gaps: ANY pattern -------------------------------------------------
mask = gl.gaps.periodic_mask(n, 10, 50)              # comb: 150 s out of 750 s
# mask = gl.gaps.random_mask(n, dt, rate_per_day=40, duration_s=600, rng=1)
# mask = gl.gaps.mask_from_intervals(n, dt, [(21000, 25000)])
gate = gl.gaps.gate_from_mask(mask, dt, taper_s=0.0) # rectangular (or tapered)
w    = gl.gaps.effective_window(gate, gl.gaps.segment_window(n, 0.05))

# --- 2. noise model: two components; lam = log10 deviations of the component
#        POWERS, so Sigma(lam) = 10^(lam_0) C_0 + 10^(lam_1) C_1, truth at 0
comps = list(gl.psd.lisa_tdi2_ae().values())         # [S_tm(f), S_oms(f)]
S_tot = gl.psd.one_sided_grid(lambda f: sum(c(f) for c in comps), n, dt)
x = gl.simulate.noise_td(S_tot, dt, np.random.default_rng(0), nch=2)

# --- 3. likelihood hierarchy ------------------------------------------------
L_full = gl.FullCovariance.from_window(w, comps, dt, f_lo, f_hi)   # dense windowed cov
L_conv = gl.DiagonalLikelihood.convolved(w, comps, dt, f_lo, f_hi) # exact diagonal
L_whit = gl.DiagonalLikelihood.whittle(w, comps, dt, f_lo, f_hi,
                                       scale_by_window_power=True) # ~ Whittle
L_td   = gl.TimeDomainExact(mask, comps, dt)                       # EXACT, no window

for L in (L_td, L_full, L_conv, L_whit):
    r = L.transform(x)                       # residual = data - template
    print(type(L).__name__, L.loglike(r, (0.0, 0.0)) - L.loglike(r, (0.05, 0.0)))

# --- 4. the same exact likelihood, matrix-free (large N) --------------------
rcg  = gl.RestrictedCG(mask, comps, dt, rtol=1e-8)   # no O(m^2) storage, no O(m^3) setup
quad, iters = rcg.quad_form((0.0, 0.0), L_td.transform(x))
# pair with a determinant: closed form from TimeDomainExact (2 components),
# or stochastic Lanczos quadrature in general
```

Each tier carries its own constant offset from its internal normalization, so
log-likelihoods are **not** comparable across tiers — only *differences*
(between parameter points, or between templates) are meaningful.

Templates from any waveform: wrap a frequency-domain callable in
`gl.Waveform(fd_func, n, dt)` and use `.td(theta)`; then
`resid = L.transform(data_td - h_td)`. For LISA MBHBs,
`gl.waveform.lisabeta_mbhb_ae(n, dt, f_lo, f_hi)` gives the (A, E) TDI-2
IMRPhenomHM waveform used in the paper — with any source parameters, any
parametrization (`to_params`) and any approximant (`wf_kw`); see
`notebooks/mbhb_gaps_demo.ipynb`.

## Reproducing the paper figures

All commands run from `paper/` with the `[pe]` extra installed.  The cached
PE chains ship in `paper/results/`, so every chain-based figure regenerates
in minutes without lisabeta and without sampling; lisabeta and the from-
scratch runs are needed only to regenerate the chains themselves.

```bash
cd paper

# -- chain-based figures, from the cached results/ (no lisabeta needed) -----
python make_figures.py        # corner_key_{A,B}, corner_noise_{A,B},
                              # corner_fullkey_{A,B}, upsilon_xi(.png/_heatmap),
                              # overview                           (~10 min)
python plot_scenC.py          # corner_key_C_full_diag_td, corner_noise_C_...,
                              # cov_colormap_C
python plot_ABC.py            # corner_key_full_ABC
python fig_cov_colormap.py    # cov_colormap_ABC

# -- the CG scaling figure (gaplike only: no lisabeta, no chains) -----------
python fig_cg_scaling.py      # times a^T Sigma_OO^-1 a, dense vs matrix-free
                              # preconditioned CG -> results/cg_scaling.json
                              # (~15-30 min; --kmax-cg/--kmax-dense to shorten)
python _mkfig.py              # -> figures/cg_scaling.{png,pdf}

# -- full regeneration of the chains (lisabeta + hours of sampling) ---------
python driver.py              # A/B x {full,diag,bare,psd}   (~80 min, 2 cores)
python scenC_run.py           # C: FD full + convolved diagonal + exact TD PE
python ensemble_noise.py      # 300-realization ensemble (Table: pseudo-true
                              # points, sandwich scatter)
```

Figure-to-script map: every `corner_*`, `upsilon_xi*` and `overview` come
from `make_figures.py`; the three standalone scripts above cover the
scenario-C corner, the A/B/C overlay corner and the covariance colormaps;
`fig_cg_scaling.py` + `_mkfig.py` produce the dense-vs-CG scaling figure.
`results/cg_scaling.json` ships with the reference timings (2 cores of an
Intel Xeon at 2.80 GHz) so `_mkfig.py` works out of the box.

## The likelihood hierarchy

| class | covariance model | cost / eval | guarantees |
|---|---|---|---|
| `TimeDomainExact` | stationary covariance restricted to observed samples | O(m) after one O(m³) setup | **exact**: scatter/width = width/true-width = 1 by construction; no window anywhere |
| `FullCovariance` | dense windowed FD covariance (band-restricted) | O(r) after one O(N_b³) setup | exact within the band-restricted, circularly-symmetric FD reduction |
| `DiagonalLikelihood.convolved` | exact diagonal of the windowed covariance | O(N_b) | unbiased widths *per bin*; ignores bin–bin correlations |
| `DiagonalLikelihood.whittle` | raw PSD (optionally × window power W₂) | O(N_b) | Whittle / "normalizing constant" approximation |
| `gaplike.cg` / `RestrictedCG` | same model as `TimeDomainExact`, matrix-free | 4 FFTs × a few hundred iterations, **no setup, no storage** | exact quadratic forms to a chosen tolerance; any number of components; determinant not included |

Noise parameters are the same throughout: `lam_k` are log10 deviations of the
component **powers** from their reference values, `Sigma(lam) = sum_k
10^(lam_k) C_k`, with the truth at `lam = 0`.

Both pencil classes (`TimeDomainExact`, `FullCovariance`) additionally assume
the covariance is linear in **exactly two** components: one simultaneous
diagonalization then makes every likelihood evaluation, Fisher matrix and
determinant closed-form. The diagonal tiers and the conjugate gradient route
accept any number of components; `gaplike.cg` is the escape hatch when the
two-component structure is not available (spline-knot spectra, extra
components) or when m is too large to factorize anything.

The paper's analytic mis-specification machinery (Fisher blocks under the
true windowed covariance, leading-order biases, MLE scatter, the Υ/Ξ
diagnostics, KL pseudo-true parameters, Godambe–White sandwich) lives in
`paper/diagnostics.py`, on top of the package.

## Honest error bars: Υ and Ξ

Two numbers summarize what an approximate likelihood does to an inference:
**Υ = true scatter / quoted width** (calibration — is the error bar honest?)
and **Ξ = quoted width / exact-analysis width** (efficiency — was all the
information used?). Honest data cannot beat the exact analysis, so
**Υ·Ξ ≥ 1**, with equality only for the exact likelihood.

<img src="assets/upsilon_xi_quadrants.gif" width="600"/>

*The four regimes as repeated experiments — every case a measured value
from Table III of the paper. A PP plot sees only Υ: the
calibrated-but-inflated analysis (top right, Ξ = 9.3) passes any PP test
while quoting error bars 9.3× wider than the data allow; the companion
`assets/upsilon_xi_pp.gif` shows exactly that. Rendered by
[`notebooks/anim_upsilon_xi.py`](notebooks/anim_upsilon_xi.py).*

[**The Υ and Ξ diagnostics**](https://gaplike.github.io/gaplike/upsilon_xi.html)
page in the docs derives both quantities, compares them with PP plots, and
explains how to read them. The
[**interactive Υ–Ξ explorer**](https://gaplike.github.io/gaplike/_static/upsilon_xi_explorer.html)
(also a single offline file, [`assets/upsilon_xi_explorer.html`](assets/upsilon_xi_explorer.html))
puts them on sliders: drag a point around the (Ξ, Υ) plane with the
forbidden Υ·Ξ < 1 region hatched out, click the paper's measured anchors to
load them, and watch the intervals and the PP curve respond.

## Repository layout

| path | content |
|---|---|
| `src/gaplike/` | the package: `gaps`, `psd`, `covariance`, `simulate`, `likelihood`, `cg`, `waveform` |
| `tests/` | unit tests (small-N brute-force exactness of every tier, CG vs dense) + machine-precision regression against the paper pipeline |
| `notebooks/exact_inference_demo.ipynb` | executable end-to-end example (**no lisabeta**): gapped two-channel LISA noise + toy chirp, joint 4-parameter PE with `TimeDomainExact` and `FullCovariance`, comparison corner |
| `notebooks/mbhb_gaps_demo.ipynb` | inject an MBHB with arbitrary parameters (lisabeta), compare five gap configurations, and set up a full 13-parameter inference ready to launch |
| `paper/` | full paper reproduction (scenarios A/B/C, PE, every figure incl. the CG scaling benchmark) — see `paper/README.md` |
| `assets/`, `notebooks/anim_upsilon_xi.py` | talk media: the gap-leakage animation, the Υ/Ξ animations and the interactive Υ–Ξ explorer |

## Tests

```bash
uv pip install -e ".[dev]"
pytest                    # ~40 s; the paper regression needs tests/data/reference.json
```

`tests/test_paper_regression.py` rebuilds the paper's three gap scenarios
purely from package primitives and reproduces the original pipeline's
windows, covariances, noise realization and all likelihood values to ~1e-9.

## Citing

If you use `gaplike`, please cite the software (see `CITATION.cff`) together
with Burke, Pozzoli & Muratore (2026), *Zurückbleiben bitte: the impact of window
functions on noise and signal parameter inference* (in prep.).

## License

MIT — see `LICENSE`.
