Consolidate likelihoods v2

NAMESPACE

@@ -15,4 +15,5 @@ import(BayesianTools)
 import(GenSA)
 import(dplyr)
 importFrom(magrittr,"%>%")
+importFrom(stats,setNames)
 useDynLib(rsofun)

NEWS.md

@@ -1,6 +1,10 @@
 # rsofun (development version)
 
 * New `run_pmodel_onestep_f_bysite()` adds single-step leaf-level simulations of the P-model.
+* `calib_sofun()` now passes parameters as a _named_ vector to cost-functions for 
+  easier processing within cost-functions. Default cost-functions were updated,
+  but currently ignore the names. This is fully backward compatible, but allows
+  to use the names in user-created cost functions.
 
 # rsofun 5.1.0
 

R/calib_sofun.R

@@ -4,10 +4,13 @@
 #' calibration of SOFUN model parameters. 
 #' 
 #' @param drivers A data frame with driver data. See \code{\link{p_model_drivers}}
-#' for a description of the data structure.
+#' for a description of the data structure. Additional columns can optionally be 
+#' provided to \code{drivers} to control e.g. the processing within a 
+#' personalized cost function.
 #' @param obs A data frame containing observational data used for model
 #'  calibration. See \code{\link{p_model_validation}} for a description of the data
-#'  structure.
+#'  structure. Additional columns can optionally be provided to \code{obs} to 
+#'  control e.g. the processing within a personalized cost function.
 #' @param settings A list containing model calibration settings. 
 #'  See the 'P-model usage' vignette for more information and examples.
 #'  \describe{
@@ -31,16 +34,16 @@
 #'  }
 #' @param optim_out A logical indicating whether the function returns the raw
 #'  output of the optimization functions (defaults to TRUE).
-#' @param ... Optional arguments passed on to the cost function specified as
-#'  \code{settings$metric}. 
-#'  . 
+#' @param ... Optional arguments, simply passed on to the cost function. 
+#' 
 #' @return A named list containing the calibrated parameter vector `par` and
 #' the output object from the optimization `mod`. For more details on this
 #' output and how to evaluate it, see \link[BayesianTools:runMCMC]{runMCMC} (also
 #' \href{https://florianhartig.github.io/BayesianTools/articles/BayesianTools.html}{this post})
 #' and \link[GenSA]{GenSA}.
 #' @export
 #' @importFrom magrittr %>%
+#' @importFrom stats setNames
 #' @import GenSA BayesianTools
 #' 
 #' @examples
@@ -153,12 +156,16 @@ calib_sofun <- function(
     }
 
     # reformat parameters
-    pars <- as.data.frame(do.call("rbind", settings$par), row.names = FALSE)
+    pars <- as.data.frame(do.call("rbind", settings$par)) 
+         # NOTE: This keeps parameters as row names. 
+         # This does not change anything to previous behavior.
+         # But this agrees better with the example data 
+         # `dput(BayesianTools::VSEMgetDefaults())`
 
     priors  <- BayesianTools::createUniformPrior(
-      unlist(pars$lower),
-      unlist(pars$upper),
-      unlist(pars$init)
+      lower = unlist(pars$lower),
+      upper = unlist(pars$upper),
+      best  = unlist(pars$init)
     )
 
     # setup the bayes run, no message forwarding is provided
@@ -168,14 +175,18 @@ calib_sofun <- function(
         do.call(
           "cost",
           list(
-            par = random_par,
+            par = setNames(random_par, rownames(pars)),
+            # NOTE: if we could make use of setup$names from within the cost
+            # function then we wouldn't need this closure (using `pars`), but it
+            # appears that BayesianTools does not pass the names into the
+            # likelihood.
             obs = obs,
             drivers = drivers
           )
         )
       },
     prior = priors,
-    names = names(settings$par)
+    names = rownames(pars)
     )    
 
     # set bt control parameters

R/cost_likelihood_biomee.R

@@ -4,7 +4,7 @@
 #' computes the log-likelihood for the biomee model fitting several target 
 #' variables for a given set of parameters.
 #' 
-#' @param par A vector containing parameter values for \code{'phiRL',
+#' @param par A named vector containing parameter values for \code{'phiRL',
 #' 'LAI_light', 'tf_base', 'par_mort'} in that order, and for the error terms
 #' corresponding to the target variables, e.g. \code{'err_GPP'} if GPP is a target. 
 #' Make sure that
@@ -34,8 +34,11 @@
 #' # BiomeE model parameter values
 #' # and the example data
 #' cost_likelihood_biomee(
-#'  par = c(3.5, 3.5, 1, 1,    # model params
-#'          0.5),              # err_GPP
+#'  par = c(phiRL = 3.5, 
+#'          LAI_light = 3.5, 
+#'          tf_base = 1, 
+#'          par_mort = 1,    # model params
+#'          err_GPP = 0.5),  # err_GPP
 #'  obs = biomee_validation,
 #'  drivers = biomee_gs_leuning_drivers,
 #'  targets = c("GPP")
@@ -52,6 +55,7 @@ cost_likelihood_biomee <- function(
   # predefine variables for CRAN check compliance
   GPP <- LAI <- Density12 <- plantC <- NULL
 
+  par <- unname(par) # reproduces previous behavior, when par was unnamed
   # Add changed model parameters to drivers, overwriting where necessary.
   drivers$params_species[[1]]$phiRL[]  <- par[1]
   drivers$params_species[[1]]$LAI_light[]  <- par[2]

R/cost_likelihood_pmodel.R

@@ -6,7 +6,7 @@
 #' observed values and with standard deviation given as an input parameter 
 #' (calibratable).
 #' 
-#' @param par A vector of values for the parameters to be calibrated, including
+#' @param par A named vector of values for the parameters to be calibrated, including
 #' a subset of model parameters (described in \code{\link{runread_pmodel_f}}),
 #' in order, and error terms
 #' for each target variable (for example \code{'gpp_err'}), in the same order as
@@ -56,8 +56,10 @@
 #' # temperature dependence of kphio 
 #' # and example data
 #' cost_likelihood_pmodel(
-#'  par = c(0.05, -0.01, 1,     # model parameters
-#'          2),                # err_gpp
+#' par = c(kphio       = 0.05, 
+#'         kphio_par_a = -0.01, 
+#'         kphio_par_b = 1,     # model parameters
+#'         err_gpp     = 2),    # err_gpp
 #'  obs = p_model_validation,
 #'  drivers = p_model_drivers,
 #'  targets = c('gpp'),
@@ -82,6 +84,7 @@ cost_likelihood_pmodel <- function(
 ){
   # predefine variables for CRAN check compliance
   sitename <- data <- gpp_mod <- NULL
+  par <- unname(par) # reproduces previous behavior, when par was unnamed
 
   ## check input parameters
   if( (length(par) + length(par_fixed)) != (9 + length(targets)) ){
@@ -95,7 +98,7 @@ cost_likelihood_pmodel <- function(
                          'beta_unitcostratio', 'rd_to_vcmax', 
                          'tau_acclim', 'kc_jmax')
 
-  if(!is.null(par_fixed)){
+  if(!is.null(par_fixed) && length(par)>0){
     params_modl <- list()
     # complete with calibrated values
     i <- 1 # start counter
@@ -107,9 +110,12 @@ cost_likelihood_pmodel <- function(
         params_modl[[par_name]] <- par_fixed[[par_name]]  # use fixed par value
       }
     }
+  }else if(length(par)==0){                # no parameters calibrated
+    params_modl <- as.list(par_fixed[calib_param_names])
+    par <- par_fixed[grepl("err_",names(par_fixed))]
   }else{
     params_modl <- as.list(par[1:9])       # all parameters calibrated
-    names(params_modl) <- calib_param_names
+    names(params_modl) <- calib_param_names# TODO: problematic, since it assumes they are in the right order
   }
 
   ## run the model

analysis/01-sensitivity-analysis.R

@@ -13,7 +13,7 @@ set.seed(432)
 
 # Define log-likelihood function
 ll_pmodel <- function(
-    par_v                 # a vector of all calibratable parameters including errors
+    par_v                 # a named vector of all calibratable parameters including errors
 ){
   rsofun::cost_likelihood_pmodel(        # likelihood cost function from package
     par_v,