1import numpy as np
  2from scipy.stats import linregress
  3
  4SEED = 3164
  5
  6
  7def correlated_field(n, beta=1.5, k0=1.0, rng=None):
  8    rng = np.random.default_rng() if rng is None else rng
  9    w = rng.normal(size=n)
 10    fk = np.fft.rfft(w)
 11    k = np.arange(len(fk), dtype=float)
 12    filt = (k + k0) ** (-beta / 2.0)
 13    g = np.fft.irfft(fk * filt, n=n)
 14    return (g - g.mean()) / (g.std() + 1e-12)
 15
 16
 17def make_mixer(n, alpha=1.3, beta=1.5, correlated=True, rng=None):
 18    rng = np.random.default_rng() if rng is None else rng
 19    g = correlated_field(n, beta=beta, rng=rng) if correlated else rng.normal(size=n)
 20    sigma = np.exp(g - 0.5 * np.var(g))
 21    signs = rng.choice([-1.0, 1.0], size=(n, n))
 22    signs = np.triu(signs, 1)
 23    signs = signs + signs.T
 24    d = np.abs(np.arange(n)[:, None] - np.arange(n)[None, :])
 25    M = signs * sigma[:, None] * sigma[None, :] * (1.0 + d) ** (-alpha)
 26    np.fill_diagonal(M, 0.0)
 27    raw_norm = np.linalg.svd(M, compute_uv=False)[0]
 28    return M / (raw_norm + 1e-12), raw_norm, sigma
 29
 30
 31def local_mixer(n, width=2, rng=None):
 32    rng = np.random.default_rng() if rng is None else rng
 33    d = np.abs(np.arange(n)[:, None] - np.arange(n)[None, :])
 34    M = rng.normal(size=(n, n)) * (d <= width) * (d > 0)
 35    return M / (np.linalg.svd(M, compute_uv=False)[0] + 1e-12)
 36
 37
 38def singular_gap(M):
 39    s = np.linalg.svd(M, compute_uv=False)
 40    return s[-1], 1.0 - s[1] / (s[0] + 1e-12)
 41
 42
 43def participation_entropy(v):
 44    p = np.asarray(v) ** 2
 45    p = p / (p.sum() + 1e-12)
 46    return float(-np.sum(p * np.log(p + 1e-15)))
 47
 48
 49def propagation_entropy(M, steps=8, trials=80, rng=None):
 50    rng = np.random.default_rng() if rng is None else rng
 51    vals = []
 52    # Absolute influence profile after repeated normalized residual propagation.
 53    A = np.eye(M.shape[0]) + 0.5 * M
 54    for _ in range(trials):
 55        x = rng.normal(size=M.shape[0])
 56        x /= np.linalg.norm(x)
 57        for _ in range(steps):
 58            x = A @ x
 59            x /= np.linalg.norm(x) + 1e-12
 60        vals.append(participation_entropy(x))
 61    return float(np.mean(vals))
 62
 63
 64def fit_scaling(ns, ys):
 65    fit = linregress(np.log(ns), np.log(np.maximum(ys, 1e-14)))
 66    return -fit.slope, fit.rvalue ** 2
 67
 68
 69def main():
 70    rng = np.random.default_rng(SEED)
 71    print('CORE MATH CHECK')
 72    n = 96
 73    M, raw_norm, _ = make_mixer(n, correlated=True, rng=rng)
 74    norm = np.linalg.svd(M, compute_uv=False)[0]
 75    print(f'raw_norm={raw_norm:.6f} normalized_norm={norm:.6f}')
 76    for c in (0.25, 0.5, 0.9):
 77        A = np.eye(n) + (c / (norm + 1e-12)) * M
 78        # The residual update has bounded one-step amplification <= 1+c.
 79        amp = np.linalg.svd(A, compute_uv=False)[0]
 80        print(f'c={c:.2f} observed_amp={amp:.6f} bound={1+c:.6f}')
 81    # Correlation check: edges sharing a node have correlated magnitudes only for the field model.
 82    def shared_edge_corr(correlated):
 83        a, _, s = make_mixer(128, correlated=correlated, rng=rng)
 84        # Remove deterministic distance effect by examining adjacent edges from a common node.
 85        x = np.abs(a[64, 1:64])
 86        y = np.abs(a[64, 65:128])
 87        return np.corrcoef(x, y)[0, 1]
 88    print(f'shared_endpoint_abs_corr iid={shared_edge_corr(False):.4f} correlated={shared_edge_corr(True):.4f}')
 89
 90    print('FINITE SIZE PROPAGATION ENTROPY')
 91    ns = [32, 48, 64, 96, 128, 192]
 92    result = {}
 93    for name in ('local', 'iid', 'correlated'):
 94        ents = []
 95        gaps = []
 96        for n in ns:
 97            vals_e, vals_g = [], []
 98            for rep in range(12):
 99                rr = np.random.default_rng(SEED + 1000 * rep + n)
100                if name == 'local':
101                    mm = local_mixer(n, rng=rr)
102                else:
103                    mm, _, _ = make_mixer(n, beta=1.5, correlated=(name == 'correlated'), rng=rr)
104                vals_e.append(propagation_entropy(mm, steps=8, trials=12, rng=rr))
105                vals_g.append(singular_gap(mm)[0])
106            ents.append(np.mean(vals_e)); gaps.append(np.mean(vals_g))
107        # Compare the two proposed entropy forms using residual sum of squares.
108        x = np.log(np.asarray(ns, float))
109        X1 = np.column_stack([x, np.ones_like(x)])
110        X2 = np.column_stack([x*x, x, np.ones_like(x)])
111        rss_log = np.sum((np.asarray(ents) - X1 @ np.linalg.lstsq(X1, ents, rcond=None)[0]) ** 2)
112        rss_log2 = np.sum((np.asarray(ents) - X2 @ np.linalg.lstsq(X2, ents, rcond=None)[0]) ** 2)
113        z, r2 = fit_scaling(ns, gaps)
114        result[name] = (ents, gaps, z, r2, rss_log, rss_log2)
115        print(f'{name:10s} entropy=' + ','.join(f'{v:.4f}' for v in ents))
116        print(f'{name:10s} gap=' + ','.join(f'{v:.6f}' for v in gaps) + f' z={z:.3f} R2={r2:.3f} RSS_log={rss_log:.6g} RSS_log2={rss_log2:.6g}')
117
118    print('STABILITY WITHOUT NORMALIZATION')
119    mm, raw_norm, _ = make_mixer(96, correlated=True, rng=np.random.default_rng(SEED))
120    raw = mm * raw_norm
121    x = np.random.default_rng(SEED + 9).normal(size=96); x /= np.linalg.norm(x)
122    for c in (0.5, 1.0):
123        y = x.copy()
124        for _ in range(30):
125            y = y + c * raw @ y
126        print(f'raw_update_c={c:.1f} final_norm={np.linalg.norm(y):.3e}')
127
128
129if __name__ == '__main__':
130    main()