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

Computation 163, the action principle Pi_mn does have, and what M4's ambiguity is
==================================================================================
comp 162 left two things.  An exact conservation result, which it said "needs an
action principle the construction does not have", and backlog M4, the divergence-free
ambiguity.  Both are addressed here, and the first sentence needs correcting: the
construction does not have an action of ITS OWN, but it does not need one, because the
identification the book already makes hands it somebody else's.

1.  THE EXACT RESULT, AND WHERE THE ACTION COMES FROM
------------------------------------------------------
\\S sec:einstein identifies Pi_mn as the CONFORMAL PART of T^(fields)_mn, the stress
tensor obtained by varying S_matter.  Diffeomorphism invariance of S_matter makes
T^(fields) covariantly conserved on shell, by Noether's second theorem.  Splitting it,

    T_mn = Pi_mn + (1/4) g_mn T,      T = g^{ab} T_ab,   Pi traceless,

and using nabla^m T_mn = 0, gives one line:

    nabla^m Pi_mn  =  -(1/4) nabla_n T                                     (EXACT)

so the traceless part is conserved exactly when the trace is covariantly constant.
Verified here on FRW with a massive scalar, on shell, component by component.

2.  AND THE COROLLARY IS THE CASE PST IS IN
--------------------------------------------
comp 158 showed the trace of Pi_mn VANISHES, for every configuration tried, because
the null threshold is the balance point of the pair-tension distribution.  T = 0 makes
the right-hand side zero identically.  So, conditional on the bridge, Pi_mn is
EXACTLY conserved on any background in the sector it actually describes, with no
(d_0/L)^2 residue at all.  The electromagnetic stress tensor is checked here to be
traceless for an ARBITRARY field strength in four dimensions, as the worked instance
of a source in that class.

    This is conditional on the bridge, which \\S sec:einstein flags as not
    established, and that condition is not hidden here.  What changes is the SHAPE of
    the open item: it is no longer "conservation is unproved", it is "conservation
    follows from an identification that is itself unproved".  comp 162's unconditional
    bound stands underneath it either way.

3.  M4: THE AMBIGUITY IS A SUPERPOTENTIAL, AND IT IS SMALLER THAN IT LOOKS
---------------------------------------------------------------------------
The freedom in a stress tensor assembled from pair contributions is the standard
Irving-Kirkwood/Hardy one: add the double divergence of a tensor with the algebraic
symmetries of Riemann.  Taking the simplest such object, H_{manb} = f (g_mn g_ab -
g_mb g_an), the shift is

    Delta_mn  =  g_mn box f  -  nabla_m nabla_n f

and three facts about it are computed, all exactly:

    (a)  it is SYMMETRIC, so it is an admissible regauging;
    (b)  its TRACE is 3 box f;
    (c)  its DIVERGENCE is exactly  -R_{ns} nabla^s f.

(c) is the M4 analogue of comp 162's result and has the same shape: zero in flat
space, curvature-suppressed otherwise.

4.  AND TRACELESSNESS CLOSES MOST OF IT
-----------------------------------------
By (b), keeping Pi traceless, which comp 158 FORCES rather than assumes, restricts
the gauge to box f = 0.  A harmonic gauge function that is localised is then zero:

    Riemannian section: zero Dirichlet data on the boundary of a region gives f = 0
    throughout, by uniqueness for the elliptic problem.  Checked discretely, and the
    operator's smallest singular value is reported so "unique" is measured rather
    than asserted.

    Lorentzian: zero Cauchy data on a spacelike slice gives f = 0 throughout the
    domain of dependence, by uniqueness for the wave equation.  Checked by explicit
    evolution with the data nonzero OUTSIDE the region, so the shrinking cone is a
    real test rather than a vacuous one.

So the ambiguity is EMPTY for any regauging confined to a region, and what survives
is a global harmonic function, which changes no integrated quantity.  M4 is not "the
tensor is undetermined"; it is "the tensor is determined up to a global harmonic
superpotential".

5.  THE UNIFYING STATEMENT
---------------------------
M3 and M4 are one phenomenon.  The pair construction is exactly conservative and
exactly gauge-rigid in flat space, and every curved-space violation is second order in
(pair scale)/(curvature radius), which is the same (d_0/L)^2 that bounds the
truncation of the projected field equations themselves (comp 154, comp 162).

A "pass": the trace identity exact on FRW; the EM trace zero for arbitrary F; the
shift symmetric with trace 3 box f and divergence exactly -R_{ns} nabla^s f; the
discrete Dirichlet solution identically zero with a nonzero smallest singular value;
and the Lorentzian cone of dependence staying at zero while the
solution outside it is order one.

WHAT THIS DOES NOT DO
---------------------
It does not establish the bridge.  It does not give Pi_mn an action of its own, and
the honest position is that it does not need one: what it needs is the identification,
and that is the open item.  The general superpotential is a full Riemann-symmetric
H_{manb}, of which only the one-scalar family is worked here; the conclusions about
trace and divergence are computed for that family and are not claimed for the general
one.
"""

import numpy as np
import sympy as sp

# ---------------------------------------------------------------------------
# a small covariant toolkit on a general metric
# ---------------------------------------------------------------------------
t, x, y, z = sp.symbols("t x y z", real=True)
X = [t, x, y, z]
a = sp.Function("a", positive=True)(t)
g = sp.diag(1, -a ** 2, -a ** 2, -a ** 2)      # FRW, spatially flat
ginv = g.inv()
N = 4


def rule(ch="="):
    print(ch * 100)


def christoffel():
    Ga = [[[0] * N for _ in range(N)] for _ in range(N)]
    for l in range(N):
        for m in range(N):
            for nu in range(N):
                Ga[l][m][nu] = sp.simplify(sum(
                    ginv[l, s] * (sp.diff(g[s, m], X[nu]) + sp.diff(g[s, nu], X[m])
                                  - sp.diff(g[m, nu], X[s]))
                    for s in range(N)) / 2)
    return Ga


GAMMA = christoffel()


def div2(T):
    """(nabla^m T_{mn})_n for a symmetric rank-2 T."""
    out = []
    for nu in range(N):
        tot = 0
        for m in range(N):
            for al in range(N):
                term = sp.diff(T[al, nu], X[m])
                term -= sum(GAMMA[s][m][al] * T[s, nu] for s in range(N))
                term -= sum(GAMMA[s][m][nu] * T[al, s] for s in range(N))
                tot += ginv[m, al] * term
        out.append(sp.simplify(tot))
    return out


def trace2(T):
    return sp.simplify(sum(ginv[i, j] * T[i, j] for i in range(N) for j in range(N)))


def hessian(f):
    df = [sp.diff(f, v) for v in X]
    H = sp.zeros(N, N)
    for i in range(N):
        for j in range(N):
            H[i, j] = sp.diff(f, X[i], X[j]) - sum(GAMMA[s][i][j] * df[s]
                                                   for s in range(N))
    return H, df


def ricci():
    Riem = [[[[0] * N for _ in range(N)] for _ in range(N)] for _ in range(N)]
    for r in range(N):
        for s in range(N):
            for m in range(N):
                for nu in range(N):
                    Riem[r][s][m][nu] = sp.simplify(
                        sp.diff(GAMMA[r][nu][s], X[m]) - sp.diff(GAMMA[r][m][s], X[nu])
                        + sum(GAMMA[r][m][l] * GAMMA[l][nu][s]
                              - GAMMA[r][nu][l] * GAMMA[l][m][s] for l in range(N)))
    Ric = sp.zeros(N, N)
    for s in range(N):
        for nu in range(N):
            Ric[s, nu] = sp.simplify(sum(Riem[r][s][r][nu] for r in range(N)))
    return Ric


# ---------------------------------------------------------------------------
# 1.  The exact trace identity
# ---------------------------------------------------------------------------
def part1_trace_identity():
    rule()
    print("  1.  nabla^m Pi_mn = -(1/4) nabla_n T, EXACTLY")
    rule()

    phi = sp.Function("phi")(t)
    m = sp.Symbol("m", positive=True)
    V = m ** 2 * phi ** 2 / 2
    dphi = [sp.diff(phi, v) for v in X]
    kin = sum(ginv[i, j] * dphi[i] * dphi[j] for i in range(N) for j in range(N))

    T = sp.zeros(N, N)
    for i in range(N):
        for j in range(N):
            T[i, j] = dphi[i] * dphi[j] - g[i, j] * (kin / 2 - V)

    boxphi = sum(ginv[i, j] * (sp.diff(phi, X[i], X[j])
                               - sum(GAMMA[s][i][j] * dphi[s] for s in range(N)))
                 for i in range(N) for j in range(N))
    eom = sp.simplify(boxphi + sp.diff(V, phi))
    phidd = sp.solve(sp.Eq(eom, 0), sp.diff(phi, t, 2))[0]
    print(f"      source: a massive scalar on FRW, eom  {eom} = 0")

    divT = [sp.simplify(d.subs(sp.diff(phi, t, 2), phidd)) for d in div2(T)]
    ok_T = all(d == 0 for d in divT)
    print(f"      nabla^m T_mn = 0 on shell                                 : {ok_T}")

    Tr = trace2(T)
    Pi = sp.zeros(N, N)
    for i in range(N):
        for j in range(N):
            Pi[i, j] = sp.simplify(T[i, j] - g[i, j] * Tr / 4)
    ok_tl = trace2(Pi) == 0
    print(f"      the traceless part Pi is traceless                        : {ok_tl}")
    print(f"      trace T = {Tr}, and it is NOT constant, so this is the hard case")

    divPi = [sp.simplify(d.subs(sp.diff(phi, t, 2), phidd)) for d in div2(Pi)]
    pred = [sp.simplify(sp.simplify(-sp.diff(Tr, v) / 4).subs(sp.diff(phi, t, 2), phidd))
            for v in X]
    ok_id = all(sp.simplify(p - q) == 0 for p, q in zip(divPi, pred))
    print(f"      nabla^m Pi_mn   = {divPi[0]}")
    print(f"      -(1/4) nabla_n T = {pred[0]}")
    print(f"      identity holds in all four components                     : {ok_id}")
    print()
    return bool(ok_T and ok_tl and ok_id)


# ---------------------------------------------------------------------------
# 2.  The corollary: a traceless source is exactly conserved
# ---------------------------------------------------------------------------
def part2_traceless_corollary():
    rule()
    print("  2.  A TRACELESS SOURCE IS CONSERVED EXACTLY, AND PST's SECTOR IS ONE")
    rule()
    print("      By part 1, T = 0 makes the right-hand side vanish identically, on any")
    print("      background.  comp 158 showed the trace of Pi_mn vanishes for every")
    print("      configuration tried, so this is the case PST is in.")
    print()
    print("      Worked instance: the electromagnetic stress tensor, for an ARBITRARY")
    print("      field strength in four dimensions.")

    F = sp.zeros(N, N)
    for i in range(N):
        for j in range(i + 1, N):
            c = sp.Function(f"F{i}{j}")(t, x, y, z)
            F[i, j] = c
            F[j, i] = -c
    Fup = ginv * F * ginv.T
    F2 = sp.simplify(sum(F[i, j] * Fup[i, j] for i in range(N) for j in range(N)))
    Tem = sp.zeros(N, N)
    for m in range(N):
        for n in range(N):
            Tem[m, n] = sp.simplify(
                sum(F[m, al] * ginv[al, be] * F[n, be]
                    for al in range(N) for be in range(N)) - g[m, n] * F2 / 4)
    tr = trace2(Tem)
    ok = sp.simplify(tr) == 0
    print(f"      trace of T^(EM) with six arbitrary component functions = {tr}")
    print(f"      identically zero, so nabla^m Pi_mn = 0 exactly            : {ok}")
    print()
    print("      Note what this changes and what it does not.  It does NOT establish")
    print("      the bridge of S sec:einstein.  It changes the shape of the open item:")
    print("      not 'conservation is unproved' but 'conservation follows from an")
    print("      identification that is itself unproved'.  comp 162's unconditional")
    print("      (d_0/L)^2 bound sits underneath either way.")
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# 3.  M4: the ambiguity is a superpotential
# ---------------------------------------------------------------------------
def part3_superpotential():
    rule()
    print("  3.  M4: THE AMBIGUITY IS A SUPERPOTENTIAL, EXACTLY CHARACTERISED")
    rule()
    print("      H_{manb} = f (g_mn g_ab - g_mb g_an) has the algebraic symmetries of")
    print("      Riemann, and the induced shift is Delta_mn = g_mn box f - nabla_m nabla_n f.")
    print()

    f = sp.Function("f")(t, x)
    H, df = hessian(f)
    boxf = sp.simplify(sum(ginv[i, j] * H[i, j] for i in range(N) for j in range(N)))
    Delta = sp.zeros(N, N)
    for i in range(N):
        for j in range(N):
            Delta[i, j] = sp.simplify(g[i, j] * boxf - H[i, j])

    ok_sym = sp.simplify(Delta - Delta.T) == sp.zeros(N, N)
    print(f"      (a) symmetric, so an admissible regauging                 : {ok_sym}")

    tr = trace2(Delta)
    ok_tr = sp.simplify(tr - 3 * boxf) == 0
    print(f"      (b) trace = 3 box f                                       : {ok_tr}")

    Ric = ricci()
    dd = div2(Delta)
    pred = [sp.simplify(-sum(Ric[nu, s] * ginv[s, k] * df[k]
                             for s in range(N) for k in range(N))) for nu in range(N)]
    ok_div = all(sp.simplify(p - q) == 0 for p, q in zip(dd, pred))
    print(f"      (c) divergence = -R_ns nabla^s f, exactly                 : {ok_div}")
    print(f"          so it is identically divergence-free in flat space, and")
    print(f"          curvature-suppressed otherwise: the same shape as comp 162.")
    print()
    return bool(ok_sym and ok_tr and ok_div)


# ---------------------------------------------------------------------------
# 4.  Tracelessness closes most of the ambiguity
# ---------------------------------------------------------------------------
def part4_locality():
    rule()
    print("  4.  TRACELESSNESS RESTRICTS THE GAUGE TO box f = 0, AND THEN LOCALITY")
    print("      KILLS IT")
    rule()
    print("      By 3(b), keeping Pi traceless, which comp 158 FORCES rather than")
    print("      assumes, requires box f = 0.  A localised harmonic gauge is zero.")
    print()

    # Riemannian section: elliptic Dirichlet problem
    n = 28
    idx = lambda i, j: i * n + j
    L = np.zeros((n * n, n * n))
    for i in range(n):
        for j in range(n):
            k = idx(i, j)
            if i in (0, n - 1) or j in (0, n - 1):
                L[k, k] = 1.0
            else:
                L[k, k] = -4.0
                for di, dj in ((1, 0), (-1, 0), (0, 1), (0, -1)):
                    L[k, idx(i + di, j + dj)] = 1.0
    sol = np.linalg.solve(L, np.zeros(n * n))
    smin = float(np.linalg.svd(L, compute_uv=False).min())
    ok_dir = np.abs(sol).max() < 1e-12 and smin > 1e-3
    print(f"      Riemannian, {n}x{n} discrete Dirichlet with zero boundary data:")
    print(f"        max |f| = {np.abs(sol).max():.3e}")
    print(f"        smallest singular value of the operator = {smin:.4f}, so the")
    print(f"        solution is unique and it is zero                       : {ok_dir}")

    # Lorentzian: zero Cauchy data does not merely stay zero, it stays zero for
    # the RIGHT reason, which is the domain of dependence.  Put nonzero data
    # everywhere OUTSIDE [0.4, 0.6] and check the shrinking cone inside it stays at
    # zero until the signal from the edge arrives, and not after.  A test that only
    # evolved zero data everywhere would be vacuous.
    nx = 2001
    dx = 1.0 / (nx - 1)
    cfl = 1.0            # the 1D wave scheme is exact at cfl = 1, so the numerical
    dt = cfl * dx        # signal speed equals the physical one and the cone is honest
    grid = np.linspace(0.0, 1.0, nx)
    u = np.where((grid < 0.4) | (grid > 0.6), np.sin(12 * np.pi * grid), 0.0)
    u_prev = u.copy()
    inside_max = 0.0
    n_steps = 150        # keeps the cone open: 150 * dx = 0.075 < 0.1
    for step in range(1, n_steps + 1):
        u_next = np.zeros(nx)
        u_next[1:-1] = (2 * u[1:-1] - u_prev[1:-1]
                        + cfl ** 2 * (u[2:] - 2 * u[1:-1] + u[:-2]))
        u_prev, u = u, u_next
        reach = step * dx                      # one cell per step at cfl = 1
        cone = (grid > 0.4 + reach + 2 * dx) & (grid < 0.6 - reach - 2 * dx)
        if cone.any():
            inside_max = max(inside_max, float(np.abs(u[cone]).max()))
    outside_amp = float(np.abs(u).max())
    ok_wave = inside_max < 1e-12 and outside_amp > 0.5
    print(f"      Lorentzian, wave equation, data nonzero only OUTSIDE [0.4, 0.6]:")
    print(f"        peak |f| inside the shrinking cone of dependence = {inside_max:.3e}")
    print(f"        peak |f| overall, so the test is not vacuous     = {outside_amp:.3f}")
    print(f"        zero data propagates as zero, by domain of dependence   : {ok_wave}")
    print()
    print("      So a regauging confined to a region is identically zero, and what")
    print("      survives is a GLOBAL harmonic function, which changes no integrated")
    print("      quantity.  M4 is not 'the tensor is undetermined'.")
    print()
    return bool(ok_dir and ok_wave)


def main():
    print()
    rule()
    print("  COMPUTATION 163: THE ACTION PRINCIPLE Pi_mn DOES HAVE, AND WHAT M4 IS")
    rule()
    print("  comp 162 said an exact result 'needs an action principle the construction")
    print("  does not have'.  It does not need one of its own: the identification in")
    print("  S sec:einstein hands it the diffeomorphism invariance of S_matter.")
    print()

    c1 = part1_trace_identity()
    c2 = part2_traceless_corollary()
    c3 = part3_superpotential()
    c4 = part4_locality()

    rule()
    print("  ASSESSMENT")
    rule()
    print(f"  1. nabla^m Pi_mn = -(1/4) nabla_n T, exact, verified on FRW     : {c1}")
    print(f"  2. a traceless source is conserved exactly; EM trace is zero    : {c2}")
    print(f"  3. the M4 shift: symmetric, trace 3 box f, div -R_ns grad^s f   : {c3}")
    print(f"  4. box f = 0 plus locality leaves nothing                       : {c4}")
    print()
    ok = all([c1, c2, c3, c4])
    print("  RESULT: " + (
        "CONSERVATION IS EXACT IN THE SECTOR Pi DESCRIBES, CONDITIONAL ON THE BRIDGE, "
        "because the traceless part of a conserved tensor obeys nabla^m Pi_mn = "
        "-(1/4) nabla_n T and comp 158 forces that trace to vanish.  M4's ambiguity is "
        "a superpotential whose trace is 3 box f, so tracelessness restricts it to "
        "harmonic gauges, and a localised harmonic gauge is identically zero.  What is "
        "open is the bridge, not the conservation and not the ambiguity"
        if ok else "MISMATCH -- revise"))
    print()


if __name__ == "__main__":
    main()
