# rm(list=ls()) ## Clear workspace. NOT in vignette
runif(1) # to force .Random.seed to be created (JN 120814)
library(rjags)
library(R2jags)
tdata<-read.table("../DATA/weeds.dat",header=TRUE)
t<-tdata$t
y<-tdata$y
rm(tdata) # cleanup

## try to drop deviance elements from coef or confidence intervals matrices,
##  but otherwise leave the object unchanged
dropdev <- function(x) {
  if (length(dim(x))==2) {
    if ("deviance" %in% rownames(x)) {
      x[rownames(x)!="deviance",]
    } else if ("deviance" %in% colnames(x)) {
      x[,colnames(x)!="deviance"]
    } else x
  } else {
    if (is.null(names(x))) x else x[names(x)!="deviance"]
  }
}



st1<-system.time(tfit_jags <- jags(model="../BUGS/weeds_bugs.txt",
      data=list(N=length(y),t=t, y=y),
      parameters.to.save=c("logb1","logb2","logb3","logsigma"),
      n.iter=10000,
      n.burnin=1000,
      n.thin=10,
      n.chains=3,
      inits=list(list(logb1=1,logb2=1,logb3=1,logsigma=1,
            .RNG.name="base::Wichmann-Hill", .RNG.seed=654321),
            list(logb1=-1,logb2=-1,logb3=-1,logsigma=-1,
            .RNG.name="base::Wichmann-Hill", .RNG.seed=321654),
            list(logb1=0,logb2=0,logb3=0,logsigma=0,
            .RNG.name="base::Wichmann-Hill", .RNG.seed=123456))))["elapsed"]
# wpred <- .... stuff to create predicted model
cat("output of jags() run\n")
tfit_jags
tfit_jags_b <- tfit_jags$BUGSoutput
library(emdbook) ## for as.mcmc.bugs
 tfit_jags_mcmc <- as.mcmc.bugs(tfit_jags_b)
library(coda)
# print(xyplot(tfit_jags_mcmc))
# x11()
# print(densityplot(tfit_jags_mcmc))
results <- list(obj=NULL,
                coef=dropdev(unlist(tfit_jags$BUGSoutput$mean)),
                confint.mcmc=dropdev(HPDinterval(lump.mcmc.list(tfit_jags_mcmc))),
                time=c(fit=st1,mcmc=st1),
                terminfo=NULL,
                convinfo=list(effsize=dropdev(effectiveSize(tfit_jags_mcmc)),
                  gelman=dropdev(gelman.diag(tfit_jags_mcmc)$psrf[,"Point est."]))
)

save("results",file="fit.RData")