setwd("F:/Ioana M/Ioana M/Articol/2025/MAx Entropy/mca")

output_folder <- "maxent_output_final"
dir.create(output_folder, showWarnings = FALSE)

library(raster)
library(dismo)
library(sf)
library(terra)

# A. Data


presence_data   <- read.csv("positive_points.csv", encoding = "UTF-8")
background_data <- read.csv("false_points.csv", encoding = "UTF-8")

colnames(presence_data)[colnames(presence_data) %in% c("xcoord","ycoord","x","y")] <- c("X","Y")
colnames(background_data)[colnames(background_data) %in% c("xcoord","ycoord","x","y")] <- c("X","Y")

p_coords  <- presence_data[, c("X", "Y")]
bg_coords <- background_data[, c("X", "Y")]

raster_files <- c(
  dem    = "Reproject.tif",
  slope  = "Slope.tif",
  aspect = "aspect.tif",
  tpi    = "TPI.tif",
  water  = "Distanta.tif"
)

predictors <- stack(raster_files)
predictors <- scale(predictors)
names(predictors) <- names(raster_files)


# B. Pearson matrix

predictors_terra <- rast(predictors)

set.seed(123)
sample_data <- as.data.frame(
  spatSample(predictors_terra, size = 10000, na.rm = TRUE)
)

cor_matrix <- cor(sample_data, method = "pearson", use = "complete.obs")

print(round(cor_matrix, 3))


# C. MODEL 


maxent_model <- maxent(
  x = predictors,
  p = p_coords,
  a = bg_coords,
  args = c(
    "betamultiplier=1.5",           
    "linear=true",                  
    "quadratic=true",               
    "product=false",
    "hinge=false",
    "threshold=false",
    "replicatetype=crossvalidate",  
    "replicates=5",
    "maximumiterations=1000",
    "responsecurves=true",
    "jackknife=true",
    "outputformat=logistic"
  ),
  path = output_folder
)



# D. Results

results <- maxent_model@results
training_auc <- mean(results[grep("Training.AUC", rownames(results)), ])
cat("Training AUC:", round(training_auc, 3), "\n")
if (any(grepl("Test.AUC", rownames(results)))) {
  test_auc <- mean(results[grep("Test.AUC", rownames(results)), ])
  cat("Test AUC:", round(test_auc, 3), "\n")
}

auc <- results[grep("Training.AUC", rownames(results)), ]


var_contrib <- results[grep("contribution", rownames(results)), ]
cat("\nContribuția variabilelor (%):\n")
print(sort(var_contrib, decreasing = TRUE))



# E. PREDICTION


prediction <- predict(maxent_model, predictors, 
                      filename = paste0(output_folder, "/archaeological_prediction.tif"), 
                      overwrite = TRUE)

threshold <- quantile(extract(prediction, p_coords), 0.1, na.rm = TRUE)
binary <- prediction > threshold
writeRaster(binary, paste0(output_folder, "/binary_suitability.tif"), overwrite = TRUE)






â