diff --git a/astra.yaml b/astra.yaml index 5e0e717a2..d93c2067c 100644 --- a/astra.yaml +++ b/astra.yaml @@ -1365,7 +1365,7 @@ analyses: psf_epoch_averaging, galaxy_pixel_weights, psf_likelihood_noise, megacam_ccd_flip, defect_fill, blend_handling, defect_weighting, epoch_masked_fraction_cut, - central_defect_veto] + central_defect_veto, defect_flagging] decisions: ngmix_seed_mode: label: Per-object RNG seeded from sky position @@ -1560,6 +1560,30 @@ analyses: options: megapipe_flip: label: Flip CCDs < 18 and 36/37 (MegaPipe orientation) + defect_flagging: + label: Adaptive exposure defects on tile objects + rationale: >- + [HARDCODED] We retain every instrument bit by OR'ing exposure flag + pixels on the object's coadd footprint (seg == obj_id), flipped to + each epoch, across epochs kept after the fraction cut and central + veto. NGMIX_IMAFLAGS_ISO invents no bit for unflagged zero-weight + or off-tile pixels. NGMIX_N_EPOCH_INTERP_ISO counts kept epochs + with an interpolated defect on that footprint, including those + unflagged pixels. Missing, empty, or stamp-edge-touching footprints + give -1 in both columns because completeness is unknown. Never-fit + rows retain the existing count convention of 0. The Napoli + collaboration decision requests adaptive exposure flags propagated + to tile objects. Stamp-wide diagnostics count a distant defect + equally for small and large objects. The coadd footprint avoids + exposure detection matching; it is not reprojected, consistent with + the existing segmentation overlay. These columns impose no cut + and do not alter the fits. + default: footprint_or_and_interp_epochs + options: + footprint_or_and_interp_epochs: + label: Footprint IMAFLAGS_ISO OR plus interpolated-epoch count + stamp_wide_only: + label: Stamp-wide interpolated count and minimum distances only defect_fill: label: Image content of defect pixels before metacal rationale: >- diff --git a/scripts/validation/masked_pixels/tile_defect_flags.md b/scripts/validation/masked_pixels/tile_defect_flags.md new file mode 100644 index 000000000..456cee404 --- /dev/null +++ b/scripts/validation/masked_pixels/tile_defect_flags.md @@ -0,0 +1,155 @@ +# Tile-object defect flags + +We propagate exposure instrument bits onto each tile object's coadd segmentation footprint, using only epochs that pass the masked-fraction cut and central veto. +The implementation does not change fits or impose a catalogue cut. +This change stacks on #925 (`feat/uberseg-seg-vignet`), which puts the int32 `SEG_VIGNET` column in the committed tile catalogue under every blend handling and makes `Tile_cat` read it. +Stacking only on #922 would leave the committed workflow without the footprint input and produce `-1` for fitted objects. + +- `NGMIX_IMAFLAGS_ISO` is the bitwise OR of exposure flag pixels on `tile_seg == obj_id`, OR'ed across kept epochs, with every instrument bit preserved. +- `NGMIX_N_EPOCH_INTERP_ISO` counts kept epochs containing at least one interpolated defect pixel on that footprint. +- Both are `-1` when segmentation is missing, the footprint is empty, or it touches a stamp edge. + We cannot establish completeness in those cases. + Never-fit final-catalogue rows retain the existing count convention of `0`. +- A zero-weight/off-tile pixel with flag `0` contributes no instrument bit, but contributes to the epoch count if interpolated. + Clean quarter-turn copies made by defect weighting do not count as defects. + +The coadd segmentation is MegaCam-flipped to each exposure's stamp orientation, without reprojection, matching the existing segmentation overlay. +The columns supplement the stamp-wide diagnostics rather than replacing the tile SExtractor's `IMAFLAGS_ISO`. + +## Fixture validation + +The synthetic path runs `prepare_postage_stamps` with multiple epochs and both flipped and unflipped CCDs. +It checks inside/outside-footprint defects, clean objects, unflagged zero-weight defects, rejected epochs, missing/empty/clipped footprints, safe float-to-integer conversion, integer FITS round trips in every shear extension, and final-catalogue values and never-fit sentinels. + +| Mutation | Observed failure | Result | +| --- | --- | --- | +| Reduce/count over the entire stamp | Outside-footprint case becomes `1028 / 2`, expected `0 / 0` | 2 assertion failures | +| Overwrite the OR on each epoch | Final clean epoch erases `1028`, leaving `0` | 2 assertion failures | +| Accumulate before either epoch cut | Rejected epochs contribute bit `64` and interpolation counts | 8 assertion failures | +| Drop the make_cat write mapping | Final-catalogue value stays `0`, expected `1028` | 1 assertion failure | +| Remove the workflow's SEG_VIGNET export | Committed detection disables the footprint | 6 assertion failures | +| Disable SEG_VIGNET in the committed ini | Committed detection disables the footprint | 6 assertion failures | +| Drop Tile_cat's SEG_VIGNET read | Ordinary clean objects receive `-1` instead of `0` | 6 assertion failures | + +All seven mutations were restored before running the suite. +`test_committed_detection_supplies_defect_flag_footprint` resolves the actual workflow and ini for data and image simulations under default, noisefill and uberseg handling, then checks the `SEG_VIGNET` load and a clean object's zero-valued diagnostics. +The runtime SIF lacks pytest; validation uses `shapepipe-develop-240b37e4.sif` on Slurm allocations under `def-mjhudson`. +The rebased branch's ngmix, make_cat, unit and PSF property suite passes 830 tests with one skipped; the workflow/DAG/params-pin suite passes 31 tests with one skipped. +The seven mutations were rerun on the rebased branch on allocation `24920131`, with source edits and pytest on the same compute node and bytecode caches cleared between mutations. +The params-pin test runs inside that image with the existing `/project/6001537/cdaley/snakemake-env/lib/python3.12/site-packages` prepended and pytest plugin autoload disabled, because the image lacks Snakemake. +Regeneration leaves the pin unchanged. +`completeness.py` counts files, not columns; this change adds no product file. + +## Flag dtype + +`split_exp.py` casts flag images to `int16`, and `get_stamps` preserves the input dtype. +On real exposure `2254826p.flag.fits.fz`, CCD HDU 1 is `uint8`, with values `0, 1, 3, 8, 10, 11`; applying the pipeline's split cast and stamp extraction produces `int16` stamps. +The real-tile check below additionally inspects all 112,734 kept flag stamps in newly regenerated SQLite stores: every stamp is `int16`. +The fixture also exercises integral `float32` flag stamps. +The new reduction validates finite, nonnegative, integral, int64-representable values before casting to `int64`, so floats cannot silently truncate or overflow. + +## Real-tile validation: 188.307 + +We regenerated inputs for tile **188.307** from the same campaign as 181.308. +It needs seven exposures rather than nine: `2608415`, `2254827`, `2608414`, `2254828`, `2254835`, `2114047`, and `2114046`. +All 21 raw exposure image/weight/flag files, the tile image and weight, and the external catalogue were already local. +Retrieval used symlinks; nothing was downloaded. +The committed detection exports `SEG_VIGNET`, and the check reads it through `Tile_cat`, without a script-only segmentation fallback. +We ran the real `prepare_postage_stamps` path with `blend_handling=uberseg` and `defect_fill=interpolate`, **without running galaxy fits, tile_ngmix, or make_cat**. +This diagnostic path stops before applying the fit-time defect-weighting scheme. +These are detections presented to epoch preparation, not a shear-selected sample. + +### Measured fractions + +Of **36,051 objects**, 374 have no kept epoch, 35,056 have a known footprint and kept epochs, and 621 have kept epochs but an edge-clipped footprint. +There are **112,734 kept epochs** out of 120,421 considered: the masked-fraction cut rejects 5,122 and the central veto rejects another 2,565. + +| Diagnostic | Count | Denominator | Fraction | +| --- | ---: | ---: | ---: | +| Interpolated defect on the footprint, `NGMIX_N_EPOCH_INTERP_ISO > 0` | 173 | 35,056 known footprints | **0.4935%** | +| Instrument bits on the footprint, `NGMIX_IMAFLAGS_ISO > 0` | 174 | 35,056 known footprints | **0.4963%** | +| Interpolation anywhere in a kept stamp | 12,943 | 35,677 objects with kept epochs | **36.2783%** | +| Unknown footprint, both new columns `-1` | 621 | 35,677 objects with kept epochs | **1.7406%** | +| Footprint touches a stamp edge, before epoch cuts | 633 | 36,051 detections | **1.7558%** | +| Empty footprint | 0 | 36,051 detections | **0%** | + +Relative to all detections, the first two rates are 0.4799% and 0.4826%, and the stamp-wide rate is 35.9019%. +All 173 interpolation-positive objects also have instrument bits; one additional object has instrument bits without footprint interpolation. +The diagnostics distinguish flagged pixels from the subset actually interpolated, so those counts need not be identical. +The 1.74% edge-sentinel share is below the proposed few-percent threshold for reconsidering the truncation policy; this tile does not justify changing it. +The script leaves no-kept-epoch array entries at `-1` to distinguish them from measured rows, whereas final-catalogue never-fit rows use `0`. +After that final-catalogue convention, a literal `NGMIX_IMAFLAGS_ISO != 0` would select **795/36,051 (2.2052%)**, including the 621 unknowns; **174** is the count with positive instrument bits. +Consumers must handle the unknown sentinel and fit status explicitly. + +| Instrument bit | Objects carrying the bit | +| --- | ---: | +| `1` | 82 | +| `2` | 24 | +| `8` | 96 | + +These counts overlap. +The exact positive OR values are `1: 55`, `3: 23`, `8: 92`, `9: 3`, and `11: 1` objects. +No other instrument bits occur on known kept footprints. + +### Poisson comparison + +The requested independent-pixel reference is the object-averaged `1 - exp(-0.0028 * A * N)`, using footprint area `A` and epoch count `N` for each known footprint. +Here the mean area is **66.77 pixels**, the median is **31 pixels**, and the mean number of kept epochs is **3.160**. +Using each object's values rather than plugging their means into the exponential gives **31.75%**, compared with the measured **0.496%** bit-bearing rate. + +The central veto protects compact footprints: **78.22%** of known footprints lie entirely inside its radius, and only **19.12 pixels** lie outside it on average. +Replacing `A` by that outside-veto area lowers the reference to **6.81%**, still well above the measured rate. +The veto is defined on interpolable defects, so this remains an explanatory reference rather than an exact conditional model for instrument bits. +For the same known objects with kept epochs, **6.05%** touch an instrument flag in any candidate epoch before either cut; the corresponding unconditioned reference is **33.00%**. +Thus the disagreement already exists before the cuts, and epoch rejection reduces the observed fraction further. + +The 0.28% rate from #887 is an exposure-level DATASEC average, not a local independent probability calibrated for this tile's detected objects. +Bad columns and cosmic rays are spatially correlated, detection is not random with respect to defects, and conditioning on the epoch cuts changes the distribution again. +In the kept stamps here, **0.1448%** of all pixels carry instrument flags, while only **0.01491%** of known-footprint pixels do. +Those measurements and the before/after-cut comparison explain why treating every footprint pixel as an independent 0.28% trial overpredicts the object rate. +We have not measured a separate full-DATASEC rate for these seven exposures, so the relative contributions of local defect density, correlations and detection selection are not isolated. + +### Spatial check + +![Kept-epoch footprint defects on tile 188.307 over the exposure defect map](tile_defect_flags_188307.png) + +Orange points are the 173 interpolation-positive objects; grey points are other detections. +The blue map is `/project/6001537/cdaley/sp-products/smk-g17-int/nflagged_smk-g17-int.hsp`, sampled on the tile at NSIDE 131072. +Its value is **flagged CCD pixels per sky pixel**, summed over exposures, not a count of flagged exposures. +The map shows the expected thin bad-column and other structured defects. +Only **23/173** affected-object centroids land directly in a positive map cell, because the flagged pixel can be displaced from the centroid within the footprint. +Every affected centroid (**173/173**) is within 8 arcsec of a positive map cell, compared with **1,644/2,000 (82.2%)** deterministic sampled controls. +This is a coarse spatial consistency check, not an independent false-positive estimate: the map combines campaign exposures and its any-touch rasterization broadens thin defects. +The grid fraction with any positive map value is 6.97%; the corresponding fraction at other object centroids is 7.36%. + +### Reproduction and resource scope + +The private run root is `/scratch/cdaley/defectflag-check/run/188307/`. +`config.yaml`, `dryrun.log` and `path-audit.txt` record a **29-job** local Snakemake plan: seven retrievals, seven splits, seven exposure PSF chains and eight tile jobs. +All outputs were audited below `/scratch/cdaley/defectflag-check/` before execution. +We mounted `/project`, `/scratch` and `/home` read-only, then overlaid only that private prefix writable. +Neither production products/run directories nor `feeder/` or `latest/` were modified. + +The coadd-weight CLASSIC vignet FITS copy was too slow on shared scratch and was stopped; it is not an input to the diagnostic path. +`EpochSnakefile`, `epoch_vignets.ini`, `epoch-dryrun.log` and `epoch-run.log` record the private one-rule restart of **only the committed MULTI-EPOCH vignetmaker section**, using the completed detection, merged headers, exposure forest and PSF store. +It completed with zero runner errors and **36,051 rows in each of the five epoch stores**. +We did not mark the original full tile-vignet stage complete. +The private completion record verifies the actual five store row counts after the writer exits successfully. + +`measure.sh` and `measure-cached.log` record the no-fitting diagnostic command. +`--cache-stores` reads those SQLite stores into RAM without modifying them; exposure WCS objects are decoded once rather than repeatedly per object. +A synthetic end-to-end check verifies identical per-object arrays with and without this cache. +The final targeted pytest rerun was interrupted by the allocation time limit, so it adds no completed-suite result beyond the earlier 830-test and 31-test runs above. +The real measurement, spatial plot and cache-equivalence smoke had already completed. + +The single `def-mjhudson` allocation was **24921492**, with 16 CPUs, 40 GiB and a two-hour limit; it has ended, and we also issued `scancel`. +`diskusage_report` showed 398/1024 GiB scratch before and 369/1024 GiB afterward; these are shared, non-transactional quota snapshots rather than the private run's incremental size. +The requested runtime lacks pytest, so we used `/project/6001537/cdaley/containers/shapepipe-develop-240b37e4.sif`. + +The measured per-object arrays, original JSON and PNG are `/scratch/cdaley/defectflag-check/188.307.{npz,json,png}`. +The committed [JSON summary](tile_defect_flags_188307.json) preserves the full input paths, counts, denominators and spatial statistics. +Run `bash /scratch/cdaley/defectflag-check/run/188307/measure-exec.sh` in an equivalent allocation to repeat the measurement after changing its job ID, or use the command in `tile_defect_flags.py` with those paths. +`--plot-only --out /scratch/cdaley/defectflag-check/188.307 --nflagged --plot ` redraws the figure from existing arrays and refreshes the spatial statistics without rereading the stores. + +Back-filling would require retained or regenerated stores, but would not require shear refits. +No back-fill has been attempted. diff --git a/scripts/validation/masked_pixels/tile_defect_flags.py b/scripts/validation/masked_pixels/tile_defect_flags.py new file mode 100644 index 000000000..bc76583dc --- /dev/null +++ b/scripts/validation/masked_pixels/tile_defect_flags.py @@ -0,0 +1,301 @@ +"""Measure footprint defect flags without fitting, using retained tile stores. + +Run inside the ShapePipe container on a compute allocation, with PYTHONPATH=src:: + + python scripts/validation/masked_pixels/tile_defect_flags.py \ + --sexcat SEXCAT --image IMAGE_SQLITE \ + --flag FLAG_SQLITE --weight WEIGHT_SQLITE --psf GALAXY_PSF_SQLITE \ + --headers HEADERS_SQLITE --background BACKGROUND_SQLITE \ + --rms BACKGROUND_RMS_SQLITE --blend-handling uberseg \ + --nflagged NFLAGGED.hsp --out /scratch/cdaley/defectflag-check/TILE \ + --plot scripts/validation/masked_pixels/tile_defect_flags-TILE.png + +All stores open read-only; nothing is written beneath an input directory. +SEXCAT must carry the committed workflow's SEG_VIGNET column. Omit +--background only for already background-subtracted images; omit --rms when the run +used stamp-estimated RMS. Pass the production blend/fill choices explicitly. +The output contains every tile row, including unavailable PSFs, zero-kept- +epoch objects and unknown footprints, distinguished by status and -1. +The Poisson expectation uses 0.0028 as the exposure DATASEC pixel flag rate +from #887. It is an unconditioned reference, not a calibration: spatially +correlated defects and the central veto violate its independence assumption. +After measuring, --plot-only --out PREFIX --nflagged MAP --plot PNG regenerates +just the figure and spatial comparison from PREFIX.npz and PREFIX.json. +""" + +import argparse +from collections import Counter +from contextlib import ExitStack +import json +from pathlib import Path +from types import SimpleNamespace + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np +from sqlitedict import SqliteDict + +from shapepipe.modules.ngmix_package.ngmix import ( + EPOCH_INTERPOLATED_DEFECT_RADIUS, Ngmix, Tile_cat, centre_distance, + prepare_postage_stamps, +) + + +def measure(args): + """Read production stamps and invoke the pipeline's epoch preparation.""" + tile = Tile_cat(args.sexcat) + if tile.seg is None: + raise ValueError("A real-tile validation requires segmentation stamps") + n = len(tile.obj_id) + data = dict(number=tile.obj_id, ra=tile.ra, dec=tile.dec) + for key in ("imaflags_iso", "n_epoch_interp_iso", "n_epoch_interp", + "n_kept", "area", "area_outside_veto", "n_candidate", + "precut_imaflags_iso"): + data[key] = np.full(n, -1, dtype=np.int64) + data["footprint_clipped"] = np.zeros(n, dtype=bool) + status, dtypes, cuts = Counter(), Counter(), Counter() + stamp_pixels = stamp_flagged = footprint_pixels = footprint_flagged = 0 + with ExitStack() as stack: + def store(path): + if path is None: + return None + if not Path(path).is_file(): + raise FileNotFoundError(path) + db = SqliteDict(path, flag="r") + if getattr(args, "cache_stores", False): + print(f"Caching {Path(path).name}", flush=True) + cached = dict(db.items()) + db.close() + print(f"Cached {len(cached)} rows", flush=True) + return cached + stack.callback(db.close) + return db + + vignet = SimpleNamespace( + gal_vign_cat=store(args.image), flag_vign_cat=store(args.flag), + weight_vign_cat=store(args.weight), psf_vign_cat=store(args.psf), + f_wcs_file=store(args.headers), bkg_vign_cat=store(args.background), + bkg_rms_vign_cat=store(args.rms), + ) + # Decode exposure-level WCS objects once rather than for every object. + vignet.f_wcs_file = dict(vignet.f_wcs_file.items()) + radius = centre_distance(tile.seg[0].shape) + for i, obj_id in enumerate(tile.obj_id): + if i % 1000 == 0: + print(f"{i}/{n}", flush=True) + footprint = tile.seg[i] == obj_id + data["area"][i] = footprint.sum() + data["footprint_clipped"][i] = ( + footprint[[0, -1], :].any() or footprint[:, [0, -1]].any() + ) + data["area_outside_veto"][i] = ( + footprint & (radius >= EPOCH_INTERPOLATED_DEFECT_RADIUS) + ).sum() + key = str(obj_id) + if key not in vignet.psf_vign_cat or key not in vignet.gal_vign_cat: + status["missing_store_row"] += 1 + continue + psf, gal = vignet.psf_vign_cat[key], vignet.gal_vign_cat[key] + if isinstance(psf, str) or isinstance(gal, str): + status["empty_store_row"] += 1 + continue + flags = vignet.flag_vign_cat[key] + data["n_candidate"][i] = len(psf) + if footprint.any() and not data["footprint_clipped"][i]: + bits = 0 + for name in psf: + own = Ngmix.MegaCamFlip(footprint, int(name.split("-")[1])) + values = np.asarray(flags[name]["VIGNET"])[own] + bits |= int(np.bitwise_or.reduce(values.astype(np.int64))) + data["precut_imaflags_iso"][i] = bits + # Reuse the fetched flag row when the production preparation reads + # it; no second SQLite fetch is needed for the pre-cut comparison. + per_object = SimpleNamespace(**(vars(vignet) | { + "flag_vign_cat": {key: flags}, + })) + stamp = prepare_postage_stamps( + per_object, obj_id, i, tile, bkg_sub=args.background is not None, + psf_obj=psf, gal_obj=gal, blend_handling=args.blend_handling, + defect_fill=args.defect_fill, + ) + data["n_kept"][i] = len(stamp.flags) + cuts.update(stamp.epoch_cuts) + for flag, epoch_seg in zip(stamp.flags, stamp.segs): + dtypes[str(flag.dtype)] += 1 + flagged_pixels = flag != 0 + stamp_pixels += flag.size + stamp_flagged += int(flagged_pixels.sum()) + if stamp.imaflags_iso >= 0: + own = epoch_seg == obj_id + footprint_pixels += int(own.sum()) + footprint_flagged += int((flagged_pixels & own).sum()) + if not stamp.flags: + status["no_kept_epoch"] += 1 + continue + for name, value in stamp.defect_diagnostics().items(): + if name in data: + data[name][i] = value + status["unknown_footprint" if stamp.imaflags_iso < 0 else "known"] += 1 + valid = data["imaflags_iso"] >= 0 + kept = data["n_kept"] > 0 + interp = data["n_epoch_interp_iso"] > 0 + flagged = data["imaflags_iso"] > 0 + wide = data["n_epoch_interp"] > 0 + summary = dict( + inputs=vars(args), n_objects=n, status=dict(status), + kept_flag_vignet_dtypes=dict(dtypes), epoch_cuts=dict(cuts), + kept_stamp_flagged_pixel_fraction=( + stamp_flagged / stamp_pixels if stamp_pixels else None + ), + kept_known_footprint_flagged_pixel_fraction=( + footprint_flagged / footprint_pixels if footprint_pixels else None + ), + footprints=dict( + n_edge_touching=int(data["footprint_clipped"].sum()), + fraction_edge_touching=float(data["footprint_clipped"].mean()), + n_empty=int((data["area"] == 0).sum()), + n_edge_touching_with_kept_epoch=int( + (data["footprint_clipped"] & kept).sum() + ), + n_unknown_with_kept_epoch=int((~valid & kept).sum()), + ), + n_interp_iso=int(interp.sum()), n_imaflags_iso=int(flagged.sum()), + n_stamp_interp=int(wide.sum()), + fractions_all_rows=dict(interp_iso=float(interp.mean()), + imaflags_iso=float(flagged.mean()), + stamp_interp=float(wide.mean())), + bits={str(1 << bit): int(((data["imaflags_iso"][valid] & (1 << bit)) != 0).sum()) + for bit in range(63) + if ((data["imaflags_iso"][valid] & (1 << bit)) != 0).any()}, + ) + if valid.any(): + summary["fractions_known_footprint"] = dict( + interp_iso=float(interp[valid].mean()), + imaflags_iso=float(flagged[valid].mean()), + ) + summary["fraction_precut_instrument_flags_known_kept"] = float( + (data["precut_imaflags_iso"][valid] > 0).mean() + ) + summary["poisson_before_cuts_known_kept"] = float(np.mean( + -np.expm1(-0.0028 * data["area"][valid] * data["n_candidate"][valid]) + )) + for area in ("area", "area_outside_veto"): + summary[f"poisson_{area}"] = float(np.mean( + -np.expm1(-0.0028 * data[area][valid] * data["n_kept"][valid]) + )) + summary["mean_area_known"] = float(data["area"][valid].mean()) + summary["median_area_known"] = float(np.median(data["area"][valid])) + summary["mean_area_outside_veto_known"] = float( + data["area_outside_veto"][valid].mean() + ) + summary["fraction_known_footprints_inside_veto"] = float( + (data["area_outside_veto"][valid] == 0).mean() + ) + summary["mean_kept_epochs_known"] = float(data["n_kept"][valid].mean()) + if kept.any(): + summary["fraction_stamp_interp_kept"] = float(wide[kept].mean()) + summary["fraction_unknown_footprint_kept"] = float((~valid)[kept].mean()) + out = Path(args.out) + out.parent.mkdir(parents=True, exist_ok=True) + np.savez_compressed(str(out) + ".npz", **data) + Path(str(out) + ".json").write_text(json.dumps(summary, indent=2) + "\n") + if args.plot: + plot(args, data, summary) + print(json.dumps(summary, indent=2)) + + +def plot(args, data, summary): + """Overlay interpolated-footprint objects on the retained nflagged map.""" + import healsparse + import hpgeom + from matplotlib.colors import LogNorm + from figstyle import INK2, ORANGE, PAPER + + if not args.nflagged: + raise ValueError("--plot requires the exposure --nflagged map") + sky = healsparse.HealSparseMap.read(args.nflagged) + ra, dec = data["ra"], data["dec"] + # Unwrap a tile crossing RA=0 before laying out its sky grid. + centre = float(np.rad2deg(np.angle(np.mean(np.exp(1j * np.deg2rad(ra)))))) + x = (ra - centre + 180) % 360 - 180 + extent = [x.min(), x.max(), dec.min(), dec.max()] + xx, yy = np.meshgrid(np.linspace(*extent[:2], 700), + np.linspace(*extent[2:], 700)) + image = sky.get_values_pos((xx + centre) % 360, yy, lonlat=True).astype(float) + image[image < 0] = np.nan + at_objects = sky.get_values_pos(ra, dec, lonlat=True) + selected = data["n_epoch_interp_iso"] > 0 + summary["map_nflagged_positive_at_interp_iso"] = int( + ((at_objects > 0) & selected).sum() + ) + summary["map_grid_fraction_positive"] = float((image > 0).mean()) + summary["map_fraction_positive_at_other_centroids"] = float( + (at_objects[~selected] > 0).mean() + ) + # A flagged footprint pixel need not share the centroid's 1.6-arcsec + # HealSparse cell. Eight arcsec covers a 51-px stamp's half-diagonal + # (6.6 arcsec) plus one map-cell half-diagonal. + def nearby(rows): + return sum(bool((sky.get_values_pix(hpgeom.query_circle( + sky.nside_sparse, ra[i], dec[i], 8 / 3600, inclusive=True, + )) > 0).any()) for i in rows) + + selected_rows = np.flatnonzero(selected) + controls = np.flatnonzero(~selected & (data["imaflags_iso"] >= 0)) + rng = np.random.default_rng(0) + controls = rng.choice(controls, size=min(len(controls), 2000), replace=False) + summary["map_positive_within_8_arcsec"] = dict( + flagged_n=len(selected_rows), flagged_matches=nearby(selected_rows), + control_n=len(controls), control_matches=nearby(controls), + ) + fig, ax = plt.subplots(figsize=(9, 8), facecolor=PAPER) + ax.set_facecolor(PAPER) + im = ax.imshow(image, origin="lower", extent=extent, cmap="Blues", + norm=LogNorm(vmin=1, vmax=max(2, np.nanmax(image))), + interpolation="nearest", aspect="auto") + ax.scatter(x[~selected], dec[~selected], s=0.6, c=INK2, alpha=.12, + linewidths=0, label="Other tile detections") + ax.scatter(x[selected], dec[selected], s=13, c=ORANGE, linewidths=.3, + edgecolors="white", label=f"Interpolated defect on footprint ({selected.sum():,})") + ax.set(xlabel=f"RA offset from {centre % 360:.3f}° [deg]", + ylabel="Declination [deg]", + title=f"{Path(args.out).name}: kept-epoch defects on object footprints") + ax.invert_xaxis() + ax.legend(loc="upper right", markerscale=2) + fig.colorbar(im, ax=ax, label="Flagged CCD pixels per sky pixel (nflagged)") + fig.tight_layout() + Path(args.plot).parent.mkdir(parents=True, exist_ok=True) + fig.savefig(args.plot, dpi=180) + plt.close(fig) + Path(args.out + ".json").write_text(json.dumps(summary, indent=2) + "\n") + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + for name in ("sexcat", "image", "flag", "weight", "psf", "headers", "out"): + parser.add_argument("--" + name, required=(name == "out")) + for name in ("background", "rms", "nflagged", "plot"): + parser.add_argument("--" + name) + parser.add_argument("--blend-handling", choices=("noisefill", "uberseg")) + parser.add_argument("--defect-fill", choices=("interpolate", "noise"), default="interpolate") + parser.add_argument("--plot-only", action="store_true") + parser.add_argument("--cache-stores", action="store_true", + help="Cache read-only stores in RAM before measuring") + args = parser.parse_args() + if args.plot_only: + if not args.plot or not args.nflagged: + parser.error("--plot-only requires --plot and --nflagged") + summary = json.loads(Path(args.out + ".json").read_text()) + with np.load(args.out + ".npz") as data: + plot(args, data, summary) + else: + for name in ("sexcat", "image", "flag", "weight", "psf", "headers", + "blend_handling"): + if getattr(args, name) is None: + parser.error(f"--{name.replace('_', '-')} is required for measurement") + measure(args) + + +if __name__ == "__main__": + main() diff --git a/scripts/validation/masked_pixels/tile_defect_flags_188307.json b/scripts/validation/masked_pixels/tile_defect_flags_188307.json new file mode 100644 index 000000000..3d7288aef --- /dev/null +++ b/scripts/validation/masked_pixels/tile_defect_flags_188307.json @@ -0,0 +1,79 @@ +{ + "inputs": { + "sexcat": "/scratch/cdaley/defectflag-check/run/188307/work/tiles/18/188.307/output/run_sp_tile_Sx/sextractor_runner/output/sexcat-188-307.fits", + "image": "/scratch/cdaley/defectflag-check/run/188307/epoch-output/run_sp_epoch_vignets/vignetmaker_runner/output/image_vignet-188-307.sqlite", + "flag": "/scratch/cdaley/defectflag-check/run/188307/epoch-output/run_sp_epoch_vignets/vignetmaker_runner/output/flag_vignet-188-307.sqlite", + "weight": "/scratch/cdaley/defectflag-check/run/188307/epoch-output/run_sp_epoch_vignets/vignetmaker_runner/output/weight_vignet-188-307.sqlite", + "psf": "/scratch/cdaley/defectflag-check/run/188307/tile-store/sp-c9b3bf5e-188.307/output/run_sp_tile_PiViVi/psfex_interp_runner/output/galaxy_psf-188-307.sqlite", + "headers": "/scratch/cdaley/defectflag-check/run/188307/tile-store/sp-c9b3bf5e-188.307/wcs/log_exp_headers-188-307.sqlite", + "out": "/scratch/cdaley/defectflag-check/188.307", + "background": "/scratch/cdaley/defectflag-check/run/188307/epoch-output/run_sp_epoch_vignets/vignetmaker_runner/output/background_vignet-188-307.sqlite", + "rms": "/scratch/cdaley/defectflag-check/run/188307/epoch-output/run_sp_epoch_vignets/vignetmaker_runner/output/background_rms_vignet-188-307.sqlite", + "nflagged": "/project/6001537/cdaley/sp-products/smk-g17-int/nflagged_smk-g17-int.hsp", + "plot": "/scratch/cdaley/defectflag-check/188.307.png", + "blend_handling": "uberseg", + "defect_fill": "interpolate", + "plot_only": false, + "cache_stores": true + }, + "n_objects": 36051, + "status": { + "no_kept_epoch": 374, + "known": 35056, + "unknown_footprint": 621 + }, + "kept_flag_vignet_dtypes": { + "int16": 112734 + }, + "epoch_cuts": { + "considered": 120421, + "masked_fraction": 5122, + "central_veto": 2565 + }, + "kept_stamp_flagged_pixel_fraction": 0.0014477094273839075, + "kept_known_footprint_flagged_pixel_fraction": 0.00014910042069069652, + "footprints": { + "n_edge_touching": 633, + "fraction_edge_touching": 0.017558458849962554, + "n_empty": 0, + "n_edge_touching_with_kept_epoch": 621, + "n_unknown_with_kept_epoch": 621 + }, + "n_interp_iso": 173, + "n_imaflags_iso": 174, + "n_stamp_interp": 12943, + "fractions_all_rows": { + "interp_iso": 0.004798757316024521, + "imaflags_iso": 0.004826495797620039, + "stamp_interp": 0.3590191672907825 + }, + "bits": { + "1": 82, + "2": 24, + "8": 96 + }, + "fractions_known_footprint": { + "interp_iso": 0.004934961204929256, + "imaflags_iso": 0.004963486992240986 + }, + "fraction_precut_instrument_flags_known_kept": 0.06053172067549065, + "poisson_before_cuts_known_kept": 0.33003992762146317, + "poisson_area": 0.3174972118531098, + "poisson_area_outside_veto": 0.06808699656730423, + "mean_area_known": 66.76725810132359, + "median_area_known": 31.0, + "mean_area_outside_veto_known": 19.11758329529895, + "fraction_known_footprints_inside_veto": 0.7822056138749429, + "mean_kept_epochs_known": 3.1601152441807394, + "fraction_stamp_interp_kept": 0.36278274518597414, + "fraction_unknown_footprint_kept": 0.017406172043613534, + "map_nflagged_positive_at_interp_iso": 23, + "map_grid_fraction_positive": 0.06971020408163266, + "map_fraction_positive_at_other_centroids": 0.07358269691733095, + "map_positive_within_8_arcsec": { + "flagged_n": 173, + "flagged_matches": 173, + "control_n": 2000, + "control_matches": 1644 + } +} diff --git a/scripts/validation/masked_pixels/tile_defect_flags_188307.png b/scripts/validation/masked_pixels/tile_defect_flags_188307.png new file mode 100644 index 000000000..f70042ce1 Binary files /dev/null and b/scripts/validation/masked_pixels/tile_defect_flags_188307.png differ diff --git a/src/shapepipe/modules/make_cat_package/make_cat.py b/src/shapepipe/modules/make_cat_package/make_cat.py index b39dd4f00..a04d4c52c 100644 --- a/src/shapepipe/modules/make_cat_package/make_cat.py +++ b/src/shapepipe/modules/make_cat_package/make_cat.py @@ -319,7 +319,10 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): observation raised, see ``ngmix.do_ngmix_metacal`` — ``NGMIX_MCAL_TYPES_FAIL``, ``NGMIX_NEIGHBOUR_FLAG`` — a blend flag set when the coadd segmentation stamp contains a non-central - footprint — and the defect diagnostics ``NGMIX_N_EPOCH_INTERP``, + footprint — the adaptive defect diagnostics + ``NGMIX_IMAFLAGS_ISO`` and ``NGMIX_N_EPOCH_INTERP_ISO`` (defined by + ``Postage_stamp.record_footprint_defects``; unknown footprints -1, + never-fit rows 0), and the stamp diagnostics ``NGMIX_N_EPOCH_INTERP``, ``NGMIX_MIN_DIST_INTERP``, ``NGMIX_MIN_DIST_NOISEFILL``, described by ``Postage_stamp.record_defects``; -1 is their absent distance, and a never-fit row carries 0 counts and -1 distances). @@ -378,6 +381,8 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): defect_columns = { "NGMIX_N_EPOCH_FAILED": "n_epoch_failed", "NGMIX_N_EPOCH_INTERP": "n_epoch_interp", + "NGMIX_IMAFLAGS_ISO": "imaflags_iso", + "NGMIX_N_EPOCH_INTERP_ISO": "n_epoch_interp_iso", "NGMIX_MIN_DIST_INTERP": "min_dist_interp", "NGMIX_MIN_DIST_NOISEFILL": "min_dist_noisefill", } @@ -409,6 +414,8 @@ def _save_ngmix_data(self, ngmix_cat_path, moments=False): self._add2dict("NGMIX_NEIGHBOUR_FLAG", np.zeros(n_obj)) self._add2dict("NGMIX_N_EPOCH_FAILED", np.zeros(n_obj)) self._add2dict("NGMIX_N_EPOCH_INTERP", np.zeros(n_obj)) + self._add2dict("NGMIX_IMAFLAGS_ISO", np.zeros(n_obj, dtype=np.int64)) + self._add2dict("NGMIX_N_EPOCH_INTERP_ISO", np.zeros(n_obj, dtype=np.int64)) self._add2dict("NGMIX_MIN_DIST_INTERP", np.full(n_obj, -1.0)) self._add2dict("NGMIX_MIN_DIST_NOISEFILL", np.full(n_obj, -1.0)) diff --git a/src/shapepipe/modules/ngmix_package/ngmix.py b/src/shapepipe/modules/ngmix_package/ngmix.py index 43995f926..04009baf5 100644 --- a/src/shapepipe/modules/ngmix_package/ngmix.py +++ b/src/shapepipe/modules/ngmix_package/ngmix.py @@ -272,6 +272,8 @@ def empty_metacal_output(): 'neighbour_flag', # defect diagnostics (Postage_stamp.record_defects) 'n_epoch_interp', + 'imaflags_iso', + 'n_epoch_interp_iso', 'min_dist_interp', 'min_dist_noisefill', 'nfev_fit', @@ -706,6 +708,8 @@ def __init__( # Defect diagnostics over the epochs kept (see # :meth:`record_defects`). self.n_epoch_interp = 0 + self.imaflags_iso = 0 + self.n_epoch_interp_iso = 0 self.min_dist_interp = np.inf self.min_dist_noisefill = np.inf self.bkg_sub = bkg_sub @@ -742,10 +746,61 @@ def record_defects(self, defect, interpolated): self.min_dist_noisefill, distance[noisefill].min() ) + def record_footprint_defects(self, flag, interpolated, seg, obj_id): + """Accumulate adaptive defect diagnostics for one kept epoch. + + @sc [decision:shape_measurement.defect_flagging,label:schema] footprint-defect-flags + ``imaflags_iso`` is the bitwise OR of all instrument bits in the + exposure flag stamp on ``seg == obj_id``, OR'ed across epochs kept + after the masked-fraction cut and central veto. ``seg`` is the coadd + SEG_VIGNET, already MegaCam-flipped to the epoch frame; it is not an + exposure detection or a reprojected isophote. ``n_epoch_interp_iso`` + counts kept epochs with any interpolated defect on that footprint. + Zero-weight/off-tile pixels with flag 0 contribute no instrument bit + (SExtractor semantics), but count if interpolated. Neither quantity + includes clean quarter-turn copies introduced by defect weighting. + Missing segmentation, an empty footprint, or a footprint touching + any stamp edge makes both quantities -1 (unknown/incomplete), not + clean; this sentinel persists across epochs. Never-fit rows in the + final catalogue carry 0, as do the other count diagnostics. + These diagnostics do not change the fit or impose a catalogue cut. + + Parameters + ---------- + flag : numpy.ndarray + Exposure instrument flags; integral floats are accepted because + vignet extraction can store flags in floating-point arrays. + Non-finite, fractional, negative or int64-overflow values on + the footprint raise ValueError rather than silently losing bits. + interpolated : numpy.ndarray of bool + Interpolated defects, as used by the fill and central veto. + seg : numpy.ndarray or None + Coadd segmentation stamp in the epoch's orientation. + obj_id : int + Tile object's segmentation label. + """ + footprint = None if seg is None else seg == obj_id + if (footprint is None or not footprint.any() + or footprint[[0, -1], :].any() + or footprint[:, [0, -1]].any()): + self.imaflags_iso = self.n_epoch_interp_iso = -1 + return + if self.imaflags_iso == -1: + return + values = np.asarray(flag)[footprint] + if (not np.isfinite(values).all() or (values < 0).any() + or (values >= 2**63).any() + or (values != np.floor(values)).any()): + raise ValueError("Non-integral or out-of-range footprint flags") + self.imaflags_iso |= int(np.bitwise_or.reduce(values.astype(np.int64))) + self.n_epoch_interp_iso += int((interpolated & footprint).any()) + def defect_diagnostics(self): """The defect diagnostic columns, -1 for an absent distance.""" return { "n_epoch_interp": self.n_epoch_interp, + "imaflags_iso": self.imaflags_iso, + "n_epoch_interp_iso": self.n_epoch_interp_iso, "min_dist_interp": ( float(self.min_dist_interp) if np.isfinite(self.min_dist_interp) else -1.0 @@ -1082,7 +1137,8 @@ def compile_results(self, results): # Defect diagnostics (Postage_stamp.record_defects), likewise # object-level and replicated. for key in ( - "n_epoch_interp", "min_dist_interp", "min_dist_noisefill" + "n_epoch_interp", "min_dist_interp", "min_dist_noisefill", + "imaflags_iso", "n_epoch_interp_iso" ): output_dict[name][key].append(results[idx][key]) # ngmix 2.x reports the solver's function-evaluation count @@ -1281,8 +1337,13 @@ def log_mean_ellipticity(self): """Log mean ellipticity from NOSHEAR HDU to the run log. Reports , with standard errors for all objects and for - objects passing the default metacal cuts (flags==0, mcal_flags==0, - 10 < SNR < 500, T/Tpsf > 0.5). + objects passing diagnostic cuts (flags==0, mcal_flags==0, + 10 <= SNR <= 500, T/T_psf_reconv >= 0.5). + + @sc [label:schema] mean-ellipticity-reconvolution-psf + Use the reconvolution PSF family, as sp_validation.calibration.metacal + does in _read_data_ngmix and _masking_gal for its relative-size cut. + This summary only reads the catalogue; it does not select saved rows. """ output_path = self.get_output_path(self._output_dir) try: @@ -1294,7 +1355,7 @@ def log_mean_ellipticity(self): mcal_flags = d['mcal_flags'] s2n = d['s2n'].astype(float) T = d['T'].astype(float) - Tpsf = d['Tpsf'].astype(float) + Tpsf = d['T_psf_reconv'].astype(float) except Exception as e: self._w_log.warning(f"Could not compute mean ellipticity: {e}") return @@ -1330,7 +1391,7 @@ def _log_stats(g1_sel, g2_sel, label): & (s2n <= 500.0) & (T / Tpsf >= 0.5) ) - _log_stats(g1[mask], g2[mask], "SNR in [10, 500], T/Tpsf > 0.5") + _log_stats(g1[mask], g2[mask], "SNR in [10, 500], T/T_psf_reconv >= 0.5") def process(self): """Process. @@ -1540,7 +1601,7 @@ def process(self): self.save_results(res_dict) saved_batch_cumul += count_batch self._w_log.info( - f"Batch-saved {count_batch} ({len(res_dict)} valid) objects," + f"Batch-saved {count_batch} objects," + f" cumul={saved_batch_cumul}" ) final_res = [] @@ -1733,6 +1794,7 @@ def prepare_postage_stamps( stamp.epoch_cuts["central_veto"] += 1 continue stamp.record_defects(defect, interpolated) + stamp.record_footprint_defects(flag_vign, interpolated, tile_seg, obj_id) # One unpickle per exposure (all CCDs), reused across this object's # epochs; the cache is per-call, so bounded by the object's exposures diff --git a/tests/module/test_make_cat.py b/tests/module/test_make_cat.py index 4a725f019..528b43ac0 100644 --- a/tests/module/test_make_cat.py +++ b/tests/module/test_make_cat.py @@ -60,6 +60,7 @@ def warning(self, message, *_args, **_kwargs): "mcal_types_fail", "neighbour_flag", "n_epoch_interp", "min_dist_interp", "min_dist_noisefill", + "imaflags_iso", "n_epoch_interp_iso", "nfev_fit", "g1", "g1_err", "g2", "g2_err", "T", "T_err", @@ -90,6 +91,7 @@ def _ngmix_row(obj_id): "n_epoch_failed": 1, "mcal_types_fail": 0, "neighbour_flag": 1, + "imaflags_iso": 1028, "n_epoch_interp_iso": 1, "n_epoch_interp": 2, "min_dist_interp": 7.5, "min_dist_noisefill": 14.25, "nfev_fit": 7, @@ -124,6 +126,7 @@ def _write_ngmix_cat(path, obj_ids): name=key, format="K" if key in ("id", "n_epoch_model", "n_epoch_failed", "mcal_types_fail", + "imaflags_iso", "n_epoch_interp_iso", "n_epoch_interp", "nfev_fit", "flags", "mcal_flags") else "D", array=np.array([row[key] for row in rows]), @@ -286,6 +289,8 @@ def test_save_ngmix_data_fills_sentinels_for_absent_objects(tmp_path): for col, key, never_fit in ( ("NGMIX_N_EPOCH_FAILED", "n_epoch_failed", 0.0), ("NGMIX_N_EPOCH_INTERP", "n_epoch_interp", 0.0), + ("NGMIX_IMAFLAGS_ISO", "imaflags_iso", 0), + ("NGMIX_N_EPOCH_INTERP_ISO", "n_epoch_interp_iso", 0), ("NGMIX_MIN_DIST_INTERP", "min_dist_interp", -1.0), ("NGMIX_MIN_DIST_NOISEFILL", "min_dist_noisefill", -1.0), ): @@ -334,6 +339,8 @@ def _metacal_result(obj_id): "n_epoch_failed": row["n_epoch_failed"], "neighbour_flag": row["neighbour_flag"], "n_epoch_interp": row["n_epoch_interp"], + "imaflags_iso": row["imaflags_iso"], + "n_epoch_interp_iso": row["n_epoch_interp_iso"], "min_dist_interp": row["min_dist_interp"], "min_dist_noisefill": row["min_dist_noisefill"], } diff --git a/tests/module/test_ngmix.py b/tests/module/test_ngmix.py index 37b3e58b5..94cf5a571 100644 --- a/tests/module/test_ngmix.py +++ b/tests/module/test_ngmix.py @@ -332,6 +332,7 @@ def _fake_metacal_result(T, T_err, T_psf, T_psf_err): "mcal_types_fail": 0, "neighbour_flag": 0, "n_epoch_interp": 0, + "imaflags_iso": 0, "n_epoch_interp_iso": 0, "min_dist_interp": -1.0, "min_dist_noisefill": -1.0, # original image PSF (psfex/mccd) family @@ -423,6 +424,45 @@ def test_compile_results_psf_families_are_unaliased(): assert noshear["g2_psf_orig"] != noshear["g2_psf_reconv"] +def test_log_mean_ellipticity_reads_compiled_reconvolution_psf(tmp_path): + """The run-log cut uses the reconvolution PSF, without changing the FITS. + + A stale Tpsf column prevents any summary; using the original PSF instead + would incorrectly admit the second object to this diagnostic selection. + """ + from pathlib import Path + from types import SimpleNamespace + + from shapepipe.modules.ngmix_package.ngmix import Ngmix + + inst = object.__new__(Ngmix) + inst._zero_point = 30.0 + inst._output_dir = str(tmp_path) + inst._file_number_string = "-000-000" + messages, warnings = [], [] + inst._w_log = SimpleNamespace(info=messages.append, warning=warnings.append) + kept = _fake_metacal_result(0.18, 0.02, 0.09, 0.001) + rejected = _fake_metacal_result(0.05, 0.01, 0.2, 0.001) + rejected["obj_id"] = 2 + rejected["noshear"]["g"] = [-0.3, 0.4] + compiled = inst.compile_results([kept, rejected]) + assert "Tpsf" not in compiled["noshear"] + inst.save_results(compiled) + path = Path(inst.get_output_path(str(tmp_path))) + before = path.read_bytes() + + inst.log_mean_ellipticity() + + assert not warnings + assert any("no cuts (N=2)" in line for line in messages) + selected = [line for line in messages if "T/T_psf_reconv" in line] + assert len(selected) == 1 + assert "(N=1)" in selected[0] + assert " = +1.0000e-02" in selected[0] + assert " = -2.0000e-02" in selected[0] + assert path.read_bytes() == before + + def test_compile_results_nan_fills_failed_fit_types(): """A failed fit type must be recorded with NaNs, not crash the tile. @@ -638,6 +678,7 @@ def test_process_counts_flagged_fits_across_batches(tmp_path, monkeypatch, flags epoch_cuts=Counter(considered=1), epoch_failures=[], jacobs=[galsim.JacobianWCS(.186, 0., 0., .186)], defect_diagnostics=lambda: dict( + imaflags_iso=0, n_epoch_interp_iso=0, n_epoch_interp=0, min_dist_interp=-1.0, min_dist_noisefill=-1.0, ), ) @@ -1334,6 +1375,7 @@ def test_process_logs_and_counts_failed_epochs(tmp_path, monkeypatch): epoch_failures=[], jacobs=[galsim.JacobianWCS(.186, 0., 0., .186)] * 2, defect_diagnostics=lambda: dict( + imaflags_iso=0, n_epoch_interp_iso=0, n_epoch_interp=0, min_dist_interp=-1.0, min_dist_noisefill=-1.0, ), ) diff --git a/tests/module/test_ngmix_defect_fill.py b/tests/module/test_ngmix_defect_fill.py index 96de7e44c..b24bf971f 100644 --- a/tests/module/test_ngmix_defect_fill.py +++ b/tests/module/test_ngmix_defect_fill.py @@ -764,6 +764,7 @@ def test_defect_diagnostics_summarise_the_kept_epochs(): ) assert stamp.defect_diagnostics() == { "n_epoch_interp": 2, + "imaflags_iso": -1, "n_epoch_interp_iso": -1, "min_dist_interp": float(EPOCH_INTERPOLATED_DEFECT_RADIUS), "min_dist_noisefill": float(EPOCH_CENTRAL_DEFECT_RADIUS), } @@ -774,6 +775,7 @@ def test_defect_diagnostics_summarise_the_kept_epochs(): psf_obj=psf_obj, gal_obj=gal_obj, ) assert stamp.defect_diagnostics() == { + "imaflags_iso": -1, "n_epoch_interp_iso": -1, "n_epoch_interp": 0, "min_dist_interp": -1.0, "min_dist_noisefill": -1.0, } diff --git a/tests/module/test_ngmix_defect_flags.py b/tests/module/test_ngmix_defect_flags.py new file mode 100644 index 000000000..53118cf5a --- /dev/null +++ b/tests/module/test_ngmix_defect_flags.py @@ -0,0 +1,109 @@ +"""Adaptive defect flags on the coadd object's footprint, over kept epochs.""" + +import numpy as np +import pytest + +from shapepipe.modules.ngmix_package.ngmix import Ngmix, prepare_postage_stamps +from tests.module.test_ngmix_defect_fill import N_STAMP, _fake_inputs + + +def _run(kind="inside", seg_kind="valid", megacam_flip=True): + seg = np.zeros((N_STAMP, N_STAMP), dtype=np.int32) + seg[23:28, 23:37] = 1 # reaches 11 px right of centre, no stamp edge + if seg_kind == "empty": + seg[:] = 0 + if seg_kind == "clipped": + seg[25, 36:] = 1 + epochs = {} + ccds = (10, 20, 21) if megacam_flip else (20, 21, 22) + for i, (ccd, bit) in enumerate(zip(ccds, (4, 1024, 0))): + flag = np.zeros(seg.shape, dtype=np.float32) + weight = np.ones(seg.shape) + pixel = (25, 35) if kind != "outside" else (35, 35) + if kind != "clean": + flag[pixel] = bit + if kind == "zero_weight" and bit: + flag[pixel] = 0 + weight[pixel] = 0 + if megacam_flip: + flag = Ngmix.MegaCamFlip(flag, ccd) + weight = Ngmix.MegaCamFlip(weight, ccd) + epochs[f"210000{i}-{ccd}"] = (flag, weight) + # A centrally vetoed epoch and a fraction-cut epoch carry otherwise + # unique bits on the footprint. Neither may contribute. + for i, fraction in enumerate((False, True), start=3): + flag = np.zeros(seg.shape, dtype=np.float32) + flag[25, 35] = 64 + if fraction: + flag[:10] = 64 + else: + flag[25, 25] = 64 + epochs[f"210000{i}-20"] = (flag, np.ones(seg.shape)) + vignet, tile, psf, gal = _fake_inputs(epochs) + tile.seg = None if seg_kind == "missing" else seg[None] + return prepare_postage_stamps( + vignet, 1, 0, tile, bkg_sub=False, psf_obj=psf, gal_obj=gal, + ) + + +@pytest.mark.decision("shape_measurement.defect_flagging") +@pytest.mark.parametrize("flip", [True, False]) +@pytest.mark.parametrize("kind, expected", [ + ("inside", (1028, 2)), ("outside", (0, 0)), + ("clean", (0, 0)), ("zero_weight", (0, 2)), +]) +def test_footprint_flags_only_count_kept_epochs(kind, expected, flip): + """Keep bit identities across epochs, ignore outside/vetoed defects. + + The final clean epoch must not overwrite the OR of earlier epochs. + A zero-weight pixel has no instrument bit but counts if interpolated. + """ + stamp = _run(kind, megacam_flip=flip) + assert len(stamp.flags) == 3 + diagnostics = stamp.defect_diagnostics() + assert (diagnostics["imaflags_iso"], + diagnostics["n_epoch_interp_iso"]) == expected + assert diagnostics["n_epoch_interp"] == (0 if kind == "clean" else 2) + + +@pytest.mark.parametrize("seg_kind", ["missing", "empty", "clipped"]) +def test_unknown_or_incomplete_footprint_is_not_clean(seg_kind): + diagnostics = _run(seg_kind=seg_kind).defect_diagnostics() + assert diagnostics["imaflags_iso"] == -1 + assert diagnostics["n_epoch_interp_iso"] == -1 + + +@pytest.mark.parametrize("value", [np.nan, np.inf, -1, 1.5, float(2**63)]) +def test_non_integral_or_unsafe_flags_raise(value): + """Do not silently truncate float flags or overflow the integer OR.""" + from shapepipe.modules.ngmix_package.ngmix import Postage_stamp + + stamp = Postage_stamp() + seg = np.zeros((5, 5), dtype=int) + seg[2, 2] = 7 + flag = np.zeros(seg.shape) + flag[2, 2] = value + with pytest.raises(ValueError, match="footprint flags"): + stamp.record_footprint_defects(flag, flag != 0, seg, 7) + + +def test_footprint_columns_survive_ngmix_fits(tmp_path): + """Every shear extension keeps integer diagnostics, including -1.""" + from astropy.io import fits + from tests.module.test_ngmix import _fake_metacal_result + + inst = object.__new__(Ngmix) + inst._zero_point = 30.0 + inst._output_dir = str(tmp_path) + inst._file_number_string = "-001-001" + rows = [] + for i, (bits, count) in enumerate(((1028, 2), (-1, -1)), start=1): + row = _fake_metacal_result(.18, .02, .09, .001) + row.update(obj_id=i, imaflags_iso=bits, n_epoch_interp_iso=count) + rows.append(row) + inst.save_results(inst.compile_results(rows)) + with fits.open(inst.get_output_path(str(tmp_path))) as hdus: + for hdu in hdus[1:]: + assert hdu.data["imaflags_iso"].tolist() == [1028, -1] + assert hdu.data["n_epoch_interp_iso"].tolist() == [2, -1] + assert hdu.data["imaflags_iso"].dtype.kind == "i" diff --git a/tests/module/test_psf_averaging_properties.py b/tests/module/test_psf_averaging_properties.py index d6799bc29..49895b3e4 100644 --- a/tests/module/test_psf_averaging_properties.py +++ b/tests/module/test_psf_averaging_properties.py @@ -260,6 +260,8 @@ def test_all_epochs_failed_raises_zero_division(flagged_specs): "NGMIX_NEIGHBOUR_FLAG": 0.0, "NGMIX_N_EPOCH_FAILED": 0.0, "NGMIX_N_EPOCH_INTERP": 0.0, + "NGMIX_IMAFLAGS_ISO": 0, + "NGMIX_N_EPOCH_INTERP_ISO": 0, "NGMIX_MIN_DIST_INTERP": -1.0, "NGMIX_MIN_DIST_NOISEFILL": -1.0, } @@ -270,6 +272,7 @@ def test_all_epochs_failed_raises_zero_division(flagged_specs): "id", "n_epoch_model", "n_epoch_failed", "mcal_types_fail", "neighbour_flag", "n_epoch_interp", "min_dist_interp", "min_dist_noisefill", "nfev_fit", + "imaflags_iso", "n_epoch_interp_iso", "g1", "g1_err", "g2", "g2_err", "T", "T_err", "flux", "flux_err", "s2n", "mag", "mag_err", "flags", "mcal_flags", "g1_psf_orig", "g2_psf_orig", "g1_err_psf_orig", "g2_err_psf_orig", @@ -280,7 +283,8 @@ def test_all_epochs_failed_raises_zero_division(flagged_specs): _INT_KEYS = { "id", "n_epoch_model", "n_epoch_failed", "mcal_types_fail", "neighbour_flag", - "n_epoch_interp", "nfev_fit", "flags", "mcal_flags" + "n_epoch_interp", "nfev_fit", "flags", "mcal_flags", + "imaflags_iso", "n_epoch_interp_iso" } _SHEAR_EXTS = ["1M", "1P", "2M", "2P", "NOSHEAR"] diff --git a/tests/module/test_psf_grammar_properties.py b/tests/module/test_psf_grammar_properties.py index 4a6442f13..d35f6a462 100644 --- a/tests/module/test_psf_grammar_properties.py +++ b/tests/module/test_psf_grammar_properties.py @@ -67,11 +67,13 @@ def info(self, *_args, **_kwargs): "id", "n_epoch_model", "n_epoch_failed", "mcal_types_fail", "neighbour_flag", "n_epoch_interp", "nfev_fit", "flags", "mcal_flags", + "imaflags_iso", "n_epoch_interp_iso", } NGMIX_KEYS = [ "id", "n_epoch_model", "n_epoch_failed", "mcal_types_fail", "neighbour_flag", "n_epoch_interp", "min_dist_interp", "min_dist_noisefill", "nfev_fit", + "imaflags_iso", "n_epoch_interp_iso", "g1", "g1_err", "g2", "g2_err", "T", "T_err", "flux", "flux_err", "s2n", "mag", "mag_err", @@ -113,6 +115,7 @@ def _base_row(obj_id): return { "id": obj_id, "n_epoch_model": 3, "n_epoch_failed": 1, "mcal_types_fail": 0, "neighbour_flag": 1, + "imaflags_iso": 1028, "n_epoch_interp_iso": 1, "n_epoch_interp": 2, "min_dist_interp": 7.5, "min_dist_noisefill": 14.25, "nfev_fit": 7, "g1": 0.10, "g1_err": 0.011, "g2": -0.20, "g2_err": 0.022, @@ -184,7 +187,8 @@ def _run_save_ngmix(ngmix_path, obj_ids, cat_size_target=None): GRAMMAR_RE = re.compile( rf"^NGMIXm?_(?:{_COMPONENT})(?:_ERR)?(?:_(?:{_OBJECT}))?_(?:{_SHEAR})$" rf"|^NGMIXm?_(?:MCAL_FLAGS|MCAL_TYPES_FAIL|N_EPOCH|NEIGHBOUR_FLAG" - rf"|N_EPOCH_FAILED|N_EPOCH_INTERP|MIN_DIST_INTERP|MIN_DIST_NOISEFILL)$" + rf"|N_EPOCH_FAILED|N_EPOCH_INTERP|MIN_DIST_INTERP|MIN_DIST_NOISEFILL" + rf"|IMAFLAGS_ISO|N_EPOCH_INTERP_ISO)$" ) # The shipped final-catalogue param files, two levels up from tests/module/. @@ -424,7 +428,8 @@ def test_param_file_ngmix_tokens_are_producible(param_path, obj_ids): r"NGMIXm?_(?:G1|G2|T|SNR|FLUX|MAG|FLAGS)(?:_ERR)?" r"(?:_(?:PSF_ORIG|PSF_RECONV))?_(?:NOSHEAR|1P|1M|2P|2M)" r"|NGMIXm?_(?:MCAL_FLAGS|MCAL_TYPES_FAIL|N_EPOCH|NEIGHBOUR_FLAG" - r"|N_EPOCH_FAILED|N_EPOCH_INTERP|MIN_DIST_INTERP|MIN_DIST_NOISEFILL)" + r"|N_EPOCH_FAILED|N_EPOCH_INTERP|MIN_DIST_INTERP|MIN_DIST_NOISEFILL" + r"|IMAFLAGS_ISO|N_EPOCH_INTERP_ISO)" # HSM: g-type, explicit PSF/STAR object, singular FLAG; the multi-epoch # sink in make_cat._save_psf_data appends a bare epoch index. r"|HSM_(?:G1|G2|T|M4_1|M4_2|RHO4)_(?:PSF|STAR)(?:_\d+)?" diff --git a/tests/workflow/test_defect_flags.py b/tests/workflow/test_defect_flags.py new file mode 100644 index 000000000..72334e49a --- /dev/null +++ b/tests/workflow/test_defect_flags.py @@ -0,0 +1,50 @@ +"""The committed detection config supplies the adaptive-flag footprint.""" + +from pathlib import Path +import re + +import numpy as np +import pytest + +from tests.workflow.harness import Campaign +from tests.workflow.test_dag import _committed_value + + +@pytest.mark.parametrize("input_type", ["data", "image_sims"]) +@pytest.mark.parametrize("blend", [None, "noisefill", "uberseg"]) +def test_committed_detection_supplies_defect_flag_footprint( + tmp_path, resolve_dag, monkeypatch, input_type, blend): + """Every handling enables SEG_VIGNET and Tile_cat consumes that column. + + This guards the shipped path, rather than assigning tile.seg by hand: + removing the workflow export, disabling the committed ini option, or + dropping the loader's SEG_VIGNET read must not silently turn an ordinary + clean object into an unknown-footprint (-1) object. The SExtractor runner + column-writing behavior is covered by test_sextractor_seg_vignet. + """ + from shapepipe.modules.ngmix_package.ngmix import Postage_stamp, Tile_cat + from tests.module.test_ngmix_uberseg import _write_tile_cat + + campaign = Campaign(tmp_path / "campaign", input_type, "psfex") + if blend is not None: + campaign.config["blend_handling"] = blend + campaign.write_config() + with resolve_dag(campaign) as dag: + job = next(job for job in dag.jobs if job.rule.name == "tile_detect") + exports = dict(re.findall(r"^export (SP_SEG_VIGNET)='([^']*)'$", + job.params.pre, re.MULTILINE)) + config_dir = (Path(__file__).parents[2] / "workflow" / "config" + / ("cfis" if input_type == "data" else "cfis_image_sims")) + enabled = _committed_value(job.shellcmd, config_dir, "SEG_VIGNET", + exports, monkeypatch) + assert enabled == "True", "committed detection disables SEG_VIGNET" + seg = np.zeros((3, 5, 5), dtype=np.int32) + seg[:, 2, 2] = [4, 5, 9] + path = _write_tile_cat(tmp_path / "tile.fits", seg) + tile = Tile_cat(str(path)) + stamp = Postage_stamp() + stamp.record_footprint_defects( + np.zeros((5, 5), dtype=np.int16), np.zeros((5, 5), dtype=bool), + None if tile.seg is None else tile.seg[0], tile.obj_id[0], + ) + assert stamp.imaflags_iso == stamp.n_epoch_interp_iso == 0 diff --git a/universes/committed.yaml b/universes/committed.yaml index 79f62f356..13b5e3503 100644 --- a/universes/committed.yaml +++ b/universes/committed.yaml @@ -56,6 +56,7 @@ analyses: psf_likelihood_noise: psf_noise_1em5 megacam_ccd_flip: megapipe_flip defect_fill: interpolate + defect_flagging: footprint_or_and_interp_epochs blend_handling: uberseg defect_weighting: fourfold_zero epoch_masked_fraction_cut: ten_percent diff --git a/workflow/config/cfis/final_cat.param b/workflow/config/cfis/final_cat.param index 1ca2774db..26da44a2e 100644 --- a/workflow/config/cfis/final_cat.param +++ b/workflow/config/cfis/final_cat.param @@ -178,6 +178,9 @@ NGMIX_NEIGHBOUR_FLAG # min distance (px) from the centre to an interpolated / noise-filled defect # pixel, -1 when none (ngmix Postage_stamp.record_defects) NGMIX_N_EPOCH_INTERP +# @sc [decision:shape_measurement.defect_flagging] +NGMIX_IMAFLAGS_ISO +NGMIX_N_EPOCH_INTERP_ISO NGMIX_MIN_DIST_INTERP NGMIX_MIN_DIST_NOISEFILL diff --git a/workflow/config/cfis_image_sims/final_cat.param b/workflow/config/cfis_image_sims/final_cat.param index 99fe70088..119477bb6 100644 --- a/workflow/config/cfis_image_sims/final_cat.param +++ b/workflow/config/cfis_image_sims/final_cat.param @@ -31,6 +31,9 @@ NGMIX_NEIGHBOUR_FLAG # min distance (px) from the centre to an interpolated / noise-filled defect # pixel, -1 when none (ngmix Postage_stamp.record_defects) NGMIX_N_EPOCH_INTERP +# @sc [decision:shape_measurement.defect_flagging] +NGMIX_IMAFLAGS_ISO +NGMIX_N_EPOCH_INTERP_ISO NGMIX_MIN_DIST_INTERP NGMIX_MIN_DIST_NOISEFILL