## This file generates one-time results for the WRITEUP.  It is somewhat slow, so I will also put the results file on svn
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)

st1 <- system.time(wfit_jags_0001boundary <- wildflower_BUGS_fit_0001boundary(dataForJags, seed = round(runif(1)*1000000)))["elapsed"]

save(wfit_jags_0001boundary, st1, file = "test_SDflatprior_101000.Rdata")

st2 <- system.time(wfit_jags_gamma_priors <- wildflower_BUGS_fit_gamma_priors(dataForJags, seed = round(runif(1)*1000000)))["elapsed"]

save(wfit_jags_gamma_priors, st2, file = "test_SDgammaprior_101000.Rdata")


st6 <- system.time(wfit_jags_0001boundary_10000<- wildflower_BUGS_fit_0001boundary(dataForJags, seed = round(runif(1)*1000000), n.chains = 3, n.iter = 11000, n.thin = 10))["elapsed"]

save(wfit_jags_0001boundary_10000, st1, file = "test_SDflatprior_10000.Rdata")

st7 <- system.time(wfit_jags_gamma_priors_10000 <- wildflower_BUGS_fit_gamma_priors(dataForJags, seed = round(runif(1)*1000000), n.chains = 3, n.iter = 11000, n.thin = 10))["elapsed"]

save(wfit_jags_gamma_priors_10000, st2, file = "test_SDgammaprior_10000.Rdata")



st3 <- system.time(wfit_jags_1000 <- wildflower_BUGS_fit(dataForJags, seed = round(runif(1)*1000000), n.chains = 3, n.iter = 2000, n.thin = 1))["elapsed"]

st4 <- system.time(wfit_jags_10000 <- wildflower_BUGS_fit(dataForJags, seed = round(runif(1)*1000000), n.chains = 3, n.iter = 11000, n.thin = 10))["elapsed"]

st5 <- system.time(wfit_jags_100000 <- wildflower_BUGS_fit(dataForJags, seed = round(runif(1)*1000000), n.chains = 3, n.iter = 101000, n.thin = 100))["elapsed"]

save(wfit_jags_1000, wfit_jags_10000, wfit_jags_100000, st3, st4, st5, file = "test_different_run_lengths.Rdata")

raw <- wfit_jags
## 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)
wfit_jags_mcmc <- as.mcmc(raw)

## 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)
coefs <- coefs[names(coefs)!="deviance"]

CI <- HPDinterval(wfit_jags_mcmc)
CI <- CI[rownames(CI)!="deviance",]
  
gresults <- list(coef=coefs,
                confint.mcmc=CI,
                time=c(fit=st1,mcmc=st1),
                convinfo=list(effsize=effectiveSize(wfit_jags_mcmc),
                  geweke=geweke.diag(wfit_jags_mcmc))
                )


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