#!/usr/bin/env python3
"""probe_selection_gapweight.py — selection-rule gap-weighting probe.

A298 found the derivative/nodal coupling carries the right (2,3,16)
direction but is monotone (no diagonal peak). A296's surviving
strength is kappa*Omega/lambda_1, which names eigenvalue data. This
tests whether a FIXED eigenvalue-gap weighting of the right-direction
coupling localizes it to the diagonal and moves magnitudes toward
(2,3,16). Small declared set of weightings, no fitting.

Base matrices (right-direction, from A298):
  D  derivative overlap   <psi_n'|x^k|psi_n'>
  R  density overlap      <psi_n| c_k pi^k x^k |psi_n>
Row weightings (row index n -> level), declared up front:
  W1  / lambda_n          (literal level-resolved strength)
  W2  / (lambda_n-lambda_1)   (gap)
  W3  * (lambda_n/lambda_1)   (inverse sign control)
  W4  / sqrt(lambda_n)
Criteria (A297/A298): T1 diagonal peak (argmax_k row n == n all n) AND
box fails; T2 normalized diagonal within 20% of (1,1.5,8). HIT iff
some declared cell meets T1 AND T2 verbatim; else NULL = sixth class,
report closest cell by L2 distance to (1,1.5,8).

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
AL = 1.0 / OM
KAP = AL**1.25
ES = 13.177
C = (2.0, 3.0, 16.0)
N = 1500
HERE = os.path.dirname(os.path.abspath(__file__))
OUT = os.path.join(HERE, "..", "..", "..", "..", "Outputs",
                   "selection_rule_gapweight.json")
TGT = [c / C[0] for c in C]


def eigfuns(potential):
    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[:3], v[:, :3] / math.sqrt(dx)


def diag_pref(M):
    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[n])) == n for n in range(3))


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)
    dx = x[1] - x[0]
    xb = x
    psib = np.column_stack([math.sqrt(2) * np.sin((n + 1) * PI * xb)
                            for n in range(3)])

    def derivM(p):
        d = np.gradient(p, dx, axis=0)
        return np.array([[float(np.sum(d[:, n]**2 * x**k) * dx)
                          for k in (1, 2, 3)] for n in range(3)])

    def densM(p):
        return np.array([[float(np.sum(p[:, n]**2 * C[k - 1] * PI**k
                                       * x**k) * dx)
                          for k in (1, 2, 3)] for n in range(3)])

    bases = {"D_derivative": (derivM(psi), derivM(psib)),
             "R_density": (densM(psi), densM(psib))}
    lam1 = lam[0]
    weightings = {
        "W1_over_lambda_n": np.array([1 / lam[n] for n in range(3)]),
        "W2_over_gap": np.array([1.0] + [1 / (lam[n] - lam1)
                                         for n in range(1, 3)]),
        "W3_times_lambda_ratio": np.array([lam[n] / lam1
                                           for n in range(3)]),
        "W4_over_sqrt_lambda": np.array([1 / math.sqrt(lam[n])
                                         for n in range(3)]),
    }

    res = {"eigenvalues": [float(l) for l in lam],
           "kappa_Omega_over_lambda1": float(KAP * OM / lam1),
           "target": TGT, "cells": {}}
    best = (1e9, None, None)
    hit = []
    for bname, (M, Mb) in bases.items():
        for wname, w in weightings.items():
            Mw = M * w[:, None]
            Mbw = Mb * w[:, None]
            peak = bool(diag_pref(Mw) and not diag_pref(Mbw))
            d = [float(Mw[n, n]) for n in range(3)]
            r = [v / d[0] for v in d]
            within20 = all(abs(r[i] - TGT[i]) <= 0.2 * TGT[i]
                           for i in range(3))
            l2 = math.sqrt(sum((r[i] - TGT[i])**2 for i in range(3)))
            cell = "%s + %s" % (bname, wname)
            res["cells"][cell] = {"diag_ratio": r, "peak": peak,
                                  "within_20pct": bool(within20),
                                  "L2_to_target": l2}
            if peak and within20:
                hit.append(cell)
            if l2 < best[0]:
                best = (l2, cell, r)
            print("  %-34s ratio (1, %.2f, %.2f) peak=%s L2=%.2f"
                  % (cell, r[1], r[2], peak, l2))

    res["VERDICT"] = ("HIT: " + ", ".join(hit) if hit else
                      "NULL: no declared gap-weighting localizes to "
                      "(2,3,16) (sixth excluded class). Closest: %s "
                      "ratio (1, %.2f, %.2f), L2=%.2f vs target (1,1.5,8)"
                      % (best[1], best[2][1], best[2][2], best[0]))
    os.makedirs(os.path.dirname(OUT), exist_ok=True)
    json.dump(res, open(OUT, "w"), indent=1)
    print("\nlambda =", [round(float(l), 2) for l in lam],
          " kappa*Omega/lambda1 =", round(KAP * OM / lam1, 4))
    print("VERDICT:", res["VERDICT"])
    print("written:", os.path.normpath(OUT))


if __name__ == "__main__":
    main()
