#' Main function to generate the Markov chain 
#' @param U: Matrix containing  the responses
#' @param X: Matrix containing  the predictors
#' @param u: Integer input. u is the dimension  of significant portion of the data 
#' @param McLength: Number of iteration of the markov chian.
#' @param Prior: two option can be used "EB" stands for emperical Bayes prior and "Improper" for the improper flat prior 
#' @param name: Name of the list variable containing all the variables of the Markov chain.
#' @export


BayesEnvelopeMC<-function(U,X,u,McLength=3500,name="DefaultName",Prior="Improper",allstep="allpossiblepair",Max_rejection=500,partition=10000){
  #oneRun_sim2_wh<-function(name,BurnIn=1100,u,Prior="Improper"){
  
  library(rstiefel)
  U=center_data_matrix(as.matrix(U))
  X=center_data_matrix(as.matrix(X))
  
  
  P_X=X%*%(solve(t(X)%*%X))%*%t(X)
  n=dim(P_X)[1]
  I=diag(dim(P_X)[1])  
  p=dim(X)[2]
  ############### Initialization ################
  true_start=F
  #print(lst    *********)
  if(!true_start){  lst=startingValue_genU(U,X,u)}
  if(true_start){  
          load("C:/Users/subhadippal/Desktop/LpSn/data/Data_RV_smallNlargeP_1.RData")
          lst=startingValue_genU1(lst)
        }
  
  gamma=lst[["gamma"]]
  gamma0=lst[["gamma0"]]
  omega=lst[["omega"]]
  omega0=lst[["omega0"]]
  O=cbind(gamma,gamma0)
  #u=dim(gamma)[2]
  r=dim(O)[2]
  
  eta=beta=list()
  
  omega_all=omega
  omega0_all=omega0
  O_all=list(O)
  
  ####################### Prior Choices ####################
  if(Prior=="EB"){
    
    e_prior=lst$beta_start
    C_prior=0*diag(p) 
    if(r>n){C_prior=diag(p)}
    #G_prior= matrix(0,nrow=r,ncol=r)
    G_prior=G_prior_EB_selection(O,omega,omega0)
    D_prior=c(omega,omega0)
    if(u==0){
      e_prior=matrix(0,nrow=r,ncol=p)
      lambda0_prior=mean(1/omega0)/(  mean(1/omega0^2)-(mean(1/omega0))^2 )
      alpha0_DF_prior=mean(1/omega0)*lambda0_prior
    }
    
    if(u==1){
      alpha_DF_prior=1/omega
      lambda_prior= 1
      lambda0_prior=mean(1/omega0)/(  mean(1/omega0^2)-(mean(1/omega0))^2 )
      alpha0_DF_prior=mean(1/omega0)*lambda0_prior
    }
    if(u==(r-1)){
      lambda_prior=mean(1/omega)/ (  mean(1/omega^2)-(mean(1/omega))^2 )
      alpha_DF_prior=mean(1/omega)*lambda_prior
      alpha0_DF_prior=1/omega0 
      lambda0_prior= 1
      
    }
    if(   (u<(r-1))*(u>1)   ){
      lambda_prior=mean(1/omega)/ (  mean(1/omega^2)-(mean(1/omega))^2 )
      lambda0_prior=mean(1/omega0)/(  mean(1/omega0^2)-(mean(1/omega0))^2 )
      alpha_DF_prior=mean(1/omega)*lambda_prior
      alpha0_DF_prior=mean(1/omega0)*lambda0_prior
      print(c(alpha_DF_prior,alpha0_DF_prior))
    }
    
    if(u==r){
      lambda_prior=mean(11/omega)/ (  mean(1/omega^2)-(mean(1/omega))^2 )
      alpha_DF_prior=mean(1/omega)*lambda_prior
    }
    
  }
  
  if(Prior=="Improper"){
    e_prior=matrix(0,nrow=r,ncol=p)
    C_prior=0*diag(p) 
    G_prior= matrix(0,nrow=r,ncol=r)
    D_prior=c(omega,omega0) # actually value of D_prior doesnot matter as G=0 and all we need is G/D[i]
    lambda_prior=0
    lambda0_prior=0
    alpha_DF_prior=-1-p/2
    alpha0_DF_prior=-1
    
  }
  
  
  
  ##########################################################
  
  #U=U-t(replicate(n=dim(U)[1],expr=(apply(U,2,"mean") ) ))
  e_tilde=(  t(U) %*% X+ e_prior%*% C_prior  ) %*% solve(t(X)%*% X+C_prior)
  G_tilde=  t(U)%*% U + e_prior%*% C_prior %*% t(e_prior) -  e_tilde %*%  (t(X)%*% X+C_prior) %*% t(e_tilde)
  
  i=1
  while(i<=McLength){
    
    if(u==0){
      #gamma=as.matrix(O[,1:u])
      gamma0=as.matrix(O[,(u+1):r])
      #A=t(gamma)%*% t(U)%*%(I-P_X)%*%U %*%gamma
      #B=t(gamma0)%*%t(U)%*%U%*%gamma0
      B=t(gamma0)%*%( t(U)%*%U)%*%gamma0 + 2*lambda0_prior
      
      omega=NULL
      omega0=UpdateOmega(omega0, diag(B)/2, alpha=(n+2*alpha0_DF_prior-1)/2,Max_rejection=Max_rejection  )
      
    }
    
    if(  (0<u)*(u<r)  ){
      gamma=as.matrix(O[,1:u])
      gamma0=as.matrix(O[,(u+1):r])
      #A=t(gamma)%*% t(U)%*%(I-P_X)%*%U %*%gamma
      A=t(gamma)%*% G_tilde  %*%gamma + 2*lambda_prior
      B=t(gamma0)%*%( t(U)%*%U  )%*%gamma0 + 2*lambda0_prior
      
      omega =UpdateOmega(omega , diag(A)/2, alpha=(n+2*alpha_DF_prior-1)/2 ,Max_rejection=Max_rejection ) ### pass the diagonal 
      omega0=UpdateOmega(omega0, diag(B)/2, alpha=(n+2*alpha0_DF_prior-1)/2,Max_rejection=Max_rejection )
    }
    
    
    if(u==r){
      gamma=as.matrix(O[,1:u])
      #gamma0=as.matrix(O[,(u+1):r])
      ##A=t(gamma)%*% t(U)%*%(I-P_X)%*%U %*%gamma
      A=t(gamma)%*% G_tilde  %*%gamma + 2*lambda_prior
      #B=t(gamma0)%*%t(U)%*%U%*%gamma0
      
      omega=UpdateOmega(omega,diag(A)/2,alpha=(n+2*alpha_DF_prior-1)/2,Max_rejection=Max_rejection ) ### pass the diagonal 
      omega0=NULL
      
    }
    
    
    O=UpdateO(O=O,G_tilde=G_tilde,G_prior=G_prior, D_prior=D_prior, U = U, omega=omega, omega0=omega0, allstep=allstep, discrete=F,MAX_rejection=Max_rejection,partition=partition)
    
    
    
    
    
    
    
    if(u>0){
      eta[[i]]=SampleEta(X, U, omega, O[,1:u],C_prior=C_prior,e_prior=e_prior)
      beta[[i]]=as.matrix( O[ , 1:u] )  %*%  eta[[i]]
    }
    
    O_all[[i]]=O  
    ##################  storing of different variables ################
    if(i==1){
      omega_all=omega  
      omega0_all=omega0_all
    }
    
    if(i>1){
      omega_all=cbind(omega_all,omega)	
      omega0_all=cbind(omega0_all,omega0)
    }
    
    
    print(i) 
    i=i+1
    
  }
  
  allList=list()
  allList$O_all=O_all
  allList$omega_all=omega_all
  allList$omega0_all=omega0_all
  allList$u=u
  allList$eta=eta
  allList$beta=beta
  allList$Prior=Prior
  allList$X=X
  allList$U=U
  allList$iterationLength=McLength
  allList$burnIn=0
  save(allList,file=paste ("output/run_",Prior, name,"u_equals_",u,".RData",sep="")) 
  return(allList)
}