API Reference

ScatteringTransforms.jl — Native Julia implementation of wavelet scattering transforms

Surfaces: 1D signals, 2D and 3D gridded fields, scattered planar points (NUFFT), and the sphere — both on a structured grid and at scattered points. Each has a monogenic variant; the gridded ones also have a localized (Mallat) field output and a multi-resolution second order. The spherical and scattered paths run dependency-free by default, with fast paths supplied by extensions.

Quick Start

Names live in submodules and are reached by their full path; alias the package for brevity.

using ScatteringTransforms: ScatteringTransforms as ST

# 1D scattering (N = signal length, J = number of octaves)
signal = randn(1024)
st = ST.Scattering1D.ScatteringTransform1D(1024, 8; Q=1, max_order=2)
coeffs = st(signal)

# 2D planar scattering (N = image size, J = scales, L = orientations)
image = randn(256, 256)
st2d = ST.Scattering2D.ScatteringTransform2D((256, 256), 4; L=8, max_order=2)
coeffs2d = st2d(image)

Implementation Notes

• FFT-based convolutions for O(N log N) performance
• Frequency-domain Morlet filter banks
• Breadth-first CSR path list, walked grouped by first-order wavelet
• Modular extensions for spherical (NUFSHT) and GPU support

References

• Mallat (2012): Group invariant scattering. Comm. Pure Appl. Math.
• Bruna & Mallat (2013): Invariant Scattering Convolution Networks. IEEE PAMI.
• Cheng & Ménard (2021): How to quantify fields or textures? A guide to the scattering transform.

ScatteringTransform1D{T,V,M,P,Tree,FB,G}

1D scattering transform: a filter bank, the admissible path tree, a spectral plan, and the workspace the cascade runs in. Every array field is a type parameter, so the same struct holds CPU, GPU or static storage.

Fields

• filter_bank: pre-computed 1D filter bank
• tree: admissible scattering paths
• groups: (j1, children) from tree, longest-first — the order the cascade walks
• max_order: maximum scattering order (1 or 2)
• plan: spectral transform plan (in-core direct sum by default; FFTW fast path if loaded)
• buffer_input: complex buffer for real→complex promotion, and multiply scratch
• buffer_signal_fft: the input spectrum, read-only for the whole cascade
• buffer_conv: inverse-transform output
• dims: the signal's length, as a 1-tuple
• buffer_u1_fft: the spectrum of the current first-order modulus, reused across that wavelet's children. The modulus itself needs no buffer of its own — the cascade writes it straight into buffer_input, which is what its transform reads.

ScatteringTransform2D{T,M,R}

2D planar scattering transform with oriented wavelets and pre-allocated workspace.

Type Parameters

• T: Real element type (Float32, Float64, ...)
• M: Complex matrix type for buffers (Matrix{Complex{T}}, CuMatrix{Complex{T}}, ...)
• R: Real matrix type for modulus buffers (Matrix{T}, ...)

Fields

• filter_bank: Pre-computed 2D filter bank
• max_order::Int: Maximum scattering order (1 or 2)
• plan: spectral transform plan (in-core direct sum by default; FFTW fast path if loaded)
• buffer_input: Complex matrix for real→complex promotion (zero alloc)
• buffer_signal_fft: Preserved copy of signal FFT (buffer_conv gets overwritten)
• buffer_conv: Complex matrix for IFFT output
• dims: the field's spatial size
• buffer_u1_fft: the spectrum of the current first-order modulus, reused across that wavelet's children. The modulus itself needs no buffer of its own — the cascade writes it straight into buffer_input, which is what its transform reads.

ScatteringTransform3D([T=Float64,] N, J; n_orient=6, max_order=2, spectral=AutoSpectralBackend())

Build a 3D volumetric scattering transform for N = (Nz, Ny, Nx) volumes over J scales and n_orient sphere directions. The element type is positional, as for zeros(T, …); omit it for Float64.

(st::ScatteringTransform3D)(volume) -> ScatteringCoefficients2D

Apply the 3D scattering transform. (Coefficients use the scales×orientations container.)

scattering_transform!(coeffs, st, signal)

In-place scattering transform. Fills pre-allocated S1/S2, returns the coefficients with S0 updated. Allocation-free; only allocates a new wrapper struct when S0 is a scalar (immutable).

scattering_transform!(coeffs, backend, st, signal)

Transform one signal on an explicit execution backend. SerialBackend runs the cascade in this task; ThreadedBackend (OhMyThreads extension) spreads the first-order wavelet groups across tasks, which is the only parallel axis available when there is a single field rather than a batch. The input transform is done once up front, so only the group loop is parallel.

scattering_transform2d!(coeffs, st, image)

In-place 2D scattering transform. Zero allocations for S1/S2 (buffers reused).

scattering_transform2d!(coeffs, backend, st, image)

Transform one image on an explicit execution backend. SerialBackend runs the cascade in this task; ThreadedBackend (OhMyThreads extension) spreads the first-order wavelet groups across tasks, which is the only parallel axis available for a single image.

scattering_transform3d!(coeffs, st, volume) -> coeffs

In-place 3D volumetric scattering transform; fills coeffs (a scales×orientations container) and returns it with S0 updated.

scattering_transform3d!(coeffs, backend, st, volume)

Transform one volume on an explicit execution backend — see the 2D counterpart.

scattering(st, x) -> ScatteringCoefficients

Non-mutating, allocation-tolerant, element-type-generic scattering transform — the autodiff-friendly counterpart of the in-place callable st(x). It composes the non-mutating Plans.forward_transform/Plans.inverse_transform with broadcast modulus and mean (no preallocated workspace, no mul!), so gradients flow through it via DifferentiationInterface (Mooncake/Zygote/Enzyme) and it accepts Dual/Float32 inputs. It returns the same coefficient container as st(x) and matches it numerically. Methods are added for the 1D/2D/3D transforms in their respective submodules.

Use st(x) (mutating, zero-alloc) for production forward passes; use scattering(st, x) when you need to differentiate the forward map (e.g. gradient-descent synthesis).

scattered_planar_scattering(x, y, ms, J; L=8, max_order=2, T=Float64,
                            spectral=SpectralBackends.AutoSpectralBackend(), period=nothing,
                            solve=false, weights=nothing, eps=nothing, maxiter=100,
                            rtol=Plans.default_solver_rtol(T, spectral, eps), damp=0,
                            nufft_nthreads=0)

Build a 2D planar scattering transform for a scalar field sampled at scattered points (x, y), using the same oriented Morlet wavelet bank as the gridded [ScatteringTransform2D] but computing the wavelet convolutions on a uniform Fourier mode grid of size ms = (m1, m2) via a nonuniform DFT: analysis maps the scattered points to the mode grid, the wavelet multiply happens there, and synthesis evaluates the filtered field back at the points. Apply it to a length-M vector of samples.

spectral selects the transform: SpectralBackends.DirectSumSpectralBackend is the in-core, dependency-free exact NUDFT (always available, O(M·prod(ms))); Plans.FINUFFTBackend and Plans.NonuniformFFTsBackend select a specific fast library (using FINUFFT / using NonuniformFFTs); SpectralBackends.NUFFTSpectralBackend takes whichever of those is loaded; and SpectralBackends.AutoSpectralBackend (the default) picks a fast library if one is loaded, else the direct sum. period is the physical domain size per axis (the Fourier period); it defaults so a uniform 0:m-1 grid reproduces the gridded FFT transform exactly. solve=false uses the fast adjoint (type-1) — exact for adequately-sampled band-limited fields, approximate on gappy/irregular data; solve=true recovers the true band-limited coefficients by least squares (slower, needed for irregular sampling). weights (length M, summing to 1) sets the quadrature for the spatial mean; the default is the uniform sample mean. eps is the NUFFT tolerance (ignored by the exact direct sum).

The solve is LSMR (Plans.lsmr_solve!), which costs one type-2 and one type-1 per iteration and keeps both ‖b - Af‖ and ‖A†(b - Af)‖ monotone, so stopping at maxiter returns the best iterate reached rather than a diverged one. rtol defaults to the accuracy of the transform underneath it — sqrt(eps(T)) for the exact direct sum, ~10·eps for a fast library, whose type-1 and type-2 are adjoints only to their own tolerance (Plans.default_solver_rtol). damp is a Tikhonov λ minimising ‖Af - b‖² + λ²‖f‖², 0 by default; with prod(ms) > M the problem is underdetermined, construction warns, and LSMR returns the minimum-norm solution unless a λ is given.

maxiter is a budget, not a target: the Krylov process terminates exactly at prod(ms) steps, so a default of 100 finishes a well-sampled band-limited field in a handful of iterations but truncates a field with a large irreducible least-squares residual. Truncation is safe — both error measures are monotone, so the returned iterate is the best one reached — but it is not converged, and an unconverged iterate is only as reproducible as the transform under it. A multi-threaded NUFFT accumulates its spreading in an order that varies between runs, and an unconverged solve amplifies that by about cond(A)²: measured at M = 400, ms = (16, 16), repeating one solve moves the coefficients by 5·10⁻⁴ at maxiter = 100 and 2·10⁻⁸ once the budget lets it converge, while nufft_nthreads = 1 gives bitwise-identical repeats at any budget. Raise maxiter for an answer that does not depend on the thread count, or pin nufft_nthreads = 1 for exact reproducibility.

nufft_nthreads sets the fast library's own thread count, and is honoured wherever the plan goes, including the per-task copies a threaded backend builds. 0 (the default) leaves it to the library, which is worth 2–3× at M ≳ 2·10⁴ on a serial backend; a per-task copy defaults to one instead, since the tasks already have the cores. Pass 1 to keep the library single-threaded everywhere — that is also the plan that executes without allocating, because a multi-threaded FFTW plan runs its parallel loop on Julia tasks.

spherical_scattering(pts_theta, pts_phi, lmax, J; max_order=2,
                     spectral=SpectralBackends.AutoSpectralBackend(),
                     rtol=SphericalCore.default_rtol(spectral), maxiter=500)

Build a spherical scattering transform for a scalar field at scattered points (θ, φ) on S², using smooth difference-of-Gaussians band-pass wavelets. spectral selects the spherical-harmonic transform: SpectralBackends.DirectSumSpectralBackend is the in-core, dependency-free exact least-squares transform (always available, O(M·(lmax+1)²)); SpectralBackends.NUFSHTSpectralBackend uses the NUFSHT fast path (needs using NUFSHT); SpectralBackends.AutoSpectralBackend (the default) picks NUFSHT if its extension is loaded, else the direct transform. Accurate analysis needs the sampling to resolve the band limit, i.e. roughly M ≳ (lmax+1)² well-distributed points.

structured_spherical_scattering(lmax, J; max_order=2,
                                spectral=SpectralBackends.AutoSpectralBackend(), T=Float64,
                                rtol=1e-8, maxiter=500)

Build a spherical scattering transform for a scalar field sampled on the structured equiangular grid (Nθ = lmax+1 colatitudes, Nφ = 2lmax+1 longitudes), using the same smooth difference-of-Gaussians band-pass wavelets as spherical_scattering. Apply it to a (Nθ, Nφ) grid of samples; obtain the grid points with structured_sphere_points. spectral selects the transform: SpectralBackends.DirectSumSpectralBackend (in-core, dependency-free) by default, or SpectralBackends.FSHTSpectralBackend (the fast exact SHT, needs using FastSphericalHarmonics); SpectralBackends.AutoSpectralBackend picks the fast path if its extension is loaded, else the direct SHT.

structured_sphere_points(lmax) -> (Θ, Φ)

Colatitudes Θ (length lmax+1) and longitudes Φ (length 2lmax+1) of the equiangular structured grid used by structured_spherical_scattering; sample a field as [f(θ, φ) for θ in Θ, φ in Φ]. In-core (no dependency); matches FastSphericalHarmonics.sph_points.

ReconstructionWorkspace(st)

Scratch for the linear wavelet inverse and for phase retrieval: the complex wavelet coefficient fields (nw of them — that is the representation's own size), the low-pass field, and a spectrum accumulator.

No conjugated filters are held. A Morlet's Fourier response is real — its analyticity is the half-plane support, not a complex value — so conj(ψ̂_λ) = ψ̂_λ and the dual frame filter is the filter itself. The bank is used directly.

wavelet_transform(st, x) -> (; wavelet, lowpass)
wavelet_transform!(ws, st, x) -> (; wavelet, lowpass)

The linear (pre-modulus) wavelet layer underlying the scattering transform: the complex wavelet coefficient fields wavelet[λ] = x ⋆ ψ_λ (the continuous wavelet transform) and the low-pass field lowpass = x ⋆ φ. Together these are an exact, invertible representation of x (see iwavelet!); taking |wavelet[λ]| and averaging is the first scattering layer.

iwavelet(st, wavelet, lowpass) -> x
iwavelet(st, wt) -> x
iwavelet!(ws, st, wavelet, lowpass) -> x

Exact inverse of wavelet_transform! via the tight-frame conjugate-filter sum x̂ = Σ_λ ψ̂_λ^* · (x̂·ψ̂_λ) + φ̂^* · (x̂·φ̂). Returns the reconstructed real field. With the tight-frame bank (Σ|ψ̂_λ|²+|φ̂|² ≡ 1) this satisfies iwavelet(st, wavelet_transform(st, x)...) ≈ x to machine precision.

The accumulation is in place: the spectrum sum is built in one buffer rather than rebuilt per wavelet, which is what makes iters rounds of phase retrieval affordable.

reconstruct_phase(st, moduli; iters=200, init=nothing, seed_lowpass=nothing) -> x

Phase retrieval from the first-order moduli moduli[λ] = |x ⋆ ψ_λ| (the real fields the scattering transform averages), via Gerchberg–Saxton alternating projections:

1. take the linear wavelet transform of the current estimate;
2. re-impose the target magnitudes on each band-pass channel (keeping the recovered phase);
3. reconstruct with the exact frame inverse iwavelet!; repeat.

The low-pass channel carries no magnitude target, so it is taken from the current estimate each iteration (or from seed_lowpass if supplied). The reconstruction is determined only up to a global sign/phase. init (defaults to a random field matched in energy to the moduli) seeds the estimate; pass one for reproducibility. Returns the real reconstructed field.

One workspace is built up front and reused across all iters, so the loop allocates nothing.

synthesize(st, target; backend, init=nothing, iters=500, lr=0.05, loss=scattering_loss) -> (; field, losses)

Reconstruct a field whose scattering coefficients match target by gradient descent (Bruna & Mallat microcanonical synthesis): starting from init (random by default), minimize loss(scattering(st, x̂), target) with Adam. target may be a precomputed coefficient container or a field (its coefficients are taken first). The gradient is obtained through DifferentiationInterface, so backend is any ADTypes backend; the synthesized result is a sample with matching statistics, not the exact original (the modulus discards local phase). Requires using DifferentiationInterface and an AD backend package.

With Enzyme, pass AutoEnzyme(; mode = Enzyme.set_runtime_activity(Enzyme.Reverse)). The cascade maps a closure over the filter bank, so each closure captures constant filters alongside the active input and Enzyme's static activity analysis cannot separate the two — runtime activity is its documented remedy for that, and a plain AutoEnzyme() raises EnzymeRuntimeActivityError instead of a gradient.

scattering_loss(c, target) -> Real

Default synthesize objective: the normalized squared error between the coefficient container c and the target, summed over the first- and (when present) second-order coefficients, ‖S₁(c)−S₁(t)‖² + ‖S₂(c)−S₂(t)‖² divided by the target energy. Differentiable in c, so it composes with scattering(st, ·) under autodiff.

MonogenicScattering{D,T,FB,Tree,P,G,CB,MW}

Monogenic scattering transform on a D-dimensional grid: an isotropic MonogenicFilterBank, the scattering path tree (strictly-increasing scale), a spectral plan, and the periodized MonogenicWorkspace the cascade runs in.

MonogenicFilterBank{D,T,A,W,R,MV}

Isotropic band-pass wavelets wavelets (one per scale/sub-octave), the D scale-free Riesz multipliers riesz, and the complementary low-pass averaging, forming a tight frame. Every container is a type parameter.

The wavelets and the low-pass are real: a radial band-pass is a real function of |k|. Only the Riesz multipliers R_d(k) = -i k_d/|k| are complex.

ComputedMonogenicFilterBank{D,T,A,R,MV}

The monogenic bank with its radial band-passes evaluated on demand instead of stored — the counterpart of FilterBanks.ComputedFilterBank2D. The Riesz multipliers stay stored: there are D of them however many scales the bank has, so they do not grow with J.

MonogenicWorkspace{T,D}

The periodized cascade plus what the monogenic amplitude needs on top of it: one real accumulator per resolution, and the Riesz-weighted band-passes R_d ψ_j periodized to each parent resolution.

R_d ψ_j is stored as one periodized filter rather than two, because periodizing the factors separately is a different function from periodizing their product. Level 1 stores nothing: there r = 1, so the two full-resolution filters fuse directly in the multiply (ScatteringCore.periodize_mul2!).

build_monogenic_bank([T=Float64,] dims::NTuple{D,Int}, J; Q=1) -> MonogenicFilterBank

Build a D-dimensional isotropic Morlet-style monogenic filter bank: J octaves × Q sub-octaves of radial band-pass wavelets (center frequency ξ_j = ξ₀·2^{-j/Q}, widths from the Lostanlen/Kymatio rule, reusing Filters.Morlet1D for the radial profile), the Riesz multipliers, and the tight-frame complementary low-pass.

riesz_multipliers(dims, ::Type{T}=Float64) -> NTuple{D, Array{Complex{T},D}}

The D Riesz-transform frequency multipliers R_d(k) = -i k_d/|k| (with R_d(0)=0) over a grid of size dims. Scale-free (a ratio of frequencies), so one set serves every wavelet scale. They satisfy Σ_d |R_d(k)|² = 1 off the DC bin.

monogenic_amplitude(m0, riesz_components) -> A

Monogenic amplitude A = √(m₀² + Σ_d m_d²) from the band-pass field m0 and the tuple/vector of Riesz component fields. Broadcast, so it is CPU/GPU/autodiff-generic.

monogenic_components(st, x, j) -> (; bandpass, riesz, amplitude, phase)

The monogenic decomposition of x band-passed by isotropic wavelet j (1-based): the band-pass field bandpass = x ⋆ ψ_j, the D Riesz component fields riesz, the monogenic amplitude √(bandpass² + Σ|riesz|²), and the local monogenic phase = atan(‖riesz‖, bandpass). The Riesz vector's direction gives the local orientation (e.g. atan(riesz[2], riesz[1]) in 2D).

spherical_monogenic_scattering(pts_theta, pts_phi, lmax, J; max_order=2)

Build a monogenic spherical scattering transform for a scalar field at scattered points (θ, φ) on S². The nonlinearity is the spherical monogenic amplitude A_j = √(U⁰_j² + |U^R_j|²), where U⁰_j is the difference-of-Gaussians band-pass and U^R_j is the spin-1 Riesz field R = ð∘(−Δ_S)^{-1/2}. The Riesz energy |U^R_j|² = |∇_S g_j|² (with g_j = (−Δ_S)^{-1/2} U⁰_j) is evaluated using only spin-0 spherical-harmonic transforms via the identity |∇_S g|² = ½ Δ_S(g²) − g Δ_S g, so no spin-weighted synthesis is required. spectral selects the spherical-harmonic transform as in spherical_scattering (dependency-free direct SH transform by default, NUFSHT fast path when loaded).

spherical_monogenic_components(st, field, j) -> (; bandpass, riesz, amplitude, phase, orientation)

Pointwise spherical monogenic decomposition of field band-passed at scale j — the S² analogue of the planar monogenic_components. Returns the band-pass field bandpass = U⁰_j, the spin-1 Riesz tangent vector riesz = (u_θ, u_φ) (U^R_j = ð∘(−Δ_S)^{-1/2} U⁰_j), the monogenic amplitude √(U⁰² + ‖U^R‖²), the phase = atan(‖U^R‖, U⁰), and the local orientation = atan(u_φ, u_θ) of the Riesz vector. st is a spherical_monogenic_scattering transform. With the in-core direct SH plan (dependency-free) the Riesz field is synthesized as the surface gradient of g = (−Δ_S)^{-1/2}U⁰; with a NUFSHT-backed plan it uses spin-weighted synthesis — the two agree to solver accuracy.

structured_spherical_monogenic_scattering(lmax, J; max_order=2,
                                          spectral=SpectralBackends.AutoSpectralBackend(),
                                          T=Float64, rtol=1e-8, maxiter=500)

Structured-grid counterpart of spherical_monogenic_scattering. spectral selects the SH transform as in structured_spherical_scattering (dependency-free direct SHT by default, fast exact SHT with using FastSphericalHarmonics).

scattering_field(st, x; subsample) -> ScatteringField{1,2}D

Localized (Mallat) scattering transform: returns the per-path low-passed, subsampled fields. Methods are added by the per-dimension transform modules.

scattering_field!(field, st, x) -> field

In-place localized scattering transform into a pre-allocated ScatteringField.

ScatteringField1D{T,A,Tree,WS}

Localized 1D scattering field. data is (M, npaths); column p is the localized field of path p at the subsampled resolution M = N ÷ s. ws is the Cascade.FieldWorkspace that fills it.

ScatteringField2D{T,A,Tree,WS}

Localized 2D scattering field. data is (My, Mx, npaths); slice [:, :, p] is the localized field of path p at subsampled resolution (My, Mx) = (Ny, Nx) ÷ s.

path_field(sf, p) -> view

Non-allocating view of path p's localized field (a vector for 1D, a matrix for 2D).

subsample_factor(sf) -> Int

The decimation factor s applied to produce the field resolution.

ScatteringCoefficients1D{T,V,M,S0}

Immutable container for 1D scattering coefficients. S0 can be scalar T (return new struct) or mutable container (update in place). Uses multiple dispatch for optimal S0 handling.

Type Parameters

• T: Element type
• V: 1D array type
• M: 2D array type
• S0: S0 storage type (T for scalar, AbstractVector{T} for mutable)

ScatteringCoefficients2D{T,V,M,S0}

Immutable container for 2D planar scattering coefficients. S0 can be scalar T or mutable container - dispatch handles both optimally.

zeroth_order(c) -> Real

Extract the zeroth-order (average / S0) scattering coefficient.

first_order(c) -> AbstractVector

Extract the first-order (S1) scattering coefficients.

second_order(c) -> AbstractMatrix

Extract the second-order (S2) scattering coefficients.

flatten1d(coeffs::ScatteringCoefficients1D{T}) -> Vector{T}

Flatten to vector: [S0; S1; vec(S2 upper triangular)]. Only includes unique S2 elements where j2 > j1 (saves ~50% space).

flatten2d(coeffs::ScatteringCoefficients2D{T}) -> Vector{T}

Flatten to vector: [S0; S1; vec(S2 upper triangular)].

flatten1d!(out, c) -> out

Zero-allocation flatten into a pre-allocated vector of length flatten_length(c).

flatten2d!(out, c) -> out

Zero-allocation flatten into a pre-allocated vector of length flatten_length(c).

flatten_length(c) -> Int

Length of the flattened coefficient vector [S0; S1; vec(S2 upper triangle)].

flat_length(n) -> Int
flat_row_s0() -> Int
flat_row_s1(j, n) -> Int
flat_row_s2(j1, j2, n) -> Int

Row layout of the flattened coefficient vector [S0; S1; vec(S2 upper triangle)] for n wavelets, as produced by flatten1d!/flatten2d!.

The batched paths write straight into flattened columns rather than filling a coefficient container first, so they need the layout as arithmetic. Defining it here keeps the one definition that flatten*! also walks — a second copy elsewhere would silently drift.

flat_eltype(c) -> Type

Element type a flattened column of c needs.

S1 and S2 are moduli and stay real whatever the input is, but S0 is the field mean and is complex for a complex field, so the flat layout is as wide as the two together.

update_S0(c, val)

Update the zeroth-order coefficient storage with val and return the coefficients.

compute_shape_sparsity(S1, S2, meta) -> (; sparsity, shape)

Reduced second-order descriptors (in the spirit of the reduced wavelet scattering transform, Allys et al. 2019; Cheng & Ménard 2021), as J × J matrices over scale pairs (j1, j2) with j2 > j1:

• sparsity (s₂₁): the orientation-averaged ratio ⟨S₂ / S₁⟩ — how much energy cascades from scale j1 to the coarser scale j2 (large for sparse/intermittent fields).
• shape (s₂₂): the anisotropy of the cascade — the normalized second angular harmonic ⟨S₂ · cos(2 Δθ)⟩ / ⟨S₂⟩ over orientation pairs, where Δθ = θ₂ − θ₁. It is ≈ 0 for statistically isotropic fields and departs from zero when the field has oriented structure.

normalized_coefficients(c) -> (; S0, s1, s2)

Amplitude-normalized coefficients: s1[j] = S1[j]/S0 and s2[j1,j2] = S2[j1,j2]/S1[j1] (zero where S1[j1] == 0).

log_coefficients(c; pad=eps) -> (; S0, logS1, logS2)

Log-coefficients log(S1 + pad), log(S2 + pad) (the small pad keeps structural zeros finite). Useful to gaussianize the heavy-tailed coefficients of intermittent fields.

scattering_batch(st::Scattering1D.ScatteringTransform1D, X) -> Matrix

Apply a 1D scattering transform to a batch of signals X of size (N, B) (signals as columns), returning a (flatten_length, B) matrix whose column b is flatten1d(st(X[:, b])). The plan and all workspace buffers are reused across the batch (only the small scalar-S0 wrapper is re-allocated per column).

scattering_batch(st::Scattering2D.ScatteringTransform2D, X) -> Matrix

Apply a 2D scattering transform to a batch of images X of size (Ny, Nx, B), returning a (flatten_length, B) matrix whose column b is flatten2d(st(X[:, :, b])). Plan and workspace are reused across the batch.

scattering_batch(st::Scattering3D.ScatteringTransform3D, X) -> Matrix

Apply a 3D scattering transform to a batch of volumes X of size (Nz, Ny, Nx, B), returning a (flatten_length, B) matrix. Plan and workspace are reused across the batch.

scattering_batch(st::SphericalCore.SphericalScattering, X) -> Matrix
scattering_batch(st::SphericalCore.SphericalMonogenicScattering, X) -> Matrix

Transform a stack of B fields sampled by the same spherical plan — X of size (field_size…, B), so (M, B) for a scattered point set and (nθ, nφ, B) on a structured grid.

Rows follow the flat layout of Coefficients.flat_length and its row accessors. The spherical cascade pairs each scale with strictly coarser ones (j2 < j1), so a pair lands in the row that layout assigns to the unordered pair, flat_row_s2(j2, j1, J).

scattering_batch!(out, st, X; workspace = nothing) -> out

In-place counterpart of scattering_batch: write the flattened coefficients of each slice of X into the preallocated (flatten_length, B) matrix out. Backend-dispatched ! methods (scattering_batch!(out, backend, st, X)) are added by the corresponding extensions.

The default transforms one slice at a time against the transform's own single-slice plan. Passing a batch_workspace instead runs the whole stack through one batched plan (Batched.batch_cascade!).

Per-slice is the default because it is faster on a CPU at every batch size measured — 1.2–1.5× serial and 4.5× threaded. The cascade is memory-bandwidth bound and issues O(nw + paths) operations against the same data, so what governs is whether that data stays resident: a slice does, a B-slice stack does not. The batched plan's amortised per-call overhead does not recover the difference. It wins on a device, where the batch is what fills the machine — which is why the GPU extension builds one.

batch_coeffs(st, T, S0T = T) -> coefficient container

A coefficient container for st whose S0 is a 1-element array, so a loop over a batch writes every field in place instead of rebuilding the struct once per slice.

S0T is given separately because S1/S2 are moduli and stay real whatever the input is, while S0 is the field mean and is complex for a complex field.

flat_rows(st) -> Int

Rows a flattened coefficient column of st occupies — the height of scattering_batch's output.

batch_workspace(st, B; spectral = ..., fft_nthreads = 1) -> Batched.BatchWorkspace

Build the batched workspace for transform st at batch size B: a spectral plan over the whole (spatial…, B) stack plus the scratch Batched.batch_cascade! runs in.

The plan is rebuilt because a batched transform is a different plan from a single-slice one — that is exactly what makes it one execution instead of B. Build it once and pass it back in through scattering_batch! when transforming many stacks of the same size.

fft_nthreads is the FFT library's own thread count, baked into the plan. It defaults to 1: the cascade is memory-bandwidth bound and the transform is only ~a quarter of a cascade step, so threading inside the FFT is capped near 1.3× by Amdahl no matter how many cores it gets, while threading over the batch or wavelet axis parallelises the whole step. Measured on 8 threads it returns 1.08–1.17× on large stacks and 0.33× on small ones, where the library's per-execution task spawn (which also allocates) exceeds the transform itself. Raise it only for a single stack large enough to be worth it, transformed with no parallelism above.

BatchWorkspace{T,D}

Preallocated state for repeated batched transforms at a fixed batch size B: a periodized cascade workspace whose arrays all carry a trailing stack axis, the real staging and reduction buffers, and the precomputed order-2 groups. Once built, a batched transform does no data-proportional allocation.

red is the spatial-reduction target, shaped (1…, 1, B) so sum! contracts every axis but the batch; that one shape serves every resolution.

invN is per level. A single 1/∏spatial would leave every decimated row short by exactly ∏r — a clean factor that no batched-versus-batched comparison can see.

batch_cascade!(out, ws, X) -> out

Write the flattened scattering coefficients of every slice of X into the columns of out.

out has Coefficients.flat_length(n) rows; row assignment goes through Coefficients.flat_row_*, the same layout flatten1d!/flatten2d! produce, so batched and per-slice results are interchangeable.

The batched form of Cascade.cascade! — same periodization and the same resolutions, with the scalar mean replaced by a per-slice reduction over the stack. Each decimation tuple gains a trailing 1 so the stack axis is never folded, and the filters carry a trailing singleton so they broadcast across it.

decimation(N, scale, α) -> NTuple{D,Int}

Per-axis decimation for a wavelet of octave scale: 2^max(scale-α, 0), reduced on each axis to the largest power of two that divides N[d].

Per axis because one scalar cannot serve every grid — (128, 100) admits (8, 4) but no scalar 8. Clamping rather than throwing keeps every N constructible; the clamp only ever reduces decimation, so it costs speed, never accuracy.

PeriodizedWorkspace{T,D}

Everything the cascade runs in: one working array and one spectrum array per live decimation, the spectral plan at each of those sizes, and the filters periodized to each resolution they are applied at.

Everything per-resolution is a vector indexed by level, with level[j] the level of wavelet j and level 1 always the full grid. Levels rather than resolution tuples because the cascade looks these up once per wavelet and once per order-2 path, and an integer index costs nothing where hashing a tuple would.

out[l] is work[l] itself when the plan inverts in place, so a resolution costs one array rather than two. spec[l] is empty except at levels that are ever a parent, where it holds Û₁ across that wavelet's children.

A/FA are not pinned to rank D: D is the spatial rank, which is what the resolutions are in, while a batched workspace's arrays carry a trailing stack axis.

build(fb, groups, dims, T, α, spectral) -> PeriodizedWorkspace

groups is the cascade work list (j1, children, pathids); it determines which resolutions are live and which (wavelet, resolution) filter pairs are ever applied, so nothing unused is built.

cascade!(S1, S2, ws, fb, groups, xfft) -> (S1, S2)

Both scattering orders in one grouped pass, each convolution produced directly on the grid its wavelet's band actually needs.

Order 1 multiplies at full resolution and periodizes the product, which is why the full-resolution filter is the one order 1 applies. Order 2 multiplies Û₁ by ψ_{j₂} already periodized to the parent's grid, then periodizes that product onto the child's. r₂ ≥ r₁ always, because admissibility makes the child strictly coarser.

group!(S1, S2, ws, fb, g, xfft) -> nothing

One first-order wavelet and all of its children. Split out because it is also the unit of parallelism: group j1 writes only S1[j1] and S2[j1, :], so tasks over groups never collide — provided each task holds its own task_copy of the workspace, whose buffers are mutable.

S2 is not cleared here; the caller does that once.

scattering_values(st, x) -> (S0, S1, S2)

The same periodized cascade as cascade!, written with allocating transforms and broadcasts so that reverse-mode backends can differentiate it. Dimension-generic; the per-dimension scattering methods wrap the result in their coefficient container.

Periodizing here is a contract, not an optimisation: synthesize differentiates this while the user reads coefficients from st(x), so at oversampling < J an unperiodized version would optimise a different function than the one being reported.

task_copy(ws, fb) -> ws′

A workspace usable concurrently with ws: fresh buffers and task-local plans, sharing the periodized filters, which are read-only. fb is the filter bank the copy will be run against.

The buffers are mutable state, so sharing a PeriodizedWorkspace across tasks is a data race — every task that runs group! needs its own.

fb is required because a computed bank's level-1 entries are all one shared scratch array. A task given its own bank (FilterBanks.task_bank) must also get a level-1 list pointing at that bank's scratch; sharing the original's would have every task refill its own array and then read another's.

with_plans(f, ws) -> ws′

ws with every plan replaced by f(plan), sharing all buffers and filters.

For instrumenting or retargeting the transforms a cascade issues without rebuilding the workspace. f must preserve Plans.inplace_inverse, which is what the shared out/work aliasing was decided from.

FieldWorkspace{T,D}

The periodized cascade plus what the localized transform S_p = (|U_p| ⋆ φ_J) ↓ s needs on top of it: φ_J periodized to every level a modulus is produced at, and one buffer and plan at the output resolution.

Built with alias_out = false, because unlike the coefficient cascade this one forward-transforms every modulus, a leaf's included.

Decimation is capped at s: a path produced coarser than the grid it is written to would have to be upsampled. The cap costs nothing at the default s = 2^(J-1), which already exceeds every wavelet's own decimation.

build_field(fb, groups, dims, T, α, s, φ; alloc, makeplan) -> FieldWorkspace

Workspace for the localized transform at output subsample factor s. φ is the full-resolution low-pass.

field_workspace(st, subsample) -> FieldWorkspace

The localized-field workspace for transform st at output factor subsample.

Duck-typed on st (filter_bank, groups, dims, plan, buffer_input, pw), so 1D, 2D and device-resident transforms all go through it: buffer_input supplies the array type and plan the plan kind, which is all that distinguishes a device build from a host one.

field_oversampling(fb, α, s) -> Int

The smallest oversampling at least α under which no wavelet decimates past s.

r = 2^max(j-α, 0) ≤ s = 2^q iff α ≥ j - q, so raising α to maxscale - q caps every resolution at the output grid without changing the dyadic family any wavelet lands on.

localize!(dst, fw, spec, level) -> dst

dst = real(ifft(spec ⋅ φ̂) ↓ (rout/r)), where spec is the spectrum of a modulus living at level.

Fused: periodizing the product onto the output grid is the decimated inverse transform, so the inverse runs at output size.

field_group!(data, fw, fb, g, xfft, p1_first) -> nothing

One first-order wavelet and all of its children, written as localized fields into data's trailing path axis. The order-1 path id is p1_first + j1 - 1; the children's come from the group.

field_cascade!(data, fw, fb, groups, xfft, root, p1_first) -> data

Every path's localized field. root is the order-0 path id, p1_first the first order-1 id.

periodize_mul!(dst, src, filt, r) -> dst
periodize!(dst, src, r) -> dst

dst[k] = (1/∏r) Σ_m src[k + m∘n] * filt[k + m∘n], where n = size(dst) and m runs over ∏r alias blocks — the spectrum of (src ⋆ filt) decimated by r per axis.

This is the whole point of the periodized cascade. For any spectrum Ŷ, with no band-limit assumption whatsoever,

ifft_N(Ŷ)[1:r:end] == ifft_M(periodize_r(Ŷ)),   M = N ./ r

so the decimated samples of a full-size inverse transform are the small inverse transform of the periodized spectrum. Transforming at full size and keeping every r-th sample is wasted work, not extra accuracy. The 1/∏r is the decimation's own factor: an inverse normalised by 1/length on a grid ∏r times smaller supplies ∏r too little.

r is per axis because a single scalar cannot serve every grid — (128, 100) admits r = (8, 4) but no scalar 8. It also fits oriented wavelets, whose support is anisotropic.

Indices are raw and 0-based, which is what makes the alias stride plain addition: these plans are unshifted, so bin k of the small grid gathers k, k+n, k+2n, … of the large one. The multiply is fused in and the 1/∏r folded into each block, so the whole thing is one pass over src and one write of dst per alias — broadcast throughout, so it carries to a device unchanged.

periodize_mul2!(dst, src, f1, f2, r) -> dst

periodize_mul! with two filters: dst[k] = (1/∏r) Σ_m src[·] * f1[·] * f2[·].

The monogenic amplitude band-passes by R_d · ψ_j, a product of two full-resolution filters. Forming that product into scratch first and then periodizing would sweep the full grid twice; fusing it costs one pass, exactly as the single-filter form does.

_modulus(z)

|z| for the scattering nonlinearity, as sqrt(abs2(z)) rather than abs(z).

Base.abs(::Complex) is hypot, which rescales by the larger component to stay exact across the whole exponent range. That costs a branch and a division per element, and the modulus is 34-46% of a cascade's runtime: hypot measures 5.1x (N=1024) to 9.2x (N=256) slower than the direct form, worth 1.22x (N=1024) to 1.62x (N=256) on the whole 2D cascade.

The two agree to one ulp — 2.2e-16 relative over samples spanning 1.7e-12 to 2.7e+11. The range hypot buys back is not reachable here: abs2 overflows only for |z| > sqrt(floatmax) (1.3e154 in Float64, 1.8e19 in Float32) and underflows below sqrt(floatmin), while a wavelet coefficient is O(‖x‖). A field near those magnitudes must be rescaled before transforming.

periodize(src, r) -> dst

Allocating periodize_mul!, for the non-mutating autodiff forward: no dst to write into and no in-place writes for a reverse-mode backend to trip over.

periodize_filter!(dst, src, r) -> dst

A filter periodized onto a coarser grid: dst[k] = Σ_m src[k + m∘n], with no 1/∏r.

Not interchangeable with rebuilding the filter on the coarse grid. The Morlet profile is a function of normalised angular frequency — kx = 2π·fftfreq(n,i-1) and k0 = 3π/(4·2^j) are both independent of n — so evaluating it on a grid of size n = N/r stretches the frequency axis by r and lands the peak at bin (N/r)·k0/2π instead of N·k0/2π: one octave too low per factor of two, i.e. a different wavelet. Periodizing the original filter is what makes the undecimated case reproduce the full cascade exactly.

No scale factor here because the 1/∏r belongs to the product being decimated, not to the filter; applying it in both places would double-count.

gaussian_lowpass!(φ, σ) -> φ
gaussian_lowpass(T, dims, J; sigma0=0.8) -> Array{T,D}

Isotropic Gaussian low-pass φ̂(k) = exp(-|k|²σ²/2) on the FFT grid of dims, with k the angular frequency (2π·fftfreq per axis) and σ = sigma0·2^J a real-space width. The scaling function of the localized (Mallat) transform S_p = (U_p ⋆ φ_J) ↓ s.

Two properties are load-bearing. φ̂(0) = 1 makes the convolution preserve a field's mean; and φ̂ must vanish on the subsampling lattice {m·N/s, m ≠ 0}, since that is exactly the condition for ⟨S_p⟩ to survive decimation and still equal the scattering coefficient ⟨U_p⟩. The lattice starts at |k| = 2π/s, so at the default s = 2^(J-1) the largest surviving term is exp(-(4π·sigma0)²/2) — 1e-22 at sigma0 = 0.8. That condition is the requirement; the Gaussian family and sigma0 are a choice with wide margin.

Distinct from a filter bank's averaging, the tight-frame complement √(max(0, 1-Σ|ψ|²)): every ψ̂ here is analytic, so that one is identically 1 across the analytic complement and fails the lattice condition outright.

plan_like(plan, prototype) -> plan′

A plan of the same kind as plan, sized for arrays like prototype.

Sizing from an array rather than from dimensions plus a backend tag is what lets one piece of code plan on host and device alike: a device plan has no spectral-backend tag to look up — spectral_backend deliberately throws on it — but it does have arrays, and AbstractFFTs dispatches on their type. Used where a cascade needs plans at resolutions not known when the transform was built.

inplace_inverse(plan) -> Bool

Whether inverse_transform!(x, plan, x) is valid. When it is, a transform points its convolution output at its multiply scratch and owns one fewer field-sized array; when it is not, the two stay distinct. Default false: the direct-sum plan alternates through its own scratch and aliasing would corrupt an odd number of axes.

FilterBank1D{T,V,W,MV}

Complete 1D filter bank for the scattering transform. Every container is a type parameter (no hardcoded Vector): V the per-filter array type (CPU/GPU/static/…), W the wavelet collection, MV the metadata collection.

Fields

• wavelets::W: wavelet filters in the Fourier domain (W<:AbstractVector{V})
• averaging::V: low-pass averaging (scaling) filter
• meta::MV: per-wavelet WaveletMeta
• J::Int: number of octaves (scales)
• Q::Int: wavelets per octave

FilterBank2D{T,M,W,MV}

Complete 2D filter bank with oriented wavelets. Containers are type parameters (no hardcoded Vector): M the per-filter matrix type, W the wavelet collection, MV the metadata collection.

Fields

• wavelets::W: oriented wavelet filters (W<:AbstractVector{M})
• averaging::M: low-pass averaging filter
• meta::MV: per-wavelet WaveletMeta
• J::Int: number of scales
• L::Int: number of orientations

FilterBank3D{T,A<:AbstractArray{Complex{T},3}}

Complete 3D oriented Morlet filter bank: J scales × n_orient sphere directions, plus a low-pass averaging filter.

WaveletMeta{T}

Concrete per-wavelet metadata: a struct rather than a NamedTuple, so the container stays concretely typed.

Fields

• scale::Int: octave index j
• q::Int: sub-octave index within the octave (1D, 0..Q-1); 0 for 2D
• orient::Int: orientation index l (2D, 0..L-1); 0 for 1D
• j_eff::T: effective log-scale used to order paths. j + q/Q in 1D, T(j) in 2D. The second-order admissibility constraint is j_eff(child) > j_eff(parent) (frequency strictly decreasing) — which for 2D means scale strictly increasing over all orientation pairs.
• center_freq::T: wavelet center frequency
• theta::T: orientation angle in radians (2D); 0 for 1D

build_filter_bank1d(N::Int, J::Int; Q::Int=1) -> FilterBank1D

Build a 1D Morlet filter bank with dyadic scales.

Arguments

• N::Int: Signal length (FFT size)
• J::Int: Maximum scale (number of octaves)
• Q::Int: Wavelets per octave (default 1 for dyadic, 8 for high Q)

Returns

• FilterBank1D: Complete filter bank with J scales

build_filter_bank1d(T, N, J; Q=1, cache=true)
build_filter_bank3d(T, N, J; n_orient=6, cache=true)

cache=false evaluates each wavelet on demand instead of storing the bank.

build_filter_bank2d(N::NTuple{2,Int}, J::Int; L::Int=8) -> FilterBank2D

Build a 2D oriented Morlet filter bank.

Arguments

• N::NTuple{2,Int}: Image dimensions (Ny, Nx)
• J::Int: Number of dyadic scales
• L::Int: Number of orientations (default 8, evenly spaced)

Returns

• FilterBank2D: Complete 2D filter bank

build_filter_bank2d(T, N, J; L=8, cache=true)

cache=false returns a ComputedFilterBank2D, which evaluates each wavelet on demand instead of holding the whole bank. Slower per multiply, and the only way a large grid fits.

build_filter_bank3d(N::NTuple{3,Int}, J::Int; n_orient::Int=6, T=Float64) -> FilterBank3D

Build a 3D oriented Morlet filter bank with J dyadic scales and n_orient near-uniform orientations on the sphere (Fibonacci spiral).

ComputedFilterBank1D{T,V,MV}
ComputedFilterBank3D{T,A,MV}

The 1D and 3D counterparts of ComputedFilterBank2D; the saving grows with dimension, since the bank is J·n_orient+1 arrays of the grid either way.

ComputedFilterBank2D{T,M,MV}

A 2D bank that keeps its wavelets' parameters rather than their samples, and evaluates one into a single scratch array when it is asked for.

A bank of J·L wavelets is J·L+1 arrays of the grid's size, and it is the largest object in a transform: 1056 MiB at 2048², J=4, L=8, Float64, against 320 MiB for the whole cascade's working set. Evaluating instead of storing takes that to two arrays — the low-pass and the scratch — at the cost of two exp per grid point per application, measured 3.2x (N=1024) to 3.8x (N=2048) on the multiply, which is 9-14% of a cascade's time.

The low-pass has to be stored: it is sqrt(1 - Σⱼ|ψⱼ|²), so recomputing it would mean evaluating the whole bank. The tight-frame rescale is a single scalar and is applied on materialisation.

The scratch is mutable state, so unlike a stored bank this one cannot be shared between tasks — task_bank gives a task its own, which is one array rather than J·L.

nwavelets(fb) -> Int

How many wavelets the bank holds.

filter_at(fb, j) -> AbstractArray

The bank's j-th wavelet in the Fourier domain.

For a stored bank this is an index. For a computed one it writes into shared scratch, so the result is valid only until the next call — every cascade here finishes one filter's multiply before asking for the next.

iscomputed(fb) -> Bool

Whether filter_at returns shared scratch rather than a stored array. A caller that needs two filters live at once, or a filter at a resolution the bank does not hold, must materialise.

task_bank(fb) -> fb′

A bank usable concurrently with fb. Stored banks are read-only and are returned as they are; a computed bank gets its own scratch, sharing everything else.

batch_views(fb, spatial) -> Vector

The bank's filters reshaped to (spatial…, 1) so they broadcast over a batch axis, built once so the batched cascade never reshapes in its inner loop.

A computed bank returns the same view nwavelets times — one view of the one scratch array. That is only sound because filter_at refills the scratch immediately before each use, which is what Batched._filter does.

Morlet1D{T<:Real}

1D Morlet wavelet in frequency domain, over normalized frequency ω ∈ [0, ½]:

Ψ(ω) = exp(-(ω-ξ)²/2σ²) - κ exp(-ω²/2σ²),   κ = exp(-(ξ/σ)²/2)

σ here is a frequency-domain width (unlike Morlet2D/Morlet3D, whose σ is a real-space width), so the zero-mean term is κ, evaluated in frequency_response! — Ψ(0) = 1·κ - κ·1 after normalising, which is what makes the wavelet admissible. The filter is analytic: zero for ω < 0.

Type Parameters

• T: Element type (Float32, Float64, etc.)

Fields

• center_freq::T: Center frequency ξ
• bandwidth::T: Frequency-domain width σ
• N::Int: Filter length (FFT size)

Morlet2D{T<:Real}

2D oriented Morlet wavelet in frequency domain.

The 2D Morlet wavelet is created by taking a 1D Morlet and rotating it to angle θ, with elliptical Gaussian envelope controlled by elongation.

Type Parameters

• T: Element type (Float32, Float64, etc.)

Fields

• center_freq::T: Center wavenumber |k₀|
• bandwidth_x::T: Bandwidth along major axis
• bandwidth_y::T: Bandwidth along minor axis (controls elongation)
• theta::T: Orientation angle in radians
• beta::T: Correction factor
• N::NTuple{2,Int}: Filter dimensions (Ny, Nx)

Morlet3D{T<:Real}

3D oriented Morlet wavelet in the frequency domain, a bump centered at k₀ n̂ for a unit direction n̂ on the sphere, with an anisotropic Gaussian envelope (std σ∥ along n̂, σ⊥ = σ∥/elongation perpendicular) and analytic on the half-space k·n̂ ≥ 0.

Fields

• center_freq::T: |k₀|
• sigma_par::T, sigma_perp::T: real-space envelope widths along / perpendicular to n̂
• direction::NTuple{3,T}: unit orientation n̂
• beta::T: zero-mean correction
• N::NTuple{3,Int}: grid dimensions

frequency_response(m::Morlet1D{T}) -> Vector{T}

Compute the frequency response Ψ(ω) of a 1D Morlet wavelet.

Returns a length-N vector with the Fourier-domain filter coefficients. The response is analytic (zero for negative frequencies) for proper wavelet transform. Element type matches the wavelet's precision.

frequency_response(m::Morlet2D{T}) -> Matrix{T}

Compute the 2D frequency response Ψ(kx, ky) of an oriented Morlet wavelet. Element type matches the wavelet's precision.

frequency_response(m::Morlet3D{T}) -> Array{T,3}

3D frequency response Ψ(kx,ky,kz) of an oriented Morlet wavelet.

frequency_response!(Ψ, m) -> Ψ

Write the response into Ψ instead of allocating it. A bank that evaluates its filters on demand rather than storing them needs exactly one array, and this is what lets it reuse that array.

fibonacci_directions(n, ::Type{T}=Float64) -> Vector{NTuple{3,T}}

n near-uniform unit directions on the sphere (Fibonacci spiral), used as 3D wavelet orientations.

ScatteringTree{IV,RV}

Flat CSR description of the scattering tree. Pure integer topology (indices into the filter bank + integer order labels) — outside the differentiable data path, so the values are integers, but the containers stay parametric (default builders produce Vector{Int} / Vector{UnitRange{Int}}, yet MArray/CuVector/Int32 storage is permitted).

Fields

• path_data: concatenated wavelet-index lists for all paths
• path_ptr: CSR offsets (length npaths+1); path p is path_data[path_ptr[p]:path_ptr[p+1]-1]
• order: scattering order of each path (0, 1, 2, …)
• by_order: by_order[o+1] is the contiguous range of path ids of order o

build_tree(j_eff, max_order) -> ScatteringTree

Enumerate all admissible paths up to max_order from the per-wavelet effective log-scales j_eff (typically [m.j_eff for m in filter_bank.meta]). Admissibility: j_eff strictly increasing along the path. Paths are laid out grouped by order (order 0, then 1, then 2, …), so by_order ranges are contiguous.

order2_groups(tree, nw) -> Vector{Tuple{Int,Vector{Int},Vector{Int}}}

Every first-order wavelet j1 paired with its admissible order-2 children j2 and the tree path ids of those (j1, j2) pairs, ordered longest-first. The path ids are what the localized-field output is indexed by; the coefficient cascade only needs (j1, j2).

The cascade walks this instead of the flat order-2 path list so that a first-order modulus is computed and transformed once per j1 and reused across its children, rather than recomputed per path. Longest-first ordering is for the parallel backends: the child count varies by roughly 3× across scales, so equal-sized chunks of a natural-order list load-imbalance badly.

order_range(tree, o) -> UnitRange

Contiguous range of path ids with scattering order o.

path_indices(tree, p) -> view

Wavelet indices of path p (empty for the order-0 root), as a non-allocating view.

npaths(tree) -> Int

Total number of paths (including the order-0 root).

AbstractScatteringPlan

Supertype for spectral transform plans. A plan implements forward_transform! and inverse_transform!; the in-core default is DirectSumPlan, with fast paths (FFTW, AbstractFFTs/device, FINUFFT, NonuniformFFTs) provided by extensions.

DirectSumPlan{T,V,D,S}

Direct-summation DFT plan: evaluates X_k = Σ_n x_n e^{-2πi kn/N} (and its inverse) by direct summation, separably over the leading D dimensions, leaving a trailing batch axis of length nbatch untouched. Memory is O(N) — a per-axis table of roots of unity and its conjugate, never an N×N matrix. scratch is nothing for D == 1 and a full-size complex array otherwise. Containers stay parametric.

forward_transform!(out, plan, x) -> out

In-place forward (fft-convention) spectral transform. Methods provided by concrete plans.

inverse_transform!(out, plan, x) -> out

In-place inverse (ifft-convention, 1/N-scaled) spectral transform.

forward_transform(plan, x) -> X̂
inverse_transform(plan, x) -> x

Non-mutating, allocating, element-type-generic spectral transforms — the autodiff-friendly counterparts of the in-place forward_transform!/inverse_transform!. They never touch the plan's preallocated scratch and never mutate their inputs, so they accept ForwardDiff.Dual/Float32 inputs and are differentiable by reverse-mode backends. Used by the non-mutating scattering(st, x) path; the in-place ! versions remain the production hot path.

DirectSumPlan implements these as dense per-axis DFT matrix-multiplies (W*x) — differentiable by every AD backend with no special rules. The matrices are built once on first use and cached on the plan, so the O(N²) build is not repeated per call.

inverse_transform(plan, x) -> x

See forward_transform.

make_plan(spectral, T, dims; nbatch=1, kwargs...) -> AbstractScatteringPlan

Build the spectral plan selected by spectral for arrays whose leading dimensions are dims and whose element type is Complex{T}, with a trailing batch axis of length nbatch.

SpectralBackends.DirectSumSpectralBackend is the dependency-free in-core default; SpectralBackends.FFTSpectralBackend requires using FFTW; SpectralBackends.AutoSpectralBackend takes the FFTW fast path when its extension is loaded and otherwise the in-core direct sum.

make_scattered_plan(spectral, x, y, ms, T; period, solve, maxiter, rtol, eps,
                    nufft_nthreads) -> AbstractScatteringPlan

Build the scattered/nonuniform planar plan selected by spectral over points (x, y) and a uniform mode grid of size ms. SpectralBackends.DirectSumSpectralBackend is the dependency-free exact NUDFT; FINUFFTBackend and NonuniformFFTsBackend select a specific fast library; SpectralBackends.NUFFTSpectralBackend takes whichever fast library is loaded, and SpectralBackends.AutoSpectralBackend falls back to the exact direct sum when neither is.

nufft_nthreads sets the fast library's own thread count (0, the default, leaves it to the library). The direct sum accepts it and ignores it, as it does eps, so a caller can pass one set of options without first knowing which backend it will get.

spectral_backend(plan) -> SpectralBackends.AbstractSpectralBackend

The tag that would rebuild plan. Each plan type declares its own, so a transform can be reconstructed faithfully on a remote worker; plans that cannot be rebuilt from a tag alone (a device-resident FFT plan needs its device too) throw rather than report a host plan.

task_local_plan(plan) -> plan

A plan equivalent to plan that is safe to use concurrently with it. Stateless plans (FFTW, AbstractFFTs) return themselves; plans carrying mutable scratch return a copy that shares their read-only tables and owns fresh scratch. Called once per task by the parallel backends.

Three obligations follow from where this is called. A method that builds rather than shares must do so under PLANNER_LOCK, because it runs inside a spawned task by construction; it must take its thread count from per_task_nthreads, since the caller has already claimed the cores; and whoever built it must close_plan! it when the task ends.

batch_width(plan) -> Int

Number of co-located fields this plan transforms per execution — 1 unless it was built for a batch.

A nonuniform plan's batch width is fixed when its guru plan is built and cannot vary per execution, so a plan is either single-field or batched. Callers use this to decide whether a stack can go through one execution per cascade step instead of one per field.

plan_points(plan) -> (x, y)

The scattered sample locations a nonuniform plan was built on, already mapped onto its 2π-periodic domain. This is what lets a transform be rebuilt on another process, where the plan itself cannot travel. Rebuilding from these requires period = (2π, 2π), since they are already scaled.

plan_analysis(plan) -> (; solve, maxiter, rtol, damp, eps, nufft_nthreads)

How a scattered-planar plan turns samples into modes, in the form its constructor takes.

Carried rather than re-derived, because none of it follows from the points: a plan rebuilt on another process without solve analyses with the plain Type-1 adjoint where the original ran a least-squares inversion, and on irregular sampling those are different transforms, not one slightly less accurate than the other. eps is nothing for the exact direct sum, which has no tolerance to honour.

nufft_guru_make(points, type, ms, iflag, ntrans, eps, T; nthreads = 0) -> guru plan
nufft_guru_setpts!(guru, x, y) -> guru
nufft_guru_exec!(guru, input, output) -> output

Creation, point assignment and execution of a FINUFFT guru plan, split so that a device binding is one method rather than a second copy of the scattered-planar plan.

Creation dispatches on the point array, since that is what decides where the transform has to run; the other two dispatch on the returned plan, so each backend's handle carries its own execution. The host methods live in the FINUFFT extension, the CUDA ones in the cuFINUFFT extension — only the NUFFT is vendor-specific, because the cascade around it is broadcasts and reductions over whatever array type the points are.

nthreads is the library's own thread count, baked into the plan, with 0 meaning "the library's default" (all cores). Measured on the scattered-planar shapes: the library's threading is worth 2.2–3.3× at M ≳ 2·10⁴ and breaks even at M ~ 500, so the default keeps it. A plan built per task defaults to one instead — see per_task_nthreads. A device binding has no CPU threads to set and ignores it.

with_fft_nthreads(f, n) -> f()

Run f with the FFT library's global thread count set to n, restoring it afterwards.

FFTW's thread count is process-global and is raised as a side effect of loading unrelated packages, so a plan built without pinning it inherits whatever was last set — and a plan built for more threads than it needs spawns (and allocates) a task per thread on every execution. Plan builders wrap construction in this so a plan's threading is a property of the plan, not of load order.

Because the count is one global, callers hold PLANNER_LOCK across this: two builds pinning it concurrently would each restore the other's value.

The default is a no-op: only the FFTW extension has a global count to set. It takes args... so the extension's fixed-arity method is strictly more specific and adds a method rather than overwriting this one — overwriting is an error during precompilation.

per_task_nthreads(requested) -> Int

Thread count for a library plan built inside a task: requested if the caller asked for one, else 1.

An explicit request is honoured here and not just at construction, or nufft_nthreads would be a keyword that silently does nothing under a threaded backend. The default is one because task_local_plan builds a plan per task, so the tasks have already claimed the cores.

Deriving the default from Sys.CPU_THREADS instead gets this wrong twice: those are logical cores, so on a hyperthreaded machine running one Julia task per physical core — already a full machine — the arithmetic hands every task a second library thread.

close_plan!(plan) -> nothing

Release the foreign-library resources plan owns, now rather than at collection. No-op by default, and safe to call more than once.

Whoever builds a plan per task must call this when the task is done. These destructors take a lock — FINUFFT installs one to serialise its FFTW planner calls — and a lock cannot be taken from a GC finalizer, so leaving them to be collected aborts the process the moment one fires inside another task's transform. Closing eagerly leaves the finalizer nothing to do, since it checks whether the C plan is already gone.

PLANNER_LOCK

Serialises plan construction across every backend in the package.

FFTW documents fftw_execute as its only thread-safe entry point, so two plan builds running at once fault inside the planner. One lock covers them all rather than one per backend, because the planner is shared far more widely than any single backend: FFTW.jl, FINUFFT, NonuniformFFTs and FastTransforms all plan through the same libfftw3, so a lock private to one of them excludes nothing. Concurrent construction is routine — task_local_plan and SphericalCore.batch_plan build inside spawned tasks — and a build racing a build of a different library was observed as a segfault in fftw_mkapiplan.

Only construction takes this lock; transforms are never serialised, so a plan per task still executes in parallel. FastTransforms additionally needs its OpenMP thread count pinned across construction and execution — see SphericalCore.with_serial_ft.

wavelet_convolve(signal_fft, filter_fft, plan)

Perform wavelet convolution via frequency-domain multiplication then inverse transform. Allocates output. For zero-allocation hot paths, use wavelet_convolve!.

wavelet_convolve!(out, signal_fft, filter_fft, plan, buffer)

Truly zero-allocation wavelet convolution.

Multiplies signal_fft .* filter_fft into buffer in-place, then applies the inverse spectral transform via Plans.inverse_transform!(out, plan, buffer) — writing directly into out.

out and buffer must both be pre-allocated complex arrays of the same size.

apply_modulus(signal)

Apply complex modulus |·| to get envelope. Allocates output. For zero-allocation hot paths, use apply_modulus!.

apply_modulus!(out, signal)

In-place modulus. Stores |signal| in pre-allocated out. Zero allocation.

modulus_mean(signal) -> Real

⟨|signal|⟩ in a single reduction. A scattering coefficient is the mean of a modulus, so the modulus field itself is never needed unless a coarser scale consumes it — this is the leaf case, which writes nothing.

modulus_mean!(out, signal) -> Real

Write |signal| into out and return ⟨|signal|⟩. Used where the modulus field is consumed downstream; the generic method is two device-friendly passes, the CPU method fuses them into one.

spatial_average(signal::AbstractArray{T}) -> T

Compute spatial average (global mean) for translation invariance. Type-stable: returns element type T.

task_workspace(st) -> st′

A transform equivalent to st that shares its read-only parts — filter bank, path tree, work list — but owns fresh buffers and a task-local spectral plan, so the two can run concurrently.

This is what lets a parallel backend give each task private scratch without duplicating the filter bank, which dominates a transform's memory (for a 256×256 J=4 L=8 transform, 16.5 MiB of the 21.5 MiB). Methods are defined per transform type.

lsmr_solve!(x, applyA!, applyAt!, b, u, t, v, w, h, hbar;
            damp, atol, btol, conlim, maxiter) -> (; istop, iters, normr, normar, normA, condA)

Minimise ‖A x − b‖² + damp²‖x‖² in place, where applyA!(dst_pts, src_modes) applies A and applyAt!(dst_modes, src_pts) applies A†.

b is read only. u/t are point-space scratch and v/w/h/hbar mode-space scratch; the two destinations t and w exist because the transforms overwrite rather than accumulate. Every array operation is copyto!, fill!, norm or a fused broadcast, so the buffers may live on a device.

istop says why it stopped: 1 the residual met btol, 2 the least-squares optimality met atol, 3 the condition estimate hit conlim, 5/6 those quantities reached the precision floor, 7 maxiter. Stopping at 7 is not a failure — both ‖r‖ and ‖A†r‖ decrease monotonically, so the iterate is the best one seen.

lsmr_solve_batched!(x, applyA!, applyAt!, b, u, t, v, w, h, hbar, work;
                    damp, atol, btol, conlim, maxiter) -> (; istop, iters, normr, normar, …)

lsmr_solve! over a stack of B right-hand sides sharing one operator.

A batched NUFFT plan's width is fixed when it is built, so no column can be transformed on its own and the stack must advance together — which turns every scalar in the recurrence into one per column. They advance on the host through the same lsmr_step the single-column path uses, and return to the device as three coefficient arrays.

A column that has stopped is frozen: its coefficients go to zero, so it contributes nothing further. It is not compacted out of the stack, because the transform width cannot shrink — removing it would save no work while costing the bookkeeping that makes a permuted, partially-retired stack correct.

istop/normr/normar describe the worst column, so a caller that checks them sees the whole stack.

lsmr_init(alpha, beta, ::Type{T}) -> LSMRState{T}

Initial state from the first bidiagonalisation pair, beta = ‖b‖ and alpha = ‖A†b/beta‖.

lsmr_step(state, alpha, beta, damp, atol, btol, ctol) -> (state, chbar, cx, ch)

Advance the scalar recurrence one iteration and return the three coefficients the vector updates need:

hbar .= h .+ chbar .* hbar
x    .= x .+ cx    .* hbar
h    .= v .+ ch    .* h

alpha/beta are this iteration's bidiagonalisation norms. damp is the Tikhonov λ, which enters only through one extra rotation, so λ = 0 costs a rotation of a zero and nothing else.

LSMRState{T}

The scalar state of an LSMR iteration: the bidiagonalisation and rotation quantities, the running estimates of ‖r‖, ‖A†r‖, ‖A‖ and cond(A), and the stopping code.

normA accumulates Σ(α² + β²), so it estimates the Frobenius norm of the bidiagonalisation built so far — bounded above by ‖A‖_F and generally above ‖A‖₂. That is what the optimality test ‖A†r‖ ≤ atol·normA·‖r‖ is scaled by, which makes it slightly conservative rather than wrong.

Immutable and advanced by lsmr_step, which touches no arrays — so the batched solver can hold one of these per column and drive them from the host while the vector work stays on the device.

BatchedLSMRWork{A2,A3,HV,SV}

Per-column bookkeeping for lsmr_solve_batched!: the two reduction targets, the three coefficient arrays the vector updates broadcast against, host mirrors of each of those five, and one LSMRState per column.

The device arrays are (1, B) over points and (1, 1, B) over modes so they broadcast against the stacks with no reshape in the loop, and they are preallocated because sum(abs2, x; dims = …) allocates on every call. Each norm is held twice — as both ranks, over the same memory — because the point stack is rank 2 and the mode stack rank 3, and a (1, 1, B) array broadcast against (M, B) would expand to (M, B, B) rather than scaling columns. Each quantity nonetheless gets its own buffer: the recurrence needs this iteration's alpha while the previous one is still live.

default_solver_rtol(::Type{T}, spectral, eps) -> T

Stopping tolerance for the scattered least-squares solve, chosen from the accuracy of the transform the solve is built on rather than from a fixed constant.

Dispatched on the spectral backend, exactly as SphericalCore.default_rtol is: the in-core direct summation evaluates an exact conjugate pair, so the arithmetic is the only limit and sqrt(Base.eps(T)) is the conventional least-squares target. A fast NUFFT's Type-1 and Type-2 are adjoints of each other only to about eps, so the operator the solver sees is Hermitian only to that accuracy and a target below it asks for something the transform cannot express — which is how a Float32 plan came to be asked for 1e-8, two orders under its own 1e-6.

Base.eps is spelled out because eps is a keyword-argument name in every plan constructor here.

default_nufft_eps(::Type{T}) -> Float64

The NUFFT tolerance a fast backend uses when the caller names none: 1e-6 for Float32, 1e-9 otherwise. One definition, because both fast backends must ask their library for the same accuracy.

warn_underdetermined(M, ms, solve, damp) -> nothing

Warn, once per session, that a solve has more modes than samples and no damping to resolve them.

With prod(ms) > M at least prod(ms) - M mode directions are constrained by no sample, so the least-squares problem has no unique solution and the answer depends on the regularisation. LSMR returns the minimum-norm iterate and its early termination is itself a regulariser, which is why this warns rather than throws — but the caller is the only one who can say whether that is the answer they wanted, or whether they should coarsen ms, add samples, or set damp for explicit Tikhonov. Both numbers are named because the ratio is what decides how ill-posed the solve is.

No λ is chosen for the caller. Sweeping λ against a dense Tikhonov reference on this operator shows why: on quasi-uniform points A is a scaled isometry and every λ across six decades returns the same coefficients to 1e-14, while on clustered points a λ small enough to be safe changes neither the residual nor ‖x‖ in the first five digits. A default that cannot be observed to do anything is worse than none, since it would still be a number the caller has to reason about.

fftw_plan(T, dims; nbatch, planning, fft_nthreads) -> AbstractScatteringPlan

Build the FFTW-backed plan. The real method lives in the FFTW extension; the definition here is a throwing stub, so an explicit FFTSpectralBackend() request dispatches straight to the extension when it is loaded and to an actionable error when it is not — with no capability lookup on either path.

abstractffts_plan(dummy_device_array; region=1:ndims(dummy_device_array))

Vendor-neutral device FFT plan constructor. Declaration only — the sole method is provided by the KernelAbstractions extension, which builds forward/inverse plans via AbstractFFTs.plan_fft / plan_ifft on dummy_device_array. Because those dispatch on the array type, the same builder yields cuFFT (CuArray), rocFFT (ROCArray), or FFTW (plain Array) plans — so the GPU scattering path is device-agnostic. region selects the transformed dimensions (e.g. (1, 2) for a batched (Ny, Nx, B) stack, leaving the batch axis untouched).

finufft_scattered_plan(x, y, ms, T; period, solve, maxiter, rtol, eps, nufft_nthreads)
nonuniformffts_scattered_plan(x, y, ms, T; period, solve, maxiter, rtol, eps, nufft_nthreads)

Fast-path scattered-planar plan constructors. The real methods live in the FINUFFT and NonuniformFFTs extensions; the definitions here are throwing stubs, so naming one of those backends explicitly costs a plain dispatch rather than a capability lookup.

See finufft_scattered_plan.

FINUFFT fast path for scattered/nonuniform planar points; requires using FINUFFT.

NonuniformFFTs.jl fast path for scattered/nonuniform planar points; requires using NonuniformFFTs.

resolve_backend(backend) -> AbstractLocalBackend

Concrete backends pass through unchanged — a request is honoured exactly or refused, never silently downgraded. AutoBackend resolves on real capability: ThreadedBackend when there is more than one thread and the OhMyThreads extension is loaded, otherwise SerialBackend.

check_available(backend) -> backend

Throw unless backend can actually execute, naming the package that would enable it. Called by the batch entry points so an unloadable request fails immediately instead of at a MethodError.

true when the OhMyThreads extension is loaded, i.e. ThreadedBackend can actually run.

true when the KernelAbstractions extension is loaded, i.e. GPUBackend can actually run.

true when the Distributed extension is loaded.

true when the MPI extension is loaded.

cascade!(S1, S2, st, signal_fft) -> (S1, S2)

Both scattering orders in one pass over the tree, grouped by first-order wavelet:

for (j1, children):  U₁ = |x ⋆ ψ_j1| ;  S1[j1] = ⟨U₁⟩
                     Û₁ = fft(U₁)     ;  S2[j1,j2] = ⟨|U₁ ⋆ ψ_j2|⟩  for each child

The first-order convolution is therefore evaluated once, not once for S1 and again for S2, and only one U₁/Û₁ pair is live at a time rather than one per wavelet. Wavelets with no admissible child skip the modulus buffer entirely, reducing to a single fused ⟨|·|⟩.

signal_fft is read only, so the caller's preserved signal spectrum survives the call.

cascade!(S1, S2, st, image_fft) -> (S1, S2)

Both scattering orders in one grouped pass — see the 1D cascade! for the scheme. Admissibility is j_eff strictly increasing, i.e. scale strictly increasing across all orientation pairs, which is what the tree encodes; same-scale different-orientation pairs are not order-2 paths.

cascade!(S1, S2, st, vol_fft) -> (S1, S2)

Both scattering orders in one grouped pass — see the 1D cascade! for the scheme.

(st::ScatteredPlanarScattering)(x) -> ScatteringCoefficients2D

Apply the scattered planar scattering transform to a length-M vector of samples at the plan's points. Allocates a coefficient container per call; use scattered_planar_scattering! to reuse one.

scattered_planar_scattering!(coeffs, st, x) -> coeffs

In-place counterpart of the callable: write into a preallocated ScatteringCoefficients2D. Allocation-free when coeffs carries a 1-element S0 (so update_S0 writes rather than rewraps).

scattered_planar_scattering_batch!(S0, S1, S2, st, X) -> (; S0, S1, S2)

Cascade over a stack of B fields sampled at the same points: X is (M, B), S0 is (B), S1 is (nw, B) and S2 is (nw, nw, B).

Every transform covers the whole stack in one execution, so the cascade issues its 1 + nparents + nw + npaths transforms once rather than once per field. st must have been built with a matching ntrans, since a guru plan's batch width is fixed when it is made.

AbstractSphericalPlan

Supertype for spherical spectral plans (scattered NUFSHT, structured SHT, …). A concrete plan must implement sphere_coeffs, sphere_apply!, and sphere_mean.

sphere_coeffs(plan, field) -> C

Spherical-harmonic analysis: the coefficients C of field (up to the plan's band limit). For a structured grid this is the fast forward SHT; for scattered points it is the exact (least-squares / CG) inversion — not the adjoint, which mis-scales the coefficients. Returned opaquely and consumed only by sphere_apply!/sphere_mean on the same plan.

sphere_coeffs!(C, plan, field) -> C

In-place sphere_coeffs: analyse field into the pre-allocated coefficient container C, which must have come from sphere_coeffs_buffer(plan). This is what lets the cascade re-analyse a first-order field without allocating a coefficient vector per scale.

sphere_coeffs_buffer(plan) -> C

A coefficient container of the right type and size for plan, suitable for sphere_coeffs!. Backends whose coefficients are not a plain vector (the FastSphericalHarmonics triangular layout, the NUFSHT dense spin layout) return their own shape.

sphere_field_buffer(plan) -> field

A grid container of the right type and size for plan, for a band synthesised at that plan's band limit. Defaults to sphere_coeffs_buffer, which is the same shape on backends whose coefficients and samples share a layout.

sphere_apply!(out, plan, C, h) -> out

Apply the per-degree multiplier h(ℓ) to (a copy of) the coefficients C and synthesise the result into out — i.e. out = Σ_{ℓm} h(ℓ) · C_{ℓm} · Y_{ℓm}. Does not mutate C.

sphere_mean(plan, field) -> scalar

Spherical average of field under the plan's sampling: an unweighted sample mean for (quasi-uniform) scattered points, the exact quadrature integral for a structured grid. Provided by each backend.

sphere_plan_at(plan, lmax) -> plan or nothing

A plan of the same kind at the lower band limit lmax, or nothing when this backend cannot be narrowed. nothing is not a defect — it means every band is synthesised at the full band limit, which is what the cascade did before.

This is the spherical form of decimation. Band-pass wavelet j is confined to degrees ≲ ℓ_j, so synthesising it at the full lmax transforms a grid whose resolution the band cannot use.

sphere_restrict!(Cdst, pdst, Csrc, psrc) -> Cdst

Copy the spherical-harmonic coefficients of degree ≤ pdst's band limit out of Csrc (laid out for psrc) into Cdst (laid out for pdst), zeroing anything Cdst holds beyond them.

Coefficient layout is backend-specific, so narrowing a band limit is a backend operation rather than an array slice.

band_lmax(lmax, J, j, headroom) -> Int

Band limit that carries wavelet j: headroom · ℓ_j, clamped to lmax.

Band j is a difference of Gaussians whose transfer above its cutoff is exp(-ℓ(ℓ+1) / 2ℓ_j(ℓ_j+1)), which decays slowly — 0.135 at 2ℓ_j, 0.011 at 3ℓ_j, 3.4e-4 at 4ℓ_j. headroom is therefore what is actually traded here, and 4 keeps the discarded tail below the transform's other approximations.

band_plans(plan, lmax, J, headroom) -> Vector or nothing

One narrowed plan per scale, or nothing if this backend cannot narrow. A scale whose band needs the full band limit reuses plan itself rather than building a duplicate.

SphericalScattering{T,P,V}

Spherical scattering transform over a backend plan::P. sigma2[k+1] is the Gaussian-transfer variance for the dyadic low-pass cutoff ℓ_k (k=0..J); band-pass wavelet j is lowpass(ℓ_j) − lowpass(ℓ_{j-1}).

SphericalMonogenicScattering{T,P}

Spherical monogenic scattering: shares the dyadic difference-of-Gaussians bands of SphericalScattering but replaces the analytic modulus with the spherical monogenic amplitude A_j = √(U⁰_j² + |∇_S g_j|²) (spin-0 Bochner identity — see monogenic_amplitude!).

SphericalWorkspace{A,C}

Scratch for one spherical cascade: the band-pass output band, the current first-order field u1, and two coefficient containers — C for the input field, C1 for the first-order field being re-analysed. Two is all the cascade ever needs, because it finishes every child of a scale before starting the next.

Build one with SphericalWorkspace(st, field) and reuse it across calls; a task that transforms concurrently needs its own (and its own Plans.task_local_plan of the spherical plan).

SphericalMonogenicWorkspace{A,C,S}

Scratch for one spherical monogenic cascade: the current amplitude field u1, a second amp for the order-2 amplitudes, two coefficient containers, and scratch — the (g, lapg, g2, lapg2) fields the Bochner identity needs (see monogenic_amplitude!).

spherical_scattering!(S1, S2, st, ws, field) -> (; S0, S1, S2)

In-place spherical scattering into pre-allocated S1/S2 using the workspace ws (see SphericalWorkspace). Allocation-free once ws exists.

Grouped by first-order scale so only one first-order field and one coefficient vector are live at a time: scale j1 is band-passed, averaged into S1[j1], analysed once, and then consumed by every coarser j2 < j1 before the next j1 begins.

spherical_scattering_batch!(S1, S2, st, ws, X) -> (; S0, S1, S2)

Cascade over a stack of B fields sampled at the same points: X is (M, B), S1 is (J, B) and S2 is (J, J, B). Every transform covers all B fields in one call, so the per-call setup is paid once per cascade step rather than B times — st.plan must have been built with a matching batch size for that to hold.

spherical_monogenic_scattering!(S1, S2, st, ws, field) -> (; S0, S1, S2)

In-place spherical monogenic scattering — the monogenic counterpart of spherical_scattering!, grouped by first-order scale so one amplitude field is live at a time rather than all J.

monogenic_amplitude!(amp, st, C, j, w) -> amp

Spherical monogenic amplitude at scale j of the field whose (already-computed) SH coefficients are C, written into amp. w is a NamedTuple of scratch fields (g, lapg, g2, lapg2).

Uses only spin-0 transforms via the Bochner/product identity: with g_j = (−Δ_S)^{-1/2} U⁰_j,

|U^R_j|² = |∇_S g_j|² = ½ Δ_S(g_j²) − g_j · Δ_S g_j,

so A_j = √(U⁰_j² + |∇_S g_j|²). The Riesz vector itself (needed for orientation/phase) requires spin-1 synthesis and is handled separately by spherical_monogenic_components.

task_local(st) -> st

A copy of the spherical transform safe to run concurrently with the original: the design matrix, Gram factor and point set are read-only and shared, while the plan's analysis/synthesis scratch is duplicated by Plans.task_local_plan. Analysis writes through that scratch, so tasks sharing one plan would overwrite each other's coefficients.

Per-degree multiplier for the difference-of-Gaussians band-pass wavelet j.

dog_sigma2(lmax, J, T) -> Vector{T}

Gaussian-transfer variances for the J+1 dyadic low-pass cutoffs ℓ_k = lmax / 2^(J-k), k=0..J (ℓ_J = lmax). Band-pass wavelet j is lowpass(ℓ_j) − lowpass(ℓ_{j-1}).

structured_grid(lmax, T) -> (Θ, Φ)

Colatitudes Θ (length lmax+1) and longitudes Φ (length 2lmax+1) of the equiangular structured grid. Matches FastSphericalHarmonics.sph_points(lmax+1).

plan_points(plan) -> (θ, φ) or nothing

The sample locations a spherical plan was built on, or nothing when the plan is defined by its band limit alone (a structured grid). This is what a transform needs to be rebuilt on another process, where the plan itself cannot travel.

plan_weights(plan) -> weights or nothing

The quadrature weights the plan averages with, or nothing for a backend that takes the unweighted sample mean.

These have to be carried, not re-derived: a structured grid's Fejér weights vary by a factor of ~5 across colatitudes, so rebuilding a plan on the same points without them silently substitutes a uniform average and shifts every coefficient.

plan_solver(plan) -> (; rtol, maxiter)

The iterative-solver settings a plan analyses with, so a rebuilt plan inverts to the same tolerance rather than to whatever the constructor defaults to.

batch_plan(plan, B) -> plan or nothing

An equivalent plan that transforms B co-located fields per call, or nothing if this backend has no batched form. Same points, band limit, weights and solver settings — only the batch width differs.

A backend that returns a plan here lets spherical_scattering_batch! issue one transform per cascade step for the whole stack instead of one per field. Returning nothing is not a defect; it means the per-field loop is the only path, and callers fall back to it.

May return plan itself when it is already B wide, so a per-task caller wants task_local_batch_plan.

task_local_batch_plan(plan, B) -> plan or nothing

A B-wide plan safe to apply while plan is applied in another task.

A spherical plan carries the scratch its analysis solves through, and batch_plan may hand back plan itself at its own width — shared across tasks that corrupts the solve into a diverging residual instead of erroring. Widening to another width already builds a fresh plan, so only the same-width case copies.

supports_batch(plan) -> Bool

Whether batch_plan can widen this plan. Asked before building anything, so a caller can choose between the batched and per-field cascades without paying for a plan it may not use.

plan_nufft(plan) -> SpectralBackends.AbstractSpectralBackend

The NUFFT backend a scattered-sphere plan resolved to, or AutoSpectralBackend() for a plan that runs none.

Plans.spectral_backend cannot express this: it reports NUFSHTSpectralBackend whether the transform is driven by FINUFFT, NonuniformFFTs, or direct summation. Carrying it separately is what lets a rebuilt transform — on a distributed worker, say — run the same transform as the original instead of silently re-resolving to whatever that process happens to have loaded.

plan_spin(plan) -> (spin-0 plan, spin-1 plan) or nothing

The spin-weighted plan pair a scattered-sphere plan carries for spherical_monogenic_components, or nothing for a plan built without one.

Built with the plan rather than per call because a spin plan is the expensive part of that decomposition — two NUFFTs plus the Wigner-d recurrence tables — while the transform it drives is one synthesis. Only the monogenic constructors ask for it, so a plain spherical transform pays nothing. Like every other buffer a plan owns, the pair is working state: concurrent use needs Plans.task_local_plan.

AnalysisNotConverged <: Exception

Thrown when an iterative analysis cannot produce a usable answer.

Not thrown merely for stopping early: LSMR decreases both ‖r‖ and ‖A†r‖ monotonically, so an iterate truncated at maxiter is the best one seen and is returned. This is for the cases where the answer means nothing — a non-finite residual, or a conditioning estimate past what the precision can represent — because those propagate into every coefficient downstream, and a silent NaN is worse than a stop.

default_rtol(spectral) -> Real

Conjugate-gradient tolerance for the scattered-sphere analysis solve under spectral.

The solve cannot resolve the field better than the transform it is built on, so the useful tolerance is backend-specific. In-core direct summation evaluates Y exactly, so the solve tolerance is the analysis accuracy and 1e-8 buys real digits. NUFSHT analysis instead floors at its own approximation error — measured against the exact in-core transform, ~2e-4 at lmax = 8, M = 500 and ~2e-5 at lmax = 16, M = 1500 — and 1e-4 reaches that floor, matching 1e-8's accuracy in up to 3× less time. 1e-2 does not: it lands about twice as far off.

with_serial_ft(f, ft_nthreads, ft_nthreads!)

Run f() with FastTransforms pinned to a single thread, reading the current thread count with ft_nthreads() and setting it with ft_nthreads!(n) (each backend passes its own accessors, since the symbols live in its own dependency).

The pin is a correctness requirement rather than a tuning choice. FastTransforms runs its butterfly and sphere transforms inside OpenMP parallel regions, and entering one from a non-root Julia task silently corrupts the result — and, during plan construction, segfaults. Pinned to one thread it takes its serial code path instead.

Reference counted, because the thread count is a process global: the first section entered records the previous value and sets one, and only the last section out restores it. Setting and restoring per section would let one section's restore re-enable OpenMP underneath another section still running in a different task — exactly the corruption the pin exists to prevent — and would leak a "previous" value of 1 whenever two sections overlapped.

make_spherical_plan(spectral, θ, φ, lmax, T; rtol, maxiter, spin) -> AbstractSphericalPlan

Build the scattered-sphere plan selected by spectral over points (θ, φ) and band limit lmax. SpectralBackends.DirectSumSpectralBackend is the dependency-free default; SpectralBackends.NUFSHTSpectralBackend (or SpectralBackends.AutoSpectralBackend once the NUFSHT extension is loaded) uses the NUFSHT fast path.

The in-core plan transforms one field per call, so ntrans is accepted and ignored rather than rejected: a caller asking for a batch still gets correct results, just not the batched transform. It ignores spin for a stronger reason — its monogenic decomposition needs no spin-weighted synthesis at all, taking the surface gradient of (−Δ_S)^{-1/2}U⁰ instead (see plan_spin).

make_structured_plan(spectral, lmax, T; rtol, maxiter) -> AbstractSphericalPlan

Build the structured-sphere plan selected by spectral. SpectralBackends.DirectSumSpectralBackend is the dependency-free default (direct SHT on the grid); SpectralBackends.FSHTSpectralBackend (or SpectralBackends.AutoSpectralBackend once the FastSphericalHarmonics extension is loaded) uses the fast exact SHT.

nusht_spherical_plan(θ, φ, lmax, T; rtol, maxiter)

Fast-path scattered-sphere plan constructor. The real method lives in the NUFSHT extension; this is its throwing stub.

fsh_structured_plan(lmax, T)

Fast-path structured-sphere plan constructor. The real method lives in the FastSphericalHarmonics extension; this is its throwing stub.

plot_coefficients(c; …) — plot scattering coefficients. Requires `using CairoMakie`.

plot_filter_bank(fb) — plot a filter bank. Requires `using CairoMakie`.