##################################################
#### Script to configure all needed variances from user inputs


source('Config//SetVariances.R')
## Bayes flags
eBayes.flag.SORO <- eBayes.flag.RFSP <- eBayes.flag.TUF <-FALSE


if(!variance.set){
	source('Code//01-01-06-Var-SORO-est.R')
}else{
    if(Ssoro == 1){
 		print('SORO: Exchangeable covariance model selected')
		SORO.vs<-set.SORO.vs
		SORO.vc<-set.SORO.vc
		Vsoro.fit<-diag(SORO.vs,Ni) + matrix(SORO.vc,Ni,Ni)
	}else if(Ssoro == 3){
		print('SORO: Spatial covariance model selected')
		SORO.v<-set.SORO.v
		SORO.l<-set.SORO.l
		Vsoro.fit<-SORO.v*exp(-(D11/SORO.l)^2)
	}else if(Ssoro == 5){
        ## Loading SORO from CSV files
        ## Matrix should be Ni x Ni = 480 x 480, Positive definite and symmetric
        if(!file.exists("Data//Vsoro.csv"))
            stop('Switch Ssoro=5, but file Data//Vsoro.csv does not exists!')
        Vsoro<-data.matrix(read.csv("Data//Vsoro.csv",header=FALSE))
		Vsoro.fit<-Vsoro
	}
    if(Sfinch){
	    Vsoro.fit<-matrix(0,Ni,Ni)
        diag(Vsoro.fit)<-Vsoro.d
    }
}

#Sanity check: is Vsoro.fit too big ??
## Re introduced sanity checks per request from Marcus, email 5th April, 2019
if(prior.switch.SORO==0 & prior.switch.SORO.Tony){
    obs.rip.var<-apply(r.s.b,2,var)
    Wsd<-obs.rip.var-diag(Vsoro.fit)
    if(min(Wsd) <= 0){
        print('The specified Vsoro values are too large: hard failure')
        print(paste('The average ripple variance across channels is',mean(obs.rip.var)))
        print('Please reset Vsoro variances to smaller values if you want to use Empirical Bayes')
        eBayes.flag.SORO <- TRUE
    }else{
        print('The specified Vsoro values are small enough')
    }
}

Vsoro<-Vsoro.fit
Vsoro.inv<-solve(Vsoro)
lower.Vs<-t(chol(Vsoro))


####################################################################
#RFSP

if(!variance.set){
	source('Code//01-01-07-Var-RFSP-est.R')
}else{
	if(Srfsp == 2){
 		print('RFSP: Exchangeable covariance model selected')
		RFSP.vs<-set.RFSP.vs2
		RFSP.vc<-set.RFSP.vc2
		RFSP.vo<-set.RFSP.vo2
		RFSP.ws<-set.RFSP.ws2
		RFSP.wc<-set.RFSP.wc2
		Vrfsp.fit<-cbind(diag(RFSP.vs,Ni) + matrix(RFSP.vc+RFSP.vo,Ni,Ni),matrix(RFSP.vo,nrow=Ni,ncol=Nd))
		Vrfsp.fit<-rbind(Vrfsp.fit,cbind(matrix(RFSP.vo,nrow=Nd,ncol=Ni),diag(RFSP.ws,Nd) + matrix(RFSP.wc+RFSP.vo,Nd,Nd)))
	}else if(Srfsp == 4){
		print('RFSP: Spatial covariance model selected')
		RFSP.vc<-set.RFSP.vc4
		RFSP.vo<-set.RFSP.vo4
		RFSP.wc<-set.RFSP.wc4
		RFSP.l<-set.RFSP.l4
		Vrfsp.fit<-cbind((RFSP.vc+RFSP.vo)*exp(-(D11/RFSP.l)^2),RFSP.vo*exp(-(D12/RFSP.l)^2))
		Vrfsp.fit<-rbind(Vrfsp.fit,cbind(RFSP.vo*exp(-(D21/RFSP.l)^2),(RFSP.wc+RFSP.vo)*exp(-(D22/RFSP.l)^2)))
	}else if(Ssoro == 5){
        ## Loading Vrfsp from CSV filr
        ## Matrix should be No x No = 534 x 534, Positive definite and symmetric
        if(!file.exists("Data//Vrfsp.csv"))
            stop('Switch Ssoro=5, but file Data//Vrfsp.csv does not exists!')
        Vrfsp<-data.matrix(read.csv("Data//Vrfsp.csv",header=FALSE))
		Vrfsp.fit<-Vrfsp
	}
    if(Sbrfsp){
        Vrfsp.fit <- 0*Vrfsp.fit
        diag(Vrfsp.fit)<-Vrfsp.d
    }
}

#Sanity check: is Vrfsp.fit too big ??
## re-introduced sanity check by request from Marcus, email 5 April, 2019
if(prior.switch.TUF==0 & prior.switch.TUF.Tony){
    obs.op.var <- c(apply(cop.J.full,2,var),apply(fop.J.full,2,var))
    Rsd <- obs.op.var-diag(Vrfsp.fit)
    if(min(Rsd) <= 0){
        print('The specified Vrfsp values are too large')
        print(paste('The average overpower variance is',mean(obs.op.var)))
        print('Please reset Vrfsp variances to smaller values if you want to use Empirical Bayes')
        eBayes.flag.RFSP <- TRUE
    }else{
        print('The specified Vrfsp values are small enough')
    }
}
Vrfsp<-Vrfsp.fit
lower.Vr<-t(chol(Vrfsp))
Vrfsp.inv<-solve(Vrfsp)



# Tuf corresponds to Critical Channel Power Error
# (see "2010_Discussion of some issues related to errors.pdf")

if(!variance.set){
	source('Code//01-01-08-Var-Tuf-est.R')
}else{
    if(Stuf == 1){
        print('TUF: exchangeable covariance model selected')
	    Vtuf = matrix(set.Vtuf.grand,Ni,Ni)                                    #Grand Mean
        Vtuf[South,South]<-Vtuf[South,South] + set.Tuf.vc1.s                   #South Common
        Vtuf[North,North]<-Vtuf[North,North] + set.Tuf.vc1.n                   #North Common
        diag(Vtuf)[South]<-diag(Vtuf)[South] + set.Tuf.vs1.s                   #South specific
        diag(Vtuf)[North]<-diag(Vtuf)[North] + set.Tuf.vs1.n                   #North specific
		Vtuf.fit<-Vtuf
	}else if(Stuf == 3){
        print('TUF: spatial covariance model selected')
        Vtuf = exp(-abs(outer(1:Ni,1:Ni,'-'))/set.Tuf.l3)                      #Spatial Correlation
        Tv = diag(c(rep(sqrt(set.Tuf.v3.s),Ni/2),                              #St dev for South and North
                    rep(sqrt(set.Tuf.v3.n),Ni/2)))
        Vtuf = Tv %*% (Vtuf %*% Tv)                                            #Variance
		Vtuf.fit<-Vtuf
	}else if(Stuf == 5){
        ## Loading Vtuf from CSV file
        ## Matrix should be Ni x Ni = 480 x 480, Positive definite and symmetric
        if(!file.exists("Data//Vtuf.csv"))
            stop('Switch Stuf=5, but file Data//Vtuf.csv does not exists!')
        Vtuf<-data.matrix(read.csv("Data//Vtuf.csv",header=FALSE))
		Vtuf.fit<-Vtuf
	}
}

## Sanity check: is Vtuf.fit too big ??
## re-introduced sanity checks as request by Marcus, email 5 April, 2019
if(prior.switch.TUF==0 & prior.switch.TUF.Tony){
    obs.ccp.var<-apply(ccp.J.full,2,var)
    Csd<-obs.ccp.var-diag(Vtuf.fit)
    if(min(Csd) <= 0){
        print('The specified Vtuf values are too large compared to the CCP data')
        print(paste('The average critical channel power variance is',mean(obs.ccp.var)))
        print('Please reset Vtuf variances to smaller values if you want to use Empirical Bayes')
        eBayes.flag.TUF <- TRUE
    }else{
        print('The specified Vtuf values are small enough')
    }
}
Vtuf<-Vtuf.fit
UVtuf = chol(Vtuf) #chol gives the U in LU
LVtuf = t(UVtuf)   #This is the lower triangle
Vtuf.inv<-solve(Vtuf)


#####################################################################
#These codes may only be needed if the superposition
#model is not sufficient

source('Code//01-01-09-Var-Zsop.R')

if(Mop == 1){
	#For no superposition, need the extra variance for perturbation
	eps<-0.01
	epsst<-0.01
	epsstst<-0.01
}

source('Code//01-01-10-Var-Zccp.R')