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::{Bounded, IntegerVariable, 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<i64, i64, Bounded>`.
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<i64, i64, Bounded>,
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<i64, i64, Bounded>;
62
63    fn target_problem(&self) -> &ILP<i64, i64, Bounded> {
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`. Malformed assignments return an extraction
72    /// error instead of fabricating an ordering.
73    fn extract_solution(
74        &self,
75        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
76    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
77        crate::rules::traits::validate_target_witness(
78            self.target_problem(),
79            target_solution,
80            |value| value.value.is_some(),
81            "target ILP assignment is infeasible",
82        )?;
83
84        Ok({
85            let m = self.num_arcs;
86            if m == 0 {
87                return Ok(Vec::new());
88            }
89
90            // Find the unique active start arc.
91            let mut current = match (0..m).find(|&a| target_solution[self.s_idx(a)] == 1) {
92                Some(a) => a,
93                None => {
94                    return Err(crate::rules::ExtractionError::invalid(
95                        "ILP witness has no active Eulerian-path start arc",
96                    ));
97                }
98            };
99
100            // Walk the active successor relation, recording each visited arc.
101            let mut order = Vec::with_capacity(m);
102            let mut visited = vec![false; m];
103            order.push(current);
104            visited[current] = true;
105
106            for _ in 1..m {
107                let next = self
108                    .pairs
109                    .iter()
110                    .enumerate()
111                    .find(|&(k, &(a, _))| a == current && target_solution[k] == 1)
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                        return Err(crate::rules::ExtractionError::invalid(format!(
122                            "ILP witness has no unvisited successor for arc {current}",
123                        )));
124                    }
125                }
126            }
127            order
128        })
129    }
130}
131
132/// Enumerate compatible ordered pairs `(a, b)` with `a != b` and
133/// `head(a) = tail(b)`. The order is `a`-major then `b`-major, matching the
134/// nested-loop construction below.
135fn compatible_pairs(arcs: &[(usize, usize)]) -> Vec<(usize, usize)> {
136    let mut pairs = Vec::new();
137    for (a, &(_, head_a)) in arcs.iter().enumerate() {
138        for (b, &(tail_b, _)) in arcs.iter().enumerate() {
139            if a != b && tail_b == head_a {
140                pairs.push((a, b));
141            }
142        }
143    }
144    pairs
145}
146
147#[crate::aggregate_reduction(ilp_feasibility)]
148impl crate::rules::AggregateReductionResult for ReductionEulerianPathToILP {}
149
150// For m arcs and P compatible distinct ordered pairs, row counts give exactly
151// 7m+6P nonzeros when m>0. P<=m(m-1), and the empty construction has zero.
152#[reduction(transform = upper_bound {
153    max_constraint_magnitude_bits = "num_arcs + 1",
154    num_vars = "3 * num_arcs + num_arcs * num_arcs",
155    num_constraints = "5 * num_arcs + 2 * num_arcs * num_arcs + 2",
156    num_nonzeros = "num_arcs + 6 * num_arcs * num_arcs",
157})]
158impl ReduceTo<ILP<i64, i64, Bounded>> for EulerianPath {
159    type Result = ReductionEulerianPathToILP;
160
161    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
162        let arcs = self.graph().arcs();
163        let m = arcs.len();
164
165        // Empty-arc instance: vacuously feasible empty ILP.
166        if m == 0 {
167            let target = ILP::new(0, Vec::new(), Vec::new(), ObjectiveSense::Minimize)
168                .map_err(Self::target_construction)?;
169            return Ok(ReductionEulerianPathToILP {
170                target,
171                pairs: Vec::new(),
172                num_arcs: 0,
173            });
174        }
175
176        let pairs = compatible_pairs(&arcs);
177        let p = pairs.len();
178        let num_vars = p + 3 * m;
179
180        // Index helpers (mirroring the struct's accessors, but we need them
181        // before the struct exists).
182        let y_idx = |k: usize| -> usize { k };
183        let s_idx = |a: usize| -> usize { p + a };
184        let e_idx = |a: usize| -> usize { p + m + a };
185        let u_idx = |a: usize| -> usize { p + 2 * m + a };
186
187        // Inverse index: for each arc `a`, list `k` indices of pairs ending in
188        // `a` (`pairs[k].1 == a`) and pairs starting at `a` (`pairs[k].0 == a`).
189        let mut incoming: Vec<Vec<usize>> = vec![Vec::new(); m];
190        let mut outgoing: Vec<Vec<usize>> = vec![Vec::new(); m];
191        for (k, &(a, b)) in pairs.iter().enumerate() {
192            outgoing[a].push(k);
193            incoming[b].push(k);
194        }
195
196        let mut constraints: Vec<LinearConstraint> = Vec::new();
197
198        // (1) Predecessor equality: s_a + sum_{(b,a) in P} y_{b,a} = 1.
199        // (2) Successor equality:   e_a + sum_{(a,b) in P} y_{a,b} = 1.
200        let m_i64 = Self::exact_i64(m, "encoding the arc order")?;
201        for a in 0..m {
202            let mut pred_terms: Vec<(usize, i64)> = vec![(s_idx(a), 1)];
203            for &k in &incoming[a] {
204                pred_terms.push((y_idx(k), 1));
205            }
206            constraints.push(LinearConstraint::eq(pred_terms, 1));
207
208            let mut succ_terms: Vec<(usize, i64)> = vec![(e_idx(a), 1)];
209            for &k in &outgoing[a] {
210                succ_terms.push((y_idx(k), 1));
211            }
212            constraints.push(LinearConstraint::eq(succ_terms, 1));
213        }
214
215        // (3) Binary upper bounds on start / end variables, and position
216        //     upper bound on `u_a`.
217        for a in 0..m {
218            constraints.push(LinearConstraint::le(vec![(s_idx(a), 1)], 1));
219            constraints.push(LinearConstraint::le(vec![(e_idx(a), 1)], 1));
220            constraints.push(LinearConstraint::le(vec![(u_idx(a), 1)], m_i64 - 1));
221        }
222
223        // (4) Binary upper bounds on successor variables.
224        // (5) Order consistency (MTZ): u_b >= u_a + 1 - m * (1 - y_{a,b})
225        //     i.e.  u_a - u_b + m * y_{a,b} <= m - 1.
226        for (k, &(a, b)) in pairs.iter().enumerate() {
227            constraints.push(LinearConstraint::le(vec![(y_idx(k), 1)], 1));
228            constraints.push(LinearConstraint::le(
229                vec![(u_idx(a), 1), (u_idx(b), -1), (y_idx(k), m_i64)],
230                m_i64 - 1,
231            ));
232        }
233
234        // (6) Unique start: sum_a s_a = 1.
235        // (7) Unique end:   sum_a e_a = 1.
236        let start_sum: Vec<(usize, i64)> = (0..m).map(|a| (s_idx(a), 1)).collect();
237        let end_sum: Vec<(usize, i64)> = (0..m).map(|a| (e_idx(a), 1)).collect();
238        constraints.push(LinearConstraint::eq(start_sum, 1));
239        constraints.push(LinearConstraint::eq(end_sum, 1));
240
241        let mut variables = vec![IntegerVariable::binary(); num_vars];
242        variables[p + 2 * m..].fill(
243            IntegerVariable::new(Some(0), Some(m_i64 - 1)).map_err(Self::target_construction)?,
244        );
245
246        let target =
247            ILP::with_variables(variables, constraints, Vec::new(), ObjectiveSense::Minimize)
248                .map_err(Self::target_construction)?;
249
250        Ok(ReductionEulerianPathToILP {
251            target,
252            pairs,
253            num_arcs: m,
254        })
255    }
256}
257
258#[cfg(feature = "example-db")]
259pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
260    use crate::topology::DirectedGraph;
261    vec![crate::example_db::specs::RuleExampleSpec {
262        id: "eulerianpath_to_ilp",
263        build: || {
264            // Canonical issue #1025 instance: V = {0,1,2},
265            // A = [(0,1), (0,1), (1,2), (2,0)] (parallel arcs a_0, a_1).
266            // Witness ordering (a_0, a_2, a_3, a_1) traces 0->1->2->0->1.
267            let source =
268                EulerianPath::new(DirectedGraph::new(3, vec![(0, 1), (0, 1), (1, 2), (2, 0)]));
269            crate::example_db::specs::rule_example_via_bounded_ilp::<_>(source)
270        },
271    }]
272}
273
274#[cfg(test)]
275#[path = "../unit_tests/rules/eulerianpath_ilp.rs"]
276mod tests;
problemreductions/rules/eulerianpath_ilp.rs

problemreductions/rules/
eulerianpath_ilp.rs
