cosmo_val dashboard: the catalogue paper's validation figures for v2

config/columns/shapepipe_v2.yaml

@@ -110,6 +110,17 @@ star:
     DEC: Declination [deg].
     HSM_G1_PSF: PSF model ellipticity g1 (HSM).
     HSM_G2_PSF: PSF model ellipticity g2 (HSM).
+    HSM_G1_STAR: Stellar ellipticity g1 (HSM).
+    HSM_G2_STAR: Stellar ellipticity g2 (HSM).
+    HSM_T_PSF: PSF model size T (HSM).
+    HSM_T_STAR: Stellar size T (HSM).
+    HSM_FLAG_PSF: HSM PSF measurement flag.
+    HSM_FLAG_STAR: HSM stellar measurement flag.
+    HSM_M4_1_PSF: PSF model fourth-moment component 1.
+    HSM_M4_2_PSF: PSF model fourth-moment component 2.
+    HSM_M4_1_STAR: Stellar fourth-moment component 1.
+    HSM_M4_2_STAR: Stellar fourth-moment component 2.
+    MAG: Star magnitude.
 
 shear:
   about: >-

scripts/dashboard/build_correlations.sbatch

@@ -0,0 +1,51 @@
+#!/bin/bash
+#SBATCH --job-name=sp-correlations
+#SBATCH --account=def-sfabbro
+#SBATCH --time=01:30:00
+#SBATCH --cpus-per-task=8
+#SBATCH --mem=128G
+#SBATCH --output=%x-%j.log
+#
+# Catalogue-paper position-shear and alpha(theta) products
+# (sp_validation.dashboard.correlations), one measurement per step so each
+# process holds only what it needs.
+#
+# Measured on the #425 snapshot (10,094 tiles, 73.0M galaxies, 5.76M valid-HSM
+# star detections), 8 cores; see MEASURED below. Pass the variables with
+# --export, e.g.
+#   sbatch --export=ALL,CATALOGUE_DIR=.../run-r425/catalogue,VERSION=SP_v2.0,\
+#     OUT_DIR=.../catpaper/products,GAIA=.../gaia_dr3_G20_ruwe1p4_unions20k.fits \
+#     scripts/dashboard/build_correlations.sbatch
+# Stage Gaia beforehand on a login node (compute nodes may lack internet):
+#   uv run scripts/dashboard/stage_gaia.py --tiles <tile lists> --out <stem>
+#
+# Measured at 10k tiles (#425, 8 cores; wall time incl. a 45 s catalogue read
+# per step; peak RSS of the Python process):
+#   stars + tiles   4 min    38 GB   (4.77M unique stars as lenses; tiles 6 s)
+#   lambda          2.5 min  15 GB   (7.3M galaxy lenses, 7.3M randoms)
+#   alpha           7 min    44 GB   (rho0 + tau0 for all and 3 magnitude bins)
+#   gaia            6 min    39 GB   (7.23M Gaia stars, 7.3M randoms)
+#   total          ~20 min
+# Extrapolated to 20.5k tiles (~2x galaxies, stars, Gaia stars and randoms;
+# TreeCorr NG/GG time scales roughly with the lens and source counts):
+# ~40-50 min and ~75-90 GB peak (alpha, stars), hence --mem=128G for headroom and 1.5 h.
+
+set -euo pipefail
+: "${CATALOGUE_DIR:?}" "${VERSION:?}" "${OUT_DIR:?}"
+SIF=${SIF:-/project/6001537/cdaley/containers/sp_validation-ddd58fe0.sif}
+SRC=${SRC:-$SLURM_SUBMIT_DIR/src}  # submit from the repository root
+module load apptainer
+
+run() {
+    apptainer exec --bind /scratch,/project,/home --env PYTHONPATH="$SRC" "$SIF" \
+        python -m sp_validation.dashboard.correlations \
+        --catalogue-dir "$CATALOGUE_DIR" --version "$VERSION" --out-dir "$OUT_DIR" \
+        --threads "$SLURM_CPUS_PER_TASK" "$@"
+}
+
+run --only stars tiles
+run --only lambda
+run --only alpha
+if [ -n "${GAIA:-}" ]; then
+    run --only gaia --gaia "$GAIA"
+fi

scripts/dashboard/build_dashboard.sh

@@ -0,0 +1,80 @@
+#!/bin/bash
+# Build every cosmo_val dashboard figure for one catalogue snapshot through Slurm.
+#
+# Submits the catalogue-level products (build_products.sbatch), the
+# position-shear products (build_correlations.sbatch) and the COSEBIs
+# scale-cut PTE grid (from the run's ξ± and covariance) in parallel, then a
+# render job that runs once all succeed. A grid already in OUT_DIR/products is
+# kept. With --after, all of it waits for that job
+# (e.g. the catalogue extraction) to finish successfully.
+#
+# Every cosmological product is checked against the blind its version is
+# declared under in the run's catalogue config: --cat-config, else RUN_DIR's
+# cat_config.yaml. With neither, the script refuses to start.
+#
+# Usage, from the repository root on a login node:
+#   scripts/dashboard/build_dashboard.sh CATALOGUE_DIR RUN_DIR VERSION OUT_DIR \
+#       [--gaia GAIA_FITS] [--after JOBID] [--cat-config CAT_CONFIG]
+# Products go to OUT_DIR/products, figures and figures.json to OUT_DIR/figures.
+# ACCOUNT (default def-sfabbro) and SPV_IMAGE (default the pinned SIF) are read
+# from the environment.
+set -euo pipefail
+catdir=$(realpath -m "${1:?catalogue directory}")
+rundir=$(realpath -m "${2:?cosmo_val run directory}")
+version=${3:?version}
+out=$(realpath -m "${4:?output directory}")
+shift 4
+gaia=""
+after=""
+catconfig=""
+while [ $# -gt 0 ]; do
+    case $1 in
+        --gaia) gaia=$(realpath "$2"); shift 2 ;;
+        --after) after=$2; shift 2 ;;
+        --cat-config) catconfig=$(realpath "$2"); shift 2 ;;
+        *) echo "unknown argument: $1" >&2; exit 2 ;;
+    esac
+done
+
+root=$(git rev-parse --show-toplevel)
+cd "$root"
+catconfig=${catconfig:-$rundir/cat_config.yaml}
+[ -e "$catconfig" ] || {
+    echo "no catalogue config at $catconfig: pass --cat-config with the config the run declared its blinds in" >&2
+    exit 2
+}
+echo "blinds declared by $catconfig"
+account=${ACCOUNT:-def-sfabbro}
+export SPV_IMAGE=${SPV_IMAGE:-/project/6001537/cdaley/containers/sp_validation-ddd58fe0.sif}
+mkdir -p "$out/products" "$out/figures" "$out/logs"
+dep=()
+[ -n "$after" ] && dep=(--dependency="afterok:$after")
+
+products=$(sbatch --parsable --account="$account" "${dep[@]}" \
+    --output="$out/logs/%x-%j.log" --export=ALL \
+    scripts/dashboard/build_products.sbatch "$catdir" "$version" "$out/products")
+correlations=$(sbatch --parsable --account="$account" "${dep[@]}" \
+    --output="$out/logs/%x-%j.log" \
+    --export=ALL,CATALOGUE_DIR="$catdir",VERSION="$version",OUT_DIR="$out/products",GAIA="$gaia",SIF="$SPV_IMAGE",SRC="$root/src" \
+    scripts/dashboard/build_correlations.sbatch)
+grid="$out/products/cosebis_grid_$version.npz"
+after_grid=""
+if [ ! -e "$grid" ]; then
+    # #425 full grid, 210 cuts: see the cosebis-grid log for timing.
+    cosebis=$(sbatch --parsable --account="$account" --job-name=dashboard-cosebis-grid \
+        --time=01:30:00 --cpus-per-task=16 --mem=32G \
+        --output="$out/logs/%x-%j.log" \
+        --wrap="module load apptainer && apptainer exec --bind /scratch,/project,/home \
+--env PYTHONPATH=$root/src $SPV_IMAGE python -m sp_validation.dashboard.cosebis_grid \
+--run-dir $rundir --version $version --cat-config $catconfig --out $grid")
+    after_grid=":$cosebis"
+fi
+render=$(sbatch --parsable --account="$account" --job-name=dashboard-render \
+    --dependency="afterok:$products:$correlations$after_grid" \
+    --time=00:30:00 --cpus-per-task=2 --mem=8G \
+    --output="$out/logs/%x-%j.log" \
+    --wrap="module load apptainer && apptainer exec --bind /scratch,/project,/home \
+--env PYTHONPATH=$root/src $SPV_IMAGE python -m sp_validation.dashboard.render \
+--run-dir $rundir --version $version --cat-config $catconfig --maps-dir $out/products \
+--correlations-dir $out/products --out $out/figures")
+echo "products $products  correlations $correlations  cosebis-grid ${cosebis:-kept}  render $render"

scripts/dashboard/build_products.sbatch

@@ -0,0 +1,38 @@
+#!/bin/bash
+#SBATCH --job-name=dashboard-products
+#SBATCH --cpus-per-task=4
+#SBATCH --mem=48G
+#SBATCH --time=04:00:00
+
+# Usage: sbatch build_products.sbatch CATALOGUE_DIR VERSION PRODUCT_DIR
+# SPV_IMAGE must name the sp_validation Apptainer image. Run from the checkout.
+# #425, 10,094 tiles: galaxy + star pass 3m42s, 7.45 GiB peak RSS
+# (Slurm step 25103890.71, 4 CPUs). Linear 20,500-tile extrapolation:
+# ~7m31s, ~15.1 GiB; quantile buffers scale with galaxy count.
+# Comprehensive pass: 59m15s, 2.99 GiB peak RSS (step 25103890.23,
+# 2 CPUs); linear 20,500-tile extrapolation ~2h00m20s, with streaming
+# memory approximately constant. Wall time is sensitive to filesystem I/O.
+set -euo pipefail
+catdir=$(realpath "${1:?catalogue directory}")
+version=${2:?version}
+out=$(realpath -m "${3:?product directory}")
+root=$(git rev-parse --show-toplevel)
+: "${SPV_IMAGE:?set SPV_IMAGE to the container image}"
+module load apptainer
+mkdir -p "$out"
+stage=$(mktemp -d "$out/.build.XXXXXXXX")
+trap 'rm -rf "$stage"' EXIT
+run() {
+    apptainer exec --bind /scratch,/project,/home \
+        --env "PYTHONPATH=$root/src" "$SPV_IMAGE" python "$@"
+}
+galaxies="$catdir/shape_catalog_cut_ngmix.fits"
+expected=$(run -c 'import fitsio,sys; print(fitsio.FITS(sys.argv[1])[1].get_nrows())' "$galaxies")
+run -m sp_validation.dashboard.field_maps \
+    --galaxies "$galaxies" --stars "$catdir/stars.hdf5" \
+    --version "$version" --out-dir "$stage"
+run -m sp_validation.dashboard.selection_hist \
+    --catalogue "$catdir/shape_catalog_comprehensive_ngmix.hdf5" \
+    --config "$catdir/config_mask.yaml" --expected "$expected" --version "$version" \
+    --out "$stage/selection_hist_${version}.npz"
+mv "$stage"/*.npz "$out/"

scripts/dashboard/stage_gaia.py

@@ -0,0 +1,149 @@
+# /// script
+# requires-python = ">=3.10"
+# dependencies = ["pyvo", "astropy", "healpy", "numpy"]
+# ///
+"""Stage Gaia DR3 stars (G < 20, RUWE < 1.4) over the UNIONS tiles.
+
+Runs on a login node (the query executes on the TAP server); compute nodes
+read the staged FITS. The region is every HEALPix level-7 (nside 128,
+nested) pixel touched by a tile, from one or more tile lists: a text file of
+``XXX.YYY`` names or a JSON manifest whose ``tiles`` mapping is keyed by name.
+Gaia ``source_id`` encodes the level-12 nested pixel in its high bits, so each
+level-7 pixel is a contiguous ``source_id`` range and the query is a union of
+indexed range scans, run synchronously in chunks. Each chunk is cached in ``<out>.parts/`` so a rerun resumes, and removed once the table is written. Writes ``<out>.fits`` and ``<out>.json`` (ADQL template,
+pixel list digest, row count).
+
+Usage::
+
+    uv run scripts/dashboard/stage_gaia.py --tiles tiles.txt MANIFEST.json \\
+        --out catpaper/aux/gaia_dr3_G20_ruwe1p4_unions20k
+"""
+
+import argparse
+import hashlib
+import json
+import shutil
+import time
+from pathlib import Path
+
+import healpy as hp
+import numpy as np
+import pyvo
+from astropy.table import Table, vstack
+
+TAP = "https://dc.zah.uni-heidelberg.de/tap"
+TABLE = "gaia.dr3lite"
+LEVEL = 7
+NSIDE = 2**LEVEL
+SOURCE_ID_PER_PIXEL = 2**35 * 4 ** (12 - LEVEL)
+TILE_HALF_DIAGONAL_DEG = 0.37  # 10000 px x 0.187"/px tiles, half-diagonal
+ADQL = (
+    "SELECT source_id, ra, dec, phot_g_mean_mag, ruwe FROM {table} "
+    "WHERE phot_g_mean_mag < 20 AND ruwe < 1.4 AND ({ranges})"
+)
+
+
+def read_tiles(path):
+    path = Path(path)
+    if path.suffix == ".json":
+        return list(json.loads(path.read_text())["tiles"])
+    return [line.strip() for line in path.read_text().split() if line.strip()]
+
+
+def tile_centre(name):
+    """CFIS tile centre (cs_util ``cfis.get_tile_coord_from_nixy``)."""
+    xi, yi = (int(part) for part in name.replace("-", ".").split("."))
+    dec = yi / 2 - 90
+    return xi / 2 / np.cos(np.deg2rad(dec)), dec
+
+
+def tile_pixels(names):
+    pixels = set()
+    radius = np.deg2rad(TILE_HALF_DIAGONAL_DEG)
+    for name in names:
+        ra, dec = tile_centre(name)
+        vec = hp.ang2vec(ra, dec, lonlat=True)
+        pixels.update(hp.query_disc(NSIDE, vec, radius, inclusive=True, nest=True))
+    return np.array(sorted(pixels))
+
+
+def ranges(pixels):
+    """Contiguous pixel runs as (first, last) pairs."""
+    breaks = np.flatnonzero(np.diff(pixels) != 1)
+    starts = np.r_[pixels[0], pixels[breaks + 1]]
+    ends = np.r_[pixels[breaks], pixels[-1]]
+    return list(zip(starts, ends))
+
+
+def where(runs):
+    return " OR ".join(
+        f"source_id BETWEEN {a * SOURCE_ID_PER_PIXEL} "
+        f"AND {(b + 1) * SOURCE_ID_PER_PIXEL - 1}"
+        for a, b in runs
+    )
+
+
+def fetch(service, pixels, attempts=4):
+    """Rows for ``pixels``; a chunk that keeps failing is fetched in halves."""
+    query = ADQL.format(table=TABLE, ranges=where(ranges(pixels)))
+    for attempt in range(attempts):
+        try:
+            return service.run_sync(query, maxrec=16_000_000).to_table()
+        except (pyvo.dal.DALServiceError, pyvo.dal.DALQueryError) as error:
+            print(f"retry {attempt + 1} ({len(pixels)} pixels): {error}", flush=True)
+            time.sleep(30)
+    if len(pixels) == 1:
+        raise RuntimeError(f"pixel {pixels[0]} failed after {attempts} attempts")
+    half = len(pixels) // 2
+    return vstack([fetch(service, pixels[:half]), fetch(service, pixels[half:])])
+
+
+def main(argv=None):
+    parser = argparse.ArgumentParser(description=__doc__.split("\n")[0])
+    parser.add_argument("--tiles", nargs="+", required=True)
+    parser.add_argument("--out", required=True, help="output path without suffix")
+    parser.add_argument("--pixels-per-query", type=int, default=300)
+    args = parser.parse_args(argv)
+
+    names = sorted({n for path in args.tiles for n in read_tiles(path)})
+    pixels = tile_pixels(names)
+    print(f"{len(names)} tiles -> {len(pixels)} level-{LEVEL} pixels", flush=True)
+    service = pyvo.dal.TAPService(TAP)
+    out = Path(args.out)
+    parts_dir = out.parent / f"{out.name}.parts"
+    parts_dir.mkdir(parents=True, exist_ok=True)
+    parts = []
+    for start in range(0, len(pixels), args.pixels_per_query):
+        chunk = pixels[start : start + args.pixels_per_query]
+        part = parts_dir / f"{chunk[0]:06d}-{chunk[-1]:06d}.fits"
+        if not part.exists():
+            rows = fetch(service, chunk)
+            rows.meta.clear()
+            rows.write(part)
+        table = Table.read(part)
+        print(f"pixels {start}-{start + len(chunk)}: {len(table)} rows", flush=True)
+        parts.append(table)
+    table = vstack(parts)
+    table.meta.clear()
+    out.parent.mkdir(parents=True, exist_ok=True)
+    table.write(out.with_suffix(".fits"), overwrite=True)
+    record = {
+        "service": TAP,
+        "table": TABLE,
+        "adql_template": ADQL.format(table=TABLE, ranges="<source_id ranges>"),
+        "healpix_level": LEVEL,
+        "nest": True,
+        "n_tiles": len(names),
+        "tile_lists": [str(Path(p).resolve()) for p in args.tiles],
+        "n_pixels": int(len(pixels)),
+        "pixels_sha256": hashlib.sha256(pixels.astype("<i8").tobytes()).hexdigest(),
+        "n_rows": int(len(table)),
+        "n_unique_source_id": int(len(np.unique(table["source_id"]))),
+    }
+    out.with_suffix(".json").write_text(json.dumps(record, indent=2) + "\n")
+    shutil.rmtree(parts_dir)
+    print(json.dumps(record, indent=2))
+
+
+if __name__ == "__main__":
+    main()