geco-bern · stineb · May 23, 2024 · Jun 19, 2024 · Aug 8, 2024 · Aug 8, 2024
diff --git a/R/cost_likelihood_biomee.R b/R/cost_likelihood_biomee.R
@@ -16,6 +16,9 @@
 #' @param targets A character vector indicating the target variables for which the
 #' optimization will be done. This should be a subset of \code{c("GPP", "LAI",
 #' "Density", "Biomass")}.
+#' @param verbose A logical specifying whether to print parameter and log-likelihood 
+#' values during each iteration of the calibration as message. Defaults to 
+#' \code{FALSE}
 #' 
 #' @return The log-likelihood of the simulated 
 #' targets by the biomee model versus the observed targets. 
@@ -46,7 +49,8 @@ cost_likelihood_biomee <- function(
   par,
   obs,
   drivers,
-  targets
+  targets,
+  verbose = FALSE
 ){
 
   # predefine variables for CRAN check compliance
@@ -108,6 +112,18 @@ cost_likelihood_biomee <- function(
   if(is.nan(ll) || is.na(ll) | ll == 0){
     ll <- -Inf
   }
+
+  if (verbose){
+    message(
+      paste("Par = ", paste(drivers$params_species, collapse = ", "))
+    )
+    message(
+      paste("Log-likelihood = ", ll)
+    )
+    message(
+      paste("--------------------------")
+    )
+  }
 
   return(ll)
 }
diff --git a/R/cost_likelihood_pmodel.R b/R/cost_likelihood_pmodel.R
@@ -28,6 +28,9 @@
 #' \code{FALSE}.
 #' @param ncores An integer specifying the number of cores used for parallel
 #' computing. Defaults to 2.
+#' @param verbose A logical specifying whether to print parameter and log-likelihood 
+#' values during each iteration of the calibration as message. Defaults to 
+#' \code{FALSE}
 #' 
 #' @return The log-likelihood of the observed target values, assuming that they
 #' are independent, normally distributed and centered on the predictions
@@ -78,7 +81,8 @@ cost_likelihood_pmodel <- function(
     targets,
     par_fixed = NULL,   # non-calibrated model parameters
     parallel = FALSE,
-    ncores = 2
+    ncores = 2,
+    verbose = FALSE
 ){
   # predefine variables for CRAN check compliance
   sitename <- data <- gpp_mod <- NULL
@@ -215,6 +219,18 @@ cost_likelihood_pmodel <- function(
 
   # trap boundary conditions
   if(is.nan(ll) | is.na(ll) | ll == 0){ll <- -Inf}
+
+  if (verbose){
+    message(
+      paste("Par = ", paste(params_modl, collapse = ", "))
+    )
+    message(
+      paste("Log-likelihood = ", ll)
+    )
+    message(
+      paste("--------------------------")
+    )
+  }
 
   return(ll)
 }

diff --git a/R/cost_rmse_biomee.R b/R/cost_rmse_biomee.R
@@ -12,6 +12,9 @@
 #' @param obs A nested data frame of observations, following the structure of \code{biomee_validation},
 #' for example.
 #' @param drivers A nested data frame of driver data, for example \code{biomee_gs_leuning_drivers}.
+#' @param verbose A logical specifying whether to print parameter and cost (RMSE) 
+#' values during each iteration of the calibration as message. Defaults to 
+#' \code{FALSE}
 #' 
 #' @return The root mean squared error (RMSE) between the observed and simulated
 #' values of \code{'GPP','LAI','Density'} and \code{'Biomass'} (all variables
@@ -38,7 +41,8 @@
 cost_rmse_biomee <- function(
     par,
     obs,
-    drivers
+    drivers,
+    verbose = FALSE
 ){
 
   # predefine variables for CRAN check compliance
@@ -85,5 +89,17 @@ cost_rmse_biomee <- function(
   ## Calculate cost (RMSE) across the N targets
   cost <- mean(dff$error_rel^2, na.rm = TRUE)
 
+  if (verbose){
+    message(
+      paste("Par = ", paste(drivers$params_species, collapse = ", "))
+    )
+    message(
+      paste("RMSE = ", cost)
+    )
+    message(
+      paste("--------------------------")
+    )
+  }
+
   return(cost)
 }
diff --git a/R/cost_rmse_pmodel.R b/R/cost_rmse_pmodel.R
@@ -27,6 +27,9 @@
 #' \code{FALSE}.
 #' @param ncores An integer specifying the number of cores used for parallel
 #' computing. Defaults to 2.
+#' @param verbose A logical specifying whether to print parameter and cost (RMSE) 
+#' values during each iteration of the calibration as message. Defaults to 
+#' \code{FALSE}
 #' 
 #' @return The root mean squared error (RMSE) between observed values and P-model
 #' predictions. The RMSE is computed for each target separately and then aggregated
@@ -75,7 +78,8 @@ cost_rmse_pmodel <- function(
     target_weights = NULL, # if using several targets, how are the individual 
                            # RMSE weighted? named vector
     parallel = FALSE,
-    ncores = 2
+    ncores = 2,
+    verbose = FALSE
 ){
 
   # predefine variables for CRAN check compliance
@@ -204,6 +208,18 @@ cost_rmse_pmodel <- function(
   }else{
     cost <- mean(rmse, na.rm = TRUE)
   }
+
+  if (verbose){
+    message(
+      paste("Par = ", paste(params_modl, collapse = ", "))
+    )
+    message(
+      paste("RMSE = ", cost)
+    )
+    message(
+      paste("--------------------------")
+    )
+  }
 
   return(cost)
 }
diff --git a/R/run_pmodel_f_bysite.R b/R/run_pmodel_f_bysite.R
@@ -221,6 +221,9 @@ run_pmodel_f_bysite <- function(
       "temp",
       "rain",
       "vpd",
+      "ppfd",
+      "netrad",
+      "fsun",
       "snow",
       "co2",
       "fapar",
@@ -282,7 +285,7 @@ run_pmodel_f_bysite <- function(
     }
 
     # Check model parameters
-    if( sum( names(params_modl) %in% c('kphio', 'kphio_par_a', 'kphio_par_b',
+    if ( sum( names(params_modl) %in% c('kphio', 'kphio_par_a', 'kphio_par_b',
                                               'soilm_thetastar', 'soilm_betao',
                                               'beta_unitcostratio', 'rd_to_vcmax', 
                                               'tau_acclim', 'kc_jmax')
@@ -302,13 +305,13 @@ run_pmodel_f_bysite <- function(
     in_netrad <- FALSE  # net radiation is currently ignored as a model forcing, but is internally simulated by SPLASH.
 
     # Check if fsun is available
-    if(! (in_ppfd & in_netrad)){
+    if (! (in_ppfd & in_netrad)){
       # fsun must be available when one of ppfd or netrad is missing
-      if(any(is.na(forcing$fsun))) continue <- FALSE
+      if (any(is.na(forcing$fsun))) continue <- FALSE
     }
   }
 
-  if(continue){
+  if (continue){
 
     # convert to matrix
     forcing <- as.matrix(forcing)

diff --git a/analysis/demo_bug.R b/analysis/demo_bug.R
@@ -3,7 +3,8 @@ library(dplyr)
 library(ggplot2)
 library(purrr)
 
-## BiomeE (original with gs-leuning) -----------
+## BiomeE -----------
+# original with gs-leuning photosynthesis
 nruns <- 3
 test_biomee_gs_leuning <- function(){
   # run the model
@@ -26,7 +27,13 @@ vec_test <- purrr::map_lgl(
 # did any simulation fail?
 if (any(!vec_test)){
   stop(
-    "At least one BiomeE-gs-leuning simulation failed."
+    paste(
+      "At least one BiomeE-gs-leuning simulation failed. Successful run? ",
+      paste(
+        vec_test,
+        collapse = ", "
+      )
+    )
   )
 } else {
   message(
@@ -57,7 +64,13 @@ vec_test <- purrr::map_lgl(
 # did any simulation fail?
 if (any(!vec_test)){
   stop(
-    "At least one BiomeEP simulation failed."
+    paste(
+      "At least one BiomeEP simulation failed. Successful run? ",
+      paste(
+        vec_test,
+        collapse = ", "
+      )
+    )
   )
 } else {
   message(

diff --git a/man/cost_likelihood_biomee.Rd b/man/cost_likelihood_biomee.Rd
diff --git a/man/cost_likelihood_pmodel.Rd b/man/cost_likelihood_pmodel.Rd
diff --git a/man/cost_rmse_biomee.Rd b/man/cost_rmse_biomee.Rd
diff --git a/man/cost_rmse_pmodel.Rd b/man/cost_rmse_pmodel.Rd
diff --git a/src/Makevars.in b/src/Makevars.in
@@ -8,6 +8,8 @@
 # -fsanitize=address: Using the Address Sanitizer https://rstudio.github.io/r-manuals/r-exts/Debugging.html#using-the-address-sanitizer
 # PKG_FFLAGS = -frecursive -fbounds-check -Wall -Wextra -Wpedantic -g -O0 -fsanitize=address
 
+# PKG_FFLAGS = -g -fsanitize=address
+
 # C objects
 C_OBJS = wrappersc.o
 

diff --git a/src/vegetation_biomee.mod.f90 b/src/vegetation_biomee.mod.f90
@@ -60,7 +60,7 @@ subroutine vegn_CNW_budget( vegn, forcing, init )
     do i = 1, vegn%n_cohorts
 
       cc => vegn%cohorts(i)
-      associate ( sp => spdata(cc%species) )
+      ! associate ( sp => spdata(cc%species) )
 
       ! increment the cohort age
       cc%age = cc%age + myinterface%dt_fast_yr
@@ -76,7 +76,7 @@ subroutine vegn_CNW_budget( vegn, forcing, init )
       cc%plabl%c%c12 = cc%plabl%c%c12 + cc%npp
       cc%plabl%n%n14 = cc%plabl%n%n14 + cc%fixedN
 
-      end associate
+      ! end associate
     enddo ! all cohorts
 
     ! update soil carbon

diff --git a/vignettes/new_cost_function.Rmd b/vignettes/new_cost_function.Rmd
@@ -61,7 +61,8 @@ pars_calib_rmse <- calib_sofun(
     tau_acclim         = 30.0,
     kc_jmax            = 0.41
   ),
-  targets = "gpp"           # define target variable GPP
+  targets = "gpp",           # define target variable GPP
+  verbose = FALSE            # set to TRUE to print parameter and cost values during each iteration of the calibration
 )
 ```
 

diff --git a/vignettes/pmodel_use.Rmd b/vignettes/pmodel_use.Rmd
@@ -211,3 +211,4 @@ ggplot(data = df_gpp_plot) +
 ```
 
 For details on the optimization settings we refer to the manuals of [GenSA](https://cran.r-project.org/package=GenSA) and [BayesianTools](https://github.com/florianhartig/BayesianTools).
+
Original file line number	Diff line number	Diff line change
Expand Up		@@ -211,3 +211,4 @@ ggplot(data = df_gpp_plot) +
		```

		For details on the optimization settings we refer to the manuals of [GenSA](https://cran.r-project.org/package=GenSA) and [BayesianTools](https://github.com/florianhartig/BayesianTools).