@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.
\nYou'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.
\nrinit_fn <- function (s1, s2, s3, e1, e2, e3, ...) {\n population <- 1e4\n S <- round(c(s1,s2,s3)*population)\n E <- round(c(e1,e2,e3)*population)\n c(S=S,E=E)\n}\n\nrinit_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 \n## 1000 1000 2000 1000 2000 3000\nstep_fn <- function (S1, S2, S3, E1, E2, E3, r, delta.t, ...) {\n dE <- rpois(n=3,lambda=r*delta.t*c(S1,S2,S3))\n S <- c(S1,S2,S3)-dE\n E <- c(E1,E2,E3)+dE\n c(S=S,E=E)\n}\n\nstep_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 \n## 88 184 272 412 516 628\nlibrary(pomp)\nsimulate(\n t0=0, times=seq(0,5),\n rinit=rinit_fn,\n rprocess=discrete_time(step_fn,delta.t=1),\n params=c(r=0.3,s1=0.1,s2=0.1,s3=0.2,e1=0.1,e2=0.2,e3=0.3),\n format=\"data.frame\"\n)## Warning: 'rmeasure' unspecified: NAs generated.\n## time .id S1 S2 S3 E1 E2 E3\n## 1 0 1 1000 1000 2000 1000 2000 3000\n## 2 1 1 708 668 1400 1292 2332 3600\n## 3 2 1 505 483 1004 1495 2517 3996\n## 4 3 1 363 338 720 1637 2662 4280\n## 5 4 1 266 246 514 1734 2754 4486\n## 6 5 1 195 170 355 1805 2830 4645\nThis 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.
\nIn 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
\ndouble *x;\nint *i;\nIn 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]\n*(p+k)\nYou 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.
\nrinit_snip <- \"\ndouble pop = 10000;\ndouble *s = &s1;\ndouble *e = &e1;\ndouble *S = &S1;\ndouble *E = &E1;\nint i;\nfor (i = 0; i < 3; i++) {\n S[i] = nearbyint(s[i]*pop);\n E[i] = nearbyint(e[i]*pop);\n}\n\"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 <- \"\ndouble pop = 10000;\ndouble *x = &s1;\ndouble *X = &S1;\nint i, j;\nfor (i = 0; i < 2; i++) {\n for (j = 0; j < 3; j++) {\n X[3*i+j] = nearbyint(x[3*i+j]*pop);\n }\n}\n\"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:
\nstep_snip <- \"\ndouble *S = &S1;\ndouble *E = &E1;\ndouble dE[3];\nint i;\nfor (i = 0; i < 3; i++) {\n dE[i] = rpois(r*dt*S[i]);\n}\nfor (i = 0; i < 3; i++) {\n S[i] -= dE[i];\n E[i] += dE[i];\n}\n\"Putting these together to generate a simulation,
\nsimulate(\n t0=0, times=seq(0,5),\n rinit=Csnippet(rinit_snip2),\n rprocess=discrete_time(Csnippet(step_snip),delta.t=1),\n params=c(r=0.3,s1=0.1,s2=0.1,s3=0.2,e1=0.1,e2=0.2,e3=0.3),\n format=\"data.frame\",\n statenames=c(sprintf(\"S%d\",1:3),sprintf(\"E%d\",1:3)),\n paramnames=c(\"r\",sprintf(\"s%d\",1:3),sprintf(\"e%d\",1:3))\n)## Warning: 'rmeasure' unspecified: NAs generated.\n## time .id S1 S2 S3 E1 E2 E3\n## 1 0 1 1000 1000 2000 1000 2000 3000\n## 2 1 1 693 702 1411 1307 2298 3589\n## 3 2 1 493 476 993 1507 2524 4007\n## 4 3 1 337 315 673 1663 2685 4327\n## 5 4 1 238 213 476 1762 2787 4524\n## 6 5 1 163 146 343 1837 2854 4657\nNote 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.
-
|
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, 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)
} |
Beta Was this translation helpful? Give feedback.
Replies: 1 comment 1 reply
-
|
@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)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)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"
)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 In these declarations, 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 ( 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 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 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))
)Note that, in the |
Beta Was this translation helpful? Give feedback.
-
|
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^_^! |
Beta Was this translation helpful? Give feedback.
@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.