API reference

AbstractComponentEmulators

Abstract type for component-level emulators that represent individual parts of the power spectrum calculation (e.g., P11, Ploop, Pct).

All concrete subtypes must implement the necessary fields to enable neural network evaluation, normalization, and postprocessing.

ComponentEmulator <: AbstractComponentEmulators

A complete emulator for a single power spectrum component, combining neural network predictions with normalization and physics-based postprocessing.

Fields

• TrainedEmulator::AbstractTrainedEmulators: The trained neural network (Lux or SimpleChains).
• kgrid::Array: Wavenumber grid on which the component is evaluated (in h/Mpc).
• InMinMax::Matrix{Float64}: Min-max normalization parameters for inputs (n_params × 2).
• OutMinMax::Array{Float64}: Min-max normalization parameters for outputs (n_k × 2).
• Postprocessing::Function: Function to apply physics transformations to raw NN output.

Details

The typical evaluation flow is:

1. Normalize input parameters using InMinMax.
2. Evaluate neural network to get normalized output.
3. Denormalize output using OutMinMax.
4. Apply postprocessing (e.g., multiply by D² for P11).

Example Postprocessing

# For linear power spectrum component
postprocess_P11 = (params, output, D, emu) -> output .* D^2

get_component(input_params, D, comp_emu::AbstractComponentEmulators)

Evaluate a component emulator to obtain power spectrum component values.

Arguments

• input_params: Array of input parameters (e.g., cosmological parameters).
• D: Growth factor at the redshift of interest.
• comp_emu::AbstractComponentEmulators: The component emulator to evaluate.

Returns

A matrix of shape (n_k, n_samples) containing the evaluated power spectrum component values on the emulator's k-grid.

Details

This function performs the full evaluation pipeline:

1. Copy input parameters to avoid mutation.
2. Apply min-max normalization to inputs.
3. Run neural network inference.
4. Denormalize network output.
5. Apply component-specific postprocessing (using D and emulator metadata).
6. Reshape to match k-grid dimensions.

The postprocessing step typically includes physics-based transformations such as scaling by powers of the growth factor.

AbstractPℓEmulators

Abstract type for complete power spectrum multipole emulators.

Concrete subtypes must combine multiple component emulators (P11, Ploop, Pct) with bias models to compute full power spectrum multipoles $P_\ell(k)$ for $\ell \in \{0, 2, 4\}$.

PℓEmulator <: AbstractPℓEmulators

Complete emulator for power spectrum multipoles in the Effective Field Theory of Large Scale Structure (EFTofLSS) framework.

Fields

• P11::ComponentEmulator: Emulator for the linear theory power spectrum component.
• Ploop::ComponentEmulator: Emulator for the one-loop corrections.
• Pct::ComponentEmulator: Emulator for the counterterm contributions.
• StochModel::Function: Function to compute stochastic (shot noise) terms.
• BiasCombination::Function: Function mapping bias parameters to linear combination weights.
• JacobianBiasCombination::Function: Analytical Jacobian of BiasCombination w.r.t. bias parameters.

Details

The power spectrum multipole is computed as:

\[P_\ell(k) = \sum_i c_i(b_1, b_2, ...) P_i(k)\]

where:

• $P_i(k)$ are the component power spectra (P11, Ploop, Pct, stochastic terms)
• $c_i(b_1, b_2, ...)$ are coefficients from the bias expansion

The BiasCombination function encodes the EFT bias model, while JacobianBiasCombination provides analytical derivatives for efficient gradient-based inference.

Example Usage

# Load pre-trained emulator
emu = trained_emulators["PyBirdmnuw0wacdm"]["0"]  # monopole

# Evaluate power spectrum
cosmology = [z, ln10As, ns, H0, ωb, ωcdm, mν, w0, wa]
bias = [b1, b2, b3, b4, b5, b6, b7, f, cϵ0, cϵ1, cϵ2]
D = 0.8  # growth factor

P0 = get_Pℓ(cosmology, D, bias, emu)

See Also

• get_Pℓ: Evaluate the power spectrum.
• get_Pℓ_jacobian: Evaluate power spectrum and its Jacobian.

get_Pℓ_jacobian(cosmology::Array, D, bs::Array, cosmoemu::AbstractPℓEmulators; stoch_kwargs...)

Compute both the power spectrum multipole $P_\ell(k)$ and its Jacobian with respect to bias parameters.

Arguments

• cosmology::Array: Array of cosmological parameters (format depends on the emulator training).
• D: Growth factor value at the redshift of interest.
• bs::Array: Array of bias parameters.
• cosmoemu::AbstractPℓEmulators: The multipole emulator containing P11, Ploop, Pct components.

Keyword Arguments

• stoch_kwargs...: Additional keyword arguments passed to the stochastic model (e.g., shot noise parameters).

Returns

A tuple (Pℓ, ∂Pℓ_∂b) where:

• Pℓ: Power spectrum multipole values evaluated on the emulator's k-grid.
• ∂Pℓ_∂b: Jacobian matrix of the power spectrum with respect to bias parameters.

Details

This function is optimized for inference workflows where both the power spectrum and its derivatives are needed (e.g., gradient-based MCMC, Fisher forecasts). It computes both quantities in a single pass, avoiding redundant neural network evaluations.

The Jacobian is computed using the analytical derivative of the bias combination function, which is significantly faster than automatic differentiation for this specific operation.

See Also

• get_Pℓ: Compute only the power spectrum without Jacobian.

_Pkμ(k, μ, Int_Mono, Int_Quad, Int_Hexa)

Reconstructs the anisotropic power spectrum $P(k, \mu)$ at a given wavenumber k and cosine of the angle to the line-of-sight μ, using its Legendre multipole moments.

Arguments

• k: The wavenumber.
• μ: The cosine of the angle to the line-of-sight.
• Int_Mono: A function or interpolant that provides the monopole moment $I_0(k)$ at wavenumber k.
• Int_Quad: A function or interpolant that provides the quadrupole moment $I_2(k)$ at wavenumber k.
• Int_Hexa: A function or interpolant that provides the hexadecapole moment $I_4(k)$ at wavenumber k.

Returns

The value of the anisotropic power spectrum $P(k, \mu)$ at the given k and μ.

Details

The anisotropic power spectrum is reconstructed as a sum of its multipole moments multiplied by the corresponding Legendre polynomials evaluated at μ. The function uses the 0th, 2nd, and 4th order Legendre polynomials.

Formula

The formula used is:

\[P(k, \mu) = I_0(k) \mathcal{L}_0(\mu) + I_2(k) \mathcal{L}_2(\mu) + I_4(k) \mathcal{L}_4(\mu)\]

where $I_l(k)$ are the multipole moments and $\mathcal{L}_l(\mu)$ are the Legendre polynomials of order $l$.

See Also

• _Legendre_0: Calculates the 0th order Legendre polynomial.
• _Legendre_2: Calculates the 2nd order Legendre polynomial.
• _Legendre_4: Calculates the 4th order Legendre polynomial.

_k_true(k_o, μ_o, q_perp, F)

Calculates the true (physical) wavenumber k from the observed wavenumber k_o and observed cosine of the angle to the line-of-sight μ_o.

This transformation accounts for anisotropic effects, likely redshift-space distortions (RSD) or anisotropic cosmological scaling, parameterized by q_perp and F.

Arguments

• k_o: The observed wavenumber (scalar).
• μ_o: The observed cosine of the angle to the line-of-sight (scalar).
• q_perp: A parameter related to perpendicular anisotropic scaling.
• F: A parameter related to parallel anisotropic scaling (often the growth rate f divided by the anisotropic scaling parameter q_parallel).

Returns

The calculated true wavenumber k (scalar).

Formula

The formula used is:

\[k = \frac{k_o}{q_\perp} \sqrt{1 + \mu_o^2 \left(\frac{1}{F^2} - 1\right)}\]

See Also

• _k_true(k_o::Array, μ_o::Array, q_perp, F): Method for arrays of observed values.
• _μ_true: Calculates the true cosine of the angle to the line-of-sight.

_k_true(k_o::Array, μ_o::Array, q_perp, F)

Calculates the true (physical) wavenumber k for arrays of observed wavenumbers k_o and observed cosines of the angle to the line-of-sight μ_o.

This method applies the transformation from observed to true wavenumber element-wise or for combinations of input arrays, accounting for anisotropic effects parameterized by q_perp and F.

Arguments

• k_o: An array of observed wavenumbers.
• μ_o: An array of observed cosines of the angle to the line-of-sight.
• q_perp: A parameter related to perpendicular anisotropic scaling.
• F: A parameter related to parallel anisotropic scaling.

Returns

A vector containing the calculated true wavenumbers k for the given input arrays.

Details

The function calculates k for pairs or combinations of values from the input arrays k_o and μ_o using a formula derived from anisotropic scaling. The calculation involves broadcasting and array operations to handle the array inputs efficiently. The result is flattened into a vector.

Formula

The underlying transformation for each pair of k_o and μ_o is:

\[k = \frac{k_o}{q_\perp} \sqrt{1 + \mu_o^2 \left(\frac{1}{F^2} - 1\right)}\]

See Also

• _k_true(k_o, μ_o, q_perp, F): Method for scalar observed values.
• _μ_true: Calculates the true cosine of the angle to the line-of-sight.

_μ_true(μ_o, F)

Calculates the true (physical) cosine of the angle to the line-of-sight μ from the observed cosine of the angle to the line-of-sight μ_o.

This transformation accounts for anisotropic effects, likely redshift-space distortions (RSD) or anisotropic cosmological scaling, parameterized by F.

Arguments

• μ_o: The observed cosine of the angle to the line-of-sight (scalar).
• F: A parameter related to parallel anisotropic scaling (often the growth rate f divided by the anisotropic scaling parameter q_parallel).

Returns

The calculated true cosine of the angle to the line-of-sight μ (scalar).

Formula

The formula used is:

\[\mu = \frac{\mu_o}{F \sqrt{1 + \mu_o^2 \left(\frac{1}{F^2} - 1\right)}}\]

See Also

• _μ_true(μ_o::Array, F): Method for an array of observed values.
• _k_true: Calculates the true wavenumber.

_μ_true(μ_o::Array, F)

Calculates the true (physical) cosine of the angle to the line-of-sight μ for an array of observed cosines of the angle to the line-of-sight μ_o.

This method applies the transformation from observed to true angle cosine element-wise, accounting for anisotropic effects parameterized by F.

Arguments

• μ_o: An array of observed cosines of the angle to the line-of-sight.
• F: A parameter related to parallel anisotropic scaling.

Returns

An array containing the calculated true cosines of the angle to the line-of-sight μ.

Details

The function calculates μ for each value in the input array μ_o using a formula derived from anisotropic scaling. Broadcasting (@.) is used to apply the calculation element-wise.

Formula

The underlying transformation for each μ_o is:

\[\mu = \frac{\mu_o}{F \sqrt{1 + \mu_o^2 \left(\frac{1}{F^2} - 1\right)}}\]

See Also

• _μ_true(μ_o, F): Method for a scalar observed value.
• _k_true: Calculates the true wavenumber.

_P_obs(k_o, μ_o, q_par, q_perp, Int_Mono, Int_Quad, Int_Hexa)

Calculates the observed power spectrum $P_{\text{obs}}(k_o, \mu_o)$ at a given observed wavenumber k_o and observed cosine of the angle to the line-of-sight μ_o.

This function transforms the observed coordinates to true (physical) coordinates, calculates the true power spectrum using provided interpolants for the multipole moments, and applies the appropriate scaling factor due to anisotropic effects.

Arguments

• k_o: The observed wavenumber.
• μ_o: The observed cosine of the angle to the line-of-sight.
• q_par: A parameter related to parallel anisotropic scaling.
• q_perp: A parameter related to perpendicular anisotropic scaling.
• Int_Mono: An interpolation function for the monopole moment $I_0(k)$ in true k.
• Int_Quad: An interpolation function for the quadrupole moment $I_2(k)$ in true k.
• Int_Hexa: An interpolation function for the hexadecapole moment $I_4(k)$ in true k.

Returns

The value of the observed power spectrum $P_{\text{obs}}(k_o, \mu_o)$.

Details

The observed coordinates $(k_o, \mu_o)$ are transformed to true coordinates $(k_t, \mu_t)$ using the _k_true and _μ_true functions, with $F = q_\parallel / q_\perp$. The true power spectrum $P(k_t, \mu_t)$ is then reconstructed using _Pkμ and the provided multipole interpolants. Finally, the result is scaled by $1 / (q_\parallel q_\perp^2)$.

Formula

The formula used is:

\[P_{\text{obs}}(k_o, \mu_o) = \frac{1}{q_\parallel q_\perp^2} P(k_t, \mu_t)\]

where

\[k_t = \text{_k_true}(k_o, \mu_o, q_\perp, F)\]

\[\mu_t = \text{_μ_true}(\mu_o, F)\]

and

\[F = q_\parallel / q_\perp\]

See Also

• _k_true: Transforms observed wavenumber to true wavenumber.
• _μ_true: Transforms observed angle cosine to true angle cosine.
• _Pkμ: Reconstructs the true power spectrum from multipole moments.
• interp_Pℓs: Creates the multipole interpolants.

interp_Pℓs(Mono_array, Quad_array, Hexa_array, k_grid)

Creates interpolation functions for the monopole, quadrupole, and hexadecapole moments of the power spectrum.

These interpolants can then be used to efficiently evaluate the multipole moments at arbitrary wavenumbers k.

Arguments

• Mono_array: An array containing the values of the monopole moment $I_0(k)$.
• Quad_array: An array containing the values of the quadrupole moment $I_2(k)$.
• Hexa_array: An array containing the values of the hexadecapole moment $I_4(k)$.
• k_grid: An array containing the corresponding wavenumber k values for the multipole arrays.

Returns

A tuple containing three interpolation functions: (Int_Mono, Int_Quad, Int_Hexa).

Details

The function uses AkimaInterpolation from the Interpolations.jl package to create the interpolants. Extrapolation is set to ExtrapolationType.Extension, which means the interpolant will use the nearest data points to extrapolate outside the provided k_grid range. Note that extrapolation can sometimes introduce errors.

See Also

• _Pkμ: Uses the interpolation functions to reconstruct the anisotropic power spectrum.

q_par_perp(z, cosmo_mcmc::AbstractCosmology, cosmo_ref::AbstractCosmology)

Calculates the parallel (q_par) and perpendicular (q_perp) Alcock-Paczynski (AP) parameters at a given redshift z, comparing a varying cosmology to a reference cosmology.

The AP parameters quantify the distortion of observed clustering due to assuming a different cosmology than the true one when converting redshifts and angles to distances.

Arguments

• z: The redshift at which to calculate the AP parameters.
• cosmo_mcmc: An AbstractCosmology struct representing the varying cosmology (e.g., from an MCMC chain).
• cosmo_ref: An AbstractCosmology struct representing the reference cosmology used for measurements.

Returns

A tuple (q_par, q_perp) containing the calculated parallel and perpendicular AP parameters at redshift z.

Details

The parallel AP parameter q_par is the ratio of the Hubble parameter in the reference cosmology to that in the varying cosmology. The perpendicular AP parameter q_perp is the ratio of the conformal angular diameter distance in the varying cosmology to that in the reference cosmology.

Formula

The formulas for the Alcock-Paczynski parameters are:

\[q_\parallel(z) = \frac{E_{\text{ref}}(z)}{E_{\text{mcmc}}(z)}\]

\[q_\perp(z) = \frac{\tilde{d}_{A,\text{mcmc}}(z)}{\tilde{d}_{A,\text{ref}}(z)}\]

where $E(z)$ is the normalized Hubble parameter and $\tilde{d}_A(z)$ is the conformal angular diameter distance.

apply_AP_check(k_input::AbstractVector, k_output::AbstractVector, Mono_array::AbstractVector, Quad_array::AbstractVector, Hexa_array::AbstractVector, q_par, q_perp)

Calculates the observed power spectrum multipole moments (monopole, quadrupole, hexadecapole) on a given observed wavenumber grid k_output, from vectors of true multipole moments provided on an input wavenumber grid k_input, using numerical integration.

This is a check version, intended for verifying results from faster methods. It is significantly slower due to the use of numerical integration over the angle μ.

Arguments

• k_input: A vector of wavenumber values on which the input true multipole moments (Mono_array, Quad_array, Hexa_array) are defined.
• k_output: A vector of observed wavenumber values at which to calculate the output observed multipoles.
• Mono_array: A vector containing the values of the true monopole moment $I_0(k)$ on the k_input grid.
• Quad_array: A vector containing the values of the true quadrupole moment $I_2(k)$ on the k_input grid.
• Hexa_array: A vector containing the values of the true hexadecapole moment $I_4(k)$ on the k_input grid.
• q_par: A parameter related to parallel anisotropic scaling.
• q_perp: A parameter related to perpendicular anisotropic scaling.

Returns

A tuple (P0_obs, P2_obs, P4_obs), where each element is a vector containing the calculated observed monopole, quadrupole, and hexadecapole moments respectively, evaluated at the wavenumbers in k_output.

Details

This method first creates interpolation functions for the true multipole moments using interp_Pℓs based on the k_input grid. It then calls the core apply_AP_check(k_grid, int_Mono, int_Quad, int_Hexa, q_par, q_perp) method, passing k_output as the grid at which to calculate the observed multipoles.

This function is a slower check implementation and should not be used in performance-critical code.

Formula

The observed multipole moments are calculated using the formula:

\[P_\ell(k_o) = (2\ell + 1) \int_{0}^1 P_{\text{obs}}(k_o, \mu_o) \mathcal{L}_\ell(\mu_o) d\mu_o\]

for $\ell \in \{0, 2, 4\}$. The observed power spectrum $P_{\text{obs}}(k_o, \mu_o)$ is calculated using _P_obs(k_o, μ_o, q_par, q_perp, int_Mono, int_Quad, int_Hexa).

See Also

• apply_AP_check(k_grid, int_Mono, int_Quad, int_Hexa, q_par, q_perp): The core method performing the integration.
• interp_Pℓs: Creates the interpolation functions for the true multipoles.
• _P_obs: Calculates the observed power spectrum.

apply_AP(k_input::AbstractVector, k_output::AbstractVector, mono::AbstractVector, quad::AbstractVector, hexa::AbstractVector, q_par, q_perp; n_GL_points=8)

Calculates the observed power spectrum multipole moments (monopole, quadrupole, hexadecapole) on a given observed wavenumber grid k_output, using arrays of true multipole moments provided on an input wavenumber grid k_input, and employing Gauss-Lobatto quadrature.

This is the standard, faster implementation for applying the Alcock-Paczynski (AP) effect to the power spectrum multipoles, designed for performance compared to the check version using generic numerical integration.

Arguments

• k_input: A vector of wavenumber values on which the input true multipole moments (mono, quad, hexa) are defined.
• k_output: A vector of observed wavenumber values at which to calculate the output observed multipoles.
• mono: A vector containing the values of the true monopole moment $I_0(k)$ on the k_input grid.
• quad: A vector containing the values of the true quadrupole moment $I_2(k)$ on the k_input grid.
• hexa: A vector containing the values of the true hexadecapole moment $I_4(k)$ on the k_input grid.
• q_par: A parameter related to parallel anisotropic scaling.
• q_perp: A parameter related to perpendicular anisotropic scaling.

Keyword Arguments

• n_GL_points: The number of Gauss-Lobatto points to use for the integration over μ. The actual number of nodes used corresponds to 2 * n_GL_points. Defaults to 8.

Returns

A tuple (P0_obs, P2_obs, P4_obs), where each element is a vector containing the calculated observed monopole, quadrupole, and hexadecapole moments respectively, evaluated at the observed wavenumbers in k_output.

Details

The function applies the AP and RSD effects by integrating the observed anisotropic power spectrum $P_{\text{obs}}(k_o, \mu_o)$ over the observed cosine of the angle to the line-of-sight $\mu_o \in [0, 1]$ (assuming symmetry for even multipoles), weighted by the corresponding Legendre polynomial $\mathcal{L}_\ell(\mu_o)$.

The process involves:

1. Determine Gauss-Lobatto nodes and weights for the interval [0, 1].
2. For each observed wavenumber k_o in the input k_output array and each μ_o node: a. Calculate the true wavenumber $k_t(k_o, \mu_o)$ using _k_true. b. Calculate the true angle cosine $\mu_t(\mu_o)$ using _μ_true. c. Interpolate the true multipole moments $I_\ell(k_t)$ using _akima_interpolation, interpolating from the k_input grid to the new k_t values. d. Calculate the true Legendre polynomials $\mathcal{L}_\ell(\mu_t)$ using _Legendre_0, _Legendre_2, _Legendre_4. e. Reconstruct the true power spectrum $P(k_t, \mu_t)$ using _Pk_recon. f. Calculate the observed power spectrum $P_{\text{obs}}(k_o, \mu_o) = P(k_t, \mu_t) / (q_\parallel q_\perp^2)$.
3. Perform the weighted sum (quadrature) over the μ_o nodes to get the observed multipoles $P_\ell(k_o)$ on the k_output grid.

This function is the standard, performant implementation for applying AP compared to the slower apply_AP_check.

Formula

The observed multipole moments are calculated using the formula:

\[P_\ell(k_o) = (2\ell + 1) \int_{0}^1 P_{\text{obs}}(k_o, \mu_o) \mathcal{L}_\ell(\mu_o) d\mu_o\]

for $\ell \in \{0, 2, 4\}$. The integral is approximated using Gauss-Lobatto quadrature.

See Also

• apply_AP_check: The slower, check version using generic numerical integration.
• _k_true: Transforms observed wavenumber to true wavenumber.
• _μ_true: Transforms observed angle cosine to true angle cosine.
• _Legendre_0, _Legendre_2, _Legendre_4: Calculate the Legendre polynomials.
• _akima_interpolation: Interpolates the true multipole moments.
• _Pk_recon: Reconstructs the true power spectrum on a grid.
• gausslobatto: Function used to get quadrature nodes and weights.

apply_AP(k_input::AbstractVector, k_output::AbstractVector, mono::AbstractMatrix, quad::AbstractMatrix, hexa::AbstractMatrix, q_par, q_perp; n_GL_points=8)

Batch version of apply_AP for processing multiple columns simultaneously using optimized matrix Akima interpolation.

This method applies the Alcock-Paczynski effect to multiple sets of multipole moments (e.g., multiple Jacobian columns or parameter variations) in a single call. It leverages the optimized matrix Akima spline implementation to interpolate all columns at once, providing significant performance improvements over column-by-column processing.

Arguments

• k_input::AbstractVector: Input wavenumber grid.
• k_output::AbstractVector: Output wavenumber grid.
• mono::AbstractMatrix: Monopole moments with shape (n_k, n_cols).
• quad::AbstractMatrix: Quadrupole moments with shape (n_k, n_cols).
• hexa::AbstractMatrix: Hexadecapole moments with shape (n_k, n_cols).
• q_par: Parallel AP parameter.
• q_perp: Perpendicular AP parameter.

Keyword Arguments

• n_GL_points::Int: Number of Gauss-Lobatto points. Default: 8.

Returns

A tuple (mono_AP, quad_AP, hexa_AP) where each is a matrix of shape (n_k_output, n_cols) containing the AP-corrected multipoles for all input columns.

Details

This optimized implementation uses matrix Akima interpolation to process all columns simultaneously rather than iterating column-by-column. For N columns, this reduces 3×N Akima calls to just 3 matrix Akima calls, providing a ~2-3× speedup.

This is particularly useful for computing Jacobians where each column represents the derivative with respect to a different parameter (typically 11 bias parameters).

Performance

For typical DESI-like scenarios (50 input k-points, 100 output k-points, 11 columns):

• Old implementation: ~6 ms (33 scalar Akima calls)
• New implementation: ~2 ms (3 matrix Akima calls)
• Speedup: ~2.5-3×

See Also

• apply_AP(k_input::AbstractVector, k_output::AbstractVector, mono::AbstractVector, quad::AbstractVector, hexa::AbstractVector, q_par, q_perp): Single-column version.
• _akima_interpolation(u::AbstractMatrix, t, t_new): Optimized matrix Akima interpolation.

_Pk_recon(mono::Matrix, quad::Matrix, hexa::Matrix, l0, l2, l4)

Reconstructs the anisotropic power spectrum $P(k, \mu)$ on a grid of wavenumbers k and cosines of the angle to the line-of-sight μ, using matrices of its Legendre multipole moments and vectors of Legendre polynomial values.

This function is designed to efficiently reconstruct the 2D power spectrum for multiple k and μ values simultaneously, assuming the multipole moments are provided as matrices (e.g., N_k x 1) and Legendre polynomials as vectors (e.g., N_μ).

Arguments

• mono: A matrix containing the monopole moment $I_0(k)$ values (expected dimensions N_k x 1).
• quad: A matrix containing the quadrupole moment $I_2(k)$ values (expected dimensions N_k x 1).
• hexa: A matrix containing the hexadecapole moment $I_4(k)$ values (expected dimensions N_k x 1).
• l0: A vector containing the 0th order Legendre polynomial $\mathcal{L}_0(\mu)$ values evaluated at the desired μ values (expected dimensions N_μ).
• l2: A vector containing the 2nd order Legendre polynomial $\mathcal{L}_2(\mu)$ values evaluated at the desired μ values (expected dimensions N_μ).
• l4: A vector containing the 4th order Legendre polynomial $\mathcal{L}_4(\mu)$ values evaluated at the desired μ values (expected dimensions N_μ).

Returns

A matrix representing the anisotropic power spectrum $P(k, \mu)$ on the N_k x N_μ grid.

Details

The function reconstructs the anisotropic power spectrum using the formula that sums the multipole moments multiplied by the corresponding Legendre polynomials. The matrix and vector operations are broadcast to calculate the result for all combinations of input k (from the rows of the moment matrices) and μ (from the elements of the Legendre polynomial vectors).

Formula

The formula used for each element $(i, j)$ of the output matrix (corresponding to the $i$-th wavenumber and $j$-th angle cosine) is:

\[P(k_i, \mu_j) = I_0(k_i) \mathcal{L}_0(\mu_j) + I_2(k_i) \mathcal{L}_2(\mu_j) + I_4(k_i) \mathcal{L}_4(\mu_j)\]

See Also

• _Pkμ: Reconstructs $P(k, \mu)$ for single k and μ.
• _Legendre_0, _Legendre_2, _Legendre_4: Calculate the Legendre polynomials.

window_convolution(W::Array{T, 4}, v::Matrix) where {T}

Applies a 4-dimensional window function or kernel W to a 2-dimensional input matrix v.

This operation performs a transformation or generalized convolution, summing over the j and l indices of the inputs to produce a 2D result indexed by i and k. This is commonly used in analyses where a 4D kernel relates input data in two dimensions to output data in another two dimensions.

Arguments

• W: A 4-dimensional array representing the window function or kernel.
• v: A 2-dimensional matrix representing the input data.

Returns

A 2-dimensional matrix representing the result of the convolution or transformation.

Details

The function implements the summation using the @tullio macro, which provides an efficient way to express tensor contractions and generalized convolutions. The operation can be thought of as applying a 4D kernel to a 2D input, resulting in a 2D output.

Formula

The operation is defined as:

\[C_{ik} = \sum_{j,l} W_{ijkl} v_{jl}\]

See Also

• window_convolution(W::AbstractMatrix, v::AbstractVector): Method for a matrix kernel and vector input.

References

• The methodology for this type of window measurement is discussed in: arXiv:1810.05051

window_convolution(W::AbstractMatrix, v::AbstractVector)

Performs matrix-vector multiplication, where the matrix W acts as a linear transformation or window applied to the vector input v.

Arguments

• W: An abstract matrix representing the linear transformation or window.
• v: An abstract vector representing the input data.

Returns

An abstract vector representing the result of the matrix-vector multiplication.

Details

This method is a direct implementation of standard matrix-vector multiplication. It applies the linear transformation defined by matrix W to the vector v.

Formula

The operation is defined as:

\[\mathbf{c} = \mathbf{W} \mathbf{v}\]

or element-wise:

\[c_i = \sum_j W_{ij} v_j\]

See Also

• window_convolution(W::Array{T, 4}, v::Matrix) where {T}: Method for a 4D kernel and matrix input.

References

• The methodology for this type of window measurement is discussed in: arXiv:1810.05051

_transformed_weights(quadrature_rule, order, a, b)

Transforms the points and weights of a standard quadrature rule from the interval [-1, 1] to a specified interval [a, b].

This is a utility function used to adapt standard quadrature rules (like Gauss-Legendre) for numerical integration over arbitrary intervals [a, b].

Arguments

• quadrature_rule: A function that takes an order and returns a tuple (points, weights) for the standard interval [-1, 1].
• order: The order of the quadrature rule (number of points).
• a: The lower bound of the target interval.
• b: The upper bound of the target interval.

Returns

A tuple (transformed_points, transformed_weights) for the interval [a, b].

Details

The transformation is applied to the standard points $x_i^{\text{std}}$ and weights $w_i^{\text{std}}$ obtained from the quadrature_rule:

• Transformed points: $x_i = \frac{b - a}{2} x_i^{\text{std}} + \frac{b + a}{2}$
• Transformed weights: $w_i = \frac{b - a}{2} w_i^{\text{std}}$

Formula

The transformation formulas are: Points: $x_i = \frac{b - a}{2} x_i^{\text{std}} + \frac{b + a}{2}$ Weights: $w_i = \frac{b - a}{2} w_i^{\text{std}}$

_Legendre_0(x)

Calculates the 0th order Legendre polynomial, $\mathcal{L}_0(x)$.

Arguments

• x: The input value (typically the cosine of an angle, -1 ≤ x ≤ 1).

Returns

The value of the 0th order Legendre polynomial evaluated at x.

Formula

The formula for the 0th order Legendre polynomial is:

\[\mathcal{L}_0(x) = 1\]

See Also

• _Legendre_2: Calculates the 2nd order Legendre polynomial.
• _Legendre_4: Calculates the 4th order Legendre polynomial.
• _Pkμ: A function that uses Legendre polynomials.

_Legendre_2(x)

Calculates the 2nd order Legendre polynomial, $\mathcal{L}_2(x)$.

Arguments

• x: The input value (typically the cosine of an angle, -1 ≤ x ≤ 1).

Returns

The value of the 2nd order Legendre polynomial evaluated at x.

Formula

The formula for the 2nd order Legendre polynomial is:

\[\mathcal{L}_2(x) = \frac{1}{2} (3x^2 - 1)\]

See Also

• _Legendre_0: Calculates the 0th order Legendre polynomial.
• _Legendre_4: Calculates the 4th order Legendre polynomial.
• _Pkμ: A function that uses Legendre polynomials.

_Legendre_4(x)

Calculates the 4th order Legendre polynomial, $\mathcal{L}_4(x)$.

Arguments

• x: The input value (typically the cosine of an angle, -1 ≤ x ≤ 1).

Returns

The value of the 4th order Legendre polynomial evaluated at x.

Formula

The formula for the 4th order Legendre polynomial is:

\[\mathcal{L}_4(x) = \frac{1}{8} (35x^4 - 30x^2 + 3)\]

See Also

• _Legendre_0: Calculates the 0th order Legendre polynomial.
• _Legendre_2: Calculates the 2nd order Legendre polynomial.
• _Pkμ: A function that uses Legendre polynomials.

_akima_interpolation(u, t, t_new)

Evaluates the one-dimensional Akima spline that interpolates the data points $(t_i, u_i)$ at new abscissae t_new.

Arguments

• u: Ordinates (function values) $u_i$ at the data nodes.
• t: Strictly increasing abscissae (knots) $t_i$ associated with u. length(t) must equal length(u).
• t_new: The query point(s) where the spline is to be evaluated.

Returns

The interpolated value(s) at t_new. A scalar input returns a scalar; a vector input returns a vector of the same length.

Details

This routine implements the original Akima piecewise-cubic method (T. Akima, 1970). On each interval $[t_j, t_{j+1}]$, a cubic polynomial is constructed. The method uses a weighted average of slopes to determine the derivative at each node, which effectively dampens oscillations without explicit shape constraints. The resulting spline is $C^1$ continuous (its first derivative is continuous) but generally not $C^2$.

Formulae

The spline on the interval $[t_j, t_{j+1}]$ is a cubic polynomial: [ Sj(w) = uj + bj w + cj w^{2} + dj w^{3}, \qquad w = t - tj ] The derivative $b_j$ at each node is determined by Akima's weighting of local slopes $m_j=(u_{j}-u_{j-1})/(t_j-t_{j-1})$: [ bj = \frac{|m{j+1}-m{j}|\,m{j-1} + |m{j-1}-m{j-2}|\,m{j}} {|m{j+1}-m{j}| + |m{j-1}-m{j-2}|} ] The remaining coefficients, ``cj$and$dj``, are found by enforcing continuity of the first derivative: [ cj = \frac{3mj - 2bj - b{j+1}}{t{j+1}-tj} ] [ dj = \frac{bj + b{j+1} - 2mj}{(t{j+1}-t_j)^2} ]

Automatic Differentiation

The implementation is free of mutation on the inputs and uses only element-wise arithmetic, making the returned value differentiable with both ForwardDiff.jl (dual numbers) and Zygote.jl (reverse-mode AD). You can therefore embed _akima_interpolation in optimization or machine-learning pipelines and back-propagate through the interpolation seamlessly.

Notes

The algorithm and numerical results are equivalent to the Akima spline in DataInterpolations.jl, but this routine is self-contained and avoids any package dependency.

_akima_interpolation(u::AbstractMatrix, t, t_new)

Akima spline interpolation for multiple data series sharing the same x-coordinates. Uses a simple comprehension-based approach that is compatible with automatic differentiation.

Arguments

• u::AbstractMatrix: Data values with shape (n_points, n_columns).
• t: X-coordinates shared by all columns.
• t_new: Query points.

Returns

Matrix of interpolated values with shape (length(t_new), n_columns).

Example

# Interpolate 11 Jacobian columns at 100 k-points
k_in = range(0.01, 0.3, length=50)
k_out = range(0.01, 0.3, length=100)
jacobian = randn(50, 11)  # 11 parameters

result = _akima_interpolation(jacobian, k_in, k_out)  # (100, 11)

_akima_slopes(u::AbstractMatrix, t)

Optimized version of _akima_slopes for matrix input where each column represents a different data series but all share the same x-coordinates t.

Performance Optimization

Computes dt = diff(t) once and reuses it for all columns, avoiding redundant computation.

Arguments

• u::AbstractMatrix: Data values with shape (n_points, n_columns).
• t: X-coordinates (same for all columns).

Returns

Matrix of slopes with shape (n_points + 3, n_columns).

_akima_coefficients(t, m::AbstractMatrix)

Optimized version of _akima_coefficients for matrix input where each column represents coefficients for a different spline series.

Performance Optimization

Computes dt = diff(t) once and reuses it for all columns.

Arguments

• t: X-coordinates.
• m::AbstractMatrix: Slopes matrix with shape (n_points + 3, n_columns).

Returns

Tuple (b, c, d) where:

• b is a matrix of shape (n_points, n_columns)
• c and d are matrices of shape (n_points - 1, n_columns)

_akima_eval(u::AbstractMatrix, t, b::AbstractMatrix, c::AbstractMatrix, d::AbstractMatrix, tq::AbstractArray)

Optimized version of _akima_eval for matrix input where each column represents a different spline series.

Performance Optimization

• Finds intervals once per query point (not per column)
• Computes polynomial weights once per query point
• Broadcasts evaluation across all columns simultaneously

This is significantly faster than calling the vector version in a loop.

Arguments

• u::AbstractMatrix: Data values with shape (n_points, n_columns).
• t: X-coordinates.
• b::AbstractMatrix, c::AbstractMatrix, d::AbstractMatrix: Spline coefficients.
• tq::AbstractArray: Query points.

Returns

Matrix of interpolated values with shape (length(tq), n_columns).

load_component_emulator(path::String; emu=LuxEmulator, k_file="k.npy", weights_file="weights.npy", inminmax_file="inminmax.npy", outminmax_file="outminmax.npy", nn_setup_file="nn_setup.json", postprocessing_file="postprocessing_file.jl")

Load a trained component emulator from disk.

Arguments

• path::String: Directory path containing the emulator files.

Keyword Arguments

• emu: Emulator type to initialize (LuxEmulator or SimpleChainsEmulator). Default: LuxEmulator.
• k_file::String: Filename for the wavenumber grid. Default: "k.npy".
• weights_file::String: Filename for neural network weights. Default: "weights.npy".
• inminmax_file::String: Filename for input normalization parameters. Default: "inminmax.npy".
• outminmax_file::String: Filename for output normalization parameters. Default: "outminmax.npy".
• nn_setup_file::String: Filename for network architecture configuration. Default: "nn_setup.json".
• postprocessing_file::String: Filename for postprocessing function. Default: "postprocessing_file.jl".

Returns

A ComponentEmulator instance ready for evaluation.

Details

This function loads all necessary files to reconstruct a trained component emulator:

1. Neural network architecture from JSON configuration.
2. Trained weights from NumPy binary format.
3. Normalization parameters for inputs and outputs.
4. Wavenumber grid.
5. Postprocessing function dynamically loaded from Julia file.

The postprocessing function is evaluated in an isolated scope to prevent namespace pollution.

Example

P11_emu = load_component_emulator("/path/to/emulator/11/")

File Structure

The expected directory structure is:

path/
├── k.npy                    # Wavenumber grid
├── weights.npy              # Neural network weights
├── inminmax.npy            # Input normalization (n_params × 2)
├── outminmax.npy           # Output normalization (n_k × 2)
├── nn_setup.json           # Network architecture
└── postprocessing_file.jl  # Postprocessing function

load_multipole_emulator(path; emu=LuxEmulator, k_file="k.npy", weights_file="weights.npy", inminmax_file="inminmax.npy", outminmax_file="outminmax.npy", nn_setup_file="nn_setup.json", postprocessing_file="postprocessing.jl", stochmodel_file="stochmodel.jl", biascombination_file="biascombination.jl", jacbiascombination_file="jacbiascombination.jl")

Load a complete power spectrum multipole emulator from disk.

Arguments

• path: Directory path containing the multipole emulator structure.

Keyword Arguments

• emu: Emulator type to initialize (LuxEmulator or SimpleChainsEmulator). Default: LuxEmulator.
• k_file::String: Filename for the wavenumber grid. Default: "k.npy".
• weights_file::String: Filename for neural network weights. Default: "weights.npy".
• inminmax_file::String: Filename for input normalization parameters. Default: "inminmax.npy".
• outminmax_file::String: Filename for output normalization parameters. Default: "outminmax.npy".
• nn_setup_file::String: Filename for network architecture configuration. Default: "nn_setup.json".
• postprocessing_file::String: Filename for postprocessing function. Default: "postprocessing.jl".
• stochmodel_file::String: Filename for stochastic model function. Default: "stochmodel.jl".
• biascombination_file::String: Filename for bias combination function. Default: "biascombination.jl".
• jacbiascombination_file::String: Filename for bias Jacobian function. Default: "jacbiascombination.jl".

Returns

A PℓEmulator instance containing all three components (P11, Ploop, Pct) and bias models.

Details

This function loads a complete multipole emulator by:

1. Loading three component emulators (P11, Ploop, Pct) from subdirectories.
2. Loading the stochastic model function (shot noise terms).
3. Loading the bias combination function (maps bias parameters to weights).
4. Loading the analytical Jacobian of the bias combination.

All functions are evaluated in isolated scopes to prevent namespace conflicts between different emulator components.

Example

# Load monopole emulator
monopole_emu = load_multipole_emulator("/path/to/artifact/0/")

# Evaluate
cosmology = [z, ln10As, ns, H0, ωb, ωcdm, mν, w0, wa]
bias = [b1, b2, b3, b4, b5, b6, b7, f, cϵ0, cϵ1, cϵ2]
D = 0.8
P0 = get_Pℓ(cosmology, D, bias, monopole_emu)

File Structure

The expected directory structure is:

path/
├── 11/                       # P11 component
│   ├── k.npy
│   ├── weights.npy
│   ├── inminmax.npy
│   ├── outminmax.npy
│   ├── nn_setup.json
│   └── postprocessing.jl
├── loop/                     # Ploop component
│   └── ... (same structure)
├── ct/                       # Pct component
│   └── ... (same structure)
├── stochmodel.jl            # Stochastic model function
├── biascombination.jl       # Bias combination function
└── jacbiascombination.jl    # Bias Jacobian function

See Also

• load_component_emulator: Load individual component emulators.
• get_Pℓ: Evaluate the loaded emulator.