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

library(stringr)

############################################################
# Definition of (mega)biome classes
############################################################

# List of megabiomes used in the reconstruction (unweighted approach)
biome=c("CLDE","TAIG","PION","CLMX","COCO","TEDE","COMX","WAMX","XERO","TUND","COST","WAST","HODE")
biome=c("TEFO","WTFO","BOFO","TUND","STEP") #DESE

############################################################
# 1. Reconstruction based on dominant (best) (mega)biome
#    (non-weighted approach)
############################################################

# Load biomisation output (best-guess megabiome per site)
OutBiomization <- read.csv2("OutLDB_Biomisation_Classe2.csv", # "OutFossil_MegaBiomization_Cooling.csv", "OutFossil_MegaBiomization_LDB.csv"
                            sep = ";", dec = ",")

# Reference structure used to define output format
Ex_Biome <- read.csv2("TEDE_LDB.csv", sep = ";", dec = ",")[, -1]

# Initialize output matrix
Output_Recon <- data.frame(matrix(NA,
                                  ncol = ncol(Ex_Biome),
                                  nrow = nrow(Ex_Biome)))

colnames(Output_Recon) <- colnames(Ex_Biome)
rownames(Output_Recon) <- rownames(Ex_Biome)

# Loop over sites
for (n in 1:nrow(OutBiomization)) {
  
  # Selected dominant biome for each site
  Best <- OutBiomization$Best[n]
  
  # Corresponding biome file
  Nom_File <- str_c(Best, "_LDB.csv")
  
  # Extract reconstructed values
  Output_Recon[n, ] <- read.csv2(Nom_File, sep = ";", dec = ",")[n, -1]
}

# Save non-weighted reconstruction
write.csv2(Output_Recon,
           "Output_Recon_CREST_Biome_LDB.csv")

############################################################
# 2. Weighted (mega)biome reconstruction
############################################################

# Load biomisation probabilities / scores
OutBiomization <- read.csv2("OutFossil_MegaBiomization_Cooling.csv",
                            sep = ";", dec = ",")


# Remove auxiliary columns not related to biome scores
OutBiomization <- OutBiomization[,  -c(1,14,16,17,18)] # -c(1,14,16,17,18) for biomes, -c(1,5,9,10) for megabiomes

# Reference structure
Ex_Biome <- read.csv2("TEDE_LDB.csv", sep = ";", dec = ",")[, -1]

# Initialize output
Output_Recon <- data.frame(matrix(NA,
                                  ncol = ncol(Ex_Biome),
                                  nrow = nrow(Ex_Biome)))

colnames(Output_Recon) <- colnames(Ex_Biome)
rownames(Output_Recon) <- rownames(Ex_Biome)

# Loop over sites
for (n in 1:nrow(OutBiomization)) {
  
  # Temporary storage for biome-weighted contributions
  data <- data.frame(matrix(NA,
                            ncol = ncol(Ex_Biome),
                            nrow = length(biome)))
  
  # Loop over biome classes
  for (m in 1:length(biome)) {
    
    Nom_File <- str_c(biome[m], "_LDB.csv")
    
    # Extract biome-specific reconstruction
    data[m, ] <- read.csv2(Nom_File, sep = ";", dec = ",")[n, -1]
    
    # Apply biomisation weight
    data[m, ] <- data[m, ] * OutBiomization[n, biome[m]]
  }
  
  # Normalized weighted mean across biomes
  Output_Recon[n, ] <- colSums(data) / sum(OutBiomization[n, ])
}

# Save weighted reconstruction
write.csv2(Output_Recon,
           "Output_Recon_CREST_Biome_Weighted_LDB.csv")

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