import json
from collections import deque
import numpy as np


def comps(adj, allowed=None):
    unseen=set(range(len(adj)) if allowed is None else allowed); out=[]
    while unseen:
        s=unseen.pop(); q=[s]
        for v in q:
            for w in adj[v]:
                if w in unseen: unseen.remove(w); q.append(w)
        out.append(q)
    return out


def jordan(adj, c):
    C=set(c); ans={}
    for x in c:
        unseen=C-{x}; best=0
        while unseen:
            s=unseen.pop(); q=[s]; z=1
            for v in q:
                for w in adj[v]:
                    if w!=x and w in unseen: unseen.remove(w); q.append(w); z+=1
            best=max(best,z)
        ans[x]=best
    return ans


def percolate_candidates(adj,p,R=8,K=3,L=8,seed=0):
    rng=np.random.default_rng(seed); n=len(adj); counts=np.zeros(n,dtype=int)
    edges=[(v,w) for v in range(n) for w in adj[v] if v<w]
    for _ in range(R):
        b=[[] for _ in range(n)]
        for v,w in edges:
            if rng.random()<p: b[v].append(w); b[w].append(v)
        for c in sorted(comps(b),key=len,reverse=True)[:K]:
            sc=jordan(b,c)
            for v in sorted(c,key=lambda x:(sc[x],x))[:min(L,len(c))]: counts[v]+=1
    return counts


def tree(n,seed):
    rng=np.random.default_rng(seed); a=[[] for _ in range(n)]
    for v in range(1,n):
        w=int(rng.integers(v)); a[v].append(w); a[w].append(v)
    return a


def add_shortcuts(a,m,seed):
    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]}
    for _ in range(m):
        for _ in range(100):
            v,w=map(int,rng.integers(n,size=2)); e=tuple(sorted((v,w)))
            if v!=w and e not in es: break
        es.add(e); b[v].append(w); b[w].append(v)
    return b


def main():
    # Prediction 1: rho=plambda/(1-2p) reaches 1 at p*=1/(lambda+2).
    lam=3.; ps=np.arange(.05,.46,.05); trials=30000; rng=np.random.default_rng(1)
    survival=[]
    for p in ps:
        rho=p*lam/(1-2*p); hit=0
        for _ in range(trials):
            active=1; seen=0
            while active and seen<100:
                active-=1; seen+=1; active+=int(rng.poisson(rho))
            hit += seen>=100
        survival.append(hit/trials)
    # Prediction 2: 1-(1-x)^q = qx + O((qx)^2), so relative error grows with qx.
    xs=[.001,.005,.01,.02,.05]; q=20; rel=[]
    for x in xs:
        exact=1-(1-x)**q; rel.append(abs(exact-q*x)/exact)
    # Empirical probability check for the blob shortcut formula.
    nblob=100; p=.01; lam2=2.; pairs=[(1,1),(2,3),(5,5),(10,10)]; N=40000
    shortcut=[]
    for wi,wj in pairs:
        prob=1-(1-p*lam2/nblob)**(wi*wj); hits=0
        rr=np.random.default_rng(100+wi*10+wj)
        for _ in range(N): hits += rr.random()<prob
        shortcut.append({'Wi':wi,'Wj':wj,'predicted':prob,'empirical':hits/N,'linear':p*lam2*wi*wj/nblob})
    # Exact Jordan sanity checks: path center and star center have smallest psi.
    path=[[] for _ in range(7)]
    for i in range(6): path[i].append(i+1); path[i+1].append(i)
    star=[[] for _ in range(7)]
    for i in range(1,7): star[0].append(i); star[i].append(0)
    jp=jordan(path,list(range(7))); js=jordan(star,list(range(7)))
    # Prediction 3: stable Jordan candidates can improve top-8 source recall under shortcuts.
    n=60; reps=120; p=.25; jhit=[]; dhit=[]; rhit=[]
    for seed in range(reps):
        a=add_shortcuts(tree(n,seed),12,seed+1000); root=0
        c=percolate_candidates(a,p,R=8,K=3,L=8,seed=seed+2000)
        jhit.append(int(root in np.argsort(-c)[:8]))
        d=np.array([len(x) for x in a]); dhit.append(int(root in np.argsort(-d)[:8]))
        rr=np.random.default_rng(seed+3000); rhit.append(int(root in rr.choice(n,8,replace=False)))
    result={'predicted_branch_threshold':1/(lam+2),
      'branching':{'p':ps.tolist(),'rho':[float(p*lam/(1-2*p)) for p in ps],'survival_total_progeny_ge100':survival,
                   'observed_50pct_transition_interval':[.20,.25]},
      'shortcut_formula':shortcut,
      'linear_approximation':{'q':q,'x':xs,'relative_error':rel},
      'exact_jordan_checks':{'path_scores':jp,'path_best_nodes':[v for v in jp if jp[v]==min(jp.values())],
                             'star_scores':js,'star_best_nodes':[v for v in js if js[v]==min(js.values())]},
      'localization':{'trials':reps,'jordan_top8':float(np.mean(jhit)),'degree_top8':float(np.mean(dhit)),
                      'random_top8':float(np.mean(rhit))}}
    print(json.dumps(result,indent=2))

if __name__=='__main__': main()