Package 'dust2'

Description

Control browser-based debugging of dust models. This help page documents three functions that can be used to control if and how the browser is enabled. You can't enter the debugger from any of these functions; it is only enabled if present in your C++ code (or if using odin2 if you have enabled it).

Usage

dust_browser_enabled(value = TRUE)

dust_browser_verbosity(level)

dust_browser_continue()

Arguments

Details

dust2 includes an extremely simple debugging system, and if you are reading this message, there's a good chance you are inside it. It is built on top of R's browser() and so all the usual tips, tricks and issues for working with this apply. We recommend setting the R option browserNLdisabled = TRUE to avoid surprises from pressing ⁠<enter>⁠.

• You can press n or c to proceed to the next enabled iteration

• You can press Q to quit the browser (this will end up as an error by the time you have control back)

These commands are established by browser, and can't be disabled. This means that if you have an variable called n you will need to work with it as (n) (i.e. in parentheses). This applies to all of browser's command variables (c, f, n, s, r and Q); please see browser for more information.

By the time the environment has been created, some variables from your model will have been copied into the environment; you can see these by running ls() and write expressions involving these objects. Changes that you make in R are not (currently) propagated back into the running system.

If you enable the debugger, you may have very many iterations to get through before control is returned back to the console. You can run dust_debug_continue() to prevent entry into the debugger until control is passed back to you; this means the time series will run to completion and then the next time you run the system the debugger will be triggered again. Alternatively, you can run dust_debug_enabled(FALSE) to disable all calls to the debugger.

Value

Both dust_browser_enabled and dust_browser_verbosity return the previous value of the option they are setting.

Description

Compile a dust2 system from a C++ input file. This function will compile the dust support around your system and return an object that you can call with no arguments to make a dust_system_generator object, suitable for using with dust functions (starting from dust_system_create).

Usage

dust_compile(
  filename,
  quiet = NULL,
  workdir = NULL,
  linking_to = NULL,
  cpp_std = NULL,
  compiler_options = NULL,
  optimisation_level = NULL,
  debug = NULL,
  skip_cache = FALSE
)

Arguments

Value

A function, which can be called with no arguments to yield a dust_system_generator function.

Description

Load example generators from dust2. These generators exist primarily for the examples and documentation and are not (yet) very interesting. The examples will likely change as the package evolves and some may be removed.

Usage

dust_example(name)

Arguments

Details

All models exist as source code in the package; to view the sir model you could write:

file.show(system.file("examples/sir.cpp", package = "dust2"))

Value

A dust_generator object, which you might pass into dust_system_create

sir

A simple SIR (Susceptible-Infected-Recovered) compartmental model. This model has parameters:

• N: total population size

• I0: initial infected population size (when using dust_system_set_state_initial)

• beta: per-contact rate of infection

• gamma: rate of recovery

• exp_noise: noise parameter used in the comparison to data

The system will have compartments S, I, R, cases_cumul and cases_inc

sirode

The same model as sir but in continuous time, deterministically

walk

A random walk in discrete time with Gaussian increments. This model has parameters:

• sd: The standard deviation of the Gaussian update (per unit time)

• len: The number of independent walks

• random_initial': Boolean, indicating if the initial position should be random (changes how dust_system_set_state_initial would initialise the system)

Examples

walk <- dust_example("walk")
walk

sys <- dust_system_create(walk, list(sd = 1), 20)
y <- dust_system_simulate(sys, 0:50)
matplot(t(y[1, , ]), type = "l", col = "#0000ff55", lty = 1,
        xlab = "Time", ylab = "Value")

Description

Create a particle filter object

Usage

dust_filter_create(
  generator,
  time_start,
  data,
  n_particles,
  n_groups = NULL,
  dt = NULL,
  ode_control = NULL,
  shared_data = FALSE,
  n_threads = 1,
  preserve_group_dimension = FALSE,
  seed = NULL
)

Arguments

Value

A dust_likelihood object, which can be used with dust_likelihood_run

Description

Prepare data for use with dust_unfilter_create or dust_filter_create. You do not have to use this function if you name your data.frame with our standard column names (i.e., time column containing the time) as it will be called within the filter functions directly. However, you can use this to validate your data separately or to use different columns than the defaults.

Usage

dust_filter_data(data, time = NULL, group = NULL)

Arguments

Value

A data.frame, with the addition of the class attribute dust_filter_data; once created you should not modify this object.

Description

Create an independent copy of a likelihood object. The new object is decoupled from the random number streams of the parent object. It is also decoupled from the state size of the parent object, so you can use this to create a new object where the system is fundamentally different but everything else is the same.

Usage

dust_likelihood_copy(obj, seed = NULL)

Arguments

Value

A new dust_likelihood object

Description

Fetch the last gradient created by running an likelihood. This errors if the last call to dust_likelihood_run did not use adjoint = TRUE. The first time you call this (after a particular set of parameters) it will trigger running the reverse model.

Usage

dust_likelihood_last_gradient(obj)

Arguments

Value

A vector (if ungrouped) or a matrix (if grouped).

Description

Get the last state from a likelihood.

Usage

dust_likelihood_last_state(obj, select_random_particle = FALSE)

Arguments

Value

An array. If ungrouped this will have dimensions state x particle, and if grouped then state x particle x group. If select_random_particle = TRUE, the second (particle) dimension will be dropped. This is the same as the state returned by dust_likelihood_last_trajectories without the time dimension but also without any state index applied (i.e., we always return all state).

Description

Fetch the last trajectories created by running a likelihood. This errors if the last call to dust_likelihood_run did not use save_trajectories = TRUE. We return the states and groups that were run via the index_state and index_group arguments to dust_likelihood_run.

Usage

dust_likelihood_last_trajectories(obj, select_random_particle = FALSE)

Arguments

Value

An array. If ungrouped this will have dimensions state x particle x time, and if grouped then state x particle x group x time. If select_random_particle = TRUE, the second (particle) dimension will be dropped.

Description

Create a monty_model from a dust_likelihood object.

Usage

dust_likelihood_monty(
  obj,
  packer,
  initial = NULL,
  domain = NULL,
  failure_is_impossible = FALSE,
  save_state = FALSE,
  save_trajectories = FALSE
)

Arguments

Value

A monty::monty_model object

Random number streams

This section is only relevant if your likelihood object is a particle filter, and therefore uses random numbers.

The short version: the seed argument that you may have passed to dust_filter_create will be ignored when using a dust_likelihood_monty model with algorithms from monty. You should generally not worry about this, it is expected.

When you initialise a filter, you provide a random number seed argument. Your filter will use n_groups * (n_particles + 1) streams (one for the filter for each group, then for each group one per particle). If you run the filter directly (e.g., with dust_likelihood_run then you will advance the state of the filter). However, if you use the filter with monty (which is why you're using dust_likelihood_monty) we will ignore this seeding.

When running MCMC with n_chains chains, we need n_chains * n_groups * (n_particles + 1) random number streams - that is enough streams for every chain to have a filter with its own set of chains. monty will look after this for us, but the upshot is that the random number state that you may have previously set when building the filter will be ignored as we need to create a series of suitable seeds.

The seeds provided by monty will start at some point in the RNG state space (2^256 possible states by default). In an MCMC, each chain will have a seed that is created by performing a "long jump", moving 2^192 steps along the chain. Then within each chain we will take a series of "jumps" (2^128 steps) for each of the streams we need across the groups, filters and particles. This ensures independence across the stochastic components of the system but also the reproducibility and predictability of the system. The initial seeding performed by monty will respond to R's RNG (i.e., it will follow set.seed) if an explicit seed is not given.

Description

Get random number generator (RNG) state from the particle filter.

Usage

dust_likelihood_rng_state(obj)

dust_likelihood_set_rng_state(obj, rng_state)

Arguments

Value

A raw vector, this could be quite long. Later we will describe how you might reseed a filter or system with this state.

Description

Compute a log likelihood based on a dynamical model, either from particle filter (created with dust_filter_create) or a deterministic model (created with dust_unfilter_create).

Usage

dust_likelihood_run(
  obj,
  pars,
  initial = NULL,
  save_trajectories = FALSE,
  adjoint = NULL,
  index_state = NULL,
  index_group = NULL
)

Arguments

Value

A vector of likelihood values, with as many elements as there are groups.

Description

Create a control object for controlling the adaptive stepper for systems of ordinary differential equations (ODEs). The returned object can be passed into a continuous-time dust model on initialisation.

Usage

dust_ode_control(
  max_steps = 10000,
  atol = 1e-06,
  rtol = 1e-06,
  step_size_min = 0,
  step_size_max = Inf,
  critical_times = NULL,
  debug_record_step_times = FALSE,
  save_history = FALSE
)

Arguments

Value

A named list of class "dust_ode_control". Do not modify this after creation.

Description

Return information about OpenMP support for this machine.

Usage

dust_openmp_support(check_compile = FALSE)

Arguments

Value

A list with information about the OpenMP support on your machine.

• The first few elements come from the OpenMP library directly: num_proc, max_threads, thread_limit; these correspond to a call to the function ⁠omp_get_<name>()⁠ in C and openmp_version which is the value of the ⁠_OPENMP⁠ macro.

• A logical has_openmp which is TRUE if it looks like runtime OpenMP support is available

• The next elements tell you about different sources that might control the number of threads allowed to run: mc.cores (from the R option with the same name), OMP_THREAD_LIMIT, OMP_NUM_THREADS, MC_CORES (from environment variables), limit_r (limit computed against R-related control variables), limit_openmp (limit computed against OpenMP-related variables) and limit the smaller of limit_r and limit_openmp

• Finally, if you specified check_compile = TRUE, the logical has_openmp_compiler will indicate if it looks like we can compile with OpenMP.

See Also

dust_openmp_threads() for setting a polite number of threads.

Examples

dust_openmp_support()

Description

Politely select a number of threads to use. See Details for the algorithm used.

Usage

dust_openmp_threads(n = NULL, action = "error")

Arguments

Details

There are two limits and we will take the smaller of the two.

The first limit comes from piggy-backing off of R's normal parallel configuration; we will use the MC_CORES environment variable and mc.cores option as a guide to how many cores you are happy to use. We take mc.cores first, then MC_CORES, which is the same behaviour as parallel::mclapply and friends.

The second limit comes from OpenMP. If you do not have OpenMP support, then we use one thread (higher numbers have no effect at all in this case). If you do have OpenMP support, we take the smallest of the number of "processors" (reported by omp_get_num_procs()) the "max threads" (reported by omp_get_max_threads() and "thread_limit" (reported by omp_get_thread_limit().

See dust_openmp_support() for the values of all the values that go into this calculation.

Value

An integer, indicating the number of threads that you can use

Examples

# Default number of threads; tries to pick something friendly,
# erring on the conservative side.
dust_openmp_threads(NULL)

# Try to pick something silly and it will be reduced for you
dust_openmp_threads(1000, action = "fix")

Description

Creates or updates the generated code for a set of dust systems in a package. The user-provided code is assumed to be in inst/dust as a series of C++ files; a file inst/dust/x.cpp will be transformed into a file src/x.cpp.

Usage

dust_package(path, quiet = NULL)

Arguments

Details

Classes used within a package must be distinct; typically these will match the filenames.

We add "cpp11 attributes" to the created functions, and will run cpp11::cpp_register() on them once the generated code has been created.

Your package needs a src/Makevars file to enable OpenMP (if your system supports it). If it is not present then a suitable Makevars will be written, containing

PKG_CXXFLAGS=$(SHLIB_OPENMP_CXXFLAGS)
PKG_LIBS=$(SHLIB_OPENMP_CXXFLAGS)

following "Writing R Extensions" (see section "OpenMP support"). If your package does contain a src/Makevars file we do not attempt to edit it but will error if it looks like it does not contain these lines or similar.

You also need to make sure that your package loads the dynamic library; if you are using roxygen, then you might create a file (say, R/zzz.R) containing

#' @useDynLib packagename, .registration = TRUE
NULL

substituting packagename for your package name as appropriate. This will create an entry in NAMESPACE.

Value

Nothing, this function is called for its side effects

Description

Compare current system state against data. This is only supported for systems that have 'compare_data' support (i.e., the system definition includes a compare_data method). The current state in the system (dust_system_state) is compared against the data provided as data.

Usage

dust_system_compare_data(sys, data)

Arguments

Value

A numeric vector with as many elements as your system has groups, corresponding to the log-likelihood of the data for each group.

Description

Create a dust system object from a system generator. This allocates a system and sets an initial set of parameters. Once created you can use other dust functions to interact with it.

Usage

dust_system_create(
  generator,
  pars = list(),
  n_particles = 1,
  n_groups = 1,
  time = 0,
  dt = NULL,
  ode_control = NULL,
  seed = NULL,
  deterministic = FALSE,
  n_threads = 1,
  preserve_particle_dimension = FALSE,
  preserve_group_dimension = FALSE
)

Arguments

Value

A dust_system object, with opaque format.

Parallelisation

Many calculations within a dust system can be parallelised straightforwardly - the most important of these is typically running the model (via dust_system_run_to_time or dust_system_simulate) but we also parallelise dust_system_set_state_initial, dust_system_compare_data and even dust_system_reorder. You need to set the number of threads for parallelism at system creation, and this number cannot be usefully larger than n_particles (or n_particles * n_groups if you have a grouped system).

Description

Return internal data from the system. This is intended for debugging only, and all formats are subject to change.

Usage

dust_system_internals(
  sys,
  include_coefficients = FALSE,
  include_history = FALSE
)

Arguments

Value

If sys is a discrete-time system, this function returns NULL, as no internal data is stored. Otherwise, for a continuous-time system we return a data.frame of statistics with one row per particle. Most of the columns are simple integers or numeric values, but dydt (the current derivative of the target function with respect to time) and step_times (times that the solver has stopped at, if debug_record_step_times is in dust_ode_control was set to TRUE) will be a list of columns, each element of which is a numeric vector. If include_coefficients is TRUE, the coefficients column exists and holds a list of coefficients (the structure of these may change over time, too).

Description

Reorder states within a system. This function is primarily used for debugging and may be removed from the interface if it is not generally useful.

Usage

dust_system_reorder(sys, index)

Arguments

Value

Nothing, called for side effects only.

Description

Fetch, and set, the random number generator (RNG) state from the system.

Usage

dust_system_rng_state(sys)

dust_system_set_rng_state(sys, rng_state)

Arguments

Value

A raw vector, this could be quite long.

See Also

You can pass the state you get back from this function as the seed object to dust_system_create and dust_system_set_rng_state

Description

Run a system, advancing time and the state by repeatedly running its update method. You can advance a system up to a time (which must be in the future).

Usage

dust_system_run_to_time(sys, time)

Arguments

Value

Nothing, called for side effects only

Description

Set system state. Takes a multidimensional array (2- or 3d depending on if the system is grouped or not). Dimensions of length 1 will be recycled as appropriate. For continuous time systems, we will initialise the solver immediately after setting state, which may cause errors if your initial state is invalid for your system. There are many ways that you can use this function to set different fractions of state (a subset of states, particles or parameter groups, recycling over any dimensions that are missing). Please see the Examples section for usage.

Usage

dust_system_set_state(
  sys,
  state,
  index_state = NULL,
  index_particle = NULL,
  index_group = NULL
)

Arguments

Value

Nothing, called for side effects only

Examples

# Consider a system with 3 particles and 1 group:
sir <- dust_example("sir")
sys <- dust_system_create(sir(), list(), n_particles = 3)
# The state for this system is packed as S, I, R, cases_cumul, cases_inc:
dust_unpack_index(sys)

# Set all particles to the same state:
dust_system_set_state(sys, c(1000, 10, 0, 0, 0))
dust_system_state(sys)

# We can set everything to different states by passing a vector
# with this shape:
m <- cbind(c(1000, 10, 0, 0, 0), c(999, 11, 0, 0, 0), c(998, 12, 0, 0, 0))
dust_system_set_state(sys, m)
dust_system_state(sys)

# Or set the state for just one state:
dust_system_set_state(sys, 1, index_state = 4)
dust_system_state(sys)

# If you want to set a different state across particles, you must
# provide a *matrix* (a vector always sets the same state into
# every particle)
dust_system_set_state(sys, rbind(c(1, 2, 3)), index_state = 4)
dust_system_state(sys)

# This will not work as it can be ambiguous what you are
# trying to do:
#> dust_system_set_state(sys, c(1, 2, 3), index_state = 4)

# State can be set for specific particles:
dust_system_set_state(sys, c(900, 100, 0, 0, 0), index_particle = 2)
dust_system_state(sys)

# And you can combine 'index_particle' with 'index_state' to set
# small rectangles of state:
dust_system_set_state(sys, matrix(c(1, 2, 3, 4), 2, 2),
                      index_particle = 2:3, index_state = 4:5)
dust_system_state(sys)

Description

Set system state from a system's initial conditions. This may depend on the current time.

Usage

dust_system_set_state_initial(sys)

Arguments

Value

Nothing, called for side effects only

Description

Set time into the system. This updates the time to the provided value but does not affect the state. You may want to call dust_system_set_state or dust_system_set_state_initial after calling this.

Usage

dust_system_set_time(sys, time)

Arguments

Value

Nothing, called for side effects only

Description

Simulate a system over a series of times, returning an array of output. This output can be quite large, so you may filter states according to some index.

Usage

dust_system_simulate(sys, times, index_state = NULL)

Arguments

Value

An array with 3 dimensions (state x particle x time) or 4 dimensions (state x particle x group x time) for a grouped system.

Description

Extract system state

Usage

dust_system_state(
  sys,
  index_state = NULL,
  index_particle = NULL,
  index_group = NULL
)

Arguments

Value

An array of system state. If your system is ungrouped (i.e., n_groups = 1 and preserve_group_dimension = FALSE), then this has two dimensions (state, particle). If grouped, this has three dimensions (state, particle, group)

See Also

dust_system_set_state() for setting state and dust_system_set_state_initial() for setting state to the system-specific initial conditions.

Description

Fetch the current time from the system.

Usage

dust_system_time(sys)

Arguments

Value

A single numeric value

See Also

dust_system_set_time

Description

Update parameters used by the system. This can be used to update a subset of parameters that do not change the extent of the system, and will be potentially faster than creating a new system object.

Usage

dust_system_update_pars(sys, pars)

Arguments

Value

Nothing, called for side effects only

Description

Create an "unfilter" object, which can be used to compute a deterministic likelihood following the same algorithm as the particle filter, but limited to a single particle. The name for this method will change in future.

Usage

dust_unfilter_create(
  generator,
  time_start,
  data,
  n_particles = 1,
  n_groups = NULL,
  dt = NULL,
  ode_control = NULL,
  shared_data = FALSE,
  n_threads = 1,
  preserve_particle_dimension = FALSE,
  preserve_group_dimension = FALSE
)

Arguments

Value

A dust_likelihood object, which can be used with dust_likelihood_run

Description

Unpack state. You will see state come out of dust2 systems in several places, for example dust_system_state, but it will typically be an unstructured vector with no names; this is not very useful! However, your model knows what each element, or group of elements "means". You can unpack your state from this unstructured array into a named list using this function. The same idea applies to the higher-dimensional arrays that you might get if your system has multiple particles, multiple parameter groups or has been run for multiple time steps.

Usage

dust_unpack_state(obj, state)

dust_unpack_index(obj)

Arguments

Value

A named list, where each element corresponds to a logical compartment.

See Also

monty::monty_packer, and within that especially documentation for ⁠$unpack()⁠, which powers this function.

Examples

sir <- dust_example("sir")
sys <- dust_system_create(sir, list(), n_particles = 10, dt = 0.25)
dust_system_set_state_initial(sys)
t <- seq(0, 100, by = 5)
y <- dust_system_simulate(sys, t)
# The result here is a 5 x 10 x 21 matrix: 5 states by 10 particles by
# 21 times.
dim(y)

# The 10 particles and 21 times (following t) are simple enough, but
# what are our 5 compartments?

# You can use dust_unpack_state() to reshape your output as a
# list:
dust_unpack_state(sys, y)

# Here, the list is named following the compartments (S, I, R,
# etc) and is a 10 x 21 matrix (i.e., the remaining dimensions
# from y)

# We could apply this to the final state, which converts a 5 x 10
# matrix of state into a 5 element list of vectors, each with
# length 10:
s <- dust_system_state(sys)
dim(s)
dust_unpack_state(sys, s)

# If you need more control, you can use 'dust_unpack_index' to map
# names to positions within the state dimension of this array
dust_unpack_index(sys)