#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<-"MAC" # WCT, MAC, AM, ABP or multiModel. Select the model to be tested
ens<-"ensemble" # choose "single" or "ensemble"
cor<-"corrected" # choose "corrected" or "uncorrected"
method<-"SLC" # choose "SMC" or "SLC"

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_A<-c("Scores","BS","BS Baseline","SS","ROC")
Summary_A<-NULL;
par(mfrow = c(3,4))
for (ev in 1:1){
	for (k in 1:12) {
		SummaryMonth_A<-NULL;
		for (l in 0:2) {
			pred_A<- as.vector(as.matrix(read.csv(paste(Directory_Communication,XXX,"_",ens,"_",cor,"_",method,"_mean_Anomalies_Start_Month_",k,"_Pred_Month_",k+l,"_Prob_Event",ev,"_SimulatedCases.csv",sep=""),sep=",",header=FALSE)))
			obs_A<-as.vector(as.matrix(read.csv(paste(Directory_Communication,XXX,"_",ens,"_",cor,"_",method,"_mean_Anomalies_Start_Month_",k,"_Pred_Month_",k+l,"_Bin_Event",ev,"_PositiveCases.csv",sep=""),sep=",",header=FALSE)))
			baseline_A<-as.vector(as.matrix(read.csv(paste(Directory_Communication,XXX,"_",ens,"_",cor,"_",method,"_mean_Anomalies_Start_Month_",k,"_Pred_Month_",k+l,"_Bin_Event",ev,"_BaselineCases.csv",sep=""),sep=",",header=FALSE)))
			
			predic<-prediction(pred_A,obs_A)
			rocScore <- performance(predic, "auc")
			bs <- performance(predic, "rmse")
			baseline <- prediction(baseline_A,obs_A)
			bsBaseline <- performance(baseline, "rmse")
			ss <- 1- bs@y.values[[1]]/bsBaseline@y.values[[1]]
			if (l==2){
				rocPlot <- performance( predic, "tpr", "fpr" )
				plot(rocPlot, main=paste("Start=",k," End=",k+l))
			}
			Skills_A<-c(bs@y.values[[1]],bsBaseline@y.values[[1]],ss,rocScore@y.values[[1]])
			SummaryMonth_A<-cbind(SummaryMonth_A,Skills_A)
		}
		SummaryMonth_A<-rbind(c(k,k+1,k+2),SummaryMonth_A)
		SummaryMonth_A<-cbind(Names_A,SummaryMonth_A)
		Summary_A<-rbind(Summary_A,SummaryMonth_A)
	}
	write.table(Summary_A,file=paste(Directory_Output,XXX,"_",ens,"_",cor,"_",method,"_Event",ev,"_mean_Anomalies_Probabilistic_Skill_Scores_3_Month_Predictions.csv",sep=""),quote=FALSE,sep=",")
}