problemreductions/rules/
stackercrane_ilp.rs

1//! Reduction from StackerCrane to ILP.
2//!
3//! One-hot position assignment for required arcs with McCormick products
4//! for consecutive-pair costs. Uses precomputed shortest-path connector
5//! distances.
6
7use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
8use crate::models::misc::StackerCrane;
9use crate::reduction;
10use crate::rules::ilp_helpers::one_hot_decode;
11use crate::rules::traits::{ReduceTo, ReductionResult};
12
13/// Result of reducing StackerCrane to ILP.
14///
15/// Variable layout (all binary):
16/// - `x_{i,p}` at index `i*m + p` for i,p in 0..m
17/// - `z_{i,j,p}` at index `m^2 + p*m^2 + i*m + j` for i,j,p in 0..m
18///
19/// Total: `m^2 + m^3` variables.
20#[derive(Debug, Clone)]
21pub struct ReductionSCToILP {
22    target: ILP<bool>,
23    num_arcs: usize,
24}
25
26impl ReductionResult for ReductionSCToILP {
27    type Source = StackerCrane;
28    type Target = ILP<bool>;
29
30    fn target_problem(&self) -> &ILP<bool> {
31        &self.target
32    }
33
34    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
35        // Decode the permutation: for each position p, find the arc a with x_{a,p} = 1
36        one_hot_decode(target_solution, self.num_arcs, self.num_arcs, 0)
37    }
38}
39
40#[reduction(
41    overhead = {
42        num_vars = "num_arcs * num_arcs + num_arcs * num_arcs * num_arcs",
43        num_constraints = "num_arcs + num_arcs + 3 * num_arcs * num_arcs * num_arcs",
44    }
45)]
46impl ReduceTo<ILP<bool>> for StackerCrane {
47    type Result = ReductionSCToILP;
48
49    fn reduce_to(&self) -> Self::Result {
50        let m = self.num_arcs();
51
52        if m == 0 {
53            return ReductionSCToILP {
54                target: ILP::new(0, vec![], vec![], ObjectiveSense::Minimize),
55                num_arcs: 0,
56            };
57        }
58
59        let num_vars = m * m + m * m * m;
60        let x_idx = |i: usize, p: usize| i * m + p;
61        let z_idx = |i: usize, j: usize, p: usize| m * m + p * m * m + i * m + j;
62
63        // Compute all-pairs shortest path distances in the mixed graph
64        let n = self.num_vertices();
65        let distances = all_pairs_shortest_paths(
66            n,
67            self.arcs(),
68            self.arc_lengths(),
69            self.edges(),
70            self.edge_lengths(),
71        );
72
73        let mut constraints = Vec::new();
74
75        // Each arc assigned to exactly one position: sum_p x_{i,p} = 1 for all i
76        for i in 0..m {
77            let terms: Vec<(usize, f64)> = (0..m).map(|p| (x_idx(i, p), 1.0)).collect();
78            constraints.push(LinearConstraint::eq(terms, 1.0));
79        }
80
81        // Each position assigned exactly one arc: sum_i x_{i,p} = 1 for all p
82        for p in 0..m {
83            let terms: Vec<(usize, f64)> = (0..m).map(|i| (x_idx(i, p), 1.0)).collect();
84            constraints.push(LinearConstraint::eq(terms, 1.0));
85        }
86
87        // McCormick linearization for z_{i,j,p} = x_{i,p} * x_{j,(p+1) mod m}
88        for p in 0..m {
89            let next_p = (p + 1) % m;
90            for i in 0..m {
91                for j in 0..m {
92                    let head_i = self.arcs()[i].1;
93                    let tail_j = self.arcs()[j].0;
94
95                    if distances[head_i][tail_j] == i64::MAX {
96                        // Infeasible pair: z_{i,j,p} = 0
97                        constraints.push(LinearConstraint::eq(vec![(z_idx(i, j, p), 1.0)], 0.0));
98                    } else {
99                        // z <= x_{i,p}
100                        constraints.push(LinearConstraint::le(
101                            vec![(z_idx(i, j, p), 1.0), (x_idx(i, p), -1.0)],
102                            0.0,
103                        ));
104                        // z <= x_{j, next_p}
105                        constraints.push(LinearConstraint::le(
106                            vec![(z_idx(i, j, p), 1.0), (x_idx(j, next_p), -1.0)],
107                            0.0,
108                        ));
109                        // z >= x_{i,p} + x_{j, next_p} - 1
110                        constraints.push(LinearConstraint::le(
111                            vec![
112                                (x_idx(i, p), 1.0),
113                                (x_idx(j, next_p), 1.0),
114                                (z_idx(i, j, p), -1.0),
115                            ],
116                            1.0,
117                        ));
118                    }
119                }
120            }
121        }
122
123        // Objective: minimize total walk length = sum_i l_i + sum_p sum_i sum_j D[head_i, tail_j] * z_{i,j,p}
124        // The constant sum_i l_i is ignored by the ILP solver (additive constant doesn't affect optimum).
125        let mut objective = Vec::new();
126        for p in 0..m {
127            for i in 0..m {
128                for j in 0..m {
129                    let head_i = self.arcs()[i].1;
130                    let tail_j = self.arcs()[j].0;
131                    let dist = distances[head_i][tail_j];
132                    if dist < i64::MAX {
133                        objective.push((z_idx(i, j, p), dist as f64));
134                    }
135                }
136            }
137        }
138
139        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize);
140
141        ReductionSCToILP {
142            target,
143            num_arcs: m,
144        }
145    }
146}
147
148/// All-pairs shortest paths via Floyd-Warshall on the mixed graph.
149fn all_pairs_shortest_paths(
150    n: usize,
151    arcs: &[(usize, usize)],
152    arc_lengths: &[i32],
153    edges: &[(usize, usize)],
154    edge_lengths: &[i32],
155) -> Vec<Vec<i64>> {
156    let mut dist = vec![vec![i64::MAX; n]; n];
157    for (i, row) in dist.iter_mut().enumerate() {
158        row[i] = 0;
159    }
160
161    // Directed arcs
162    for (&(u, v), &length) in arcs.iter().zip(arc_lengths) {
163        let cost = i64::from(length);
164        if cost < dist[u][v] {
165            dist[u][v] = cost;
166        }
167    }
168
169    // Undirected edges (both directions)
170    for (&(u, v), &length) in edges.iter().zip(edge_lengths) {
171        let cost = i64::from(length);
172        if cost < dist[u][v] {
173            dist[u][v] = cost;
174        }
175        if cost < dist[v][u] {
176            dist[v][u] = cost;
177        }
178    }
179
180    // Floyd-Warshall
181    for via in 0..n {
182        for src in 0..n {
183            if dist[src][via] == i64::MAX {
184                continue;
185            }
186            for dst in 0..n {
187                if dist[via][dst] == i64::MAX {
188                    continue;
189                }
190                let through = dist[src][via] + dist[via][dst];
191                if through < dist[src][dst] {
192                    dist[src][dst] = through;
193                }
194            }
195        }
196    }
197
198    dist
199}
200
201#[cfg(feature = "example-db")]
202pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
203    vec![crate::example_db::specs::RuleExampleSpec {
204        id: "stackercrane_to_ilp",
205        build: || {
206            // Simple: 3 vertices, 2 arcs, 1 edge
207            let source =
208                StackerCrane::new(3, vec![(0, 1), (2, 0)], vec![(1, 2)], vec![1, 1], vec![1]);
209            crate::example_db::specs::rule_example_via_ilp::<_, bool>(source)
210        },
211    }]
212}
213
214#[cfg(test)]
215#[path = "../unit_tests/rules/stackercrane_ilp.rs"]
216mod tests;
problemreductions/rules/stackercrane_ilp.rs

problemreductions/rules/
stackercrane_ilp.rs
