Merge pull request #98 from mrc-ide/mrc-5865

Import and adapt docs from dust1

DESCRIPTION

@@ -27,6 +27,7 @@ Suggests:
     callr,
     cpp11,
     decor,
+    fs,
     glue,
     knitr,
     mockery,

R/browser.R

@@ -1,5 +1,5 @@
 ##' Control browser-based debugging of dust models.  This help page documents
-##' three funtions that can be used to control if and how the browser
+##' 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).
@@ -9,7 +9,7 @@
 ##' 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 presssing `<enter>`.
+##' 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

README.md

@@ -1,11 +1,30 @@
-# dust2
+# dust2 <img src='man/figures/logo.png' align="right" height="139" />
 
 <!-- badges: start -->
 [![Project Status: Concept – Minimal or no implementation has been done yet, or the repository is only intended to be a limited example, demo, or proof-of-concept.](https://www.repostatus.org/badges/latest/concept.svg)](https://www.repostatus.org/#concept)
 [![R-CMD-check](https://github.com/mrc-ide/dust2/actions/workflows/R-CMD-check.yaml/badge.svg?branch=main)](https://github.com/mrc-ide/dust2/actions/workflows/R-CMD-check.yaml)
 [![codecov.io](https://codecov.io/github/mrc-ide/dust2/coverage.svg?branch=main)](https://codecov.io/github/mrc-ide/dust2?branch=main)
 <!-- badges: end -->
 
+The `dust2` package provides an engine for running dynamical systems in discrete or continuous time and where the processes are stochastic or deterministic.  We focus on Markov models where the problem reduces to describing how model state changes as a function of its current state (and possibly time) but without reference to where it has come from.  Superficially, the problem is not very hard (see `vignette("design")`), but `dust2` takes care of many practical and bookkeeping details such as:
+
+* Running systems in parallel on multi-core machines, even those involving random numbers
+* Providing useful verbs for efficiently working with groups of simulations (different parameters, starting conditions or stochastic realisations)
+* Comparing simulations to time-series of data, and implementing sequential Monte Carlo methods such as a bootstrap particle filter
+
+## Get started
+
+* `vignette("design")` describes the problems `dust` tries to solve ([read on package website](https://mrc-ide.github.io/dust/articles/design.html))
+* `vignette("dust2")` describes `dust` by example, showing two simple systems and the methods that can drive them ([read on package website](https://mrc-ide.github.io/dust/articles/dust.html))
+* [`odin2`](https://mrc-ide.github.io/odin2) is the way most `dust2` systems are written
+* The [odin & monty book](https://mrc-ide.github.io/odin-monty) shows how the package works in context
+* `vignette("writing")` shows how to write a `dust2` system by hand in C++
+* `vignette("packaging")` describes how to package a system for easy reuse
+
+## Roadmap
+
+This package is a ground-up rewrite of [`dust`](https://mrc-ide.github.io/dust) and will eventually become version 2.0.0 of `dust`, which we will then release to CRAN.  It exists separately for now to facilitate development and use alongside the original `dust`, and is being developed in parallel with [`odin2`](https://mrc-ide.github.io/odin2) and [`monty`](https://mrc-ide.github.io/monty) (previously [`mcstate`](https://mrc-ide.github.io/mcstate)).  Some of the functionality here was originally found in `mcstate` and some of the previous version of `dust` can now be found in `monty` (e.g., the random number library).
+
 ## Installation
 
 To install `dust2`:

_pkgdown.yml

@@ -64,8 +64,16 @@ reference:
       - dust_package
 
 articles:
+  - title: Introduction
+    navbar: ~
+    contents:
+      - dust2
+      - design
   - title: Details
     navbar: Details
     contents:
+      - data
       - details
+      - packaging
       - periodic
+      - writing

inst/WORDLIST

@@ -1,21 +1,35 @@
+Blackmann
 CMD
+CMRG
+CRAN's
 Dormand
 HPC
+Kalman
+L'Ecuyer
 Makevars
+Mersenne
 ODEs
 OpenMP
+Poisson
 R's
+SEIR
+Vigna
 codecov
 cpp
 etc
 io
 lfloor
 monty
+odin
 parallelisable
 pkgload
 resettings
 rfloor
 rng
 roxygen
+stochasticity
 th
+truthy
 unfilter
+unitless
+xoshiro

vignettes/data.Rmd

@@ -0,0 +1,145 @@
+---
+title: "Comparing dust systems to data"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Comparing dust systems to data}
+  %\VignetteEngine{knitr::rmarkdown}
+  %\VignetteEncoding{UTF-8}
+---
+
+```{r, include = FALSE}
+source("support.R")
+set.seed(1)
+```
+One of our aims with `dust` was to enable the creation of fast particle filters.  This vignette outlines the steps in implementing the comparison directly as part of the model, and compares the result against a known deterministic result.
+
+```{r}
+library(dust2)
+```
+
+We start with a simple example, a model of volatility.  For completeness (and because the maths around this will need adding later!) we show the full code here:
+
+```{r echo = FALSE, results = "asis"}
+volatility_code <- readLines("examples/volatility.cpp")
+cc_output(volatility_code)
+```
+
+```{r}
+volatility <- dust_compile("examples/volatility.cpp", quiet = TRUE)
+```
+
+To demonstrate the approach, we simulate some data from the model itself.  We have to include the simulation of the observation process here which adds an additional sample from the normal distribution.
+
+These are the parameters we will use:
+
+```{r}
+pars <- list(
+  # Generation process
+  alpha = 0.91,
+  sigma = 1,
+  # Observation process
+  gamma = 1,
+  tau = 1)
+```
+
+```{r}
+data <- local({
+  sys <- dust_system_create(volatility(), pars)
+  dust_system_set_state_initial(sys)
+  times <- seq(1, 100, by = 1)
+  observed <- drop(dust_system_simulate(sys, times)) + rnorm(length(times))
+  data.frame(time = times, observed = observed)
+})
+head(data)
+plot(observed ~ time, data, type = "o", pch = 19, las = 1)
+```
+
+Now, we construct a particle filter:
+
+
+```{r}
+filter <- dust_filter_create(volatility(), 0, data, n_particles = 1000)
+filter
+```
+
+Running the particle filter simulates the process on all $10^3$ particles and compares at each timestep the simulated data with your observed data using the provided comparison function. It returns the log-likelihood:
+
+```{r}
+dust_likelihood_run(filter, pars)
+```
+
+This is stochastic and each time you run it, the estimate will differ:
+
+```{r}
+dust_likelihood_run(filter, pars)
+```
+
+In this case the model is simple enough that we can use a [Kalman Filter](https://en.wikipedia.org/wiki/Kalman_filter) to calculate the likelihood exactly:
+
+```{r volatility_kalman}
+kalman_filter <- function(pars, data) {
+  y <- data$observed
+
+  mu <- 0
+  s <- 1
+  log_likelihood <- 0
+
+  for (t in seq_along(y)) {
+    mu <- pars$alpha * mu
+    s <- pars$alpha^2 * s + pars$sigma^2
+    m <- pars$gamma * mu
+
+    S <- pars$gamma^2 * s + pars$tau^2
+    K <- pars$gamma * s / S
+
+    mu <- mu + K * (y[t] - m)
+    s <- s - pars$gamma * K * s
+
+    log_likelihood <- log_likelihood + dnorm(y[t], m, sqrt(S), log = TRUE)
+  }
+
+  log_likelihood
+}
+ll_k <- kalman_filter(pars, data)
+ll_k
+```
+
+Unlike the particle filter the Kalman filter is deterministic:
+
+```{r}
+kalman_filter(pars, data)
+```
+
+```{r}
+ll <- replicate(200, dust_likelihood_run(filter, pars))
+hist(ll, col = "steelblue3")
+abline(v = ll_k, col = "red", lty = 2, lwd = 2)
+```
+
+As the number of particles used changes, the variance of this estimate will change.  For example, here is a filter with only 100 particles (1/10th the number as in the previous).  We plot the range of observed likelihoods from the larger filter as orange vertical lines.
+
+```{r}
+filter2 <- dust_filter_create(volatility(), 0, data, n_particles = 100)
+ll2 <- replicate(200, dust_likelihood_run(filter2, pars))
+hist(ll2, col = "steelblue3")
+abline(v = ll_k, col = "red", lty = 2, lwd = 2)
+abline(v = range(ll), col = "orange", lty = 3, lwd = 2)
+```
+
+If you run a particle filter with `save_history = TRUE`, it will record the (filtered) trajectories, which you can then extract with `dust_filter_last_history()`:
+
+```{r}
+dust_likelihood_run(filter, pars, save_history = TRUE)
+history <- dust_likelihood_last_history(filter)
+dim(history)
+```
+
+This is a `n_state` (here 1) x `n_particles` (1000) x `n_time` (100) 3d array, but we will drop the first rank of this for plotting, and transpose so that time is the first axis:
+
+```{r}
+matplot(data$time, t(drop(history)), xlab = "Time", ylab = "Value",
+        las = 1, type = "l", lty = 1, col = "#00000002")
+points(observed ~ time, data, col = "red", pch = 19)
+```
+
+Here you can see the sampled trajectories fitting the observed data, with fewer extant trajectories the further back in time you go (fewer, thicker, lines due to the sampling process).