library(R2jags)
source("wildflower_BUGS_funs.R")

sessionInfo()

rawData <- read.csv("../DATA/wildflower_pods_complete.csv",
                    header=TRUE, row.names = 1)

dataForJags <- setup_wildflower_data(rawData)

set.seed(81571)

## Use this for full runs
st1 <- system.time(wfit_jags <- wildflower_BUGS_fit(dataForJags,
                                                    seed = round(runif(1)*1000000)))["elapsed"]

## Use this for quick test runs
##st1 <- system.time(wfit_jags <- wildflower_BUGS_fit(dataForJags, seed = round(runif(1)*1000000), n.iter = 2000, n.thin = 1))["elapsed"]

## Do some transformations to match parameterizations with ADMB and R.
## I will do this for the full MCMC output.  The BUGSoutput includes 3
## copies of the chains in different formats. I determined that the
## sims.matrix version seems to be what as.mcmc uses, so I'll
## transform that one.

## 1. intercept1 is the reference intercept.  Make all others relative to it.
F.int <- wfit_jags$BUGSoutput$sims.matrix[,'intercept[1]'] 
for(i in 2:5) {
  label <- paste('intercept[',i,']',sep='')
  wfit_jags$BUGSoutput$sims.matrix[,label] <- wfit_jags$BUGSoutput$sims.matrix[,label] - F.int
}

## Now it is as if we had sampled with non-adjusted data and with contrasts as done for R and ADMB

wfit_jags_mcmc <- as.mcmc(wfit_jags)

## coefs <- wfit_jags$BUGSoutput$mean ## Can't use this because it was calculated before the transformations of the posterior sample
coefs <- apply(wfit_jags$BUGSoutput$sims.matrix, 2, mean)
fixef.names <- c(paste("intercept[",1:5,"]",sep=""),"slope")
coefs.fix <- coefs[fixef.names]
vnamebase <- c("plantIntercept","plantSlope","yearIntercept")
if (any(grepl("Precision",names(coefs)))) {
  vname_ext <- "Precision"
  vpow <- -1  ## power to get from coef to variance estimate
} else {
  vname_ext <- "SD"
  vpow <- 2
}
vnames <- paste(vnamebase,vname_ext,sep="")
vmcmc <- wfit_jags_mcmc[,vnames]^vpow
colnames(vmcmc) <- paste(vnamebase,"Var",sep="")
coefs.vc <- colMeans(vmcmc)

CI <- HPDinterval(wfit_jags_mcmc)
CI.fix <- CI[fixef.names,]

CI.vc <- HPDinterval(vmcmc)

results <- list(coef=c(coefs.fix,coefs.vc),
                confint.mcmc=rbind(CI.fix,CI.vc),
                time=c(fit=as.numeric(st1),mcmc=as.numeric(st1)),
                convinfo=list(effsize=effectiveSize(wfit_jags_mcmc),
                  geweke=geweke.diag(wfit_jags_mcmc))
                )


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