Simulating a Weak-Lensing Galaxy Catalog

Author

NumCosmo developers

Purpose

Builds a mock weak-lensing galaxy catalog from the position, redshift, and shape factor calculators described in The Galaxy Population & Weak-Lensing Framework. Read that page first for the vocabulary used here: Factor, Pop, Obs, and the per-galaxy Data.

The orchestrator that wires the three factors together has not been built yet, so this tutorial threads the per-galaxy data between them by hand. That also shows exactly what each factor needs from its neighbours.

Setting Up the Lens and the Models

We need a lens (a single NFW halo at a fixed redshift and mass, observed through its projected surface mass density) and three population models: the redshift distribution of the source sample, its photo-\(z\) error kernel, and the intrinsic-ellipticity distribution of the sources.

import numpy as np
from IPython.display import HTML
from numcosmo_py import Nc, Ncm
from numcosmo_py.catalog import catalog_to_table

Ncm.cfg_init()


def show_table(table, float_format="%.4f", max_rows=10):
    """Render an astropy Table as a compact HTML table."""
    return HTML(table.to_pandas().to_html(float_format=float_format, max_rows=max_rows))


cosmo = Nc.HICosmoDEXcdm.new()
dist = Nc.Distance.new(100.0)

hms = Nc.HaloCMParam.new(Nc.HaloMassSummaryMassDef.MEAN, 200.0)
hms.param_set_by_name("log10MDelta", 14.0)
density_profile = Nc.HaloDensityProfileNFW.new(hms)

halo_position = Nc.HaloPosition.new(dist)
halo_position.param_set_by_name("z", 0.3)
halo_position.prepare(cosmo)

surface_mass_density = Nc.WLSurfaceMassDensity.new(dist)

pop_z = Nc.GalaxyRedshiftPopLSSTSRD.new_y1_source()
obs_z = Nc.GalaxyRedshiftObsGauss.new()

pop_shape = Nc.GalaxyShapePopGauss.new()
pop_shape["sigma"] = 0.3

All seven models are collected into one NcmMSet: every calculator below reads whichever of these it needs by resolving them from this same mset, rather than holding a reference of its own.

mset = Ncm.MSet.empty_new()
for model in (cosmo, density_profile, halo_position, surface_mass_density,
              pop_z, obs_z, pop_shape):
    mset.set(model)

Building the Three Factors

Each factor is a plain calculator (not a model, not in the mset): a flat footprint for the position, the composed true-\(z \otimes\) photo-\(z\) scheme for the redshift, and the variance-add marginalization for the shape.

position_factor = Nc.GalaxyPositionFactorFlat.new(-0.2, 0.2, -0.2, 0.2)
redshift_factor = Nc.GalaxyRedshiftFactorComposed.new(0.0, 20.0)
shape_factor = Nc.GalaxyShapeFactorVarAdd.new(Nc.GalaxyWLObsEllipConv.TRACE_DET)

Each factor’s data_new allocates the per-galaxy fragment for that factor. The shape factor’s fragment is built from the position and redshift fragments, since the shape likelihood needs the projected radius (from ra/dec) and the source redshift z (see the framework page’s Z \(\Rightarrow\) P \(\Rightarrow\) S pipeline) — this one call is the whole “orchestration” this tutorial needs to do by hand:

pos_data = Nc.GalaxyPositionFactorData.new(position_factor, mset)
z_data = Nc.GalaxyRedshiftFactorData.new(redshift_factor, mset)
shape_data = Nc.GalaxyShapeFactorData.new(
    shape_factor, mset, pos_data, z_data
)

Generating the Catalog

Every factor exposes the columns its fragment needs through required_columns; concatenating them (skipping repeats) gives exactly the columns a storage catalog must provide, with no factor needing to know the others’ column names:

columns = []
for cols in (
    Nc.GalaxyPositionFactorData.required_columns(pos_data),
    Nc.GalaxyRedshiftFactorData.required_columns(z_data),
    Nc.GalaxyShapeFactorData.required_columns(shape_data),
):
    for col in cols:
        if col not in columns:
            columns.append(col)

print(columns)
['ra', 'dec', 'z', 'zp', 'sigma0', 'epsilon_int_1', 'epsilon_int_2', 'epsilon_obs_1', 'epsilon_obs_2', 'std_noise', 'c1', 'c2', 'm']

The storage object is NcGalaxyWLObs — a plain named-column catalog, agnostic to which calculators fill it in.

n_galaxies = 8
wl_obs = Nc.GalaxyWLObs.new(
    Nc.GalaxyWLObsEllipConv.TRACE_DET,
    Nc.WLEllipticityFrame.CELESTIAL,
    n_galaxies,
    columns,
)

For each galaxy we draw a position and a redshift, fix the noise/calibration inputs (std_noise, c1, c2, m — these stay plain scalars on the shape factor’s data, not a pluggable model, for now), draw the shape, and write the three fragments into their row. gen is the same code path integ uses below for the likelihood — see Generation and Likelihood Share One Engine.

rng = Ncm.RNG.seeded_new(None, 123)

for i in range(n_galaxies):
    position_factor.gen(mset, pos_data, rng)
    redshift_factor.gen(mset, z_data, rng)

    shape_factor.data_set(
        shape_data, 0.0, 0.0, 0.03, 0.0, 0.0, 0.0,
        Nc.WLEllipticityFrame.CELESTIAL,
    )
    shape_factor.gen(mset, shape_data, rng)

    pos_data.write_row(wl_obs, i)
    z_data.write_row(wl_obs, i)
    shape_data.write_row(wl_obs, i)

NcGalaxyWLObs derives from NcmCatalog, so the generic catalog_to_table helper (used in the halo/cluster mock tutorial) works here unchanged:

Code 
show_table(catalog_to_table(wl_obs))
ra dec z zp sigma0 epsilon_int_1 epsilon_int_2 epsilon_obs_1 epsilon_obs_2 std_noise c1 c2 m
0 0.0852 0.0786 0.3102 0.3102 0.0000 0.4265 -0.0163 0.4787 -0.0188 0.0300 0.0000 0.0000 0.0000
1 0.0319 0.0739 0.1700 0.1700 0.0000 -0.4682 0.0472 -0.4322 0.0201 0.0300 0.0000 0.0000 0.0000
2 0.0376 -0.0408 0.2626 0.2626 0.0000 -0.4247 -0.0081 -0.3701 -0.0827 0.0300 0.0000 0.0000 0.0000
3 0.0961 0.0898 0.1354 0.1354 0.0000 -0.2198 0.3891 -0.1855 0.3490 0.0300 0.0000 0.0000 0.0000
4 0.1108 -0.1632 0.3993 0.3993 0.0000 0.2055 0.4641 0.1693 0.4415 0.0300 0.0000 0.0000 0.0000
5 0.0774 0.1574 1.2337 1.2337 0.0000 -0.0089 0.5179 0.0256 0.5060 0.0300 0.0000 0.0000 0.0000
6 0.1667 0.1465 0.8541 0.8541 0.0000 -0.0158 0.1423 -0.0976 0.1831 0.0300 0.0000 0.0000 0.0000
7 -0.0869 0.0412 1.0196 1.0196 0.0000 -0.1114 -0.9073 -0.0664 -0.8821 0.0300 0.0000 0.0000 0.0000
Table 1: Mock weak-lensing galaxy catalog.

Reading a Row Back

read_row is the inverse of write_row: it re-populates a fragment from one row of the catalog. This is exactly what a likelihood evaluation does when loading an observed catalog rather than a generated one — a fragment does not care whether the values it is holding were drawn or measured.

pos_data_readback = Nc.GalaxyPositionFactorData.new(position_factor, mset)
pos_data_readback.read_row(wl_obs, 0)

print(f"row 0: ra={pos_data_readback.ra:.4f}, dec={pos_data_readback.dec:.4f}")
row 0: ra=0.0852, dec=0.0786

Evaluating the Shape Likelihood Factor

The same shape_factor that generated the mock also evaluates the per-galaxy likelihood factor \(p(\epsilon_\mathrm{obs} \mid z)\) as a function of the source redshift. prepare_data_array resolves the per-galaxy geometry caches described on the framework page (projected radius, surface-density optimization structs) once, and integ returns the reusable integrand callback:

shape_factor.prepare_data_array(mset, [shape_data], True, True)
integrand = shape_factor.integ(mset, False)

z_grid = np.linspace(0.05, 1.5, 5)
for z in z_grid:
    print(f"z={z:.2f}  ->  p(eps_obs | z) = {integrand.eval(z, shape_data):.4e}")
z=0.05  ->  p(eps_obs | z) = 2.3663e-02
z=0.41  ->  p(eps_obs | z) = 2.3470e-02
z=0.78  ->  p(eps_obs | z) = 2.3232e-02
z=1.14  ->  p(eps_obs | z) = 2.3148e-02
z=1.50  ->  p(eps_obs | z) = 2.3106e-02

In a full likelihood this integrand would be combined with the redshift factor’s own \(p(z_p, z)\) and integrated over \(z\) once per galaxy — the single outer integral the orchestrator (once built) will perform automatically.

Summary

We built a lens and three population models, wired one factor per observable (position, redshift, shape), generated a mock catalog by calling each factor’s gen and writing the results into a shared NcGalaxyWLObs, read a row back to confirm the round trip, and evaluated the shape factor’s own likelihood integrand for the very galaxies it had just generated. See The Galaxy Population & Weak-Lensing Framework for why the pieces are shaped this way.

API Reference

The catalog_to_table helper lives in the numcosmo_py.catalog subpackage.