Skip to main content

problemreductions/rules/
acyclicpartition_ilp.rs

1//! Reduction from AcyclicPartition to ILP<i32>.
2//!
3//! One-hot assignment x_{v,c}, McCormick same-class indicators s_{t,c},
4//! crossing flags y_t, class ordering o_c, vertex-order copies p_v.
5//! See the paper entry for the full formulation.
6
7use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
8use crate::models::graph::AcyclicPartition;
9use crate::reduction;
10use crate::rules::ilp_helpers::mccormick_product;
11use crate::rules::traits::{ReduceTo, ReductionResult};
12
13#[derive(Debug, Clone)]
14pub struct ReductionAcyclicPartitionToILP {
15    target: ILP<i32>,
16    n: usize,
17}
18
19impl ReductionResult for ReductionAcyclicPartitionToILP {
20    type Source = AcyclicPartition<i32>;
21    type Target = ILP<i32>;
22
23    fn target_problem(&self) -> &ILP<i32> {
24        &self.target
25    }
26
27    /// One-hot decode: for each vertex v, output the unique c with x_{v,c} = 1.
28    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
29        let n = self.n;
30        (0..n)
31            .map(|v| {
32                (0..n)
33                    .find(|&c| target_solution[v * n + c] == 1)
34                    .unwrap_or(0)
35            })
36            .collect()
37    }
38}
39
40#[reduction(
41    overhead = {
42        num_vars = "num_vertices * num_vertices + num_arcs * num_vertices + num_arcs + 2 * num_vertices",
43        num_constraints = "num_vertices + num_vertices + num_arcs * num_vertices + num_arcs + 1 + 2 * num_vertices + 2 * num_vertices * num_vertices + num_arcs",
44    }
45)]
46impl ReduceTo<ILP<i32>> for AcyclicPartition<i32> {
47    type Result = ReductionAcyclicPartitionToILP;
48
49    fn reduce_to(&self) -> Self::Result {
50        let n = self.num_vertices();
51        let arcs = self.graph().arcs();
52        let m = arcs.len();
53
54        // Variable indices:
55        // x_{v,c} : v*n + c                          [0, n^2)
56        // s_{t,c} : n^2 + t*n + c                    [n^2, n^2 + m*n)
57        // y_t     : n^2 + m*n + t                     [n^2 + m*n, n^2 + m*n + m)
58        // o_c     : n^2 + m*n + m + c                 [n^2 + m*n + m, n^2 + m*n + m + n)
59        // p_v     : n^2 + m*n + m + n + v              [n^2 + m*n + m + n, n^2 + m*n + m + 2n)
60        let x_idx = |v: usize, c: usize| -> usize { v * n + c };
61        let s_idx = |t: usize, c: usize| -> usize { n * n + t * n + c };
62        let y_idx = |t: usize| -> usize { n * n + m * n + t };
63        let o_idx = |c: usize| -> usize { n * n + m * n + m + c };
64        let p_idx = |v: usize| -> usize { n * n + m * n + m + n + v };
65
66        let num_vars = n * n + m * n + m + 2 * n;
67        let mut constraints = Vec::new();
68        let big_m = n as f64;
69
70        // 1) Assignment: Σ_c x_{v,c} = 1  for each vertex v
71        for v in 0..n {
72            let terms: Vec<(usize, f64)> = (0..n).map(|c| (x_idx(v, c), 1.0)).collect();
73            constraints.push(LinearConstraint::eq(terms, 1.0));
74        }
75
76        // 2) Weight bound: Σ_v w_v * x_{v,c} ≤ B  for each class c
77        for c in 0..n {
78            let terms: Vec<(usize, f64)> = self
79                .vertex_weights()
80                .iter()
81                .enumerate()
82                .map(|(v, &w)| (x_idx(v, c), w as f64))
83                .collect();
84            constraints.push(LinearConstraint::le(terms, *self.weight_bound() as f64));
85        }
86
87        // 3) McCormick: s_{t,c} = x_{u_t,c} * x_{v_t,c}
88        for (t, &(u, v)) in arcs.iter().enumerate() {
89            for c in 0..n {
90                constraints.extend(mccormick_product(s_idx(t, c), x_idx(u, c), x_idx(v, c)));
91            }
92        }
93
94        // 4) Crossing: y_t + Σ_c s_{t,c} = 1
95        for t in 0..m {
96            let mut terms: Vec<(usize, f64)> = vec![(y_idx(t), 1.0)];
97            for c in 0..n {
98                terms.push((s_idx(t, c), 1.0));
99            }
100            constraints.push(LinearConstraint::eq(terms, 1.0));
101        }
102
103        // 5) Cost bound: Σ_t cost(a_t) * y_t ≤ K
104        let cost_terms: Vec<(usize, f64)> = self
105            .arc_costs()
106            .iter()
107            .enumerate()
108            .map(|(t, &c)| (y_idx(t), c as f64))
109            .collect();
110        constraints.push(LinearConstraint::le(cost_terms, *self.cost_bound() as f64));
111
112        // 6) Order bounds: 0 ≤ o_c ≤ n-1, 0 ≤ p_v ≤ n-1
113        for c in 0..n {
114            constraints.push(LinearConstraint::ge(vec![(o_idx(c), 1.0)], 0.0));
115            constraints.push(LinearConstraint::le(vec![(o_idx(c), 1.0)], (n - 1) as f64));
116        }
117        for v in 0..n {
118            constraints.push(LinearConstraint::ge(vec![(p_idx(v), 1.0)], 0.0));
119            constraints.push(LinearConstraint::le(vec![(p_idx(v), 1.0)], (n - 1) as f64));
120        }
121
122        // 7) Link p_v to o_c: p_v - o_c ≤ (n-1)(1 - x_{v,c}) and o_c - p_v ≤ (n-1)(1 - x_{v,c})
123        for v in 0..n {
124            for c in 0..n {
125                // p_v - o_c + (n-1)*x_{v,c} ≤ n-1
126                constraints.push(LinearConstraint::le(
127                    vec![
128                        (p_idx(v), 1.0),
129                        (o_idx(c), -1.0),
130                        (x_idx(v, c), (n - 1) as f64),
131                    ],
132                    (n - 1) as f64,
133                ));
134                // o_c - p_v + (n-1)*x_{v,c} ≤ n-1
135                constraints.push(LinearConstraint::le(
136                    vec![
137                        (o_idx(c), 1.0),
138                        (p_idx(v), -1.0),
139                        (x_idx(v, c), (n - 1) as f64),
140                    ],
141                    (n - 1) as f64,
142                ));
143            }
144        }
145
146        // 8) DAG ordering: p_{v_t} - p_{u_t} ≥ 1 - n * Σ_c s_{t,c}
147        //    i.e., p_{v_t} - p_{u_t} + n * Σ_c s_{t,c} ≥ 1
148        for (t, &(u, v)) in arcs.iter().enumerate() {
149            let mut terms = vec![(p_idx(v), 1.0), (p_idx(u), -1.0)];
150            for c in 0..n {
151                terms.push((s_idx(t, c), big_m));
152            }
153            constraints.push(LinearConstraint::ge(terms, 1.0));
154        }
155
156        let target = ILP::new(num_vars, constraints, vec![], ObjectiveSense::Minimize);
157
158        ReductionAcyclicPartitionToILP { target, n }
159    }
160}
161
162#[cfg(feature = "example-db")]
163pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
164    use crate::export::SolutionPair;
165    use crate::topology::DirectedGraph;
166    vec![crate::example_db::specs::RuleExampleSpec {
167        id: "acyclicpartition_to_ilp",
168        build: || {
169            let source = AcyclicPartition::new(
170                DirectedGraph::new(4, vec![(0, 1), (1, 2), (2, 3)]),
171                vec![1, 1, 1, 1],
172                vec![1, 1, 1],
173                3,
174                2,
175            );
176            let reduction: ReductionAcyclicPartitionToILP =
177                crate::rules::ReduceTo::<ILP<i32>>::reduce_to(&source);
178            let ilp_sol = crate::solvers::ILPSolver::new()
179                .solve(reduction.target_problem())
180                .expect("ILP should be solvable");
181            let extracted = reduction.extract_solution(&ilp_sol);
182            crate::example_db::specs::rule_example_with_witness::<_, ILP<i32>>(
183                source,
184                SolutionPair {
185                    source_config: extracted,
186                    target_config: ilp_sol,
187                },
188            )
189        },
190    }]
191}
192
193#[cfg(test)]
194#[path = "../unit_tests/rules/acyclicpartition_ilp.rs"]
195mod tests;