[Pomp-commits] r1201 - in pkg: pomp pomp/R pompExamples/inst/examples

[noreply at r-forge.r-project.org]

Author: kingaa
Date: 2015-06-11 13:56:42 +0200 (Thu, 11 Jun 2015)
New Revision: 1201

Modified:
   pkg/pomp/DESCRIPTION
   pkg/pomp/R/builder.R
   pkg/pomp/R/pomp.R
   pkg/pompExamples/inst/examples/ebola.R
Log:
- fix bug with Csnippet initializer
- update ebola example to use this facility

Modified: pkg/pomp/DESCRIPTION
===================================================================
--- pkg/pomp/DESCRIPTION	2015-06-11 11:56:30 UTC (rev 1200)
+++ pkg/pomp/DESCRIPTION	2015-06-11 11:56:42 UTC (rev 1201)
@@ -2,7 +2,7 @@
 Type: Package
 Title: Statistical Inference for Partially Observed Markov Processes
 Version: 0.67-1
-Date: 2015-06-10
+Date: 2015-06-11
 Authors at R: c(person(given=c("Aaron","A."),family="King",
 		role=c("aut","cre"),email="kingaa at umich.edu"),
 	  person(given=c("Edward","L."),family="Ionides",role=c("aut")),

Modified: pkg/pomp/R/builder.R
===================================================================
--- pkg/pomp/R/builder.R	2015-06-11 11:56:30 UTC (rev 1200)
+++ pkg/pomp/R/builder.R	2015-06-11 11:56:42 UTC (rev 1201)
@@ -338,7 +338,7 @@
                dmeasure= "__pomp_dmeasure",
                step.fn="__pomp_stepfn",
                skeleton="__pomp_skelfn",
-               skeleton="__pomp_initializer",
+               initializer="__pomp_initializer",
                fromEstimationScale="__pomp_par_trans",
                toEstimationScale="__pomp_par_untrans",
                rprior="__pomp_rprior",

Modified: pkg/pomp/R/pomp.R
===================================================================
--- pkg/pomp/R/pomp.R	2015-06-11 11:56:30 UTC (rev 1200)
+++ pkg/pomp/R/pomp.R	2015-06-11 11:56:42 UTC (rev 1201)
@@ -144,22 +144,7 @@
   
   ## handle initializer
   default.init <- missing(initializer) || (is(initializer,"pomp.fun") && initializer at mode == pompfunmode$undef)
-
-  if (default.init) {
-    initializer <- pomp.fun()
-  } else {
-    initializer <- pomp.fun(
-                            f=initializer,
-                            PACKAGE=PACKAGE,
-                            proto=quote(initializer(params,t0,...)),
-                            slotname="initializer",
-                            libname=libname,
-                            statenames=statenames,
-                            paramnames=paramnames,
-                            obsnames=obsnames,
-                            covarnames=covarnames
-                            )
-  }
+  if (default.init) initializer <- pomp.fun()
   
   ## default rprocess & dprocess
   if (missing(rprocess))
@@ -220,6 +205,21 @@
     libname <- ''
   }
   
+  ## handle initializer
+  if (!default.init) {
+    initializer <- pomp.fun(
+                            f=initializer,
+                            PACKAGE=PACKAGE,
+                            proto=quote(initializer(params,t0,...)),
+                            slotname="initializer",
+                            libname=libname,
+                            statenames=statenames,
+                            paramnames=paramnames,
+                            obsnames=obsnames,
+                            covarnames=covarnames
+                            )
+  }
+  
   ## handle rprocess
   if (is(rprocess,"pompPlugin")) {
     rprocess <- plugin.handler(

Modified: pkg/pompExamples/inst/examples/ebola.R
===================================================================
--- pkg/pompExamples/inst/examples/ebola.R	2015-06-11 11:56:30 UTC (rev 1200)
+++ pkg/pompExamples/inst/examples/ebola.R	2015-06-11 11:56:42 UTC (rev 1201)
@@ -163,7 +163,19 @@
   DN_IR = gamma * I;
 ')
 
+initlzr <- Csnippet("
+  int j;
+  double m, *E;
+  E = &E1;
+  m = N/(S_0+E_0+I_0+R_0);
+  S = nearbyint(m*S_0);
+  I = nearbyint(m*I_0);
+  R = nearbyint(m*R_0);
+  m = nearbyint(m*E_0/nstageE);
+  for (j = 0; j < nstageE; j++) E[j] = m;
+")
 
+
 ebolaModel <- function (country=c("Guinea", "SierraLeone", "Liberia", "WestAfrica"),
                         data = NULL,
                         timestep = 0.01, nstageE = 3L,
@@ -219,26 +231,18 @@
        params=theta,
        globals=globs,
        obsnames=c("cases","deaths"),
-       statenames=c("S","E1","I","R","N_EI","N_IR"),
+       statenames=c("S",sprintf("E%d",1:nstageE),"I","R","N_EI","N_IR"),
        zeronames=if (type=="raw") c("N_EI","N_IR") else character(0),
        paramnames=c("N","R0","alpha","gamma","rho","k","cfr",
          "S_0","E_0","I_0","R_0"),
-       nstageE=nstageE,
        dmeasure=if (least.sq) dObsLS else dObs,
        rmeasure=if (least.sq) rObsLS else rObs,
        rprocess=discrete.time.sim(step.fun=rSim,delta.t=timestep),
        skeleton=skel,
        skeleton.type="vectorfield",
-       parameter.transform=partrans,
-       parameter.inv.transform=paruntrans,
-       initializer=function (params, t0, nstageE, ...) {
-         all.state.names <- c("S",paste0("E",1:nstageE),"I","R","N_EI","N_IR")
-         comp.names <- c("S",paste0("E",1:nstageE),"I","R")
-         x0 <- setNames(numeric(length(all.state.names)),all.state.names)
-         frac <- c(params["S_0"],rep(params["E_0"]/nstageE,nstageE),params["I_0"],params["R_0"])
-         x0[comp.names] <- round(params["N"]*frac/sum(frac))
-         x0
-       }
+       fromEstimationScale=partrans,
+       toEstimationScale=paruntrans,
+       initializer=initlzr
        ) -> po
 }