Skip to main content

problemreductions/rules/
boundedcomponentspanningforest_ilp.rs

1//! Reduction from BoundedComponentSpanningForest to `ILP<i64, i64, Bounded>`.
2//!
3//! Assign every vertex to one of K components, bound weight, certify
4//! connectivity inside each used component via flow.
5//! See the paper entry for the full formulation.
6
7use crate::models::algebraic::{Bounded, IntegerVariable, LinearConstraint, ObjectiveSense, ILP};
8use crate::models::graph::BoundedComponentSpanningForest;
9use crate::reduction;
10use crate::rules::ilp_helpers::one_hot_decode_rows;
11use crate::rules::traits::{ReduceTo, ReductionResult};
12use crate::topology::{Graph, SimpleGraph};
13
14#[derive(Debug, Clone)]
15pub struct ReductionBCSFToILP {
16    target: ILP<i64, i64, Bounded>,
17    n: usize,
18    k: usize,
19}
20
21impl ReductionResult for ReductionBCSFToILP {
22    type Source = BoundedComponentSpanningForest<SimpleGraph, i64>;
23    type Target = ILP<i64, i64, Bounded>;
24
25    fn target_problem(&self) -> &ILP<i64, i64, Bounded> {
26        &self.target
27    }
28
29    /// One-hot decode: for each vertex v, output the unique component c with x_{v,c} = 1.
30    fn extract_solution(
31        &self,
32        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
33    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
34        crate::rules::traits::validate_target_witness(
35            self.target_problem(),
36            target_solution,
37            |value| value.value.is_some(),
38            "target ILP assignment is infeasible",
39        )?;
40
41        one_hot_decode_rows(target_solution, self.n, self.k, 0)
42    }
43}
44
45#[crate::aggregate_reduction(ilp_feasibility)]
46impl crate::rules::AggregateReductionResult for ReductionBCSFToILP {}
47
48// Assignment/weight/size/root/product/conservation rows contribute
49// <=16nk+6k terms; arc capacities and conservation add <=14mk.
50// Zero weights, n=0/1 capacity coefficients, and self-loop flows only cancel.
51#[reduction(transform = {
52    exact {
53        num_vars = "3 * num_vertices * max_components + 2 * max_components + 2 * num_edges * max_components",
54        num_constraints = "num_vertices + 5 * max_components + 6 * num_vertices * max_components + 6 * num_edges * max_components",
55    },
56    upper_bound {
57        max_constraint_magnitude_bits = "max_weight_bits + num_vertices",
58        num_nonzeros = "16 * num_vertices * max_components + 6 * max_components + 14 * num_edges * max_components",
59    },
60})]
61impl ReduceTo<ILP<i64, i64, Bounded>> for BoundedComponentSpanningForest<SimpleGraph, i64> {
62    type Result = ReductionBCSFToILP;
63
64    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
65        let n = self.num_vertices();
66        let edges = self.graph().edges();
67        let m = edges.len();
68        let k = self.max_components();
69
70        let x_idx = |v: usize, c: usize| -> usize { v * k + c };
71        let u_idx = |c: usize| -> usize { n * k + c };
72        let r_idx = |v: usize, c: usize| -> usize { n * k + k + v * k + c };
73        let s_idx = |c: usize| -> usize { 2 * n * k + k + c };
74        let b_idx = |v: usize, c: usize| -> usize { 2 * n * k + 2 * k + v * k + c };
75        let f_idx =
76            |i: usize, eta: usize, c: usize| -> usize { 3 * n * k + 2 * k + (i * 2 + eta) * k + c };
77
78        let num_vars = 3 * n * k + 2 * k + 2 * m * k;
79        let n_i64 = Self::exact_i64(n, "encoding the vertex count")?;
80        let mut constraints = Vec::new();
81        let weights = self.weights();
82        let max_weight = *self.max_weight();
83
84        // 1) Assignment: sum_c x_{v,c} = 1 for each vertex v
85        for v in 0..n {
86            let terms: Vec<(usize, i64)> = (0..k).map(|c| (x_idx(v, c), 1)).collect();
87            constraints.push(LinearConstraint::eq(terms, 1));
88        }
89
90        // 2) Weight bound: sum_v w_v * x_{v,c} <= B for each component c
91        for c in 0..k {
92            let terms: Vec<(usize, i64)> = weights
93                .iter()
94                .enumerate()
95                .map(|(vertex, &weight)| (x_idx(vertex, c), weight))
96                .collect();
97            constraints.push(LinearConstraint::le(terms, max_weight));
98        }
99
100        // 3) Size: s_c = sum_v x_{v,c}
101        for c in 0..k {
102            let mut terms: Vec<(usize, i64)> = vec![(s_idx(c), -1)];
103            for v in 0..n {
104                terms.push((x_idx(v, c), 1));
105            }
106            constraints.push(LinearConstraint::eq(terms, 0));
107        }
108
109        // 4) Nonempty indicator: u_c <= s_c and s_c <= n * u_c
110        for c in 0..k {
111            constraints.push(LinearConstraint::le(vec![(u_idx(c), 1), (s_idx(c), -1)], 0));
112            constraints.push(LinearConstraint::le(
113                vec![(s_idx(c), 1), (u_idx(c), -n_i64)],
114                0,
115            ));
116        }
117
118        // 5) Root selection: sum_v r_{v,c} = u_c and r_{v,c} <= x_{v,c}
119        for c in 0..k {
120            let mut terms: Vec<(usize, i64)> = (0..n).map(|v| (r_idx(v, c), 1)).collect();
121            terms.push((u_idx(c), -1));
122            constraints.push(LinearConstraint::eq(terms, 0));
123
124            for v in 0..n {
125                constraints.push(LinearConstraint::le(
126                    vec![(r_idx(v, c), 1), (x_idx(v, c), -1)],
127                    0,
128                ));
129            }
130        }
131
132        // 6) Product linearization: b_{v,c} = s_c * r_{v,c}
133        for v in 0..n {
134            for c in 0..k {
135                // b <= s_c
136                constraints.push(LinearConstraint::le(
137                    vec![(b_idx(v, c), 1), (s_idx(c), -1)],
138                    0,
139                ));
140                // b <= n * r
141                constraints.push(LinearConstraint::le(
142                    vec![(b_idx(v, c), 1), (r_idx(v, c), -n_i64)],
143                    0,
144                ));
145                // b >= s - n*(1-r) => b - s - n*r >= -n
146                constraints.push(LinearConstraint::ge(
147                    vec![(b_idx(v, c), 1), (s_idx(c), -1), (r_idx(v, c), -n_i64)],
148                    -n_i64,
149                ));
150                // b >= 0
151                constraints.push(LinearConstraint::ge(vec![(b_idx(v, c), 1)], 0));
152            }
153        }
154
155        // 7) Flow capacity: 0 <= f_{i,eta,c} <= (n-1)*x_{u_i,c} and <= (n-1)*x_{v_i,c}
156        let cap = n_i64 - 1;
157        for (i, &(u_e, v_e)) in edges.iter().enumerate() {
158            for eta in 0..2usize {
159                for c in 0..k {
160                    constraints.push(LinearConstraint::ge(vec![(f_idx(i, eta, c), 1)], 0));
161                    constraints.push(LinearConstraint::le(
162                        vec![(f_idx(i, eta, c), 1), (x_idx(u_e, c), -cap)],
163                        0,
164                    ));
165                    constraints.push(LinearConstraint::le(
166                        vec![(f_idx(i, eta, c), 1), (x_idx(v_e, c), -cap)],
167                        0,
168                    ));
169                }
170            }
171        }
172
173        // 8) Flow conservation: net_flow(v,c) = b_{v,c} - x_{v,c}
174        for v in 0..n {
175            for c in 0..k {
176                let mut terms: Vec<(usize, i64)> = Vec::new();
177
178                for (i, &(u_e, v_e)) in edges.iter().enumerate() {
179                    if u_e == v {
180                        terms.push((f_idx(i, 0, c), 1));
181                        terms.push((f_idx(i, 1, c), -1));
182                    }
183                    if v_e == v {
184                        terms.push((f_idx(i, 0, c), -1));
185                        terms.push((f_idx(i, 1, c), 1));
186                    }
187                }
188
189                terms.push((b_idx(v, c), -1));
190                terms.push((x_idx(v, c), 1));
191                constraints.push(LinearConstraint::eq(terms, 0));
192            }
193        }
194
195        let mut variables = vec![IntegerVariable::binary(); num_vars];
196        variables[2 * n * k + k..3 * n * k + 2 * k]
197            .fill(IntegerVariable::new(Some(0), Some(n_i64)).map_err(Self::target_construction)?);
198        variables[3 * n * k + 2 * k..].fill(
199            IntegerVariable::new(Some(0), Some(cap.max(0))).map_err(Self::target_construction)?,
200        );
201
202        let target = ILP::with_variables(variables, constraints, vec![], ObjectiveSense::Minimize)
203            .map_err(Self::target_construction)?;
204        Ok(ReductionBCSFToILP { target, n, k })
205    }
206}
207
208#[cfg(feature = "example-db")]
209pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
210    use crate::export::SolutionPair;
211    vec![crate::example_db::specs::RuleExampleSpec {
212        id: "boundedcomponentspanningforest_to_ilp",
213        build: || {
214            let source = BoundedComponentSpanningForest::new(
215                SimpleGraph::new(4, vec![(0, 1), (1, 2), (2, 3)]),
216                vec![1, 2, 2, 1],
217                2,
218                4,
219            );
220            let reduction: ReductionBCSFToILP =
221                crate::rules::ReduceTo::<ILP<i64, i64, Bounded>>::reduce_to(&source)
222                    .expect("reduction should succeed");
223            let ilp_sol = crate::solvers::ILPSolver::new()
224                .solve(reduction.target_problem())
225                .expect("ILP should be solvable");
226            let extracted = reduction.extract_solution(&ilp_sol).unwrap();
227            crate::example_db::specs::rule_example_with_witness::<_, ILP<i64, i64, Bounded>>(
228                source,
229                SolutionPair {
230                    source_config: serde_json::json!(extracted),
231                    target_config: serde_json::json!(ilp_sol),
232                },
233            )
234        },
235    }]
236}
237
238#[cfg(test)]
239#[path = "../unit_tests/rules/boundedcomponentspanningforest_ilp.rs"]
240mod tests;