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

Computation 199 -- the conformal part of the pair sum reaches exactly the traceless
tensors that obey the null energy condition, so the identification has no obstruction
===================================================================================
WHAT THIS DOES.  comp 164 found the one obstruction to the conformal identification of
sec:einstein: every term of eq:Pimunu is a single-signed weight times a dyad of the
pair separation, so the pair sum is semidefinite, while a static magnetic or electric
field is metric-traceless and indefinite (tension along the field lines, pressure
across them).  No pair sum equals it, and the only rescue was a universal constant
shift that nothing fixes.  That obstruction has stood since v32.74 as the first
structural item.

THE MOVE.  The identification is a statement about CONFORMAL PARTS, and the book
defines the conformal part of a tensor as the tensor minus a quarter of its metric
trace times the metric.  Apply that definition to the pair sum itself.  Term by term,

        Delta_m Delta_n  ->  Delta_m Delta_n - (s/4) g_mn,        s = Delta . Delta,

which is traceless identically, with no balance condition, and which is INDEFINITE:
a spacelike bond now carries pressure along itself and tension across itself.  Two
families of bonds lying in the plane transverse to a direction sum to exactly the
magnetic stress tensor with the field along that direction, and the electric tensor
is the same construction in the other transverse plane.  So the conformal part of the
pair sum reaches what comp 164 showed the pair sum cannot, and nothing in eq:Pimunu
changed: the conformal part of Pi_mn is Pi_mn - (1/4) g_mn tr Pi, and comp 158 showed
tr Pi = 0 in balance, so on every balanced cloud, which is every result the book rests
on here, the two objects coincide.

THE THEOREM.  Let R be the closed cone generated by the traceless dyads
Delta Delta - (s/4) g over all Delta, with the positive coefficients the forced
weight W(s) = -|s| and the bond kernel supply.  Then

        R  =  { T symmetric, metric-traceless :  T(k, k) >= 0 for every null k },

the metric-traceless tensors obeying the NULL ENERGY CONDITION.

  Inclusion.  For null k the metric term drops, so (Delta Delta - (s/4) g)(k, k)
  = (Delta . k)^2 >= 0: every generator obeys the condition, hence so does the cone.

  Equality.  The dual of R, paired against traceless tensors, is the set of
  contravariant matrices Y that are positive semidefinite and metric-traceless (take Y
  traceless without loss, then <generator, Y> = Y(Delta, Delta) >= 0 for all Delta is
  semidefiniteness).  That set is a spectrahedron, the semidefinite cone cut by one
  linear equation, and its normalised slice has one more; by Pataki's bound an extreme
  point of rank r of a spectrahedron with m equations has r(r+1)/2 <= m, so with m = 2
  every extreme ray has rank one, Y = k k with g(k, k) = 0, a null dyad.  The slice is
  compact, so every Y is a positive combination of null dyads, the dual of R is the
  cone of null dyads, and R, being closed, is the dual of that: the null energy
  condition on traceless tensors.

WHAT THIS CLOSES.  Every classical field obeys the null energy condition, radiation,
null waves, static electric and magnetic fields among them, so the conformal part of
the field stress tensor is reachable by the conformal part of the pair sum for every
field this theory contains.  The obstruction comp 164 found is gone, and with comp 196
having shown the magnitude is the state's, the identification is left with no gap: it
is an identification, of the kind Appendix C.6 lists, and no longer a derivation the
structural front owes.  What is excluded is exactly what a null-energy violation is:
a traceless tensor with negative energy density along some null direction, such as a
pure anisotropic stress with zero energy density, and nothing the book contains is of
that kind.

WHAT THIS COSTS, STATED RATHER THAN HIDDEN.  Out of balance the pair sum has a trace,
and by comp 158 that trace is minus the summed pair tension relative to the null
threshold, the balance scalar.  The identification is therefore of conformal parts:
the conformal part of Pi_mn is the conformal part of T^(fields)_mn, and the trace of
an unbalanced pair sum is a scalar the identification does not use.  comp 159 had
already met that scalar: the worldline sector of the pair sum is exactly kinetic,
w = beta^2 / 3, dust at rest, and carries a few per cent of the density for
homogeneous matter, which is the share by which the projection there is conformal
rather than exactly traceless.  comp 163's exact conservation is the balanced case, where the
trace is constant; eq:trace-conservation is the general statement and it is an identity
for the conformal part of any conserved tensor.  And the forced weight is load-bearing
in a new way: a weight that changed sign across the light cone would flip the
spacelike generators and the cone would reach null-energy violations, so comp 158's
W(s) = -|s| is what confines the pair sum to physical matter.

CHECKS.  (1) The per-term trace vanishes identically (symbolic); the conformal part of
the sum is the sum of conformal parts; the trace of the sum is sum of |s| s.
(2) Explicit clouds: pairs on concentric regular polygons in the x-y plane give
diag(u, u, u, -u) to machine precision after the conformal projection (magnetic, field
along z), the same in the y-z plane gives the electric tensor along x; in comp 158's
balanced regime (uniform cloud, hard cutoff, growing window) the trace is below a part
in a thousand of the scale and the projection changes the pair sum by that much, so
comps 158 to 161 are untouched; comp 164's six shapes are tabulated by their trace
fraction, balanced or not, and every one obeys the null energy condition.  (3) Inclusion,
symbolic: (Delta . k)^2 for generic Delta and null k.  (4) Equality, numerical: random
positive semidefinite traceless Y decompose into sampled null dyads with a residual
that falls with the sampling density; random traceless tensors obeying the null
energy condition decompose into sampled generators with a residual that falls the same
way, and random traceless tensors violating it do not, the residual staying above a
fixed fraction of the violation, a null dyad being the certificate.  The rank-two
case of the dual, where Pataki's bound bites, is done by hand: such a Y spans a
Lorentzian plane and equals a non-negative combination of the dyads of that plane's
two null directions, exactly, on 200 random draws.  (5) The magnetic tensor sits on
the boundary of the cone, the condition saturating along the field lines, and is
reached exactly.  (6) Tamper: a signed weight W = -s reaches the zero-energy
anisotropic stress that the forced weight excludes, and gives a cloud strung along
the time axis a negative energy density.
"""

import numpy as np
import sympy as sp
from scipy.optimize import linprog

ETA = np.diag([1.0, -1.0, -1.0, -1.0])


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


# ---------------------------------------------------------------------------
# the pair sum, eq:Pimunu with comp 164's conventions, and its conformal part
# ---------------------------------------------------------------------------
def pi_of(pts, cut=4.0, signed=False):
    """Covariant Pi_mn = sum |s| Delta_m Delta_n over pairs with |s| <= cut.
    signed=True is Part 6's tamper: weight s instead of |s|."""
    P = np.zeros((4, 4))
    for i in range(len(pts)):
        D = pts[i] - pts[i + 1:]
        if len(D) == 0:
            continue
        s = np.einsum("pi,ij,pj->p", D, ETA, D)
        keep = np.abs(s) <= cut
        D, s = D[keep], s[keep]
        if len(s) == 0:
            continue
        Dl = D @ ETA
        w = -s if signed else np.abs(s)
        P += np.einsum("p,pm,pn->mn", w, Dl, Dl)
    return P


def g_trace(T):
    """metric trace of a covariant tensor."""
    return float(np.einsum("mn,mn->", np.linalg.inv(ETA), T))


def conformal(T):
    return T - 0.25 * g_trace(T) * ETA


def nec_min(T, n=4000, seed=0):
    """min over null k = (1, n) of T(k, k), n on the unit sphere (covariant T, contravariant k)."""
    rng = np.random.default_rng(seed)
    v = rng.normal(size=(n, 3))
    v /= np.linalg.norm(v, axis=1)[:, None]
    k = np.column_stack([np.ones(n), v])
    return float(np.einsum("pi,ij,pj->p", k, T, k).min())


# ---------------------------------------------------------------------------
# Part 1: the per-term projection
# ---------------------------------------------------------------------------
def part1():
    rule()
    print("  1.  THE CONFORMAL PART, TERM BY TERM")
    rule()
    d0, d1, d2, d3 = sp.symbols("d0 d1 d2 d3", real=True)
    g = sp.diag(1, -1, -1, -1)
    D = sp.Matrix([d0, d1, d2, d3])
    Dl = g * D
    s = (D.T * g * D)[0]
    G = Dl * Dl.T - s / 4 * g
    tr = sp.simplify((g.inv() * G).trace())
    ok1 = tr == 0
    print(f"      tr_g [Delta Delta - (s/4) g] = {tr}, identically, for every Delta            : {ok1}")
    # mixed components in the rest frame of a timelike Delta and along a spacelike one
    Gt = (g.inv() * G.subs({d0: 1, d1: 0, d2: 0, d3: 0}))
    Gs = (g.inv() * G.subs({d0: 0, d1: 0, d2: 0, d3: 1}))
    print(f"      timelike unit Delta, mixed components T^m_n = diag{tuple(Gt[i, i] for i in range(4))}: rho = 3/4, p = 1/4 isotropic, radiation-like")
    print(f"      spacelike unit Delta along z, mixed            = diag{tuple(Gs[i, i] for i in range(4))}: rho = 1/4, tension 1/4 across, pressure 3/4 along")
    ok2 = tuple(Gs[i, i] for i in range(4)) == (sp.Rational(1, 4), sp.Rational(1, 4), sp.Rational(1, 4), sp.Rational(-3, 4))
    # linearity: conformal part of the sum is the sum of conformal parts
    rng = np.random.default_rng(1)
    pts = rng.normal(size=(60, 4))
    P = pi_of(pts)
    # recompute the sum of per-term conformal parts
    Q = np.zeros((4, 4))
    for i in range(len(pts)):
        Dm = pts[i] - pts[i + 1:]
        sv = np.einsum("pi,ij,pj->p", Dm, ETA, Dm)
        keep = np.abs(sv) <= 4.0
        Dm, sv = Dm[keep], sv[keep]
        Dl = Dm @ ETA
        Q += np.einsum("p,pm,pn->mn", np.abs(sv), Dl, Dl) - 0.25 * np.sum(np.abs(sv) * sv) * ETA
    ok3 = np.allclose(conformal(P), Q, atol=1e-9 * abs(P).max())
    print(f"      conformal part of the sum = sum of per-term conformal parts (60 points)       : {ok3}")
    ok4 = abs(g_trace(P) - sum(abs(sv) * sv for i in range(len(pts)) for sv in [np.einsum('pi,ij,pj->p', pts[i] - pts[i + 1:], ETA, pts[i] - pts[i + 1:])] for sv in sv[np.abs(sv) <= 4.0])) < 1e-9 * abs(P).max()
    print(f"      tr_g Pi = sum |s| s, minus comp 158's summed (T - tau_null): the balance scalar : {ok4}")
    print()
    return ok1 and ok2 and ok3 and ok4


# ---------------------------------------------------------------------------
# Part 2: explicit clouds reach the magnetic and electric tensors
# ---------------------------------------------------------------------------
def planar_cloud(rng, n, axes):
    pts = np.zeros((n, 4))
    pts[:, axes] = rng.normal(size=(n, 2))
    return pts


def ring_cloud(axes, rings=(0.5, 1.0, 1.5, 2.0), per_ring=24):
    """points on concentric regular polygons in one spatial plane: isotropic in the plane by symmetry."""
    pts = []
    for r in rings:
        for j in range(per_ring):
            a = 2 * np.pi * j / per_ring
            q = np.zeros(4)
            q[axes[0]], q[axes[1]] = r * np.cos(a), r * np.sin(a)
            pts.append(q)
    return np.array(pts)


def uniform_cloud(rng, n, R):
    """comp 158's cloud: t uniform in [-R, R], space uniform in the ball of radius R."""
    t = rng.uniform(-R, R, n)
    v = rng.normal(size=(n, 3))
    v /= np.linalg.norm(v, axis=1)[:, None]
    d = rng.uniform(0, 1, n) ** (1 / 3) * R
    return np.column_stack([t, v * d[:, None]])


def part2():
    rule()
    print("  2.  THE MAGNETIC AND ELECTRIC TENSORS, FROM EXPLICIT PAIR CLOUDS")
    rule()
    rng = np.random.default_rng(20261010)
    ok = True
    for name, axes, pattern in (("magnetic, field along z: pairs in the x-y plane", (1, 2), np.array([1, 1, 1, -1.0])),
                                ("electric, field along x: pairs in the y-z plane", (2, 3), np.array([1, -1, 1, 1.0]))):
        pts = ring_cloud(axes)
        P = pi_of(pts, cut=1e9)
        C = conformal(P)
        u = C[0, 0]
        target = u * np.diag(pattern)
        off = np.abs(C - np.diag(np.diag(C))).max() / u
        diag_dev = np.abs(np.diag(C) - np.diag(target)).max() / u
        print(f"      {name} ({len(pts)} points on concentric regular polygons)")
        print(f"        conformal part / u = diag{tuple(np.round(np.diag(C) / u, 6).tolist())},  target diag{tuple(pattern.tolist())}; max off-diagonal {off:.1e}, max diagonal deviation {diag_dev:.1e}")
        print(f"        trace of the raw sum / u = {g_trace(P) / u:+.3f} (unbalanced: every pair spacelike); raw sum min eigenvalue / u = {np.linalg.eigvalsh(P).min() / u:+.3f} (semidefinite, comp 164)")
        ok &= diag_dev < 1e-9 and off < 1e-9 and np.linalg.eigvalsh(C).min() < -0.9 * u
    # exact two-bond construction
    ex = conformal(np.diag([0, 1, 0, 0.0])) + conformal(np.diag([0, 0, 1, 0.0]))
    ok_ex = np.allclose(ex, 0.5 * np.diag([1, 1, 1, -1.0]))
    print(f"      one x-bond plus one y-bond, exactly: {np.diag(ex)} = (1/2) diag(1, 1, 1, -1)                       : {ok_ex}")
    # balanced cloud, comp 158's regime: uniform cloud, hard cutoff |s| <= 0.5, window growing
    print("      comp 158's balanced regime (uniform cloud, 1500 points, hard cutoff |s| <= 0.5): the trace fraction and the change the projection makes")
    frac = rel = None
    for R in (2.0, 4.0, 8.0, 16.0):
        P = pi_of(uniform_cloud(rng, 1500, R), cut=0.5)
        frac = abs(g_trace(P)) / abs(P).max()
        rel = np.abs(conformal(P) - P).max() / abs(P).max()
        print(f"        window {R:5.1f}: |tr Pi| / scale = {frac:.1e},  |conformal part - Pi| / scale = {rel:.1e}")
    print("        so on a balanced cloud the conformal part is the pair sum itself, and every result of comps 158 to 161 is untouched")
    ok &= frac < 2e-3 and rel < 2e-3
    # comp 164's six shapes, classified
    shapes = [
        ("isotropic gaussian", lambda: rng.normal(size=(250, 4)) * 1.5),
        ("flattened to a plane", lambda: rng.normal(size=(250, 4)) * np.array([1.5, 1.5, 1.5, 0.05])),
        ("stretched along time", lambda: rng.normal(size=(250, 4)) * np.array([6.0, 0.4, 0.4, 0.4])),
        ("two clumps", lambda: np.vstack([rng.normal(size=(125, 4)) * 0.4, rng.normal(size=(125, 4)) * 0.4 + np.array([0, 3, 0, 0])])),
        ("laid on a light cone", lambda: (lambda v: np.column_stack([np.linalg.norm(v, axis=1), v]))(rng.normal(size=(250, 3)))),
        ("collapsed onto a line", lambda: np.column_stack([rng.normal(size=250) * 0.02, np.linspace(-3, 3, 250), rng.normal(size=250) * 0.02, rng.normal(size=250) * 0.02])),
    ]
    print("      comp 164's six shapes: trace fraction of the raw sum, and the minimum eigenvalue of the conformal part")
    for name, gen in shapes:
        P = pi_of(gen())
        C = conformal(P)
        print(f"        {name:22s} tr/scale {g_trace(P) / abs(P).max():+.3f}   min eig(conformal)/scale {np.linalg.eigvalsh(C).min() / abs(C).max():+.3f}   NEC min {nec_min(C) / abs(C).max():+.3f}")
    print("      balanced shapes keep the projection trivial; unbalanced ones make it indefinite, and every one obeys the null energy condition")
    print()
    return bool(ok and ok_ex)


# ---------------------------------------------------------------------------
# Part 3: inclusion, symbolic
# ---------------------------------------------------------------------------
def part3():
    rule()
    print("  3.  INCLUSION: EVERY GENERATOR OBEYS THE NULL ENERGY CONDITION")
    rule()
    d0, d1, d2, d3, th, ph = sp.symbols("d0 d1 d2 d3 theta phi", real=True)
    g = sp.diag(1, -1, -1, -1)
    D = sp.Matrix([d0, d1, d2, d3])
    Dl = g * D
    s = (D.T * g * D)[0]
    G = Dl * Dl.T - s / 4 * g
    k = sp.Matrix([1, sp.sin(th) * sp.cos(ph), sp.sin(th) * sp.sin(ph), sp.cos(th)])
    null = sp.simplify((k.T * g * k)[0])
    val = sp.simplify((k.T * G * k)[0] - ((Dl.T * k)[0]) ** 2)
    ok = null == 0 and val == 0
    print(f"      k = (1, n) is null: k.k = {null};  G(k, k) - (Delta . k)^2 = {val}, so G(k, k) = (Delta . k)^2 >= 0        : {ok}")
    print("      hence the cone the generators span lies inside the null-energy traceless tensors")
    print()
    return ok


# ---------------------------------------------------------------------------
# Part 4: equality, numerically
# ---------------------------------------------------------------------------
def sym_basis():
    """basis of symmetric 4x4 matrices as vectors (10 components)."""
    idx = [(i, j) for i in range(4) for j in range(i, 4)]
    return idx


IDX = sym_basis()


def vec(T):
    return np.array([T[i, j] for i, j in IDX])


def lp_residual(targets, generators):
    """min L1 residual of a nonnegative combination of generator vectors matching a target."""
    Gm = np.array([vec(G) for G in generators]).T      # 10 x N
    out = []
    for T in targets:
        t = vec(T)
        N = Gm.shape[1]
        # variables: w (N) >= 0, r+ (10) >= 0, r- (10) >= 0 ; Gm w + r+ - r- = t ; min sum(r+ + r-)
        c = np.concatenate([np.zeros(N), np.ones(20)])
        A = np.hstack([Gm, np.eye(10), -np.eye(10)])
        res = linprog(c, A_eq=A, b_eq=t, bounds=[(0, None)] * (N + 20), method="highs")
        out.append(res.fun / np.abs(t).max())
    return np.array(out)


def null_dyads(n, rng):
    v = rng.normal(size=(n, 3))
    v /= np.linalg.norm(v, axis=1)[:, None]
    k = np.column_stack([np.ones(n), v])
    return [np.outer(a, a) for a in k]


def generators(n, rng, rapidity=2.0, signed=False):
    """traceless dyads for unit timelike and spacelike Delta, boosts up to the given rapidity."""
    out = []
    for _ in range(n // 2):
        v = rng.normal(size=3)
        v /= np.linalg.norm(v)
        chi = rng.uniform(0, rapidity)
        Dt = np.array([np.cosh(chi), *(np.sinh(chi) * v)])          # unit timelike, s = +1
        w = rng.normal(size=3)
        w -= v * (w @ v)
        w /= np.linalg.norm(w)
        Ds = np.array([np.sinh(chi), *(np.cosh(chi) * w)]) if rng.uniform() < 0.5 else np.array([0, *w])   # unit spacelike, s = -1
        for Dv, sgn in ((Dt, 1.0), (Ds, -1.0)):
            Dl = ETA @ Dv
            G = np.outer(Dl, Dl) - 0.25 * sgn * ETA
            out.append(-G if (signed and sgn < 0) else G)
    return out


def random_psd_traceless(rng):
    A = rng.normal(size=(4, 4))
    Z = A @ A.T
    d = Z[1, 1] + Z[2, 2] + Z[3, 3] - Z[0, 0]
    if d > 0:
        Z[0, 0] += d
    else:
        Z[1:, 1:] += (-d / 3) * np.eye(3)
    return Z / abs(Z).max()


def random_traceless(rng):
    A = rng.normal(size=(4, 4))
    T = (A + A.T) / 2
    return conformal(T)


def part4():
    rule()
    print("  4.  EQUALITY: THE DUAL CONE IS THE NULL DYADS, AND THE PRIMAL CONE IS THE NULL-ENERGY TENSORS")
    rule()
    rng = np.random.default_rng(7)
    ok = True
    # (a) dual: random PSD traceless Y into null dyads, residual falls with the sampling density
    Ys = [random_psd_traceless(rng) for _ in range(20)]
    print("      (a) positive semidefinite, metric-traceless Y as positive combinations of sampled null dyads")
    prev = None
    res_a = []
    for n in (50, 200, 800):
        r = lp_residual(Ys, null_dyads(n, rng))
        res_a.append(r.max())
        print(f"          {n:4d} null directions: worst relative L1 residual over 20 random Y = {r.max():.2e}")
    ok &= res_a[-1] < 2e-2 and res_a[-1] < res_a[0]
    # (b) primal: traceless tensors classified by the null energy condition
    Ts, cls = [], []
    e00 = np.zeros((4, 4)); e00[0, 0] = 1.0
    for j in range(40):
        Z = random_traceless(rng)
        m = nec_min(Z)                      # shifting the energy density by lam moves T(k, k) by exactly lam
        lam = -m + (0.15 if j % 2 == 0 else -0.6) * abs(Z).max()
        T = conformal(Z + lam * e00)
        T /= abs(T).max()
        Ts.append(T); cls.append(nec_min(T) > 0.05)
    cls = np.array(cls)
    viol = np.array([-nec_min(T) for T in Ts])        # the size of the violation, zero or negative on the obeying set
    print(f"      (b) random traceless tensors: {cls.sum()} obey the null energy condition with margin, {(~cls).sum()} violate it with margin")
    res_b = []
    for n in (100, 400, 1600):
        r = lp_residual(Ts, generators(n, rng))
        ratio = (r[~cls] / viol[~cls]).min()
        res_b.append((r[cls].max(), r[~cls].min(), ratio))
        print(f"          {n:4d} generators: worst residual on the obeying set {r[cls].max():.2e};  on the violating set the residual is at least {ratio:.2f} times the violation (min residual {r[~cls].min():.2e})")
    print("          a null dyad k k is a certificate: the cone pairs non-negatively with it while a violator pairs to T(k, k) < 0, so the residual cannot fall below a fixed fraction of the violation")
    ok &= res_b[-1][0] < 3e-2 and res_b[-1][0] < res_b[0][0] and res_b[-1][2] > 0.1
    # (c) the two named tensors
    mag = np.diag([1, 1, 1, -1.0])
    ele = np.diag([1, -1, 1, 1.0])
    stress = np.diag([0, 1, -1, 0.0])       # zero energy density, anisotropic stress: traceless, violates the condition
    neg = -mag
    r = lp_residual([mag, ele, stress, neg], generators(1600, rng))
    print(f"      (c) magnetic {r[0]:.2e}, electric {r[1]:.2e} (reached, on the boundary: NEC min {nec_min(mag):+.2f} along the field);  zero-energy anisotropic stress {r[2]:.2e}, negative-energy radiation {r[3]:.2e} (not reached; NEC min {nec_min(stress):+.2f}, {nec_min(neg):+.2f})")
    ok &= r[0] < 3e-2 and r[1] < 3e-2 and r[2] > 0.1 and r[3] > 0.1
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# Part 5: Pataki's bound, checked on the slice
# ---------------------------------------------------------------------------
def part5():
    rule()
    print("  5.  THE EXTREME RAYS OF THE DUAL CONE ARE RANK ONE, HENCE NULL DYADS")
    rule()
    # Pataki: extreme points of {X >= 0, <A_i, X> = b_i, i = 1..m} have rank r with r(r+1)/2 <= m.
    # Here m = 2 (metric trace zero, Euclidean trace one), so r = 1: X = w w^T, and metric trace zero makes w null.
    m = 2
    ranks = [r for r in range(1, 5) if r * (r + 1) // 2 <= m]
    print(f"      spectrahedron with m = {m} equations: admissible extreme ranks {ranks}")
    # numerical witness: maximise random linear functionals over the slice (an SDP solved by sampling + LP over null dyads
    # would be circular); instead check that random rank-2 PSD traceless Y are never extreme: they split into null dyads
    rng = np.random.default_rng(11)
    ok = ranks == [1]
    # the rank-two case by hand, which is where the bound bites: Y = a a + b b with g(a, a) + g(b, b) = 0 spans a plane
    # whose induced metric has discriminant g(a, b)^2 - g(a, a) g(b, b) = g(a, b)^2 + g(a, a)^2 >= 0, so the plane is
    # Lorentzian and holds two null directions k1, k2; g-tracelessness kills the cross term k1 k2 + k2 k1 because
    # g(k1, k2) is not zero, and positivity makes the two coefficients non-negative: Y = c1 k1 k1 + c2 k2 k2 exactly.
    worst = 0.0
    cmin = np.inf
    for _ in range(200):
        a, b = rng.normal(size=4), rng.normal(size=4)
        f = np.sqrt((a[1:] @ a[1:] + b[1:] @ b[1:]) / (a[0] ** 2 + b[0] ** 2))
        a[0] *= f; b[0] *= f
        Y = np.outer(a, a) + np.outer(b, b)
        gaa, gab, gbb = a @ ETA @ a, a @ ETA @ b, b @ ETA @ b
        disc = gab ** 2 - gaa * gbb
        roots = [(-gab + sg * np.sqrt(disc)) / gaa if abs(gaa) > 1e-12 else None for sg in (1, -1)]
        ks = [a * r + b for r in roots]                    # the two null directions of the plane
        M = np.array([[np.outer(k, k)[i, j] for k in ks] for i, j in IDX])
        c, *_ = np.linalg.lstsq(M, vec(Y), rcond=None)
        resid = np.abs(M @ c - vec(Y)).max() / abs(Y).max()
        worst = max(worst, resid)
        cmin = min(cmin, c.min())
    print(f"      200 random rank-two members: each equals c1 k1 k1 + c2 k2 k2 on the two null directions of its plane, worst residual {worst:.1e}, smallest coefficient {cmin:+.2e}   : {worst < 1e-9 and cmin > -1e-9}")
    ok &= worst < 1e-9 and cmin > -1e-9
    print()
    return bool(ok)


# ---------------------------------------------------------------------------
# Part 6: tamper, the signed weight
# ---------------------------------------------------------------------------
def part6():
    rule()
    print("  6.  TAMPER: A SIGNED WEIGHT W = -s REACHES NULL-ENERGY VIOLATIONS")
    rule()
    rng = np.random.default_rng(5)
    stress = np.diag([0, 1, -1, 0.0])
    r_forced = lp_residual([stress], generators(1600, rng))[0]
    r_signed = lp_residual([stress], generators(1600, rng, signed=True))[0]
    print(f"      zero-energy anisotropic stress: residual with the forced weight {r_forced:.2e}, with the signed weight {r_signed:.2e}")
    # and on a cloud: with W = -s the timelike pairs carry a negative weight, so a cloud strung along the time axis
    # has a negative energy density after the projection
    pts = np.column_stack([np.linspace(-2, 2, 200), rng.normal(size=(200, 3)) * 0.02])
    C = conformal(pi_of(pts, cut=1e9, signed=True))
    print(f"      a cloud along the time axis under the signed weight: energy density / scale = {C[0, 0] / abs(C).max():+.3f}, NEC min / scale = {nec_min(C) / abs(C).max():+.3f}")
    ok = r_forced > 0.1 and r_signed < 3e-2 and C[0, 0] < 0
    print(f"      so comp 158's single-signed weight is what confines the pair sum to physical matter             : {ok}")
    print()
    return bool(ok)


def main():
    rule()
    print("  COMPUTATION 199.  The conformal part of the pair sum reaches exactly the null-energy traceless tensors")
    rule()
    results = [
        ("per-term conformal projection: traceless identically, linear, trace of the sum is the balance scalar", part1()),
        ("explicit clouds reach the magnetic and electric tensors; balanced clouds are unchanged", part2()),
        ("inclusion: every generator obeys the null energy condition (symbolic)", part3()),
        ("equality: dual cone is the null dyads, primal cone is the null-energy traceless tensors (numerical)", part4()),
        ("Pataki: extreme rays of the dual are rank one, hence null dyads", part5()),
        ("tamper: a signed weight reaches null-energy violations, so the forced weight is load-bearing", part6()),
    ]
    rule("-")
    for i, (name, r) in enumerate(results, 1):
        print(f"  {i}. {name}: {'PASS' if r else 'FAIL'}")
    allok = all(r for _, r in results)
    rule("-")
    print("  VERDICT:", "PASS" if allok else "FAIL")
    print("""
  The conformal part of the pair sum reaches exactly the metric-traceless tensors that
  obey the null energy condition.  Static electric and magnetic fields are reached, by
  bonds in the transverse plane, and nothing a classical field produces is excluded.
  The obstruction comp 164 found is gone; with comp 196 the magnitude is the state's;
  the conformal identification has no gap left and leaves the structural front as an
  identification.  The cost is a reading: the identification is of conformal parts, and
  an unbalanced pair sum carries a trace that the identification does not use.""")
    return allok


if __name__ == "__main__":
    main()