#################################################################
#' Calculate DIC for a Markov chian object of the class BayesEnvelope
#' 

Calculate_DIC<-function(MC,Option="Gelman"){
  if(!(Option %in% c("Gelman","Spiegel")))  {Option="Gelman"; print("methods can only be 'Gelman' or 'Spiegel', for all other choices default is Gelman method")} 
  U=MC$U
  X=MC$X
  #X=as.matrix(data[,6:7])
  #U=as.matrix(data[,1:5])
  D=D_vect(MC,X,U)
  if(Option=="Spiegel"){
    theta_bar=calculate_theta_bar(MC)
    DIC=2*mean(D)- calculateD(theta_bar,X,U)
    
  }
  if(Option=="Gelman"){
    DIC=mean(D)+.5*var(D)
    
  }
  return(DIC)
  
}
 ################################################

tr<-function(A){return(sum(diag(A)))}

calculateD<-function(theta,X,U){
  
  n=dim(X)[1]
  r=dim(U)[2]
  
  u=ifelse(is.null(theta[["gamma"]]),  0,  dim(as.matrix(theta[["gamma"]]))[2])
  
  if(u==0){ 
    Gamma0=theta$gamma0
    Omega0=diag(theta$omega0,nrow=(r-u))   
    x=  tr( U%*%(Gamma0)  %*%  solve(Omega0)  %*%  t(U%*% Gamma0) )
    x=  x  +   n*(   sum(log(theta$omega0))   )
  }
  
  if((0<u)*(u<r)){
    Gamma=theta$gamma
    Gamma0=theta$gamma0
    Omega=diag(theta$omega,nrow=u)
    Omega0=diag(theta$omega0,nrow=(r-u))
    Eta=theta$eta
    x=tr(   (U%*%(Gamma)-X%*%t(Eta))  %*%  solve(Omega)  %*%   t(U%*%(Gamma)-X%*%t(Eta))    )    
    x=x+   tr( U%*%(Gamma0)  %*% solve(Omega0)  %*% t(U%*% Gamma0) )
    x=x+n*(   sum(log(theta$omega))   + sum(log(theta$omega0)))
  }
  
  if(u==r){
    Gamma=theta$gamma
    Omega=diag(theta$omega,nrow=u)
    Eta=theta$eta
    x=tr(      (U%*%(Gamma)-X%*%t(Eta))  %*%  solve(Omega)  %*%   t(U%*%(Gamma)-X%*%t(Eta))    )    
    x=x+n*(        sum(log(theta$omega))     ) 
    
  }
  return(x)
}


D_vect<-function(MC,X,U){
  D=0
  
  len=length(MC[[1]])
  u=MC$u
  r=dim(U)[2]
  for( i in 1:len ){
    theta=list()
    if(u<r){  theta$gamma0=((MC[["O_all"]])[[i]])[,(u+1):r]  }
    
    theta$omega=pick_ith_element(MC$omega_all,i)
    theta$omega0=pick_ith_element(MC$omega0_all,i)
    if(u>0){
      theta$eta=MC[["eta"]][[i]]
      theta$gamma=((MC[["O_all"]])[[i]])[,1:u] 
    }
    D[i]=calculateD(theta,X,U)
  }
  
  return(D)
  
}

#################################################################
pick_ith_element<-function(M,i){
  if(is.vector(M)){return(M[i])}
  if(is.matrix(M)){return(M[,i])}
  return(NULL)
  
}


###################################################################

Generate_MC_fromAllList<-function(allList, burnIn=1000){
  if(burnIn<1){MC=allList}
  if(burnIn>=1){
  len=length(allList[[1]])  
  #if(burnIn<1){burnIn=1}
  j=burnIn
  burnIn=burnIn+1
  MC=list()
  MC$O_all=allList$O_all[burnIn:len]
  MC$omega_all=allList$omega_all[,burnIn:len]
  MC$omega0_all=allList$omega0_all[,burnIn:len]
  MC$u=allList$u
  MC$eta=allList$eta[(burnIn-1):(len-1)]
  
  MC$beta=allList$beta[(burnIn-1):(len-1)]
  MC$Prior=allList$Prior
  MC$X=allList$X
  MC$U=allList$U
  MC$iterationLength=allList$McLength-burnIn+1
  MC$burnIn=burnIn-1
  }
  return(MC)
  
  
}


#######################################################################
calculate_theta_bar<-function(MC){
  u=MC$u
  r=dim((MC$O_all)[[1]])[1]
  theta=list()
  
  if(u==0){
    O_mean=Mean_list(MC$O_all)
    theta$gamma=NULL
    theta$gamma0=(O_mean)[,(u+1):r]
    OmegaAVG=apply(rbind(MC$omega_all,MC$omega0_all),1,mean)
    theta$omega=NULL
    theta$omega0=OmegaAVG[(u+1):r]
    theta$eta=Mean_list(MC$eta)
  }
  
  if((u>0)*(u<r)){
    O_mean=Mean_list(MC$O_all)
    theta$gamma=(O_mean)[,1:u]
    theta$gamma0=(O_mean)[,(u+1):r]
    OmegaAVG=apply(rbind(MC$omega_all,MC$omega0_all),1,mean)
    theta$omega=OmegaAVG[1:u]
    theta$omega0=OmegaAVG[(u+1):r]
    theta$eta=Mean_list(MC$eta)
  }
  if(u==r){
    O_mean=Mean_list(MC$O_all)
    theta$gamma=(O_mean)[,1:u]
    theta$gamma0=NULL
    OmegaAVG=apply(rbind(MC$omega_all,MC$omega0_all),1,mean)
    theta$omega=OmegaAVG[1:u]
    theta$omega0=NULL
    theta$eta=Mean_list(MC$eta)
  }
  
  
  
  
  return(theta)  
}

########################################################################

Mean_list<-function(lst){
  n=length(lst)
  sum=lst[[1]]
  for(i in 2:n){
    sum=sum+lst[[i]]
    
  }
  return(sum/n)
  
  
}



Var_list<-function(lst){
  n=length(lst)
  all=as.vector(lst[[1]])
  for(i in 2:n){
    all=rbind(all,as.vector(lst[[i]]))
    
  }
  return(apply(all,2,var))
  
  
}