# This script is used to analyse the C profiles generated by the Differential Evolution algorithm

rm(list=ls(all = TRUE))
cat("\014")

library("ggplot2")
library("reshape2")
library("gridExtra")
library("readxl")
library("corrplot")
library("RColorBrewer")
library("gridGraphics")
library("grid")
library("dplyr")
library("scales")
library("ggpubr")
library("GGally")
library("patchwork")

source("f_createClasses.R")
source("f_load_measured_data.R")

# ---------------------------
# The measured data is loaded
# ---------------------------

createInputClass()
par <- new("inputObject")
par@evalDepth <- .7
out <- loadMeasuredData("Hainich")
measuredData_soil <- out$measuredData_soil
rm(out)

# ------------------------------------------
# The locations of the data files to be used
# -----------------------------------------

# These files contain the model output for *all* parameter combinations
# that have been tested by the DE algorithm
dir_C <- "Analysis DE output/outputForSeries_C.RData"
dir_C_d13C <- "Analysis DE output/outputForSeries_C_d13C.RData"
dir_C_d14C <- "Analysis DE output/outputForSeries_C_d14C.RData"
dir_C_d13C_d14C <- "Analysis DE output/outputForSeries_C_d13C_d14C.RData"

# ------------------------------------------------------------------------------
# Defining the fraction of best runs to be retained
# ------------------------------------------------------------------------------

Qthreshold <- 0.10

# ------------------------------------------------------------------------------
# C only calibration
# ------------------------------------------------------------------------------

# The data is loaded
load(dir_C)

C <- as.data.frame(t(outputForSeries$Ctot))
d13C <- as.data.frame(t(outputForSeries$d13C))
d14C <- as.data.frame(t(outputForSeries$d14C))
depth <- outputForSeries$midDepth
error <- as.data.frame(outputForSeries$error)
names(error) <- "error"

relError_OC_POC <- as.data.frame(outputForSeries$relError_OC_POC)
names(relError_OC_POC) <- "relError_OC_POC"

relError_OC_MAOC <- as.data.frame(outputForSeries$relError_OC_MAOC)
names(relError_OC_MAOC) <- "relError_OC_MAOC"

relError_d13C_POC <- as.data.frame(outputForSeries$relError_d13C_POC)
names(relError_d13C_POC) <- "relError_d13C_POC"

relError_d13C_MAOC <- as.data.frame(outputForSeries$relError_d13C_MAOC)
names(relError_d13C_MAOC) <- "relError_d13C_MAOC"

relError_d14C_POC <- as.data.frame(outputForSeries$relError_d14C_POC)
names(relError_d14C_POC) <- "relError_d14C_POC"

relError_d14C_MAOC <- as.data.frame(outputForSeries$relError_d14C_MAOC)
names(relError_d14C_MAOC) <- "relError_d14C_MAOC"

POC_perc <- as.data.frame(t(outputForSeries$POC_perc))
bulkOC_perc <- as.data.frame(t(outputForSeries$bulkOC_perc))

# The distribution of the errors is plotted
e1 <- ggplot(error, aes(error)) +
  geom_histogram(binwidth = 0.1) +
  scale_x_continuous(limits = c(0,50)) +
  scale_y_continuous(limits = c(0,1000))

# Results with a high error are removed
# errorTol_C <- 0.05
# n <- which(error > errorTol_C)

# The X % lowest errors are retained
q_C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
n <- which(error > q_C) # The rows with errors larger than this quantile

C <- subset(C,select = -n)
d13C <- subset(d13C,select = -n)
d14C <- subset(d14C,select = -n)
error_opt <- error[-n,1]

# The errors for the different evaluated criteria
error_OC_POC_opt <- relError_OC_POC[-n,1]
error_OC_MAOC_opt <- relError_OC_MAOC[-n,1]
error_d13C_POC_opt <- relError_d13C_POC[-n,1]
error_d13C_MAOC_opt <- relError_d13C_MAOC[-n,1]
error_d14C_POC_opt <- relError_d14C_POC[-n,1]
error_d14C_MAOC_opt <- relError_d14C_MAOC[-n,1]

POC_perc <- subset(POC_perc,select = -n)
bulkOC_perc <- subset(bulkOC_perc,select = -n)

# # The total amounts of POC and bulkOC are calcultaed
POC_perc_tot <- unname(colSums(POC_perc))
bulkOC_perc_tot <- unname(colSums(bulkOC_perc))

# # The average errors for OC, d13C and d14C
avgError_OC_Conly <- mean(error_OC_POC_opt) + mean(error_OC_MAOC_opt)
avgError_d13C_Conly <- mean(error_d13C_POC_opt) + mean(error_d13C_MAOC_opt)
avgError_d14C_Conly <- mean(error_d14C_POC_opt) + mean(error_d14C_MAOC_opt)
totalError_Conly <- sum(avgError_OC_Conly + avgError_d13C_Conly + avgError_d14C_Conly)

# The errors per run are weighted based on the simulated amounts of POC and bulkOC
indErrors_OC_Conly <- ((POC_perc_tot * error_OC_POC_opt) + (bulkOC_perc_tot * error_OC_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d13C_Conly <- ((POC_perc_tot * error_d13C_POC_opt) + (bulkOC_perc_tot * error_d13C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d14C_Conly <- ((POC_perc_tot * error_d14C_POC_opt) + (bulkOC_perc_tot * error_d14C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)

# The average errors are calculated
avgError_OC_Conly <- mean(indErrors_OC_Conly)
avgError_d13C_Conly <- mean(indErrors_d13C_Conly)
avgError_d14C_Conly <- mean(indErrors_d14C_Conly)
totalError_Conly <- avgError_OC_Conly + avgError_d13C_Conly + avgError_d14C_Conly

# Q1 and Q3 are calculated
Q.OC.Conly <- c(summary(indErrors_OC_Conly)[["1st Qu."]], summary(indErrors_OC_Conly)[["3rd Qu."]])
Q.d13C.Conly <- c(summary(indErrors_d13C_Conly)[["1st Qu."]], summary(indErrors_d13C_Conly)[["3rd Qu."]])
Q.d14C.Conly <- c(summary(indErrors_d14C_Conly)[["1st Qu."]], summary(indErrors_d14C_Conly)[["3rd Qu."]])

# The distribution of the errors is plotted
# ggplot() +
#   geom_histogram(data = error, aes(error), binwidth = 0.005) +
#   geom_histogram(data = as.data.frame(error[-n,]), aes(error[-n,]), binwidth = 0.005, color = "green")
# 
# e1 <- ggplot(as.data.frame(error[-n,]), aes(error), color = "blue") +
#   geom_histogram(binwidth = 0.005)

# Unrealistic d13C and d14C simulations are removed
# r <- which(d13C > -24, arr.ind = T)
# a <- unique(r[,2])
# d13C[,a] <- NA

# r <- which(d14C > 200, arr.ind = T)
# a <- unique(r[,2])
# d14C[,a] <- NA

# The average C, d13C and d14C depth profiles are calculated
df_avg <- data.frame(depth = depth,
                     avg_C = rowMeans(C,na.rm = T),
                     avg_d13C = rowMeans(d13C,na.rm = T),
                     avg_d14C = rowMeans(d14C,na.rm = T))

# The optimal profiles are stored
rowNum_opt <- which(error_opt == min(error_opt))
df_opt <- data.frame("depth" = depth,
                     "opt_C" = C[,rowNum_opt],
                     "opt_d13C" = d13C[,rowNum_opt],
                     "opt_d14C" = d14C[,rowNum_opt])

# The depth is added to the data frames
C$depth = depth
d13C$depth = depth
d14C$depth = depth

# The data frames are reformatted for plotting
df_C <- melt(C, id.vars = "depth")
df_d13C <- melt(d13C, id.vars = "depth")
df_d14C <- melt(d14C, id.vars = "depth")

# A function to plot the C depth profiles
plot_C <- function(df_C, df_avg, df_opt, measuredData_soil, color_lines, color_avg, par){

  p <- ggplot(df_C, aes(x = value, y = depth, colour = variable)) +
    geom_path() +
    scale_y_continuous(expand = c(0,0), trans = "reverse", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_x_continuous(expand = c(0,0), position = "top", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_colour_manual(labels = unique(df_C$variable),
                        values = rep(color_lines,length(unique(df_C$variable)))) +
    coord_cartesian(xlim = c(0, 10), ylim = c(df_avg$depth[9], 0)) +
    geom_path(data = df_avg, aes(x = avg_C, y = depth), size = 2, colour = color_avg) +
    # geom_path(data = df_opt, aes(x = opt_C, y = depth), size = 1, colour = "black", linetype = "dotted") +
    geom_point(data = measuredData_soil, aes(x = fPOC_conc + MAOC_conc, y = Depth), color = "black") +
    
    geom_errorbarh(data = measuredData_soil, aes(y = Depth, xmin = (fPOC_conc +  MAOC_conc) - sqrt(fPOC_conc_stDev^2 + MAOC_conc_stDev^2), xmax = (fPOC_conc +  MAOC_conc) + sqrt(fPOC_conc_stDev^2 + MAOC_conc_stDev^2), height = .01), color = "black", inherit.aes = FALSE) +
    
    theme_classic() +
    theme(legend.position = "none",
          panel.border = element_rect(color = "black",
                                      fill = NA,
                                      size = 1),
          # axis.line = element_line(colour = "black", 
          #                          size = 0),
          axis.text=element_text(size = 16),
          axis.title=element_text(size = 16),
          plot.margin = margin(0, 5, 0, 0, "mm")) + # #top, right, bottom, left
    labs(x = "C (%)")
  
  return(p)
  
}

# A function to plot the d13C depth profiles
plot_d13C <- function(df_d13C, df_avg, df_opt, measuredData_soil, color_lines, color_avg, par){
  
  p <- ggplot(df_d13C, aes(x = value, y = depth, colour = variable)) + 
    geom_path() +
    scale_y_continuous(expand = c(0,0), trans = "reverse", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_x_continuous(expand = c(0,0), position = "top", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_colour_manual(labels = unique(df_C$variable),
                        values = rep(color_lines,length(unique(df_C$variable)))) +
    coord_cartesian(xlim = c(-28, -24), ylim = c(df_avg$depth[9], 0)) +
    geom_path(data = df_avg, aes(x = avg_d13C, y = depth), size = 2, colour = color_avg) +
    geom_point(data = measuredData_soil, aes(x = total_d13C, y = Depth), color = "black") +
    # geom_path(data = df_opt, aes(x = opt_d13C, y = depth), size = 2, colour = "black", linetype = "dotted") +
    
    geom_errorbarh(data = measuredData_soil, aes(y = Depth, xmin = ((fPOC_d13C*fPOC_conc +  MAOC_d13C*MAOC_conc)/(fPOC_conc +  MAOC_conc)) - sqrt(fPOC_d13C_stDev^2 + MAOC_d13C_stDev^2), xmax = ((fPOC_d13C*fPOC_conc +  MAOC_d13C*MAOC_conc)/(fPOC_conc +  MAOC_conc)) + sqrt(fPOC_d13C_stDev^2 + MAOC_d13C_stDev^2), height = .01), color = "black", inherit.aes = FALSE) +
    
    theme_classic() +
    theme(legend.position = "none",
          panel.border = element_rect(color = "black",
                                      fill = NA,
                                      size = 1),
          # axis.line = element_line(colour = "black", 
          #                          size = 0),
          axis.text=element_text(size = 16),
          axis.title=element_text(size = 16),
          plot.margin = margin(0, 5, 0, 0, "mm")) + # #top, right, bottom, left
    labs(x = expression(paste(delta^{13}, "C (\u2030)")))
  
  return(p)
  
}

# A function to plot the d14C depth profiles
plot_d14C <- function(df_d14C, df_avg, df_opt, measuredData_soil, color_lines, color_avg, par){
  
  p <- ggplot(df_d14C, aes(x = value, y = depth, colour = variable)) + 
    geom_path() +
    scale_y_continuous(expand = c(0,0), trans = "reverse", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_x_continuous(expand = c(0,0), position = "top", sec.axis = sec_axis(~ ., labels = NULL)) +
    coord_cartesian(xlim = c(-500, 200), ylim = c(df_avg$depth[9], 0)) +
    scale_colour_manual(labels = unique(df_C$variable),
                        values = rep(color_lines,length(unique(df_C$variable)))) +
    geom_path(data = df_avg, aes(x = avg_d14C, y = depth), size = 2, colour = color_avg) +
    geom_point(data = measuredData_soil, aes(x = total_d14C, y = Depth), color = "black") +
    # geom_path(data = df_opt, aes(x = opt_d14C, y = depth), size = 2, colour = "black", linetype = "dotted") +
    
    geom_errorbarh(data = measuredData_soil, aes(y = Depth, xmin = ((fPOC_d14C*fPOC_conc +  MAOC_d14C*MAOC_conc)/(fPOC_conc +  MAOC_conc)) - sqrt(fPOC_d14C_stDev^2 + MAOC_d14C_stDev^2), xmax = ((fPOC_d14C*fPOC_conc +  MAOC_d14C*MAOC_conc)/(fPOC_conc +  MAOC_conc)) + sqrt(fPOC_d14C_stDev^2 + MAOC_d14C_stDev^2), height = .01), color = "black", inherit.aes = FALSE) +
    
    theme_classic() +
    theme(legend.position = "none",
          panel.border = element_rect(color = "black",
                                      fill = NA,
                                      size = 1),
          # axis.line = element_line(colour = "black", 
          #                          size = 0),
          axis.text=element_text(size = 16),
          axis.title=element_text(size = 16),
          plot.margin = margin(0, 5, 0, 0, "mm")) + # #top, right, bottom, left
    labs(x = expression(paste(Delta^{14}, "C (\u2030)")))
  
  return(p)
  
}

# --------
# Plotting
# --------

# A function to round the error to 2 number after the comma
round2dec <- function(num,dec){
  return(round(num*10^dec)/(10^dec))
}

p1 <- plot_C(df_C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p2 <- plot_d13C(df_d13C, df_avg, df_opt, measuredData_soil, color_lines = "grey75", color_avg = "grey40", par)
p3 <- plot_d14C(df_d14C, df_avg, df_opt, measuredData_soil, color_lines = "grey75", color_avg = "grey40", par)

p1 <- p1 + annotate("label", x = 9.7, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic') +
  annotate("text", x = 9.7, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_OC_Conly,3), " [", round2dec(Q.OC.Conly[1],3), ";", round2dec(Q.OC.Conly[2],3),"]"))

p2 <- p2 + annotate("label", x = -27.9, y = 0.57, label = "Calibrated for C", hjust = 0, size = 4, fontface = 'italic') +
  annotate("text", x = -27.9, y = 0.5, hjust = 0, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d13C_Conly,3), " [", round2dec(Q.d13C.Conly[1],3), ";", round2dec(Q.d13C.Conly[2],3),"]"))

p3 <- p3 + annotate("label", x = 180, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic') +
  annotate("text", x = 180, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d14C_Conly,3), " [", round2dec(Q.d14C.Conly[1],3), ";", round2dec(Q.d14C.Conly[2],3),"]"))


# p2 <- p2 + annotate("label", x = -27.9, y = 0.57, label = paste("Calibrated for C; ", "\u03F5\u0304", " = ", round2dec(avgError_d13C_Conly,3)), hjust = 0, size = 4, fontface = 'italic')
# p3 <- p3 + annotate("label", x = -170, y = 0.57, label = paste("Calibrated for C; ", "\u03F5\u0304", " = ", round2dec(avgError_d14C_Conly,3)), hjust = 0, size = 4, fontface = 'italic')

rm(list = c("df_C", "df_d13C", "df_d14C", "df_avg", "df_opt"))

# grid.arrange(p1,p2,p3, nrow = 1)

# ------------------------------------------------------------------------------
# C and d13C
# ------------------------------------------------------------------------------

# The data is loaded
load(dir_C_d13C)

C <- as.data.frame(t(outputForSeries$Ctot))
d13C <- as.data.frame(t(outputForSeries$d13C))
d14C <- as.data.frame(t(outputForSeries$d14C))
depth <- outputForSeries$midDepth
error <- as.data.frame(outputForSeries$error)
names(error) <- "error"

relError_OC_POC <- as.data.frame(outputForSeries$relError_OC_POC)
names(relError_OC_POC) <- "relError_OC_POC"

relError_OC_MAOC <- as.data.frame(outputForSeries$relError_OC_MAOC)
names(relError_OC_MAOC) <- "relError_OC_MAOC"

relError_d13C_POC <- as.data.frame(outputForSeries$relError_d13C_POC)
names(relError_d13C_POC) <- "relError_d13C_POC"

relError_d13C_MAOC <- as.data.frame(outputForSeries$relError_d13C_MAOC)
names(relError_d13C_MAOC) <- "relError_d13C_MAOC"

relError_d14C_POC <- as.data.frame(outputForSeries$relError_d14C_POC)
names(relError_d14C_POC) <- "relError_d14C_POC"

relError_d14C_MAOC <- as.data.frame(outputForSeries$relError_d14C_MAOC)
names(relError_d14C_MAOC) <- "relError_d14C_MAOC"

POC_perc <- as.data.frame(t(outputForSeries$POC_perc))
bulkOC_perc <- as.data.frame(t(outputForSeries$bulkOC_perc))

# # Results with a high error are removed
# errorTol_C_d13C <- 0.4
# n <- which(error > errorTol_C_d13C)

# The X % lowest errors are retained
q_C_d13C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
n <- which(error > q_C_d13C) # The rows with errors larger than this quantile

C <- subset(C,select = -n)
d13C <- subset(d13C,select = -n)
d14C <- subset(d14C,select = -n)
error_opt <- error[-n,1]

# The errors for the different evaluated criteria
error_OC_POC_opt <- relError_OC_POC[-n,1]
error_OC_MAOC_opt <- relError_OC_MAOC[-n,1]
error_d13C_POC_opt <- relError_d13C_POC[-n,1]
error_d13C_MAOC_opt <- relError_d13C_MAOC[-n,1]
error_d14C_POC_opt <- relError_d14C_POC[-n,1]
error_d14C_MAOC_opt <- relError_d14C_MAOC[-n,1]

POC_perc <- subset(POC_perc,select = -n)
bulkOC_perc <- subset(bulkOC_perc,select = -n)

# The total amounts of POC and bulkOC are calcultaed
POC_perc_tot <- unname(colSums(POC_perc))
bulkOC_perc_tot <- unname(colSums(bulkOC_perc))

# The errors per run are weighted based on the simulated amounts of POC and bulkOC
indErrors_OC_C_d13C <- ((POC_perc_tot * error_OC_POC_opt) + (bulkOC_perc_tot * error_OC_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d13C_C_d13C <- ((POC_perc_tot * error_d13C_POC_opt) + (bulkOC_perc_tot * error_d13C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d14C_C_d13C <- ((POC_perc_tot * error_d14C_POC_opt) + (bulkOC_perc_tot * error_d14C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)

# The average errors are calculated
avgError_OC_C_d13C <- mean(indErrors_OC_C_d13C)
avgError_d13C_C_d13C <- mean(indErrors_d13C_C_d13C)
avgError_d14C_C_d13C <- mean(indErrors_d14C_C_d13C)
totalError_C_d13C <- avgError_OC_C_d13C + avgError_d13C_C_d13C + avgError_d14C_C_d13C

# Q1 and Q3 are calculated
Q.OC.C_d13C <- c(summary(indErrors_OC_C_d13C)[["1st Qu."]], summary(indErrors_OC_C_d13C)[["3rd Qu."]])
Q.d13C.C_d13C <- c(summary(indErrors_d13C_C_d13C)[["1st Qu."]], summary(indErrors_d13C_C_d13C)[["3rd Qu."]])
Q.d14C.C_d13C <- c(summary(indErrors_d14C_C_d13C)[["1st Qu."]], summary(indErrors_d14C_C_d13C)[["3rd Qu."]])

# The distribution of the errors is plotted
# ggplot() +
#   geom_histogram(data = error, aes(error), binwidth = 0.005) +
#   geom_histogram(data = as.data.frame(error[-n,]), aes(error[-n,]), binwidth = 0.005, color = "green")

# Unrealistic d13C and d14C simulations are removed
# r <- which(d13C > -24, arr.ind = T)
# a <- unique(r[,2])
# d13C[,a] <- NA

# r <- which(d14C > 200, arr.ind = T)
# a <- unique(r[,2])
# d14C[,a] <- NA

# The average C, d13C and d14C depth profiles are calculated
df_avg <- data.frame(depth = depth,
                     avg_C = rowMeans(C,na.rm = T),
                     avg_d13C = rowMeans(d13C,na.rm = T),
                     avg_d14C = rowMeans(d14C,na.rm = T))

# The optimal profiles are stored
rowNum_opt <- which(error_opt == min(error_opt))
df_opt <- data.frame("depth" = depth,
                     "opt_C" = C[,rowNum_opt],
                     "opt_d13C" = d13C[,rowNum_opt],
                     "opt_d14C" = d14C[,rowNum_opt])

# The depth is added to the data frames
C$depth = depth
d13C$depth = depth
d14C$depth = depth

# The data frames are reformatted for plotting
df_C <- melt(C, id.vars = "depth")
df_d13C <- melt(d13C, id.vars = "depth")
df_d14C <- melt(d14C, id.vars = "depth")

# Plotting
p4 <- plot_C(df_C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p5 <- plot_d13C(df_d13C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p6 <- plot_d14C(df_d14C, df_avg, df_opt, measuredData_soil, color_lines = "grey75", color_avg = "grey40", par)

p4 <- p4 + annotate("label", x = 9.7, y = 0.57, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = 9.7, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_OC_C_d13C,3), " [", round2dec(Q.OC.C_d13C[1],3), ";", round2dec(Q.OC.C_d13C[2],3),"]"))

p5 <- p5 + annotate("label", x = -27.9, y = 0.57, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 0, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = -27.9, y = 0.5, hjust = 0, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d13C_C_d13C,3), " [", round2dec(Q.d13C.C_d13C[1],3), ";", round2dec(Q.d13C.C_d13C[2],3),"]"))

p6 <- p6 + annotate("label", x = 180, y = 0.57, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 0.5) +
  annotate("text", x = 180, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d14C_C_d13C,3), " [", round2dec(Q.d14C.C_d13C[1],3), ";", round2dec(Q.d14C.C_d13C[2],3),"]"))

# p4 <- p4 + annotate("label", x = 3.5, y = 0.56, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 0, size = 4, fontface = 'italic')
# # p4 <- p4 + annotate("label", x = 3.5, y = 0.56, label = expression(paste("Calibrated based on C and ", delta^{13}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p5 <- p5 + annotate("label", x = -27.9, y = 0.56, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p6 <- p6 + annotate("label", x = -310, y = 0.56, label = expression(paste("Calibrated for C and ", delta^{13}, "C")), hjust = 0, size = 4, fontface = 'italic')

rm(list = c("df_C", "df_d13C", "df_d14C"))

# grid.arrange(p4,p5,p6, nrow = 1)

# ------------------------------------------------------------------------------
# C and d14C
# ------------------------------------------------------------------------------

# The data is loaded
load(dir_C_d14C)

C <- as.data.frame(t(outputForSeries$Ctot))
d13C <- as.data.frame(t(outputForSeries$d13C))
d14C <- as.data.frame(t(outputForSeries$d14C))
depth <- outputForSeries$midDepth
error <- as.data.frame(outputForSeries$error)
names(error) <- "error"

relError_OC_POC <- as.data.frame(outputForSeries$relError_OC_POC)
names(relError_OC_POC) <- "relError_OC_POC"

relError_OC_MAOC <- as.data.frame(outputForSeries$relError_OC_MAOC)
names(relError_OC_MAOC) <- "relError_OC_MAOC"

relError_d13C_POC <- as.data.frame(outputForSeries$relError_d13C_POC)
names(relError_d13C_POC) <- "relError_d13C_POC"

relError_d13C_MAOC <- as.data.frame(outputForSeries$relError_d13C_MAOC)
names(relError_d13C_MAOC) <- "relError_d13C_MAOC"

relError_d14C_POC <- as.data.frame(outputForSeries$relError_d14C_POC)
names(relError_d14C_POC) <- "relError_d14C_POC"

relError_d14C_MAOC <- as.data.frame(outputForSeries$relError_d14C_MAOC)
names(relError_d14C_MAOC) <- "relError_d14C_MAOC"

POC_perc <- as.data.frame(t(outputForSeries$POC_perc))
bulkOC_perc <- as.data.frame(t(outputForSeries$bulkOC_perc))

# # Results with a high error are removed
# errorTol_C_d14C <- 0.8
# n <- which(error > errorTol_C_d14C)

# The X % lowest errors are retained
q_C_d14C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
n <- which(error > q_C_d14C) # The rows with errors larger than this quantile

C <- subset(C,select = -n)
d13C <- subset(d13C,select = -n)
d14C <- subset(d14C,select = -n)
error_opt <- error[-n,1]

# The errors for the different evaluated criteria
error_OC_POC_opt <- relError_OC_POC[-n,1]
error_OC_MAOC_opt <- relError_OC_MAOC[-n,1]
error_d13C_POC_opt <- relError_d13C_POC[-n,1]
error_d13C_MAOC_opt <- relError_d13C_MAOC[-n,1]
error_d14C_POC_opt <- relError_d14C_POC[-n,1]
error_d14C_MAOC_opt <- relError_d14C_MAOC[-n,1]

POC_perc <- subset(POC_perc,select = -n)
bulkOC_perc <- subset(bulkOC_perc,select = -n)

# The total amounts of POC and bulkOC are calcultaed
POC_perc_tot <- unname(colSums(POC_perc))
bulkOC_perc_tot <- unname(colSums(bulkOC_perc))

# The errors per run are weighted based on the simulated amounts of POC and bulkOC
indErrors_OC_C_d14C <- ((POC_perc_tot * error_OC_POC_opt) + (bulkOC_perc_tot * error_OC_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d13C_C_d14C <- ((POC_perc_tot * error_d13C_POC_opt) + (bulkOC_perc_tot * error_d13C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d14C_C_d14C <- ((POC_perc_tot * error_d14C_POC_opt) + (bulkOC_perc_tot * error_d14C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)

# The average errors are calculated
avgError_OC_C_d14C <- mean(indErrors_OC_C_d14C)
avgError_d13C_C_d14C <- mean(indErrors_d13C_C_d14C)
avgError_d14C_C_d14C <- mean(indErrors_d14C_C_d14C)
totalError_C_d14C <- avgError_OC_C_d14C + avgError_d13C_C_d14C + avgError_d14C_C_d14C

# Q1 and Q3 are calculated
Q.OC.C_d14C <- c(summary(indErrors_OC_C_d14C)[["1st Qu."]], summary(indErrors_OC_C_d14C)[["3rd Qu."]])
Q.d13C.C_d14C <- c(summary(indErrors_d13C_C_d14C)[["1st Qu."]], summary(indErrors_d13C_C_d14C)[["3rd Qu."]])
Q.d14C.C_d14C <- c(summary(indErrors_d14C_C_d14C)[["1st Qu."]], summary(indErrors_d14C_C_d14C)[["3rd Qu."]])

# The distribution of the errors is plotted
# ggplot() +
#   geom_histogram(data = error, aes(error), binwidth = 0.005) +
#   geom_histogram(data = as.data.frame(error[-n,]), aes(error[-n,]), binwidth = 0.005, color = "green")

# Unrealistic d13C and d14C simulations are removed
# r <- which(d13C > -24, arr.ind = T)
# a <- unique(r[,2])
# d13C[,a] <- NA

# r <- which(d14C > 200, arr.ind = T)
# a <- unique(r[,2])
# d14C[,a] <- NA

# The average C, d13C and d14C depth profiles are calculated
df_avg <- data.frame(depth = depth,
                     avg_C = rowMeans(C,na.rm = T),
                     avg_d13C = rowMeans(d13C,na.rm = T),
                     avg_d14C = rowMeans(d14C,na.rm = T))

# The optimal profiles are stored
rowNum_opt <- which(error_opt == min(error_opt))
df_opt <- data.frame("depth" = depth,
                     "opt_C" = C[,rowNum_opt],
                     "opt_d13C" = d13C[,rowNum_opt],
                     "opt_d14C" = d14C[,rowNum_opt])

# The depth is added to the data frames
C$depth = depth
d13C$depth = depth
d14C$depth = depth

# The data frames are reformatted for plotting
df_C <- melt(C, id.vars = "depth")
df_d13C <- melt(d13C, id.vars = "depth")
df_d14C <- melt(d14C, id.vars = "depth")

# Plotting
p7 <- plot_C(df_C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p8 <- plot_d13C(df_d13C, df_avg, df_opt, measuredData_soil, color_lines = "grey75", color_avg = "grey40", par)
p9 <- plot_d14C(df_d14C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)

p7 <- p7 + annotate("label", x = 9.7, y = 0.57, label = expression(paste("Calibrated for C and ", Delta^{14}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = 9.7, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_OC_C_d14C,3), " [", round2dec(Q.OC.C_d14C[1],3), ";", round2dec(Q.OC.C_d14C[2],3),"]"))

p8 <- p8 + annotate("label", x = -27.9, y = 0.57, label = expression(paste("Calibrated for C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = -27.9, y = 0.5, hjust = 0, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d13C_C_d14C,3), " [", round2dec(Q.d13C.C_d14C[1],3), ";", round2dec(Q.d13C.C_d14C[2],3),"]"))

p9 <- p9 + annotate("label", x = 180, y = 0.57, label = expression(paste("Calibrated for C and ", Delta^{14}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = 180, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d14C_C_d14C,3), " [", round2dec(Q.d14C.C_d14C[1],3), ";", round2dec(Q.d14C.C_d14C[2],3),"]"))

# p7 <- p7 + annotate("label", x = 3.5, y = 0.56, label = expression(paste("Calibrated based on C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p8 <- p8 + annotate("label", x = -27.9, y = 0.56, label = expression(paste("Calibrated based on C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p9 <- p9 + annotate("label", x = -310, y = 0.56, label = expression(paste("Calibrated based on C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')

rm(list = c("df_C", "df_d13C", "df_d14C"))

# grid.arrange(p7,p8,p9, nrow = 1)

# ------------------------------------------------------------------------------
# C, d13C, d14C
# ------------------------------------------------------------------------------

# The data is loaded
load(dir_C_d13C_d14C)

C <- as.data.frame(t(outputForSeries$Ctot))
d13C <- as.data.frame(t(outputForSeries$d13C))
d14C <- as.data.frame(t(outputForSeries$d14C))
depth <- outputForSeries$midDepth
error <- as.data.frame(outputForSeries$error)
names(error) <- "error"

relError_OC_POC <- as.data.frame(outputForSeries$relError_OC_POC)
names(relError_OC_POC) <- "relError_OC_POC"

relError_OC_MAOC <- as.data.frame(outputForSeries$relError_OC_MAOC)
names(relError_OC_MAOC) <- "relError_OC_MAOC"

relError_d13C_POC <- as.data.frame(outputForSeries$relError_d13C_POC)
names(relError_d13C_POC) <- "relError_d13C_POC"

relError_d13C_MAOC <- as.data.frame(outputForSeries$relError_d13C_MAOC)
names(relError_d13C_MAOC) <- "relError_d13C_MAOC"

relError_d14C_POC <- as.data.frame(outputForSeries$relError_d14C_POC)
names(relError_d14C_POC) <- "relError_d14C_POC"

relError_d14C_MAOC <- as.data.frame(outputForSeries$relError_d14C_MAOC)
names(relError_d14C_MAOC) <- "relError_d14C_MAOC"

POC_perc <- as.data.frame(t(outputForSeries$POC_perc))
bulkOC_perc <- as.data.frame(t(outputForSeries$bulkOC_perc))

# # Results with a high error are removed
# errorTol_C_d13C_d14C <- 1.2
# n <- which(error > errorTol_C_d13C_d14C)

# The X % lowest errors are retained
q_C_d13C_d14C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
n <- which(error > q_C_d13C_d14C) # The rows with errors larger than this quantile

C <- subset(C,select = -n)
d13C <- subset(d13C,select = -n)
d14C <- subset(d14C,select = -n)
error_opt <- error[-n,1]

# The errors for the different evaluated criteria
error_OC_POC_opt <- relError_OC_POC[-n,1]
error_OC_MAOC_opt <- relError_OC_MAOC[-n,1]
error_d13C_POC_opt <- relError_d13C_POC[-n,1]
error_d13C_MAOC_opt <- relError_d13C_MAOC[-n,1]
error_d14C_POC_opt <- relError_d14C_POC[-n,1]
error_d14C_MAOC_opt <- relError_d14C_MAOC[-n,1]

POC_perc <- subset(POC_perc,select = -n)
bulkOC_perc <- subset(bulkOC_perc,select = -n)

# The total amounts of POC and bulkOC are calcultaed
POC_perc_tot <- unname(colSums(POC_perc))
bulkOC_perc_tot <- unname(colSums(bulkOC_perc))

# The errors per run are weighted based on the simulated amounts of POC and bulkOC
indErrors_OC_C_d13C_d14C <- ((POC_perc_tot * error_OC_POC_opt) + (bulkOC_perc_tot * error_OC_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d13C_C_d13C_d14C <- ((POC_perc_tot * error_d13C_POC_opt) + (bulkOC_perc_tot * error_d13C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)
indErrors_d14C_C_d13C_d14C <- ((POC_perc_tot * error_d14C_POC_opt) + (bulkOC_perc_tot * error_d14C_MAOC_opt)) / (POC_perc_tot + bulkOC_perc_tot)

# The average errors are calculated
avgError_OC_C_d13C_d14C <- mean(indErrors_OC_C_d13C_d14C)
avgError_d13C_C_d13C_d14C <- mean(indErrors_d13C_C_d13C_d14C)
avgError_d14C_C_d13C_d14C <- mean(indErrors_d14C_C_d13C_d14C)
totalError_C_d13C_d14C <- avgError_OC_C_d13C_d14C + avgError_d13C_C_d13C_d14C + avgError_d14C_C_d13C_d14C

# Q1 and Q3 are calculated
Q.OC.C_d13C_d14C <- c(summary(indErrors_OC_C_d13C_d14C)[["1st Qu."]], summary(indErrors_OC_C_d13C_d14C)[["3rd Qu."]])
Q.d13C.C_d13C_d14C <- c(summary(indErrors_d13C_C_d13C_d14C)[["1st Qu."]], summary(indErrors_d13C_C_d13C_d14C)[["3rd Qu."]])
Q.d14C.C_d13C_d14C <- c(summary(indErrors_d14C_C_d13C_d14C)[["1st Qu."]], summary(indErrors_d14C_C_d13C_d14C)[["3rd Qu."]])

# The distribution of the errors is plotted
# ggplot() +
#   geom_histogram(data = error, aes(error), binwidth = 0.005) +
#   geom_histogram(data = as.data.frame(error[-n,]), aes(error[-n,]), binwidth = 0.005, color = "green")

# Unrealistic d13C and d14C simulations are removed
# r <- which(d13C > -24, arr.ind = T)
# a <- unique(r[,2])
# d13C[,a] <- NA

# r <- which(d14C > 200, arr.ind = T)
# a <- unique(r[,2])
# d14C[,a] <- NA

# The average C, d13C and d14C depth profiles are calculated
df_avg <- data.frame(depth = depth,
                     avg_C = rowMeans(C,na.rm = T),
                     avg_d13C = rowMeans(d13C,na.rm = T),
                     avg_d14C = rowMeans(d14C,na.rm = T))

# The optimal profiles are stored
rowNum_opt <- which(error_opt == min(error_opt))
df_opt <- data.frame("depth" = depth,
                     "opt_C" = C[,rowNum_opt],
                     "opt_d13C" = d13C[,rowNum_opt],
                     "opt_d14C" = d14C[,rowNum_opt])

# The depth is added to the data frames
C$depth = depth
d13C$depth = depth
d14C$depth = depth

# The data frames are reformatted for plotting
df_C <- melt(C, id.vars = "depth")
df_d13C <- melt(d13C, id.vars = "depth")
df_d14C <- melt(d14C, id.vars = "depth")

# Plotting
p10 <- plot_C(df_C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p11 <- plot_d13C(df_d13C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p12 <- plot_d14C(df_d14C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)

p10 <- p10 + annotate("label", x = 9.7, y = 0.57, label = expression(paste("Calibrated for C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = 9.7, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_OC_C_d13C_d14C,3), " [", round2dec(Q.OC.C_d13C_d14C[1],3), ";", round2dec(Q.OC.C_d13C_d14C[2],3),"]"))

p11 <- p11 + annotate("label", x = -27.9, y = 0.57, label = expression(paste("Calibrated for C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = -27.9, y = 0.5, hjust = 0, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d13C_C_d13C_d14C,3), " [", round2dec(Q.d13C.C_d13C_d14C[1],3), ";", round2dec(Q.d13C.C_d13C_d14C[2],3),"]"))

p12 <- p12 + annotate("label", x = 180, y = 0.57, label = expression(paste("Calibrated for C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 1, size = 4, fontface = 'italic', alpha = 1) +
  annotate("text", x = 180, y = 0.5, hjust = 1, "label" = paste0("\u03F5\u0304", " = ", round2dec(avgError_d14C_C_d13C_d14C,3), " [", round2dec(Q.d14C.C_d13C_d14C[1],3), ";", round2dec(Q.d14C.C_d13C_d14C[2],3),"]"))

# p10 <- p10 + annotate("label", x = 1.9, y = 0.56, label = expression(paste("Calibrated based on C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p11 <- p11 + annotate("label", x = -27.9, y = 0.56, label = expression(paste("Calibrated based on C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
# p12 <- p12 + annotate("label", x = -400, y = 0.56, label = expression(paste("Calibrated based on C, ", delta^{13}, "C and ", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic', alpha = .75)

removeAxisLabel <- function(plotObj, ax){
  if(ax == "x"){
    plotObj <- plotObj +
      theme(axis.title.x = element_blank())
  } else if(ax == "y"){
    plotObj <- plotObj +
      theme(axis.title.y = element_blank())
  } else if (ax == "xy"){
    plotObj <- plotObj +
      theme(axis.title.x = element_blank(),
            axis.title.y = element_blank())
  }
  return(plotObj)
}

t <- 0.6
r <- 0.5
b <- 0
l <- 0.5

p1 <- p1 + theme(plot.margin = margin(0.3,r,b,l, "cm"))
p2 <- removeAxisLabel(p2,"y") + theme(plot.margin = margin(0,r,b,l, "cm"))
p3 <- removeAxisLabel(p3,"y") + theme(plot.margin = margin(0,r,b,l, "cm"))
p4 <- removeAxisLabel(p4,"x") + theme(plot.margin = margin(t,r,b,l, "cm"))
p5 <- removeAxisLabel(p5,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))
p6 <- removeAxisLabel(p6,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))
p7 <- removeAxisLabel(p7,"x") + theme(plot.margin = margin(t,r,b,l, "cm"))
p8 <- removeAxisLabel(p8,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))
p9 <- removeAxisLabel(p9,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))
p10 <- removeAxisLabel(p10,"x") + theme(plot.margin = margin(t,r,b,l, "cm"))
p11 <- removeAxisLabel(p11,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))
p12 <- removeAxisLabel(p12,"xy") + theme(plot.margin = margin(t,r,b,l, "cm"))

rm(list = c("df_C", "df_d13C", "df_d14C"))

# All plots are shown together
quartz(title = "Depth profiles optimal parameter values", width = 32*0.3937, height = 32*0.3937)
ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12, ncol = 3, nrow = 4, labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)", "(j)", "(k)", "(l)"))

# g <- ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12, ncol = 3, nrow = 4, labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)", "(j)", "(k)", "(l)"), font.label=list(color="black",size=18, face = "plain"))
# ggsave("Figures for manuscript/Figure_calibrationConstraints.png", g, device = "png", dpi = 300, units = "cm", width = 32, height = 32)

# grid.arrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11,p12, nrow = 4)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# The turnover times of microbes are are plotted
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

# ------------------------------------------------------------------------------
# C only calibration
# ------------------------------------------------------------------------------

# Defining some functions

get_TT_data <- function(outputForSeries, Qthreshold){
  
  death_rhizo_day <- as.data.frame(t(outputForSeries$death_R_day))
  death_bulk_day <- as.data.frame(t(outputForSeries$death_B_day))
  error <- as.data.frame(outputForSeries$error)
  names(error) <- "error"
  
  # The X % lowest errors are retained
  q_C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
  n <- which(error > q_C) # The rows with errors larger than this quantile
  
  death_rhizo_day <- subset(death_rhizo_day,select = -n)
  death_bulk_day <- subset(death_bulk_day,select = -n)
  
  # The average turnover times are calculated
  df_avg <- data.frame(depth = depth,
                       avg_rhizo = rowMeans(death_rhizo_day,na.rm = T),
                       avg_bulk = rowMeans(death_bulk_day,na.rm = T))
  
  # The depth is added to the data frames
  death_rhizo_day$depth = depth
  death_bulk_day$depth = depth
  
  # The data frames are reformatted for plotting
  df_death_rhizo_day <- melt(death_rhizo_day, id.vars = "depth")
  df_death_bulk_day <- melt(death_bulk_day, id.vars = "depth")
  
  return(list(df_death_rhizo_day, df_death_bulk_day, df_avg))
  
}

# A function to plot the microbial turnover times - rhizosphere
plot_TT_rhizo <- function(df, df_avg, color_lines, color_avg, xMax){
  
  p <- ggplot(df, aes(x = value, y = depth, colour = variable)) +
    geom_path() +
    scale_y_continuous(expand = c(0,0), trans = "reverse", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_x_continuous(expand = c(0,0), position = "top", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_colour_manual(labels = unique(df$variable),
                        values = rep(color_lines,length(unique(df$variable)))) +
    coord_cartesian(xlim = c(0, xMax), ylim = c(df_avg$depth[9], 0))+
    geom_path(data = df_avg, aes(x = avg_rhizo, y = depth), size = 2, colour = color_avg) +
    theme_classic() +
    theme(legend.position = "none",
          panel.border = element_rect(color = "black",
                                      fill = NA,
                                      size = 1),
          # axis.line = element_line(colour = "black", 
          #                          size = 0),
          axis.text=element_text(size = 16),
          axis.title=element_text(size = 16),
          plot.margin = margin(0, 5, 0, 0, "mm")) + # #top, right, bottom, left
    labs(x = "Turnover time of rhizosphere microbes (days)")
  
  return(p)
  
}

# A function to plot the microbial turnover times - bulk soil
plot_TT_bulk <- function(df, df_avg, color_lines, color_avg, xMax){
  
  p <- ggplot(df, aes(x = value, y = depth, colour = variable)) +
    geom_path() +
    scale_y_continuous(expand = c(0,0), trans = "reverse", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_x_continuous(expand = c(0,0), position = "top", sec.axis = sec_axis(~ ., labels = NULL)) +
    scale_colour_manual(labels = unique(df$variable),
                        values = rep(color_lines,length(unique(df$variable)))) +
    coord_cartesian(xlim = c(0, xMax), ylim = c(df_avg$depth[9], 0))+
    geom_path(data = df_avg, aes(x = avg_bulk, y = depth), size = 2, colour = color_avg) +
    theme_classic() +
    theme(legend.position = "none",
          panel.border = element_rect(color = "black",
                                      fill = NA,
                                      size = 1),
          # axis.line = element_line(colour = "black", 
          #                          size = 0),
          axis.text=element_text(size = 16),
          axis.title=element_text(size = 16),
          plot.margin = margin(0, 5, 0, 0, "mm")) + # #top, right, bottom, left
    labs(x = "Turnover time of bulk soil microbes (days)")
  
  return(p)
  
}

# The data is loaded
load(dir_C)

out_fun <- get_TT_data(outputForSeries, Qthreshold)
df_death_rhizo_day <- out_fun[[1]]
df_death_bulk_day <- out_fun[[2]]
df_avg <- out_fun[[3]]

p_rhizo_C <- plot_TT_rhizo(df_death_rhizo_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)
p_bulk_C <- plot_TT_bulk(df_death_bulk_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)

p_rhizo_C <- p_rhizo_C + annotate("label", x = 7.5, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic')
p_bulk_C <- p_bulk_C + annotate("label", x = 2e4, y = 0.57, label = "Calibrated for C", hjust = 0, size = 4, fontface = 'italic')

# ------------------------------------------------------------------------------
# Calibration for C and d13C
# ------------------------------------------------------------------------------

load(dir_C_d13C)

out_fun <- get_TT_data(outputForSeries, Qthreshold)
df_death_rhizo_day <- out_fun[[1]]
df_death_bulk_day <- out_fun[[2]]
df_avg <- out_fun[[3]]

p_rhizo_C_d13C <- plot_TT_rhizo(df_death_rhizo_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)
p_bulk_C_d13C <- plot_TT_bulk(df_death_bulk_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)

p_rhizo_C_d13C <- p_rhizo_C_d13C + annotate("label", x = 7.5, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic')
p_bulk_C_d13C <- p_bulk_C_d13C + annotate("label", x = 2e4, y = 0.57, label = "Calibrated for C", hjust = 0, size = 4, fontface = 'italic')

# ------------------------------------------------------------------------------
# Calibration for C and d14C
# ------------------------------------------------------------------------------

load(dir_C_d14C)

out_fun <- get_TT_data(outputForSeries, Qthreshold)
df_death_rhizo_day <- out_fun[[1]]
df_death_bulk_day <- out_fun[[2]]
df_avg <- out_fun[[3]]

p_rhizo_C_d14C <- plot_TT_rhizo(df_death_rhizo_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)
p_bulk_C_d14C <- plot_TT_bulk(df_death_bulk_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)

p_rhizo_C_d14C <- p_rhizo_C_d14C + annotate("label", x = 7.5, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic')
p_bulk_C_d14C <- p_bulk_C_d14C + annotate("label", x = 2e4, y = 0.57, label = "Calibrated for C", hjust = 0, size = 4, fontface = 'italic')

# ------------------------------------------------------------------------------
# Calibration for C, d13C and d14C
# ------------------------------------------------------------------------------

load(dir_C_d13C_d14C)

out_fun <- get_TT_data(outputForSeries, Qthreshold)
df_death_rhizo_day <- out_fun[[1]]
df_death_bulk_day <- out_fun[[2]]
df_avg <- out_fun[[3]]

p_rhizo_C_d13C_d14C <- plot_TT_rhizo(df_death_rhizo_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)
p_bulk_C_d13C_d14C <- plot_TT_bulk(df_death_bulk_day, df_avg, color_lines = "deepskyblue3", color_avg = "dodgerblue4", xMax = 10)

p_rhizo_C_d13C_d14C <- p_rhizo_C_d13C_d14C + annotate("label", x = 7.5, y = 0.57, label = "Calibrated for C", hjust = 1, size = 4, fontface = 'italic')
p_bulk_C_d13C_d14C <- p_bulk_C_d13C_d14C + annotate("label", x = 2e4, y = 0.57, label = "Calibrated for C", hjust = 0, size = 4, fontface = 'italic')

# ---------------------
# All plots are grouped
# ---------------------

quartz(title = "Turnover times microbes", width = 32*0.3937, height = 32*0.3937)
ggarrange(p_rhizo_C, p_bulk_C, p_rhizo_C_d13C, p_bulk_C_d13C, p_rhizo_C_d14C, p_bulk_C_d14C, p_rhizo_C_d13C_d14C, p_bulk_C_d13C_d14C, ncol = 2, nrow = 4, labels = c("(A)", "(B)", "(C)", "(D)", "(E)", "(F)", "(G)", "(H)"))

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# The ranges in parameter values are plotted
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

# # The DE ouptput it loaded, to plot the optimal values
load("GA output/Hainich/17. New calibration/C/Soil/DE_out.RData")
outDE_C <- data.frame("bestmem" = outDE$optim$bestmem)
load("GA output/Hainich/17. New calibration/C_d13C/Soil/DE_out.RData")
outDE_C_d13C <- data.frame("bestmem" = outDE$optim$bestmem)
load("GA output/Hainich/17. New calibration/C_d14C/Soil/DE_out.RData")
outDE_C_d14C <- data.frame("bestmem" = outDE$optim$bestmem)
load("GA output/Hainich/17. New calibration/C_d13C_d14C/Soil/DE_out.RData")
outDE_C_d13C_d14C <- data.frame("bestmem" = outDE$optim$bestmem)

rm(outDE)

# The best values are combined in a data.frame
df_bestMem <- data.frame("C" = outDE_C$bestmem,
                         "C_d13C" = outDE_C_d13C$bestmem,
                         "C_d14C" = outDE_C_d14C$bestmem,
                         "C_d13C_d14C" = outDE_C_d13C_d14C$bestmem)

rownames(df_bestMem) <- rownames(outDE_C)

# The data for the unique parameter values is loaded and reformatted
load(dir_C)
output_C <- outputForSeries

load(dir_C_d13C)
output_C_d13C <- outputForSeries

load(dir_C_d14C)
output_C_d14C <- outputForSeries

load(dir_C_d13C_d14C)
output_C_d13C_d14C <- outputForSeries

rm(outputForSeries)

# Only parameter values for runs with an error smaller than the threshold are retained
n_C <- which(output_C$error < q_C)
n_C_d13C <- which(output_C_d13C$error < q_C_d13C)
n_C_d14C <- which(output_C_d14C$error < q_C_d14C)
n_C_d13C_d14C <- which(output_C_d13C_d14C$error < q_C_d13C_d14C)

paramValues_optimal_C <- as.data.frame(output_C$paramValues[n_C,])
paramValues_optimal_C_d13C <- as.data.frame(output_C_d13C$paramValues[n_C_d13C,])
paramValues_optimal_C_d14C <- as.data.frame(output_C_d14C$paramValues[n_C_d14C,])
paramValues_optimal_C_d13C_d14C <- as.data.frame(output_C_d13C_d14C$paramValues[n_C_d13C_d14C,])

# The errors are added
paramValues_optimal_C$error <- output_C$error[n_C]
paramValues_optimal_C_d13C$error <- output_C_d13C$error[n_C_d13C]
paramValues_optimal_C_d14C$error <- output_C_d14C$error[n_C_d14C]
paramValues_optimal_C_d13C_d14C$error <- output_C_d13C_d14C$error[n_C_d13C_d14C]

# ----------------------------------------------------------
# A similar data.frame with all tested parameters is created
# ----------------------------------------------------------

allTestedParam_C <- as.data.frame(output_C$paramValues)
allTestedParam_C$error <- NA
allTestedParam_C$optimal <- 0
paramValues_optimal_C$optimal <- 1
allParam_C <- rbind(allTestedParam_C,paramValues_optimal_C)

allTestedParam_C_d13C <- as.data.frame(output_C_d13C$paramValues)
allTestedParam_C_d13C$error <- NA
allTestedParam_C_d13C$optimal <- 0
paramValues_optimal_C_d13C$optimal <- 1
allParam_C_d13C <- rbind(allTestedParam_C_d13C,paramValues_optimal_C_d13C)

allTestedParam_C_d14C <- as.data.frame(output_C_d14C$paramValues)
allTestedParam_C_d14C$error <- NA
allTestedParam_C_d14C$optimal <- 0
paramValues_optimal_C_d14C$optimal <- 1
allParam_C_d14C <- rbind(allTestedParam_C_d14C,paramValues_optimal_C_d14C)

allTestedParam_C_d13C_d14C <- as.data.frame(output_C_d13C_d14C$paramValues)
allTestedParam_C_d13C_d14C$error <- NA
allTestedParam_C_d13C_d14C$optimal <- 0
paramValues_optimal_C_d13C_d14C$optimal <- 1
allParam_C_d13C_d14C <- rbind(allTestedParam_C_d13C_d14C,paramValues_optimal_C_d13C_d14C)

# A function to group the parameters
groupPerParameter <- function(data_C, data_d13C, data_d14C, data_d13C_d14C){
  
  # The data for this variable is isolated
  subset_C <- as.data.frame(data_C[,1])
  subset_C$ID <- "C"
  subset_C$error <- data_C$error
  subset_C$optimal <- data_C$optimal
  names(subset_C) <- c("value", "ID", "error", "optimal")
  
  subset_C_d13C <- as.data.frame(data_d13C[,1])
  subset_C_d13C$ID <- "C_d13C"
  subset_C_d13C$error <- data_d13C$error
  subset_C_d13C$optimal <- data_d13C$optimal
  names(subset_C_d13C) <- c("value", "ID", "error", "optimal")
  
  subset_C_d14C <- as.data.frame(data_d14C[,1])
  subset_C_d14C$ID <- "C_d14C"
  subset_C_d14C$error <- data_d14C$error
  subset_C_d14C$optimal <- data_d14C$optimal
  names(subset_C_d14C) <- c("value", "ID", "error", "optimal")
  
  subset_C_d13C_d14C <- as.data.frame(data_d13C_d14C[,1])
  subset_C_d13C_d14C$ID <- "C_d13C_d14C"
  subset_C_d13C_d14C$error <- data_d13C_d14C$error
  subset_C_d13C_d14C$optimal <- data_d13C_d14C$optimal
  names(subset_C_d13C_d14C) <- c("value", "ID", "error", "optimal")

  # The data is combined into a single data frame
  df <- data.frame("ID" = c(subset_C$ID, subset_C_d13C$ID, subset_C_d14C$ID, subset_C_d13C_d14C$ID),
                   "value" = c(subset_C$value, subset_C_d13C$value, subset_C_d14C$value, subset_C_d13C_d14C$value),
                   "error" = c(subset_C$error, subset_C_d13C$error, subset_C_d14C$error, subset_C_d13C_d14C$error),
                   "optimal" = c(subset_C$optimal, subset_C_d13C$optimal, subset_C_d14C$optimal, subset_C_d13C_d14C$optimal))
  
  # The ID's are converted to factors to determine the order in which they are plotted
  df$ID <- factor(df$ID, levels = c("C_d13C_d14C", "C_d14C", "C_d13C", "C"))
  
  return(df)

}

# The data for every parameter is extracted and re-formatted using the function
df_VmaxD_root_R <- groupPerParameter(allParam_C %>% select(VmaxD_root_R_stage1, error, optimal), allParam_C_d13C %>% select(VmaxD_root_R_stage1, error, optimal), allParam_C_d14C %>% select(VmaxD_root_R_stage1, error, optimal), allParam_C_d13C_d14C %>% select(VmaxD_root_R_stage1, error, optimal))
df_VmaxU_BioAv <- groupPerParameter(allParam_C %>% select(VmaxU_BioAv_stage1, error, optimal), allParam_C_d13C %>% select(VmaxU_BioAv_stage1, error, optimal), allParam_C_d14C %>% select(VmaxU_BioAv_stage1, error, optimal), allParam_C_d13C_d14C %>% select(VmaxU_BioAv_stage1, error, optimal))
df_Db0 <- groupPerParameter(allParam_C %>% select(Db0_stage1, error, optimal), allParam_C_d13C %>% select(Db0_stage1, error, optimal), allParam_C_d14C %>% select(Db0_stage1, error, optimal), allParam_C_d13C_d14C %>% select(Db0_stage1, error, optimal))
df_betaRoots <- groupPerParameter(allParam_C %>% select(betaRoots, error, optimal), allParam_C_d13C %>% select(betaRoots, error, optimal), allParam_C_d14C %>% select(betaRoots, error, optimal), allParam_C_d13C_d14C %>% select(betaRoots, error, optimal))
df_Km_ads <- groupPerParameter(allParam_C %>% select(Km_ads, error, optimal), allParam_C_d13C %>% select(Km_ads, error, optimal), allParam_C_d14C %>% select(Km_ads, error, optimal), allParam_C_d13C_d14C %>% select(Km_ads, error, optimal))
df_Km_depol_M <- groupPerParameter(allParam_C %>% select(Km_depol_M_stage1, error, optimal), allParam_C_d13C %>% select(Km_depol_M_stage1, error, optimal), allParam_C_d14C %>% select(Km_depol_M_stage1, error, optimal), allParam_C_d13C_d14C %>% select(Km_depol_M_stage1, error, optimal))
df_kDes_init <- groupPerParameter(allParam_C %>% select(kDes_init, error, optimal), allParam_C_d13C %>% select(kDes_init, error, optimal), allParam_C_d14C %>% select(kDes_init, error, optimal), allParam_C_d13C_d14C %>% select(kDes_init, error, optimal))
df_VmaxD_M <- groupPerParameter(allParam_C %>% select(VmaxD_M_stage1, error, optimal), allParam_C_d13C %>% select(VmaxD_M_stage1, error, optimal), allParam_C_d14C %>% select(VmaxD_M_stage1, error, optimal), allParam_C_d13C_d14C %>% select(VmaxD_M_stage1, error, optimal))
df_Vmax_ads <- groupPerParameter(allParam_C %>% select(Vmax_ads, error, optimal), allParam_C_d13C %>% select(Vmax_ads, error, optimal), allParam_C_d14C %>% select(Vmax_ads, error, optimal), allParam_C_d13C_d14C %>% select(Vmax_ads, error, optimal))

# The averages of the parameters that have a limited effect on model outputs
# avg_Db_efold <- mean(df_Db_eFold_depth$value)
# avg_vMax <- mean(df_Vmax_ads$value)
# avg_avectionRate <- mean(df_advectionRate_polyC$value)

# The colors for plotting are defined
c_C <- "#80962a"
c_C_d13C <- "#b85c16"
c_C_d14C <- "#bb005e"
c_C_d13C_d14C <- "#4a32a8"
colors <- c(c_C, c_C_d13C, c_C_d14C, c_C_d13C_d14C)

# A function to find the range in parameter values for all calibrations
findRange <- function(df){
  
  range_C <- findRange_single(df, "C")
  range_C_d13C <- findRange_single(df, "C_d13C")
  range_C_d14C <- findRange_single(df, "C_d14C")
  range_C_d13C_d14C <- findRange_single(df, "C_d13C_d14C")
  
  # The lowest value
  low <- min(c(range_C, range_C_d13C, range_C_d14C, range_C_d13C_d14C))
  
  out_C <- reduceDec(range_C/low,100)
  out_C_d13C <- reduceDec(range_C_d13C/low,100)
  out_C_d14C <- reduceDec(range_C_d14C/low,100)
  out_C_d13C_d14C <- reduceDec(range_C_d13C_d14C/low,100)
  
  return(as.list(c(out_C, out_C_d13C, out_C_d14C, out_C_d13C_d14C)))
  
}

# A function to find the range for a single calibration
findRange_single <- function(df, scen){
  
  nRow <- which(df$ID == scen)
  minValue <- min(df$value[nRow])
  maxValue <- max(df$value[nRow])
  range <- maxValue - minValue
  
  return(range)
  
}

df <- df_betaRoots
df_opt <- df[df$optimal == 1,]

summ <- summary(df_opt[df_opt$ID == "C","value"])
median <- summ[["Median"]]
IQ <- summ[["3rd Qu."]] - summ[["1st Qu."]]




# A function to find the median and interquartile range for the optimal parameters
find_IQ_med <- function(df){
  
  df_opt <- df[df$optimal == 1,]
  
  out_C <- find_IQ_med_single(df_opt,"C")
  out_C_d13C <- find_IQ_med_single(df_opt,"C_d13C")
  out_C_d14C <- find_IQ_med_single(df_opt,"C_d14C")
  out_C_d13C_d14C <- find_IQ_med_single(df_opt,"C_d13C_d14C")
  
  out <- c()
  out[["C.med"]] <- out_C[["med"]]
  out[["C.IQ"]] <- out_C[["IQ"]]
  out[["C_d13C.med"]] <- out_C_d13C[["med"]]
  out[["C_d13C.IQ"]] <- out_C_d13C[["IQ"]]
  out[["C_d14C.med"]] <- out_C_d14C[["med"]]
  out[["C_d14C.IQ"]] <- out_C_d14C[["IQ"]]
  out[["C_d13C_d14C.med"]] <- out_C_d13C_d14C[["med"]]
  out[["C_d13C_d14C.IQ"]] <- out_C_d13C_d14C[["IQ"]]

  return(out)
  
}

# A function to find the range for a single calibration
find_IQ_med_single <- function(df, scen){
  
  df_opt <- df[df$optimal == 1,]
  
  summ <- summary(df_opt[df_opt$ID == scen,"value"])
  median <- summ[["Median"]]
  IQ <- summ[["3rd Qu."]] - summ[["1st Qu."]]
  
  return(c("med" = median, "IQ" = IQ))
  
}

# A function to reduce the number of decimals
reduceDec <- function(input, factor){
  
  a <- input * factor
  a <- round(a)
  a <- a / factor
  
  return(a)
  
}

# A function to add the relative ranges to the graphs
addRanges <- function(plotObj, ranges, x){
  
  plotObj <- plotObj + 
    annotate("text", x = x, y = 1, label = as.character(ranges[[4]]), hjust = 0, size = 4, fontface = 'italic') +
    annotate("text", x = x, y = 2, label = as.character(ranges[[3]]), hjust = 0, size = 4, fontface = 'italic') +
    annotate("text", x = x, y = 3, label = as.character(ranges[[2]]), hjust = 0, size = 4, fontface = 'italic') +
    annotate("text", x = x, y = 4, label = as.character(ranges[[1]]), hjust = 0, size = 4, fontface = 'italic')
  
  return(plotObj)
  
}

# The relative range in parameter values is obtained
ranges_VmaxD_root_R <- findRange(df_VmaxD_root_R)
ranges_VmaxU_BioAv <- findRange(df_VmaxU_BioAv)
ranges_Db0 <- findRange(df_Db0)
ranges_betaRoots <- findRange(df_betaRoots)
ranges_Km_ads <- findRange(df_Km_ads)
ranges_Km_depol_M <- findRange(df_Km_depol_M)
ranges_kDes_init <- findRange(df_kDes_init)
ranges_VmaxD_M <- findRange(df_VmaxD_M)
ranges_Vmax_ads <- findRange(df_Vmax_ads)

# Find the median and interquartile ranges
IQ_med_VmaxD_root_R <- find_IQ_med(df_VmaxD_root_R)
IQ_med_VmaxU_BioAv <- find_IQ_med(df_VmaxU_BioAv)
IQ_med_Db0 <- find_IQ_med(df_Db0)
IQ_med_betaRoots <- find_IQ_med(df_betaRoots)
IQ_med_Vmax_ads <- find_IQ_med(df_Vmax_ads)
IQ_med_Km_ads <- find_IQ_med(df_Km_ads)
IQ_med_VmaxD_M <- find_IQ_med(df_VmaxD_M)
IQ_med_Km_depol_M <- find_IQ_med(df_Km_depol_M)
IQ_med_kDes_init <- find_IQ_med(df_kDes_init)

# A function to create the plotting object (boxplots)
createPlotObj <- function(dataset, allJitterColor, optimalJitterColor, boxplotAlpha, xMin, xMax, title, colors, bestMem){
  
  p <- ggplot() +
    # geom_violin(scale = "width",adjust = 1) +
    geom_jitter(data = dataset[dataset$optimal == 0,], aes(y = ID, x = value, fill = ID), shape=16, height = 0.3, cex=1, color = allJitterColor) +
    geom_jitter(data = dataset[dataset$optimal == 1,], aes(y = ID, x = value, fill = ID), shape=16, height = 0.3, cex=1, color = optimalJitterColor) +
    geom_boxplot(data = dataset[dataset$optimal == 1,], aes(y = ID, x = value, fill = ID), alpha = boxplotAlpha) +
    # geom_point(data = data.frame("x" = as.numeric(bestMem), "y" = c("C", "C_d13C", "C_d14C", "C_d13C_d14C")), aes(x=x, y=y), color = "#045a8d", size = 4, shape = 17) +
    # geom_dotplot(binaxis='x', stackdir='center', dotsize=1, binwidth = 0.00005) +
    # xlim(c(xMin, xMax)) +
    scale_x_continuous(expand = c(0,0), limits = c(xMin, xMax), sec.axis = sec_axis(~ .)) +
    scale_y_discrete(labels = c(expression(paste("C,", delta^{13}, "C and ", Delta^{14}, "C")), expression(paste("C and ", Delta^{14}, "C")), expression(paste("C and ", delta^{13}, "C")),"C")) +
    theme_classic() +
    guides(y.sec = "axis", ) +
    theme(axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          plot.title = element_text(hjust = 0.5),
          legend.position = "none",
          axis.text.y.right = element_blank(),
          plot.margin=unit(c(.5,.5,.5,.5),"cm")) +
    ggtitle(title) +
    scale_fill_manual(values = colors)
  
  return(p)
  
}

createPlotObj_medIQ <- function(dataset, allJitterColor, optimalJitterColor, boxplotAlpha, xMin, xMax, title, colors, bestMem, addMed = 0, IQ_med = NA, numDec = NA){
  
  xVal1 <- xMax + (xMax - xMin) * 2
  
  p <- ggplot() +
    geom_jitter(data = dataset[dataset$optimal == 0,], aes(y = ID, x = value, fill = ID), shape=16, height = 0.3, cex=1, color = allJitterColor) +
    geom_jitter(data = dataset[dataset$optimal == 1,], aes(y = ID, x = value, fill = ID), shape=16, height = 0.3, cex=1, color = optimalJitterColor) +
    geom_boxplot(data = dataset[dataset$optimal == 1,], aes(y = ID, x = value, fill = ID), alpha = boxplotAlpha) +
    scale_x_continuous(expand = c(0,0), limits = c(xMin, xVal1), sec.axis = sec_axis(~ .)) +
    scale_y_discrete(labels = c(expression(paste("C,", delta^{13}, "C and ", Delta^{14}, "C")), expression(paste("C and ", Delta^{14}, "C")), expression(paste("C and ", delta^{13}, "C")),"C")) +
    guides(y.sec = "axis", ) +
    theme(axis.title.x=element_blank(),
          axis.title.y=element_blank(),
          plot.title = element_text(hjust = 0.5),
          axis.text = element_text(size = 10),
          legend.position = "none",
          axis.text.y.right = element_blank(),
          panel.background = element_rect(fill = "white",
                                          colour = "white"),
          panel.border = element_rect(colour = "black", fill=NA, size=.5)) +
    ggtitle(title) +
    scale_fill_manual(values = colors)
  
  if(addMed == 1){
    
    xVal2 <- xMax + (xMax - xMin) / 20
    
    textObj <- c(paste0(round2dec(IQ_med$C.med,numDec), "\n(", round2dec(IQ_med$C.IQ,numDec), ")"),
                 paste0(round2dec(IQ_med$C_d13C.med,numDec), "\n(", round2dec(IQ_med$C_d13C.IQ,numDec), ")"),
                 paste0(round2dec(IQ_med$C_d14C.med,numDec), "\n(", round2dec(IQ_med$C_d14C.IQ,numDec), ")"),
                 paste0(round2dec(IQ_med$C_d13C_d14C.med,numDec), "\n(", round2dec(IQ_med$C_d13C_d14C.IQ,numDec), ")"))
    
    p <- p +
      geom_text(aes(x = xVal2,
                    y = unique(dataset$ID),
                    label = textObj),
                hjust = 0,
                size = 3.5) +
      coord_cartesian(clip = 'off',
                      xlim = c(xMin,xMax)) +
      theme(plot.margin = unit(c(0,4,1,1), "lines"))
  }
  
  return(p)
  
}

boxplotAlpha <- 0.8

c_allJitter <- "grey80"
c_optimalJitter <- "black"

# p1 <- createPlotObj(dataset = df_VmaxD_root_R, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 1, title = expression(italic(paste(V[maxD]) )), colors = colors, bestMem = df_bestMem["VmaxD_root_R_stage1",])
# p2 <- createPlotObj(dataset = df_VmaxU_BioAv, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 1, title = expression(italic(paste(V[maxU_BioAv]) )), colors = colors, bestMem = df_bestMem["VmaxU_BioAv_stage1",])
# p3 <- createPlotObj(dataset = df_Db0, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 7.5e-5, title = expression(italic(D[b](0))), colors = colors, bestMem = df_bestMem["Db0_stage1",])
# p4 <- createPlotObj(dataset = df_betaRoots, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.85, xMax = 0.93, title = expression(italic(beta[r])), colors = colors, bestMem = df_bestMem["betaRoots",])
# p5 <- createPlotObj(dataset = df_Vmax_ads, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.1, xMax = 1e3, title = expression(italic(V["max,ads"])), colors = colors, bestMem = df_bestMem["Vmax_ads",])
# p6 <- createPlotObj(dataset = df_Km_ads, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.001, xMax = 1, title = expression(italic(K["m,ads"])), colors = colors, bestMem = df_bestMem["Km_ads",])
# p7 <- createPlotObj(dataset = df_VmaxD_M, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.1, xMax = 1e3, title = expression(italic(V["max,poly-b"])), colors = colors, bestMem = df_bestMem["VmaxD_M_stage1",])
# p8 <- createPlotObj(dataset = df_Km_depol_M, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha , xMin = 1e-6, xMax = 1, title = expression(italic(K["m,depol-M"])), colors = colors, bestMem = df_bestMem["Km_depol_M_stage1",])
# p9 <- createPlotObj(dataset = df_kDes_init, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 1, title = expression(italic(k["deprotect"](0))), colors = colors, bestMem = df_bestMem["kDes_init",])

# Plotting with the median and interquartile range included
p1 <- createPlotObj_medIQ(dataset = df_VmaxD_root_R, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.5, xMax = 1, title = expression(italic(paste(V["max,POC-r"]) )), colors = colors, bestMem = df_bestMem["VmaxD_root_R_stage1",], addMed = 1, IQ_med = IQ_med_VmaxD_root_R, numDec = 2)
p2 <- createPlotObj_medIQ(dataset = df_VmaxU_BioAv, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 1, title = expression(italic(paste(V["maxU,mic-r"]) )), colors = colors, bestMem = df_bestMem["VmaxU_BioAv_stage1",], addMed = 1, IQ_med = IQ_med_VmaxU_BioAv, numDec = 2)
p3 <- createPlotObj_medIQ(dataset = df_Db0, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 3.5e-7, xMax = 1e-4, title = expression(italic(D[b](0))), colors = colors, bestMem = df_bestMem["Db0_stage1",], addMed = 1, IQ_med = IQ_med_Db0, numDec = 7)
p4 <- createPlotObj_medIQ(dataset = df_betaRoots, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.85, xMax = 0.97, title = expression(italic(beta[r])), colors = colors, bestMem = df_bestMem["betaRoots",], addMed = 1, IQ_med = IQ_med_betaRoots, numDec = 3)
p5 <- createPlotObj_medIQ(dataset = df_Vmax_ads, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.1, xMax = 1e3, title = expression(italic(V["max,ads"])), colors = colors, bestMem = df_bestMem["Vmax_ads",], addMed = 1, IQ_med = IQ_med_Vmax_ads, numDec = 0)
p6 <- createPlotObj_medIQ(dataset = df_Km_ads, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.001, xMax = 1, title = expression(italic(K["m,ads"])), colors = colors, bestMem = df_bestMem["Km_ads",], addMed = 1, IQ_med = IQ_med_Km_ads, numDec = 3)
p7 <- createPlotObj_medIQ(dataset = df_VmaxD_M, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0.1, xMax = 1e3, title = expression(italic(V["max,DOC-b"])), colors = colors, bestMem = df_bestMem["VmaxD_M_stage1",], addMed = 1, IQ_med = IQ_med_VmaxD_M, numDec = 0)
p8 <- createPlotObj_medIQ(dataset = df_Km_depol_M, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha , xMin = 1e-6, xMax = 1, title = expression(italic(K["m,DOC-b"])), colors = colors, bestMem = df_bestMem["Km_depol_M_stage1",], addMed = 1, IQ_med = IQ_med_Km_depol_M, numDec = 2)
p9 <- createPlotObj_medIQ(dataset = df_kDes_init, allJitterColor = c_allJitter, optimalJitterColor = c_optimalJitter, boxplotAlpha, xMin = 0, xMax = 1, title = expression(italic(k["deprotect"](0))), colors = colors, bestMem = df_bestMem["kDes_init",], addMed = 1, IQ_med = IQ_med_kDes_init, numDec = 2)


# All plots are shown together
quartz(title = "Depth profiles optimal parameter values", width = 32*0.3937, height = 32*0.3937)
ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,ncol = 3,nrow = 3, labels = c("(A)", "(B)", "(C)", "(D)", "(E)", "(F)", "(G)", "(H)", "(I)"))

# g <- ggarrange(p1,p2,p3,p4,p5,p6,p7,p8,p9,ncol = 3,nrow = 3, labels = c("(a)", "(b)", "(c)", "(d)", "(e)", "(f)", "(g)", "(h)", "(i)"), font.label=list(color="black",size=18, face = "plain"))
# ggsave("Figures for manuscript/Figure_ranges_optimalValues.png", g, device = "png", dpi = 300, units = "cm", width = 32, height = 32)

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Plotting the error of every parameter value
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

df_Db0$relError <- df_Db0$error / max(df_Db0$error)


ggplot() + 
  geom_point(data = subset(df_Db0, df_Db0$ID == "C"), aes(x = value, y = relError/max(relError)), color = "red") +
  # geom_ma(data = subset(df_Db0, df_Db0$ID == "C"), aes(x = value, y = relError/max(relError)), color = "red", ma_fun = SMA, n = 30)
  geom_line(data = subset(df_Db0, df_Db0$ID == "C"), aes(x = value[order(relError)], y = rollmean(relError/max(relError),21 , fill = NA)), color = "blue") +
  geom_point(data = subset(df_Db0, df_Db0$ID == "C_d13C"), aes(x = value, y = relError/max(relError)), color = "blue") +
  geom_point(data = subset(df_Db0, df_Db0$ID == "C_d14C"), aes(x = value, y = relError/max(relError)), color = "green") +
  geom_point(data = subset(df_Db0, df_Db0$ID == "C_d13C_d14C"), aes(x = value, y = relError/max(relError)), color = "red")
  

ggplot() + 
  geom_point(data = df_Db0, aes(x = value, y = relError, color = ID))


a <- subset(df_Db0, df_Db0$ID == "C")
relError <- a$relError / max(a$relError)

y  <-  rollmean(relError/max(relError),21 , fill = NA)
x <- a$value


# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# Correlations between parameter values are plotted
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------

library("RColorBrewer")

# The colnames are changed for the corrplots
# IDs <- c("$D[b](0)", "$D[b_depth]", "$beta[r]", "$V[max_ads]", "$K[m_ads]", "$V[maxD_M]", "$K[m_depol_M]", "$k[deprotect](0)")
IDs <- c("$V[max_POC-r]",
         "$V[maxU_mic-r]",
         "$D[b](0)",
         "$beta[r]",
         "$K[m_ads]",
         "$K[m_DOC-b]",
         "$k[deprotect](0)",
         "$V[max_DOC-b]",
         "$V[max_ads]"
         )

colnames(paramValues_optimal_C) <- IDs
colnames(paramValues_optimal_C_d13C) <- IDs
colnames(paramValues_optimal_C_d14C) <- IDs
colnames(paramValues_optimal_C_d13C_d14C) <- IDs

# A function to create the correlation plots
makeCorrPlot <- function(inData){
  
  testRes = cor.mtest(inData, conf.level = 0.95)
  plotObj <- corrplot(cor(inData),
                       p.mat = testRes$p,
                       diag = FALSE,
                       insig='blank',
                       method = "square",
                       type="lower", 
                       order="alphabet",
                       tl.col = 'black',
                       tl.srt = 45,
                       # col = COL2("BrBG"),
                      col = COL2("RdBu"),
                       # col=brewer.pal(n=8, name="RdYlBu"),
                       addCoef.col = "black")$corrPos -> p1
  text(p1$x, p1$y, round(p1$corr, 2))
  
  return(plotObj)
  
}

# ---------------------------
# Plots for C and C_d13C_d14C
# ---------------------------

quartz(title = "Correlation plots", width = 32*0.3937, height = 16*0.3937)
par(mfrow = c(1, 2))

makeCorrPlot(paramValues_optimal_C[,seq(1,9)])
mtext("(A) Calibrated based on C", side=3, line=0, cex=1.5, las=1, adj = 0.05, padj = -2)

makeCorrPlot(paramValues_optimal_C_d13C_d14C[,seq(1,9)])
mtext(expression(paste("(B) Calibrated based on C and ", delta^{13}, "C and ", Delta^{14}, "C")), side=3, line=0, cex=1.5, las=1, adj = 0.05, padj = -1)
par(mfrow = c(1, 1)) #To clear layout

# quartz.save("Figures for manuscript/Figure_corrPlot.png", type="png")

# ---------------------------
# Plots for C_d13C and C_d14C
# ---------------------------

quartz(title = "Correlation plots", width = 32*0.3937, height = 16*0.3937)
par(mfrow = c(1, 2))

makeCorrPlot(paramValues_optimal_C_d13C[,seq(1,9)])
mtext(expression(paste("(A) Calibrated based on C and ", delta^{13}, "C")), side=3, line=0, cex=1.5, las=1, adj = 0.05, padj = -1)

makeCorrPlot(paramValues_optimal_C_d14C[,seq(1,9)])
mtext(expression(paste("(B) Calibrated based on C and ", Delta^{14}, "C")), side=3, line=0, cex=1.5, las=1, adj = 0.05, padj = -1)
par(mfrow = c(1, 1)) #To clear layout

# quartz.save("Figures for manuscript/Supplement_Figure_correlation_d13C_d14C.png", type="png")

dev.off()

# ------------------------------------------------------------------------------
# The correlation plots with data points
# ------------------------------------------------------------------------------

library(PerformanceAnalytics)
# chart.Correlation(paramValues_optimal_C[,seq(1,10)], histogram=TRUE, pch=19)

# --------------------------
# Some functions are defined
# --------------------------

panel.cor <- function(x, y, digits = 2, prefix = "", cex.cor, ...) {
  usr <- par("usr")
  on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r <- abs(cor(x, y, use = "complete.obs"))
  txt <- format(c(r, 0.123456789), digits = digits)[1]
  txt <- paste(prefix, txt, sep = "")
  if (missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex =  cex.cor * (1 + r) / 2)
}

panel.hist <- function(x, ...) {
  usr <- par("usr")
  on.exit(par(usr))
  par(usr = c(usr[1:2], 0, 1.5) )
  h <- hist(x, plot = FALSE)
  breaks <- h$breaks
  nB <- length(breaks)
  y <- h$counts
  y <- y/max(y)
  rect(breaks[-nB], 0, breaks[-1], y, col = "white", ...)
}

# ---------------------
# The plots are created
# ---------------------

quartz(title = "Correlation plots", width = 32*0.3937, height = 32*0.3937)
par(mfrow = c(2, 2))

# x <- paramValues_optimal_C[,seq(1,9)]
x <- paramValues_optimal_C_d13C_d14C[,seq(1,9)]

pairs(x,
      upper.panel = panel.cor,
      diag.panel  = panel.hist,
      lower.panel = panel.smooth)






my_fn <- function(data, mapping, ...){
  p <- ggplot(data = data, mapping = mapping) + 
    geom_point() + 
    # geom_smooth(method=loess, fill="red", color="red", ...) +
    geom_smooth(method=lm, fill="blue", color="blue", ...)
  p
}

ggpairs(x,lower = list(continuous = my_fn))





# -----------------------
# Let's try something new
# -----------------------

# IDs_ggpairs <- c("Db0", "VmaxD_POC", "VmaxD_poly", "Vmax_ads", "k_deprotect(0)", "k_mic_R", "k_mic_b", "beta_r")
IDs_ggpairs <- IDs

# Some changes are made for plotting purposes
allParam_C$optimal[allParam_C$optimal == 0] <- "All data"
allParam_C$optimal[allParam_C$optimal == 1] <- "Optimal data"
allParam_C_d13C$optimal[allParam_C_d13C$optimal == 0] <- "All data"
allParam_C_d13C$optimal[allParam_C_d13C$optimal == 1] <- "Optimal data"
allParam_C_d14C$optimal[allParam_C_d14C$optimal == 0] <- "All data"
allParam_C_d14C$optimal[allParam_C_d14C$optimal == 1] <- "Optimal data"
allParam_C_d13C_d14C$optimal[allParam_C_d13C_d14C$optimal == 0] <- "All data"
allParam_C_d13C_d14C$optimal[allParam_C_d13C_d14C$optimal == 1] <- "Optimal data"

# C only
quartz(title = "Correlation plots: C only", width = 38*0.3937, height = 32*0.3937)
ggpairs(data = allParam_C[,c(1:9)],
        columnLabels = IDs_ggpairs,
        aes(color = as.character(allParam_C$optimal)),
        lower = list(combo = wrap("facethist", bins = 30)),
        diag = list(continuous = wrap("densityDiag", alpha = 0.5)))

# C, d13C
quartz(title = "Correlation plots: C and d13C", width = 38*0.3937, height = 32*0.3937)
ggpairs(data = allParam_C_d13C[,c(1:8)],
        columnLabels = IDs_ggpairs,
        aes(color = as.character(allParam_C_d13C$optimal)),
        lower = list(combo = wrap("facethist", bins = 30)),
        diag = list(continuous = wrap("densityDiag", alpha = 0.5)))

# C, d14C
quartz(title = "Correlation plots: C and d14C", width = 38*0.3937, height = 32*0.3937)
ggpairs(data = allParam_C_d14C[,c(1:8)],
        columnLabels = IDs_ggpairs,
        aes(color = as.character(allParam_C_d14C$optimal)),
        lower = list(combo = wrap("facethist", bins = 30)),
        diag = list(continuous = wrap("densityDiag", alpha = 0.5)))
  
# C, d13C and d14C
quartz(title = "Correlation plots: C, d13C and d14C", width = 38*0.3937, height = 32*0.3937)
ggpairs(data = allParam_C_d13C_d14C[,c(1:8)],
        columnLabels = IDs_ggpairs,
        aes(color = as.character(allParam_C_d13C_d14C$optimal)),
        lower = list(combo = wrap("facethist", bins = 30)),
        diag = list(continuous = wrap("densityDiag", alpha = 0.5)))


# ------------------------------------------------------------------------------
# The turnover time of mic-r and mic-b is plotted
# ------------------------------------------------------------------------------

# ----------------------------------------------
# The data is loaded: depth profiles of microbes
# ----------------------------------------------

rm(list=ls(all = TRUE))
cat("\014")

dir_C_d13C_d14C <- "Analysis DE output/7. Long calibration/outputForSeries_C_d13C_d14C.RData"
load(dir_C_d13C_d14C)

C <- as.data.frame(t(outputForSeries$Ctot))
d13C <- as.data.frame(t(outputForSeries$d13C))
d14C <- as.data.frame(t(outputForSeries$d14C))
depth <- outputForSeries$midDepth
error <- as.data.frame(outputForSeries$error)
names(error) <- "error"

# The distribution of the errors is plotted
# ggplot(error, aes(error)) + 
#   geom_histogram(binwidth = 0.01) + 
#   xlim(c(0,5))

# # Results with a high error are removed
# errorTol_C_d13C_d14C <- 1.2
# n <- which(error > errorTol_C_d13C_d14C)

# The X % lowest errors are retained
q_C_d13C_d14C <- quantile(error$error,Qthreshold) # The error value which marks the Qthreshold
n <- which(error > q_C_d13C_d14C) # The rows with errors larger than this quantile

C <- subset(C,select = -n)
d13C <- subset(d13C,select = -n)
d14C <- subset(d14C,select = -n)
error_opt <- error[-n,1]

# The average C, d13C and d14C depth profiles are calculated
df_avg <- data.frame(depth = depth,
                     avg_C = rowMeans(C),
                     avg_d13C = rowMeans(d13C),
                     avg_d14C = rowMeans(d14C))

# The optimal profiles are stored
rowNum_opt <- which(error_opt == min(error_opt))
df_opt <- data.frame("depth" = depth,
                     "opt_C" = C[,rowNum_opt],
                     "opt_d13C" = d13C[,rowNum_opt],
                     "opt_d14C" = d14C[,rowNum_opt])

# The depth is added to the data frames
C$depth = depth
d13C$depth = depth
d14C$depth = depth

# The data frames are reformatted for plotting
df_C <- melt(C, id.vars = "depth")
df_d13C <- melt(d13C, id.vars = "depth")
df_d14C <- melt(d14C, id.vars = "depth")

# Plotting
p10 <- plot_C(df_C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p11 <- plot_d13C(df_d13C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)
p12 <- plot_d14C(df_d14C, df_avg, df_opt, measuredData_soil, color_lines = "deepskyblue3", color_avg = "dodgerblue4", par)

p10 <- p10 + annotate("label", x = 1.2, y = 0.56, label = expression(paste("Calibrated based on C and ", delta^{13}, "C and", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
p11 <- p11 + annotate("label", x = -27.9, y = 0.56, label = expression(paste("Calibrated based on C and ", delta^{13}, "C and", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')
p12 <- p12 + annotate("label", x = -430, y = 0.56, label = expression(paste("Calibrated based on C and ", delta^{13}, "C and", Delta^{14}, "C")), hjust = 0, size = 4, fontface = 'italic')

rm(list = c("df_C", "df_d13C", "df_d14C"))


















