"""HDF5 reader for summary statistics, reconstructing astropy Quantities."""
from __future__ import annotations
from pathlib import Path
import h5py
import numpy as np
from astropy.cosmology import FlatLambdaCDM
[docs]
class SummaryStatReader:
"""Read summary statistics from an HDF5 file written by :class:`SummaryStatWriter`.
Parameters
----------
path : str or Path
Path to the HDF5 file.
Examples
--------
>>> reader = SummaryStatReader("results.h5")
>>> groups = reader.list_statistics()
>>> result = reader.read_esd("esd/DESI_BGS_KiDS")
>>> result["delta_sigma"] # ndarray in M_sun/pc^2
>>> result["cov"] # covariance matrix
>>> result["cosmology"] # FlatLambdaCDM instance
>>> reader.close()
"""
def __init__(self, path: str | Path) -> None:
self._path = Path(path)
self._file = h5py.File(self._path, "r")
# ------------------------------------------------------------------
# Navigation
# ------------------------------------------------------------------
[docs]
def list_statistics(self) -> list[str]:
"""Return all leaf group paths (summary-statistic groups) in the file.
Returns
-------
paths : list of str
"""
paths: list[str] = []
def _visitor(name: str, obj: h5py.Group | h5py.Dataset) -> None:
if isinstance(obj, h5py.Group) and "cosmology" in obj:
paths.append(name)
elif isinstance(obj, h5py.Group) and "ell_eff" in obj:
paths.append(name)
self._file.visititems(_visitor)
return sorted(paths)
# ------------------------------------------------------------------
# Typed readers
# ------------------------------------------------------------------
[docs]
def read_lf(self, group: str) -> dict:
"""Read a luminosity function group.
Returns
-------
dict with keys: M_centres, bin_edges, phi, phi_err, cov,
cosmology (FlatLambdaCDM), attrs (dict).
"""
return self._read_group(group, ["M_centres", "bin_edges", "phi", "phi_err", "cov"])
[docs]
def read_smf(self, group: str) -> dict:
"""Read a stellar mass function group."""
return self._read_group(group, ["log10mstar_centres", "bin_edges", "phi", "phi_err", "cov"])
[docs]
def read_twopcf(self, group: str) -> dict:
"""Read a two-point correlation function group."""
return self._read_group(group, ["sep_centres", "bin_edges", "xi", "cov"])
[docs]
def read_multipoles(self, group: str) -> dict:
"""Read a multipoles group, returning xi0, xi2, xi4 if present."""
keys = ["s_centres", "bin_edges", "cov"]
result = self._read_group(group, keys)
g = self._file[group]
for ell in (0, 2, 4):
key = f"xi{ell}"
if key in g:
result[key] = np.array(g[key])
return result
[docs]
def read_pk(self, group: str) -> dict:
"""Read a power spectrum group, returning pk0, pk2, pk4 if present."""
keys = ["k_centres", "bin_edges", "cov"]
result = self._read_group(group, keys)
g = self._file[group]
for ell in (0, 2, 4):
key = f"pk{ell}"
if key in g:
result[key] = np.array(g[key])
return result
[docs]
def read_cl(self, group: str) -> dict:
"""Read an angular power spectrum group."""
return self._read_group(group, ["ell_eff", "bin_edges", "cl", "nl", "cov"])
[docs]
def read_esd(self, group: str) -> dict:
"""Read an excess surface mass density group."""
return self._read_group(group, ["rp_centres", "bin_edges", "delta_sigma", "cov"])
[docs]
def cosmology(self, group: str) -> FlatLambdaCDM:
"""Reconstruct the astropy cosmology stored in a group.
Parameters
----------
group : str
HDF5 group path containing a ``cosmology`` subgroup.
Returns
-------
FlatLambdaCDM
"""
g = self._file[group]["cosmology"]
H0 = float(g["H0"][()])
Om0 = float(g["Om0"][()])
Ob0_val = float(g["Ob0"][()])
name = g.attrs.get("name", None) or None
kwargs: dict = {"H0": H0, "Om0": Om0, "name": name}
if not np.isnan(Ob0_val):
kwargs["Ob0"] = Ob0_val
return FlatLambdaCDM(**kwargs)
# ------------------------------------------------------------------
# Context manager
# ------------------------------------------------------------------
[docs]
def close(self) -> None:
"""Close the HDF5 file."""
self._file.close()
def __enter__(self) -> SummaryStatReader:
return self
def __exit__(self, *args) -> None:
self.close()
# ------------------------------------------------------------------
# Private helpers
# ------------------------------------------------------------------
def _read_group(self, group: str, keys: list[str]) -> dict:
g = self._file[group]
result: dict = {}
for key in keys:
if key in g:
result[key] = np.array(g[key])
# Metadata attributes
result["attrs"] = dict(g.attrs)
# Cosmology sub-group
if "cosmology" in g:
result["cosmology"] = self.cosmology(group)
return result