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


############Environmental data###################
library (raster)

biogeo7<-raster("biogeo07_30s.asc")
biogeo10<-raster("biogeo10_30s.asc")
biogeo12<-raster("biogeo12_30s.asc")
biogeo13<-raster("biogeo13_30s.asc")
biogeo17<-raster("biogeo17_30s.asc")

###########Stack predictors###########################
predictors<-stack(biogeo7,biogeo10,biogeo12,biogeo13,biogeo17)
plot(predictors)


occ <- read.csv('anemona_MARSPEC.csv', header=TRUE, sep=',')[,-1]


###################################################  
#### bg con sesgo
library (dismo)
# get the file names
files <- list.files(path='xxxxxxxxxxxxxxxxxxxxx',  pattern='asc',  full.names=TRUE )
####  we use the first file to create a RasterLayer
mask <- raster(files)
plot (mask)
set.seed(1963)
bg <- randomPoints(mask, 1000, prob=T)
plot (bg)
plot (mask)
points(bg)
########################################################################  

dir.create('ENMeval_ENERO') # creates a folder to store the results
 
###k-1 Jackknife
jack <- get.jackknife(occ, bg)
plot(bg, col='gray', legend=FALSE)
points(occ, pch=21, bg=jack$occ.grp)  # note that colors are repeated here


library (ENMeval)

enmeval_results <- ENMevaluate(occ, predictors, occ.grp = jack$occ.grp, 
            bg.grp = jack$bg.grp, RMvalues = seq(1, 3, 0.5), 
            fc = c("L", "LQ"), method = 'jackknife',
            clamp = TRUE, rasterPreds = TRUE, 
            parallel = TRUE, numCores = 12, algorithm='maxent.jar', progbar = TRUE )

eval2.par <- ENMevaluate(occ, predictors, bg, method='jackknife',
            RMvalues=c(1,2), fc=c('L','LQ'), parallel=TRUE, numCores = 12, algorithm='maxent.jar')


data(enmeval_results)
enmeval_results


enmeval_results@results### See table of evaluation metrics

### Plot prediction with lowest AICc
plot(enmeval_results@predictions[[which (enmeval_results@results$delta.AICc == 0) ]])
points(enmeval_results@occ.pts, pch=21, bg=enmeval_results@occ.grp)
plot(enmeval_results@predictions)

#Maps
maps <- enmeval_results@predictions
maps



library (knitr)  ###beautiful table # unpack the results data frame, the list of models, and the RasterStack of raw predictions
evalTbl <- enmeval_results@results
kable (evalTbl)   
eval.plot(evalTbl)


evalMods <- enmeval_results@models
evalMods

evalPreds <- enmeval_results@predictions
evalPreds 
plot (evalPreds)


# Results

enmeval_results
plot(enmeval_results@predictions[[which(enmeval_results@results$delta.AICc==0)]], main="Best model")


aic.opt <- enmeval_results@models[[which(enmeval_results@results$delta.AICc==0)]]
aic.opt
aic.opt@results
var.importance(aic.opt)


#### curves
maxmod <- enmeval_results@models[[which(enmeval_results@results$delta.AICc==0)]]
maxmod 


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

head(enmeval_results@results)

# parsimonious model
aicmods <- which(enmeval_results@results$AICc == min(na.omit(enmeval_results@results$AICc)))
enmeval_results@results[aicmods,]


###Plot ENMeval models comparisons:
par(mfrow=c(2,3))
eval.plot(enmeval_results@results)
eval.plot(enmeval_results@results, 'Mean.AUC', legend = F)
eval.plot(enmeval_results@results, 'Mean.AUC.DIFF', legend = F)
eval.plot(enmeval_results@results, 'Mean.OR10', legend = F)
eval.plot(enmeval_results@results, 'Mean.ORmin', legend = F)
plot(enmeval_results@results$Mean.AUC, enmeval_results@results$delta.AICc, bg=enmeval_results@results$features, pch=21, cex= enmeval_results@results$rm/2, xlab = "Mean.AUC", ylab = 'delta.AICc', cex.lab = 1.5)
legend("topright", legend=unique(enmeval_results@results$features), pt.bg=enmeval_results@results$features, pch=21)
mtext("Circle size proportional to regularization multiplier value")

#####most parsimonious model to build a new Maxent model.

aicmods <- which(enmeval_results@results$AICc == min(na.omit(enmeval_results@results$AICc)))[1] # AIC model
aicmods <- enmeval_results@results[aicmods,]
FC_best <- as.character(aicmods$features[1]) # Get FCs from the AIC model
rm_best <- aicmods$rm # Get RM from the AIC model

# Built the model
maxent.args <- make.args(RMvalues = rm_best, fc = FC_best) # Define arguments

mx_best <- maxent(predictors, occ, args=maxent.args[[1]],
                  path = 'I:/MS_pinochet_NOV_2018/SoloMARSPEC/ENMeval_ENERO', overight=T)
r_best <- predict(mx_best, predictors, overwrite=TRUE, progress = 'text')


plot (r_best)