--- title: "odin functions" author: "Rich FitzJohn" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true vignette: > %\VignetteIndexEntry{odin functions} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ``` {r include = FALSE} path <- system.file("functions.json", package = "odin", mustWork = TRUE) functions <- jsonlite::fromJSON(path, FALSE) vcapply <- odin:::vcapply render <- function(x) { render1 <- function(el) { if (!is.null(el$description)) { ret <- sprintf("**%s**: %s", el$title, el$description) } else { ret <- sprintf("**%s**", el$title) } if (!is.null(el$examples)) { examples <- vcapply(el$examples, function(ex) sprintf("`%s` → `%s`", ex[[1]], ex[[2]])) ret <- sprintf("%s (e.g., %s)", ret, paste(examples, collapse = "; ")) } ret } value <- vapply(x, render1, "", USE.NAMES = FALSE) paste(sprintf("* `%s` -- %s\n", names(x), value), collapse = "\n") } ``` `odin` supports many functions that you'd expect to see for constructing differential equation models; primarily mathematical functions available through R's "Rmath" library. These include all mathematical operations, and many more obscure mathematical functions. Special support is provided for working with arrays. A further set of functions is available for working discrete time stochastic models. ## Basic operators ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$basic)) ``` Because general programming is not supported in `odin` and because every line must contain an assignment, instead of writing ```r if (mycondition) { a <- true_value } else { a <- false_value } ``` instead write ```r a <- if (mycondition) true_value else false_value ``` (this works in normal R too!) ## Array support There are a group of functions for interacting with arrays ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$arrays)) ``` When working with arrays, use generally implies a "for loop" in the generated C code. For example, in the example in [the main package vignette](odin.html) the derivatives are computed as ```r deriv(y[]) <- r[i] * y[i] * (1 - sum(ay[i, ])) ``` The indexes on the right hand side can be one of `i`, `j`, `k`, `l` `i5`, `i6`, `i7` or `i8` corresponding to the index on the *left hand side* being iterated over (`odin` supports arrays up to 8 dimensions). The left-hand-side here contains no explicit entry (`y[]`) which is equivalent to `y[1:length(y)]`, which expands (approximately) to the "for loop" ```r for (i in 1:length(y)) { deriv(y[i]) <- r[i] * y[i] * (1 - sum(ay[i, ])) } ``` (except slightly different, and in C). Similarly, the expression ```r ay[, ] <- a[i, j] * y[j] ``` involves loops over two dimensions (`ay[, ]` becomes `ay[1:dim(ay, 1), 1:dim(ay, 2)]` and so the loop becomes ```r for (i in 1:dim(ay, 1)) { for (j in 1:dim(ay, 2)) { ay[i, j] <- a[i, j] * y[j] } } ``` Due to constraints with using C, few things can be used as an index; in particular the following will not work: ```r idx <- if (t > 5) 2 else 1 x <- vec[idx] ``` (or where `idx` is some general odin variable as the result of a different assignment). You must use `as.integer` to cast this to integer immediately before indexing: ```r idx <- if (t > 5) 2 else 1 x <- vec[as.integer(idx)] ``` This will *truncate* the value (same behaviour as `truncate`) so be warned if passing in things that may be approximately integer - you may want to use `as.integer(round(x))` in that case. The interpolation functions are described in more detail in the [main package vignette](odin.html) ## Operators A number of logical-returning operators exist, primarily to support the `if` statement; all the usual comparison operators exist (though not vectorised `|` or `&`). ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$operators)) ``` Be wary of strict equality with `==` or `!=` as numbers may be floating point numbers, which have some surprising properties for the uninitiated, for example ``` {r } sqrt(3)^2 == 3 ``` ## Mathematical operators ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$maths)) ``` The exact for `%%` and `%/%` for floating point numbers and signed numbers are complicated - please see `?Arithmetic`. The rules for operators in `odin` are exactly those in R as the same underlying functions are used. Similarly, for the differences between `round`, `floor`, `ceiling` and `truncate`, see the help page `?round`. Note that R's behaviour for rounding away from 0.5 is exactly followed and that this slightly changed behaviour at version 4.0.0 All the usual trig functions are also available: ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$trig)) ``` ## Stochastic models For discrete time stochastic models, all of R's normal stochastic distribution functions are available: ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$stochastic)) ``` With random number functions we can write: ```r x <- runif(10, 20) ``` which will generate a random number from the uniform distribution. If you write: ```r x[] <- runif(10, 20) ``` then each element of `x` will be filled with a *different* random number drawn from this distribution (which is generally what you want). Random numbers are considered to be *time varying* which means they will automatically generate each time step, so if you write ```r x <- rnorm(0, 10) update(y[]) <- y[i] + x ``` then at each time step, each element of `y` will be updated by the same random number from a normal distribution with a mean of zero and a standard deviation of 10 - the number will change each time step but be the same for each element of `y` in the example above. In addition, two functions that are vector returning and require some care to use: ``` {r echo = FALSE, results = "asis"} writeLines(render(functions$inplace)) ``` Both these functions require a vector input (of probabilities for `rmultinom` and of counts for `rmhyper`) and return a vector the same length. So the expression ```r y[] <- rmultinom(10, p) ``` will produce a vector `y` of samples from the multinomial distribution with parameters `size = 10` (so after wards `sum(y)` is 10) and probabilities `p`. It is very important that `y` and `p` have the same size. At the moment it is not possible to use expressions like ```r y[1, ] <- rmultinom(10, p[i, ]) ``` but this is planned for implementation in the future. A full example of using `rmultinom` is given in the [discrete models](discrete.html) vignette.