# This program take as an entrie the predictions of all the model ensembles and give back the predictions of the multi-model for each Sy,Sm,L, the variation on Sy being due to the cross validation.

ens <- "ensemble"
cor <- "corrected" # choose \corrected\" or \"uncorrected\""
method <- "SMC"

m<- 24 # Number of years


Directory_Communication <- "/Volumes/shared_data/paulb/Results/Communication/Verification/"
Directory_Output <- "/Volumes/shared_data/paulb/Results/Results/MultiModel/" #Directory where you want your output to be stocked.

if (cor == "uncorrected") {correct <- cor}
if (cor == "corrected") {correct <- paste(cor, "_", method, sep = "")}

par(mfrow = c(3,4))

for (k in 1:12) {
	for (l in 0:2) {
		WeightedPrediction<-NULL
		for (y in 1:m){ # The year we want to avoid (for cross validation)
			pred <- NULL
			obs <- NULL
			Test<-NULL
		
			# Computation of F (m x (9*4))
			for (XXX in c("WCT", "MAC", "AM", "ABP")) {
			
				# We need to formate the entry vector into a m x 9 matrix
				predic <- t(matrix(as.vector(as.matrix(read.csv(paste(Directory_Communication, XXX, 
					"_", ens, "_", correct, "_MonthlyMean_Start_Month_", k, "_Pred_Month_", k + l, 
					"_SimulatedCases.csv", sep = ""), sep = ",", header = FALSE))),ncol=m))
							
				Test<-c(Test,predic[y,]) # The forecast of the model we are going to weight
				predic <-predic[-y,] # Put asside year y from the computation of the weight
				pred <- cbind(pred, predic)
			
			}
		
			# Computation of A (m x 1)
			obs <- as.vector(as.matrix(read.csv(paste(Directory_Communication, "MonthlyMean_Start_Month_", k, "_Pred_Month_", k + l,"_PositiveCases.csv", sep = ""), sep = ",", header = FALSE)))
			obs <- obs[-y] # Put asside year y
				
			s<-svd(pred) # Compute the SVD of F
			
			if (l==2 & y==24){ # Plot the explained variance vs total variance when year 24 is out
				plot(seq(1,m-1,by=1), s$d^2/sum(s$d^2), type="h", col=c(2,rep(1,m-1)),ylab="Normalized eigenvalues",xlab="Eigenvectors",main=paste("Sm=",k,", L=",l+1,sep="")) # Plot the normalized eigenvalues
			}
			
			#sv<-s$d[s$d^2/sum(s$d^2)>0.01] # Keep only the sv^2>0.01*Total Variance
			sv<-s$d
			p<-length(sv) # Number of dimention kept
			V<-s$v[,1:p] # Extract corresponding matrix V
			U<-s$u # Matrix of the ev
		
			Sigma<-rbind(diag(sv),matrix(data=0,nrow=(m-1-p),ncol=p)) # Simgma_p
			Ts<-U %*% Sigma # Time series
		
			C<-t(Ts) %*% obs # Matrix of correlations
			EV<-diag(sv^2) # EV=t(Ts)Ts
			
			WeightedPrediction[y] <- Test %*% V %*% solve(EV) %*% C
		}
		write.table(WeightedPrediction, file = paste(Directory_Communication, "MultiModel3", "_", "single", "_", correct, "_MonthlyMean_Start_Month_", k, "_Pred_Month_", k + l, "_SimulatedCases.csv", sep = ""), quote = FALSE, sep = ",", row.names = FALSE, col.names = FALSE)
	}
}