API reference

High level

solvephase.retrieve

retrieve(images, pupil, wavelength, *, sampling=None, pixel_scale=None, diversity=None, basis=36, zonal_refinement=False, method='auto', loss='poisson', read_noise=0.0, weights=None, fit_background=True, fit_tilt=False, spectral_weights=None, oversample=1, diameter=None, start=None, max_iter=None, device='cpu', precision=None)

Retrieve a pupil wavefront from focal-plane images in one call.

Parameters:

Name	Type	Description	Default
images	Any	(K, my, mx) images (or one (my, mx) image), ideally in photo-electrons, background included.	required
pupil	Pupil | int	A :class:~solvephase.Pupil, or an int n for an n x n circular pupil of diameter (default 1 m).	required
wavelength	float | Sequence[float]	Wavelength in metres, or an array sampling the band.	required
sampling	float | None	Detector sampling in pixels per lambda/D, or pixel scale in radians (give one).	None
pixel_scale	float | None	Detector sampling in pixels per lambda/D, or pixel scale in radians (give one).	None
diversity	Any	Known per-image aberrations: (K, ny, nx) OPD maps in metres; a list of K numbers, read as defocus (Noll 4) RMS in metres per image; or {"noll": j, "amounts": [...]}. Required for more than one image.	None
basis	Basis | int | str | None	Wavefront parameterization: number of Zernike modes (default 36, tip/tilt upward), a :class:~solvephase.Basis, or "zonal".	36
zonal_refinement	bool	After the modal solve, refine pixel by pixel (L-BFGS).	False
method	str	"auto" (robust multi-start: Levenberg-Marquardt from a flat start and from Gerchberg-Saxton, best kept, then a maximum-likelihood polish), "lm", "lbfgs", "gs" (Gerchberg-Saxton/Misell only) or "gs+lm".	'auto'
loss	str | Loss	Final data term: "poisson" (default; images in photo-electrons), "gaussian" or "amplitude".	'poisson'
read_noise	float	Read noise in electrons (Poisson loss).	0.0
weights	Any	Per-pixel weights; zero masks bad pixels.	None
fit_background	bool	Fit a constant background per image.	True
fit_tilt	bool	Fit a registration tip/tilt per image relative to the first.	False
spectral_weights	Sequence[float] | None	Relative photon flux per wavelength for broadband data.	None
oversample	int	Model pixels per detector pixel per axis (pixel integration).	1
diameter	float | None	Diameter in metres when pupil is an int.	None
start	Any	Initial guess: OPD map, coefficients or a previous :class:Result.	None
max_iter	int | None	Iteration limit of the final stage.	None
device	BackendLike	"cpu", "gpu" or "auto"; "single" or "double".	'cpu'
precision	BackendLike	"cpu", "gpu" or "auto"; "single" or "double".	'cpu'

Returns:

Type	Description
Result	result.opd (metres), result.coefficients, model images, flux/background, convergence history. result.extra["stages"] records each stage's loss and time.

Examples:

>>> import solvephase as sp
>>> pupil = sp.Pupil.circular(64, 8.0, obscuration=0.14)
>>> result = sp.retrieve(images, pupil, 1.6e-6, sampling=2.0,
...                      diversity=[0.0, 0.4e-6])
>>> result.opd, result.coefficients

solvephase.Result dataclass

Result(method, opd, phase, amplitude, pupil, wavelength, coefficients=None, basis=None, model_images=None, flux=None, background=None, tilts=None, loss=math.nan, history=list(), times=list(), n_iter=0, converged=False, message='', elapsed=0.0, device='cpu', extra=dict())

Outcome of a phase retrieval.

Array fields are backend arrays (CuPy for GPU solves) until you call :meth:to_numpy. Unknown or not-applicable fields are None.

Attributes:

Name	Type	Description
method	str	Algorithm name.
opd	Any	(ny, nx) retrieved OPD in metres, zero outside the pupil.
phase	Any	The same wavefront as phase in radians at :attr:wavelength.
amplitude	Any	(ny, nx) pupil amplitude (the known pupil unless it was fitted).
pupil	Pupil | None	The pupil the wavefront is defined on.
wavelength	float	Reference wavelength of :attr:phase, in metres.
coefficients	ndarray | None	Modal coefficients (RMS OPD in metres per unit-RMS mode) when a modal basis was used.
basis	Basis | None	The basis :attr:coefficients refer to.
model_images	Any	(K, my, mx) model of the data at the solution, in data units.
flux, background		Per-channel fitted flux (data units, full PSF) and background (data units per pixel).
tilts	ndarray | None	(K, 2) per-channel (tip, tilt) OPD RMS in metres added to the common wavefront (registration), when fitted.
loss	float	Final objective value.
history	list[float]	Objective value per iteration.
times	list[float]	Seconds since the start per :attr:history entry.
n_iter	int	Iterations performed.
converged	bool	Whether a convergence criterion stopped the solver.
message	str	Why it stopped.
elapsed	float	Total wall-clock seconds, including setup.
device	str	"cpu" or "gpu".
extra	dict[str, Any]	Algorithm-specific outputs (object estimate, support, error metrics ...).

to_numpy

to_numpy()

Copy with every array field on the host.

rms

rms(remove='piston')

RMS of the retrieved OPD in metres over the pupil (see :func:~solvephase.rms).

residuals

residuals(images)

images - model_images on the host.

summary

summary()

A short human-readable report.

save

save(path)

Save arrays and metadata to a compressed .npz file.

plot

plot(images=None, *, unit='nm', path=None)

Plot the OPD, and data/model/residual images if images is given.

Needs matplotlib (pip install 'solvephase[plot]'). Returns the figure.

Optics

aocore.Pupil dataclass

Pupil(amplitude, pitch, diameter, obscuration=0.0, segments=None, name='custom')

A sampled entrance pupil.

Attributes:

Name	Type	Description
amplitude	ndarray	(ny, nx) real amplitude transmission, >= 0 (host NumPy).
pitch	float	Pixel pitch in metres.
diameter	float	Reference diameter in metres. It sets lambda / D and the unit radius of Zernike polynomials.
obscuration	float	Central obscuration as a fraction of diameter (used by annular Zernikes; 0 for an unobscured pupil).
segments	ndarray | None	Optional (ny, nx) integer label image: 0 outside, 1..S for each segment of a segmented aperture.
name	str	A short descriptive label.

shape property

shape

Grid shape (ny, nx).

mask property

mask

Boolean support: pixels with non-zero amplitude.

n_valid property

n_valid

Number of pixels in :attr:mask.

area property

area

Transmitted area sum(amplitude^2) * pitch^2 in square metres.

n_segments property

n_segments

Number of labelled segments (0 for a monolithic pupil).

coordinates

coordinates()

Return (y, x) coordinate grids in metres.

lambda_over_d

lambda_over_d(wavelength)

wavelength / diameter in radians.

with_amplitude

with_amplitude(amplitude)

Copy with a different amplitude map on the same grid.

padded

padded(n)

Embed the pupil centred in a larger n x n grid (same pitch).

downsampled

downsampled(factor)

Average factor x factor pixel blocks (pitch grows by factor).

from_array classmethod

from_array(amplitude, *, diameter=None, pitch=None, obscuration=0.0, segments=None, name='custom')

Wrap an existing amplitude map.

Give diameter (the grid spans it across the largest illuminated extent if pitch is omitted), pitch, or both. With neither, the grid width is taken as 1 m.

circular classmethod

circular(n, diameter=1.0, *, obscuration=0.0, spiders=0, spider_width=0.0, spider_angle=0.0, fill=1.0, supersample=4, name=None)

Circular, optionally obscured pupil with radial spider vanes.

Parameters:

Name	Type	Description	Default
n	int	Grid size in pixels (n x n).	required
diameter	float	Outer diameter in metres.	1.0
obscuration	float	Central obscuration diameter as a fraction of diameter.	0.0
spiders	int	Number of radial vanes, evenly spaced from spider_angle.	0
spider_width	float	Vane width in metres.	0.0
spider_angle	float	Angle of the first vane in degrees, counter-clockwise from +x.	0.0
fill	float	Fraction of the grid width the diameter spans (<= 1).	1.0
supersample	int	Sub-samples per pixel axis for anti-aliased edges (1 = binary).	4

segmented_hexagonal classmethod

segmented_hexagonal(n, *, rings, segment_size, gap=0.0, missing_center=True, obscuration=0.0, spiders=0, spider_width=0.0, spider_angle=90.0, orientation=0.0, fill=1.0, supersample=4, name='segmented')

Hexagonally segmented aperture (Keck-, JWST-, ELT-like).

Parameters:

Name	Type	Description	Default
rings	int	Number of segment rings around the centre.	required
segment_size	float	Segment flat-to-flat width in metres.	required
gap	float	Gap between neighbouring segments in metres.	0.0
missing_center	bool	Leave out the central segment.	True
obscuration	float	Circular central obscuration as a fraction of the circumscribed diameter.	0.0
orientation	float	Rotation of the segment lattice in degrees.	0.0

Notes

The returned pupil's :attr:~Pupil.segments labels the segments 1..S, ordered by ring and then angle; :attr:~Pupil.diameter is the circumscribed diameter.

keck classmethod

keck(n=256, *, supersample=4)

Keck-like primary: 36 pointy-top hexagonal segments, 10.95 m across corners.

The segment pitch is 10.95 / 7 m (1.8 m corner to corner) with 3 mm gaps; six 26 mm support arms at 30 + 60 k degrees and a 2.57 m circular central obscuration. :attr:diameter is the conventional 10.95 m (the circumscribed circle is 11.19 m). This follows the geometry of makewfs' Keck HAKA example, without its fitted hexagonal secondary shadow and offset. Good for algorithm development, not instrument modelling.

jwst classmethod

jwst(n=256, *, supersample=4)

JWST-like primary: 18 hexagonal segments of 1.32 m with three struts.

Approximate geometry (7 mm gaps, 0.74 m central obscuration, three 0.1 m struts) for algorithm development.

vlt classmethod

vlt(n=256, *, supersample=4)

VLT Unit Telescope: 8.0 m with a 1.12 m obscuration and four vanes.

hst classmethod

hst(n=256, *, supersample=4)

Hubble-like: 2.4 m, 0.33 obscuration and four 2.6 cm vanes.

solvephase.Basis

Basis(pupil, modes, *, labels=None, variances=None, kind='custom')

A linear wavefront parameterization over a :class:~solvephase.pupil.Pupil.

Build one with a constructor (:meth:zernike, :meth:kl, :meth:zonal, ...), then use :meth:synthesize and :meth:fit. Solvers move the basis to their device with :meth:on.

Attributes:

Name	Type	Description
pupil		The pupil whose valid pixels the modes are defined on.
modes		(n_modes, n_valid) host matrix, or None for the zonal basis.
labels		One short name per mode.
variances		Optional prior variance of each coefficient in m^2 (KL eigenvalues), used for maximum a posteriori regularization.
kind		"zernike", "kl", "zonal", "segments", "dm", ...

is_zonal property

is_zonal

Whether this is the identity (pixel) basis.

n_modes property

n_modes

Number of coefficients.

mode_maps

mode_maps(indices=None)

Return modes as (k, ny, nx) maps (zero outside the pupil).

synthesize

synthesize(coefficients)

OPD map(s) B c for coefficient vector(s) (..., n_modes) (host).

fit

fit(opd, *, weights=None, rcond=1e-10)

Least-squares coefficients of an OPD map over the pupil (host).

weights defaults to the pupil intensity amplitude**2.

project

project(opd)

Best approximation of an OPD map within the span of the basis.

on

on(backend)

The (n_modes, n_valid) mode matrix on backend (cached).

concatenate

concatenate(*others)

Stack the modes of several bases defined on the same pupil.

subset

subset(indices)

A basis made of some of these modes.

orthonormalized

orthonormalized()

Orthonormal over the pupil (intensity-weighted, unit RMS), same span and order.

zonal classmethod

zonal(pupil)

One coefficient per illuminated pixel (OPD in metres).

from_maps classmethod

from_maps(pupil, maps, *, labels=None, variances=None, kind='custom')

Wrap (n_modes, ny, nx) mode maps (values outside the pupil are ignored).

zernike classmethod

zernike(pupil, n_modes, *, start=2, ordering='noll', annular=None, orthonormalize=False)

Zernike polynomials (via :class:aobasis.ZernikeBasisGenerator).

Parameters:

Name	Type	Description	Default
n_modes	int	Number of modes returned.	required
start	int	First index in ordering (Noll: 1 = piston, 2 = tip, 4 = defocus). The default 2 skips piston, which phase retrieval cannot see.	2
ordering	str	"noll", "ansi" or "fringe".	'noll'
annular	bool | None	Use annular Zernikes orthonormal over the obscured pupil. Defaults to True when the pupil has an :attr:~Pupil.obscuration.	None
orthonormalize	bool	Re-orthonormalize over the sampled pupil (spiders, gaps, edge pixels), which improves solver conditioning. Coefficients then belong to the orthonormalized modes, not the textbook polynomials.	False

Notes

Mode labels are Z<j> in the chosen ordering, and each mode has unit RMS over the continuous (annular) disk, so a coefficient is RMS OPD in metres.

from_aobasis classmethod

from_aobasis(pupil, generator, n_modes, *, kind=None, **generate_kwargs)

Evaluate any :mod:aobasis generator class at the pupil's pixels.

generator is a generator class (e.g. aobasis.FourierBasisGenerator) or a factory f(positions) -> generator; it is constructed with the pupil's pixel positions in metres. Keywords go to generate.

fourier classmethod

fourier(pupil, n_modes, *, orthonormalize=True)

Sine/cosine modes (via :class:aobasis.FourierBasisGenerator), piston-free.

kl classmethod

kl(pupil, n_modes, *, r0=0.15, outer_scale=30.0, r0_wavelength=5e-07, max_points=2500)

Karhunen-Loeve modes of von Karman turbulence (via :mod:aobasis).

The covariance is diagonalized on a coarse grid of at most max_points pupil samples and the modes are interpolated (cubic) to the full pupil, then re-orthonormalized. :attr:variances holds each coefficient's turbulence variance in m^2 (OPD), so the basis can act as a Bayesian prior.

Parameters:

Name	Type	Description	Default
r0	float	Fried parameter in metres at r0_wavelength.	0.15
outer_scale	float	Outer scale in metres (numpy.inf for Kolmogorov).	30.0

segments classmethod

segments(pupil, terms=('piston', 'tip', 'tilt'))

Per-segment piston, tip and tilt for a segmented pupil.

Each mode is non-zero on one segment and unit RMS over it (piston = 1, tip/tilt are ramps in x/y about the segment centroid). Requires :attr:Pupil.segments.

influence_functions classmethod

influence_functions(pupil, influence, *, m2c=None, labels=None)

Deformable-mirror basis: coefficients are DM commands.

Parameters:

Name	Type	Description	Default
influence	Any	(n_actuators, ny, nx) influence-function OPD maps in metres per unit command (zero outside the pupil is not required).	required
m2c	Any	Optional (n_actuators, n_modes) modal-to-command matrix (e.g. from :mod:aobasis); coefficients are then modal commands.	None

gaussian_dm classmethod

gaussian_dm(pupil, actuators_across, *, coupling=0.15, stroke=1e-06)

Square-grid DM with Gaussian influence functions (via :mod:aobasis).

Actuators are on a actuators_across-wide square grid spanning the pupil diameter (rim actuators on the edge), keeping those within one pitch of the pupil. A unit command gives stroke metres of OPD at the actuator.

solvephase.FocalPlaneModel

FocalPlaneModel(pupil, wavelength, image_shape, *, pixel_scale=None, sampling=None, weights=None, reference_wavelength=None, diversity=None, oversample=1, offset=0.0, method='auto', device='cpu', precision=None)

Pupil-to-detector forward model for one or more diversity channels.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The entrance pupil.	required
wavelength	float | Sequence[float]	Wavelength in metres, or an array of wavelengths sampling a band.	required
image_shape	Any	Detector image shape (my, mx) (an int for square images).	required
pixel_scale	float | None	Detector pixel scale in radians. Give this or sampling.	None
sampling	float | None	Detector pixels per lambda_ref / D (2 = Nyquist), an alternative to pixel_scale.	None
weights	Sequence[float] | None	Relative spectral weights (photon flux) of the wavelengths; normalized to sum to 1. Uniform by default.	None
reference_wavelength	float | None	Wavelength at which phase is expressed. Defaults to the weight-averaged wavelength.	None
diversity	Any	Known OPD added in each channel, (K, ny, nx) in metres (or a sequence of maps; None entries mean zero). None gives a single channel with no diversity.	None
oversample	int	Model each detector pixel with oversample x oversample samples and integrate (sum) them, which models the pixel's response.	1
offset	Any	Detector-window offset in pixels (dy, dx): pixel j samples the angle (j - (m - 1) / 2 + offset) * pixel_scale, so the optical axis falls on pixel (m - 1) / 2 - offset. offset=-0.5 puts it on pixel m / 2 of an even-sized image (the HCIPy/fftshift convention).	0.0
method	str	Propagation engine: "auto", "fft" or "mft".	'auto'
device	BackendLike	Backend selection, see :func:solvephase.get_backend.	'cpu'
precision	BackendLike	Backend selection, see :func:solvephase.get_backend.	'cpu'

n_channels property

n_channels

Number of diversity channels K.

n_wavelengths property

n_wavelengths

Number of wavelength samples L.

sampling property

sampling

Detector pixels per lambda_ref / D.

nyquist_sampled property

nyquist_sampled

Whether the shortest wavelength has at least 2 pixels per lambda/D.

opd_to_phase

opd_to_phase(opd)

Convert OPD in metres to phase in radians at the reference wavelength.

phase_to_opd

phase_to_opd(phase)

Convert reference-wavelength phase in radians to OPD in metres.

forward

forward(phase, amplitude=None)

Model images for a phase map.

Parameters:

Name	Type	Description	Default
phase	Any	(ny, nx) phase in radians at the reference wavelength (shared by all channels), or (K, ny, nx) per channel. Backend array.	required
amplitude	Any	Optional (ny, nx) pupil amplitude; defaults to the pupil's.	None

Returns:

Type	Description
ForwardState	state.images holds the (K, my, mx) normalized images.

images

images(opd=None, *, amplitude=None)

Convenience: normalized (K, my, mx) images for an OPD map in metres.

Accepts NumPy or backend arrays; returns a backend array.

vjp

vjp(state, grad_images, *, amplitude=False)

Back-propagate dL/d images to the phase (and optionally amplitude).

Parameters:

Name	Type	Description	Default
state	ForwardState	The :class:ForwardState of the matching :meth:forward call.	required
grad_images	Any	(K, my, mx) real gradient of a scalar with respect to state.images.	required
amplitude	bool	Also return the amplitude gradient.	False

Returns:

Type	Description
(grad_phase, grad_amplitude)	(K, ny, nx) per-channel phase gradient (sum over channels for a shared phase), and (ny, nx) amplitude gradient or None.

jvp

jvp(state, directions, *, per_channel=False)

Image derivatives along phase directions.

Parameters:

Name	Type	Description	Default
state	ForwardState	The :class:ForwardState at the linearization point.	required
directions	Any	(P, ny, nx) phase directions in radians (applied to every channel), or with per_channel=True (P, K, ny, nx).	required

Returns:

Type	Description
``(P, K, my, mx)`` array of ``d images / d direction``.

solvephase.zernike_diversity

zernike_diversity(pupil, noll, amounts, *, annular=None)

Diversity OPD maps (K, ny, nx): one Zernike mode at several RMS amplitudes.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The pupil.	required
noll	int	Noll index of the mode (4 = defocus, 5/6 = astigmatism).	required
amounts	Sequence[float]	RMS OPD of the mode for each channel, in metres (0 for an in-focus channel).	required

Examples:

Two images, in focus and with 0.5 waves RMS of defocus at 1.6 um:

>>> div = zernike_diversity(pupil, 4, [0.0, 0.8e-6])

Focal-plane solvers

solvephase.FocalPlaneProblem

FocalPlaneProblem(model, images, *, basis=None, loss='gaussian', weights=None, read_noise=0.0, fit_flux=True, flux=None, fit_background=False, background=0.0, fit_tilt=False, fit_amplitude=False, regularization=0.0, prior=False, smoothness=0.0)

Phase retrieval from one or more focal-plane images.

Parameters:

Name	Type	Description	Default
model	FocalPlaneModel	The forward model (pupil, wavelengths, sampling, diversity channels).	required
images	Any	Measured images (K, my, mx) (or (my, mx) for one channel), in any linear units (ADU, photo-electrons). Use photo-electrons with the Poisson loss.	required
basis	Basis | int | None	Wavefront parameterization. None means zonal (one value per pupil pixel); an int n means Basis.zernike(pupil, n).	None
loss	str | Loss	"gaussian" (weighted least squares), "poisson" (maximum likelihood for photon-counting data) or "amplitude"; or a :class:~solvephase.losses.Loss.	'gaussian'
weights	Any	Per-pixel weights (K, my, mx), broadcastable; zero masks a pixel. For the Gaussian loss, inverse variances make it a likelihood (see :func:noise_weights). Defaults to 1.	None
read_noise	float	Read noise in data units, used by the Poisson loss.	0.0
fit_flux	bool	Fit each channel's total flux. Otherwise flux is used.	True
flux	Any	Fixed per-channel flux (scalar or (K,)) when fit_flux is False; by default the data sum.	None
fit_background	bool	Fit a constant background per channel. Otherwise background is used.	False
background	Any	Fixed per-channel background (data units per pixel). Default 0.	0.0
fit_tilt	bool	Fit a tip/tilt per channel relative to channel 0 (image registration); requires at least two channels.	False
fit_amplitude	bool	Also retrieve the pupil amplitude (zonal, on the pupil support).	False
regularization	float	Tikhonov weight on the phase coefficients (radians^-2).	0.0
prior	bool	Add the Gaussian prior 0.5 * sum c^2 / var from :attr:Basis.variances (e.g. a KL basis), i.e. maximum a posteriori.	False
smoothness	float	Weight of a first-difference smoothness penalty on a zonal phase.	0.0

n_params property

n_params

Length of the internal parameter vector.

n_channels property

n_channels

Number of image channels.

initial

initial(start=None)

Internal starting vector.

start may be None (flat wavefront), an OPD map in metres (ny, nx), a coefficient vector in metres, or a :class:Result (its OPD, amplitude, flux and background are reused).

value

value(x)

Objective value at x (data term + regularization).

objective

objective(x)

Objective value and its gradient with respect to the internal vector x.

gauss_newton

gauss_newton(x, *, chunk=None)

Objective, gradient and Gauss-Newton (Fisher) matrix at x.

Builds the image Jacobian column-block by column-block with forward-mode derivatives; chunk directions are propagated at once (default: as many as fit in about 64 MB of propagation workspace). Not available with fit_amplitude or a zonal basis (too many parameters); use L-BFGS there.

result

result(x, *, method, opt=None, elapsed=0.0, extra=None)

Package an internal vector as a :class:~solvephase.result.Result.

solvephase.solve

solve(problem, *, method='auto', start=None, max_iter=None, callback=None, **options)

Minimize a :class:FocalPlaneProblem.

Parameters:

Name	Type	Description	Default
problem	FocalPlaneProblem	The problem to solve.	required
method	str	"lbfgs" (default for zonal or amplitude problems), "lm" (Levenberg-Marquardt; default for modal problems with at most 300 parameters), "adam", or "auto".	'auto'
start	Any	Starting point, see :meth:FocalPlaneProblem.initial.	None
max_iter	int | None	Iteration limit (defaults: 300 for L-BFGS, 50 for LM, 1000 for Adam).	None
callback	Callable[[int, Any, float], bool | None] | None	callback(iteration, x, f); return True to stop.	None
**options	Any	Passed to the optimizer (ftol, gtol, memory, damping, learning_rate ...).	{}

solvephase.gerchberg_saxton

gerchberg_saxton(model, images, *, iterations=200, start=None, basis=None, background=0.0, weights=None, momentum=0.0, tol=1e-06, check_every=10, unwrap=True, callback=None)

Gerchberg-Saxton / Misell phase retrieval with a known pupil amplitude.

Parameters:

Name	Type	Description	Default
model	FocalPlaneModel	Forward model. Broadband models are run at the reference wavelength (a warning is issued when the fractional bandwidth exceeds 5%).	required
images	Any	(K, my, mx) measured images (any linear units).	required
iterations	int	Maximum number of iterations.	200
start	Any	Initial OPD map in metres (default flat).	None
basis	Basis | int | None	Optional basis (or number of Zernike modes) the unwrapped phase is fitted to; the returned OPD is then the modal fit.	None
background	Any	Per-channel background subtracted before taking square roots.	0.0
weights	Any	Per-pixel weights; zero marks unmeasured pixels, whose model value is kept.	None
momentum	float	Over-relaxation beta of the accelerated GS update theta <- P(theta) + beta (P(theta) - P(theta_prev)) (0 = plain GS).	0.0
tol	float	Stop when the relative change of the error metric over check_every iterations is below tol.	1e-06
check_every	int	Error-metric interval (each check synchronizes the GPU).	10
unwrap	bool	Unwrap the final phase (otherwise the OPD is the wrapped phase).	True
callback	Callable[[int, Any, float], bool | None] | None	callback(iteration, phase, error); return True to stop.	None

Returns:

Type	Description
Result	extra["wrapped_phase"] holds the wrapped pupil phase and history the normalized amplitude error sum (|E| - sqrt(d))^2 / sum d over measured pixels.

solvephase.noise_weights

noise_weights(images, read_noise=0.0, *, gain=1.0, floor=1.0)

Inverse-variance weights 1 / (gain * max(d, 0) + read_noise^2) for a Gaussian loss.

images are in photo-electrons (or data units with Poisson gain); floor bounds the variance from below.

Algorithm families

solvephase.algorithms.phase_diversity

Phase diversity with an unknown extended object.

K images d_k = o * s_k + n_k of the same unknown scene o are recorded through the unknown wavefront plus a known diversity per channel (defocus, astigmatism ...). Gonsalves (Opt. Eng. 21, 829, 1982) showed that the object can be eliminated analytically. In Fourier space, with D_k = F d_k and the optical transfer functions S_k = F s_k, the least-squares object for fixed S_k is the multi-frame Wiener estimate

O = sum_k D_k conj(S_k) / (sum_k |S_k|^2 + gamma),

and substituting it back leaves a metric of the wavefront alone (Paxman, Schulz & Fienup, JOSA A 9, 1072, 1992):

L = sum_f [ sum_k |D_k|^2 - |sum_k D_k conj(S_k)|^2 / (sum_k |S_k|^2 + gamma) ].

gamma regularizes the division (the inverse signal-to-noise power ratio of a Wiener filter). Following Löfdahl & Scharmer (A&AS 107, 243, 1994) the data are apodized against periodic-convolution edge effects and the sum is restricted to spatial frequencies inside the diffraction cutoff, where the transfer functions carry information. Mugnier, Blanc & Idier (Adv. Imaging Electron Phys. 141, 1, 2006) review the method.

:class:PhaseDiversityProblem evaluates L and its exact gradient with respect to the modal coefficients: dL/dS_k follows from the envelope theorem (O is optimal, so its own variation drops out), the real FFT is back-propagated to dL/d psf_k, and :meth:~solvephase.focal.FocalPlaneModel.vjp takes it to the pupil phase. The metric is evaluated in the equivalent form sum_k |D_k - O S_k|^2 + gamma |O|^2, a sum of non-negative terms that keeps its accuracy in single precision. :func:phase_diversity minimizes it with L-BFGS, optionally coarse-to-fine over the number of modes.

Internal parameters are those of :class:~solvephase.retrieval.FocalPlaneProblem: phase coefficients in radians at the reference wavelength (reported in metres), a relative zonal amplitude, and per-channel tip/tilt for registration.

PhaseDiversityProblem

PhaseDiversityProblem(model, images, *, basis=None, regularization='auto', object_regularization='auto', window='tukey', window_alpha=0.25, frequency_mask='diffraction', fit_amplitude=False, fit_tilt=False, equalize_flux=True)

Phase diversity of an unknown extended object, with the object eliminated.

Parameters:

Name	Type	Description	Default
model	FocalPlaneModel	Forward model with one channel per image. Its image_shape must be the data shape; its normalized PSFs (sum ~1) are the convolution kernels. Sample at Nyquist or better at the shortest wavelength.	required
images	Any	(K, my, mx) images of the same scene, K >= 2, in any linear units (a constant background is absorbed by the object). Channels must be registered to a pixel or better (see fit_tilt).	required
basis	Basis | int | None	Wavefront parameterization. None means zonal; an int n means Basis.zernike(pupil, n, start=4) (defocus upwards), because a common tip/tilt only moves the unknown object and is not measurable.	None
regularization	float | str	gamma of the metric, added to sum_k |S_k|^2 (dimensionless: the OTFs are 1 at zero frequency). It only needs to stabilize the division where every OTF vanishes; larger values bias the wavefront towards sharper PSFs. "auto" is the measured noise power over the peak in-band data power (the inverse of the best per-frequency SNR, squared).	'auto'
object_regularization	float | str	gamma of the reported object estimate (a Wiener filter, which needs more damping than the metric). "auto" is the noise power over the mean in-band data power.	'auto'
window	Any	Apodization: each image minus its edge level (the mean weighted by 1 - window) is multiplied by the window and the level added back, so the frame is periodic and a uniform background passes untouched. "tukey" (cosine-tapered edges, flat centre), "hann", None or a (my, mx) array. Use None only for scenes whose images already fade to a uniform level at the frame edges.	'tukey'
window_alpha	float	Tapered fraction of each axis for the Tukey window.	0.25
frequency_mask	str | float | None	"diffraction" sums the metric over spatial frequencies inside the diffraction cutoff D / lambda_min; a float a uses a times that radius; None uses every frequency. Zero frequency (the total flux, which carries no wavefront information) is always excluded.	'diffraction'
fit_amplitude	bool	Also retrieve the pupil amplitude (zonal, on the pupil support).	False
fit_tilt	bool	Fit a tip/tilt per channel relative to channel 0 (image registration).	False
equalize_flux	bool	Scale every image to the mean total flux, so unequal exposures do not bias the metric.	True

Notes

The objective is the reduced metric divided by the noise power per frequency, so it is a chi-square-like number of order (K - 1) x (in-band frequencies) at the solution when the noise power can be measured outside the cutoff (otherwise it is divided by the mean in-band data power).

n_params property

n_params

Length of the internal parameter vector.

n_channels property

n_channels

Number of images K.

with_basis

with_basis(basis)

The same data and settings with another basis (used for coarse-to-fine solves).

initial

initial(start=None)

Internal starting vector.

start may be None (flat wavefront), an OPD map in metres (ny, nx), a coefficient vector in metres, or a :class:Result (its OPD, and its amplitude and tilts when fitted, are reused).

value

value(x)

Reduced metric at x (normalized by the noise power).

objective

objective(x)

Reduced metric and its gradient with respect to the internal vector x.

object_estimate

object_estimate(x)

Wiener object estimate (my, mx) at x on the data grid, in data units.

Uses :attr:object_regularization. The estimate is registered with the optical axis, band-limited to the frequency mask (plus the mean), and carries the apodization window and the flux equalization.

result

result(x, *, opt=None, elapsed=0.0, method='phase_diversity')

Package an internal vector as a :class:~solvephase.result.Result.

model_images is the model of the apodized, flux-equalized data (extra["apodized_images"]). extra also holds "object" (:meth:object_estimate), "otf" ((K, my, mx // 2 + 1) real-FFT transfer functions, 1 at zero frequency for an unbounded detector), "regularization" and "object_regularization" (the gammas), "noise_power", "window" and "frequency_mask" (real-FFT layout).

phase_diversity

phase_diversity(model, images, *, basis=None, regularization='auto', object_regularization='auto', window='tukey', window_alpha=0.25, frequency_mask='diffraction', fit_amplitude=False, fit_tilt=False, equalize_flux=True, start=None, coarse_to_fine=None, max_iter=300, callback=None, **options)

Retrieve a wavefront from phase-diverse images of an unknown extended object.

Builds a :class:PhaseDiversityProblem (see it for basis, regularization, object_regularization, window, window_alpha, frequency_mask, fit_amplitude, fit_tilt and equalize_flux) and minimizes the reduced metric with L-BFGS.

Parameters:

Name	Type	Description	Default
model	FocalPlaneModel	Forward model with one channel per image (diversity included).	required
images	Any	(K, my, mx) images of the same scene.	required
start	Any	Starting point: None (flat), an OPD map or coefficient vector in metres, or a :class:~solvephase.result.Result.	None
coarse_to_fine	Sequence[int] | None	Increasing mode counts, e.g. (5, 10): solve with the first n modes of basis, then continue with more, ending with the full basis. Widens the capture range for large aberrations. Needs a modal basis.	None
max_iter	int	L-BFGS iteration limit per stage.	300
callback	Callable[[int, Any, float], bool | None] | None	callback(iteration, x, f); return True to stop the stage.	None
**options	Any	Passed to :func:solvephase.optimize.lbfgs (ftol, gtol, memory ...).	{}

Returns:

Type	Description
Result	method="phase_diversity"; coefficients are RMS OPD in metres and extra["object"] is the Wiener object estimate on the data grid, in data units (see :meth:PhaseDiversityProblem.result). A common tip/tilt is not measurable and stays at its starting value.

Examples:

>>> pupil = Pupil.circular(64)
>>> div = zernike_diversity(pupil, 4, [0.0, 0.29 * 1.6e-6])  # 1 wave PV
>>> model = FocalPlaneModel(pupil, 1.6e-6, 128, sampling=2.0, diversity=div)
>>> result = phase_diversity(model, images, basis=21)

solvephase.algorithms.lift

LIFT: low-order wavefront sensing from one astigmatic focal-plane image.

LIFT (LInearized Focal-plane Technique; Meimon, Fusco & Mugnier, Opt. Lett. 35, 3036, 2010; Plantet et al., Opt. Express 21, 16337, 2013) estimates a handful of low-order modes (tip, tilt, defocus, astigmatism, coma ...) from a single, possibly very faint, focal-plane image. A known astigmatism is added in the optical path. Without it a single in-focus image cannot tell the sign of the even part of the wavefront (defocus, astigmatism, spherical...): the twin phi(x) -> -phi(-x) gives the same image. The astigmatism moves that twin far from the true solution (for the total phase phi + phi_d the twin needs phi_d -> -phi_d), so a local maximum-likelihood search started near zero lands on the true wavefront.

The estimator is maximum likelihood under photon and read noise, computed by iterated linearization of the PSF around the current estimate (Gauss-Newton, with the noise variance model + read_noise^2 as weights). In solvephase this is exactly :func:solvephase.solve with method="lm" on a :class:~solvephase.FocalPlaneProblem with a Poisson (shifted-Poisson) loss, a modal basis and a fitted flux; this module only sets that problem up, keeps the model for repeated estimates (:class:LIFT), and computes the Cramer-Rao lower bound of the modal coefficients (:func:lift_crlb), which LIFT approaches at high flux.

Units: coefficients are RMS OPD in metres of unit-RMS Zernike modes (Noll order, starting at tip by default); the astigmatism is RMS OPD in metres; images are in photo-electrons (read noise in electrons RMS).

LIFT

LIFT(pupil, wavelength, image_shape, *, sampling=None, pixel_scale=None, n_modes=10, basis=None, astigmatism=None, astigmatism_mode=5, read_noise=0.0, loss='poisson', fit_flux=True, fit_background=False, background=0.0, max_iter=20, resolve_twin=True, oversample=1, offset=0.0, device='cpu', precision=None)

A reusable LIFT wavefront sensor: the model is built once, estimates are repeated.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The entrance pupil.	required
wavelength	Any	Wavelength in metres (or an array sampling a band, see :class:~solvephase.FocalPlaneModel).	required
image_shape	Any	Detector image shape (my, mx).	required
sampling	float | None	Detector pixels per lambda / D or pixel scale in radians (give one).	None
pixel_scale	float | None	Detector pixels per lambda / D or pixel scale in radians (give one).	None
n_modes	int	Number of Zernike modes estimated, Noll order from tip (Noll 2). Ignored when basis is given.	10
basis	Basis | None	Any modal :class:~solvephase.Basis instead of Zernikes.	None
astigmatism	float | None	RMS OPD of the known astigmatism in metres (default wavelength / 10, see :func:default_astigmatism).	None
astigmatism_mode	int	Noll index of the diversity astigmatism (5 or 6).	5
read_noise	float	Read-noise standard deviation in electrons.	0.0
loss	str | Loss	"poisson" (maximum likelihood under photon + read noise; the default) or "gaussian" (weighted least squares with variances from the image, :func:~solvephase.noise_weights, as in the original LIFT); or a :class:~solvephase.Loss.	'poisson'
fit_flux	bool	Nuisance parameters, see :class:~solvephase.FocalPlaneProblem.	True
fit_background	bool	Nuisance parameters, see :class:~solvephase.FocalPlaneProblem.	True
background	bool	Nuisance parameters, see :class:~solvephase.FocalPlaneProblem.	True
max_iter	int	Gauss-Newton (Levenberg-Marquardt) iteration limit per estimate.	20
resolve_twin	bool	After convergence, compare the solution with its exact twin (same image, see :meth:twin) and report the one with the smaller wavefront RMS. Only active when the twin is exact: centro-symmetric pupil, centred image window and a basis closed under the twin (Zernikes are).	True
oversample	int	Detector pixel integration and optical-axis offset, see :class:~solvephase.FocalPlaneModel.	1
offset	int	Detector pixel integration and optical-axis offset, see :class:~solvephase.FocalPlaneModel.	1
device	BackendLike	Backend selection.	'cpu'
precision	BackendLike	Backend selection.	'cpu'

Examples:

>>> sensor = LIFT(pupil, 1.6e-6, 32, sampling=2, n_modes=10)
>>> for frame in frames:
...     result = sensor.estimate(frame)
...     dm_command -= gain * result.coefficients

backend property

backend

The model's :class:~solvephase.Backend.

n_modes property

n_modes

Number of estimated modal coefficients.

twin

twin(coefficients)

The twin coefficients (metres) that produce exactly the same image.

Raises ValueError when the twin is not exact for this sensor (an asymmetric pupil, an off-centre window, or a basis not closed under phi(x) -> -phi(-x)).

problem

problem(image)

The :class:~solvephase.FocalPlaneProblem LIFT solves for image.

estimate

estimate(image, start=None, *, max_iter=None, **options)

Estimate the modal coefficients from one image.

Parameters:

Name	Type	Description	Default
image	Any	(my, mx) image in photo-electrons (NumPy or backend array).	required
start	Any	Starting point: None (flat), coefficients in metres, an OPD map or a previous :class:~solvephase.Result (warm start in a loop).	None
max_iter	int | None	Override the iteration limit.	None
**options	Any	Passed to the Levenberg-Marquardt solver (ftol, damping).	{}

Returns:

Type	Description
Result	method="lift"; coefficients in metres RMS per mode, the known astigmatism is not included in opd.

expected_image

expected_image(coefficients=None, *, photons=1.0)

Noise-free image photons * PSF + background (backend array).

coefficients are in metres (default: zero wavefront).

fisher

fisher(coefficients=None, *, photons)

Fisher information matrix at coefficients (metres) for photons detected.

Covers the modal coefficients and, when fitted, the nuisance parameters (log-flux, background), in that order, in SI units for the coefficients (m^-2). Pixel noise is shifted Poisson: variance model + read_noise^2.

crlb

crlb(coefficients=None, *, photons)

Cramer-Rao lower bound: covariance (n_modes, n_modes) of the coefficients in m^2.

Nuisance parameters (flux, background) are marginalized: the bound is the coefficient block of the inverse Fisher matrix.

lift_diversity

lift_diversity(pupil, wavelength, astigmatism=None, *, noll=5)

The known LIFT astigmatism as an OPD map (ny, nx) in metres.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The pupil.	required
wavelength	float	Wavelength in metres (sets the default amplitude).	required
astigmatism	float | None	RMS OPD of the astigmatism in metres. None gives :func:default_astigmatism (wavelength / 10 RMS).	None
noll	int	Noll index of the astigmatism mode, 5 (oblique) or 6 (vertical).	5

default_astigmatism

default_astigmatism(wavelength)

Default LIFT astigmatism: wavelength / 10 RMS OPD (about half a wave peak to valley).

A unit-RMS Zernike astigmatism spans 2 sqrt 6 peak to valley over the disk, so this is 0.49 waves PV. It is close to the 170 nm RMS at 1.6 um of VLT/IRLOS (Kuznetsov et al. 2024). The Cramer-Rao bound is flat between about 0.05 and 0.15 waves RMS, while the capture range grows with the astigmatism (the twin solution sits at -a_astig - 2 d).

lift

lift(image, pupil, wavelength, *, sampling=None, pixel_scale=None, n_modes=10, basis=None, astigmatism=None, astigmatism_mode=5, read_noise=0.0, loss='poisson', fit_background=False, background=0.0, start=None, max_iter=20, device='cpu', precision=None, **options)

Estimate low-order modes from one astigmatic focal-plane image (LIFT).

Parameters:

Name	Type	Description	Default
image	Any	(my, mx) image in photo-electrons, taken with the known astigmatism in the beam, background subtracted (or set background / fit_background).	required
pupil	Pupil	The optics, see :class:~solvephase.FocalPlaneModel.	required
wavelength	Pupil	The optics, see :class:~solvephase.FocalPlaneModel.	required
sampling	Pupil	The optics, see :class:~solvephase.FocalPlaneModel.	required
pixel_scale	Pupil	The optics, see :class:~solvephase.FocalPlaneModel.	required
n_modes	int	Zernike modes from tip (Noll 2), or any modal basis.	10
basis	int	Zernike modes from tip (Noll 2), or any modal basis.	10
astigmatism	float | None	Known astigmatism, RMS OPD in metres (default wavelength / 10).	None
astigmatism_mode	int	5 (oblique) or 6 (vertical astigmatism).	5
read_noise	float	Read noise in electrons RMS.	0.0
loss	str | Loss	"poisson" (default, maximum likelihood) or "gaussian".	'poisson'
start	Any	Optional starting point (coefficients, OPD or Result).	None
max_iter	int	Gauss-Newton iteration limit.	20
device	BackendLike	Backend selection.	'cpu'
precision	BackendLike	Backend selection.	'cpu'
**options	Any	Passed to the Levenberg-Marquardt solver.	{}

Returns:

Type	Description
Result	With coefficients (metres RMS per mode) and opd (without the known astigmatism).
For repeated estimates with the same optics use :class:`LIFT`, which
builds the model once.

lift_crlb

lift_crlb(pupil, wavelength, image_shape, coefficients=None, *, photons, read_noise=0.0, sampling=None, pixel_scale=None, n_modes=10, basis=None, astigmatism=None, astigmatism_mode=5, background=0.0, fit_flux=True, device='cpu', precision=None)

Cramer-Rao lower bound on the LIFT modal coefficients.

The Fisher matrix J^T diag(1 / (m + read_noise^2)) J is built from the model Jacobian J (forward-mode derivatives of the expected image m with respect to the coefficients and the flux), at the true coefficients. Its inverse bounds the covariance of any unbiased estimator.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The optics.	required
wavelength	Pupil	The optics.	required
image_shape	Pupil	The optics.	required
sampling	Pupil	The optics.	required
pixel_scale	Pupil	The optics.	required
coefficients	Any	True coefficients in metres (default: zero wavefront).	None
photons	float	Total detected photo-electrons in the PSF.	required
read_noise	float	Read noise in electrons RMS.	0.0
n_modes	int	As for :func:lift.	10
basis	int	As for :func:lift.	10
astigmatism	int	As for :func:lift.	10
astigmatism_mode	int	As for :func:lift.	10
background	float	Known background in electrons per pixel (adds photon noise).	0.0
fit_flux	bool	Whether the flux is an unknown (marginalized) parameter.	True
device	BackendLike	Backend selection (double precision is advisable for the inverse).	'cpu'
precision	BackendLike	Backend selection (double precision is advisable for the inverse).	'cpu'

Returns:

Type	Description
ndarray	(n_modes, n_modes) covariance bound in m^2; the square root of its diagonal is the smallest achievable RMS error per mode.

solvephase.algorithms.fast_furious

Fast & Furious: sequential small-phase focal-plane wavefront sensing.

Fast & Furious (F&F; Keller et al., Proc. SPIE 8447, 2012; Korkiakoski et al., Appl. Opt. 53, 4565, 2014; on sky: Bos et al., A&A 639, A52, 2020) measures a pupil-plane wavefront from focal-plane images in a closed loop with a deformable mirror (DM). It builds on the weak-phase solution of Gonsalves (Opt. Lett. 26, 684, 2001).

For a weak phase phi = phi_e + phi_o (even and odd parts under x -> -x) and a real, centro-symmetric pupil amplitude A, the focal field is E = a + i v - y to first order, with a = F{A}, v = F{A phi_e} (real, even) and i y = F{A phi_o} (y real, odd), so the image is

p = a^2 + v^2 + y^2 - 2 a y.

• Its odd part p_o = -2 a y gives the odd phase directly.
• Its even part p_e = a^2 + v^2 + y^2 gives |v| but not its sign.
• The sign comes from the previous image, taken with a known phase difference phi_d (the DM change between the two frames): p_e,prev - p_e = 2 v v_d + v_d^2 + 2 y y_d + y_d^2.

The sign is taken from that difference and the magnitude from the current image alone (Korkiakoski et al. 2014, Eq. 19), which keeps the noise of the two-image difference out of the estimate.

Implemented with the modifications of Korkiakoski et al. (2014):

• images are normalized to the energy of the unaberrated PSF a^2 over the detector window (Eq. 12), and a scaled a^2 is added so the image peak matches max a^2 (Eq. 13), an improved first-order model that accounts for the Strehl loss;
• the division by a is regularized, y = -a p_o / (2 a^2 + epsilon) (Eq. 16; their y has the opposite sign convention);
• |v| = |p_e - a^2 - y^2|^(1/2) (absolute value inside the root, Eq. 19);
• the focal-plane field v + i y is filtered by a concave parabola before the inverse transform (Eq. 21), damping noisy high spatial frequencies.

Each step costs one forward FFT (of the DM change) and one inverse FFT, all on the backend. Units: OPD in metres; images in any linear unit (they are normalized).

The first frame has no predecessor. :class:FastAndFurious then estimates the odd part and, with first_even=True, the even part with positive signs everywhere: a valid even phase whose correction is the even diversity the next frame needs (with only an odd correction the DM change would have no even part and the signs could never be resolved).

FastAndFurious

FastAndFurious(pupil, wavelength, image_shape, *, sampling=None, pixel_scale=None, epsilon=0.0001, filter_radius=None, strehl_compensation=True, first_even=True, device='cpu', precision=None)

Sequential F&F wavefront sensor with its own previous-image state.

Parameters:

Name	Type	Description	Default
pupil	Pupil	The pupil. F&F assumes a centro-symmetric amplitude (circular, annular, symmetric spiders).	required
wavelength	float	Monochromatic wavelength in metres.	required
image_shape	Any	Detector image shape (my, mx); the optical axis is at the window centre.	required
sampling	float | None	Detector pixels per lambda / D or pixel scale in radians (give one). wavelength / (pupil.pitch * pixel_scale) must be an integer FFT size at least as large as the pupil and image grids, e.g. sampling=2 with a pupil filling its grid.	None
pixel_scale	float | None	Detector pixels per lambda / D or pixel scale in radians (give one). wavelength / (pupil.pitch * pixel_scale) must be an integer FFT size at least as large as the pupil and image grids, e.g. sampling=2 with a pupil filling its grid.	None
epsilon	float	Regularization of the division by a (Korkiakoski et al. Eq. 16), relative to the unaberrated peak max a^2 = 1. Korkiakoski et al. recommend 50-500 times the per-pixel noise level of the normalized image. Smaller values recover more of the faint-ring signal (higher spatial frequencies) but amplify noise there.	0.0001
filter_radius	float | None	Radius of the parabolic focal-plane filter max(0, 1 - (r / R)^2) in lambda / D. Default: the largest circle inside the detector window.	None
strehl_compensation	bool	Add the scaled a^2 term of Korkiakoski et al. Eq. (13).	True
first_even	bool	On the first frame, estimate the even part with positive signs (see the module docstring). With False the first estimate is odd only.	True
device	BackendLike	Backend selection.	'cpu'
precision	BackendLike	Backend selection.	'cpu'

Examples:

>>> ff = FastAndFurious(pupil, 1.6e-6, 64, sampling=2)
>>> change = None
>>> for _ in range(20):
...     estimate = ff.step(camera.read(), change)
...     change = -0.5 * estimate  # wavefront change made by the DM
...     dm.apply(change)

sampling property

sampling

Detector pixels per lambda / D.

has_previous property

has_previous

Whether a previous image is stored (the even-part signs can be resolved).

reset

reset()

Forget the previous image (e.g. after opening the loop).

step

step(image, dm_change_opd=None)

Process one image and return the estimated pupil OPD.

Parameters:

Name	Type	Description	Default
image	Any	(my, mx) focal-plane image (any linear unit, background subtracted), NumPy or backend array.	required
dm_change_opd	Any	(ny, nx) change of the wavefront OPD in metres made between the previous image and this one (the DM correction applied after the previous step, as the camera sees it). Ignored on the first frame; None means no change, so the even-part signs are not resolved and the even estimate is zero.	None

Returns:

Type	Description
Backend array ``(ny, nx)``: the estimated wavefront OPD in metres at
this image, zero outside the pupil. Apply ``-gain`` times it with the
DM to correct. The odd and even parts are kept in
attr:`last_odd_opd` and :attr:`last_even_opd`.

ClosedLoopResult dataclass

ClosedLoopResult(residual_rms, strehl, residual_opd, dm_opd, step_times, history=list(), n_iter=0, converged=False, message='', elapsed=0.0, device='cpu')

Record of :func:simulate_closed_loop.

Attributes:

Name	Type	Description
residual_rms	ndarray	(n_iter + 1,) piston-removed, intensity-weighted RMS residual OPD in metres before each frame (entry 0 is the open-loop aberration).
strehl	ndarray	Marechal Strehl ratio of each :attr:residual_rms.
residual_opd	Any	Final residual OPD map in metres (backend array).
dm_opd	Any	Final DM OPD in metres (backend array; the correction it applies).
step_times	ndarray	Wall-clock seconds of each F&F step (synchronized on the GPU).
history	list[float]	:attr:residual_rms as a list.
n_iter	int	Loop iterations run.
converged	bool	Whether the residual fell below tol.
message	str	Why the loop stopped.
elapsed	float	Total wall-clock seconds.
device	str	"cpu" or "gpu".

simulate_closed_loop

simulate_closed_loop(ff, true_opd, n_iterations=20, *, gain=0.5, leak=1.0, photons=None, read_noise=0.0, seed=None, basis=None, tol=None, reset=True)

Run F&F in a simulated closed loop against a static aberration.

Each iteration images the residual true_opd + dm (monochromatic, :class:~solvephase.FocalPlaneModel with the sensor's sampling), adds noise, calls :meth:FastAndFurious.step with the previous DM change, and updates the DM with the leaky integrator of Korkiakoski et al. (Eq. 22), dm <- leak * dm - gain * estimate.

Parameters:

Name	Type	Description	Default
ff	FastAndFurious	The sensor (its state is reset first when reset).	required
true_opd	Any	(ny, nx) static aberration in metres (e.g. NCPA).	required
n_iterations	int	Number of frames.	20
gain	float	Integrator gain and leak (leak = 1 is a pure integrator).	0.5
leak	float	Integrator gain and leak (leak = 1 is a pure integrator).	0.5
photons	float | None	Detected photons per frame (full PSF) for Poisson noise; None gives noise-free images.	None
read_noise	float	Gaussian read noise per pixel in photo-electrons (with photons).	0.0
seed	Any	Seed or :class:numpy.random.Generator for the noise (host draws).	None
basis	Basis | None	Optional DM modal basis: each estimate is projected onto it (a DM cannot make arbitrary shapes). None is a zonal DM at pupil resolution.	None
tol	float | None	Stop once the residual RMS (metres) is below this.	None
reset	bool	Call :meth:FastAndFurious.reset first.	True

Returns:

Type	Description
ClosedLoopResult

solvephase.algorithms.cdi

Coherent diffraction imaging (CDI): iterative projection algorithms.

A compact object is illuminated coherently and only the modulus of its far-field (Fourier) transform is measured. When the diffraction pattern is oversampled (the object occupies less than half of the reconstruction window along each axis, Miao, Sayre & Chapman, JOSA A 15, 1662, 1998) the phase is fixed by the measured modulus together with an object-domain constraint (a support, and optionally realness or positivity). The algorithms here alternate projections onto the two constraint sets:

• P_M - the Fourier-modulus projection: keep the phase of the current transform and impose the measured modulus on measured pixels;
• P_S - the object-domain projection: zero outside the support and, inside it, the nearest complex/real/non-negative value (and amplitude bounds).

Implemented update rules (x is the iterate, y = P_M x, R = 2 P - I):

• ER, error reduction x <- P_S P_M x (Gerchberg & Saxton, Optik 35, 237, 1972; Fienup, Opt. Lett. 3, 27, 1978; Appl. Opt. 21, 2758, 1982);
• HIO, hybrid input-output (Fienup 1982): y where it satisfies the object constraints, x - beta y elsewhere;
• DM, the difference map (Elser, JOSA A 20, 40, 2003) with gamma_S = -1/beta and gamma_M = 1/beta;
• RAAR, relaxed averaged alternating reflections (Luke, Inverse Problems 21, 37, 2005);
• RRR, relax-reflect-reflect (Elser, Lan & Bendory, SIAM J. Imaging Sci. 11, 2429, 2018);
• ASR, averaged successive reflections = Douglas-Rachford (Bauschke, Combettes & Luke, JOSA A 19, 1334, 2002);
• HPR, hybrid projection-reflection (Bauschke, Combettes & Luke, JOSA A 20, 1025, 2003);
• OSS, oversampling smoothness (Rodriguez, Xu, Chen, Zou & Miao, J. Appl. Cryst. 46, 312, 2013): HIO whose out-of-support region is low-pass filtered by a Gaussian of decreasing width;
• shrinkwrap support refinement (Marchesini et al., Phys. Rev. B 68, 140101(R), 2003), usable with any of the above.

Algorithms are chained with a schedule ("hio:400,er:100"), and many random starts run simultaneously as a leading batch dimension, so each iteration is one batched FFT pair. The far-field arrays are centred (DC at n // 2); they are shifted once on entry and once on exit, never inside the loop. Sub-pixel registration in :func:align_object follows Guizar-Sicairos, Thurman & Fienup, Opt. Lett. 33, 156 (2008).

CDI_ALGORITHMS module-attribute

CDI_ALGORITHMS = ('er', 'hio', 'dm', 'raar', 'rrr', 'asr', 'hpr', 'oss')

Names accepted in a :func:cdi schedule.

DEFAULT_BETA module-attribute

DEFAULT_BETA = {'er': 1.0, 'hio': 0.9, 'dm': 0.9, 'raar': 0.98, 'rrr': 0.5, 'asr': 1.0, 'hpr': 0.9, 'oss': 0.9}

Per-algorithm beta used when neither the stage nor :func:cdi sets one.

ER and ASR have no parameter (their entry is informational).

ShrinkwrapConfig dataclass

ShrinkwrapConfig(every=20, threshold=0.2, sigma=3.0, sigma_min=1.5, decay=0.99, start=0, stop=None)

Shrinkwrap support refinement (Marchesini et al. 2003).

Every :attr:every iterations the support becomes G_sigma * |object| > threshold * max(G_sigma * |object|), where G_sigma is a Gaussian of standard deviation sigma pixels and |object| is the modulus of the current data-consistent image P_M x. sigma starts at :attr:sigma and is multiplied by :attr:decay after each update until it reaches :attr:sigma_min.

Attributes:

Name	Type	Description
every	int	Iterations between support updates.
threshold	float	Fraction of the maximum of the blurred modulus kept in the support.
sigma	float	Initial Gaussian standard deviation in pixels.
sigma_min	float	Smallest standard deviation in pixels.
decay	float	Factor applied to sigma after every update (<= 1).
start	int	No update happens before this (global) iteration.
stop	int | None	No update happens after this iteration (None: never stop). Useful to let a final ER stage refine with a fixed support.

CDIResult dataclass

CDIResult(object, support, objects, history, support_history, times, n_iter, converged, message, elapsed, device, best_start, start_errors, modulus_error, support_error, schedule)

Outcome of :func:cdi.

Array fields are backend arrays (CuPy for GPU solves) until you call :meth:to_numpy. All images are centred like the inputs.

Attributes:

Name	Type	Description
object	Any	(ny, nx) complex object estimate of the best start (the support projection of the final data-consistent image).
support	Any	(ny, nx) boolean support of the best start (after shrinkwrap).
objects	Any	(starts, ny, nx) final estimates of every start (for averaging or phase-retrieval transfer functions).
history	list[float]	Fourier-modulus error of the best start at each check; the last entry is the final :attr:modulus_error.
support_history	list[float]	Object-domain constraint violation of the best start at the same points.
times	list[float]	Seconds since the start at each :attr:history entry.
n_iter	int	Iterations performed.
converged	bool	Whether the modulus error reached tol (or the callback stopped).
message	str	Why the run stopped.
elapsed	float	Total wall-clock seconds, including setup.
device	str	"cpu" or "gpu".
best_start	int	Index of the start with the lowest final modulus error.
start_errors	ndarray	Final modulus error of every start (host array).
modulus_error	float	Final normalized Fourier-modulus error of :attr:object, sqrt(sum_meas (|F o| - A)^2 / sum_meas A^2).
support_error	float	Final normalized constraint violation ||u - P_S u|| / ||u|| of the image u :attr:object was projected from.
schedule	list[tuple[str, int, float]]	The parsed schedule [(name, iterations, beta), ...].

to_numpy

to_numpy()

Copy with every array field on the host.

SimulatedCDI dataclass

SimulatedCDI(magnitudes, intensity, object, support)

Synthetic CDI data from :func:simulate_cdi (centred backend arrays).

Attributes:

Name	Type	Description
magnitudes	Any	(ny, nx) far-field modulus: |F o| without noise, or the square root of the Poisson counts.
intensity	Any	(ny, nx) far-field intensity: |F o|^2 or the Poisson counts.
object	Any	(ny, nx) zero-padded true object, scaled so that |F o|^2 is the expected intensity.
support	Any	(ny, nx) tight support |o| > 0 of the padded object.

parse_schedule

parse_schedule(schedule, beta=None)

Normalize a CDI schedule to [(algorithm, iterations, beta), ...].

Parameters:

Name	Type	Description	Default
schedule	ScheduleLike	Either a string "hio:400,er:100" (an optional third field sets the stage's beta: "raar:300:0.8") or a sequence of (name, n) / (name, n, beta) tuples. Names are those in :data:CDI_ALGORITHMS ("dr" is an alias of "asr").	required
beta	float | None	Feedback/relaxation parameter of stages that do not set their own; None uses :data:DEFAULT_BETA.	None

Returns:

Type	Description
list of tuple	One (name, iterations, beta) per non-empty stage.

cdi

cdi(magnitudes, support, *, schedule='hio:500,er:100', starts=1, constraint='complex', bounds=None, shrinkwrap=None, measured=None, intensity=False, beta=None, initial=None, seed=None, device='cpu', precision=None, check_every=10, tol=1e-06, oss_stages=10, callback=None)

Reconstruct an object from its far-field modulus by iterative projections.

Parameters:

Name	Type	Description	Default
magnitudes	Any	(ny, nx) measured far-field modulus |F o|, centred (DC at (ny // 2, nx // 2), i.e. fftshift-ed), any linear units. Non-finite values are treated as unmeasured.	required
support	Any	(ny, nx) boolean object support, centred like the object. Use :func:autocorrelation_support with shrinkwrap when it is unknown.	required
schedule	ScheduleLike	Algorithms to run in order, e.g. "hio:400,er:100" or [("raar", 300, 0.8), ("er", 50)]; see :func:parse_schedule.	'hio:500,er:100'
starts	int	Number of independent random starts run simultaneously as a batch. The start with the lowest final modulus error is returned.	1
constraint	str	Object-domain constraint inside the support: "complex" (support only), "real" or "positive" (real and non-negative).	'complex'
bounds	tuple[float | None, float | None] | None	Optional (min, max) bounds on the object modulus inside the support (either may be None), in object units.	None
shrinkwrap	ShrinkwrapLike	Support refinement: None/False (fixed support), True (defaults), a dict of :class:ShrinkwrapConfig fields, or a config. Refined supports stay inside the initial support (unless it is all True), so pass a generous one, e.g. :func:autocorrelation_support.	None
measured	Any	Optional (ny, nx) boolean mask, centred, False where the modulus is unknown (beamstop, detector gaps); the modulus is left free there.	None
intensity	bool	If True, magnitudes holds intensities |F o|^2 (negative values are clipped to zero before the square root).	False
beta	float | None	Feedback (HIO, HPR, OSS) / relaxation (DM, RAAR, RRR) parameter of stages that do not set their own. None uses the per-algorithm :data:DEFAULT_BETA (0.9; RAAR 0.98; RRR 0.5). ER and ASR ignore it.	None
initial	Any	Optional starting object, (ny, nx) (shared by all starts) or (starts, ny, nx), centred. Default: random Fourier phases with the measured modulus, restricted to the support.	None
seed	Any	Seed or :class:numpy.random.Generator for the random starts (drawn on the host, so CPU and GPU runs start identically).	None
device	BackendLike	Backend selection, see :func:solvephase.get_backend.	'cpu'
precision	BackendLike	Backend selection, see :func:solvephase.get_backend.	'cpu'
check_every	int	Iterations between error evaluations (each one synchronizes the GPU).	10
tol	float	Stop when the best start's modulus error is at or below tol.	1e-06
oss_stages	int	Number of filter widths an OSS stage steps through (its iterations are split evenly); at the end of each the best iterate is kept.	10
callback	Callable[[int, Any, float], bool | None] | None	callback(iteration, object, error) at every check with the centred estimate of the currently best start; return True to stop.	None

Returns:

Type	Description
CDIResult	The best start's object and support, error histories and per-start final errors.

Notes

Each iteration of ER, HIO, RAAR, RRR, ASR and HPR costs one batched FFT pair; DM and OSS cost two. An error check adds one forward FFT.

autocorrelation_support

autocorrelation_support(intensity, threshold=0.04, *, measured=None, device=None)

Loose initial support from the thresholded autocorrelation.

The inverse Fourier transform of the far-field intensity is the object's autocorrelation, whose support is twice the object's extent (the difference set S - S); thresholded, it is the usual starting support for shrinkwrap.

Parameters:

Name	Type	Description	Default
intensity	Any	(ny, nx) centred far-field intensity |F o|^2 (square magnitudes first). Non-finite values count as unmeasured.	required
threshold	float	Fraction of the autocorrelation maximum kept.	0.04
measured	Any	Optional centred boolean mask, False where the intensity is unknown (set to zero before transforming).	None
device	BackendLike	Backend; by default the backend of intensity.	None

Returns:

Type	Description
array	(ny, nx) centred boolean support on the backend.

align_object

align_object(estimate, reference, *, subpixel=True, upsample=100, twin=True, scale=True)

Register a CDI reconstruction to a reference, removing the trivial ambiguities.

A far-field modulus does not change under a global phase factor, a translation, or the twin image conj(o(-r)). This finds the translation (integer cross-correlation peak refined to 1/upsample pixel by an upsampled matrix DFT, Guizar-Sicairos et al. 2008), the twin choice and the complex factor c that minimize ||c * shift(est) - ref|| / ||ref||.

Parameters:

Name	Type	Description	Default
estimate	Any	(ny, nx) complex or real images of the same shape.	required
reference	Any	(ny, nx) complex or real images of the same shape.	required
subpixel	bool	Refine the translation below one pixel (applied as a Fourier phase ramp, i.e. a circular shift).	True
upsample	int	Sub-pixel refinement factor.	100
twin	bool	Also try the twin conj(est[::-1, ::-1]) and keep the better one.	True
scale	bool	If True c is any complex number (phase and scale); if False only a global phase |c| = 1 is removed.	True

Returns:

Name	Type	Description
aligned	Any	c * shift(est) (or its twin), on the backend of estimate.
error	float	||aligned - ref|| / ||ref||.

simulate_cdi

simulate_cdi(obj, oversampling=2.0, *, shape=None, photons=None, seed=None, device='cpu', precision=None)

Zero-pad an object and compute its centred far-field diffraction data.

Parameters:

Name	Type	Description	Default
obj	Any	(my, mx) complex or real object.	required
oversampling	float	Linear oversampling per axis: the window is round(oversampling * m) pixels along each axis. Unique recovery needs Miao's oversampling ratio (window area / object area) above 2, i.e. oversampling above sqrt(2) for an object filling its box.	2.0
shape	tuple[int, int] | None	Explicit (ny, nx) window instead of oversampling.	None
photons	float | None	If given, the expected total photon count: the intensity is scaled to it and Poisson noise is drawn (on the host).	None
seed	Any	Seed or :class:numpy.random.Generator for the noise.	None
device	BackendLike	Backend of the returned arrays.	'cpu'
precision	BackendLike	Backend of the returned arrays.	'cpu'

Returns:

Type	Description
SimulatedCDI

solvephase.algorithms.wirtinger

Phase retrieval for generic linear measurement models: the Wirtinger-flow family.

Recover x (up to a global phase) from y = |A x|^2 for a known linear operator A (:mod:solvephase.operators). The algorithms are the gradient methods catalogued by PhasePack (Chandra et al., Asilomar 2017), each started from a spectral-type initialization:

Initializations * spectral - leading eigenvector of A^H diag(y) A (Netrapalli, Jain & Sanghavi, NeurIPS 2013; Candès, Li & Soltanolkotabi, IEEE Trans. Inf. Theory 61, 1985, 2015); * truncated spectral - drops measurements above alpha_y^2 mean(y) (Chen & Candès, Comm. Pure Appl. Math. 70, 822, 2017); * orthogonality-promoting - the m/6 largest amplitudes, unweighted (Wang, Giannakis & Eldar, IEEE Trans. Inf. Theory 64, 773, 2018); * weighted maximal correlation - the 3m/13 largest amplitudes weighted by psi^(1/2) (Wang, Giannakis, Saad & Chen, IEEE Trans. Signal Process. 66, 2818, 2018); * optimal preprocessing - T(y) = (y - 1) / (y + sqrt(delta) - 1) with y normalized to unit mean and delta = m / n (Luo, Alghamdi & Lu, IEEE Trans. Signal Process. 67, 2347, 2019); * random.

Algorithms * "wf" - Wirtinger Flow on the intensity loss with the step schedule mu_t = min(1 - exp(-t / tau0), mu_max) (Candès et al. 2015; tau0 = 330; mu_max defaults to 0.1 because the paper's 0.2 diverges or stalls on many Gaussian instances with n >= 128); * "twf" - Truncated Wirtinger Flow on the Poisson likelihood (Chen & Candès 2017); * "taf" - Truncated Amplitude Flow (Wang, Giannakis & Eldar 2018); * "raf" - Reweighted Amplitude Flow (Wang, Giannakis, Saad & Chen 2018); * "lbfgs" - L-BFGS (:func:solvephase.optimize.lbfgs) on the amplitude or intensity least-squares loss, using a real view of the complex vector.

Scaling. The literature states every rule for the unit-variance complex Gaussian model (E |A_kj|^2 = 1). Internally the operator is rescaled to that convention, A -> sqrt(s) A and y -> s y with s = m n / ||A||_F^2, so the published step sizes and truncation thresholds apply unchanged to any operator (for coded diffraction with a unitary FFT this reproduces the unnormalized-FFT model of the CDP papers). Measurement vectors are assumed to have nearly equal norms, which holds for the Gaussian and CDP models; the row-norm normalizations of the TAF/RAF initializations are therefore dropped.

Whitening. When the operator reports G = diag(A^H A) (:meth:~solvephase.operators.LinearOperator.gram_diagonal; every built-in operator does), the power method runs on G^{-1/2} Y G^{-1/2} and maps the eigenvector back with G^{1/2}. For the Gaussian model G ~ m I and nothing changes; for coded diffraction with octanary masks it removes eigenvectors localized on pixels where most masks are strong, which otherwise beat the signal eigenvector for images of about 10^4 pixels or more.

All iterations stay on the backend device; the residual is transferred to the host only every check_every iterations.

GenericResult dataclass

GenericResult(x, history, times, n_iter, converged, message, elapsed, device, method, init, options=dict())

Outcome of :func:wirtinger.

Attributes:

Name	Type	Description
x	Any	Recovered vector, shape operator.in_shape (backend array), in the units of the true signal (the global phase is arbitrary).
history	list[float]	Relative amplitude residual || |A x| - sqrt(y) || / || sqrt(y) || at the start and at every check.
times	list[float]	Wall-clock seconds since the call started, for each history entry (the first entry includes the initialization).
n_iter	int	Iterations performed.
converged	bool	Whether a tolerance (rather than the iteration limit) stopped the run.
message	str	Why the run stopped.
elapsed	float	Total wall-clock seconds, including the initialization.
device	str	"cpu" or "gpu".
method	str	Algorithm used.
init	str	Initialization used ("given" for a user-supplied start).
options	dict[str, float]	Method parameters used (step and the method's options).

residual property

residual

Final relative amplitude residual.

to_numpy

to_numpy()

A copy with :attr:x on the host as a NumPy array.

relative_error

relative_error(x, truth)

Distance to the truth modulo the global phase.

min_phi ||x - exp(i phi) truth|| / ||truth||, with the optimal phase exp(i phi) = truth^H x / |truth^H x|. Works for NumPy or CuPy input (computed on the device of x).

init_spectral

init_spectral(operator, measurements, *, magnitudes=False, iterations=100, tol=1e-06, seed=None)

Spectral initialization (Netrapalli et al. 2013; Candès et al. 2015).

Leading eigenvector of Y = (1/m) A^H diag(y) A by the power method, scaled to the norm estimate ||x||^2 ~ mean(y) (in the unit-variance convention; see the module docstring).

Parameters:

Name	Type	Description	Default
operator	LinearOperator	The measurement operator A.	required
measurements	Any	Intensities y = |A x|^2 (or magnitudes |A x| with magnitudes=True), shape operator.out_shape.	required
magnitudes	bool	Whether measurements are magnitudes rather than intensities.	False
iterations	int	Maximum power iterations.	100
tol	float	Relative change of the Rayleigh quotient that stops the power method.	1e-06
seed	Any	Seed of the random power-method start (host RNG).	None

Returns:

Type	Description
array	Initial estimate of shape operator.in_shape on the operator's backend.

init_truncated_spectral

init_truncated_spectral(operator, measurements, *, magnitudes=False, alpha_y=3.0, iterations=100, tol=1e-06, seed=None)

Truncated spectral initialization (Chen & Candès 2017).

As :func:init_spectral but keeping only y_k <= alpha_y^2 mean(y), which removes the heavy tail that biases the plain spectral estimate. Parameters are those of :func:init_spectral plus alpha_y (dimensionless, default 3 as in the paper).

init_orthogonality_promoting

init_orthogonality_promoting(operator, measurements, *, magnitudes=False, fraction=1.0 / 6.0, iterations=100, tol=1e-06, seed=None)

Orthogonality-promoting initialization (Wang, Giannakis & Eldar 2018).

Leading eigenvector of sum_{k in S} a_k a_k^* over the ceil(fraction m) largest amplitudes (the measurement vectors most aligned with x), scaled to sqrt(mean(y)). Parameters are those of :func:init_spectral plus fraction (default 1/6, as in TAF).

init_weighted_correlation

init_weighted_correlation(operator, measurements, *, magnitudes=False, fraction=3.0 / 13.0, exponent=0.5, iterations=100, tol=1e-06, seed=None)

Weighted maximal correlation initialization (Wang, Giannakis, Saad & Chen 2018).

Leading eigenvector of sum_{k in S} psi_k^exponent a_k a_k^* over the ceil(fraction m) largest amplitudes psi = sqrt(y), scaled to sqrt(mean(y)). Defaults (3/13 and 1/2) follow the RAF paper.

init_optimal_spectral

init_optimal_spectral(operator, measurements, *, magnitudes=False, iterations=200, tol=1e-06, seed=None)

Spectral initialization with optimal preprocessing (Luo, Alghamdi & Lu 2019).

Leading (largest algebraic) eigenvector of (1/m) A^H diag(T(y)) A with T(y) = (y - 1) / (y + sqrt(delta) - 1), y normalized to unit mean and delta = m / n > 1. T is negative for small y, so the power method runs on the matrix shifted by ||A||^2 / (m (sqrt(delta) - 1)) (||A|| from :meth:~solvephase.operators.LinearOperator.norm_estimate), which slows it; the default iterations is therefore larger. Asymptotically this weighting minimizes the number of measurements needed for a non-trivial correlation with x.

init_random

init_random(operator, measurements, *, magnitudes=False, seed=None)

Random complex Gaussian start scaled to the norm estimate sqrt(mean(y)).

wirtinger

wirtinger(operator, measurements, *, method='raf', init='auto', magnitudes=False, iterations=1000, tol=1e-07, check_every=10, step=None, loss='amplitude', options=None, power_iterations=None, seed=None, callback=None)

Recover x from y = |A x|^2 with a Wirtinger-flow-family solver.

Parameters:

Name	Type	Description	Default
operator	LinearOperator	Measurement operator A (:mod:solvephase.operators); its backend sets the device and precision.	required
measurements	Any	Intensities y = |A x|^2, shape operator.out_shape (host or device array). Negative values (noise) are clipped to zero.	required
method	str	"raf" (Reweighted Amplitude Flow, default), "taf" (Truncated Amplitude Flow), "twf" (Truncated Wirtinger Flow), "wf" (Wirtinger Flow) or "lbfgs" (L-BFGS on the loss).	'raf'
init	Any	"auto" (the method's own initialization, see :data:DEFAULT_INIT), "spectral", "truncated", "orthogonal", "weighted", "optimal", "random", or an array of shape operator.in_shape to start from.	'auto'
magnitudes	bool	Whether measurements are magnitudes |A x| instead of intensities.	False
iterations	int	Maximum iterations (each costs one forward and one adjoint; an L-BFGS iteration may cost more through its line search).	1000
tol	float	Stop when the relative amplitude residual falls below tol (exact fit), or changes by less than tol (relative) between checks (stationary point, e.g. with noisy data). It is raised to 8 eps of the working precision (about 1e-6 in single precision).	1e-07
check_every	int	Iterations between convergence checks (the only host syncs).	10
step	float | None	Step size mu overriding the method default (mu_max for WF).	None
loss	str	"amplitude" or "intensity" least squares, for method="lbfgs".	'amplitude'
options	dict[str, float] | None	Method parameters overriding :data:DEFAULT_OPTIONS: tau0 (WF); alpha_lb, alpha_ub, alpha_h (TWF); gamma (TAF); beta (RAF); memory (L-BFGS).	None
power_iterations	int | None	Power-method iterations of the initialization (default 100, 200 for "optimal").	None
seed	Any	Seed for the initialization (host RNG; CPU and GPU agree).	None
callback	Callback | None	callback(iteration, x, residual) at each check; returning True stops the run.	None

Returns:

Type	Description
GenericResult	The estimate (global phase arbitrary; compare with :func:relative_error) and the residual history.

solvephase.operators

Generic linear measurement operators with exact adjoints.

The generic phase-retrieval solvers (:mod:solvephase.algorithms.wirtinger) recover x from y = |A x|^2 for any linear A. This module supplies the operators:

• :class:MatrixOperator - a dense complex matrix, with :meth:MatrixOperator.gaussian for the i.i.d. complex Gaussian model used in the Wirtinger-flow literature (Candès, Li & Soltanolkotabi, IEEE Trans. Inf. Theory 61, 1985, 2015);
• :class:CodedDiffractionOperator - coded diffraction patterns (CDP): L random modulation masks followed by a unitary FFT (Candès, Li & Soltanolkotabi, Appl. Comput. Harmon. Anal. 39, 277, 2015);
• :class:OversampledFourierOperator - zero padding followed by a unitary FFT (the classical Fourier phase-retrieval model).

Every operator maps (..., *in_shape) to (..., *out_shape) (leading batch axes are allowed), works in the precision of its :class:Backend, and ships :meth:~LinearOperator.adjoint, its exact conjugate transpose.

LinearOperator

Bases: ABC

A linear map A between complex arrays, with its exact adjoint.

Subclasses set :attr:backend, :attr:in_shape and :attr:out_shape and implement :meth:forward and :meth:adjoint.

in_size property

in_size

Number of unknowns n (complex entries of the input).

out_size property

out_size

Number of measurements m (entries of the output).

forward abstractmethod

forward(x)

Apply A: (..., *in_shape) to (..., *out_shape).

adjoint abstractmethod

adjoint(y)

Apply A^H: (..., *out_shape) to (..., *in_shape).

frobenius_norm_squared

frobenius_norm_squared(*, probes=32, seed=0)

trace(A^H A) = ||A||_F^2.

||A||_F^2 / n is the mean energy gain E ||A x||^2 / ||x||^2 for a random x; the solvers use it to put every operator on the scale of the unit-variance Gaussian model. The base implementation is a Hutchinson estimate from probes random phase vectors (relative error about 1 / sqrt(probes) at worst); concrete operators override it with the exact value.

gram_diagonal

gram_diagonal()

diag(A^H A) as a real array of shape :attr:in_shape, or None.

The squared norm of each column of A. Spectral initializations use it to whiten the power method (see :mod:solvephase.algorithms.wirtinger). The base class returns None (unknown); every built-in operator returns the exact value.

norm_estimate

norm_estimate(*, iterations=100, tol=1e-06, seed=0, check_every=5)

Spectral norm ||A||_2 by power iteration on A^H A.

Parameters:

Name	Type	Description	Default
iterations	int	Maximum power iterations.	100
tol	float	Stop when the estimate changes by less than tol (relative) between checks.	1e-06
seed	Any	Seed of the random start (host RNG, so CPU and GPU agree).	0
check_every	int	Iterations between scalar convergence checks (device syncs).	5

Returns:

Type	Description
float	The estimate, which approaches ||A||_2 from below.

MatrixOperator

MatrixOperator(matrix, *, backend=None)

Bases: LinearOperator

A dense complex matrix A of shape (m, n): y = A x.

Row k of A is the conjugated measurement vector a_k^*, so y_k = |a_k^* x|^2. forward accepts (..., n) and returns (..., m).

Parameters:

Name	Type	Description	Default
matrix	Any	(m, n) array (host or device); stored in the backend's complex precision.	required
backend	BackendLike	Target backend (default CPU, double precision).	None

gaussian classmethod

gaussian(m, n, *, seed=None, backend=None)

i.i.d. circularly-symmetric complex Gaussian matrix, E |A_kj|^2 = 1.

Entries are (g1 + i g2) / sqrt(2) with g1, g2 ~ N(0, 1), drawn on the host from seed (Candès et al. 2015: a_k ~ N(0, I/2) + i N(0, I/2), so E y_k = ||x||^2).

frobenius_norm_squared

frobenius_norm_squared(*, probes=32, seed=0)

Exact ||A||_F^2 (probes and seed are ignored).

gram_diagonal

gram_diagonal()

Squared column norms sum_k |A_kj|^2.

CodedDiffractionOperator

CodedDiffractionOperator(masks, *, backend=None)

Bases: LinearOperator

Coded diffraction patterns: y_l = |F(d_l * x)|^2 for l < L.

Each of the L masks d_l multiplies the (ny, nx) image before a unitary 2-D FFT F (batched over the masks). forward maps (..., ny, nx) to (..., L, ny, nx); the adjoint is sum_l conj(d_l) * F^H y_l. Candès, Li & Soltanolkotabi (2015) prove recovery from L = O(log^4 n) octanary masks; L = 6-8 suffices in practice.

Parameters:

Name	Type	Description	Default
masks	Any	(L, ny, nx) complex mask array.	required
backend	BackendLike	Target backend.	None

random classmethod

random(shape, n_masks, *, kind='octanary', seed=None, backend=None)

Random masks drawn on the host from seed.

Parameters:

Name	Type	Description	Default
shape	Any	Image shape (ny, nx) (or an int for a square image).	required
n_masks	int	Number of masks L.	required
kind	str	"octanary" - d = b1 * b2 with b1 uniform on {1, -1, i, -i} and b2 = 1/sqrt(2) with probability 4/5, sqrt(3) with probability 1/5 (Candès et al. 2015, eq. 1.9; E |d|^2 = 1); "uniform" - unit-modulus masks with phases uniform on [0, 2 pi).	'octanary'

frobenius_norm_squared

frobenius_norm_squared(*, probes=32, seed=0)

Exact ||A||_F^2 = sum |d|^2 (probes and seed are ignored).

gram_diagonal

gram_diagonal()

sum_l |d_l|^2: A^H A is diagonal for a unitary FFT.

OversampledFourierOperator

OversampledFourierOperator(in_shape, oversampling=2.0, *, out_shape=None, backend=None)

Bases: LinearOperator

Zero padding followed by a unitary 2-D FFT.

The in_shape image is placed in the corner of an out_shape grid and transformed, so A^H A = I (an isometry) and A^H crops the inverse FFT. forward maps (..., *in_shape) to (..., *out_shape).

Fourier magnitudes alone do not determine an image up to a global phase: shifts and the conjugate flip x(-r)* give the same data. The Wirtinger-flow solvers therefore rarely succeed with this operator by itself; use it with support/positivity solvers or add masks (:class:CodedDiffractionOperator).

Parameters:

Name	Type	Description	Default
in_shape	Any	Image shape (ny, nx).	required
oversampling	Any	Output size per input size (scalar or (y, x) pair, >= 1); ignored when out_shape is given.	2.0
out_shape	Any	Explicit FFT grid shape, at least in_shape.	None
backend	BackendLike	Target backend.	None

frobenius_norm_squared

frobenius_norm_squared(*, probes=32, seed=0)

Exact ||A||_F^2 = n (an isometry; arguments are ignored).

gram_diagonal

gram_diagonal()

All ones (A^H A = I).

solvephase.algorithms.tie

Transport-of-intensity (TIE) phase retrieval from a defocus stack.

In the paraxial regime the axial change of intensity of a monochromatic beam is fixed by its transverse phase (Teague, JOSA 73, 1434, 1983)::

-k dI/dz = div(I grad(phi)),    k = 2 pi / lambda.

Given intensities I(z_j) near the plane z = 0 the solver estimates the axial derivative dI/dz and the in-focus intensity I0, then inverts the elliptic equation for phi (radians at z = 0):

• Axial derivative. "central" uses the nearest planes on each side of z = 0 (and z = 0 itself when measured): a central difference for a symmetric pair, the three-point Lagrange derivative otherwise. "polyfit" fits a polynomial of degree order in z to every plane by least squares, pixel by pixel (Savitzky-Golay style), and takes its slope at z = 0; with more planes than coefficients it averages noise, with order = N - 1 it is the exact N-point stencil whose truncation error falls as a higher power of the plane spacing (Waller et al., Opt. Express 18, 12552, 2010). The measured z = 0 plane is the in-focus intensity when present; otherwise the interpolant (or fit) value at z = 0 is used.
• Spectral solvers. "fft" assumes a periodic field (real FFTs); "dct" assumes zero normal slope at the window edge and expands the field in the cosine series of its mirror extension (type-II DCT; gradients are type-II DST series), which avoids the wrap-around coupling of opposite edges for non-periodic fields (Zuo, Chen & Asundi, Opt. Express 22, 9220, 2014). Both differentiate exactly and use either the uniform-intensity approximation laplacian(phi) = -k dI/dz / mean(I0) (uniform=True) or Teague's auxiliary function psi with grad(psi) = I grad(phi): two Poisson solves and a division by the clipped intensity (Gureyev & Nugent, JOSA A 13, 1670, 1996; Paganin & Nugent, PRL 80, 2586, 1998).
• Masked exact solver. Teague's construction assumes I grad(phi) is curl free. "pcg" drops that assumption and solves the five-point discretization of div(I grad(phi)) = -k dI/dz inside a support mask (default: where I0 exceeds intensity_floor of its maximum) by conjugate gradients preconditioned with the DCT Poisson solver. Edges that leave the mask carry no flux, so it is the right solver for an illuminated aperture with dark surroundings, where the edge signal carries the boundary slopes (curvature sensing, Roddier, Appl. Opt. 27, 1223, 1988).
• Regularization. The first Poisson solve of "fft"/"dct" uses the Tikhonov filter lam / (lam**2 + alpha) on the eigenvalues lam of -laplacian with alpha = regularization * lam_min**2: the lowest non-zero spatial frequency of the grid is attenuated by exactly 1 / (1 + regularization). Raise it to suppress the low-frequency "cloud" artefacts that noise produces (Zuo et al., Opt. Lasers Eng. 135, 106187, 2020). Teague's second solve is bounded and left unregularized.

Piston is not measurable and is removed (mean over the mask, or over each connected region of the mask for "pcg").

TIEResult dataclass

TIEResult(phase, opd, intensity, didz, mask, wavelength, method, derivative, uniform, history=list(), n_iter=0, converged=True, message='', elapsed=0.0, device='cpu')

Outcome of a transport-of-intensity retrieval.

Array fields are backend arrays (CuPy for GPU solves) until you call :meth:to_numpy.

Attributes:

Name	Type	Description
phase	Any	(ny, nx) phase at z = 0 in radians, piston removed, zero outside mask.
opd	Any	The same wavefront as optical path difference in metres (phase * wavelength / (2 pi)).
intensity	Any	(ny, nx) in-focus intensity I0 used by the solver (measured or estimated), in data units.
didz	Any	(ny, nx) axial intensity derivative estimate, data units per metre.
mask	Any	Boolean support the phase is defined on.
wavelength	float	Wavelength in metres.
method	str	Solver: "fft", "dct" or "pcg".
derivative	str	Axial derivative estimator: "central" or "polyfit".
uniform	bool	Whether the uniform-intensity approximation was used.
history	list[float]	Relative residual norm per convergence check ("pcg" only).
n_iter	int	Conjugate-gradient iterations (0 for the direct solvers).
converged	bool	Whether the solve finished (always true for direct solvers).
message	str	Why it stopped.
elapsed	float	Wall-clock seconds, including the derivative estimate.
device	str	"cpu" or "gpu".

to_numpy

to_numpy()

Copy with every array field on the host.

simulate_defocus_stack

simulate_defocus_stack(field, distances, pitch, wavelength, *, paraxial=False, device=None, precision=None)

Intensities of a complex field propagated to several distances.

Uses :class:~solvephase.propagation.AngularSpectrumPropagator (periodic boundaries: pad the field if light must not wrap around).

Parameters:

Name	Type	Description	Default
field	Any	(ny, nx) complex field at z = 0 (amplitude, not intensity).	required
distances	float | Sequence[float]	Propagation distances in metres (positive downstream); 0 returns |field|**2.	required
pitch	float	Sample pitch in metres.	required
wavelength	float	Wavelength in metres.	required
paraxial	bool	Use the Fresnel transfer function instead of the exact angular spectrum.	False
device	BackendLike	Backend; defaults to the backend of field.	None
precision	BackendLike	Backend; defaults to the backend of field.	None

Returns:

Type	Description
``(Z, ny, nx)`` intensity stack on the backend (``(ny, nx)`` for a scalar
distance).

tie

tie(intensities, distances, *, pitch, wavelength, method='dct', uniform=False, derivative='central', order=None, in_focus=None, regularization=1e-06, intensity_floor=0.001, mask=None, tol=1e-06, max_iter=2000, check_every=10, device=None, precision=None)

Retrieve the phase at z = 0 from intensities at known defocus.

Solves -k dI/dz = div(I grad phi); see the module notes and the transport-of-intensity guide page.

Parameters:

Name	Type	Description	Default
intensities	Any	(N, ny, nx) intensity images (NumPy or CuPy), any consistent units, N >= 2.	required
distances	Sequence[float] | Any	(N,) defocus of each image in metres, relative to the plane where the phase is wanted (positive downstream). May include 0.	required
pitch	float | tuple[float, float]	Pixel pitch in metres, scalar or (y, x).	required
wavelength	float	Wavelength in metres.	required
method	TIEMethod	"dct" (Neumann boundaries, default), "fft" (periodic) or "pcg" (exact non-uniform solve inside mask, Neumann).	'dct'
uniform	bool	Use the uniform-intensity approximation I = mean(I0) (one Poisson solve) instead of the non-uniform solution.	False
derivative	TIEDerivative	"central" (nearest planes about z = 0) or "polyfit" (least-squares polynomial over all planes).	'central'
order	int | None	Polynomial degree for "polyfit"; default min(N - 1, 3).	None
in_focus	Any	Optional (ny, nx) in-focus intensity; overrides the measured z = 0 plane and the estimate.	None
regularization	float	Dimensionless Tikhonov weight for "fft"/"dct": the lowest non-zero frequency is attenuated by 1 / (1 + regularization). Ignored by "pcg".	1e-06
intensity_floor	float	Fraction of max(I0) below which the intensity is clipped: when dividing by it ("fft"/"dct") and in the "pcg" weights, whose condition number it bounds. The default "pcg" mask is where I0 exceeds it.	0.001
mask	Any	Optional (ny, nx) boolean support. The phase is zero outside it and its piston is removed over it. For "pcg" the equation is solved only inside it.	None
tol	float	"pcg" stopping rule: relative residual, iteration limit, and how often the residual is checked (each check syncs a GPU).	1e-06
max_iter	float	"pcg" stopping rule: relative residual, iteration limit, and how often the residual is checked (each check syncs a GPU).	1e-06
check_every	float	"pcg" stopping rule: relative residual, iteration limit, and how often the residual is checked (each check syncs a GPU).	1e-06
device	BackendLike	Backend; defaults to the backend of intensities.	None
precision	BackendLike	Backend; defaults to the backend of intensities.	None

Returns:

Type	Description
TIEResult	Phase (radians) and OPD (metres) on the backend.

Losses and optimizers

solvephase.losses

Data-fidelity terms comparing model images to measured images.

Each loss returns its value, its gradient with respect to the model, and a non-negative per-pixel curvature (the expected Hessian, i.e. Fisher information) used by Gauss-Newton/Levenberg-Marquardt. All operations run on the backend of the input arrays; values are 0-d backend arrays so the caller decides when to synchronize.

weights (same shape as the data, or broadcastable) scales each pixel's contribution; zero excludes a pixel (bad pixels, saturation, beyond a field stop).

Loss

Bases: ABC

A data-fidelity term L(model, data).

evaluate abstractmethod

evaluate(backend, model, data, weights)

Return (value, d value / d model, curvature) on the backend.

GaussianLoss

Bases: Loss

Weighted least squares 0.5 * sum w (model - data)^2.

With weights = 1 / variance this is the Gaussian negative log-likelihood. :func:noise_weights gives a photon + read-noise variance.

PoissonLoss

PoissonLoss(read_noise=0.0, floor=None)

Bases: Loss

Poisson negative log-likelihood (Paxman, Schulz & Fienup 1992).

Evaluated as the deviance sum w [ (m' - d') + d' log(d' / m') ] with m' = m + s, d' = d + s and s = read_noise^2 (in the data's units squared), the shifted-Poisson approximation that also covers Gaussian read noise. The deviance differs from the negative log-likelihood by a data-only constant, is non-negative, and is about chi^2 / 2 at a good fit.

Below a small positive floor m0 (default 1e-6 max|d|) the loss continues as its second-order Taylor expansion about m0, so value and gradient stay finite and consistent when a line search probes negative model values.

Parameters:

Name	Type	Description	Default
read_noise	float	Read-noise standard deviation in data units (e.g. photo-electrons).	0.0
floor	float | None	Explicit model floor m0 in data units.	None

AmplitudeLoss

AmplitudeLoss(floor=None)

Bases: Loss

Least squares on square-root intensities, 0.5 * sum w (sqrt m - sqrt d)^2.

The amplitude metric of Fienup's iterative transform algorithms; it variance-stabilizes photon noise and converges fast from poor starts. Below a small positive floor (default 1e-6 max|d|) it continues as its second-order Taylor expansion, like :class:PoissonLoss.

make_loss

make_loss(loss, **kwargs)

Return a :class:Loss from a name ("gaussian", "poisson", "amplitude").

solvephase.optimize

Device-resident optimizers.

Parameters stay on the backend (CPU or GPU) for the whole solve; only scalars (objective values and directional derivatives) cross to the host for line search decisions. The objective is a callable fun(x) -> (f, g) returning a Python float and a backend gradient array.

• :func:lbfgs - limited-memory BFGS with a strong-Wolfe line search (Nocedal & Wright, Numerical Optimization, 2nd ed., Alg. 7.4, 3.5, 3.6).
• :func:levenberg_marquardt - damped Gauss-Newton for problems with few parameters, with Nielsen's damping update.
• :func:adam - first-order Adam (Kingma & Ba 2015), for very large or noisy objectives.

OptimizeResult dataclass

OptimizeResult(x, fun, n_iter, n_eval, converged, message, history=list(), times=list())

Outcome of an optimizer run.

Attributes:

Name	Type	Description
x	Any	Final parameter vector (backend array).
fun	float	Final objective value.
n_iter	int	Iterations taken.
n_eval	int	Objective evaluations.
converged	bool	Whether a convergence test (rather than the iteration limit) stopped the run.
message	str	Why the run stopped.
history	list[float]	Objective value after each iteration (the first entry is the start).
times	list[float]	Wall-clock seconds since the start at each history entry.

lbfgs

lbfgs(fun, x0, backend, *, max_iter=200, memory=10, ftol=1e-10, gtol=1e-10, max_eval=None, callback=None)

Minimize with limited-memory BFGS and a strong-Wolfe line search.

Parameters:

Name	Type	Description	Default
fun	Objective	fun(x) -> (f, g) with f a Python float and g a backend array.	required
x0	Any	Starting point (backend array, any shape; treated as a flat vector).	required
max_iter	int	Iteration limit.	200
memory	int	Number of correction pairs kept.	10
ftol	float	Stop when the relative decrease of f over an iteration is below ftol.	1e-10
gtol	float	Stop when max |g| falls below gtol times its initial value.	1e-10
max_eval	int | None	Limit on objective evaluations (default 10 * max_iter).	None
callback	Callback | None	callback(iteration, x, f); returning True stops the run.	None

levenberg_marquardt

levenberg_marquardt(system, value, x0, backend, *, max_iter=50, ftol=1e-10, xtol=1e-10, damping=0.001, callback=None)

Damped Gauss-Newton (Levenberg-Marquardt) for small parameter counts.

Parameters:

Name	Type	Description	Default
system	Callable[[Any], tuple[float, Any, Any]]	system(x) -> (f, g, H): objective, gradient and a positive semi-definite Hessian approximation (Gauss-Newton/Fisher), as backend arrays except f.	required
value	Callable[[Any], float]	value(x) -> f for trial steps.	required
damping	float	Initial damping relative to the largest diagonal of H.	0.001

Notes

The step solves (H + mu diag(H)) dx = -g (Marquardt scaling) on the host in double precision; mu follows Nielsen's update.

adam

adam(fun, x0, backend, *, max_iter=500, learning_rate=0.05, beta1=0.9, beta2=0.999, eps=1e-08, ftol=0.0, callback=None)

First-order Adam with bias correction.

learning_rate is in parameter units per step (internal parameters are radians of phase, so ~0.01-0.1 is typical).

Shared primitives (aocore)

The backend, pupils, propagators, wavefront metrics and phase unwrapping come from aocore, the stack's shared core. solvephase re-exports them unchanged (solvephase.Pupil, solvephase.FocalPlanePropagator, solvephase.rms, ...).

aocore.propagation

Linear optical propagators with exact adjoints.

All propagators map a complex field of shape (..., ny, nx) to (..., my, mx) and provide :meth:~Propagator.adjoint, the exact conjugate transpose of :meth:~Propagator.forward. The adjoint is what makes analytic gradients (and therefore fast optimization) possible.

Conventions

• Arrays are (y, x). Coordinates sit on pixel centres and are centred: pixel i of n is at (i - (n - 1) / 2) * pitch.
• Fraunhofer propagation uses the kernel exp(-2 pi i x . alpha / lambda), so a pupil OPD ramp OPD = a * x moves the image to angle +a (towards increasing column index). See CONVENTIONS.md, rule 3.1.
• Normalization is unitary: when the output window captures all the light, sum |E|^2 == sum |u|^2 (Parseval). A finite window loses light; it is never renormalized.

Propagator

Bases: ABC

A linear map between two sampled planes with an exact adjoint.

forward abstractmethod

forward(field)

Propagate (..., *in_shape) to (..., *out_shape).

adjoint abstractmethod

adjoint(field)

Apply the conjugate transpose: (..., *out_shape) to (..., *in_shape).

FFTPropagator

FFTPropagator(in_shape, out_shape, n_fft, *, offset=0.0, backend=None)

Bases: Propagator

Zero-padded FFT between centred grids.

Equivalent to :class:MFTPropagator with samples = n_fft but costs one n_fft-sized FFT. The output window of out_shape pixels is centred on the optical axis (shifted by offset output pixels); fractional offsets are exact, applied as an input phase ramp.

Parameters:

Name	Type	Description	Default
in_shape	Any	Input (pupil) grid shape.	required
out_shape	Any	Output (image) window shape; each axis must be <= n_fft.	required
n_fft	Any	FFT size per axis, >= the input size. One output pixel is 1 / n_fft cycles per input pixel.	required
offset	Any	Output-window centre offset in output pixels, (dy, dx).	0.0

MFTPropagator

MFTPropagator(in_shape, out_shape, samples, *, offset=0.0, stack=False, backend=None)

Bases: Propagator

Matrix Fourier transform with arbitrary sampling (Soummer et al. 2007).

E = A_y u A_x^T with A[j, i] = exp(-2 pi i (i - c_in)(j - c_out) / N) / sqrt(N), where N (samples) is the equivalent FFT size and need not be an integer. Several N values (one per wavelength) give a stacked propagator mapping (..., L, ny, nx) to (..., L, my, mx).

Parameters:

Name	Type	Description	Default
in_shape	Any	Input and output grid shapes.	required
out_shape	Any	Input and output grid shapes.	required
samples	Any	Equivalent FFT size N = 1 / (input pitch x output pitch) in dimensionless units: a scalar or (y, x) pair, or with stack=True an (L,) or (L, 2) array, one entry per wavelength.	required
offset	Any	Output-window centre offset in output pixels, (dy, dx).	0.0

FocalPlanePropagator

FocalPlanePropagator(pupil_shape, pupil_pitch, wavelengths, pixel_scale, out_shape, *, offset=0.0, method='auto', backend=None)

Bases: Propagator

Fraunhofer propagation from a pupil to a detector, at one or more wavelengths.

Maps a pupil field (..., L, ny, nx) (one slice per wavelength) to an image field (..., L, my, mx). With a single wavelength the wavelength axis is still present (L = 1).

The engine is chosen per call to :meth:build: a zero-padded FFT when the equivalent FFT size N = lambda / (pupil pitch x pixel scale) is an integer at every wavelength and the FFT is cheaper, otherwise a stacked matrix Fourier transform, which handles any sampling (including undersampled detectors) and broadband light exactly.

Parameters:

Name	Type	Description	Default
pupil_shape	Any	(ny, nx) of the pupil grid.	required
pupil_pitch	float	Pupil pixel pitch in metres.	required
wavelengths	Any	Wavelength(s) in metres.	required
pixel_scale	float	Detector (or model) pixel scale in radians.	required
out_shape	Any	(my, mx) of the image.	required
offset	Any	Image-window centre offset in pixels (dy, dx).	0.0
method	str	"auto", "fft" or "mft".	'auto'

n_wavelengths property

n_wavelengths

Number of wavelength slices L.

AngularSpectrumPropagator

AngularSpectrumPropagator(shape, pitch, wavelength, distances, *, paraxial=False, backend=None)

Bases: Propagator

Scalar angular-spectrum (Fresnel/near-field) propagation over a distance.

E_z = IFFT[ FFT[u] * H ] with H = exp(i 2 pi z sqrt(1/lambda^2 - f^2)) and evanescent components removed. Shapes are unchanged. Several distances give a stack (..., Z, ny, nx) from a field broadcast against it.

Parameters:

Name	Type	Description	Default
shape	Any	(ny, nx) of the field.	required
pitch	float	Sample pitch in metres.	required
wavelength	float	Wavelength in metres.	required
distances	float | Sequence[float]	Propagation distance(s) in metres (positive downstream).	required
paraxial	bool	Use the Fresnel (paraxial) transfer function exp(-i pi lambda z f^2) (global phase dropped).	False

centered_coordinates

centered_coordinates(n, pitch=1.0, offset=0.0, *, dtype=None, backend=None)

Pixel-centre coordinates (i - (n - 1) / 2 + offset) * pitch (CONVENTIONS 1.2).

The values are computed in float64 and then cast, so a float32 grid equals the float32 cast of the float64 one on either device.

Parameters:

Name	Type	Description	Default
n	int	Number of pixels along the axis.	required
pitch	float	Pixel pitch (any unit; metres for pupil planes).	1.0
offset	float	Shift of the grid in pixels.	0.0
dtype	Any	Floating-point dtype of the result. Defaults to the backend's working precision (:attr:Backend.real_dtype), which is float64 for the default CPU backend.	None
backend	BackendLike	Where to build the array: a :class:~aocore.Backend or a device name (see :func:~aocore.get_backend). Default: host NumPy, double precision.	None

aocore.metrics

Wavefront error metrics that respect phase-retrieval ambiguities.

Phase retrieval cannot see some components of the wavefront: piston always, absolute tip/tilt when the image position is free, and for a single in-focus image of a centro-symmetric pupil the twin phi(x) -> -phi(-x). These helpers compare an estimate with a reference after removing exactly those components, weighted by pupil intensity.

remove_modes

remove_modes(opd, pupil, remove='piston')

Subtract the intensity-weighted least-squares fit of the named modes.

remove is "piston", "tip", "tilt", "tiptilt" (piston + tip + tilt), "segment_piston", a list of these, or None.

rms

rms(opd, pupil, remove='piston')

Intensity-weighted RMS over the pupil after :func:remove_modes (CONVENTIONS 4.1).

pupil is a :class:~aocore.Pupil or a pupil amplitude array.

rms_unweighted

rms_unweighted(opd, mask=None)

Plain quadratic mean of opd over mask, piston included (CONVENTIONS 4.1 variant).

Returns sqrt(mean(opd[mask] ** 2)): every pixel inside the mask counts equally, whatever the pupil transmission, and nothing is removed, so a constant OPD c gives |c|. This is the figure that packages such as makewfs report for a residual OPD that is already piston-free. Use :func:rms for the convention's intensity-weighted, piston-removed RMS.

Works on NumPy and CuPy arrays without copying the OPD to the host; the sum is accumulated in float64 on the array's device, and only the scalar result crosses to the host.

Parameters:

Name	Type	Description	Default
opd	Any	OPD map (any shape; the result is in its units). NumPy or CuPy.	required
mask	Pupil | Any	Pixels to include: a boolean array, an amplitude-like array whose non-zero pixels count, or a :class:~aocore.Pupil (its :attr:~aocore.Pupil.mask). Same shape as opd. None uses the whole array. Values outside the mask are ignored, even if they are NaN.	None

rms_tiptilt_removed

rms_tiptilt_removed(opd, pupil)

Intensity-weighted RMS with piston, tip and tilt removed (CONVENTIONS 4.1 variant).

The residual after subtracting the weighted least-squares fit of c0 + c1 x + c2 y, with weights |pupil|^2: sqrt(sum a^2 r^2 / sum a^2). Equal to rms(opd, pupil, "tiptilt"), but it runs on the array's device: NumPy and CuPy inputs stay where they are, the reductions run in float64 there, and only a handful of scalars (the 2x2 normal equations and the result) cross to the host.

Parameters:

Name	Type	Description	Default
opd	Any	2-D OPD map (the result is in its units). NumPy or CuPy.	required
pupil	Pupil | Any	A :class:~aocore.Pupil or a pupil amplitude array of the same shape. Pixels with zero amplitude do not count (their OPD may be NaN).	required

strehl_from_rms

strehl_from_rms(rms_opd, wavelength)

Marechal approximation exp(-(2 pi rms / lambda)^2) (alias of strehl_marechal).

wavefront_error

wavefront_error(estimate, truth, pupil, *, remove='piston', allow_twin=False)

RMS difference between two OPD maps over the pupil, modulo ambiguities.

Parameters:

Name	Type	Description	Default
estimate	Any	OPD maps (same units; the result is in those units).	required
truth	Any	OPD maps (same units; the result is in those units).	required
remove	RemoveSpec	Modes removed from the difference before the RMS (see :func:remove_modes).	'piston'
allow_twin	bool	Also compare against the twin -estimate(-x) and return the smaller error (for single-image retrieval without diversity). Only meaningful for pupils symmetric under 180-degree rotation.	False

aocore.unwrap

Least-squares phase unwrapping over arbitrary (obscured, segmented) pupils.

Iterative-transform algorithms return phase wrapped into (-pi, pi]. The weighted least-squares unwrapper of Ghiglia & Romero (JOSA A 11, 107, 1994) finds the smooth phase whose finite differences best match the wrapped differences of the input, using only pixel pairs that are both inside the mask. It is solved by conjugate gradients preconditioned with the unweighted Neumann Poisson solver (a DCT), on CPU or GPU.

Least-squares unwrapping is exact when the true phase changes by less than pi between neighbouring pixels. Spiders and segment gaps split a pupil into disconnected regions whose relative 2 pi multiples the data cannot fix. With bridge=True (default) each region's multiple of 2 pi is chosen so that the regions join a smooth low-order polynomial surface; that assumes the true wavefront is continuous across the gap, which is right for spider vanes but is a guess for independently pistoned mirror segments.

wrap

wrap(phase)

Wrap phase into (-pi, pi].

unwrap_phase

unwrap_phase(wrapped, mask=None, *, weights=None, tol=1e-08, max_iter=500, bridge=True, bridge_order=4, device=None)

Weighted least-squares phase unwrapping.

Parameters:

Name	Type	Description	Default
wrapped	Any	(ny, nx) wrapped phase in radians (NumPy or CuPy).	required
mask	Any	Boolean support; pixels outside are ignored and returned as 0. Defaults to everything.	None
weights	Any	Optional per-pixel reliability in [0, 1] (e.g. pupil amplitude); a pair is weighted by the smaller of its two pixel weights.	None
tol	float	Relative residual for the conjugate-gradient stop.	1e-08
max_iter	int	Conjugate-gradient iteration limit.	500
bridge	bool	Choose the relative 2 pi multiples of disconnected regions so they join a smooth surface (see module notes).	True
bridge_order	int	Total degree of the polynomial surface used for bridging.	4
device	BackendLike	Backend override; defaults to the backend of wrapped.	None

Returns:

Type	Description
Unwrapped phase on the same backend, zero outside ``mask``. Within each
connected region its mean matches the wrapped input's circular mean.

aocore.backend

Array backends: NumPy/SciPy on the CPU and CuPy on CUDA GPUs.

Numerical routines in the AO stack take their array namespace from a :class:Backend, so the same algorithm runs on either device. A backend fixes three things:

• device — "cpu" (NumPy arrays, :mod:scipy.fft with all cores) or "gpu" (CuPy arrays, cuFFT/cuBLAS);
• precision — "double" (float64/complex128) or "single" (float32/complex64). CPU defaults to double, GPU to single;
• the unitary FFT used by every propagator.

Arrays returned by a GPU solver stay on the GPU. :func:to_numpy is the one explicit host boundary.

Backend dataclass

Backend(device, precision)

One array namespace, device and floating-point precision.

Use :func:get_backend rather than constructing this directly.

Attributes:

Name	Type	Description
device	str	"cpu" or "gpu".
precision	str	"double" or "single".

xp property

xp

The array namespace: :mod:numpy or :mod:cupy.

is_gpu property

is_gpu

Whether arrays live on a CUDA device.

real_dtype property

real_dtype

float64 for double precision, float32 for single.

complex_dtype property

complex_dtype

complex128 for double precision, complex64 for single.

eps property

eps

Machine epsilon of :attr:real_dtype.

asarray

asarray(value, dtype=None)

Move value onto this backend.

dtype may be a NumPy dtype, "real" or "complex" (this backend's working precision), or None to promote real input to :attr:real_dtype and complex input to :attr:complex_dtype. Integer and boolean arrays keep their type when dtype is None.

zeros

zeros(shape, dtype='real')

Zero-filled array in working precision ("real"/"complex").

empty

empty(shape, dtype='real')

Uninitialized array in working precision ("real"/"complex").

fft2

fft2(array, *, axes=(-2, -1))

Unitary (norm="ortho") two-dimensional FFT, uncentred.

ifft2

ifft2(array, *, axes=(-2, -1))

Unitary inverse of :meth:fft2 (also its adjoint).

rfft2

rfft2(array, *, axes=(-2, -1))

Unnormalized real-input 2-D FFT (used for convolutions).

irfft2

irfft2(array, shape, *, axes=(-2, -1))

Inverse of :meth:rfft2 for an output of shape.

blas_limit

blas_limit(work)

Context limiting CPU BLAS threads for products of work multiply-adds.

A no-op on the GPU. Entering and leaving costs a few microseconds.

dot

dot(a, b)

Re(sum(conj(a) * b)) over all elements, as a Python float.

On the CPU this avoids BLAS level-1 calls: OpenBLAS threads ddot above ~10k elements, and those threads stall (~1 ms per call) whenever the cores are busy. Use it for every scalar product in iterative loops.

to_numpy

to_numpy(array)

Copy an array to host NumPy (no-op on the CPU).

scalar

scalar(value)

Convert a 0-d array to a Python float (a device sync on GPU).

synchronize

synchronize()

Block until queued device work finishes (no-op on the CPU).

random

random(seed)

A host :class:numpy.random.Generator.

Random draws (initial guesses, noise) are made on the host and moved to the device so a seed gives the same numbers on CPU and GPU.

gpu_available cached

gpu_available()

Return whether CuPy is importable and sees at least one CUDA device.

get_backend

get_backend(device='cpu', precision=None)

Return the :class:Backend for a device and precision.

Parameters:

Name	Type	Description	Default
device	BackendLike	"cpu", "gpu" (CuPy; raises if unavailable), "auto" (GPU when available, else CPU), or an existing :class:Backend, which is returned unchanged unless precision overrides it.	'cpu'
precision	str | None	"double" or "single". Defaults to double on CPU, single on GPU.	None

to_numpy

to_numpy(array)

Return array as a host NumPy array (copying from the GPU if needed).

Non-array values (lists, scalars) go through :func:numpy.asarray.

backend_of

backend_of(array, precision=None)

Infer a backend from an array (CuPy arrays give the GPU backend).

Without precision, single precision is chosen for float32/complex64 input and double otherwise.

Simulation

solvephase.simulate

Synthetic data for tests, tutorials and benchmarks.

For physically detailed data, see :mod:solvephase.interop: atmospheric OPD from pyturb and detector frames from getframes.

random_aberration

random_aberration(pupil, rms, *, n_modes=36, power=1.0, start=2, seed=None)

A random Zernike aberration with a power-law spectrum.

Coefficient j (counting from start in Noll order) has standard deviation proportional to j^-power; the map is scaled to exactly rms metres RMS over the pupil (piston removed).

Parameters:

Name	Type	Description	Default
pupil	Pupil	The pupil.	required
rms	float	Target RMS OPD in metres.	required
n_modes	int	Number of Zernike modes.	36
power	float	Spectral slope (0 = white, ~1 typical of polished optics).	1.0
start	int	First Noll index (2 includes tip/tilt, 4 starts at defocus).	2
seed	Any	Seed or :class:numpy.random.Generator.	None

Returns:

Type	Description
``(ny, nx)`` OPD map in metres.

simulate_images

simulate_images(model, opd=None, *, photons=1000000.0, background=0.0, read_noise=0.0, noise=True, amplitude=None, seed=None)

Simulate photon-counting focal-plane images for a model.

Parameters:

Name	Type	Description	Default
model	FocalPlaneModel	The forward model.	required
opd	Any	Pupil OPD in metres (default flat).	None
photons	float | Any	Photons collected per channel from the whole PSF (scalar or (K,)).	1000000.0
background	float | Any	Background photons per pixel per channel.	0.0
read_noise	float	Gaussian read noise standard deviation in electrons.	0.0
noise	bool	Add Poisson and read noise; False returns the expectation.	True
amplitude	Any	Optional pupil amplitude override.	None
seed	Any	Seed or :class:numpy.random.Generator.	None

Returns:

Type	Description
``(K, my, mx)`` host array in photo-electrons (quantum efficiency 1).