"""Catalogue dataclasses for galaxy, shape, and photo-z calibration data."""
from __future__ import annotations
from dataclasses import dataclass
import numpy as np
from astropy.cosmology import FlatLambdaCDM
[docs]
@dataclass
class GalaxyCatalogue:
"""Observed and derived properties of a galaxy or random catalogue.
Parameters
----------
ra : array_like
Right ascension [degrees].
dec : array_like
Declination [degrees].
redshift : array_like
Spectroscopic or photometric redshift (dimensionless).
weight : array_like, optional
Combined completeness and systematic weight (w_comp * w_sys).
Defaults to an array of ones.
mag : array_like, optional
Observed magnitudes, shape (N, n_bands).
flux : array_like, optional
Emission line fluxes [erg/s/cm^2], shape (N, n_lines).
abs_mag : array_like, optional
Absolute magnitudes, shape (N,) or (N, n_bands).
log10_mstar : array_like, optional
Stellar mass log10(M_* / M_sun).
log10_sfr : array_like, optional
Star formation rate log10(SFR / (M_sun yr^-1)).
log10_ssfr : array_like, optional
Specific SFR log10(sSFR / yr^-1).
mag_limit : float, optional
Apparent magnitude limit of the survey (for V_max calculation).
flux_limit : float, optional
Flux limit [erg/s/cm^2] of the survey (for V_max calculation).
"""
ra: np.ndarray
dec: np.ndarray
redshift: np.ndarray
weight: np.ndarray | None = None
mag: np.ndarray | None = None
flux: np.ndarray | None = None
abs_mag: np.ndarray | None = None
log10_mstar: np.ndarray | None = None
log10_sfr: np.ndarray | None = None
log10_ssfr: np.ndarray | None = None
mag_limit: float | None = None
flux_limit: float | None = None
def __post_init__(self) -> None:
self.ra = np.asarray(self.ra, dtype=np.float64)
self.dec = np.asarray(self.dec, dtype=np.float64)
self.redshift = np.asarray(self.redshift, dtype=np.float64)
n = self.ra.size
if self.dec.size != n or self.redshift.size != n:
raise ValueError("ra, dec, redshift must have the same length.")
if self.weight is None:
self.weight = np.ones(n, dtype=np.float64)
else:
self.weight = np.asarray(self.weight, dtype=np.float64)
if self.weight.size != n:
raise ValueError("weight must have the same length as ra.")
for attr in ("mag", "flux", "abs_mag", "log10_mstar", "log10_sfr", "log10_ssfr"):
val = getattr(self, attr)
if val is not None:
arr = np.asarray(val, dtype=np.float64)
if arr.shape[0] != n:
raise ValueError(f"{attr} first dimension must match ra length ({n}).")
setattr(self, attr, arr)
@property
def n(self) -> int:
"""Number of objects."""
return self.ra.size
[docs]
def subsample(self, n_max: int) -> GalaxyCatalogue:
"""Return a random subsample of at most n_max objects.
Parameters
----------
n_max : int
Maximum number of objects to keep.
Returns
-------
GalaxyCatalogue
Subsampled catalogue (same object if n <= n_max).
"""
if self.n <= n_max:
return self
idx = np.sort(np.random.default_rng().choice(self.n, size=n_max, replace=False))
def _s(arr):
return arr[idx] if arr is not None else None
return GalaxyCatalogue(
ra=self.ra[idx],
dec=self.dec[idx],
redshift=self.redshift[idx],
weight=self.weight[idx],
mag=_s(self.mag),
flux=_s(self.flux),
abs_mag=_s(self.abs_mag),
log10_mstar=_s(self.log10_mstar),
log10_sfr=_s(self.log10_sfr),
log10_ssfr=_s(self.log10_ssfr),
mag_limit=self.mag_limit,
flux_limit=self.flux_limit,
)
[docs]
def comoving_distance(self, cosmo: FlatLambdaCDM) -> np.ndarray:
"""Comoving distance to each object [Mpc].
Parameters
----------
cosmo : FlatLambdaCDM
Astropy cosmology instance.
Returns
-------
chi : ndarray, shape (N,)
Comoving distances [Mpc].
"""
return cosmo.comoving_distance(self.redshift).to("Mpc").value
[docs]
def vmax(
self,
z_min: float,
z_max: float,
area_sr: float,
cosmo: FlatLambdaCDM,
z_max_individual: np.ndarray | None = None,
) -> np.ndarray:
"""Comoving V_max for each object [Mpc^3].
V_max_i = (area_sr / 4pi) * [V_c(min(z_max, z_max_i)) - V_c(z_min)]
Parameters
----------
z_min : float
Minimum redshift of the survey volume.
z_max : float
Maximum redshift of the survey volume.
area_sr : float
Survey solid angle [sr].
cosmo : FlatLambdaCDM
Astropy cosmology instance.
z_max_individual : array_like, optional
Maximum redshift at which each object would still enter the sample
(e.g. from the mag/flux limit). Defaults to z_max for all objects.
Returns
-------
vmax : ndarray, shape (N,)
V_max values [Mpc^3].
"""
if z_max_individual is None:
z_max_individual = np.full(self.n, z_max)
else:
z_max_individual = np.asarray(z_max_individual, dtype=np.float64)
z_eff = np.minimum(z_max, z_max_individual)
vc_max = cosmo.comoving_volume(z_eff).to("Mpc3").value
vc_min = cosmo.comoving_volume(z_min).to("Mpc3").value
return (area_sr / (4.0 * np.pi)) * (vc_max - vc_min)
[docs]
@dataclass
class ShapeCatalogue:
"""Weak-lensing shape catalogue.
Parameters
----------
ra : array_like
Right ascension [degrees].
dec : array_like
Declination [degrees].
z_phot : array_like
Photometric redshift (dimensionless).
e1, e2 : array_like
Ellipticity components (dimensionless).
w : array_like
Per-galaxy lensing weight.
m : array_like
Multiplicative shear bias per galaxy (dimensionless).
e_rms : array_like
Per-galaxy rms ellipticity (dimensionless).
R11, R22 : array_like
Diagonal elements of the shear response matrix (dimensionless).
"""
ra: np.ndarray
dec: np.ndarray
z_phot: np.ndarray
e1: np.ndarray
e2: np.ndarray
w: np.ndarray
m: np.ndarray
e_rms: np.ndarray
R11: np.ndarray
R22: np.ndarray
def __post_init__(self) -> None:
arrays = [
self.ra,
self.dec,
self.z_phot,
self.e1,
self.e2,
self.w,
self.m,
self.e_rms,
self.R11,
self.R22,
]
names = ["ra", "dec", "z_phot", "e1", "e2", "w", "m", "e_rms", "R11", "R22"]
n = None
for name, arr in zip(names, arrays):
arr = np.asarray(arr, dtype=np.float64)
object.__setattr__(self, name, arr)
if n is None:
n = arr.size
elif arr.size != n:
raise ValueError(f"ShapeCatalogue: {name} length {arr.size} != {n}.")
self._n = n
@property
def n(self) -> int:
"""Number of source galaxies."""
return self._n
[docs]
def mean_m(self) -> float:
"""Weighted mean multiplicative shear bias <m>."""
return float(np.average(self.m, weights=self.w))
[docs]
def mean_R(self) -> tuple[float, float]:
"""Weighted mean shear response (R11, R22)."""
R11_mean = float(np.average(self.R11, weights=self.w))
R22_mean = float(np.average(self.R22, weights=self.w))
return R11_mean, R22_mean
[docs]
@dataclass
class PhotoZCalibTable:
"""Photo-z calibration table mapping z_phot to z_true with weights.
Used for photo-z bias correction in galaxy-galaxy lensing (dsigma pipeline).
Parameters
----------
z_phot : array_like
Photometric redshift.
z_true : array_like
True (spectroscopic) redshift.
w : array_like
Per-object weight.
w_sys : array_like, optional
Systematic weight. Defaults to ones.
"""
z_phot: np.ndarray
z_true: np.ndarray
w: np.ndarray
w_sys: np.ndarray | None = None
def __post_init__(self) -> None:
self.z_phot = np.asarray(self.z_phot, dtype=np.float64)
self.z_true = np.asarray(self.z_true, dtype=np.float64)
self.w = np.asarray(self.w, dtype=np.float64)
n = self.z_phot.size
if self.z_true.size != n or self.w.size != n:
raise ValueError("z_phot, z_true, w must have the same length.")
if self.w_sys is None:
self.w_sys = np.ones(n, dtype=np.float64)
else:
self.w_sys = np.asarray(self.w_sys, dtype=np.float64)
if self.w_sys.size != n:
raise ValueError("w_sys must have the same length as z_phot.")
@property
def n(self) -> int:
"""Number of calibration objects."""
return self.z_phot.size
[docs]
def n_of_z(
self,
z_bins: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
"""Weighted n(z) histogram of the true redshifts.
Parameters
----------
z_bins : array_like
Bin edges for z_true, shape (n_bins+1,).
Returns
-------
z_centres : ndarray, shape (n_bins,)
n_z : ndarray, shape (n_bins,)
Normalised n(z) (sums to 1).
"""
z_bins = np.asarray(z_bins, dtype=np.float64)
counts, _ = np.histogram(self.z_true, bins=z_bins, weights=self.w * self.w_sys)
z_centres = 0.5 * (z_bins[:-1] + z_bins[1:])
norm = counts.sum()
if norm > 0:
counts = counts / norm
return z_centres, counts