import json, time
import numpy as np
from scipy.optimize import nnls


def certify(a, b, retained, tol=1e-9):
    A = np.asarray(a, float); B = np.asarray(b, float)
    R = list(retained); out = {}
    for j in range(len(A)):
        if j in R: continue
        lam, _ = nnls(A[R].T, A[j])
        resid = np.linalg.norm(A[j] - A[R].T @ lam)
        offset = float(lam @ B[R] - B[j])
        out[j] = (resid <= tol and offset <= tol, resid, offset, lam)
    return out


def reduce_constraints(a, b, tol=1e-8):
    A=np.asarray(a,float); B=np.asarray(b,float)
    ang=np.mod(np.arctan2(A[:,1], A[:,0]),2*np.pi)
    order=np.argsort(ang)
    gaps=np.diff(np.r_[ang[order],ang[order[0]]+2*np.pi])
    if gaps.max() > np.pi + tol:
        k=int(np.argmax(gaps)); R=[int(order[k]),int(order[(k+1)%len(order)])]
    else:
        R=[int(order[0]),int(order[len(order)//3]),int(order[2*len(order)//3])]
    R=list(dict.fromkeys(R))
    while True:
        c=certify(A,B,R,tol)
        bad=[j for j,v in c.items() if not v[0]]
        if not bad: return sorted(R), c
        R.append(bad[0])


def feasible(x,a,b,tol=1e-8):
    return bool(np.all(np.asarray(a)@x <= np.asarray(b)+tol))


def project(x0,a,b):
    """2-D Euclidean projection onto intersection of halfspaces."""
    A=np.asarray(a,float); B=np.asarray(b,float); x=np.asarray(x0,float)
    if feasible(x,A,B): return x.copy()
    best=None; bestv=np.inf
    # Candidates are projections onto one boundary and intersections of two.
    for i in range(len(A)):
        den=A[i]@A[i]; y=x-A[i]*(A[i]@x-B[i])/den
        if feasible(y,A,B,1e-7) and np.sum((y-x)**2)<bestv: best,bestv=y,np.sum((y-x)**2)
    for i in range(len(A)):
        for j in range(i):
            M=np.vstack((A[i],A[j]))
            if abs(np.linalg.det(M))<1e-12: continue
            y=np.linalg.solve(M,np.array([B[i],B[j]]))
            if feasible(y,A,B,1e-7) and np.sum((y-x)**2)<bestv: best,bestv=y,np.sum((y-x)**2)
    if best is None: raise RuntimeError('infeasible or numerical projection failure')
    return best


def make_cone(m, spread=1.2, margin=0.05, rng=None):
    rng=np.random.default_rng(0) if rng is None else rng
    theta=np.linspace(-spread,spread,m)+rng.normal(0,spread/(10*m),m)
    A=np.c_[np.cos(theta),np.sin(theta)]
    # b=1 for extreme rows; interior rows are deliberately redundant.
    b=np.ones(m)*3.0
    b[0]=b[-1]=1.0
    return A,b


def make_circle(m, margin=.05):
    th=np.linspace(0,2*np.pi,m,endpoint=False)
    return np.c_[np.cos(th),np.sin(th)],np.ones(m)*(1+margin)


def run(seed=7):
    rng=np.random.default_rng(seed); rows=[]
    # Prediction 1: adding constraints inside a pointed cone does not change exact set.
    for m in [4,8,16,32,64]:
        A,B=make_cone(m,rng=rng); R,c=reduce_constraints(A,B)
        errs=[]; agree=0
        for _ in range(50):
            x=rng.uniform(-2,2,2); yf=project(x,A,B); yr=project(x,A[R],B[R])
            errs.append(np.linalg.norm(yf-yr)); agree += feasible(x,A,B)==feasible(x,A[R],B[R])
        rows.append({'test':'cone_scaling','m':m,'retained':len(R),'fraction':len(R)/m,'max_projection_error':max(errs),'decision_agreement':agree/50})
    # Prediction 2: crossing spread=pi/2 loses the two-extreme initialization.
    for spread in [0.5, 1.2, 1.55, 1.7, 2.4]:
        A,B=make_cone(32,spread=spread,rng=rng); R,c=reduce_constraints(A,B)
        rows.append({'test':'halfplane_transition','spread':spread,'predicted_two_extremes':spread < np.pi/2,'retained':len(R),'fraction':len(R)/32,'max_resid':max([v[1] for v in c.values()] or [0.0])})
    # Prediction 3: full-circle normals cannot be represented by two extremes; all facets remain.
    for m in [8,16,32]:
        A,B=make_circle(m); R,c=reduce_constraints(A,B)
        rows.append({'test':'angular_transition','m':m,'retained':len(R),'fraction':len(R)/m,'all_certified':all(v[0] for v in c.values())})
    # Prediction 4: nominal feasible point is accepted identically, infeasible is projected identically.
    A,B=make_cone(64,rng=rng); R,c=reduce_constraints(A,B)
    nominal=np.array([0.,0.]); xbad=np.array([2.,.7])
    rows.append({'test':'nominal_preservation','nominal_full':feasible(nominal,A,B),'nominal_reduced':feasible(nominal,A[R],B[R]),'bad_projection_error':float(np.linalg.norm(project(xbad,A,B)-project(xbad,A[R],B[R]))),'certified':all(v[0] for v in c.values())})
    # Direct containment check: random points satisfying R must satisfy every omitted row.
    A,B=make_cone(48,rng=rng); R,c=reduce_constraints(A,B); safe=0; false_safe=0
    for _ in range(10000):
        x=rng.uniform(-4,4,2)
        if feasible(x,A[R],B[R]):
            safe += 1; false_safe += int(not feasible(x,A,B))
    rows.append({'test':'farkas_containment','retained':len(R),'retained_feasible_samples':safe,'false_safe_samples':false_safe,'max_certificate_residual':max([v[1] for v in c.values()] or [0.0]),'max_offset_violation':max([v[2] for v in c.values()] or [0.0])})
    # Runtime comparison, same projection routine and inputs.
    timing=[]
    for m in [16,32,64,128]:
        A,B=make_cone(m,rng=rng); R,_=reduce_constraints(A,B); xs=rng.uniform(-2,2,(20,2))
        t=time.perf_counter(); [project(x,A,B) for x in xs]; tf=time.perf_counter()-t
        t=time.perf_counter(); [project(x,A[R],B[R]) for x in xs]; tr=time.perf_counter()-t
        timing.append({'m':m,'retained':len(R),'full_sec':tf,'reduced_sec':tr,'speedup':tf/tr})
    return {'rows':rows,'timing':timing}

if __name__=='__main__':
    result=run(); print(json.dumps(result,indent=2))