From fb80566a212fecbd36a7ba3877bcac5e83b8532b Mon Sep 17 00:00:00 2001 From: edknock Date: Thu, 7 Nov 2024 11:55:42 +0000 Subject: [PATCH 1/5] add vaccination example for arrays --- arrays.qmd | 209 ++++++++++++++++++++++++++++++++++------------------- 1 file changed, 133 insertions(+), 76 deletions(-) diff --git a/arrays.qmd b/arrays.qmd index ff35f05..0bdb89c 100644 --- a/arrays.qmd +++ b/arrays.qmd @@ -5,30 +5,67 @@ source("common.R") ``` -::: {.callout-warning} -These examples come from the 2023 teaching slides, and have not been turned into actual content yet. -::: +## When scalars are just not enough -When scalars are just not enough +The aim of this section is to show you how you use odin's array syntax, via some motivating examples. The examples here are necessarily a bit longer than in the previous sections because we will generally need a few more moving parts. -The aim of this section is to show you how you use odin's array syntax, via some motivating examples. We may need to break this up into several sections! The examples here are necessarily a bit longer than in the previous sections because we will generally need a few more moving parts. +## Example: an age structured SIR model {#sec-stochastic-age} -* Need a good motivating example for an array - previously we used an age structured model -* show example code +Let's take the simple stochastic SIR model from @sec-stochastic-sir and add age structure to it, with heterogeneous mixing between age groups -``` -update(S) <- S - n_SI -update(I) <- I + n_SI - n_IR -update(R) <- R + n_IR -``` - -``` -update(S[]) <- S[i] - n_SI[i] -update(I[]) <- I[i] + n_SI[i] - n_IR[i] -update(R[]) <- R[i] + n_IR[i] +```r +odin2::odin({ + # Equations for transitions between compartments by age group + update(S[]) <- S[i] - n_SI[i] + update(I[]) <- I[i] + n_SI[i] - n_IR[i] + update(R[]) <- R[i] + n_IR[i] + + # Individual probabilities of transition: + p_SI[] <- 1 - exp(-lambda[i] * dt) # S to I + p_IR <- 1 - exp(-gamma * dt) # I to R + + # Calculate force of infection + ## age-structured contact matrix: m[i, j] is mean number of contacts an + ## individual in group i has with an individual in group j per time unit + m <- parameter() + ## here s_ij[i, j] gives the mean number of contacts and individual in group + ## i will have with the currently infectious individuals of group j + s_ij[, ] <- m[i, j] * I[j] + ## lambda[i] is the total force of infection on an individual in group i + lambda[] <- beta * sum(s_ij[i, ]) + + # Draws from binomial distributions for numbers changing between + # compartments: + n_SI[] <- Binomial(S[i], p_SI[i]) + n_IR[] <- Binomial(I[i], p_IR) + + initial(S[]) <- S0[i] + initial(I[]) <- I0[i] + initial(R[]) <- 0 + + # User defined parameters - default in parentheses: + S0 <- parameter() + I0 <- parameter() + beta <- parameter(0.0165) + gamma <- parameter(0.1) + + # Dimensions of arrays + N_age <- parameter() + dim(S0) <- N_age + dim(I0) <- N_age + dim(S) <- N_age + dim(I) <- N_age + dim(R) <- N_age + dim(n_SI) <- N_age + dim(n_IR) <- N_age + dim(p_SI) <- N_age + dim(m) <- c(N_age, N_age) + dim(s_ij) <- c(N_age, N_age) + dim(lambda) <- N_age +}) ``` -The odin code +In the odin code above ``` update(S[]) <- S[i] - n_SI[i] @@ -42,68 +79,88 @@ for (int i = 0; i < S_length; ++i) { } ``` -* Don't use index variables (`i`, `j`, `k`, etc) on the left hand side -* Use multiple lines to account for boundary conditions -* Can crash the program if out of bounds (we have plans for this) - -## An age structured SIR model - -```r -# Equations for transitions between compartments by age group -update(S[]) <- S[i] - n_SI[i] -update(I[]) <- I[i] + n_SI[i] - n_IR[i] -update(R[]) <- R[i] + n_IR[i] - -# Individual probabilities of transition: -p_SI[] <- 1 - exp(-lambda[i] * dt) # S to I -p_IR <- 1 - exp(-gamma * dt) # I to R +so there is an implicit indexing by `i` on the LHS in that equation of the odin code, and the generated code will then loop over all values of `i`. -# Force of infection -m <- parameter() # age-structured contact matrix -s_ij[, ] <- m[i, j] * I[j] -lambda[] <- beta * sum(s_ij[i, ]) -# Draws from binomial distributions for numbers changing between -# compartments: -n_SI[] <- Binomial(S[i], p_SI[i]) -n_IR[] <- Binomial(I[i], p_IR) +## Another example: an age-structured SIR model with vaccination +Let's take the above model and additionally add some vaccination to it. -initial(S[]) <- S_ini[i] -initial(I[]) <- I_ini[i] -initial(R[]) <- 0 - -# User defined parameters - default in parentheses: -S_ini <- parameter() -I_ini <- parameter() -beta <- parameter(0.0165) -gamma <- parameter(0.1) - -# Dimensions of arrays -N_age <- parameter() -dim(S_ini) <- N_age -dim(I_ini) <- N_age -dim(S) <- N_age -dim(I) <- N_age -dim(R) <- N_age -dim(n_SI) <- N_age -dim(n_IR) <- N_age -dim(p_SI) <- N_age -dim(m) <- c(N_age, N_age) -dim(s_ij) <- c(N_age, N_age) -dim(lambda) <- N_age +``` +odin2::odin({ + # Equations for transitions between compartments by age group + update(S[, ]) <- new_S[i, j] + update(I[, ]) <- I[i, j] + n_SI[i, j] - n_IR[i, j] + update(R[, ]) <- R[i, j] + n_IR[i, j] + + # Individual probabilities of transition: + p_SI[, ] <- 1 - exp(-rel_susceptibility[i, j] * lambda[i] * dt) # S to I + p_IR <- 1 - exp(-gamma * dt) # I to R + p_vax[, ] <- 1 - exp(-eta[i, j] * dt) + + # Force of infection + m <- parameter() # age-structured contact matrix + s_ij[, ] <- m[i, j] * sum(I[j, ]) + lambda[] <- beta * sum(s_ij[i, ]) + + # Draws from binomial distributions for numbers changing between + # compartments: + n_SI[, ] <- Binomial(S[i, j], p_SI[i, j]) + n_IR[, ] <- Binomial(I[i, j], p_IR) + + # Nested binomial draw for vaccination in S + # Assume you cannot move vaccine class and get infected in same step + n_S_vax[, ] <- Binomial(S[i, j] - n_SI[i, j], p_vax[i, j]) + new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, N_vax] + new_S[, 2:N_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1]) + + initial(S[, ]) <- S0[i, j] + initial(I[, ]) <- I0[i, j] + initial(R[, ]) <- 0 + + # User defined parameters - default in parentheses: + S0 <- parameter() + I0 <- parameter() + beta <- parameter(0.0165) + gamma <- parameter(0.1) + eta <- parameter() + rel_susceptibility[, ] <- parameter() + + # Dimensions of arrays + N_age <- parameter() + N_vax <- parameter() + dim(S0) <- c(N_age, N_vax) + dim(I0) <- c(N_age, N_vax) + dim(S) <- c(N_age, N_vax) + dim(I) <- c(N_age, N_vax) + dim(R) <- c(N_age, N_vax) + dim(n_SI) <- c(N_age, N_vax) + dim(n_IR) <- c(N_age, N_vax) + dim(p_SI) <- c(N_age, N_vax) + dim(m) <- c(N_age, N_age) + dim(s_ij) <- c(N_age, N_age) + dim(lambda) <- N_age + dim(eta) <- c(N_age, N_vax) + dim(rel_susceptibility) <- c(N_age, N_vax) + dim(p_vax) <- c(N_age, N_vax) + dim(n_S_vax) <- c(N_age, N_vax) + dim(new_S) <- c(N_age, N_vax) +}) ``` -Relevant changes from the above: - +We see we can use multiple lines to deal with boundary conditions ``` -m <- parameter() # age-structured contact matrix -s_ij[, ] <- m[i, j] * I[j] -lambda[] <- beta * sum(s_ij[i, ]) -p_SI[] <- 1 - exp(-lambda[i] * dt) -update(S[]) <- S[i] - n_SI[i] -N_age <- user() -dim(S) <- N_age -dim(m) <- c(N_age, N_age) +new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, N_vax] +new_S[, 2:N_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1] ``` - -## Vaccination +which we could also write as +``` +new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] +new_S[, 1] <- new_S[i, j] + n_S_vax[i, N_vax] +new_S[, 2:N_vax] <- new_S[i, j] + n_S_vax[i, j - 1] +``` +or another alternative way of writing this would be to use `if else` +``` +new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + + (if (j == 1) n_S_vax[i, N_vax] else n_S_vax[i, j - 1]) +``` +Note that in odin, an `if` always requires an `else`! \ No newline at end of file From 57ef95ee706a7cf38f9fad8180745cae8cd30f3b Mon Sep 17 00:00:00 2001 From: edknock Date: Thu, 7 Nov 2024 15:08:10 +0000 Subject: [PATCH 2/5] more explanation of array syntax --- arrays.qmd | 73 +++++++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 70 insertions(+), 3 deletions(-) diff --git a/arrays.qmd b/arrays.qmd index 0bdb89c..dd26819 100644 --- a/arrays.qmd +++ b/arrays.qmd @@ -93,7 +93,7 @@ odin2::odin({ update(R[, ]) <- R[i, j] + n_IR[i, j] # Individual probabilities of transition: - p_SI[, ] <- 1 - exp(-rel_susceptibility[i, j] * lambda[i] * dt) # S to I + p_SI[, ] <- 1 - exp(-rel_susceptibility[j] * lambda[i] * dt) # S to I p_IR <- 1 - exp(-gamma * dt) # I to R p_vax[, ] <- 1 - exp(-eta[i, j] * dt) @@ -140,7 +140,7 @@ odin2::odin({ dim(s_ij) <- c(N_age, N_age) dim(lambda) <- N_age dim(eta) <- c(N_age, N_vax) - dim(rel_susceptibility) <- c(N_age, N_vax) + dim(rel_susceptibility) <- c(N_vax) dim(p_vax) <- c(N_age, N_vax) dim(n_S_vax) <- c(N_age, N_vax) dim(new_S) <- c(N_age, N_vax) @@ -163,4 +163,71 @@ or another alternative way of writing this would be to use `if else` new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + (if (j == 1) n_S_vax[i, N_vax] else n_S_vax[i, j - 1]) ``` -Note that in odin, an `if` always requires an `else`! \ No newline at end of file +Note that in odin, an `if` always requires an `else`! + +## Indexing + +We have seen in the above examples uses of 1d and 2d arrays, making use of index variables `i` and `j`. Note that while we never use these index variables on the LHS, it is implicit that on the LHS we are indexing by `i` for 1d arrays and `i` and `j` for 2d arrays! Use of these index variables on the RHS corresponds to the indexing on the LHS. + +Of course you may want to use higher-dimensional arrays! We currently supporting up to 8 dimensions, with index variables `i`, `j`, `k`, `l`, `i5`, `i6`, `i7` and `i8`. + +Be careful with array equations! It's possible that in an array equation you end up going out of bounds with an array used on the RHS. This can crash the program when running. + +## Declaring dimensions +You must declare the dimensions of all arrays within your `odin` code using a `dim` equation. + +Dimensions can be hardcoded: +``` +dim(y) <- 10 +``` +or dimensions themselves can be parameters: +``` +dim(x) <- c(nr, 5) +nr <- parameter() +``` +You may have array parameters, these can have dimensions explicitly declared +``` +dim(m) <- c(N_age, N_age) +``` +or they can also have their dimensions detected when the parameters are passed into the system - this still needs a `dim` equation where you must explicitly state the rank: +``` +dim(m) <- parameter(rank = 2) +``` + +## Semantics of random number draws + +Stochastic functions are called for each element in an array they are assigned to, at each time. So here: + +``` +x[] <- Normal(0, 1) +``` + +`x` will be filled with each element having a different draw from a standard normal. In contrast, in: + +``` +a <- Normal(0, 1) +x[] <- a +``` + +`x` will be a vector where every element is the same, the result of a *single* draw from a standard normal. + + +## Summing over arrays +Frequently, you will want to take a sum over an array, or part of an array, using `sum`. To sum over all elements of an array, use `sum()` with the name of the array you would like to sum over: + +```r +dim(x) <- 10 +x_tot <- sum(x) +``` + +If `m` is a matrix you can compute the sums over the second column by writing: + +```r +m1_tot <- sum(m[, 2]) +``` + +This partial sum approach is frequently used within implicit loops: + +``` +m_col_totals[] <- sum(m[, i]) +``` \ No newline at end of file From ed156d287af153566735e9ed78686cb0aa246f3f Mon Sep 17 00:00:00 2001 From: edknock Date: Thu, 7 Nov 2024 15:35:58 +0000 Subject: [PATCH 3/5] some formatting --- arrays.qmd | 41 +++++++++++++++++++++++------------------ 1 file changed, 23 insertions(+), 18 deletions(-) diff --git a/arrays.qmd b/arrays.qmd index dd26819..23da684 100644 --- a/arrays.qmd +++ b/arrays.qmd @@ -9,12 +9,17 @@ source("common.R") The aim of this section is to show you how you use odin's array syntax, via some motivating examples. The examples here are necessarily a bit longer than in the previous sections because we will generally need a few more moving parts. +```{r} +library(odin2) +library(dust2) +``` + ## Example: an age structured SIR model {#sec-stochastic-age} Let's take the simple stochastic SIR model from @sec-stochastic-sir and add age structure to it, with heterogeneous mixing between age groups -```r -odin2::odin({ +```{r} +sir_age <- odin({ # Equations for transitions between compartments by age group update(S[]) <- S[i] - n_SI[i] update(I[]) <- I[i] + n_SI[i] - n_IR[i] @@ -67,13 +72,13 @@ odin2::odin({ In the odin code above -``` +```r update(S[]) <- S[i] - n_SI[i] ``` becomes (approximately) -``` +```r for (int i = 0; i < S_length; ++i) { update_S[i] = S[i] + n_SI[i]; } @@ -85,8 +90,8 @@ so there is an implicit indexing by `i` on the LHS in that equation of the odin ## Another example: an age-structured SIR model with vaccination Let's take the above model and additionally add some vaccination to it. -``` -odin2::odin({ +```{r} +sir_age_vax <- odin({ # Equations for transitions between compartments by age group update(S[, ]) <- new_S[i, j] update(I[, ]) <- I[i, j] + n_SI[i, j] - n_IR[i, j] @@ -111,7 +116,7 @@ odin2::odin({ # Assume you cannot move vaccine class and get infected in same step n_S_vax[, ] <- Binomial(S[i, j] - n_SI[i, j], p_vax[i, j]) new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, N_vax] - new_S[, 2:N_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1]) + new_S[, 2:N_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1] initial(S[, ]) <- S0[i, j] initial(I[, ]) <- I0[i, j] @@ -147,19 +152,19 @@ odin2::odin({ }) ``` -We see we can use multiple lines to deal with boundary conditions -``` +We see we can use multiple lines to deal with boundary conditions: +```r new_S[, 1] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, N_vax] new_S[, 2:N_vax] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + n_S_vax[i, j - 1] ``` which we could also write as -``` +```r new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] new_S[, 1] <- new_S[i, j] + n_S_vax[i, N_vax] new_S[, 2:N_vax] <- new_S[i, j] + n_S_vax[i, j - 1] ``` or another alternative way of writing this would be to use `if else` -``` +```r new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + (if (j == 1) n_S_vax[i, N_vax] else n_S_vax[i, j - 1]) ``` @@ -177,20 +182,20 @@ Be careful with array equations! It's possible that in an array equation you end You must declare the dimensions of all arrays within your `odin` code using a `dim` equation. Dimensions can be hardcoded: -``` +```r dim(y) <- 10 ``` or dimensions themselves can be parameters: -``` +```r dim(x) <- c(nr, 5) nr <- parameter() ``` You may have array parameters, these can have dimensions explicitly declared -``` +```r dim(m) <- c(N_age, N_age) ``` or they can also have their dimensions detected when the parameters are passed into the system - this still needs a `dim` equation where you must explicitly state the rank: -``` +```r dim(m) <- parameter(rank = 2) ``` @@ -198,13 +203,13 @@ dim(m) <- parameter(rank = 2) Stochastic functions are called for each element in an array they are assigned to, at each time. So here: -``` +```r x[] <- Normal(0, 1) ``` `x` will be filled with each element having a different draw from a standard normal. In contrast, in: -``` +```r a <- Normal(0, 1) x[] <- a ``` @@ -228,6 +233,6 @@ m1_tot <- sum(m[, 2]) This partial sum approach is frequently used within implicit loops: -``` +```r m_col_totals[] <- sum(m[, i]) ``` \ No newline at end of file From f4a3412b63821aeef1a89c7df4fd3164bb8b3848 Mon Sep 17 00:00:00 2001 From: edknock <47318334+edknock@users.noreply.github.com> Date: Thu, 7 Nov 2024 16:39:56 +0000 Subject: [PATCH 4/5] Apply Rich's suggestions Co-authored-by: Rich FitzJohn --- arrays.qmd | 3 +++ 1 file changed, 3 insertions(+) diff --git a/arrays.qmd b/arrays.qmd index 23da684..e0ac0bb 100644 --- a/arrays.qmd +++ b/arrays.qmd @@ -88,6 +88,7 @@ so there is an implicit indexing by `i` on the LHS in that equation of the odin ## Another example: an age-structured SIR model with vaccination + Let's take the above model and additionally add some vaccination to it. ```{r} @@ -179,6 +180,7 @@ Of course you may want to use higher-dimensional arrays! We currently supporting Be careful with array equations! It's possible that in an array equation you end up going out of bounds with an array used on the RHS. This can crash the program when running. ## Declaring dimensions + You must declare the dimensions of all arrays within your `odin` code using a `dim` equation. Dimensions can be hardcoded: @@ -217,6 +219,7 @@ x[] <- a `x` will be a vector where every element is the same, the result of a *single* draw from a standard normal. + ## Summing over arrays Frequently, you will want to take a sum over an array, or part of an array, using `sum`. To sum over all elements of an array, use `sum()` with the name of the array you would like to sum over: From 8f8e0410fe6e52443cc50dde54e00d3338c5adc2 Mon Sep 17 00:00:00 2001 From: edknock Date: Fri, 8 Nov 2024 12:07:36 +0000 Subject: [PATCH 5/5] add 2 group example --- arrays.qmd | 212 +++++++++++++++++++++++++++++++++++++---------------- 1 file changed, 148 insertions(+), 64 deletions(-) diff --git a/arrays.qmd b/arrays.qmd index e0ac0bb..ca464bc 100644 --- a/arrays.qmd +++ b/arrays.qmd @@ -14,9 +14,139 @@ library(odin2) library(dust2) ``` -## Example: an age structured SIR model {#sec-stochastic-age} -Let's take the simple stochastic SIR model from @sec-stochastic-sir and add age structure to it, with heterogeneous mixing between age groups +## Example: a two-group age structured SIR model + +Let's take the simple stochastic SIR model from @sec-stochastic-sir and add age structure to it by having two groups (children and adults), with heterogeneous mixing between the groups. + +```{r} +sir_2 <- odin({ + # Equations for transitions between compartments by age group + update(S[]) <- S[i] - n_SI[i] + update(I[]) <- I[i] + n_SI[i] - n_IR[i] + update(R[]) <- R[i] + n_IR[i] + + # Individual probabilities of transition: + p_SI[] <- 1 - exp(-lambda[i] * dt) # S to I + p_IR <- 1 - exp(-gamma * dt) # I to R + + # Calculate force of infection + ## age-structured contact matrix: m[i, j] is mean number of contacts an + ## individual in group i has with an individual in group j per time unit + m <- parameter() + ## here s_ij[i, j] gives the mean number of contacts and individual in group + ## i will have with the currently infectious individuals of group j + s_ij[, ] <- m[i, j] * I[j] + ## lambda[i] is the total force of infection on an individual in group i + lambda[] <- beta * (s_ij[i, 1] + s_ij[i, 2]) + + # Draws from binomial distributions for numbers changing between + # compartments: + n_SI[] <- Binomial(S[i], p_SI[i]) + n_IR[] <- Binomial(I[i], p_IR) + + initial(S[]) <- S0[i] + initial(I[]) <- I0[i] + initial(R[]) <- 0 + + # User defined parameters - default in parentheses: + S0 <- parameter() + I0 <- parameter() + beta <- parameter(0.0165) + gamma <- parameter(0.1) + + # Dimensions of arrays + dim(S0) <- 2 + dim(I0) <- 2 + dim(S) <- 2 + dim(I) <- 2 + dim(R) <- 2 + dim(n_SI) <- 2 + dim(n_IR) <- 2 + dim(p_SI) <- 2 + dim(m) <- c(2, 2) + dim(s_ij) <- c(2, 2) + dim(lambda) <- 2 +}) +``` + +In the odin code above + +```r +update(S[]) <- S[i] - n_SI[i] +``` + +becomes (approximately) + +```r +for (int i = 0; i < S_length; ++i) { + update_S[i] = S[i] + n_SI[i]; +} +``` + +so there is an implicit indexing by `i` on the LHS in that equation of the odin code, and the generated code will then produce a for loop over all values of `i`. + +## Declaring dimensions + +We can see in our first example that it is necessary to declare the dimensions of all arrays within your `odin` code using a `dim` equation. + +Dimensions can be hardcoded: +```r +dim(S) <- 2 +``` +or can be generalised so that the dimensions themselves become parameters: +```r +dim(S) <- N_age +N_age <- parameter() +``` +You may have array parameters, these can have dimensions explicitly declared +```r +dim(m) <- c(N_age, N_age) +``` +or they can also have their dimensions detected when the parameters are passed into the system - this still needs a `dim` equation where you must explicitly state the rank: +```r +dim(m) <- parameter(rank = 2) +``` + +## Summing over arrays + +Frequently, you will want to take a sum over an array, or part of an array, using `sum`. To sum over all elements of an array, use `sum()` with the name of the array you would like to sum over: + +```r +dim(x) <- 10 +x_tot <- sum(x) +``` + +If `m` is a matrix you can compute the sums over the second column by writing: + +```r +m1_tot <- sum(m[, 2]) +``` + +This partial sum approach is frequently used within implicit loops: + +```r +m_col_totals[] <- sum(m[, i]) +``` + +In our example we calculated the force of infection by summing over the two age groups +```r +lambda[] <- beta * (s_ij[i, 1] + s_ij[i, 2]) +``` +so we could have used `sum` to write this as +```r +lambda[] <- beta * sum(s_ij[i, ]) +``` +and this equation has been generalised to any number of age groups now! + +::: {.callout-note} +The same syntax applies to other functions that reduce arrays: `min`, `max` and `prod`. +::: + + +## Example: a generalised age structured SIR model {#sec-stochastic-age} + +Let's now take the two group SIR model and generalise it to `N_age` age groups ```{r} sir_age <- odin({ @@ -70,24 +200,29 @@ sir_age <- odin({ }) ``` -In the odin code above - +::: {.callout-note} +With the equations ```r -update(S[]) <- S[i] - n_SI[i] +s_ij[, ] <- m[i, j] * I[j] +lambda[] <- beta * sum(s_ij[i, ]) ``` - -becomes (approximately) - +what we are really doing is matrix multiplication, and in R this would be ```r -for (int i = 0; i < S_length; ++i) { - update_S[i] = S[i] + n_SI[i]; -} +lambda <- beta * (m %*% I) ``` +We are aiming to support matrix multiplication in future to help simplify this code. +::: + + +## Indexing -so there is an implicit indexing by `i` on the LHS in that equation of the odin code, and the generated code will then loop over all values of `i`. +We have seen in the above examples uses of 1d and 2d arrays, making use of index variables `i` and `j`. Note that while we never use these index variables on the LHS, it is implicit that on the LHS we are indexing by `i` for 1d arrays and `i` and `j` for 2d arrays! Use of these index variables on the RHS corresponds to the indexing on the LHS. +Of course you may want to use higher-dimensional arrays! We currently supporting up to 8 dimensions, with index variables `i`, `j`, `k`, `l`, `i5`, `i6`, `i7` and `i8`. + +Be careful with array equations! It's possible that in an array equation you end up going out of bounds with an array used on the RHS. This can crash the program when running. -## Another example: an age-structured SIR model with vaccination +## Example: an age-structured SIR model with vaccination Let's take the above model and additionally add some vaccination to it. @@ -171,35 +306,6 @@ new_S[, ] <- S[i, j] - n_SI[i, j] - n_S_vax[i, j] + ``` Note that in odin, an `if` always requires an `else`! -## Indexing - -We have seen in the above examples uses of 1d and 2d arrays, making use of index variables `i` and `j`. Note that while we never use these index variables on the LHS, it is implicit that on the LHS we are indexing by `i` for 1d arrays and `i` and `j` for 2d arrays! Use of these index variables on the RHS corresponds to the indexing on the LHS. - -Of course you may want to use higher-dimensional arrays! We currently supporting up to 8 dimensions, with index variables `i`, `j`, `k`, `l`, `i5`, `i6`, `i7` and `i8`. - -Be careful with array equations! It's possible that in an array equation you end up going out of bounds with an array used on the RHS. This can crash the program when running. - -## Declaring dimensions - -You must declare the dimensions of all arrays within your `odin` code using a `dim` equation. - -Dimensions can be hardcoded: -```r -dim(y) <- 10 -``` -or dimensions themselves can be parameters: -```r -dim(x) <- c(nr, 5) -nr <- parameter() -``` -You may have array parameters, these can have dimensions explicitly declared -```r -dim(m) <- c(N_age, N_age) -``` -or they can also have their dimensions detected when the parameters are passed into the system - this still needs a `dim` equation where you must explicitly state the rank: -```r -dim(m) <- parameter(rank = 2) -``` ## Semantics of random number draws @@ -217,25 +323,3 @@ x[] <- a ``` `x` will be a vector where every element is the same, the result of a *single* draw from a standard normal. - - - -## Summing over arrays -Frequently, you will want to take a sum over an array, or part of an array, using `sum`. To sum over all elements of an array, use `sum()` with the name of the array you would like to sum over: - -```r -dim(x) <- 10 -x_tot <- sum(x) -``` - -If `m` is a matrix you can compute the sums over the second column by writing: - -```r -m1_tot <- sum(m[, 2]) -``` - -This partial sum approach is frequently used within implicit loops: - -```r -m_col_totals[] <- sum(m[, i]) -``` \ No newline at end of file