problemreductions/rules/
sumofsquarespartition_ilp.rs

1//! Reduction from SumOfSquaresPartition to ILP (Integer Linear Programming).
2//!
3//! The objective Σ_g (Σ_{i ∈ g} s_i)^2 is quadratic, so we linearize using McCormick:
4//!
5//! Variables:
6//! - x_{i,g}: binary, element i in group g (index: i * K + g)
7//! - z_{i,j,g}: binary product for x_{i,g} * x_{j,g} (index: n*K + (i*n + j) * K + g)
8//!
9//! Constraints:
10//! - Σ_g x_{i,g} = 1 for each element i (assignment)
11//! - McCormick for each (i,j,g):
12//!   z_{i,j,g} ≤ x_{i,g}, z_{i,j,g} ≤ x_{j,g}, z_{i,j,g} ≥ x_{i,g} + x_{j,g} - 1
13//!
14//! Note: Σ_g (Σ_i s_i * x_{i,g})^2 = Σ_g Σ_{i,j} s_i * s_j * x_{i,g} * x_{j,g}
15//!       which equals Σ_g Σ_{i,j} s_i * s_j * z_{i,j,g} after linearization.
16//!
17//! Objective: Minimize Σ_g Σ_{i,j} s_i * s_j * z_{i,j,g}
18
19use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
20use crate::models::misc::SumOfSquaresPartition;
21use crate::reduction;
22use crate::rules::ilp_helpers::mccormick_product;
23use crate::rules::traits::{ReduceTo, ReductionResult};
24
25/// Result of reducing SumOfSquaresPartition to ILP.
26///
27/// Variable layout:
28/// - x_{i,g} at index i * K + g  (i ∈ 0..n, g ∈ 0..K)
29/// - z_{i,j,g} at index n*K + (i*n + j) * K + g  (i,j ∈ 0..n, g ∈ 0..K)
30///
31/// Total: n*K + n^2*K variables.
32#[derive(Debug, Clone)]
33pub struct ReductionSSPToILP {
34    target: ILP<bool>,
35    num_elements: usize,
36    num_groups: usize,
37}
38
39impl ReductionSSPToILP {
40    fn x_var(&self, i: usize, g: usize) -> usize {
41        i * self.num_groups + g
42    }
43
44    fn z_var(&self, i: usize, j: usize, g: usize) -> usize {
45        let n = self.num_elements;
46        let k = self.num_groups;
47        n * k + (i * n + j) * k + g
48    }
49}
50
51impl ReductionResult for ReductionSSPToILP {
52    type Source = SumOfSquaresPartition;
53    type Target = ILP<bool>;
54
55    fn target_problem(&self) -> &ILP<bool> {
56        &self.target
57    }
58
59    /// Extract solution: for each element i, find the unique group g where x_{i,g} = 1.
60    fn extract_solution(
61        &self,
62        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
63    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
64        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
65
66        crate::rules::ilp_helpers::one_hot_decode_rows(
67            target_solution,
68            self.num_elements,
69            self.num_groups,
70            0,
71        )
72    }
73}
74
75#[reduction(transform = {
76    exact {
77        num_vars = "num_elements * num_groups + num_elements^2 * num_groups",
78        num_constraints = "num_elements + 3 * num_elements^2 * num_groups",
79        num_nonzeros = "7 * num_elements^2 * num_groups",
80    },
81    upper_bound {
82        max_constraint_magnitude_bits = "2",
83    },
84})]
85impl ReduceTo<ILP<bool>> for SumOfSquaresPartition {
86    type Result = ReductionSSPToILP;
87
88    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
89        let n = self.num_elements();
90        let k = self.num_groups();
91        let num_vars = n * k + n * n * k;
92
93        let result = ReductionSSPToILP {
94            target: ILP::empty(),
95            num_elements: n,
96            num_groups: k,
97        };
98
99        let mut constraints = Vec::new();
100
101        // Assignment constraints: for each element i, Σ_g x_{i,g} = 1
102        for i in 0..n {
103            let terms: Vec<(usize, i64)> = (0..k).map(|g| (result.x_var(i, g), 1)).collect();
104            constraints.push(LinearConstraint::eq(terms, 1));
105        }
106
107        // McCormick linearization for z_{i,j,g} = x_{i,g} * x_{j,g}
108        for i in 0..n {
109            for j in 0..n {
110                for g in 0..k {
111                    let z = result.z_var(i, j, g);
112                    let xi = result.x_var(i, g);
113                    let xj = result.x_var(j, g);
114
115                    constraints.extend(mccormick_product(z, xi, xj));
116                }
117            }
118        }
119
120        // Objective: Minimize Σ_g Σ_{i,j} s_i * s_j * z_{i,j,g}
121        let sizes = self.sizes();
122        let mut objective: Vec<(usize, i64)> = Vec::new();
123        for i in 0..n {
124            for j in 0..n {
125                for g in 0..k {
126                    let product = sizes[i].checked_mul(sizes[j]).ok_or_else(|| {
127                        crate::rules::ReductionError::integer_overflow::<
128                            SumOfSquaresPartition,
129                            ILP<bool>,
130                        >("multiplying two partition element sizes")
131                    })?;
132                    let coeff = product;
133                    if coeff != 0 {
134                        objective.push((result.z_var(i, j, g), coeff));
135                    }
136                }
137            }
138        }
139
140        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Minimize)
141            .map_err(Self::target_construction)?;
142
143        Ok(ReductionSSPToILP {
144            target,
145            num_elements: n,
146            num_groups: k,
147        })
148    }
149}
150
151#[cfg(feature = "example-db")]
152pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
153    vec![crate::example_db::specs::RuleExampleSpec {
154        id: "sumofsquarespartition_to_ilp",
155        build: || {
156            // 4 elements [1, 2, 3, 4], K=2 groups
157            // Group {1,4}: sum=5, Group {2,3}: sum=5 → sos = 25+25 = 50
158            let source = SumOfSquaresPartition::new(vec![1, 2, 3, 4], 2);
159            crate::example_db::specs::rule_example_via_ilp::<_, bool>(source)
160        },
161    }]
162}
163
164#[cfg(test)]
165#[path = "../unit_tests/rules/sumofsquarespartition_ilp.rs"]
166mod tests;
problemreductions/rules/sumofsquarespartition_ilp.rs

problemreductions/rules/
sumofsquarespartition_ilp.rs
