"""Check of the cross identities behind the harmonicity of the tubular half-density (Section 16).

For F = Ber(L_S), L_S A = S_+ A - A S_- on Hom(V_-, V_+), and the flat supertrace Laplacians
Delta_+ and Delta_- of the two normal blocks, the proof in Section 16 shows

    Delta_+ log F = - |grad_- log F|^2,      Delta_- log F = - |grad_+ log F|^2,

so that F^{-1} (Delta_+ + Delta_-) F = Delta log F + |grad log F|^2 = 0.
Equivalently, F Delta_+ F = |grad_+ F|^2 - |grad_- F|^2 and F Delta_- F = |grad_- F|^2 - |grad_+ F|^2.
Here |grad_+ f|^2 = (Delta_+(f^2) - 2 f Delta_+ f)/2.  The script checks both identities numerically for every Grassmann coefficient at random
complex base points, using the finite Grassmann jet algebra with floating-point complex
coefficients of tube_super_check.py.
Run: python3 -I tube_cross_identity.py   (from this directory, or with this directory on sys.path).
"""
import os
import random
import sys

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import tube_super_check as sc  # noqa: E402


def block_laplacian(f, evars, ovars, par, blk, kappa):
    """Flat supertrace Laplacian of one normal block (blk = '+' or '-')."""
    res = sc.SA.const(0)
    eidx = {key: k for k, key in enumerate(evars)}
    oidx = {key: k for k, key in enumerate(ovars)}
    n = len(par)
    for i in range(n):
        for j in range(n):
            if par[i] == par[j]:
                sgn = 1 if par[i] == 0 else -1
                res = res + f.d_even(eidx[(blk, j, i)]).d_even(eidx[(blk, i, j)]) * sgn
            elif par[i] == 0:
                res = res + f.d_odd(oidx[(blk, i, j)]).d_odd(oidx[(blk, j, i)]) * kappa
    return res


def check(par_p, par_m, x0, s=0.75, kappa=-2.0):
    a, b = -1j * s, 1j * s
    evars, ovars, ent, Sp, Sm = sc.build(par_p, par_m, a, b, x0)
    M, parH = sc.L_matrix(Sp, Sm, par_p, par_m)
    F = sc.ber(M, parH)
    Lp = lambda f: block_laplacian(f, evars, ovars, par_p, '+', kappa)
    Lm = lambda f: block_laplacian(f, evars, ovars, par_m, '-', kappa)
    LpF, LmF = Lp(F), Lm(F)
    gp = (Lp(F * F) - F * LpF * 2) * 0.5
    gm = (Lm(F * F) - F * LmF * 2) * 0.5
    e_plus = (F * LpF - (gp - gm)).at_base().maxabs()
    e_minus = (F * LmF - (gm - gp)).at_base().maxabs()
    size = max((gp.at_base()).maxabs(), (gm.at_base()).maxabs())
    print(f'V+ = {par_p.count(0)}|{par_p.count(1)}, V- = {par_m.count(0)}|{par_m.count(1)}: '
          f'size of |grad F|^2 terms {size:.3g}; '
          f'max coeff of F Delta_+ F - (|grad_+F|^2 - |grad_-F|^2) = {e_plus:.2g}; '
          f'of F Delta_- F - (|grad_-F|^2 - |grad_+F|^2) = {e_minus:.2g}')
    return max(e_plus, e_minus)


if __name__ == '__main__':
    random.seed(7)
    worst = 0.0
    for pp, pm in [([0, 1], [0, 1]), ([0], [0, 1, 1]), ([0, 1, 1], [0]), ([0, 0, 1], [0]),
                   ([0, 0], [0, 1]), ([0, 1], [0, 0, 1])]:
        evars, ovars = sc.setup(pp, pm)
        x0 = {key: complex(random.uniform(-.2, .2), random.uniform(-.2, .2)) for key in evars}
        worst = max(worst, check(pp, pm, x0))
    print(f'worst coefficient error: {worst:.2g}')