rm(list = ls())

#Comparison function
#With this code we perform the binary comparisons. 
#First to perform the control and then to evaluate the transferability.

# Comparison function between two binary maps
#item{0 no habitat in neither two times}
#'           \item{1 habitat loss}
#'           \item{10 habitat gain}
#'           \item{11 preserved habitat} 

#Comparation binaries (compare two binaries)
reportFuture <- function(raster_present_bin, raster_future_bin){
  reclass_matrix <- matrix(c(0, 1, 0, 10), ncol=2)
  offset_raster <- raster::reclassify(raster_future_bin, reclass_matrix)
  result <- raster_present_bin + offset_raster
    return(list(stats=raster::freq(result, useNA="no"),
              raster_change=result))
}

#-------------------------------GEOGRAPHIC COMPARISON
reportChange <- function(values){
  change_stats<-as.data.frame(t(values))
  colnames(change_stats)<-c("V0","V1","V10","V11")
  change_stats$IO <-(change_stats$V11/(change_stats$V1+change_stats$V11))
  change_stats$TSS <-(((change_stats$V11/(change_stats$V1+change_stats$V11))) + 
                        ((change_stats$V0/(change_stats$V0+change_stats$V10)))) -1
  #FPR (Overfitting)
  change_stats$FPR <-((change_stats$V10/(change_stats$V10+change_stats$V0)))
    #FNR (Overestimation)
  change_stats$FNR <-((change_stats$V1/(change_stats$V11+change_stats$V1)))
  return(cbind.data.frame(specie,change_stats[2,]))
  }

#----------------------------Analysis
setwd("C:/requisito/scrip")
#install.packages("devtools")
library("devtools")
#install_github("danlwarren/ENMTools") #seleccionar la opcion 3
library("ENMTools")
library("rgdal", quietly = TRUE)
library("readr", quietly = TRUE)
library("raster", quietly = TRUE)
library("dismo")
library("rowr")

#-------------------DATA: ALL COMPARASION
#folders
setwd<-"C:/requisito/sp_todas/thoatr"
dir.create("C:/requisito/modelos_tvar/thoatr/analisis")
outputFolder<-"C:/requisito/modelos_tvar/thoatr/analisis"

#data base
base<-read.csv("C:/requisito/sp_todas/thoatr/interR/p2_train.csv")

#species name
specie<-as.character(base[1,1])
specie

#---------------------------------MODELS--------
#GLM (CCHANGE MODEL AND SPECIE)
#Read binaries models
glm_Bp2p2<-raster("C:/requisito/modelos_tvar/thoatr/glmp2p2bin_q10.tif") #same time data and predictors (binary: auto 1)
glm_Bp2p3<-raster("C:/requisito/modelos_tvar/thoatr/glmp2p2_3bin_q10.tif") #data one time, predictors other time (binary: cross 1)

glm_Bp3p3<-raster("C:/requisito/modelos_tvar/thoatr/glmp3p3bin_q10.tif") #same time data and predictors (binary: auto 2)
glm_Bp3p2<-raster("C:/requisito/modelos_tvar/thoatr/glmp3p3_2bin_q10.tif") #data one time, predictors other time (binary: cross 2)

#---------------------------------
#------------------------Geographic comparation: cross-temporal approach
#p2p3 (Control model 1)
proyeccion_glmp2p3 <- reportFuture(glm_Bp2p2, glm_Bp2p3)
#writeRaster(proyeccion_glmp2p3$raster_change,'C:/requisito/modelos_tvar/thoatr/analisis/glmp2p2_3.tif', overwrite=TRUE)
plot(proyeccion_glmp2p3$raster_change)

#p3p2 (Control model 2)
proyeccion_glmp3p2 <- reportFuture(glm_Bp3p3, glm_Bp3p2)
#writeRaster(proyeccion_glmp3p2$raster_change,"C:/requisito/modelos_tvar/thoatr/analisis/glmp3p3_2.tif", overwrite=TRUE)
plot(proyeccion_glmp3p2$raster_change)

#glm
change_glmp2p3<-reportChange(proyeccion_glmp2p3$stats)
change_glmp3p2<-reportChange(proyeccion_glmp3p2$stats)
head(cambio_glmp2p3)

cross1<-cbind.data.frame(change_glmp2p3, Time="Cross1", model="glm")
cross2<-cbind.data.frame(change_glmp3p2, Time="Cross2", model="glm")

analysis_cross<-cbind.data.frame (cross1, cross2)
#write.csv (analysis_cross,"C:/requisito/comp_s/glmp3p2_thoatr.csv")

#-------------------------------Transferability Comparison
#---------------------------------------Binaries
# Forecast
glm_Bp2p2<-raster("C:/requisito/sp_todas/thoatr/interR/glmp2p2bin_q10.tif") # Evaluation
glm_Bp2p3<-raster("C:/requisito/sp_todas/thoatr/interR/glmp2p3bin_q10.tif") # Transference

# Hindcast
glm_Bp3p3<-raster("C:/requisito/sp_todas/thoatr/interR/glmp3p3bin_q10.tif") #Evaluation
glm_Bp3p2<-raster("C:/requisito/sp_todas/thoatr/interR/glmp3p2bin_q10.tif") #Transference

#---------------------------Geographic comparation: transferability comparison
# Forecast
proyeccion_glmp2p3 <- reportFuture(glm_Bp2p2, glm_Bp2p3)
#writeRaster(proyeccion_glmp2p3$raster_change,'C:/requisito/sp_todas/thoatr/interR/analisis/glmp2p3.tif', overwrite=TRUE)
plot(proyeccion_glmp2p3$raster_change)

# Hindcast
proyeccion_glmp3p2 <- reportFuture(glm_Bp3p3, glm_Bp3p2)
#writeRaster(proyeccion_glmp3p2$raster_change,"C:/requisito/sp_todas/thoatr/interR/analisis/glmp3p2.tif", overwrite=TRUE)
plot(proyeccion_glmp3p2$raster_change)

# transferability report
change_glmp2p3<-reportChange(proyeccion_glmp2p3$stats) #Forecast
change_glmp3p2<-reportaChange(proyeccion_glmp3p2$stats) #Hindcast
head(cambio_glmp2p3)

forecast<-cbind.data.frame(change_glmp2p3, Time="Forecast", model="glm")
hindcast<-cbind.data.frame(change_glmp3p2, Time="Hindcast", model="glm")

analysis_transferability<-cbind.data.frame (forecast, hindcast)
#write.csv (transferability,"C:/requisito/comp_s/glmp3p2_thoatr.csv")