# =========================================================
# R version 4.3.2 (2023-10-31)  #
# Operating System: arm64    #
# Code for statistical analysis   #
# Supplement to: Gabriel Fénisse*, Manuel Chevalier, Odile Peyron, David Vincent Bekaert, Pierre-Henri Blard: Advancing Last Glacial Maximum paleoclimate reconstructions in Europe using pollen data: a multi-method (mega)biomization approach, EGUsphere [preprint], https://doi.org/10.5194/egusphere-2026-1590, 2026. #
# Contact: *Gabriel Fénisse (gabriel.fenisse@univ-lorraine.fr)[Centre de Recherches Pétrographiques et Géochimiques, Vandoeuvre-les-Nancy, 2025] #
# =========================================================

# Set up environment
# MAT transfer function reconstruction - SQ chord
# Please, install these packages if you want to process it - rioja and dplyr packages (without analogue package).
library(rioja)#1.0-6 version
library(dplyr)#1.1.4 version

# Data initialization process
# Initialization of file names
# File name for fossil pollen data
nameFile_AssF <- "./CoolingPollen.csv" # "./CoolingPollen.csv", "./AssFLDB-2.csv", "./AssFLGP.csv", "./AssFLGP_2.csv"
# File name for fossil pollen taxa identification
nameFile_Pollen <- "./Name_CoolingPollen_Harmonized.txt" # "./Name_CoolingPollen_Harmonized.txt", "./NameLDB_Harmonized.csv", "./NameLGP_Harmonized.csv", "./NameLGP2_Harmonized.csv"
# Choice 1 - Age model available (for LDB and LGP)
Ages <- read.csv2("AgeModel_LDB.txt", sep=',', dec='.') # "AgeModel_LDB.txt", "AgeModel_LGP.txt", "AgeModel_LGP2.txt"
AgeSelect <- Ages$X2
# Choice 2 - If no age model is available  
AssF <- read.csv2(nameFile_AssF, sep=";", dec=',')
AgeSelect <- seq_len(nrow(AssF))  # Sequential sample indices
# Match ages with fossil pollen data
AssF$Ages <- AgeSelect

# Variable initialization process
# Choice 1 - Biomisation selection
site <- "Cooling"   # "Cooling", "LDB", "LGP", "LGP2"
biome <- c("CLDE","TAIG","PION","CLMX","COCO","TEDE","COMX","WAMX","XERO","TUND","COST","WAST","CODE","HODE") # Without AQU and ANTH
# Automatic file generation
biomeSelect <- paste0("./", site, "/", biome, "_", site, ".csv")
biomeSelect
# Choice 2 - Megabiome selection
site <- "Cooling" # "Cooling", "LDB", "LGP", "LGP2"
# Folder names for megabiomes
folder_site <- paste0(site, "_M")
# File suffix exceptions
suffix_site <- site
# Megabiomes
biome <- c("TEFO","WTFO","BOFO","TUND","STEP")#"DESE"
# Automatic file generation
biomeSelect <- paste0("./", folder_site, "/", biome, "_", suffix_site, ".csv")
biomeSelect

# Initialisation of climate variables
VariableClim  <- rep(c(1:3),times=length(unique(biome)))
RMSEP  <-data.frame(matrix(NA,ncol=3,nrow=length(VariableClim)))
Variable  <- c("TANN","MTCO","MTWA") #Check that the climate variable names match the modern training set file

# =========================================================
# Climate reconstruction using MAT
# =========================================================
# ---------------------------------------------------------
# Load modern biome-climate calibration dataset
# ---------------------------------------------------------

OutBiomization <- read.csv2(
  "ModernClimate_Synthesis_Megabiome.csv", #"ModernClimate_Synthesis_Biome.csv", "ModernClimate_Synthesis_Megabiome.csv"
  sep = ";",
  dec = "."
)

# Remove unsampled sites
NotSampled <- OutBiomization$SampleName[
  is.na(OutBiomization$TANN)
]

OutBiomization <- subset(
  OutBiomization,
  !(Place %in% NotSampled)
)

# ---------------------------------------------------------
# Load calibration and fossil datasets (EMPD2)
# ---------------------------------------------------------

AssCal <- read.csv2(
  "EMPD2-Set.txt",
  sep = ",",
  dec = "."
)

NameAssCal <- read.csv2(
  "NamesAssCal.txt",
  sep = "\t",
  dec = "."
)

AssF <- read.csv2(
  nameFile_AssF,
  sep = ";",
  dec = ","
)

NameF <- read.csv2(
  nameFile_Pollen,
  sep = ";",
  dec = "."
)

# ---------------------------------------------------------
# Fossil pollen formatting
# ---------------------------------------------------------

rownames(AssF) <- make.names(AgeSelect, unique = TRUE)

AssF <- AssF[, -1]
AssF <- AssF * 100

colnames(AssF) <- NameF$Names

# Merge duplicated taxa
for (i in unique(NameF$Names)) {
  cols <- colnames(AssF)[colnames(AssF) == i]
  AssF[, i] <- rowSums(
    AssF[, cols, drop = FALSE]
  )
}

AssF <- AssF[, unique(NameF$Names)]

# ---------------------------------------------------------
# Modern calibration formatting
# ---------------------------------------------------------

AssCal <- subset(
  AssCal,
  AssCal$Place %in% OutBiomization$Place
)

rownames(AssCal) <- AssCal$Place
colnames(AssCal) <- NameAssCal$Names

AssCal <- AssCal[
  ,
  !colnames(AssCal) %in%
    grep("Delete", colnames(AssCal), value = TRUE)
]

# Merge duplicated taxa
for (i in names(AssF)) {
  cols <- grep(i, colnames(AssCal), value = TRUE)
  AssCal[, i] <- rowSums(
    AssCal[, cols, drop = FALSE]
  )
}

AssCal <- AssCal[, names(AssF)]

# Convert to percentages
AssCal <- sweep(
  AssCal,
  1,
  rowSums(AssCal),
  "/"
) * 100

# Save clean datasets
AssCalInt <- AssCal
AssFInt <- AssF

# =========================================================
# Function for MAT reconstruction
# =========================================================

run_MAT <- function(selected_biome, variable_index) {
  
  # Select biome
  clim <- subset(
    OutBiomization,
    Biome == selected_biome
  )
  
  modern_clim <- as.numeric(
    clim[, variable_index + 5]
  )
  
  # Select calibration spectra
  AssCal <- subset(
    AssCalInt,
    rownames(AssCalInt) %in% clim$Place
  )
  
  # Remove absent taxa
  TaxaSelect <- colnames(AssCal)[
    colSums(AssCal) == 0
  ]
  
  AssCal <- AssCal[
    ,
    !colnames(AssCal) %in% TaxaSelect
  ]
  
  AssF <- AssFInt[
    ,
    !colnames(AssFInt) %in% TaxaSelect
  ]
  
  # MAT calibration
  modMAT <- MAT(
    AssCal,
    modern_clim,
    dist.method = "sq.chord",
    lean = FALSE,
    k = 10
  )
  
  # Optimal analogs
  analogs <- which.min(
    performance(modMAT)$object[1:10, 1]
  )
  
  # Final model
  modMAT <- MAT(
    AssCal,
    modern_clim,
    dist.method = "sq.chord",
    lean = FALSE,
    k = analogs
  )
  
  # Prediction
  predMAT <- predict(
    modMAT,
    newdata = AssF,
    sse = TRUE,
    lean = FALSE,
    k = analogs,
    nboot = 200
  )
  
  # Output
  data.frame(
    Age = AgeSelect,
    Reconstruction = predMAT$fit[, 1],
    SEP = predMAT$SEP.boot[, 1]
  )
}

# =========================================================
# Main loop
# =========================================================

results <- list()

for (b in seq_along(biome)) {
  
  cat("\n")
  cat("Biome:", biome[b], "\n")
  
  reconstructions <- list()
  
  for (h in seq_along(Variable)) {
    
    cat("Variable:", Variable[h], "\n")
    
    reconstructions[[Variable[h]]] <- run_MAT(
      selected_biome = biome[b],
      variable_index = h
    )
  }
  
  climateRecon_MAT <- data.frame(
    Age = AgeSelect,
    TANN = reconstructions[[1]]$Reconstruction,
    TANN_SEP = reconstructions[[1]]$SEP,
    MTCO = reconstructions[[2]]$Reconstruction,
    MTCO_SEP = reconstructions[[2]]$SEP,
    MTWA = reconstructions[[3]]$Reconstruction,
    MTWA_SEP = reconstructions[[3]]$SEP
  )
  
  # Export
  write.csv2(
    climateRecon_MAT,
    biomeSelect[b],
    row.names = FALSE
  )
  
  results[[biome[b]]] <- climateRecon_MAT
}
# =========================================================
# END
# =========================================================
