Parameter calibration (P-Model) to multiple targets (GPP, bigD13C)

#    226: Setup c)      Delta^'13 C*',VJ
#    231: Setup h)      Delta^'13 C*',VJ, GPP

setup_rsofun_calibration <- function(scenario) {
  ## Load forcing and targets data from disk (here now from package) ----
  bigD13C_vj_gpp_drivers <- rsofun::pmodel_drivers
  bigD13C_vj_gpp_obs     <- rsofun::pmodel_validation

  # ## Read test-train split from disk ----
  # read_csv(here::here("data/01_test_train_split.csv"))
  # Not needed here for vignette

  ## Preprocess observation data (gpp) ----
  ## # no additinal QC needed

  ## Apply test-train split to data ----
  train_drivers <- bigD13C_vj_gpp_drivers
  train_obs     <- bigD13C_vj_gpp_obs

  test_drivers <- bigD13C_vj_gpp_drivers |> slice(0)
  test_obs     <- bigD13C_vj_gpp_obs     |> slice(0)

  ## Setup the settings for the different calibration scenarios ----
  ## Define default parameter
  default_par_fixed <- list( # fix parameter value from previous calibration
    kphio              = 0.04998,    # value from Stocker et al. 2020
    kphio_par_a        = 0.0,        # 0 corresponds to no temperature dependency of kphio (ORG setup in Stocker et al. 2020)
    kphio_par_b        = 1.0,        #
    soilm_thetastar    = 0.6 * 240,  # to recover paper setup with soil moisture stress
    soilm_betao        = 0.01,       # 1 corresponds to no reduction, 0 to full reduction at theta==0
    beta_unitcostratio = 146.0,      # value from Stocker et al. 2020
    rd_to_vcmax        = 0.014,      # value from Atkin et al. 2015 for C3 herbaceous
    tau_acclim         = 14.0,       # value from Liu et al. 2024
    kc_jmax            = 0.41        # value from Stocker et al. 2024 (citing Wang et al. 2017)
  )

  ## Define parameters to estimate and their priors
  if (scenario %in% c(223, 226, 231)) { # 231 is like 223 but using posteriors from 226
    par_to_estimate <- list(
      kphio           = list(lower = 0.02, upper = 0.15, init = 0.05),
      kphio_par_a     = list(lower = -0.004, upper = -0.001, init = -0.0025),
      kphio_par_b     = list(lower = 10, upper = 30, init = 20),
      soilm_thetastar = list(lower = 1, upper = 250, init = 40),
      beta_unitcostratio = as.list(c(lower = 0.1, upper = 3.0, init = 1.0) * 146.0),
      # truncated normal, with ~14 days as mean, taken from Liu et al.
      # 2024, Nat.Plants and Mäkelä et al. 2004, Tree Phys.
      tau_acclim      = list(mean = 14, sd = 8, lower = 0.01, upper = 40),
      kc_jmax         = as.list(0.41 * c(lower = 0.1, upper = 3.0, init = 1.0)),
      err_gpp         = list(lower = 0.01, upper = 3, init = 0.8),
      err_bigD13C     = list(lower = 0.01, upper = 3, init = 0.8)
    )
    if (scenario %in% c(223, 226, 231)) {
      par_to_estimate$err_bigD13C     <- list(lower = 0.01, upper = 15, init = 0.8)
    }
  } else {
    stop(sprintf("Unsupported scenario: %d", scenario))
  }

  if (scenario %in% c(231)) {   # 231 is like 223 but using posteriors from 226
    # for beta_unitcostratio and kc_jmax as priors
    # Below was done iteratively
    # calib_scen226 <- readr::read_rds(file.path(rsofun_doc_output_path, "data","calibrations","out_calib__scen226_DREAMzs-100000-0iter_8x3chains_on_CPU8x1_continued.rds"))
    # # i) extract samples as a data.frame
    # burnins_scen226 <- 30000
    # samples_scen226 <- getSample(calib_scen226$mod, thin = 1, start = burnins_scen226) %>% as.data.frame()
    #
    # # ii) fit normal and lognormal distributions for each parameter
    # param_normals_scen226 <- lapply(setNames(names(samples_scen226), names(samples_scen226)), function(p) {
    #   list(mean = mean(samples_scen226[[p]]),
    #       sd   = sd(  samples_scen226[[p]]))
    # })[c('beta_unitcostratio', 'kc_jmax')] # only keep these
    #
    # # then pass on these as prior for these
    # par_to_estimate$beta_unitcostratio <- param_normals_scen226$beta_unitcostratio
    # par_to_estimate$kc_jmax            <- param_normals_scen226$kc_jmax

    par_to_estimate$beta_unitcostratio <- list(mean    = 207.86, sd    = 6.79)
    par_to_estimate$kc_jmax            <- list(mean    = 0.4244, sd    = 0.0217)
    if (scenario %in% c(231)) {
      # truncated normal:
      par_to_estimate$beta_unitcostratio <- list(mean    = 207.86, sd    = 6.79,
        lower = 207.86 - 3 * 6.79,
        upper = 207.86 + 3 * 6.79)
      # normal:
      par_to_estimate$kc_jmax            <- list(mean    = 0.4244, sd    = 0.0217)
    }
  }
  # Remove parameters that are defined to be estimated from default_par_fixed
  par_to_fix <- default_par_fixed[!(names(default_par_fixed) %in% names(par_to_estimate))]

  ## Setup the data (drivers and obs) for the three calibration scenarios ----

  # Subset different combination sites to define target variables.
  # For easier handling do this in combined drivobs-object.
  drivobs_train_bigD13C_vj_gpp <- dplyr::inner_join(
    train_drivers,
    train_obs,
    by = join_by(sitename))

  drivobs_test_bigD13C_vj_gpp <- dplyr::inner_join(
    test_drivers,
    test_obs,
    by = join_by(sitename))

  all_potential_targets <- c("bigD13C" = "cornwell", "gpp" = "fluxnet", "le" = "fluxnet")
  if (scenario %in% c(223, 224, 225, 229, 230, 231)) { # GPP and traits data
    targets_to_keep <- c("gpp", "bigD13C")
  } else if (scenario %in% c(226)) {       # only traits data, either both, or vj only, or bigD13C only
    targets_to_keep <- c("bigD13C")
    # } else if (scenario %in% c(220)) {  # only GPP data from FR-Pue
    #   targets_to_keep <- c("gpp")
    #   sites_to_keep <- "FR-Pue"
  } else {
    stop(sprintf("Unsupported scenario: %d", scenario))
  }
  stopifnot(all(targets_to_keep %in% names(all_potential_targets)))

  targets_to_request <- all_potential_targets[targets_to_keep]
  
  drivobs_train <- drivobs_train_bigD13C_vj_gpp |>
    # keep only sites with targets we want
    dplyr::filter(source %in% unname(targets_to_request))

  drivobs_test <- drivobs_test_bigD13C_vj_gpp ## for the test data set keep all

  ## return ---
  return(list(
    drivobs_train = drivobs_train,
    drivobs_test = drivobs_test,
    targets = targets_to_request,
    par_fixed = par_to_fix,
    par = par_to_estimate
  ))
}


setup_rsofun_calibration(226)
#> $drivobs_train
#> # A tibble: 12 × 6
#>    sitename                params_siml      site_info forcing  source   data    
#>    <chr>                   <list>           <list>    <list>   <chr>    <list>  
#>  1 lon_+010.52_lat_+051.08 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  2 lon_+011.10_lat_+048.30 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  3 lon_+112.58_lat_+023.13 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  4 lon_+145.13_lat_-005.83 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  5 lon_+146.13_lat_-032.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  6 lon_+148.30_lat_-036.10 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  7 lon_+153.00_lat_-026.85 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  8 lon_-079.10_lat_+035.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  9 lon_-116.45_lat_+047.16 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 10 lon_-119.82_lat_+034.50 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 11 lon_-122.98_lat_+038.40 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 12 lon_-149.61_lat_+063.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 
#> $drivobs_test
#> # A tibble: 0 × 6
#> # ℹ 6 variables: sitename <chr>, params_siml <list>, site_info <list>,
#> #   forcing <list>, source <chr>, data <list>
#> 
#> $targets
#>    bigD13C 
#> "cornwell" 
#> 
#> $par_fixed
#> $par_fixed$soilm_betao
#> [1] 0.01
#> 
#> $par_fixed$rd_to_vcmax
#> [1] 0.014
#> 
#> 
#> $par
#> $par$kphio
#> $par$kphio$lower
#> [1] 0.02
#> 
#> $par$kphio$upper
#> [1] 0.15
#> 
#> $par$kphio$init
#> [1] 0.05
#> 
#> 
#> $par$kphio_par_a
#> $par$kphio_par_a$lower
#> [1] -0.004
#> 
#> $par$kphio_par_a$upper
#> [1] -0.001
#> 
#> $par$kphio_par_a$init
#> [1] -0.0025
#> 
#> 
#> $par$kphio_par_b
#> $par$kphio_par_b$lower
#> [1] 10
#> 
#> $par$kphio_par_b$upper
#> [1] 30
#> 
#> $par$kphio_par_b$init
#> [1] 20
#> 
#> 
#> $par$soilm_thetastar
#> $par$soilm_thetastar$lower
#> [1] 1
#> 
#> $par$soilm_thetastar$upper
#> [1] 250
#> 
#> $par$soilm_thetastar$init
#> [1] 40
#> 
#> 
#> $par$beta_unitcostratio
#> $par$beta_unitcostratio$lower
#> [1] 14.6
#> 
#> $par$beta_unitcostratio$upper
#> [1] 438
#> 
#> $par$beta_unitcostratio$init
#> [1] 146
#> 
#> 
#> $par$tau_acclim
#> $par$tau_acclim$mean
#> [1] 14
#> 
#> $par$tau_acclim$sd
#> [1] 8
#> 
#> $par$tau_acclim$lower
#> [1] 0.01
#> 
#> $par$tau_acclim$upper
#> [1] 40
#> 
#> 
#> $par$kc_jmax
#> $par$kc_jmax$lower
#> [1] 0.041
#> 
#> $par$kc_jmax$upper
#> [1] 1.23
#> 
#> $par$kc_jmax$init
#> [1] 0.41
#> 
#> 
#> $par$err_gpp
#> $par$err_gpp$lower
#> [1] 0.01
#> 
#> $par$err_gpp$upper
#> [1] 3
#> 
#> $par$err_gpp$init
#> [1] 0.8
#> 
#> 
#> $par$err_bigD13C
#> $par$err_bigD13C$lower
#> [1] 0.01
#> 
#> $par$err_bigD13C$upper
#> [1] 15
#> 
#> $par$err_bigD13C$init
#> [1] 0.8
setup_rsofun_calibration(226)$drivobs_train
#> # A tibble: 12 × 6
#>    sitename                params_siml      site_info forcing  source   data    
#>    <chr>                   <list>           <list>    <list>   <chr>    <list>  
#>  1 lon_+010.52_lat_+051.08 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  2 lon_+011.10_lat_+048.30 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  3 lon_+112.58_lat_+023.13 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  4 lon_+145.13_lat_-005.83 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  5 lon_+146.13_lat_-032.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  6 lon_+148.30_lat_-036.10 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  7 lon_+153.00_lat_-026.85 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  8 lon_-079.10_lat_+035.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#>  9 lon_-116.45_lat_+047.16 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 10 lon_-119.82_lat_+034.50 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 11 lon_-122.98_lat_+038.40 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
#> 12 lon_-149.61_lat_+063.97 <tibble [1 × 2]> <tibble>  <tibble> cornwell <tibble>
setup_rsofun_calibration(231)
#> $drivobs_train
#> # A tibble: 20 × 6
#>    sitename                params_siml       site_info forcing  source  data    
#>    <chr>                   <list>            <list>    <list>   <chr>   <list>  
#>  1 CH-Dav                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  2 CZ-BK1                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  3 DK-Sor                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  4 FI-Hyy                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  5 FR-Pue                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  6 GF-Guy                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  7 US-Ha1                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  8 US-PFa                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  9 lon_+010.52_lat_+051.08 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 10 lon_+011.10_lat_+048.30 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 11 lon_+112.58_lat_+023.13 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 12 lon_+145.13_lat_-005.83 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 13 lon_+146.13_lat_-032.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 14 lon_+148.30_lat_-036.10 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 15 lon_+153.00_lat_-026.85 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 16 lon_-079.10_lat_+035.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 17 lon_-116.45_lat_+047.16 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 18 lon_-119.82_lat_+034.50 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 19 lon_-122.98_lat_+038.40 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 20 lon_-149.61_lat_+063.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 
#> $drivobs_test
#> # A tibble: 0 × 6
#> # ℹ 6 variables: sitename <chr>, params_siml <list>, site_info <list>,
#> #   forcing <list>, source <chr>, data <list>
#> 
#> $targets
#>        gpp    bigD13C 
#>  "fluxnet" "cornwell" 
#> 
#> $par_fixed
#> $par_fixed$soilm_betao
#> [1] 0.01
#> 
#> $par_fixed$rd_to_vcmax
#> [1] 0.014
#> 
#> 
#> $par
#> $par$kphio
#> $par$kphio$lower
#> [1] 0.02
#> 
#> $par$kphio$upper
#> [1] 0.15
#> 
#> $par$kphio$init
#> [1] 0.05
#> 
#> 
#> $par$kphio_par_a
#> $par$kphio_par_a$lower
#> [1] -0.004
#> 
#> $par$kphio_par_a$upper
#> [1] -0.001
#> 
#> $par$kphio_par_a$init
#> [1] -0.0025
#> 
#> 
#> $par$kphio_par_b
#> $par$kphio_par_b$lower
#> [1] 10
#> 
#> $par$kphio_par_b$upper
#> [1] 30
#> 
#> $par$kphio_par_b$init
#> [1] 20
#> 
#> 
#> $par$soilm_thetastar
#> $par$soilm_thetastar$lower
#> [1] 1
#> 
#> $par$soilm_thetastar$upper
#> [1] 250
#> 
#> $par$soilm_thetastar$init
#> [1] 40
#> 
#> 
#> $par$beta_unitcostratio
#> $par$beta_unitcostratio$mean
#> [1] 207.86
#> 
#> $par$beta_unitcostratio$sd
#> [1] 6.79
#> 
#> $par$beta_unitcostratio$lower
#> [1] 187.49
#> 
#> $par$beta_unitcostratio$upper
#> [1] 228.23
#> 
#> 
#> $par$tau_acclim
#> $par$tau_acclim$mean
#> [1] 14
#> 
#> $par$tau_acclim$sd
#> [1] 8
#> 
#> $par$tau_acclim$lower
#> [1] 0.01
#> 
#> $par$tau_acclim$upper
#> [1] 40
#> 
#> 
#> $par$kc_jmax
#> $par$kc_jmax$mean
#> [1] 0.4244
#> 
#> $par$kc_jmax$sd
#> [1] 0.0217
#> 
#> 
#> $par$err_gpp
#> $par$err_gpp$lower
#> [1] 0.01
#> 
#> $par$err_gpp$upper
#> [1] 3
#> 
#> $par$err_gpp$init
#> [1] 0.8
#> 
#> 
#> $par$err_bigD13C
#> $par$err_bigD13C$lower
#> [1] 0.01
#> 
#> $par$err_bigD13C$upper
#> [1] 15
#> 
#> $par$err_bigD13C$init
#> [1] 0.8
setup_rsofun_calibration(231)$drivobs_train
#> # A tibble: 20 × 6
#>    sitename                params_siml       site_info forcing  source  data    
#>    <chr>                   <list>            <list>    <list>   <chr>   <list>  
#>  1 CH-Dav                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  2 CZ-BK1                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  3 DK-Sor                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  4 FI-Hyy                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  5 FR-Pue                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  6 GF-Guy                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  7 US-Ha1                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  8 US-PFa                  <tibble [1 × 12]> <tibble>  <tibble> fluxnet <tibble>
#>  9 lon_+010.52_lat_+051.08 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 10 lon_+011.10_lat_+048.30 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 11 lon_+112.58_lat_+023.13 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 12 lon_+145.13_lat_-005.83 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 13 lon_+146.13_lat_-032.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 14 lon_+148.30_lat_-036.10 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 15 lon_+153.00_lat_-026.85 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 16 lon_-079.10_lat_+035.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 17 lon_-116.45_lat_+047.16 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 18 lon_-119.82_lat_+034.50 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 19 lon_-122.98_lat_+038.40 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>
#> 20 lon_-149.61_lat_+063.97 <tibble [1 × 2]>  <tibble>  <tibble> cornwe… <tibble>

run_mcmc_rsofun <- function(
    curr_calibration_scenario,
    # MCMC setup:
    iterations = 6,
    burnin = 0,
    n_chains_independent    = 3, # number of independent chains (called 'nrChains' in runMCMC)
    n_chains_within_sampler = 3, # number of internal chains to be sampled (called 'startValue' in runMCMC, at least 2 for DEzs)

    # parallelization:
    n_parallel_independent  = 3, # number of cores for parallelization of independent chains     https://cran.r-project.org/web/packages/BayesianTools/vignettes/InterfacingAModel.html#running-several-mcmcs-in-parallel
    outpath = "data",
    suffix_str = "",
    sampler = "DREAMzs"
    ) {
  # Setup simulation model
  res <- setup_rsofun_calibration(scenario = curr_calibration_scenario)
  # res$drivobs_train
  # res$drivobs_test
  # res$par_fixed
  # res$par
  # res$targets

  # Load loglikelihood (e.g. if using a tailor-made one, instead of package-provided):
  # source("R/calibration_helpers.R", echo = FALSE)
  # source("R/cost_likelihood_pmodel.R", echo = FALSE)

  # Setup MCMC          # uses drivobs_train, not drivobs_test
  driver_to_use_for_mcmc <- select(res$drivobs_train, sitename, params_siml, site_info, forcing)
  obs_to_use_for_mcmc    <- select(res$drivobs_train, sitename, source, data)

  calib_sofun_settings <- list(
    method = "BayesianTools",
    metric = rsofun::cost_likelihood_pmodel,
    control = list(
      sampler_runMCMC = sampler,
      settings_runMCMC = list(
        burnin     = burnin,                 # 10000,
        iterations = iterations,             # 50000,
        nrChains   = n_chains_independent,   # number of independent chains to be sampled
        startValue = n_chains_within_sampler # number of internal chains to be sampled
      ),
      n_parallel_nrChains = n_parallel_independent
    ),
    par = res$par
  )

  # setup output directories for parallel calibration
  logpath <- file.path(outpath, "calibrations", paste0("out_calib_", suffix_str, ".rds.log.txt"))
  rds_path <- gsub(".log.txt", "", logpath)
  dir.create(path = dirname(logpath), showWarnings = FALSE)

  # Run calibration in parallel
  cat(paste0(Sys.time(), ": start sampling of ", paste0("out_calib_", suffix_str, ".rds")))
  out_calib <- calib_sofun(
    drivers   = driver_to_use_for_mcmc,
    obs       = obs_to_use_for_mcmc,
    settings_calib = calib_sofun_settings,
    logpath   = logpath,
    # other arguments for the cost function
    targets = res$targets,
    par_fixed = res$par_fixed
  )
  ## Store result object to file: ----
  ## This enables easy continuation of MCMC sampling
  saveRDS(out_calib, file = rds_path, compress = "xz")
  cat(paste0(Sys.time(), ": end   sampling of ", out_calib$name,
    ". \nWritten *.rds-output to: ", rds_path))

  # return only performance results as tibble
  timings <- tibble(
    setup       = curr_calibration_scenario,
    name        = gsub(".rds.log.txt", "", out_calib$name),
    resultfile  = rds_path,
    # sampling options:
    sampler        = calib_sofun_settings$control$sampler_runMCMC,
    burnin         = burnin,
    iterations     = iterations,
    n_chains       = n_chains_independent,
    n_chains_inner = n_chains_within_sampler,
    # performance results:
    cores          = n_parallel_independent,
    walltime       = out_calib$walltime,
    # logfile        = out_calib$logpath
  )
  return(timings)
}

dir.create("vignettes/files/param_calib_multitarget/", showWarnings = FALSE)

df_timings226 <- run_mcmc_rsofun(226, suffix_str = "setup_226",
  outpath = "vignettes/files/param_calib_multitarget/",
  # mcmc setup:
  sampler = "DREAMzs",
  iterations = 600,
  burnin = 0,
  n_chains_independent      = 8,
  n_chains_within_sampler   = 3,
  # parallelization
  n_parallel_independent    = 1
)
# 2026-07-09 10:39:19.622351: start sampling of out_calib_setup_226.rds
# runMCMC terminated after 4.39500000000407seconds00 . Current logp  -389.1087 -421.1468 -396.6275 Please wait!
# runMCMC terminated after 3.95499999998719seconds00 . Current logp  -377.7371 -379.5134 -377.3347 Please wait!
# ...
# runMCMC terminated after 4.23500000000058seconds00 . Current logp  -397.9886 -379.6501 -382.6354 Please wait!

df_timings231 <- run_mcmc_rsofun(231, suffix_str = "setup_231",
  outpath = "vignettes/files/param_calib_multitarget/",
  # mcmc setup:
  sampler = "DREAMzs",
  iterations = 600,
  burnin = 0,
  n_chains_independent      = 8,
  n_chains_within_sampler   = 3,
  # parallelization
  n_parallel_independent    = 1
)
# 2026-07-09 10:39:53.032407: start sampling of out_calib_setup_231.rds
# Running DREAM-MCMC, chain  1 iteration 300 of 600 . Current logp  -14645.97 -14571.83 -15471.88 Please wait!
# runMCMC terminated after 16.8090000000084seconds00 . Current logp  -14616.55 -14462.62 -14529.8 Please wait!
# Running DREAM-MCMC, chain  2 iteration 300 of 600 . Current logp  -16760.01 -15828.73 -14948.44 Please wait!
# runMCMC terminated after 16.4020000000019seconds00 . Current logp  -14955.26 -14941.47 -14441.98 Please wait!
# ...
# runMCMC terminated after 17.8410000000003seconds00 . Current logp  -15023.04 -14854.75 -14265.61 Please wait!

df_timings231p <- run_mcmc_rsofun(231, suffix_str = "setup_231_parallel",
  outpath = "vignettes/files/param_calib_multitarget/",
  # mcmc setup:
  sampler = "DREAMzs",
  iterations = 600,
  burnin = 0,
  n_chains_independent      = 8,
  n_chains_within_sampler   = 3,
  # parallelization
  n_parallel_independent    = 8     # here now we parallelize the 8 chains to 8 cores
)
# 2026-07-09 10:42:24.395628: start sampling of out_calib_setup_231_parallel.rds
# Writing MCMC sampling log to: files/param_calib_multitarget//calibrations/out_calib_setup_231_parallel.rds.log.txt
# ...
# 2026-07-09 10:43:07.689861: end   sampling of out_calib_setup_231_parallel.rds.log.txt.
# Written *.rds-output to: files/param_calib_multitarget//calibrations/out_calib_setup_231_parallel.rds

# calib_timings <- bind_rows(
#   df_timings226,
#   df_timings231,
#   df_timings231p)
knitr::kable(calib_timings |>
  dplyr::mutate(name = gsub(".rds.log.txt", "", name)) |>
  dplyr::mutate(resultfile = gsub("vignettes/files/param_calib_multitarget//", "", resultfile)) |>
  dplyr::select(-"resultfile", "resultfile")
)

# This loads MCMC chains that might have been generated with a separate R script
out_calib_226 <- readRDS("files/param_calib_multitarget/calibrations/out_calib_setup_226.rds")
out_calib_231p <- readRDS("files/param_calib_multitarget/calibrations/out_calib_setup_231_parallel.rds")

rsofun_symbol_parname_description <- dplyr::tribble(
  ~variable,            ~label,                                                      ~Symbol_R,
  "kphio",              "italic(phi)[0]^'*'*' (mol mol-\u00B9)'",               "italic(phi)[0]^'*'",
  "kphio_par_a",        "italic(a)[italic(phi)]*' (\u00B0C-\u00B2)'",           "italic(a)[italic(phi)]",
  "kphio_par_b",        "italic(b)[italic(phi)]*' (\u00B0C)'",                       "italic(b)[italic(phi)]",
  "soilm_thetastar",    "italic(theta)^'*'*' ('*'mm'*')'",                           "italic(theta)^'*'",
  "soilm_betao",        "italic(beta)[0]*' (-)'",                                    "italic(beta)[0]",
  "beta_unitcostratio", "italic(beta)*' (-)'",                                       "italic(beta)",
  "rd_to_vcmax",        "italic(b)[0]*' (-)'",                                       "italic(b)[0]",
  "tau_acclim",         "italic(tau)*' (days)'",                                     "italic(tau)",
  "kc_jmax",            "italic(c)^'*'*' (-)'",                                      "italic(c)^'*'",
  "err_gpp",            "italic(sigma)['GPP']*' (gC m -\u00B2 s-\u00B9)'", "italic(sigma)['GPP']",
  "err_bigD13C",        "italic(sigma)[Delta]*' ('*'\u2030'*')'",                    "italic(sigma)[Delta]"
)
label_vec_short <- setNames(rsofun_symbol_parname_description$Symbol_R, rsofun_symbol_parname_description$variable)
custom_labeller_variable <- function(labels, multi_line = TRUE) { # adapted from label_parsed
  replaced_labels <- left_join(
    select(labels, variable), # NOTE: this has variable hardcoded
    select(rsofun_symbol_parname_description, variable, label),
    by = join_by(variable))
  # print(tibble(replaced_labels))
  replaced_labels <- replaced_labels |> select(-variable)
  replaced_labels <- label_value(replaced_labels, multi_line = multi_line)
  # print(replaced_labels)
  lapply(unname(replaced_labels), lapply, function(values) {
    c(parse(text = as.character(values)))
  })
}

# Define plotting function
plot_mcmc_trace <- function(out_calib, nr_internal_chains, burnin_to_skip, burnin_to_skip_gelman = burnin_to_skip, dont_thin = FALSE, end = NULL) {
  x <- out_calib$mod
  title <- basename(out_calib$name)

  curr_iter <- x[[1]]$settings$iterations
  if (dont_thin || curr_iter < 10000) {
    curr_thin <- 1
  } else {
    curr_thin <- floor(curr_iter / 10000)
  }

  xsample <- BayesianTools::getSample(x, coda = T, thin = curr_thin, start = burnin_to_skip, end = end)

  # nr_internal_chains will have same color
  dat_to_plot <- lapply(xsample, function(single_chain) {
    as_tibble(single_chain) |> mutate(iteration = burnin_to_skip + curr_thin * (1:n()))
  }) |>
    dplyr::bind_rows(.id = "chain_id") |>
    tidyr::pivot_longer(-c(iteration, chain_id), names_to = "variable") |>
    # mark inner and outer chains (assumes DEzs or DREAMzs):
    mutate(outerChain = as.factor(ceiling(as.numeric(chain_id) / 3)),
      innerChain = (as.numeric(chain_id) + 2) %% 3 + 1,
      innerChain_str = letters[innerChain],
      chain_id_str = paste0(outerChain, letters[innerChain])) #|>
  # fix order: in order of appearance
  # mutate(variable = forcats::as_factor(variable))

  pl <- ggplot(dat_to_plot,
    aes(x = iteration, y = value, color = outerChain, linetype = innerChain_str)) + geom_line() +
    # geom_rug(sides = "r") +
    theme_classic() +
    # scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
    facet_wrap(~variable,  nrow = 2, scales = "free_y", labeller = custom_labeller_variable) +
    theme(
      legend.position = "bottom", strip.text = element_text(size = 12)
    ) +
    labs(y = "", color = "chain", linetype = "internal\nchains")

  # add Gelman Diagnostics
  get_gelman_diag <- function(mcmc, burnin_to_skip, end) {
    gelman_df <- BayesianTools::gelmanDiagnostics(mcmc, start = burnin_to_skip, end = end)
    psrf_values <- gelman_df$psrf[, "Point est."]
    # psrf_strings <- paste0(substr(names(psrf_values),1,5), "..=", sprintf("%.2f", psrf_values))
    psrf_strings <- paste0("'*", unname(label_vec_short[names(psrf_values)]), "*'=", sprintf("%.2f", psrf_values))
    psrf_string <- paste0(psrf_strings, collapse = ", ")

    subtitle <- sprintf("'GelmanDiagnostics: mpsrf=%.1f; psrf:%s'",
      gelman_df$mpsrf,
      psrf_string)
  }
  subtitle <- tryCatch(get_gelman_diag(x, burnin_to_skip_gelman + 1, end = end), error = function(e) {
    e
  }) # unsure why min burnin of 1 is needed
  pl <- pl + ggtitle(title, subtitle = parse(text = subtitle))

  pl <- pl + geom_vline(xintercept = burnin_to_skip_gelman, color = "red", linetype = "dashed")

  return(pl)
}

plot_mcmc_trace(out_calib_226, nr_internal_chains = 3, burnin_to_skip = 0, burnin_to_skip_gelman = 150)

plot_mcmc_trace(out_calib_231p, nr_internal_chains = 3, burnin_to_skip = 0, burnin_to_skip_gelman = 150)

plot_prior_posterior_density <- function(x, burnin_to_skip) {
  require(BayesianTools)
  require(dplyr)
  require(tidyr)
  require(ggplot2)

  # Get matrices of prior and posterior samples
  posteriorMat <- BayesianTools::getSample(x, parametersOnly = TRUE, start = burnin_to_skip)
  priorMat <-  BayesianTools:::getSetup(x)$prior$sampler(10000) # nPriorDraws = 10000

  # Parameter names
  parNames <- colnames(posteriorMat)
  # rename columns priorMat
  colnames(priorMat) <- parNames

  # Create data frame for plotting
  df_plot <- rbind(
    data.frame(priorMat,     par_estimation = "prior"),
    data.frame(posteriorMat, par_estimation = "posterior")
  ) |>
    # dplyr::mutate(par_estimation = forcats::fct_inorder(par_estimation)) |> # order by appearance
    tibble()
  levels(df_plot$par_estimation)
  # Plot with facet wrap
  gg <- df_plot |>
    tidyr::pivot_longer(-c(par_estimation), names_to = "variable") |>
    # dplyr::mutate(variable = forcats::fct_inorder(variable)) |> # order by appearance
    ggplot(
      aes(x = value, fill = par_estimation)
    ) +
    geom_density(alpha = 0.8) +
    theme_classic() +
    facet_wrap(~variable,  nrow = 2, scales = "free", labeller = custom_labeller_variable) +
    # facet_wrap( ~ variable , nrow = 2, scales = "free") +
    theme(
      legend.position = "bottom",
      axis.title.x = element_text(""),
      axis.ticks.y = element_blank(),
      axis.text.y = element_blank(),
    ) +
    scale_fill_manual(NULL, values = c("posterior" = "#29a274ff",
      "prior" = "#777055ff")) # GECO colors

  return(gg)
}

library(cowplot)
p1 <- plot_prior_posterior_density(out_calib_226$mod,  burnin_to_skip = 200) +  ggtitle(out_calib_226$name)
#> Loading required package: BayesianTools
#> Loading required package: tidyr
p2 <- plot_prior_posterior_density(out_calib_231p$mod, burnin_to_skip = 20) +  ggtitle(out_calib_231p$name)
cowplot::plot_grid(p1, p2, ncol = 1)

BayesianTools::correlationPlot(out_calib_231p$mod)

# Samples posterior parameters (MAP or N=n_samples) and runs
# model predictions
#
# This allows to do two things:
# A) make predictions at sites/dates/... where we have observations
# B) make predictions also  on dates/... where we have no observations, e.g. for continuous time series
# we do B) once {by setting `return_continuous_timeseries = TRUE`} and then derive A) from it.
run_prediction_rsofun <- function(
    mcmc_posterior,
    prediction = c("both", "test", "train"),
    burnin_to_skip = 0,
    n_samples = 100,      # if n_samples == 1, use MAP
    n_cores = NULL) {

  if (length(prediction) == 1 && (prediction %in% c("both", "test", "train"))) {
    # as expected
  } else {
    stop("Provide prediction as either: 'both', 'test' or 'train'")
  }
  n_cores <- ifelse(is.null(n_cores), 1, n_cores)
  stopifnot(is(mcmc_posterior$mod, "mcmcSamplerList"))

  curr_calibration_scenario <- as.integer(gsub(
    ".*_setup_([0-9]*)_.*",  # NOTE: hardcoded expected format of filename
    "\\1",
    mcmc_posterior$name))
  stopifnot(!is.na(curr_calibration_scenario)) # Catches if expected format is inaccurate

  # Setup simulation model
  res <- setup_rsofun_calibration(scenario = curr_calibration_scenario)

  # Set random seed for reproducibility
  set.seed(2023)

  # Sample parameters from MCMC posterior
  # Evaluation of the uncertainty coming from the model parameters' uncertainty
  if (n_samples > 1) {
    samples_par <- getSample(
      mcmc_posterior$mod,
      thin = 1,
      start = burnin_to_skip, numSamples = n_samples
    ) |>
      as.data.frame() |>
      # Add sample IDs
      dplyr::mutate(posterior_sample_id = 1:n()) |>
      tidyr::nest(.by = posterior_sample_id, .key = "pars")
  } else {
    # mcmc_posterior$par # these are already precomputed...
    # but more robust to recompute:
    samples_par <- BayesianTools::MAP(mcmc_posterior$mod)$parametersMAP |>
      as.list() |> as_tibble() |>
      # Add sample IDs
      dplyr::mutate(posterior_sample_id = 0L) |> # posterior_sample_id == 0 means MAP
      tidyr::nest(.by = posterior_sample_id, .key = "pars")
  }

  # Setup prediction
  predict_sofun_settings <- list(n_cores = n_cores)

  if (prediction == "both") {
    curr_driver <- bind_rows(res$drivobs_train, res$drivobs_test) |>
      select(sitename, params_siml, site_info, forcing)
    curr_obs    <- bind_rows(res$drivobs_train, res$drivobs_test) |>
      select(sitename, source, data)
  } else if (prediction == "train") {
    curr_driver <- select(res$drivobs_train, sitename, params_siml, site_info, forcing)
    curr_obs    <- select(res$drivobs_train, sitename, source, data)
  } else if (prediction == "test") {
    curr_driver <- select(res$drivobs_test, sitename, params_siml, site_info, forcing)
    curr_obs    <- select(res$drivobs_test, sitename, source, data)
    stopifnot(nrow(curr_obs) > 0)
  }

  # Run prediction
  df_pred_vs_obs <- predict_sofun_parallelized(
    drivers     = curr_driver,
    obs         = curr_obs,
    settings    = predict_sofun_settings,
    par         = samples_par,
    par_fixed   = res$par_fixed
  )

  return(df_pred_vs_obs)
}

# Runs the requested model predictions and returns results
# Requested model predictions are defined by data.frame `par_df` containing
# model parameter sets
predict_sofun_parallelized <- function(
    drivers,
    obs,
    settings,
    targets,
    par_fixed,
    par_df = NULL) {
  # Set number of cores if not specified
  if (is.null(settings$n_cores)) {
    settings$n_cores <- min(detectCores() - 1, 20) # at most 20
  }
  # ensure not more than needed
  settings$n_cores <- min(nrow(par_df), settings$n_cores)

  # Function to run prediction for a single parameter set
  run_pmodel_single_prediction <- function(par, par_fixed, targets, drivers, obs) {
    # Function that runs the P-model for a sample of parameters
    # but does not add the observation error

    # Taken from cost_likelihood_pmodel()
    stopifnot(nrow(obs) > 0)     # ensure some observation data are provided
    stopifnot(nrow(drivers) > 0) # ensure some driver data are provided

    # A) Include current parameters ----
    stopifnot(length(intersect(names(par), names(par_fixed))) == 0) # no overlap
    params_modl <- c(par, par_fixed)
    
    # B1) append column targets to obs data.frame: (code copy from cost_likelihood)
    #   This ensures that variant that requires targets as function argument can be 
    #   run with the code below (which was written with targets as a column of obs defined as:
    #       \item{targets}{A single one-row tibble with a column for each target
    #                      containing TRUE or FALSE (or NA_logical)}
    targets = c( 
      "gpp"     = "fluxnet",   # NOTE: hardcoding TODO: remove
      "le"      = "fluxnet",   # NOTE: hardcoding TODO: remove
      "bigD13C" = "cornwell",  # NOTE: hardcoding TODO: remove
      "vj"      = "cornwell")  # NOTE: hardcoding TODO: remove
    targets_df <- tidyr::tibble(target = names(targets), source = unname(targets)) |> # or use enframe(..., name = "target", value = "source")
      dplyr::group_by(source) |> dplyr::summarise(targets = unique(list(target)))
    obs <- dplyr::left_join(obs, targets_df, by = dplyr::join_by("source")) |> dplyr::select(-"source")
  
    # B,C) Run model and bring together with observed ----
    df_pred_vs_obs <- rsofun:::get_mod_obs_pmodel(
      drivers,
      obs,
      params_modl,
      parallel = FALSE,
      ncores = 1,
      return_continuous_timeseries = TRUE)

    # D) (DON'T) Sample error model ----
    # NOTE: sampling is not done here, but can optionally be done before plotting
    #       Here we rename to clarify that no error model has yet been applied.
    df_pred_vs_obs <- df_pred_vs_obs |>
      # clarify name of model output (containing not yet any error model term)
      rename(mod_no_err = mod) |>
      relocate(c(mod_no_err, err_par_sd), .after = last_col())

    return(df_pred_vs_obs)
  }

  # Run the P-model predictions for each set of parameters
  ## Prepare parallelization:
  parallel_flag <- (settings$n_cores > 1 && nrow(par_df) > 1)
  # setup helpers for conditional parallelization
  do_if <- function(df, cond, f) { # enable conditional lines in dplyr piping
    if (cond) f(df) else df
  }
  cl <- if (parallel_flag) {
    multidplyr::new_cluster(settings$n_cores) |>
      multidplyr::cluster_library(
        packages = c("dplyr", "tidyr", "purrr", "rsofun")
      )
  }
  ## Run predictions:
  df_model_predictions <- par_df |>
    do_if(parallel_flag, function(df) multidplyr::partition(df, cl)) %>%
    dplyr::mutate("sim" = purrr::map(
      pars,
      ~ run_pmodel_single_prediction(
        par = .x,
        par_fixed,
        targets,
        drivers,
        obs
      )
    )) |>
    do_if(parallel_flag, function(df) dplyr::collect(df))

  return(df_model_predictions)
}

apply_bias_correction_and_sample_error <- function(df_pred, N_sample_error) {
  set.seed(1982)
  df_pred |>
    # lazy_dt() |>   # if needed, uses lazy data.table and dtplyr for speed
    # a) repeat lines: once for each sampled error
    # following line is basically a cross_join: with
    slice(rep(1:n(), each = N_sample_error)) |>
    mutate(error_sample_id = as.integer(rep(1:N_sample_error, n() / N_sample_error))) |>
    arrange(posterior_sample_id) |>
    # b) sample the error:
    group_by(err_par_sd, err_par_bias) |>
    mutate(err_sample                = rnorm(n(),   sd = err_par_sd),
      mod_biasremoved_no_err    = mod_no_err - err_par_bias,
      mod_biasremoved_with_err  = mod_no_err - err_par_bias + err_sample) |>
    ungroup() |>
    as_tibble() |> # to access results of lazy dtplyr-computation as normal tibble()
    # keep output light: i.e. remove unneded columns:
    select(-err_sample) |>
    mutate(across(where(is.character), as.factor))
}

# get MAP and run model for MAP parameter set
df_predict_MAP_train_231p <- run_prediction_rsofun(
  mcmc_posterior = out_calib_231p,
  prediction     = "train",
  burnin_to_skip = 200,
  n_samples      = 1, # n_samples == 1, requests MAP
  n_cores        = 1) # this could be increased e.g. to 4

# get parameter samples and run model for thes parameter sets
df_predict_train_231p <- run_prediction_rsofun(
  mcmc_posterior = out_calib_231p,
  prediction     = "train",
  burnin_to_skip = 200,
  n_samples      = 20,
  n_cores        = 1) # this could be increased e.g. to 4

df_predict <- bind_rows(
  # df_predict_test      |> mutate(is_train0_test1 = 1L, is_MAP = FALSE),
  df_predict_train_231p     |> mutate(is_train0_test1 = 0L, is_MAP = FALSE),
  # df_predict_MAP_test  |> mutate(is_train0_test1 = 1L, is_MAP = TRUE),
  df_predict_MAP_train_231p |> mutate(is_train0_test1 = 0L, is_MAP = TRUE)
)

# Free memory:
# rm(df_predict_train_231p,
#    df_predict_MAP_train_231p)

# df_predict_params<- df_predict |> select(posterior_sample_id, is_train0_test1, is_MAP, pars) |> unnest(pars)
df_predict_gpp     <- df_predict |> select(posterior_sample_id, is_train0_test1, is_MAP, sim)  |> unnest(sim) |> filter(target == "gpp")     |> unnest(obs_metadata)
df_predict_bigD13C <- df_predict |> select(posterior_sample_id, is_train0_test1, is_MAP, sim)  |> unnest(sim) |> filter(target == "bigD13C") |> unnest(obs_metadata)

# apply bias-correction and simulate structural error
N_sample_error <- 3
df_predict_231p_gpp_sampled <- apply_bias_correction_and_sample_error(
  df_predict_gpp,
  N_sample_error = N_sample_error)
df_predict_231p_bigD13C_sampled <- apply_bias_correction_and_sample_error(
  df_predict_bigD13C,
  N_sample_error = N_sample_error)

# Save predictions to plot with separate script:
# FOR STORAGE REASONS unnest the obs_metadata and save different targets separately:
out_prediction_template <- file.path("vignettes/files/param_calib_multitarget/",
                                      "predictions", "pred_231p_XXX.rds")
dir.create(dirname(out_prediction_template), showWarnings = FALSE)
fname_out_gpp     <- gsub("_XXX",
                          paste0(    "_gpp_sampled_N", N_sample_error, "errors"),
                          out_prediction_template)
fname_out_bigD13C <- gsub("_XXX",
                          paste0("_bigD13C_sampled_N", N_sample_error, "errors"),
                          out_prediction_template)
df_predict_231p_gpp_sampled     |>
  select(posterior_sample_id,error_sample_id,is_train0_test1,is_MAP,sitename,target,obs,   date,err_par_bias,err_par_sd,mod_no_err,mod_biasremoved_no_err,mod_biasremoved_with_err) |>
  saveRDS(fname_out_gpp,     compress = "xz")
df_predict_231p_bigD13C_sampled |>
  select(posterior_sample_id,error_sample_id,is_train0_test1,is_MAP,sitename,target,obs,id,     err_par_bias,err_par_sd,mod_no_err,mod_biasremoved_no_err,mod_biasremoved_with_err) |>
  saveRDS(fname_out_bigD13C, compress = "xz")

# Load predictions for plotting with separate script:
df_predict_231p_gpp_sampled     <- readRDS("files/param_calib_multitarget//predictions/pred_231p_gpp_sampled_N3errors.rds")      # NOTE: for rendering remove prefix "vignettes/"
df_predict_231p_bigD13C_sampled <- readRDS("files/param_calib_multitarget//predictions/pred_231p_bigD13C_sampled_N3errors.rds")  # NOTE: for rendering remove prefix "vignettes/"
# df_predict_231p_bigD13C_sampled <- readr::read_rds(paste0("/storage/scratch/giub_geco/fbernhard/rsofun_doc_outputs/data/predictions/out_predict_",n_post,"_30000burnin__out_calib__scen222_DREAMzs-100000-0iter_8x3chains_on_CPU8x1_continued.rds_vj_sampled",n_err,".rds"))
# df_predict_231p_gpp_sampled     <- readr::read_rds(paste0("/storage/scratch/giub_geco/fbernhard/rsofun_doc_outputs/data/predictions/out_predict_",n_post,"_30000burnin__out_calib__scen222_DREAMzs-100000-0iter_8x3chains_on_CPU8x1_continued.rds_gpp_sampled",n_err,".rds"))

# add transparency to colors
t_col <- function(color, percent = 50, name = NULL) {
  #      color = color name
  #    percent = % transparency
  #       name = an optional name for the color
  rgb.val <- col2rgb(color) # Get RGB values for named color
  ## Make new color using input color as base and alpha set by transparency
  t.col <- rgb(
    rgb.val[1], rgb.val[2], rgb.val[3],
    max = 255,
    alpha = (100 - percent) * 255 / 100,
    names = name)
  invisible(t.col) ## Save the color
}

## copied from sofunCalVal package (and adapted to analyse_modobs3)
analyse_modobs3 <- function(df,
                            mod,
                            obs,
                            type = "points",
                            filnam = NA,
                            relative = FALSE,
                            lower_xlim = 0,
                            use_factor = NULL,
                            shortsubtitle = FALSE,
                            plot_subtitle = TRUE,
                            plot_linmod = TRUE,
                            pal = "batlowW",
                            ...) {
  require(ggplot2)
  require(dplyr)

  # if (identical(filnam, NA)) filnam <- "analyse_modobs.pdf"

  ## rename to 'mod' and 'obs' and remove rows with NA in mod or obs
  df <- df %>%
    as_tibble() %>%
    ungroup() %>%
    dplyr::select(dplyr::all_of(c(mod = mod, obs = obs))) %>%
    tidyr::drop_na(mod, obs)

  ## get linear regression (coefficients)
  linmod <- lm(obs ~ mod, data = df)

  ## construct metrics table using the 'yardstick' library
  df_metrics <- df %>%
    yardstick::metrics(obs, mod) %>%
    dplyr::bind_rows(tibble(.metric = "n", .estimator = "standard", .estimate = summarise(df, numb = n()) %>% unlist())) %>%
    dplyr::bind_rows(tibble(.metric = "slope", .estimator = "standard", .estimate = coef(linmod)[2])) %>%
    # dplyr::bind_rows( tibble( .metric = "nse",      .estimator = "standard", .estimate = hydroGOF::NSE( obs, mod, na.rm=TRUE ) ) ) %>%
    dplyr::bind_rows(tibble(.metric = "mean_obs", .estimator = "standard", .estimate = summarise(df, mean = mean(obs, na.rm = TRUE)) %>% unlist())) %>%
    dplyr::bind_rows(tibble(
      .metric = "prmse", .estimator = "standard",
      .estimate = dplyr::filter(., .metric == "rmse") %>% dplyr::select(.estimate) %>% unlist() /
        dplyr::filter(., .metric == "mean_obs") %>%
          dplyr::select(.estimate) %>%
          unlist()
    )) %>%
    dplyr::bind_rows(tibble(
      .metric = "pmae", .estimator = "standard",
      .estimate = dplyr::filter(., .metric == "mae") %>% dplyr::select(.estimate) %>% unlist() /
        dplyr::filter(., .metric == "mean_obs") %>%
          dplyr::select(.estimate) %>%
          unlist()
    )) %>%
    dplyr::bind_rows(tibble(.metric = "bias", .estimator = "standard", .estimate = dplyr::summarise(df, mean((mod - obs), na.rm = TRUE)) %>% unlist())) %>%
    dplyr::bind_rows(tibble(.metric = "pbias", .estimator = "standard", .estimate = dplyr::summarise(df, mean((mod - obs) / obs, na.rm = TRUE)) %>% unlist()))

  rsq_val <- df_metrics %>%
    dplyr::filter(.metric == "rsq") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()
  rmse_val <- df_metrics %>%
    dplyr::filter(.metric == "rmse") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()
  mae_val <- df_metrics %>%
    dplyr::filter(.metric == "mae") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()
  bias_val <- df_metrics %>%
    dplyr::filter(.metric == "bias") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()
  slope_val <- df_metrics %>%
    dplyr::filter(.metric == "slope") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()
  n_val <- df_metrics %>%
    dplyr::filter(.metric == "n") %>%
    dplyr::select(.estimate) %>%
    unlist() %>%
    unname()

  if (relative) {
    rmse_val <- rmse_val / mean(df$obs, na.rm = TRUE)
    bias_val <- bias_val / mean(df$obs, na.rm = TRUE)
  }

  rsq_lab <- format(rsq_val, digits = 2)
  rmse_lab <- format(rmse_val, digits = 3)
  mae_lab <- format(mae_val, digits = 3)
  bias_lab <- format(bias_val, digits = 3)
  slope_lab <- format(slope_val, digits = 3)
  n_lab <- format(n_val, digits = 3)

  results <- tibble(rsq = rsq_val, rmse = rmse_val, mae = mae_val, bias = bias_val, slope = slope_val, n = n_val)

  if (shortsubtitle) {
    subtitle <- bquote(italic(R)^2 == .(rsq_lab) ~ ~
      RMSE == .(rmse_lab))
  } else {
    subtitle <- bquote(italic(R)^2 == .(rsq_lab) ~ ~
      RMSE == .(rmse_lab) ~ ~
      bias == .(bias_lab) ~ ~
      slope == .(slope_lab) ~ ~
      italic(N) == .(n_lab))
  }

  if (type == "hex") {

    upper_xlim <- round(max(quantile(df$mod, 0.9999), quantile(df$obs, 0.9999)))
    stopifnot(is.null(lower_xlim) || (is.numeric(lower_xlim) && length(lower_xlim) == 1))
    lower_xlim <- ifelse(is.null(lower_xlim),
      round(min(quantile(df$mod, 0.0001), quantile(df$obs, 0.0001))),
      lower_xlim)

    ## ggplot hexbin
    gg <- df %>%
      ggplot2::ggplot(aes(x = mod, y = obs)) +
      geom_hex(bins = 50, show.legend = FALSE) +
      geom_abline(intercept = 0, slope = 1, linetype = "dotted") +
      # geom_hline(yintercept = 0, linetype = "dotted") +
      # geom_vline(xintercept = 0, linetype = "dotted") +
      coord_fixed() +
      xlim(lower_xlim, upper_xlim) +
      ylim(lower_xlim, upper_xlim) +
      theme_classic() +
      labs(x = mod, y = obs)

    if (pal == "batlowW") {
      gg <- gg + khroma::scale_fill_batlowW(trans = "log", reverse = TRUE)
    } else if (pal == "davos") {
      gg <- gg + khroma::scale_fill_davos(trans = "log", reverse = TRUE)
    }

    if (plot_subtitle) gg <- gg + labs(subtitle = subtitle)
    if (plot_linmod) gg <- gg + geom_smooth(method = "lm", color = "red", linewidth = 0.5, se = FALSE)

    if (!identical(filnam, NA)) {
      ggsave(filnam, width = 5, height = 5)
    }
  }

  return(list(df_metrics = df_metrics, gg = gg, linmod = linmod, results = results))
}

plot_predobs_gpp_timeseries3 <- function(ts_to_plot) {
  # separate obs
  df_tsplot_gpp_obs <- ts_to_plot |>
    select(sitename, target, date, obs, Scenario, dataset) |>
    distinct()

  # compute stats of sampled distributions before plotting them
  tibble_to_plot <- ts_to_plot |>
    filter(target == "gpp") |>
    # dtplyr::lazy_dt() |> # THIS WAS NEEDED FOR THE FULL DATA.FRAME
    group_by(Scenario, dataset, parameters, sitename, target, date, model_output_type, y_facet) |>
    summarise( # mod_no_err_p50 = quantile(mod_no_err, 0.5),
      modelled_p50 = quantile(modelled, 0.5),
      modelled_p95 = quantile(modelled, 0.95),
      modelled_p05 = quantile(modelled, 0.05)
    ) |>
    as_tibble() |>
    # and bind back obsevations
    left_join(
      df_tsplot_gpp_obs,
      by = join_by(sitename, target, date, Scenario, dataset))

  n_sites <- tibble_to_plot$sitename |> unique() |> length()

  pl_timeseries_gpp <- ggplot(
    data = tibble_to_plot,
    mapping = aes(x = date, y = modelled_p50)) +
    # Observations underneath (following Cameron 2022)
    geom_point(
      data = function(df) df |> select(-parameters, -model_output_type, -y_facet, -starts_with("modelled")) |> distinct(),
      # data = df_tsplot_gpp_obs, # variant 1, but does not allow %+%-replacement of underlying data
      mapping = aes(y = obs), color = "black", shape = 4, alpha = 0.5, size = 0.5) +
    # Structural uncertainty (including error model), a.k.a prediction band
    geom_ribbon(
      alpha = 0.5,
      data = function(df) df |> filter(model_output_type == "with struct. uncert."),
      mapping = aes(ymin = modelled_p05, ymax = modelled_p95, fill = "Post.+Error")) +
    # Parametric uncertainty (without error model, only parameter sampling), a.k.a confidence band
    geom_ribbon(
      alpha = 0.5,
      data = function(df) df |> filter(model_output_type == "rsofun"),
      mapping = aes(ymin = modelled_p05, ymax = modelled_p95, fill = "Posterior")) +
    geom_line(
      data = function(df) df |> filter(model_output_type == "rsofun"),
      mapping = aes(y = modelled_p50, color = "Posterior")) +
    # layout
    facet_wrap(~sitename, scales = "free_x") +
    theme_classic() + theme(legend.position = "bottom") +
    labs(
      x = NULL, #' Date',
      y = expression(paste("GPP (g C m"^-2, "s"^-1, ")"))
    ) +
    scale_fill_manual(NULL, aesthetics = c("colour", "fill"),
      breaks = c("Post.+Error",
        "Posterior"),
      values = c("Posterior"   = t_col("#29a274ff"),
        "Post.+Error" = t_col("#777055ff"))
    ) +
    scale_x_date(date_breaks = "12 months", date_labels = "%Y-%m")
}

# Time series of GPP (Figure C5)
dflong_gpp_train <- df_predict_231p_gpp_sampled |>
  mutate(Scenario = "231") |>
  # select only training sites
  filter(is_train0_test1 == 0) |>
  select(posterior_sample_id, error_sample_id, Scenario, is_train0_test1, is_MAP, sitename, target,
    obs,                 date, # these are target specific observation_metadata
    mod_no_err, mod_biasremoved_no_err, mod_biasremoved_with_err) |>
  # pivot the model_output_types to long
  tidyr::pivot_longer(c(mod_no_err, mod_biasremoved_no_err, mod_biasremoved_with_err),
    names_to = "model_output_type", values_to = "modelled") |>
  mutate(model_output_type = factor(
    model_output_type,
    levels = c("mod_no_err", "mod_biasremoved_no_err", "mod_biasremoved_with_err"),
    labels = c("rsofun",    "bias-corrected",        "with struct. uncert."))) |>
  # derive column `parameters` ("MAP" or "Posterior") from `is_MAP`
  mutate(is_MAP = factor(ifelse(is_MAP, "MAP", "Posterior"))) |>
  rename(parameters = is_MAP) |>
  # derive column `dataset` ("train" or "test") from column `is_train0_test1`
  mutate(is_train0_test1 = factor(ifelse(is_train0_test1 == 1, "test", "train"))) |>
  rename(dataset = is_train0_test1)

df_B3_timeseries <- dflong_gpp_train |>
  # remove the bias-corrected values for gpp since we did not fit a bias
  filter(!(model_output_type %in% c("bias-corrected"))) |>
  # select what to plot and how to name it
  mutate(y_facet = case_when(
    Scenario == "231" & model_output_type == "rsofun" &               parameters == "MAP"       ~ "MAP",
    Scenario == "231" & model_output_type == "rsofun" &               parameters == "Posterior" ~ "Posterior",
    Scenario == "231" & model_output_type == "with struct. uncert." & parameters == "Posterior" ~ "Post.+Error",
    # all else is not plotted
    TRUE ~ "remove") |> factor(levels = c("MAP", "Posterior", "Post.+Error"))) |>
  filter(y_facet != "remove")

# Plot raw predictions
pl_timeseries_gpp <- plot_predobs_gpp_timeseries3(df_B3_timeseries)
pl_timeseries_gpp

# fake output since prediction was run previous to vignette rendering
knitr::include_graphics("files/param_calib_multitarget/predictions/pl_timeseries_gpp.png", # NOTE: for rendering remove prefix "vignettes/"
  dpi = 100)

gpp_labs <- function(xNULL = FALSE) {
  list(labs(x = ifelse(xNULL, " ", "Predicted GPP (g C m-\u00B2 s-\u00B9)"),
    y = "Observed GPP (g C m-\u00B2 s-\u00B9)"))
}
bigD13C_labs <- function(xNULL = FALSE) {
  list(labs(x = ifelse(xNULL, " ", "Predicted Δ (\u2030)"),
    y = "Observed Δ (\u2030)"))
}

# compute scatters (and skills) for MAP of test set:
list_of_scenarios_to_loop_over <- list(
  s231 = list("bigD13C" = df_predict_231p_bigD13C_sampled,  "gpp" = df_predict_231p_gpp_sampled)
)

parameter_set <- "MAP"
list_of_scatters <- lapply(list_of_scenarios_to_loop_over, function(list_of_targets) {
  lapply(list_of_targets, function(df_target_prediction) {
    curr_target <- first(df_target_prediction$target)

    # (MAP, Posterior, Posterior+Error)
    # MAP:             is posterior_sample_id==0
    #                  e.g. filter(df_gpp, is_train0_test1 == 1,     is_MAP, error_sample_id==1)
    # Posterior:       just take one error sampling, but all (~25 posterior samples)
    #                  e.g. filter(df_gpp, is_train0_test1 == 1,             error_sample_id==1),
    # Posterior+Error: take all errors
    #                  e.g. filter(df_gpp, is_train0_test1 == 1,             error_sample_id>=1)
    df <- df_target_prediction %>%
      {
        if (parameter_set == "MAP") filter(., is_MAP) else .
      } %>% filter(error_sample_id == 1) # , is_train0_test1 == 1 # TODO: why was this removed?
    mod <- "mod_biasremoved_no_err"
    lower_xlim <- ifelse(curr_target == "gpp", list(0), list(NULL))[[1]]

    if (nrow(df) == 0) {
      return(list(gg = ggplot() + theme_void()))
    } else {
      return(analyse_modobs3(df, mod = mod, obs = "obs", type = "hex"))
    }
  })
})

### single comparison plot: ----
mark_as_target <- theme(panel.background = element_rect(fill = t_col("darkgreen", 80)))
pl_scatter_comparison <- cowplot::plot_grid(
  list_of_scatters$s231$bigD13C$gg + mark_as_target + bigD13C_labs(),
  list_of_scatters$s231$gpp$gg     + mark_as_target + gpp_labs()
)
pl_scatter_comparison

# fake output since prediction was run previous to vignette rendering
knitr::include_graphics("files/param_calib_multitarget/predictions/pl_scatter_comparison.png", # NOTE: for rendering remove prefix "vignettes/"
  dpi = 150)