R/build_model.R

In multiscape: Multi-Action, Multi-Objective Spatial Planning

#' @include internal.R
#'
#' @title Build optimization model from Problem
#'
#' @description
#' Materializes (builds) the optimization model using the current state of the `Problem` object:
#' prepared data tables, stored objective settings, and stored constraints (e.g., targets).
#'
#' @param x Problem object (class "Problem") created with create_problem().
#'
#' @return Updated `Problem` object with model pointer and model snapshot.
#' @keywords internal
.pa_build_model <- function(x) {

  stopifnot(inherits(x, "Problem"))
  stopifnot(!is.null(x$data$pu))
  stopifnot(!is.null(x$data$dist_features))

  # ------------------------------------------------------------
  # local helpers (keep minimal here; you can move to internal.R)
  # ------------------------------------------------------------
  .pa_abort <- function(...) stop(paste0(...), call. = FALSE)
  .pa_has_rows <- function(df) !is.null(df) && inherits(df, "data.frame") && nrow(df) > 0

  .pa_log_init <- function(x) {
    x$data$meta <- x$data$meta %||% list()
    x$data$meta$build_log <- x$data$meta$build_log %||% list(
      dropped_status3 = 0L,
      dropped_missing_internal = 0L,
      dropped_bad_cost = 0L,
      filtered_effects_rows = 0L,
      filtered_profit_rows = 0L
    )
    x
  }
  .pa_log_add <- function(x, key, n) {
    x$data$meta$build_log[[key]] <- as.integer(x$data$meta$build_log[[key]] %||% 0L) + as.integer(n)
    x
  }
  .pa_emit_build_warnings <- function(x) {
    log <- x$data$meta$build_log %||% list()
    msgs <- character(0)

    if (isTRUE(log$dropped_status3 > 0L)) {
      msgs <- c(msgs, paste0("Dropped ", log$dropped_status3, " dist_actions rows with status==3 (locked-out)."))
    }
    if (isTRUE(log$dropped_missing_internal > 0L)) {
      msgs <- c(msgs, paste0("Dropped ", log$dropped_missing_internal, " dist_actions rows with missing internal_pu/internal_action."))
    }
    if (isTRUE(log$dropped_bad_cost > 0L)) {
      msgs <- c(msgs, paste0("Dropped ", log$dropped_bad_cost, " dist_actions rows with missing/non-finite cost."))
    }
    if (isTRUE(log$filtered_effects_rows > 0L)) {
      msgs <- c(msgs, paste0("Filtered ", log$filtered_effects_rows, " dist_effects rows not matching feasible dist_actions."))
    }
    if (isTRUE(log$filtered_profit_rows > 0L)) {
      msgs <- c(msgs, paste0("Filtered ", log$filtered_profit_rows, " dist_profit rows not matching feasible dist_actions."))
    }

    if (length(msgs) > 0) {
      warning(paste(msgs, collapse = "\n"), call. = FALSE, immediate. = TRUE)
    }
    x
  }

  .pa_check_single_objective <- function(args) {
    if (isTRUE(args$mo_mode)) return(invisible(TRUE))

    if (is.null(args)) return(invisible(TRUE))

    if (!is.null(args$objectives)) {
      nobj <- if (is.data.frame(args$objectives)) nrow(args$objectives) else length(args$objectives)
      if (is.finite(nobj) && nobj > 1) {
        .pa_abort(
          "Multiple objectives detected (", nobj, "). multiscape builds a single-objective MILP.\n",
          "Use a multiobjective method for epsilon-constraint/AUGMECON/interactive methods (e.g., set_method_weighted)."
        )
      }
    }

    if (is.vector(args$model_type) && length(args$model_type) > 1) {
      .pa_abort(
        "Multiple model_type values detected: ",
        paste(args$model_type, collapse = ", "),
        ". Choose exactly one objective."
      )
    }

    invisible(TRUE)
  }

  # ------------------------------------------------------------
  # init build log and args
  # ------------------------------------------------------------
  x <- .pa_log_init(x)

  input_format <- x$data$meta$input_format %||% "new"

  if (is.null(x$data$model_args)) x$data$model_args <- list()
  .pa_check_single_objective(x$data$model_args)

  # objective required for new pipeline (legacy can default)
  if (is.null(x$data$model_args$model_type)) {
    if (identical(input_format, "legacy")) {
      x$data$model_args$model_type <- "minimizeCosts"
    } else {
      .pa_abort(
        "No objective configured. Use add_objective_*() (new pipeline)."
      )
    }
  }

  # defaults for later printing / adapters
  if (is.null(x$data$model_args$budget))   x$data$model_args$budget   <- 0
  if (is.null(x$data$model_args$blm))      x$data$model_args$blm      <- 0
  if (is.null(x$data$model_args$curve))    x$data$model_args$curve    <- 1L
  if (is.null(x$data$model_args$segments)) x$data$model_args$segments <- 3L

  # ---- NEW: initialize + set needs flags (single-objective safe, MO-safe)
  x <- .pa_build_model_set_needs_from_objective(x)

  # ------------------------------------------------------------
  # implicit simple conservation model
  # ------------------------------------------------------------
  # If the user supplied planning units, features and dist_features, but did
  # not supply actions or effects, create the classical conservation-planning
  # model internally:
  #
  #   one implicit conservation action per planning unit
  #   effect = 0
  #   action outcome = reference feature amount
  #
  # This must happen before ensure_tables(), prepare_tables(), targets, and
  # objective validation.
  x <- .pa_build_model_ensure_implicit_conservation(x)


  # ------------------------------------------------------------
  # legacy adapters
  # ------------------------------------------------------------
  # if (identical(input_format, "legacy")) {
  #   x <- .pa_build_model_apply_legacy_adapters(x)
  # }

  # ------------------------------------------------------------
  # ensure tables exist (allow "no actions" case)
  # ------------------------------------------------------------
  x <- .pa_build_model_ensure_tables(x)

  # ------------------------------------------------------------
  # post-adapter validation
  # ------------------------------------------------------------
  x <- .pa_build_model_validate_pipeline_state(x, input_format = input_format)

  # ------------------------------------------------------------
  # prepare model-ready tables (filters + joins)
  # ------------------------------------------------------------
  x <- .pa_build_model_prepare_tables(x)
  x <- .pa_prepare_joint_effects_model(x)

  x <- .pa_build_model_validate_locked_in_action_feasibility(x)
  x <- .pa_validate_action_cardinality_model(x)
  .pa_validate_action_relations_specs(x$data$constraints$action_relations,
                                     x$data$pu, x$data$actions)

  # ------------------------------------------------------------
  # early validation: objective dependencies
  # ------------------------------------------------------------
  x <- .pa_build_model_validate_objective_requirements(x)

  # ------------------------------------------------------------
  # build C++ model pointer + base vars + minimal constraints
  # ------------------------------------------------------------
  x <- .pa_build_model_build_cpp_core(x)

  # ------------------------------------------------------------
  # prepare auxiliary variables/constraints required by needs
  # ------------------------------------------------------------
  x <- .pa_build_model_prepare_needs_cpp(x)
  x <- .pa_build_joint_effects(x)
  x <- .pa_build_effect_losses(x)

  # ------------------------------------------------------------
  # objective (C++ side)
  # ------------------------------------------------------------
  x <- .pa_build_model_set_objective_cpp(x)

  # ------------------------------------------------------------
  # apply stored targets/constraints
  # ------------------------------------------------------------
  x <- .pa_build_model_apply_constraints(x)

  # ------------------------------------------------------------
  # store metadata + snapshot + warnings
  # ------------------------------------------------------------
  x <- .pa_build_model_finalize(x)
  x <- .pa_emit_build_warnings(x)

  x$data$has_model <- TRUE
  x
}

# -------------------------------------------------------------------------
# Helpers (recommended to live in internal.R, but ok here as noRd)
# -------------------------------------------------------------------------

.pa_build_model_ensure_tables <- function(x) {

  if (is.null(x$data$dist_actions)) {
    x$data$dist_actions <- data.frame(
      pu = integer(0),
      action = character(0),
      cost = numeric(0),
      status = integer(0),
      internal_pu = integer(0),
      internal_action = integer(0),
      stringsAsFactors = FALSE
    )
  }
  if (is.null(x$data$actions)) {
    x$data$actions <- data.frame(
      id = character(0),
      internal_id = integer(0),
      stringsAsFactors = FALSE
    )
  }
  x
}

.pa_build_model_validate_pipeline_state <- function(x, input_format) {

  mtype_check <- x$data$model_args$model_type %||% "minimizeCosts"

  if (!identical(input_format, "legacy") && identical(mtype_check, "minimizeCosts")) {
    has_targets <- !is.null(x$data$targets) &&
      inherits(x$data$targets, "data.frame") &&
      nrow(x$data$targets) > 0L

    area_constraints <- x$data$constraints$area %||% NULL
    has_selection_area <- !is.null(area_constraints) &&
      inherits(area_constraints, "data.frame") &&
      nrow(area_constraints) > 0L &&
      all(c("sense", "value") %in% names(area_constraints)) &&
      any(
        area_constraints$sense %in% c("min", "equal") &
          is.finite(area_constraints$value) &
          area_constraints$value > 0,
        na.rm = TRUE
      )

    pu <- x$data$pu %||% NULL
    has_locked_in_pu <- !is.null(pu) &&
      inherits(pu, "data.frame") &&
      "locked_in" %in% names(pu) &&
      any(as.logical(pu$locked_in), na.rm = TRUE)

    dist_actions <- x$data$dist_actions %||% NULL
    has_locked_in_action <- !is.null(dist_actions) &&
      inherits(dist_actions, "data.frame") &&
      "status" %in% names(dist_actions) &&
      any(dist_actions$status %in% c(1L, 2L), na.rm = TRUE)

    cardinality <- x$data$constraints$action_cardinality
    has_selection_cardinality <- is.data.frame(cardinality) && nrow(cardinality) > 0L &&
      any(cardinality$sense %in% c("min", "equal") & cardinality$count > 0L)

    has_selection_requirement <- any(c(
      has_targets,
      has_selection_area,
      has_locked_in_pu,
      has_locked_in_action,
      has_selection_cardinality
    ))

    if (!isTRUE(has_selection_requirement)) {
      warning(
        paste0(
          "The minimum-cost problem has no feature targets, positive minimum/equality area or action-count constraint, ",
          "or locked-in decisions. The all-zero solution may therefore be optimal. ",
          "Add a selection requirement if an empty solution is not intended."
        ),
        call. = FALSE
      )
    }
  }

  x
}
.pa_build_model_prepare_tables <- function(x) {

  .pa_abort <- function(...) stop(paste0(...), call. = FALSE)
  .has_rows <- function(df) !is.null(df) && inherits(df, "data.frame") && nrow(df) > 0

  da <- x$data$dist_actions
  if (is.null(da)) da <- data.frame()

  if (!inherits(da, "data.frame")) {
    .pa_abort("dist_actions must be a data.frame (or NULL).")
  }

  n0 <- nrow(da)

  if (n0 > 0) {

    # status filter
    if ("status" %in% names(da)) {
      keep <- da$status != 3L
      dropped <- sum(!keep, na.rm = TRUE)
      if (dropped > 0) x <- .pa_log_add(x, "dropped_status3", dropped)
      da <- da[keep, , drop = FALSE]
    }

    # internal ids
    if (all(c("internal_pu", "internal_action") %in% names(da))) {
      keep <- !is.na(da$internal_pu) & !is.na(da$internal_action)
      dropped <- sum(!keep)
      if (dropped > 0) x <- .pa_log_add(x, "dropped_missing_internal", dropped)
      da <- da[keep, , drop = FALSE]
    }

    # costs
    if ("cost" %in% names(da)) {
      da$cost <- as.numeric(da$cost)
      keep <- is.finite(da$cost) & !is.na(da$cost)
      dropped <- sum(!keep)
      if (dropped > 0) x <- .pa_log_add(x, "dropped_bad_cost", dropped)
      da <- da[keep, , drop = FALSE]
    }

    # normalize types
    if ("pu" %in% names(da)) da$pu <- as.integer(da$pu)
    if ("action" %in% names(da)) da$action <- as.character(da$action)
    if ("internal_pu" %in% names(da)) da$internal_pu <- as.integer(da$internal_pu)
    if ("internal_action" %in% names(da)) da$internal_action <- as.integer(da$internal_action)

    # deterministic order
    if (all(c("internal_pu","internal_action") %in% names(da)) && nrow(da) > 0) {
      da <- da[order(da$internal_pu, da$internal_action), , drop = FALSE]
    }
  }

  # internal_row
  da$internal_row <- seq_len(nrow(da))

  # validate internal id ranges (if present)
  if (.has_rows(da) && all(c("internal_pu","internal_action") %in% names(da))) {
    n_pu <- nrow(x$data$pu)
    if (any(da$internal_pu < 1L | da$internal_pu > n_pu, na.rm = TRUE)) {
      .pa_abort("dist_actions contains internal_pu indices out of range (1..n_pu).")
    }
    if (!is.null(x$data$actions) && inherits(x$data$actions, "data.frame") && nrow(x$data$actions) > 0) {
      n_a <- nrow(x$data$actions)
      if (any(da$internal_action < 1L | da$internal_action > n_a, na.rm = TRUE)) {
        .pa_abort("dist_actions contains internal_action indices out of range (1..n_actions).")
      }
    }
  }

  # duplicates check (pu, action)
  if (.has_rows(da) && all(c("pu","action") %in% names(da))) {
    key <- paste(da$pu, da$action, sep = "||")
    if (anyDuplicated(key)) {
      .pa_abort(
        "dist_actions has duplicated (pu, action) pairs. Please aggregate/resolve duplicates before building the model."
      )
    }
  }

  x$data$dist_actions_model <- da

  # dist_effects model-ready
  de <- x$data$dist_effects

  if (is.null(de)) {
    de <- data.frame(
      pu = integer(0),
      action = character(0),
      feature = integer(0),
      benefit = numeric(0),
      loss = numeric(0),
      internal_pu = integer(0),
      internal_action = integer(0),
      internal_feature = integer(0),
      stringsAsFactors = FALSE
    )
  }

  if (!inherits(de, "data.frame")) {
    .pa_abort("dist_effects must be a data.frame (or NULL).")
  }

  if (.has_rows(de) && .has_rows(da)) {

    n_before <- nrow(de)

    if ("pu" %in% names(de))      de$pu      <- as.integer(de$pu)
    if ("action" %in% names(de))  de$action  <- as.character(de$action)
    if ("feature" %in% names(de)) de$feature <- as.integer(de$feature)

    de <- dplyr::inner_join(
      de,
      da[, c("pu", "action"), drop = FALSE],
      by = c("pu", "action")
    )

    n_after <- nrow(de)
    dropped <- n_before - n_after
    if (dropped > 0) x <- .pa_log_add(x, "filtered_effects_rows", dropped)

    pu_map <- x$data$pu[, c("id", "internal_id")]
    act_map <- x$data$actions[, c("id", "internal_id")]
    feat_map <- x$data$features[, c("id", "internal_id")]

    if (!("internal_pu" %in% names(de))) {
      de$internal_pu <- pu_map$internal_id[match(de$pu, pu_map$id)]
    } else {
      de$internal_pu <- as.integer(de$internal_pu)
    }

    if (!("internal_action" %in% names(de))) {
      de$internal_action <- act_map$internal_id[match(de$action, act_map$id)]
    } else {
      de$internal_action <- as.integer(de$internal_action)
    }

    if (!("internal_feature" %in% names(de))) {
      if (!("feature" %in% names(de))) {
        .pa_abort("dist_effects must contain column 'feature' (feature id) to derive internal_feature.")
      }
      de$internal_feature <- feat_map$internal_id[match(de$feature, feat_map$id)]
    } else {
      de$internal_feature <- as.integer(de$internal_feature)
    }

    if (anyNA(de$internal_pu)) {
      .pa_abort("dist_effects has pu ids not found in x$data$pu$id (cannot derive internal_pu).")
    }
    if (anyNA(de$internal_action)) {
      .pa_abort("dist_effects has action ids not found in x$data$actions$id (cannot derive internal_action).")
    }
    if (anyNA(de$internal_feature)) {
      .pa_abort("dist_effects has feature ids not found in x$data$features$id (cannot derive internal_feature).")
    }
  }

  if (.has_rows(de)) {
    de <- .pa_add_feature_labels(
      df = de,
      features_df = x$data$features,
      feature_col = "feature",
      internal_feature_col = "internal_feature",
      out_col = "feature_name"
    )

    de <- .pa_add_action_labels(
      df = de,
      actions_df = x$data$actions,
      action_col = "action",
      internal_action_col = "internal_action",
      out_col = "action_name"
    )
  }

  x$data$dist_effects_model <- de



  # dist_benefit_model derived (now has internal_* too)
  db <- NULL
  if (.has_rows(de) && "benefit" %in% names(de)) {
    db <- de
    db$benefit <- as.numeric(db$benefit)
  }
  x$data$dist_benefit_model <- db

  # dist_profit model-ready (MUST carry internal_pu/internal_action)
  # dist_profit model-ready (MUST carry internal_pu/internal_action)
  dp <- x$data$dist_profit
  if (.has_rows(dp) && .has_rows(da)) {

    n_before <- nrow(dp)

    # tipos base
    if ("pu" %in% names(dp)) dp$pu <- as.integer(dp$pu)
    if ("action" %in% names(dp)) dp$action <- as.character(dp$action)

    # filtra a (pu, action) factibles
    dp <- dplyr::inner_join(
      dp,
      da[, c("pu", "action"), drop = FALSE],
      by = c("pu", "action")
    )

    n_after <- nrow(dp)
    dropped <- n_before - n_after
    if (dropped > 0) x <- .pa_log_add(x, "filtered_profit_rows", dropped)

    # --- NEW: asegurar columnas internas (igual que en dist_effects)
    pu_map  <- x$data$pu[, c("id", "internal_id")]
    act_map <- x$data$actions[, c("id", "internal_id")]

    pu_map$id  <- as.integer(pu_map$id)
    act_map$id <- as.character(act_map$id)

    if (!("internal_pu" %in% names(dp))) {
      dp$internal_pu <- pu_map$internal_id[match(dp$pu, pu_map$id)]
    } else {
      dp$internal_pu <- as.integer(dp$internal_pu)
    }

    if (!("internal_action" %in% names(dp))) {
      dp$internal_action <- act_map$internal_id[match(dp$action, act_map$id)]
    } else {
      dp$internal_action <- as.integer(dp$internal_action)
    }

    if (anyNA(dp$internal_pu)) {
      .pa_abort("dist_profit has pu ids not found in x$data$pu$id (cannot derive internal_pu).")
    }
    if (anyNA(dp$internal_action)) {
      .pa_abort("dist_profit has action ids not found in x$data$actions$id (cannot derive internal_action).")
    }

    # opcional: orden determinista
    dp <- dp[order(dp$internal_pu, dp$internal_action), , drop = FALSE]
  }

  x$data$dist_profit_model <- dp


  x
}

.pa_build_model_validate_objective_requirements <- function(x) {

  .pa_abort <- function(...) stop(paste0(...), call. = FALSE)
  .has_rows <- function(df) !is.null(df) && inherits(df, "data.frame") && nrow(df) > 0

  args <- x$data$model_args %||% list()

  if (is.null(args$model_type) || !nzchar(as.character(args$model_type)[1])) {
    .pa_abort("No active objective has been defined in x$data$model_args$model_type.")
  }

  mtype <- as.character(args$model_type)[1]
  oargs <- args$objective_args %||% list()

  if (.pa_uses_aggregate_effects(x) && mtype %in% c("maximizeBenefits", "minimizeLosses")) {
    .pa_joint_objective_terms(x, if (mtype == "minimizeLosses") "loss" else oargs$benefit_col %||% "effect",
                             oargs$actions, oargs$features)
    return(x)
  }

  has_actions_model <- .has_rows(x$data$dist_actions_model)
  has_effects_model <- .has_rows(x$data$dist_effects_model)
  has_targets <- !is.null(x$data$targets) &&
    inherits(x$data$targets, "data.frame") &&
    nrow(x$data$targets) > 0


  needs_actions <- mtype %in% c(
    "maximizeBenefits",
    "minimizeLosses",
    "maximizeProfit",
    "maximizeNetProfit",
    "minimizeActionFragmentation",
    "minimizeInterventionFragmentation",
    "minimizeInterventionImpact"
  )

  if (needs_actions && !isTRUE(has_actions_model)) {
    .pa_abort(
      "Objective '", mtype, "' requires actions, but no actions are available.\n",
      "Run add_actions() (and ensure feasible dist_actions rows exist) before solve()."
    )
  }

  # ------------------------------------------------------------
  # Actions + targets over features require effects
  # ------------------------------------------------------------
  if (isTRUE(has_actions_model) && isTRUE(has_targets) && !isTRUE(has_effects_model)) {
    .pa_abort(
      "This problem includes explicit actions and feature targets, but no action effects were provided.\n",
      "Because targets are defined on features, multiscape needs to know how each explicit action contributes to each feature.\n",
      "Run add_effects() after add_actions() before solve().\n",
      "If you intended to use the simple conservation-planning model, do not call add_actions(); ",
      "multiscape will then create an implicit conservation action automatically."
    )
  }

  if (identical(mtype, "maximizeBenefits")) {

    if (!.has_rows(x$data$dist_effects_model)) {
      .pa_abort(
        "Objective 'maximizeBenefits' requires effects, but dist_effects is empty.\n",
        "In explicit action-based models, run add_effects() after add_actions().\n",
        "In the simple conservation model, omit both add_actions() and add_effects() ",
        "so multiscape can create implicit conservation effects."
      )
    }

    bcol <- as.character(oargs$benefit_col %||% "effect")[1]

    if (!(bcol %in% names(x$data$dist_effects_model))) {
      .pa_abort(
        "Objective 'maximizeBenefits' requires column '", bcol, "' in dist_effects."
      )
    }

    de <- x$data$dist_effects_model

    # Apply the same action/feature filters that are later used by
    # .pa_build_model_set_objective_cpp().
    acts <- oargs$actions %||% NULL
    if (!is.null(acts)) {
      de <- de[de$internal_action %in% as.integer(acts), , drop = FALSE]
    }

    feats <- oargs$features %||% NULL
    if (!is.null(feats)) {
      de <- de[de$internal_feature %in% as.integer(feats), , drop = FALSE]
    }

    if (!.has_rows(de)) {
      .pa_abort("Objective 'maximizeBenefits' has no matching effect rows for its action/feature subset.")
    }
    bb <- as.numeric(de[[bcol]])
    if (anyNA(bb) || any(!is.finite(bb))) {
      .pa_abort("Objective 'maximizeBenefits' coefficient column contains NA or non-finite values.")
    }
    # Negative and identically zero net objectives are valid.
  }

  if (identical(mtype, "minimizeLosses")) {
    if (!.has_rows(x$data$dist_effects_model)) {
      .pa_abort(
        "Objective 'minimizeLosses' requires effects/losses, but dist_effects is empty.\n",
        "Run add_losses() / add_effects() after add_actions()."
      )
    }
    lcol <- as.character(oargs$loss_col %||% "loss")[1]
    if (!(lcol %in% names(x$data$dist_effects_model))) {
      .pa_abort(
        "Objective 'minimizeLosses' requires column '", lcol, "' in dist_effects."
      )
    }
  }

  if (mtype %in% c("maximizeProfit", "maximizeNetProfit")) {
    if (!.has_rows(x$data$dist_profit_model)) {
      .pa_abort(
        "Objective '", mtype, "' requires dist_profit, but dist_profit is empty.\n",
        "Run add_profit() after add_actions()."
      )
    }
    pcol <- as.character(oargs$profit_col %||% "profit")[1]
    if (!(pcol %in% names(x$data$dist_profit_model))) {
      .pa_abort(
        "Objective '", mtype, "' requires column '", pcol, "' in dist_profit."
      )
    }
  }

  # ------------------------------------------------------------
  # Fragmentation objectives: require spatial relation
  # ------------------------------------------------------------
  if (mtype %in% c(
    "minimizeFragmentation",
    "minimizeActionFragmentation",
    "minimizeInterventionFragmentation"
  )) {

    rel_name <- as.character(oargs$relation_name %||% "boundary")[1]
    rels <- x$data$spatial_relations

    if (is.null(rels) || !is.list(rels) || is.null(rels[[rel_name]])) {
      .pa_abort(
        "Objective '", mtype, "' requires a spatial relation named '", rel_name, "'.\n",
        "But x$data$spatial_relations[['", rel_name, "']] is missing."
      )
    }

    rel <- rels[[rel_name]]
    if (!inherits(rel, "data.frame") || nrow(rel) == 0) {
      .pa_abort("Spatial relation '", rel_name, "' exists but is empty or not a data.frame.")
    }

    need_cols <- c("internal_pu1", "internal_pu2", "weight")
    if (!all(need_cols %in% names(rel))) {
      .pa_abort(
        "Spatial relation '", rel_name, "' must contain columns: ",
        paste(need_cols, collapse = ", "), "."
      )
    }

    n_pu <- nrow(x$data$pu)
    if (anyNA(rel$internal_pu1) || anyNA(rel$internal_pu2)) {
      .pa_abort("Spatial relation '", rel_name, "' contains NA in internal_pu1/internal_pu2.")
    }
    if (any(rel$internal_pu1 < 1L | rel$internal_pu1 > n_pu) ||
        any(rel$internal_pu2 < 1L | rel$internal_pu2 > n_pu)) {
      .pa_abort("Spatial relation '", rel_name, "' has internal PU indices out of range (1..n_pu).")
    }

    ww <- as.numeric(rel$weight)
    # if (any(!is.finite(ww)) || any(ww < 0)) {
    #   .pa_abort("Spatial relation '", rel_name, "' has non-finite or negative weights.")
    # }

    if (identical(mtype, "minimizeActionFragmentation")) {
      if (is.null(x$data$actions) || !inherits(x$data$actions, "data.frame") || nrow(x$data$actions) == 0) {
        .pa_abort("Objective 'minimizeActionFragmentation' requires x$data$actions to exist (non-empty).")
      }
      if (!is.null(oargs$actions_to_use)) {
        a <- as.integer(oargs$actions_to_use)
        if (anyNA(a)) .pa_abort("actions_to_use contains NA.")
        if (any(a < 1L | a > nrow(x$data$actions))) {
          .pa_abort("actions_to_use has internal_action ids out of range (1..n_actions).")
        }
      }
      # if (!is.null(oargs$action_weights)) {
      #   aw <- as.numeric(oargs$action_weights)
      #   if (any(!is.finite(aw)) || any(aw < 0)) {
      #     .pa_abort("action_weights must be finite and >= 0.")
      #   }
      # }
    }
  }

  # ------------------------------------------------------------
  # Intervention impact: does NOT require spatial relation
  # ------------------------------------------------------------
  if (identical(mtype, "minimizeInterventionImpact")) {

    if (!.has_rows(x$data$dist_features)) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact' requires dist_features, but dist_features is empty."
      )
    }

    icol <- as.character(oargs$impact_col %||% "amount")[1]
    if (!(icol %in% names(x$data$dist_features))) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact' requires column '", icol, "' in dist_features."
      )
    }

    if (is.null(x$data$actions) || !inherits(x$data$actions, "data.frame") || nrow(x$data$actions) == 0) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact' requires x$data$actions to exist (non-empty)."
      )
    }

    acts <- oargs$actions %||% NULL
    if (is.null(acts)) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact' requires `actions=` in objective_args.\n",
        "Example: add_objective_min_intervention_impact(actions = c('restoration'))."
      )
    }

    act_subset <- tryCatch(
      .pa_resolve_action_subset(x, subset = acts),
      error = function(e) .pa_abort(conditionMessage(e))
    )

    if (!inherits(act_subset, "data.frame") || nrow(act_subset) == 0) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact': action subset resolved to zero actions."
      )
    }

    feats <- oargs$features %||% NULL
    if (!is.null(feats)) {
      feat_subset <- tryCatch(
        .pa_resolve_feature_subset(x, features = feats),
        error = function(e) .pa_abort(conditionMessage(e))
      )

      if (!inherits(feat_subset, "data.frame") || nrow(feat_subset) == 0) {
        .pa_abort(
          "Objective 'minimizeInterventionImpact': feature subset resolved to zero features."
        )
      }
    }
  }

  x
}

.pa_build_model_validate_locked_in_action_feasibility <- function(x) {

  stopifnot(inherits(x, "Problem"))

  pu <- x$data$pu
  da <- x$data$dist_actions_model

  if (is.null(pu) || !inherits(pu, "data.frame") || nrow(pu) == 0) {
    return(x)
  }

  if (!("locked_in" %in% names(pu)) || !any(pu$locked_in, na.rm = TRUE)) {
    return(x)
  }

  locked_in_pu <- as.integer(pu$internal_id[isTRUE(pu$locked_in) | (!is.na(pu$locked_in) & pu$locked_in)])

  if (length(locked_in_pu) == 0) return(x)

  feasible_pu <- integer(0)
  if (!is.null(da) && inherits(da, "data.frame") && nrow(da) > 0 && "internal_pu" %in% names(da)) {
    feasible_pu <- unique(as.integer(da$internal_pu))
  }

  bad <- setdiff(locked_in_pu, feasible_pu)

  if (length(bad) > 0) {
    bad_ext <- pu$id[match(bad, pu$internal_id)]
    stop(
      "Some locked-in planning units have no feasible actions, but the model requires w_i = 1 ",
      "to imply at least one action.\n",
      "Affected PU ids: ", paste(utils::head(bad_ext, 20), collapse = ", "),
      if (length(bad_ext) > 20) paste0(" ... (", length(bad_ext), " total)") else "",
      "\nFix this by either:\n",
      "- adding at least one feasible action to those units,\n",
      "- removing locked_in from those units, or\n",
      "- making 'Conservation' an explicit feasible action there.",
      call. = FALSE
    )
  }

  x
}

.pa_build_model_build_cpp_core <- function(x) {


  .pa_abort <- function(...) stop(paste0(...), call. = FALSE)

  if (!exists("rcpp_new_optimization_problem", mode = "function")) {
    .pa_abort("Missing rcpp_new_optimization_problem() in the package.")
  }
  if (!exists("rcpp_add_base_variables", mode = "function")) {
    .pa_abort("Missing rcpp_add_base_variables() in the package.")
  }

  # ---- NEW: read needs flags (default-safe)
  args  <- x$data$model_args %||% list()
  needs <- args$needs %||% list()
  need_z <- isTRUE(needs$z)

  # Reserve for the actual core rather than a million slots in every
  # temporary MO model. Native vectors grow when auxiliaries are added.
  core_columns <- nrow(x$data$pu) + nrow(x$data$dist_actions_model)
  op <- rcpp_new_optimization_problem(
    nrow = as.integer(max(4096, 2 * core_columns)),
    ncol = as.integer(max(4096, core_columns)),
    ncell = as.integer(max(100000, 4 * core_columns))
  )

  # registry placeholder for future MO updates (constraint/objective IDs)
  x$data$model_registry <- list(
    cons = list(),
    vars = list(),
    obj_templates = list(),
    objective = list()
  )

  idx <- rcpp_add_base_variables(
    op,
    pu_data            = x$data$pu,
    dist_actions_data  = x$data$dist_actions_model,
    dist_features_data = x$data$dist_features,
    add_z              = need_z            # <- NEW
  )

  # structural constraints
  if (!is.null(x$data$dist_actions_model) && nrow(x$data$dist_actions_model) > 0) {

    if (!exists("rcpp_add_linking_x_le_w", mode = "function")) {
      .pa_abort("Missing rcpp_add_linking_x_le_w() in the package.")
    }
    res_locks <- rcpp_add_linking_x_le_w(op, x$data$dist_actions_model)
    x$data$model_registry$cons$x_le_w <- res_locks

    # NEW: w_i <= sum_a x_ia
    if (!exists("rcpp_add_linking_w_le_sum_x", mode = "function")) {
      .pa_abort("Missing rcpp_add_linking_w_le_sum_x() in the package.")
    }
    res_locks <- rcpp_add_linking_w_le_sum_x(op, x$data$dist_actions_model)
    x$data$model_registry$cons$w_le_sum_x <- res_locks
  }

  # only add z <= w if z exists
  if (isTRUE(need_z) && exists("rcpp_add_linking_z_le_w", mode = "function")) {
    res_locks <- rcpp_add_linking_z_le_w(op, x$data$dist_features)
    x$data$model_registry$cons$z_le_w <- res_locks
  }

  if (exists("rcpp_add_pu_locks", mode = "function")) {
    res_locks <- rcpp_add_pu_locks(op, x$data$pu)
    x$data$model_registry$cons$pu_locks <- res_locks
  }

  if (!is.null(x$data$dist_actions_model) && nrow(x$data$dist_actions_model) > 0 &&
      exists("rcpp_add_action_locks", mode = "function")) {

    res_locks <- rcpp_add_action_locks(op, x$data$dist_actions_model)
    x$data$model_registry$cons$action_locks <- res_locks
  }

  x$data$model_ptr   <- op
  x$data$model_index <- idx

  x
}


.pa_build_model_set_objective_cpp <- function(x) {

  .pa_abort <- function(...) stop(paste0(...), call. = FALSE)

  .pa_prepare_relation_model <- function(rel) {
    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)
    rel <- rel[order(rel$internal_pu1, rel$internal_pu2), , drop = FALSE]
    rel$internal_edge <- seq_len(nrow(rel))
    rel
  }

  op <- x$data$model_ptr
  if (is.null(op)) .pa_abort("model_ptr is NULL. Build core model first.")

  args  <- x$data$model_args %||% list()
  mtype <- args$model_type %||% "minimizeCosts"
  oargs <- args$objective_args %||% list()

  # ------------------------------------------------------------
  # Decide modelsense for the atomic objective we are activating
  # ------------------------------------------------------------
  modelsense <- if (mtype %in% c("minimizeCosts",
                                 "minimizeLosses",
                                 "minimizeFragmentation",
                                 "minimizeActionFragmentation",
                                 "minimizeInterventionFragmentation",
                                 "minimizeInterventionImpact")) "min" else "max"

  if (!exists("rcpp_reset_objective", mode = "function")) {
    .pa_abort("Missing rcpp_reset_objective() in the package.")
  }
  if (identical(mtype, "maximizeBenefits")) modelsense <- oargs$effect_sense %||% "max"
  rcpp_reset_objective(op, modelsense)


  # ------------------------------------------------------------
  # ONE atomic objective (prepare + add)
  # ------------------------------------------------------------
  if (identical(mtype, "minimizeCosts")) {

    if (!exists("rcpp_prepare_objective_min_cost", mode = "function")) .pa_abort("Missing rcpp_prepare_objective_min_cost().")
    if (!exists("rcpp_add_objective_min_cost",     mode = "function")) .pa_abort("Missing rcpp_add_objective_min_cost().")

    # Keep native row indices: only the cost contribution is scoped, never
    # the feasible decisions or the global planning-unit cost component.
    da_cost <- x$data$dist_actions_model
    if (!is.null(oargs$actions)) {
      action_ids <- .pa_resolve_action_subset(x, oargs$actions)$internal_id
      da_cost <- da_cost[da_cost$internal_action %in% action_ids, , drop = FALSE]
    }

    rcpp_prepare_objective_min_cost(
      op,
      pu_data = x$data$pu,
      dist_actions_data = da_cost,
      include_pu_cost = isTRUE(oargs$include_pu_cost %||% TRUE),
      include_action_cost = isTRUE(oargs$include_action_cost %||% TRUE),
      block_name = "objective_min_cost",
      tag = as.character(oargs$tag %||% "")[1]
    )

    res <- rcpp_add_objective_min_cost(
      op,
      pu_data = x$data$pu,
      dist_actions_data = da_cost,
      include_pu_cost = isTRUE(oargs$include_pu_cost %||% TRUE),
      include_action_cost = isTRUE(oargs$include_action_cost %||% TRUE),
      weight = 1.0
    )

    objective_id <- "min_cost"

  } else if (identical(mtype, "maximizeBenefits")) {

    if (.pa_uses_aggregate_effects(x)) {
      vector <- .pa_joint_objective_vector(x, oargs$benefit_col %||% "effect",
                                           oargs$actions, oargs$features)
      rcpp_model_set_objective_vector(op, vector, modelsense)
      x$data$model_args$modelsense <- modelsense
      x$data$model_args$objective_id <- if (modelsense == "min") "min_effect" else (x$data$model_args$objective_id %||% "max_effect")
      x$data$model_registry$objective <- list(type = mtype, id = x$data$model_args$objective_id, joint_effects = TRUE)
      return(x)
    }

    if (!exists("rcpp_prepare_objective_max_benefit", mode = "function")) {
      .pa_abort("Missing rcpp_prepare_objective_max_benefit().")
    }
    if (!exists("rcpp_add_objective_max_benefit", mode = "function")) {
      .pa_abort("Missing rcpp_add_objective_max_benefit().")
    }

    de <- x$data$dist_effects_model
    if (is.null(de) || !inherits(de, "data.frame")) {
      .pa_abort("Missing x$data$dist_effects_model for objective 'maximizeBenefits'.")
    }

    bcol <- as.character(oargs$benefit_col %||% "effect")[1]

    if (is.na(bcol) || !nzchar(bcol)) {
      .pa_abort("`benefit_col` must be a non-empty string.")
    }

    for (nm in c("internal_pu", "internal_action", "internal_feature")) {
      if (!(nm %in% names(de))) {
        .pa_abort("dist_effects_model must contain column '", nm, "'.")
      }
    }

    if (!(bcol %in% names(de))) {
      .pa_abort(
        "dist_effects_model must contain column '", bcol,
        "' for objective 'maximizeBenefits'."
      )
    }

    # filter action subset if requested
    acts <- oargs$actions %||% NULL
    if (!is.null(acts)) {
      de <- de[de$internal_action %in% as.integer(acts), , drop = FALSE]
    }

    # filter feature subset if requested
    feats <- oargs$features %||% NULL
    if (!is.null(feats)) {
      de <- de[de$internal_feature %in% as.integer(feats), , drop = FALSE]
    }

    # C++ expects a column named 'benefit'; supply signed effects in a local
    # copy while preserving the canonical positive/loss reporting components.
    de$benefit <- as.numeric(de[[bcol]])

    if (anyNA(de$benefit) || any(!is.finite(de$benefit))) {
      .pa_abort(
        "Objective 'maximizeBenefits' uses column '", bcol,
        "', but it contains NA or non-finite values."
      )
    }

    prep <- rcpp_prepare_objective_max_benefit(
      op,
      dist_actions_data = x$data$dist_actions_model,
      dist_effects_data = de,
      block_name = "objective_max_benefit",
      tag = as.character(oargs$tag %||% "")[1]
    )

    coef_x <- as.numeric(prep$coef_x)

    res <- rcpp_add_objective_max_benefit(
      op,
      coef_x = coef_x,
      weight = 1.0,
      block_name = "objective_max_benefit",
      tag = as.character(oargs$tag %||% "")[1]
    )

    rcpp_model_set_objective_vector(op, rcpp_optimization_problem_as_list(op)$obj, modelsense)
    objective_id <- if (modelsense == "min") "min_effect" else (x$data$model_args$objective_id %||% "max_effect")

  } else if (identical(mtype, "minimizeLosses")) {

    if (.pa_uses_aggregate_effects(x)) {
      vector <- .pa_joint_objective_vector(x, "loss", oargs$actions, oargs$features)
      rcpp_model_set_objective_vector(op, vector, "min")
      x$data$model_args$modelsense <- "min"
      x$data$model_args$objective_id <- "min_loss"
      x$data$model_registry$objective <- list(type = mtype, id = "min_loss", joint_effects = TRUE)
      return(x)
    }

    if (!exists("rcpp_prepare_objective_min_loss", mode = "function")) {
      .pa_abort("Missing rcpp_prepare_objective_min_loss().")
    }
    if (!exists("rcpp_add_objective_min_loss", mode = "function")) {
      .pa_abort("Missing rcpp_add_objective_min_loss().")
    }

    de <- x$data$dist_effects_model
    if (is.null(de) || !inherits(de, "data.frame")) {
      .pa_abort("Missing x$data$dist_effects_model for objective 'minimizeLosses'.")
    }

    for (nm in c("internal_pu", "internal_action", "internal_feature", "loss")) {
      if (!(nm %in% names(de))) {
        .pa_abort("dist_effects_model must contain column '", nm, "'.")
      }
    }

    # filter action subset if requested
    acts <- oargs$actions %||% NULL
    if (!is.null(acts)) {
      de <- de[de$internal_action %in% as.integer(acts), , drop = FALSE]
    }

    # filter feature subset if requested
    feats <- oargs$features %||% NULL
    if (!is.null(feats)) {
      de <- de[de$internal_feature %in% as.integer(feats), , drop = FALSE]
    }

    prep <- rcpp_prepare_objective_min_loss(
      x = op,
      dist_actions_data = x$data$dist_actions_model,
      dist_effects_data = de,
      block_name = "objective_min_loss",
      tag = as.character(oargs$tag %||% "")[1]
    )

    res <- rcpp_add_objective_min_loss(
      x = op,
      coef_x = as.numeric(prep$coef_x),
      weight = 1.0,
      block_name = "objective_min_loss",
      tag = as.character(oargs$tag %||% "")[1]
    )

    objective_id <- "min_loss"

  } else if (identical(mtype, "maximizeProfit")) {

    if (!exists("rcpp_prepare_objective_max_profit", mode = "function")) .pa_abort("Missing rcpp_prepare_objective_max_profit().")
    if (!exists("rcpp_add_objective_max_profit",     mode = "function")) .pa_abort("Missing rcpp_add_objective_max_profit().")

    pcol <- as.character(oargs$profit_col %||% "profit")[1]

    rcpp_prepare_objective_max_profit(
      op,
      dist_actions_data = x$data$dist_actions_model,
      dist_profit_data = x$data$dist_profit_model,
      profit_col = pcol,
      block_name = "objective_max_profit",
      tag = as.character(oargs$tag %||% "")[1]
    )

    res <- rcpp_add_objective_max_profit(
      op,
      dist_actions_data = x$data$dist_actions_model,
      dist_profit_data = x$data$dist_profit_model,
      profit_col = pcol,
      weight = 1.0
    )

    objective_id <- "max_profit"

  } else if (identical(mtype, "maximizeNetProfit")) {

    if (!exists("rcpp_prepare_objective_max_net_profit", mode = "function")) .pa_abort("Missing rcpp_prepare_objective_max_net_profit().")
    if (!exists("rcpp_add_objective_max_net_profit",     mode = "function")) .pa_abort("Missing rcpp_add_objective_max_net_profit().")

    pcol <- as.character(oargs$profit_col %||% "profit")[1]

    rcpp_prepare_objective_max_net_profit(
      op,
      pu_data = x$data$pu,
      dist_actions_data = x$data$dist_actions_model,
      dist_profit_data = x$data$dist_profit_model,
      profit_col = pcol,
      include_pu_cost = isTRUE(oargs$include_pu_cost %||% TRUE),
      include_action_cost = isTRUE(oargs$include_action_cost %||% TRUE),
      block_name = "objective_max_net_profit",
      tag = as.character(oargs$tag %||% "")[1]
    )

    res <- rcpp_add_objective_max_net_profit(
      op,
      pu_data = x$data$pu,
      dist_actions_data = x$data$dist_actions_model,
      dist_profit_data = x$data$dist_profit_model,
      profit_col = pcol,
      include_pu_cost = isTRUE(oargs$include_pu_cost %||% TRUE),
      include_action_cost = isTRUE(oargs$include_action_cost %||% TRUE),
      weight = 1.0
    )

    objective_id <- "max_net_profit"


  } else if (identical(mtype, "minimizeActionFragmentation")) {

    if (!exists("rcpp_add_objective_min_fragmentation_actions", mode = "function")) {
      .pa_abort("Missing rcpp_add_objective_min_fragmentation_actions().")
    }

    rel_name  <- as.character(oargs$relation_name %||% "boundary")[1]
    rel       <- x$data$spatial_relations[[rel_name]]
    if (is.null(rel)) .pa_abort("Missing spatial relation: ", rel_name)

    rel_model <- .pa_prepare_relation_model(rel)
    x$data$spatial_relations_model <- x$data$spatial_relations_model %||% list()
    x$data$spatial_relations_model[[rel_name]] <- rel_model

    subset_actions_raw <- oargs$actions_to_use %||% oargs$actions %||% NULL

    act_subset <- NULL
    actions_to_use <- NULL
    if (!is.null(subset_actions_raw)) {
      act_subset <- .pa_resolve_action_subset(x, subset = subset_actions_raw)
      actions_to_use <- as.integer(act_subset$internal_id)
    }

    aw_vec <- .pa_action_weights_vector(
      actions_df = x$data$actions,
      action_weights = oargs$action_weights %||% NULL,
      subset_actions = actions_to_use,
      default_weight = 1
    )

    res <- rcpp_add_objective_min_fragmentation_actions(
      op,
      dist_actions_data = x$data$dist_actions_model,
      relation_data     = rel_model,
      actions_to_use    = actions_to_use,
      action_weights    = aw_vec,
      weight_multiplier = as.numeric(oargs$weight_multiplier %||% 1)[1]
    )

    objective_id <- "min_action_fragmentation"

  } else if (identical(mtype, "minimizeFragmentation")) {

    if (!exists("rcpp_add_objective_min_fragmentation", mode = "function")) .pa_abort("Missing rcpp_add_objective_min_fragmentation().")

    rel_name  <- as.character(oargs$relation_name %||% "boundary")[1]
    rel       <- x$data$spatial_relations[[rel_name]]
    if (is.null(rel)) .pa_abort("Missing spatial relation: ", rel_name)

    rel_model <- x$data$spatial_relations_model[[rel_name]] %||% .pa_prepare_relation_model(rel)
    x$data$spatial_relations_model <- x$data$spatial_relations_model %||% list()
    x$data$spatial_relations_model[[rel_name]] <- rel_model

    # IMPORTANT: prepare of aux vars (idempotent)
    if (exists("rcpp_prepare_objective_min_fragmentation", mode = "function")) {
      rcpp_prepare_objective_min_fragmentation(op, rel_model)
    }

    res <- rcpp_add_objective_min_fragmentation(
      op,
      relation_data = rel_model,
      weight = 1.0 ,
      weight_multiplier = as.numeric(oargs$weight_multiplier %||% 1)[1]
    )

    objective_id <- "min_fragmentation"

  } else if (identical(mtype, "minimizeInterventionImpact")) {

    if (!exists("rcpp_prepare_objective_min_intervention_impact", mode = "function")) {
      .pa_abort("Missing rcpp_prepare_objective_min_intervention_impact().")
    }
    if (!exists("rcpp_add_objective_min_intervention_impact", mode = "function")) {
      .pa_abort("Missing rcpp_add_objective_min_intervention_impact().")
    }

    icol <- as.character(oargs$impact_col %||% "amount")[1]

    # ---- feature subset
    feats <- oargs$features %||% NULL
    feats_internal <- integer(0)
    if (!is.null(feats)) {
      m <- match(feats, x$data$features$id)
      if (anyNA(m)) {
        .pa_abort(
          "minimizeInterventionImpact: some features not in x$data$features$id: ",
          paste(feats[is.na(m)], collapse = ", ")
        )
      }
      feats_internal <- as.integer(x$data$features$internal_id[m])
    }

    # ---- action subset (REQUIRED for option 2)
    acts <- oargs$actions %||% NULL
    if (is.null(acts)) {
      .pa_abort(
        "Objective 'minimizeInterventionImpact' now requires `actions=`.\n",
        "This objective is defined over intervention-group variables linked to a subset of actions."
      )
    }

    act_subset <- .pa_resolve_action_subset(x, subset = acts)
    acts_internal <- as.integer(act_subset$internal_id)

    if (length(acts_internal) == 0) {
      .pa_abort("minimizeInterventionImpact: action subset resolved to zero actions.")
    }

    # optional identifier for the action subset/group
    subset_key <- .pa_subset_to_string(act_subset$id)

    prep <- rcpp_prepare_objective_min_intervention_impact(
      op,
      pu_data = x$data$pu,
      dist_actions_data = x$data$dist_actions_model,
      dist_features_data = x$data$dist_features,
      actions_to_use = acts_internal,
      subset_key = subset_key,
      impact_col = icol,
      features_to_use = feats_internal,
      internal_feature_col = "internal_feature",
      block_name = "objective_min_intervention_impact",
      tag = as.character(oargs$tag %||% "")[1]
    )

    res <- rcpp_add_objective_min_intervention_impact(
      op,
      pu_data = x$data$pu,
      dist_actions_data = x$data$dist_actions_model,
      dist_features_data = x$data$dist_features,
      actions_to_use = acts_internal,
      subset_key = subset_key,
      impact_col = icol,
      features_to_use = feats_internal,
      internal_feature_col = "internal_feature",
      weight = 1.0
    )

    x$data$model_registry$vars$u_intervention_impact <- prep
    objective_id <- "min_intervention_impact"
  } else {
    .pa_abort("Unknown model_type: ", mtype)
  }

  x$data$model_args$modelsense   <- modelsense
  x$data$model_args$objective_id <- objective_id

  x$data$model_registry$objective <- list(type = mtype, id = objective_id)
  x$data$model_registry$objective$cpp <- res

  x
}

.pa_build_model_apply_constraints <- function(x) {

  if (!is.null(x$data$targets) &&
      inherits(x$data$targets, "data.frame") &&
      nrow(x$data$targets) > 0) {

    if (!exists(".pa_apply_targets_if_present", mode = "function")) {
      stop(".pa_apply_targets_if_present() is missing.", call. = FALSE)
    }

    x <- .pa_apply_targets_if_present(x)
  }

  x <- .pa_apply_action_max_per_pu_default(x)
  x <- .pa_apply_action_cardinality_if_present(x)
  x <- .pa_apply_action_relations_if_present(x)

  if (exists(".pa_apply_area_constraints_if_present", mode = "function")) {
    x <- .pa_apply_area_constraints_if_present(x)
  }

  if (exists(".pa_apply_budget_constraints_if_present", mode = "function")) {
    x <- .pa_apply_budget_constraints_if_present(x)
  }

  if (exists(".pa_apply_group_area_constraints_if_present", mode = "function")) {
    x <- .pa_apply_group_area_constraints_if_present(x)
  }

  x
}


.pa_build_model_finalize <- function(x) {

  if (!exists(".pa_refresh_model_snapshot", mode = "function")) {
    stop(".pa_refresh_model_snapshot() is missing.", call. = FALSE)
  }

  x <- .pa_refresh_model_snapshot(x)
  x
}


.pa_apply_action_max_per_pu_default <- function(x) {

  stopifnot(inherits(x, "Problem"))

  da <- x$data$dist_actions_model
  if (is.null(da) || !inherits(da, "data.frame") || nrow(da) == 0) {
    return(x)
  }

  if (!exists("rcpp_add_action_max_per_pu", mode = "function")) {
    stop("Missing rcpp_add_action_max_per_pu().", call. = FALSE)
  }

  # Only explicit TOTAL maxima or equalities replace the implicit maximum.
  # Do not pass an empty filtered vector to C++: there it means all units.
  specs <- x$data$constraints$action_cardinality
  covered <- .pa_action_cardinality_covered_pu(specs)
  default_pu <- integer()
  if (length(covered)) {
    default_pu <- unique(as.integer(da$internal_pu[!da$pu %in% covered]))
    if (!length(default_pu)) return(x)
  }

  res <- rcpp_add_action_max_per_pu(
    x$data$model_ptr,
    dist_actions_data = da,
    max_per_pu = 1L,
    internal_pu_ids = default_pu,
    internal_action_ids = integer()
  )

  x$data$model_registry$cons$action_max_per_pu <- res
  x
}


#' Ensure implicit conservation model when actions/effects are omitted
#'
#' @noRd
.pa_build_model_ensure_implicit_conservation <- function(x) {
  stopifnot(inherits(x, "Problem"))

  if (.pa_needs_implicit_conservation_model(x)) {
    x <- .pa_add_implicit_conservation_model(x)
  }

  x
}


.pa_needs_implicit_conservation_model <- function(x) {
  stopifnot(inherits(x, "Problem"))

  has_rows <- function(df) {
    !is.null(df) && inherits(df, "data.frame") && nrow(df) > 0
  }

  has_actions <- has_rows(x$data$actions) || has_rows(x$data$dist_actions)
  has_effects <- has_rows(x$data$dist_effects)

  !has_actions && !has_effects
}


.pa_add_implicit_conservation_model <- function(x,
                                                action_id = "conservation",
                                                action_name = "Conservation") {
  stopifnot(inherits(x, "Problem"))

  has_rows <- function(df) {
    !is.null(df) && inherits(df, "data.frame") && nrow(df) > 0
  }

  has_actions <- has_rows(x$data$actions) || has_rows(x$data$dist_actions)
  has_effects <- has_rows(x$data$dist_effects)

  if (has_actions || has_effects) {
    return(x)
  }

  if (is.null(x$data$pu) ||
      is.null(x$data$features) ||
      is.null(x$data$dist_features)) {
    stop(
      "Implicit conservation mode requires planning units, features, and dist_features.",
      call. = FALSE
    )
  }

  if (!all(c("id", "internal_id") %in% names(x$data$pu))) {
    stop(
      "Implicit conservation mode requires x$data$pu to contain 'id' and 'internal_id'.",
      call. = FALSE
    )
  }

  if (!all(c("id", "internal_id") %in% names(x$data$features))) {
    stop(
      "Implicit conservation mode requires x$data$features to contain 'id' and 'internal_id'.",
      call. = FALSE
    )
  }

  if (!all(c("pu", "feature", "amount") %in% names(x$data$dist_features))) {
    stop(
      "Implicit conservation mode requires dist_features columns 'pu', 'feature', and 'amount'.",
      call. = FALSE
    )
  }

  actions <- data.frame(
    id = action_id,
    name = action_name,
    action_set = "conservation",
    stringsAsFactors = FALSE
  )

  x <- add_actions(
    x = x,
    actions = actions,
    include_pairs = NULL,
    exclude_pairs = NULL,
    cost = 0
  )

  effects <- data.frame(
    pu = x$data$dist_features$pu,
    action = action_id,
    feature = x$data$dist_features$feature,
    effect = 0,
    stringsAsFactors = FALSE
  )

  x <- add_effects(
    x = x,
    effects = effects
  )

  x$data$meta <- x$data$meta %||% list()
  x$data$meta$model_mode <- "simple_conservation"
  x$data$meta$implicit_actions <- TRUE
  x$data$meta$implicit_effects <- TRUE
  x$data$meta$implicit_action_id <- action_id
  x$data$meta$model_dirty <- TRUE

  x$data$effects_meta <- x$data$effects_meta %||% list()
  x$data$effects_meta$input <- "implicit_conservation"
  x$data$effects_meta$effect_type <- "after"
  x$data$effects_meta$component <- "amount_after"
  x$data$effects_meta$action_id <- action_id
  x$data$effects_meta$action_name <- action_name

  x
}