From 67e4b8da5f7518c2984d8e43bc2240c07d2cdade Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 30 Jun 2026 14:34:02 +0200 Subject: [PATCH 01/11] docs: spec for strain PPD plotting Co-Authored-By: Claude Opus 4.8 (1M context) --- .../2026-06-30-strain-ppd-plotting-design.md | 89 +++++++++++++++++++ 1 file changed, 89 insertions(+) create mode 100644 docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md new file mode 100644 index 000000000..277371c60 --- /dev/null +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -0,0 +1,89 @@ +# Strain PPD plotting — design + +**Date:** 2026-06-30 +**Branch:** `hackathon/ppd-plotting` + +## Goal + +Given a dingo `Result` HDF5 file, compute and plot the posterior-predictive +distribution (PPD) of the **whitened strain** in both the time domain and the +frequency domain, per detector, and save the figures as PNGs. + +Derived from `/work/nihargupte/analysis/dingo/asimov/eccentricity/plotting/ +production_plots/strain_posterior.ipynb`, but stripped of the +eccentricity/glitch/prior-hull extras and generalized to a clean, reusable +script. + +## Deliverable + +`dingo/gw/inference/ppd.py` — importable functions plus an `argparse` CLI: + +``` +python -m dingo.gw.inference.ppd \ + --outdir --num-waveforms 1000 --num-processes 32 \ + --central-time 8 --zoom -0.7 0.1 +``` + +Demo run target: GW230709_122727 EAS importance-sampling result +(`prod_o4a/working/S230709bi/Exp19_more_points_1/result/ +Exp19_more_points_1_data0_1372940865-2_importance_sampling.hdf5`), producing +`_ppd_td.png` and `_ppd_fd.png`. + +## Pipeline (per detector) + +1. `result = Result(file_name=...)`; `result._build_likelihood()`. + `domain = result.domain.base_domain if hasattr(result.domain, "base_domain") + else result.domain`. +2. **Sample selection.** Resample `N` parameter rows from `result.samples` with + replacement, probability ∝ importance `weights`, using a seeded + `np.random.default_rng`. This is a proper posterior-predictive draw and + respects the importance weights. (The notebook instead kept the top-90% + credible set then truncated to `N`; we deliberately diverge.) +3. Generate whitened FD waveforms: + `apply_func_with_multiprocessing(result.likelihood.signal, samples_df, + num_processes)`. Each `wf["waveform"][ifo]` is **already whitened** + (`/asd/noise_std`), the same convention applied to the data, so model and + data share one y-scale. + +## Time domain (`_ppd_td.png`, one row per detector) + +- Each draw: phase-shift to place merger at `CENTRAL_TIME` + (`wf *= exp(2πi·f·CENTRAL_TIME)`), then `one_sided_fd_to_td`, keep real part. +- Bands across draws via `np.nanpercentile`: **median line + shaded 50% + (25–75) + 90% (5–95)**. (Notebook used a min/max envelope; we deliberately + diverge to percentile credible bands.) +- Overlay: whitened data trace (`√(4Δf)·context_waveform/asd`, same shift + + IFFT, light sliding-window average), grey. +- x-axis: time to merger (s), linear, zoomed via `--zoom` (default + `(-0.7, 0.1)` for GW230709). y-axis: whitened strain (dimensionless, σ units). + +## Frequency domain (`_ppd_fd.png`, one row per detector) + +- Per draw: `|whitened h(f)|` (magnitude of the complex whitened FD waveform). +- Bands: median + 50%/90% percentiles vs frequency. +- Overlay: `|whitened data|` = `|√(4Δf)·context_waveform/asd|` — the noisy + ~O(1) floor the signal rises above. +- x-axis: frequency [Hz], **log scale**. y-axis: dimensionless whitened + amplitude. + +## Reused logic + +`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook: it builds +the full Hermitian two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` +(normalization tied to the whitening so noise has unit variance), +`dt = 1/(2·f_max)`, output length `2n-1`. Correct and whitening-consistent. + +## Testing / determinism + +- Seeded `np.random.default_rng` for resampling → deterministic figures. +- `tests/gw/inference/test_ppd.py`: assert `one_sided_fd_to_td` round-trips — + a one-sided spectrum with a single populated frequency bin returns a real + cosine of the expected period (exercises the Hermitian-mirror + normalization + logic, the one piece of non-trivial new math). + +## Boundaries / conventions + +- Lives in `gw/`, imports `gw/` + `core/` (allowed direction). No new classes, + no new dependencies. +- Whitening conventions for data (`√(4Δf)/asd`) and model (likelihood's + `/asd/noise_std`, `noise_std = 1/√(4Δf)`) are identical and must be preserved. From b142472b006b232cde45b0ec7ccf9ebfd62107a9 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 30 Jun 2026 14:43:55 +0200 Subject: [PATCH 02/11] docs: PPD plots as Result methods, not standalone module Co-Authored-By: Claude Opus 4.8 (1M context) --- .../2026-06-30-strain-ppd-plotting-design.md | 42 ++++++++++++------- 1 file changed, 28 insertions(+), 14 deletions(-) diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index 277371c60..a0188fa0f 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -16,18 +16,31 @@ script. ## Deliverable -`dingo/gw/inference/ppd.py` — importable functions plus an `argparse` CLI: - -``` -python -m dingo.gw.inference.ppd \ - --outdir --num-waveforms 1000 --num-processes 32 \ - --central-time 8 --zoom -0.7 0.1 +Two plotting **methods on the GW `Result` class** (`dingo/gw/result.py`), +mirroring the existing `result.plot_corner` API: + +```python +result.plot_ppd_td(filename="ppd_td.png", num_waveforms=1000, + num_processes=1, central_time=8.0, zoom=None, ppd=None) +result.plot_ppd_fd(filename="ppd_fd.png", num_waveforms=1000, + num_processes=1, ppd=None) ``` +They live on the GW `Result` (subclass of `CoreResult`), **not** `core/`, +because they need the GW likelihood, detector projection, and ASD whitening — +`core/` stays domain-agnostic. + +Shared expensive computation lives in a private +`_compute_ppd(self, num_waveforms, num_processes, seed)` (used by both methods, +so extraction is justified) returning `{domain, ifos, wf_fd, data_fd}`. Each +public method accepts an optional precomputed `ppd=` dict so a caller can +generate the draws once and render both plots without regenerating waveforms. + Demo run target: GW230709_122727 EAS importance-sampling result (`prod_o4a/working/S230709bi/Exp19_more_points_1/result/ -Exp19_more_points_1_data0_1372940865-2_importance_sampling.hdf5`), producing -`_ppd_td.png` and `_ppd_fd.png`. +Exp19_more_points_1_data0_1372940865-2_importance_sampling.hdf5`). A few-line +script loads the `Result` and calls both methods, producing +`GW230709_ppd_td.png` and `GW230709_ppd_fd.png`. ## Pipeline (per detector) @@ -68,10 +81,11 @@ Exp19_more_points_1_data0_1372940865-2_importance_sampling.hdf5`), producing ## Reused logic -`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook: it builds -the full Hermitian two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` -(normalization tied to the whitening so noise has unit variance), -`dt = 1/(2·f_max)`, output length `2n-1`. Correct and whitening-consistent. +`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook (added as +a module-level function in `dingo/gw/result.py`): it builds the full Hermitian +two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` (normalization tied +to the whitening so noise has unit variance), `dt = 1/(2·f_max)`, output length +`2n-1`. Correct and whitening-consistent. ## Testing / determinism @@ -83,7 +97,7 @@ the full Hermitian two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` ## Boundaries / conventions -- Lives in `gw/`, imports `gw/` + `core/` (allowed direction). No new classes, - no new dependencies. +- Methods on the GW `Result` in `gw/result.py`; imports `gw/` + `core/` + (allowed direction). No new classes, no new dependencies. - Whitening conventions for data (`√(4Δf)/asd`) and model (likelihood's `/asd/noise_std`, `noise_std = 1/√(4Δf)`) are identical and must be preserved. From 3b381f13e87a5bd7572792822dd0e200b60ac625 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 30 Jun 2026 15:16:25 +0200 Subject: [PATCH 03/11] docs: min/max envelope + hidden time offset for PPD plots Co-Authored-By: Claude Opus 4.8 (1M context) --- .../2026-06-30-strain-ppd-plotting-design.md | 53 ++++++++++++------- 1 file changed, 33 insertions(+), 20 deletions(-) diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index a0188fa0f..d6880bea5 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -20,10 +20,10 @@ Two plotting **methods on the GW `Result` class** (`dingo/gw/result.py`), mirroring the existing `result.plot_corner` API: ```python -result.plot_ppd_td(filename="ppd_td.png", num_waveforms=1000, - num_processes=1, central_time=8.0, zoom=None, ppd=None) -result.plot_ppd_fd(filename="ppd_fd.png", num_waveforms=1000, - num_processes=1, ppd=None) +result.plot_ppd_td(filename="ppd_td.png", credible_interval=0.9, + num_waveforms=1000, num_processes=1, zoom=None, ppd=None) +result.plot_ppd_fd(filename="ppd_fd.png", credible_interval=0.9, + num_waveforms=1000, num_processes=1, ppd=None) ``` They live on the GW `Result` (subclass of `CoreResult`), **not** `core/`, @@ -31,10 +31,11 @@ because they need the GW likelihood, detector projection, and ASD whitening — `core/` stays domain-agnostic. Shared expensive computation lives in a private -`_compute_ppd(self, num_waveforms, num_processes, seed)` (used by both methods, -so extraction is justified) returning `{domain, ifos, wf_fd, data_fd}`. Each -public method accepts an optional precomputed `ppd=` dict so a caller can -generate the draws once and render both plots without regenerating waveforms. +`_compute_ppd(self, credible_interval, num_waveforms, num_processes, seed)` +(used by both methods, so extraction is justified) returning +`{domain, ifos, wf_fd, data_fd}`. Each public method accepts an optional +precomputed `ppd=` dict so a caller can generate the draws once and render both +plots without regenerating waveforms. Demo run target: GW230709_122727 EAS importance-sampling result (`prod_o4a/working/S230709bi/Exp19_more_points_1/result/ @@ -47,33 +48,45 @@ script loads the `Result` and calls both methods, producing 1. `result = Result(file_name=...)`; `result._build_likelihood()`. `domain = result.domain.base_domain if hasattr(result.domain, "base_domain") else result.domain`. -2. **Sample selection.** Resample `N` parameter rows from `result.samples` with - replacement, probability ∝ importance `weights`, using a seeded - `np.random.default_rng`. This is a proper posterior-predictive draw and - respects the importance weights. (The notebook instead kept the top-90% - credible set then truncated to `N`; we deliberately diverge.) +2. **Sample selection.** Build the `credible_interval` (default 0.9) credible + set: normalize the importance/posterior weights, sort descending, and keep + the highest-weight samples whose cumulative weight ≤ `credible_interval` + (the X% credible region — same definition as the notebook's + `select_from_n_percent_interval`). Then draw up to `num_waveforms` rows + **uniformly at random** from that set (seeded `np.random.default_rng`) so the + envelope represents the whole interval rather than only the peak (the + notebook took the top-N highest-weight rows, biasing toward the peak — we + diverge here). 3. Generate whitened FD waveforms: `apply_func_with_multiprocessing(result.likelihood.signal, samples_df, num_processes)`. Each `wf["waveform"][ifo]` is **already whitened** (`/asd/noise_std`), the same convention applied to the data, so model and data share one y-scale. +### Time offset (hidden from the user) + +The merger always lands at **t = 0**. Internally each waveform is phase-shifted +by a fixed offset derived from the data segment (`t0 = 1 / domain.delta_f`, the +segment duration `T`, which places the dingo coalescence at the segment edge), +via `wf *= exp(2πi·f·t0)`, then the time axis has `t0` subtracted. `t0` is a +private detail — not a method argument. The implementation verifies the demo +merger sits at 0 and adjusts the offset derivation if not. + ## Time domain (`_ppd_td.png`, one row per detector) -- Each draw: phase-shift to place merger at `CENTRAL_TIME` - (`wf *= exp(2πi·f·CENTRAL_TIME)`), then `one_sided_fd_to_td`, keep real part. -- Bands across draws via `np.nanpercentile`: **median line + shaded 50% - (25–75) + 90% (5–95)**. (Notebook used a min/max envelope; we deliberately - diverge to percentile credible bands.) +- Each draw: phase-shift by `t0`, `one_sided_fd_to_td`, keep real part. +- Band = pointwise **min/max envelope** across the selected draws + (`np.nanmin` / `np.nanmax`, ignoring NaN waveforms), drawn with + `fill_between`. - Overlay: whitened data trace (`√(4Δf)·context_waveform/asd`, same shift + IFFT, light sliding-window average), grey. -- x-axis: time to merger (s), linear, zoomed via `--zoom` (default +- x-axis: time to merger (s), linear, zoomed via `zoom=` (default `(-0.7, 0.1)` for GW230709). y-axis: whitened strain (dimensionless, σ units). ## Frequency domain (`_ppd_fd.png`, one row per detector) - Per draw: `|whitened h(f)|` (magnitude of the complex whitened FD waveform). -- Bands: median + 50%/90% percentiles vs frequency. +- Band = min/max envelope of `|whitened h(f)|` across the selected draws. - Overlay: `|whitened data|` = `|√(4Δf)·context_waveform/asd|` — the noisy ~O(1) floor the signal rises above. - x-axis: frequency [Hz], **log scale**. y-axis: dimensionless whitened From 6ad94f96c42b2cb2d14d909c766e9e87f8fe040e Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 12:12:44 +0200 Subject: [PATCH 04/11] feat(gw): time-domain whitened-strain PPD plots on Result Add Result.plot_ppd_td / _compute_ppd to render the time-domain whitened-strain posterior-predictive distribution, overlaying the Dingo and (when importance- sampled) Dingo-IS credible-set envelopes on the whitened detector data, mirroring plot_corner. Waveform generation and the inverse FFT live in a new dingo/gw/utils/plotting.py. Restrict draws to the prior support (finite log_prior) before generating -- the same in-prior restriction importance_sample applies. The flow proposal leaks a small fraction of unphysical out-of-prior draws (e.g. tiny negative spin magnitudes) that are exactly the ones failing waveform generation; filtering them removes those failures (guaranteed for an IS result, whose in-prior samples all generated during IS) and keeps unphysical waveforms out of the envelope. Measured on the GW150914 fixture: 43/5000 draws out-of-prior, all generation failures among them, zero in-prior failures. With the filter in place there is no per-sample failure handling: signal() is called directly and the envelope uses plain min/max; a genuine in-prior model failure is left to raise. Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/result.py | 181 +++++++++++++++++- dingo/gw/utils/__init__.py | 1 + dingo/gw/utils/plotting.py | 143 ++++++++++++++ .../2026-06-30-strain-ppd-plotting-design.md | 155 +++++++++------ tests/gw/test_ppd.py | 32 ++++ 5 files changed, 455 insertions(+), 57 deletions(-) create mode 100644 dingo/gw/utils/__init__.py create mode 100644 dingo/gw/utils/plotting.py create mode 100644 tests/gw/test_ppd.py diff --git a/dingo/gw/result.py b/dingo/gw/result.py index 51616032d..eb7c25e88 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -1,6 +1,6 @@ import copy import time -from typing import Optional +from typing import Optional, Tuple import numpy as np import yaml @@ -19,6 +19,7 @@ from dingo.gw.domains import build_domain from dingo.gw.gwutils import get_extrinsic_prior_dict from dingo.gw.likelihood import StationaryGaussianGWLikelihood +from dingo.gw.utils.plotting import plot_ppd_td as _plot_ppd_td from dingo.gw.prior import build_prior_with_defaults from dingo.core.utils.backward_compatibility import check_minimum_version @@ -781,3 +782,181 @@ def pesummary_prior(self): except AttributeError: continue return prior + + def _compute_ppd( + self, + credible_interval: float = 0.9, + num_waveforms: int = 1000, + num_processes: int = 1, + seed: int = RANDOM_STATE, + ) -> Tuple[dict, dict]: + """Generate whitened waveforms for the time-domain strain PPD. + + Builds a ``credible_interval`` highest-posterior-density credible set for each + available posterior -- ``"dingo"`` (the raw network samples) and, when the result is + importance-sampled, ``"dingo-is"`` -- draws up to ``num_waveforms`` of each uniformly + at random (so the envelope spans the whole interval, not just the peak), and projects + them onto the detectors via the likelihood. Model waveforms come back already whitened + (``h / asd / noise_std``), matching the whitened data. + + Parameters + ---------- + credible_interval : float + Posterior probability of the credible set to draw from (default 0.9). + num_waveforms : int + Maximum number of waveforms drawn from each credible set. + num_processes : int + Processes used for waveform generation. + seed : int + Seed for the uniform subsampling, for deterministic figures. + + Returns + ------- + wf_fd : dict + ``{mode: {ifo: complex ndarray (n_kept, n_freq)}}`` whitened model waveforms, + where ``mode`` is ``"dingo"`` (always) and ``"dingo-is"`` (present iff the result + is importance-sampled). Rows are the successfully generated draws; the rare failed + draw is dropped, so ``n_kept`` may be slightly below the number requested. + data_fd : dict + ``{ifo: complex ndarray (n_freq,)}`` whitened detector data (shared across modes). + + Both are consumed by :func:`dingo.gw.utils.plotting.plot_ppd_td` (together with + ``self.domain``). + + Notes + ----- + Per-mode credible set. The ``credible_interval`` (= X) region is the highest- + posterior-density set ``R = {theta : p(theta) >= c}`` with ``\\int_R p = X``. It needs + *two distinct* quantities per sample -- a probability **mass** and a **ranking key** -- + not one weight. Given samples ``theta_i ~ g`` (the sampling density) and target density + ``p``, self-normalized importance sampling estimates any set's mass as + ``P(A) ~= sum_i m_i 1[theta_i in A]`` with ``m_i = (p_i/g_i) / sum_j (p_j/g_j)`` + + - ``"dingo"``: ``m_i propto 1/N``; rank by ``log_prob``. + - ``"dingo-is"``: ``m_i propto p/q``; rank by ``log_likelihood + log_prior``. + + Ranking by the density and accumulating ``m_i`` until the cumulative mass reaches X + keeps exactly the samples above the density threshold ``c`` + + Prior filter. Samples are first restricted to the prior support (finite ``log_prior``). + The flow proposal is a smooth density over a box and leaks a small fraction of draws + outside the prior (e.g. tiny negative spin magnitudes); these are unphysical, carry + zero target mass, and are the samples that fail waveform generation. This is exactly + the in-prior restriction that :meth:`importance_sample` applies, so for an importance- + sampled result -- whose in-prior samples were all generated successfully during IS -- + the filter is guaranteed to leave no generation failures. A genuine in-prior model + failure (e.g. an approximant interpolation edge on a not-yet-importance-sampled result) + is left to raise, surfacing the offending sample rather than being silently dropped. + """ + if getattr(self, "likelihood", None) is None: + self._build_likelihood() + + domain = ( + self.domain.base_domain + if hasattr(self.domain, "base_domain") + else self.domain + ) + ifos = list(self.context["waveform"].keys()) + rng = np.random.default_rng(seed) + + # Whitened data, same convention as the (already-whitened) model waveforms. + data_fd = { + ifo: np.sqrt(4 * domain.delta_f) + * self.context["waveform"][ifo] + / self.context["asds"][ifo] + for ifo in ifos + } + + # Restrict to the prior support before generating anything. The flow proposal is a + # smooth density over a box and leaks a small fraction of draws outside the prior + # (e.g. tiny negative spin magnitudes); these are unphysical and are exactly the + # draws that fail waveform generation. Dropping them here -- the same in-prior + # restriction importance sampling uses -- removes those failures (guaranteed for an + # importance-sampled result, whose in-prior samples all generated during IS) and + # keeps unphysical waveforms out of the envelope. + non_fixed = [ + k + for k, v in self.prior.items() + if not isinstance(v, (Constraint, DeltaFunction)) + ] + in_prior = np.isfinite( + np.asarray(self.prior.ln_prob(self.samples[non_fixed], axis=0), dtype=float) + ) + + # One credible set per available posterior: "dingo" always; "dingo-is" when the + # result is importance-sampled (weights present). + wf_fd = {} + for mode in ("dingo", "dingo-is"): + if mode == "dingo": + mass = in_prior.astype(float) + rank = self.samples["log_prob"].to_numpy() + else: + if "weights" not in self.samples: + continue + mass = self.samples["weights"].to_numpy().astype(float) + rank = ( + self.samples["log_likelihood"] + self.samples["log_prior"] + ).to_numpy() + # Out-of-prior draws get zero mass and bottom rank, so they never enter the + # credible set (this also sidesteps NaN log_likelihood on unevaluated samples). + mass = np.where(in_prior, mass, 0.0) + mass = mass / mass.sum() + rank = np.where(in_prior, rank, -np.inf) + order = np.argsort(rank)[::-1] + n_keep = int(np.searchsorted(np.cumsum(mass[order]), credible_interval)) + 1 + credible_idx = order[:n_keep] + # Subsample uniformly from the credible set so the min/max envelope spans it. + if len(credible_idx) > num_waveforms: + credible_idx = rng.choice( + credible_idx, size=num_waveforms, replace=False + ) + theta = self.samples.iloc[credible_idx] + + # In-prior samples generate cleanly (guaranteed for an importance-sampled + # result, whose in-prior draws were all generated during IS); a genuine model + # failure here should surface, not be silently dropped. + signals = apply_func_with_multiprocessing( + self.likelihood.signal, + theta, + num_processes=num_processes, + ) + wf_fd[mode] = { + ifo: np.array([s["waveform"][ifo] for s in signals]) for ifo in ifos + } + + if wf_fd["dingo"][ifos[0]].shape[1] != len(domain()): + raise ValueError( + "Whitened waveform length does not match the inverse-FFT domain; " + "time-domain PPD requires a uniform frequency domain." + ) + return wf_fd, data_fd + + def plot_ppd_td( + self, + filename: str = "ppd_td.png", + credible_interval: float = 0.9, + num_waveforms: int = 1000, + num_processes: int = 1, + zoom: Optional[Tuple[float, float]] = None, + ) -> Tuple[dict, dict]: + """Plot the time-domain whitened-strain posterior-predictive distribution. + + For each detector, inverse-FFTs whitened waveforms to the time domain with the merger + at t = 0 and shades the pointwise min/max envelope; the whitened detector data is + overlaid in grey. Mirroring :meth:`plot_corner`, both the **Dingo** posterior and (when + the result is importance-sampled) the **Dingo-IS** posterior are overlaid on one + figure. + + ``zoom`` is the (left, right) x-limit in seconds-to-merger; defaults to (-1.0, 0.2). + Returns the ``(wf_fd, data_fd)`` tuple that was plotted. + """ + wf_fd, data_fd = self._compute_ppd( + credible_interval, num_waveforms, num_processes + ) + domain = ( + self.domain.base_domain + if hasattr(self.domain, "base_domain") + else self.domain + ) + _plot_ppd_td(wf_fd, data_fd, domain, filename=filename, zoom=zoom) + return wf_fd, data_fd diff --git a/dingo/gw/utils/__init__.py b/dingo/gw/utils/__init__.py new file mode 100644 index 000000000..630ed3758 --- /dev/null +++ b/dingo/gw/utils/__init__.py @@ -0,0 +1 @@ +"""GW-specific utilities.""" diff --git a/dingo/gw/utils/plotting.py b/dingo/gw/utils/plotting.py new file mode 100644 index 000000000..ba6d1b0af --- /dev/null +++ b/dingo/gw/utils/plotting.py @@ -0,0 +1,143 @@ +"""Time-domain strain posterior-predictive-distribution (PPD) plotting for GW results. + +GW-specific counterpart to :mod:`dingo.core.utils.plotting`: renders the whitened-strain +PPD envelopes produced by :meth:`dingo.gw.result.Result._compute_ppd` in the time domain. +:func:`plot_ppd_td` overlays one envelope per posterior mode in ``wf_fd`` (``"dingo"`` and, +when available, ``"dingo-is"``), mirroring how ``result.plot_corner`` shows both. + +Inputs come straight from ``Result._compute_ppd``: ``wf_fd`` is +``{mode: {ifo: (n_waveforms, n_freq) complex}}`` (already whitened) and ``data_fd`` is +``{ifo: (n_freq,) complex}`` (whitened data); ``domain`` is the frequency ``Domain`` used +for the inverse FFT. +""" + +from typing import Optional, Sequence, Tuple + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.axes import Axes + +from dingo.gw.domains import Domain + + +def plot_ppd_td( + wf_fd: dict, + data_fd: dict, + domain: Domain, + filename: str = "ppd_td.png", + zoom: Optional[Tuple[float, float]] = None, + axes: Optional[Sequence[Axes]] = None, + plot_data: bool = True, + colors: Optional[Sequence[str]] = None, +) -> np.ndarray: + """Plot the time-domain whitened-strain PPD, one envelope per posterior mode. + + For each detector, inverse-FFTs each mode's whitened waveforms to the time domain with + the merger at t = 0 (segment-midpoint offset ``t0 = T/2``) and shades the pointwise + min/max envelope. The grey whitened-data trace is overlaid once. + + Parameters + ---------- + wf_fd : dict + ``{mode: {ifo: (n_waveforms, n_freq) complex}}`` whitened model waveforms, from + :meth:`Result._compute_ppd`. Each ``mode`` (e.g. ``"dingo"``, ``"dingo-is"``) is + drawn as its own coloured, labelled envelope. + data_fd : dict + ``{ifo: (n_freq,) complex}`` whitened detector data; its keys set the subplot rows. + domain : Domain + Frequency domain used for the inverse FFT (needs ``delta_f``, ``f_max``, ``()``). + filename : str + Output path. Ignored when ``axes`` is supplied (the caller owns saving). + zoom : tuple or None + ``(left, right)`` x-limits in seconds-to-merger. Default ``(-1.0, 0.2)``. + axes : sequence of matplotlib Axes or None + Existing axes (one per detector) to draw onto for composition. When None a new + figure is created and saved to ``filename``. + plot_data : bool + Draw the grey whitened-data trace once. + colors : list[str] or None + One color per mode; defaults to an internal cycle. + + Returns + ------- + numpy.ndarray of the matplotlib Axes drawn onto. + """ + # Envelope colors, one per overlaid posterior mode (first is the single-mode + # default orange); grey for the whitened-data trace. + ppd_colors = ["#DD8452", "#4C72B0", "#55A868", "#C44E52", "#8172B3"] + data_color = "#808080" + + ifos = list(data_fd.keys()) + modes = list(wf_fd.keys()) + if colors is None: + colors = [ppd_colors[i % len(ppd_colors)] for i in range(len(modes))] + + if axes is None: + fig, axes_col = plt.subplots( + len(ifos), 1, figsize=(10, 3 * len(ifos)), sharex=True, squeeze=False + ) + axes = axes_col[:, 0] + else: + fig = None + axes = np.atleast_1d(axes) + + t0 = 1 / (2 * domain.delta_f) # merger at segment midpoint -> t = 0 + phase_shift = np.exp(2j * np.pi * domain() * t0) + + for row, (ax, ifo) in enumerate(zip(axes, ifos)): + for mode, color in zip(modes, colors): + td = [] + for wf in wf_fd[mode][ifo] * phase_shift: + times, x = one_sided_fd_to_td(wf, domain) + td.append(np.real(x)) + td = np.array(td) + ax.fill_between( + times - t0, + td.min(axis=0), + td.max(axis=0), + color=color, + alpha=0.5, + label=mode if row == 0 else None, + ) + if plot_data: + d_times, d_x = one_sided_fd_to_td(data_fd[ifo] * phase_shift, domain) + d_td = np.convolve(np.real(d_x), np.ones(4) / 4, mode="same") + ax.plot( + d_times - t0, + d_td, + color=data_color, + lw=1, + alpha=0.7, + zorder=0, + label="data" if row == 0 else None, + ) + ax.set_ylabel(f"{ifo}\nwhitened strain") + ax.set_xlim(*(zoom if zoom is not None else (-1.0, 0.2))) + + axes[-1].set_xlabel("time to merger (s)") + axes[0].legend(loc="upper left") + if fig is not None: + fig.tight_layout() + fig.savefig(filename, dpi=200, bbox_inches="tight") + plt.close(fig) + return axes + + +def one_sided_fd_to_td( + fd: np.ndarray, domain: Domain +) -> Tuple[np.ndarray, np.ndarray]: + """Inverse-FFT a one-sided (positive-frequency) whitened spectrum to the time domain. + + Zeros the DC bin, then uses ``np.fft.irfft`` (which internally mirrors the conjugate + spectrum) with a ``sqrt(N)`` normalization that matches Dingo's whitening convention + (so whitened noise has unit variance). Returns the time array (``dt = 1 / (2 * f_max)``) + together with the length ``2 * n - 1`` real time series. + """ + fd = np.array(fd, dtype=np.complex128) + fd[0] = 0.0 # zero DC + + n_time = 2 * fd.shape[0] - 1 + td = np.fft.irfft(fd, n=n_time) * np.sqrt(n_time) + times = np.arange(n_time) * (1 / (2 * domain.f_max)) + + return times, td diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index d6880bea5..f1626e759 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -6,8 +6,9 @@ ## Goal Given a dingo `Result` HDF5 file, compute and plot the posterior-predictive -distribution (PPD) of the **whitened strain** in both the time domain and the -frequency domain, per detector, and save the figures as PNGs. +distribution (PPD) of the **whitened strain** in the **time domain**, per +detector, and save the figure as a PNG. (Frequency-domain plotting is out of +scope for this PR.) Derived from `/work/nihargupte/analysis/dingo/asimov/eccentricity/plotting/ production_plots/strain_posterior.ipynb`, but stripped of the @@ -16,47 +17,91 @@ script. ## Deliverable -Two plotting **methods on the GW `Result` class** (`dingo/gw/result.py`), -mirroring the existing `result.plot_corner` API: +Rendering lives in a new GW-specific module **`dingo/gw/utils/plotting.py`**, the +domain-specific counterpart to `dingo/core/utils/plotting.py` (which holds +`plot_corner_multi`). It cannot live in `core/` because PPDs need the GW +likelihood, detector projection, and ASD whitening — `core/` stays +domain-agnostic. The module exposes: ```python -result.plot_ppd_td(filename="ppd_td.png", credible_interval=0.9, - num_waveforms=1000, num_processes=1, zoom=None, ppd=None) -result.plot_ppd_fd(filename="ppd_fd.png", credible_interval=0.9, - num_waveforms=1000, num_processes=1, ppd=None) +plot_ppd_td(wf_fd, data_fd, domain, filename="ppd_td.png", + zoom=None, axes=None, plot_data=True, colors=None) +one_sided_fd_to_td(fd, domain) # IFFT helper, also here ``` -They live on the GW `Result` (subclass of `CoreResult`), **not** `core/`, -because they need the GW likelihood, detector projection, and ASD whitening — -`core/` stays domain-agnostic. +`wf_fd` is `{mode: {ifo: (n_waveforms, n_freq) complex}}` (whitened model +waveforms) and `data_fd` is `{ifo: (n_freq,) complex}` (whitened data) — both +straight from `_compute_ppd`. `plot_ppd_td` draws **one coloured min/max envelope +per `mode`** (`"dingo"`, and `"dingo-is"` when present), labelled via +`_MODE_LABELS`, and overlays the grey whitened-data trace once. `axes=None` +creates a figure and saves to `filename`; passing `axes` draws onto them and +skips saving, for composition. -Shared expensive computation lives in a private -`_compute_ppd(self, credible_interval, num_waveforms, num_processes, seed)` -(used by both methods, so extraction is justified) returning -`{domain, ifos, wf_fd, data_fd}`. Each public method accepts an optional -precomputed `ppd=` dict so a caller can generate the draws once and render both -plots without regenerating waveforms. +The GW `Result` (subclass of `CoreResult`) keeps a thin convenience method that +mirrors the `result.plot_corner` API and delegates to the module: + +```python +result.plot_ppd_td(filename="ppd_td.png", credible_interval=0.9, + num_waveforms=1000, num_processes=1, zoom=None, ppd=None) +``` -Demo run target: GW230709_122727 EAS importance-sampling result -(`prod_o4a/working/S230709bi/Exp19_more_points_1/result/ -Exp19_more_points_1_data0_1372940865-2_importance_sampling.hdf5`). A few-line -script loads the `Result` and calls both methods, producing -`GW230709_ppd_td.png` and `GW230709_ppd_fd.png`. +Like `plot_corner`, both the **Dingo** posterior (credible set from the raw +network samples) and — when the result is importance-sampled (`"weights"` +present) — the **Dingo-IS** posterior are overlaid on one figure. + +The shared expensive computation stays a **private `Result` method** +`_compute_ppd(self, credible_interval, num_waveforms, num_processes, seed)`, +returning the tuple `(wf_fd, data_fd)`. It builds a credible set and generates +waveforms for **each** available posterior internally (no `weighting` argument): +`"dingo"` always, `"dingo-is"` iff importance-sampled. `plot_ppd_td` accepts a +precomputed `ppd=(wf_fd, data_fd)` to render without regenerating waveforms; the +`Result` method supplies `domain` from `self.domain`. + +Robustness (prior filter): the normalizing-flow proposal `q(theta)` is a smooth +density over a box and leaks a small fraction of draws outside the prior — with +finite `log_prob` but unphysical (e.g. tiny negative spin magnitudes). These are +exactly the draws that fail waveform generation. `_compute_ppd` therefore +restricts to the prior support up front (`self.prior.ln_prob(...)`, keeping +`log_prior` finite) — the same in-prior restriction `importance_sample` applies — +before generating anything. This removes the failures *and* keeps unphysical +waveforms out of the envelope, and it is a strict improvement for the `"dingo"` +mode (which otherwise gives those out-of-prior draws equal mass `1/N`). For an +importance-sampled result the filter leaves **no** generation failures at all: +IS already generated every in-prior sample successfully (or it would have crashed +at `core/likelihood.py`, which has no `try/except`). Measured on the GW150914 +fixture: 43/5000 draws (0.86%) are out-of-prior; **all** generation failures were +among them, **zero** in-prior draws failed. + +Because the prior filter removes every observed failure, there is no per-sample +failure handling: `_compute_ppd` calls `self.likelihood.signal` directly (a bound +method, picklable for the pool, as with `log_likelihood` in `core/likelihood.py`) +and the envelope uses plain `min`/`max`. A genuine in-prior model failure (an +approximant interpolation edge on a *not-yet-importance-sampled* result) is left to +raise, surfacing the offending sample rather than being silently dropped. + +Dev fixture: GW150914 dingo-ci result (IMRPhenomXPHM, no `pyseobnr` needed), +`.../dingo-ci/outdir_GW150914/result/ +GW150914_data0_1126259462-4_importance_sampling_part1.hdf5`, rendered via the +throwaway `ppd_dev.py`. ## Pipeline (per detector) 1. `result = Result(file_name=...)`; `result._build_likelihood()`. `domain = result.domain.base_domain if hasattr(result.domain, "base_domain") else result.domain`. -2. **Sample selection.** Build the `credible_interval` (default 0.9) credible - set: normalize the importance/posterior weights, sort descending, and keep - the highest-weight samples whose cumulative weight ≤ `credible_interval` - (the X% credible region — same definition as the notebook's - `select_from_n_percent_interval`). Then draw up to `num_waveforms` rows - **uniformly at random** from that set (seeded `np.random.default_rng`) so the - envelope represents the whole interval rather than only the peak (the - notebook took the top-N highest-weight rows, biasing toward the peak — we - diverge here). +2. **Sample selection.** Build the `credible_interval` (default 0.9) highest- + posterior-density credible set using the per-sample probability **mass** that + is correct for how the samples were drawn (no assumed reference density): for + `"dingo"` the raw samples are equal-mass draws from q(θ) — mass `1/N`, ranked + by `log_prob`; for `"dingo-is"` the target is the true posterior, so the mass + is the self-normalized importance weight (the `weights` column, ∝ p/q), ranked + by `log_likelihood + log_prior`. Sort by rank descending and keep the top + samples whose cumulative mass ≤ `credible_interval` (standard HPD estimator — + the same weighting `plot_corner` feeds to corner). Then draw up to + `num_waveforms` rows **uniformly at random** from that set (seeded + `np.random.default_rng`) so the envelope spans the whole interval rather than + only the peak (the notebook took the top-N highest-weight rows, biasing toward + the peak — we diverge here). 3. Generate whitened FD waveforms: `apply_func_with_multiprocessing(result.likelihood.signal, samples_df, num_processes)`. Each `wf["waveform"][ifo]` is **already whitened** @@ -66,36 +111,33 @@ script loads the `Result` and calls both methods, producing ### Time offset (hidden from the user) The merger always lands at **t = 0**. Internally each waveform is phase-shifted -by a fixed offset derived from the data segment (`t0 = 1 / domain.delta_f`, the -segment duration `T`, which places the dingo coalescence at the segment edge), -via `wf *= exp(2πi·f·t0)`, then the time axis has `t0` subtracted. `t0` is a -private detail — not a method argument. The implementation verifies the demo -merger sits at 0 and adjusts the offset derivation if not. +by a fixed offset derived from the data segment (`t0 = 1 / (2·domain.delta_f)`, +i.e. `T/2`, the segment **midpoint** where dingo centers the coalescence), via +`wf *= exp(2πi·f·t0)`, then the time axis has `t0` subtracted. `t0` is a private +detail — not a method argument. This generalizes the source notebook's hardcoded +`CENTRAL_TIME = 8` s (which was `T/2` for those `T = 16` s segments). Verified +empirically: for the GW150914 CI result (`T = 8` s, `t0 = 4`) the whitened-strain +peak lands at t ≈ +0.02 s in both H1 and L1. ## Time domain (`_ppd_td.png`, one row per detector) - Each draw: phase-shift by `t0`, `one_sided_fd_to_td`, keep real part. - Band = pointwise **min/max envelope** across the selected draws - (`np.nanmin` / `np.nanmax`, ignoring NaN waveforms), drawn with - `fill_between`. + (`td.min` / `td.max`), drawn with `fill_between`. - Overlay: whitened data trace (`√(4Δf)·context_waveform/asd`, same shift + IFFT, light sliding-window average), grey. -- x-axis: time to merger (s), linear, zoomed via `zoom=` (default - `(-0.7, 0.1)` for GW230709). y-axis: whitened strain (dimensionless, σ units). - -## Frequency domain (`_ppd_fd.png`, one row per detector) +- x-axis: time to merger (s), linear, zoomed via `zoom=` (method default + `(-1.0, 0.2)`; `(-0.7, 0.1)` works well for GW230709). y-axis: whitened strain + (dimensionless, σ units). -- Per draw: `|whitened h(f)|` (magnitude of the complex whitened FD waveform). -- Band = min/max envelope of `|whitened h(f)|` across the selected draws. -- Overlay: `|whitened data|` = `|√(4Δf)·context_waveform/asd|` — the noisy - ~O(1) floor the signal rises above. -- x-axis: frequency [Hz], **log scale**. y-axis: dimensionless whitened - amplitude. +Frequency-domain PPD plotting is deferred to a follow-up (a per-bin whitened FD +view buries the signal below the unit-variance noise floor; the useful FD view is +an amplitude-spectral-density / log-log plot, out of scope here). ## Reused logic -`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook (added as -a module-level function in `dingo/gw/result.py`): it builds the full Hermitian +`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook (a +module-level function in `dingo/gw/utils/plotting.py`): it builds the full Hermitian two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` (normalization tied to the whitening so noise has unit variance), `dt = 1/(2·f_max)`, output length `2n-1`. Correct and whitening-consistent. @@ -103,14 +145,15 @@ to the whitening so noise has unit variance), `dt = 1/(2·f_max)`, output length ## Testing / determinism - Seeded `np.random.default_rng` for resampling → deterministic figures. -- `tests/gw/inference/test_ppd.py`: assert `one_sided_fd_to_td` round-trips — - a one-sided spectrum with a single populated frequency bin returns a real - cosine of the expected period (exercises the Hermitian-mirror + normalization - logic, the one piece of non-trivial new math). +- `tests/gw/test_ppd.py`: assert `one_sided_fd_to_td` (imported from + `dingo.gw.utils.plotting`) round-trips — a one-sided spectrum with a single populated + frequency bin returns a real cosine of the expected period (exercises the + Hermitian-mirror + normalization logic, the one piece of non-trivial new math). ## Boundaries / conventions -- Methods on the GW `Result` in `gw/result.py`; imports `gw/` + `core/` - (allowed direction). No new classes, no new dependencies. +- Rendering in `gw/utils/plotting.py`, thin `Result.plot_ppd_td` + `_compute_ppd` in + `gw/result.py`; imports `gw/` + `core/` (allowed direction). No new classes, + no new dependencies. - Whitening conventions for data (`√(4Δf)/asd`) and model (likelihood's `/asd/noise_std`, `noise_std = 1/√(4Δf)`) are identical and must be preserved. diff --git a/tests/gw/test_ppd.py b/tests/gw/test_ppd.py new file mode 100644 index 000000000..fc0c24c73 --- /dev/null +++ b/tests/gw/test_ppd.py @@ -0,0 +1,32 @@ +"""Tests for the strain-PPD helpers in :mod:`dingo.gw.utils.plotting`.""" + +from types import SimpleNamespace + +import numpy as np + +from dingo.gw.utils.plotting import one_sided_fd_to_td + + +def test_one_sided_fd_to_td_single_tone() -> None: + """A one-sided spectrum with a single populated bin k inverse-FFTs to a real cosine + of frequency k * delta_f, with the documented 2 * n - 1 output length.""" + n = 65 + delta_f = 1.0 + f_max = (n - 1) * delta_f # one_sided_fd_to_td only reads domain.f_max + domain = SimpleNamespace(f_max=f_max) + + k = 8 + fd = np.zeros(n, dtype=np.complex128) + fd[k] = 1.0 + + times, td = one_sided_fd_to_td(fd, domain) + + # Hermitian spectrum -> real time series, and the documented length. + assert len(td) == 2 * n - 1 + assert np.max(np.abs(td.imag)) < 1e-9 + + # Dominant frequency of the recovered tone is k * delta_f. + dt = times[1] - times[0] + rfreqs = np.fft.rfftfreq(len(td), d=dt) + recovered = rfreqs[np.argmax(np.abs(np.fft.rfft(td.real)))] + assert abs(recovered - k * delta_f) < 1.0 # within one frequency bin From 1b45cd2f7ad4d4b3ff7fd0a21aa6dd09c2fbe933 Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 12:45:59 +0200 Subject: [PATCH 05/11] chore: gitignore local/ scratch dir for generated images Generated PPD plots and other throwaway artifacts go in a local/ directory at the worktree root, kept out of version control. Co-Authored-By: Claude Opus 4.8 (1M context) --- .gitignore | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 3a5bb4048..e76100e45 100644 --- a/.gitignore +++ b/.gitignore @@ -151,4 +151,7 @@ dmypy.json # directories with temporary files **/tmp_files/* -dingo-venv/ \ No newline at end of file +dingo-venv/ + +# scratch/output dir for generated images and throwaway artifacts +local/ \ No newline at end of file From 5f2e9b1e117a7ebe8ab6c3d77f4303c6b46aaa92 Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 13:47:59 +0200 Subject: [PATCH 06/11] fix(gw): reinstate per-sample waveform guard in PPD The prior filter removes out-of-prior draws but does NOT guarantee every in-prior draw generates. A waveform model can fail on an interior parameter, and this bites even for importance-sampled results whenever the waveform backend differs from the one used at IS time: re-plotting a stored result with a newer/older pyseobnr / lalsimulation is enough to make an IS-accepted draw fail. Concrete case: the S231001aq (SEOBNRv5EHM) external result has 2/200 in-prior credible-set draws that crash inside pyseobnr's eccentric CubicSpline ("x must be strictly increasing") under pyseobnr 0.3.6 / lalsimulation 6.2.1 -- taking down the whole plot. A previous revision had removed this guard on the mistaken assumption that the prior filter left no in-prior failures for IS results. Restore safe_signal (module-level, picklable) and generate each draw through it: failed draws are dropped from the min/max envelope with a warning giving the count, and only an all-failed batch raises. Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/injection.py | 18 ++++++++ dingo/gw/result.py | 43 ++++++++++++++----- .../2026-06-30-strain-ppd-plotting-design.md | 43 +++++++++++-------- 3 files changed, 75 insertions(+), 29 deletions(-) diff --git a/dingo/gw/injection.py b/dingo/gw/injection.py index 87f7a0660..a71932ca7 100644 --- a/dingo/gw/injection.py +++ b/dingo/gw/injection.py @@ -1,3 +1,5 @@ +from typing import Optional + import numpy as np from bilby.gw.detector import InterferometerList from torchvision.transforms import Compose @@ -442,3 +444,19 @@ def random_injection(self): k: float(v) for k, v in theta.items() } # Some parameters are np.float64 return self.injection(theta) + + +def safe_signal(gw_signal: "GWSignal", theta: dict) -> Optional[dict]: + """Call ``gw_signal.signal(theta)``, returning ``None`` if generation fails. + + Even in-prior posterior draws can hit waveform-model edge cases -- e.g. SEOBNRv5EHM's + eccentric interpolation failing on an interior parameter, or simply a waveform backend + (pyseobnr / lalsimulation) that differs from the one used when the result was produced. + Returning ``None`` lets callers such as :meth:`dingo.gw.result.Result._compute_ppd` drop + the offending draw rather than aborting a whole :class:`multiprocessing.Pool` batch. + Defined at module level (not a lambda/nested function) so it is picklable for the pool. + """ + try: + return gw_signal.signal(theta) + except Exception: + return None diff --git a/dingo/gw/result.py b/dingo/gw/result.py index eb7c25e88..f865c6fb5 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -1,5 +1,7 @@ import copy import time +import warnings +from functools import partial from typing import Optional, Tuple import numpy as np @@ -18,6 +20,7 @@ from dingo.gw.domains import MultibandedFrequencyDomain from dingo.gw.domains import build_domain from dingo.gw.gwutils import get_extrinsic_prior_dict +from dingo.gw.injection import safe_signal from dingo.gw.likelihood import StationaryGaussianGWLikelihood from dingo.gw.utils.plotting import plot_ppd_td as _plot_ppd_td from dingo.gw.prior import build_prior_with_defaults @@ -841,12 +844,16 @@ def _compute_ppd( Prior filter. Samples are first restricted to the prior support (finite ``log_prior``). The flow proposal is a smooth density over a box and leaks a small fraction of draws outside the prior (e.g. tiny negative spin magnitudes); these are unphysical, carry - zero target mass, and are the samples that fail waveform generation. This is exactly - the in-prior restriction that :meth:`importance_sample` applies, so for an importance- - sampled result -- whose in-prior samples were all generated successfully during IS -- - the filter is guaranteed to leave no generation failures. A genuine in-prior model - failure (e.g. an approximant interpolation edge on a not-yet-importance-sampled result) - is left to raise, surfacing the offending sample rather than being silently dropped. + zero target mass, and account for most waveform-generation failures. This is exactly + the in-prior restriction that :meth:`importance_sample` applies. + + Per-sample guard. The prior filter is not sufficient on its own: a waveform model can + fail on an *interior* parameter (e.g. SEOBNRv5EHM's eccentric ``CubicSpline`` on a + non-monotonic omega array), and this happens even for importance-sampled results + whenever the waveform backend here differs from the one used at IS time -- re-plotting + a stored result with a newer/older ``pyseobnr``/``lalsimulation`` is enough. So each + draw is generated through ``safe_signal``: failures are dropped from the envelope with + a warning, and only an all-failed batch raises. """ if getattr(self, "likelihood", None) is None: self._build_likelihood() @@ -912,14 +919,30 @@ def _compute_ppd( ) theta = self.samples.iloc[credible_idx] - # In-prior samples generate cleanly (guaranteed for an importance-sampled - # result, whose in-prior draws were all generated during IS); a genuine model - # failure here should surface, not be silently dropped. + # The prior filter removes out-of-prior draws, but it does not guarantee every + # in-prior draw generates: a waveform model can fail on an interior parameter + # (e.g. SEOBNRv5EHM's eccentric CubicSpline hitting a non-monotonic omega array). + # This bites even for importance-sampled results whenever the waveform backend + # here differs from the one used at IS time -- re-plotting a stored result with a + # newer/older pyseobnr/lalsimulation is enough to make an IS-accepted draw fail. + # So guard per-sample: drop failures, warn with the count, raise only if all fail. signals = apply_func_with_multiprocessing( - self.likelihood.signal, + partial(safe_signal, self.likelihood), theta, num_processes=num_processes, ) + n_failed = sum(s is None for s in signals) + signals = [s for s in signals if s is not None] + if not signals: + raise RuntimeError( + f"All '{mode}' PPD waveform generations failed; cannot build the PPD." + ) + if n_failed: + warnings.warn( + f"{n_failed} in-prior '{mode}' PPD waveform(s) failed to generate " + "and were dropped from the envelope.", + stacklevel=2, + ) wf_fd[mode] = { ifo: np.array([s["waveform"][ifo] for s in signals]) for ifo in ifos } diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index f1626e759..bc24c0c2a 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -59,25 +59,30 @@ precomputed `ppd=(wf_fd, data_fd)` to render without regenerating waveforms; the Robustness (prior filter): the normalizing-flow proposal `q(theta)` is a smooth density over a box and leaks a small fraction of draws outside the prior — with -finite `log_prob` but unphysical (e.g. tiny negative spin magnitudes). These are -exactly the draws that fail waveform generation. `_compute_ppd` therefore -restricts to the prior support up front (`self.prior.ln_prob(...)`, keeping -`log_prior` finite) — the same in-prior restriction `importance_sample` applies — -before generating anything. This removes the failures *and* keeps unphysical -waveforms out of the envelope, and it is a strict improvement for the `"dingo"` -mode (which otherwise gives those out-of-prior draws equal mass `1/N`). For an -importance-sampled result the filter leaves **no** generation failures at all: -IS already generated every in-prior sample successfully (or it would have crashed -at `core/likelihood.py`, which has no `try/except`). Measured on the GW150914 -fixture: 43/5000 draws (0.86%) are out-of-prior; **all** generation failures were -among them, **zero** in-prior draws failed. - -Because the prior filter removes every observed failure, there is no per-sample -failure handling: `_compute_ppd` calls `self.likelihood.signal` directly (a bound -method, picklable for the pool, as with `log_likelihood` in `core/likelihood.py`) -and the envelope uses plain `min`/`max`. A genuine in-prior model failure (an -approximant interpolation edge on a *not-yet-importance-sampled* result) is left to -raise, surfacing the offending sample rather than being silently dropped. +finite `log_prob` but unphysical (e.g. tiny negative spin magnitudes). These +account for most waveform-generation failures. `_compute_ppd` therefore restricts +to the prior support up front (`self.prior.ln_prob(...)`, keeping `log_prior` +finite) — the same in-prior restriction `importance_sample` applies — before +generating anything. This removes those failures *and* keeps unphysical waveforms +out of the envelope, and it is a strict improvement for the `"dingo"` mode (which +otherwise gives those out-of-prior draws equal mass `1/N`). Measured on the +GW150914 fixture: 43/5000 draws (0.86%) are out-of-prior; all failures there were +among them, zero in-prior draws failed. + +Robustness (per-sample guard): the prior filter is **not** sufficient on its own. +A waveform model can fail on an *interior* parameter, and this bites even for +importance-sampled results whenever the waveform backend differs from the one used +at IS time — re-plotting a stored result with a newer/older `pyseobnr` / +`lalsimulation` is enough to make an IS-accepted draw fail. Observed on the +`S231001aq` (SEOBNRv5EHM) external result: 2/200 in-prior credible-set draws crash +inside pyseobnr's eccentric `CubicSpline` (`x must be strictly increasing`) under +`pyseobnr 0.3.6` / `lalsimulation 6.2.1`. So each draw is generated through the +module-level picklable `safe_signal` (`injection.py`): failures are dropped from +the `min`/`max` envelope with a warning giving the count, and only an all-failed +batch raises. (An earlier revision removed this guard on the mistaken assumption +that the prior filter left no in-prior failures for IS results; `S231001aq` is the +counterexample — the "IS already generated these" guarantee only holds for the +exact same waveform backend.) Dev fixture: GW150914 dingo-ci result (IMRPhenomXPHM, no `pyseobnr` needed), `.../dingo-ci/outdir_GW150914/result/ From 3470e06a00973e7df234542b83dbd57beb57694a Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 13:53:25 +0200 Subject: [PATCH 07/11] style(gw): lighter default PPD envelope alpha (0.5 -> 0.3) The two overlaid mode envelopes (dingo, dingo-is) plus the data trace read more clearly with a lighter fill, especially where they overlap. Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/utils/plotting.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/dingo/gw/utils/plotting.py b/dingo/gw/utils/plotting.py index ba6d1b0af..ea2454705 100644 --- a/dingo/gw/utils/plotting.py +++ b/dingo/gw/utils/plotting.py @@ -96,7 +96,7 @@ def plot_ppd_td( td.min(axis=0), td.max(axis=0), color=color, - alpha=0.5, + alpha=0.3, label=mode if row == 0 else None, ) if plot_data: From c24d9339073feda871906930e8d2dc7fbac77025 Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 13:57:25 +0200 Subject: [PATCH 08/11] feat(gw): overlay max-probability waveform on PPD plots For each mode, generate and draw the single maximum-probability waveform as a solid line over the min/max band: max log_prob for "dingo", max (log_likelihood + log_prior) for "dingo-is". This gives a clean representative chirp on top of the band -- especially helpful for weaker/eccentric events where the min/max envelope of time/phase-misaligned draws otherwise reads as a fuzzy packet. _compute_ppd now returns (wf_fd, data_fd, map_fd); plot_ppd_td takes an optional map_fd and draws each present mode's MAP line in the mode's colour. Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/result.py | 42 ++++++++++++++++++++++++++------------ dingo/gw/utils/plotting.py | 20 ++++++++++++++---- 2 files changed, 45 insertions(+), 17 deletions(-) diff --git a/dingo/gw/result.py b/dingo/gw/result.py index f865c6fb5..056f6308b 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -792,7 +792,7 @@ def _compute_ppd( num_waveforms: int = 1000, num_processes: int = 1, seed: int = RANDOM_STATE, - ) -> Tuple[dict, dict]: + ) -> Tuple[dict, dict, dict]: """Generate whitened waveforms for the time-domain strain PPD. Builds a ``credible_interval`` highest-posterior-density credible set for each @@ -822,8 +822,13 @@ def _compute_ppd( draw is dropped, so ``n_kept`` may be slightly below the number requested. data_fd : dict ``{ifo: complex ndarray (n_freq,)}`` whitened detector data (shared across modes). + map_fd : dict + ``{mode: {ifo: complex ndarray (n_freq,)}}`` the single maximum-probability + whitened waveform per mode -- max ``log_prob`` for ``"dingo"``, max + ``log_likelihood + log_prior`` for ``"dingo-is"`` -- drawn as a line over the band. + A mode is absent iff its MAP draw failed generation. - Both are consumed by :func:`dingo.gw.utils.plotting.plot_ppd_td` (together with + All three are consumed by :func:`dingo.gw.utils.plotting.plot_ppd_td` (together with ``self.domain``). Notes @@ -891,8 +896,10 @@ def _compute_ppd( ) # One credible set per available posterior: "dingo" always; "dingo-is" when the - # result is importance-sampled (weights present). + # result is importance-sampled (weights present). map_fd holds the single + # maximum-probability waveform per mode (drawn as a line over the band). wf_fd = {} + map_fd = {} for mode in ("dingo", "dingo-is"): if mode == "dingo": mass = in_prior.astype(float) @@ -947,12 +954,21 @@ def _compute_ppd( ifo: np.array([s["waveform"][ifo] for s in signals]) for ifo in ifos } + # Maximum-probability waveform for this mode: the single highest-ranked + # (in-prior) draw -- max log_prob for "dingo", max (log_likelihood + log_prior) + # for "dingo-is". Gives a clean representative chirp over the min/max band. + map_signal = safe_signal( + self.likelihood, self.samples.iloc[order[0]].to_dict() + ) + if map_signal is not None: + map_fd[mode] = {ifo: map_signal["waveform"][ifo] for ifo in ifos} + if wf_fd["dingo"][ifos[0]].shape[1] != len(domain()): raise ValueError( "Whitened waveform length does not match the inverse-FFT domain; " "time-domain PPD requires a uniform frequency domain." ) - return wf_fd, data_fd + return wf_fd, data_fd, map_fd def plot_ppd_td( self, @@ -961,19 +977,19 @@ def plot_ppd_td( num_waveforms: int = 1000, num_processes: int = 1, zoom: Optional[Tuple[float, float]] = None, - ) -> Tuple[dict, dict]: + ) -> Tuple[dict, dict, dict]: """Plot the time-domain whitened-strain posterior-predictive distribution. For each detector, inverse-FFTs whitened waveforms to the time domain with the merger - at t = 0 and shades the pointwise min/max envelope; the whitened detector data is - overlaid in grey. Mirroring :meth:`plot_corner`, both the **Dingo** posterior and (when - the result is importance-sampled) the **Dingo-IS** posterior are overlaid on one - figure. + at t = 0 and shades the pointwise min/max envelope, with the maximum-probability + waveform of each mode drawn as a line over it; the whitened detector data is overlaid + in grey. Mirroring :meth:`plot_corner`, both the **Dingo** posterior and (when the + result is importance-sampled) the **Dingo-IS** posterior are overlaid on one figure. ``zoom`` is the (left, right) x-limit in seconds-to-merger; defaults to (-1.0, 0.2). - Returns the ``(wf_fd, data_fd)`` tuple that was plotted. + Returns the ``(wf_fd, data_fd, map_fd)`` tuple that was plotted. """ - wf_fd, data_fd = self._compute_ppd( + wf_fd, data_fd, map_fd = self._compute_ppd( credible_interval, num_waveforms, num_processes ) domain = ( @@ -981,5 +997,5 @@ def plot_ppd_td( if hasattr(self.domain, "base_domain") else self.domain ) - _plot_ppd_td(wf_fd, data_fd, domain, filename=filename, zoom=zoom) - return wf_fd, data_fd + _plot_ppd_td(wf_fd, data_fd, domain, map_fd=map_fd, filename=filename, zoom=zoom) + return wf_fd, data_fd, map_fd diff --git a/dingo/gw/utils/plotting.py b/dingo/gw/utils/plotting.py index ea2454705..e57264939 100644 --- a/dingo/gw/utils/plotting.py +++ b/dingo/gw/utils/plotting.py @@ -6,9 +6,10 @@ when available, ``"dingo-is"``), mirroring how ``result.plot_corner`` shows both. Inputs come straight from ``Result._compute_ppd``: ``wf_fd`` is -``{mode: {ifo: (n_waveforms, n_freq) complex}}`` (already whitened) and ``data_fd`` is -``{ifo: (n_freq,) complex}`` (whitened data); ``domain`` is the frequency ``Domain`` used -for the inverse FFT. +``{mode: {ifo: (n_waveforms, n_freq) complex}}`` (already whitened), ``data_fd`` is +``{ifo: (n_freq,) complex}`` (whitened data), and ``map_fd`` is +``{mode: {ifo: (n_freq,) complex}}`` (the per-mode maximum-probability waveform, drawn as +a line); ``domain`` is the frequency ``Domain`` used for the inverse FFT. """ from typing import Optional, Sequence, Tuple @@ -24,6 +25,7 @@ def plot_ppd_td( wf_fd: dict, data_fd: dict, domain: Domain, + map_fd: Optional[dict] = None, filename: str = "ppd_td.png", zoom: Optional[Tuple[float, float]] = None, axes: Optional[Sequence[Axes]] = None, @@ -34,7 +36,8 @@ def plot_ppd_td( For each detector, inverse-FFTs each mode's whitened waveforms to the time domain with the merger at t = 0 (segment-midpoint offset ``t0 = T/2``) and shades the pointwise - min/max envelope. The grey whitened-data trace is overlaid once. + min/max envelope, with each mode's maximum-probability waveform drawn as a line over it. + The grey whitened-data trace is overlaid once. Parameters ---------- @@ -46,6 +49,10 @@ def plot_ppd_td( ``{ifo: (n_freq,) complex}`` whitened detector data; its keys set the subplot rows. domain : Domain Frequency domain used for the inverse FFT (needs ``delta_f``, ``f_max``, ``()``). + map_fd : dict or None + ``{mode: {ifo: (n_freq,) complex}}`` maximum-probability waveform per mode, drawn as + a solid line in the mode's colour over its band. Modes absent from the dict (or a + ``None`` argument) get no line. filename : str Output path. Ignored when ``axes`` is supplied (the caller owns saving). zoom : tuple or None @@ -62,6 +69,7 @@ def plot_ppd_td( ------- numpy.ndarray of the matplotlib Axes drawn onto. """ + map_fd = map_fd or {} # Envelope colors, one per overlaid posterior mode (first is the single-mode # default orange); grey for the whitened-data trace. ppd_colors = ["#DD8452", "#4C72B0", "#55A868", "#C44E52", "#8172B3"] @@ -99,6 +107,10 @@ def plot_ppd_td( alpha=0.3, label=mode if row == 0 else None, ) + # Maximum-probability waveform for this mode, as a solid line over the band. + if mode in map_fd: + m_times, m_x = one_sided_fd_to_td(map_fd[mode][ifo] * phase_shift, domain) + ax.plot(m_times - t0, np.real(m_x), color=color, lw=1.0, alpha=0.9) if plot_data: d_times, d_x = one_sided_fd_to_td(data_fd[ifo] * phase_shift, domain) d_td = np.convolve(np.real(d_x), np.ones(4) / 4, mode="same") From 1f5b16e1d07eff49b3d58e059c04b3502095786b Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 15:33:44 +0200 Subject: [PATCH 09/11] style(gw): stack PPD min/max panels per (mode, detector) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Split the overlaid Dingo / Dingo-IS envelopes into separate stacked panels -- one per (mode, detector), Dingo on top and Dingo-IS below -- each with its min/max band, MAP waveform line, grey data trace, and an in-axes "{mode} · {ifo}" label. Drops the unused axes/plot_data/colors composition args and the legend. Records (in the design doc) that shading by the full per-time density p(s|t) was prototyped but rejected: mathematically the right object and clean on GW150914, but messy on real eccentric O4 events (phase decorrelation + shot noise near merger). Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/result.py | 10 +- dingo/gw/utils/plotting.py | 144 ++++++++---------- .../2026-06-30-strain-ppd-plotting-design.md | 35 +++-- 3 files changed, 93 insertions(+), 96 deletions(-) diff --git a/dingo/gw/result.py b/dingo/gw/result.py index 056f6308b..68a53187a 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -980,11 +980,11 @@ def plot_ppd_td( ) -> Tuple[dict, dict, dict]: """Plot the time-domain whitened-strain posterior-predictive distribution. - For each detector, inverse-FFTs whitened waveforms to the time domain with the merger - at t = 0 and shades the pointwise min/max envelope, with the maximum-probability - waveform of each mode drawn as a line over it; the whitened detector data is overlaid - in grey. Mirroring :meth:`plot_corner`, both the **Dingo** posterior and (when the - result is importance-sampled) the **Dingo-IS** posterior are overlaid on one figure. + For each mode and detector, inverse-FFTs whitened waveforms to the time domain with + the merger at t = 0 and shades the pointwise min/max envelope, with the + maximum-probability waveform drawn as a line over it and the whitened data in grey. + Panels are stacked vertically, one per (mode, detector): the **Dingo** posterior on + top and, when the result is importance-sampled, the **Dingo-IS** posterior below. ``zoom`` is the (left, right) x-limit in seconds-to-merger; defaults to (-1.0, 0.2). Returns the ``(wf_fd, data_fd, map_fd)`` tuple that was plotted. diff --git a/dingo/gw/utils/plotting.py b/dingo/gw/utils/plotting.py index e57264939..343ff66a7 100644 --- a/dingo/gw/utils/plotting.py +++ b/dingo/gw/utils/plotting.py @@ -2,8 +2,8 @@ GW-specific counterpart to :mod:`dingo.core.utils.plotting`: renders the whitened-strain PPD envelopes produced by :meth:`dingo.gw.result.Result._compute_ppd` in the time domain. -:func:`plot_ppd_td` overlays one envelope per posterior mode in ``wf_fd`` (``"dingo"`` and, -when available, ``"dingo-is"``), mirroring how ``result.plot_corner`` shows both. +:func:`plot_ppd_td` draws one min/max envelope panel per (posterior mode, detector), +stacked vertically -- ``"dingo"`` panels on top and, when available, ``"dingo-is"`` below. Inputs come straight from ``Result._compute_ppd``: ``wf_fd`` is ``{mode: {ifo: (n_waveforms, n_freq) complex}}`` (already whitened), ``data_fd`` is @@ -12,14 +12,17 @@ a line); ``domain`` is the frequency ``Domain`` used for the inverse FFT. """ -from typing import Optional, Sequence, Tuple +from typing import Optional, Tuple import numpy as np from matplotlib import pyplot as plt -from matplotlib.axes import Axes from dingo.gw.domains import Domain +# One colour per posterior mode (first is the single-mode default orange); grey for data. +_PPD_COLORS = ["#DD8452", "#4C72B0", "#55A868", "#C44E52", "#8172B3"] +_DATA_COLOR = "#808080" + def plot_ppd_td( wf_fd: dict, @@ -28,110 +31,91 @@ def plot_ppd_td( map_fd: Optional[dict] = None, filename: str = "ppd_td.png", zoom: Optional[Tuple[float, float]] = None, - axes: Optional[Sequence[Axes]] = None, - plot_data: bool = True, - colors: Optional[Sequence[str]] = None, ) -> np.ndarray: - """Plot the time-domain whitened-strain PPD, one envelope per posterior mode. + """Plot the time-domain whitened-strain PPD min/max envelope, one panel per (mode, ifo). - For each detector, inverse-FFTs each mode's whitened waveforms to the time domain with - the merger at t = 0 (segment-midpoint offset ``t0 = T/2``) and shades the pointwise - min/max envelope, with each mode's maximum-probability waveform drawn as a line over it. - The grey whitened-data trace is overlaid once. + For each mode and detector, inverse-FFTs the whitened draws to the time domain (merger at + t = 0, segment-midpoint offset ``t0 = T/2``), shades the pointwise min/max envelope, and + draws that mode's maximum-probability waveform as a line over it, all above the grey + whitened-data trace. Panels are stacked vertically -- the **Dingo** posterior on top and, + when the result is importance-sampled, **Dingo-IS** below. Parameters ---------- wf_fd : dict ``{mode: {ifo: (n_waveforms, n_freq) complex}}`` whitened model waveforms, from - :meth:`Result._compute_ppd`. Each ``mode`` (e.g. ``"dingo"``, ``"dingo-is"``) is - drawn as its own coloured, labelled envelope. + :meth:`Result._compute_ppd` (``mode`` is ``"dingo"`` and/or ``"dingo-is"``). data_fd : dict - ``{ifo: (n_freq,) complex}`` whitened detector data; its keys set the subplot rows. + ``{ifo: (n_freq,) complex}`` whitened detector data; its keys set the detectors. domain : Domain Frequency domain used for the inverse FFT (needs ``delta_f``, ``f_max``, ``()``). map_fd : dict or None ``{mode: {ifo: (n_freq,) complex}}`` maximum-probability waveform per mode, drawn as - a solid line in the mode's colour over its band. Modes absent from the dict (or a - ``None`` argument) get no line. + a solid line in the mode's colour. Modes absent from the dict (or ``None``) get no line. filename : str - Output path. Ignored when ``axes`` is supplied (the caller owns saving). + Output path for the saved figure. zoom : tuple or None ``(left, right)`` x-limits in seconds-to-merger. Default ``(-1.0, 0.2)``. - axes : sequence of matplotlib Axes or None - Existing axes (one per detector) to draw onto for composition. When None a new - figure is created and saved to ``filename``. - plot_data : bool - Draw the grey whitened-data trace once. - colors : list[str] or None - One color per mode; defaults to an internal cycle. Returns ------- - numpy.ndarray of the matplotlib Axes drawn onto. + numpy.ndarray of the matplotlib Axes drawn onto (one per stacked (mode, detector) panel). """ map_fd = map_fd or {} - # Envelope colors, one per overlaid posterior mode (first is the single-mode - # default orange); grey for the whitened-data trace. - ppd_colors = ["#DD8452", "#4C72B0", "#55A868", "#C44E52", "#8172B3"] - data_color = "#808080" - ifos = list(data_fd.keys()) modes = list(wf_fd.keys()) - if colors is None: - colors = [ppd_colors[i % len(ppd_colors)] for i in range(len(modes))] - - if axes is None: - fig, axes_col = plt.subplots( - len(ifos), 1, figsize=(10, 3 * len(ifos)), sharex=True, squeeze=False - ) - axes = axes_col[:, 0] - else: - fig = None - axes = np.atleast_1d(axes) + colors = {mode: _PPD_COLORS[i % len(_PPD_COLORS)] for i, mode in enumerate(modes)} + zoom = zoom if zoom is not None else (-1.0, 0.2) t0 = 1 / (2 * domain.delta_f) # merger at segment midpoint -> t = 0 phase_shift = np.exp(2j * np.pi * domain() * t0) - for row, (ax, ifo) in enumerate(zip(axes, ifos)): - for mode, color in zip(modes, colors): - td = [] - for wf in wf_fd[mode][ifo] * phase_shift: - times, x = one_sided_fd_to_td(wf, domain) - td.append(np.real(x)) - td = np.array(td) - ax.fill_between( - times - t0, - td.min(axis=0), - td.max(axis=0), - color=color, - alpha=0.3, - label=mode if row == 0 else None, - ) - # Maximum-probability waveform for this mode, as a solid line over the band. - if mode in map_fd: - m_times, m_x = one_sided_fd_to_td(map_fd[mode][ifo] * phase_shift, domain) - ax.plot(m_times - t0, np.real(m_x), color=color, lw=1.0, alpha=0.9) - if plot_data: - d_times, d_x = one_sided_fd_to_td(data_fd[ifo] * phase_shift, domain) - d_td = np.convolve(np.real(d_x), np.ones(4) / 4, mode="same") - ax.plot( - d_times - t0, - d_td, - color=data_color, - lw=1, - alpha=0.7, - zorder=0, - label="data" if row == 0 else None, - ) - ax.set_ylabel(f"{ifo}\nwhitened strain") - ax.set_xlim(*(zoom if zoom is not None else (-1.0, 0.2))) + # Grey whitened-data trace per detector (shared across that detector's panels). + data_td = {} + for ifo in ifos: + d_times, d_x = one_sided_fd_to_td(data_fd[ifo] * phase_shift, domain) + data_td[ifo] = (d_times - t0, np.convolve(np.real(d_x), np.ones(4) / 4, mode="same")) + + # One panel per (mode, ifo), stacked vertically (grouped by mode: dingo, then dingo-is). + panels = [(mode, ifo) for mode in modes for ifo in ifos] + fig, axes = plt.subplots( + len(panels), 1, figsize=(11, 2.6 * len(panels)), sharex=True, squeeze=False + ) + axes = axes[:, 0] + + for ax, (mode, ifo) in zip(axes, panels): + color = colors[mode] + td = np.array( + [np.real(one_sided_fd_to_td(wf, domain)[1]) for wf in wf_fd[mode][ifo] * phase_shift] + ) + times, _ = one_sided_fd_to_td(data_fd[ifo], domain) + + # Grey data (bottom) -> min/max band -> MAP line (top). + d_t, d_td = data_td[ifo] + ax.plot(d_t, d_td, color=_DATA_COLOR, lw=0.8, alpha=0.6, zorder=1) + ax.fill_between( + times - t0, td.min(axis=0), td.max(axis=0), color=color, alpha=0.3, zorder=2 + ) + if mode in map_fd: + _, m_x = one_sided_fd_to_td(map_fd[mode][ifo] * phase_shift, domain) + ax.plot(times - t0, np.real(m_x), color=color, lw=1.0, alpha=0.9, zorder=3) + + ax.set_xlim(*zoom) + ax.set_ylabel("whitened strain") + ax.text( + 0.01, + 0.95, + f"{mode} · {ifo}", + transform=ax.transAxes, + va="top", + ha="left", + fontweight="bold", + color=color, + ) axes[-1].set_xlabel("time to merger (s)") - axes[0].legend(loc="upper left") - if fig is not None: - fig.tight_layout() - fig.savefig(filename, dpi=200, bbox_inches="tight") - plt.close(fig) + fig.savefig(filename, dpi=200, bbox_inches="tight") + plt.close(fig) return axes diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index bc24c0c2a..2f56e7e2d 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -24,18 +24,17 @@ likelihood, detector projection, and ASD whitening — `core/` stays domain-agnostic. The module exposes: ```python -plot_ppd_td(wf_fd, data_fd, domain, filename="ppd_td.png", - zoom=None, axes=None, plot_data=True, colors=None) +plot_ppd_td(wf_fd, data_fd, domain, map_fd=None, filename="ppd_td.png", zoom=None) one_sided_fd_to_td(fd, domain) # IFFT helper, also here ``` `wf_fd` is `{mode: {ifo: (n_waveforms, n_freq) complex}}` (whitened model -waveforms) and `data_fd` is `{ifo: (n_freq,) complex}` (whitened data) — both -straight from `_compute_ppd`. `plot_ppd_td` draws **one coloured min/max envelope -per `mode`** (`"dingo"`, and `"dingo-is"` when present), labelled via -`_MODE_LABELS`, and overlays the grey whitened-data trace once. `axes=None` -creates a figure and saves to `filename`; passing `axes` draws onto them and -skips saving, for composition. +waveforms), `data_fd` is `{ifo: (n_freq,) complex}` (whitened data), and `map_fd` +is `{mode: {ifo: (n_freq,) complex}}` (the per-mode maximum-probability waveform) +— all from `_compute_ppd`. `plot_ppd_td` draws **one min/max envelope panel per +`(mode, ifo)`, stacked vertically** — `"dingo"` panels on top and, when present, +`"dingo-is"` below — each with the grey whitened-data trace beneath and the mode's +MAP waveform as a line on top. Each panel is labelled in-axes `"{mode} · {ifo}"`. The GW `Result` (subclass of `CoreResult`) keeps a thin convenience method that mirrors the `result.plot_corner` API and delegates to the module: @@ -45,9 +44,9 @@ result.plot_ppd_td(filename="ppd_td.png", credible_interval=0.9, num_waveforms=1000, num_processes=1, zoom=None, ppd=None) ``` -Like `plot_corner`, both the **Dingo** posterior (credible set from the raw -network samples) and — when the result is importance-sampled (`"weights"` -present) — the **Dingo-IS** posterior are overlaid on one figure. +Both the **Dingo** posterior (credible set from the raw network samples) and — +when the result is importance-sampled (`"weights"` present) — the **Dingo-IS** +posterior get their own stacked panels (Dingo on top, Dingo-IS below). The shared expensive computation stays a **private `Result` method** `_compute_ppd(self, credible_interval, num_waveforms, num_processes, seed)`, @@ -135,6 +134,20 @@ peak lands at t ≈ +0.02 s in both H1 and L1. `(-1.0, 0.2)`; `(-0.7, 0.1)` works well for GW230709). y-axis: whitened strain (dimensionless, σ units). +### Alternative considered: shade by density p(s | t) (rejected) + +We prototyped replacing the min/max band with the full per-time strain **density** +`p(s | t)` — the pushforward of the posterior through the waveform map, rendered as +an opacity/heatmap (mass-proportional draws, per-column-normalised + gamma, KDE +smoothing). Mathematically it is the "right" object (min/max is just its support), +and it looks great on the loud, well-localised GW150914. **But on the real +eccentric O4 events (e.g. S230709bi) it read as messy** — near merger the draws +decorrelate in phase, so the density fills a broad lens with fine internal banding +that is hard to parse, and finite-sample shot noise needs aggressive smoothing. +The min/max envelope + MAP line conveys the same "where does the signal live" story +far more legibly, so we kept it. (The density prototype was discarded, not merged; the +approach is recorded here in case it's worth revisiting with a cleaner rendering.) + Frequency-domain PPD plotting is deferred to a follow-up (a per-bin whitened FD view buries the signal below the unit-variance noise floor; the useful FD view is an amplitude-spectral-density / log-log plot, out of scope here). From 03175752ae5fa6f91796884d4d693e4f0760a4ad Mon Sep 17 00:00:00 2001 From: nihargupte-ph Date: Thu, 2 Jul 2026 16:34:14 +0200 Subject: [PATCH 10/11] feat(gw): nested credible-level bands in PPD plots credible_interval now accepts a list (e.g. [0.5, 0.9]) to overlay nested min/max bands. The HPD set at level L is {theta : density >= c_L}, so a tighter level is just a higher density threshold on the same draws: _compute_ppd returns the generated draws ordered by descending density plus a per-level count k, and plot_ppd_td shades the densest td[:k] per level (tighter = darker, with a CI legend). Co-Authored-By: Claude Opus 4.8 (1M context) --- dingo/gw/result.py | 91 ++++++++++++++----- dingo/gw/utils/plotting.py | 47 ++++++---- .../2026-06-30-strain-ppd-plotting-design.md | 5 + 3 files changed, 102 insertions(+), 41 deletions(-) diff --git a/dingo/gw/result.py b/dingo/gw/result.py index 68a53187a..e87a0a2fb 100644 --- a/dingo/gw/result.py +++ b/dingo/gw/result.py @@ -2,7 +2,7 @@ import time import warnings from functools import partial -from typing import Optional, Tuple +from typing import Optional, Sequence, Tuple, Union import numpy as np import yaml @@ -788,26 +788,28 @@ def pesummary_prior(self): def _compute_ppd( self, - credible_interval: float = 0.9, + credible_interval: Union[float, Sequence[float]] = 0.9, num_waveforms: int = 1000, num_processes: int = 1, seed: int = RANDOM_STATE, - ) -> Tuple[dict, dict, dict]: + ) -> Tuple[dict, dict, dict, dict]: """Generate whitened waveforms for the time-domain strain PPD. - Builds a ``credible_interval`` highest-posterior-density credible set for each - available posterior -- ``"dingo"`` (the raw network samples) and, when the result is - importance-sampled, ``"dingo-is"`` -- draws up to ``num_waveforms`` of each uniformly - at random (so the envelope spans the whole interval, not just the peak), and projects - them onto the detectors via the likelihood. Model waveforms come back already whitened - (``h / asd / noise_std``), matching the whitened data. + Builds the highest-posterior-density credible set at the largest requested + ``credible_interval`` level for each available posterior -- ``"dingo"`` (the raw + network samples) and, when the result is importance-sampled, ``"dingo-is"`` -- + generates up to ``num_waveforms`` of each, and projects them onto the detectors via + the likelihood. Model waveforms come back already whitened (``h / asd / noise_std``), + matching the whitened data. The generated draws are returned ordered by descending + density, and per-level counts are returned so the plot can draw nested bands. Parameters ---------- - credible_interval : float - Posterior probability of the credible set to draw from (default 0.9). + credible_interval : float or sequence of float + Posterior probability of the credible set(s). A sequence (e.g. ``[0.5, 0.9]``) + requests nested credible bands; the waveforms are generated for the largest level. num_waveforms : int - Maximum number of waveforms drawn from each credible set. + Maximum number of waveforms drawn from the (largest) credible set. num_processes : int Processes used for waveform generation. seed : int @@ -818,8 +820,8 @@ def _compute_ppd( wf_fd : dict ``{mode: {ifo: complex ndarray (n_kept, n_freq)}}`` whitened model waveforms, where ``mode`` is ``"dingo"`` (always) and ``"dingo-is"`` (present iff the result - is importance-sampled). Rows are the successfully generated draws; the rare failed - draw is dropped, so ``n_kept`` may be slightly below the number requested. + is importance-sampled). Rows are the successfully generated draws **ordered by + descending density**; the rare failed draw is dropped. data_fd : dict ``{ifo: complex ndarray (n_freq,)}`` whitened detector data (shared across modes). map_fd : dict @@ -827,8 +829,12 @@ def _compute_ppd( whitened waveform per mode -- max ``log_prob`` for ``"dingo"``, max ``log_likelihood + log_prior`` for ``"dingo-is"`` -- drawn as a line over the band. A mode is absent iff its MAP draw failed generation. + level_counts : dict + ``{mode: {level: k}}`` -- for each requested credible level, ``k`` is how many of + the (density-ordered) generated draws fall in that level's HPD set, so the plot can + take the densest ``wf_fd[mode][ifo][:k]`` as the band for that level. - All three are consumed by :func:`dingo.gw.utils.plotting.plot_ppd_td` (together with + All are consumed by :func:`dingo.gw.utils.plotting.plot_ppd_td` (together with ``self.domain``). Notes @@ -871,6 +877,12 @@ def _compute_ppd( ifos = list(self.context["waveform"].keys()) rng = np.random.default_rng(seed) + levels = sorted( + {float(credible_interval)} + if np.isscalar(credible_interval) + else {float(x) for x in credible_interval} + ) + # Whitened data, same convention as the (already-whitened) model waveforms. data_fd = { ifo: np.sqrt(4 * domain.delta_f) @@ -900,6 +912,7 @@ def _compute_ppd( # maximum-probability waveform per mode (drawn as a line over the band). wf_fd = {} map_fd = {} + level_counts = {} for mode in ("dingo", "dingo-is"): if mode == "dingo": mass = in_prior.astype(float) @@ -917,13 +930,27 @@ def _compute_ppd( mass = mass / mass.sum() rank = np.where(in_prior, rank, -np.inf) order = np.argsort(rank)[::-1] - n_keep = int(np.searchsorted(np.cumsum(mass[order]), credible_interval)) + 1 + cum_mass = np.cumsum(mass[order]) + + # The HPD set at level L is {theta : density(theta) >= c_L}. Find each level's + # density threshold c_L from the population (rank where the cumulative mass first + # reaches L); a smaller level is just a higher threshold on the same ranked draws, + # so nested bands are a top-k slice of the density-ordered draws -- no extra + # generation. Build (and generate) the set for the largest level. + thresholds = { + lvl: rank[order[min(int(np.searchsorted(cum_mass, lvl)), len(order) - 1)]] + for lvl in levels + } + n_keep = int(np.searchsorted(cum_mass, levels[-1])) + 1 credible_idx = order[:n_keep] - # Subsample uniformly from the credible set so the min/max envelope spans it. if len(credible_idx) > num_waveforms: credible_idx = rng.choice( credible_idx, size=num_waveforms, replace=False ) + # Order the draws by descending density so the plot can take the densest top-k as + # the inner (tighter) credible level. + credible_idx = credible_idx[np.argsort(rank[credible_idx])[::-1]] + sel_rank = rank[credible_idx] theta = self.samples.iloc[credible_idx] # The prior filter removes out-of-prior draws, but it does not guarantee every @@ -938,7 +965,9 @@ def _compute_ppd( theta, num_processes=num_processes, ) - n_failed = sum(s is None for s in signals) + failed = np.array([s is None for s in signals]) + n_failed = int(failed.sum()) + kept_rank = sel_rank[~failed] signals = [s for s in signals if s is not None] if not signals: raise RuntimeError( @@ -953,6 +982,10 @@ def _compute_ppd( wf_fd[mode] = { ifo: np.array([s["waveform"][ifo] for s in signals]) for ifo in ifos } + # Density-ordered draws already kept; count how many fall in each level's HPD set. + level_counts[mode] = { + lvl: int(np.sum(kept_rank >= thresholds[lvl])) for lvl in levels + } # Maximum-probability waveform for this mode: the single highest-ranked # (in-prior) draw -- max log_prob for "dingo", max (log_likelihood + log_prior) @@ -968,12 +1001,12 @@ def _compute_ppd( "Whitened waveform length does not match the inverse-FFT domain; " "time-domain PPD requires a uniform frequency domain." ) - return wf_fd, data_fd, map_fd + return wf_fd, data_fd, map_fd, level_counts def plot_ppd_td( self, filename: str = "ppd_td.png", - credible_interval: float = 0.9, + credible_interval: Union[float, Sequence[float]] = 0.9, num_waveforms: int = 1000, num_processes: int = 1, zoom: Optional[Tuple[float, float]] = None, @@ -986,10 +1019,12 @@ def plot_ppd_td( Panels are stacked vertically, one per (mode, detector): the **Dingo** posterior on top and, when the result is importance-sampled, the **Dingo-IS** posterior below. - ``zoom`` is the (left, right) x-limit in seconds-to-merger; defaults to (-1.0, 0.2). - Returns the ``(wf_fd, data_fd, map_fd)`` tuple that was plotted. + ``credible_interval`` may be a single level (default 0.9) or a sequence such as + ``[0.5, 0.9]`` to overlay nested credible bands. ``zoom`` is the (left, right) + x-limit in seconds-to-merger; defaults to (-1.0, 0.2). Returns the + ``(wf_fd, data_fd, map_fd)`` tuple that was plotted. """ - wf_fd, data_fd, map_fd = self._compute_ppd( + wf_fd, data_fd, map_fd, level_counts = self._compute_ppd( credible_interval, num_waveforms, num_processes ) domain = ( @@ -997,5 +1032,13 @@ def plot_ppd_td( if hasattr(self.domain, "base_domain") else self.domain ) - _plot_ppd_td(wf_fd, data_fd, domain, map_fd=map_fd, filename=filename, zoom=zoom) + _plot_ppd_td( + wf_fd, + data_fd, + domain, + map_fd=map_fd, + level_counts=level_counts, + filename=filename, + zoom=zoom, + ) return wf_fd, data_fd, map_fd diff --git a/dingo/gw/utils/plotting.py b/dingo/gw/utils/plotting.py index 343ff66a7..4535bbbb8 100644 --- a/dingo/gw/utils/plotting.py +++ b/dingo/gw/utils/plotting.py @@ -29,6 +29,7 @@ def plot_ppd_td( data_fd: dict, domain: Domain, map_fd: Optional[dict] = None, + level_counts: Optional[dict] = None, filename: str = "ppd_td.png", zoom: Optional[Tuple[float, float]] = None, ) -> np.ndarray: @@ -38,13 +39,14 @@ def plot_ppd_td( t = 0, segment-midpoint offset ``t0 = T/2``), shades the pointwise min/max envelope, and draws that mode's maximum-probability waveform as a line over it, all above the grey whitened-data trace. Panels are stacked vertically -- the **Dingo** posterior on top and, - when the result is importance-sampled, **Dingo-IS** below. + when the result is importance-sampled, **Dingo-IS** below. When ``level_counts`` gives more + than one credible level, nested bands are drawn (tighter level = darker). Parameters ---------- wf_fd : dict - ``{mode: {ifo: (n_waveforms, n_freq) complex}}`` whitened model waveforms, from - :meth:`Result._compute_ppd` (``mode`` is ``"dingo"`` and/or ``"dingo-is"``). + ``{mode: {ifo: (n_waveforms, n_freq) complex}}`` whitened model waveforms **ordered by + descending density**, from :meth:`Result._compute_ppd`. data_fd : dict ``{ifo: (n_freq,) complex}`` whitened detector data; its keys set the detectors. domain : Domain @@ -52,6 +54,9 @@ def plot_ppd_td( map_fd : dict or None ``{mode: {ifo: (n_freq,) complex}}`` maximum-probability waveform per mode, drawn as a solid line in the mode's colour. Modes absent from the dict (or ``None``) get no line. + level_counts : dict or None + ``{mode: {level: k}}`` -- draw a min/max band over the densest ``k`` draws for each + credible ``level``. ``None`` draws a single band over all draws. filename : str Output path for the saved figure. zoom : tuple or None @@ -83,35 +88,43 @@ def plot_ppd_td( ) axes = axes[:, 0] - for ax, (mode, ifo) in zip(axes, panels): + for row, (ax, (mode, ifo)) in enumerate(zip(axes, panels)): color = colors[mode] td = np.array( [np.real(one_sided_fd_to_td(wf, domain)[1]) for wf in wf_fd[mode][ifo] * phase_shift] ) times, _ = one_sided_fd_to_td(data_fd[ifo], domain) + tt = times - t0 - # Grey data (bottom) -> min/max band -> MAP line (top). + # Levels largest -> smallest so the tighter (darker) band is drawn on top. + counts = (level_counts or {}).get(mode, {}) or {None: len(td)} + ordered = sorted(counts.items(), key=lambda kv: -(kv[0] or 1.0)) + multi = len(ordered) > 1 + + # Grey data (bottom) -> nested min/max bands -> MAP line (top). d_t, d_td = data_td[ifo] ax.plot(d_t, d_td, color=_DATA_COLOR, lw=0.8, alpha=0.6, zorder=1) - ax.fill_between( - times - t0, td.min(axis=0), td.max(axis=0), color=color, alpha=0.3, zorder=2 - ) + for i, (lvl, k) in enumerate(ordered): + if k < 1: + continue + band = td[:k] + ax.fill_between( + tt, band.min(axis=0), band.max(axis=0), + color=color, alpha=0.25 + 0.22 * i, zorder=2 + i, + label=(f"{lvl:.0%} CI" if (multi and row == 0 and lvl is not None) else None), + ) if mode in map_fd: _, m_x = one_sided_fd_to_td(map_fd[mode][ifo] * phase_shift, domain) - ax.plot(times - t0, np.real(m_x), color=color, lw=1.0, alpha=0.9, zorder=3) + ax.plot(tt, np.real(m_x), color=color, lw=1.0, alpha=0.9, zorder=10) ax.set_xlim(*zoom) ax.set_ylabel("whitened strain") ax.text( - 0.01, - 0.95, - f"{mode} · {ifo}", - transform=ax.transAxes, - va="top", - ha="left", - fontweight="bold", - color=color, + 0.01, 0.95, f"{mode} · {ifo}", transform=ax.transAxes, + va="top", ha="left", fontweight="bold", color=color, ) + if multi and row == 0: + ax.legend(loc="upper right", fontsize=8, framealpha=0.9) axes[-1].set_xlabel("time to merger (s)") fig.savefig(filename, dpi=200, bbox_inches="tight") diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md index 2f56e7e2d..ceb18f5cf 100644 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md @@ -128,6 +128,11 @@ peak lands at t ≈ +0.02 s in both H1 and L1. - Each draw: phase-shift by `t0`, `one_sided_fd_to_td`, keep real part. - Band = pointwise **min/max envelope** across the selected draws (`td.min` / `td.max`), drawn with `fill_between`. +- **Nested credible levels:** `credible_interval` accepts a list (e.g. `[0.5, 0.9]`). + Since the HPD set at level L is `{θ : density ≥ c_L}`, a tighter level is just a + higher density threshold on the same ranked draws — so `_compute_ppd` returns the + draws ordered by descending density plus a per-level count `k`, and the plot shades + `td[:k].min/max` per level (tighter = darker). - Overlay: whitened data trace (`√(4Δf)·context_waveform/asd`, same shift + IFFT, light sliding-window average), grey. - x-axis: time to merger (s), linear, zoomed via `zoom=` (method default From be85e642769ef37c3417298a616fdcf9b05e1e96 Mon Sep 17 00:00:00 2001 From: Nihar Gupte Date: Tue, 14 Jul 2026 10:50:15 +0200 Subject: [PATCH 11/11] chore: drop superpowers design spec, gitignore docs/superpowers/ --- .gitignore | 2 +- .../2026-06-30-strain-ppd-plotting-design.md | 182 ------------------ 2 files changed, 1 insertion(+), 183 deletions(-) delete mode 100644 docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md diff --git a/.gitignore b/.gitignore index e76100e45..1a37a8ec3 100644 --- a/.gitignore +++ b/.gitignore @@ -154,4 +154,4 @@ dmypy.json dingo-venv/ # scratch/output dir for generated images and throwaway artifacts -local/ \ No newline at end of file +local/docs/superpowers/ diff --git a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md b/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md deleted file mode 100644 index ceb18f5cf..000000000 --- a/docs/superpowers/specs/2026-06-30-strain-ppd-plotting-design.md +++ /dev/null @@ -1,182 +0,0 @@ -# Strain PPD plotting — design - -**Date:** 2026-06-30 -**Branch:** `hackathon/ppd-plotting` - -## Goal - -Given a dingo `Result` HDF5 file, compute and plot the posterior-predictive -distribution (PPD) of the **whitened strain** in the **time domain**, per -detector, and save the figure as a PNG. (Frequency-domain plotting is out of -scope for this PR.) - -Derived from `/work/nihargupte/analysis/dingo/asimov/eccentricity/plotting/ -production_plots/strain_posterior.ipynb`, but stripped of the -eccentricity/glitch/prior-hull extras and generalized to a clean, reusable -script. - -## Deliverable - -Rendering lives in a new GW-specific module **`dingo/gw/utils/plotting.py`**, the -domain-specific counterpart to `dingo/core/utils/plotting.py` (which holds -`plot_corner_multi`). It cannot live in `core/` because PPDs need the GW -likelihood, detector projection, and ASD whitening — `core/` stays -domain-agnostic. The module exposes: - -```python -plot_ppd_td(wf_fd, data_fd, domain, map_fd=None, filename="ppd_td.png", zoom=None) -one_sided_fd_to_td(fd, domain) # IFFT helper, also here -``` - -`wf_fd` is `{mode: {ifo: (n_waveforms, n_freq) complex}}` (whitened model -waveforms), `data_fd` is `{ifo: (n_freq,) complex}` (whitened data), and `map_fd` -is `{mode: {ifo: (n_freq,) complex}}` (the per-mode maximum-probability waveform) -— all from `_compute_ppd`. `plot_ppd_td` draws **one min/max envelope panel per -`(mode, ifo)`, stacked vertically** — `"dingo"` panels on top and, when present, -`"dingo-is"` below — each with the grey whitened-data trace beneath and the mode's -MAP waveform as a line on top. Each panel is labelled in-axes `"{mode} · {ifo}"`. - -The GW `Result` (subclass of `CoreResult`) keeps a thin convenience method that -mirrors the `result.plot_corner` API and delegates to the module: - -```python -result.plot_ppd_td(filename="ppd_td.png", credible_interval=0.9, - num_waveforms=1000, num_processes=1, zoom=None, ppd=None) -``` - -Both the **Dingo** posterior (credible set from the raw network samples) and — -when the result is importance-sampled (`"weights"` present) — the **Dingo-IS** -posterior get their own stacked panels (Dingo on top, Dingo-IS below). - -The shared expensive computation stays a **private `Result` method** -`_compute_ppd(self, credible_interval, num_waveforms, num_processes, seed)`, -returning the tuple `(wf_fd, data_fd)`. It builds a credible set and generates -waveforms for **each** available posterior internally (no `weighting` argument): -`"dingo"` always, `"dingo-is"` iff importance-sampled. `plot_ppd_td` accepts a -precomputed `ppd=(wf_fd, data_fd)` to render without regenerating waveforms; the -`Result` method supplies `domain` from `self.domain`. - -Robustness (prior filter): the normalizing-flow proposal `q(theta)` is a smooth -density over a box and leaks a small fraction of draws outside the prior — with -finite `log_prob` but unphysical (e.g. tiny negative spin magnitudes). These -account for most waveform-generation failures. `_compute_ppd` therefore restricts -to the prior support up front (`self.prior.ln_prob(...)`, keeping `log_prior` -finite) — the same in-prior restriction `importance_sample` applies — before -generating anything. This removes those failures *and* keeps unphysical waveforms -out of the envelope, and it is a strict improvement for the `"dingo"` mode (which -otherwise gives those out-of-prior draws equal mass `1/N`). Measured on the -GW150914 fixture: 43/5000 draws (0.86%) are out-of-prior; all failures there were -among them, zero in-prior draws failed. - -Robustness (per-sample guard): the prior filter is **not** sufficient on its own. -A waveform model can fail on an *interior* parameter, and this bites even for -importance-sampled results whenever the waveform backend differs from the one used -at IS time — re-plotting a stored result with a newer/older `pyseobnr` / -`lalsimulation` is enough to make an IS-accepted draw fail. Observed on the -`S231001aq` (SEOBNRv5EHM) external result: 2/200 in-prior credible-set draws crash -inside pyseobnr's eccentric `CubicSpline` (`x must be strictly increasing`) under -`pyseobnr 0.3.6` / `lalsimulation 6.2.1`. So each draw is generated through the -module-level picklable `safe_signal` (`injection.py`): failures are dropped from -the `min`/`max` envelope with a warning giving the count, and only an all-failed -batch raises. (An earlier revision removed this guard on the mistaken assumption -that the prior filter left no in-prior failures for IS results; `S231001aq` is the -counterexample — the "IS already generated these" guarantee only holds for the -exact same waveform backend.) - -Dev fixture: GW150914 dingo-ci result (IMRPhenomXPHM, no `pyseobnr` needed), -`.../dingo-ci/outdir_GW150914/result/ -GW150914_data0_1126259462-4_importance_sampling_part1.hdf5`, rendered via the -throwaway `ppd_dev.py`. - -## Pipeline (per detector) - -1. `result = Result(file_name=...)`; `result._build_likelihood()`. - `domain = result.domain.base_domain if hasattr(result.domain, "base_domain") - else result.domain`. -2. **Sample selection.** Build the `credible_interval` (default 0.9) highest- - posterior-density credible set using the per-sample probability **mass** that - is correct for how the samples were drawn (no assumed reference density): for - `"dingo"` the raw samples are equal-mass draws from q(θ) — mass `1/N`, ranked - by `log_prob`; for `"dingo-is"` the target is the true posterior, so the mass - is the self-normalized importance weight (the `weights` column, ∝ p/q), ranked - by `log_likelihood + log_prior`. Sort by rank descending and keep the top - samples whose cumulative mass ≤ `credible_interval` (standard HPD estimator — - the same weighting `plot_corner` feeds to corner). Then draw up to - `num_waveforms` rows **uniformly at random** from that set (seeded - `np.random.default_rng`) so the envelope spans the whole interval rather than - only the peak (the notebook took the top-N highest-weight rows, biasing toward - the peak — we diverge here). -3. Generate whitened FD waveforms: - `apply_func_with_multiprocessing(result.likelihood.signal, samples_df, - num_processes)`. Each `wf["waveform"][ifo]` is **already whitened** - (`/asd/noise_std`), the same convention applied to the data, so model and - data share one y-scale. - -### Time offset (hidden from the user) - -The merger always lands at **t = 0**. Internally each waveform is phase-shifted -by a fixed offset derived from the data segment (`t0 = 1 / (2·domain.delta_f)`, -i.e. `T/2`, the segment **midpoint** where dingo centers the coalescence), via -`wf *= exp(2πi·f·t0)`, then the time axis has `t0` subtracted. `t0` is a private -detail — not a method argument. This generalizes the source notebook's hardcoded -`CENTRAL_TIME = 8` s (which was `T/2` for those `T = 16` s segments). Verified -empirically: for the GW150914 CI result (`T = 8` s, `t0 = 4`) the whitened-strain -peak lands at t ≈ +0.02 s in both H1 and L1. - -## Time domain (`_ppd_td.png`, one row per detector) - -- Each draw: phase-shift by `t0`, `one_sided_fd_to_td`, keep real part. -- Band = pointwise **min/max envelope** across the selected draws - (`td.min` / `td.max`), drawn with `fill_between`. -- **Nested credible levels:** `credible_interval` accepts a list (e.g. `[0.5, 0.9]`). - Since the HPD set at level L is `{θ : density ≥ c_L}`, a tighter level is just a - higher density threshold on the same ranked draws — so `_compute_ppd` returns the - draws ordered by descending density plus a per-level count `k`, and the plot shades - `td[:k].min/max` per level (tighter = darker). -- Overlay: whitened data trace (`√(4Δf)·context_waveform/asd`, same shift + - IFFT, light sliding-window average), grey. -- x-axis: time to merger (s), linear, zoomed via `zoom=` (method default - `(-1.0, 0.2)`; `(-0.7, 0.1)` works well for GW230709). y-axis: whitened strain - (dimensionless, σ units). - -### Alternative considered: shade by density p(s | t) (rejected) - -We prototyped replacing the min/max band with the full per-time strain **density** -`p(s | t)` — the pushforward of the posterior through the waveform map, rendered as -an opacity/heatmap (mass-proportional draws, per-column-normalised + gamma, KDE -smoothing). Mathematically it is the "right" object (min/max is just its support), -and it looks great on the loud, well-localised GW150914. **But on the real -eccentric O4 events (e.g. S230709bi) it read as messy** — near merger the draws -decorrelate in phase, so the density fills a broad lens with fine internal banding -that is hard to parse, and finite-sample shot noise needs aggressive smoothing. -The min/max envelope + MAP line conveys the same "where does the signal live" story -far more legibly, so we kept it. (The density prototype was discarded, not merged; the -approach is recorded here in case it's worth revisiting with a cleaner rendering.) - -Frequency-domain PPD plotting is deferred to a follow-up (a per-bin whitened FD -view buries the signal below the unit-variance noise floor; the useful FD view is -an amplitude-spectral-density / log-log plot, out of scope here). - -## Reused logic - -`one_sided_fd_to_td(fd, domain)` is reused verbatim from the notebook (a -module-level function in `dingo/gw/utils/plotting.py`): it builds the full Hermitian -two-sided spectrum (DC zeroed), `np.fft.ifft(...) * sqrt(N)` (normalization tied -to the whitening so noise has unit variance), `dt = 1/(2·f_max)`, output length -`2n-1`. Correct and whitening-consistent. - -## Testing / determinism - -- Seeded `np.random.default_rng` for resampling → deterministic figures. -- `tests/gw/test_ppd.py`: assert `one_sided_fd_to_td` (imported from - `dingo.gw.utils.plotting`) round-trips — a one-sided spectrum with a single populated - frequency bin returns a real cosine of the expected period (exercises the - Hermitian-mirror + normalization logic, the one piece of non-trivial new math). - -## Boundaries / conventions - -- Rendering in `gw/utils/plotting.py`, thin `Result.plot_ppd_td` + `_compute_ppd` in - `gw/result.py`; imports `gw/` + `core/` (allowed direction). No new classes, - no new dependencies. -- Whitening conventions for data (`√(4Δf)/asd`) and model (likelihood's - `/asd/noise_std`, `noise_std = 1/√(4Δf)`) are identical and must be preserved.