From 2168738e6bd95e6caf4e0586f78bfe76c5358d88 Mon Sep 17 00:00:00 2001 From: Claude Date: Thu, 6 Aug 2026 13:42:13 +0000 Subject: [PATCH] Add cr_method dial (per-frame deepCR route) and PSF star-pass decoupling Community-reported (#61, #62, @samlange04): driz_cr's blotted-median reference reads systematically low on steep gradients, flagging genuine deflector-core flux as cosmic rays (~37% core flux loss measured on SLACS ACS/WFC F814W), and PSF star finding on that same CR-rejected mosaic holes star cores before DAOStarFinder sees them (344 -> 599 usable stars when rebuilt from a no-CR pass). - TargetSpec.cr_method ("driz_cr" default | "deepcr"): the per-frame route writes deepCR masks (the machinery frame_products already uses -- #61's reporter used LACosmic; deepCR is the same per-frame, no-median-reference class without a new dependency) into each exposure's DQ as the AstroDrizzle CR bit (4096), then plain weighted-mean drizzles with median=blot=driz_cr=False and resetbits=0. resetbits=0 is #61's trap -- AstroDrizzle's default resetbits=4096 clears exactly the bit the masks were written into -- and is pinned by unit test, along with a regression test that the default driz_cr route and the single-exposure branch are unchanged. The default deliberately stays "driz_cr": the flip is human-gated on SLACS validation with a tuned-driz_cr comparison arm (driz_cr_scale/snr 1.5/1.2 per STScI guidance), so #61/#62 stay open. - TargetSpec.psf_star_pass ("auto" default | "science" | "no_cr"): Tier-1/1b star finding + stamp extraction are decoupled from the shipped science mosaic. "auto" never adds a drizzle (the science mosaic is recorded, with the reason it is or is not the least-CR-rejected pass); "no_cr" is the explicit opt-in second AstroDrizzle pass that treats CR DQ flags as good (final_bits | 4096) without clearing them (resetbits=0), preserving the science flags for frame products. The reduction.json psf block records star_source_pass in every case so the coupling cannot silently regress. - Unsupported combinations (non-astrodrizzle backends, instruments without a deepCR model, no_cr + psf_from_frames) fail fast in reduce_target before any download; both dials are documented in docs/design/hst_acs_pipeline.md stages 3 and 5 as deviations-with- justification, default unchanged. Refs #61, #62 (both stay open pending the SLACS validation arms). Co-Authored-By: Claude Fable 5 Claude-Session: https://claude.ai/code/session_01UBpwwfzbewbMLhE137h231 --- autoreduce/drizzle/combine.py | 131 ++++++++++++++++++ autoreduce/pipeline.py | 88 +++++++++++- autoreduce/target.py | 35 +++++ docs/design/hst_acs_pipeline.md | 36 +++++ test_autoreduce/test_psf_and_package.py | 123 ++++++++++++++++ .../test_target_and_instruments.py | 100 +++++++++++++ 6 files changed, 508 insertions(+), 5 deletions(-) diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py index 3318fad..cc9ff7a 100644 --- a/autoreduce/drizzle/combine.py +++ b/autoreduce/drizzle/combine.py @@ -13,9 +13,16 @@ from pathlib import Path from typing import Dict, List, Tuple +import numpy as np + from ..instruments import InstrumentAdapter from ..target import TargetSpec +# AstroDrizzle's cosmic-ray DQ bit. The cr_method="deepcr" route writes its +# per-frame masks into this same bit, so every downstream consumer of CR +# flags (final drizzle bit masking, frame products) reads one convention. +CR_DQ_BIT = 4096 + def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int) -> Dict: """Assemble the AstroDrizzle keyword set; pure function, unit-testable.""" @@ -36,9 +43,87 @@ def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures median=multi, blot=multi, ) + if multi and spec.cr_method == "deepcr": + # Per-frame route (issue #61): CR masks are already in the DQ arrays + # (apply_per_frame_cr_masks), so the stack rejection is off and the + # combine is a plain weighted mean. resetbits=0 is load-bearing — + # AstroDrizzle's default resetbits (4096) clears exactly the DQ bit + # the masks were written into, silently producing an unmasked mosaic. + kwargs.update(driz_cr=False, median=False, blot=False, resetbits=0) return kwargs +def star_pass_kwargs_for( + spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int +) -> Dict: + """ + Keyword set for the dedicated PSF-star drizzle pass + (``psf_star_pass="no_cr"``, issue #62): the science kwargs with every + stack CR-rejection step off and previously written CR DQ flags treated + as good (``final_bits`` gains the CR bit). ``resetbits=0`` so the + science pass's DQ flags survive in the input exposures — frame products + and re-runs still see them. Pure function, unit-testable. + """ + kwargs = drizzle_kwargs_for(spec, adapter, n_exposures) + kwargs.update( + driz_cr=False, + median=False, + blot=False, + resetbits=0, + final_bits=int(kwargs.get("final_bits", 0)) | CR_DQ_BIT, + ) + return kwargs + + +def dq_with_cr_flags(dq: np.ndarray, cr_mask: np.ndarray, cr_bit: int = CR_DQ_BIT) -> np.ndarray: + """ + One chip's DQ array with the CR bit rewritten from ``cr_mask``: the bit + is cleared everywhere first (so re-runs are idempotent, mirroring + AstroDrizzle's own resetbits behaviour) then set where the mask flags a + cosmic ray. Pure function; preserves the DQ dtype. + """ + dq = np.asarray(dq) + cr_mask = np.asarray(cr_mask, dtype=bool) + if dq.shape != cr_mask.shape: + raise ValueError(f"DQ/mask shape mismatch: {dq.shape} vs {cr_mask.shape}") + cleared = (dq & ~cr_bit).astype(dq.dtype) + return np.where(cr_mask, cleared | cr_bit, cleared).astype(dq.dtype) + + +def apply_per_frame_cr_masks(exposures: List[Path], adapter: InstrumentAdapter) -> Dict: + """ + The ``cr_method="deepcr"`` masking step (issue #61): per-frame deepCR + cosmic-ray masks written into each exposure's DQ arrays as the + AstroDrizzle CR bit, replacing driz_cr's blotted-median stack rejection + (which reads systematically low on steep gradients and flags genuine + core flux). Mutates the exposures' DQ in place — exactly as driz_cr + itself does — idempotently (the CR bit is cleared before rewriting). + Returns the provenance fragment. + """ + from astropy.io import fits + + from ..package import cosmic_rays as cr_mod + + masker = cr_mod.masker_for(adapter.key) + n_cr_pixels = {} + for path in exposures: + with fits.open(path, mode="update") as hdul: + total = 0 + for hdu in hdul: + if hdu.name != "SCI": + continue + mask = masker(hdu.data) + dq_hdu = hdul["DQ", hdu.ver] + dq_hdu.data = dq_with_cr_flags(dq_hdu.data, mask) + total += int(mask.sum()) + n_cr_pixels[Path(path).name] = total + return { + **cr_mod.cr_method_record(adapter.key), + "cr_bit": CR_DQ_BIT, + "n_cr_pixels": n_cr_pixels, + } + + def combine( exposures: List[Path], spec: TargetSpec, @@ -75,6 +160,9 @@ def combine( # 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() + per_frame_cr = None + if len(exposures) > 1 and spec.cr_method == "deepcr": + per_frame_cr = apply_per_frame_cr_masks(exposures, adapter) kwargs = drizzle_kwargs_for(spec, adapter, len(exposures)) with chdir_scratch(output_dir) as output_dir: astrodrizzle.AstroDrizzle( @@ -95,6 +183,9 @@ def _one(suffix: str) -> Path: sci = _one("_sci.fits") wht = _one("_wht.fits") + tail = {"cr_method": spec.cr_method} + if per_frame_cr is not None: + tail["per_frame_cr"] = per_frame_cr provenance = combine_provenance( spec, adapter, @@ -102,5 +193,45 @@ def _one(suffix: str) -> Path: fits.getdata(wht), kwargs_key="drizzle_kwargs", kwargs={k: kwargs[k] for k in sorted(kwargs)}, + tail=tail, ) return sci, wht, provenance + + +def combine_star_pass( + exposures: List[Path], + spec: TargetSpec, + adapter: InstrumentAdapter, + output_dir: Path, +) -> Tuple[Path, Dict]: + """ + Drizzle the dedicated CR-flag-ignoring PSF-star mosaic + (``psf_star_pass="no_cr"``, issue #62) onto the same grid/geometry as + the science combine; returns (sci_path, kwargs). AstroDrizzle backend + only — the caller gates on ``adapter.combine_backend``. + """ + if adapter.combine_backend != "astrodrizzle": + raise ValueError( + "combine_star_pass supports the astrodrizzle backend only " + f"(instrument {adapter.key!r})" + ) + + from drizzlepac import astrodrizzle + + from ._common import chdir_scratch + + output_name = f"{spec.name}_{spec.filter_name}_starpass".lower() + kwargs = star_pass_kwargs_for(spec, adapter, len(exposures)) + with chdir_scratch(output_dir) as output_dir: + astrodrizzle.AstroDrizzle( + input=[str(p) for p in exposures], + output=output_name, + **kwargs, + ) + hits = sorted(glob.glob(f"{output_dir / output_name}*_sci.fits")) + if len(hits) != 1: + raise FileNotFoundError( + f"expected exactly one star-pass _sci.fits for " + f"{output_dir / output_name}, got {hits}" + ) + return Path(hits[0]), kwargs diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index 4a6078c..6ef1e86 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -25,6 +25,7 @@ from .drizzle.diagnostics import check_weight_uniformity from .instruments import InstrumentAdapter from .noise import rms as rms_mod +from .package import cosmic_rays as cr_mod from .package import cutout as cutout_mod from .package import frames as frames_mod from .package import provenance as provenance_mod @@ -481,6 +482,57 @@ def _noise(ctx: _StageContext, sci, wht, exptime: float) -> np.ndarray: return noise +def _star_pass_image(ctx: _StageContext, sci, header): + """ + The mosaic that PSF star finding + stamp extraction run on (issue #62), + decoupled from the shipped science mosaic: driz_cr's blotted-median + rejection reads systematically low on star cores, holing PSF stars + before DAOStarFinder sees them. Returns (image, header, provenance) — + provenance always records which pass fed the stars, so the coupling + cannot silently regress. + + ``psf_star_pass="auto"`` (default) never adds a drizzle: it uses the + science mosaic, which is already the least-CR-rejected pass at zero + cost when no stack rejection ran (single-exposure branch, per-frame + cr_method, non-astrodrizzle backends) — and otherwise is simply the + only pass available without doubling combine time. The explicit + ``"no_cr"`` opt-in buys the dedicated CR-flag-ignoring pass with a + second full AstroDrizzle run. + """ + spec, adapter = ctx.spec, ctx.adapter + single = bool(ctx.record["drizzle"].get("single_exposure_branch")) + if spec.psf_star_pass == "no_cr" and not single: + from astropy.io import fits + + sci_path, kwargs = combine_mod.combine_star_pass( + ctx.exposures, spec, adapter, ctx.work_dir + ) + with fits.open(sci_path) as hdul: + star_sci = hdul[0].data.astype(float) + star_header = hdul[0].header.copy() + return star_sci, star_header, { + "star_source_pass": "no_cr_drizzle", + "star_pass_kwargs": {k: kwargs[k] for k in sorted(kwargs)}, + } + if single: + reason = "single-exposure branch: no CR rejection ran" + elif getattr(adapter, "combine_backend", None) != "astrodrizzle": + reason = ( + f"combine backend {adapter.combine_backend!r}: the dedicated " + "star pass applies to the astrodrizzle path only" + ) + elif spec.cr_method != "driz_cr": + reason = f"cr_method={spec.cr_method!r}: per-frame CR masks, no stack rejection" + elif spec.psf_star_pass == "science": + reason = "psf_star_pass='science': pinned to the shipped mosaic" + else: + reason = ( + "psf_star_pass='auto': no cheaper less-CR-rejected pass exists; " + "opt into 'no_cr' for a dedicated star drizzle" + ) + return sci, header, {"star_source_pass": "science", "star_source_reason": reason} + + def _psf(ctx: _StageContext, sci, header, noise=None): from astropy.io import fits from astropy.wcs import WCS @@ -521,7 +573,8 @@ def _psf(ctx: _StageContext, sci, header, noise=None): ) ctx.record["psf"] = diag return psf, psf_full - target_xy = WCS(header).world_to_pixel_values(spec.ra, spec.dec) + star_sci, star_header, star_pass_prov = _star_pass_image(ctx, sci, header) + target_xy = WCS(star_header).world_to_pixel_values(spec.ra, spec.dec) selection = stars_mod.StarSelection() max_single_exptime = adapter.max_single_exposure_seconds( [fits.getheader(p) for p in ctx.exposures] @@ -532,7 +585,7 @@ def _psf(ctx: _StageContext, sci, header, noise=None): else None ) stars = stars_mod.find_stars( - sci, + star_sci, selection, target_xy=(float(target_xy[0]), float(target_xy[1])), peak_max=peak_max, @@ -548,18 +601,21 @@ def _psf(ctx: _StageContext, sci, header, noise=None): "psf_backend='starred' needs the noise map (STARRED weights " "stars by per-pixel noise); pipeline did not pass it" ) + # Stamps come from the star pass; the science-pass noise map is a + # faithful weight there (same exposures/geometry — the passes differ + # only in which pixels CR rejection removed). psf, psf_full, psf_diag = starred_mod.build_starred_epsf( - sci, noise, stars, spec.psf_shape, spec.psf_full_shape + star_sci, noise, stars, spec.psf_shape, spec.psf_full_shape ) else: psf, psf_full, psf_diag = epsf_mod.build_epsf( - sci, stars, spec.psf_shape, spec.psf_full_shape + star_sci, stars, spec.psf_shape, spec.psf_full_shape ) if adapter.observatory == "keck": # Tier B: in-field ePSF. Usable, but an AO PSF from field stars at a # different anisoplanatic angle is still provisional by contract. psf_diag = {"psf_provisional": True, **psf_diag} - ctx.record["psf"] = psf_diag + ctx.record["psf"] = {**star_pass_prov, **psf_diag} return psf, psf_full @@ -705,6 +761,28 @@ def reduce_target( "psf_from_frames supports HST and JWST only " f"(instrument {spec.instrument!r})" ) + backend = getattr(adapter, "combine_backend", None) + if spec.cr_method != "driz_cr" and backend != "astrodrizzle": + raise ValueError( + f"cr_method={spec.cr_method!r} supports the HST astrodrizzle " + f"path only (instrument {spec.instrument!r})" + ) + if spec.cr_method == "deepcr" and spec.instrument not in cr_mod.DEEPCR_MODELS: + raise ValueError( + f"cr_method='deepcr' needs a registered deepCR model; none for " + f"{spec.instrument!r} (known: {sorted(cr_mod.DEEPCR_MODELS)}; " + "wfc3_ir cosmic rays are already per-frame flagged by calwf3 " + "ramp fitting)" + ) + if spec.psf_star_pass == "no_cr" and ( + backend != "astrodrizzle" or spec.psf_from_frames + ): + raise ValueError( + "psf_star_pass='no_cr' builds a second AstroDrizzle pass for " + "mosaic star finding — HST astrodrizzle path without " + f"psf_from_frames only (instrument {spec.instrument!r}, " + f"psf_from_frames={spec.psf_from_frames})" + ) if spec.inject_image and adapter.domain != "visibility" and ( observatory, getattr(adapter, "combine_backend", None) ) not in ( diff --git a/autoreduce/target.py b/autoreduce/target.py index 040df57..be3ae95 100644 --- a/autoreduce/target.py +++ b/autoreduce/target.py @@ -42,6 +42,20 @@ class TargetSpec: final_pixfrac: float = 0.8 final_kernel: str = "square" + # Cosmic-ray rejection route for multi-exposure AstroDrizzle combines + # (issue #61). "driz_cr" (STScI default): blotted-median reference + + # driz_cr — on steep gradients (galaxy cores, stars) the sub-pixel- + # dithered median reads systematically low, so genuine core flux can be + # flagged as CR. "deepcr": per-frame deepCR masks written into each + # exposure's DQ arrays (the AstroDrizzle CR bit, 4096), then a plain + # weighted-mean drizzle with median/blot/driz_cr off and resetbits=0 + # (the default resetbits, 4096, would clear exactly the bit the masks + # were written into). #61's reporter used LACosmic; the pipeline's + # established per-frame CR machinery is deepCR (package/cosmic_rays.py), + # reused here instead of adding a dependency. Default deliberately + # unchanged pending SLACS validation (see hst_acs_pipeline.md stage 3). + cr_method: str = "driz_cr" + # PSF products (design doc stage 5). psf_shape: Tuple[int, int] = (21, 21) psf_full_shape: Tuple[int, int] = (61, 61) @@ -64,6 +78,18 @@ class TargetSpec: # when psf_from_frames or the Keck tier-A path own the PSF stage. psf_backend: str = "epsf" + # Which drizzle pass feeds mosaic PSF star finding + stamp extraction + # (issue #62): driz_cr's core rejection can hole PSF-star cores before + # DAOStarFinder sees them, so the star pass is decoupled from the + # shipped science mosaic. "auto" (default): the least-CR-rejected pass + # that costs no extra drizzle — today that is always the science mosaic + # (recorded with the reason); "science": pin star finding to the shipped + # mosaic; "no_cr": build a dedicated CR-flag-ignoring drizzle for the + # stars — an explicit opt-in, because it is a second full AstroDrizzle + # pass (HST astrodrizzle path only). Provenance records which pass fed + # the stars either way. + psf_star_pass: str = "auto" + # Alignment: residual (pixels) above which TweakReg refinement triggers. alignment_tolerance_pix: float = 0.1 @@ -135,6 +161,15 @@ def __post_init__(self): raise ValueError(f"dec out of range: {self.dec}") if not 0.0 < self.final_pixfrac <= 1.0: raise ValueError(f"final_pixfrac must be in (0, 1]: {self.final_pixfrac}") + if self.cr_method not in ("driz_cr", "deepcr"): + raise ValueError( + f"cr_method must be 'driz_cr' or 'deepcr': {self.cr_method!r}" + ) + if self.psf_star_pass not in ("auto", "science", "no_cr"): + raise ValueError( + "psf_star_pass must be 'auto', 'science' or 'no_cr': " + f"{self.psf_star_pass!r}" + ) if self.alma_width < 0: raise ValueError( f"alma_width must be >= 0 (0 = collapse the spw): {self.alma_width}" diff --git a/docs/design/hst_acs_pipeline.md b/docs/design/hst_acs_pipeline.md index 9fd19c4..c869f19 100644 --- a/docs/design/hst_acs_pipeline.md +++ b/docs/design/hst_acs_pipeline.md @@ -114,6 +114,27 @@ median-combine baseline), final drizzle of all exposures onto one grid. Undrizzled artifacts (`_single_sci`, masks) stay in the transient cache; only the mosaic + weight map proceed. +**Cosmic-ray rejection — `TargetSpec.cr_method` (issue #61).** The default +stays the STScI flow above: `driz_cr` against the blotted-median stack. On +steep gradients (deflector cores, PSF stars) that median reference reads +systematically low — sub-pixel dither shifts smear the peak — so driz_cr can +flag genuine core flux as cosmic rays; #61 measured ~37% deflector-core flux +loss on SLACS ACS/WFC F814W at the pipeline thresholds, worse with looser +ones. `cr_method="deepcr"` is the documented per-frame alternative: deepCR +masks — the per-frame CR machinery frame products already use; #61's reporter +used L.A.Cosmic (van Dokkum 2001, per-frame Laplacian), and deepCR is the +same no-median-reference class without adding a dependency — are written into +each exposure's DQ as the AstroDrizzle CR bit (4096), followed by a plain +weighted-mean drizzle (`median=blot=driz_cr=False`) with **`resetbits=0`**. +That last keyword is load-bearing: AstroDrizzle's default `resetbits=4096` +clears exactly the DQ bit the masks were written into, silently producing an +unmasked mosaic that can score flawlessly against itself (#61's trap, pinned +by unit test). **The default flip is deliberately human-gated on SLACS +validation** with two comparison arms: the per-frame route vs a *tuned* +driz_cr (raising `driz_cr_scale`/`driz_cr_snr` per STScI guidance — published +reprocessing used 1.5/1.2 against the pipeline default 1.2/0.7). Until that +lands, `driz_cr` remains the default and the dial is opt-in. + **The SLACS-V caveat (literature finding):** [Bolton et al. 2008](https://arxiv.org/abs/0805.1931) did **not** drizzle the F814W snapshot data — "the 'drizzle' re-sampling algorithm … is not well suited to @@ -206,6 +227,21 @@ mosaic. Never pair a native-frame PSF with a drizzled image. under-concentrated — and the photutils ePSF wins there (consistent with the #35 adversarial undersampling result). Prefer STARRED for well-sampled and crowded/few-star fields; keep photutils (or Tier 2) for undersampled SW. +- **Star-pass decoupling — `TargetSpec.psf_star_pass` (issue #62):** Tier-1/1b + star finding *and* stamp extraction draw from the least-CR-rejected drizzle + pass, decoupled from the shipped science mosaic: driz_cr's blotted-median + rejection holes star cores before DAOStarFinder sees them (#62 measured + 344→599 usable stars, +74%, rebuilding from a no-CR pass, rescuing 4 + lens/filter pairs from the model-PSF fallback). `"auto"` (default) never + adds a drizzle — the science mosaic is used and the choice + reason + recorded (it is genuinely the least-rejected pass on the single-exposure + branch and the per-frame `cr_method="deepcr"` route); `"no_cr"` is the + explicit opt-in that drizzles a second, CR-flag-ignoring star pass + (`final_bits | 4096`; `resetbits=0` so the science pass's DQ flags survive + for frame products) onto the same grid — same kernel/pixfrac/scale/rot, so + the drizzled-PSF invariant holds; `"science"` pins star finding to the + shipped mosaic. The `reduction.json` `psf` block records + `star_source_pass` in every case, so the coupling cannot silently regress. - **Tier 2 — TinyTim + focus model** (fallback; SLACS elliptical snapshot fields are typically star-poor): model PSFs raytraced per exposure with TinyTim, focus ("breathing") estimated by matching whatever stars exist, diff --git a/test_autoreduce/test_psf_and_package.py b/test_autoreduce/test_psf_and_package.py index 7d17dcb..8079772 100644 --- a/test_autoreduce/test_psf_and_package.py +++ b/test_autoreduce/test_psf_and_package.py @@ -283,3 +283,126 @@ def test_registered_ratios_recovers_known_shift_and_scale(): with pytest.raises(ValueError, match="shape mismatch"): registered_ratios(new_data[:100], new_noise[:100], ref_data, ref_noise) + + +class TestStarPassDecoupling: + """_psf's star-source pass selection (issue #62): star finding is + decoupled from the shipped science mosaic, and provenance always + records which pass fed the stars.""" + + def _ctx(self, spec, single_exposure, work_dir=None): + from pathlib import Path + + from autoreduce import instruments + from autoreduce.pipeline import _StageContext + + ctx = _StageContext( + spec=spec, + adapter=instruments.get("acs_wfc"), + cache=None, + out_dir=Path("."), + work_dir=work_dir or Path("."), + ) + ctx.record["drizzle"] = {"single_exposure_branch": single_exposure} + return ctx + + def _spec(self, **overrides): + from autoreduce.target import TargetSpec + + return TargetSpec(name="lens", ra=2.0, dec=-0.1, **overrides) + + def test_auto_default_uses_science_mosaic_and_records_why(self): + from autoreduce.pipeline import _star_pass_image + + ctx = self._ctx(self._spec(), single_exposure=False) + sci, header = np.zeros((5, 5)), object() + star_sci, star_header, prov = _star_pass_image(ctx, sci, header) + assert star_sci is sci and star_header is header + assert prov["star_source_pass"] == "science" + assert "no_cr" in prov["star_source_reason"] + + def test_single_exposure_branch_is_already_least_rejected(self): + from autoreduce.pipeline import _star_pass_image + + ctx = self._ctx(self._spec(psf_star_pass="no_cr"), single_exposure=True) + sci, header = np.zeros((5, 5)), object() + star_sci, _, prov = _star_pass_image(ctx, sci, header) + # "no_cr" on the single-exposure branch never drizzles again: the + # science mosaic had no CR rejection to begin with. + assert star_sci is sci + assert prov["star_source_pass"] == "science" + assert "single-exposure" in prov["star_source_reason"] + + def test_deepcr_route_records_per_frame_reason(self): + from autoreduce.pipeline import _star_pass_image + + ctx = self._ctx(self._spec(cr_method="deepcr"), single_exposure=False) + _, _, prov = _star_pass_image(ctx, np.zeros((5, 5)), object()) + assert prov["star_source_pass"] == "science" + assert "deepcr" in prov["star_source_reason"] + + def test_no_cr_opt_in_builds_the_dedicated_pass(self, monkeypatch, tmp_path): + from astropy.io import fits + + from autoreduce import pipeline as pipeline_mod + from autoreduce.pipeline import _star_pass_image + + star_path = tmp_path / "lens_f814w_starpass_sci.fits" + fits.PrimaryHDU(np.full((4, 4), 7.0, dtype=np.float32)).writeto(star_path) + calls = {} + + def fake_star_pass(exposures, spec, adapter, output_dir): + calls["output_dir"] = output_dir + return star_path, {"resetbits": 0, "driz_cr": False} + + monkeypatch.setattr( + pipeline_mod.combine_mod, "combine_star_pass", fake_star_pass + ) + ctx = self._ctx( + self._spec(psf_star_pass="no_cr"), + single_exposure=False, + work_dir=tmp_path, + ) + science = np.zeros((5, 5)) + star_sci, star_header, prov = _star_pass_image(ctx, science, object()) + assert calls["output_dir"] == tmp_path + assert star_sci is not science + assert star_sci[0, 0] == pytest.approx(7.0) + assert prov["star_source_pass"] == "no_cr_drizzle" + assert prov["star_pass_kwargs"]["resetbits"] == 0 + + +class TestCrDialFailFast: + """reduce_target rejects unsupported dial/instrument combinations + before any download (issues #61/#62).""" + + def _reduce(self, tmp_path, **spec_overrides): + from autoreduce.pipeline import reduce_target + from autoreduce.target import TargetSpec + + spec = TargetSpec(name="x", ra=2.0, dec=-0.1, **spec_overrides) + return reduce_target( + spec, cache_root=tmp_path / "cache", output_root=tmp_path / "out" + ) + + def test_deepcr_needs_the_astrodrizzle_backend(self, tmp_path): + with pytest.raises(ValueError, match="astrodrizzle"): + self._reduce(tmp_path, instrument="nircam_lw", cr_method="deepcr") + + def test_deepcr_needs_a_registered_model(self, tmp_path): + # wfc3_ir is astrodrizzle-combined but has no deepCR model — its + # cosmic rays are already per-frame flagged by calwf3 ramp fitting. + with pytest.raises(ValueError, match="deepCR model"): + self._reduce(tmp_path, instrument="wfc3_ir", cr_method="deepcr") + + def test_no_cr_star_pass_needs_the_astrodrizzle_backend(self, tmp_path): + with pytest.raises(ValueError, match="psf_star_pass"): + self._reduce(tmp_path, instrument="nircam_lw", psf_star_pass="no_cr") + + def test_no_cr_star_pass_conflicts_with_psf_from_frames(self, tmp_path): + # psf_from_frames never finds stars on the mosaic; a silent no-op of + # an explicitly requested second drizzle would hide the mistake. + with pytest.raises(ValueError, match="psf_star_pass"): + self._reduce( + tmp_path, psf_star_pass="no_cr", psf_from_frames=True + ) diff --git a/test_autoreduce/test_target_and_instruments.py b/test_autoreduce/test_target_and_instruments.py index 97d6147..0770177 100644 --- a/test_autoreduce/test_target_and_instruments.py +++ b/test_autoreduce/test_target_and_instruments.py @@ -41,6 +41,21 @@ def test_dec_bounds(self): with pytest.raises(ValueError): TargetSpec(name="x", ra=0.0, dec=91.0) + def test_cr_dials_default_to_current_behaviour(self): + # The default flip to a per-frame route is human-gated on SLACS + # validation (#61); "auto" adds no drizzle pass (#62). + spec = TargetSpec(name="x", ra=0.0, dec=0.0) + assert spec.cr_method == "driz_cr" + assert spec.psf_star_pass == "auto" + + def test_invalid_cr_method_rejected(self): + with pytest.raises(ValueError, match="cr_method"): + TargetSpec(name="x", ra=0.0, dec=0.0, cr_method="lacosmic") + + def test_invalid_psf_star_pass_rejected(self): + with pytest.raises(ValueError, match="psf_star_pass"): + TargetSpec(name="x", ra=0.0, dec=0.0, psf_star_pass="always") + class TestInstrumentRegistry: def test_acs_wfc_registered(self): @@ -102,3 +117,88 @@ def test_drizzle_kwargs_single_vs_multi_exposure(): assert single["final_wht_type"] == "IVM" with pytest.raises(ValueError): drizzle_kwargs_for(spec, adapter, 0) + + +class TestCrMethodDrizzleKwargs: + """The cr_method routes through drizzle_kwargs_for (issue #61).""" + + def test_driz_cr_route_unchanged(self): + # Regression: the default route is byte-identical to the pre-dial + # behaviour — no resetbits key, so AstroDrizzle's own default rules. + from autoreduce.drizzle.combine import drizzle_kwargs_for + + spec = TargetSpec(name="x", ra=0.0, dec=0.0) + kwargs = drizzle_kwargs_for(spec, instruments.get("acs_wfc"), 4) + assert kwargs["driz_cr"] and kwargs["median"] and kwargs["blot"] + assert "resetbits" not in kwargs + + def test_deepcr_route_is_plain_mean_with_resetbits_zero(self): + # The #61 trap: AstroDrizzle's default resetbits=4096 clears exactly + # the DQ bit the per-frame masks were written into. + from autoreduce.drizzle.combine import drizzle_kwargs_for + + spec = TargetSpec(name="x", ra=0.0, dec=0.0, cr_method="deepcr") + kwargs = drizzle_kwargs_for(spec, instruments.get("acs_wfc"), 4) + assert not (kwargs["driz_cr"] or kwargs["median"] or kwargs["blot"]) + assert kwargs["resetbits"] == 0 + + def test_single_exposure_branch_unaffected_by_cr_method(self): + from autoreduce.drizzle.combine import drizzle_kwargs_for + + adapter = instruments.get("acs_wfc") + default = drizzle_kwargs_for(TargetSpec(name="x", ra=0.0, dec=0.0), adapter, 1) + deepcr = drizzle_kwargs_for( + TargetSpec(name="x", ra=0.0, dec=0.0, cr_method="deepcr"), adapter, 1 + ) + assert default == deepcr + assert "resetbits" not in deepcr + + def test_invalid_cr_method_raises_at_spec_construction(self): + with pytest.raises(ValueError, match="cr_method"): + TargetSpec(name="x", ra=0.0, dec=0.0, cr_method="median") + + def test_star_pass_kwargs_ignore_cr_flags_without_clearing_them(self): + # psf_star_pass="no_cr" (#62): plain mean, prior CR DQ flags treated + # as good via final_bits, never cleared (frame products still need + # the science pass's flags in the inputs). + from autoreduce.drizzle.combine import CR_DQ_BIT, star_pass_kwargs_for + + spec = TargetSpec(name="x", ra=0.0, dec=0.0) + kwargs = star_pass_kwargs_for(spec, instruments.get("acs_wfc"), 4) + assert not (kwargs["driz_cr"] or kwargs["median"] or kwargs["blot"]) + assert kwargs["resetbits"] == 0 + assert kwargs["final_bits"] & CR_DQ_BIT + + +class TestDqCrFlagWrite: + """The pure DQ update behind cr_method='deepcr' (issue #61).""" + + def test_sets_clears_and_preserves_dtype(self): + import numpy as np + + from autoreduce.drizzle.combine import CR_DQ_BIT, dq_with_cr_flags + + dq = np.array([[0, CR_DQ_BIT], [16, CR_DQ_BIT | 16]], dtype=np.int16) + mask = np.array([[True, False], [False, True]]) + out = dq_with_cr_flags(dq, mask) + assert out.dtype == np.int16 + # Set where masked; stale CR bits cleared; other bits untouched. + assert out.tolist() == [[CR_DQ_BIT, 0], [16, CR_DQ_BIT | 16]] + + def test_idempotent_on_rerun(self): + import numpy as np + + from autoreduce.drizzle.combine import dq_with_cr_flags + + dq = np.array([[0, 32], [8192, 0]], dtype=np.int16) + mask = np.array([[True, True], [False, False]]) + once = dq_with_cr_flags(dq, mask) + assert dq_with_cr_flags(once, mask).tolist() == once.tolist() + + def test_shape_mismatch_is_loud(self): + import numpy as np + + from autoreduce.drizzle.combine import dq_with_cr_flags + + with pytest.raises(ValueError, match="shape"): + dq_with_cr_flags(np.zeros((2, 2), np.int16), np.zeros((3, 3), bool))