improve testing

tests/issue222.R

@@ -84,4 +84,66 @@ try(
     )
 )
 
+library(pomp)
+
+simulate(
+  t0=-1/52,
+  times=seq(0,10,by=1/52),
+  params=c(
+    gamma = 26, mu = 0.02, iota = 0.1,
+    Beta = 400, Beta_sd = 0.01,
+    rho = 0.6, k = 100,
+    pop = 1e6,
+    S_0 = 25/400, I_0 = 0.001, R_0 = 375/400
+  ),
+  rprocess = euler(
+    step.fun=Csnippet(r"{
+      double dW = rgammawn(Beta_sd,dt);
+      double rate[6];
+      double trans[6];
+      rate[0] = mu*pop;
+      rate[1] = Beta*(I+iota)/pop*dW/dt;
+      rate[2] = mu;
+      rate[3] = gamma;
+      rate[4] = mu;
+      rate[5] = mu;
+      trans[0] = rate[0]*dt;
+      eeulermultinom(2,S,&rate[1],dt,&trans[1]);
+      eeulermultinom(2,I,&rate[3],dt,&trans[3]);
+      eeulermultinom(1,R,&rate[5],dt,&trans[5]);
+      // balance equations
+      S += trans[0] - trans[1] - trans[2];
+      I += trans[1] - trans[3] - trans[4];
+      R += trans[3] - trans[5];
+      C += trans[3];
+      W += (dW-dt)/Beta_sd;}"
+    ),
+    delta.t=1/365
+  ),
+  rmeasure = Csnippet(r"{
+      reports = rnbinom_mu(k,rho*C);}"
+  ),
+  rinit=Csnippet(r"{
+      double m = pop/(S_0+I_0+R_0);
+      S = m*S_0;
+      I = m*I_0;
+      R = m*R_0;
+      C = 0;
+      W = 0;}"
+  ),
+  partrans=parameter_trans(
+    logit=c("rho"),
+    log=c("gamma","mu","k","Beta","Beta_sd","iota"),
+    barycentric=c("S_0","I_0","R_0")
+  ),
+  obsnames = "reports",
+  accumvars = c("W","C"),
+  statenames = c("S","I","R","C","W"),
+  paramnames = c(
+    "pop","rho","k","gamma","mu","Beta","Beta_sd","iota",
+    "S_0","I_0","R_0"
+  )
+) |>
+  plot()
+
 dev.off()

tests/issue222.Rout.save

@@ -106,6 +106,68 @@ Error : in 'simulate': argument is missing, with no default
 Error : in 'simulate': The argument is not recognized.
 Use the 'userdata' argument to supply extra objects to basic model components. See '?userdata'.
 > 
+> library(pomp)
+> 
+> simulate(
++   t0=-1/52,
++   times=seq(0,10,by=1/52),
++   params=c(
++     gamma = 26, mu = 0.02, iota = 0.1,
++     Beta = 400, Beta_sd = 0.01,
++     rho = 0.6, k = 100,
++     pop = 1e6,
++     S_0 = 25/400, I_0 = 0.001, R_0 = 375/400
++   ),
++   rprocess = euler(
++     step.fun=Csnippet(r"{
++       double dW = rgammawn(Beta_sd,dt);
++       double rate[6];
++       double trans[6];
++       rate[0] = mu*pop;
++       rate[1] = Beta*(I+iota)/pop*dW/dt;
++       rate[2] = mu;
++       rate[3] = gamma;
++       rate[4] = mu;
++       rate[5] = mu;
++       trans[0] = rate[0]*dt;
++       eeulermultinom(2,S,&rate[1],dt,&trans[1]);
++       eeulermultinom(2,I,&rate[3],dt,&trans[3]);
++       eeulermultinom(1,R,&rate[5],dt,&trans[5]);
++       // balance equations
++       S += trans[0] - trans[1] - trans[2];
++       I += trans[1] - trans[3] - trans[4];
++       R += trans[3] - trans[5];
++       C += trans[3];
++       W += (dW-dt)/Beta_sd;}"
++     ),
++     delta.t=1/365
++   ),
++   rmeasure = Csnippet(r"{
++       reports = rnbinom_mu(k,rho*C);}"
++   ),
++   rinit=Csnippet(r"{
++       double m = pop/(S_0+I_0+R_0);
++       S = m*S_0;
++       I = m*I_0;
++       R = m*R_0;
++       C = 0;
++       W = 0;}"
++   ),
++   partrans=parameter_trans(
++     logit=c("rho"),
++     log=c("gamma","mu","k","Beta","Beta_sd","iota"),
++     barycentric=c("S_0","I_0","R_0")
++   ),
++   obsnames = "reports",
++   accumvars = c("W","C"),
++   statenames = c("S","I","R","C","W"),
++   paramnames = c(
++     "pop","rho","k","gamma","mu","Beta","Beta_sd","iota",
++     "S_0","I_0","R_0"
++   )
++ ) |>
++   plot()
+> 
 > dev.off()
 null device 
           1 

tests/mif2_index.R

@@ -7,10 +7,12 @@ gompertz() -> gompertz
 set.seed(1481104436)
 
 gompertz |>
-  mif2(Nmif=4,Np=1000,
+  mif2(
+    Nmif=4,Np=1000,
     .indices=seq.int(1000),
     rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-    cooling.fraction.50=0.5) |>
+    cooling.fraction.50=0.5
+  ) |>
   slot("indices") -> idx
 stopifnot(
   length(idx)==1000,
@@ -20,16 +22,24 @@ stopifnot(
 set.seed(962724905)
 
 gompertz |>
-  mif2(Nmif=4,Np=100,
+  mif2(
+    Nmif=4,Np=100,
     .indices=as.list(seq.int(100)),
     rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-    cooling.fraction.50=0.5) |>
+    cooling.fraction.50=0.5
+  ) |>
   slot("indices") -> idx
 stopifnot(
   length(idx)==100,
   class(idx)=="list"
 )
 
-try(mif2(gompertz,Nmif=1,Np=100,
-  .indices=1:5,rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-  cooling.fraction=0.5))
+try(
+  mif2(
+    gompertz,
+    Nmif=1,Np=100,
+    .indices=1:5,
+    rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
+    cooling.fraction.50=0.5
+  )
+)

tests/mif2_index.Rout.save

@@ -26,10 +26,12 @@ Type 'q()' to quit R.
 > set.seed(1481104436)
 > 
 > gompertz |>
-+   mif2(Nmif=4,Np=1000,
++   mif2(
++     Nmif=4,Np=1000,
 +     .indices=seq.int(1000),
 +     rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-+     cooling.fraction.50=0.5) |>
++     cooling.fraction.50=0.5
++   ) |>
 +   slot("indices") -> idx
 > stopifnot(
 +   length(idx)==1000,
@@ -39,19 +41,26 @@ Type 'q()' to quit R.
 > set.seed(962724905)
 > 
 > gompertz |>
-+   mif2(Nmif=4,Np=100,
++   mif2(
++     Nmif=4,Np=100,
 +     .indices=as.list(seq.int(100)),
 +     rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-+     cooling.fraction.50=0.5) |>
++     cooling.fraction.50=0.5
++   ) |>
 +   slot("indices") -> idx
 > stopifnot(
 +   length(idx)==100,
 +   class(idx)=="list"
 + )
 > 
-> try(mif2(gompertz,Nmif=1,Np=100,
-+   .indices=1:5,rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
-+   cooling.fraction=0.5))
-Error : in 'mif2': The argument 'cooling.fraction'is not recognized.
-Use the 'userdata' argument to supply extra objects to basic model components. See '?userdata'.
+> try(
++   mif2(
++     gompertz,
++     Nmif=1,Np=100,
++     .indices=1:5,
++     rw.sd=rw_sd(r=0.02,K=0.02,sigma=0.02),
++     cooling.fraction.50=0.5
++   )
++ )
+Error : in 'mif2': '.indices' has improper length.
 >