Working tomographic n(z): SOM calibration in sp_validation, read by cosmo_val as a file

config/columns/shapepipe_v2.yaml

@@ -28,6 +28,7 @@ galaxy:
     XWIN_WORLD: SExtractor windowed right ascension [deg]; tile-overlap cut.
     YWIN_WORLD: SExtractor windowed declination [deg]; tile-overlap cut.
     TILE_ID: Tile the object was detected on.
+    NUMBER: Object number within its tile.
     FLAGS: SExtractor extraction flags.
     MAG_AUTO: SExtractor magnitude.
     N_EPOCH: Number of exposures covering the object.

cosmo_val/cat_config.yaml

@@ -1292,3 +1292,67 @@ GLASS_mock_validation:
     e2_col: e2
     path: /n17data/UNIONS/WL/v1.6.x/unions_shapepipe_star_2024_v1.6.a.fits
   patch_number: 150
+SP_v2.0:
+  host: nibi
+  blind: unions
+  subdir: /project/6001537/cdaley/cosmo-val-v2-tomo/catalogue
+  pipeline: SP
+  colour: violet
+  getdist_colour: 0.0, 0.5, 1.0
+  ls: solid
+  marker: '*'
+  patch_number: 100
+  cov_th:
+    A: 4689.662867153573
+    n_e: 6.86756988279185
+    sigma_e: 0.28908294714031374
+    sigma_e_convention: per_component
+    n_psf: 0.580817039505617
+    tomo:
+      A: 3814.775466932069
+      mask: /project/6001537/cdaley/cosmo-val-v2-tomo/catalogue/footprint-tomo-nside4096-ring.fits
+      n_e:
+      - 0.7427992141190386
+      - 0.9163698955600097
+      - 1.006977964461924
+      - 1.3252290285029567
+      sigma_e:
+      - 0.28814703635003425
+      - 0.2851918267345306
+      - 0.28830578863062367
+      - 0.27332974092856765
+  mask: /project/6001537/cdaley/cosmo-val-v2-tomo/catalogue/footprint-nside4096-ring.fits
+  shear:
+    path: shape_catalog_bound9_tomo.fits
+    R: 1.0
+    ra_col: RA
+    dec_col: Dec
+    w_col: w_des
+    e1_col: e1
+    e2_col: e2
+    e1_PSF_col: e1_PSF
+    e2_PSF_col: e2_PSF
+    # Calibration recipe: cosmo_val/nz/v2_bound9.yaml
+    redshift_path: /project/def-mjhudson/cdaley/tomo-nz-v2-full/products/nz_v2_bound9_full_tomo.txt
+    tomo_bin_col: TOM_BIN_ID
+  psf:
+    path: stars.hdf5
+    ra_col: RA
+    dec_col: DEC
+    e1_PSF_col: HSM_G1_PSF
+    e2_PSF_col: HSM_G2_PSF
+    e1_star_col: HSM_G1_STAR
+    e2_star_col: HSM_G2_STAR
+    PSF_size: HSM_T_PSF
+    star_size: HSM_T_STAR
+    PSF_flag: HSM_FLAG_PSF
+    star_flag: HSM_FLAG_STAR
+    mask_config: ../config/calibration/mask_v2.0.yaml
+    mask_cuts:
+    - MASK_8_Manual
+  star:
+    path: psf_catalog_ngmix.fits
+    ra_col: RA
+    dec_col: Dec
+    e1_col: e1
+    e2_col: e2

cosmo_val/nz/README.md

@@ -0,0 +1,40 @@
+# Provisional v2 n(z)
+
+`v2_bound9.yaml` records the four-bin SOM calibration recipe for the full-area bound-9 selection.
+The output has z in column zero and four densities whose integrals are the retained shear-weight fractions; their sum has unit integral.
+The n(z) describes only `TOM_BIN_ID > 0` objects, so a non-tomographic run over all objects pairs an all-object sample with a gold-only n(z).
+Source weights come from the 22,200-tile selection (191,016,178 rows), rather than the 19,940-tile catalogue measured by cosmo_val (170.6M rows).
+cosmo_val reads the saved file through `shear.redshift_path` and never runs calibration or loads photometry.
+The `SP_v2.0` entry uses Nibi catalogue and n(z) paths and declares `host: nibi`.
+The candide n(z) is `/n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/calibration/nz_tomo.txt`.
+At another site, set `shear.redshift_path` to the local copy.
+
+This EGS-only calibration is for validation and null tests.
+It has no independent-field variance or model of spectroscopic selection and does not support parameter constraints.
+
+Run from this checkout inside the sp_validation container on a Slurm allocation:
+
+```bash
+python scripts/prepare_nz_photometry.py cosmo_val/nz/v2_bound9.yaml
+python scripts/calibrate_nz.py cosmo_val/nz/v2_bound9.yaml --workers 16
+```
+
+Photometry preparation is a separate step; reuse its cache for subsequent calibrations.
+The calibration uses frozen SOM cells, groups and gold masks, and computes group weights over the full source selection.
+`--rebuild-model` retrains the SOM, Ward groups and gold masks from the inputs.
+`--output DIRECTORY` selects an output directory; `--figures DIRECTORY` renders diagnostics.
+The frozen-delivery reproduction gate lives in `scratch/cdaley/nz_som_reproduction/`.
+
+The calibration sample, `inputs/matched_calibration.npz`, contains the v1.4.6.3 sample's fiducial spectroscopic `Z_B` as `z`, DR6 BPZ `Z_B` as `photoz`, and detected ugriz magnitudes as `mag`.
+The selected shear Parquet and its provenance are `inputs/bound9_full_selection.parquet` and `inputs/selection.json` in the analysis directory named by the recipe.
+Photometry matching uses the nearest DR6 position within 0.8 arcsec in the same tile and preserves source row order.
+ShapePipe `NUMBER` identifies the source within its tile.
+BPZ stays float64 to preserve values on bin edges.
+
+The recipe fixes a 101×101, 60-epoch batch SOM of all ten ugriz colours plus r magnitude, standardized per feature, and 2,000 Ward groups.
+BPZ intervals [0.10, 0.37, 0.55, 0.76, 1.20] include lower edges and exclude upper edges.
+Each gold group has in-bin spectra and passes the source-BPZ versus mean-spec-z threshold `max(0.4, 5 NMAD)` on the 600-tile reference sample.
+Groups empty in that reference skip the BPZ-versus-spec-z check.
+The source shear-weight sum divided by the number of in-bin calibrators weights each spectrum in its group.
+Linear deposition on [0, 5] with Δz=0.01 conserves mass and first moment under trapezoidal integration.
+Canonical z takes precedence; detected z2 fills missing canonical z after the offset in the recipe.

cosmo_val/nz/v2_bound9.yaml

@@ -0,0 +1,31 @@
+name: provisional_v2_bound9_DR6_BPZ_SOM_EGS_full_z_z2
+limitation: EGS-only provisional calibration for validation and null tests; no independent-field variance or spectroscopic-selection model.
+inputs:
+  # Fiducial spectroscopic Z_B and matched DR6 BPZ/ugriz; 11,114 EGS spectra.
+  calibrators: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/inputs/matched_calibration.npz
+  model: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/inputs/som_assignment_model.npz
+  gold_reference: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/inputs/gold_reference_600tiles.npz
+photometry:
+  selection: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/inputs/bound9_full_selection.parquet
+  dr6_root: /n17data/UNIONS/WL/photometry/UNIONS_DR6
+  output: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/cache/matched_full.npz
+som:
+  side: 101
+  epochs: 60
+  seed: 42
+groups: 2000
+edges: [0.10, 0.37, 0.55, 0.76, 1.20]
+gold_qc:
+  floor: 0.4
+  nmad_factor: 5
+# Canonical z takes precedence over z2.
+z2_offset: -0.1104888916015625
+grid: {min: 0, max: 5, step: 0.01}
+output: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/calibration
+diagnostics:
+  published_nz: /n17data/sguerrini/UNIONS/WL/nz/v1.4.6.3/nz_SP_v1.4.6.3_B.txt
+  tile_nz: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_first/products/nz_v2_bound9_tomo.txt
+  full_label: Full area, z + z2
+  tile_label: 600-tile reference
+  published_label: Published v1.4.6.3
+  source_label: 600-tile reference

scratch/cdaley/nz_som_reproduction/README.md

@@ -0,0 +1,13 @@
+# Full-area SOM n(z) reproduction
+
+Run inside the sp_validation container on an allocated candide compute node:
+
+```bash
+python scratch/cdaley/nz_som_reproduction/reproduce.py cosmo_val/nz/v2_bound9.yaml scratch/nz-reproduction --workers 16
+python scratch/cdaley/nz_som_reproduction/reproduce.py cosmo_val/nz/v2_bound9.yaml scratch/nz-rebuild --workers 16 --rebuild-model
+```
+
+Both commands compare the output against `reference.yaml` and assert zero `TOM_BIN_ID` differences across every source row and per-bin mean redshift differences below 1e-4.
+The rebuild command also checks exact equality of the SOM code, feature normalization, groups, edges and gold masks.
+`verify.py reference.yaml OUTPUT` can check an existing output independently.
+The small cosmo_val reader and SACC round-trip check is in `test_nz_calibration.py`.

scratch/cdaley/nz_som_reproduction/figures.py

@@ -0,0 +1,65 @@
+"""Render calibration diagnostics and saved-delivery z2 coverage."""
+
+import argparse
+import json
+from pathlib import Path
+
+import matplotlib
+
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import yaml
+
+from sp_validation.nz_diagnostics import COLORS, render
+
+matplotlib.rcParams.update({"text.usetex": False, "font.family": "DejaVu Sans"})
+
+parser = argparse.ArgumentParser(description=__doc__)
+parser.add_argument("recipe")
+parser.add_argument("output")
+parser.add_argument("figures")
+args = parser.parse_args()
+recipe = yaml.safe_load(Path(args.recipe).read_text())
+output = Path(args.output)
+figures = Path(args.figures)
+render(
+    dict(
+        recipe["diagnostics"],
+        calibration_arrays=output / "calibration_arrays.npz",
+        model=output / "som_assignment_model.npz",
+        summary=output / "summary.json",
+        source_catalogue=recipe["inputs"]["gold_reference"],
+        side=recipe["som"]["side"],
+        workers=16,
+    ),
+    figures,
+)
+# This comparison uses saved delivery counts; production runs only z+z2.
+delivery = Path("/n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/products")
+retention = json.loads((delivery / "summary.json").read_text())["retention"]
+stages = ["all", "complete_photometry", "bpz_bounds", "gold"]
+x = np.arange(len(stages))
+fig, ax = plt.subplots(figsize=(8, 4.8), layout="constrained")
+width = 0.34
+for j, (key, label, color) in enumerate(
+    (("zonly", "Canonical z only", COLORS[2]), ("both", "z + z2", COLORS[0]))
+):
+    bars = ax.bar(
+        x + (j - 0.5) * width,
+        [100 * retention[key][stage] for stage in stages],
+        width,
+        label=label,
+        color=color,
+    )
+    ax.bar_label(bars, fmt="%.1f%%", fontsize=8, padding=2)
+ax.set(
+    xticks=x,
+    xticklabels=["All", "Complete\nphotometry", "BPZ\nbounds", "Gold"],
+    ylabel="Retained bound-9 shear weight (%)",
+    ylim=(0, 110),
+    title="Coverage with and without z2 (saved delivery)",
+)
+ax.legend()
+fig.savefig(figures / "coverage_z2.png", dpi=100, bbox_inches="tight")
+plt.close(fig)

scratch/cdaley/nz_som_reproduction/reference.yaml

@@ -0,0 +1,3 @@
+reference:
+  nz: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/products/nz_v2_bound9_full_tomo.txt
+  tom_bins: /n17data/cdaley/unions/analyses/shear_tomographic/nz_v2_full/products/tom_bin_full.parquet

scratch/cdaley/nz_som_reproduction/reproduce.py

@@ -0,0 +1,33 @@
+"""Reproduce the full-area delivery, optionally rebuilding the assignment model."""
+
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+import yaml
+from verify import verify
+
+from sp_validation.nz_pipeline import run
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("recipe")
+    parser.add_argument("output")
+    parser.add_argument("--workers", type=int, default=16)
+    parser.add_argument("--rebuild-model", action="store_true")
+    args = parser.parse_args()
+    config = yaml.safe_load(Path(args.recipe).read_text())
+    run(config, args.output, workers=args.workers, rebuild=args.rebuild_model)
+    reference = yaml.safe_load(Path(__file__).with_name("reference.yaml").read_text())
+    verify(reference, args.output)
+    if args.rebuild_model:
+        new = np.load(Path(args.output) / "som_assignment_model.npz")
+        old = np.load(config["inputs"]["model"])
+        keys = ["code", "center", "scale", "group_map", "edges", "gold"]
+        result = {key: bool(np.array_equal(new[key], old[key])) for key in keys}
+        result["max_code_difference"] = float(abs(new["code"] - old["code"]).max())
+        result["gold_mismatches"] = int(np.count_nonzero(new["gold"] != old["gold"]))
+        (Path(args.output) / "rebuild.json").write_text(json.dumps(result, indent=2))
+        print(json.dumps(result, indent=2), flush=True)
+        assert all(result[key] for key in keys), result

scratch/cdaley/nz_som_reproduction/verify.py

@@ -0,0 +1,60 @@
+"""Check the frozen full-area delivery against a calibration output directory."""
+
+import argparse
+import json
+from pathlib import Path
+
+import numpy as np
+import pyarrow.parquet as pq
+import yaml
+
+
+def verify(config, output):
+    """Compare delivered full-area n(z) and every sidecar bin ID."""
+    output = Path(output)
+    reference = config["reference"]
+    actual = np.loadtxt(output / "nz_tomo.txt")
+    expected = np.loadtxt(reference["nz"])
+    z = actual[:, 0]
+    means = np.trapezoid(z[:, None] * actual[:, 1:], z, axis=0)
+    means /= np.trapezoid(actual[:, 1:], z, axis=0)
+    ref_mean = np.trapezoid(expected[:, :1] * expected[:, 1:], expected[:, 0], axis=0)
+    ref_mean /= np.trapezoid(expected[:, 1:], expected[:, 0], axis=0)
+    error = abs(means - ref_mean)
+    assert np.all(error < 1e-4), error
+    mismatches, count = 0, 0
+    # Both sidecars preserve the selected-catalogue row order.
+    a = pq.ParquetFile(output / "tom_bin.parquet")
+    b = pq.ParquetFile(reference["tom_bins"])
+    assert a.metadata.num_rows == b.metadata.num_rows
+    for aa, bb in zip(
+        a.iter_batches(columns=["TILE_ID", "NUMBER", "TOM_BIN_ID"]),
+        b.iter_batches(columns=["TILE_ID", "NUMBER", "TOM_BIN_ID"]),
+    ):
+        assert aa.column(0).equals(bb.column(0))
+        assert aa.column(1).equals(bb.column(1))
+        mismatches += np.count_nonzero(
+            aa.column(2).to_numpy() != bb.column(2).to_numpy()
+        )
+        count += len(aa)
+    assert count == a.metadata.num_rows, (count, a.metadata.num_rows)
+    assert mismatches == 0, mismatches
+    result = dict(
+        means=means.tolist(),
+        reference_means=ref_mean.tolist(),
+        absolute_mean_errors=error.tolist(),
+        rows=count,
+        tom_bin_mismatches=int(mismatches),
+        max_density_difference=float(abs(actual - expected).max()),
+    )
+    (output / "verification.json").write_text(json.dumps(result, indent=2) + "\n")
+    print(json.dumps(result, indent=2), flush=True)
+    return result
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("reference")
+    parser.add_argument("output")
+    args = parser.parse_args()
+    verify(yaml.safe_load(Path(args.reference).read_text()), args.output)

scripts/calibrate_nz.py

@@ -0,0 +1,6 @@
+"""Run a saved SOM n(z) recipe outside cosmo_val (see cosmo_val/nz/README.md)."""
+
+from sp_validation.nz_pipeline import main
+
+if __name__ == "__main__":
+    main()

scripts/prepare_nz_photometry.py

@@ -0,0 +1,24 @@
+#!/usr/bin/env python
+"""Prepare the DR6 positional photometry cache used by SOM n(z)."""
+
+import argparse
+from pathlib import Path
+
+import yaml
+
+from sp_validation.nz_photometry import prepare_nz_photometry
+
+
+def main():
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("recipe", help="YAML recipe containing a photometry section")
+    args = parser.parse_args()
+    photometry = yaml.safe_load(Path(args.recipe).read_text())["photometry"]
+    selection, dr6_root, output = (
+        photometry[key] for key in ("selection", "dr6_root", "output")
+    )
+    prepare_nz_photometry(selection, dr6_root, output)
+
+
+if __name__ == "__main__":
+    main()