Opal.jl

ALL, WARMUP, TRAIN, WASHOUT, SPREAD, DA
OBSALL, OBSWARMUP, OBSTRAIN, OBSWASHOUT, OBSSPREAD, OBSDA

Integer phase tags used to index a TimeStencils descriptor or a StencilArray. The OBS* variants select the same time window but at the observation cadence (dt_obs) rather than the model cadence (dt).

ℓ2(n::Int) -> AbstractMatrix

Returns the $n \times n$ Euclidean ($\ell^2$) distance matrix for a one-dimensional grid of n equally-spaced points. Also available as the alias euclidean.

struct AdaptiveFilter{F<:KalmanFilter} <: KalmanFilter

Wraps an inner Kalman filter F and estimates the process-noise covariance Q and observation-noise covariance R on-line using an innovation-based EM-like update.

At each step, the cached transition Jacobian F, observation Jacobian H, prior/posterior covariances, and innovations from the previous step are used to form new estimates $\hat Q$ and $\hat R$, which are then blended into the running estimates via an exponential moving average with weight step.

Fields:

• filter: the underlying KalmanFilter (can be any concrete subtype);
• last_posterior: copy of the posterior from the previous analysis step, needed for the adaptive update;
• step: EMA step size $\alpha \in (0,1]$ controlling how fast $Q$ and $R$ adapt;
• cache: pre-allocated workspace for Jacobians, covariances, and innovation history.

Construct via AdaptiveFilter(filter; step=0.1).

struct AddBias{A<:Number} <: BiasStyle

Appends a constant value as an extra component to the input vector, increasing its length by one.

struct AdjointProblem{A,B}

Encapsulates the pieces needed to identify unknown parameters – or, for an unparameterised StateMap (no ParamSpace attached), the initial condition – of a PDE- or ODE-constrained observation model via gradient-based optimisation with automatic differentiation.

The parameter-to-state map, observation operator, and state-to-loss map are precomposed into a single μ_obs_to_ℓ(μ,obs) -> Real closure at construction time (see build_loss); they are not stored as separate fields.

Fields:

• μ_to_u::A: the state map (an ODEStateMap or PDEStateMap); its own bounds and initial_condition (derived from its attached pspace, or nothing if unparameterised) supply the optimisation's box constraints and default initial guess – there is no separate pspace field here anymore;
• μ_obs_to_ℓ::B: the precomposed (μ,obs) -> ℓ loss closure.

Construct via AdjointProblem(μ_to_u, u_to_obs, args...), where args... is forwarded to build_loss (e.g. obs_noise for the StateToObservationMap case).

struct AlgebraicModel{T,A<:AbstractMatrix{T}} <: LinearModel
  matrix::A
end

Standard implementation of a LinearModel. The field matrix represents the constant Jacobian of the model itself.

struct BiasAwareFilter{A<:KalmanFilter} <: KalmanFilter

Wraps an inner Kalman filter A and corrects for systematic observation bias using a RecurrentNeuralNetwork trained online.

During a warm-up phase (first maxiter steps) the inner filter runs as-is, accumulating innovations to train the bias model. Once the counter exceeds maxiter (isaware returns true) the analysis step adjusts the innovation by the current bias estimate before computing the Kalman gain and state update.

Fields:

• filter: the underlying KalmanFilter (can be any concrete subtype);
• bias_model: a RecurrentNeuralNetwork that maps the current innovation to a bias vector;
• regularisation: penalty coefficient γ applied to the bias Jacobian in the modified gain;
• awareness: [IterCounter] tracking how many analysis steps have been performed;
• cache: pre-allocated workspace for Jacobians and intermediate matrices.

Construct via BiasAwareFilter(f, bias_model; γ=10, maxiter=0) where f is any KalmanFilter instance.

abstract type BiasStyle end

Base type for bias-appending strategies within a Modifier pipeline. Concrete subtypes: NoBias, AddBias.

struct BickelLevina <: TaperFunction end

Hard-thresholding taper: returns 1 for ‖x‖ ≤ 1 and 0 elsewhere.

struct Cai <: TaperFunction end

Piecewise-linear taper: 1 for |x| ≤ 1/2, linearly decaying to 0 at |x| = 1.

struct CalibratedFilter{A<:Filter,B<:Calibration} <: Filter

Wraps any Filter with a Calibration that corrects the observation error estimate at each assimilation step. All filter interface methods (get_prior, transition!, analyse!, etc.) delegate to the inner filter; only the observation update step additionally calls the calibration.

Construct via:

CalibratedFilter(filter, calibration)

struct CatArray{T,N,D} <: AbstractArray{T,N}

A block-structured array that stores its data in a flat Array{T,N} together with D pointer vectors (ptrs) that delimit the blocks along each blocked dimension.

Construct with vcatarray (blocks stacked along dimension 1) or dcatarray (blocks placed on the block diagonal of a matrix). Individual blocks are accessed via BlockArrays.Block indices or via nblocks.

struct ConstrainedLaw{A,B,N} <: Law{N}

Wraps a Law together with a box constraint (ConstrainTo) so that samples and resampled particles are projected back into the feasible set after each update step.

Construct via ConstrainTo(lb, ub) applied to a law (see also ConstrainTo):

prior = Particle(ConstrainTo(lb, ub), particles, weights)

abstract type Filter end

Type reserved for filter operators for data processes, such as the Kalman filter (KF). A filter should be given as input an transition operator, a observation operator, and a prior distribution, which should be representative of the distribution of the initial sample. The prior should be updated iteratively, assimilating observations at multiple different time instants, as well as the last bayesian estimate of the (unknown) state variable. Check the function loop for suitable iterative updates of the prior. Check here for more details on Kalman filtering. Subtypes:

• KalmanFilter

struct DAResults
  state_history::History
  obs_measures::ResultsTable
end

const DEnKF{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} = EnsembleKalmanFilter{A,B,C,D,E,F,DEnKFStrategy}

Implements the DEnKF. Construct by passing an Ensemble with strategy=DEnKFStrategy() as the prior.

struct DEnKFStrategy <: EnsembleStyle end

Trait for ensembles mimicking the DEnKF (deterministic EnKF) method:

• run the forecast step on each ensemble member (see forecast!);
• compute the Kalman gain $K$ as usual (see kalman_gain!);
• compute the ensemble innovations $ỹ$ (see innovation!);
• update the ensemble mean according to the formula:

\[μ ← μ + K ⋅ mean(ỹ),\]

where $μ$ is the ensemble mean;

• update the ensemble anomaly according to the formula:

\[A ← (I - rac{1}{2} K⋅H) A \]

where $H$ is the jacobian of the observation model, evaluated in the forecasted ensemble mean, and $A$ is the ensemble anomaly. This is the so-called deterministic approximation of DEnKF;

• update the ensemble $E$ according to the formula:

\[E[:,i] = A[:,i] + μ \]

for every $i = 1,...,nₑ$. The covariance $Σ$ is then estimated via its ensemble counterpart:

\[Σ = A ⋅ Aᵀ / (nₑ - 1)\]

abstract type DataAugmentation <: DataTransformation end

Base type for input-augmentation strategies applied to training data before feeding into the reservoir. Concrete subtypes: NoAugmentation, ScaledAugmentation.

abstract type DataRegularisation <: DataTransformation end

Base type for data-regularisation strategies (e.g. additive noise) applied to training inputs. Concrete subtypes: NoRegularisation, AdditiveNoiseRegularisation.

struct EchoStateNetwork <: RecurrentNeuralNetwork

Echo State Network (ESN) / reservoir computing model.

The hidden state evolves as

\[s_{t+1} = (1 - \alpha)\, s_t + \alpha \, f(\rho W s_t + \sigma W_{\mathrm{in}} x_t')\]

where $f$ is the activation function, $\rho$ is the spectral radius, $\sigma$ is the input scaling, and $\alpha$ is the leak rate. The readout $y_t = W_{\mathrm{out}} s_t'$ is linear and fitted by ridge regression via TrainRecurrentNeuralNetwork.

Fields:

• activation: element-wise activation (default tanh);
• state: mutable hidden-state vector $s_t$;
• initial_state: initial hidden-state vector $s_0$;
• weights: reservoir weight matrix (sparse);
• weights_in: input weight matrix;
• weights_out_T: transposed readout matrix (fitted during training);
• modifier_in, modifier_state: Modifier pipelines for inputs and states;
• radius: mutable spectral-radius reference $\rho$;
• scaling: mutable input-scaling reference $\sigma$;
• leak: leak rate $\alpha$.

Construct via EchoStateNetwork(ninput, nstate, noutput; kwargs...).

const EnKF{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} = EnsembleKalmanFilter{A,B,C,D,E,F,EnKFStrategy}

Implements the standard EnKF. Construct by passing an Ensemble with strategy=EnKFStrategy() (the default) as the prior.

struct EnKFStrategy <: EnsembleStyle end

Trait for ensembles mimicking the EnKF method:

• run the forecast step on each ensemble member (see forecast!);
• compute the Kalman gain $K$ as usual (see kalman_gain!);
• compute the ensemble innovations $ỹ$ (see innovation!);
• update the ensemble according to the formula:

\[E ← E + K ⋅ ỹ \]

where $E$ is the ensemble matrix. The mean $μ$ and covariance $Σ$ are then estimated via their ensemble counterparts:

\[μ = ∑ᵢ E[:,i] / nₑ Σ = ∑ᵢ (E[:,i] - μ)⋅(E[:,i] - μ)ᵀ / (nₑ - 1)\]

const EnSRKF{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} = EnsembleKalmanFilter{A,B,C,D,E,F,EnSRKFStrategy}

Ensemble Square-Root Kalman Filter. Deterministic (no perturbed observations), so the ensemble spread is updated via a square-root anomaly formula instead of the stochastic EnKF perturbation.

Construct by passing an Ensemble with strategy=EnSRKFStrategy() as the prior:

prior = Ensemble(values; strategy=EnSRKFStrategy())
f = KalmanFilter(transition, observation, prior)

struct EnSRKFStrategy <: EnsembleStyle end

Trait for ensembles using the Ensemble Square-Root Kalman Filter (EnSRKF) method. The filter is deterministic: observations are not perturbed, so the ensemble spread is updated exactly via a square-root factorisation of the innovation covariance rather than by the stochastic perturbation of the standard EnKF.

The analysis proceeds as follows:

• run the forecast step on each ensemble member (see forecast!);
• compute the observation anomaly $S = H A_f$ where $H$ is the Jacobian of the observation model and $A_f$ the forecast anomaly;
• assemble the innovation covariance

\[C = (n_e - 1)R + S S^{\top}\]

and eigendecompose it as $C = \Phi \Lambda \Phi^{\top}$;

• compute the mean Kalman gain

\[K = P_{xy} C^{-1}\]

where $P_{xy}$ is the state–observation cross-covariance;

• form $E = \Lambda^{-1/2} \Phi^{\top} S$ (shape $m \times n_e$) and compute its SVD $E = U \Sigma V^{\top}$;
• update the ensemble mean

\[\mu_a = \mu_f + K \tilde{y}\]

where $\tilde{y} = z - H \mu_f$ is the mean innovation;

• update the ensemble anomaly

\[A_a = A_f \cdot V \sqrt{I - \Sigma^{\top}\Sigma}\, V^{\top}\]

• rebuild each ensemble member as $E[:,i] = \mu_a + A_a[:,i]$.

struct Ensemble{C<:EnsembleStyle,A<:AbstractMatrix,B<:AbstractVector} <: SecondMoment
  values::A
  mean::B 
  anomaly::A
  strategy::C
end

This SecondMoment distribution represents the ensemble needed to run the EnsembleKalmanFilter. Fields:

• values: $n × n_e$-dimensional matrix storing the ensemble values;
• mean: $n$-dimensional vector representing the ensemble mean;
• anomaly: $n × n_e$-dimensional matrix storing the ensemble anomaly;
• strategy: trait of type EnsembleStyle which determines the update type of the ensemble.

struct EnsembleKalmanFilter{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law,G<:EnsembleStyle} <: KalmanFilter

Implements an Ensemble Kalman Filter. In particular:

• instead of propagating a single probability distribution as in a Kalman filter, we do so for several different (ensemble) distributions.
• the explicit update of the state covariance matrix is not required. Indeed, the variability of the state is implicitly encoded in the ensemble's spread.
• the remaining steps are equivalent to a standard Kalman Filter.

The EnsembleStyle type parameter G is inferred from the strategy of the prior Ensemble at construction time. Subtypes:

• EnKF
• DEnKF
• EnSRKF

abstract type EnsembleStyle end

Trait specifying how the ensemble covariance of an Ensemble distribution should be updated. The reason why this is kept as a parameter is that, in ensemble filtering strategies, computing the ensemble covariance with the usual formula

\[Σ = ∑ᵢ (E[:,i] - μ)*(E[:,i] - μ)ᵀ / (nₑ - 1)\]

is generally expensive, and thus alternative strategies are sought. Here, $E$ are the ensemble members. Subtypes:

• EnKFStrategy
• DEnKFStrategy

const FirstMoment = Law{1}

Type reserved for distributions characterised only by their first moment, i.e. the mean, accessed via the function mean.

struct GaspariCohn <: TaperFunction end

Fifth-order piecewise-polynomial taper (Gaspari & Cohn 1999). Compactly supported on |x| ≤ 1, smooth at the origin, and non-negative everywhere. The standard default for ensemble covariance localisation.

struct GaussianTaper <: TaperFunction end

Gaussian taper exp(-x²/2). Globally supported (never exactly zero), so it damps but does not hard-threshold long-range correlations.

struct GenericFirstMoment{A<:AbstractVector} <: FirstMoment
  mean::A 
end

Most basic implementation of a FirstMoment distribution.

struct GenericModel{F<:FType} <: NonlinearModel
  form::F
end

Standard implementation of a NonlinearModel. The field form represents the function or Gridap Map characterising the model itself.

struct ImportanceSampling <: ResamplingStyle

Resampling strategy for Particle filters that uses systematic importance resampling without kernel smoothing (SIR filter). Particles are resampled when the effective sample size drops below the threshold stored in the Particle prior.

struct InflationFilter{A<:KalmanFilter,B<:InflationModel} <: KalmanFilter

Wraps an inner Kalman filter A and inflates the forecast covariance before computing the Kalman gain, using the inflation strategy B.

Covariance inflation counteracts ensemble collapse and filter divergence by artificially broadening the prior. Two strategies are available:

• MultInflation: multiplies the covariance by a fixed scalar ρ;
• NLLInflation: adaptively tunes ρ at each step by minimising the negative log-likelihood of the innovation.

Fields:

• filter: the underlying KalmanFilter;
• inflation: an InflationModel that supplies (and optionally optimises) the inflation factor;
• cache: stashed prior/obs-prior and optimisation workspace for the adaptive case.

Construct via InflationFilter(filter; inflation=MultInflation()) or from raw transition/observation/prior arguments (which automatically includes localisation).

abstract type InflationModel end

Base type for covariance-inflation strategies used by InflationFilter. An InflationModel supplies a scalar factor ρ by which the forecast covariance is multiplied before the Kalman gain is computed.

Subtypes:

• MultInflation: constant multiplicative inflation;
• NLLInflation: adaptive inflation tuned by negative log-likelihood.

struct KalmanFilter{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} <: Filter
  transition::A
  observation::B
  prior::C
  obs_prior::D
  noise::E
  obs_noise::F
  cache::KalmanCache
end

Kalman filter subtype implementing a Kalman filter procedure. Fields:

• transition: Model representing the transition operator;
• observation: Model representing the observation operator;
• prior: Law representing the probability distribution for the state;
• obs_prior: Law representing the probability distribution for the observation;
• noise: Law representing the probability distribution for the process (state) noise;
• obs_noise: Law representing the probability distribution for the observation noise;
• cache: cached object allowing for efficient in-place operations.

struct KrigingCalibration <: UnbiasedCalibration

Calibration strategy based on Kriging that estimates a bias correction and inflated observation-error covariance from reduced-basis snapshot data.

At each assimilation step the Kriging system is solved using the precomputed variogram fit to produce a weight vector λ; the correction is then applied to the current observation error estimate.

Construct via:

KrigingCalibration(observation, fesnaps, rbsnaps; variogram=ParametricSphere{2}())
KrigingCalibration(observation, fesnaps, rbsnaps, ts)  # restrict to DA phase of ts

abstract type Law{N} end

Type representing a probability distribution characterised by N moments. Subtypes:

• FirstMoment
• SecondMoment

struct Linear <: Linearity end

Trait used by LinearModel.

const LinearModel = Model{Linear}

Models that are fully characterised by an $m × n$-dimensional Jacobian matrix $J$, i.e.

\[x ↦ J⋅x\]

represents the action of a LinearModel on an $n$-dimensional vector $x$. If the input is a distribution $d$, then:

• if $d$ is a FirstMoment distribution with mean $μ$, then the output is a FirstMoment

with mean $J⋅μ$;

• $d$ is a SecondMoment distribution with mean $μ$ and covariance $Σ$, then the output

is a SecondMoment with mean $J⋅μ$, and covariance $J⋅P⋅Jᵀ$.

abstract type Linearity end

A Model trait that allows certain optimizations during the evaluation of the models themselves.

struct LocalisationFilter{A<:KalmanFilter} <: KalmanFilter

Wraps an inner Kalman filter A and applies covariance localisation after each forecast step via a TaperModel.

Localisation suppresses spurious long-range correlations that arise in ensemble-based filters when the ensemble size is small relative to the state dimension. After the transition step the sample covariance (or anomaly matrix) is element-wise multiplied by a distance-decay taper function.

Fields:

• filter: the underlying KalmanFilter;
• taper: a TaperModel that defines the taper function and the distance matrix;
• cache: pre-allocated workspace for the taper computation.

Construct via LocalisationFilter(filter; taper=GaspariCohn(), npoints=n) or by passing transition, observation, and prior directly.

struct MemoryModel{A<:Linearity} <: Model{A}

A Model wrapper that pre-allocates and reuses the evaluation cache across calls.

Unlike a bare model — which allocates a fresh cache on every evaluate! call — a MemoryModel holds a persistent cache and prior that are reused. This makes it suitable for models (e.g. ODEModel) whose cache encapsulates integrator state that must be preserved between steps.

Fields:

• model: the wrapped Model;
• prior: the input distribution, updated in-place on each call;
• cache: the persistent evaluation cache returned by return_cache(model, prior).

Construct via MemoryModel(model, prior).

abstract type Model{A<:Linearity} <: Map end

Type used for operator-like quantities, such as functions or Gridap Maps. For performance reasons, we distinguish models depending on their Linearity trait. To evaluate a Model a in a point $x$ — usually an $n$-dimensional vector or a scalar — simply call

a(x)

which automatically triggers the function evaluate(a,x), originally defined in Gridap, and overwritten here a few times. For in-place evaluations, use instead the syntax

evaluate!(cache,a,x)

where cache = return_cache(a,x) is a suitable cached object.

The main characteristic of a Model is that it may also be evaluated in a probability distribution. Given an input Law prior, the output

posteriori = a(priori)

returns another distribution posteriori, which should be thought of the propagation of priori through the model a. The type of Model and input Law determine the expression of posterior.

struct Modifier{A<:NormaliseStyle,B<:TransformStyle,C<:BiasStyle} <: DataTransformation

Composable pre-processing pipeline applied to reservoir inputs or states. Sequentially applies normalisation → transformation → bias appending.

Construct via Modifier(; normalisation, transformation, bias) or use DoNotModify for the identity pipeline.

struct MultInflation <: InflationModel

Constant multiplicative covariance inflation. The forecast covariance is scaled by the fixed factor ρ at every analysis step.

Fields:

• ρ: inflation factor (default 1.0, i.e. no inflation).

Construct via MultInflation(ρ) or MultInflation() for the identity.

struct NLLInflation <: InflationModel

Adaptive multiplicative covariance inflation that tunes the factor ρ at every analysis step by minimising the negative log-likelihood of the innovation.

The optimal ρ is found by a scalar line search over [lower, upper]. The iteration stops when the reduction in NLL falls below tolerance.

Fields:

• bounds: (lower, upper) box constraint on ρ (defaults: (1e-3, 10.0));
• tolerance: stopping criterion for the inner NLL optimisation (default 1e-1);
• ρ: mutable reference to the current inflation factor (reset to -1.0 each step).

Construct via NLLInflation(; lower=1e-3, upper=10.0, tolerance=1e-1).

const NLLInflationKalmanFilter{A<:KalmanFilter} = InflationFilter{A,<:NLLInflation}

abstract type NeuralNetwork <: Map end

Base type for all neural-network models. Subtypes must implement get_parameters (returns the trainable weight matrices) and extend train / forecast. The main concrete subtype is RecurrentNeuralNetwork (and in particular EchoStateNetwork).

struct NoAugmentation <: DataAugmentation end

Identity augmentation: returns the input unchanged.

struct NoBias <: BiasStyle end

Identity bias: returns the input unchanged (no extra component appended).

struct NoNormalisation <: NormaliseStyle end

Identity normalisation: returns the input unchanged.

struct NoRegularisation <: DataRegularisation end

Identity regularisation: returns the input unchanged.

struct NoTransformation <: TransformStyle end

Identity transformation: returns the input unchanged.

struct Nonlinear <: Linearity end

Trait used by NonlinearModel.

const NonlinearModel = Model{Nonlinear}

Models that are in general characterised by either a function, or a Gridap Map. Denoting such function/Map by by f, the action of a NonlinearModel on an $n$-dimensional vector $x$ is simply defined by

\[x ↦ f(x)\]

where the output is an $m$-dimensional vector, or a scalar. If the input is a distribution $d$, then:

• if $d$ is a FirstMoment distribution with mean $μ$, then the output is a FirstMoment

with mean f(μ);

• $d$ is a SecondMoment distribution with mean $μ$ and covariance $Σ$, then the output

is a SecondMoment with mean f(μ), and covariance whose definition depends on the types of boht a and $d$.

struct NormalLaw{A<:AbstractVector,B<:AbstractMatrix} <: SecondMoment
  mean::A 
  covariance::B
end

Normal distribution of mean mean and covariance covariance; a SecondMoment distribution defaults to a NormalLaw.

struct Normalisation{A<:AbstractVector} <: NormaliseStyle

Element-wise normalisation: divides each component x[i] by factor[i].

abstract type NormaliseStyle end

Base type for normalisation strategies within a Modifier pipeline. Concrete subtypes: NoNormalisation, Normalisation.

const ODEAdjointProblem = AdjointProblem{<:ODEStateMap}

An AdjointProblem built from an ODE (ODEWrapper) state map; dispatches to the Optimization.jl/Fminbox(BFGS()) optimise method.

struct ODEModel{A} <: DifferentialModel

Wraps an ODEWrapper (an ODE integrator) as a NonlinearModel.

Each call to evaluate! advances the integrator by one step from the current state, making it suitable as the transition model in a filter operating on ODE-governed systems. Construct via Model(sol::ODEWrapper).

const PDEAdjointProblem = AdjointProblem{<:PDEStateMap}

An AdjointProblem built from a PDE (AffineFEStateMap/ NonlinearFEStateMap) state map; dispatches to the Optim.jl/Fminbox(LBFGS()) optimise method.

struct ParametricSphere{N} <: VariogramModel

Parametric spherical variogram model with N free parameters, used inside KrigingCalibration to fit the spatial covariance structure of the observation errors from snapshot data.

• ParametricSphere{1}(): unit-sill spherical model (one range parameter);
• ParametricSphere{2}(): two-parameter model (sill + range).

struct Particle{A,B,C,D} <: FirstMoment

Weighted particle cloud representing a probability distribution for use with ParticleFilter.

Fields:

• particles: $n \times N_s$ matrix; each column is one particle;
• weights: $N_s$-vector of normalised importance weights;
• mean: $n$-vector of the particle mean;
• strategy: ResamplingStrategy that controls when and how resampling occurs.

Construct via:

Particle(particles, weights, ImportanceSampling())  # SIR
Particle(particles, weights, RegularisedSampling()) # RPF (default)
Particle(particles, weights)                        # same as RegularisedSampling

Pass nthreshold as a keyword to override the effective-sample-size threshold.

struct ParticleFilter{...} <: KalmanFilter

Sequential Monte Carlo filter using a particles approximation of the state distribution. The resampling strategy is determined by the ResamplingStrategy embedded in the Particle prior; by default this is RegularisedSampling, which runs the Regularised Particle Filter (RPF) of Berry & Sauer (2002).

Construct via the generic interface:

KalmanFilter(transition, observation, prior::Particle; Q, R, noise, obs_noise)

Fields:

• transition: state transition Model;
• observation: observation Model;
• prior: Particle prior distribution;
• obs_prior: observation-space particles distribution;
• noise: process noise law (default: zero-mean NormalLaw);
• obs_noise: observation noise law.

struct RTS <: Smoother

Rauch–Tung–Striebel (RTS) smoother.

Given the sequence of filter posteriors produced by loop, runs a backward pass that refines each estimate using information from all future observations. For a linear-Gaussian model the RTS smoother is the exact fixed-interval smoother.

Use via smooth_loop to obtain smoothed histories directly, or call smoothen! on an existing filter history.

abstract type RecurrentNeuralNetwork <: NeuralNetwork end

Base type for recurrent neural networks. Subtypes must implement get_state (internal hidden state) and get_output (current readout). The main concrete subtype is EchoStateNetwork.

struct RecycleValidation{A<:TrainMethod,B<:UpdateRule} <: TrainMethod

Training method that wraps another TrainMethod and tunes hyper-parameters via cross-validation on held-out forecast windows (the "recycle" trick).

For each candidate in updates the network is trained, then its closed-loop forecast error over each window in windows is evaluated with loss. The candidate with the lowest total loss is selected, and a final training pass is performed with it.

Fields:

• method: inner training method (e.g. TrainRecurrentNeuralNetwork);
• updates: an UpdateRule iterable of hyper-parameter candidates;
• windows: tuple of index ranges defining the validation forecast windows;
• loss: scoring function (true, predicted) -> Real (default log10RMSE).
• iterations: number of iterations for local refinement step after the grid search.

Construct via RecycleValidation(method, ranges...; Nfolds, Ntrain, Nvalidation, loss, iterations).

struct RegularisedParticleMetadata <: ParticleMetadata

Pre-allocated scratch workspace for the Regularised Particle Filter resampling step.

Fields:

• particles_scratch: n_x × n_s buffer used for systematic resampling;
• Sk: n_x × n_x buffer for the weighted empirical state covariance;
• Dk: n_x × n_x buffer for the lower Cholesky factor of Sk;
• c: n_x buffer for sampling from the Epanechnikov kernel.

struct RegularisedSampling <: ResamplingStyle

Resampling strategy for Particle filters that applies the optimal-bandwidth Epanechnikov kernel perturbation after systematic resampling (Regularised Particle Filter, RPF). This is the default resampling style and prevents particle collapse in low-noise settings while preserving the correct second moment.

struct RidgeRegression <: GridapType

Tikhonov-regularised (ridge) least-squares solver.

Solves $\min_X \|AX - B\|^2 + \lambda \|X\|^2$ via the normal equations $(A A^\top + \lambda I) X = A B^\top$. The regularisation parameter λ is stored as a mutable reference so it can be updated by RecycleValidation during hyper-parameter search.

Construct via RidgeRegression(λ::Real).

const SecondMoment = Law{2}

Type reserved for distributions characterised by their first two moments, i.e. mean and covariance, accessed via the functions mean and cov.

mutable struct SecondOrderResultsTable <: ResultsTable
  innovation_means::AbstractVector{<:AbstractVector{<:Real}}
  innovation_stds::AbstractVector{<:AbstractVector{<:Real}}
  innovation_nis::AbstractVector{<:Real}
  innovation_rmse::AbstractVector{<:Real}
end

Stores innovation diagnostics collected step-by-step during loop. Access via result.table on the returned DAResults.

Fields:

• innovation_means: mean innovation vector at each step (length-m vectors)
• innovation_stds: diagonal std of the observation prior at each step
• innovation_nis: Normalized Innovation Squared at each step
• innovation_rmse: RMS of the mean innovation at each step

struct SigmaPoints{A<:AbstractVector,B<:AbstractMatrix,C<:AbstractMatrix,D<:AbstractVector,E<:Real} <: SecondMoment
  mean::A 
  covariance::B
  points::C 
  weights_mean::D 
  weights_cov::D
  λ::E
end

This SecondMoment distribution represents the sigma points needed to run the UnscentedKalmanFilter. Fields:

• points: $n × (2*L + 1)$-dimensional matrix storing the values of the sigma points;
• mean: $n$-dimensional vector representing the weighted mean of points;
• covariance: $n × n$-dimensional matrix representing the weighted covariance of points;
• weights_mean: $(2*L + 1)$-dimensional vector storing weights for mean;
• weights_cov: $(2*L + 1)$-dimensional vector storing weights for covariance;
• λ: real value used in the update of the sigma points points.

In an Unscented Transformation, two things occur in an iterative fashion:

• the field points is updated in-place via a call to sigma_points!;
• the fields mean and covariance are updated in-place via a call to update!.

struct StateToObservationMap{A<:Linearity} <: Model{A}
  model::Model{A}
  u_to_obs_ids::AbstractVector
end

Wraps a Model model (typically a LinearModel selection/observation operator) and scatters its output back into a state-sized vector, placing model(u)[k] at u_to_obs_ids[k] and zero everywhere else.

Two evaluation modes are supported:

• evaluate(a,u): the raw scatter of model(u);
• evaluate(a,u,obs,W): the scatter of W*(model(u) - obs), i.e. a weighted residual against observations obs – this is the form build_loss uses to assemble the parameter-identification loss.

Custom rrules are defined for both modes below; both linearise model (via jac) at the point of evaluation, consistent with how build_loss always uses it. u_to_obs_ids defaults to _find_u_to_obs_ids when omitted, which requires model::LinearModel and errors if it is not injective.

struct StencilArray{A<:AbstractArray,B}

Wraps an array together with a TimeStencils descriptor and a phase tag so that it can be re-indexed at a different phase (with automatic restriction or expansion).

Construct via to_stencil; index with a phase constant to extract or rephase the data:

sa = to_stencil(x, ts, TRAIN)
sa[DA]    # restrict (or expand) x to the DA phase

struct TaperModel <: NonlinearModel

Applies element-wise covariance localisation to a sample covariance or ensemble anomaly.

Given a covariance matrix $\Sigma$, the localised matrix is

\[\Sigma_{ij}^{\text{loc}} = f\!\left(d_{ij}/r\right) \cdot \Sigma_{ij}\]

where $f$ is a TaperFunction, $d_{ij}$ is the pre-computed distance between indices $i$ and $j$, and $r$ is the localisation radius (mutable, may be optimised).

The result is projected back to the nearest positive-semidefinite matrix via SVD to ensure it remains a valid covariance.

Fields:

• taper: a TaperFunction instance;
• distance: symmetric matrix of pairwise distances;
• radius: mutable reference to the localisation radius.

Construct via TaperModel(npoints; taper=GaspariCohn(), distance=ℓ2, radius=1.0).

struct TimeStencils

Descriptor of the time axis partitioned into named phases used by the RC/DA pipeline.

Fields:

• dt: model integration time step;
• dt_obs: observation time step (defaults to dt);
• all_window, warmup_window, train_window, washout_window, spread_window, da_window: (t_start, t_end) tuples for each phase.

Construct via keyword arguments:

ts = TimeStencils(dt=0.01, dt_obs=0.1, t_warmup=10.0, t_train=50.0, t_wash=5.0, t_spread=5.0, t_da=20.0)

Index with phase constants to get the corresponding time stencil:

ts[DA]      # time steps in the DA phase
ts[TRAIN]   # time steps in the training phase

struct TrainRecurrentNeuralNetwork <: TrainMethod

Training strategy for RecurrentNeuralNetwork subtypes.

The readout weights are fitted by ridge regression on the reservoir states collected during an open-loop run. Optional data augmentation and regularisation are applied before fitting, and a forget period discards the initial transient.

Fields:

• solver: a RidgeRegression solver (holds the Tikhonov parameter λ);
• augmentation: a DataAugmentation applied to inputs before training;
• regularisation: a DataRegularisation applied to reservoir states;
• forget: number of initial steps to discard before regression.

Construct via TrainRecurrentNeuralNetwork(; augmentation=..., regularisation=..., forget=0, λ=1e-16).

abstract type TransformStyle end

Base type for nonlinear input-transformation strategies within a Modifier. Concrete subtypes: NoTransformation, T₁, T₂, T₃.

struct TransientPDEModel{A<:ODESolution} <: DifferentialModel

Wraps a Gridap/GridapROMs ODESolution (a parametric or non-parametric PDE solution object) as a NonlinearModel.

Each call to evaluate! advances the PDE solver by one time step, making it suitable as the transition model in a filter for PDE-governed dynamical systems. Construct via Model(sol::ODESolution).

TransientPDEStateMap(step,U,V,P,u0,grid,pspace,output_ids=eachindex(grid); p=nothing, kwargs...)

Build a PDEStateMap for a transient PDE. p is only used when pspace is nothing (an unparameterised map): it is the fixed physical parameter combined with the varying u0 at each step (see evaluate) – the PDE analogue of ODEWrapper's prob.p in the ODEStateMap{Nothing} case. When pspace is an actual ParamSpace, p is unused (the physical parameter is the thing being varied instead) and should be left at its default.

struct T₁ <: TransformStyle end

Transformation that squares odd-indexed components: x[i] = i odd ? x[i]^2 : x[i].

struct T₂ <: TransformStyle end

Transformation that replaces odd-indexed components (after the first) with the product of the two preceding components: x[i] = i > 1 && i odd ? x[i-1]*x[i-2] : x[i].

struct T₃ <: TransformStyle end

Transformation that replaces interior odd-indexed components with the product of their neighbours: x[i] = 1 < i < n && i odd ? x[i+1]*x[i-1] : x[i].

struct UniformLaw{A<:AbstractVector,B<:AbstractMatrix} <: SecondMoment
  mean::A 
  covariance::B
  lower_bound::A
  upper_bound::A
end

Uniform distribution of mean mean and covariance covariance.

struct UnscentedKalmanFilter{A<:Model,B<:Model,C<:Law,D<:Law,E<:Law,F<:Law} <: KalmanFilter

Implements the unscented transform. This model is analogous to a Kalman filter procedure for nonlinear transition/observation models, handling nonlinearities via sigma points (see SigmaPoints). These are a set of interpolation points used to approximate the first and second moments of a transitioned/observed distribution via a weighted linear combination. In particular:

• for the forecasting step: we update the sigma points from the analysed distribution of the previous iteration; evaluate the transition model at each sigma point; update mean and covariance via weighted linear combinations (see SigmaPoints);
• for the analysis step: we evaluate the observation model at each sigma point; update mean and covariance; compute innovation, Kalman gain, and state posterior as in a standard Kalman filter.

The prior C and observation prior D must both be (or wrap) SigmaPoints. Construct via KalmanFilter(transition, observation, prior) with a SigmaPoints (or ConstrainedSigmaPoints) prior.

abstract type UpdateRule end

Base type for hyper-parameter search strategies used by RecycleValidation. Concrete subtypes iterate over candidate parameter sets:

• UpdateRule(stencils...): grid-search over Cartesian product of ranges (NetworkUpdate);
• UpdateRule(tikhonov, ranges...): joint grid over ridge λ and network parameters (NetworkAndTikhonovUpdate).

struct VariationalMethod{A<:StateMap,B} <: DAMethod

Variational data-assimilation method (4D-Var / weak-constraint 4D-Var style).

Minimises a cost function that balances fidelity to observations (obs_to_ℓ) against departure from the background (back_noise) over a time window defined by μ_to_u.

Fields:

• μ_to_u: StateMap from control parameters to model state trajectories;
• u_to_obs: StateToObservationMap projecting the trajectory to observation space;
• obs_to_ℓ: loss function (ỹ, μ) -> ℓ measuring observation misfit;
• back_noise, obs_noise: background and observation error SecondMoment distributions.

Construct via:

VariationalMethod(μ_to_u, u_to_obs, obs_to_ℓ; B=..., R=...)

DoNotModify() -> Modifier

Returns the identity Modifier: no normalisation, no transformation, no bias.

InnovationACF(t::ResultsTable;maxlag=20) -> AbstractVector

Autocorrelation function of the innovation-norm time series from table. Under a well-specified filter innovations should be white noise, so ACF[lag > 0] ≈ 0.

NEES(true_values::AbstractVector,d::Law) -> Real

Normalized Estimation Error Squared. Under a well-calibrated filter, the expected value is 1. Values > 1 indicate underestimated uncertainty; < 1 indicate overestimated uncertainty.

NEES(true_values::AbstractMatrix,history::History) -> AbstractVector

NIS(true_values::AbstractVector,d::Law) -> Real

Normalized Innovation Squared. Under a consistent filter, the expected value is 1.

NIS(true_values::AbstractMatrix,history::History) -> AbstractVector

NLL(true_values::AbstractVector,d::Law) -> Real

Computes the Negative Log Likelihood between the true data true_values and the first moment of the distribution d.

NLL(true_values::AbstractMatrix,history::History) -> AbstractVector

Computes the Negative Log Likelihood between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

NRMSE(true_values::AbstractVector,d::Law) -> Real

Computes the Normalised Root Mean Square Error between the true data true_values and the first moment of the distribution d. This is equal to the Root Mean Square Error divided by the standard deviation of d.

NRMSE(true_values::AbstractMatrix,history::History) -> AbstractVector

Computes the Normalised Root Mean Square Error between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

Noise(Σ::AbstractMatrix) -> NormalLaw

Constructs a zero-mean NormalLaw with covariance Σ. Intended for specifying process noise $\theta$ and observation noise $\eta$ in a filter.

NovoaEchoStateNetwork(args...; kwargs...) -> EchoStateNetwork

Variant of EchoStateNetwork using the sparse reservoir initialisation from Novoa et al., with one non-zero entry per row in weights_in. Accepts the same arguments as EchoStateNetwork.

RMSE(true_values::AbstractVector,d::Law) -> Real

Computes the Root Mean Square Error between the true data true_values and the first moment of the distribution d.

RMSE(true_values::AbstractMatrix,history::History) -> AbstractVector

Computes the Root Mean Square Error between the true data true_values and history, the historical distributions obtained by running the Kalman iterations.

RankHistogram(true_values::AbstractMatrix,history::AbstractVector{<:Ensemble}) -> AbstractVector

Rank (Talagrand) histogram. For each time step and state variable, counts the rank of the true value within the sorted ensemble. A flat histogram indicates a well-calibrated ensemble; U-shaped indicates underdispersion; dome-shaped indicates overdispersion. Returns a normalised frequency vector of length ne + 1.

SpreadSkillRatio(true_values::AbstractVector,d::Law) -> Real

Ratio of mean ensemble spread to RMSE. Values ≈ 1 indicate a well-calibrated ensemble. Values < 1 indicate underdispersion; > 1 indicate overdispersion.

SpreadSkillRatio(true_values::AbstractMatrix,history::History) -> AbstractVector

ad_compatible(a) -> Function

Wraps a coefficient-building function a, where a(μ) returns a spatial closure x -> ... for a numeric parameter vector μ, so that the result also works when called with a parameter CellField μh (e.g. an FEFunction on a ConstantFESpace, or a dual-valued CellField produced internally by Gridap's automatic differentiation w.r.t. μh's DOFs). The returned function aμ(μh) builds a proper lazy CellField via Operation, evaluated at the true quadrature points of whichever Measure it ends up composed into – it never extracts μh into a raw numeric vector, which would break under AD (compare to test/AD.jl's equivalent, hand-written Operation-based pattern).

Example

a(μ) = x -> 1 + exp(-x[1]/sum(μ))
aμ = ad_compatible(a)
afe(u,v,μh) = ∫(aμ(μh)*∇(v)⋅∇(u))dΩ

analyse!(posterior::Law,f::Filter,args...) -> Law

In-place execution of the analysis step of a Kalman filter algorithm. This step consists of the following operations:

• Running the observation model on the forecasted distribution (see observation!)
• Computing the innovation (see innovation!)
• Computing the Kalman gain (see kalman_gain!)
• Updating the forecasted distribution by accounting for the Kalman gain (see update!)

To run an iteration of the Kalman filter, one must run the forecasting step in forecast! prior to the analysis one. If no optional argument is provided, the analysis is not performed.

build_linear_observation_model(
  ids::AbstractVector,
  obs_ids::AbstractVector = ids
) -> AlgebraicModel

Builds a linear observation model (selection matrix $H$) that extracts the components at indices obs_ids from a state vector of length length(ids).

obs_ids must already be indices into the full state vector (offset by the caller if the observed quantity is a sub-block of a larger vector). Returns an AlgebraicModel wrapping the sparse (length(obs_ids) × length(ids)) matrix.

build_observations(
  a::Model,
  x::AbstractMatrix,
  [obs_noise::Law],
  args::Function...
) -> AbstractMatrix

Applies the observation model a column-by-column to the state matrix x (each column is one time step) and collects the results into an (m × T) observation matrix, where m = dimension(a(x[:,1])) and T = size(x, 2).

Optional positional Function arguments b are called as b(x[:,k]) and their output is added to the observation at step k (useful for adding a bias term). If obs_noise::Law is provided, a random draw from it is added to every column.

build_prior(state::AbstractVector; kwargs...) -> FirstMoment
build_prior(states::AbstractMatrix; nsamples=1, kwargs...) -> Ensemble
build_prior(state, noise::SecondMoment; nsamples=1, kwargs...) -> SecondMoment or Ensemble
build_prior(states, c::AbstractConstraint; kwargs...) -> ConstrainedLaw

Constructs the initial prior distribution from a state snapshot or ensemble matrix.

• A Vector with no noise returns a FirstMoment;
• a Matrix returns an Ensemble;
• adding a noise::SecondMoment draws nsamples perturbations and sets the covariance;
• adding a constraint c wraps the result in a ConstrainedLaw.

AbstractBlockVector/AbstractBlockMatrix inputs produce block-structured joint distributions.

dimension(d::Law) -> Int

Dimension of the (vector) space on which the distribution is defined.

draw(d::SecondMoment) -> AbstractVector 
draw(d::SecondMoment,nsamples::Int) -> AbstractMatrix

Draws a $n$-dimensional random vector from the distribution d, where $n$ represents the dimension of d (see distribution). If an integer nsamples is also provided, the output will be an $n × nsamples$ - dimensional matrix. γ represents a scaling factor for the draw.

equispaced_windows(nobs, nwindows=1) -> NTuple of UnitRange
equispaced_windows(grid, nwindows=1)
equispaced_windows(ts::TimeStencils, nwindows=1)

Partition nobs observation indices (or the DA phase of ts) into nwindows nearly equal contiguous windows. Returns a tuple of UnitRanges.

execute(a::Model, prior::Law, stencil::AbstractVector) -> History
execute(f::Filter, stencil::AbstractVector) -> History

Runs the model a (or the transition model of filter f) forward for length(stencil) steps starting from prior, recording a copy of the distribution at each step. Returns the history as a vector of distributions.

forecast!(posterior::Law,f::Filter) -> Law

In-place execution of the forecast step of a Kalman filter algorithm. This step consists of the following operations:

• Running the transition model on the prior distribution (see transition!)

To complete a single iteration of the Kalman filter, one must run the analysis step in analyse! following the forecast one.

forecast(network::NeuralNetwork, stencil; kwargs...) -> AbstractArray

Runs network in closed-loop (autoregressive) mode for length(stencil) steps: the output at each step is fed back as the input for the next. Returns the sequence of outputs.

forecasted_history(args...) -> History

Alias for execute. Runs the model forward and returns the full history of distributions. Accepts the same arguments as execute.

forecasted_law(args...) -> Law

Returns only the final distribution from forecasted_history.

geostrophic(n::Int) -> AbstractMatrix
geostrophic(n::Dims{2}) -> AbstractMatrix

Returns the pairwise Euclidean distance matrix for a 2-D Cartesian grid of shape (n, n) (or n when given as a Dims{2} tuple), laid out in column-major order. Intended for geophysical applications where the state is defined on a regular 2-D domain.

get_noise(f::Filter) -> Law

Fetches the transition noise from the filter f. This noise, denoted by $θ$, is such that

\[xₙ₊₁ := F(xₙ) + θ\]

where ${xₖ}ₖ$ is the state process, and $F$ is the transition model.

get_observation_model(f::Filter) -> Model

Fetches the observation model from the filter f. This model, denoted by $H$, is such that

\[yₙ := H(xₙ) + η\]

where ${xₖ}ₖ$ is the state process, ${yₖ}ₖ$ is the observed process, and $η$ is a (usually Gaussian) random variable.

get_observation_noise(f::Filter) -> Law

Fetches the observation noise from the filter f. This noise, denoted by $η$, is such that

\[yₙ := H(xₙ) + η\]

where ${xₖ}ₖ$ is the state process, ${yₖ}ₖ$ is the observed process, and $H$ is the observation model.

get_observation_prior(f::Filter) -> Law

Fetches the distribution of the observation variable from the filter f.

get_prior(f::Filter) -> Law

Fetches the distribution of the state variable from the filter f.

get_transition_model(f::Filter) -> Model

Fetches the transition model from the filter f. This model, denoted by $F$, is such that

\[xₙ₊₁ := F(xₙ) + θ\]

where ${xₖ}ₖ$ is the state process, and $θ$ is a (usually Gaussian) random variable.

identify_parameter(
  ad::AdjointProblem,
  obs::AbstractMatrix,
  args...;
  x0 = initial_condition(ad.μ_to_u),
  iterations = 1000,
  show_trace = true,
  kwargs...
)

Identify the vector μ that best explains observations obs, for an AdjointProblem built from an ODE (ODEWrapper) state map, by minimising the loss precomposed into ad.μ_obs_to_ℓ (see build_loss). μ is the physical parameter for a ParamSpace-attached ODEStateMap, or the initial condition for an unparameterised one – see initial_condition.

Gradients are obtained via Zygote AD (Optimization.AutoZygote()). Optimisation uses Fminbox(BFGS()) (via OptimizationOptimJL) with box constraints from bounds(ad.μ_to_u) passed as lb/ub ((nothing,nothing) for an unparameterised map, i.e. unconstrained) – OptimizationPolyalgorithms.PolyOpt() was tried first, but despite checking prob.lb/prob.ub internally, OptimizationBase rejects it outright as an algorithm that "does not support box constraints", so it can't be used here.

Arguments

• obs: observed data, in the shape expected by ad.μ_obs_to_ℓ;
• x0: initial guess (defaults to initial_condition(ad.μ_to_u));
• iterations: forwarded to Optimization.solve as maxiters;
• show_trace: if true, the loss value is displayed at every callback;
• kwargs...: additional keyword arguments forwarded to Optimization.solve.

Returns the minimiser as a plain Array. Also reachable as optimise(ad::ODEAdjointProblem, obs, args...; p=initial_condition(ad.μ_to_u), kwargs...), which forwards to this function as x0=p – the p spelling is what VariationalMethod's loop/optimise machinery passes.

innovation!(f::Filter,posterior::Law,z::InType) -> InType

Given a distribution posterior and an observation z, returns the innovation $ỹ$ such that

\[ỹ = z - yₙ = z - [H(xᶠₙ) + η]\]

where $yₙ$ represents the observation forecasted by the filter f, and $xᶠₙ$ is distributed according to posterior.

jac(a::Union{Function,Map},x) -> AbstractMatrix

Returns a's Jacobian matrix evaluated in x.

joint_law(d::Law...) -> Law

Given a list of marginal distributions d = (d1,...,dn), returns their joint distribution.

kalman_gain!(f::Filter,posterior::Law) -> AbstractMatrix

In-place computation of the Kalman gain $K$ according to the formula

\[K = Σxy ⋅ Σy⁻¹\]

where $Σxy$ and $Σy$ are the state-observation and observation covariance matrices, respectively. $K$ is storede as a cached object in f, while the covariances $Σxy$ and $Σy$ can be computed by suitably accessing the transition and observation distributions via get_prior and get_observation_prior, respectively.

linearise(a::Model,x::InType) -> LinearModel

Linearizes a model a around $x$. If a is a LinearModel, it returns a itself.

loop(f::Filter,obs::AbstractArray) -> DAResults

Given a filter f and observations obs, iteratively runs the forecast-analyse paradigm typical of a Kalman filter, producing a DAResults that acts as a vector of posterior distributions and also carries a ResultsTable of innovation diagnostics.

If a slice of obs along its last dimension is all-NaN, the analysis step is skipped (forecast only) and NaN entries are pushed to the table. Use expand to embed sparse observations into a fine grid before passing to loop.

mixed_cov!(Σ::AbstractMatrix,f::Filter,posterior::Law) -> AbstractMatrix

In-place computation of the state-observation "mixed" covariance Σ.

mixed_cov!(cache,a::SigmaPoints,b::SigmaPoints) -> AbstractMatrix

In-place computation of the covariance between the SigmaPoints distributions a and b. These two distributions should have the same $L$ (i.e. the same number of sigma points) and $λ$ parameters, which also implies that they share the same mean/covariance weights. The formula used here is:

\[Σ = ∑ᵢ₌₁²ᴸ⁺¹ w[i] ⋅ (χᵃ[:,i] - μᵃ) ⋅ (χᵇ[:,i] - μᵇ)\]

where $χᵃ$ and $μᵃ$ are the sigma points and their mean for a, $χᵇ$ and $μᵇ$ are the sigma points and their mean for b, and $w$ are the covariance weights of either a or b.

observation!(f::Filter,posterior::Law) -> Law

In-place application of the observation model stored in f on the distribution posterior, which represents the posterior distribution of the state variable. In essence, denoting by $F$ the transition model, if posterior represents the forecasted distribution of the state variable xᶠₙ at the $n$th iteration, this step runs the observation model

\[yₙ := H(xᶠₙ) + η\]

overwriting the result yₙ in the distribution of the observed variable stored in f, and accessible through get_observation_prior. This function should be run during the analysis step in a Kalman filter algorithm.

observe(a::Model, d::Law) -> AbstractArray

Applies the model a to the state (mean or ensemble matrix) of the distribution d and returns the result without modifying d. Equivalent to a(get_state(d)).

optimise(
  ad::PDEAdjointProblem,
  obs::AbstractVecOrMat,
  args...;
  p = initial_condition(ad.μ_to_u),
  iterations = 1000,
  x_abstol = 1e-12,
  x_reltol = 1e-6,
  show_trace = true,
  kwargs...
) -> AbstractVector

Identify the vector μ that best explains observations obs, for an AdjointProblem built from a PDE state map, by minimising the loss precomposed into ad.μ_obs_to_ℓ (see build_loss) – typically the weighted least-squares residual

\[\ell(\mu) = \|R^{-1/2}(\mathcal{H}(u(\mu)) - y)\|_{\Omega}^2\]

where $u(\mu)$ solves the parametric PDE, $\mathcal{H}$ is the observation operator supplied at construction, and the norm is the finite-element $L^2(\Omega)$ norm. μ is the physical parameter for a ParamSpace-attached PDEStateMap, or the initial condition for an unparameterised one – see initial_condition.

Gradients are obtained via Zygote AD, which differentiates through the custom rrules on AffineFEStateMap, StateParamMap, and LinearModel. Optimization uses Fminbox(LBFGS()) with box constraints from bounds(ad.μ_to_u), falling back to unconstrained LBFGS() when those bounds are (nothing,nothing) (the unparameterised case).

Arguments

• obs: observed data in the observation space;
• p: initial guess (defaults to initial_condition(ad.μ_to_u));
• iterations, x_abstol, x_reltol, show_trace: forwarded to Optim.Options;
• kwargs...: additional keyword arguments forwarded to Optim.Options.

Returns the minimiser (Optim.minimizer(res)) as a plain vector.

predicted_history(args...) -> DAResults

Alias for loop. Runs the filter forward and returns the full assimilation history (posteriors + innovation table). Accepts the same arguments as loop.

predicted_law(args...) -> Law

Returns only the final posterior distribution from predicted_history.

sigma_points!(dcache::SigmaPoints,d::SigmaPoints;kwargs...) -> AbstractMatrix

In-place update of the sigma points according to the formula:

\[χₖ₊₁[:,1] = μₖ χₖ₊₁[:,2:L+1] = μₖ + √((L + λ)Pₖ) χₖ₊₁[:,L+2:2L+1] = μₖ - √((L + λ)Pₖ)\]

where $μₖ$ and $Pₖ$ are the previous mean and covariance fields. The output $χₖ₊₁$ overwrites the field points.

sigma_points(d::SecondMoment;kwargs...) -> AbstractMatrix

Given an input distribution d, computes the sigma points $χ$ according to the formula:

\[χ[:,1] = μ χ[:,2:L+1] = μ + √((L + λ)Σ) χ[:,L+2:2L+1] = μ - √((L + λ)Σ)\]

where μ and Σ are the mean and covariance of d, respectively. The variables $L$ and $λ$ may be passed as keyword arguments, and assume the following default values:

\[L = dimension(d) λ = 3 - L\]

smooth_loop(f::Filter, obs::AbstractArray; smoother=RTS()) -> DAResults

Combined forward-filter and backward-smoother pass.

Runs loop on f with observations obs, then immediately applies smoothen! to refine the history. Returns a DAResults whose state_history contains the smoothed posteriors.

smoothen!(
  smooth_history::History,
  smoother::Smoother,
  filter::Filter,
  history::History
)

In-place backward smoothing pass. Given the forward-filter posteriors history (from loop), overwrites smooth_history with the smoothed distributions produced by smoother.

The two-argument form smoothen!(smooth_history, filter, history; smoother=RTS()) also works and defaults to the RTS smoother.

state_blocks(d::Law) -> Tuple

Splits the state vector of d into its constituent blocks, as returned by BlockArrays.blocks. Used in variational methods to separate control parameters from physical-state components.

transition!(posterior::Law,f::Filter) -> Law

In-place application of the transition model stored in f on the distribution posterior, which represents the posterior distribution of the state variable. In essence, denoting by $F$ the transition model, if posterior represents the distribution of the state variable $xₙ$ at the $n$th iteration, this step runs the transition model

\[xₙ₊₁ := F(xₙ) + θ\]

overwriting the result $xₙ₊₁$ in posterior. This function should be run during the forecast step in a Kalman filter algorithm.

update!(cache,d::SigmaPoints) -> SigmaPoints

Update the mean and covariance of the sigma points d.

update!(d::Ensemble) -> Ensemble

Update the mean and anomaly of the ensemble d.

update!(posterior::Law,f::Filter) -> Law

In-place update of the distribution posterior through the action of Kalman gain matrix cached in f. Denoting by $K$ the Kalman gain computed by running kalman_gain!, and by $ỹ$ the innovation computed by running innovation!, if posterior represents the forecasted distribution of the state variable $xᶠₙ$ at the $n$th iteration, this step runs the formula

\[xᵃₙ := xᶠₙ + K ⋅ ỹ\]

and overwrites the analysed distribution of the state variable $xᵃₙ$ in posterior.

update!(table::ResultsTable,f::Filter)

Extracts the current innovation and observation-prior std from f and appends NIS, RMSE, mean and std to table. Called automatically inside loop.

warmup!(f::Filter, stencil::AbstractVector)
warmup!(a::MemoryModel, prior::Law, stencil::AbstractVector)

Advances the model or filter forward for length(stencil) steps, updating the internal state in-place without recording the history. Used to spin up memory-based models (e.g. EchoStateNetwork, ODEModel) before the assimilation window.

ℓ1(n::Int) -> AbstractMatrix

Returns the $n \times n$ periodic $\ell^1$ (Manhattan) distance matrix for a one-dimensional grid of n equally-spaced points, with wrap-around at the boundary.

train(method::TrainMethod, network::NeuralNetwork, x, y; kwargs...) -> states

Fits network on training data (x, y) using method. Allocates the training cache, runs the open-loop pass, and solves for the readout weights in-place. Returns the collected reservoir states.