library(runjags)
library(R2jags)
library(gtools)
library(gdata)

  num.iso = 1  # the number of isotopes
  N = 20       # the number of data points to generate (# of consumers)
  # Z is a matrix of ~ Normal(0,1) numbers to create datasets
  Z = matrix(rnorm(N*num.iso,0,1),N, num.iso)
  # if you generate a small sample, you may want to standardize:
  #Z = Z /(sd(Z)); Z = Z - mean(Z)
  p.contrib = c(0.6,0.3,0.1) # this is a vector of diet proportions
  u.mid = 5  # u.mid = arbitrary location of the isotope to use for simulations
  num.prey = length(p.contrib) # number of sources

  # These are all MCMC parameters...same as the ones we used
  mcmc.chainLength <- as.integer(70000)  # burn-in plus post-burn
  mcmc.burn <- as.integer(50000)
  mcmc.thin = 1;
  mcmc.adapt = 5000
  mcmc.chains = 3       # needs to be at least 2 for DIC
  dic.samples = 5000
  alpha=rep(1,num.prey)  # alpha is the hyperparameters for the Dirichlet

  u.steps = seq(0,10,0.1) # this is the vector of increments to explore over the simulations
  u.steps = u.steps[-1] # we'll ignore 0 separation for now
  
  # note: right now, this only stores the quantiles for the 1st proportion...others could be stored too
  mixSIR.results = matrix(0, length(u.steps),3) # stores quantiles for mixSIR model
  bayes.results = matrix(0, length(u.steps),3)  # stores quantiles for bayesian model
    
  for(k in 1:length(u.steps)) {
    u.dif = u.steps[k]
    
    # u is the vector of means for this isotope
    u = matrix(0, num.prey, num.iso)
    u[,1] = c(u.mid-u.dif, u.mid, u.mid+u.dif)
    # sigma2 is the vector of variances...
    sigma2 = matrix(0, num.prey, num.iso)
    sigma2[,1] = c(1.15,1.5,1.25)
    
    # this series of loops creates a new dataset, X.  Note that only 1 Z was drawn above, and we're using the same deviations over and over
    X = matrix(0, N, num.iso)
    for(i in 1:N) {
      for(j in 1:num.prey) {
      	  for(ii in 1:num.iso) {
           X[i,ii] = X[i,ii] + p.contrib[j] * (u[j,ii] + sqrt(sigma2[j,ii])*Z[i,ii])
        }      
      }
    }
    
    # This block runs the mixSIR-type model
    jags.data = list("u","sigma2", "N", "num.prey", "num.iso","X", "alpha")
    jags.params=c("p")
    model.loc = "mixSIR model.txt"  
    jags.mixSIRmodel = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = FALSE)

    # store the median and 95% credible intervals
    mixSIR.results[k,]=quantile(jags.mixSIRmodel$BUGSoutput$sims.matrix[,1],c(0.025,0.5,0.975))   
    
    # this block sets up the parameters for the bayesian code.  Samp.size is the vector of sample sizes collected for each isotope.  I used a value of 3, because it's about the smallest I've found in some published studies.  I have experimented with larger values...10, 20, 100
    samp.size = matrix(3, num.prey, num.iso)
    sigma.v0 = matrix(0, num.prey, num.iso)
    u.loc = u
    u.scale = matrix(0, num.prey, num.iso)
    sigma.scale2 = matrix(0, num.prey, num.iso) # scale param for var
    sigma.v0.n = samp.size + sigma.v0
    sigma2.n = (sigma.v0*sigma.scale2 + (samp.size-1)*sigma2 + ((u.scale*samp.size)*(u - u.loc)^2)/(u.scale + samp.size))/sigma.v0.n 
    
    # This block runs the fully bayesian model
    jags.data = list("u","sigma2", "N", "num.prey", "num.iso","X", "alpha", "samp.size", "sigma.v0.n", "sigma2.n", "u.loc", "u.scale")
    jags.params=c("p")
    model.loc = "bayes model.txt"  
    jags.newModel = jags(jags.data, inits = NULL, parameters.to.save= jags.params, model.file=model.loc, n.chains = mcmc.chains, n.burnin = mcmc.burn, n.thin = mcmc.thin, n.iter = mcmc.chainLength, DIC = FALSE)
    # again, store the median and 95% intervals
    bayes.results[k,]=quantile(jags.newModel$BUGSoutput$sims.matrix[,1],c(0.025,0.5,0.975))   

        
    if(k==50) {
    	 sir.data.for.plot = jags.mixSIRmodel$BUGSoutput$sims.matrix
    	 bayes.data.for.plot = jags.newModel$BUGSoutput$sims.matrix
    }
  }