research: plan the Keck NIRC2 AO reduction workflow (SHARP context)

[Jammy2211]

Overview

Plan how PyAutoReduce reduces Keck adaptive-optics (NIRC2) imaging of strong lenses into modeling-ready datasets for PyAutoLens/PyAutoGalaxy, alongside the existing HST (ACS/WFC3) and JWST (NIRCam) phases. The scientific grounding comes from the SHARP programme (Strong-lensing at High Angular Resolution Programme) papers, which are the reference practice for AO strong-lens imaging. Deliverable is a design/plan document (docs/design/keck_ao.md), not code — implementation phases become follow-up prompts.

This run is --auto at effective level supervised (research cap). Plan below per the autonomy contract; judgment forks will be batched as comments here.

Plan

• Read the SHARP survey papers (Fassnacht, Chen, Vegetti, McKean et al.) plus the standard Keck/NIRC2 reduction references (NIRC2 distortion solutions, the KAI pipeline) and extract what AO imaging demands of a lensing reduction — above all the time-variable, poorly-known AO PSF.
• Survey existing NIRC2 reduction tooling and decide wrap-vs-implement per PyAutoReduce's "default pipelines first" boundary.
• Map the AO reduction onto PyAutoReduce's stage architecture (acquire/align/drizzle/noise/psf/package + instruments/ adapters); identify genuinely new stages (sky subtraction, PSF estimation strategy) vs per-stage deltas.
• Choose a validation anchor: a SHARP lens with public Keck Observatory Archive (KOA) NIRC2 data, with defined quality/parity checks.
• Write docs/design/keck_ao.md (per-stage deltas vs the HST design, following jwst.md) and update docs/design/roadmap.md.
• Enumerate implementation phases as follow-up PyAutoMind prompts.

Sequencing note: PyAutoReduce is currently claimed by refactor-post-phase3 (#8, in-place). The research legs are read-only and proceed now, with findings and the draft design content landing on this issue; the repo-write leg (keck_ao.md + roadmap + PR) is parked until that task ships. If the refactor changes the stage/module layout, the design doc adopts the post-refactor names before writing.

Detailed implementation plan

Work Classification

Library (PyAutoReduce) — design doc only, no source changes in this task.

Affected Repositories

• PyAutoReduce (primary)

Branch Survey

Repository	Current Branch	Dirty?
./PyAutoReduce	refactor/post-phase3-consolidation	clean (claimed in-place by refactor-post-phase3)

Suggested branch: feature/keck-ao-reduction-plan (created only for the write leg, in a task worktree — the main checkout stays with the refactor)
Worktree root: ~/Code/PyAutoLabs-wt/keck-ao-reduction-plan/ (deferred until the write leg unblocks)

Implementation Steps

1. Literature pass — SHARP series papers (the AO-imaging lens-modelling and H0 papers carry the most reduction detail, esp. PSF handling); NIRC2 geometric-distortion solutions (Yelda et al. 2010, Service et al. 2016); the KAI (Keck AO Imaging) pipeline as the community reduction baseline. Extract per-stage requirements: calibration (darks/flats/bad-pixel), sky subtraction, distortion correction, frame registration + coaddition, noise-map construction, PSF strategy.
2. Requirements → architecture mapping — table of PyAutoReduce stages vs AO needs: acquire (KOA archive access, e.g. PyKOA/astroquery vs MAST), align/drizzle (dewarp + shift-and-add vs drizzle of dithered frames), noise (read noise + Poisson + sky-dominated regime), psf (the AO-specific problem: field-star PSF, PSF-as-free-model, reconstruction — what the output contract psf.fits/psf_full.fits should mean when the PSF is uncertain), package (unchanged al.Imaging.from_fits contract).
3. Wrap-vs-implement per stage — deviation from community-standard tooling must be justified by a lensing requirement (repo boundary rule).
4. Validation anchor — pick a SHARP lens with public KOA NIRC2 narrow-camera data; define checks (astrometric parity vs HST imaging of the same lens, photometric sanity, noise-map statistics, PSF sanity/fit quality).
5. Write the design doc — docs/design/keck_ao.md mirroring jwst.md's "per-stage deltas vs the HST design" pattern + roadmap update. (Parked leg — needs the PyAutoReduce claim.)
6. Follow-up prompts — phase decomposition for implementation (e.g. instruments/nirc2 adapter, sky stage, PSF stage extension).

Key Files

• docs/design/keck_ao.md — the deliverable (new)
• docs/design/roadmap.md — ground-based/AO section update
• docs/design/hst_acs_pipeline.md — the baseline design the deltas are written against
• docs/design/jwst.md — the per-stage-delta pattern to follow
• autoreduce/instruments/ — adapter pattern the plan must respect

Original Prompt

Click to expand starting prompt

Plan the Keck adaptive-optics (NIRC2 AO) data-reduction workflow

Type: research
Target: PyAutoReduce
Difficulty: medium
Autonomy: supervised
Priority: normal
Status: formalised

Original request: "For PyAutoReduce can you plan out how we reduce keck-AO data, read SHARP papers for scientific context."

Plan the Keck adaptive-optics (NIRC2 AO) imaging data-reduction workflow for PyAutoReduce — how raw Keck AO data becomes modeling-ready datasets for PyAutoLens/PyAutoGalaxy, alongside the existing HST (ACS/WFC3) and JWST (NIRCam) phases. Read the SHARP (Strong-lensing High Angular Resolution Programme) papers for scientific context on AO strong-lens imaging — PSF handling (AO PSFs are time-variable and poorly known, unlike HST/JWST), sky subtraction, distortion correction, and coaddition. Deliverable is a design/plan document, not code.

[Jammy2211]

Research complete — findings + draft docs/design/keck_ao.md

The literature and tooling passes are done. Below: (1) what the SHARP papers and the NIRC2 tooling landscape actually say, (2) the full draft design doc (per-stage deltas vs the HST design, jwst.md pattern), (3) proposed phasing, and (4) three batched decisions (recommendations stated; defaulting to them unless overridden). The repo write leg stays parked on refactor-post-phase3 (#8).

1. Findings

The SHARP series (what each paper contributes to reduction design)

Paper	System / role	Reduction-relevant content
SHARP I — Lagattuta et al. 2012	B1938+666 IR Einstein ring, K′+H, narrow camera	Canonical extended-source AO lens observation: LGS-AO with R=15 tip-tilt star at 18″; K′ 14,760 s / H 6,600 s; reduction "following Auger et al. 2011"; PSF strategy: multiple PSF-star observations taken at different times, selected between by Bayesian evidence during lens modelling
SHARP II — Hsueh et al. 2016	edge-on disc flux anomaly	modelling-side; confirms pipeline description
SHARP III — Chen et al. 2016	RXJ1131−1231, first AO time-delay cosmology	Wide camera (39.7 mas/pix, 40″×40″), 61 × (6 coadd × 10 s) = 3,660 s, tip-tilt star R=15.8 at 54.5″; the iterative PSF-reconstruction method for lensed AGN: initialize from any star (or an AGN image), decompose PSF into analytic-Gaussian global + fine-structure correction grid grown gradually (×1.2 per step), curvature regularization with evidence-optimized strength, nested-loop convergence; PSF uncertainty propagated into the noise map
SHARP IV — Hsueh et al. 2017	B0712+472	modelling-side
Chen et al. 2019 — A SHARP view of H0LiCOW	H0 from 3 AO lenses	The mature pipeline statement: "flat-field correction, sky subtraction, correction of the optical distortion, coaddition"; narrow 9.942 mas/pix resampled to 10 mas (wide 39.686 → 40 mas); 2×2 binning to 20 mas for modelling efficiency; 8 PSF-reconstruction iterations for PG1115 → FWHM 0.07″; AO models fit "down to the noise level"
SHARP VIII — Chen et al. 2022	J0924+0219	same practice, most recent

Supporting reduction craft (the Auger-method lineage, e.g. Marshall et al. 2007, SWELLS Treu et al. 2011): dithered exposures; sky from temporally adjacent object-masked frames, masks updated iteratively; B-spline residual background removal; geometric distortion corrected and frames combined with drizzle (Fruchter & Hook 2002).

The AO-specific physics the design must respect: the PSF is time-variable and spatially variable, there is usually no star inside the ~10″ narrow-camera field, and K′ sky varies on minutes timescales — so PSF and sky handling, trivial for HST/JWST, become first-class stages.

Tooling landscape (the wrap-vs-implement input)

• KAI (Keck AO Imaging DRP, Lu group) is the community-standard NIRC2/OSIRIS imaging pipeline — but it is Python 2.7 + IRAF/PyRAF, 32-bit. It is not a viable wrap target for a modern stack. This is the decisive difference from HST (drizzlepac) and JWST (jwst): there is no maintained, importable default pipeline for NIRC2 AO imaging.
• keckcode (Fassnacht) — the SHARP team's own Python code; research-grade, MIT-licensed; useful as a reference implementation to read, not a dependency.
• Distortion solutions are published lookup tables: Yelda et al. 2010 (pre-servicing) and Service et al. 2016 (post-2015-04-13), distributed via the Keck dewarp page. Observation date selects the solution.
• Archive access: PyKOA (NExScI, TAP-based) gives programmatic access to KOA raw level-0 science + calibration frames (level-1 quick-look products are explicitly not science-grade). All SHARP data are past the 18-month proprietary window.

Wrap-vs-implement verdict: implement the NIRC2 stages natively (numpy/astropy/photutils/scipy — every operation is a simple array op: dark/flat/sky/dewarp/coadd), justified under the "default pipelines first" rule by the absence of any maintained default; use the drizzle PyPI package (STScI-maintained — the same resampling engine inside JWST resample) for the dewarp+combine step, which keeps the drizzled-PSF invariant and the Casertano correlated-noise machinery intact. Cross-check against KAI outputs is optional (py2.7 conda env cost) — the SHARP-paper numbers and internal closures are the primary validation.

2. Draft docs/design/keck_ao.md

# Keck NIRC2 AO — per-stage deltas vs the HST design

Phase 4 (planned). The first ground-based instrument, which forces two new
seams: an **acquire backend** (KOA/PyKOA vs MAST) and two genuinely new
stages — **calibrate** and **sky** — that HST/JWST get for free from their
level-2 products. Scientific reference practice: the SHARP programme
(Lagattuta et al. 2012; Chen et al. 2016, 2019, 2022; reduction lineage
Auger et al. 2011). There is no maintained community pipeline to wrap (KAI
is Python 2.7 + IRAF), so stages are implemented natively and validated
against SHARP published numbers.

| Stage | Delta vs HST |
|-------|--------------|
| spec | new per-target fields: `camera` (narrow 9.942 mas / wide 39.686 mas), `ao_mode` (LGS/NGS + tip-tilt star), `psf_star` (coords or "in-field" or "none"), KOA program/semester IDs |
| acquire | **KOA via PyKOA** (TAP), not MAST: raw level-0 science frames of the target **plus the night's calibrations** (darks matched to ITIME×COADDS, dome/sky flats, bad-pixel inputs) **plus PSF-star frames** from the same program/night. Reference files = the two published distortion solutions (Yelda 2010 pre-2015-04-13, Service 2016 after), vendored + versioned in `reduction.json`. Cache semantics unchanged |
| calibrate | **new stage** (ground-based only): DN→e⁻ via gain, dark subtraction, flat-field, bad-pixel mask. Frames are ITIME×COADDS; read noise scales 1/√COADDS. Pure numpy/astropy |
| sky | **new stage**, the big NIR delta: running sky = median of N temporally adjacent, object-masked frames (K′ sky varies on minutes); iterative object-mask update; optional B-spline residual background (Auger-method). Per-frame sky level recorded — it feeds the noise stage |
| align+combine | dewarp (distortion lookup tables) **through the `drizzle` package** onto the final grid — distortion handled as the pixel mapping, exactly as drizzlepac does for ACS. Registration by cross-correlation on the target (header dither offsets are approximate). Final scales: **10 mas** narrow / **40 mas** wide (SHARP convention); optional 2×2 binning to 20 mas as a packaging dial (Chen 2019). Sigma-clip outlier rejection at combine handles cosmic rays/bad pixels |
| noise | construct from first principles (no pipeline ERR exists): per-frame variance = sky + dark + read (per-coadd-corrected) + source Poisson, propagated through the coadd coverage/weight maps, × the same Casertano R (drizzle correlates pixels identically). Blank-sky closure check as for HST/JWST. Chen-2016's PSF-uncertainty noise term is a *modelling-stage* extension, out of reduction scope |
| psf | **redefined contract for AO — the PSF is provisional, never final.** Tier A (default, extended-source lenses): reduce the PSF-star exposures through the *identical* pipeline onto the same grid; ship **all candidates** (`psf_0.fits`, `psf_1.fits`, …) because SHARP I selects between epochs by Bayesian evidence *in modelling*. Tier B: in-field ePSF (existing photutils machinery) when the field has stars (wide camera). Tier C (lensed AGN): initial PSF only; iterative reconstruction (Chen 2016) / STARRED / PSFr is modelling-stage. `reduction.json` gains `psf_provisional: true`, per-candidate FWHM/epoch/tip-tilt-distance, and an AO-conditions block (Strehl estimate, seeing, airmass) |
| package | unchanged `al.Imaging.from_fits` contract; BUNIT e⁻/s; headers intact |

## Validation anchor — B1938+666 (SHARP I)

K′ narrow-camera imaging of the canonical IR Einstein ring (the Vegetti et
al. 2012 substructure-detection system). Acceptance bar is
**"internally consistent + parity with published SHARP measurements"** (the
JWST-phase precedent — there is no legacy PyAuto Keck dataset to diff):

1. Internal closures: blank-sky vs noise-map prediction; coverage/WHT
   uniformity over the cutout; masked-pixel policy.
2. PSF core FWHM in the ~65–70 mas range reported by SHARP.
3. Astrometric parity of ring features vs HST imaging of the same system.
4. End-to-end PyAutoLens fit reproduces the published Einstein radius
   within statistical errors — science invariance, as for SLACS.

Stretch anchor (AGN/wide-camera regime): RXJ1131−1231 (SHARP III).

## Non-goals

PSF reconstruction during modelling (Chen 2016 / STARRED / PSFr — a
PyAutoLens-side concern); OSIRIS imager; NGS-only archival programs beyond
the anchor; spectroscopy.

3. Proposed phasing (follow-up prompts, issued one at a time per queue discipline)

• K1 — acquire + calibrate + sky: KOA/PyKOA adapter behind an acquire-backend seam, calibrate + sky stages, on B1938+666 raw frames.
• K2 — dewarp + combine + noise: distortion-as-drizzle-mapping, registration, coadd, first-principles noise map + closures.
• K3 — PSF tiers + package + validation: PSF-candidate products, provenance extensions, the four-point anchor validation incl. the PyAutoLens fit.

4. Batched decisions (supervised checkpoint — proceeding on the recommendations unless redirected)

1. Combine backend: drizzle package (recommended) vs pure-interpolation dewarp+stack (SHARP's own pipeline style). Drizzle keeps the correlated-noise (R) and drizzled-PSF machinery shared with HST/JWST; the SHARP-parity counterargument is that their products were interpolation-combined, but our acceptance bar is internal consistency + published-number parity, not pixel-level reproduction.
2. Validation anchor: B1938+666 (recommended) — extended-ring regime matches PyAutoLens/PyAutoGalaxy's primary use and the PSF-star strategy; RXJ1131 as stretch. Alternative: start with RXJ1131 wide-camera (easier PSF from in-field stars, but AGN-dominated and the wrong regime for the main user base).
3. Sequencing: write keck_ao.md + roadmap update after refactor-post-phase3 ships (recommended), adopting post-refactor module names; K1 implementation follows as its own prompt. Alternative: take a parallel docs-only claim now (new file, low collision risk) — declined by default since the refactor is actively consolidating stage layout.

Sources beyond those linked above: KAI docs, PyKOA NIRC2 guide.

[Jammy2211]

Design delivered: docs/design/keck_ao.md merged to main in #12 (44a8d1a), implementing the K1–K3 plan drafted here with the user-adopted decisions (drizzle backend, B1938+666 anchor). Validated on the SHARP I discovery data — closure 0.883, PSF FWHM 72 mas. Follow-ups (astrometric/Einstein-radius acceptance, wide camera, CR rejection at combine) are tracked as keck_ao.md open items for their own prompts. Closing.