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

Computation 192 -- The time function the slab is cut in: isotropy makes Phi(tau-hat)
hypersurface-orthogonal, its level sets are the Cauchy surfaces, and the slab's
normalisation is the proper 3-volume, invariant under relabelling the time
========================================================================================
THE ITEM.  comp 191 left the slab-wise CAR limit resting on two things: A3's
equidistribution rate for the emergence-determined cloud, and "the time function in
which the slab is cut, in which that cloud must be 4-volume equidistributed", which the
book had not constructed.  comp 191's part 6 showed why the second is not cosmetic: a
cloud equidistributed in space but not in the time coordinate returns the time
marginal's density at the slice rather than the slice integral.  This file constructs
the time function from what the book already has, proves that the construction does
not depend on the one choice left in it, and shows that the normalisation comp 191
assumed is the metric's own.

WHAT THE BOOK ALREADY HAS, in order.  (1) P3 fixes the modal-sublimation direction
tau-hat in V_7 globally at threshold crossing.  (2) The time-orientation passage of
Chapter 5 maps it, by permutation equivariance and Mosco continuity, to a continuous
nowhere-vanishing timelike field on M = R x S^3, "Phi(tau-hat) = d/dt along the R
factor".  (3) The Cosmological Principle theorem makes the large-scale geometry
homogeneous and isotropic, hence of Friedmann-Robertson-Walker form, eq:frw-metric,
ds^2 = dt^2 - a(t)^2 dOmega_3^2.  (4) A3 says the realised cloud is equidistributed
against the 4-volume form sqrt(-g) d^4x = a(t)^3 dt dOmega_3.  Nothing below is new
input; each part is one of these four, made to do the work it was already doing.

PART 1, THE FOLIATION IS FORCED.  A timelike unit field u is tangent to a foliation by
hypersurfaces if and only if it is hypersurface-orthogonal, u ^ du = 0 (Frobenius),
equivalently its twist omega_mu = eps_{mu nu rho sigma} u^nu nabla^rho u^sigma vanishes.
The twist is a spatial vector field, orthogonal to u.  Isotropy about every point, which
the Cosmological Principle supplies, means no spatial vector field is invariant unless it
is zero: so the twist vanishes, Phi(tau-hat) is hypersurface-orthogonal, and its
orthogonal hypersurfaces ARE the level sets t = const.  Checked symbolically on
eq:frw-metric: u = dt has u ^ du = 0 identically for every a(t); and the check
discriminates, since a rotating field u = dt + Omega r^2 dphi (the Goedel-type tamper,
which breaks isotropy) has u ^ du = 2 Omega r dr ^ dphi ^ dt, not zero, so no foliation
and no slab.  Part 1 also averages a random spatial vector over SO(3) at a point and
gets zero, which is the isotropy step made numerical rather than asserted.

PART 2, THE LEVEL SETS ARE COMPACT CAUCHY SURFACES.  {t} x S^3 is compact, and it is
Cauchy if every inextendible causal curve meets it exactly once.  Along a causal curve
dt/dlambda = g(u, xdot) is nonzero (u timelike, xdot causal, both future-directed), so t
is strictly monotone and no causal curve meets a level set twice; and t runs over the
whole range of the expansion, so each curve meets each level set once.  Checked by
integrating null geodesics on the book's instantiation-dominated solution a = exp(H t),
and on a matter-plus-instantiation solution, from random initial points and directions:
each crosses each of a set of level sets exactly once, with the crossing counted from
the sign changes of t - t_k.  comp 4's entry asserted the Cauchy surfaces; here the
assertion is verified for the time function being constructed.

PART 3, THE NORMALISATION IS THE PROPER 3-VOLUME.  A3 equidistributes the cloud against
a(t)^3 dt dOmega_3, so the cloud's time marginal is proportional to a(t)^3 and NOT
uniform in t; by comp 191's part 6 a uniform normalisation would then return
a(t)^3 / <a^3> times the slice value.  The slab S(t_0) = {|t - t_0| <= h/2} has proper
3-volume a(t_0)^3 V_3 per slice and proper thickness h (t is proper time along u), so
the covariant normalisation is nu = lambda a(t_0)^3 V_3 h with lambda = N / integral
a^3 dt dOmega_3 the 4-density.  Measured on a Halton cloud pushed through the cumulative
distribution of a^3: the covariant estimator lands on the slice integral at every t_0
tested while the uniform one is off by exactly a(t_0)^3 / <a^3>, from 0.4 to 2.3 across
the range.  This is comp 191's tamper read the right way round: the modulation is the
scale factor, and the metric absorbs it.

PART 4, THE ONE CHOICE LEFT DOES NOT MATTER.  The foliation is fixed by Part 1; the
labelling of its leaves, t -> f(t) with f monotone, is not.  Relabel: the slab of
coordinate thickness h' in f-time is the slab of thickness h'/f'(t_0) in t, and the
proper-volume normalisation, which is built from the metric and so from f', absorbs the
factor.  Measured with f(t) = t + 0.1 sin(2 pi t), monotone with f' between 0.37 and
1.63, on the same cloud: the covariant estimator is unchanged to the discrepancy level
while the uniform one moves by 1/f'(t_0).  (A first draft used amplitude 0.3, for which
f' changes sign, so f is not a relabelling at all; the check caught it, since a
non-monotone f maps one f-slab onto several t-slabs and nothing is invariant.)  So the time function is defined up to the one
freedom that carries no content, and the slab-wise CAR limit is invariant under it.

WHAT THIS CLOSES.  The time function comp 191 asked for is cosmic time: the parameter
along Phi(tau-hat) whose level sets are the isotropy-orthogonal slices, with the slab
normalised by the proper 3-volume of the slice.  With comp 191, the slab-wise CAR limit
holds for the emergence-determined cloud given exactly the conditions the Mosco chain
already carries, A1 to A4 with A3's rate supplied by P3, and nothing else; and with
comp 185, W4 follows.  That is the same standing equicoercivity A7 has, "a theorem given
A1, A2, A4 and A3", and the register counts A7 closed.  So the item is closed on the
same terms, and the Appendix C row says so in those words.

CHECKS.
  1.  u = dt on eq:frw-metric: u ^ du = 0 symbolically for symbolic a(t); the rotating
      tamper gives a nonzero 3-form; the SO(3) average of a random spatial vector is 0
      to 1e-12.
  2.  Null geodesics on a = exp(H t) and on a matter-plus-instantiation a(t): 200 curves,
      each crossing each of 6 level sets exactly once.
  3.  Covariant slab estimator within 1% of the slice integral at 5 slices of the
      a^3-weighted Halton cloud; uniform estimator within 2% of a(t_0)^3/<a^3> times it.
  4.  Under t -> t + 0.1 sin(2 pi t): covariant estimator unchanged within 1%; uniform
      estimator moves by 1/f'(t_0) within 2%.
A "pass" is all four.
"""
import math
import numpy as np
import sympy as sp

CHUNK = 1 << 17


# --- part 1: Frobenius -------------------------------------------------------------
def wedge_u_du(metric, coords, u_cov):
    """u ^ du as the antisymmetric 3-tensor (u ^ du)_{abc}, for a covector u on coords."""
    n = len(coords)
    du = [[sp.diff(u_cov[b], coords[a]) - sp.diff(u_cov[a], coords[b]) for b in range(n)] for a in range(n)]
    out = {}
    for a in range(n):
        for b in range(a + 1, n):
            for c in range(b + 1, n):
                val = (u_cov[a] * du[b][c] - u_cov[b] * du[a][c] + u_cov[c] * du[a][b])
                out[(a, b, c)] = sp.simplify(val)
    return out

def part1():
    print("\n  (1) THE FOLIATION IS FORCED: u = dt is hypersurface-orthogonal, a rotating field is not")
    t, chi, th, ph = sp.symbols('t chi theta phi', real=True)
    a = sp.Function('a')(t)
    coords = [t, chi, th, ph]
    # FRW on R x S^3 (k = +1): ds^2 = dt^2 - a^2 (dchi^2 + sin^2 chi (dtheta^2 + sin^2 theta dphi^2))
    u = [sp.Integer(1), 0, 0, 0]                       # u = dt, the covector of Phi(tau-hat) = d/dt
    w = wedge_u_du(None, coords, u)
    frw_zero = all(v == 0 for v in w.values())
    print(f"      FRW, symbolic a(t): every component of u ^ du is 0: {frw_zero}")
    # the tamper: a rotating frame field u = dt + Omega r^2 dphi in cylindrical coordinates
    tt, r, phi, z, Om = sp.symbols('t r phi z Omega', real=True)
    w2 = wedge_u_du(None, [tt, r, phi, z], [sp.Integer(1), 0, Om * r ** 2, 0])
    nonzero = {k: v for k, v in w2.items() if v != 0}
    print(f"      rotating tamper u = dt + Omega r^2 dphi: nonzero components {nonzero}")
    tamper_ok = nonzero == {(0, 1, 2): sp.simplify(2 * Om * r)} or (len(nonzero) == 1 and sp.simplify(list(nonzero.values())[0] - 2 * Om * r) == 0)
    # isotropy kills a spatial vector: average over SO(3) (Haar, via random rotations)
    rng = np.random.default_rng(192)
    v = rng.normal(size=3)
    acc = np.zeros(3)
    n_rot = 20000
    for _ in range(n_rot):
        q = rng.normal(size=4); q /= np.linalg.norm(q)
        w0, x, y, z_ = q
        R = np.array([[1 - 2 * (y * y + z_ * z_), 2 * (x * y - z_ * w0), 2 * (x * z_ + y * w0)],
                      [2 * (x * y + z_ * w0), 1 - 2 * (x * x + z_ * z_), 2 * (y * z_ - x * w0)],
                      [2 * (x * z_ - y * w0), 2 * (y * z_ + x * w0), 1 - 2 * (x * x + y * y)]])
        acc += R @ v
    mean = acc / n_rot
    # Monte Carlo over Haar: the exact average is 0, the sample mean is O(n^-1/2)
    iso_ok = np.linalg.norm(mean) < 5 / math.sqrt(n_rot)
    print(f"      SO(3) average of a random spatial vector over {n_rot} rotations: |mean| = {np.linalg.norm(mean):.1e} (exact 0; sampling floor {1/math.sqrt(n_rot):.1e})")
    print("      so the twist of Phi(tau-hat), a spatial vector invariant under isotropy, is zero, and the level sets of t are the foliation")
    return frw_zero and tamper_ok and iso_ok


# --- part 2: Cauchy surfaces -------------------------------------------------------
def null_geodesic_crossings(a_of_t, t0, t1, levels, n_curves, rng):
    """Radial null geodesics on R x S^3 in FRW: chi(t) = chi_0 +- int dt / a(t).  The
    crossing count of each level set is read from the sign changes of t(lambda) - t_k
    along the curve, with lambda an affine parameter, dt/dlambda = E / a(t)."""
    ok = True
    for _ in range(n_curves):
        ts = t0 + (t1 - t0) * rng.random()
        E = 1.0 + rng.random()
        sgn = 1.0 if rng.random() < 0.5 else -1.0
        # integrate dt/dlambda = E / a (future-directed), both ways from the start to cover the range
        def march(direction):
            t, lam, out = ts, 0.0, [ts]
            dl = 1e-3
            while t0 - 1e-9 <= t <= t1 + 1e-9 and lam < 1e4:
                t += direction * dl * E / a_of_t(t); lam += dl; out.append(t)
            return out
        path = march(-1.0)[::-1] + march(+1.0)[1:]
        path = np.array(path)
        for tk in levels:
            s = np.sign(path - tk)
            crossings = int(np.sum(s[1:] * s[:-1] < 0))
            if crossings != 1:
                ok = False
    return ok

def part2():
    print("\n  (2) THE LEVEL SETS ARE COMPACT CAUCHY SURFACES: every null geodesic crosses each once")
    rng = np.random.default_rng(1920)
    H = 1.0
    levels = np.linspace(0.3, 2.7, 6)
    ok1 = null_geodesic_crossings(lambda t: math.exp(H * t), 0.0, 3.0, levels, 100, rng)
    # matter plus instantiation: a(t) = (sinh(3Ht/2))^(2/3), the LCDM-type solution
    ok2 = null_geodesic_crossings(lambda t: math.sinh(1.5 * H * t) ** (2 / 3), 0.05, 3.0, levels, 100, rng)
    print(f"      a = exp(Ht), 100 null geodesics, 6 level sets: each crossed exactly once: {ok1}")
    print(f"      a = sinh(3Ht/2)^(2/3) (matter + instantiation), 100 null geodesics: each crossed exactly once: {ok2}")
    print("      (dt/dlambda = E/a > 0 along every future-directed causal curve, so t is monotone and the crossing is unique)")
    return ok1 and ok2


# --- parts 3 and 4: the slab on an a^3-weighted cloud ------------------------------
def radical_inverse(idx, base):
    x = np.zeros(idx.shape); f = 1.0 / base; i = idx.copy()
    while i.any():
        x += f * (i % base); i //= base; f /= base
    return x

def F(t, x):
    """comp 186's and comp 191's smearing integrand, even about each slice value in t
    once shifted, so the bias is the same on every slice."""
    f = 1.0 + np.sin(2 * np.pi * x[:, 0])
    g = 2.0 + np.cos(2 * np.pi * x[:, 1])
    return f * g

TARGET = 2.0          # the slice integral of F over the unit spatial cube, at every t

def a3_cloud(N, H, t0, t1):
    """Halton cloud in the unit 4-box, with the time coordinate pushed through the CDF of
    a(t)^3 = exp(3Ht) on [t0, t1]: A3's 4-volume equidistribution in coordinates (t, x)."""
    idx = np.arange(1, N + 1, dtype=np.int64)
    u = radical_inverse(idx, 2)
    x = np.column_stack([radical_inverse(idx, b) for b in (3, 5, 7)])
    k = 3 * H
    t = np.log(np.exp(k * t0) + u * (np.exp(k * t1) - np.exp(k * t0))) / k
    return t, x

def part3():
    print("\n  (3) THE NORMALISATION IS THE PROPER 3-VOLUME: covariant versus uniform, on an a^3-weighted cloud")
    H, t0, t1, N, h = 1.0, 0.0, 1.0, 1 << 20, 0.02
    t, x = a3_cloud(N, H, t0, t1)
    k = 3 * H
    four_vol = (np.exp(k * t1) - np.exp(k * t0)) / k          # integral of a^3 dt over the unit cube
    lam = N / four_vol
    mean_a3 = four_vol / (t1 - t0)
    print(f"      {'t_0':>6} {'a^3':>8} {'covariant':>11} {'uniform':>9} {'a^3/<a^3>':>10} {'ratio':>7}")
    ok = True
    for tc in (0.1, 0.3, 0.5, 0.7, 0.9):
        m = np.abs(t - tc) <= h / 2
        s = float(np.sum(F(t[m], x[m])))
        cov = s / (lam * math.exp(k * tc) * h)       # proper 3-volume of the unit cube is a^3, times proper thickness h
        uni = s / ((N / (t1 - t0)) * h)              # the uniform normalisation comp 191 used on a flat cloud
        pred = math.exp(k * tc) / mean_a3
        ok &= abs(cov / TARGET - 1) < 0.01 and abs(uni / (TARGET * pred) - 1) < 0.02
        print(f"      {tc:6.2f} {math.exp(k*tc):8.3f} {cov:11.4f} {uni:9.4f} {pred:10.3f} {uni/TARGET/pred:7.3f}")
    print(f"      the covariant estimator lands on {TARGET} at every slice; the uniform one is off by a(t_0)^3/<a^3>: {ok}")
    return ok

def part4():
    print("\n  (4) RELABELLING THE LEAVES, t -> f(t) = t + 0.1 sin(2 pi t), monotone: the covariant estimator does not move")
    H, t0, t1, N, h = 1.0, 0.0, 1.0, 1 << 20, 0.02
    t, x = a3_cloud(N, H, t0, t1)
    k = 3 * H
    lam = N / ((np.exp(k * t1) - np.exp(k * t0)) / k)
    f = lambda s: s + 0.1 * np.sin(2 * np.pi * s)        # monotone: f' = 1 + 0.2 pi cos(2 pi s) in [0.37, 1.63]
    fp = lambda s: 1 + 0.1 * 2 * np.pi * np.cos(2 * np.pi * s)
    tf = f(t)
    print(f"      {'t_0':>6} {'f(t_0)':>7} {'f(t_0)':>8} {'cov in t':>9} {'cov in f':>9} {'uni in t':>9} {'uni in f':>9} {'ratio':>7}")
    ok = True
    for tc in (0.1, 0.3, 0.5, 0.7, 0.9):
        m = np.abs(t - tc) <= h / 2
        cov_t = float(np.sum(F(t[m], x[m]))) / (lam * math.exp(k * tc) * h)
        uni_t = float(np.sum(F(t[m], x[m]))) / ((N / (t1 - t0)) * h)
        # in f-time: the slab of coordinate thickness h about f(t_c); its proper thickness is h / f'(t_c),
        # and the metric in f-time carries g_ff = 1/f'^2, so the covariant normalisation uses sqrt(g_ff) h = h / f'
        mf = np.abs(tf - f(tc)) <= h / 2
        cov_f = float(np.sum(F(t[mf], x[mf]))) / (lam * math.exp(k * tc) * h / fp(tc))
        uni_f = float(np.sum(F(t[mf], x[mf]))) / ((N / (t1 - t0)) * h)
        ok &= abs(cov_f / cov_t - 1) < 0.01 and abs(uni_f / uni_t / (1 / fp(tc)) - 1) < 0.02
        print(f"      {tc:6.2f} {f(tc):7.3f} {fp(tc):8.3f} {cov_t:9.4f} {cov_f:9.4f} {uni_t:9.4f} {uni_f:9.4f} {uni_f/uni_t*fp(tc):7.3f}")
    print(f"      covariant unchanged within 1%, uniform moved by 1/f'(t_0): {ok}")
    return ok


def main():
    print("=" * 100)
    print("  Computation 192 -- The time function: forced foliation, Cauchy level sets, proper-volume slab")
    print("=" * 100)
    checks = [
        ("Phi(tau-hat) = d/dt is hypersurface-orthogonal on FRW, the rotating tamper is not, isotropy kills the twist", part1()),
        ("the level sets of t are compact Cauchy surfaces: each null geodesic crosses each once", part2()),
        ("the slab normalisation is the proper 3-volume; the uniform one is off by a^3/<a^3>", part3()),
        ("relabelling the leaves leaves the covariant estimator unchanged", part4()),
    ]
    print("\n" + "=" * 100)
    for name, ok in checks:
        print(f"  {'ok  ' if ok else 'FAIL'}  {name}")
    print("=" * 100)
    print("  The time function is cosmic time along Phi(tau-hat), forced by isotropy up to a relabelling that")
    print("  carries no content; its level sets are the Cauchy surfaces; the slab is normalised by the proper")
    print("  3-volume.  With comps 185 and 191 the slab-wise CAR limit, and W4 with it, is a theorem given")
    print("  A1 to A4 with A3 supplied by P3: the standing of equicoercivity A7.")
    print("  PASS" if all(ok for _, ok in checks) else "  FAIL")

main()