Nothing
## Prompt user for confirmation to proceed
ask_for_confirmation = function(n_combinations, force) {
if (force) {
return()
}
if (!interactive()) { # nocov start
logger("Quitting run_parameter_scan() in non-interactive context with
`force == FALSE`.", level = INFO)
rlang::interrupt()
}
print(sprintf("Number of parameter combinations: %s", n_combinations))
response = prompt_user("Continue? [Y/n] ")
if (length(response) > 0 && tolower(response) %in% c("x", "n", "c", "q")) {
logger("Quitting as requested by user.", level = INFO)
rlang::interrupt()
} # nocov end
}
#' Parameter Scan
#'
#' Run ModVege for a different sets of parameters.
#'
#' @note
#' Special care has to be taken in the creation of the *param_values* argument.
#' It's possible to choose values that do not allow for any valid combination.
#' Confer [create_combinations()].
#'
#' @param environment Either a [ModvegeEnvironment] instance with all the site,
#' management and weather inputs expected by [ModvegeSite]`$run()` **or** a
#' string representing the name of a config file to read in order to
#' generate the [ModvegeEnvironment] with [read_config()]. Note that, in
#' the latter case, only the first found configuration is used if there are
#' more than one valid uncommented lines in the config file.
#' @param param_values A named list where each key stands for a ModVege
#' parameter, i.e. a member of [ModvegeParameters]`$parameter_names`. Each
#' list entry then has to be a vector containing the allowed values for the
#' respective parameter. All possible allowed combinations of these
#' parameter values are then generated and fed into a ModVege run.
#' @param force Boolean. By default (`force = FALSE`), the function first
#' counts the number of parameter combinations that need to be run and asks
#' the user, if it should proceed. This can be suppressed by letting `force
#' = TRUE`.
#' @param outfilename String. If nonempty, the results are stored as an `rds`
#' file with filename *outfilename* using the [saveRDS()] function.
#'
#' @return results A list containing an entry for each supplied parameter set
#' in *param_values*. Each entry is itself a list containing the following
#' keys:
#' \describe{
#' \item{params}{The parameter set that was used to run ModVege for this
#' entry.}
#' \item{data}{A list containing for each simulated year a [ModvegeSite]
#' object which was run for the respective year and therefore carries the
#' respective results.}
#' }
#'
#' @seealso [ModvegeParameters], [saveRDS()], [create_combinations()]
#'
#' @examples
#' env = create_example_environment()
#' # We're creating a trivial list of parameters to explore here in order to
#' # prevent the example from requiring a long time to execute. See
#' # [create_combinations()] for more realistic uses of param_values.
#' param_values = list(w_FGA = c(0, 1), w_FGB = c(0, 1))
#' run_parameter_scan(env, param_values, force = TRUE)
#'
#' @md
#' @export
run_parameter_scan = function(environment, param_values, force = FALSE,
outfilename = "") {
start_time = Sys.time()
# Read environment from config file, if *environment* is a string.
if (is.character(environment)) {
# Only take first run from whatever config is specified
environment = read_config(environment)[[1]]
}
# Ensure all supplied parameters are valid
environment$parameters$check_parameters(names(param_values),
check_for_completeness = FALSE)
# Create all parameter combinations
parameter_sets = create_combinations(param_values)
n_combinations = length(parameter_sets)
ask_for_confirmation(n_combinations, force)
param_names = names(parameter_sets[[1]])
results = list()
for (i in 1:n_combinations) {
logger(sprintf("Parameter combination %i out of %i.", i, n_combinations),
level = INFO
)
environment$parameters$set_parameters(parameter_sets[[i]])
mvs = growR_run_loop(c(environment), output_dir = "")
results[[i]] = list(params = parameter_sets[[i]],
data = mvs[[1]])
}
logger("Completed parameter scan.", level = INFO)
logger(paste0("Time used: ", Sys.time() - start_time), level = INFO)
if (outfilename != "") { # nocov start
question = "About to write results to `%s`. Continue? [Y/n] "
response = prompt_user(sprintf(question, outfilename))
if (!response %in% c("y", "Y", "")) {
logger(sprintf("Storing results as `%s`.", outfilename))
}
saveRDS(results, outfilename)
} # nocov end
return(results)
}
#' Analyze results of a parameter scan
#'
#' @param parameter_scan_results String or List. If a string, it is
#' interpreted as the name of a `rds` file that contains the results of a
#' parameter scan which is then loaded using [readRDS()]. Otherwise, it
#' should be the output of [run_parameter_scan()] directly.
#' @param datafile Name or path to a file containing measured data. The model
#' results in *parameter_scan_results* are compared to the data therein. If
#' empty, the site is inferred from the [ModvegeSite] objects in
#' *parameter_scan_results* and a corresponding data file is searched for
#' in `getOption("growR.data_dir", default = "data").
#' @param smooth_interval Int. Number of days over which the variable `dBM`
#' is smoothened. Should be set to make experimental data and simulated
#' data to be as comparable as possible.
#' @return analyzed A list with five keys: `dBM`, `cBM`, `cBM_end`, `metrics`
#' and `params`.
#' \describe{
#' \item{dBM}{A data.frame with `1 + n_params + n_metrics` columns
#' where each row represents a different parameter combination.
#' The first column (`n`) gives the row number and is used to identify a
#' parameter combination. The subsequent `n_params` columns give the
#' values of the parameters used in this combination. The final `n_metrics`
#' columns give the resulting performance score of the model run with
#' these parameters for each metric applied to model variable `dBM`.}
#' \item{cBM}{A data.frame of same format as for the key *dBM*. The first
#' `n_params + 1` columns are identical to the data.frame in *dBM*. The
#' difference is that the final `n_metrics` columns give performance
#' scores with respect to the model variable `cBM`.}
#' \item{cBM_end}{A data.frame analogous to *dBM* and *cBM*, only this time
#' the last `n_metrics` columns give performance scores with respect to
#' the variable `cBM_end`, which is the final value of `cBM`, i.e. the
#' cumulative grown biomass at the end of the year.}
#' \item{params}{A vector containing the names of the scanned parameters.
#' These are also the column names of columns `2:(n_params+1)` in
#' *results*.}
#' \item{metrics}{A vector containing the names of the employed
#' performance metrics. These are also the column names of the last
#' `n_metrics` columns in *results*.}
#' }
#'
#' @seealso [run_parameter_scan()], [readRDS()]
#'
#' @examples
#' # There needs to be data available with which the model is to be compared.
#' # For this example, use data provided by the package.
#' path = system.file("extdata", package = "growR")
#' datafile = file.path(path, "posieux1.csv")
#'
#' # We also use example parameter scan data provided by the package.
#' # In the real world, you would generally create your own data using
#' # `run_parameter_scan()`.
#' analyze_parameter_scan(parameter_scan_example, datafile = datafile)
#'
#' @md
#' @export
analyze_parameter_scan = function(parameter_scan_results, datafile = "",
smooth_interval = 28) {
supported_variables = c("dBM", "cBM", "cBM_end")
metrics_to_use = c("bias", "MAE", "RMSE")
# Read results from .Rds file, if applicable
if (is.character(parameter_scan_results)) {
results = readRDS(parameter_scan_results)
} else {
results = parameter_scan_results
}
n_combinations = length(results)
first = results[[1]]$data
years = sapply(first, function(mv) {mv$year})
n_years = length(years)
# Load measured data for detected site
if (datafile == "") {
site = results[[1]]$data[[1]]$site_name
measured_data = load_data_for_sites(c(site))[[1]]
} else {
measured_data = load_measured_data(datafile)[[1]]
}
# Reduce to relevant years
relevant_data = measured_data[measured_data$year %in% years, ]
# Calculate measured cBM_end
measured_cBM_ends = c()
for (year in years) {
mask2 = measured_data$year == year
measured_cBM_end = sum(measured_data[mask2, ]$cBM)
measured_cBM_ends = c(measured_cBM_ends, measured_cBM_end)
}
# Assemble metrics.
for (combination in 1:n_combinations) {
modvegesites = results[[combination]]$data
# Collect year curves into one vector and smoothen dBM
dBM = c()
cBM = c()
cBM_end = c()
for (i_year in 1:n_years) {
mv = modvegesites[[i_year]]
mask = relevant_data$DOY[relevant_data$year == years[i_year]]
cBM = c(cBM, mv$cBM[mask])
cBM_end = c(cBM_end, mv$cBM[length(mv$cBM)])
smoothed = box_smooth(mv$dBM, box_width = smooth_interval)
dBM = c(dBM, smoothed[mask])
}
results[[combination]][["cBM"]] = list()
results[[combination]][["dBM"]] = list()
results[[combination]][["cBM_end"]] = list()
# Employ performance metrics
for (metric in metrics_to_use) {
# cBM
m_cBM = metric_map[[metric]][["func"]](cBM, relevant_data$cBM)
results[[combination]][["cBM"]][[metric]] = m_cBM
# dBM
m_dBM = metric_map[[metric]][["func"]](dBM, relevant_data$dBM)
results[[combination]][["dBM"]][[metric]] = m_dBM
# cBM_end
m_cBM_end = metric_map[[metric]][["func"]](cBM_end, measured_cBM_ends)
results[[combination]][["cBM_end"]][[metric]] = m_cBM_end
}
}
# Reformat results to a data.frame
params = names(results[[1]][["params"]])
return_list = list(metrics = metrics_to_use, params = params)
for (variable in supported_variables) {
analyzed = data.frame(n = 1:n_combinations)
for (i_param in 1:length(params)) {
param_name = params[[i_param]]
analyzed[[param_name]] = extract_results(results, "params", i_param)
}
for (i_metric in 1:length(metrics_to_use)) {
metric_name = metrics_to_use[[i_metric]]
analyzed[[metric_name]] = extract_results(results, variable, i_metric)
}
return_list[[variable]] = analyzed
}
return(return_list)
}
extract_results = function(results, key, index) {
sapply(results, function(r) {r[[key]][[index]]})
}
#' Plot Parameter Scan Results
#'
#' @description
#' This class facilitates interactive visual analysis of parameter scan
#' results.
#'
#' @seealso [plot_parameter_scan()]
#'
#' @md
#' @export
PscanPlotter = R6Class(
"PscanPlotter",
public = list(
#-Public-fields---------------------------------------------------------------
#' @field analyzed List, as output by [analyze_parameter_scan()].
analyzed = NULL,
#' @field params Vector of names of scanned parameters.
params = NULL,
#' @field metrics Vector of names of model performance metrics to use.
metrics = NULL,
#' @field n_params Number of scanned parameters.
n_params = NULL,
#' @field n_metrics Number of performance metrics to apply.
n_metrics = NULL,
#' @field res data.frame holding parameter scan results. It contains
#' `n_params + n_metrics + 1` columns: one column for each scanned
#' parameter, one for each employed metric and an additional column
#' (name `n`) to give each parameter combination (i.e. each row) an
#' identifying number.
res = NULL,
#' @field n_combinations Number of rowns in `res`.
n_combinations = NULL,
#' @field sorted List containing copies of `res` which are each sorted by
#' a different performance metric. List keys are the values in
#' `self$metrics`.
sorted = NULL,
#' @field selection Vector of integers corresponding to the ID number of
#' combinations (column `n` in `self$res`) that is to be highlighted.
selection = NULL,
#-Public-methods--------------------------------------------------------------
#' @description
#' Construct and set up a [PscanPlotter] instance.
#'
#' @param analyzed List; Output of [analyze_parameter_scan()].
#' @param variable Str; Name of variable in *analyzed* to visualize. Can
#' be changed later with `set_variable()`. Allowed values are the keys
#' in *analyzed* except for `params` and `metrics`.
#'
#' @seealso [analyze_parameter_scan()]
#'
#' @md
initialize = function(analyzed, variable = "dBM") {
self$analyzed = analyzed
# Get the number of metrics and the number of parameters from the data
self$params = analyzed[["params"]]
self$n_params = length(self$params)
self$metrics = analyzed[["metrics"]]
self$n_metrics = length(self$metrics)
self$set_variable(variable)
},
#' @description Choose which variable to visualize.
#'
#' @param variable Chosen variable name. One of "dBM", "cBM", "cBM_end"
set_variable = function(variable) {
choices = private$get_variable_options()
if (!variable %in% choices) {
msg = "Invalid variable `%s`. Choose one of: `%s`."
stop(sprintf(msg, variable, paste(choices, collapse = "`, `")))
}
private$variable = variable
self$res = self$analyzed[[variable]]
private$sort_results()
},
#' @description Enter analysis loop.
#'
#' This plots the analysis results and enters a simple command-line
#' interface through which more insights can be gathered. Particularly,
#' it allows highlighting specific parameter combinations, either by
#' their index number or by selecting the best performers according to a
#' given metric.
#'
#' @seealso `PscanPlotter$plot()`
#' @md
analyze = function() {
private$quit_requested = FALSE
self$plot()
while (!private$quit_requested) {
private$ask_input()
}
cat("[PSP]Goodbye!\n")
},
#' @description Plot parameter scan results.
#'
#' For every combination of scanned parameter and metric, a subplot is
#' generated in which the parameter values are plotted against
#' performance score in that metric for every parameter combination.
#'
#' The result of this is static. Use `PscanPlotter$analyze()`
#' for an interactive version.
#'
#' @md
plot = function() {
# Create new graphics device
if (is.null(private$device_number)) {
# dev.new()
device_number = dev.cur()
private$device_number = ifelse(device_number == 1, 2, device_number)
} else {
dev.set(private$device_number)
}
# Set up the plot window
oldpar = private$set_par()
on.exit(par(oldpar))
for (i_metric in 1:self$n_metrics) {
for (i_param in 1:self$n_params) {
private$draw_points(i_metric, i_param)
}
}
# Force creation of Rplots.pdf when run from Rscript
if (!interactive()) {
dev.off()
}
},
#' @description Print information on selected parameter combinations.
#'
#' The parameter values and performance scores of all combinations
#' referred to by the integers in *selection* are printed to console.
#'
#' @param selection Vector of integers representing IDs of parameter
#' combinations (i.e. column `n` in `self$res`).
#'
print_info = function(selection) {
cat(sprintf("Performance scores for variable `%s`.\n", private$variable))
ncol = length(self$res)
base_string = paste("|", rep("%6.5s", ncol), collapse = " ")
base_string = paste(base_string, "|\n")
cat(do.call(sprintf, c(base_string, as.list(colnames(self$res)))))
for (s in selection) {
values = self$res[s, ]
cat(do.call(sprintf, c(base_string, as.list(c(values)))))
}
}
),
private = list(
#-Private-fields--------------------------------------------------------------
quit_requested = FALSE,
base_prompt = "\n[PSP] %s \n>>> ",
base_instructions = "Select combinations to highlight or enter command.
`?` for help. `q` to quit.",
help_message = "You can highlight a specific parameter combination and
print its parameter set and performance scores by issuing the index of the
combination or a comma-separated list of indices.
The following additional commands are available:
? Show this help message.
q Quit the parameter scan analyzer session.
b Highlight best performers. You will be queried to select a metric.
The 5 best scoring parameter combinations for the selected metric will
then be highlighted.
v Select model variable for which performance is to be evaluated.
",
## @field variable String representing the key in *self.analyzed* for
## which performance is plotted.
variable = NULL,
device_number = NULL,
#-Private-methods-------------------------------------------------------------
## Update self$sorted
sort_results = function() {
# Create different orderings for every metric
sorted = list()
for (metric in self$metrics) {
performance = abs(self$res[[metric]] - metric_map[[metric]]$target)
sorted[[metric]] = self$res[order(performance), ]
}
self$sorted = sorted
},
## Prompt user for input at the "main" level.
ask_input = function() {
input = prompt_user(sprintf(private$base_prompt,
private$base_instructions))
# Quit
if (input == "q") {
private$quit_requested = TRUE
return()
# Help
} else if (input == "?") {
cat(private$help_message)
# Standard: highlight specific combinations
} else if (private$parse_numbers(input)) {
self$plot()
self$print_info(self$selection)
# Highlight *n* best performers
} else if (input == "b") {
private$query_best_performers()
} else if (input == "v") {
private$query_variable()
} else {
private$input_not_understood(input)
}
},
## Ask for user input for displaying best scoring combinations.
query_best_performers = function() {
# First, ask user to select metric
msg = "Select metric for which to highlight best performers:"
input = menu(self$metrics, title = msg)
if (input == 0) {
return()
} else {
metric = self$metrics[[input]]
}
## Continue by querying number of points to highlight
n_to_highlight = 5
best_performers = self$sorted[[metric]][1:n_to_highlight,]
self$selection = best_performers$n
# Go!
self$plot()
self$print_info(self$selection)
},
## Return the allowed variables, i.e. those variable names that are
## present in self$analyzed.
get_variable_options = function() {
return(setdiff(names(self$analyzed), c("params", "metrics")))
},
query_variable = function() {
msg = "Select variable for which to show performance:"
choices = private$get_variable_options()
selection = menu(choices, title = msg)
if (selection == 0) {
return()
} else {
# Input understood as integer
variable = choices[selection]
}
self$set_variable(variable)
self$plot()
},
## Try parsing a user input of potentially comma-separated numbers
parse_numbers = function(input) {
# Check for comma
if (grepl(",", input)) {
numbers = strsplit(input, ",")[[1]]
self$selection = as.integer(numbers)
return(TRUE)
} else {
number = as.integer(input)
if (is.na(number)) {
return(FALSE)
} else {
self$selection = number
return(TRUE)
}
}
},
## Handle unrecognized user input.
input_not_understood = function(input) {
cat(sprintf("Your input `%s` was not understood.\n", input))
},
## Plot points into subplot at *row*, *col*.
## Rows correspond to metrics and columns to scan parameters.
draw_points = function(row, col) {
# Get data to plot
param_name = self$params[[col]]
metric_name = self$metrics[[row]]
param_values = self$res[[param_name]]
metric_values = self$res[[metric_name]]
# Select correct plot
# graphics::par(mfg = c(row, col))
graphics::plot(param_values, metric_values,
xlab = "", ylab = "",
ylim = metric_map[[metric_name]][["limits"]],
axes = FALSE,
frame.plot = TRUE
)
# Highlight points
graphics::lines(param_values[self$selection],
metric_values[self$selection],
type = "p", col = "red")
# Add ticks with or without labels
xticklabels = ifelse(row == self$n_metrics, TRUE, FALSE)
yticklabels = ifelse(col == 1, TRUE, FALSE)
graphics::Axis(side = 1, labels = xticklabels)
graphics::Axis(side = 2, labels = yticklabels)
if (row == 1 & col == 1) { #as.integer(self$n_params/2 + 0.5)) {
## Global title
graphics::title(private$variable, adj = 0)
}
if (row == self$n_metrics) {
graphics::mtext(param_name, side = 1, line = 3)
}
if (col == 1) {
graphics::mtext(metric_name, side = 2, line = 3)#, outer = TRUE
}
},
## Set up plot par() and return original par() for resetting.
set_par = function() {
oldpar = par(no.readonly = TRUE)
graphics::par(mfrow = c(self$n_metrics, self$n_params),
mar = c(1.0, 1.0, 1.0, 1.0),
oma = c(4, 4, 1, 1)
)
return(oldpar)
}
)
)
#' Plot Parameter Scan Results
#'
#' Visualize the results of a parameter scan and allow interactive inspection
#' of model performance in parameter space.
#'
#' @details
#' Under the hood this function just creates a [PscanPlotter] object and calls its
#' `analyze` method.
#'
#' @param analyzed List; Output of [analyze_parameter_scan()].
#' @param variable Str; Name of variable in *analyzed* to visualize. Can
#' be changed later with `PscanPlotter$set_variable()`. Allowed values are
#' the keys in *analyzed* except for `params` and `metrics`.
#' @param interactive boolean; Toggle between just creating a static plot
#' (`interactive = FALSE`) or entering a small, interactive analysis
#' setting (`interactive = TRUE`, default).
#' @return A [PscanPlotter] object.
#'
#' @seealso [analyze_parameter_scan()], [PscanPlotter]`$analyze()`
#'
#' @examplesIf interactive()
#' # There needs to be data available with which the modle is to be compared.
#' # For this example, use data provided by the package.
#' path = system.file("extdata", package = "growR")
#' datafile = file.path(path, "posieux1.csv")
#'
#' # Analyze example output of `run_parameter_scan()`.
#' results = analyze_parameter_scan(parameter_scan_example, datafile = datafile)
#' # The following plots the results.
#' psp = plot_parameter_scan(results, interactive = FALSE)
#'
#' # The interactive session can still be entered later from the returned
#' # PscanPlotter object
#' psp$analyze()
#'
#' @md
#' @export
plot_parameter_scan = function(analyzed, variable = "dBM", interactive = TRUE) {
PSP = PscanPlotter$new(analyzed, variable = variable)
if (interactive) {
PSP$analyze()
} else {
PSP$plot()
}
return(PSP)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.