fix: key the drizzle DQ bits on exposure count, and add a guard that can see the defect (#65)

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(

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,
+    }

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)),
     )
 )

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."""

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)),
     )
 )

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)),
     )
 )

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

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}"