source("../R/shobbs.R")
source("../R/shobbslik.R")
source("../R/shobbslikn.R")


library(optimx)
library(plyr)
library(R2admb)
library(R2jags)

## print(shobbs.lik(log(c(2,5,3,1))))
## ansR5 <- optimx(log(c(2,5,3,1)),shobbs.lik, shobbs.lg, control=list(all.methods=TRUE))

simpars <- c(logb1=0.6,logb2=1.6,logb3=1.15,logSigma=-0.8)
simfun <- function(p=simpars,t=1:12) {
  x <- exp(p[1:3])
  m <- 100.0*x[1]/(1+x[2]*10.*exp(-0.1*x[3]*t))
  rnorm(length(m),m,sd=exp(p[4]))
}

start0 <- log(c(logb1=2,logb2=5,logb3=3,logSigma=1))
nsim <- 100
## L <- list()
## for (i in 1:nsim) {
##   s <- simfun()
##   ans <- optimx(start0,fn=shobbs.lik, gr=shobbs.lg, y=s,
##                 control=list(all.methods=TRUE))
##   ## hessian ignored??
  
##   L[[i]] <- data.frame(method=paste("R",unlist(ans$method),sep="_"),
##                        .n=i,
##                        elapsed=unlist(ans$xtime),
##                        do.call(rbind,ans$par),
##                        matrix(NA,ncol=8,nrow=length(ans$par)))
## }
## L2 <- do.call(rbind,L)


file.copy("../ADMB/weed.tpl",".")
system("admb weed.tpl")
## alternately:
L2  <- rdply(nsim,
             {
               s <- simfun()
               ans <- optimx(start0,fn=shobbs.lik, gr=shobbs.lg, y=s,
                             control=list(all.methods=TRUE))
               Rvals <- data.frame(method=paste("R",unlist(ans$method),sep="_"),
                                   elapsed=unlist(ans$xtime),
                                   do.call(rbind,ans$par),
                                   matrix(NA,ncol=8,nrow=length(ans$par)))
               ###
               write_dat("weed",list(noObs=12,d=cbind(t=1:12,y=s)))
               st2 <- system.time(system("./weed"))
               ADMB_ans <- read_pars("weed")
               CI <- with(ADMB_ans,coefficients+ c(-1.96,1.96) %o% se)
               ADMBvals <- data.frame(method="ADMB",matrix(c(st2["elapsed"],coef(ADMB_ans),CI),nrow=1))
               names(ADMBvals) <- names(Rvals)
               ###
               y <- s
               t <- 1:12
               ## Note that n.iter=10000 and n.burning=1000 was not sufficient.
               st3 <- system.time(tfit_jags <- jags(model="../BUGS/weeds_bugs.txt",
                                 data=list(N=length(y),t=t, y=y),
                                 parameters.to.save=c("b1","b2","b3","sigma"),
                                 n.iter=30000,
                                 n.burnin=5000,
                                 n.thin=10,
                                 n.chains=3,
                                 inits=list(list(b1=1,b2=1,b3=1,
                                   .RNG.name="base::Wichmann-Hill", .RNG.seed=654321),
                                   list(b1=2,b2=2,b3=2,
                                        .RNG.name="base::Wichmann-Hill", .RNG.seed=321654),
                                   list(b1=0,b2=0,b3=0,
                                        .RNG.name="base::Wichmann-Hill", .RNG.seed=123456))
                                 ))
               bb <- tfit_jags$BUGSoutput
               bbs <- bb$summary[c(1:3,5),c("mean","2.5%","97.5%")] ## omit deviance
               BUGSvals <- data.frame(method="BUGS",matrix(c(st3["elapsed"],bbs[,"mean"],t(bbs[,-1])),nrow=1))
               names(BUGSvals) <- names(Rvals)
               ###
               rbind(Rvals,ADMBvals,BUGSvals)
             },
             .progress="text")



save("L2","simpars",file="weeds_simdat.RData")