1import json, time
 2import numpy as np
 3from scipy.optimize import minimize
 4
 5SEED = 7
 6RNG = np.random.default_rng(SEED)
 7
 8def entropy(A, pi):
 9    x = np.clip(A, 1e-9, 1 - 1e-9)
10    return float(np.sum(pi[:, None] * pi[None, :] * (x*np.log(x) + (1-x)*np.log(1-x))))
11
12def edge_density(A, pi):
13    return float(np.einsum('i,j,ij->', pi, pi, A))
14
15def triangle_density(A, pi):
16    return float(np.einsum('i,j,k,ij,jk,ki->', pi, pi, pi, A, A, A))
17
18def wedge_density(A, pi):
19    return float(np.einsum('i,j,k,ij,jk->', pi, pi, pi, A, A))
20
21def optimize_entropy(pi, edge_target, triangle_target=None, seed=0):
22    m = len(pi); rng = np.random.default_rng(seed)
23    x0 = np.clip(edge_target + rng.normal(0, .04, (m, m)), .03, .97).ravel()
24    constraints = [{'type':'eq', 'fun': lambda x: edge_density(x.reshape(m,m), pi)-edge_target}]
25    if triangle_target is not None:
26        constraints.append({'type':'eq', 'fun': lambda x: triangle_density(x.reshape(m,m), pi)-triangle_target})
27    result = minimize(lambda x: entropy(x.reshape(m,m), pi), x0, method='SLSQP',
28                      bounds=[(1e-5, .99999)]*(m*m), constraints=constraints,
29                      options={'ftol':1e-12, 'maxiter':1000})
30    return result.x.reshape(m,m), result
31
32def exact_vs_sampled_motif():
33    pi = np.array([.2, .5, .3])
34    A = np.array([[.15,.4,.7],[.25,.6,.35],[.8,.3,.55]])
35    exact = {'edge': edge_density(A,pi), 'wedge': wedge_density(A,pi), 'triangle': triangle_density(A,pi)}
36    out = []
37    for n in [100, 300, 1000, 3000, 10000, 30000]:
38        errs = []
39        for rep in range(40):
40            z = RNG.choice(3, size=(n,3), p=pi)
41            errs.append(np.mean(A[z[:,0],z[:,1]]*A[z[:,1],z[:,2]]*A[z[:,2],z[:,0]])-exact['triangle'])
42        out.append({'samples':n, 'rmse':float(np.sqrt(np.mean(np.square(errs))))})
43    return exact, out
44
45def entropy_sweep():
46    # Prediction 1: with only edge density, maximum entropy is homogeneous A=p.
47    # Prediction 2: entropy decreases monotonically as a non-extremal target triangle
48    # is forced away from its homogeneous value, while edge density remains fixed.
49    pi = np.array([.25,.35,.40]); p=.30
50    homogeneous_triangle=p**3
51    rows=[]
52    for tri in [homogeneous_triangle, .030, .040, .050, .060, .075, .090]:
53        A,res=optimize_entropy(pi,p,tri,seed=11)
54        rows.append({'target_triangle':tri, 'achieved_triangle':triangle_density(A,pi),
55                     'edge_error':abs(edge_density(A,pi)-p), 'entropy':entropy(A,pi),
56                     'heterogeneity':float(np.sqrt(np.sum(pi[:,None]*pi[None,:]*(A-p)**2))),
57                     'success':bool(res.success)})
58    one=[]
59    for m in [2,3,4,5]:
60        pi=np.ones(m)/m
61        A,res=optimize_entropy(pi,p,None,seed=m)
62        one.append({'m':m,'max_abs_deviation':float(np.max(np.abs(A-p))),
63                    'weighted_rmse':float(np.sqrt(np.sum(pi[:,None]*pi[None,:]*(A-p)**2))),
64                    'success':bool(res.success)})
65    return {'homogeneous_triangle':homogeneous_triangle,'triangle_sweep':rows,'edge_only_block_sweep':one}
66
67def propagation_benchmark():
68    # Compare dense relation tensor application with the proposed block factorization.
69    rng=np.random.default_rng(19); n=1800; m=8; r=3; d=32
70    P=rng.dirichlet(np.ones(m), size=n).astype(np.float64)
71    V=rng.normal(size=(n,d)).astype(np.float64)
72    A=rng.uniform(.1,.9,size=(r,m,m)).astype(np.float64)
73    # Dense kernel is n*n*r storage and application; use n=1800 to stay modest.
74    W=np.einsum('vi,kij,wj->kvw',P,A,P)
75    dense_bytes=W.nbytes
76    t0=time.perf_counter(); dense=np.einsum('kvw,wd->vd',W,V); dense_time=time.perf_counter()-t0
77    t0=time.perf_counter()
78    S=np.einsum('wj,wd->jd',P,V)
79    block=np.einsum('vi,kij,jd->vd',P,A,S)
80    block_time=time.perf_counter()-t0
81    rel_err=float(np.max(np.abs(dense-block))/(np.max(np.abs(dense))+1e-12))
82    block_bytes=P.nbytes+A.nbytes
83    return {'n':n,'m':m,'relations':r,'hidden':d,'dense_bytes':dense_bytes,
84            'block_bytes':block_bytes,'memory_reduction':dense_bytes/block_bytes,
85            'dense_seconds':dense_time,'block_seconds':block_time,'max_relative_error':rel_err}
86
87def main():
88    exact, sampling=exact_vs_sampled_motif()
89    result={'seed':SEED,'motif_exact':exact,'motif_sampling_rmse':sampling,
90            'entropy_sweep':entropy_sweep(),'propagation':propagation_benchmark()}
91    with open('results.json','w') as f: json.dump(result,f,indent=2)
92    print(json.dumps(result,indent=2))
93
94if __name__=='__main__': main()