import math,json,csv,itertools,sys,hashlib
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
P=Path(__file__).resolve().parent
gamma=.5772156649015328606

def e1(a):
    if not math.isfinite(a) or not 0<a<=1:raise ValueError('Series domain 0<a<=1')
    term=1.;total=0.
    for k in range(1,101):
        term*=(-a)/k;total+=term/k
        if abs(term/k)<1e-18:break
    return -gamma-math.log(a)-total

def quadrature(a,n):
    upper=math.log(50/a);step=upper/n
    vals=[math.exp(-a*math.exp(i*step)) for i in range(n+1)]
    return step/3*(vals[0]+vals[-1]+4*math.fsum(vals[1:-1:2])+2*math.fsum(vals[2:-1:2]))

def name(t,f):return {(0,0):'neither',(1,0):'for_only',(0,1):'against_only',(1,1):'both'}[t,f]
states=[(0,0),(1,0),(0,1),(1,1)]
unary=[dict(state=name(*a),support_for=a[0],support_against=a[1],negation=name(a[1],a[0])) for a in states]
binary=[]
for a,b in itertools.product(states,repeat=2):
    binary.append(dict(left=name(*a),right=name(*b),conjunction=name(a[0]&b[0],a[1]|b[1]),disjunction=name(a[0]|b[0],a[1]&b[1]),merge=name(a[0]|b[0],a[1]|b[1])))
arc=[dict(angle_degrees=k,angle_radians=k*math.pi/180,radius=1,arc_length=k*math.pi/180,chord_length=2*math.sin(k*math.pi/360)) for k in range(181)]
geo=[];convergence=[]
for k in range(121):
    r=10**(k/40);a=1/r**2;v=math.pi*e1(a);cap=math.pi*r*r*(-math.expm1(-a))
    geo.append(dict(top_radius=r,top_height=math.exp(-a),top_slice_volume=v,cap_volume=cap,cylindrical_clip_volume=v+cap))
for r in [1,2,10,100,1000]:
    for n in [250,500,1000,2000]:convergence.append(dict(radius=r,steps=n,quadrature_volume=math.pi*quadrature(1/r**2,n),series_volume=math.pi*e1(1/r**2),absolute_error=math.pi*abs(quadrature(1/r**2,n)-e1(1/r**2))))
def save(name,rows):
    with(P/name).open('w',newline='',encoding='utf8')as f:
        w=csv.DictWriter(f,fieldnames=list(rows[0]));w.writeheader();w.writerows(rows)
for n,r in [('unary-reports.csv',unary),('binary-reports.csv',binary),('arc-chord.csv',arc),('finite-geometry.csv',geo),('quadrature-check.csv',convergence)]:save(n,r)
assert all(r['chord_length']<=r['arc_length']+1e-15 for r in arc)
assert arc[0]['chord_length']==arc[0]['arc_length']==0 and arc[-1]['chord_length']==2
assert all(b['top_slice_volume']>a['top_slice_volume'] for a,b in zip(geo,geo[1:]))
for a in states:
    assert name(a[0],a[1])==name(*a)
    for b in states:
        assert (a[0]|b[0],a[1]|b[1])==(b[0]|a[0],b[1]|a[1])
        # Truth tables independently checked against sets of supporting assertions.
        row=next(r for r in binary if r['left']==name(*a) and r['right']==name(*b))
        assert row['conjunction']==name(int(all([a[0],b[0]])),int(any([a[1],b[1]])))
        assert row['disjunction']==name(int(any([a[0],b[0]])),int(all([a[1],b[1]])))
        assert row['merge']==name(int(any([a[0],b[0]])),int(any([a[1],b[1]])))
assert max(r['absolute_error']for r in convergence if r['steps']==2000)<1e-9
legacy_paths=list((P.parent.parent/'prototypes').glob('*axis_geometry_results.json'))
parity=[];hashes={}
if legacy_paths:
    path=legacy_paths[0];old=json.loads(path.read_text());hashes['reference_results']=hashlib.sha256(path.read_bytes()).hexdigest();src=path.with_name(path.name.replace('_results.json','.py'));hashes['reference_generator']=hashlib.sha256(src.read_bytes()).hexdigest()
    for r in old['results_if_revolved_about_y']['truncated_volumes']:
        val=math.pi*e1(1/r['radius_cap']**2);parity.append(dict(radius=r['radius_cap'],series=val,reference=r['volume'],difference=val-r['volume']))
    assert max(abs(r['difference'])for r in parity)<1e-10
if parity:save('reference-comparison.csv',parity)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(1,2,figsize=(10,4.4),layout='constrained');theta=np.linspace(0,2*math.pi/3,201)
ax[0].plot(np.cos(theta),np.sin(theta),label='Allowed arc',color='#7957a8');ax[0].plot([1,-.5],[0,math.sqrt(3)/2],label='Interior chord',color='#278c89');ax[0].set_aspect('equal');ax[0].set(xlabel='Horizontal coordinate',ylabel='Vertical coordinate',title='A  Two paths between the same endpoints');ax[0].legend(fontsize=8)
ax[1].plot([r['angle_degrees']for r in arc],[r['arc_length']for r in arc],label='Arc');ax[1].plot([r['angle_degrees']for r in arc],[r['chord_length']for r in arc],label='Chord');ax[1].set(xlabel='Angle (degrees)',ylabel='Length / circle radius',title='B  Geometry does not establish physical access');ax[1].legend();fig.suptitle('Unit circle: path length depends on the permitted route')
for ext in ['png','svg','pdf']:fig.savefig(P/f'arc-and-chord.{ext}')
plt.close(fig)
fig,ax=plt.subplots(1,2,figsize=(11,4.5),layout='constrained');ys=np.linspace(.001,math.exp(-.25),300);rs=1/np.sqrt(-np.log(ys))
ax[0].fill_betweenx(ys,-rs,rs,alpha=.35,label='Top-sliced solid section');ax[0].fill_betweenx([math.exp(-.25),1],-2,2,alpha=.3,color='#7957a8',label='Extra cylindrical-clip cap');ax[0].axhline(1,ls='--',color='gray');ax[0].set(xlabel='Horizontal coordinate (dimensionless)',ylabel='Height (dimensionless)',title='A  R = 2: two different truncations');ax[0].legend(fontsize=8)
for field,label in [('top_slice_volume','Stop at h = exp(−1/R²)'),('cylindrical_clip_volume','Clip full region at cylindrical radius R')]:ax[1].semilogx([r['top_radius']for r in geo],[r[field]for r in geo],label=label)
ax[1].set(xlabel='Cutoff radius R (dimensionless)',ylabel='Volume (dimensionless)',title='B  Bounded height, growing volume');ax[1].legend(fontsize=8);ax[1].grid(alpha=.2);fig.suptitle('Geometric idealization: no mass, energy or propulsion model')
for ext in ['png','svg','pdf']:fig.savefig(P/f'cutoff-geometry.{ext}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,quadrature_2000_max_error=max(r['absolute_error']for r in convergence if r['steps']==2000),reference_max_error=max((abs(r['difference'])for r in parity),default=None),arc_limits=True,monotone_volume=True,rows=dict(unary=4,binary=16,arc=181,geometry=121,quadrature=20,reference=4)),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),baseline_sha256=hashes,dimensionless=True,logic='Recorded support-for/support-against assertions, not probabilities',radius_range=[1,1000],radius_points=121,series_domain='0<a<=1 only',quadrature='Composite Simpson after t=a exp(s), upper t=50; neglected tail below exp(-50)/50',geometry='Top-slice pi E1(1/R²); cylindrical clip adds pi R²(1-exp(-1/R²))'),indent=2)+'\n')
print('PASS',json.loads((P/'validation.json').read_text()))