"""Excess surface mass density ΔΣ(r_p) for galaxy-galaxy lensing.
Provides:
- Production measurement via the dsigma pipeline: :func:`delta_sigma`
The :func:`delta_sigma` function follows the same steps used for DES Y3,
HSC, and KiDS DR4 measurements:
1. Lens-source pair precomputation with optional photo-z n(z) calibration.
2. Random-point subtraction to suppress additive systematics.
3. Adaptive jackknife resampling for covariance estimation.
4. Stacking with survey-specific shear corrections.
Use the survey loaders in :mod:`sum_stat.lensing.surveys` to prepare the
source and n(z) tables, then pass them directly to :func:`delta_sigma`.
"""
from __future__ import annotations
import numpy as np
from astropy.cosmology import FlatLambdaCDM
from astropy.table import Table
from ..catalogue import GalaxyCatalogue
def _adaptive_n_jackknife(n_lenses: int) -> int:
"""Return jackknife region count scaled to the lens sample size."""
if n_lenses < 100:
return 10
if n_lenses < 500:
return 50
return 100
def _galaxy_catalogue_to_dsigma_table(cat: GalaxyCatalogue) -> Table:
"""Convert a :class:`GalaxyCatalogue` to a dsigma-formatted lens table."""
import dsigma.helpers as dh
raw = Table(
{
"ra": cat.ra,
"dec": cat.dec,
"z": cat.redshift,
"w_sys": cat.weight,
}
)
return dh.dsigma_table(raw, "lens")
[docs]
def delta_sigma(
lens: GalaxyCatalogue,
table_s: Table,
rp_bins: np.ndarray,
cosmo: FlatLambdaCDM,
*,
rand: GalaxyCatalogue | None = None,
table_n: Table | None = None,
lens_source_cut: float = 0.1,
n_jackknife: int | None = None,
comoving: bool = True,
corrections: dict | None = None,
) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Excess surface mass density ΔΣ(r_p) via the dsigma pipeline.
Follows the same steps as the DES Y3, HSC, and KiDS DR4 measurement
scripts: precompute lens-source pair sums, subtract randoms, assign
jackknife regions adaptively, and stack with survey-specific corrections.
Parameters
----------
lens : GalaxyCatalogue
Lens catalogue (ra, dec, redshift, weight required).
table_s : Table
dsigma-formatted source shape table. Prepare with the survey-specific
loaders in :mod:`sum_stat.lensing.surveys` (e.g.
:func:`~sum_stat.lensing.surveys.load_des_sources`).
rp_bins : ndarray, shape (n_bins+1,)
Projected separation bin edges [Mpc].
cosmo : FlatLambdaCDM
Cosmology for Σ_crit computation.
rand : GalaxyCatalogue, optional
Random-point catalogue used for additive-bias subtraction. Required
when ``corrections`` contains ``random_subtraction=True``.
table_n : Table, optional
Photo-z calibration n(z) table (columns ``z`` and ``n``). Used by
dsigma to correct for photo-z dilution when ``lens_source_cut`` is
applied. Returned by the survey loaders alongside ``table_s``.
lens_source_cut : float
Minimum redshift separation Δz = z_s − z_l required to include a
lens-source pair. Typical values: 0.1 (DES, KiDS), 0.3 (HSC).
n_jackknife : int, optional
Number of jackknife regions. When *None* (default), the count is
chosen adaptively: 10 for N_lens < 100, 50 for N_lens < 500, 100
otherwise — matching the reference scripts.
comoving : bool
Use comoving coordinates for Σ_crit (default True).
corrections : dict, optional
Boolean correction flags forwarded to
``dsigma.stacking.excess_surface_density``. Survey-specific defaults:
* **DES Y3** ::
{'matrix_shear_response_correction': True,
'random_subtraction': True}
* **HSC** ::
{'scalar_shear_response_correction': True,
'shear_responsivity_correction': True,
'selection_bias_correction': True,
'boost_correction': False,
'random_subtraction': True}
* **KiDS DR4** ::
{'scalar_shear_response_correction': True,
'random_subtraction': True}
Returns
-------
rp_centres : ndarray, shape (n_bins,)
Projected separation bin centres [Mpc] (geometric mean of bin edges).
ds : ndarray, shape (n_bins,)
ΔΣ(r_p) [M_sun/pc^2].
cov : ndarray, shape (n_bins, n_bins)
Jackknife covariance [(M_sun/pc^2)^2].
Raises
------
ValueError
If ``corrections`` requests random subtraction but ``rand`` is None.
ValueError
If fewer than 10 lenses remain after the precomputation filter.
Performance
-----------
~30 s/call (N_lens=1000, N_src=3000, n_bins=10, n_jk=100, CPU x86-64)
Examples
--------
DES Y3 measurement:
>>> from sum_stat.lensing.surveys import load_des_sources
>>> from sum_stat import GalaxyCatalogue
>>> from astropy.cosmology import Planck15
>>> import numpy as np
>>>
>>> table_s, table_n = load_des_sources('des_y3.hdf5')
>>> lens = GalaxyCatalogue(ra=..., dec=..., redshift=...)
>>> rand = GalaxyCatalogue(ra=..., dec=..., redshift=...)
>>> rp_bins = 10 ** np.arange(-2.0, 1.6, 0.4)
>>> rp, ds, cov = delta_sigma(
... lens, table_s, rp_bins, Planck15,
... rand=rand, table_n=table_n,
... corrections={'matrix_shear_response_correction': True,
... 'random_subtraction': True},
... )
"""
import dsigma.jackknife as djk
import dsigma.precompute as dp
import dsigma.stacking as ds_stack
if corrections is None:
corrections = {}
if corrections.get("random_subtraction", False) and rand is None:
raise ValueError("corrections['random_subtraction']=True requires a rand catalogue.")
rp_bins = np.asarray(rp_bins, dtype=np.float64)
rp_centres = np.sqrt(rp_bins[:-1] * rp_bins[1:])
# Convert GalaxyCatalogue → dsigma astropy Table
table_l = _galaxy_catalogue_to_dsigma_table(lens)
table_r = _galaxy_catalogue_to_dsigma_table(rand) if rand is not None else None
# Precompute lens-source and random-source pair sums
precompute_kwargs = dict(
cosmology=cosmo,
comoving=comoving,
lens_source_cut=lens_source_cut,
progress_bar=False,
)
if table_n is not None:
precompute_kwargs["table_n"] = table_n
dp.precompute(table_l, table_s, rp_bins, **precompute_kwargs)
if table_r is not None:
dp.precompute(table_r, table_s, rp_bins, **precompute_kwargs)
# Remove lenses/randoms that had no nearby source
table_l = table_l[np.sum(table_l["sum 1"], axis=1) > 0]
if table_r is not None:
table_r = table_r[np.sum(table_r["sum 1"], axis=1) > 0]
if len(table_l) < 10:
raise ValueError(
f"Only {len(table_l)} lenses survived the source-match filter; "
"too few for a meaningful measurement."
)
# Adaptive jackknife regions
n_jk = n_jackknife if n_jackknife is not None else _adaptive_n_jackknife(len(table_l))
centers = djk.compute_jackknife_fields(table_l, n_jk, weights=np.sum(table_l["sum 1"], axis=1))
if table_r is not None:
djk.compute_jackknife_fields(table_r, centers)
# Build stacking kwargs, injecting the random table when required
stack_kwargs = dict(corrections)
stack_kwargs["return_table"] = True
if table_r is not None:
stack_kwargs["table_r"] = table_r
# Stack and compute jackknife covariance
result = ds_stack.excess_surface_density(table_l, **stack_kwargs)
ds_vals = np.array(result["ds"])
stack_kwargs["return_table"] = False
cov = np.array(
djk.jackknife_resampling(ds_stack.excess_surface_density, table_l, **stack_kwargs)
)
return rp_centres, ds_vals, cov