How to use arrays when constructing an age-structured model

[Chaofan13]

Currently, I am attempting to construct an age-structured model and I have encountered a technical challenge. Typically, outside the pomp framework, I would organize the population state as an array, for example, S <- c(S1, S2, S3, S4). However, when translating this into the rprocess within pomp, I find that I cannot directly use S but must specify each state individually, such as S1, S2, S3, and S4. Consequently, my rprocess function ends up with an extensive list of individual state variables. How can I handle this within the pomp framework in C snippets and R？

SIR_step <- function (
  S1, S2, S3, S4, S5,
  E1, E2, E3, E4, E5,
  I1, I2, I3, I4, I5,
  V1, V2, V3, V4, V5,
  R1, R2, R3, R4, R5,
  inci1, inci2, inci3, inci4, inci5,
  no.vac1, no.vac2, no.vac3, no.vac4, no.vac5,
  DALY1, DALY2, DALY3, DALY4, DALY5,
  Death1, Death2, Death3, Death4, Death5,
  cost_d1, cost_d2, cost_d3, cost_d4, cost_d5,
  cost_indi1, cost_indi2, cost_indi3, cost_indi4, cost_indi5,
  b1, b2, b3, b4, b5,
  seas_m, ...
) {
  vac_coverage <-  0.2
  birth_rate <- 0.2
  death_r <- rep(0.2,5)
  age_migration_rates <- rep(0.2,5)
  vac_eff <- 0.20
  waifw <- WAIFW_matrix(as.numeric(c(b1,b2,b3,b4,b5)))

  alpha <- 0.1
  phi <- 0.1

  S <- c(S1,S2,S3,S4,S5)
  E <- c(E1,E2,E3,E4,E5)
  R <- c(R1,R2,R3,E4,R5)
  I <- c(I1,I2,I3,I4,I5)
  V <- c(V1,V2,V3,V4,V5)

  E <- E %*% aging_matrix(age_migration_rates,death_r)
  I <- I %*% aging_matrix(age_migration_rates,death_r)
  R <- R %*% aging_matrix(age_migration_rates,death_r)
  V <- V %*% aging_matrix(age_migration_rates,death_r)

  N <- S + I + R + E + V
  Birth <- c(sum(N)*birth_rate/52,0,0,0,0)

  S <- S + Birth*(1 - vac_coverage*vac_eff)
  V <- V + Birth*vac_coverage*vac_eff

  dS <- numeric(num_groups)
  dE <- numeric(num_groups)
  dI <- numeric(num_groups)
  dR <- numeric(num_groups)
  beta_effective <- waifw %*% t(I / N)

  inci <- rpois(5,(1 - exp(-beta_effective * 1/52))*S)
  dE <- inci - rpois(5, (1 - exp(-alpha * 1/52))*E)
  dI <- rpois(5, (1 - exp(-alpha * 1/52))*E) - rpois(5, (1 - exp(-phi * 1/52))*I)
  dR <- rpois(5, (1 - exp(-phi * 1/52))*I)

  S <- S + inci
  E <- E + dE
  I <- I + dI
  R <- R + dR
  c(S=S,E=E,I=I,R=R,V=V)
}

rinit_SIR <- function (
  s1, s2, s3, s4, s5,
  repo, e, ...
) {
  population <- c(200,400,500,600,700)
  incidence <- (1000 * repo / 100000) * sum(population)
  age_inci_propotion <- c(0.25, 0.25, 0.25, 0.15, 0.05)
  E <- incidence * age_inci_propotion * e
  I <- incidence* age_inci_propotion
  V <- rep(0, 5)
  S <- (population - E - I - V) * c(s1,s2,s3,s4,s5)
  R <- population - S - E - I - V
  c(S=S,E=E,I=I,R=R,V=V)
}

[kingaa]

@Chaofan13, your example very well demonstrates the inconvenience of having to treat every state variable, parameter, observable, and covariate by name when you supply basic model components to pomp as R functions. This illustrates one more advantage of the C snippet interface.

You've given us a lot to work with here, but it's a bit too much really. I'll illustrate how to write a C snippet to accomplish something similar, but simpler. Let's suppose we want to translate the following basic model components—which are written as R functions—into C snippets.

rinit_fn <- function (s1, s2, s3, e1, e2, e3, ...) {
  population <- 1e4
  S <- round(c(s1,s2,s3)*population)
  E <- round(c(e1,e2,e3)*population)
  c(S=S,E=E)
}

rinit_fn(s1=0.1,s2=0.1,s3=0.2,e1=0.1,e2=0.2,e3=0.3)

##   S1   S2   S3   E1   E2   E3 
## 1000 1000 2000 1000 2000 3000

step_fn <- function (S1, S2, S3, E1, E2, E3, r, delta.t, ...) {
  dE <- rpois(n=3,lambda=r*delta.t*c(S1,S2,S3))
  S <- c(S1,S2,S3)-dE
  E <- c(E1,E2,E3)+dE
  c(S=S,E=E)
}

step_fn(S1=100,S2=200,S3=300,E1=400,E2=500,E3=600,r=0.1,delta.t=1)

##  S1  S2  S3  E1  E2  E3 
##  88 184 272 412 516 628

library(pomp)
simulate(
  t0=0, times=seq(0,5),
  rinit=rinit_fn,
  rprocess=discrete_time(step_fn,delta.t=1),
  params=c(r=0.3,s1=0.1,s2=0.1,s3=0.2,e1=0.1,e2=0.2,e3=0.3),
  format="data.frame"
)

## Warning: 'rmeasure' unspecified: NAs generated.

##   time .id   S1   S2   S3   E1   E2   E3
## 1    0   1 1000 1000 2000 1000 2000 3000
## 2    1   1  708  668 1400 1292 2332 3600
## 3    2   1  505  483 1004 1495 2517 3996
## 4    3   1  363  338  720 1637 2662 4280
## 5    4   1  266  246  514 1734 2754 4486
## 6    5   1  195  170  355 1805 2830 4645

This is a very simple stochastic dynamical system. It's not very interesting, but hopefully you can see that it illustrates the issue you ask about.

In C, one can easily work with arrays, provided their entries are stored in contiguous memory locations. This is accomplished using pointers and pointer arithmetic. We need to know what data type is stored in the array and we need a pointer to the first element in the array. Both of these things are furnished by each of the declarations

double *x;
int *i;

In these declarations, x and y are declared to be pointers to arrays of floating-point and integer numbers, respectively.

The C language defines pointer arithmetic in such a way that one can easily access any given element. In addition to the usual arithmetic operations, pointer arithmetic has the "address" and "dereference" operators. The address operator (&) extracts a pointer to the given variable. Thus, if x is a variable, &x is the address of x in the computer's memory. The dereference operator (*) is the inverse of &: if p is a pointer to some memory (i.e., an address), then *p is the number stored at that memory location. Moreover, if p is a pointer to some memory location, then for any integer k, p+k points to the location k units farther along in the memory. In particular, in C, the following expressions are equivalent, and each returns the kth element of the array that begins at memory location p:

p[k]
*(p+k)

You can read more about pointers here, for example. Now, when a pomp C snippet is executed, the named parameters, state variables, etc. are provided by name. Assuming that you have stored variables belonging to an array contiguously and in order, you can access them from the C snippet using pointers. In our simplified example, we might do something like the following.

rinit_snip <- "
double pop = 10000;
double *s = &s1;
double *e = &e1;
double *S = &S1;
double *E = &E1;
int i;
for (i = 0; i < 3; i++) {
  S[i] = nearbyint(s[i]*pop);
  E[i] = nearbyint(e[i]*pop);
}
"

More generally, we can use multi-dimensional arrays. Suppose for instance that we wish to treat the S and E variables, respectively, as the first and second columns of a 2-D array. Then the following would work for an rinit C snippet:

rinit_snip2 <- "
double pop = 10000;
double *x = &s1;
double *X = &S1;
int i, j;
for (i = 0; i < 2; i++) {
  for (j = 0; j < 3; j++) {
    X[3*i+j] = nearbyint(x[3*i+j]*pop);
  }
}
"

Note that here I have assumed that the 2-D arrays X and x are in column-major order. This is the usual convention in scientific computing. When this C snippet is executed, it will place the S and I variables in their proper places, but, very importantly, it assumes that the parameters s1,e1,... are stored in a very particular order. It would be your responsibility, as the user, to make sure that you put the parameters in their proper locations at all times!

A corresponding rprocess snippet might look like this:

step_snip <- "
double *S = &S1;
double *E = &E1;
double dE[3];
int i;
for (i = 0; i < 3; i++) {
  dE[i] = rpois(r*dt*S[i]);
}
for (i = 0; i < 3; i++) {
  S[i] -= dE[i];
  E[i] += dE[i];
}
"

Putting these together to generate a simulation,

simulate(
  t0=0, times=seq(0,5),
  rinit=Csnippet(rinit_snip2),
  rprocess=discrete_time(Csnippet(step_snip),delta.t=1),
  params=c(r=0.3,s1=0.1,s2=0.1,s3=0.2,e1=0.1,e2=0.2,e3=0.3),
  format="data.frame",
  statenames=c(sprintf("S%d",1:3),sprintf("E%d",1:3)),
  paramnames=c("r",sprintf("s%d",1:3),sprintf("e%d",1:3))
)

## Warning: 'rmeasure' unspecified: NAs generated.

##   time .id   S1   S2   S3   E1   E2   E3
## 1    0   1 1000 1000 2000 1000 2000 3000
## 2    1   1  693  702 1411 1307 2298 3589
## 3    2   1  493  476  993 1507 2524 4007
## 4    3   1  337  315  673 1663 2685 4327
## 5    4   1  238  213  476 1762 2787 4524
## 6    5   1  163  146  343 1837 2854 4657

Note that, in the params argument, we have been careful to put the parameters in the order that the C snippets expect! Having done that, however, our rinit snippet ensures that the state variables are in the order that the rprocess snippet expects.

[Chaofan13]

Dear Dr. King

Thank you so much for your comprehensive and patient explanations, they have been incredibly enlightening and helpful, especially for someone like myself who is new to this field. Your approach not only clarified my doubts but also deepened my understanding of C++. I feel so fortunate to have had the opportunity to learn from you.

Thank you once again for your invaluable guidance^_^!