"""Resolve the open degree-3 boxes near the face F = {th12+th13+th23 = 2 pi} (pattern 33,33,34).
Step 1 (monotone push): for each open box B, certify d/dt G(x+t n) > 0, n=(1,1,1), on the prism swept
  from B to F (affine-arithmetic gradient enclosures on sub-boxes; no integrals). Moving along n keeps
  the triangle inequalities, so G(x) <= G(x*) with x* = x + t* n on F.
Step 2 (face bound): cover the face patch {x*: x in B} by parallelogram cells in coordinates
  (a,b) = (th12, th13), th23 = 2 pi - a - b. Each cell: G <= G(x1) + max over the cell's 4 vertices of
  gc.(y-x1) + sum_m gr_m max|y_m-x1_m|, with x1 strictly inside D near the cell, the gradient enclosed over
  the box hull of the cell and x1, and G(x1) in Arb. Cells are subdivided adaptively.
Usage: python deep3.py  (reads open_orig.json; writes deep3_result.json)"""
import json, sys, time, itertools, numpy as np
import verify as v
v.init()
import grad3aa as ga, grad as g0
from flint import arb
TWO_PI=2*np.pi; PAT=['33','33','34']; V=0.1529003
rs=v.rvals(PAT)
d=json.load(open("leaves_33,33,34_0.15290030000000002.json"))
items=[("E",e) for e in d["empties"]]+[("L",l) for l in d["leaves"]]
openk=json.load(open("open_orig.json"))
def grad_enc(lo,hi):
    a=ga.grad_box_aa(rs,lo,hi); b=g0.grad_box(rs,lo,hi)
    out=[]
    for x,y in zip(a,b):
        if x.is_finite() and y.is_finite():
            try: out.append(x.intersection(y))
            except Exception: out.append(x)
        elif x.is_finite(): out.append(x)
        elif y.is_finite(): out.append(y)
        else: out.append(None)
    return out
# ---------- step 1
def push_ok(lo,hi,depth=0,stats=None):
    lo=list(lo); hi=list(hi)
    if sum(lo)>TWO_PI: return True                      # no point of D strictly below F here
    gg=grad_enc(lo,hi)
    if all(x is not None for x in gg):
        sdir=gg[0]+gg[1]+gg[2]
        if sdir>0: return True
    if depth>=14: stats["fail"].append((lo,hi)); return False
    m=int(np.argmax([h-l for l,h in zip(lo,hi)])); mid=(lo[m]+hi[m])/2
    h1=list(hi); h1[m]=mid; l2=list(lo); l2[m]=mid
    a=push_ok(lo,h1,depth+1,stats); b=push_ok(l2,hi,depth+1,stats); return a and b
# ---------- step 2
def face_pt(a,b): return [a,b,TWO_PI-a-b]
def cell_ok(alo,ahi,blo,bhi):
    verts=[face_pt(a,b) for a in (alo,ahi) for b in (blo,bhi)]
    ac,bc=(alo+ahi)/2,(blo+bhi)/2; delta=1e-7
    x1=[t-delta for t in face_pt(ac,bc)]
    if not v.strictly_inD(x1): return False,None
    lo=[min([p[m] for p in verts]+[x1[m]]) for m in range(3)]; hi=[max([p[m] for p in verts]+[x1[m]]) for m in range(3)]
    gg=grad_enc(lo,hi)
    if any(x is None for x in gg): return False,None
    gc=[float(x.mid()) for x in gg]; gr=[abs(x-arb(c)).upper() for x,c in zip(gg,gc)]
    lin=None
    for p in verts:
        t=sum((arb(gc[m])*(arb(p[m])-arb(x1[m])) for m in range(3)),arb(0))
        lin=t if lin is None else lin.union(t)
    err=sum((gr[m]*max((abs(arb(p[m])-arb(x1[m])) for p in verts),key=lambda z:float(z.upper())).upper() for m in range(3)),arb(0))
    # cheap float screen before the Arb integral
    est=float(lin.upper())+float(err.upper())
    for tol in (arb(2)**-30,arb(2)**-40):
        Gv=v.G_rig(rs,x1,tol)
        if not Gv.is_finite(): continue
        ub=Gv+lin.upper()+err
        if ub<arb(V): return True,float(ub.upper())
    return False,float(ub.upper()) if Gv.is_finite() else None
def face_cover(rect,stats,depth=0):
    alo,ahi,blo,bhi=rect
    # skip cells whose face points all violate a triangle inequality or exceed pi (outside D)
    ok,ub=cell_ok(alo,ahi,blo,bhi)
    if ok: stats["cells"]+=1; stats["max_ub"]=max(stats["max_ub"],ub); return True
    if depth>=16: stats["face_fail"].append((rect,ub)); return False
    if ahi-alo>=bhi-blo: m=(alo+ahi)/2; return face_cover((alo,m,blo,bhi),stats,depth+1) & face_cover((m,ahi,blo,bhi),stats,depth+1)
    m=(blo+bhi)/2; return face_cover((alo,ahi,blo,m),stats,depth+1) & face_cover((alo,ahi,m,bhi),stats,depth+1)
if __name__=="__main__":
    t0=time.time(); res=dict(boxes=len(openk),push={},face=None)
    rects=[]
    for k in openk:
        o=items[k][1]; lo,hi=o["lo"],o["hi"]
        tmax=max(0.0,(TWO_PI-sum(lo))/3)
        elo=list(lo); ehi=[h+tmax for h in hi]
        st=dict(fail=[]); ok=push_ok(elo,ehi,0,st)
        res["push"][k]=dict(ok=ok,fail=len(st["fail"])); print("push",k,ok,len(st["fail"]),f"{time.time()-t0:.0f}s",flush=True)
        # face patch bounding rectangle in (a,b)
        alo=(2*lo[0]-hi[1]-hi[2])/3+TWO_PI/3; ahi=(2*hi[0]-lo[1]-lo[2])/3+TWO_PI/3
        blo=(2*lo[1]-hi[0]-hi[2])/3+TWO_PI/3; bhi=(2*hi[1]-lo[0]-lo[2])/3+TWO_PI/3
        rects.append((alo,ahi,blo,bhi))
    A0=min(r[0] for r in rects); A1=max(r[1] for r in rects); B0=min(r[2] for r in rects); B1=max(r[3] for r in rects)
    print("face rectangle",A0,A1,B0,B1,flush=True)
    st=dict(cells=0,max_ub=-1.0,face_fail=[])
    okf=face_cover((A0,A1,B0,B1),st)
    res["face"]=dict(ok=okf,rect=(A0,A1,B0,B1),cells=st["cells"],max_ub=st["max_ub"],fail=st["face_fail"][:20],nfail=len(st["face_fail"]))
    res["seconds"]=round(time.time()-t0)
    json.dump(res,open("deep3_result.json","w"),default=str)
    print("RESULT push all ok:",all(p["ok"] for p in res["push"].values()),"face ok:",okf,"cells",st["cells"],"max_ub",st["max_ub"],"fails",len(st["face_fail"]),f"{time.time()-t0:.0f}s",flush=True)