prysm v0.22#

Released (eventually)

This release, the first in several years, marks the beginning of a split between the core library which will have better stability between releases, and a new prysm.x-perimental (shorthand, x or x/) directory for new modules which have not had the ~10 years of development that the core library has.

A basic summary of the changelog is:

  • Expanded propagation module, particularly for multi-plane problems such as coronagraphs.

  • New experimental modules for optimization, polarization, and massive expansion of the raytracing module.

  • Vast expansion of the polynomials module, including analytical derivatives and sums+derivatives using Clenshaw’s method (from Forbes’ Q polynomial papers) for extremely efficient computation of surfaces and their partial derivatives, essential for high-performance raytracing.

  • Adjoint / gradient backpropagation routines for a large subset of the API, allowing gradient-based optimization of a wide range of problems in optics such as:

      • Phase Retrieval

      • Coronagraph mask optimization

      • In raytracing: Using backpropagation to obtain partial derivatives of a

      cost function with respect to tolerances and other parameters, similar to wavefront differential tolerancing (which is also implemented).

  • Several new features and bug fixes in several core modules. See the release notes for each.

The new x/ modules are detailed in the release notes below.

Core library#

Polynomials#

A breaking change has been made by renaming xxx_sequence to _seq, to be consistent with using _der for derivatives.

Utilities to orthogonalize and normalize modes over arbitrary apertures with special routines for annular apertures:

  • distort_annular_grid() - deforms an ordinary (r, t) coordinate mesh to represent a new unit domain for using conventional orthogonal polynomial sets. For example, vanilla zernikes can be generated on the return from this function, and they will be orthogonal over the annulus.

  • orthogonalize_modes() - using a QR factorization (more stable gramm-schmidt style process) approach to generate numerically orthogonalized modes.

  • normalize_modes() - normalize modes so that they are orthonormal or otherwise (e.g., to=’std’).

See Ins and Outs of Polynomials for usage examples

Rich XY polynomial capability has been introduced:

Partial derivatives of the XY monomials have been added. Three named callables plus their sequence variants are provided rather than one with an axis= kwarg, so the chain rule reads naturally at call sites:

Additionally, Laguerre polynomials have been introduced, which can be used for generating Laguerre-Gaussian beams:

Derivatives of the Dickson polynomials of the first and second kind are now available via the standard differentiated three-term recurrence:

New functions for cartesian derivatives of Zernikes, as well as using Clenshaw’s method to compute their sum and x- and y- derivatives in one sweep:

Derivatives of Q polynomials, as for Zernikes:

All of the _seq polynomial functions have been revised. Previously, they returned generators to allow weighted sums of extremely high order expansions to be computed in a reduced memory footprint. This lead to the most common usage being basis = array(list(xxx_sequence())). This benefit has been more theoretical than practical, e.g. a 1000-term expansion at 4096x4096 requires much less than 1GB of memory to hold the dense coefficient array. Now equivalent usage is basis = xxx_seq(), which returns the dense array of shape (K,N,M) directly (K=num modes, (N,M) = spatial dimensionality).

Propagation#

A significant breaking change has been made to how fixed sampling / DFT based propagations are done. The previous API used a global matrix DFT or CZT object and the user passed arguments that implicitly define a grid to propagation functions. This created significant complexity with backpropagation. The API now requires two lines of code to do the equivalent, but usage is significantly simplified for models that repeatedly propagate between the same samplings. The previous usage looked like:

WAS.py#
1from prysm.propagation import wavefront as WF
2
3wf = WF(...)
4efoc = wf.focus_fixed_sampling(efl, focal_dx, Nfocal)
5epup = efoc.unfocus_fixed_sampling(efl, wf.dx, wf.data.shape)

The new usage looks like:

IS.py#
1from prysm.propagation import Wavefront as WF
2
3wf = WF(...)
4mdft = wf.prepare_executor(efl, focal_dx, Nfocal) # can add shifts
5efoc = wf.focus_dft(mdft)
6epup = efoc.unfocus_dft(mdft)

This change significantly simplifies logic around backpropagation.

Functions have been renamed:

New features:

prysm.propagation has been split into a few files. This should cause no change for users.

The versions of propagation functions that are not methods of Wavefront previously accepted Wavefront or Richdata instances, with optional arguments such as dx. They now operate only on arrays; users should use one API or the other, not both at once.

fttools#

The fttools module has been reworked, and greatly simplified. The module was created to serve the needs of the propagation module, which resulted in an indirect form of input to the DFT, where instead of passing the spatial and frequency grids, abstractions over them were passed. In the new module, the matrix DFT and chirp Z transform are implemented directly. Most users do not interact with the fttools module directly. The previous API signature was

from prysm.fttools import mdft
mdft.dft2(self, ary, Q, samples_out, shift=(0, 0))

The new API is:

from prysm.fttools import MDFT
mdft = MDFT(x, y, fx, fy)
mdft() # call it
mdft.adjoint() # adjoint operator for backpropagation

The CZT class works in the same way.

Segmented Systems#

Bayer#

  • prysm.bayer.wb_scale has been renamed to wb_postscale()

  • wb_postscale() now has additional safe and saturation kwargs for colorimetrically correct handling of saturation

  • wb_prescale() now has additional safe and saturation kwargs for colorimetrically correct handling of saturation

    New:

  • demosaic_deinterlace() for deinterlace-style demosaicing, which cuts resolution by a factor of two but imparts no blur or color channel crosstalk

  • assemble_superresolved() to assemble a trichromatic image from super-resolved Bayer color planes

PSF#

New:

  • airydisk_efield(); new function that returns the complex electric field of an Airy disk.

i/o#

New:

Detectors#

new apply_lut() function, and associated kwarg lut when initializing a Detector instance. This adds the capability to simulate detector nonlinearity that is homogeneous over the whole array via arbitrary mapping

mathops#

The global prysm.conf.config.precision setting now accepts any real floating dtype-like value, not only the previous 32 and 64 integer selectors.

from prysm.conf import config
import numpy as np

config.precision = np.float16
config.precision = "float32"
config.precision = np.dtype("float64")

Convenience methods to swap the backend to Cupy, pytorch, apple MLX, or MKL-FFT:

There is no special support for either library beyond these functions.

The math operations shim now exposes prysm.mathops.linalg.

  • array_to_true_numpy() converts arrays from the above libraries into numpy arrays for plotting, serialization, etc.

  • the plot2d() method of RichData now has an extend keyword argument, which controls the extension of the colorbar beyond the color limits

Interferograms#

  • pad() now accepts either samples to pad by a number of pixels, or shape to pad to an explicit output shape. The fill value is now passed separately as value

  • slope() returns x, y, and non-directional slope as RichData instances

  • save_zygo_dat() saves interferograms as Zygo .dat files

Geometry#

New shapes available:

Coordinates#

New function chebygauss_quadrature_xy() for generating optimally spaced spiral sampling patterns

A suite of functions has been added for working with homographies, a type of 3D transformation that can used for projecting surfaces onto inclined planes and other geometric manipulations:

Thin Films#

The API has been changed and the user passes separate thickness and index arrays. This makes the module easier to use and less error prone by swapping the position of the index and thickness variables. The substrate also is now handled separately, and no longer needs to be passed as a very thick film layer. This second change avoids subtle etalon effect errors that could occur.

WAS.py#
1stack = [
2    (n_MgF2, .150),
3    (n_C7980, 10_000),
4]
5r, _ = thinfilm.multilayer_stack_rt(stack, wvl, 'p')
IS.py#
1indices = [n_MgF2]
2thicknesses = [.150]
3# note last arg is substrate index
4r, _ = thinfilm.multilayer_stack_rt(indices, thicknesses, wvl, 'p', n_C7980)

snell_aor() previously took a degrees kwarg, while the sibling functions brewsters_angle() and critical_angle() took deg. All three now use deg consistently.

Thin Lens#

Several new functions have been added:

Refractive#

  • internal_transmission() computes internal transmission from slab thickness, extinction coefficient, and wavelength

eXperimental Modules#

A total of seven new x/ modules have been introduced in this release. The most substantial is the new optimization module, as well as upgrades to raytracing. Additional modules include polarization, fibers, several types of interferometers and the shack-hartmann wavefront sensor, as well as deformable mirrors.

x/optym#

This module is for (gradient-based) optimization with numerous optimizers. The API is different to scipy.optimize. In scipy, you call minimize in a fire-and-forget way. The only way you can interact during optimization is via a callback function. In optym, optimizers only know how to step() and advance one iteration. At the lowest level of usage, you could just call` step in a loop and access the optimizer and its internal state, and thing you are optimizing, on every iteration. You could for example update a plot, or decide to stop if something in your model got astray.

The package also has a type called Governor, which inspects the optimizer and its internal state on each iteration, and decides when to stop. For example, one governor checks that the function has stopped decreasing to within a tolerance. Another checks that the gradient has become ~0.

All of the optimizers are GPU-compatible and 32-bit compatible. Many are able to run at 1000’s of iterations per second on relatively inexpensive problems.

The basic building blocks are the user providing a minimal function cost(x) -> (f, g) where f is the scalar cost function, and g is d/dx of f. This can be provided as a more formal Problem object which can be more computationally efficient with certain optimizers. This is driven to the minimum value of f using an optimizer. Optimizers have a simple API: the user calls .step() -> (x, f, g) to iterate one step, returning the value of the parameter vector, cost, and gradient before the step. The new x vector is then accessible at optimizer.x. In most cases, you will not step the optimizer directly; instead calling run_until() to run until a stopping condition is met, or run_for() to run for a fixed number of steps. run_until can be given any Governor to define one or more stopping condition. When multiple governors are given, the first of them to call stop will end iteration.

In addition to API level documentation that describes each of these items in detail, there are several new optimization-related tutorials and how-tos:

A strong Wolfe conditions line searcher is available via ls_strong_wolfe().

Optimizers#

LBFGSB is a wrapper around the low level scipy L-BFGS-B driver, in prysm’s interface. It wraps the Fortran optimizer for older versions of scipy, or the new C driver for newer versions of scipy. PrysmLBFGSB is a bespoke, pure python/numpy/numpy-like implementation which is competitive in speed on large problems, and slower on very small problems due to python overheads. It supports reduced precision such as 32-bit, and can run on the GPU if the backend is swapped to CuPy. The LBFGSB driver is CPU and 64-bit only. Damped Least Squares supports all lagrangian constraints, equality, greater than, less than. All other optimizers support box constraints, l < x < u.

The point of equality is approximately a thousand variables. For smaller problems, PrysmLBFGSB is a factor of a few (2-5x) slower than the Scipy driver due to python overhead. This may be reduced when using the python JIT PYTHON_JIT=1 python3.15 my_script.py.

The PrysmLBFGSB uses some different methods of solving for the cauchy point and then performing subspace minimization. It appears to be able to converge more reliably and in many cases in fewer iterations than the scipy L-BFGS-B driver.

Governors#

These are what trigger the end of iteration with run_until().

Two are special or meta, and combine the others:

  • AnyGovernor; # stop when any of the other governors are met

  • AllGovernor; # stop when all of the other governors are met

Ordinary governors:

The convergence may be plotted using plot_convergence().

Benchmark and sample problems#

Problem forms of these come pre-populated with f, g, fg, h, and hvp methods. You can use them with optimizers, or use them as examples of building a Problem.

Activation functions, discretizers, cost/loss functions, and operators are differentiable building blocks for building optimization problems:

Activation functions and discretizers#

Cost or loss functions#

operators#

x/polarization#

New module for Jones calculus and other polarization calculations. Included is an adapter that generalizes all routines within the propagation module to propagation of Jones states. A tutorial is available at Intro to Polarization and a How-To at Polarized Propagation.

Jones Vectors#

Jones Matrices#

Conversion to Mueller matrices and simple data reduction with Pauli spin matrices:

x/fibers#

New module with routines to parametrically study cylindrical step index fibers and wavesguides. Contains functions to identify the \(LP_{\ell{}m}\) modes of single and multi-mode fibers as well as evaluate them numerically. Also contains the mode overlap integral used to model coupling of complex E-fields into fibers and waveguides.

The main user-facing routines are:

The tutorial on Single-Mode Fibers covers basic usage. Multimode fibers have the same as single mode except compute_LP_modes() is used instead of smf_mode_field().

x/psi, x/pdi, x/sri, x/shack_hartmann#

These four modules are for the modeling of Shack-Hartmann sensors and two types of pinhole and/or fiber/waveguide based interferometers. Extensive phase shifting interferometry (PSI) reconstruction capability is included, both of wavefront phase as well as complex E-field.

  • Forward modeling of Shack Hartmann wavefront sensors using shack_hartmann() and the propagation module

  • Forward modeling of Phase Shifting Point Diffraction Interferometers, aka Medecki interferometers using PSPDI and the routines and consants of x/psi

  • Forward modeling of Self-Referenced Interferometers (SRIs), which use a fiber to generate the reference wave using light from the input port using SelfReferencedInterferometer

  • PSI routines:

All of the PSI routines revolve around the Scheme type, which contains the phase shifts and associated sine and cosine weights

    • psi_accumulate() for accumulating the sums of de groot’s formalism, an essential intermediate step in full complex E-field reconstruction and differential reconstruction

x/dm#

  • copy() method to clone a DM, when e.g. the two DMs in a system are the same

  • new Nout parameter that controls the amount of padding or cropping of the natural model resolution is done. The behavior here is similar to PROPER

  • the forward model of the DM is now differentiable. render_adjoint() applies the adjoint of render()

  • rotation definitions have been changed, and a related bug that would cause a transposition of the DM surface for some rotations fixed.

x/raytracing#

The sequential raytracing module has been substantially reworked. The surface API is now built around explicit sag-shape objects and a simpler Surface container instead of the previous factory/metaclass-heavy implementation. Built-in shapes include planes, spheres, conics, off-axis conics, even aspheres, Forbes Q2D, Zernike, XY, Chebyshev, Jacobi, toroid, and biconic sags, with consistent F / FFp height-and-slope contracts. The Spencer-and-Murty tracer now carries richer status information through apertures, interaction failures, gratings, and inactive rays.

A new Prescription wrapper provides an ergonomic sequential-system interface over the low-level surface list. It stores fields, wavelengths, EPD, ambient index, and reference wavelength; can build simple refractive sequences; can solve and set the paraxial image plane or target focal length; and exposes one-call launch, trace, RMS spot, focus optimization, transverse ray fan, wave aberration fan, and layout/fan plotting helpers.

The launch and first-order layers have grown into reusable user APIs: Field describes finite or angular object fields, Sampling constructs chief-ray, fan, cross, rectangular, hexagonal, and radial-spiral pupil samplings, and launch creates the corresponding ray bundles. The new paraxial solver reports system matrix quantities, image distance, effective focal length, front/back focal length, and a FirstOrderProperties bundle.

Analysis, design, and tolerancing tools have been added around the tracer. Analysis routines now cover transverse ray aberration, wavefront OPD and Zernike fitting, distortion, field curvature, axial color, and lateral color. The design layer introduces variable factories for curvature, radius, conic constant, indexed coefficients, position, and air/glass thickness; operands for RMS spot radius, ray height, boresight, EFL, BFL, paraxial image distance, wavefront RMS, Zernike coefficients, distortion, and field curvature; and a problem wrapper compatible with x/optym least-squares solvers. The tolerance layer adds perturbation distributions, finite-difference sensitivity tables, and Monte Carlo tolerance runs.

Zemax .zmx and Code V .seq readers can now construct prysm surface prescriptions, including support code for glass-name lookup through the optional refractiveindex.info SQLite database. The materials module also provides air, vacuum, glass, lookup, and load_material_db helpers, and parsing recognizes mirror/ambient special cases.

Plotting has been expanded from simple ray-path views to prescription-aware optical layout plots, transverse ray aberration plots, wave-aberration fan plots, spot diagrams, and lens-outline construction for refractive groups. Geometric helpers for ray-plane, ray-sphere, ray-conic, line-intersection, and closest-approach calculations have been factored into small modules so surface intersection and analysis code share the same math.

Performance Optimizations#

  • angular_spectrum_transfer_function() has been optimized. The new runtime is approximately the square root of that of the old. For example, on a 1024x1024 array, in version 0.21 this function took 31 ms on a desktop. It now takes 4 ms for the same array size and output

  • prysm.propagation.Wavefront.intensity() was made ~60% faster.

  • Matrix DFT and CTZ transforms for asymmetric axis sizes and rectangular arrays have been optimized by a factor of 2-8x depending on how rectangular and how asymmetric the array size is.

  • rectangle() has been optimized when the rotation angle is zero

  • rectangle() has been optimized when the coordinates are exactly square/cartesian (not rotated)

  • read_zygo_dat() now only performs big/little endian conversions on phase arrays when necessary (little endian systems), which creates a slight performance enhancement for big endian systems, such as apple silicon

Bug Fixes#

  • The sign of thin_lens() was incorrect, requiring a propagation by the negative of the focal length to go to the focus. The sign has been swapped; (wf * thin_lens(f,...)).free_space(f) now goes to the focus

  • mtf_from_psf() as well as the ptf and otf functions used the wrong pixel as the origin for normalization, when array sizes were odd. This has been fixed

  • a bug in scipy.special.factorial2 has been fixed in a recent version. Like all respectable software, prysm depended on that bug. Q2D polynomials would return NaN for m=1, n=0 (Q-coma) with scipy’s bugfix. This has been corrected within prysm in this version, and Q-coma is no longer destined for NaN

  • prysm.polynomials.zernike.barplot() and barplot_magnitudes() now apply axis labels to the correct axis when plotting on a figure with multiple axes

  • fixed a bug in prysm.psf.encircled_energy() where x,y axes were double meshgrided

  • Non-square arrays now work properly in read_codev_gridint() and write_codev_gridint()

  • strip_latcal() no longer fails to work properly

  • demosaic_malvar() had a scaling error for the blue channel corrected.

  • prysm.polynomials.cheby functions now are scaled properly when called on N dimensional coordinate arrays

  • A latent bug in xy_seq() was also fixed: the internal monomial table was generated from dickson1_seq() with alpha=0, but the Dickson convention gives \(D_0 = 2\) rather than the monomial \(x^0 = 1\). This meant xy_seq([(0,0)], x, y) returned an array of 4 while xy(0, 0, x, y) returned 1. Both forms now return 1 consistently.

  • Some bugs in edge cases with clenshaw’s method for sums of Q polynomials were fixed. Notably, indexing beyond the end of the array, requiring the user to pass an extra zero coefficient at the end of the vector.

  • Chirp Z Transforms for nonsquare or shifted arrays have been fixed. Previously the chirp was not properly aligned to the shifted coordinates, resulting in incorrect results.

Code Removal and other Breaking Changes#

Packaging has moved from setuptools configuration in setup.py, setup.cfg, and MANIFEST.in to pyproject.toml using hatchling. The minimum supported Python version is now 3.10. The top-level prysm.__version__ attribute has also been removed; use importlib.metadata.version('prysm') instead.

Project linting has moved to Ruff.

Numerous features related to MTF benches have been removed. The code was extremely old, had incomplete test coverage, and is rarely used:

  • prysm.io.read_trioptics_mtfvfvf()

  • prysm.io.read_trioptics_mtf_vs_field()

  • prysm.io.read_trioptics_mtf()

  • the entire mtf_utils module

  • sample Trioptics mht and txt files

Within the geometry module, all functions now use homogeneous names of x, y, r, and t for arguments. The circle() and truecircle() routines have had some of their arguments renamed.

The following functions have been removed from the polynomials submodule:

  • separable_2d_sequence

  • mode_1d_to_2d

  • sum_of_xy_modes

They assumed strict separability of the two axes, with no cross terms. This can be acheived by having terms where only m or n is positive in the new XY routines. In general, suppressing cross terms artificially is not intended and the functions have been removed to avoid confusion.

The degredations module has been modernized, and two bugs have been fixed in doing so. The magnitude of jitter now matches more common modern formalisms, and is twice as large for the same “scale” parameter previously. The smear parametrization has been modified from (mag,ang) to (mag x, mag y). Pass width=0 or height=0 for monodirectional smear. This also corrects a bug, in which only the diagonal elements of the transfer function were corectly populated with sinc() when rotation != 0 previously.

Encircled Energy related functions have been moved from the psf module to the otf module:

prysm.io.read_zygo_dat() was reworked to improve code reuse with the new write function. In doing so, some of the nesting in the dictionary representation of the metadata has become flat or unnested. The reading of phase and intensity is unchanged.