Skip to main content

problemreductions/rules/
consecutiveonesmatrixaugmentation_ilp.rs

1//! Reduction from ConsecutiveOnesMatrixAugmentation to ILP.
2//!
3//! Choose a column permutation and, for each row, choose an interval that will
4//! become its consecutive block of 1s. Flips are needed only for zeros inside
5//! that interval.
6
7use crate::models::algebraic::{
8    ConsecutiveOnesMatrixAugmentation, LinearConstraint, ObjectiveSense, ILP,
9};
10use crate::reduction;
11use crate::rules::ilp_helpers::{one_hot_assignment_constraints, one_hot_decode};
12use crate::rules::traits::{ReduceTo, ReductionResult};
13
14#[derive(Debug, Clone)]
15pub struct ReductionCOMAToILP {
16    target: ILP<bool>,
17    num_cols: usize,
18}
19
20impl ReductionResult for ReductionCOMAToILP {
21    type Source = ConsecutiveOnesMatrixAugmentation;
22    type Target = ILP<bool>;
23
24    fn target_problem(&self) -> &ILP<bool> {
25        &self.target
26    }
27
28    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
29        one_hot_decode(target_solution, self.num_cols, self.num_cols, 0)
30    }
31}
32
33#[reduction(
34    overhead = {
35        num_vars = "num_cols * num_cols + 5 * num_rows * num_cols",
36        num_constraints = "num_cols + num_cols + num_rows * num_cols + 2 * num_rows + num_rows + 3 * num_rows * num_cols + 4 * num_rows * num_cols + 1",
37    }
38)]
39impl ReduceTo<ILP<bool>> for ConsecutiveOnesMatrixAugmentation {
40    type Result = ReductionCOMAToILP;
41
42    fn reduce_to(&self) -> Self::Result {
43        let m = self.num_rows();
44        let n = self.num_cols();
45
46        // Variable layout (all binary):
47        // x_{c,p}: n^2 at [0, n^2)
48        // a_{r,p}: m*n at [n^2, n^2 + m*n)
49        // l_{r,p}: m*n at [n^2 + m*n, n^2 + 2*m*n)
50        // u_{r,p}: m*n at [n^2 + 2*m*n, n^2 + 3*m*n)
51        // h_{r,p}: m*n at [n^2 + 3*m*n, n^2 + 4*m*n)
52        // f_{r,p}: m*n at [n^2 + 4*m*n, n^2 + 5*m*n)
53        let x_off = 0;
54        let a_off = n * n;
55        let l_off = n * n + m * n;
56        let u_off = n * n + 2 * m * n;
57        let h_off = n * n + 3 * m * n;
58        let f_off = n * n + 4 * m * n;
59        let num_vars = n * n + 5 * m * n;
60
61        let mut constraints = Vec::new();
62
63        // One-hot permutation assignment
64        constraints.extend(one_hot_assignment_constraints(n, n, x_off));
65
66        // a_{r,p} = sum_c A_{r,c} * x_{c,p}
67        for r in 0..m {
68            for p in 0..n {
69                let a_idx = a_off + r * n + p;
70                let mut terms = vec![(a_idx, 1.0)];
71                for c in 0..n {
72                    if self.matrix()[r][c] {
73                        terms.push((x_off + c * n + p, -1.0));
74                    }
75                }
76                constraints.push(LinearConstraint::eq(terms, 0.0));
77            }
78        }
79
80        // Per-row interval constraints
81        for r in 0..m {
82            let beta_r: f64 = if self.matrix()[r].iter().any(|&v| v) {
83                1.0
84            } else {
85                0.0
86            };
87
88            // sum_p l_{r,p} = beta_r
89            let l_terms: Vec<(usize, f64)> = (0..n).map(|p| (l_off + r * n + p, 1.0)).collect();
90            constraints.push(LinearConstraint::eq(l_terms, beta_r));
91
92            // sum_p u_{r,p} = beta_r
93            let u_terms: Vec<(usize, f64)> = (0..n).map(|p| (u_off + r * n + p, 1.0)).collect();
94            constraints.push(LinearConstraint::eq(u_terms, beta_r));
95
96            // sum_p p*l_{r,p} <= sum_p p*u_{r,p} + (n-1)*(1 - beta_r)
97            // => sum_p p*l_{r,p} - sum_p p*u_{r,p} <= (n-1)*(1 - beta_r)
98            let mut order_terms = Vec::new();
99            for p in 0..n {
100                order_terms.push((l_off + r * n + p, p as f64));
101                order_terms.push((u_off + r * n + p, -(p as f64)));
102            }
103            constraints.push(LinearConstraint::le(
104                order_terms,
105                (n as f64 - 1.0) * (1.0 - beta_r),
106            ));
107
108            for p in 0..n {
109                let h_idx = h_off + r * n + p;
110                let a_idx = a_off + r * n + p;
111                let f_idx = f_off + r * n + p;
112
113                // h_{r,p} <= sum_{q=0}^{p} l_{r,q}
114                let l_prefix: Vec<(usize, f64)> =
115                    (0..=p).map(|q| (l_off + r * n + q, -1.0)).collect();
116                let mut h_le_l = vec![(h_idx, 1.0)];
117                h_le_l.extend(l_prefix);
118                constraints.push(LinearConstraint::le(h_le_l, 0.0));
119
120                // h_{r,p} <= sum_{q=p}^{n-1} u_{r,q}
121                let u_suffix: Vec<(usize, f64)> =
122                    (p..n).map(|q| (u_off + r * n + q, -1.0)).collect();
123                let mut h_le_u = vec![(h_idx, 1.0)];
124                h_le_u.extend(u_suffix);
125                constraints.push(LinearConstraint::le(h_le_u, 0.0));
126
127                // h_{r,p} >= sum_{q=0}^{p} l_{r,q} + sum_{q=p}^{n-1} u_{r,q} - 1
128                let mut h_ge_terms = vec![(h_idx, 1.0)];
129                for q in 0..=p {
130                    h_ge_terms.push((l_off + r * n + q, -1.0));
131                }
132                for q in p..n {
133                    h_ge_terms.push((u_off + r * n + q, -1.0));
134                }
135                constraints.push(LinearConstraint::ge(h_ge_terms, -1.0));
136
137                // a_{r,p} <= h_{r,p}
138                constraints.push(LinearConstraint::le(vec![(a_idx, 1.0), (h_idx, -1.0)], 0.0));
139
140                // h_{r,p} <= a_{r,p} + f_{r,p}
141                constraints.push(LinearConstraint::le(
142                    vec![(h_idx, 1.0), (a_idx, -1.0), (f_idx, -1.0)],
143                    0.0,
144                ));
145
146                // f_{r,p} <= h_{r,p}
147                constraints.push(LinearConstraint::le(vec![(f_idx, 1.0), (h_idx, -1.0)], 0.0));
148
149                // f_{r,p} + a_{r,p} <= 1
150                constraints.push(LinearConstraint::le(vec![(f_idx, 1.0), (a_idx, 1.0)], 1.0));
151            }
152        }
153
154        // Augmentation budget: sum f_{r,p} <= K
155        let mut budget_terms = Vec::new();
156        for r in 0..m {
157            for p in 0..n {
158                budget_terms.push((f_off + r * n + p, 1.0));
159            }
160        }
161        constraints.push(LinearConstraint::le(budget_terms, self.bound() as f64));
162
163        let target = ILP::new(num_vars, constraints, vec![], ObjectiveSense::Minimize);
164        ReductionCOMAToILP {
165            target,
166            num_cols: n,
167        }
168    }
169}
170
171#[cfg(feature = "example-db")]
172pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
173    use crate::export::SolutionPair;
174    vec![crate::example_db::specs::RuleExampleSpec {
175        id: "consecutiveonesmatrixaugmentation_to_ilp",
176        build: || {
177            let source = ConsecutiveOnesMatrixAugmentation::new(
178                vec![vec![true, false, true], vec![false, true, true]],
179                1,
180            );
181            // Identity permutation [0,1,2]:
182            // Row 0: [1,0,1] needs 1 flip (the middle 0), cost=1
183            // Row 1: [0,1,1] needs 0 flips, cost=0
184            // Total = 1 <= 1
185            let reduction: ReductionCOMAToILP = ReduceTo::<ILP<bool>>::reduce_to(&source);
186            let ilp_solver = crate::solvers::ILPSolver::new();
187            let target_config = ilp_solver
188                .solve(reduction.target_problem())
189                .expect("ILP should be solvable");
190            let extracted = reduction.extract_solution(&target_config);
191            crate::example_db::specs::rule_example_with_witness::<_, ILP<bool>>(
192                source,
193                SolutionPair {
194                    source_config: extracted,
195                    target_config,
196                },
197            )
198        },
199    }]
200}
201
202#[cfg(test)]
203#[path = "../unit_tests/rules/consecutiveonesmatrixaugmentation_ilp.rs"]
204mod tests;