1import json, math, random, time
  2import numpy as np
  3from coherence_sampler import LaplacianColumns, coherent_sample, uniform_sample
  4
  5
  6def weighted_graph(n, edges):
  7    a=[[] for _ in range(n)]
  8    for i,j,w in edges:
  9        a[i].append((j,w)); a[j].append((i,w))
 10    return a
 11
 12def two_cliques(k, alpha=0.0):
 13    n=2*k; e=[]
 14    for off in (0,k):
 15        for i in range(off,off+k):
 16            for j in range(i+1,off+k): e.append((i,j,1.0))
 17    if alpha:
 18        # One weak bridge keeps the block structure explicit.
 19        e.append((k-1,k,alpha))
 20    return weighted_graph(n,e)
 21
 22def complete(n):
 23    return weighted_graph(n,[(i,j,1.0) for i in range(n) for j in range(i+1,n)])
 24
 25def dense_laplacian(adj):
 26    n=len(adj); L=np.zeros((n,n))
 27    for i,nbrs in enumerate(adj):
 28        L[i,i]=sum(w for _,w in nbrs)
 29        for j,w in nbrs: L[i,j]-=w
 30    return L
 31
 32def coverage(sel,k):
 33    return len(set(i//k for i in sel))
 34
 35def mean_uniform_coverage(n,k,batch,reps=3000):
 36    vals=[]
 37    for r in range(reps): vals.append(coverage(uniform_sample(n,batch,random.Random(10000+r)),k))
 38    return float(np.mean(vals))
 39
 40def verify_sparse_math():
 41    adj=two_cliques(5,0.23); c=LaplacianColumns(adj); L=dense_laplacian(adj)
 42    max_col=max(np.max(np.abs(np.array([c.column(i).get(j,0.) for j in range(10)])-L[:,i])) for i in range(10))
 43    max_dot=max(abs(c.dot(i,j)-float(L[:,i]@L[:,j])) for i in range(10) for j in range(10))
 44    return {"max_column_error":float(max_col),"max_inner_product_error":float(max_dot)}
 45
 46def prediction_disconnected_sweep():
 47    # Prediction 1: at zero inter-block coupling, cross-block signatures have exactly zero coherence.
 48    out=[]
 49    for k in [3,5,8,12]:
 50        c=LaplacianColumns(two_cliques(k,0.0))
 51        cross=max(c.coherence(i,j) for i in range(k) for j in range(k,2*k))
 52        within=min(c.coherence(i,j) for i in range(k) for j in range(1,k))
 53        sel=coherent_sample(c,list(range(2*k)),2,random.Random(7))
 54        out.append({"k":k,"observed_cross_coherence":cross,"predicted":0.0,
 55                    "selected_blocks":coverage(sel,k),"predicted_blocks":2,
 56                    "within_coherence_min":within})
 57    return out
 58
 59def prediction_pool_sweep():
 60    # Prediction 2: with two disconnected equal blocks and an all-node candidate pool,
 61    # a batch of two selects both blocks; smaller random pools can fail only when
 62    # the pool itself misses a block.
 63    k=10; c=LaplacianColumns(two_cliques(k,0.0)); n=2*k; b=2
 64    out=[]
 65    for m in [2,3,4,6,10,15,20]:
 66        vals=[]; selected_both=0
 67        for r in range(1000):
 68            rng=random.Random(200000+m*1000+r)
 69            cand=rng.sample(range(n),m)
 70            s=coherent_sample(c,cand,b,rng)
 71            vals.append(coverage(s,k)); selected_both += coverage(s,k)==2
 72        pool_has_both=1-(2*math.comb(k,m) / math.comb(2*k,m) if m<=k else 0)
 73        out.append({"candidate_pool":m,"observed_mean_blocks":float(np.mean(vals)),
 74                    "observed_both_fraction":selected_both/1000,
 75                    "predicted_pool_has_both":float(pool_has_both)})
 76    return out
 77
 78def prediction_coupling_sweep():
 79    # Prediction 3: the cross-block column coherence is zero at alpha=0 and
 80    # increases smoothly with bridge weight; selection's block coverage falls
 81    # only once alpha is large enough to make signatures less distinct.
 82    out=[]; k=8
 83    for alpha in [0,1e-4,1e-3,1e-2,0.05,0.1,0.25,0.5,1.0,2.0]:
 84        c=LaplacianColumns(two_cliques(k,alpha))
 85        cross=max(c.coherence(i,j) for i in range(k) for j in range(k,2*k))
 86        # The maximum is attained by the two bridge endpoints.  For a k-clique
 87        # joined by one edge of weight alpha, direct substitution in c(i,S)
 88        # gives this exact prediction.
 89        d=(k-1)+alpha
 90        predicted=2*alpha*d/(d*d+(k-1)+alpha*alpha)
 91        vals=[]
 92        for r in range(200):
 93            vals.append(coverage(coherent_sample(c,list(range(2*k)),2,random.Random(3000+r)),k))
 94        uniform=[]
 95        for r in range(200):
 96            uniform.append(coverage(uniform_sample(2*k,2,random.Random(9000+r)),k))
 97        out.append({"alpha":alpha,"max_cross_coherence":float(cross),
 98                    "predicted_bridge_coherence":float(predicted),
 99                    "coherent_mean_blocks":float(np.mean(vals)),
100                    "uniform_mean_blocks":float(np.mean(uniform))})
101    return out
102
103def estimator_variance():
104    # A transparent small loss-vector experiment: empirical inclusion p is
105    # estimated from repeated candidate pools, then HT weights are evaluated.
106    n=20; k=10; c=LaplacianColumns(two_cliques(k,0.0)); B=4; M=10; R=3000
107    counts=np.zeros(n)
108    for r in range(R):
109        rng=random.Random(4000+r); cand=rng.sample(range(n),M)
110        for i in coherent_sample(c,cand,B,rng): counts[i]+=1
111    p=counts/R; loss=np.linspace(0.2,1.2,n)
112    estimates=[]; uni=[]
113    for r in range(2000):
114        rng=random.Random(8000+r); cand=rng.sample(range(n),M); s=coherent_sample(c,cand,B,rng)
115        estimates.append(sum(loss[i]/max(p[i],1/R) for i in s)/n)
116        u=uniform_sample(n,B,rng); uni.append(n*sum(loss[i] for i in u)/(B*n))
117    truth=float(np.mean(loss))
118    return {"true_mean_loss":truth,"coherent_HT_mean":float(np.mean(estimates)),
119            "coherent_HT_std":float(np.std(estimates)),"uniform_mean":float(np.mean(uni)),
120            "uniform_std":float(np.std(uni)),"min_estimated_p":float(np.min(p))}
121
122def main():
123    t=time.time()
124    result={"sparse_math":verify_sparse_math(),
125            "disconnected_sweep":prediction_disconnected_sweep(),
126            "candidate_pool_sweep":prediction_pool_sweep(),
127            "coupling_sweep":prediction_coupling_sweep(),
128            "estimator":estimator_variance(),
129            "runtime_sec":time.time()-t}
130    with open("results.json","w") as f: json.dump(result,f,indent=2)
131    print(json.dumps(result,indent=2))
132if __name__=="__main__": main()