From 23f62bdf158434182fc0640558798a588b31b74c Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Fri, 21 Nov 2025 10:41:13 +0100 Subject: [PATCH 01/16] DEoptim and graDiEnt packages integration --- DESCRIPTION | 4 +- NAMESPACE | 5 + R/global_optim_function.R | 529 ++++++++++++++++++++++ R/optim_switch.R | 91 +++- R/wrap_DE.R | 100 ++++ R/wrap_graDiEnt.R | 136 ++++++ man/plot_global_optim.Rd | 53 +++ man/plot_valuesVSit_2D_go.Rd | 79 ++++ man/plot_valuesVSit_go.Rd | 72 +++ man/post_treat_global_optim.Rd | 55 +++ man/summary_global_optim.Rd | 52 +++ man/wrap_DEoptim.Rd | 33 ++ man/wrap_graDiEnt.Rd | 25 + tests/testthat/test_estim_param_DEoptim.R | 176 +++++++ tests/testthat/test_estim_param_GEoptim.R | 179 ++++++++ 15 files changed, 1563 insertions(+), 26 deletions(-) create mode 100644 R/global_optim_function.R create mode 100644 R/wrap_DE.R create mode 100644 R/wrap_graDiEnt.R create mode 100644 man/plot_global_optim.Rd create mode 100644 man/plot_valuesVSit_2D_go.Rd create mode 100644 man/plot_valuesVSit_go.Rd create mode 100644 man/post_treat_global_optim.Rd create mode 100644 man/summary_global_optim.Rd create mode 100644 man/wrap_DEoptim.Rd create mode 100644 man/wrap_graDiEnt.Rd create mode 100644 tests/testthat/test_estim_param_DEoptim.R create mode 100644 tests/testthat/test_estim_param_GEoptim.R diff --git a/DESCRIPTION b/DESCRIPTION index a2467948..84e7ed9f 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -36,9 +36,11 @@ Imports: Matrix, BayesianTools, crayon, + DEoptim, doParallel, dplyr, ggplot2, + graDiEnt, lhs, lifecycle, lubridate, @@ -76,4 +78,4 @@ Encoding: UTF-8 Language: en-US LazyData: true Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index 165dc471..fec65822 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -14,6 +14,8 @@ export(likelihood_log_ciidn_corr) export(plot_estimVSinit) export(plot_valuesVSit) export(plot_valuesVSit_2D) +export(plot_valuesVSit_2D_go) +export(plot_valuesVSit_go) export(test_wrapper) importFrom(BayesianTools,MAP) importFrom(BayesianTools,applySettingsDefault) @@ -24,6 +26,7 @@ importFrom(BayesianTools,gelmanDiagnostics) importFrom(BayesianTools,getSample) importFrom(BayesianTools,marginalPlot) importFrom(BayesianTools,runMCMC) +importFrom(DEoptim,DEoptim) importFrom(dplyr,"%>%") importFrom(dplyr,filter) importFrom(dplyr,select) @@ -45,6 +48,8 @@ importFrom(ggplot2,scale_y_log10) importFrom(ggplot2,theme) importFrom(ggplot2,xlim) importFrom(ggplot2,ylim) +importFrom(graDiEnt,GetAlgoParams) +importFrom(graDiEnt,optim_SQGDE) importFrom(purrr,modify) importFrom(rlang,.data) importFrom(stats,setNames) diff --git a/R/global_optim_function.R b/R/global_optim_function.R new file mode 100644 index 00000000..8785169d --- /dev/null +++ b/R/global_optim_function.R @@ -0,0 +1,529 @@ +#' @title Summarizes results of global optim methods +#' +#' @inheritParams estim_param +#' +#' @param optim_results Results list returned by global optim method wrappers +#' +#' @return Prints results of global optim methods +#' +summary_global_optim<- function(optim_options, param_info, optim_results, out_dir) { + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + est_values <- optim_results$est_values + ind_min_crit <- optim_results$ind_min_crit + min_crit_value <- optim_results$min_crit_value + + cat(paste( + "\nList of observed variables used:", + paste(optim_results$obs_var_list, collapse = ", ") + )) + + # Display of parameters values for the candidate which has the + # smallest criterion + for (ipar in 1:nb_params) { + cat(paste( + "\nEstimated value for", param_names[ipar], ": ", + format(est_values[ind_min_crit, ipar], + scientific = FALSE, + digits = 2, nsmall = 0 + ) + )) + } + cat(paste( + "\nMinimum value of the criterion:", + format(min_crit_value, scientific = FALSE, digits = 2, nsmall = 0) + )) + cat(paste( + "\nComplementary graphs and results can be found in ", out_dir, + "\n" + )) +} + + +#' @title Post-treat results of global optim methods +#' +#' @inheritParams estim_param +#' +#' @param optim_results Results list returned by global optim method wrappers +#' @param crit_options List containing several arguments given to `estim_param` +#' function: `param_names`, `obs_list`, `model_function`, +#' `model_options`, `param_info`, `transform_obs`, `transform_sim` +#' that must be passed to main_crit function by the methods wrappers. +#' +#' @return Updated results of global optim method +#' +post_treat_global_optim <- function(optim_options, param_info, optim_results, + crit_options) { + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + info_crit_list <- crit_options$info_crit_list + + # Recompute final value of minimized criterion + # (just to check it is correct and to get the observation list used) + info_final <- main_crit( + param_values = optim_results$final_values, + crit_options = c(crit_options, return_detailed_info = TRUE) + ) + if (info_final$crit != optim_results$min_crit_value) { + stop(paste( + "Internal error: incoherent computation of minimum criterion value. \nValue obtained in method wrapper:", + optim_results$min_crit_value, "\nValue obtained afterwards:", + info_final$crit + )) + } + optim_results$forced_param_values <- info_final$forced_param_values + + sapply(info_crit_list, function(x) { + final_info_crit <- x( + obs_list = info_final$obs_intersect, + crit = info_final$crit, + param_nb = nb_params + ) + optim_results[x()$name] <<- final_info_crit + }) + + return(optim_results) +} + + +#' @title Generate plots for global optim methods +#' +#' @inheritParams estim_param +#' +#' @param optim_results Results list returned by global optim method wrappers +#' +#' @return Returns the list of plots + save them in a pdf file. +#' +#' @keywords internal +#' +plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) { + bounds <- get_params_bounds(param_info) + init_values <- optim_results$init_values + est_values <- optim_results$est_values + #crit_values <- optim_results$crit_values + p_all <- list() + + # ValuesVSit plot + + tryCatch( + { + grDevices::pdf( + file = file.path(out_dir, "ValuesVSit_globaloptim2.pdf"), + width = 9, height = 9 + ) + }, + error = function(cond) { + filename <- paste0("ValuesVSit_new.pdf") + warning( + "Error trying to create ", out_dir, + "/ValuesVSit.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", + filename + ) + grDevices::pdf( + file = file.path(out_dir, filename), + width = 9, height = 9 + ) + } + ) + + # DEoptim :we should construct a dataframe df "individu × itération" + if (!is.null(optim_results$DE)) { + + DE <- optim_results$DE + param_names <- get_params_names(param_info) + + # storepopfrom = 1 should be activated id wrap_DE + storepop <- DE$member$storepop + bestvalit <- DE$member$bestvalit + + if (length(storepop) > 0) { + itermax <- length(storepop) + NP <- nrow(storepop[[1]]) + + df_list <- vector("list", itermax) + + for (it in seq_len(itermax)) { + pop_it <- as.data.frame(storepop[[it]]) + colnames(pop_it) <- param_names + + pop_it$iter <- it # it id + pop_it$ind <- seq_len(NP) # id individual + # ------------------------------------------------ + ### To see with Samuel + # compute criterion per individual (To proprely plot the last figure of ValueVSit) (I have to add in parameters crit_options) + #crit_pop <- apply(storepop[[it]], 1, function(par) { + # main_crit(par, c(crit_options, return_detailed_info = FALSE))$crit + #}) + #pop_it$crit <- crit_pop + # ------------------------------------------------ + pop_it$crit <- bestvalit[it] + df_list[[it]] <- pop_it + } + + df_DE <- dplyr::bind_rows(df_list) + + # for plot_valuesVSit, we need eval col + df_DE$eval <- df_DE$iter + + p <- plot_valuesVSit_go( + df = df_DE, + param_info = param_info, + iter_or_eval = "iter", + crit_log = FALSE, + ind_label = "begin" + ) + + for (plot in p) print(plot) + grDevices::dev.off() + p_all$valuesVSit_DE <- p + } else { + warning("DE$member$storepop is empty : verify storepopfrom in DEoptim.control") + grDevices::dev.off() + } + } + + # ValuesVSit_2D plot + + if (!is.null(est_values)) { + + # Build a data.frame with final parameters + df_2d <- as.data.frame(est_values) + colnames(df_2d) <- get_params_names(param_info) + df_2d$eval <- seq_len(nrow(df_2d)) + df_2d$iter <- df_2d$eval + + tryCatch( + { + grDevices::pdf( + file = file.path(out_dir, "ValuesVSit_2D.pdf"), + width = 9, height = 9 + ) + }, + error = function(cond) { + filename <- paste0("ValuesVSit_2D_new.pdf") + warning( + "Error trying to create ", out_dir, + "/ValuesVSit_2D.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", + filename + ) + grDevices::pdf( + file = file.path(out_dir, filename), + width = 9, height = 9 + ) + } + ) + + + p <- plot_valuesVSit_2D_go( + df = df_2d, + param_info = param_info, + iter_or_eval = "eval", + fill = "eval", + crit_log = FALSE, + lines = FALSE, + rep_label = "end" + ) + + for (plot in p) { + print(plot) + } + grDevices::dev.off() + p_all$valuesVSit_2D <- p + } + + + + return(p_all) +} + + +#' @title Create plots of parameters and criterion values per iteration or +#' evaluation number +#' +#' @inheritParams estim_param +#' @param df Data.frame containing values of each individual and iteration parameters values (one column per +#' estimated parameter), criterion (crit column), individual index (ind), +#' iteration number (iter) and evaluation number (eval) +#' (similar to params_and_crit). +#' See Details section for comments about the difference between evaluations +#' and iterations. +#' @param iter_or_eval Values of the x axis: "iter" for iteration number, +#' "eval" for evaluation number +#' @param crit_log If TRUE, consider criterion values in log scale +#' @param ind_label Indicate if labels for the individual number must be +#' plotted at both beginning and end of lines ("begin_end"), +#' only at the beginning ("begin") or only at the end ("end") +#' +#' @return A named list containing one plot per parameter and a plot for the +#' criterion. +#' +#' @details Evaluation means evaluation of the criterion from proposed values of +#' the parameters by the parameter estimation algorithm. +#' An iteration is reached when an evaluation lead to a better value of the +#' criterion than the previously obtained values. +#' There are thus more evaluations than iterations. The criterion decreases when +#' iteration number increases while it is not the case when evaluation number +#' increases. +#' +#' @importFrom ggplot2 ggplot aes_string theme element_text geom_point +#' scale_color_gradient2 geom_line geom_label aes labs scale_y_log10 +#' @importFrom dplyr select filter %>% +#' +#' @export +#' +plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), + crit_log = TRUE, + ind_label = c("begin_end", "begin", "end")) { + param_names <- get_params_names(param_info) + bounds <- get_params_bounds(param_info) + + lab <- "evaluations" + if (iter_or_eval[1] == "iter") { + df <- filter(df, !is.na(.data$iter)) + lab <- "iterations" + } + trans <- "identity" + mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) + if (crit_log) { + if (all(df$crit > 0)) { + trans <- "log10" + mid <- (max(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 + } + } + tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) + tmp[tmp == Inf | tmp == -Inf] <- NA + minvalue <- apply(tmp, 2, min, na.rm = TRUE) + maxvalue <- apply(tmp, 2, max, na.rm = TRUE) + minvalue <- minvalue - 0.05 * (maxvalue - minvalue) + maxvalue <- maxvalue + 0.05 * (maxvalue - minvalue) + + p <- list() + + for (param_name in param_names) { + p[[param_name]] <- ggplot(df, aes( + x = !!rlang::sym(iter_or_eval[1]), + y = !!rlang::sym(param_name), + color = crit + )) + + labs( + title = paste0( + "Evolution of ", param_name, + " \n in function of the minimization ", lab + ), + y = param_name, + x = paste(lab, "number"), + fill = "Criterion" + ) + + theme(plot.title = element_text(hjust = 0.5)) + + geom_point(alpha = 0.5) + + scale_color_gradient2( + midpoint = mid, low = "blue", mid = "yellow", + high = "red", space = "Lab", trans = trans + ) + + for (iind in unique(df$ind)) { + p[[param_name]] <- p[[param_name]] + + geom_line(data = filter(df, ind == iind)) + if (ind_label[1] == "begin_end" || ind_label[1] == "begin") { + p[[param_name]] <- p[[param_name]] + + geom_label(aes(label = ind), + data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), + size = 3 + ) + } + if (ind_label[1] == "begin_end" || ind_label[1] == "end") { + p[[param_name]] <- p[[param_name]] + + geom_label(aes(label = ind), + data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)), + size = 3 + ) + } + } + ylim(minvalue[param_name], maxvalue[param_name]) + } + + df$ind <- as.factor(df$ind) + p[["criterion"]] <- ggplot(df, aes_string( + x = iter_or_eval[1], y = "crit", + color = "ind" + )) + + labs( + title = paste0( + "Evolution of the minimized criterion \n in function of the minimization ", + lab + ), + y = "Minimized criterion", + x = paste(lab, "number"), + fill = "Individual" + ) + + theme(plot.title = element_text(hjust = 0.5)) + + geom_point(alpha = 0.5) + + for (iind in unique(df$ind)) { + p[["criterion"]] <- p[["criterion"]] + + geom_line(data = filter(df, ind == iind)) + + geom_label(aes(label = ind), + data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), + size = 3 + ) + + geom_label(aes(label = ind), + data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)) + ) + } + + if (crit_log) { + p[["criterion"]] <- p[["criterion"]] + scale_y_log10() + } + + return(p) +} + + +#' @title Create 2D plots of parameters values evolution per iteration or +#' evaluation number +#' +#' @inheritParams estim_param +#' @param df Data.frame containing values of parameters (one column per +#' estimated parameter), criterion (crit column), repetition number (rep), +#' iteration number (iter) and evaluation number (eval) +#' (similar to params_and_crit). +#' See Details section for comments about the difference between evaluations +#' and iterations. +#' @param iter_or_eval "iter" for plotting the values for each iteration, +#' "eval" for plotting the values for each evaluation +#' @param fill If "crit", colours the points and lines in function of the +#' minimized criterion value, if "rep" colours in function of the +#' repetition number. +#' @param crit_log If TRUE, consider criterion values in log scale +#' @param lines If TRUE add lines between points of a same repetition +#' @param rep_label Indicate if labels for the repetition number must be plotted +#' at both beginning and end of lines ("begin_end"), only at the beginning +#' ("begin") or only at the end ("end") +#' +#' @return A list containing one plot per parameter pair. +#' +#' @details Evaluation means evaluation of the criterion from proposed values of +#' the parameters by the parameter estimation algorithm. +#' An iteration is reached when an evaluation lead to a better value of the +#' criterion than the previously obtained values. +#' There are thus more evaluations than iterations. The criterion decreases when +#' iteration number increases while it is not the case when evaluation number +#' increases. +#' +#' @importFrom ggplot2 ggplot aes_string theme element_text geom_point labs +#' xlim ylim geom_path scale_y_log10 +#' @importFrom dplyr select filter %>% +#' +#' @export +#' +plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "iter"), + fill = c("crit", "rep"), crit_log = TRUE, + lines = FALSE, + rep_label = c("begin_end", "begin", "end")) { + param_names <- get_params_names(param_info) + if (length(param_names) <= 1) { + return() + } + bounds <- get_params_bounds(param_info) + + lab <- "evaluations" + if (iter_or_eval[1] == "iter") { + df <- filter(df, !is.na(.data$iter)) + lab <- "iterations" + } + if ("rep" %in% names(df)) { + df$rep <- as.factor(df$rep) + } + has_crit <- "crit" %in% names(df) + trans <- "identity" + mid <- NA_real_ + if (fill[1] == "crit" && has_crit) { + mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) + if (crit_log) { + if (all(df$crit > 0)) { + trans <- "log10" + mid <- (max(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 + } + } + } else { + crit_log <- FALSE + } + tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) + # -.data$ avoid NOTES on check ... + tmp[tmp == Inf | tmp == -Inf] <- NA + minvalue <- apply(tmp, 2, min, na.rm = TRUE) + maxvalue <- apply(tmp, 2, max, na.rm = TRUE) + minvalue <- minvalue - 0.05 * (maxvalue - minvalue) + maxvalue <- maxvalue + 0.05 * (maxvalue - minvalue) + + p <- list() + + df_pairs <- utils::combn(param_names, 2) + + for (ipair in seq_len(ncol(df_pairs))) { + p[[ipair]] <- ggplot(df, aes_string( + x = df_pairs[1, ipair], + y = df_pairs[2, ipair], color = fill[1] + )) + + labs( + title = paste0( + "Evolution of ", df_pairs[1, ipair], " and ", + df_pairs[2, ipair], " \n in function of the minimization ", + lab + ), + y = paste("Estimated value for", df_pairs[2, ipair]), + x = paste("Estimated value for", df_pairs[1, ipair]), + fill = "Criterion" + ) + + theme(plot.title = element_text(hjust = 0.5)) + + geom_point(alpha = 0.5) + + if (fill[1] == "crit" && has_crit){ + p[[ipair]] <- p[[ipair]] + + scale_color_gradient2( + midpoint = mid, low = "blue", mid = "yellow", + high = "red", space = "Lab", trans = trans + ) + } else if (fill[1] == "eval") { + # palette continue simple pour eval + p[[ipair]] <- p[[ipair]] + + scale_color_gradient2( + low = "lightblue", + high = "darkblue" + ) + } + + if (lines) { + for (irep in unique(df$rep)) { + p[[ipair]] <- p[[ipair]] + + geom_path(data = filter(df, rep == irep)) + 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)), + 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)), + size = 3 + ) + } + } + } + xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) + ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]]) + } + + return(p) +} diff --git a/R/optim_switch.R b/R/optim_switch.R index 2e19e32b..1da98bce 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -23,7 +23,11 @@ optim_switch <- function(...) { flag_error <- FALSE on.exit({ - res$forced_param_values <- crit_options$forced_param_values + #### + if (!is.null(crit_options)) { + res$forced_param_values <- crit_options$forced_param_values + } + ##### if (!is.null(res$final_values)) { res$forced_param_values <- compute_eq_const( res$forced_param_values, @@ -68,7 +72,8 @@ optim_switch <- function(...) { } } - if (!is.null(crit_options$out_dir) & length(res) > 0) { + #if (!is.null(crit_options$out_dir) & length(res) > 0) { + if (!is.null(crit_options) && !is.null(crit_options$out_dir) && length(res) > 0) { save(res, file = file.path( crit_options$out_dir, "optim_results.Rdata" @@ -91,6 +96,19 @@ optim_switch <- function(...) { optim_method <- arguments$optim_method optim_options <- arguments$optim_options wrap_args <- within(arguments, rm("optim_method")) + ########################## + # + # param_info <- NULL + #crit_options <- NULL + + # if ("param_info" %in% names(arguments)) { + # param_info <- arguments$param_info + # } + # if ("crit_options" %in% names(arguments)) { + # crit_options <- arguments$crit_options + # } + ########################## + if (nargs() > 2) { param_info <- arguments$param_info crit_options <- arguments$crit_options @@ -112,7 +130,7 @@ optim_switch <- function(...) { optim_results = res, crit_options = crit_options ) - res$plots <- plot_frequentist( + res$plots <- plot_frequentist( optim_options = optim_options, param_info = param_info, optim_results = res, @@ -132,7 +150,8 @@ optim_switch <- function(...) { res$obs_situation_list <- .croptEnv$obs_situation_list res$plots <- plot_bayesian( optim_options = optim_options, - param_info = param_info, optim_results = res, + param_info = param_info, + optim_results = res, out_dir = crit_options$out_dir ) summary_bayesian( @@ -141,32 +160,54 @@ optim_switch <- function(...) { out_dir = crit_options$out_dir ) } - } else if (optim_method == "optim") { - res <- do.call(wrap_optim, wrap_args) + } else if (optim_method == "DEoptim") { + res <- do.call(wrap_DEoptim, wrap_args) if (nargs() > 2) { res$obs_var_list <- .croptEnv$obs_var_list res$obs_situation_list <- .croptEnv$obs_situation_list if (arguments$crit_options$info_level >= 1) { res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) } - res <- post_treat_frequentist( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - crit_options = crit_options - ) - res$plots <- plot_frequentist( - optim_options = optim_options, - param_info = param_info, optim_results = res, - out_dir = crit_options$out_dir - ) - summary_frequentist( - optim_options = optim_options, param_info = param_info, - optim_results = res, - out_dir = crit_options$out_dir - ) - } - } else { + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + summary_global_optim( + optim_options = optim_options, param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) +} + } else if (optim_method == "graDiEnt") { + res <- do.call(wrap_graDiEnt, wrap_args) + if (nargs() > 2) { + res$obs_var_list <- .croptEnv$obs_var_list + res$obs_situation_list <- .croptEnv$obs_situation_list + if (arguments$crit_options$info_level >= 1) { + res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) + } + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + + }}else { flag_unknown_method <- TRUE }, error = function(cond) { @@ -179,7 +220,7 @@ optim_switch <- function(...) { flag_error <- TRUE stop(paste0( "Unknown method ", optim_method, - ", please choose between nloptr.simplex, BayesianTools.dreamzs and optim." + ", please choose between nloptr.simplex, BayesianTools.dreamzs, DEoptim and graDiEnt." )) } diff --git a/R/wrap_DE.R b/R/wrap_DE.R new file mode 100644 index 00000000..975db799 --- /dev/null +++ b/R/wrap_DE.R @@ -0,0 +1,100 @@ +#' @title A wrapper for DEoptim package +#' +#' @inheritParams optim_switch +#' +#' @return prints, graphs and a list containing: +#' `final_values`, the vector of estimated values for optimized parameters +#' for the repetition that lead to the lowest value of the criterion +#' `init_values`, the vector of initial values for optimized parameters +#' `est_values`, data frame of values of all population +#' `ind_min_crit`, the index of the repetition that lead to the lowest value +#' of the criterion +#' `DE`, the data structure returned by DE +#' `crit_values` the vector of criterion values for the all population +#' @importFrom DEoptim DEoptim +#' @keywords internal + +wrap_DEoptim <- function(optim_options, param_info, crit_options) { + if(is.null((ranseed <- optim_options$ranseed))) { + ranseed = NULL + } + if (is.null((reltol <- optim_options$reltol))) { + optim_options$reltol <- 1e-10 + } + if (is.null((NP <- optim_options$NP))) { + stop( + "Optim_options$NP is mandatory, please define it (e.g., 10 * number of parameters)" + ) + } + if (is.null((trace <- optim_options$trace))) { + optim_options$trace <- FALSE + } + optim_options$ranseed = NULL # ranseed is set to NULL because DEoptim.controll doesn't regonise it + optim_options$storepopfrom <- 1 + optim_options$storepopfreq <- 1 + control <- do.call(DEoptim::DEoptim.control, optim_options) + algorithm <- "DEoptim" + + # return requested information if only optim_options is given in argument + if (nargs() == 1 & methods::hasArg(optim_options)) { + return(list( + package = "DEoptim", family = "Global", + method = "DEoptim", init_values_nb = control$NP + )) + } + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + + crit_options$tot_max_eval <- control$NP * control$itermax + bounds <- get_params_bounds(param_info) + init_values <- get_init_values(param_info) + control$initialpop <- as.matrix(init_values) + + + start_time <- Sys.time() + + if (!is.null(ranseed)) set.seed(ranseed) + fn_de <- function(x) { + main_crit(x, crit_options = crit_options) + } + DE <- tryCatch( + DEoptim::DEoptim( + fn = fn_de, + lower = bounds$lb, + upper = bounds$ub, + control = control + ), + error = function(e) { warning(sprintf("DEoptim failed: %s", e$message)); NULL } + ) + elapsed <- Sys.time() - start_time + + # Verify criterion value + if (is.na( DE$optim$bestval)) { + stop( + "The value of the criterion is NA." + ) + } + + + # Get the estimated values ==> matrix of NP * nb_params + est_values <- NULL + est_values <- as.matrix(DE$member$pop) + colnames(est_values) <- param_names + if (ncol(est_values) == nb_params) { + colnames(est_values) <- param_names + } + + # Final solution + final_values <- DE$optim$bestmem + names(final_values) <- param_names + + res <- list( + final_values = final_values, + init_values = init_values, + est_values = est_values, + min_crit_value = DE$optim$bestval, + DE = DE + ) + + return(res) +} diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R new file mode 100644 index 00000000..def4231c --- /dev/null +++ b/R/wrap_graDiEnt.R @@ -0,0 +1,136 @@ +#' @title A wrapper for graDiEnt package +#' +#' @inheritParams optim_switch +#' +#' @return list with final_values, init_values, est_values, GE, min_crit_value +#' @importFrom graDiEnt optim_SQGDE GetAlgoParams +#' @keywords internal +wrap_graDiEnt <- function(optim_options, param_info, crit_options) { + if(is.null((ranseed <- optim_options$ranseed))) { + ranseed = NULL + } + if (is.null((reltol <- optim_options$reltol))) { + optim_options$converge_crit <- "stdev" + } + if (is.null((n_params <- optim_options$n_params))) { + stop( + "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)" + ) + } + if (is.null((trace <- optim_options$return_trace))) { + optim_options$return_trace <- TRUE + } + optim_options$ranseed = NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it + algorithm <- "graDiEnt" + + # return requested information if only optim_options is given in argument + if (nargs() == 1 && methods::hasArg(optim_options)) { + init_nb <- if (!is.null(optim_options$n_particles)) { + optim_options$n_particles + } else { + NA_integer_ + } + return(list( + package = "graDiEnt", family = "Global", + method = "SQGDE", init_values_nb = init_nb + )) + } + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + bounds <- get_params_bounds(param_info) + init_values <- get_init_values(param_info) + + control_params <- do.call(graDiEnt::GetAlgoParams,optim_options) + + crit_options$tot_max_eval <- control_params$n_particles * control_params$n_iter + + start_time <- Sys.time() + + if (!is.null(ranseed)) set.seed(ranseed) + range_bounds <- bounds$ub - bounds$lb + ObjFun <- function(u) { + u <- as.numeric(u) + names(u) <- param_names + x <- bounds$lb + u * range_bounds + names(x) <- param_names + main_crit(x, crit_options = crit_options) + } + SQGDE <- tryCatch( + graDiEnt::optim_SQGDE( + ObjFun = ObjFun, + control_params = control_params + ), + error = function(e) { warning(sprintf("GraDiEnt failed: %s", e$message)); NULL } + ) + elapsed <- Sys.time() - start_time + + # Verify criterion value + if (is.null(SQGDE) || is.null(SQGDE$solution)) { + stop( + "The value of the criterion is NA." + ) + } + u_sol <- SQGDE$solution + names(u_sol) <- param_names + + final_values <- bounds$lb + u_sol * range_bounds + names(final_values) <- param_names + min_crit_value <- SQGDE$solution_weight + if (is.null(min_crit_value) ) { + min_crit_value <- main_crit(final_values, crit_options) + } + + + + # Get the estimated values ==> matrix of 1 * nb_params + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) + colnames(est_values) <- param_names + if (ncol(est_values) == nb_params) { + colnames(est_values) <- param_names + } + + # Final solution + if (!is.null(SQGDE$particles_trace)) { + tr <- SQGDE$particles_trace + if (length(dim(tr)) != 3) { + warning("Unexpected dimension for SQGDE$particles_trace; est_values set to final_values only.") + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) + } else { + n_it <- dim(tr)[1] + last_pop <- tr[n_it, , ] # matrice n_particles × n_params (en u) + est_u <- as.matrix(last_pop) + colnames(est_u) <- param_names + + # transformation vers l’échelle physique + est_values <- sweep(est_u, 2, range_bounds, `*`) + est_values <- sweep(est_values, 2, bounds$lb, `+`) + colnames(est_values) <- param_names + } + } else { + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) + } + + + res <- list( + final_values = final_values, + init_values = init_values, + est_values = est_values, + min_crit_value = min_crit_value, + SQGDE = SQGDE + ) + + return(res) +} + diff --git a/man/plot_global_optim.Rd b/man/plot_global_optim.Rd new file mode 100644 index 00000000..2f2d32f5 --- /dev/null +++ b/man/plot_global_optim.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/global_optim_function.R +\name{plot_global_optim} +\alias{plot_global_optim} +\title{Generate plots for global optim methods} +\usage{ +plot_global_optim(optim_options, param_info, optim_results, out_dir) +} +\arguments{ +\item{optim_options}{List of options of the parameter estimation method, containing: +\itemize{ +\item \code{ranseed} Set random seed so that each execution of estim_param give the same +results when using the same seed. If you want randomization, set it to NULL, +otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed) +\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex} +and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods. +\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead. +}} + +\item{param_info}{Information on the parameters to estimate. +Either +a list containing: +\itemize{ +\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), +\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) +\item \code{init_values}, a data.frame containing initial +values to test for the parameters (optional, if not provided, or if less values +than number of repetitions of the minimization are provided, the, or part +of the, initial values will be randomly generated using LHS sampling within +parameter bounds). +} + +or a named list containing for each parameter: +\itemize{ +\item \code{sit_list}, list the groups of situations for which the current estimated +parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} +for an example), +\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), +\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). +\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) +}} + +\item{optim_results}{Results list returned by global optim method wrappers} + +\item{out_dir}{Path to the directory where the optimization results will be written. (optional, default to \code{getwd()})} +} +\value{ +Returns the list of plots + save them in a pdf file. +} +\description{ +Generate plots for global optim methods +} +\keyword{internal} diff --git a/man/plot_valuesVSit_2D_go.Rd b/man/plot_valuesVSit_2D_go.Rd new file mode 100644 index 00000000..0ee03ced --- /dev/null +++ b/man/plot_valuesVSit_2D_go.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/global_optim_function.R +\name{plot_valuesVSit_2D_go} +\alias{plot_valuesVSit_2D_go} +\title{Create 2D plots of parameters values evolution per iteration or +evaluation number} +\usage{ +plot_valuesVSit_2D_go( + df, + param_info, + iter_or_eval = c("eval", "iter"), + fill = c("crit", "rep"), + crit_log = TRUE, + lines = FALSE, + rep_label = c("begin_end", "begin", "end") +) +} +\arguments{ +\item{df}{Data.frame containing values of parameters (one column per +estimated parameter), criterion (crit column), repetition number (rep), +iteration number (iter) and evaluation number (eval) +(similar to params_and_crit). +See Details section for comments about the difference between evaluations +and iterations.} + +\item{param_info}{Information on the parameters to estimate. +Either +a list containing: +\itemize{ +\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), +\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) +\item \code{init_values}, a data.frame containing initial +values to test for the parameters (optional, if not provided, or if less values +than number of repetitions of the minimization are provided, the, or part +of the, initial values will be randomly generated using LHS sampling within +parameter bounds). +} + +or a named list containing for each parameter: +\itemize{ +\item \code{sit_list}, list the groups of situations for which the current estimated +parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} +for an example), +\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), +\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). +\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) +}} + +\item{iter_or_eval}{"iter" for plotting the values for each iteration, +"eval" for plotting the values for each evaluation} + +\item{fill}{If "crit", colours the points and lines in function of the +minimized criterion value, if "rep" colours in function of the +repetition number.} + +\item{crit_log}{If TRUE, consider criterion values in log scale} + +\item{lines}{If TRUE add lines between points of a same repetition} + +\item{rep_label}{Indicate if labels for the repetition number must be plotted +at both beginning and end of lines ("begin_end"), only at the beginning +("begin") or only at the end ("end")} +} +\value{ +A list containing one plot per parameter pair. +} +\description{ +Create 2D plots of parameters values evolution per iteration or +evaluation number +} +\details{ +Evaluation means evaluation of the criterion from proposed values of +the parameters by the parameter estimation algorithm. +An iteration is reached when an evaluation lead to a better value of the +criterion than the previously obtained values. +There are thus more evaluations than iterations. The criterion decreases when +iteration number increases while it is not the case when evaluation number +increases. +} diff --git a/man/plot_valuesVSit_go.Rd b/man/plot_valuesVSit_go.Rd new file mode 100644 index 00000000..9034cbbd --- /dev/null +++ b/man/plot_valuesVSit_go.Rd @@ -0,0 +1,72 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/global_optim_function.R +\name{plot_valuesVSit_go} +\alias{plot_valuesVSit_go} +\title{Create plots of parameters and criterion values per iteration or +evaluation number} +\usage{ +plot_valuesVSit_go( + df, + param_info, + iter_or_eval = c("iter", "eval"), + crit_log = TRUE, + ind_label = c("begin_end", "begin", "end") +) +} +\arguments{ +\item{df}{Data.frame containing values of each individual and iteration parameters values (one column per +estimated parameter), criterion (crit column), individual index (ind), +iteration number (iter) and evaluation number (eval) +(similar to params_and_crit). +See Details section for comments about the difference between evaluations +and iterations.} + +\item{param_info}{Information on the parameters to estimate. +Either +a list containing: +\itemize{ +\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), +\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) +\item \code{init_values}, a data.frame containing initial +values to test for the parameters (optional, if not provided, or if less values +than number of repetitions of the minimization are provided, the, or part +of the, initial values will be randomly generated using LHS sampling within +parameter bounds). +} + +or a named list containing for each parameter: +\itemize{ +\item \code{sit_list}, list the groups of situations for which the current estimated +parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} +for an example), +\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), +\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). +\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) +}} + +\item{iter_or_eval}{Values of the x axis: "iter" for iteration number, +"eval" for evaluation number} + +\item{crit_log}{If TRUE, consider criterion values in log scale} + +\item{ind_label}{Indicate if labels for the individual number must be +plotted at both beginning and end of lines ("begin_end"), +only at the beginning ("begin") or only at the end ("end")} +} +\value{ +A named list containing one plot per parameter and a plot for the +criterion. +} +\description{ +Create plots of parameters and criterion values per iteration or +evaluation number +} +\details{ +Evaluation means evaluation of the criterion from proposed values of +the parameters by the parameter estimation algorithm. +An iteration is reached when an evaluation lead to a better value of the +criterion than the previously obtained values. +There are thus more evaluations than iterations. The criterion decreases when +iteration number increases while it is not the case when evaluation number +increases. +} diff --git a/man/post_treat_global_optim.Rd b/man/post_treat_global_optim.Rd new file mode 100644 index 00000000..16900d51 --- /dev/null +++ b/man/post_treat_global_optim.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/global_optim_function.R +\name{post_treat_global_optim} +\alias{post_treat_global_optim} +\title{Post-treat results of global optim methods} +\usage{ +post_treat_global_optim(optim_options, param_info, optim_results, crit_options) +} +\arguments{ +\item{optim_options}{List of options of the parameter estimation method, containing: +\itemize{ +\item \code{ranseed} Set random seed so that each execution of estim_param give the same +results when using the same seed. If you want randomization, set it to NULL, +otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed) +\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex} +and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods. +\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead. +}} + +\item{param_info}{Information on the parameters to estimate. +Either +a list containing: +\itemize{ +\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), +\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) +\item \code{init_values}, a data.frame containing initial +values to test for the parameters (optional, if not provided, or if less values +than number of repetitions of the minimization are provided, the, or part +of the, initial values will be randomly generated using LHS sampling within +parameter bounds). +} + +or a named list containing for each parameter: +\itemize{ +\item \code{sit_list}, list the groups of situations for which the current estimated +parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} +for an example), +\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), +\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). +\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) +}} + +\item{optim_results}{Results list returned by global optim method wrappers} + +\item{crit_options}{List containing several arguments given to \code{estim_param} +function: \code{param_names}, \code{obs_list}, \code{model_function}, +\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim} +that must be passed to main_crit function by the methods wrappers.} +} +\value{ +Updated results of global optim method +} +\description{ +Post-treat results of global optim methods +} diff --git a/man/summary_global_optim.Rd b/man/summary_global_optim.Rd new file mode 100644 index 00000000..6c3f65a8 --- /dev/null +++ b/man/summary_global_optim.Rd @@ -0,0 +1,52 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/global_optim_function.R +\name{summary_global_optim} +\alias{summary_global_optim} +\title{Summarizes results of global optim methods} +\usage{ +summary_global_optim(optim_options, param_info, optim_results, out_dir) +} +\arguments{ +\item{optim_options}{List of options of the parameter estimation method, containing: +\itemize{ +\item \code{ranseed} Set random seed so that each execution of estim_param give the same +results when using the same seed. If you want randomization, set it to NULL, +otherwise set it to a number of your choice (e.g. 1234) (optional, default to NULL, which means random seed) +\item specific options depending on the method used. Click on the links to see examples with the \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_simple_case.html}{simplex} +and \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_DREAM.html}{DreamZS} methods. +\item \code{out_dir} \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}} Definition of \code{out_dir} in \code{optim_options} is no longer supported, use the new argument \code{out_dir} of \code{estim_param} instead. +}} + +\item{param_info}{Information on the parameters to estimate. +Either +a list containing: +\itemize{ +\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), +\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) +\item \code{init_values}, a data.frame containing initial +values to test for the parameters (optional, if not provided, or if less values +than number of repetitions of the minimization are provided, the, or part +of the, initial values will be randomly generated using LHS sampling within +parameter bounds). +} + +or a named list containing for each parameter: +\itemize{ +\item \code{sit_list}, list the groups of situations for which the current estimated +parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} +for an example), +\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), +\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). +\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) +}} + +\item{optim_results}{Results list returned by global optim method wrappers} + +\item{out_dir}{Path to the directory where the optimization results will be written. (optional, default to \code{getwd()})} +} +\value{ +Prints results of global optim methods +} +\description{ +Summarizes results of global optim methods +} diff --git a/man/wrap_DEoptim.Rd b/man/wrap_DEoptim.Rd new file mode 100644 index 00000000..34e31468 --- /dev/null +++ b/man/wrap_DEoptim.Rd @@ -0,0 +1,33 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wrap_DE.R +\name{wrap_DEoptim} +\alias{wrap_DEoptim} +\title{A wrapper for DEoptim package} +\usage{ +wrap_DEoptim(optim_options, param_info, crit_options) +} +\arguments{ +\item{optim_options}{see description in estim_param} + +\item{param_info}{see description in estim_param} + +\item{crit_options}{List containing several arguments given to \code{estim_param} +function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function}, +\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir} +that must be passed to main_crit function by the methods wrappers.} +} +\value{ +prints, graphs and a list containing: +\code{final_values}, the vector of estimated values for optimized parameters +for the repetition that lead to the lowest value of the criterion +\code{init_values}, the vector of initial values for optimized parameters +\code{est_values}, data frame of values of all population +\code{ind_min_crit}, the index of the repetition that lead to the lowest value +of the criterion +\code{DE}, the data structure returned by DE +\code{crit_values} the vector of criterion values for the all population +} +\description{ +A wrapper for DEoptim package +} +\keyword{internal} diff --git a/man/wrap_graDiEnt.Rd b/man/wrap_graDiEnt.Rd new file mode 100644 index 00000000..b0e3224f --- /dev/null +++ b/man/wrap_graDiEnt.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wrap_graDiEnt.R +\name{wrap_graDiEnt} +\alias{wrap_graDiEnt} +\title{A wrapper for graDiEnt package} +\usage{ +wrap_graDiEnt(optim_options, param_info, crit_options) +} +\arguments{ +\item{optim_options}{see description in estim_param} + +\item{param_info}{see description in estim_param} + +\item{crit_options}{List containing several arguments given to \code{estim_param} +function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function}, +\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir} +that must be passed to main_crit function by the methods wrappers.} +} +\value{ +list with final_values, init_values, est_values, GE, min_crit_value +} +\description{ +A wrapper for graDiEnt package +} +\keyword{internal} diff --git a/tests/testthat/test_estim_param_DEoptim.R b/tests/testthat/test_estim_param_DEoptim.R new file mode 100644 index 00000000..4f62982d --- /dev/null +++ b/tests/testthat/test_estim_param_DEoptim.R @@ -0,0 +1,176 @@ +context("Test the estim_param function using a toy model") +library(tibble) + +# Define a toy model and its wrapper +toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) { + # Simulate biomass and yield with a simple logistic model + # + # Arguments + # - rB: growth rate + # - Bmin: minimum biomass + # - Bmax: maximum biomass + # - h: harvest rate + # - Bini: initial biomass + # + # Value + # A list containing the simulated biomass and yield + # + biom <- numeric(nb_days) + biom[1] <- Bini + for (i in 2:nb_days) { + biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax) + } + yield <- h * biom + return(list(biom = biom, yield = yield)) +} + +toymodel_wrapper <- function(param_values = NULL, situation, + model_options, ...) { + # A wrapper for toymodel + # + # Arguments + # - param_values: (optional) a named vector or a tibble containing the + # values of the toymodel parameters to force. + # - situation: Vector of situations names for which results must be + # returned. + # In this case, the names of the situations are coded as "year_suffix" + # - model_options: a tibble containing the begin and end date of simulation + # for each situation. + # - ...: mandatory since CroptimizR will give additional arguments not used + # here + # + # Value: + # A named list of tibble per situation. + # Each tibble contains columns: + # - Date (POSIXct dates of simulated results), + # - One column per simulated variable (biomass and yield) + # + results <- list( + sim_list = setNames( + vector("list", length(situation)), + nm = situation + ), + error = FALSE + ) + attr(results$sim_list, "class") <- "cropr_simulation" + + # Remove dummy in param_values, just here for specific tests + if (!is.null(param_values) && "dummy" %in% names(param_values)) { + param_values <- param_values[!names(param_values) %in% "dummy"] + } + + for (sit in situation) { + # Retrieve begin and end date from situation name + begin_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(begin_date) + end_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(end_date) + + date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day") + + if (!all(names(param_values) %in% c( + "rB", "Bmin", "Bmax", "h", "Bini" + ))) { + warning( + paste( + "Unknown parameters in param_values:", + paste(names(param_values), collapse = ",") + ) + ) + results$error <- TRUE + return(results) + } + + # Call the toymodel with varying arguments depending on what is given in + # param_values + res <- do.call( + "toymodel", + c(nb_days = length(date_sequence), as.list(param_values)) + ) + + # Fill the result variable + results$sim_list[[sit]] <- + dplyr::tibble( + Date = date_sequence, + biomass = res$biom, + yield = res$yield + ) + } + return(results) +} + +# To make these function accessible to the test environment +.GlobalEnv$toymodel <- toymodel +.GlobalEnv$toymodel_wrapper <- toymodel_wrapper + +# Use setup() to define shared data for all tests in this file +setup({ + # These variables will be available in all tests + .GlobalEnv$model_options <- tibble::tibble( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")), + end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18")) + ) + .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7) + + # Generate synthetic observations + tmp <- toymodel_wrapper( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + param_values = param_true_values, + model_options = model_options + ) + + .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) { + return(dplyr::filter(x, dplyr::row_number() %% 10 == 0)) + }) +}) + +# For cleaning up after tests +teardown({ + if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv) + if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv) + if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv) +}) + + + +# Test estim_param with DEoptim +test_that("estim_param 1 step OLS criterion", { + param_info <- list( + rB = list(lb = 0, ub = 1,init_values = 0.2), + h = list(lb = 0, ub = 1,init_values = 0.2), + Bmax = list(lb = 5, ub = 15,init_values = 6) + ) + + optim_options <- list( + ranseed = 1234, + itermax = 60, + NP = 30 + ) + + res <- estim_param( + obs_list = obs_synth, + model_function = toymodel_wrapper, + model_options = model_options, + optim_method = "DEoptim", + crit_function = crit_ols, + optim_options = optim_options, + param_info = param_info, + obs_var = c("biomass", "yield"), + situation = c("sit1_2000", "sit1_2001", "sit2_2003"), + out_dir = tempdir() + ) + + expect_equal(res$final_values[["rB"]], + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 + ) + expect_equal(res$final_values[["h"]], + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 + ) + expect_equal(res$final_values[["Bmax"]], + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 + ) +}) diff --git a/tests/testthat/test_estim_param_GEoptim.R b/tests/testthat/test_estim_param_GEoptim.R new file mode 100644 index 00000000..6b3361a3 --- /dev/null +++ b/tests/testthat/test_estim_param_GEoptim.R @@ -0,0 +1,179 @@ +context("Test the estim_param function using a toy model") +library(tibble) + +# Define a toy model and its wrapper +toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) { + # Simulate biomass and yield with a simple logistic model + # + # Arguments + # - rB: growth rate + # - Bmin: minimum biomass + # - Bmax: maximum biomass + # - h: harvest rate + # - Bini: initial biomass + # + # Value + # A list containing the simulated biomass and yield + # + biom <- numeric(nb_days) + biom[1] <- Bini + for (i in 2:nb_days) { + biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax) + } + yield <- h * biom + return(list(biom = biom, yield = yield)) +} + +toymodel_wrapper <- function(param_values = NULL, situation, + model_options, ...) { + # A wrapper for toymodel + # + # Arguments + # - param_values: (optional) a named vector or a tibble containing the + # values of the toymodel parameters to force. + # - situation: Vector of situations names for which results must be + # returned. + # In this case, the names of the situations are coded as "year_suffix" + # - model_options: a tibble containing the begin and end date of simulation + # for each situation. + # - ...: mandatory since CroptimizR will give additional arguments not used + # here + # + # Value: + # A named list of tibble per situation. + # Each tibble contains columns: + # - Date (POSIXct dates of simulated results), + # - One column per simulated variable (biomass and yield) + # + results <- list( + sim_list = setNames( + vector("list", length(situation)), + nm = situation + ), + error = FALSE + ) + attr(results$sim_list, "class") <- "cropr_simulation" + + # Remove dummy in param_values, just here for specific tests + if (!is.null(param_values) && "dummy" %in% names(param_values)) { + param_values <- param_values[!names(param_values) %in% "dummy"] + } + + for (sit in situation) { + # Retrieve begin and end date from situation name + begin_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(begin_date) + end_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(end_date) + + date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day") + + if (!all(names(param_values) %in% c( + "rB", "Bmin", "Bmax", "h", "Bini" + ))) { + warning( + paste( + "Unknown parameters in param_values:", + paste(names(param_values), collapse = ",") + ) + ) + results$error <- TRUE + return(results) + } + + # Call the toymodel with varying arguments depending on what is given in + # param_values + res <- do.call( + "toymodel", + c(nb_days = length(date_sequence), as.list(param_values)) + ) + + # Fill the result variable + results$sim_list[[sit]] <- + dplyr::tibble( + Date = date_sequence, + biomass = res$biom, + yield = res$yield + ) + } + return(results) +} + +# To make these function accessible to the test environment +.GlobalEnv$toymodel <- toymodel +.GlobalEnv$toymodel_wrapper <- toymodel_wrapper + +# Use setup() to define shared data for all tests in this file +setup({ + # These variables will be available in all tests + .GlobalEnv$model_options <- tibble::tibble( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")), + end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18")) + ) + .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7) + + # Generate synthetic observations + tmp <- toymodel_wrapper( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + param_values = param_true_values, + model_options = model_options + ) + + .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) { + return(dplyr::filter(x, dplyr::row_number() %% 10 == 0)) + }) +}) + +# For cleaning up after tests +teardown({ + if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv) + if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv) + if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv) +}) + + + +# Test estim_param with DEoptim +test_that("estim_param 1 step OLS criterion", { + param_info <- list( + rB = list(lb = 0, ub = 1), + h = list(lb = 0, ub = 1), + Bmax = list(lb = 5, ub = 15) + ) + + optim_options <- list( + ranseed = 1234, + n_params = length(param_info), + n_iter = 60, + n_particles = 30, + n_diff = 2, + return_trace = TRUE + ) + + res <- estim_param( + obs_list = obs_synth, + model_function = toymodel_wrapper, + model_options = model_options, + optim_method = "graDiEnt", + crit_function = crit_ols, + optim_options = optim_options, + param_info = param_info, + obs_var = c("biomass", "yield"), + situation = c("sit1_2000", "sit1_2001", "sit2_2003"), + out_dir = tempdir() + ) + + expect_equal(res$final_values[["rB"]], + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 + ) + expect_equal(res$final_values[["h"]], + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 + ) + expect_equal(res$final_values[["Bmax"]], + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 + ) +}) From cbe27c6254a1c9020c6d00f72918146f307c2610 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Mon, 24 Nov 2025 02:27:46 +0100 Subject: [PATCH 02/16] global_optim_function update update of global_optim_function in order to suit to global optim wraps --- R/global_optim_function.R | 381 ++++++++++------------ R/optim_switch.R | 6 + R/wrap_DE.R | 27 +- R/wrap_graDiEnt.R | 42 ++- tests/testthat/test_estim_param_DEoptim.R | 1 + 5 files changed, 240 insertions(+), 217 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index 8785169d..7eb9394a 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -9,37 +9,35 @@ summary_global_optim<- function(optim_options, param_info, optim_results, out_dir) { param_names <- get_params_names(param_info) nb_params <- length(param_names) - est_values <- optim_results$est_values - ind_min_crit <- optim_results$ind_min_crit + final_values <- optim_results$final_values min_crit_value <- optim_results$min_crit_value + est_values <- matrix(final_values, nrow = 1, dimnames = list(NULL, param_names)) cat(paste( "\nList of observed variables used:", paste(optim_results$obs_var_list, collapse = ", ") - )) - + )) # Display of parameters values for the candidate which has the # smallest criterion for (ipar in 1:nb_params) { cat(paste( "\nEstimated value for", param_names[ipar], ": ", - format(est_values[ind_min_crit, ipar], + format(est_values[1, ipar], scientific = FALSE, digits = 2, nsmall = 0 - ) + ) )) - } + } cat(paste( "\nMinimum value of the criterion:", format(min_crit_value, scientific = FALSE, digits = 2, nsmall = 0) - )) + )) cat(paste( "\nComplementary graphs and results can be found in ", out_dir, "\n" )) } - #' @title Post-treat results of global optim methods #' #' @inheritParams estim_param @@ -52,8 +50,9 @@ summary_global_optim<- function(optim_options, param_info, optim_results, out_di #' #' @return Updated results of global optim method #' + post_treat_global_optim <- function(optim_options, param_info, optim_results, - crit_options) { + crit_options) { param_names <- get_params_names(param_info) nb_params <- length(param_names) info_crit_list <- crit_options$info_crit_list @@ -63,27 +62,25 @@ post_treat_global_optim <- function(optim_options, param_info, optim_results, info_final <- main_crit( param_values = optim_results$final_values, crit_options = c(crit_options, return_detailed_info = TRUE) - ) + ) if (info_final$crit != optim_results$min_crit_value) { stop(paste( "Internal error: incoherent computation of minimum criterion value. \nValue obtained in method wrapper:", optim_results$min_crit_value, "\nValue obtained afterwards:", info_final$crit )) - } + } optim_results$forced_param_values <- info_final$forced_param_values - sapply(info_crit_list, function(x) { final_info_crit <- x( obs_list = info_final$obs_intersect, crit = info_final$crit, param_nb = nb_params - ) + ) optim_results[x()$name] <<- final_info_crit - }) - + }) return(optim_results) -} + } #' @title Generate plots for global optim methods @@ -96,145 +93,92 @@ post_treat_global_optim <- function(optim_options, param_info, optim_results, #' #' @keywords internal #' -plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) { - bounds <- get_params_bounds(param_info) - init_values <- optim_results$init_values +plot_global_optim <- function(optim_options, param_info, optim_results, out_dir){ est_values <- optim_results$est_values - #crit_values <- optim_results$crit_values + trace_df <- optim_results$trace_df + p_all <- list() # ValuesVSit plot - - tryCatch( - { - grDevices::pdf( - file = file.path(out_dir, "ValuesVSit_globaloptim2.pdf"), - width = 9, height = 9 - ) - }, - error = function(cond) { - filename <- paste0("ValuesVSit_new.pdf") - warning( - "Error trying to create ", out_dir, - "/ValuesVSit.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", - filename + if (!is.null(trace_df)) { + tryCatch( + { + grDevices::pdf( + file = file.path(out_dir, "ValuesVSit_globaloptim.pdf"), + width = 9, height = 9 + ) + }, + error = function(cond) { + filename <- paste0("ValuesVSit_new.pdf") + warning( + "Error trying to create ", out_dir, + "/ValuesVSit.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", + filename + ) + grDevices::pdf( + file = file.path(out_dir, filename), + width = 9, height = 9 + ) + } ) - grDevices::pdf( - file = file.path(out_dir, filename), - width = 9, height = 9 + # we should construct a dataframe df "individu × itération" + p <- plot_valuesVSit_go( + df = trace_df, + param_info = param_info, + iter_or_eval = "iter", + crit_log = FALSE, + ind_label = "begin" ) + for (plot in p) print(plot) + grDevices::dev.off() + p_all$valuesVSit <- p } - ) - - # DEoptim :we should construct a dataframe df "individu × itération" - if (!is.null(optim_results$DE)) { - - DE <- optim_results$DE - param_names <- get_params_names(param_info) - - # storepopfrom = 1 should be activated id wrap_DE - storepop <- DE$member$storepop - bestvalit <- DE$member$bestvalit - - if (length(storepop) > 0) { - itermax <- length(storepop) - NP <- nrow(storepop[[1]]) - - df_list <- vector("list", itermax) - - for (it in seq_len(itermax)) { - pop_it <- as.data.frame(storepop[[it]]) - colnames(pop_it) <- param_names - - pop_it$iter <- it # it id - pop_it$ind <- seq_len(NP) # id individual - # ------------------------------------------------ - ### To see with Samuel - # compute criterion per individual (To proprely plot the last figure of ValueVSit) (I have to add in parameters crit_options) - #crit_pop <- apply(storepop[[it]], 1, function(par) { - # main_crit(par, c(crit_options, return_detailed_info = FALSE))$crit - #}) - #pop_it$crit <- crit_pop - # ------------------------------------------------ - pop_it$crit <- bestvalit[it] - df_list[[it]] <- pop_it - } - - df_DE <- dplyr::bind_rows(df_list) - - # for plot_valuesVSit, we need eval col - df_DE$eval <- df_DE$iter - - p <- plot_valuesVSit_go( - df = df_DE, - param_info = param_info, - iter_or_eval = "iter", - crit_log = FALSE, - ind_label = "begin" - ) - - for (plot in p) print(plot) - grDevices::dev.off() - p_all$valuesVSit_DE <- p - } else { - warning("DE$member$storepop is empty : verify storepopfrom in DEoptim.control") - grDevices::dev.off() - } - } # ValuesVSit_2D plot if (!is.null(est_values)) { - - # Build a data.frame with final parameters - df_2d <- as.data.frame(est_values) - colnames(df_2d) <- get_params_names(param_info) - df_2d$eval <- seq_len(nrow(df_2d)) - df_2d$iter <- df_2d$eval - - tryCatch( - { - grDevices::pdf( - file = file.path(out_dir, "ValuesVSit_2D.pdf"), - width = 9, height = 9 - ) - }, - error = function(cond) { - filename <- paste0("ValuesVSit_2D_new.pdf") - warning( - "Error trying to create ", out_dir, - "/ValuesVSit_2D.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", - filename - ) - grDevices::pdf( - file = file.path(out_dir, filename), - width = 9, height = 9 - ) - } - ) - - - p <- plot_valuesVSit_2D_go( - df = df_2d, - param_info = param_info, - iter_or_eval = "eval", - fill = "eval", - crit_log = FALSE, - lines = FALSE, - rep_label = "end" - ) - - for (plot in p) { - print(plot) - } - grDevices::dev.off() - p_all$valuesVSit_2D <- p - } - - + # Build a data.frame with final parameters + df_2d <- as.data.frame(est_values) + colnames(df_2d) <- get_params_names(param_info) + df_2d$eval <- seq_len(nrow(df_2d)) + df_2d$iter <- df_2d$eval + tryCatch( + { + grDevices::pdf( + file = file.path(out_dir, "ValuesVSit_2D.pdf"), + width = 9, height = 9 + ) + }, + error = function(cond) { + filename <- paste0("ValuesVSit_2D_new.pdf") + warning( + "Error trying to create ", out_dir, + "/ValuesVSit_2D.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", + filename + ) + grDevices::pdf( + file = file.path(out_dir, filename), + width = 9, height = 9 + ) + } + ) + p <- plot_valuesVSit_2D_go( + df = df_2d, + param_info = param_info, + iter_or_eval = "eval", + fill = "eval", + crit_log = FALSE, + lines = FALSE, + rep_label = "end" + ) + for (plot in p) + print(plot) + grDevices::dev.off() + p_all$valuesVSit_2D <- p + } return(p_all) -} + } #' @title Create plots of parameters and criterion values per iteration or @@ -271,9 +215,10 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) #' #' @export #' + plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), - crit_log = TRUE, - ind_label = c("begin_end", "begin", "end")) { + crit_log = TRUE, + ind_label = c("begin_end", "begin", "end")) { param_names <- get_params_names(param_info) bounds <- get_params_bounds(param_info) @@ -282,18 +227,21 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), df <- filter(df, !is.na(.data$iter)) lab <- "iterations" } + have_crit <- "crit" %in% names(df) && !all(is.na(df$crit)) trans <- "identity" - mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) - if (crit_log) { + mid <- NA_real_ + if (crit_log && have_crit) { if (all(df$crit > 0)) { trans <- "log10" mid <- (max(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.") + 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 + } + } else if (!have_crit) { crit_log <- FALSE - } - } + } tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) tmp[tmp == Inf | tmp == -Inf] <- NA minvalue <- apply(tmp, 2, min, na.rm = TRUE) @@ -304,27 +252,31 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), p <- list() for (param_name in param_names) { - p[[param_name]] <- ggplot(df, aes( - x = !!rlang::sym(iter_or_eval[1]), - y = !!rlang::sym(param_name), - color = crit - )) + + p[[param_name]] <- ggplot( + df, aes(x = !!rlang::sym(iter_or_eval[1]), + y = !!rlang::sym(param_name), + color = crit + ) + ) + labs( title = paste0( "Evolution of ", param_name, " \n in function of the minimization ", lab - ), - y = param_name, - x = paste(lab, "number"), + ), + y = param_name, + x = paste(lab, "number"), fill = "Criterion" - ) + + ) + theme(plot.title = element_text(hjust = 0.5)) + - geom_point(alpha = 0.5) + - scale_color_gradient2( - midpoint = mid, low = "blue", mid = "yellow", - high = "red", space = "Lab", trans = trans - ) + geom_point(alpha = 0.5) + if (have_crit) { + p[[param_name]] <- p[[param_name]] + + scale_color_gradient2( + midpoint = mid, low = "blue", mid = "yellow", + high = "red", space = "Lab", trans = trans + ) + } for (iind in unique(df$ind)) { p[[param_name]] <- p[[param_name]] + geom_line(data = filter(df, ind == iind)) @@ -334,7 +286,7 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), size = 3 ) - } + } if (ind_label[1] == "begin_end" || ind_label[1] == "end") { p[[param_name]] <- p[[param_name]] + geom_label(aes(label = ind), @@ -342,20 +294,20 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), size = 3 ) } - } + } ylim(minvalue[param_name], maxvalue[param_name]) - } + } df$ind <- as.factor(df$ind) p[["criterion"]] <- ggplot(df, aes_string( x = iter_or_eval[1], y = "crit", color = "ind" - )) + + )) + labs( title = paste0( "Evolution of the minimized criterion \n in function of the minimization ", lab - ), + ), y = "Minimized criterion", x = paste(lab, "number"), fill = "Individual" @@ -373,14 +325,14 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), geom_label(aes(label = ind), data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)) ) - } + } if (crit_log) { p[["criterion"]] <- p[["criterion"]] + scale_y_log10() - } + } return(p) -} + } #' @title Create 2D plots of parameters values evolution per iteration or @@ -421,9 +373,9 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), #' @export #' plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "iter"), - fill = c("crit", "rep"), crit_log = TRUE, - lines = FALSE, - rep_label = c("begin_end", "begin", "end")) { + fill = c("crit", "rep"), crit_log = TRUE, + lines = FALSE, + rep_label = c("begin_end", "begin", "end")) { param_names <- get_params_names(param_info) if (length(param_names) <= 1) { return() @@ -436,28 +388,28 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite lab <- "iterations" } if ("rep" %in% names(df)) { - df$rep <- as.factor(df$rep) + df$rep <- as.factor(df$rep) + } + has_crit <- "crit" %in% names(df) + trans <- "identity" + mid <- NA_real_ + if (fill[1] == "crit" && has_crit) { + mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) + if (crit_log) { + if (all(df$crit > 0)) { + trans <- "log10" + mid <- (max(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 } - has_crit <- "crit" %in% names(df) - trans <- "identity" - mid <- NA_real_ - if (fill[1] == "crit" && has_crit) { - mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) - if (crit_log) { - if (all(df$crit > 0)) { - trans <- "log10" - mid <- (max(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 - } - } - } else { - crit_log <- FALSE - } + } + }else { + crit_log <- FALSE + } tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) - # -.data$ avoid NOTES on check ... + # check ... tmp[tmp == Inf | tmp == -Inf] <- NA minvalue <- apply(tmp, 2, min, na.rm = TRUE) maxvalue <- apply(tmp, 2, max, na.rm = TRUE) @@ -487,19 +439,18 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite geom_point(alpha = 0.5) if (fill[1] == "crit" && has_crit){ - p[[ipair]] <- p[[ipair]] + - scale_color_gradient2( - midpoint = mid, low = "blue", mid = "yellow", - high = "red", space = "Lab", trans = trans - ) - } else if (fill[1] == "eval") { - # palette continue simple pour eval - p[[ipair]] <- p[[ipair]] + - scale_color_gradient2( - low = "lightblue", - high = "darkblue" - ) - } + p[[ipair]] <- p[[ipair]] + + scale_color_gradient2( + midpoint = mid, low = "blue", mid = "yellow", + high = "red", space = "Lab", trans = trans + ) + } else if (fill[1] == "eval") { + p[[ipair]] <- p[[ipair]] + + scale_color_gradient2( + low = "lightblue", + high = "darkblue" + ) + } if (lines) { for (irep in unique(df$rep)) { @@ -508,15 +459,15 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite 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)), - size = 3 + 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)), - size = 3 + data = filter(df, rep == irep) %>% filter(eval == max(.data$eval)), + size = 3 ) } } @@ -526,4 +477,4 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite } return(p) -} + } diff --git a/R/optim_switch.R b/R/optim_switch.R index 1da98bce..1528262c 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -206,6 +206,12 @@ optim_switch <- function(...) { optim_results = res, out_dir = crit_options$out_dir ) + summary_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) }}else { flag_unknown_method <- TRUE diff --git a/R/wrap_DE.R b/R/wrap_DE.R index 975db799..22628889 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -88,12 +88,37 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { final_values <- DE$optim$bestmem names(final_values) <- param_names + trace_df <- NULL + if (!is.null(DE$member$storepop)) { + storepop <- DE$member$storepop + bestvalit <- DE$member$bestvalit + itermax <- length(storepop) + NP <- nrow(storepop[[1]]) + df_list <- vector("list", itermax) + eval_counter <- 0 + for (it in seq_len(itermax)) { + pop_it <- as.data.frame(storepop[[it]]) + colnames(pop_it) <- param_names + pop_it$ind <- seq_len(NP) # + pop_it$iter <- it + pop_it$crit <- bestvalit[it] + idx <- seq_len(NP) + pop_it$eval <- eval_counter + idx + eval_counter <- eval_counter + NP + pop_it$rep <- 1L + pop_it$method <- "DEoptim" + df_list[[it]] <- pop_it + } + trace_df <- dplyr::bind_rows(df_list) + } + res <- list( final_values = final_values, init_values = init_values, est_values = est_values, min_crit_value = DE$optim$bestval, - DE = DE + DE = DE, + trace_df = trace_df ) return(res) diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index def4231c..a7f23060 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -122,13 +122,53 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { ) } + trace_df <- NULL + if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) + { + tr <- SQGDE$particles_trace + n_it <- dim(tr)[1] + n_pop <- dim(tr)[2] + df_list <- vector("list", n_it) + eval_counter <- 0 + for (it in seq_len(n_it)) { + pop_u <- tr[it, , ] + pop_u <- as.matrix(pop_u) + colnames(pop_u) <- param_names + + # transformation + pop_x <- sweep(pop_u, 2, range_bounds, `*`) + pop_x <- sweep(pop_x, 2, bounds$lb, `+`) + + df <- as.data.frame(pop_x) + + df$ind <- seq_len(n_pop) + df$iter <- it + + crit_pop <- apply(pop_x, 1, function(par) { + names(par) <- param_names + main_crit(par, crit_options = crit_options) + }) + df$crit <- crit_pop + + idx <- seq_len(n_pop) + df$eval <- eval_counter + idx + eval_counter <- eval_counter + n_pop + + df$method <- "graDiEnt" + + df_list[[it]] <- df + } + trace_df <- dplyr::bind_rows(df_list) + } + res <- list( final_values = final_values, init_values = init_values, est_values = est_values, min_crit_value = min_crit_value, - SQGDE = SQGDE + SQGDE = SQGDE, + trace_df = trace_df ) return(res) diff --git a/tests/testthat/test_estim_param_DEoptim.R b/tests/testthat/test_estim_param_DEoptim.R index 4f62982d..7f752239 100644 --- a/tests/testthat/test_estim_param_DEoptim.R +++ b/tests/testthat/test_estim_param_DEoptim.R @@ -160,6 +160,7 @@ test_that("estim_param 1 step OLS criterion", { situation = c("sit1_2000", "sit1_2001", "sit2_2003"), out_dir = tempdir() ) + str(res$trace_df) expect_equal(res$final_values[["rB"]], param_true_values[["rB"]], From eedbb736724e65e6ee0e58973ead85977e07b6b8 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Mon, 24 Nov 2025 15:32:17 +0100 Subject: [PATCH 03/16] Update of global optim graphics --- R/global_optim_function.R | 98 +++++++++++++++++++++------------------ R/wrap_DE.R | 4 +- 2 files changed, 54 insertions(+), 48 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index 7eb9394a..a51e4158 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -99,83 +99,84 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) p_all <- list() - # ValuesVSit plot + # ValuesVSit + if (!is.null(trace_df)) { tryCatch( { grDevices::pdf( file = file.path(out_dir, "ValuesVSit_globaloptim.pdf"), width = 9, height = 9 - ) - }, + ) + }, error = function(cond) { filename <- paste0("ValuesVSit_new.pdf") warning( "Error trying to create ", out_dir, "/ValuesVSit.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", filename - ) + ) grDevices::pdf( file = file.path(out_dir, filename), width = 9, height = 9 ) - } - ) - # we should construct a dataframe df "individu × itération" + } + ) + + # trace_df : params, ind, iter, crit, eval, rep, method p <- plot_valuesVSit_go( df = trace_df, param_info = param_info, iter_or_eval = "iter", - crit_log = FALSE, + crit_log = FALSE, ind_label = "begin" - ) + ) for (plot in p) print(plot) grDevices::dev.off() p_all$valuesVSit <- p - } + } + - # ValuesVSit_2D plot + # ValuesVSit_2D + + if (!is.null(trace_df)) { + df_2d <- trace_df - if (!is.null(est_values)) { - # Build a data.frame with final parameters - df_2d <- as.data.frame(est_values) - colnames(df_2d) <- get_params_names(param_info) - df_2d$eval <- seq_len(nrow(df_2d)) - df_2d$iter <- df_2d$eval tryCatch( { grDevices::pdf( file = file.path(out_dir, "ValuesVSit_2D.pdf"), width = 9, height = 9 ) - }, + }, error = function(cond) { filename <- paste0("ValuesVSit_2D_new.pdf") warning( "Error trying to create ", out_dir, "/ValuesVSit_2D.pdf file. It is maybe opened in a pdf viewer and locked. It will be created under the name ", filename - ) + ) grDevices::pdf( file = file.path(out_dir, filename), width = 9, height = 9 ) - } - ) + } + ) + p <- plot_valuesVSit_2D_go( - df = df_2d, - param_info = param_info, + df = df_2d, + param_info = param_info, iter_or_eval = "eval", - fill = "eval", - crit_log = FALSE, - lines = FALSE, - rep_label = "end" - ) - for (plot in p) - print(plot) + fill = "eval", + crit_log = FALSE, + lines = FALSE, + rep_label = "end" + ) + + for (plot in p) print(plot) grDevices::dev.off() p_all$valuesVSit_2D <- p - } + } return(p_all) } @@ -230,18 +231,21 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), have_crit <- "crit" %in% names(df) && !all(is.na(df$crit)) trans <- "identity" mid <- NA_real_ - if (crit_log && have_crit) { - if (all(df$crit > 0)) { - trans <- "log10" - mid <- (max(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 - } - } else if (!have_crit) { - crit_log <- FALSE - } + if (have_crit) { + mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) + if (crit_log) { + if (all(df$crit > 0)) { + trans <- "log10" + mid <- (max(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 + } + } + } else { + crit_log <- FALSE + } tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) tmp[tmp == Inf | tmp == -Inf] <- NA minvalue <- apply(tmp, 2, min, na.rm = TRUE) @@ -294,7 +298,8 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), size = 3 ) } - } + } + p[[param_name]] <- p[[param_name]] + ylim(minvalue[param_name], maxvalue[param_name]) } @@ -472,8 +477,9 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite } } } - xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) - ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]]) + p[[ipair]] <- p[[ipair]] + + xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) + + ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]]) } return(p) diff --git a/R/wrap_DE.R b/R/wrap_DE.R index 22628889..94f9f0c3 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -99,13 +99,13 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { for (it in seq_len(itermax)) { pop_it <- as.data.frame(storepop[[it]]) colnames(pop_it) <- param_names - pop_it$ind <- seq_len(NP) # + pop_it$ind <- seq_len(NP) pop_it$iter <- it pop_it$crit <- bestvalit[it] idx <- seq_len(NP) pop_it$eval <- eval_counter + idx eval_counter <- eval_counter + NP - pop_it$rep <- 1L + pop_it$rep <- 1 pop_it$method <- "DEoptim" df_list[[it]] <- pop_it } From d35e6966757c003d83dd21ca3fbda8d2470ec0b9 Mon Sep 17 00:00:00 2001 From: sbuis Date: Mon, 24 Nov 2025 17:41:11 +0100 Subject: [PATCH 04/16] Fix: print of final estimated values of parameters --- R/global_optim_function.R | 166 +++++++++++++++++++------------------- 1 file changed, 83 insertions(+), 83 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index a51e4158..8ba6a259 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -6,32 +6,32 @@ #' #' @return Prints results of global optim methods #' -summary_global_optim<- function(optim_options, param_info, optim_results, out_dir) { +summary_global_optim <- function(optim_options, param_info, optim_results, out_dir) { param_names <- get_params_names(param_info) nb_params <- length(param_names) - final_values <- optim_results$final_values + final_values <- optim_results$final_values min_crit_value <- optim_results$min_crit_value est_values <- matrix(final_values, nrow = 1, dimnames = list(NULL, param_names)) cat(paste( "\nList of observed variables used:", paste(optim_results$obs_var_list, collapse = ", ") - )) + )) # Display of parameters values for the candidate which has the # smallest criterion for (ipar in 1:nb_params) { cat(paste( "\nEstimated value for", param_names[ipar], ": ", - format(est_values[1, ipar], - scientific = FALSE, - digits = 2, nsmall = 0 - ) + format(final_values[[ipar]], + scientific = FALSE, + digits = 2, nsmall = 0 + ) )) - } + } cat(paste( "\nMinimum value of the criterion:", format(min_crit_value, scientific = FALSE, digits = 2, nsmall = 0) - )) + )) cat(paste( "\nComplementary graphs and results can be found in ", out_dir, "\n" @@ -52,7 +52,7 @@ summary_global_optim<- function(optim_options, param_info, optim_results, out_di #' post_treat_global_optim <- function(optim_options, param_info, optim_results, - crit_options) { + crit_options) { param_names <- get_params_names(param_info) nb_params <- length(param_names) info_crit_list <- crit_options$info_crit_list @@ -62,25 +62,25 @@ post_treat_global_optim <- function(optim_options, param_info, optim_results, info_final <- main_crit( param_values = optim_results$final_values, crit_options = c(crit_options, return_detailed_info = TRUE) - ) + ) if (info_final$crit != optim_results$min_crit_value) { stop(paste( "Internal error: incoherent computation of minimum criterion value. \nValue obtained in method wrapper:", optim_results$min_crit_value, "\nValue obtained afterwards:", info_final$crit )) - } + } optim_results$forced_param_values <- info_final$forced_param_values sapply(info_crit_list, function(x) { final_info_crit <- x( obs_list = info_final$obs_intersect, crit = info_final$crit, param_nb = nb_params - ) + ) optim_results[x()$name] <<- final_info_crit - }) + }) return(optim_results) - } +} #' @title Generate plots for global optim methods @@ -93,10 +93,9 @@ post_treat_global_optim <- function(optim_options, param_info, optim_results, #' #' @keywords internal #' -plot_global_optim <- function(optim_options, param_info, optim_results, out_dir){ +plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) { est_values <- optim_results$est_values - trace_df <- optim_results$trace_df - + trace_df <- optim_results$trace_df p_all <- list() # ValuesVSit @@ -125,11 +124,11 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) # trace_df : params, ind, iter, crit, eval, rep, method p <- plot_valuesVSit_go( - df = trace_df, - param_info = param_info, + df = trace_df, + param_info = param_info, iter_or_eval = "iter", - crit_log = FALSE, - ind_label = "begin" + crit_log = FALSE, + ind_label = "begin" ) for (plot in p) print(plot) grDevices::dev.off() @@ -164,13 +163,13 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) ) p <- plot_valuesVSit_2D_go( - df = df_2d, - param_info = param_info, + df = df_2d, + param_info = param_info, iter_or_eval = "eval", - fill = "eval", - crit_log = FALSE, - lines = FALSE, - rep_label = "end" + fill = "eval", + crit_log = FALSE, + lines = FALSE, + rep_label = "end" ) for (plot in p) print(plot) @@ -179,7 +178,7 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) } return(p_all) - } +} #' @title Create plots of parameters and criterion values per iteration or @@ -233,19 +232,19 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), mid <- NA_real_ if (have_crit) { mid <- (max(df$crit) - min(df$crit)) / 2 + min(df$crit) - if (crit_log) { - if (all(df$crit > 0)) { - trans <- "log10" - mid <- (max(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 - } - } - } else { - crit_log <- FALSE - } + if (crit_log) { + if (all(df$crit > 0)) { + trans <- "log10" + mid <- (max(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 + } + } + } else { + crit_log <- FALSE + } tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) tmp[tmp == Inf | tmp == -Inf] <- NA minvalue <- apply(tmp, 2, min, na.rm = TRUE) @@ -257,20 +256,21 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), for (param_name in param_names) { p[[param_name]] <- ggplot( - df, aes(x = !!rlang::sym(iter_or_eval[1]), - y = !!rlang::sym(param_name), - color = crit - ) - ) + + df, aes( + x = !!rlang::sym(iter_or_eval[1]), + y = !!rlang::sym(param_name), + color = crit + ) + ) + labs( title = paste0( "Evolution of ", param_name, " \n in function of the minimization ", lab - ), - y = param_name, - x = paste(lab, "number"), + ), + y = param_name, + x = paste(lab, "number"), fill = "Criterion" - ) + + ) + theme(plot.title = element_text(hjust = 0.5)) + geom_point(alpha = 0.5) @@ -280,39 +280,39 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), midpoint = mid, low = "blue", mid = "yellow", high = "red", space = "Lab", trans = trans ) - } + } for (iind in unique(df$ind)) { p[[param_name]] <- p[[param_name]] + geom_line(data = filter(df, ind == iind)) if (ind_label[1] == "begin_end" || ind_label[1] == "begin") { p[[param_name]] <- p[[param_name]] + geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), - size = 3 + data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), + size = 3 ) - } + } if (ind_label[1] == "begin_end" || ind_label[1] == "end") { p[[param_name]] <- p[[param_name]] + geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)), - size = 3 + data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)), + size = 3 ) } } p[[param_name]] <- p[[param_name]] + - ylim(minvalue[param_name], maxvalue[param_name]) - } + ylim(minvalue[param_name], maxvalue[param_name]) + } df$ind <- as.factor(df$ind) p[["criterion"]] <- ggplot(df, aes_string( x = iter_or_eval[1], y = "crit", color = "ind" - )) + + )) + labs( title = paste0( "Evolution of the minimized criterion \n in function of the minimization ", lab - ), + ), y = "Minimized criterion", x = paste(lab, "number"), fill = "Individual" @@ -324,20 +324,20 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), p[["criterion"]] <- p[["criterion"]] + geom_line(data = filter(df, ind == iind)) + geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), - size = 3 + data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), + size = 3 ) + geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)) + data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)) ) - } + } if (crit_log) { p[["criterion"]] <- p[["criterion"]] + scale_y_log10() - } + } return(p) - } +} #' @title Create 2D plots of parameters values evolution per iteration or @@ -377,10 +377,10 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), #' #' @export #' -plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "iter"), - fill = c("crit", "rep"), crit_log = TRUE, - lines = FALSE, - rep_label = c("begin_end", "begin", "end")) { +plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "iter"), + fill = c("crit", "rep"), crit_log = TRUE, + lines = FALSE, + rep_label = c("begin_end", "begin", "end")) { param_names <- get_params_names(param_info) if (length(param_names) <= 1) { return() @@ -404,13 +404,13 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite 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 } - } - }else { + } + } else { crit_log <- FALSE } tmp <- rbind(bounds$lb, bounds$ub, select(df, all_of(param_names))) @@ -443,7 +443,7 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite theme(plot.title = element_text(hjust = 0.5)) + geom_point(alpha = 0.5) - if (fill[1] == "crit" && has_crit){ + if (fill[1] == "crit" && has_crit) { p[[ipair]] <- p[[ipair]] + scale_color_gradient2( midpoint = mid, low = "blue", mid = "yellow", @@ -455,7 +455,7 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite low = "lightblue", high = "darkblue" ) - } + } if (lines) { for (irep in unique(df$rep)) { @@ -464,23 +464,23 @@ plot_valuesVSit_2D_go <- function(df, param_info, iter_or_eval = c("eval", "ite 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)), - size = 3 + 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)), - size = 3 + data = filter(df, rep == irep) %>% filter(eval == max(.data$eval)), + size = 3 ) } } } p[[ipair]] <- p[[ipair]] + - xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) + - ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]]) + xlim(minvalue[df_pairs[1, ipair]], maxvalue[df_pairs[1, ipair]]) + + ylim(minvalue[df_pairs[2, ipair]], maxvalue[df_pairs[2, ipair]]) } return(p) - } +} From cf6ce30089fd3190ac2075d528d3be96f3796b31 Mon Sep 17 00:00:00 2001 From: sbuis Date: Mon, 24 Nov 2025 17:42:47 +0100 Subject: [PATCH 05/16] format: apply styler formatting --- R/optim_switch.R | 84 +++++++++++++++++++++++------------------------ R/wrap_graDiEnt.R | 39 ++++++++++++---------- 2 files changed, 63 insertions(+), 60 deletions(-) diff --git a/R/optim_switch.R b/R/optim_switch.R index 1528262c..3438dc46 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -72,7 +72,7 @@ optim_switch <- function(...) { } } - #if (!is.null(crit_options$out_dir) & length(res) > 0) { + # if (!is.null(crit_options$out_dir) & length(res) > 0) { if (!is.null(crit_options) && !is.null(crit_options$out_dir) && length(res) > 0) { save(res, file = file.path( crit_options$out_dir, @@ -99,7 +99,7 @@ optim_switch <- function(...) { ########################## # # param_info <- NULL - #crit_options <- NULL + # crit_options <- NULL # if ("param_info" %in% names(arguments)) { # param_info <- arguments$param_info @@ -130,7 +130,7 @@ optim_switch <- function(...) { optim_results = res, crit_options = crit_options ) - res$plots <- plot_frequentist( + res$plots <- plot_frequentist( optim_options = optim_options, param_info = param_info, optim_results = res, @@ -168,25 +168,25 @@ optim_switch <- function(...) { if (arguments$crit_options$info_level >= 1) { res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) } - res <- post_treat_global_optim( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - crit_options = crit_options - ) - res$plots <- plot_global_optim( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - out_dir = crit_options$out_dir - ) - summary_global_optim( - optim_options = optim_options, param_info = param_info, - optim_results = res, - out_dir = crit_options$out_dir - ) -} - } else if (optim_method == "graDiEnt") { + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + summary_global_optim( + optim_options = optim_options, param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + } + } else if (optim_method == "graDiEnt") { res <- do.call(wrap_graDiEnt, wrap_args) if (nargs() > 2) { res$obs_var_list <- .croptEnv$obs_var_list @@ -194,26 +194,26 @@ optim_switch <- function(...) { if (arguments$crit_options$info_level >= 1) { res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) } - res <- post_treat_global_optim( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - crit_options = crit_options - ) - res$plots <- plot_global_optim( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - out_dir = crit_options$out_dir - ) - summary_global_optim( - optim_options = optim_options, - param_info = param_info, - optim_results = res, - out_dir = crit_options$out_dir - ) - - }}else { + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + summary_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + } + } else { flag_unknown_method <- TRUE }, error = function(cond) { diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index a7f23060..e9a30724 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -6,8 +6,8 @@ #' @importFrom graDiEnt optim_SQGDE GetAlgoParams #' @keywords internal wrap_graDiEnt <- function(optim_options, param_info, crit_options) { - if(is.null((ranseed <- optim_options$ranseed))) { - ranseed = NULL + if (is.null((ranseed <- optim_options$ranseed))) { + ranseed <- NULL } if (is.null((reltol <- optim_options$reltol))) { optim_options$converge_crit <- "stdev" @@ -20,7 +20,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { if (is.null((trace <- optim_options$return_trace))) { optim_options$return_trace <- TRUE } - optim_options$ranseed = NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it + optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it algorithm <- "graDiEnt" # return requested information if only optim_options is given in argument @@ -40,7 +40,8 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { bounds <- get_params_bounds(param_info) init_values <- get_init_values(param_info) - control_params <- do.call(graDiEnt::GetAlgoParams,optim_options) + control_params <- do.call(graDiEnt::GetAlgoParams, optim_options) + # control_params$init_center <- as.vector(init_values[1, ]) crit_options$tot_max_eval <- control_params$n_particles * control_params$n_iter @@ -57,15 +58,18 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { } SQGDE <- tryCatch( graDiEnt::optim_SQGDE( - ObjFun = ObjFun, + ObjFun = ObjFun, control_params = control_params ), - error = function(e) { warning(sprintf("GraDiEnt failed: %s", e$message)); NULL } + error = function(e) { + warning(sprintf("GraDiEnt failed: %s", e$message)) + NULL + } ) elapsed <- Sys.time() - start_time # Verify criterion value - if (is.null(SQGDE) || is.null(SQGDE$solution)) { + if (is.null(SQGDE) || is.null(SQGDE$solution)) { stop( "The value of the criterion is NA." ) @@ -76,18 +80,18 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { final_values <- bounds$lb + u_sol * range_bounds names(final_values) <- param_names min_crit_value <- SQGDE$solution_weight - if (is.null(min_crit_value) ) { - min_crit_value <- main_crit(final_values, crit_options) + if (is.null(min_crit_value)) { + min_crit_value <- main_crit(final_values, crit_options) } # Get the estimated values ==> matrix of 1 * nb_params - est_values <- matrix( - final_values, - nrow = 1L, - dimnames = list(NULL, param_names) - ) + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) colnames(est_values) <- param_names if (ncol(est_values) == nb_params) { colnames(est_values) <- param_names @@ -104,9 +108,9 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { dimnames = list(NULL, param_names) ) } else { - n_it <- dim(tr)[1] - last_pop <- tr[n_it, , ] # matrice n_particles × n_params (en u) - est_u <- as.matrix(last_pop) + n_it <- dim(tr)[1] + last_pop <- tr[n_it, , ] # matrice n_particles × n_params (en u) + est_u <- as.matrix(last_pop) colnames(est_u) <- param_names # transformation vers l’échelle physique @@ -173,4 +177,3 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { return(res) } - From aed85c9505f37ba972781454fac5d009c108f792 Mon Sep 17 00:00:00 2001 From: sbuis Date: Mon, 24 Nov 2025 16:53:25 +0000 Subject: [PATCH 06/16] Style code (GHA) --- R/FwdRegAgMIP.R | 4 --- R/bayesian_functions.R | 1 - R/information_criteria.R | 2 -- R/param_info_functions.R | 2 -- R/wrap_BayesianTools.R | 1 - R/wrap_DE.R | 31 +++++++++++---------- R/wrap_graDiEnt.R | 14 ++++------ tests/testthat/test-get_params.R | 3 -- tests/testthat/test-init_values_functions.R | 2 -- tests/testthat/test_estim_param_DEoptim.R | 19 ++++++------- tests/testthat/test_estim_param_GEoptim.R | 20 ++++++------- vignettes/Designing_a_model_wrapper.Rmd | 4 +-- 12 files changed, 43 insertions(+), 60 deletions(-) diff --git a/R/FwdRegAgMIP.R b/R/FwdRegAgMIP.R index 0aa3e8fc..4ecc5a10 100644 --- a/R/FwdRegAgMIP.R +++ b/R/FwdRegAgMIP.R @@ -96,8 +96,6 @@ select_param_FwdRegAgMIP <- function(oblig_param_list, add_param_list, crt_list, } - - #' @title Post-treat the results of the Forward Selection algorithm proposed in #' AgMIP calibration phaseIII protocol #' @@ -198,8 +196,6 @@ summary_FwdRegAgMIP <- function(param_selection_steps, } - - #' @title Save the results of the Forward Selection algorithm proposed in #' AgMIP calibration phaseIII protocol #' diff --git a/R/bayesian_functions.R b/R/bayesian_functions.R index f169053a..c7f8cc43 100644 --- a/R/bayesian_functions.R +++ b/R/bayesian_functions.R @@ -32,7 +32,6 @@ summary_bayesian <- function(optim_options, param_info, optim_results, out_dir) } - #' @title Generate plots for bayesian methods #' #' @inheritParams estim_param diff --git a/R/information_criteria.R b/R/information_criteria.R index e6c8346a..0aaf4af4 100644 --- a/R/information_criteria.R +++ b/R/information_criteria.R @@ -25,7 +25,6 @@ AIC <- function(obs_list, crit_value, param_nb) { } - #' @title Computes AICc for ordinary least squares #' #' @inheritParams estim_param @@ -62,7 +61,6 @@ AICc <- function(obs_list, crit_value, param_nb) { } - #' @title Computes BIC for ordinary least squares #' #' @inheritParams estim_param diff --git a/R/param_info_functions.R b/R/param_info_functions.R index 4e80b551..95566c0b 100644 --- a/R/param_info_functions.R +++ b/R/param_info_functions.R @@ -71,7 +71,6 @@ get_params_names <- function(param_info, short_list = FALSE) { } - #' @title Extract bounds from parameter information #' #' @inheritParams estim_param @@ -418,7 +417,6 @@ get_init_values <- function(param_info) { } - #' @title Complete initial values in param_info by sampling values from bounds taking into account inequality constraints between parameters #' #' @inheritParams estim_param diff --git a/R/wrap_BayesianTools.R b/R/wrap_BayesianTools.R index 6045182b..a84ade2c 100644 --- a/R/wrap_BayesianTools.R +++ b/R/wrap_BayesianTools.R @@ -37,7 +37,6 @@ wrap_BayesianTools <- function(optim_options, param_info, crit_options) { } - crit_options$tot_max_eval <- optim_options$iterations + optim_options$startValue - optim_options$iterations %% optim_options$startValue diff --git a/R/wrap_DE.R b/R/wrap_DE.R index 94f9f0c3..9ec5455d 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -15,9 +15,9 @@ #' @keywords internal wrap_DEoptim <- function(optim_options, param_info, crit_options) { - if(is.null((ranseed <- optim_options$ranseed))) { - ranseed = NULL - } + if (is.null((ranseed <- optim_options$ranseed))) { + ranseed <- NULL + } if (is.null((reltol <- optim_options$reltol))) { optim_options$reltol <- 1e-10 } @@ -29,7 +29,7 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { if (is.null((trace <- optim_options$trace))) { optim_options$trace <- FALSE } - optim_options$ranseed = NULL # ranseed is set to NULL because DEoptim.controll doesn't regonise it + optim_options$ranseed <- NULL # ranseed is set to NULL because DEoptim.controll doesn't regonise it optim_options$storepopfrom <- 1 optim_options$storepopfreq <- 1 control <- do.call(DEoptim::DEoptim.control, optim_options) @@ -59,17 +59,20 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { } DE <- tryCatch( DEoptim::DEoptim( - fn = fn_de, + fn = fn_de, lower = bounds$lb, upper = bounds$ub, control = control ), - error = function(e) { warning(sprintf("DEoptim failed: %s", e$message)); NULL } + error = function(e) { + warning(sprintf("DEoptim failed: %s", e$message)) + NULL + } ) elapsed <- Sys.time() - start_time # Verify criterion value - if (is.na( DE$optim$bestval)) { + if (is.na(DE$optim$bestval)) { stop( "The value of the criterion is NA." ) @@ -90,27 +93,27 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { trace_df <- NULL if (!is.null(DE$member$storepop)) { - storepop <- DE$member$storepop + storepop <- DE$member$storepop bestvalit <- DE$member$bestvalit itermax <- length(storepop) - NP <- nrow(storepop[[1]]) - df_list <- vector("list", itermax) + NP <- nrow(storepop[[1]]) + df_list <- vector("list", itermax) eval_counter <- 0 for (it in seq_len(itermax)) { pop_it <- as.data.frame(storepop[[it]]) colnames(pop_it) <- param_names - pop_it$ind <- seq_len(NP) + pop_it$ind <- seq_len(NP) pop_it$iter <- it pop_it$crit <- bestvalit[it] idx <- seq_len(NP) pop_it$eval <- eval_counter + idx eval_counter <- eval_counter + NP - pop_it$rep <- 1 + pop_it$rep <- 1 pop_it$method <- "DEoptim" df_list[[it]] <- pop_it - } - trace_df <- dplyr::bind_rows(df_list) } + trace_df <- dplyr::bind_rows(df_list) + } res <- list( final_values = final_values, diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index e9a30724..0cbf6eb8 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -85,7 +85,6 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { } - # Get the estimated values ==> matrix of 1 * nb_params est_values <- matrix( final_values, @@ -127,10 +126,9 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { } trace_df <- NULL - if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) - { + if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) { tr <- SQGDE$particles_trace - n_it <- dim(tr)[1] + n_it <- dim(tr)[1] n_pop <- dim(tr)[2] df_list <- vector("list", n_it) eval_counter <- 0 @@ -145,13 +143,13 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df <- as.data.frame(pop_x) - df$ind <- seq_len(n_pop) + df$ind <- seq_len(n_pop) df$iter <- it crit_pop <- apply(pop_x, 1, function(par) { names(par) <- param_names main_crit(par, crit_options = crit_options) - }) + }) df$crit <- crit_pop idx <- seq_len(n_pop) @@ -161,9 +159,9 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df$method <- "graDiEnt" df_list[[it]] <- df - } - trace_df <- dplyr::bind_rows(df_list) } + trace_df <- dplyr::bind_rows(df_list) + } res <- list( diff --git a/tests/testthat/test-get_params.R b/tests/testthat/test-get_params.R index bd4f6b04..a3f1d739 100644 --- a/tests/testthat/test-get_params.R +++ b/tests/testthat/test-get_params.R @@ -39,7 +39,6 @@ test_that("get_params_bounds", { }) - # Test get_init_values prior_1 <- list( init_values = c(dlaimax = 0.001, durvieF = 200), @@ -112,7 +111,6 @@ test_that("get_init_values", { }) - # Test get_params_names prior_1 <- list( lb = c(dlaimax = 0.0005, durvieF = 50), @@ -176,7 +174,6 @@ test_that("get_params_names", { }) - # Test get_params_per_sit sg <- list( p1 = list(sit_list = list( diff --git a/tests/testthat/test-init_values_functions.R b/tests/testthat/test-init_values_functions.R index 004a337c..7ea35be7 100644 --- a/tests/testthat/test-init_values_functions.R +++ b/tests/testthat/test-init_values_functions.R @@ -124,8 +124,6 @@ test_that("set_init_values: Case 1 with groups of situations per parameter", { }) - - # Case 2 with groups of situations per parameter param_info <- list() diff --git a/tests/testthat/test_estim_param_DEoptim.R b/tests/testthat/test_estim_param_DEoptim.R index 7f752239..32e13ce9 100644 --- a/tests/testthat/test_estim_param_DEoptim.R +++ b/tests/testthat/test_estim_param_DEoptim.R @@ -133,13 +133,12 @@ teardown({ }) - # Test estim_param with DEoptim test_that("estim_param 1 step OLS criterion", { param_info <- list( - rB = list(lb = 0, ub = 1,init_values = 0.2), - h = list(lb = 0, ub = 1,init_values = 0.2), - Bmax = list(lb = 5, ub = 15,init_values = 6) + rB = list(lb = 0, ub = 1, init_values = 0.2), + h = list(lb = 0, ub = 1, init_values = 0.2), + Bmax = list(lb = 5, ub = 15, init_values = 6) ) optim_options <- list( @@ -163,15 +162,15 @@ test_that("estim_param 1 step OLS criterion", { str(res$trace_df) expect_equal(res$final_values[["rB"]], - param_true_values[["rB"]], - tolerance = param_true_values[["rB"]] * 1e-2 + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 ) expect_equal(res$final_values[["h"]], - param_true_values[["h"]], - tolerance = param_true_values[["h"]] * 1e-2 + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 ) expect_equal(res$final_values[["Bmax"]], - param_true_values[["Bmax"]], - tolerance = param_true_values[["Bmax"]] * 1e-2 + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 ) }) diff --git a/tests/testthat/test_estim_param_GEoptim.R b/tests/testthat/test_estim_param_GEoptim.R index 6b3361a3..09c88a80 100644 --- a/tests/testthat/test_estim_param_GEoptim.R +++ b/tests/testthat/test_estim_param_GEoptim.R @@ -144,10 +144,10 @@ test_that("estim_param 1 step OLS criterion", { optim_options <- list( ranseed = 1234, - n_params = length(param_info), - n_iter = 60, - n_particles = 30, - n_diff = 2, + n_params = length(param_info), + n_iter = 60, + n_particles = 30, + n_diff = 2, return_trace = TRUE ) @@ -165,15 +165,15 @@ test_that("estim_param 1 step OLS criterion", { ) expect_equal(res$final_values[["rB"]], - param_true_values[["rB"]], - tolerance = param_true_values[["rB"]] * 1e-2 + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 ) expect_equal(res$final_values[["h"]], - param_true_values[["h"]], - tolerance = param_true_values[["h"]] * 1e-2 + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 ) expect_equal(res$final_values[["Bmax"]], - param_true_values[["Bmax"]], - tolerance = param_true_values[["Bmax"]] * 1e-2 + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 ) }) diff --git a/vignettes/Designing_a_model_wrapper.Rmd b/vignettes/Designing_a_model_wrapper.Rmd index 859ddc26..89f0557e 100644 --- a/vignettes/Designing_a_model_wrapper.Rmd +++ b/vignettes/Designing_a_model_wrapper.Rmd @@ -67,9 +67,7 @@ Here's a header for a basic version of a model wrapper: #' o `error`: an error code indicating if at least one simulation ended #' with an error. -model_wrapper <- function(param_values = NULL, situation, model_options, ...) { - -} +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. From aaa33b980023a1c880f27e14bd944538dc02b717 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Fri, 5 Dec 2025 13:02:28 +0100 Subject: [PATCH 07/16] Code cleanup in wrap_graDiEnt --- R/wrap_graDiEnt.R | 88 +++++++++++++++++++++-------------------------- 1 file changed, 39 insertions(+), 49 deletions(-) diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index 0cbf6eb8..03646ff2 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -5,44 +5,40 @@ #' @return list with final_values, init_values, est_values, GE, min_crit_value #' @importFrom graDiEnt optim_SQGDE GetAlgoParams #' @keywords internal +#' wrap_graDiEnt <- function(optim_options, param_info, crit_options) { - if (is.null((ranseed <- optim_options$ranseed))) { - ranseed <- NULL - } - if (is.null((reltol <- optim_options$reltol))) { - optim_options$converge_crit <- "stdev" - } + if(is.null((ranseed <- optim_options$ranseed))) { + ranseed = NULL + } if (is.null((n_params <- optim_options$n_params))) { stop( "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)" ) - } + } if (is.null((trace <- optim_options$return_trace))) { optim_options$return_trace <- TRUE - } + } optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it algorithm <- "graDiEnt" # return requested information if only optim_options is given in argument if (nargs() == 1 && methods::hasArg(optim_options)) { - init_nb <- if (!is.null(optim_options$n_particles)) { - optim_options$n_particles - } else { - NA_integer_ - } - return(list( + init_nb <- if (!is.null(optim_options$n_particles)) { + optim_options$n_particles + } else { + NA_integer_ + } + return(list( package = "graDiEnt", family = "Global", method = "SQGDE", init_values_nb = init_nb - )) - } + )) + } param_names <- get_params_names(param_info) nb_params <- length(param_names) bounds <- get_params_bounds(param_info) init_values <- get_init_values(param_info) - control_params <- do.call(graDiEnt::GetAlgoParams, optim_options) - # control_params$init_center <- as.vector(init_values[1, ]) - + control_params <- do.call(graDiEnt::GetAlgoParams,optim_options) crit_options$tot_max_eval <- control_params$n_particles * control_params$n_iter start_time <- Sys.time() @@ -58,13 +54,12 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { } SQGDE <- tryCatch( graDiEnt::optim_SQGDE( - ObjFun = ObjFun, + ObjFun = ObjFun, control_params = control_params ), error = function(e) { warning(sprintf("GraDiEnt failed: %s", e$message)) - NULL - } + NULL } ) elapsed <- Sys.time() - start_time @@ -74,27 +69,26 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { "The value of the criterion is NA." ) } + u_sol <- SQGDE$solution names(u_sol) <- param_names final_values <- bounds$lb + u_sol * range_bounds names(final_values) <- param_names min_crit_value <- SQGDE$solution_weight - if (is.null(min_crit_value)) { + if (is.null(min_crit_value) ) { min_crit_value <- main_crit(final_values, crit_options) } - - # Get the estimated values ==> matrix of 1 * nb_params est_values <- matrix( final_values, nrow = 1L, dimnames = list(NULL, param_names) - ) + ) colnames(est_values) <- param_names if (ncol(est_values) == nb_params) { colnames(est_values) <- param_names - } + } # Final solution if (!is.null(SQGDE$particles_trace)) { @@ -107,28 +101,28 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { dimnames = list(NULL, param_names) ) } else { - n_it <- dim(tr)[1] - last_pop <- tr[n_it, , ] # matrice n_particles × n_params (en u) - est_u <- as.matrix(last_pop) + n_it <- dim(tr)[1] + last_pop <- tr[n_it, , ] + est_u <- as.matrix(last_pop) colnames(est_u) <- param_names - # transformation vers l’échelle physique + # transformation est_values <- sweep(est_u, 2, range_bounds, `*`) est_values <- sweep(est_values, 2, bounds$lb, `+`) colnames(est_values) <- param_names - } - } else { - est_values <- matrix( + } + } else { + est_values <- matrix( final_values, nrow = 1L, dimnames = list(NULL, param_names) - ) - } - + ) + } trace_df <- NULL - if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) { + if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) + { tr <- SQGDE$particles_trace - n_it <- dim(tr)[1] + n_it <- dim(tr)[1] n_pop <- dim(tr)[2] df_list <- vector("list", n_it) eval_counter <- 0 @@ -143,27 +137,22 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df <- as.data.frame(pop_x) - df$ind <- seq_len(n_pop) + df$ind <- seq_len(n_pop) df$iter <- it - crit_pop <- apply(pop_x, 1, function(par) { names(par) <- param_names main_crit(par, crit_options = crit_options) - }) + } + ) df$crit <- crit_pop - idx <- seq_len(n_pop) df$eval <- eval_counter + idx eval_counter <- eval_counter + n_pop - df$method <- "graDiEnt" - df_list[[it]] <- df - } + } trace_df <- dplyr::bind_rows(df_list) - } - - + } res <- list( final_values = final_values, init_values = init_values, @@ -175,3 +164,4 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { return(res) } + From 48822f3e87317a11fc320925decce8da4da231f3 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Sun, 14 Dec 2025 20:45:29 +0100 Subject: [PATCH 08/16] parma::cmaes package --- NAMESPACE | 2 + R/optim_switch.R | 58 ++++++- R/wrap_cmaes.R | 215 ++++++++++++++++++++++++ man/wrap_cmaes.Rd | 25 +++ tests/testthat/test_estim_param_cmaes.R | 180 ++++++++++++++++++++ 5 files changed, 474 insertions(+), 6 deletions(-) create mode 100644 R/wrap_cmaes.R create mode 100644 man/wrap_cmaes.Rd create mode 100644 tests/testthat/test_estim_param_cmaes.R diff --git a/NAMESPACE b/NAMESPACE index fec65822..c5ef7162 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -50,6 +50,8 @@ importFrom(ggplot2,xlim) importFrom(ggplot2,ylim) importFrom(graDiEnt,GetAlgoParams) importFrom(graDiEnt,optim_SQGDE) +importFrom(parma,cmaes) +importFrom(parma,cmaes.control) importFrom(purrr,modify) importFrom(rlang,.data) importFrom(stats,setNames) diff --git a/R/optim_switch.R b/R/optim_switch.R index 3438dc46..58b21b2a 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -16,11 +16,13 @@ #' optim_switch <- function(...) { + message("optim_switch nargs = ", nargs()) + message("names(arguments) = ", paste(names(list(...)), collapse = ", ")) + # Store and save results at the end of the current optimization step + res <- list() + flag_error <- FALSE if (nargs() > 2) { - # Initialize res - res <- list() - flag_error <- FALSE on.exit({ #### @@ -213,13 +215,57 @@ optim_switch <- function(...) { out_dir = crit_options$out_dir ) } - } else { + } else if (optim_method == "cmaes") { + + + if (is.null(optim_options$n_particles) || is.na(optim_options$n_particles)) { + if (!is.null(optim_options$ctrl$options$PopSize) && + !is.na(optim_options$ctrl$options$PopSize)) { + optim_options$n_particles <- as.integer(optim_options$ctrl$options$PopSize) + } else { + optim_options$n_particles <- 30L + } + } + + wrap_args$optim_options <- optim_options + + res <- do.call(wrap_cmaes, wrap_args) + + if (nargs() > 2) { + res$obs_var_list <- .croptEnv$obs_var_list + res$obs_situation_list <- .croptEnv$obs_situation_list + if (arguments$crit_options$info_level >= 1) { + res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) + } + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + summary_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + }} + +else { flag_unknown_method <- TRUE }, error = function(cond) { - warning(cond) - flag_error <<- TRUE + message("=== ERROR CAUGHT IN optim_switch ===") + traceback(20) + stop(cond) } + ) if (flag_unknown_method) { diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R new file mode 100644 index 00000000..67a74e44 --- /dev/null +++ b/R/wrap_cmaes.R @@ -0,0 +1,215 @@ +#' @title A wrapper for parma::cmaes package (CMA-ES) +#' +#' @inheritParams optim_switch +#' +#' @return list with final_values, init_values, est_values, GE, min_crit_value +#' @importFrom parma cmaes cmaes.control +#' @keywords internal +#' +wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { + if(is.null((ranseed <- optim_options$ranseed))) { + ranseed = NULL + } + n_params <- optim_options$n_params + if (is.null(n_params)) { + n_params <- length(param_info) + } + if (is.null((trace <- optim_options$return_trace))) { + optim_options$return_trace <- TRUE + } + optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it + algorithm <- "parma_cmaes" + + if (is.null(param_info) || is.null(crit_options)) { + init_nb <- if (!is.null(optim_options$n_particles)) optim_options$n_particles else NA_integer_ + return(list( + package = "parma", family = "Global", + method = "cmaes", init_values_nb = init_nb + )) + } + param_names <- get_params_names(param_info, short_list = TRUE) + nb_params <- length(param_names) + bounds <- get_params_bounds(param_info) + init_values <- get_init_values(param_info) + + control_params <- optim_options$ctrl + + # base ctrl (toujours valide) + ctrl <- parma::cmaes.control() + + # si l'utilisateur a déjà fourni un ctrl complet (déjà construit) + if (!is.null(control_params) && !is.null(control_params$options) && !is.null(control_params$CMA) && + !is.list(control_params$options)) { + # (cas rare) on ignore, car options devrait être une liste + } + + # si l'utilisateur a fourni une liste ctrl = list(options=..., CMA=...) + if (!is.null(control_params) && is.list(control_params) && + (!is.null(control_params$options) || !is.null(control_params$CMA))) { + + if (!is.null(control_params$options)) { + for (nm in names(control_params$options)) { + ctrl$options[[nm]] <- control_params$options[[nm]] + } + } + + if (!is.null(control_params$CMA)) { + for (nm in names(control_params$CMA)) { + ctrl$CMA[[nm]] <- control_params$CMA[[nm]] + } + } + + # si l'utilisateur a fourni directement une liste d'options (sans options/CMA) + } else if (!is.null(control_params) && is.list(control_params)) { + for (nm in names(control_params)) { + ctrl$options[[nm]] <- control_params[[nm]] + } + } + + if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){ + crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter) + } + + start_time <- Sys.time() + + if (!is.null(ranseed)) set.seed(ranseed) + trace_env <- new.env(parent = emptyenv()) + trace_env$eval <- 0L + trace_env$crit <- numeric(0) + trace_env$X <- list() + keep_trace <- isTRUE(optim_options$return_trace) + ObjFun <- function(x) { + x <- as.numeric(x) + names(x) <- param_names + val <- main_crit(x, crit_options = crit_options) + + if (keep_trace) { + trace_env$eval <- trace_env$eval + 1L + trace_env$crit[trace_env$eval] <- val + trace_env$X[[trace_env$eval]] <- x + } + return(val) + } + insigma <- optim_options$insigma + if (is.null(insigma)) { + # default value if it is not provided + insigma <- rep(0.3, nb_params) + } + if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) + + # 🟩 parma::cmaes expects a numeric vector for 'pars' (length = nb_params) + x0 <- init_values + + # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) + if (is.data.frame(x0)) x0 <- as.matrix(x0) + + if (is.matrix(x0)) { + # take the first row as starting point (you can also use colMeans) + x0 <- x0[1, , drop = TRUE] + } + + x0 <- as.numeric(x0) + names(x0) <- param_names + + # safety: if still wrong length, fallback to mid-bounds + if (length(x0) != nb_params) { + x0 <- (bounds$lb + bounds$ub) / 2 + names(x0) <- param_names + } + + parma_res <- tryCatch( + parma::cmaes( + pars = x0, + fun = ObjFun, + lower = bounds$lb, + upper = bounds$ub, + insigma = insigma, + ctrl = ctrl + ), + error = function(e) { + warning(sprintf("parma::cmaes failed: %s", e$message)) + NULL + } + ) + if (isTRUE(optim_options$debug)) { + message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) + } + + + elapsed <- Sys.time() - start_time + + if (is.null(parma_res)) { + stop("parma::cmaes failed (object NULL).") + } + + sol <- NULL + if (!is.null(parma_res$bestever$x)) { + sol <- parma_res$bestever$x + } else if (!is.null(parma_res$par)) { + sol <- parma_res$par + } else if (!is.null(parma_res$xbest)) { + sol <- parma_res$xbest + } + if (is.null(sol)) { + stop("Impossible to find a solution in the object sent by parma::cmaes().") + } + final_values <- as.numeric(sol) + names(final_values) <- param_names + + # IMPORTANT: make post_treat_global_optim() happy (strict !=) + min_crit_value <- main_crit(final_values, crit_options = crit_options) + + # Keep CMAES native objectives for debugging (do NOT use for post_treat) + cmaes_bestever_f <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f else NA_real_ + cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ + + + # Get the estimated values ==> matrix of 1 * nb_params + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) + colnames(est_values) <- param_names + if (ncol(est_values) == nb_params) { + colnames(est_values) <- param_names + } + + # Final solution + trace_df <- NULL + if (keep_trace && trace_env$eval > 0L) { + Xmat <- do.call(rbind, trace_env$X) + colnames(Xmat) <- param_names + + trace_df <- as.data.frame(Xmat) + trace_df$crit <- trace_env$crit + trace_df$eval <- seq_len(trace_env$eval) + lambda <- ctrl$options$PopSize + if (is.null(lambda) || is.na(lambda)) + lambda <- 30L + lambda <- as.integer(lambda) + + trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES + trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory + trace_df$rep <- 1L + trace_df$method <- "cmaes" + } + + res <- list( + final_values = final_values, + init_values = init_values, + est_values = est_values, + min_crit_value = min_crit_value, + CMAES = parma_res, + trace_df = trace_df, + cmaes_bestever_f = cmaes_bestever_f, + cmaes_objective = cmaes_objective, + counteval = if (!is.null(parma_res$counteval)) parma_res$counteval else NA_integer_, + stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, + out_dir = crit_options$out_dir + + ) + + return(res) +} + diff --git a/man/wrap_cmaes.Rd b/man/wrap_cmaes.Rd new file mode 100644 index 00000000..56ad4462 --- /dev/null +++ b/man/wrap_cmaes.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wrap_parma.R +\name{wrap_cmaes} +\alias{wrap_cmaes} +\title{A wrapper for parma::cmaes package (CMA-ES)} +\usage{ +wrap_cmaes(optim_options, param_info, crit_options) +} +\arguments{ +\item{optim_options}{see description in estim_param} + +\item{param_info}{see description in estim_param} + +\item{crit_options}{List containing several arguments given to \code{estim_param} +function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function}, +\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir} +that must be passed to main_crit function by the methods wrappers.} +} +\value{ +list with final_values, init_values, est_values, GE, min_crit_value +} +\description{ +A wrapper for parma::cmaes package (CMA-ES) +} +\keyword{internal} diff --git a/tests/testthat/test_estim_param_cmaes.R b/tests/testthat/test_estim_param_cmaes.R new file mode 100644 index 00000000..9d3cac4d --- /dev/null +++ b/tests/testthat/test_estim_param_cmaes.R @@ -0,0 +1,180 @@ +context("Test the estim_param function using a toy model") +library(tibble) + +# Define a toy model and its wrapper +toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) { + # Simulate biomass and yield with a simple logistic model + # + # Arguments + # - rB: growth rate + # - Bmin: minimum biomass + # - Bmax: maximum biomass + # - h: harvest rate + # - Bini: initial biomass + # + # Value + # A list containing the simulated biomass and yield + # + biom <- numeric(nb_days) + biom[1] <- Bini + for (i in 2:nb_days) { + biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax) + } + yield <- h * biom + return(list(biom = biom, yield = yield)) +} + +toymodel_wrapper <- function(param_values = NULL, situation, + model_options, ...) { + # A wrapper for toymodel + # + # Arguments + # - param_values: (optional) a named vector or a tibble containing the + # values of the toymodel parameters to force. + # - situation: Vector of situations names for which results must be + # returned. + # In this case, the names of the situations are coded as "year_suffix" + # - model_options: a tibble containing the begin and end date of simulation + # for each situation. + # - ...: mandatory since CroptimizR will give additional arguments not used + # here + # + # Value: + # A named list of tibble per situation. + # Each tibble contains columns: + # - Date (POSIXct dates of simulated results), + # - One column per simulated variable (biomass and yield) + # + results <- list( + sim_list = setNames( + vector("list", length(situation)), + nm = situation + ), + error = FALSE + ) + attr(results$sim_list, "class") <- "cropr_simulation" + + # Remove dummy in param_values, just here for specific tests + if (!is.null(param_values) && "dummy" %in% names(param_values)) { + param_values <- param_values[!names(param_values) %in% "dummy"] + } + + for (sit in situation) { + # Retrieve begin and end date from situation name + begin_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(begin_date) + end_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(end_date) + + date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day") + + if (!all(names(param_values) %in% c( + "rB", "Bmin", "Bmax", "h", "Bini" + ))) { + warning( + paste( + "Unknown parameters in param_values:", + paste(names(param_values), collapse = ",") + ) + ) + results$error <- TRUE + return(results) + } + + # Call the toymodel with varying arguments depending on what is given in + # param_values + res <- do.call( + "toymodel", + c(nb_days = length(date_sequence), as.list(param_values)) + ) + + # Fill the result variable + results$sim_list[[sit]] <- + dplyr::tibble( + Date = date_sequence, + biomass = res$biom, + yield = res$yield + ) + } + return(results) +} + +# To make these function accessible to the test environment +.GlobalEnv$toymodel <- toymodel +.GlobalEnv$toymodel_wrapper <- toymodel_wrapper + +# Use setup() to define shared data for all tests in this file +setup({ + # These variables will be available in all tests + .GlobalEnv$model_options <- tibble::tibble( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")), + end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18")) + ) + .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7) + + # Generate synthetic observations + tmp <- toymodel_wrapper( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + param_values = param_true_values, + model_options = model_options + ) + + .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) { + return(dplyr::filter(x, dplyr::row_number() %% 10 == 0)) + }) +}) + +# For cleaning up after tests +teardown({ + if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv) + if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv) + if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv) +}) + + + +# Test estim_param with CMA-ES (parma) +test_that("estim_param 1 step OLS criterion", { + param_info <- list( + rB = list(lb = 0, ub = 1), + h = list(lb = 0, ub = 1), + Bmax = list(lb = 5, ub = 15) + ) + + optim_options <- list( + ranseed = 1234, + return_trace = TRUE, + insigma = 0.2, + ctrl = list( + options = list(MaxIter = 200, PopSize = 30, Trace = 0), + CMA = list(active = 1) + ) + ) + + res <- estim_param( + obs_list = obs_synth, + model_function = toymodel_wrapper, + model_options = model_options, + optim_method = "cmaes", + crit_function = crit_ols, + optim_options = optim_options, + param_info = param_info, + obs_var = c("biomass", "yield"), + situation = c("sit1_2000", "sit1_2001", "sit2_2003"), + out_dir = tempdir() + ) + + expect_equal(res$final_values[["rB"]], + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 + ) + expect_equal(res$final_values[["h"]], + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 + ) + expect_equal(res$final_values[["Bmax"]], + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 + ) +}) From 886febb34b4407bda29d8b093d7bb680213880c3 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Tue, 16 Dec 2025 15:53:46 +0100 Subject: [PATCH 09/16] Harmonize global optim plot - Remove lines between points - Color points by criterion - Unify CMAES / DE / graDiEnt plotting --- R/global_optim_function.R | 31 +++++++++++++++---------------- R/wrap_DE.R | 7 ++++++- R/wrap_cmaes.R | 27 +++++++++++---------------- R/wrap_graDiEnt.R | 1 + 4 files changed, 33 insertions(+), 33 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index 8ba6a259..5ed10e9a 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -166,7 +166,7 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) df = df_2d, param_info = param_info, iter_or_eval = "eval", - fill = "eval", + fill = "crit", crit_log = FALSE, lines = FALSE, rep_label = "end" @@ -282,8 +282,6 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), ) } for (iind in unique(df$ind)) { - p[[param_name]] <- p[[param_name]] + - geom_line(data = filter(df, ind == iind)) if (ind_label[1] == "begin_end" || ind_label[1] == "begin") { p[[param_name]] <- p[[param_name]] + geom_label(aes(label = ind), @@ -304,9 +302,10 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), } df$ind <- as.factor(df$ind) - p[["criterion"]] <- ggplot(df, aes_string( - x = iter_or_eval[1], y = "crit", - color = "ind" + p[["criterion"]] <- ggplot(df, aes( + x = !!rlang::sym(iter_or_eval[1]), + y = crit, + color = crit )) + labs( title = paste0( @@ -315,20 +314,16 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), ), y = "Minimized criterion", x = paste(lab, "number"), - fill = "Individual" + color = "Criterion" ) + theme(plot.title = element_text(hjust = 0.5)) + - geom_point(alpha = 0.5) + geom_point(alpha = 0.6) - for (iind in unique(df$ind)) { + if (have_crit) { p[["criterion"]] <- p[["criterion"]] + - geom_line(data = filter(df, ind == iind)) + - geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), - size = 3 - ) + - geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)) + scale_color_gradient2( + midpoint = mid, low = "blue", mid = "yellow", + high = "red", space = "Lab", trans = trans ) } @@ -336,6 +331,10 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), p[["criterion"]] <- p[["criterion"]] + scale_y_log10() } + if (crit_log) { + p[["criterion"]] <- p[["criterion"]] + scale_y_log10() + } + return(p) } diff --git a/R/wrap_DE.R b/R/wrap_DE.R index 9ec5455d..d0d3344f 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -104,7 +104,12 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { colnames(pop_it) <- param_names pop_it$ind <- seq_len(NP) pop_it$iter <- it - pop_it$crit <- bestvalit[it] + crit_pop <- apply(pop_it[, param_names, drop = FALSE], 1, function(par) { + par <- as.numeric(par) + names(par) <- param_names + main_crit(par, crit_options = crit_options) + }) + pop_it$crit <- crit_pop idx <- seq_len(NP) pop_it$eval <- eval_counter + idx eval_counter <- eval_counter + NP diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R index 67a74e44..499d67d2 100644 --- a/R/wrap_cmaes.R +++ b/R/wrap_cmaes.R @@ -18,32 +18,28 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { optim_options$return_trace <- TRUE } optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it - algorithm <- "parma_cmaes" + algorithm <- "cmaes" if (is.null(param_info) || is.null(crit_options)) { - init_nb <- if (!is.null(optim_options$n_particles)) optim_options$n_particles else NA_integer_ + return(list( package = "parma", family = "Global", - method = "cmaes", init_values_nb = init_nb + method = "cmaes", init_values_nb = init_nb <- if (!is.null(optim_options$ctrl$options$PopSize)) optim_options$ctrl$options$PopSize else NA_integer_ )) } - param_names <- get_params_names(param_info, short_list = TRUE) + param_names <- get_params_names(param_info) nb_params <- length(param_names) bounds <- get_params_bounds(param_info) init_values <- get_init_values(param_info) control_params <- optim_options$ctrl - # base ctrl (toujours valide) + ctrl <- parma::cmaes.control() - # si l'utilisateur a déjà fourni un ctrl complet (déjà construit) - if (!is.null(control_params) && !is.null(control_params$options) && !is.null(control_params$CMA) && - !is.list(control_params$options)) { - # (cas rare) on ignore, car options devrait être une liste - } - # si l'utilisateur a fourni une liste ctrl = list(options=..., CMA=...) + + # if ctrl = list(options=..., CMA=...) existe if (!is.null(control_params) && is.list(control_params) && (!is.null(control_params$options) || !is.null(control_params$CMA))) { @@ -59,7 +55,6 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } } - # si l'utilisateur a fourni directement une liste d'options (sans options/CMA) } else if (!is.null(control_params) && is.list(control_params)) { for (nm in names(control_params)) { ctrl$options[[nm]] <- control_params[[nm]] @@ -97,21 +92,21 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) - # 🟩 parma::cmaes expects a numeric vector for 'pars' (length = nb_params) + #because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) x0 <- init_values - # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) + # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) our case data frame if (is.data.frame(x0)) x0 <- as.matrix(x0) if (is.matrix(x0)) { - # take the first row as starting point (you can also use colMeans) + # take the first row as starting point x0 <- x0[1, , drop = TRUE] } x0 <- as.numeric(x0) names(x0) <- param_names - # safety: if still wrong length, fallback to mid-bounds + # safety: if still initial-value isn(t provided) if (length(x0) != nb_params) { x0 <- (bounds$lb + bounds$ub) / 2 names(x0) <- param_names diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index 03646ff2..a69d5026 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -149,6 +149,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df$eval <- eval_counter + idx eval_counter <- eval_counter + n_pop df$method <- "graDiEnt" + df$rep <- 1L df_list[[it]] <- df } trace_df <- dplyr::bind_rows(df_list) From 6bfe42237f33377c2a0088172ef652979740783c Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Wed, 17 Dec 2025 08:42:30 +0100 Subject: [PATCH 10/16] Harmonize global optimization visualizations --- R/global_optim_function.R | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index 5ed10e9a..8fe58f0f 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -281,22 +281,8 @@ plot_valuesVSit_go <- function(df, param_info, iter_or_eval = c("iter", "eval"), high = "red", space = "Lab", trans = trans ) } - for (iind in unique(df$ind)) { - if (ind_label[1] == "begin_end" || ind_label[1] == "begin") { - p[[param_name]] <- p[[param_name]] + - geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == min(.data$eval)), - size = 3 - ) - } - if (ind_label[1] == "begin_end" || ind_label[1] == "end") { - p[[param_name]] <- p[[param_name]] + - geom_label(aes(label = ind), - data = filter(df, ind == iind) %>% filter(eval == max(.data$eval)), - size = 3 - ) - } - } + + p[[param_name]] <- p[[param_name]] + ylim(minvalue[param_name], maxvalue[param_name]) } From 070071af140d6eed8dbc5653838772f48444802c Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Wed, 17 Dec 2025 07:50:55 +0000 Subject: [PATCH 11/16] Style code (GHA) --- R/optim_switch.R | 13 +-- R/wrap_cmaes.R | 101 +++++++++++------------- R/wrap_graDiEnt.R | 66 ++++++++-------- tests/testthat/test_estim_param_cmaes.R | 13 ++- 4 files changed, 90 insertions(+), 103 deletions(-) diff --git a/R/optim_switch.R b/R/optim_switch.R index 58b21b2a..1f3e0e9b 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -23,7 +23,6 @@ optim_switch <- function(...) { res <- list() flag_error <- FALSE if (nargs() > 2) { - on.exit({ #### if (!is.null(crit_options)) { @@ -215,12 +214,10 @@ optim_switch <- function(...) { out_dir = crit_options$out_dir ) } - } else if (optim_method == "cmaes") { - - + } else if (optim_method == "cmaes") { if (is.null(optim_options$n_particles) || is.na(optim_options$n_particles)) { if (!is.null(optim_options$ctrl$options$PopSize) && - !is.na(optim_options$ctrl$options$PopSize)) { + !is.na(optim_options$ctrl$options$PopSize)) { optim_options$n_particles <- as.integer(optim_options$ctrl$options$PopSize) } else { optim_options$n_particles <- 30L @@ -255,9 +252,8 @@ optim_switch <- function(...) { optim_results = res, out_dir = crit_options$out_dir ) - }} - -else { + } + } else { flag_unknown_method <- TRUE }, error = function(cond) { @@ -265,7 +261,6 @@ else { traceback(20) stop(cond) } - ) if (flag_unknown_method) { diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R index 499d67d2..32e9a482 100644 --- a/R/wrap_cmaes.R +++ b/R/wrap_cmaes.R @@ -7,8 +7,8 @@ #' @keywords internal #' wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { - if(is.null((ranseed <- optim_options$ranseed))) { - ranseed = NULL + if (is.null((ranseed <- optim_options$ranseed))) { + ranseed <- NULL } n_params <- optim_options$n_params if (is.null(n_params)) { @@ -21,7 +21,6 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { algorithm <- "cmaes" if (is.null(param_info) || is.null(crit_options)) { - return(list( package = "parma", family = "Global", method = "cmaes", init_values_nb = init_nb <- if (!is.null(optim_options$ctrl$options$PopSize)) optim_options$ctrl$options$PopSize else NA_integer_ @@ -38,11 +37,9 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { ctrl <- parma::cmaes.control() - # if ctrl = list(options=..., CMA=...) existe if (!is.null(control_params) && is.list(control_params) && - (!is.null(control_params$options) || !is.null(control_params$CMA))) { - + (!is.null(control_params$options) || !is.null(control_params$CMA))) { if (!is.null(control_params$options)) { for (nm in names(control_params$options)) { ctrl$options[[nm]] <- control_params$options[[nm]] @@ -54,16 +51,15 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { ctrl$CMA[[nm]] <- control_params$CMA[[nm]] } } - } else if (!is.null(control_params) && is.list(control_params)) { for (nm in names(control_params)) { ctrl$options[[nm]] <- control_params[[nm]] } } - if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){ + if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)) { crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter) - } + } start_time <- Sys.time() @@ -74,25 +70,25 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { trace_env$X <- list() keep_trace <- isTRUE(optim_options$return_trace) ObjFun <- function(x) { - x <- as.numeric(x) - names(x) <- param_names - val <- main_crit(x, crit_options = crit_options) - - if (keep_trace) { - trace_env$eval <- trace_env$eval + 1L - trace_env$crit[trace_env$eval] <- val - trace_env$X[[trace_env$eval]] <- x - } - return(val) + x <- as.numeric(x) + names(x) <- param_names + val <- main_crit(x, crit_options = crit_options) + + if (keep_trace) { + trace_env$eval <- trace_env$eval + 1L + trace_env$crit[trace_env$eval] <- val + trace_env$X[[trace_env$eval]] <- x } + return(val) + } insigma <- optim_options$insigma if (is.null(insigma)) { - # default value if it is not provided - insigma <- rep(0.3, nb_params) + # default value if it is not provided + insigma <- rep(0.3, nb_params) } if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) - #because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) + # because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) x0 <- init_values # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) our case data frame @@ -127,8 +123,8 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } ) if (isTRUE(optim_options$debug)) { - message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) - } + message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) + } elapsed <- Sys.time() - start_time @@ -139,15 +135,15 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { sol <- NULL if (!is.null(parma_res$bestever$x)) { - sol <- parma_res$bestever$x - } else if (!is.null(parma_res$par)) { - sol <- parma_res$par - } else if (!is.null(parma_res$xbest)) { - sol <- parma_res$xbest - } + sol <- parma_res$bestever$x + } else if (!is.null(parma_res$par)) { + sol <- parma_res$par + } else if (!is.null(parma_res$xbest)) { + sol <- parma_res$xbest + } if (is.null(sol)) { - stop("Impossible to find a solution in the object sent by parma::cmaes().") - } + stop("Impossible to find a solution in the object sent by parma::cmaes().") + } final_values <- as.numeric(sol) names(final_values) <- param_names @@ -156,7 +152,7 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { # Keep CMAES native objectives for debugging (do NOT use for post_treat) cmaes_bestever_f <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f else NA_real_ - cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ + cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ # Get the estimated values ==> matrix of 1 * nb_params @@ -173,22 +169,23 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { # Final solution trace_df <- NULL if (keep_trace && trace_env$eval > 0L) { - Xmat <- do.call(rbind, trace_env$X) - colnames(Xmat) <- param_names - - trace_df <- as.data.frame(Xmat) - trace_df$crit <- trace_env$crit - trace_df$eval <- seq_len(trace_env$eval) - lambda <- ctrl$options$PopSize - if (is.null(lambda) || is.na(lambda)) - lambda <- 30L - lambda <- as.integer(lambda) - - trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES - trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory - trace_df$rep <- 1L - trace_df$method <- "cmaes" - } + Xmat <- do.call(rbind, trace_env$X) + colnames(Xmat) <- param_names + + trace_df <- as.data.frame(Xmat) + trace_df$crit <- trace_env$crit + trace_df$eval <- seq_len(trace_env$eval) + lambda <- ctrl$options$PopSize + if (is.null(lambda) || is.na(lambda)) { + lambda <- 30L + } + lambda <- as.integer(lambda) + + trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES + trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory + trace_df$rep <- 1L + trace_df$method <- "cmaes" + } res <- list( final_values = final_values, @@ -198,13 +195,11 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { CMAES = parma_res, trace_df = trace_df, cmaes_bestever_f = cmaes_bestever_f, - cmaes_objective = cmaes_objective, + cmaes_objective = cmaes_objective, counteval = if (!is.null(parma_res$counteval)) parma_res$counteval else NA_integer_, - stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, + stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, out_dir = crit_options$out_dir - ) return(res) } - diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index a69d5026..c5552fa2 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -7,38 +7,38 @@ #' @keywords internal #' wrap_graDiEnt <- function(optim_options, param_info, crit_options) { - if(is.null((ranseed <- optim_options$ranseed))) { - ranseed = NULL - } + if (is.null((ranseed <- optim_options$ranseed))) { + ranseed <- NULL + } if (is.null((n_params <- optim_options$n_params))) { stop( "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)" ) - } + } if (is.null((trace <- optim_options$return_trace))) { optim_options$return_trace <- TRUE - } + } optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it algorithm <- "graDiEnt" # return requested information if only optim_options is given in argument if (nargs() == 1 && methods::hasArg(optim_options)) { - init_nb <- if (!is.null(optim_options$n_particles)) { - optim_options$n_particles - } else { - NA_integer_ - } - return(list( + init_nb <- if (!is.null(optim_options$n_particles)) { + optim_options$n_particles + } else { + NA_integer_ + } + return(list( package = "graDiEnt", family = "Global", method = "SQGDE", init_values_nb = init_nb - )) - } + )) + } param_names <- get_params_names(param_info) nb_params <- length(param_names) bounds <- get_params_bounds(param_info) init_values <- get_init_values(param_info) - control_params <- do.call(graDiEnt::GetAlgoParams,optim_options) + control_params <- do.call(graDiEnt::GetAlgoParams, optim_options) crit_options$tot_max_eval <- control_params$n_particles * control_params$n_iter start_time <- Sys.time() @@ -59,7 +59,8 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { ), error = function(e) { warning(sprintf("GraDiEnt failed: %s", e$message)) - NULL } + NULL + } ) elapsed <- Sys.time() - start_time @@ -76,7 +77,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { final_values <- bounds$lb + u_sol * range_bounds names(final_values) <- param_names min_crit_value <- SQGDE$solution_weight - if (is.null(min_crit_value) ) { + if (is.null(min_crit_value)) { min_crit_value <- main_crit(final_values, crit_options) } # Get the estimated values ==> matrix of 1 * nb_params @@ -84,11 +85,11 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { final_values, nrow = 1L, dimnames = list(NULL, param_names) - ) + ) colnames(est_values) <- param_names if (ncol(est_values) == nb_params) { colnames(est_values) <- param_names - } + } # Final solution if (!is.null(SQGDE$particles_trace)) { @@ -101,28 +102,27 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { dimnames = list(NULL, param_names) ) } else { - n_it <- dim(tr)[1] + n_it <- dim(tr)[1] last_pop <- tr[n_it, , ] - est_u <- as.matrix(last_pop) + est_u <- as.matrix(last_pop) colnames(est_u) <- param_names # transformation est_values <- sweep(est_u, 2, range_bounds, `*`) est_values <- sweep(est_values, 2, bounds$lb, `+`) colnames(est_values) <- param_names - } - } else { - est_values <- matrix( + } + } else { + est_values <- matrix( final_values, nrow = 1L, dimnames = list(NULL, param_names) - ) - } + ) + } trace_df <- NULL - if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) - { + if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) { tr <- SQGDE$particles_trace - n_it <- dim(tr)[1] + n_it <- dim(tr)[1] n_pop <- dim(tr)[2] df_list <- vector("list", n_it) eval_counter <- 0 @@ -137,13 +137,12 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df <- as.data.frame(pop_x) - df$ind <- seq_len(n_pop) + df$ind <- seq_len(n_pop) df$iter <- it crit_pop <- apply(pop_x, 1, function(par) { names(par) <- param_names main_crit(par, crit_options = crit_options) - } - ) + }) df$crit <- crit_pop idx <- seq_len(n_pop) df$eval <- eval_counter + idx @@ -151,9 +150,9 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { df$method <- "graDiEnt" df$rep <- 1L df_list[[it]] <- df - } - trace_df <- dplyr::bind_rows(df_list) } + trace_df <- dplyr::bind_rows(df_list) + } res <- list( final_values = final_values, init_values = init_values, @@ -165,4 +164,3 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { return(res) } - diff --git a/tests/testthat/test_estim_param_cmaes.R b/tests/testthat/test_estim_param_cmaes.R index 9d3cac4d..d8c560a4 100644 --- a/tests/testthat/test_estim_param_cmaes.R +++ b/tests/testthat/test_estim_param_cmaes.R @@ -133,7 +133,6 @@ teardown({ }) - # Test estim_param with CMA-ES (parma) test_that("estim_param 1 step OLS criterion", { param_info <- list( @@ -166,15 +165,15 @@ test_that("estim_param 1 step OLS criterion", { ) expect_equal(res$final_values[["rB"]], - param_true_values[["rB"]], - tolerance = param_true_values[["rB"]] * 1e-2 + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 ) expect_equal(res$final_values[["h"]], - param_true_values[["h"]], - tolerance = param_true_values[["h"]] * 1e-2 + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 ) expect_equal(res$final_values[["Bmax"]], - param_true_values[["Bmax"]], - tolerance = param_true_values[["Bmax"]] * 1e-2 + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 ) }) From 67f34a6479f5f36f77c6dc67ca3f427f7177783d Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Wed, 17 Dec 2025 09:47:38 +0100 Subject: [PATCH 12/16] parma added to Description --- DESCRIPTION | 1 + man/wrap_cmaes.Rd | 4 ++-- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2f2b35ff..1bfb3d1d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -46,6 +46,7 @@ Imports: lubridate, methods, nloptr, + parma, purrr, rlang, rstudioapi, diff --git a/man/wrap_cmaes.Rd b/man/wrap_cmaes.Rd index 56ad4462..deb10b0a 100644 --- a/man/wrap_cmaes.Rd +++ b/man/wrap_cmaes.Rd @@ -1,10 +1,10 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/wrap_parma.R +% Please edit documentation in R/wrap_cmaes.R \name{wrap_cmaes} \alias{wrap_cmaes} \title{A wrapper for parma::cmaes package (CMA-ES)} \usage{ -wrap_cmaes(optim_options, param_info, crit_options) +wrap_cmaes(optim_options, param_info = NULL, crit_options) } \arguments{ \item{optim_options}{see description in estim_param} From fb2c0304212786c00e29e3b21222c7b45b029302 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Wed, 17 Dec 2025 18:35:00 +0100 Subject: [PATCH 13/16] herit problems solved --- R/global_optim_function.R | 3 ++- R/multi_steps_functions.R | 4 ---- man/plot_valuesVSit_go.Rd | 23 +---------------------- 3 files changed, 3 insertions(+), 27 deletions(-) diff --git a/R/global_optim_function.R b/R/global_optim_function.R index 8fe58f0f..e57e161c 100644 --- a/R/global_optim_function.R +++ b/R/global_optim_function.R @@ -184,7 +184,8 @@ plot_global_optim <- function(optim_options, param_info, optim_results, out_dir) #' @title Create plots of parameters and criterion values per iteration or #' evaluation number #' -#' @inheritParams estim_param +#' @param param_info List describing the parameters to estimate + #' @param df Data.frame containing values of each individual and iteration parameters values (one column per #' estimated parameter), criterion (crit column), individual index (ind), #' iteration number (iter) and evaluation number (eval) diff --git a/R/multi_steps_functions.R b/R/multi_steps_functions.R index 967a2520..61f741ce 100644 --- a/R/multi_steps_functions.R +++ b/R/multi_steps_functions.R @@ -1,7 +1,5 @@ #' @title Summarizes results of multi-step procedure #' -#' @inheritParams estim_param -#' #' @param results_multi_step Results of the multi_step procedure as returned by post_treat_multi_step #' #' @param path_results Folder path where results of the multi-step optimization process can be found @@ -45,8 +43,6 @@ summary_multi_step <- function(results_multi_step, path_results) { #' @title Post-treat results of multi-step procedure #' -#' @inheritParams estim_param -#' #' @param optim_results_list List of results returned for each step of the parameter estimation process #' #' @return List of estimated and forced parameters values diff --git a/man/plot_valuesVSit_go.Rd b/man/plot_valuesVSit_go.Rd index 9034cbbd..7fa7e99b 100644 --- a/man/plot_valuesVSit_go.Rd +++ b/man/plot_valuesVSit_go.Rd @@ -21,28 +21,7 @@ iteration number (iter) and evaluation number (eval) See Details section for comments about the difference between evaluations and iterations.} -\item{param_info}{Information on the parameters to estimate. -Either -a list containing: -\itemize{ -\item \code{ub} and \code{lb}, named vectors of upper and lower bounds (-Inf and Inf can be used if init_values is provided), -\item \code{default}, named vectors of default values (optional, corresponding parameters are set to these values when the parameter is part of the \code{candidate_param} list and when it is not estimated ; these values are also used as first initial values when the parameters are estimated) -\item \code{init_values}, a data.frame containing initial -values to test for the parameters (optional, if not provided, or if less values -than number of repetitions of the minimization are provided, the, or part -of the, initial values will be randomly generated using LHS sampling within -parameter bounds). -} - -or a named list containing for each parameter: -\itemize{ -\item \code{sit_list}, list the groups of situations for which the current estimated -parameter must take different values (see \href{https://sticsrpacks.github.io/CroptimizR/articles/Parameter_estimation_Specific_and_Varietal.html}{here} -for an example), -\item \code{ub} and \code{lb}, vectors of upper and lower bounds (one value per group), -\item \code{init_values}, the list of initial values per group (data.frame, one column per group, optional). -\item \code{default}, vector of default values per group (optional, the parameter is set to its default value when it is part of the \code{candidate_param} list and when it is not estimated ; the default value is also used as first initial value when the parameter is estimated) -}} +\item{param_info}{List describing the parameters to estimate} \item{iter_or_eval}{Values of the x axis: "iter" for iteration number, "eval" for evaluation number} From 165a6e67a323c94dc82434f36a50550c017a93ad Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Sat, 10 Jan 2026 21:44:31 +0100 Subject: [PATCH 14/16] Update wrappers and evaluations counting --- NAMESPACE | 1 + R/wrap_DE.R | 54 +++++++++++---------- R/wrap_cmaes.R | 118 +++++++++++++++++++++++++--------------------- R/wrap_graDiEnt.R | 60 +++++++++++------------ man/wrap_cmaes.Rd | 4 +- 5 files changed, 124 insertions(+), 113 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index c5ef7162..dff958e2 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -17,6 +17,7 @@ export(plot_valuesVSit_2D) export(plot_valuesVSit_2D_go) export(plot_valuesVSit_go) export(test_wrapper) +export(wrap_cmaes) importFrom(BayesianTools,MAP) importFrom(BayesianTools,applySettingsDefault) importFrom(BayesianTools,correlationPlot) diff --git a/R/wrap_DE.R b/R/wrap_DE.R index d0d3344f..acb34e31 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -53,9 +53,18 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { start_time <- Sys.time() + .trace_env <- new.env(parent = emptyenv()) + .trace_env$x_list <- list() + .trace_env$crit_list <- list() + .trace_env$k <- 0L + if (!is.null(ranseed)) set.seed(ranseed) fn_de <- function(x) { - main_crit(x, crit_options = crit_options) + val <- main_crit(x, crit_options = crit_options) + .trace_env$k <- .trace_env$k + 1L + .trace_env$x_list[[.trace_env$k]] <- as.numeric(x) + .trace_env$crit_list[[.trace_env$k]] <- val + return(val) } DE <- tryCatch( DEoptim::DEoptim( @@ -94,31 +103,26 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { trace_df <- NULL if (!is.null(DE$member$storepop)) { storepop <- DE$member$storepop - bestvalit <- DE$member$bestvalit - itermax <- length(storepop) - NP <- nrow(storepop[[1]]) - df_list <- vector("list", itermax) - eval_counter <- 0 - for (it in seq_len(itermax)) { - pop_it <- as.data.frame(storepop[[it]]) - colnames(pop_it) <- param_names - pop_it$ind <- seq_len(NP) - pop_it$iter <- it - crit_pop <- apply(pop_it[, param_names, drop = FALSE], 1, function(par) { - par <- as.numeric(par) - names(par) <- param_names - main_crit(par, crit_options = crit_options) - }) - pop_it$crit <- crit_pop - idx <- seq_len(NP) - pop_it$eval <- eval_counter + idx - eval_counter <- eval_counter + NP - pop_it$rep <- 1 - pop_it$method <- "DEoptim" - df_list[[it]] <- pop_it + itermax <- length(storepop) + NP <- nrow(storepop[[1]]) + + pop_all <- do.call(rbind, lapply(seq_len(itermax), function(it) { + m <- as.data.frame(storepop[[it]]) + colnames(m) <- param_names + m$ind <- seq_len(NP) + m$iter <- it + m + })) + crit_vec <- as.numeric(unlist(.trace_env$crit_list)) + n <- min(nrow(pop_all), length(crit_vec)) + pop_all <- pop_all[seq_len(n), , drop = FALSE] + pop_all$crit <- crit_vec[seq_len(n)] + + pop_all$eval <- seq_len(n) + pop_all$rep <- 1L + pop_all$method <- "DEoptim" + trace_df <- pop_all } - trace_df <- dplyr::bind_rows(df_list) - } res <- list( final_values = final_values, diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R index 32e9a482..f80247c7 100644 --- a/R/wrap_cmaes.R +++ b/R/wrap_cmaes.R @@ -1,13 +1,16 @@ -#' @title A wrapper for parma::cmaes package (CMA-ES) +#' Wrap CMA-ES using parma #' #' @inheritParams optim_switch #' #' @return list with final_values, init_values, est_values, GE, min_crit_value #' @importFrom parma cmaes cmaes.control #' @keywords internal -#' +#' @export wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { - if (is.null((ranseed <- optim_options$ranseed))) { + if (!requireNamespace("parma", quietly = TRUE)) { + stop("the package 'parma' (not installed) becouse it required a dependencies that doesn't exist in the cluster") + } + if(is.null((ranseed <- optim_options$ranseed))) { ranseed <- NULL } n_params <- optim_options$n_params @@ -21,6 +24,7 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { algorithm <- "cmaes" if (is.null(param_info) || is.null(crit_options)) { + return(list( package = "parma", family = "Global", method = "cmaes", init_values_nb = init_nb <- if (!is.null(optim_options$ctrl$options$PopSize)) optim_options$ctrl$options$PopSize else NA_integer_ @@ -37,9 +41,11 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { ctrl <- parma::cmaes.control() + # if ctrl = list(options=..., CMA=...) existe if (!is.null(control_params) && is.list(control_params) && - (!is.null(control_params$options) || !is.null(control_params$CMA))) { + (!is.null(control_params$options) || !is.null(control_params$CMA))) { + if (!is.null(control_params$options)) { for (nm in names(control_params$options)) { ctrl$options[[nm]] <- control_params$options[[nm]] @@ -51,15 +57,17 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { ctrl$CMA[[nm]] <- control_params$CMA[[nm]] } } + + # si l'utilisateur a fourni directement une liste d'options (sans options/CMA) } else if (!is.null(control_params) && is.list(control_params)) { for (nm in names(control_params)) { ctrl$options[[nm]] <- control_params[[nm]] } } - if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)) { + if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){ crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter) - } + } start_time <- Sys.time() @@ -70,25 +78,25 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { trace_env$X <- list() keep_trace <- isTRUE(optim_options$return_trace) ObjFun <- function(x) { - x <- as.numeric(x) - names(x) <- param_names - val <- main_crit(x, crit_options = crit_options) - - if (keep_trace) { - trace_env$eval <- trace_env$eval + 1L - trace_env$crit[trace_env$eval] <- val - trace_env$X[[trace_env$eval]] <- x - } + x <- as.numeric(x) + names(x) <- param_names + val <- main_crit(x, crit_options = crit_options) + + if (keep_trace) { + trace_env$eval <- trace_env$eval + 1L + trace_env$crit[trace_env$eval] <- val + trace_env$X[[trace_env$eval]] <- x + } return(val) - } + } insigma <- optim_options$insigma if (is.null(insigma)) { - # default value if it is not provided - insigma <- rep(0.3, nb_params) + # default value if it is not provided + insigma <- rep(0.3, nb_params) } if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) - # because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) + #because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) x0 <- init_values # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) our case data frame @@ -102,7 +110,7 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { x0 <- as.numeric(x0) names(x0) <- param_names - # safety: if still initial-value isn(t provided) + if (length(x0) != nb_params) { x0 <- (bounds$lb + bounds$ub) / 2 names(x0) <- param_names @@ -123,8 +131,8 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } ) if (isTRUE(optim_options$debug)) { - message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) - } + message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) + } elapsed <- Sys.time() - start_time @@ -135,27 +143,28 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { sol <- NULL if (!is.null(parma_res$bestever$x)) { - sol <- parma_res$bestever$x - } else if (!is.null(parma_res$par)) { - sol <- parma_res$par - } else if (!is.null(parma_res$xbest)) { - sol <- parma_res$xbest - } + sol <- parma_res$bestever$x + } else if (!is.null(parma_res$par)) { + sol <- parma_res$par + } else if (!is.null(parma_res$xbest)) { + sol <- parma_res$xbest + } if (is.null(sol)) { - stop("Impossible to find a solution in the object sent by parma::cmaes().") - } + stop("Impossible to find a solution in the object sent by parma::cmaes().") + } final_values <- as.numeric(sol) names(final_values) <- param_names - # IMPORTANT: make post_treat_global_optim() happy (strict !=) - min_crit_value <- main_crit(final_values, crit_options = crit_options) + # for post_treat_global_optim() + min_crit_value <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f + else main_crit(final_values, crit_options) + - # Keep CMAES native objectives for debugging (do NOT use for post_treat) cmaes_bestever_f <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f else NA_real_ - cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ + cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ - # Get the estimated values ==> matrix of 1 * nb_params + # Getting the estimated values est_values <- matrix( final_values, nrow = 1L, @@ -169,23 +178,22 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { # Final solution trace_df <- NULL if (keep_trace && trace_env$eval > 0L) { - Xmat <- do.call(rbind, trace_env$X) - colnames(Xmat) <- param_names - - trace_df <- as.data.frame(Xmat) - trace_df$crit <- trace_env$crit - trace_df$eval <- seq_len(trace_env$eval) - lambda <- ctrl$options$PopSize - if (is.null(lambda) || is.na(lambda)) { - lambda <- 30L - } - lambda <- as.integer(lambda) - - trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES - trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory - trace_df$rep <- 1L - trace_df$method <- "cmaes" - } + Xmat <- do.call(rbind, trace_env$X) + colnames(Xmat) <- param_names + + trace_df <- as.data.frame(Xmat) + trace_df$crit <- trace_env$crit + trace_df$eval <- seq_len(trace_env$eval) + lambda <- ctrl$options$PopSize + if (is.null(lambda) || is.na(lambda)) + lambda <- 30L + lambda <- as.integer(lambda) + + trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES + trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory + trace_df$rep <- 1L + trace_df$method <- "cmaes" + } res <- list( final_values = final_values, @@ -195,11 +203,13 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { CMAES = parma_res, trace_df = trace_df, cmaes_bestever_f = cmaes_bestever_f, - cmaes_objective = cmaes_objective, + cmaes_objective = cmaes_objective, counteval = if (!is.null(parma_res$counteval)) parma_res$counteval else NA_integer_, - stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, + stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, out_dir = crit_options$out_dir + ) return(res) } + diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index c5552fa2..9dc99f8b 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -7,6 +7,7 @@ #' @keywords internal #' wrap_graDiEnt <- function(optim_options, param_info, crit_options) { + if (is.null((ranseed <- optim_options$ranseed))) { ranseed <- NULL } @@ -43,6 +44,11 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { start_time <- Sys.time() + .trace_env <- new.env(parent = emptyenv()) + .trace_env$par_list <- list() + .trace_env$crit_list <- list() + .trace_env$k <- 0L + if (!is.null(ranseed)) set.seed(ranseed) range_bounds <- bounds$ub - bounds$lb ObjFun <- function(u) { @@ -50,7 +56,11 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { names(u) <- param_names x <- bounds$lb + u * range_bounds names(x) <- param_names - main_crit(x, crit_options = crit_options) + val <- main_crit(x, crit_options = crit_options) + .trace_env$k <- .trace_env$k + 1L + .trace_env$par_list[[.trace_env$k]] <- x + .trace_env$crit_list[[.trace_env$k]] <- val + return(val) } SQGDE <- tryCatch( graDiEnt::optim_SQGDE( @@ -120,39 +130,25 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { ) } trace_df <- NULL - if (!is.null(SQGDE$particles_trace) && length(dim(SQGDE$particles_trace)) == 3) { - tr <- SQGDE$particles_trace - n_it <- dim(tr)[1] - n_pop <- dim(tr)[2] - df_list <- vector("list", n_it) - eval_counter <- 0 - for (it in seq_len(n_it)) { - pop_u <- tr[it, , ] - pop_u <- as.matrix(pop_u) - colnames(pop_u) <- param_names + if (.trace_env$k > 0L) { + n_pop <- control_params$n_particles + k <- .trace_env$k + pars_mat <- do.call(rbind, lapply(.trace_env$par_list, function(v) as.numeric(v))) + colnames(pars_mat) <- param_names + crit_vec <- as.numeric(unlist(.trace_env$crit_list)) - # transformation - pop_x <- sweep(pop_u, 2, range_bounds, `*`) - pop_x <- sweep(pop_x, 2, bounds$lb, `+`) + n_it_est <- ceiling(k / n_pop) + iter <- rep(seq_len(n_it_est), each = n_pop)[seq_len(k)] + ind <- rep(seq_len(n_pop), times = n_it_est)[seq_len(k)] - df <- as.data.frame(pop_x) - - df$ind <- seq_len(n_pop) - df$iter <- it - crit_pop <- apply(pop_x, 1, function(par) { - names(par) <- param_names - main_crit(par, crit_options = crit_options) - }) - df$crit <- crit_pop - idx <- seq_len(n_pop) - df$eval <- eval_counter + idx - eval_counter <- eval_counter + n_pop - df$method <- "graDiEnt" - df$rep <- 1L - df_list[[it]] <- df - } - trace_df <- dplyr::bind_rows(df_list) - } + trace_df <- as.data.frame(pars_mat) + trace_df$ind <- ind + trace_df$iter <- iter + trace_df$crit <- crit_vec + trace_df$eval <- seq_len(k) + trace_df$method <- "graDiEnt" + trace_df$rep <- 1L + } res <- list( final_values = final_values, init_values = init_values, diff --git a/man/wrap_cmaes.Rd b/man/wrap_cmaes.Rd index deb10b0a..e038fdbf 100644 --- a/man/wrap_cmaes.Rd +++ b/man/wrap_cmaes.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/wrap_cmaes.R \name{wrap_cmaes} \alias{wrap_cmaes} -\title{A wrapper for parma::cmaes package (CMA-ES)} +\title{Wrap CMA-ES using parma} \usage{ wrap_cmaes(optim_options, param_info = NULL, crit_options) } @@ -20,6 +20,6 @@ that must be passed to main_crit function by the methods wrappers.} list with final_values, init_values, est_values, GE, min_crit_value } \description{ -A wrapper for parma::cmaes package (CMA-ES) +Wrap CMA-ES using parma } \keyword{internal} From 60853d9f54d83205c8d07fbe713d68dc0fdbe679 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Sat, 10 Jan 2026 20:47:10 +0000 Subject: [PATCH 15/16] Style code (GHA) --- R/wrap_DE.R | 18 ++++---- R/wrap_cmaes.R | 105 +++++++++++++++++++++++----------------------- R/wrap_graDiEnt.R | 9 ++-- 3 files changed, 65 insertions(+), 67 deletions(-) diff --git a/R/wrap_DE.R b/R/wrap_DE.R index acb34e31..a21d6482 100644 --- a/R/wrap_DE.R +++ b/R/wrap_DE.R @@ -54,7 +54,7 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { start_time <- Sys.time() .trace_env <- new.env(parent = emptyenv()) - .trace_env$x_list <- list() + .trace_env$x_list <- list() .trace_env$crit_list <- list() .trace_env$k <- 0L @@ -62,7 +62,7 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { fn_de <- function(x) { val <- main_crit(x, crit_options = crit_options) .trace_env$k <- .trace_env$k + 1L - .trace_env$x_list[[.trace_env$k]] <- as.numeric(x) + .trace_env$x_list[[.trace_env$k]] <- as.numeric(x) .trace_env$crit_list[[.trace_env$k]] <- val return(val) } @@ -103,26 +103,26 @@ wrap_DEoptim <- function(optim_options, param_info, crit_options) { trace_df <- NULL if (!is.null(DE$member$storepop)) { storepop <- DE$member$storepop - itermax <- length(storepop) - NP <- nrow(storepop[[1]]) + itermax <- length(storepop) + NP <- nrow(storepop[[1]]) pop_all <- do.call(rbind, lapply(seq_len(itermax), function(it) { m <- as.data.frame(storepop[[it]]) colnames(m) <- param_names - m$ind <- seq_len(NP) + m$ind <- seq_len(NP) m$iter <- it m - })) + })) crit_vec <- as.numeric(unlist(.trace_env$crit_list)) n <- min(nrow(pop_all), length(crit_vec)) pop_all <- pop_all[seq_len(n), , drop = FALSE] pop_all$crit <- crit_vec[seq_len(n)] - pop_all$eval <- seq_len(n) - pop_all$rep <- 1L + pop_all$eval <- seq_len(n) + pop_all$rep <- 1L pop_all$method <- "DEoptim" trace_df <- pop_all - } + } res <- list( final_values = final_values, diff --git a/R/wrap_cmaes.R b/R/wrap_cmaes.R index f80247c7..0fb278cf 100644 --- a/R/wrap_cmaes.R +++ b/R/wrap_cmaes.R @@ -10,7 +10,7 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { if (!requireNamespace("parma", quietly = TRUE)) { stop("the package 'parma' (not installed) becouse it required a dependencies that doesn't exist in the cluster") } - if(is.null((ranseed <- optim_options$ranseed))) { + if (is.null((ranseed <- optim_options$ranseed))) { ranseed <- NULL } n_params <- optim_options$n_params @@ -24,7 +24,6 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { algorithm <- "cmaes" if (is.null(param_info) || is.null(crit_options)) { - return(list( package = "parma", family = "Global", method = "cmaes", init_values_nb = init_nb <- if (!is.null(optim_options$ctrl$options$PopSize)) optim_options$ctrl$options$PopSize else NA_integer_ @@ -41,11 +40,9 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { ctrl <- parma::cmaes.control() - # if ctrl = list(options=..., CMA=...) existe if (!is.null(control_params) && is.list(control_params) && - (!is.null(control_params$options) || !is.null(control_params$CMA))) { - + (!is.null(control_params$options) || !is.null(control_params$CMA))) { if (!is.null(control_params$options)) { for (nm in names(control_params$options)) { ctrl$options[[nm]] <- control_params$options[[nm]] @@ -65,9 +62,9 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } } - if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){ + if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)) { crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter) - } + } start_time <- Sys.time() @@ -78,25 +75,25 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { trace_env$X <- list() keep_trace <- isTRUE(optim_options$return_trace) ObjFun <- function(x) { - x <- as.numeric(x) - names(x) <- param_names - val <- main_crit(x, crit_options = crit_options) - - if (keep_trace) { - trace_env$eval <- trace_env$eval + 1L - trace_env$crit[trace_env$eval] <- val - trace_env$X[[trace_env$eval]] <- x - } - return(val) + x <- as.numeric(x) + names(x) <- param_names + val <- main_crit(x, crit_options = crit_options) + + if (keep_trace) { + trace_env$eval <- trace_env$eval + 1L + trace_env$crit[trace_env$eval] <- val + trace_env$X[[trace_env$eval]] <- x } + return(val) + } insigma <- optim_options$insigma if (is.null(insigma)) { - # default value if it is not provided - insigma <- rep(0.3, nb_params) + # default value if it is not provided + insigma <- rep(0.3, nb_params) } if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) - #because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) + # because parma::cmaes expects a numeric vector for 'pars' (length = nb_params) x0 <- init_values # init_values can be a matrix/data.frame with several rows (nb_values x nb_params) our case data frame @@ -131,8 +128,8 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { } ) if (isTRUE(optim_options$debug)) { - message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) - } + message("parma_res fields: ", paste(names(parma_res), collapse = ", ")) + } elapsed <- Sys.time() - start_time @@ -143,25 +140,28 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { sol <- NULL if (!is.null(parma_res$bestever$x)) { - sol <- parma_res$bestever$x - } else if (!is.null(parma_res$par)) { - sol <- parma_res$par - } else if (!is.null(parma_res$xbest)) { - sol <- parma_res$xbest - } + sol <- parma_res$bestever$x + } else if (!is.null(parma_res$par)) { + sol <- parma_res$par + } else if (!is.null(parma_res$xbest)) { + sol <- parma_res$xbest + } if (is.null(sol)) { - stop("Impossible to find a solution in the object sent by parma::cmaes().") - } + stop("Impossible to find a solution in the object sent by parma::cmaes().") + } final_values <- as.numeric(sol) names(final_values) <- param_names # for post_treat_global_optim() - min_crit_value <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f - else main_crit(final_values, crit_options) + min_crit_value <- if (!is.null(parma_res$bestever$f)) { + parma_res$bestever$f + } else { + main_crit(final_values, crit_options) + } cmaes_bestever_f <- if (!is.null(parma_res$bestever$f)) parma_res$bestever$f else NA_real_ - cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ + cmaes_objective <- if (!is.null(parma_res$objective)) parma_res$objective else NA_real_ # Getting the estimated values @@ -178,22 +178,23 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { # Final solution trace_df <- NULL if (keep_trace && trace_env$eval > 0L) { - Xmat <- do.call(rbind, trace_env$X) - colnames(Xmat) <- param_names - - trace_df <- as.data.frame(Xmat) - trace_df$crit <- trace_env$crit - trace_df$eval <- seq_len(trace_env$eval) - lambda <- ctrl$options$PopSize - if (is.null(lambda) || is.na(lambda)) - lambda <- 30L - lambda <- as.integer(lambda) - - trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES - trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory - trace_df$rep <- 1L - trace_df$method <- "cmaes" - } + Xmat <- do.call(rbind, trace_env$X) + colnames(Xmat) <- param_names + + trace_df <- as.data.frame(Xmat) + trace_df$crit <- trace_env$crit + trace_df$eval <- seq_len(trace_env$eval) + lambda <- ctrl$options$PopSize + if (is.null(lambda) || is.na(lambda)) { + lambda <- 30L + } + lambda <- as.integer(lambda) + + trace_df$iter <- ((trace_df$eval - 1L) %/% lambda) + 1L # génération CMA-ES + trace_df$ind <- ((trace_df$eval - 1L) %% lambda) + 1L # individu dans la pop # single trajectory + trace_df$rep <- 1L + trace_df$method <- "cmaes" + } res <- list( final_values = final_values, @@ -203,13 +204,11 @@ wrap_cmaes <- function(optim_options, param_info = NULL, crit_options) { CMAES = parma_res, trace_df = trace_df, cmaes_bestever_f = cmaes_bestever_f, - cmaes_objective = cmaes_objective, + cmaes_objective = cmaes_objective, counteval = if (!is.null(parma_res$counteval)) parma_res$counteval else NA_integer_, - stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, + stopflag = if (!is.null(parma_res$stopflag)) parma_res$stopflag else NA_character_, out_dir = crit_options$out_dir - ) return(res) } - diff --git a/R/wrap_graDiEnt.R b/R/wrap_graDiEnt.R index 9dc99f8b..282652f4 100644 --- a/R/wrap_graDiEnt.R +++ b/R/wrap_graDiEnt.R @@ -7,7 +7,6 @@ #' @keywords internal #' wrap_graDiEnt <- function(optim_options, param_info, crit_options) { - if (is.null((ranseed <- optim_options$ranseed))) { ranseed <- NULL } @@ -45,7 +44,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { start_time <- Sys.time() .trace_env <- new.env(parent = emptyenv()) - .trace_env$par_list <- list() + .trace_env$par_list <- list() .trace_env$crit_list <- list() .trace_env$k <- 0L @@ -58,7 +57,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { names(x) <- param_names val <- main_crit(x, crit_options = crit_options) .trace_env$k <- .trace_env$k + 1L - .trace_env$par_list[[.trace_env$k]] <- x + .trace_env$par_list[[.trace_env$k]] <- x .trace_env$crit_list[[.trace_env$k]] <- val return(val) } @@ -139,7 +138,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { n_it_est <- ceiling(k / n_pop) iter <- rep(seq_len(n_it_est), each = n_pop)[seq_len(k)] - ind <- rep(seq_len(n_pop), times = n_it_est)[seq_len(k)] + ind <- rep(seq_len(n_pop), times = n_it_est)[seq_len(k)] trace_df <- as.data.frame(pars_mat) trace_df$ind <- ind @@ -148,7 +147,7 @@ wrap_graDiEnt <- function(optim_options, param_info, crit_options) { trace_df$eval <- seq_len(k) trace_df$method <- "graDiEnt" trace_df$rep <- 1L - } + } res <- list( final_values = final_values, init_values = init_values, From bddd5dfa3f6d2da3f12e64a7270675db851df173 Mon Sep 17 00:00:00 2001 From: nisrinemouhrim Date: Fri, 27 Feb 2026 11:34:39 +0100 Subject: [PATCH 16/16] Add rgenoud : global method --- NAMESPACE | 1 + R/optim_switch.R | 35 +++- R/wrap_parma.R | 136 ++++++++++++++ R/wrap_rgenoud.R | 218 ++++++++++++++++++++++ man/wrap_cmaes.Rd | 10 +- man/wrap_rgenoud.Rd | 37 ++++ tests/testthat/test_estim_param_DEoptim.R | 4 +- tests/testthat/test_estim_param_GEoptim.R | 3 +- tests/testthat/test_estim_param_rgenoud.R | 194 +++++++++++++++++++ 9 files changed, 625 insertions(+), 13 deletions(-) create mode 100644 R/wrap_parma.R create mode 100644 R/wrap_rgenoud.R create mode 100644 man/wrap_rgenoud.Rd create mode 100644 tests/testthat/test_estim_param_rgenoud.R diff --git a/NAMESPACE b/NAMESPACE index dff958e2..3ac4784a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -54,6 +54,7 @@ importFrom(graDiEnt,optim_SQGDE) importFrom(parma,cmaes) importFrom(parma,cmaes.control) importFrom(purrr,modify) +importFrom(rgenoud,genoud) importFrom(rlang,.data) importFrom(stats,setNames) importFrom(tibble,tibble) diff --git a/R/optim_switch.R b/R/optim_switch.R index 1f3e0e9b..dd05c51b 100644 --- a/R/optim_switch.R +++ b/R/optim_switch.R @@ -214,18 +214,35 @@ optim_switch <- function(...) { out_dir = crit_options$out_dir ) } - } else if (optim_method == "cmaes") { - if (is.null(optim_options$n_particles) || is.na(optim_options$n_particles)) { - if (!is.null(optim_options$ctrl$options$PopSize) && - !is.na(optim_options$ctrl$options$PopSize)) { - optim_options$n_particles <- as.integer(optim_options$ctrl$options$PopSize) - } else { - optim_options$n_particles <- 30L + }else if (optim_method == "genoud") { + res <- do.call(wrap_rgenoud, wrap_args) + if (nargs() > 2) { + res$obs_var_list <- .croptEnv$obs_var_list + res$obs_situation_list <- .croptEnv$obs_situation_list + if (arguments$crit_options$info_level >= 1) { + res$params_and_crit <- dplyr::bind_rows(.croptEnv$params_and_crit) } + res <- post_treat_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + crit_options = crit_options + ) + res$plots <- plot_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) + summary_global_optim( + optim_options = optim_options, + param_info = param_info, + optim_results = res, + out_dir = crit_options$out_dir + ) } - + } else if (optim_method == "cmaes") { wrap_args$optim_options <- optim_options - res <- do.call(wrap_cmaes, wrap_args) if (nargs() > 2) { diff --git a/R/wrap_parma.R b/R/wrap_parma.R new file mode 100644 index 00000000..43920b06 --- /dev/null +++ b/R/wrap_parma.R @@ -0,0 +1,136 @@ +#' @title A wrapper for parma::cmaes package (CMA-ES) +#' +#' @inheritParams optim_switch +#' +#' @return list with final_values, init_values, est_values, GE, min_crit_value +#' @importFrom parma cmaes cmaes.control +#' @keywords internal +#' +wrap_cmaes <- function(optim_options, param_info, crit_options) { + if(is.null((ranseed <- optim_options$ranseed))) { + ranseed = NULL + } + if (is.null((n_params <- optim_options$n_params))) { + stop( + "Optim_options$n_params is mandatory, please define it (e.g., 10 * number of parameters)" + ) + } + if (is.null((trace <- optim_options$return_trace))) { + optim_options$return_trace <- TRUE + } + optim_options$ranseed <- NULL # ranseed is set to NULL because getAlgoParams doesn't regonise it + algorithm <- "parma_cmaes" + + # return requested information if only optim_options is given in argument + if (nargs() == 1 && methods::hasArg(optim_options)) { + init_nb <- if (!is.null(optim_options$n_particles)) { + optim_options$n_particles + } else { + NA_integer_ + } + return(list( + package = "parma", family = "Global", + method = "cmaes", init_values_nb = init_nb + )) + } + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + bounds <- get_params_bounds(param_info) + init_values <- get_init_values(param_info) + + control_params <- optim_options$ctrl + if (is.null(control_params )) { + control_params <- parma::cmaes.control() + } else{ + if (is.null(control_params$options) && is.null(control_params$CMA)) { + ctrl <- parma::cmaes.control(options = control_params, CMA = list()) + }else { + ctrl <- do.call(parma::cmaes.control, control_params) + } + } + if (!is.null(ctrl$options$PopSize) && !is.null(ctrl$options$MaxIter)){ + crit_options$tot_max_eval <- as.integer(ctrl$options$PopSize * ctrl$options$MaxIter) + } + + start_time <- Sys.time() + + if (!is.null(ranseed)) set.seed(ranseed) + ObjFun <- function(x){ + x <- as.numeric(x) + names(x) <- param_names + main_crit(x, crit_options = crit_options) + } + insigma <- optim_options$insigma + if (is.null(insigma)) { + # default value if it is not provided + insigma <- rep(0.3, nb_params) + } + if (length(insigma) == 1L) insigma <- rep(insigma, nb_params) + parma_res <- tryCatch( + parma::cmaes( + pars = init_values, + fun = ObjFun, + lower = bounds$lb, + upper = bounds$ub, + insigma = insigma, + ctrl = ctrl + ), + error = function(e) { + warning(sprintf("parma::cmaes failed: %s", e$message)) + NULL + } + ) + elapsed <- Sys.time() - start_time + + if (is.null(parma_res)) { + stop("parma::cmaes failed (object NULL).") + } + + sol <- NULL + if (!is.null(parma_res$par)) sol <- parma_res$par + if (is.null(sol) && !is.null(parma_res$bestever$x)) + sol <- parma_res$bestever$x + if (is.null(sol) && !is.null(parma_res$xbest)) + sol <- parma_res$xbest + if (is.null(sol)) { + stop("Impossible to find a solution in the object sent by parma::cmaes().") + } + + final_values <- as.numeric(sol) + names(final_values) <- param_names + + min_crit_value <- NULL + if (!is.null(parma_res$value)) + min_crit_value <- parma_res$value + if (is.null(min_crit_value) && !is.null(parma_res$bestever$f)) + min_crit_value <- parma_res$bestever$f + if (is.null(min_crit_value)) { + min_crit_value <- main_crit(final_values, crit_options = crit_options) + } + + # Get the estimated values ==> matrix of 1 * nb_params + est_values <- matrix( + final_values, + nrow = 1L, + dimnames = list(NULL, param_names) + ) + colnames(est_values) <- param_names + if (ncol(est_values) == nb_params) { + colnames(est_values) <- param_names + } + + # Final solution + trace_df <- NULL + + res <- list( + final_values = final_values, + init_values = init_values, + est_values = est_values, + min_crit_value = min_crit_value, + CMAES = parma_res, + trace_df = trace_df + ) + + return(res) +} + diff --git a/R/wrap_rgenoud.R b/R/wrap_rgenoud.R new file mode 100644 index 00000000..c01efda1 --- /dev/null +++ b/R/wrap_rgenoud.R @@ -0,0 +1,218 @@ +#' @title A wrapper for rgenoud package (single criterion) +#' @inheritParams optim_switch +#' @importFrom rgenoud genoud +#' @return A list containing: +#' `final_values`: vector of estimated parameter values +#' corresponding to the best solution. +#' `init_values`: vector or matrix of initial parameter values. +#' `est_values`: matrix of estimated values (final solution). +#' `min_crit_value`: minimum value of the criterion. +#' `trace_df`: data.frame containing parameters, criterion, +#' iteration, evaluation and individual indices . +#' `obs_var_list`: list of observed variables used by the criterion. +#' `elapsed`: computation time. +#' `genoud.pro`: an output file called located in the temporary +#' directory provided by tempdir, containing solutions in each +#' generation (criterion and parameter values ) +#' +#' @keywords internal + +wrap_rgenoud <- function(optim_options, param_info, crit_options) { + if (nargs() == 1 && methods::hasArg(optim_options)) { + init_TP <- if (!is.null(optim_options$pop.size)) + optim_options$pop.size + else + NA_integer_ + return(list( + package = "rgenoud", family = "Global", + method = "genoud", init_values_nb = init_TP + )) + } + if (is.null(optim_options$pop.size)) { + stop("optim_options$pop.size is mandatory, please define it (e.g., 10 * number of parameters)") + } + + ranseed <- optim_options$ranseed + optim_options$ranseed <- NULL + + param_names <- get_params_names(param_info) + nb_params <- length(param_names) + + bounds <- get_params_bounds(param_info) + init_values <- get_init_values(param_info) + + Domains <- cbind(bounds$lb, bounds$ub) + + # read final population from genoud.pro to construct est_value + read_last_population_from_pro <- function(project_path, nb_params, param_names) { + if (is.null(project_path) || !file.exists(project_path)) + return(NULL) + L <- readLines(project_path, warn = FALSE) + gen_idx <- grep("^Generation:\\s*[0-9]+", L) + if (length(gen_idx) == 0) + return(NULL) + start <- gen_idx[length(gen_idx)] # last generation + block <- L[start:length(L)] + rows <- grep("^\\s*[0-9]+\\s+", block, value = TRUE) + if (length(rows) == 0) + return(NULL) + tok <- strsplit(trimws(rows), "\\s+") + mat <- do.call(rbind, lapply(tok, function(x) as.numeric(x))) + if (ncol(mat) < (2 + nb_params)) #(je vérifie si ce n est pas 3+nb_params) + return(NULL) + crit <- mat[,2] + pars <- mat[, 3:(2+nb_params), drop = FALSE] # ici aussi je dois vérifier!! + colnames(pars) <- param_names + list(pars = pars, crit = crit) + } + + starting.values <- init_values + if (is.vector(starting.values)) { + starting.values <- matrix(starting.values, nrow = 1) + } + if (is.data.frame(starting.values)) { + starting.values <- as.matrix(starting.values) + } + if (!is.null(ranseed)) { + set.seed(ranseed) + if (is.null(optim_options$unif.seed)) + optim_options$unif.seed <- sample.int(.Machine$integer.max, 1) + if (is.null(optim_options$int.seed)) + optim_options$int.seed <- sample.int(.Machine$integer.max, 1) + } + # project file + if (is.null(optim_options$project.path)) { + optim_options$project.path <- file.path( + crit_options$out_dir, + "genoud.pro" + ) + } + # tot_max_eval : genoud could do more (local search!!) + crit_options$tot_max_eval <- optim_options$pop.size * optim_options$max.generations + + # trace env + .trace_env <- new.env(parent = emptyenv()) + .trace_env$x_list <- list() + .trace_env$crit_list <- list() + .trace_env$k <- 0 + + fn_gen <- function(x) { + val <- main_crit(x, crit_options = crit_options) + .trace_env$k <- .trace_env$k + 1 + .trace_env$x_list[[.trace_env$k]] <- as.numeric(x) + .trace_env$crit_list[[.trace_env$k]] <- val + val + } + + start_time <- Sys.time() + + GENOUD <- tryCatch( + do.call( + rgenoud::genoud, + c( + list( + fn = fn_gen, + nvars = nb_params, + max = FALSE, + Domains = Domains, + starting.values = starting.values + ), + optim_options + ) + ), + error = function(e) { + warning(sprintf("genoud failed: %s", e$message)) + NULL + } + ) + + elapsed <- Sys.time() - start_time + + if (is.null(GENOUD)) { + return(NULL) + } + + # final solution + final_values <- as.numeric(GENOUD$par) + names(final_values) <- param_names + min_crit_value <- GENOUD$value + + # est_values : the best individual / final population en amont !! + est_values <- matrix(final_values, nrow = 1, dimnames = list(NULL, param_names)) + # trace_df reconstruction from saved file + trace_df <- NULL + if (.trace_env$k > 0L) { + X <- do.call(rbind, .trace_env$x_list) # k x nb_params + critv <- as.numeric(unlist(.trace_env$crit_list)) + eval <- seq_along(critv) + popsize_used <- optim_options$pop.size + + + # iter = batch (pas génération genoud ), ind = index dans le batch + iter <- ((eval - 1L) %/% popsize_used) + 1L #à ver + ind <- ((eval - 1L) %% popsize_used) + 1L + + colnames(X) <- param_names + trace_df <- data.frame( + iter = as.integer(iter), + ind = as.integer(ind), + as.data.frame(X), + crit = critv, + eval = as.integer(eval), + rep = 1L, + method = "rgenoud" + ) + # est_values from genoud.pro (final gener) + pro_pop <- read_last_population_from_pro( + project_path = optim_options$project.path, + nb_params = nb_params, + param_names = param_names + ) + if (!is.null(pro_pop)) { + est_values <- as.matrix(pro_pop$pars) + } + } + #obs_var_list + info_final <- tryCatch( + main_crit( + param_values = final_values, + crit_options = c(crit_options, return_detailed_info = TRUE) + ), + error = function(e) NULL + ) + + obs_var_list <- NULL + if (!is.null(info_final)) { + if (!is.null(info_final$obs_var_list)) + obs_var_list <- info_final$obs_var_list + else if (!is.null(info_final$obs_intersect)) + obs_var_list <- names(info_final$obs_intersect) + } + # Stop info + stop_reason <- "HARD MAXIMUM GENERATION LIMIT HIT" + if (!is.null(GENOUD$hard.generation.limit) && isTRUE(GENOUD$hard.generation.limit)) { + stop_reason <- "hard maximum generation limit hit" + } + if (!is.null(GENOUD$message) && nzchar(GENOUD$message)) { + stop_reason <- GENOUD$message + } + generations_run <- NA_integer_ + if (!is.null(GENOUD$generations)) + generations_run <- as.integer(GENOUD$generations) + + res <- list( + final_values = final_values, + init_values = init_values, + est_values = est_values, + min_crit_value = min_crit_value, + GENOUD = GENOUD, + trace_df = trace_df, + obs_var_list = obs_var_list, + elapsed = elapsed, + stop_reason = stop_reason, + generations_run = generations_run, + project_path = optim_options$project.path + ) + + return(res) +} diff --git a/man/wrap_cmaes.Rd b/man/wrap_cmaes.Rd index e038fdbf..6e8cb8a2 100644 --- a/man/wrap_cmaes.Rd +++ b/man/wrap_cmaes.Rd @@ -1,10 +1,12 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/wrap_cmaes.R +% Please edit documentation in R/wrap_cmaes.R, R/wrap_parma.R \name{wrap_cmaes} \alias{wrap_cmaes} \title{Wrap CMA-ES using parma} \usage{ -wrap_cmaes(optim_options, param_info = NULL, crit_options) +wrap_cmaes(optim_options, param_info, crit_options) + +wrap_cmaes(optim_options, param_info, crit_options) } \arguments{ \item{optim_options}{see description in estim_param} @@ -17,9 +19,13 @@ function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model that must be passed to main_crit function by the methods wrappers.} } \value{ +list with final_values, init_values, est_values, GE, min_crit_value + list with final_values, init_values, est_values, GE, min_crit_value } \description{ Wrap CMA-ES using parma + +A wrapper for parma::cmaes package (CMA-ES) } \keyword{internal} diff --git a/man/wrap_rgenoud.Rd b/man/wrap_rgenoud.Rd new file mode 100644 index 00000000..46329331 --- /dev/null +++ b/man/wrap_rgenoud.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/wrap_rgenoud.R +\name{wrap_rgenoud} +\alias{wrap_rgenoud} +\title{A wrapper for rgenoud package (single criterion)} +\usage{ +wrap_rgenoud(optim_options, param_info, crit_options) +} +\arguments{ +\item{optim_options}{see description in estim_param} + +\item{param_info}{see description in estim_param} + +\item{crit_options}{List containing several arguments given to \code{estim_param} +function: \code{param_names}, \code{obs_list}, \code{crit_function}, \code{model_function}, +\code{model_options}, \code{param_info}, \code{transform_obs}, \code{transform_sim}, \code{out_dir} +that must be passed to main_crit function by the methods wrappers.} +} +\value{ +A list containing: +\code{final_values}: vector of estimated parameter values +corresponding to the best solution. +\code{init_values}: vector or matrix of initial parameter values. +\code{est_values}: matrix of estimated values (final solution). +\code{min_crit_value}: minimum value of the criterion. +\code{trace_df}: data.frame containing parameters, criterion, +iteration, evaluation and individual indices . +\code{obs_var_list}: list of observed variables used by the criterion. +\code{elapsed}: computation time. +\code{genoud.pro}: an output file called located in the temporary +directory provided by tempdir, containing solutions in each +generation (criterion and parameter values ) +} +\description{ +A wrapper for rgenoud package (single criterion) +} +\keyword{internal} diff --git a/tests/testthat/test_estim_param_DEoptim.R b/tests/testthat/test_estim_param_DEoptim.R index 32e13ce9..a653cdf3 100644 --- a/tests/testthat/test_estim_param_DEoptim.R +++ b/tests/testthat/test_estim_param_DEoptim.R @@ -143,7 +143,9 @@ test_that("estim_param 1 step OLS criterion", { optim_options <- list( ranseed = 1234, - itermax = 60, + itermax = 50, + reltol = 3.929623e-04, + steptol = 4, NP = 30 ) diff --git a/tests/testthat/test_estim_param_GEoptim.R b/tests/testthat/test_estim_param_GEoptim.R index 09c88a80..b0db6540 100644 --- a/tests/testthat/test_estim_param_GEoptim.R +++ b/tests/testthat/test_estim_param_GEoptim.R @@ -145,9 +145,10 @@ test_that("estim_param 1 step OLS criterion", { optim_options <- list( ranseed = 1234, n_params = length(param_info), - n_iter = 60, + n_iter = 50, n_particles = 30, n_diff = 2, + stop_check = 3, return_trace = TRUE ) diff --git a/tests/testthat/test_estim_param_rgenoud.R b/tests/testthat/test_estim_param_rgenoud.R new file mode 100644 index 00000000..86359997 --- /dev/null +++ b/tests/testthat/test_estim_param_rgenoud.R @@ -0,0 +1,194 @@ +context("Test the estim_param function using a toy model") +library(tibble) + +# Define a toy model and its wrapper +toymodel <- function(nb_days, rB = 0.1, Bmin = 0.1, Bmax = 10, h = 0.5, Bini = 0.01) { + # Simulate biomass and yield with a simple logistic model + # + # Arguments + # - rB: growth rate + # - Bmin: minimum biomass + # - Bmax: maximum biomass + # - h: harvest rate + # - Bini: initial biomass + # + # Value + # A list containing the simulated biomass and yield + # + biom <- numeric(nb_days) + biom[1] <- Bini + for (i in 2:nb_days) { + biom[i] <- biom[i - 1] + rB * (biom[i - 1] + Bmin) * (1 - biom[i - 1] / Bmax) + } + yield <- h * biom + return(list(biom = biom, yield = yield)) +} +clip <- function(x, lb, ub) pmin(pmax(x, lb), ub) + +toymodel_wrapper <- function(param_values = NULL, situation, + model_options, ...) { + # A wrapper for toymodel + # + # Arguments + # - param_values: (optional) a named vector or a tibble containing the + # values of the toymodel parameters to force. + # - situation: Vector of situations names for which results must be + # returned. + # In this case, the names of the situations are coded as "year_suffix" + # - model_options: a tibble containing the begin and end date of simulation + # for each situation. + # - ...: mandatory since CroptimizR will give additional arguments not used + # here + # + # Value: + # A named list of tibble per situation. + # Each tibble contains columns: + # - Date (POSIXct dates of simulated results), + # - One column per simulated variable (biomass and yield) + # + results <- list( + sim_list = setNames( + vector("list", length(situation)), + nm = situation + ), + error = FALSE + ) + attr(results$sim_list, "class") <- "cropr_simulation" + + # Remove dummy in param_values, just here for specific tests + if (!is.null(param_values) && "dummy" %in% names(param_values)) { + param_values <- param_values[!names(param_values) %in% "dummy"] + } + + for (sit in situation) { + # Retrieve begin and end date from situation name + begin_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(begin_date) + end_date <- dplyr::filter(model_options, situation == sit) %>% + dplyr::select(end_date) + + date_sequence <- seq(from = begin_date[[1]], to = end_date[[1]], by = "day") + + if (!all(names(param_values) %in% c( + "rB", "Bmin", "Bmax", "h", "Bini" + ))) { + warning( + paste( + "Unknown parameters in param_values:", + paste(names(param_values), collapse = ",") + ) + ) + results$error <- TRUE + return(results) + } + + # --- force parameters into bounds (important for Nelder-Mead) --- + if (!is.null(param_values)) { + if ("rB" %in% names(param_values)) param_values["rB"] <- clip(param_values["rB"], 0, 1) + if ("h" %in% names(param_values)) param_values["h"] <- clip(param_values["h"], 0, 1) + if ("Bmax" %in% names(param_values)) param_values["Bmax"] <- clip(param_values["Bmax"], 5, 15) + } + # --------------------------------------------------------------- + + res <- do.call( + "toymodel", + c(nb_days = length(date_sequence), as.list(param_values)) + ) + + # Fill the result variable + results$sim_list[[sit]] <- + dplyr::tibble( + Date = date_sequence, + biomass = res$biom, + yield = res$yield + ) + } + return(results) +} + +# To make these function accessible to the test environment +.GlobalEnv$toymodel <- toymodel +.GlobalEnv$toymodel_wrapper <- toymodel_wrapper + +# Use setup() to define shared data for all tests in this file +setup({ + # These variables will be available in all tests + .GlobalEnv$model_options <- tibble::tibble( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + begin_date = as.Date(c("2000-05-01", "2001-05-12", "2003-05-05", "2004-05-15")), + end_date = as.Date(c("2000-11-05", "2001-11-20", "2003-11-15", "2004-11-18")) + ) + .GlobalEnv$param_true_values <- c(rB = 0.08, h = 0.55, Bmax = 7) + + # Generate synthetic observations + tmp <- toymodel_wrapper( + situation = c("sit1_2000", "sit1_2001", "sit2_2003", "sit2_2004"), + param_values = param_true_values, + model_options = model_options + ) + + .GlobalEnv$obs_synth <- lapply(tmp$sim_list, function(x) { + return(dplyr::filter(x, dplyr::row_number() %% 10 == 0)) + }) +}) + +# For cleaning up after tests +teardown({ + if (exists("model_options", envir = .GlobalEnv)) rm(model_options, envir = .GlobalEnv) + if (exists("param_true_values", envir = .GlobalEnv)) rm(param_true_values, envir = .GlobalEnv) + if (exists("obs_synth", envir = .GlobalEnv)) rm(obs_synth, envir = .GlobalEnv) +}) + + +# Test estim_param with DEoptim +test_that("estim_param 1 step OLS criterion", { + param_info <- list( + rB = list(lb = 0, ub = 1, init_values = 0.2), + h = list(lb = 0, ub = 1, init_values = 0.2), + Bmax = list(lb = 5, ub = 15, init_values = 6) + ) + + + optim_options <- list( + ranseed = 1234, + pop.size = 30, + max.generations = 20, + wait.generations = 30, + BFGS = TRUE, + optim.method = "Nelder-Mead", + boundary.enforcement = 1, + gradient.check = FALSE, + control = list(maxit = 1, reltol = 1e-06), + solution.tolerance = 1e-3, + hard.generation.limit = FALSE, + print.level = 3, + P9 = 0 + ) + + res <- estim_param( + obs_list = obs_synth, + model_function = toymodel_wrapper, + model_options = model_options, + optim_method = "genoud", + crit_function = crit_ols, + optim_options = optim_options, + param_info = param_info, + obs_var = c("biomass", "yield"), + situation = c("sit1_2000", "sit1_2001", "sit2_2003"), + out_dir = tempdir() + ) + str(res$trace_df) + + expect_equal(res$final_values[["rB"]], + param_true_values[["rB"]], + tolerance = param_true_values[["rB"]] * 1e-2 + ) + expect_equal(res$final_values[["h"]], + param_true_values[["h"]], + tolerance = param_true_values[["h"]] * 1e-2 + ) + expect_equal(res$final_values[["Bmax"]], + param_true_values[["Bmax"]], + tolerance = param_true_values[["Bmax"]] * 1e-2 + ) +})