WFS-in-the-loop: ADU calibration, input-referred read noise, faster gain stage, closed-loop guide (2.6.0)

AGENTS.md

@@ -38,8 +38,8 @@ removals. A JOSS paper + citation remain a post-2.0 follow-up.
 | `backend.py` | Optional NumPy/CuPy array and RNG boundary, explicit host conversion, and backend convolution. NumPy is the reference/default. Owns the stack's `device` vocabulary (`"cpu"`, `"gpu"`, `"gpu:N"`, `"auto"`; parsed by `_parse_device`, resolved by `get_backend`) and `precision` vocabulary (`resolve_precision`: `"single"`/`"double"`, aliases `"float32"`/`"float64"`). A GPU `ArrayBackend` carries a `device_id`; `activate()` makes it current. Kept local rather than `aocore.Backend`: getframes draws per-frame noise from a device-native (cuRAND) stream, which `aocore.Backend.random` (host-side `Generator`) does not provide. |
 | `_cuda.py` | Private fused CuPy `ElementwiseKernel`s for the launch-bound GPU hot path: the photo + total expectation (`fused_expectation`, called from `simulate_frame`) and the whole readout (`fused_readout`, called from `digitize`). Each mirrors the separate CuPy operations it replaces op for op (same order, working precision, Python scalars rounded to it, `--fmad=false`), so seeded GPU frames are bit-identical to the unfused path, which remains the fallback for any operand it cannot take exactly. Imported only on the GPU backend; CuPy is imported lazily. |
 | `frame.py` | `Frame` container: a NumPy array (ADU) plus metadata; array-like; optional FITS export. |
-| `camera.py` | `Camera`, the main user-facing object. Orchestrates config + scene + noise into `Frame`s. Holds the RNG and high-level methods (`dark_frame`, `dark_series`, reset-correlated `nondestructive_series`, `correlated_double_sample`, `expose`, `observe`, `*_series`, `master_*`). Reset-correlated readout has one core, the private `_ramp_reads`, which walks a ramp on an arbitrary per-read interval pattern; `nondestructive_series` drives it uniformly and `correlated_double_sample` drives it as pedestal-then-signal. Put new ramp-readout modes there rather than in a second loop — the interval-scaled bias/settling/avalanche terms are easy to get subtly wrong twice. Every public method that touches device arrays is decorated with `backend._on_device`, which runs it (each step, for generator methods) with the camera's CUDA device current; decorate new ones too, or `device="gpu:N"` silently allocates on the wrong card. |
-| `calibrate.py` | Master-frame builders (`combine`) and `calibrate` reduction — the raw → reduced → truth loop (phase 1.1). |
+| `camera.py` | `Camera`, the main user-facing object. Orchestrates config + scene + noise into `Frame`s. Holds the RNG and high-level methods (`dark_frame`, `dark_series`, reset-correlated `nondestructive_series`, `correlated_double_sample`, `expose`, `observe`, `*_series`, `master_*`, plus `calibration` and the input-referred `effective_read_noise_e`). Reset-correlated readout has one core, the private `_ramp_reads`, which walks a ramp on an arbitrary per-read interval pattern; `nondestructive_series` drives it uniformly and `correlated_double_sample` drives it as pedestal-then-signal. Put new ramp-readout modes there rather than in a second loop — the interval-scaled bias/settling/avalanche terms are easy to get subtly wrong twice. Every public method that touches device arrays is decorated with `backend._on_device`, which runs it (each step, for generator methods) with the camera's CUDA device current; decorate new ones too, or `device="gpu:N"` silently allocates on the wrong card. |
+| `calibrate.py` | Master-frame builders (`combine`), `calibrate` reduction — the raw → reduced → truth loop (phase 1.1) — and `Calibration`, the ADU → input-referred electrons → photons inverse of the mean signal chain (mean master dark, `gain_e_per_adu / em_gain`, nonlinearity inverse, optional flat, scalar/curve QE) used for WFS frames. It works on NumPy and CuPy frames; keep its inverse in step with the forward model in `noise.py`. |
 | `observation.py` | `Observation` / `ObservationTruth` / `Pointing`: the time-series driver, jitter/drift/dither, per-frame truth (phase 1.2). |
 | `spectral.py` | Opt-in spectral mode: `QE`, `SED` (relative or absolute via `from_flux_density`), `Spectrum`, `SpectralBandpass`, effective-QE folding, transmission-product helpers (`product`, `from_file`/`from_product`), optional `astropy.units` coercion. |
 | `scene/` | The scene/optics layer: `Scene`, `Source` hierarchy (`PointSource`, `ExtendedSource`, `UniformIllumination`, `Catalog`; point/extended sources accept a `flux_sed` absolute SED), PSFs (`GaussianPSF`/`MoffatPSF`/`AiryPSF`/`ArrayPSF`/`EllipticalGaussianPSF`), `Telescope` (+ `Vignetting`/`RadialDistortion`), `Bandpass` (Vega `johnson` / AB `ab` ugriz·Gaia·2MASS + `Extinction`) in `photometry.py`, `Thermal` graybody background in `thermal.py`, `WCSInfo`, `LightCurve`. Renders a photon-rate map; no randomness. |
@@ -211,7 +211,7 @@ make the smallest change that fits the existing patterns; and if a change spans
 layers, respect the one-way data flow. When something is genuinely done and
 verified, say so plainly; when a step was skipped or a test fails, say *that*.
 
-## Performance notes (GPU hot path)
+## Performance notes (hot path)
 
 - Small GPU frames are launch-bound (each CuPy call costs 10–20 µs of host time on a
   slow core), so per-frame work counts CuPy calls, not FLOPs. Large frames, and the
@@ -229,6 +229,15 @@ verified, say so plainly; when a step was skipped or a test fails, say *that*.
   `backend.asarray` for a value that may be a scalar, and pass scalar distribution
   parameters through the `_CuPyGenerator` wrapper (it reuses a cached device copy),
   so a GPU frame does not upload a fresh 0-d array per call.
+- **CPU gain stage.** NumPy's `Generator.gamma` returns `Gamma(shape=0)` as exactly
+  zero *without consuming a draw* but still pays its per-element dispatch, so
+  `apply_gain_stage` samples only the lit pixels (C order) when they are a minority.
+  This is bit-identical only because zero-shape elements consume nothing;
+  `tests/test_wfs_calibration.py::test_gain_stage_matches_full_gamma_draw_bit_for_bit`
+  guards it against the installed NumPy. Do not extend the trick to Poisson: there a
+  zero mean also consumes nothing, but every pixel carries dark current. An exact
+  vectorised replica of NumPy's small-rate Poisson (for CIC) was tried and measured
+  no faster than `rng.poisson`, so it was not kept.
 - Identity arithmetic on a whole frame (`+ 0`, `/ 1`, `* 1` by a *Python* scalar) is
   skipped via `noise._is_scalar_value`; it cannot change any ADU (at most the sign
   of a zero, which rounding erases). Keep the test to Python scalars: reading a 0-d

CHANGELOG.md

@@ -6,6 +6,68 @@ to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
 
+## [2.6.0] - 2026-10-10
+
+### Added
+
+- **`Calibration`: ADU to input-referred electrons and incident photons**
+  (`gf.Calibration`, `Camera.calibration`). `Calibration.from_camera(cam,
+  exposure, n_darks=64, seed=...)` averages exposure-matched dark frames into a
+  mean master dark on the camera's device and precision; `to_electrons` applies
+  `L^-1[(ADU - dark) * gain_e_per_adu / em_gain] / flat`, inverting the
+  configured nonlinearity (closed form or Newton) by default; `to_photons`
+  divides by the scalar QE `expose` applies, a `quantum_efficiency`, the
+  `qe_curve` at `wavelength_nm`, or its effective value over a `bandpass`;
+  `saturated` flags pixels whose charge reached the input well, the output
+  register or the ADC ceiling. It works on NumPy and CuPy frames, takes an
+  optional unit-mean flat, and handles binned frames. Tests check that the mean
+  calibrated flux matches the incident photons to 0.5% for EMCCD, qCMOS, sCMOS,
+  eAPD, CCD and CMOS presets (issue #10).
+- **`Camera.effective_read_noise_e(binning=..., binning_mode=...)`**: the
+  per-frame read noise referred to the input (electrons collected in the
+  pixel): output read noise (the quadrature mean of a per-pixel map), kTC/reset
+  noise and avalanche-scaled noise, divided by the EM/avalanche gain. The
+  docstring and the new guide give the excess-noise-factor SNR.
+- **Guide: [A WFS camera in a closed loop](docs/guides/wfs-closed-loop.md)**:
+  exposure vs. frame period, ADU calibration and its accuracy, input-referred
+  noise and the excess noise factor, recommended EMCCD/qCMOS/sCMOS/eAPD/CCD
+  presets, and keeping small frames cheap.
+- `benchmarks/bench_small_frames.py`: seeded 256x256 WFS frames for one preset
+  per detector class, with before/after comparison and a SHA-256 digest of
+  seeded output for bit-identity checks; results in
+  `benchmarks/small-frame-results.md`.
+
+### Performance
+
+- **Photon-starved EMCCD and eAPD frames are faster on the CPU, with
+  bit-identical seeded output.** NumPy returns `Gamma(shape=0)` as exactly zero
+  without consuming a draw, but still visits every such pixel; the gain stage
+  now samples only the pixels that hold charge when they are a minority (a
+  wavefront sensor's dark background), and skips the `* 1` shape scaling at
+  `F = sqrt(2)`. Same draws, same order, same generator state afterwards;
+  seeded frames of every preset, in both precisions, are bit-identical to 2.5.0.
+  On one Neoverse-N1 core (256x256 Shack-Hartmann spots, 1 ms, in-process A/B):
+  the Gamma draw 1.43 → 0.78 ms; whole frames `andor_ocam2k` 6.57 → 5.95 ms
+  (1.10x), `leonardo_saphira` 4.79 → 4.29 ms (1.12x), `first_light_imaging_cred_one`
+  7.16 → 6.63 ms (1.08x); presets without a gain stage are unchanged. The CPU
+  chain is otherwise at its sampler floor (about 85% of a frame is NumPy's
+  Poisson, Gamma and Gaussian draws). An exact vectorised replica of NumPy's
+  small-rate Poisson for the CIC draw was tried and measured no faster, so it
+  was not kept. See `benchmarks/small-frame-results.md`.
+
+### Documented
+
+- The 1--2% calibrated-flux deficit seen in the dithered Shack-Hartmann
+  prototype is not a deterministic offset of the signal chain: averaged over
+  300 frames the prototype recovers 1.002 (OCAM2K), 0.999 (ORCA-Quest 2) and
+  1.001 (generic sCMOS) of the incident flux, while a single 288x288 frame
+  scatters by 1.2%, 1.0% and 3.3%. The guide lists the deterministic offsets a
+  pipeline can introduce (median master dark, scalar vs. spectral QE,
+  nonlinearity, saturation, unflattened fixed patterns) and their remedies.
+- The Nüvü HNü 240 / 128 presets cap the *amplified* charge at their 500 / 800
+  e- `full_well_e` (no `output_full_well_e`), so at EM gain 1000 a pixel holding
+  one electron saturates; the guide warns about it (#11).
+
 ## [2.5.0] - 2026-10-07
 
 ### Added

README.md

@@ -127,6 +127,10 @@ for the methodology.
 - **Calibration & ground truth** — master bias/dark/flat builders and a
   `calibrate` reduction that closes the
   **[raw → reduced → truth loop](https://jacotay7.github.io/getframes/guides/calibration/)**.
+- **Wavefront-sensor cameras** — `Calibration` turns ADU back into
+  input-referred electrons and photons, `Camera.effective_read_noise_e` gives
+  the input-referred noise, and a guide covers
+  **[a WFS camera in a closed loop](https://jacotay7.github.io/getframes/guides/wfs-closed-loop/)**.
 - **Observations** — `Observation` drives time series with jitter, drift, dither
   and persistence, carrying per-frame
   **[truth](https://jacotay7.github.io/getframes/guides/time-series/)**.

benchmarks/bench_small_frames.py

@@ -0,0 +1,192 @@
+# SPDX-License-Identifier: MIT
+"""Benchmark small, frequent wavefront-sensor frames on the CPU.
+
+A closed-loop wavefront-sensor study calls :meth:`getframes.Camera.expose` once
+per loop step on a small frame, so per-frame latency, not throughput on large
+detectors, is what limits it. This script times seeded ``expose`` calls on a
+256x256 Shack-Hartmann-like spot pattern (16x16 subapertures of 16x16 pixels,
+one Gaussian spot each) for one preset of each WFS detector class, at a
+photon-starved and a bright flux, in two call styles:
+
+* ``default``: ``expose(rate, t, seed=s, include_truth=False)``;
+* ``workspace``: the same with a reused :class:`getframes.DetectorWorkspace` and
+  a caller-owned ``out`` frame.
+
+It also records a SHA-256 digest of a few seeded frames per configuration, so two
+revisions can be compared for bit-identical output (``--compare``). Timings are
+process CPU time per frame over interleaved, warmed blocks (the chain is
+single-threaded); they are machine dependent and are not part of the test gate.
+Each case reports the median and the best (minimum) block; on a shared host,
+compare the best blocks, which other load can only inflate.
+
+Run from the repository root::
+
+    python benchmarks/bench_small_frames.py --output small-frame-results.json
+    python benchmarks/bench_small_frames.py --compare before.json after.json
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import platform
+import statistics
+import time
+from pathlib import Path
+from typing import Any
+
+import numpy as np
+
+import getframes as gf
+
+SHAPE = (256, 256)
+SUBAPERTURE_PX = 16
+EXPOSURE_S = 1e-3
+#: One preset per detector class a WFS is built from.
+PRESETS = (
+    ("EMCCD", "andor_ocam2k"),
+    ("EMCCD", "andor_ixon_ultra_888"),
+    ("qCMOS", "hamamatsu_orca_quest_2"),
+    ("sCMOS", "generic_scmos"),
+    ("eAPD", "first_light_imaging_cred_one"),
+)
+#: Photons per subaperture per frame: photon-starved and bright.
+FLUXES = (100.0, 2000.0)
+DIGEST_SEEDS = (0, 1, 2)
+
+
+def spot_pattern(photons_per_subaperture: float) -> np.ndarray[Any, np.dtype[np.float64]]:
+    """A 256x256 photon map of one Gaussian spot (sigma 1.5 px) per subaperture."""
+    coordinates = np.arange(SUBAPERTURE_PX) - (SUBAPERTURE_PX - 1) / 2
+    spot = np.exp(-(coordinates[:, None] ** 2 + coordinates[None, :] ** 2) / (2 * 1.5**2))
+    spot *= photons_per_subaperture / spot.sum()
+    reps = (SHAPE[0] // SUBAPERTURE_PX, SHAPE[1] // SUBAPERTURE_PX)
+    return np.tile(spot, reps)
+
+
+def make_camera(preset: str) -> gf.Camera:
+    return gf.Camera.from_preset(preset).with_config(resolution=SHAPE, roi=None)
+
+
+def _block_ms(function: Any, calls: int) -> float:
+    # Process CPU time, not wall time: the chain is single-threaded, and CPU time
+    # is far less sensitive to other load on a shared host.
+    start = time.process_time_ns()
+    for index in range(calls):
+        function(index)
+    return (time.process_time_ns() - start) / calls / 1e6
+
+
+def bench_case(preset: str, flux: float, calls: int, repetitions: int) -> dict[str, Any]:
+    camera = make_camera(preset)
+    rate = spot_pattern(flux) / EXPOSURE_S
+    workspace = gf.DetectorWorkspace()
+    out = np.empty(SHAPE, dtype=np.uint32)
+
+    def default(seed: int) -> None:
+        camera.expose(rate, EXPOSURE_S, seed=seed, include_truth=False)
+
+    def reused(seed: int) -> None:
+        camera.expose(
+            rate, EXPOSURE_S, seed=seed, include_truth=False, workspace=workspace, out=out
+        )
+
+    for seed in range(10):
+        default(seed)
+        reused(seed)
+    samples: dict[str, list[float]] = {"default": [], "workspace": []}
+    variants = [("default", default), ("workspace", reused)]
+    for repetition in range(repetitions):
+        for name, function in variants if repetition % 2 == 0 else variants[::-1]:
+            samples[name].append(_block_ms(function, calls))
+
+    digest = hashlib.sha256()
+    for seed in DIGEST_SEEDS:
+        frame = camera.expose(rate, EXPOSURE_S, seed=seed)
+        digest.update(np.ascontiguousarray(frame.data).tobytes())
+        assert frame.truth is not None
+        digest.update(np.ascontiguousarray(frame.truth.mean_electrons).tobytes())
+    return {
+        "preset": preset,
+        "photons_per_subaperture": flux,
+        "default_ms": statistics.median(samples["default"]),
+        "workspace_ms": statistics.median(samples["workspace"]),
+        "default_best_ms": min(samples["default"]),
+        "workspace_best_ms": min(samples["workspace"]),
+        "default_samples_ms": samples["default"],
+        "workspace_samples_ms": samples["workspace"],
+        "seeded_digest": digest.hexdigest(),
+    }
+
+
+def run(calls: int, repetitions: int) -> dict[str, Any]:
+    results = []
+    for sensor, preset in PRESETS:
+        for flux in FLUXES:
+            case = bench_case(preset, flux, calls, repetitions)
+            case["class"] = sensor
+            results.append(case)
+            print(
+                f"{sensor:6s} {preset:30s} {flux:7.0f} ph/subap  "
+                f"default {case['default_ms']:6.2f} ms  workspace {case['workspace_ms']:6.2f} ms",
+                flush=True,
+            )
+    return {
+        "schema_version": 1,
+        "getframes": gf.__version__,
+        "numpy": np.__version__,
+        "platform": platform.platform(),
+        "processor": platform.processor() or platform.machine(),
+        "shape": list(SHAPE),
+        "exposure_s": EXPOSURE_S,
+        "calls_per_block": calls,
+        "repetitions": repetitions,
+        "results": results,
+    }
+
+
+def compare(before_path: Path, after_path: Path) -> str:
+    """A Markdown table of before/after best blocks and whether seeded output matches."""
+    before = json.loads(before_path.read_text(encoding="utf-8"))
+    after = json.loads(after_path.read_text(encoding="utf-8"))
+    lines = [
+        f"Before: getframes {before['getframes']}; after: getframes {after['getframes']} "
+        f"({after['processor']}, NumPy {after['numpy']}, {after['shape'][0]}x{after['shape'][1]}).",
+        "",
+        "| Class | Preset | Photons/subap | Before (ms, best) | After (ms, best) | Speed-up | "
+        "After, workspace+out (ms, best) | Seeded output |",
+        "| --- | --- | ---: | ---: | ---: | ---: | ---: | --- |",
+    ]
+    old = {(r["preset"], r["photons_per_subaperture"]): r for r in before["results"]}
+    for row in after["results"]:
+        ref = old.get((row["preset"], row["photons_per_subaperture"]))
+        if ref is None:
+            continue
+        same = "identical" if ref["seeded_digest"] == row["seeded_digest"] else "CHANGED"
+        lines.append(
+            f"| {row['class']} | {row['preset']} | {row['photons_per_subaperture']:.0f} | "
+            f"{ref['default_best_ms']:.2f} | {row['default_best_ms']:.2f} | "
+            f"{ref['default_best_ms'] / row['default_best_ms']:.2f}x | "
+            f"{row['workspace_best_ms']:.2f} | {same} |"
+        )
+    return "\n".join(lines) + "\n"
+
+
+def main() -> None:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--calls", type=int, default=40)
+    parser.add_argument("--repetitions", type=int, default=7)
+    parser.add_argument("--output", type=Path)
+    parser.add_argument("--compare", nargs=2, type=Path, metavar=("BEFORE", "AFTER"))
+    arguments = parser.parse_args()
+    if arguments.compare:
+        print(compare(*arguments.compare), end="")
+        return
+    payload = run(arguments.calls, arguments.repetitions)
+    if arguments.output is not None:
+        arguments.output.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8")
+
+
+if __name__ == "__main__":
+    main()