Weight wfs_input_opd_rms_m by the pupil and remove piston (2.0.0)

AGENTS.md

@@ -87,6 +87,16 @@ write a failing integration test/design note, use the conditional gates in
 - Intensities, not fields, are summed over incoherent wavelengths, modulation
   points, finite-source samples, and sodium slices.
 - Cropping reports lost flux; it does not renormalize it away.
+- Wavefront metrics follow aocore CONVENTIONS 4.1. Frame metadata
+  `wfs_input_opd_rms_m` is the pupil-intensity-weighted, piston-removed RMS of
+  the input OPD on the input grid (weights from `WavefrontSensor._rms_weights`:
+  the analytic pupil evaluated on `input.shape`, or a custom mask area-averaged
+  from the engine's `configured_pupil`); `wfs_input_opd_rms_unweighted_m` is the
+  whole-grid quadratic mean with piston kept. Both are reduced on the device and
+  cross in the one `backend.scalars` batch, so do not swap in `aocore.rms` or
+  `aocore.rms_unweighted` there: they return host floats and would add a
+  synchronization each. Any other RMS-like key must say its variant in its
+  name (`_unweighted`, `_tiptilt_removed`).
 - The intended top-level API is `load_config`, `WavefrontSensor`, and `simulate`.
   Keep other implementation objects out of `makewfs.__init__` unless an API review
   explicitly accepts them.
@@ -153,8 +163,13 @@ Follow the target layout in `ROADMAP.md`:
   rotation and rectangular grids on the selected backend, none of which
   `aocore.Pupil` models. `backend.ArrayBackend` takes a dtype per array, explicit
   FFT workers and `ndimage` helpers, which `aocore.Backend` does not, and its
-  centred FFTs keep the `fftshift` convention. `sampling.block_sum` wraps
-  `aocore.block_sum` but keeps its own factor-two fast path.
+  centred FFTs keep the `fftshift` convention; `ArrayBackend.centered_coordinates`
+  builds aocore's coordinates directly on the device. `sampling.block_sum` is a
+  thin wrapper over `aocore.block_sum` (no local fast path since aocore 0.1.3,
+  which measured at least as fast at the SH call sites).
+  `sampling.area_rebin` is exact-overlap area averaging between grids of the
+  same extent and any shape ratio, which `aocore.block_sum` (integer factors
+  only) does not cover.
 - `sensors/` contains deterministic ideal optical engines and no camera noise.
   `_shack_hartmann_cuda.py` is a private first-use-JIT execution plan for exact
   compatible CUDA geometries; `shack_hartmann.py` remains the readable physics

CHANGELOG.md

@@ -4,6 +4,65 @@ All notable changes to `makewfs` are documented here.
 
 ## [Unreleased]
 
+## [2.0.0] - 2026-10-07
+
+### Breaking
+
+- **`wfs_input_opd_rms_m` is now the pupil-weighted, piston-removed RMS.**
+  It follows aocore CONVENTIONS.md 4.1: the RMS of the input OPD in metres,
+  weighted by the intensity of the pupil the optics use and with the
+  intensity-weighted mean (piston) removed. Before 2.0 it was the quadratic
+  mean over the whole input grid, with every pixel counted equally, pixels
+  outside the pupil included and piston kept. For the same wavefront the new
+  value is usually smaller; a pure piston now reports 0. The weights are the
+  configured pupil's intensity (the amplitude the engines use, squared) on the
+  input grid: an analytic pupil is evaluated on `input.shape` with the
+  configured `numerics.pupil_supersampling`, the same map
+  `WavefrontSensor.pupil_illumination()` returns; a custom mask, which exists
+  only on the engine's pupil grid, is area-averaged onto the input grid.
+  Phase input is converted to OPD first, and `expose_integrated` reports the
+  RMS of the mean OPD, as before. The value is still reduced on the device and
+  crosses to the host in the same single batched transfer as the captured
+  photon rate. A pupil with no transmission on the input grid is now rejected
+  when the sensor is built.
+- **`makewfs.provenance.metadata`**, an internal helper, takes the two
+  reduced RMS values as required `opd_rms_m` and `opd_rms_unweighted_m`
+  arguments and no longer accepts `opd_m`.
+
+See "Migrating to 2.0" in the stability guide.
+
+### Added
+
+- **`wfs_input_opd_rms_unweighted_m`** frame metadata keeps the 1.x quantity,
+  the unweighted RMS over the whole input grid with piston included, so no
+  information is lost. Read it wherever the old number is still wanted.
+- Both sensor engines expose `configured_pupil`, the pupil amplitude on their
+  own pupil grid, and `makewfs.sampling.area_rebin` area-averages a map onto
+  another grid of the same extent, for any shape ratio.
+
+### Changed
+
+- The `closed_loop_injection.py` and `showcase.py` examples take the residual
+  RMS they plot from each frame's `wfs_input_opd_rms_m` instead of a
+  whole-grid `np.std`, so the numbers they print are pupil-weighted. The
+  checked-in `makewfs_showcase.webp` was not regenerated.
+- **Requires `aocore>=0.1.3,<0.2`.** `makewfs.sampling.block_sum` drops its
+  own factor-two shortcut and delegates every factor to `aocore.block_sum`,
+  whose strided CPU adds and single CuPy kernel measured at least as fast on
+  the spot stacks the Shack-Hartmann actually bins (for example
+  `(400, 16, 16)` float32: CPU 112 vs 173 us, Quadro P620 65 vs 94 us;
+  `(3600, 12, 12)` float64: CPU 0.90 vs 1.27 ms, GPU 98 vs 173 us). The
+  additions happen in a different order, so Shack-Hartmann float32 images
+  differ from 1.2.0 by float32 rounding (at most about 1e-7 relative) and
+  float64 ones by about 2e-16; pyramid images are bit-for-bit unchanged.
+  Pixel-centre coordinates for the input, pupil and DFT detector grids are
+  now built on the selected device by `aocore.centered_coordinates`
+  (`ArrayBackend.centered_coordinates`), with identical values.
+- The input-RMS tests check the unweighted key against
+  `aocore.rms_unweighted`. The per-frame path keeps its own device reduction,
+  because the aocore functions return host floats and would each add a
+  synchronization.
+
 ## [1.2.0] - 2026-10-07
 
 - **Fixed: phase input was reported as metres.** With

CITATION.cff

@@ -5,7 +5,7 @@ type: software
 authors:
   - family-names: Taylor
     given-names: Jacob
-version: 1.2.0
+version: 2.0.0
 date-released: 2026-10-07
 license: MIT
 repository-code: "https://github.com/jacotay7/makewfs"

docs/concepts.md

@@ -31,6 +31,35 @@ in single precision.
 When an input grid does not divide into the configured lenslets, OPD is
 resampled on physical coordinates. Wrapped phase is never interpolated.
 
+## Input wavefront RMS in frame metadata
+
+Every frame records two RMS values of the input wavefront, both in OPD metres
+(phase input is converted at `input.reference_wavelength_m` first). Both
+describe the input array as given, without `input.static_opd_path`; for
+`expose_integrated` they describe the mean of the temporal samples.
+
+| Key | Definition |
+| --- | --- |
+| `wfs_input_opd_rms_m` | Pupil-weighted, piston-removed RMS, the `rms` of aocore CONVENTIONS 4.1: `sqrt(sum a^2 (opd - <opd>_a)^2 / sum a^2)`, where `<opd>_a` is the intensity-weighted mean. |
+| `wfs_input_opd_rms_unweighted_m` | `sqrt(mean(opd^2))` over every pixel of the input grid, with piston included and pixels outside the pupil counted. This was `wfs_input_opd_rms_m` before 2.0. |
+
+The weight `a^2` is the intensity of the pupil the optics use, on the input
+grid. The engines use the configured pupil as a field amplitude, so the
+intensity is its square. An analytic pupil (diameter, obscuration, spiders,
+segment gaps, rotation) is evaluated on `input.shape` with the same
+`numerics.pupil_supersampling`; this is the map
+`WavefrontSensor.pupil_illumination()` returns, computed in float64. A
+`telescope.custom_mask_path` mask exists only on the engine's own pupil grid
+(`lenslets_across_pupil` times the pupil samples per lenslet for a
+Shack-Hartmann, `pixels_across_pupil` for a pyramid), so its intensity is
+area-averaged onto the input grid: each input
+pixel takes the mean over the mask cells it overlaps, which is exact when the
+two shapes match. A piston, or OPD outside the pupil, therefore changes only
+the unweighted value.
+
+Both values are reduced on the selected device and cross to the host together
+with the captured photon rate in one transfer.
+
 ## Closed-loop use
 
 The package intentionally stops at the detector image. A downstream controller

docs/stability.md

@@ -25,3 +25,20 @@ The physical per-frame contract is a wavefront array plus configuration (and an
 optional detector seed). `out` controls storage lifetime only and cannot change
 the simulated result. Atmosphere, reconstruction, controllers, and detector
 physics remain outside the package boundary.
+
+## Migrating to 2.0
+
+makewfs 2.0 changes the meaning of one frame-metadata key; the API and the
+configuration schema are unchanged.
+
+- `wfs_input_opd_rms_m` is now the pupil-weighted, piston-removed RMS of
+  aocore CONVENTIONS 4.1 (see
+  [Concepts](concepts.md#input-wavefront-rms-in-frame-metadata)). For the same
+  wavefront it is usually smaller than before: it no longer includes piston or
+  any OPD outside the pupil.
+- The old quantity, the unweighted RMS over the whole input grid with piston
+  included, is still recorded as `wfs_input_opd_rms_unweighted_m`. Code that
+  needs the 1.x number reads that key instead.
+- `makewfs.provenance.metadata`, an internal helper, now takes both RMS values
+  as required `opd_rms_m` and `opd_rms_unweighted_m` arguments and no longer
+  accepts `opd_m`.

examples/closed_loop_injection.py

@@ -5,7 +5,9 @@
 next *residual* wavefront back to ``WavefrontSensor.expose``. The residual is a
 low-order aberration (defocus + astigmatism) so the Shack-Hartmann spots shift
 differently across the pupil and the relaxation toward the flat reference is
-visible, and a convergence panel tracks the residual RMS.
+visible, and a convergence panel tracks the residual RMS that each frame
+records in its ``wfs_input_opd_rms_m`` metadata: pupil-weighted with piston
+removed, as aocore CONVENTIONS 4.1 defines ``rms``.
 """
 
 from __future__ import annotations
@@ -19,7 +21,7 @@
 
 
 def _low_order_residual(shape: tuple[int, int], amplitude_m: float) -> np.ndarray:
-    """Return a defocus + astigmatism OPD normalized to ``amplitude_m`` RMS."""
+    """Return a defocus + astigmatism OPD with ``amplitude_m`` RMS over the grid."""
     yy, xx = np.mgrid[: shape[0], : shape[1]]
     x = (xx - (shape[1] - 1) / 2) / (shape[1] / 2)
     y = (yy - (shape[0] - 1) / 2) / (shape[0] / 2)
@@ -56,8 +58,9 @@ def main() -> None:
     residual_rms_nm: list[float] = []
     for step in range(args.steps):
         residual = _low_order_residual(shape, amplitude_m=200e-9) * (0.6**step)
-        residual_rms_nm.append(float(np.std(residual) * 1e9))
-        frames.append(np.asarray(sensor.expose(residual, seed=0)).astype(np.float64))
+        frame = sensor.expose(residual, seed=0)
+        residual_rms_nm.append(frame.metadata["wfs_input_opd_rms_m"] * 1e9)
+        frames.append(np.asarray(frame).astype(np.float64))
 
     columns = min(4, args.steps)
     shown = sorted({0, 1, args.steps // 2, args.steps - 1})[:columns]
@@ -91,7 +94,7 @@ def main() -> None:
     convergence.plot(range(args.steps), residual_rms_nm, "o-")
     convergence.set_title("external loop convergence")
     convergence.set_xlabel("loop step")
-    convergence.set_ylabel("residual OPD RMS (nm)")
+    convergence.set_ylabel("residual OPD RMS over the pupil (nm)")
     convergence.grid(True, alpha=0.3)
 
     frame_rms = [float(np.sqrt(np.mean(diff**2))) for diff in diffs]

examples/showcase.py

@@ -263,16 +263,17 @@ def benchmark_panel(panel: Panel, sync: Callable[[], None]) -> float:
 def collect_frames(panel: Panel) -> np.ndarray:
     """``(N_FRAMES, *output_shape)`` detector ADU on the host, plus the OPD rms."""
     import getframes
-    import pyturb
 
     sensor = panel.build_sensor()
     atmosphere = panel.build_atmosphere()
     out = np.empty((N_FRAMES, *panel.output_shape), dtype=np.float32)
     rms = []
     for index in range(N_FRAMES):
         residual = panel.residual(atmosphere, index)
-        rms.append(float(np.std(pyturb.to_numpy(residual))) * 1e9)
         frame = sensor.expose(residual, seed=index)
+        # Pupil-weighted, piston-removed (aocore CONVENTIONS 4.1), reduced on
+        # the sensor's device.
+        rms.append(frame.metadata["wfs_input_opd_rms_m"] * 1e9)
         out[index] = getframes.to_numpy(frame.data).astype(np.float32)
     panel.badge_extra = f"residual {np.mean(rms):.0f} nm rms"
     return np.clip(out - panel.bias_adu, 0.0, None)

pyproject.toml

@@ -25,7 +25,7 @@ classifiers = [
     "Typing :: Typed",
 ]
 dependencies = [
-    "aocore>=0.1.2,<0.2",
+    "aocore>=0.1.3,<0.2",
     "numpy>=1.23",
     "scipy>=1.10",
     "getframes>=2.2.0",

src/makewfs/__about__.py

@@ -1,3 +1,3 @@
 """Package version metadata."""
 
-__version__ = "1.2.0"
+__version__ = "2.0.0"

src/makewfs/api.py

@@ -16,10 +16,11 @@
 from .detector import DetectorAdapter
 from .provenance import metadata as build_metadata
 from .pupil import make_pupil
+from .sampling import area_rebin
 from .sensors.base import OpticalResult, SensorEngine
 from .sensors.pyramid import PyramidEngine
 from .sensors.shack_hartmann import ShackHartmannEngine
-from .wavefront import iter_phase_samples
+from .wavefront import grid_rms, iter_phase_samples, pupil_rms, pupil_weights
 
 
 class WavefrontSensor:
@@ -53,10 +54,12 @@ def __init__(self, config: WFSConfig, *, _backend: ArrayBackend | None = None) -
             launched_rate=0.0,
             captured_rate=0.0,
             opd_rms_m=0.0,
+            opd_rms_unweighted_m=0.0,
             seed=None,
             source_states=self.engine.source_states,
             file_digests=self.engine.file_digests,
         )
+        self._rms_weights = pupil_weights(self._input_pupil_intensity(), backend=self.backend)
 
     @classmethod
     def from_toml(cls, path: str | Path) -> WavefrontSensor:
@@ -66,16 +69,50 @@ def from_toml(cls, path: str | Path) -> WavefrontSensor:
     def _render(self, wavefront: ArrayLike) -> OpticalResult:
         return self.engine.render(cast(NDArray[np.float64], wavefront))
 
-    def _opd_rms(self, opd: Any) -> Any:
-        """Reduce OPD RMS on-device before the batched metadata crossing."""
-        return self.backend.sqrt(self.backend.mean(opd**2))
+    def _input_pupil_intensity(self) -> Any:
+        """Return the pupil intensity the optics use, on the input grid.
+
+        An analytic pupil is evaluated on ``input.shape`` from the same
+        telescope model and ``numerics.pupil_supersampling`` the engine uses on
+        its own grid, as :meth:`pupil_illumination` does; both grids span
+        ``input.grid_extent_m``. A custom mask exists only on the engine's
+        configured pupil grid, so its intensity is area-averaged from there
+        onto the input grid, which is exact when the two shapes match. The
+        engines use the pupil as a field amplitude, so the intensity is its
+        square.
+        """
+        if self.config.telescope.custom_mask_path is None:
+            amplitude = make_pupil(
+                self.config.telescope,
+                self.config.input.shape,
+                self.config.input.grid_extent_m,
+                supersampling=self.config.numerics.pupil_supersampling,
+                backend=self.backend,
+                dtype=np.float64,
+            )
+            return amplitude * amplitude
+        amplitude = self.backend.asarray(self.engine.configured_pupil, dtype=np.float64)
+        return area_rebin(amplitude * amplitude, self.config.input.shape, backend=self.backend)
+
+    def _opd_rms(self, opd: Any) -> tuple[Any, Any]:
+        """Reduce both input-OPD RMS values on the device.
+
+        Returns the pupil-weighted, piston-removed RMS and the unweighted
+        whole-grid RMS as device scalars, for the one batched metadata
+        crossing in ``backend.scalars``.
+        """
+        return (
+            pupil_rms(opd, self._rms_weights, backend=self.backend),
+            grid_rms(opd, backend=self.backend),
+        )
 
     def _frame_metadata(
         self,
         *,
         launched_rate: float,
         captured_rate: float,
         opd_rms_m: float,
+        opd_rms_unweighted_m: float,
         seed: int | None,
     ) -> dict[str, Any]:
         """Copy cached static provenance and fill the per-frame values."""
@@ -85,6 +122,7 @@ def _frame_metadata(
                 "wfs_launched_photons_s": float(launched_rate),
                 "wfs_captured_photons_s": float(captured_rate),
                 "wfs_input_opd_rms_m": float(opd_rms_m),
+                "wfs_input_opd_rms_unweighted_m": float(opd_rms_unweighted_m),
                 "wfs_seed": seed if seed is not None else "internal",
             }
         )
@@ -182,18 +220,26 @@ def expose(
         seed: int | None = None,
         out: Any | None = None,
     ) -> Any:
-        """Render one wavefront into optional caller-owned detector storage."""
+        """Render one wavefront into optional caller-owned detector storage.
+
+        Besides provenance and timings, the frame metadata records the input
+        wavefront's ``wfs_input_opd_rms_m``, its RMS in OPD metres weighted by
+        the pupil intensity with piston removed (aocore CONVENTIONS 4.1), and
+        ``wfs_input_opd_rms_unweighted_m``, its RMS over the whole input grid
+        with piston kept.
+        """
         total_start = perf_counter()
         optical_start = total_start
         result = self._render(wavefront)
-        captured_rate, opd_rms = self.backend.scalars(
-            result.captured_rate_per_s, self._opd_rms(result.opd_m)
+        captured_rate, opd_rms, opd_rms_unweighted = self.backend.scalars(
+            result.captured_rate_per_s, *self._opd_rms(result.opd_m)
         )
         optical_elapsed = perf_counter() - optical_start
         frame_metadata = self._frame_metadata(
             launched_rate=result.launched_rate_per_s,
             captured_rate=captured_rate,
             opd_rms_m=opd_rms,
+            opd_rms_unweighted_m=opd_rms_unweighted,
             seed=seed,
         )
         detector_start = perf_counter()
@@ -244,13 +290,14 @@ def expose_integrated(
             if not samples:
                 raise ValueError("phase_samples must contain at least one sample")
             result = batched(samples)
-            captured_rate, opd_rms = self.backend.scalars(
-                result.captured_rate_per_s, self._opd_rms(result.opd_m)
+            captured_rate, opd_rms, opd_rms_unweighted = self.backend.scalars(
+                result.captured_rate_per_s, *self._opd_rms(result.opd_m)
             )
             frame_metadata = self._frame_metadata(
                 launched_rate=result.launched_rate_per_s,
                 captured_rate=captured_rate,
                 opd_rms_m=opd_rms,
+                opd_rms_unweighted_m=opd_rms_unweighted,
                 seed=seed,
             )
             frame_metadata["wfs_temporal_samples"] = len(samples)
@@ -305,13 +352,14 @@ def expose_integrated(
             None if spectral_rate_sum is None else spectral_rate_sum / sample_count
         )
         average_opd = opd_sum / sample_count
-        captured_rate, opd_rms = self.backend.scalars(
-            captured / sample_count, self._opd_rms(average_opd)
+        captured_rate, opd_rms, opd_rms_unweighted = self.backend.scalars(
+            captured / sample_count, *self._opd_rms(average_opd)
         )
         frame_metadata = self._frame_metadata(
             launched_rate=launched,
             captured_rate=captured_rate,
             opd_rms_m=opd_rms,
+            opd_rms_unweighted_m=opd_rms_unweighted,
             seed=seed,
         )
         frame_metadata["wfs_temporal_samples"] = sample_count