problemreductions/rules/eulerianpath_ilp.rs
1//! Reduction from EulerianPath to ILP (Integer Linear Programming).
2//!
3//! Encodes the directed Eulerian-trail witness structure as an ILP feasibility
4//! instance with integer variables. Given a directed multigraph
5//! `D = (V, A)` with arc occurrences `A = {a_1, ..., a_m}`:
6//!
7//! - For every compatible ordered pair `(a, b) in P` (where
8//! `P = { (a, b) : a != b and head(a) = tail(b) }`), introduce an integer
9//! successor variable `y_{a,b}` (intended `0/1`: `1` iff `b` immediately
10//! follows `a` in the trail).
11//! - For every arc `a in A`, introduce integer variables `s_a`, `e_a`
12//! (`0/1`: `s_a = 1` iff `a` is first, `e_a = 1` iff `a` is last) and an
13//! integer position variable `u_a` (intended value `0..m-1`).
14//! - The predecessor and successor equalities, together with the unique start
15//! and unique end constraints and Miller--Tucker--Zemlin-style ordering
16//! constraints, force any feasible solution to encode a directed trail that
17//! uses every arc occurrence exactly once.
18//!
19//! The empty-arc instance (`m = 0`) maps to the empty ILP with no variables
20//! and no constraints, which is vacuously feasible.
21//!
22//! References: Ebert, "Computing Eulerian trails," IPL 28(2):93--97 (1988);
23//! Bang-Jensen and Gutin, *Digraphs: Theory, Algorithms and Applications*,
24//! 2nd ed., Springer (2009).
25
26use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
27use crate::models::graph::EulerianPath;
28use crate::reduction;
29use crate::rules::traits::{ReduceTo, ReductionResult};
30
31/// Result of reducing EulerianPath to `ILP<i32>`.
32///
33/// Variable layout (all in the non-negative integer domain, with explicit
34/// upper bounds enforcing the intended `0/1` and `0..m-1` ranges):
35/// - `y_{a,b}` at index `k` for the `k`-th compatible pair in `pairs` order
36/// (a single sweep over `(a, b)` with `a, b in 0..m`, `a != b`, in
37/// row-major order),
38/// - `s_a` at index `p + a` for `a in 0..m`,
39/// - `e_a` at index `p + m + a` for `a in 0..m`,
40/// - `u_a` at index `p + 2 * m + a` for `a in 0..m`,
41///
42/// where `p = pairs.len()` is the number of compatible ordered pairs.
43#[derive(Debug, Clone)]
44pub struct ReductionEulerianPathToILP {
45 target: ILP<i32>,
46 /// Compatible ordered pairs `(a, b)` in the order their `y_{a,b}` variables
47 /// appear in the ILP, for `m > 0`. Empty when `m = 0`.
48 pairs: Vec<(usize, usize)>,
49 /// Number of arc occurrences in the source instance.
50 num_arcs: usize,
51}
52
53impl ReductionEulerianPathToILP {
54 fn s_idx(&self, a: usize) -> usize {
55 self.pairs.len() + a
56 }
57}
58
59impl ReductionResult for ReductionEulerianPathToILP {
60 type Source = EulerianPath;
61 type Target = ILP<i32>;
62
63 fn target_problem(&self) -> &ILP<i32> {
64 &self.target
65 }
66
67 /// Decode an ILP assignment into a source arc ordering.
68 ///
69 /// Reads the unique active start arc (`s_a = 1`) and walks the active
70 /// successor relation (`y_{a,b} = 1`) one step at a time, producing an arc
71 /// permutation of length `m`. If the assignment is malformed (no start,
72 /// no successor mid-walk, or revisits an arc) we fall back to the identity
73 /// ordering `0..m` in release builds; debug builds trip a
74 /// `debug_assert!` to surface the caller bug. Callers must independently
75 /// check feasibility on the source side via
76 /// `EulerianPath::is_valid_solution`.
77 fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
78 let m = self.num_arcs;
79 if m == 0 {
80 return Vec::new();
81 }
82 let fallback: Vec<usize> = (0..m).collect();
83
84 // Find the unique active start arc.
85 let mut current = match (0..m)
86 .find(|&a| target_solution.get(self.s_idx(a)).copied().unwrap_or(0) == 1)
87 {
88 Some(a) => a,
89 None => {
90 debug_assert!(
91 false,
92 "EulerianPath -> ILP extract_solution: malformed assignment, no active start arc (expected exactly one s_a = 1)",
93 );
94 return fallback;
95 }
96 };
97
98 // Walk the active successor relation, recording each visited arc.
99 let mut order = Vec::with_capacity(m);
100 let mut visited = vec![false; m];
101 order.push(current);
102 visited[current] = true;
103
104 for _ in 1..m {
105 let next = self
106 .pairs
107 .iter()
108 .enumerate()
109 .find(|&(k, &(a, _))| {
110 a == current && target_solution.get(k).copied().unwrap_or(0) == 1
111 })
112 .map(|(_, &(_, b))| b);
113
114 match next {
115 Some(b) if !visited[b] => {
116 order.push(b);
117 visited[b] = true;
118 current = b;
119 }
120 _ => {
121 debug_assert!(
122 false,
123 "EulerianPath -> ILP extract_solution: malformed assignment at arc {} (expected exactly one active successor y_{{{},b}} = 1 leading to an unvisited arc)",
124 current,
125 current,
126 );
127 return fallback;
128 }
129 }
130 }
131 order
132 }
133}
134
135/// Enumerate compatible ordered pairs `(a, b)` with `a != b` and
136/// `head(a) = tail(b)`. The order is `a`-major then `b`-major, matching the
137/// nested-loop construction below.
138fn compatible_pairs(arcs: &[(usize, usize)]) -> Vec<(usize, usize)> {
139 let mut pairs = Vec::new();
140 for (a, &(_, head_a)) in arcs.iter().enumerate() {
141 for (b, &(tail_b, _)) in arcs.iter().enumerate() {
142 if a != b && tail_b == head_a {
143 pairs.push((a, b));
144 }
145 }
146 }
147 pairs
148}
149
150#[reduction(
151 overhead = {
152 num_vars = "3 * num_arcs + num_arcs * num_arcs",
153 num_constraints = "5 * num_arcs + 2 * num_arcs * num_arcs + 2",
154 }
155)]
156impl ReduceTo<ILP<i32>> for EulerianPath {
157 type Result = ReductionEulerianPathToILP;
158
159 fn reduce_to(&self) -> Self::Result {
160 let arcs = self.graph().arcs();
161 let m = arcs.len();
162
163 // Empty-arc instance: vacuously feasible empty ILP.
164 if m == 0 {
165 let target = ILP::new(0, Vec::new(), Vec::new(), ObjectiveSense::Minimize);
166 return ReductionEulerianPathToILP {
167 target,
168 pairs: Vec::new(),
169 num_arcs: 0,
170 };
171 }
172
173 let pairs = compatible_pairs(&arcs);
174 let p = pairs.len();
175 let num_vars = p + 3 * m;
176
177 // Index helpers (mirroring the struct's accessors, but we need them
178 // before the struct exists).
179 let y_idx = |k: usize| -> usize { k };
180 let s_idx = |a: usize| -> usize { p + a };
181 let e_idx = |a: usize| -> usize { p + m + a };
182 let u_idx = |a: usize| -> usize { p + 2 * m + a };
183
184 // Inverse index: for each arc `a`, list `k` indices of pairs ending in
185 // `a` (`pairs[k].1 == a`) and pairs starting at `a` (`pairs[k].0 == a`).
186 let mut incoming: Vec<Vec<usize>> = vec![Vec::new(); m];
187 let mut outgoing: Vec<Vec<usize>> = vec![Vec::new(); m];
188 for (k, &(a, b)) in pairs.iter().enumerate() {
189 outgoing[a].push(k);
190 incoming[b].push(k);
191 }
192
193 let mut constraints: Vec<LinearConstraint> = Vec::new();
194
195 // (1) Predecessor equality: s_a + sum_{(b,a) in P} y_{b,a} = 1.
196 // (2) Successor equality: e_a + sum_{(a,b) in P} y_{a,b} = 1.
197 for a in 0..m {
198 let mut pred_terms: Vec<(usize, f64)> = vec![(s_idx(a), 1.0)];
199 for &k in &incoming[a] {
200 pred_terms.push((y_idx(k), 1.0));
201 }
202 constraints.push(LinearConstraint::eq(pred_terms, 1.0));
203
204 let mut succ_terms: Vec<(usize, f64)> = vec![(e_idx(a), 1.0)];
205 for &k in &outgoing[a] {
206 succ_terms.push((y_idx(k), 1.0));
207 }
208 constraints.push(LinearConstraint::eq(succ_terms, 1.0));
209 }
210
211 // (3) Binary upper bounds on start / end variables, and position
212 // upper bound on `u_a`.
213 for a in 0..m {
214 constraints.push(LinearConstraint::le(vec![(s_idx(a), 1.0)], 1.0));
215 constraints.push(LinearConstraint::le(vec![(e_idx(a), 1.0)], 1.0));
216 constraints.push(LinearConstraint::le(
217 vec![(u_idx(a), 1.0)],
218 (m as f64) - 1.0,
219 ));
220 }
221
222 // (4) Binary upper bounds on successor variables.
223 // (5) Order consistency (MTZ): u_b >= u_a + 1 - m * (1 - y_{a,b})
224 // i.e. u_a - u_b + m * y_{a,b} <= m - 1.
225 for (k, &(a, b)) in pairs.iter().enumerate() {
226 constraints.push(LinearConstraint::le(vec![(y_idx(k), 1.0)], 1.0));
227 constraints.push(LinearConstraint::le(
228 vec![(u_idx(a), 1.0), (u_idx(b), -1.0), (y_idx(k), m as f64)],
229 (m as f64) - 1.0,
230 ));
231 }
232
233 // (6) Unique start: sum_a s_a = 1.
234 // (7) Unique end: sum_a e_a = 1.
235 let start_sum: Vec<(usize, f64)> = (0..m).map(|a| (s_idx(a), 1.0)).collect();
236 let end_sum: Vec<(usize, f64)> = (0..m).map(|a| (e_idx(a), 1.0)).collect();
237 constraints.push(LinearConstraint::eq(start_sum, 1.0));
238 constraints.push(LinearConstraint::eq(end_sum, 1.0));
239
240 let target = ILP::new(num_vars, constraints, Vec::new(), ObjectiveSense::Minimize);
241
242 ReductionEulerianPathToILP {
243 target,
244 pairs,
245 num_arcs: m,
246 }
247 }
248}
249
250#[cfg(feature = "example-db")]
251pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
252 use crate::topology::DirectedGraph;
253 vec![crate::example_db::specs::RuleExampleSpec {
254 id: "eulerianpath_to_ilp",
255 build: || {
256 // Canonical issue #1025 instance: V = {0,1,2},
257 // A = [(0,1), (0,1), (1,2), (2,0)] (parallel arcs a_0, a_1).
258 // Witness ordering (a_0, a_2, a_3, a_1) traces 0->1->2->0->1.
259 let source =
260 EulerianPath::new(DirectedGraph::new(3, vec![(0, 1), (0, 1), (1, 2), (2, 0)]));
261 crate::example_db::specs::rule_example_via_ilp::<_, i32>(source)
262 },
263 }]
264}
265
266#[cfg(test)]
267#[path = "../unit_tests/rules/eulerianpath_ilp.rs"]
268mod tests;