Repository navigation
Conversation
Co-Authored-By: Claude Sonnet 5.5 <noreply@anthropic.com>
Collaborator
Author
Agentic Review ReportSubject: rule Structural and semantic
Quality
Feature testI used the
FindingsNone open. All findings from earlier rounds are resolved:
|
Collaborator
Author
Verification certificateVerdict: verified. The construction is correct and the proof holds; no counterexample or disagreement was found. Type chain. Source Construction. One vertex per pair of contacts (i, j), numbered Checks (all 0 failures).
Scripts follow (each collapsed). Proof (Typst)#set document(title: "MaximumContactMapOverlap to MaximumClique")
#set page(margin: 2cm)
#set heading(numbering: "1.")
= MaximumContactMapOverlap $arrow.r$ MaximumClique⟨SimpleGraph, One⟩
*Provenance.* The reduction of contact-map overlap to maximum clique on a graph of compatible residue pairings is due to Strickland, Barnes and Sokol (Operations Research 53(3), 389--402, 2005; DOI 10.1287/opre.1040.0189). The full text could not be obtained (only the abstract and secondary summaries were reachable), so the construction and proof below are *derived here*, cf. Strickland et al.; the cited paper is background only. Rahrig, Leontis and Zirbel (Bioinformatics 26(21), 2010) describe the same order-preserving alignment graph for the two-contact case.
*Type resolution.* Source: `MaximumContactMapOverlap` has `type Value = Max<i64>` (`src/models/graph/maximum_contact_map_overlap.rs`), configuration `Vec<usize>` of length $n_1$, `0` unmatched, `j+1` image $j$, nonzero entries strictly increasing; the all-unmatched vector is always feasible. Target: `MaximumClique<SimpleGraph, One>` has `type Value = Max<W::Sum>` with `<One as WeightElement>::Sum = i64` (`src/models/graph/maximum_clique.rs`, `src/types.rs`), i.e. `Max<i64>`, configuration `Vec<bool>`, value = clique cardinality. `Max<i64>` $->$ `Max<i64>`: witness-compatible. Redundancy: `pred path MaximumContactMapOverlap MaximumClique/SimpleGraph/One` reports no path; the only outgoing edge of the source is to `ILP/bool/i64/general`, and no path leads back from ILP to MaximumClique. The new edge is not dominated. This is an optimization rule.
*Theorem.* Let $I = (n_1, E_1, n_2, E_2)$ be an instance with contact lists $E_1 = [e_1, dots, e_(m_1)]$, $E_2 = [d_1, dots, d_(m_2)]$ (each contact $(u,v)$ with $u<v$). The graph $G_I$ below has $m_1 m_2$ vertices, its maximum clique size equals the optimum of $I$, and every maximum clique decodes by the extraction below to an optimal alignment.
== Construction
+ *Vertices.* For each $i in [m_1]$, $j in [m_2]$, with $e_i = (u,v)$ and $d_j = (a,b)$, create one vertex $p_(i j)$ with assignment set $A_(i j) = {(u,a), (v,b)}$. Order vertices by $i$ first, then $j$; index $p_(i j) = i m_2 + j$. Every vertex has weight 1.
+ *Edges.* Distinct vertices $p, q$ are adjacent iff for every $(s,t) in A_p$ and $(s', t') in A_q$: $"sgn"(s - s') = "sgn"(t - t')$, i.e. $s<s' arrow.l.r.double t<t'$, $s=s' arrow.l.r.double t=t'$, $s>s' arrow.l.r.double t>t'$.
+ *Extraction.* Given a vertex set $K$ that is a clique, start with $f = (0, dots, 0)$ of length $n_1$; for every $p in K$ and every $(s,t) in A_p$ set $f[s] = t+1$. Return $f$. If $K$ is not a clique, extraction reports an error.
Notation: for a feasible alignment $f$ the _score_ $sigma(f)$ is the number of $(u,v) in E_1$ with $f(u), f(v)$ both matched and $(f(u), f(v)) in E_2$ (as $f$ is increasing, $f(u)<f(v)$).
== Overhead
#table(columns: 3,
[target parameter], [formula], [contract],
[`num_vertices`], [$m_1 m_2$ (`num_contacts_1 * num_contacts_2`)], [exact],
[`num_edges`], [$<= binom(m_1 m_2, 2)$; sharper: $<= 2 binom(m_1,2) binom(m_2,2)$], [upper bound],
)
The vertex count is exact by construction. For the edge bound: two vertices $p_(i j), p_(i j')$ with $j != j'$ share the source contact $e_i = (u,v)$ so they both assign $u$, to $a$ and to $a'$ respectively; adjacency needs $a = a'$ and (from $v$) $b = b'$, giving $d_j = d_(j')$, contradiction. Symmetrically, let $i != i'$ with $e_i = (u,v)$, $e_(i') = (u',v')$ and $d_j = (a,b)$. If $u != u'$, the assignments $(u,a)$ and $(u',a)$ have $u != u'$ but equal images, violating equality preservation; if $u = u'$ then $v != v'$ and $(v,b),(v',b)$ violate it in the same way. So only pairs differing in both indices can be edges, and their number is $binom(m_1 m_2, 2) - m_1 binom(m_2, 2) - m_2 binom(m_1, 2) = 2 binom(m_1,2) binom(m_2,2)$. The bound $binom(N, 2)$, $N = m_1 m_2$, is attained only when $N <= 1$, so it is not exact for general inputs and is declared as an upper bound. If $m_1 = 0$ or $m_2 = 0$ the graph is empty. Construction time and output size are $O(N + N^2)$.
== Correctness
*Lemma 1 (a vertex is a locally increasing injection).* $A_(i j) = {(u,a),(v,b)}$ has $u<v$ and $a<b$.
*(⇒) Source alignment gives a clique.* Let $f$ be a feasible alignment with score $sigma$. Let $S = {(u,v) in E_1 : f(u), f(v) "matched", (f(u),f(v)) in E_2}$, $|S| = sigma$. For each $e_i in S$ let $j(i)$ be the index of the target contact $(f(u), f(v))$; it is unique since $E_2$ has no duplicates. The vertices $p_(i, j(i))$ for $e_i in S$ are distinct (their first index differs), and $A_(i, j(i)) subset.eq {(s, f(s))}$, the graph of $f$. Take two of them and pairs $(s,t),(s',t')$ from their assignment sets: all four are of the form $(s, f(s))$. As $f$ is a function, $s = s' arrow.double t = t'$; as $f$ is injective, $t = t' arrow.double s = s'$; as $f$ is strictly increasing on its matched positions, $s<s' arrow.l.r.double t<t'$. All three comparisons agree, so the vertices are adjacent. Thus $S$ yields a clique of size $sigma$, and $omega(G_I) >= "OPT"(I)$.
*(⇐) Clique gives an alignment.* Let $K$ be a clique and $R = union.big_(p in K) A_p$. (i) $R$ is a function: if $(s,t),(s,t') in R$, either they come from the same vertex (then $t=t'$ because $s$ occurs once in an assignment set, by $u<v$) or from distinct adjacent vertices, where $s=s$ forces $t=t'$. (ii) $R$ is injective and strictly increasing: if $(s,t),(s',t') in R$ with $s<s'$, adjacency (or Lemma 1 within one vertex) gives $t<t'$. Hence the extracted $f$, with $f[s]=t+1$ exactly for $(s,t) in R$, is a feasible alignment (nonzero entries strictly increasing; $t<n_2$ so $t+1 <= n_2$). Each $p_(i j) in K$ certifies the source contact $e_i$ as preserved: $f(u)=a$, $f(v)=b$, and $(a,b)=d_j in E_2$. Distinct clique vertices certify distinct source contacts: if $p_(i j), p_(i j') in K$ then, as shown above, $d_j = d_(j')$ so $j=j'$. Therefore $sigma(f) >= |K|$ and $"OPT"(I) >= omega(G_I)$.
Combining both directions: $omega(G_I) = "OPT"(I)$. Moreover, for a maximum clique $K$, $"OPT" >= sigma(f) >= |K| = omega = "OPT"$, so $f$ is optimal. For an arbitrary (non-maximum) clique only the inequality $sigma(f) >= |K|$ holds: pairwise compatibility can force extra preserved contacts that were not selected (e.g. source and target both the triangle on positions $0,1,2$; the vertices $p_(0 0)$ and $p_(1 1)$ extract $f = [1,2,3]$ with score 3 but the clique has size 2). Hence the result mapping is the identity only at optimum, and no value map is claimed for non-optimal witnesses. Also $omega <= min(m_1, m_2)$ since distinct clique vertices have distinct source contacts and (by the same argument with roles swapped) distinct target contacts.
== Solution extraction
Described as step 3 of the construction. It validates the clique once (every pair of selected vertices is an edge, selection length $= m_1 m_2$), then writes $t+1$ at source position $s$ for each assignment. Special cases: the empty clique yields the all-unmatched vector; an empty contact set yields an empty graph and the all-unmatched vector of length $n_1$; $n_1 > n_2$ needs no special case because alignments are partial; contacts sharing a source endpoint must agree on its image (equality preservation); disjoint source contacts cannot share a target endpoint; crossing assignments are non-adjacent (strict-inequality preservation); no assumption about connectedness or isolated positions is used. A selection that is not a clique, or of the wrong length, is rejected.
== YES example
$n_1 = n_2 = 4$, $E_1 = [(0,2),(1,3),(0,3)]$, $E_2 = [(0,1),(2,3),(0,3)]$ (both are paths on four vertices, but the orders differ). Vertices $p_(i j)$, index $3i+j$:
#table(columns: 3,
[vertex], [pair $(i,j)$], [assignments],
[0], [(0,0)], [$0 -> 0, 2 -> 1$],
[1], [(0,1)], [$0 -> 2, 2 -> 3$],
[2], [(0,2)], [$0 -> 0, 2 -> 3$],
[3], [(1,0)], [$1 -> 0, 3 -> 1$],
[4], [(1,1)], [$1 -> 2, 3 -> 3$],
[5], [(1,2)], [$1 -> 0, 3 -> 3$],
[6], [(2,0)], [$0 -> 0, 3 -> 1$],
[7], [(2,1)], [$0 -> 2, 3 -> 3$],
[8], [(2,2)], [$0 -> 0, 3 -> 3$],
)
Overhead: 9 vertices ($= 3 dot 3$), 2 edges $<= binom(9,2) = 36$ and $<= 2 binom(3,2) binom(3,2) = 18$. The edge list is $[(0,8),(4,8)]$: e.g. $p_0, p_8$ assign ${0->0, 2->1}$ and ${0->0, 3->3}$; the comparisons $0=0$ with $0=0$, $0<3$ with $0<3$, $2>0$ with $1>0$, $2<3$ with $1<3$ all agree. Vertices $0$ and $4$ are not adjacent: $2 -> 1$ and $1 -> 2$ reverse the order of source positions $1<2$. The clique ${0,8}$ has size 2 and extracts $f = [1,0,2,4]$: position 0 to 0, position 1 unmatched, position 2 to 1, position 3 to 3; it preserves $(0,2) -> (0,1)$ and $(0,3) -> (0,3)$, score 2 $=$ optimum. The suboptimal cliques ${0}$ and ${4}$ extract $[1,0,2,0]$ and $[0,3,0,4]$, each with score 1; the empty clique extracts $[0,0,0,0]$ with score 0.
== NO example
Same instance with threshold 3 (all three source contacts preserved). Target side: the maximum clique has size 2 since the only edges are $(0,8)$ and $(4,8)$, which form no triangle, so no clique of size 3 exists. Source side: preserving three contacts needs positions $0,1,2,3$ all matched (contacts $(0,2),(1,3),(0,3)$ cover every position); an increasing injection of four positions into four positions is the identity $[1,2,3,4]$, which preserves only $(0,3)$ (contacts $(0,2)$ and $(1,3)$ map to $(0,2)$ and $(1,3)$, absent from $E_2$). So the optimum is 2 and threshold 3 is infeasible on both sides, although the two contact graphs are isomorphic as unordered graphs (an order-ignoring construction would report 3).
Construction script (Python)#!/usr/bin/env python3
"""Constructor verification: MaximumContactMapOverlap -> MaximumClique<SimpleGraph, One>.
Source: Max<i64>, config = Vec<usize> of length n1, 0 = unmatched, j+1 = image j,
nonzero entries strictly increasing (all-unmatched always feasible).
Target: Max<i64> (One weights), config = Vec<bool>; a clique has value = its size.
Seven sections (see the audit at the end); the last line is
TOTAL: <checks> checks, <failures> failures
"""
import itertools
import json
import os
import random
import sys
import sympy as sp
from hypothesis import HealthCheck, given, settings, strategies as st
SEED = 20261103
random.seed(SEED)
COUNTS = {}
FAILS = {}
FAIL_MSGS = []
def check(section, cond, msg=""):
COUNTS[section] = COUNTS.get(section, 0) + 1
if not cond:
FAILS[section] = FAILS.get(section, 0) + 1
if len(FAIL_MSGS) < 20:
FAIL_MSGS.append(f"[{section}] {msg}")
# ---------------------------------------------------------------- reduction
def canon(contacts):
return sorted((min(u, v), max(u, v)) for u, v in contacts)
def cmp(x, y):
return (x > y) - (x < y)
def reduce(n1, E1, n2, E2):
"""Return (candidates, adjacency) where candidates[p] = (i, j, assignments)."""
E1, E2 = list(E1), list(E2)
cands = []
for i, (u, v) in enumerate(E1):
for j, (a, b) in enumerate(E2):
cands.append((i, j, ((u, a), (v, b))))
N = len(cands)
edges = set()
for p in range(N):
for q in range(p + 1, N):
ok = all(
cmp(s, s2) == cmp(t, t2)
for (s, t) in cands[p][2]
for (s2, t2) in cands[q][2]
)
if ok:
edges.add((p, q))
return cands, edges
def extract(n1, n2, cands, edges, sel):
"""sel: set of selected vertices; returns None if not a clique."""
sel = sorted(sel)
for x, y in itertools.combinations(sel, 2):
if (x, y) not in edges:
return None
cfg = [0] * n1
for p in sel:
for (s, t) in cands[p][2]:
cfg[s] = t + 1
return cfg
def extract_bits(n1, n2, cands, edges, bits):
"""Target-shaped witness (Vec<bool> of length N); None if wrong length or not a clique."""
if len(bits) != len(cands):
return None
return extract(n1, n2, cands, edges, {p for p, b in enumerate(bits) if b})
# ------------------------------------------------------- source semantics
def feasible(cfg, n2):
nz = [c for c in cfg if c != 0]
return all(0 <= c <= n2 for c in cfg) and all(a < b for a, b in zip(nz, nz[1:]))
def score(cfg, E1, E2):
S2 = set(E2)
cnt = 0
for (u, v) in E1:
if cfg[u] and cfg[v]:
a, b = cfg[u] - 1, cfg[v] - 1
if (min(a, b), max(a, b)) in S2:
cnt += 1
return cnt
def all_alignments(n1, n2):
"""Every feasible config, by choosing the matched subset and an increasing image."""
out = []
for k in range(0, min(n1, n2) + 1):
for src in itertools.combinations(range(n1), k):
for img in itertools.combinations(range(n2), k):
cfg = [0] * n1
for s, t in zip(src, img):
cfg[s] = t + 1
out.append(cfg)
return out
_ALIGN_CACHE = {}
def alignments(n1, n2):
if (n1, n2) not in _ALIGN_CACHE:
_ALIGN_CACHE[(n1, n2)] = all_alignments(n1, n2)
return _ALIGN_CACHE[(n1, n2)]
def max_clique_size(N, edges):
adj = [set() for _ in range(N)]
for x, y in edges:
adj[x].add(y)
adj[y].add(x)
best = 0
def bk(R, P, X):
nonlocal best
if not P and not X:
best = max(best, R)
return
if R + len(P) <= best:
return
u = max(P | X, key=lambda w: len(adj[w] & P))
for v in list(P - adj[u]):
bk(R + 1, P & adj[v], X & adj[v])
P = P - {v}
X = X | {v}
bk(0, set(range(N)), set())
return best
def all_cliques(N, edges):
"""Yield every clique (as a tuple) including the empty one."""
adj = [set() for _ in range(N)]
for x, y in edges:
adj[x].add(y)
adj[y].add(x)
def rec(cur, cand):
yield tuple(cur)
for v in sorted(cand):
yield from rec(cur + [v], {w for w in cand if w > v and w in adj[v]})
yield from rec([], set(range(N)))
# ------------------------------------------------------------ instance sets
def contact_sets(n, max_m=None):
pairs = list(itertools.combinations(range(n), 2))
top = len(pairs) if max_m is None else min(max_m, len(pairs))
for m in range(top + 1):
for c in itertools.combinations(pairs, m):
yield list(c)
def instances_exhaustive():
# all contact-set pairs through 3 vertices (both sides), n >= 1
for n1 in range(1, 4):
for n2 in range(1, 4):
for E1 in contact_sets(n1):
for E2 in contact_sets(n2):
yield n1, E1, n2, E2
# up to 5 vertices, at most 2 contacts on each side
for n1 in range(1, 6):
for n2 in range(1, 6):
if n1 <= 3 and n2 <= 3:
continue
for E1 in contact_sets(n1, 2):
for E2 in contact_sets(n2, 2):
yield n1, E1, n2, E2
def sampled_instances(count, nmax=7, mmax=6):
for _ in range(count):
n1, n2 = random.randint(2, nmax), random.randint(2, nmax)
P1 = list(itertools.combinations(range(n1), 2))
P2 = list(itertools.combinations(range(n2), 2))
m1, m2 = random.randint(0, min(mmax, len(P1))), random.randint(0, min(mmax, len(P2)))
yield n1, sorted(random.sample(P1, m1)), n2, sorted(random.sample(P2, m2))
# ------------------------------------------------------ per-instance checks
def check_instance(n1, E1, n2, E2, full_cliques=True):
m1, m2 = len(E1), len(E2)
cands, edges = reduce(n1, E1, n2, E2)
N = len(cands)
als = alignments(n1, n2)
scores = {tuple(c): score(c, E1, E2) for c in als}
opt = max(scores.values())
# --- section 4: size vs formula
check("4 size", N == m1 * m2, f"vertices {N} != {m1*m2}")
check("4 size", len(edges) <= N * (N - 1) // 2, "edges above bound")
# --- section 5: well-formedness
check("5 structure", all(x < y < N for x, y in edges), "edge not a normalised pair")
for (i, j, A) in cands:
(s1, t1), (s2, t2) = A
check("5 structure", s1 < s2 and t1 < t2, "candidate not locally increasing")
check("5 structure", 0 <= t1 < t2 < n2 and 0 <= s1 < s2 < n1, "candidate out of range")
index = {(c[0], c[1]) for c in cands}
check("5 structure", len(index) == N, "duplicate candidate")
check("5 structure", [(c[0], c[1]) for c in cands] ==
[(i, j) for i in range(m1) for j in range(m2)], "candidate order")
# largest clique bounded by min(m1, m2)
mc = max_clique_size(N, edges)
check("5 structure", mc <= min(m1, m2), "clique exceeds min(m1,m2)")
# --- section 2: optimum preserved both ways
# forward: every alignment yields a clique of size = score
for c in als:
sel = []
for ci, (u, v) in enumerate(E1):
if c[u] and c[v]:
a, b = c[u] - 1, c[v] - 1
if (a, b) in set(E2):
sel.append(ci * m2 + E2.index((a, b)))
ok = all((x, y) in edges for x, y in itertools.combinations(sorted(sel), 2))
check("2 forward/backward", ok and len(sel) == scores[tuple(c)],
f"forward clique fails {n1,E1,n2,E2,c}")
check("2 forward/backward", mc == opt, f"opt {opt} != max clique {mc} on {n1,E1,n2,E2}")
# --- section 3: enumerator self-test, rejection of non-cliques and wrong lengths
if N <= 10:
brute = set()
for mask in range(1 << N):
sel = [p for p in range(N) if mask >> p & 1]
if all((x, y) in edges for x, y in itertools.combinations(sel, 2)):
brute.add(tuple(sel))
enum = list(all_cliques(N, edges))
check("3 extraction", set(enum) == brute and len(enum) == len(brute),
"clique enumerator disagrees with brute force")
nonedges = [(x, y) for x in range(N) for y in range(x + 1, N) if (x, y) not in edges]
for (x, y) in nonedges[:25]:
bits = [p in (x, y) for p in range(N)]
check("3 extraction", extract_bits(n1, n2, cands, edges, bits) is None,
"non-clique accepted")
bits2 = [True] * N
check("3 extraction", extract_bits(n1, n2, cands, edges, bits2) is None,
"all-ones non-clique accepted")
for bad in (N + 1, max(N - 1, 0)) if N > 0 else (1,):
check("3 extraction", extract_bits(n1, n2, cands, edges, [False] * bad) is None,
"wrong length accepted")
check("3 extraction", extract_bits(n1, n2, cands, edges, [False] * N) == [0] * n1,
"empty selection")
# --- section 3: extraction from every clique
if full_cliques:
for K in all_cliques(N, edges):
cfg = extract(n1, n2, cands, edges, K)
check("3 extraction", cfg is not None, "clique rejected")
check("3 extraction", feasible(cfg, n2), f"extracted infeasible {cfg}")
sc = scores[tuple(cfg)]
check("3 extraction", sc >= len(K), f"score {sc} < size {len(K)}")
if len(K) == mc:
check("3 extraction", sc == opt, "maximum clique decodes to non-optimal")
# non-edges are rejected
else:
for K in [()]:
check("3 extraction", extract(n1, n2, cands, edges, K) == [0] * n1, "empty clique")
return mc, opt
# ------------------------------------------------------------------ sections
def section1():
"""Tie the declared overhead formulas to the construction's own index structure."""
m1, m2 = sp.symbols("m1 m2", integer=True, nonnegative=True)
declared_v = m1 * m2
declared_e = m1 * m2 * (m1 * m2 - 1) / 2
sharp = 2 * sp.binomial(m1, 2) * sp.binomial(m2, 2)
# Symbolic: pairs of vertices minus pairs sharing a contact index equals the sharp bound
shares_i = m1 * sp.binomial(m2, 2)
shares_j = m2 * sp.binomial(m1, 2)
check("1 symbolic",
sp.simplify(sp.expand_func(declared_e - shares_i - shares_j - sharp)) == 0,
"C(N,2) - sharing pairs != sharp bound")
P = [(0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3), (0, 4), (1, 4)]
for a_ in range(0, 7):
for b_ in range(0, 7):
c, e = reduce(5, P[:a_], 5, P[:b_])
vs = {m1: a_, m2: b_}
# measured |V| and |E| against the declared formulas, from the built target
check("1 symbolic", len(c) == int(declared_v.subs(vs)), f"|V| {a_},{b_}")
check("1 symbolic", len(e) <= declared_e.subs(vs), f"|E| bound {a_},{b_}")
check("1 symbolic", len(e) <= sharp.subs(vs), f"|E| sharp {a_},{b_}")
# the construction's own index structure: vertices differing in both indices
both = sum(1 for p in range(len(c)) for q in range(p + 1, len(c))
if c[p][0] != c[q][0] and c[p][1] != c[q][1])
check("1 symbolic", both == sharp.subs(vs), f"index pairs {a_},{b_}")
check("1 symbolic", all(c[p][0] != c[q][0] and c[p][1] != c[q][1] for p, q in e),
"edge between vertices sharing an index")
def section_exhaustive():
n = 0
for inst in instances_exhaustive():
n += 1
check_instance(*inst)
return n
def section_sampled():
for inst in sampled_instances(400):
check_instance(*inst, full_cliques=True)
# larger: only max clique vs optimum plus all-alignment forward check
for inst in sampled_instances(60, nmax=7, mmax=12):
check_instance(*inst, full_cliques=False)
def section_edge_exactness():
"""Section 4: the edge count is only a bound; show it can be strict and attained."""
c, e = reduce(4, [(0, 2), (1, 3), (0, 3)], 4, [(0, 1), (2, 3), (0, 3)])
check("4 size", len(c) == 9 and len(e) == 2 and len(e) < 36, "issue example size")
# complete case: m1 = m2 = 1, N = 1, bound 0 attained
c, e = reduce(2, [(0, 1)], 2, [(0, 1)])
check("4 size", len(c) == 1 and len(e) == 0, "single candidate")
# empty side -> empty graph
c, e = reduce(3, [], 3, [(0, 1)])
check("4 size", len(c) == 0 and len(e) == 0, "empty E1")
c, e = reduce(3, [(0, 1)], 3, [])
check("4 size", len(c) == 0 and len(e) == 0, "empty E2")
# attained bound: identical single-chain maps give compatible cliques; find max edges ratio
# Derivation: two candidates sharing a source contact or a target contact are never
# adjacent, so edges <= C(N,2) - m1*C(m2,2) - m2*C(m1,2) = 2*C(m1,2)*C(m2,2).
# Hence C(N,2) is attained only for N <= 1; test strictness and the sharper bound.
attained = tight = 0
for n1, E1, n2, E2 in instances_exhaustive():
if n1 > 3 or n2 > 3:
break
c, e = reduce(n1, E1, n2, E2)
N, m1, m2 = len(c), len(E1), len(E2)
cn = lambda x: x * (x - 1) // 2
sharper = 2 * cn(m1) * cn(m2)
check("4 size", len(e) <= sharper, "sharper bound violated")
check("4 size", N <= 1 or len(e) < cn(N), "C(N,2) attained with N >= 2")
check("4 size", sharper == cn(N) - m1 * cn(m2) - m2 * cn(m1), "identity")
attained += len(e) == cn(N)
tight += len(e) == sharper
return attained, tight
def section_special_cases():
# n1 > n2: no special case, alignment partial
n1, E1, n2, E2 = 5, [(0, 4), (1, 2), (3, 4)], 2, [(0, 1)]
mc, opt = check_instance(n1, E1, n2, E2)
check("5 structure", mc == opt == 1)
# shared source endpoints must agree
c, e = reduce(3, [(0, 1), (0, 2)], 3, [(0, 1), (0, 2), (1, 2)])
pid = {(x[0], x[1]): p for p, x in enumerate(c)}
check("5 structure", (pid[(0, 0)], pid[(1, 2)]) not in e and (pid[(0, 0)], pid[(1, 0)]) not in e
and (pid[(0, 0)], pid[(1, 1)]) in e)
# disjoint source contacts cannot reuse target endpoints
c, e = reduce(4, [(0, 1), (2, 3)], 3, [(0, 1), (1, 2)])
pid = {(x[0], x[1]): p for p, x in enumerate(c)}
check("5 structure", (pid[(0, 0)], pid[(1, 1)]) not in e)
# crossing assignments excluded
c, e = reduce(4, [(0, 2), (1, 3)], 4, [(0, 3), (1, 2)])
pid = {(x[0], x[1]): p for p, x in enumerate(c)}
check("5 structure", (pid[(0, 0)], pid[(1, 1)]) not in e)
# witness-value caveat: nonmaximum clique decodes to higher score
n1, E1, n2, E2 = 3, [(0, 1), (1, 2), (0, 2)], 3, [(0, 1), (1, 2), (0, 2)]
c, e = reduce(n1, E1, n2, E2)
p01 = [p for p, x in enumerate(c) if (x[0], x[1]) == (0, 0)][0]
p12 = [p for p, x in enumerate(c) if (x[0], x[1]) == (1, 1)][0]
cfg = extract(n1, n2, c, e, {p01, p12})
check("3 extraction", cfg == [1, 2, 3], f"cfg {cfg}")
check("3 extraction", score(cfg, E1, E2) == 3 > 2, "caveat: score should exceed clique size")
check("3 extraction", max_clique_size(len(c), e) == 3)
# non-clique selection is rejected, wrong-order pair is not an edge
nonedge = next((x, y) for x in range(len(c)) for y in range(x + 1, len(c)) if (x, y) not in e)
check("3 extraction", extract(3, 3, c, e, set(nonedge)) is None, "non-clique accepted")
# an unordered (graph-isomorphism) construction would overestimate: mutation test
ex = (4, [(0, 2), (1, 3), (0, 3)], 4, [(0, 1), (2, 3), (0, 3)])
check("7 NO", mutant_unordered_opt(*ex) == 3, "mutant should overestimate (detects order loss)")
def mutant_unordered_opt(n1, E1, n2, E2):
"""Deliberately wrong: injective partial maps with no order restriction."""
S2 = set(E2)
best = 0
for k in range(0, min(n1, n2) + 1):
for src in itertools.combinations(range(n1), k):
for img in itertools.permutations(range(n2), k):
f = dict(zip(src, img))
sc = sum(1 for u, v in E1 if u in f and v in f
and (min(f[u], f[v]), max(f[u], f[v])) in S2)
best = max(best, sc)
return best
def section_yes_no():
n1, E1, n2, E2 = 4, [(0, 2), (1, 3), (0, 3)], 4, [(0, 1), (2, 3), (0, 3)]
c, e = reduce(n1, E1, n2, E2)
expect_assign = [
((0, 0), (2, 1)), ((0, 2), (2, 3)), ((0, 0), (2, 3)),
((1, 0), (3, 1)), ((1, 2), (3, 3)), ((1, 0), (3, 3)),
((0, 0), (3, 1)), ((0, 2), (3, 3)), ((0, 0), (3, 3)),
]
check("6 YES", [x[2] for x in c] == expect_assign, "nine candidates")
check("6 YES", sorted(e) == [(0, 8), (4, 8)], f"edges {sorted(e)}")
check("6 YES", len(c) == 9 and len(e) == 2)
cfg = extract(n1, n2, c, e, {0, 8})
check("6 YES", cfg == [1, 0, 2, 4], f"cfg {cfg}")
check("6 YES", feasible(cfg, n2) and score(cfg, E1, E2) == 2)
check("6 YES", extract(n1, n2, c, e, {0, 4}) is None, "{0,4} is not a clique")
check("6 YES", extract(n1, n2, c, e, {0}) == [1, 0, 2, 0])
check("6 YES", extract(n1, n2, c, e, {4}) == [0, 3, 0, 4])
check("6 YES", extract(n1, n2, c, e, set()) == [0, 0, 0, 0])
check("6 YES", max_clique_size(9, e) == 2)
check("6 YES", max(score(a, E1, E2) for a in alignments(4, 4)) == 2)
# NO: threshold 3 unreachable on both sides
check("7 NO", max_clique_size(9, e) < 3, "target has a clique of size 3")
check("7 NO", all(score(a, E1, E2) < 3 for a in alignments(4, 4)), "source reaches 3")
check("7 NO", not any(len(K) >= 3 for K in all_cliques(9, e)), "3-clique exists")
# all-three-contacts alignment would have to be total => identity
tot = [a for a in alignments(4, 4) if all(a)]
check("7 NO", tot == [[1, 2, 3, 4]], "total alignments")
check("7 NO", score([1, 2, 3, 4], E1, E2) == 1)
# unordered graphs are isomorphic paths yet overlap is 2
check("7 NO", mutant_unordered_opt(n1, E1, n2, E2) == 3)
# infeasible target config: wrong-length / out-of-range / non-increasing configs
check("7 NO", not feasible([2, 1, 0, 0], 4) and not feasible([1, 5, 0, 0], 4))
check("7 NO", not feasible([1, 1, 0, 0], 4))
def section_benchmark():
root = os.path.expanduser("<benchmark>/MaximumContactMapOverlap")
if not os.path.isdir(root):
return 0
n = 0
for fn in sorted(os.listdir(root)):
if not fn.endswith(".json"):
continue
d = json.load(open(os.path.join(root, fn)))["data"]
n1, n2 = d["num_vertices_1"], d["num_vertices_2"]
E1, E2 = canon(d["contacts_1"]), canon(d["contacts_2"])
cands, edges = reduce(n1, E1, n2, E2)
N = len(cands)
mc = max_clique_size(N, edges)
opt = max(score(a, E1, E2) for a in alignments(n1, n2))
check("2 forward/backward", N == len(E1) * len(E2) and mc == opt,
f"{fn}: clique {mc} vs opt {opt}")
check("4 size", len(edges) <= N * (N - 1) // 2)
n += 1
return n
# ---------------------------------------------------------------- hypothesis
@st.composite
def instance_small(draw):
n1, n2 = draw(st.integers(1, 6)), draw(st.integers(1, 6))
P1 = list(itertools.combinations(range(n1), 2))
P2 = list(itertools.combinations(range(n2), 2))
E1 = draw(st.lists(st.sampled_from(P1), unique=True, max_size=4)) if P1 else []
E2 = draw(st.lists(st.sampled_from(P2), unique=True, max_size=4)) if P2 else []
return n1, sorted(E1), n2, sorted(E2)
@st.composite
def instance_dense_equal(draw):
n = draw(st.integers(3, 6))
P = list(itertools.combinations(range(n), 2))
E1 = draw(st.lists(st.sampled_from(P), unique=True, min_size=1, max_size=5))
# target = noisy order-preserving copy of source under a random increasing map
img = sorted(draw(st.lists(st.integers(0, n + 1), unique=True, min_size=n, max_size=n)))
E2 = {(img[u], img[v]) for u, v in E1}
extra = draw(st.lists(st.sampled_from(list(itertools.combinations(range(n + 2), 2))),
unique=True, max_size=2))
E2 |= set(extra)
return n, sorted(E1), n + 2, sorted(E2)
def section_hypothesis():
prof = dict(max_examples=250, deadline=None, derandomize=True,
suppress_health_check=list(HealthCheck))
@settings(**prof)
@given(instance_small())
def t1(inst):
check_instance(*inst, full_cliques=True)
@settings(**prof)
@given(instance_dense_equal())
def t2(inst):
n1, E1, n2, E2 = inst
mc, opt = check_instance(*inst, full_cliques=False)
# planted copy: optimum at least |E1|
check("2 forward/backward", opt >= len(E1), "planted alignment not found")
t1()
t2()
# ---------------------------------------------------------------------- main
def main():
section1()
n_ex = section_exhaustive()
section_sampled()
attained = section_edge_exactness()
section_special_cases()
section_yes_no()
n_benchmark = section_benchmark()
section_hypothesis()
print(f"exhaustive instances: {n_ex}; benchmark instances: {n_benchmark}; "
f"instances attaining C(N,2): {attained[0]}, attaining 2C(m1,2)C(m2,2): {attained[1]}")
names = ["1 symbolic", "2 forward/backward", "3 extraction", "4 size",
"5 structure", "6 YES", "7 NO"]
for s in names:
print(f" section {s:<20} checks={COUNTS.get(s, 0):>9} failures={FAILS.get(s, 0)}")
if COUNTS.get(s, 0) == 0:
FAILS[s] = FAILS.get(s, 0) + 1 # an empty section is a failure
for m in FAIL_MSGS:
print("FAIL", m)
total = sum(COUNTS.values())
fails = sum(FAILS.values())
print(f"TOTAL: {total} checks, {fails} failures")
return 1 if fails else 0
if __name__ == "__main__":
sys.exit(main())Independent check script (Python)#!/usr/bin/env python3
"""Independent check for MaximumContactMapOverlap -> MaximumClique<SimpleGraph, One>.
Written from the proof alone (the constructor script was not read before
this file was written). Implements the construction, the extraction and both
feasibility checkers from scratch, then attacks every claim of the proof:
vertex count, both edge bounds, the "edges only between pairs differing in both
indices" lemma, (=>) and (<=) directions, omega <= min(m1, m2), extraction
validity, sigma(f) >= |K| for every clique, optimality of every maximum clique,
rejection of non-cliques and wrong lengths, both worked examples and the
triangle remark.
Usage: python3 independent check.py (full run, prints audit and TOTAL line)
python3 independent check.py --quick (reduced run, used by mutants.py)
The module exposes `run(reduce_fn, extract_fn, quick)` so that mutants.py can
inject broken constructions.
"""
from __future__ import annotations
import itertools
import json
import os
import random
import sys
from collections import Counter
from math import comb
BENCH = "<benchmark>/MaximumContactMapOverlap"
SEED = 1205
CLIQUE_CAP = 200_000
SHARP_FACTOR = 2 # the proof's sharper edge bound is SHARP_FACTOR * C(m1,2) * C(m2,2)
# --------------------------------------------------------------------------
# Construction and extraction (from the proof, independently written)
# --------------------------------------------------------------------------
def sgn(x: int) -> int:
return (x > 0) - (x < 0)
def assignment_sets(inst):
"""A_{ij} = {(u,a),(v,b)} for e_i=(u,v), d_j=(a,b), index i*m2+j."""
_n1, E1, _n2, E2 = inst
return [((u, a), (v, b)) for (u, v) in E1 for (a, b) in E2]
def compatible(Ap, Aq) -> bool:
for (s, t) in Ap:
for (s2, t2) in Aq:
if sgn(s - s2) != sgn(t - t2):
return False
return True
def reduce(inst):
"""Return (num_vertices, frozenset of edges (p,q) with p<q)."""
A = assignment_sets(inst)
N = len(A)
edges = set()
for p in range(N):
for q in range(p + 1, N):
if compatible(A[p], A[q]):
edges.add((p, q))
return N, frozenset(edges)
class ExtractionError(ValueError):
pass
def extract(inst, target, sel):
"""Selection (list of bools of length N) -> source config; error if not a clique."""
n1, E1, n2, E2 = inst
N, edges = target
if len(sel) != N:
raise ExtractionError("wrong length")
K = [p for p, x in enumerate(sel) if x]
for p, q in itertools.combinations(K, 2):
if (p, q) not in edges:
raise ExtractionError("not a clique")
A = assignment_sets(inst)
f = [0] * n1
for p in K:
for (s, t) in A[p]:
f[s] = t + 1
return f
# --------------------------------------------------------------------------
# Independent oracles
# --------------------------------------------------------------------------
def source_feasible(inst, f) -> bool:
n1, _E1, n2, _E2 = inst
if len(f) != n1:
return False
if any((not isinstance(x, int)) or x < 0 or x > n2 for x in f):
return False
nz = [x for x in f if x]
return all(a < b for a, b in zip(nz, nz[1:]))
def score(inst, f) -> int:
_n1, E1, _n2, E2 = inst
S2 = set(E2)
c = 0
for (u, v) in E1:
if f[u] and f[v] and (f[u] - 1, f[v] - 1) in S2:
c += 1
return c
def alignments(n1, n2):
"""All increasing partial injections, as configs."""
for k in range(min(n1, n2) + 1):
for src in itertools.combinations(range(n1), k):
for dst in itertools.combinations(range(n2), k):
f = [0] * n1
for s, t in zip(src, dst):
f[s] = t + 1
yield f
def adjacency(target):
N, edges = target
adj = [0] * N
for p, q in edges:
adj[p] |= 1 << q
adj[q] |= 1 << p
return adj
def all_cliques(target):
"""Every clique (including empty) as a sorted tuple."""
N, _ = target
adj = adjacency(target)
out = []
def ext(cur, cand):
out.append(tuple(cur))
if len(out) > CLIQUE_CAP:
raise RuntimeError("clique enumeration exceeded cap (graph far denser than the bound allows)")
while cand:
low = cand & -cand
v = low.bit_length() - 1
cand ^= low
ext(cur + [v], cand & adj[v])
ext([], (1 << N) - 1)
return out
def is_clique(target, K) -> bool:
_N, edges = target
return all((p, q) in edges for p, q in itertools.combinations(sorted(K), 2))
def omega_bruteforce(target) -> int:
"""Subset brute force, independent of all_cliques (only for tiny N)."""
N, _ = target
best = 0
for r in range(N, 0, -1):
if r <= best:
break
for K in itertools.combinations(range(N), r):
if is_clique(target, K):
return r
return best
# --------------------------------------------------------------------------
# Check machinery
# --------------------------------------------------------------------------
class Audit:
def __init__(self):
self.count = Counter()
self.fail = Counter()
self.first = {}
def check(self, name, ok, ctx=None):
self.count[name] += 1
if not ok:
self.fail[name] += 1
self.first.setdefault(name, ctx)
return ok
def safe(fn, *a):
try:
return fn(*a), None
except ExtractionError as e:
return None, e
except Exception as e: # a mutant may crash; a crash is a failure, not an abort
return None, e
def check_instance(audit, inst, reduce_fn, extract_fn, rng, full=True):
n1, E1, n2, E2 = inst
m1, m2 = len(E1), len(E2)
try:
target = reduce_fn(inst)
except Exception as e:
audit.check("reduce_runs", False, (inst, repr(e)))
return
audit.check("reduce_runs", True)
N, edges = target
# --- overhead / structure
audit.check("vertices_exact_m1m2", N == m1 * m2, inst)
audit.check("edges_le_C(N,2)", len(edges) <= comb(N, 2), inst)
audit.check("edges_le_2C(m1,2)C(m2,2)", len(edges) <= SHARP_FACTOR * comb(m1, 2) * comb(m2, 2), inst)
audit.check("C(N,2)_not_attained_for_N>=2", N < 2 or len(edges) < comb(N, 2), inst)
audit.check("edges_well_formed",
all(0 <= p < q < N for p, q in edges), inst)
if m2:
audit.check("edges_differ_in_both_indices",
all(p // m2 != q // m2 and p % m2 != q % m2 for p, q in edges), inst)
# --- (=>): every alignment yields a clique of size sigma(f)
idx2 = {d: j for j, d in enumerate(E2)}
opt = 0
for f in alignments(n1, n2):
sig = score(inst, f)
opt = max(opt, sig)
if not full and rng.random() > 0.3:
continue
K = [i * m2 + idx2[(f[u] - 1, f[v] - 1)]
for i, (u, v) in enumerate(E1)
if f[u] and f[v] and (f[u] - 1, f[v] - 1) in idx2]
audit.check("forward_alignment_gives_clique",
len(set(K)) == sig and all(k < N for k in K) and is_clique(target, K),
(inst, f))
# --- cliques
try:
cl = all_cliques(target)
except RuntimeError as e:
audit.check("clique_enumeration_tractable", False, (inst, str(e)))
return
audit.check("clique_enumeration_tractable", True)
omega = max(len(K) for K in cl)
audit.check("omega_equals_OPT", omega == opt, (inst, omega, opt))
audit.check("omega_le_min(m1,m2)", omega <= min(m1, m2), inst)
if N <= 12:
audit.check("omega_matches_subset_bruteforce", omega_bruteforce(target) == omega, inst)
for K in cl:
if m2:
audit.check("clique_distinct_source_and_target_contacts",
len({p // m2 for p in K}) == len(K) and len({p % m2 for p in K}) == len(K),
(inst, K))
sel = [False] * N
for p in K:
sel[p] = True
f, err = safe(extract_fn, inst, target, sel)
if not audit.check("extract_clique_accepted", err is None, (inst, K, repr(err))):
continue
if len(K) == 0:
audit.check("empty_clique_all_unmatched", f == [0] * n1, (inst, f))
feas = source_feasible(inst, f)
audit.check("extract_feasible", feas, (inst, K, f))
if not feas:
continue
sig = score(inst, f)
audit.check("extract_score_ge_clique_size", sig >= len(K), (inst, K, f))
if len(K) == omega:
audit.check("max_clique_extracts_optimum", sig == opt, (inst, K, f, opt))
# --- rejection
_, err = safe(extract_fn, inst, target, [False] * (N + 1))
audit.check("extract_rejects_long", err is not None, inst)
if N:
_, err = safe(extract_fn, inst, target, [False] * (N - 1))
audit.check("extract_rejects_short", err is not None, inst)
noncliques = []
for p in range(N):
for q in range(p + 1, N):
if (p, q) not in edges:
noncliques.append((p, q))
rng.shuffle(noncliques)
for p, q in noncliques[:5]:
sel = [False] * N
sel[p] = sel[q] = True
# widen with a random extra vertex sometimes
if N > 2 and rng.random() < 0.5:
sel[rng.randrange(N)] = True
_, err = safe(extract_fn, inst, target, sel)
audit.check("extract_rejects_nonclique", err is not None, (inst, p, q))
# --------------------------------------------------------------------------
# Instance families
# --------------------------------------------------------------------------
def contact_sets(n):
pairs = list(itertools.combinations(range(n), 2))
for mask in range(1 << len(pairs)):
yield tuple(pairs[k] for k in range(len(pairs)) if mask >> k & 1)
def exhaustive(quick):
"""All contact-set pairs with n1,n2 <= 4; plus n=5 on one side x n<=3 on the other."""
lim = 3 if quick else 4
for n1 in range(lim + 1):
for n2 in range(lim + 1):
for E1 in contact_sets(n1):
for E2 in contact_sets(n2):
yield (n1, E1, n2, E2)
if quick:
return
for n5 in [5]:
for n_small in range(4):
for E5 in contact_sets(n5):
for Es in contact_sets(n_small):
yield (n5, E5, n_small, Es)
yield (n_small, Es, n5, E5)
def shuffled(inst, rng):
n1, E1, n2, E2 = inst
E1, E2 = list(E1), list(E2)
rng.shuffle(E1)
rng.shuffle(E2)
return (n1, tuple(E1), n2, tuple(E2))
def load_benchmark():
out = []
if not os.path.isdir(BENCH):
return out
for name in sorted(os.listdir(BENCH)):
d = json.load(open(os.path.join(BENCH, name)))["data"]
out.append((name, (d["num_vertices_1"], tuple(tuple(sorted(e)) for e in d["contacts_1"]),
d["num_vertices_2"], tuple(tuple(sorted(e)) for e in d["contacts_2"]))))
return out
# --------------------------------------------------------------------------
# Worked examples
# --------------------------------------------------------------------------
YES = (4, ((0, 2), (1, 3), (0, 3)), 4, ((0, 1), (2, 3), (0, 3)))
TRIANGLE = (3, ((0, 1), (0, 2), (1, 2)), 3, ((0, 1), (0, 2), (1, 2)))
def check_examples(audit, reduce_fn, extract_fn):
inst = YES
t = reduce_fn(inst)
N, edges = t
audit.check("yes_num_vertices_9", N == 9)
table = [((0, 0), (2, 1)), ((0, 2), (2, 3)), ((0, 0), (2, 3)),
((1, 0), (3, 1)), ((1, 2), (3, 3)), ((1, 0), (3, 3)),
((0, 0), (3, 1)), ((0, 2), (3, 3)), ((0, 0), (3, 3))]
audit.check("yes_assignment_table", assignment_sets(inst) == table)
audit.check("yes_edges_[(0,8),(4,8)]", sorted(edges) == [(0, 8), (4, 8)], sorted(edges))
audit.check("yes_bounds_36_18", comb(9, 2) == 36 and SHARP_FACTOR * comb(3, 2) ** 2 == 18 and len(edges) <= 18)
audit.check("yes_0_4_not_adjacent", (0, 4) not in edges)
def ex(K):
sel = [p in K for p in range(N)]
f, err = safe(extract_fn, inst, t, sel)
return f
f = ex({0, 8})
audit.check("yes_clique_0_8_extracts_[1,0,2,4]", f == [1, 0, 2, 4], f)
audit.check("yes_score_2", f is not None and score(inst, f) == 2)
opt = max(score(inst, g) for g in alignments(4, 4))
audit.check("yes_optimum_2", opt == 2, opt)
audit.check("yes_clique_0", ex({0}) == [1, 0, 2, 0] and score(inst, [1, 0, 2, 0]) == 1)
audit.check("yes_clique_4", ex({4}) == [0, 3, 0, 4] and score(inst, [0, 3, 0, 4]) == 1)
audit.check("yes_empty", ex(set()) == [0, 0, 0, 0])
# oracle self-tests (so a broken oracle cannot hide a broken construction)
audit.check("oracle_feasibility_sanity",
source_feasible(inst, [0, 0, 0, 0]) and source_feasible(inst, [1, 2, 3, 4])
and not source_feasible(inst, [2, 1, 0, 0]) and not source_feasible(inst, [1, 1, 0, 0])
and not source_feasible(inst, [1, 5, 0, 0]) and not source_feasible(inst, [1, 0, 0]))
audit.check("oracle_clique_enumeration_sanity",
len(all_cliques((3, frozenset({(0, 1), (0, 2), (1, 2)})))) == 8
and len(all_cliques(t)) == 1 + N + len(edges))
# NO example: threshold 3
cl = all_cliques(t) # cap raises on a blown-up target; run() records it as a failure
audit.check("no_no_3_clique", max(len(K) for K in cl) == 2)
audit.check("no_no_alignment_score_3", all(score(inst, g) < 3 for g in alignments(4, 4)))
audit.check("no_identity_scores_1", score(inst, [1, 2, 3, 4]) == 1)
# order-ignoring overlap would be 3 (the two contact graphs are isomorphic)
best_unordered = 0
for perm in itertools.permutations(range(4)):
f = [p + 1 for p in perm]
s = sum(1 for (u, v) in inst[1] if tuple(sorted((perm[u], perm[v]))) in set(inst[3]))
best_unordered = max(best_unordered, s)
audit.check("no_order_ignoring_would_give_3", best_unordered == 3)
# triangle remark: {p00, p11} extracts [1,2,3] with score 3 > 2
tt = reduce_fn(TRIANGLE)
Nt = tt[0]
sel = [p in (0, 4) for p in range(Nt)]
f, err = safe(extract_fn, TRIANGLE, tt, sel)
audit.check("triangle_p00_p11_extracts_[1,2,3]_score_3",
err is None and f == [1, 2, 3] and score(TRIANGLE, f) == 3, (f, err))
# --------------------------------------------------------------------------
# Hypothesis strategies aimed at the special cases
# --------------------------------------------------------------------------
def run_hypothesis(audit, reduce_fn, extract_fn, max_examples):
from hypothesis import given, settings, strategies as st, HealthCheck
rng = random.Random(SEED + 7)
common = settings(max_examples=max_examples, deadline=None, derandomize=True,
suppress_health_check=list(HealthCheck), database=None)
def side(draw, n, max_m, kind):
pairs = list(itertools.combinations(range(n), 2))
if not pairs:
return ()
if kind == "dense":
E = draw(st.lists(st.sampled_from(pairs), unique=True, max_size=max_m))
elif kind == "star": # many contacts sharing endpoints (equality preservation)
hubs = draw(st.lists(st.integers(0, n - 1), min_size=1, max_size=2, unique=True))
E = [p for p in pairs if p[0] in hubs or p[1] in hubs]
E = draw(st.permutations(E))[:max_m]
elif kind == "arcs": # crossing / nested arcs, protein-like backbone
E = [(i, i + 1) for i in range(n - 1)] if draw(st.booleans()) else []
extra = draw(st.lists(st.sampled_from(pairs), unique=True, max_size=3))
E = list(dict.fromkeys(E + extra))[:max_m]
else:
E = []
return tuple(draw(st.permutations(E)))
@st.composite
def cmo(draw, kinds, n_range):
n1 = draw(st.integers(*n_range))
n2 = draw(st.integers(*n_range))
k1 = draw(st.sampled_from(kinds))
k2 = draw(st.sampled_from(kinds))
return (n1, side(draw, n1, 7, k1), n2, side(draw, n2, 7, k2))
@st.composite
def lopsided(draw): # n1 > n2, tiny target, or one empty side
n2 = draw(st.integers(0, 3))
n1 = draw(st.integers(n2, 7))
k1 = draw(st.sampled_from(["dense", "star", "arcs", "empty"]))
k2 = draw(st.sampled_from(["dense", "star", "empty"]))
inst = (n1, side(draw, n1, 7, k1), n2, side(draw, n2, 3, k2))
if draw(st.booleans()):
inst = (inst[2], inst[3], inst[0], inst[1])
return inst
seen = Counter()
for name, strat in [("dense", cmo(["dense"], (0, 7))),
("shared_endpoint", cmo(["star"], (2, 7))),
("crossing_nested", cmo(["arcs", "dense"], (2, 7))),
("lopsided_or_empty", lopsided())]:
@common
@given(strat)
def prop(inst):
seen[name] += 1
check_instance(audit, inst, reduce_fn, extract_fn, rng)
prop()
return seen
# --------------------------------------------------------------------------
def run(reduce_fn=reduce, extract_fn=extract, quick=False, verbose=True):
audit = Audit()
rng = random.Random(SEED)
sections = {}
def section(name, fn):
before = sum(audit.count.values())
fn()
sections[name] = sum(audit.count.values()) - before
n_ex = Counter()
def s_exh():
for inst in exhaustive(quick):
n_ex["instances"] += 1
check_instance(audit, inst, reduce_fn, extract_fn, rng, full=not quick)
section("exhaustive", s_exh)
def s_shuf():
r = random.Random(SEED + 1)
base = [(a, E, b, F) for a in range(5) for b in range(5)
for E in contact_sets(a) for F in contact_sets(b) if len(E) >= 2 or len(F) >= 2]
for inst in r.sample(base, 300 if quick else 1500):
n_ex["shuffled"] += 1
check_instance(audit, shuffled(inst, r), reduce_fn, extract_fn, rng)
section("shuffled_contact_order", s_shuf)
def s_rand():
r = random.Random(SEED + 2)
for _ in range(30 if quick else 400):
n1, n2 = r.randint(5, 7), r.randint(5, 7)
p1, p2 = r.random(), r.random()
E1 = tuple(e for e in itertools.combinations(range(n1), 2) if r.random() < p1 * 0.5)[:8]
E2 = tuple(e for e in itertools.combinations(range(n2), 2) if r.random() < p2 * 0.5)[:8]
n_ex["random"] += 1
check_instance(audit, shuffled((n1, E1, n2, E2), r), reduce_fn, extract_fn, rng, full=False)
section("random_n5to7", s_rand)
def s_benchmark():
for _name, inst in load_benchmark()[: (5 if quick else 100)]:
n_ex["benchmark"] += 1
check_instance(audit, inst, reduce_fn, extract_fn, rng, full=False)
section("benchmark", s_benchmark)
def s_examples():
try:
check_examples(audit, reduce_fn, extract_fn)
except Exception as e:
audit.check("examples_run", False, repr(e))
section("worked_examples", s_examples)
def s_overhead():
# the proof's count: pairs differing in both indices = C(N,2) - m1 C(m2,2) - m2 C(m1,2)
# = 2 C(m1,2) C(m2,2); counted directly, and symbolically
for m1 in range(0, 9):
for m2 in range(0, 9):
idx = [(i, j) for i in range(m1) for j in range(m2)]
both = sum(1 for (i, j), (k, l) in itertools.combinations(idx, 2) if i != k and j != l)
N = m1 * m2
audit.check("overhead_pair_count_identity",
both == comb(N, 2) - m1 * comb(m2, 2) - m2 * comb(m1, 2)
== SHARP_FACTOR * comb(m1, 2) * comb(m2, 2))
audit.check("overhead_C(N,2)_strictly_above_sharp_iff_N>=2",
(comb(N, 2) > SHARP_FACTOR * comb(m1, 2) * comb(m2, 2)) == (N >= 2))
import sympy as sp
a, b = sp.symbols("m1 m2", integer=True, nonnegative=True)
C2 = lambda x: x * (x - 1) / 2
audit.check("overhead_symbolic_identity",
sp.expand(C2(a * b) - a * C2(b) - b * C2(a) - SHARP_FACTOR * C2(a) * C2(b)) == 0)
section("overhead_formulas", s_overhead)
hyp = {}
def s_hyp():
hyp.update(run_hypothesis(audit, reduce_fn, extract_fn, 40 if quick else 300))
section("hypothesis", s_hyp)
if verbose:
print("=== independent check: MaximumContactMapOverlap -> MaximumClique<SimpleGraph,One> ===")
print("instances:", dict(n_ex), "hypothesis:", dict(hyp))
print("--- checks per section")
for k, v in sections.items():
print(f" {k:28s} {v}")
print("--- checks per claim (count / failures)")
for k in sorted(audit.count):
print(f" {k:48s} {audit.count[k]:8d} {audit.fail[k]:6d}")
for k, ctx in audit.first.items():
print(" FIRST FAILURE", k, str(ctx)[:300])
print(f"TOTAL: {sum(audit.count.values())} checks, {sum(audit.fail.values())} failures")
return audit
if __name__ == "__main__":
a = run(quick="--quick" in sys.argv)
sys.exit(1 if sum(a.fail.values()) else 0)Cross-comparison script (Python)#!/usr/bin/env python3
"""Cross-compare MaximumContactMapOverlap -> MaximumClique<SimpleGraph, One> targets.
1. Script vs script: on every shared instance, the constructor's target
(check_a.py) and the independent check's target (check_b.py)
must be identical graphs. The documented isomorphism is the identity: vertex
p_{ij} has index i*m2 + j in both, i over contacts_1 and j over contacts_2 in
their stored (canonicalized, list) order. On a subset, both scripts' optimum
oracles must also agree (max clique = brute-force alignment optimum).
2. With --pred <binary>: each pred instance is written in the source model's JSON
format, reduced with `pred reduce --via <the single direct path>`, and the
bundle's target must be MaximumClique/SimpleGraph/One with unit weights and the
same vertex count and edge set as both scripts' targets.
Exit status 0 only if everything agrees; non-zero on any mismatch or pred failure.
python3 cross_compare.py # scripts only
python3 cross_compare.py --pred path/to/pred # scripts + pred
"""
from __future__ import annotations
import argparse
import contextlib
import importlib.util
import io
import itertools
import json
import os
import random
import subprocess
import sys
import tempfile
HERE = os.path.dirname(os.path.abspath(__file__))
SOURCE = "MaximumContactMapOverlap"
TARGET = "MaximumClique/SimpleGraph/One"
TARGET_VARIANT = {"graph": "SimpleGraph", "weight": "One"}
SEED = 12051
def load(name, path):
spec = importlib.util.spec_from_file_location(name, path)
mod = importlib.util.module_from_spec(spec)
with contextlib.redirect_stdout(io.StringIO()):
spec.loader.exec_module(mod)
return mod
con = load("constructor_cc", os.path.join(HERE, "constructor.py"))
adv = load("independent check_cc", os.path.join(HERE, "independent check", "independent check.py"))
def canonical(raw):
"""What the source model stores: (min, max) per contact, list order kept."""
return tuple((min(u, v), max(u, v)) for u, v in raw)
def con_target(inst):
n1, E1, n2, E2 = inst
cands, edges = con.reduce(n1, list(E1), n2, list(E2))
return len(cands), frozenset(edges)
def adv_target(inst):
return adv.reduce(inst)
# ------------------------------------------------------------- shared instances
def shared_instances():
"""(label, raw instance) - raw contacts may be non-canonical or shuffled."""
out = []
# constructor's exhaustive family
for n1, E1, n2, E2 in con.instances_exhaustive():
out.append(("constructor-exhaustive", (n1, tuple(E1), n2, tuple(E2))))
# independent check's exhaustive family
for inst in adv.exhaustive(False):
out.append(("independent check-exhaustive", inst))
# shuffled order and reversed orientation (exercise the stored order)
rng = random.Random(SEED)
small = [i for _, i in out if len(i[1]) >= 2 and len(i[3]) >= 2]
for inst in rng.sample(small, 1500):
n1, E1, n2, E2 = inst
E1, E2 = list(E1), list(E2)
rng.shuffle(E1)
rng.shuffle(E2)
E1 = [(v, u) if rng.random() < 0.5 else (u, v) for u, v in E1]
E2 = [(v, u) if rng.random() < 0.5 else (u, v) for u, v in E2]
out.append(("shuffled-reversed", (n1, tuple(E1), n2, tuple(E2))))
# benchmark (contacts as stored in the files, unsorted)
benchmark = os.path.expanduser("<benchmark>/MaximumContactMapOverlap")
if os.path.isdir(benchmark):
for fn in sorted(os.listdir(benchmark)):
d = json.load(open(os.path.join(benchmark, fn)))["data"]
out.append(("benchmark", (d["num_vertices_1"], tuple(map(tuple, d["contacts_1"])),
d["num_vertices_2"], tuple(map(tuple, d["contacts_2"])))))
# worked examples and special cases of the proof
special = [
adv.YES, adv.TRIANGLE,
(0, (), 0, ()), (0, (), 3, ((0, 2),)), (3, ((1, 2),), 0, ()), # empty sides
(5, ((0, 4), (1, 2), (3, 4)), 2, ((0, 1),)), # n1 > n2
(3, ((0, 1), (0, 2)), 3, ((0, 1), (0, 2), (1, 2))), # shared source endpoint
(4, ((0, 1), (2, 3)), 3, ((0, 1), (1, 2))), # disjoint vs shared target
(4, ((0, 2), (1, 3)), 4, ((0, 3), (1, 2))), # crossing vs nested
(6, ((5, 0), (4, 1), (3, 2)), 6, ((2, 3), (1, 4), (0, 5))), # reversed orientation
(7, ((0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6)), 7,
((0, 1), (1, 2), (2, 3), (3, 4), (4, 5), (5, 6))), # identical backbones
]
out += [("special", s) for s in special]
return out
def pred_subset(instances):
rng = random.Random(SEED + 1)
keep = [(l, i) for l, i in instances if l in ("benchmark", "special")]
keep += rng.sample([(l, i) for l, i in instances if l == "shuffled-reversed"], 150)
keep += rng.sample([(l, i) for l, i in instances if l.endswith("exhaustive")], 150)
return keep
# ---------------------------------------------------------------------- pred
def run(cmd, timeout=60):
p = subprocess.run(cmd, capture_output=True, text=True, timeout=timeout)
return p.returncode, p.stdout, p.stderr
def find_route(pred, workdir):
code, out, err = run([pred, "path", SOURCE, TARGET, "--limit", "all", "--quiet", "--json"])
if code != 0:
return None, f"pred path failed (exit {code}): {err.strip()[:300]}"
try:
batch = json.loads(out)
except json.JSONDecodeError as e:
return None, f"pred path printed no JSON: {e}"
direct = [p for p in batch.get("paths", []) if p.get("steps") == 1]
if len(direct) != 1:
return None, f"expected one direct {SOURCE} -> {TARGET} path, found {len(direct)}"
route = os.path.join(workdir, "route.json")
with open(route, "w") as f:
json.dump(direct[0], f, sort_keys=True)
return route, None
def source_json(inst):
n1, E1, n2, E2 = inst
return {"type": SOURCE, "variant": {},
"data": {"num_vertices_1": n1, "contacts_1": [list(e) for e in E1],
"num_vertices_2": n2, "contacts_2": [list(e) for e in E2]}}
def parse_target(bundle):
t = bundle["target"]
problems = []
if t.get("type") != "MaximumClique":
problems.append(f"target type {t.get('type')}")
if t.get("variant") != TARGET_VARIANT:
problems.append(f"target variant {t.get('variant')}")
g = t["data"]["graph"]
N = g["num_vertices"]
raw = [tuple(e) for e in g["edges"]]
edges = frozenset((min(e), max(e)) for e in raw)
if len(edges) != len(raw):
problems.append("duplicate edges")
if any(a == b for a, b in raw):
problems.append("self-loop")
w = t["data"].get("weights")
if w is not None and (len(w) != N or any(x != 1 for x in w)):
problems.append(f"weights not unit/length N: {w[:10]}")
return (N, edges), problems
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--pred", help="path to the pred binary")
ap.add_argument("--keep", help="directory to keep the written instances and bundles")
args = ap.parse_args()
instances = shared_instances()
by = {}
mism = []
opt_checked = 0
rng = random.Random(SEED + 2)
for label, raw in instances:
inst = (raw[0], canonical(raw[1]), raw[2], canonical(raw[3]))
a, c = adv_target(inst), con_target(inst)
by.setdefault(label, [0, 0])
by[label][0] += 1
if a != c:
by[label][1] += 1
mism.append((label, raw, "independent check != constructor", a, c))
continue
# optimum agreement on a subset (both scripts' independent oracles)
if label in ("special", "benchmark") or rng.random() < 0.05:
om_adv = max(len(K) for K in adv.all_cliques(a))
om_con = con.max_clique_size(*c)
opt_adv = max(adv.score(inst, f) for f in adv.alignments(inst[0], inst[2]))
opt_con = max(con.score(f, list(inst[1]), list(inst[3]))
for f in con.alignments(inst[0], inst[2]))
opt_checked += 1
if not (om_adv == om_con == opt_adv == opt_con):
by[label][1] += 1
mism.append((label, raw, "optimum disagreement", (om_adv, om_con), (opt_adv, opt_con)))
print("script vs script (independent check target == constructor target, identity labelling):")
for label, (n, bad) in by.items():
print(f" {label:24s} {n:7d} instances, {bad} mismatches")
print(f" optimum agreement (2 clique oracles, 2 alignment oracles) on {opt_checked} instances")
pred_bad = 0
pred_n = 0
if args.pred:
tmp = None if args.keep else tempfile.TemporaryDirectory(prefix="cmo_cc_")
work = args.keep or tmp.name
os.makedirs(work, exist_ok=True)
route, err = find_route(args.pred, work)
if err:
print(f"PRED FAILURE: {err}")
return 2
for k, (label, raw) in enumerate(pred_subset(instances)):
inst = (raw[0], canonical(raw[1]), raw[2], canonical(raw[3]))
expect = adv_target(inst)
assert expect == con_target(inst)
src = os.path.join(work, f"{k:04d}-src.json")
dst = os.path.join(work, f"{k:04d}-bundle.json")
with open(src, "w") as f:
json.dump(source_json(raw), f)
code, _out, perr = run([args.pred, "reduce", src, "--via", route, "-o", dst, "-q"])
pred_n += 1
if code != 0:
pred_bad += 1
mism.append((label, raw, f"pred reduce exit {code}: {perr.strip()[:200]}", None, expect))
continue
try:
got, problems = parse_target(json.load(open(dst)))
except Exception as e: # malformed bundle
pred_bad += 1
mism.append((label, raw, f"unreadable bundle: {e!r}", None, expect))
continue
if problems or got != expect:
pred_bad += 1
mism.append((label, raw, f"pred target differs {problems}",
(got[0], sorted(got[1])[:12]), (expect[0], sorted(expect[1])[:12])))
print(f"pred reduce vs both scripts: {pred_n} instances, {pred_bad} mismatches"
+ (f" (instances kept in {work})" if args.keep else ""))
if tmp:
tmp.cleanup()
else:
print("pred comparison skipped (no --pred given)")
for m in mism[:15]:
print("MISMATCH", m[0], m[1], m[2], "got", m[3], "expected", m[4])
total = sum(n for n, _ in by.values())
print(f"CROSS-COMPARE: {total} shared instances, {pred_n} via pred, {len(mism)} mismatches")
return 1 if mism else 0
if __name__ == "__main__":
sys.exit(main()) |
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #1235 +/- ##
========================================
Coverage 97.19% 97.20%
========================================
Files 1146 1148 +2
Lines 151558 151879 +321
========================================
+ Hits 147308 147632 +324
+ Misses 4250 4247 -3 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Fixes #1205
Adds the reduction
MaximumContactMapOverlap→MaximumClique<SimpleGraph, One>(order-preserving association graph): one vertex per pair of contacts, an edge between two pairs that differ in both coordinates and keep the order of the residue assignments consistent. Optimal cliques decode to optimal alignments.src/rules/maximumcontactmapoverlap_maximumclique.rs, registered insrc/rules/mod.rswith a canonical example.num_vertices = num_contacts_1 * num_contacts_2(exact);num_edgesis the upper boundm1(m1-1)m2(m2-1)/2.src/unit_tests/rules/maximumcontactmapoverlap_maximumclique.rs(closed-loop, rejection paths, overflow boundary).maximumcontactmapoverlap_ilp.rs: one-line disambiguation, since the source now has twoReduceToimpls.Local gate:
make checkandmake paperpass.🤖 Generated with Claude Code