#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))
baseline<-as.matrix(read.csv(paste(Directory_Communication,XXX,"_Mean_Anomalies_Baseline.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){
		ACor<- verify(obs[i,], predCor[i,], baseline = baseline[i,],  frcst.type = "cont", obs.type = "cont")
		SkillsCor<-c(A$MAE,A$ME,A$MSE,A$MSE.baseline,A$MSE.pers,A$SS.baseline,A$ACC,A$bias)
	}
}
RMSE_COR<-RMSE_COR / 24
ACOR_COR<-ACOR_COR / 24
RMSE_UNCOR<-RMSE_UNCOR / 24
ACOR_UNCOR<-ACOR_UNCOR / 24

#Plot
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.new()
plot.window(c(1,12), c(min(cbind(ACOR_COR,ACOR_UNCOR)),max(cbind(ACOR_COR,ACOR_UNCOR))))
axis(1, at=seq(1,12,by=1))
#axis(2, las=1, at=seq(1,max(cbind(RMSE_COR,RMSE_UNCOR))))
box()
lines(ACOR_COR, type="o", col=2, pch=1)
lines(ACOR_UNCOR, type="o", col=1, pch=1)
title(main=paste("RMSE for seasonial lead time and yearly prevision by a 9 members ensemble of the ",XXX," model"), col.main="black", font.main=1)
title(xlab="Lead Time, from April to Mars")
title(ylab=paste("RMSE")) 
legend("topright", c("Unecorrected","SMC corrected"), cex=0.8, col=c(1,2), pch=c(1,1), horiz=TRUE);