# R code developed by:
# International Research Institute for Climate and Society
# Columbia University in the City of New York, USA
# Last update: Jully 3, 2013
# --------------------------------------------------------------------------------
# Dependancies : (see each function for further details)
# 		-- Files "Directory_Communication/XXX_ens_correct_MonthlyMean_Anomalies_Start_Month_k_Pred_Month_k+l_Prob_Eventev_SimulatedCases.csv" for all 1=<k=<12 and 0=<l=<2 generated by Scilab function Verification_MonthlyMean.sce
# 		-- Files "MonthlyMean_Anomalies_Start_Month_k_Pred_Month_k+l_Bin_Eventev_PositiveCases.csv" "MonthlyMean_Anomalies_Start_Month_k_Pred_Month_k+l_Bin_Event_ev_BaselineCases.csv" and  for all 1=<k=<12 and 0=<l=<2 generated by Scilab function ActualCases_And_Climatology.sce (and called by funcion Verification_MonthlyMean.sce)
# 		-- Package "ROCR"
# --------------------------------------------------------------------------------
# Results :
#		-- Write the probalistic scores of the selected model in file "Directory_Output/XXX_ens_correct_Eventev_MonthlyMean_Anomalies_Probabilistic_Skill_Scores_3_Month_Predictions.csv"
# 		-- (Optional) Plot the ROC diagrams for some configurations
# --------------------------------------------------------------------------------
# Parameters :
#		-- Directory_Communication the directory where the files generated by Scilab are
# 		-- Directory_Output the directory where the files containing the scores will be written
# 		-- XXX the model
# 		-- ens="ensemble"
# 		-- correct="corrected_method" or "uncorected" indicate the way the results have been corrected (or not)
# --------------------------------------------------------------------------------
# Actions :
#		-- Load the data
# 		-- Compute the scores
# 		-- Save them
# --------------------------------------------------------------------------------

ProbalisticScores <- function(Directory_Communication, Directory_Output, XXX, ens, correct) {
	Names_A <- c("Scores", "BS", "BS Baseline", "SS", "ROC") # Scores that will be computed
	Summary_A <- NULL
	#dev.new()
	#par(mfrow = c(3,4))
for (ev in 1:2) {
		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, "_", correct, "_MonthlyMean_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, "MonthlyMean_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, "MonthlyMean_Anomalies_Start_Month_", 
					k, "_Pred_Month_", k + l, "_Bin_Event", ev, "_BaselineCases.csv", sep = ""), 
					sep = ",", header = FALSE)))

				# Use functions "prediction" and "performance" of package ROCR (see documentation for further details)
				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]]
				
				# Un-comment these lines if you want to plot the ROC diagrams
				#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, "_", correct, 
			"_", "Event", ev, "_MonthlyMean_Anomalies_Probabilistic_Skill_Scores_3_Month_Predictions.csv", 
			sep = ""), quote = FALSE, sep = ",")
	}
}