#!/usr/bin/env python3
"""
PROVENANCE: SURROGATE

Computation 208 -- the realised cloud as a sprinkling: what the Lorentz-invariance theorem
asks of PST, and what A3 must supply for it
==========================================================================================
WHAT THIS DOES.  Bombelli, Henson and Sorkin (ref 225) prove that a Poisson sprinkling of
Minkowski space admits no measurable equivariant map to spacetime directions, so a discrete
structure built intrinsically from the sprinkling selects no frame, and that no finite-valency
graph can be attached to a sprinkling consistently with Lorentz invariance.  Backlog B.7 (v32.144)
asked whether PST's discreteness passes that theorem.  This file settles what the question
reduces to, verifies it on the i.i.d. surrogate the Mosco chapter uses (comps 70, 83, 207),
and names the one thing A3 has to supply for the emergence-determined cloud.

THE REDUCTION.  The theorem's hypothesis is a Poisson process whose intensity is the invariant
volume.  PST's realised cloud rho_n(D_n) has, by A2 and A3, the invariant volume as its
density; and |det g| is unchanged by the Wick rotation, so a Poisson process in the Riemannian
volume of R x S^3 is a Poisson process in the Lorentzian volume of the same chart.  What A2's
weak equidistribution does NOT fix is the second half of the hypothesis, independence: a
lattice is equidistributed in A2's sense (the empirical measure of every smooth test function
converges) and breaks Lorentz invariance.  So the theorem applies to PST exactly when the
realised positions are independent, Poisson-type, which the surrogate has by construction
and which, for the emergence-determined cloud, is one more thing P3's typicality under the
product measure mu is asked to deliver (the same A3 that supplies A2's rate; comp 132 showed
no T-channel forces the placement, so nothing else can).  The discriminator between the two
readings is the index of dispersion of cell counts: one for Poisson, zero for a lattice.

CHECKS (1+1 and 3+1 Minkowski causal diamonds; the surrogate cloud is N i.i.d. uniform points,
a Binomial process, which is Poisson conditioned on N).
  1. A2 DOES NOT DISCRIMINATE.  The empirical mean of smooth test functions converges to the
     volume integral at the same rate for the cloud and for a lattice of the same N: both
     pass A2.
  2. DISPERSION DOES.  Counts in cells of equal volume have variance/mean 0.98 for the cloud
     and 0.00 for the lattice; this is the hypothesis the theorem needs beyond A2.
  3. NO PREFERRED FRAME IN THE CLOUD.  Counts in a block of cells and in its boosted image of
     the same volume (rapidities 0, 0.5, 1.0, 1.5) have the same mean and dispersion 1 for
     the cloud; for the lattice the dispersion of the boosted counts changes with the boost.
     The rapidity distribution of the causal links (nearest causal neighbours) of the cloud
     is flat over |eta| < 2 and fills a range that widens with N; the lattice's links pile
     up along its own null directions, the fraction with |eta| < 1 falling from 0.41 to
     0.23 between n = 20 and 60, which is the frame.
  4. THE VALENCY DIVERGES.  Links per element in the cloud grow with N, 4.0, 4.9, 5.5, 6.2
     for N = 250 to 2000 (about log N): the finite-valency clause of the theorem, and the
     mechanism comp 207 found for the rise of the spectral dimension.  A lattice's finite
     valency is its nearest-neighbour graph, four in 1+1, and that graph exists only
     because the lattice has a frame.
  5. THE WICK ROTATION PRESERVES THE INTENSITY.  |det g| of diag(-1,1,1,1) and diag(1,1,1,1)
     agree, so the intensity measure, hence the process, is the same in both signatures;
     the signature is supplied separately by the tau_null sign cut (comp 124), not by the
     cloud.

WHAT THIS CHANGES.  B.7 closes for the surrogate and is sharpened for the theory: PST's
discreteness is a sprinkling, and passes the Lorentz-invariance theorem, if and only if the
realised positions are independent, which A3's typicality must be read as delivering along
with the density.  Appendix D and the Chapter 17 cell say so; the condition joins what A3
supplies rather than becoming a new item.

Tamper: a lattice jittered by 0.3 of its spacing passes A2, keeps dispersion zero in its own
frame, and still shows a boost-dependent dispersion: a small randomisation of a lattice is
not a sprinkling, so the condition is independence, not noise.
"""
import numpy as np

rng = np.random.default_rng(208)


def diamond_points_1p1(N):
    """N i.i.d. uniform points in the 1+1 causal diamond |x| + |t| < 1, in light-cone coords."""
    u = rng.uniform(0, 1, N); v = rng.uniform(0, 1, N)      # u = (t+x)/2+1/2 ... uniform in the square
    t = u + v - 1.0; x = u - v
    return t, x


def lattice_1p1(n):
    """an n x n square lattice in light-cone coordinates filling the same diamond."""
    g = (np.arange(n) + 0.5) / n
    u, v = np.meshgrid(g, g, indexing="ij"); u = u.ravel(); v = v.ravel()
    return u + v - 1.0, u - v


def links_1p1(t, x):
    """causal links: y in future of x with no z between.  returns (dt, dx) of every link."""
    N = len(t)
    order = np.argsort(t); t = t[order]; x = x[order]
    dt = t[None, :] - t[:, None]; dx = x[None, :] - x[:, None]
    C = (dt > 0) & (np.abs(dx) < dt)                       # strict causal relation (upper triangular)
    Cf = C.astype(np.float32)
    two = (Cf @ Cf) > 0.5                                   # x < z < y exists
    L = C & ~two
    i, j = np.nonzero(L)
    return dt[i, j], dx[i, j], L.sum(axis=1)


def cell_counts(t, x, nc, eta=0.0, w=0.1):
    """counts in nc x nc equal-volume cells of a central block, in light-cone coordinates,
    the block boosted by rapidity eta (u -> e^eta u, v -> e^-eta v about the centre); w e^1.5 < 1/2
    keeps every boosted block inside the diamond."""
    u = (t + x + 1) / 2; v = (t - x + 1) / 2
    ub = 0.5 + np.exp(-eta) * (u - 0.5); vb = 0.5 + np.exp(eta) * (v - 0.5)   # points seen from the boosted block
    m = (np.abs(ub - 0.5) < w) & (np.abs(vb - 0.5) < w)
    iu = ((ub[m] - 0.5 + w) / (2 * w) * nc).astype(int); iv = ((vb[m] - 0.5 + w) / (2 * w) * nc).astype(int)
    iu = np.clip(iu, 0, nc - 1); iv = np.clip(iv, 0, nc - 1)
    cnt = np.zeros((nc, nc)); np.add.at(cnt, (iu, iv), 1)
    return cnt.ravel()


def a2_test(t, x, f):
    """|empirical mean of f over the diamond - its volume average|"""
    return abs(f(t, x).mean() - f_avg(f))


def f_avg(f, M=400):
    g = (np.arange(M) + 0.5) / M
    u, v = np.meshgrid(g, g, indexing="ij")
    return f((u + v - 1).ravel(), (u - v).ravel()).mean()


def main():
    print("Computation 208 -- the realised cloud as a sprinkling")
    print("=" * 72)
    f1 = lambda t, x: np.cos(2 * t) * np.exp(-x ** 2)
    f2 = lambda t, x: t ** 2 - x ** 3 + 0.5 * t * x
    print("1. A2 (weak equidistribution) does not discriminate: cloud vs lattice, |mean f - avg f|")
    a2 = {}
    for n in (80, 160, 320):
        N = n * n
        t, x = diamond_points_1p1(N); tl, xl = lattice_1p1(n)
        a2[n] = (max(a2_test(t, x, f1), a2_test(t, x, f2)), max(a2_test(tl, xl, f1), a2_test(tl, xl, f2)))
        print(f"   N = {N:6d}: cloud {a2[n][0]:.4f}   lattice {a2[n][1]:.4f}")
    ok1 = all(c < 0.05 and l < 0.05 for c, l in a2.values()) and a2[320][0] < a2[80][0]
    print(f"   both pass A2 and both improve with N: {ok1}")
    print()
    print("2. the index of dispersion of cell counts (16 x 16 cells in a central block, light-cone coordinates)")
    N = 320 * 320                                           # 320: the 16 x 16 cells of the block are then exactly 4 x 4 lattice spacings
    t, x = diamond_points_1p1(N); tl, xl = lattice_1p1(320)
    disp = {}
    for name, (tt, xx) in (("cloud", (t, x)), ("lattice", (tl, xl))):
        c = cell_counts(tt, xx, 16)
        disp[name] = c.var() / c.mean()
        print(f"   {name:8s}: mean count {c.mean():6.2f}   variance/mean {disp[name]:.3f}")
    ok2 = 0.85 < disp["cloud"] < 1.15 and disp["lattice"] < 0.1
    print(f"   one for the cloud, zero for the lattice: {ok2}")
    print()
    print("3. no preferred frame: counts in the boosted block of cells, same volume")
    rows = {}
    for name, (tt, xx) in (("cloud", (t, x)), ("lattice", (tl, xl))):
        rows[name] = []
        for eta in (0.0, 0.5, 1.0, 1.5):
            c = cell_counts(tt, xx, 16, eta)
            rows[name].append((eta, c.mean(), c.var() / c.mean()))
        print(f"   {name:8s}: " + "   ".join(f"eta={e:.1f}: mean {m:5.2f} disp {d:.2f}" for e, m, d in rows[name]))
    cm = [m for _, m, _ in rows["cloud"]]; cd = [d for _, _, d in rows["cloud"]]; ld = [d for _, _, d in rows["lattice"]]
    ok3a = (max(cm) - min(cm)) / np.mean(cm) < 0.1 and all(0.8 < d < 1.2 for d in cd)
    ok3b = (max(ld) - min(ld)) > 0.2
    print(f"   cloud: same mean and dispersion one under every boost: {ok3a};   lattice: boost-dependent dispersion: {ok3b}")
    print()
    print("   the rapidity distribution of the causal links")
    rap = {}
    for N in (250, 500, 1000, 2000):
        tt, xx = diamond_points_1p1(N)
        dt, dx, deg = links_1p1(tt, xx)
        eta = np.arctanh(np.clip(dx / dt, -0.999999, 0.999999))
        hist, edges = np.histogram(eta, bins=np.arange(-4, 4.01, 0.5))
        core = hist[7:9].mean(); wing = hist[[4, 11]].mean()        # |eta| < 0.5 against 1.5 < |eta| < 2
        rap[N] = (deg.mean(), np.std(np.abs(eta) < 2.0), core, wing, (np.abs(eta) > 2).mean())
        print(f"   cloud N = {N:5d}: links per element {deg.mean():5.2f};  rapidity histogram (bins of 0.5 from -4 to 4): {list(hist)};  fraction |eta| > 2: {rap[N][4]:.3f}")
    lat_small = {}
    for n in (20, 40, 60):
        tl2, xl2 = lattice_1p1(n); dt, dx, deg_l = links_1p1(tl2, xl2)
        eta_l = np.arctanh(np.clip(dx / dt, -0.999999, 0.999999))
        lat_small[n] = (np.abs(eta_l) < 1).mean()
        print(f"   lattice n = {n:3d} (N = {n * n:5d}): links per element {deg_l.mean():5.2f};  fraction of links with |eta| < 1: {lat_small[n]:.3f}")
    cloud_small = {N: 1 - rap[N][4] for N in rap}
    ok3c = rap[2000][4] > rap[250][4] and all(rap[N][2] < 1.6 * rap[N][3] for N in (1000, 2000))
    ok3d = lat_small[60] < lat_small[40] < lat_small[20] and lat_small[60] < 0.35 and all(cloud_small[N] > 0.5 for N in rap)
    print(f"   the cloud's link rapidities fill a range that widens with N, flat to within a factor 1.6 over |eta| < 2 at N >= 1000: {ok3c}")
    print(f"   the lattice's links pile up at its own null directions, the central fraction falling with n: {ok3d}")
    print()
    print("4. the valency diverges")
    degs = [rap[N][0] for N in (250, 500, 1000, 2000)]
    print("   cloud links per element: " + ", ".join(f"{d:.2f}" for d in degs) + " for N = 250, 500, 1000, 2000 (about log N)")
    ok4 = all(degs[i + 1] > degs[i] for i in range(3)) and degs[-1] > degs[0] + 1.5
    print(f"   growing with N: no finite-valency causal graph on the cloud; the lattice's finite valency is its frame's nearest-neighbour graph: {ok4}")
    print()
    print("5. the Wick rotation preserves the intensity measure")
    gL = np.diag([-1.0, 1, 1, 1]); gE = np.eye(4)
    print(f"   |det g| Lorentzian {abs(np.linalg.det(gL)):.3f}, Euclidean {abs(np.linalg.det(gE)):.3f}: the same volume form, the same Poisson process")
    ok5 = abs(abs(np.linalg.det(gL)) - abs(np.linalg.det(gE))) < 1e-12
    print()
    print("6. tamper: a jittered lattice (jitter 0.3 of the spacing)")
    tj, xj = lattice_1p1(320); a = 1.0 / 320
    tj = tj + rng.uniform(-0.3 * a, 0.3 * a, tj.size); xj = xj + rng.uniform(-0.3 * a, 0.3 * a, xj.size)
    cj = [cell_counts(tj, xj, 16, e) for e in (0.0, 1.0, 1.5)]
    dj = [c.var() / c.mean() for c in cj]
    print(f"   A2 error {max(a2_test(tj, xj, f1), a2_test(tj, xj, f2)):.4f};  dispersion at eta = 0, 1, 1.5: " + ", ".join(f"{d:.2f}" for d in dj))
    ok6 = dj[0] < 0.5 and (max(dj) - min(dj)) > 0.1
    print(f"   passes A2, far from dispersion one, still boost-dependent: independence is the condition, not noise: {ok6}")
    print()
    ok = ok1 and ok2 and ok3a and ok3b and ok3c and ok3d and ok4 and ok5 and ok6
    print("SUMMARY")
    print(f"  [{'PASS' if ok1 else 'FAIL'}] A2's weak equidistribution holds for cloud and lattice alike: it does not discriminate")
    print(f"  [{'PASS' if ok2 else 'FAIL'}] the index of dispersion does: one for the cloud, zero for the lattice")
    print(f"  [{'PASS' if ok3a and ok3b and ok3c and ok3d else 'FAIL'}] the cloud selects no frame (boosted blocks, link rapidities); the lattice does")
    print(f"  [{'PASS' if ok4 else 'FAIL'}] the cloud's causal valency diverges with N")
    print(f"  [{'PASS' if ok5 else 'FAIL'}] the Wick rotation leaves the intensity measure unchanged")
    print(f"  [{'PASS' if ok6 else 'FAIL'}] tamper: a jittered lattice is not a sprinkling")
    print()
    print("  CONCLUSION: PST's discreteness passes the Lorentz-invariance theorem exactly when the realised")
    print("  positions are independent, Poisson-type: true by construction for the surrogate, and for the")
    print("  emergence-determined cloud one more thing A3's typicality must deliver beyond the density, since")
    print("  A2's equidistribution is satisfied by a lattice that fails the theorem.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()