"""LAB01: exact three-dimensional Clifford and graded Fierz calculations.

Run: python spinor_algebra.py
Change one input: python spinor_algebra.py --flip-orientation
Requires SymPy 1.13.3; no network, data files, random sampling or plotting.
Conventions: arXiv:1212.3388v3 Appendix A.1; canonical S03 and T09 lessons.
"""
import argparse
import itertools
import json
import sympy as s


def add(*polys):
    result = {}
    for poly in polys:
        for mask, coefficient in poly.items():
            result[mask] = s.expand(result.get(mask, 0)+coefficient)
    return {m:c for m,c in result.items() if c != 0}


def scale(number, poly):
    return add({m:number*c for m,c in poly.items()})


def wedge(left, right):
    """Masks encode increasing products theta1 ... theta8.

    A repeated bit gives zero. Count crossings between a left generator
    and every lower-index right generator to restore increasing order.
    Coefficients are ordinary commuting SymPy expressions.
    """
    out = {}
    for a, ca in left.items():
        for b, cb in right.items():
            if a & b:
                continue
            crossings = sum((b & ((1 << i)-1)).bit_count()
                            for i in range(8) if a & (1 << i))
            mask = a | b
            out[mask] = out.get(mask, 0)+(-1)**crossings*ca*cb
    return add(out)


def bilinear(u, matrix, v):
    return add(*(scale(matrix[a,b], wedge(u[a],v[b]))
                 for a in range(2) for b in range(2)))


def run(flip=False):
    gamma = [s.diag(1,-1), s.Matrix([[0,-1],[-1,0]]),
             s.Matrix([[0,s.I],[-s.I,0]])]
    if flip:
        gamma[2] = -gamma[2]
    C = s.Matrix([[0,-1],[1,0]])
    eye = s.eye(2)
    checks = {}

    def verify(group, condition):
        if not condition:
            raise AssertionError(group)
        checks[group] = checks.get(group, 0)+1

    for a,b in itertools.product(range(3),repeat=2):
        verify('clifford', gamma[a]*gamma[b]+gamma[b]*gamma[a] == 2*int(a==b)*eye)
    oriented = gamma[0]*gamma[1]*gamma[2] == s.I*eye
    if flip:
        # This intentional changed-input run diagnoses the mismatch.
        assert not oriented
        return {'clifford_checks':checks['clifford'],
                'orientation_matches_epsilon123_plus_one':False,
                'triple_product':'-i I',
                'interpretation':'Clifford alone does not fix orientation.'}
    verify('orientation', oriented)
    for a,b in itertools.product(range(3),repeat=2):
        verify('orientation', gamma[a]*gamma[b] == int(a==b)*eye+s.I*sum(
            (s.LeviCivita(a,b,c)*gamma[c] for c in range(3)),s.zeros(2)))
    # Independent component tensor form of completeness: source (A.4).
    lowered = [g*C.T for g in gamma]
    for a,b,c,d in itertools.product(range(2),repeat=4):
        verify('lowered_gamma_completeness', sum(g[a,b]*g[c,d] for g in lowered)
               == -C[a,c]*C[b,d]-C[a,d]*C[b,c])
    generators = [{1<<i:s.Integer(1)} for i in range(8)]
    for a,b in itertools.product(range(8),repeat=2):
        verify('grassmann_anticommutation', not add(wedge(generators[a],generators[b]),
                                                   wedge(generators[b],generators[a])))
    for a,b,c in itertools.product(range(8),repeat=3):
        verify('grassmann_associativity',
               wedge(wedge(generators[a],generators[b]),generators[c]) ==
               wedge(generators[a],wedge(generators[b],generators[c])))
    for parities in itertools.product([0,1],repeat=4):
        spinors = []
        for j,parity in enumerate(parities):
            spinors.append(generators[2*j:2*j+2] if parity else
                           [{0:s.Symbol(f'z{j}_{a}')} for a in range(2)])
        u,v,w,x = spinors
        scalar = lambda a,b: bilinear(a,C,b)
        vector = lambda a,b: [bilinear(a,C*g,b) for g in gamma]
        contract = lambda a,b: add(*(wedge(aa,bb) for aa,bb in zip(a,b)))
        S = wedge(scalar(u,v),scalar(w,x))
        V = contract(vector(u,v),vector(w,x))
        Sp = wedge(scalar(u,x),scalar(w,v))
        Vp = contract(vector(u,x),vector(w,v))
        pv,pw,px = parities[1:]
        sign = (-1)**(pv*pw+px*(pv+pw))
        verify('graded_fierz', not add(S,scale(-s.Rational(sign,2),add(Sp,Vp))))
        verify('graded_fierz', not add(V,scale(-s.Rational(sign,2),add(scale(3,Sp),scale(-1,Vp)))))
    psi,chi = generators[:2],generators[2:4]
    scalar = bilinear(psi,C,chi)
    vector = [bilinear(psi,C*g,chi) for g in gamma]
    scalar2 = wedge(scalar,scalar)
    vector2 = add(*(wedge(j,j) for j in vector))
    verify('quartic_example',scalar2 == {15:-2})
    verify('quartic_example',vector2 == {15:6})
    verify('quartic_example',not add(vector2,scale(3,scalar2)))
    # A Hermitian commuting projector is a different operation.
    zeta = s.Matrix([1, 1+2*s.I])
    n = s.simplify((zeta.conjugate().T*zeta)[0])
    eta = s.Matrix([s.simplify((zeta.conjugate().T*g*zeta)[0]/n) for g in gamma])
    rho = s.simplify(zeta*zeta.conjugate().T/n)
    verify('commuting_projector',s.simplify(rho*rho-rho) == s.zeros(2))
    verify('commuting_projector',eta.dot(eta) == 1)
    verify('commuting_projector',s.simplify(rho-(eye+sum((eta[a]*gamma[a] for a in range(3)),s.zeros(2)))/2) == s.zeros(2))
    P = s.Matrix([(zeta.T*C*g*zeta)[0] for g in gamma])
    verify('commuting_projector',s.simplify(P.dot(P)) == 0)
    return {'arithmetic':'exact integers, rationals and symbolic polynomials',
            'checks':checks,'total':sum(checks.values()),'parity_assignments':16,
            'coefficient_of_theta1_theta2_theta3_theta4':{'scalar_squared':-2,'vector_squared':6},
            'commuting_spinor':[str(z) for z in zeta], 'unit_direction':[str(z) for z in eta],
            'scope':'Pointwise algebra; no Killing equation or quantum composite renormalization.'}


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--flip-orientation',action='store_true',help='Reverse gamma3 only and diagnose the orientation mismatch.')
    args = parser.parse_args()
    print(json.dumps(run(args.flip_orientation),indent=2))