library(siar)
library("BRugs")
library("R2WinBUGS")
################################
#  DEFINE MODEL PARAMETERS
################################
p.contrib = c(0.7,0.2,0.1)  # vector of proportions following Jackson

N = 10 # this is the number of data points to generate
num.prey = length(p.contrib)
num.iso = 2

sigmaW = 0.5 # following Jackson
sigmaResid = 0 # this can be anything, Jackson used (0,1)

u = matrix(0,num.prey,num.iso)  # This is the matrix of means
u[,1] = c(-10,0,10) # means for isotope 1
u[,2] = c(0, 17.32, 0) # means for isotope 2
u = u
sources.sd = matrix(sigmaW,num.prey,num.iso)
sources.dat = cbind(u[,1],sources.sd[,1],u[,2],sources.sd[,2])

sources = (sourcesdemo[1:3,])
sources[,2:5] = sources.dat
nameVec =c("Min","Q1","Median","Mean","Q3","Max","Mode")
outputNames=c(paste(nameVec,".1",sep=""),paste(nameVec,".2",sep=""),paste(nameVec,".3",sep=""),paste(nameVec,".4",sep=""))
numSims = 1000
summaryResults = matrix(0, numSims, 28)
colnames(summaryResults)=outputNames

for(ii in 1:numSims) {
  # X is the random data matrix for simulations
  X = matrix(0, N, num.iso)
  Z = matrix(rnorm(N*num.iso,0,1),N, num.iso)
  for(i in 1:N) {
    #for(j in 1:num.iso) {
      # loop over isotopes, draw a random source value
    #  X[i,j] = rnorm(num.prey,u[,j],rep(sigmaW,num.prey))%*%p + rnorm(1,0,sigmaResid)
    #}
    for(j in 1:num.prey) {
       X[i,1] = X[i,1] + p.contrib[j] * (u[j,1] + sources.sd[j,1]*Z[i,1])
       X[i,2] = X[i,2] + p.contrib[j] * (u[j,2] + sources.sd[j,2]*Z[i,2])       
    }
    #rand.sources = matrix(rnorm(num.prey*num.iso,c(u),sigmaW),num.prey,num.iso)
    #sources.dat[,1] = rand.sources[,1]
    #sources.dat[,3] = rand.sources[,2]
    #X[i,] = t(t(rand.sources)%*%p) + rnorm(num.iso,0,sigmaResid)
  }
  # add resid variance down here
  X = X + matrix(rnorm(N*num.iso,0,sigmaResid),N,num.iso)

  # fit the model using BUGS
  chainLength <- 5000;  # post-burn
  burn <- 5000;
  numPrey = num.prey
  numIsotopes=num.iso
  alpha=rep(1,num.prey)
  vars = sources.sd^2
  data = list("u","N","numPrey","numIsotopes","X","alpha","vars");
    parameters <- c("p","sigma")
    model <- 0;
    model <- bugs (data, inits=NULL, parameters, DIC = TRUE, "bugsMix_noResid.txt", n.chains=1, n.burnin = burn, 
       n.thin = 1, n.iter=(burn+chainLength), program="openbugs", debug = FALSE)
    try(attach.bugs (model, overwrite=TRUE), silent = TRUE)
    modes = 0
    for(i in 1:3) {
      d = density(p[,i])
      modes[i] = d$x[which(d$y==max(d$y))]
    }
    summaryResults[ii,1:21] = as.numeric(c(summary(p[,1]),modes[1],summary(p[,2]),modes[2],summary(p[,3]),modes[3]))
#  model = siarmcmcdirichlet(X, sources.dat)
#  # for each parameter calculate the median, posterior mode, and quantiles
#  modes = 0
#  for(i in 1:4) {
#    d = density(model$output[,i])
#    modes[i] = d$x[which(d$y==max(d$y))]
#  }
#  summaryResults[ii,] = c(as.numeric(c(summary(model$output[,1]),modes[1])),as.numeric(c(summary(model$output[,2]),modes[2])),as.numeric(c(summary(model$output[,3]),modes[3])),as.numeric(c(summary(model$output[,4]),modes[4])))

  write.table(summaryResults[1:ii,],"bugs_out.txt",row.names=F,col.names=T)
}