Calculations

The workhorse of this package, split into focused files under src/Calculations/

append_helmholtz_rotational_divergent_rows(sums, counts, distance_bins)

Append Helmholtz-derived rotational/divergent rows to native six-row single-pass sums and counts. Rows 1 through 6 are copied unchanged; rows 7 and 8 are rotational and divergent second-order components.

auxiliary_joint2d!(sums, counts, sf_type, x, u, distance_bins, value_bins)

Single-type joint 2D batch (serial); output (n_dist, n_val, batch…).

auxiliary_shared_positions!(sums, counts, x_mat, u_batch, sf_type, distance_bins; strip_width=32)

Fixed geometry batch: x is (N_dims, N), u has trailing batch dims.

auxiliary_varying_positions!(sums, counts, x_batch, u_batch, sf_type, distance_bins)

Varying geometry batch: x and u have matching trailing batch dims.

calculate_structure_function!(sums, counts, sf, x, fields, distance_bins; backend, kwargs...)

Accumulate a multi-field's pairs into sums/counts on backend.

calculate_structure_function(sf, x, fields, distance_bins[, count_eltype]; kwargs...)

Structure function of a multi-field: several quantities sampled at the same points, swept together in one pass, so a mixed moment reads every field at each pair.

calculate_structure_function(sf, x, u, nodes::HarmonicNodes, spectral_backend; distance_metric, weights, valid, output_type)

The kernel-binned structure function of u sampled at the points x of a sphere, at the nodes' separations, by spherical harmonic pseudo-coefficients (see harmonic_sweep!). x is (lon, lat) in the angle unit of distance_metric, which must be spherical; u is a (D, cells...) field of (east, north[, radial]) components, a (1, cells...) scalar, or a multi-field, with prod(cells) == size(x, 2). weights default to one per point; on a grid the cell measure makes the statistic an area average.

The result's distance is the HarmonicNodes object itself, one value per node, and its counts are the kernel-weighted pair counts, floating point.

calculate_structure_function(sf, schedule::ScatteredModesSchedule, u, distance_bins, spectral_backend; weights, backend, output_type, verbose)

The soft-binned structure function of scattered points by non-uniform FFT. u is (D, N) over the N points of schedule, or a multi-field over them; spectral_backend is a provider tag, NonuniformFFTsSpectralBackend or FINUFFTSpectralBackend. Counts are Float64, the kernel-weighted pair mass, and the result's distance is a ModeBinEdges carrying the schedule. See ScatteredModesSchedule for what is and is not exact.

calculate_structure_function_2d_batch!(sums, counts, sf_type, x, u, distance_bins, value_bins; backend=..., ...)

Batch 2D joint histograms over (N_dims, N_points, T); outputs (n_dist, n_val, T).

calculate_structure_function_batch!(sums, counts, sf_type, x, u, distance_bins; backend=..., workspace=nothing, ...)

Batch structure functions over the third dimension of matrix inputs (N_dims, N_points, T). Host sums, counts must have shape (NB, T) where NB = length(distance_bins) - 1.

GPU: fully implemented via GPUBackend when KernelAbstractions is loaded. CPU backends: not yet implemented (use a loop over t with calculate_structure_function!).

calculate_structure_function_tensor!(sums, counts, order, x, u, distance_bins[, axis_bins]; second_axis, backend, distance_metric)

Accumulate the rank-order increment moment tensor of a point list into sums (D, …, D, n_bins[, n_axis], auxiliary...) and counts (n_bins[, n_axis], auxiliary...) on backend, the in-place form of calculate_structure_function_tensor. With axis_bins the tensor is binned jointly in separation and in the angle second_axis reads from each pair.

calculate_structure_function_tensor(order, x, u, distance_bins, axis_bins; second_axis, backend, output_type, count_eltype, distance_metric)

The increment moment tensor joint in separation and the angle second_axis reads from each pair's direction, sums of shape (D, …, D, n_bins, n_axis): the tensor resolved by direction. One field (no auxiliary axes), a flat metric, the CPU backends.

calculate_structure_functions_single_pass!(sums, counts, x, u, distance_bins; backend=CB.SerialBackend(), kwargs...)

Accumulate into pre-allocated $(6, n_bins)$ buffers using the requested execution backend.

calculate_structure_functions_single_pass(x, u, distance_bins; backend=CB.AutoBackend(),
                                          output_type=StructureFunction, kwargs...)

Compute the six native invariant structure functions (S2, L2, T2, S3, L3, L1T2) in one pair pass, returned as a NamedTuple keyed by invariant. Each entry is a single-operator result of the requested output_type (default the averaged StructureFunction; pass StructureFunctionSumsAndCounts for the raw sums+counts). For point-field input a :helmholtz entry (a HelmholtzDecomposition2D) is included.

calculate_structure_functions_single_pass_2d!(sums_3d, counts_3d, x, u, distance_bins, value_bins; backend, distance_metric, kwargs...)

Accumulate the six invariants' joint distance × value histograms of a point list into sums_3d and counts_3d, (6, n_bins, n_val) each, on backend; the in-place form of calculate_structure_functions_single_pass_2d. value_bins is one edge vector for every invariant or a tuple of six.

calculate_structure_functions_single_pass_2d(x, u, distance_bins, value_bins; backend=CB.AutoBackend(),
                                             output_type=StructureFunction2DSumsAndCounts, kwargs...)

Compute the six invariant 2D joint structure-function histograms in one pass, returned as a NamedTuple keyed by invariant (S2, L2, T2, S3, L3, L1T2). Each entry is a StructureFunction2DSumsAndCounts view into the stacked accumulator (the 2D joint histogram has no averaged form, so output_type must be StructureFunction2DSumsAndCounts).

calculate_structure_functions_single_pass_2d_batch!(sums, counts, x, u, distance_bins, value_bins; backend=..., ...)

Batch six invariant distance × value joint histograms over (N_dims, N_points, T); outputs (6, NB, n_val, T). Pass shared bin types or NTuple{6,...}; use Tuple(v...) if you have a length-6 vector of bin objects.

calculate_structure_functions_single_pass_batch!(sums, counts, x, u, distance_bins; backend=..., ...)

Batch six invariant 1D distance histograms over (N_dims, N_points, T); outputs (6, NB, T).

cell_measure(grid) -> Vector

The measure — length, area or volume — of every cell of grid, in the cell order the gridded entries use. As weights = cell_measure(grid) it makes a structure function an area average. Supplied by the FlowGeometries extension.

covariance(result, variance) -> (separations, C)

Covariance function from a second-order structure function and the field's variance.

A structure function is twice the variogram, D(r) = 2[C(0) - C(r)], so

\[C(r) = C(0) - D(r)/2\]

with C(0) the variance. variance must be supplied and cannot be recovered from D, which is blind to it — the same blindness that makes k = 0 unavailable to isotropic_spectrum.

Valid only if the field is second-order stationary. A field that is merely intrinsically stationary has a variogram but need not have a finite variance, and then no covariance exists to compute; supplying a number anyway produces a curve with no referent.

covariance_matrix(points, separations, C; metric, check_posdef = true, posdef_rtol = 1e-6)

Covariance matrix over points, evaluating the covariance function (separations, C) at each pair distance by linear interpolation and holding it constant outside the sampled range.

A covariance matrix must be positive semi-definite, and one built this way need not be, so it is checked. Two different things trip the check and the message distinguishes them: a covariance function that is not a valid kernel at all, and one that is valid but sampled too coarsely — interpolating a kernel does not preserve positive-definiteness, and the error falls as the square of the separation spacing. check_posdef = false returns the matrix regardless.

Loop-over-slice gold reference for batch parity.

enstrophy_flux(operator::VectorDotSFType, separations, SF_Au, wavenumbers)
enstrophy_flux(result, wavenumbers)

Interscale enstrophy flux of a two-dimensional isotropic flow from the velocity's advective structure function SF_Au = ⟨δu · δ𝓐_u⟩, 𝓐_u = u·∇u, sampled at separations up to R:

Π^ω_K = (K³/2) ∫₀^R SF_Au(r) [J₃(Kr) - J₁(Kr)]/2 dr + (K²/2) [ SF_Au(R) J₂(KR) + (1/K) SF_Au'(R) J₁(KR) ]

The J₁ relation on ⟨δω δ𝓐_ω⟩ = -∇² SF_Au integrated by parts twice; [J₃ - J₁]/2 = -J₂'. The slope SF_Au'(R) is the one-sided difference of the last two samples, so the boundary term is sensitive to noise in the last bins. operator must name two distinct fields, the velocity and its advection.

fit_flux(sf, k_edges, method; W = nothing) -> (k, F, ξ, ε, covariance, flux_covariance)

The piecewise-constant spectral energy flux that reproduces the third-order structure function sf = ⟨δu_L ‖δu‖²⟩ (S3SFType) of a two-dimensional isotropic flow through FluxForwardModel. F is the flux at the bin centres k (positive towards small scales), ξ the injection density on the bins and ε the flux towards large scales below the first bin; covariance is over (ε, ξ…) and flux_covariance = G C Gᵀ over F.

fit_helmholtz_spectra(L2, T2, k_edges, method; W = nothing) -> (k, E, B, covariance)

The gradient (E) and curl (B) shell spectra on the bins k_edges that reproduce the longitudinal and transverse second-order structure functions of a two-dimensional field through HelmholtzForwardModel. L2 and T2 share their bins; W covers the stacked values [D_LL; D_TT] over the bins both hold (or over all bins of each). The covariance is over the stacked [E; B].

fit_spectrum(sf, k_edges, method, ::Val{D}; W = nothing) -> (k, E, covariance, …)

The shell spectrum on the wavenumber bins k_edges that reproduces the second-order trace sf through SpectrumForwardModel, fitted by method. W is the data covariance: the variance of each bin's value (a vector over the bins holding a value, or over all bins), or a full matrix. Returns the bin centres k, the fitted E on the bins and the posterior covariance (nothing for a method that gives none). A SegmentedPowerLaw returns E evaluated at k, and in addition the fitted parameters, their covariance, the segment breakpoints and converged.

The matrix taking a flux model's (ε, ξ…) to the flux at the bin centres.

The matrix of a forward model.

gpu_calculate_structure_function(...)

GPU-accelerated structure function calculation. Requires loading KernelAbstractions.jl to activate the GPUExt extension. The backend can be KernelAbstractions.CPU() (for testing parity) or any GPU backend like CUDABackend() from CUDA.jl.

Device histogram buffers are UInt32; count_eltype (default UInt32) selects the host count type after download. Pass workspace=GPUSFWorkspace(...) to reuse device histogram buffers across repeated calls (see GPUSFWorkspace).

This stub exists so the extension can legally extend this function.

gpu_calculate_structure_function!(output_sums, output_counts, ...)

In-place GPU structure function reduction. Requires the GPUExt extension. Accumulates into caller-owned output_sums and output_counts (same contract as serial_calculate_structure_function! / threaded_calculate_structure_function!). Reuses host buffers only; pass workspace=GPUSFWorkspace(...) to reuse device histogram buffers across calls.

gpu_calculate_structure_function_2d(sf_type, backend, x_mat, u_mat, distance_bins, value_bins; kwargs...)

GPU 2D joint structure function (distance × SF value histogram) for one sf_type. Requires loading KernelAbstractions.jl to activate the GPUExt extension. Device counts are UInt32; count_eltype selects the host matrix type after download.

gpu_calculate_structure_function_2d_batch(sf, backend, x, u, distance_bins, value_bins; kwargs...)

The value-binned joint histogram of a field with auxiliary axes on a device, one histogram per auxiliary slice; supplied by the KernelAbstractions extension.

gpu_calculate_structure_functions_single_pass_2d(backend, x, u, distance_bins, value_bins; kwargs...)

Six invariant native distance × value joint histograms on a KernelAbstractions backend. Requires loading KernelAbstractions.jl to activate the GPUExt extension. Device counts are UInt32; count_eltype selects the host array type after download.

gpu_calculate_structure_functions_single_pass_2d!(sums, counts, backend, x, u, distance_bins, value_bins; kwargs...)

In-place GPU 2D single-pass accumulation. Requires the GPUExt extension.

gridded_spectrum(u, schedule, ::Val{D}, spectral_backend; valid, weights, taper, missing_lags) -> (wavenumbers, density)

Spectral density of a gridded field, by transforming its structure function over the whole lag space.

No direction is averaged over and no separation is binned, so this carries none of the angular assumption isotropic_spectrum makes. It reads the field only through its structure function and the variance of its held cells, so it accepts a field with cells missing, which is the case where the field's own transform is meaningless while the pair average is still unbiased.

On a complete periodic grid the result is the field's own spectrum to round-off. With cells missing, or on a bounded direction, it is an estimate: the exact transform of the exact masked structure function, equal to the complete field's spectrum in expectation, with a statistical error set by the pair count behind each lag. A bounded direction is padded so the transform of the lags |h| < n is their linear transform; wavenumbers then has the padded length along it. taper weights the lags (see AbstractTaper) and missing_lags says what to do with a lag no held pair names (see AbstractMissingLagPolicy). weights, one per cell, makes every pair sum and pair count a w_k w_kp-weighted one and the variance a weighted variance, so the estimate is that of the weighted structure function.

wavenumbers is one angular-wavenumber vector per grid direction; density is the Dg-dimensional array over those, on the same convention as isotropic_spectrum — integrating it over d^D k returns the variance of the held cells, so sum(density) * prod(step) does exactly, with step the wavenumber spacing of each direction. The k = 0 value is the sum of the weighted autocovariance over the lags: zero to round-off on a complete periodic grid, where no pair sees the mean, and the estimator's own low-wavenumber value otherwise.

harmonic_spectra(x, u, lmax, spectral_backend; distance_metric, weights, valid) -> (l, C) or (l, EE, BB, EB)

Pseudo-spectra of the masked, weighted field: C̃_l = (1/(2l+1)) Σ_m |ũ_lm|² for a (1, N) scalar; for a (2, N) tangent vector in (east, north), the gradient (E) and curl (B) spectra and their cross spectrum from the spin-1 coefficients ₁U_lm of U = u_θ + i u_φ and the spin-(−1) coefficients ₋₁Ū_lm of its conjugate, which follow from the first as ₋₁Ū_lm = (−1)^{m+1} conj(₁U_{l,−m}):

Φ_lm = (₁U_lm − ₋₁Ū_lm) / (2√(l(l+1))),   Ψ_lm = −i (₁U_lm + ₋₁Ū_lm) / (2√(l(l+1))),
C^E_l = l(l+1) (1/(2l+1)) Σ_m |Φ_lm|²,   C^B_l likewise from Ψ,   C^{EB}_l from Φ conj(Ψ),

so that Σ_l (2l+1)/(4π) (C^E_l + C^B_l) is the mean square of u on a complete sphere with exact quadrature weights. On a masked or unevenly sampled sphere these are the pseudo-spectra of the window and the field together, the input to the kernel-binned statistics.

harmonic_sweep!(sums, counts, sf, geometry, x, weights, data, nodes, ::Val{D}, ::Val{V}, ::Val{K}, spectral_backend; valid)

Accumulate the kernel-binned pair statistic of sf at the nodes' separations into sums, and the kernel-weighted pair count into counts, from spherical-harmonic pseudo-coefficients of the masked, weighted field. x is (lon, lat) per point in the metric's angle unit; data is the packed field with vector fields in (east, north[, radial]) components; weights one weight per point; valid which points hold a datum. spectral_backend names what computes the pseudo-coefficients: the direct sum, O(N lmax²), or a fast spherical harmonic transform from an extension.

For a spin-s monomial F at one end and a spin-s′ monomial G at the other, rotated into the pair's geodesic frame,

Σ_{ij} w_i w_j F̄_i conj(Ḡ_j) K^{ss′}(γ_ij, β) = (−1)^s Σ_l (2l+1)/(4π) b_l X^{FG}_l d^l_{ss′}(β),
K^{ss′}(γ, β) = (1/16π²) Σ_l (2l+1) b_l d^l_{ss′}(γ) d^l_{ss′}(β),

with X^{FG}_l = (1/(2l+1)) Σ_m F̃_lm conj(G̃_lm) the cross pseudo-spectrum, so every polynomial operator is a finite sum of such series and the counts are the mask's own Legendre series. The kernel tends to (1/8π²) δ(cos γ − cos β) as lmax grows, so the ratio sums ./ counts tends to the hard-binned pair average.

An operator odd in a scalar increment is refused: the kernel sum runs over both readings of every pair, so such a moment is identically zero here.

helmholtz_decompose_2d(distance_bins, sums, counts)

Run the 2D isotropic Helmholtz decomposition from native single-pass sums/counts. Rows 2 and 3 must contain L2SF and T2SF.

helmholtz_decompose_2d(distance_bins, L2_sums, L2_counts, T2_sums, T2_counts)

Run the 2D isotropic Helmholtz decomposition using the trapezoidal rule over binned longitudinal/transverse second-order structure functions. This implements the cumulative integral equations described by Lindborg (JFM 2015) and Bühler, Callies, and Ferrari (JFM 2014):

$D_rot(r) = D_TT(r) + I(r)$, $D_div(r) = D_LL(r) - I(r)$, where $I(r) = ∫_0^r [D_TT(s) - D_LL(s)]/s ds$.

The quadrature starts at the first bin midpoint, so $I$ omits $∫_0^{r_1}$. That segment is D_LL(r_1) for an $r^{2/3}$ inertial range, so the decomposition is quantitative only for $r ≫ r_1$; choose a first bin well below the scales of interest.

helmholtz_spectra(L2, T2, wavenumbers; asymptote)

Rotational and divergent spectral densities of a two-dimensional field from its longitudinal and transverse second-order structure functions, binned on the same edges.

The sum and the difference of the two projections transform separately,

\[P_E + P_B = -\frac{1}{4π} ∫_0^∞ (D_{LL} + D_{TT} - a)\, J_0(kr)\, r\, dr, \qquad P_E - P_B = \frac{1}{4π} ∫_0^∞ (D_{LL} - D_{TT})\, J_2(kr)\, r\, dr,\]

with E the divergent (gradient) part and B the rotational (curl) part, on the convention that a density integrates over d²k to the variance. The first line is isotropic_spectrum of the trace and needs the trace's large-separation limit a, which defaults to its largest value; the second needs none, since D_LL − D_TT decays on its own. Both kernels need Bessels.

Its error is the truncation of the Hankel integrals at the last separation, which is a different error from the one the real-space route through helmholtz_decompose_2d carries.

helmholtz_spectra(L2, T2, geometry::SphericalGeometry, lmax; variance) -> (l, E, B)

Gradient (E) and curl (B) angular power spectra, l = 1:lmax, of a tangent vector field on a sphere, from its longitudinal and transverse second-order structure functions binned in r = R σ.

In the geodesic frame ξ₊ = C_LL + C_TT = Σ_l (2l+1)/(4π) (C^E_l + C^B_l) d^l_{11}(σ) and ξ₋ = C_LL − C_TT = −Σ_l (2l+1)/(4π) (C^E_l − C^B_l) d^l_{1,−1}(σ), and the Wigner functions are orthogonal with ∫₀^π d^l_{ss'} d^{l'}_{ss'} sin σ dσ = 2δ_{ll'}/(2l+1), so

\[C^E_l - C^B_l = π ∫_0^π (D_{LL} - D_{TT})\, d^l_{1,-1}\, \sin σ\, dσ, \qquad C^E_l + C^B_l = 2π ∫_0^π \Big(⟨‖u‖²⟩ - \frac{D_{LL} + D_{TT}}{2}\Big)\, d^l_{11}\, \sin σ\, dσ .\]

The difference needs no constant because C_LL(0) = C_TT(0); the sum needs the field's mean square ⟨‖u‖²⟩, which a structure function cannot supply, so variance is required. Bins are integrated in closed form as in isotropic_spectrum on a sphere, nodes by their quadrature weights.

helmholtz_spectra(h, wavenumbers; rotational_asymptote, divergent_asymptote)

Rotational and divergent kinetic-energy spectra from a two-dimensional Helmholtz decomposition.

Each component of the decomposition is the share of the second-order trace carried by one of the two fields, so each transforms to a spectrum by the same route the trace does — in two dimensions, with the J₀ kernel, hence Bessels. Returns a NamedTuple of the two densities, on the convention isotropic_spectrum documents.

Bins holding no pair are dropped. The asymptotes default to each component's own largest value.

independent_pair_variance(joint::StructureFunction2DSumsAndCounts) -> Vector

The variance of each distance bin's mean pair value under the assumption that the pairs are independent: the variance of the pair values in the bin over the bin's count. Read from a value-binned joint histogram (InvariantValueAxis), each value cell contributing its own mean, so the spread inside a value cell is not seen and the result is a lower bound. A bin holding no pair is NaN. This is the data covariance W the fits take.

isotropic_spectrum(result, geometry::SphericalGeometry, lmax) -> (l, C)

Angular power spectrum C_l, l = 1:lmax, of a scalar or of the trace on a sphere, from a second-order structure function binned in separation r = R σ.

The isotropic correlation is C(σ) = Σ_l (2l+1)/(4π) C_l P_l(cos σ) and D(σ) = 2[C(0) − C(σ)], so by the orthogonality of the Legendre polynomials

\[C_l = -π ∫_0^π D(σ)\, P_l(\cos σ)\, \sin σ\, dσ, \qquad l ≥ 1,\]

with no variance needed: the constant C(0) is orthogonal to every P_l but P_0, which is the one degree lost — the sphere's k = 0. Each bin's value is taken as constant across the bin and the kernel is integrated over the bin in closed form, so the only error is that of the binning itself. Every bin must hold a value: the integral runs over the whole sphere. A kernel-binned result on HarmonicNodes integrates by the nodes' quadrature weights.

isotropic_spectrum(operator, separations, values, wavenumbers, ::Val{D}; asymptote)

Power spectral density at each of wavenumbers, from a second-order structure function sampled at separations.

D(r) = 2[C(0) - C(r)], so the transform of D differs from that of -2C only by a constant, whose transform is confined to k = 0. Every returned wavenumber must therefore be nonzero, and the k = 0 mode is not recoverable — no variance argument would help, because D does not carry it.

asymptote is the large-separation limit of the structure function, subtracted so the integrand decays; it defaults to the largest value supplied.

Normalised so that integrating the density over d^D k returns the field's variance, which fixes the convention: ∫₀^∞ shell_spectrum(...) dk == var(u).

isotropic_spectrum(result, wavenumbers, ::Val{D}; asymptote)

Spectral density from a structure function result, taking the operator, the separations and the values from the result itself.

The abscissa is the bin representative of result's edges. Bins holding no pair are dropped rather than carried as NaN, which would otherwise propagate through the quadrature into every wavenumber.

The transform averages over the directions of the separation, so it assumes the pairs behind each bin sample direction uniformly. Scattered points do; a rectilinear grid does not, and on gridded data the separations available at a given r are biased toward the lattice axes.

256-aligned joint compile width (capped at max). Implemented in GPU extension.

Exact joint histogram cell count n_dist × n_val. Implemented in GPU extension.

Compile-time joint 2D shared-histogram width SF_GPU_MAX_2D_HIST (4096). Implemented in GPU extension.

marginalize_sp2d_then_append_helmholtz_rows(sums_6, counts_6, distance_bins)

Marginalize six native invariant 2D joint histograms to 1D, then append rotational/divergent rows via append_helmholtz_rotational_divergent_rows. Returns $(sums_8, counts_8)$.

nufft_half_support(tag::NonuniformFFTsSpectralBackend) -> Int

The kernel half-support M at which NonuniformFFTs' backwards Kaiser–Bessel kernel, of shape β = π γ M (2 − 1/σ) (Potts & Steidl 2003, eq. 5.12) at its oversampling σ = 2, reaches the tag's tolerance: its aliasing error is ε = exp(−π M √((2 − 1/σ)² − 1/σ²)) = exp(−2π M √(1 − 1/σ)), so M = ⌈−ln ε / (2π √(1 − 1/σ))⌉; γ ≤ 1 makes the realised error slightly better than the estimate. The tolerance is floored at eps(Float64), below which nothing is left to resolve, and M at 2. 1e-7 gives 4, 1e-12 gives 7, 1e-15 gives 8.

nufft_monomial_transforms(tag, schedule::ScatteredModesSchedule, data, valid, weights, keys, ::Val{Pm}; to = identity) -> Vector

Type-1 non-uniform FFTs of the masked, weighted monomials keys of the packed point field onto the schedule's mode grid, one array per key in the real-to-complex half-spectrum layout (M₁ ÷ 2 + 1, M₂, …), each weighted by the schedule's taper; the arrays in the family to returns. tag is a provider tag, NonuniformFFTsSpectralBackend or FINUFFTSpectralBackend; supplied by the provider's extension.

Release device buffers held by a GPUSFWorkspace (optional explicit free).

Zero device histogram buffers in a GPUSFWorkspace before the next launch.

Zero every accumulator held by a CPUSFWorkspace; the drivers do this per call.

segmented_spectrum(p, k, edges) -> E(k)

The continuous piecewise power law with parameters p = (b₁, α₁, …, α_S) on the segments edges: b_s k^{α_s} on segment s, b_{s+1} = b_s k_s^{α_s - α_{s+1}} at each breakpoint. k outside the segments takes the nearest end segment's law.

select_segments(sf, k_edges, S_range, ::Val{D}; W = nothing, slope_bounds = (-4.0, 1.0))
    -> (segments, misfits, fits)

Fit a SegmentedPowerLaw with each number of segments in S_range and report the one minimising the mean relative misfit mean(|S₂^{fit} - S₂| / S₂), with every misfit and every fit.

serial_calculate_structure_function!(sums, counts, sf, x, fields, distance_bins; kwargs...)

Accumulate a multi-field's pairs into the 1-D distance histogram.

A field of one vector field and no scalars is the array path — it forwards there, so a bare u and Fields(vectors = (u,)) give the same answer through the same kernel. Anything carrying more than one field takes this sweep, which builds each pair's whole increment and hands it to the operator; it is a scalar loop, where the single-field path is a SIMD compute/scatter split.

serial_calculate_structure_function(sf, x, u, distance_bins[, value_bins][, count_eltype]; kwargs...)

The point-list structure function on the serial CPU backend, returning the raw StructureFunctionSumsAndCounts (or the joint StructureFunction2DSumsAndCounts with value_bins). Takes the arguments of calculate_structure_function without backend; the public entry picks the output type.

serial_calculate_structure_function_tensor!(sums, counts, order, shape, x, u, distance_bins; distance_metric, axis)

The serial pair loop behind calculate_structure_function_tensor!: every pair's increment in its pair frame, its order-fold outer product added to the bin of the pair's separation (and of its angle when axis names one).

serial_calculate_structure_functions_single_pass!(sums, counts, x, u, distance_bins)

Six-invariant-type single-pass 1D batch (serial).

serial_calculate_structure_functions_single_pass_2d!(sums, counts, x, u, distance_bins, value_bins)

Six-invariant-type SP2D batch (serial); output (6, n_dist, n_val, batch…).

serial_calculate_structure_functions_single_pass_2d(x, u, distance_bins, value_bins, sums_3d, counts_3d; kwargs...)

Zero sums_3d and counts_3d, (6, n_bins, n_val) each, and accumulate the six invariants' joint distance × value histograms into them through calculate_structure_functions_single_pass_2d!.

shell_average(wavenumbers, density, edges)

Radially bin a Dg-dimensional spectral density onto edges, returning (midpoints, E) with E the shell-integrated spectrum: summing E over the shells returns what summing density over the wavenumber cells returns.

shell_spectrum(P, wavenumbers, ::Val{D})

Shell-integrated spectrum E(k) = Ω_D k^(D-1) P(k) from a power spectral density.

E integrates over k to the variance the density integrates over d^D k to, which is the form the inertial-range scaling laws are stated in.

spectral_flux(L3, S3, wavenumbers)

The L3 route on two result objects sharing their distance bins, over the bins both hold.

spectral_flux(operator, separations, values, wavenumbers)

Interscale flux at each of wavenumbers, from an advective structure function sampled at separations.

Π_K = -(K/2) ∫₀^∞ SF_A(r) J₁(Kr) dr

SF_A is the advective structure function averaged over the directions of the separation, so the same uniform-direction assumption isotropic_spectrum carries applies here. The relation itself assumes no isotropy of the flow, which is what these estimators are for. The integral is the trapezoid rule over the samples from the origin, where the integrand vanishes, to the last separation.

The kernel's first peak sets which separations carry the most weight at a given wavenumber: J₁ peaks at Kr ≈ 1.84, so the flux at K is reported on mostly by separations near 1.84/K. That is a statement about weighting, not about where Π itself is largest — the explicit factor of K means |Π| keeps growing with K for a fixed feature.

This is the J₁ relation, which takes the advective structure function and assumes no isotropy of the flow. The companion relations on third-order structure functions are the methods on S3SFType, L3SFType and MixedSFType{1,0,2}, each with the boundary term its integration by parts leaves at the last separation.

spectral_flux(::MixedSFType{1,0,2}, separations, S, wavenumbers)

Interscale flux of the variance of a scalar θ (enstrophy, for the vorticity) in a two-dimensional isotropic flow from S = ⟨δu_L (δθ)²⟩, sampled at separations up to R:

Π_K = -(K²/4) ∫₀^R S(r) J₂(Kr) dr - (K/4) S(R) J₁(KR)

The J₁ relation on ⟨δθ δ𝓐_θ⟩ = (1/2r) d(r S)/dr integrated by parts, boundary term kept.

spectral_flux(::L3SFType, separations, L3, S3, wavenumbers)

Interscale energy flux of a two-dimensional isotropic flow from L3 = ⟨δu_L³⟩, with S3 = ⟨δu_L ‖δu‖²⟩ on the same separations for the boundary term:

Π_K = -(K³/12) ∫₀^R L3(r) J₃(Kr) r dr - (K²/12) [ R L3(R) J₂(KR) + (3/K) S3(R) J₁(KR) ]

The S3 relation integrated by parts once more with S3 = (1/3r²) d(r³ L3)/dr, which holds for an isotropic incompressible two-dimensional flow. Both boundary terms are kept. J₃ peaks at Kr ≈ 4.20.

spectral_flux(::S3SFType, separations, S3, wavenumbers)

Interscale energy flux of a two-dimensional isotropic flow from S3 = ⟨δu_L ‖δu‖²⟩, sampled at separations up to R:

Π_K = -(K²/4) ∫₀^R S3(r) J₂(Kr) dr - (K/4) S3(R) J₁(KR)

The J₁ relation integrated by parts with SF_A = (1/2r) d(r S3)/dr, the isotropic divergence of ⟨δu ‖δu‖²⟩. The boundary term at R does not vanish at any finite R and is kept. J₂ peaks at Kr ≈ 3.05, so the flux at K is reported on mostly by separations near 3.05/K.

threaded_calculate_structure_function!(sums, counts, sf, x, u, distance_bins[, value_bins]; kwargs...)

The in-place form of threaded_calculate_structure_function, accumulating into sums and counts; supplied by the OhMyThreads extension.

threaded_calculate_structure_function(sf, x, u, distance_bins[, value_bins][, count_eltype]; kwargs...)

The point-list structure function on the threaded CPU backend, returning the raw sums and counts. Takes the arguments of calculate_structure_function without backend; supplied by the OhMyThreads extension, so using OhMyThreads is required.

tradeoff_curve(model, y, W, priors) -> (misfit, norm)

For each prior in priors (each a prior variance shared by every value, or a vector of them), the regularised least-squares fit of y through model under the data covariance W, reporting the W-normalised misfit ‖W^{-1/2}(H x̂ - y)‖² and ‖x̂‖₂. The curve is returned for the caller to choose from; nothing is chosen.

AbstractFitMethod

How values on wavenumber bins are fitted to a structure function through a forward model.

AbstractForwardModel

A linear map y = H x from values on wavenumber bins to a structure function at fixed separations. forward_matrix returns H; model.r the separations, model.k_edges the bin edges.

CPUSFWorkspace{kind}(x, u, distance_bins[, value_bins]; backend = SerialBackend(), count_eltype = UInt32)

Reusable CPU scratch for the batch drivers: the batch-leading working copies of x/u, one accumulator per task, and the full-width reduction accumulator. Pass to a batch entry point as workspace = ... to make repeated calls allocation-free.

Without it a batch call allocates one accumulator per task, tens of MiB at large B, on every call; the GC pauses that follow land on some calls and widen the spread of call times.

kind matches GPUSFWorkspace: :sf1d, :joint2d, :single_pass or :single_pass_2d. It is a type parameter, so the accumulator rank is fixed at construction and every field is concretely typed.

Built from the same arguments as the call it serves — including backend, which fixes how many accumulators it holds — and checked against them on use: a workspace built for a different shape, element type or backend is a hard error, never a silent reallocation. Composes with BatchLeading — when the input is already batch-leading the transpose buffers have length zero and are never touched.

Not thread-safe: one workspace serves one call at a time, exactly like GPUSFWorkspace.

FINUFFTSpectralBackend(; tolerance = 1e-12)

The non-uniform FFT tag that names FINUFFT.jl as the provider — cuFINUFFT for points on a CUDA device; using FINUFFT supplies it, with using CUDA for the device. tolerance is the relative accuracy the transforms are asked for.

FluxForwardModel(r, k_edges)

The third-order structure function S3 = ⟨δu_L ‖δu‖²⟩ of a two-dimensional isotropic flow at the separations r from a spectral energy flux that is piecewise constant in wavenumber: x = (ε, ξ₁, …, ξ_n) with F(k) = -ε + Σ_j ξ_j Δk_j H(k - k_j), k_j the bin centres and Δk_j the bin widths, so ε is the flux towards large scales below the first bin and ξ_j the injection density (flux gained per unit wavenumber) in bin j. From S3(r) = -4r ∫_0^∞ F(k) J_2(kr)/k\,dk,

\[S3(r_i) = 2 ε r_i - Σ_j 4\,\frac{ξ_j Δk_j}{k_j} J_1(k_j r_i) ,\]

exactly for that F. F at the bin centres, just past each jump, is G x with G[i, 1] = -1, G[i, j + 1] = Δk_j for j ≤ i; flux_matrix returns G. The sign is that of spectral_flux: positive towards small scales. Needs Bessels.

GPUSFWorkspace

Reusable device histogram buffers and cached distance-bin edge arrays for GPU structure-function launches. Construct once per (backend, distance_bins[, value_bins]) configuration; pass to gpu_calculate_structure_function(!) or slice drivers to avoid per-call KA.zeros allocation and repeated edge uploads.

Kinds (kind type parameter)

• :sf1d — 1D distance histogram, out_sums_dev/out_cnts_dev of rank 1
• :joint2d — distance × value joint histogram, rank 2 (exact n_dist × n_val smem by default; override via joint2d_compile_cells)
• :single_pass — six invariant 1D distance histograms (6, NB)
• :single_pass_2d — six invariant distance × value joint histograms (6, NB, n_val); on-chip modes (:shared, :typeplane) flush shared histograms directly to out_*; :direct uses lazy.partition_sums_dev / lazy.partition_counts_dev plus merge.

kind is a type parameter, so the histogram rank is fixed at construction and every field is concretely typed — matching CPUSFWorkspace. Use the matching constructor overload; reset_histogram!(ws) zeroes device outputs before each launch. The constructors live in StructureFunctionsKernelAbstractionsExt and need using KernelAbstractions.

HelmholtzForwardModel(r, k_edges)

The longitudinal and transverse second-order structure functions of a two-dimensional field, stacked as [D_LL; D_TT] at the separations r, from the gradient (E) and curl (B) shell spectra stacked as [E_E; E_B] on the bins k_edges:

\[D_{LL} + D_{TT} = 2 ∫ (E_E + E_B)\,[1 - J_0(kr)]\,dk, \qquad D_{LL} - D_{TT} = 2 ∫ (E_E - E_B)\, J_2(kr)\,dk ,\]

the relations helmholtz_spectra inverts. Needs Bessels.

NonNegativeLeastSquares()

Least squares with every fitted value constrained non-negative, by the Lawson–Hanson active-set algorithm on the whitened system. For a flux that is the monotone-flux model, which cannot resolve a sink; for a spectrum it is the constraint that a spectral density is non-negative. It returns no covariance.

NonuniformFFTsSpectralBackend(; tolerance = 1e-12)

The non-uniform FFT tag that names NonuniformFFTs.jl as the provider, on the CPU or on the device the points live on; using NonuniformFFTs supplies it. tolerance is the relative accuracy asked of the transforms; the provider's kernel half-support follows from it by nufft_half_support, and is also the least number of modes a direction of the schedule may have.

RegularizedLeastSquares(prior)

Regularised least squares with a Gaussian prior on the fitted values: for y = H x + e, e of covariance W and a prior covariance P,

\[\hat x = (Hᵀ W^{-1} H + P^{-1})^{-1} Hᵀ W^{-1} y, \qquad C_{xx} = (Hᵀ W^{-1} H + P^{-1})^{-1} .\]

prior is the prior variance of each value (a vector), a full prior covariance matrix, or nothing for no prior, which is ordinary weighted least squares. The data covariance W is required by the fits taking this method; independent_pair_variance supplies one from a joint histogram, and tradeoff_curve shows how the misfit and the norm of x̂ trade against the prior.

ScatteredModesSchedule(points, r_max, modes; taper = NoTaper())

Scattered points (Dg, N) in the box that is their extent padded by r_max along every direction, so that no pair within r_max wraps, mapped to [0, 2π)^Dg, with a periodic mode grid of modes per direction.

The transform engine takes the forward transforms of the masked, weighted monomials as type-1 non-uniform FFTs onto that grid, so every lag-space product is the pair sum Σ_ij w_i w_j m_i m_j μ_i ν_j K(θ_j − θ_i − 2πh/M) with the periodic kernel K(φ) = Π_d (1/M_d) Σ_k b_k e^{i k φ_d} of the mode set in place of a delta at lag h: a soft bin of width about L_d/M_d, whose sidelobes taper trades for width (NoTaper() is the Dirichlet kernel; GaussianTaper(σ), σ a length, weights mode k by exp(−|k|²σ²/2)). The statistic is not exact: it converges to the hard-binned pair sum as the modes grow, and no finite mode set reproduces it. Its pair count is the kernel-weighted pair mass, a floating-point number, so results carry ModeBinEdges and need floating-point counts. Every point's pair with itself deposits half a unit of that mass on the smallest lags, as the harmonic route's self term does; bins that start above the kernel's width exclude it. The accuracy of the transforms themselves is the provider's, set on its tag: NonuniformFFTsSpectralBackend or FINUFFTSpectralBackend.

SegmentedPowerLaw(segments; slope_bounds = (-4.0, 1.0))

A shell spectrum that is a power law on each of segments wavenumber ranges, log-uniform between the smallest and the largest wavenumber asked for, continuous at the breakpoints, with slopes held in slope_bounds and a non-negative amplitude: S + 1 parameters (b₁, α₁, …, α_S). Fitted by bounded Levenberg–Marquardt through LsqFit (using LsqFit), on the relative residual (S₂^{fit} - S₂)/S₂ when no data covariance is given.

SpectrumForwardModel(::Val{D}, r, k_edges)

The second-order trace at the separations r from a shell spectrum E(k) constant on each of the wavenumber bins k_edges, in D dimensions:

\[S_2(r_i) = 2 ∫_0^∞ E(k)\,[1 - \mathrm{kernel}_D(k r_i)]\,dk = Σ_j H_{ij} E_j, \qquad H_{ij} = 2 ∫_{k_j}^{k_{j+1}} [1 - \mathrm{kernel}_D(k r_i)]\,dk ,\]

kernel_D the angular average of cos(k·r) (isotropic_kernel; J₀ in two dimensions, which needs Bessels). E integrates over k to the variance, the convention of shell_spectrum. Each bin is integrated by Gauss–Legendre quadrature.