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(
73    transform = 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    },
77    unavailable = {
78        num_nonzeros = "the exact target parameter is not represented by this reduction's symbolic transform",
79    }
80)]
81impl ReduceTo<ILP<bool>> for MaximumContactMapOverlap {
82    type Result = ReductionCMOToILP;
83
84    fn reduce_to(&self) -> Result<Self::Result, crate::rules::ReductionError> {
85        let n1 = self.num_vertices_1();
86        let n2 = self.num_vertices_2();
87        let contacts_1 = self.contacts_1();
88        let contacts_2 = self.contacts_2();
89
90        let num_x = n1 * n2;
91        let x_idx = |i: usize, j: usize| -> usize { i * n2 + j };
92
93        // y-variables: one per (contact in E_1) x (contact in E_2). The
94        // canonicalized contacts already satisfy i < k and j < l.
95        let num_y = contacts_1.len() * contacts_2.len();
96        let num_vars = num_x + num_y;
97        let y_idx = |seq: usize| -> usize { num_x + seq };
98
99        let mut constraints: Vec<LinearConstraint> = Vec::new();
100
101        // Row constraints: each residue of G_1 maps to at most one residue of G_2.
102        for i in 0..n1 {
103            let terms: Vec<(usize, i64)> = (0..n2).map(|j| (x_idx(i, j), 1)).collect();
104            constraints.push(LinearConstraint::le(terms, 1));
105        }
106
107        // Column constraints: each residue of G_2 receives at most one residue of G_1.
108        for j in 0..n2 {
109            let terms: Vec<(usize, i64)> = (0..n1).map(|i| (x_idx(i, j), 1)).collect();
110            constraints.push(LinearConstraint::le(terms, 1));
111        }
112
113        // Order-preservation: for i < k in V_1 and j >= l in V_2,
114        // forbid x_(i,j) + x_(k,l) <= 1. This rules out crossings (j > l)
115        // as well as equal-image matches (j == l).
116        for i in 0..n1 {
117            for k in (i + 1)..n1 {
118                for j in 0..n2 {
119                    for l in 0..=j {
120                        constraints.push(LinearConstraint::le(
121                            vec![(x_idx(i, j), 1), (x_idx(k, l), 1)],
122                            1,
123                        ));
124                    }
125                }
126            }
127        }
128
129        // Linking constraints for every contact pair: y_(i,k,j,l) <= x_(i,j) and
130        // y_(i,k,j,l) <= x_(k,l). Because each y has positive objective
131        // coefficient and there is no negative coupling, an optimum sets y to 1
132        // exactly when both endpoint-match variables are selected.
133        let mut seq = 0usize;
134        for &(i, k) in contacts_1 {
135            for &(j, l) in contacts_2 {
136                let yv = y_idx(seq);
137                constraints.push(LinearConstraint::le(vec![(yv, 1), (x_idx(i, j), -1)], 0));
138                constraints.push(LinearConstraint::le(vec![(yv, 1), (x_idx(k, l), -1)], 0));
139                seq += 1;
140            }
141        }
142
143        // Objective: maximize the number of preserved contacts.
144        let objective: Vec<(usize, i64)> = (0..num_y).map(|s| (y_idx(s), 1)).collect();
145
146        let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Maximize)
147            .map_err(Self::target_construction)?;
148
149        Ok(ReductionCMOToILP {
150            target,
151            num_vertices_1: n1,
152            num_vertices_2: n2,
153        })
154    }
155}
156
157#[cfg(feature = "example-db")]
158pub(crate) fn canonical_rule_example_specs() -> Vec<crate::example_db::specs::RuleExampleSpec> {
159    vec![crate::example_db::specs::RuleExampleSpec {
160        id: "maximumcontactmapoverlap_to_ilp",
161        build: || {
162            // Canonical instance from issue #1043:
163            //   G_1: n_1 = 4, E_1 = {{0,2}, {1,3}}
164            //   G_2: n_2 = 5, E_2 = {{0,3}, {1,4}, {0,2}}
165            // Optimal CMO alignment 0->0, 1->1, 2->3, 3->4 preserves 2 contacts.
166            let source = MaximumContactMapOverlap::new(
167                4,
168                vec![(0, 2), (1, 3)],
169                5,
170                vec![(0, 3), (1, 4), (0, 2)],
171            );
172            crate::example_db::specs::rule_example_via_ilp::<_, bool>(source)
173        },
174    }]
175}
176
177#[cfg(test)]
178#[path = "../unit_tests/rules/maximumcontactmapoverlap_ilp.rs"]
179mod tests;
problemreductions/rules/maximumcontactmapoverlap_ilp.rs

problemreductions/rules/
maximumcontactmapoverlap_ilp.rs
