Refactor/projection and validation gaps

R/get_recommended_interventions.R

@@ -30,6 +30,86 @@ 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
+#' 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.
+#' @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 +665,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,

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,

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,
@@ -81,11 +82,32 @@ 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.
+  # 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
+  # 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) {

R/validate_inputs.R

@@ -348,12 +348,47 @@ 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 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,
+    # 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

tests/testthat/test-minimize-and-bounds.R

@@ -90,23 +90,33 @@ 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, 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
+# 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,