# =========================================================
# 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-CREST
##R_Package
library(crestr) # 1.3.1.9
library(rio) # 1.1.0

# 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

# =========================================================
# CREST reconstruction - optimized version
# =========================================================

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

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

df <- read.csv2(
  "EMPD2-Set-Pollen_Normalized.txt",
  sep = "\t",
  dec = "."
)
df <- df[, colnames(df) != "ID"]

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

colnames(df) <- NameAssCal$Names

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

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

# ---------------------------
# Preprocess fossil 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)]

# remove empty taxa
TaxaSelect <- colnames(AssF)[colSums(AssF) == 0]
AssF <- AssF[, !colnames(AssF) %in% TaxaSelect]

# ---------------------------
# Loop by biome
# ---------------------------

for (b in seq_along(biome)) {
  
  cat("\nBiome:", biome[b], "\n")
  
  # ---------------------------
  # Filter metadata ONCE per biome
  # ---------------------------
  
  meta_b <- subset(metadata, Biome == biome[b])
  
  NotSampled <- meta_b$Place[is.na(meta_b$TANN)]
  meta_b <- subset(meta_b, !(Place %in% NotSampled))
  
  df_b <- subset(df, Place %in% meta_b$Place)
  
  # remove taxa not in fossil
  df_b <- df_b[, c("Place", intersect(names(AssF), colnames(df_b)))]
  
  # ---------------------------
  # Build climate space
  # ---------------------------
  
  coordinates_id <- 2:3
  climate_id <- c(6, 7, 8)
  
  climate_names <- colnames(meta_b)[climate_id]
  
  climate_space <- meta_b[, c(coordinates_id, climate_id)]
  colnames(climate_space)[1:2] <- c("longitude", "latitude")
  
  # ---------------------------
  # Build distributions WITHOUT rbind loop
  # ---------------------------
  
  taxa_names <- intersect(names(AssF), colnames(df_b))
  
  distributions_list <- lapply(taxa_names, function(tax) {
    
    w <- which(df_b[, tax] > 0)
    
    if (length(w) == 0) return(NULL)
    
    data.frame(
      taxonid = tax,
      ProxyName = tax,
      longitude = meta_b[w, coordinates_id[1]],
      latitude = meta_b[w, coordinates_id[2]],
      climate_space[w, climate_id-3],
      weight = df_b[w, tax]
    )
  })
  
  distributions <- do.call(rbind, distributions_list)
  
  # ---------------------------
  # Reduce fossil dataset
  # ---------------------------
  
  AssF_b <- AssF[, intersect(colnames(AssF), distributions$ProxyName)]
  AssF_b <- na.omit(AssF_b)
  
  AssF_b$Ages <- AgeSelect
  AssF_b=AssF_b[,c('Ages',unique(distributions$ProxyName))]
  
  # ---------------------------
  # CREST model
  # ---------------------------
  
  rcnstrctn <- crest.set_modern_data(
    distributions,
    weight = TRUE,
    minGridCells = 20,
    df = AssF_b,
    climate_space = climate_space,
    climate = climate_names
  )
  
  rcnstrctn <- crest.calibrate(
    rcnstrctn,
    geoWeighting = TRUE,
    climateSpaceWeighting = TRUE,
    shape = rep("normal", length(climate_names)),
    bin_width = rep(2, length(climate_names)),
    npoints = 2000,
    verbose = TRUE
  )
  
  # ---------------------------
  # Remove bad taxa
  # ---------------------------
  
  TaxaSelect <- data.frame()
  
  for (pol in seq_along(rcnstrctn$inputs$taxa.name)) {
    
    if (
      anyNA(rcnstrctn$modelling$pdfs[[pol]]) ||
      nrow(unique(rcnstrctn$modelling$distributions[[pol]])) < 5
    ) {
      TaxaSelect <- rbind(
        TaxaSelect,
        rcnstrctn$inputs$taxa.name[pol]
      )
    }
  }
  
  if (nrow(TaxaSelect) > 0) {
    rcnstrctn <- excludeTaxa(
      rcnstrctn,
      taxa = TaxaSelect[, 1],
      climate = rcnstrctn$parameters$climate
    )
  }
  
  # ---------------------------
  # Reconstruction
  # ---------------------------
  
  rcnstrctn <- crest.reconstruct(
    rcnstrctn,
    presenceThreshold = 0,
    taxWeight = "normalisation",
    uncertainties = c(0.5, 0.95),
    verbose = TRUE
  )
  
  # ---------------------------
  # Export
  # ---------------------------
  
  data <- data.frame(
    TANN = rcnstrctn$reconstructions$TANN$optima,
    TANN_u = rcnstrctn$reconstructions$TANN$uncertainties,
    MTCO = rcnstrctn$reconstructions$MTCO$optima,
    MTCO_u = rcnstrctn$reconstructions$MTCO$uncertainties,
    MTWA = rcnstrctn$reconstructions$MTWA$optima,
    MTWA_u = rcnstrctn$reconstructions$MTWA$uncertainties
  )
  
  write.csv2(data, biomeSelect[b], row.names = FALSE)
}

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

