feat(gw): time-domain whitened-strain PPD plots on Result#393
feat(gw): time-domain whitened-strain PPD plots on Result#393nihargupte-ph wants to merge 11 commits into
Conversation
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
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) <noreply@anthropic.com>
lorenzopompili00
left a comment
There was a problem hiding this comment.
Thanks @nihargupte-ph! I left a few small comments.
Also, would it be possible to have an example of the figure that gets generated? It might also be nice to have this incorporated in the documentation.
| 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), |
There was a problem hiding this comment.
Just noting that this is a different statistic than what bilby (and most LVK reconstruction plots?) show, so I would avoid labelling it, e.g., "90% CI". The method plot_interferometer_waveform_posterior in bilby plots the pointwise 90% band: at each time sample, the central 90% of posterior waveforms. Here we plot the min/max envelope of draws from the 90% HPD set in parameter space, which is a different quantity. Both are reasonable things to plot, but with the same "90% CI" label someone comparing against a bilby plot of the same event might conclude the posteriors disagree when it's only the band definition. I would suggest labelling it "90% HPD" or "90% HPD envelope" (or similar), and maybe mentioning the difference with Bilby in the docstring.
| 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")) |
There was a problem hiding this comment.
Is this needed? Why not just plot the raw time domain data (as done by Bilby)?
| return axes | ||
|
|
||
|
|
||
| def one_sided_fd_to_td( |
There was a problem hiding this comment.
The docstring says the sqrt N normalization makes whitened noise have unit variance, but it comes out at sqrt 2, while and bilby's convention is exactly 1, so our y-axis is sqrt 2 larger than a bilby "whitened strain" plot of the same data. Since data and model go through the same function the plot shape is fine, it's just a scale issue.
A quick check from Claude:
import numpy as np
rng = np.random.default_rng(0)
n = 4097 # one-sided bins, f_max = (n-1)*df
fd = rng.normal(size=(500, n)) + 1j * rng.normal(size=(500, n))
# DINGO one_sided_fd_to_td
fd_d = fd.copy(); fd_d[:, 0] = 0
n_time = 2 * n - 1
td_dingo = np.fft.irfft(fd_d, n=n_time, axis=-1) * np.sqrt(n_time)
# Bilby get_whitened_time_series_from_whitened_frequency_series (full mask)
mask = np.ones(n, dtype=bool); mask[0] = False
fd_b = fd.copy(); fd_b[:, ~mask] = 0
td_bilby = np.fft.irfft(fd_b, axis=-1) * mask.sum()**0.5 / mask.mean()
print(td_dingo.std()) # 1.414
print(td_bilby.std()) # 1.000
print(td_dingo.std() / td_bilby.std()) # ≈ sqrt(2)More broadly, this makes me think if one could reuse more of the bilby functionality instead of reimplementing things. For example, there is a method get_whitened_time_series_from_whitened_frequency_series which already does what we need, and we already use the bilby Interferometer class so this wouldn't add extra dependencies.
| 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. |
There was a problem hiding this comment.
This docstring feels a bit too specific? I'm not sure there's any realistic case where using a different version of pyseobnr / lalsimulation would have waveform failing that used to generate tbh
| """Plot the time-domain whitened-strain posterior-predictive distribution. | ||
|
|
||
| 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 |
There was a problem hiding this comment.
Merger at t=0 reads very generic: what does "merger" mean exactly? Is it the trigger time, or the time at which the waveform peaks? Is it defined at the geocenter, or at the detector? The docstring (and label on the plot) could be more specific.
| 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 |
There was a problem hiding this comment.
Again, the comment seems overly specific
| 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 |
There was a problem hiding this comment.
Again, the comment seems overly specific
There was a problem hiding this comment.
Would it make sense to combine this with dingo.gw.gwutils?
Summary
Adds time-domain whitened-strain posterior predictive distribution (PPD) plotting to
dingo.gw.result.Result.Resultmethods to generate strain PPD plots (not a standalone module)dingo/gw/utils/plotting.py; injection support indingo/gw/injection.pydocs/specs/2026-06-30-strain-ppd-plotting-design.mdtests/gw/test_ppd.py.gitignores alocal/scratch dir for generated imagesTest plan
tests/gw/test_ppd.py🤖 Generated with Claude Code