1import json, time
  2import numpy as np
  3from scipy.optimize import nnls
  4
  5
  6def certify(a, b, retained, tol=1e-9):
  7    A = np.asarray(a, float); B = np.asarray(b, float)
  8    R = list(retained); out = {}
  9    for j in range(len(A)):
 10        if j in R: continue
 11        lam, _ = nnls(A[R].T, A[j])
 12        resid = np.linalg.norm(A[j] - A[R].T @ lam)
 13        offset = float(lam @ B[R] - B[j])
 14        out[j] = (resid <= tol and offset <= tol, resid, offset, lam)
 15    return out
 16
 17
 18def reduce_constraints(a, b, tol=1e-8):
 19    A=np.asarray(a,float); B=np.asarray(b,float)
 20    ang=np.mod(np.arctan2(A[:,1], A[:,0]),2*np.pi)
 21    order=np.argsort(ang)
 22    gaps=np.diff(np.r_[ang[order],ang[order[0]]+2*np.pi])
 23    if gaps.max() > np.pi + tol:
 24        k=int(np.argmax(gaps)); R=[int(order[k]),int(order[(k+1)%len(order)])]
 25    else:
 26        R=[int(order[0]),int(order[len(order)//3]),int(order[2*len(order)//3])]
 27    R=list(dict.fromkeys(R))
 28    while True:
 29        c=certify(A,B,R,tol)
 30        bad=[j for j,v in c.items() if not v[0]]
 31        if not bad: return sorted(R), c
 32        R.append(bad[0])
 33
 34
 35def feasible(x,a,b,tol=1e-8):
 36    return bool(np.all(np.asarray(a)@x <= np.asarray(b)+tol))
 37
 38
 39def project(x0,a,b):
 40    """2-D Euclidean projection onto intersection of halfspaces."""
 41    A=np.asarray(a,float); B=np.asarray(b,float); x=np.asarray(x0,float)
 42    if feasible(x,A,B): return x.copy()
 43    best=None; bestv=np.inf
 44    # Candidates are projections onto one boundary and intersections of two.
 45    for i in range(len(A)):
 46        den=A[i]@A[i]; y=x-A[i]*(A[i]@x-B[i])/den
 47        if feasible(y,A,B,1e-7) and np.sum((y-x)**2)<bestv: best,bestv=y,np.sum((y-x)**2)
 48    for i in range(len(A)):
 49        for j in range(i):
 50            M=np.vstack((A[i],A[j]))
 51            if abs(np.linalg.det(M))<1e-12: continue
 52            y=np.linalg.solve(M,np.array([B[i],B[j]]))
 53            if feasible(y,A,B,1e-7) and np.sum((y-x)**2)<bestv: best,bestv=y,np.sum((y-x)**2)
 54    if best is None: raise RuntimeError('infeasible or numerical projection failure')
 55    return best
 56
 57
 58def make_cone(m, spread=1.2, margin=0.05, rng=None):
 59    rng=np.random.default_rng(0) if rng is None else rng
 60    theta=np.linspace(-spread,spread,m)+rng.normal(0,spread/(10*m),m)
 61    A=np.c_[np.cos(theta),np.sin(theta)]
 62    # b=1 for extreme rows; interior rows are deliberately redundant.
 63    b=np.ones(m)*3.0
 64    b[0]=b[-1]=1.0
 65    return A,b
 66
 67
 68def make_circle(m, margin=.05):
 69    th=np.linspace(0,2*np.pi,m,endpoint=False)
 70    return np.c_[np.cos(th),np.sin(th)],np.ones(m)*(1+margin)
 71
 72
 73def run(seed=7):
 74    rng=np.random.default_rng(seed); rows=[]
 75    # Prediction 1: adding constraints inside a pointed cone does not change exact set.
 76    for m in [4,8,16,32,64]:
 77        A,B=make_cone(m,rng=rng); R,c=reduce_constraints(A,B)
 78        errs=[]; agree=0
 79        for _ in range(50):
 80            x=rng.uniform(-2,2,2); yf=project(x,A,B); yr=project(x,A[R],B[R])
 81            errs.append(np.linalg.norm(yf-yr)); agree += feasible(x,A,B)==feasible(x,A[R],B[R])
 82        rows.append({'test':'cone_scaling','m':m,'retained':len(R),'fraction':len(R)/m,'max_projection_error':max(errs),'decision_agreement':agree/50})
 83    # Prediction 2: crossing spread=pi/2 loses the two-extreme initialization.
 84    for spread in [0.5, 1.2, 1.55, 1.7, 2.4]:
 85        A,B=make_cone(32,spread=spread,rng=rng); R,c=reduce_constraints(A,B)
 86        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])})
 87    # Prediction 3: full-circle normals cannot be represented by two extremes; all facets remain.
 88    for m in [8,16,32]:
 89        A,B=make_circle(m); R,c=reduce_constraints(A,B)
 90        rows.append({'test':'angular_transition','m':m,'retained':len(R),'fraction':len(R)/m,'all_certified':all(v[0] for v in c.values())})
 91    # Prediction 4: nominal feasible point is accepted identically, infeasible is projected identically.
 92    A,B=make_cone(64,rng=rng); R,c=reduce_constraints(A,B)
 93    nominal=np.array([0.,0.]); xbad=np.array([2.,.7])
 94    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())})
 95    # Direct containment check: random points satisfying R must satisfy every omitted row.
 96    A,B=make_cone(48,rng=rng); R,c=reduce_constraints(A,B); safe=0; false_safe=0
 97    for _ in range(10000):
 98        x=rng.uniform(-4,4,2)
 99        if feasible(x,A[R],B[R]):
100            safe += 1; false_safe += int(not feasible(x,A,B))
101    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])})
102    # Runtime comparison, same projection routine and inputs.
103    timing=[]
104    for m in [16,32,64,128]:
105        A,B=make_cone(m,rng=rng); R,_=reduce_constraints(A,B); xs=rng.uniform(-2,2,(20,2))
106        t=time.perf_counter(); [project(x,A,B) for x in xs]; tf=time.perf_counter()-t
107        t=time.perf_counter(); [project(x,A[R],B[R]) for x in xs]; tr=time.perf_counter()-t
108        timing.append({'m':m,'retained':len(R),'full_sec':tf,'reduced_sec':tr,'speedup':tf/tr})
109    return {'rows':rows,'timing':timing}
110
111if __name__=='__main__':
112    result=run(); print(json.dumps(result,indent=2))