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

Computation 216 -- the black-hole interior on the theory's current lengths, and the
central singularity dissolved by the fibre bound
==============================================================================
WHAT THIS DOES.  Two things about the one region the projected equations do not reach.

(A) THE NUMBERS WERE ON THE OLD LENGTH.  comp 154 evaluated the book's validity condition
(k d_0)^2 << 1 on Schwarzschild with d_0 = 2.3 nm, the Casimir bound of its day, and the
book's black-hole chapter still prints those values: (k d_0)^2 = 2.1e-24 at a solar-mass
horizon, r_break = 38 micrometres, 1.3e-8 of r_s, 2e-24 of the volume, a mass floor of
2.9e18 kg.  Since comp 201 the coherence length is d_0 = sqrt 2 hbar c / m_h = 2.2 am, nine
orders smaller, and the audit of v32.149 (section A) caught the Casimir chapter and the
Conclusion but not this chapter.  Re-evaluated on d_0 = 2.2 am: (k d_0)^2 = 1.9e-42 at the
horizon, r_break = 37 picometres = 1.2e-14 r_s, enclosing 1.9e-42 of the interior volume,
mass floor 2.8e9 kg; on the grain l_* = hbar c / M_* = 1.5e-19 m: 9.4e-45, 6.2 pm, 1.9e8 kg
(the grain's floor was already right in the book, so the chapter mixed one length's
numbers with the other's).  Every conclusion of comp 154 survives, by wider margins: the
uncovered region is the singularity's neighbourhood, not the horizon.  With d_0 = 2.3 nm
the old numbers are reproduced exactly (tamper i), so the discrepancy is the length and
nothing else.

(B) THE SINGULARITY HAS NO SUBSTRATE COUNTERPART.  The book declines to transfer the
initial-singularity argument ("loci are not primitive") to the black hole's centre, and
says the argument has not been made.  It can be made in a different form, from comp 212.
A spacetime point is a fibre of rho, and a fibre is a clique of pairwise-null properties,
which in V_7 holds at most seven members (the Gram rank), fewer for an obtuse null angle,
and is a singleton almost surely on the theory's measure.  So the realised content of a
collapsing body, n properties with n far above seven, cannot occupy one point: some pair
must fail to be null, and its proper distance |T - tau_null|^{1/4} (eq. metric) is the
region's diameter.  How far from null that pair must be is a packing bound.  With the
G2-invariant pair tension a function of the angle (T = |v_a + v_b| = sqrt(2 + 2 cos) as
comp 212's exemplar) and tau_null at cos theta_0 = c_0:
   Rankin:  among n unit vectors some pair has cos >= -1/(n-1), so if c_0 < 0 the
            deviation is at least |c_0| - 1/(n-1)  (> 0 for n > 1 + 1/|c_0|);
   Welch:   among n > 7 unit vectors in R^7 some pair has |cos| >= sqrt((n-7)/(7(n-1))),
            so the deviation is at least that minus |c_0|  (1/7 at n = 8 and c_0 = 0,
            achieved by the simplex).
Part B evaluates both bounds and, as upper bounds on the same minimum, the best
configurations a soft-max descent on the sphere finds for n = 4 to 50 at the two
placements comp 212 uses (the 5 percent quantile, c_0 = -0.62, and the median, c_0 = 0).
At the median seven properties can be exactly null (an orthonormal frame; the descent
reaches 7e-4), eight cannot (Welch's 1/7 = 0.143, found to 0.1429), and the deviation
grows with n: 0.25, 0.28, 0.39, 0.45, 0.58 at n = 10, 14, 20, 30, 50; at the obtuse
placement only two can.  In the proper distance the diameter of a region holding n
properties is at least 0.56 (n = 8) to 0.80 (n = 50) at the median and 0.76 to 1.05 at the
obtuse placement, in the units of the exemplar tension, and it saturates: the pair
tension is bounded (0 <= T <= 2), so the single-pair distance is bounded by (2 - tau)^{1/4},
which is the local character of eq. metric and not a statement about M at large.

WHAT THIS CHANGES.  The chapter's numbers move to the current lengths.  The open item
"no result on the central singularity" becomes narrower: the singularity as a point
carrying the collapsed content is excluded on the substrate side, by the fibre bound and
the packing bounds, with the same strength as comp 212's fibre result and no more; what
replaces it, a region of n properties whose minimal diameter the bounds give in local
units and whose dynamics the projected equations do not reach, the theory does not
describe.  The statistical entropy, the information question and the horizon's status
as an event horizon are untouched.  The uncovered region holds about 1e25 grains for a
solar mass, 1e-42 of the hole's volume.

Tampers: (i) d_0 = 2.3 nm reproduces every printed number; (ii) the simplex attains
Welch at n = 8, so the bound is sharp where it can be checked; (iii) at c_0 = -0.62 the
descent's n = 4 and 7 values equal Rankin's bound to 1e-3, so the bound is attained
there too (the regular simplex, which embeds up to n = 8).

Runs in about a minute; numpy only.
"""
import math
import numpy as np

G_N, C_L, HBAR_C = 6.674e-11, 2.998e8, 1.97327e-7     # SI; hbar c in eV m
M_SUN = 1.989e30
D0_OLD, D0_NEW = 2.3e-9, 2.2e-18                        # the Casimir bound of comp 154's day; comp 201's value
L_GRAIN = HBAR_C / 1.285e12                             # hbar c / M_*, M_* = 1.285 TeV
rng = np.random.default_rng(216)


def horizon_numbers(d, M=M_SUN):
    rs = 2 * G_N * M / C_L ** 2
    supp = math.sqrt(12) * (d / rs) ** 2
    rb = (12 ** 0.25 * d) ** (2 / 3) * rs ** (1 / 3)
    floor = 12 ** 0.25 * d * C_L ** 2 / (2 * G_N)
    return dict(rs=rs, supp=supp, rb=rb, frac=rb / rs, vol=(rb / rs) ** 3, floor=floor)


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


def null_cos(q, n=4000):
    v = directions(n); G = (v @ v.T)[np.triu_indices(n, 1)]
    return float(np.quantile(G, q))


def min_deviation(n, c0, restarts=10, iters=2000, beta0=5.0):
    """Upper bound on min over n unit vectors in R^7 of max_pairs |cos - c0|: soft-max descent on the sphere."""
    best = np.inf; iu = np.triu_indices(n, 1)
    for _ in range(restarts):
        v = directions(n); lr = 0.05
        for it in range(iters):
            beta = beta0 * (1 + it / 300)
            D = v @ v.T - c0; a = np.abs(D[iu]); w = np.exp(beta * (a - a.max())); w /= w.sum()
            W = np.zeros((n, n)); W[iu] = w * np.sign(D[iu]); W = W + W.T
            g = W @ v; g -= (g * v).sum(1, keepdims=True) * v
            v -= lr * g; v /= np.linalg.norm(v, axis=1, keepdims=True)
        best = min(best, float(np.abs(v @ v.T - c0)[iu].max()))
    return best


def rankin(n, c0):
    return max(0.0, -c0 - 1 / (n - 1)) if c0 < 0 else 0.0


def welch(n, d=7):
    return math.sqrt((n - d) / (d * (n - 1))) if n > d else 0.0


def main():
    print("=" * 100)
    print("Computation 216: the interior on the current lengths, and the singularity dissolved by the fibre bound")
    print("=" * 100)
    checks = {}
    print("PART A.  comp 154's horizon numbers on the three lengths (one solar mass)")
    rows = {}
    for name, d in (("d_0 = 2.3 nm (comp 154's input, the old Casimir bound)", D0_OLD), ("d_0 = 2.2 am (comp 201)", D0_NEW), ("grain l_* = hbar c / M_*", L_GRAIN)):
        h = horizon_numbers(d); rows[name] = h
        print(f"  {name:52s}: (k d)^2 at r_s = {h['supp']:.1e}; r_break = {h['rb']:.2e} m = {h['frac']:.1e} r_s; volume fraction {h['vol']:.1e}; mass floor {h['floor']:.1e} kg")
    old = rows["d_0 = 2.3 nm (comp 154's input, the old Casimir bound)"]; new = rows["d_0 = 2.2 am (comp 201)"]; gr = rows["grain l_* = hbar c / M_*"]
    checks['old'] = abs(old['supp'] / 2.1e-24 - 1) < 0.05 and abs(old['rb'] / 38e-6 - 1) < 0.02 and abs(old['floor'] / 2.9e18 - 1) < 0.02
    checks['new'] = abs(new['supp'] / 1.9e-42 - 1) < 0.05 and abs(new['rb'] / 37e-12 - 1) < 0.02 and abs(new['floor'] / 2.8e9 - 1) < 0.03 and abs(gr['floor'] / 1.9e8 - 1) < 0.03
    grains = (4 / 3) * math.pi * (new['rb'] / L_GRAIN) ** 3
    print(f"  the printed 2.1e-24, 38 um, 2.9e18 kg are the 2.3 nm row: {checks['old']}; on 2.2 am they read 1.9e-42, 37 pm, 2.8e9 kg (grain floor 1.9e8 kg as printed): {checks['new']}")
    print(f"  the uncovered region holds about {grains:.1e} grains (radius {new['rb'] / L_GRAIN:.1e} grains), {new['vol']:.0e} of the hole's volume")
    print()
    print("PART B.  the minimal diameter of a region holding n properties, on the theory's measure")
    okB = True; table = {}
    for q in (0.05, 0.5):
        c0 = null_cos(q); tau = math.sqrt(2 + 2 * c0)
        print(f"  tau_null at quantile {q:.2f}: cos theta_0 = {c0:+.3f}, tau_null = {tau:.3f}; a fibre holds at most {7 if c0 > -1 / 6 else int(1 + 1 / abs(c0))}")
        print(f"      n   Rankin  Welch   found   -> deviation in T   diameter |T - tau|^(1/4)")
        for n in (2, 4, 7, 8, 10, 14, 20, 30, 50):
            lb = max(rankin(n, c0), welch(n) - abs(c0))      # some pair has |cos| >= Welch, so |cos - c0| >= Welch - |c0|
            found = min_deviation(n, c0, restarts=6 if n > 20 else 10)
            dT = found / tau; table[(q, n)] = (lb, found, dT ** 0.25)
            print(f"     {n:3d}   {rankin(n, c0):.3f}   {max(0.0, welch(n) - abs(c0)):.3f}   {found:.4f}        {dT:.4f}             {dT ** 0.25:.3f}")
            okB &= found >= lb - 2e-3
    med7 = table[(0.5, 7)][1]; med8 = table[(0.5, 8)][1]; obt2 = table[(0.05, 2)][1]; obt4 = table[(0.05, 4)][1]
    checks['bounds'] = okB
    checks['fibre'] = med7 < 1e-2 and abs(med8 - 1 / 7) < 3e-3 and obt2 < 1e-2 and obt4 > 0.25
    grow = all(table[(0.5, a)][1] < table[(0.5, b)][1] for a, b in ((8, 10), (10, 14), (14, 20), (20, 30), (30, 50)))
    checks['grows'] = grow
    print(f"  at the median seven can be exactly null and eight cannot (Welch 1/7): {checks['fibre']}; the deviation grows with n: {grow}; every found value respects its lower bound: {okB}")
    print(f"  the single-pair distance saturates at (2 - tau)^(1/4) = {(2 - math.sqrt(2)) ** 0.25:.3f} (median), the local character of eq. metric")
    print()
    print("TAMPERS")
    t1 = checks['old']
    print(f"  (i)   d_0 = 2.3 nm reproduces the printed numbers: {t1}")
    t2 = abs(table[(0.5, 8)][1] - welch(8)) < 3e-3
    print(f"  (ii)  the simplex attains Welch at n = 8: found {table[(0.5, 8)][1]:.4f} against 1/7 = {welch(8):.4f}: {t2}")
    c0 = null_cos(0.05); t3 = abs(table[(0.05, 4)][1] - rankin(4, c0)) < 3e-3 and abs(table[(0.05, 7)][1] - rankin(7, c0)) < 3e-3
    print(f"  (iii) at the obtuse placement Rankin is attained at n = 4 and 7: found {table[(0.05, 4)][1]:.4f}, {table[(0.05, 7)][1]:.4f} against {rankin(4, c0):.4f}, {rankin(7, c0):.4f}: {t3}")
    checks['tampers'] = t1 and t2 and t3
    print()
    ok = all(checks.values())
    print("SUMMARY")
    print(f"  [{'PASS' if checks['old'] and checks['new'] else 'FAIL'}] the chapter's horizon numbers are the 2.3 nm values; on 2.2 am the uncovered region is 37 pm, 1.9e-42 of the volume, floor 2.8e9 kg")
    print(f"  [{'PASS' if checks['fibre'] else 'FAIL'}] seven properties can share a point at the median placement and eight cannot; two at the obtuse one")
    print(f"  [{'PASS' if checks['bounds'] and checks['grows'] else 'FAIL'}] the minimal diameter of n properties is bounded below by Rankin and Welch and grows with n")
    print(f"  [{'PASS' if checks['tampers'] else 'FAIL'}] tampers rejected")
    print()
    print("  CONCLUSION: the region the projected equations do not reach is six orders smaller than printed, and")
    print("  the singularity at its centre has no substrate counterpart: a point carries at most seven properties,")
    print("  so a collapsed body's content occupies a region whose minimal diameter the packing bounds give in")
    print("  local units.  What that region does, the theory does not say; nothing else about the hole changes.")
    print("  ALL CHECKS PASSED" if ok else "  SOME CHECKS FAILED")


if __name__ == "__main__":
    main()