Fix #186: Add MinimumMultiwayCut -> QUBO reduction

docs/paper/reductions.typ

@@ -3657,6 +3657,54 @@ where $P$ is a penalty weight large enough that any constraint violation costs m
   _Solution extraction._ Discard slack variables: return $bold(z)[0..n]$.
 ]
 
+#let mwc_qubo = load-example("MinimumMultiwayCut", "QUBO")
+#let mwc_qubo_sol = mwc_qubo.solutions.at(0)
+#let mwc_qubo_edges = mwc_qubo.source.instance.graph.inner.edges.map(e => (e.at(0), e.at(1)))
+#let mwc_qubo_weights = mwc_qubo.source.instance.edge_weights
+#let mwc_qubo_terminals = mwc_qubo.source.instance.terminals
+#let mwc_qubo_n = mwc_qubo.source.instance.graph.inner.nodes.len()
+#let mwc_qubo_k = mwc_qubo_terminals.len()
+#let mwc_qubo_nq = mwc_qubo_n * mwc_qubo_k
+#let mwc_qubo_alpha = mwc_qubo_weights.fold(0, (a, w) => a + w) + 1
+#let mwc_qubo_cut_indices = mwc_qubo_sol.source_config.enumerate().filter(((i, v)) => v == 1).map(((i, _)) => i)
+#let mwc_qubo_cut_cost = mwc_qubo_cut_indices.fold(0, (a, i) => a + mwc_qubo_weights.at(i))
+#reduction-rule("MinimumMultiwayCut", "QUBO",
+  example: true,
+  example-caption: [$n = #mwc_qubo_n$ vertices, $k = #mwc_qubo_k$ terminals $T = {#mwc_qubo_terminals.map(str).join(", ")}$, $|E| = #mwc_qubo_edges.len()$ edges],
+  extra: [
+    *Step 1 -- Source instance.* The canonical graph has $n = #mwc_qubo_n$ vertices, $m = #mwc_qubo_edges.len()$ edges with weights $(#mwc_qubo_weights.map(str).join(", "))$, and $k = #mwc_qubo_k$ terminals $T = {#mwc_qubo_terminals.map(str).join(", ")}$.
+
+    *Step 2 -- Introduce binary variables.* Assign $k = #mwc_qubo_k$ indicator variables per vertex: $x_(u,t) = 1$ means vertex $u$ belongs to terminal $t$'s component. This gives $n k = #mwc_qubo_n times #mwc_qubo_k = #mwc_qubo_nq$ QUBO variables:
+    $ underbrace(x_(0,0) x_(0,1) x_(0,2), "vertex 0") #h(4pt) underbrace(x_(1,0) x_(1,1) x_(1,2), "vertex 1") #h(4pt) dots.c #h(4pt) underbrace(x_(4,0) x_(4,1) x_(4,2), "vertex 4") $
+
+    *Step 3 -- Penalty coefficient.* $alpha = 1 + sum_(e in E) w(e) = 1 + #mwc_qubo_weights.map(str).join(" + ") = #mwc_qubo_alpha$.
+
+    *Step 4 -- Build $H_A$ (constraints).* One-hot: diagonal entries $Q_(u k+t, u k+t) = -#mwc_qubo_alpha$, off-diagonal $Q_(u k+s, u k+t) = #(2 * mwc_qubo_alpha)$ within each vertex's group. Terminal pinning: for each terminal vertex $t_i$, the wrong-position diagonal entries $Q_(t_i k+s, t_i k+s) += #mwc_qubo_alpha$ for $s != i$, effectively canceling the one-hot incentive for those positions.\
+
+    *Step 5 -- Build $H_B$ (cut cost).* For each edge $(u,v)$ with weight $w$ and each pair $s != t$, add $w$ to $Q_(u k+s, v k+t)$. For example, edge $(0,1)$ with weight $2$ contributes $2$ to positions $(x_(0,0), x_(1,1))$, $(x_(0,0), x_(1,2))$, $(x_(0,1), x_(1,0))$, $(x_(0,1), x_(1,2))$, $(x_(0,2), x_(1,0))$, and $(x_(0,2), x_(1,1))$.\
+
+    *Step 6 -- Verify a solution.* The QUBO ground state $bold(x) = (#mwc_qubo_sol.target_config.map(str).join(", "))$ decodes to the partition: vertex 0 in component 0, vertices 1--3 in component 1, vertex 4 in component 2. Cut edges: $\{#mwc_qubo_cut_indices.map(i => "(" + str(mwc_qubo_edges.at(i).at(0)) + "," + str(mwc_qubo_edges.at(i).at(1)) + ")").join(", ")\}$ with total weight #mwc_qubo_cut_indices.map(i => str(mwc_qubo_weights.at(i))).join(" + ") $= #mwc_qubo_cut_cost$ #sym.checkmark.
+  ],
+)[
+  The multiway cut problem requires a partition of vertices into $k$ components — one per terminal — minimizing the total weight of edges crossing components. The penalty method (@sec:penalty-method) encodes two constraints as QUBO penalties: (1) each vertex belongs to exactly one component (one-hot), and (2) each terminal is pinned to its own component. The cut-cost Hamiltonian counts edge weight across distinct components. Reference: @Heidari2022.
+][
+  _Construction._ Given $G = (V, E)$ with $n = |V|$, edge weights $w: E -> RR_(>0)$, and $k$ terminals $T = {t_0, ..., t_(k-1)}$. Introduce $n k$ binary variables $x_(u,t) in {0,1}$ (indexed by $u dot k + t$), where $x_(u,t) = 1$ means vertex $u$ is in terminal $t$'s component. Let $alpha = 1 + sum_(e in E) w(e)$.
+
+  The QUBO Hamiltonian is $H = H_A + H_B$ where:
+  $ H_A = alpha (sum_(u in V) (1 - sum_(t=0)^(k-1) x_(u,t))^2 + sum_(i=0)^(k-1) sum_(s != i) x_(t_i, s)) $
+  The first term is a _one-hot constraint_ ensuring each vertex is assigned to exactly one component. The second term _pins_ each terminal $t_i$ to position $i$ by penalizing any other assignment. Expanding the one-hot term using $x^2 = x$:
+  $ Q_(u k+t, u k+t) = -alpha, quad Q_(u k+s, u k+t) = 2 alpha quad (s < t) $
+  Terminal pinning adds $alpha$ to the diagonal $Q_(t_i k+s, t_i k+s)$ for $s != i$, canceling the one-hot incentive.
+
+  The cut-cost Hamiltonian:
+  $ H_B = sum_((u,v) in E) sum_(s != t) w(u,v) dot x_(u,s) dot x_(v,t) $
+  counts the total weight of edges whose endpoints lie in different components.
+
+  _Correctness._ ($arrow.r.double$) A valid multiway cut with cost $C$ maps to a QUBO solution with $H_A = 0$ (valid partition with correct terminal pinning) and $H_B = C$. ($arrow.l.double$) If $H_A > 0$, the penalty $alpha > sum_e w(e)$ exceeds the entire cut-cost range, so any QUBO minimizer has $H_A = 0$, encoding a valid partition. Among valid partitions, $H_B$ equals the cut cost, and the minimizer achieves the minimum multiway cut.
+
+  _Solution extraction._ For each vertex $u$, find terminal position $t$ with $x_(u,t) = 1$. For each edge $(u,v)$, output 1 (cut) if $u$ and $v$ are in different components, 0 otherwise.
+]
+
 #let qubo_ilp = load-example("QUBO", "ILP")
 #let qubo_ilp_sol = qubo_ilp.solutions.at(0)
 #reduction-rule("QUBO", "ILP",

docs/paper/references.bib

@@ -798,6 +798,14 @@ @article{papadimitriou1982
   doi     = {10.1145/322307.322309}
 }
 
+@techreport{Heidari2022,
+  author      = {Heidari, Shahrokh and Dinneen, Michael J. and Delmas, Patrice},
+  title       = {An Equivalent {QUBO} Model to the Minimum Multi-Way Cut Problem},
+  institution = {Centre for Discrete Mathematics and Theoretical Computer Science, University of Auckland},
+  year        = {2022},
+  number      = {CDMTCS-565}
+}
+
 @phdthesis{booth1975,
   author  = {Booth, Kellogg S.},
   title   = {{PQ}-Tree Algorithms},

src/rules/minimummultiwaycut_qubo.rs

@@ -0,0 +1,163 @@
+//! Reduction from MinimumMultiwayCut to QUBO.
+//!
+//! Variable mapping: k*n binary variables x_{u,t} for each vertex u and
+//! terminal position t. x_{u,t} = 1 means vertex u is assigned to terminal t's
+//! component. Variable index: u * k + t.
+//!
+//! QUBO Hamiltonian: H = H_A + H_B
+//!
+//! H_A enforces valid partition (one-hot per vertex) and terminal pinning.
+//! H_B encodes the cut cost objective.
+//!
+//! Reference: Heidari, Dinneen & Delmas (2022).
+
+use crate::models::algebraic::QUBO;
+use crate::models::graph::MinimumMultiwayCut;
+use crate::reduction;
+use crate::rules::traits::{ReduceTo, ReductionResult};
+use crate::topology::{Graph, SimpleGraph};
+
+/// Result of reducing MinimumMultiwayCut to QUBO.
+#[derive(Debug, Clone)]
+pub struct ReductionMinimumMultiwayCutToQUBO {
+    target: QUBO<f64>,
+    num_vertices: usize,
+    num_terminals: usize,
+    edges: Vec<(usize, usize)>,
+}
+
+impl ReductionResult for ReductionMinimumMultiwayCutToQUBO {
+    type Source = MinimumMultiwayCut<SimpleGraph, i32>;
+    type Target = QUBO<f64>;
+
+    fn target_problem(&self) -> &Self::Target {
+        &self.target
+    }
+
+    /// Decode one-hot assignment: for each vertex find its terminal, then
+    /// for each edge check if endpoints are in different terminals.
+    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
+        let k = self.num_terminals;
+        let n = self.num_vertices;
+
+        // For each vertex, find which terminal position it is assigned to
+        let assignments: Vec<usize> = (0..n)
+            .map(|u| {
+                (0..k)
+                    .find(|&t| target_solution[u * k + t] == 1)
+                    .unwrap_or(0)
+            })
+            .collect();
+
+        // For each edge, output 1 (cut) if endpoints differ, 0 (keep) otherwise
+        self.edges
+            .iter()
+            .map(|&(u, v)| {
+                if assignments[u] != assignments[v] {
+                    1
+                } else {
+                    0
+                }
+            })
+            .collect()
+    }
+}
+
+#[reduction(overhead = { num_vars = "num_terminals * num_vertices" })]
+impl ReduceTo<QUBO<f64>> for MinimumMultiwayCut<SimpleGraph, i32> {
+    type Result = ReductionMinimumMultiwayCutToQUBO;
+
+    fn reduce_to(&self) -> Self::Result {
+        let n = self.num_vertices();
+        let k = self.num_terminals();
+        let edges = self.graph().edges();
+        let edge_weights = self.edge_weights();
+        let terminals = self.terminals();
+        let nq = n * k;
+
+        // Penalty: sum of all edge weights + 1
+        let alpha: f64 = edge_weights.iter().map(|&w| (w as f64).abs()).sum::<f64>() + 1.0;
+
+        let mut matrix = vec![vec![0.0f64; nq]; nq];
+
+        // Helper: add value to upper-triangular position
+        let mut add_upper = |i: usize, j: usize, val: f64| {
+            let (lo, hi) = if i <= j { (i, j) } else { (j, i) };
+            matrix[lo][hi] += val;
+        };
+
+        // H_A: one-hot constraint per vertex
+        // (1 - sum_t x_{u,t})^2 = 1 - sum_t x_{u,t} + 2 * sum_{s<t} x_{u,s} * x_{u,t}
+        // (using x^2 = x for binary variables)
+        for u in 0..n {
+            // Diagonal: -alpha for each terminal position
+            for s in 0..k {
+                add_upper(u * k + s, u * k + s, -alpha);
+            }
+            // Off-diagonal within same vertex: +2*alpha for each pair
+            for s in 0..k {
+                for t in (s + 1)..k {
+                    add_upper(u * k + s, u * k + t, 2.0 * alpha);
+                }
+            }
+        }
+
+        // H_A: terminal pinning
+        // For each terminal vertex, penalize assignment to wrong position
+        for (t_pos, &t_vertex) in terminals.iter().enumerate() {
+            for s in 0..k {
+                if s != t_pos {
+                    add_upper(t_vertex * k + s, t_vertex * k + s, alpha);
+                }
+            }
+        }
+
+        // H_B: cut cost
+        // For each edge (u,v) with weight w, for each pair of distinct
+        // terminal positions s != t: add w to Q[u*k+s, v*k+t]
+        for (edge_idx, &(u, v)) in edges.iter().enumerate() {
+            let w = edge_weights[edge_idx] as f64;
+            for s in 0..k {
+                for t in 0..k {
+                    if s != t {
+                        add_upper(u * k + s, v * k + t, w);
+                    }
+                }
+            }
+        }
+
+        ReductionMinimumMultiwayCutToQUBO {
+            target: QUBO::from_matrix(matrix),
+            num_vertices: n,
+            num_terminals: k,
+            edges,
+        }
+    }
+}
+
+#[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: "minimummultiwaycut_to_qubo",
+        build: || {
+            use crate::models::algebraic::QUBO;
+            use crate::models::graph::MinimumMultiwayCut;
+            use crate::topology::SimpleGraph;
+            let graph = SimpleGraph::new(5, vec![(0, 1), (1, 2), (2, 3), (3, 4), (0, 4), (1, 3)]);
+            let source = MinimumMultiwayCut::new(graph, vec![0, 2, 4], vec![2, 3, 1, 2, 4, 5]);
+            crate::example_db::specs::rule_example_with_witness::<_, QUBO<f64>>(
+                source,
+                SolutionPair {
+                    source_config: vec![1, 0, 0, 1, 1, 0],
+                    target_config: vec![1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 1],
+                },
+            )
+        },
+    }]
+}
+
+#[cfg(test)]
+#[path = "../unit_tests/rules/minimummultiwaycut_qubo.rs"]
+mod tests;

src/rules/mod.rs

@@ -24,6 +24,7 @@ mod maximumindependentset_triangular;
 pub(crate) mod maximummatching_maximumsetpacking;
 mod maximumsetpacking_casts;
 pub(crate) mod maximumsetpacking_qubo;
+pub(crate) mod minimummultiwaycut_qubo;
 pub(crate) mod minimumvertexcover_maximumindependentset;
 pub(crate) mod minimumvertexcover_minimumsetcovering;
 pub(crate) mod sat_circuitsat;
@@ -94,6 +95,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::Ru
     specs.extend(maximumindependentset_maximumsetpacking::canonical_rule_example_specs());
     specs.extend(maximummatching_maximumsetpacking::canonical_rule_example_specs());
     specs.extend(maximumsetpacking_qubo::canonical_rule_example_specs());
+    specs.extend(minimummultiwaycut_qubo::canonical_rule_example_specs());
     specs.extend(minimumvertexcover_maximumindependentset::canonical_rule_example_specs());
     specs.extend(minimumvertexcover_minimumsetcovering::canonical_rule_example_specs());
     specs.extend(sat_circuitsat::canonical_rule_example_specs());

src/unit_tests/rules/minimummultiwaycut_qubo.rs

@@ -0,0 +1,100 @@
+use super::*;
+use crate::solvers::BruteForce;
+use crate::traits::Problem;
+use crate::types::SolutionSize;
+
+#[test]
+fn test_minimummultiwaycut_to_qubo_closed_loop() {
+    // 5 vertices, terminals {0,2,4}, 6 edges with weights [2,3,1,2,4,5]
+    let graph = SimpleGraph::new(5, vec![(0, 1), (1, 2), (2, 3), (3, 4), (0, 4), (1, 3)]);
+    let source = MinimumMultiwayCut::new(graph, vec![0, 2, 4], vec![2, 3, 1, 2, 4, 5]);
+
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&source);
+    let qubo = reduction.target_problem();
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_all_best(qubo);
+
+    assert!(!qubo_solutions.is_empty(), "QUBO solver found no solutions");
+
+    // All QUBO optimal solutions should extract to valid source solutions with cost 8
+    for sol in &qubo_solutions {
+        let extracted = reduction.extract_solution(sol);
+        let metric = source.evaluate(&extracted);
+        assert_eq!(metric, SolutionSize::Valid(8));
+    }
+}
+
+#[test]
+fn test_minimummultiwaycut_to_qubo_small() {
+    // 3 vertices, 2 terminals {0,2}, edges [(0,1),(1,2)] with weights [1,1]
+    let graph = SimpleGraph::new(3, vec![(0, 1), (1, 2)]);
+    let source = MinimumMultiwayCut::new(graph, vec![0, 2], vec![1, 1]);
+
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&source);
+    let qubo = reduction.target_problem();
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_all_best(qubo);
+
+    assert!(!qubo_solutions.is_empty(), "QUBO solver found no solutions");
+
+    // All solutions should extract to valid cuts
+    for sol in &qubo_solutions {
+        let extracted = reduction.extract_solution(sol);
+        let metric = source.evaluate(&extracted);
+        // With 2 terminals and path 0-1-2, minimum cut is 1 (cut either edge)
+        assert_eq!(metric, SolutionSize::Valid(1));
+    }
+}
+
+#[test]
+fn test_minimummultiwaycut_to_qubo_sizes() {
+    // 5 vertices, 3 terminals => QUBO has k*n = 3*5 = 15 variables
+    let graph = SimpleGraph::new(5, vec![(0, 1), (1, 2), (2, 3), (3, 4), (0, 4), (1, 3)]);
+    let source = MinimumMultiwayCut::new(graph, vec![0, 2, 4], vec![2, 3, 1, 2, 4, 5]);
+
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&source);
+    assert_eq!(reduction.target_problem().num_variables(), 15);
+}
+
+#[test]
+fn test_minimummultiwaycut_to_qubo_terminal_pinning() {
+    // Verify that in all QUBO optimal solutions, each terminal vertex is
+    // assigned to its own terminal position.
+    let graph = SimpleGraph::new(5, vec![(0, 1), (1, 2), (2, 3), (3, 4), (0, 4), (1, 3)]);
+    let terminals = vec![0, 2, 4];
+    let source = MinimumMultiwayCut::new(graph, terminals.clone(), vec![2, 3, 1, 2, 4, 5]);
+
+    let reduction = ReduceTo::<QUBO<f64>>::reduce_to(&source);
+    let qubo = reduction.target_problem();
+
+    let solver = BruteForce::new();
+    let qubo_solutions = solver.find_all_best(qubo);
+
+    let k = terminals.len();
+    for sol in &qubo_solutions {
+        for (t_pos, &t_vertex) in terminals.iter().enumerate() {
+            // Terminal vertex should be assigned to its own position
+            assert_eq!(
+                sol[t_vertex * k + t_pos],
+                1,
+                "Terminal {} at position {} should be 1",
+                t_vertex,
+                t_pos
+            );
+            // And not assigned to any other position
+            for s in 0..k {
+                if s != t_pos {
+                    assert_eq!(
+                        sol[t_vertex * k + s],
+                        0,
+                        "Terminal {} at position {} should be 0",
+                        t_vertex,
+                        s
+                    );
+                }
+            }
+        }
+    }
+}