Installation#

This guide explains how to install the required components to run the Makefile and draw the GWTC-5.0 distribution. It is recommended to use uv to create and manage the virtual environment, which avoids dependency conflicts.

Note

The model describes neutron stars and the primary black hole mass distribution as a broken power law between minimum and maximum masses, with two Gaussian peaks at \(\sim 30~M_\odot\) and \(\sim 9~M_\odot\) (posterior median; see Population Model (GWTC-5.0): FullPop for the full table). The file includes the posterior values for mass, spin, and merger-rate hyperparameters.

Using these hyperparameters, the model can be used to generate synthetic CBC distributions under the GWTC-5.0 FullPop population model. Here, we draw a distribution of one million samples to be used with the observing-scenarios pipeline together with ligo.skymap to simulate observing campaigns for upcoming runs.

Requirements

You will need the following:

gwpopulation

Python >= 3.11

Environment setup with uv
curl -LsSf https://astral.sh/uv/install.sh | sh
uv sync

Read the hyperparams file#

Below we provide a small script to read the FullPop (GWTC-5.0) result file with bilby and extract the MAP and posterior-median samples.

get_hyperparams

Show script
  1"""Build data/derived/hyperparams_map.csv and the joint mass-rate grid
  2caches from the GWTC-4.0/GWTC-5.0 source files. Missing sources are
  3fetched automatically (LIGO DCC for GWTC-4.0, Zenodo for GWTC-5.0). Run
  4directly to (re)build the cache; downstream code only reads the small
  5cached outputs, never the raw files."""
  6
  7import json
  8import os
  9import tarfile
 10import urllib.request
 11from pathlib import Path
 12
 13import numpy as np
 14import pandas as pd
 15from astropy.utils.data import download_file
 16from bilby.core.result import read_in_result
 17
 18#: Anchor paths on this file's own location, not the caller's cwd.
 19SCRIPT_DIR = Path(__file__).resolve().parent
 20REPO_ROOT = SCRIPT_DIR.parent.parent
 21DATA_DIR = REPO_ROOT / "data"
 22DERIVED_DIR = DATA_DIR / "derived"
 23
 24# --- Remote sources -----------------------------------------------------
 25ZENODO_RECORD_GWTC5 = "20292639"
 26#: Only available inside this archive, not as an individual download.
 27ZENODO_ARCHIVE_GWTC5 = "popsummary_files.tar.gz"
 28ZENODO_FILENAME_GWTC5 = (
 29    "production_1_mass_NotchFilterBinnedPairingMassDistribution_"
 30    "redshift_powerlaw_mag_iid_spin_magnitude_gaussian_tilt_"
 31    "iid_spin_orientation_popsummary.h5"
 32)
 33
 34# GWTC-4.0: AllCBC_FullPop.h5 (from Zenodo, inside analyses_AllCBC.tar) is
 35# the source actually used to build the published rate table. Its
 36# per-category rate_bns/rate_nsbh/rate_bbh/rate_full columns match the
 37# paper's Table 2 at the MAP row (rate_full = 130.28 vs 130 published).
 38# The DCC bilby Result file below has the same hyperparameter posterior
 39# (log_likelihood matches row-for-row) but a *different*, non-reproducible
 40# `rate` column. compute_rate_posterior() in gwpopulation_pipe draws it
 41# from a random gamma distribution, so two independent post-processing
 42# runs of the same posterior don't agree, and it never had rate_bns/
 43# rate_nsbh/rate_bbh/rate_full at all. Kept only as a fallback/for tests.
 44DCC_BASE_URL = "https://dcc.ligo.org/LIGO-T2500311/public"
 45DCC_FILENAME_GWTC4 = (
 46    "baseline5_widesigmachi2_mass_NotchFilterBinnedPairingMassDistribution_"
 47    "redshift_powerlaw_mag_iid_spin_magnitude_gaussian_tilt_"
 48    "iid_spin_orientation_result.hdf5"
 49)
 50
 51ZENODO_RECORD_GWTC4 = "16911563"
 52ZENODO_ARCHIVE_GWTC4 = "analyses_AllCBC.tar"
 53ZENODO_FILENAME_GWTC4 = "AllCBC_FullPop.h5"
 54
 55#: PixelPop's popsummary file, from the same GWTC-5.0 archive as FullPop's.
 56#: Unlike FullPop, it carries no scalar `rate` column. The rate lives in
 57#: the `joint_pixelpop_rate` grid and has to be integrated (see
 58#: compute_rate_summary_pixelpop).
 59ZENODO_MEMBER_PIXELPOP = "all_cbc_varcut1_popsummary.h5"
 60ZENODO_FILENAME_PIXELPOP = "pixelpop_popsummary.h5"
 61
 62
 63def _verify_size(path: Path, expected_size: int, source_desc: str) -> None:
 64    """Raise OSError (and delete the file) if its size doesn't match
 65    expected_size. Catches truncated/incomplete downloads."""
 66    actual_size = path.stat().st_size
 67    if actual_size != expected_size:
 68        path.unlink()
 69        raise OSError(
 70            f"Incomplete download of {source_desc}: expected {expected_size} "
 71            f"bytes, got {actual_size}. Deleted the truncated file at {path}, "
 72            f"re-run to retry."
 73        )
 74
 75
 76def download_from_zenodo(record_id: str, filename: str, dest_path: Path) -> Path:
 77    """
 78    Download `filename` from a public Zenodo record via the official API
 79    (queries the record's file list rather than guessing a download URL),
 80    unless it already exists at `dest_path`.
 81    """
 82    dest_path = Path(dest_path)
 83    if dest_path.exists():
 84        return dest_path
 85
 86    api_url = f"https://zenodo.org/api/records/{record_id}"
 87    with urllib.request.urlopen(api_url, timeout=30) as response:
 88        record = json.load(response)
 89
 90    matches = [f for f in record["files"] if f["key"] == filename]
 91    if not matches:
 92        available = [f["key"] for f in record["files"]]
 93        raise FileNotFoundError(
 94            f"{filename!r} not found in Zenodo record {record_id}. "
 95            f"Available files: {available}"
 96        )
 97
 98    download_url = matches[0]["links"]["self"]
 99    expected_size = matches[0]["size"]
100    print(f"Downloading {filename} from Zenodo record {record_id} ...")
101    cached_path = Path(download_file(download_url, cache=True, show_progress=True))
102    _verify_size(cached_path, expected_size, f"{filename!r} (Zenodo {record_id})")
103
104    dest_path.parent.mkdir(parents=True, exist_ok=True)
105    cached_path.replace(dest_path)
106    print(f"Saved to {dest_path}")
107    return dest_path
108
109
110def download_and_extract_from_zenodo_tarball(
111    record_id: str, archive_name: str, member_filename: str, dest_path: Path
112) -> Path:
113    """Download `archive_name` (a .tar or .tar.gz, compression is
114    auto-detected) from a public Zenodo record and extract only the member
115    ending in `member_filename` to `dest_path`. No-op if `dest_path` already
116    exists. The archive itself can be several GB, so this is meant as a
117    one-time step, not a routine call."""
118    dest_path = Path(dest_path)
119    if dest_path.exists():
120        return dest_path
121
122    api_url = f"https://zenodo.org/api/records/{record_id}"
123    with urllib.request.urlopen(api_url, timeout=30) as response:
124        record = json.load(response)
125
126    matches = [f for f in record["files"] if f["key"] == archive_name]
127    if not matches:
128        available = [f["key"] for f in record["files"]]
129        raise FileNotFoundError(
130            f"{archive_name!r} not found in Zenodo record {record_id}. "
131            f"Available files: {available}"
132        )
133
134    download_url = matches[0]["links"]["self"]
135    expected_size = matches[0]["size"]
136    size_gb = expected_size / 1e9
137    print(
138        f"Downloading {archive_name} ({size_gb:.1f} GB) from Zenodo record {record_id} ..."
139    )
140    cached_path = Path(download_file(download_url, cache=True, show_progress=True))
141    _verify_size(cached_path, expected_size, f"{archive_name!r} (Zenodo {record_id})")
142
143    print(f"Looking for a member ending in {member_filename!r} ...")
144    with tarfile.open(cached_path) as tar:
145        member_matches = [
146            m for m in tar.getmembers() if m.name.endswith(member_filename)
147        ]
148        if not member_matches:
149            raise FileNotFoundError(
150                f"No member ending in {member_filename!r} inside {archive_name}"
151            )
152        member = member_matches[0]
153        print(f"Extracting {member.name} ...")
154        tar.extract(member, path=dest_path.parent, filter="data")
155        extracted_path = dest_path.parent / member.name
156
157    dest_path.parent.mkdir(parents=True, exist_ok=True)
158    extracted_path.replace(dest_path)
159    print(f"Saved to {dest_path}")
160    return dest_path
161
162
163def download_from_dcc(filename: str, dest_path: Path) -> Path:
164    """
165    Download `filename` from the public LIGO DCC directory for T2500311,
166    unless it already exists at `dest_path`.
167    """
168    dest_path = Path(dest_path)
169    if dest_path.exists():
170        return dest_path
171
172    url = f"{DCC_BASE_URL}/{filename}"
173    head_request = urllib.request.Request(url, method="HEAD")
174    with urllib.request.urlopen(head_request, timeout=30) as response:
175        expected_size = int(response.headers["Content-Length"])
176
177    print(f"Downloading {filename} from LIGO DCC (T2500311) ...")
178    cached_path = Path(download_file(url, cache=True, show_progress=True))
179    _verify_size(cached_path, expected_size, f"{filename!r} (LIGO DCC T2500311)")
180
181    dest_path.parent.mkdir(parents=True, exist_ok=True)
182    cached_path.replace(dest_path)
183    print(f"Saved to {dest_path}")
184    return dest_path
185
186
187def to_rst(df, title="Hyperparameters of the FullPop-4.0 model"):
188    widths = [max(len(str(x)) for x in df[col]) for col in df.columns]
189    widths = [max(w, len(col)) for w, col in zip(widths, df.columns)]
190
191    def hline(sep="-"):
192        return "+" + "+".join(sep * (w + 2) for w in widths) + "+"
193
194    def row(cells):
195        return (
196            "|" + "|".join(f" {str(c).ljust(w)} " for c, w in zip(cells, widths)) + "|"
197        )
198
199    # build the grid table body (no directive yet)
200    body = [hline("-"), row(df.columns), hline("=")]
201    for _, r in df.iterrows():
202        body.append(row(r))
203        body.append(hline())
204
205    # indent body so it belongs to the table directive
206    indented = ["   " + line for line in body]  # 3 spaces is conventional
207
208    lines = [f".. table:: {title}", ""]
209    lines.extend(indented)
210    return "\n".join(lines)
211
212
213# --- in  LaTeX ---
214def to_latex(df, caption="Hyperparameters of the FullPop-4.0 model"):
215    out = []
216    out.append(r"\begin{table}[ht]")
217    out.append(r"\centering")
218    out.append(rf"\caption{{{caption}}}")
219    out.append(r"\begin{tabular}{lll}")
220    out.append(r"\hline")
221    out.append(r"Parameter & Description & Value\\")
222    out.append(r"\hline")
223    for _, row in df.iterrows():
224        out.append(rf"{row['Parameter']} & {row['Description']} & {row['Value']} \\")
225    out.append(r"\hline")
226    out.append(r"\end{tabular}")
227    out.append(r"\end{table}")
228    return "\n".join(out)
229
230
231# --- paper-style LaTeX table (booktabs, sectioned) -----------------------
232def to_latex_paper(
233    row: pd.Series,
234    caption: str = "Hyperparameters of the FullPop-4.0 Distribution Model",
235    label: str = "tab:hyperparams",
236) -> str:
237    """Build the sectioned booktabs-style hyperparameter table used in the
238    paper (Mass Distribution / Pairing Function / Spin Distribution), from
239    one row of `row` (a MAP hyperparameter Series, e.g.
240    hyperparams_map.loc["GWTC-5.0 (popsummary)"])."""
241
242    def v(key, decimals):
243        return f"{row[key]:.{decimals}f}"
244
245    return "\n".join(
246        [
247            r"% !TEX root = ../main.tex",
248            r"\begin{table*}[htb]",
249            r"\renewcommand\arraystretch{1.08}",
250            r"\setlength{\tabcolsep}{16pt}",
251            r"\centering",
252            rf"\caption{{{caption}}}",
253            rf"\label{{{label}}}",
254            r"\begin{tabular}{llr}",
255            r"\toprule",
256            r"\textbf{Parameter} & \textbf{Description} & \textbf{Posterior Value} \\",
257            r"\midrule",
258            r"\multicolumn{3}{c}{\textit{Mass Distribution}} \\",
259            r"\midrule",
260            rf"$m_{{\mathrm{{min,NS}}}}$ & Minimum neutron star mass ($M_\odot$) & ${v('NSmin', 2)}$ \\",
261            rf"$m_{{\mathrm{{max,NS}}}} \equiv \gamma_{{\mathrm{{low}},1}}$ & Maximum neutron star mass ($M_\odot$) & ${v('NSmax', 2)}$ \\",
262            rf"$m_{{\mathrm{{min,BH}}}} \equiv \gamma_{{\mathrm{{high}},1}}$ & Minimum black hole mass ($M_\odot$) & ${v('BHmin', 2)}$ \\",
263            rf"$\gamma_{{\mathrm{{low}},2}}$  & Lower boundary of pair-instability gap ($M_\odot$) & ${v('UPPERmin', 2)}$ \\",
264            rf"$\gamma_{{\mathrm{{high}},2}}$ & Upper boundary of pair-instability gap ($M_\odot$) & ${v('UPPERmax', 2)}$ \\",
265            rf"$m_{{\mathrm{{max,BH}}}}$ & Maximum black hole mass ($M_\odot$) & ${v('BHmax', 1)}$ \\",
266            r"\midrule",
267            rf"$\alpha_1$ & Power-law exponent for masses below $m_{{\mathrm{{max,NS}}}}$ & ${v('alpha_1', 2)}$ \\",
268            rf"$\alpha_{{\mathrm{{dip}}}}$ & Power-law exponent within the NS–BH mass gap & ${v('alpha_dip', 2)}$ \\",
269            rf"$\alpha_2$ & Power-law exponent for masses above $m_{{\mathrm{{min,BH}}}}$ & ${v('alpha_2', 2)}$ \\",
270            r"\midrule",
271            rf"$\mu_{{\mathrm{{peak}},1}}$ & Mean of primary Gaussian peak ($M_\odot$) & ${v('mu1', 2)}$ \\",
272            rf"$\sigma_{{\mathrm{{peak}},1}}$ & Std. dev. of primary Gaussian peak ($M_\odot$) & ${v('sig1', 2)}$ \\",
273            rf"$c_1$ & Mixing fraction of primary Gaussian peak & ${v('mix1', 2)}$ \\",
274            rf"$\mu_{{\mathrm{{peak}},2}}$ & Mean of secondary Gaussian peak ($M_\odot$) & ${v('mu2', 2)}$ \\",
275            rf"$\sigma_{{\mathrm{{peak}},2}}$ & Std. dev. of secondary Gaussian peak ($M_\odot$) & ${v('sig2', 2)}$ \\",
276            rf"$c_2$ & Mixing fraction of secondary Gaussian peak & ${v('mix2', 2)}$ \\",
277            r"\midrule",
278            rf"$A_1$ & Depth of primary mass gap suppression & ${v('A', 3)}$ \\",
279            rf"$A_2$ & Depth of pair-instability gap suppression & ${v('A2', 3)}$ \\",
280            "",
281            r"\midrule",
282            rf"$\eta_0$ & Sharpness of low-mass truncation & ${v('n0', 0)}$ \\",
283            rf"$\eta_1$ & Sharpness at $m_{{\mathrm{{max,NS}}}}$ & ${v('n1', 0)}$ \\",
284            rf"$\eta_2$ & Sharpness at $m_{{\mathrm{{min,BH}}}}$ & ${v('n2', 0)}$ \\",
285            rf"$\eta_3$ & Sharpness at $\gamma_{{\mathrm{{low}},2}}$ & ${v('n3', 0)}$ \\",
286            rf"$\eta_4$ & Sharpness at $\gamma_{{\mathrm{{high}},2}}$ & ${v('n4', 0)}$ \\",
287            rf"$\eta_5$ & Sharpness of high-mass truncation & ${v('n5', 2)}$ \\",
288            "",
289            r"\midrule",
290            r"\multicolumn{3}{c}{\textit{Pairing Function}} \\",
291            r"\midrule",
292            rf"$m_{{\mathrm{{break}}}}$ & Pairing function break mass ($M_\odot$) & ${v('mbreak', 1)}$ \\",
293            rf"$\beta_1$ & Pairing power-law index for $m_2 < m_{{\mathrm{{break}}}}$ & ${v('beta_pair_1', 2)}$ \\",
294            rf"$\beta_2$ & Pairing power-law index for $m_2 \geq m_{{\mathrm{{break}}}}$ & ${v('beta_pair_2', 2)}$ \\",
295            "",
296            r"\midrule",
297            r"\multicolumn{3}{c}{\textit{Spin Distribution}} \\",
298            r"\midrule",
299            rf"$\mu_{{\chi}}$ & Mean of spin magnitude Gaussian component & ${v('mu_chi', 3)}$ \\",
300            rf"$\sigma_{{\chi}}$ & Std. dev. of spin magnitude Gaussian component & ${v('sigma_chi', 2)}$ \\",
301            rf"$a_{{\mathrm{{max}}}}$ & Maximum spin magnitude & ${v('amax', 0)}$ \\",
302            r"\midrule",
303            rf"$\xi_{{\mathrm{{spin}}}}$ & Fraction of BHs in preferentially aligned component & ${v('xi_spin', 2)}$ \\",
304            rf"$\sigma_{{\mathrm{{spin}}}}$ & Width of preferentially aligned component & ${v('sigma_spin', 3)}$ \\",
305            r"\bottomrule",
306            r"\end{tabular}",
307            r"\end{table*}",
308            "",
309        ]
310    )
311
312
313def _get_map_sample(hyperparams) -> pd.Series:
314    """
315    Return MAP sample if prior is informative, else ML sample.
316    """
317    post = hyperparams.posterior.copy()
318    if "log_prior" in post and post["log_prior"].nunique() > 1:
319        score = post.log_likelihood + post.log_prior
320        return post.iloc[np.argmax(score)]
321    else:
322        return post.iloc[np.argmax(post.log_likelihood)]
323
324
325def load_hyperparams_map(path: str) -> pd.Series:
326    """Load the MAP hyperparameter sample from a FullPop result file,
327    dispatching on format: popsummary (via PopulationResult) or bilby
328    Result. Only the hyperparameter table is loaded, not the rate grids."""
329    try:
330        from popsummary.popresult import PopulationResult
331
332        popresult = PopulationResult(path)
333        df = pd.DataFrame(
334            popresult.get_hyperparameter_samples(),
335            columns=popresult.get_metadata("hyperparameters"),
336        )
337        return df.iloc[(df.log_likelihood + df.log_prior).idxmax()]
338    except (KeyError, OSError):
339        hyperparams = read_in_result(path)
340        return _get_map_sample(hyperparams)
341
342
343def extract_hyperparams_map(
344    sources: dict[str, str], cache_csv: str, force: bool = False
345) -> pd.DataFrame:
346    """
347    Return a DataFrame of MAP hyperparameters, one row per label in
348    `sources` (label -> file path). Labels already present in `cache_csv`
349    are reused instead of reopening their (possibly multi-GB) source file;
350    pass force=True to recompute everything.
351    """
352    cached = pd.DataFrame()
353    if not force and os.path.exists(cache_csv):
354        cached = pd.read_csv(cache_csv, index_col=0)
355        sources = {
356            label: path for label, path in sources.items() if label not in cached.index
357        }
358
359    if not sources:
360        return cached
361
362    rows = {label: load_hyperparams_map(path) for label, path in sources.items()}
363    new_df = pd.DataFrame(rows).T
364    df = pd.concat([cached, new_df]) if not cached.empty else new_df
365    df.index.name = "label"
366    os.makedirs(os.path.dirname(cache_csv) or ".", exist_ok=True)
367    df.to_csv(cache_csv)
368    return df
369
370
371def load_hyperparams_median(path: str) -> pd.Series:
372    """Load the posterior median of every hyperparameter from a FullPop
373    result file, dispatching on format like load_hyperparams_map. Matches
374    GWTC-5.0's own convention of reporting posterior medians rather than
375    a single MAP sample (see compute_rate_summary_fullpop's docstring)."""
376    try:
377        from popsummary.popresult import PopulationResult
378
379        popresult = PopulationResult(path)
380        df = pd.DataFrame(
381            popresult.get_hyperparameter_samples(),
382            columns=popresult.get_metadata("hyperparameters"),
383        )
384        return df.median(numeric_only=True)
385    except (KeyError, OSError):
386        hyperparams = read_in_result(path)
387        return hyperparams.posterior.median(numeric_only=True)
388
389
390def extract_hyperparams_median(
391    sources: dict[str, str], cache_csv: str, force: bool = False
392) -> pd.DataFrame:
393    """
394    Return a DataFrame of posterior-median hyperparameters, one row per
395    label in `sources` (label -> file path). Same caching pattern as
396    extract_hyperparams_map: labels already present in `cache_csv` are
397    reused instead of reopening their (possibly multi-GB) source file;
398    pass force=True to recompute everything.
399    """
400    cached = pd.DataFrame()
401    if not force and os.path.exists(cache_csv):
402        cached = pd.read_csv(cache_csv, index_col=0)
403        sources = {
404            label: path for label, path in sources.items() if label not in cached.index
405        }
406
407    if not sources:
408        return cached
409
410    rows = {label: load_hyperparams_median(path) for label, path in sources.items()}
411    new_df = pd.DataFrame(rows).T
412    df = pd.concat([cached, new_df]) if not cached.empty else new_df
413    df.index.name = "label"
414    os.makedirs(os.path.dirname(cache_csv) or ".", exist_ok=True)
415    df.to_csv(cache_csv)
416    return df
417
418
419#: The three small joint-mass rate grids present in both catalogues'
420#: popsummary files (600x600 each). Excludes
421#: primary_mass_secondary_mass_joint_full_posterior, which only GWTC-5.0
422#: carries and is ~7 GB once loaded.
423JOINT_GRID_NAMES = (
424    "primary_mass_secondary_mass_joint_median",
425    "primary_mass_secondary_mass_joint_ppd",
426    "primary_mass_secondary_mass_joint_uncertainty",
427)
428
429
430def load_joint_grids(
431    path: str, grid_names: tuple[str, ...] = JOINT_GRID_NAMES
432) -> dict[str, tuple]:
433    """
434    Load the small 2D joint-mass rate grids from a popsummary file.
435
436    Returns {grid_name: (m1_positions, m2_positions, rates)}. Grid names
437    that don't exist in this particular file are silently skipped (e.g. a
438    future catalogue might not ship all three). Only ever touches the
439    named grids, never primary_mass_secondary_mass_joint_full_posterior.
440    """
441    from popsummary.popresult import PopulationResult
442
443    popresult = PopulationResult(path)
444    grids = {}
445    for name in grid_names:
446        try:
447            (m1_pos, m2_pos), rates = popresult.get_rates_on_grids(name)
448        except KeyError:
449            continue
450        grids[name] = (
451            np.asarray(m1_pos).ravel(),
452            np.asarray(m2_pos).ravel(),
453            np.asarray(rates),
454        )
455    return grids
456
457
458def extract_joint_grids(
459    sources: dict[str, str],
460    cache_dir: str,
461    grid_names: tuple[str, ...] = JOINT_GRID_NAMES,
462    force: bool = False,
463) -> dict[str, str]:
464    """
465    Save the small joint-mass grids for each label in `sources` (label ->
466    popsummary file path) to one compressed .npz per label under
467    `cache_dir`. Skips labels whose .npz already exists unless force=True.
468    Returns {label: npz_path}.
469    """
470    os.makedirs(cache_dir, exist_ok=True)
471    safe = str.maketrans({c: "_" for c in " ()."})
472    out_paths = {}
473    for label, path in sources.items():
474        npz_path = os.path.join(cache_dir, f"{label.translate(safe)}_grids.npz")
475        out_paths[label] = npz_path
476        if not force and os.path.exists(npz_path):
477            continue
478        grids = load_joint_grids(path, grid_names)
479        arrays = {}
480        for name, (m1_pos, m2_pos, rates) in grids.items():
481            arrays[f"{name}__m1"] = m1_pos
482            arrays[f"{name}__m2"] = m2_pos
483            arrays[f"{name}__rates"] = rates
484        np.savez_compressed(npz_path, **arrays)
485    return out_paths
486
487
488def compute_rate_summary_fullpop(path: str, rate_column: str = "rate") -> dict:
489    """Summarise a FullPop-style scalar rate column, giving the median and
490    5%/95% credible interval (the convention GWTC-5.0 itself uses, since
491    its published rates come from `np.median` over the posterior, not a
492    single MAP sample, see compute_rate_summary_pixelpop for why
493    PixelPop has no MAP/ML point at all here), plus the MAP point for
494    reference (matches GWTC-4.0's older, MAP-based convention)."""
495    from popsummary.popresult import PopulationResult
496
497    popresult = PopulationResult(path)
498    df = pd.DataFrame(
499        popresult.get_hyperparameter_samples(),
500        columns=popresult.get_metadata("hyperparameters"),
501    )
502    rate = df[rate_column]
503    lower_5, median, upper_95 = np.percentile(rate, [5, 50, 95])
504    map_row = df.iloc[(df.log_likelihood + df.log_prior).idxmax()]
505    return {
506        "map": map_row[rate_column],
507        "median": median,
508        "lower_5": lower_5,
509        "upper_95": upper_95,
510    }
511
512
513def compute_rate_summary_pixelpop(
514    path: str, grid_name: str = "joint_pixelpop_rate"
515) -> dict:
516    """Summarise PixelPop's total rate by integrating its per-sample 2D
517    joint mass-rate grid over the whole domain (log_mass_1 x log_mass_2).
518
519    PixelPop's popsummary file has no scalar `rate` column at all. The
520    rate only exists as a density on this grid. The grid's `positions`
521    array (n_edges x n_edges) is one larger per axis than its `rates`
522    array (n_bins x n_bins = (n_edges-1) x (n_edges-1)). `positions` are
523    bin edges, `rates` are the density at bin centres. Reshaping
524    `rates` as (n_samples, n_bins, n_bins) and summing weighted by the
525    (uniform) bin widths reproduces the file's own 1D marginals to
526    within 0.01%, verified before trusting this integration.
527
528    Also unlike FullPop, there is no `log_prior` column here (PixelPop's
529    binned-GP smoothing hyperparameters aren't sampled with an
530    informative prior in the same sense), so there is no true MAP, only
531    a maximum-likelihood (ML) point, returned as `best_fit_ml`.
532    """
533    import h5py
534
535    with h5py.File(path, "r") as f:
536        g = f["posterior/rates_on_grids"]
537        edges = np.unique(g["log_mass_1/positions"][:])
538        n_bins = len(edges) - 1
539        d_logm = np.diff(edges).mean()
540
541        rates = g[f"{grid_name}/rates"][:]
542        n_samples = rates.shape[0]
543        rates_3d = rates.reshape(n_samples, n_bins, n_bins)
544        total_rate = rates_3d.sum(axis=(1, 2)) * d_logm * d_logm
545
546        names = list(f.attrs["hyperparameters"])
547        idx = {n: i for i, n in enumerate(names)}
548        log_likelihood = f["posterior/hyperparameter_samples"][:, idx["log_likelihood"]]
549
550    best_i = np.argmax(log_likelihood)
551    lower_5, median, upper_95 = np.percentile(total_rate, [5, 50, 95])
552    return {
553        "best_fit_ml": total_rate[best_i],
554        "median": median,
555        "lower_5": lower_5,
556        "upper_95": upper_95,
557    }
558
559
560#: label -> (path, kind, rate_column). kind picks which compute_rate_summary_*
561#: function to use; rate_column is only meaningful for kind="fullpop".
562RATE_SOURCES = {
563    "GWTC-4.0 FullPop": (
564        str(DATA_DIR / "raw" / ZENODO_FILENAME_GWTC4),
565        "fullpop",
566        "rate_full",
567    ),
568    "GWTC-5.0 FullPop": (
569        str(DATA_DIR / "raw" / ZENODO_FILENAME_GWTC5),
570        "fullpop",
571        "rate",
572    ),
573    "GWTC-5.0 PixelPop": (
574        str(DATA_DIR / "raw" / ZENODO_FILENAME_PIXELPOP),
575        "pixelpop",
576        None,
577    ),
578}
579
580
581#: Per-class merger rates as published in GWTC-5.0's results paper
582#: (https://arxiv.org/abs/2605.27226), Table 2, FullPop and PixelPop
583#: rows. Static citations, not derived from any local file, kept here so
584#: they exist in exactly one place instead of being copy-pasted into every
585#: script that needs them (population_stats.py, detection_rate.ipynb, ...).
586PUBLISHED_RATES_TABLE2 = [
587    {
588        "catalog": "GWTC-5.0 FullPop",
589        "population": "BNS",
590        "lower": 15.4,
591        "mid": 59.3,
592        "upper": 154.7,
593    },
594    {
595        "catalog": "GWTC-5.0 FullPop",
596        "population": "NSBH",
597        "lower": 6.7,
598        "mid": 14.2,
599        "upper": 26.2,
600    },
601    {
602        "catalog": "GWTC-5.0 FullPop",
603        "population": "BBH",
604        "lower": 27.5,
605        "mid": 36.0,
606        "upper": 47.1,
607    },
608    {
609        "catalog": "GWTC-5.0 PixelPop",
610        "population": "BNS",
611        "lower": 5.2,
612        "mid": 23.4,
613        "upper": 78.1,
614    },
615    {
616        "catalog": "GWTC-5.0 PixelPop",
617        "population": "NSBH",
618        "lower": 7.0,
619        "mid": 15.9,
620        "upper": 32.8,
621    },
622    {
623        "catalog": "GWTC-5.0 PixelPop",
624        "population": "BBH",
625        "lower": 28.4,
626        "mid": 37.5,
627        "upper": 49.4,
628    },
629]
630
631
632def write_published_rates_table2(cache_csv: str) -> pd.DataFrame:
633    """Write PUBLISHED_RATES_TABLE2 to `cache_csv`. Always overwrites,
634    since these are static citations, not the output of a slow
635    computation, so there's no cache to preserve across runs."""
636    df = pd.DataFrame(PUBLISHED_RATES_TABLE2)
637    os.makedirs(os.path.dirname(cache_csv) or ".", exist_ok=True)
638    df.to_csv(cache_csv, index=False)
639    return df
640
641
642def extract_rate_summary(
643    rate_sources: dict[str, tuple[str, str, str | None]],
644    cache_csv: str,
645    force: bool = False,
646) -> pd.DataFrame:
647    """Return a DataFrame of rate summaries, one row per label in
648    `rate_sources`. Labels already present in `cache_csv` are reused
649    unless force=True."""
650    cached = pd.DataFrame()
651    if not force and os.path.exists(cache_csv):
652        cached = pd.read_csv(cache_csv, index_col=0)
653        rate_sources = {
654            label: v for label, v in rate_sources.items() if label not in cached.index
655        }
656
657    if not rate_sources:
658        return cached
659
660    rows = {}
661    for label, (path, kind, rate_column) in rate_sources.items():
662        if kind == "fullpop":
663            rows[label] = compute_rate_summary_fullpop(path, rate_column)
664        elif kind == "pixelpop":
665            rows[label] = compute_rate_summary_pixelpop(path)
666        else:
667            raise ValueError(f"unknown kind {kind!r} for label {label!r}")
668
669    new_df = pd.DataFrame(rows).T
670    df = pd.concat([cached, new_df]) if not cached.empty else new_df
671    df.index.name = "label"
672    os.makedirs(os.path.dirname(cache_csv) or ".", exist_ok=True)
673    df.to_csv(cache_csv)
674    return df
675
676
677PARAMS_INFO = {
678    "alpha_1": (
679        r":math:`\alpha_1`",
680        r"Power-law exponent for masses below :math:`m_{\mathrm{max,NS}}`",
681    ),
682    "alpha_2": (
683        r":math:`\alpha_2`",
684        r"Power-law exponent for masses above :math:`m_{\mathrm{min,BH}}`",
685    ),
686    "alpha_dip": (r":math:`\alpha_d`", r"Power-law exponent within the NS–BH mass gap"),
687    "NSmin": (
688        r":math:`m_{\mathrm{min,NS}}`",
689        r"Minimum neutron star mass (:math:`M_\odot`)",
690    ),
691    "NSmax": (
692        r":math:`\gamma_{\mathrm{low},1}`",
693        r"Maximum neutron star mass (:math:`M_\odot`)",
694    ),
695    "BHmin": (
696        r":math:`\gamma_{\mathrm{high},1}`",
697        r"Minimum black hole mass (:math:`M_\odot`)",
698    ),
699    "BHmax": (
700        r":math:`m_{\mathrm{max,BH}}`",
701        r"Maximum black hole mass (:math:`M_\odot`)",
702    ),
703    "A": (
704        r":math:`\mathrm{A}`",
705        r"Depth of primary mass gap suppression",
706    ),
707    "UPPERmin": (
708        r":math:`\gamma_{\mathrm{low},2}`",
709        r"Lower boundary of pair-instability gap (:math:`M_\odot`)",
710    ),
711    "UPPERmax": (
712        r":math:`\gamma_{\mathrm{high},2}`",
713        r"Upper boundary of pair-instability gap (:math:`M_\odot`)",
714    ),
715    "mu1": (
716        r":math:`\mu_{\mathrm{peak},1}`",
717        r"Mean of primary Gaussian peak (:math:`M_\odot`)",
718    ),
719    "sig1": (
720        r":math:`\sigma_{\mathrm{peak},1}`",
721        r"Std. dev. of primary Gaussian peak (:math:`M_\odot`)",
722    ),
723    "mix1": (r":math:`\mathrm{c}_1`", r"Mixing fraction of primary Gaussian peak"),
724    "mu2": (
725        r":math:`\mu_{\mathrm{peak},2}`",
726        r"Mean of secondary Gaussian peak (:math:`M_\odot`)",
727    ),
728    "sig2": (
729        r":math:`\sigma_{\mathrm{peak},2}`",
730        r"Std. dev. of secondary Gaussian peak (:math:`M_\odot`)",
731    ),
732    "mix2": (
733        r":math:`\mathrm{c}_2`",
734        r"Mixing fraction of secondary Gaussian peak",
735    ),
736    "absolute_mmin": (r":math:`m_\mathrm{abs,min}`", r"Absolute minimum truncation"),
737    "absolute_mmax": (r":math:`m_\mathrm{abs,max}`", r"Absolute maximum truncation"),
738    "n0": (
739        r":math:`\eta_0`",
740        r"Sharpness of low-mass truncation",
741    ),
742    "n5": (
743        r":math:`\eta_5`",
744        r"Sharpness of high-mass truncation",
745    ),
746    "n1": (
747        r":math:`\eta_1`",
748        r"Sharpness at :math:`m_{\mathrm{max,NS}}`",
749    ),
750    "n2": (
751        r":math:`\eta_2`",
752        r"Sharpness at :math:`m_{\mathrm{min,BH}}`",
753    ),
754    "n3": (
755        r":math:`\eta_3`",
756        r"Sharpness at :math:`\gamma_{\mathrm{low},2}`",
757    ),
758    "n4": (
759        r":math:`\eta_4`",
760        r"Sharpness at :math:`\gamma_{\mathrm{high},2}`",
761    ),
762}
763
764
765# label -> source file, downloaded into data/raw/ under its official name
766# if missing. Add an entry here (e.g. a future GWTC-6.0 file) and only
767# that one gets fetched/loaded on the next run.
768SOURCES = {
769    "GWTC-4.0 (popsummary)": str(DATA_DIR / "raw" / ZENODO_FILENAME_GWTC4),
770    "GWTC-5.0 (popsummary)": str(DATA_DIR / "raw" / ZENODO_FILENAME_GWTC5),
771}
772
773
774def main(force: bool = False) -> None:
775    """Rebuild data/derived/ from data/raw/, fetching whatever's missing
776    from DCC/Zenodo. force=False (default) reuses existing cache entries;
777    force=True re-extracts everything from data/raw/, downloading first
778    if needed."""
779    download_and_extract_from_zenodo_tarball(
780        ZENODO_RECORD_GWTC4,
781        ZENODO_ARCHIVE_GWTC4,
782        ZENODO_FILENAME_GWTC4,
783        Path(SOURCES["GWTC-4.0 (popsummary)"]),
784    )
785    download_and_extract_from_zenodo_tarball(
786        ZENODO_RECORD_GWTC5,
787        ZENODO_ARCHIVE_GWTC5,
788        ZENODO_FILENAME_GWTC5,
789        Path(SOURCES["GWTC-5.0 (popsummary)"]),
790    )
791    download_and_extract_from_zenodo_tarball(
792        ZENODO_RECORD_GWTC5,
793        ZENODO_ARCHIVE_GWTC5,
794        ZENODO_MEMBER_PIXELPOP,
795        Path(RATE_SOURCES["GWTC-5.0 PixelPop"][0]),
796    )
797
798    hyperparams_csv = DERIVED_DIR / "hyperparams_map.csv"
799    hyperparams_map = extract_hyperparams_map(
800        SOURCES, cache_csv=str(hyperparams_csv), force=force
801    )
802    print(f"Hyperparameters MAP cache: {hyperparams_csv}")
803
804    hyperparams_median_csv = DERIVED_DIR / "hyperparams_median.csv"
805    hyperparams_median = extract_hyperparams_median(
806        SOURCES, cache_csv=str(hyperparams_median_csv), force=force
807    )
808    print(f"Hyperparameters median cache: {hyperparams_median_csv}")
809
810    rate_summary_csv = DERIVED_DIR / "rate_summary.csv"
811    extract_rate_summary(RATE_SOURCES, cache_csv=str(rate_summary_csv), force=force)
812    print(f"Rate summary cache: {rate_summary_csv}")
813
814    published_rates_csv = DERIVED_DIR / "published_rates_table2.csv"
815    write_published_rates_table2(str(published_rates_csv))
816    print(f"Published Table 2 rates cache: {published_rates_csv}")
817
818    # Joint mass-rate grids. Every current SOURCES entry is a popsummary
819    # file, so all of them carry grids.
820    grid_cache_paths = extract_joint_grids(
821        SOURCES, cache_dir=str(DERIVED_DIR), force=force
822    )
823    for label, npz_path in grid_cache_paths.items():
824        print(f"Joint grids cache [{label}]: {npz_path}")
825
826    # Table below is built from the GWTC-5.0 MAP.
827    maxp_samp = hyperparams_map.loc["GWTC-5.0 (popsummary)"]
828
829    rows = []
830    for key, (param, desc) in PARAMS_INFO.items():
831        if key in maxp_samp:
832            val = maxp_samp[key]
833            rows.append((f"{param}", desc, f"{val:.3g}"))
834
835    df = pd.DataFrame(rows, columns=["Parameter", "Description", "Value"])
836
837    title = "Hyperparameters of the FullPop model (GWTC-5.0)"
838
839    # Saved next to this script, wherever it's invoked from.
840    with open(SCRIPT_DIR / "hyperparams_table.rst", "w") as f:
841        f.write(to_rst(df, title=title))
842
843    with open(SCRIPT_DIR / "hyperparams_table.tex", "w") as f:
844        f.write(to_latex(df, caption=title))
845
846    # Paper-style (booktabs, sectioned) table, built straight from the MAP
847    # row rather than the simplified `df` above. Caption says "MAP"
848    # explicitly -- it used not to, which is exactly what caused this to
849    # be mistaken for the posterior-median table down the line.
850    with open(SCRIPT_DIR / "hyperparams_table_gwtc5.tex", "w") as f:
851        f.write(
852            to_latex_paper(
853                maxp_samp,
854                caption="Hyperparameters of the FullPop Distribution Model (GWTC-5.0), MAP.",
855            )
856        )
857
858    # Same table, posterior median instead of MAP -- GWTC-5.0's own
859    # convention (see load_hyperparams_median's docstring). Kept as a
860    # separate file/caption rather than replacing the MAP one above: both
861    # are legitimate summaries of the same posterior, and collapsing them
862    # into one file under one name is what caused the confusion this is
863    # fixing.
864    median_samp = hyperparams_median.loc["GWTC-5.0 (popsummary)"]
865    with open(SCRIPT_DIR / "hyperparams_table_gwtc5_median.tex", "w") as f:
866        f.write(
867            to_latex_paper(
868                median_samp,
869                caption="Hyperparameters of the FullPop Distribution Model (GWTC-5.0), posterior median.",
870            )
871        )
872
873
874if __name__ == "__main__":
875    main()

Run the Pipeline#

Running the pipeline

Two equivalent options are available - choose one:

Run commands without activating the environment explicitly:

uv run make

Activate the .venv created by uv and run the pipeline:

source .venv/bin/activate
make