From 3fb77b22b609c7aa99f865b4dd5ddd963ce0c2dd Mon Sep 17 00:00:00 2001
From: zzylol <50204836+zzylol@users.noreply.github.com>
Date: Sat, 31 Oct 2026 17:38:13 +0000
Subject: [PATCH 1/4] =?UTF-8?q?feat(plan-selection):=20solve=20the=20selec?=
 =?UTF-8?q?tion=20MILP's=20independent=20parts=20apart=20(#1095=20=C2=A75,?=
 =?UTF-8?q?=203a)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The model splits into the connected components of the graph joining
each row's variables; each part is solved alone with its share of the
deadline, and the joined solution goes through the one completion,
decode, realization and plan_cost check.
---
 crates/plan-selection/src/milp/model.rs       | 463 +++++++++++++++---
 crates/plan-selection/src/milp/model/tests.rs |  89 +++-
 crates/plan-selection/src/select.rs           |   6 +-
 3 files changed, 469 insertions(+), 89 deletions(-)

diff --git a/crates/plan-selection/src/milp/model.rs b/crates/plan-selection/src/milp/model.rs
index f148955be..2c2755795 100644
--- a/crates/plan-selection/src/milp/model.rs
+++ b/crates/plan-selection/src/milp/model.rs
@@ -80,6 +80,9 @@ pub(crate) struct SelectionMilpOptions {
     /// Write the accuracy rows (row 8); off when the facade ignores
     /// accuracy guarantees.
     pub(crate) enforce_accuracy: bool,
+    /// Solve each independent part alone ([`independent_parts`]); off only
+    /// to compare with one solve in tests.
+    pub(crate) split_into_parts: bool,
 }
 
 /// Which part of the model a variable or row belongs to, for the size
@@ -262,6 +265,8 @@ pub(crate) struct SelectionMilp<'a> {
     /// producer in the other slot is not justified by `r`.
     pub(crate) slot_links:
         BTreeMap<(usize, PhysicalCandidateId, usize), Vec<(ResultSlot, Variable)>>,
+    /// [`SelectionMilpOptions::split_into_parts`].
+    pub(crate) split_into_parts: bool,
 }
 
 /// Variables and rows of one model.
@@ -304,6 +309,8 @@ pub(crate) struct SelectionMilpSolution {
     /// slot it reads.
     pub(crate) slot_reads: BTreeMap<(usize, PhysicalCandidateId, usize), ResultSlot>,
     pub(crate) size: SelectionMilpSize,
+    /// Per independent part, how its solve went, in part order.
+    pub(crate) parts: Vec<PartSolve>,
     pub(crate) build_ms: u64,
     /// Building the warm start, HiGHS's first incumbent.
     pub(crate) warm_start_ms: u64,
@@ -489,6 +496,7 @@ pub(crate) fn build_selection_milp<'a>(
         producer_gives: HashMap::new(),
         window_cache: HashMap::new(),
         slot_links: BTreeMap::new(),
+        split_into_parts: options.split_into_parts,
     };
     choice_variables(&mut milp);
     structural_rows(&mut milp);
@@ -1654,8 +1662,154 @@ fn accuracy_rows(milp: &mut SelectionMilp<'_>) {
     }
 }
 
-/// Solves `milp` once with HiGHS, with the time left before `deadline` as
-/// its time limit, from the warm start.
+/// One independent part of the model (#1095 §5, 3a): a connected component
+/// of the graph joining the variables of each row, or all components
+/// without a query together. Parts share no variable and no row, and the
+/// objective is a sum over variables, so solving each alone and joining
+/// their solutions solves the whole model.
+#[derive(Debug, Clone, PartialEq, Eq)]
+pub(crate) struct IndependentPart {
+    /// The queries whose output references it covers, ascending.
+    pub(crate) queries: Vec<usize>,
+    /// Its columns (variables in creation order), ascending.
+    pub(crate) columns: Vec<usize>,
+    /// Its rows, ascending.
+    pub(crate) rows: Vec<usize>,
+}
+
+/// The independent parts of `milp`, in order of their first query (the
+/// part without a query last). Queries fall into one part when they may
+/// share a shareable candidate, a run of `I`, a source's raw retention or
+/// anything else a row joins.
+pub(crate) fn independent_parts(milp: &SelectionMilp<'_>) -> Vec<IndependentPart> {
+    let column_of: HashMap<Variable, usize> = milp
+        .builder
+        .variables
+        .iter_variables_with_def()
+        .enumerate()
+        .map(|(column, (variable, _))| (variable, column))
+        .collect();
+    let count = column_of.len();
+    // Union-find over the columns, joined along each row.
+    let mut parent: Vec<usize> = (0..count).collect();
+    fn root(parent: &mut [usize], mut node: usize) -> usize {
+        while parent[node] != node {
+            parent[node] = parent[parent[node]];
+            node = parent[node];
+        }
+        node
+    }
+    for row in &milp.builder.rows {
+        let mut first: Option<usize> = None;
+        for (variable, _) in row.expression().linear_coefficients() {
+            let column = root(&mut parent, column_of[&variable]);
+            match first {
+                None => first = Some(column),
+                Some(known) if known != column => parent[column] = known,
+                Some(_) => {}
+            }
+        }
+    }
+    let output_refs = output_references(&milp.physical.query_roots);
+    let mut query_of_component: BTreeMap<usize, usize> = BTreeMap::new();
+    for (reference, &(query, _)) in output_refs.iter().enumerate() {
+        for out in &milp.inputs.output_candidates[reference] {
+            if let Some(variable) = milp.choice[query].get(out) {
+                let component = root(&mut parent, column_of[variable]);
+                query_of_component.entry(component).or_insert(query);
+            }
+        }
+    }
+    // The part of each component: its first query's, or the last.
+    let part_of = |component: usize| query_of_component.get(&component).copied();
+    let mut parts: BTreeMap<Option<usize>, IndependentPart> = BTreeMap::new();
+    let mut part_of_column = Vec::with_capacity(count);
+    for column in 0..count {
+        let key = part_of(root(&mut parent, column));
+        part_of_column.push(key);
+        parts
+            .entry(key)
+            .or_insert_with(|| IndependentPart {
+                queries: Vec::new(),
+                columns: Vec::new(),
+                rows: Vec::new(),
+            })
+            .columns
+            .push(column);
+    }
+    for (reference, &(query, _)) in output_refs.iter().enumerate() {
+        let Some(variable) = milp.inputs.output_candidates[reference]
+            .iter()
+            .find_map(|out| milp.choice[query].get(out))
+        else {
+            continue;
+        };
+        let key = part_of_column[column_of[variable]];
+        let queries = &mut parts.get_mut(&key).expect("its column's part").queries;
+        if !queries.contains(&query) {
+            queries.push(query);
+        }
+    }
+    for (index, row) in milp.builder.rows.iter().enumerate() {
+        // A row without a variable is a constant check: kept in the first
+        // part.
+        let key = row
+            .expression()
+            .linear_coefficients()
+            .next()
+            .map(|(variable, _)| part_of_column[column_of[&variable]])
+            .unwrap_or_else(|| *parts.keys().next().expect("a part"));
+        parts.get_mut(&key).expect("a part").rows.push(index);
+    }
+    // `None` sorts first; the part without a query goes last.
+    let mut parts: Vec<(Option<usize>, IndependentPart)> = parts.into_iter().collect();
+    parts.sort_by_key(|(key, _)| key.unwrap_or(usize::MAX));
+    parts
+        .into_iter()
+        .map(|(_, mut part)| {
+            part.queries.sort_unstable();
+            part
+        })
+        .collect()
+}
+
+/// How one part's solve went, for the report.
+#[derive(Debug, Clone, PartialEq, serde::Serialize)]
+pub(crate) struct PartSolve {
+    pub(crate) queries: usize,
+    pub(crate) variables: usize,
+    pub(crate) binaries: usize,
+    pub(crate) rows: usize,
+    pub(crate) status: String,
+    pub(crate) solve_ms: u64,
+    /// HiGHS's dual bound on the part's objective (cost per second).
+    pub(crate) dual_bound_per_second: f64,
+}
+
+/// One part's model ready for HiGHS: its variables, rows and objective
+/// terms, renamed to its own columns when the model has several parts.
+struct PartModel {
+    variables: ProblemVariables,
+    rows: Vec<Constraint>,
+    objective_terms: Vec<(Variable, f64)>,
+    /// Per column of the part, the column in the whole model.
+    columns: Vec<usize>,
+    queries: usize,
+}
+
+/// A part solved: its columns' values, the HiGHS objective and dual bound
+/// (cost per second), and the status.
+struct SolvedPart {
+    status: SelectionMilpStatus,
+    columns: Vec<f64>,
+    objective_per_second: f64,
+    dual_bound_per_second: f64,
+    report: PartSolve,
+}
+
+/// Solves `milp` with HiGHS, one solve per independent part
+/// ([`independent_parts`]) when `split_into_parts`, each part from the warm
+/// start with its share of the time left before `deadline` (by variables).
 pub(crate) fn solve_selection_milp(
     milp: SelectionMilp<'_>,
     deadline: Option<Instant>,
@@ -1672,21 +1826,18 @@ pub(crate) fn solve_selection_milp(
     let start = warm_start::warm_start(&milp, &completion, Some(search_budget));
     let warm_start_ms = warm_started.elapsed().as_millis() as u64;
     let warm_start_cost_per_second = start.as_ref().map(|start| start.cost_per_second);
-    let scale = objective_scale(
-        &milp.objective_terms,
-        start.as_ref().map(|start| start.values.as_slice()),
-    );
-    let objective: Expression = milp
-        .objective_terms
-        .iter()
-        .map(|&(variable, coefficient)| coefficient * scale * variable)
-        .sum();
+    let parts = if milp.split_into_parts {
+        independent_parts(&milp)
+    } else {
+        Vec::new()
+    };
     let SelectionMilp {
         builder,
         choice,
         selected,
         slot_links,
         request_cost_per_second,
+        objective_terms,
         ..
     } = milp;
     let ModelBuilder {
@@ -1697,100 +1848,150 @@ pub(crate) fn solve_selection_milp(
         .enumerate()
         .map(|(column, (variable, _))| (variable, column))
         .collect();
-    let problem = variables
-        .minimise(objective.clone())
-        .using(highs)
-        .with_all(rows);
-    let mut model = problem
-        .try_into_inner()
-        .map_err(|error| SelectionMilpError::Solver(error.to_string()))?;
-    model.set_option("mip_rel_gap", 0.0);
-    model.set_option("mip_abs_gap", 0.0);
-    model.set_option("output_flag", false);
-    // HiGHS's MIP search is mostly serial; a thread per core costs
-    // gigabytes of address space on a large host (#1036 under 8 GB).
-    model.set_option("threads", HIGHS_THREADS);
-    let started_with_plan = start.is_some();
-    if let Some(start) = start {
-        let mut columns = vec![0.0; column_of.len()];
+    let column_count = column_of.len();
+    let start_columns: Option<Vec<f64>> = start.map(|start| {
+        let mut columns = vec![0.0; column_count];
         for (variable, value) in start.values {
             columns[column_of[&variable]] = value;
         }
-        model
-            .try_set_solution(Some(&columns), None, None, None)
-            .map_err(|status| SelectionMilpError::Solver(format!("warm start: {status:?}")))?;
-    }
-    let mut only_read_back = false;
-    if let Some(deadline) = deadline {
-        let seconds = deadline
-            .saturating_duration_since(Instant::now())
-            .as_secs_f64();
-        if seconds <= 0.0 {
-            if !started_with_plan {
-                return Err(SelectionMilpError::DeadlineReached);
+        columns
+    });
+    let models: Vec<PartModel> = if parts.len() <= 1 {
+        vec![PartModel {
+            variables,
+            rows,
+            objective_terms,
+            columns: (0..column_count).collect(),
+            queries: parts.first().map_or(0, |part| part.queries.len()),
+        }]
+    } else {
+        let definitions: Vec<_> = variables
+            .iter_variables_with_def()
+            .map(|(_, definition)| definition.clone())
+            .collect();
+        let mut part_of_column = vec![0; column_count];
+        let mut renamed: Vec<Option<Variable>> = vec![None; column_count];
+        let mut models: Vec<PartModel> = parts
+            .iter()
+            .enumerate()
+            .map(|(index, part)| {
+                let mut variables = ProblemVariables::new();
+                for &column in &part.columns {
+                    part_of_column[column] = index;
+                    renamed[column] = Some(variables.add(definitions[column].clone()));
+                }
+                PartModel {
+                    variables,
+                    rows: Vec::with_capacity(part.rows.len()),
+                    objective_terms: Vec::new(),
+                    columns: part.columns.clone(),
+                    queries: part.queries.len(),
+                }
+            })
+            .collect();
+        let rename = |variable: Variable| renamed[column_of[&variable]].expect("in a part");
+        for (index, part) in parts.iter().enumerate() {
+            for &row in &part.rows {
+                let row = &rows[row];
+                let mut expression = Expression::from(row.expression().constant());
+                for (variable, coefficient) in row.expression().linear_coefficients() {
+                    expression.add_mul(coefficient, rename(variable));
+                }
+                models[index].rows.push(if row.is_equality() {
+                    constraint!(expression == 0.0)
+                } else {
+                    constraint!(expression <= 0.0)
+                });
             }
-            only_read_back = true;
         }
-        model.set_option("time_limit", seconds.max(MOMENT_TO_READ_THE_START_S));
-    }
-    let started = Instant::now();
-    let solved = model
-        .try_solve()
-        .map_err(|status| SelectionMilpError::Solver(format!("{status:?}")))?;
-    let solve_ms = started.elapsed().as_millis() as u64;
-    use highs::{HighsModelStatus, HighsSolutionStatus};
-    let status = match solved.status() {
-        HighsModelStatus::Optimal if only_read_back => SelectionMilpStatus::TimeLimit,
-        HighsModelStatus::Optimal => SelectionMilpStatus::Optimal,
-        HighsModelStatus::Infeasible | HighsModelStatus::UnboundedOrInfeasible => {
-            return Err(SelectionMilpError::Infeasible)
+        for (variable, coefficient) in objective_terms {
+            let column = column_of[&variable];
+            models[part_of_column[column]]
+                .objective_terms
+                .push((rename(variable), coefficient));
         }
-        HighsModelStatus::ReachedTimeLimit
-        | HighsModelStatus::ReachedInterrupt
-        | HighsModelStatus::ReachedIterationLimit
-        | HighsModelStatus::ReachedMemoryLimit
-        | HighsModelStatus::ReachedSolutionLimit => SelectionMilpStatus::TimeLimit,
-        other => return Err(SelectionMilpError::Solver(format!("{other:?}"))),
+        models
     };
-    if solved.primal_solution_status() != HighsSolutionStatus::Feasible {
-        return Err(SelectionMilpError::NoSolution {
-            status: format!("{:?}", solved.status()),
+    // Smallest first: time a small part leaves goes to the larger ones.
+    let mut order: Vec<usize> = (0..models.len()).collect();
+    order.sort_by_key(|&index| models[index].columns.len());
+    let mut remaining_columns: usize = models.iter().map(|model| model.columns.len()).sum();
+    let mut models: Vec<Option<PartModel>> = models.into_iter().map(Some).collect();
+    let mut columns = vec![0.0; column_count];
+    let mut solved_parts: Vec<Option<SolvedPart>> = (0..models.len()).map(|_| None).collect();
+    let started = Instant::now();
+    for index in order {
+        let model = models[index].take().expect("solved once");
+        let share = model.columns.len().max(1) as f64 / remaining_columns.max(1) as f64;
+        remaining_columns -= model.columns.len();
+        let part_deadline = deadline.map(|deadline| {
+            Instant::now()
+                + deadline
+                    .saturating_duration_since(Instant::now())
+                    .mul_f64(share)
+        });
+        let start = start_columns.as_ref().map(|start| {
+            model
+                .columns
+                .iter()
+                .map(|&column| start[column])
+                .collect::<Vec<f64>>()
         });
+        let part_columns = model.columns.clone();
+        let solved = solve_part(model, start, part_deadline)?;
+        for (&column, &value) in part_columns.iter().zip(&solved.columns) {
+            columns[column] = value;
+        }
+        solved_parts[index] = Some(solved);
     }
+    let solve_ms = started.elapsed().as_millis() as u64;
+    let solved_parts: Vec<SolvedPart> = solved_parts
+        .into_iter()
+        .map(|part| part.expect("every part solved"))
+        .collect();
+    let status = if solved_parts
+        .iter()
+        .all(|part| part.status == SelectionMilpStatus::Optimal)
+    {
+        SelectionMilpStatus::Optimal
+    } else {
+        SelectionMilpStatus::TimeLimit
+    };
+    let highs_objective: f64 = solved_parts
+        .iter()
+        .map(|part| part.objective_per_second)
+        .sum();
+    let dual_bound: f64 = solved_parts
+        .iter()
+        .map(|part| part.dual_bound_per_second)
+        .sum();
     let relative_gap = match status {
         SelectionMilpStatus::Optimal => 0.0,
-        SelectionMilpStatus::TimeLimit => solved.mip_gap(),
+        SelectionMilpStatus::TimeLimit => {
+            (highs_objective - dual_bound) / highs_objective.abs().max(f64::MIN_POSITIVE)
+        }
     };
-    let dual_bound = solved
-        .double_info_value(c"mip_dual_bound")
-        .unwrap_or(f64::NAN);
-    let solution = solved.get_solution();
-    let columns = solution.columns();
     // The objective of the integral plan: the incumbent's binaries,
     // rounded, with every continuous variable at its cheapest value (the
     // solver's own value carries its tolerances, relative to the scaled
     // objective). At the deadline, HiGHS's incumbent objective (an upper
     // bound) when the completion does not apply.
-    let objective_value = match completion.complete_columns(columns) {
-        Ok(completed) => completed.cost_per_second * scale,
+    let objective_value = match completion.complete_columns(&columns) {
+        Ok(completed) => completed.cost_per_second,
         Err(reason) => {
             warm_start::debug(|| format!("completion of the incumbent failed: {reason}"));
-            solved.objective_value()
+            highs_objective
         }
     };
     warm_start::debug(|| {
-        format!(
-            "objective: HiGHS {:e}, completed {:e}",
-            solved.objective_value() / scale,
-            objective_value / scale
-        )
+        format!("objective: HiGHS {highs_objective:e}, completed {objective_value:e}")
     });
     let on = |variable: Variable| columns[column_of[&variable]] > 0.5;
     Ok(SelectionMilpSolution {
         status,
-        cost_per_second: objective_value / scale,
+        cost_per_second: objective_value,
         request_cost_per_second,
-        cost_lower_bound_per_second: dual_bound / scale,
+        cost_lower_bound_per_second: dual_bound,
         relative_gap,
         chosen: choice
             .iter()
@@ -1817,6 +2018,7 @@ pub(crate) fn solve_selection_milp(
             })
             .collect(),
         size,
+        parts: solved_parts.into_iter().map(|part| part.report).collect(),
         build_ms,
         warm_start_ms,
         warm_start_cost_per_second,
@@ -1824,6 +2026,111 @@ pub(crate) fn solve_selection_milp(
     })
 }
 
+/// Solves one part once with HiGHS, from `start` (its columns' values) when
+/// given, with the time left before `deadline` as its time limit.
+fn solve_part(
+    model: PartModel,
+    start: Option<Vec<f64>>,
+    deadline: Option<Instant>,
+) -> Result<SolvedPart, SelectionMilpError> {
+    let PartModel {
+        variables,
+        rows,
+        objective_terms,
+        columns: part_columns,
+        queries,
+    } = model;
+    let start_values: Option<Vec<(Variable, f64)>> = start.as_ref().map(|start| {
+        variables
+            .iter_variables_with_def()
+            .map(|(variable, _)| variable)
+            .zip(start.iter().copied())
+            .collect()
+    });
+    let scale = objective_scale(&objective_terms, start_values.as_deref());
+    let objective: Expression = objective_terms
+        .iter()
+        .map(|&(variable, coefficient)| coefficient * scale * variable)
+        .sum();
+    let binaries = variables
+        .iter_variables_with_def()
+        .filter(|(_, definition)| definition.is_integer())
+        .count();
+    let row_count = rows.len();
+    let problem = variables.minimise(objective).using(highs).with_all(rows);
+    let mut model = problem
+        .try_into_inner()
+        .map_err(|error| SelectionMilpError::Solver(error.to_string()))?;
+    model.set_option("mip_rel_gap", 0.0);
+    model.set_option("mip_abs_gap", 0.0);
+    model.set_option("output_flag", false);
+    // HiGHS's MIP search is mostly serial; a thread per core costs
+    // gigabytes of address space on a large host (#1036 under 8 GB).
+    model.set_option("threads", HIGHS_THREADS);
+    let started_with_plan = start.is_some();
+    if let Some(start) = &start {
+        model
+            .try_set_solution(Some(start), None, None, None)
+            .map_err(|status| SelectionMilpError::Solver(format!("warm start: {status:?}")))?;
+    }
+    let mut only_read_back = false;
+    if let Some(deadline) = deadline {
+        let seconds = deadline
+            .saturating_duration_since(Instant::now())
+            .as_secs_f64();
+        if seconds <= 0.0 {
+            if !started_with_plan {
+                return Err(SelectionMilpError::DeadlineReached);
+            }
+            only_read_back = true;
+        }
+        model.set_option("time_limit", seconds.max(MOMENT_TO_READ_THE_START_S));
+    }
+    let started = Instant::now();
+    let solved = model
+        .try_solve()
+        .map_err(|status| SelectionMilpError::Solver(format!("{status:?}")))?;
+    let solve_ms = started.elapsed().as_millis() as u64;
+    use highs::{HighsModelStatus, HighsSolutionStatus};
+    let status = match solved.status() {
+        HighsModelStatus::Optimal if only_read_back => SelectionMilpStatus::TimeLimit,
+        HighsModelStatus::Optimal => SelectionMilpStatus::Optimal,
+        HighsModelStatus::Infeasible | HighsModelStatus::UnboundedOrInfeasible => {
+            return Err(SelectionMilpError::Infeasible)
+        }
+        HighsModelStatus::ReachedTimeLimit
+        | HighsModelStatus::ReachedInterrupt
+        | HighsModelStatus::ReachedIterationLimit
+        | HighsModelStatus::ReachedMemoryLimit
+        | HighsModelStatus::ReachedSolutionLimit => SelectionMilpStatus::TimeLimit,
+        other => return Err(SelectionMilpError::Solver(format!("{other:?}"))),
+    };
+    if solved.primal_solution_status() != HighsSolutionStatus::Feasible {
+        return Err(SelectionMilpError::NoSolution {
+            status: format!("{:?}", solved.status()),
+        });
+    }
+    let dual_bound = solved
+        .double_info_value(c"mip_dual_bound")
+        .unwrap_or(f64::NAN)
+        / scale;
+    Ok(SolvedPart {
+        status,
+        columns: solved.get_solution().columns().to_vec(),
+        objective_per_second: solved.objective_value() / scale,
+        dual_bound_per_second: dual_bound,
+        report: PartSolve {
+            queries,
+            variables: part_columns.len(),
+            binaries,
+            rows: row_count,
+            status: format!("{status:?}"),
+            solve_ms,
+            dual_bound_per_second: dual_bound,
+        },
+    })
+}
+
 /// HiGHS's threads (its `threads` option).
 const HIGHS_THREADS: i32 = 4;
 
diff --git a/crates/plan-selection/src/milp/model/tests.rs b/crates/plan-selection/src/milp/model/tests.rs
index 861a11a2e..60e329da4 100644
--- a/crates/plan-selection/src/milp/model/tests.rs
+++ b/crates/plan-selection/src/milp/model/tests.rs
@@ -64,7 +64,7 @@ pub(crate) fn through_stages(queries: &[Query<'_>], raw_data_retained: bool) ->
 /// [`through_stages`] over metric `m` declared non-negative (so relative
 /// errors compose through quotients, and Pass 1 keeps DDSketch).
 pub(crate) fn through_stages_non_negative(queries: &[Query<'_>]) -> Fixture {
-    through_stages_over(queries, true, false, true)
+    through_stages_over(queries, &["m"], true, false, true)
 }
 
 /// [`through_stages`], with rows arriving (100 per second) when
@@ -74,12 +74,14 @@ pub(crate) fn through_stages_with(
     raw_data_retained: bool,
     ingesting: bool,
 ) -> Fixture {
-    through_stages_over(queries, raw_data_retained, ingesting, false)
+    through_stages_over(queries, &["m"], raw_data_retained, ingesting, false)
 }
 
-/// The fixture builder behind [`through_stages`].
+/// The fixture builder behind [`through_stages`], over the metrics
+/// `metrics`, each with 100 series sampled every second.
 fn through_stages_over(
     queries: &[Query<'_>],
+    metrics: &[&str],
     raw_data_retained: bool,
     ingesting: bool,
     values_non_negative: bool,
@@ -115,9 +117,15 @@ fn through_stages_over(
             entry
         })
         .collect();
-    let mut data = DataWorkload::sampled_every(["m".to_string()], DurationMs(1_000));
-    data.sources.get_mut("m").expect("m").series_count = Some(100);
-    data.sources.get_mut("m").expect("m").values_non_negative = values_non_negative;
+    let mut data = DataWorkload::sampled_every(
+        metrics.iter().map(|metric| metric.to_string()),
+        DurationMs(1_000),
+    );
+    for metric in metrics {
+        let source = data.sources.get_mut(*metric).expect("declared");
+        source.series_count = Some(100);
+        source.values_non_negative = values_non_negative;
+    }
     if ingesting {
         data.arrival = asap_types::workload::DataArrival::Ingesting {
             rows_per_second: 100.0,
@@ -176,15 +184,17 @@ pub(crate) fn solve_and_check(
     cost_inputs: &CostInputs,
     restrict: impl FnOnce(&mut SelectionMilp<'_>),
 ) -> Result<Checked, SelectionMilpError> {
-    solve_and_check_by(fixture, cost_inputs, restrict, None)
+    solve_and_check_by(fixture, cost_inputs, restrict, None, true)
 }
 
-/// [`solve_and_check`] with HiGHS's time limit `limit`.
+/// [`solve_and_check`] with HiGHS's time limit `limit`, and the
+/// independent parts solved apart when `split_into_parts`.
 fn solve_and_check_by(
     fixture: &Fixture,
     cost_inputs: &CostInputs,
     restrict: impl FnOnce(&mut SelectionMilp<'_>),
     limit: Option<Duration>,
+    split_into_parts: bool,
 ) -> Result<Checked, SelectionMilpError> {
     let inputs = prepare_milp_inputs(
         &fixture.physical,
@@ -199,6 +209,7 @@ fn solve_and_check_by(
         cost_inputs,
         SelectionMilpOptions {
             enforce_accuracy: true,
+            split_into_parts,
         },
     )?;
     assert!(milp.inexact.is_empty(), "inexact terms: {:?}", milp.inexact);
@@ -642,6 +653,7 @@ fn a_phantom_reader_does_not_satisfy_minimum_query_readers() {
     let inputs = prepare_milp_inputs(&fixture.physical, &fixture.entries, &fixture.data, &costs);
     let options = SelectionMilpOptions {
         enforce_accuracy: true,
+        split_into_parts: true,
     };
     let mut milp = build_selection_milp(
         &fixture.physical,
@@ -775,6 +787,7 @@ fn cost_models_2_and_3_are_unsupported() {
             prepare_milp_inputs(&fixture.physical, &fixture.entries, &fixture.data, &costs);
         let options = SelectionMilpOptions {
             enforce_accuracy: true,
+            split_into_parts: true,
         };
         assert!(matches!(
             build_selection_milp(
@@ -825,6 +838,7 @@ fn accuracy_rows_follow_the_composition_on_a_ratio() {
             &costs,
             SelectionMilpOptions {
                 enforce_accuracy: true,
+                split_into_parts: true,
             },
         )
         .unwrap();
@@ -911,8 +925,14 @@ fn the_warm_start_is_a_feasible_plan_costed_exactly() {
 // time ends with a validated incumbent, reported as the deadline.
 #[test]
 fn a_deadline_yields_the_warm_starts_plan() {
-    let checked =
-        solve_and_check_by(&two_windows(), &cost_inputs(), |_| {}, Some(Duration::ZERO)).unwrap();
+    let checked = solve_and_check_by(
+        &two_windows(),
+        &cost_inputs(),
+        |_| {},
+        Some(Duration::ZERO),
+        true,
+    )
+    .unwrap();
     assert_eq!(checked.solution.status, SelectionMilpStatus::TimeLimit);
     let start = checked
         .solution
@@ -956,3 +976,52 @@ fn a_widened_scan_reading_its_own_memo_group_plans_exactly() {
         assert_eq!(checked.solution.status, SelectionMilpStatus::Optimal);
     }
 }
+
+// #1095 §5 (3a): queries over two metrics share no candidate, run of `I` or
+// source, so the model splits into one part per metric; solving the parts
+// apart gives the plan and cost of one solve of the whole model.
+#[test]
+fn independent_query_sets_solve_apart_to_the_same_plan() {
+    let queries = [
+        query("sum_over_time(m[5m])", 60_000),
+        query("sum_over_time(m[1h])", 60_000),
+        query("sum(n) + max(n)", 30_000),
+    ];
+    let fixture = through_stages_over(&queries, &["m", "n"], true, false, false);
+    let costs = cost_inputs();
+    let inputs = prepare_milp_inputs(&fixture.physical, &fixture.entries, &fixture.data, &costs);
+    let options = SelectionMilpOptions {
+        enforce_accuracy: true,
+        split_into_parts: true,
+    };
+    let milp = build_selection_milp(
+        &fixture.physical,
+        &inputs,
+        &fixture.entries,
+        &costs,
+        options,
+    )
+    .unwrap();
+    let parts = independent_parts(&milp);
+    let queries_per_part: Vec<Vec<usize>> = parts.iter().map(|part| part.queries.clone()).collect();
+    assert_eq!(queries_per_part, vec![vec![0, 1], vec![2]]);
+    let columns: usize = parts.iter().map(|part| part.columns.len()).sum();
+    let rows: usize = parts.iter().map(|part| part.rows.len()).sum();
+    assert_eq!(columns, milp.builder.variables.len());
+    assert_eq!(rows, milp.builder.rows.len());
+
+    let apart = solve_and_check_by(&fixture, &costs, |_| {}, None, true).unwrap();
+    let whole = solve_and_check_by(&fixture, &costs, |_| {}, None, false).unwrap();
+    assert_eq!(apart.solution.parts.len(), 2);
+    assert_eq!(whole.solution.parts.len(), 1);
+    assert_eq!(apart.solution.status, SelectionMilpStatus::Optimal);
+    assert_eq!(whole.solution.status, SelectionMilpStatus::Optimal);
+    let tolerance = 1e-9 * whole.solution.cost_per_second.abs();
+    assert!(
+        (apart.solution.cost_per_second - whole.solution.cost_per_second).abs() <= tolerance,
+        "apart {:e} != whole {:e}",
+        apart.solution.cost_per_second,
+        whole.solution.cost_per_second
+    );
+    assert_eq!(apart.dag.candidate_of_node, whole.dag.candidate_of_node);
+}
diff --git a/crates/plan-selection/src/select.rs b/crates/plan-selection/src/select.rs
index 966be802d..877341c6b 100644
--- a/crates/plan-selection/src/select.rs
+++ b/crates/plan-selection/src/select.rs
@@ -98,7 +98,10 @@ pub(crate) fn select_plan_with_milp(
         &inputs,
         entries,
         cost_inputs,
-        SelectionMilpOptions { enforce_accuracy },
+        SelectionMilpOptions {
+            enforce_accuracy,
+            split_into_parts: true,
+        },
     )
     .map_err(selection_error)?;
     let build_ms = construction_started.elapsed().as_millis() as u64;
@@ -336,6 +339,7 @@ fn report_run(
         "plan_cost_per_second": plan_cost_per_second,
         "dual_bound_per_second": solution.cost_lower_bound_per_second,
         "relative_gap": solution.relative_gap,
+        "parts": solution.parts,
     });
     use std::io::Write;
     if let Ok(mut file) = std::fs::OpenOptions::new()

From 2f8500d951c8c742bfb21513d8024e61f9a4ff6e Mon Sep 17 00:00:00 2001
From: zzylol <50204836+zzylol@users.noreply.github.com>
Date: Sat, 31 Oct 2026 17:43:24 +0000
Subject: [PATCH 2/4] =?UTF-8?q?feat(plan-selection):=20LP=20relaxation=20d?=
 =?UTF-8?q?iagnosis=20behind=20the=20lp-diagnosis=20feature=20(#1095=20?=
 =?UTF-8?q?=C2=A76.1)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

With ASAP_STAGE3_LP_DIAGNOSIS=<path>, a build with the feature appends,
per solve, the LP relaxation and root-node bounds, the LP with each row
family dropped, the LP with cover rows on configuration and window
terms, and the objective split by variable family at the LP, the warm
start and the incumbent. Measurement only; production builds leave it out.
---
 crates/devtools/Cargo.toml                    |   5 +
 crates/plan-selection/Cargo.toml              |   4 +
 .../plan-selection/src/milp/lp_diagnosis.rs   | 314 ++++++++++++++++++
 crates/plan-selection/src/milp/mod.rs         |   6 +-
 crates/plan-selection/src/milp/model.rs       |  41 +++
 5 files changed, 369 insertions(+), 1 deletion(-)
 create mode 100644 crates/plan-selection/src/milp/lp_diagnosis.rs

diff --git a/crates/devtools/Cargo.toml b/crates/devtools/Cargo.toml
index 65a65831e..8c71ffb2f 100644
--- a/crates/devtools/Cargo.toml
+++ b/crates/devtools/Cargo.toml
@@ -32,6 +32,11 @@ serde_json = "1"
 serde_yaml = "0.9"
 tokio = { version = "1", features = ["rt", "macros", "rt-multi-thread"] }
 
+[features]
+# Each Stage 3 selection MILP solve also reports where its LP relaxation
+# loses its bound, with `ASAP_STAGE3_LP_DIAGNOSIS=<path>` (measurement only).
+lp-diagnosis = ["asap-plan-selection/lp-diagnosis"]
+
 [dev-dependencies]
 asap-types = { path = "../types", features = ["test-support"] }
 asap-physical-asap-optimization = { path = "../physical-asap-optimization", features = ["test-support"] }
diff --git a/crates/plan-selection/Cargo.toml b/crates/plan-selection/Cargo.toml
index 193554f1f..557755d49 100644
--- a/crates/plan-selection/Cargo.toml
+++ b/crates/plan-selection/Cargo.toml
@@ -25,6 +25,10 @@ highs = "=2.4.0"
 # `reference_pricer`: the whole-workload selector over enumerated plans, for
 # the developer tools' reference pricing (`asap-devtools`).
 reference-pricer = []
+# `lp-diagnosis`: with `ASAP_STAGE3_LP_DIAGNOSIS=<path>`, each selection MILP
+# solve also reports where its LP relaxation loses its bound (measurement
+# only, for the developer tools).
+lp-diagnosis = []
 
 [dev-dependencies]
 asap-types = { path = "../types", features = ["test-support"] }
diff --git a/crates/plan-selection/src/milp/lp_diagnosis.rs b/crates/plan-selection/src/milp/lp_diagnosis.rs
new file mode 100644
index 000000000..5f2c58f7d
--- /dev/null
+++ b/crates/plan-selection/src/milp/lp_diagnosis.rs
@@ -0,0 +1,314 @@
+//! #1095 §6.1: where the selection MILP's LP relaxation loses its bound
+//! (measurement only, behind the `lp-diagnosis` feature). With
+//! `ASAP_STAGE3_LP_DIAGNOSIS=<path>`, each solve appends one JSON line:
+//! - `lp`: the LP relaxation of the model as built, and `root`: HiGHS's
+//!   bound after its root node (presolve, cuts), both from the warm start;
+//! - `lp_without`: the LP with one row family dropped, per family: how much
+//!   of the LP bound that family's rows hold up;
+//! - `lp_tightened`: the LP with valid rows added: `Σ terms ≥ run` for each
+//!   run charged only on its configuration or window terms (at integral
+//!   binaries one of them is 1);
+//! - `objective_by_family`: the objective split by the family of each
+//!   charged variable, at the LP's solution, at the warm start and at the
+//!   incumbent, so `incumbent − lp` splits into one share per family;
+//! - `runs`: the charged `p` and `x` variables with the largest difference
+//!   between the incumbent and the LP, with their values in both.
+use std::collections::BTreeMap;
+use std::time::Instant;
+
+use good_lp::solvers::highs::highs;
+use good_lp::{constraint, variable, Constraint, Expression, ProblemVariables, SolverModel};
+
+use super::model::SelectionMilp;
+
+/// The environment variable naming the file the diagnosis appends to.
+const DIAGNOSIS_ENV: &str = "ASAP_STAGE3_LP_DIAGNOSIS";
+
+/// The runs listed per solve.
+const RUNS_SHOWN: usize = 40;
+
+/// The time limit of each LP and of the root node.
+const LIMIT_S: f64 = 300.0;
+
+/// One LP (or root node) solved.
+#[derive(Debug, Clone, serde::Serialize)]
+struct Bound {
+    /// The bound, cost per second (`None` when HiGHS gave none).
+    per_second: Option<f64>,
+    status: String,
+    ms: u64,
+}
+
+/// The LP solves, before the MILP's.
+pub(crate) struct Diagnosis {
+    path: String,
+    lp: Bound,
+    root: Bound,
+    lp_without: BTreeMap<&'static str, Bound>,
+    lp_tightened: Bound,
+    tightening_rows: usize,
+    /// The LP's columns.
+    lp_columns: Vec<f64>,
+    /// The objective, by column.
+    objective: Vec<(usize, f64)>,
+    column_families: Vec<&'static str>,
+    /// Per `p` and `x` column, what it stands for.
+    labels: BTreeMap<usize, String>,
+}
+
+impl Diagnosis {
+    /// The LP solves of `milp` from `start` (columns), when the diagnosis is
+    /// asked for.
+    pub(crate) fn run(milp: &SelectionMilp<'_>, start: Option<&[f64]>) -> Option<Self> {
+        let path = std::env::var(DIAGNOSIS_ENV).ok()?;
+        let definitions: Vec<_> = milp
+            .builder
+            .variables
+            .iter_variables_with_def()
+            .map(|(variable, definition)| (variable, definition.clone()))
+            .collect();
+        let column_of: std::collections::HashMap<_, usize> = definitions
+            .iter()
+            .enumerate()
+            .map(|(column, (variable, _))| (*variable, column))
+            .collect();
+        let objective: Vec<(usize, f64)> = milp
+            .objective_terms
+            .iter()
+            .map(|&(variable, coefficient)| (column_of[&variable], coefficient))
+            .collect();
+        let scale = scale_of(&objective, start);
+        // The relaxation: the same variables (same indices), none integral.
+        let relaxed = || {
+            let mut variables = ProblemVariables::new();
+            for (_, definition) in &definitions {
+                variables.add(
+                    variable()
+                        .min(definition.get_min())
+                        .max(definition.get_max()),
+                );
+            }
+            variables
+        };
+        let objective_expression: Expression = milp
+            .objective_terms
+            .iter()
+            .map(|&(variable, coefficient)| coefficient * scale * variable)
+            .sum();
+        let solve = |variables: ProblemVariables, rows: Vec<Constraint>, root: bool| {
+            let started = Instant::now();
+            let problem = variables
+                .minimise(objective_expression.clone())
+                .using(highs)
+                .with_all(rows);
+            let mut model = problem.try_into_inner().expect("a HiGHS model");
+            model.set_option("output_flag", false);
+            model.set_option("threads", 4);
+            model.set_option("time_limit", LIMIT_S);
+            if root {
+                // The root node, then stop: its bound is the root bound.
+                model.set_option("mip_max_nodes", 1);
+                model.set_option("mip_rel_gap", 0.0);
+                model.set_option("mip_abs_gap", 0.0);
+                if let Some(start) = start {
+                    let _ = model.try_set_solution(Some(start), None, None, None);
+                }
+            }
+            let Ok(solved) = model.try_solve() else {
+                return (
+                    Bound {
+                        per_second: None,
+                        status: "solver error".into(),
+                        ms: started.elapsed().as_millis() as u64,
+                    },
+                    Vec::new(),
+                );
+            };
+            let value = if root {
+                solved.double_info_value(c"mip_dual_bound").ok()
+            } else if solved.status() == highs::HighsModelStatus::Optimal {
+                Some(solved.objective_value())
+            } else {
+                None
+            };
+            let columns = if root {
+                Vec::new()
+            } else {
+                solved.get_solution().columns().to_vec()
+            };
+            (
+                Bound {
+                    per_second: value.filter(|v| v.is_finite()).map(|v| v / scale),
+                    status: format!("{:?}", solved.status()),
+                    ms: started.elapsed().as_millis() as u64,
+                },
+                columns,
+            )
+        };
+        let rows = &milp.builder.rows;
+        let families = &milp.builder.row_families;
+        let (lp, lp_columns) = solve(relaxed(), rows.clone(), false);
+        let (root, _) = solve(milp.builder.variables.clone(), rows.clone(), true);
+        let mut lp_without = BTreeMap::new();
+        let present: std::collections::BTreeSet<&'static str> = families.iter().copied().collect();
+        for family in present {
+            let kept: Vec<Constraint> = rows
+                .iter()
+                .zip(families)
+                .filter(|(_, &f)| f != family)
+                .map(|(row, _)| row.clone())
+                .collect();
+            lp_without.insert(family, solve(relaxed(), kept, false).0);
+        }
+        let tightening = cover_rows(milp);
+        let tightening_rows = tightening.len();
+        let mut tightened = rows.clone();
+        tightened.extend(tightening);
+        let (lp_tightened, _) = solve(relaxed(), tightened, false);
+        let candidate = |c: asap_types::stages::PhysicalCandidateId| {
+            let data = &milp.physical.candidates[c.0];
+            format!(
+                "c{} {} {:?} {:?}",
+                c.0,
+                milp.physical.physical_node(c).operator.kind_name(),
+                data.timing,
+                data.materialization
+            )
+        };
+        let mut labels = BTreeMap::new();
+        for (context, choices) in milp.choice.iter().enumerate() {
+            let place = match milp.inputs.contexts[context].kind {
+                crate::groups::ContextKind::Query(query) => format!("q{query}"),
+                crate::groups::ContextKind::Ingestion => "I".to_string(),
+                crate::groups::ContextKind::Shared(m) => format!("m{}", m.0),
+            };
+            for (&c, variable) in choices {
+                labels.insert(column_of[variable], format!("p[{place}] {}", candidate(c)));
+            }
+        }
+        for (&m, variable) in &milp.selected {
+            labels.insert(column_of[variable], format!("x {}", candidate(m)));
+        }
+        Some(Self {
+            labels,
+            path,
+            lp,
+            root,
+            lp_without,
+            lp_tightened,
+            tightening_rows,
+            lp_columns,
+            objective,
+            column_families: milp.builder.column_families.clone(),
+        })
+    }
+
+    /// Appends the line, with the warm start's and the incumbent's columns
+    /// (completed) and the MILP's results.
+    pub(crate) fn write(
+        self,
+        start: Option<&[f64]>,
+        incumbent: Option<&[f64]>,
+        objective_per_second: f64,
+        dual_bound_per_second: f64,
+        status: &str,
+    ) {
+        let by_family = |columns: Option<&[f64]>| -> Option<BTreeMap<&'static str, f64>> {
+            let columns = columns.filter(|columns| !columns.is_empty())?;
+            let mut split: BTreeMap<&'static str, f64> = BTreeMap::new();
+            for &(column, coefficient) in &self.objective {
+                *split.entry(self.column_families[column]).or_default() +=
+                    coefficient * columns[column];
+            }
+            Some(split)
+        };
+        // Per labeled column: (LP charge, incumbent charge, LP value,
+        // incumbent value).
+        let mut runs: Vec<(String, f64, f64, f64, f64)> = Vec::new();
+        if let Some(incumbent) = incumbent.filter(|columns| !columns.is_empty()) {
+            let mut charge: BTreeMap<usize, f64> = BTreeMap::new();
+            for &(column, coefficient) in &self.objective {
+                if self.labels.contains_key(&column) {
+                    *charge.entry(column).or_default() += coefficient;
+                }
+            }
+            for (column, coefficient) in charge {
+                let (lp, at) = (self.lp_columns[column], incumbent[column]);
+                runs.push((
+                    self.labels[&column].clone(),
+                    coefficient * lp,
+                    coefficient * at,
+                    lp,
+                    at,
+                ));
+            }
+            runs.sort_by(|a, b| (b.2 - b.1).abs().total_cmp(&(a.2 - a.1).abs()));
+            runs.truncate(RUNS_SHOWN);
+        }
+        let line = serde_json::json!({
+            "runs": runs,
+            "variables": self.column_families.len(),
+            "status": status,
+            "objective_per_second": objective_per_second,
+            "dual_bound_per_second": dual_bound_per_second,
+            "lp": self.lp,
+            "root": self.root,
+            "lp_without": self.lp_without,
+            "lp_tightened": self.lp_tightened,
+            "tightening_rows": self.tightening_rows,
+            "objective_by_family": {
+                "lp": by_family(Some(&self.lp_columns)),
+                "warm_start": by_family(start),
+                "incumbent": by_family(incumbent),
+            },
+        });
+        use std::io::Write;
+        if let Ok(mut file) = std::fs::OpenOptions::new()
+            .create(true)
+            .append(true)
+            .open(&self.path)
+        {
+            let _ = writeln!(file, "{line}");
+        }
+    }
+}
+
+/// `Σ terms ≥ run` for each run whose terms are all configuration or
+/// window terms (none on the run's own indicator): at integral binaries a
+/// used run has one term at 1.
+fn cover_rows(milp: &SelectionMilp<'_>) -> Vec<Constraint> {
+    let mut rows = Vec::new();
+    for (&(context, candidate), terms) in &milp.run_terms_cache {
+        if terms.is_empty() || terms.iter().any(|term| term.plain) {
+            continue;
+        }
+        let Some(run) = milp.run_indicator(context, candidate) else {
+            continue;
+        };
+        let sum: Expression = terms.iter().map(|term| term.indicator.clone()).sum();
+        rows.push(constraint!(sum >= run));
+    }
+    rows
+}
+
+/// The objective's scale for HiGHS, as the MILP's solve sets it.
+fn scale_of(objective: &[(usize, f64)], start: Option<&[f64]>) -> f64 {
+    let largest = objective
+        .iter()
+        .map(|(_, coefficient)| coefficient.abs())
+        .fold(0.0, f64::max);
+    if !(largest > 0.0 && largest.is_finite()) {
+        return 1.0;
+    }
+    let start_cost: f64 = start.map_or(0.0, |start| {
+        objective
+            .iter()
+            .map(|&(column, coefficient)| coefficient * start[column])
+            .sum()
+    });
+    if start_cost > 0.0 && start_cost.is_finite() {
+        (1_000.0 / start_cost).min(1e7 / largest)
+    } else {
+        1_000.0 / largest
+    }
+}
diff --git a/crates/plan-selection/src/milp/mod.rs b/crates/plan-selection/src/milp/mod.rs
index c0cf8b822..7c0e9ab2f 100644
--- a/crates/plan-selection/src/milp/mod.rs
+++ b/crates/plan-selection/src/milp/mod.rs
@@ -12,11 +12,15 @@
 //! - [`completion`]: the exact objective of a plan given as binaries, the
 //!   cheapest values of the continuous variables found without a solver;
 //! - [`decode`]: the solution as a `DagOfSelectedCandidates`, one node per
-//!   candidate run.
+//!   candidate run;
+//! - `lp_diagnosis` (feature `lp-diagnosis`, measurement only): where the
+//!   LP relaxation loses its bound.
 //!
 //! [`crate::select`] runs them, realizes the plan once and checks it.
 pub(crate) mod completion;
 pub(crate) mod cost_terms;
 pub(crate) mod decode;
+#[cfg(feature = "lp-diagnosis")]
+pub(crate) mod lp_diagnosis;
 pub(crate) mod model;
 pub(crate) mod warm_start;
diff --git a/crates/plan-selection/src/milp/model.rs b/crates/plan-selection/src/milp/model.rs
index 2c2755795..eadc12b6f 100644
--- a/crates/plan-selection/src/milp/model.rs
+++ b/crates/plan-selection/src/milp/model.rs
@@ -125,6 +125,9 @@ pub(crate) struct ModelBuilder {
     pub(crate) scopes: BTreeMap<Scope, FamilySize>,
     /// Where new variables and rows are counted.
     pub(crate) scope: Scope,
+    /// The family of each variable, in creation order.
+    #[cfg(feature = "lp-diagnosis")]
+    pub(crate) column_families: Vec<&'static str>,
 }
 
 impl ModelBuilder {
@@ -136,10 +139,16 @@ impl ModelBuilder {
             families: BTreeMap::new(),
             scopes: BTreeMap::new(),
             scope: Scope::Global,
+            #[cfg(feature = "lp-diagnosis")]
+            column_families: Vec::new(),
         }
     }
 
     fn count(&mut self, family: &'static str, size: FamilySize) {
+        #[cfg(feature = "lp-diagnosis")]
+        if size.variables > 0 {
+            self.column_families.push(family);
+        }
         self.families.entry(family).or_default().add(size);
         self.scopes.entry(self.scope).or_default().add(size);
     }
@@ -1826,6 +1835,25 @@ pub(crate) fn solve_selection_milp(
     let start = warm_start::warm_start(&milp, &completion, Some(search_budget));
     let warm_start_ms = warm_started.elapsed().as_millis() as u64;
     let warm_start_cost_per_second = start.as_ref().map(|start| start.cost_per_second);
+    #[cfg(feature = "lp-diagnosis")]
+    let diagnosis = {
+        let column_of: HashMap<Variable, usize> = milp
+            .builder
+            .variables
+            .iter_variables_with_def()
+            .enumerate()
+            .map(|(column, (variable, _))| (variable, column))
+            .collect();
+        let start_columns = start.as_ref().map(|start| {
+            let mut columns = vec![0.0; column_of.len()];
+            for &(variable, value) in &start.values {
+                columns[column_of[&variable]] = value;
+            }
+            columns
+        });
+        super::lp_diagnosis::Diagnosis::run(&milp, start_columns.as_deref())
+            .map(|diagnosis| (diagnosis, start_columns))
+    };
     let parts = if milp.split_into_parts {
         independent_parts(&milp)
     } else {
@@ -1986,6 +2014,19 @@ pub(crate) fn solve_selection_milp(
     warm_start::debug(|| {
         format!("objective: HiGHS {highs_objective:e}, completed {objective_value:e}")
     });
+    #[cfg(feature = "lp-diagnosis")]
+    if let Some((diagnosis, start_columns)) = diagnosis {
+        let incumbent = completion.complete_columns(&columns).ok();
+        diagnosis.write(
+            start_columns.as_deref(),
+            incumbent
+                .as_ref()
+                .map(|completed| completed.columns.as_slice()),
+            objective_value,
+            dual_bound,
+            &format!("{status:?}"),
+        );
+    }
     let on = |variable: Variable| columns[column_of[&variable]] > 0.5;
     Ok(SelectionMilpSolution {
         status,

From bf4e240e70ffc04232eab8481947a8d02f48ccb6 Mon Sep 17 00:00:00 2001
From: zzylol <50204836+zzylol@users.noreply.github.com>
Date: Sat, 31 Oct 2026 17:56:16 +0000
Subject: [PATCH 3/4] =?UTF-8?q?perf(plan-selection):=20tighter=20big-M=20b?=
 =?UTF-8?q?ounds=20for=20context=20execution=20and=20accuracy=20shares=20(?=
 =?UTF-8?q?#1095=20=C2=A76.1)?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Ē[m] in V ≥ E − Ē (x[m] − z) is the sum over the context's slots of its
costliest run, not over every candidate; N̄_q in the accuracy row is the
number of places (context, slot) q can depend on an approximate value
at, not the number of runs. Both stay valid upper bounds, so the
objective still equals plan_cost at Optimal.
---
 crates/plan-selection/src/milp/cost_terms.rs | 29 ++++++++++++++------
 crates/plan-selection/src/milp/model.rs      | 10 ++++++-
 2 files changed, 30 insertions(+), 9 deletions(-)

diff --git a/crates/plan-selection/src/milp/cost_terms.rs b/crates/plan-selection/src/milp/cost_terms.rs
index 3be9a9ecb..99b12ceb8 100644
--- a/crates/plan-selection/src/milp/cost_terms.rs
+++ b/crates/plan-selection/src/milp/cost_terms.rs
@@ -63,7 +63,9 @@ use asap_types::workload::bounds::source_name;
 use asap_types::workload::{QueryRecurrence, RepeatedDemand};
 use good_lp::{constraint, Expression, IntoAffineExpression, Variable};
 
-use crate::groups::{CandidateUsage, ContextKind, RunKind, UsageAndOutput, UsageContext};
+use crate::groups::{
+    CandidateUsage, ContextKind, ResultSlot, RunKind, UsageAndOutput, UsageContext,
+};
 
 use super::model::{Reader, Scope, SelectionMilp, SelectionMilpError};
 
@@ -985,8 +987,10 @@ pub(crate) fn encode_objective(milp: &mut SelectionMilp<'_>) -> Result<(), Selec
         }
     };
     // Per query-time shareable context: the execution cost of its runs
-    // other than `m`'s plain own term, and the most it can be.
-    let mut context_execution: BTreeMap<usize, (Expression, f64)> = BTreeMap::new();
+    // other than `m`'s plain own term, and per slot (`None`: `m`'s own run)
+    // the most it can be.
+    let mut context_execution: BTreeMap<usize, (Expression, BTreeMap<Option<ResultSlot>, f64>)> =
+        BTreeMap::new();
     // `H[m]` per query-time shareable context.
     let mut held_time: HashMap<usize, (Variable, f64)> = HashMap::new();
     let runs = milp.runs();
@@ -996,6 +1000,9 @@ pub(crate) fn encode_objective(milp: &mut SelectionMilp<'_>) -> Result<(), Selec
         let candidate_data = &milp.physical.candidates[candidate.0];
         let (timing, materialization) = (candidate_data.timing, candidate_data.materialization);
         let own = milp.inputs.contexts[context].kind == ContextKind::Shared(candidate);
+        // The most this run adds to its context's execution cost: all its
+        // terms at 1.
+        let mut run_execution = 0.0;
         for term in &terms {
             let numbers = term.numbers;
             charge(&term.indicator, numbers.standing_cost);
@@ -1008,11 +1015,11 @@ pub(crate) fn encode_objective(milp: &mut SelectionMilp<'_>) -> Result<(), Selec
                             charge(&component.into(), numbers.execution_cost * rate);
                         }
                     } else {
-                        let (sum, most) = context_execution
+                        let (sum, _) = context_execution
                             .entry(context)
-                            .or_insert_with(|| (Expression::from(0.0), 0.0));
+                            .or_insert_with(|| (Expression::from(0.0), BTreeMap::new()));
                         *sum += numbers.execution_cost * term.indicator.clone();
-                        *most += numbers.execution_cost;
+                        run_execution += numbers.execution_cost;
                     }
                 }
             }
@@ -1093,11 +1100,17 @@ pub(crate) fn encode_objective(milp: &mut SelectionMilp<'_>) -> Result<(), Selec
                 ContextKind::Ingestion => {}
             }
         }
+        if let Some((_, most_per_slot)) = context_execution.get_mut(&context) {
+            let slot = (!own).then(|| milp.inputs.slot[candidate.0]);
+            let entry = most_per_slot.entry(slot).or_default();
+            *entry = entry.max(run_execution);
+        }
     }
     // A context's execution at its varying rate: `V ≥ E − Ē (x[m] − z)` per
-    // component `z`.
+    // component `z`, `Ē` the most `E` can be: per slot, its costliest run.
     milp.builder.scope = Scope::SharedContexts;
-    for (context, (sum, most)) in context_execution {
+    for (context, (sum, most_per_slot)) in context_execution {
+        let most: f64 = most_per_slot.values().sum();
         let ContextRate::Components(components) = &rates[context] else {
             continue;
         };
diff --git a/crates/plan-selection/src/milp/model.rs b/crates/plan-selection/src/milp/model.rs
index eadc12b6f..d481a1295 100644
--- a/crates/plan-selection/src/milp/model.rs
+++ b/crates/plan-selection/src/milp/model.rs
@@ -1579,6 +1579,9 @@ fn accuracy_rows(milp: &mut SelectionMilp<'_>) {
         let exact_target = milp.entries[query].requirements.accuracy.is_exact();
         // (candidate, dep, parts) per run `q` may depend on.
         let mut counted: Vec<(PhysicalCandidateId, Expression, usize)> = Vec::new();
+        // Where each counted run gives its value: (context, slot), `None`
+        // for a shareable candidate's own run. One value per place.
+        let mut places: BTreeSet<(usize, Option<ResultSlot>)> = BTreeSet::new();
         for &(context, candidate) in &approximate_runs {
             let Some(weight) = inputs.accuracy_weights[query].get(&candidate).copied() else {
                 continue;
@@ -1621,6 +1624,10 @@ fn accuracy_rows(milp: &mut SelectionMilp<'_>) {
                 AccuracyWeight::Unusable(_) => false,
             };
             if usable {
+                places.insert(match inputs.contexts[context].kind {
+                    ContextKind::Shared(own) if own == candidate => (context, None),
+                    _ => (context, Some(inputs.slot[candidate.0])),
+                });
                 counted.push((
                     candidate,
                     dep,
@@ -1642,7 +1649,8 @@ fn accuracy_rows(milp: &mut SelectionMilp<'_>) {
             milp.builder
                 .row("8 accuracy (pair)", constraint!(one + other <= 1.0));
         }
-        let most = counted.len() as f64;
+        // The most values `q` can depend on together.
+        let most = places.len() as f64;
         if counted.iter().all(|&(_, _, parts)| parts as f64 >= most) {
             continue;
         }

From 55ed1d099a6ac225fe6685acd633ec1b90962fe5 Mon Sep 17 00:00:00 2001
From: zzylol <50204836+zzylol@users.noreply.github.com>
Date: Sat, 31 Oct 2026 19:03:39 +0000
Subject: [PATCH 4/4] fix(plan-selection): review fixes for the independent
 parts and the LP diagnosis (#1095)

- The LP diagnosis pushes the MILP's deadline back by its own time.
- A part without binaries reports its LP optimum as its dual bound
  (HiGHS leaves mip_dual_bound at 0 for an LP).
- The comment on a row without variables names the part it goes to.
- A test pins the tighter context execution bound against the old one
  (same plan and cost); a test-only switch rebuilds the old bounds.
- lp_diagnosis reuses objective_scale and HIGHS_THREADS.
---
 crates/plan-selection/Cargo.toml              |  3 +-
 crates/plan-selection/src/milp/cost_terms.rs  | 10 ++-
 .../plan-selection/src/milp/lp_diagnosis.rs   | 43 ++++------
 crates/plan-selection/src/milp/model.rs       | 48 ++++++++---
 crates/plan-selection/src/milp/model/tests.rs | 80 +++++++++++++++++++
 5 files changed, 146 insertions(+), 38 deletions(-)

diff --git a/crates/plan-selection/Cargo.toml b/crates/plan-selection/Cargo.toml
index 557755d49..ac64258a5 100644
--- a/crates/plan-selection/Cargo.toml
+++ b/crates/plan-selection/Cargo.toml
@@ -27,7 +27,8 @@ highs = "=2.4.0"
 reference-pricer = []
 # `lp-diagnosis`: with `ASAP_STAGE3_LP_DIAGNOSIS=<path>`, each selection MILP
 # solve also reports where its LP relaxation loses its bound (measurement
-# only, for the developer tools).
+# only, for the developer tools; a diagnosed run's timings are not
+# comparable with an undiagnosed one's).
 lp-diagnosis = []
 
 [dev-dependencies]
diff --git a/crates/plan-selection/src/milp/cost_terms.rs b/crates/plan-selection/src/milp/cost_terms.rs
index 99b12ceb8..dbc2483ab 100644
--- a/crates/plan-selection/src/milp/cost_terms.rs
+++ b/crates/plan-selection/src/milp/cost_terms.rs
@@ -1110,7 +1110,15 @@ pub(crate) fn encode_objective(milp: &mut SelectionMilp<'_>) -> Result<(), Selec
     // component `z`, `Ē` the most `E` can be: per slot, its costliest run.
     milp.builder.scope = Scope::SharedContexts;
     for (context, (sum, most_per_slot)) in context_execution {
-        let most: f64 = most_per_slot.values().sum();
+        let most: f64 = if super::model::tests_loose_big_m() {
+            // Before: every term's cost.
+            sum.clone()
+                .linear_coefficients()
+                .map(|(_, cost)| cost)
+                .sum()
+        } else {
+            most_per_slot.values().sum()
+        };
         let ContextRate::Components(components) = &rates[context] else {
             continue;
         };
diff --git a/crates/plan-selection/src/milp/lp_diagnosis.rs b/crates/plan-selection/src/milp/lp_diagnosis.rs
index 5f2c58f7d..c946d61b4 100644
--- a/crates/plan-selection/src/milp/lp_diagnosis.rs
+++ b/crates/plan-selection/src/milp/lp_diagnosis.rs
@@ -13,13 +13,17 @@
 //!   incumbent, so `incumbent − lp` splits into one share per family;
 //! - `runs`: the charged `p` and `x` variables with the largest difference
 //!   between the incumbent and the LP, with their values in both.
+//!
+//! The deadline is pushed back by the diagnosis's time, so the MILP keeps
+//! its budget; a diagnosed run's wall time (and the facade's later rounds)
+//! is still not comparable with an undiagnosed one.
 use std::collections::BTreeMap;
 use std::time::Instant;
 
 use good_lp::solvers::highs::highs;
 use good_lp::{constraint, variable, Constraint, Expression, ProblemVariables, SolverModel};
 
-use super::model::SelectionMilp;
+use super::model::{objective_scale, SelectionMilp, HIGHS_THREADS};
 
 /// The environment variable naming the file the diagnosis appends to.
 const DIAGNOSIS_ENV: &str = "ASAP_STAGE3_LP_DIAGNOSIS";
@@ -28,7 +32,7 @@ const DIAGNOSIS_ENV: &str = "ASAP_STAGE3_LP_DIAGNOSIS";
 const RUNS_SHOWN: usize = 40;
 
 /// The time limit of each LP and of the root node.
-const LIMIT_S: f64 = 300.0;
+const LP_TIME_LIMIT_S: f64 = 300.0;
 
 /// One LP (or root node) solved.
 #[derive(Debug, Clone, serde::Serialize)]
@@ -77,7 +81,14 @@ impl Diagnosis {
             .iter()
             .map(|&(variable, coefficient)| (column_of[&variable], coefficient))
             .collect();
-        let scale = scale_of(&objective, start);
+        let start_values: Option<Vec<_>> = start.map(|start| {
+            definitions
+                .iter()
+                .map(|(variable, _)| *variable)
+                .zip(start.iter().copied())
+                .collect()
+        });
+        let scale = objective_scale(&milp.objective_terms, start_values.as_deref());
         // The relaxation: the same variables (same indices), none integral.
         let relaxed = || {
             let mut variables = ProblemVariables::new();
@@ -103,8 +114,8 @@ impl Diagnosis {
                 .with_all(rows);
             let mut model = problem.try_into_inner().expect("a HiGHS model");
             model.set_option("output_flag", false);
-            model.set_option("threads", 4);
-            model.set_option("time_limit", LIMIT_S);
+            model.set_option("threads", HIGHS_THREADS);
+            model.set_option("time_limit", LP_TIME_LIMIT_S);
             if root {
                 // The root node, then stop: its bound is the root bound.
                 model.set_option("mip_max_nodes", 1);
@@ -290,25 +301,3 @@ fn cover_rows(milp: &SelectionMilp<'_>) -> Vec<Constraint> {
     }
     rows
 }
-
-/// The objective's scale for HiGHS, as the MILP's solve sets it.
-fn scale_of(objective: &[(usize, f64)], start: Option<&[f64]>) -> f64 {
-    let largest = objective
-        .iter()
-        .map(|(_, coefficient)| coefficient.abs())
-        .fold(0.0, f64::max);
-    if !(largest > 0.0 && largest.is_finite()) {
-        return 1.0;
-    }
-    let start_cost: f64 = start.map_or(0.0, |start| {
-        objective
-            .iter()
-            .map(|&(column, coefficient)| coefficient * start[column])
-            .sum()
-    });
-    if start_cost > 0.0 && start_cost.is_finite() {
-        (1_000.0 / start_cost).min(1e7 / largest)
-    } else {
-        1_000.0 / largest
-    }
-}
diff --git a/crates/plan-selection/src/milp/model.rs b/crates/plan-selection/src/milp/model.rs
index d481a1295..ee987788f 100644
--- a/crates/plan-selection/src/milp/model.rs
+++ b/crates/plan-selection/src/milp/model.rs
@@ -1650,7 +1650,11 @@ fn accuracy_rows(milp: &mut SelectionMilp<'_>) {
                 .row("8 accuracy (pair)", constraint!(one + other <= 1.0));
         }
         // The most values `q` can depend on together.
-        let most = places.len() as f64;
+        let most = if tests_loose_big_m() {
+            counted.len() as f64
+        } else {
+            places.len() as f64
+        };
         if counted.iter().all(|&(_, _, parts)| parts as f64 >= most) {
             continue;
         }
@@ -1768,8 +1772,9 @@ pub(crate) fn independent_parts(milp: &SelectionMilp<'_>) -> Vec<IndependentPart
         }
     }
     for (index, row) in milp.builder.rows.iter().enumerate() {
-        // A row without a variable is a constant check: kept in the first
-        // part.
+        // A row without a variable is a constant check: kept in the part
+        // without a query when there is one (`None` sorts first), else in
+        // the first query's part.
         let key = row
             .expression()
             .linear_coefficients()
@@ -1843,6 +1848,10 @@ pub(crate) fn solve_selection_milp(
     let start = warm_start::warm_start(&milp, &completion, Some(search_budget));
     let warm_start_ms = warm_started.elapsed().as_millis() as u64;
     let warm_start_cost_per_second = start.as_ref().map(|start| start.cost_per_second);
+    // The diagnosis's LP solves do not count against the deadline: it is
+    // pushed back by their time.
+    #[cfg(feature = "lp-diagnosis")]
+    let diagnosis_started = Instant::now();
     #[cfg(feature = "lp-diagnosis")]
     let diagnosis = {
         let column_of: HashMap<Variable, usize> = milp
@@ -1862,6 +1871,8 @@ pub(crate) fn solve_selection_milp(
         super::lp_diagnosis::Diagnosis::run(&milp, start_columns.as_deref())
             .map(|diagnosis| (diagnosis, start_columns))
     };
+    #[cfg(feature = "lp-diagnosis")]
+    let deadline = deadline.map(|deadline| deadline + diagnosis_started.elapsed());
     let parts = if milp.split_into_parts {
         independent_parts(&milp)
     } else {
@@ -2159,10 +2170,17 @@ fn solve_part(
             status: format!("{:?}", solved.status()),
         });
     }
-    let dual_bound = solved
-        .double_info_value(c"mip_dual_bound")
-        .unwrap_or(f64::NAN)
-        / scale;
+    // HiGHS leaves `mip_dual_bound` at 0 for a part without binaries (an
+    // LP): its optimal objective is its bound.
+    let dual_bound = if binaries > 0 {
+        solved
+            .double_info_value(c"mip_dual_bound")
+            .unwrap_or(f64::NAN)
+    } else if status == SelectionMilpStatus::Optimal {
+        solved.objective_value()
+    } else {
+        f64::NAN
+    } / scale;
     Ok(SolvedPart {
         status,
         columns: solved.get_solution().columns().to_vec(),
@@ -2181,7 +2199,7 @@ fn solve_part(
 }
 
 /// HiGHS's threads (its `threads` option).
-const HIGHS_THREADS: i32 = 4;
+pub(crate) const HIGHS_THREADS: i32 = 4;
 
 /// HiGHS's time limit when the warm start used up the time to the deadline.
 const MOMENT_TO_READ_THE_START_S: f64 = 1e-3;
@@ -2192,7 +2210,7 @@ const WARM_START_MOST_TIME: std::time::Duration = std::time::Duration::from_secs
 /// The factor the objective is multiplied by for HiGHS: the warm start's
 /// cost becomes about 1,000, its largest coefficient at most 10⁷; without
 /// a warm start cost, the largest coefficient becomes 1,000.
-fn objective_scale(terms: &[(Variable, f64)], start: Option<&[(Variable, f64)]>) -> f64 {
+pub(crate) fn objective_scale(terms: &[(Variable, f64)], start: Option<&[(Variable, f64)]>) -> f64 {
     let largest = terms
         .iter()
         .map(|(_, coefficient)| coefficient.abs())
@@ -2212,5 +2230,17 @@ fn objective_scale(terms: &[(Variable, f64)], start: Option<&[(Variable, f64)]>)
     }
 }
 
+/// Whether this thread builds the big-M bounds as before #1095 step 1
+/// (tests only: to compare the tightened bounds with the looser ones).
+#[cfg(test)]
+pub(crate) fn tests_loose_big_m() -> bool {
+    tests::LOOSE_BIG_M.with(|loose| loose.get())
+}
+
+#[cfg(not(test))]
+pub(crate) fn tests_loose_big_m() -> bool {
+    false
+}
+
 #[cfg(test)]
 pub(crate) mod tests;
diff --git a/crates/plan-selection/src/milp/model/tests.rs b/crates/plan-selection/src/milp/model/tests.rs
index 60e329da4..ba87120eb 100644
--- a/crates/plan-selection/src/milp/model/tests.rs
+++ b/crates/plan-selection/src/milp/model/tests.rs
@@ -1025,3 +1025,83 @@ fn independent_query_sets_solve_apart_to_the_same_plan() {
     );
     assert_eq!(apart.dag.candidate_of_node, whole.dag.candidate_of_node);
 }
+
+thread_local! {
+    /// Set by a test to build the big-M bounds as before #1095 step 1
+    /// ([`super::tests_loose_big_m`]).
+    pub(crate) static LOOSE_BIG_M: std::cell::Cell<bool> = const { std::cell::Cell::new(false) };
+}
+
+/// Per row of family `family` in `fixture`'s model (built with the loose
+/// big-M bounds when `loose`), the largest coefficient on an `x` variable.
+fn big_m_per_row(fixture: &Fixture, family: &str, loose: bool) -> Vec<f64> {
+    LOOSE_BIG_M.with(|cell| cell.set(loose));
+    let costs = aws_fargate_cost_inputs();
+    let inputs = prepare_milp_inputs(&fixture.physical, &fixture.entries, &fixture.data, &costs);
+    let milp = build_selection_milp(
+        &fixture.physical,
+        &inputs,
+        &fixture.entries,
+        &costs,
+        SelectionMilpOptions {
+            enforce_accuracy: true,
+            split_into_parts: true,
+        },
+    )
+    .unwrap();
+    LOOSE_BIG_M.with(|cell| cell.set(false));
+    let selected: std::collections::HashSet<Variable> = milp.selected.values().copied().collect();
+    milp.builder
+        .rows
+        .iter()
+        .zip(&milp.builder.row_families)
+        .filter(|(_, &f)| f == family)
+        .map(|(row, _)| {
+            row.expression()
+                .linear_coefficients()
+                .filter(|(variable, _)| selected.contains(variable))
+                .map(|(_, coefficient)| coefficient.abs())
+                .fold(0.0, f64::max)
+        })
+        .collect()
+}
+
+// #1095 step 1: `Ē[m]` in `V ≥ E − Ē (x[m] − z)` is the sum over the
+// context's slots of its costliest run. With a rollup read every 20 s and
+// every 30 s, a shared context has slots with several runs, so the bound is
+// below the old sum over every run; the plan and its cost are unchanged.
+#[test]
+fn a_tighter_context_execution_bound_keeps_the_plan() {
+    let fixture = through_stages(
+        &[
+            query("sum by (job) (sum_over_time(m[10m]))", 20_000),
+            query("sum(sum_over_time(m[10m]))", 30_000),
+        ],
+        true,
+    );
+    let family = "V (context execution at a rate component)";
+    let tight = big_m_per_row(&fixture, family, false);
+    let loose = big_m_per_row(&fixture, family, true);
+    assert_eq!(tight.len(), loose.len());
+    assert!(tight
+        .iter()
+        .zip(&loose)
+        .all(|(t, l)| *t <= l * (1.0 + 1e-12)));
+    assert!(tight.iter().zip(&loose).any(|(t, l)| *t < l * (1.0 - 1e-9)));
+
+    let tightened = solve_and_check(&fixture, &aws_fargate_cost_inputs(), |_| {}).unwrap();
+    LOOSE_BIG_M.with(|cell| cell.set(true));
+    let looser = solve_and_check(&fixture, &aws_fargate_cost_inputs(), |_| {});
+    LOOSE_BIG_M.with(|cell| cell.set(false));
+    let looser = looser.unwrap();
+    assert_eq!(tightened.solution.status, SelectionMilpStatus::Optimal);
+    assert_eq!(looser.solution.status, SelectionMilpStatus::Optimal);
+    let tolerance = 1e-9 * looser.solution.cost_per_second.abs();
+    assert!(
+        (tightened.solution.cost_per_second - looser.solution.cost_per_second).abs() <= tolerance
+    );
+    assert_eq!(
+        tightened.dag.candidate_of_node,
+        looser.dag.candidate_of_node
+    );
+}