From 49fb0b1680e8eeec5317ebf60656bc15c2e7e181 Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 11:27:40 +0100 Subject: [PATCH 1/4] fix: pass HDUList to WCS in footprint filter (HST CPDIS lookup-table distortion) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Phase-3 regression: covers_target broke the ACS path — FLC headers carry lookup-table distortion astropy can only resolve with fobj; JWST cal files don't, so the four-band integration never tripped it. Caught by the refactor witness baseline. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WNhuuvZLgZjHvWBGdhUst3 --- autoreduce/acquire/footprint.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/autoreduce/acquire/footprint.py b/autoreduce/acquire/footprint.py index 32cbcad..c6ac6ee 100644 --- a/autoreduce/acquire/footprint.py +++ b/autoreduce/acquire/footprint.py @@ -32,7 +32,9 @@ def covers_target(path: Path, ra: float, dec: float, margin_arcsec: float) -> bo for hdu in hdul: if hdu.name != "SCI" or hdu.data is None: continue - wcs = WCS(hdu.header, naxis=2) + # fobj: HST headers carry lookup-table distortion (CPDIS/D2IM) + # that astropy can only resolve with the open HDUList in hand. + wcs = WCS(hdu.header, fobj=hdul, naxis=2) ny, nx = hdu.data.shape[-2], hdu.data.shape[-1] x, y = wcs.world_to_pixel_values(ra, dec) if not (np.isfinite(x) and np.isfinite(y)): From 7b34bf4de8bc07c76ffd163e7a51bbe1f51674cf Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 11:35:57 +0100 Subject: [PATCH 2/4] =?UTF-8?q?fix:=20widen=20ePSF=20extraction=20window?= =?UTF-8?q?=20(+20)=20=E2=80=94=20photutils=202.3=20recentering=20crash?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Drifting stars at maxiters=10 hit a malformed overlap_slices inside photutils 2.3 when the fitting region overruns the +10 window; +20 restores its designed warn-and-skip handling. Broke the ACS path on main (star set changed by the phase-2/3 selection fixes); JWST ePSFs shift marginally with the wider window — disclosed on issue #8. edge_margin follows the window (46). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WNhuuvZLgZjHvWBGdhUst3 --- autoreduce/psf/epsf.py | 7 +++++-- autoreduce/psf/stars.py | 4 ++-- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/autoreduce/psf/epsf.py b/autoreduce/psf/epsf.py index 9e066df..e23744e 100644 --- a/autoreduce/psf/epsf.py +++ b/autoreduce/psf/epsf.py @@ -58,8 +58,11 @@ def build_epsf( positions = Table( {"x": stars_table["xcentroid"], "y": stars_table["ycentroid"]} ) - # Extraction window comfortably larger than the extended kernel. - size = max(psf_full_shape) + 10 + # Extraction window comfortably larger than the extended kernel — the + # +20 pad gives EPSFBuilder's recentering iterations room; with only +10, + # drifting stars trip a photutils 2.3 internal crash (malformed + # overlap_slices) instead of its designed warn-and-skip handling. + size = max(psf_full_shape) + 20 if size % 2 == 0: size += 1 stars = extract_stars(NDData(sci), positions, size=size) diff --git a/autoreduce/psf/stars.py b/autoreduce/psf/stars.py index 11d0f24..7680e34 100644 --- a/autoreduce/psf/stars.py +++ b/autoreduce/psf/stars.py @@ -23,9 +23,9 @@ class StarSelection: round_limit: float = 0.3 saturation_fraction: float = 0.7 # of adapter.saturation_dn, in counts min_separation_pix: float = 25.0 - # Must exceed half the ePSF extraction window (psf_full 61 + 10 pad -> 36) + # Must exceed half the ePSF extraction window (psf_full 61 + 20 pad -> 41) # or edge stars pass selection only to be dropped at extraction. - edge_margin_pix: int = 36 + edge_margin_pix: int = 46 exclusion_radius_pix: float = 50.0 # around the target itself From 869dfb52a0742123df373dda3c434048c764fcbb Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 11:54:11 +0100 Subject: [PATCH 3/4] refactor: post-phase-3 consolidation (behaviour-preserving) - reduce_target decomposed into per-stage functions with a _StageContext; record layout and key order unchanged. - autoreduce/validation/parity.py extracts sub-pixel registration + ratio statistics; slacs0008 and cosmos_web_ring scripts thinned onto it (registered_ratios is a superset of both originals). - drizzle/_common.py: chdir_scratch context manager + combine_provenance (canonical key order preserved); both backends rewritten onto them. - reject_crowded vectorized (broadcasted distance matrix); randomized equivalence test against the original loop. 85 tests. Witness byte-comparison vs _baseline_refactor follows. Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WNhuuvZLgZjHvWBGdhUst3 --- autoreduce/drizzle/_common.py | 61 ++++++++++ autoreduce/drizzle/combine.py | 40 +++---- autoreduce/drizzle/jwst_combine.py | 52 ++++----- autoreduce/pipeline.py | 145 ++++++++++++++++-------- autoreduce/psf/stars.py | 11 +- autoreduce/validation/__init__.py | 7 ++ autoreduce/validation/parity.py | 80 +++++++++++++ scripts/reduce_cosmos_web_ring.py | 53 +-------- scripts/reduce_slacs0008.py | 46 +------- test_autoreduce/test_psf_and_package.py | 54 +++++++++ 10 files changed, 344 insertions(+), 205 deletions(-) create mode 100644 autoreduce/drizzle/_common.py create mode 100644 autoreduce/validation/__init__.py create mode 100644 autoreduce/validation/parity.py diff --git a/autoreduce/drizzle/_common.py b/autoreduce/drizzle/_common.py new file mode 100644 index 0000000..b1001f9 --- /dev/null +++ b/autoreduce/drizzle/_common.py @@ -0,0 +1,61 @@ +""" +Shared combine-backend plumbing: the scratch-directory discipline and the +provenance fragment both backends emit identically. +""" + +import os +from contextlib import contextmanager +from pathlib import Path +from typing import Dict, List + +import numpy as np + +from ..instruments import InstrumentAdapter +from ..noise.rms import casertano_r +from ..target import TargetSpec +from .diagnostics import check_weight_uniformity + + +@contextmanager +def chdir_scratch(output_dir: Path): + """Resolve, create and chdir into a backend scratch dir; always restore.""" + output_dir = Path(output_dir).resolve() + output_dir.mkdir(parents=True, exist_ok=True) + cwd = os.getcwd() + os.chdir(output_dir) + try: + yield output_dir + finally: + os.chdir(cwd) + + +def combine_provenance( + spec: TargetSpec, + adapter: InstrumentAdapter, + exposures: List[Path], + wht: np.ndarray, + kwargs_key: str, + kwargs: Dict, + head: Dict = None, + tail: Dict = None, +) -> Dict: + """ + The provenance dict every combine backend records, assembled in the + canonical key order (kept stable — reduction.json is byte-compared by the + refactor witnesses). `head`/`tail` carry backend-specific extras. + """ + out = dict(head or {}) + out.update( + { + "n_exposures": len(exposures), + "exposures": [Path(p).name for p in exposures], + "single_exposure_branch": len(exposures) == 1, + kwargs_key: kwargs, + "correlated_noise_factor": casertano_r( + spec.final_pixfrac, adapter.scale_ratio(spec.final_scale) + ), + "weight_uniformity": check_weight_uniformity(wht), + } + ) + out.update(tail or {}) + return out diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py index fef6b42..8494fc5 100644 --- a/autoreduce/drizzle/combine.py +++ b/autoreduce/drizzle/combine.py @@ -15,8 +15,6 @@ from ..instruments import InstrumentAdapter from ..target import TargetSpec -from .diagnostics import check_weight_uniformity -from ..noise.rms import casertano_r def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int) -> Dict: @@ -66,28 +64,21 @@ def combine( from astropy.io import fits from drizzlepac import astrodrizzle - output_dir = Path(output_dir) - output_dir.mkdir(parents=True, exist_ok=True) + from ._common import chdir_scratch, combine_provenance + # Drizzlepac lowercases output filenames internally, which breaks absolute - # paths containing capitals on case-sensitive filesystems — so chdir into - # the work dir and pass a relative, already-lowercase output root. This + # paths containing capitals on case-sensitive filesystems — so run inside + # the scratch dir with a relative, already-lowercase output root. This # also keeps AstroDrizzle's cwd scratch files contained. output_name = f"{spec.name}_{spec.filter_name}".lower() - output_root = str(output_dir / output_name) - kwargs = drizzle_kwargs_for(spec, adapter, len(exposures)) - import os - - cwd = os.getcwd() - os.chdir(output_dir) - try: + with chdir_scratch(output_dir) as output_dir: astrodrizzle.AstroDrizzle( input=[str(p) for p in exposures], output=output_name, **kwargs, ) - finally: - os.chdir(cwd) + output_root = str(output_dir / output_name) def _one(suffix: str) -> Path: hits = sorted(glob.glob(f"{output_root}*{suffix}")) @@ -100,15 +91,12 @@ def _one(suffix: str) -> Path: sci = _one("_sci.fits") wht = _one("_wht.fits") - wht_data = fits.getdata(wht) - provenance = { - "n_exposures": len(exposures), - "exposures": [Path(p).name for p in exposures], - "single_exposure_branch": len(exposures) == 1, - "drizzle_kwargs": {k: kwargs[k] for k in sorted(kwargs)}, - "correlated_noise_factor": casertano_r( - spec.final_pixfrac, adapter.scale_ratio(spec.final_scale) - ), - "weight_uniformity": check_weight_uniformity(wht_data), - } + provenance = combine_provenance( + spec, + adapter, + exposures, + fits.getdata(wht), + kwargs_key="drizzle_kwargs", + kwargs={k: kwargs[k] for k in sorted(kwargs)}, + ) return sci, wht, provenance diff --git a/autoreduce/drizzle/jwst_combine.py b/autoreduce/drizzle/jwst_combine.py index 589d047..73cc4ac 100644 --- a/autoreduce/drizzle/jwst_combine.py +++ b/autoreduce/drizzle/jwst_combine.py @@ -18,8 +18,6 @@ from ..instruments import InstrumentAdapter from ..target import TargetSpec -from .diagnostics import check_weight_uniformity -from ..noise.rms import casertano_r def combine( @@ -34,23 +32,18 @@ def combine( from jwst.associations.lib.rules_level3_base import DMS_Level3_Base from jwst.pipeline import Image3Pipeline - # Resolve before the chdir below: relative paths would dangle afterwards. - output_dir = Path(output_dir).resolve() - output_dir.mkdir(parents=True, exist_ok=True) - product_name = f"{spec.name}_{spec.filter_name}".lower() - - asn = asn_from_list( - [str(p) for p in exposures], rule=DMS_Level3_Base, product_name=product_name - ) - asn_path = output_dir / f"{product_name}_asn.json" - _, serialized = asn.dump(format="json") - asn_path.write_text(serialized) + from ._common import chdir_scratch, combine_provenance - import os + product_name = f"{spec.name}_{spec.filter_name}".lower() + with chdir_scratch(output_dir) as output_dir: + asn = asn_from_list( + [str(p) for p in exposures], rule=DMS_Level3_Base, + product_name=product_name, + ) + asn_path = output_dir / f"{product_name}_asn.json" + _, serialized = asn.dump(format="json") + asn_path.write_text(serialized) - cwd = os.getcwd() - os.chdir(output_dir) - try: Image3Pipeline.call( str(asn_path), output_dir=str(output_dir), @@ -68,8 +61,6 @@ def combine( }, }, ) - finally: - os.chdir(cwd) i2d = output_dir / f"{product_name}_i2d.fits" if not i2d.exists(): @@ -92,22 +83,19 @@ def combine( fits.PrimaryHDU(data, header=header).writeto(path, overwrite=True) paths[name] = path - provenance = { - "backend": "jwst_image3", - "n_exposures": len(exposures), - "exposures": [Path(p).name for p in exposures], - "single_exposure_branch": len(exposures) == 1, - "resample_kwargs": { + provenance = combine_provenance( + spec, + adapter, + exposures, + wht, + kwargs_key="resample_kwargs", + kwargs={ "pixel_scale": spec.final_scale, "pixfrac": spec.final_pixfrac, "kernel": spec.final_kernel, "rotation": 0.0, }, - "correlated_noise_factor": casertano_r( - spec.final_pixfrac, adapter.scale_ratio(spec.final_scale) - ), - "weight_uniformity": check_weight_uniformity(wht), - "err_path": str(paths["err"]), - "i2d_path": str(i2d), - } + head={"backend": "jwst_image3"}, + tail={"err_path": str(paths["err"]), "i2d_path": str(i2d)}, + ) return paths["sci"], paths["wht"], provenance diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index 444d6e1..0087ca8 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -8,8 +8,9 @@ photutils) are imported inside stages so the package imports without them. """ +from dataclasses import dataclass, field from pathlib import Path -from typing import Dict, Optional +from typing import Dict, List, Optional import numpy as np @@ -20,6 +21,8 @@ from .acquire import mast as mast_mod from .align import diagnostics as align_mod from .drizzle import combine as combine_mod +from .drizzle.diagnostics import check_weight_uniformity +from .instruments import InstrumentAdapter from .noise import rms as rms_mod from .package import cutout as cutout_mod from .package import provenance as provenance_mod @@ -28,24 +31,22 @@ from .target import TargetSpec -def reduce_target( - spec: TargetSpec, - cache_root: Path, - output_root: Path, - size_cap_bytes: Optional[int] = None, - evict_when_done: bool = False, -) -> Dict: - """Run the full pipeline for one target; returns the provenance record.""" - adapter = instruments.get(spec.instrument) - cache = cache_mod.ExposureCache(Path(cache_root), size_cap_bytes=size_cap_bytes) - out_dir = Path(output_root) / spec.name - out_dir.mkdir(parents=True, exist_ok=True) - work_dir = out_dir / "work" - work_dir.mkdir(exist_ok=True) +@dataclass +class _StageContext: + """Everything the stages share; `record` is the growing provenance.""" + + spec: TargetSpec + adapter: InstrumentAdapter + cache: cache_mod.ExposureCache + out_dir: Path + work_dir: Path + record: Dict = field(default_factory=dict) + exposures: List[Path] = field(default_factory=list) - record: Dict = {"target": spec.as_dict(), "instrument": adapter.key} - # -- acquire ------------------------------------------------------------ +def _acquire(ctx: _StageContext) -> None: + """Download (or reuse) exposures, sync references, footprint-filter.""" + spec, adapter, cache = ctx.spec, ctx.adapter, ctx.cache crds_mod.configure_environment(cache.references_dir, adapter) exposures = cache.exposures_for(spec.name) downloaded = False @@ -81,7 +82,8 @@ def reduce_target( exposures, skipped = footprint_mod.filter_to_target( exposures, spec.ra, spec.dec, margin_arcsec=cutout_extent + 15.0 ) - record["acquire"] = { + ctx.exposures = exposures + ctx.record["acquire"] = { "n_exposures": len(exposures), "exposures": [Path(p).name for p in exposures], "n_skipped_off_target": len(skipped), @@ -89,20 +91,23 @@ def reduce_target( "references_synced": refs_synced, } - # -- align ---------------------------------------------------------------- - record["align"] = { - "wcs_solutions": align_mod.wcs_solution_names(exposures), + +def _align(ctx: _StageContext) -> None: + ctx.record["align"] = { + "wcs_solutions": align_mod.wcs_solution_names(ctx.exposures), "tweakreg_run": False, # a-priori WCS accepted by default (stage 2) } - # -- drizzle --------------------------------------------------------------- - sci_path, wht_path, drizzle_prov = combine_mod.combine( - exposures, spec, adapter, work_dir - ) - record["drizzle"] = drizzle_prov +def _combine(ctx: _StageContext): + """Run the backend combine; load the mosaic; return (sci, header, wht, exptime).""" from astropy.io import fits + sci_path, wht_path, drizzle_prov = combine_mod.combine( + ctx.exposures, ctx.spec, ctx.adapter, ctx.work_dir + ) + ctx.record["drizzle"] = drizzle_prov + with fits.open(sci_path) as hdul: sci = hdul[0].data.astype(float) header = hdul[0].header.copy() @@ -111,14 +116,18 @@ def reduce_target( # Only the HST noise construction divides by exposure time; the JWST path # reads propagated ERR and needs no exptime — record it if present, but # never hard-fail a reduction that doesn't use it. - if adapter.combine_backend != "jwst_image3" and ( + if ctx.adapter.combine_backend != "jwst_image3" and ( exptime is None or exptime <= 0 ): raise ValueError(f"mosaic header carries no positive EXPTIME: {exptime}") - exptime = float(exptime) if exptime else 0.0 + return sci, header, wht, float(exptime) if exptime else 0.0 + + +def _noise(ctx: _StageContext, sci, wht, exptime: float) -> np.ndarray: + drizzle_prov = ctx.record["drizzle"] + if ctx.adapter.combine_backend == "jwst_image3": + from astropy.io import fits - # -- noise ----------------------------------------------------------------- - if adapter.combine_backend == "jwst_image3": from .noise import jwst_rms as jwst_rms_mod err = fits.getdata(drizzle_prov["err_path"]).astype(float) @@ -127,7 +136,7 @@ def reduce_target( sci, correlated_noise_factor=drizzle_prov["correlated_noise_factor"], ) - record["noise"] = { + ctx.record["noise"] = { "recipe": "R * ERR (propagated by calwebb_image3 resample)", "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], "exptime": float(exptime), @@ -140,7 +149,7 @@ def reduce_target( exptime=float(exptime), correlated_noise_factor=drizzle_prov["correlated_noise_factor"], ) - record["noise"] = { + ctx.record["noise"] = { "recipe": "R * sqrt(max(sci,0)/exptime + 1/wht)", "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], "exptime": float(exptime), @@ -148,10 +157,14 @@ def reduce_target( sci[np.isfinite(noise)] ), } + return noise + - # -- psf ------------------------------------------------------------------- +def _psf(ctx: _StageContext, sci, header): + from astropy.io import fits from astropy.wcs import WCS + spec, adapter = ctx.spec, ctx.adapter target_xy = WCS(header).world_to_pixel_values(spec.ra, spec.dec) selection = stars_mod.StarSelection() if adapter.observatory == "hst": @@ -160,7 +173,7 @@ def reduce_target( # full well by the longest single-exposure time — never the mosaic # total. max_single_exptime = max( - float(fits.getheader(p).get("EXPTIME", 0.0)) for p in exposures + float(fits.getheader(p).get("EXPTIME", 0.0)) for p in ctx.exposures ) if max_single_exptime <= 0.0: raise ValueError("no exposure carries a positive EXPTIME header") @@ -184,9 +197,14 @@ def reduce_target( psf, psf_full, psf_diag = epsf_mod.build_epsf( sci, stars, spec.psf_shape, spec.psf_full_shape ) - record["psf"] = psf_diag + ctx.record["psf"] = psf_diag + return psf, psf_full - # -- package ---------------------------------------------------------------- + +def _package(ctx: _StageContext, sci, header, wht, noise, psf, psf_full) -> None: + from astropy.io import fits + + spec, out_dir = ctx.spec, ctx.out_dir data_cut, data_header, center_xy = cutout_mod.make_cutout( sci, header, spec.ra, spec.dec, spec.cutout_shape ) @@ -205,16 +223,16 @@ def reduce_target( rms_mod.assert_finite_within(noise_cut, f"{spec.name} cutout") cutout_mod.write_fits(data_cut, data_header, out_dir / "data.fits") cutout_mod.write_fits(noise_cut, noise_header, out_dir / "noise_map.fits") - record["bad_pixel_policy"] = mask_diag + ctx.record["bad_pixel_policy"] = mask_diag # The mosaic-wide WHT uniformity mixes coverage tiers across the full # union footprint; the science verdict belongs to the cutout region. - from .drizzle.diagnostics import check_weight_uniformity - wht_cut, _, _ = cutout_mod.make_cutout( wht, header, spec.ra, spec.dec, spec.cutout_shape ) - record["drizzle"]["weight_uniformity_cutout"] = check_weight_uniformity(wht_cut) + ctx.record["drizzle"]["weight_uniformity_cutout"] = check_weight_uniformity( + wht_cut + ) fits.PrimaryHDU(psf.astype(np.float32)).writeto( out_dir / "psf.fits", overwrite=True @@ -222,7 +240,7 @@ def reduce_target( fits.PrimaryHDU(psf_full.astype(np.float32)).writeto( out_dir / "psf_full.fits", overwrite=True ) - record["package"] = { + ctx.record["package"] = { "products": ["data.fits", "noise_map.fits", "psf.fits", "psf_full.fits"], "cutout_shape": list(spec.cutout_shape), "pixel_scale": spec.final_scale, @@ -230,12 +248,45 @@ def reduce_target( "data_units": str(header.get("BUNIT", "unknown")), } - provenance_mod.write_reduction_json(out_dir, record) - # -- evict -------------------------------------------------------------------- - cache.mark_completed(spec.name) +def _evict(ctx: _StageContext, evict_when_done: bool) -> None: + ctx.cache.mark_completed(ctx.spec.name) if evict_when_done: - cache.evict(spec.name) - cache.enforce_cap() + ctx.cache.evict(ctx.spec.name) + ctx.cache.enforce_cap() + + +def reduce_target( + spec: TargetSpec, + cache_root: Path, + output_root: Path, + size_cap_bytes: Optional[int] = None, + evict_when_done: bool = False, +) -> Dict: + """Run the full pipeline for one target; returns the provenance record.""" + adapter = instruments.get(spec.instrument) + cache = cache_mod.ExposureCache(Path(cache_root), size_cap_bytes=size_cap_bytes) + out_dir = Path(output_root) / spec.name + out_dir.mkdir(parents=True, exist_ok=True) + work_dir = out_dir / "work" + work_dir.mkdir(exist_ok=True) + + ctx = _StageContext( + spec=spec, + adapter=adapter, + cache=cache, + out_dir=out_dir, + work_dir=work_dir, + record={"target": spec.as_dict(), "instrument": adapter.key}, + ) + + _acquire(ctx) + _align(ctx) + sci, header, wht, exptime = _combine(ctx) + noise = _noise(ctx, sci, wht, exptime) + psf, psf_full = _psf(ctx, sci, header) + _package(ctx, sci, header, wht, noise, psf, psf_full) + provenance_mod.write_reduction_json(out_dir, ctx.record) + _evict(ctx, evict_when_done) - return record + return ctx.record diff --git a/autoreduce/psf/stars.py b/autoreduce/psf/stars.py index 7680e34..4179b46 100644 --- a/autoreduce/psf/stars.py +++ b/autoreduce/psf/stars.py @@ -31,13 +31,10 @@ class StarSelection: def reject_crowded(x: np.ndarray, y: np.ndarray, min_separation: float) -> np.ndarray: """Boolean mask keeping sources with no neighbour within min_separation.""" - keep = np.ones(len(x), dtype=bool) - for i in range(len(x)): - d2 = (x - x[i]) ** 2 + (y - y[i]) ** 2 - d2[i] = np.inf - if (d2 < min_separation**2).any(): - keep[i] = False - return keep + n = len(x) + d2 = (x[:, None] - x[None, :]) ** 2 + (y[:, None] - y[None, :]) ** 2 + d2[np.arange(n), np.arange(n)] = np.inf + return ~(d2 < min_separation**2).any(axis=1) def reject_edges( diff --git a/autoreduce/validation/__init__.py b/autoreduce/validation/__init__.py new file mode 100644 index 0000000..b289d2e --- /dev/null +++ b/autoreduce/validation/__init__.py @@ -0,0 +1,7 @@ +""" +Validation utilities shared by the integration/acceptance scripts: sub-pixel +registration and reference-parity statistics (the SLACS-parity method, reused +verbatim for every instrument since phase 1). +""" + +from .parity import registered_ratios, subpixel_offset diff --git a/autoreduce/validation/parity.py b/autoreduce/validation/parity.py new file mode 100644 index 0000000..a65b2cf --- /dev/null +++ b/autoreduce/validation/parity.py @@ -0,0 +1,80 @@ +""" +Sub-pixel registration + parity statistics against a reference dataset. + +Extracted from the three per-instrument integration scripts (which had +triplicated it); behaviour matches the phase-3 version — bounds-guarded +parabolic peak refinement, bright-pixel data ratios, masked-pixel-excluded +noise ratios. +""" + +from typing import Dict, Tuple + +import numpy as np + +# Pixels at/above this noise value are masked-by-noise products (see +# noise.rms.MASKED_NOISE_VALUE) or their shift-interpolation bleed. +_MASKED_EXCLUSION_THRESHOLD = 1.0e6 + + +def subpixel_offset(a: np.ndarray, b: np.ndarray) -> Tuple[float, float]: + """(dy, dx) shift of `b` relative to `a`: FFT cross-correlation + parabolic peak.""" + a0 = np.nan_to_num(a - np.nanmedian(a)) + b0 = np.nan_to_num(b - np.nanmedian(b)) + corr = np.fft.fftshift( + np.fft.irfft2(np.fft.rfft2(a0) * np.conj(np.fft.rfft2(b0)), s=a0.shape) + ) + peak = np.unravel_index(np.argmax(corr), corr.shape) + + def parabolic(axis): + prev = list(peak); prev[axis] -= 1 + nxt = list(peak); nxt[axis] += 1 + if min(prev[axis], 0) < 0 or nxt[axis] >= corr.shape[axis]: + return 0.0 + cm, c0, cp = corr[tuple(prev)], corr[peak], corr[tuple(nxt)] + denom = cm - 2 * c0 + cp + return 0.0 if denom == 0 else 0.5 * (cm - cp) / denom + + return ( + peak[0] - a.shape[0] // 2 + parabolic(0), + peak[1] - a.shape[1] // 2 + parabolic(1), + ) + + +def registered_ratios( + new_data: np.ndarray, + new_noise: np.ndarray, + ref_data: np.ndarray, + ref_noise: np.ndarray, + bright_sigma: float = 10.0, +) -> Dict: + """Register `new` onto `ref` (sub-pixel) and report parity statistics.""" + from scipy.ndimage import shift as nd_shift + + if new_data.shape != ref_data.shape: + raise ValueError( + f"shape mismatch: new {new_data.shape} vs reference {ref_data.shape}" + ) + + dy, dx = subpixel_offset(ref_data, new_data) + new_data_r = nd_shift(np.nan_to_num(new_data), (dy, dx), order=3) + new_noise_r = nd_shift(np.nan_to_num(new_noise), (dy, dx), order=1) + + bright = ref_data > bright_sigma * np.nanmedian(ref_noise) + data_ratio = new_data_r[bright] / ref_data[bright] + # Exclude masked-by-noise pixels and their shift-interpolation bleed. + valid = new_noise_r < _MASKED_EXCLUSION_THRESHOLD + noise_ratio = np.where(valid, new_noise_r / ref_noise, np.nan) + return { + "offset": [float(dy), float(dx)], + "n_bright": int(bright.sum()), + "data_ratio_median": float(np.nanmedian(data_ratio)), + "data_ratio_16_84": [ + float(np.nanpercentile(data_ratio, 16)), + float(np.nanpercentile(data_ratio, 84)), + ], + "noise_ratio_median": float(np.nanmedian(noise_ratio)), + "noise_ratio_16_84": [ + float(np.nanpercentile(noise_ratio, 16)), + float(np.nanpercentile(noise_ratio, 84)), + ], + } diff --git a/scripts/reduce_cosmos_web_ring.py b/scripts/reduce_cosmos_web_ring.py index 81220b1..7b21aa5 100644 --- a/scripts/reduce_cosmos_web_ring.py +++ b/scripts/reduce_cosmos_web_ring.py @@ -16,8 +16,6 @@ import sys from pathlib import Path -import numpy as np - REPO = Path(__file__).resolve().parent.parent sys.path.insert(0, str(REPO)) @@ -56,64 +54,19 @@ def spec_for(band: str) -> TargetSpec: ) -def subpixel_offset(a, b): - a0 = np.nan_to_num(a - np.nanmedian(a)) - b0 = np.nan_to_num(b - np.nanmedian(b)) - corr = np.fft.fftshift( - np.fft.irfft2(np.fft.rfft2(a0) * np.conj(np.fft.rfft2(b0)), s=a0.shape) - ) - peak = np.unravel_index(np.argmax(corr), corr.shape) - - def parabolic(axis): - prev = list(peak); prev[axis] -= 1 - nxt = list(peak); nxt[axis] += 1 - if min(prev[axis], 0) < 0 or nxt[axis] >= corr.shape[axis]: - return 0.0 - cm, c0, cp = corr[tuple(prev)], corr[peak], corr[tuple(nxt)] - denom = cm - 2 * c0 + cp - return 0.0 if denom == 0 else 0.5 * (cm - cp) / denom - - return ( - peak[0] - a.shape[0] // 2 + parabolic(0), - peak[1] - a.shape[1] // 2 + parabolic(1), - ) - - def compare(band: str, out_dir: Path) -> dict: from astropy.io import fits - from scipy.ndimage import shift as nd_shift + + from autoreduce.validation import registered_ratios new_data = fits.getdata(out_dir / "data.fits").astype(float) new_noise = fits.getdata(out_dir / "noise_map.fits").astype(float) demo_dir = DEMO_ROOT / band demo_data = fits.getdata(demo_dir / "data.fits").astype(float) demo_noise = fits.getdata(demo_dir / "noise_map.fits").astype(float) - - if new_data.shape != demo_data.shape: - raise ValueError( - f"{band}: shape mismatch new {new_data.shape} vs demo {demo_data.shape}" - ) - - dy, dx = subpixel_offset(demo_data, new_data) - new_data_r = nd_shift(np.nan_to_num(new_data), (dy, dx), order=3) - new_noise_r = nd_shift(np.nan_to_num(new_noise), (dy, dx), order=1) - - bright = demo_data > 10 * np.nanmedian(demo_noise) - data_ratio = new_data_r[bright] / demo_data[bright] - # Exclude masked-by-noise pixels (1e8) and their shift-interpolation - # bleed from the parity statistics. - valid = new_noise_r < 1.0e6 - noise_ratio = np.where(valid, new_noise_r / demo_noise, np.nan) return { "band": band, - "offset": [float(dy), float(dx)], - "n_bright": int(bright.sum()), - "data_ratio_median": float(np.nanmedian(data_ratio)), - "noise_ratio_median": float(np.nanmedian(noise_ratio)), - "noise_ratio_16_84": [ - float(np.nanpercentile(noise_ratio, 16)), - float(np.nanpercentile(noise_ratio, 84)), - ], + **registered_ratios(new_data, new_noise, demo_data, demo_noise), } diff --git a/scripts/reduce_slacs0008.py b/scripts/reduce_slacs0008.py index 9254798..6e2a723 100644 --- a/scripts/reduce_slacs0008.py +++ b/scripts/reduce_slacs0008.py @@ -31,33 +31,13 @@ ) -def subpixel_offset(a: np.ndarray, b: np.ndarray): - """(dy, dx) shift of b relative to a: FFT cross-correlation + parabolic peak.""" - a0 = np.nan_to_num(a - np.nanmedian(a)) - b0 = np.nan_to_num(b - np.nanmedian(b)) - corr = np.fft.irfft2(np.fft.rfft2(a0) * np.conj(np.fft.rfft2(b0)), s=a0.shape) - corr = np.fft.fftshift(corr) - peak = np.unravel_index(np.argmax(corr), corr.shape) - - def parabolic(idx, axis): - c0 = corr[peak] - prev = list(peak); prev[axis] -= 1 - nxt = list(peak); nxt[axis] += 1 - cm, cp = corr[tuple(prev)], corr[tuple(nxt)] - denom = cm - 2 * c0 + cp - return 0.0 if denom == 0 else 0.5 * (cm - cp) / denom - - dy = peak[0] - a.shape[0] // 2 + parabolic(peak, 0) - dx = peak[1] - a.shape[1] // 2 + parabolic(peak, 1) - return dy, dx - - def main(): record = reduce_target(SPEC, cache_root=CACHE_ROOT, output_root=OUTPUT_ROOT) print(json.dumps(record["drizzle"]["weight_uniformity"], indent=2)) from astropy.io import fits - from scipy.ndimage import shift as nd_shift + + from autoreduce.validation import registered_ratios out_dir = OUTPUT_ROOT / SPEC.name new_data = fits.getdata(out_dir / "data.fits").astype(float) @@ -65,29 +45,9 @@ def main(): legacy_data = fits.getdata(LEGACY_DIR / "data.fits").astype(float) legacy_noise = fits.getdata(LEGACY_DIR / "noise_map.fits").astype(float) - dy, dx = subpixel_offset(legacy_data, new_data) - print(f"[parity] sub-pixel offset (legacy vs new): dy={dy:.3f}, dx={dx:.3f}") - new_data_r = nd_shift(np.nan_to_num(new_data), (dy, dx), order=3) - new_noise_r = nd_shift(np.nan_to_num(new_noise), (dy, dx), order=1) - - bright = legacy_data > 10 * np.nanmedian(legacy_noise) - data_ratio = new_data_r[bright] / legacy_data[bright] - noise_ratio = new_noise_r / legacy_noise - summary = { "n_exposures": record["acquire"]["n_exposures"], - "offset": [float(dy), float(dx)], - "n_bright": int(bright.sum()), - "data_ratio_median": float(np.nanmedian(data_ratio)), - "data_ratio_16_84": [ - float(np.nanpercentile(data_ratio, 16)), - float(np.nanpercentile(data_ratio, 84)), - ], - "noise_ratio_median": float(np.nanmedian(noise_ratio)), - "noise_ratio_16_84": [ - float(np.nanpercentile(noise_ratio, 16)), - float(np.nanpercentile(noise_ratio, 84)), - ], + **registered_ratios(new_data, new_noise, legacy_data, legacy_noise), "correlated_noise_factor_applied": record["noise"]["correlated_noise_factor"], "psf_diagnostics": record["psf"], } diff --git a/test_autoreduce/test_psf_and_package.py b/test_autoreduce/test_psf_and_package.py index d196615..0f48fc3 100644 --- a/test_autoreduce/test_psf_and_package.py +++ b/test_autoreduce/test_psf_and_package.py @@ -145,3 +145,57 @@ def test_mast_query_hygiene(): assert is_direct_observation("j9op01010", "10886") assert not is_direct_observation("hst_skycell-p1322x03y02_acs_wfc_f814w_all", "--") assert not is_direct_observation("j9op01010", "--") + + +def test_reject_crowded_matches_reference_loop(): + """Randomized equivalence vs the original O(N^2) loop implementation.""" + from autoreduce.psf.stars import reject_crowded + + def reference(x, y, min_separation): + keep = np.ones(len(x), dtype=bool) + for i in range(len(x)): + d2 = (x - x[i]) ** 2 + (y - y[i]) ** 2 + d2[i] = np.inf + if (d2 < min_separation**2).any(): + keep[i] = False + return keep + + rng = np.random.default_rng(7) + for n in (0, 1, 2, 50, 300): + x = rng.uniform(0, 500, n) + y = rng.uniform(0, 500, n) + for sep in (1.0, 25.0, 100.0): + assert ( + reject_crowded(x, y, sep) == reference(x, y, sep) + ).all(), (n, sep) + + +def test_registered_ratios_recovers_known_shift_and_scale(): + from scipy.ndimage import shift as nd_shift + + from autoreduce.validation import registered_ratios + + rng = np.random.default_rng(3) + ref_data = rng.normal(0.0, 0.01, (120, 120)) + yy, xx = np.mgrid[0:120, 0:120] + ref_data += 8.0 * np.exp(-(((xx - 60) ** 2 + (yy - 60) ** 2) / (2 * 3.0**2))) + ref_noise = np.full((120, 120), 0.01) + + new_data = 1.5 * nd_shift(ref_data, (1.25, -0.75), order=3) + new_noise = 2.0 * ref_noise + out = registered_ratios(new_data, new_noise, ref_data, ref_noise) + # The offset is the shift applied to `new` to register it onto `ref` — + # the negative of new's displacement. + assert out["offset"][0] == pytest.approx(-1.25, abs=0.15) + assert out["offset"][1] == pytest.approx(0.75, abs=0.15) + assert out["data_ratio_median"] == pytest.approx(1.5, rel=0.05) + assert out["noise_ratio_median"] == pytest.approx(2.0, rel=0.05) + + # Masked-by-noise pixels are excluded from the noise statistics. + new_noise_masked = new_noise.copy() + new_noise_masked[5, 5] = 1.0e8 + out2 = registered_ratios(new_data, new_noise_masked, ref_data, ref_noise) + assert out2["noise_ratio_median"] == pytest.approx(2.0, rel=0.05) + + with pytest.raises(ValueError, match="shape mismatch"): + registered_ratios(new_data[:100], new_noise[:100], ref_data, ref_noise) From 493bc30a733d65885f26fdf5b0a8df73c31125cf Mon Sep 17 00:00:00 2001 From: Jammy2211 Date: Thu, 9 Jul 2026 12:12:55 +0100 Subject: [PATCH 4/4] Review cosmetics: dead imports pruned, edge-margin comment explains headroom Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01WNhuuvZLgZjHvWBGdhUst3 --- autoreduce/drizzle/jwst_combine.py | 1 - autoreduce/psf/stars.py | 4 ++-- scripts/reduce_slacs0008.py | 1 - 3 files changed, 2 insertions(+), 4 deletions(-) diff --git a/autoreduce/drizzle/jwst_combine.py b/autoreduce/drizzle/jwst_combine.py index 73cc4ac..a1cd760 100644 --- a/autoreduce/drizzle/jwst_combine.py +++ b/autoreduce/drizzle/jwst_combine.py @@ -10,7 +10,6 @@ stage (noise, psf, package) is backend-agnostic. """ -import json from pathlib import Path from typing import Dict, List, Tuple diff --git a/autoreduce/psf/stars.py b/autoreduce/psf/stars.py index 4179b46..3214022 100644 --- a/autoreduce/psf/stars.py +++ b/autoreduce/psf/stars.py @@ -23,8 +23,8 @@ class StarSelection: round_limit: float = 0.3 saturation_fraction: float = 0.7 # of adapter.saturation_dn, in counts min_separation_pix: float = 25.0 - # Must exceed half the ePSF extraction window (psf_full 61 + 20 pad -> 41) - # or edge stars pass selection only to be dropped at extraction. + # Must exceed half the ePSF extraction window (psf_full 61 + 20 pad -> + # half-window 41), plus headroom for EPSFBuilder recentering drift. edge_margin_pix: int = 46 exclusion_radius_pix: float = 50.0 # around the target itself diff --git a/scripts/reduce_slacs0008.py b/scripts/reduce_slacs0008.py index 6e2a723..08539e1 100644 --- a/scripts/reduce_slacs0008.py +++ b/scripts/reduce_slacs0008.py @@ -12,7 +12,6 @@ import sys from pathlib import Path -import numpy as np REPO = Path(__file__).resolve().parent.parent sys.path.insert(0, str(REPO))