#library(verification,lib.loc="/home/paulb/R/x86_64-redhat-linux-gnu-library/2.15");
source("/Volumes/shared_data/paulb/R/Functions/verify_mod.R") #I modified the function "verify" of the package to add the computation of ACC

#Parameters of the program you can change
########################################################################
XXX="WCT" # WCT, MAC, AM, ABP or multiModel. Select the model to be tested
Directory_Data<-"/Volumes/shared_data/paulb/Data/Kenya/" #Directory where the parameters are defined
Directory_Communication<-"/Volumes/shared_data/paulb/Results/Communication/Verification/"
Directory_Output<-"/Volumes/shared_data/paulb/Results/Results/Skills/" #Directory where you want your output to be stocked. 

###################
# Lauch on scilab :
# exec('~/Scilab/infectedVerification/infected_XXX') and
# exec('~/Scilab/infectedVerification/Climatology.sci')
###################
# Computation of the scores for the monthly mean values
Names<-c("Scores","MAE","ME","MSE","MSE baseline","MSE persistence","SS baseline","ACC","bias")

RMSE_COR<-mat.or.vec(12,1)
ACOR_COR<-mat.or.vec(12,1)
RMSE_UNCOR<-mat.or.vec(12,1)
ACOR_UNCOR<-mat.or.vec(12,1)

obs<-as.matrix(read.csv(paste(Directory_Communication,XXX,"_Mean_Anomalies_SimulatedCases.csv",sep=""),sep=",",header=FALSE))

for (k in 1:9) { # Each member of the ensemble
	
	predCor<-as.matrix(read.csv(paste(Directory_Communication,XXX,"_Corrected_Mean_Anomalies_Answer_",k,"_SimulatedCases.csv",sep=""),sep=",",header=FALSE))
	predUncor<-as.matrix(read.csv(paste(Directory_Communication,XXX,"_Mean_Anomalies_Answer_",k,"_SimulatedCases.csv",sep=""),sep=",",header=FALSE))
	
	for (i in 1:12){
		RMSE_COR[i]<-RMSE_COR[i]+mean((predCor[i,]-obs[i,])^2)
		ACOR_COR[i]<-ACOR_COR[i]+(sum(predCor[i,]*obs[i,]))/(sqrt(sum(predCor[i,]^2))*sqrt(sum(obs[i,]^2)))
		RMSE_UNCOR[i]<-RMSE_UNCOR[i]+mean((predUncor[i,]-obs[i,])^2)
		ACOR_UNCOR[i]<-ACOR_UNCOR[i]+(sum(predUncor[i,]*obs[i,]))/(sqrt(sum(predUncor[i,]^2))*sqrt(sum(obs[i,]^2)))
	}
}
RMSE_COR<-RMSE_COR / 24
ACOR_COR<-ACOR_COR / 24
RMSE_UNCOR<-RMSE_UNCOR / 24
ACOR_UNCOR<-ACOR_UNCOR / 24

#Plot
#Plot
source("/Volumes/shared_data/paulb/R/Functions/Plot_Comparison.R")
par(mfrow = c(2,1))

Plot_Comparison(RMSE_UNCOR,RMSE_COR,"Unecorrected","SMC corrected","Lead Time, from April to Mars","RMSE",paste("RMSE for seasonial lead time and yearly prevision by a 9 members ensemble of the ",XXX," model"))

Plot_Comparison(ACOR_UNCOR,ACOR_COR,"Unecorrected","SMC corrected","Lead Time, from April to Mars","ACOR",paste("ACOR for seasonial lead time and yearly prevision by a 9 members ensemble of the ",XXX," model"))