# =========================================================
# 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 - SQ chord - Manual CV
library(crestr)# 1.3.1.9
library(rio)# 1.1.4 version

#Biomisation process - Load input data 
biome=c("CLDE","TAIG","PION","CLMX","COCO","TEDE","COMX","WAMX","XERO","TUND","COST","WAST","HODE","CODE")
biomeSelect=c("./Cooling/CLDE_Cooling.csv","./Cooling/TAIG_Cooling.csv","./Cooling/PION_Cooling.csv","./Cooling/CLMX_Cooling.csv","./Cooling/COCO_Cooling.csv","./Cooling/TEDE_Cooling.csv","./Cooling/COMX_Cooling.csv","./Cooling/WAMX_Cooling.csv","./Cooling/XERO_Cooling.csv","./Cooling/TUND_Cooling.csv","./Cooling/COST_Cooling.csv","./Cooling/WAST_Cooling.csv","./Cooling/HODE_Cooling.csv","./Cooling/CODE_Cooling.csv")

#Megabiomisation process - Load input data 
biome=c("TEFO","WTFO","BOFO","TUND","STEP")#"DESE"
biomeSelect=c("./Cooling_M/TEFO_Cooling.csv","./Cooling_M/WTFO_Cooling.csv","./Cooling_M/BOFO_Cooling.csv","./Cooling_M/TUND_Cooling.csv","./Cooling_M/STEP_Cooling.csv","./Cooling_M/DESE_Cooling.csv")

#File names
nameFile_AssF="./CoolingPollen_copie.csv"
nameFile_Pollen="./Name_CoolingPollen_Harmonized.txt"
AgeSelect=seq(1,nrow(AssF))# List

# =====================================================
# 0. LOAD DATA
# =====================================================

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

#To change
metadata <- metadata[,!colnames(metadata) %in% c("TANN", "MTWA")]

df <- read.csv2(
  "EMPD2-Set-Pollen_Normalized.txt",
  sep = "\t",
  dec = ".",
  stringsAsFactors = FALSE
)
colnames(df)=NameAssCal$Names
NotSampled=metadata$Place[is.na(metadata$TANN)]
metadata=subset(metadata,!(metadata$Place %in% NotSampled))
df=subset(df,(df$Place %in% metadata$Place))
Place=df$Place

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

colnames(df) <- NameAssCal$Names

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

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

# =====================================================
# 1. FOSSIL PREPROCESSING
# =====================================================

AgeSelect=AssF$SiteName
rownames(AssF)=make.names(AgeSelect,unique=TRUE)#Age
AssF=AssF[,-c(1)]
AssF=AssF*100
colnames(AssF)=NameF$Names

AssF=AssF[,which(colnames(AssF) %in% c(NameF$Names))]
for (i in unique(NameF$Names))
{AssF[,i]=rowSums(AssF[grep(i,colnames(AssF),value=TRUE)])}
AssF=AssF[,which(colnames(AssF) %in% c(unique(NameF$Names)))]

TaxaSelect=colnames(AssF)[colSums(AssF)==0]
AssF=AssF[,!colnames(AssF) %in% c(TaxaSelect)]
NameF=TaxaSelect

for (i in names(AssF))
{df[,i]=rowSums(df[grep(i,colnames(df),value=TRUE)])}
df=df[,which(colnames(df) %in% c(names(AssF)))]
df=df[,!colnames(df) %in% grep("Delete",colnames(df),value=TRUE)]

AssF$Ages=AgeSelect
df$Place=Place

# =====================================================
# 2. BIOME LOOP STORAGE
# =====================================================

biome_list <- unique(metadata$Biome)

pred_matrix <- NULL

all_obs <- NULL

pred_matrix <- matrix(
  NA,
  nrow = nrow(metadata),   
  ncol = length(biome)
)

colnames(pred_matrix) <- biome
rownames(pred_matrix) <- metadata$Place

# =====================================================
# 3. LOOP OVER BIOMES (CREST CALIBRATION)
# =====================================================

for (b in seq_along(biome_list)) {
  
  cat("\nBiome:", biome_list[b], "\n")
  
  # =====================================================
  # CREST SET MODERN DATA
  # =====================================================
  
  full_df <- merge(metadata, df, on=c('place'), all=TRUE)
  metadata_id    <- c(1:2,3:4)
  coordinates_id <- 2:3
  climate_id     <- c(6)
  climate_names  <- colnames(full_df)[climate_id]
  taxa_names     <- colnames(full_df)[7:length(full_df)]
  
  print("Attention BIOME")
  full_df <- subset(full_df,(full_df$Biome == biome[b]))
  print(biome[b])
  
  site_idx <- full_df$Place[full_df$Biome == biome_list[b]]
  
  distributions <- full_df[FALSE, c(1,2,coordinates_id,climate_id,2)]
  for(tax in taxa_names) {
    w <- which(full_df[, tax] > 0)
    if (tax  %in% names(AssF)&((length(w)!=0)))
    {{distributions <- rbind(distributions,
                             cbind( 'taxonid' = rep(tax, length(w)),
                                    'ProxyName' = rep(tax, length(w)),
                                    full_df[w, coordinates_id],
                                    'MTCO'=full_df[w, climate_id],
                                    'weight'=full_df[w, tax]
                             )
    )
    }}}  
  
  colnames(distributions)[1:4]=c('taxonid', 'ProxyName','longitude', 'latitude')
  head(distributions)
  
  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))]
  
  climate_space <- full_df[, c(coordinates_id, climate_id)]
  colnames(climate_space)[1:2]=c('longitude', 'latitude')
  head(climate_space)
  
  full_df[,1:6] <- NULL
  full_df$Ages <- 1:nrow(full_df)
  full_df=full_df[,c("Ages",unique(distributions$ProxyName))]
  
  # =====================================================
  # CALIBRATION
  # =====================================================
  
  rcnstrctn <- crest.set_modern_data(distributions,weight=TRUE,minGridCells=10,df=full_df,climate_space=climate_space,climate=climate_names)
  
  rcnstrctn <- crest.calibrate(
    rcnstrctn,
    geoWeighting = TRUE,
    climateSpaceWeighting = TRUE,
    shape=c('normal'),
    bin_width =c(2),
    npoints=2000,
    verbose=TRUE)
  
  # =====================================================
  # RECONSTRUCTION (ALL SAMPLES)
  # =====================================================
  
  #Taxa select before rcnstrct
  TaxaSelect=data.frame(matrix(NA,ncol=1,nrow=length(rcnstrctn$inputs$taxa.name)))
  for (pol in 1:length(rcnstrctn$inputs$taxa.name))
  {if((anyNA(list(rcnstrctn$modelling$pdfs[pol]),recursive=TRUE)==TRUE))
  {TaxaSelect[pol,1]=rcnstrctn$inputs$taxa.name[pol]}
    if (nrow(unique(data.frame(rcnstrctn$modelling$distributions[pol])))<5)
    {print(rcnstrctn$inputs$taxa.name[pol])
      TaxaSelect[pol,1]=rcnstrctn$inputs$taxa.name[pol]}}
  TaxaSelect=na.omit(TaxaSelect)
  rcnstrctn <- excludeTaxa(rcnstrctn,taxa=TaxaSelect[,1], climate=rcnstrctn$parameters$climate)
  
  # Prediction
  pred <- crest.reconstruct(rcnstrctn,presenceThreshold = 0,taxWeight = "normalisation",uncertainties = c(0.5,0.95),verbose=TRUE)
  
  pred_matrix[site_idx,biome[b]] <- pred$reconstruction$MTCO$optima$mean
}

# =====================================================
# 4. MULTI-BIOME WEIGHTED FUSION
# =====================================================

final_predictions <- rep(
  NA,
  nrow(pred_matrix)
)

biome_scores <- read.csv2(
  "OutModern_Megabiomization_Cooling.csv", #"OutModern_Biomization_Cooling.csv", "OutModern_Megabiomization_Cooling.csv"
  sep = ";",
  dec = ".",
  header = TRUE,
  stringsAsFactors = FALSE
)

for (i in 1:nrow(pred_matrix)) {
  
  weights <- as.numeric(
    biome_scores[i, biome]
  )
  
  weights <- weights / sum(weights, na.rm = TRUE)
  
  final_predictions[i] <-
    weighted.mean(
      pred_matrix[i, ],
      weights,
      na.rm = TRUE
    )
}

# =====================================================
# 5. VALIDATION
# =====================================================

valid_idx <- which(
  !is.na(final_predictions)
)

allClim_select=metadata$MTCO
correlation <- cor(
  final_predictions[valid_idx],
  allClim_select[valid_idx]
)

rmse <- sqrt(mean(
  (final_predictions[valid_idx] -
     allClim_select[valid_idx])^2
))

plot(final_predictions[valid_idx],
     allClim_select[valid_idx])
print(correlation)
print(rmse)

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