From 17d98c4fb8d388b9bbd34f162f950ffca14a7c03 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Wed, 5 Aug 2026 17:24:00 +0000 Subject: [PATCH 1/3] Make the restart selection callable without a solver, and fix two gaps The numerical optimizer's choice among its restarts was written inline inside a function nested in optimize_cost_nlcoptim(). Two of its steps, projecting the chosen restart onto the intervention bounds and recomputing its cost there, only run when every restart has left the box, which no two-component problem does. The regression suite had to reach them through a three-component configuration whose eleven restarts each miss the bounds by a unit or so in the last place, so a solver converging more tightly would have stopped the fallback firing and those tests would have kept passing while covering nothing. That could not be guarded from the outside: the projection leaves a component exactly on a bound and so does a solver that converged there, so the returned value cannot distinguish them. The filter, the fallback, the choice, the projection and the cost recomputation are one decision, so they move together into select_restart_within_bounds(), which takes the restart results and the bounds and needs no solver, no model and no data. The body is unchanged. Its unit tests cover all restarts in the box, some, none, a tie in cost, a failed restart among good ones, every restart failed, a single restart, and a winner genuinely outside the bounds. Removing the projection now fails thirteen of them, where before it was caught only by the fragile fixture, which is kept for the integration path. Two smaller gaps found while reviewing the same code: rec_int_processor() passed character(0) as the set of coefficient names already claimed by other terms, so at that call site the anchored-name fallback excluded nothing and the guard against a covariate being read as a fixed effect could not fire. It now builds the list the same way get_confidence_set() does, which meant threading additional_covariates through from lago_optimization(). The fallback turns out not to be reachable from any model the package fits, since the term mapping is only unavailable for an intercept-only fit, but the argument was wrong and a mutation to it survived the whole suite. Center weights within a thousandth of summing to 1 were accepted and then used as they were. They multiply the per-center outcomes and are summed, so weights summing to 0.999 scaled every reported outcome by 0.999, including the outcome the recommendation is chosen against. They are now renormalised where they are validated. The tolerance is unchanged, since it is documented and states that the input was meant to be a set of weights; renormalising is exact for a caller whose weights already sum to 1, and silent, because warning would fire on rounded input the documentation invites. The suite goes from 629 expectations to 712. --- R/get_recommended_interventions.R | 132 +++-- R/lago_optimization.R | 1 + R/rec_int_processor.R | 23 +- R/validate_inputs.R | 33 +- man/select_restart_within_bounds.Rd | 67 +++ tests/testthat/test-minimize-and-bounds.R | 41 +- tests/testthat/test-outcome-internals.R | 639 ++++++++++++++++++++++ 7 files changed, 881 insertions(+), 55 deletions(-) create mode 100644 man/select_restart_within_bounds.Rd diff --git a/R/get_recommended_interventions.R b/R/get_recommended_interventions.R index ae555aa..67e9ceb 100644 --- a/R/get_recommended_interventions.R +++ b/R/get_recommended_interventions.R @@ -30,6 +30,84 @@ unachievable_goal_message <- function(lower_outcome_goal) { } +#' select_restart_within_bounds +#' +#' @description Internal function that picks the recommended intervention out +#' of the points a multi-start numerical optimizer converged to, and returns it +#' together with the cost of the point being recommended. +#' +#' @details The filter, the fallback, the choice, the projection and the cost +#' recomputation are one decision and are deliberately kept together: each +#' step's correctness depends on the one before it. +#' +#' costs holds the cost each restart converged to, so the best restart is the +#' cheapest one. Every restart satisfies the outcome constraint, which +#' NlcOptim::solnl() enforces through confun, so cost is the only thing to +#' choose on among the ones that are actually implementable. +#' +#' solnl() treats the box as a soft constraint and will step a little outside it +#' to buy a lower objective, so the cheapest restart is systematically the one +#' furthest outside the bounds: selecting on cost alone selects for the +#' violation. Restarts that left the box are therefore dropped before the +#' comparison. If every restart left it they are all kept, so a recommendation +#' is still returned, and the survivor is brought back onto the box. +#' +#' The chosen restart can still sit a solver tolerance outside the box, when +#' every restart did. A recommendation has to be implementable, so it is +#' projected onto the bounds and its cost recomputed at the value actually being +#' recommended rather than reported from the point the solver stopped at. The +#' projection can only raise the cost of a monotone cost function, so reporting +#' the solver's cost understates what the recommendation costs. +#' +#' @param restart_points A numeric matrix with one column per restart and one +#' row per intervention component, holding the point each restart converged to. +#' @param costs A numeric vector, one entry per column of restart_points, with +#' the cost that restart converged to. NA marks a restart whose optimization +#' failed. Not every entry may be NA: the caller refuses that case with its own +#' message before calling this. +#' @param lower_bounds A numeric vector. The lower bounds of the intervention +#' components. +#' @param upper_bounds A numeric vector. The upper bounds of the intervention +#' components. +#' @param cost_fun A function of one numeric vector returning the total cost of +#' that intervention. +#' +#' @return A list with: +#' - int_components: the chosen intervention, projected onto the bounds. +#' - rec_int_cost: cost_fun() evaluated at int_components. +#' +#' @keywords internal +select_restart_within_bounds <- function(restart_points, + costs, + lower_bounds, + upper_bounds, + cost_fun) { + in_box <- apply( + restart_points, 2, + function(x) { + all(x >= lower_bounds) && + all(x <= upper_bounds) + } + ) + valid_indices <- which(!is.na(costs) & in_box) + if (length(valid_indices) == 0) { + valid_indices <- which(!is.na(costs)) + } + min_position <- valid_indices[which.min(costs[valid_indices])] + int_components <- restart_points[, min_position] + + int_components <- pmin( + pmax(int_components, lower_bounds), + upper_bounds + ) + + list( + int_components = int_components, + rec_int_cost = cost_fun(int_components) + ) +} + + #' get_recommended_interventions #' #' @description Internal function that calculates the LAGO recommended @@ -585,44 +663,26 @@ get_recommended_interventions <- function( )) } - # cost_results holds the cost each restart converged to, so the best - # restart is the cheapest one. Every restart satisfies the outcome - # constraint, which solnl() enforces through confun, so cost is the only - # thing to choose on among the ones that are actually implementable. - # - # solnl() treats the box as a soft constraint and will step a little - # outside it to buy a lower objective, so the cheapest restart is - # systematically the one furthest outside the bounds: selecting on cost - # alone selects for the violation. Restarts that left the box are - # therefore dropped before the comparison. If every restart left it they - # are all kept, so a recommendation is still returned, and the result is - # brought back onto the box below. - in_box <- apply( - results_int_components, 2, - function(x) { - all(x >= intervention_lower_bounds) && - all(x <= intervention_upper_bounds) - } - ) - valid_indices <- which(!is.na(cost_results) & in_box) - if (length(valid_indices) == 0) { - valid_indices <- which(!is.na(cost_results)) - } - min_position <- valid_indices[which.min(cost_results[valid_indices])] - int_components <- results_int_components[, min_position] - - # The chosen restart can still sit a solver tolerance outside the box, - # when every restart did. A recommendation has to be implementable, so - # it is brought back onto the bounds and its cost recomputed at the - # value actually being recommended rather than reported from the point - # the solver stopped at. - int_components <- pmin( - pmax(int_components, intervention_lower_bounds), - intervention_upper_bounds + # Choosing among the restarts, and making the winner implementable, is + # one decision and lives in select_restart_within_bounds(): the in-box + # filter, the keep-everything fallback, the cheapest-survivor choice, + # the projection onto the bounds and the cost recomputation at the + # projected point. It is a package-level internal rather than inline + # here so that each of those steps can be tested with a hand-built + # restart matrix, with no solver, model or data in between. The last two + # steps only do anything when EVERY restart left the box, which is + # exceptional from the outside, so they are not otherwise reachable from + # a test. + selected <- select_restart_within_bounds( + restart_points = results_int_components, + costs = cost_results, + lower_bounds = intervention_lower_bounds, + upper_bounds = intervention_upper_bounds, + cost_fun = cost_obj_fun ) - rec_int_cost <- cost_obj_fun(int_components) + rec_int_cost <- selected$rec_int_cost - est_rec_int <- int_components + est_rec_int <- selected$int_components est_reachable_outcome <- get_est_reachable_outcome( x = est_rec_int, diff --git a/R/lago_optimization.R b/R/lago_optimization.R index 415208d..4a13b40 100644 --- a/R/lago_optimization.R +++ b/R/lago_optimization.R @@ -393,6 +393,7 @@ lago_optimization <- function( data = data, model = if (lower_outcome_goal) new_model else model, center_characteristics = center_characteristics, + additional_covariates = additional_covariates, include_center_effects = include_center_effects, include_time_effects = include_time_effects, include_interaction_terms = include_interaction_terms, diff --git a/R/rec_int_processor.R b/R/rec_int_processor.R index a859cdf..d8af126 100644 --- a/R/rec_int_processor.R +++ b/R/rec_int_processor.R @@ -2,6 +2,7 @@ rec_int_processor <- function( data, model, center_characteristics, + additional_covariates = NULL, include_center_effects, include_time_effects, include_interaction_terms, @@ -83,9 +84,27 @@ rec_int_processor <- function( # coefficient that is not a center effect at all. # The fallback for a model whose term mapping could not be rebuilt, which # outcome_model_fitting() does not produce, is the name search this used to - # do, anchored so that only a name beginning with the term is matched. + # do, anchored so that only a name beginning with the term is matched, and + # restricted to the coefficients no other block below claims for itself. + # + # That exclusion list is the point of named_predictors, and it has to be the + # names THIS function looks up on its own account, exactly as + # get_confidence_set() passes its own: the intervention components, the + # additional covariates, the center characteristics and the intercept. Passing + # character(0) excluded nothing, so on the fallback a covariate or + # characteristic named center_size was claimed as a center dummy and + # period_flag as a period dummy. Reading center_size as a center dummy + # adds its coefficient to all_center_lvl_effects and makes that vector one + # longer than center_weights_for_outcome_goal, which shifts every predicted + # outcome and recycles the weights: the #68 defect, one layer down. + named_predictors <- gsub("`", "", c( + "(Intercept)", intervention_components, additional_covariates, + center_characteristics + )) fixed_effect_coefs <- function(term) { - fixed_effect_coef_names(term, coef_mapping, names(all_coefs), character(0)) + fixed_effect_coef_names( + term, coef_mapping, names(all_coefs), named_predictors + ) } if (include_center_effects) { diff --git a/R/validate_inputs.R b/R/validate_inputs.R index b3268ba..80337be 100644 --- a/R/validate_inputs.R +++ b/R/validate_inputs.R @@ -348,12 +348,43 @@ validate_inputs <- function( # check if values of center_weights_for_outcome_goal sum up to 1 if (!is.null(center_weights_for_outcome_goal)) { - if (abs(sum(center_weights_for_outcome_goal) - 1) >= 0.001) { + weights_sum <- sum(center_weights_for_outcome_goal) + if (abs(weights_sum - 1) >= 0.001) { stop(paste( "values in center_weights_for_outcome_goal must", "sum up to 1." )) } + # The tolerance above says the input was MEANT to be a set of weights. It + # does not make it one: the weights multiply the per-center outcomes and are + # then summed, so a set summing to 0.999 scales every reported outcome by + # 0.999 and one summing to 1.001 by 1.001. On the documented example weights + # that is a bias of 6.7e-4 on the reported outcome, in the direction of the + # error and on every outcome the run reports, including the goal comparison + # the recommendation is chosen against. + # + # Renormalising here rather than tightening the tolerance: the tolerance is + # documented and callers rely on it, so tightening it would turn accepted + # input into a hard error. Renormalising is EXACTLY a no-op for a compliant + # caller -- weights already summing to 1 divide by 1 and come back bit for + # bit unchanged, since x / 1 is exact in floating point -- and removes the + # bias for everyone else. + # + # Done once, here, where the weights are validated, so nothing downstream + # needs to know: this is the only place lago_optimization() obtains them, + # whether from the caller, from the center sample sizes, or from a single + # named center, and the value returned from here is what every optimizer, + # the shrinking method and the confidence set are all given. + # + # Silently, deliberately. The correction is at most 0.1% of a weight and is + # what the caller already asked for by passing something the tolerance + # accepts as a set of weights. Warning on it would fire on rounded input + # that the documentation invites -- c(0.333, 0.333, 0.334) is the natural + # way to write three equal weights -- so it would be noise on correct + # usage while the actual behaviour is now right. A caller who did not intend + # weights is still refused, by the check above. + center_weights_for_outcome_goal <- + center_weights_for_outcome_goal / weights_sum } # check if include_time_effects is logical diff --git a/man/select_restart_within_bounds.Rd b/man/select_restart_within_bounds.Rd new file mode 100644 index 0000000..db59d05 --- /dev/null +++ b/man/select_restart_within_bounds.Rd @@ -0,0 +1,67 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/get_recommended_interventions.R +\name{select_restart_within_bounds} +\alias{select_restart_within_bounds} +\title{select_restart_within_bounds} +\usage{ +select_restart_within_bounds( + restart_points, + costs, + lower_bounds, + upper_bounds, + cost_fun +) +} +\arguments{ +\item{restart_points}{A numeric matrix with one column per restart and one +row per intervention component, holding the point each restart converged to.} + +\item{costs}{A numeric vector, one entry per column of restart_points, with +the cost that restart converged to. NA marks a restart whose optimization +failed. Not every entry may be NA: the caller refuses that case with its own +message before calling this.} + +\item{lower_bounds}{A numeric vector. The lower bounds of the intervention +components.} + +\item{upper_bounds}{A numeric vector. The upper bounds of the intervention +components.} + +\item{cost_fun}{A function of one numeric vector returning the total cost of +that intervention.} +} +\value{ +A list with: +- int_components: the chosen intervention, projected onto the bounds. +- rec_int_cost: cost_fun() evaluated at int_components. +} +\description{ +Internal function that picks the recommended intervention out +of the points a multi-start numerical optimizer converged to, and returns it +together with the cost of the point being recommended. +} +\details{ +The filter, the fallback, the choice, the projection and the cost +recomputation are one decision and are deliberately kept together: each +step's correctness depends on the one before it. + +costs holds the cost each restart converged to, so the best restart is the +cheapest one. Every restart satisfies the outcome constraint, which +NlcOptim::solnl() enforces through confun, so cost is the only thing to +choose on among the ones that are actually implementable. + +solnl() treats the box as a soft constraint and will step a little outside it +to buy a lower objective, so the cheapest restart is systematically the one +furthest outside the bounds: selecting on cost alone selects for the +violation. Restarts that left the box are therefore dropped before the +comparison. If every restart left it they are all kept, so a recommendation +is still returned, and the survivor is brought back onto the box. + +The chosen restart can still sit a solver tolerance outside the box, when +every restart did. A recommendation has to be implementable, so it is +projected onto the bounds and its cost recomputed at the value actually being +recommended rather than reported from the point the solver stopped at. The +projection can only raise the cost of a monotone cost function, so reporting +the solver's cost understates what the recommendation costs. +} +\keyword{internal} diff --git a/tests/testthat/test-minimize-and-bounds.R b/tests/testthat/test-minimize-and-bounds.R index a766152..398469e 100644 --- a/tests/testthat/test-minimize-and-bounds.R +++ b/tests/testthat/test-minimize-and-bounds.R @@ -90,23 +90,32 @@ bb_config <- function(outcome_goal, # no center_characteristics argument for that reason and so cannot use # bb_config(). # -# HOW FRAGILE THIS IS. Reaching the fallback requires EVERY restart to stop -# outside the box, and most of them miss it by one to six units in the last -# place, around 1e-16 on a bound of 1. A solver that converges more tightly, -# which a new NlcOptim could reasonably ship, would land some of those restarts -# exactly on the bound instead. The fallback would then stop firing and the two -# tests below would keep passing while no longer covering the projection or the -# cost recomputation at all. They would not fail, they would quietly stop -# testing anything, and nothing here detects that. +# WHAT THIS FIXTURE DOES AND DOES NOT COVER. Reaching the fallback requires +# EVERY restart to stop outside the box, and most of them miss it by one to six +# units in the last place, around 1e-16 on a bound of 1. A solver that converges +# more tightly, which a new NlcOptim could reasonably ship, lands some of those +# restarts exactly on the bound instead: with solnl() at tolX = 1e-8, +# tolFun = tolCon = 1e-9, two restarts come back to a violation of exactly 0, +# the fallback stops firing, and the projection is unreachable from here again. +# It would not fail, it would quietly stop exercising that branch, and nothing +# in this file detects it. Nor can it be detected from the outside: the +# projection puts a component exactly on a bound, and so does a solver that +# converged there, so the returned value cannot distinguish the two. # -# It cannot be detected from the outside either: the projection puts a component -# exactly on a bound, and so does a solver that converged there, so the returned -# value cannot distinguish the two. Covering these two lines durably means -# lifting the projection out of the function nested inside -# optimize_cost_nlcoptim() to somewhere a test can call it directly, with no -# solver in between. That is a source change and is deliberately not made here. -# Until then, treat this fixture as covering those lines TODAY and not as a -# guarantee that it still will after a solver upgrade. +# That is why the projection and the cost recomputation no longer DEPEND on this +# fixture. They live in select_restart_within_bounds(), and +# test-outcome-internals.R calls it directly with a hand-built restart matrix, +# so the out-of-box fallback, the projection and the recomputation are each +# pinned with no solver, model or data in between. Verified: deleting the +# projection, or reporting the solver's cost instead of the cost at the +# projected point, fails those unit tests even with the tighter tolerances above +# applied, i.e. even when this fixture has stopped reaching the fallback at all. +# +# What this fixture still adds, and why it stays: it is the only case in the +# suite where the INTEGRATION path reaches that branch, so it is what shows the +# restart selection is wired into optimize_cost_nlcoptim() correctly and that a +# real solnl() run can leave every restart outside the box. Treat it as covering +# the integration, not as the guarantee for those two lines. bb_three_component_config <- function() { list( data = BB_data, diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index 4647cb7..9bcac65 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -26,6 +26,14 @@ # own: the search is anchored, and it skips the names the assembly can # supply for itself. The mapping is available for every model the suite # fits, so the fallback never runs and neither guard is tested end to end. +# - select_restart_within_bounds() projects the chosen restart onto the +# intervention bounds and recomputes its cost there. Both only do anything +# when EVERY restart left the box, which needs solnl() to stop a tolerance +# outside every one of its bounds on every restart. The three-component +# integration fixture in test-minimize-and-bounds.R does reach it, but only +# because most of its restarts miss the box by around 1e-16, so a solver +# that converged more tightly would disarm it silently. Called directly with +# a hand-built restart matrix there is no solver to depend on. # # getFromNamespace() is how test-helpers.R reaches the other internals, and is # used here for the same reason. @@ -248,3 +256,634 @@ test_that("the fixed-effect fallback is anchored and skips the named predictors" 0 ) }) + + +test_that("both callers pass the names they have already claimed", { + # The guard above is a PARAMETER, so it has two halves: the helper filtering + # against the list, and the caller building a list to filter against. The + # helper half is covered above. This is the caller half, which + # rec_int_processor() used to defeat by passing character(0): a list that + # excludes nothing makes the helper's exclusion a no-op, so the two are not + # independent guards and covering one does not cover the other. + # + # The exclusion list each caller owes is the names it looks up on its own + # account: the intercept, the intervention components, the additional + # covariates and the center characteristics. Built here the way both callers + # build it, over the coefficient names glm() actually produces for a fit with + # real center and period dummies alongside a center_size covariate and a + # period_flag covariate. + fixed_effect_coef_names <- getFromNamespace( + "fixed_effect_coef_names", "LAGO" + ) + term_coef_names <- getFromNamespace("term_coef_names", "LAGO") + + bb <- as.data.frame(BB_data) + bb$center <- factor(rep_len(paste0("c", 1:3), nrow(bb))) + bb$period <- factor(rep_len(1:3, nrow(bb))) + bb$center_size <- bb$staff_nurse + bb$period_flag <- bb$distance_10 + model <- suppressWarnings(glm( + pp3_oxytocin_mother ~ center + period + coaching_updt + + launch_duration + center_size + period_flag, + data = bb, family = binomial() + )) + model_coef_names <- names(coef(model)) + # the fixture has to contain the confusable names, or there is nothing for an + # empty exclusion list to wrongly claim + expect_true(all(c("centerc2", "centerc3", "period2", "period3", + "center_size", "period_flag") %in% model_coef_names)) + + # the list, built exactly as both call sites build it + named_predictors <- gsub("`", "", c( + "(Intercept)", c("coaching_updt", "launch_duration"), + c("center_size", "period_flag"), NULL + )) + + # with the list the callers owe, only the real dummies are claimed + expect_equal( + fixed_effect_coef_names( + "center", NULL, model_coef_names, named_predictors + ), + c("centerc2", "centerc3") + ) + expect_equal( + fixed_effect_coef_names( + "period", NULL, model_coef_names, named_predictors + ), + c("period2", "period3") + ) + + # and EMPTIED, which is what rec_int_processor() used to pass, the covariates + # are claimed as dummies. This is the assertion that fails if either call site + # regresses to character(0). + expect_equal( + fixed_effect_coef_names( + "center", NULL, model_coef_names, character(0) + ), + c("centerc2", "centerc3", "center_size") + ) + expect_equal( + fixed_effect_coef_names( + "period", NULL, model_coef_names, character(0) + ), + c("period2", "period3", "period_flag") + ) + # the two disagree, so the argument is load-bearing rather than incidental + expect_false(identical( + fixed_effect_coef_names( + "center", NULL, model_coef_names, named_predictors + ), + fixed_effect_coef_names("center", NULL, model_coef_names, character(0)) + )) + + # WHY it matters, spelled out in the units rec_int_processor() works in. + # all_center_lvl_effects is the intercept followed by one entry per center + # dummy, and it is averaged against center_weights_for_outcome_goal, which has + # one entry per center. An extra claimed coefficient makes it one longer than + # the weights, so the weights recycle and every predicted outcome shifts. + all_coefs <- coef(model) + effects_of <- function(exclusions) { + dummies <- fixed_effect_coef_names( + "center", NULL, model_coef_names, exclusions + ) + intercept <- all_coefs[["(Intercept)"]] + c(intercept, all_coefs[dummies] + intercept) + } + n_centers <- length(levels(bb$center)) + expect_length(effects_of(named_predictors), n_centers) + expect_length(effects_of(character(0)), n_centers + 1) + + # a center characteristic is excluded on the same footing as a covariate, so + # the list covers both ways center_size can enter the model + as_characteristic <- gsub("`", "", c( + "(Intercept)", c("coaching_updt", "launch_duration"), NULL, "center_size" + )) + expect_equal( + fixed_effect_coef_names( + "center", NULL, model_coef_names, as_characteristic + ), + c("centerc2", "centerc3") + ) + + # For the record on reachability, which is why this is a + # documentation-of-intent guard rather than a live defect: the fallback runs + # only when term_coef_names() returns NULL, and it does not for this model or + # for any model outcome_model_fitting() builds. The mapping resolves the + # dummies exactly and the fallback is never consulted. + mapping <- term_coef_names(model) + expect_false(is.null(mapping)) + expect_equal(mapping$center, c("centerc2", "centerc3")) + expect_equal( + fixed_effect_coef_names( + "center", mapping, model_coef_names, character(0) + ), + c("centerc2", "centerc3") + ) +}) + + +test_that("rec_int_processor() itself excludes the covariates it names", { + # The test above builds the exclusion list the way the callers build it, which + # pins what the list must BE but cannot observe a caller that stops passing + # one: it never runs the caller. This does, through rec_int_processor(), on a + # model whose term mapping has been removed so the fallback is the code path + # actually taken. + # + # The pulesa data has 16 clinics, so a correct run resolves 15 center dummies + # and all_center_lvl_effects has 16 entries, matching the 16 weights. With the + # exclusion list emptied the center_size covariate is claimed as a 16th dummy, + # all_center_lvl_effects becomes 17 long, the weights recycle against it and + # the estimated outcome changes. That difference is what this asserts. + rec_int_processor <- getFromNamespace("rec_int_processor", "LAGO") + term_coef_names <- getFromNamespace("term_coef_names", "LAGO") + + pulesa <- as.data.frame(main_pulesa_data) + pulesa$center <- pulesa$Clinic + # a covariate whose own name begins with "center", i.e. the #68 shape + pulesa$center_size <- 5 + 0.01 * seq_len(nrow(pulesa)) + model <- glm( + Proportions ~ center + AccessMedicines + AccessBPMachines + center_size, + data = pulesa, family = gaussian() + ) + n_centers <- length(levels(pulesa$Clinic)) + expect_equal(n_centers, 16) + + # the fallback is only reached with no mapping, and the mapping is available + # for this model as fitted, so it is removed on purpose here. Emptying + # term.labels is what makes term_coef_names() return NULL. + no_mapping <- model + attr(no_mapping$terms, "term.labels") <- character(0) + expect_false(is.null(term_coef_names(model))) + expect_true(is.null(term_coef_names(no_mapping))) + + run <- function(fitted) { + suppressWarnings(suppressMessages(rec_int_processor( + data = pulesa, + model = fitted, + center_characteristics = NULL, + additional_covariates = "center_size", + include_center_effects = TRUE, + include_time_effects = FALSE, + include_interaction_terms = FALSE, + main_components = NULL, + intervention_components = c("AccessMedicines", "AccessBPMachines"), + optimization_method = "grid_search", + optimization_grid_search_step_size = c(5, 0.5), + link = "identity", + center_weights_for_outcome_goal = rep(1 / n_centers, n_centers), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + intervention_lower_bounds = c(0, 0), + intervention_upper_bounds = c(10, 1), + outcome_goal = 0.6, + center_characteristics_optimization_values = NULL, + time_effect_optimization_value = NULL, + lower_outcome_goal = FALSE, + prev_recommended_interventions = NULL, + shrinkage_threshold = 0.25, + power_goal = NULL, + power_goal_approach = "unconditional", + num_centers_in_next_stage = NULL, + patients_per_center_in_next_stage = NULL, + outcome_name = "Proportions" + ))) + } + + # taking the fallback must give the SAME answer as resolving through the + # mapping: the fallback is a reconstruction of the mapping's result, so the + # two agreeing is the whole requirement on it. They agree only while the + # exclusion list is passed. Emptied, the fallback claims center_size as a 16th + # dummy and the two diverge. + via_mapping <- run(model) + via_fallback <- run(no_mapping) + expect_identical( + via_fallback$est_outcome_goal, via_mapping$est_outcome_goal + ) + expect_identical(via_fallback$rec_int, via_mapping$rec_int) + expect_identical(via_fallback$rec_int_cost, via_mapping$rec_int_cost) + + # and the value itself, so this is not two wrong numbers agreeing. Verified by + # hand below. + expect_equal(via_mapping$est_outcome_goal, -1.28046742730091, + tolerance = 1e-12 + ) + + # the OTHER call site, get_confidence_set(), on the same model and the same + # fallback. It assembles its prediction columns from the resolved names and + # then checks them against the model's coefficients, so an emptied exclusion + # list there does not report a wrong number: the extra claimed dummy makes the + # block one column too wide and the coefficient check refuses the model. That + # this SUCCEEDS is therefore the assertion, and it is what fails if that call + # site stops passing its list. + cs <- suppressWarnings(suppressMessages(get_confidence_set( + predictors_data = pulesa[, c( + "center", "AccessMedicines", "AccessBPMachines", "center_size" + ), drop = FALSE], + include_center_effects = TRUE, + center_weights_for_outcome_goal = rep(1 / n_centers, n_centers), + additional_covariates = "center_size", + intervention_components = c("AccessMedicines", "AccessBPMachines"), + outcome_data = pulesa$Proportions, + fitted_model = no_mapping, + link = "identity", + outcome_goal = 0.6, + outcome_type = "continuous", + intervention_lower_bounds = c(0, 0), + intervention_upper_bounds = c(10, 1), + confidence_set_grid_step_size = c(5, 0.5), + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + rec_int = c(5, 0.5) + ))) + expect_equal( + as.numeric(cs$rec_int_ci), + c(-2.414, -0.262), + tolerance = 1e-3 + ) +}) + + +# --------------------------------------------------------------------------- +# select_restart_within_bounds(): the numerical optimizer's restart selection +# --------------------------------------------------------------------------- + +# a linear total cost, the shape create_cost_function() builds from +# cost_list_of_vectors = list(c(0, 2), c(0, 5)). Linear and increasing, which is +# what makes projecting a below-bound component RAISE the cost, so reporting the +# solver's cost understates the recommendation rather than merely differing from +# it. +srb_cost <- function(x) sum(c(2, 5) * x) +srb_lower <- c(1, 1) +srb_upper <- c(10, 5) + +# the cost function get_recommended_interventions() actually assembles: one +# closure per component from create_cost_function(), reduced with mapply(). Same +# values as srb_cost() on a two-component intervention, but it ERRORS on a +# zero-length intervention instead of returning 0, which is used below where +# that distinction is the point. +srb_cost_as_assembled <- local({ + create_cost_function <- function(coeffs) { + function(x) { + sum(sapply(seq_along(coeffs), function(i) coeffs[i] * x^(i - 1))) + } + } + cost_functions <- lapply(list(c(0, 2), c(0, 5)), create_cost_function) + function(x) sum(mapply(function(f, x) f(x), cost_functions, x)) +}) + +# restart columns and their costs together, so a test cannot accidentally pair a +# cost with the wrong restart: the cost of every restart IS srb_cost() at it, +# which is what solnl() converges to and reports in cost_results. +srb_restarts <- function(...) { + points <- cbind(...) + list(points = points, costs = apply(points, 2, srb_cost)) +} + + +test_that("out-of-box restarts are dropped and the cheapest is taken", { + # The first three steps, on restarts that need no projection. solnl() steps a + # little outside the box to buy a lower objective, so the cheapest restart is + # systematically the one furthest outside the bounds: choosing on cost alone + # chooses the violation. That is why the in-box filter runs FIRST and the + # comparison is only over the survivors. + select_restart_within_bounds <- getFromNamespace( + "select_restart_within_bounds", "LAGO" + ) + + # ALL in box. The filter keeps everything, so the answer is just the cheapest, + # and the projection has nothing to move: the chosen point comes back + # unchanged, bit for bit. + all_in <- srb_restarts(c(2, 2), c(3, 1), c(5, 4)) + expect_equal(all_in$costs, c(14, 11, 30)) + chosen <- select_restart_within_bounds( + all_in$points, all_in$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) + + # SOME in box, and the cheapest restart overall is one of the ones that left + # it: column 1 costs 6, which is cheaper than every survivor, and it is + # discarded because its first component is 0.5 against a lower bound of 1. + # Without the filter it would win, and its projection to c(1, 1) would then be + # recommended at a cost of 7 rather than the genuine in-box optimum of 11. + some_in <- srb_restarts(c(0.5, 1), c(3, 1), c(4, 2)) + expect_equal(some_in$costs, c(6, 11, 18)) + chosen <- select_restart_within_bounds( + some_in$points, some_in$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) + # and it really did pass over a cheaper number + expect_lt(min(some_in$costs), chosen$rec_int_cost) + + # an UPPER-bound violation is dropped just the same as a lower-bound one, so + # the filter is not one-sided + upper_out <- srb_restarts(c(3, 1), c(10, 5.5), c(10.5, 1)) + chosen <- select_restart_within_bounds( + upper_out$points, upper_out$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + + # a SINGLE restart, in box: the degenerate case of the same three steps, and + # the one that a which.min() over an empty selection would break on + single <- srb_restarts(c(3, 1)) + chosen <- select_restart_within_bounds( + single$points, single$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) +}) + + +test_that("a tie in restart cost resolves to the first such restart", { + # which.min() takes the first minimum, so two restarts of equal cost resolve + # to the earlier column. This is not arbitrary: the restarts are ordered by + # their start point along the box diagonal, so the tie-break is deterministic + # and reproducible rather than dependent on the order NlcOptim happened to + # return. Pinning it is what makes a future change of reduction (e.g. to + # which() plus sample(), or to the LAST minimum) visible. + select_restart_within_bounds <- getFromNamespace( + "select_restart_within_bounds", "LAGO" + ) + + # c(3, 1) and c(1, 1.8) both cost 11, exactly, and both are in the box + tie <- srb_restarts(c(3, 1), c(1, 1.8), c(5, 4)) + expect_identical(tie$costs[[1]], tie$costs[[2]]) + expect_identical(tie$costs, c(11, 11, 30)) + chosen <- select_restart_within_bounds( + tie$points, tie$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) + + # the tie-break happens AFTER the filter, not before it: put an equally cheap + # restart outside the box in front of the winner and the winner is unchanged, + # because the out-of-box one is not in the comparison at all + tie_out <- srb_restarts(c(0.5, 2), c(3, 1), c(5, 4)) + expect_identical(tie_out$costs, c(11, 11, 30)) + chosen <- select_restart_within_bounds( + tie_out$points, tie_out$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) +}) + + +test_that("a restart whose optimization failed is not selected", { + # solnl() is wrapped in tryCatch() and a restart that errors leaves NA in + # cost_results with its column of the restart matrix left at whatever it was + # initialised to, i.e. all zeros. An NA cost must therefore exclude the + # restart outright: which.min() ignores NA, but `!is.na(costs)` is what stops + # the all-zeros column from being read as a legitimate in-box point of unknown + # cost. + select_restart_within_bounds <- getFromNamespace( + "select_restart_within_bounds", "LAGO" + ) + + # column 1 is the zeros a failed restart leaves behind. It is OUT of the box + # (0 < lower bound 1) and its cost is NA. + with_na <- srb_restarts(c(0, 0), c(3, 1), c(5, 4)) + costs <- c(NA_real_, with_na$costs[[2]], with_na$costs[[3]]) + chosen <- select_restart_within_bounds( + with_na$points, costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) + + # the harder version: the NA is on a restart that IS in the box, and would be + # the cheapest if its cost were read as 0. The in-box filter alone does not + # exclude it, so this is the NA guard on its own. + na_in_box <- srb_restarts(c(1, 1), c(3, 1), c(5, 4)) + chosen <- select_restart_within_bounds( + na_in_box$points, c(NA_real_, 11, 30), srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(3, 1)) + expect_identical(chosen$rec_int_cost, 11) + + # and when the ONLY in-box restart is the failed one, the fallback keeps the + # out-of-box ones rather than returning nothing, and the cheapest of THOSE is + # projected. c(0.5, 1) costs 6 unprojected and c(1, 1) costs 7. + only_na_in_box <- srb_restarts(c(2, 2), c(0.5, 1), c(11, 2)) + expect_identical(only_na_in_box$costs, c(14, 6, 32)) + chosen <- select_restart_within_bounds( + only_na_in_box$points, c(NA_real_, 6, 32), + srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(1, 1)) + expect_identical(chosen$rec_int_cost, 7) + + # EVERY restart failing selects nothing at all: both the filter and the + # fallback are `which()` over an all-NA condition, so there is no column to + # index and no intervention to project. get_recommended_interventions() + # refuses that case with its own "Numerical optimization failed to find a + # solution" message BEFORE calling this, which is why there is no stop() here. + # This pins that the refusal upstream is load-bearing: reached anyway, the + # cost function the optimizer assembles errors out, so a deleted upstream + # refusal surfaces as an opaque failure rather than as a recommendation. + all_failed <- srb_restarts(c(2, 2), c(3, 1), c(5, 4)) + expect_error( + select_restart_within_bounds( + all_failed$points, rep(NA_real_, 3), srb_lower, srb_upper, + srb_cost_as_assembled + ), + "invalid 'type'" + ) + # the selection itself is what is empty, independently of the cost function: + # with a cost function that tolerates a zero-length argument, the returned + # intervention has no components at all. + chosen <- select_restart_within_bounds( + all_failed$points, rep(NA_real_, 3), srb_lower, srb_upper, srb_cost + ) + expect_length(chosen$int_components, 0) + # and the two cost functions agree on every non-degenerate case, so that + # substitution is not smuggling in different arithmetic + expect_identical(srb_cost(c(3, 1)), srb_cost_as_assembled(c(3, 1))) + expect_identical(srb_cost(c(4, 1)), srb_cost_as_assembled(c(4, 1))) +}) + + +test_that("with no restart in the box the winner is projected and recosted", { + # THE POINT OF THE EXTRACTION. Both the projection and the cost recomputation + # are only reachable when every restart left the box, which from the outside + # needs solnl() to stop a tolerance outside every bound on every restart. Here + # it is just an argument. + # + # A recommendation the user's own bounds forbid is not a recommendation, so + # the winner is brought back onto the box; and the cost has to be that of the + # point being recommended, not of the point the solver stopped at. For an + # increasing cost function the projection can only RAISE the cost, so + # reporting the solver's number understates what the recommendation costs. + select_restart_within_bounds <- getFromNamespace( + "select_restart_within_bounds", "LAGO" + ) + + # none of the three is in the box: column 1 is below the first lower bound, + # column 2 below the second, column 3 above the first upper bound. + none_in <- srb_restarts(c(0.9, 2), c(4, 0.5), c(11, 2)) + expect_identical(none_in$costs, c(11.8, 10.5, 32)) + # the precondition this test is about, asserted rather than assumed + expect_false(any(apply( + none_in$points, 2, + function(x) all(x >= srb_lower) && all(x <= srb_upper) + ))) + + chosen <- select_restart_within_bounds( + none_in$points, none_in$costs, srb_lower, srb_upper, srb_cost + ) + + # the fallback kept them all, so the cheapest of the three wins: column 2, at + # 10.5. Its second component is 0.5, below its lower bound of 1, and the + # PROJECTION is what puts it on the bound. Exactly on it, not near it. + expect_identical(chosen$int_components, c(4, 1)) + expect_identical(chosen$int_components[[2]], srb_lower[[2]]) + + # the RECOMPUTATION: the reported cost is srb_cost() at the projected point, + # 2*4 + 5*1 = 13, and NOT the 10.5 the solver stopped at. The two differ by + # 2.5, so reporting the solver's cost is a 19% understatement of what the + # recommendation actually costs. + expect_identical(chosen$rec_int_cost, 13) + expect_identical(chosen$rec_int_cost, srb_cost(chosen$int_components)) + expect_false(isTRUE(all.equal(chosen$rec_int_cost, none_in$costs[[2]]))) + expect_gt(chosen$rec_int_cost, none_in$costs[[2]]) + + # the same on a SINGLE restart that violates both bounds in opposite + # directions, so the projection is a pmax on one component and a pmin on the + # other in the same call. c(0.25, 7) projects to c(1, 5) and costs 27, against + # the 35.5 the solver reported. + both_ways <- srb_restarts(c(0.25, 7)) + expect_identical(both_ways$costs, 35.5) + chosen <- select_restart_within_bounds( + both_ways$points, both_ways$costs, srb_lower, srb_upper, srb_cost + ) + expect_identical(chosen$int_components, c(1, 5)) + expect_identical(chosen$rec_int_cost, 27) + expect_identical(chosen$rec_int_cost, srb_cost(chosen$int_components)) + # here the projection LOWERS the cost, since the binding violation is on the + # upper bound. The requirement is not a direction, it is that the cost belongs + # to the returned point. + expect_lt(chosen$rec_int_cost, both_ways$costs) + + # the returned intervention is inside the box in every one of these, which is + # the invariant lago_optimization() reports to the user + for (case in list(none_in, both_ways)) { + result <- select_restart_within_bounds( + case$points, case$costs, srb_lower, srb_upper, srb_cost + ) + expect_true(all(result$int_components >= srb_lower)) + expect_true(all(result$int_components <= srb_upper)) + expect_identical(result$rec_int_cost, srb_cost(result$int_components)) + } +}) + + +# --------------------------------------------------------------------------- +# validate_inputs(): the center weights are renormalised, not merely checked +# --------------------------------------------------------------------------- + +test_that("center weights are renormalised exactly when they need it", { + # The tolerance says the input was MEANT to be a set of weights; it does not + # make it one. The weights multiply the per-center outcomes and are summed, so + # a set summing to 1 - d scales EVERY reported outcome by 1 - d: the relative + # bias is exactly sum(w) - 1, on the estimated outcome, on the goal comparison + # the recommendation is chosen against, and on the confidence set. + # + # Renormalising here rather than tightening the tolerance keeps documented + # input accepted. The requirement is two-sided: a no-op for a compliant + # caller, and unbiased for everyone else. + validate_inputs <- getFromNamespace("validate_inputs", "LAGO") + get_outcome <- getFromNamespace("get_outcome", "LAGO") + + pulesa <- as.data.frame(main_pulesa_data) + pulesa$center <- pulesa$Clinic + vi <- function(weights) { + suppressWarnings(suppressMessages(validate_inputs( + data = pulesa, + outcome_name = "Proportions", + outcome_type = "continuous", + intervention_components = c("AccessMedicines", "AccessBPMachines"), + intervention_lower_bounds = c(0, 0), + intervention_upper_bounds = c(10, 1), + outcome_goal = 0.6, + outcome_goal_intention = "maximize", + power_goal = NULL, + power_goal_approach = "unconditional", + cost_list_of_vectors = list(c(0, 1), c(0, 1)), + include_center_effects = TRUE, + center_weights_for_outcome_goal = weights + )))$center_weights_for_outcome_goal + } + + # COMPLIANT: bit-identical, not merely equal. x / 1 is exact in floating + # point, so a caller already summing to 1 is unaffected by construction. + n_centers <- length(levels(pulesa$Clinic)) + compliant <- rep(1 / n_centers, n_centers) + expect_identical(sum(compliant), 1) + expect_identical(vi(compliant), compliant) + + # NON-COMPLIANT but inside the tolerance: 15 weights of 0.0624 and one of + # 0.0632 sum to 0.9992, which the 0.001 check accepts. It comes back + # renormalised, and to the same value as dividing by the sum by hand. + raw <- c(rep(0.0624, n_centers - 1), 0.0632) + expect_equal(sum(raw), 0.9992) + expect_lt(abs(sum(raw) - 1), 0.001) + corrected <- vi(raw) + expect_identical(corrected, raw / sum(raw)) + expect_identical(sum(corrected), 1) + expect_false(identical(corrected, raw)) + + # the bias this removes, on the estimated outcome, computed by hand. The + # identity link makes get_outcome() a plain weighted sum of the per-center + # linear predictors, so scaling the weights scales the result: the relative + # error is exactly sum(w) - 1 and nothing else. + center_effects <- c(0.30, 0.45, 0.55, 0.40, 0.52, 0.35, 0.48, 0.42, + 0.38, 0.50, 0.44, 0.36, 0.46, 0.41, 0.53, 0.39) + beta <- c(0.05, 0.02, 0.10) + int_vector <- c(1, 4, 0.8) + eta <- center_effects + sum(beta * int_vector) - beta[1] + biased <- get_outcome(raw, center_effects, beta, int_vector, 0, 0, "identity") + unbiased <- get_outcome( + corrected, center_effects, beta, int_vector, 0, 0, "identity" + ) + # both hand derivations, independent of the package + expect_equal(biased, sum(raw * eta), tolerance = 1e-14) + expect_equal(unbiased, sum((raw / sum(raw)) * eta), tolerance = 1e-14) + # the relative bias IS sum(w) - 1, i.e. -8e-4 here + expect_equal((biased - unbiased) / unbiased, sum(raw) - 1, tolerance = 1e-9) + expect_equal((biased - unbiased) / unbiased, -8e-4, tolerance = 1e-9) + # and the unbiased value is the one a weighted mean must give: bracketed by + # the values it averages, which the biased one is free to leave + expect_gte(unbiased, min(eta)) + expect_lte(unbiased, max(eta)) + + # the tolerance is still the check that the input was MEANT to be weights, so + # something that is not is refused rather than silently rescaled. 16 weights + # of 0.05 sum to 0.8. + expect_error( + vi(rep(0.05, n_centers)), + "must sum up to 1" + ) + + # the correction is SILENT: no message and no warning. It is at most 0.1% of a + # weight and is what the caller already asked for by passing something the + # tolerance accepts. Warning would fire on rounded input the documentation + # invites, e.g. three weights written as c(0.333, 0.333, 0.334). + messages <- character(0) + warnings <- character(0) + withCallingHandlers( + suppressWarnings(suppressMessages(vi(raw))), + message = function(m) { + messages <<- c(messages, conditionMessage(m)) + invokeRestart("muffleMessage") + }, + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + expect_length(grep("weight", messages, ignore.case = TRUE), 0) + expect_length(grep("weight", warnings, ignore.case = TRUE), 0) + + # the DEFAULT weights, which validate_inputs() derives from the center sample + # sizes rather than taking from the caller, already sum to 1 and so are also + # unaffected. The renormalisation sits after every branch that can produce + # them, so it covers the caller's weights, the sample-size default and the + # single-named-center indicator alike. + expect_identical(sum(vi(NULL)), 1) +}) From 954505575614b3f91cfb00f4a989de09c4129524 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Wed, 5 Aug 2026 18:12:18 +0000 Subject: [PATCH 2/3] Correct the claim about which way the projection moves a cost The new internal's documentation said the projection can only raise the cost of a monotone cost function, so that reporting the solver's cost understates what the recommendation costs. That is wrong in both halves. Projection moves a component up to a lower bound or down to an upper one, so on an increasing cost function it can move the cost either way, and a test in the same file demonstrates it: c(0.25, 7) projects to c(1, 5) and costs 27 against the 35.5 the solver reported. The cost function need not be monotone at all, since create_cost_function() builds whatever polynomial the caller's coefficients describe. The claim appeared in the roxygen, in the generated manual page and in two test comments. All four now say the recomputed cost can be on either side of the solver's, which is the reason it is recomputed rather than adjusted. Two smaller corrections in the same files. The note on the three-component fixture said two restarts return an exactly-zero bound violation under tighter solver tolerances; the count depends on how the objective is built, so it now says several, which is what the argument needs. And an assertion on a recorded outcome value was described as verified by hand below, with no derivation following: it is a fixture recorded from the fixed tree, and says so. --- R/get_recommended_interventions.R | 6 ++++-- man/select_restart_within_bounds.Rd | 6 ++++-- tests/testthat/test-minimize-and-bounds.R | 5 +++-- tests/testthat/test-outcome-internals.R | 22 +++++++++++++--------- 4 files changed, 24 insertions(+), 15 deletions(-) diff --git a/R/get_recommended_interventions.R b/R/get_recommended_interventions.R index 67e9ceb..811ec6c 100644 --- a/R/get_recommended_interventions.R +++ b/R/get_recommended_interventions.R @@ -56,8 +56,10 @@ unachievable_goal_message <- function(lower_outcome_goal) { #' every restart did. A recommendation has to be implementable, so it is #' projected onto the bounds and its cost recomputed at the value actually being #' recommended rather than reported from the point the solver stopped at. The -#' projection can only raise the cost of a monotone cost function, so reporting -#' the solver's cost understates what the recommendation costs. +#' recomputed cost can be either side of the one the solver reported, since +#' projection moves a component up to a lower bound or down to an upper one, so +#' the solver's cost is not a bound on the recommendation's cost in either +#' direction. #' #' @param restart_points A numeric matrix with one column per restart and one #' row per intervention component, holding the point each restart converged to. diff --git a/man/select_restart_within_bounds.Rd b/man/select_restart_within_bounds.Rd index db59d05..8e1751d 100644 --- a/man/select_restart_within_bounds.Rd +++ b/man/select_restart_within_bounds.Rd @@ -61,7 +61,9 @@ The chosen restart can still sit a solver tolerance outside the box, when every restart did. A recommendation has to be implementable, so it is projected onto the bounds and its cost recomputed at the value actually being recommended rather than reported from the point the solver stopped at. The -projection can only raise the cost of a monotone cost function, so reporting -the solver's cost understates what the recommendation costs. +recomputed cost can be either side of the one the solver reported, since +projection moves a component up to a lower bound or down to an upper one, so +the solver's cost is not a bound on the recommendation's cost in either +direction. } \keyword{internal} diff --git a/tests/testthat/test-minimize-and-bounds.R b/tests/testthat/test-minimize-and-bounds.R index 398469e..5989b16 100644 --- a/tests/testthat/test-minimize-and-bounds.R +++ b/tests/testthat/test-minimize-and-bounds.R @@ -95,8 +95,9 @@ bb_config <- function(outcome_goal, # units in the last place, around 1e-16 on a bound of 1. A solver that converges # more tightly, which a new NlcOptim could reasonably ship, lands some of those # restarts exactly on the bound instead: with solnl() at tolX = 1e-8, -# tolFun = tolCon = 1e-9, two restarts come back to a violation of exactly 0, -# the fallback stops firing, and the projection is unreachable from here again. +# tolFun = tolCon = 1e-9, several restarts come back to a violation of exactly +# 0, which is enough for the fallback to stop firing and the projection to be +# unreachable from here again. # It would not fail, it would quietly stop exercising that branch, and nothing # in this file detects it. Nor can it be detected from the outside: the # projection puts a component exactly on a bound, and so does a solver that diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index 9bcac65..d6f16b0 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -461,8 +461,10 @@ test_that("rec_int_processor() itself excludes the covariates it names", { expect_identical(via_fallback$rec_int, via_mapping$rec_int) expect_identical(via_fallback$rec_int_cost, via_mapping$rec_int_cost) - # and the value itself, so this is not two wrong numbers agreeing. Verified by - # hand below. + # and the value itself, so this is not two wrong numbers agreeing. It is a + # fixture, recorded from the fixed tree rather than derived here: the two + # assertions above are what pin the behaviour, and this one only holds them to + # a number that was checked once. expect_equal(via_mapping$est_outcome_goal, -1.28046742730091, tolerance = 1e-12 ) @@ -506,10 +508,11 @@ test_that("rec_int_processor() itself excludes the covariates it names", { # --------------------------------------------------------------------------- # a linear total cost, the shape create_cost_function() builds from -# cost_list_of_vectors = list(c(0, 2), c(0, 5)). Linear and increasing, which is -# what makes projecting a below-bound component RAISE the cost, so reporting the -# solver's cost understates the recommendation rather than merely differing from -# it. +# cost_list_of_vectors = list(c(0, 2), c(0, 5)). Linear and increasing, so +# projecting a component UP to a lower bound raises the cost and projecting one +# DOWN to an upper bound lowers it. The recomputed cost is therefore not on a +# predictable side of the solver's, which is why it is recomputed rather than +# adjusted. srb_cost <- function(x) sum(c(2, 5) * x) srb_lower <- c(1, 1) srb_upper <- c(10, 5) @@ -707,9 +710,10 @@ test_that("with no restart in the box the winner is projected and recosted", { # # A recommendation the user's own bounds forbid is not a recommendation, so # the winner is brought back onto the box; and the cost has to be that of the - # point being recommended, not of the point the solver stopped at. For an - # increasing cost function the projection can only RAISE the cost, so - # reporting the solver's number understates what the recommendation costs. + # point being recommended, not of the point the solver stopped at. Which way + # the cost then moves depends on which bound was crossed, so the solver's + # number is not usable as an estimate of it: the test below projects a + # component down to an upper bound and the cost falls, from 35.5 to 27. select_restart_within_bounds <- getFromNamespace( "select_restart_within_bounds", "LAGO" ) From 51b6e8ff9a80fe7420cf8867c659d17956b79f74 Mon Sep 17 00:00:00 2001 From: Ante Bing Date: Wed, 5 Aug 2026 18:25:57 +0000 Subject: [PATCH 3/3] Test that the caller forwards additional_covariates, and correct two claims Review found that deleting the single line in lago_optimization() that forwards additional_covariates to rec_int_processor() left the whole suite green. The existing test drives rec_int_processor() directly, so it pins what the callee does with the argument and never runs the caller, which is the same structural gap one level up from the one this branch set out to close. The fallback the argument protects needs a model whose term mapping is missing, and lago_optimization() never fits one, so it cannot be reached by running an optimization: the test instead replaces rec_int_processor() for one call and asserts what it received. Deleting the forwarding line now fails it. Two claims in comments were wrong. The note on renormalising center weights said it is exactly a no-op for a compliant caller, which does not hold on the default weights: a vector of center sizes divided by its own total need not sum to exactly 1, and about one such vector in four hundred is a unit in the last place away, so those runs move by that much. It now says so, and says the movement is the correction working rather than an error introduced. The note in rec_int_processor() said the term mapping is always available for models this package fits, making the fallback unreachable. A model fitted with model = FALSE whose data has left scope cannot rebuild its mapping, and the exported get_confidence_set() takes whatever model the caller passes, so the fallback is reachable and the guard on it is not decoration. --- R/rec_int_processor.R | 11 ++++--- R/validate_inputs.R | 12 ++++--- tests/testthat/test-outcome-internals.R | 43 +++++++++++++++++++++++++ 3 files changed, 58 insertions(+), 8 deletions(-) diff --git a/R/rec_int_processor.R b/R/rec_int_processor.R index d8af126..918bba1 100644 --- a/R/rec_int_processor.R +++ b/R/rec_int_processor.R @@ -82,10 +82,13 @@ rec_int_processor <- function( # center-level effects would shift every predicted outcome by its # coefficient, and the first-element lookups downstream would then read a # coefficient that is not a center effect at all. - # The fallback for a model whose term mapping could not be rebuilt, which - # outcome_model_fitting() does not produce, is the name search this used to - # do, anchored so that only a name beginning with the term is matched, and - # restricted to the coefficients no other block below claims for itself. + # The fallback for a model whose term mapping could not be rebuilt is the name + # search this used to do, anchored so that only a name beginning with the term + # is matched, and restricted to the coefficients no other block below claims + # for itself. It is reachable: a model fitted with model = FALSE whose data + # has since left scope cannot rebuild its mapping, and the exported + # get_confidence_set() accepts any model the caller passes, so the fallback is + # what a caller in that position gets rather than dead code. # # That exclusion list is the point of named_predictors, and it has to be the # names THIS function looks up on its own account, exactly as diff --git a/R/validate_inputs.R b/R/validate_inputs.R index 80337be..f50a307 100644 --- a/R/validate_inputs.R +++ b/R/validate_inputs.R @@ -365,10 +365,14 @@ validate_inputs <- function( # # Renormalising here rather than tightening the tolerance: the tolerance is # documented and callers rely on it, so tightening it would turn accepted - # input into a hard error. Renormalising is EXACTLY a no-op for a compliant - # caller -- weights already summing to 1 divide by 1 and come back bit for - # bit unchanged, since x / 1 is exact in floating point -- and removes the - # bias for everyone else. + # input into a hard error. Renormalising leaves weights that sum to exactly + # 1 bit for bit unchanged, since dividing by 1 is exact, and removes the + # bias for everyone else. It is not quite a no-op on the defaults computed + # below: a vector of center sizes divided by its own total need not sum to + # exactly 1, and about one such vector in four hundred sums to one unit in + # the last place away from it, so those runs move by the same amount. That + # is the correction doing its job on a genuinely non-unit sum, not an error + # introduced here. # # Done once, here, where the weights are validated, so nothing downstream # needs to know: this is the only place lago_optimization() obtains them, diff --git a/tests/testthat/test-outcome-internals.R b/tests/testthat/test-outcome-internals.R index d6f16b0..4f134ef 100644 --- a/tests/testthat/test-outcome-internals.R +++ b/tests/testthat/test-outcome-internals.R @@ -502,6 +502,49 @@ test_that("rec_int_processor() itself excludes the covariates it names", { ) }) +test_that("lago_optimization() passes additional_covariates on to the processor", { + # The test above runs rec_int_processor() directly, so it pins what the callee + # does with the argument but not that its caller supplies it. Deleting the one + # line in lago_optimization() that forwards additional_covariates leaves that + # test green, because nothing else drives the caller. + # + # The fallback the argument protects needs a model whose term mapping is + # missing, which lago_optimization() never fits, so this cannot be reached + # end to end by running an optimization. What can be checked is the forwarding + # itself: rec_int_processor() is replaced for the duration of one call and + # asked what it received. + seen <- new.env(parent = emptyenv()) + real <- getFromNamespace("rec_int_processor", "LAGO") + spy <- function(...) { + args <- list(...) + seen$additional_covariates <- args$additional_covariates + do.call(real, args) + } + + testthat::with_mocked_bindings( + { + suppressWarnings(suppressMessages(lago_optimization( + data = BB_data, + outcome_name = "pp3_oxytocin_mother", + outcome_type = "binary", + glm_family = "binomial", + intervention_components = c("coaching_updt", "launch_duration"), + additional_covariates = "birth_volume_100", + intervention_lower_bounds = c(1, 1), + intervention_upper_bounds = c(40, 5), + cost_list_of_vectors = list(c(0, 1.7), c(0, 8)), + outcome_goal = 0.85, + include_confidence_set = FALSE, + quiet = TRUE + ))) + }, + rec_int_processor = spy, + .package = "LAGO" + ) + + expect_identical(seen$additional_covariates, "birth_volume_100") +}) + # --------------------------------------------------------------------------- # select_restart_within_bounds(): the numerical optimizer's restart selection