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

Computation 211 -- maximum entropy supplies the independence: the realised cloud is a
sprinkling at the grain, and A3's rate touches only the scales above the constraint
=====================================================================================
WHAT THIS DOES.  comp 208 reduced the Lorentz invariance of PST's discreteness to one
condition beyond A2's density: that the realised positions be independent, Poisson-type
(index of dispersion one), which a lattice fails while passing A2; it named that condition
as one more thing P3's typicality must deliver.  comp 209 showed the explicit placement
cannot test it.  This file tests the clause itself: P3's maximum-entropy sublimation, read
as a law on the realised cloud rather than as a placement formula.  Two facts, one classical
and one measured, settle what that clause delivers.

  1. THE THEOREM.  Among all laws of N labelled positions with the same one-point marginals,
     the product law has the largest entropy: H(X_1, ..., X_N) <= sum_i H(X_i) with equality
     iff the positions are independent (subadditivity; the deficit is the relative entropy
     to the product law).  So "maximum entropy given the density" IS independence.  Part 1
     verifies this exactly on K = 8 cells and N = 16 points over the exchangeable family
     p(n) ~ multinomial(n) exp(-J sum_k n_k^2): the labelled-configuration entropy is
     N log K = 33.27 at J = 0 and strictly below it on both sides, the dispersion of the
     cell counts at J = 0 is the binomial 1 - 1/K = 0.875, and a lattice-like law (J > 0)
     and a clustered one (J < 0) pay for their dispersion in entropy.  The entropy that is
     maximised is that of the labelled cloud; the entropy of the occupancy vector alone is
     not maximal at J = 0 (it is maximal for a uniform law on compositions), which is why
     the theorem is about positions and not about counts.

  2. THE CONSTRAINED CLOUD.  A3 asks the realised cloud to be equidistributed at a rate:
     the empirical mean of every smooth test function matches its volume average to within
     a margin.  The maximum-entropy law under that constraint is the product law conditioned
     on the constraint set, i.e. the uniform law on the clouds that satisfy it.  Part 2
     samples it: N = 4000 points in the 1+1 causal diamond (light-cone coordinates u, v
     uniform on the unit square, as in comp 208), the constraint |mean f| <= eps for the
     eight lowest Fourier modes cos(pi j u) cos(pi k v), j, k <= 2, with eps = 0.5 and 0.1
     times the i.i.d. scale 1/sqrt(N) (the i.i.d. cloud violates the tighter one), sampled
     by a greedy descent into the set followed by single-point Metropolis moves rejected
     outside it, thirty independent clouds per eps.  Measured over the clouds: the index of
     dispersion of cell counts at cell sizes 0.01 to 0.5, and in a central block of small
     cells boosted by rapidities 0 to 1.5.
     RESULT: the dispersion is one to within 0.01 at every cell size up to 0.05 for every
     eps and 0.92 to 1.02 in the boosted small cells at every rapidity; it falls only as the
     cells approach the constraint's half-wavelength 0.25 (cell 0.10: 1.01, 0.94, 0.89; cell
     0.25: 0.95, 0.67, 0.54; cell 0.50: 0.75, 0.36, 0.24 for i.i.d. and the two eps, the
     i.i.d. values being the binomial 1 - 1/cells).  The constraint makes the cloud
     hyperuniform at the Mosco scale and leaves it Poisson at the grain: A3 and the
     sprinkling hypothesis are not in tension, they live at different scales, and maximum
     entropy delivers both at once.

  3. THE LINKS.  On one constrained cloud of N = 1000 (eps = 0.1/sqrt N) the causal links
     have the sprinkling's signature from comp 208: 5.6 links per element, a rapidity
     distribution flat to within a factor 1.2 between |eta| < 0.5 and 1.5 < |eta| < 2, and
     no pile-up at any null direction.

WHAT THIS CHANGES.  B.7 closes.  The independence comp 208 named is not a further thing
A3 must be read as supplying by hand: it is what maximum entropy means once the clause is
read as a law on the realised cloud (the theorem), and the rate A3 states does not disturb
it at the grain (the measurement).  PST's discreteness is a sprinkling in the sense the
Bombelli--Henson--Sorkin theorem needs, with the density fixed by A2 and the
independence by the maximum-entropy clause, both on P3's footing.  What stays open is
B.6's question of what the fibres of rho are, which this file does not touch: it places
labelled sites, and a point of M that is a fibre of several sites is a thinning of this
cloud, which preserves the Poisson property.

TAMPERS.  (a) A jittered lattice (jitter 0.3 of the spacing) satisfies the same constraints
six times inside the tight eps and has dispersion 0.35 and 0.14 at cells 0.02 and 0.05, and
0.30, 0.38, 1.08 in the boosted small cells: A2 plus the rate does not give independence and
does not remove the frame; the clause does.  (b) Conditioning with eps a tenth of the i.i.d.
scale is a strong conditioning (the greedy phase moves about 150 of 4000 points) and still
leaves the small-cell dispersion at one: the set the cloud is conditioned into is large and
its small-scale statistics are the product law's.

Runs in about half a minute with numpy.
"""
import math
import numpy as np

rng = np.random.default_rng(211)


# ---------------------------------------------------------------- part 1: the theorem on cells
def compositions(N, K):
    out = []
    def rec(prefix, rem, slots):
        if slots == 1:
            out.append(prefix + [rem]); return
        for k in range(rem + 1):
            rec(prefix + [k], rem - k, slots - 1)
    rec([], N, K)
    return np.array(out, dtype=float)


def part1(K=8, N=16):
    C = compositions(N, K)
    logmult = np.array([math.lgamma(N + 1) - sum(math.lgamma(n + 1) for n in c) for c in C])
    S2 = (C ** 2).sum(1)
    rows = []
    for J in (-0.1, -0.03, 0.0, 0.03, 0.1, 0.3):
        w = logmult - N * math.log(K) - J * S2
        w -= w.max(); p = np.exp(w); p /= p.sum()
        H_counts = -(p * np.log(p)).sum()
        H_full = H_counts + (p * logmult).sum()            # entropy of the labelled cloud
        mean = (p[:, None] * C).sum(0); var = (p[:, None] * (C - mean) ** 2).sum(0)
        rows.append((J, H_full, H_counts, (var / mean).mean()))
    return rows, N * math.log(K), 1 - 1 / K


# ---------------------------------------------------------------- part 2: the constrained cloud
FAM = [(j, k) for j in range(3) for k in range(3) if (j, k) != (0, 0)]


def F(u, v):
    return np.stack([np.cos(np.pi * j * u) * np.cos(np.pi * k * v) for j, k in FAM], axis=1)


def F1(u, v):
    """F for one point as a flat vector: the same products in the same order (bit-identical
    to F(np.array([u]), np.array([v]))[0], checked on 200,000 points), seven times faster."""
    out = np.empty(len(FAM))
    for i, (j, k) in enumerate(FAM):
        out[i] = np.cos(np.pi * j * u) * np.cos(np.pi * k * v)
    return out


def constrained_cloud(N, eps, sweeps):
    """the uniform law on {clouds : |mean f| <= eps for every f in FAM}; eps = None is i.i.d."""
    u = rng.uniform(0, 1, N); v = rng.uniform(0, 1, N)
    if eps is None:
        return u, v, 0
    S = F(u, v).sum(0); lim = eps * N
    viol = lambda S: np.maximum(np.abs(S) - lim, 0).sum()
    moved = 0
    while viol(S) > 0:                                      # greedy descent into the constraint set
        i = rng.integers(N); un, vn = rng.uniform(), rng.uniform()
        Sn = S - F1(u[i], v[i]) + F1(un, vn)
        if viol(Sn) <= viol(S):
            u[i], v[i], S = un, vn, Sn; moved += 1
    for _ in range(sweeps):                                 # Metropolis, uniform within the set
        idx = rng.integers(0, N, N); UN = rng.uniform(0, 1, N); VN = rng.uniform(0, 1, N)
        for i, un, vn in zip(idx, UN, VN):
            Sn = S - F1(u[i], v[i]) + F1(un, vn)
            if np.all(np.abs(Sn) <= lim):
                u[i], v[i], S = un, vn, Sn
    return u, v, moved


def jittered_lattice(n, jitter):
    g = (np.arange(n) + 0.5) / n
    u, v = np.meshgrid(g, g, indexing="ij"); u = u.ravel(); v = v.ravel()
    u = (u + rng.uniform(-jitter / n, jitter / n, u.size)) % 1.0
    v = (v + rng.uniform(-jitter / n, jitter / n, v.size)) % 1.0
    return u, v


def disp_cells(u, v, s):
    nc = int(round(1 / s))
    iu = np.minimum((u * nc).astype(int), nc - 1); iv = np.minimum((v * nc).astype(int), nc - 1)
    cnt = np.zeros((nc, nc)); np.add.at(cnt, (iu, iv), 1)
    return cnt.var() / cnt.mean()


def disp_block(u, v, nc, eta, w):
    """counts in nc x nc cells of a central block of half-width w, the block boosted by rapidity eta"""
    ub = 0.5 + np.exp(-eta) * (u - 0.5); vb = 0.5 + np.exp(eta) * (v - 0.5)
    m = (np.abs(ub - 0.5) < w) & (np.abs(vb - 0.5) < w)
    iu = np.clip(((ub[m] - 0.5 + w) / (2 * w) * nc).astype(int), 0, nc - 1)
    iv = np.clip(((vb[m] - 0.5 + w) / (2 * w) * nc).astype(int), 0, nc - 1)
    cnt = np.zeros((nc, nc)); np.add.at(cnt, (iu, iv), 1)
    return cnt.var() / cnt.mean()


def links_1p1(t, x):
    order = np.argsort(t); t = t[order]; x = x[order]
    dt = t[None, :] - t[:, None]; dx = x[None, :] - x[:, None]
    Cr = (dt > 0) & (np.abs(dx) < dt)
    two = (Cr.astype(np.float32) @ Cr.astype(np.float32)) > 0.5
    L = Cr & ~two
    i, j = np.nonzero(L)
    return dt[i, j], dx[i, j], L.sum(axis=1)


SIZES = (0.01, 0.02, 0.05, 0.1, 0.25, 0.5)
ETAS = (0.0, 0.5, 1.0, 1.5)


def part2(N=4000, R=30, sweeps=6):
    scale = 1 / math.sqrt(N)
    table = {}
    for eps in (None, 0.5 * scale, 0.1 * scale):
        d = np.zeros((R, len(SIZES))); b = np.zeros((R, len(ETAS))); moved = []
        for r in range(R):
            u, v, mv = constrained_cloud(N, eps, sweeps); moved.append(mv)
            d[r] = [disp_cells(u, v, s) for s in SIZES]
            b[r] = [disp_block(u, v, 8, e, 0.08) for e in ETAS]      # cells of 0.02 in a block of half-width 0.08
        table[eps] = (d.mean(0), d.std(0) / math.sqrt(R), b.mean(0), b.std(0) / math.sqrt(R), np.mean(moved))
    return table, scale


def main():
    print("Computation 211 -- maximum entropy supplies the independence")
    print("=" * 72)
    rows, Hmax, binom = part1()
    print("1. the theorem on K = 8 cells, N = 16 labelled points, p(n) ~ multinomial(n) exp(-J sum n_k^2)")
    for J, Hf, Hc, disp in rows:
        print(f"   J = {J:+.2f}: labelled-cloud entropy {Hf:7.4f}   (counts alone {Hc:6.3f})   dispersion of cell counts {disp:.3f}")
    print(f"   N log K = {Hmax:.4f}: the product law's entropy, reached at J = 0 and nowhere else; binomial dispersion 1 - 1/K = {binom:.3f}")
    H0 = [Hf for J, Hf, _, _ in rows if J == 0][0]
    ok1 = abs(H0 - Hmax) < 1e-9 and all(Hf < H0 - 1e-6 for J, Hf, _, _ in rows if J != 0) and abs([d for J, _, _, d in rows if J == 0][0] - binom) < 1e-9
    print(f"   entropy maximal exactly at independence, dispersion binomial there: {ok1}")
    print()
    table, scale = part2()
    print(f"2. the constrained cloud: N = 4000, eight lowest modes constrained to |mean f| <= eps, 30 clouds each (i.i.d. scale 1/sqrt N = {scale:.4f})")
    print("   cell size:        " + "  ".join(f"{s:>6.2f}" for s in SIZES) + "   |  boosted 0.02-cells, eta = " + ", ".join(f"{e:.1f}" for e in ETAS))
    for eps, (dm, de, bm, be, mv) in table.items():
        tag = "i.i.d.     " if eps is None else f"eps={eps:.4f}"
        print(f"   {tag}: " + "  ".join(f"{m:6.3f}" for m in dm) + "   |  " + "  ".join(f"{m:.2f}" for m in bm) + f"   (greedy moves {mv:.0f}; s.e. of each cell-size entry <= {de.max():.3f})")
    keys = list(table)
    small = [i for i, s in enumerate(SIZES) if s <= 0.05]
    ok2a = all(abs(table[k][0][i] - 1.0) < 0.05 for k in keys for i in small)
    ok2b = table[keys[2]][0][-1] < 0.4 and table[keys[1]][0][-1] < table[keys[0]][0][-1] - 0.2 and abs(table[keys[0]][0][-1] - 0.75) < 0.1
    ok2c = all(abs(table[k][2][j] - 1.0) < 0.12 for k in keys for j in range(len(ETAS)))
    print(f"   dispersion one (within 0.05) at every cell size <= 0.05 for every eps: {ok2a}")
    print(f"   suppressed only toward the constraint scale (cell 0.5: i.i.d. {table[keys[0]][0][-1]:.2f}, the binomial 1 - 1/4; eps {table[keys[1]][0][-1]:.2f} and {table[keys[2]][0][-1]:.2f}): {ok2b}")
    print(f"   boosted small cells at one (within 0.12) under every rapidity for every eps: {ok2c}")
    print()
    print("3. the causal links of one constrained cloud (N = 1000, eps = 0.1/sqrt N) against a lattice")
    u, v, _ = constrained_cloud(1000, 0.1 / math.sqrt(1000), 6)
    t = u + v - 1; x = u - v
    dt, dx, deg = links_1p1(t, x)
    eta = np.arctanh(np.clip(dx / dt, -0.999999, 0.999999))
    hist, _ = np.histogram(eta, bins=np.arange(-4, 4.01, 0.5)); hist = [int(h) for h in hist]
    core = (hist[7] + hist[8]) / 2; wing = (hist[4] + hist[11]) / 2
    ul, vl = jittered_lattice(32, 0.0); dtl, dxl, degl = links_1p1(ul + vl - 1, ul - vl)
    etal = np.arctanh(np.clip(dxl / dtl, -0.999999, 0.999999))
    print(f"   cloud: links per element {deg.mean():.2f}; rapidity histogram (bins of 0.5, -4 to 4): {hist}; |eta| < 0.5 against 1.5 < |eta| < 2: {core:.0f} vs {wing:.0f}")
    print(f"   lattice n = 32: links per element {degl.mean():.2f}; fraction of links with |eta| < 1: {(np.abs(etal) < 1).mean():.2f} (cloud {(np.abs(eta) < 1).mean():.2f})")
    ok3 = core < 1.6 * wing and deg.mean() > 4.0 and (np.abs(etal) < 1).mean() < 0.5
    print(f"   flat rapidities and a growing valency on the constrained cloud, the lattice's pile-up absent: {ok3}")
    print()
    print("4. tamper: a jittered lattice (jitter 0.3 of the spacing) under the same constraints")
    ul, vl = jittered_lattice(64, 0.3)
    viol = np.abs(F(ul, vl).mean(0)).max()
    dl = [disp_cells(ul, vl, s) for s in (0.02, 0.05)]; bl = [disp_block(ul, vl, 8, e, 0.08) for e in (0.0, 1.0, 1.5)]
    print(f"   max |mean f| = {viol:.5f} (inside the tight eps = {0.1 / math.sqrt(4096):.5f}); dispersion at cells 0.02, 0.05: {dl[0]:.2f}, {dl[1]:.2f}; boosted small cells eta = 0, 1, 1.5: " + ", ".join(f"{b:.2f}" for b in bl))
    ok4 = viol < 0.1 / math.sqrt(4096) and dl[0] < 0.5 and (max(bl) - min(bl)) > 0.2
    print(f"   passes the constraint, fails independence, selects a frame: the clause, not the rate, is what gives independence: {ok4}")
    print()
    ok = ok1 and ok2a and ok2b and ok2c and ok3 and ok4
    print("SUMMARY")
    print(f"  [{'PASS' if ok1 else 'FAIL'}] maximum entropy given the density is independence (exact on cells)")
    print(f"  [{'PASS' if ok2a and ok2b else 'FAIL'}] the constrained maximum-entropy cloud is Poisson at the grain and hyperuniform only at the constraint scale")
    print(f"  [{'PASS' if ok2c else 'FAIL'}] no frame at the grain under boosts")
    print(f"  [{'PASS' if ok3 else 'FAIL'}] the causal links carry the sprinkling's signature")
    print(f"  [{'PASS' if ok4 else 'FAIL'}] tamper: the rate alone does not give independence")
    print()
    print("  CONCLUSION: read as a law on the realised cloud, P3's maximum-entropy clause supplies the")
    print("  independence comp 208 asked for, by the subadditivity of entropy, and the rate A3 states")
    print("  conditions only the modes above the Mosco scale, leaving the grain a Poisson sprinkling in")
    print("  every frame.  B.7 closes on P3's footing; the fibres of rho (B.6) are a thinning of this cloud.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()