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..ac64258a5 100644
--- a/crates/plan-selection/Cargo.toml
+++ b/crates/plan-selection/Cargo.toml
@@ -25,6 +25,11 @@ 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; a diagnosed run's timings are not
+# comparable with an undiagnosed one's).
+lp-diagnosis = []
 
 [dev-dependencies]
 asap-types = { path = "../types", features = ["test-support"] }
diff --git a/crates/plan-selection/src/milp/cost_terms.rs b/crates/plan-selection/src/milp/cost_terms.rs
index 3be9a9ecb..dbc2483ab 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,25 @@ 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 = 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
new file mode 100644
index 000000000..c946d61b4
--- /dev/null
+++ b/crates/plan-selection/src/milp/lp_diagnosis.rs
@@ -0,0 +1,303 @@
+//! #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.
+//!
+//! 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::{objective_scale, SelectionMilp, HIGHS_THREADS};
+
+/// 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 LP_TIME_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 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();
+            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", 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);
+                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
+}
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 f148955be..ee987788f 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
@@ -122,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 {
@@ -133,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);
     }
@@ -262,6 +274,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 +318,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 +505,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);
@@ -1562,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;
@@ -1604,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,
@@ -1625,7 +1649,12 @@ 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 = if tests_loose_big_m() {
+            counted.len() as f64
+        } else {
+            places.len() as f64
+        };
         if counted.iter().all(|&(_, _, parts)| parts as f64 >= most) {
             continue;
         }
@@ -1654,8 +1683,155 @@ 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 part
+        // without a query when there is one (`None` sorts first), else in
+        // the first query's 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 +1848,43 @@ 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();
+    // 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
+            .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))
+    };
+    #[cfg(feature = "lp-diagnosis")]
+    let deadline = deadline.map(|deadline| deadline + diagnosis_started.elapsed());
+    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 +1895,163 @@ 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}")
     });
+    #[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,
-        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 +2078,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,8 +2086,120 @@ 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()),
+        });
+    }
+    // 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(),
+        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;
+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;
@@ -1836,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())
@@ -1856,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 861a11a2e..ba87120eb 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,132 @@ 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);
+}
+
+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
+    );
+}
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()