Skip to main content

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