From c514eec6d322d114ffcce2bd636ed7187b38e539 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 7 Aug 2026 20:27:41 +0000 Subject: [PATCH 1/2] Add a local weight-deficit guard for the science region (#65 leg 2) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Both existing coverage guards only answer the *total*-loss question, so a stripe of reduced-but-finite IVM weight through a deflector core ships silently: `mask_isolated_bad_pixels` tests `~isfinite | <= 0`, so a degraded pixel is never a candidate — not for the clustering check, not for the 1.5" `protect_radius_arcsec` protection whose entire purpose is to guarantee the lens itself is clean; and `weight_uniformity` is a global RMS/median (the slacs0008 spike measured 0.066 against a 0.2 limit) that a few columns cannot move. `check_local_weight_deficit` interrogates the same 1.5" science region and reports, as fractions of the cutout median, the region median weight and the worst row/column median. Both axes are tested because a detector-column defect lands on an image row or column depending on the frame's orientation on the sky. Recorded in reduction.json beside weight_uniformity_cutout. The 0.9 limit is provisional and the verdict is RECORDED, NEVER RAISED: losing one exposure of N leaves weight (N-1)/N, so 0.9 detects a single lost exposure for any N <= 9 — the regime where sqrt(N/(N-1)) noise inflation is visible — but it is uncalibrated against real reductions, and a fatal guard at an uncalibrated limit would refuse datasets that are fine. The leg-1 control test calibrates it. Leg 2 lands first deliberately: it is the detector for leg 1, so the control test gets an objective pass/fail instead of an eyeball judgement. The test suite pins the blindness itself — one synthetic striped map that the two pre-existing guards both pass and the new one catches. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_014po4zMWnxACBWxatVBMn2f --- autoreduce/drizzle/diagnostics.py | 117 ++++++++++++++++++++++++ autoreduce/pipeline.py | 9 +- docs/design/hst_acs_pipeline.md | 30 ++++++ test_autoreduce/test_psf_and_package.py | 113 +++++++++++++++++++++++ 4 files changed, 268 insertions(+), 1 deletion(-) diff --git a/autoreduce/drizzle/diagnostics.py b/autoreduce/drizzle/diagnostics.py index dcf3e4d..a9205ce 100644 --- a/autoreduce/drizzle/diagnostics.py +++ b/autoreduce/drizzle/diagnostics.py @@ -30,3 +30,120 @@ def check_weight_uniformity(wht: np.ndarray) -> dict: "limit": WEIGHT_UNIFORMITY_LIMIT, "acceptable": value <= WEIGHT_UNIFORMITY_LIMIT, } + + +# The science region the local diagnostic below interrogates — deliberately +# the same 1.5" as `noise.rms.mask_isolated_bad_pixels`'s +# `protect_radius_arcsec`, because the two guards answer the same question +# ("is the lens itself clean?") about different failure modes. +SCIENCE_RADIUS_ARCSEC = 1.5 + +# A line whose weight falls below this fraction of the cutout median is +# flagged. Derivation: losing one exposure of N along a column leaves IVM +# weight (N-1)/N, so 0.9 detects a single lost exposure for any N <= 9 — +# exactly the few-exposure regime where the sqrt(N/(N-1)) noise inflation +# is visible (x1.41 at N=2 ... x1.06 at N=9). PROVISIONAL: uncalibrated +# against real reductions, which is why the verdict is recorded and never +# raised (issue #65 leg 2). +LOCAL_WEIGHT_DEFICIT_LIMIT = 0.9 + + +def local_weight_deficit( + wht: np.ndarray, + center_xy, + pixel_scale: float, + radius_arcsec: float = SCIENCE_RADIUS_ARCSEC, +) -> dict: + """ + Local coverage deficit inside the science region (issue #65 leg 2). + + `weight_uniformity` is a *global* RMS/median over the whole cutout, so a + handful of degraded columns cannot move it (the slacs0008 spike measured + 0.066 against a 0.2 limit); `mask_isolated_bad_pixels` only sees weight + that has gone non-finite or non-positive. Neither can see the failure + mode this reports: coverage that is finite and positive but materially + *reduced* along a line running through the lens — the ACS/WFC stripe + class, where rejecting hot/warm pixels on column-organised CCD defects + drops one exposure's contribution along a column rather than zeroing it. + + Reports, as fractions of the cutout's median positive weight: + + - ``science_median_ratio`` — median weight inside `radius_arcsec`. Catches + a uniformly-degraded science region. + - ``min_line_ratio`` — the worst row/column median, over lines crossing + the science region. Catches the stripe: a single degraded column barely + shifts the region median but drives this down to (N-1)/N. + + Both axes are tested because a detector-column defect lands on an image + column only for a north-up frame at the detector's own orientation; after + drizzling to a common sky frame it can run either way. + """ + wht = np.asarray(wht, dtype=float) + if wht.ndim != 2: + raise ValueError(f"weight map must be 2D, got shape {wht.shape}") + good = np.isfinite(wht) & (wht > 0.0) + if not good.any(): + raise ValueError("weight map has no positive pixels — empty coverage") + median_all = float(np.median(wht[good])) + + cx, cy = center_xy + ys, xs = np.indices(wht.shape) + r_arcsec = np.hypot(ys - cy, xs - cx) * pixel_scale + science = r_arcsec < radius_arcsec + if not science.any(): + raise ValueError( + f"science region (r < {radius_arcsec}\") falls outside the weight " + f"map — centre {center_xy} is off the cutout" + ) + + # Zero/non-finite weight inside the science region is a *total* loss, not + # a deficit; report it as ratio 0 rather than dropping it from the median, + # so a hole cannot masquerade as clean coverage. + science_weights = np.where(good, wht, 0.0)[science] + science_median_ratio = float(np.median(science_weights) / median_all) + + min_line_ratio = np.inf + min_line_axis, min_line_index = None, None + for axis, name in ((0, "column"), (1, "row")): + # Lines are indexed along the *other* axis: axis=0 collapses rows, so + # each entry is one image column. + crossed = np.flatnonzero(science.any(axis=axis)) + for index in crossed: + line = science.take(index, axis=1 - axis) + values = np.where(good, wht, 0.0).take(index, axis=1 - axis)[line] + ratio = float(np.median(values) / median_all) + if ratio < min_line_ratio: + min_line_ratio, min_line_axis, min_line_index = ratio, name, int(index) + + return { + "science_median_ratio": science_median_ratio, + "min_line_ratio": min_line_ratio, + "min_line_axis": min_line_axis, + "min_line_index": min_line_index, + "radius_arcsec": radius_arcsec, + "n_science_pixels": int(science.sum()), + } + + +def check_local_weight_deficit( + wht: np.ndarray, + center_xy, + pixel_scale: float, + radius_arcsec: float = SCIENCE_RADIUS_ARCSEC, +) -> dict: + """ + Compute the local deficit and its verdict for the provenance record. + + **Records, never raises.** The limit is provisional (see + `LOCAL_WEIGHT_DEFICIT_LIMIT`) and has not been calibrated against real + reductions, so a fatal guard here would refuse datasets that are very + likely fine. `acceptable=False` in `reduction.json` is the signal; the + control test for issue #65 leg 1 is what calibrates the limit. + """ + stats = local_weight_deficit(wht, center_xy, pixel_scale, radius_arcsec) + worst = min(stats["science_median_ratio"], stats["min_line_ratio"]) + return { + **stats, + "limit": LOCAL_WEIGHT_DEFICIT_LIMIT, + "acceptable": worst >= LOCAL_WEIGHT_DEFICIT_LIMIT, + } diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py index 6ef1e86..756a45d 100644 --- a/autoreduce/pipeline.py +++ b/autoreduce/pipeline.py @@ -22,7 +22,7 @@ from .acquire import quality as quality_mod from .align import diagnostics as align_mod from .drizzle import combine as combine_mod -from .drizzle.diagnostics import check_weight_uniformity +from .drizzle.diagnostics import check_local_weight_deficit, check_weight_uniformity from .instruments import InstrumentAdapter from .noise import rms as rms_mod from .package import cosmic_rays as cr_mod @@ -651,6 +651,13 @@ def _package(ctx: _StageContext, sci, header, wht, noise, psf, psf_full) -> None ctx.record["drizzle"]["weight_uniformity_cutout"] = check_weight_uniformity( wht_cut ) + # ...and the cutout-wide statistic still averages the lens away: a stripe + # of reduced-but-finite coverage through the deflector core moves neither + # it nor the bad-pixel policy above (issue #65 leg 2). Recorded, not + # raised — the limit is provisional until the leg-1 control test. + ctx.record["drizzle"]["local_weight_deficit"] = check_local_weight_deficit( + wht_cut, center_xy=center_xy, pixel_scale=spec.final_scale + ) fits.PrimaryHDU(psf.astype(np.float32)).writeto( out_dir / "psf.fits", overwrite=True diff --git a/docs/design/hst_acs_pipeline.md b/docs/design/hst_acs_pipeline.md index c869f19..2fc3086 100644 --- a/docs/design/hst_acs_pipeline.md +++ b/docs/design/hst_acs_pipeline.md @@ -292,6 +292,36 @@ PSFs and `reduction.json`. Optionally emit the auxiliary modeling-prep files (`info.json` skeleton) but leave scientific annotations (positions, extra galaxies) to the modeling workflow — reduction ends at the dataset. +**Coverage guards, and the gap between them (issue #65 leg 2).** Two guards +run over the packaged cutout, and until leg 2 both answered only the +*total*-loss question: + +- `noise.rms.mask_isolated_bad_pixels` masks isolated dead/rejected pixels and + fails loudly on structured clusters, on more than 0.5% of the cutout, or on + any bad pixel within `protect_radius_arcsec` (1.5″) of the target — but it + tests `~isfinite(noise) | noise <= 0`, so a pixel whose coverage is merely + *reduced* is never even a candidate. +- `drizzle.diagnostics.weight_uniformity` is a global RMS/median over the + cutout against a 0.2 limit; the slacs0008 spike measured 0.066, and a + handful of degraded columns cannot shift it. + +Between them sits the ACS/WFC stripe class: finite, positive, but materially +reduced IVM weight along a line through the deflector core. It passes both — +including the 1.5″ protection whose whole purpose is to guarantee the lens +itself is clean. `drizzle.diagnostics.check_local_weight_deficit` closes that +gap: inside the same 1.5″ radius it reports the science-region median weight +and the worst row/column median, each as a fraction of the cutout median, and +records them in `reduction.json` beside `weight_uniformity_cutout`. Both axes +are tested because a detector-column defect maps onto an image row or column +depending on the frame's orientation on the sky. + +The limit (0.9) is **provisional and the verdict is recorded, never raised**: +losing one exposure of N leaves weight (N-1)/N, so 0.9 detects a single lost +exposure for any N ≤ 9 — the regime where the `sqrt(N/(N-1))` noise inflation +is visible — but it has not been calibrated against real reductions, and a +fatal guard at an uncalibrated limit would refuse datasets that are fine. The +leg-1 control test is what calibrates it. + ## Validation — SLACS parity study End-to-end on 2–3 SLACS lenses (e.g. `slacs0008-0004` plus one well-behaved diff --git a/test_autoreduce/test_psf_and_package.py b/test_autoreduce/test_psf_and_package.py index 8079772..2cab0f7 100644 --- a/test_autoreduce/test_psf_and_package.py +++ b/test_autoreduce/test_psf_and_package.py @@ -200,6 +200,119 @@ def test_weight_uniformity_diagnostic(): weight_uniformity(np.zeros((5, 5))) +class TestLocalWeightDeficit: + """ + The local coverage guard (issue #65 leg 2) — the detector for the ACS + stripe class that both pre-existing guards pass silently. + """ + + # 0.05"/pix, so the 1.5" science radius is 30 px around the centre. + PIXEL_SCALE = 0.05 + CENTER = (60.0, 60.0) + + def _striped(self, n_exposures=4): + """ + Uniform coverage with one column through the lens core down to + (N-1)/N — one exposure lost along a CCD column, finite and positive + throughout, exactly the reported failure. + """ + wht = np.full((121, 121), 100.0) + wht[:, 60] = 100.0 * (n_exposures - 1) / n_exposures + return wht + + def test_flat_coverage_is_clean(self): + from autoreduce.drizzle.diagnostics import check_local_weight_deficit + + verdict = check_local_weight_deficit( + np.full((121, 121), 100.0), self.CENTER, self.PIXEL_SCALE + ) + assert verdict["acceptable"] + assert verdict["science_median_ratio"] == pytest.approx(1.0) + assert verdict["min_line_ratio"] == pytest.approx(1.0) + + def test_stripe_through_the_core_is_caught(self): + from autoreduce.drizzle.diagnostics import check_local_weight_deficit + + verdict = check_local_weight_deficit( + self._striped(4), self.CENTER, self.PIXEL_SCALE + ) + assert not verdict["acceptable"] + # One lost exposure of four: the column sits at 3/4 of the median... + assert verdict["min_line_ratio"] == pytest.approx(0.75) + assert verdict["min_line_axis"] == "column" + assert verdict["min_line_index"] == 60 + # ...while the science region as a whole barely notices it, which is + # precisely why a region-median test alone would not do. + assert verdict["science_median_ratio"] == pytest.approx(1.0) + + def test_the_existing_guards_are_blind_to_the_same_map(self): + # The regression this leg exists for: without the diagnostic above, + # a striped reduction ships clean. + from autoreduce.drizzle.diagnostics import check_weight_uniformity + from autoreduce.noise.rms import mask_isolated_bad_pixels + + wht = self._striped(4) + assert check_weight_uniformity(wht)["acceptable"] + + # The bad-pixel policy sees the noise map, so mirror the stripe into + # noise space: reduced IVM weight -> elevated but finite noise. + noise = np.full((121, 121), 1.0) + noise[:, 60] = np.sqrt(4.0 / 3.0) + _, _, diag = mask_isolated_bad_pixels( + np.zeros_like(noise), + noise, + center_xy=self.CENTER, + pixel_scale=self.PIXEL_SCALE, + ) + # No pixel is non-finite or <= 0, so nothing is even a candidate — + # not the clustering check, not the 1.5" lens-core protection. + assert diag["n_masked_pixels"] == 0 + + def test_row_oriented_stripe_is_caught_too(self): + # A detector column lands on an image row after drizzling to a sky + # frame at the right orientation, so both axes are tested. + from autoreduce.drizzle.diagnostics import check_local_weight_deficit + + wht = np.full((121, 121), 100.0) + wht[60, :] = 50.0 + verdict = check_local_weight_deficit(wht, self.CENTER, self.PIXEL_SCALE) + assert not verdict["acceptable"] + assert verdict["min_line_axis"] == "row" + assert verdict["min_line_ratio"] == pytest.approx(0.5) + + def test_zero_coverage_inside_the_region_reads_as_total_loss(self): + from autoreduce.drizzle.diagnostics import check_local_weight_deficit + + wht = np.full((121, 121), 100.0) + wht[:, 60] = 0.0 + verdict = check_local_weight_deficit(wht, self.CENTER, self.PIXEL_SCALE) + assert verdict["min_line_ratio"] == pytest.approx(0.0) + assert not verdict["acceptable"] + + def test_uniformly_degraded_science_region_is_caught(self): + # The complementary failure: no stripe, but the whole lens region + # sits below the cutout median — min_line_ratio and the region + # median agree here. + from autoreduce.drizzle.diagnostics import check_local_weight_deficit + + wht = np.full((121, 121), 100.0) + # The 1.5" radius is 30 px here, so the region spans 30..90. + wht[30:91, 30:91] = 70.0 + verdict = check_local_weight_deficit(wht, self.CENTER, self.PIXEL_SCALE) + assert not verdict["acceptable"] + assert verdict["science_median_ratio"] == pytest.approx(0.7) + + def test_empty_and_off_cutout_inputs_raise(self): + from autoreduce.drizzle.diagnostics import local_weight_deficit + + with pytest.raises(ValueError, match="coverage"): + local_weight_deficit(np.zeros((10, 10)), (5.0, 5.0), self.PIXEL_SCALE) + with pytest.raises(ValueError, match="off the cutout"): + local_weight_deficit( + np.full((10, 10), 1.0), (500.0, 500.0), self.PIXEL_SCALE + ) + + def test_provenance_record(tmp_path): import json From 52342b7bfbd8a293a3a190819278b43cab91300d Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 7 Aug 2026 20:30:49 +0000 Subject: [PATCH 2/2] Key the drizzle DQ bits on exposure count, per MDRIZTAB (#65 leg 1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `drizzle_kwargs_for` never set `final_bits`/`driz_sep_bits` and no adapter did either — the strings appeared nowhere in the repo — so every reduction inherited drizzlepac's package default `final_bits="0"` (drizzlepac/pars/astrodrizzle.cfg:101): no DQ bit treated as good, every flagged pixel rejected, including the hot/warm/blob pixels calacs/calwf3 have already corrected. That is an unintended inheritance, not a chosen deviation, and it produced structured defects on both HST instruments — zero-coverage holes on WFC3/IR (PJ011646, DQ 512 blobs at the same detector pixels in all five exposures) and high-noise stripes through deflector cores on ACS/WFC F814W SLACS-gold. Adapters now carry STScI's own MDRIZTAB rows as (min_exposures, driz_sep_bits, final_bits), with MDRIZTAB's semantics (the last row whose min_exposures <= N): acs_wfc, wfc3_uvis: (1, 65535, 65535), (2, 336, 336) wfc3_ir: (1, 65535, 65535), (2, 65535, 528), (4, 528, 528) 336 = 16+64+256 (hot, warm, saturated); 528 = 512+16 (blob, hot). Two properties of the table are load-bearing. The value is genuinely N-dependent — single-exposure data uses 65535, because with one exposure there is nothing to fill a masked pixel with, so a flat constant would be wrong; this is also an independent explanation for why the legacy SLACS SNAP maps look clean. And the two columns DIFFER for wfc3_ir at N=2-3, which is why rows carry both rather than one bits value per exposure count. TargetSpec.final_bits / driz_sep_bits override the table at every N. reduction.json records each value AND its source (adapter_mdriztab / target_spec / unset) so datasets stay re-derivable as the tables move. Non-AstroDrizzle backends (jwst_image3, nirc2_native) declare no table and emit no keyword at all — absent, not 0, since 0 is the bug. mdriztab=True is deliberately NOT used: it would import final_scale, final_pixfrac, final_kernel and final_rot too, silently reverting the justified lensing deviations in hst_acs_pipeline.md stage 3. star_pass_kwargs_for picks this up for free — its `int(kwargs.get("final_bits", 0)) | CR_DQ_BIT` was the fingerprint of the missing dial and now ORs onto a real value. CONTROL TEST NOT RUN: it needs drizzlepac, the CRDS cache and archive data, none of which exist in the session that wrote this. Re-drizzle a striped SLACS target at 0 vs the MDRIZTAB value and diff the weight/noise maps — leg 2's diagnostic scores it — before these defaults reach a release. If the stripes do not move, revert them rather than shipping the dial anyway. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_014po4zMWnxACBWxatVBMn2f --- autoreduce/drizzle/combine.py | 43 +++++++- autoreduce/instruments/acs_wfc.py | 7 ++ autoreduce/instruments/adapter.py | 37 +++++++ autoreduce/instruments/wfc3_ir.py | 10 ++ autoreduce/instruments/wfc3_uvis.py | 4 + autoreduce/target.py | 20 ++++ docs/design/hst_acs_pipeline.md | 44 ++++++++ docs/design/wfc3.md | 51 +++++++++ .../test_target_and_instruments.py | 103 ++++++++++++++++++ 9 files changed, 318 insertions(+), 1 deletion(-) diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py index cc9ff7a..1a8d7a8 100644 --- a/autoreduce/drizzle/combine.py +++ b/autoreduce/drizzle/combine.py @@ -43,6 +43,18 @@ def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures median=multi, blot=multi, ) + # Which DQ bits count as good (issue #65). Precedence: TargetSpec + # override > adapter's MDRIZTAB-derived table > emit nothing (the + # non-AstroDrizzle backends, which have no such keyword). The keys are + # deliberately absent rather than 0 when unset — 0 is drizzlepac's own + # "no bit is good" default and was the bug. + bits = adapter.dq_bits_for(n_exposures) + if bits is not None: + kwargs.update(driz_sep_bits=bits[0], final_bits=bits[1]) + if spec.driz_sep_bits is not None: + kwargs["driz_sep_bits"] = spec.driz_sep_bits + if spec.final_bits is not None: + kwargs["final_bits"] = spec.final_bits 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 @@ -53,6 +65,32 @@ def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures return kwargs +def dq_bits_provenance( + spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int +) -> Dict: + """ + Where each bits value came from, so a dataset stays re-derivable after + the adapter tables or the user's spec change (issue #65). Pure function. + """ + table = adapter.dq_bits_for(n_exposures) + record = {"n_exposures": n_exposures} + for i, name in enumerate(("driz_sep_bits", "final_bits")): + override = getattr(spec, name) + if override is not None: + record[name] = override + record[f"{name}_source"] = "target_spec" + elif table is not None: + record[name] = table[i] + record[f"{name}_source"] = "adapter_mdriztab" + else: + record[name] = None + record[f"{name}_source"] = "unset" + if adapter.dq_bits_rows is not None: + applicable = [r for r in adapter.dq_bits_rows if r[0] <= n_exposures] + record["adapter_row_min_exposures"] = applicable[-1][0] + return record + + def star_pass_kwargs_for( spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int ) -> Dict: @@ -183,7 +221,10 @@ def _one(suffix: str) -> Path: sci = _one("_sci.fits") wht = _one("_wht.fits") - tail = {"cr_method": spec.cr_method} + tail = { + "cr_method": spec.cr_method, + "dq_bits": dq_bits_provenance(spec, adapter, len(exposures)), + } if per_frame_cr is not None: tail["per_frame_cr"] = per_frame_cr provenance = combine_provenance( diff --git a/autoreduce/instruments/acs_wfc.py b/autoreduce/instruments/acs_wfc.py index c9ea305..9ffa393 100644 --- a/autoreduce/instruments/acs_wfc.py +++ b/autoreduce/instruments/acs_wfc.py @@ -23,5 +23,12 @@ "final_rot": 0.0, }, saturation_dn=80_000.0, + # MDRIZTAB acs/37g1550cj_mdz.fits (issue #65). 336 = 16 + 64 + 256: + # hot, warm and saturated pixels are treated as usable once several + # exposures overlap, because the dark subtraction has already + # corrected them; rejecting them instead removes pixels that are + # partly column-organised on an aged CCD (trap columns, CTE trails), + # reducing IVM weight along those columns and striping the noise map. + dq_bits_rows=((1, 65535, 65535), (2, 336, 336)), ) ) diff --git a/autoreduce/instruments/adapter.py b/autoreduce/instruments/adapter.py index 57073f1..4aa6014 100644 --- a/autoreduce/instruments/adapter.py +++ b/autoreduce/instruments/adapter.py @@ -56,6 +56,43 @@ class InstrumentAdapter: # injected MEAN never depends on it. None where unused (HST e-/s # frames are already electron-referred). e_per_dn: Optional[float] = None + # Which DQ bits AstroDrizzle should treat as GOOD, as a function of + # exposure count (issue #65). Rows are + # ``(min_exposures, driz_sep_bits, final_bits)`` in ascending order and + # carry MDRIZTAB's own semantics: the applicable row is the last one whose + # ``min_exposures <= n``. Mirrors STScI's shipped MDRIZTAB reference files + # rather than enabling ``mdriztab=True``, which would also import + # final_scale/final_pixfrac/final_kernel/final_rot and silently revert the + # deliberate lensing deviations in hst_acs_pipeline.md stage 3. + # + # None = emit no bits keywords at all, leaving the backend's own default + # untouched — correct for the non-AstroDrizzle backends (jwst_image3, + # nirc2_native), which have no such keyword. + dq_bits_rows: Optional[Tuple[Tuple[int, int, int], ...]] = None + + def dq_bits_for(self, n_exposures: int) -> Optional[Tuple[int, int]]: + """ + ``(driz_sep_bits, final_bits)`` for this exposure count, or None when + the instrument declares no table. + + The value is genuinely N-dependent: single-exposure data uses 65535 + (every bit good) because with one exposure there is nothing to fill a + masked pixel with, so the standard recipe keeps flagged pixels rather + than punching holes. A flat constant would be wrong. + """ + if self.dq_bits_rows is None: + return None + if n_exposures < 1: + raise ValueError(f"need at least one exposure, got {n_exposures}") + applicable = [row for row in self.dq_bits_rows if row[0] <= n_exposures] + if not applicable: + raise ValueError( + f"instrument {self.key!r} has no DQ-bits row covering " + f"{n_exposures} exposure(s); rows start at " + f"{self.dq_bits_rows[0][0]}" + ) + _, driz_sep_bits, final_bits = applicable[-1] + return driz_sep_bits, final_bits def ground_detector(self): """The detector constants, loud when a ground stage needs them.""" diff --git a/autoreduce/instruments/wfc3_ir.py b/autoreduce/instruments/wfc3_ir.py index 5ea229a..11c2a1a 100644 --- a/autoreduce/instruments/wfc3_ir.py +++ b/autoreduce/instruments/wfc3_ir.py @@ -30,5 +30,15 @@ }, saturation_dn=78_000.0, recommended_final_scale=0.065, + # MDRIZTAB wfc3/3562021pi_mdz.fits (issue #65). 528 = 512 + 16 — the + # blob bit plus hot pixels. Blobs are detector-fixed IR channel + # features that calwf3 flags but does not remove; rejecting them on + # snapshot data with tiny dithers punches structured zero-coverage + # holes in the mosaic (PJ011646, 5 exposures, a 123-px hole at + # r = 5.3"). The IR rows are the reason this is a table and not a + # pair: the two bits columns DIFFER at N = 2-3, where the separate + # (median-building) drizzle still keeps every bit while the final + # drizzle already drops to 528. + dq_bits_rows=((1, 65535, 65535), (2, 65535, 528), (4, 528, 528)), ) ) diff --git a/autoreduce/instruments/wfc3_uvis.py b/autoreduce/instruments/wfc3_uvis.py index 2596507..6a7781d 100644 --- a/autoreduce/instruments/wfc3_uvis.py +++ b/autoreduce/instruments/wfc3_uvis.py @@ -26,5 +26,9 @@ }, saturation_dn=63_000.0, recommended_final_scale=0.0396, + # MDRIZTAB wfc3/2ck18260i_mdz.fits (issue #65) — identical rows to + # ACS/WFC: 336 = 16 + 64 + 256 (hot, warm, saturated) once two or + # more exposures overlap. + dq_bits_rows=((1, 65535, 65535), (2, 336, 336)), ) ) diff --git a/autoreduce/target.py b/autoreduce/target.py index be3ae95..33d43d3 100644 --- a/autoreduce/target.py +++ b/autoreduce/target.py @@ -42,6 +42,19 @@ class TargetSpec: final_pixfrac: float = 0.8 final_kernel: str = "square" + # DQ bits AstroDrizzle treats as GOOD (issue #65). None = use the + # adapter's exposure-count-keyed table, which mirrors STScI's MDRIZTAB + # reference files; an explicit value overrides it at every N. Set both + # deliberately: `driz_sep_bits` governs the separate drizzle that builds + # the CR-rejection median (inert on the single-exposure branch), while + # `final_bits` governs the shipped mosaic. Leaving them None is the + # normal case — the dial exists because the *previous* behaviour was an + # unintended inheritance of drizzlepac's package default `final_bits=0`, + # which treats no bit as good and rejects every flagged pixel, including + # the hot/warm/blob pixels calacs/calwf3 have already corrected. + final_bits: Optional[int] = None + driz_sep_bits: Optional[int] = None + # 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- @@ -161,6 +174,13 @@ 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}") + for bits_name in ("final_bits", "driz_sep_bits"): + bits = getattr(self, bits_name) + if bits is not None and (not isinstance(bits, int) or bits < 0): + raise ValueError( + f"{bits_name} must be a non-negative int bitmask or None " + f"(None = the adapter's MDRIZTAB-derived default): {bits!r}" + ) if self.cr_method not in ("driz_cr", "deepcr"): raise ValueError( f"cr_method must be 'driz_cr' or 'deepcr': {self.cr_method!r}" diff --git a/docs/design/hst_acs_pipeline.md b/docs/design/hst_acs_pipeline.md index 2fc3086..e9fbb0d 100644 --- a/docs/design/hst_acs_pipeline.md +++ b/docs/design/hst_acs_pipeline.md @@ -114,6 +114,50 @@ 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. +**DQ bits — which flagged pixels count as good (issue #65).** Not a deviation: +this is the pipeline *returning* to STScI practice. No bits keyword was +emitted anywhere until this landed, so every reduction inherited drizzlepac's +package default `final_bits = "0"` (`drizzlepac/pars/astrodrizzle.cfg:101`), +which treats **no** DQ bit as good and rejects every flagged pixel — including +the hot, warm and saturated pixels that `calacs`/`calwf3` have already +corrected. The adapters now mirror STScI's MDRIZTAB reference files, keyed on +exposure count: + +| Detector | Reference file | `numimages` | `driz_sep_bits` | `final_bits` | +|---|---|---|---|---| +| ACS/WFC | `acs/37g1550cj_mdz.fits` | 1 | 65535 | 65535 | +| ACS/WFC | " | ≥2 | 336 | **336** | +| WFC3/UVIS | `wfc3/2ck18260i_mdz.fits` | 1 | 65535 | 65535 | +| WFC3/UVIS | " | ≥2 | 336 | **336** | + +`336 = 16 + 64 + 256` — hot, warm, saturated. `TargetSpec.final_bits` / +`driz_sep_bits` override the table at every N; both the value and its source +(`adapter_mdriztab` / `target_spec` / `unset`) are recorded in +`reduction.json`, so existing datasets stay re-derivable as the tables move. + +*Why it matters for lensing.* Rejecting hot (16) and warm (64) pixels on top +of the genuine bad columns (4, 128) removes pixels that are partly +**column-organised** on an aged ACS CCD — trap columns and CTE trails — so the +IVM weight along those columns is *reduced* rather than zeroed. Noise then +rises by `sqrt(N/(N-1))` per lost exposure: ×1.41 at N=2, ×1.22 at N=3, ×1.15 +at N=4, ×1.08 at N=7. Few-exposure targets stripe visibly and many-exposure +ones do not, which is exactly the pattern reported on the F814W SLACS-gold +noise maps (2026-08-04) and absent from the legacy SLACS reductions. The +single-exposure row is an independent part of that story: SLACS SNAP data is +the N=1 regime, where the standard recipe keeps every flagged pixel. + +**Not `mdriztab=True`.** MDRIZTAB carries the whole parameter set — +`final_scale`, `final_pixfrac`, `final_kernel`, `final_rot` — so enabling it +would silently revert the deviations in the table above. Only the bits columns +are mirrored; our explicit kwargs stand. + +**Control test — not yet run.** Re-drizzle one striped SLACS target at the old +`0` and at the MDRIZTAB value and diff the weight and noise maps; the stage-6 +local weight-deficit diagnostic gives that comparison an objective pass/fail. +If the stripes do not move, the cause is elsewhere — exposure count, dither +geometry, or genuine bad columns — and these defaults must be reverted rather +than kept. Run it before these values reach a release. + **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 diff --git a/docs/design/wfc3.md b/docs/design/wfc3.md index 3a658b7..1ced340 100644 --- a/docs/design/wfc3.md +++ b/docs/design/wfc3.md @@ -35,6 +35,57 @@ scale, R = 1.5; comparisons account for whether R is applied.) | saturation | ~78 ke- effective full well | | psf | same tiers; STScI focus-diverse ePSF grids exist for IR when tier 2 lands | +## DQ bits — the blob bit, and why the IR table has three rows (issue #65) + +Until this landed, no bits keyword was emitted anywhere in the pipeline, so +every reduction inherited drizzlepac's package default `final_bits = "0"` +(`drizzlepac/pars/astrodrizzle.cfg:101`) — **no** DQ bit treated as good, and +every flagged pixel rejected. That is far more aggressive than STScI's own +practice, and on the IR channel it has a specific, structural consequence. + +**Blobs (DQ 512).** IR-channel blobs are detector-fixed features — shadows +cast by particulate contamination on the channel-select mechanism mirror — +which `calwf3` *flags* but does not remove. They are fixed in detector +coordinates, so a dither pattern of a few pixels moves them barely at all: the +same detector pixels are flagged in every exposure of the visit. Reject them +and the mosaic gets a structured **zero-coverage hole** rather than the +speckle that a large dither would produce. PJ011646 (program 14653, F160W, +5 exposures, 2–6 px dithers) failed packaging on exactly this — a single +123-px hole at r = 5.3″, DQ 512 at the same detector pixels in all five +exposures — while a trusted external reduction of the same data has none. +STScI's own MDRIZTAB passes the blob bit, so under standard practice that hole +could not have occurred. + +The adapter now carries STScI's rows, read from the shipped reference file: + +| Reference file | `numimages` | `driz_sep_bits` | `final_bits` | +|---|---|---|---| +| `wfc3/3562021pi_mdz.fits` | 1 | 65535 | 65535 | +| " | 2 | 65535 | **528** | +| " | ≥4 | 528 | 528 | + +`528 = 512 + 16` — blob plus hot. Two things about this table are load-bearing: + +- **The two columns differ at N = 2–3.** The separate (median-building) + drizzle still keeps every bit while the final drizzle is already at 528, so + the adapter stores rows of `(min_exposures, driz_sep_bits, final_bits)` + rather than one bits value per exposure count. UVIS and ACS/WFC keep the two + columns equal at every N, but the IR channel does not, and a single-value + design would silently misreport one of them. +- **Single-exposure data uses 65535 — every bit good.** With one exposure + there is nothing to fill a masked pixel with, so the standard recipe keeps + flagged pixels rather than punching holes. Any fix here has to be N-aware; a + flat constant would be wrong in both directions. + +Rows follow MDRIZTAB's own semantics: the applicable row is the last one whose +`numimages` does not exceed the actual exposure count. + +**Not `mdriztab=True`.** Enabling MDRIZTAB wholesale would import the entire +parameter set — `final_scale`, `final_pixfrac`, `final_kernel`, `final_rot` — +and silently revert the deliberate, justified lensing deviations in +`hst_acs_pipeline.md` stage 3 (0.05″/pix, pixfrac 0.8, north-up). The bits +columns are mirrored onto the adapters instead, and our explicit kwargs stand. + ## Coverage audit vs `ajshajib/hst-lens` (the checklist, not the architecture) Their three notebooks (Download / IR / UVIS) cover: archive download, diff --git a/test_autoreduce/test_target_and_instruments.py b/test_autoreduce/test_target_and_instruments.py index 0770177..dc030b9 100644 --- a/test_autoreduce/test_target_and_instruments.py +++ b/test_autoreduce/test_target_and_instruments.py @@ -119,6 +119,106 @@ def test_drizzle_kwargs_single_vs_multi_exposure(): drizzle_kwargs_for(spec, adapter, 0) +class TestDqBitsDial: + """ + The exposure-count-keyed DQ-bits dial (issue #65 leg 1). Values mirror + STScI's MDRIZTAB reference files; before this, no bits keyword was ever + emitted and every reduction inherited drizzlepac's package default + `final_bits="0"` — no bit treated as good, every flagged pixel rejected. + """ + + def _bits(self, key, n, **spec_kwargs): + from autoreduce.drizzle.combine import drizzle_kwargs_for + + spec = TargetSpec(name="x", ra=0.0, dec=0.0, **spec_kwargs) + kwargs = drizzle_kwargs_for(spec, instruments.get(key), n) + return kwargs["driz_sep_bits"], kwargs["final_bits"] + + @pytest.mark.parametrize("key", ["acs_wfc", "wfc3_uvis"]) + def test_acs_like_rows(self, key): + # Single exposure: every bit good — there is nothing to fill a masked + # pixel with, so the standard recipe keeps flagged pixels. + assert self._bits(key, 1) == (65535, 65535) + # N >= 2: 336 = 16 + 64 + 256 (hot, warm, saturated). + assert self._bits(key, 2) == (336, 336) + assert self._bits(key, 7) == (336, 336) + assert 336 == 16 | 64 | 256 + + def test_wfc3_ir_rows_differ_between_the_two_columns(self): + # The IR table is why this is a table and not a pair: at N = 2-3 the + # separate drizzle keeps every bit while the final drizzle is already + # at 528 = 512 + 16 (blob + hot). + assert self._bits("wfc3_ir", 1) == (65535, 65535) + assert self._bits("wfc3_ir", 2) == (65535, 528) + assert self._bits("wfc3_ir", 3) == (65535, 528) + assert self._bits("wfc3_ir", 4) == (528, 528) + assert 528 == 512 | 16 + + def test_pj011646_would_not_have_holed(self): + # The regression this leg exists for: PJ011646 was F160W with five + # exposures, so it lands on the numimages >= 4 row, where STScI + # passes exactly the blob bit (512) that punched the 123-px hole. + _, final_bits = self._bits("wfc3_ir", 5) + assert final_bits & 512 + + def test_target_spec_overrides_the_adapter_at_every_n(self): + assert self._bits("acs_wfc", 4, final_bits=0) == (336, 0) + assert self._bits("acs_wfc", 1, driz_sep_bits=8) == (8, 65535) + assert self._bits("wfc3_ir", 5, final_bits=65535, driz_sep_bits=65535) == ( + 65535, + 65535, + ) + + def test_non_hst_adapters_emit_no_bits_keywords(self): + # jwst_image3 / nirc2_native have no such keyword; the key must be + # ABSENT rather than 0, since 0 is drizzlepac's "no bit is good". + from autoreduce.drizzle.combine import drizzle_kwargs_for + + for key in ("nircam_sw", "nirc2_narrow"): + adapter = instruments.get(key) + assert adapter.dq_bits_for(4) is None + kwargs = drizzle_kwargs_for( + TargetSpec(name="x", ra=0.0, dec=0.0), adapter, 4 + ) + assert "final_bits" not in kwargs + assert "driz_sep_bits" not in kwargs + + def test_invalid_bits_raise_at_spec_construction(self): + with pytest.raises(ValueError, match="final_bits"): + TargetSpec(name="x", ra=0.0, dec=0.0, final_bits=-1) + with pytest.raises(ValueError, match="driz_sep_bits"): + TargetSpec(name="x", ra=0.0, dec=0.0, driz_sep_bits="336") + + def test_zero_exposures_still_raises(self): + with pytest.raises(ValueError, match="at least one exposure"): + instruments.get("acs_wfc").dq_bits_for(0) + + def test_provenance_records_value_and_source(self): + from autoreduce.drizzle.combine import dq_bits_provenance + + adapter = instruments.get("wfc3_ir") + default = dq_bits_provenance(TargetSpec(name="x", ra=0.0, dec=0.0), adapter, 5) + assert default["final_bits"] == 528 + assert default["final_bits_source"] == "adapter_mdriztab" + assert default["adapter_row_min_exposures"] == 4 + assert default["n_exposures"] == 5 + + overridden = dq_bits_provenance( + TargetSpec(name="x", ra=0.0, dec=0.0, final_bits=65535), adapter, 5 + ) + assert overridden["final_bits"] == 65535 + assert overridden["final_bits_source"] == "target_spec" + # A partial override leaves the other column on the adapter default. + assert overridden["driz_sep_bits_source"] == "adapter_mdriztab" + + unset = dq_bits_provenance( + TargetSpec(name="x", ra=0.0, dec=0.0), instruments.get("nircam_sw"), 4 + ) + assert unset["final_bits"] is None + assert unset["final_bits_source"] == "unset" + assert "adapter_row_min_exposures" not in unset + + class TestCrMethodDrizzleKwargs: """The cr_method routes through drizzle_kwargs_for (issue #61).""" @@ -168,6 +268,9 @@ def test_star_pass_kwargs_ignore_cr_flags_without_clearing_them(self): assert not (kwargs["driz_cr"] or kwargs["median"] or kwargs["blot"]) assert kwargs["resetbits"] == 0 assert kwargs["final_bits"] & CR_DQ_BIT + # It ORs onto whatever the science pass resolved, so it inherits the + # bits dial (#65) rather than the old `.get("final_bits", 0)` zero. + assert kwargs["final_bits"] == 336 | CR_DQ_BIT class TestDqCrFlagWrite: