####Function to run the sampler#####
##SPECIFYING PRIORS######
#########################
##if priors = NULL => uses randomly generated priors
##else: priors is a list of following objects:
        ##MuBeta:= prior mean for betas & intercepts; 
        ##SigmaBeta:= prior variance for betas & intercept;
        ##MuAlpha:= prior mean for alpha; 
        ##SigmaAlpha:= prior variance for alpha;
        ##MuZ; VarZ;
        ##PriorA; PriorB 
##TUNING PARAMETERS#######
##########################
##if tune = NULL => uses auto tuning
##else: tune is a list of following objects:
       ##tuneAlpha = 0.9
       ##tuneBeta = array, dim=c(PP,KK)
       ##tuneInt = vec, len = KK
       ##tuneZ =  list( vec(len = nn[x]])) length of list = KK
##############################################################          


#library(MASS)
HLSM=function(X, Y, initialVals = NULL, priors = NULL, tune = NULL,
	tuneIn = TRUE, TT = NULL,dd, niter,intervention)
{
    #X and Y are provided as list. 
    nn = sapply(1:length(X),function(x) nrow(X[[x]]))
    KK = length(X)
    PP = dim(X[[1]])[3]
    XX = unlist(X)
    YY = unlist(Y)
    YY[which(is.na(YY))] = 0
    XX[which(is.na(XX))] = 0

    #Priors

    if(is.null(priors)){
        MuBeta= rep(0,(PP+1)) 
        SigmaBeta = rep(1,(PP+1)) 
        MuAlpha=0 
        SigmaAlpha = 100 
        MuZ = c(0,0)
        VarZ = c(100,100)
        PriorA = 4.5
        PriorB = 17.5
    }else{
	if(class(priors) != 'list')(stop("priors must be of class list, if not NULL"))
	MuBeta = priors$MuBeta
	SigmaBeta = priors$SigmaBeta
	MuAlpha = priors$MuAlpha
	SigmaAlpha = priors$SigmaAlpha
	MuZ = priors$MuZ
	VarZ = priors$VarZ
	PriorA = priors$PriorA
	PriorB = priors$PriorB
  }
##starting values
    if(is.null(initialVals)){
        Z0 = array(0, dim = c(sum(nn),dd))
        cc = 1
        for(i in 1:KK){  
            cc1 = (cc-1)+nn[i]
            ZZ = t(replicate(nn[i],rnorm(dd,0,1)))
            ZZ[1,]=c(0,0)
            ZZ[2,2]=0
            if(ZZ[2,1] < ZZ[1,1]){
                ZZ[2,1] = -1*ZZ[2,1]}
            ZZ[3,2] = abs(ZZ[3,2])
            Z0[cc:cc1,] = ZZ
            cc = cc+nn[i]  
        }
        Z0 = unlist(Z0)
        beta0 = replicate(KK,rnorm(PP,0,1))
        intercept0  = rnorm(KK, 0,1)
        if(intervention == 1){    alpha0=rnorm(1, 0, 1) } 
    print("Starting Values Set")
    }else{
	if(class(initialVals)!= 'list')(stop("initialVals must be of class list, if not NULL"))
	Z0 = initialVals$ZZ
	beta0 = initialVals$beta
	intercept0 = initialVals$intercept
	if(intervention == 1){ alpha0 = initialVals$alpha}
	}

    if(intervention == 0){
        alpha0 = 0
        TT = rep(0, KK)
    }
	
###tuning parameters#####
    if(is.null(tune)){
            a.number = 5
            tuneAlpha = 0.9
            tuneBeta = array(1,dim=c(PP,KK))
            tuneInt = rep(0.2,KK)
            tuneZ =  lapply(1:KK,function(x) rep(1.2,nn[x]))          
            } else{
		if(class(tune) != 'list')(stop("tune must be of class list, if not NULL"))
                 a.number = 1
                 tuneAlpha = tune$tuneAlpha
                 tuneBeta = tune$tuneBeta
                 tuneInt = tune$tuneInt
                 tuneZ = tune$tuneZ
          }       
  
###Tuning the Sampler####
    do.again = 1
    tuneX = 1
    if(tuneIn == TRUE){
    while(do.again ==1){
        print('Tuning the Sampler')
        for(counter in 1:a.number){
            rslt = MCMCfunction(nn=nn,PP=PP,KK=KK,dd=dd,XX = XX,YY = YY,ZZ = Z0,TT = TT,beta = beta0 ,intercept = intercept0,
		alpha = alpha0,MuAlpha = MuAlpha,SigmaAlpha = SigmaAlpha,MuBeta = MuBeta,SigmaBeta = SigmaBeta,MuZ = MuZ,
		VarZ = VarZ,tuneBetaAll = tuneBeta, tuneInt = tuneInt, tuneAlpha = tuneAlpha,tuneZAll = unlist(tuneZ),
		niter = 200,PriorA = PriorA, PriorB = PriorB, intervention = intervention)

        tuneAlpha = adjust.my.tune(tuneAlpha, rslt$acc$alpha,1)
        tuneZ = lapply(1:KK,function(x)adjust.my.tune(tuneZ[[x]], rslt$acc$Z[[x]], 2))
        tuneBeta = array(sapply(1:KK,function(x)adjust.my.tune(tuneBeta[,x], rslt$acc$beta[,x],1)),dim = c(PP,KK))
        tuneInt = sapply(1:KK,function(x)adjust.my.tune(tuneInt[x],rslt$acc$intercept[x], 1))
        print(paste('TuneDone = ',tuneX))
        tuneX = tuneX+1
    }
    extreme = lapply(1:KK,function(x)which.suck(rslt$acc$Z[[x]],2))
    do.again = max(sapply(extreme, length)) > 5
 
}
    print("Tuning is finished")  
}
    rslt = MCMCfunction(nn=nn,PP=PP,KK=KK,dd= dd,XX = XX,YY = YY,ZZ = Z0,TT = TT,beta = beta0 ,intercept = intercept0,
		alpha = alpha0,MuAlpha = MuAlpha,SigmaAlpha = SigmaAlpha,MuBeta = MuBeta,SigmaBeta = SigmaBeta,MuZ = MuZ,
		VarZ = VarZ,tuneBetaAll = tuneBeta, tuneInt = tuneInt, tuneAlpha = tuneAlpha,tuneZAll = unlist(tuneZ),
		niter = niter,PriorA = PriorA, PriorB = PriorB,intervention = intervention)

    rslt$call = match.call()
    class(rslt) = 'HLSM'
    rslt
}