From e66db05b20ef97ee701418a46a1c8ec803f1b91a Mon Sep 17 00:00:00 2001 From: "Documenter.jl" Date: Wed, 18 Dec 2024 12:50:14 +0000 Subject: [PATCH] build based on 238af83 --- dev/.documenter-siteinfo.json | 2 +- dev/index.html | 53 +++++++++++++++++++++++----------- dev/objects.inv | Bin 868 -> 1033 bytes dev/search_index.js | 2 +- 4 files changed, 38 insertions(+), 19 deletions(-) diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index e8f3d8a..5862c8c 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.11.1","generation_timestamp":"2024-11-16T10:29:32","documenter_version":"1.8.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.11.2","generation_timestamp":"2024-12-18T12:50:08","documenter_version":"1.8.0"}} \ No newline at end of file diff --git a/dev/index.html b/dev/index.html index a03a712..dfa9de6 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,46 +1,57 @@ -ParticleDA.jl · ParticleDA
GitHub

ParticleDA.jl

ParticleDA.jl is a Julia package to perform data assimilation with particle filters, distributed using MPI.

Installation

To install the latest stable version of the package, open the Julia REPL, enter the package manager with ], then run the command

add ParticleDA.jl

If you plan to develop the package (make changes, submit pull requests, etc), in the package manager mode run this command

dev ParticleDA

This will automatically clone the repository to your local directory ~/.julia/dev/ParticleDA.

You can exit from the package manager mode by pressing CTRL + C or, alternatively, the backspace key when there is no input in the prompt.

Usage

After installing the package, you can start using it in Julia's REPL with

using ParticleDA

The run the particle filter you can use the function run_particle_filter:

ParticleDA.run_particle_filter — Function
run_particle_filter(
+ParticleDA.jl · ParticleDA
GitHub
ParticleDA.jl
ParticleDA.jl is a Julia package to perform data assimilation with particle filters, distributed using MPI.
Installation
To install the latest stable version of the package, open the Julia REPL, enter the package manager with ], then run the command
add ParticleDA.jl
If you plan to develop the package (make changes, submit pull requests, etc), in the package manager mode run this command
dev ParticleDA
This will automatically clone the repository to your local directory ~/.julia/dev/ParticleDA.
You can exit from the package manager mode by pressing CTRL + C or, alternatively, the backspace key when there is no input in the prompt.
Usage
After installing the package, you can start using it in Julia's REPL with
using ParticleDA
The run the particle filter you can use the function run_particle_filter:
ParticleDA.run_particle_filter — Function
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.
source
To simulate observations from the model (which can be used to for example test the filtering algorithms) you can use the function simulate_observations_from_model
ParticleDA.simulate_observations_from_model — Function
simulate_observations_from_model(
+) -> 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.
source
To simulate observations from the model (which can be used to for example test the filtering algorithms) you can use the function simulate_observations_from_model
ParticleDA.simulate_observations_from_model — Function
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.
source
The filter_type argument to run_particle_filter should be a concrete subtype of the ParticleFilter abstract type.
The summary_stat_type argument to run_particle_filter should be a concrete subtype of the AbstractSummaryStat abstract type.
ParticleDA.AbstractCustomReductionSummaryStat — Type
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.
source
ParticleDA.AbstractSumReductionSummaryStat — Type
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.
source
ParticleDA.MeanAndVarSummaryStat — Type
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.
source
ParticleDA.MeanSummaryStat — Type
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.
source
ParticleDA.NaiveMeanAndVarSummaryStat — Type
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.
source
ParticleDA.NaiveMeanSummaryStat — Type
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.
source
The next section details how to write the interface between the model and the particle filter.
Interfacing the model
The model needs to define a custom data structure and a few functions, that will be used by run_particle_filter:
  • a custom structure which holds the data about the model.  This will be used to dispatch the methods to be defined, listed below;
  • an initialisation function with the following signature:
    init(params_dict::Dict, n_tasks::Integer) -> model
    with params_dict a dictionary with the parameters of the model and n_tasks an integer specifying the maximum number tasks (coroutines) parallelisable operations will be scheduled over. This initialisation function should create an instance of the model data structure and return it. The value of n_tasks can be used to create  task-specific buffers for writing to when computing the model updates to avoid  reallocating data structures on each function call. As tasks may be run in parallel  over multiple threads, any buffers used in functions called within tasks should be unique to the task; to facilitate this functions in the model interface (see below) which may be called within tasks scheduled in parallel are passed a task_index  argument which is a integer index in 1:n_tasks which is guaranteed to be unique to a particular task and so can be used to index in to task specific buffers.
  • The model needs to extend the following methods, using the model type for dispatch:
ParticleDA.get_state_dimension — Function
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.
source
ParticleDA.get_observation_dimension — Function
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.
source
ParticleDA.get_state_eltype — Function
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.
source
ParticleDA.get_observation_eltype — Function
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.
source
ParticleDA.sample_initial_state! — Function
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.
source
ParticleDA.update_state_deterministic! — Function
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.
source
ParticleDA.update_state_stochastic! — Function
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.
source
ParticleDA.sample_observation_given_state! — Function
ParticleDA.sample_observation_given_state!(
+) -> 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.
source
The filter_type argument to run_particle_filter should be a concrete subtype of the ParticleFilter abstract type.
The summary_stat_type argument to run_particle_filter should be a concrete subtype of the AbstractSummaryStat abstract type.
ParticleDA.AbstractCustomReductionSummaryStat — Type
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.
source
ParticleDA.AbstractSumReductionSummaryStat — Type
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.
source
ParticleDA.MeanAndVarSummaryStat — Type
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.
source
ParticleDA.MeanSummaryStat — Type
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.
source
ParticleDA.NaiveMeanAndVarSummaryStat — Type
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.
source
ParticleDA.NaiveMeanSummaryStat — Type
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.
source
The next section details how to write the interface between the model and the particle filter.
Interfacing the model
The model needs to define a custom data structure and a few functions, that will be used by run_particle_filter:
  • a custom structure which holds the data about the model.  This will be used to dispatch the methods to be defined, listed below;
  • an initialisation function with the following signature:
    init(params_dict::Dict, n_tasks::Integer) -> model
    with params_dict a dictionary with the parameters of the model and n_tasks an integer specifying the maximum number tasks (coroutines) parallelisable operations will be scheduled over. This initialisation function should create an instance of the model data structure and return it. The value of n_tasks can be used to create  task-specific buffers for writing to when computing the model updates to avoid  reallocating data structures on each function call. As tasks may be run in parallel  over multiple threads, any buffers used in functions called within tasks should be unique to the task; to facilitate this functions in the model interface (see below) which may be called within tasks scheduled in parallel are passed a task_index  argument which is a integer index in 1:n_tasks which is guaranteed to be unique to a particular task and so can be used to index in to task specific buffers.
  • The model needs to extend the following methods, using the model type for dispatch:
ParticleDA.get_state_dimension — Function
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.
source
ParticleDA.get_observation_dimension — Function
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.
source
ParticleDA.get_state_eltype — Function
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.
source
ParticleDA.get_observation_eltype — Function
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.
source
ParticleDA.sample_initial_state! — Function
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.
source
ParticleDA.update_state_deterministic! — Function
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.
source
ParticleDA.update_state_stochastic! — Function
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.
source
ParticleDA.sample_observation_given_state! — Function
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.
source
ParticleDA.get_log_density_observation_given_state — Function
ParticleDA.get_log_density_observation_given_state(
+)
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.
source
ParticleDA.get_log_density_observation_given_state — Function
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.
source
ParticleDA.write_model_metadata — Function
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.
source
  • Optionally, if the model has additive Gaussian observation and state noise, it may also extend the following methods, again using the model type for dispatch, to allow using the more statistically efficient OptimalFilter for filtering
ParticleDA.get_observation_mean_given_state! — Function
ParticleDA.get_observation_mean_given_state!(
+) -> 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.
source
ParticleDA.write_model_metadata — Function
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.
source
  • Optionally, if the model has additive Gaussian observation and state noise, it may also extend the following methods, again using the model type for dispatch, to allow using the more statistically efficient OptimalFilter for filtering
ParticleDA.get_observation_mean_given_state! — Function
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.
source
ParticleDA.get_covariance_state_noise — Function
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.
source
ParticleDA.get_covariance_observation_noise — Function
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.
source
ParticleDA.get_covariance_state_observation_given_previous_state — Function
ParticleDA.get_covariance_state_observation_given_previous_state(
+)
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.
source
ParticleDA.get_covariance_state_noise — Function
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.
source
ParticleDA.get_covariance_observation_noise — Function
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.
source
ParticleDA.get_covariance_state_observation_given_previous_state — Function
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.
source
ParticleDA.get_covariance_observation_observation_given_previous_state — Function
ParticleDA.get_covariance_observation_observation_given_previous_state(
+) -> 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.
source
ParticleDA.get_covariance_observation_observation_given_previous_state — Function
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.
source
Other function you may want to extend
ParticleDA.get_state_indices_correlated_to_observations — Function
ParticleDA.get_state_indices_correlated_to_observations(
+) -> 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.
source
Functions are provided to run tests to check a model correctly implements the expected interfaces:
ParticleDA.run_unit_tests_for_generic_model_interface — Function
run_unit_tests_for_generic_model_interface(model, seed)
Run tests for model model correctly implementing generic model interface with random seed seed.
Tests that all required methods are defined for model and that they have expected behaviour.
source
ParticleDA.run_tests_for_optimal_proposal_model_interface — Function
run_tests_for_optimal_proposal_model_interface(
+    model, seed, estimate_bound_constant, estimate_n_samples
+)
Run tests for model model correctly implementing locally optimal proposal model interface with random seed seed, constant factor for checks of functions exhibiting expected Monte Carlo error scaling estimate_bound_constant and number of samples to use in checks of Monte Carlo error scaling estimate_n_samples.
Tests that all additional methods required to use locally optimal proposal filter are defined for model and that they have expected behaviour.
source
ParticleDA.run_tests_for_convergence_of_filter_estimates_against_kalman_filter — Function
run_tests_for_convergence_of_filter_estimates_against_kalman_filter(
+    filter_type,
+    init_model,
+    model_parameters_dict,
+    seed,
+    n_time_step,
+    n_particles,
+    mean_rmse_bound_constant,
+    log_var_rmse_bound_constant,
+)
Run tests to check convergence of estimates produced using filter with type filter_type for (linear-Gaussian) model initialised by init_model and with parameters model_parameter_dict against ground-truth values computed using a Kalman filter, running filtering for n_time_step and with a sequence of n_particles ensemble sizes and with random number generator seeded with seed. Estimates of filtering distribution means and (log) variances are checked to show the expected Monte Carlo error (1 / sqrt(n_particle)) scaling with constant factors mean_rmse_bound_constant and log_var_rmse_bound_constant used in checks - that is error < bound_constant / sqrt(n_particle).
source
Other function you may want to extend
ParticleDA.get_state_indices_correlated_to_observations — Function
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)`.
source
ParticleDA.get_initial_state_mean! — Function
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.
source
ParticleDA.write_snapshot — Function
ParticleDA.write_snapshot(
+) -> 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)`.
source
ParticleDA.get_initial_state_mean! — Function
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.
source
ParticleDA.write_snapshot — Function
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.
source
ParticleDA.write_state — Function
ParticleDA.write_state(
+)
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.
source
ParticleDA.write_state — Function
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.
source
ParticleDA.get_covariance_initial_state — Function
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.
source
ParticleDA.init_filter — Function
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.
source
ParticleDA.get_covariance_initial_state — Function
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.
source
ParticleDA.init_filter — Function
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.
source
ParticleDA.sample_proposal_and_compute_log_weights! — Function
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.
source
ParticleDA.write_observation — Function
ParticleDA.write_observation(
+)
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.
source
ParticleDA.write_observation — Function
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.
source
ParticleDA.write_weights — Function
ParticleDA.write_weights(
+)
Write the observations at time index time_index represented by the vector observation to the HDF5 file file for the model represented by model.
source
ParticleDA.write_weights — Function
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.
source
Input parameters
You can store the input parameters in an YAML file with the following structure
filter:
+)
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.
source

Input parameters

You can store the input parameters in an YAML file with the following structure

filter:
   key1: value1
 
 model:
@@ -49,7 +60,7 @@
     key3: value3
   model_name2:
     key4: value4
-    key5: value5

The parameters under filter are related to the particle filter, under model you can specify the parameters for different models.

The particle filter parameters are saved in the following data structure:

ParticleDA.FilterParameters — Type
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.
source

Example: estimating the state in a tsunami model

A full example of a model interfacing ParticleDA is available in test/models/llw2d.jl. This model represents a simple two dimensional simulation of tsunami dynamics and is partly based on the tsunami data assimilation code by Takuto Maeda. The particle filter can be run with observations simulated from the model as follows

# Load ParticleDA
+    key5: value5

The parameters under filter are related to the particle filter, under model you can specify the parameters for different models.

The particle filter parameters are saved in the following data structure:

ParticleDA.FilterParameters — Type
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.
source

Example: estimating the state in a tsunami model

A full example of a model interfacing ParticleDA is available in test/models/llw2d.jl. This model represents a simple two dimensional simulation of tsunami dynamics and is partly based on the tsunami data assimilation code by Takuto Maeda. The particle filter can be run with observations simulated from the model as follows

# Load ParticleDA
 using ParticleDA
 
 # Save some variables for later use
@@ -75,7 +86,7 @@
 # to use NaiveMeanAndVarSummaryStat instead
 final_states, final_statistics = run_particle_filter(
   LLW2d.init, input_file, observation_file, OptimalFilter, MeanAndVarSummaryStat
-)

Parameters of the test model

The Input parameters section shows how to pass the input parameters for the filter and the model. For the model included in the test suite, called llw2d, you can set the following parameters:

Main.LLW2d.LLW2dModelParameters — Type
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)
source

Observation station coordinates

The coordinates of the observation stations can be set in two different ways.

  1. Provide the coordinates in an input file. Set the parameter station_filename to the name of your input file. The input file is in plain text, the format is one row per station containing x, y - coordinates in metres. Here is a simple example with two stations

    # Comment lines starting with '#' will be ignored by the code
+)

    Parameters of the test model

    The Input parameters section shows how to pass the input parameters for the filter and the model. For the model included in the test suite, called llw2d, you can set the following parameters:

    Main.LLW2d.LLW2dModelParameters — Type
    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)
    source

    Observation station coordinates

    The coordinates of the observation stations can be set in two different ways.

    1. Provide the coordinates in an input file. Set the parameter station_filename to the name of your input file. The input file is in plain text, the format is one row per station containing x, y - coordinates in metres. Here is a simple example with two stations

      # Comment lines starting with '#' will be ignored by the code
 # This file contains two stations: at [1km, 1km] and [2km, 2km]
 1.0e3, 1.0e3
 2.0e3, 2.0e3

    2. Provide parameters for an equispaced rectilinear grid of observation stations. The values of these parameters should then be set:

      • n_stations_x : Number of observation stations in the x direction
      • n_stations_y : Number of observation stations in the y direction
      • station_distance_x : Distance between stations in the x direction (m)
      • station_distance_y : Distance between stations in the y direction (m)
      • station_boundary_x : Distance between bottom left edge of box and first station in the x direction (m)
      • station_boundary_y : Distance between bottom left edge of box and first station in the y direction (m)

      As an example, one could set

      n_stations_x=2,
@@ -83,4 +94,12 @@
 station_distance_x=1.0e3,
 station_distance_y=1.0e3,
 station_boundary_x=10.0e3,
-station_boundary_y=10.0e3,

      to generate 4 stations at [10km, 10km], [10km, 11km], [11km, 10km] and [11km, 11km].

    Output

    If the filter parameter verbose is set to true, run_particle_filter will produce an HDF5 file in the run directory. The file name is particle_da.h5 by default (this is configurable using the output_filename filter parameter). The file contains the summary statistics of the estimate state distribution (by default the mean and variance), particle weights, parameters used, and other metadata at each time step observation were assimilated. To read the output file, use the HDF5 library.

    A basic plotting tool is provided in a Jupyter notebook. This is intended as a template to build more sophisticated postprocessing tools, but can be used for some simple analysis. Set the variable timestamp in the third cell to plot different time slices from the output file. More functionality may be added as the package develops.

    Running in parallel

    The particle state update is parallelised using both MPI and threading. According to our preliminary tests both methods work well at small scale. To use the threading, set the environment variable JULIA_NUM_THREADS to the number of threads you want to use before starting Julia and then call the run_particle_filter function normally. You can check the number of threads julia has available by calling in Julia's REPL

    Threads.nthreads()

    To use the MPI parallelisation, write a Julia script that calls the run_particle_filter function for the relevant model and observations and run it in an Unix shell with

    mpirun -np <your_number_of_processes> julia <your_julia_script>

    Note that the parallel performance may vary depending on the performance of the algorithm. In general, a degeneracy of the particle weights will lead to poor load balance and parallel performance. See this issue for more details.

    Testing

    We have a basic test suite for ParticleDA.jl. You can run the tests by entering the package manager mode in Julia's REPL with ] and running the command

    test ParticleDA

    License

    The ParticleDA.jl package is licensed under the MIT "Expat" License.

+station_boundary_y=10.0e3,

to generate 4 stations at [10km, 10km], [10km, 11km], [11km, 10km] and [11km, 11km].

Output

If the filter parameter verbose is set to true, run_particle_filter will produce an HDF5 file in the run directory. The file name is particle_da.h5 by default (this is configurable using the output_filename filter parameter). The file contains the summary statistics of the estimate state distribution (by default the mean and variance), particle weights, parameters used, and other metadata at each time step observation were assimilated. To read the output file, use the HDF5 library.

A basic plotting tool is provided in a Jupyter notebook. This is intended as a template to build more sophisticated postprocessing tools, but can be used for some simple analysis. Set the variable timestamp in the third cell to plot different time slices from the output file. More functionality may be added as the package develops.

Running in parallel

The particle state update is parallelised using both MPI and threading. According to our preliminary tests both methods work well at small scale. To use the threading, set the environment variable JULIA_NUM_THREADS to the number of threads you want to use before starting Julia and then call the run_particle_filter function normally. You can check the number of threads julia has available by calling in Julia's REPL

Threads.nthreads()

To use the MPI parallelisation, write a Julia script that calls the run_particle_filter function for the relevant model and observations and run it in an Unix shell with

mpirun -np <your_number_of_processes> julia <your_julia_script>

Note that the parallel performance may vary depending on the performance of the algorithm. In general, a degeneracy of the particle weights will lead to poor load balance and parallel performance. See this issue for more details.

Kalman filters

To allow validation that the filter implementations give the expected results when applied to linear-Gaussian state space models, dense and matrix-free Kalman filter implementations are available in the ParticleDA.Kalman module

ParticleDA.Kalman.KalmanFilter — Type

Kalman filter for linear Gaussian state space models.

Explicitly constructs state transition and observation matrices. Assumes state transition matrix is time-invariant.

source
ParticleDA.Kalman.MatrixFreeKalmanFilter — Type

Kalman filter for linear Gaussian state space models using matrix-free updates.

Applies model state transition and observation functions directly to perform covariance updates. This is more memory efficient and allows for time varying state transitions but may be slower compared to using explicit matrix-matrix multiplies.

source
ParticleDA.Kalman.run_kalman_filter — Function
run_kalman_filter(
+    model, observation_sequence[, filter_type=KalmanFilter]
+)

Run Kalman filter on a linear-Gaussian state space model model with observations observation_sequence. The filter_type argument can be used to set the implementation used for the filtering updates.

source
ParticleDA.Kalman.lmult_by_observation_matrix! — Function
lmult_by_observation_matrix!(
+    output_matrix::AbstractMatrix, rhs_matrix::AbstractMatrix, model
+)

Compute output_matrix = observation_matrix * rhs_matrix.

source
ParticleDA.Kalman.pre_and_postmultiply_by_state_transition_matrix! — Function
pre_and_postmultiply_by_state_transition_matrix!(
+    state_covar::Matrix, filter::MatrixFreeKalmanFilter, time_index::Int
+)

Compute state_covar = transition_matrix * state_covar * transition_matrix'.

source
ParticleDA.Kalman.pre_and_postmultiply_by_observation_matrix! — Function
pre_and_postmultiply_by_state_transition_matrix!(
+    state_covar::Matrix, filter::MatrixFreeKalmanFilter, time_index::Int
+)

Compute observation_covar = observation_matrix * state_covar * observation_matrix'.

source
ParticleDA.Kalman.lmult_by_state_transition_matrix! — Function
lmult_by_state_transition_matrix!(matrix::AbstractMatrix, model, time_index)

Compute matrix = state_transition_matrix * matrix.

source

Testing

We have a basic test suite for ParticleDA.jl. You can run the tests by entering the package manager mode in Julia's REPL with ] and running the command

test ParticleDA

License

The ParticleDA.jl package is licensed under the MIT "Expat" License.

diff --git a/dev/objects.inv b/dev/objects.inv index dcd511e7e1c0c08d1403a1abfcccdca7a40142b2..778b67b0703e881d5e2fccaa993e2eb8be2af3b2 100644 GIT binary patch delta 921 zcmV;K17`f>28jreet%cVQsXcXz2_@5HE?x^cLshDtW*F96UHFn{~StmCR0@O8hm%xtLR zT4C$rccHZ@xuTBsUAOQk;ZY0kb5B9BUSYRp7=CqF>XP+KIf9X1Apj!}Mgcs#;{{W8 zh|T-P`TB;40beUxK&c`0aMTtm#Ejv-ZKN!{LuFOp%dFvZR`ib_7{HkqG) znrA!MM?^=1kbh4wKwKe;l0i!GR>3PF_T=|y0b_$#jVi(^?sS7C-Qj1d9CjMBtz)X} zAE>uNdn*K12&|B|LQbP-HG623@&!Q|zO959;5A`<=6g)Lm7BVRx=5&wl>XEuR~$+& zKAZHSk0eUiNQ57hZ`guS_1H*B&F&GBI}ydZZxnaLivzH#4XK!7h8K}s!qY&HRhA#} zkl|dFe}AaSGfHg4ubTrAr_r@xgsPI%+B)BLycWCnHSa`^yXc!xKm@X_+l)_V-CzT^J{U-MFoyEggP_8cS z?*3kt@pIw1$=VP5$I;W+Yfiq%bdm-7)0mGo`Iv35|GJG*XtiMm){hsrRXRhGE^y2e zT7Lmq)=gT6vBO(QWvG+T72JlrlJ5q@E?48ClyFtTRog;s#xQ@2DT&jao6Ix`Zs+q( zAx^`efxM3sJf+hnuC~ri6jo!grjstON0Y%GtJ9x(J@bO{j%L4QVdgC7m!p$=S5156 vl{U@TNUvD-XndJq)^#_Hket%ZUZrd;rz56Rzv^g?B(pzrs1_Vf&BI$AzCL(P!;3AMCC;s~m zNhy{rMV8x(?L_a*^5$_UWrciU0mYSYCF2*2{Q!@u)bh*n{7El<1z5}UT*IT(5%K;5 z~+?p7Q(lGA^Y8zYONv)>T z<5@eXlryLM`U0lX@`^{P3}720HYvHX95{yBf>EPk8~WMDlmv|>2W$YL#v5sMLvMHWNHGHY>P zTN@jqHlnSx9?`YI#L8Pt*UH~PZRM>5yR_mr_=nmU@Vf_UmObe!+*Ohd{a;st#3+&^ zg@{%$$|};Ud>2`-WfOB2>+JHPpWxZ6h_q0Jjl;v32Y=OsG_a!bu5tl%0ag^HGh{|t zXtyq&rCM(=OqDh>#GuYk&S5|g)ili1HcU4>r`sm18MM=UF<;+zlUlg>7fzAyaDXRT>N>wul+Z}-+ zA2xGHqD~#EW8A@M>MaOzVYz{mcFLYI2z}tNjQ>w>#hB4GiWWu=7&|TT`YAGw|%{S0d@m|J@TomYfm-$%Izahh7&N?)5k^lez diff --git a/dev/search_index.js b/dev/search_index.js index ce5bfa3..91b28a6 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"#ParticleDA.jl","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.jl is a Julia package to perform data assimilation with particle filters, distributed using MPI.","category":"page"},{"location":"#Installation","page":"ParticleDA.jl","title":"Installation","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To install the latest stable version of the package, open the Julia REPL, enter the package manager with ], then run the command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"add ParticleDA.jl","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"If you plan to develop the package (make changes, submit pull requests, etc), in the package manager mode run this command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"dev ParticleDA","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"This will automatically clone the repository to your local directory ~/.julia/dev/ParticleDA.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"You can exit from the package manager mode by pressing CTRL + C or, alternatively, the backspace key when there is no input in the prompt.","category":"page"},{"location":"#Usage","page":"ParticleDA.jl","title":"Usage","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"After installing the package, you can start using it in Julia's REPL with","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"using ParticleDA","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The run the particle filter you can use the function run_particle_filter:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"run_particle_filter","category":"page"},{"location":"#ParticleDA.run_particle_filter","page":"ParticleDA.jl","title":"ParticleDA.run_particle_filter","text":"run_particle_filter(\n init_model,\n input_file_path,\n observation_file_path,\n filter_type=BootstrapFilter,\n summary_stat_type=MeanAndVarSummaryStat;\n rng=Random.TaskLocalRNG()\n) -> Tuple{Matrix, Union{NamedTuple, Nothing}}\n\nRun 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.\n\nReturns 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.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To simulate observations from the model (which can be used to for example test the filtering algorithms) you can use the function simulate_observations_from_model","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"simulate_observations_from_model","category":"page"},{"location":"#ParticleDA.simulate_observations_from_model","page":"ParticleDA.jl","title":"ParticleDA.simulate_observations_from_model","text":"simulate_observations_from_model(\n init_model, input_file_path, output_file_path; rng=Random.TaskLocalRNG()\n) -> Matrix\n\nSimulate 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.\n\nThe 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.\n\nThe simulated observation sequence is returned as a matrix with columns corresponding to the observation vectors at each time step.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The filter_type argument to run_particle_filter should be a concrete subtype of the ParticleFilter abstract type.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.ParticleFilter\nBootstrapFilter\nOptimalFilter","category":"page"},{"location":"#ParticleDA.ParticleFilter","page":"ParticleDA.jl","title":"ParticleDA.ParticleFilter","text":"ParticleFilter\n\nAbstract type for particle filters. Currently implemented subtypes are:\n\nBootstrapFilter\nOptimalFilter\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.BootstrapFilter","page":"ParticleDA.jl","title":"ParticleDA.BootstrapFilter","text":"BootstrapFilter <: ParticleFilter\n\nSingleton type BootstrapFilter. This can be used as argument of run_particle_filter to select the bootstrap filter.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.OptimalFilter","page":"ParticleDA.jl","title":"ParticleDA.OptimalFilter","text":"OptimalFilter <: ParticleFilter\n\nSingleton type OptimalFilter. This can be used as argument of run_particle_filter to select the optimal proposal filter (for conditionally linear-Gaussian models).\n\n\n\n\n\n","category":"type"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The summary_stat_type argument to run_particle_filter should be a concrete subtype of the AbstractSummaryStat abstract type.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.AbstractSummaryStat\nParticleDA.AbstractCustomReductionSummaryStat\nParticleDA.AbstractSumReductionSummaryStat\nParticleDA.MeanAndVarSummaryStat\nParticleDA.MeanSummaryStat\nParticleDA.NaiveMeanAndVarSummaryStat\nParticleDA.NaiveMeanSummaryStat","category":"page"},{"location":"#ParticleDA.AbstractSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractSummaryStat","text":"AbstractSummaryStat\n\nAbstract 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.\n\nSee also: AbstractSumReductionSummaryStat, AbstractCustomReductionSummaryStat.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.AbstractCustomReductionSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractCustomReductionSummaryStat","text":"AbstractCustomReductionSummaryStat <: AbstractSummaryStat\n\nAbstract 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.AbstractSumReductionSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractSumReductionSummaryStat","text":"AbstractSumReductionSummaryStat <: AbstractSummaryStat\n\nAbstract 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.MeanAndVarSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.MeanAndVarSummaryStat","text":"MeanAndVarSummaryStat <: AbstractCustomReductionSummaryStat\n\nCustom 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.MeanSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.MeanSummaryStat","text":"MeanSummaryStat <: AbstractCustomReductionSummaryStat\n\nCustom 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.NaiveMeanAndVarSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.NaiveMeanAndVarSummaryStat","text":"NaiveMeanAndVarSummaryStat <: AbstractSumReductionSummaryStat\n\nSum 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.NaiveMeanSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.NaiveMeanSummaryStat","text":"NaiveMeanSummaryStat <: AbstractSumReductionSummaryStat\n\nSum 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.\n\n\n\n\n\n","category":"type"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The next section details how to write the interface between the model and the particle filter.","category":"page"},{"location":"#Interfacing-the-model","page":"ParticleDA.jl","title":"Interfacing the model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The model needs to define a custom data structure and a few functions, that will be used by run_particle_filter:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"a custom structure which holds the data about the model. This will be used to dispatch the methods to be defined, listed below;\nan initialisation function with the following signature:\ninit(params_dict::Dict, n_tasks::Integer) -> model\nwith params_dict a dictionary with the parameters of the model and n_tasks an integer specifying the maximum number tasks (coroutines) parallelisable operations will be scheduled over. This initialisation function should create an instance of the model data structure and return it. The value of n_tasks can be used to create task-specific buffers for writing to when computing the model updates to avoid reallocating data structures on each function call. As tasks may be run in parallel over multiple threads, any buffers used in functions called within tasks should be unique to the task; to facilitate this functions in the model interface (see below) which may be called within tasks scheduled in parallel are passed a task_index argument which is a integer index in 1:n_tasks which is guaranteed to be unique to a particular task and so can be used to index in to task specific buffers.\nThe model needs to extend the following methods, using the model type for dispatch:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_state_dimension\nParticleDA.get_observation_dimension\nParticleDA.get_state_eltype\nParticleDA.get_observation_eltype\nParticleDA.sample_initial_state!\nParticleDA.update_state_deterministic!\nParticleDA.update_state_stochastic!\nParticleDA.sample_observation_given_state!\nParticleDA.get_log_density_observation_given_state\nParticleDA.write_model_metadata","category":"page"},{"location":"#ParticleDA.get_state_dimension","page":"ParticleDA.jl","title":"ParticleDA.get_state_dimension","text":"ParticleDA.get_state_dimension(model) -> Integer\n\nReturn the positive integer dimension of the state vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_observation_dimension","page":"ParticleDA.jl","title":"ParticleDA.get_observation_dimension","text":"ParticleDA.get_observation_dimension(model) -> Integer\n\nReturn the positive integer dimension of the observation vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_state_eltype","page":"ParticleDA.jl","title":"ParticleDA.get_state_eltype","text":"ParticleDA.get_state_eltype(model) -> Type\n\nReturn the element type of the state vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_observation_eltype","page":"ParticleDA.jl","title":"ParticleDA.get_observation_eltype","text":"ParticleDA.get_observation_eltype(model) -> Type\n\nReturn the element type of the observation vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_initial_state!","page":"ParticleDA.jl","title":"ParticleDA.sample_initial_state!","text":"ParticleDA.sample_initial_state!(state, model, rng, task_index=1)\n\nSample 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.update_state_deterministic!","page":"ParticleDA.jl","title":"ParticleDA.update_state_deterministic!","text":"ParticleDA.update_state_deterministic!(state, model, time_index, task_index=1)\n\nApply 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.update_state_stochastic!","page":"ParticleDA.jl","title":"ParticleDA.update_state_stochastic!","text":"ParticleDA.update_state_stochastic!(state, model, rng, task_index=1)\n\nApply 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_observation_given_state!","page":"ParticleDA.jl","title":"ParticleDA.sample_observation_given_state!","text":"ParticleDA.sample_observation_given_state!(\n observation, state, model, rng, task_index=1\n)\n\nSimulate noisy observations of the state state of model described by model and write to observation array using rng to generate any random draws.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_log_density_observation_given_state","page":"ParticleDA.jl","title":"ParticleDA.get_log_density_observation_given_state","text":"ParticleDA.get_log_density_observation_given_state(\n observation, state, model, task_index=1\n) -> Real\n\nReturn 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_model_metadata","page":"ParticleDA.jl","title":"ParticleDA.write_model_metadata","text":"ParticleDA.write_model_metadata(file::HDF5.File, model)\n\nWrite metadata for with the model described by model to the HDF5 file file.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Optionally, if the model has additive Gaussian observation and state noise, it may also extend the following methods, again using the model type for dispatch, to allow using the more statistically efficient OptimalFilter for filtering","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_observation_mean_given_state!\nParticleDA.get_covariance_state_noise\nParticleDA.get_covariance_observation_noise\nParticleDA.get_covariance_state_observation_given_previous_state\nParticleDA.get_covariance_observation_observation_given_previous_state","category":"page"},{"location":"#ParticleDA.get_observation_mean_given_state!","page":"ParticleDA.jl","title":"ParticleDA.get_observation_mean_given_state!","text":"ParticleDA.get_observation_mean_given_state!(\n observation_mean, state, model, task_index=1\n)\n\nCompute the mean of the multivariate normal distribution on the observations given the current state and write to the first argument.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_state_noise","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_state_noise","text":"ParticleDA.get_covariance_state_noise(model, i, j) -> Real\n\nReturn covariance cov(U[i], U[j]) between components of the zero-mean Gaussian state noise vector U.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_observation_noise","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_observation_noise","text":"ParticleDA.get_covariance_observation_noise(model, i, j) -> Real\n\nReturn covariance cov(V[i], V[j]) between components of the zero-mean Gaussian observation noise vector V.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_state_observation_given_previous_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_state_observation_given_previous_state","text":"ParticleDA.get_covariance_state_observation_given_previous_state(\n model, i, j\n) -> Real\n\nReturn 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.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_observation_observation_given_previous_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_observation_observation_given_previous_state","text":"ParticleDA.get_covariance_observation_observation_given_previous_state(\n model, i, j\n) -> Real\n\nReturn 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.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#Other-function-you-may-want-to-extend","page":"ParticleDA.jl","title":"Other function you may want to extend","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_state_indices_correlated_to_observations\nParticleDA.get_initial_state_mean!\nParticleDA.write_snapshot\nParticleDA.write_state\nParticleDA.get_covariance_initial_state\nParticleDA.init_filter\nParticleDA.sample_proposal_and_compute_log_weights!\nParticleDA.write_observation\nParticleDA.write_weights\nParticleDA.get_initial_state_mean","category":"page"},{"location":"#ParticleDA.get_state_indices_correlated_to_observations","page":"ParticleDA.jl","title":"ParticleDA.get_state_indices_correlated_to_observations","text":"ParticleDA.get_state_indices_correlated_to_observations(\n model\n) -> AbstractVector{Int}\n\nReturn 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)`.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_initial_state_mean!","page":"ParticleDA.jl","title":"ParticleDA.get_initial_state_mean!","text":"ParticleDA.get_initial_state_mean!(state_mean, model)\n\nCompute the mean of the multivariate normal distribution on the initial state and write to the first argument.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_snapshot","page":"ParticleDA.jl","title":"ParticleDA.write_snapshot","text":"ParticleDA.write_snapshot(\n output_filename, model, filter_data, states, time_index, save_states\n)\n\nWrite 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_state","page":"ParticleDA.jl","title":"ParticleDA.write_state","text":"ParticleDA.write_state(\n file::HDF5.File,\n state::AbstractVector{T},\n time_index::Int,\n group_name::String,\n model\n)\n\nWrite 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_initial_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_initial_state","text":"ParticleDA.get_covariance_initial_state(model, i, j) -> Real\n\nReturn covariance cov(X[i], X[j]) between components of the initial state vector X.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.init_filter","page":"ParticleDA.jl","title":"ParticleDA.init_filter","text":"ParticleDA.init_filter(filter_params, model, nprt_per_rank, ::T) -> NamedTuple\n\nInitialise 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.\n\nNew filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_proposal_and_compute_log_weights!","page":"ParticleDA.jl","title":"ParticleDA.sample_proposal_and_compute_log_weights!","text":"ParticleDA.sample_proposal_and_compute_log_weights!(\n states, log_weights, observation, model, filter_data, ::T, rng\n)\n\nSample 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.\n\nNew filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_observation","page":"ParticleDA.jl","title":"ParticleDA.write_observation","text":"ParticleDA.write_observation(\n file::HDF5.File,\n observation::AbstractVector,\n time_index::Int,\n model\n)\n\nWrite the observations at time index time_index represented by the vector observation to the HDF5 file file for the model represented by model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_weights","page":"ParticleDA.jl","title":"ParticleDA.write_weights","text":"ParticleDA.write_weights(\n file::HDF5.File,\n weights::AbstractVector{T},\n time_index::Int,\n model\n)\n\nWrite the particle weights at time index time_index represented by the vector weights to the HDF5 file file for the model represented by model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_initial_state_mean","page":"ParticleDA.jl","title":"ParticleDA.get_initial_state_mean","text":"ParticleDA.get_initial_state_mean(model)\n\nCompute the mean of the multivariate normal distribution on the initial state.\n\n\n\n\n\n","category":"function"},{"location":"#Input-parameters","page":"ParticleDA.jl","title":"Input parameters","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"You can store the input parameters in an YAML file with the following structure","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"filter:\n key1: value1\n\nmodel:\n model_name1:\n key2: value2\n key3: value3\n model_name2:\n key4: value4\n key5: value5","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The parameters under filter are related to the particle filter, under model you can specify the parameters for different models.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The particle filter parameters are saved in the following data structure:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.FilterParameters","category":"page"},{"location":"#ParticleDA.FilterParameters","page":"ParticleDA.jl","title":"ParticleDA.FilterParameters","text":"FilterParameters()\n\nParameters for ParticleDA run. Keyword arguments:\n\nmaster_rank : Id of MPI rank that performs serial computations\nverbose::Bool : Flag to control whether to write output\noutput_filename::String : Name of output file\nnprt::Int : Number of particles for particle filter\nenable_timers::Bool : Flag to control run time measurements\nparticle_save_time_indices: Set of time indices to save particles at\nseed: Seed to initialise state of random number generator used for filtering\nn_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.\n\n\n\n\n\n","category":"type"},{"location":"#Example:-estimating-the-state-in-a-tsunami-model","page":"ParticleDA.jl","title":"Example: estimating the state in a tsunami model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"A full example of a model interfacing ParticleDA is available in test/models/llw2d.jl. This model represents a simple two dimensional simulation of tsunami dynamics and is partly based on the tsunami data assimilation code by Takuto Maeda. The particle filter can be run with observations simulated from the model as follows","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"# Load ParticleDA\nusing ParticleDA\n\n# Save some variables for later use\ntest_dir = joinpath(dirname(pathof(ParticleDA)), \"..\", \"test\")\nllw2d_src = joinpath(test_dir, \"models\", \"llw2d.jl\")\ninput_file = joinpath(test_dir, \"integration_test_1.yaml\")\nobservation_file = tempname()\n\n# Instantiate the test environment\nusing Pkg\nPkg.activate(test_dir)\nPkg.instantiate()\n\n# Include the sample model source code and load it\ninclude(llw2d_src)\nusing .LLW2d\n\n# Simulate observations from the model to use \nsimulate_observations_from_model(LLW2d.init, input_file, observation_file)\n\n# Run the (optimal proposal) particle filter with simulated observations computing the\n# mean and variance of the particle ensemble. On non-Intel architectures you may need\n# to use NaiveMeanAndVarSummaryStat instead\nfinal_states, final_statistics = run_particle_filter(\n LLW2d.init, input_file, observation_file, OptimalFilter, MeanAndVarSummaryStat\n)","category":"page"},{"location":"#Parameters-of-the-test-model","page":"ParticleDA.jl","title":"Parameters of the test model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The Input parameters section shows how to pass the input parameters for the filter and the model. For the model included in the test suite, called llw2d, you can set the following parameters:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"LLW2d.LLW2dModelParameters","category":"page"},{"location":"#Main.LLW2d.LLW2dModelParameters","page":"ParticleDA.jl","title":"Main.LLW2d.LLW2dModelParameters","text":"LLW2dModelParameters()\n\nParameters for the linear long wave two-dimensional (LLW2d) model. Keyword arguments:\n\nnx::Int : Number of grid points in the x direction\nny::Int : Number of grid points in the y direction\nx_length::AbstractFloat : Domain size in metres in the x direction\ny_length::AbstractFloat : Domain size in metres in the y direction\ndx::AbstractFloat : Distance in metres between grid points in the x direction\ndy::AbstractFloat : Distance in metrtes between grid points in the y direction\nstation_filename::String : Name of input file for station coordinates\nn_stations_x::Int : Number of observation stations in the x direction (if using regular grid)\nn_stations_y::Int : Number of observation stations in the y direction (if using regular grid)\nstation_distance_x::Float : Distance in metres between stations in the x direction (if using regular grid)\nstation_distance_y::Float : Distance in metres between stations in the y direction (if using regular grid)\nstation_boundary_x::Float : Distance in metres between bottom left edge of box and first station in the x direction (if using regular grid)\nstation_boundary_y::Float : Distance in metres between bottom left edge of box and first station in the y direction (if using regular grid)\nn_integration_step::Int : Number of sub-steps to integrate the forward model per time step\ntime_step::AbstractFloat : Time step length in seconds\npeak_position::Vector{AbstractFloat} : The `[x, y] coordinates in metres of the initial wave peak\npeak_height::AbstractFloat : The height in metres of the initial wave peak\nsource_size::AbstractFloat : Cutoff distance in metres from the peak for the initial wave\nbathymetry_setup::AbstractFloat : Bathymetry set-up\nlambda::AbstractFloat : Length scale for Matérn covariance kernel in background noise\nnu::AbstractFloat : Smoothess parameter for Matérn covariance kernel in background noise\nsigma::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in background noise\nlambda_initial_state::AbstractFloat : Length scale for Matérn covariance kernel in initial state of particles\nnu_initial_state::AbstractFloat : Smoothess parameter for Matérn covariance kernel in initial state of particles\nsigma_initial_state::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in initial state of particles\npadding::Int : Min padding for circulant embedding gaussian random field generator\nprimes::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.\nuse_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.\nabsorber_thickness_fraction::Float : Thickness of absorber for sponge absorbing boundary conditions, fraction of grid size\nboundary_damping::Float : Damping for boundaries\ncutoff_depth::Float : Shallowest water depth\nobs_noise_std::Vector: Standard deviations of noise added to observations of the true state\nobserved_state_var_indices::Vector: Vector containing the indices of the observed state variables (1: height, 2: velocity x-component, 3: velocity y-component)\n\n\n\n\n\n","category":"type"},{"location":"#Observation-station-coordinates","page":"ParticleDA.jl","title":"Observation station coordinates","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The coordinates of the observation stations can be set in two different ways.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Provide the coordinates in an input file. Set the parameter station_filename to the name of your input file. The input file is in plain text, the format is one row per station containing x, y - coordinates in metres. Here is a simple example with two stations\n# Comment lines starting with '#' will be ignored by the code\n# This file contains two stations: at [1km, 1km] and [2km, 2km]\n1.0e3, 1.0e3\n2.0e3, 2.0e3\nProvide parameters for an equispaced rectilinear grid of observation stations. The values of these parameters should then be set:\nn_stations_x : Number of observation stations in the x direction\nn_stations_y : Number of observation stations in the y direction\nstation_distance_x : Distance between stations in the x direction (m)\nstation_distance_y : Distance between stations in the y direction (m)\nstation_boundary_x : Distance between bottom left edge of box and first station in the x direction (m)\nstation_boundary_y : Distance between bottom left edge of box and first station in the y direction (m)\nAs an example, one could set\nn_stations_x=2,\nn_stations_y=2,\nstation_distance_x=1.0e3,\nstation_distance_y=1.0e3,\nstation_boundary_x=10.0e3,\nstation_boundary_y=10.0e3,\nto generate 4 stations at [10km, 10km], [10km, 11km], [11km, 10km] and [11km, 11km].","category":"page"},{"location":"#Output","page":"ParticleDA.jl","title":"Output","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"If the filter parameter verbose is set to true, run_particle_filter will produce an HDF5 file in the run directory. The file name is particle_da.h5 by default (this is configurable using the output_filename filter parameter). The file contains the summary statistics of the estimate state distribution (by default the mean and variance), particle weights, parameters used, and other metadata at each time step observation were assimilated. To read the output file, use the HDF5 library.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"A basic plotting tool is provided in a Jupyter notebook. This is intended as a template to build more sophisticated postprocessing tools, but can be used for some simple analysis. Set the variable timestamp in the third cell to plot different time slices from the output file. More functionality may be added as the package develops.","category":"page"},{"location":"#Running-in-parallel","page":"ParticleDA.jl","title":"Running in parallel","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The particle state update is parallelised using both MPI and threading. According to our preliminary tests both methods work well at small scale. To use the threading, set the environment variable JULIA_NUM_THREADS to the number of threads you want to use before starting Julia and then call the run_particle_filter function normally. You can check the number of threads julia has available by calling in Julia's REPL","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Threads.nthreads()","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To use the MPI parallelisation, write a Julia script that calls the run_particle_filter function for the relevant model and observations and run it in an Unix shell with","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"mpirun -np julia ","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Note that the parallel performance may vary depending on the performance of the algorithm. In general, a degeneracy of the particle weights will lead to poor load balance and parallel performance. See this issue for more details.","category":"page"},{"location":"#Testing","page":"ParticleDA.jl","title":"Testing","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"We have a basic test suite for ParticleDA.jl. You can run the tests by entering the package manager mode in Julia's REPL with ] and running the command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"test ParticleDA","category":"page"},{"location":"#License","page":"ParticleDA.jl","title":"License","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The ParticleDA.jl package is licensed under the MIT \"Expat\" License.","category":"page"}] +[{"location":"#ParticleDA.jl","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.jl is a Julia package to perform data assimilation with particle filters, distributed using MPI.","category":"page"},{"location":"#Installation","page":"ParticleDA.jl","title":"Installation","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To install the latest stable version of the package, open the Julia REPL, enter the package manager with ], then run the command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"add ParticleDA.jl","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"If you plan to develop the package (make changes, submit pull requests, etc), in the package manager mode run this command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"dev ParticleDA","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"This will automatically clone the repository to your local directory ~/.julia/dev/ParticleDA.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"You can exit from the package manager mode by pressing CTRL + C or, alternatively, the backspace key when there is no input in the prompt.","category":"page"},{"location":"#Usage","page":"ParticleDA.jl","title":"Usage","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"After installing the package, you can start using it in Julia's REPL with","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"using ParticleDA","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The run the particle filter you can use the function run_particle_filter:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"run_particle_filter","category":"page"},{"location":"#ParticleDA.run_particle_filter","page":"ParticleDA.jl","title":"ParticleDA.run_particle_filter","text":"run_particle_filter(\n init_model,\n input_file_path,\n observation_file_path,\n filter_type=BootstrapFilter,\n summary_stat_type=MeanAndVarSummaryStat;\n rng=Random.TaskLocalRNG()\n) -> Tuple{Matrix, Union{NamedTuple, Nothing}}\n\nRun 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.\n\nReturns 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.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To simulate observations from the model (which can be used to for example test the filtering algorithms) you can use the function simulate_observations_from_model","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"simulate_observations_from_model","category":"page"},{"location":"#ParticleDA.simulate_observations_from_model","page":"ParticleDA.jl","title":"ParticleDA.simulate_observations_from_model","text":"simulate_observations_from_model(\n init_model, input_file_path, output_file_path; rng=Random.TaskLocalRNG()\n) -> Matrix\n\nSimulate 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.\n\nThe 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.\n\nThe simulated observation sequence is returned as a matrix with columns corresponding to the observation vectors at each time step.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The filter_type argument to run_particle_filter should be a concrete subtype of the ParticleFilter abstract type.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.ParticleFilter\nBootstrapFilter\nOptimalFilter","category":"page"},{"location":"#ParticleDA.ParticleFilter","page":"ParticleDA.jl","title":"ParticleDA.ParticleFilter","text":"ParticleFilter\n\nAbstract type for particle filters. Currently implemented subtypes are:\n\nBootstrapFilter\nOptimalFilter\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.BootstrapFilter","page":"ParticleDA.jl","title":"ParticleDA.BootstrapFilter","text":"BootstrapFilter <: ParticleFilter\n\nSingleton type BootstrapFilter. This can be used as argument of run_particle_filter to select the bootstrap filter.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.OptimalFilter","page":"ParticleDA.jl","title":"ParticleDA.OptimalFilter","text":"OptimalFilter <: ParticleFilter\n\nSingleton type OptimalFilter. This can be used as argument of run_particle_filter to select the optimal proposal filter (for conditionally linear-Gaussian models).\n\n\n\n\n\n","category":"type"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The summary_stat_type argument to run_particle_filter should be a concrete subtype of the AbstractSummaryStat abstract type.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.AbstractSummaryStat\nParticleDA.AbstractCustomReductionSummaryStat\nParticleDA.AbstractSumReductionSummaryStat\nParticleDA.MeanAndVarSummaryStat\nParticleDA.MeanSummaryStat\nParticleDA.NaiveMeanAndVarSummaryStat\nParticleDA.NaiveMeanSummaryStat","category":"page"},{"location":"#ParticleDA.AbstractSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractSummaryStat","text":"AbstractSummaryStat\n\nAbstract 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.\n\nSee also: AbstractSumReductionSummaryStat, AbstractCustomReductionSummaryStat.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.AbstractCustomReductionSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractCustomReductionSummaryStat","text":"AbstractCustomReductionSummaryStat <: AbstractSummaryStat\n\nAbstract 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.AbstractSumReductionSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.AbstractSumReductionSummaryStat","text":"AbstractSumReductionSummaryStat <: AbstractSummaryStat\n\nAbstract 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.MeanAndVarSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.MeanAndVarSummaryStat","text":"MeanAndVarSummaryStat <: AbstractCustomReductionSummaryStat\n\nCustom 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.MeanSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.MeanSummaryStat","text":"MeanSummaryStat <: AbstractCustomReductionSummaryStat\n\nCustom 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.NaiveMeanAndVarSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.NaiveMeanAndVarSummaryStat","text":"NaiveMeanAndVarSummaryStat <: AbstractSumReductionSummaryStat\n\nSum 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.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.NaiveMeanSummaryStat","page":"ParticleDA.jl","title":"ParticleDA.NaiveMeanSummaryStat","text":"NaiveMeanSummaryStat <: AbstractSumReductionSummaryStat\n\nSum 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.\n\n\n\n\n\n","category":"type"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The next section details how to write the interface between the model and the particle filter.","category":"page"},{"location":"#Interfacing-the-model","page":"ParticleDA.jl","title":"Interfacing the model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The model needs to define a custom data structure and a few functions, that will be used by run_particle_filter:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"a custom structure which holds the data about the model. This will be used to dispatch the methods to be defined, listed below;\nan initialisation function with the following signature:\ninit(params_dict::Dict, n_tasks::Integer) -> model\nwith params_dict a dictionary with the parameters of the model and n_tasks an integer specifying the maximum number tasks (coroutines) parallelisable operations will be scheduled over. This initialisation function should create an instance of the model data structure and return it. The value of n_tasks can be used to create task-specific buffers for writing to when computing the model updates to avoid reallocating data structures on each function call. As tasks may be run in parallel over multiple threads, any buffers used in functions called within tasks should be unique to the task; to facilitate this functions in the model interface (see below) which may be called within tasks scheduled in parallel are passed a task_index argument which is a integer index in 1:n_tasks which is guaranteed to be unique to a particular task and so can be used to index in to task specific buffers.\nThe model needs to extend the following methods, using the model type for dispatch:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_state_dimension\nParticleDA.get_observation_dimension\nParticleDA.get_state_eltype\nParticleDA.get_observation_eltype\nParticleDA.sample_initial_state!\nParticleDA.update_state_deterministic!\nParticleDA.update_state_stochastic!\nParticleDA.sample_observation_given_state!\nParticleDA.get_log_density_observation_given_state\nParticleDA.write_model_metadata","category":"page"},{"location":"#ParticleDA.get_state_dimension","page":"ParticleDA.jl","title":"ParticleDA.get_state_dimension","text":"ParticleDA.get_state_dimension(model) -> Integer\n\nReturn the positive integer dimension of the state vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_observation_dimension","page":"ParticleDA.jl","title":"ParticleDA.get_observation_dimension","text":"ParticleDA.get_observation_dimension(model) -> Integer\n\nReturn the positive integer dimension of the observation vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_state_eltype","page":"ParticleDA.jl","title":"ParticleDA.get_state_eltype","text":"ParticleDA.get_state_eltype(model) -> Type\n\nReturn the element type of the state vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_observation_eltype","page":"ParticleDA.jl","title":"ParticleDA.get_observation_eltype","text":"ParticleDA.get_observation_eltype(model) -> Type\n\nReturn the element type of the observation vector which is assumed to be fixed for all time steps.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_initial_state!","page":"ParticleDA.jl","title":"ParticleDA.sample_initial_state!","text":"ParticleDA.sample_initial_state!(state, model, rng, task_index=1)\n\nSample 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.update_state_deterministic!","page":"ParticleDA.jl","title":"ParticleDA.update_state_deterministic!","text":"ParticleDA.update_state_deterministic!(state, model, time_index, task_index=1)\n\nApply 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.update_state_stochastic!","page":"ParticleDA.jl","title":"ParticleDA.update_state_stochastic!","text":"ParticleDA.update_state_stochastic!(state, model, rng, task_index=1)\n\nApply 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_observation_given_state!","page":"ParticleDA.jl","title":"ParticleDA.sample_observation_given_state!","text":"ParticleDA.sample_observation_given_state!(\n observation, state, model, rng, task_index=1\n)\n\nSimulate noisy observations of the state state of model described by model and write to observation array using rng to generate any random draws.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_log_density_observation_given_state","page":"ParticleDA.jl","title":"ParticleDA.get_log_density_observation_given_state","text":"ParticleDA.get_log_density_observation_given_state(\n observation, state, model, task_index=1\n) -> Real\n\nReturn 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.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_model_metadata","page":"ParticleDA.jl","title":"ParticleDA.write_model_metadata","text":"ParticleDA.write_model_metadata(file::HDF5.File, model)\n\nWrite metadata for with the model described by model to the HDF5 file file.\n\nThis method is intended to be extended by the user with the above signature, specifying the type of model.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Optionally, if the model has additive Gaussian observation and state noise, it may also extend the following methods, again using the model type for dispatch, to allow using the more statistically efficient OptimalFilter for filtering","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_observation_mean_given_state!\nParticleDA.get_covariance_state_noise\nParticleDA.get_covariance_observation_noise\nParticleDA.get_covariance_state_observation_given_previous_state\nParticleDA.get_covariance_observation_observation_given_previous_state","category":"page"},{"location":"#ParticleDA.get_observation_mean_given_state!","page":"ParticleDA.jl","title":"ParticleDA.get_observation_mean_given_state!","text":"ParticleDA.get_observation_mean_given_state!(\n observation_mean, state, model, task_index=1\n)\n\nCompute the mean of the multivariate normal distribution on the observations given the current state and write to the first argument.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_state_noise","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_state_noise","text":"ParticleDA.get_covariance_state_noise(model, i, j) -> Real\n\nReturn covariance cov(U[i], U[j]) between components of the zero-mean Gaussian state noise vector U.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_observation_noise","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_observation_noise","text":"ParticleDA.get_covariance_observation_noise(model, i, j) -> Real\n\nReturn covariance cov(V[i], V[j]) between components of the zero-mean Gaussian observation noise vector V.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_state_observation_given_previous_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_state_observation_given_previous_state","text":"ParticleDA.get_covariance_state_observation_given_previous_state(\n model, i, j\n) -> Real\n\nReturn 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.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_observation_observation_given_previous_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_observation_observation_given_previous_state","text":"ParticleDA.get_covariance_observation_observation_given_previous_state(\n model, i, j\n) -> Real\n\nReturn 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.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Functions are provided to run tests to check a model correctly implements the expected interfaces:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.run_unit_tests_for_generic_model_interface\nParticleDA.run_tests_for_optimal_proposal_model_interface\nParticleDA.run_tests_for_convergence_of_filter_estimates_against_kalman_filter","category":"page"},{"location":"#ParticleDA.run_unit_tests_for_generic_model_interface","page":"ParticleDA.jl","title":"ParticleDA.run_unit_tests_for_generic_model_interface","text":"run_unit_tests_for_generic_model_interface(model, seed)\n\nRun tests for model model correctly implementing generic model interface with random seed seed.\n\nTests that all required methods are defined for model and that they have expected behaviour.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.run_tests_for_optimal_proposal_model_interface","page":"ParticleDA.jl","title":"ParticleDA.run_tests_for_optimal_proposal_model_interface","text":"run_tests_for_optimal_proposal_model_interface(\n model, seed, estimate_bound_constant, estimate_n_samples\n)\n\nRun tests for model model correctly implementing locally optimal proposal model interface with random seed seed, constant factor for checks of functions exhibiting expected Monte Carlo error scaling estimate_bound_constant and number of samples to use in checks of Monte Carlo error scaling estimate_n_samples.\n\nTests that all additional methods required to use locally optimal proposal filter are defined for model and that they have expected behaviour.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.run_tests_for_convergence_of_filter_estimates_against_kalman_filter","page":"ParticleDA.jl","title":"ParticleDA.run_tests_for_convergence_of_filter_estimates_against_kalman_filter","text":"run_tests_for_convergence_of_filter_estimates_against_kalman_filter(\n filter_type,\n init_model,\n model_parameters_dict,\n seed,\n n_time_step,\n n_particles,\n mean_rmse_bound_constant,\n log_var_rmse_bound_constant,\n)\n\nRun tests to check convergence of estimates produced using filter with type filter_type for (linear-Gaussian) model initialised by init_model and with parameters model_parameter_dict against ground-truth values computed using a Kalman filter, running filtering for n_time_step and with a sequence of n_particles ensemble sizes and with random number generator seeded with seed. Estimates of filtering distribution means and (log) variances are checked to show the expected Monte Carlo error (1 / sqrt(n_particle)) scaling with constant factors mean_rmse_bound_constant and log_var_rmse_bound_constant used in checks - that is error < bound_constant / sqrt(n_particle).\n\n\n\n\n\n","category":"function"},{"location":"#Other-function-you-may-want-to-extend","page":"ParticleDA.jl","title":"Other function you may want to extend","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.get_state_indices_correlated_to_observations\nParticleDA.get_initial_state_mean!\nParticleDA.write_snapshot\nParticleDA.write_state\nParticleDA.get_covariance_initial_state\nParticleDA.init_filter\nParticleDA.sample_proposal_and_compute_log_weights!\nParticleDA.write_observation\nParticleDA.write_weights\nParticleDA.get_initial_state_mean","category":"page"},{"location":"#ParticleDA.get_state_indices_correlated_to_observations","page":"ParticleDA.jl","title":"ParticleDA.get_state_indices_correlated_to_observations","text":"ParticleDA.get_state_indices_correlated_to_observations(\n model\n) -> AbstractVector{Int}\n\nReturn 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)`.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_initial_state_mean!","page":"ParticleDA.jl","title":"ParticleDA.get_initial_state_mean!","text":"ParticleDA.get_initial_state_mean!(state_mean, model)\n\nCompute the mean of the multivariate normal distribution on the initial state and write to the first argument.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_snapshot","page":"ParticleDA.jl","title":"ParticleDA.write_snapshot","text":"ParticleDA.write_snapshot(\n output_filename, model, filter_data, states, time_index, save_states\n)\n\nWrite 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_state","page":"ParticleDA.jl","title":"ParticleDA.write_state","text":"ParticleDA.write_state(\n file::HDF5.File,\n state::AbstractVector{T},\n time_index::Int,\n group_name::String,\n model\n)\n\nWrite 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_covariance_initial_state","page":"ParticleDA.jl","title":"ParticleDA.get_covariance_initial_state","text":"ParticleDA.get_covariance_initial_state(model, i, j) -> Real\n\nReturn covariance cov(X[i], X[j]) between components of the initial state vector X.\n\nThis 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.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.init_filter","page":"ParticleDA.jl","title":"ParticleDA.init_filter","text":"ParticleDA.init_filter(filter_params, model, nprt_per_rank, ::T) -> NamedTuple\n\nInitialise 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.\n\nNew filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.sample_proposal_and_compute_log_weights!","page":"ParticleDA.jl","title":"ParticleDA.sample_proposal_and_compute_log_weights!","text":"ParticleDA.sample_proposal_and_compute_log_weights!(\n states, log_weights, observation, model, filter_data, ::T, rng\n)\n\nSample 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.\n\nNew filter implementations should extend this method specifying T as the appropriate singleton type for the new filter.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_observation","page":"ParticleDA.jl","title":"ParticleDA.write_observation","text":"ParticleDA.write_observation(\n file::HDF5.File,\n observation::AbstractVector,\n time_index::Int,\n model\n)\n\nWrite the observations at time index time_index represented by the vector observation to the HDF5 file file for the model represented by model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.write_weights","page":"ParticleDA.jl","title":"ParticleDA.write_weights","text":"ParticleDA.write_weights(\n file::HDF5.File,\n weights::AbstractVector{T},\n time_index::Int,\n model\n)\n\nWrite the particle weights at time index time_index represented by the vector weights to the HDF5 file file for the model represented by model.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.get_initial_state_mean","page":"ParticleDA.jl","title":"ParticleDA.get_initial_state_mean","text":"ParticleDA.get_initial_state_mean(model)\n\nCompute the mean of the multivariate normal distribution on the initial state.\n\n\n\n\n\n","category":"function"},{"location":"#Input-parameters","page":"ParticleDA.jl","title":"Input parameters","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"You can store the input parameters in an YAML file with the following structure","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"filter:\n key1: value1\n\nmodel:\n model_name1:\n key2: value2\n key3: value3\n model_name2:\n key4: value4\n key5: value5","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The parameters under filter are related to the particle filter, under model you can specify the parameters for different models.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The particle filter parameters are saved in the following data structure:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.FilterParameters","category":"page"},{"location":"#ParticleDA.FilterParameters","page":"ParticleDA.jl","title":"ParticleDA.FilterParameters","text":"FilterParameters()\n\nParameters for ParticleDA run. Keyword arguments:\n\nmaster_rank : Id of MPI rank that performs serial computations\nverbose::Bool : Flag to control whether to write output\noutput_filename::String : Name of output file\nnprt::Int : Number of particles for particle filter\nenable_timers::Bool : Flag to control run time measurements\nparticle_save_time_indices: Set of time indices to save particles at\nseed: Seed to initialise state of random number generator used for filtering\nn_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.\n\n\n\n\n\n","category":"type"},{"location":"#Example:-estimating-the-state-in-a-tsunami-model","page":"ParticleDA.jl","title":"Example: estimating the state in a tsunami model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"A full example of a model interfacing ParticleDA is available in test/models/llw2d.jl. This model represents a simple two dimensional simulation of tsunami dynamics and is partly based on the tsunami data assimilation code by Takuto Maeda. The particle filter can be run with observations simulated from the model as follows","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"# Load ParticleDA\nusing ParticleDA\n\n# Save some variables for later use\ntest_dir = joinpath(dirname(pathof(ParticleDA)), \"..\", \"test\")\nllw2d_src = joinpath(test_dir, \"models\", \"llw2d.jl\")\ninput_file = joinpath(test_dir, \"integration_test_1.yaml\")\nobservation_file = tempname()\n\n# Instantiate the test environment\nusing Pkg\nPkg.activate(test_dir)\nPkg.instantiate()\n\n# Include the sample model source code and load it\ninclude(llw2d_src)\nusing .LLW2d\n\n# Simulate observations from the model to use \nsimulate_observations_from_model(LLW2d.init, input_file, observation_file)\n\n# Run the (optimal proposal) particle filter with simulated observations computing the\n# mean and variance of the particle ensemble. On non-Intel architectures you may need\n# to use NaiveMeanAndVarSummaryStat instead\nfinal_states, final_statistics = run_particle_filter(\n LLW2d.init, input_file, observation_file, OptimalFilter, MeanAndVarSummaryStat\n)","category":"page"},{"location":"#Parameters-of-the-test-model","page":"ParticleDA.jl","title":"Parameters of the test model","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The Input parameters section shows how to pass the input parameters for the filter and the model. For the model included in the test suite, called llw2d, you can set the following parameters:","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"LLW2d.LLW2dModelParameters","category":"page"},{"location":"#Main.LLW2d.LLW2dModelParameters","page":"ParticleDA.jl","title":"Main.LLW2d.LLW2dModelParameters","text":"LLW2dModelParameters()\n\nParameters for the linear long wave two-dimensional (LLW2d) model. Keyword arguments:\n\nnx::Int : Number of grid points in the x direction\nny::Int : Number of grid points in the y direction\nx_length::AbstractFloat : Domain size in metres in the x direction\ny_length::AbstractFloat : Domain size in metres in the y direction\ndx::AbstractFloat : Distance in metres between grid points in the x direction\ndy::AbstractFloat : Distance in metrtes between grid points in the y direction\nstation_filename::String : Name of input file for station coordinates\nn_stations_x::Int : Number of observation stations in the x direction (if using regular grid)\nn_stations_y::Int : Number of observation stations in the y direction (if using regular grid)\nstation_distance_x::Float : Distance in metres between stations in the x direction (if using regular grid)\nstation_distance_y::Float : Distance in metres between stations in the y direction (if using regular grid)\nstation_boundary_x::Float : Distance in metres between bottom left edge of box and first station in the x direction (if using regular grid)\nstation_boundary_y::Float : Distance in metres between bottom left edge of box and first station in the y direction (if using regular grid)\nn_integration_step::Int : Number of sub-steps to integrate the forward model per time step\ntime_step::AbstractFloat : Time step length in seconds\npeak_position::Vector{AbstractFloat} : The `[x, y] coordinates in metres of the initial wave peak\npeak_height::AbstractFloat : The height in metres of the initial wave peak\nsource_size::AbstractFloat : Cutoff distance in metres from the peak for the initial wave\nbathymetry_setup::AbstractFloat : Bathymetry set-up\nlambda::AbstractFloat : Length scale for Matérn covariance kernel in background noise\nnu::AbstractFloat : Smoothess parameter for Matérn covariance kernel in background noise\nsigma::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in background noise\nlambda_initial_state::AbstractFloat : Length scale for Matérn covariance kernel in initial state of particles\nnu_initial_state::AbstractFloat : Smoothess parameter for Matérn covariance kernel in initial state of particles\nsigma_initial_state::AbstractFloat : Marginal standard deviation for Matérn covariance kernel in initial state of particles\npadding::Int : Min padding for circulant embedding gaussian random field generator\nprimes::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.\nuse_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.\nabsorber_thickness_fraction::Float : Thickness of absorber for sponge absorbing boundary conditions, fraction of grid size\nboundary_damping::Float : Damping for boundaries\ncutoff_depth::Float : Shallowest water depth\nobs_noise_std::Vector: Standard deviations of noise added to observations of the true state\nobserved_state_var_indices::Vector: Vector containing the indices of the observed state variables (1: height, 2: velocity x-component, 3: velocity y-component)\n\n\n\n\n\n","category":"type"},{"location":"#Observation-station-coordinates","page":"ParticleDA.jl","title":"Observation station coordinates","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The coordinates of the observation stations can be set in two different ways.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Provide the coordinates in an input file. Set the parameter station_filename to the name of your input file. The input file is in plain text, the format is one row per station containing x, y - coordinates in metres. Here is a simple example with two stations\n# Comment lines starting with '#' will be ignored by the code\n# This file contains two stations: at [1km, 1km] and [2km, 2km]\n1.0e3, 1.0e3\n2.0e3, 2.0e3\nProvide parameters for an equispaced rectilinear grid of observation stations. The values of these parameters should then be set:\nn_stations_x : Number of observation stations in the x direction\nn_stations_y : Number of observation stations in the y direction\nstation_distance_x : Distance between stations in the x direction (m)\nstation_distance_y : Distance between stations in the y direction (m)\nstation_boundary_x : Distance between bottom left edge of box and first station in the x direction (m)\nstation_boundary_y : Distance between bottom left edge of box and first station in the y direction (m)\nAs an example, one could set\nn_stations_x=2,\nn_stations_y=2,\nstation_distance_x=1.0e3,\nstation_distance_y=1.0e3,\nstation_boundary_x=10.0e3,\nstation_boundary_y=10.0e3,\nto generate 4 stations at [10km, 10km], [10km, 11km], [11km, 10km] and [11km, 11km].","category":"page"},{"location":"#Output","page":"ParticleDA.jl","title":"Output","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"If the filter parameter verbose is set to true, run_particle_filter will produce an HDF5 file in the run directory. The file name is particle_da.h5 by default (this is configurable using the output_filename filter parameter). The file contains the summary statistics of the estimate state distribution (by default the mean and variance), particle weights, parameters used, and other metadata at each time step observation were assimilated. To read the output file, use the HDF5 library.","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"A basic plotting tool is provided in a Jupyter notebook. This is intended as a template to build more sophisticated postprocessing tools, but can be used for some simple analysis. Set the variable timestamp in the third cell to plot different time slices from the output file. More functionality may be added as the package develops.","category":"page"},{"location":"#Running-in-parallel","page":"ParticleDA.jl","title":"Running in parallel","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The particle state update is parallelised using both MPI and threading. According to our preliminary tests both methods work well at small scale. To use the threading, set the environment variable JULIA_NUM_THREADS to the number of threads you want to use before starting Julia and then call the run_particle_filter function normally. You can check the number of threads julia has available by calling in Julia's REPL","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Threads.nthreads()","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To use the MPI parallelisation, write a Julia script that calls the run_particle_filter function for the relevant model and observations and run it in an Unix shell with","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"mpirun -np julia ","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"Note that the parallel performance may vary depending on the performance of the algorithm. In general, a degeneracy of the particle weights will lead to poor load balance and parallel performance. See this issue for more details.","category":"page"},{"location":"#Kalman-filters","page":"ParticleDA.jl","title":"Kalman filters","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"To allow validation that the filter implementations give the expected results when applied to linear-Gaussian state space models, dense and matrix-free Kalman filter implementations are available in the ParticleDA.Kalman module","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"ParticleDA.Kalman.KalmanFilter\nParticleDA.Kalman.MatrixFreeKalmanFilter\nParticleDA.Kalman.run_kalman_filter\nParticleDA.Kalman.lmult_by_observation_matrix!\nParticleDA.Kalman.pre_and_postmultiply_by_state_transition_matrix!\nParticleDA.Kalman.pre_and_postmultiply_by_observation_matrix!\nParticleDA.Kalman.lmult_by_state_transition_matrix!","category":"page"},{"location":"#ParticleDA.Kalman.KalmanFilter","page":"ParticleDA.jl","title":"ParticleDA.Kalman.KalmanFilter","text":"Kalman filter for linear Gaussian state space models.\n\nExplicitly constructs state transition and observation matrices. Assumes state transition matrix is time-invariant.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.Kalman.MatrixFreeKalmanFilter","page":"ParticleDA.jl","title":"ParticleDA.Kalman.MatrixFreeKalmanFilter","text":"Kalman filter for linear Gaussian state space models using matrix-free updates.\n\nApplies model state transition and observation functions directly to perform covariance updates. This is more memory efficient and allows for time varying state transitions but may be slower compared to using explicit matrix-matrix multiplies.\n\n\n\n\n\n","category":"type"},{"location":"#ParticleDA.Kalman.run_kalman_filter","page":"ParticleDA.jl","title":"ParticleDA.Kalman.run_kalman_filter","text":"run_kalman_filter(\n model, observation_sequence[, filter_type=KalmanFilter]\n)\n\nRun Kalman filter on a linear-Gaussian state space model model with observations observation_sequence. The filter_type argument can be used to set the implementation used for the filtering updates.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.Kalman.lmult_by_observation_matrix!","page":"ParticleDA.jl","title":"ParticleDA.Kalman.lmult_by_observation_matrix!","text":"lmult_by_observation_matrix!(\n output_matrix::AbstractMatrix, rhs_matrix::AbstractMatrix, model\n)\n\nCompute output_matrix = observation_matrix * rhs_matrix.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.Kalman.pre_and_postmultiply_by_state_transition_matrix!","page":"ParticleDA.jl","title":"ParticleDA.Kalman.pre_and_postmultiply_by_state_transition_matrix!","text":"pre_and_postmultiply_by_state_transition_matrix!(\n state_covar::Matrix, filter::MatrixFreeKalmanFilter, time_index::Int\n)\n\nCompute state_covar = transition_matrix * state_covar * transition_matrix'.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.Kalman.pre_and_postmultiply_by_observation_matrix!","page":"ParticleDA.jl","title":"ParticleDA.Kalman.pre_and_postmultiply_by_observation_matrix!","text":"pre_and_postmultiply_by_state_transition_matrix!(\n state_covar::Matrix, filter::MatrixFreeKalmanFilter, time_index::Int\n)\n\nCompute observation_covar = observation_matrix * state_covar * observation_matrix'.\n\n\n\n\n\n","category":"function"},{"location":"#ParticleDA.Kalman.lmult_by_state_transition_matrix!","page":"ParticleDA.jl","title":"ParticleDA.Kalman.lmult_by_state_transition_matrix!","text":"lmult_by_state_transition_matrix!(matrix::AbstractMatrix, model, time_index)\n\nCompute matrix = state_transition_matrix * matrix.\n\n\n\n\n\n","category":"function"},{"location":"#Testing","page":"ParticleDA.jl","title":"Testing","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"We have a basic test suite for ParticleDA.jl. You can run the tests by entering the package manager mode in Julia's REPL with ] and running the command","category":"page"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"test ParticleDA","category":"page"},{"location":"#License","page":"ParticleDA.jl","title":"License","text":"","category":"section"},{"location":"","page":"ParticleDA.jl","title":"ParticleDA.jl","text":"The ParticleDA.jl package is licensed under the MIT \"Expat\" License.","category":"page"}] }