feat: HST/ACS reduction pipeline (phase 1)

.gitignore

@@ -13,6 +13,8 @@ venv/
 *.fits.gz
 prototypes/output/
 prototypes/cache/
+scripts/output/
+scripts/cache/
 mastDownload/
 
 # Notebooks

autoreduce/__init__.py

@@ -16,3 +16,6 @@
     __version__ = _version("autoreduce")
 except PackageNotFoundError:
     __version__ = "0.0.dev0"
+
+from .target import TargetSpec
+from .pipeline import reduce_target

autoreduce/acquire/cache.py

@@ -0,0 +1,149 @@
+"""
+The transient exposure cache (design doc stage 1).
+
+Full-frame exposures are transient: download per target, reduce, package,
+evict. A manifest records what came from where so eviction never costs
+reproducibility. CRDS reference files live under the same root but are the
+one component never evicted per target — they are shared across targets.
+"""
+
+import json
+import shutil
+import time
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Dict, List, Optional
+
+MANIFEST_NAME = "cache_manifest.json"
+REFERENCES_DIRNAME = "crds"
+
+
+@dataclass
+class ExposureCache:
+    """Size-capped per-target exposure storage with a provenance manifest."""
+
+    root: Path
+    size_cap_bytes: Optional[int] = None  # None = uncapped
+
+    def __post_init__(self):
+        self.root = Path(self.root)
+        self.root.mkdir(parents=True, exist_ok=True)
+
+    # -- manifest -----------------------------------------------------------
+
+    @property
+    def manifest_path(self) -> Path:
+        return self.root / MANIFEST_NAME
+
+    def read_manifest(self) -> Dict:
+        if self.manifest_path.exists():
+            manifest = json.loads(self.manifest_path.read_text())
+            if "targets" not in manifest:
+                raise ValueError(
+                    f"{self.manifest_path} is not an ExposureCache manifest "
+                    f"(keys: {sorted(manifest)}); refusing to guess — point the "
+                    f"cache at a fresh directory (spike-era caches are not "
+                    f"compatible)"
+                )
+            return manifest
+        return {"targets": {}}
+
+    def _write_manifest(self, manifest: Dict) -> None:
+        self.manifest_path.write_text(json.dumps(manifest, indent=2))
+
+    # -- per-target lifecycle ------------------------------------------------
+
+    def target_dir(self, target_name: str) -> Path:
+        return self.root / target_name
+
+    def record_download(
+        self, target_name: str, files: List[str], source: str
+    ) -> None:
+        """Register downloaded exposures so a re-run can re-fetch deterministically."""
+        manifest = self.read_manifest()
+        manifest["targets"][target_name] = {
+            "files": sorted(str(f) for f in files),
+            "source": source,
+            "downloaded_at": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()),
+            "evicted": False,
+        }
+        self._write_manifest(manifest)
+
+    def exposures_for(self, target_name: str) -> List[Path]:
+        entry = self.read_manifest()["targets"].get(target_name)
+        if entry is None or entry["evicted"]:
+            return []
+        paths = [Path(f) for f in entry["files"]]
+        missing = [p for p in paths if not p.exists()]
+        if missing:
+            raise FileNotFoundError(
+                f"cache manifest lists exposures that are gone (not via evict): "
+                f"{[str(m) for m in missing]}"
+            )
+        return paths
+
+    def evict(self, target_name: str) -> None:
+        """Drop a target's exposures; the manifest keeps the provenance."""
+        manifest = self.read_manifest()
+        entry = manifest["targets"].get(target_name)
+        if entry is None:
+            raise KeyError(f"no cache entry for target {target_name!r}")
+        target_dir = self.target_dir(target_name)
+        if target_dir.exists():
+            shutil.rmtree(target_dir)
+        entry["evicted"] = True
+        self._write_manifest(manifest)
+
+    # -- size cap -------------------------------------------------------------
+
+    def size_bytes(self) -> int:
+        """Total evictable payload (excludes the shared CRDS references)."""
+        total = 0
+        for path in self.root.rglob("*"):
+            if (
+                path.is_file()
+                and REFERENCES_DIRNAME not in path.parts
+                and path.name != MANIFEST_NAME
+            ):
+                total += path.stat().st_size
+        return total
+
+    def enforce_cap(self) -> List[str]:
+        """
+        Evict oldest completed targets until under the cap. Returns the
+        evicted target names. Targets are eligible only once marked evictable
+        (their products written) via ``mark_completed``.
+        """
+        if self.size_cap_bytes is None:
+            return []
+        manifest = self.read_manifest()
+        evicted: List[str] = []
+        entries = sorted(
+            (
+                (name, e)
+                for name, e in manifest["targets"].items()
+                if not e["evicted"] and e.get("completed", False)
+            ),
+            key=lambda item: item[1]["downloaded_at"],
+        )
+        for name, _ in entries:
+            if self.size_bytes() <= self.size_cap_bytes:
+                break
+            self.evict(name)
+            evicted.append(name)
+        return evicted
+
+    def mark_completed(self, target_name: str) -> None:
+        """Products for this target are written; its exposures may be evicted."""
+        manifest = self.read_manifest()
+        entry = manifest["targets"].get(target_name)
+        if entry is None:
+            raise KeyError(f"no cache entry for target {target_name!r}")
+        entry["completed"] = True
+        self._write_manifest(manifest)
+
+    # -- CRDS references -------------------------------------------------------
+
+    @property
+    def references_dir(self) -> Path:
+        return self.root / REFERENCES_DIRNAME

autoreduce/acquire/crds.py

@@ -0,0 +1,66 @@
+"""
+CRDS reference-file sync (design doc stage 1, spike finding).
+
+AstroDrizzle's IVM weighting resolves calibration files through the
+adapter's reference environment variable (``jref$`` for ACS), so best
+references must exist locally before the drizzle stage. References are
+shared across targets and are never evicted.
+"""
+
+import os
+import subprocess
+import sys
+from pathlib import Path
+from typing import List
+
+from ..instruments import InstrumentAdapter
+
+CRDS_SERVER_URL = "https://hst-crds.stsci.edu"
+
+
+def configure_environment(references_root: Path, adapter: InstrumentAdapter) -> dict:
+    """
+    Set the CRDS variables for this process. Must run before drizzlepac is
+    imported anywhere in the process. Returns the mapping applied.
+
+    Deliberately overrides any inherited CRDS_PATH/jref: the pipeline is a
+    pure function of the target spec plus the archive, so its reference files
+    live in *its* cache, not wherever the shell environment happens to point.
+    """
+    env = {
+        "CRDS_SERVER_URL": CRDS_SERVER_URL,
+        "CRDS_PATH": str(references_root),
+        adapter.reference_env_key: str(
+            Path(references_root) / adapter.crds_reference_subpath
+        )
+        + "/",
+    }
+    os.environ.update(env)
+    return env
+
+
+def references_present(references_root: Path, adapter: InstrumentAdapter) -> bool:
+    """True if the instrument's reference directory exists and is non-empty."""
+    ref_dir = Path(references_root) / adapter.crds_reference_subpath
+    return ref_dir.is_dir() and any(ref_dir.iterdir())
+
+
+def sync_best_references(exposures: List[Path]) -> None:
+    """Fetch + assign best references for the exposures (network)."""
+    if not exposures:
+        raise ValueError("no exposures to sync references for")
+    cmd = [
+        sys.executable,
+        "-m",
+        "crds.bestrefs",
+        "--files",
+        *[str(p) for p in exposures],
+        "--sync-references=1",
+        "--update-bestrefs",
+    ]
+    result = subprocess.run(cmd, capture_output=True, text=True)
+    if result.returncode != 0:
+        tail = "\n".join(
+            result.stdout.splitlines()[-5:] + result.stderr.splitlines()[-5:]
+        )
+        raise RuntimeError(f"crds.bestrefs failed (exit {result.returncode}):\n{tail}")

autoreduce/acquire/mast.py

@@ -0,0 +1,100 @@
+"""
+MAST acquisition (design doc stage 1).
+
+Query hygiene (spike finding): plain coordinate queries also match HAP
+skycell products, whose member lists re-reference the same exposures many
+times over and pull in neighbouring pointings. We therefore keep only
+*direct* calibration-level observations (numeric proposal IDs, obs_id not a
+``hst_skycell`` product) and optionally filter by proposal, then download the
+adapter's calibrated exposure products.
+"""
+
+from pathlib import Path
+from typing import List, Optional, Sequence
+
+from ..instruments import InstrumentAdapter
+
+
+def is_direct_observation(obs_id: str, proposal_id: str) -> bool:
+    """True for a direct program observation, False for HAP skycell products."""
+    if str(obs_id).startswith("hst_skycell"):
+        return False
+    proposal = str(proposal_id).strip()
+    return proposal.isdigit()
+
+
+def select_observations(
+    obs_table,
+    proposal_ids: Optional[Sequence[str]] = None,
+):
+    """Filter a MAST observation table to direct program observations."""
+    keep = []
+    for row in obs_table:
+        if not is_direct_observation(row["obs_id"], row["proposal_id"]):
+            continue
+        if proposal_ids is not None and str(row["proposal_id"]) not in set(
+            str(p) for p in proposal_ids
+        ):
+            continue
+        keep.append(row)
+    return keep
+
+
+def query_exposures(
+    ra: float,
+    dec: float,
+    adapter: InstrumentAdapter,
+    filter_name: str,
+    radius: str = "0.5 arcmin",
+    proposal_ids: Optional[Sequence[str]] = None,
+):
+    """Query MAST for direct observations of the target. Network."""
+    from astropy.coordinates import SkyCoord
+    from astroquery.mast import Observations
+
+    coord = SkyCoord(ra, dec, unit="deg")
+    obs = Observations.query_criteria(
+        coordinates=coord,
+        radius=radius,
+        obs_collection="HST",
+        instrument_name=adapter.mast_instrument_name,
+        filters=filter_name,
+        dataproduct_type="image",
+    )
+    selected = select_observations(obs, proposal_ids=proposal_ids)
+    if not selected:
+        raise LookupError(
+            f"no direct {adapter.mast_instrument_name} {filter_name} observations "
+            f"at ({ra}, {dec}) within {radius}"
+            + (f" for proposals {list(proposal_ids)}" if proposal_ids else "")
+        )
+    return selected
+
+
+def download_exposures(
+    observations,
+    adapter: InstrumentAdapter,
+    download_dir: Path,
+) -> List[Path]:
+    """Download the calibrated exposure products for the observations. Network."""
+    from astropy.table import vstack
+    from astroquery.mast import Observations
+
+    products = vstack([Observations.get_product_list(row) for row in observations])
+    calibrated = Observations.filter_products(
+        products,
+        productSubGroupDescription=[adapter.calibrated_suffix],
+        mrp_only=False,
+    )
+    if len(calibrated) == 0:
+        raise LookupError(
+            f"observations carry no {adapter.calibrated_suffix} products"
+        )
+    Observations.download_products(calibrated, download_dir=str(download_dir))
+    suffix = f"_{adapter.calibrated_suffix.lower()}.fits"
+    paths = sorted(set(Path(download_dir).rglob(f"*{suffix}")))
+    if not paths:
+        raise FileNotFoundError(
+            f"download reported success but no *{suffix} files under {download_dir}"
+        )
+    return list(paths)

autoreduce/align/diagnostics.py

@@ -0,0 +1,36 @@
+"""
+Alignment (design doc stage 2): trust the MAST a-priori WCS by default;
+TweakReg refinement is a *triggered* fallback, not a default step.
+
+The trigger diagnostic compares each exposure's WCS-predicted position of
+the brightest compact source near the target against the stack consensus;
+sub-tolerance scatter means the a-priori solutions are good enough for
+drizzling and TweakReg is skipped.
+"""
+
+from pathlib import Path
+from typing import Dict, List
+
+
+def wcs_solution_names(exposures: List[Path]) -> Dict[str, str]:
+    """Record which WCS solution each exposure carries (provenance)."""
+    from astropy.io import fits
+
+    names = {}
+    for path in exposures:
+        with fits.open(path) as hdul:
+            header = hdul["SCI", 1].header
+            names[Path(path).name] = header.get("WCSNAME", "unknown")
+    return names
+
+
+def run_tweakreg(exposures: List[Path]) -> None:
+    """Relative alignment refinement. Only called when the trigger demands."""
+    from drizzlepac import tweakreg
+
+    tweakreg.TweakReg(
+        [str(p) for p in exposures],
+        interactive=False,
+        updatehdr=True,
+        shiftfile=False,
+    )