#!/usr/bin/env python3
"""verify_P312.py -- Verifier for Addendum 312 (P18-T2 sub-result (i),
renormalization sector + remaining-sector survey).

Recomputes from scratch and asserts that the renormalization generator
R = r d/dr + Delta_psi of P18 S3.6 is FORCED, INCLUDING the additive
Delta_psi that A304 flagged as a free c-number. The dilation generator
D = r d/dr is the unique scale-transformation generator up to the affine
freedom a*D + b*I (eigenfunctions r^p, eigenvalue p). The additive shift b
is fixed by self-adjointness in the geometric measure: under integration by
parts in measure r^k dr the adjoint is D* = -D - (k+1) I, so the skew-
symmetric shift is c = (k+1)/2. For the boundary-field radial measure
r^{d-2} dr (d=4, weight r^2) this is c = 3/2 = (d-1)/2 = Delta_psi.

This is a verified per-category uniqueness result for a SECOND sector of
P18-T2 sub-result (i), on top of A304's layer cycle T_cycle. It also carries
a machine check of the fiber-Laplacian survey item (unique 2nd-order
rotation-invariant operator on S^1, up to scale). The other survey items
(Dirac^2, boundary Laplacian, V_self/rho/M formula-uniqueness) are argued in
the addendum, not asserted here. Sub-result (i) is ADVANCED, not closed;
A294's re-typing of the whole P18-T2 stands; (ii)=A310 and (iii)=A302 are
unchanged.

  S1  Dilation generator: eigenfunctions and affine uniqueness  - checks 1-3
  S2  Skew constant in the geometric measure forces Delta_psi    - checks 4-6
  S3  Survey: fiber Laplacian on S^1 (machine-checkable item)     - checks 7-8

Copyright Léon Fernando Vlegels -- CC BY 4.0
"""
import sys
import math

import numpy as np

PI = math.pi
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}")


D_DIM = 4
BOUNDARY_DIM = D_DIM - 1            # 3
FIELD_K = D_DIM - 2                 # 2  -> boundary-field radial measure r^2
BULK_K = D_DIM - 1                  # 3  -> full B^4 radial measure r^3
DELTA_PSI = BOUNDARY_DIM / 2.0      # 3/2


def Dapply(f, r):
    """D = r d/dr via central differences."""
    return r * np.gradient(f, r)


def bump(r, c0, wd):
    x = (r - c0) / wd
    out = np.zeros_like(r)
    m = np.abs(x) < 1
    out[m] = np.exp(-1.0 / (1.0 - x[m] ** 2))
    return out


def inner(f, g, r, k):
    h = r[1] - r[0]
    return float(np.sum(f * g * r ** k) * h)


def skew_constant(r, k, fs, gs):
    """Least-squares c making (D + c) skew-symmetric in measure r^k dr."""
    A, B = [], []
    for f in fs:
        Df = Dapply(f, r)
        for g in gs:
            Dg = Dapply(g, r)
            A.append(inner(f, Dg, r, k) + inner(Df, g, r, k))
            B.append(inner(f, g, r, k))
    A = np.array(A)
    B = np.array(B)
    return -float(np.sum(A * B) / (2.0 * np.sum(B * B)))


r = np.linspace(0.2, 6.0, 4000)

# S1: dilation eigenfunctions and affine uniqueness -------------------------
print("S1  Dilation generator D = r d/dr: eigenfunctions and affine uniqueness")

# (1) D r^p = p r^p for p = 1, 2 (eigenfunctions powers, eigenvalue p)
errs = []
for p in (1.0, 2.0):
    f = r ** p
    rel = np.abs(Dapply(f, r)[3:-3] - p * f[3:-3]) / (np.abs(p * f[3:-3]) + 1e-12)
    errs.append(float(np.max(rel)))
check(1, "D r^p = p r^p for p=1,2 (max rel resid %.2e)" % max(errs),
      max(errs) < 5e-3)

# (2) D r^p = p r^p for a non-integer p (p = 1/2)
f = r ** 0.5
rel_half = float(np.max(np.abs(Dapply(f, r)[3:-3] - 0.5 * f[3:-3]) /
                        (np.abs(0.5 * f[3:-3]) + 1e-12)))
check(2, "D r^p = p r^p for p=1/2 (eigenvalue = scaling dimension; resid %.2e)"
      % rel_half, rel_half < 5e-3)

# (3) affine uniqueness: a*D + b*I still has r^p eigenfunctions, eig a*p+b
a_s, b_s, p = 2.0, -0.5, 2.0
f = r ** p
lhs = a_s * Dapply(f, r) + b_s * f
rhs = (a_s * p + b_s) * f
aff_rel = float(np.max(np.abs(lhs[3:-3] - rhs[3:-3]) /
                       (np.abs(rhs[3:-3]) + 1e-12)))
check(3, "any scale-covariant 1st-order op is a*D + b*I (eig a*p+b; resid %.2e)"
      % aff_rel, aff_rel < 5e-3)

# S2: skew constant forces Delta_psi ----------------------------------------
print("S2  Skew-symmetry in the geometric measure forces Delta_psi = 3/2")
fs = [bump(r, 2.0, 1.0), bump(r, 3.0, 1.2), bump(r, 4.0, 0.8)]
gs = [bump(r, 2.5, 1.1), bump(r, 3.5, 0.9), bump(r, 2.2, 1.3)]

c_field = skew_constant(r, FIELD_K, fs, gs)   # k=2 -> 3/2
c_bulk = skew_constant(r, BULK_K, fs, gs)     # k=3 -> 2

# (4) boundary-field measure r^{d-2}=r^2 gives skew constant c = 3/2
check(4, "skew constant of (D+c) in measure r^2 dr is c = 3/2 (got %.6f, tol 1e-3)"
      % c_field, abs(c_field - DELTA_PSI) < 1e-3)

# (5) general law c=(k+1)/2: bulk r^3 gives c=2, distinct from the r^2 value
check(5, "general law c=(k+1)/2: bulk r^3 measure gives c=2 (got %.6f, tol 1e-3)"
      % c_bulk, abs(c_bulk - (BULK_K + 1) / 2.0) < 1e-3 and abs(c_bulk - 2.0) < 1e-3)

# (6) Delta_psi = 3/2 = (d-1)/2, d=4 -- the forced scaling dimension
check(6, "Delta_psi = 3/2 = (d-1)/2 with d=4 (forced by self-adjointness)",
      abs(DELTA_PSI - 1.5) < 1e-12
      and abs(DELTA_PSI - BOUNDARY_DIM / 2.0) < 1e-12
      and abs(c_field - DELTA_PSI) < 1e-3)

# S3: survey -- fiber Laplacian on S^1 (machine-checkable item) --------------
print("S3  Survey: fiber Laplacian on S^1 is the unique invariant 2nd-order op")
m = 256
phi = np.linspace(0.0, 2 * PI, m, endpoint=False)
dphi = phi[1] - phi[0]
Lap1 = np.zeros((m, m))
for i in range(m):
    Lap1[i, (i - 1) % m] = 1.0 / dphi ** 2
    Lap1[i, i] = -2.0 / dphi ** 2
    Lap1[i, (i + 1) % m] = 1.0 / dphi ** 2

s1_ok = True
for k in (1, 2, 3):
    fk = np.cos(k * phi)
    lam_disc = -(2.0 / dphi ** 2) * (1.0 - math.cos(k * dphi))
    rel = float(np.max(np.abs(Lap1 @ fk - lam_disc * fk)) /
                (np.max(np.abs(lam_disc * fk)) + 1e-12))
    if rel > 1e-8:
        s1_ok = False
# eigenvalue approaches the Laplace-Beltrami spectrum -k^2
lam1 = -(2.0 / dphi ** 2) * (1.0 - math.cos(dphi))
check(7, "Laplace-Beltrami on S^1: e^{ik phi} are eigenfunctions, eigenvalue "
      "-> -k^2 (lam_1=%.5f ~ -1)" % lam1,
      s1_ok and abs(lam1 + 1.0) < 1e-3)

Shift = np.zeros((m, m))
for i in range(m):
    Shift[i, (i + 1) % m] = 1.0
commutator = float(np.linalg.norm(Lap1 @ Shift - Shift @ Lap1))
check(8, "the operator is rotation-invariant: [Lap, rotation] = 0 (norm %.2e), "
      "so it is the unique 2nd-order invariant op up to scale" % commutator,
      commutator < 1e-9)

print(f"\n{'='*60}\nRESULT: {PASS} PASS / {FAIL} FAIL")
sys.exit(0 if FAIL == 0 else 1)