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

[jacotay7]

Closes #10

Summary

• gf.Calibration (and Camera.calibration(...)): Calibration.from_camera(cam, exposure, n_darks=64, seed=...) averages seeded, exposure-matched dark frames into a mean master dark on the camera's device and in its precision.
  • to_electrons(frame) computes L^-1[(ADU - dark) * gain_e_per_adu / em_gain] / flat. It inverts the configured nonlinearity by default: closed form for the single-parameter model, Newton for the polynomial one.
  • to_photons(frame) divides by one of: the scalar QE expose applies (the default), quantum_efficiency=, the qe_curve at wavelength_nm=, or effective_qe over a bandpass=.
  • saturated(frame) flags pixels at the input-well, output-register or ADC ceiling.
  • It works on NumPy and CuPy frames (a host calibration applied to a device frame computes on the device), takes an optional unit-mean flat, handles binned frames, and accepts user-supplied darks for real data.
• Camera.effective_read_noise_e(binning=, binning_mode=) sums read, kTC/reset and avalanche-scaled noise in quadrature and divides by the EM/avalanche gain. For a per-pixel sCMOS map it uses the quadrature mean of the map. The docstring and guide give the excess-noise-factor SNR S / sqrt(F²(S+D) + σ_eff²). OCAM2K: 0.36 e-.
• Faster gain stage, bit-identical output. NumPy's Generator.gamma returns 0 for a zero shape without consuming a draw, but it still visits every such pixel. The CPU gain stage now gathers only the lit pixels when they are a minority. The draws, their order and the generator state afterwards are unchanged.
• Guide: "A WFS camera in a closed loop" (docs/guides/wfs-closed-loop.md, in the nav): exposure vs. frame period, ADU calibration and its accuracy per detector class, input-referred noise and ENF, recommended EMCCD/qCMOS/sCMOS/eAPD/CCD presets with effective read noise, and how to keep frames cheap.
• benchmarks/bench_small_frames.py plus benchmarks/small-frame-results.{md,json}. Version 2.6.0, CHANGELOG dated 2026-10-10. I bumped __about__.py because the task asked for it; AGENTS.md normally leaves the version to the release step.

Evidence

Calibration accuracy

tests/test_wfs_calibration.py uses a 64x64 Shack-Hartmann spot pattern with 1000 photons per spot, averaged over 150 frames, with a 64-frame mean dark. Mean calibrated flux vs. incident photons:

• Within 0.5%: andor_ocam2k, andor_ixon_ultra_888 (EMCCD), hamamatsu_orca_quest_2 (qCMOS), generic_scmos (sCMOS), generic_eapd, leonardo_saphira (eAPD), generic_ccd, generic_cmos.
• Within 1%: first_light_imaging_cred_one, with a flat and a per-frame background from unlit pixels.

Further tests cover:

• the flat removing PRNU on a 9-pixel spot;
• linearisation removing the nonlinearity bias, for both models;
• a median master dark reading an EMCCD photon-starved flux ~8% high;
• the QE options, saturation flags, seeded and single-precision darks, on-chip binning, validation, and the effective read noise;
• the effective read noise against measured bias-frame noise (sCMOS, 3%).

Cause of the 1–2% flux offset

It is not a deterministic bias in the signal chain. I reran the prototype's own SHMakewfs camera path: 288x288, nph=100 per subaperture, 64-frame mean dark, /em_gain/QE.

• Over 300 frames it recovers 1.0020 ± 0.0007 (OCAM2K), 0.9993 ± 0.0006 (ORCA-Quest 2) and 1.0010 ± 0.0019 (generic sCMOS) of the incident flux.
• A single frame scatters by 1.2%, 1.0% and 3.3% respectively. The main contributors are EM excess noise and read noise summed over 83k pixels.
• So a 1–2% deficit seen on a few frames is per-frame scatter.

Every stage of the chain is linear in expectation: Poisson, Gamma mean nG, zero-mean read noise, and rounding with ≥0.5 ADU of read noise. The guide lists the deterministic offsets a pipeline can introduce, and Calibration handles each one:

• median master dark: biased by the CIC and hot-pixel mean. from_camera always averages.
• scalar vs. spectral QE: the ORCA-Quest 2 scalar is 0.85, but its curve gives 0.50 at 700 nm. Use wavelength_nm= or bandpass=.
• nonlinearity: 0.67% at 2/3 of full well on generic_scmos. It is inverted by default.
• saturation: the OCAM2K output register caps a pixel at 467 input e- at gain 600. saturated() flags these pixels.
• unflattened PRNU / gain maps: pass flat=.

Found along the way

• nuvu_hnu_240 / nuvu_hnu_128_omega saturate at one electron. Their full_well_e of 500/800 e- also caps the amplified charge, because output_full_well_e is unset. Calibrated flux comes out at 3.4% of the truth. The Nüvü datasheet gives no well depth, so I did not invent a value. Opened Nüvü EMCCD presets saturate at one input electron (full_well_e doubles as the post-gain ceiling) #11; the guide carries a warning, and the tests use andor_ixon_ultra_888 instead.
• The C-RED One's 69 ADU per-frame common mode is about 6 input e- on every pixel and dominates a frame's summed flux. It is zero-mean; the guide explains how to remove it.

Benchmark (before/after)

Setup: 256x256 spot frames, 1 ms exposure, seeded, include_truth=False, one pinned Ampere Neoverse-N1 core (cfl-test-bench, NumPy 2.5.3), process CPU time.

In-process A/B (old and new gain stage alternated, best of 15 blocks):

Preset	ph/subap	2.5.0	2.6.0
andor_ocam2k	100	6.57 ms	5.95 ms	1.10x
andor_ocam2k	2000	7.36	6.91	1.07x
leonardo_saphira	100	4.79	4.29	1.12x
first_light_imaging_cred_one	100	7.16	6.63	1.08x
generic_scmos (control, no gain stage)	100	4.32	4.32	1.00x

• The Gamma draw alone: 1.43 → 0.78 ms (1.82x) with 14% of pixels lit.
• Cross-process: OCAM2K 6.49 → 5.84 ms. With workspace=/out= it is 5.05 ms, 1.28x vs. the 2.5.0 default call.
• Seeded output is bit-identical. I compared SHA-256 digests of frame and truth for every benchmark case, and for all 22 presets × {double, single} × two fluxes × three seeds.

Reading this honestly:

• About 85% of a frame is NumPy's Poisson/Gamma/Gaussian samplers. Under the "same config + inputs + seed → same frame" guarantee they cannot be swapped or threaded, so the CPU chain is now at its sampler floor.
• A frame here is 4–7 ms of CPU, not 40 ms. On the i7-10700 the prototype's whole noisy-measurement step took 14 ms, of which expose was about 5 ms.
• I also wrote an exact, bit-identical vectorised replica of NumPy's small-rate Poisson for the CIC draw. It measured no faster on x86 or Arm, so it is not included; AGENTS.md records this.
• For a large speed-up, use device="gpu".

Quality gate

• ruff check . ✔, ruff format --check . ✔, mypy ✔ (strict, 28 files).
• pytest -m "not gpu": 616 passed.
• mkdocs build --strict ✔. The guide's code blocks were executed.
• GPU (pytest -m gpu, Quadro P620): 24 passed, including the two new calibration tests (device frames, CPU/GPU statistical parity, host-calibration-on-device-frame). The pre-existing test_gpu_supports_full_detector_artifact_chain_and_sensor_families hits OOM in the full GPU run because the shared 2 GB card has ~1.75 GB in use by other processes. It passes when run alone on both this branch and main.

🤖 Generated with Claude Code

https://claude.ai/code/session_01NhyUAKxoAf5ymTWJx5WiiF