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

Computation 209 -- the classical-scaling placement of the signed-tension metric has no
Euclidean geometry, so the independence of realised positions cannot be tested on it
=======================================================================================
WHAT THIS DOES.  comp 208 reduced the Lorentz invariance of PST's discreteness to one
condition: that the realised positions be independent, Poisson-type, which A2's weak
equidistribution does not fix and which A3's typicality is asked to supply.  The natural
next step was to measure it on an emergence-determined cloud rather than on the surrogate.
The one explicit placement the book has is classical scaling (the multidimensional-scaling
reframing of the realisation map, sec:lorentzian-derivation-explicit) of the signed-tension
metric eq:metric on mu-typical distinction data: directions v(a) i.i.d. uniform on S^6 (as
comp 132 draws them), pairwise tension T({a,b}) = |v(a) + v(b)|, tau_null a low quantile
of the tensions, d(a,b)^2 = |T({a,b}) - tau_null|.  This file builds that placement and
finds that it carries no geometry in which independence could be measured.

THE GRAM FORM HAS THE WRONG SIGN.  With T = sqrt(2 + 2 v_a.v_b) and tau_null below every
T that occurs (a low quantile), d^2 = T - tau_null = (sqrt2 - tau_null) + (sqrt2/2) v_a.v_b
+ O((v_a.v_b)^2).  Classical scaling forms B = -(1/2) J d^2 J, and J kills the constant, so
    B = -(sqrt2/4) G + (higher order),   G the centred Gram matrix of the directions.
G is positive of rank seven (the direction sphere), so B's signal is NEGATIVE of rank
seven: the seven direction axes appear as negative eigenvalues, the kind classical scaling
reads as timelike, and the positive part of B, which is what an embedding into a
Euclidean M would use, is the higher-order remainder.

CHECKS.
  1. SEVEN NEGATIVE EIGENVALUES PROPORTIONAL TO N.  At N = 200, 400, 800, 1600 the seven
     most negative eigenvalues of B are -6.7 to -11.3, -15 to -22, -33 to -44, -67 to -86,
     matching -(sqrt2/4) times the Gram eigenvalues to within the higher-order terms, and
     the eighth is an order of magnitude smaller.
  2. THE EUCLIDEAN PART IS NOISE GROWING LIKE sqrt N.  The largest positive eigenvalue is
     1.5, 2.0, 3.0, 4.9 at the same N: a ratio to the negative signal falling like
     N^(-1/2), so no Euclidean cloud emerges as N grows.
  3. THE DISPERSION TEST HAS NO POWER HERE.  Positions taken from the leading positive
     eigenvectors, cells of equal reference mass, index of dispersion of the counts: for
     i.i.d. directions 0.87 and 1.20 (N = 400, 600; the value expected for independent
     positions with a reference four times the sample is 1.25, the noise about 0.3), and
     for a design evened on the sphere by Riesz-energy descent until the nearest-neighbour
     spread is 0.03 instead of 0.19, 1.09 and 1.65: the evened design reads no lower than
     the i.i.d. one, because the coordinates the test reads are noise.  (Under the flipped
     sign of check 4 the same evened design reads 0.58 against 1.27.)
  4. WITH THE SIGN FLIPPED THE GEOMETRY RETURNS.  Taking d^2 = tau - T (aligned directions
     near, the opposite of eq:metric) makes B's signal +(sqrt2/4) G: seven positive
     eigenvalues proportional to N, the realised positions a pointwise function of each
     site's own direction, and the dispersion test regains its power.  This is a control
     on the method, not a proposal: eq:metric's sign is the book's.

WHAT THIS CHANGES.  Comp 209 does not close B.7; it shows that the explicit placement cannot
be the test of it.  Classical scaling of the signed-tension metric on typical data places
the direction sphere on the negative side of the form and leaves nothing on the positive
side but noise, which is the quantitative form of the book's own remark that the Gram
construction cannot supply the signature (it supplies seven negative directions, not one).
The realisation into M is therefore not this scaling, and the independence of realised
positions stays where comp 208 left it, as A3's, with the explicit placement ruled out as
its instrument.

Tamper: a cut above the bulk of the tensions (the 95 percent quantile) puts d^2 = tau - T
for most pairs and reverses the signal to positive rank seven.  So the sign is set by which
side of the cut the bulk sits on; eq:metric reads T > tau_null as spacelike, which puts the
bulk on the T - tau_null side, and that is the negative case.
"""
import numpy as np

rng = np.random.default_rng(209)


def directions(n):
    v = rng.normal(size=(n, 7)); return v / np.linalg.norm(v, axis=1, keepdims=True)


def even_design(v, steps=300):
    """Riesz-energy descent on S^6: an evened, non-typical design."""
    v = v.copy(); n = len(v)
    for _ in range(steps):
        diff = v[:, None, :] - v[None, :, :]
        r2 = (diff ** 2).sum(-1) + np.eye(n)
        F = (diff / r2[..., None] ** 2).sum(1)
        F -= (F * v).sum(1, keepdims=True) * v
        v = v + 0.1 * F / (np.linalg.norm(F, axis=1, keepdims=True).max() + 1e-12)
        v /= np.linalg.norm(v, axis=1, keepdims=True)
    return v


def nn_spread(v):
    G = v @ v.T; np.fill_diagonal(G, -2); nn = np.arccos(np.clip(G.max(1), -1, 1))
    return nn.std() / nn.mean()


def gram_form(v, sign=+1, q=0.05):
    """B = -(1/2) J d^2 J for d^2 = sign * (T - tau_null) made non-negative as eq:metric does."""
    N = len(v)
    T = np.linalg.norm(v[:, None, :] + v[None, :, :], axis=-1)
    tau = np.quantile(T[np.triu_indices(N, 1)], q if sign > 0 else 1 - q)
    d2 = np.abs(T - tau); np.fill_diagonal(d2, 0)
    J = np.eye(N) - 1.0 / N
    return -0.5 * J @ d2 @ J


def kmeans(X, k, it=30):
    C = X[rng.choice(len(X), k, replace=False)].copy()
    for _ in range(it):
        lab = np.argmin(((X[:, None, :] - C[None, :, :]) ** 2).sum(-1), axis=1)
        for j in range(k):
            if (lab == j).any(): C[j] = X[lab == j].mean(0)
    return C


def dispersion(v_test, v_ref, sign=+1, dim=4, k=40):
    """index of dispersion of the test cloud's counts in cells of equal reference mass, the
    positions taken from the leading positive eigenvectors of the joint Gram form."""
    v = np.vstack([v_test, v_ref]); B = gram_form(v, sign)
    from scipy.linalg import eigh as _eigh
    n = B.shape[0]
    w, V = _eigh(B, subset_by_index=[n - dim, n - 1])     # the top dim eigenpairs only (v32.162; np.linalg.eigh took all)
    idx = np.argsort(w)[::-1][:dim]
    X = V[:, idx] * np.sqrt(np.clip(w[idx], 0, None)); nt = len(v_test)
    C = kmeans(X[nt:], k)
    lab = lambda Y: np.argmin(((Y[:, None, :] - C[None, :, :]) ** 2).sum(-1), axis=1)
    n = np.bincount(lab(X[:nt]), minlength=k); e = nt * np.bincount(lab(X[nt:]), minlength=k) / len(v_ref)
    return float(np.mean((n - e) ** 2 / e))


def main():
    print("Computation 209 -- the classical-scaling placement of the signed-tension metric")
    print("=" * 78)
    print("1. the Gram form's signal is negative of rank seven and grows like N; the Euclidean part like sqrt N")
    neg7, pos1 = {}, {}
    for N in (200, 400, 800, 1600):
        v = directions(N); B = gram_form(v); w = np.sort(np.linalg.eigvalsh(B))
        Gc = (v - v.mean(0)) @ (v - v.mean(0)).T; wg = np.sort(np.linalg.eigvalsh(Gc))[::-1][:7]
        neg7[N] = w[:7]; pos1[N] = w[-1]
        print(f"   N = {N:4d}: seven most negative {np.round(w[:7], 1)}, eighth {w[7]:.1f};  -(sqrt2/4) Gram {np.round(-np.sqrt(2) / 4 * wg, 1)};  largest positive {w[-1]:.1f}")
    ok1 = all(neg7[N][6] < 5 * np.sort(np.linalg.eigvalsh(gram_form(directions(N))))[7] for N in (400, 1600)) and \
          all(abs(neg7[2 * N].mean() / neg7[N].mean() - 2) < 0.3 for N in (200, 400, 800))
    ok2 = all(pos1[2 * N] / pos1[N] < 1.8 for N in (200, 400, 800)) and pos1[1600] / abs(neg7[1600].mean()) < pos1[200] / abs(neg7[200].mean())
    print(f"   seven negative eigenvalues an order above the eighth and doubling with N: {ok1};  the positive part grows slower than N and its ratio to the signal falls: {ok2}")
    print()
    print("2. the dispersion test on the positive-side coordinates has no power (book sign)")
    res = {}
    for nt, nr in ((400, 1600), (600, 2400)):
        vt = directions(nt); vr = directions(nr); ve = even_design(directions(nt))
        res[nt] = (dispersion(vt, vr), dispersion(ve, vr), nn_spread(vt), nn_spread(ve))
        print(f"   N = {nt}, reference {nr}: i.i.d. index {res[nt][0]:.2f} (independent positions expect {1 + nt / nr:.2f}), evened design index {res[nt][1]:.2f};  nearest-neighbour spread {res[nt][2]:.2f} against {res[nt][3]:.2f}")
    ok3 = all(0.6 < r[0] < 2.0 and r[1] > r[0] - 0.3 and r[3] < 0.3 * r[2] for r in res.values())
    print(f"   the evened design reads no lower than the i.i.d. one although it is far more regular: the test has no object here: {ok3}")
    print()
    print("3. control, the sign flipped (aligned directions near): the geometry returns and the test regains power")
    v = directions(800); w = np.sort(np.linalg.eigvalsh(gram_form(v, -1)))
    print(f"   flipped sign, N = 800: seven largest positive {np.round(w[-7:], 1)}, eighth {w[-8]:.1f};  most negative {w[0]:.1f}")
    vt = directions(400); vr = directions(1600); ve = even_design(directions(400))
    di, de = dispersion(vt, vr, -1), dispersion(ve, vr, -1)
    print(f"   flipped sign, N = 400: i.i.d. index {di:.2f}, evened design index {de:.2f}")
    ok4 = w[-7] > 5 * w[-8] and de < di - 0.3
    print(f"   seven positive eigenvalues stand out and the evened design is now detected: {ok4}")
    print()
    print("4. tamper: a cut above the bulk reverses the signal, so the sign is set by the side the bulk sits on")
    v = directions(400); w = np.sort(np.linalg.eigvalsh(gram_form(v, +1, q=0.95)))
    ok5 = w[-7] > 5 * abs(w[-8]) and w[0] > -0.3 * w[-1]
    print(f"   cut at the 95 percent quantile: seven largest positive {np.round(w[-7:], 1)}, eighth {w[-8]:.1f}, most negative {w[0]:.1f}: reversed to positive rank seven: {ok5}")
    ok = ok1 and ok2 and ok3 and ok4 and ok5
    print()
    print("SUMMARY")
    print(f"  [{'PASS' if ok1 else 'FAIL'}] the Gram form of the signed-tension metric is negative of rank seven, proportional to N")
    print(f"  [{'PASS' if ok2 else 'FAIL'}] its Euclidean part is noise: no cloud emerges as N grows")
    print(f"  [{'PASS' if ok3 else 'FAIL'}] the dispersion test cannot tell an evened design from an i.i.d. one on it")
    print(f"  [{'PASS' if ok4 else 'FAIL'}] control: with the sign flipped the geometry and the test's power return")
    print(f"  [{'PASS' if ok5 else 'FAIL'}] tamper: a cut above the bulk reverses the sign; the book's cut is the negative case")
    print()
    print("  CONCLUSION: classical scaling of the signed-tension metric on typical data places the direction")
    print("  sphere on the negative side of the form and leaves noise on the positive side, so it is not a")
    print("  realisation into a Euclidean M and cannot test whether realised positions are independent;")
    print("  that condition stays A3's, with the explicit placement ruled out as its instrument.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()