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::{mccormick_product, 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(
35        &self,
36        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
37    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
38        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
39
40        Ok({
41            // Decode the permutation: for each position p, find the arc a with x_{a,p} = 1
42            one_hot_decode(target_solution, self.num_arcs, self.num_arcs, 0)?
43        })
44    }
45}
46
47#[reduction(
48    transform = upper_bound {
49        num_vars = "num_arcs * num_arcs + num_arcs * num_arcs * num_arcs",
50        num_constraints = "num_arcs + num_arcs + 4 * num_arcs * num_arcs * num_arcs",
51    },
52    unavailable = {
53        num_nonzeros = "the exact target parameter is not represented by this reduction's symbolic transform",
54    }
55)]
56impl ReduceTo<ILP<bool>> for StackerCrane {
57    type Result = ReductionSCToILP;
58
59    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
60        let m = self.num_arcs();
61
62        if m == 0 {
63            return Ok(ReductionSCToILP {
64                target: ILP::new(0, vec![], vec![], ObjectiveSense::Minimize)
65                    .map_err(Self::target_construction)?,
66                num_arcs: 0,
67            });
68        }
69
70        let num_vars = m * m + m * m * m;
71        let x_idx = |i: usize, p: usize| i * m + p;
72        let z_idx = |i: usize, j: usize, p: usize| m * m + p * m * m + i * m + j;
73
74        // Compute all-pairs shortest path distances in the mixed graph
75        let n = self.num_vertices();
76        let distances = all_pairs_shortest_paths(
77            n,
78            self.arcs(),
79            self.arc_lengths(),
80            self.edges(),
81            self.edge_lengths(),
82        );
83
84        let mut constraints = Vec::new();
85
86        // Each arc assigned to exactly one position: sum_p x_{i,p} = 1 for all i
87        for i in 0..m {
88            let terms: Vec<(usize, i64)> = (0..m).map(|p| (x_idx(i, p), 1)).collect();
89            constraints.push(LinearConstraint::eq(terms, 1));
90        }
91
92        // Each position assigned exactly one arc: sum_i x_{i,p} = 1 for all p
93        for p in 0..m {
94            let terms: Vec<(usize, i64)> = (0..m).map(|i| (x_idx(i, p), 1)).collect();
95            constraints.push(LinearConstraint::eq(terms, 1));
96        }
97
98        // McCormick linearization for z_{i,j,p} = x_{i,p} * x_{j,(p+1) mod m}
99        for p in 0..m {
100            let next_p = (p + 1) % m;
101            for i in 0..m {
102                for j in 0..m {
103                    let head_i = self.arcs()[i].1;
104                    let tail_j = self.arcs()[j].0;
105
106                    // The product relation is required for every pair, even
107                    // unreachable ones: z = 0 alone does not forbid adjacency.
108                    constraints.extend(mccormick_product(
109                        z_idx(i, j, p),
110                        x_idx(i, p),
111                        x_idx(j, next_p),
112                    ));
113                    if distances[head_i][tail_j] == i64::MAX {
114                        constraints.push(LinearConstraint::eq(vec![(z_idx(i, j, p), 1)], 0));
115                    }
116                }
117            }
118        }
119
120        // Objective: minimize total walk length = sum_i l_i + sum_p sum_i sum_j D[head_i, tail_j] * z_{i,j,p}
121        // The constant sum_i l_i is ignored by the ILP solver (additive constant doesn't affect optimum).
122        let mut objective = Vec::new();
123        for p in 0..m {
124            for i in 0..m {
125                for j in 0..m {
126                    let head_i = self.arcs()[i].1;
127                    let tail_j = self.arcs()[j].0;
128                    let dist = distances[head_i][tail_j];
129                    if dist < i64::MAX {
130                        objective.push((z_idx(i, j, p), dist));
131                    }
132                }
133            }
134        }
135
136        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize)
137            .map_err(Self::target_construction)?;
138
139        Ok(ReductionSCToILP {
140            target,
141            num_arcs: m,
142        })
143    }
144}
145
146/// All-pairs shortest paths via Floyd-Warshall on the mixed graph.
147fn all_pairs_shortest_paths(
148    n: usize,
149    arcs: &[(usize, usize)],
150    arc_lengths: &[i64],
151    edges: &[(usize, usize)],
152    edge_lengths: &[i64],
153) -> Vec<Vec<i64>> {
154    let mut dist = vec![vec![i64::MAX; n]; n];
155    for (i, row) in dist.iter_mut().enumerate() {
156        row[i] = 0;
157    }
158
159    // Directed arcs
160    for (&(u, v), &length) in arcs.iter().zip(arc_lengths) {
161        let cost = length;
162        if cost < dist[u][v] {
163            dist[u][v] = cost;
164        }
165    }
166
167    // Undirected edges (both directions)
168    for (&(u, v), &length) in edges.iter().zip(edge_lengths) {
169        let cost = length;
170        if cost < dist[u][v] {
171            dist[u][v] = cost;
172        }
173        if cost < dist[v][u] {
174            dist[v][u] = cost;
175        }
176    }
177
178    // Floyd-Warshall
179    for via in 0..n {
180        for src in 0..n {
181            if dist[src][via] == i64::MAX {
182                continue;
183            }
184            for dst in 0..n {
185                if dist[via][dst] == i64::MAX {
186                    continue;
187                }
188                let through = dist[src][via] + dist[via][dst];
189                if through < dist[src][dst] {
190                    dist[src][dst] = through;
191                }
192            }
193        }
194    }
195
196    dist
197}
198
199#[cfg(feature = "example-db")]
200pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
201    vec![crate::example_db::specs::RuleExampleSpec {
202        id: "stackercrane_to_ilp",
203        build: || {
204            // Simple: 3 vertices, 2 arcs, 1 edge
205            let source =
206                StackerCrane::new(3, vec![(0, 1), (2, 0)], vec![(1, 2)], vec![1, 1], vec![1]);
207            crate::example_db::specs::rule_example_via_ilp::<_, bool>(source)
208        },
209    }]
210}
211
212#[cfg(test)]
213#[path = "../unit_tests/rules/stackercrane_ilp.rs"]
214mod tests;
problemreductions/rules/stackercrane_ilp.rs

problemreductions/rules/
stackercrane_ilp.rs
