#!/usr/bin/env python3 """Exact checks for transport, Hamilton projection, Hopf-Lax and eikonal examples. Only standard-library rational arithmetic; no sampling is used to select minima. Use --output to write a JSON result. Without it, print results only. """ from fractions import Fraction as F from itertools import product from collections import Counter from pathlib import Path import json,argparse def piecewise(knots,slopes,value_at_first=F(0)): """Continuous piecewise affine function, including two unbounded outer pieces.""" knots=tuple(map(F,knots));slopes=tuple(map(F,slopes)) if len(slopes)!=len(knots)+1 or any(a>=b for a,b in zip(knots,knots[1:])):raise ValueError('ordered knots and one slope per piece required') if not knots:raise ValueError('at least one knot required') bs=[F(value_at_first)-slopes[0]*knots[0]] for i,k in enumerate(knots):bs.append(bs[-1]+(slopes[i]-slopes[i+1])*k) return [(None if i==0 else knots[i-1],None if i==len(knots) else knots[i],m,b) for i,(m,b) in enumerate(zip(slopes,bs))] def value(p,x): x=F(x) for lo,hi,m,b in p: if (lo is None or x>=lo) and (hi is None or x<=hi):return m*x+b raise AssertionError('pieces do not cover R') def hopf(p,t,x): t=F(t);x=F(x) if t<=0:raise ValueError('positive time required') candidates=[] for lo,hi,m,b in p: y=x-t*m if lo is not None:y=max(y,lo) if hi is not None:y=min(y,hi) candidates.append((m*y+b+(x-y)**2/(2*t),y)) best=min(v for v,y in candidates) return best,sorted({y for v,y in candidates if v==best}) def plus(t,x):return x*x/(2*t) if abs(x)<=t else abs(x)-t/2 def minus(t,x):return -abs(x)-t/2 def run(): counts=Counter() def ck(k,b): counts[k]+=1 if not b:raise AssertionError((k,counts[k])) # Transport: u=sum coefficient*x^j*E^k, E=e^-t, symbolic coefficient residual. terms={(2,3):F(1),(1,0):F(1,2),(1,2):F(-1,2)} residual={key:(-key[1]+key[0]+1)*v for key,v in terms.items()};residual[(1,0)]-=1 ck('transport_polynomial_identity',all(v==0 for v in residual.values())) for E,x in product([F(i,20) for i in range(1,21)],[F(i,7) for i in range(-21,22)]): u=E**3*x*x+x*(1-E**2)/2;ut=-3*E**3*x*x+x*E**2;ux=2*E**3*x+(1-E**2)/2 ck('transport_residual',ut+x*ux+u==x) xi=E*x;along=E*xi*xi+xi*(1/E-E)/2;ck('transport_inverse_label',u==along) ck('transport_capstone',(F(1,2)**3*4+1*(1-F(1,2)**2))==F(5,4)) # Rational unit-circle parameters encode exact Hamilton phases; no trig roundoff. for q,xi in product([F(i,30) for i in range(31)],[F(i,3) for i in range(-15,16)]): c=(1-q*q)/(1+q*q);s=2*q/(1+q*q);X=xi*c;P=-xi*s;U=-xi*xi*s*c/2 ck('hamilton_energy',X*X+P*P==xi*xi) ck('hamilton_action_derivative',-xi*xi*(c*c-s*s)/2==P*P-(X*X+P*P)/2) ck('hamilton_label_gradient',-xi*s*c==P*c) if c: ux=-X*s/c;ut=-X*X/(2*c*c);ck('hamilton_PDE',ut+(X*X+ux*ux)/2==0 and ux==P and -X*X*s/(2*c)==U) else:ck('projection_collapse',X==0 and P==-xi) pplus=piecewise([0],[-1,1]);pminus=piecewise([0],[1,-1]);models=[pplus,pminus,piecewise([-3,-1,0,2],[F(1,2),-1,2,-2,F(1,3)]),piecewise([-2,1],[-2,0,1]),piecewise([0],[0,0],F(3,7))] times=[F(1,10),F(1,3),F(1,2),F(1),F(3,2),F(3)];xs=[F(i,8) for i in range(-32,33)] for p in models: L=max(abs(m) for lo,hi,m,b in p) for t,x in product(times,xs): u,ys=hopf(p,t,x);ck('all_minimizer_displacement',all(abs(x-y)<=L*t for y in ys));ck('initial_uniform_bound',value(p,x)-L*L*t/2<=u<=value(p,x)) for y in ys:ck('minimizer_value',value(p,y)+(x-y)**2/(2*t)==u) # Independent evaluation on nearby competitors supplements exact piecewise minima. for z in [x+F(j,3) for j in range(-4,5)]:ck('competitor_lower_bound',u<=value(p,z)+(x-z)**2/(2*t)) h=F(2,7);v,_=hopf(p,t,x+h);ck('spatial_Lipschitz',abs(u-v)<=L*abs(h)) w,_=hopf(p,t+F(1,5),x);ck('temporal_Lipschitz',abs(w-u)<=L*L/10) for s,t,x in product(times[:4],times[:4],xs[::4]): whole,ys=hopf(p,s+t,x) for y in ys: z=(t*y+s*x)/(s+t);ck('two_leg_cost',((z-y)**2/(2*s)+(x-z)**2/(2*t))==(x-y)**2/(2*(s+t))) previous,_=hopf(p,s,z);ck('exact_semigroup',previous+(x-z)**2/(2*t)==whole) for t,x in product(times,xs[::4]): shifted=[(lo,hi,m,b+F(3,11)) for lo,hi,m,b in p];ck('constant_shift_contraction',hopf(shifted,t,x)[0]-hopf(p,t,x)[0]==F(3,11)) for t,x in product(times,xs): up,yp=hopf(pplus,t,x);um,ym=hopf(pminus,t,x) expectedplus=max(x-t,F(0)) if x>=0 else min(x+t,F(0));expectedminus=[x+t] if x>0 else ([x-t] if x<0 else [-t,t]) ck('absolute_value_envelopes',up==plus(t,x) and um==minus(t,x)) ck('all_absolute_value_labels',yp==[expectedplus] and ym==expectedminus) if abs(x)=0 else F(-1) ck('positive_envelope_PDE',ut+ux*ux/2==0) for p in [F(i,30) for i in range(-30,31)]: ck('downward_cusp_subtest',F(-1,2)+p*p/2<=0) ck('eikonal_upper_test',abs(p)-1<=0) ck('upward_HJ_cusp_rejected',F(-1,2)+F(0)**2/2<0);ck('negative_tent_rejected',abs(F(0))-1<0) for t in times: ck('positive_join_derivatives',(-t*t/(2*t*t),t/t)==(F(-1,2),F(1))) # Domain may be disconnected; nearest point and short path remain componentwise valid. intervals=[(F(-3),F(-1)),(F(0),F(2)),(F(4),F(7))];boundary=[z for pair in intervals for z in pair] def distance(x):return min(abs(x-a) for a in boundary) for lo,hi in intervals: for i in range(1,30): x=lo+(hi-lo)*F(i,30);d=distance(x);near=[a for a in boundary if abs(x-a)==d] for a in near: e=(x-a)/d for h in [d/4,d/2,3*d/4]:ck('nearest_boundary_direction',distance(x-h*e)==d-h and lo