Skip to content

Performance and precision

Construct a WavefrontSensor once. Pupil masks, source quadrature, FFT geometry, pyramid masks, and the getframes.Camera adapter are built at construction; photon_rate and expose are the warm paths.

The CPU implementation uses SciPy batched FFTs and supports fft_workers plus float32/complex64 or float64/complex128 optical arithmetic. Float32 is useful for throughput and memory; float64 is the reference path. Benchmark cold construction separately from warm frames and record Python, NumPy, SciPy, and hardware versions.

The array-reference Shack--Hartmann path integrates the common two-times oversampling case with direct strided sums. Temporally integrated exposures batch their fields and avoid a caller-side optical render for every sample. These are mathematically equivalent allocation/reduction optimizations.

First-use-JIT Shack--Hartmann execution

Compatible sampled-DFT Shack--Hartmann configurations on CUDA automatically use a shape-specialized compiled executor. The persistent sensor generates a CUDA C kernel for its precision, temporal sample count, lenslet geometry, detector sampling, field stop, margin, wavelength count, and detector-owned charge diffusion footprint. CuPy compiles it on the first launch and caches the binary both in the process and on disk. A later process can therefore reuse the disk cache, but a new geometry or isolated cache pays the compilation cost again.

The kernel fuses phase formation, separable DFT propagation, field-stop masking, temporal intensity averaging, pixel integration, mosaic assembly, and spectral and total photon-rate accumulation; the executor then performs one captured-flux reduction. When charge diffusion is configured, its focal-plane convolution and native-pixel integration are composed once at plan construction into one exact linear operator. This retains the zero-padded focal boundary while removing the repeated convolution traversal from every lenslet, wavelength, and temporal sample. Wavelength/source states still launch and accumulate in their declared order, so incoherent intensity and spectral-QE inputs are preserved.

The readable array implementation remains the reference and automatic fallback for CPU execution, FFT-resolved spot geometries, continuous optical Gaussian blur, measured native-pixel optical kernels, and geometries exceeding portable CUDA block/shared-memory limits. Compiled renders allocate independent public rate outputs on every call; cached plan storage is immutable and never aliases a returned result.

An isolated-cache Quadro P620 benchmark of the physically complete HAKA example (57x57 lenslets, 4x4 detector pixels, four temporal samples, eight wavelengths, 4x focal sampling, and the measured 9x9 OCAM2K charge-diffusion kernel) measured 24.85 ms compiled versus 569.67 ms reference p50: 22.93x faster and a 95.64% latency reduction. The cold first use took 3.176 s. Photon-rate, spectral-rate, and captured-flux disagreements were respectively 4.85e-8, 9.74e-8, and 4.71e-8 relative. Reproduce the matched alternating-order measurement with:

CUPY_CACHE_DIR=/tmp/makewfs-isolated-cache \
  python benchmarks/benchmark_compiled_sh_executor.py \
  --haka --frames 20 --temporal-samples 4

The versioned record is benchmarks/haka-compiled-sh-executor-quadro-p620.json. It is an owner-isolated hardware measurement, not an end-to-end or portable latency guarantee.

In the array fallback on GPU, Shack--Hartmann source states that resolve to the same FFT size are propagated in one device batch. State accumulation order and the CPU reference remain unchanged. A configured wavelength-scaled field stop keeps states sequential because its mask differs with sampling. On a Quadro P620, the eight-state, four-temporal-sample HAKA-class path groups the first five wavelengths and leaves the other three single. A matched, alternating-order 64-frame benchmark measured 37.19 ms sequential versus 36.39 ms grouped median optics time, a 2.16% reduction. The maximum relative ideal photon-rate difference was 5.1e-8. Reproduce it with:

python benchmarks/benchmark_sh_state_batching.py path/to/haka-resolved-wfs.toml \
  --frames 64 --temporal-samples 4

The versioned Quadro P620 record is benchmarks/haka-sh-state-batching-quadro-p620.json. It is a local matched measurement, not a portable latency guarantee.

Install the optional CUDA extra, then set numerics.device = "gpu" in the WFS TOML:

python -m pip install 'makewfs[gpu]'
python -m pytest -q -m gpu

Sensor allocations, reductions, FFTs, interpolation, blur hooks, photon-rate maps, detector stochastic samples, truth, and ADU remain on the device. frame.data is the zero-copy CuPy interface; np.asarray(frame), FITS output, plotting, or getframes.to_numpy intentionally transfer to the host. CPU and GPU use independent random streams, so optical arrays are compared numerically while detector parity is validated statistically and seeded repetition is checked on each backend.

CPU versus GPU reference throughput

The August 2026 development reference used an AMD Ryzen 9 9950X3D and NVIDIA RTX 5090 with Python 3.12, NumPy 2.2.6, SciPy 1.16.3, CuPy 14.1.1, getframes 2.1.1, and pyturb 1.0.0. Each row constructs one persistent sensor, performs an untimed end-to-end warm-up, and records 100 frames. A zero-OPD array enters on the selected device; optics, detector stochastic samples, truth, and ADU remain there. Every frame receives a distinct deterministic detector-noise seed. CUDA is synchronized around timed regions; construction and host transfers are excluded from frames/s.

“Work samples” is source states multiplied by modulation points: it is 9 for the three-wavelength by three-range LGS case and 8 or 32 for the modulated pyramid cases.

`shack_hartmann_quadrature_9sample.toml` crosses three wavelengths with three beacon
ranges to reach nine work samples. That spectral axis is a deliberate quadrature
*load*, not sodium physics: a sodium beacon returns the 589 nm D2 line broadened
by roughly 0.003 nm, some three orders of magnitude narrower than the config's
585–595 nm sweep. Spot elongation comes from the range spread and the off-axis
launch, not from colour. Use `examples/showcase.py`, which samples a
monochromatic beacon across the sodium layer's depth, as the physical reference.
Configuration Sensor Output Work samples CPU (frames/s) GPU (frames/s) Speedup
shack_hartmann_20x20_float32.toml SH 160x160 1 126.5 2,041.8 16.13x
shack_hartmann_60x60_float64.toml SH 360x360 1 17.9 899.3 50.11x
shack_hartmann_quadrature_9sample.toml SH 64x64 9 186.0 780.4 4.20x
pyramid_40_float32.toml pyramid 54x54 1 3,526.4 1,604.7 0.46x
pyramid_60_mod8_float32.toml pyramid 80x80 8 659.0 1,645.9 2.50x
pyramid_80_mod32_float64.toml pyramid 108x108 32 24.3 923.6 37.94x

Higher frames/s is better. Static source/range geometry, modulation phasors, interpolation grids, normalizers, and monochromatic spectral views are cached. SH uses an intensity-only centered FFT that removes an irrelevant input permutation; both sensors use native orthonormal FFT scaling, a fixed illuminated piston reference, and batched metadata scalar transfers. The GPU detector improvements described in getframes are included in every end-to-end row.

Persistent Shack--Hartmann engines also cache half-sample ramps, sampled-DFT kernels, field-stop masks, and backend blur kernels. This targets geometry-heavy field-stop and arbitrary-sampling cases; the dominant plain large-FFT path keeps the same transform workload and should not be expected to change materially.

Pyramid behavior still shows the GPU crossover clearly: launch overhead makes the tiny unmodulated case slower, eight modulation points provide a 2.50x gain, and the 32-point float64 case reaches 37.9x. The nine-sample quadrature case has nine incoherent source states but a small 64x64 detector, so it gains only 4.20x.

The rendered snapshot and raw JSON record exact configuration paths, command, revision, dirty-checkout flag, environment, timings, and detector-only rates. The checked-in snapshot records a specific local environment and is not a cross-hardware performance guarantee.

The repository benchmark runner separates cold construction, warm optical frames, end-to-end frames, and detector-only frames:

python benchmarks/run.py --frames 10 --output benchmark-results.json
python benchmarks/run.py --device both --frames 100 \
  --config benchmarks/configs/shack_hartmann_20x20_float32.toml \
  --config benchmarks/configs/shack_hartmann_60x60_float64.toml \
  --config benchmarks/configs/shack_hartmann_quadrature_9sample.toml \
  --config benchmarks/configs/pyramid_40_float32.toml \
  --config benchmarks/configs/pyramid_60_mod8_float32.toml \
  --config benchmarks/configs/pyramid_80_mod32_float64.toml \
  --output benchmarks/device-results.json
python benchmarks/render_device_table.py benchmarks/device-results.json \
  --output benchmarks/device-results.md

Benchmark JSON includes source-state count, output shape, exact invocation, revision state, environment, methodology, and per-frame optical, end-to-end, and detector-only timings. CUDA is synchronized around timed regions.

Each returned detector frame also records wfs_optical_render_s, wfs_detector_expose_s, and wfs_total_expose_s. These are lightweight diagnostics for a closed-loop driver, not scheduling guarantees or controller abstractions.

Sequential pipeline owners may provide a backend-native uint32 out= array to WavefrontSensor.expose() or expose_integrated(). The detector adapter uses a reusable getframes.DetectorWorkspace only for this explicit path, preserving independent frame lifetimes for ordinary calls. The returned Frame.data is the exact destination and remains valid only until the caller reuses it. On the local i7-10700, the physical OCAM2K 240×240 detector / 228×228 ROI case improved 1.205× and cut traced peak allocation by 49.5%. Isolated Quadro P620 detector timing and a three-pair full AO-loop trial were throughput-neutral, so the supported GPU graph does not enable the destination automatically.

Representative Shack–Hartmann configurations are provided under benchmarks/configs/:

python benchmarks/run.py --representative --frames 3 --output benchmark-results.json

They cover 20×20 float32 and 60×60 float64 lenslet grids, a nine-sample spectrally and range-resolved SH source, and 40/60/80-pixel pyramid cases with 1/8/32 modulation samples while keeping the minimal SH and pyramid smoke cases in the default benchmark run.

The repository also keeps a dated reference snapshot with hardware, Python, dependency, precision, source-state, construction, warm-optics, and detector columns in benchmarks/reference-table.md. Regenerate it on a new machine with:

python benchmarks/run.py --representative --frames 3 --measure-memory \
  --output benchmarks/reference-results.json
python benchmarks/render_table.py benchmarks/reference-results.json --output benchmarks/reference-table.md

For call-level hotspot inspection, write a compact cumulative cProfile report:

python benchmarks/profile_warm.py --frames 3 --output profile-results.json

CI also runs benchmarks/check_regression.py on the representative report. It checks generous same-run ratios for 60x60 versus 20x20 SH, the nine-state quadrature LGS case versus minimal SH, and 32-sample versus unmodulated pyramid optics. These are order-of-magnitude regression alarms, not portable latency promises; inspect the full JSON report for machine-specific performance work.