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.

BatchedSyntheticObservations(observations; weights)

Return a collection of observations with weights, where observations is a Vector or Tuple of SyntheticObservations. weights are unity by default.

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 free parameters and 0 dependent parameters
├── names: (:ν, :κ)
└── priors:
    ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
    └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)

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 free parameters and 1 dependent parameter
├── names: (:ν, :κ)
├── priors:
│   ├── ν => Normal{Float64}(μ=0.0001, σ=1.0e-5)
│   └── κ => Normal{Float64}(μ=0.001, σ=1.0e-5)
└── dependent parameters:
    └── 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.

BatchedInverseProblem(batched_ip; weights=Tuple(1 for o in batched_ip))

Return a collection of observations with weights, where observations is a Vector or Tuple of SyntheticObservations. weights are unity by default.

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,
               initialize_with_observations = true)

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::InverseProblem, parameter_ensemble)

Initialize ip.simulation with parameter_ensemble and run it forward. Output is stored in ip.time_series_collector. forward_run can also be called with one parameter set.

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,
                        pseudo_Δt = 1.0,
                        resampler = Resampler(),
                        unconstrained_parameters = nothing,
                        forward_map_output = nothing,
                        mark_failed_particles = NormExceedsMedian(1e9),
                        ensemble_kalman_process = Inversion(),
                        Nensemble = Nensemble(inverse_problem),
                        tikhonov = false)

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.

• pseudo_Δt: The pseudo time-step; Default: 1.0.

• resampler: controls particle resampling procedure. See Resampler.

• unconstrained_parameters: Default: nothing.

• forward_map_output: Default: nothing.

• mark_failed_particles: The particle failure condition. Default: NormExceedsMedian(1e9).

• ensemble_kalman_process: The Ensemble Kalman process. Default: Inversion().

• tikhonov: Whether to incorporate prior information in the EKI objective via Tikhonov regularization. See Chada et al. "Tikhonov Regularization Within Ensemble Kalman Inversion." SIAM J. Numer. Anal. 2020.

struct NormExceedsMedian{T}

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.

ObjectiveLossThreshold(multiple = 4.0;
                       baseline = nanmedian,
                       distance = median_absolute_deviation)

Return a failure criterion that defines failure for particle k as

Φₖ > baseline(Φ) + multiple * distance(Φ)

where Φ is the objective loss function.

By default, baseline = nanmedian, the distance is the median absolute deviation, and multiple = 4.0.

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

Convenience function for running pseudo_step! multiple times with the same argument. Iterates the ensemble Kalman inversion dynamic forward by iterations given current state eki.

Keyword arguments

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

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

• kwargs (NamedTuple): Keyword arguments to be passed to pseudo_step! at each iteration. (Defaults: see pseudo_step!)

Return

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

pseudo_step!(eki::EnsembleKalmanInversion; 
             pseudo_Δt = nothing,
             pseudo_stepping = nothing,
             covariance_inflation = 0.0,
             momentum_parameter = 0.0)

Step forward X = eki.unconstrained_parameters using y = eki.mapped_observations, Γy = eki.noise_covariance, and G = eki.forward_map_output.

Keyword arguments

• pseudo_Δt (Float64 or Nothing): Pseudo time-step. If pseudo_Δt is nothing, the time step is set according to the algorithm specified by the pseudo_stepping scheme; If pseudo_Δt is a Float64, pseudo_stepping is ignored. (Default: nothing)

• pseudo_stepping (Float64): Scheme for selecting a time step if pseudo_Δt is nothing. (Default: eki.pseudo_stepping)

• covariance_inflation: (Default: 0.)

• momentum_parameter: (Default: 0.)

ConstantConvergence(; convergence_ratio=0.7)

Return the ConstantConvergence pseudo-stepping scheme with target convergence_ratio. With ConstantConvergence, the ensemble Kalman inversion (EKI) pseudo step size is adjusted such that the $n$-th root of the determinant of the parameter covariance is decreased by convergence_ratio (e.g., for convergence_ratio = 0.7, by 70%) after one EKI iteration, where $n$ is the number of parameters.