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

Computation 141 -- The directional marginal over v is FORCED by G_2; only independence is adopted
=================================================================================================
Sharpens Computation 127, which established that P2's irreducible primitive is the
non-balance of the elementary directions, |Sum_{a in D} v(a)| > 0, and closed with
"PST fixes NO measure over v: D -> V_7 ... so 'generic' has no substrate referent and
non-balance is adopted, not derived."

That verdict is narrower than it reads, and this computation locates exactly what is
still adopted.  A measure over v factors into a DIRECTIONAL marginal on S^6 and a RADIAL
marginal on [0, inf).  The directional one is not free: G_2 acts TRANSITIVELY on S^6, so
there is exactly ONE G_2-invariant probability measure there.  The radial one turns out
not to matter at all.  What is left adopted is neither -- it is INDEPENDENCE across
properties, and independence is precisely what the ungroundedness argument asserts on
structural grounds: at substrate level nothing exists that could coordinate the values.

Computation 127 check (2) demonstrated the G_2 refutation "with random SO(7) maps (which
contain the G_2 action)".  This one builds G_2 itself, as the automorphism group of the
octonion cross product, so the refutation is confirmed with genuine automorphisms rather
than with a group that merely contains them.

CHECKS
------
  (1) G_2 IS CONSTRUCTED, NOT PROXIED.  The derivations of the octonion cross product on
      V_7 = Im(O) form a 14-dimensional algebra -- the dimension of g_2 -- and their
      exponentials are orthogonal, of determinant +1, and preserve the cross product.
  (2) THE BALANCED LOCUS IS G_2-INVARIANT (confirms 127 check (2), now with real G_2).
      g(Sum v) = Sum g(v(a)) by linearity, so a balanced assignment stays balanced under
      every automorphism.  No invariance argument can exclude an invariant set, which is
      why the Hurwitz/G_2 EXCLUSION route is refuted and must not be re-attempted.
  (3) G_2 IS TRANSITIVE ON S^6.  The orbit through a unit vector is 6-dimensional, so
      S^6 = G_2/SU(3) and carries a UNIQUE G_2-invariant probability measure.  The
      directional marginal is therefore derived, not adopted -- the same shape of
      argument as the Bernoulli measure on P(D), which is fixed by uniqueness under an
      invariance rather than chosen.
  (4) THE RADIAL MARGINAL IS IRRELEVANT.  With directions drawn from the measure of (3),
      independently across properties, P(Sum v = 0) = 0 for EVERY radial law that puts no
      atom at zero -- unit, uniform, exponential, heavy-tailed, or wildly unequal.  No
      continuity assumption on the radii is needed.
  (5) THE ONE EXCEPTION IS DEGENERATE, AND SELF-CLOSING.  If a property may carry no
      directional content at all, v(a) = 0, with probability p, then balance requires
      every property to draw zero and P(balanced) = p^n exactly.  For p < 1 this vanishes
      as |D| grows, so the exception closes by itself in the large-substrate limit; only
      p = 1 survives, which is not a substrate with directional content.
  (6) INDEPENDENCE IS THE LOAD-BEARING INPUT.  A law that coordinates the values across
      properties can make balance CERTAIN; an independent law makes it null.  This is the
      whole of what remains adopted.

CONCLUSION
----------
"PST fixes no measure over v" overstates the gap.  The directional marginal is forced by
G_2, the radial marginal is provably irrelevant, and what is still assumed is independence
across properties -- not a measure choice but a structural claim, and exactly the claim
that exact cancellation would need a coordinating structure which the substrate does not
provide.  The ungroundedness argument does not replace the genericity argument; it supplies
its missing premise.  The narrowed open question therefore becomes: is independence over v
grounded at substrate level, or is IT the irreducible primitive?
"""
import numpy as np
from scipy.linalg import expm

np.random.seed(141)

# ---------------------------------------------------------------- octonions
# Cross product on Im(O) from the Fano-plane triples.
TRIPLES = [(1, 2, 3), (1, 4, 5), (1, 7, 6), (2, 4, 6), (2, 5, 7), (3, 4, 7), (3, 6, 5)]
STRUCT = np.zeros((7, 7, 7))
for (a, b, c) in TRIPLES:
    i, j, k = a - 1, b - 1, c - 1
    for (x, y, z), s in (((i, j, k), 1), ((j, k, i), 1), ((k, i, j), 1),
                         ((j, i, k), -1), ((k, j, i), -1), ((i, k, j), -1)):
        STRUCT[x, y, z] = s


def cross(u, v):
    return np.einsum('i,j,ijk->k', u, v, STRUCT)


E = np.eye(7)

# --------------------------------------------- (1) build g_2 = Der(Im(O), x)
so7 = []
for i in range(7):
    for j in range(i + 1, 7):
        A = np.zeros((7, 7))
        A[i, j], A[j, i] = 1, -1
        so7.append(A)
so7 = np.array(so7)

rows = []
for A in so7:
    r = []
    for a in range(7):
        for b in range(7):
            r.append(A @ cross(E[a], E[b]) - cross(A @ E[a], E[b]) - cross(E[a], A @ E[b]))
    rows.append(np.concatenate(r))
_, sv, Vt = np.linalg.svd(np.array(rows).T)
g2 = Vt[np.sum(sv > 1e-9):]                      # coefficient vectors spanning g_2

print("(1) G_2 constructed as Der(Im(O), x)")
print(f"    dim = {g2.shape[0]}   (g_2 has dimension 14)")
assert g2.shape[0] == 14


def random_G2(scale=1.0):
    c = g2.T @ np.random.randn(g2.shape[0]) * scale
    return expm(np.tensordot(c, so7, axes=(0, 0)))


g = random_G2()
orth = np.allclose(g.T @ g, np.eye(7), atol=1e-10)
det1 = np.isclose(np.linalg.det(g), 1.0)
perr = max(np.linalg.norm(g @ cross(E[a], E[b]) - cross(g @ E[a], g @ E[b]))
           for a in range(7) for b in range(7))
print(f"    orthogonal={orth}  det=+1 {det1}  preserves the cross product "
      f"(max err {perr:.1e})")
assert orth and det1 and perr < 1e-9

# ------------------------------------- (2) the balanced locus is G_2-invariant
worst = 0.0
for n in (2, 3, 5, 12, 40):
    for _ in range(200):
        V = np.random.randn(n, 7)
        V[-1] = -V[:-1].sum(axis=0)
        worst = max(worst, np.linalg.norm((V @ random_G2().T).sum(axis=0)))
print("\n(2) balanced stays balanced under genuine G_2 automorphisms")
print(f"    worst |Sum g(v)| over 1000 trials = {worst:.1e}")
print("    -> an invariance argument cannot exclude an invariant set;")
print("       the Hurwitz/G_2 EXCLUSION route is refuted (do not re-attempt).")
assert worst < 1e-9

# --------------------------------------------- (3) transitivity on S^6
tang = [np.tensordot(c, so7, axes=(0, 0)) @ E[0] for c in g2]
orbit_dim = np.linalg.matrix_rank(np.array(tang), tol=1e-9)
print("\n(3) G_2 orbit through a unit vector")
print(f"    dim = {orbit_dim}   (S^6 has dimension 6)")
print("    -> transitive, so S^6 carries a UNIQUE G_2-invariant probability")
print("       measure. The directional marginal is DERIVED, not adopted.")
assert orbit_dim == 6


# Parts 4 to 6 each drew an (m, n, 7) array of directions in one go: 90 MB at
# m = 200000, n = 8, and 180 MB with the normalised copy beside it.  The browser
# runtime this book publishes to cannot be relied on for that, so they are drawn a
# block at a time.  The numbers are the same numbers in the same order: splitting a
# randn request does not change the legacy stream, and where the radii are drawn
# AFTER the directions, the generator is wound forward, the radii taken, and then
# rewound with get_state/set_state so the directions repeat exactly.
BLOCK = 20_000


def directions(m, n):
    x = np.random.randn(m, n, 7)
    return x / np.linalg.norm(x, axis=2, keepdims=True)


def resultant_norms(m, n, rad=None, r=None):
    """|Sum_i r_i u_i| per row, u uniform on S^6, without holding all m rows."""
    s1 = None
    if r is None:
        s0 = np.random.get_state()
        for i in range(0, m, BLOCK):
            np.random.randn(min(BLOCK, m - i), n, 7)   # wind past the directions
        r = rad(m, n)
        s1 = np.random.get_state()
        np.random.set_state(s0)
    out = np.empty(m)
    for i in range(0, m, BLOCK):
        k = min(BLOCK, m - i)
        u = directions(k, n)
        out[i:i + k] = np.linalg.norm((u * r[i:i + k, :, None]).sum(axis=1), axis=1)
    if s1 is not None:
        np.random.set_state(s1)
    return out


# ------------------------------- (4) the radial marginal does not matter
print("\n(4) P(|Sum v| < eps) with directions uniform on S^6, for radial laws")
print("    that put no atom at zero.  eps = 1e-9 throughout.")
M = 200_000
RADIAL = {
    "all radii = 1":        lambda m, n: np.ones((m, n)),
    "U(0.5, 1.5)":          lambda m, n: np.random.uniform(.5, 1.5, (m, n)),
    "Exp(1), shifted off 0": lambda m, n: 0.1 + np.random.exponential(1, (m, n)),
    "|Cauchy|, heavy tail": lambda m, n: 0.1 + np.abs(np.random.standard_cauchy((m, n))),
    "10^U(-3,3), 6 decades": lambda m, n: 10 ** np.random.uniform(-3, 3, (m, n)),
}
for name, rad in RADIAL.items():
    hits = []
    for n in (2, 3, 8):
        S = resultant_norms(M, n, rad=rad)
        hits.append(int((S < 1e-9).sum()))
    print(f"    {name:<24} n=2,3,8 -> {hits} balanced out of {M} each")
    assert sum(hits) == 0

# -------------------------- (5) the degenerate exception, and its closure
print("\n(5) if v(a) = 0 is allowed with probability p, balance needs ALL of them")
for p in (0.3, 0.5):
    for n in (2, 3, 5):
        m = 400_000
        r = np.where(np.random.rand(m, n) < p, 0.0,
                     np.random.uniform(.5, 1.5, (m, n)))
        S = resultant_norms(m, n, r=r)
        got, pred = float(np.mean(S < 1e-12)), p ** n
        print(f"    p={p} n={n}: measured {got:.4f}   p^n = {pred:.4f}")
        assert abs(got - pred) < 0.01
print("    and p^n -> 0 as |D| grows, for every p < 1:")
for p in (0.5, 0.9, 0.99):
    print(f"      p={p}: " + "  ".join(f"|D|={n}: {p**n:.1e}"
                                       for n in (6, 20, 100, 1000)))
print("    -> the exception closes by itself; only p=1 survives, which is")
print("       not a substrate carrying directional content.")

# ------------------------------------- (6) independence is what does the work
m, n = 200_000, 8
c_hits = i_hits = 0
for _i in range(0, m, BLOCK):
    V = np.random.randn(min(BLOCK, m - _i), n, 7)
    coord = V - V.mean(axis=1, keepdims=True)    # a law that coordinates the values
    c_hits += int((np.linalg.norm(coord.sum(axis=1), axis=1) < 1e-10).sum())
    i_hits += int((np.linalg.norm(V.sum(axis=1), axis=1) < 1e-10).sum())
p_coord, p_indep = c_hits / m, i_hits / m
print("\n(6) is independence load-bearing?")
print(f"    coordinated law : P(balanced) = {p_coord:.4f}")
print(f"    independent law : P(balanced) = {p_indep:.4f}")
assert p_coord > 0.99 and p_indep == 0.0

print("\n" + "=" * 72)
print("RESULT: the directional marginal is forced by G_2 transitivity, the radial")
print("marginal is irrelevant, and what remains adopted is INDEPENDENCE across")
print("properties -- not a measure choice but the structural claim that nothing at")
print("substrate level coordinates the values.  The ungroundedness argument supplies")
print("the genericity argument's missing premise rather than replacing it.")
print("Narrowed question: is independence over v grounded, or is IT the primitive?")