From eb1b7095715927d825ff946d4b9bb54cfd9fa97b Mon Sep 17 00:00:00 2001 From: Cail Daley Date: Sat, 31 Oct 2026 01:32:39 +0100 Subject: [PATCH] Mask PSF stars and galaxies through one mask-config path A mask config's optional `maps` names a sky map per mask column (path or arc:/vos: URL, fetched once with vcp into ~/.cache/sp_validation/masks). get_masks_from_config reads a cut's column from the catalogue, or from the map at the objects' positions when the catalogue lacks it. catalogue_selection, the AND of those cuts, is what scripts/masking.py and the rho/tau star selection both evaluate. A cat_config psf block names its mask config (mask_config, relative to the catalogue config) and the cuts stars take (mask_cuts). ShapePipe's star catalogue has no MASK_* columns, so with mask_v2.0.yaml's MASK_8_Manual -> DR6 bit-8 map, stars inside the manual mask (M33's HII regions and clusters, shapepipe#993) are dropped under the same cut galaxies take from their column. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01JK7o9GrF3X5AAqQo1X5UmX --- config/calibration/mask_v2.0.yaml | 9 + scripts/masking.py | 23 ++- src/sp_validation/cosmo_val/core.py | 7 + .../cosmo_val/psf_systematics.py | 37 +++-- src/sp_validation/masks.py | 112 ++++++++++++- src/sp_validation/tests/test_masks.py | 154 +++++++++++++++++- 6 files changed, 317 insertions(+), 25 deletions(-) diff --git a/config/calibration/mask_v2.0.yaml b/config/calibration/mask_v2.0.yaml index ae8a7176..268d333b 100644 --- a/config/calibration/mask_v2.0.yaml +++ b/config/calibration/mask_v2.0.yaml @@ -32,6 +32,15 @@ params: sky_regions: False verbose: True +# Sky maps for mask columns a catalogue lacks: a cut on such a column reads +# the map (True = masked) at each object's position. Galaxy catalogues carry +# MASK_8_Manual; ShapePipe's PSF-star catalogue does not, so stars selected +# with this config (a cat_config psf block's mask_config/mask_cuts) take the +# manual mask from the DR6 bit-8 map. Values are paths or arc:/vos: URLs, +# fetched once into ~/.cache/sp_validation/masks (sp_validation.masks.map_path). +maps: + MASK_8_Manual: arc:home/mhudson/masks/mask_ugriz_nside131072_n8.hsp + # Masks ## Cuts on catalogue columns, in the v2 grammar (sp_validation.grammar) dat: diff --git a/scripts/masking.py b/scripts/masking.py index 20a4e890..78a3dd99 100644 --- a/scripts/masking.py +++ b/scripts/masking.py @@ -7,7 +7,7 @@ import yaml from sp_validation.io import Catalogue -from sp_validation.masks import apply_condition, catalogue_cuts +from sp_validation.masks import catalogue_cuts, catalogue_selection # ------------------------- # Spatially-structured cuts: these define the survey footprint. @@ -51,6 +51,8 @@ def apply_masks(data, mask_config, footprint_only=False): Dictionary parsed from the YAML mask configuration file. Expected structure: - mask_config["dat"] : list of cuts applied to `data` + - mask_config["maps"] : optional sky maps for cut columns + `data` lacks (see sp_validation.masks.MapColumns) - mask_config["metacal"] : derived-quantity parameters (e.g. relative size limits) @@ -69,17 +71,14 @@ def apply_masks(data, mask_config, footprint_only=False): False indicates the object is masked (removed). """ - # Initialize mask - mask = np.ones(len(data), dtype=bool) - - for cut in catalogue_cuts(mask_config): - col = cut["col_name"] - if footprint_only and col not in SPATIAL_CUTS: - continue - kind = cut["kind"] - value = cut["value"] - - mask &= apply_condition(data[col], kind, value, cut.get("ignore_bits", 0)) + cuts = None + if footprint_only: + cuts = [ + cut["col_name"] + for cut in catalogue_cuts(mask_config) + if cut["col_name"] in SPATIAL_CUTS + ] + mask = catalogue_selection(mask_config, data, cuts=cuts) # --- metacal relative size (skip for footprint-only) --- if not footprint_only: diff --git a/src/sp_validation/cosmo_val/core.py b/src/sp_validation/cosmo_val/core.py index 5e5576ca..7113fd51 100644 --- a/src/sp_validation/cosmo_val/core.py +++ b/src/sp_validation/cosmo_val/core.py @@ -406,6 +406,13 @@ def resolve_paths_for_version(ver): section["path"] = ( str(path) if path.is_absolute() else str(subdir / path) ) + # A psf block's mask config is a config file, not data: relative + # to the catalogue config, not to subdir. + psf = cc[ver].get("psf") + if isinstance(psf, dict) and psf.get("mask_config"): + psf["mask_config"] = str( + self.catalog_config_path.parent / psf["mask_config"] + ) resolve_paths_for_version("nz") processed = {"nz"} diff --git a/src/sp_validation/cosmo_val/psf_systematics.py b/src/sp_validation/cosmo_val/psf_systematics.py index ef0e004b..c649ae6f 100644 --- a/src/sp_validation/cosmo_val/psf_systematics.py +++ b/src/sp_validation/cosmo_val/psf_systematics.py @@ -12,6 +12,7 @@ import matplotlib.ticker as mticker import numpy as np import treecorr +import yaml from cs_util import plots as cs_plots from shear_psf_leakage import leakage from shear_psf_leakage import plots as psfleak_plots @@ -20,6 +21,7 @@ from .. import blinding, sacc_io from ..io import open_entry +from ..masks import catalogue_selection from ..rho_tau import ( get_rho_tau_w_cov, get_samples, @@ -358,16 +360,31 @@ def _get_galaxy_mask(self, ver, tomo_bin_id): return gal_mask def _get_star_mask(self, ver): - cat_star = open_entry(self.cc[ver]["psf"]) - PSF_flag = self.cc[ver]["psf"].get("PSF_flag") - star_flag = self.cc[ver]["psf"].get("star_flag") - if PSF_flag is not None: - if star_flag is not None: - star_mask = (cat_star[PSF_flag] == 0) & (cat_star[star_flag] == 0) - else: - star_mask = cat_star[PSF_flag] == 0 - else: - star_mask = np.ones(len(cat_star), dtype=bool) + """Return the stars that enter rho/tau for version ``ver``. + + A star is kept when the psf block's ``PSF_flag`` and ``star_flag`` + columns, where configured, are zero, and it passes the block's + ``mask_config``: a mask config file (relative to the catalogue + config) whose ``mask_cuts`` (column names; default all) select the + stars as they select galaxies. Star catalogues carry no mask columns, + so those cuts read the config's ``maps`` at ``ra_col``/``dec_col``. + """ + entry = self.cc[ver]["psf"] + cat_star = open_entry(entry) + star_mask = np.ones(len(cat_star), dtype=bool) + for key in ("PSF_flag", "star_flag"): + if entry.get(key) is not None: + star_mask &= cat_star[entry[key]] == 0 + if entry.get("mask_config"): + with open(entry["mask_config"]) as f: + mask_config = yaml.safe_load(f) + star_mask &= catalogue_selection( + mask_config, + cat_star, + cuts=entry.get("mask_cuts"), + ra_col=entry["ra_col"], + dec_col=entry["dec_col"], + ) return star_mask def set_params_rho_tau(self, ver, params, params_psf): diff --git a/src/sp_validation/masks.py b/src/sp_validation/masks.py index 6c3ca0f3..a48b532f 100644 --- a/src/sp_validation/masks.py +++ b/src/sp_validation/masks.py @@ -7,6 +7,11 @@ :Author: Martin Kilbinger """ +import functools +import os +import subprocess +from pathlib import Path + import healsparse as hsp import numpy as np from astropy.io import fits @@ -14,6 +19,8 @@ _KIND_ALIASES = {"smaller_equal": "less_equal"} +REMOTE_SCHEMES = ("arc:", "vos:") + _KIND_OPS = { "equal": lambda a, v: a == v, "not_equal": lambda a, v: a != v, @@ -64,6 +71,67 @@ def apply_condition(array, kind, value, ignore_bits=0): return _KIND_OPS[kind](array, value) +def map_path(source): + """Return a local path to the sky map ``source`` names. + + ``source`` is a path, or an ``arc:``/``vos:`` URL on CANFAR. A URL is + fetched once with ``vcp`` into the mask cache, ``$SP_VALIDATION_MASK_CACHE`` + or ``~/.cache/sp_validation/masks``, mirroring its remote path, and read + from there afterwards. + """ + source = str(source) + if not source.startswith(REMOTE_SCHEMES): + return source + cache = Path( + os.environ.get( + "SP_VALIDATION_MASK_CACHE", Path.home() / ".cache/sp_validation/masks" + ) + ) + local = cache / source.split(":", 1)[1].lstrip("/") + if not local.exists(): + local.parent.mkdir(parents=True, exist_ok=True) + part = local.with_name(f"{local.name}.part{os.getpid()}") + subprocess.run(["vcp", source, str(part)], check=True) + os.replace(part, local) + return str(local) + + +@functools.lru_cache(maxsize=4) +def read_map(source): + """Return the healsparse map ``source`` names (see :func:`map_path`).""" + return hsp.HealSparseMap.read(map_path(source)) + + +class MapColumns: + """A catalogue whose missing columns are read from sky maps. + + ``maps`` takes a column name to a boolean healsparse map (True = masked). + A column the catalogue has is read from the catalogue; a column it lacks + is the map's value at each object's ``ra_col``/``dec_col``. So one cut, + e.g. ``MASK_8_Manual equal False``, masks galaxies from their column and + stars, whose catalogue has no mask columns, from the map. + """ + + def __init__(self, cat, maps, ra_col="RA", dec_col="Dec"): + self._cat = cat + self._maps = maps + self._ra_col = ra_col + self._dec_col = dec_col + + def __len__(self): + return len(self._cat) + + def __getitem__(self, col): + if col not in self._maps: + return self._cat[col] + try: + return self._cat[col] + except (KeyError, ValueError): + ra, dec = (np.asarray(self._cat[c]) for c in (self._ra_col, self._dec_col)) + values = read_map(self._maps[col]).get_values_pos(ra, dec, lonlat=True) + return np.asarray(values, dtype=bool) + + def correlation_matrix(masks, confidence_level=0.9): n_key = len(masks) print(n_key) @@ -351,11 +419,15 @@ def catalogue_cuts(config): return config.get("dat") or [] -def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False): +def get_masks_from_config( + config, dat, masks_to_apply=None, verbose=False, ra_col="RA", dec_col="Dec" +): """Get Masks From Config. Return the masks of the config's ``dat`` cut list (``catalogue_cuts``), - evaluated on ``dat``. + evaluated on ``dat``. A cut on a column ``dat`` lacks reads it from the + config's ``maps`` (column name -> healsparse map path or ``arc:``/``vos:`` + URL), at ``dat``'s ``ra_col``/``dec_col`` (see :class:`MapColumns`). Parameters ---------- @@ -368,6 +440,8 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False): masks to apply exclusively; if `None` (default), use all masks verbose : bool, optional verbose output if ``True``; default is ``False`` + ra_col, dec_col : str, optional + position columns for map lookups; default ``RA``, ``Dec`` Returns ------- @@ -377,6 +451,9 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False): list of indices for given mask column name (label) """ + if config.get("maps"): + dat = MapColumns(dat, config["maps"], ra_col=ra_col, dec_col=dec_col) + # List to store all mask objects masks = [] @@ -402,3 +479,34 @@ def get_masks_from_config(config, dat, masks_to_apply=None, verbose=False): masks.append(my_mask) return masks, labels + + +def catalogue_selection(config, dat, cuts=None, ra_col="RA", dec_col="Dec"): + """Return the objects of ``dat`` that pass a mask config's cuts. + + The one place a mask config is evaluated on a catalogue, galaxies and PSF + stars alike: the AND of :func:`get_masks_from_config`'s masks. + + Parameters + ---------- + config : dict + mask config (``dat`` cut list, optional ``maps``) + dat : mapping + catalogue, indexable by column name + cuts : list of str, optional + column names of the cuts to apply; default all + ra_col, dec_col : str, optional + position columns for map lookups + + Returns + ------- + numpy.ndarray + boolean, True for objects kept + """ + masks, _ = get_masks_from_config( + config, dat, masks_to_apply=cuts, ra_col=ra_col, dec_col=dec_col + ) + keep = np.ones(len(dat), dtype=bool) + for my_mask in masks: + keep &= np.asarray(my_mask._mask, dtype=bool) + return keep diff --git a/src/sp_validation/tests/test_masks.py b/src/sp_validation/tests/test_masks.py index 60adc12a..4c25b02d 100644 --- a/src/sp_validation/tests/test_masks.py +++ b/src/sp_validation/tests/test_masks.py @@ -3,7 +3,9 @@ Covers ``sp_validation.masks.apply_condition`` — the single ``kind``/``value`` evaluator shared by the object-selection ``Mask`` class and the cosmo_inference footprint builder (issue #181) — plus an equivalence check -against ``Mask.apply()`` itself. +against ``Mask.apply()`` itself, and the one catalogue-masking path +(``catalogue_selection``) that galaxies and PSF stars share, with mask +columns a catalogue lacks read from the config's sky ``maps``. :Author: cdaley @@ -12,15 +14,22 @@ import importlib.util from pathlib import Path +import healpy as hp +import healsparse as hsp import numpy as np import numpy.testing as npt import pytest +import yaml +from astropy.table import Table +from sp_validation import masks from sp_validation.masks import ( Mask, apply_condition, catalogue_cuts, + catalogue_selection, get_masks_from_config, + map_path, ) pytestmark = pytest.mark.fast @@ -153,3 +162,146 @@ def test_masking_script_refuses_dat_ext(): with pytest.raises(ValueError, match="dat_ext"): masking.apply_masks(dat, config, footprint_only=True) assert "IMAFLAGS_ISO" in masking.SPATIAL_CUTS + + +_MANUAL_CUT = { + "col_name": "MASK_8_Manual", + "label": "manual mask", + "kind": "equal", + "value": False, +} + + +def _manual_map(path, *positions): + """Write a bool healsparse map, True at ``positions`` (ra, dec in deg).""" + nside = 1024 + sky_map = hsp.HealSparseMap.make_empty(32, nside, dtype=np.bool_) + for ra, dec in positions: + sky_map[hp.ang2pix(nside, ra, dec, lonlat=True, nest=True)] = True + sky_map.write(str(path)) + return str(path) + + +def _mask_config(tmp_path): + """A galaxy-style mask config: a FLAGS cut, the manual cut, a bit-8 map.""" + config = { + "maps": {"MASK_8_Manual": _manual_map(tmp_path / "manual.hsp", (10.1, 20.1))}, + "dat": [ + {"col_name": "FLAGS", "label": "SE", "kind": "less_equal", "value": 2}, + _MANUAL_CUT, + ], + } + (tmp_path / "mask.yaml").write_text(yaml.safe_dump(config)) + return config + + +def _star_entry(tmp_path, **mask): + stars = Table( + dict( + RA=[10.0, 10.1, 10.2, 10.3], + DEC=[20.0, 20.1, 20.2, 20.3], + HSM_FLAG_PSF=[0, 0, 0, 1], + HSM_FLAG_STAR=[0, 0, 0, 0], + ) + ) + tmp_path.mkdir(exist_ok=True) + stars.write(tmp_path / "stars.fits") + return dict( + path=str(tmp_path / "stars.fits"), + ra_col="RA", + dec_col="DEC", + PSF_flag="HSM_FLAG_PSF", + star_flag="HSM_FLAG_STAR", + **mask, + ) + + +def _star_mask(entry): + from sp_validation.cosmo_val.psf_systematics import PSFSystematicsMixin + + stub = PSFSystematicsMixin() + stub.cc = {"v": {"psf": entry}} + return stub._get_star_mask("v") + + +def test_star_inside_the_manual_map_is_dropped(tmp_path): + """A star catalogue has no MASK_8_Manual; the config's map supplies it.""" + _mask_config(tmp_path) + unmasked = _star_mask(_star_entry(tmp_path / "a")) + masked = _star_mask( + _star_entry( + tmp_path / "b", + mask_config=str(tmp_path / "mask.yaml"), + mask_cuts=["MASK_8_Manual"], + ) + ) + np.testing.assert_array_equal(unmasked, [True, True, True, False]) + np.testing.assert_array_equal(masked, [True, False, True, False]) + + +def test_galaxy_column_wins_over_the_map(tmp_path): + """A catalogue that has the mask column is cut on it, map or no map.""" + config = _mask_config(tmp_path) + galaxies = np.zeros( + 4, dtype=[("RA", "f8"), ("Dec", "f8"), ("FLAGS", "i2"), ("MASK_8_Manual", "?")] + ) + galaxies["RA"], galaxies["Dec"] = [10.0, 10.1, 10.2, 10.3], [20.0, 20.1, 20.2, 20.3] + galaxies["FLAGS"] = [0, 0, 4, 0] + galaxies["MASK_8_Manual"] = [False, False, False, True] + expected = [True, True, False, False] + np.testing.assert_array_equal(catalogue_selection(config, galaxies), expected) + del config["maps"] + np.testing.assert_array_equal(catalogue_selection(config, galaxies), expected) + + +def test_galaxy_selection_unchanged_on_the_v2_config(tmp_path): + """mask_v2.0.yaml's cuts select the same galaxies as evaluating each cut + directly, its bit-8 map unread because the catalogue has the column.""" + config_path = Path(__file__).resolve().parents[3] / "config/calibration" + config = yaml.safe_load((config_path / "mask_v2.0.yaml").read_text()) + rng = np.random.default_rng(1) + galaxies = { + cut["col_name"]: np.where(rng.random(200) < 0.9, cut["value"], not cut["value"]) + for cut in config["dat"] + if isinstance(cut["value"], bool) + } + galaxies.update( + FLAGS=rng.choice([0, 1, 16, 18, 4], 200), + N_EPOCH=rng.choice([1, 2, 3, 3], 200), + mag=rng.uniform(10, 35, 200), + NGMIX_MCAL_TYPES_FAIL=rng.choice([0, 0, 0, 1], 200), + NGMIX_N_EPOCH=rng.choice([0, 1, 2, 2], 200), + NGMIX_G1_PSF_ORIG_NOSHEAR=rng.choice([-10.0, 0.1, 0.1, 0.1], 200), + NGMIX_G2_PSF_ORIG_NOSHEAR=rng.choice([-10.0, 0.1, 0.1, 0.1], 200), + RA=rng.uniform(0, 1, 200), + Dec=rng.uniform(0, 1, 200), + ) + galaxies = Table(galaxies) + direct = np.ones(200, dtype=bool) + for cut in config["dat"]: + direct &= apply_condition( + galaxies[cut["col_name"]], + cut["kind"], + cut["value"], + cut.get("ignore_bits", 0), + ) + assert "MASK_8_Manual" in config["maps"] + assert 0 < direct.sum() < 200 + np.testing.assert_array_equal(catalogue_selection(config, galaxies), direct) + + +def test_remote_map_is_fetched_once_into_the_cache(tmp_path, monkeypatch): + calls = [] + + def fake_run(cmd, check): + calls.append(cmd) + Path(cmd[2]).write_text("map") + + monkeypatch.setenv("SP_VALIDATION_MASK_CACHE", str(tmp_path)) + monkeypatch.setattr(masks.subprocess, "run", fake_run) + url = "arc:home/someone/masks/m.hsp" + local = tmp_path / "home/someone/masks/m.hsp" + assert map_path(url) == str(local) == map_path(url) + assert local.read_text() == "map" + assert len(calls) == 1 and calls[0][:2] == ["vcp", url] + assert map_path("/a/local/map.hsp") == "/a/local/map.hsp"