## ----setup, include = FALSE--------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(rsofun) library(dplyr) library(ggplot2) ## ----run GenSA calibration, eval = FALSE-------------------------------------- # # Define calibration settings and parameter ranges from previous work # settings_rmse <- list( # method = 'GenSA', # minimizes the RMSE # metric = cost_rmse_pmodel, # our cost function returning the RMSE # control = list( # control parameters for optimizer GenSA # maxit = 100), # par = list( # bounds for the parameter space # kphio = list(lower=0.02, upper=0.2, init=0.05) # ) # ) # # # Calibrate the model and optimize the free parameters using # # demo datasets # pars_calib_rmse <- calib_sofun( # # calib_sofun arguments: # drivers = p_model_drivers, # obs = p_model_validation, # settings = settings_rmse, # # extra arguments passed to the cost function: # par_fixed = list( # fix all other parameters # kphio_par_a = 0.0, # set to zero to disable temperature-dependence # # of kphio, setup ORG # kphio_par_b = 1.0, # soilm_thetastar = 0.6 * 240, # to recover paper setup with soil moisture stress # soilm_betao = 0.0, # beta_unitcostratio = 146.0, # rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous # tau_acclim = 30.0, # kc_jmax = 0.41 # ), # targets = "gpp" # define target variable GPP # ) ## ----include = FALSE---------------------------------------------------------- # fake output since calibration isn't run pars_calib_rmse <- list( par = c(kphio = 0.03580938038563619), mod = list( value = 1.20014048299052467, par = c(kphio = 0.03580938038563619), trace.mat = matrix( c(1, 5230, 4.241255082597149, 1.20014048299052467), nrow = 1L, ncol = 4L, dimnames = list(NULL, c("nb.steps", "temperature", "function.value", "current.minimum")) ), counts = 286L ) ) ## ----------------------------------------------------------------------------- pars_calib_rmse ## ----run Bayesian calibration, eval = FALSE----------------------------------- # # Define calibration settings # settings_likelihood <- list( # method = 'BayesianTools', # metric = cost_likelihood_pmodel, # our cost function # control = list( # optimization control settings for # sampler = 'DEzs', # BayesianTools::runMCMC # settings = list( # burnin = 1500, # iterations = 3000 # )), # par = list( # kphio = list(lower = 0, upper = 0.2, init = 0.05), # kphio_par_a = list(lower = -0.5, upper = 0.5, init = -0.1), # kphio_par_b = list(lower = 10, upper = 40, init =25), # err_gpp = list(lower = 0.1, upper = 4, init = 0.8) # ) # ) # # # Calibrate the model and optimize the free parameters using # # demo datasets # pars_calib_likelihood <- calib_sofun( # # calib_sofun arguments: # drivers = p_model_drivers, # obs = p_model_validation, # settings = settings_likelihood, # # extra arguments passed ot the cost function: # par_fixed = list( # fix all other parameters # soilm_thetastar = 0.6 * 240, # to recover paper setup with soil moisture stress # soilm_betao = 0.0, # beta_unitcostratio = 146.0, # rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous # tau_acclim = 30.0, # kc_jmax = 0.41 # ), # targets = "gpp" # ) ## ----eval = FALSE, include = FALSE-------------------------------------------- # # store output since calibration isn't run # saveRDS(pars_calib_likelihood, "files/new_cost_function.Rmd__pars_calib_likelihood.RDS", compress = "xz") ## ----include = FALSE---------------------------------------------------------- # fake output since calibration isn't run pars_calib_likelihood <- readRDS("files/new_cost_function.Rmd__pars_calib_likelihood.RDS") ## ----------------------------------------------------------------------------- pars_calib_likelihood ## ----run concatenated calibration, eval = FALSE------------------------------- # # Define calibration settings for two targets # settings_joint_likelihood <- list( # method = "BayesianTools", # metric = cost_likelihood_pmodel, # control = list( # sampler = "DEzs", # settings = list( # burnin = 1500, # kept artificially low # iterations = 3000 # )), # par = list(kc_jmax = list(lower = 0.2, upper = 0.8, init = 0.41), # uniform priors # err_gpp = list(lower = 0.001, upper = 4, init = 1), # err_vcmax25 = list(lower = 0.000001, upper = 0.0001, init = 0.00001)) # ) # # # Run the calibration on the concatenated data # par_calib_join <- calib_sofun( # drivers = rbind(p_model_drivers, # p_model_drivers_vcmax25), # obs = rbind(p_model_validation, # p_model_validation_vcmax25), # # settings = settings_joint_likelihood, # # arguments for the cost function # par_fixed = list( # fix parameter value from previous calibration # kphio = 0.041, # kphio_par_a = 0.0, # kphio_par_b = 16, # soilm_thetastar = 0.6 * 240, # to recover paper setup with soil moisture stress # soilm_betao = 0.0, # beta_unitcostratio = 146.0, # rd_to_vcmax = 0.014, # value from Atkin et al. 2015 for C3 herbaceous # tau_acclim = 30.0 # ), # targets = c('gpp', 'vcmax25') # ) ## ----eval = FALSE, include = FALSE-------------------------------------------- # # store output since calibration isn't run # saveRDS(par_calib_join, "files/new_cost_function.Rmd__par_calib_join.RDS", compress = "xz") ## ----include = FALSE---------------------------------------------------------- # fake output since calibration isn't run par_calib_join <- readRDS("files/new_cost_function.Rmd__par_calib_join.RDS") ## ----------------------------------------------------------------------------- par_calib_join ## ----eval = FALSE------------------------------------------------------------- # # Define the custom cost function to be used # cost_mae <- function(par, obs, drivers, my_own_message){ # # Your code # browser() # can facilitate the development, remove afterwards # } # # # Define calibration settings and parameter ranges # settings_mae <- list( # method = 'GenSA', # metric = cost_mae, # directly uses the custom cost function # control = list( # maxit = 100 # ), # par = list( # soilm_thetastar = list(lower=0.0, upper=3000, init=0.6*240), # soilm_betao = list(lower=0, upper=1, init=0.2) # ) # ) # # # Calibrate the model and optimize the free parameters # pars_calib_mae <- calib_sofun( # drivers = p_model_drivers, # obs = p_model_validation, # settings = settings_mae, # # optional arguments if needed in the cost function # my_own_message = "Hi from inside the cost_mae function." # ) ## ----eval = FALSE------------------------------------------------------------- # cost_mae <- function(par, obs, drivers, my_own_message){ # # # Set values for the list of calibrated and non-calibrated model parameters # params_modl <- list( # kphio = 0.09423773, # kphio_par_a = 0.0, # kphio_par_b = 25, # soilm_thetastar = par[["soilm_thetastar"]], # soilm_betao = par[["soilm_betao"]], # beta_unitcostratio = 146.0, # rd_to_vcmax = 0.014, # tau_acclim = 30.0, # kc_jmax = 0.41 # ) # # # Run the model # df <- runread_pmodel_f( # drivers, # par = params_modl, # makecheck = TRUE, # parallel = FALSE # ) # # # Your code to compute the cost # print(my_own_message) # useless, but showcases how to use additional arguments # browser() # can facilitate the development, remove afterwards # } ## ----define custom cost function, eval = TRUE--------------------------------- cost_mae <- function(par, obs, drivers){ # Set values for the list of calibrated and non-calibrated model parameters params_modl <- list( kphio = 0.09423773, kphio_par_a = 0.0, kphio_par_b = 25, soilm_thetastar = par[["soilm_thetastar"]], soilm_betao = par[["soilm_betao"]], beta_unitcostratio = 146.0, rd_to_vcmax = 0.014, tau_acclim = 30.0, kc_jmax = 0.41 ) # Set values for the list of calibrated and non-calibrated model parameters # Run the model df <- runread_pmodel_f( drivers = drivers, par = params_modl, makecheck = FALSE, parallel = FALSE ) # Clean model output to compute cost df <- df %>% dplyr::select(sitename, data) %>% tidyr::unnest(data) # Clean validation data to compute cost obs <- obs %>% dplyr::select(sitename, data) %>% tidyr::unnest(data) %>% dplyr::rename('gpp_obs' = 'gpp') # rename for later # Left join model output with observations by site and date df <- dplyr::left_join(df, obs, by = c('sitename', 'date')) # Compute mean absolute error cost <- mean(abs(df$gpp - df$gpp_obs), na.rm = TRUE) # Return the computed cost return(cost) # browser() # can facilitate the development, remove afterwards } ## ----run custom calibration, eval=FALSE--------------------------------------- # # Define calibration settings and parameter ranges # settings_mae <- list( # method = 'GenSA', # metric = cost_mae, # our cost function # control = list( # maxit = 100), # par = list( # soilm_thetastar = list(lower=0.0, upper=3000, init=0.6*240), # soilm_betao = list(lower=0, upper=1, init=0.2) # ) # ) # # # Calibrate the model and optimize the free parameters # pars_calib_mae <- calib_sofun( # drivers = p_model_drivers, # obs = p_model_validation, # settings = settings_mae # ) ## ----eval = FALSE, include = FALSE-------------------------------------------- # # store output since calibration isn't run # saveRDS(pars_calib_mae, "files/new_cost_function.Rmd__pars_calib_mae.RDS", compress = "xz") ## ----include = FALSE---------------------------------------------------------- # fake output since calibration isn't run pars_calib_mae <- readRDS("files/new_cost_function.Rmd__pars_calib_mae.RDS") ## ----------------------------------------------------------------------------- pars_calib_mae