library("BRugs")
library("R2WinBUGS")

###################################
# STEP 1: SET PARAMETERS
###################################
numClones<-50;  # must be at least 2
chainLength <- 2000;  # min 2000, any statistic = function of sample size
burn <- 4000; # min 4000

###################################
# STEP 2: READ IN THE DATA
###################################
d <- read.csv("dat.csv",header=TRUE)

# replace 0s with NAs
for(j in 1:dim(d)[2]) {
   ind <- which(d[,j] < 0);
   d[ind,j] <- NA;
   d[-ind,j] <- log(d[-ind,j])
}
N <- dim(d)[2]; # number of sites
numY <- dim(d)[1];  # number of years
DATA <- matrix(0,numY,N);
for(i in 1:N) DATA[,i] <- d[,i];

#####################################################################
# USER NEEDS TO CHANGE THE FLAGS (F2, F3) AND THE GROUP VECTORS BELOW
# To summarize, there are N sites being monitored, and there may be M subpopulations (1 <= M <= N).
# The vector 'groups' assigns sites to subpopulations, and the remaining vectors allow for parameters
# to be shared among subpopulations.  For example, if 2 subpopulations have the same process error vars,
# but different growth rates and observation errors, U.groups = c(1,2), Q.groups = c(1,1), R.groups = c(1,2)
#####################################################################
f1 <- 4; # FLAG FOR B (1, 2, 3, 4), interaction matrix
#case 1 % no interactions (identity matrix)
#case 2 % diagonal matrix, equal
#case 3 % diagonal matrix, unequal
#case 4 % unconstrained interaction matrix
f2 <- 1; # FLAG FOR Q (1, 2, 3), process error
f3 <- 2; # FLAG FOR R (1, 2, 3), observation error
#case 1 %unconstrained    
#case 2 %diagonal covariance matrix, variances can be equal or not
#case 3 %equal variances and covariances; the same but not independent

# THESE NEXT 4 VECTORS CAN BE MODIFIED FOR DIFFERENT STRUCTURES
groups <- c(1,2,3,4);  # this is an Nx1 vector describing which site belongs to which group
M <- max(groups);      # M is the clusters / groups / signals, DO NOT CHANGE THIS LINE
if(M==1) f2 = 2;       # this line prevents an error from being thrown
U.groups <- rep(1,M);  # this is an Mx1 vector describing which group-specific growth rates are unique
#Q.groups <- seq(1,M);
Q.groups <- rep(1,M);# this is an Mx1 vector describing which group-specific process errors are unique
#R.groups <- seq(1,N);
R.groups <- rep(1,N);  # this is an Mx1 vector describing which group-specific observation errors are unique

############################################################
#  USER DOES NOT CHANGE THESE VALUES                                                          
############################################################
U.G <- max(U.groups);  # this is the number of unique Us being estimated
Q.G <- max(Q.groups);  # this is the number of unique Qs being estimated
R.G <- max(R.groups);  # this is the number of unique Rs being estimated
source("writeModel.R");
# note: B can't be in the parameter list if it is fixed
if(f1 == 1) parameters <- c("U","Q","R","A","states1")
if(f1 > 1) parameters <- c("B","U","Q","R","A","states1")
mod <- bugs(data,inits=NULL,parameters,DIC=TRUE,"model.txt",n.chains=1,n.burnin=burn, 
   n.thin=1,n.iter=(burn+chainLength),program="openbugs",debug=TRUE)
attach.bugs(mod, overwrite=TRUE)
dev <- deviance/numClones;
nPar <- mean(dev)-min(dev)

#bestDev <- which(dev==min(dev))

# calculate the deviance at the mean posterior
pD <- mod$pD;
Dbar <- mean(deviance);
Dhat <- (Dbar - pD)/numClones;
Dhat <- (mod$DIC - 2*pD)/numClones;

z<- rep(0,N);
z[1] <- min(dev);
z[2] <- mean(dev);
z[3] <- nPar;
z[4] <- pD;
z[5] <- Dbar;
z[6] <- Dhat;
z[7] <- mod$DIC;

# calculate the distribution of the states
mat <- matrix(0,dim(d)[1],dim(d)[2]);
for(i in 1:dim(d)[1]) {
   for(j in 1:dim(d)[2]) {
      mat[i,j]<-(states1[1,i,j]);
   }
}
if(f2 == 1) { 
   # unconstrained
   Qmat <- matrix(0,dim(d)[2],dim(d)[2]);
   if(length(Q) == chainLength) {
      den <- density(Q);
	Qmat[1,1] <- den$x[which(den$y==max(den$y))];
   }
   if(length(Q) > chainLength) {
      for(i in 1:U.G) {
         for(j in 1:U.G) {
      	den <- density(Q[,i,j]);
		Qmat[i,j] <- den$x[which(den$y==max(den$y))];
         }
      }
   }
}
if(f2 == 2) { 
   # constrained
   Qmat <- matrix(0,dim(d)[2],dim(d)[2]);
   if(length(Q) == chainLength) {
      den <- density(Q);
	Qmat[1,1] <- den$x[which(den$y==max(den$y))];
   }
   if(length(Q) > chainLength) {
      for(i in 1:U.G) {
         den <- density(Q[,i]);
	   Qmat[i,i] <- den$x[which(den$y==max(den$y))];
      }
   }
}
if(f3 == 1) { 
   # unconstrained
   Rmat <- matrix(0,dim(d)[2],dim(d)[2]);
   for(i in 1:dim(d)[2]) {
      for(j in 1:dim(d)[2]) {
         den <- density(R[,i,j]);
	   Rmat[i,j] <- den$x[which(den$y==max(den$y))];
      }
   }
}
if(f3 == 2) { 
   # constrained
   Rmat <- matrix(0,dim(d)[2],dim(d)[2]);
   for(i in 1:dim(d)[2]) {
	den <- density(R[,i]);
	Rmat[i,i] <- den$x[which(den$y==max(den$y))];
   }
}