ProCreations/repro-compact-conformal-subgraphs

1"""Corollary 1: parametric max-flow computation of the whole nested sequence,2certified against brute force, then timed. Replaces 2^n enumeration."""3import time, json, math, numpy as np4from collections import deque5EPS = 1e-9; INF = float("inf")6 7class Dinic:8    def __init__(s, N):9        s.N=N; s.to=[]; s.cap=[]; s.head=[-1]*N; s.nxt=[]; s.radj=[[] for _ in range(N)]; s.ops=010    def add(s,u,v,c):11        s.to.append(v); s.cap.append(c); s.nxt.append(s.head[u]); s.head[u]=len(s.to)-112        s.radj[v].append((u,len(s.to)-1))13        s.to.append(u); s.cap.append(0.0); s.nxt.append(s.head[v]); s.head[v]=len(s.to)-114        s.radj[u].append((v,len(s.to)-1))15    def bfs(s,src,snk):16        s.lv=[-1]*s.N; s.lv[src]=0; q=deque([src])17        while q:18            u=q.popleft(); e=s.head[u]19            while e!=-1:20                s.ops+=1; v=s.to[e]21                if s.cap[e]>EPS and s.lv[v]<0: s.lv[v]=s.lv[u]+1; q.append(v)22                e=s.nxt[e]23        return s.lv[snk]>=024    def dfs(s,u,snk,f):25        if u==snk: return f26        while s.it[u]!=-1:27            e=s.it[u]; v=s.to[e]; s.ops+=128            if s.cap[e]>EPS and s.lv[v]==s.lv[u]+1:29                d=s.dfs(v,snk,min(f,s.cap[e]))30                if d>EPS: s.cap[e]-=d; s.cap[e^1]+=d; return d31            s.it[u]=s.nxt[e]32        return 0.033    def maxflow(s,src,snk):34        fl=0.035        while s.bfs(src,snk):36            s.it=s.head[:]37            while True:38                f=s.dfs(src,snk,INF)39                if f<=EPS: break40                fl+=f41        return fl42    def reach_t(s,snk):43        seen=[False]*s.N; seen[snk]=True; q=deque([snk])44        while q:45            v=q.popleft()46            for u,eid in s.radj[v]:47                if not seen[u] and s.cap[eid]>EPS: seen[u]=True; q.append(u)48        return seen49 50def solve(n, edges, w, lam, stats):51    """maximal minimizer of |S| - lam*w(S), via max-closure min-cut."""52    m=len(edges); s=n+m; t=n+m+1; g=Dinic(n+m+2)53    for j,e in enumerate(edges):54        g.add(s,n+j,lam*float(w[j]))55        for v in e: g.add(n+j,v,INF)56    for v in range(n): g.add(v,t,1.0)57    g.maxflow(s,t); stats["flows"]+=1; stats["ops"]+=g.ops58    can=g.reach_t(t)59    return frozenset(v for v in range(n) if not can[v])60 61def wsum(S,edges,w): return float(sum(w[j] for j,e in enumerate(edges) if all(v in S for v in e)))62 63def sequence(n, edges, w, stats, lam_max=None):64    if lam_max is None: lam_max=2.0*(n+1)/max(1e-9,min(w))+1.065    lo,hi=0.0,lam_max66    Slo=solve(n,edges,w,lo,stats); Shi=solve(n,edges,w,hi,stats)67    out={Slo,Shi}68    def rec(a,b,Sa,Sb,d=0):69        if Sa==Sb or d>60: return70        wa,ca=wsum(Sa,edges,w),len(Sa); wb,cb=wsum(Sb,edges,w),len(Sb)71        if abs(wb-wa)<1e-12: return72        lam=(cb-ca)/(wb-wa)73        if not (a+1e-12 < lam < b-1e-12): return74        S=solve(n,edges,w,lam,stats); out.add(S)75        if S!=Sa and S!=Sb: rec(a,lam,Sa,S,d+1); rec(lam,b,S,Sb,d+1)76    rec(lo,hi,Slo,Shi)77    return sorted(out,key=len)78 79# ---- brute force reference (the logbook's original primitive) --------------80def brute(n, edges, w):81    """exact reference: walk the lower envelope of |S| - lam*w(S) over all 2^n masks.82    Only O(gamma) probes are needed, so this stays exact without the O(4^n) crossing scan."""83    M=1<<n; masks=np.arange(M)84    bits=np.zeros(M,dtype=np.int64)85    for v in range(n): bits+=(masks>>v)&186    bym=np.zeros(M)87    for j,e in enumerate(edges):88        em=089        for v in e: em|=1<<v90        bym+=w[j]*((masks&em)==em)91    def argmax_set(lam):92        obj=bits-lam*bym; best=obj.min()93        return int(np.bitwise_or.reduce(masks[obj<=best+2e-9]))94    bps=[]; lam=0.095    for _ in range(4*n+8):96        u=argmax_set(lam); w0=bym[u]; b0=bits[u]97        gt=bym>w0+1e-1298        if not gt.any(): break99        cand=(bits[gt]-b0)/(bym[gt]-w0)100        cand=cand[cand>lam+1e-9]101        if cand.size==0: break102        lam=float(cand.min()); bps.append(lam)103    probes=[0.0]+bps104    for a,b in zip(bps,bps[1:]): probes.append((a+b)/2)105    probes.append((bps[-1] if bps else 0.0)+1.0)106    seq=[]107    for lam in sorted(set(probes)):108        u=argmax_set(lam); S=frozenset(v for v in range(n) if u>>v&1)109        if not seq or S!=seq[-1]: seq.append(S)110    return sorted(set(seq),key=len)111 112 113def rand_hg(rng,n,m):114    import math as _m115    kmax=min(4,n)116    avail=sum(_m.comb(n,k) for k in range(2,kmax+1))   # never request more edges than exist117    m=min(m,avail); E=set(); guard=0118    while len(E)<m and guard<200*m+1000:119        guard+=1120        k=int(rng.integers(2,min(5,n+1))); E.add(tuple(sorted(int(x) for x in rng.choice(n,k,replace=False))))121    E=sorted(E); return E, rng.uniform(0.4,2.5,len(E))122 123if __name__=="__main__":124    rng=np.random.default_rng(5); ok=0; bad=0; tot=0125    for trial in range(150):126        n=int(rng.integers(4,11)); m=int(rng.integers(2,min(14,2**n)))127        E,w=rand_hg(rng,n,m); st={"flows":0,"ops":0}128        a=sequence(n,E,w,st); b=brute(n,E,w)129        tot+=1130        if set(a)==set(b): ok+=1131        else:132            bad+=1133            if bad<=3: print("  MISMATCH n=%d m=%d\n    flow %s\n    brute %s"%(n,m,[sorted(x) for x in a],[sorted(x) for x in b]))134    print("  certification: %d/%d instances match brute force exactly"%(ok,tot))135