#These are the comands to execute in R to make a Sobol sensitivity analysis of the diferents models
#Sets of parameters are alreay defined for each model. Be sure to have the files of parameters definitions in the correct directory. Format : XXX_Parameters.csv.
#Type of files :
#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_Output<-"/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.

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

#Reading the Parameters
Pa<-read.csv(paste(Directory_Entrie,XXX,"_Parameters_Mu.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("/data/paulb/R/Functions/Random_Param.R") #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_Output
write.table(x$X,file=paste(Directory_Output,XXX,"_Sobol_ParamToBeTested_",n,"_",nb,"_Mu.csv",sep=""),quote=FALSE,sep=",",row.name=FALSE,col.name=FALSE)
save.image(file=paste(Directory_Output,XXX,"_ParamToBeTested_",n,"_",nb,"_UNFINISHED_Mu.RData",sep=""))

##################################################################
##################################################################
##Execute Scilab Functions infected_XXX.sci and 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_Output,XXX,"_ComputedResults_",n,"_",nb,"_NEW.csv",sep=""),header = FALSE,sep = ",")));
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, execute these lines

########################################################################
t<-3300 #Change here the final time
Seasons<-c("winter","spring","summer","Autumn");
Y <-as.matrix(read.csv(paste(Directory_Output,XXX,"_Climatology_ComputedResults_",n,"_",nb,"_NEW.csv",sep=""),header = FALSE,sep = ","));
Xx<-list()
for (k in 1:4) {
	xx <- x
	y<-Y[k,]
	tell(xx, y);
	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],"_NEW.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],"_NEW.csv",sep=""),quote=FALSE,sep=",")#Save the Sobol Total Indices
	print(xx)
	Xx[[k]]<-list()
	Xx[[k]]<-xx	
}
	save.image(file=paste(Directory_Results,XXX,"_Sobol_Indices_Seasons_",n,"_",nb,"_NEW.RData",sep=""))
	source("/data/paulb/R/Functions/plot_sobol_seasons.R")
	source("/data/paulb/R/Functions/nodeplot.R")
	plot_sobol_seasons(Xx)
	
#Plot Locally
load("/Volumes/shared_data/paulb/Results/Results/AM_Sobol_Indices_Seasons_1000_10000_NEW.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)
########################################################################

#If output at all times, execute these lines

########################################################################
t<-9000-2 #Change here the final time
All_Sobol_First<-NULL;
All_Sobol_Total<-NULL;

Y <-as.matrix(read.csv(paste(Directory_Output,XXX,"_ComputedResults.csv",sep=""),header = FALSE,sep = ","));

for (k in 1:t) {
	x <- sobol2002(model = NULL, X1, X2, nboot = nb)
	y<-Y[,1]
	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
	All_Sobol_First<-rbind(All_Sobol_First,x$S);
	All_Sobol_Total<-rbind(All_Sobol_Total,x$T);
}
write.table(All_Sobol_First,file=paste(Directory_Output,XXX,"_Sobol_Indices_First_Order_All_Times.csv",sep=""),quote=FALSE,sep=",")#Save the Sobol First Indices
write.table(All_Sobol_Total,file=paste(Directory_Output,XXX,"_Sobol_Indices_Total_All_Times.csv",sep=""),quote=FALSE,sep=",")#Save the Sobol Total Indices

source("plotAll.r")
indice<-"Total"
names<-c("r","alpha")
plotAll(All_Sobol_Total,p,names,indice)
########################################################################