#!/usr/bin/env python3
"""Independent Python checks for the OMEGA geometry/operator bridge.

This is not an execution of the Julia file or a proof-assistant certificate.
Dependencies: NumPy, SymPy (their actual versions print in the execution log).
Run: python verify_reference.py ; python -O verify_reference.py
No Python assert statements are used. Failures raise explicit exceptions.
"""
from __future__ import annotations
import ast
import math
import sys
from pathlib import Path
from typing import Callable
import numpy as np
import sympy as sp

GROUPS: list[str] = []
CHECKS = 0

def require(ok: bool, message: str) -> None:
    global CHECKS
    CHECKS += 1
    if not bool(ok):
        raise ValueError(message)

def close(a, b, *, atol: float = 2e-12, rtol: float = 2e-11) -> bool:
    return bool(np.linalg.norm(np.asarray(a)-np.asarray(b)) <=
                atol+rtol*max(np.linalg.norm(a),np.linalg.norm(b)))

def group(name: str, fn: Callable[[], None]) -> None:
    fn()
    GROUPS.append(name)
    print(f"PASS {len(GROUPS):02d} {name}")

def zero(x: sp.Expr, message: str) -> None:
    require(sp.simplify(sp.expand(x)) == 0, message)

def geometry(a: float = 1.0, epsilon: float = math.pi/648000) -> dict:
    require(math.isfinite(a) and a > 0, "invalid scale")
    require(math.isfinite(epsilon) and 0 < epsilon < math.pi, "invalid seed")
    phi = (1+math.sqrt(5))/2
    s,c = math.sin(epsilon/2),math.cos(epsilon/2)
    z = np.array([-phi**2-s,-1-s,-s,s,1+s,phi**2+s])
    X = a*np.vstack((z,np.full(6,c)))
    R,r = a*np.linalg.norm(z)/math.sqrt(6),a*c
    return dict(a=a,epsilon=epsilon,phi=phi,s=s,c=c,z=z,X=X,G=X.T@X,
                R=R,r=r,bearings=np.arctan2(z,c))

def operator(R: float,r: float,center: complex = 0j,angle: float = 0.0) -> np.ndarray:
    require(all(math.isfinite(t) for t in (R,r,center.real,center.imag,angle))
            and R >= r > 0,"invalid operator parameters")
    d=math.sqrt(R-r)*math.sqrt(R+r)
    return center*np.eye(2,dtype=complex)+np.exp(1j*angle)*np.array([[d,2*r],[0,-d]])

def boundary(T: np.ndarray, theta: float):
    U=np.exp(-1j*theta)*T
    vals,vecs=np.linalg.eigh((U+U.conj().T)/2)
    v=vecs[:,-1]
    return np.vdot(v,T@v),vals[-1],v

def witness(R: float,r: float,w: complex) -> np.ndarray:
    d=math.sqrt(R-r)*math.sqrt(R+r)
    x,y=w.real,w.imag
    t=(x/R)**2+(y/r)**2
    require(t<=1+256*np.finfo(float).eps,"outside ellipse")
    X,Y,Z=np.array([r*x/R**2,y/r,d*x/R**2])+math.sqrt(max(0,1-t))*np.array([d/R,0,-r/R])
    if Z>=0:
        v1=math.sqrt((1+Z)/2)
        return np.array([v1,(X+1j*Y)/(2*v1)])
    v2=math.sqrt((1-Z)/2)
    return np.array([(X-1j*Y)/(2*v2),v2])

def blocks(R: float,r: float,n: int,center: complex=0j,angle: float=0.0):
    require(1<=n<=512,"invalid block count")
    B=center*np.eye(2*n,dtype=complex)
    T=operator(R,r)
    for j in range(n):
        B[2*j:2*j+2,2*j:2*j+2]+=np.exp(1j*angle)*math.ldexp(1.0,-j)*T
    return B

def pval(coeffs,z):
    ans=np.zeros_like(coeffs[0],dtype=complex)
    for C in reversed(coeffs):
        ans=ans*z+C
    return ans

def pmat(coeffs,T):
    ans=np.zeros((len(coeffs[0])*len(T),)*2,dtype=complex)
    power=np.eye(len(T),dtype=complex)
    for C in coeffs:
        ans+=np.kron(C,power)
        power=power@T
    return ans

def opnorm(T):
    return float(np.linalg.norm(T,2))

def main() -> None:
    print("OMEGA geometry/operator bridge: independent Python reference v1.0")
    print("Python",sys.version.split()[0],"NumPy",np.__version__,"SymPy",sp.__version__)
    print("Optimization mode:",sys.flags.optimize)
    print("Julia execution: NOT RUN. This is a separate implementation.")
    t=sp.symbols('s',real=True);c=sp.symbols('c',real=True)
    phi=(1+sp.sqrt(5))/2
    zs=sp.Matrix([-phi**2-t,-1-t,-t,t,1+t,phi**2+t])
    ones=sp.ones(6,1)
    X=sp.Matrix.vstack(zs.T,c*ones.T)
    rsq=t*t+sp.Rational(2,3)*(phi+2)*t+phi**2
    def golden():
        zero(phi**2-phi-1,"golden")
        zero(1+phi**4-3*phi**2,"fourth power")
    group("exact golden and fourth-power identities",golden)
    def gram():
        G=X.T*X-c*c*(ones*ones.T)-zs*zs.T
        for x in G: zero(x,"Gram entry")
        zero(sum(zs),"zero signed sum")
    group("exact 36-entry Gram factorization and signed sum",gram)
    def moment():
        M=X*X.T/6-sp.diag(rsq,c*c)
        for x in M: zero(x,"moment entry")
        zero(rsq-(1-t*t)-(2*t*t+sp.Rational(2,3)*(phi+2)*t+phi),"focal polynomial")
    group("exact moment and focal polynomials",moment)
    def bloch():
        r,d=sp.symbols('r d',positive=True)
        L=sp.Matrix([[r,0,d],[0,r,0]])
        for x in L*L.T-sp.diag(r*r+d*d,r*r):zero(x,"Bloch ellipse")
        T=sp.Matrix([[d,2*r],[0,-d]])
        for x in T*T-d*d*sp.eye(2):zero(x,"T square")
        q=sp.symbols('q')
        zero((T.T*T).charpoly(q).as_expr()-(q*q-(2*d*d+4*r*r)*q+d**4),"norm roots")
    group("exact Bloch projection and operator polynomial",bloch)
    def recovery_identity():
        q,h=sp.symbols('q h')
        zero((q+1)*t*t+h*t+phi**2-q-
             ((t*t+h*t+phi**2)-q*(1-t*t)),"inverse quadratic")
        # Positive derivative numerator for q(s) on 0 <= s < 1.
        h=sp.Rational(2,3)*(phi+2)
        expected=h+2*(1+phi**2)*t+h*t*t
        zero(sp.diff(rsq/(1-t*t),t)*(1-t*t)**2-expected,"aspect derivative")
    group("exact inverse-family and monotonicity identities",recovery_identity)

    k=geometry();R,r=k['R'],k['r'];T=operator(R,r)
    max_errors={}
    def source_nodes():
        # Selected values transcribed from the inspected website JSON.
        expected=np.array([-2.6180364128183005,-1.0000024240684056,
            -2.424068405545306e-6,2.424068405545306e-6,
            1.0000024240684056,2.6180364128183005])
        require(close(k['z'],expected,atol=2e-15,rtol=2e-15),"source x")
        require(close(k['c'],0.9999999999970619,atol=2e-15,rtol=2e-15),"source y")
        max_errors['source_node_max_abs']=float(np.max(np.abs(k['z']-expected)))
    group("published unit-node transcription",source_nodes)
    def full_curve():
        for j in (2,3):require(close(np.linalg.norm(k['X'][:,j]),1),"seed radius")
        require(close(np.diff(k['z']),[k['phi'],1,2*k['s'],1,k['phi']]),"gaps")
        for angle,j in ((k['epsilon']/2,3),(2*math.pi-k['epsilon']/2,2)):
            require(close([math.sin(angle),math.cos(angle)],k['X'][:,j]),"arc endpoint")
        for sign,j in ((-1,0),(1,5)):
            require(close([sign*(k['s']+k['phi']**2),k['c']],k['X'][:,j]),"foot")
    group("seed radii, retained-arc endpoints, feet and gaps",full_curve)
    def angles():
        G=k['G'];betas=k['bearings']
        for i in range(6):
            for j in range(i+1,6):
                require(close(math.cos(betas[j]-betas[i]),G[i,j]/math.sqrt(G[i,i]*G[j,j])),"angle")
        require(close(betas[3],k['epsilon']/2,atol=2e-18,rtol=2e-12),"theta0")
        require(close(betas[4],math.pi/4+k['epsilon']/4),"theta1")
    group("all fifteen center angles and exact-form bearing identities",angles)
    def ranks():
        require(np.linalg.matrix_rank(k['X'])==2,"circle-centered rank")
        require(np.linalg.matrix_rank(k['X']-k['X'].mean(axis=1,keepdims=True))==1,"centroid rank")
    group("circle-centered rank two versus centroid-centered rank one",ranks)
    def placement():
        E=np.zeros((9,2));E[[0,2],0]=1/math.sqrt(2)
        E[1,1]=1/math.sqrt(2);E[8,1]=-1/math.sqrt(2)
        A=np.arange(9).reshape(-1,1)/10
        Y=np.column_stack((k['X'],[0,-1],[-(k['s']+0.3),k['c']]))
        P=A+E@Y
        require(close((P-A).T@(P-A),Y.T@Y),"ambient Gram")
    group("isometric placement of marks, arc and foot samples",placement)
    def ball_image():
        L=k['X']/math.sqrt(6);M=L@L.T
        require(close(M,np.diag([R*R,r*r])),"ellipse matrix")
        require(close(np.linalg.svd(L,compute_uv=False),[R,r]),"singular axes")
        for theta in np.linspace(0,2*math.pi,129):
            y=np.array([R*math.cos(theta),r*math.sin(theta)])
            q=L.T@np.linalg.solve(M,y)
            require(close(L@q,y) and close(np.linalg.norm(q),1),"constructive image")
    group("ellipse singular values and constructive ball preimages",ball_image)
    def affine_chart():
        D=np.diag([R/k['a'],r/k['a']]);H=np.diag([(k['a']/R)**2,(k['a']/r)**2])
        Y=D@k['X']
        require(close(Y.T@H@Y,k['G']),"metric Gram recovery")
        for theta in np.linspace(k['epsilon']/2,2*math.pi-k['epsilon']/2,129):
            y=D@(k['a']*np.array([math.sin(theta),math.cos(theta)]))
            require(close((y[0]/R)**2+(y[1]/r)**2,1),"ellipse arc image")
        x=D@np.array([1.,1.]);y=D@np.array([1.,-1.])
        require(abs(x@y)>.1,"ordinary angles not invariant")
        require(close(x@H@y,0),"weighted angle recovery")
    group("full-kernel affine ellipse chart and angle-recovery metric",affine_chart)
    def scales():
        for a in (.1,1,3.25,10):
            kk=geometry(a)
            require(close([kk['R'],kk['r']],a*np.array([R,r])),"axis scaling")
            require(close(kk['G'],a*a*k['G']),"Gram scaling")
            require(close(operator(kk['R'],kk['r']),a*T),"operator scaling")
    group("scale covariance of geometry, ellipse and operator",scales)
    def inverse_family():
        err=0.0
        for epsilon in (math.pi/648000,.2,1.2):
            for a in (.3,1,4):
                kk=geometry(a,epsilon);q=(kk['R']/kk['r'])**2;ph=kk['phi']
                h=2*(ph+2)/3
                s0=2*(q-ph**2)/(h+math.sqrt(h*h+4*(q+1)*(q-ph**2)))
                a0=kk['r']/math.sqrt((1-s0)*(1+s0))
                require(close(s0,kk['s'],atol=3e-15,rtol=2e-12),"s inverse")
                require(close(a0,a),"a inverse")
                err=max(err,abs(s0-kk['s']))
        max_errors['inverse_s_max_abs']=err
    group("shape/scale recovery within the declared golden-step family",inverse_family)
    def witnesses():
        err=0.0
        for AA,bb in ((R,r),(3,2),(1,1)):
            A=operator(AA,bb)
            for radius in (0,.2,.8,1):
                for theta in np.linspace(0,2*math.pi,129):
                    w=radius*(AA*math.cos(theta)+1j*bb*math.sin(theta))
                    v=witness(AA,bb,w)
                    require(close(np.linalg.norm(v),1),"unit vector")
                    w2=np.vdot(v,A@v);require(close(w2,w),"numerical witness")
                    err=max(err,abs(w2-w))
        max_errors['witness_max_abs']=float(err)
    group("1548 sampled point-to-unit-vector witnesses",witnesses)
    def supports():
        ew=eh=0.0
        for AA,bb,zeta,angle in ((R,r,0j,0),(5,1,2-3j,.37),(1,1,-1+2j,1.1)):
            A=operator(AA,bb,zeta,angle)
            for theta in np.linspace(0,2*math.pi,721):
                w,h,_=boundary(A,float(theta));u=theta-angle
                q=math.hypot(AA*math.cos(u),bb*math.sin(u))
                e=zeta+np.exp(1j*angle)*(AA*AA*math.cos(u)/q+1j*bb*bb*math.sin(u)/q)
                hs=float(np.real(np.exp(-1j*theta)*zeta))+q
                require(close(w,e) and close(h,hs),"boundary/support")
                ew=max(ew,abs(w-e));eh=max(eh,abs(h-hs))
        max_errors['boundary_max_abs']=float(ew);max_errors['support_max_abs']=float(eh)
    group("2163 support/boundary comparisons including translations",supports)
    def similarity():
        d=math.sqrt(R-r)*math.sqrt(R+r);S=np.array([[1,-r/R],[0,d/R]])
        require(close(T@T,d*d*np.eye(2)),"square")
        require(close(opnorm(T),R+r),"operator norm")
        require(close(T@S,S@np.diag([d,-d])),"diagonalization")
        require(close(np.linalg.cond(S),math.sqrt((R+r)/(R-r))),"similarity constant")
    group("operator square, norm and explicit diagonalization",similarity)
    def regression():
        A=operator(1,1,3+0j);bad=np.zeros((3,3),complex);bad[:2,:2]=A
        good=bad.copy();good[2,2]=3
        require(close(boundary(bad,math.pi)[1],0),"zero tail enlarges range")
        require(close(boundary(good,math.pi)[1],-2),"center tail preserves range")
    group("translated-ellipse zero-tail regression and correction",regression)
    def block_support():
        for n in (1,2,5,9):
            B=blocks(R,r,n,2-3j,.37)
            for theta in np.linspace(0,2*math.pi,65):
                expected=np.real(np.exp(-1j*theta)*(2-3j))+math.hypot(R*math.cos(theta-.37),r*math.sin(theta-.37))
                require(close(boundary(B,float(theta))[1],expected),"block numerical range")
            B1=blocks(R,r,n+1);B0=np.zeros_like(B1);B0[:2*n,:2*n]=blocks(R,r,n)
            require(close(opnorm(B1-B0),math.ldexp(R+r,-n)),"truncation norm")
    group("finite block supports and exact-form tail-error checks",block_support)
    def polynomial():
        coeffs=[np.array([[1,0],[0,2]],complex),np.array([[.2,1j/3],[-.1,.4]]),
                np.array([[.1,-.2],[.3,.05]])]
        d=math.sqrt(R-r)*math.sqrt(R+r);kappa=math.sqrt((R+r)/(R-r))
        for n in (1,3,6):
            A=blocks(R,r,n)
            b=max(opnorm(pval(coeffs,sign*math.ldexp(d,-j))) for j in range(n) for sign in (-1,1))
            require(opnorm(pmat(coeffs,A))<=kappa*b+2e-11,"matrix polynomial bound")
    group("matrix-polynomial similarity bound without a Crouzeix oracle",polynomial)
    def infinite_product_sections():
        # Independent block variations have an N-dimensional Gram matrix.
        # The all-N conclusion is proved in the note, not by this finite test.
        n=12;delta=k['X'][:,5]-k['X'][:,0]
        D=np.zeros((2*n,n))
        for j in range(n):D[2*j:2*j+2,j]=math.ldexp(1.,-j)*delta
        require(np.linalg.matrix_rank(D)==n,"independent product directions")
        weights=np.array([math.ldexp(1.,-j) for j in range(n)])
        require(close(D.T@D,np.diag(weights**2*np.dot(delta,delta))),"product metric")
    group("twelve independent directions in a finite product section",infinite_product_sections)
    def precision():
        # 80-digit evaluation, still not a certified enclosure.
        eps=sp.pi/648000;ss=sp.sin(eps/2);cc=sp.cos(eps/2)
        RR=sp.sqrt((ss**2+(1+ss)**2+(phi**2+ss)**2)/3)
        require(close(float(sp.N(RR,80)),R,atol=2e-15,rtol=2e-15),"high precision R")
        require(close(float(sp.N(cc,80)),r,atol=2e-15,rtol=2e-15),"high precision r")
        print("REFERENCE_80_DIGIT R",sp.N(RR,70))
        print("REFERENCE_80_DIGIT r",sp.N(cc,70))
    group("80-digit canonical-value cross-check",precision)
    def validation():
        for action in (lambda:geometry(-1),lambda:geometry(1,0),
                       lambda:operator(1,2),lambda:operator(1,1,complex(float('inf'),0)),
                       lambda:witness(1,1,2),lambda:blocks(1,1,0)):
            try:action()
            except ValueError:pass
            else:raise ValueError("invalid input accepted")
        tree=ast.parse(Path(__file__).read_text(encoding='utf-8'))
        require(not any(isinstance(n,ast.Assert) for n in ast.walk(tree)),"assert present")
    group("invalid-input rejection and zero Python assert nodes",validation)
    print("GROUPS",len(GROUPS),"PASS; explicit require calls",CHECKS)
    print("CANONICAL phi",format(k['phi'],'.17g'))
    print("CANONICAL s",format(k['s'],'.17g'))
    print("CANONICAL R",format(R,'.17g'),"r",format(r,'.17g'))
    print("CANONICAL d",format(math.sqrt(R-r)*math.sqrt(R+r),'.17g'))
    print("CANONICAL norm(T)",format(opnorm(T),'.17g'))
    print("CANONICAL kappa",format(math.sqrt((R+r)/(R-r)),'.17g'))
    for key,value in max_errors.items():print("MAX_ERROR",key,format(value,'.17g'))
    print("No Julia runtime, Lean build, infinite simulation or physical experiment is claimed.")

if __name__=='__main__':
    main()