MeteoModels.jl

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.

const DEnKF{A<:Model,B<:Model} = EnsembleKalmanFilter{A,B,<:Ensemble{DEnKFUpdate},<:Ensemble}

Implements the DEnKF. Simply requires a specialization of the update! function.

struct DEnKFUpdate <: NonstandardCovUpdate 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 + 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ₑ$.

abstract type Distribution{M} end

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

• FirstMoment
• SecondMoment

const EnKF{A<:Model,B<:Model} = EnsembleKalmanFilter{A,B,<:Ensemble{EnKFUpdate},<:Ensemble}

Implements the standard EnKF. Simply requires a specialization of the update! function.

struct EnKFUpdate <: NonstandardCovUpdate 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:

\[ensemble = ensemble + K ⋅ ỹ + θ\]

where $θ$ is an $n × nₑ$-dimensional (usually Gaussian) random matrix. This term represents an inflation to add to the ensemble to prevent the ensemble spread from collapsing after just a few EnKF iterations.

abstract type EnsembleCovStyle 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

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

is generally expensive, and thus alternative strategies are sought. Subtypes:

• StandardCovUpdate
• NonstandardCovUpdate

const EnsembleKalmanFilter{A<:Model,B<:Model,C<:Ensemble,D<:Ensemble} = KalmanFilter{A,B,C,D}

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.
• other than the different treatment of the state's covariance, the other steps (transition, observation, innovation and Kalman gain) are equivalent to a standard Kalman Filter.

Subtypes:

• EnKF
• DEnKF

const ExtendedKalmanFilter{C<:Distribution,D<:Distribution} = 
  KalmanFilter{<:StochasticLinearisedModel,<:StochasticLinearisedModel,C,D}

Implements the Extended Kalman Filter (EKF). In particular:

• for the forecasting step: we linearise the transition operator around the analysed state of the

previous iteration;

• for the analysis step: we linearise the observation operator around the forecasted state of the

current iteration. The remaining scheme is equivalent to that of a standard Kalman Filter. From an implementation standpoint, an ExtendedKalmanFilter simply requires the transition and observation models to both be StochasticLinearisedModel.

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

const FirstMoment = Distribution{1}

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

abstract type FunctionFilter <: Filter end

Subtype reserved for filters whose transition and observation models are time-dependent functions. In practice, the only difference with a standard Filter is that the models of a FunctionFilter must be explicitly evaluated at each iteration before running the Kalman iterations (see loop for more details).

struct FunctionKalmanFilter{A<:Function,B<:Function,C<:Distribution,D<:Distribution} <: FunctionFilter
  transition::A 
  observation::B
  prior::C
  obs_prior::D
  cache::KalmanCache
end

Filter subtype implementing a Kalman filter procedure. Fields:

• transition: Real -> Model function representing the transition operator. The real input it receives

could be, for example, the time instant of the current Kalman iteration. This field should be evaluated at each iteration to successfully run the Kalman iterations, e.g. via loop;

• observation: Real -> Model function representing the observation operator. The real input it receives

could be, for example, the time instant of the current Kalman iteration. This field should be evaluated at each iteration to successfully run the Kalman iterations, e.g. via loop;

• prior: Distribution representing the probability distribution for the state;
• obs_prior: Distribution representing the probability distribution for the observation;
• cache: cached object allowing for efficient in-place operations.

struct GenericFirstMoment{T,A<:AbstractVector{T}} <: 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 GenericSecondMoment{T,A<:AbstractVector{T},B<:AbstractMatrix{T}} <: SecondMoment
  mean::A 
  covariance::B
end

Most basic implementation of a SecondMoment distribution.

struct KalmanFilter{A<:Model,B<:Model,C<:Distribution,D<:Distribution} <: Filter
  transition::A 
  observation::B
  prior::C
  obs_prior::D
  cache::KalmanCache
end

Filter subtype implementing a Kalman filter procedure. Fields:

• transition: Model representing the transition operator;
• observation: Model representing the observation operator;
• prior: Distribution representing the probability distribution for the state;
• obs_prior: Distribution representing the probability distribution for the observation;
• cache: cached object allowing for efficient in-place operations.

struct Linear <: ModelStyle 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 $P$, then the output

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

struct LinearisedModel{T,A<:AbstractMatrix{T},F<:FType} <: LinearModel
  form::F
  cache::A
end

Type reserved for (generally nonlinear) a function or Gridap Map form that is linearised around some point $x$ (to be later specified). The $x$-dependent Jacobian should be stored in-place in the field cache.

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

Type used for operator-like quantities, such as functions or Gridap Maps. For performance reasons, we distinguish models depending on their ModelStyle 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 Distribution 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 Distribution determine the expression of posterior.

abstract type ModelStyle end

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

struct Nonlinear <: ModelStyle 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 $P$, then the output

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

abstract type NonstandardCovUpdate <: EnsembleCovStyle end

Trait used for ensembles whose covariance is generally not directly incorporated in the filtering procedure. Subtypes:

• EnKFUpdate
• DEnKFUpdate

const SecondMoment = Distribution{2}

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

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

This SecondMoment distribution represents the sigma points needed to run the UnscentedTransform. 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 StandardCovUpdate <: EnsembleCovStyle end

Standard computation of the ensemble covariance, according to the formula:

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

where $μ$ is the $n$-dimensional ensemble mean, and ensemble is the $n × nₑ$ ensemble matrix. This formula is highly expensive, depending on the value of $n$, and should be used only for the observations ensemble.

struct StochasticModel{A<:ModelStyle,B<:Model{A},C<:Distribution,D<:NoiseStrategy} <: Model{A}
  model::B
  noise::C
  strategy::D
end

Models characterised by an underlying deterministic Model model, and a stochastic noise component, as specified by the field noise. Usually, noise is a SecondMoment distribution with zero mean and a certain covariance Q. The field strategy determines how the stochastic component is added to the deterministic component. Suppose that

\[θ ∼ SecondMoment(η,R),\]

where $θ$ is the output distribution such that

\[θ = model(d)\]

for a given input distribution

\[d ∼ SecondMoment(μ,P),\]

Then if:

• strategy::Default (default): we augment $μ ← μ + mean(noise)$, and $P ← P + cov(noise)$;
• strategy::Additive: we augment $μ ← μ + mean(noise) + ω$, and $P ← P + cov(noise)$, where

$ω$ is a random vector drawn according to noise.

const UnscentedTransform{A<:Model,B<:Model} = KalmanFilter{A,B,<:SigmaPoints,<:SigmaPoints}

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

• for the forecasting step: we update the expression of the sigma points based on the analysed sigma points

from the previous iteration; we evaluate the state prior distribution in the sigma points, and update its mean and covariance according to weighted linear combinations (see SigmaPoints);

• for the analysis step: we evaluate the obvervation prior distribution in the sigma points, and

update its mean and covariance according to the same weighted linear combinations; we compute innovation, Kalman gain, and state posterior distribution as in a standard Kalman Filter scheme. From an implementation standpoint, an UnscentedTransform simply requires the transition and observation priors to both be SigmaPoints.

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

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 Kalman gain (see kalman_gain!)
• Computing the innovation (see innovation!)
• Updating the forecasted distribution by accounting for the Kalman gain (see update!)

To run a single iteration of the Kalman filter, one must run the forecasting step in forecast! prior to the analysis one.

dimension(a::Model) -> Int

If a is a Model encoding an operator from $Rᵐ$ to $Rⁿ$, it returns the integer $n$.

dimension(d::Distribution) -> Int

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

dimension(a::Model) -> Int

If a is a Model encoding an operator from $Rᵐ$ to $Rⁿ$, it returns the integer $m$.

draw(d::Distribution) -> AbstractVector 
draw(d::Distribution,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.

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

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.

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. Though $H$ need not necessarily be stochastic, the standard implementation is such that $H$ is a StochasticModel.

get_observation_prior(f::Filter) -> Distribution

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

get_prior(f::Filter) -> Distribution

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. Though $F$ need not necessarily be stochastic, the standard implementation is such that $F$ is a StochasticModel.

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

Given an observation z, returns the innovation ỹ such that

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

where $yₙ$ represents the observation forecasted by the filter f.

jac(a::Model,x::InType) -> AbstractMatrix

Returns a's Jacobian matrix evaluated in $x$.

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

In-place computation of the Kalman gain K according to the formula $K = Pxy * Pyy⁻¹$ where Pxy and Pyy are the state-observation and observation covariance matrices, respectively. K is storede as a cached object in f, while the covariances Pxy and Pyy 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 -> AbstractVector{<:Distribution}

Given a filter f and a list of observations obs, iteratively runs the forecast-analyse paradigm typical of a Kalman filter, producing a list of posterior distributions of the state variable. In practice, one iteration of the loop consists of one call to forecast!, followed by one to analyse!. The posterior resulting from each analysis is then fed as the prior distribution to the next forecast step.

mixed_cov!(P::AbstractMatrix,f::Filter,posterior::Distribution) -> AbstractMatrix

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

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:

\[P = ∑ᵢ₌₁²ᴸ⁺¹ 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::Distribution) -> Distribution

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.

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 + λ)P) χ[:,L+2:2L+1] = μ - √((L + λ)P)\]

where μ and P 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\]

similar_distribution(d::Distribution,dim::Int=dimension(d)) -> Distribution

Returns a distribution of same type as d, with a possibly different dimension specified by the optional argument dim.

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

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::Ensemble) -> Ensemble

Update the mean, anomaly and covariance of the ensemble d.

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

Update the mean and covariance of the sigma points d.

update!(posterior::Distribution,f::Filter,args...) -> Distribution

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.

visualize(
  history::AbstractVector{<:Distribution},
  grid=eachindex(history);
  index::Int=1
)

visualize(
  true_values::AbstractMatrix,
  history::AbstractVector{<:Distribution},
  grid=eachindex(history);
  index::Int=1
)

Plot the historical distributions obtained by running the Kalman iterations, for e.g. via the function loop. The primary estimator is the mean of the distributions. If the distributions feature a second moment (i.e. they are equipped by a variance) then a confidence interval is also drawn. The true data true_values, if known, may be provided, and will be plotted on the same figure.