"""Symbolic check of the convention of t:grass-change (Section 3): rows of the Jacobian supermatrix
are the old coordinates, left derivatives, Ber J = det(A - B D^{-1} C)/det D.  The transposed
arrangement is shown to fail."""
# Check Berezin's change of variables with J_{ab} = (left) derivative of new coordinate b with respect to old coordinate a,
import sys, os; sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
# Ber J = det(A - B D^{-1} C)/det D, A = dy/dx, B = (d eta_b/dx), C = (dy/d xi_a), D = (d eta_b/ d xi_a).
import sympy as sp
from grassmann import G, func_even, det, inv_even
x, Y = sp.symbols('x Y', real=True)
xi1, xi2 = G.gen(1), G.gen(2)
# change of variables: body map x -> 2x + x^3 (diffeo of R), nilpotent shifts, odd mixing depending on x
y = G.const(2 * x + x**3) + sp.sin(x) * xi1 * xi2
eta1 = (1 + x**2) * xi1 + sp.cos(x) * xi2
eta2 = x * xi1 + 3 * xi2
# Jacobian blocks
A = y.apply(lambda e: sp.diff(e, x))
B = [eta1.apply(lambda e: sp.diff(e, x)), eta2.apply(lambda e: sp.diff(e, x))]     # row (1 x 2), odd
C = [y.dxi(1), y.dxi(2)]                                                           # column (2 x 1), odd
D = [[eta1.dxi(1), eta2.dxi(1)], [eta1.dxi(2), eta2.dxi(2)]]                       # D[a][b] = d eta_b / d xi_a
detD = det(D)
idet = inv_even(detD)
Dinv = [[D[1][1] * idet, -D[0][1] * idet], [-D[1][0] * idet, D[0][0] * idet]]
BDC = G()
for a in range(2):
    for b in range(2):
        BDC = BDC + B[a] * Dinv[a][b] * C[b]
Ber = (A - BDC) * idet
Ber = Ber.simplify()
print('Ber =', Ber.d)
# test functions f(Y, eta) = g0(Y) + g1(Y) eta1 + g2(Y) eta2 + g12(Y) eta1 eta2 with Gaussian decay
g0 = sp.exp(-Y**2) * (1 + Y); g12 = sp.exp(-Y**2 / 2) * Y**2; g1 = sp.exp(-Y**2); g2 = Y * sp.exp(-Y**2)
lhs = sp.integrate(g12, (Y, -sp.oo, sp.oo))   # Berezin: coefficient of eta1 eta2
F = func_even(g0, y, Y) + func_even(g1, y, Y) * eta1 + func_even(g2, y, Y) * eta2 + func_even(g12, y, Y) * eta1 * eta2
top = (F * Ber).coeff((1, 2))
top = sp.simplify(top)
rhs = sp.Integral(top, (x, -sp.oo, sp.oo)).evalf(30)
print('lhs', sp.N(lhs, 20), 'rhs', rhs)
# with the transposed (rows = new coordinates, left derivatives) arrangement the result should fail
BDC2 = G()
Bt = [y.dxi(1), y.dxi(2)]; Ct = [eta1.apply(lambda e: sp.diff(e, x)), eta2.apply(lambda e: sp.diff(e, x))]
Dt = [[D[0][0], D[1][0]], [D[0][1], D[1][1]]]
Dtinv = [[Dt[1][1] * idet, -Dt[0][1] * idet], [-Dt[1][0] * idet, Dt[0][0] * idet]]
for a in range(2):
    for b in range(2):
        BDC2 = BDC2 + Bt[a] * Dtinv[a][b] * Ct[b]
Ber2 = ((A - BDC2) * idet).simplify()
top2 = sp.simplify((F * Ber2).coeff((1, 2)))
print('transposed arrangement:', sp.Integral(top2, (x, -sp.oo, sp.oo)).evalf(30))