dust: Create a dust model from a C++ input file
In mrc-ide/dust: Iterate Multiple Realisations of Stochastic Models

View source: R/interface.R

dust	R Documentation

Create a dust model from a C++ input file

Description

Create a dust model from a C++ input file. This function will compile the dust support around your model and return an object that can be used to work with the model (see the Details below, and dust_generator).

Usage

dust(
  filename,
  quiet = FALSE,
  workdir = NULL,
  gpu = FALSE,
  real_type = NULL,
  linking_to = NULL,
  cpp_std = NULL,
  compiler_options = NULL,
  optimisation_level = NULL,
  skip_cache = FALSE
)

Arguments

`filename`	The path to a single C++ file
`quiet`	Logical, indicating if compilation messages from `pkgbuild` should be displayed. Error messages will be displayed on compilation failure regardless of the value used.
`workdir`	Optional working directory to use. If `NULL` uses a temporary directory. By using a different directory of your choosing you can see the generated code.
`gpu`	Logical, indicating if we should generate GPU code. This requires a considerable amount of additional software installed (CUDA toolkit and drivers) as well as a CUDA-compatible GPU. If `TRUE`, then we call dust_cuda_options with no arguments. Alternatively, call that function and pass the value here (e.g, `gpu = dust::dust_cuda_options(debug = TRUE)`). Note that due to the use of the `⁠__syncwarp()⁠` primitive this may require a GPU with compute version 70 or higher.
`real_type`	Optionally, a string indicating a substitute type to swap in for your model's `real_type` declaration. If given, then we replace the string `⁠using real_type = (double\|float)⁠` with the given type. This is primarily intended to be used as `⁠gpu = TRUE, real_type = "float"⁠` in order to create model for the GPU that will use 32 bit `floats` (rather than 64 bit doubles, which are much slower). For CPU models decreasing precision of your real type will typically just decrease precision for no additional performance.
`linking_to`	Optionally, a character vector of additional packages to add to the `DESCRIPTION`'s `LinkingTo` field. Use this when your model pulls in C++ code that is packaged within another package's header-only library.
`cpp_std`	The C++ standard to use, if you need to set one explicitly. See the section "Using C++ code" in "Writing R extensions" for the details of this, and how it interacts with the R version currently being used. For R 4.0.0 and above, C++11 will be used; as dust depends on at least this version of R you will never need to specify a version this low. Sensible options are `C++14`, `C++17`, etc, depending on the features you need and what your compiler supports.
`compiler_options`	A character vector of additional options to pass through to the C++ compiler. These will be passed through without any shell quoting or validation, so check the generated commands and outputs carefully in case of error. Note that R will apply these before anything in your personal `Makevars`.
`optimisation_level`	A shorthand way of specifying common compiler options that control optimisation level. By default (`NULL`) no options are generated from this, and the optimisation level will depend on your user `Makevars` file. Valid options are `none` which disables optimisation (`-O0`), which will be faster to compile but much slower, `standard` which enables standard level of optimisation (`-O2`), useful if your Makevars/pkgload configuration is disabling optimisation, or `max` (`-O3` and `--fast-math`) which enables some slower-to-compile and potentially unsafe optimisations. These options are applied after `compiler_options` and may override options provided there. Note that as for `compiler_options`, R will apply these before anything in your personal `Makevars`
`skip_cache`	Logical, indicating if the cache of previously compiled models should be skipped. If `TRUE` then your model will not be looked for in the cache, nor will it be added to the cache after compilation.

Value

A dust_generator object based on your source files

Input requirements

Your input dust model must satisfy a few requirements.

Define some class that implements your model (below model is assumed to be the class name)
That class must define a type internal_type (so model::internal_type) that contains its internal data that the model may change during execution (i.e., that is not shared between particles). If no such data is needed, you can do ⁠using internal_type = dust::no_internal;⁠ to indicate this.
We also need a type shared_type that contains constant internal data is shared between particles (e.g., dimensions, arrays that are read but not written). If no such data is needed, you can do ⁠using share_type = dust::no_shared;⁠ to indicate this.
That class must also include a type alias that describes the model's floating point type, real_type. Most models can include ⁠using real_type = double;⁠ in their public section.
The class must also include a type alias that describes the model's data type. If your model does not support data, then write ⁠using data_type = dust::no_data;⁠, which disables the compare_data and set_data methods. Otherwise see vignette("data") for more information.
The class must have a constructor that accepts ⁠const dust::pars_type<model>& pars⁠ for your type model. This will have elements shared and internal which you can assign into your model if needed.
The model must have a method size() returning size_t which returns the size of the system. This size may depend on values in your initialisation object but is constant within a model run.
The model must have a method initial (which may not be const), taking a time step number (size_t) and returning a ⁠std::vector<real_type>⁠ of initial state for the model.
The model must have a method update taking arguments:
- ⁠size_t time⁠: the time step number
- ⁠const double * state⁠: the state at the beginning of the time step
- ⁠dust::rng_state_type<real_type>& rng_state⁠: the dust random number generator state - this must be a reference, as it will be modified as random numbers are drawn
- double *state_next: the end state of the model (to be written to by your function)

Your update function is the core here and should update the state of the system - you are expected to update all values of state on return.

It is very important that none of the functions in the class use the R API in any way as these functions will be called in parallel.

You must also provide a data/parameter-wrangling function for producing an object of type ⁠dust::pars_type<model>⁠ from an R list. We use cpp11 for this. Your function will look like:

namespace dust {
template <>
dust::pars_type<model> dust_pars<model>(cpp11::list pars) {
  // ...
  return dust::pars_type<model>(shared, internal);
}
}

With the body interacting with pars to create an object of type model::shared_type and model::internal_type before returning the dust::pars_type object. This function will be called in serial and may use anything in the cpp11 API. All elements of the returned object must be standard C/C++ (e.g., STL) types and not cpp11/R types. If your model uses only shared or internal, you may use the single-argument constructor overload to dust::pars_type which is equivalent to using dust::no_shared or dust::no_internal for the missing argument.

Your model may provided a template specialisation ⁠dust::dust_info<model>()⁠ returning a cpp11::sexp for returning arbitrary information back to the R session:

namespace dust {
template <>
cpp11::sexp dust_info<model>(const dust::pars_type<sir>& pars) {
  return cpp11::wrap(...);
}
}

What you do with this is up to you. If not present then the info() method on the created object will return NULL. Potential use cases for this are to return information about variable ordering, or any processing done while accepting the pars object used to create the pars fed into the particles.

Configuring your model

You can optionally use C++ pseudo-attributes to configure the generated code. Currently we support two attributes:

⁠[[dust::class(classname)]]⁠ will tell dust the name of your target C++ class (in this example classname). You will need to use this if your file uses more than a single class, as otherwise will try to detect this using extremely simple heuristics.
⁠[[dust::name(modelname)]]⁠ will tell dust the name to use for the class in R code. For technical reasons this must be alphanumeric characters only (sorry, no underscore) and must not start with a number. If not included then the C++ type name will be used (either specified with ⁠[[dust::class()]]⁠ or detected).

Error handling

Your model should only throw exceptions as a last resort. One such last resort exists already if rbinom is given invalid inputs to prevent an infinite loop. If an error is thrown, all particles will complete their current run, and then the error will be rethrown - this is required by our parallel processing design. Once this happens though the state of the system is "inconsistent" as it contains particles that have run for different lengths of time. You can extract the state of the system at the point of failure (which may help with debugging) but you will be unable to continue running the object until either you reset it (with ⁠$update_state()⁠). An error will be thrown otherwise.

Things are worse on a GPU; if an error is thrown by the RNG code (happens in rbinom when given impossible inputs such as negative sizes, probabilities less than zero or greater than 1) then we currently use CUDA's ⁠__trap()⁠ function which will require a process restart to be able to use anything that uses the GPU again, covering all methods in the class. However, this is preferable to the infinite loop that would otherwise be caused.

Examples


# dust includes a couple of very simple examples
filename <- system.file("examples/walk.cpp", package = "dust")

# This model implements a random walk with a parameter coming from
# R representing the standard deviation of the walk
writeLines(readLines(filename))

# The model can be compiled and loaded with dust::dust(filename)
# but it's faster in this example to use the prebuilt version in
# the package
model <- dust::dust_example("walk")

# Print the object and you can see the methods that it provides
model

# Create a model with standard deviation of 1, initial time step zero
# and 30 particles
obj <- model$new(list(sd = 1), 0, 30)
obj

# Curent state is all zero
obj$state()

# Current time is also zero
obj$time()

# Run the model up to time step 100
obj$run(100)

# Reorder/resample the particles:
obj$reorder(sample(30, replace = TRUE))

# See the state again
obj$state()

mrc-ide/dust documentation built on Oct. 14, 2024, 12:57 a.m.