problemreductions/rules/
maximumcontactmapoverlap_ilp.rs

1//! Reduction from MaximumContactMapOverlap to ILP (Integer Linear Programming).
2//!
3//! Binary mapping variables `x_(i,j)` indicate that residue `i in V_1` is
4//! aligned to residue `j in V_2`. Row and column inequalities encode a partial
5//! injective alignment. Order-preservation is enforced by forbidding crossings
6//! and equal-image matches: for every `i < k in V_1` and every `j >= l in V_2`,
7//! we add `x_(i,j) + x_(k,l) <= 1`. For every pair of contacts
8//! `({i,k} in E_1, {j,l} in E_2)` with `i < k` and `j < l` we introduce a
9//! binary `y_(i,k,j,l)` linked by `y <= x_(i,j)` and `y <= x_(k,l)`. The ILP
10//! objective is `max sum y_(i,k,j,l)`, which equals the number of preserved
11//! contacts under the alignment.
12//!
13//! This is a direct ILP rendering of the polyhedral formulation studied by
14//! Andonov, Malod-Dognin, and Yanev (J. Comput. Biol., 2011) and by
15//! Xie and Sahinidis (J. Comput. Biol., 2007).
16
17use crate::models::algebraic::{LinearConstraint, ObjectiveSense, ILP};
18use crate::models::graph::MaximumContactMapOverlap;
19use crate::reduction;
20use crate::rules::traits::{ReduceTo, ReductionResult};
21
22/// Result of reducing MaximumContactMapOverlap to ILP.
23///
24/// Variable layout (all binary):
25/// - `x_(i,j)` at index `i * n2 + j` for `i in V_1`, `j in V_2`
26/// - `y_(i,k,j,l)` for each contact pair from `E_1 x E_2` (with `i < k` and
27///   `j < l` enforced by the contact canonicalization), indexed sequentially
28///   after the `x` block in the order they are enumerated by the constructor.
29#[derive(Debug, Clone)]
30pub struct ReductionCMOToILP {
31    target: ILP<bool>,
32    num_vertices_1: usize,
33    num_vertices_2: usize,
34}
35
36impl ReductionResult for ReductionCMOToILP {
37    type Source = MaximumContactMapOverlap;
38    type Target = ILP<bool>;
39
40    fn target_problem(&self) -> &ILP<bool> {
41        &self.target
42    }
43
44    /// Extract the CMO configuration from the ILP assignment.
45    ///
46    /// For each source residue `i in V_1`, find the unique `j` with
47    /// `x_(i,j) = 1` and encode it as `j + 1` (CMO's `bot` is `0`); if no
48    /// `x_(i,*)` is selected, the residue is left unmatched (`0`).
49    fn extract_solution(
50        &self,
51        target_solution: &<Self::Target as crate::traits::Problem>::Solution,
52    ) -> crate::rules::ExtractionResult<<Self::Source as crate::traits::Problem>::Solution> {
53        crate::rules::traits::validate_target_solution(self.target_problem(), target_solution)?;
54
55        let n2 = self.num_vertices_2;
56        (0..self.num_vertices_1)
57            .map(|residue| {
58                let mut selected =
59                    (0..n2).filter(|&mapped| target_solution[residue * n2 + mapped] == 1);
60                match (selected.next(), selected.next()) {
61                    (Some(mapped), None) => Ok(mapped + 1),
62                    (None, _) => Ok(0),
63                    (Some(_), Some(_)) => Err(crate::rules::ExtractionError::invalid(format!(
64                        "source residue {residue} maps to multiple target residues"
65                    ))),
66                }
67            })
68            .collect()
69    }
70}
71
72#[reduction(transform = {
73    exact {
74        num_vars = "num_vertices_1 * num_vertices_2 + num_contacts_1 * num_contacts_2",
75        num_constraints = "num_vertices_1 + num_vertices_2 + num_vertices_1 * (num_vertices_1 - 1) / 2 * num_vertices_2 * (num_vertices_2 + 1) / 2 + 2 * num_contacts_1 * num_contacts_2",
76        num_nonzeros = "2 * num_vertices_1 * num_vertices_2 + num_vertices_1 * (num_vertices_1 - 1) * num_vertices_2 * (num_vertices_2 + 1) / 2 + 4 * num_contacts_1 * num_contacts_2",
77    },
78    upper_bound {
79        max_constraint_magnitude_bits = "2",
80    },
81})]
82impl ReduceTo<ILP<bool>> for MaximumContactMapOverlap {
83    type Result = ReductionCMOToILP;
84
85    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
86        let n1 = self.num_vertices_1();
87        let n2 = self.num_vertices_2();
88        let contacts_1 = self.contacts_1();
89        let contacts_2 = self.contacts_2();
90
91        let num_x = n1 * n2;
92        let x_idx = |i: usize, j: usize| -> usize { i * n2 + j };
93
94        // y-variables: one per (contact in E_1) x (contact in E_2). The
95        // canonicalized contacts already satisfy i < k and j < l.
96        let num_y = contacts_1.len() * contacts_2.len();
97        let num_vars = num_x + num_y;
98        let y_idx = |seq: usize| -> usize { num_x + seq };
99
100        let mut constraints: Vec<LinearConstraint> = Vec::new();
101
102        // Row constraints: each residue of G_1 maps to at most one residue of G_2.
103        for i in 0..n1 {
104            let terms: Vec<(usize, i64)> = (0..n2).map(|j| (x_idx(i, j), 1)).collect();
105            constraints.push(LinearConstraint::le(terms, 1));
106        }
107
108        // Column constraints: each residue of G_2 receives at most one residue of G_1.
109        for j in 0..n2 {
110            let terms: Vec<(usize, i64)> = (0..n1).map(|i| (x_idx(i, j), 1)).collect();
111            constraints.push(LinearConstraint::le(terms, 1));
112        }
113
114        // Order-preservation: for i < k in V_1 and j >= l in V_2,
115        // forbid x_(i,j) + x_(k,l) <= 1. This rules out crossings (j > l)
116        // as well as equal-image matches (j == l).
117        for i in 0..n1 {
118            for k in (i + 1)..n1 {
119                for j in 0..n2 {
120                    for l in 0..=j {
121                        constraints.push(LinearConstraint::le(
122                            vec![(x_idx(i, j), 1), (x_idx(k, l), 1)],
123                            1,
124                        ));
125                    }
126                }
127            }
128        }
129
130        // Linking constraints for every contact pair: y_(i,k,j,l) <= x_(i,j) and
131        // y_(i,k,j,l) <= x_(k,l). Because each y has positive objective
132        // coefficient and there is no negative coupling, an optimum sets y to 1
133        // exactly when both endpoint-match variables are selected.
134        let mut seq = 0usize;
135        for &(i, k) in contacts_1 {
136            for &(j, l) in contacts_2 {
137                let yv = y_idx(seq);
138                constraints.push(LinearConstraint::le(vec![(yv, 1), (x_idx(i, j), -1)], 0));
139                constraints.push(LinearConstraint::le(vec![(yv, 1), (x_idx(k, l), -1)], 0));
140                seq += 1;
141            }
142        }
143
144        // Objective: maximize the number of preserved contacts.
145        let objective: Vec<(usize, i64)> = (0..num_y).map(|s| (y_idx(s), 1)).collect();
146
147        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Maximize)
148            .map_err(Self::target_construction)?;
149
150        Ok(ReductionCMOToILP {
151            target,
152            num_vertices_1: n1,
153            num_vertices_2: n2,
154        })
155    }
156}
157
158#[cfg(feature = "example-db")]
159pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
160    vec![crate::example_db::specs::RuleExampleSpec {
161        id: "maximumcontactmapoverlap_to_ilp",
162        build: || {
163            // Canonical instance from issue #1043:
164            //   G_1: n_1 = 4, E_1 = {{0,2}, {1,3}}
165            //   G_2: n_2 = 5, E_2 = {{0,3}, {1,4}, {0,2}}
166            // Optimal CMO alignment 0->0, 1->1, 2->3, 3->4 preserves 2 contacts.
167            let source = MaximumContactMapOverlap::new(
168                4,
169                vec![(0, 2), (1, 3)],
170                5,
171                vec![(0, 3), (1, 4), (0, 2)],
172            );
173            crate::example_db::specs::rule_example_via_ilp::<_, bool>(source)
174        },
175    }]
176}
177
178#[cfg(test)]
179#[path = "../unit_tests/rules/maximumcontactmapoverlap_ilp.rs"]
180mod tests;
problemreductions/rules/maximumcontactmapoverlap_ilp.rs

problemreductions/rules/
maximumcontactmapoverlap_ilp.rs
