R/funcs_run-optimization.R

In cath-parkinson/mm.reoptimise: Run and Compare Budget Optimisation in Different Scenarios

#' Run optimization
#'
#' This function takes a prepared scenario object, and runs the optimization algorithm,
#' to get the final results.The results are appeneded to the scenario object,
#' plus 4 key columns are populated in curves_full \code{optim_spend},
#' \code{response_at_optim}, \code{same_units} and \code{net_response_at_optim}
#'
#' @param scenario A named list representing a scenario list-object.
#' @param constant Possibility to change the constant used.
#'
#' @return An updated scenario list-object.
#'
#' @export
run_optimization <- function(scenario, constant = NULL){

   # stop errors in R package check
   param2 <- NULL

   # Auto detect which type of optimization is needed
   check_for_net <- get_scenario_budget(scenario) %>% is.na() # if NA then apply "net"
   check_for_pct <- scenario$curves_filtered %>% dplyr::pull(param2) %>% max() # if less than one apply "alloc_pct"

   if(check_for_net){ optim_type <- "net" } else
   if(check_for_pct < 1){ optim_type <- "alloc_pct" } else { optim_type <- "alloc" }

   # Default constants based on optim_type
   # If constant was supplied by user, we use that one
   if (is.null(constant)){
      if(optim_type == "alloc")     constant = 1e-6
      if(optim_type == "alloc_pct") constant = 1e10
      if(optim_type == "net")       constant = 1e-6
   }

   # Create data to be used in the optimisation
   scenario$data_optim <- get_data_optim(scenario)

   # Create constraints strings
   scenario$constraints_strings <- get_constraints_strings(scenario,
                                                           optim_type,
                                                           constant)

   # Run optimisation based on optim_type
   scenario <- get_optim_results(scenario, optim_type, constant)

   return(scenario)
}


get_optim_results <- function(scenario, optim_type, constant){

      data_optim <- scenario$data_optim

      .x0 <- rep(0, nrow(scenario$data_optim))
      .opts <- get_optim_opts(scenario, optim_type)
      .eval_f <- get_eval_f(optim_type)
      .eval_const <- get_eval_const(scenario, optim_type)

      results <- nloptr::nloptr(
            x0 = .x0,
            eval_f = .eval_f,
            opts = .opts,
            lb = rep(0, nrow(data_optim)),
            eval_g_ineq = .eval_const$.eval_g_ineq,
            eval_jac_g_ineq = .eval_const$.eval_jac_g_ineq,
            eval_g_eq = .eval_const$.eval_g_eq,
            eval_jac_g_eq = .eval_const$.eval_jac_g_eq,
            curve_equations = data_optim$equation,
            constant = constant,
            p1 = data_optim$param1,
            p2 = data_optim$param2,
            p3 = data_optim$param3,
            p4 = data_optim$param4,
            p5 = data_optim$param5,
            p6 = data_optim$param6,
            p7 = data_optim$param7,
            p8 = data_optim$param8,
            p9 = data_optim$param9
      )

      if (optim_type == "alloc_pct") constant = 1

      results$optim_spend <- tibble::tibble(
            alloc.unit_id = data_optim$alloc.unit_id,
            optim_spend = results$solution / constant
      )

      scenario$optim_results <- results

      scenario$curves_full <- add_optim_results(scenario$curves_full,
                                                scenario$optim_results$optim_spend)

      return(scenario)
}