problemreductions/rules/
consecutiveonessubmatrix_ilp.rs

1//! Reduction from ConsecutiveOnesSubmatrix to ILP.
2//!
3//! Select exactly K columns, permute only those selected columns, and require
4//! every row to have a single consecutive block within the chosen submatrix.
5//! The output is the column-selection bits s_c.
6
7use crate::models::algebraic::{ConsecutiveOnesSubmatrix, LinearConstraint, ObjectiveSense, ILP};
8use crate::reduction;
9use crate::rules::traits::{ReduceTo, ReductionResult};
10
11#[derive(Debug, Clone)]
12pub struct ReductionCOSToILP {
13    target: ILP<bool>,
14    num_cols: usize,
15}
16
17impl ReductionResult for ReductionCOSToILP {
18    type Source = ConsecutiveOnesSubmatrix;
19    type Target = ILP<bool>;
20
21    fn target_problem(&self) -> &ILP<bool> {
22        &self.target
23    }
24
25    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
26        // Output the selection bits s_c (first num_cols variables)
27        target_solution[..self.num_cols].to_vec()
28    }
29}
30
31#[reduction(
32    overhead = {
33        num_vars = "num_cols + num_cols * bound + 5 * num_rows * bound",
34        num_constraints = "1 + num_cols + bound + num_rows * bound + 2 * num_rows + num_rows + 3 * num_rows * bound + 4 * num_rows * bound",
35    }
36)]
37impl ReduceTo<ILP<bool>> for ConsecutiveOnesSubmatrix {
38    type Result = ReductionCOSToILP;
39
40    fn reduce_to(&self) -> Self::Result {
41        let m = self.num_rows();
42        let n = self.num_cols();
43        let k = self.bound() as usize;
44
45        // Variable layout (all binary):
46        // s_c: n vars at [0, n)  — column selection
47        // x_{c,p}: n*K vars at [n, n + n*K)  — column c placed at position p in [0..K)
48        // a_{r,p}: m*K vars at [n + n*K, n + n*K + m*K)  — value at row r, position p
49        // l_{r,p}: m*K vars  — left boundary
50        // u_{r,p}: m*K vars  — right boundary
51        // h_{r,p}: m*K vars  — inside interval
52        // f_{r,p}: m*K vars  — flip indicator (not used for budget, but needed for C1P)
53        let s_off = 0;
54        let x_off = n;
55        let a_off = n + n * k;
56        let l_off = a_off + m * k;
57        let u_off = l_off + m * k;
58        let h_off = u_off + m * k;
59        let f_off = h_off + m * k;
60        let num_vars = f_off + m * k;
61
62        let mut constraints = Vec::new();
63
64        // sum_c s_c = K
65        let s_terms: Vec<(usize, f64)> = (0..n).map(|c| (s_off + c, 1.0)).collect();
66        constraints.push(LinearConstraint::eq(s_terms, k as f64));
67
68        // sum_p x_{c,p} = s_c for all c
69        for c in 0..n {
70            let mut terms: Vec<(usize, f64)> = (0..k).map(|p| (x_off + c * k + p, 1.0)).collect();
71            terms.push((s_off + c, -1.0));
72            constraints.push(LinearConstraint::eq(terms, 0.0));
73        }
74
75        // sum_c x_{c,p} = 1 for all p in {0, ..., K-1}
76        for p in 0..k {
77            let terms: Vec<(usize, f64)> = (0..n).map(|c| (x_off + c * k + p, 1.0)).collect();
78            constraints.push(LinearConstraint::eq(terms, 1.0));
79        }
80
81        // a_{r,p} = sum_c A_{r,c} * x_{c,p}
82        for r in 0..m {
83            for p in 0..k {
84                let a_idx = a_off + r * k + p;
85                let mut terms = vec![(a_idx, 1.0)];
86                for c in 0..n {
87                    if self.matrix()[r][c] {
88                        terms.push((x_off + c * k + p, -1.0));
89                    }
90                }
91                constraints.push(LinearConstraint::eq(terms, 0.0));
92            }
93        }
94
95        // C1P interval constraints on the K-position permuted submatrix
96        for r in 0..m {
97            // beta_r = 1 if row r has at least one 1 in the original matrix
98            // (among any column, not just selected ones — the ILP will determine)
99            // We use beta_r = 1 for rows that have any 1, to allow intervals
100            let beta_r: f64 = if self.matrix()[r].iter().any(|&v| v) {
101                1.0
102            } else {
103                0.0
104            };
105
106            // sum_p l_{r,p} = beta_r
107            let l_terms: Vec<(usize, f64)> = (0..k).map(|p| (l_off + r * k + p, 1.0)).collect();
108            constraints.push(LinearConstraint::eq(l_terms, beta_r));
109
110            // sum_p u_{r,p} = beta_r
111            let u_terms: Vec<(usize, f64)> = (0..k).map(|p| (u_off + r * k + p, 1.0)).collect();
112            constraints.push(LinearConstraint::eq(u_terms, beta_r));
113
114            // sum_p p*l_{r,p} <= sum_p p*u_{r,p} + (K-1)*(1 - beta_r)
115            if k > 0 {
116                let mut order_terms = Vec::new();
117                for p in 0..k {
118                    order_terms.push((l_off + r * k + p, p as f64));
119                    order_terms.push((u_off + r * k + p, -(p as f64)));
120                }
121                constraints.push(LinearConstraint::le(
122                    order_terms,
123                    (k as f64 - 1.0).max(0.0) * (1.0 - beta_r),
124                ));
125            }
126
127            for p in 0..k {
128                let h_idx = h_off + r * k + p;
129                let a_idx = a_off + r * k + p;
130                let f_idx = f_off + r * k + p;
131
132                // h_{r,p} <= sum_{q=0}^{p} l_{r,q}
133                let mut h_le_l = vec![(h_idx, 1.0)];
134                for q in 0..=p {
135                    h_le_l.push((l_off + r * k + q, -1.0));
136                }
137                constraints.push(LinearConstraint::le(h_le_l, 0.0));
138
139                // h_{r,p} <= sum_{q=p}^{K-1} u_{r,q}
140                let mut h_le_u = vec![(h_idx, 1.0)];
141                for q in p..k {
142                    h_le_u.push((u_off + r * k + q, -1.0));
143                }
144                constraints.push(LinearConstraint::le(h_le_u, 0.0));
145
146                // h_{r,p} >= sum_{q=0}^{p} l_{r,q} + sum_{q=p}^{K-1} u_{r,q} - 1
147                let mut h_ge_terms = vec![(h_idx, 1.0)];
148                for q in 0..=p {
149                    h_ge_terms.push((l_off + r * k + q, -1.0));
150                }
151                for q in p..k {
152                    h_ge_terms.push((u_off + r * k + q, -1.0));
153                }
154                constraints.push(LinearConstraint::ge(h_ge_terms, -1.0));
155
156                // a_{r,p} <= h_{r,p}  — every 1 must be inside the interval
157                constraints.push(LinearConstraint::le(vec![(a_idx, 1.0), (h_idx, -1.0)], 0.0));
158
159                // For C1P (no augmentation): the interval must exactly cover the 1s
160                // h_{r,p} <= a_{r,p} + f_{r,p} — position inside interval but 0 costs a flip
161                constraints.push(LinearConstraint::le(
162                    vec![(h_idx, 1.0), (a_idx, -1.0), (f_idx, -1.0)],
163                    0.0,
164                ));
165
166                // f_{r,p} <= h_{r,p}
167                constraints.push(LinearConstraint::le(vec![(f_idx, 1.0), (h_idx, -1.0)], 0.0));
168
169                // f_{r,p} + a_{r,p} <= 1
170                constraints.push(LinearConstraint::le(vec![(f_idx, 1.0), (a_idx, 1.0)], 1.0));
171            }
172        }
173
174        // No augmentation allowed: sum f_{r,p} = 0
175        // This is the key difference from COMA: C1P requires zero flips
176        let mut flip_terms = Vec::new();
177        for r in 0..m {
178            for p in 0..k {
179                flip_terms.push((f_off + r * k + p, 1.0));
180            }
181        }
182        if !flip_terms.is_empty() {
183            constraints.push(LinearConstraint::eq(flip_terms, 0.0));
184        }
185
186        let target = ILP::new(num_vars, constraints, vec![], ObjectiveSense::Minimize);
187        ReductionCOSToILP {
188            target,
189            num_cols: n,
190        }
191    }
192}
193
194#[cfg(feature = "example-db")]
195pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
196    use crate::export::SolutionPair;
197    vec![crate::example_db::specs::RuleExampleSpec {
198        id: "consecutiveonessubmatrix_to_ilp",
199        build: || {
200            // Tucker matrix (3x4), K=3
201            let source = ConsecutiveOnesSubmatrix::new(
202                vec![
203                    vec![true, true, false, true],
204                    vec![true, false, true, true],
205                    vec![false, true, true, false],
206                ],
207                3,
208            );
209            let reduction: ReductionCOSToILP = ReduceTo::<ILP<bool>>::reduce_to(&source);
210            let ilp_solver = crate::solvers::ILPSolver::new();
211            let target_config = ilp_solver
212                .solve(reduction.target_problem())
213                .expect("ILP should be solvable");
214            let extracted = reduction.extract_solution(&target_config);
215            crate::example_db::specs::rule_example_with_witness::<_, ILP<bool>>(
216                source,
217                SolutionPair {
218                    source_config: extracted,
219                    target_config,
220                },
221            )
222        },
223    }]
224}
225
226#[cfg(test)]
227#[path = "../unit_tests/rules/consecutiveonessubmatrix_ilp.rs"]
228mod tests;
problemreductions/rules/consecutiveonessubmatrix_ilp.rs

problemreductions/rules/
consecutiveonessubmatrix_ilp.rs
