ParticleDA.jl

run_particle_filter(
    init_model,
    input_file_path,
    observation_file_path,
    filter_type=BootstrapFilter,
    summary_stat_type=MeanAndVarSummaryStat;
    rng=Random.TaskLocalRNG()
) -> Tuple{Matrix, Union{NamedTuple, Nothing}}

Run particle filter. init_model is the function which initialise the model, input_file_path is the path to the YAML file with the input parameters. observation_file_path is the path to the HDF5 file containing the observation sequence to perform filtering for. filter_type is the particle filter type to use. See ParticleFilter for the possible values. summary_stat_type is a type specifying the summary statistics of the particles to compute at each time step. See AbstractSummaryStat for the possible values. rng is a random number generator to use to generate random variates while filtering - a seeded random number generator may be specified to ensure reproducible results. If running with multiple threads a thread-safe generator such as Random.TaskLocalRNG (the default) must be used.

Returns a tuple containing the state particles representing an estimate of the filtering distribution at the final observation time (with each particle a column of the returned matrix) and a named tuple containing the estimated summary statistics of this final filtering distribution. If running on multiple ranks using MPI, the returned states array will correspond only to the particles local to this rank and the summary statistics will be returned only on the master rank with all other ranks returning nothing for their second return value.

simulate_observations_from_model(
    init_model, input_file_path, output_file_path; rng=Random.TaskLocalRNG()
) -> Matrix

Simulate observations from the state space model initialised by the init_model function with parameters specified by the model key in the input YAML file at input_file_path and save the simulated observation and state sequences to a HDF5 file at output_file_path. rng is a random number generator to use to generate random variates while simulating from the model - a seeded random number generator may be specified to ensure reproducible results.

The input YAML file at input_file_path should have a simulate_observations key with value a dictionary with keys seed and n_time_step corresponding to respectively the number of time steps to generate observations for from the model and the seed to use to initialise the state of the random number generator used to simulate the observations.

The simulated observation sequence is returned as a matrix with columns corresponding to the observation vectors at each time step.

ParticleFilter

Abstract type for particle filters. Currently implemented subtypes are:

• BootstrapFilter
• OptimalFilter

BootstrapFilter <: ParticleFilter

Singleton type BootstrapFilter. This can be used as argument of run_particle_filter to select the bootstrap filter.

OptimalFilter <: ParticleFilter

Singleton type OptimalFilter. This can be used as argument of run_particle_filter to select the optimal proposal filter (for conditionally linear-Gaussian models).

AbstractSummaryStat

Abstract type for summary statistics of particle ensemble. Concrete subtypes can be passed as the filter_type argument to run_particle_filter to specify which summary statistics to record and how they are computed.

See also: AbstractSumReductionSummaryStat, AbstractCustomReductionSummaryStat.

AbstractCustomReductionSummaryStat <: AbstractSummaryStat

Abstract type for summary statistics computed using custom MPI reductions. Allows greater flexibility in computing statistics which can support more numerically stable implementations, but at a cost of not being compatible with all CPU architectures. In particular, MPI.jl does not currently support custom operators on Power PC and ARM architecures.

AbstractSumReductionSummaryStat <: AbstractSummaryStat

Abstract type for summary statistics computed using standard MPI sum reductions. Compatible with a wider range of CPU architectures but may require less numerically stable implementations.

MeanAndVarSummaryStat <: AbstractCustomReductionSummaryStat

Custom reduction based summary statistic type which computes the means and variances o the particle ensemble for each state dimension. On CPU architectures which do not support custom reductions NaiveMeanAndVarSummaryStat can be used instead.

MeanSummaryStat <: AbstractCustomReductionSummaryStat

Custom reduction based summary statistic type which computes the means of the particle ensemble for each state dimension. On CPU architectures which do not support custom reductions NaiveMeanSummaryStat can be used instead.

NaiveMeanAndVarSummaryStat <: AbstractSumReductionSummaryStat

Sum reduction based summary statistic type which computes the means and variances of the particle ensemble for each state dimension. The mean and variance are computed by directly accumulating the sums of the particle values, the squared particle values and number of particles on each rank, with the variance computed as the scaled difference between the sum of the squares and square of the sums. This 'naive' implementation avoids custom MPI reductions but can be numerically unstable for large ensembles or state components with large values. If custom reductions are supported by the CPU architecture in use the more numerically stable MeanAndVarSummaryStat should be used instead.

NaiveMeanSummaryStat <: AbstractSumReductionSummaryStat

Sum reduction based summary statistic type which computes the means of the particle ensemble for each state dimension. The mean is computed by directly accumulating the sums of the particle values and number of particles on each rank. If custom reductions are supported by the CPU architecture in use the more numerically stable MeanSummaryStat should be used instead.

ParticleDA.get_state_dimension(model) -> Integer

Return the positive integer dimension of the state vector which is assumed to be fixed for all time steps.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.get_observation_dimension(model) -> Integer

Return the positive integer dimension of the observation vector which is assumed to be fixed for all time steps.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.get_state_eltype(model) -> Type

Return the element type of the state vector which is assumed to be fixed for all time steps.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.get_observation_eltype(model) -> Type

Return the element type of the observation vector which is assumed to be fixed for all time steps.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.sample_initial_state!(state, model, rng, task_index=1)

Sample value for state vector from its initial distribution for model described by model using random number generator rng to generate random draws and writing to state argument.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.update_state_deterministic!(state, model, time_index, task_index=1)

Apply the deterministic component of the state time update at discrete time index time_index for the model described by model for the state vector state writing the updated state back to the state argument.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.update_state_stochastic!(state, model, rng, task_index=1)

Apply the stochastic component of the state time update for the model described by model for the state vector state, using random number generator rng to generate random draws and writing the updated state back to the state argument.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.sample_observation_given_state!(
    observation, state, model, rng, task_index=1
)

Simulate noisy observations of the state state of model described by model and write to observation array using rng to generate any random draws.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.get_log_density_observation_given_state(
    observation, state, model, task_index=1
) -> Real

Return the logarithm of the probability density of an observation vector observation given a state vector state for the model associated with model. Any additive terms that are constant with respect to the state may be neglected.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.write_model_metadata(file::HDF5.File, model)

Write metadata for with the model described by model to the HDF5 file file.

This method is intended to be extended by the user with the above signature, specifying the type of model.

ParticleDA.get_observation_mean_given_state!(
    observation_mean, state, model, task_index=1
)

Compute the mean of the multivariate normal distribution on the observations given the current state and write to the first argument.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering conditionally Gaussian models with the optimal proposal filter implementation in OptimalFilter.

ParticleDA.get_covariance_state_noise(model, i, j) -> Real

Return covariance cov(U[i], U[j]) between components of the zero-mean Gaussian state noise vector U.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering conditionally Gaussian models with the optimal proposal filter implementation in OptimalFilter.

ParticleDA.get_covariance_observation_noise(model, i, j) -> Real

Return covariance cov(V[i], V[j]) between components of the zero-mean Gaussian observation noise vector V.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering conditionally Gaussian models with the optimal proposal filter implementation in OptimalFilter.

ParticleDA.get_covariance_state_observation_given_previous_state(
    model, i, j
) -> Real

Return the covariance cov(X[i], Y[j]) between components of the state vector X = F(x) + U and observation vector Y = H * X + V where H is the linear observation operator, F the (potentially non-linear) forward operator describing the deterministic state dynamics, U is a zero-mean Gaussian state noise vector, V is a zero-mean Gaussian observation noise vector and x is the state at the previous observation time.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering conditionally Gaussian models with the optimal proposal filter implementation in OptimalFilter.

ParticleDA.get_covariance_observation_observation_given_previous_state(
    model, i, j
) -> Real

Return covariance cov(Y[i], Y[j]) between components of the observation vector Y = H * (F(x) + U) + V where H is the linear observation operator, F the (potentially non-linear) forward operator describing the deterministic state dynamics, U is a zero-mean Gaussian state noise vector, V is a zero-mean Gaussian observation noise vector and x is the state at the previous observation time.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering conditionally Gaussian models with the optimal proposal filter implementation in OptimalFilter.

ParticleDA.get_state_indices_correlated_to_observations(
    model
) -> AbstractVector{Int}

Return the vector containing the indices of the state vector X which at least one of the observations Yare correlated to, that isi ∈ 1:getstatedimension(model)such thatcov(X[i], Y[j]) > 0for at least onej ∈ 1:getobservationdimension(model). This is used to avoid needing to compute and store zero covariance terms. Defaults to returning all state indices which will always give correct results but will be inefficient if there are zero blocks incov(X, Y)`.

ParticleDA.get_initial_state_mean!(state_mean, model)

Compute the mean of the multivariate normal distribution on the initial state and write to the first argument.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering linear-Gaussian models with the Kalman filter for testing.

ParticleDA.write_snapshot(
    output_filename, model, filter_data, states, time_index, save_states
)

Write a snapshot of the model and filter states to the HDF5 file output_filename for the model and filters described by model and filter_data respectively at time index time_index, optionally saving the current ensemble of state particles represented by the two-dimensional array states (first axis state component, second particle index) if save_states == true. time_index == 0 corresponds to the initial model and filter states before any updates and non-time dependent model data will be written out when called with this value of time_index.

ParticleDA.write_state(
    file::HDF5.File,
    state::AbstractVector{T},
    time_index::Int,
    group_name::String,
    model
)

Write the model state at time index time_index represented by the vector state to the HDF5 file file with group_name for the model represented by model.

ParticleDA.get_covariance_initial_state(model, i, j) -> Real

Return covariance cov(X[i], X[j]) between components of the initial state vector X.

This method is intended to be extended by the user with the above signature, specifying the type of model. Only required for filtering linear-Gaussian models with the Kalman filter for testing.

ParticleDA.init_filter(filter_params, model, nprt_per_rank, ::T) -> NamedTuple

Initialise any data structures required by filter of type T, with filtering specific parameters specified by filter_params, state space model to perform filtering with described by model and nprt_per_rank particles per MPI rank.

New filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.

ParticleDA.sample_proposal_and_compute_log_weights!(
    states, log_weights, observation, model, filter_data, ::T, rng
)

Sample new values for two-dimensional array of state vectors states from proposal distribution writing in-place to states array and compute logarithm of unnormalized particle weights writing to log_weights given observation vector observation, with state space model described by model, named tuple of filter specific data structures filter_data, filter type T and random number generator rng used to generate any random draws.

New filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.

ParticleDA.write_observation(
    file::HDF5.File,
    observation::AbstractVector,
    time_index::Int,
    model
)

Write the observations at time index time_index represented by the vector observation to the HDF5 file file for the model represented by model.

ParticleDA.write_weights(
    file::HDF5.File,
    weights::AbstractVector{T},
    time_index::Int,
    model
)

Write the particle weights at time index time_index represented by the vector weights to the HDF5 file file for the model represented by model.

ParticleDA.get_initial_state_mean(model)

Compute the mean of the multivariate normal distribution on the initial state.

FilterParameters()

Parameters for ParticleDA run. Keyword arguments:

• master_rank : Id of MPI rank that performs serial computations
• verbose::Bool : Flag to control whether to write output
• output_filename::String : Name of output file
• nprt::Int : Number of particles for particle filter
• enable_timers::Bool : Flag to control run time measurements
• particle_save_time_indices: Set of time indices to save particles at
• seed: Seed to initialise state of random number generator used for filtering
• n_tasks: Number of tasks to use for running parallelisable operations. Positive integers indicate the number of tasks directly, while the absolute value of negative integers indicate the number of task to use per-thread (as reported by Threads.nthreads()). Using multiple tasks per thread will improve the ability of the scheduler to balance load across threads but potentially increase overheads. If simulation of the model being filtered use multiple threads then it may be beneficial to set the n_tasks = 1 to avoid too much contention between threads.

LLW2dModelParameters()

Parameters for the linear long wave two-dimensional (LLW2d) model. Keyword arguments:

• nx::Int : Number of grid points in the x direction
• ny::Int : Number of grid points in the y direction
• x_length::AbstractFloat : Domain size in metres in the x direction
• y_length::AbstractFloat : Domain size in metres in the y direction
• dx::AbstractFloat : Distance in metres between grid points in the x direction
• dy::AbstractFloat : Distance in metrtes between grid points in the y direction
• station_filename::String : Name of input file for station coordinates
• n_stations_x::Int : Number of observation stations in the x direction (if using regular grid)
• n_stations_y::Int : Number of observation stations in the y direction (if using regular grid)
• station_distance_x::Float : Distance in metres between stations in the x direction (if using regular grid)
• station_distance_y::Float : Distance in metres between stations in the y direction (if using regular grid)
• station_boundary_x::Float : Distance in metres between bottom left edge of box and first station in the x direction (if using regular grid)
• station_boundary_y::Float : Distance in metres between bottom left edge of box and first station in the y direction (if using regular grid)
• n_integration_step::Int : Number of sub-steps to integrate the forward model per time step
• time_step::AbstractFloat : Time step length in seconds
• peak_position::Vector{AbstractFloat} : The `[x, y] coordinates in metres of the initial wave peak
• peak_height::AbstractFloat : The height in metres of the initial wave peak
• source_size::AbstractFloat : Cutoff distance in metres from the peak for the initial wave
• bathymetry_setup::AbstractFloat : Bathymetry set-up
• lambda::AbstractFloat : Length scale for Matérn covariance kernel in background noise
• nu::AbstractFloat : Smoothess parameter for Matérn covariance kernel in background noise
• sigma::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in background noise
• lambda_initial_state::AbstractFloat : Length scale for Matérn covariance kernel in initial state of particles
• nu_initial_state::AbstractFloat : Smoothess parameter for Matérn covariance kernel in initial state of particles
• sigma_initial_state::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in initial state of particles
• padding::Int : Min padding for circulant embedding gaussian random field generator
• primes::Int: Whether the size of the minimum circulant embedding of the covariance matrix can be written as a product of small primes (2, 3, 5 and 7). Default is true.
• use_peak_initial_state_mean::Bool: Whether to set mean of initial height field to a wave peak (true) or to all zeros (false). In both cases the initial mean of the other state variables is zero.
• absorber_thickness_fraction::Float : Thickness of absorber for sponge absorbing boundary conditions, fraction of grid size
• boundary_damping::Float : Damping for boundaries
• cutoff_depth::Float : Shallowest water depth
• obs_noise_std::Vector: Standard deviations of noise added to observations of the true state
• observed_state_var_indices::Vector: Vector containing the indices of the observed state variables (1: height, 2: velocity x-component, 3: velocity y-component)