#!/usr/bin/env python3
"""probe_selection_familycount.py — selection-rule family-count bound.

A297-A300 pursued the localizing factor of OI-287-1: why family n draws
the degree-n coefficient of the density rho. This probe asks a different,
boundary question: how many families does that rule permit?

The density is fixed by rho'(x) = 48 pi^3 x^2 + 6 pi^2 x + 2 pi, so
rho(x) = 16 pi^3 x^3 + 3 pi^2 x^2 + 2 pi x -- a CUBIC, with exactly three
nonzero coefficients (2 pi, 3 pi^2, 16 pi^3), the (2,3,16) target. The
index-matching rule (family n <-> degree-n coefficient) therefore has a
fourth term only if rho has a degree-4 coefficient. It does not.

The fourth eigenstate exists as an operator mode (node count 3), but the
coefficient it would draw, c_4, is identically zero. So the rule bounds
the family count at deg(rho) = 3 with no extra postulate. This matches
P20's topological 'no fourth family' (Z_3 = pi_1(S^3/Z_3) has three
eigenspaces) from the lepton-coupling side rather than the topology side.

HIT (supports no-fourth-family): c_4 = 0 exactly AND families 1-3 keep
their nonzero (1, 1.12, 7.50) lead AND family 4 coupling = 0.
NULL/tension: c_4 != 0, or family 4 coupling is a sensible nonzero value.

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

import numpy as np

PI = math.pi
OM = 4 * PI**3 + PI**2 + PI
ES = 13.177
N = 1500
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "selection_rule_familycount.json")

# rho(x) = 16 pi^3 x^3 + 3 pi^2 x^2 + 2 pi x : coefficients c_k in units pi^k
C = {1: 2.0, 2: 3.0, 3: 16.0, 4: 0.0}  # deg-4 term absent in a cubic


def eigfuns(potential, m=4):
    dx = 1.0 / N
    x = np.linspace(dx, 1 - dx, N - 1)
    off = -1.0 / dx**2 * np.ones(N - 2)
    H = np.diag(2.0 / dx**2 + potential(x)) + np.diag(off, 1) \
        + np.diag(off, -1)
    w, v = np.linalg.eigh(H)
    return x, w[:m], v[:, :m] / math.sqrt(dx)


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


def main():
    rho_p = lambda x: 48 * PI**3 * x**2 + 6 * PI**2 * x + 2 * PI
    Vpot = lambda x: rho_p(x)**2 / (2 * OM**2) + ES * x**2 * (1 - x)**2
    x, lam, psi = eigfuns(Vpot, 4)
    dx = x[1] - x[0]
    nu = [node_count(psi[:, n]) for n in range(4)]

    # family n (n=1..4) draws degree-n coefficient, A299 form (/ lambda_n)
    def diag_entry(n):
        k = n
        ck = C.get(k, 0.0)
        return float(np.sum(psi[:, n - 1]**2 * ck * PI**k * x**k) * dx) \
            / lam[n - 1]

    d = [diag_entry(n) for n in (1, 2, 3, 4)]
    r = [v / d[0] for v in d]

    c4 = C[4]
    fam4_zero = abs(r[3]) < 1e-12
    fam123_lead = (abs(r[1] - 1.12) < 0.02 and abs(r[2] - 7.5) < 0.05)
    hit = (c4 == 0.0) and fam4_zero and fam123_lead

    res = {"eigenvalues": [float(l) for l in lam], "node_counts": nu,
           "density": "rho = 16 pi^3 x^3 + 3 pi^2 x^2 + 2 pi x (cubic)",
           "coeffs_c_k": {str(k): C[k] for k in (1, 2, 3, 4)},
           "diag_raw": d, "diag_normalized": r,
           "family_count_bound": 3, "deg_rho": 3,
           "VERDICT": ("HIT: c_4 = 0 exactly; family 4 coupling = 0; "
                       "families 1-3 keep the (1, %.2f, %.2f) lead. The "
                       "selection rule bounds the family count at "
                       "deg(rho) = 3, matching P20's Z_3 no-fourth-family "
                       "from the coupling side." % (r[1], r[2]))
           if hit else
           ("NULL/tension: c_4 = %.3g, family-4 coupling = %.3g; the rule "
            "does not cleanly bound the family count at 3." % (c4, r[3]))}

    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    json.dump(res, open(OUT, "w"), indent=1)
    print("eigenvalues:", [round(float(l), 2) for l in lam])
    print("node counts:", nu, " (0,1,2,3 -> 4th mode exists)")
    print("c_k:", {k: C[k] for k in (1, 2, 3, 4)}, " -> c_4 =", c4)
    print("diagonal normalized:", [round(v, 4) for v in r])
    print("VERDICT:", res["VERDICT"])
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
