y = as.matrix(read.table('y.txt', header = F))
M = as.matrix(read.table('M.txt', header = F, sep = ','))

L = 826
R = 2
nIter = 25000
burnIn = 5000

BayesM = function(nIter, burnIn, Y, X, R, L){
  #Initialization of variables
  delta = 1e-3
  sigma2r = array(1, dim = c(1, R))*0.01
  alphaPlus = array(1, dim = c(R, 1))/R
  
  #---------------#
  # Place holders #
  #---------------#
  Trho = array(0, dim = c(nIter, R))
  TalphaPlus = array(0, dim = c(nIter, R))
  Tsigma2r = array(0, dim = c(nIter, R))
  Tdelta = array(0, dim = c(nIter, R))
  
  #-------------------------------------------------------#
  # Sig and TX and rho are parameters for the random walk #
  #-------------------------------------------------------#
  sig = array(1, dim = c(1, R))*5e-2
  Tx = array(1, dim = c(R, 1))*0.3
  rho = array(0, dim = c(1, R))
  
  t1 = proc.time() # Start timing
  
  for (iter in 1:nIter){
    tmp = sample(R)
    r = tmp[1]
    comp.r = setdiff(tmp, r)
    for (k in comp.r){
      # Choosing at random one of the components #
      
      ind.k = which(comp.r == k)
      Sk = sigma2r[comp.r]
      Sk = Sk[-ind.k]
      Muk = as.matrix(X[, comp.r])
      Muk = Muk[,-ind.k]
      alphak = alphaPlus[comp.r]
      alphak = alphak[-ind.k]
      if (length(alphak) == 0 & length(Sk) == 0 & length(Muk) == 0){
        alphak = 0
        Sk = 0
        Muk = 0
      }
      
      if (Tx[k] > 0.4 && (iter - 1)%%100 == 0){
        sig[k] = sig[k]*5
      } else if (Tx[k] < 0.3 && (iter - 1)%%100 == 0){
        sig[k] = sig[k]/5  
      } 
      
      # Proposed value
      alpha = alphaPlus[k] + sqrt(sig[k])*rnorm(1)
      
      # If constraint is fulfilled, proceed #
      if(alpha > 0 & alpha < (1 - sum(alphak))){
        
        #-----------------------------#
        # Random-walk Metropolis step #
        #-----------------------------#
        alphaStar = alpha
        
        Mualpha = X[,1:(R - 1)]*alphaPlus[1:(R - 1)] + X[,R]*(1-sum(alphaPlus[1:(R - 1)]))
        C.alpha = sigma2r[1:(R - 1)]*(alphaPlus[1:(R - 1)]^2) + sigma2r[R]*(1-sum(alphaPlus[1:(R - 1)]))^2
        
        MualphaStar = X[,k]*alphaStar + Muk*alphak + X[,r]*(1 - (alphaStar + sum(alphak)))
        C.alphaStar = sigma2r[k]*(alphaStar^2) + Sk*(alphak^2) + sigma2r[r]*(1 - (alphaStar + sum(alphak)))^2
        
        # Difference between the logarithms of the distributions #
        d = 0.5*((max(svd(Y - MualphaStar)$d)^2/C.alphaStar) - max(svd(Y - Mualpha)$d)^2/C.alpha) + 
          (L/2)*log(C.alphaStar/C.alpha)
        
      } else {
        d = Inf
      }
      
      # Acceptance / reject #
      if (d < 0.0 | exp(-d) > runif(1)){
        alphaPlus[k] = alphaStar
        rho[k] = 1
        
      } else {
        rho[k] = 0
      }
      
      alphaPlus.out = array(0, dim = c(R, 1))
      alphaPlus.out[comp.r] = alphaPlus[comp.r]
      alphaPlus.out[r] = 1 - sum(alphaPlus[comp.r])
      alphaPlus.out = t(alphaPlus.out)
      
      if (iter %% 100 == 0){
        Tx = colMeans(Trho[(iter - 99):iter,])
      }
      
      E = (max(svd(Y - X%*%t(alphaPlus.out))$d)^2 + 2*sum(alphaPlus.out^2)*delta)/(2*sum(alphaPlus.out^2))
      sig2inv = rgamma(1, L/2 + 1, E)
      sigma2r = rep(1/sig2inv, R)
      
      a = R
      b = R/sigma2r[1]
      delta = rgamma(1, a, b)
      
      TalphaPlus[iter, ] = alphaPlus.out
      Trho[iter, ] = t(rho)
    }
    
    Tsigma2r[iter,] = sigma2r  
    Tdelta[iter,] = delta
    
    telliter = 1000
    if(iter %% telliter == 0) cat(paste("Iteration", iter,"completed\n"))
  }
  t2 = proc.time() # Stop the clock
  ttime = t2 - t1 # Total time
  
  return(list(Time = ttime, Alpha = TalphaPlus[-c(1:burnIn),], 
              Sigma = Tsigma2r[-c(1:burnIn),], Delta = Tdelta[-c(1:burnIn),], Rho = Trho))
}

m1 = BayesM(nIter, burnIn, Y = y, X = M, R = 2, L = 826)