y = array(readtable("C:/Documents/Master/y.txt", header = false))
M = array(readtable("C:/Documents/Master/M.txt", header = false))


nIter = 25000
X = M
Y = y
L = size(M,1)
R = size(M,2)

#Initial values
    delta = 1e-3
    sigma2r = ones(Float64, 1, R).*0.01
    alphaPlus = ones(Float64, R, 1)./R

#Place holders
    Trho = zeros(Int64, nIter, R)
    TalphaPlus = zeros(Float64, nIter, R)
    Tsigma2r = zeros(Float64, nIter, R)
    Tdelta = zeros(Float64, nIter, R)

#Parameters for random walk
    sig = ones(Float64, 1, R).*5e-2
    Tx = ones(Float64, R, 1).*0.3
    rho = zeros(Int64, 1, R)

    for iter in 1:nIter
        tmp = randperm(R)
        r = tmp[1]
        comp_r = setdiff(tmp, r)
        comp_r = comp_r[1]

        for k in comp_r
            ind_k = find(comp_r == k)
            Sk = sigma2r[comp_r]
            #splice!(Sk, ind_k[1])
            Muk = X[:;comp_r]
            Muk = Muk[:,1:size(Muk,2).!= ind_k[1]]
            alphak = alphaPlus[comp_r]
            #splice!(alphak, ind_k[1])

            #if isempty(alphak) &&  isempty(Sk) && isempty(Muk)
                alphak = 0.0
                Sk = 0.0
                Muk = 0.0
            #end

            if Tx[k] .> 0.4 && (iter - 1)%100 == 0
                sig[k] = sig[k].*5
            elseif Tx[k] .< 0.3 && (iter - 1)%100 == 0
                sig[k] = sig[k]./5
            end

            alpha = alphaPlus[k] + sqrt(sig[k]).*rand(Normal())

            if alpha .> 0.0 && alpha .< (1.0 - sum(alphak))
                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

                d = 0.5.*((norm(Y - MualphaStar).^2)./C_alphaStar) - (norm(Y - Mualpha)^2)./C_alpha + (L ./2).*log(C_alphaStar ./C_alpha)

            else
                d[1] = Inf
            end

            if d[1] .< 0.0 || exp(-d[1]) .> rand()
                alphaPlus[k] = alphaStar
                rho[k] = 1
            else
                rho[k] = 0
            end

            alphaPlus_out = zeros(Float64, R, 1)
            alphaPlus_out[comp_r] = alphaPlus[comp_r]
            alphaPlus_out[r] = 1 - sum(alphaPlus[comp_r])
            #alphaPlus_out = alphaPlus_out'

            if iter%100 == 0
                Tx = mean(Trho[(iter - 99):iter;:],2)
            end

            E = (norm(Y - X*alphaPlus_out)^2 + 2 .*sum(alphaPlus_out.^2).*delta)./(2 .*sum(alphaPlus_out.^2))
            sig2inv = rand(Gamma(L./2, 1 ./E))
            sigma2 = 1 ./sig2inv
            sigma2r = kron(sigma2, ones(Float64, 1, R))

            a = R
            b = R./sigma2r[1]
            delta = rand(Gamma(a, 1 ./b))

            TalphaPlus[iter;:] = alphaPlus_out'
            Trho[iter;:] = rho

        end

        Tsigma2r[iter;:] = sigma2r
        Tdelta[iter;:] = delta

        if iter%1000 == 0
            println(iter, " Iterations completed\n")
        end
    end