diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index 79b4fbdea..a412a3a29 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -5330,6 +5330,53 @@ The following reductions to Integer Linear Programming are straightforward formu _Solution extraction._ For each edge $e$ at index $"idx"$, read $x_e = x^*_(k n + "idx")$. The source configuration is $"config"[e] = x_e$ (1 = cut, 0 = keep). ] +#let st_ilp = load-example("SteinerTree", "ILP") +#let st_ilp_sol = st_ilp.solutions.at(0) +#let st_edges = st_ilp.source.instance.graph.edges +#let st_weights = st_ilp.source.instance.edge_weights +#let st_terminals = st_ilp.source.instance.terminals +#let st_root = st_terminals.at(0) +#let st_non_root_terminals = range(1, st_terminals.len()).map(i => st_terminals.at(i)) +#let st_selected_edge_indices = st_ilp_sol.source_config.enumerate().filter(((i, v)) => v == 1).map(((i, _)) => i) +#let st_selected_edges = st_selected_edge_indices.map(i => st_edges.at(i)) +#let st_cost = st_selected_edge_indices.map(i => st_weights.at(i)).sum() + +#reduction-rule("SteinerTree", "ILP", + example: true, + example-caption: [Canonical Steiner tree instance ($n = #st_ilp.source.instance.graph.num_vertices$, $m = #st_edges.len()$, $|T| = #st_terminals.len()$)], + extra: [ + *Step 1 -- Choose a root and one commodity per remaining terminal.* The canonical source instance has terminals $T = {#st_terminals.map(t => $v_#t$).join(", ")}$. The reduction fixes the first terminal as root $r = v_#st_root$ and creates one flow commodity for each remaining terminal: $v_#st_non_root_terminals.at(0)$ and $v_#st_non_root_terminals.at(1)$. + + *Step 2 -- Count the variables from the source edge order.* The first #st_edges.len() target variables are the edge selectors $bold(y) = (#st_ilp_sol.target_config.slice(0, st_edges.len()).map(str).join(", "))$, one per source edge in the order #st_edges.enumerate().map(((i, e)) => [$e_#i = (#(e.at(0)), #(e.at(1)))$]).join(", "). The remaining #(st_ilp.target.instance.num_vars - st_edges.len()) variables are directed flow indicators: $2 m (|T| - 1) = 2 times #st_edges.len() times #st_non_root_terminals.len() = #(st_ilp.target.instance.num_vars - st_edges.len())$. + + *Step 3 -- Count the constraints commodity-by-commodity.* Each non-root terminal contributes one flow-conservation equality per vertex and two capacity inequalities per source edge. For this fixture that is $#st_ilp.source.instance.graph.num_vertices times #st_non_root_terminals.len() = #(st_ilp.source.instance.graph.num_vertices * st_non_root_terminals.len())$ equalities plus $#(2 * st_edges.len()) times #st_non_root_terminals.len() = #(2 * st_edges.len() * st_non_root_terminals.len())$ inequalities, totaling #st_ilp.target.instance.constraints.len() constraints. + + *Step 4 -- Read the canonical witness pair.* The source witness selects edges ${#st_selected_edges.map(e => $(v_#(e.at(0)), v_#(e.at(1)))$).join(", ")}$, so $bold(y)$ already encodes the Steiner tree. In the target witness, the commodity for $v_2$ routes along $v_0 arrow v_1 arrow v_2$, while the commodity for $v_4$ routes along $v_0 arrow v_1 arrow v_3 arrow v_4$. Every flow 1-entry therefore sits under a selected edge variable #sym.checkmark + + *Step 5 -- Verify the objective end-to-end.* The selected-edge prefix is $bold(y) = (#st_ilp_sol.target_config.slice(0, st_edges.len()).map(str).join(", "))$, matching the source witness $(#st_ilp_sol.source_config.map(str).join(", "))$. The ILP objective is #st_selected_edge_indices.map(i => $#(st_weights.at(i))$).join($+$) $= #st_cost$, exactly the Steiner tree optimum stored in the fixture. + + *Multiplicity:* The fixture stores one canonical witness. Other optimal Steiner trees could yield different feasible ILP witnesses, but every valid witness still exposes the source solution in the first $m$ variables. + ], +)[ + The rooted multi-commodity flow formulation @wong1984steiner @kochmartin1998steiner introduces one binary selector $y_e$ for each source edge and, for every non-root terminal $t$, one binary flow variable on each directed source edge. Flow conservation sends one unit from the root to each terminal, while the linking inequalities $f^t_(u,v) <= y_e$ ensure that every used flow arc is backed by a selected source edge. The resulting binary ILP has $m + 2 m (k - 1)$ variables and $n (k - 1) + 2 m (k - 1)$ constraints. +][ + _Construction._ Given an undirected weighted graph $G = (V, E, w)$ with strictly positive edge weights, terminals $T = {t_0, dots, t_(k-1)}$, and root $r = t_0$, introduce binary edge selectors $y_e in {0,1}$ for every $e in E$. For each non-root terminal $t in T backslash {r}$ and each directed copy of an undirected edge $(u, v) in E$, introduce a binary flow variable $f^t_(u,v) in {0,1}$. The target objective is + $ min sum_(e in E) w_e y_e. $ + For every commodity $t$ and vertex $v$, enforce flow conservation: + $ sum_(u : (u, v) in A) f^t_(u,v) - sum_(u : (v, u) in A) f^t_(v,u) = b_(t,v), $ + where $A$ contains both orientations of every undirected edge, $b_(t,v) = -1$ at the root $v = r$, $b_(t,v) = 1$ at the sink $v = t$, and $b_(t,v) = 0$ otherwise. For every commodity $t$ and undirected edge $e = {u, v}$, add the capacity-linking inequalities + $ f^t_(u,v) <= y_e quad "and" quad f^t_(v,u) <= y_e. $ + Binary flow variables suffice because any Steiner tree yields a unique simple root-to-terminal path for each commodity, so every commodity can be realized as a 0/1 path indicator. + + _Correctness._ ($arrow.r.double$) If $S subset.eq E$ is a Steiner tree, set $y_e = 1$ exactly for $e in S$. For each non-root terminal $t$, the unique path from $r$ to $t$ inside the tree defines a binary flow assignment satisfying the conservation equations, and every used arc lies on a selected edge, so all linking inequalities hold. The ILP objective equals $sum_(e in S) w_e$. ($arrow.l.double$) Any feasible ILP solution with edge selector set $Y = {e in E : y_e = 1}$ supports one unit of flow from $r$ to every non-root terminal, so the selected edges contain a connected subgraph spanning all terminals. Because all edge weights are strictly positive, any cycle in the selected subgraph has positive total cost; the optimizer therefore never includes redundant edges, so the selected subgraph is already a Steiner tree. Therefore an optimal ILP solution induces a minimum-cost Steiner tree. + + _Variable mapping._ The first $m$ ILP variables are the source-edge indicators $y_0, dots, y_(m-1)$ in source edge order. For terminal $t_p$ with $p in {1, dots, k-1}$, the next block of $2 m$ variables stores the directed arc indicators $f^(t_p)_(u,v)$ and $f^(t_p)_(v,u)$ for each source edge $(u, v)$. + + _Solution extraction._ Read the first $m$ target variables as the source edge-selection vector. Since those coordinates are exactly the $y_e$ variables, the extracted source configuration is valid whenever the selected subgraph is pruned to its Steiner tree witness. + + _Remark._ Zero-weight edges are excluded because they allow degenerate optimal ILP solutions containing redundant cycles at no cost; following the convention of practical solvers (e.g., SCIP-Jack @kochmartin1998steiner), such edges should be contracted before applying the reduction. +] + == Unit Disk Mapping #reduction-rule("MaximumIndependentSet", "KingsSubgraph")[ diff --git a/docs/paper/references.bib b/docs/paper/references.bib index 9fe022a25..dc2a63d5b 100644 --- a/docs/paper/references.bib +++ b/docs/paper/references.bib @@ -1188,6 +1188,28 @@ @article{schaefer1978 doi = {10.1145/800133.804350} } +@article{wong1984steiner, + author = {R. T. Wong}, + title = {A Dual Ascent Approach for Steiner Tree Problems on a Directed Graph}, + journal = {Mathematical Programming}, + volume = {28}, + number = {3}, + pages = {271--287}, + year = {1984}, + doi = {10.1007/BF02612335} +} + +@article{kochmartin1998steiner, + author = {Thorsten Koch and Alexander Martin}, + title = {Solving Steiner Tree Problems in Graphs to Optimality}, + journal = {Networks}, + volume = {32}, + number = {3}, + pages = {207--232}, + year = {1998}, + doi = {10.1002/(SICI)1097-0037(199810)32:3<207::AID-NET5>3.0.CO;2-O} +} + @inproceedings{stockmeyer1973, author = {Larry J. Stockmeyer and Albert R. Meyer}, title = {Word Problems Requiring Exponential Time: Preliminary Report}, diff --git a/src/rules/mod.rs b/src/rules/mod.rs index 90cf69544..4b55c9aae 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -82,6 +82,8 @@ pub(crate) mod qubo_ilp; #[cfg(feature = "ilp-solver")] pub(crate) mod sequencingtominimizeweightedcompletiontime_ilp; #[cfg(feature = "ilp-solver")] +pub(crate) mod steinertree_ilp; +#[cfg(feature = "ilp-solver")] pub(crate) mod travelingsalesman_ilp; pub use graph::{ @@ -138,6 +140,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec, + num_edges: usize, +} + +impl ReductionResult for ReductionSteinerTreeToILP { + type Source = SteinerTree; + type Target = ILP; + + fn target_problem(&self) -> &ILP { + &self.target + } + + fn extract_solution(&self, target_solution: &[usize]) -> Vec { + target_solution[..self.num_edges].to_vec() + } +} + +#[reduction( + overhead = { + num_vars = "num_edges + 2 * num_edges * (num_terminals - 1)", + num_constraints = "num_vertices * (num_terminals - 1) + 2 * num_edges * (num_terminals - 1)", + } +)] +impl ReduceTo> for SteinerTree { + type Result = ReductionSteinerTreeToILP; + + fn reduce_to(&self) -> Self::Result { + assert!( + self.edge_weights().iter().all(|&weight| weight > 0), + "SteinerTree -> ILP requires strictly positive edge weights (zero-weight edges should be contracted beforehand)" + ); + + let n = self.num_vertices(); + let m = self.num_edges(); + let root = self.terminals()[0]; + let non_root_terminals = &self.terminals()[1..]; + let edges = self.graph().edges(); + let num_vars = m + 2 * m * non_root_terminals.len(); + let num_constraints = n * non_root_terminals.len() + 2 * m * non_root_terminals.len(); + let mut constraints = Vec::with_capacity(num_constraints); + + let edge_var = |edge_idx: usize| edge_idx; + let flow_var = |terminal_pos: usize, edge_idx: usize, dir: usize| -> usize { + m + terminal_pos * 2 * m + 2 * edge_idx + dir + }; + + for (terminal_pos, &terminal) in non_root_terminals.iter().enumerate() { + for vertex in 0..n { + let mut terms = Vec::new(); + for (edge_idx, &(u, v)) in edges.iter().enumerate() { + if v == vertex { + terms.push((flow_var(terminal_pos, edge_idx, 0), 1.0)); + terms.push((flow_var(terminal_pos, edge_idx, 1), -1.0)); + } + if u == vertex { + terms.push((flow_var(terminal_pos, edge_idx, 0), -1.0)); + terms.push((flow_var(terminal_pos, edge_idx, 1), 1.0)); + } + } + + let rhs = if vertex == root { + -1.0 + } else if vertex == terminal { + 1.0 + } else { + 0.0 + }; + constraints.push(LinearConstraint::eq(terms, rhs)); + } + } + + for terminal_pos in 0..non_root_terminals.len() { + for edge_idx in 0..m { + let selector = edge_var(edge_idx); + constraints.push(LinearConstraint::le( + vec![(flow_var(terminal_pos, edge_idx, 0), 1.0), (selector, -1.0)], + 0.0, + )); + constraints.push(LinearConstraint::le( + vec![(flow_var(terminal_pos, edge_idx, 1), 1.0), (selector, -1.0)], + 0.0, + )); + } + } + + let objective: Vec<(usize, f64)> = self + .edge_weights() + .iter() + .enumerate() + .map(|(edge_idx, &weight)| (edge_var(edge_idx), weight as f64)) + .collect(); + + let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize); + + ReductionSteinerTreeToILP { + target, + num_edges: m, + } + } +} + +#[cfg(feature = "example-db")] +pub(crate) fn canonical_rule_example_specs() -> Vec { + use crate::export::SolutionPair; + + vec![crate::example_db::specs::RuleExampleSpec { + id: "steinertree_to_ilp", + build: || { + let source = SteinerTree::new( + SimpleGraph::new( + 5, + vec![(0, 1), (1, 2), (1, 3), (3, 4), (0, 3), (3, 2), (2, 4)], + ), + vec![2, 2, 1, 1, 5, 5, 6], + vec![0, 2, 4], + ); + crate::example_db::specs::rule_example_with_witness::<_, ILP>( + source, + SolutionPair { + source_config: vec![1, 1, 1, 1, 0, 0, 0], + target_config: vec![ + 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, + 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, + ], + }, + ) + }, + }] +} + +#[cfg(test)] +#[path = "../unit_tests/rules/steinertree_ilp.rs"] +mod tests; diff --git a/src/unit_tests/rules/steinertree_ilp.rs b/src/unit_tests/rules/steinertree_ilp.rs new file mode 100644 index 000000000..25420eaa9 --- /dev/null +++ b/src/unit_tests/rules/steinertree_ilp.rs @@ -0,0 +1,97 @@ +use super::*; +use crate::models::algebraic::{ObjectiveSense, ILP}; +use crate::models::graph::SteinerTree; +use crate::rules::ReduceTo; +use crate::solvers::{BruteForce, ILPSolver}; +use crate::topology::SimpleGraph; +use crate::traits::Problem; +use crate::types::SolutionSize; + +fn canonical_instance() -> SteinerTree { + let graph = SimpleGraph::new( + 5, + vec![(0, 1), (1, 2), (1, 3), (3, 4), (0, 3), (3, 2), (2, 4)], + ); + SteinerTree::new(graph, vec![2, 2, 1, 1, 5, 5, 6], vec![0, 2, 4]) +} + +#[test] +fn test_reduction_creates_expected_ilp_shape() { + let problem = canonical_instance(); + let reduction: ReductionSteinerTreeToILP = ReduceTo::>::reduce_to(&problem); + let ilp = reduction.target_problem(); + + assert_eq!(ilp.num_vars, 35); + assert_eq!(ilp.constraints.len(), 38); + assert_eq!(ilp.sense, ObjectiveSense::Minimize); + assert_eq!( + ilp.objective, + vec![ + (0, 2.0), + (1, 2.0), + (2, 1.0), + (3, 1.0), + (4, 5.0), + (5, 5.0), + (6, 6.0), + ] + ); +} + +#[test] +fn test_steinertree_to_ilp_closed_loop() { + let problem = canonical_instance(); + let reduction: ReductionSteinerTreeToILP = ReduceTo::>::reduce_to(&problem); + let ilp = reduction.target_problem(); + + let bf = BruteForce::new(); + let ilp_solver = ILPSolver::new(); + let best_source = bf.find_all_best(&problem); + let ilp_solution = ilp_solver.solve(ilp).expect("ILP should be solvable"); + let extracted = reduction.extract_solution(&ilp_solution); + + assert_eq!(problem.evaluate(&best_source[0]), SolutionSize::Valid(6)); + assert_eq!(problem.evaluate(&extracted), SolutionSize::Valid(6)); + assert!(problem.is_valid_solution(&extracted)); +} + +#[test] +fn test_solution_extraction_reads_edge_selector_prefix() { + let problem = canonical_instance(); + let reduction: ReductionSteinerTreeToILP = ReduceTo::>::reduce_to(&problem); + + let target_solution = vec![ + 1, 1, 1, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, + 0, 0, 0, 0, 0, + ]; + + assert_eq!( + reduction.extract_solution(&target_solution), + vec![1, 1, 1, 1, 0, 0, 0] + ); +} + +#[test] +fn test_solve_reduced_uses_new_rule() { + let problem = canonical_instance(); + let solution = ILPSolver::new() + .solve_reduced(&problem) + .expect("solve_reduced should find the Steiner tree via ILP"); + assert_eq!(problem.evaluate(&solution), SolutionSize::Valid(6)); +} + +#[test] +#[should_panic(expected = "SteinerTree -> ILP requires strictly positive edge weights")] +fn test_reduction_rejects_negative_weights() { + let graph = SimpleGraph::new(3, vec![(0, 1), (1, 2), (0, 2)]); + let problem = SteinerTree::new(graph, vec![1, -2, 3], vec![0, 1]); + let _ = ReduceTo::>::reduce_to(&problem); +} + +#[test] +#[should_panic(expected = "SteinerTree -> ILP requires strictly positive edge weights")] +fn test_reduction_rejects_zero_weights() { + let graph = SimpleGraph::new(3, vec![(0, 1), (1, 2), (0, 2)]); + let problem = SteinerTree::new(graph, vec![0, 0, 0], vec![0, 1]); + let _ = ReduceTo::>::reduce_to(&problem); +}