Add vignette

.Rbuildignore

@@ -4,6 +4,9 @@
 # reason, there should not be any blank comment lines, because they would
 # exclude any files with `#` in the filename, possibly causing confusing issues.
 
+^Rhistory$
+# History files
+
 ^.*\.Rproj$
 # Designates the directory as an RStudio Project.
 

DESCRIPTION

@@ -8,14 +8,19 @@ Authors@R: c(
            comment = c(ORCID = "0000-0002-4912-9783")),
     person("BioCroField authors", role = "cph")
     )
-License: MIT + file LICENSE
-Encoding: UTF-8
-LazyData: true
 Depends:
     R (>= 3.6.0)
 Imports:
     BioCro (>= 3.2.0)
 Suggests:
+    dfoptim,
+    lattice,
+    knitr,
+    rmarkdown,
     testthat (>= 3.0.0)
+VignetteBuilder: knitr
+License: MIT + file LICENSE
+Encoding: UTF-8
+LazyData: true
 URL: https://github.com/BioCro/BioCroField, https://biocro.github.io/BioCroValidation/
 Config/testthat/edition: 3

NAMESPACE

@@ -1,2 +1,4 @@
+export(bounds_table)
 export(objective_function)
+export(update_model)
 export(write_model)

NEWS.md

@@ -37,7 +37,10 @@ be directly added to this file to describe the related changes.
 
 - Added 2002 and 2005 SoyFACE biomass data.
 
-- Added the `objective_function` function.
+- Added several new functions: `objective_function`, `update_model`, and
+  `bounds_table`
+
+- Added a vignette illustrating how to perform a model parameterization.
 
 # Changes in BioCroValidation Version 0.1.0
 

R/bounds_table.R

@@ -0,0 +1,154 @@
+# A helper function for checking the bounds list for mistakes; if an issue is
+# found, this function will throw an error; otherwise it will be silent with no
+# return value.
+check_bounds_list <- function(bounds_list, independent_args) {
+    # Must be a list of named elements
+    if (!is.list(bounds_list) | is.null(names(bounds_list))) {
+        stop('`bounds_list` must be a list of named elements')
+    }
+
+    # Must contain all elements in independent_args
+    missing_element <- sapply(names(independent_args), function(x) {
+        !x %in% names(bounds_list)
+    })
+
+    if (any(missing_element)) {
+        msg <- paste0(
+            'The following elements were included in ',
+            '`independent_args` but not `bounds_list`: ',
+            paste(names(independent_args)[missing_element], collapse = ', ')
+        )
+        stop(msg)
+    }
+
+    # Each element must have length 2
+    length_two <- sapply(bounds_list, function(x) {
+        xlen <- length(x)
+
+        if (is.finite(xlen)) {
+            length(x) == 2
+        } else {
+            FALSE
+        }
+    })
+
+    if (any(!length_two)) {
+        msg <- paste0(
+            'The following elements of `bounds_list` do not have a length of 2: ',
+            paste(names(bounds_list)[!length_two], collapse = ', ')
+        )
+        stop(msg)
+    }
+
+    # Each element must be numeric
+    not_numeric <- sapply(bounds_list, function(x) {!is.numeric(x)})
+
+    if (any(not_numeric)) {
+        msg <- paste0(
+            'The following elements of `bounds_list` are not numeric: ',
+            paste(names(bounds_list)[not_numeric], collapse = ', ')
+        )
+        stop(msg)
+    }
+
+    return(invisible(NULL))
+}
+
+# A helper function for checking the initial guess for mistakes; if an issue is
+# found, this function will throw an error or a warning; otherwise it will be
+# silent with no return value.
+check_initial_ind_arg_values <- function(
+    independent_args,
+    lbounds,
+    ubounds,
+    initial_ind_arg_values
+)
+{
+    # Check the length
+    if (length(initial_ind_arg_values) != length(independent_args)) {
+        stop('`initial_ind_arg_values` must have the same length as `independent_args`')
+    }
+
+    # Check to make sure the initial values are not outside the bounds
+    outside_bounds <- sapply(seq_along(initial_ind_arg_values), function(i) {
+        initial_ind_arg_values[i] < lbounds[i] | initial_ind_arg_values[i] > ubounds[i]
+    })
+
+    if (any(outside_bounds)) {
+        msg <- paste0(
+            'The initial values for the following arguments lie outside the bounds: ',
+            paste(names(independent_args)[outside_bounds], collapse = ', ')
+        )
+        stop(msg)
+    }
+
+    # Check to see if any initial values are on the bounds
+    eps <- sqrt(.Machine$double.eps)
+
+    on_bounds <- sapply(seq_along(initial_ind_arg_values), function(i) {
+        abs(initial_ind_arg_values[i] - lbounds[i]) < eps |
+            abs(initial_ind_arg_values[i] - ubounds[i]) < eps
+    })
+
+    if (any(on_bounds)) {
+        msg <- paste0(
+            'The initial values for the following arguments lie on the ',
+            'bounds, which can be problematic for some optimizers: ',
+            paste(names(independent_args)[on_bounds], collapse = ', ')
+        )
+        warning(msg)
+    }
+
+    return(invisible(NULL))
+}
+
+bounds_table <- function(
+    independent_args,
+    bounds_list,
+    initial_ind_arg_values = NULL
+)
+{
+    # Check the bounds_list
+    check_bounds_list(bounds_list, independent_args)
+
+    # Get an ordering for the elements of `bounds_list` so they match the order
+    # of elements in `independent_args`; note that this will also exclude any
+    # elements of `bounds_list` that are not included in `independent_args`.
+    ordering <- match(
+        names(independent_args),
+        names(bounds_list)
+    )
+
+    bounds_list <- bounds_list[ordering]
+
+    # Get the lower and upper bounds
+    lbounds <- sapply(bounds_list, min)
+    ubounds <- sapply(bounds_list, max)
+
+    # Form the bounds table
+    bounds_table <- data.frame(
+        arg_name = names(independent_args),
+        lower = lbounds,
+        upper = ubounds,
+        stringsAsFactors = FALSE
+    )
+
+    # Include initial values in the table if they were provided
+    if (!is.null(initial_ind_arg_values)) {
+        # Check the values
+        check_initial_ind_arg_values(
+            independent_args,
+            lbounds,
+            ubounds,
+            initial_ind_arg_values
+        )
+
+        # Include the values
+        bounds_table$initial_value <- initial_ind_arg_values
+    }
+
+    # Remove row names and return the table
+    rownames(bounds_table) <- NULL
+
+    bounds_table
+}

R/objective_function.R

@@ -23,7 +23,8 @@ objective_function <- function(
     check_args_to_vary(
         independent_args,
         dependent_arg_function,
-        data_driver_pairs
+        data_driver_pairs,
+        verbose_startup
     )
 
     # Get the model runners

R/objective_function_helpers.R

@@ -62,6 +62,13 @@ get_model_runner <- function(
                     c(x, as.numeric(dependent_arg_function(x_for_dep_arg_func)))
                 }
 
+                if (any(!is.finite(x_for_partial))) {
+                    stop(
+                        'At least one independent or dependent argument ',
+                        'value is not finite'
+                    )
+                }
+
                 initial_res <- partial_func(x_for_partial)
 
                 if (is.null(post_process_function)) {
@@ -272,13 +279,17 @@ one_error <- function(
     normalization
 )
 {
-    qw <- if (predicted < observed) {
-        quantity_weight[1] # Underprediction
+    if (!is.finite(predicted)) {
+        NA
     } else {
-        quantity_weight[2] # Overprediction
-    }
+        qw <- if (predicted < observed) {
+            quantity_weight[1] # Underprediction
+        } else {
+            quantity_weight[2] # Overprediction
+        }
 
-    (observed - predicted)^2 * qw * ddp_weight * var_weight / normalization
+        (observed - predicted)^2 * qw * ddp_weight * var_weight / normalization
+    }
 }
 
 # Helping function for returning a failure value
@@ -348,8 +359,14 @@ error_from_res <- function(
 
     # Return the sum of the penalty and error terms, or the individual errors
     if (return_terms) {
+        error_terms_by_quantity <- as.list(tapply(
+            errors,
+            long_form_data_table[['quantity_name']],
+            sum
+        ))
+
         list(
-            least_squares_term = error_sum,
+            least_squares_terms = error_terms_by_quantity,
             extra_penalty = penalty
         )
     } else {
@@ -386,7 +403,18 @@ get_obj_fun <- function(
     regularization_method
 )
 {
-    function(x, lambda = 0, return_terms = FALSE) {
+    function(x, lambda = 0, return_terms = FALSE, debug_mode = FALSE) {
+        if (debug_mode) {
+            msg <- paste0(
+                '\nTime: ',
+                Sys.time(),
+                '    Independent argument values: ',
+                paste(x, collapse = ', '),
+                '\n'
+            )
+            cat(msg)
+        }
+
         errors <- lapply(seq_along(model_runners), function(i) {
             runner <- model_runners[[i]]
             res    <- runner(x)
@@ -405,15 +433,36 @@ get_obj_fun <- function(
         reg_penalty <- regularization_penalty(x, regularization_method, lambda)
 
         if (return_terms) {
-            list(
+            error_metric_terms <- list(
                 terms_from_data_driver_pairs = stats::setNames(
                     errors,
                     names(model_runners)
                 ),
                 regularization_penalty = reg_penalty
             )
+
+            if (debug_mode) {
+                cat(paste0('Time: ', Sys.time()), '    Error metric terms: ')
+                utils::str(error_metric_terms)
+                cat('\n')
+            }
+
+            error_metric_terms
         } else {
-            sum(as.numeric(errors)) + reg_penalty
+            error_metric <- sum(as.numeric(errors)) + reg_penalty
+
+            if (debug_mode) {
+                msg <- paste0(
+                    'Time: ',
+                    Sys.time(),
+                    '    Error metric: ',
+                    error_metric,
+                    '\n'
+                )
+                cat(msg)
+            }
+
+            error_metric
         }
     }
 }

R/objective_function_input_checks.R

@@ -148,7 +148,8 @@ check_data_driver_pairs <- function(base_model_definition, data_driver_pairs) {
 check_args_to_vary <- function(
     independent_args,
     dependent_arg_function,
-    data_driver_pairs
+    data_driver_pairs,
+    verbose
 )
 {
     # Make sure the independent arguments have names
@@ -168,6 +169,15 @@ check_args_to_vary <- function(
         }
     }
 
+    # Print argument names, if necessary
+    if (verbose) {
+        cat('\nThe independent arguments and their initial values:\n\n')
+        utils::str(independent_args)
+
+        cat('\nThe dependent arguments and their initial values:\n\n')
+        utils::str(dependent_arg_function(independent_args))
+    }
+
     # Make sure no drivers were specified
     arg_names <- get_full_arg_names(independent_args, dependent_arg_function)
 
@@ -235,8 +245,19 @@ check_runner_results <- function(
         stop(msg)
     }
 
+    # Make sure each runner produces a data frame
+    is_df <- sapply(initial_runner_res, is.data.frame)
+
+    if (any(!is_df)) {
+        msg <- paste(
+            'Some runners did not produce data frames:',
+            paste(names(initial_runner_res)[!is_df], collapse = ', ')
+        )
+        stop(msg)
+    }
+
     # Make sure each runner produces the necessary columns in its output
-    expected_columns <- as.character(full_data_definitions)
+    expected_columns <- c('time', as.character(full_data_definitions))
 
     missing_columns <- lapply(initial_runner_res, function(res) {
         expected_columns[!expected_columns %in% colnames(res)]
@@ -402,12 +423,15 @@ check_obj_fun <- function(obj_fun, initial_ind_arg_values, verbose) {
     initial_error_terms <-
         obj_fun(as.numeric(initial_ind_arg_values), return_terms = TRUE)
 
+    initial_error <- sum(unlist(initial_error_terms))
+
     if (verbose) {
         cat('\nThe initial error metric terms:\n\n')
         utils::str(initial_error_terms)
-    }
 
-    initial_error <- sum(unlist(initial_error_terms))
+        cat('\nThe initial error metric value:\n\n')
+        print(initial_error)
+    }
 
     if (!is.finite(initial_error)) {
         stop(