re-run polio; update everything else

bobdouma · Apr 20, 2016 · a55ba70 · a55ba70
1 parent c741011
commit a55ba70
Show file tree

Hide file tree

Showing 18 changed files with 440 additions and 492 deletions.
diff --git a/measles/measles.R b/measles/measles.R
@@ -182,11 +182,11 @@ library(doParallel)
 registerDoParallel()
 
 set.seed(998468235L,kind="L'Ecuyer")
-mcopts <- list(preschedule=FALSE,set.seed=TRUE)
 
 foreach(i=1:4,
         .packages="pomp",
-        .options.multicore=mcopts) %dopar% {
+        .options.multicore=list(set.seed=TRUE)
+        ) %dopar% {
   pfilter(m1,Np=10000,params=theta)
 } -> pfs
 logmeanexp(sapply(pfs,logLik),se=TRUE)

diff --git a/measles/measles.Rmd b/measles/measles.Rmd
@@ -996,11 +996,11 @@ library(doParallel)
 registerDoParallel()
 
 set.seed(998468235L,kind="L'Ecuyer")
-mcopts <- list(preschedule=FALSE,set.seed=TRUE)
 
 foreach(i=1:4,
         .packages="pomp",
-        .options.multicore=mcopts) %dopar% {
+        .options.multicore=list(set.seed=TRUE)
+        ) %dopar% {
   pfilter(m1,Np=10000,params=theta)
 } -> pfs
 logmeanexp(sapply(pfs,logLik),se=TRUE)

diff --git a/measles/measles.html b/measles/measles.html
@@ -825,11 +825,11 @@ <h4>Constructing the <code>pomp</code> object</h4>
 registerDoParallel()
 
 set.seed(998468235L,kind=&quot;L'Ecuyer&quot;)
-mcopts &lt;- list(preschedule=FALSE,set.seed=TRUE)
 
 foreach(i=1:4,
         .packages=&quot;pomp&quot;,
-        .options.multicore=mcopts) %dopar% {
+        .options.multicore=list(set.seed=TRUE)
+        ) %dopar% {
   pfilter(m1,Np=10000,params=theta)
 } -&gt; pfs
 logmeanexp(sapply(pfs,logLik),se=TRUE)</code></pre>

diff --git a/mif/mif.R b/mif/mif.R
@@ -1,6 +1,7 @@
 library(pomp)
 options(stringsAsFactors=FALSE)
-stopifnot(packageVersion("pomp")>="1.4")
+stopifnot(packageVersion("pomp")>="1.4.7")
+set.seed(557976883)
 bsflu_data <- read.table("http://kingaa.github.io/short-course/stochsim/bsflu_data.txt")
 statenames <- c("S","I","R1","R2")
 paramnames <- c("Beta","mu_I","rho","mu_R1","mu_R2")
@@ -67,8 +68,6 @@ library(foreach)
 library(doParallel)
 
 registerDoParallel()
-
-set.seed(2036049659,kind="L'Ecuyer")
 stew(file="pf.rda",{
   t_pf <- system.time(
     pf <- foreach(i=1:10,.packages='pomp',

diff --git a/mif/mif.Rmd b/mif/mif.Rmd
@@ -51,7 +51,8 @@ options(keep.source=TRUE,encoding="UTF-8")
 ```{r prelims,include=FALSE,purl=TRUE,cache=FALSE}
 library(pomp)
 options(stringsAsFactors=FALSE)
-stopifnot(packageVersion("pomp")>="1.4")
+stopifnot(packageVersion("pomp")>="1.4.7")
+set.seed(557976883)
 ```
 
 ----------------------------
@@ -383,7 +384,7 @@ Let's treat $\mu_{R_1}$ and  $\mu_{R_2}$ as known, and fix these parameters at t
 We will estimate $\beta$, $\mu_I$, and $\rho$.
 
 It will be helpful to parallelize most of the computations.
-Most machines nowadays have multiple cores and using this computational capacity involves only 
+Most machines nowadays have multiple cores and using this computational capacity involves only:
 
 i. the following lines of code to let R know you plan to use multiple processors; 
 i. the parallel for loop provided by the **foreach** package; and
@@ -394,14 +395,11 @@ library(foreach)
 library(doParallel)
 
 registerDoParallel()
-
-set.seed(2036049659,kind="L'Ecuyer")
 ```
 
 The first two lines load the **foreach** and **doParallel** packages, the latter being a "backend" for the **foreach** package.
 The next line tells **foreach** that we will use the **doParallel** backend.
 By default, **R** will guess how many concurrent **R** processes to run on this single multicore machine.
-The final line sets up a parallel random number generator and fixes its seed.
 
 ### Running a particle filter.
 

diff --git a/mif/mif.html b/mif/mif.html
diff --git a/parest/parest.R b/parest/parest.R
@@ -1,3 +1,6 @@
+library(pomp)
+stopifnot(packageVersion("pomp")>"1.4.5")
+set.seed(1173489184)
 f <- seq(0,1,length=100)
 R0 <- -log(1-f)/f
 plot(f~R0,type='l',xlab=expression(R[0]),ylab="fraction infected",bty='l')
@@ -161,64 +164,8 @@ plot(I~time,data=out,type='l',log='y')
 plot(R~time,data=out,type='l',log='y')
 plot(I~S,data=out,log='xy',pch='.',cex=0.5)
 par(op) 
-rain <- read.csv(paste0(baseurl,"data/dacca_rainfall.csv"))
-rain$time <- with(rain,year+(month-1)/12)
-
-plot(rainfall~time,data=rain,type='l')
-
-interpol <- with(rain,approxfun(time,rainfall,rule=2,method='constant'))
-
-data.frame(time=seq(from=1920,to=1930,by=1/365)) -> smoothed.rain
-smoothed.rain$rainfall <- sapply(smoothed.rain$time,interpol)
-
-plot(rainfall~time,data=rain,col='black',cex=2,pch=16,log='')
-lines(rainfall~time,data=smoothed.rain,col='red')
-
-
-rain.sir.model <- function (t, x, params) {
-  a <- params["a"]
-  b <- params["b"]
-  mu <- params["mu"]
-  gamma <- params["gamma"]
-  R <- interpol(t)
-  beta <- a*R/(b+R)
-  dSdt <- mu*(1-x[1])-beta*x[1]*x[2]
-  dIdt <- beta*x[1]*x[2]-(mu+gamma)*x[2]
-  dRdt <- gamma*x[2]-mu*x[3]
-  list(c(dSdt,dIdt,dRdt))
-}
-
-params <- c(a=500,b=50,mu=1/50,gamma=365/13)
-xstart <- c(S=0.07,I=0.00039,R=0.92961)
-times <- seq(from=1920,to=1930,by=7/365)
-out <- as.data.frame(
-                     ode(
-                         func=rain.sir.model,
-                         y=xstart,
-                         times=times,
-                         parms=params
-                         )
-                     )
-
-op <- par(fig=c(0,1,0,1),mfrow=c(2,2),
-          mar=c(3,3,1,1),mgp=c(2,1,0))
-plot(S~time,data=out,type='l',log='y')
-plot(I~time,data=out,type='l',log='y')
-plot(R~time,data=out,type='l',log='y')
-plot(I~S,data=out,log='xy',pch='.',cex=1)
-par(op) 
-
-
-fit <- loess(log1p(rainfall)~time,data=rain,span=0.05)
-
-data.frame(time=seq(from=1920,to=1930,by=1/365)) -> smoothed.rain
-smoothed.rain$rainfall <- expm1(predict(fit,newdata=smoothed.rain))
-
-plot(rainfall~time,data=rain,col='black',cex=2,pch=16,log='')
-lines(rainfall~time,data=smoothed.rain,col='red')
-
-url <- paste0(baseurl,"data/niamey.csv")
-niamey <- read.csv(url,comment.char="#")
+# niamey <- read.csv("http://kingaa.github.io/parest/niamey.csv",comment.char="#")
+niamey <- read.csv("niamey.csv",comment.char="#")
 plot(measles~biweek,data=niamey,type='n')
 lines(measles~biweek,data=subset(niamey,community=="A"),col=1)
 lines(measles~biweek,data=subset(niamey,community=="B"),col=2)
@@ -286,7 +233,7 @@ f <- function (beta, X.0) {
 grid <- expand.grid(beta=seq(from=100,to=300,length=50),
                     X.0=seq(from=4000,to=20000,length=50))
 grid$SSE <- with(grid,mapply(f,beta,X.0))
-require(lattice)
+library(lattice)
 contourplot(sqrt(SSE)~beta+X.0,data=grid,cuts=30)
 ## ?optim
 dat <- subset(niamey,community=="A")

diff --git a/pfilter/pfilter.R b/pfilter/pfilter.R
@@ -6,7 +6,7 @@ library(grid)
 library(plyr)
 library(reshape2)
 options(stringsAsFactors=FALSE)
-set.seed(594709947L)
+set.seed(594709947)
 base_url <- "http://kingaa.github.io/short-course/"
 url <- paste0(base_url,"stochsim/bsflu_data.txt")
 bsflu <- read.table(url)

diff --git a/pfilter/pfilter.Rmd b/pfilter/pfilter.Rmd
@@ -69,7 +69,7 @@ library(grid)
 library(plyr)
 library(reshape2)
 options(stringsAsFactors=FALSE)
-set.seed(594709947L)
+set.seed(594709947)
 ```
 
 ## Objectives

diff --git a/pfilter/pfilter.html b/pfilter/pfilter.html
diff --git a/polio/global_search.rda b/polio/global_search.rda
diff --git a/polio/local_search.rda b/polio/local_search.rda
diff --git a/polio/polio.R b/polio/polio.R
@@ -189,6 +189,7 @@ pf1 %>%
 
 ## ----local_search--------------------------------------------------------
 stew(file="local_search.rda",{
+  w1 <- getDoParWorkers()
   t1 <- system.time({
     m1 <- foreach(i=1:90,
                   .packages='pomp', .combine=rbind,
@@ -250,6 +251,7 @@ polio_box <- rbind(
 
 ## ----global_search-------------------------------------------------------
 stew(file="global_search.rda",{
+  w2 <- getDoParWorkers()
   t2 <- system.time({
     m2 <- foreach(i=1:400,.packages='pomp',.combine=rbind,
                   .options.multicore=list(set.seed=TRUE)
@@ -310,37 +312,37 @@ library(reshape2)
 library(magrittr)
 
 bake(file="profile_rho.rds",{
-    params %>% 
-        subset(loglik>max(loglik)-20,
-               select=-c(loglik,loglik.se,rho)) %>% 
-        melt(id=NULL) %>% 
-        daply(~variable,function(x)range(x$value)) -> box
-    
-    starts <- profileDesign(rho=seq(0.01,0.025,length=30),
-                            lower=box[,1],upper=box[,2],
-                            nprof=10)
-
-    foreach(start=iter(starts,"row"),
-            .combine=rbind,
-            .options.multicore=list(set.seed=TRUE),
-            .options.mpi=list(seed=290860873,chunkSize=1)
-            ) %dopar% {
-                mf <- mif2(polio,
-                           start=unlist(start),
-                           Np=2000,
-                           Nmif=300,
-                           cooling.type="geometric",
-                           cooling.fraction.50=0.5,
-                           transform=TRUE,
-                           rw.sd=rw.sd(
-                               b1=0.02, b2=0.02, b3=0.02, b4=0.02, b5=0.02, b6=0.02,
-                               psi=0.02, tau=0.02, sigma_dem=0.02, sigma_env=0.02,
-                               IO_0=ivp(0.2), SO_0=ivp(0.2)
-                           ))
-                mf <- mif2(mf,Np=5000,Nmif=100,cooling.fraction.50=0.1)
-                ll <- logmeanexp(replicate(10,logLik(pfilter(mf,Np=5000))),se=TRUE)
-                data.frame(as.list(coef(mf)),loglik=ll[1],loglik.se=ll[2])
-            }
+  params %>% 
+    subset(loglik>max(loglik)-20,
+           select=-c(loglik,loglik.se,rho)) %>% 
+    melt(id=NULL) %>% 
+    daply(~variable,function(x)range(x$value)) -> box
+
+  starts <- profileDesign(rho=seq(0.01,0.025,length=30),
+                          lower=box[,1],upper=box[,2],
+                          nprof=10)
+  
+  foreach(start=iter(starts,"row"),
+          .combine=rbind,
+          .options.multicore=list(set.seed=TRUE),
+          .options.mpi=list(seed=290860873,chunkSize=1)
+  ) %dopar% {
+    mf <- mif2(polio,
+               start=unlist(start),
+               Np=2000,
+               Nmif=300,
+               cooling.type="geometric",
+               cooling.fraction.50=0.5,
+               transform=TRUE,
+               rw.sd=rw.sd(
+                 b1=0.02, b2=0.02, b3=0.02, b4=0.02, b5=0.02, b6=0.02,
+                 psi=0.02, tau=0.02, sigma_dem=0.02, sigma_env=0.02,
+                 IO_0=ivp(0.2), SO_0=ivp(0.2)
+               ))
+    mf <- mif2(mf,Np=5000,Nmif=100,cooling.fraction.50=0.1)
+    ll <- logmeanexp(replicate(10,logLik(pfilter(mf,Np=5000))),se=TRUE)
+    data.frame(as.list(coef(mf)),loglik=ll[1],loglik.se=ll[2])
+  }
 }) -> m3
 
 ## ----save_profile_rho----------------------------------------------------
@@ -375,7 +377,7 @@ params %>%
 
 bake(file="sims.rds",seed=398906785,
      simulate(polio,nsim=2000,as.data.frame=TRUE,include.data=TRUE)
-     ) -> sims
+) -> sims
 ddply(sims,~sim,summarize,zeros=sum(cases==0)) -> num_zeros
 
 num_zeros %>%
@@ -413,7 +415,7 @@ print(pl2,vp=viewport(x=0.8,y=0.75,height=0.4,width=0.2))
 sims %>%
   subset(sim!="data") %>%
   summarize(imports=coef(polio,"psi")*mean(SO+SB1+SB2+SB3+SB4+SB5+SB6)/12
-            ) -> imports
+  ) -> imports
 
 ## ----check_class_m3,eval=F-----------------------------------------------
 ## class(m3)