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

`+`

–**Plus**: Both infix (`a + b`

) and prefix (`+a`

) versions supported (e.g.,`1 + 2`

→`3`

)`-`

–**Minus**: Both infix (`a - b`

) and prefix (`-a`

) versions supported (e.g.,`10 - 1`

→`9`

)`*`

–**Multiply**: Multiply two numbers together (e.g.,`2 * 6`

→`12`

)`/`

–**Divide**: Divide two numbers (e.g.,`12 / 6`

→`2`

)`^`

–**Power**: Raise the first number to the power of the second. Either number may be a floating point number (e.g.,`2.3 ^ 1.2`

→`2.716898`

)`(`

–**Parenthesis**: Group expressions together (e.g.,`(1 + 5) * 2`

→`12`

)`if`

–**Conditional**: Inline conditional statement. This takes a form slightly different to typically seen in R with the result of the statement directly assigned (e.g.,`if (9 > 10) 1 else 2`

→`2`

)

Because general programming is not supported in `odin`

and
because every line must contain an assignment, instead of writing

```
if (mycondition) {
a <- true_value
} else {
a <- false_value
}
```

instead write

`a <- if (mycondition) true_value else false_value`

(this works in normal R too!)

There are a group of functions for interacting with arrays

`sum`

–**Sum**: Compute sum over an array, or over some dimension of an array`length`

–**Length**: Total length of an array`dim`

–**Length of one of an array’s dimensions**: If an array`x`

has 10 rows and 20 columns, then`dim(x, 1)`

is`10`

and`dim(x, 2)`

is`20`

. Note that this differs from`dim`

in R and`dim(x, i)`

is closer to`dim(x)[[i]]`

`[`

–**Subset an array**: See below`interpolate`

–**Interpolate an array over time**: See below

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 the derivatives are computed as

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”

(except slightly different, and in C).

Similarly, the expression

`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

Due to constraints with using C, few things can be used as an index; in particular the following will not work:

```
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:

```
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

A number of logical-returning operators exist, primarily to support
the `if`

statement; all the usual comparison operators exist
(though not vectorised `|`

or `&`

).

`>`

–**Greater than**(e.g.,`1 > 2`

→`FALSE`

)`<`

–**Less than**(e.g.,`1 < 2`

→`TRUE`

)`>=`

–**Greater than or equal to**(e.g.,`1 >= 2`

→`FALSE`

)`<=`

–**Less than or equal to**(e.g.,`1 <= 2`

→`TRUE`

)`==`

–**Is exactly equal to**(e.g.,`1 == 1`

→`TRUE`

)`!=`

–**Is not exactly equal to**(e.g.,`1 != 2`

→`TRUE`

)`&&`

–**Boolean AND**(e.g.,`(1 == 1) && (2 > 1)`

→`TRUE`

)`||`

–**Boolean OR**(e.g.,`(1 == 1) && (2 > 1)`

→`TRUE`

)

Be wary of strict equality with `==`

or `!=`

as
numbers may be floating point numbers, which have some surprising
properties for the uninitiated, for example

`sqrt(3)^2 == 3`

`## [1] FALSE`

`%%`

–**Modulo**: Finds the remainder after division of one number by another (e.g.,`123 %% 100`

→`23`

)`%/%`

–**Integer divide**: Different to floating point division, effectively the full number of times one number divides into another (e.g.,`20 %/% 7`

→`2`

)`abs`

–**Absolute value**(e.g.,`abs(-1)`

→`1`

)`sign`

–**Sign function**: Returns the sign of its argument as either -1, 0 or 1, which may be useful for multiplying by another argument (e.g.,`sign(-100)`

→`-1`

)`round`

–**Round a number**(e.g.,`round(1.23)`

→`1`

;`round(1.23, 1)`

→`1.2`

)`floor`

–**Floor of a number**: Largest integer not greater than the provided number (e.g.,`floor(6.5)`

→`6`

)`ceiling`

–**Ceiling of a number**: Smallest integer not less than the provided number (e.g.,`ceiling(6.5)`

→`7`

)`trunc`

–**Truncate a number**: Round a number towards zero`max`

–**Maximum**: Returns maximum of all arguments given (e.g.,`max(2, 6, 1)`

→`6`

)`min`

–**Minimum**(e.g.,`min(2, 6, 1)`

→`1`

)`exp`

–**Exponential function**(e.g.,`exp(1)`

→`2.718282`

)`expm1`

–**Computes exp(x) - 1 accurately for small |x|**(e.g.,`exp(1)`

→`1.718282`

)`log`

–**Logarithmic function**(e.g.,`log(1)`

→`0`

)`log2`

–**Logarithmic function in base 2**(e.g.,`log2(1024)`

→`10`

)`log10`

–**Logarithmic function in base 10**(e.g.,`log10(1000)`

→`3`

)`log1p`

–**Computes log(x + 1) accurately for small |x|**(e.g.,`log1p(1)`

→`0.6931472`

)`sqrt`

–**Square root function**(e.g.,`sqrt(4)`

→`2`

)`beta`

–**Beta function**(e.g.,`beta(3, 5)`

→`0.00952381`

)`lbeta`

–**Log beta function**(e.g.,`lbeta(3, 5)`

→`-4.65396`

)`choose`

–**Binomial coefficients**(e.g.,`choose(60, 3)`

→`34220`

)`lchoose`

–**Log binomial coefficients**(e.g.,`choose(60, 3)`

→`10.44057`

)`gamma`

–**Gamma function**(e.g.,`gamma(10)`

→`362880`

)`lgamma`

–**Log gamma function**(e.g.,`lgamma(10)`

→`12.80183`

)

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:

`cos`

–**Cosine function**`sin`

–**Sine function**`tan`

–**Tangent function**`acos`

–**Arc-cosine function**`asin`

–**Arc-sin function**`atan`

–**Arc-tangent function**`atan2`

–**Two-arg arc-tangent function**`cosh`

–**Hyperbolic cosine function**`sinh`

–**Hyperbolic sine function**`tanh`

–**Hyperbolic tangent function**`acosh`

–**Hyperbolic arc-cosine function**`asinh`

–**Hyperbolic arc-sine function**`atanh`

–**Hyperbolic arc-tangent function**

For discrete time stochastic models, all of R’s normal stochastic distribution functions are available:

`unif_rand`

–**Standard uniform distribution**: Sample from the uniform distribution on [0, 1] - more efficient than but equivalent to runif(0, 1)`norm_rand`

–**Standard normal distribution**: Sample from the standard normal distribution - more efficient than but equivalent to rnorm(0, 1)`exp_rand`

–**Standard exponential distribution**: Sample from the exponential distribution with rate 1 - more efficient than but equivalent to rexp(1)`rbeta`

–**Beta distribution**: With parameters shape1 and shape2 (see`?rbeta`

for details)`rbinom`

–**Binomial distribution**: With parameters`size`

(number of trials) and`prob`

(probability of success)`rcauchy`

–**Cauchy distribution**: With parameters`location`

and`scale`

`rchisq`

–**Chi-Squared distribution**: With parameter`df`

`rexp`

–**Exponential distribution**: With parameter`rate`

`rf`

–**F-distribution**: With parameter`df1`

and `df2`rgamma`

–**Gamma distribution**: With parameters`shape`

and`rate`

`rgeom`

–**Geometric distribution**: Distribution with parameters`prob`

`rhyper`

–**Hypergeometric distribution**: With parameters`m`

(the number of white balls in the urn),`n`

(the number of black balls in the urn) and`k`

(the number of balls drawn from the urn)`rlogis`

–**Logistic distribution**: With parameters`location`

and`scale`

`rlnorm`

–**Log-normal distribution**: With parameters`meanlog`

and`sdlog`

`rnbinom`

–**Negative binomial distribution**: With parameters`size`

,`prob`

and`mu`

`rnorm`

–**Normal distribution**: With parameters`mean`

and`sd`

`rpois`

–**Poisson distribution**: With parameter`lambda`

`rt`

–**Student’s t distribution**: With parameter`df`

`runif`

–**uniform distribution**: With parameters`min`

and`max`

`rweibull`

–**Weibull distribution**: With parameters`shape`

and`scale`

`rwilcox`

–**Wilcoxon rank sum statistic distribution**: With parameters`n`

and`m`

`rsignrank`

–**Wilcoxon signed rank statistic distribution**: With parameter`n`

With random number functions we can write:

`x <- runif(10, 20)`

which will generate a random number from the uniform distribution. If you write:

`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

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:

`rmultinom`

–**multinomial distribution**: The first parameter is the number of samples and the second is the per-class probability and must be a vector`rmhyper`

–**Multivariate hypergeometric distribution**: The first parameter is the number of samples and the second is the per-class count and must be a vector

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

`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

`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 vignette.