Nothing
#' @include internal.R
# -------------------------------------------------------------------------
# Utilities
# -------------------------------------------------------------------------
# -------------------------------------------------------------------------
# Promote Problem -> Problem (internal; never exposed to user)
# -------------------------------------------------------------------------
# .pamo_as_mo <- function(x) {
# if (inherits(x, "Problem")) return(x)
#
# if (inherits(x, "Problem")) {
# obj <- pproto(NULL, Problem, base = x)
#
# # fuerza clase S3 para que UseMethod("solve") encuentre solve.Problem
# cls <- class(obj)
# if (is.null(cls)) cls <- character()
# class(obj) <- unique(c("Problem", cls))
#
# return(obj)
# }
#
# stop("Expected a Problem or a Problem.", call. = FALSE)
# }
# -------------------------------------------------------------------------
# Atomic objective registry accessors (stored in multiscape::Problem)
# Registry: x$data$objectives[[alias]] = list(...)
# -------------------------------------------------------------------------
.pamo_get_specs <- function(x) {
stopifnot(inherits(x, "Problem"))
specs <- x$data$objectives %||% list()
if (!is.list(specs)) specs <- list()
specs
}
# internal: get ONE atomic objective spec by alias
.pamo_get_objective_spec <- function(x, alias) {
stopifnot(inherits(x, "Problem"))
alias <- as.character(alias)[1]
if (is.na(alias) || !nzchar(alias)) {
stop("alias must be a non-empty string.", call. = FALSE)
}
specs <- .pamo_get_specs(x)
sp <- specs[[alias]]
if (is.null(sp)) {
stop(
"Objective alias not found: '", alias, "'.\n",
"Tip: register objectives via add_objective_* (alias='...') before calling set_method_*().",
call. = FALSE
)
}
# defensive minimal validation (clearer errors later)
if (is.null(sp$objective_id)) stop("Objective '", alias, "' has no objective_id.", call. = FALSE)
if (is.null(sp$sense)) stop("Objective '", alias, "' has no sense.", call. = FALSE)
sp
}
# internal: get MANY objective specs in order (and validate duplicates)
.pamo_get_objective_specs <- function(x, aliases) {
stopifnot(inherits(x, "Problem"))
aliases <- as.character(aliases)
if (length(aliases) == 0L) stop("aliases must have length > 0.", call. = FALSE)
if (anyNA(aliases) || any(!nzchar(aliases))) {
stop("aliases must be non-empty strings.", call. = FALSE)
}
if (anyDuplicated(aliases) != 0L) {
dups <- unique(aliases[duplicated(aliases)])
stop("Duplicated objective aliases: ", paste(dups, collapse = ", "), call. = FALSE)
}
specs <- lapply(aliases, function(a) .pamo_get_objective_spec(x, a))
names(specs) <- aliases
specs
}
# internal: validate registry sanity
.pamo_validate_objectives <- function(x) {
stopifnot(inherits(x, "Problem"))
specs <- .pamo_get_specs(x)
if (length(specs) == 0) {
stop("No objectives registered. Use add_objective_* (with alias=...) first.", call. = FALSE)
}
al <- names(specs)
if (anyNA(al) || any(!nzchar(al))) stop("Objective registry has invalid aliases.", call. = FALSE)
if (anyDuplicated(al)) stop("Duplicated objective aliases are not allowed.", call. = FALSE)
invisible(TRUE)
}
# -------------------------------------------------------------------------
# Objective IR (intermediate representation)
# Each IR carries:
# - sense: "min"|"max"
# - terms: list of atomic terms (still indivisible objectives)
# - objective_id, objective_args: to be able to re-activate the objective in a single-objective model
# -------------------------------------------------------------------------
.pamo_objvec_action_boundary_cut <- function(base_superset, term) {
stopifnot(inherits(base_superset, "Problem"))
op <- base_superset$data$model_ptr
if (is.null(op)) {
stop("Superset model is not built (missing model_ptr).", call. = FALSE)
}
ml <- base_superset$data$model_list %||% NULL
if (is.null(ml) || !is.list(ml)) {
base_superset <- .pa_refresh_model_snapshot(base_superset)
ml <- base_superset$data$model_list %||% NULL
}
if (is.null(ml) || !is.list(ml)) {
stop("Model snapshot is missing in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
m0 <- .pa_model_from_ptr(
op,
args = base_superset$data$model_args %||% list(),
drop_triplets = TRUE
)
n_col <- as.integer(
ml$ncol %||%
ml$n_col %||%
length(m0$obj)
)
n_x <- as.integer(ml$n_x %||% NA_integer_)
x_offset <- as.integer(ml$x_offset %||% NA_integer_)
y0 <- as.integer(ml$y_action_offset %||% NA_integer_)
n_y <- as.integer(ml$n_y_action %||% NA_integer_)
if (!is.finite(n_col) || is.na(n_col) || n_col <= 0L) {
stop("Invalid model dimension (ncol) in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
if (!is.finite(n_x) || is.na(n_x) || n_x <= 0L) {
stop("Invalid n_x in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
if (!is.finite(x_offset) || is.na(x_offset) || x_offset < 0L) {
stop("Invalid x_offset in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
if (!is.finite(y0) || is.na(y0) || y0 < 0L) {
stop("Invalid y_action_offset in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
if (!is.finite(n_y) || is.na(n_y) || n_y <= 0L) {
stop("Invalid n_y_action in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
v <- numeric(n_col)
# ---- relation
rel_name <- as.character(term$relation_name %||% "boundary")[1]
rel_model <- base_superset$data$spatial_relations_model[[rel_name]] %||%
base_superset$data$spatial_relations[[rel_name]] %||%
NULL
if (is.null(rel_model) || !inherits(rel_model, "data.frame")) {
stop("Missing relation '", rel_name, "' in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
for (nm in c("internal_pu1", "internal_pu2", "weight")) {
if (!nm %in% names(rel_model)) {
stop("Relation '", rel_name, "' must contain column '", nm, "'.", call. = FALSE)
}
}
# ---- dist_actions model-ready
da <- base_superset$data$dist_actions_model %||% NULL
if (is.null(da) || !inherits(da, "data.frame")) {
stop("dist_actions_model is missing in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
for (nm in c("internal_row", "internal_pu", "internal_action")) {
if (!nm %in% names(da)) {
stop("dist_actions_model must contain column '", nm, "'.", call. = FALSE)
}
}
# ---- infer global number of actions
all_actions <- sort(unique(as.integer(da$internal_action)))
all_actions <- all_actions[is.finite(all_actions) & !is.na(all_actions)]
if (length(all_actions) == 0L) {
stop("Could not infer internal actions in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
n_actions_global <- max(all_actions)
# ---- selected actions
act_raw <- term$actions %||% NULL
if (is.null(act_raw) || length(act_raw) == 0L) {
act_int <- all_actions
} else {
numeric_ids <- suppressWarnings(as.integer(act_raw))
act_int <- if (all(!is.na(numeric_ids)) && all(numeric_ids %in% all_actions)) {
numeric_ids
} else {
.pa_resolve_action_subset(base_superset, act_raw)$internal_id
}
act_int <- sort(unique(act_int[is.finite(act_int) & !is.na(act_int)]))
}
if (length(act_int) == 0L) {
return(v)
}
# ---- action weights
aw <- term$action_weights %||% NULL
if (is.null(aw)) {
aw_map <- stats::setNames(rep(1, length(act_int)), act_int)
} else {
aw <- as.numeric(aw)
if (length(aw) == length(all_actions)) {
tmp_map <- stats::setNames(aw, all_actions)
aw_map <- tmp_map[as.character(act_int)]
} else if (length(aw) == length(act_int)) {
aw_map <- stats::setNames(aw, act_int)
} else {
stop(
"action_weights length must match either all actions or the selected action subset.",
call. = FALSE
)
}
}
weight_multiplier <- as.numeric(term$weight_multiplier %||% 1)[1]
if (!is.finite(weight_multiplier)) weight_multiplier <- 1
n_pu <- nrow(base_superset$data$pu)
if (!is.finite(n_pu) || is.na(n_pu) || n_pu <= 0L) {
stop("Invalid n_pu in .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
ip1 <- as.integer(rel_model$internal_pu1)
ip2 <- as.integer(rel_model$internal_pu2)
ww <- as.numeric(rel_model$weight)
directed <- if ("directed" %in% names(rel_model)) unique(as.logical(rel_model$directed)) else FALSE
if (length(directed) != 1L || is.na(directed)) stop("Invalid directed metadata in spatial relation.", call. = FALSE)
support_a <- pmin(ip1[ip1 != ip2], ip2[ip1 != ip2])
support_b <- pmax(ip1[ip1 != ip2], ip2[ip1 != ip2])
support <- unique(data.frame(a = support_a, b = support_b))
support <- support[order(support$a, support$b), , drop = FALSE]
k_edges <- nrow(support)
if (k_edges * n_actions_global != n_y) stop("Mismatch in y_action layout.", call. = FALSE)
support_key <- paste(support$a, support$b, sep = "::")
edge_index <- stats::setNames(seq_len(k_edges), support_key)
key2 <- function(i, a) paste0(i, "::", a)
x_map <- stats::setNames(x_offset + as.integer(da$internal_row) - 1L,
key2(da$internal_pu, da$internal_action))
for (act in act_int) {
awi <- as.numeric(aw_map[[as.character(act)]])
if (!is.finite(awi) || awi == 0) next
for (r in seq_along(ip1)) {
i <- ip1[r]; j <- ip2[r]; wij <- ww[r]
xi0 <- unname(x_map[key2(i, act)])[1]
if (i == j) {
if (directed) stop("Directed spatial relations cannot contain self-arcs.", call. = FALSE)
if (is.finite(xi0)) v[xi0 + 1L] <- v[xi0 + 1L] + weight_multiplier * awi * wij
next
}
if (is.finite(xi0)) v[xi0 + 1L] <- v[xi0 + 1L] + weight_multiplier * awi * wij
if (!directed) {
xj0 <- unname(x_map[key2(j, act)])[1]
if (is.finite(xj0)) v[xj0 + 1L] <- v[xj0 + 1L] + weight_multiplier * awi * wij
}
e <- unname(edge_index[paste(min(i, j), max(i, j), sep = "::")])
bcol0 <- y0 + (e - 1L) * n_actions_global + (act - 1L)
v[bcol0 + 1L] <- v[bcol0 + 1L] - weight_multiplier * awi * wij * if (directed) 1 else 2
}
}
v}
.pamo_objective_to_ir <- function(x, spec) {
stopifnot(inherits(x, "Problem"))
stopifnot(is.list(spec), !is.null(spec$objective_id))
id <- as.character(spec$objective_id)[1]
a <- spec$objective_args %||% list()
sense <- as.character(spec$sense %||% NA_character_)[1]
if (is.na(sense) || !nzchar(sense)) sense <- NA_character_
.c1 <- function(z, default = NULL) {
if (is.null(z)) return(default)
as.character(z)[1]
}
.n1 <- function(z, default = NULL) {
if (is.null(z)) return(default)
as.numeric(z)[1]
}
.l1 <- function(z, default = FALSE) {
if (is.null(z)) return(default)
isTRUE(z)
}
.chr <- function(z) {
if (is.null(z)) return(NULL)
z <- as.character(z)
z <- unique(z[!is.na(z) & nzchar(z)])
if (length(z) == 0) return(NULL)
z
}
# ------------------------------------------------------------------
# min_cost
# ------------------------------------------------------------------
if (identical(id, "min_cost")) {
inc_pu <- .l1(a$include_pu_cost, TRUE)
inc_act <- .l1(a$include_action_cost, TRUE)
actions <- .chr(a$actions)
feats <- .chr(a$features)
terms <- list()
if (inc_pu) {
terms <- c(terms, list(list(
type = "pu_cost",
features = feats
)))
}
if (inc_act) {
# Cost setters store external ids. Resolve them before the shared IR
# evaluator, where numeric selectors denote internal column ids.
action_internal <- if (is.null(actions)) NULL else
as.character(.pa_resolve_action_subset(x, actions)$internal_id)
terms <- c(terms, list(list(
type = "action_cost",
actions = action_internal
)))
}
a$include_pu_cost <- inc_pu
a$include_action_cost <- inc_act
a$actions <- actions
a$features <- feats
return(list(
sense = "min",
terms = terms,
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# max_benefit
# ------------------------------------------------------------------
if (id %in% c("max_benefit", "max_effect", "min_effect")) {
bcol <- .c1(a$benefit_col, "effect")
actions <- .chr(a$actions)
feats <- .chr(a$features)
a$effect_sense <- if (id == "min_effect") "min" else "max"
a$benefit_col <- bcol
a$actions <- actions
a$features <- feats
return(list(
sense = if (id == "min_effect") "min" else "max",
terms = list(list(
type = "benefit",
benefit_col = bcol,
actions = actions,
features = feats
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# max_profit
# ------------------------------------------------------------------
if (identical(id, "max_profit")) {
pcol <- .c1(a$profit_col, "profit")
actions <- .chr(a$actions)
a$profit_col <- pcol
a$actions <- actions
return(list(
sense = "max",
terms = list(list(
type = "profit",
profit_col = pcol,
actions = actions
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# max_net_profit
# ------------------------------------------------------------------
if (identical(id, "max_net_profit")) {
pcol <- .c1(a$profit_col, "profit")
inc_pu <- .l1(a$include_pu_cost, TRUE)
inc_act <- .l1(a$include_action_cost, TRUE)
actions <- .chr(a$actions)
feats <- .chr(a$features)
a$profit_col <- pcol
a$include_pu_cost <- inc_pu
a$include_action_cost <- inc_act
a$actions <- actions
a$features <- feats
return(list(
sense = "max",
terms = list(list(
type = "net_profit",
profit_col = pcol,
include_pu_cost = inc_pu,
include_action_cost = inc_act,
actions = actions,
features = feats
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# min_fragmentation
# ------------------------------------------------------------------
if (identical(id, "min_fragmentation")) {
rel <- .c1(a$relation_name, "boundary")
mul <- .n1(a$weight_multiplier, 1)
feats <- .chr(a$features)
a$relation_name <- rel
a$weight_multiplier <- mul
a$features <- feats
return(list(
sense = "min",
terms = list(list(
type = "boundary_cut",
relation_name = rel,
weight_multiplier = mul,
features = feats
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# min_action_fragmentation
# ------------------------------------------------------------------
if (identical(id, "min_action_fragmentation")) {
rel <- .c1(a$relation_name, "boundary")
mul <- .n1(a$weight_multiplier, 1)
actions <- .chr(a$actions)
aw <- a$action_weights %||% NULL
if (!is.null(aw)) aw <- as.numeric(aw)
a$relation_name <- rel
a$weight_multiplier <- mul
a$actions <- actions
a$action_weights <- aw
return(list(
sense = "min",
terms = list(list(
type = "action_boundary_cut",
relation_name = rel,
weight_multiplier = mul,
actions = actions,
action_weights = aw
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# min_loss
# ------------------------------------------------------------------
if (identical(id, "min_loss")) {
lcol <- .c1(a$loss_col, "loss")
actions <- .chr(a$actions)
feats <- .chr(a$features)
a$loss_col <- lcol
a$actions <- actions
a$features <- feats
return(list(
sense = "min",
terms = list(list(
type = "loss",
actions = actions,
features = feats
)),
objective_id = id,
objective_args = a
))
}
# ------------------------------------------------------------------
# min_intervention_impact
# ------------------------------------------------------------------
if (identical(id, "min_intervention_impact")) {
icol <- .c1(a$impact_col, "amount")
feats <- .chr(a$features)
actions <- .chr(a$actions)
a$impact_col <- icol
a$features <- feats
a$actions <- actions
return(list(
sense = "min",
terms = list(list(
type = "intervention_impact",
impact_col = icol,
features = feats,
actions = actions
)),
objective_id = id,
objective_args = a
))
}
stop("Unknown objective_id in .pamo_objective_to_ir(): ", id, call. = FALSE)
}
# -------------------------------------------------------------------------
# Cloning base Problem safely for MO runs
# (Avoid copying externalptr / built model pointer)
# -------------------------------------------------------------------------
.pamo_clone_base <- function(base) {
stopifnot(inherits(base, "Problem"))
b <- pproto(NULL, base)
# copia superficial del contenedor
b$data <- base$data
# independizar solo lo mutable
b$data$model_args <- unserialize(serialize(base$data$model_args %||% list(), NULL))
b$data$meta <- unserialize(serialize(base$data$meta %||% list(), NULL))
# limpiar estado de modelo para este run
b$data$model_ptr <- NULL
b$data$has_model <- FALSE
b$data$model_list <- NULL
b$data$model_index <- NULL
b$data$runtime_updates <- NULL
b$data$mo_cache <- NULL
b$data$meta$model_dirty <- TRUE
b
}
# -------------------------------------------------------------------------
# Activate an IR as a single-objective config in multiscape::Problem
# (Used by the "rebuild + pad" objective vector strategy)
# -------------------------------------------------------------------------
.pamo_activate_ir_as_single_objective <- function(x, ir) {
stopifnot(inherits(x, "Problem"))
stopifnot(is.list(ir), !is.null(ir$objective_id))
id <- as.character(ir$objective_id)[1]
a <- ir$objective_args %||% list()
map <- list(
min_cost = "minimizeCosts",
max_benefit = "maximizeBenefits",
max_effect = "maximizeBenefits",
min_effect = "maximizeBenefits",
max_profit = "maximizeProfit",
min_loss = "minimizeLosses",
max_net_profit = "maximizeNetProfit",
min_fragmentation = "minimizeFragmentation",
min_action_fragmentation = "minimizeActionFragmentation",
min_intervention_impact = "minimizeInterventionImpact",
custom = "custom"
)
mt <- map[[id]]
if (is.null(mt)) stop("No model_type mapping for objective_id: ", id, call. = FALSE)
if (is.null(x$data$model_args) || !is.list(x$data$model_args)) x$data$model_args <- list()
x$data$model_args$model_type <- mt
x$data$model_args$objective_id <- id
x$data$model_args$objective_args <- a
if (is.null(x$data$meta) || !is.list(x$data$meta)) x$data$meta <- list()
x$data$meta$model_dirty <- TRUE
x
}
# -------------------------------------------------------------------------
# Superset model selection + objective-vector extraction
# Strategy:
# - Build ONE "superset" model using the objective that likely introduces most aux vars
# - For each objective, build its own model and extract model$obj
# - Pad with zeros to match superset length
# This avoids dealing with offsets/variable mapping for now.
# -------------------------------------------------------------------------
.pamo_prepare_relation_model <- function(rel) {
stopifnot(is.data.frame(rel), nrow(rel) > 0)
need <- c("internal_pu1", "internal_pu2", "weight")
if (!all(need %in% names(rel))) {
stop(
"Spatial relation must contain columns: ",
paste(need, collapse = ", "),
call. = FALSE
)
}
rel <- rel[, c("internal_pu1","internal_pu2","weight",
intersect(names(rel), c("directed","distance","source","relation_name"))),
drop = FALSE]
rel$internal_pu1 <- as.integer(rel$internal_pu1)
rel$internal_pu2 <- as.integer(rel$internal_pu2)
rel$weight <- as.numeric(rel$weight)
if (anyNA(rel$internal_pu1) || anyNA(rel$internal_pu2)) {
stop("Relation has NA in internal_pu1/internal_pu2.", call. = FALSE)
}
if (any(!is.finite(rel$weight)) || any(rel$weight < 0)) {
stop("Relation has non-finite or negative weights.", call. = FALSE)
}
rel <- rel[order(rel$internal_pu1, rel$internal_pu2), , drop = FALSE]
rel$internal_edge <- seq_len(nrow(rel))
rel
}
.pamo_prepare_superset_model <- function(base, ir_list) {
stopifnot(inherits(base, "Problem"))
stopifnot(is.list(ir_list), length(ir_list) > 0)
spec <- .pamo_compile_superset_spec(base, ir_list)
# elegir un objetivo "semilla" válido para que el core del modelo se materialice.
# preferimos min_cost si existe; si no, el primero.
pick_seed_ir <- function(ir_list) {
idx_cost <- which(vapply(
ir_list,
function(ir) identical(ir$objective_id %||% "", "min_cost"),
logical(1)
))
if (length(idx_cost) > 0L) return(ir_list[[idx_cost[1]]])
ir_list[[1]]
}
ir_seed <- pick_seed_ir(ir_list)
b <- .pamo_clone_base(base)
b <- .pamo_activate_ir_as_single_objective(b, ir_seed)
b$data$model_args <- b$data$model_args %||% list()
b$data$model_args$mo_mode <- TRUE
# la clave del superset: needs estructurales, independientes del método
b$data$model_args$needs <- utils::modifyList(
b$data$model_args$needs %||% list(),
spec$needs
)
if (!is.null(spec$relation_name)) {
b$data$model_args$needs$relation_name <- spec$relation_name
}
b <- .pa_build_model(b)
op <- b$data$model_ptr
if (is.null(op)) {
stop("Superset build failed: model_ptr is NULL.", call. = FALSE)
}
# sanity check mÃnimo: si hay objetivos de representación, el modelo debe tener z
if (isTRUE(spec$needs$z)) {
op_list <- .pa_model_from_ptr(
op,
args = b$data$model_args %||% list(),
drop_triplets = TRUE
)
n_z <- as.integer(op_list$n_z %||% 0L)
if (n_z <= 0L) {
stop(
"MO superset requires z variables, but the built model has n_z = 0.",
call. = FALSE
)
}
}
# sanity check PU fragmentation
if (isTRUE(spec$needs$y_pu)) {
op_list <- .pa_model_from_ptr(
op,
args = b$data$model_args %||% list(),
drop_triplets = TRUE
)
n_y_pu <- as.integer(op_list$n_y_pu %||% 0L)
if (n_y_pu <= 0L) {
stop(
"MO superset requires PU fragmentation auxiliaries, but the built model has n_y_pu = 0.",
call. = FALSE
)
}
}
b
}
.pamo_model_obj_length <- function(x) {
stopifnot(inherits(x, "Problem"))
if (is.null(x$data$model_ptr)) stop("Problem has no model_ptr; build model first.", call. = FALSE)
m <- .pa_model_from_ptr(
x$data$model_ptr,
args = x$data$model_args %||% list(),
drop_triplets = TRUE
)
length(m$obj)
}
# -------------------------------------------------------------------------
# Optional normalization helpers (for obj vectors)
# NOTE: this is normalization of objective *vectors* (not weights).
# If you don't want it now, keep but don't use.
# -------------------------------------------------------------------------
.pamo_normalize_vec <- function(v, method = c("max", "l1", "l2")) {
method <- match.arg(method)
v <- as.numeric(v)
if (!any(is.finite(v))) return(v)
if (method == "max") {
s <- max(abs(v), na.rm = TRUE)
if (is.finite(s) && s > 0) return(v / s)
return(v)
}
if (method == "l1") {
s <- sum(abs(v), na.rm = TRUE)
if (is.finite(s) && s > 0) return(v / s)
return(v)
}
if (method == "l2") {
s <- sqrt(sum(v^2, na.rm = TRUE))
if (is.finite(s) && s > 0) return(v / s)
return(v)
}
v
}
# ---------------------------------------------------------
# Internal: solve weighted
# ---------------------------------------------------------
.pamo_solve_weighted <- function(x, ...) {
stopifnot(inherits(x, "Problem"))
data <- x$data %||% list()
method <- x$data$method %||% list()
# method <- .pamo_normalize_legacy_method_runs(method)
aliases <- as.character(method$aliases %||% character(0))
normalize_weights <- isTRUE(method$normalize_weights)
objective_scaling <- isTRUE(method$objective_scaling)
dots <- list(...)
verbose <- .pamo_cli_enabled(x, dots)
gap_limit <- dots$gap_limit %||% NULL
time_limit <- dots$time_limit %||% NULL
.pamo_cli_header("Weighted method", verbose = verbose)
.pamo_cli_step("Preparing weighted design.", verbose = verbose)
if (length(aliases) == 0L) {
stop("Weighted method: missing aliases.", call. = FALSE)
}
design_df <- .pamo_compile_weighted_runs(
x = x,
method = method
)
n_runs <- nrow(design_df)
.pamo_cli_step(
"Design contains {.val {n_runs}} run{?s}.",
verbose = verbose
)
scales <- NULL
if (objective_scaling) {
.pamo_cli_step(
"Computing payoff ranges for {.field objective_scaling = TRUE}.",
verbose = verbose
)
scales <- .pamo_compute_weighted_ranges(
x = x,
aliases = aliases,
verbose = verbose,
gap_limit = gap_limit,
time_limit = time_limit
)
}
solutions <- vector("list", n_runs)
value_mat <- matrix(NA_real_, n_runs, length(aliases))
colnames(value_mat) <- paste0("value_", aliases)
status <- character(n_runs)
runtime <- numeric(n_runs)
gap <- numeric(n_runs)
for (r in seq_len(n_runs)) {
w_r <- as.numeric(design_df[r, paste0("weight_", aliases), drop = TRUE])
one <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = aliases,
weights = w_r,
normalize_weights = normalize_weights,
objective_scaling = objective_scaling,
scales = scales,
gap_limit = gap_limit,
time_limit = time_limit
)
)
solutions[r] <- list(one$solution)
status[r] <- as.character(one$status %||% NA_character_)
runtime[r] <- as.numeric(one$runtime %||% NA_real_)
gap[r] <- as.numeric(one$gap %||% NA_real_)
alias_values <- stats::setNames(
vapply(
aliases,
function(a) .pamo_eval_alias_on_solution(x, one$solution, a),
numeric(1)
),
aliases
)
value_mat[r, ] <- unname(as.numeric(alias_values))
if (!is.null(solutions[[r]]) && inherits(solutions[[r]], "Solution")) {
solutions[[r]]$solution$alias_values <- alias_values
solutions[[r]]$meta$run_id <- design_df$run_id[r]
solutions[[r]]$method$type <- "weighted"
solutions[[r]]$method$weights <- stats::setNames(as.numeric(w_r), aliases)
}
}
runs <- data.frame(
run_id = design_df$run_id,
status = status,
runtime = runtime,
gap = gap,
stringsAsFactors = FALSE
)
runs <- cbind(
runs,
as.data.frame(value_mat, stringsAsFactors = FALSE)
)
summary_set <- .pamo_bind_solution_summaries(
solutions = solutions,
run_ids = design_df$run_id
)
solutions <- .pa_strip_nested_solution_problems(solutions)
.pamo_cli_done("Weighted solve finished.", verbose = verbose)
out <- pproto(
NULL, SolutionSet,
problem = x,
solution = list(
design = design_df,
runs = runs,
solutions = solutions
),
summary = summary_set,
diagnostics = list(
n_design = nrow(design_df),
n_runs = nrow(runs),
n_solutions = length(solutions),
status_summary = .pa_solutionset_status_summary(runs),
runtime_range = if ("runtime" %in% names(runs)) .pa_solutionset_range_text(runs$runtime, digits = 3) else "none",
gap_range = if ("gap" %in% names(runs)) .pa_solutionset_range_text(runs$gap, digits = 4) else "none"
),
method = method,
meta = list(
call = match.call()
),
name = "solset"
)
.pa_finalize_solution_ids(out)
}
.pamo_solve_epsilon_constraint <- function(x, ...) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% NULL
if (is.null(method) || !is.list(method)) {
stop("epsilon_constraint: x$data$method is missing or invalid.", call. = FALSE)
}
data <- x$data %||% list()
primary <- as.character(method$primary %||% NA_character_)[1]
aliases <- as.character(method$aliases %||% character(0))
constrained <- as.character(method$constrained %||% character(0))
# mode <- as.character(method$mode %||% "manual")[1]
# pass-through solver controls / progress / verbosity
dots <- list(...)
gap_limit <- dots$gap_limit %||% NULL
time_limit <- dots$time_limit %||% NULL
verbose <- .pamo_cli_enabled(x, dots)
# progress control
progress <- dots$progress %||% dots$verbose %||% interactive()
progress <- isTRUE(progress)
.pamo_cli_header("Epsilon-constraint method", verbose = verbose)
if (is.na(primary) || !nzchar(primary)) {
stop("epsilon_constraint: missing primary objective.", call. = FALSE)
}
if (length(aliases) == 0L) {
stop("epsilon_constraint: aliases are missing.", call. = FALSE)
}
if (length(constrained) == 0L) {
stop("epsilon_constraint: no constrained objectives were defined.", call. = FALSE)
}
# if (!mode %in% c("manual", "auto")) {
# stop("epsilon_constraint: unknown mode '", mode, "'.", call. = FALSE)
# }
.pamo_cli_step(
"Primary objective: {.val {primary}}. Secondary objectives: {.val {paste(constrained, collapse = ', ')}}.",
verbose = verbose
)
# ---------------------------------------------------------
# Build or recover epsilon design
# ---------------------------------------------------------
# if (identical(mode, "manual")) {
# .pamo_cli_step("Using user-supplied epsilon grid.", verbose = verbose)
#
# design_df <- method$runs
#
# if (is.null(design_df) || !inherits(design_df, "data.frame") || nrow(design_df) == 0L) {
# stop("epsilon_constraint (manual): x$data$method$runs is missing/empty.", call. = FALSE)
# }
#
# } else {
# .pamo_cli_step("Building automatic epsilon grid.", verbose = verbose)
#
# if (!exists(".pamo_build_auto_epsilon_runs", mode = "function")) {
# stop(
# "epsilon_constraint (auto): missing internal helper .pamo_build_auto_epsilon_runs().",
# call. = FALSE
# )
# }
#
# design_df <- method$runs %||% NULL
# if (is.null(design_df)) {
# design_df <- .pamo_build_auto_epsilon_runs(
# x,
# gap_limit = gap_limit,
# time_limit = time_limit,
# verbose = verbose
# )
# }
#
# if (is.null(design_df) || !inherits(design_df, "data.frame") || nrow(design_df) == 0L) {
# stop("epsilon_constraint (auto): generated epsilon runs are empty.", call. = FALSE)
# }
#
# x$data$method$runs <- design_df
# }
design_df <- .pamo_compile_epsilon_runs(
x = x,
method = method,
gap_limit = gap_limit,
time_limit = time_limit,
verbose = verbose
)
if (!("run_id" %in% names(design_df))) {
design_df$run_id <- seq_len(nrow(design_df))
}
eps_cols <- paste0("eps_", constrained)
miss_cols <- setdiff(eps_cols, names(design_df))
if (length(miss_cols) > 0L) {
stop(
"epsilon_constraint: design is missing epsilon columns: ",
paste(miss_cols, collapse = ", "),
call. = FALSE
)
}
.pamo_cli_step(
"Grid ready with {.val {nrow(design_df)}} run{?s}.",
verbose = verbose
)
# ---------------------------------------------------------
# Solve runs
# ---------------------------------------------------------
n_runs <- nrow(design_df)
solutions <- vector("list", n_runs)
value_mat <- matrix(NA_real_, n_runs, length(aliases))
colnames(value_mat) <- paste0("value_", aliases)
status <- character(n_runs)
runtime <- numeric(n_runs)
gap <- numeric(n_runs)
messages <- character(n_runs)
control <- method$control %||% list()
stop_on_infeasible <- isTRUE(
control$stop_on_infeasible %||%
method$stop_on_infeasible %||%
FALSE
)
stop_on_no_solution <- isTRUE(
control$stop_on_no_solution %||%
method$stop_on_no_solution %||%
FALSE
)
stop_on_error <- isTRUE(
control$stop_on_error %||%
method$stop_on_error %||%
TRUE
)
progress_id <- NULL
if (progress && n_runs > 1L) {
progress_id <- NULL
if (progress && n_runs > 1L) {
progress_id <- cli::cli_progress_bar(
name = "Solving runs",
total = n_runs,
clear = FALSE,
format = "{cli::pb_name}{cli::pb_bar} {cli::pb_current}/{cli::pb_total} | ETA: {cli::pb_eta}"
)
}
}
for (r in seq_len(n_runs)) {
eps_r <- as.list(design_df[r, eps_cols, drop = FALSE])
names(eps_r) <- constrained
one <- .pamo_solve_one_or_keep_failed(
x = x,
spec = list(
type = "epsilon_constraint",
primary = primary,
eps = eps_r,
gap_limit = gap_limit,
time_limit = time_limit
),
stop_on_infeasible = stop_on_infeasible,
stop_on_no_solution = stop_on_no_solution,
stop_on_error = stop_on_error
)
solutions[r] <- list(one$solution)
status[r] <- as.character(one$status %||% NA_character_)
runtime[r] <- as.numeric(one$runtime %||% NA_real_)
gap[r] <- as.numeric(one$gap %||% NA_real_)
if (!is.null(solutions[[r]]) && inherits(solutions[[r]], "Solution")) {
alias_values <- stats::setNames(
vapply(
aliases,
function(a) .pamo_eval_alias_on_solution(x, one$solution, a),
numeric(1)
),
aliases
)
value_mat[r, ] <- unname(as.numeric(alias_values))
solutions[[r]]$solution$alias_values <- alias_values
solutions[[r]]$meta$run_id <- design_df$run_id[r]
solutions[[r]]$method$type <- "epsilon_constraint"
solutions[[r]]$method$primary_alias <- primary
solutions[[r]]$method$eps <- stats::setNames(
as.numeric(design_df[r, eps_cols, drop = TRUE]),
constrained
)
} else {
value_mat[r, ] <- NA_real_
}
if (!is.null(progress_id)) {
cli::cli_progress_update(id = progress_id, set = r)
}
}
if (!is.null(progress_id)) {
cli::cli_progress_done(id = progress_id)
}
# ---------------------------------------------------------
# Runs summary: outputs only
# ---------------------------------------------------------
runs <- data.frame(
run_id = design_df$run_id,
status = status,
runtime = runtime,
gap = gap,
message = messages,
stringsAsFactors = FALSE
)
runs <- cbind(
runs,
as.data.frame(value_mat, stringsAsFactors = FALSE)
)
n_infeasible <- sum(
runs$status %in% c("infeasible", "infeasible_or_unbounded"),
na.rm = TRUE
)
n_no_solution <- sum(
runs$status %in% c("no_solution"),
na.rm = TRUE
)
n_failed <- sum(
runs$status %in% c("failed"),
na.rm = TRUE
)
if (n_infeasible > 0L) {
warning(
n_infeasible, " of ", nrow(runs),
" epsilon-constraint runs were infeasible and were kept in the SolutionSet with missing objective values. ",
"Feasible runs were retained.",
call. = FALSE,
immediate. = TRUE
)
}
if (n_no_solution > 0L) {
warning(
n_no_solution, " of ", nrow(runs),
" epsilon-constraint runs did not return a solution vector and were kept in the SolutionSet with missing objective values. ",
"This may indicate infeasibility, numerical issues, or solver termination before finding a feasible solution.",
call. = FALSE,
immediate. = TRUE
)
}
if (n_failed > 0L) {
warning(
n_failed, " of ", nrow(runs),
" epsilon-constraint runs failed for reasons other than infeasibility/no-solution.",
call. = FALSE,
immediate. = TRUE
)
}
# ---------------------------------------------------------
# Extras: method-specific metadata
# ---------------------------------------------------------
extras <- list()
extras$epsilon_grid <- design_df
ext <- attr(design_df, "extremes", exact = TRUE)
if (!is.null(ext)) {
extras$extremes <- ext
}
if (!is.null(method$meta)) {
extras$method_meta <- method$meta
}
summary_set <- .pamo_bind_solution_summaries(
solutions = solutions,
run_ids = design_df$run_id
)
solutions <- .pa_strip_nested_solution_problems(solutions)
.pamo_cli_done("Epsilon-constraint solve finished.", verbose = verbose)
ext <- attr(design_df, "extremes", exact = TRUE)
ranges_df <- NULL
if (!is.null(ext)) {
primary_vals <- c(
as.numeric(ext$primary_at_primary %||% NA_real_),
as.numeric(ext$primary_at_secondary %||% NA_real_)
)
primary_vals <- primary_vals[is.finite(primary_vals)]
secondary_vals <- c(
as.numeric(ext$secondary_min %||% NA_real_),
as.numeric(ext$secondary_max %||% NA_real_)
)
secondary_vals <- secondary_vals[is.finite(secondary_vals)]
ranges_df <- data.frame(
alias = c(primary, constrained[1]),
role = c("primary", "secondary"),
min = c(
if (length(primary_vals)) min(primary_vals) else NA_real_,
if (length(secondary_vals)) min(secondary_vals) else NA_real_
),
max = c(
if (length(primary_vals)) max(primary_vals) else NA_real_,
if (length(secondary_vals)) max(secondary_vals) else NA_real_
),
source = c("anchors", "auto_grid"),
stringsAsFactors = FALSE
)
}
out <- pproto(
NULL, SolutionSet,
problem = x,
solution = list(
design = design_df,
runs = runs,
solutions = solutions,
frontier = list(
primary = primary,
secondary = constrained,
ranges = ranges_df,
extremes = ext
)
),
summary = summary_set,
diagnostics = list(
n_design = nrow(design_df),
n_runs = nrow(runs),
n_solutions = sum(vapply(solutions, inherits, logical(1), "Solution")),
n_optimal = sum(runs$status == "optimal", na.rm = TRUE),
n_infeasible = sum(runs$status %in% c("infeasible", "infeasible_or_unbounded"), na.rm = TRUE),
n_no_solution = sum(runs$status == "no_solution", na.rm = TRUE),
n_failed = sum(runs$status == "failed", na.rm = TRUE),
status_summary = .pa_solutionset_status_summary(runs),
runtime_range = if ("runtime" %in% names(runs)) .pa_solutionset_range_text(runs$runtime, digits = 3) else "none",
gap_range = if ("gap" %in% names(runs)) .pa_solutionset_range_text(runs$gap, digits = 4) else "none"
),
method = method,
meta = list(
call = match.call()
),
name = "solset"
)
.pa_finalize_solution_ids(out)
}
.pamo_build_auto_epsilon_runs <- function(x,
gap_limit = NULL,
time_limit = NULL,
verbose = FALSE) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% list()
primary <- as.character(method$primary %||% NA_character_)[1]
aliases <- as.character(method$aliases %||% character(0))
constrained <- as.character(method$constrained %||% character(0))
n_points <- as.integer(method$n_points %||% NA_integer_)[1]
include_extremes <- isTRUE(method$include_extremes)
if (is.na(primary) || !nzchar(primary)) {
stop("Auto epsilon mode: missing primary objective.", call. = FALSE)
}
if (length(aliases) == 0L) {
stop("Auto epsilon mode: aliases are missing.", call. = FALSE)
}
if (!primary %in% aliases) {
stop("Auto epsilon mode: primary must be included in aliases.", call. = FALSE)
}
if (length(constrained) == 0L) {
stop("Auto epsilon mode requires at least one constrained objective.", call. = FALSE)
}
if (!is.finite(n_points) || is.na(n_points) || n_points < 2L) {
stop("Auto epsilon mode requires n_points >= 2.", call. = FALSE)
}
if (length(aliases) != 2L) {
stop(
"Auto epsilon mode currently supports exactly 2 objectives.\n",
"You provided ", length(aliases), " objective(s): ",
paste(aliases, collapse = ", "), ".\n",
"Use mode='manual' for 3+ objectives.",
call. = FALSE
)
}
if (length(constrained) != 1L) {
stop(
"Auto epsilon mode currently expects exactly 1 constrained objective.",
call. = FALSE
)
}
secondary <- constrained[1]
.pamo_cli_step(
"Computing automatic epsilon bounds for {.val {secondary}} against {.val {primary}}.",
verbose = verbose
)
ext <- .pamo_compute_epsilon_extremes_2obj(
x = x,
primary = primary,
secondary = secondary,
gap_limit = gap_limit,
time_limit = time_limit
)
sec_min <- as.numeric(ext$secondary_min)[1]
sec_max <- as.numeric(ext$secondary_max)[1]
if (!is.finite(sec_min) || !is.finite(sec_max)) {
stop("Auto epsilon mode: computed epsilon bounds are not finite.", call. = FALSE)
}
if (sec_min > sec_max) {
tmp <- sec_min
sec_min <- sec_max
sec_max <- tmp
}
.pamo_cli_step(
"Secondary range: [{.val {format(sec_min, digits = 6)}}, {.val {format(sec_max, digits = 6)}}].",
verbose = verbose
)
eps_vals <- seq(from = sec_min, to = sec_max, length.out = n_points)
if (!isTRUE(include_extremes)) {
if (length(eps_vals) <= 2L) {
stop(
"include_extremes=FALSE requires n_points >= 3.",
call. = FALSE
)
}
eps_vals <- eps_vals[2:(length(eps_vals) - 1L)]
}
if (length(eps_vals) == 0L) {
stop("Auto epsilon mode produced an empty epsilon grid.", call. = FALSE)
}
eps_col <- paste0("eps_", secondary)
design_df <- data.frame(
run_id = seq_along(eps_vals),
stringsAsFactors = FALSE
)
design_df[[eps_col]] <- as.numeric(eps_vals)
attr(design_df, "extremes") <- ext
attr(design_df, "primary_alias") <- primary
attr(design_df, "secondary_alias") <- secondary
attr(design_df, "include_extremes") <- include_extremes
attr(design_df, "n_points_requested") <- n_points
.pamo_cli_done(
"Automatic epsilon grid built with {.val {nrow(design_df)}} run{?s}.",
verbose = verbose
)
design_df
}
# -------------------------------------------------------------------------
# "Custom objectives" (DEV/advanced) constructor + registry add
# (This keeps the path you started, but stores in base$data$objectives)
# -------------------------------------------------------------------------
# .pamo_objective <- function(alias, sense = c("min", "max"), build, eval, meta = list()) {
# sense <- match.arg(sense)
# if (!is.character(alias) || length(alias) != 1 || !nzchar(alias)) {
# stop("objective `alias` must be a non-empty string.", call. = FALSE)
# }
# if (!is.function(build)) stop("objective `build` must be a function.", call. = FALSE)
# if (!is.function(eval)) stop("objective `eval` must be a function.", call. = FALSE)
#
# structure(
# list(alias = alias, sense = sense, build = build, eval = eval, meta = meta),
# class = "pa_objective"
# )
# }
#
# add_objective <- function(x, objective) {
# #x <- .pamo_as_mo(x)
# x <- .pa_clone_data(x)
#
# if (!inherits(objective, "pa_objective")) {
# stop("add_objective() expects an objective of class 'pa_objective'.", call. = FALSE)
# }
#
# # store it as an atomic spec in the multiscape registry
# if (is.null(x$data$objectives) || !is.list(x$data$objectives)) {
# x$data$objectives <- list()
# }
# if (!is.null(x$data$objectives[[objective$alias]])) {
# stop("Objective alias already exists: '", objective$alias, "'.", call. = FALSE)
# }
#
# x$data$objectives[[objective$alias]] <- list(
# alias = objective$alias,
# objective_id = "custom",
# model_type = "custom",
# objective_args = list(
# build = objective$build,
# eval = objective$eval,
# meta = objective$meta %||% list()
# ),
# sense = objective$sense
# )
#
# x
# }
.pamo_objvec_from_ir <- function(base_superset, ir) {
stopifnot(inherits(base_superset, "Problem"))
op <- base_superset$data$model_ptr
if (is.null(op)) {
stop("Superset model is not built (missing model_ptr).", call. = FALSE)
}
# ------------------------------------------------------------
# local helpers
# ------------------------------------------------------------
.chr <- function(z) {
if (is.null(z)) return(NULL)
z <- as.character(z)
z <- unique(z[!is.na(z) & nzchar(z)])
if (length(z) == 0) return(NULL)
z
}
.resolve_feature_internal <- function(features) {
if (is.null(features)) return(integer(0))
feat_df <- base_superset$data$features
if (is.null(feat_df) || !inherits(feat_df, "data.frame") || nrow(feat_df) == 0) {
stop("Cannot resolve feature subset: x$data$features is missing/empty.", call. = FALSE)
}
if (!all(c("id", "internal_id") %in% names(feat_df))) {
stop("x$data$features must contain 'id' and 'internal_id'.", call. = FALSE)
}
vals_chr <- .chr(features)
vals_int <- suppressWarnings(as.integer(vals_chr))
# 1) try internal ids first
if (length(vals_int) > 0 && all(!is.na(vals_int)) &&
all(vals_int %in% feat_df$internal_id)) {
return(as.integer(unique(vals_int)))
}
# 2) then external ids
idx <- match(vals_chr, as.character(feat_df$id))
if (all(!is.na(idx))) {
return(as.integer(feat_df$internal_id[idx]))
}
# 3) then names if available
if ("name" %in% names(feat_df)) {
idx2 <- match(vals_chr, as.character(feat_df$name))
if (all(!is.na(idx2))) {
return(as.integer(feat_df$internal_id[idx2]))
}
}
bad <- vals_chr[
!(vals_chr %in% as.character(feat_df$id)) &
!(vals_chr %in% as.character(feat_df$internal_id)) &
!("name" %in% names(feat_df) && vals_chr %in% as.character(feat_df$name))
]
stop(
"Unknown features in objective subset: ",
paste(unique(bad), collapse = ", "),
call. = FALSE
)
}
.resolve_action_internal <- function(actions) {
if (is.null(actions)) return(integer(0))
act_df <- base_superset$data$actions
if (is.null(act_df) || !inherits(act_df, "data.frame") || nrow(act_df) == 0) {
stop("Cannot resolve action subset: x$data$actions is missing/empty.", call. = FALSE)
}
if (!all(c("id", "internal_id") %in% names(act_df))) {
stop("x$data$actions must contain 'id' and 'internal_id'.", call. = FALSE)
}
vals_chr <- .chr(actions)
vals_int <- suppressWarnings(as.integer(vals_chr))
# 1) try internal ids first
if (length(vals_int) > 0 && all(!is.na(vals_int)) &&
all(vals_int %in% act_df$internal_id)) {
return(as.integer(unique(vals_int)))
}
# 2) otherwise fall back to the standard resolver (ids / action_set / etc.)
out <- .pa_resolve_action_subset(base_superset, subset = vals_chr)
as.integer(out$internal_id)
}
.subset_dist_actions <- function(actions = NULL) {
da <- base_superset$data$dist_actions_model
if (is.null(da) || !inherits(da, "data.frame")) {
stop("dist_actions_model is missing.", call. = FALSE)
}
act_int <- .resolve_action_internal(actions)
if (length(act_int) == 0) return(da)
da[da$internal_action %in% act_int, , drop = FALSE]
}
.subset_dist_effects <- function(df, actions = NULL, features = NULL) {
if (is.null(df) || !inherits(df, "data.frame")) return(df)
act_int <- .resolve_action_internal(actions)
feat_int <- .resolve_feature_internal(features)
out <- df
if (length(act_int) > 0 && "internal_action" %in% names(out)) {
out <- out[out$internal_action %in% act_int, , drop = FALSE]
}
if (length(feat_int) > 0 && "internal_feature" %in% names(out)) {
out <- out[out$internal_feature %in% feat_int, , drop = FALSE]
}
out
}
.subset_dist_profit <- function(df, actions = NULL) {
if (is.null(df) || !inherits(df, "data.frame")) return(df)
act_int <- .resolve_action_internal(actions)
if (length(act_int) == 0) return(df)
df[df$internal_action %in% act_int, , drop = FALSE]
}
.subset_dist_features <- function(df, features = NULL) {
if (is.null(df) || !inherits(df, "data.frame")) return(df)
feat_int <- .resolve_feature_internal(features)
if (length(feat_int) == 0) return(df)
df[df$internal_feature %in% feat_int, , drop = FALSE]
}
# Convención MO: motor en MIN; luego en R se invierte signo para objetivos max
rcpp_reset_objective(op, "min")
terms <- ir$terms %||% list()
for (t in terms) {
type <- t$type
if (identical(type, "pu_cost")) {
feats <- .chr(t$features)
if (!is.null(feats)) {
stop(
"Subset by features is not currently supported for 'pu_cost'. ",
"PU costs are attached to w variables, not directly to feature-specific contributions.",
call. = FALSE
)
}
rcpp_add_objective_min_cost(
op,
pu_data = base_superset$data$pu,
dist_actions_data = base_superset$data$dist_actions_model,
include_pu_cost = TRUE,
include_action_cost = FALSE,
weight = 1.0
)
} else if (identical(type, "action_cost")) {
da_sub <- .subset_dist_actions(actions = t$actions)
rcpp_add_objective_min_cost(
op,
pu_data = base_superset$data$pu,
dist_actions_data = da_sub,
include_pu_cost = FALSE,
include_action_cost = TRUE,
weight = 1.0
)
} else if (identical(type, "boundary_cut")) {
feats <- .chr(t$features)
if (!is.null(feats)) {
stop(
"Subset by features is not currently supported for 'boundary_cut'. ",
"PU fragmentation is defined over selected units, not over feature-specific selections.",
call. = FALSE
)
}
rel_name <- as.character(t$relation_name %||% "boundary")[1]
rel_model <- base_superset$data$spatial_relations_model[[rel_name]] %||%
base_superset$data$spatial_relations[[rel_name]]
if (is.null(rel_model)) {
stop("Missing relation '", rel_name, "'.", call. = FALSE)
}
rcpp_prepare_objective_min_fragmentation(op, rel_model)
rcpp_add_objective_min_fragmentation(
op,
relation_data = rel_model,
weight_multiplier = as.numeric(t$weight_multiplier %||% 1)[1]
)
} else if (identical(type, "benefit")) {
if (.pa_uses_aggregate_effects(base_superset)) {
current <- rcpp_optimization_problem_as_list(op)$obj
vector <- .pa_joint_objective_vector(base_superset, t$benefit_col %||% "effect", t$actions, t$features)
rcpp_model_set_objective_vector(op, as.numeric(current) + vector, "min")
next
}
de_sub <- .subset_dist_effects(
df = base_superset$data$dist_effects_model,
actions = t$actions,
features = t$features
)
da_sub <- .subset_dist_actions(actions = t$actions)
if (is.null(de_sub) || !inherits(de_sub, "data.frame")) {
stop("dist_effects_model is missing for benefit objective.", call. = FALSE)
}
for (nm in c("internal_pu", "internal_action", "benefit")) {
if (!(nm %in% names(de_sub))) {
stop("dist_effects_model must contain column '", nm, "' for benefit objective.", call. = FALSE)
}
}
# Canonical positive/loss components remain unchanged for reporting.
bcol <- t$benefit_col %||% "effect"
if (!bcol %in% names(de_sub)) stop("Missing benefit coefficient column '", bcol, "'.", call. = FALSE)
de_sub$benefit <- as.numeric(de_sub[[bcol]])
prep <- rcpp_prepare_objective_max_benefit(
x = op,
dist_actions_data = da_sub,
dist_effects_data = de_sub,
block_name = "prepare_max_benefit",
tag = "from=.pamo_objvec_from_ir"
)
rcpp_add_objective_max_benefit(
x = op,
coef_x = as.numeric(prep$coef_x),
weight = 1.0,
block_name = "objective_add_max_benefit",
tag = "from=.pamo_objvec_from_ir"
)
} else if (identical(type, "loss")) {
if (.pa_uses_aggregate_effects(base_superset)) {
current <- rcpp_optimization_problem_as_list(op)$obj
vector <- .pa_joint_objective_vector(base_superset, "loss", t$actions, t$features)
rcpp_model_set_objective_vector(op, as.numeric(current) + vector, "min")
next
}
de_sub <- .subset_dist_effects(
df = base_superset$data$dist_effects_model,
actions = t$actions,
features = t$features
)
da_sub <- .subset_dist_actions(actions = t$actions)
if (is.null(de_sub) || !inherits(de_sub, "data.frame")) {
stop("dist_effects_model is missing for loss objective.", call. = FALSE)
}
for (nm in c("internal_pu", "internal_action", "loss")) {
if (!(nm %in% names(de_sub))) {
stop("dist_effects_model must contain column '", nm, "' for loss objective.", call. = FALSE)
}
}
prep <- rcpp_prepare_objective_min_loss(
x = op,
dist_actions_data = da_sub,
dist_effects_data = de_sub,
block_name = "prepare_min_loss",
tag = "from=.pamo_objvec_from_ir"
)
rcpp_add_objective_min_loss(
x = op,
coef_x = as.numeric(prep$coef_x),
weight = 1.0,
block_name = "objective_add_min_loss",
tag = "from=.pamo_objvec_from_ir"
)
} else if (identical(type, "profit")) {
pcol <- as.character(t$profit_col %||% "profit")[1]
dp_sub <- .subset_dist_profit(
df = base_superset$data$dist_profit_model,
actions = t$actions
)
da_sub <- .subset_dist_actions(actions = t$actions)
rcpp_add_objective_max_profit(
op,
dist_actions_data = da_sub,
dist_profit_data = dp_sub,
profit_col = pcol,
weight = 1.0
)
} else if (identical(type, "net_profit")) {
pcol <- as.character(t$profit_col %||% "profit")[1]
inc_pu <- isTRUE(t$include_pu_cost %||% TRUE)
inc_act <- isTRUE(t$include_action_cost %||% TRUE)
feats <- .chr(t$features)
if (!is.null(feats) && isTRUE(inc_pu)) {
stop(
"Subset by features is not currently supported for 'net_profit' when include_pu_cost=TRUE. ",
"PU costs are attached to w variables and cannot yet be attributed to feature subsets.",
call. = FALSE
)
}
da_sub <- .subset_dist_actions(actions = t$actions)
dp_sub <- .subset_dist_profit(
df = base_superset$data$dist_profit_model,
actions = t$actions
)
rcpp_add_objective_max_net_profit(
op,
pu_data = base_superset$data$pu,
dist_actions_data = da_sub,
dist_profit_data = dp_sub,
profit_col = pcol,
include_pu_cost = inc_pu,
include_action_cost = inc_act,
weight = 1.0
)
} else if (identical(type, "action_boundary_cut")) {
rel_name <- as.character(t$relation_name %||% "boundary")[1]
rel_model <- base_superset$data$spatial_relations_model[[rel_name]] %||%
base_superset$data$spatial_relations[[rel_name]]
if (is.null(rel_model)) {
stop("Missing relation '", rel_name, "'.", call. = FALSE)
}
act_int <- .resolve_action_internal(t$actions)
rcpp_prepare_objective_min_fragmentation_actions(
op,
dist_actions_data = base_superset$data$dist_actions_model,
relation_data = rel_model,
actions_to_use = if (length(act_int) > 0) act_int else NULL,
block_name = "fragmentation_actions",
tag = "from=.pamo_objvec_from_ir"
)
term2 <- t
term2$actions <- if (length(act_int) > 0) as.integer(act_int) else NULL
v_term <- .pamo_objvec_action_boundary_cut(base_superset, term2)
if (length(v_term) == 0L) {
stop("Empty vector returned by .pamo_objvec_action_boundary_cut().", call. = FALSE)
}
if (!exists("rcpp_model_set_objective_vector", mode = "function")) {
m0 <- .pa_model_from_ptr(
op,
args = base_superset$data$model_args %||% list(),
drop_triplets = TRUE
)
obj0 <- as.numeric(m0$obj)
if (length(obj0) != length(v_term)) {
stop("Current model objective length does not match action_boundary_cut term vector.", call. = FALSE)
}
rcpp_model_set_objective_vector(
x = op,
obj = obj0 + v_term,
model_sense = "min"
)
} else {
m0 <- .pa_model_from_ptr(
op,
args = base_superset$data$model_args %||% list(),
drop_triplets = TRUE
)
obj0 <- as.numeric(m0$obj)
if (length(obj0) != length(v_term)) {
stop("Current model objective length does not match action_boundary_cut term vector.", call. = FALSE)
}
rcpp_model_set_objective_vector(
x = op,
obj = obj0 + v_term,
model_sense = "min"
)
}
} else if (identical(type, "intervention_impact")) {
icol <- as.character(t$impact_col %||% "amount")[1]
df_sub <- .subset_dist_features(
df = base_superset$data$dist_features,
features = t$features
)
feat_int <- .resolve_feature_internal(t$features)
act_int <- .resolve_action_internal(t$actions)
if (length(act_int) == 0) {
stop(
"intervention_impact requires a non-empty action subset.",
call. = FALSE
)
}
subset_key <- .pa_subset_to_string(as.character(t$actions))
rcpp_add_objective_min_intervention_impact(
op,
pu_data = base_superset$data$pu,
dist_features_data = df_sub,
dist_actions_data = base_superset$data$dist_actions_model,
subset_key = subset_key,
impact_col = icol,
features_to_use = as.integer(feat_int),
actions_to_use = as.integer(act_int),
internal_feature_col = "internal_feature",
weight = 1.0,
weight_multiplier = as.numeric(t$weight_multiplier %||% 1)[1],
block_name = "objective_add_min_intervention_impact",
tag = "from=.pamo_objvec_from_ir"
)
} else if (identical(type, "custom")) {
stop("custom objectives are not supported in weighted yet.", call. = FALSE)
} else {
stop("Unknown term type in .pamo_objvec_from_ir(): ", type, call. = FALSE)
}
}
m <- .pa_model_from_ptr(
op,
args = base_superset$data$model_args %||% list(),
drop_triplets = TRUE
)
as.numeric(m$obj)
}
.pamo_apply_weighted_objective <- function(base,
ir_list,
weights,
objective_scaling = FALSE,
scales = NULL) {
stopifnot(inherits(base, "Problem"))
stopifnot(is.list(ir_list), length(ir_list) > 0)
weights <- as.numeric(weights)
if (length(weights) != length(ir_list)) {
stop("weights length must match ir_list length.", call. = FALSE)
}
# 1) construir superset UNA vez (incluye prepares necesarios)
base2 <- .pamo_prepare_superset_model(base, ir_list)
# 2) construir objvec para cada IR sobre el MISMO superset (sin rebuild)
objvecs <- vector("list", length(ir_list))
for (i in seq_along(ir_list)) {
v <- .pamo_objvec_from_ir(base2, ir_list[[i]])
# convención: todo se resuelve como MIN en el motor;
# si IR es "max", flipa signo acá.
if (identical(ir_list[[i]]$sense, "max")) v <- -v
objvecs[[i]] <- v
}
# 3) objective scaling by payoff range
if (isTRUE(objective_scaling)) {
if (is.null(scales) || length(scales) != length(objvecs)) {
stop(
"objective_scaling = TRUE requires one positive scale per objective.",
call. = FALSE
)
}
objvecs <- Map(
function(v, s) {
s <- as.numeric(s)[1]
if (!is.finite(s) || s <= 0) {
warning("Encountered non-positive scale in weighted objective scaling; leaving objective unscaled.",
call. = FALSE)
return(v)
}
v / s
},
objvecs,
as.list(scales)
)
}
# 4) combinar pesos -> objetivo final
obj_w <- Reduce(`+`, Map(`*`, objvecs, as.list(weights)))
# 5) IMPORTANTÃSIMO: inyectar runtime update para que multiscape lo use en solve()
base2$data$runtime_updates <- list(
obj = obj_w,
modelsense = "min"
)
# 6) asegurar que multiscape NO reconstruya el modelo y no te borre runtime_updates
base2$data$meta$model_dirty <- FALSE
base2$data$has_model <- TRUE
# 7) guarda cache para evaluación posterior (opcional pero muy útil)
base2$data$mo_cache <- list(
ir_list = ir_list,
weights = weights,
objective_scaling = objective_scaling,
scales = scales,
objvecs = objvecs,
obj_weighted = obj_w
)
base2
}
# ---------------------------------------------------------
# Internal: solve a single run using multiscape as engine
# ---------------------------------------------------------
.pamo_solve_one <- function(x, spec) {
if (!inherits(x, "Problem")) stop("Problem$base must be a Problem.", call. = FALSE)
base <- .pamo_clone_base(x)
if (identical(spec$type, "weighted")) {
aliases <- as.character(spec$aliases)
weights <- as.numeric(spec$weights)
objective_scaling <- isTRUE(spec$objective_scaling)
scales <- spec$scales %||% NULL
specs <- .pamo_get_objective_specs(x, aliases)
ir_list <- lapply(specs, function(sp) .pamo_objective_to_ir(base, sp))
base <- .pamo_apply_weighted_objective(
base = base,
ir_list = ir_list,
weights = weights,
objective_scaling = objective_scaling,
scales = scales
)
} else if (identical(spec$type, "epsilon_constraint")) {
primary <- as.character(spec$primary)[1]
eps_list <- spec$eps %||% list()
eps_tol <- spec$eps_tol %||% list()
sec_aliases <- names(eps_list)
if (is.na(primary) || !nzchar(primary)) stop("epsilon_constraint: primary is invalid.", call. = FALSE)
if (length(sec_aliases) == 0) stop("epsilon_constraint: eps list is empty.", call. = FALSE)
# ---- build IRs for (primary + constrained) so we can build ONE superset
all_aliases <- unique(c(primary, sec_aliases))
specs_all <- .pamo_get_objective_specs(x, all_aliases)
ir_all <- lapply(specs_all, function(sp) .pamo_objective_to_ir(base, sp))
# ---- 0) build superset once (includes prepares)
base <- .pamo_prepare_superset_model(base, ir_all)
# ---- 1) add epsilon constraints for secondaries
for (a in sec_aliases) {
sp_sec <- specs_all[[a]]
ir_sec <- ir_all[[a]]
eps_val <- as.numeric(eps_list[[a]])[1]
if (!is.finite(eps_val)) {
stop("epsilon_constraint: eps for '", a, "' must be finite.", call. = FALSE)
}
tol_a <- as.numeric(eps_tol[[a]] %||% 0)[1]
if (!is.finite(tol_a) || tol_a < 0) {
stop("epsilon_constraint: eps_tol for '", a, "' must be finite and >= 0.", call. = FALSE)
}
sec_sense <- as.character(sp_sec$sense %||% "min")[1]
eps_relaxed <- if (identical(sec_sense, "max")) {
eps_val - tol_a
} else {
eps_val + tol_a
}
base <- .pamo_apply_epsilon_constraint(
base = base,
ir = ir_sec,
eps = eps_relaxed,
sense = sec_sense,
name = paste0("eps_", a),
block_name = "epsilon_constraint"
)
}
# ---- 2) set primary objective (as runtime update, like weighted)
sp_primary <- specs_all[[primary]]
ir_primary <- ir_all[[primary]]
base <- .pamo_apply_single_objective(
base = base,
ir = ir_primary,
sense = as.character(sp_primary$sense %||% "min")[1]
)
} else if (identical(spec$type, "augmecon")) {
primary <- as.character(spec$primary)[1]
eps_list <- spec$eps %||% list()
eps_tol <- spec$eps_tol %||% list()
sec_aliases <- names(eps_list)
sec_ranges <- as.numeric(spec$secondary_ranges %||% numeric(0))
augmentation <- as.numeric(spec$augmentation %||% NA_real_)[1]
slack_upper_bound <- as.numeric((x$data$method %||% list())$slack_upper_bound %||% 1e6)[1]
if (is.na(primary) || !nzchar(primary)) {
stop("augmecon: primary is invalid.", call. = FALSE)
}
if (length(sec_aliases) == 0L) {
stop("augmecon: eps list is empty.", call. = FALSE)
}
if (!is.finite(augmentation) || augmentation <= 0) {
stop("augmecon: augmentation must be a finite positive number.", call. = FALSE)
}
if (!is.finite(slack_upper_bound) || slack_upper_bound <= 0) {
stop("augmecon: slack_upper_bound must be a finite positive number.", call. = FALSE)
}
if (length(sec_ranges) != length(sec_aliases)) {
stop(
"augmecon: secondary_ranges must have one value per secondary objective.\n",
"length(secondary_ranges) = ", length(sec_ranges),
", length(sec_aliases) = ", length(sec_aliases), ".",
call. = FALSE
)
}
if (any(!is.finite(sec_ranges)) || any(sec_ranges <= 0)) {
stop("augmecon: secondary_ranges must be finite and > 0.", call. = FALSE)
}
names(sec_ranges) <- sec_aliases
# ---- build IRs for primary + secondaries
all_aliases <- unique(c(primary, sec_aliases))
specs_all <- .pamo_get_objective_specs(x, all_aliases)
ir_all <- lapply(specs_all, function(sp) .pamo_objective_to_ir(base, sp))
# ---- build one superset
base <- .pamo_prepare_superset_model(base, ir_all)
# ---- add explicit slack variables (one per secondary objective)
base <- .pamo_add_augmecon_slacks(
base = base,
secondary_aliases = sec_aliases,
ub = slack_upper_bound
)
slack_cols <- base$data$mo_cache$augmecon$slack_cols_0based[sec_aliases]
if (length(slack_cols) != length(sec_aliases)) {
stop("augmecon: failed to recover slack columns for all secondary objectives.", call. = FALSE)
}
# ---- add equality constraints with slack: g_j(x) + s_j = eps_j
for (j in seq_along(sec_aliases)) {
a <- sec_aliases[j]
sp_sec <- specs_all[[a]]
ir_sec <- ir_all[[a]]
eps_val <- as.numeric(eps_list[[a]])[1]
if (!is.finite(eps_val)) {
stop("augmecon: eps for '", a, "' must be finite.", call. = FALSE)
}
tol_a <- as.numeric(eps_tol[[a]] %||% 0)[1]
if (!is.finite(tol_a) || tol_a < 0) {
stop("augmecon: eps_tol for '", a, "' must be finite and >= 0.", call. = FALSE)
}
base <- .pamo_apply_augmecon_constraint(
base = base,
ir = ir_sec,
eps = eps_val + tol_a,
slack_col = as.integer(slack_cols[[a]]),
sense = as.character(sp_sec$sense %||% "min")[1],
name = paste0("augmecon_", a),
block_name = "augmecon_constraint"
)
}
# ---- primary objective + normalized slack rewards
ir_primary <- ir_all[[primary]]
slack_weights <- augmentation / sec_ranges
base <- .pamo_apply_augmecon_objective(
base = base,
ir_primary = ir_primary,
slack_cols = as.integer(unname(slack_cols)),
slack_weights = as.numeric(slack_weights)
)
} else {
stop("Unsupported spec$type in .pamo_solve_one: ", spec$type, call. = FALSE)
}
gap_limit <- spec$gap_limit %||% NULL
time_limit <- spec$time_limit %||% NULL
# base$data$method <- NULL
# base$data$results <- NULL
# base$data$runtime_updates <- NULL
# MO scalar objectives can contain very small augmentation/weight terms
# next to monetary costs. Request numerically precise solves while keeping
# the user's explicit solver parameters authoritative.
base$data$runtime_updates <- base$data$runtime_updates %||% list()
base$data$runtime_updates$mo_numeric <- TRUE
out <- .pa_solve_single_problem(
base,
gap_limit = gap_limit,
time_limit = time_limit
)
if (identical(spec$type, "augmecon") && identical(out$diagnostics$status_code, 0L)) {
out <- .pamo_refine_augmecon(base, out, spec)
}
if (spec$type %in% c("epsilon_constraint", "augmecon")) {
values <- stats::setNames(vapply(names(spec$eps), function(a)
.pamo_eval_alias_on_solution(base, out, a), numeric(1)), names(spec$eps))
senses <- vapply(specs_all[names(values)], function(s) s$sense, character(1))
effective <- out$diagnostics$solver_args$effective_solver_params %||% list()
feasibility_tol <- as.numeric(effective$FeasibilityTol %||%
effective$primalTolerance %||% 1e-7)
out$diagnostics$epsilon_checks <- .pamo_check_secondary_values(
values, senses, spec$eps, spec$eps_tol, feasibility_tol
)
out$solution$alias_values <- values
}
.pamo_extract_solution(out)
}
.pamo_refine_augmecon <- function(base, original, spec) {
cache <- base$data$mo_cache$augmecon
primary <- .pamo_eval_alias_on_solution(base, original, spec$primary)
original_secondary <- stats::setNames(vapply(names(spec$eps), function(a)
.pamo_eval_alias_on_solution(base, original, a), numeric(1)), names(spec$eps))
senses <- vapply(.pamo_get_objective_specs(base, names(spec$eps)),
function(s) s$sense, character(1))
ranges <- as.numeric(spec$secondary_ranges)
levels <- as.numeric(unlist(spec$eps, use.names = FALSE))
relaxation <- vapply(names(spec$eps), function(a)
as.numeric(spec$eps_tol[[a]] %||% 0), numeric(1))
score <- function(p, secondary) {
canon <- if (cache$ir_primary$sense == "max") -p else p
slack <- ifelse(senses == "max", secondary - levels - relaxation,
levels + relaxation - secondary)
canon - spec$augmentation * sum(slack / ranges)
}
original_score <- score(primary, original_secondary)
settings <- .pa_get_solve_args(base, gap_limit = spec$gap_limit, time_limit = spec$time_limit)
effective <- original$diagnostics$solver_args$effective_solver_params %||% list()
limit <- as.numeric(effective$TimeLimit %||% effective$sec %||% settings$time_limit)
remaining <- limit - original$diagnostics$runtime
if (!is.finite(remaining)) remaining <- settings$time_limit
if (remaining <= 0) {
original$diagnostics$augmentation_refinement <- "time_limit"
return(original)
}
# Refine only on the attained primary level. Removing monetary coefficients
# from this pass makes the normalized slack criterion numerically visible.
# The first pass still solves the user's complete augmented scalar objective.
v <- .pamo_objvec_from_ir(base, cache$ir_primary)
nz <- which(v != 0)
rcpp_add_linear_constraint(base$data$model_ptr, as.integer(nz - 1L),
as.numeric(v[nz]), "=", primary, name = "augmecon_primary_fix",
block_name = "augmecon_primary_fix")
base <- .pa_refresh_model_snapshot(base)
reward <- numeric(length(v))
reward[cache$slack_cols_0based + 1L] <- -1 / ranges
base$data$runtime_updates <- list(obj = reward, modelsense = "min", mo_numeric = TRUE,
start = original$solution$vector)
base$data$meta$model_dirty <- FALSE
# Raw solver parameters also obey the remaining per-run budget.
if (!is.null(settings$solver_params$TimeLimit)) base$data$solve_args$solver_params$TimeLimit <- remaining
if (!is.null(settings$solver_params$sec)) base$data$solve_args$solver_params$sec <- as.character(remaining)
t0 <- Sys.time()
refined <- tryCatch(.pa_solve_single_problem(base, gap_limit = spec$gap_limit,
time_limit = remaining), error = function(e) e)
elapsed <- as.numeric(difftime(Sys.time(), t0, units = "secs"))
original$diagnostics$runtime <- original$diagnostics$runtime + elapsed
if (inherits(refined, "error")) {
original$diagnostics$augmentation_refinement <- conditionMessage(refined)
return(original)
}
p <- .pamo_eval_alias_on_solution(base, refined, spec$primary)
secondary <- stats::setNames(vapply(names(spec$eps), function(a)
.pamo_eval_alias_on_solution(base, refined, a), numeric(1)), names(spec$eps))
candidate_score <- score(p, secondary)
precision <- 32 * .Machine$double.eps * max(1, abs(primary), abs(original_score))
if (abs(p - primary) > precision || candidate_score > original_score + precision) {
original$diagnostics$augmentation_refinement <- "original_retained"
return(original)
}
# Retain the first pass's global status/gap and settings. The refinement is
# optimal on a restricted primary level, not a new global proof by itself.
first_diagnostics <- original$diagnostics
first_diagnostics$augmentation_refinement <- list(
status_code = refined$diagnostics$status_code, primary = primary,
original_objective = original_score, refined_objective = candidate_score,
solver_args = refined$diagnostics$solver_args)
refined$diagnostics <- first_diagnostics
refined$solution$objective <- candidate_score
refined
}
.pamo_check_secondary_values <- function(values, senses, eps, eps_tol = NULL,
feasibility_tol = 1e-9) {
aliases <- names(values)
levels <- as.numeric(unlist(eps[aliases], use.names = FALSE))
relaxation <- vapply(aliases, function(a) as.numeric(eps_tol[[a]] %||% 0), numeric(1))
limits <- ifelse(senses == "max", levels - relaxation, levels + relaxation)
violations <- pmax(0, ifelse(senses == "max", limits - values, values - limits))
# Allow floating-point accumulation/cancellation, not the old 1e-6 default
# solver tolerance. Explicit solver tolerances remain authoritative.
tolerance <- pmax(feasibility_tol, 1e-12 * pmax(1, abs(values), abs(limits)))
valid <- is.finite(values) & violations <= tolerance
if (any(!valid)) {
bad <- aliases[!valid][1]
stop("Solver returned a numerically invalid solution: secondary objective '",
bad, "' violates its epsilon bound after evaluating the selected decisions.",
call. = FALSE)
}
data.frame(alias = aliases, sense = unname(senses), value = unname(values),
epsilon = levels, limit = limits, violation = unname(violations),
tolerance = unname(tolerance), valid = unname(valid))
}
.pamo_apply_epsilon_constraint <- function(base, ir, eps, sense = c("min","max"),
name = "", block_name = "epsilon_constraint", tag = "") {
stopifnot(inherits(base, "Problem"))
sense <- match.arg(sense)
if (is.null(base$data$model_ptr)) stop("Model not built (model_ptr is NULL).", call. = FALSE)
eps <- as.numeric(eps)[1]
if (!is.finite(eps)) stop("eps must be finite.", call. = FALSE)
v <- .pamo_objvec_from_ir(base, ir)
if (length(v) == 0) stop("epsilon objvec is empty.", call. = FALSE)
# canonical: <=
if (identical(sense, "max")) {
v <- -v
eps <- -eps
}
idx <- which(v != 0)
# add row: sum(v[j]*x[j]) <= eps
rcpp_add_linear_constraint(
base$data$model_ptr,
j0 = as.integer(idx - 1L),
x = as.numeric(v[idx]),
sense = "<=",
rhs = as.numeric(eps),
name = as.character(name %||% ""),
block_name = block_name,
tag = tag
)
base$data$meta$model_dirty <- TRUE
base <- .pa_refresh_model_snapshot(base)
base
}
# ---------------------------------------------------------
# Internal: extract (minimal) results from multiscape Solution
# ---------------------------------------------------------
.pamo_extract_solution <- function(out) {
objval <- tryCatch(out$solution$objective, error = function(e) NA_real_)
gap <- tryCatch(out$diagnostics$gap, error = function(e) NA_real_)
runtime <- tryCatch(out$diagnostics$runtime, error = function(e) NA_real_)
status <- tryCatch(getStatus(out), error = function(e) "unknown")
list(
solution = out,
status = status,
runtime = runtime,
gap = gap,
objval = objval
)
}
.pamo_solver_error_status <- function(e) {
msg <- conditionMessage(e)
if (grepl(
"status:\\s*infeasible_or_unbounded|\\binfeasible_or_unbounded\\b",
msg,
ignore.case = TRUE
)) {
return("infeasible_or_unbounded")
}
if (grepl(
"status:\\s*infeasible|\\binfeasible\\b",
msg,
ignore.case = TRUE
)) {
return("infeasible")
}
if (grepl(
"empty solution vector|no solution vector|solution vector is empty|objval:\\s*NA|time limit before first feasible",
msg,
ignore.case = TRUE
)) {
return("no_solution")
}
NA_character_
}
.pamo_failed_run <- function(status = "failed", message = NULL) {
list(
solution = NULL,
status = status,
runtime = NA_real_,
gap = NA_real_,
objval = NA_real_,
message = as.character(message %||% NA_character_)
)
}
.pamo_solve_one_or_keep_failed <- function(
x,
spec,
stop_on_infeasible = FALSE,
stop_on_no_solution = FALSE,
stop_on_error = TRUE
) {
tryCatch(
.pamo_solve_one(x = x, spec = spec),
error = function(e) {
msg <- conditionMessage(e)
status <- .pamo_solver_error_status(e)
if (identical(status, "infeasible") ||
identical(status, "infeasible_or_unbounded")) {
if (isTRUE(stop_on_infeasible)) {
stop(e)
}
return(.pamo_failed_run(
status = status,
message = msg
))
}
if (identical(status, "no_solution")) {
if (isTRUE(stop_on_no_solution)) {
stop(e)
}
return(.pamo_failed_run(
status = status,
message = msg
))
}
if (isTRUE(stop_on_error)) {
stop(e)
}
.pamo_failed_run(
status = "failed",
message = msg
)
}
)
}
#' @keywords internal
.mo_abort <- function(...) stop(paste0(...), call. = FALSE)
#' @keywords internal
.mo_get_solution_from <- function(x, run = NULL, solution_id = NULL) {
if (inherits(x, "Solution")) {
return(x)
}
if (!inherits(x, "SolutionSet")) {
.mo_abort("x must be a SolutionSet object returned by solve().")
}
if (!is.null(run) && !is.null(solution_id)) {
.mo_abort("Use either `run` or `solution_id`, not both.")
}
sols <- x$solution$solutions %||% NULL
if (is.null(sols) || !is.list(sols) || length(sols) == 0L) {
.mo_abort("No stored solutions found in this SolutionSet.")
}
is_sol <- vapply(
sols,
inherits,
logical(1),
"Solution"
)
sols <- sols[is_sol]
if (length(sols) == 0L) {
.mo_abort("No stored valid solutions found in this SolutionSet.")
}
# Ensure names exist. New convention: names are character versions of
# numeric solution ids, which match run ids.
if (is.null(names(sols)) || any(!nzchar(names(sols)))) {
ids <- vapply(seq_along(sols), function(i) {
sid <- sols[[i]]$meta$solution_id %||%
sols[[i]]$meta$run_id %||%
i
sid <- as.integer(sid)[1]
if (is.na(sid) || sid < 1L) {
as.character(i)
} else {
as.character(sid)
}
}, character(1))
names(sols) <- ids
}
# ------------------------------------------------------------
# Case 1: selected by solution_id
# ------------------------------------------------------------
if (!is.null(solution_id)) {
sid_int <- suppressWarnings(as.integer(solution_id))[1]
if (is.na(sid_int) || sid_int < 1L) {
.mo_abort("`solution_id` must be a positive integer.")
}
sid_chr <- as.character(sid_int)
if (!sid_chr %in% names(sols)) {
.mo_abort(
"No stored solution found for solution_id = ", sid_int, ". ",
"Available solution ids are: ",
paste(names(sols), collapse = ", "),
"."
)
}
return(sols[[sid_chr]])
}
# ------------------------------------------------------------
# Case 2: selected by run_id
# ------------------------------------------------------------
if (!is.null(run)) {
run <- as.integer(run)[1]
if (!is.finite(run) || is.na(run) || run < 1L) {
.mo_abort("`run` must be a positive integer run id.")
}
runs <- x$solution$runs %||% NULL
if (
!is.null(runs) &&
inherits(runs, "data.frame") &&
all(c("run_id", "solution_id") %in% names(runs))
) {
idx <- match(run, as.integer(runs$run_id))
if (is.na(idx)) {
.mo_abort("No run found with run_id = ", run, ".")
}
sid <- suppressWarnings(as.integer(runs$solution_id[idx]))[1]
if (is.na(sid) || sid < 1L) {
status_msg <- ""
if ("status" %in% names(runs)) {
status_msg <- paste0(" Run status is '", runs$status[idx], "'.")
}
.mo_abort(
"No stored solution is available for run_id = ",
run,
".",
status_msg
)
}
sid_chr <- as.character(sid)
if (!sid_chr %in% names(sols)) {
.mo_abort(
"Run ",
run,
" points to solution_id = ",
sid,
", but that solution is not stored."
)
}
return(sols[[sid_chr]])
}
# Fallback for old objects without a valid run table.
stored_run_ids <- vapply(seq_along(sols), function(i) {
rid <- sols[[i]]$meta$run_id %||% NA_integer_
as.integer(rid)[1]
}, integer(1))
idx <- match(run, stored_run_ids)
if (is.na(idx)) {
.mo_abort("No stored solution is available for run_id = ", run, ".")
}
return(sols[[idx]])
}
# ------------------------------------------------------------
# Case 3: default = first stored solution
# ------------------------------------------------------------
sols[[1]]
}
# -------------------------------------------------------------------------
# Internal: extract raw decision vector from a multiscape Solution
# -------------------------------------------------------------------------
.pamo_get_solution_vector <- function(sol) {
if (!is.null(sol$solution$vector) && is.numeric(sol$solution$vector) && length(sol$solution$vector) > 0) {
return(as.numeric(sol$solution$vector))
}
# fallback defensivo por si entra una lista rara interna
candidates <- c("vector", "solution", "sol", "x", "best_solution", "solution_vector")
d <- sol$solution %||% sol
for (nm in candidates) {
v <- d[[nm]] %||% NULL
if (is.numeric(v) && length(v) > 0) {
return(as.numeric(v))
}
}
stop(
"Could not extract the raw decision vector from the solution object.\n",
"Expected x$solution$vector or an equivalent numeric field.",
call. = FALSE
)
}
.pamo_eval_boundary_cut_on_solution <- function(x, solution, term) {
stopifnot(inherits(x, "Problem"))
rel_name <- as.character(term$relation_name %||% "boundary")[1]
rel <- x$data$spatial_relations[[rel_name]] %||% NULL
if (is.null(rel) || !inherits(rel, "data.frame")) {
stop("Missing spatial relation '", rel_name, "' while evaluating boundary_cut.", call. = FALSE)
}
for (nm in c("internal_pu1", "internal_pu2", "weight")) {
if (!nm %in% names(rel)) {
stop("Spatial relation '", rel_name, "' must contain column '", nm, "'.", call. = FALSE)
}
}
sol_vec <- .pamo_get_solution_vector(solution)
n_pu <- nrow(x$data$pu)
if (!is.finite(n_pu) || is.na(n_pu) || n_pu <= 0L) {
stop("Could not determine n_pu while evaluating boundary_cut.", call. = FALSE)
}
if (length(sol_vec) < n_pu) {
stop(
"Solution vector is shorter than n_pu while evaluating boundary_cut.\n",
"length(sol_vec) = ", length(sol_vec), ", n_pu = ", n_pu, ".",
call. = FALSE
)
}
# Phase 1/2A assumption: w occupies the first n_pu decision variables
w <- as.numeric(sol_vec[seq_len(n_pu)])
w_bin <- as.numeric(w > 0.5)
ip1 <- as.integer(rel$internal_pu1)
ip2 <- as.integer(rel$internal_pu2)
ww <- as.numeric(rel$weight)
ok <- is.finite(ip1) & is.finite(ip2) & is.finite(ww)
ip1 <- ip1[ok]
ip2 <- ip2[ok]
ww <- ww[ok]
if (length(ip1) == 0L) return(0)
directed <- if ("directed" %in% names(rel)) unique(as.logical(rel$directed)) else FALSE
if (length(directed) != 1L || is.na(directed)) stop("Invalid directed metadata in spatial relation.", call. = FALSE)
val <- 0
for (r in seq_along(ip1)) {
i <- ip1[r]; j <- ip2[r]; wij <- ww[r]
if (i == j) {
if (directed) stop("Directed spatial relations cannot contain self-arcs.", call. = FALSE)
val <- val + wij * w_bin[i]
} else if (directed) {
val <- val + wij * w_bin[i] * (1 - w_bin[j])
} else {
val <- val + wij * abs(w_bin[i] - w_bin[j])
}
}
mult <- as.numeric(term$weight_multiplier %||% 1)[1]
if (!is.finite(mult)) mult <- 1
as.numeric(mult * val)
}
# -------------------------------------------------------------------------
# Internal: evaluate one registered objective alias on a solved solution
# Uses the same IR/objective-vector machinery, so evaluation is consistent
# with the MO engine.
# -------------------------------------------------------------------------
.pamo_get_alias_value_from_solution <- function(x, solution, alias) {
stopifnot(inherits(x, "Problem"))
.pamo_eval_alias_on_solution(x, solution, alias)
}
.pamo_eval_alias_on_solution <- function(x, solution, alias) {
stopifnot(inherits(x, "Problem"))
alias <- as.character(alias)[1]
if (is.na(alias) || !nzchar(alias)) {
stop("alias must be a non-empty string.", call. = FALSE)
}
spec <- .pamo_get_objective_spec(x, alias)
ir <- .pamo_objective_to_ir(x, spec)
terms <- ir$terms %||% list()
# Fase 1:
# si el objetivo es exactamente min_fragmentation (PU boundary cut),
# lo evaluamos directamente desde w sin reconstruir auxiliares y_pu.
if (length(terms) == 1L && identical(terms[[1]]$type %||% "", "boundary_cut")) {
val <- .pamo_eval_boundary_cut_on_solution(x, solution, terms[[1]])
return(as.numeric(val))
}
if (length(terms) == 1L && identical(terms[[1]]$type %||% "", "intervention_impact")) {
val <- .pamo_eval_intervention_impact_on_solution(x, solution, terms[[1]])
return(as.numeric(val))
}
if (length(terms) == 1L && identical(terms[[1]]$type %||% "", "action_boundary_cut")) {
return(as.numeric(.pamo_eval_action_boundary_cut_on_solution(x, solution, terms[[1]])))
}
# Standard selection costs need no new MILP or auxiliary-variable snapshot.
# Use the solved model's exact variable layout and canonical cost tables.
cost_terms <- length(terms) > 0L && all(vapply(terms, function(t) {
t$type %in% c("pu_cost", "action_cost") && is.null(t$features)
}, logical(1)))
if (cost_terms) {
problem <- solution$problem %||% x
ml <- problem$data$model_list
da <- problem$data$dist_actions_model %||% problem$data$dist_actions
pu <- problem$data$pu
values <- .pamo_get_solution_vector(solution)
val <- 0
for (t in terms) {
if (identical(t$type, "pu_cost")) {
w <- values[seq_len(nrow(pu))]
selected <- which(w > 0.5)
val <- val + sum(pu$cost[selected])
} else {
ids <- suppressWarnings(as.integer(t$actions))
action_ids <- if (length(ids) && all(!is.na(ids)) && all(ids %in% problem$data$actions$internal_id)) {
ids
} else {
.pa_resolve_action_subset(problem, t$actions)$internal_id
}
rows <- which(da$internal_action %in% action_ids)
x_offset <- as.integer(ml$x_offset %||% nrow(pu))
internal_row <- da$internal_row %||% seq_len(nrow(da))
xv <- values[x_offset + internal_row[rows]]
selected <- which(xv > 0.5)
val <- val + sum(da$cost[rows[selected]])
}
}
return(as.numeric(val))
}
# resto de objetivos: evaluación por objvec
base_eval <- .pamo_prepare_superset_model(x, list(ir))
# Scoped loss auxiliaries need not occupy the same columns when rebuilt for
# evaluation. Evaluate ecological changes from the atomic action decisions,
# not from a prefix of the solution's auxiliary columns.
if (.pa_uses_aggregate_effects(base_eval) && length(terms) == 1L &&
terms[[1]]$type %in% c("benefit", "loss")) {
term <- terms[[1]]
return(.pa_eval_aggregate_effect_objective(base_eval,
.pamo_get_solution_vector(solution), term$type, term$actions, term$features))
}
obj_vec <- .pamo_objvec_from_ir(base_eval, ir)
sol_vec <- .pamo_get_solution_vector(solution)
n_obj <- length(obj_vec)
n_sol <- length(sol_vec)
if (n_obj > n_sol) {
stop(
"Objective vector length does not match solution vector length when evaluating alias '", alias, "'.\n",
"length(obj_vec) = ", n_obj, ", length(sol_vec) = ", n_sol, ".",
call. = FALSE
)
}
# permitir evaluar objetivos 'base' sobre soluciones de un superset
sol_use <- as.numeric(sol_vec[seq_len(n_obj)])
# Evaluate the implemented binary decisions, as the spatial summaries and
# legacy impact evaluator do. Solver integrality noise must not become an
# ecological amount or an automatic epsilon endpoint.
binary <- base_eval$data$model_list$vtype == "B"
sol_use[binary] <- as.numeric(sol_use[binary] > 0.5)
val_engine <- sum(as.numeric(obj_vec) * sol_use, na.rm = TRUE)
as.numeric(val_engine)
}
.pamo_add_alias_upper_bound_constraint <- function(base_eval, x, alias, rhs, tol = 0, name = NULL) {
stopifnot(inherits(base_eval, "Problem"))
stopifnot(inherits(x, "Problem"))
alias <- as.character(alias)[1]
rhs <- as.numeric(rhs)[1]
tol <- as.numeric(tol)[1]
if (!is.finite(rhs)) {
stop("rhs must be finite in .pamo_add_alias_upper_bound_constraint().", call. = FALSE)
}
if (!is.finite(tol) || tol < 0) {
stop("tol must be a finite non-negative number.", call. = FALSE)
}
spec <- .pamo_get_objective_spec(x, alias)
ir <- .pamo_objective_to_ir(x, spec)
obj_vec <- .pamo_objvec_from_ir(base_eval, ir)
op_list <- .pa_model_from_ptr(
base_eval$data$model_ptr,
args = base_eval$data$model_args %||% list(),
drop_triplets = TRUE
)
n_model <- as.integer(op_list$ncol %||% length(obj_vec))
if (length(obj_vec) != n_model) {
stop(
"Cannot add alias bound for '", alias, "': objective vector length does not match model dimension.\n",
"length(obj_vec) = ", length(obj_vec), ", model dimension = ", n_model, ".",
call. = FALSE
)
}
nz <- which(abs(obj_vec) > 0)
if (length(nz) == 0L) {
rcpp_add_linear_constraint(base_eval$data$model_ptr, integer(), numeric(), "<=",
as.numeric(rhs + tol), name = name %||% paste0("eps_bound_", alias))
return(.pa_refresh_model_snapshot(base_eval))
}
if (is.null(name)) {
name <- paste0("eps_bound_", alias)
}
base_eval <- .pa_add_linear_constraint(
base_eval,
var_index_0based = as.integer(nz - 1L),
coeff = as.numeric(obj_vec[nz]),
sense = "<=",
rhs = as.numeric(rhs + tol),
name = name
)
base_eval
}
# -------------------------------------------------------------------------
# Internal: compute 2-objective extreme points for auto epsilon mode
#
# Returns:
# - primary_at_primary
# - secondary_at_primary
# - primary_at_secondary
# - secondary_at_secondary
# - secondary_min
# - secondary_max
# -------------------------------------------------------------------------
.pamo_compute_epsilon_extremes_2obj <- function(x, primary, secondary, gap_limit = NULL, time_limit = NULL) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% list()
do_lexi <- isTRUE(method$lexicographic)
lexi_tol <- as.numeric(method$lexicographic_tol %||% 0)[1]
if (!is.finite(lexi_tol) || lexi_tol < 0) lexi_tol <- 0
specs_all <- .pamo_get_specs(x)
obj_alias <- names(specs_all)
if (!primary %in% obj_alias) {
stop("primary alias not found: '", primary, "'.", call. = FALSE)
}
if (!secondary %in% obj_alias) {
stop("secondary alias not found: '", secondary, "'.", call. = FALSE)
}
# ------------------------------------------------------------
# Caso simple (no lexicográfico)
# ------------------------------------------------------------
if (!isTRUE(do_lexi)) {
# Extreme 1: optimize primary
sol_primary <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = primary,
weights = 1,
objective_scaling = FALSE,
scales = NULL
)
)
primary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary$solution, primary)
secondary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary$solution, secondary)
# Extreme 2: optimize secondary
sol_secondary <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = secondary,
weights = 1,
normalize = FALSE,
gap_limit = gap_limit,
time_limit = time_limit
)
)
primary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary$solution, primary)
secondary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary$solution, secondary)
vals_secondary <- c(secondary_at_primary, secondary_at_secondary)
if (any(!is.finite(vals_secondary))) {
stop(
"Could not compute finite extreme values for the secondary objective.\n",
"secondary_at_primary = ", secondary_at_primary, "\n",
"secondary_at_secondary = ", secondary_at_secondary,
call. = FALSE
)
}
return(list(
primary = primary,
secondary = secondary,
primary_at_primary = as.numeric(primary_at_primary),
secondary_at_primary = as.numeric(secondary_at_primary),
primary_at_secondary = as.numeric(primary_at_secondary),
secondary_at_secondary = as.numeric(secondary_at_secondary),
secondary_min = min(vals_secondary),
secondary_max = max(vals_secondary),
vals_secondary = vals_secondary,
meta = list(
lexicographic = FALSE,
lexicographic_tol = 0
)
))
}
# ------------------------------------------------------------
# Caso lexicográfico
# ------------------------------------------------------------
# 1) resolver primary puro
sol_primary_stage1 <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = primary,
weights = 1,
normalize = FALSE,
gap_limit = gap_limit,
time_limit = time_limit
)
)
primary_star <- .pamo_eval_alias_on_solution(x, sol_primary_stage1$solution, primary)
# 2) resolver secondary sujeto a primary <= primary_star + tol
sol_primary_lexi <- .pamo_solve_one(
x = x,
spec = list(
type = "epsilon_constraint",
primary = secondary,
eps = stats::setNames(list(primary_star), primary),
gap_limit = gap_limit,
time_limit = time_limit,
eps_tol = stats::setNames(list(lexi_tol), primary)
)
)
primary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary_lexi$solution, primary)
secondary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary_lexi$solution, secondary)
# 3) resolver secondary puro
sol_secondary_stage1 <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = secondary,
weights = 1,
normalize = FALSE,
gap_limit = gap_limit,
time_limit = time_limit
)
)
secondary_star <- .pamo_eval_alias_on_solution(x, sol_secondary_stage1$solution, secondary)
# 4) resolver primary sujeto a secondary <= secondary_star + tol
sol_secondary_lexi <- .pamo_solve_one(
x = x,
spec = list(
type = "epsilon_constraint",
primary = primary,
eps = stats::setNames(list(secondary_star), secondary),
gap_limit = gap_limit,
time_limit = time_limit,
eps_tol = stats::setNames(list(lexi_tol), secondary)
)
)
primary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary_lexi$solution, primary)
secondary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary_lexi$solution, secondary)
vals_secondary <- c(secondary_at_primary, secondary_at_secondary)
if (any(!is.finite(vals_secondary))) {
stop(
"Could not compute finite extreme values for the secondary objective.\n",
"secondary_at_primary = ", secondary_at_primary, "\n",
"secondary_at_secondary = ", secondary_at_secondary,
call. = FALSE
)
}
list(
primary = primary,
secondary = secondary,
primary_at_primary = as.numeric(primary_at_primary),
secondary_at_primary = as.numeric(secondary_at_primary),
primary_at_secondary = as.numeric(primary_at_secondary),
secondary_at_secondary = as.numeric(secondary_at_secondary),
secondary_min = min(vals_secondary),
secondary_max = max(vals_secondary),
vals_secondary = vals_secondary,
meta = list(
lexicographic = TRUE,
lexicographic_tol = lexi_tol
)
)
}
.pamo_compile_superset_spec <- function(base, ir_list) {
stopifnot(inherits(base, "Problem"))
stopifnot(is.list(ir_list), length(ir_list) > 0)
all_terms <- unlist(
lapply(ir_list, function(ir) ir$terms %||% list()),
recursive = FALSE
)
all_types <- unique(vapply(all_terms, `[[`, character(1), "type"))
need_z <- "representation" %in% all_types
need_y_pu <- "boundary_cut" %in% all_types
need_y_act <- "action_boundary_cut" %in% all_types
need_y_int <- "intervention_boundary_cut" %in% all_types
need_u_int <- "intervention_impact" %in% all_types
if (isTRUE(need_y_int)) {
stop(
"Phase 1 superset does not yet support 'intervention_boundary_cut' objectives.\n",
"Please remove intervention fragmentation from the MO method for now.",
call. = FALSE
)
}
# ------------------------------------------------------------
# Spatial fragmentation terms: at most one relation_name in phase 1
# ------------------------------------------------------------
spatial_terms <- Filter(
function(t) t$type %in% c("boundary_cut", "action_boundary_cut", "intervention_boundary_cut"),
all_terms
)
rel_names <- unique(vapply(
spatial_terms,
function(t) as.character(t$relation_name %||% "boundary")[1],
character(1)
))
rel_names <- rel_names[!is.na(rel_names) & nzchar(rel_names)]
if (length(rel_names) > 1L) {
stop(
"Phase 1 superset supports only one relation_name across fragmentation objectives.\n",
"Found: ", paste(rel_names, collapse = ", "),
call. = FALSE
)
}
# ------------------------------------------------------------
# intervention_impact: phase 1 supported ONLY if all terms
# share the same action subset
# ------------------------------------------------------------
u_int_actions <- NULL
if (isTRUE(need_u_int)) {
int_terms <- Filter(function(t) identical(t$type, "intervention_impact"), all_terms)
norm_actions <- function(x) {
x <- as.character(x %||% character(0))
x <- unique(x[!is.na(x) & nzchar(x)])
sort(x)
}
action_keys <- vapply(
int_terms,
function(t) paste(norm_actions(t$actions), collapse = "||"),
character(1)
)
uniq_keys <- unique(action_keys)
if (length(uniq_keys) != 1L) {
pretty <- vapply(
int_terms,
function(t) {
aa <- norm_actions(t$actions)
if (length(aa) == 0L) "<empty>" else paste(aa, collapse = ", ")
},
character(1)
)
stop(
"Phase 1 superset supports 'intervention_impact' only when all such objectives ",
"share the same action subset.\n",
"Found subsets: ", paste(unique(pretty), collapse = " | "),
call. = FALSE
)
}
u_int_actions <- norm_actions(int_terms[[1]]$actions)
if (length(u_int_actions) == 0L) {
stop(
"Phase 1 superset requires a non-empty action subset for 'intervention_impact'.",
call. = FALSE
)
}
}
list(
needs = list(
z = isTRUE(need_z),
y_pu = isTRUE(need_y_pu),
y_action = isTRUE(need_y_act),
y_intervention = FALSE,
u_intervention = isTRUE(need_u_int),
u_intervention_actions = u_int_actions,
effect_loss_scopes = Filter(function(t) identical(t$type, "loss"), all_terms)
),
relation_name = if (length(rel_names) == 1L) rel_names[1] else NULL
)
}
.pamo_apply_single_objective <- function(base, ir, sense = c("min", "max")) {
stopifnot(inherits(base, "Problem"))
sense <- match.arg(sense)
if (is.null(base$data$model_ptr)) {
stop("Model not built (model_ptr is NULL).", call. = FALSE)
}
v <- .pamo_objvec_from_ir(base, ir)
# El motor trabaja siempre en minimización
if (identical(sense, "max")) {
v <- -as.numeric(v)
} else {
v <- as.numeric(v)
}
# Si existe un setter C++ del objetivo, úsalo
if (exists("rcpp_model_set_objective_vector", mode = "function")) {
rcpp_model_set_objective_vector(
x = base$data$model_ptr,
obj = v,
model_sense = "min"
)
} else {
# fallback temporal: deja la actualización pendiente
base$data$runtime_updates <- list(
obj = v,
modelsense = "min"
)
}
base$data$meta <- base$data$meta %||% list()
base$data$meta$model_dirty <- FALSE
base$data$has_model <- TRUE
base
}
.pamo_eval_action_boundary_cut_on_solution <- function(x, solution, term) {
stopifnot(inherits(x, "Problem"))
rel_name <- as.character(term$relation_name %||% "boundary")[1]
sol_problem <- solution$problem %||% NULL
rel <- sol_problem$data$spatial_relations_model[[rel_name]] %||%
sol_problem$data$spatial_relations[[rel_name]] %||%
x$data$spatial_relations_model[[rel_name]] %||%
x$data$spatial_relations[[rel_name]] %||%
NULL
if (is.null(rel) || !inherits(rel, "data.frame")) {
stop("Missing spatial relation '", rel_name, "' while evaluating action_boundary_cut.", call. = FALSE)
}
for (nm in c("internal_pu1", "internal_pu2", "weight")) {
if (!nm %in% names(rel)) {
stop("Spatial relation '", rel_name, "' must contain column '", nm, "'.", call. = FALSE)
}
}
da <- sol_problem$data$dist_actions_model %||%
sol_problem$data$dist_actions %||%
x$data$dist_actions_model %||%
x$data$dist_actions %||%
NULL
if (is.null(da) || !inherits(da, "data.frame")) {
stop("Missing dist_actions data while evaluating action_boundary_cut.", call. = FALSE)
}
for (nm in c("internal_pu", "internal_action")) {
if (!nm %in% names(da)) {
stop("dist_actions data must contain column '", nm, "'.", call. = FALSE)
}
}
sol_vec <- .pamo_get_solution_vector(solution)
ml <- sol_problem$data$model_list %||% x$data$model_list %||% NULL
n_pu <- as.integer(ml$n_pu %||% nrow(x$data$pu))
if (!is.finite(n_pu) || is.na(n_pu) || n_pu <= 0L) {
stop("Could not determine n_pu while evaluating action_boundary_cut.", call. = FALSE)
}
n_x <- as.integer(ml$n_x %||% nrow(da))
if (!is.finite(n_x) || is.na(n_x) || n_x <= 0L) {
stop("Could not determine n_x while evaluating action_boundary_cut.", call. = FALSE)
}
# IMPORTANT:
# fallback must follow the structural layout w | x | ...
x0 <- as.integer(ml$x_offset %||% n_pu)
if (!is.finite(x0) || is.na(x0) || x0 < 0L) {
stop("Could not determine x_offset while evaluating action_boundary_cut.", call. = FALSE)
}
if (!"internal_row" %in% names(da)) {
if (nrow(da) != n_x) {
stop(
"dist_actions data does not contain 'internal_row', and nrow(dist_actions) != n_x.\n",
"nrow(dist_actions) = ", nrow(da), ", n_x = ", n_x, ".",
call. = FALSE
)
}
da$internal_row <- seq_len(nrow(da))
}
if (any(!is.finite(da$internal_row)) || any(is.na(da$internal_row))) {
stop("dist_actions$internal_row contains NA/non-finite values.", call. = FALSE)
}
if (any(da$internal_row < 1L) || any(da$internal_row > n_x)) {
stop("dist_actions$internal_row is out of valid range 1..n_x.", call. = FALSE)
}
if (length(sol_vec) < (x0 + n_x)) {
stop(
"Solution vector is too short for action boundary evaluation.\n",
"length(sol_vec) = ", length(sol_vec),
", required >= ", x0 + n_x, ".",
call. = FALSE
)
}
all_actions <- sort(unique(as.integer(da$internal_action)))
all_actions <- all_actions[is.finite(all_actions) & !is.na(all_actions)]
act_raw <- term$actions %||% NULL
if (is.null(act_raw) || length(act_raw) == 0L) {
act_int <- all_actions
} else {
numeric_ids <- suppressWarnings(as.integer(act_raw))
act_int <- if (all(!is.na(numeric_ids)) && all(numeric_ids %in% all_actions)) {
numeric_ids
} else {
.pa_resolve_action_subset(sol_problem %||% x, act_raw)$internal_id
}
act_int <- sort(unique(act_int))
}
if (length(act_int) == 0L) return(0)
aw <- term$action_weights %||% NULL
if (is.null(aw)) {
aw_map <- stats::setNames(rep(1, length(act_int)), act_int)
} else {
aw <- as.numeric(aw)
if (length(aw) == length(all_actions)) {
tmp_map <- stats::setNames(aw, all_actions)
aw_map <- tmp_map[as.character(act_int)]
} else if (length(aw) == length(act_int)) {
aw_map <- stats::setNames(aw, act_int)
} else {
stop(
"action_weights length must match either all actions or the selected action subset.",
call. = FALSE
)
}
}
da_sub <- da[, c("internal_row", "internal_pu", "internal_action"), drop = FALSE]
da_sub <- da_sub[da_sub$internal_action %in% act_int, , drop = FALSE]
if (nrow(da_sub) == 0L) return(0)
x_vals <- sol_vec[x0 + da_sub$internal_row]
x_bin <- as.numeric(x_vals > 0.5)
ip1 <- as.integer(rel$internal_pu1)
ip2 <- as.integer(rel$internal_pu2)
ww <- as.numeric(rel$weight)
ok <- is.finite(ip1) & is.finite(ip2) & is.finite(ww)
ip1 <- ip1[ok]
ip2 <- ip2[ok]
ww <- ww[ok]
directed <- if ("directed" %in% names(rel)) unique(as.logical(rel$directed)) else FALSE
if (length(directed) != 1L || is.na(directed)) stop("Invalid directed metadata in spatial relation.", call. = FALSE)
val <- 0
for (act in act_int) {
awi <- as.numeric(aw_map[[as.character(act)]])
if (!is.finite(awi) || awi == 0) next
x_i <- numeric(n_pu)
idx_act <- da_sub$internal_action == act
if (any(idx_act)) x_i[da_sub$internal_pu[idx_act]] <- x_bin[idx_act]
for (r in seq_along(ip1)) {
i <- ip1[r]; j <- ip2[r]; wij <- ww[r]
if (i == j) {
if (directed) stop("Directed spatial relations cannot contain self-arcs.", call. = FALSE)
val <- val + awi * wij * x_i[i]
} else if (directed) {
val <- val + awi * wij * x_i[i] * (1 - x_i[j])
} else {
val <- val + awi * wij * abs(x_i[i] - x_i[j])
}
}
}
mult <- as.numeric(term$weight_multiplier %||% 1)[1]
if (!is.finite(mult)) mult <- 1
as.numeric(mult * val)
}
.pamo_eval_intervention_impact_on_solution <- function(x, solution, term) {
stopifnot(inherits(x, "Problem"))
icol <- as.character(term$impact_col %||% "amount")[1]
feats <- as.character(term$features %||% character(0))
acts <- as.character(term$actions %||% character(0))
distf <- x$data$dist_features %||% NULL
if (is.null(distf) || !inherits(distf, "data.frame") || nrow(distf) == 0) {
stop("dist_features is missing or empty while evaluating intervention_impact.", call. = FALSE)
}
if (!(icol %in% names(distf))) {
stop("impact column '", icol, "' not found in dist_features.", call. = FALSE)
}
# filtrar features si corresponde
if (length(feats) > 0) {
feat_subset <- .pa_resolve_feature_subset(x, features = feats)
feat_ids <- unique(as.character(feat_subset$id))
distf <- distf[as.character(distf$feature) %in% feat_ids, , drop = FALSE]
}
if (nrow(distf) == 0) return(0)
# resolver subset de acciones
act_subset <- .pa_resolve_action_subset(x, subset = acts)
if (!inherits(act_subset, "data.frame") || nrow(act_subset) == 0) {
stop("intervention_impact resolved to zero actions.", call. = FALSE)
}
act_ids <- unique(as.character(act_subset$id))
# tabla de acciones de la solución
act_tbl <- solution$summary$actions %||% NULL
if (is.null(act_tbl) || !inherits(act_tbl, "data.frame")) {
stop("No actions table found in solution while evaluating intervention_impact.", call. = FALSE)
}
if (!("selected" %in% names(act_tbl))) {
stop("Solution actions table has no 'selected' column.", call. = FALSE)
}
# columnas esperadas: pu, action
pu_col <- if ("pu" %in% names(act_tbl)) "pu" else if ("internal_pu" %in% names(act_tbl)) "internal_pu" else NULL
ac_col <- if ("action" %in% names(act_tbl)) "action" else if ("id" %in% names(act_tbl)) "id" else NULL
if (is.null(pu_col) || is.null(ac_col)) {
stop(
"Could not identify PU/action columns in solution$summary$actions while evaluating intervention_impact.",
call. = FALSE
)
}
act_sel <- act_tbl[
act_tbl$selected > 0.5 & as.character(act_tbl[[ac_col]]) %in% act_ids,
,
drop = FALSE
]
if (nrow(act_sel) == 0) return(0)
pu_hit <- unique(act_sel[[pu_col]])
# sumar impacto por PU intervenida
val <- sum(
distf[distf$pu %in% pu_hit, icol, drop = TRUE],
na.rm = TRUE
)
as.numeric(val)
}
.pamo_bind_solution_summaries <- function(solutions, run_ids = NULL) {
if (is.null(solutions) || length(solutions) == 0L) {
return(list(
pu = NULL,
actions = NULL,
features = NULL,
targets = NULL
))
}
if (is.null(run_ids)) {
run_ids <- seq_along(solutions)
}
if (length(run_ids) != length(solutions)) {
stop("run_ids length must match solutions length.", call. = FALSE)
}
bind_one <- function(name) {
out <- vector("list", length(solutions))
for (i in seq_along(solutions)) {
sol <- solutions[[i]]
if (!inherits(sol, "Solution")) next
tb <- sol$summary[[name]] %||% NULL
if (is.null(tb) || !inherits(tb, "data.frame") || nrow(tb) == 0L) next
tb$run_id <- run_ids[i]
out[[i]] <- tb
}
out <- Filter(Negate(is.null), out)
if (length(out) == 0L) return(NULL)
# usar tu helper de fill si las columnas no coinciden exactamente
acc <- out[[1]]
if (length(out) > 1L) {
for (j in 2:length(out)) {
acc <- .pa_rbind_fill(acc, out[[j]])
}
}
acc
}
list(
pu = bind_one("pu"),
actions = bind_one("actions"),
features = bind_one("features"),
targets = bind_one("targets")
)
}
.pamo_eval_alias_canonical_on_solution <- function(x, solution, alias) {
stopifnot(inherits(x, "Problem"))
alias <- as.character(alias)[1]
if (is.na(alias) || !nzchar(alias)) {
stop("alias must be a non-empty string.", call. = FALSE)
}
val <- .pamo_eval_alias_on_solution(x, solution, alias)
spec <- .pamo_get_objective_spec(x, alias)
sense <- as.character(spec$sense %||% NA_character_)[1]
if (!sense %in% c("min", "max")) {
stop("Objective alias '", alias, "' has invalid sense.", call. = FALSE)
}
if (identical(sense, "max")) {
return(as.numeric(-val))
}
as.numeric(val)
}
.pamo_compute_weighted_ranges <- function(x,
aliases,
tol = 1e-9,
verbose = FALSE,
gap_limit = NULL,
time_limit = NULL) {
stopifnot(inherits(x, "Problem"))
aliases <- as.character(aliases)
if (length(aliases) == 0L) {
stop("aliases must be non-empty in .pamo_compute_weighted_ranges().", call. = FALSE)
}
.pamo_cli_step(
"Computing payoff table for {.val {length(aliases)}} objective{?s}.",
verbose = verbose
)
payoff <- matrix(
NA_real_,
nrow = length(aliases),
ncol = length(aliases),
dimnames = list(aliases, aliases)
)
for (i in seq_along(aliases)) {
a_i <- aliases[i]
.pamo_cli_step(
"Anchor {.val {i}}/{.val {length(aliases)}}: optimizing {.val {a_i}}.",
verbose = verbose
)
one <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = a_i,
weights = 1,
objective_scaling = FALSE,
scales = NULL,
gap_limit = gap_limit,
time_limit = time_limit
)
)
sol_i <- one$solution
if (is.null(sol_i) || !inherits(sol_i, "Solution")) {
stop(
"Failed to compute scaling range anchor for objective '", a_i, "'.",
call. = FALSE
)
}
vals_i <- vapply(
aliases,
function(a_j) .pamo_eval_alias_canonical_on_solution(x, sol_i, a_j),
numeric(1)
)
payoff[i, ] <- vals_i
}
ideal <- apply(payoff, 2, min, na.rm = TRUE)
nadir <- apply(payoff, 2, max, na.rm = TRUE)
ranges <- nadir - ideal
names(ranges) <- aliases
bad <- !is.finite(ranges) | (ranges <= tol)
if (any(bad)) {
.pamo_cli_warn(
"Some objectives had near-zero payoff range: {.val {paste(names(ranges)[bad], collapse = ', ')}}. Using scale = 1.",
verbose = verbose
)
ranges[bad] <- 1
}
.pamo_cli_done("Payoff ranges computed.", verbose = verbose)
ranges
}
################################################################################
######## AUGMECON ##############################################################
################################################################################
.pamo_apply_augmecon_objective <- function(base,
ir_primary,
slack_cols,
slack_weights) {
stopifnot(inherits(base, "Problem"))
slack_cols <- as.integer(slack_cols)
slack_weights <- as.numeric(slack_weights)
if (length(slack_cols) == 0L) {
stop("slack_cols must contain at least one slack column.", call. = FALSE)
}
if (length(slack_cols) != length(slack_weights)) {
stop(
"slack_cols and slack_weights must have the same length.",
call. = FALSE
)
}
if (any(!is.finite(slack_cols)) || any(is.na(slack_cols)) || any(slack_cols < 0L)) {
stop("slack_cols must contain valid 0-based column indices.", call. = FALSE)
}
if (any(!is.finite(slack_weights)) || any(is.na(slack_weights)) || any(slack_weights <= 0)) {
stop("slack_weights must contain finite positive values.", call. = FALSE)
}
v_primary <- .pamo_objvec_from_ir(base, ir_primary)
# canonical minimization form
if (identical(ir_primary$sense, "max")) {
v_primary <- -v_primary
}
obj <- as.numeric(v_primary)
# add AUGMECON slack contribution:
# minimize primary - sum_j weight_j * slack_j
for (i in seq_along(slack_cols)) {
j <- slack_cols[i] + 1L
if (j > length(obj)) {
stop(
"Slack column ", slack_cols[i], " is outside the model dimension (length(obj) = ",
length(obj), ").",
call. = FALSE
)
}
obj[j] <- obj[j] - slack_weights[i]
}
base$data$runtime_updates <- list(
obj = obj,
modelsense = "min"
)
base$data$meta <- base$data$meta %||% list()
base$data$meta$model_dirty <- FALSE
base$data$has_model <- TRUE
base$data$mo_cache <- base$data$mo_cache %||% list()
base$data$mo_cache$augmecon <- utils::modifyList(
base$data$mo_cache$augmecon %||% list(),
list(
ir_primary = ir_primary,
slack_cols_0based = slack_cols,
slack_weights = slack_weights,
obj_augmented = obj
)
)
base
}
.pamo_solve_augmecon <- function(x, ...) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% NULL
if (is.null(method) || !is.list(method)) {
stop("augmecon: x$data$method is missing or invalid.", call. = FALSE)
}
primary <- as.character(method$primary %||% NA_character_)[1]
aliases <- as.character(method$aliases %||% character(0))
secondary <- as.character(method$secondary %||% character(0))
# grid <- method$grid %||% NULL
# n_points <- as.integer(method$n_points %||% NA_integer_)[1]
augmentation <- as.numeric(method$augmentation %||% NA_real_)[1]
dots <- list(...)
gap_limit <- dots$gap_limit %||% NULL
time_limit <- dots$time_limit %||% NULL
verbose <- .pamo_cli_enabled(x, dots)
progress <- dots$progress %||% dots$verbose %||% interactive()
progress <- isTRUE(progress)
.pamo_cli_header("AUGMECON method", verbose = verbose)
if (is.na(primary) || !nzchar(primary)) {
stop("augmecon: missing primary objective.", call. = FALSE)
}
if (length(aliases) == 0L) {
stop("augmecon: aliases are missing.", call. = FALSE)
}
if (length(secondary) == 0L) {
stop("augmecon: secondary objectives are missing.", call. = FALSE)
}
if (!is.finite(augmentation) || is.na(augmentation) || augmentation <= 0) {
stop("augmecon: augmentation must be a finite positive number.", call. = FALSE)
}
secondary_ranges <- vapply(
secondary,
function(a) {
rr <- .pamo_compute_secondary_range_against_primary(
x = x,
primary = primary,
secondary = a,
gap_limit = gap_limit,
time_limit = time_limit
)
rng <- as.numeric(rr$secondary_max)[1] - as.numeric(rr$secondary_min)[1]
if (!is.finite(rng) || rng <= 0) rng <- 1
rng
},
numeric(1)
)
.pamo_cli_step(
"Primary objective: {.val {primary}}. Secondary objectives: {.val {paste(secondary, collapse = ', ')}}.",
verbose = verbose
)
# design
# design_df <- method$runs %||% NULL
# if (is.null(design_df)) {
# if (is.null(grid)) {
# .pamo_cli_step("Building automatic epsilon grid.", verbose = verbose)
#
# design_df <- .pamo_build_auto_augmecon_runs(
# x = x,
# gap_limit = gap_limit,
# time_limit = time_limit,
# verbose = verbose
# )
# } else {
# .pamo_cli_step("Using user-supplied epsilon grid.", verbose = verbose)
#
# design_df <- .pamo_build_manual_augmecon_runs(
# x = x,
# verbose = verbose
# )
# }
# x$data$method$runs <- design_df
# }
design_df <- .pamo_compile_augmecon_runs(
x = x,
method = method,
gap_limit = gap_limit,
time_limit = time_limit,
verbose = verbose
)
if (!("run_id" %in% names(design_df))) {
design_df$run_id <- seq_len(nrow(design_df))
}
eps_cols <- paste0("eps_", secondary)
miss_cols <- setdiff(eps_cols, names(design_df))
if (length(miss_cols) > 0L) {
stop(
"augmecon: design is missing epsilon columns: ",
paste(miss_cols, collapse = ", "),
call. = FALSE
)
}
.pamo_cli_step(
"Design ready with {.val {nrow(design_df)}} run{?s}. Augmentation = {.val {augmentation}}.",
verbose = verbose
)
n_runs <- nrow(design_df)
solutions <- vector("list", n_runs)
value_mat <- matrix(NA_real_, n_runs, length(aliases))
colnames(value_mat) <- paste0("value_", aliases)
status <- character(n_runs)
runtime <- numeric(n_runs)
gap <- numeric(n_runs)
messages <- character(n_runs)
stop_on_infeasible <- isTRUE(method$stop_on_infeasible %||% FALSE)
stop_on_no_solution <- isTRUE(method$stop_on_no_solution %||% FALSE)
stop_on_error <- isTRUE(method$stop_on_error %||% TRUE)
progress_id <- NULL
if (progress && n_runs > 1L) {
progress_id <- NULL
if (progress && n_runs > 1L) {
progress_id <- cli::cli_progress_bar(
name = "Solving runs",
total = n_runs,
clear = FALSE,
format = "{cli::pb_name}{cli::pb_bar} {cli::pb_current}/{cli::pb_total} | ETA: {cli::pb_eta}"
)
}
}
for (r in seq_len(n_runs)) {
eps_r <- as.list(design_df[r, eps_cols, drop = FALSE])
names(eps_r) <- secondary
one <- .pamo_solve_one_or_keep_failed(
x = x,
spec = list(
type = "augmecon",
primary = primary,
eps = eps_r,
eps_tol = stats::setNames(as.list(rep(0, length(secondary))), secondary),
augmentation = augmentation,
secondary_ranges = secondary_ranges,
gap_limit = gap_limit,
time_limit = time_limit
),
stop_on_infeasible = stop_on_infeasible,
stop_on_no_solution = stop_on_no_solution,
stop_on_error = stop_on_error
)
solutions[r] <- list(one$solution)
status[r] <- as.character(one$status %||% NA_character_)
runtime[r] <- as.numeric(one$runtime %||% NA_real_)
gap[r] <- as.numeric(one$gap %||% NA_real_)
messages[r] <- as.character(one$message %||% NA_character_)
if (!is.null(solutions[[r]]) && inherits(solutions[[r]], "Solution")) {
alias_values <- stats::setNames(
vapply(
aliases,
function(a) .pamo_eval_alias_on_solution(x, one$solution, a),
numeric(1)
),
aliases
)
solutions[[r]]$solution$alias_values <- alias_values
solutions[[r]]$meta$run_id <- design_df$run_id[r]
solutions[[r]]$method$type <- "augmecon"
solutions[[r]]$method$primary_alias <- primary
solutions[[r]]$method$eps <- stats::setNames(
as.numeric(design_df[r, eps_cols, drop = TRUE]),
secondary
)
solutions[[r]]$method$augmentation <- augmentation
value_mat[r, ] <- unname(as.numeric(alias_values))
} else {
value_mat[r, ] <- NA_real_
}
if (!is.null(progress_id)) {
cli::cli_progress_update(id = progress_id, set = r)
}
}
if (!is.null(progress_id)) {
cli::cli_progress_done(id = progress_id)
}
runs <- data.frame(
run_id = design_df$run_id,
status = status,
runtime = runtime,
gap = gap,
message = messages,
stringsAsFactors = FALSE
)
runs <- cbind(runs, as.data.frame(value_mat, stringsAsFactors = FALSE))
summary_set <- .pamo_bind_solution_summaries(
solutions = solutions,
run_ids = design_df$run_id
)
solutions <- .pa_strip_nested_solution_problems(solutions)
.pamo_cli_done("AUGMECON solve finished.", verbose = verbose)
ext <- attr(design_df, "extremes", exact = TRUE)
ranges_list <- list()
# rango del primario a partir de los anchors de todos los secundarios
primary_vals <- numeric(0)
if (!is.null(ext) && length(ext) > 0L) {
primary_vals <- unlist(
lapply(ext, function(e) {
c(
as.numeric(e$primary_at_primary %||% NA_real_),
as.numeric(e$primary_at_secondary %||% NA_real_)
)
}),
use.names = FALSE
)
primary_vals <- primary_vals[is.finite(primary_vals)]
}
ranges_list[[1]] <- data.frame(
alias = primary,
role = "primary",
min = if (length(primary_vals)) min(primary_vals) else NA_real_,
max = if (length(primary_vals)) max(primary_vals) else NA_real_,
source = "anchors",
stringsAsFactors = FALSE
)
if (!is.null(ext) && length(ext) > 0L) {
for (s in secondary) {
e <- ext[[s]]
if (is.null(e)) next
ranges_list[[length(ranges_list) + 1L]] <- data.frame(
alias = s,
role = "secondary",
min = as.numeric(e$secondary_min %||% NA_real_)[1],
max = as.numeric(e$secondary_max %||% NA_real_)[1],
source = "auto_grid",
stringsAsFactors = FALSE
)
}
}
ranges_df <- do.call(rbind, ranges_list)
rownames(ranges_df) <- NULL
n_infeasible <- sum(
runs$status %in% c("infeasible", "infeasible_or_unbounded"),
na.rm = TRUE
)
n_no_solution <- sum(
runs$status %in% c("no_solution"),
na.rm = TRUE
)
n_failed <- sum(
runs$status %in% c("failed"),
na.rm = TRUE
)
if (n_infeasible > 0L) {
warning(
n_infeasible, " of ", nrow(runs),
" AUGMECON runs were infeasible and were kept in the SolutionSet with missing objective values. ",
"Feasible runs were retained.",
call. = FALSE,
immediate. = TRUE
)
}
if (n_no_solution > 0L) {
warning(
n_no_solution, " of ", nrow(runs),
" AUGMECON runs did not return a solution vector and were kept in the SolutionSet with missing objective values. ",
"This may indicate infeasibility, numerical issues, or solver termination before finding a feasible solution.",
call. = FALSE,
immediate. = TRUE
)
}
if (n_failed > 0L) {
warning(
n_failed, " of ", nrow(runs),
" AUGMECON runs failed for reasons other than infeasibility/no-solution.",
call. = FALSE,
immediate. = TRUE
)
}
out <- pproto(
NULL, SolutionSet,
problem = x,
solution = list(
design = design_df,
runs = runs,
solutions = solutions,
frontier = list(
primary = primary,
secondary = secondary,
ranges = ranges_df,
extremes = ext
)
),
summary = summary_set,
diagnostics = list(
n_design = nrow(design_df),
n_runs = nrow(runs),
n_solutions = sum(vapply(solutions, inherits, logical(1), "Solution")),
n_optimal = sum(runs$status == "optimal", na.rm = TRUE),
n_infeasible = sum(runs$status %in% c("infeasible", "infeasible_or_unbounded"), na.rm = TRUE),
n_no_solution = sum(runs$status == "no_solution", na.rm = TRUE),
n_failed = sum(runs$status == "failed", na.rm = TRUE),
status_summary = .pa_solutionset_status_summary(runs),
runtime_range = if ("runtime" %in% names(runs)) .pa_solutionset_range_text(runs$runtime, digits = 3) else "none",
gap_range = if ("gap" %in% names(runs)) .pa_solutionset_range_text(runs$gap, digits = 4) else "none"
),
method = method,
meta = list(
call = match.call()
),
name = "solset"
)
.pa_finalize_solution_ids(out)
}
.pamo_compute_secondary_range_against_primary <- function(
x,
primary,
secondary,
gap_limit = NULL,
time_limit = NULL
) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% list()
do_lexi <- isTRUE(method$lexicographic)
lexi_tol <- as.numeric(method$lexicographic_tol %||% 0)[1]
if (!is.finite(lexi_tol) || lexi_tol < 0) lexi_tol <- 0
primary <- as.character(primary)[1]
secondary <- as.character(secondary)[1]
if (is.na(primary) || !nzchar(primary)) {
stop("`primary` must be a non-empty string.", call. = FALSE)
}
if (is.na(secondary) || !nzchar(secondary)) {
stop("`secondary` must be a non-empty string.", call. = FALSE)
}
if (identical(primary, secondary)) {
stop("`primary` and `secondary` must be different aliases.", call. = FALSE)
}
specs_all <- .pamo_get_specs(x)
obj_alias <- names(specs_all)
if (!primary %in% obj_alias) {
stop("Primary alias not found: '", primary, "'.", call. = FALSE)
}
if (!secondary %in% obj_alias) {
stop("Secondary alias not found: '", secondary, "'.", call. = FALSE)
}
# ------------------------------------------------------------
# Case 1: non-lexicographic anchors
# ------------------------------------------------------------
if (!isTRUE(do_lexi)) {
# Anchor A: optimize primary
sol_primary <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = primary,
weights = 1,
objective_scaling = FALSE,
scales = NULL
)
)
if (is.null(sol_primary$solution) || !inherits(sol_primary$solution, "Solution")) {
stop("Failed to compute primary anchor solution for '", primary, "'.", call. = FALSE)
}
primary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary$solution, primary)
secondary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary$solution, secondary)
# Anchor B: optimize secondary
sol_secondary <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = secondary,
weights = 1,
objective_scaling = FALSE,
scales = NULL
)
)
if (is.null(sol_secondary$solution) || !inherits(sol_secondary$solution, "Solution")) {
stop("Failed to compute secondary anchor solution for '", secondary, "'.", call. = FALSE)
}
primary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary$solution, primary)
secondary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary$solution, secondary)
vals_secondary <- c(secondary_at_primary, secondary_at_secondary)
if (any(!is.finite(vals_secondary))) {
stop(
"Could not compute finite anchor values for secondary objective '", secondary, "'.\n",
"secondary_at_primary = ", secondary_at_primary, "\n",
"secondary_at_secondary = ", secondary_at_secondary,
call. = FALSE
)
}
return(list(
primary = primary,
secondary = secondary,
primary_at_primary = as.numeric(primary_at_primary),
secondary_at_primary = as.numeric(secondary_at_primary),
primary_at_secondary = as.numeric(primary_at_secondary),
secondary_at_secondary = as.numeric(secondary_at_secondary),
secondary_min = min(vals_secondary),
secondary_max = max(vals_secondary),
# solution_primary = sol_primary$solution,
# solution_secondary = sol_secondary$solution,
meta = list(
lexicographic = FALSE,
lexicographic_tol = 0
)
))
}
# ------------------------------------------------------------
# Case 2: lexicographic anchors
# ------------------------------------------------------------
# Stage 1A: optimize primary
sol_primary_stage1 <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = primary,
weights = 1,
objective_scaling = FALSE,
scales = NULL
)
)
if (is.null(sol_primary_stage1$solution) || !inherits(sol_primary_stage1$solution, "Solution")) {
stop("Failed to compute stage-1 primary solution for '", primary, "'.", call. = FALSE)
}
primary_star <- .pamo_eval_alias_on_solution(x, sol_primary_stage1$solution, primary)
# Stage 2A: optimize secondary subject to primary <= primary_star + tol
sol_primary_lexi <- .pamo_solve_one(
x = x,
spec = list(
type = "epsilon_constraint",
primary = secondary,
eps = stats::setNames(list(primary_star), primary),
eps_tol = stats::setNames(list(lexi_tol), primary)
)
)
if (is.null(sol_primary_lexi$solution) || !inherits(sol_primary_lexi$solution, "Solution")) {
stop("Failed to compute lexicographic primary anchor for '", primary, "'.", call. = FALSE)
}
primary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary_lexi$solution, primary)
secondary_at_primary <- .pamo_eval_alias_on_solution(x, sol_primary_lexi$solution, secondary)
# Stage 1B: optimize secondary
sol_secondary_stage1 <- .pamo_solve_one(
x = x,
spec = list(
type = "weighted",
aliases = secondary,
weights = 1,
objective_scaling = FALSE,
scales = NULL
)
)
if (is.null(sol_secondary_stage1$solution) || !inherits(sol_secondary_stage1$solution, "Solution")) {
stop("Failed to compute stage-1 secondary solution for '", secondary, "'.", call. = FALSE)
}
secondary_star <- .pamo_eval_alias_on_solution(x, sol_secondary_stage1$solution, secondary)
# Stage 2B: optimize primary subject to secondary <= secondary_star + tol
sol_secondary_lexi <- .pamo_solve_one(
x = x,
spec = list(
type = "epsilon_constraint",
primary = primary,
eps = stats::setNames(list(secondary_star), secondary),
eps_tol = stats::setNames(list(lexi_tol), secondary)
)
)
if (is.null(sol_secondary_lexi$solution) || !inherits(sol_secondary_lexi$solution, "Solution")) {
stop("Failed to compute lexicographic secondary anchor for '", secondary, "'.", call. = FALSE)
}
primary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary_lexi$solution, primary)
secondary_at_secondary <- .pamo_eval_alias_on_solution(x, sol_secondary_lexi$solution, secondary)
vals_secondary <- c(secondary_at_primary, secondary_at_secondary)
if (any(!is.finite(vals_secondary))) {
stop(
"Could not compute finite anchor values for secondary objective '", secondary, "'.\n",
"secondary_at_primary = ", secondary_at_primary, "\n",
"secondary_at_secondary = ", secondary_at_secondary,
call. = FALSE
)
}
list(
primary = primary,
secondary = secondary,
primary_at_primary = as.numeric(primary_at_primary),
secondary_at_primary = as.numeric(secondary_at_primary),
primary_at_secondary = as.numeric(primary_at_secondary),
secondary_at_secondary = as.numeric(secondary_at_secondary),
secondary_min = min(vals_secondary),
secondary_max = max(vals_secondary),
# solution_primary = sol_primary_lexi$solution,
# solution_secondary = sol_secondary_lexi$solution,
meta = list(
lexicographic = TRUE,
lexicographic_tol = lexi_tol
)
)
}
.pamo_build_auto_augmecon_runs <- function(x,
gap_limit = NULL,
time_limit = NULL,
verbose = FALSE) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% list()
primary <- as.character(method$primary %||% NA_character_)[1]
aliases <- as.character(method$aliases %||% character(0))
secondary <- as.character(method$secondary %||% character(0))
n_points <- as.integer(method$n_points %||% NA_integer_)[1]
include_extremes <- isTRUE(method$include_extremes)
if (is.na(primary) || !nzchar(primary)) {
stop("Auto AUGMECON mode: missing primary objective.", call. = FALSE)
}
if (length(aliases) == 0L) {
stop("Auto AUGMECON mode: aliases are missing.", call. = FALSE)
}
if (!primary %in% aliases) {
stop("Auto AUGMECON mode: primary must be included in aliases.", call. = FALSE)
}
if (length(secondary) == 0L) {
stop("Auto AUGMECON mode requires at least one secondary objective.", call. = FALSE)
}
if (!all(secondary %in% aliases)) {
stop(
"Auto AUGMECON mode: all secondary objectives must be included in aliases.\n",
"Missing: ",
paste(setdiff(secondary, aliases), collapse = ", "),
call. = FALSE
)
}
if (!is.finite(n_points) || is.na(n_points) || n_points < 2L) {
stop("Auto AUGMECON mode requires n_points >= 2.", call. = FALSE)
}
.pamo_cli_step(
"Computing automatic ranges for {.val {length(secondary)}} secondary objective{?s}.",
verbose = verbose
)
# ------------------------------------------------------------
# build one epsilon sequence per secondary objective
# ------------------------------------------------------------
grid_list <- vector("list", length(secondary))
names(grid_list) <- secondary
extremes <- vector("list", length(secondary))
names(extremes) <- secondary
for (i in seq_along(secondary)) {
s <- secondary[i]
.pamo_cli_step(
"Secondary {.val {i}}/{.val {length(secondary)}}: {.val {s}}.",
verbose = verbose
)
ext_s <- .pamo_compute_secondary_range_against_primary(
x = x,
primary = primary,
secondary = s,
gap_limit = gap_limit,
time_limit = time_limit
)
s_min <- as.numeric(ext_s$secondary_min)[1]
s_max <- as.numeric(ext_s$secondary_max)[1]
if (!is.finite(s_min) || !is.finite(s_max)) {
stop(
"Auto AUGMECON mode: computed non-finite bounds for secondary objective '", s, "'.",
call. = FALSE
)
}
if (s_min > s_max) {
tmp <- s_min
s_min <- s_max
s_max <- tmp
}
.pamo_cli_step(
"Range for {.val {s}}: [{.val {format(s_min, digits = 6)}}, {.val {format(s_max, digits = 6)}}].",
verbose = verbose
)
vals <- seq(from = s_min, to = s_max, length.out = n_points)
if (!isTRUE(include_extremes)) {
if (length(vals) <= 2L) {
stop(
"include_extremes = FALSE requires n_points >= 3.",
call. = FALSE
)
}
vals <- vals[2:(length(vals) - 1L)]
}
if (length(vals) == 0L) {
stop(
"Auto AUGMECON mode produced an empty grid for secondary objective '", s, "'.",
call. = FALSE
)
}
grid_list[[s]] <- as.numeric(vals)
extremes[[s]] <- ext_s
}
# ------------------------------------------------------------
# Cartesian product of epsilon levels
# ------------------------------------------------------------
design_df <- expand.grid(
grid_list,
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
names(design_df) <- paste0("eps_", names(grid_list))
design_df$run_id <- seq_len(nrow(design_df))
# reorder columns: run_id first
design_df <- design_df[, c("run_id", setdiff(names(design_df), "run_id")), drop = FALSE]
attr(design_df, "extremes") <- extremes
attr(design_df, "primary_alias") <- primary
attr(design_df, "secondary_aliases") <- secondary
attr(design_df, "include_extremes") <- include_extremes
attr(design_df, "n_points_requested") <- n_points
.pamo_cli_done(
"Automatic AUGMECON grid built with {.val {nrow(design_df)}} run{?s}.",
verbose = verbose
)
design_df
}
# .pamo_build_manual_augmecon_runs <- function(x, verbose = FALSE) {
# stopifnot(inherits(x, "Problem"))
#
# method <- x$data$method %||% list()
# grid <- method$grid %||% NULL
# secondary <- as.character(method$secondary %||% character(0))
#
# if (is.null(grid) || !is.list(grid) || length(grid) == 0L) {
# stop("Manual AUGMECON mode requires a non-empty named grid list.", call. = FALSE)
# }
#
# if (length(secondary) == 0L) {
# stop("Manual AUGMECON mode: secondary objectives are missing.", call. = FALSE)
# }
#
# if (!all(secondary %in% names(grid))) {
# stop(
# "Manual AUGMECON mode: grid is missing secondary objectives: ",
# paste(setdiff(secondary, names(grid)), collapse = ", "),
# call. = FALSE
# )
# }
#
# .pamo_cli_step(
# "Building manual AUGMECON grid for {.val {length(secondary)}} secondary objective{?s}.",
# verbose = verbose
# )
#
# grid <- grid[secondary]
#
# for (i in seq_along(secondary)) {
# s <- secondary[i]
# vals <- grid[[s]]
#
# .pamo_cli_step(
# "Secondary {.val {i}}/{.val {length(secondary)}}: {.val {s}} with {.val {length(vals)}} epsilon level{?s}.",
# verbose = verbose
# )
# }
#
# design_df <- expand.grid(
# grid,
# KEEP.OUT.ATTRS = FALSE,
# stringsAsFactors = FALSE
# )
#
# names(design_df) <- paste0("eps_", secondary)
# design_df$run_id <- seq_len(nrow(design_df))
# design_df <- design_df[, c("run_id", setdiff(names(design_df), "run_id")), drop = FALSE]
#
# attr(design_df, "primary_alias") <- as.character(method$primary %||% NA_character_)[1]
# attr(design_df, "secondary_aliases") <- secondary
# attr(design_df, "manual_grid") <- TRUE
#
# .pamo_cli_done(
# "Manual AUGMECON grid built with {.val {nrow(design_df)}} run{?s}.",
# verbose = verbose
# )
#
# design_df
# }
.pamo_cli_enabled <- function(x, dots = list()) {
# prioridad:
# 1) verbose pasado a solve(...)
# 2) solver verbose si existe
# 3) interactive()
v <- dots$verbose %||%
x$data$solver$verbose %||%
x$data$solver_args$verbose %||%
interactive()
isTRUE(v)
}
.pamo_cli_header <- function(..., verbose = TRUE) {
if (isTRUE(verbose)) {
cli::cli_h2(..., .envir = parent.frame())
}
}
.pamo_cli_step <- function(..., verbose = TRUE) {
if (isTRUE(verbose)) {
cli::cli_alert_info(..., .envir = parent.frame())
}
}
.pamo_cli_done <- function(..., verbose = TRUE) {
if (isTRUE(verbose)) {
cli::cli_alert_success(..., .envir = parent.frame())
}
}
.pamo_cli_warn <- function(..., verbose = TRUE) {
if (isTRUE(verbose)) {
cli::cli_alert_warning(..., .envir = parent.frame())
}
}
.pamo_add_augmecon_slacks <- function(base,
secondary_aliases,
ub = 1e6,
prefix = "aug_slack") {
stopifnot(inherits(base, "Problem"))
secondary_aliases <- as.character(secondary_aliases)
secondary_aliases <- secondary_aliases[!is.na(secondary_aliases) & nzchar(secondary_aliases)]
if (length(secondary_aliases) == 0L) {
stop("secondary_aliases must contain at least one objective alias.", call. = FALSE)
}
ub <- as.numeric(ub)[1]
if (!is.finite(ub) || ub <= 0) {
stop("`ub` must be a single finite positive number.", call. = FALSE)
}
if (is.null(base$data$model_ptr)) {
stop("Model not built (model_ptr is NULL).", call. = FALSE)
}
# snapshot actual del modelo
m0 <- .pa_model_from_ptr(
base$data$model_ptr,
args = base$data$model_args %||% list(),
drop_triplets = TRUE
)
n_old <- length(m0$obj)
if (!is.finite(n_old) || n_old <= 0L) {
stop("Could not determine current number of model columns.", call. = FALSE)
}
k <- length(secondary_aliases)
# necesitamos una helper de bajo nivel para añadir columnas al modelo ya construido
if (!exists("rcpp_model_add_columns", mode = "function")) {
stop(
"AUGMECON with explicit slacks requires a low-level helper `rcpp_model_add_columns()`.\n",
"This helper is not available yet, so slack variables cannot be appended to the built model.",
call. = FALSE
)
}
slack_names <- paste0(prefix, "_", secondary_aliases)
# Añadir k columnas continuas:
# obj = 0 por ahora (el objetivo aumentado se define después)
# lb = 0
# ub = ub
# vtype = "C"
#
# Se asume que rcpp_model_add_columns():
# - modifica el modelo en sitio
# - acepta vectores de obj/lb/ub/vtype/names
# - deja el modelo consistente para luego añadir restricciones
rcpp_model_add_columns(
x = base$data$model_ptr,
obj = rep(0, k),
lb = rep(0, k),
ub = rep(ub, k),
vtype = rep("C", k),
names = slack_names
)
# refrescar snapshot/model_list tras añadir columnas
base <- .pa_refresh_model_snapshot(base)
m1 <- .pa_model_from_ptr(
base$data$model_ptr,
args = base$data$model_args %||% list(),
drop_triplets = TRUE
)
n_new <- length(m1$obj)
if (!identical(n_new, n_old + k)) {
stop(
"Failed to append AUGMECON slack variables.\n",
"Expected ", n_old + k, " columns after insertion, got ", n_new, ".",
call. = FALSE
)
}
# columnas nuevas en Ãndice 0-based
slack_cols_0based <- seq.int(from = n_old, length.out = k)
names(slack_cols_0based) <- secondary_aliases
# guardar metadata útil
base$data$mo_cache <- base$data$mo_cache %||% list()
base$data$mo_cache$augmecon <- utils::modifyList(
base$data$mo_cache$augmecon %||% list(),
list(
secondary_aliases = secondary_aliases,
slack_names = slack_names,
slack_cols_0based = slack_cols_0based,
slack_cols_1based = slack_cols_0based + 1L,
slack_upper_bound = ub
)
)
base$data$meta <- base$data$meta %||% list()
base$data$meta$model_dirty <- FALSE
base$data$has_model <- TRUE
base
}
.pamo_apply_augmecon_constraint <- function(base, ir, eps, slack_col, sense = c("min", "max"),
name = "", block_name = "augmecon_constraint", tag = "") {
stopifnot(inherits(base, "Problem"))
sense <- match.arg(sense)
if (is.null(base$data$model_ptr)) {
stop("Model not built (model_ptr is NULL).", call. = FALSE)
}
eps <- as.numeric(eps)[1]
if (!is.finite(eps)) {
stop("eps must be finite.", call. = FALSE)
}
slack_col <- as.integer(slack_col)[1]
if (!is.finite(slack_col) || is.na(slack_col) || slack_col < 0L) {
stop("slack_col must be a valid 0-based column index.", call. = FALSE)
}
v <- .pamo_objvec_from_ir(base, ir)
if (length(v) == 0L) {
stop("AUGMECON objective vector is empty.", call. = FALSE)
}
# canonical minimization form
if (identical(sense, "max")) {
v <- -v
eps <- -eps
}
# build equality: g(x) + s = eps
v2 <- as.numeric(v)
m <- .pa_model_from_ptr(
base$data$model_ptr,
args = base$data$model_args %||% list(),
drop_triplets = TRUE
)
n_model <- length(m$obj)
if (length(v2) != n_model) {
stop(
"AUGMECON constraint vector length does not match model dimension.\n",
"length(v2) = ", length(v2), ", model dimension = ", n_model, ".",
call. = FALSE
)
}
if (slack_col + 1L > n_model) {
stop("slack_col is outside model dimension.", call. = FALSE)
}
v2[slack_col + 1L] <- v2[slack_col + 1L] + 1
idx <- which(abs(v2) > 1e-12)
if (length(idx) == 0L) {
stop("AUGMECON constraint vector has no non-zero coefficients.", call. = FALSE)
}
rcpp_add_linear_constraint(
base$data$model_ptr,
j0 = as.integer(idx - 1L),
x = as.numeric(v2[idx]),
sense = "=",
rhs = as.numeric(eps),
name = as.character(name %||% ""),
block_name = block_name,
tag = tag
)
base$data$meta$model_dirty <- TRUE
base <- .pa_refresh_model_snapshot(base)
base
}
.pamo_compile_problem <- function(x, method_name = NULL, ...) {
stopifnot(inherits(x, "Problem"))
method <- x$data$method %||% NULL
if (is.null(method) || !is.list(method) || length(method) == 0L) {
stop(
"Multi-objective compilation requires a valid x$data$method configuration.",
call. = FALSE
)
}
if (is.null(method_name)) {
method_name <- as.character(method$type %||% method$name %||% NA_character_)[1]
}
if (is.na(method_name) || !nzchar(method_name)) {
stop(
"Invalid multi-objective method configuration: missing method name.",
call. = FALSE
)
}
aliases <- switch(
method_name,
weighted = as.character(method$aliases %||% character(0)),
epsilon_constraint = unique(c(
as.character(method$primary %||% character(0)),
as.character(method$constrained %||% character(0))
)),
augmecon = unique(c(
as.character(method$primary %||% character(0)),
as.character(method$secondary %||% character(0))
)),
stop(
"Unknown/unsupported multi-objective method: '", method_name, "'.",
call. = FALSE
)
)
aliases <- aliases[!is.na(aliases) & nzchar(aliases)]
if (length(aliases) == 0L) {
stop(
"Could not determine the active objective aliases for multi-objective compilation.",
call. = FALSE
)
}
specs <- .pamo_get_objective_specs(x, aliases)
ir_list <- lapply(specs, function(sp) .pamo_objective_to_ir(x, sp))
out <- .pamo_prepare_superset_model(x, ir_list)
out$data$meta <- out$data$meta %||% list()
out$data$meta$model_dirty <- FALSE
out$data$has_model <- TRUE
out$data$mo_cache <- out$data$mo_cache %||% list()
out$data$mo_cache$compiled_method <- method_name
out$data$mo_cache$compiled_aliases <- aliases
out
}
# -------------------------------------------------------------------------
# Weighted run compiler
# -------------------------------------------------------------------------
.pamo_compile_weighted_runs <- function(x, method) {
stopifnot(inherits(x, "Problem"))
aliases <- as.character(method$aliases %||% character(0))
if (length(aliases) == 0L) {
stop("Weighted method: missing objective aliases.", call. = FALSE)
}
runs <- method$runs %||% NULL
.pamo_check_run_design(runs)
if (.pamo_is_run_manual(runs)) {
return(.pamo_compile_weighted_manual_runs(
values = runs$values,
aliases = aliases,
normalize_weights = isTRUE(method$normalize_weights)
))
}
if (.pamo_is_run_grid(runs)) {
return(.pamo_compile_weighted_grid_runs(
aliases = aliases,
n = runs$n,
include_extremes = runs$include_extremes,
normalize_weights = isTRUE(method$normalize_weights)
))
}
stop("Unsupported weighted run design.", call. = FALSE)
}
.pamo_compile_weighted_manual_runs <- function(values,
aliases,
normalize_weights = TRUE) {
values <- as.data.frame(values, stringsAsFactors = FALSE)
weight_cols <- paste0("weight_", aliases)
miss <- setdiff(weight_cols, names(values))
if (length(miss) > 0L) {
stop(
"Manual weighted runs are missing columns: ",
paste(miss, collapse = ", "),
call. = FALSE
)
}
if (!("run_id" %in% names(values))) {
values$run_id <- seq_len(nrow(values))
}
out <- values[, c("run_id", weight_cols), drop = FALSE]
.pamo_validate_weighted_design(
out,
aliases = aliases,
normalize_weights = normalize_weights
)
}
.pamo_compile_weighted_grid_runs <- function(aliases,
n,
include_extremes = TRUE,
normalize_weights = TRUE) {
aliases <- as.character(aliases)
k <- length(aliases)
if (k < 2L) {
stop("Weighted run grids require at least two objectives.", call. = FALSE)
}
n <- as.integer(n)[1]
if (!is.finite(n) || is.na(n) || n < 2L) {
stop("`n` must be an integer >= 2.", call. = FALSE)
}
# 2 objectives: simple line grid
if (k == 2L) {
w1 <- seq(1, 0, length.out = n)
w2 <- 1 - w1
W <- cbind(w1, w2)
colnames(W) <- paste0("weight_", aliases)
if (!isTRUE(include_extremes)) {
keep <- W[, 1] > 0 & W[, 1] < 1
W <- W[keep, , drop = FALSE]
}
out <- data.frame(
run_id = seq_len(nrow(W)),
W,
check.names = FALSE,
stringsAsFactors = FALSE
)
return(.pamo_validate_weighted_design(
out,
aliases = aliases,
normalize_weights = normalize_weights
))
}
# 3+ objectives: regular simplex grid
# n controls the number of levels per weight, so denominator is n - 1.
denom <- n - 1L
comp <- .pamo_simplex_integer_compositions(total = denom, parts = k)
W <- comp / denom
colnames(W) <- paste0("weight_", aliases)
if (!isTRUE(include_extremes)) {
# Drop pure-objective vertices such as (1,0,0,...).
W <- W[apply(W, 1, max) < 1, , drop = FALSE]
}
out <- data.frame(
run_id = seq_len(nrow(W)),
W,
check.names = FALSE,
stringsAsFactors = FALSE
)
.pamo_validate_weighted_design(
out,
aliases = aliases,
normalize_weights = normalize_weights
)
}
.pamo_validate_weighted_design <- function(design,
aliases,
normalize_weights = TRUE,
tol = 1e-10) {
aliases <- as.character(aliases)
weight_cols <- paste0("weight_", aliases)
for (nm in weight_cols) {
design[[nm]] <- as.numeric(design[[nm]])
}
W <- as.matrix(design[, weight_cols, drop = FALSE])
if (any(!is.finite(W) | is.na(W))) {
stop("Weighted runs contain non-finite or missing weights.", call. = FALSE)
}
# if (any(W < -tol)) {
# stop("Weighted runs cannot contain negative weights.", call. = FALSE)
# }
row_sum <- rowSums(W)
# if (any(row_sum <= tol)) {
# stop("Each weighted run must contain at least one positive weight.", call. = FALSE)
# }
if (isTRUE(normalize_weights)) {
W <- sweep(W, 1, row_sum, "/")
}
design[, weight_cols] <- W
# Keep run_id first
design <- design[, c("run_id", weight_cols), drop = FALSE]
rownames(design) <- NULL
design
}
.pamo_simplex_integer_compositions <- function(total, parts) {
total <- as.integer(total)[1]
parts <- as.integer(parts)[1]
if (parts == 1L) {
return(matrix(total, nrow = 1L, ncol = 1L))
}
out <- list()
idx <- 1L
rec <- function(prefix, remaining, p_left) {
if (p_left == 1L) {
out[[idx]] <<- c(prefix, remaining)
idx <<- idx + 1L
return(invisible(NULL))
}
for (v in seq.int(remaining, 0L)) {
rec(c(prefix, v), remaining - v, p_left - 1L)
}
invisible(NULL)
}
rec(integer(0), total, parts)
mat <- do.call(rbind, out)
storage.mode(mat) <- "integer"
mat
}
# -------------------------------------------------------------------------
# Epsilon-constraint run compiler
# -------------------------------------------------------------------------
.pamo_compile_epsilon_runs <- function(x,
method,
gap_limit = NULL,
time_limit = NULL,
verbose = FALSE) {
stopifnot(inherits(x, "Problem"))
runs <- method$runs %||% NULL
.pamo_check_run_design(runs)
constrained <- as.character(method$constrained %||% character(0))
if (length(constrained) == 0L) {
stop("epsilon_constraint: missing constrained objective(s).", call. = FALSE)
}
if (.pamo_is_run_manual(runs)) {
return(.pamo_compile_epsilon_manual_runs(
values = runs$values,
constrained = constrained
))
}
if (.pamo_is_run_grid(runs)) {
# Reuse the existing automatic builder by passing n/include_extremes
# through the method object.
x2 <- x
x2$data$method <- method
x2$data$method$n_points <- runs$n
x2$data$method$include_extremes <- runs$include_extremes
return(.pamo_build_auto_epsilon_runs(
x = x2,
gap_limit = gap_limit,
time_limit = time_limit,
verbose = verbose
))
}
stop("Unsupported epsilon-constraint run design.", call. = FALSE)
}
.pamo_compile_epsilon_manual_runs <- function(values, constrained) {
values <- as.data.frame(values, stringsAsFactors = FALSE)
eps_cols <- paste0("eps_", constrained)
miss <- setdiff(eps_cols, names(values))
if (length(miss) > 0L) {
stop(
"Manual epsilon runs are missing columns: ",
paste(miss, collapse = ", "),
call. = FALSE
)
}
if (!("run_id" %in% names(values))) {
values$run_id <- seq_len(nrow(values))
}
out <- values[, c("run_id", eps_cols), drop = FALSE]
for (nm in eps_cols) {
out[[nm]] <- as.numeric(out[[nm]])
}
if (any(!is.finite(as.matrix(out[, eps_cols, drop = FALSE])))) {
stop("Manual epsilon runs contain non-finite epsilon values.", call. = FALSE)
}
rownames(out) <- NULL
out
}
# -------------------------------------------------------------------------
# AUGMECON run compiler
# -------------------------------------------------------------------------
.pamo_compile_augmecon_runs <- function(x,
method,
gap_limit = NULL,
time_limit = NULL,
verbose = FALSE) {
stopifnot(inherits(x, "Problem"))
runs <- method$runs %||% NULL
.pamo_check_run_design(runs)
secondary <- as.character(method$secondary %||% character(0))
if (length(secondary) == 0L) {
stop("augmecon: missing secondary objective(s).", call. = FALSE)
}
if (.pamo_is_run_manual(runs)) {
return(.pamo_compile_augmecon_manual_runs(
values = runs$values,
secondary = secondary
))
}
if (.pamo_is_run_grid(runs)) {
# Reuse the existing automatic builder by passing n/include_extremes
# through the method object.
x2 <- x
x2$data$method <- method
x2$data$method$n_points <- runs$n
x2$data$method$include_extremes <- runs$include_extremes
return(.pamo_build_auto_augmecon_runs(
x = x2,
gap_limit = gap_limit,
time_limit = time_limit,
verbose = verbose
))
}
stop("Unsupported AUGMECON run design.", call. = FALSE)
}
.pamo_compile_augmecon_manual_runs <- function(values, secondary) {
values <- as.data.frame(values, stringsAsFactors = FALSE)
eps_cols <- paste0("eps_", secondary)
miss <- setdiff(eps_cols, names(values))
if (length(miss) > 0L) {
stop(
"Manual AUGMECON runs are missing columns: ",
paste(miss, collapse = ", "),
call. = FALSE
)
}
if (!("run_id" %in% names(values))) {
values$run_id <- seq_len(nrow(values))
}
out <- values[, c("run_id", eps_cols), drop = FALSE]
for (nm in eps_cols) {
out[[nm]] <- as.numeric(out[[nm]])
}
if (any(!is.finite(as.matrix(out[, eps_cols, drop = FALSE])))) {
stop("Manual AUGMECON runs contain non-finite epsilon values.", call. = FALSE)
}
rownames(out) <- NULL
out
}
.pa_deprecate_arg <- function(old, new) {
warning(
"`", old, "` is deprecated. Use `", new, "` instead.",
call. = FALSE
)
}
.pamo_eps_to_manual_df <- function(eps, constrained) {
if (is.numeric(eps) && !is.null(names(eps))) {
eps_list <- as.list(eps)
eps_list <- lapply(eps_list, function(v) c(as.numeric(v)[1]))
} else if (is.list(eps) && length(eps) > 0L && !is.null(names(eps))) {
eps_list <- lapply(eps, function(v) as.numeric(v))
} else {
stop(
"`eps` must be a named numeric vector or a named list of numeric vectors.",
call. = FALSE
)
}
miss <- setdiff(constrained, names(eps_list))
if (length(miss) > 0L) {
stop(
"`eps` must include all constrained objectives. Missing: ",
paste(miss, collapse = ", "),
call. = FALSE
)
}
extra <- setdiff(names(eps_list), constrained)
if (length(extra) > 0L) {
stop(
"`eps` contains aliases that are not constrained objectives: ",
paste(extra, collapse = ", "),
call. = FALSE
)
}
bad_empty <- names(eps_list)[
vapply(eps_list, length, integer(1)) == 0L
]
if (length(bad_empty) > 0L) {
stop(
"`eps` contains empty vectors for: ",
paste(bad_empty, collapse = ", "),
call. = FALSE
)
}
bad_nonfinite <- names(eps_list)[
vapply(
eps_list,
function(v) any(!is.finite(v)),
logical(1)
)
]
if (length(bad_nonfinite) > 0L) {
stop(
"`eps` contains non-finite values for: ",
paste(bad_nonfinite, collapse = ", "),
call. = FALSE
)
}
grid <- expand.grid(
eps_list[constrained],
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
if (nrow(grid) == 0L) {
stop("Empty epsilon grid.", call. = FALSE)
}
names(grid) <- paste0("eps_", names(grid))
rownames(grid) <- NULL
grid
}
.pamo_augmecon_grid_to_manual_df <- function(grid, secondary) {
if (!is.list(grid) || length(grid) == 0L) {
stop("`grid` must be a non-empty named list.", call. = FALSE)
}
if (is.null(names(grid)) ||
anyNA(names(grid)) ||
any(!nzchar(names(grid)))) {
stop(
"`grid` must be a fully named list.",
call. = FALSE
)
}
if (anyDuplicated(names(grid)) != 0L) {
stop(
"`grid` names must not contain duplicates.",
call. = FALSE
)
}
missing_secondary <- setdiff(secondary, names(grid))
unknown_secondary <- setdiff(names(grid), secondary)
if (length(missing_secondary) > 0L) {
stop(
"`grid` is missing secondary objective(s): ",
paste(missing_secondary, collapse = ", "),
".",
call. = FALSE
)
}
if (length(unknown_secondary) > 0L) {
stop(
"`grid` contains unknown secondary objective(s): ",
paste(unknown_secondary, collapse = ", "),
".",
call. = FALSE
)
}
# Preserve the order of the secondary objectives
grid <- grid[secondary]
valid <- vapply(
grid,
function(z) {
is.numeric(z) &&
length(z) > 0L &&
!anyNA(z) &&
all(is.finite(z))
},
logical(1)
)
if (any(!valid)) {
bad <- names(grid)[!valid]
stop(
"Every element of `grid` must be a non-empty finite numeric vector. ",
"Invalid objective(s): ",
paste(bad, collapse = ", "),
".",
call. = FALSE
)
}
out <- expand.grid(
grid,
KEEP.OUT.ATTRS = FALSE,
stringsAsFactors = FALSE
)
names(out) <- paste0("eps_", names(out))
rownames(out) <- NULL
out
}
#' Get run ids for stored run-level results
#'
#' @noRd
.pa_get_stored_run_ids <- function(x) {
if (!inherits(x, "SolutionSet")) {
stop("x must be a SolutionSet object.", call. = FALSE)
}
sols <- x$solution$solutions %||% NULL
if (is.null(sols) || !is.list(sols) || length(sols) == 0L) {
return(integer(0))
}
nms <- names(sols)
if (!is.null(nms) && all(nzchar(nms))) {
out <- suppressWarnings(as.integer(nms))
if (!anyNA(out)) {
return(out)
}
}
# Fallback for older objects.
out <- vapply(seq_along(sols), function(i) {
sol <- sols[[i]]
as.integer(sol$meta$run_id %||% i)
}, integer(1))
out
}
#' Make sequential solution ids
#'
#' @noRd
.pa_make_solution_ids <- function(n, prefix = "s") {
n <- as.integer(n)[1]
if (is.na(n) || n < 0L) {
stop("`n` must be a non-negative integer.", call. = FALSE)
}
if (n == 0L) {
return(character(0))
}
paste0(prefix, seq_len(n))
}
#' Finalize solution ids in a SolutionSet
#'
#' @noRd
.pa_finalize_solution_ids <- function(x) {
if (!inherits(x, "SolutionSet")) {
stop("x must be a SolutionSet object.", call. = FALSE)
}
runs <- x$solution$runs %||% NULL
sols <- x$solution$solutions %||% NULL
if (is.null(runs) || !inherits(runs, "data.frame")) {
return(x)
}
if (!("run_id" %in% names(runs))) {
runs$run_id <- seq_len(nrow(runs))
}
runs$run_id <- as.integer(runs$run_id)
if (anyNA(runs$run_id) || any(runs$run_id < 1L)) {
stop(
"Run table contains invalid `run_id` values.",
call. = FALSE
)
}
# Public/internal convention:
# - solution_id is numeric.
# - solution_id matches run_id when a stored solution exists.
# - runs without stored solutions have solution_id = NA.
runs$solution_id <- NA_integer_
if (is.null(sols) || !is.list(sols) || length(sols) == 0L) {
x$solution$runs <- runs
x$solution$solutions <- list()
return(x)
}
is_sol <- vapply(
sols,
inherits,
logical(1),
"Solution"
)
sols_valid <- sols[is_sol]
if (length(sols_valid) == 0L) {
x$solution$runs <- runs
x$solution$solutions <- list()
return(x)
}
original_positions <- which(is_sol)
stored_run_ids <- vapply(seq_along(sols_valid), function(i) {
sol_i <- sols_valid[[i]]
rid <- sol_i$meta$run_id %||% NULL
if (is.null(rid)) {
original_i <- original_positions[i]
nm <- names(sols)[original_i] %||% NA_character_
rid <- suppressWarnings(as.integer(nm))
if (is.na(rid)) {
rid <- original_i
}
}
rid <- as.integer(rid)[1]
if (is.na(rid) || rid < 1L) {
NA_integer_
} else {
rid
}
}, integer(1))
keep <- !is.na(stored_run_ids) & stored_run_ids %in% runs$run_id
if (!all(keep)) {
sols_valid <- sols_valid[keep]
stored_run_ids <- stored_run_ids[keep]
}
if (length(sols_valid) == 0L) {
x$solution$runs <- runs
x$solution$solutions <- list()
return(x)
}
if (anyDuplicated(stored_run_ids) > 0L) {
duplicated_runs <- unique(stored_run_ids[duplicated(stored_run_ids)])
stop(
"Multiple stored solutions point to the same run_id: ",
paste(duplicated_runs, collapse = ", "),
".",
call. = FALSE
)
}
solution_ids <- as.integer(stored_run_ids)
for (i in seq_along(sols_valid)) {
rid <- stored_run_ids[i]
sid <- solution_ids[i]
idx <- match(rid, runs$run_id)
if (!is.na(idx)) {
runs$solution_id[idx] <- sid
}
sols_valid[[i]]$meta <- sols_valid[[i]]$meta %||% list()
sols_valid[[i]]$meta$run_id <- rid
sols_valid[[i]]$meta$solution_id <- sid
}
# List names must be character in R, but the stored metadata and run table
# keep solution_id as integer.
names(sols_valid) <- as.character(solution_ids)
x$solution$runs <- runs
x$solution$solutions <- sols_valid
# Add/update solution_id in summary tables when they have run_id.
if (!is.null(x$summary) && is.list(x$summary)) {
run_map <- runs[, c("run_id", "solution_id"), drop = FALSE]
for (nm in names(x$summary)) {
tab <- x$summary[[nm]]
if (
inherits(tab, "data.frame") &&
"run_id" %in% names(tab)
) {
# Replace old solution_id values such as "s1", "s2", ...
if ("solution_id" %in% names(tab)) {
tab$solution_id <- NULL
}
tab <- merge(
tab,
run_map,
by = "run_id",
all.x = TRUE,
sort = FALSE
)
first <- intersect(c("run_id", "solution_id"), names(tab))
rest <- setdiff(names(tab), first)
tab <- tab[, c(first, rest), drop = FALSE]
rownames(tab) <- NULL
x$summary[[nm]] <- tab
}
}
}
x
}
# -------------------------------------------------------------------------
# Internal run-design helpers
# -------------------------------------------------------------------------
.pamo_is_run_design <- function(x) {
inherits(x, "RunDesign")
}
.pamo_is_run_grid <- function(x) {
inherits(x, "RunGrid")
}
.pamo_is_run_manual <- function(x) {
inherits(x, "RunManual")
}
.pamo_check_run_design <- function(runs) {
if (missing(runs) || is.null(runs)) {
stop(
"`runs` must be created with `set_runs_grid()` or `set_runs_manual()`.",
call. = FALSE
)
}
if (!.pamo_is_run_design(runs)) {
stop(
"`runs` must be created with `set_runs_grid()` or `set_runs_manual()`.",
call. = FALSE
)
}
if (.pamo_is_run_grid(runs)) {
if (
!identical(runs$type, "grid") ||
length(runs$n) != 1L ||
is.na(runs$n) ||
!is.finite(runs$n) ||
runs$n < 2L ||
runs$n != floor(runs$n) ||
!is.logical(runs$include_extremes) ||
length(runs$include_extremes) != 1L ||
is.na(runs$include_extremes)
) {
stop(
"Invalid RunGrid object.",
call. = FALSE
)
}
return(invisible(TRUE))
}
if (.pamo_is_run_manual(runs)) {
if (
!identical(runs$type, "manual") ||
!inherits(runs$values, "data.frame") ||
nrow(runs$values) == 0L
) {
stop(
"Invalid RunManual object.",
call. = FALSE
)
}
return(invisible(TRUE))
}
stop(
"Unsupported RunDesign subclass.",
call. = FALSE
)
}
#' Validate and resolve a multi-objective control object
#'
#' @param control A control object created with \code{mo_control()}, or
#' \code{NULL} to use the default settings.
#'
#' @return A validated object of class \code{MOControl}.
#'
#' @noRd
.pamo_check_mo_control <- function(control) {
if (is.null(control)) {
return(set_runs_control())
}
if (!inherits(control, "MOControl") && !inherits(control, "RunsControl")) {
stop(
"`control` must be created with set_runs_control().",
call. = FALSE
)
}
required <- c(
"stop_on_infeasible",
"stop_on_no_solution",
"stop_on_error"
)
missing <- setdiff(required, names(control))
if (length(missing) > 0L) {
stop(
"Invalid control object. Missing field(s): ",
paste(missing, collapse = ", "),
".",
call. = FALSE
)
}
for (nm in required) {
value <- control[[nm]]
if (
!is.logical(value) ||
length(value) != 1L ||
is.na(value)
) {
stop(
"Invalid control object. `", nm, "` must be TRUE or FALSE.",
call. = FALSE
)
}
}
structure(
list(
stop_on_infeasible = isTRUE(control$stop_on_infeasible),
stop_on_no_solution = isTRUE(control$stop_on_no_solution),
stop_on_error = isTRUE(control$stop_on_error)
),
class = c("RunsControl", "MOControl", "list")
)
}
#' Drop Problem objects stored inside run-level Solution objects
#'
#' @description
#' Internal helper used to reduce the memory footprint of SolutionSet objects.
#' The full Problem is stored once at SolutionSet$problem, so keeping an
#' additional copy inside each run-level Solution is unnecessary.
#'
#' @noRd
.pa_strip_nested_solution_problems <- function(solutions) {
if (is.null(solutions)) {
return(solutions)
}
if (!is.list(solutions)) {
return(solutions)
}
for (i in seq_along(solutions)) {
if (inherits(solutions[[i]], "Solution")) {
solutions[[i]]$problem <- NULL
}
}
solutions
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.