Simulating a Weak-Lensing Galaxy Catalog

Author

NumCosmo developers

Purpose

This tutorial builds a mock weak-lensing galaxy catalog directly from the per-observation factor calculators described in The Galaxy Population & Weak-Lensing Framework: one for the sky position, one for the redshift, one for the shape. Reading that page first will make the vocabulary here (Factor, Pop, Obs, per-galaxy Data) much less cryptic; this page is deliberately about using the pieces, not about why they are shaped this way.

Because the orchestrator that wires the three factors automatically has not been built yet (see the framework page’s status note), this tutorial threads the per-galaxy data between them by hand — which is also the best way to see exactly what each factor needs from its neighbors.

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 ⇒ P ⇒ 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 freshly 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.