# R code developed by:
# International Research Institute for Climate and Society
# Columbia University in the City of New York, USA
# Last update: June 26, 2013
# --------------------------------------------------------------------------------
# Dependancies : (see each function for further details)
#		-- function "bootstat.R" from package "sensitivity"
#		-- modified function "src.R" issued from package "sensitivity"
#		-- R function "Random_Param" to randomly generate sets of parameters
#		-- files "XXX_Parameters.csv" (format described below)
#		-- R functions "plot_SRC_seasons.R" and "nodeplot.R" to plot the results
#		-- Scilab program "Compute_Sol.sce"
# --------------------------------------------------------------------------------
# Functions :
#		-- Read the parameters and their ranges in files "XXX_Parameters.csv"
#		-- Generate two random matrix of parameters, each parameter following a certain distribution (see option of funtion "Random_Param")
#		-- Generate the results to be tested by the Scilab codes with the function "sobol2002" (of the library "sensitivity")
#		-- Analyse these results with function "sobol2002" and plot the Sobol indices
# --------------------------------------------------------------------------------

#Format of files XXX_Parameters.csv :
#X            min max base unit
#parameter 1   x   x   x    u
#...
#parameter p   x   x   x    u


#Parameters of the program you can change
########################################################################
XXX="ABP" #WCT, MAC, AM or ABP. Select the model to be tested
Dist="Norm99" #Norm95 or Norm99 or Unif. Select the distribution you want for the parameters. Norm99 is the normal distribution with less than 1% chance of being out of the defined interval.
Directory_Entrie<-"/data/paulb/Data/Parameters/Entries_Sensitivity/" #Directory where the parameters are defined
Directory_Communication<-"/data/paulb/Results/Communication/" #Directory where you want your output to be stocked. Has to match with the entry  directory of Scilab.
Directory_Results<-"/data/paulb/Results/Results/" #Directory where you want your results to be stocked.
Directory_R_Data<-"/data/paulb/R/Data/" #Directory where the R objects will be stocked
Directory_R_Functions <- "/data/paulb/R/Functions"

n <- 1000 # total cost of (p + 2) * n model evaluations
nb <- 10000 #nboot
########################################################################

#Reading the Parameters
Pa<-read.csv(paste(Directory_Entrie,XXX,"_Parameters.csv",sep=""),sep=";"); 
param_min<-Pa$min # Definition of the parameters
param_max<-Pa$max
p<-nrow(Pa);

source(paste(Directory_R_Functions,"Sensitivity/Random_Param.R",sep="")) #see comments in file "Random_Param.r"
X<-data.frame(Random_Param(n,param_min,param_max,Dist))
write.table(X,file=paste(Directory_Communication,XXX,"_SRC_ParamToBeTested_",n,"_",nb,".csv",sep=""),quote=FALSE,sep=",",row.name=FALSE,col.name=FALSE)
save.image(file=paste(Directory_R_Data,XXX,"_SRC_ParamToBeTested_",n,"_",nb,"_UNFINISHED.RData",sep=""))

##################################################################
##################################################################
##Execute Scilab Function Compute_Sol.sce ##
##################################################################
##################################################################

#If output at several times, execute these lines

########################################################################
source(paste(Directory_R_Functions,"Sensitivity/src.R",sep=""))
source(paste(Directory_R_Functions,"Sensitivity/bootstats.R",sep="")) # Only these two function of the package sensitivity are used. src.R has been slightly modified  (see comments in the function) 
library(boot)
Seasons<-c("winter","spring","summer","Autumn");
Y <-as.matrix(read.csv(paste(Directory_Communication,XXX,"_SRC_Climatology_ComputedResults_",n,"_",nb,".csv",sep=""),header = FALSE,sep = ",")); # Read the answers given by Scilab. Has to match with the output of function Compute_Sol.sce
Xx<-list()
for (k in 1:4) {
	x <- src(X, Y[k,], nboot = nb) # compute the SRC indices
	rownames(x$SRC)<-Pa$X; # Give their names to the parameters.
	print(x);
	write.table(x$SRC,file=paste(Directory_Results,XXX,"_SRC_Indices_",n,"_",nb,"_",Seasons[k],".csv",sep=""),quote=FALSE,sep=",")#Save the SRC
	Xx[[k]]<-list()
	Xx[[k]]<-x
}
save.image(file=paste(Directory_Results,XXX,"_SRC_Indices_Seasons_",n,"_",nb,".RData",sep=""))
source(paste(Directory_R_Functions,"Sensitivity/plot_SRC_seasons.R",sep="")) # See comment on this function
source(paste(Directory_R_Functions,"Sensitivity/nodeplot.R",sep="")) # Unmodified function from package "sensitivity"
plot_SRC_seasons(Xx)
	
#Plot locally
load("/Volumes/shared_data/paulb/Results/Results/Sensitivity/ABP_SRC_Indices_Seasons_1000_10000.RData")
source("/Volumes/shared_data/paulb/R/Functions/Sensitivity/plot_SRC_seasons.R")
source("/Volumes/shared_data/paulb/R/Functions/Sensitivity/nodeplot.R")
plot_SRC_seasons(Xx)
########################################################################

#If output at all times, execute these lines

########################################################################
t=9000-2
Y <-as.matrix(read.csv(paste(Directory_Communication,XXX,"_SRC_AllTimes_ComputedResults_",n,"_",nb,".csv",sep=""),header = FALSE,sep = ","));
All_SRC<-NULL;
#WCT:
names(X)<-Pa$X
for (k in 1:t) {
	x <- src(X, Y[k,], nboot = 100)#compute the SRC indices
	All_SRC<-rbind(All_SRC,x$SRC);
}
write.table(All_SRC,file=paste(Directory_Results,XXX,"_SRC_AllTimes_Indices_Seasons_",n,"_",nb,".RData",sep=""),quote=FALSE,sep=",")#Save the SRC coeficients

#Plot the results
source(paste(Directory_R_Functions,"Sensitivity/PlotAllSquared.r",sep="")) # See comment on this function
names<-Pa$X #Choose here the parameters to be plotted
xlim<-c(0,t) #Choose the range of times for your plot
ylim<-c(0,1) #Choose the range for y
plotAllSquared(All_SRC,q,names,xlim,ylim)

########################################################################