Fix #123: [Rule] SteinerTree to ILP

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")[

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},

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<crate::example_db::specs::Ru
         specs.extend(qubo_ilp::canonical_rule_example_specs());
         specs
             .extend(sequencingtominimizeweightedcompletiontime_ilp::canonical_rule_example_specs());
+        specs.extend(steinertree_ilp::canonical_rule_example_specs());
         specs.extend(travelingsalesman_ilp::canonical_rule_example_specs());
     }
     specs

src/rules/steinertree_ilp.rs

@@ -0,0 +1,156 @@
+//! Reduction from SteinerTree to ILP (Integer Linear Programming).
+//!
+//! Uses the standard rooted multi-commodity flow formulation:
+//! - Variables: edge selectors `y_e` plus directed flow variables `f^t_(u,v)`
+//!   for each non-root terminal `t`
+//! - Constraints: flow conservation for each commodity and capacity linking
+//!   `f^t_(u,v) <= y_e`
+//! - Objective: minimize the total weight of selected edges
+
+use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
+use crate::models::graph::SteinerTree;
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+use crate::topology::{Graph, SimpleGraph};
+
+/// Result of reducing SteinerTree to ILP.
+///
+/// Variable layout (all binary):
+/// - `y_e` for each undirected source edge `e` (indices `0..m`)
+/// - `f^t_(u,v)` and `f^t_(v,u)` for each non-root terminal `t` and each source edge
+///   `(u, v)` (indices `m..m + 2m(k-1)`)
+#[derive(Debug, Clone)]
+pub struct ReductionSteinerTreeToILP {
+    target: ILP<bool>,
+    num_edges: usize,
+}
+
+impl ReductionResult for ReductionSteinerTreeToILP {
+    type Source = SteinerTree<SimpleGraph, i32>;
+    type Target = ILP<bool>;
+
+    fn target_problem(&self) -> &ILP<bool> {
+        &self.target
+    }
+
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        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<ILP<bool>> for SteinerTree<SimpleGraph, i32> {
+    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<crate::example_db::specs::RuleExampleSpec> {
+    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<bool>>(
+                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;

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<SimpleGraph, i32> {
+    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::<ILP<bool>>::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::<ILP<bool>>::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::<ILP<bool>>::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::<ILP<bool>>::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::<ILP<bool>>::reduce_to(&problem);
+}