1import json
  2from collections import deque
  3import numpy as np
  4
  5
  6def comps(adj, allowed=None):
  7    unseen=set(range(len(adj)) if allowed is None else allowed); out=[]
  8    while unseen:
  9        s=unseen.pop(); q=[s]
 10        for v in q:
 11            for w in adj[v]:
 12                if w in unseen: unseen.remove(w); q.append(w)
 13        out.append(q)
 14    return out
 15
 16
 17def jordan(adj, c):
 18    C=set(c); ans={}
 19    for x in c:
 20        unseen=C-{x}; best=0
 21        while unseen:
 22            s=unseen.pop(); q=[s]; z=1
 23            for v in q:
 24                for w in adj[v]:
 25                    if w!=x and w in unseen: unseen.remove(w); q.append(w); z+=1
 26            best=max(best,z)
 27        ans[x]=best
 28    return ans
 29
 30
 31def percolate_candidates(adj,p,R=8,K=3,L=8,seed=0):
 32    rng=np.random.default_rng(seed); n=len(adj); counts=np.zeros(n,dtype=int)
 33    edges=[(v,w) for v in range(n) for w in adj[v] if v<w]
 34    for _ in range(R):
 35        b=[[] for _ in range(n)]
 36        for v,w in edges:
 37            if rng.random()<p: b[v].append(w); b[w].append(v)
 38        for c in sorted(comps(b),key=len,reverse=True)[:K]:
 39            sc=jordan(b,c)
 40            for v in sorted(c,key=lambda x:(sc[x],x))[:min(L,len(c))]: counts[v]+=1
 41    return counts
 42
 43
 44def tree(n,seed):
 45    rng=np.random.default_rng(seed); a=[[] for _ in range(n)]
 46    for v in range(1,n):
 47        w=int(rng.integers(v)); a[v].append(w); a[w].append(v)
 48    return a
 49
 50
 51def add_shortcuts(a,m,seed):
 52    rng=np.random.default_rng(seed); b=[x[:] for x in a]; n=len(a); es={tuple(sorted((v,w))) for v in range(n) for w in a[v]}
 53    for _ in range(m):
 54        for _ in range(100):
 55            v,w=map(int,rng.integers(n,size=2)); e=tuple(sorted((v,w)))
 56            if v!=w and e not in es: break
 57        es.add(e); b[v].append(w); b[w].append(v)
 58    return b
 59
 60
 61def main():
 62    # Prediction 1: rho=plambda/(1-2p) reaches 1 at p*=1/(lambda+2).
 63    lam=3.; ps=np.arange(.05,.46,.05); trials=30000; rng=np.random.default_rng(1)
 64    survival=[]
 65    for p in ps:
 66        rho=p*lam/(1-2*p); hit=0
 67        for _ in range(trials):
 68            active=1; seen=0
 69            while active and seen<100:
 70                active-=1; seen+=1; active+=int(rng.poisson(rho))
 71            hit += seen>=100
 72        survival.append(hit/trials)
 73    # Prediction 2: 1-(1-x)^q = qx + O((qx)^2), so relative error grows with qx.
 74    xs=[.001,.005,.01,.02,.05]; q=20; rel=[]
 75    for x in xs:
 76        exact=1-(1-x)**q; rel.append(abs(exact-q*x)/exact)
 77    # Empirical probability check for the blob shortcut formula.
 78    nblob=100; p=.01; lam2=2.; pairs=[(1,1),(2,3),(5,5),(10,10)]; N=40000
 79    shortcut=[]
 80    for wi,wj in pairs:
 81        prob=1-(1-p*lam2/nblob)**(wi*wj); hits=0
 82        rr=np.random.default_rng(100+wi*10+wj)
 83        for _ in range(N): hits += rr.random()<prob
 84        shortcut.append({'Wi':wi,'Wj':wj,'predicted':prob,'empirical':hits/N,'linear':p*lam2*wi*wj/nblob})
 85    # Exact Jordan sanity checks: path center and star center have smallest psi.
 86    path=[[] for _ in range(7)]
 87    for i in range(6): path[i].append(i+1); path[i+1].append(i)
 88    star=[[] for _ in range(7)]
 89    for i in range(1,7): star[0].append(i); star[i].append(0)
 90    jp=jordan(path,list(range(7))); js=jordan(star,list(range(7)))
 91    # Prediction 3: stable Jordan candidates can improve top-8 source recall under shortcuts.
 92    n=60; reps=120; p=.25; jhit=[]; dhit=[]; rhit=[]
 93    for seed in range(reps):
 94        a=add_shortcuts(tree(n,seed),12,seed+1000); root=0
 95        c=percolate_candidates(a,p,R=8,K=3,L=8,seed=seed+2000)
 96        jhit.append(int(root in np.argsort(-c)[:8]))
 97        d=np.array([len(x) for x in a]); dhit.append(int(root in np.argsort(-d)[:8]))
 98        rr=np.random.default_rng(seed+3000); rhit.append(int(root in rr.choice(n,8,replace=False)))
 99    result={'predicted_branch_threshold':1/(lam+2),
100      'branching':{'p':ps.tolist(),'rho':[float(p*lam/(1-2*p)) for p in ps],'survival_total_progeny_ge100':survival,
101                   'observed_50pct_transition_interval':[.20,.25]},
102      'shortcut_formula':shortcut,
103      'linear_approximation':{'q':q,'x':xs,'relative_error':rel},
104      'exact_jordan_checks':{'path_scores':jp,'path_best_nodes':[v for v in jp if jp[v]==min(jp.values())],
105                             'star_scores':js,'star_best_nodes':[v for v in js if js[v]==min(js.values())]},
106      'localization':{'trials':reps,'jordan_top8':float(np.mean(jhit)),'degree_top8':float(np.mean(dhit)),
107                      'random_top8':float(np.mean(rhit))}}
108    print(json.dumps(result,indent=2))
109
110if __name__=='__main__': main()