Skip to content

Changelog

All notable changes to solvephase are documented here. The project follows Semantic Versioning.

[Unreleased]

[0.2.0] - 2026-10-07

Performance

Small GPU problems were bound by host-side launch overhead (15-30 us per CuPy operation on an Arm host), not by the GPU. Results are unchanged: bitwise identical to 0.1.0 on CPU and GPU on the reference machine (Ampere Neoverse-N1, 12 pinned cores, NumPy 2.5.3/SciPy 1.18.1; RTX 4060 with CuPy 14.2.0). Speed-ups below are medians of 0.1.0 and 0.2.0 run alternately, three times each (the Benchmarks page, docs/benchmarks.md, has the tables).

  • Fast & Furious on the GPU replays each step from a CUDA graph: one copy in, one graph launch and one copy out per frame, running the same kernels. A step takes 0.11 ms at 64² (was 1.41 ms, 12x) and 0.27 ms at 256² (5.8x), so the GPU now beats the CPU on every size.
  • Fused GPU kernels (solvephase._kernels) for the focal-plane model, its reverse- and forward-mode derivatives, the Poisson, Gaussian and amplitude losses, the Gerchberg-Saxton projection, the Fast & Furious step and the phase-diversity metric. On the RTX 4060: objective + gradient 1.5-1.8x, Levenberg-Marquardt 1.4-1.7x, LIFT 1.8x, phase diversity 1.4-1.5x, retrieve() 1.2x, Gerchberg-Saxton 1.1-1.2x.
  • L-BFGS reuses the scalar products it already has, fuses its vector updates and evaluates real scalar products on the GPU with less overhead: coded-diffraction L-BFGS 1.2-1.3x on the GPU.
  • FocalPlaneProblem keeps device copies of its flux, background and scaling vectors (no host-to-device copy per evaluation), caches the Gauss-Newton direction maps and assembles gradients and Jacobian blocks without per-channel loops.
  • CPU: unit phasors are built from cos/sin instead of a complex exp (same values, about 1.4x faster), Gerchberg-Saxton computes its diversity phasors once instead of every iteration, and the Poisson and amplitude losses skip their Taylor terms when no pixel is below the floor. Gerchberg-Saxton 1.2x, focal-plane objective and LM up to 1.1x on 12 Neoverse-N1 cores.

Fixed

  • The CDI CPU/GPU parity test failed on aarch64. After 50 iterations it compared the objects to atol=1e-9, but RAAR and DM are chaotic and amplify the FFT libraries' last-bit differences (about 1e5x between iterations 10 and 40 on an Arm Neoverse-N1 with an RTX 4060). The test now compares the objects after a short run (atol=1e-10) and the error histories after the long one (rtol=1e-6). The implementations already agreed.

Added

  • Validation on real Keck/NIRC2 data (validation/nirc2.py). Daytime focus-diversity calibrations of the Keck AO bench, with known coma and trefoil patterns injected on the Xinetics DM, are retrieved with FocalPlaneProblem.
  • The injected patterns are recovered with correlations of 0.945–0.965. One DM gain fits all six runs (3 % scatter), and strength ratios come out within 1 % of the commanded values.
  • The IDL sharpening correction is predicted with no free parameter (correlation 0.90, amplitude within 1 %).
  • The fitted DM gain (652 nm OPD/V) and pupil size (9.75 actuator pitches) agree with the Keck AO software's Xinetics parameters (600 nm/V; a 10.0-pitch control aperture).
  • The data are not distributed; the script reads them from --data or SOLVEPHASE_NIRC2_DATA. The report and figures are in validation/artifacts.

[0.1.0] - 2026-10-07

First release.

Focal-plane phase retrieval

  • FocalPlaneModel: pupil-to-detector image formation for any number of diversity channels.
  • Broadband light (exact per-wavelength sampling).
  • Any detector sampling, including undersampled, with pixel integration (oversample).
  • Reverse-mode (vjp) and forward-mode (jvp) derivatives, verified against finite differences.
  • Propagators with exact adjoints:
  • zero-padded FFT;
  • matrix Fourier transform, stacked over wavelengths;
  • angular spectrum (Fresnel).

An automatic engine choice picks FFT or MFT per model. The model matches HCIPy to machine precision. - FocalPlaneProblem + solve: maximum-likelihood retrieval. - Losses: Poisson deviance with read noise, weighted Gaussian, or the amplitude metric. - Nuisance fits: per-image flux, background and registration tilt, and an optional zonal pupil amplitude. - Regularization: Tikhonov, a KL-turbulence prior (MAP), and zonal smoothness. - Optimizers: Levenberg-Marquardt with an exact Gauss-Newton/Fisher matrix, L-BFGS with a strong-Wolfe line search, or Adam. All are device-resident. - Single precision reaches the double-precision noise floor, because likelihoods and normal equations accumulate in float64. - gerchberg_saxton: Gerchberg-Saxton and the multi-plane Misell algorithm, with exact modulus projection on the full FFT grid. - unwrap_phase: weighted least-squares unwrapping over obscured, spidered and segmented pupils (DCT-preconditioned CG), with 2π bridging of disconnected regions. - retrieve: one-call retrieval with a robust default strategy. - Amplitude-metric Levenberg-Marquardt from a flat start and from Gerchberg-Saxton; the better fit is kept. - Then a Poisson maximum-likelihood polish and optional zonal refinement.

Algorithm families

  • Phase diversity with an unknown extended object (phase_diversity, PhaseDiversityProblem).
  • Gonsalves/Paxman reduced metric with an exact gradient.
  • Apodization and a diffraction-cutoff frequency mask.
  • Automatic noise and regularization estimates.
  • Wiener object estimate.
  • LIFT (lift, LIFT, lift_crlb): maximum-likelihood low-order sensing from one astigmatic image, with a twin check and Cramér-Rao bounds.
  • Fast & Furious (FastAndFurious, simulate_closed_loop): sequential weak-phase focal-plane sensing with the Korkiakoski et al. (2014) modifications. It takes one FFT pair per step.
  • Coherent diffraction imaging (cdi).
  • Algorithms: ER, HIO, DM, RAAR, RRR, ASR/Douglas-Rachford, HPR and OSS.
  • Constraints: real, positive, amplitude bounds; beamstop masks.
  • Chainable schedules, shrinkwrap, and batched multi-start (one batched FFT per iteration).
  • Helpers: autocorrelation_support, align_object (ambiguity-aware) and simulate_cdi.
  • Generic measurement models (solvephase.algorithms.wirtinger, solvephase.operators).
  • Algorithms: Wirtinger flow, truncated Wirtinger flow, truncated and reweighted amplitude flows, and an L-BFGS solver.
  • Initializations: spectral, truncated, orthogonality-promoting and optimal-preprocessing spectral.
  • Operators: Gaussian matrix, coded diffraction and oversampled Fourier, each with an exact adjoint.
  • Transport of intensity (tie, simulate_defocus_stack).
  • Solvers: FFT and DCT Poisson, uniform intensity and Teague's non-uniform solution, and an exact masked conjugate-gradient solver.
  • Axial derivative: central differences or a multi-plane polynomial fit.

Optics and integration

  • Pupil: anti-aliased circular, annular and spidered pupils; segmented hexagonal apertures with segment labels; Keck, JWST, VLT and Hubble presets.
  • Basis:
  • Zernike (annular), KL (with prior variances), Fourier and any other aobasis generator;
  • segment piston/tip/tilt;
  • Gaussian and measured DM influence functions;
  • zonal.
  • solvephase.interop: pyturb turbulence (turbulence_opd) and getframes detector frames (expose). expose converts ADU to electrons and masks saturated pixels.
  • Ambiguity-aware metrics: rms, wavefront_error (with twin handling) and strehl_from_rms.
  • A solvephase command line with info and retrieve.
  • CPU threading tuned for contended and hyper-threaded machines:
  • FFT worker counts scale with the transform size.
  • BLAS threads are limited, per call, for matrix Fourier transforms and Gauss-Newton products (via threadpoolctl). This made broadband gradients up to 10x faster.
  • Scalar products avoid threaded level-1 BLAS.

Comparison

  • benchmarks/methods.py runs every method on one reference problem. All focal-plane methods see the same wavefront. It reports:
  • the data each method needs;
  • accuracy at a standard SNR;
  • warm CPU and GPU time;
  • capture range.

Its output feeds the comparison tables on the "Choosing an algorithm" docs page, which also has a decision flowchart. - The README showcase races every focal-plane method on one clock, with a live error-versus-time chart, plus a gallery of CDI, coded diffraction and TIE. - The Keck preset now follows the documented geometry: 10.95 m across the corners, pointy-top segments, six 26 mm support arms, and a 2.57 m central obscuration.

Shared core

  • The backend, pupils, propagators, wavefront metrics and phase unwrapping now live in aocore, the AO stack's shared core, and solvephase re-exports them unchanged.
  • tests/test_conformance.py runs aocore's convention checks against solvephase: image centring, tilt direction, unit flux and the RMS definition.

Quality

  • More than 280 tests: adjoints, gradients, recovery accuracy, Cramér-Rao efficiency, and CPU/GPU parity.
  • A validation suite against theory and HCIPy.
  • Benchmark and comparison suites with versioned artifacts.
  • An mkdocs documentation site whose code examples are executed in CI.