Skip to main content

problemreductions/rules/
rootedtreestorageassignment_ilp.rs

1//! Reduction from RootedTreeStorageAssignment to ILP (Integer Linear Programming).
2//!
3//! Uses parent indicators p_{v,u}, depth variables d_v, ancestor indicators
4//! a_{u,v}, transitive-closure helpers h_{u,v,w}, and per-subset gadgets
5//! (top/bottom selectors, pair selectors, endpoint depths, extension costs).
6
7use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
8use crate::models::set::RootedTreeStorageAssignment;
9use crate::reduction;
10use crate::rules::traits::{ReduceTo, ReductionResult};
11
12// Index helpers
13
14fn idx_p(n: usize, v: usize, u: usize) -> usize {
15    v * n + u
16}
17
18fn idx_d(n: usize, v: usize) -> usize {
19    n * n + v
20}
21
22fn idx_a(n: usize, u: usize, v: usize) -> usize {
23    n * n + n + u * n + v
24}
25
26fn idx_h(n: usize, u: usize, v: usize, w: usize) -> usize {
27    2 * n * n + n + (u * n + v) * n + w
28}
29
30fn idx_t(n: usize, r: usize, s: usize, u: usize) -> usize {
31    let _ = r;
32    n * n * n + 2 * n * n + n + s * n + u
33}
34
35fn idx_b(n: usize, r: usize, s: usize, v: usize) -> usize {
36    n * n * n + 2 * n * n + n + r * n + s * n + v
37}
38
39fn idx_m(n: usize, r: usize, s: usize, u: usize, v: usize) -> usize {
40    n * n * n + 2 * n * n + n + 2 * r * n + s * n * n + u * n + v
41}
42
43fn idx_big_t(n: usize, r: usize, s: usize) -> usize {
44    n * n * n + 2 * n * n + n + 2 * r * n + r * n * n + s
45}
46
47fn idx_big_b(n: usize, r: usize, s: usize) -> usize {
48    n * n * n + 2 * n * n + n + 2 * r * n + r * n * n + r + s
49}
50
51fn idx_c(n: usize, r: usize, s: usize) -> usize {
52    n * n * n + 2 * n * n + n + 2 * r * n + r * n * n + 2 * r + s
53}
54
55fn total_vars(n: usize, r: usize) -> usize {
56    n * n * n + 2 * n * n + n + r * (n * n + 2 * n + 3)
57}
58
59#[derive(Debug, Clone)]
60pub struct ReductionRTSAToILP {
61    target: ILP<i32>,
62    n: usize,
63}
64
65impl ReductionResult for ReductionRTSAToILP {
66    type Source = RootedTreeStorageAssignment;
67    type Target = ILP<i32>;
68
69    fn target_problem(&self) -> &ILP<i32> {
70        &self.target
71    }
72
73    /// Decode parent array from one-hot parent indicators p_{v,u}.
74    fn extract_solution(&self, target_solution: &[usize]) -> Vec<usize> {
75        let n = self.n;
76        (0..n)
77            .map(|v| {
78                (0..n)
79                    .find(|&u| target_solution[idx_p(n, v, u)] == 1)
80                    .unwrap_or(v)
81            })
82            .collect()
83    }
84}
85
86#[reduction(
87    overhead = {
88        num_vars = "universe_size * universe_size * universe_size + 2 * universe_size * universe_size + universe_size + num_subsets * (universe_size * universe_size + 2 * universe_size + 3)",
89        num_constraints = "universe_size * universe_size * universe_size + universe_size * universe_size + universe_size * universe_size + num_subsets * universe_size * universe_size",
90    }
91)]
92impl ReduceTo<ILP<i32>> for RootedTreeStorageAssignment {
93    type Result = ReductionRTSAToILP;
94
95    fn reduce_to(&self) -> Self::Result {
96        let n = self.universe_size();
97        let subsets = self.subsets();
98        let bound = self.bound();
99
100        // Nontrivial subsets (size >= 2)
101        let nontrivial: Vec<usize> = (0..subsets.len())
102            .filter(|&k| subsets[k].len() >= 2)
103            .collect();
104        let r = nontrivial.len();
105
106        if n == 0 {
107            return ReductionRTSAToILP {
108                target: ILP::new(0, vec![], vec![], ObjectiveSense::Minimize),
109                n,
110            };
111        }
112
113        let nv = total_vars(n, r);
114        let big_m = n as f64;
115        let big_m_depth = (n - 1) as f64;
116
117        let mut constraints = Vec::new();
118
119        // === Rooted-tree constraints ===
120
121        // Σ_u p_{v,u} = 1  ∀ v
122        for v in 0..n {
123            let terms: Vec<(usize, f64)> = (0..n).map(|u| (idx_p(n, v, u), 1.0)).collect();
124            constraints.push(LinearConstraint::eq(terms, 1.0));
125        }
126
127        // Σ_v p_{v,v} = 1 (exactly one root)
128        let root_terms: Vec<(usize, f64)> = (0..n).map(|v| (idx_p(n, v, v), 1.0)).collect();
129        constraints.push(LinearConstraint::eq(root_terms, 1.0));
130
131        // p_{v,u} binary: upper bound p_{v,u} <= 1
132        for v in 0..n {
133            for u in 0..n {
134                constraints.push(LinearConstraint::le(vec![(idx_p(n, v, u), 1.0)], 1.0));
135            }
136        }
137
138        // d_v <= (n-1)(1 - p_{v,v})  ∀ v  (root has depth 0)
139        for v in 0..n {
140            constraints.push(LinearConstraint::le(
141                vec![(idx_d(n, v), 1.0), (idx_p(n, v, v), big_m_depth)],
142                big_m_depth,
143            ));
144        }
145
146        // d_v >= 0  ∀ v
147        for v in 0..n {
148            constraints.push(LinearConstraint::ge(vec![(idx_d(n, v), 1.0)], 0.0));
149        }
150
151        // d_v <= n-1  ∀ v
152        for v in 0..n {
153            constraints.push(LinearConstraint::le(vec![(idx_d(n, v), 1.0)], big_m_depth));
154        }
155
156        // For u != v: d_v - d_u >= 1 - n(1 - p_{v,u})
157        //             d_v - d_u <= 1 + n(1 - p_{v,u})
158        for v in 0..n {
159            for u in 0..n {
160                if u != v {
161                    // d_v - d_u + n*p_{v,u} >= 1 - n + n = 1
162                    // => d_v - d_u + n*p_{v,u} >= 1 - n*(1 - p_{v,u})
163                    // Rewrite: d_v - d_u + n*p_{v,u} >= 1 - n + n*p_{v,u} ... no.
164                    // Original: d_v - d_u >= 1 - n(1 - p_{v,u})
165                    // => d_v - d_u + n - n*p_{v,u} >= 1
166                    // => d_v - d_u - n*p_{v,u} >= 1 - n
167                    constraints.push(LinearConstraint::ge(
168                        vec![
169                            (idx_d(n, v), 1.0),
170                            (idx_d(n, u), -1.0),
171                            (idx_p(n, v, u), -big_m),
172                        ],
173                        1.0 - big_m,
174                    ));
175
176                    // d_v - d_u <= 1 + n(1 - p_{v,u})
177                    // => d_v - d_u - n + n*p_{v,u} <= 1
178                    // => d_v - d_u + n*p_{v,u} <= 1 + n
179                    constraints.push(LinearConstraint::le(
180                        vec![
181                            (idx_d(n, v), 1.0),
182                            (idx_d(n, u), -1.0),
183                            (idx_p(n, v, u), big_m),
184                        ],
185                        1.0 + big_m,
186                    ));
187                }
188            }
189        }
190
191        // === Ancestor relation ===
192
193        // a_{v,v} = 1  ∀ v
194        for v in 0..n {
195            constraints.push(LinearConstraint::eq(vec![(idx_a(n, v, v), 1.0)], 1.0));
196        }
197
198        // h_{u,v,v} = 0  ∀ u,v
199        for u in 0..n {
200            for v in 0..n {
201                constraints.push(LinearConstraint::eq(vec![(idx_h(n, u, v, v), 1.0)], 0.0));
202            }
203        }
204
205        // For u != v: a_{u,v} = Σ_w h_{u,v,w}
206        for u in 0..n {
207            for v in 0..n {
208                if u != v {
209                    let mut terms = vec![(idx_a(n, u, v), -1.0)];
210                    for w in 0..n {
211                        terms.push((idx_h(n, u, v, w), 1.0));
212                    }
213                    constraints.push(LinearConstraint::eq(terms, 0.0));
214                }
215            }
216        }
217
218        // h_{u,v,w} <= p_{v,w}  ∀ u,v,w with w != v
219        // h_{u,v,w} <= a_{u,w}  ∀ u,v,w with w != v
220        // h_{u,v,w} >= p_{v,w} + a_{u,w} - 1  ∀ u,v,w with w != v
221        for u in 0..n {
222            for v in 0..n {
223                for w in 0..n {
224                    if w != v {
225                        constraints.push(LinearConstraint::le(
226                            vec![(idx_h(n, u, v, w), 1.0), (idx_p(n, v, w), -1.0)],
227                            0.0,
228                        ));
229                        constraints.push(LinearConstraint::le(
230                            vec![(idx_h(n, u, v, w), 1.0), (idx_a(n, u, w), -1.0)],
231                            0.0,
232                        ));
233                        constraints.push(LinearConstraint::ge(
234                            vec![
235                                (idx_h(n, u, v, w), 1.0),
236                                (idx_p(n, v, w), -1.0),
237                                (idx_a(n, u, w), -1.0),
238                            ],
239                            -1.0,
240                        ));
241                    }
242                }
243            }
244        }
245
246        // Binary bounds for a, h
247        for u in 0..n {
248            for v in 0..n {
249                constraints.push(LinearConstraint::le(vec![(idx_a(n, u, v), 1.0)], 1.0));
250                for w in 0..n {
251                    constraints.push(LinearConstraint::le(vec![(idx_h(n, u, v, w), 1.0)], 1.0));
252                }
253            }
254        }
255
256        // === Subset gadgets ===
257        for (s, &orig_k) in nontrivial.iter().enumerate() {
258            let subset = &subsets[orig_k];
259            let subset_size = subset.len();
260
261            // Top selectors: Σ_{u ∈ S} t_{s,u} = 1, t_{s,u} = 0 for u ∉ S
262            let top_terms: Vec<(usize, f64)> =
263                subset.iter().map(|&u| (idx_t(n, r, s, u), 1.0)).collect();
264            constraints.push(LinearConstraint::eq(top_terms, 1.0));
265            for u in 0..n {
266                if !subset.contains(&u) {
267                    constraints.push(LinearConstraint::eq(vec![(idx_t(n, r, s, u), 1.0)], 0.0));
268                }
269                // Binary bound
270                constraints.push(LinearConstraint::le(vec![(idx_t(n, r, s, u), 1.0)], 1.0));
271            }
272
273            // Bottom selectors: Σ_{v ∈ S} b_{s,v} = 1, b_{s,v} = 0 for v ∉ S
274            let bot_terms: Vec<(usize, f64)> =
275                subset.iter().map(|&v| (idx_b(n, r, s, v), 1.0)).collect();
276            constraints.push(LinearConstraint::eq(bot_terms, 1.0));
277            for v in 0..n {
278                if !subset.contains(&v) {
279                    constraints.push(LinearConstraint::eq(vec![(idx_b(n, r, s, v), 1.0)], 0.0));
280                }
281                constraints.push(LinearConstraint::le(vec![(idx_b(n, r, s, v), 1.0)], 1.0));
282            }
283
284            // Pair selectors (McCormick): m_{s,u,v} = t_{s,u} * b_{s,v}
285            for u in 0..n {
286                for v in 0..n {
287                    constraints.push(LinearConstraint::le(
288                        vec![(idx_m(n, r, s, u, v), 1.0), (idx_t(n, r, s, u), -1.0)],
289                        0.0,
290                    ));
291                    constraints.push(LinearConstraint::le(
292                        vec![(idx_m(n, r, s, u, v), 1.0), (idx_b(n, r, s, v), -1.0)],
293                        0.0,
294                    ));
295                    constraints.push(LinearConstraint::ge(
296                        vec![
297                            (idx_m(n, r, s, u, v), 1.0),
298                            (idx_t(n, r, s, u), -1.0),
299                            (idx_b(n, r, s, v), -1.0),
300                        ],
301                        -1.0,
302                    ));
303                    constraints.push(LinearConstraint::le(vec![(idx_m(n, r, s, u, v), 1.0)], 1.0));
304                }
305            }
306
307            // Path condition: m_{s,u,v} <= a_{u,v} (top is ancestor of bottom)
308            for u in 0..n {
309                for v in 0..n {
310                    constraints.push(LinearConstraint::le(
311                        vec![(idx_m(n, r, s, u, v), 1.0), (idx_a(n, u, v), -1.0)],
312                        0.0,
313                    ));
314                }
315            }
316
317            // Every subset element w lies on the chain:
318            // m_{s,u,v} <= a_{u,w} and m_{s,u,v} <= a_{w,v}  ∀ w ∈ S, u, v
319            for &w in subset {
320                for u in 0..n {
321                    for v in 0..n {
322                        constraints.push(LinearConstraint::le(
323                            vec![(idx_m(n, r, s, u, v), 1.0), (idx_a(n, u, w), -1.0)],
324                            0.0,
325                        ));
326                        constraints.push(LinearConstraint::le(
327                            vec![(idx_m(n, r, s, u, v), 1.0), (idx_a(n, w, v), -1.0)],
328                            0.0,
329                        ));
330                    }
331                }
332            }
333
334            // Endpoint depths: T_s, B_s
335            // T_s - d_u <= (n-1)(1 - t_{s,u})  and  d_u - T_s <= (n-1)(1 - t_{s,u})
336            for &u in subset {
337                constraints.push(LinearConstraint::le(
338                    vec![
339                        (idx_big_t(n, r, s), 1.0),
340                        (idx_d(n, u), -1.0),
341                        (idx_t(n, r, s, u), big_m_depth),
342                    ],
343                    big_m_depth,
344                ));
345                constraints.push(LinearConstraint::le(
346                    vec![
347                        (idx_d(n, u), 1.0),
348                        (idx_big_t(n, r, s), -1.0),
349                        (idx_t(n, r, s, u), big_m_depth),
350                    ],
351                    big_m_depth,
352                ));
353            }
354            // B_s - d_v <= (n-1)(1 - b_{s,v})  and  d_v - B_s <= (n-1)(1 - b_{s,v})
355            for &v in subset {
356                constraints.push(LinearConstraint::le(
357                    vec![
358                        (idx_big_b(n, r, s), 1.0),
359                        (idx_d(n, v), -1.0),
360                        (idx_b(n, r, s, v), big_m_depth),
361                    ],
362                    big_m_depth,
363                ));
364                constraints.push(LinearConstraint::le(
365                    vec![
366                        (idx_d(n, v), 1.0),
367                        (idx_big_b(n, r, s), -1.0),
368                        (idx_b(n, r, s, v), big_m_depth),
369                    ],
370                    big_m_depth,
371                ));
372            }
373
374            // Depth bounds for T_s, B_s
375            constraints.push(LinearConstraint::ge(vec![(idx_big_t(n, r, s), 1.0)], 0.0));
376            constraints.push(LinearConstraint::le(
377                vec![(idx_big_t(n, r, s), 1.0)],
378                big_m_depth,
379            ));
380            constraints.push(LinearConstraint::ge(vec![(idx_big_b(n, r, s), 1.0)], 0.0));
381            constraints.push(LinearConstraint::le(
382                vec![(idx_big_b(n, r, s), 1.0)],
383                big_m_depth,
384            ));
385
386            // Extension cost: c_s = B_s - T_s + 1 - |S|
387            // => c_s - B_s + T_s = 1 - |S|
388            constraints.push(LinearConstraint::eq(
389                vec![
390                    (idx_c(n, r, s), 1.0),
391                    (idx_big_b(n, r, s), -1.0),
392                    (idx_big_t(n, r, s), 1.0),
393                ],
394                1.0 - subset_size as f64,
395            ));
396
397            // c_s >= 0
398            constraints.push(LinearConstraint::ge(vec![(idx_c(n, r, s), 1.0)], 0.0));
399        }
400
401        // Total cost bound: Σ c_s <= K
402        if r > 0 {
403            let cost_terms: Vec<(usize, f64)> = (0..r).map(|s| (idx_c(n, r, s), 1.0)).collect();
404            constraints.push(LinearConstraint::le(cost_terms, bound as f64));
405        }
406
407        let target = ILP::new(nv, constraints, vec![], ObjectiveSense::Minimize);
408        ReductionRTSAToILP { target, n }
409    }
410}
411
412#[cfg(feature = "example-db")]
413pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
414    use crate::export::SolutionPair;
415    vec![crate::example_db::specs::RuleExampleSpec {
416        id: "rootedtreestorageassignment_to_ilp",
417        build: || {
418            let source = RootedTreeStorageAssignment::new(3, vec![vec![0, 1], vec![1, 2]], 1);
419            let reduction: ReductionRTSAToILP = ReduceTo::<ILP<i32>>::reduce_to(&source);
420            let target_config = {
421                let ilp_solver = crate::solvers::ILPSolver::new();
422                ilp_solver
423                    .solve(reduction.target_problem())
424                    .expect("ILP should be solvable")
425            };
426            let source_config = reduction.extract_solution(&target_config);
427            crate::example_db::specs::rule_example_with_witness::<_, ILP<i32>>(
428                source,
429                SolutionPair {
430                    source_config,
431                    target_config,
432                },
433            )
434        },
435    }]
436}
437
438#[cfg(test)]
439#[path = "../unit_tests/rules/rootedtreestorageassignment_ilp.rs"]
440mod tests;