#!/usr/bin/env python3
"""probe_p18t2_coeff_rigidity.py — P18-T2 sub-result (iii), local form.

P18-T2 (assembly uniqueness) needs three sub-results; A294 parked all
three. The third -- the four coefficients (alpha, gamma, zeta, beta) are
the only spectrum-matching assignment -- is the one amenable to a sharp
local test: spectral rigidity.

Treat (alpha, gamma, zeta, beta) as FREE assembly weights on their four
sectors of P18's radial operator (S3, S4.3), holding the fixed sectors
(kinetic, self-lensing) at canonical:

  O_{l,k} = -d^2/dr^2 - (3/r) d/dr + l(l+2)/r^2
            + V_self(r)                      [fixed sector]
            + alpha * rho(r)                 [density coupling]
            + zeta  * (r d/dr)               [renormalization generator]
            + gamma * Re(omega^k)            [layer-cycle, omega=e^{2pi i/3}]
            + beta  * (mu_n/mu_0)            [moment hierarchy, P_n diagonal]

rho(r) = 16 pi^3 r^3 + 3 pi^2 r^2 + 2 pi r;  mu_n = int_0^1 x^n rho dx.
Canonical: alpha=1/137.036, gamma=3/4, zeta=alpha^{5/4}, beta=3pi/20.

The spectrum map lambda_{n,k}(alpha,gamma,zeta,beta) = [n-th radial
eigenvalue with density-weight alpha, renorm-weight zeta] + gamma Re(w^k)
+ beta (mu_n/mu_0). gamma and beta enter analytically; alpha and zeta are
finite-differenced through the operator. The Jacobian J = d lambda_i /
d c_j (rows = modes (n,k), cols = the four coefficients) is formed at the
canonical point and its singular values inspected.

HIT (local uniqueness; partial support for sub-result iii): J has
numerical rank 4, sigma_min/sigma_max > 1e-3 -- no flat direction, so the
canonical coefficients are locally the unique spectrum-match.
OBSTRUCTION (against uniqueness): sigma_min/sigma_max < 1e-6 -- a
coefficient combination leaves the spectrum invariant to first order;
report the flat direction. Pre-registered watch: alpha (density) vs beta
(moment) may be degenerate, since mu_n/mu_0 and <rho>_n are both
density-weighted.

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

import numpy as np

PI = math.pi
N = 600
NMODES = 6
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "p18t2_coeff_rigidity.json")

A0 = 1.0 / 137.036
G0 = 3.0 / 4.0
Z0 = A0 ** 1.25
B0 = 3.0 * PI / 20.0
ESELF = 13.177
M0 = 137.036
R0 = A0 ** 1.25


def rho(r):
    return 16 * PI**3 * r**3 + 3 * PI**2 * r**2 + 2 * PI * r


def radial_eigs(alpha, zeta, l=0, m=NMODES):
    """Lowest m real eigenvalues (by real part) of the radial operator
    with density-weight alpha and renorm-weight zeta. Fixed sectors
    (kinetic, V_self) at canonical."""
    dr = 1.0 / N
    r = np.linspace(dr, 1 - dr, N - 1)
    n = len(r)
    d2 = (np.diag(-2.0 * 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)
    Vself = (ESELF / M0**2) * (1 - np.exp(-r / R0))
    H = (-d2
         - np.diag(3.0 / r) @ d1
         + np.diag(l * (l + 2) / r**2)
         + np.diag(Vself)
         + alpha * np.diag(rho(r))
         + zeta * np.diag(r) @ d1)
    w = np.linalg.eigvals(H)
    w = np.sort(w.real[np.argsort(w.real)])  # by real part
    return np.sort(w.real)[:m]


def main():
    # moments mu_n / mu_0
    dr = 1.0 / N
    x = np.linspace(dr, 1 - dr, N - 1)
    rh = rho(x)
    mu = [float(np.sum(x**nn * rh) * dr) for nn in range(NMODES)]
    mu_ratio = [m / mu[0] for m in mu]

    # base spectrum + finite differences for alpha, zeta
    base = radial_eigs(A0, Z0)
    eps = 1e-6
    dA = (radial_eigs(A0 + eps, Z0) - radial_eigs(A0 - eps, Z0)) / (2 * eps)
    dZ = (radial_eigs(A0, Z0 + eps) - radial_eigs(A0, Z0 - eps)) / (2 * eps)
    reomega = [math.cos(2 * PI * k / 3) for k in range(3)]  # 1, -1/2, -1/2

    rows = []
    for k in range(3):
        for nidx in range(NMODES):
            rows.append([dA[nidx], reomega[k], dZ[nidx], mu_ratio[nidx]])
    J = np.array(rows)  # (3*NMODES) x 4, cols = [alpha, gamma, zeta, beta]

    # column-normalize so singular values compare responses, not scales
    norms = np.linalg.norm(J, axis=0)
    Jn = J / norms
    sv = np.linalg.svd(Jn, compute_uv=False)
    ratio = float(sv[-1] / sv[0])

    # which pair is most collinear? cosine between normalized columns
    cols = Jn
    names = ["alpha", "gamma", "zeta", "beta"]
    cos = {}
    for a in range(4):
        for b in range(a + 1, 4):
            c = abs(float(cols[:, a] @ cols[:, b]))
            cos["%s.%s" % (names[a], names[b])] = c

    rank4 = ratio > 1e-3
    obstruction = ratio < 1e-6
    res = {"canonical": {"alpha": A0, "gamma": G0, "zeta": Z0, "beta": B0},
           "mu_ratio": mu_ratio,
           "dlambda_dalpha": list(map(float, dA)),
           "dlambda_dzeta": list(map(float, dZ)),
           "singular_values_normalized": list(map(float, sv)),
           "sigma_min_over_max": ratio,
           "column_cosines": cos,
           "VERDICT": ("HIT: Jacobian full rank 4 (sigma_min/sigma_max="
                       "%.2e > 1e-3). The four assembly coefficients are "
                       "LOCALLY the unique spectrum-matching assignment; "
                       "no flat direction. Partial support for P18-T2 "
                       "sub-result (iii): local coefficient uniqueness." % ratio)
           if rank4 else
           ("OBSTRUCTION: flat direction (sigma_min/sigma_max=%.2e). A "
            "coefficient combination leaves the spectrum invariant; "
            "coefficient uniqueness fails locally." % ratio
            if obstruction else
            "INCONCLUSIVE: sigma_min/sigma_max=%.2e in the gray zone "
            "[1e-6, 1e-3]; near-degenerate, report and refine." % ratio)}
    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    json.dump(res, open(OUT, "w"), indent=1)
    print("canonical (alpha,gamma,zeta,beta) = (%.4g, %.3g, %.4g, %.4g)"
          % (A0, G0, Z0, B0))
    print("mu_n/mu_0      :", [round(v, 4) for v in mu_ratio])
    print("dlambda/dalpha :", [round(v, 3) for v in dA])
    print("dlambda/dzeta  :", [round(v, 4) for v in dZ])
    print("normalized singular values:", [round(float(s), 4) for s in sv])
    print("sigma_min/sigma_max = %.3e" % ratio)
    print("most collinear column pair:",
          max(cos, key=cos.get), "cos=%.4f" % max(cos.values()))
    print("VERDICT:", res["VERDICT"])
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
