Fix critical/high findings from whole-codebase audit

CHANGELOG.md

@@ -6,6 +6,38 @@ match the `VERSION` file and `v*` git tags.
 
 ## [Unreleased]
 
+### Changed
+
+- **`--auto` no longer applies one preset's settings to every target.** A setting only some target
+  presets define was blended over those presets alone, so e.g. the galaxy-only coarse chroma pass,
+  the nebula anisotropic diffusion and `pre_gradient_removal` ran on star fields and clusters too.
+  Presets that don't list a setting now vote for its default. **Expect different `--auto` output**
+  on most targets, and re-save any `_config.toml` you rely on.
+
+### Fixed
+
+- Darks whose exposure differs from the lights are now scaled on the default parallel path (pool
+  workers never received `dark_exptime` and subtracted the dark unscaled, silently).
+- Hierarchical runs: sessions with darks no longer lose their `info.json` metadata (Bayer pattern,
+  WCS, GPS, target name) to a variable-shadowing bug.
+- `--photometry` and `--color-calibrate` work again: the stacked cube's 3-axis WCS made Gaia
+  projection fail silently, and the colour-calibration scale was always exactly 1.
+- With `--auto` on 15+ frames, `--stack-method median`/`linear_fit`/`ivw`/`wavelet` are honoured
+  instead of silently becoming an unrejected patch-weighted mean.
+- Seeded affine star matching failed for any shift over ~5 px (sign error), leaving frames on
+  translation-only registration whenever the blind matcher also failed.
+- Resuming from a phase-1 or phase-2 checkpoint now runs target inference and the `--auto` advisor.
+- Blind PSF estimation (used by `--deconvolve` under `--auto`) no longer flattens the PSF; the
+  broken "blind RL refinement" was removed.
+- `--transient-detect`: S_corr has unit variance when the two epochs' flux scales differ (the ZOGY
+  kernel flux powers were swapped), so thresholds mean what they say.
+- Origin `DATE-OBS` is local time: the `TIMEZONE` keyword is now applied (and kept in the stacked
+  header), fixing the `--fix-atmospheric-dispersion` zenith angle and light-curve MJD/airmass.
+  A below-horizon zenith angle now declines instead of correcting the wrong way.
+- `--stream`: a pixel seeded by a single frame no longer freezes on that frame's value.
+- `--use-gpu`: the next target in the same process no longer reuses the previous target's GPU-side
+  calibration masters.
+
 ## [2.3.0] - 2026-09-23
 
 ### Added

ext/astro_native/Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "astro_native"
-version = "0.36.0"
+version = "0.37.0"
 edition = "2021"
 description = "Native (Rust) hot-path kernels for OriginStack: stacking combine, etc."
 

ext/astro_native/pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "astro_native"
-version = "0.36.0"
+version = "0.37.0"
 description = "Native Rust hot-path kernels for OriginStack"
 requires-python = ">=3.10"
 classifiers = ["Programming Language :: Rust"]

ext/astro_native/src/lib.rs

@@ -569,25 +569,34 @@
 /// borders uncovered). Running the accept-test against that fabricated state
 /// would reject the first real sample forever; instead initialize directly
 /// from it.
+///
+/// Below `ONLINE_CLIP_MIN_SAMPLES` accepted samples there is no spread to
+/// clip against (one sample has M2 = 0, so every later sample failed the
+/// 1e-6 test and the pixel froze on a single frame), so samples are accepted
+/// unconditionally until the state holds that many.
 #[inline]
 fn fold_pixel(mean: f64, m2: f64, n_acc: f64, x: f64, sigma: f64) -> (f64, f64, f64, bool) {
     if n_acc <= 0.0 {
         return (x, 0.0, 1.0, true);
     }
-    let var_est = m2 / n_acc;
-    let std_est = var_est.max(1e-12).sqrt();
-    if (x - mean).abs() <= sigma * std_est {
-        let n_acc_new = n_acc + 1.0;
-        let delta = x - mean;
-        let new_mean = mean + delta / n_acc_new;
-        let delta2 = x - new_mean;
-        let new_m2 = m2 + delta * delta2;
-        (new_mean, new_m2, n_acc_new, true)
-    } else {
-        (mean, m2, n_acc, false)
+    if n_acc >= ONLINE_CLIP_MIN_SAMPLES {
+        let std_est = (m2 / n_acc).max(1e-12).sqrt();
+        if (x - mean).abs() > sigma * std_est {
+            return (mean, m2, n_acc, false);
+        }
     }
+    let n_acc_new = n_acc + 1.0;
+    let delta = x - mean;
+    let new_mean = mean + delta / n_acc_new;
+    let delta2 = x - new_mean;
+    let new_m2 = m2 + delta * delta2;
+    (new_mean, new_m2, n_acc_new, true)
 }
 
+/// Accepted samples a pixel needs before the online sigma-clip starts
+/// rejecting. Mirrored by `ONLINE_CLIP_MIN_SAMPLES` in src/stacking.py.
+const ONLINE_CLIP_MIN_SAMPLES: f64 = 3.0;
+
 /// Per-pixel online sigma-clip: a MAD-rejected burn-in window (first `k =
 /// min(burn_in, n)` samples) seeds a running (mean, M2) Welford state; each
 /// remaining sample is tested against that running estimate before being
@@ -1382,7 +1391,7 @@
 }
 
 #[pyfunction]
 fn dwt2_native<'py>(
     py: Python<'py>,
     img: PyReadonlyArray2<'py, f64>,
     dec_lo: PyReadonlyArray1<'py, f64>,
@@ -1462,7 +1471,7 @@
 }
 
 #[pyfunction]
 fn idwt2_native<'py>(
     py: Python<'py>,
     ca: PyReadonlyArray2<'py, f64>,
     ch: PyReadonlyArray2<'py, f64>,
@@ -3190,7 +3199,7 @@
 /// reference wherever the values are f32-representable.
 #[pyfunction]
 #[pyo3(signature = (channel, emission_mask, patch_size, masked_frac_thresh, sky_ref, sky_std))]
 fn dbe_sample_patches<'py>(
     py: Python<'py>,
     channel: PyReadonlyArray2<'py, f32>,
     emission_mask: PyReadonlyArray2<'py, f32>,

src/auto_settings.py

@@ -503,13 +503,12 @@ def _apply_dynamic_settings(
 ) -> List[str]:
     """Continuous replacement for _apply_target_settings(): every parameter
     is blended across all 8 presets weighted by _blend_weights, instead of
-    looking up one bucket's fixed table. Numeric parameters get a weighted
-    average (renormalized over just the presets that define that
-    parameter); boolean/string parameters -- which can't be fractionally
-    blended -- take the value from whichever contributing preset has the
-    single highest weight (nearest-neighbor-in-signal-space), so the
-    decision still shifts continuously with the signals even though the
-    final choice at any instant is binary.
+    looking up one bucket's fixed table. A preset that does not list a
+    parameter contributes the baseline value args already holds. Numeric
+    parameters get a weighted average; boolean/string parameters -- which
+    can't be fractionally blended -- take the value with the largest total
+    weight, so the decision still shifts continuously with the signals even
+    though the final choice at any instant is binary.
     """
     changes: List[str] = []
     _explicit = getattr(args, '_explicit_cli_dests', set())
@@ -532,18 +531,33 @@ def _set(attr: str, val: object) -> None:
                 seen.add(attr)
                 all_attrs.append(attr)
 
+    # A preset that does not list an attr means "leave it at the baseline",
+    # not "abstain": it votes for the value args held before this blend, with
+    # its own weight. Renormalizing over only the listing presets applied a
+    # one-preset value at full strength to every target -- galaxy's coarse
+    # chroma pass (chroma_nr_large_sigma 0 -> 50) on a 97%-star-field blend.
+    # A baseline of None (attr absent from args) has nothing to vote for, so
+    # that attr falls back to the listing presets alone.
     for attr in all_attrs:
-        contributors = [(t, val) for t, rows in _TARGET_SETTINGS.items()
-                        for a, val in rows if a == attr]
-        w_sum = sum(weights.get(t, 0.0) for t, _ in contributors)
+        listed = {t: val for t, rows in _TARGET_SETTINGS.items()
+                  for a, val in rows if a == attr}
+        sample_val = next(iter(listed.values()))
+        baseline = getattr(args, attr, None)
+        votes = [(weights.get(t, 0.0), listed[t] if t in listed else baseline)
+                 for t in _TARGET_SETTINGS
+                 if t in listed or baseline is not None]
+        w_sum = sum(w for w, _ in votes)
         if w_sum <= 0:
             continue
-        sample_val = contributors[0][1]
         if isinstance(sample_val, (bool, str)):
-            best_t, best_val = max(contributors, key=lambda tv: weights.get(tv[0], 0.0))
-            _set(attr, best_val)
+            # Largest total weight per value, not the single heaviest preset:
+            # six presets that leave a flag off outweigh one that turns it on.
+            totals: Dict[object, float] = {}
+            for w, val in votes:
+                totals[val] = totals.get(val, 0.0) + w
+            _set(attr, max(totals, key=totals.get))
         else:
-            blended = sum(weights.get(t, 0.0) * val for t, val in contributors) / w_sum
+            blended = sum(w * float(val) for w, val in votes) / w_sum
             if isinstance(sample_val, int):
                 blended = int(round(blended))
             else:

src/cli.py

@@ -6,7 +6,6 @@
 import logging
 import math
 import os
-import re
 import sys
 import tempfile
 import threading
@@ -27,6 +26,7 @@
 from src.models import Config, ProcessingStats
 from src.pipeline import stack_target
 from src.utils import (
+    TZ_OFFSET_RE,
     disable_astropy_network,
     format_time,
     print_header,
@@ -676,10 +676,6 @@
 _ROTATION_SPLIT_THRESHOLD_DEG = 3.0
 
 
-# A FITS TIMEZONE value this code can append to DATE-OBS: '-0700', '+05:30'.
-# Anything else ('PDT', 'US/Pacific', '-07:00 (PDT)') would produce an
-# unparseable timestamp and silently disable the rotation prediction.
-_TZ_OFFSET_RE = re.compile(r'[+-]\d{2}:?\d{2}')
 
 
 def _predict_rotation_spread(subdirs: List[str]) -> Optional[float]:
@@ -728,7 +724,7 @@
         h0, h1 = lights[0].header or {}, lights[-1].header or {}
 
         offset = str(h0.get('TIMEZONE', '') or '').strip()
-        if offset and not _TZ_OFFSET_RE.fullmatch(offset):
+        if offset and not TZ_OFFSET_RE.fullmatch(offset):
             return _give_up(d, f"TIMEZONE {offset!r} is not a +/-HHMM offset")
 
         ra, dec = math.degrees(si.ra_rad), math.degrees(si.dec_rad)
@@ -1050,9 +1046,11 @@
                     safe_print(f"    Bias:  pedestal={b_med:.1f} ADU  "
                                f"noise={b_std:.1f} ADU  -> {b_quality}")
                 if masters.get('dark') is not None:
-                    d = masters['dark']
-                    dark_med = float(np.median(d))
-                    dark_peak = float(d.max())
+                    # not `d`: that is the enclosing loop's target directory,
+                    # read again below for args._input_directory
+                    dk = masters['dark']
+                    dark_med = float(np.median(dk))
+                    dark_peak = float(dk.max())
                     dark_et = masters.get('dark_exptime')
                     dark_hdr = frames['dark'][0].header if frames.get('dark') else {}
                     dark_temp_c = dark_hdr.get('CCD-TEMP')
@@ -2642,7 +2640,7 @@
                       "— rerun the same command to resume.")
         raise SystemExit(130)  # 128 + SIGINT, standard convention
     except Exception as e:
         print(f'ERROR: {str(e)}', file=sys.stderr)
         import traceback
         traceback.print_exc()
         raise SystemExit(1)

src/color_calibrate.py

@@ -47,7 +47,7 @@
 def _field_radius_deg(header) -> float:
     """Estimate the field-of-view radius in degrees from WCS."""
     try:
-        wcs = WCS(header)
+        wcs = WCS(header).celestial  # (3, H, W) cube header -> 2-axis
         naxis1 = int(header.get("NAXIS1", 0))
         naxis2 = int(header.get("NAXIS2", 0))
         if naxis1 == 0 or naxis2 == 0:
@@ -236,6 +236,21 @@
     return fluxes[0], fluxes[1], fluxes[2]
 
 
+def _channel_correction(meas: np.ndarray, expected: np.ndarray) -> float:
+    """Multiplicative factor that brings a channel's measured star fluxes to
+    the catalogue's expected ones: median(expected / measured).
+
+    ``apply_photometric_calibration`` multiplies the image by this, so a
+    channel reading too bright gets a factor below 1. Units (ADU vs. the
+    catalogue's 10^(-0.4 m)) are common to all three channels and cancel in
+    the callers' mean normalisation. This used to divide each star's ratio by
+    the channel's own median ratio before taking the median -- always 1.0 --
+    so --color-calibrate never changed the image.
+    """
+    ratio = np.maximum(expected, 1e-30) / np.maximum(meas, 1e-30)
+    return float(np.median(ratio))
+
+
 def fit_channel_scales_spcc(img: np.ndarray, header, catalog,
                             channel_response=None,
                             verbose: bool = False) -> Tuple[float, float, float]:
@@ -262,7 +277,7 @@
     valid = np.all(np.isfinite(fluxes) & (fluxes > 0), axis=1)
     if valid.sum() < 10:
         if verbose:
             print(f"  [SPCC] Too few valid stars ({valid.sum()}) — skipping")
         return 1.0, 1.0, 1.0
     fluxes = fluxes[valid]
 
@@ -313,14 +328,9 @@
     flux_r_expected[fallback_mask] = 10.0 ** (-0.4 * fallback_r[fallback_mask])
     flux_b_expected[fallback_mask] = 10.0 ** (-0.4 * fallback_b[fallback_mask])
 
-    def _robust_ratio(meas: np.ndarray, expected: np.ndarray) -> float:
-        ratio = meas / np.maximum(expected, 1e-30)
-        ratio = ratio / np.median(ratio)
-        return float(np.median(ratio))
-
-    scale_r = _robust_ratio(fluxes[:, 0], flux_r_expected)
-    scale_g = _robust_ratio(fluxes[:, 1], flux_g_expected)
-    scale_b = _robust_ratio(fluxes[:, 2], flux_b_expected)
+    scale_r = _channel_correction(fluxes[:, 0], flux_r_expected)
+    scale_g = _channel_correction(fluxes[:, 1], flux_g_expected)
+    scale_b = _channel_correction(fluxes[:, 2], flux_b_expected)
 
     mean_scale = (scale_r + scale_g + scale_b) / 3.0
     if mean_scale > 0:
@@ -403,15 +413,10 @@
         flux_b_expected = 10.0 ** (-0.4 * k)
 
     # Compute per-star measured ratios vs expected ratios
-    # scale_R = median(flux_R_measured / flux_R_expected) normalised to G
-    def _robust_ratio(meas: np.ndarray, expected: np.ndarray) -> float:
-        ratio = meas / np.maximum(expected, 1e-30)
-        ratio = ratio / np.median(ratio)  # normalise so green ≈ 1
-        return float(np.median(ratio))
-
-    scale_r = _robust_ratio(fluxes[:, 0], flux_r_expected)
-    scale_g = _robust_ratio(fluxes[:, 1], flux_g_expected)
-    scale_b = _robust_ratio(fluxes[:, 2], flux_b_expected)
+    # Per-channel correction = median(expected / measured); mean-normalised below.
+    scale_r = _channel_correction(fluxes[:, 0], flux_r_expected)
+    scale_g = _channel_correction(fluxes[:, 1], flux_g_expected)
+    scale_b = _channel_correction(fluxes[:, 2], flux_b_expected)
 
     # Normalise so that mean(scale) = 1 (preserve overall brightness)
     mean_scale = (scale_r + scale_g + scale_b) / 3.0

src/difference_imaging.py

@@ -276,9 +276,13 @@ def zogy(new: np.ndarray, ref: np.ndarray,
     del d_hat, pd_hat
 
     # --- Noise terms for S_corr (eq. 26-30) ---
-    # Matched-filter kernels for each input image.
-    kn_hat = fr * fn ** 2 * np.conj(pn_hat) * abs_pr2 / denom
-    kr_hat = fn * fr ** 2 * np.conj(pr_hat) * abs_pn2 / denom
+    # Matched-filter kernels for each input image: S = kn*N - kr*R, which is
+    # F_D * D (x) P_D expanded (eq. 28-29). The flux powers were once swapped,
+    # invisible at flux_new == flux_ref == 1 but mis-scaling the variance
+    # whenever the epochs differ (S_corr std 0.4-0.65 on pure noise at a
+    # flux ratio of 3 or 0.3, so a "5 sigma" cut really sat at 8-12 sigma).
+    kn_hat = fn * fr ** 2 * np.conj(pn_hat) * abs_pr2 / denom
+    kr_hat = fr * fn ** 2 * np.conj(pr_hat) * abs_pn2 / denom
     del abs_pn2, abs_pr2, denom, pn_hat, pr_hat
 
     # Variance maps live on the padded grid too. The padding is sky, so it