#function to generate HTML file

getFCs <- function(html) {
  htmlRead <- readLines(html)
  featureTypes <- htmlRead[grep("Feature types", htmlRead)]
  substr(featureTypes, start=21, stop=nchar(featureTypes)-4)
}

library (dismo)
library (raster)
library (rgdal)
library (ENMeval)

lista_variables <- list.files(path="    ",pattern='*.asc', full.names=TRUE) 
variables <- stack(lista_variables)
plot(variables)
presencia<-read.table("    ",header=T, sep=',')
presencia<-presencia[,2:3]
names(presencia)<-c("x","y")
npresencias<-nrow(presencia)
tabla_presencia<-extract(variables, presencia)
head(tabla_presencia)
areaE <-readOGR(dsn="C:/Users/EVOL-ECOL/Desktop/PseudoausencesRhinella/AREA_estudio",layer="ANDES")
mask = raster("C:/Users/EVOL-ECOL/Desktop/PseudoausencesRhinella/mascara/bio_1.asc")
e <- extent(-81.13291, -63.08904, -55.71303, 11.21101)
bg <- randomPoints(mask, 1000, ext=e, excludep=TRUE)
plot(!is.na(mask), legend=FALSE)
plot(e, add=TRUE, col='red')
points(bg, cex=0.5)
tabla_background<-extract(variables, bg)
pb<-c(rep(1, nrow(tabla_presencia)), rep(0, nrow(tabla_background)))
head(pb)
tabla_pb<-data.frame(cbind(pb, rbind(tabla_presencia, tabla_background)))
head(tabla_pb)
enmeval <-ENMevaluate(presencia, variables, bg.coords= bg, RMvalues = seq(1, 10, 0.5), 
                      fc = c("L", "LQ", "H", "LQH", "LQHP", "LQHPT"),
                      n.bg = 5000, method = 'randomkfold',  kfolds=5, 
                      clamp = F, rasterPreds = TRUE, 
                      parallel = TRUE, numCores = 12, progbar = TRUE, algorithm='maxent.jar')

enmeval@results
library (knitr)
Tbl <- enmeval@results
kable (Tbl)   
eval.plot(Tbl)
plot(enmeval@predictions[[which (enmeval@results$delta.AICc == 0) ]])
points(enmeval@occ.pts, pch=21, bg=enmeval@occ.grp)
plot(enmeval@predictions)
enmeval
plot(enmeval@predictions[[which(enmeval@results$delta.AICc==0)]], main="Best model")
aic.opt <- enmeval@models[[which(enmeval@results$delta.AICc==0)]]
aic.opt
aic.opt@results
var.importance(aic.opt)
maxmod <- enmeval@models[[which(enmeval@results$delta.AICc==0)]]
maxmod 
aicmods <- which(enmeval@results$AICc == min(na.omit(enmeval@results$AICc)))
enmeval@results[aicmods,]
BestModel <- enmeval@results[aicmods,]
kable (BestModel)   


aicmods <- which(enmeval@results$AICc == min(na.omit(enmeval@results$AICc)))[1] 
aicmods <- enmeval@results[aicmods,]
FC_best <- as.character(aicmods$features[1]) 
rm_best <- aicmods$rm 


maxent.args <- make.args(RMvalues = rm_best, fc = FC_best) 

mx_best <- maxent(variables, presencia, args=maxent.args[[1]],
                  path = '       ', overight=T)
r_best <- predict(mx_best, variables, overwrite=TRUE, progress = 'text')
plot (r_best)
var.importance(mx_best)
writeRaster(r_best, filename="    ", format="ascii", overwrite=TRUE)

library (ecospat)

datostest <- read.csv("    ", header=TRUE, sep=",")
boyce  <- ecospat.boyce (r_best, datostest ) 
boyce$Spearman.cor


library(kuenm)


model <- raster("   ")
ind_data <- read.csv("   ")
thres <- 5 
rand_perc <- 50 
iterac <- 500 
p_roc <- kuenm_proc(occ.test = ind_data, model = model, threshold = thres,
                    rand.percent = rand_perc, iterations = iterac)
p_roc$pROC_summary  
p_roc$pROC_results