# Temporary solution 
# Once: source(paste0(ddmore:::DDMORE.checkConfiguration(),"/BlackBoxWinBugs/ddmorewinbugs.R"))

mdlfile<-"Simeoni_PAGE_Estimation_BUGS.mdl"

BUGS<-run.winbugs(mdlfile)
#BUGSSO<-run.winbugs(mdlfile, subfolder="WINBUGSSO")
#BUGSCoda<-run.winbugs(mdlfile, subfolder="WINBUGSCoda", LoadCoda=TRUE)

getPopulationParameters(BUGS)$Bayesian
getPopulationParameters(BUGS, what="estimates")$Bayesian # o other options,  what= precisions, intervals, all
# it is possible also to omitt $Bayesian
getIndividualParameters(BUGS)




######################  Example for coda (possible instructions   #################

##to read only one chain
#coda_out <- read.coda(output.file="output1.txt", index.file="outputIndex.txt")

##to read more chains (3 in this case) use read.coda.interactive()
##coda_out <- read.coda.interactive()
## Enter in the order:
#  #outputIndex.txt
#  #output1.txt  
#  #output2.txt  
#  #output2.txt

## summary statistics
#summary(coda_out)

## trace and density plots
#plot(coda_out)

## gelman-rubin diagnostics (if at least 2 chains)
##gelman.diag(coda_out)
##gelman.plot(coda_out)


## Marginal densities correlation for each chain

##install.packages("IDPmisc")
#library(IDPmisc)

#panel.hist <- function(x, ...)
#{
#  usr <- par("usr"); on.exit(par(usr))
#  par(usr = c(usr[1:2], 0, 1.5) )
#  h <- hist(x, plot = FALSE)
#  breaks <- h$breaks; nB <- length(breaks)
#  y <- h$counts; y <- y/max(y)
#  rect(breaks[-nB], 0, breaks[-1], y, col="blue4", ...)
#}

#panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
#{
#  usr <- par("usr"); on.exit(par(usr))
#  par(usr = c(0, 1, 0, 1))
#  r <- abs(cor(x, y, method = "spearman"))
#  txt <- format(c(r, 0.123456789), digits=digits)[1]
#  txt <- paste(prefix, txt, sep="")
#  if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
#  text(0.5, 0.5, txt, cex = cex * r)
#}

#betterPairs <- function(YourData){
#  return(pairs(YourData, lower.panel=function(...) {par(new=TRUE);ipanel.smooth(...)}, diag.panel=panel.hist, upper.panel=panel.cor))
#}


#betterPairs(data.frame(coda_out$output1.txt))
#betterPairs(data.frame(coda_out$output2.txt))
#betterPairs(data.frame(coda_out$output3.txt))
