From 0a5709027fae8b9363629a0ebc31440a97582eb0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Vezy?= Date: Wed, 11 Dec 2024 11:03:34 +0100 Subject: [PATCH 1/4] Create flint.yaml --- .github/workflows/flint.yaml | 35 +++++++++++++++++++++++++++++++++++ 1 file changed, 35 insertions(+) create mode 100644 .github/workflows/flint.yaml diff --git a/.github/workflows/flint.yaml b/.github/workflows/flint.yaml new file mode 100644 index 0000000..871d84c --- /dev/null +++ b/.github/workflows/flint.yaml @@ -0,0 +1,35 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + release: + types: [published] + workflow_dispatch: + +name: flint + +jobs: + flint: + runs-on: macOS-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: flint-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-r@v2 + + - name: Install flint + run: install.packages("flint", repos = c("https://etiennebacher.r-universe.dev/", getOption("repos"))) + shell: Rscript {0} + + - name: Run flint + run: flint::lint() + shell: Rscript {0} From 1df18563a4a8736e27a015b0adf549611b80dc65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Vezy?= Date: Wed, 11 Dec 2024 14:20:17 +0100 Subject: [PATCH 2/4] Create style.yaml --- .github/workflows/style.yaml | 85 ++++++++++++++++++++++++++++++++++++ 1 file changed, 85 insertions(+) create mode 100644 .github/workflows/style.yaml diff --git a/.github/workflows/style.yaml b/.github/workflows/style.yaml new file mode 100644 index 0000000..0fd50fe --- /dev/null +++ b/.github/workflows/style.yaml @@ -0,0 +1,85 @@ +# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples +# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help +on: + pull_request: + paths: + [ + "**.[rR]", + "**.[qrR]md", + "**.[rR]markdown", + "**.[rR]nw", + "**.[rR]profile", + ] + workflow_dispatch: + +name: style + +permissions: read-all + +jobs: + style: + runs-on: ubuntu-latest + permissions: + contents: write + env: + GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + steps: + - name: Checkout repo + uses: actions/checkout@v4 + with: + fetch-depth: 0 + + - name: Setup R + uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - name: Install dependencies + uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::styler, any::roxygen2 + needs: styler + + - name: Enable styler cache + run: styler::cache_activate() + shell: Rscript {0} + + - name: Determine cache location + id: styler-location + run: | + cat( + "location=", + styler::cache_info(format = "tabular")$location, + "\n", + file = Sys.getenv("GITHUB_OUTPUT"), + append = TRUE, + sep = "" + ) + shell: Rscript {0} + + - name: Cache styler + uses: actions/cache@v4 + with: + path: ${{ steps.styler-location.outputs.location }} + key: ${{ runner.os }}-styler-${{ github.sha }} + restore-keys: | + ${{ runner.os }}-styler- + ${{ runner.os }}- + + - name: Style + run: styler::style_pkg() + shell: Rscript {0} + + - name: Commit and push changes + run: | + if FILES_TO_COMMIT=($(git diff-index --name-only ${{ github.sha }} \ + | egrep --ignore-case '\.(R|[qR]md|Rmarkdown|Rnw|Rprofile)$')) + then + git config --local user.name "$GITHUB_ACTOR" + git config --local user.email "$GITHUB_ACTOR@users.noreply.github.com" + git commit ${FILES_TO_COMMIT[*]} -m "Style code (GHA)" + git pull --ff-only + git push origin + else + echo "No changes to commit." + fi From 928fd90dd6a3b44dd60d88a0d8f97118a757e1c7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Vezy?= Date: Wed, 11 Dec 2024 14:20:32 +0100 Subject: [PATCH 3/4] Delete flint.yaml --- .github/workflows/flint.yaml | 35 ----------------------------------- 1 file changed, 35 deletions(-) delete mode 100644 .github/workflows/flint.yaml diff --git a/.github/workflows/flint.yaml b/.github/workflows/flint.yaml deleted file mode 100644 index 871d84c..0000000 --- a/.github/workflows/flint.yaml +++ /dev/null @@ -1,35 +0,0 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help -on: - push: - branches: [main, master] - pull_request: - branches: [main, master] - release: - types: [published] - workflow_dispatch: - -name: flint - -jobs: - flint: - runs-on: macOS-latest - # Only restrict concurrency for non-PR jobs - concurrency: - group: flint-${{ github.event_name != 'pull_request' || github.run_id }} - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - permissions: - contents: write - steps: - - uses: actions/checkout@v4 - - - uses: r-lib/actions/setup-r@v2 - - - name: Install flint - run: install.packages("flint", repos = c("https://etiennebacher.r-universe.dev/", getOption("repos"))) - shell: Rscript {0} - - - name: Run flint - run: flint::lint() - shell: Rscript {0} From 1fc867a8e1f0c51dd0adfa5059c2252af01cddf3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?R=C3=A9mi=20Vezy?= Date: Wed, 11 Dec 2024 14:20:51 +0100 Subject: [PATCH 4/4] Style pkg --- R/FwdRegAgMIP.R | 8 +- R/bayesian_functions.R | 12 +- R/compute_eq_const.R | 33 +++--- R/estim_param.R | 26 ++-- R/filter_obs.R | 1 - R/frequentist_functions.R | 28 +++-- R/intersect_sim_obs.R | 5 +- R/is_data.R | 14 +-- R/ls_criteria.R | 25 ++-- R/main_crit.R | 38 +++--- R/obsSim_consistency.R | 63 +++++----- R/optim_switch.R | 6 +- R/param_info_functions.R | 3 +- R/test_wrapper.R | 1 - tests/spelling.R | 9 +- tests/testthat/test-compute_eq_const.R | 16 ++- tests/testthat/test-filter_obs.R | 1 - tests/testthat/test-filter_params.R | 18 ++- tests/testthat/test-get_params.R | 28 +++-- tests/testthat/test-init_values_functions.R | 29 +++-- tests/testthat/test-intersect_sim_obs.R | 12 +- tests/testthat/test-isdata.R | 21 ++-- tests/testthat/test-likelihood.R | 16 ++- tests/testthat/test-ls_criteria.R | 54 ++++++--- tests/testthat/test-obsSim_consistency.R | 65 ++++++---- tests/testthat/test-optim_switch.R | 8 +- tests/testthat/test-sample_params.R | 12 +- .../AgMIP_Calibration_Phenology_protocol.Rmd | 109 ++++++++++------- ...psimX_parameter_estimation_simple_case.Rmd | 111 ++++++++++-------- vignettes/CroptimizR.Rmd | 4 +- vignettes/Designing_a_model_wrapper.Rmd | 82 ++++++------- vignettes/Parameter_estimation_DREAM.Rmd | 30 +++-- ...meter_estimation_Specific_and_Varietal.Rmd | 27 ++--- .../Parameter_estimation_simple_case.Rmd | 33 ++---- 34 files changed, 544 insertions(+), 404 deletions(-) diff --git a/R/FwdRegAgMIP.R b/R/FwdRegAgMIP.R index 1c73759..e69acdd 100644 --- a/R/FwdRegAgMIP.R +++ b/R/FwdRegAgMIP.R @@ -47,7 +47,6 @@ select_param_FwdRegAgMIP <- function(oblig_param_list, add_param_list, crt_list, return(res) } else if (crt_list[length(crt_list)] == add_param_list[length(add_param_list)]) { - # we tested all parameters if (crt_info_crit < min(prev_info_crit)) { res$selected <- TRUE @@ -56,7 +55,6 @@ select_param_FwdRegAgMIP <- function(oblig_param_list, add_param_list, crt_list, } return(res) } else if (length(crt_list) == length(oblig_param_list)) { - # we only tested so far the obligatory parameters res$selected <- TRUE res$next_candidates <- c(oblig_param_list, add_param_list[1]) @@ -69,7 +67,6 @@ select_param_FwdRegAgMIP <- function(oblig_param_list, add_param_list, crt_list, add_param_list[which(add_param_list == crt_list[length(crt_list)]) + 1] ) } else { - # Replace the last candidate parameter by the next candidate res$selected <- FALSE res$next_candidates <- c( @@ -118,7 +115,8 @@ post_treat_FwdRegAgMIP <- function(optim_results, crit_options, crt_list, ## Store the results per step v_init <- as.vector( - t(optim_results$init_values[optim_results$ind_min_crit, ])) + t(optim_results$init_values[optim_results$ind_min_crit, ]) + ) names(v_init) <- names(optim_results$init_values) info_new_step <- setNames( tibble::tibble( @@ -138,7 +136,7 @@ post_treat_FwdRegAgMIP <- function(optim_results, crit_options, crt_list, info_crit_func()$name, "Selected step" ) ) -param_selection_steps <- dplyr::bind_rows(param_selection_steps, info_new_step) + param_selection_steps <- dplyr::bind_rows(param_selection_steps, info_new_step) ind_min_infocrit <- which.min(param_selection_steps[[info_crit_func()$name]]) param_selection_steps[, "Selected step"] <- "" param_selection_steps[ind_min_infocrit, "Selected step"] <- "X" diff --git a/R/bayesian_functions.R b/R/bayesian_functions.R index f73d3d0..cab464c 100644 --- a/R/bayesian_functions.R +++ b/R/bayesian_functions.R @@ -52,7 +52,8 @@ plot_bayesian <- function(optim_options, param_info, optim_results) { nb_chains <- length(out$chain) nb_iterations <- nrow(optim_results$post_sample) / nb_chains - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "iterAndDensityPlots.pdf"), width = 9, height = 9 @@ -81,7 +82,8 @@ plot_bayesian <- function(optim_options, param_info, optim_results) { } ) - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "marginalPlots.pdf"), width = 9, height = 9 @@ -110,7 +112,8 @@ plot_bayesian <- function(optim_options, param_info, optim_results) { ) if (nb_params >= 2) { - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "correlationPlots.pdf"), width = 9, height = 9 @@ -145,7 +148,8 @@ plot_bayesian <- function(optim_options, param_info, optim_results) { # an error if (is.null(optim_options$thin)) optim_options$thin <- 1 if (nb_iterations >= (optim_options$thin + 50)) { - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "gelmanDiagPlots.pdf"), width = 9, height = 9 diff --git a/R/compute_eq_const.R b/R/compute_eq_const.R index ef7d9f9..ed3ec3a 100644 --- a/R/compute_eq_const.R +++ b/R/compute_eq_const.R @@ -11,47 +11,48 @@ #' @keywords internal #' compute_eq_const <- function(forced_param_values, param_values) { - comp_forced_values <- NULL is_vector <- is.vector(param_values) if (!is.null(forced_param_values)) { - param_values <- tibble::tibble(!!!param_values) param_values$situation <- NULL - nrows <- max(1,seq_len(nrow(param_values))) - comp_forced_values <- matrix(ncol = length(forced_param_values), - nrow = nrows) + nrows <- max(1, seq_len(nrow(param_values))) + comp_forced_values <- matrix( + ncol = length(forced_param_values), + nrow = nrows + ) colnames(comp_forced_values) <- names(forced_param_values) # Backticks are added here and in the following to handle parameters names # including special characters expr_ls <- - lapply(names(forced_param_values), function(x) paste0("`",x,"`","<-", - forced_param_values[[x]])) + lapply(names(forced_param_values), function(x) { + paste0( + "`", x, "`", "<-", + forced_param_values[[x]] + ) + }) names(expr_ls) <- names(forced_param_values) for (irow in 1:nrows) { - for (par in names(param_values)) { - eval(parse(text = paste0("`",par,"`","<-",param_values[[irow, par]]))) + eval(parse(text = paste0("`", par, "`", "<-", param_values[[irow, par]]))) } for (par in names(forced_param_values)) { eval(parse(text = expr_ls[[par]])) - eval(parse(text = paste0("comp_forced_values[irow,\"",par,"\"] <- ", - "`",par,"`"))) + eval(parse(text = paste0( + "comp_forced_values[irow,\"", par, "\"] <- ", + "`", par, "`" + ))) } - } if (is_vector) { - comp_forced_values <- comp_forced_values[1,] + comp_forced_values <- comp_forced_values[1, ] } else { comp_forced_values <- tibble::as_tibble(comp_forced_values) } - } return(comp_forced_values) - } - diff --git a/R/estim_param.R b/R/estim_param.R index 13e2f3c..36afbbf 100644 --- a/R/estim_param.R +++ b/R/estim_param.R @@ -198,9 +198,8 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, CroptimizR::BIC, CroptimizR::AICc, CroptimizR::AIC ), - weight=NULL, + weight = NULL, var_names = lifecycle::deprecated()) { - # Managing parameter names changes between versions: if (rlang::has_name(optim_options, "path_results")) { lifecycle::deprecate_warn("0.5.0", "estim_param(optim_options = 'is deprecated, use `out_dir` instead of `path_results`')") @@ -299,12 +298,14 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, "The following parameters are defined both in forced_param_values and param_info arguments of estim_param function while they should not (a parameter cannot be both forced and estimated except if it is part of the `candidate` parameters):", - paste(tmp,collapse = ","), + paste(tmp, collapse = ","), "\n They will be removed from forced_param_values." ) forced_param_values <- - forced_param_values[setdiff(names(forced_param_values), - setdiff(param_names, candidate_param))] + forced_param_values[setdiff( + names(forced_param_values), + setdiff(param_names, candidate_param) + )] } } @@ -360,8 +361,8 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, crt_candidates <- oblig_param_list if (length(crt_candidates) == 0) crt_candidates <- candidate_param[[1]] # in case there are only candidates ... count <- 1 - param_selection_steps<-NULL - tmp <- optim_switch(optim_method=optim_method,optim_options=optim_options) + param_selection_steps <- NULL + tmp <- optim_switch(optim_method = optim_method, optim_options = optim_options) # Parameter selection loop while (!is.null(crt_candidates)) { @@ -422,7 +423,7 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, forced_param_values = forced_param_values_tmp, info_level = info_level, info_crit_list = info_crit_list, - weight=weight + weight = weight ) ## Run the estimation @@ -432,14 +433,13 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, ) ## In case no results, there was an error during the estimation process => stop - if (length(res_tmp)==0) { + if (length(res_tmp) == 0) { stop("There was an error during the parameter estimation process. Please check warnings and messages displayed above and/or by running warnings().") } ## The following is done only if parameter selection is activated if (!is.null(candidate_param)) { - ### Update results in param_selection_steps param_selection_steps <- post_treat_FwdRegAgMIP( res_tmp, crit_options, @@ -465,8 +465,10 @@ estim_param <- function(obs_list, crit_function = crit_log_cwss, model_function, # Print and store results of parameter estimation steps if parameter selection was activated if (!is.null(candidate_param)) { - summary_FwdRegAgMIP(param_selection_steps, info_crit_list, path_results_ORI, - res) + summary_FwdRegAgMIP( + param_selection_steps, info_crit_list, path_results_ORI, + res + ) save_results_FwdRegAgMIP(param_selection_steps, path_results_ORI) res$param_selection_steps <- param_selection_steps } diff --git a/R/filter_obs.R b/R/filter_obs.R index 18c03db..86cdc77 100644 --- a/R/filter_obs.R +++ b/R/filter_obs.R @@ -41,7 +41,6 @@ filter_obs <- function(obs_list, var = NULL, situation = NULL, dates = NULL, include = FALSE, var_names = lifecycle::deprecated(), sit_names = lifecycle::deprecated()) { - # Managing parameter names changes between versions: if (lifecycle::is_present(sit_names)) { lifecycle::deprecate_warn("0.5.0", "filter_obs(sit_names)", "filter_obs(situation)") diff --git a/R/frequentist_functions.R b/R/frequentist_functions.R index af1c43d..65ca138 100644 --- a/R/frequentist_functions.R +++ b/R/frequentist_functions.R @@ -107,7 +107,8 @@ plot_frequentist <- function(optim_options, param_info, optim_results) { # EstimatedVSinit plot - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "EstimatedVSinit.pdf"), width = 9, height = 9 @@ -127,7 +128,8 @@ plot_frequentist <- function(optim_options, param_info, optim_results) { } ) - tryCatch({ + tryCatch( + { p <- plot_estimVSinit( init_values, est_values, crit_values, bounds$lb, bounds$ub @@ -155,7 +157,8 @@ plot_frequentist <- function(optim_options, param_info, optim_results) { # ValuesVSit plot - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "ValuesVSit.pdf"), width = 9, height = 9 @@ -184,7 +187,8 @@ plot_frequentist <- function(optim_options, param_info, optim_results) { # ValuesVSit_2D plot - tryCatch({ + tryCatch( + { grDevices::pdf( file = file.path(path_results, "ValuesVSit_2D.pdf"), width = 9, height = 9 @@ -284,7 +288,7 @@ plot_estimVSinit <- function(init_values, est_values, crit, lb, ub, theme(plot.title = element_text(hjust = 0.5)) if (bubble) { - p[[param_name]] <- p[[param_name]] + geom_point(alpha = 0.5, color = "red") + p[[param_name]] <- p[[param_name]] + geom_point(alpha = 0.5, color = "red") } p[[param_name]] <- p[[param_name]] + @@ -374,7 +378,7 @@ plot_valuesVSit <- function(df, param_info, iter_or_eval = c("iter", "eval"), if (all(df$crit > 0)) { trans <- "log10" mid <- (max(log10(df$crit)) - - min(log10(df$crit))) / 2 + min(log10(df$crit)) + min(log10(df$crit))) / 2 + min(log10(df$crit)) } else { warning("The criterion takes negative values, log transformation will not be done.") crit_log <- FALSE @@ -438,8 +442,10 @@ plot_valuesVSit <- function(df, param_info, iter_or_eval = c("iter", "eval"), color = "rep" )) + labs( - title = paste0("Evolution of the minimized criterion \n in function of the minimization ", - lab), + title = paste0( + "Evolution of the minimized criterion \n in function of the minimization ", + lab + ), y = "Minimized criterion", x = paste(lab, "number"), fill = "Repetition" @@ -526,7 +532,7 @@ plot_valuesVSit_2D <- function(df, param_info, iter_or_eval = c("eval", "iter"), if (all(df$crit > 0)) { trans <- "log10" mid <- (max(log10(df$crit)) - - min(log10(df$crit))) / 2 + min(log10(df$crit)) + min(log10(df$crit))) / 2 + min(log10(df$crit)) } else { warning("The criterion takes negative values, log transformation will not be done.") crit_log <- FALSE @@ -578,14 +584,14 @@ plot_valuesVSit_2D <- function(df, param_info, iter_or_eval = c("eval", "iter"), if (rep_label[1] == "begin_end" || rep_label[1] == "begin") { p[[ipair]] <- p[[ipair]] + geom_label(aes(label = rep), - data = filter(df, rep == irep) %>% filter(eval == min(.data$eval)), + data = filter(df, rep == irep) %>% filter(eval == min(.data$eval)), size = 3 ) } if (rep_label[1] == "begin_end" || rep_label[1] == "end") { p[[ipair]] <- p[[ipair]] + geom_label(aes(label = rep), - data = filter(df, rep == irep) %>% filter(eval == max(.data$eval)), + data = filter(df, rep == irep) %>% filter(eval == max(.data$eval)), size = 3 ) } diff --git a/R/intersect_sim_obs.R b/R/intersect_sim_obs.R index 8569eec..8c8cb66 100644 --- a/R/intersect_sim_obs.R +++ b/R/intersect_sim_obs.R @@ -9,7 +9,6 @@ #' @keywords internal #' intersect_sim_obs <- function(sim_list, obs_list) { - # Intersect situations situations <- intersect(names(sim_list), names(obs_list)) if (length(situations) == 0) { @@ -79,11 +78,11 @@ intersect_sim_obs <- function(sim_list, obs_list) { } sim_list <- sapply(situations, - function(x) sim_list[[x]][is.element(sim_list[[x]]$Date, list_dates[[x]]), ], + function(x) sim_list[[x]][is.element(sim_list[[x]]$Date, list_dates[[x]]), ], simplify = F ) obs_list <- sapply(situations, - function(x) obs_list[[x]][is.element(obs_list[[x]]$Date, list_dates[[x]]), ], + function(x) obs_list[[x]][is.element(obs_list[[x]]$Date, list_dates[[x]]), ], simplify = F ) diff --git a/R/is_data.R b/R/is_data.R index 6fbda8d..5b0f868 100644 --- a/R/is_data.R +++ b/R/is_data.R @@ -47,9 +47,8 @@ #' @keywords internal #' is.data <- function(data_list) { - # Check data_list format - if (length(data_list)==0 || !is.list(data_list) || !all(sapply( + if (length(data_list) == 0 || !is.list(data_list) || !all(sapply( data_list, function(x) is.data.frame(x) ))) { @@ -74,11 +73,12 @@ is.data <- function(data_list) { ))) { warning(paste( "Incorrect format, Date column include replicated dates for situations", - paste(names(data_list)[sapply( - data_list, - function(x) length(unique(x$Date)) < length(x$Date) - )], - collapse = "," + paste( + names(data_list)[sapply( + data_list, + function(x) length(unique(x$Date)) < length(x$Date) + )], + collapse = "," ) )) return(FALSE) diff --git a/R/ls_criteria.R b/R/ls_criteria.R index 1873cdf..cb4207b 100644 --- a/R/ls_criteria.R +++ b/R/ls_criteria.R @@ -102,8 +102,9 @@ crit_ols <- function(sim_list, obs_list) { #' @export #' @rdname ls_criteria crit_wls <- function(sim_list, obs_list, weight) { - - if(!nargs()) return("ls") # return criterion type (ls, log-ls, likelihood, log-likelihood) if no argument given + if (!nargs()) { + return("ls") + } # return criterion type (ls, log-ls, likelihood, log-likelihood) if no argument given var_list <- unique(unlist(lapply(obs_list, function(x) colnames(x)))) var_list <- setdiff(var_list, "Date") @@ -118,20 +119,24 @@ crit_wls <- function(sim_list, obs_list, weight) { res <- res[id_not_na] sigma <- weight(obs[id_not_na], var) - if (any(sigma==0)) { - stop(paste("Error in crit_wls: weight is zero for variable",var, - ". The wls criterion takes Inf value.", - "Please handle this case in the function given in weight argument of estim_param.")) + if (any(sigma == 0)) { + stop(paste( + "Error in crit_wls: weight is zero for variable", var, + ". The wls criterion takes Inf value.", + "Please handle this case in the function given in weight argument of estim_param." + )) } if (any(is.na(sigma))) { - stop(paste("Error in crit_wls: weight is NA for variable",var, - ". The wls criterion takes NA value.", - "Please handle this case in the function given in weight argument of estim_param.")) + stop(paste( + "Error in crit_wls: weight is NA for variable", var, + ". The wls criterion takes NA value.", + "Please handle this case in the function given in weight argument of estim_param." + )) } sz <- length(res) - result <- result + (res/sigma) %*% (res/sigma) + result <- result + (res / sigma) %*% (res / sigma) } return(as.numeric(result)) diff --git a/R/main_crit.R b/R/main_crit.R index 2d295ff..ae0ddd3 100644 --- a/R/main_crit.R +++ b/R/main_crit.R @@ -20,11 +20,13 @@ #' main_crit <- function(param_values, crit_options) { on.exit({ - if (exists("obs_sim_list") && !is.null(obs_sim_list$obs_list)) { .croptEnv$obs_var_list <- unique( - unlist(lapply(obs_sim_list$obs_list, - function(x) setdiff(names(x),c("Date"))))) + unlist(lapply( + obs_sim_list$obs_list, + function(x) setdiff(names(x), c("Date")) + )) + ) } if (crit_options$info_level >= 1 && !is.null(crit_options$tot_max_eval)) { @@ -47,7 +49,7 @@ main_crit <- function(param_values, crit_options) { .croptEnv$params_and_crit[[.croptEnv$eval_count]] <- dplyr::bind_cols( .croptEnv$params_and_crit[[.croptEnv$eval_count]], - forced_param_values=tibble::tibble(!!!forced_param_values) + forced_param_values = tibble::tibble(!!!forced_param_values) ) } @@ -56,7 +58,7 @@ main_crit <- function(param_values, crit_options) { ## should be changed for more robust test later ... if ((.croptEnv$eval_count == 1) || - (crit_options$irep > tail(.croptEnv$params_and_crit[[.croptEnv$eval_count - 1]]$rep,1))) { + (crit_options$irep > tail(.croptEnv$params_and_crit[[.croptEnv$eval_count - 1]]$rep, 1))) { eval <- 1 iter <- NA .croptEnv$last_iter <- 0 @@ -66,7 +68,7 @@ main_crit <- function(param_values, crit_options) { .croptEnv$last_iter <- iter } } else { - eval <- tail(.croptEnv$params_and_crit[[.croptEnv$eval_count - 1]]$eval,1) + 1 + eval <- tail(.croptEnv$params_and_crit[[.croptEnv$eval_count - 1]]$eval, 1) + 1 iter <- NA if (!is.na(crit) && (is.na(.croptEnv$last_crit) || crit < .croptEnv$last_crit)) { @@ -322,7 +324,7 @@ main_crit <- function(param_values, crit_options) { } if (!is.null(transform_var)) { model_results$sim_list <- lapply(model_results$sim_list, function(x) { - for (var in intersect(names(x),names(transform_var))) { + for (var in intersect(names(x), names(transform_var))) { x[var] <- transform_var[[var]](x[var]) } return(x) @@ -332,7 +334,7 @@ main_crit <- function(param_values, crit_options) { } # Check results, return NA if incorrect if (is.null(model_results) || - (!is.null(model_results$error) && model_results$error)) { + (!is.null(model_results$error) && model_results$error)) { warning("Error in transformation of simulation results.") return(crit <- NA) } @@ -365,7 +367,7 @@ main_crit <- function(param_values, crit_options) { } if (!is.null(transform_var)) { obs_list <- lapply(obs_list, function(x) { - for (var in intersect(names(x),names(transform_var))) { + for (var in intersect(names(x), names(transform_var))) { x[var] <- transform_var[[var]](x[var]) } return(x) @@ -395,8 +397,8 @@ main_crit <- function(param_values, crit_options) { # check presence of Inf/NA in simulated results where obs is not NA if (is_sim_inf_or_na(obs_sim_list$sim_list, obs_sim_list$obs_list, param_values)) { - warning("\nThe optimized criterion is set to NA.\n") - return(crit <- NA) + warning("\nThe optimized criterion is set to NA.\n") + return(crit <- NA) } # Filter reserved columns that should not be taken into account in the computation of the criterion @@ -425,16 +427,17 @@ main_crit <- function(param_values, crit_options) { ) # Compute criterion value - potential_arglist <- list(sim_list=obs_sim_list$sim_list, - obs_list=obs_sim_list$obs_list, - weight=crit_options$weight) + potential_arglist <- list( + sim_list = obs_sim_list$sim_list, + obs_list = obs_sim_list$obs_list, + weight = crit_options$weight + ) # Test weight function is well defined if (!is.null(crit_options$weight)) { - var_list_tmp <- names(obs_sim_list$obs_list[[1]]) var_tmp <- setdiff(var_list_tmp, "Date")[1] - obsvec_tmp <- obs_sim_list$obs_list[[1]][,var_tmp] + obsvec_tmp <- obs_sim_list$obs_list[[1]][, var_tmp] tryCatch( w <- crit_options$weight(na.omit(obsvec_tmp), var_tmp), error = function(cond) { @@ -449,12 +452,11 @@ main_crit <- function(param_values, crit_options) { stop("Caught an error while testing argument weight: \n it must be function that returns a numeric value (or vector of).") } - if (length(w)!=1 & length(w)!=length(na.omit(obsvec_tmp))) { + if (length(w) != 1 & length(w) != length(na.omit(obsvec_tmp))) { stop("Caught an error while testing argument weight: \n it must be a function that returns a single value or a vector of values of size the size of the vector of observed values given as first argument.") } - } crit <- do.call(crit_function, args = potential_arglist[names(formals(crit_function))]) diff --git a/R/obsSim_consistency.R b/R/obsSim_consistency.R index f2106ff..9557784 100644 --- a/R/obsSim_consistency.R +++ b/R/obsSim_consistency.R @@ -11,7 +11,6 @@ #' @keywords internal #' make_obsSim_consistent <- function(sim_list, obs_list) { - # Check homogeneity of types of obs and sim Date columns # ------------------------------------------------------ @@ -38,14 +37,14 @@ make_obsSim_consistent <- function(sim_list, obs_list) { if (!all(isPosixObs | isDateObs)) { stop(paste( - "Error: incorrect format for Date column of observation list for situations", + "Error: incorrect format for Date column of observation list for situations", paste(names(obs_list)[!(isPosixObs | isDateObs)], collapse = ", "), ". \n Date columns should be of type Date or POSIXct." )) } if (!all(isPosixSim | isDateSim)) { stop(paste( - "Error: incorrect format for Date column of simulation list for situations", + "Error: incorrect format for Date column of simulation list for situations", paste(names(sim_list)[!(isPosixSim | isDateSim)], collapse = ", "), ". \n Date columns should be of type Date or POSIXct." )) @@ -139,7 +138,6 @@ make_obsSim_consistent <- function(sim_list, obs_list) { #' @keywords internal #' check_obsSim_consistency <- function(sim_list, obs_list) { - # Check homogeneity of types of obs and sim Date columns # ------------------------------------------------------ @@ -162,14 +160,14 @@ check_obsSim_consistency <- function(sim_list, obs_list) { if (!all(isPosixObs | isDateObs)) { stop(paste( - "Error: incorrect format for Date column of observation list for situations", + "Error: incorrect format for Date column of observation list for situations", paste(names(obs_list)[!(isPosixObs | isDateObs)], collapse = ", "), ". \n Date columns should be of type Date or POSIXct." )) } if (!all(isPosixSim | isDateSim)) { stop(paste( - "Error: incorrect format for Date column of simulation list for situations", + "Error: incorrect format for Date column of simulation list for situations", paste(names(sim_list)[!(isPosixSim | isDateSim)], collapse = ", "), ". \n Date columns should be of type Date or POSIXct." )) @@ -226,21 +224,19 @@ check_obsSim_consistency <- function(sim_list, obs_list) { #' @keywords internal #' is_sim_inf_or_na <- function(sim_list, obs_list, param_values) { - res <- FALSE mess <- "" # check presence of Inf/NA in simulated results where obs is not NA for (sit in names(sim_list)) { - var_list <- lapply(names(sim_list[[sit]]), function(x) { - if (any(is.infinite(sim_list[[sit]][!is.na(obs_list[[sit]][,x]),][[x]])) || - any(is.na(sim_list[[sit]][!is.na(obs_list[[sit]][,x]),][[x]]))) { + if (any(is.infinite(sim_list[[sit]][!is.na(obs_list[[sit]][, x]), ][[x]])) || + any(is.na(sim_list[[sit]][!is.na(obs_list[[sit]][, x]), ][[x]]))) { return(list(sim_list[[sit]]$Date[(is.infinite(sim_list[[sit]][[x]]) & - !is.na(obs_list[[sit]][,x])) | - (is.na(sim_list[[sit]][[x]]) & - !is.na(obs_list[[sit]][,x]))])) - } else { + !is.na(obs_list[[sit]][, x])) | + (is.na(sim_list[[sit]][[x]]) & + !is.na(obs_list[[sit]][, x]))])) + } else { return(NULL) } }) @@ -248,27 +244,32 @@ is_sim_inf_or_na <- function(sim_list, obs_list, param_values) { names(var_list) <- names(sim_list[[sit]]) if (!is.null(unlist(var_list))) { - mess <- paste(mess, - paste("\n for situation ",sit, - ", variable(s):", - paste(sapply(names(var_list), function(x) { - if (!is.null(var_list[[x]])) { - paste("\n o",x," at date(s) ", - paste(var_list[[x]][[1]], collapse = ", ")) - } else { - return("") - } - }),collapse="") - )) + mess <- paste( + mess, + paste( + "\n for situation ", sit, + ", variable(s):", + paste(sapply(names(var_list), function(x) { + if (!is.null(var_list[[x]])) { + paste( + "\n o", x, " at date(s) ", + paste(var_list[[x]][[1]], collapse = ", ") + ) + } else { + return("") + } + }), collapse = "") + ) + ) res <- TRUE } - } if (res) { - warning("\n\nSimulated values are NA or Inf," ,mess, ",\nwhen using input parameters: ", - paste(names(param_values), collapse = ", "), ", with values: ", - paste(param_values, collapse = ", ")) + warning( + "\n\nSimulated values are NA or Inf,", mess, ",\nwhen using input parameters: ", + paste(names(param_values), collapse = ", "), ", with values: ", + paste(param_values, collapse = ", ") + ) } return(res) - } diff --git a/R/optim_switch.R b/R/optim_switch.R index fcc26aa..fc8010c 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -16,20 +16,16 @@ #' optim_switch <- function(...) { - # Store and save results at the end of the current optimization step if (nargs() > 2) { - # Initialize res res <- list() flag_error <- FALSE on.exit({ - - res$forced_param_values = crit_options$forced_param_values + res$forced_param_values <- crit_options$forced_param_values if (exists(".croptEnv")) { - # Save results even in case parameter estimation crash res$obs_var_list <- .croptEnv$obs_var_list if (exists(".croptEnv$obs_var_list")) { diff --git a/R/param_info_functions.R b/R/param_info_functions.R index d2415a0..da43940 100644 --- a/R/param_info_functions.R +++ b/R/param_info_functions.R @@ -384,7 +384,6 @@ get_init_values <- function(param_info) { # Case of simultaneous estimation of varietal and specific parameters } else if (is.list(param_info[[1]])) { - # check if colnames were set to params_names, if not set them # and handle translation if necessary for (i in 1:length(param_info)) { @@ -398,7 +397,7 @@ get_init_values <- function(param_info) { } } else { if (!is.null(param_info[[i]]$sit_list)) { - param_info[[i]]$init_values <- data.frame(t(rep(NA,length(param_info[[i]]$sit_list)))) + param_info[[i]]$init_values <- data.frame(t(rep(NA, length(param_info[[i]]$sit_list)))) } else { param_info[[i]]$init_values <- data.frame(NA) } diff --git a/R/test_wrapper.R b/R/test_wrapper.R index b4d1ec5..74e98d1 100644 --- a/R/test_wrapper.R +++ b/R/test_wrapper.R @@ -40,7 +40,6 @@ test_wrapper <- function(model_function, model_options, param_values, situation, var = NULL, sit_names = lifecycle::deprecated(), var_names = lifecycle::deprecated()) { - # Managing parameter names changes between versions: if (lifecycle::is_present(sit_names)) { lifecycle::deprecate_warn("0.5.0", "test_wrapper(sit_names)", "test_wrapper(situation)") diff --git a/tests/spelling.R b/tests/spelling.R index 33ef2c7..13f77d9 100644 --- a/tests/spelling.R +++ b/tests/spelling.R @@ -1,3 +1,6 @@ -if (requireNamespace("spelling", quietly = TRUE)) - spelling::spell_check_test(vignettes = TRUE, error = FALSE, - skip_on_cran = TRUE) +if (requireNamespace("spelling", quietly = TRUE)) { + spelling::spell_check_test( + vignettes = TRUE, error = FALSE, + skip_on_cran = TRUE + ) +} diff --git a/tests/testthat/test-compute_eq_const.R b/tests/testthat/test-compute_eq_const.R index 4512860..97bc41b 100644 --- a/tests/testthat/test-compute_eq_const.R +++ b/tests/testthat/test-compute_eq_const.R @@ -35,9 +35,15 @@ test_that("correct format", { test_that("correct content", { expect_equal(v_new_forced_param_values, c(p3 = 100, p4 = 1001)) expect_equal(v_new_forced_param_values, c(p3 = 100, p4 = 1001)) - expect_equal(t_new_forced_param_values, tibble::tibble(p3 = c(100, 100, 100), - p4 = c(1001, 2001, 3001))) - expect_equal(t_new_forced_param_values_v2, - tibble::tibble(p3 = c(100), - p4 = c(1001))) + expect_equal(t_new_forced_param_values, tibble::tibble( + p3 = c(100, 100, 100), + p4 = c(1001, 2001, 3001) + )) + expect_equal( + t_new_forced_param_values_v2, + tibble::tibble( + p3 = c(100), + p4 = c(1001) + ) + ) }) diff --git a/tests/testthat/test-filter_obs.R b/tests/testthat/test-filter_obs.R index 6ef101c..23c8651 100644 --- a/tests/testthat/test-filter_obs.R +++ b/tests/testthat/test-filter_obs.R @@ -39,7 +39,6 @@ test_that("filter_obs filters-out correctly", { filter_obs(obs_list, var = c("var2")), "No observations found in situation\\(s\\) sit2, sit3", ) - }) diff --git a/tests/testthat/test-filter_params.R b/tests/testthat/test-filter_params.R index 414e5f5..5654e9f 100644 --- a/tests/testthat/test-filter_params.R +++ b/tests/testthat/test-filter_params.R @@ -28,10 +28,16 @@ param_info2$durvieF <- list( param_info2_res <- param_info2["durvieF"] test_that("filter_param_info", { - expect_equal(eval( - parse(text = "CroptimizR:::filter_param_info(param_info1, \"durvieF\")")), - param_info1_res) - expect_equal(eval( - parse(text = "CroptimizR:::filter_param_info(param_info2, \"durvieF\")")), - param_info2_res) + expect_equal( + eval( + parse(text = "CroptimizR:::filter_param_info(param_info1, \"durvieF\")") + ), + param_info1_res + ) + expect_equal( + eval( + parse(text = "CroptimizR:::filter_param_info(param_info2, \"durvieF\")") + ), + param_info2_res + ) }) diff --git a/tests/testthat/test-get_params.R b/tests/testthat/test-get_params.R index 2f9071a..36eccdf 100644 --- a/tests/testthat/test-get_params.R +++ b/tests/testthat/test-get_params.R @@ -75,8 +75,8 @@ prior_3$durvieF <- list( ) ), init_values = data.frame(c(200, 300), c(250, 350)), - lb = 50, ub = 400 - ) + lb = 50, ub = 400 +) prior_4 <- list( init_values = c(dlaimax = 0.001), @@ -159,7 +159,8 @@ test_that("get_params_names", { ) expect_equal( eval( - parse(text = "CroptimizR:::get_params_names(prior_2, short_list=TRUE)")), + parse(text = "CroptimizR:::get_params_names(prior_2, short_list=TRUE)") + ), c("dlaimax", "durvieF") ) expect_equal( @@ -168,7 +169,8 @@ test_that("get_params_names", { ) expect_equal( eval( - parse(text = "CroptimizR:::get_params_names(prior_3, short_list=TRUE)")), + parse(text = "CroptimizR:::get_params_names(prior_3, short_list=TRUE)") + ), c("durvieF", "dlaimax") ) }) @@ -193,10 +195,16 @@ sg <- list( vec <- c(1, 2, 3, 4, 5) names(vec) <- eval(parse(text = "CroptimizR:::get_params_names(sg)")) test_that("get_params_per_sit", { - expect_equal(eval( - parse(text = "CroptimizR:::get_params_per_sit(sg, \"sit2\", vec)")), - c(p1 = 1, p2 = 3, p3 = 4)) - expect_equal(eval( - parse(text = "CroptimizR:::get_params_per_sit(sg, \"sit4\", vec)")), - c(p1 = 2, p2 = 3, p3 = 5)) + expect_equal( + eval( + parse(text = "CroptimizR:::get_params_per_sit(sg, \"sit2\", vec)") + ), + c(p1 = 1, p2 = 3, p3 = 4) + ) + expect_equal( + eval( + parse(text = "CroptimizR:::get_params_per_sit(sg, \"sit4\", vec)") + ), + c(p1 = 2, p2 = 3, p3 = 5) + ) }) diff --git a/tests/testthat/test-init_values_functions.R b/tests/testthat/test-init_values_functions.R index a86c54a..004a337 100644 --- a/tests/testthat/test-init_values_functions.R +++ b/tests/testthat/test-init_values_functions.R @@ -48,10 +48,13 @@ eval(parse(text = "init_values <- CroptimizR:::get_init_values(param_info)")) param_info_new <- within(param_info, rm(init_values)) eval( parse( -text = - "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)")) + text = + "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)" + ) +) eval(parse( - text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)")) + text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)" +)) test_that("set_init_values: simple case 1", { expect_identical(init_values, init_values_new) @@ -70,9 +73,11 @@ eval(parse(text = "init_values <- CroptimizR:::get_init_values(param_info)")) param_info_new <- within(param_info, rm(init_values)) eval(parse( text = - "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)")) + "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)" +)) eval(parse( - text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)")) + text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)" +)) test_that("set_init_values: simple case 2", { expect_identical(init_values, init_values_new) @@ -107,10 +112,12 @@ param_info_new <- param_info param_info_new$dlaimax$init_values <- NULL param_info_new$durvieF$init_values <- NULL eval(parse( -text = - "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)")) + text = + "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)" +)) eval(parse( - text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)")) + text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)" +)) test_that("set_init_values: Case 1 with groups of situations per parameter", { expect_identical(init_values, init_values_new) @@ -150,9 +157,11 @@ param_info_new$dlaimax$init_values <- NULL param_info_new$durvieF$init_values <- NULL eval(parse( text = - "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)")) + "param_info_new <- CroptimizR:::set_init_values(param_info, init_values)" +)) eval(parse( - text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)")) + text = "init_values_new <- CroptimizR:::get_init_values(param_info_new)" +)) test_that("set_init_values: Case 2 with groups of situations per parameter", { expect_identical(init_values, init_values_new) diff --git a/tests/testthat/test-intersect_sim_obs.R b/tests/testthat/test-intersect_sim_obs.R index 88fe7a8..b2e9297 100644 --- a/tests/testthat/test-intersect_sim_obs.R +++ b/tests/testthat/test-intersect_sim_obs.R @@ -35,10 +35,12 @@ sim_list <- list(sit1 = data.frame( var1 = rep(2, 4), var2 = c(NA, NA, NA, 4) )) res <- suppressWarnings(eval(parse( - text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)"))) + text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)" +))) test_that("intersect_sim_obs", { expect_warning(eval(parse( - text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)")), "dates") + text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)" + )), "dates") expect_equal(res, NA) }) obs_list <- list(sit1 = data.frame( @@ -53,9 +55,11 @@ sim_list <- list(sit1 = data.frame( var1 = rep(2, 4), var2 = c(NA, NA, NA, 4) )) res <- suppressWarnings(eval(parse( - text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)"))) + text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)" +))) test_that("intersect_sim_obs", { expect_warning(eval(parse( - text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)")), "*variables") + text = "CroptimizR:::intersect_sim_obs(obs_list,sim_list)" + )), "*variables") expect_equal(res, NA) }) diff --git a/tests/testthat/test-isdata.R b/tests/testthat/test-isdata.R index 0cc1c38..4c76def 100644 --- a/tests/testthat/test-isdata.R +++ b/tests/testthat/test-isdata.R @@ -31,9 +31,11 @@ obs_list3 <- list( test_that("is.obs", { expect_true(eval(parse(text = "CroptimizR:::is.obs(obs_list1)"))) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.obs(obs_list2)")))) + eval(parse(text = "CroptimizR:::is.obs(obs_list2)")) + )) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.obs(obs_list3)")))) + eval(parse(text = "CroptimizR:::is.obs(obs_list3)")) + )) }) @@ -74,9 +76,11 @@ sim_list3 <- list( test_that("is.sim", { expect_true(CroptimizR:::is.sim(sim_list1)) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.sim(sim_list2)")))) + eval(parse(text = "CroptimizR:::is.sim(sim_list2)")) + )) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.sim(sim_list3)")))) + eval(parse(text = "CroptimizR:::is.sim(sim_list3)")) + )) }) @@ -127,9 +131,12 @@ data_list4 <- list( test_that("is.data", { expect_true(CroptimizR:::is.data(data_list1)) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.data(data_list2)")))) + eval(parse(text = "CroptimizR:::is.data(data_list2)")) + )) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.data(data_list3)")))) + eval(parse(text = "CroptimizR:::is.data(data_list3)")) + )) expect_false(suppressWarnings( - eval(parse(text = "CroptimizR:::is.data(data_list4)")))) + eval(parse(text = "CroptimizR:::is.data(data_list4)")) + )) }) diff --git a/tests/testthat/test-likelihood.R b/tests/testthat/test-likelihood.R index 208575b..205a639 100644 --- a/tests/testthat/test-likelihood.R +++ b/tests/testthat/test-likelihood.R @@ -13,20 +13,24 @@ sim_list2 <- list( test_that("likelihood_log_ciidn", { expect_equal(likelihood_log_ciidn(sim_list, sim_list), Inf) - expect_equal(likelihood_log_ciidn(obs_list, sim_list), - log((4^ (-4) * 4^ (-2.5)))) + expect_equal( + likelihood_log_ciidn(obs_list, sim_list), + log((4^(-4) * 4^(-2.5))) + ) expect_equal( likelihood_log_ciidn(obs_list2, sim_list2), - log((4^ (-4) * 17^ (-7 / 2))) + log((4^(-4) * 17^(-7 / 2))) ) }) test_that("likelihood_log_ciidn_corr", { expect_equal(likelihood_log_ciidn_corr(sim_list, sim_list), Inf) - expect_equal(likelihood_log_ciidn_corr(obs_list, sim_list), - log((4^ (-5 / 2)))) + expect_equal( + likelihood_log_ciidn_corr(obs_list, sim_list), + log((4^(-5 / 2))) + ) expect_equal( likelihood_log_ciidn_corr(obs_list2, sim_list2), - log(((4 + 13 / 2)^ (-3))) + log(((4 + 13 / 2)^(-3))) ) }) diff --git a/tests/testthat/test-ls_criteria.R b/tests/testthat/test-ls_criteria.R index 16a5162..f7aeddd 100644 --- a/tests/testthat/test-ls_criteria.R +++ b/tests/testthat/test-ls_criteria.R @@ -18,24 +18,50 @@ test_that("crit_ols", { }) test_that("crit_wls", { - expect_equal(crit_wls(sim_list, sim_list, - function(...) { return(Inf) }), 0) - expect_equal(crit_wls(sim_list2, obs_list2, - function(...) { return(1) }), - crit_ols(obs_list2, sim_list2)) - expect_equal(crit_wls(sim_list2, obs_list2, - function(...) { return(Inf) }), 0) - expect_equal(crit_wls(sim_list2, obs_list2, - function(obs, ...) { return(obs) }), 10) - expect_error(crit_wls(sim_list2, obs_list2, - function(obs, ...) { return(0) })) - expect_error(crit_wls(sim_list2, obs_list2, - function(obs, ...) { return(NA) })) + expect_equal(crit_wls( + sim_list, sim_list, + function(...) { + return(Inf) + } + ), 0) + expect_equal( + crit_wls( + sim_list2, obs_list2, + function(...) { + return(1) + } + ), + crit_ols(obs_list2, sim_list2) + ) + expect_equal(crit_wls( + sim_list2, obs_list2, + function(...) { + return(Inf) + } + ), 0) + expect_equal(crit_wls( + sim_list2, obs_list2, + function(obs, ...) { + return(obs) + } + ), 10) + expect_error(crit_wls( + sim_list2, obs_list2, + function(obs, ...) { + return(0) + } + )) + expect_error(crit_wls( + sim_list2, obs_list2, + function(obs, ...) { + return(NA) + } + )) }) test_that("crit_log_cwss", { expect_equal(crit_log_cwss(sim_list, sim_list), -Inf) expect_equal(crit_log_cwss(sim_list, obs_list), log(2)) - expect_equal(crit_log_cwss(sim_list2, obs_list2), log((17 / 3)^ (3 / 2))) + expect_equal(crit_log_cwss(sim_list2, obs_list2), log((17 / 3)^(3 / 2))) }) test_that("crit_log_cwss_corr", { expect_equal(crit_log_cwss_corr(sim_list, sim_list), -Inf) diff --git a/tests/testthat/test-obsSim_consistency.R b/tests/testthat/test-obsSim_consistency.R index 4b6c232..147fb52 100644 --- a/tests/testthat/test-obsSim_consistency.R +++ b/tests/testthat/test-obsSim_consistency.R @@ -15,7 +15,8 @@ obs_list <- list( # Check if nothing is captured if sim and obs lists are identical res <- eval(parse( - text = "CroptimizR:::make_obsSim_consistent(obs_list, obs_list)")) + text = "CroptimizR:::make_obsSim_consistent(obs_list, obs_list)" +)) test_that("obs and sim are identical", { expect_identical(res$obs_list, obs_list) expect_identical(res$sim_list, obs_list) @@ -28,10 +29,13 @@ sim_list <- obs_list sim_list[[1]]$var1 <- as.character(sim_list[[1]]$var1) sim_list[[3]]$var2 <- as.character(sim_list[[3]]$var2) test_that("Capture non-consistent variable types", { - expect_error(eval(parse( - text = - "CroptimizR:::check_obsSim_consistency(sim_list, obs_list)")), - "different types") + expect_error( + eval(parse( + text = + "CroptimizR:::check_obsSim_consistency(sim_list, obs_list)" + )), + "different types" + ) }) # Check if it captures that some Date columns are of non expected type @@ -39,9 +43,12 @@ sim_list <- obs_list sim_list[[1]]$Date <- as.character(sim_list[[1]]$Date) sim_list[[3]]$Date <- as.character(sim_list[[3]]$Date) test_that("Capture unexpected type for Date columns", { - expect_error(eval( - parse(text = "CroptimizR:::check_obsSim_consistency(sim_list, obs_list)")), - "incorrect format") + expect_error( + eval( + parse(text = "CroptimizR:::check_obsSim_consistency(sim_list, obs_list)") + ), + "incorrect format" + ) }) # Check if it corrects non consistent types for sim and obs Dates @@ -49,7 +56,8 @@ sim_list <- obs_list sim_list[[1]]$Date <- as.Date(sim_list[[1]]$Date) sim_list[[3]]$Date <- as.Date(sim_list[[3]]$Date) res <- eval(parse( - text = "CroptimizR:::make_obsSim_consistent(sim_list, obs_list)")) + text = "CroptimizR:::make_obsSim_consistent(sim_list, obs_list)" +)) test_that("Check correction of dates types", { expect_true(lubridate::is.POSIXct(res$sim_list[[1]]$Date)) expect_true(lubridate::is.POSIXct(res$sim_list[[3]]$Date)) @@ -58,9 +66,10 @@ test_that("Check correction of dates types", { # Check if is_sim_inf_or_na return FALSE when it must sim_list <- obs_list -param_values <- c(p1=1.0, p2=2.0) +param_values <- c(p1 = 1.0, p2 = 2.0) res <- eval(parse( - text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)")) + text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)" +)) test_that("Check is_sim_inf_or_na return FALSE when sim is not Inf neither NA when there is a corresponding observed value", { expect_false(res) }) @@ -68,21 +77,27 @@ test_that("Check is_sim_inf_or_na return FALSE when sim is not Inf neither NA wh # Check if is_sim_inf_or_na return TRUE when it must, missing value for one variable, one situation and one date sim_list <- obs_list sim_list$sit1$var1[[1]] <- Inf -param_values <- c(p1=1.0, p2=2.0) +param_values <- c(p1 = 1.0, p2 = 2.0) test_that("Check is_sim_inf_or_na return TRUE when sim is Inf or NA when there is a corresponding observed value, case 1", { - expect_warning(eval(parse( - text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)")), - "sit1.*var1.*2009-11-30") + expect_warning( + eval(parse( + text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)" + )), + "sit1.*var1.*2009-11-30" + ) }) # Check if is_sim_inf_or_na return TRUE when it must, missing values for several dates sim_list <- obs_list sim_list$sit3$var2 <- NA -param_values <- c(p1=1.0, p2=2.0) +param_values <- c(p1 = 1.0, p2 = 2.0) test_that("Check is_sim_inf_or_na return TRUE when sim is Inf or NA when there is a corresponding observed value, case 1", { - expect_warning(eval(parse( - text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)")), - "sit3.*var2.*2010-10-03.*2010-10-04") + expect_warning( + eval(parse( + text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)" + )), + "sit3.*var2.*2010-10-03.*2010-10-04" + ) }) @@ -91,10 +106,12 @@ sim_list <- obs_list sim_list$sit1$var1[[1]] <- Inf sim_list$sit3$var1[[2]] <- NA sim_list$sit3$var2 <- NA -param_values <- c(p1=1.0, p2=2.0) +param_values <- c(p1 = 1.0, p2 = 2.0) test_that("Check is_sim_inf_or_na return TRUE when sim is Inf or NA when there is a corresponding observed value, case 1", { - expect_warning(eval(parse( - text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)")), - "sit1.*var1.*2009-11-30.*sit3.*var1.*2010-10-04.*var2.*2010-10-03.*2010-10-04") + expect_warning( + eval(parse( + text = "CroptimizR:::is_sim_inf_or_na(sim_list, obs_list, param_values)" + )), + "sit1.*var1.*2009-11-30.*sit3.*var1.*2010-10-04.*var2.*2010-10-03.*2010-10-04" + ) }) - diff --git a/tests/testthat/test-optim_switch.R b/tests/testthat/test-optim_switch.R index 835e930..6638cad 100644 --- a/tests/testthat/test-optim_switch.R +++ b/tests/testthat/test-optim_switch.R @@ -1,9 +1,11 @@ test_that("Returns an error when unknown method is given", { expect_error( - eval(parse(text = - "CroptimizR:::optim_switch(optim_method=\"unknown_method\", + eval(parse( + text = + "CroptimizR:::optim_switch(optim_method=\"unknown_method\", optim_options=NULL, param_info=NULL, - crit_options=NULL)")), + crit_options=NULL)" + )), "Unknown method" ) }) diff --git a/tests/testthat/test-sample_params.R b/tests/testthat/test-sample_params.R index 4b57999..d0604c4 100644 --- a/tests/testthat/test-sample_params.R +++ b/tests/testthat/test-sample_params.R @@ -1,10 +1,14 @@ bounds_list <- list(lb = c(0, 1, 2), ub = c(1, 2, 3)) -res1 <- eval(parse(text = - "CroptimizR:::sample_params(bounds_list, 5, seed = 1234)")) +res1 <- eval(parse( + text = + "CroptimizR:::sample_params(bounds_list, 5, seed = 1234)" +)) Sys.sleep(2) res2 <- eval(parse(text = "CroptimizR:::sample_params(bounds_list, 5)")) -res3 <- eval(parse(text = - "CroptimizR:::sample_params(bounds_list, 5, seed = 1234)")) +res3 <- eval(parse( + text = + "CroptimizR:::sample_params(bounds_list, 5, seed = 1234)" +)) test_that("sample_params", { expect_false(isTRUE(all.equal(res1, res2))) diff --git a/vignettes/AgMIP_Calibration_Phenology_protocol.Rmd b/vignettes/AgMIP_Calibration_Phenology_protocol.Rmd index 5dc511c..531c771 100644 --- a/vignettes/AgMIP_Calibration_Phenology_protocol.Rmd +++ b/vignettes/AgMIP_Calibration_Phenology_protocol.Rmd @@ -65,14 +65,15 @@ The STICS crop model is used in this example. Initialization steps required by t ## Plotting the observations ```{r setup_initializations, eval=params$eval_auto_test, message=FALSE, results=FALSE, warning=FALSE, include=FALSE} - # DEFINE THE PATH TO YOUR LOCALLY INSTALLED VERSION OF JAVASTICS javastics_path <- path_to_JavaStics # Download the example USMs and define the path to the JavaStics workspace # (JavaStics XML input files): -data_dir <- file.path(SticsRFiles::download_data(example_dirs = "study_case_1", - stics_version = "V9.0")) +data_dir <- file.path(SticsRFiles::download_data( + example_dirs = "study_case_1", + stics_version = "V9.0" +)) javastics_workspace_path <- file.path(data_dir, "XmlFiles") @@ -83,22 +84,29 @@ model_wrapper <- stics_wrapper stics_inputs_path <- file.path(data_dir, "TxtFiles") dir.create(stics_inputs_path, showWarnings = FALSE) -gen_usms_xml2txt(javastics = javastics_path, - workspace = javastics_workspace_path, - out_dir = stics_inputs_path, verbose = TRUE) +gen_usms_xml2txt( + javastics = javastics_path, + workspace = javastics_workspace_path, + out_dir = stics_inputs_path, verbose = TRUE +) ``` ```{r eval=params$eval_auto_test, results='hide', message=FALSE, warning=FALSE, include=FALSE} - # Set the model options (see '? stics_wrapper_options' for details) -model_options <- stics_wrapper_options(javastics = javastics_path, - workspace = stics_inputs_path, - parallel = TRUE, cores = 4) +model_options <- stics_wrapper_options( + javastics = javastics_path, + workspace = stics_inputs_path, + parallel = TRUE, cores = 4 +) # Run the model on all situations found in stics_inputs_path -sim_res <- stics_wrapper(model_options = model_options, - param_values = c(stlevamf = 350, stamflax = 650, - tdmin = 7, tdmax = 26), - var = c("iamfs", "ilaxs")) +sim_res <- stics_wrapper( + model_options = model_options, + param_values = c( + stlevamf = 350, stamflax = 650, + tdmin = 7, tdmax = 26 + ), + var = c("iamfs", "ilaxs") +) set.seed(1234) obs_list <- lapply(sim_res$sim_list, function(x) { tmp <- x[nrow(x), c("Date", "iamfs", "ilaxs")] @@ -131,8 +139,10 @@ Let's plot the values of these observations: ```{r eval=params$eval_auto_vignette, message=FALSE, warning=FALSE} obs_df <- dplyr::bind_rows(obs_list) %>% tidyr::pivot_longer( - cols = c("iamfs", "ilaxs"), names_to = "stage") -ggplot(obs_df, aes(x = stage, y = value)) + geom_boxplot() + cols = c("iamfs", "ilaxs"), names_to = "stage" +) +ggplot(obs_df, aes(x = stage, y = value)) + + geom_boxplot() ``` The "end juvenile" stage occurred around day 180 after sowing, while the "maximum LAI" stage occurred around day 230 after sowing. @@ -148,7 +158,7 @@ In this example, we will consider 5 parameters of the STICS crop model: `stlevam param_info <- list( lb = c(stlevamf = 100, stamflax = 300, tdmin = 4, stressdev = 0, tdmax = 25), ub = c(stlevamf = 500, stamflax = 800, tdmin = 8, stressdev = 1, tdmax = 32) - ) +) ``` Note that `-Inf` or `Ìnf` can be used for, respectively, lower and upper bounds. In that case, initial values must @@ -211,20 +221,21 @@ By default, the information criterion used for the parameter selection is the BI Note, however, that other criteria can be used (AIC and AICc) using the `info_crit_func` argument. ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} -res <- estim_param(obs_list = obs_list, - crit_function = crit_ols, - model_function = model_wrapper, - model_options = model_options, - optim_options = optim_options, - forced_param_values = forced_param_values, - candidate_param = candidate_param, - param_info = param_info) +res <- estim_param( + obs_list = obs_list, + crit_function = crit_ols, + model_function = model_wrapper, + model_options = model_options, + optim_options = optim_options, + forced_param_values = forced_param_values, + candidate_param = candidate_param, + param_info = param_info +) ``` At the end of the execution, some information such as the selected step number, the list of selected parameters, their estimated values and some stats about the execution time are printed in the R console: ```{r eval=FALSE, echo=TRUE} - ... # ---------------------- # End of parameter selection process @@ -273,9 +284,11 @@ load(file.path("ResultsAgmipPhenology", "optim_results.Rdata")) ```{r eval=params$eval_manual_vignette, message=FALSE, warning=FALSE} param_values <- c(res$final_values, res$forced_param_values) -sim_after_optim <- model_wrapper(param_values = param_values, - model_options = model_options, - var = c("iamfs", "ilaxs")) +sim_after_optim <- model_wrapper( + param_values = param_values, + model_options = model_options, + var = c("iamfs", "ilaxs") +) ``` Then, the simulated versus observed values can be plotted using the CroPlotR plot function: @@ -296,13 +309,17 @@ knitr::include_graphics("ResultsAgmipPhenology/all_situations.png") MSE and its components can be calcuted for each variable using the CroPlotR summary function: ```{r eval=FALSE, warning=FALSE, message=FALSE} -summary(sim_after_optim$sim_list, obs = obs_list, - stats = c("MSE", "Bias2", "SDSD", "LCS")) +summary(sim_after_optim$sim_list, + obs = obs_list, + stats = c("MSE", "Bias2", "SDSD", "LCS") +) ``` ```{r eval=params$eval_manual_vignette, warning=FALSE, message=FALSE, echo=FALSE} -stats <- summary(sim_after_optim$sim_list, obs = obs_list, - stats = c("MSE", "Bias2", "SDSD", "LCS")) +stats <- summary(sim_after_optim$sim_list, + obs = obs_list, + stats = c("MSE", "Bias2", "SDSD", "LCS") +) save(stats, file = file.path(params$result_path, "stats.Rdata")) ``` @@ -315,16 +332,26 @@ print(stats) ```{r move_results, eval=params$eval_manual_vignette, include=FALSE} # Move the files produced since the temp. folder is removed after Rmd execution file.copy(file.path(optim_options$out_dir, "param_selection_steps.csv"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) file.copy(file.path(optim_options$out_dir, "param_selection_steps.Rdata"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) file.copy(file.path(optim_options$out_dir, "optim_results.Rdata"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) file.copy(file.path(optim_options$out_dir, "all_situations.png"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) dir.create(file.path(params$result_path, "results_all_steps")) -file.copy(from = file.path(optim_options$out_dir, "results_all_steps"), - to = file.path(params$result_path, "results_all_steps"), - recursive = TRUE, - overwrite = TRUE) +file.copy( + from = file.path(optim_options$out_dir, "results_all_steps"), + to = file.path(params$result_path, "results_all_steps"), + recursive = TRUE, + overwrite = TRUE +) ``` diff --git a/vignettes/ApsimX_parameter_estimation_simple_case.Rmd b/vignettes/ApsimX_parameter_estimation_simple_case.Rmd index 7e40da9..f01ebd9 100644 --- a/vignettes/ApsimX_parameter_estimation_simple_case.Rmd +++ b/vignettes/ApsimX_parameter_estimation_simple_case.Rmd @@ -35,7 +35,6 @@ The parameter estimation is performed using the Nelder-Mead simplex method imple ```{r setup_initializations, message=FALSE, results=FALSE, warning=FALSE} - # Install and load the needed libraries if (!require("CroptimizR")) { devtools::install_github("SticsRPacks/CroptimizR@*release") @@ -79,11 +78,10 @@ apsimx_path <- "D:\\Home\\sbuis\\Documents\\OUTILS-INFORMATIQUE\\APSIM2021.01.20 ```{r message=FALSE, warning=FALSE} - sit_name <- "GattonRowSpacingRowSpace25cm" - # among "GattonRowSpacingRowSpace25cm", - # "GattonRowSpacingRowSpace50cm", - # "GattonRowSpacingRowSpaceN0" +# among "GattonRowSpacingRowSpace25cm", +# "GattonRowSpacingRowSpace50cm", +# "GattonRowSpacingRowSpaceN0" var_name <- c("Wheat.Leaf.LAI") # or "Wheat.AboveGround.Wt" ``` @@ -96,11 +94,10 @@ In this case, the argument `param_values` of the wrapper is not set: the values ```{r results='hide', message=FALSE, warning=FALSE} - - # Set the model options (see '? apsimx_wrapper_options' for details) files_path <- system.file(file.path("extdata", "apsimx_files"), - package = "ApsimOnR") + package = "ApsimOnR" +) apsimx_file <- file.path(files_path, "template.apsimx") # Setting met files path @@ -113,13 +110,15 @@ obs_files_path <- files_path predicted_table_name <- "DailyReport" observed_table_name <- "Observed" -model_options <- apsimx_wrapper_options(apsimx_path = apsimx_path, - apsimx_file = apsimx_file, - variable_names = var_name, - predicted_table_name = predicted_table_name, - met_files_path = met_files_path, - observed_table_name = observed_table_name, - obs_files_path = obs_files_path) +model_options <- apsimx_wrapper_options( + apsimx_path = apsimx_path, + apsimx_file = apsimx_file, + variable_names = var_name, + predicted_table_name = predicted_table_name, + met_files_path = met_files_path, + observed_table_name = observed_table_name, + obs_files_path = obs_files_path +) # Run the model (on all situations found in the apsimx_file) sim_before_optim <- apsimx_wrapper(model_options = model_options) @@ -136,10 +135,12 @@ We only keep observations for situation `sit_name` and variable `var_name` (`obs # simulation ran before optimization. But they may be loaded using the original # xlsx data file (from the files_path) -obs_list <- read_apsimx_output(sim_before_optim$db_file_name, - model_options$observed_table_name, - model_options$variable_names, - names(sim_before_optim$sim_list)) +obs_list <- read_apsimx_output( + sim_before_optim$db_file_name, + model_options$observed_table_name, + model_options$variable_names, + names(sim_before_optim$sim_list) +) obs_list <- filter_obs(obs_list, situation = sit_name, include = TRUE) ``` @@ -156,10 +157,16 @@ Initial values for the minimization can also be provided in `param_info` (see `? ```{r message=FALSE, warning=FALSE} # 2 parameters here: ExtinctionCoeff and RUE, of bounds [0.4,0.6] and [1.4,1.6] param_info <- - list(lb = c(.Simulations.Replacements.Wheat.Leaf.ExtinctionCoeff.VegetativePhase.FixedValue = 0.4, - .Simulations.Replacements.Wheat.Leaf.Photosynthesis.RUE.FixedValue = 1.4), - ub = c(.Simulations.Replacements.Wheat.Leaf.ExtinctionCoeff.VegetativePhase.FixedValue = 0.6, - .Simulations.Replacements.Wheat.Leaf.Photosynthesis.RUE.FixedValue = 1.6)) + list( + lb = c( + .Simulations.Replacements.Wheat.Leaf.ExtinctionCoeff.VegetativePhase.FixedValue = 0.4, + .Simulations.Replacements.Wheat.Leaf.Photosynthesis.RUE.FixedValue = 1.4 + ), + ub = c( + .Simulations.Replacements.Wheat.Leaf.ExtinctionCoeff.VegetativePhase.FixedValue = 0.6, + .Simulations.Replacements.Wheat.Leaf.Photosynthesis.RUE.FixedValue = 1.6 + ) + ) ``` @@ -174,24 +181,23 @@ The number of repetitions `nb_rep` is advised to be set at least to 5, while 10 ```{r message=FALSE, warning=FALSE} - optim_options <- list() optim_options$nb_rep <- 7 # Number of repetitions of the minimization - # (each time starting with different initial - # values for the estimated parameters) +# (each time starting with different initial +# values for the estimated parameters) optim_options$maxeval <- 500 # Maximum number of evaluations of the - # minimized criteria +# minimized criteria optim_options$xtol_rel <- 1e-03 # Tolerance criterion between two iterations - # (threshold for the relative difference of - # parameter values between the 2 previous - # iterations) +# (threshold for the relative difference of +# parameter values between the 2 previous +# iterations) optim_options$out_dir <- getwd() # path where to store the results - # (graph and Rdata) +# (graph and Rdata) optim_options$ranseed <- 1234 # set random seed so that each execution give the - # same results. If you want randomization, - # don't set it. +# same results. If you want randomization, +# don't set it. ``` @@ -203,12 +209,13 @@ Same for crit_function: a value is set by default (`crit_log_cwss`, see `? crit_ ```{r message=FALSE, warning=FALSE} - -res <- estim_param(obs_list = obs_list, - model_function = apsimx_wrapper, - model_options = model_options, - optim_options = optim_options, - param_info = param_info) +res <- estim_param( + obs_list = obs_list, + model_function = apsimx_wrapper, + model_options = model_options, + optim_options = optim_options, + param_info = param_info +) ``` The estimated values of the parameters are the following: @@ -227,7 +234,6 @@ Complementary graphs and data are stored in the folder which path is given in `o ```{r eval=TRUE, echo=FALSE, out.width = '45%'} - knitr::include_graphics("ResultsSimpleCaseApsim/estimInit_ExtinctionCoeff.png") knitr::include_graphics("ResultsSimpleCaseApsim/estimInit_RUE.png") @@ -242,9 +248,10 @@ In this case, the `param_values` argument is set so that estimated values of the ```{r results='hide', message=FALSE, warning=FALSE} - -sim_after_optim <- apsimx_wrapper(param_values = res$final_values, - model_options = model_options) +sim_after_optim <- apsimx_wrapper( + param_values = res$final_values, + model_options = model_options +) ``` @@ -262,8 +269,10 @@ p2 <- p[[sit_name]] + labs(title = "After Optimization") + p <- grid.arrange(grobs = list(p1, p2), nrow = 1, ncol = 2) # Save the graph -ggsave(file.path(optim_options$out_dir, - paste0("sim_obs_plots", ".png")), plot = p) +ggsave(file.path( + optim_options$out_dir, + paste0("sim_obs_plots", ".png") +), plot = p) ``` @@ -278,9 +287,15 @@ knitr::include_graphics("ResultsSimpleCaseApsim/sim_obs_plots.png") ```{r move_results, include=FALSE} # Move the files produced since the temp. folder is removed after Rmd execution file.copy(file.path(optim_options$out_dir, "EstimatedVSinit.pdf"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) file.copy(file.path(optim_options$out_dir, "optim_results.Rdata"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) file.copy(file.path(optim_options$out_dir, "sim_obs_plots.png"), - params$result_path, overwrite = TRUE) + params$result_path, + overwrite = TRUE +) ``` diff --git a/vignettes/CroptimizR.Rmd b/vignettes/CroptimizR.Rmd index 2709b4b..07b0035 100644 --- a/vignettes/CroptimizR.Rmd +++ b/vignettes/CroptimizR.Rmd @@ -90,7 +90,7 @@ Basically, to use CroptimizR on their own case, users must provide: The parameter estimation process automatically runs the given crop model on a set of observed situations by forcing the estimated parameters to the values proposed by the selected algorithm. It will then compute the selected goodness-of-fit criterion using the simulated values and corresponding observations of the model output variables (see Figure 1). ```{r, out.width='85%', echo=FALSE, fig.align='center', fig.cap='Figure 1: Schematic representation of the parameter estimation process.'} -knitr::include_graphics('GettingStarted/main_scheme.png') +knitr::include_graphics("GettingStarted/main_scheme.png") ``` In output, CroptimizR generates specific statistics and graphs depending on the @@ -115,7 +115,7 @@ Observed values of model output variables have to be provided under the shape of Let's consider the observations are, as often, available under the shape of a text file as displayed in Figure 2. ```{r, out.width='50%', echo=FALSE, fig.align='center', fig.cap='Figure 2: Example of a text file defining observations.'} -knitr::include_graphics('GettingStarted/obs_file.png') +knitr::include_graphics("GettingStarted/obs_file.png") ``` In this example, observations are available for two situations, called here `sit1` and `sit2`, three variables for situation `sit1` (`LAI`, `Biomass` and `Yield`) and two for `sit2` (`Biomass` and `Yield`). Observations of `Yield` where only provided for the last date. diff --git a/vignettes/Designing_a_model_wrapper.Rmd b/vignettes/Designing_a_model_wrapper.Rmd index 388a135..2d1cca9 100644 --- a/vignettes/Designing_a_model_wrapper.Rmd +++ b/vignettes/Designing_a_model_wrapper.Rmd @@ -64,14 +64,13 @@ Here's a header for a basic version of a model wrapper: #' @return A list containing: #' o `sim_list`: a `named list` (names = situations names) of data.frames #' (or tibbles) of simulated output values (one element of the list per -#` simulated situation) +# ` simulated situation) #' o `error`: an error code indicating if at least one simulation ended #' with an error. model_wrapper <- function(param_values = NULL, situation, model_options, ...) { } - ``` **Each argument detailed here must be defined for any CroptimizR model wrapper. You can use this header to develop yours. @@ -111,13 +110,12 @@ The data.frames must have one column called `Date` containing the simulations da For example, if `situation` is `c("situation1","situation2","situation3"), sim_list should look like: ```{r eval=FALSE} - sim_list$situation1 #> A tibble: *** x *** #> Date var1 var2 var3 var4 ... #> -#> 1 1994-10-17 0 2.53 4.80 150 +#> 1 1994-10-17 0 2.53 4.80 150 #> 2 1994-10-18 0 2.31 4.66 150 #> 3 1994-10-19 0 4.55 4.44 150 #> @@ -141,8 +139,6 @@ sim_list$situation3 #> 1 1996-10-17 0 2.41 4.81 142 #> 2 1996-10-18 0 3.03 4.71 142 #> 3 1996-10-19 0 5.10 4.47 142 - - ``` ** Note that there is not for the moment any dedicated data structure handled in CroptimizR and CroPlotR for storing simulated variables that do no depend on time (e.g. such as the julian day of a given phenological stage). You must thus include them in the data.frames described here-before. For this variables, it is strongly advised to replicate their (unique) values for all rows of the data.frames (see e.g. var4 in the data.frames described here-before). This makes it possible to compare these simulated values to the corresponding observed values regardless of the dates with which they are associated in the data.frames including the observations.** @@ -159,13 +155,15 @@ If any simulation failed for any reason, use the R function `warning` to print a A typical pseudo-code implementation of a basic wrapper function is thus: ```{r eval=FALSE} - model_wrapper <- function(param_values = NULL, situation, model_options, ...) { # Initializations # param_names <- names(param_values) - results <- list(sim_list = setNames(vector("list", length(situation)), - nm = situation), - error = FALSE) + results <- list( + sim_list = setNames(vector("list", length(situation)), + nm = situation + ), + error = FALSE + ) attr(results$sim_list, "class") <- "cropr_simulation" for (situation in situation) { @@ -181,7 +179,6 @@ model_wrapper <- function(param_values = NULL, situation, model_options, ...) { } return(results) } - ``` ### Test your wrapper @@ -214,15 +211,18 @@ If all tests succeed, it returns: ```{r eval=TRUE, echo=FALSE} cat(crayon::green( - "Test the wrapper returns outputs in expected format ...\n")) + "Test the wrapper returns outputs in expected format ...\n" +)) cat(crayon::green("... OK\n")) cat("\n") cat(crayon::green( -"Test the wrapper gives identical results when running with same inputs ...\n")) + "Test the wrapper gives identical results when running with same inputs ...\n" +)) cat(crayon::green("... OK\n")) cat("\n") cat(crayon::green( -"Test the wrapper gives different results when running with different inputs ...\n")) + "Test the wrapper gives different results when running with different inputs ...\n" +)) cat(crayon::green("... OK\n")) ``` @@ -251,17 +251,15 @@ If you want your model wrapper to be used for simultaneous estimation of specifi Depending on the application (simple case or simultaneous specific and varietal parameters estimation), the `estim_param` function will thus pass to the model wrapper either a named vector or a data.frame, including or not the `situation` column. Both shapes have thus to be handled in the model wrapper in this case. This may be done as in the following piece of code: ```{r, eval=FALSE} - # convert param_values in a tibble param_values <- tibble::tibble(!!!param_values) # loop on the situations to simulate for (sit in situation_list) { - # extract the parameters values to apply for the given situation params <- NULL if (!is.null(param_values)) { - if (! "situation" %in% names(param_values)) { + if (!"situation" %in% names(param_values)) { params <- param_values } else { params <- dplyr::filter(param_values, situation == sit) %>% @@ -270,7 +268,6 @@ for (sit in situation_list) { } # run the model with parameters values defined by params - } ``` @@ -293,7 +290,7 @@ test_parallel <- function(cores_nb = 1, pa = FALSE, max_it = 5) { # pa: flag to switch from case with pre-allocated list to case with - #returned statement + # returned statement # Launching the cluster cl <- parallel::makeCluster(cores_nb) @@ -305,9 +302,9 @@ test_parallel <- function(cores_nb = 1, # Parallel loop out <- foreach(i = 1:max_it) %dopar% { if (pa) { - out_pa[[i]] <- i # store results in a pre-allocted list + out_pa[[i]] <- i # store results in a pre-allocted list } else { - return(i) # return the results + return(i) # return the results } } @@ -332,7 +329,6 @@ out_pa <- test_parallel(cores_nb, TRUE) out out_pa - ``` * Outputs selection @@ -351,7 +347,6 @@ It is advised to define them as optionnal arguments with default value equal to ### Example of a basic wrapper for an LAI toy model ```{r, eval=TRUE} - lai_toymodel <- function(year, max_lai = 8, julday_incslope = 100, inc_slope = 5, julday_decslope = 200, dec_slope = 2) { # Simulate lai for a single crop over 2 years @@ -368,12 +363,14 @@ lai_toymodel <- function(year, max_lai = 8, julday_incslope = 100, # - lai: vector of simulated lai # - dates: vector of dates (POSIXct) for which lai is computed - end_day <- format(as.Date(paste0(year + 1, "-12-31"), format = "%Y-%m-%d", - origin = paste0(year, "-01-01")), "%j") + end_day <- format(as.Date(paste0(year + 1, "-12-31"), + format = "%Y-%m-%d", + origin = paste0(year, "-01-01") + ), "%j") jul_days <- 1:as.numeric(end_day) lai <- max_lai * (1 / (1 + exp((julday_incslope - jul_days) / inc_slope)) - - 1 / (1 + exp((julday_decslope - jul_days) / dec_slope))) + 1 / (1 + exp((julday_decslope - jul_days) / dec_slope))) lai[lai < 0] <- 0 dates <- @@ -381,7 +378,8 @@ lai_toymodel <- function(year, max_lai = 8, julday_incslope = 100, as.character( as.Date( jul_days, - origin = paste0(year, "-01-01")) + origin = paste0(year, "-01-01") + ) ), format = "%Y-%m-%d", tz = "UTC" ) @@ -392,7 +390,6 @@ lai_toymodel <- function(year, max_lai = 8, julday_incslope = 100, laitm_simple_wrapper <- function(param_values = NULL, situation, model_options, ...) { - # A basic wrapper for lai_toymodel # # Arguments @@ -430,7 +427,6 @@ laitm_simple_wrapper <- function(param_values = NULL, situation, attr(results$sim_list, "class") <- "cropr_simulation" for (sit in situation) { - # Retrieve year, emergence and crop_duration from situation name tmp <- stringr::str_split(sit, "_") year <- as.numeric(tmp[[1]][[1]]) @@ -453,8 +449,7 @@ laitm_simple_wrapper <- function(param_values = NULL, situation, "dec_slope", "julday_incslope", "julday_decslope" - ) - )) { + ))) { warning( paste( "Unknown parameters in param_values:", @@ -535,10 +530,11 @@ if (!require("CroPlotR")) { devtools::install_github("SticsRPacks/CroPlotR@*release") library("CroPlotR") } -tmp <- laitm_simple_wrapper(situation = c("2005_a", "2006_b"), - param_values = res$final_values) +tmp <- laitm_simple_wrapper( + situation = c("2005_a", "2006_b"), + param_values = res$final_values +) plot(tmp$sim_list, obs = obs_synth, type = "scatter") - ``` ### Another version taking into account the different shapes param_values can take @@ -548,7 +544,6 @@ This one can be used to perform simultaneous estimation of specific and varietal ```{r, eval=TRUE} laitm_simple_wrapper_v2 <- function(param_values = NULL, situation, model_options, ...) { - # A basic wrapper for lai_toymodel # # Arguments @@ -586,7 +581,7 @@ laitm_simple_wrapper_v2 <- function(param_values = NULL, situation, ) attr(results$sim_list, "class") <- "cropr_simulation" - param_values <- tibble::tibble(!!!param_values) # convert param_values + param_values <- tibble::tibble(!!!param_values) # convert param_values # in a tibble for (sit in situation) { @@ -609,7 +604,7 @@ laitm_simple_wrapper_v2 <- function(param_values = NULL, situation, # extract the parameters values to apply for the given situation params <- NULL if (!is.null(param_values)) { - if (! "situation" %in% names(param_values)) { + if (!"situation" %in% names(param_values)) { params <- param_values } else { params <- @@ -623,8 +618,7 @@ laitm_simple_wrapper_v2 <- function(param_values = NULL, situation, "dec_slope", "julday_incslope", "julday_decslope" - ) - )) { + ))) { warning( paste( "Unknown parameters in param_values:", @@ -710,10 +704,12 @@ param_info <- list( ) optim_options <- list(nb_rep = 5, maxeval = 100, xtol_rel = 1e-2) -res <- estim_param(obs_synth, crit_function = crit_ols, - model_function = laitm_simple_wrapper_v2, - optim_options = optim_options, - param_info = param_info) +res <- estim_param(obs_synth, + crit_function = crit_ols, + model_function = laitm_simple_wrapper_v2, + optim_options = optim_options, + param_info = param_info +) res$final_values # Plot the simulations obtained with the optimized values of the parameters VS diff --git a/vignettes/Parameter_estimation_DREAM.Rmd b/vignettes/Parameter_estimation_DREAM.Rmd index 45e7c6f..d40bff7 100644 --- a/vignettes/Parameter_estimation_DREAM.Rmd +++ b/vignettes/Parameter_estimation_DREAM.Rmd @@ -48,7 +48,6 @@ library(CroptimizR) ```{r setup_initializations, eval=params$eval_auto_test, message=FALSE, results=FALSE, warning=FALSE, echo=FALSE} - # DEFINE THE PATH TO YOUR LOCALLY INSTALLED VERSION OF JAVASTICS javastics_path <- path_to_JavaStics @@ -141,12 +140,12 @@ The main basic options are set here. The reader can refer to `? DREAMzs` (see do ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} optim_options <- list() optim_options$iterations <- 10000 # Total number of iterations - # (=> optim_options$iterations/ - # optim_options$startValue - # iterations per chain) +# (=> optim_options$iterations/ +# optim_options$startValue +# iterations per chain) optim_options$startValue <- 3 # Number of markov chains optim_options$out_dir <- data_dir # path where to store the results - # (graph and Rdata) +# (graph and Rdata) optim_options$ranseed <- 1234 # seed for random numbers ``` @@ -234,11 +233,14 @@ In this case, we obtain: ```{r eval=params$eval_auto_vignette, echo=FALSE, out.width = '45%'} knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/Gelman_dlaimax.PNG") + "ResultsSpecificVarietalDREAM/N10000/Gelman_dlaimax.PNG" +) knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/Gelman_durvieF1.PNG") + "ResultsSpecificVarietalDREAM/N10000/Gelman_durvieF1.PNG" +) knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/Gelman_durvieF2.PNG") + "ResultsSpecificVarietalDREAM/N10000/Gelman_durvieF2.PNG" +) ``` Figure 1: Gelman diagnostic. @@ -249,11 +251,14 @@ marginalPlots.pdf shows estimation of the prior (blue) and posterior (pink) dens ```{r eval=params$eval_auto_vignette, echo=FALSE, out.width = '45%'} knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/densities_dlaimax.PNG") + "ResultsSpecificVarietalDREAM/N10000/densities_dlaimax.PNG" +) knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/densities_durvieF1.PNG") + "ResultsSpecificVarietalDREAM/N10000/densities_durvieF1.PNG" +) knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/densities_durvieF2.PNG") + "ResultsSpecificVarietalDREAM/N10000/densities_durvieF2.PNG" +) ``` Figure 2: Prior and posterior densities. @@ -263,7 +268,8 @@ correlationPlots.pdf gives information about the correlation between parameters: ```{r eval=params$eval_auto_vignette, echo=FALSE, out.width = '75%'} knitr::include_graphics( - "ResultsSpecificVarietalDREAM/N10000/correlation_plot.PNG") + "ResultsSpecificVarietalDREAM/N10000/correlation_plot.PNG" +) ``` Figure 3: Correlation plot. diff --git a/vignettes/Parameter_estimation_Specific_and_Varietal.Rmd b/vignettes/Parameter_estimation_Specific_and_Varietal.Rmd index e3e8c0a..198f7ea 100644 --- a/vignettes/Parameter_estimation_Specific_and_Varietal.Rmd +++ b/vignettes/Parameter_estimation_Specific_and_Varietal.Rmd @@ -53,7 +53,6 @@ library(gridExtra) ```{r setup_initializations, eval=params$eval_auto_test, message=FALSE, results=FALSE, warning=FALSE, echo=FALSE} - # DEFINE THE PATH TO YOUR LOCALLY INSTALLED VERSION OF JAVASTICS javastics_path <- path_to_JavaStics @@ -138,7 +137,6 @@ param_info$durvieF <- list( ## Set options for the model ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} - model_options <- stics_wrapper_options( javastics = javastics_path, @@ -151,7 +149,6 @@ model_options <- ## Set options for the parameter estimation method ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} - optim_options <- list() optim_options$nb_rep <- 7 # Number of repetitions of the minimization # (each time starting with different initial @@ -212,7 +209,6 @@ Figure 1: plots of estimated vs initial values of parameters dlaimax and durvieF A piece of code to run the model with the values of the parameters before and after the optimization and to create a couple of plots using the [CroPlotR](https://github.com/SticsRPacks/CroPlotR) package to check if the calibration reduced the difference between simulations and observations: ```{r eval=params$eval_manual_vignette, message=FALSE, warning=FALSE} - # Run the model without and with forcing the optimized values of the parameters sim_before_optim <- stics_wrapper(model_options = model_options) sim_after_optim <- stics_wrapper( @@ -230,11 +226,12 @@ p <- plot( p1 <- p[[1]] + - labs(title = - paste( - "Before Optimization, \n RMSE=", - round(stats$RMSE, digits = 2) - ) + labs( + title = + paste( + "Before Optimization, \n RMSE=", + round(stats$RMSE, digits = 2) + ) ) + theme(plot.title = element_text(size = 9, hjust = 0.5)) @@ -249,11 +246,12 @@ p <- p2 <- p[[1]] + - labs(title = - paste( - "After Optimization, \n RMSE=", - round(stats$RMSE, digits = 2) - ) + labs( + title = + paste( + "After Optimization, \n RMSE=", + round(stats$RMSE, digits = 2) + ) ) + theme(plot.title = element_text(size = 9, hjust = 0.5)) @@ -266,7 +264,6 @@ ggsave( ), plot = p ) - ``` This gives: diff --git a/vignettes/Parameter_estimation_simple_case.Rmd b/vignettes/Parameter_estimation_simple_case.Rmd index a5ac701..5c37103 100644 --- a/vignettes/Parameter_estimation_simple_case.Rmd +++ b/vignettes/Parameter_estimation_simple_case.Rmd @@ -60,7 +60,6 @@ if (!require("SticsRPacks")) { ``` ```{r setup_initializations, eval=params$eval_auto_test, message=FALSE, results=FALSE, warning=FALSE} - # DEFINE THE PATH TO YOUR LOCALLY INSTALLED VERSION OF JAVASTICS javastics_path <- path_to_JavaStics @@ -109,7 +108,6 @@ Here model parameters values are read in the model input files. ```{r eval=params$eval_auto_test, results='hide', message=FALSE, warning=FALSE} - # Set the model options (see '? stics_wrapper_options' for details) model_options <- stics_wrapper_options( @@ -135,11 +133,10 @@ In variables and parameters names, "(\*)" must be replaced by "_\*" to be handle ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} - -sit_name <- "bo96iN+" # can be a vector of situation names if you want to - # consider several, e.g. c("bo96iN+","bou00t1") -var_name <- "lai_n" # can be a vector of variable names if you want to - # consider several, e.g. c("lai_n","masec_n") +sit_name <- "bo96iN+" # can be a vector of situation names if you want to +# consider several, e.g. c("bo96iN+","bou00t1") +var_name <- "lai_n" # can be a vector of variable names if you want to +# consider several, e.g. c("lai_n","masec_n") obs_list <- get_obs(javastics_workspace_path, usm = sit_name) obs_list <- filter_obs(obs_list, var = var_name, include = TRUE) ``` @@ -178,19 +175,19 @@ The number of repetitions is advised to be set to at least 5, while 10 is a reas ```{r eval=params$eval_auto_test, message=FALSE, warning=FALSE} optim_options <- list() optim_options$nb_rep <- 7 # Number of repetitions of the minimization - # (each time starting with different initial - # values for the estimated parameters) +# (each time starting with different initial +# values for the estimated parameters) optim_options$maxeval <- 500 # Maximum number of evaluations of the - # minimized criteria +# minimized criteria optim_options$xtol_rel <- 1e-03 # Tolerance criterion between two iterations - # (threshold for the relative difference of - # parameter values between the 2 previous - # iterations) +# (threshold for the relative difference of +# parameter values between the 2 previous +# iterations) optim_options$out_dir <- data_dir # path where to store the results - # (graph and Rdata) +# (graph and Rdata) optim_options$ranseed <- 1234 # set random seed so that each execution give the - # same results - # If you want randomization, don't set it. +# same results +# If you want randomization, don't set it. ``` @@ -291,7 +288,6 @@ This data is also stored in the `optim_options$out_dir` folder (file optim_resul ```{r eval=params$eval_auto_vignette, echo=FALSE, out.width = '45%'} - knitr::include_graphics("ResultsSimpleCase/estimInit_dlaimax.PNG") knitr::include_graphics("ResultsSimpleCase/estimInit_durvieF.PNG") @@ -303,13 +299,11 @@ Figure 1: plots of estimated vs initial values of parameters dlaimax and durvieF * the ValuesVSit.pdf file contains: ```{r eval=params$eval_auto_vignette, echo=FALSE, out.width = '45%', fig.show = "hold", fig.align="default"} - knitr::include_graphics("ResultsSimpleCase/critVSit.PNG") knitr::include_graphics("ResultsSimpleCase/dlaimaxVSit.PNG") knitr::include_graphics("ResultsSimpleCase/durvieFVSit.PNG") - ``` Figure 2: The top left figure displays the evolution of the minimized criterion value in function of the iterations of the minimization. Colors and labels represent the repetition number. The top right figure displays the evolution of the value of the dlaimax parameter in function of the iterations of the minimization. Labels represent the repetition number and colors the value of the minimized criterion. The figure at the bottom is the same but for the durvieF parameter. @@ -379,7 +373,6 @@ ggsave( ), plot = p ) - ```