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

Computation 195 -- Maximum entropy supplies non-degeneracy: under the symmetry package
the maximum-entropy law on T's contingent content has a density, both exceptional sets
are null under it, and the one thing the argument needs is that P3's typicality clause
reach the content of T and not only the configuration
=========================================================================================
THE ITEM.  comp 193 left the structural front with one condition, the non-degeneracy of
T's contingent content (t, v): the magnitude t must avoid the tie set, closed of
codimension one in the admissible cone, and the directions v must avoid the balance
set, closed of codimension seven in (S^6)^|D|.  Both faces fail only on closed
nowhere-dense sets and both fail with positive probability only under DISCRETE content.
sec:non-balance-status says what would close the condition: "an argument about
correlation rather than about direction", since comp 152's all-parallel law mu_pm is
invariant under everything the postulates assert (relabelling of D, per-property
complementation, the diagonal G_2 action) and still balances 31% of the time at
|D| = 6.  comp 193 said the same thing from the other side: one argument, that P3's
maximum-entropy sublimation instantiates non-degenerate content, would close both faces.

THE ARGUMENT, which is a theorem of probability and not of PST.  Fix the invariant
reference measure nu on the content space: the uniform product on (S^6)^|D| for the
directions, Lebesgue measure on the admissible cone (with a scale constraint, since the
cone is unbounded) for the magnitudes.  The Gibbs variational principle says that for
every law P the relative entropy D(P || nu) is non-negative with equality only at
P = nu, and that under any finite set of moment constraints the maximum-entropy law is
an exponential family with a density against nu.  So the maximum-entropy law on the
content, with or without constraints of that kind, is ABSOLUTELY CONTINUOUS, and every
closed nowhere-dense set, the tie set and the balance set among them, is null under
it.  mu_pm is on the other side of the same principle: it is supported on a
seven-dimensional subset of a 6|D|-dimensional space, so its relative entropy is
+infinity, and it is the FARTHEST law from maximum entropy that the symmetry package
admits, not a competitor to it.  The symmetry package does not force non-degeneracy
(comp 152); the symmetry package plus maximum entropy does.

WHAT THE ARGUMENT NEEDS, stated exactly.  P3's clause as the book states it reads:
"among the substrate configurations that cross the threshold, the realised one carries
no unforced density structure beyond what T determines", a maximum-entropy selection
of the CONFIGURATION, consistent with the tension data, that supplies A3.  The content
(t, v) of T is that tension data: the input to the selection, not its output.  So the
argument above closes the condition if and only if the typicality clause is read as
reaching the contingent content of T as well as the configuration.  The book does not
make that reading on its own; this file records that it is the one thing missing, that
it is a reading of P3 and not a new postulate, and that under it the condition would
have exactly the standing A3 has.

CHECKS.
  1.  Directions.  A one-parameter family inside the symmetry package: v(a) drawn from
      the von Mises-Fisher law of concentration kappa about +-u, u uniform, signs
      i.i.d. fair (the signs comp 178 showed complementation forces).  Its entropy,
      in closed form from the Bessel functions, decreases strictly in kappa, with its
      maximum at kappa = 0, the uniform product.  Exact balance has probability zero at
      every finite kappa (a steep power law, slope between 6 and 7 at r = 0.5 to 0.75 on a
      million draws at kappa = 0 and 4, approaching the codimension 7 as r -> 0, which
      comp 193 measured at 6.7 on four million draws), and the binomial
      mass C(|D|, |D|/2) 2^-|D| only at kappa = infinity, which is mu_pm.
  2.  Magnitudes.  The maximum-entropy law on the admissible cone with a scale
      constraint has exponential increments along inclusion (Gibbs); 200 draws on three
      graphs carry only the forced atom.  A lattice of increments with L levels is the
      discrete tamper: its tie fraction falls as L grows and is 1 at L = 1, which is T_ex.
  3.  The reference measure is load-bearing: with the fourteen octonion units as the
      only directions, comp 152's quantised scope, the maximum-entropy law is uniform on
      a finite set and balances with probability 1/14 at |D| = 2 and positively at every
      even |D|.  The continuum of directions P1 supplies is what makes the exceptional
      set null.
  4.  The two faces under one law: a joint maximum-entropy draw of (t, v) satisfies both
      faces on every one of 500 draws, with the distance to degeneracy of comp 193
      positive each time.
A "pass" is all four.
"""
import itertools
import math
import numpy as np
from scipy.special import iv, gammaln


# --- directions: von Mises-Fisher on S^6 --------------------------------------------
P = 7                                  # ambient dimension of V_7
LOG_AREA_S6 = math.log(16 * math.pi ** 3 / 15)

def vmf_entropy(kappa):
    """differential entropy of the vMF law on S^6 against the surface measure, exact."""
    if kappa == 0.0:
        return LOG_AREA_S6
    nu = P / 2 - 1                                                  # 5/2
    log_C = nu * math.log(kappa) - (P / 2) * math.log(2 * math.pi) - math.log(iv(nu, kappa))
    A = iv(P / 2, kappa) / iv(nu, kappa)                            # E[u . x]
    return -log_C - kappa * A

def sample_vmf(kappa, mu, n, rng):
    """Wood's algorithm for vMF on S^{P-1}; kappa = 0 gives the uniform law."""
    if kappa == 0.0:
        v = rng.normal(size=(n, P)); return v / np.linalg.norm(v, axis=1, keepdims=True)
    b = (-2 * kappa + math.sqrt(4 * kappa ** 2 + (P - 1) ** 2)) / (P - 1)
    x0 = (1 - b) / (1 + b)
    c = kappa * x0 + (P - 1) * math.log(1 - x0 ** 2)
    out = np.empty((n, P)); k = 0
    while k < n:
        m = n - k
        z = rng.beta((P - 1) / 2, (P - 1) / 2, size=m)
        w = (1 - (1 + b) * z) / (1 - (1 - b) * z)
        u = rng.uniform(size=m)
        acc = kappa * w + (P - 1) * np.log(1 - x0 * w) - c >= np.log(u)
        w = w[acc]
        if w.size == 0:
            continue
        t = rng.normal(size=(w.size, P)); t -= np.outer(t @ mu, mu); t /= np.linalg.norm(t, axis=1, keepdims=True)
        v = np.sqrt(1 - w[:, None] ** 2) * t + w[:, None] * mu
        take = min(v.shape[0], n - k); out[k:k + take] = v[:take]; k += take
    return out

def draw_family(kappa, n, draws, rng):
    """draws clouds of n directions: v(a) ~ vMF(kappa, eps_a u), u uniform, eps fair."""
    out = np.empty((draws, n, P))
    for d in range(draws):
        u = rng.normal(size=P); u /= np.linalg.norm(u)
        eps = rng.choice([-1.0, 1.0], size=n)
        if np.isinf(kappa):
            out[d] = eps[:, None] * u
        else:
            vs = sample_vmf(kappa, u, n, rng)
            out[d] = eps[:, None] * vs
    return out


def part1():
    print("\n  (1) DIRECTIONS: entropy is maximal at the uniform product, and balance is null at every finite kappa")
    rng = np.random.default_rng(195)
    kappas = [0.0, 0.5, 1.0, 2.0, 4.0, 8.0, 16.0]
    H = [vmf_entropy(k) for k in kappas]
    mono = all(H[i] > H[i + 1] for i in range(len(H) - 1))
    print("      kappa:   " + "".join(f"{k:9.1f}" for k in kappas))
    print("      entropy: " + "".join(f"{h:9.4f}" for h in H) + f"   (log area S^6 = {LOG_AREA_S6:.4f}; strictly decreasing: {mono})")
    n = 8; ok = mono
    print(f"      {'kappa':>8} {'exact balance':>14} {'P(|sum|<0.5)':>13} {'P(|sum|<0.75)':>14} {'slope':>7}")
    for k in (0.0, 4.0, np.inf):
        s = np.concatenate([np.linalg.norm(draw_family(k, n, 200000, rng).sum(axis=1), axis=1) for _ in range(5 if not np.isinf(k) else 1)])
        exact = (s < 1e-9).mean()
        p1, p2 = (s < 0.5).mean(), (s < 0.75).mean()
        slope = math.log(p2 / p1) / math.log(1.5) if p1 > 0 else float('nan')
        if np.isinf(k):
            ok &= abs(exact - math.comb(n, n // 2) / 2 ** n) < 0.01
        else:
            ok &= exact == 0.0 and slope > 5.5
        print(f"      {str(k):>8} {exact:14.4f} {p1:13.6f} {p2:13.6f} {slope:7.2f}")
    print(f"      binomial at kappa = inf is {math.comb(n, n // 2) / 2 ** n:.4f}; mu_pm is the kappa = inf point, relative entropy +inf")
    return ok


# --- magnitudes: exponential increments on the admissible cone -----------------------
def bits(n):
    m = np.arange(1 << n); return np.array([(m >> a) & 1 for a in range(n)], dtype=np.int8)

def independent_mask(n, edges, bt):
    bad = np.zeros(1 << n, dtype=np.int8)
    for u, v in edges:
        bad |= bt[u] & bt[v]
    return bad == 0

def maxent_t(n, indep, rng, levels=None):
    """t = 0 on independent sets; elsewhere the maximum over C minus one property plus an
    increment: exponential (the maximum-entropy increment under a mean constraint) or,
    for the tamper, uniform on a lattice of `levels` values."""
    N = 1 << n; t = np.zeros(N)
    for C in range(N):
        if indep[C]:
            continue
        best = max(t[C & ~(1 << a)] for a in range(n) if C >> a & 1)
        inc = rng.exponential(1.0) if levels is None else rng.integers(1, levels + 1) / levels
        t[C] = best + inc
    return t

def unforced_atom(t, indep):
    vals, counts = np.unique(np.round(t[~indep], 12), return_counts=True)
    return counts.max() / t.size

def part2():
    print("\n  (2) MAGNITUDES: exponential increments carry only the forced atom; a lattice of increments ties")
    rng = np.random.default_rng(1950)
    ok = True
    for n, edges in ((6, [(i, (i + 1) % 6) for i in range(6)]), (8, [e for e in itertools.combinations(range(8), 2) if rng.random() < 0.35]), (8, list(itertools.combinations(range(4), 2)))):
        bt = bits(n); indep = independent_mask(n, edges, bt)
        worst = max(unforced_atom(maxent_t(n, indep, rng), indep) for _ in range(200))
        ok &= worst <= 1 / (1 << n) + 1e-15
        print(f"      |D| = {n}, {len(edges):2d} distinctions: largest unforced atom over 200 maximum-entropy draws = {worst:.4f} (one configuration is {1/(1<<n):.4f})")
    n = 8; edges = [(i, (i + 1) % 8) for i in range(8)]; bt = bits(n); indep = independent_mask(n, edges, bt)
    print(f"      {'levels L':>9} {'mean tie mass':>14}")
    prev = None
    for L in (1, 2, 4, 16, 256):
        m = float(np.mean([unforced_atom(maxent_t(n, indep, rng, levels=L), indep) for _ in range(50)]))
        ok &= (prev is None or m <= prev + 1e-12) and (L > 1 or m > 0.2)
        prev = m
        print(f"      {L:9d} {m:14.4f}" + ("   (L = 1 is T_ex: every increment 1, the edge count)" if L == 1 else ""))
    return ok


def part3():
    print("\n  (3) THE REFERENCE MEASURE IS LOAD-BEARING: quantised directions balance under maximum entropy")
    rng = np.random.default_rng(1951)
    units = np.vstack([np.eye(P), -np.eye(P)])                    # the fourteen octonion imaginary units
    ok = True
    print(f"      {'|D|':>4} {'balance, uniform on 14 units':>29} {'exact':>8}")
    for n in (2, 4, 6):
        idx = rng.integers(0, 14, size=(400000, n))
        s = units[idx].sum(axis=1)
        frac = (np.abs(s).sum(axis=1) < 1e-9).mean()
        exact = 1 / 14 if n == 2 else None
        if exact is not None:
            ok &= abs(frac - exact) < 0.005
        ok &= frac > 0
        print(f"      {n:4d} {frac:29.4f} {('1/14' if n == 2 else '> 0'):>8}")
    print("      maximum entropy on a finite set is uniform on it and has no density: the continuum P1 supplies does the work")
    return ok


def part4():
    print("\n  (4) THE TWO FACES UNDER ONE MAXIMUM-ENTROPY LAW: 500 joint draws")
    rng = np.random.default_rng(1952)
    n = 8; edges = [(i, (i + 1) % 8) for i in range(8)]; bt = bits(n); indep = independent_mask(n, edges, bt)
    dmin = np.inf; ok = True
    for _ in range(500):
        t = maxent_t(n, indep, rng); v = sample_vmf(0.0, None, n, rng)
        s = np.sort(t[~indep]); dt = np.min(np.diff(s)) / math.sqrt(2)
        dv = float(np.linalg.norm(v.sum(axis=0)) / math.sqrt(n))
        d = min(dt, dv); dmin = min(dmin, d); ok &= d > 0
    print(f"      smallest distance to degeneracy over 500 draws: {dmin:.3e} (> 0 on every draw: {ok})")
    return ok


def main():
    print("=" * 100)
    print("  Computation 195 -- Maximum entropy supplies non-degeneracy, if P3's typicality reaches the content of T")
    print("=" * 100)
    checks = [
        ("directions: entropy maximal at the uniform product, balance null at every finite kappa, binomial only at mu_pm", part1()),
        ("magnitudes: exponential increments carry only the forced atom, lattice increments tie", part2()),
        ("quantised directions balance under maximum entropy: the continuum reference measure is load-bearing", part3()),
        ("both faces hold on every joint maximum-entropy draw", part4()),
    ]
    print("\n" + "=" * 100)
    for name, ok in checks:
        print(f"  {'ok  ' if ok else 'FAIL'}  {name}")
    print("=" * 100)
    print("  Under maximum entropy relative to the invariant reference measure the content of T has a density,")
    print("  so the tie set and the balance set are null and both faces of the condition hold almost surely.")
    print("  mu_pm is the farthest law from maximum entropy the package admits.  What the argument needs is")
    print("  that P3's typicality clause reach the content of T and not only the configuration: a reading of")
    print("  P3, not a new postulate, and the book has not made it.")
    print("  PASS" if all(ok for _, ok in checks) else "  FAIL")

main()