diff --git a/AGENTS.md b/AGENTS.md index 4682772..0414f4a 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -62,6 +62,13 @@ pupil OPD/phase + static config because it is part of WFS image formation, but it never predicts LGS return flux or evolves the sodium layer. +- Stack conventions and generic optics primitives belong to + [`aocore`](https://github.com/jacotay7/aocore), a core dependency. Its + `CONVENTIONS.md` is the contract, and §9 says to import the primitives rather + than re-implement them. Today makewfs takes `centered_coordinates`, + `ARCSEC_TO_RAD`/`RAD_TO_ARCSEC`, `opd_to_phase`/`phase_to_opd` and + `block_sum` from it. Change a convention in aocore, never here. + Do not copy sibling physics for convenience. If their public API is insufficient, write a failing integration test/design note, use the conditional gates in `ROADMAP.md`, and make the smallest change in the owning repository. @@ -141,6 +148,13 @@ Follow the target layout in `ROADMAP.md`: - `config.py` parses and validates; it does not propagate optics. - `wavefront.py`, `pupil.py`, and `sampling.py` hold shared numerical rules. + Some aocore look-alikes stay local on purpose; each says why in its module + docstring. `pupil.make_pupil` builds wedge spiders, square segment gaps, + 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. - `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 diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d0fe23..89351e1 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,24 @@ All notable changes to `makewfs` are documented here. ## [Unreleased] +- **Changed: generic optics helpers now come from + [aocore](https://github.com/jacotay7/aocore).** `aocore>=0.1.2,<0.2` moves + from the `dev` extra to a core dependency. Under CONVENTIONS.md §9, makewfs + now imports these primitives instead of keeping its own copies: + `centered_coordinates` for the pupil, field-stop and DFT detector grids; + `ARCSEC_TO_RAD`/`RAD_TO_ARCSEC` for source angles and the Shack-Hartmann + plate scale; `opd_to_phase`/`phase_to_opd` for phase input and the sensor + phasors; and `block_sum` for pixel integration. No public name changes. + `makewfs.sampling.block_sum` keeps its signature, error messages and + factor-two fast path, and delegates other factors to aocore. The behaviour + already followed the conventions. Float32 renders of the example + configurations are bit-for-bit unchanged, and float64 renders differ only by + floating-point rounding (at most about 6e-16 relative), from the changed + operation order in the angle and phase conversions and, for oversampling + factors above two, in pixel binning. The pupil masks, the + `ArrayBackend` and its `fftshift`-convention centred FFTs stay local + because aocore has no equivalent for them; the module docstrings say why. + - **Tests: conformance with the AO stack conventions.** `tests/test_conformance.py` runs [aocore](https://github.com/jacotay7/aocore)'s checks against both sensors. Shack-Hartmann spots centre between pixels for a flat wavefront and diff --git a/README.md b/README.md index 46aaa56..67dc28f 100644 --- a/README.md +++ b/README.md @@ -22,7 +22,9 @@ quadrature, and sodium-layer LGS geometry — everything from the pupil to the A The package owns the wavefront-sensor optics and nothing else. It deliberately reuses [`pyturb`](https://github.com/jacotay7/pyturb) for atmospheric OPD and [`getframes`](https://github.com/jacotay7/getframes) for detector response and -noise; neither model is reimplemented here. It runs on NumPy by default and +noise; neither model is reimplemented here. Shared conventions and generic optics +primitives (coordinates, angle constants, OPD/phase conversion, pixel binning) +come from [`aocore`](https://github.com/jacotay7/aocore). It runs on NumPy by default and switches to CUDA (via CuPy) with a single configuration field. ## Install diff --git a/docs/units-and-coordinates.md b/docs/units-and-coordinates.md index d65faf8..0f4f96d 100644 --- a/docs/units-and-coordinates.md +++ b/docs/units-and-coordinates.md @@ -1,14 +1,24 @@ # Units and coordinates +`makewfs` follows the AO stack conventions defined by +[`aocore`](https://github.com/jacotay7/aocore) (`CONVENTIONS.md`), and +`tests/test_conformance.py` runs aocore's checks against both sensors. Pixel +centres, arcsecond/radian constants, OPD/phase conversion and flux-conserving +pixel binning come from `aocore` rather than local copies. + - Input arrays are `(y, x)` NumPy arrays. - OPD is metres. Phase is radians and requires `input.reference_wavelength_m`. - `input.grid_extent_m` is the physical width/height represented by the input array; coordinates are centred on pixel centres. -- `source.field_angle_arcsec` is `[x, y]` on-sky angle. +- `source.field_angle_arcsec` is `[x, y]` on-sky angle, converted with + `aocore.ARCSEC_TO_RAD`. - `lgs_launch_position_m` is `[x, y]` in the entrance-pupil plane. - Optical output is photons/s/native detector pixel. Detector output is ADU. -The centred FFT helpers use unitary normalization. Piston removal is a numerical +The centred FFT helpers use unitary normalization. They keep the `fftshift` +convention (zero frequency on pixel `n // 2`) internally; the detector images +the sensors produce follow the stack's centring, with the optical axis at pixel +`(n - 1) / 2`, which the conformance tests check. Piston removal is a numerical stabilization only; it does not alter intensity. Detector cropping reports lost flux rather than renormalizing the image. diff --git a/pyproject.toml b/pyproject.toml index 99afa51..e40e49f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -25,6 +25,7 @@ classifiers = [ "Typing :: Typed", ] dependencies = [ + "aocore>=0.1.2,<0.2", "numpy>=1.23", "scipy>=1.10", "getframes>=2.2.0", @@ -33,7 +34,6 @@ dependencies = [ [project.optional-dependencies] dev = [ - "aocore>=0.1.1", "build>=1.0", "mypy>=1.8", "pytest>=7.0", diff --git a/src/makewfs/api.py b/src/makewfs/api.py index 1ab6a00..9aa05a9 100644 --- a/src/makewfs/api.py +++ b/src/makewfs/api.py @@ -8,6 +8,7 @@ from typing import Any, cast import numpy as np +from aocore import RAD_TO_ARCSEC from numpy.typing import ArrayLike, NDArray from .backend import ArrayBackend, cpu_backend, cupy_backend @@ -20,8 +21,6 @@ from .sensors.shack_hartmann import ShackHartmannEngine from .wavefront import iter_phase_samples -_ARCSEC_PER_RADIAN = 206264.80624709636 - class WavefrontSensor: """Configured wavefront-sensor facade. @@ -125,7 +124,7 @@ def subaperture_plate_scale_arcsec(self) -> float: sampling = engine._spot_sampling(self.config.sensor.wavelength_m) subaperture_m = self.config.input.grid_extent_m / engine.n_lenslets lambda_over_d_rad = self.config.sensor.wavelength_m / subaperture_m - return float(lambda_over_d_rad / sampling * _ARCSEC_PER_RADIAN) + return float(lambda_over_d_rad / sampling * RAD_TO_ARCSEC) def subaperture_field_of_view_arcsec(self) -> float: """Return the on-sky field of view one subaperture's window spans. diff --git a/src/makewfs/backend.py b/src/makewfs/backend.py index 1a62f3c..a7c64e0 100644 --- a/src/makewfs/backend.py +++ b/src/makewfs/backend.py @@ -5,6 +5,14 @@ NumPy directly. That boundary is deliberately private today; it gives a future CuPy implementation one place to provide array creation, reductions, FFT, and interpolation semantics without changing the optical equations. + +It is deliberately separate from ``aocore.Backend``. The sensors choose a dtype +per array rather than one precision per backend, pass explicit FFT worker +counts and ``overwrite_x``, and need SciPy/CuPy ``ndimage`` helpers. The +centred FFTs here use the ``fftshift`` convention (zero frequency on pixel +``n // 2``). The engines own the conversion to CONVENTIONS 1.3 centring: for +example ``sampling.spot_intensity`` applies a half-sample phase ramp for even +detectors, and ``tests/test_conformance.py`` checks both sensors. """ from __future__ import annotations diff --git a/src/makewfs/pupil.py b/src/makewfs/pupil.py index 01acb6a..57e05fb 100644 --- a/src/makewfs/pupil.py +++ b/src/makewfs/pupil.py @@ -1,4 +1,12 @@ -"""Analytic and custom entrance-pupil masks.""" +"""Analytic and custom entrance-pupil masks. + +These stay local rather than using ``aocore.Pupil``: the configured geometry +has no aocore equivalent. Spiders are angular wedges of ``width_fraction`` of +pi, segment gaps lie on a square grid, the whole pattern rotates with +``pupil_rotation_deg``, grids may be rectangular, and masks are built on the +selected array backend in the working precision. Pixel centres follow +CONVENTIONS 1.2 through ``wavefront._coordinates``. +""" from __future__ import annotations diff --git a/src/makewfs/sampling.py b/src/makewfs/sampling.py index 2d50914..2245516 100644 --- a/src/makewfs/sampling.py +++ b/src/makewfs/sampling.py @@ -8,6 +8,8 @@ from typing import Any, cast import numpy as np +from aocore import block_sum as _aocore_block_sum +from aocore import centered_coordinates from numpy.typing import NDArray from .backend import ArrayBackend, centered_fft_intensity, cpu_backend @@ -142,7 +144,13 @@ def crop_center(array: NDArray[Any], shape: tuple[int, int]) -> NDArray[Any]: def block_sum( array: NDArray[Any], factor: int, *, backend: ArrayBackend | None = None ) -> NDArray[Any]: - """Sum square pixel blocks while preserving total flux.""" + """Sum square pixel blocks while preserving total flux. + + A thin wrapper over ``aocore.block_sum``, which owns flux-conserving + binning (CONVENTIONS 9). It keeps this module's error messages and a + factor-two fast path. ``backend`` is accepted for compatibility; the + reduction runs on the array's own namespace. + """ if factor < 1: raise ValueError("factor must be positive") if factor == 1: @@ -161,10 +169,7 @@ def block_sum( + array[..., 1::2, 0::2] + array[..., 1::2, 1::2], ) - resolved = backend or cpu_backend() - reshaped = array.reshape((*array.shape[:-2], height // factor, factor, width // factor, factor)) - reduced_x = resolved.sum(reshaped, axis=-1) - return cast(NDArray[Any], resolved.sum(reduced_x, axis=-2)) + return cast(NDArray[Any], _aocore_block_sum(array, factor)) @dataclass @@ -245,9 +250,8 @@ def build( coordinate = backend.arange(nfft, dtype=np.float64) plan.half_sample = backend.exp(-1j * math.pi * coordinate / nfft) else: - detector_coordinate = ( - backend.arange(high_resolution_pixels, dtype=np.float64) - - (high_resolution_pixels - 1) / 2.0 + detector_coordinate = backend.asarray( + centered_coordinates(high_resolution_pixels), dtype=np.float64 ) / (sampling * oversampling) pupil_coordinate = backend.arange(samples_per_lenslet, dtype=np.float64) kernel = backend.exp( @@ -259,12 +263,11 @@ def build( ) plan.dft_kernel = backend.astype(kernel, plan.field_dtype) if field_stop_radius_lambda_over_d is not None: - coordinates = backend.arange(high_resolution_pixels, dtype=np.float64) + coordinates = backend.asarray( + centered_coordinates(high_resolution_pixels), dtype=np.float64 + ) y, x = backend.meshgrid(coordinates, coordinates, indexing="ij") - radius_lambda_over_d = backend.hypot( - x - (high_resolution_pixels - 1) / 2.0, - y - (high_resolution_pixels - 1) / 2.0, - ) / (oversampling * sampling) + radius_lambda_over_d = backend.hypot(x, y) / (oversampling * sampling) plan.field_stop_mask = radius_lambda_over_d <= field_stop_radius_lambda_over_d return plan diff --git a/src/makewfs/sensors/pyramid.py b/src/makewfs/sensors/pyramid.py index 167e772..05d951a 100644 --- a/src/makewfs/sensors/pyramid.py +++ b/src/makewfs/sensors/pyramid.py @@ -6,6 +6,7 @@ from typing import Any, cast import numpy as np +from aocore import opd_to_phase from numpy.typing import NDArray from ..backend import ArrayBackend, centered_fft2, centered_ifft2, complex_dtype, real_dtype @@ -153,7 +154,7 @@ def _base_phasor( total_opd = internal + self._field_angle_opd[state_index] piston = total_opd[self._piston_index] relative_opd = total_opd - piston - phase = 2.0 * math.pi * relative_opd / state.wavelength_m + phase = opd_to_phase(relative_opd, state.wavelength_m) return cast( NDArray[Any], self.backend.asarray(self.backend.exp(1j * phase), dtype=self._complex_dtype), diff --git a/src/makewfs/sensors/shack_hartmann.py b/src/makewfs/sensors/shack_hartmann.py index 5366ea5..b6ea4e0 100644 --- a/src/makewfs/sensors/shack_hartmann.py +++ b/src/makewfs/sensors/shack_hartmann.py @@ -7,6 +7,7 @@ from typing import Any, cast import numpy as np +from aocore import opd_to_phase from numpy.typing import NDArray from ..backend import ArrayBackend, complex_dtype, real_dtype @@ -379,7 +380,7 @@ def _field( # it makes the exact physical piston invariance survive finite precision # in exp(i*phase) and has no effect on the intensity. relative_opd = total_opd - piston - phase = 2.0 * math.pi * relative_opd / state.wavelength_m + phase = opd_to_phase(relative_opd, state.wavelength_m) return cast( NDArray[Any], self.backend.asarray( @@ -416,7 +417,7 @@ def _fields( [self.source_states[index].wavelength_m for index in state_group] * samples, dtype=self._real_dtype, )[:, None, None] - phase = 2.0 * math.pi * (total_opd - piston) / wavelengths + phase = opd_to_phase(total_opd - piston, wavelengths) return cast( NDArray[Any], self.backend.asarray( diff --git a/src/makewfs/source.py b/src/makewfs/source.py index 3ee527b..92f7d14 100644 --- a/src/makewfs/source.py +++ b/src/makewfs/source.py @@ -7,6 +7,7 @@ from pathlib import Path import numpy as np +from aocore import ARCSEC_TO_RAD from numpy.typing import NDArray from .config import WFSConfig @@ -33,8 +34,8 @@ def _normalised(values: tuple[float, ...], count: int) -> NDArray[np.float64]: def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]: """Return ``(x angle, y angle, weight)`` Gaussian angular quadrature.""" source = config.source - base_x = math.radians(source.field_angle_arcsec[0] / 3600.0) - base_y = math.radians(source.field_angle_arcsec[1] / 3600.0) + base_x = source.field_angle_arcsec[0] * ARCSEC_TO_RAD + base_y = source.field_angle_arcsec[1] * ARCSEC_TO_RAD if source.angular_kernel_path is not None: data = np.atleast_2d( np.asarray(np.loadtxt(Path(source.angular_kernel_path)), dtype=np.float64) @@ -48,8 +49,8 @@ def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]: weights = data[:, 2] / np.sum(data[:, 2]) return [ ( - base_x + math.radians(float(x) / 3600.0), - base_y + math.radians(float(y) / 3600.0), + base_x + float(x) * ARCSEC_TO_RAD, + base_y + float(y) * ARCSEC_TO_RAD, float(weight), ) for (x, y), weight in zip(data[:, :2], weights) @@ -60,7 +61,7 @@ def _angular_states(config: WFSConfig) -> list[tuple[float, float, float]]: # Three-point Gauss-Hermite is exact through fourth order for a normal # expectation. For higher configured orders use an evenly weighted grid; # this keeps the rule deterministic without requiring an optional package. - sigma = math.radians(source.angular_fwhm_arcsec / 3600.0) / 2.3548200450309493 + sigma = source.angular_fwhm_arcsec * ARCSEC_TO_RAD / 2.3548200450309493 if order == 3: nodes = np.array([-math.sqrt(3.0), 0.0, math.sqrt(3.0)]) one_d = np.array([1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0]) diff --git a/src/makewfs/wavefront.py b/src/makewfs/wavefront.py index ae4d2a8..8dc6ae7 100644 --- a/src/makewfs/wavefront.py +++ b/src/makewfs/wavefront.py @@ -7,6 +7,7 @@ from typing import Any, cast import numpy as np +from aocore import centered_coordinates, phase_to_opd from numpy.typing import ArrayLike, NDArray from .backend import ArrayBackend, cpu_backend @@ -20,11 +21,16 @@ def _coordinates( backend: ArrayBackend | None = None, dtype: object = np.float64, ) -> tuple[NDArray[np.float64], NDArray[np.float64]]: - """Return centered physical ``(x, y)`` coordinates for an array shape.""" + """Return centered physical ``(x, y)`` coordinates for an array shape. + + Pixel centres follow ``aocore.centered_coordinates`` (CONVENTIONS 1.2). The + unit-pitch offsets are exact in either precision; scaling to metres happens + on the backend in ``dtype``. + """ resolved = backend or cpu_backend() height, width = shape - x = (resolved.arange(width, dtype=dtype) - (width - 1) / 2.0) * extent_m / width - y = (resolved.arange(height, dtype=dtype) - (height - 1) / 2.0) * extent_m / height + x = resolved.asarray(centered_coordinates(width), dtype=dtype) * extent_m / width + y = resolved.asarray(centered_coordinates(height), dtype=dtype) * extent_m / height xx, yy = resolved.meshgrid(x, y) return xx, yy @@ -91,7 +97,7 @@ def opd( converted = self.backend.asarray(array, dtype=np.float64) if self.config.input.quantity == "phase": assert self.config.input.reference_wavelength_m is not None - converted = converted * self.config.input.reference_wavelength_m / (2.0 * np.pi) + converted = phase_to_opd(converted, self.config.input.reference_wavelength_m) if self.static_opd is not None: converted = converted + self.static_opd if not self.backend.scalar(self.backend.all(self.backend.isfinite(converted))): diff --git a/tests/test_conformance.py b/tests/test_conformance.py index b5b4e4a..19a8045 100644 --- a/tests/test_conformance.py +++ b/tests/test_conformance.py @@ -3,11 +3,10 @@ from __future__ import annotations import numpy as np -import pytest +from aocore import ARCSEC_TO_RAD, conformance from makewfs import WavefrontSensor, load_config -conformance = pytest.importorskip("aocore.conformance") SH_CONFIG = "examples/configs/shack_hartmann_minimal.toml" PYRAMID_CONFIG = "examples/configs/pyramid_minimal.toml" @@ -19,8 +18,6 @@ def _sensor(path: str): def test_shack_hartmann_spots_centre_between_pixels_and_follow_tilt() -> None: - from aocore import ARCSEC_TO_RAD - config, sensor, image = _sensor(SH_CONFIG) shape = config.input.shape conformance.check_image_centring(image, pupil_shape=shape)