
#library:------------------------------------------
library(rgrass7)

library(caret)
library(randomForest)
library(e1071)

library(vcd)


#define input variables------------------------------------------

args <- commandArgs(trailingOnly = TRUE)

typevar=args[1]
year=args[2]
layers=args[3]
train=args[4]
valid=args[5]
optimise=args[6]
workdir=args[7]
seasons=args[8]

layers=strsplit(layers, ",")[[1]]

set.seed(975)

#kappa and confusion function-----------------------------------
fun.kappa=function(validate,prediction){
	#confusion
	Kappa_confusion2<-confusionMatrix(prediction,validate ,dnn = c("Predic.", "Reference"))
	#kappa	
	Kappaest<-Kappa(Kappa_confusion2$table)
	Kappaint<-confint(Kappaest)
	#calculate comision, omision
	Kappa_confusion=as.data.frame.matrix(Kappa_confusion2$table)
	Kappa_confusion[,"Total"]<-NA
	Kappa_confusion["Total",]<-NA
	Kappa_confusion[,"ExU"]<-NA
	Kappa_confusion["ExP",]<-NA
	Kappa_confusion[,"CE"]<-NA
	Kappa_confusion["OE",]<-NA
	Kappa_confusion[,"condK"]<-NA
	F_correct<-0
	for (exp in 1:(length(Kappa_confusion[,1])-3)){
		F_correct<-Kappa_confusion[exp,exp]+F_correct
		Kappa_confusion[exp,"Total"]<-sum(Kappa_confusion[exp,1:(length(Kappa_confusion[,1])-3)])
		Kappa_confusion["Total",exp]<-sum(Kappa_confusion[1:(length(Kappa_confusion[,1])-3),exp])
		Kappa_confusion[exp,"ExU"]<-round(Kappa_confusion[exp,exp]/Kappa_confusion[exp,"Total"]*100,digits=3)
		Kappa_confusion["ExP",exp]<-round(Kappa_confusion[exp,exp]/Kappa_confusion["Total",exp]*100,digits=3)
		Kappa_confusion[exp,"CE"]<-round(((Kappa_confusion[exp,"Total"]-Kappa_confusion[exp,exp])/Kappa_confusion[exp,"Total"])*100,digits=3)
		Kappa_confusion["OE",exp]<-round(((Kappa_confusion["Total",exp]-Kappa_confusion[exp,exp])/Kappa_confusion["Total",exp])*100,digits=3)
	}
	Kappa_confusion["Total","Total"]<-sum(Kappa_confusion[1:(length(Kappa_confusion[,1])-3),"Total"])
	#Accuracy:
	F_global<-F_correct/Kappa_confusion["Total","Total"]*100
	#Conditional kappa
	for (exp in 1:(length(Kappa_confusion[,1])-3)){
		Kappa_confusion[exp,"condK"]=round(((Kappa_confusion["Total","Total"]*Kappa_confusion[exp,exp])-(Kappa_confusion[exp,"Total"]*Kappa_confusion["Total",exp])) / ((Kappa_confusion["Total","Total"]*Kappa_confusion[exp,"Total"])-(Kappa_confusion[exp,"Total"]*Kappa_confusion["Total",exp]))*100,digits=3)
	}
	for (exp in 1:(length(Kappa_confusion[1,])-4)){
		Kappa_confusion["condK",exp]=round(((Kappa_confusion["Total","Total"]*Kappa_confusion[exp,exp])-(Kappa_confusion["Total",exp]*Kappa_confusion[exp,"Total"])) / ((Kappa_confusion["Total","Total"]*Kappa_confusion["Total",exp])-(Kappa_confusion["Total",exp,]*Kappa_confusion[exp,"Total"]))*100,digits=3)
	}
	#results
	tablak<-data.frame(t(rep(NA,5)));tablak<-tablak[-1,]
	names(tablak)<-c("Accuracy","Kappa","Int1K","Int2K","ASE");
	tablak["Resum",]<-c(round(F_global,digit=4),round(as.vector(Kappaest$Unweighted[1]),digit=4),round(as.vector(Kappaint[1]),digit=4),round(as.vector(Kappaint[3]),digit=4),round(as.vector(Kappaest$Unweighted[2]),digit=4));

	#export
	return=list(Kappa_confusion,tablak)
	return(return)
}


#input data------------------------------------------

#raster
RASTER <- readRAST(c(train,valid,layers))
RASTER[[train]] <- as.factor(RASTER[[train]])
RASTER[[valid]] <- as.factor(RASTER[[valid]])
#
datRASTER <- data.frame(RASTER@data)
datcal <- datRASTER[!is.na(datRASTER[,train]), c(train,layers)]
#
datval <- datRASTER[!is.na(datRASTER[,valid]), c(valid,layers)]


#define formula for models----------------------------

formulaag <- as.formula(paste(train," ~ ",paste0(layers,collapse="+")))

#define matrices for results and folder to results----------------------------

#folder to save results:
if(dir.exists(paste0(workdir,"/clasifresults"))==F){dir.create(paste0(workdir,"/clasifresults"))}

#summary exists?, if not create
if(file.exists(paste0(workdir,"/summaryCAL.Rdata"))==F){
	tableA <- data.frame(typevar=character(0),clasif=character(0),NumVar=integer(0),param1=double(0),param2=double(0),Accuracy=double(0),Kappa=double(0),Optim=integer(0),stringsAsFactors=FALSE)
}else{
	load(paste0(workdir,"/summaryCAL.Rdata"))
}
if(file.exists(paste0(workdir,"/summaryVAL.Rdata"))==F){
	tableK <- data.frame(typevar=character(0),clasif=character(0),Accuracy=double(0),Kappa=double(0),Int1k=double(0),Int2k=double(0),ASE=double(0),Optim=integer(0),stringsAsFactors=FALSE)
}else{
	load(paste0(workdir,"/summaryVAL.Rdata"))
}

#define for calibration---------------------------

if(optimise==1){
	controlObject <- trainControl(method="repeatedcv",repeats=5,number=10)
}


############################################################
#RANDOMFOREST
#############################################################

if(optimise==1){
	#fit
	fitRF <- train(formulaag,data=datcal,method="rf",ntree=500,trControl=controlObject,fit=F,importance=T,tuneGrid=data.frame(.mtry=c(1:15)),na.action=na.omit)

	con2tol=oneSE(fitRF$results,metric="Kappa",num=10,maximize=T)
	RFparam=data.frame(ntree=500,mtry=fitRF$results[con2tol,1:5][["mtry"]])
	#fit
	fit2RF <- randomForest(formulaag,data=datcal,importance=T,na.action=na.omit,mtry=RFparam[["mtry"]],ntree=RFparam[["ntree"]])

	#plot optim
	png(paste0(workdir,"/clasifresults/clasif_rf",typevar,year,seasons,"_OPTIM.png"),width=640,height=480)
		plot(fitRF$results$mtry,fitRF$results$Kappa,type="l",col="lightblue",ylab="Kappa (Repeated Cross-Validation)",xlab="Randomly Selectred Predictors (mtry)")
		grid(col="grey")
		points(fitRF$results$mtry,fitRF$results$Kappa,pch=16,col="lightblue")
		points(RFparam$mtry,subset(fitRF$results$Kappa,fitRF$results$mtry==RFparam$mtry),col="red",pch=16)
		text(RFparam$mtry,subset(fitRF$results$Kappa,fitRF$results$mtry==RFparam$mtry),col="darkred",labels=paste("mtry: ",RFparam$mtry,"\n ntree: ",RFparam$ntree,"\n Kappa: ",round(subset(fitRF$results$Kappa,fitRF$results$mtry==RFparam$mtry),4),sep=""),pos=1,offset=1.5)
	dev.off()

}else{
	#fit
	fit2RF <- randomForest(formulaag,data=datcal,importance=T,na.action=na.omit)
}

#predict and export to GRASS
datRASTER$pred2RF <- predict(fit2RF,newdata=datRASTER[layers],na.action=na.exclude)
RASTER@data$predRF=as.integer(datRASTER$pred2RF)
#
writeRAST(RASTER, paste0("clasif_rf",typevar,year,seasons), zcol = "predRF", NODATA=NULL, ignore.stderr = FALSE,overwrite=TRUE)

#results------------------------------------------------------------

#plot mtree
png(paste0(workdir,"/clasifresults/clasif_rf",typevar,year,seasons,"_ntree.png"),width=640,height=480)
	plot(fit2RF,main=paste0("RF ntree, ",typevar,", ",year))	
	grid(col="grey")
	legend("topright", colnames(fit2RF$err.rate),col=1:length(ncol(fit2RF$err.rate)),cex=0.8,fill=1:ncol(fit2RF$err.rate))	
dev.off()
#plot importance
png(paste0(workdir,"/clasifresults/clasif_rf",typevar,year,seasons,"_importance.png"),width=640,height=480)
	rfImp1=importance(fit2RF,2)
	rownames(rfImp1)=substring(rownames(fit2RF$importance),2)


par(las=2)
par(mar=c(8,8,1,1)) # adjust as needed

	plot(sort(rfImp1),type="l",lty=3,col="black",ylab=colnames(rfImp1),xlab="",main=paste0("RF importance, ",typevar,", ",year),xaxt="n")
	points(sort(rfImp1),col="blue",pch=16,cex=1)
	axis(1,labels=rownames(rfImp1)[order(rfImp1)],at=1:length(order(rfImp1)),cex.axis=1,las=2,xlab="")
	grid(col="grey")
dev.off()

#results:
datcal$RFpredicted <- predict(fit2RF,newdata=datcal[layers],na.action=na.exclude)
adjustRF <- confusionMatrix(datcal[,train],datcal$RFpredicted ,dnn = c("Predic.", "Reference"))
#
tableA[nrow(tableA)+1,] <- c(typevar,"RF",length(layers),fit2RF[["ntree"]],fit2RF[["mtry"]],round(adjustRF$overall[["Accuracy"]],4),round(adjustRF$overall[["Kappa"]],4),optimise)

#validation-----------------------------------------------------------
datval$RFpred <- predict(fit2RF, datval[layers], na.action = na.exclude)

kappadata=fun.kappa(validate=datval[,valid],prediction=datval[,"RFpred"])
#
tableK[nrow(tableK)+1,]<-c(typevar,"RF",kappadata[[2]][[1]],kappadata[[2]][[2]],kappadata[[2]][[3]],kappadata[[2]][[4]],kappadata[[2]][[5]],optimise)

#export confusion table:
write.table(kappadata[[1]],file=paste0(workdir,"/clasifresults/confusion_rf",typevar,year,seasons,".txt"),sep=";")

############################################################
#SVM
#############################################################

if(optimise==1){

	#ajust
	fitSVM=train(formulaag,data=datcal,method="svmRadial",tuneLength=10,trControl=controlObject,fit=F,na.action=na.omit)

	con2tol=oneSE(fitSVM$results,metric="Kappa",num=10,maximize=T)
	SVMparam=data.frame(C=fitSVM$results[con2tol,1:6][["C"]],sigma=fitSVM$results[con2tol,1:6][["sigma"]])
	#fit
	fit2SVM <- svm(formulaag,data=datcal,kernel="radial",fitted=T,na.action=na.omit,gamma =SVMparam$sigma,cost=SVMparam$C)

	#plot optim
	png(paste0(workdir,"/clasifresults/clasif_svm",typevar,year,seasons,"_OPTIM.png"),width=640,height=480)
		plot(fitSVM$results$C,fitSVM$results$Kappa,type="l",col="lightblue",xlab="Cost",ylab="Kappa (Repeated Cross-Validation)")
		grid(col="grey")
		points(fitSVM$results$C,fitSVM$results$Kappa,pch=16,col="lightblue")
		points(SVMparam$C,subset(fitSVM$results$Kappa,fitSVM$results$C==SVMparam$C),col="red",pch=16)
		text(SVMparam$C,subset(fitSVM$results$Kappa,fitSVM$results$C==SVMparam$C),col="darkred",labels=paste("C: ",SVMparam$C,"\n Sigma: ",round(SVMparam$sigma,4),"\n Kappa: ",round(subset(fitSVM$results$Kappa,fitSVM$results$C==SVMparam$C),4),sep=""),pos=1,offset=1.5)
	dev.off()

}else{
	#fit
	fit2SVM <- svm(formulaag,data=datcal,kernel="radial",fitted=T,na.action=na.omit)
}

#predict and export to GRASS
datRASTER$pred2SVM <- predict(fit2SVM,newdata=datRASTER[layers],na.action=na.exclude)
RASTER@data$predSVM=as.integer(datRASTER$pred2SVM)
#
writeRAST(RASTER, paste0("clasif_svm",typevar,year,seasons), zcol = "predSVM", NODATA=NULL, ignore.stderr = FALSE,overwrite=TRUE)

#results------------------------------------------------------------
datcal$SVMpredicted <- predict(fit2SVM,newdata=datcal[layers],na.action=na.exclude)
adjustSVM <- confusionMatrix(datcal[,train],datcal$SVMpredicted ,dnn = c("Predic.", "Refeference"))
#
tableA[nrow(tableA)+1,] <- c(typevar,"SVM",length(layers),round(fit2SVM[["cost"]],3),round(fit2SVM[["gamma"]],4),round(adjustSVM$overall[["Accuracy"]],4),round(adjustSVM$overall[["Kappa"]],4),optimise)

#validation-----------------------------------------------------------
datval$SVMpred <- predict(fit2SVM, datval[layers], na.action = na.exclude)

kappadata=fun.kappa(validate=datval[,valid],prediction=datval[,"SVMpred"])
#
tableK[nrow(tableK)+1,]<-c(typevar,"SVM",kappadata[[2]][[1]],kappadata[[2]][[2]],kappadata[[2]][[3]],kappadata[[2]][[4]],kappadata[[2]][[5]],optimise)

#export confusion table:
write.table(kappadata[[1]],file=paste0(workdir,"/clasifresults/confusion_svm",typevar,year,seasons,".txt"),sep=";")

############################################################
#EXPORT SUMMARY DATA

write.table(tableA,file=paste0(workdir,"/summaryCAL",year,seasons,".txt"),sep=";",,row.names=F)
write.table(tableK,file=paste0(workdir,"/summaryVAL",year,seasons,".txt"),sep=";",,row.names=F)

save(tableA,tableK,file=paste0(workdir,"/summaryCAL.Rdata"))
save(tableA,tableK,file=paste0(workdir,"/summaryVAL.Rdata"))
