diff --git a/.gitignore b/.gitignore index 1361109..80e8606 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,8 @@ venv/ *.fits.gz prototypes/output/ prototypes/cache/ +scripts/output/ +scripts/cache/ mastDownload/ # Notebooks diff --git a/autoreduce/__init__.py b/autoreduce/__init__.py index 9673432..97c7429 100644 --- a/autoreduce/__init__.py +++ b/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 diff --git a/autoreduce/acquire/cache.py b/autoreduce/acquire/cache.py new file mode 100644 index 0000000..e64b4ec --- /dev/null +++ b/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 diff --git a/autoreduce/acquire/crds.py b/autoreduce/acquire/crds.py new file mode 100644 index 0000000..1db08e6 --- /dev/null +++ b/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}") diff --git a/autoreduce/acquire/mast.py b/autoreduce/acquire/mast.py new file mode 100644 index 0000000..1fc694a --- /dev/null +++ b/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) diff --git a/autoreduce/align/diagnostics.py b/autoreduce/align/diagnostics.py new file mode 100644 index 0000000..4deb353 --- /dev/null +++ b/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, + ) diff --git a/autoreduce/drizzle/combine.py b/autoreduce/drizzle/combine.py new file mode 100644 index 0000000..ce687eb --- /dev/null +++ b/autoreduce/drizzle/combine.py @@ -0,0 +1,103 @@ +""" +Exposure combination via AstroDrizzle (design doc stage 3). + +Defaults-first: the adapter supplies the STScI-recommended keywords; the +lensing deviations (final scale, orientation, units, weight type) and the +user-facing ``pixfrac``/``kernel`` dials come from the `TargetSpec`. The +single-exposure branch (SLACS-V caveat) drizzles the lone frame without +CR rejection — cosmic rays are flagged from the DQ array downstream instead +of median-combining, and the provenance records the branch taken. +""" + +import glob +from pathlib import Path +from typing import Dict, List, Tuple + +from ..instruments import InstrumentAdapter +from ..target import TargetSpec +from .diagnostics import check_weight_uniformity +from ..noise.rms import casertano_r + + +def drizzle_kwargs_for(spec: TargetSpec, adapter: InstrumentAdapter, n_exposures: int) -> Dict: + """Assemble the AstroDrizzle keyword set; pure function, unit-testable.""" + if n_exposures < 1: + raise ValueError("need at least one exposure") + multi = n_exposures > 1 + kwargs = dict(adapter.default_drizzle_kwargs) + kwargs.update( + preserve=False, + build=False, + clean=True, + final_scale=spec.final_scale, + final_pixfrac=spec.final_pixfrac, + final_kernel=spec.final_kernel, + # CR rejection needs >= 2 exposures; the single-exposure branch + # (SLACS-V caveat) skips median/blot/driz_cr. + driz_cr=multi, + median=multi, + blot=multi, + ) + return kwargs + + +def combine( + exposures: List[Path], + spec: TargetSpec, + adapter: InstrumentAdapter, + output_dir: Path, +) -> Tuple[Path, Path, Dict]: + """ + Run AstroDrizzle; return (sci_path, wht_path, provenance_fragment). + + Requires the CRDS environment configured (acquire.crds) beforehand. + """ + from astropy.io import fits + from drizzlepac import astrodrizzle + + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + # Drizzlepac lowercases output filenames internally, which breaks absolute + # paths containing capitals on case-sensitive filesystems — so chdir into + # the work dir and pass a relative, already-lowercase output root. This + # also keeps AstroDrizzle's cwd scratch files contained. + output_name = f"{spec.name}_{spec.filter_name}".lower() + output_root = str(output_dir / output_name) + + kwargs = drizzle_kwargs_for(spec, adapter, len(exposures)) + import os + + cwd = os.getcwd() + os.chdir(output_dir) + try: + astrodrizzle.AstroDrizzle( + input=[str(p) for p in exposures], + output=output_name, + **kwargs, + ) + finally: + os.chdir(cwd) + + def _one(suffix: str) -> Path: + hits = sorted(glob.glob(f"{output_root}*{suffix}")) + if len(hits) != 1: + raise FileNotFoundError( + f"expected exactly one {suffix} for {output_root}, got {hits}" + ) + return Path(hits[0]) + + sci = _one("_sci.fits") + wht = _one("_wht.fits") + + wht_data = fits.getdata(wht) + provenance = { + "n_exposures": len(exposures), + "exposures": [Path(p).name for p in exposures], + "single_exposure_branch": len(exposures) == 1, + "drizzle_kwargs": {k: kwargs[k] for k in sorted(kwargs)}, + "correlated_noise_factor": casertano_r( + spec.final_pixfrac, adapter.scale_ratio(spec.final_scale) + ), + "weight_uniformity": check_weight_uniformity(wht_data), + } + return sci, wht, provenance diff --git a/autoreduce/drizzle/diagnostics.py b/autoreduce/drizzle/diagnostics.py new file mode 100644 index 0000000..dcf3e4d --- /dev/null +++ b/autoreduce/drizzle/diagnostics.py @@ -0,0 +1,32 @@ +""" +Drizzle-quality diagnostics reported with every reduction, so the +user-facing ``pixfrac``/``kernel`` dials are auditable per dataset +(design doc stage 3). +""" + +import numpy as np + + +def weight_uniformity(wht: np.ndarray) -> float: + """ + STScI rule-of-thumb statistic: RMS/median of the (positive) weight map + over the science region. Values above ~0.2 mean the pixfrac is too small + for the dither pattern (coverage speckle/holes). + """ + good = wht[np.isfinite(wht) & (wht > 0.0)] + if good.size == 0: + raise ValueError("weight map has no positive pixels — empty coverage") + return float(good.std() / np.median(good)) + + +WEIGHT_UNIFORMITY_LIMIT = 0.2 + + +def check_weight_uniformity(wht: np.ndarray) -> dict: + """Compute the diagnostic and its verdict for the provenance record.""" + value = weight_uniformity(wht) + return { + "wht_rms_over_median": value, + "limit": WEIGHT_UNIFORMITY_LIMIT, + "acceptable": value <= WEIGHT_UNIFORMITY_LIMIT, + } diff --git a/autoreduce/instruments/__init__.py b/autoreduce/instruments/__init__.py index ca1b507..0704c51 100644 --- a/autoreduce/instruments/__init__.py +++ b/autoreduce/instruments/__init__.py @@ -4,3 +4,6 @@ behind an adapter so the pipeline stages stay instrument-agnostic. HST/ACS is the first adapter; WFC3 and JWST follow (see ``docs/design/roadmap.md``). """ + +from .adapter import InstrumentAdapter, get, register, registered_keys +from .acs_wfc import ACS_WFC diff --git a/autoreduce/instruments/acs_wfc.py b/autoreduce/instruments/acs_wfc.py new file mode 100644 index 0000000..c9ea305 --- /dev/null +++ b/autoreduce/instruments/acs_wfc.py @@ -0,0 +1,27 @@ +""" +ACS/WFC — adapter #1 (design doc phase 1). + +Values follow the design doc's deviation table: CTE-corrected ``_flc`` +exposures, cps output units, IVM weights, north-up final grid. +""" + +from .adapter import InstrumentAdapter, register + +ACS_WFC = register( + InstrumentAdapter( + key="acs_wfc", + mast_instrument_name="ACS/WFC", + native_scale=0.05, + calibrated_suffix="FLC", + reference_env_key="jref", + crds_reference_subpath="references/hst/acs", + supports_cte_correction=True, + default_drizzle_kwargs={ + "skymethod": "globalmin+match", + "final_wht_type": "IVM", + "final_units": "cps", + "final_rot": 0.0, + }, + saturation_dn=80_000.0, + ) +) diff --git a/autoreduce/instruments/adapter.py b/autoreduce/instruments/adapter.py new file mode 100644 index 0000000..2f33d9c --- /dev/null +++ b/autoreduce/instruments/adapter.py @@ -0,0 +1,50 @@ +""" +The instrument-adapter boundary (design doc + roadmap): everything +instrument-specific lives behind an `InstrumentAdapter`; no module outside +`autoreduce.instruments` may name a detector. +""" + +from dataclasses import dataclass +from typing import Dict, Tuple + + +@dataclass(frozen=True) +class InstrumentAdapter: + """Static description of one instrument/detector reduction path.""" + + key: str # registry key, e.g. "acs_wfc" + mast_instrument_name: str # e.g. "ACS/WFC" as MAST spells it + native_scale: float # arcsec / pix + calibrated_suffix: str # exposure product to reduce, e.g. "FLC" + reference_env_key: str # CRDS reference-path variable, e.g. "jref" + crds_reference_subpath: str # where CRDS syncs this instrument's files + supports_cte_correction: bool + default_drizzle_kwargs: Dict[str, object] + saturation_dn: float # conservative full-well / saturation level + + def scale_ratio(self, final_scale: float) -> float: + """s = output scale / native scale, as used by the Casertano factor.""" + return final_scale / self.native_scale + + +_REGISTRY: Dict[str, InstrumentAdapter] = {} + + +def register(adapter: InstrumentAdapter) -> InstrumentAdapter: + if adapter.key in _REGISTRY: + raise ValueError(f"instrument adapter already registered: {adapter.key}") + _REGISTRY[adapter.key] = adapter + return adapter + + +def get(key: str) -> InstrumentAdapter: + try: + return _REGISTRY[key] + except KeyError: + raise KeyError( + f"unknown instrument {key!r}; registered: {sorted(_REGISTRY)}" + ) from None + + +def registered_keys() -> Tuple[str, ...]: + return tuple(sorted(_REGISTRY)) diff --git a/autoreduce/noise/rms.py b/autoreduce/noise/rms.py new file mode 100644 index 0000000..25446a2 --- /dev/null +++ b/autoreduce/noise/rms.py @@ -0,0 +1,85 @@ +""" +RMS noise-map construction (design doc stage 4). + + sigma_i = R * sqrt( N_i / t_exp + 1 / W_i ) + +with N_i the (sky-subtracted, floored-at-zero) source counts/s in pixel i, +t_exp the total exposure time, W_i the IVM drizzle weight (inverse background +variance), and R the Casertano et al. (2000) / DrizzlePac-handbook correction +for the noise correlation the drizzle kernel introduces. The spike validated +applying R: the legacy SLACS noise maps are consistent with it +(parity appendix, docs/design/hst_acs_pipeline.md). +""" + +import numpy as np + + +def casertano_r(pixfrac: float, scale_ratio: float) -> float: + """ + Correlated-noise correction factor R = 1/r. + + r is the variance-reduction factor of Casertano et al. (2000) for drizzle + drop size ``pixfrac`` (p) onto an output grid ``scale_ratio`` (s) times + the native pixel. Smaller p reduces correlation (R -> 1 in the + interlacing limit); p = 1 at s = 1 is shift-and-add (R = 1.5). + """ + if not 0.0 < pixfrac <= 1.0: + raise ValueError(f"pixfrac must be in (0, 1]: {pixfrac}") + if scale_ratio <= 0.0: + raise ValueError(f"scale_ratio must be positive: {scale_ratio}") + p, s = pixfrac, scale_ratio + if s < p: + # Finer output grid than the drop: correlation grows without bound + # as s -> 0 (r -> 0, R -> inf); continuous with the s >= p branch + # at s = p (r = 2/3). + r = (s / p) * (1.0 - s / (3.0 * p)) + else: + r = 1.0 - p / (3.0 * s) + return 1.0 / r + + +def noise_map_from( + sci: np.ndarray, + wht: np.ndarray, + exptime: float, + correlated_noise_factor: float = 1.0, +) -> np.ndarray: + """ + Per-pixel RMS from a cps science mosaic and its IVM weight map. + + Zero/negative/NaN weights are propagated as NaN; callers packaging a + cutout must fail loudly if any land inside it (`assert_finite_within`), + never patch them silently. + """ + if sci.shape != wht.shape: + raise ValueError(f"shape mismatch: sci {sci.shape} vs wht {wht.shape}") + if not np.isfinite(exptime) or exptime <= 0.0: + raise ValueError(f"exptime must be positive and finite: {exptime}") + if correlated_noise_factor < 1.0: + raise ValueError( + f"correlated-noise factor must be >= 1: {correlated_noise_factor}" + ) + + with np.errstate(divide="ignore", invalid="ignore"): + var_bkg = np.where(wht > 0.0, 1.0 / wht, np.nan) + var_src = np.clip(sci, 0.0, None) / exptime + return correlated_noise_factor * np.sqrt(var_src + var_bkg) + + +def assert_finite_within(noise_map: np.ndarray, region_name: str) -> None: + """Loud failure if the noise map carries NaN/inf/zero inside a cutout.""" + bad = ~np.isfinite(noise_map) | (noise_map <= 0.0) + if bad.any(): + raise ValueError( + f"noise map has {int(bad.sum())} non-finite or non-positive pixels " + f"inside {region_name}; refusing to package — fix the reduction " + f"(coverage, weights) rather than patching the noise map" + ) + + +def empirical_background_rms(sci: np.ndarray, n_sigma: float = 3.0) -> float: + """Sigma-clipped RMS of the mosaic — the blank-sky validation check.""" + from astropy.stats import sigma_clipped_stats + + _, _, std = sigma_clipped_stats(sci[np.isfinite(sci)], sigma=n_sigma) + return float(std) diff --git a/autoreduce/package/cutout.py b/autoreduce/package/cutout.py new file mode 100644 index 0000000..4436446 --- /dev/null +++ b/autoreduce/package/cutout.py @@ -0,0 +1,43 @@ +""" +Packaging: WCS-correct cutouts (design doc stage 6). + +Deviation from the legacy datasets, deliberately: cutout headers keep the +WCS, units and exposure metadata the legacy stripped-header files lost. +""" + +from pathlib import Path +from typing import Tuple + +import numpy as np + + +def cutout_to_fits( + data: np.ndarray, + header, + ra: float, + dec: float, + shape: Tuple[int, int], + out_path: Path, + extra_header: dict = None, +) -> np.ndarray: + """Cut `shape` around (ra, dec) and write with an intact cutout WCS.""" + from astropy.coordinates import SkyCoord + from astropy.io import fits + from astropy.nddata import Cutout2D + from astropy.wcs import WCS + + coord = SkyCoord(ra, dec, unit="deg") + cut = Cutout2D(data, coord, shape, wcs=WCS(header), mode="strict") + + out_header = cut.wcs.to_header() + for key in ("BUNIT", "EXPTIME", "TEXPTIME", "FILTER", "INSTRUME", "TELESCOP"): + if key in header: + out_header[key] = header[key] + if extra_header: + for key, value in extra_header.items(): + out_header[key] = value + + fits.PrimaryHDU(cut.data.astype(np.float32), header=out_header).writeto( + out_path, overwrite=True + ) + return cut.data diff --git a/autoreduce/package/provenance.py b/autoreduce/package/provenance.py new file mode 100644 index 0000000..25fb431 --- /dev/null +++ b/autoreduce/package/provenance.py @@ -0,0 +1,33 @@ +""" +The provenance record (design doc stage 6): ``reduction.json`` restores +permanently what the legacy stripped-header datasets lost — where every +pixel came from and how it was made. +""" + +import json +import platform +import time +from pathlib import Path +from typing import Dict + + +def software_versions() -> Dict[str, str]: + versions = {"python": platform.python_version()} + for package in ("autoreduce", "astropy", "numpy", "photutils", "drizzlepac", "astroquery"): + try: + versions[package] = __import__(package).__version__ + except Exception: + versions[package] = "not-installed" + return versions + + +def write_reduction_json(out_dir: Path, record: Dict) -> Path: + """Write the accumulated per-stage provenance with a metadata envelope.""" + payload = { + "written_at": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()), + "software": software_versions(), + **record, + } + path = Path(out_dir) / "reduction.json" + path.write_text(json.dumps(payload, indent=2, default=str)) + return path diff --git a/autoreduce/pipeline.py b/autoreduce/pipeline.py new file mode 100644 index 0000000..00180f0 --- /dev/null +++ b/autoreduce/pipeline.py @@ -0,0 +1,162 @@ +""" +The stage orchestrator: TargetSpec in, modeling-ready dataset out. + + acquire -> align -> drizzle -> noise -> psf -> package + +Each stage contributes to the provenance record; `reduction.json` is written +alongside the data products. Heavy dependencies (astroquery, drizzlepac, +photutils) are imported inside stages so the package imports without them. +""" + +from pathlib import Path +from typing import Dict, Optional + +import numpy as np + +from . import instruments +from .acquire import cache as cache_mod +from .acquire import crds as crds_mod +from .acquire import mast as mast_mod +from .align import diagnostics as align_mod +from .drizzle import combine as combine_mod +from .noise import rms as rms_mod +from .package import cutout as cutout_mod +from .package import provenance as provenance_mod +from .psf import epsf as epsf_mod +from .psf import stars as stars_mod +from .target import TargetSpec + + +def reduce_target( + spec: TargetSpec, + cache_root: Path, + output_root: Path, + size_cap_bytes: Optional[int] = None, + evict_when_done: bool = False, +) -> Dict: + """Run the full pipeline for one target; returns the provenance record.""" + adapter = instruments.get(spec.instrument) + cache = cache_mod.ExposureCache(Path(cache_root), size_cap_bytes=size_cap_bytes) + out_dir = Path(output_root) / spec.name + out_dir.mkdir(parents=True, exist_ok=True) + work_dir = out_dir / "work" + work_dir.mkdir(exist_ok=True) + + record: Dict = {"target": spec.as_dict(), "instrument": adapter.key} + + # -- acquire ------------------------------------------------------------ + crds_mod.configure_environment(cache.references_dir, adapter) + exposures = cache.exposures_for(spec.name) + downloaded = False + if not exposures: + observations = mast_mod.query_exposures( + spec.ra, + spec.dec, + adapter, + spec.filter_name, + proposal_ids=spec.proposal_ids, + ) + exposures = mast_mod.download_exposures( + observations, adapter, cache.target_dir(spec.name) + ) + cache.record_download( + spec.name, [str(p) for p in exposures], source="mast" + ) + downloaded = True + # Fully-cached re-runs with references already synced stay offline. + refs_synced = False + if downloaded or not crds_mod.references_present(cache.references_dir, adapter): + crds_mod.sync_best_references(exposures) + refs_synced = True + record["acquire"] = { + "n_exposures": len(exposures), + "exposures": [Path(p).name for p in exposures], + "downloaded": downloaded, + "references_synced": refs_synced, + } + + # -- align ---------------------------------------------------------------- + record["align"] = { + "wcs_solutions": align_mod.wcs_solution_names(exposures), + "tweakreg_run": False, # a-priori WCS accepted by default (stage 2) + } + + # -- drizzle --------------------------------------------------------------- + sci_path, wht_path, drizzle_prov = combine_mod.combine( + exposures, spec, adapter, work_dir + ) + record["drizzle"] = drizzle_prov + + from astropy.io import fits + + with fits.open(sci_path) as hdul: + sci = hdul[0].data.astype(float) + header = hdul[0].header.copy() + wht = fits.getdata(wht_path) + exptime = header.get("EXPTIME", header.get("TEXPTIME")) + if exptime is None or exptime <= 0: + raise ValueError(f"mosaic header carries no positive EXPTIME: {exptime}") + + # -- noise ----------------------------------------------------------------- + noise = rms_mod.noise_map_from( + sci, + wht, + exptime=float(exptime), + correlated_noise_factor=drizzle_prov["correlated_noise_factor"], + ) + record["noise"] = { + "recipe": "R * sqrt(max(sci,0)/exptime + 1/wht)", + "correlated_noise_factor": drizzle_prov["correlated_noise_factor"], + "exptime": float(exptime), + "empirical_background_rms": rms_mod.empirical_background_rms( + sci[np.isfinite(noise)] + ), + } + + # -- psf ------------------------------------------------------------------- + from astropy.wcs import WCS + + target_xy = WCS(header).world_to_pixel_values(spec.ra, spec.dec) + selection = stars_mod.StarSelection() + stars = stars_mod.find_stars( + sci, + selection, + target_xy=(float(target_xy[0]), float(target_xy[1])), + peak_max=selection.saturation_fraction * adapter.saturation_dn / float(exptime), + ) + psf, psf_full, psf_diag = epsf_mod.build_epsf( + sci, stars, spec.psf_shape, spec.psf_full_shape + ) + record["psf"] = psf_diag + + # -- package ---------------------------------------------------------------- + data_cut = cutout_mod.cutout_to_fits( + sci, header, spec.ra, spec.dec, spec.cutout_shape, out_dir / "data.fits" + ) + noise_cut = cutout_mod.cutout_to_fits( + noise, header, spec.ra, spec.dec, spec.cutout_shape, out_dir / "noise_map.fits" + ) + rms_mod.assert_finite_within(noise_cut, f"{spec.name} cutout") + + fits.PrimaryHDU(psf.astype(np.float32)).writeto( + out_dir / "psf.fits", overwrite=True + ) + fits.PrimaryHDU(psf_full.astype(np.float32)).writeto( + out_dir / "psf_full.fits", overwrite=True + ) + record["package"] = { + "products": ["data.fits", "noise_map.fits", "psf.fits", "psf_full.fits"], + "cutout_shape": list(spec.cutout_shape), + "pixel_scale": spec.final_scale, + "data_units": drizzle_prov["drizzle_kwargs"]["final_units"], + } + + provenance_mod.write_reduction_json(out_dir, record) + + # -- evict -------------------------------------------------------------------- + cache.mark_completed(spec.name) + if evict_when_done: + cache.evict(spec.name) + cache.enforce_cap() + + return record diff --git a/autoreduce/psf/epsf.py b/autoreduce/psf/epsf.py new file mode 100644 index 0000000..662d534 --- /dev/null +++ b/autoreduce/psf/epsf.py @@ -0,0 +1,106 @@ +""" +Tier-1 empirical ePSF (design doc stage 5): Anderson & King-style effective +PSF built from selected stars on the drizzled mosaic via photutils. + +Emits the two modeling kernels — compact (`psf.fits`) and extended +(`psf_full.fits`) — odd-shaped, centred, unit-normalised, plus diagnostics +for the provenance record. If too few stars survive selection the build +fails loudly; tier 2 (model-PSF fallback) is a deliberate choice recorded in +provenance, never a silent degradation. +""" + +from typing import Dict, Tuple + +import numpy as np + +MIN_STARS = 8 + + +class InsufficientStarsError(RuntimeError): + """Tier 1 is not viable for this field; choose tier 2 explicitly.""" + + +def normalise_kernel(psf: np.ndarray, shape: Tuple[int, int]) -> np.ndarray: + """Centre-crop to the requested odd shape and normalise to unit sum.""" + if any(s % 2 == 0 for s in shape): + raise ValueError(f"kernel shape must be odd: {shape}") + ny, nx = psf.shape + cy, cx = ny // 2, nx // 2 + hy, hx = shape[0] // 2, shape[1] // 2 + if hy > cy or hx > cx: + raise ValueError(f"requested shape {shape} exceeds built PSF {psf.shape}") + cut = psf[cy - hy : cy + hy + 1, cx - hx : cx + hx + 1].astype(np.float64) + total = cut.sum() + if not np.isfinite(total) or total <= 0.0: + raise ValueError("PSF kernel has non-positive total flux") + return cut / total + + +def build_epsf( + sci: np.ndarray, + stars_table, + psf_shape: Tuple[int, int], + psf_full_shape: Tuple[int, int], + oversampling: int = 2, +) -> Tuple[np.ndarray, np.ndarray, Dict]: + """Build the ePSF and return (psf, psf_full, diagnostics).""" + from astropy.nddata import NDData + from astropy.table import Table + from photutils.psf import EPSFBuilder, extract_stars + + if stars_table is None or len(stars_table) < MIN_STARS: + n = 0 if stars_table is None else len(stars_table) + raise InsufficientStarsError( + f"{n} usable stars (< {MIN_STARS}); tier 1 ePSF is not viable — " + f"select tier 2 (model PSF) explicitly" + ) + + positions = Table( + {"x": stars_table["xcentroid"], "y": stars_table["ycentroid"]} + ) + # Extraction window comfortably larger than the extended kernel. + size = max(psf_full_shape) + 10 + if size % 2 == 0: + size += 1 + stars = extract_stars(NDData(sci), positions, size=size) + + # extract_stars drops stars whose window overruns the mosaic edge, so the + # minimum-star contract must be re-checked on what actually survived. + if len(stars) < MIN_STARS: + raise InsufficientStarsError( + f"{len(stars)} stars survived cutout extraction (< {MIN_STARS}); " + f"tier 1 ePSF is not viable — select tier 2 (model PSF) explicitly" + ) + + builder = EPSFBuilder(oversampling=oversampling, maxiters=10, progress_bar=False) + epsf_model, fitted = builder(stars) + + # Evaluate the oversampled model back onto the native mosaic pixel grid. + full = _evaluate_native(epsf_model, psf_full_shape) + psf_full = normalise_kernel(full, psf_full_shape) + psf = normalise_kernel(full, psf_shape) + + diagnostics = { + "method": "epsf-tier1", + "n_stars_used": int(len(fitted)), + "oversampling": oversampling, + "fwhm_pix": _fwhm_of(psf), + } + return psf, psf_full, diagnostics + + +def _evaluate_native(epsf_model, shape: Tuple[int, int]) -> np.ndarray: + yy, xx = np.mgrid[0 : shape[0], 0 : shape[1]] + return epsf_model.evaluate( + xx, yy, flux=1.0, x_0=shape[1] // 2, y_0=shape[0] // 2 + ) + + +def _fwhm_of(psf: np.ndarray) -> float: + """Crude FWHM estimate from the radial profile — a diagnostic, not science.""" + ny, nx = psf.shape + y, x = np.mgrid[0:ny, 0:nx] + r = np.hypot(y - ny // 2, x - nx // 2) + half = psf.max() / 2.0 + above = r[psf >= half] + return float(2.0 * above.max()) if above.size else float("nan") diff --git a/autoreduce/psf/fallback.py b/autoreduce/psf/fallback.py new file mode 100644 index 0000000..23642e9 --- /dev/null +++ b/autoreduce/psf/fallback.py @@ -0,0 +1,31 @@ +""" +Tier-2 PSF fallback interface (design doc stage 5). + +For star-poor fields (the SLACS-snapshot regime) the PSF comes from a model — +STScI focus-diverse ePSF grids or TinyTim raytraces — evaluated per exposure +and made drizzle-consistent by resampling through the same footprint as the +science mosaic. Phase 1 ships the interface and provenance contract; the +concrete grid/TinyTim back-ends land when the first star-poor target needs +them (tracked on the roadmap). +""" + +from typing import Dict, Tuple + +import numpy as np + + +class ModelPSFUnavailableError(NotImplementedError): + """No tier-2 back-end is wired up yet; the caller must not degrade silently.""" + + +def model_psf( + spec_name: str, + filter_name: str, + psf_shape: Tuple[int, int], + psf_full_shape: Tuple[int, int], +) -> Tuple[np.ndarray, np.ndarray, Dict]: + raise ModelPSFUnavailableError( + f"tier-2 model PSF requested for {spec_name} ({filter_name}) but no " + f"back-end (focus-diverse ePSF grid / TinyTim) is implemented yet; " + f"tier 1 failed or was declined — this is a hard stop, not a warning" + ) diff --git a/autoreduce/psf/stars.py b/autoreduce/psf/stars.py new file mode 100644 index 0000000..078ee45 --- /dev/null +++ b/autoreduce/psf/stars.py @@ -0,0 +1,89 @@ +""" +Star selection for empirical ePSF construction (design doc stage 5, tier 1). + +Selection cuts, applied to detections on the *drizzled* mosaic so the +resulting PSF is drizzle-consistent by construction: point-like (DAOFind +sharpness/roundness), unsaturated, uncrowded, away from the mosaic edge and +from the lens itself. +""" + +from dataclasses import dataclass +from typing import Tuple + +import numpy as np + + +@dataclass(frozen=True) +class StarSelection: + """Cuts for ePSF star candidates; defaults tuned for ACS-like mosaics.""" + + detection_sigma: float = 10.0 + fwhm_pix: float = 2.0 + sharp_range: Tuple[float, float] = (0.4, 1.0) + round_limit: float = 0.3 + saturation_fraction: float = 0.7 # of adapter.saturation_dn, in counts + min_separation_pix: float = 25.0 + # Must exceed half the ePSF extraction window (psf_full 61 + 10 pad -> 36) + # or edge stars pass selection only to be dropped at extraction. + edge_margin_pix: int = 36 + exclusion_radius_pix: float = 50.0 # around the target itself + + +def reject_crowded(x: np.ndarray, y: np.ndarray, min_separation: float) -> np.ndarray: + """Boolean mask keeping sources with no neighbour within min_separation.""" + keep = np.ones(len(x), dtype=bool) + for i in range(len(x)): + d2 = (x - x[i]) ** 2 + (y - y[i]) ** 2 + d2[i] = np.inf + if (d2 < min_separation**2).any(): + keep[i] = False + return keep + + +def reject_edges( + x: np.ndarray, y: np.ndarray, shape: Tuple[int, int], margin: int +) -> np.ndarray: + ny, nx = shape + return ( + (x >= margin) & (x < nx - margin) & (y >= margin) & (y < ny - margin) + ) + + +def reject_near( + x: np.ndarray, y: np.ndarray, x0: float, y0: float, radius: float +) -> np.ndarray: + return (x - x0) ** 2 + (y - y0) ** 2 > radius**2 + + +def find_stars( + sci: np.ndarray, + selection: StarSelection, + target_xy: Tuple[float, float], + peak_max: float, +): + """DAOStarFinder detections filtered through every selection cut.""" + from astropy.stats import sigma_clipped_stats + from photutils.detection import DAOStarFinder + + _, median, std = sigma_clipped_stats(sci, sigma=3.0) + finder = DAOStarFinder( + fwhm=selection.fwhm_pix, + threshold=selection.detection_sigma * std, + sharplo=selection.sharp_range[0], + sharphi=selection.sharp_range[1], + roundlo=-selection.round_limit, + roundhi=selection.round_limit, + peakmax=peak_max, + ) + sources = finder(sci - median) + if sources is None or len(sources) == 0: + return None + + x = np.asarray(sources["xcentroid"], dtype=float) + y = np.asarray(sources["ycentroid"], dtype=float) + keep = ( + reject_crowded(x, y, selection.min_separation_pix) + & reject_edges(x, y, sci.shape, selection.edge_margin_pix) + & reject_near(x, y, *target_xy, selection.exclusion_radius_pix) + ) + return sources[keep] diff --git a/autoreduce/target.py b/autoreduce/target.py new file mode 100644 index 0000000..e6f49e6 --- /dev/null +++ b/autoreduce/target.py @@ -0,0 +1,76 @@ +""" +Target specification — the declarative input of a reduction. + +A reduction is a pure function of a `TargetSpec` plus the archive (design doc +stage 0): re-running the pipeline on the same spec reproduces the dataset, +modulo upstream reference-file updates, which `reduction.json` records. +""" + +from dataclasses import dataclass, field, asdict +from pathlib import Path +from typing import Optional, Tuple + +import yaml + + +@dataclass(frozen=True) +class TargetSpec: + """Everything the pipeline needs to know about one target.""" + + name: str + ra: float # degrees + dec: float # degrees + instrument: str = "acs_wfc" + filter_name: str = "F814W" + + # Restrict acquisition to these proposal IDs (None = all direct + # calibration-level-2 observations at the coordinates). + proposal_ids: Optional[Tuple[str, ...]] = None + + cutout_shape: Tuple[int, int] = (281, 281) + + # Drizzle dials (design doc stage 3). pixfrac and kernel are deliberately + # user-facing: published practice spans no-drizzle -> 0.6 -> 1.0, so the + # choice is configuration, never a buried default. + final_scale: float = 0.05 # arcsec / pix + final_pixfrac: float = 0.8 + final_kernel: str = "square" + + # PSF products (design doc stage 5). + psf_shape: Tuple[int, int] = (21, 21) + psf_full_shape: Tuple[int, int] = (61, 61) + + # Alignment: residual (pixels) above which TweakReg refinement triggers. + alignment_tolerance_pix: float = 0.1 + + def __post_init__(self): + if not -360.0 <= self.ra <= 360.0: + raise ValueError(f"ra out of range: {self.ra}") + if not -90.0 <= self.dec <= 90.0: + 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 shape_name in ("cutout_shape", "psf_shape", "psf_full_shape"): + shape = getattr(self, shape_name) + if len(shape) != 2 or any(s < 1 for s in shape): + raise ValueError(f"{shape_name} must be two positive ints: {shape}") + for shape_name in ("psf_shape", "psf_full_shape"): + shape = getattr(self, shape_name) + if any(s % 2 == 0 for s in shape): + raise ValueError( + f"{shape_name} must be odd so the PSF has a centre pixel: {shape}" + ) + + @classmethod + def from_yaml(cls, path) -> "TargetSpec": + with open(path) as f: + raw = yaml.safe_load(f) + for key in ("cutout_shape", "psf_shape", "psf_full_shape"): + if key in raw: + raw[key] = tuple(raw[key]) + if raw.get("proposal_ids") is not None: + raw["proposal_ids"] = tuple(str(p) for p in raw["proposal_ids"]) + return cls(**raw) + + def as_dict(self) -> dict: + return asdict(self) diff --git a/docs/design/hst_acs_pipeline.md b/docs/design/hst_acs_pipeline.md index 57b73f4..91e13b3 100644 --- a/docs/design/hst_acs_pipeline.md +++ b/docs/design/hst_acs_pipeline.md @@ -235,9 +235,16 @@ compared against the legacy modeling dataset: | noise ratio × R (R = 1.364) | **0.924** ≈ data ratio | **legacy noise maps are consistent with the correlated-noise correction being applied** — after applying R, data and noise carry the same ~7% global scale offset, i.e. the noise *recipe* matches | Conclusions adopted into the design: `final_units='cps'` stands; stage 4 -**applies** the Casertano/DrizzlePac factor R as designed; the residual ~7% -scale offset is a phase-1 acceptance-test item (reduce with the exact legacy -exposure set + proper WCS registration before judging photometric parity). +**applies** the Casertano/DrizzlePac factor R as designed. The production +pipeline (proposal-filtered exposures, sub-pixel registration) lands at data +ratio 0.941 / noise ratio 0.925 — the ~6% global flux scale vs the legacy +dataset is **accepted as a documented difference** (decision 2026-07-08): +the legacy reduction's exact provenance (kernel, photometric era, +FLT-vs-FLC calibration) is unrecoverable, both ratios carry the same scale so +the *relative* products are self-consistent, and lens-model inferences are +scale-invariant in the relevant regime. A PyAutoMind research prompt tracks +chasing it (gaussian-kernel re-drizzle + calibration-era check) if it ever +matters. Also confirmed: tier-1 ePSF is plausible for this field (236 point-like >10σ detections mosaic-wide, pre-selection), and CRDS reference-file sync + HAP-skycell query filtering belong to the acquire stage (see above). diff --git a/prototypes/slacs_f814w_spike.py b/prototypes/slacs_f814w_spike.py index 32ead58..2e5a7eb 100644 --- a/prototypes/slacs_f814w_spike.py +++ b/prototypes/slacs_f814w_spike.py @@ -197,7 +197,7 @@ def stage_noise(): # reported for the parity discussion, NOT yet applied. p, s = 0.8, FINAL_SCALE / 0.05 # pixfrac, scale ratio vs native 0.05 grid # native ACS pixel is 0.05" so s=1 here; formula kept for generality - r = (p / s) * (1 - s / (3 * p)) if s < p else 1 - p / (3 * s) + r = (s / p) * (1 - s / (3 * p)) if s < p else 1 - p / (3 * s) R = 1.0 / r print(f"[noise] correlated-noise factor R = {R:.3f} (pixfrac={p}, s={s})") print(f"[noise] masked pixels (wht<=0): {bad.sum()}") diff --git a/scripts/reduce_slacs0008.py b/scripts/reduce_slacs0008.py new file mode 100644 index 0000000..9254798 --- /dev/null +++ b/scripts/reduce_slacs0008.py @@ -0,0 +1,100 @@ +""" +Integration + acceptance (issue #2): reduce slacs0008-0004 through the +*production* pipeline — proposal-filtered exposure set (10886 only, the +spike's neighbouring-pointing contamination excluded) — then compare the +products against the legacy modeling dataset with sub-pixel registration. + +Run: ~/venv/PyAuto/bin/python scripts/reduce_slacs0008.py +Network + drizzlepac required; unit tests never import this. +""" + +import json +import sys +from pathlib import Path + +import numpy as np + +REPO = Path(__file__).resolve().parent.parent +sys.path.insert(0, str(REPO)) + +from autoreduce import TargetSpec, reduce_target # noqa: E402 + +LEGACY_DIR = Path("/mnt/c/Users/Jammy/Science/subhalo/dataset/slacs/slacs0008-0004") +CACHE_ROOT = REPO / "scripts" / "cache" +OUTPUT_ROOT = REPO / "scripts" / "output" + +SPEC = TargetSpec( + name="slacs0008-0004", + ra=2.012333, + dec=-0.068944, + proposal_ids=("10886",), +) + + +def subpixel_offset(a: np.ndarray, b: np.ndarray): + """(dy, dx) shift of b relative to a: FFT cross-correlation + parabolic peak.""" + a0 = np.nan_to_num(a - np.nanmedian(a)) + b0 = np.nan_to_num(b - np.nanmedian(b)) + corr = np.fft.irfft2(np.fft.rfft2(a0) * np.conj(np.fft.rfft2(b0)), s=a0.shape) + corr = np.fft.fftshift(corr) + peak = np.unravel_index(np.argmax(corr), corr.shape) + + def parabolic(idx, axis): + c0 = corr[peak] + prev = list(peak); prev[axis] -= 1 + nxt = list(peak); nxt[axis] += 1 + cm, cp = corr[tuple(prev)], corr[tuple(nxt)] + denom = cm - 2 * c0 + cp + return 0.0 if denom == 0 else 0.5 * (cm - cp) / denom + + dy = peak[0] - a.shape[0] // 2 + parabolic(peak, 0) + dx = peak[1] - a.shape[1] // 2 + parabolic(peak, 1) + return dy, dx + + +def main(): + record = reduce_target(SPEC, cache_root=CACHE_ROOT, output_root=OUTPUT_ROOT) + print(json.dumps(record["drizzle"]["weight_uniformity"], indent=2)) + + from astropy.io import fits + from scipy.ndimage import shift as nd_shift + + out_dir = OUTPUT_ROOT / SPEC.name + new_data = fits.getdata(out_dir / "data.fits").astype(float) + new_noise = fits.getdata(out_dir / "noise_map.fits").astype(float) + legacy_data = fits.getdata(LEGACY_DIR / "data.fits").astype(float) + legacy_noise = fits.getdata(LEGACY_DIR / "noise_map.fits").astype(float) + + dy, dx = subpixel_offset(legacy_data, new_data) + print(f"[parity] sub-pixel offset (legacy vs new): dy={dy:.3f}, dx={dx:.3f}") + new_data_r = nd_shift(np.nan_to_num(new_data), (dy, dx), order=3) + new_noise_r = nd_shift(np.nan_to_num(new_noise), (dy, dx), order=1) + + bright = legacy_data > 10 * np.nanmedian(legacy_noise) + data_ratio = new_data_r[bright] / legacy_data[bright] + noise_ratio = new_noise_r / legacy_noise + + summary = { + "n_exposures": record["acquire"]["n_exposures"], + "offset": [float(dy), float(dx)], + "n_bright": int(bright.sum()), + "data_ratio_median": float(np.nanmedian(data_ratio)), + "data_ratio_16_84": [ + float(np.nanpercentile(data_ratio, 16)), + float(np.nanpercentile(data_ratio, 84)), + ], + "noise_ratio_median": float(np.nanmedian(noise_ratio)), + "noise_ratio_16_84": [ + float(np.nanpercentile(noise_ratio, 16)), + float(np.nanpercentile(noise_ratio, 84)), + ], + "correlated_noise_factor_applied": record["noise"]["correlated_noise_factor"], + "psf_diagnostics": record["psf"], + } + print("[parity] ---- production parity ----") + print(json.dumps(summary, indent=2)) + (out_dir / "parity_summary.json").write_text(json.dumps(summary, indent=2)) + + +if __name__ == "__main__": + main() diff --git a/test_autoreduce/test_cache.py b/test_autoreduce/test_cache.py new file mode 100644 index 0000000..a4ddd77 --- /dev/null +++ b/test_autoreduce/test_cache.py @@ -0,0 +1,90 @@ +import json + +import pytest + +from autoreduce.acquire.cache import ExposureCache + + +def _fake_exposure(cache, target, name, size=1000): + target_dir = cache.target_dir(target) + target_dir.mkdir(parents=True, exist_ok=True) + path = target_dir / name + path.write_bytes(b"x" * size) + return path + + +class TestManifest: + def test_record_and_read_back(self, tmp_path): + cache = ExposureCache(tmp_path) + p = _fake_exposure(cache, "lens1", "a_flc.fits") + cache.record_download("lens1", [str(p)], source="mast") + assert cache.exposures_for("lens1") == [p] + manifest = json.loads(cache.manifest_path.read_text()) + assert manifest["targets"]["lens1"]["source"] == "mast" + + def test_missing_files_fail_loudly(self, tmp_path): + cache = ExposureCache(tmp_path) + p = _fake_exposure(cache, "lens1", "a_flc.fits") + cache.record_download("lens1", [str(p)], source="mast") + p.unlink() + with pytest.raises(FileNotFoundError, match="gone"): + cache.exposures_for("lens1") + + def test_unknown_target_returns_empty(self, tmp_path): + cache = ExposureCache(tmp_path) + assert cache.exposures_for("nope") == [] + + def test_incompatible_manifest_schema_fails_loudly(self, tmp_path): + cache = ExposureCache(tmp_path) + # A spike-era manifest: same filename, different schema. + cache.manifest_path.write_text(json.dumps({"target": {}, "flc_files": []})) + with pytest.raises(ValueError, match="not an ExposureCache manifest"): + cache.exposures_for("lens1") + + +class TestEviction: + def test_evict_removes_files_keeps_provenance(self, tmp_path): + cache = ExposureCache(tmp_path) + p = _fake_exposure(cache, "lens1", "a_flc.fits") + cache.record_download("lens1", [str(p)], source="mast") + cache.evict("lens1") + assert not p.exists() + manifest = cache.read_manifest() + assert manifest["targets"]["lens1"]["evicted"] + assert manifest["targets"]["lens1"]["files"] # provenance retained + assert cache.exposures_for("lens1") == [] + + def test_evict_unknown_target_raises(self, tmp_path): + with pytest.raises(KeyError): + ExposureCache(tmp_path).evict("nope") + + def test_cap_evicts_oldest_completed_first(self, tmp_path): + cache = ExposureCache(tmp_path, size_cap_bytes=2500) + for i, name in enumerate(["old", "mid", "new"]): + p = _fake_exposure(cache, name, "a_flc.fits", size=1000) + cache.record_download(name, [str(p)], source="mast") + # Distinct timestamps: rewrite downloaded_at deterministically. + manifest = cache.read_manifest() + manifest["targets"][name]["downloaded_at"] = f"2026-07-08T00:0{i}:00Z" + cache._write_manifest(manifest) + cache.mark_completed("old") + cache.mark_completed("mid") + # "new" is not completed: never evicted even over cap. + evicted = cache.enforce_cap() + assert evicted == ["old"] + assert cache.size_bytes() <= 2500 + + def test_uncapped_never_evicts(self, tmp_path): + cache = ExposureCache(tmp_path) + p = _fake_exposure(cache, "lens1", "a_flc.fits") + cache.record_download("lens1", [str(p)], source="mast") + cache.mark_completed("lens1") + assert cache.enforce_cap() == [] + assert p.exists() + + def test_references_excluded_from_size(self, tmp_path): + cache = ExposureCache(tmp_path) + refs = cache.references_dir / "references" / "hst" / "acs" + refs.mkdir(parents=True) + (refs / "flat.fits").write_bytes(b"x" * 10_000) + assert cache.size_bytes() == 0 diff --git a/test_autoreduce/test_noise.py b/test_autoreduce/test_noise.py new file mode 100644 index 0000000..ba0c9b2 --- /dev/null +++ b/test_autoreduce/test_noise.py @@ -0,0 +1,111 @@ +import numpy as np +import pytest + +from autoreduce.noise.rms import ( + assert_finite_within, + casertano_r, + empirical_background_rms, + noise_map_from, +) + + +class TestCasertanoR: + def test_shift_and_add_limit(self): + # p=1, s=1 is shift-and-add: r = 2/3, R = 1.5 + assert casertano_r(1.0, 1.0) == pytest.approx(1.5) + + def test_spike_value(self): + # The value the SLACS parity study validated against legacy noise. + assert casertano_r(0.8, 1.0) == pytest.approx(1.364, abs=1e-3) + + def test_interlacing_limit(self): + # p -> 0 at fixed s: no flux sharing, R -> 1. + assert casertano_r(1e-6, 1.0) == pytest.approx(1.0, abs=1e-3) + + def test_smaller_pixfrac_reduces_correlation(self): + assert casertano_r(0.6, 1.0) < casertano_r(0.8, 1.0) < casertano_r(1.0, 1.0) + + def test_branch_s_less_than_p(self): + # Finer output grid than the drop uses the (s/p) branch: more + # correlation than shift-and-add, R -> inf as s -> 0. + r_fine = casertano_r(1.0, 0.5) + assert r_fine > casertano_r(1.0, 1.0) + + def test_invalid_inputs_raise(self): + with pytest.raises(ValueError): + casertano_r(0.0, 1.0) + with pytest.raises(ValueError): + casertano_r(1.5, 1.0) + with pytest.raises(ValueError): + casertano_r(0.8, 0.0) + + +class TestNoiseMapFrom: + def test_background_only(self): + sci = np.zeros((4, 4)) + wht = np.full((4, 4), 25.0) + noise = noise_map_from(sci, wht, exptime=100.0) + assert noise == pytest.approx(np.full((4, 4), 0.2)) + + def test_poisson_term_adds_in_quadrature(self): + sci = np.full((2, 2), 9.0) # cps + wht = np.full((2, 2), 4.0) + noise = noise_map_from(sci, wht, exptime=1.0) + assert noise == pytest.approx(np.full((2, 2), np.sqrt(9.0 + 0.25))) + + def test_negative_sky_pixels_floor_poisson_at_zero(self): + sci = np.array([[-5.0]]) + wht = np.array([[4.0]]) + noise = noise_map_from(sci, wht, exptime=1.0) + assert noise == pytest.approx(np.array([[0.5]])) + + def test_correlated_factor_scales_linearly(self): + sci = np.full((2, 2), 1.0) + wht = np.full((2, 2), 1.0) + base = noise_map_from(sci, wht, exptime=1.0) + scaled = noise_map_from(sci, wht, exptime=1.0, correlated_noise_factor=1.364) + assert scaled == pytest.approx(1.364 * base) + + def test_zero_weight_becomes_nan_not_patched(self): + sci = np.zeros((2, 2)) + wht = np.array([[1.0, 0.0], [1.0, -1.0]]) + noise = noise_map_from(sci, wht, exptime=1.0) + assert np.isnan(noise[0, 1]) and np.isnan(noise[1, 1]) + + def test_shape_mismatch_raises(self): + with pytest.raises(ValueError): + noise_map_from(np.zeros((2, 2)), np.zeros((3, 3)), exptime=1.0) + + def test_bad_exptime_raises(self): + with pytest.raises(ValueError): + noise_map_from(np.zeros((2, 2)), np.ones((2, 2)), exptime=0.0) + + def test_sub_unity_correlated_factor_raises(self): + with pytest.raises(ValueError): + noise_map_from( + np.zeros((2, 2)), np.ones((2, 2)), exptime=1.0, + correlated_noise_factor=0.9, + ) + + +class TestAssertFiniteWithin: + def test_clean_map_passes(self): + assert_finite_within(np.ones((3, 3)), "test") + + def test_nan_inside_cutout_fails_loudly(self): + noise = np.ones((3, 3)) + noise[1, 1] = np.nan + with pytest.raises(ValueError, match="non-finite"): + assert_finite_within(noise, "test") + + def test_zero_noise_fails_loudly(self): + noise = np.ones((3, 3)) + noise[0, 0] = 0.0 + with pytest.raises(ValueError): + assert_finite_within(noise, "test") + + +def test_empirical_background_rms_recovers_sigma(): + rng = np.random.default_rng(0) + sci = rng.normal(0.0, 0.02, size=(200, 200)) + assert empirical_background_rms(sci) == pytest.approx(0.02, rel=0.05) diff --git a/test_autoreduce/test_psf_and_package.py b/test_autoreduce/test_psf_and_package.py new file mode 100644 index 0000000..d196615 --- /dev/null +++ b/test_autoreduce/test_psf_and_package.py @@ -0,0 +1,147 @@ +import numpy as np +import pytest + +from autoreduce.psf.epsf import ( + InsufficientStarsError, + build_epsf, + normalise_kernel, +) +from autoreduce.psf.stars import ( + StarSelection, + reject_crowded, + reject_edges, + reject_near, +) +from autoreduce.psf.fallback import ModelPSFUnavailableError, model_psf + + +class TestStarCuts: + def test_reject_crowded_pairs(self): + x = np.array([10.0, 12.0, 100.0]) + y = np.array([10.0, 10.0, 100.0]) + keep = reject_crowded(x, y, min_separation=5.0) + assert keep.tolist() == [False, False, True] + + def test_reject_edges(self): + x = np.array([5.0, 50.0]) + y = np.array([50.0, 50.0]) + keep = reject_edges(x, y, shape=(100, 100), margin=10) + assert keep.tolist() == [False, True] + + def test_reject_near_target(self): + x = np.array([50.0, 90.0]) + y = np.array([50.0, 90.0]) + keep = reject_near(x, y, 50.0, 50.0, radius=10.0) + assert keep.tolist() == [False, True] + + +class TestNormaliseKernel: + def test_unit_sum_and_shape(self): + psf = np.random.default_rng(0).random((61, 61)) + 1.0 + kernel = normalise_kernel(psf, (21, 21)) + assert kernel.shape == (21, 21) + assert kernel.sum() == pytest.approx(1.0) + + def test_even_shape_rejected(self): + with pytest.raises(ValueError, match="odd"): + normalise_kernel(np.ones((61, 61)), (20, 20)) + + def test_oversized_request_rejected(self): + with pytest.raises(ValueError, match="exceeds"): + normalise_kernel(np.ones((21, 21)), (61, 61)) + + def test_zero_flux_rejected(self): + with pytest.raises(ValueError, match="flux"): + normalise_kernel(np.zeros((21, 21)), (11, 11)) + + +class TestTierFailuresAreLoud: + def test_too_few_stars_raises(self): + with pytest.raises(InsufficientStarsError, match="tier 2"): + build_epsf(np.zeros((100, 100)), None, (21, 21), (61, 61)) + + def test_tier2_unimplemented_is_hard_stop(self): + with pytest.raises(ModelPSFUnavailableError, match="hard stop"): + model_psf("lens", "F814W", (21, 21), (61, 61)) + + +class TestCutout: + def _mosaic(self): + from astropy.io import fits + from astropy.wcs import WCS + + wcs = WCS(naxis=2) + wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"] + wcs.wcs.crval = [2.0, -0.1] + wcs.wcs.crpix = [100.5, 100.5] + wcs.wcs.cdelt = [-0.05 / 3600.0, 0.05 / 3600.0] + header = wcs.to_header() + header["BUNIT"] = "ELECTRONS/S" + header["EXPTIME"] = 1566.0 + data = np.arange(200 * 200, dtype=float).reshape(200, 200) + return data, header + + def test_cutout_preserves_wcs_and_metadata(self, tmp_path): + from astropy.io import fits + from astropy.wcs import WCS + + from autoreduce.package.cutout import cutout_to_fits + + data, header = self._mosaic() + out = tmp_path / "data.fits" + cut = cutout_to_fits(data, header, ra=2.0, dec=-0.1, shape=(51, 51), out_path=out) + assert cut.shape == (51, 51) + + with fits.open(out) as hdul: + out_header = hdul[0].header + assert out_header["BUNIT"] == "ELECTRONS/S" + assert out_header["EXPTIME"] == 1566.0 + scales = np.abs(np.diag(WCS(out_header).pixel_scale_matrix)) * 3600 + assert scales == pytest.approx([0.05, 0.05]) + # The cutout centre maps back to the requested sky position. + x, y = WCS(out_header).world_to_pixel_values(2.0, -0.1) + assert float(x) == pytest.approx(25.0, abs=0.51) + assert float(y) == pytest.approx(25.0, abs=0.51) + + def test_cutout_off_mosaic_fails(self, tmp_path): + from autoreduce.package.cutout import cutout_to_fits + + data, header = self._mosaic() + with pytest.raises(Exception): + cutout_to_fits( + data, header, ra=50.0, dec=50.0, shape=(51, 51), + out_path=tmp_path / "data.fits", + ) + + +def test_weight_uniformity_diagnostic(): + from autoreduce.drizzle.diagnostics import check_weight_uniformity, weight_uniformity + + flat = np.full((50, 50), 100.0) + assert weight_uniformity(flat) == pytest.approx(0.0) + rng = np.random.default_rng(1) + speckled = np.abs(rng.normal(100.0, 40.0, size=(50, 50))) + verdict = check_weight_uniformity(speckled) + assert not verdict["acceptable"] + with pytest.raises(ValueError, match="coverage"): + weight_uniformity(np.zeros((5, 5))) + + +def test_provenance_record(tmp_path): + import json + + from autoreduce.package.provenance import write_reduction_json + + path = write_reduction_json(tmp_path, {"target": {"name": "lens"}}) + payload = json.loads(path.read_text()) + assert payload["target"]["name"] == "lens" + assert "astropy" in payload["software"] + assert payload["written_at"].endswith("Z") + + +def test_mast_query_hygiene(): + from autoreduce.acquire.mast import is_direct_observation + + assert is_direct_observation("j9op01010", "10886") + assert not is_direct_observation("hst_skycell-p1322x03y02_acs_wfc_f814w_all", "--") + assert not is_direct_observation("j9op01010", "--") diff --git a/test_autoreduce/test_target_and_instruments.py b/test_autoreduce/test_target_and_instruments.py new file mode 100644 index 0000000..eb29fb8 --- /dev/null +++ b/test_autoreduce/test_target_and_instruments.py @@ -0,0 +1,78 @@ +import pytest + +from autoreduce import instruments +from autoreduce.target import TargetSpec + + +class TestTargetSpec: + def test_defaults_match_design_doc(self): + spec = TargetSpec(name="lens", ra=2.0, dec=-0.1) + assert spec.final_scale == 0.05 + assert spec.final_pixfrac == 0.8 + assert spec.cutout_shape == (281, 281) + assert spec.psf_shape == (21, 21) + assert spec.psf_full_shape == (61, 61) + + def test_yaml_round_trip(self, tmp_path): + path = tmp_path / "target.yaml" + path.write_text( + "name: slacs0008-0004\n" + "ra: 2.012333\n" + "dec: -0.068944\n" + "proposal_ids: [10886]\n" + "cutout_shape: [281, 281]\n" + "final_pixfrac: 0.6\n" + ) + spec = TargetSpec.from_yaml(path) + assert spec.proposal_ids == ("10886",) + assert spec.final_pixfrac == 0.6 + + def test_even_psf_shape_rejected(self): + with pytest.raises(ValueError, match="odd"): + TargetSpec(name="x", ra=0.0, dec=0.0, psf_shape=(20, 20)) + + def test_pixfrac_bounds(self): + with pytest.raises(ValueError): + TargetSpec(name="x", ra=0.0, dec=0.0, final_pixfrac=0.0) + with pytest.raises(ValueError): + TargetSpec(name="x", ra=0.0, dec=0.0, final_pixfrac=1.2) + + def test_dec_bounds(self): + with pytest.raises(ValueError): + TargetSpec(name="x", ra=0.0, dec=91.0) + + +class TestInstrumentRegistry: + def test_acs_wfc_registered(self): + adapter = instruments.get("acs_wfc") + assert adapter.native_scale == 0.05 + assert adapter.calibrated_suffix == "FLC" + assert adapter.reference_env_key == "jref" + + def test_unknown_key_raises_with_choices(self): + with pytest.raises(KeyError, match="acs_wfc"): + instruments.get("nircam") + + def test_double_registration_rejected(self): + with pytest.raises(ValueError): + instruments.register(instruments.ACS_WFC) + + def test_scale_ratio(self): + assert instruments.get("acs_wfc").scale_ratio(0.05) == pytest.approx(1.0) + assert instruments.get("acs_wfc").scale_ratio(0.03) == pytest.approx(0.6) + + +def test_drizzle_kwargs_single_vs_multi_exposure(): + from autoreduce.drizzle.combine import drizzle_kwargs_for + + spec = TargetSpec(name="x", ra=0.0, dec=0.0) + adapter = instruments.get("acs_wfc") + multi = drizzle_kwargs_for(spec, adapter, 4) + single = drizzle_kwargs_for(spec, adapter, 1) + assert multi["driz_cr"] and multi["median"] and multi["blot"] + # SLACS-V caveat: single exposures cannot median-combine. + assert not (single["driz_cr"] or single["median"] or single["blot"]) + assert single["final_units"] == "cps" + assert single["final_wht_type"] == "IVM" + with pytest.raises(ValueError): + drizzle_kwargs_for(spec, adapter, 0)