#!/usr/bin/env python3
"""verify_P330.py -- Verifier for Addendum 330 (P18-T2 sub-result (iii), the
GLOBAL and HIGHER-ORDER coefficient-uniqueness question; the Heart 2
grounding/orientation probe for the third sub-result).

Copyright Léon Fernando Vlegels -- CC BY 4.0

A302 established the FIRST-ORDER LOCAL form of sub-result (iii): treating
(alpha, gamma, zeta, beta) as FREE assembly weights on their four sectors of
P18's radial operator (P18 l.363-372), the spectrum-response Jacobian at the
canonical point has full rank 4 (sigma_min/sigma_max = 9.3e-3 on l=0,
strengthening to 3.5e-2 with l<=2), so the coefficients are LOCALLY the unique
spectrum-matching assignment; the density (alpha) and renormalization (zeta)
sectors are ~99.86% spectrally collinear -- the SOFT DIRECTION where uniqueness
is fragile. A302 is explicitly silent on global injectivity and second order.

This probe RE-GROUNDS A302 (S1: exact reproduction of the rank-4 + collinearity
finding) and takes the first concrete steps on the global/higher-order frontier:

  S1  Re-ground A302              - checks 1-3   (rank 4; soft dir; strengthens)
  S2  Free-vs-tied resolution     - checks 4-7   (the corpus tie zeta=alpha^5/4
                                                  removes the soft PLANE)
  S3  Second-order (Hessian)      - checks 8-10  (residual-Hessian is positive-
                                                  definite => strict local min;
                                                  strengthens with sectors)
  S4  Global second-preimage      - checks 11-13 (FREE: a near-flat zeta valley
                                                  is the GLOBAL fragility; TIED:
                                                  a sharp isolated well)

THE CENTRAL FINDING (free-vs-tied). A302 treated all four coefficients as FREE
weights -- the hardest version, and the one P18 line 376 supports ("the four
coefficients ... are pinned by the references", each INDEPENDENTLY). P18 also
fixes zeta = alpha^{5/4} (l.348, l.358, l.579, l.608), a FUNCTIONAL tie of zeta
to alpha. The soft direction A302 named is, in raw coordinates, almost pure
'move zeta alone' (the raw near-null right-singular vector of the alpha-zeta
block is (0.014, 0.99990)): zeta is spectrally cheap (||dlambda/dzeta|| is ~71x
smaller than ||dlambda/dalpha||). The corpus tie FORBIDS moving zeta alone --
under zeta=alpha^{5/4}, dzeta/dalpha = (5/4) alpha^{1/4}, so any motion in zeta
drags alpha along, and alpha has the large response. The tie direction (1,
dz/da) is TRANSVERSE to the raw flat direction (|overlap| = 0.356), and the
tied single-direction spectral response is ~820x the free flat-direction
response. So: the alpha-zeta soft PLANE that carries A302's fragility EXISTS
only in the free-weight version; under the corpus tie the plane collapses to a
1-D curve and the fragility is dissolved. The honest verdict: sub-result (iii)
for the corpus-TIED family is far less fragile than A302's free-weight bound
suggests; the free-weight global question is the genuinely hard one, and S4
exhibits its fragility directly (a near-flat zeta valley) and the tie's cure (a
sharp well along the tied curve).
"""
import sys
import math

import numpy as np

PASS = FAIL = 0


def check(n, desc, cond):
    global PASS, FAIL
    ok = bool(cond)
    PASS += ok
    FAIL += (not ok)
    print(f"  [{'PASS' if ok else 'FAIL'}] {n:>2}. {desc}")


PI = math.pi
N = 600
NM = 6
A0 = 1.0 / 137.036
Z0 = A0 ** 1.25          # the corpus tie value zeta = alpha^{5/4} (P18 l.348)
ES = 13.177
M0 = 137.036
R0 = A0 ** 1.25
DZ_DA = 1.25 * A0 ** 0.25  # d zeta / d alpha along the tie zeta=alpha^{5/4}


def rho(r):
    # P18/P03 density: 16 pi^3 r^3 + 3 pi^2 r^2 + 2 pi r (A302's exact form)
    return 16 * PI**3 * r**3 + 3 * PI**2 * r**2 + 2 * PI * r


def eigs(alpha, zeta, l, m=NM, nmesh=N):
    """Lowest m eigenvalues of P18's radial operator with free density coupling
    alpha and free renormalization generator zeta, angular momentum l. The
    kinetic, self-lensing sectors are held canonical (A302's construction,
    reproduced exactly at nmesh=N)."""
    dr = 1.0 / nmesh
    r = np.linspace(dr, 1 - dr, nmesh - 1)
    n = len(r)
    d2 = (np.diag(-2 * np.ones(n)) + np.diag(np.ones(n - 1), 1)
          + np.diag(np.ones(n - 1), -1)) / dr**2
    d1 = (np.diag(np.ones(n - 1), 1) - np.diag(np.ones(n - 1), -1)) / (2 * dr)
    Vs = (ES / M0**2) * (1 - np.exp(-r / R0))
    H = (-d2 - np.diag(3 / r) @ d1 + np.diag(l * (l + 2) / r**2)
         + np.diag(Vs) + alpha * np.diag(rho(r)) + zeta * np.diag(r) @ d1)
    return np.sort(np.linalg.eigvals(H).real)[:m]


# moment ratios and family character (A302's fixed sectors)
dr = 1.0 / N
x = np.linspace(dr, 1 - dr, N - 1)
mu = [float(np.sum(x**k * rho(x)) * dr) for k in range(NM)]
mur = [m / mu[0] for m in mu]
eps = 1e-6
reom = [math.cos(2 * PI * k / 3) for k in range(3)]


def col_alpha(l):
    return (eigs(A0 + eps, Z0, l) - eigs(A0 - eps, Z0, l)) / (2 * eps)


def col_zeta_central(l):
    return (eigs(A0, Z0 + eps, l) - eigs(A0, Z0 - eps, l)) / (2 * eps)


def jac_free(sectors):
    """A302's free-weight response Jacobian, columns (alpha,gamma,zeta,beta)."""
    rows = []
    for l in sectors:
        dA = col_alpha(l)
        dZ = col_zeta_central(l)
        for k in range(3):
            for ni in range(NM):
                rows.append([dA[ni], reom[k], dZ[ni], mur[ni]])
    return np.array(rows)


# ----------------------------------------------------------------------
print("S1  Re-ground A302 (first-order local rank-4 + the soft direction)")

J0 = jac_free([0])
Jn0 = J0 / np.linalg.norm(J0, axis=0)
sv0 = np.linalg.svd(Jn0, compute_uv=False)
caz0 = abs(float(Jn0[:, 0] @ Jn0[:, 2]))
r0 = sv0[-1] / sv0[0]
J1 = jac_free([0, 1]); Jn1 = J1 / np.linalg.norm(J1, axis=0)
sv1 = np.linalg.svd(Jn1, compute_uv=False); r1 = sv1[-1] / sv1[0]
J2 = jac_free([0, 1, 2]); Jn2 = J2 / np.linalg.norm(J2, axis=0)
sv2 = np.linalg.svd(Jn2, compute_uv=False); r2 = sv2[-1] / sv2[0]

check(1, "A302 reproduced: free-weight Jacobian rank exactly 4, sigma_min/"
      "sigma_max=%.2e>1e-3 (l=0)" % r0,
      np.sum(sv0 > 1e-6) == 4 and r0 > 1e-3)
check(2, "A302 reproduced: alpha-zeta spectrally collinear ~0.99+ (soft dir): "
      "%.4f" % caz0, caz0 > 0.99)
check(3, "A302 reproduced: local uniqueness strengthens with sectors %.2e<%.2e"
      "<%.2e" % (r0, r1, r2), r0 < r1 < r2)

# ----------------------------------------------------------------------
print("S2  Free-vs-tied: the corpus tie zeta=alpha^{5/4} removes the soft PLANE")

# raw (un-normalized) alpha-zeta block over l=0,1,2; the soft direction in
# RAW coordinates is the near-null right-singular vector.
A = np.concatenate([col_alpha(l) for l in [0, 1, 2]])
Zc = np.concatenate([col_zeta_central(l) for l in [0, 1, 2]])
Jraw = np.column_stack([A, Zc])
_, Sr, Vtr = np.linalg.svd(Jraw)
soft_raw = Vtr[-1]                      # raw near-flat direction in (alpha,zeta)
soft_raw = soft_raw / np.linalg.norm(soft_raw)
na, nz = np.linalg.norm(A), np.linalg.norm(Zc)

check(4, "raw soft direction is almost pure 'move zeta alone' (|zeta-comp|"
      ">0.95): %s" % np.round(soft_raw, 4),
      abs(soft_raw[1]) > 0.95 and abs(soft_raw[0]) < 0.1)
check(5, "zeta is spectrally cheap: ||dlam/dalpha||/||dlam/dzeta|| > 10 (got "
      "%.1f)" % (na / nz), na / nz > 10)

# the tie direction (1, dz/da); is it transverse to the flat direction?
tievec = np.array([1.0, DZ_DA]); tievec = tievec / np.linalg.norm(tievec)
overlap = abs(float(tievec @ soft_raw))
check(6, "the tie direction (1, dzeta/dalpha) is TRANSVERSE to the soft "
      "direction (|overlap|<0.5): %.3f" % overlap, overlap < 0.5)

# tied single-direction response vs free flat-direction response
tied_resp = np.concatenate([col_alpha(l) + DZ_DA * col_zeta_central(l)
                            for l in [0, 1, 2]])
flat_resp = Jraw @ soft_raw
ratio = np.linalg.norm(tied_resp) / np.linalg.norm(flat_resp)
check(7, "tied-curve spectral response >> free flat-direction response (ratio"
      ">100): %.0f -- the soft PLANE collapses to a stiff curve" % ratio,
      ratio > 100)

# ----------------------------------------------------------------------
print("S3  Second-order: the residual-Hessian is positive-definite (strict min)")

# Residual R(c)=sum_modes (lambda_i(c)-lambda_i(c*))^2. At the root c*, R=0 and
# grad R=0, and the Hessian is exactly 2 J^T J (Gauss-Newton exact at a zero
# residual). Positive-definiteness => strict local minimum (second-order
# rigidity), upgrading A302's non-degenerate critical point.
H_free = 2 * J2.T @ J2
w_free = np.linalg.eigh(H_free)[0]
soft_eig = w_free[0]

check(8, "FREE residual-Hessian is positive-definite (all eig>0): %s -> strict"
      " local min, second-order rigidity" % np.round(w_free, 2),
      bool(np.all(w_free > 0)))

# second-order rigidity strengthens with sector coverage, mirroring first order
soft_curv = []
for secs in ([0], [0, 1], [0, 1, 2]):
    Js = jac_free(secs)
    soft_curv.append(float(np.linalg.eigh(2 * Js.T @ Js)[0][0]))
check(9, "second-order rigidity strengthens with sectors %.3f<%.3f<%.3f"
      % tuple(soft_curv), soft_curv[0] < soft_curv[1] < soft_curv[2])

# TIED residual-Hessian (coords alpha,gamma,beta; zeta slaved). Still PD, and
# the alpha-zeta soft plane is GONE (only 3 coords; no flat zeta direction).
rows_t = []
for l in [0, 1, 2]:
    dAt = col_alpha(l) + DZ_DA * col_zeta_central(l)
    for k in range(3):
        for ni in range(NM):
            rows_t.append([dAt[ni], reom[k], mur[ni]])
Jt = np.array(rows_t)
w_tied = np.linalg.eigh(2 * Jt.T @ Jt)[0]
check(10, "TIED residual-Hessian (alpha,gamma,beta) positive-definite, soft "
      "plane removed: softest curv %.3f (all eig>0: %s)"
      % (w_tied[0], bool(np.all(w_tied > 0))), bool(np.all(w_tied > 0)))

# ----------------------------------------------------------------------
print("S4  Global second-preimage: FREE zeta-valley vs TIED sharp well")

# Residual to the canonical spectrum with gamma,beta re-optimized (they enter as
# mode-structured additive shifts). Coarser mesh for the scan (the global claim
# is structural, not a sub-e-12 numeric).
NS = 200
g0, b0 = 0.75, 3 * PI / 20
drs = 1.0 / NS
xs = np.linspace(drs, 1 - drs, NS - 1)
mus = [float(np.sum(xs**k * rho(xs)) * drs) for k in range(NM)]
murs = [m / mus[0] for m in mus]
target = {l: eigs(A0, Z0, l, nmesh=NS) for l in [0, 1, 2]}


def resid(alpha, zeta):
    base, G, B, T = [], [], [], []
    for l in [0, 1, 2]:
        e = eigs(alpha, zeta, l, nmesh=NS)
        for k in range(3):
            for ni in range(NM):
                base.append(e[ni]); G.append(reom[k]); B.append(murs[ni])
                T.append(target[l][ni] + g0 * reom[k] + b0 * murs[ni])
    base = np.array(base); D = np.column_stack([np.array(G), np.array(B)])
    T = np.array(T)
    coef, *_ = np.linalg.lstsq(D, T - base, rcond=None)
    return float(np.sqrt(np.mean((base + D @ coef - T)**2)))

r_canon = resid(A0, Z0)
# FREE: scan zeta at fixed alpha -- the near-flat valley (the global fragility)
r_z2 = resid(A0, 2 * Z0)
r_z4 = resid(A0, 4 * Z0)
# TIED: scan along the tied curve zeta=alpha^{5/4} -- a sharp isolated well
r_a15 = resid(1.5 * A0, (1.5 * A0) ** 1.25)
r_a20 = resid(2.0 * A0, (2.0 * A0) ** 1.25)

check(11, "canonical point is the global min: residual ~machine zero (%.2e)"
      % r_canon, r_canon < 1e-6)
check(12, "FREE: a NEAR-FLAT zeta valley carries the global fragility -- zeta "
      "to 4x raises residual only to %.2e (<<the tied rise)" % r_z4,
      r_z4 < 5e-3 and r_z2 < r_z4)
check(13, "TIED: the curve zeta=alpha^{5/4} is a SHARP isolated well -- alpha "
      "to 1.5x/2x raises residual to %.2e/%.2e (>>10x the free zeta valley)"
      % (r_a15, r_a20),
      r_a15 > 10 * r_z4 and r_a20 > r_a15)

# ----------------------------------------------------------------------
print(f"\n{'='*64}\nRESULT: {PASS} PASS / {FAIL} FAIL")
print("VERDICT: GROUNDED + STEP. A302's first-order local rank-4 reproduced; "
      "the free-vs-tied framing RESOLVED -- the alpha-zeta soft plane that "
      "carries A302's fragility exists only in the FREE-weight version; the "
      "corpus tie zeta=alpha^{5/4} (P18 l.348) collapses it to a stiff curve "
      "(tied response ~820x the free flat-direction response), dissolving the "
      "fragility. Second order: the residual-Hessian is positive-definite "
      "(strict local min), strengthening with sectors -- A302 upgraded from "
      "first to second order. Global: the FREE version has a near-flat zeta "
      "valley (the genuine global risk); the TIED curve is a sharp isolated "
      "well. Sub-result (iii) for the corpus-tied family is far less fragile "
      "than A302's free-weight bound; the FREE-weight global injectivity (the "
      "second-preimage scan / a structural monotonicity argument) is the "
      "warranted fan-out. Global uniqueness stays OPEN; no number changes.")
sys.exit(0 if FAIL == 0 else 1)
