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(
26        &self,
27        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
28    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
29        crate::rules::traits::validate_target_witness(
30            self.target_problem(),
31            target_solution,
32            |value| value.value.is_some(),
33            "target ILP assignment is infeasible",
34        )?;
35
36        Ok({
37            // Output the selection bits s_c (first num_cols variables)
38            target_solution[..self.num_cols]
39                .iter()
40                .map(|&value| value == 1)
41                .collect()
42        })
43    }
44}
45
46#[crate::aggregate_reduction(ilp_feasibility)]
47impl crate::rules::AggregateReductionResult for ReductionCOSToILP {}
48
49#[reduction(transform = upper_bound {
50    max_constraint_magnitude_bits = "bound + num_cols + 1",
51    num_vars = "num_cols + num_cols * bound + 5 * num_rows * bound",
52    num_constraints = "2 + num_cols + bound + 3 * num_rows + 8 * num_rows * bound",
53    num_nonzeros = "(num_cols + num_cols * bound + 5 * num_rows * bound) * (2 + num_cols + bound + 3 * num_rows + 8 * num_rows * bound)",
54})]
55impl ReduceTo<ILP<bool>> for ConsecutiveOnesSubmatrix {
56    type Result = ReductionCOSToILP;
57
58    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
59        let m = self.num_rows();
60        let n = self.num_cols();
61        let k = self.bound() as usize;
62
63        // Variable layout (all binary):
64        // s_c: n vars at [0, n)  — column selection
65        // x_{c,p}: n*K vars at [n, n + n*K)  — column c placed at position p in [0..K)
66        // a_{r,p}: m*K vars at [n + n*K, n + n*K + m*K)  — value at row r, position p
67        // l_{r,p}: m*K vars  — left boundary
68        // u_{r,p}: m*K vars  — right boundary
69        // h_{r,p}: m*K vars  — inside interval
70        // f_{r,p}: m*K vars  — flip indicator (not used for budget, but needed for C1P)
71        let s_off = 0;
72        let x_off = n;
73        let a_off = n + n * k;
74        let l_off = a_off + m * k;
75        let u_off = l_off + m * k;
76        let h_off = u_off + m * k;
77        let f_off = h_off + m * k;
78        let num_vars = f_off + m * k;
79
80        let mut constraints = Vec::new();
81
82        // sum_c s_c = K
83        let s_terms: Vec<(usize, i64)> = (0..n).map(|c| (s_off + c, 1)).collect();
84        constraints.push(LinearConstraint::eq(
85            s_terms,
86            <Self as ReduceTo<ILP<bool>>>::exact_i64(k, "encoding the selected column count")?,
87        ));
88
89        // sum_p x_{c,p} = s_c for all c
90        for c in 0..n {
91            let mut terms: Vec<(usize, i64)> = (0..k).map(|p| (x_off + c * k + p, 1)).collect();
92            terms.push((s_off + c, -1));
93            constraints.push(LinearConstraint::eq(terms, 0));
94        }
95
96        // sum_c x_{c,p} = 1 for all p in {0, ..., K-1}
97        for p in 0..k {
98            let terms: Vec<(usize, i64)> = (0..n).map(|c| (x_off + c * k + p, 1)).collect();
99            constraints.push(LinearConstraint::eq(terms, 1));
100        }
101
102        // a_{r,p} = sum_c A_{r,c} * x_{c,p}
103        for r in 0..m {
104            for p in 0..k {
105                let a_idx = a_off + r * k + p;
106                let mut terms = vec![(a_idx, 1)];
107                for c in 0..n {
108                    if self.matrix()[r][c] {
109                        terms.push((x_off + c * k + p, -1));
110                    }
111                }
112                constraints.push(LinearConstraint::eq(terms, 0));
113            }
114        }
115
116        // C1P interval constraints on the K-position permuted submatrix
117        for r in 0..m {
118            // A row has either no interval (all selected entries are zero), or
119            // exactly one pair of left and right boundaries. Existing a <= h
120            // constraints force the latter whenever a selected entry is one.
121            let l_terms: Vec<(usize, i64)> = (0..k).map(|p| (l_off + r * k + p, 1)).collect();
122            constraints.push(LinearConstraint::le(l_terms.clone(), 1));
123
124            let u_terms: Vec<(usize, i64)> = (0..k).map(|p| (u_off + r * k + p, 1)).collect();
125            let mut boundary_count_terms = l_terms;
126            boundary_count_terms.extend(u_terms.into_iter().map(|(index, _)| (index, -1)));
127            constraints.push(LinearConstraint::eq(boundary_count_terms, 0));
128
129            // The left boundary cannot occur after the right boundary
130            if k > 0 {
131                let mut order_terms = Vec::new();
132                for p in 0..k {
133                    let p_i64 = <Self as ReduceTo<ILP<bool>>>::exact_i64(
134                        p,
135                        "encoding a selected-column position",
136                    )?;
137                    order_terms.push((l_off + r * k + p, p_i64));
138                    order_terms.push((u_off + r * k + p, -p_i64));
139                }
140                constraints.push(LinearConstraint::le(order_terms, 0));
141            }
142
143            for p in 0..k {
144                let h_idx = h_off + r * k + p;
145                let a_idx = a_off + r * k + p;
146                let f_idx = f_off + r * k + p;
147
148                // h_{r,p} <= sum_{q=0}^{p} l_{r,q}
149                let mut h_le_l = vec![(h_idx, 1)];
150                for q in 0..=p {
151                    h_le_l.push((l_off + r * k + q, -1));
152                }
153                constraints.push(LinearConstraint::le(h_le_l, 0));
154
155                // h_{r,p} <= sum_{q=p}^{K-1} u_{r,q}
156                let mut h_le_u = vec![(h_idx, 1)];
157                for q in p..k {
158                    h_le_u.push((u_off + r * k + q, -1));
159                }
160                constraints.push(LinearConstraint::le(h_le_u, 0));
161
162                // h_{r,p} >= sum_{q=0}^{p} l_{r,q} + sum_{q=p}^{K-1} u_{r,q} - 1
163                let mut h_ge_terms = vec![(h_idx, 1)];
164                for q in 0..=p {
165                    h_ge_terms.push((l_off + r * k + q, -1));
166                }
167                for q in p..k {
168                    h_ge_terms.push((u_off + r * k + q, -1));
169                }
170                constraints.push(LinearConstraint::ge(h_ge_terms, -1));
171
172                // a_{r,p} <= h_{r,p}  — every 1 must be inside the interval
173                constraints.push(LinearConstraint::le(vec![(a_idx, 1), (h_idx, -1)], 0));
174
175                // For C1P (no augmentation): the interval must exactly cover the 1s
176                // h_{r,p} <= a_{r,p} + f_{r,p} — position inside interval but 0 costs a flip
177                constraints.push(LinearConstraint::le(
178                    vec![(h_idx, 1), (a_idx, -1), (f_idx, -1)],
179                    0,
180                ));
181
182                // f_{r,p} <= h_{r,p}
183                constraints.push(LinearConstraint::le(vec![(f_idx, 1), (h_idx, -1)], 0));
184
185                // f_{r,p} + a_{r,p} <= 1
186                constraints.push(LinearConstraint::le(vec![(f_idx, 1), (a_idx, 1)], 1));
187            }
188        }
189
190        // No augmentation allowed: sum f_{r,p} = 0
191        // This is the key difference from COMA: C1P requires zero flips
192        let mut flip_terms = Vec::new();
193        for r in 0..m {
194            for p in 0..k {
195                flip_terms.push((f_off + r * k + p, 1));
196            }
197        }
198        if !flip_terms.is_empty() {
199            constraints.push(LinearConstraint::eq(flip_terms, 0));
200        }
201
202        let target = ILP::new(num_vars, constraints, vec![], ObjectiveSense::Minimize)
203            .map_err(Self::target_construction)?;
204        Ok(ReductionCOSToILP {
205            target,
206            num_cols: n,
207        })
208    }
209}
210
211#[cfg(feature = "example-db")]
212pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
213    use crate::export::SolutionPair;
214    vec![crate::example_db::specs::RuleExampleSpec {
215        id: "consecutiveonessubmatrix_to_ilp",
216        build: || {
217            // Tucker matrix (3x4), K=3
218            let source = ConsecutiveOnesSubmatrix::new(
219                vec![
220                    vec![true, true, false, true],
221                    vec![true, false, true, true],
222                    vec![false, true, true, false],
223                ],
224                3,
225            );
226            let reduction: ReductionCOSToILP =
227                ReduceTo::<ILP<bool>>::reduce_to(&source).expect("reduction should succeed");
228            let ilp_solver = crate::solvers::ILPSolver::new();
229            let target_config = ilp_solver
230                .solve(reduction.target_problem())
231                .expect("ILP should be solvable");
232            let extracted = reduction.extract_solution(&target_config).unwrap();
233            crate::example_db::specs::rule_example_with_witness::<_, ILP<bool>>(
234                source,
235                SolutionPair {
236                    source_config: serde_json::json!(extracted),
237                    target_config: serde_json::to_value(target_config)
238                        .expect("solution serialization must succeed"),
239                },
240            )
241        },
242    }]
243}
244
245#[cfg(test)]
246#[path = "../unit_tests/rules/consecutiveonessubmatrix_ilp.rs"]
247mod tests;
problemreductions/rules/consecutiveonessubmatrix_ilp.rs

problemreductions/rules/
consecutiveonessubmatrix_ilp.rs
