# ------------------------------------------------------------------
# 🔹 Load Environmental Variables (Pre-processed and Masked)
# ------------------------------------------------------------------
current_clim_m <- rast("Cholera/Climate/Mask/current_clim_m.nc")

# ------------------------------------------------------------------
# 🔹 Summarize Raster Layers to Detect Issues (NA, Zero Values, Outliers)
# This helps identify missing or invalid data before proceeding with the model
# ------------------------------------------------------------------
summary(current_clim_m)

# ------------------------------------------------------------------
# 🔹 Identify and Remove Problematic Variables
# Some layers may contain excessive NA values, zero values, or unrelated variables.
# We remove these to avoid errors or bias in the model.

# ------------------------------------------------------------------
filtered_clim_m <- current_clim_m[[-c(5, 7, 11, 17, 18, 20)]]

# Print the names of the remaining raster layers to verify the selection
names(filtered_clim_m)

# ------------------------------------------------------------------
# 🔹 Scale Environmental Variables (Standardization)
# Scaling transforms each variable to have a mean of 0 and standard deviation of 1.
# This ensures variables with different ranges (e.g., temperature vs. nutrient levels)
# do not disproportionately influence the model.
# ------------------------------------------------------------------
scaled_clim_m <- lapply(1:nlyr(filtered_clim_m), function(i) {
  raster::raster(scale(filtered_clim_m[[i]]))  # Standardize each raster layer separately
})

# Convert the list of scaled layers to a RasterStack for compatibility with modeling
current_clim_scale <- raster::stack(scaled_clim_m)

# ------------------------------------------------------------------
# 🔹 Visualize Scaled Environmental Variables
# This step ensures that all layers were successfully scaled and are ready for analysis.
# ------------------------------------------------------------------
plot(current_clim_scale)

# Load occurrence data
cholera <- read.csv("Cholera/data/Occurrence_crop.csv")

# Add species column
cholera$species <- "Cholera"

# Select relevant columns
cholera <- cholera %>% 
  dplyr::select("x", "y", "species")

# Thin the occurrence data to reduce spatial autocorrelation
thinned_cholera <- thin(loc.data = cholera,
                        lat.col = "y",
                        long.col = "x",
                        spec.col = "species",
                        thin.par = 10,  # Minimum distance threshold in km
                        reps = 100,     # Number of repetitions
                        locs.thinned.list.return = TRUE,
                        write.files = TRUE,
                        max.files = 1,
                        out.dir = "Cholera/thinning",
                        out.base = "cholera_thinned", 
                        write.log.file = TRUE)

# Load the thinned dataset
thinned_cholera_data <- read.csv("Cholera/thinning/cholera_thinned_thin1.csv")

# Convert thinned data to spatial vector
thinned_cholera_vec <- vect(thinned_cholera_data, geom = c("x", "y"),
                            crs = "+proj=longlat +datum=WGS84 +no_defs", 
                            keepgeom = FALSE)

# Extract climate data at occurrence locations
thinned_cholera_clim <- terra::extract(rast(current_clim_scale), 
                                       thinned_cholera_vec,
                                       xy = TRUE) %>%
  na.omit() %>%
  dplyr::select("x", "y") %>%
  mutate(species = "Cholera")

# ------------------------------------------------------------------
# 🔹 Compute Correlation Matrix for Environmental Variables
# ------------------------------------------------------------------
# Convert the RasterStack to a DataFrame before computing correlation
M <- cor(as.data.frame(rast(current_clim_scale)))

# ------------------------------------------------------------------
# 🔹 Visualize the Correlation Matrix as a Heatmap
# This helps identify highly correlated variables, which may cause redundancy in modeling.
# 'hc.order = TRUE' orders variables hierarchically for better interpretation.
# 'type = "lower"' shows only the lower triangle of the correlation matrix.
# 'lab = TRUE' adds correlation values inside the heatmap.
# ------------------------------------------------------------------
ggcorrplot(M, 
           hc.order = TRUE, 
           type = "lower", 
           lab = TRUE, 
           lab_size = 2,
           title = "Correlation Heatmap of Environmental Variables")

# ------------------------------------------------------------------
# 🔹 Compute Variance Inflation Factor (VIF) to Detect Multicollinearity
# ------------------------------------------------------------------
# VIF measures how much a variable is correlated with other predictor variables.
# High VIF (>10) indicates multicollinearity, which can lead to unreliable model predictions.
# We set a threshold (th = 5) to remove highly correlated variables.
# ------------------------------------------------------------------
current_clim_vif <- vifstep(rast(current_clim_scale), th = 5)

# Print the VIF results
current_clim_vif

# ------------------------------------------------------------------
# 🔹 Check Selected Variables After Removing Highly Correlated Ones
# This step ensures that only independent variables are retained for modeling.
# ------------------------------------------------------------------
selected_variables <- current_clim_vif@results$Variables

# Print the final list of selected variables
selected_variables

# Set a random seed for reproducibility
set.seed(1234)

# ------------------------------------------------------------------
# 🔹 Split Data into Training (70%) and Testing (30%) Sets
# ------------------------------------------------------------------
# This step ensures that the model is trained on a subset of data (70%)
# and tested on the remaining 30% to evaluate its predictive performance.
# 
# 'method = "random"' randomly selects training and testing points.
# 'longitude' and 'latitude' specify spatial coordinates.
# 'train_proportion = 0.70' sets the training data percentage.
# 'save = FALSE' prevents saving the split dataset as a file.
# 'name = "occ"' assigns a name to the occurrence dataset.
# ------------------------------------------------------------------
data_split <- ellipsenm::split_data(thinned_cholera_clim,
                                    method = "random",
                                    longitude = "x",
                                    latitude = "y",
                                    train_proportion = 0.70,
                                    save = FALSE,
                                    name = "occ")

# Define the predictor variables
all_variables <- c(current_clim_vif@results$Variables)

# Function to generate three possible combinations
generate_combinations <- function(vars) {
  comb_list <- list()
  for (m in 3) {
    combs <- combn(vars, m, simplify = FALSE)
    comb_list <- c(comb_list, combs)
  }
  return(comb_list)
}

# Generate variable sets
all_combinations <- generate_combinations(all_variables)
all_combinations

# Prepare sets for ellipsoidal models
variable_sets <- prepare_sets(raster::stack(current_clim_scale), 
                              all_combinations)

# ------------------------------------------------------------------
# 🔹 Define Ellipsoid Calibration Method
# ------------------------------------------------------------------
# 'covmat' (covariance matrix) is a method used to build ellipsoid niche models.
# It assumes the ecological niche can be represented by a multidimensional
# covariance matrix of environmental variables.
# ------------------------------------------------------------------
methods <- c("covmat")

# ------------------------------------------------------------------
# 🔹 Run Ellipsoid Model Calibration
# ------------------------------------------------------------------
# This step calibrates the model based on training data and selected variables.
# 
# - 'data = data_split' -> Uses the previously split training dataset.
# - 'species = "species"' -> Specifies the species being modeled.
# - 'longitude' and 'latitude' -> Define spatial coordinates.
# - 'variables = variable_sets' -> Environmental variables used in the model.
# - 'methods = methods' -> Uses the covariance matrix approach ('covmat').
# - 'level = 95' -> 95% confidence level for the ellipsoid.
# - 'selection_criteria = "S_OR_P"' -> Selects the best model based on suitability.
# - 'error = 5' -> Error threshold for calibration.
# - 'iterations = 100' -> Number of model iterations.
# - 'percentage = 50' -> 50% of occurrences used for testing.
# - 'overwrite = TRUE' -> Allows overwriting existing outputs.
# - 'parallel = TRUE' -> Runs computations in parallel to speed up processing.
# ------------------------------------------------------------------
calib <- ellipsoid_calibration(data = data_split, 
                               species = "species", 
                               longitude = "x", 
                               latitude = "y", 
                               variables = variable_sets, 
                               methods = methods, 
                               level = 99, 
                               selection_criteria = "S_OR_P", 
                               error = 5, 
                               iterations = 1000, 
                               percentage = 70,
                               overwrite = TRUE,
                               parallel = TRUE)

# ------------------------------------------------------------------
# 🔹 View Model Summary
# This provides an overview of the calibration results, including model performance.
# ------------------------------------------------------------------
summary(calib)

# ------------------------------------------------------------------
# 🔹 Extract the Selected Variable Set
# This retrieves the environmental variables chosen during calibration.
# ------------------------------------------------------------------
selected_vars <- variable_sets$variable_sets[calib@selected_parameters$Variable_set][1]

# Flatten the list and extract only the values
selected_vars <- unname(unlist(selected_vars))

# Display the final calibrated variable set
selected_vars

# ------------------------------------------------------------------
# 🔹 Extract Climate Data at Cholera Occurrence Locations
# ------------------------------------------------------------------
# This extracts environmental variable values from 'current_clim_scale' 
# at the locations where Cholera occurrences were recorded.
#
# - 'rast(current_clim_scale)' -> Converts raster data into a format compatible with terra.
# - 'thinned_cholera_vec' -> The vector of Cholera occurrence locations.
# - 'xy = TRUE' -> Retains spatial coordinates in the extracted dataset.
# - 'na.omit()' -> Removes rows with missing values to ensure clean data input.
# ------------------------------------------------------------------
cholera_vec_clim <- terra::extract(rast(current_clim_scale), 
                                   thinned_cholera_vec, xy = TRUE) %>%
  na.omit()

# ------------------------------------------------------------------
# 🔹 Generate Hypervolume Model Using All Environmental Variables
# ------------------------------------------------------------------
# The hypervolume model represents the ecological niche of Cholera.
#
# - 'cholera_vec_clim[, all_variables]' -> Uses all selected environmental variables.
# - 'name = "Cholera"' -> Assigns a label to the hypervolume model.
# - 'hypervolume_svm()' -> Uses Support Vector Machine (SVM) method to construct the hypervolume.
# ------------------------------------------------------------------
hv_cholera <- hypervolume_svm(cholera_vec_clim[, all_variables],
                              name = "Cholera",
                              samples.per.point = 100)

# ------------------------------------------------------------------
# 🔹 Plot the Hypervolume Model
# Enhances visualization with custom aesthetics for better interpretation.
# ------------------------------------------------------------------
plot(hv_cholera, 
     cex.legend = 1.5,    # Increases legend size
     contour.lwd = 3,     # Adjusts contour line thickness
     cex.axis = 1.5,      # Enlarges axis labels
     cex.names = 1.5,     # Enlarges variable names
     cex.centroid = 3)    # Increases centroid size

# ------------------------------------------------------------------
# 🔹 Generate Hypervolume Model Using Calibrated Variables
# ------------------------------------------------------------------
# This hypervolume model is built using only the variables selected 
# from the VIF analysis, reducing redundancy and improving efficiency.
#
# - 'cholera_vec_clim[, selected_vars]' -> Uses only the calibrated environmental variables.
# - 'name = "Cholera calibration"' -> Assigns a specific label for this model.
# ------------------------------------------------------------------
hv_cholera_calib <- hypervolume_svm(cholera_vec_clim[, selected_vars],
                                    name = "Cholera calibration",
                                    samples.per.point = 100)

# ------------------------------------------------------------------
# 🔹 Plot the Hypervolume Model with Calibrated Variables
# ------------------------------------------------------------------
plot(hv_cholera_calib,
     cex.legend = 1.5,    # Increases legend size
     contour.lwd = 3,     # Adjusts contour line thickness
     cex.axis = 1.5,      # Enlarges axis labels
     cex.names = 1.5,     # Enlarges variable names
     cex.centroid = 3)    # Increases centroid size

# Project Cholera hypervolume within the defined M region
hv_cholera_projected <- hypervolume_project(
  hv_cholera,
  rasters = raster::stack(current_clim_scale[[all_variables]]),
  type = 'inclusion',
  fast.or.accurate = 'fast'
)

# Plot projected hypervolume for M region
plot(hv_cholera_projected)
maps::map("world", add = TRUE)

# Project Cholera hypervolume with calibrated factors within the defined M region
hv_cholera_calib_projected <- hypervolume_project(
  hv_cholera_calib,
  rasters = raster::stack(current_clim_scale[[selected_vars]]),
  type = 'inclusion',
  fast.or.accurate = 'fast'
)

# Plot projected hypervolume for M region
plot(hv_cholera_calib_projected)
maps::map("world", add = TRUE)

# Load global environmental conditions
current_clim_world <- rast("Cholera/Climate/Current/current_clim.nc")

# Plot global climate data
plot(current_clim_world)

# Identify and remove problematic variables (NA values and irrelevant variables)
# Removing layers
filtered_clim_world <- current_clim_world[[-c(5, 7, 11, 17, 18, 20)]]

# Print the names of the remaining raster layers
names(filtered_clim_world)

# Scale environmental variables (standardization)
scaled_clim_world <- lapply(1:nlyr(filtered_clim_world), function(i) {
  raster::raster(scale(filtered_clim_world[[i]]))  # Scale each layer separately
})

# Convert list to RasterStack for analysis
current_clim_world_scale <- raster::stack(scaled_clim_world)

# Project Cholera hypervolume with all environmental factors globally
hv_cholera_world_projected <- hypervolume_project(
  hv_cholera,
  rasters = raster::stack(current_clim_world_scale[[all_variables]]),
  type = 'inclusion',
  fast.or.accurate = 'fast'
)

# Plot projected hypervolume for global conditions
plot(hv_cholera_world_projected)
maps::map("world", add = TRUE)

# Project Cholera hypervolume with calibrated environmental factors globally
hv_cholera_world_calib_projected <- hypervolume_project(
  hv_cholera_calib,
  rasters = raster::stack(current_clim_world_scale[[selected_vars]]),
  type = 'inclusion',
  fast.or.accurate = 'fast'
)

# Plot projected hypervolume for global conditions
plot(hv_cholera_world_calib_projected)
maps::map("world", add = TRUE)

# Load future climate projections for SSP119 (2050s)
future_ssp119 <- rast("Cholera/Climate/Future_2050_2060/future_clim_ssp119_2050s.nc")

# Remove unwanted raster layers (e.g., NA values or unrelated variables)
# Removed layers: 5 (Cloud Cover), 18 (Sea Water Speed), 16 (Sea Ice Thickness)
future_ssp119 <- future_ssp119[[-c(5, 18, 16)]]

# Load future climate projections for SSP585 (2050s)
future_ssp585 <- rast("Cholera/Climate/Future_2050_2060/future_clim_ssp585_2050s.nc")

# Apply the same filtering process for SSP585
future_ssp585 <- future_ssp585[[-c(5, 18, 16)]]

# ------------------------------------------------------------------
# 🔹 Scale environmental variables (standardization) for SSP119
# Each raster layer is standardized individually
# Scaling transforms the data to have mean = 0 and standard deviation = 1
# This prevents bias in the model due to different variable ranges
# ------------------------------------------------------------------
future_clim_ssp119_scale <- lapply(1:nlyr(future_ssp119), function(i) {
  raster::raster(scale(future_ssp119[[i]]))  # Convert SpatRaster layer to Raster after scaling
})

# Convert list of scaled layers to a RasterStack
future_clim_ssp119_stack <- raster::stack(future_clim_ssp119_scale)

# ------------------------------------------------------------------
# 🔹 Scale environmental variables (standardization) for SSP585
# Apply the same standardization process to SSP585 layers
# ------------------------------------------------------------------
future_clim_ssp585_scale <- lapply(1:nlyr(future_ssp585), function(i) {
  raster::raster(scale(future_ssp585[[i]]))  # Standardize each raster layer separately
})

# Convert list of scaled layers to a RasterStack
future_clim_ssp585_stack <- raster::stack(future_clim_ssp585_scale)

# ------------------------------------------------------------------
# 🔹 Visualizing Scaled Climate Data
# ------------------------------------------------------------------
# SSP119
plot(future_clim_ssp119_stack)

# SSP585
plot(future_clim_ssp585_stack)

# ------------------------------------------------------------------
# 🔹 Train the Ecological Niche Model for Cholera
# ------------------------------------------------------------------
# This model predicts the environmental suitability for Cholera
# using an ellipsoid niche approach based on climatic variables.
#
# - 'thinned_cholera_data' -> Cholera occurrence dataset (after thinning).
# - 'species = "species"' -> Defines the target species being modeled.
# - 'longitude' & 'latitude' -> Specify spatial coordinates of occurrences.
# - 'raster_layers' -> Current climate data for training (selected variables only).
# - 'method = "covmat"' -> Uses the covariance matrix method for modeling.
# - 'level = 95' -> Sets the confidence level at 95% for niche space.
# - 'replicates = 100' -> Runs 100 bootstrap replicates for robust estimates.
# - 'replicate_type = "bootstrap"' -> Uses bootstrapping to assess model stability.
# - 'bootstrap_percentage = 70' -> 70% of occurrence data is used for training.
# ------------------------------------------------------------------
cholera_mod <- tryCatch(ellipsoid_model(
  data = thinned_cholera_data, 
  species = "species", 
  longitude = "x",
  latitude = "y", 
  raster_layers = raster::stack(current_clim_world_scale[[selected_vars]]),
  method = "covmat", 
  level = 95, 
  replicates = 10, 
  replicate_type = "bootstrap",
  bootstrap_percentage = 70, 

  # ------------------------------------------------------------------
  # 🔹 Project the Model onto Future Climate Scenarios
  # This step projects the trained niche model onto future climate conditions.
  #
  # - 'future_ssp119' -> Future projection under SSP1-1.9 (low-emission scenario).
  # - 'future_ssp585' -> Future projection under SSP5-8.5 (high-emission scenario).
  # - 'raster::stack()' -> Ensures compatibility of projection data format.
  # ------------------------------------------------------------------
  projection_variables = list(
    future_ssp119 = raster::stack(future_clim_ssp119_stack[[selected_vars]]),
    future_ssp585 = raster::stack(future_clim_ssp585_stack[[selected_vars]])
  ),

  # ------------------------------------------------------------------
  # 🔹 Model Output Settings
  # - 'prvariables_format = "GTiff"' -> Saves the projected variables as GeoTIFF files.
  # - 'prediction = "suitability"' -> Outputs a continuous suitability map.
  # - 'return_numeric = TRUE' -> Returns numeric model output values.
  # - 'tolerance = 1e-60' -> Sets numerical tolerance to avoid computational errors.
  # - 'format = "GTiff"' -> Saves all model outputs as GeoTIFF files.
  # - 'color_palette = viridis::magma' -> Uses a color scheme for visual clarity.
  # - 'output_directory' -> Directory where model results are saved.
  # - 'overwrite = TRUE' -> Overwrites previous results if they exist.
  # ------------------------------------------------------------------
  prvariables_format = "GTiff", 
  prediction = "suitability", 
  return_numeric = TRUE,
  tolerance = 1e-60, 
  format = "GTiff", 
  color_palette = viridis::magma,
  output_directory = "Cholera/Models/ellipsenm_cholera", 
  overwrite = TRUE
))

# ------------------------------------------------------------------
# 🔹 View the Model Summary
# This provides an overview of the trained model, including performance metrics.
# ------------------------------------------------------------------
summary(cholera_mod)

# Load the projected suitability map (Current Climate)
cholera_suitability_current <- rast("Cholera/Models/ellipsenm_cholera/mean_suitability_calibration_Cholera.tif")

# Plot current suitability map
plot(cholera_suitability_current, 
     main = "Cholera Suitability - Current Climate", 
     range = c(0, 1))
maps::map("world", add = TRUE)

# Load the projected suitability map for SSP119
cholera_suitability_future_ssp119 <- rast("Cholera/Models/ellipsenm_cholera/mean_suitability_future_ssp119_Cholera.tif")

# Plot future suitability map (SSP119)
plot(cholera_suitability_future_ssp119, 
     main = "Cholera Suitability - SSP119 (2050s)", 
     range = c(0, 1))
maps::map("world", add = TRUE)

# Load the projected suitability map for SSP585
cholera_suitability_future_ssp585 <- rast("Cholera/Models/ellipsenm_cholera/mean_suitability_future_ssp585_Cholera.tif")

# Plot future suitability map (SSP585)
plot(cholera_suitability_future_ssp585, 
     main = "Cholera Suitability - SSP585 (2050s)", 
     range = c(0, 1))
maps::map("world", add = TRUE)