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).

dust(
filename,
quiet = FALSE,
workdir = NULL,
gpu = FALSE,
real_type = 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.

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 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 step: the step number

• const double * state: the state at the beginning of the 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).

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. ## See also dust_generator for a description of the class of created objects, and dust_example() for some pre-built examples. If you want to just generate the code and load it yourself with pkgload::load_all or some other means, see dust_generate) ## 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)) #> class walk { #> public: #> using real_type = double; #> using data_type = dust::no_data; #> using internal_type = dust::no_internal; #> using rng_state_type = dust::random::generator<real_type>; #> #> struct shared_type { #> real_type sd; #> }; #> #> walk(const dust::pars_type<walk>& pars) : shared(pars.shared) { #> } #> #> size_t size() const { #> return 1; #> } #> #> std::vector<real_type> initial(size_t step) { #> std::vector<real_type> ret = {0}; #> return ret; #> } #> #> void update(size_t step, const real_type * state, rng_state_type& rng_state, #> real_type * state_next) { #> state_next[0] = state[0] + #> dust::random::normal<real_type>(rng_state, 0, shared->sd); #> } #> #> private: #> dust::shared_ptr<walk> shared; #> }; #> #> namespace dust { #> #> template <> #> dust::pars_type<walk> dust_pars<walk>(cpp11::list pars) { #> walk::real_type sd = cpp11::as_cpp<walk::real_type>(pars["sd"]); #> return dust::pars_type<walk>(walk::shared_type{sd}); #> } #> #> } # 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 #> <dust> object generator #> Public: #> initialize: function (pars, step, n_particles, n_threads = 1L, seed = NULL, #> name: function () #> param: function () #> run: function (step_end) #> simulate: function (step_end) #> set_index: function (index) #> index: function () #> n_threads: function () #> n_state: function () #> n_particles: function () #> n_particles_each: function () #> shape: function () #> update_state: function (pars = NULL, state = NULL, step = NULL, set_initial_state = NULL) #> state: function (index = NULL) #> step: function () #> reorder: function (index) #> resample: function (weights) #> info: function () #> pars: function () #> rng_state: function (first_only = FALSE, last_only = FALSE) #> set_rng_state: function (rng_state) #> has_openmp: function () #> has_gpu_support: function (fake_gpu = FALSE) #> has_compare: function () #> real_size: function () #> rng_algorithm: function () #> uses_gpu: function (fake_gpu = FALSE) #> n_pars: function () #> set_n_threads: function (n_threads) #> set_data: function (data, shared = FALSE) #> compare_data: function () #> filter: function (step_end = NULL, save_trajectories = FALSE, step_snapshot = NULL, #> gpu_info: function () #> Private: #> pars_: NULL #> pars_multi_: NULL #> index_: NULL #> info_: NULL #> n_threads_: NULL #> n_particles_: NULL #> n_particles_each_: NULL #> shape_: NULL #> ptr_: NULL #> gpu_config_: NULL #> methods_: NULL #> param_: NULL #> reload_: NULL #> Parent env: <environment: namespace:dust> #> Locked objects: TRUE #> Locked class: FALSE #> Portable: TRUE # Create a model with standard deviation of 1, initial step zero # and 30 particles obj <- model$new(list(sd = 1), 0, 30)
obj
#> <dust>
#>   Public:
#>     compare_data: function ()
#>     filter: function (step_end = NULL, save_trajectories = FALSE, step_snapshot = NULL,
#>     gpu_info: function ()
#>     has_compare: function ()
#>     has_gpu_support: function (fake_gpu = FALSE)
#>     has_openmp: function ()
#>     index: function ()
#>     info: function ()
#>     initialize: function (pars, step, n_particles, n_threads = 1L, seed = NULL,
#>     n_pars: function ()
#>     n_particles: function ()
#>     n_particles_each: function ()
#>     n_state: function ()
#>     n_threads: function ()
#>     name: function ()
#>     param: function ()
#>     pars: function ()
#>     real_size: function ()
#>     reorder: function (index)
#>     resample: function (weights)
#>     rng_algorithm: function ()
#>     rng_state: function (first_only = FALSE, last_only = FALSE)
#>     run: function (step_end)
#>     set_data: function (data, shared = FALSE)
#>     set_index: function (index)
#>     set_rng_state: function (rng_state)
#>     shape: function ()
#>     simulate: function (step_end)
#>     state: function (index = NULL)
#>     step: function ()
#>     update_state: function (pars = NULL, state = NULL, step = NULL, set_initial_state = NULL)
#>     uses_gpu: function (fake_gpu = FALSE)
#>   Private:
#>     gpu_config_: NULL
#>     index_: NULL
#>     info_: NULL
#>     methods_: list
#>     n_particles_: 30
#>     n_particles_each_: 30
#>     param_: NULL
#>     pars_: list
#>     pars_multi_: FALSE
#>     ptr_: externalptr
#>     shape_: 30

# Curent state is all zero
obj$state() #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14] #> [1,] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 #> [,15] [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] [,24] [,25] [,26] #> [1,] 0 0 0 0 0 0 0 0 0 0 0 0 #> [,27] [,28] [,29] [,30] #> [1,] 0 0 0 0 # Current step is also zero obj$step()
#> [1] 0

# Run the model up to step 100
obj$run(100) #> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] #> [1,] 8.575403 -3.856778 -5.096966 14.6775 -15.44855 6.08033 4.434121 -4.934855 #> [,9] [,10] [,11] [,12] [,13] [,14] [,15] #> [1,] 6.561592 21.90629 -6.113029 -0.4245573 3.757224 -6.83798 -4.262336 #> [,16] [,17] [,18] [,19] [,20] [,21] [,22] [,23] #> [1,] 5.786386 -5.858829 -22.99371 -5.320212 6.420181 -30.38696 3.76536 20.639 #> [,24] [,25] [,26] [,27] [,28] [,29] [,30] #> [1,] -10.21101 -2.108581 -7.354048 -10.74424 5.05436 -0.9596275 -7.513779 # Reorder/resample the particles: obj$reorder(sample(30, replace = TRUE))

# See the state again
obj\$state()
#>         [,1]     [,2]      [,3]      [,4]      [,5]      [,6]    [,7]      [,8]
#> [1,] 14.6775 4.434121 -4.934855 -5.096966 -4.262336 -4.934855 14.6775 -10.74424
#>         [,9]     [,10]     [,11]   [,12]     [,13]      [,14]   [,15]   [,16]
#> [1,] 14.6775 -7.354048 -10.74424 3.76536 -5.858829 -0.9596275 5.05436 3.76536
#>          [,17]      [,18]    [,19]     [,20]     [,21]     [,22]     [,23]
#> [1,] -2.108581 -0.4245573 4.434121 -6.113029 -5.858829 -7.354048 -5.320212
#>          [,24]    [,25]     [,26]     [,27]   [,28]     [,29]    [,30]
#> [1,] -22.99371 3.757224 -3.856778 -4.262336 5.05436 -10.74424 6.561592