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

Computation 131 -- the A6 weak (H^1-tested) discrepancy is BLIND to the null-cone
vanishing-separation set (codim >= 1); only the codim-0 bulk node-density binds (open_research §1 A6)
=================================================================================
Backs the A2-core characterization (cycle 7, 27/27 survived). The A6 bias coefficient
    c_u = |D_n|^{-1} sum_a |grad u(rho_n(a))|^2 / (|K|^{-1} int_K |grad u|^2)
is a sum over SINGLE nodes against the empirical node measure. The Lorentzian directed
metric the signed-tension metric makes near-null pairs (s^2 -> 0) -- vanishing STRONG (pair) separation --
but the vanishing-separation locus is a codimension >= 1 relation on pair-space D x D, so it
does not perturb the single-node measure: c_u is blind to it. The strong (separation/star)
channel therefore does NOT bind the weak discrepancy; the residual is purely the codim-0
node-density channel (a positive-volume void/excess), which the matched count A1 does not force.

This script realises the contrast on a 3-cube with the Lorentzian form s^2 = -dx0^2+dx1^2+dx2^2:

  (1) STRONG min-separation collapses to 0 once null pairs are present -- in BOTH a uniform
      and a void cloud -- since near-null pairs are generic (codim >= 1, but the MINIMUM hits 0).
  (2) The WEAK defect |c_u - 1| is UNCHANGED by injecting null pairs at bulk-typical positions
      (the paired nodes sit at volume-typical locations, so the single-node average is unmoved):
      a uniform cloud keeps c_u -> 1 even with strong separation driven to 0.
  (3) Only a codim-0 VOID biases c_u: |c_u - 1| stays ~0 for the void-free cloud and grows with
      the void volume, independent of the null-cone pairs.

A "pass": (1) min|s| -> ~0 with injected null pairs; (2) |c_u-1| essentially identical with and
without the null-pair injection (weak survives the null cone); (3) |c_u-1| ~ 0 uniform, large for
the void -- the codim-0 density is the sole binding channel.
"""

from __future__ import annotations
import numpy as np

CONT = 3 * np.pi ** 2 / 8.0   # (1/|K|) int_{[0,1]^3} |grad sin(pi x)sin(pi y)sin(pi z)|^2


def gradf3(P):
    s = np.sin(np.pi * P)
    c = np.cos(np.pi * P)
    sx, sy, sz = s[:, 0], s[:, 1], s[:, 2]
    cx, cy, cz = c[:, 0], c[:, 1], c[:, 2]
    return np.pi ** 2 * (cx ** 2 * sy ** 2 * sz ** 2
                         + sx ** 2 * cy ** 2 * sz ** 2
                         + sx ** 2 * sy ** 2 * cz ** 2)


def c_u(P):
    return gradf3(P).mean() / CONT


def lorentz_min_sep(P, rng, sample=400000):
    """min over sampled pairs of |s|, s^2 = -dx0^2 + dx1^2 + dx2^2 (signature -,+,+)."""
    N = len(P)
    i = rng.integers(0, N, sample)
    j = rng.integers(0, N, sample)
    m = i != j
    d = P[i] - P[j]
    s2 = -d[:, 0] ** 2 + d[:, 1] ** 2 + d[:, 2] ** 2
    return np.sqrt(np.abs(s2[m]).min())


def jittered_grid3(n, rng):
    xs = (np.arange(n) + 0.5) / n
    X, Y, Z = np.meshgrid(xs, xs, xs)
    P = np.column_stack([X.ravel(), Y.ravel(), Z.ravel()])
    P += rng.uniform(-0.3, 0.3, P.shape) / n
    return np.clip(P, 1e-6, 1 - 1e-6)


def void3(n, r, rng):
    P = jittered_grid3(n, rng)
    d = np.sum((P - 0.5) ** 2, axis=1)
    return P[d > r * r]


def inject_null_pairs(P, K, delta, rng):
    """add K exact null pairs: a uniform, partner a+(delta,delta,0) => s^2 = -delta^2+delta^2 = 0."""
    a = rng.uniform(0.1, 0.9, (K, 3))
    b = np.clip(a + np.array([delta, delta, 0.0]), 1e-6, 1 - 1e-6)
    return np.vstack([P, a, b])


def main():
    print("=" * 100)
    print("  Computation 131 -- weak A6 discrepancy is blind to the null cone (codim>=1); only codim-0 density binds")
    print("=" * 100)
    print()
    rng = np.random.default_rng(20260628)
    n = 18

    uni = jittered_grid3(n, rng)
    uni_null = inject_null_pairs(uni, 60, 0.02, rng)
    vd = void3(n, 0.22, rng)

    print("  (1) STRONG Lorentzian min-separation min|s|, s^2 = -dx0^2+dx1^2+dx2^2:")
    m_uni = lorentz_min_sep(uni, rng)
    m_uni_null = lorentz_min_sep(uni_null, rng)
    m_void = lorentz_min_sep(vd, rng)
    print(f"      uniform                 : min|s| = {m_uni:.5f}")
    print(f"      uniform + 60 null pairs : min|s| = {m_uni_null:.5e}   (collapsed by the null cone)")
    print(f"      void                    : min|s| = {m_void:.5f}")
    p1 = m_uni_null < 1e-6
    print(f"      -> injected null pairs drive strong separation to ~0: {p1}")
    print()

    print("  (2) WEAK defect |c_u - 1| is UNCHANGED by the null-pair injection (weak survives the null cone):")
    d_uni = abs(c_u(uni) - 1.0)
    d_uni_null = abs(c_u(uni_null) - 1.0)
    print(f"      uniform                 : |c_u-1| = {d_uni:.5f}")
    print(f"      uniform + 60 null pairs : |c_u-1| = {d_uni_null:.5f}   (strong sep = 0, yet weak unmoved)")
    p2 = abs(d_uni_null - d_uni) < 0.01 and d_uni_null < 0.02
    print(f"      -> the codim>=1 null-cone pairs do not perturb the single-node measure: {p2}")
    print()

    print("  (3) only the codim-0 VOID biases c_u (the sole binding channel; cf. the stronger 2D void in Comp 130):")
    print(f"      {'void r':>7} {'|c_u-1|':>10}")
    ds = []
    for r in (0.0, 0.10, 0.15, 0.20, 0.22):
        d = abs(c_u(void3(n, r, rng)) - 1.0)
        ds.append(d)
        print(f"      {r:>7.2f} {d:>10.5f}")
    p3 = ds[0] < 0.005 and ds[-1] > 0.012 and ds[-1] > 5 * ds[0]
    print(f"      -> defect ~0 with no void, rising with void volume to {ds[-1]:.3f} ({ds[-1]/ds[0]:.0f}x baseline): {p3}")
    print()

    print("=" * 100)
    print("  ASSESSMENT")
    print("=" * 100)
    print(f"  (1) null pairs collapse STRONG separation to ~0                  : {p1}")
    print(f"  (2) WEAK c_u unchanged by null pairs (survives the null cone)    : {p2}")
    print(f"  (3) only the codim-0 void biases c_u (density is the binding channel): {p3}")
    ok = p1 and p2 and p3
    print()
    msg = ("CONFIRMS the A2-core characterization -- the weak (H^1-tested) A6 discrepancy is blind to the "
           "codim>=1 null-cone vanishing-separation set; the residual is purely the codim-0 bulk node-density "
           "(void/excess) equidistribution, which matched scaling does not force") if ok else "MISMATCH -- revise"
    print("  RESULT: " + msg)


if __name__ == "__main__":
    main()
