#!/usr/bin/env python3
# Copyright Léon Fernando Vlegels -- CC BY 4.0
"""probe_selection_hopf.py -- Probe for OI-287-1 (A317). THE KEYSTONE PROBE.

FIT-VS-FORCED probe. A316 closed the family-to-mode route in the REAL radial
reduction and named the real blocker: the three charged-lepton families are the
three COMPLEX cube-root eigenspaces of the layer-cycle (P20: omega^k, k=0,1,2,
omega = exp(2 pi i / 3)), but the real radial reduction folds omega and omega^2
together. Concretely, the real layer-cycle shift Re(omega^k) gamma =
gamma*(1, -1/2, -1/2) takes EQUAL values on k=1 and k=2, so families 2 and 3
collapse to the same radial operator and the real reduction structurally cannot
produce three distinct magnitude entries. A316 verdict NULL closest L2=0.399,
mu=1.222 (the pure-radial 4D operator). A316's redirection: move to a COMPLEX
realization where omega^k act as genuine distinct phases.

THE CONSTRUCTION THIS PROBE BUILDS (the point of the probe).

P18 S2.2 puts the master operator on sections over B^4 with boundary S^3,
Hopf-fibered S^1 -> S^3 -> S^2, fiber coordinate psi in [0, 2pi). On the Hopf
fiber a field decomposes into fiber-momentum sectors e^{i m psi}, m integer
(P18 S3.3 fiber Laplacian Delta_{S^1} = d^2/d psi^2, eigenvalues of -Delta_{S^1}
are m^2, m in Z). The layer-cycle T_cycle acts as psi -> psi + 2pi/3 on the
Hopf fiber (P18 S3.5 / S2.5), so on the m-sector it multiplies by
e^{i m 2pi/3} = omega^m. Therefore the omega^k eigenspace of T_cycle is exactly
{ m : m == k (mod 3) }, and its lowest representative is m = k. CRUCIALLY the
fiber momentum m enters the reduced radial operator as a genuine m-DEPENDENT
term, and m = 0, 1, 2 give THREE DISTINCT values -- unlike Re(omega^k), which
gave only two. This is exactly the complex structure the real reduction discarded:
the fiber-Laplacian eigenvalue m^2 separates m=1 from m=2 (1 vs 4), where
Re(omega) = Re(omega^2) = -1/2 could not.

BASE RADIAL OPERATOR (held fixed, verbatim from A311/A313/A316 C_4Dself):
  O = -d^2/dr^2 - (3/r) d/dr + V_eff(r), self-adjoint in r^3 dr,
  V_eff(r) = rho'(r)^2/(2 Omega^2) + E_self r^2 (1-r)^2,
  rho'(r) = 48 pi^3 r^2 + 6 pi^2 r + 2 pi, Omega = 4 pi^3 + pi^2 + pi,
  E_self = 13.177, rho(r) = 16 pi^3 r^3 + 3 pi^2 r^2 + 2 pi r,
  C = (2, 3, 16) -> target normalized (1, 1.5, 8). N = 1500.
  Eigensolver gen_eigh (numpy Cholesky reduction; scipy unavailable), the
  4D-radial assembly radial4d_matrices, and the box control
  sqrt(2) sin((n+1) pi r) in r^3 dr are reused VERBATIM from A313/A316.

THE FIBER-MOMENTUM TERM (two principled corpus-faithful forms, BOTH declared):
  (F-lap)  The fiber-Laplacian eigenvalue as a FLAT ADDITIVE term: the operator
           for family k gains + m_k^2 (the -Delta_{S^1} eigenvalue, P18 S3.3),
           m_k = k. This is the fiber as an internal S^1 of fixed circumference:
           its Laplacian contributes a constant energy m^2 to the radial sector.
           O_m = -d^2/dr^2 - (3/r) d/dr + V_eff(r) + m^2.
  (F-Hopf) The Hopf-connection / fiber-on-the-cone form: in the cone metric
           ds^2 = dr^2 + r^2 dOmega_3^2 (P18 S2.1) the Hopf fiber S^1 sits inside
           the radius-r boundary S^3, so the fiber circle scales with r and its
           Laplacian eigenvalue m^2 enters with the 1/r^2 weight that every
           angular/fiber mode on a cone carries (the same 1/r^2 that the angular
           term l(l+2)/r^2 carries in P18 S2.2). O_m = -d^2/dr^2 - (3/r) d/dr +
           m^2/r^2 + V_eff(r). This is the geometrically faithful reading: the
           fiber momentum couples through the Hopf/cone structure as an
           m-dependent centrifugal term, m_k = k. (For m=0 both forms reduce to
           the pure-radial C_4Dself, the A316 node-count ground state -> the
           construction contains the A316 closest map as its k=0 family.)

FAMILY-TO-STATE MAP (CANONICAL, declared before computing):
  family k (k = 0, 1, 2; i.e. lepton generation n = k + 1) -> the GROUND radial
  state of the m = k fiber-momentum sector (the omega^k eigenspace). Three
  GENUINELY DISTINCT operators (m = 0, 1, 2 give distinct m^2 = 0, 1, 4 or
  m^2/r^2 terms), hence three distinct eigenstates and three distinct overlaps:
  NO degeneracy. This is the corpus family definition (P20: families ARE the
  omega^k eigenspaces), realized complexly so omega^k are distinct phases.

The diagonal entry n (n = k+1) is the degree-n density overlap of the m=k ground
state, divided by its eigenvalue (A299/A311 normalization), all in r^3 dr:
  R[n][k] = <psi_{m_n} | C[k-1] pi^k r^k | psi_{m_n}>_{r^3}
            / <psi_{m_n}|psi_{m_n}>_{r^3} / lambda_{m_n},  C = (2,3,16).
Report normalized diagonal, mu = ratio[1], tau = ratio[2], L2 to (1,1.5,8),
valid_peak (column-standardized argmax of row n == n for all n AND the box
control fails the same test in the same r^3 measure).

PRE-REGISTERED VERDICT (declared before any number is computed):
  HIT / FORCED (CLOSES the magnitude piece of OI-287-1) iff the CANONICAL complex
    map (family k -> ground state of the m = k fiber-momentum sector, the omega^k
    eigenspace) reproduces (1, 1.5, 8) to within 20% on all three entries AND
    restores a valid column-standardized diagonal peak, with the fiber term in
    its corpus-faithful form (F-lap or F-Hopf, both declared) and NO fitting.
    This would be the derivation the whole A296-A317 arc has sought.
  STRENGTHENED-PARTIAL iff the complex realization produces THREE DISTINCT entries
    (breaking the A316 degeneracy -- itself a real structural result) and moves mu
    toward 1.5 (mu in [1.3, 1.7]) but does not fully hit or peak. Report which
    fiber form (F-lap vs F-Hopf) and the diagonal.
  NULL iff the complex realization does not improve on the real-reduction closest
    (L2 = 0.399, mu = 1.222) -- then even the genuine omega^k phase does not
    supply the magnitude, the residual is elsewhere, and this is a major
    redirection. Report clearly.

CRITICAL HONESTY: the canonical map is family k -> m = k sector (the omega^k
eigenspace, the corpus family definition). A HIT requires THAT map, not a tuned
m-assignment. Both fiber forms are declared up front; if a non-canonical
m-assignment is needed to hit, that is a FIT and is flagged as such. tau and peak
are reported, not just mu. The first deliverable is simply: does the complex
realization BREAK the A316 degeneracy and give three distinct entries? Confirmed
explicitly below. No NaN; the generalized-eigenproblem weight r^3 > 0 on (0,1).
"""
import os
import json
import math

import numpy as np


def gen_eigh(A, Bm):
    """Generalized symmetric-definite eig A v = lambda W v, W SPD, via numpy
    Cholesky reduction (scipy unavailable). Verbatim from A313/A316.
        W = L L^T ; C = L^{-1} A L^{-T} ; C y = lambda y ; v = L^{-T} y.
    """
    L = np.linalg.cholesky(Bm)
    Linv = np.linalg.inv(L)
    Cmat = Linv @ A @ Linv.T
    Cmat = 0.5 * (Cmat + Cmat.T)
    w, y = np.linalg.eigh(Cmat)
    v = Linv.T @ y
    return w, v


PI = math.pi
OMEGA = 4 * PI**3 + PI**2 + PI            # alpha^{-1}
E_SELF = 13.177
C = (2.0, 3.0, 16.0)                     # raw target coefficients
TARGET = [c / C[0] for c in C]           # normalized: (1, 1.5, 8)
N = 1500

# CANONICAL family-to-state map: family k (k=0,1,2) -> m = k fiber-momentum
# sector (the omega^k eigenspace, P20). Lowest representative of {m == k mod 3}.
CANONICAL_M = [0, 1, 2]
# fiber forms, both declared up front
FIBER_FORMS = ["F_lap", "F_Hopf"]

HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "selection_rule_hopf.json")


def rho_p(x):
    return 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI


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


def Vpot(x):
    return rho_p(x)**2 / (2 * OMEGA**2) + E_SELF * x**2 * (1 - x)**2


def grid():
    dx = 1.0 / N
    x = np.linspace(dx, 1 - dx, N - 1)
    return x, dx


def radial4d_matrices(x, dx, V_on_grid, weight):
    """Symmetric 4D-radial operator A and mass W, self-adjoint in weight dr.
    Verbatim from A313/A316 (radial4d_matrices)."""
    n = N - 1
    w_mid = 0.5 * (weight[:-1] + weight[1:])
    w_mid_left0 = max(0.5 * (0.0 + weight[0]), 0.0)
    A = np.zeros((n, n))
    for i in range(n):
        wl = w_mid[i - 1] if i > 0 else w_mid_left0
        wr = w_mid[i] if i + 1 < n else 0.5 * (weight[i] + weight[i])
        A[i, i] = (wl + wr) / dx**2
        if i + 1 < n:
            A[i, i + 1] = -w_mid[i] / dx**2
            A[i + 1, i] = -w_mid[i] / dx**2
    A += np.diag(weight * V_on_grid)
    A = 0.5 * (A + A.T)
    return A, np.diag(weight)


def solve_fiber_sector(m, form, n_modes=2):
    """Solve the m-fiber-momentum-sector operator in r^3 dr.

    F-lap : O_m = -d^2/dr^2 - (3/r)d/dr + V_eff(r) + m^2          (flat additive)
    F-Hopf: O_m = -d^2/dr^2 - (3/r)d/dr + m^2/r^2 + V_eff(r)      (cone/Hopf 1/r^2)

    Returns x, dx, eigenvalues, r^3-normalized eigenvectors (lowest n_modes).
    The GROUND mode (index 0) is the family-k state for m = k.
    """
    x, dx = grid()
    w_arr = x**3
    if form == "F_lap":
        fiber_term = float(m**2) * np.ones_like(x)
    elif form == "F_Hopf":
        fiber_term = (float(m**2) / x**2) if m > 0 else np.zeros_like(x)
    else:
        raise ValueError("unknown fiber form: %s" % form)
    A, Wm = radial4d_matrices(x, dx, Vpot(x) + fiber_term, w_arr)
    w, v = gen_eigh(A, Wm)
    order = np.argsort(w)
    w = w[order]
    v = v[:, order]
    vv = v[:, :n_modes].copy()
    for k in range(n_modes):
        nrm = math.sqrt(float(np.sum(vv[:, k]**2 * w_arr) * dx))
        vv[:, k] = vv[:, k] / nrm
    return x, dx, w[:n_modes], vv


def overlap_row(col, x, dx, w_arr, lam):
    """Degree-k normalized density overlap row for one mode, divided by lam
    (A299/A311 normalization), measure = r^3 (w_arr)."""
    norm = float(np.sum(col**2 * w_arr) * dx)
    return [float(np.sum(col**2 * C[k - 1] * PI**k * x**k * w_arr) * dx)
            / norm / lam for k in (1, 2, 3)]


def node_count(col):
    s = np.sign(col[np.abs(col) > 1e-9])
    return int(np.sum(s[1:] != s[:-1]))


def diag_peaks(M):
    """Column-standardized argmax of column k == k for all k."""
    Z = np.zeros_like(M)
    for k in range(3):
        c = M[:, k]
        Z[:, k] = (c - c.mean()) / c.std() if c.std() > 0 else 0.0
    return all(int(np.argmax(Z[:, k])) == k for k in range(3))


def norm_diag(M):
    d = [M[n, n] for n in range(3)]
    return [v / d[0] for v in d]


def l2(ratio):
    return math.sqrt(sum((ratio[i] - TARGET[i])**2 for i in range(3)))


def box_states(x):
    psi_box = np.column_stack([math.sqrt(2.0) * np.sin((n + 1) * PI * x)
                               for n in range(3)])
    lam_box = np.array([(n + 1)**2 * PI**2 for n in range(3)])
    return psi_box, lam_box


def box_matrix(psi_box, lam_box, x, dx, w_arr):
    M = np.zeros((3, 3))
    for n in range(3):
        M[n, :] = overlap_row(psi_box[:, n], x, dx, w_arr, lam_box[n])
    return M


def build_canonical(form, x, dx, w_arr):
    """Build the canonical complex map for one fiber form: family k -> ground
    state of m=k sector. Returns the 3x3 B-matrix, eigenvalues, node counts."""
    states = np.zeros((N - 1, 3))
    lams = np.zeros(3)
    nodes = [0, 0, 0]
    for k in range(3):
        m = CANONICAL_M[k]
        _, _, w_m, v_m = solve_fiber_sector(m, form, n_modes=2)
        states[:, k] = v_m[:, 0]      # GROUND mode of the m=k sector
        lams[k] = w_m[0]
        nodes[k] = node_count(v_m[:, 0])
    M = np.zeros((3, 3))
    for n in range(3):
        M[n, :] = overlap_row(states[:, n], x, dx, w_arr, lams[n])
    return M, lams, nodes


def evaluate_form(form, x, dx, w_arr, psi_box, lam_box):
    M, lams, nodes = build_canonical(form, x, dx, w_arr)
    B_box = box_matrix(psi_box, lam_box, x, dx, w_arr)
    ratio = norm_diag(M)
    peak = bool(diag_peaks(M))
    box_peak = bool(diag_peaks(B_box))
    within20 = all(abs(ratio[i] - TARGET[i]) <= 0.20 * TARGET[i]
                   for i in range(3))
    # degeneracy guard (the A316 failure mode): two families share BOTH an
    # eigenvalue and an entire overlap row -> collapse to the same state.
    degenerate = False
    for a in range(3):
        for b in range(a + 1, 3):
            same_eig = abs(lams[a] - lams[b]) < 1e-6
            same_row = all(abs(M[a, q] - M[b, q]) < 1e-6 for q in range(3))
            if same_eig and same_row:
                degenerate = True
    # three distinct entries? (the first deliverable: degeneracy broken)
    diag_entries = [M[n, n] for n in range(3)]
    distinct = (abs(diag_entries[0] - diag_entries[1]) > 1e-6
                and abs(diag_entries[0] - diag_entries[2]) > 1e-6
                and abs(diag_entries[1] - diag_entries[2]) > 1e-6)
    return {
        "fiber_form": form,
        "canonical_m_assignment": CANONICAL_M,
        "diagonal_ratio": [round(float(v), 4) for v in ratio],
        "mu": round(float(ratio[1]), 4),
        "tau": round(float(ratio[2]), 4),
        "peaks_diagonal": peak,
        "box_control_peaks": box_peak,
        "valid_peak": bool(peak and not box_peak),
        "within_20pct": bool(within20),
        "hits_target": bool(within20 and peak and not box_peak),
        "L2_to_target": round(float(l2(ratio)), 4),
        "eigenvalues": [round(float(v), 4) for v in lams],
        "node_counts": nodes,
        "degenerate_selection": degenerate,
        "three_distinct_entries": bool(distinct),
        "mu_in_partial_band": bool(1.3 <= ratio[1] <= 1.7),
    }


def main():
    x, dx = grid()
    w_arr = x**3
    psi_box, lam_box = box_states(x)

    # the m=0 sector ground state must reproduce the A316 node-count ground
    # (C_4Dself ground), the closest real-reduction map: mu_diag entry sanity.
    # (m=0 is identical in both forms -> the pure-radial 4D operator.)
    forms = {}
    for form in FIBER_FORMS:
        forms[form] = evaluate_form(form, x, dx, w_arr, psi_box, lam_box)

    # A316 real-reduction closest reference (the bar to beat)
    A316_closest_L2 = 0.399
    A316_closest_mu = 1.222

    # the canonical map's outcome per form
    canonical_hits = [f for f in FIBER_FORMS if forms[f]["hits_target"]]
    closest_form = min(FIBER_FORMS, key=lambda f: forms[f]["L2_to_target"])

    # degeneracy broken? (BOTH forms should give three distinct entries since
    # m=0,1,2 are genuinely distinct operators)
    degeneracy_broken = all(forms[f]["three_distinct_entries"]
                            and not forms[f]["degenerate_selection"]
                            for f in FIBER_FORMS)

    # improvement over the A316 real-reduction closest?
    improves = any(forms[f]["L2_to_target"] < A316_closest_L2 - 1e-9
                   for f in FIBER_FORMS)
    moves_mu = any(abs(forms[f]["mu"] - 1.5) < abs(A316_closest_mu - 1.5) - 1e-9
                   for f in FIBER_FORMS)
    in_band = [f for f in FIBER_FORMS if forms[f]["mu_in_partial_band"]]

    # NaN guard
    nan_found = False
    for f in FIBER_FORMS:
        for val in forms[f]["diagonal_ratio"] + forms[f]["eigenvalues"]:
            if val != val:
                nan_found = True

    # ---- VERDICT (canonical map only counts for FORCED) ----------------
    if canonical_hits:
        verdict = "FORCED"
    elif degeneracy_broken and (in_band or (moves_mu and improves)):
        verdict = "STRENGTHENED_PARTIAL"
    elif degeneracy_broken and (moves_mu or improves):
        verdict = "STRENGTHENED_PARTIAL"
    else:
        verdict = "NULL"

    summary = {
        "addendum": "A317",
        "addresses": "OI-287-1",
        "builds_on": ["A311", "A313", "A316"],
        "thesis": ("the three families are the omega^k eigenspaces of the "
                   "layer-cycle (psi -> psi + 2pi/3 on the Hopf fiber); on the "
                   "Hopf fibration the omega^k eigenspace is the fiber-momentum "
                   "sector m == k (mod 3), m enters the radial operator as a "
                   "genuine m-dependent fiber-Laplacian term (m^2 or m^2/r^2), "
                   "and m=0,1,2 give THREE DISTINCT operators -- the complex "
                   "structure the real reduction (Re(omega^k)) discarded. "
                   "Canonical map: family k -> ground state of the m=k sector."),
        "target_normalized": TARGET,
        "constants": {"OMEGA_alpha_inv": OMEGA, "E_self": E_SELF},
        "canonical_m_assignment": CANONICAL_M,
        "fiber_forms_declared": FIBER_FORMS,
        "forms": forms,
        "F_lap_diagonal": forms["F_lap"]["diagonal_ratio"],
        "F_lap_mu": forms["F_lap"]["mu"],
        "F_lap_tau": forms["F_lap"]["tau"],
        "F_lap_L2": forms["F_lap"]["L2_to_target"],
        "F_lap_valid_peak": forms["F_lap"]["valid_peak"],
        "F_lap_three_distinct": forms["F_lap"]["three_distinct_entries"],
        "F_Hopf_diagonal": forms["F_Hopf"]["diagonal_ratio"],
        "F_Hopf_mu": forms["F_Hopf"]["mu"],
        "F_Hopf_tau": forms["F_Hopf"]["tau"],
        "F_Hopf_L2": forms["F_Hopf"]["L2_to_target"],
        "F_Hopf_valid_peak": forms["F_Hopf"]["valid_peak"],
        "F_Hopf_three_distinct": forms["F_Hopf"]["three_distinct_entries"],
        "A316_degeneracy_broken": bool(degeneracy_broken),
        "A316_real_reduction_closest_L2": A316_closest_L2,
        "A316_real_reduction_closest_mu": A316_closest_mu,
        "improves_on_A316_closest_L2": bool(improves),
        "moves_mu_toward_1.5_vs_A316": bool(moves_mu),
        "forms_with_mu_in_partial_band": in_band,
        "canonical_map_hits": sorted(canonical_hits),
        "closest_form": closest_form,
        "closest_form_diagonal": forms[closest_form]["diagonal_ratio"],
        "closest_form_L2": forms[closest_form]["L2_to_target"],
        "nan_found": nan_found,
        "closes_OI_287_1_magnitude": bool(verdict == "FORCED"),
        "verdict": verdict,
    }

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    with open(OUT, "w") as f:
        json.dump(summary, f, indent=2)

    print("=" * 74)
    print("A317 PROBE -- OI-287-1: the COMPLEX Hopf-fibered realization. Does the")
    print("   genuine omega^k fiber-momentum phase (m=0,1,2) break the A316")
    print("   degeneracy and supply the magnitude (1,1.5,8) with a peak?")
    print("=" * 74)
    print("target (normalized):", TARGET)
    print("canonical map: family k -> ground state of m=k fiber sector "
          "(omega^k eigenspace)")
    print(f"A316 real-reduction closest: L2={A316_closest_L2}, mu={A316_closest_mu}")
    print("-" * 74)
    for form in FIBER_FORMS:
        v = forms[form]
        print(f"{form:8s} diag={v['diagonal_ratio']} mu={v['mu']} tau={v['tau']} "
              f"peak={v['valid_peak']} <=20%={v['within_20pct']} "
              f"HIT={v['hits_target']} L2={v['L2_to_target']}")
        print(f"         eig={v['eigenvalues']} nodes={v['node_counts']} "
              f"3distinct={v['three_distinct_entries']} "
              f"degen={v['degenerate_selection']}")
    print("-" * 74)
    print(f"A316 DEGENERACY BROKEN (three distinct entries, both forms): "
          f"{degeneracy_broken}")
    print(f"improves on A316 closest L2 (<{A316_closest_L2}): {improves}")
    print(f"moves mu toward 1.5 vs A316 (1.222): {moves_mu}")
    print(f"forms with mu in [1.3,1.7]: {in_band}")
    print(f"canonical map hits (FORCED-eligible): {sorted(canonical_hits)}")
    print(f"closest form: {closest_form} diag={forms[closest_form]['diagonal_ratio']} "
          f"L2={forms[closest_form]['L2_to_target']}")
    print(f"NaN found: {nan_found}")
    print(f"\nVERDICT: {verdict}")
    print(f"A316 degeneracy broken: {degeneracy_broken} "
          f"(three distinct entries in the complex realization)")
    if verdict == "FORCED":
        print("  -> the CANONICAL complex map (family k -> m=k omega^k eigenspace)")
        print("     reproduces (1,1.5,8) within 20% with a valid peak, fiber term")
        print("     corpus-faithful, NO fitting. The magnitude piece of OI-287-1")
        print("     CLOSES. This is the derivation the arc sought.")
    elif verdict == "STRENGTHENED_PARTIAL":
        print("  -> the complex realization BREAKS the A316 degeneracy: three")
        print("     genuinely distinct entries from m=0,1,2 (real result). It moves")
        print("     mu and/or improves L2 vs the real-reduction closest, but does")
        print("     NOT fully hit (1,1.5,8) with a peak under the canonical map.")
    else:
        print("  -> the complex realization does NOT improve on the real-reduction")
        print("     closest (L2=0.399, mu=1.222). Even the genuine omega^k phase")
        print("     does not supply the magnitude; the residual is elsewhere. A")
        print("     major redirection.")
    print(f"wrote {OUT}")


if __name__ == "__main__":
    main()
