Public Documentation

Main module for ParameterEstimocean.jl – a Julia software for parameter estimation for Oceananigans models using Ensemble Kalman Inversion.

Transformation(; time=nothing, space=nothing, normalization=nothing)

Return a transformation that is applied on the observation. Examples include slicing the data or multiplying with weight factors to make the loss function putting more weight in particular regions of the domain or particular times. Also, we can denote a normalization procedure applied to the data after the space- and time- transformations.

Slicing is prescribed as SpaceIndices and TimeIndices. For example,

Transformation(time = TimeIndices(4:10))

only keeps time instances 4 to 10 from the observations. Similarly,

Transformation(space = SpaceIndices(x=:, y=1:10, z=2:2:20))

does not affect the x dimension of the data, but slices the observations in y and z as prescribed.

Keyword Arguments

• time: The time transformation either as a TimeIndices or as an AbstractVector of weights of same size as observations.times. If nothing is given, then, by default, the transformation ignores the first snapshot (initial state). (Default: nothing)

• space: The space trasformation either as a SpaceIndices or as an AbstractArray of weights of same size as a snapshot of the observations. (Default: nothing)

• normalization: The normalization that is applied to the data after space and time transformations have been applied first. (Default: nothing)

ZScore(data)

Return the ZScore normalization of a FieldTimeSeries after computing its mean and its variance.

SyntheticObservations(path=nothing;
                      field_names,
                      forward_map_names = field_names,
                      transformation = Transformation(),
                      times = nothing,
                      field_time_serieses = nothing,
                      architecture = CPU(),
                      regrid = nothing)

Return a time series of synthetic observations generated by Oceananigans.jl's simulations gridded as Oceananigans.jl fields.

struct FreeParameters{N, P, D}

A container for free parameters that includes the parameter names and their corresponding prior distributions.

• names

  free parameters

• priors

  prior distributions for free parameters

• dependent_parameters

  dependent parameters

FreeParameters(priors; names = Symbol.(keys(priors)), dependent_parameters=NamedTuple())

Return named FreeParameters with priors. Free parameter names are inferred from the keys of priors if not provided. Optionally, dependent_parameters are prescribed as a NamedTuple whose keys are the names of "additional" parameters, and whose values are functions that return those parameters given a vector of free parameters in names.

Example

julia> using Distributions, ParameterEstimocean

julia> priors = (ν = Normal(1e-4, 1e-5), κ = Normal(1e-3, 1e-5))
(ν = Normal{Float64}(μ=0.0001, σ=1.0e-5), κ = Normal{Float64}(μ=0.001, σ=1.0e-5))

julia> free_parameters = FreeParameters(priors)
FreeParameters with 2 parameters
├── names: (:ν, :κ)
├── priors: Dict{Symbol, Any}
│   ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
│   └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)
└── dependent parameters: Dict{Symbol, Any}

julia> c(p) = p.ν + p.κ # compute a third dependent parameter `c` as a function of `ν` and `κ`
c (generic function with 1 method)

julia> free_parameters_with_a_dependent = FreeParameters(priors, dependent_parameters=(; c))
FreeParameters with 2 parameters and 1 dependent parameter
├── names: (:ν, :κ)
├── priors: Dict{Symbol, Any}
│   ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
│   └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)
└── dependent parameters: Dict{Symbol, Any}
    └── c => c (generic function with 1 method)

ScaledLogitNormal([FT=Float64;] bounds=(0, 1), mass=0.5, interval=nothing, μ=nothing, σ=nothing)

Return a ScaledLogitNormal distribution with compact support within bounds.

interval is an optional 2-element tuple or Array. When specified, the parameters μ and σ of the underlying Normal distribution are calculated so that mass fraction of the probability density lies within interval.

If interval is not specified, then μ=0 and σ=1 by default.

Notes

ScaledLogitNormal is a four-parameter distribution generated by the transformation

\[Y = L + (U - L) / [1 + \exp(X)],\]

of the normally-distributed variate $X ∼ 𝒩(μ, σ)$. The four parameters governing the distribution of $Y$ are thus

• $L$: lower bound (0 for the LogitNormal distribution)
• $U$: upper bound (1 for the LogitNormal distribution)
• $μ$: mean of the underlying Normal distribution
• $σ²$: variance of the underlying Normal distribution

lognormal(; mean, std)

Return Lognormal distribution parameterized by the distribution mean and standard deviation std.

Notes

A variate X is LogNormal distributed if

\[\log(X) ∼ 𝒩(μ, σ²) ,\]

where $𝒩(μ, σ²)$ is the Normal distribution with mean $μ$ and variance $σ²$.

The mean and variance $s²$ (where $s$ is the standard deviation or std) are related to the parameters $μ$ and $σ²$ via

\[ m = \exp(μ + σ² / 2),\]

\[s² = [\exp(σ²) - 1] m².\]

These formula allow us to calculate $μ$ and $σ$ given $m$ and $s²$, since rearranging the formula for $s²$ gives

\[\exp(σ²) = m² / s² + 1\]

which then yields

\[σ = \sqrt{\log(m² / s² + 1)}.\]

We then find that

\[μ = \log(m) - σ² / 2 .\]

See also wikipedia.

struct ConcatenatedOutputMap

Forward map transformation of simulation output to the concatenated vectors of the simulation output.

InverseProblem(observations,
               simulation,
               free_parameters;
               output_map = ConcatenatedOutputMap(),
               time_series_collector = nothing,
               initialize_simulation = nothingfunction)

Return an InverseProblem.

forward_map(ip, parameters)

Run ip.simulation forward with parameters and return the data, transformed into an array format expected by EnsembleKalmanProcesses.jl.

forward_run!(ip, parameters)

Initialize ip.simulation with parameters and run it forward. Output is stored in ip.time_series_collector.

observation_map(ip::InverseProblem)

Transform and return ip.observations appropriate for ip.output_map.

observation_map_variance_across_time(map::ConcatenatedOutputMap, observation::SyntheticObservations)

Return an array of size (Nensemble, Ny * Nz * Nfields, Ny * Nz * Nfields) that stores the covariance of each element of the observation map measured across time, for each ensemble member, where Nensemble is the ensemble size, Ny is either the number of grid elements in y or the batch size, Nz is the number of grid elements in the vertical, and Nfields is the number of fields in observation.

EnsembleKalmanInversion(inverse_problem;
                        noise_covariance = 1,
                        pseudo_stepping = nothing,
                        resampler = Resampler(),
                        unconstrained_parameters = nothing,
                        forward_map_output = nothing,
                        mark_failed_particles = NormExceedsMedian(1e9),
                        ensemble_kalman_process = Inversion())

Return an object that finds local minima of the inverse problem:

\[y = G(θ) + η,\]

for the parameters $θ$, where $y$ is a vector of observations (often normalized), $G(θ)$ is a forward map that predicts the observations, and $η ∼ 𝒩(0, Γ_y)$ is zero-mean random noise with a noise_covariance matrix $Γ_y$ representing uncertainty in the observations.

The "forward map output" G is model output mapped to the space of inverse_problem.observations.

(For more details on the Ensemble Kalman Inversion algorithm refer to the EnsembleKalmanProcesses.jl Documentation.)

Positional Arguments

• inverse_problem (InverseProblem): Represents an inverse problem representing the comparison between synthetic observations generated by Oceananigans.jl and model predictions, also generated by Oceananigans.jl.

Keyword Arguments

• noise_covariance (Number or AbstractMatrix): Covariance matrix representing observational uncertainty. noise_covariance::Number is converted to a scaled identity matrix.

• pseudo_stepping: The pseudo time-stepping scheme for stepping EKI forward.

• resampler: controls particle resampling procedure. See Resampler.

• unconstrained_parameters: Default: nothing.

• forward_map_output: Default: nothing.

• ensemble_kalman_process: Process type defined by EnsembleKalmanProcesses.jl. Default: Inversion().

NormExceedsMedian(minimum_relative_norm = 1e9)

The particle failure condition. A particle is marked "failed" if the forward map norm is larger than minimum_relative_norm times more than the median value of the ensemble. By default minimum_relative_norm = 1e9.

Return a BitVector indicating whether the norm of the forward map for a given particle exceeds the median by mrn.minimum_relative_norm.

iterate!(eki::EnsembleKalmanInversion;
         iterations = 1,
         pseudo_Δt = eki.pseudo_Δt,
         pseudo_stepping = eki.pseudo_stepping,
         show_progress = true)

Iterate the ensemble Kalman inversion problem eki forward by iterations.

Keyword arguments

• iterations (Int): Number of iterations to run. (Default: 1)

• pseudo_Δt (Float64): Pseudo time-step. When convegence_rate is specified, this is an initial guess for finding an adaptive time-step. (Default: eki.pseudo_Δt)

• pseudo_stepping (Float64): Ensemble convergence rate for adaptive time-stepping. (Default: eki.pseudo_stepping)

• show_progress (Boolean): Whether to show a progress bar. (Default: true)

Return

• best_parameters: the ensemble mean of all parameter values after the last iteration.