# =========================================================
# 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] #
# =========================================================

#Biomisation-WAPLS
##R_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"  # "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

# =========================================================
# WAPLS climate reconstruction (optimized version)
# =========================================================

# ---------------------------
# Load static data ONCE
# ---------------------------

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


NotSampled <- OutBiomization$SampleName[
  is.na(OutBiomization$TANN)
]

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

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 preprocessing ONCE
# ---------------------------

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

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

colnames(AssF) <- NameF$Names

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

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

AssFInt <- AssF
Age <- AgeSelect

# ---------------------------
# Calibration preprocessing ONCE
# ---------------------------

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 taxa (modern)
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

AssCalInt <- AssCal

# ---------------------------
# Output storage
# ---------------------------

m <- 1
results <- list()

# =========================================================
# MAIN LOOP
# =========================================================

for (b in seq_along(biome)) {
  
  for (h in seq_along(Variable)) {
    
    # ---------------------------
    # Extract climate variable
    # ---------------------------
    
    clim <- subset(
      OutBiomization,
      Biome == biome[b]
    )
    
    modern_clim <- as.numeric(clim[, h + 5])
    
    # ---------------------------
    # Subset calibration
    # ---------------------------
    
    AssCal_b <- AssCalInt[
      rownames(AssCalInt) %in% clim$Place,
      ,
      drop = FALSE
    ]
    
    # Remove empty taxa
    TaxaSelect <- colnames(AssCal_b)[
      colSums(AssCal_b) == 0
    ]
    
    AssCal_b <- AssCal_b[
      ,
      !colnames(AssCal_b) %in% TaxaSelect
    ]
    
    AssF_b <- AssFInt[
      ,
      !colnames(AssFInt) %in% TaxaSelect
    ]
    
    # ---------------------------
    # WAPLS model
    # ---------------------------
    
    fit <- WAPLS(
      sqrt(AssCal_b),
      modern_clim,
      npls = 2
    )
    
    cross <- crossval(fit)
    TT <- rand.t.test(cross)
    
    nb_comp <- 1
    
    if (TT[2, 6] <= -5 & TT[2, 7] <= 0.05) {
      nb_comp <- 2
    }
    
    cat("Biome:", biome[b], "\n")
    cat("Variable:", Variable[h], "\n")
    cat("Components:", nb_comp, "\n")
    
    RMSEP[m, ] <- TT[nb_comp, ]
    m <- m + 1
    
    pred <- predict(
      fit,
      sqrt(AssF_b),
      npls = nb_comp,
      sse = TRUE
    )
    
    # ---------------------------
    # Store results
    # ---------------------------
    
    results[[paste(b, h, sep = "_")]] <- data.frame(
      Age = Age,
      value = pred$fit[, 1],
      error = pred$SEP.boot[, 1]
      #biome = biome[b],
      #variable = Variable[h]
    )
    
    # ---------------------------
    # Export per variable
    # ---------------------------
    
    if (h == 1) ReconTANN <- results[[paste(b, h, sep = "_")]]
    if (h == 2) ReconMTCO <- results[[paste(b, h, sep = "_")]]
    if (h == 3) ReconMTWA <- results[[paste(b, h, sep = "_")]]
    
    if (h == 3) {
      
      climateRecon_WAPLS <- data.frame(
        ReconTANN,
        ReconMTCO,
        ReconMTWA
      )
      
      write.csv2(
        climateRecon_WAPLS,
        biomeSelect[b],
        row.names = FALSE
      )
    }
  }
}

# =========================================================
# END
# =========================================================

