diff --git a/docs/paper/reductions.typ b/docs/paper/reductions.typ index 8e0462b7a..3eb66f288 100644 --- a/docs/paper/reductions.typ +++ b/docs/paper/reductions.typ @@ -14022,6 +14022,103 @@ The following reductions to Integer Linear Programming are straightforward formu _Solution extraction._ For each source residue $i in V_1$, find the unique $j$ with $x_(i,j) = 1$ and output $j + 1$; if no such $j$ exists, output $0$ (the CMO sentinel for an unmatched residue). ] +#let cmo_clique = load-example( + "MaximumContactMapOverlap", + "MaximumClique", + target-variant: (graph: "SimpleGraph", weight: "One"), +) +#let cmo_clique_sol = cmo_clique.solutions.at(0) +#let cmo_clique_inst = cmo_clique.source.instance +#let cmo_clique_edges = cmo_clique.target.instance.graph.edges +#let cmo_clique_chosen = cmo_clique_sol.target_config.enumerate().filter(((i, x)) => x).map(((i, x)) => i) +#reduction-rule("MaximumContactMapOverlap", "MaximumClique", + example: true, + example-target-variant: (graph: "SimpleGraph", weight: "One"), + example-caption: [$|V_1| = #cmo_clique_inst.num_vertices_1$, $|E_1| = #cmo_clique_inst.contacts_1.len()$, $|V_2| = #cmo_clique_inst.num_vertices_2$, $|E_2| = #cmo_clique_inst.contacts_2.len()$], + extra: [ + #pred-commands( + "pred create --example " + rule-spec(cmo_clique) + " -o source.json", + "pred reduce source.json --via route.json -o bundle.json", + "pred solve bundle.json", + "pred evaluate source.json --config " + cli-config(cmo_clique_sol.source_config), + ) + #{ + let m1 = cmo_clique_inst.contacts_1.len() + let m2 = cmo_clique_inst.contacts_2.len() + let n = m1 * m2 + assert(cmo_clique.target.instance.graph.num_vertices == n) + // Recompute the assignments and the compatibility graph from the source instance. + let assign = () + for c1 in cmo_clique_inst.contacts_1 { + for c2 in cmo_clique_inst.contacts_2 { + assign.push(((c1.at(0), c2.at(0)), (c1.at(1), c2.at(1)))) + } + } + let sgn(x, y) = if x < y { -1 } else if x == y { 0 } else { 1 } + let adjacent(p, q) = p.all(a => q.all(b => sgn(a.at(0), b.at(0)) == sgn(a.at(1), b.at(1)))) + let computed = () + for p in range(n) { + for q in range(p + 1, n) { + if adjacent(assign.at(p), assign.at(q)) { computed.push((p, q)) } + } + } + assert(computed == cmo_clique_edges.map(e => (e.at(0), e.at(1)))) + // Cliques by exhaustive enumeration of vertex subsets. + let has-edge(p, q) = computed.contains((calc.min(p, q), calc.max(p, q))) + let sizes = () + for mask in range(calc.pow(2, n)) { + let sel = range(n).filter(i => calc.rem(calc.quo(mask, calc.pow(2, i)), 2) == 1) + if sel.all(p => sel.all(q => p == q or has-edge(p, q))) { sizes.push(sel.len()) } + } + let omega = calc.max(..sizes) + assert(omega == 2 and sizes.filter(k => k >= 3).len() == 0) + // The witness: a clique whose extraction preserves the listed contacts. + assert(cmo_clique_chosen.all(p => cmo_clique_chosen.all(q => p == q or has-edge(p, q)))) + let f = cmo_clique_sol.source_config + let preserved = cmo_clique_inst.contacts_1.filter(c => f.at(c.at(0)) > 0 and f.at(c.at(1)) > 0 and cmo_clique_inst.contacts_2.contains((f.at(c.at(0)) - 1, f.at(c.at(1)) - 1))) + assert(preserved.len() == omega and omega == cmo_clique_chosen.len()) + // Non-adjacent pair of the text: assignments reverse the order of source positions. + assert(not has-edge(0, 4)) + // Threshold three: all contacts preserved needs every position matched, i.e. the identity. + let covered = cmo_clique_inst.contacts_1.flatten().dedup().sorted() + assert(covered == range(cmo_clique_inst.num_vertices_1) and cmo_clique_inst.num_vertices_1 == cmo_clique_inst.num_vertices_2) + let n1 = cmo_clique_inst.num_vertices_1 + let identity = range(1, n1 + 1) + assert(f.len() == n1) + // Both contact maps are paths: n1 vertices, n1 - 1 contacts, degrees 1,1,2,...,2. + let degrees(inst-contacts, nv) = range(nv).map(v => inst-contacts.filter(c => c.contains(v)).len()).sorted() + let path-degrees = (1, 1) + (2,) * (n1 - 2) + assert(cmo_clique_inst.num_vertices_2 == n1 and m1 == n1 - 1 and m2 == n1 - 1) + assert(degrees(cmo_clique_inst.contacts_1, n1) == path-degrees and degrees(cmo_clique_inst.contacts_2, n1) == path-degrees) + // The assignments of vertices 0 and 4 that disagree on order. + let reversed = () + for a in assign.at(0) { + for b in assign.at(4) { + if sgn(a.at(0), b.at(0)) != sgn(a.at(1), b.at(1)) { reversed.push((a, b)) } + } + } + assert(reversed.len() == 1) + let (rev-a, rev-b) = reversed.at(0) + assert(rev-a.at(0) != rev-b.at(0) and rev-a.at(1) != rev-b.at(1)) + let identity-score = cmo_clique_inst.contacts_1.filter(c => cmo_clique_inst.contacts_2.contains((c.at(0), c.at(1)))).len() + assert(identity-score == 1 and m1 == 3) + [ + The $#m1 dot #m2 = #n$ vertices pair the contacts $E_1 = #cmo_clique_inst.contacts_1.map(c => $(#c.at(0), #c.at(1))$).join(", ")$ with $E_2 = #cmo_clique_inst.contacts_2.map(c => $(#c.at(0), #c.at(1))$).join(", ")$ in that order, vertex $i |E_2| + j$ pairing the $i$-th source contact with the $j$-th target contact. As unordered graphs both maps are paths on #n1 vertices, yet the compatibility graph has only #computed.len() edges, #computed.map(e => $(#e.at(0), #e.at(1))$).join(" and "): vertices $0$ and $4$ are not adjacent because their assignments #assign.at(0).map(a => $#a.at(0) -> #a.at(1)$).join(", ") and #assign.at(4).map(a => $#a.at(0) -> #a.at(1)$).join(", ") reverse the order of source positions $#calc.min(rev-a.at(0), rev-b.at(0)) < #calc.max(rev-a.at(0), rev-b.at(0))$ (#rev-a.at(0) $->$ #rev-a.at(1) and #rev-b.at(0) $->$ #rev-b.at(1)). The clique $\{#cmo_clique_chosen.map(str).join(", ")\}$ has size #cmo_clique_chosen.len() and extracts to the alignment $(#fmt-values(f))$, which preserves the contacts #preserved.map(c => $(#c.at(0), #c.at(1)) -> (#(f.at(c.at(0)) - 1), #(f.at(c.at(1)) - 1))$).join(" and "), so its score equals the clique size. The maximum clique size and the optimum are both #omega. *NO instance* (threshold $#(m1)$, all contacts preserved): every position $0, dots, #(cmo_clique_inst.num_vertices_1 - 1)$ is an endpoint of a source contact, so preserving all of them forces all positions to be matched, and the only increasing injection of #cmo_clique_inst.num_vertices_1 positions into #cmo_clique_inst.num_vertices_2 positions is the identity $(#fmt-values(identity))$, which preserves only #identity-score contact. Correspondingly the target has no clique of size $#m1$ (its #computed.len() edges form no triangle). + ] + } + ], +)[ + Build an order-preserving association graph: one unit-weight vertex for every pair of a source contact and a target contact, adjacent exactly when their proposed residue assignments are mutually consistent. A clique is then a partial injective order-preserving alignment, and the maximum clique size equals the maximum contact map overlap. The clique formulation of contact map overlap is due to @StricklandBarnesSokol2005CMO#footnote[The construction here was derived independently for the definitions of this library; the full text of the cited article was not compared.]. The target has $|E_1| |E_2|$ vertices and at most $|E_1|(|E_1| - 1)|E_2|(|E_2| - 1) slash 2$ edges. +][ + _Construction._ Let $G_1 = (V_1, E_1)$, $G_2 = (V_2, E_2)$ with contacts $(u, v)$, $u < v$, in stored order $e_0, dots, e_(m_1 - 1)$ and $d_0, dots, d_(m_2 - 1)$. For each $0 <= i < m_1$, $0 <= j < m_2$ with $e_i = (u, v)$ and $d_j = (a, b)$, create a vertex $p_(i j)$ with index $i m_2 + j$, weight $1$ and assignment set $A_(i j) = {(u, a), (v, b)}$. Distinct vertices $p, q$ are adjacent iff for all $(s, t) in A_p$ and $(s', t') in A_q$: $s < s' arrow.l.r.double t < t'$, $s = s' arrow.l.r.double t = t'$ and $s > s' arrow.l.r.double t > t'$. If either contact set is empty the graph is empty. + + _Correctness._ ($arrow.r.double$) Let $f$ be an alignment (increasing partial injection) preserving the contacts $S subset.eq E_1$. For $e_i in S$ let $j(i)$ index the target contact $(f(u), f(v))$. All assignments of the vertices $p_(i, j(i))$ lie on the graph of $f$, and a function that is injective and increasing on its matched positions preserves $<$, $=$ and $>$; hence these $|S|$ vertices form a clique. ($arrow.l.double$) Let $K$ be a clique and $R = union.big_(p in K) A_p$. Equality preservation makes $R$ a function, and strict-inequality preservation makes it injective and increasing, so $R$ is an alignment. Each $p_(i j) in K$ certifies that $e_i$ is preserved, and two clique vertices cannot share a source contact: they would assign the same images to its endpoints, hence share the target contact. So the score of $R$ is at least $|K|$. Therefore the maximum clique size equals the optimum, and the alignment extracted from a maximum clique is optimal. A clique that is not maximum can extract to an alignment scoring more than its size, because compatible assignments may imply further preserved contacts (if both maps are the triangle on positions $0, 1, 2$, the vertices pairing $(0, 1)$ with $(0, 1)$ and $(1, 2)$ with $(1, 2)$ already force the third contact). Hence no value map for non-optimal witnesses is claimed. + + _Solution extraction._ Reject a selection that is not a clique. Start from the all-unmatched vector of length $|V_1|$ and, for every selected vertex and each of its assignments $(s, t)$, write $t + 1$ at position $s$. + + _Overhead._ The vertex count $|E_1||E_2|$ is exact. Vertices sharing a source contact or a target contact are never adjacent, so at most $2 binom(|E_1|, 2) binom(|E_2|, 2)$ of the $binom(|E_1||E_2|, 2)$ vertex pairs are edges; the edge count is not exact in general. +] + #let mewkc_ilp = load-example( "MaximumEdgeWeightedKClique", "ILP", diff --git a/docs/paper/references.bib b/docs/paper/references.bib index 0dda81fd2..18f9f32ef 100644 --- a/docs/paper/references.bib +++ b/docs/paper/references.bib @@ -2352,3 +2352,14 @@ @article{een2006 pages = {1--26}, doi = {10.3233/SAT190014}, } + +@article{StricklandBarnesSokol2005CMO, + author = {Dawn M. Strickland and Earl Barnes and Joel S. Sokol}, + title = {Optimal Protein Structure Alignment Using Maximum Cliques}, + journal = {Operations Research}, + volume = {53}, + number = {3}, + pages = {389--402}, + year = {2005}, + doi = {10.1287/opre.1040.0189} +} diff --git a/src/rules/maximumcontactmapoverlap_ilp.rs b/src/rules/maximumcontactmapoverlap_ilp.rs index cbcd46917..d8d3e80f0 100644 --- a/src/rules/maximumcontactmapoverlap_ilp.rs +++ b/src/rules/maximumcontactmapoverlap_ilp.rs @@ -145,7 +145,7 @@ impl ReduceTo> for MaximumContactMapOverlap { let objective: Vec<(usize, i64)> = (0..num_y).map(|s| (y_idx(s), 1)).collect(); let target = ILP::new(num_vars, constraints, objective, ObjectiveSense::Maximize) - .map_err(Self::target_construction)?; + .map_err(>>::target_construction)?; Ok(ReductionCMOToILP { target, diff --git a/src/rules/maximumcontactmapoverlap_maximumclique.rs b/src/rules/maximumcontactmapoverlap_maximumclique.rs new file mode 100644 index 000000000..04993b927 --- /dev/null +++ b/src/rules/maximumcontactmapoverlap_maximumclique.rs @@ -0,0 +1,165 @@ +//! Reduction from MaximumContactMapOverlap to MaximumClique (unit weights) via an +//! order-preserving association graph. +//! +//! One target vertex is created for every pair (contact of `G_1`, contact of `G_2`); +//! it proposes the two residue assignments `u -> a` and `v -> b`. Two vertices are +//! adjacent when every pair of proposed assignments compares identically on both +//! sides (`<`, `=`, `>`), so a clique is a set of mutually consistent assignments: +//! a partial injective order-preserving alignment. A maximum clique has as many +//! vertices as the optimal alignment has preserved contacts. +//! +//! Background: Strickland, Barnes and Sokol, "Optimal Protein Structure Alignment +//! Using Maximum Cliques", Operations Research 53(3), 2005, formulate contact-map +//! overlap through maximum cliques. The construction here was derived for this +//! library's definitions and is not taken from that article. + +use crate::models::graph::{MaximumClique, MaximumContactMapOverlap}; +use crate::reduction; +use crate::rules::traits::{ReduceTo, ReductionResult}; +use crate::topology::SimpleGraph; +use crate::types::One; + +/// Result of reducing MaximumContactMapOverlap to MaximumClique. +/// +/// Target vertex `i * num_contacts_2 + j` stands for the pair (`i`-th contact of +/// `G_1`, `j`-th contact of `G_2`), in the order the contacts are stored. +#[derive(Debug, Clone)] +pub struct ReductionCMOToClique { + target: MaximumClique, + num_vertices_1: usize, + /// Per target vertex, its two assignments `(source residue, target residue)`. + assignments: Vec<[(usize, usize); 2]>, +} + +impl ReductionResult for ReductionCMOToClique { + type Source = MaximumContactMapOverlap; + type Target = MaximumClique; + + fn target_problem(&self) -> &Self::Target { + &self.target + } + + /// Write `t + 1` at source position `s` for every assignment `(s, t)` of a selected + /// vertex; every other position stays unmatched. + /// + /// The result of a maximum clique is an optimal alignment. A clique that is not + /// maximum may extract to an alignment preserving more contacts than the clique has + /// vertices, because compatible assignments can imply further preserved contacts. + fn extract_solution( + &self, + target_solution: &::Solution, + ) -> crate::rules::ExtractionResult<::Solution> { + crate::rules::traits::validate_target_witness( + self.target_problem(), + target_solution, + |value| value.0.is_some(), + "selected vertices are not a clique", + )?; + + let mut alignment = vec![0; self.num_vertices_1]; + for (pair, _) in self + .assignments + .iter() + .zip(target_solution) + .filter(|&(_, &selected)| selected) + { + for &(source, image) in pair { + alignment[source] = image + 1; + } + } + Ok(alignment) + } +} + +/// Two assignment pairs are compatible when all cross comparisons agree on both sides. +fn compatible(p: &[(usize, usize); 2], q: &[(usize, usize); 2]) -> bool { + p.iter() + .all(|&(s, t)| q.iter().all(|&(s2, t2)| s.cmp(&s2) == t.cmp(&t2))) +} + +/// Number of target vertices, an overflow error when it does not fit in `usize`. +fn pair_count( + num_contacts_1: usize, + num_contacts_2: usize, +) -> Result { + num_contacts_1.checked_mul(num_contacts_2).ok_or_else(|| { + crate::rules::ReductionError::integer_overflow::< + MaximumContactMapOverlap, + MaximumClique, + >("counting contact pairs") + }) +} + +#[reduction(transform = { + exact { + num_vertices = "num_contacts_1 * num_contacts_2", + }, + upper_bound { + num_edges = "num_contacts_1 * (num_contacts_1 - 1) * num_contacts_2 * (num_contacts_2 - 1) / 2", + }, +})] +impl ReduceTo> for MaximumContactMapOverlap { + type Result = ReductionCMOToClique; + + fn reduce_to(&self) -> Result { + let contacts_1 = self.contacts_1(); + let contacts_2 = self.contacts_2(); + let num_vertices = pair_count(contacts_1.len(), contacts_2.len())?; + + let mut assignments = Vec::new(); + for &(u, v) in contacts_1 { + for &(a, b) in contacts_2 { + assignments.push([(u, a), (v, b)]); + } + } + debug_assert_eq!(assignments.len(), num_vertices); + + let mut edges = Vec::new(); + for (p, left) in assignments.iter().enumerate() { + for (offset, right) in assignments[p + 1..].iter().enumerate() { + if compatible(left, right) { + edges.push((p, p + 1 + offset)); + } + } + } + + let graph = SimpleGraph::try_new(num_vertices, edges) + .map_err(>>::target_construction)?; + Ok(ReductionCMOToClique { + target: MaximumClique::new(graph, vec![One; num_vertices]), + num_vertices_1: self.num_vertices_1(), + assignments, + }) + } +} + +#[cfg(feature = "example-db")] +pub(crate) fn canonical_rule_example_specs() -> Vec { + use crate::export::SolutionPair; + + vec![crate::example_db::specs::RuleExampleSpec { + id: "maximumcontactmapoverlap_to_maximumclique", + build: || { + let source = MaximumContactMapOverlap::new( + 4, + vec![(0, 2), (1, 3), (0, 3)], + 4, + vec![(0, 1), (2, 3), (0, 3)], + ); + let mut clique = vec![false; 9]; + clique[0] = true; + clique[8] = true; + crate::example_db::specs::rule_example_with_witness::<_, MaximumClique>( + source, + SolutionPair { + source_config: serde_json::json!(vec![1, 0, 2, 4]), + target_config: serde_json::json!(clique), + }, + ) + }, + }] +} + +#[cfg(test)] +#[path = "../unit_tests/rules/maximumcontactmapoverlap_maximumclique.rs"] +mod tests; diff --git a/src/rules/mod.rs b/src/rules/mod.rs index 7606f4923..0c8bdadec 100644 --- a/src/rules/mod.rs +++ b/src/rules/mod.rs @@ -240,6 +240,7 @@ pub(crate) mod maximumclique_ilp; pub(crate) mod maximumcokplex_ilp; pub(crate) mod maximumcommonedgesubgraph_ilp; pub(crate) mod maximumcontactmapoverlap_ilp; +pub(crate) mod maximumcontactmapoverlap_maximumclique; pub(crate) mod maximumdomaticnumber_ilp; pub(crate) mod maximumedgeweightedkclique_ilp; pub(crate) mod maximumleafspanningtree_ilp; @@ -431,6 +432,7 @@ pub(crate) fn canonical_rule_example_specs() -> Vec MaximumContactMapOverlap { + MaximumContactMapOverlap::new( + 4, + vec![(0, 2), (1, 3), (0, 3)], + 4, + vec![(0, 1), (2, 3), (0, 3)], + ) +} + +fn reduce(source: &MaximumContactMapOverlap) -> ReductionCMOToClique { + ReduceTo::>::reduce_to(source).expect("reduction succeeds") +} + +fn edge_set(reduction: &ReductionCMOToClique) -> Vec<(usize, usize)> { + let mut edges: Vec<_> = reduction + .target_problem() + .graph() + .edges() + .into_iter() + .map(|(u, v)| (u.min(v), u.max(v))) + .collect(); + edges.sort_unstable(); + edges +} + +fn selection(len: usize, chosen: &[usize]) -> Vec { + (0..len).map(|i| chosen.contains(&i)).collect() +} + +fn optimum(source: &MaximumContactMapOverlap) -> i64 { + let witness = BruteForce::new().solve(source).unwrap().unwrap(); + source.evaluate(&witness).unwrap().0.unwrap() +} + +/// Checks the optimum equality, extraction from every clique and from every maximum clique. +fn check_instance(source: &MaximumContactMapOverlap) { + let reduction = reduce(source); + let target = reduction.target_problem(); + let n = target.num_vertices(); + let opt = optimum(source); + let mut best = 0; + let mut cliques = Vec::new(); + for mask in 0u32..(1 << n) { + let config: Vec = (0..n).map(|i| mask >> i & 1 == 1).collect(); + if let Max(Some(size)) = target.evaluate(&config).unwrap() { + best = best.max(size); + cliques.push((config, size)); + } + } + assert_eq!(best, opt, "{source:?}"); + for (config, size) in cliques { + let alignment = reduction.extract_solution(&config).unwrap(); + let score = source + .evaluate(&alignment) + .unwrap() + .0 + .expect("feasible alignment"); + assert!(score >= size, "{source:?}: score {score} < clique {size}"); + if size == best { + assert_eq!(score, opt, "{source:?}"); + } + } +} + +fn contact_sets(n: usize, max_len: usize) -> Vec> { + let pairs: Vec<_> = (0..n) + .flat_map(|u| ((u + 1)..n).map(move |v| (u, v))) + .collect(); + (0u32..(1 << pairs.len())) + .map(|mask| { + pairs + .iter() + .enumerate() + .filter(|(i, _)| mask >> i & 1 == 1) + .map(|(_, &pair)| pair) + .collect::>() + }) + .filter(|set| set.len() <= max_len) + .collect() +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_closed_loop() { + let source = example(); + let reduction = reduce(&source); + assert_eq!(reduction.target_problem().num_vertices(), 9); + assert_eq!(edge_set(&reduction), vec![(0, 8), (4, 8)]); + assert_optimization_round_trip_from_optimization_target( + &source, + &reduction, + "MaximumContactMapOverlap->MaximumClique closed loop", + ); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_canonical_example() { + let source = example(); + let reduction = reduce(&source); + let len = 9; + assert_eq!( + reduction + .extract_solution(&selection(len, &[0, 8])) + .unwrap(), + vec![1, 0, 2, 4] + ); + assert_eq!( + reduction.extract_solution(&selection(len, &[0])).unwrap(), + vec![1, 0, 2, 0] + ); + assert_eq!( + reduction.extract_solution(&selection(len, &[4])).unwrap(), + vec![0, 3, 0, 4] + ); + assert_eq!( + reduction.extract_solution(&selection(len, &[])).unwrap(), + vec![0; 4] + ); + // Candidates 0 and 4 assign crossing images to source positions 1 and 2. + assert!(!edge_set(&reduction).contains(&(0, 4))); + // The optimum is two; the unordered contact graphs are isomorphic paths but the + // order forbids preserving all three contacts. + assert_eq!(optimum(&source), 2); + let max_clique = BruteForce::new() + .solve(reduction.target_problem()) + .unwrap() + .unwrap(); + assert_eq!( + reduction.target_problem().evaluate(&max_clique).unwrap(), + Max(Some(2)) + ); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_all_small_pairs() { + // Every pair of contact sets on up to three residues per side. + for n1 in 1..=3 { + for n2 in 1..=3 { + for e1 in contact_sets(n1, usize::MAX) { + for e2 in contact_sets(n2, usize::MAX) { + check_instance(&MaximumContactMapOverlap::new(n1, e1.clone(), n2, e2)); + } + } + } + } +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_sparse_pairs_to_five_residues() { + // Every pair of contact sets with at most two contacts per side, up to five residues + // (source brute force is limited to alignment spaces of at most 1300 configurations). + for n1 in 1..=5usize { + for n2 in 1..=5usize { + if (n2 + 1).pow(n1 as u32) > 1300 || (n1 <= 3 && n2 <= 3) { + continue; + } + for e1 in contact_sets(n1, 2) { + for e2 in contact_sets(n2, 2) { + check_instance(&MaximumContactMapOverlap::new(n1, e1.clone(), n2, e2)); + } + } + } + } +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_empty_contacts() { + for source in [ + MaximumContactMapOverlap::new(3, vec![], 3, vec![(0, 1)]), + MaximumContactMapOverlap::new(3, vec![(0, 1)], 3, vec![]), + MaximumContactMapOverlap::new(0, vec![], 0, vec![]), + ] { + let reduction = reduce(&source); + assert_eq!(reduction.target_problem().num_vertices(), 0); + assert_eq!(reduction.target_problem().num_edges(), 0); + assert_eq!( + reduction.extract_solution(&vec![]).unwrap(), + vec![0; source.num_vertices_1()] + ); + } +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_more_source_than_target_vertices() { + let source = MaximumContactMapOverlap::new(5, vec![(0, 4), (1, 2), (3, 4)], 2, vec![(0, 1)]); + check_instance(&source); + assert_eq!(optimum(&source), 1); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_compatibility_rules() { + // Shared source endpoint: contacts (0,1),(0,2) onto (0,1),(0,2),(1,2). + let source = + MaximumContactMapOverlap::new(3, vec![(0, 1), (0, 2)], 3, vec![(0, 1), (0, 2), (1, 2)]); + let reduction = reduce(&source); + let edges = edge_set(&reduction); + let vertex = |i: usize, j: usize| i * 3 + j; + // Images of position 0 agree: 0->0 in both. + assert!(edges.contains(&(vertex(0, 0), vertex(1, 1)))); + // Images of position 0 disagree: 0->0 versus 0->1. + assert!(!edges.contains(&(vertex(0, 0), vertex(1, 2)))); + // Vertices of one source contact are never adjacent, nor those of one target contact. + assert!(!edges.contains(&(vertex(0, 0), vertex(0, 1)))); + assert!(!edges.contains(&(vertex(0, 0), vertex(1, 0)))); + + // Disjoint source contacts cannot reuse a target endpoint. + let source = MaximumContactMapOverlap::new(4, vec![(0, 1), (2, 3)], 3, vec![(0, 1), (1, 2)]); + assert!(!edge_set(&reduce(&source)).contains(&(0, 3))); + + // Crossing assignments are excluded. + let source = MaximumContactMapOverlap::new(4, vec![(0, 2), (1, 3)], 4, vec![(0, 3), (1, 2)]); + assert!(!edge_set(&reduce(&source)).contains(&(0, 3))); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_vertices_follow_stored_contact_order() { + // Contacts keep the order they were given in; vertex `i * m2 + j` pairs the `i`-th + // source contact with the `j`-th target contact. + let source = MaximumContactMapOverlap::new(4, vec![(2, 3), (0, 1)], 4, vec![(2, 3), (0, 1)]); + let reduction = reduce(&source); + let len = 4; + assert_eq!( + reduction.extract_solution(&selection(len, &[0])).unwrap(), + vec![0, 0, 3, 4] + ); + assert_eq!( + reduction.extract_solution(&selection(len, &[1])).unwrap(), + vec![0, 0, 1, 2] + ); + assert_eq!( + reduction.extract_solution(&selection(len, &[2])).unwrap(), + vec![3, 4, 0, 0] + ); + // (2,3)->(2,3) with (0,1)->(0,1) is the identity alignment. + assert_eq!( + reduction + .extract_solution(&selection(len, &[0, 3])) + .unwrap(), + vec![1, 2, 3, 4] + ); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_threshold_three_is_infeasible() { + let source = example(); + // Threshold three is infeasible on both sides. + assert!(optimum(&source) < 3); + let reduction = reduce(&source); + assert!( + BruteForce::new() + .find_all_witnesses(reduction.target_problem()) + .unwrap() + .iter() + .all(|w| w.iter().filter(|&&b| b).count() < 3) + ); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_extraction_rejects_malformed_witnesses() { + let reduction = reduce(&example()); + // Vertices 0 and 4 are not adjacent. + assert!(reduction.extract_solution(&selection(9, &[0, 4])).is_err()); + // All vertices selected. + assert!(reduction.extract_solution(&vec![true; 9]).is_err()); + // Wrong lengths. + assert!(reduction.extract_solution(&vec![false; 8]).is_err()); + assert!(reduction.extract_solution(&vec![false; 10]).is_err()); + assert!(reduction.extract_solution(&vec![]).is_err()); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_nonmaximum_clique_scores_above_its_size() { + // Both maps are the triangle on positions 0, 1, 2. The vertices pairing contact + // (0,1) with (0,1) and (1,2) with (1,2) form a clique of size two whose alignment + // also preserves the contact (0,2). + let triangle = vec![(0, 1), (0, 2), (1, 2)]; + let source = MaximumContactMapOverlap::new(3, triangle.clone(), 3, triangle); + let reduction = reduce(&source); + let witness = selection(9, &[0, 2 * 3 + 2]); + // contacts are stored as (0,1),(0,2),(1,2): index 2 is (1,2) + let alignment = reduction.extract_solution(&witness).unwrap(); + assert_eq!(alignment, vec![1, 2, 3]); + let size = reduction.target_problem().evaluate(&witness).unwrap(); + assert_eq!(size, Max(Some(2))); + let score = source.evaluate(&alignment).unwrap(); + assert_eq!(score, Max(Some(3))); + assert!(score.0 > size.0); + // At the optimum the two sides agree. + assert_eq!(optimum(&source), 3); +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_parameter_predictions() { + for source in [ + example(), + MaximumContactMapOverlap::new(3, vec![(0, 1)], 3, vec![]), + MaximumContactMapOverlap::new(4, vec![(0, 2), (1, 3)], 5, vec![(0, 3), (1, 4), (0, 2)]), + ] { + let reduction = reduce(&source); + crate::rules::test_helpers::assert_parameter_predictions(&source, &reduction); + let m = source.num_contacts_1() * source.num_contacts_2(); + assert_eq!(reduction.target_problem().num_vertices(), m); + } +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_pair_count_boundaries() { + assert_eq!(pair_count(0, usize::MAX).unwrap(), 0); + assert_eq!(pair_count(usize::MAX, 1).unwrap(), usize::MAX); + assert_eq!(pair_count(1 << 32, 1 << 31).unwrap(), 1 << 63); + for (a, b) in [(usize::MAX, 2), (1 << 32, 1 << 32)] { + assert!(matches!( + pair_count(a, b), + Err(crate::rules::ReductionError::IntegerOverflow { .. }) + )); + } +} + +#[test] +fn test_maximumcontactmapoverlap_to_maximumclique_edge_bound_holds() { + let mut seen = HashSet::new(); + for e1 in contact_sets(4, 4) { + for e2 in contact_sets(4, 4) { + let source = MaximumContactMapOverlap::new(4, e1.clone(), 4, e2.clone()); + let reduction = reduce(&source); + let (m1, m2) = (e1.len(), e2.len()); + let bound = m1 * m1.saturating_sub(1) * m2 * m2.saturating_sub(1) / 2; + assert!(reduction.target_problem().num_edges() <= bound); + seen.insert((m1, m2)); + } + } + assert!(seen.len() > 10); +}