# 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)
#		-- package "sensitivity" (available on CRAN)
#		-- R function "Random_Param" to randomly generate sets of parameters
#		-- files "XXX_Parameters.csv" (format described below)
#		-- R functions "plot_sobol_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

#Package "sensitivity" needed (available on CRAN)
library(sensitivity,lib.loc="/home/paulb/R/x86_64-redhat-linux-gnu-library/2.15");

#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/Sensitivity/" #Directory where you want your output to be stocked. Has to match with the entry directory of Scilab.
Directory_Results<-"/data/paulb/Results/Results/Sensitivity/" #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 (number of bootstrap replicates)
########################################################################

#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);

#Generate the vectors of parameters to be tested by Scilab
source(paste(Directory_R_Functions,"Sensitivity/Random_Param.R",sep="")) #see comments in file "Random_Param.r"
X1 <- data.frame(Random_Param(n,param_min,param_max,Dist)) 
X2 <- data.frame(Random_Param(n,param_min,param_max,Dist)) # The method of sobol requires 2 samples
x <- sobol2002(model = NULL, X1, X2, nboot = nb) #Lauch the function sobol2002. Result : x$X is a matrix, each line is a set of parameters to be tested by scilab.

# Write the sets of parameters to be tested in Directory_Communication
write.table(x$X,file=paste(Directory_Communication,XXX,"_Sobol_ParamToBeTested_",n,"_",nb,".csv",sep=""),quote=FALSE,sep=",",row.name=FALSE,col.name=FALSE)
save.image(file=paste(Directory_R_Data,XXX,"_ParamToBeTested_",n,"_",nb,"_UNFINISHED.RData",sep="")) # Save the R objects generated

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

#If only one time for the output, execute these lines

########################################################################
t<-9000 #Change here the final time
y <-as.vector(as.matrix(read.csv(paste(Directory_Communication,XXX,"_ComputedResults_",n,"_",nb,".csv",sep=""),header = FALSE,sep = ","))); # Read the answers given by Scilab
tell(x, y);
rownames(x$S)<-Pa$X; #Give their names to the parameters.
rownames(x$T)<-Pa$X; #Give their names to the parameters.
write.table(x$S,file=paste(Directory_Results,XXX,"_Sobol_Indices_First_Order_time_",t,"_",n,"_",nb,".csv",sep=""),quote=FALSE,sep=",")
write.table(x$T,file=paste(Directory_Results,XXX,"_Sobol_Indices_Total_time_",t,"_",n,"_",nb,".csv",sep=""),quote=FALSE,sep=",")
print(x) #Display the results
plot(x)
########################################################################

#If output at several times (seasons), execute these lines

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

t<-3300 #Change here the final time
Seasons<-c("winter","spring","summer","Autumn");
Y <-as.matrix(read.csv(paste(Directory_Communication,XXX,"_Climatology_ComputedResults_",n,"_",nb,".csv",sep=""),header = FALSE,sep = ",")); # Read the answers given by Scilab
Xx<-list()
for (k in 1:4) {
	xx <- x
	y<-Y[k,]
	tell(xx, y); # See documentation of package sensitivity for explanation on function "tell"
	rownames(xx$S)<-Pa$X; #Give their names to the parameters.
	rownames(xx$T)<-Pa$X; #Give their names to the parameters
	write.table(xx$S,file=paste(Directory_Results,XXX,"_Sobol_Indices_First_Order_",n,"_",nb,"_",Seasons[k],".csv",sep=""),quote=FALSE,sep=",")#Save the Sobol First Indices
	write.table(xx$T,file=paste(Directory_Results,XXX,"_Sobol_Indices_Total_",n,"_",nb,"_",Seasons[k],".csv",sep=""),quote=FALSE,sep=",")#Save the Sobol Total Indices
	print(xx)
	Xx[[k]]<-list()
	Xx[[k]]<-xx	
}
	save.image(file=paste(Directory_R_Data,XXX,"_Sobol_Indices_Seasons_",n,"_",nb,".RData",sep=""))
	source(paste(Directory_R_Functions,"Sensitivity/plot_sobol_seasons.R",sep="")) # See comment on this functions
	source(paste(Directory_R_Functions,"Sensitivity/nodeplot.R",sep="")) # Unmodified function from package "sensitivity"
	plot_sobol_seasons(Xx)
	
#Plot Locally
load("/Volumes/shared_data/paulb/Results/Results/ABP_Sobol_Indices_Seasons_1000_10000.RData")
source("/Volumes/shared_data/paulb/R/Functions/plot_sobol_seasons.R")
source("/Volumes/shared_data/paulb/R/Functions/nodeplot.R")
plot_sobol_seasons(Xx)