import math,json,csv,statistics,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

def water(Q,c,r,e):
 if not all(math.isfinite(v)for v in [Q,c,r,e])or Q<0 or e<0 or not 0<=c<=1 or not 0<=r<=1:raise ValueError('Invalid ledger input')
 recovered=Q*c*r;unc=Q*(1-c);reject=Q*c*(1-r)
 return dict(throughput_kg_day=Q,collection_fraction=c,recovery_fraction=r,energy_kWh_per_processed_kg=e,recovered_kg_day=recovered,uncollected_kg_day=unc,reject_kg_day=reject,makeup_kg_day=unc+reject,energy_kWh_day=e*Q*c)
def constituent(F,w,y):
 if not all(math.isfinite(v)for v in [F,w,y])or F<0 or not 0<=w<=1 or not 0<=y<=1:raise ValueError('Invalid constituent input')
 return dict(feed_kg=F,feed_constituent_fraction=w,retention_fraction=y,incoming_constituent_kg=F*w,retained_constituent_kg=F*w*y,residual_constituent_kg=F*w*(1-y))
def reading_only(values):
 if len(values)!=3 or any(v is not None and (not math.isfinite(v) or not 0<=v<=2) for v in values):raise ValueError('Three normalized readings required')
 valid=[v for v in values if v is not None]
 median=statistics.median(valid)if len(valid)>=2 else None
 return median,None if median is None else median>=1,len(valid)<3
wr=[water(100,i/20,j/20,.02)for i in range(21)for j in range(21)]
cr=[constituent(100,i/20,j/20)for i in range(21)for j in range(21)]
cases=[('All below threshold',.5,[.5,.5,.5]),('All above threshold',1.5,[1.5,1.5,1.5]),('One high dissent',.5,[.5,.5,1.5]),('One low dissent',1.5,[1.5,1.5,.5]),('Shared low readings miss',1.5,[.5,.5,1.5]),('Shared high readings alarm',.5,[1.5,1.5,.5]),('One missing mixed readings',1.5,[None,.5,1.5]),('One missing both low',1.5,[None,.5,.5]),('Two missing one high',1.5,[None,None,1.5]),('All missing',1.5,[None,None,None]),('Nominal stuck indistinguishable',.5,[.5,.5,.5])]
sr=[]
for label,truth,values in cases:
 median,alarm,degraded=reading_only(values);sr.append(dict(case=label,truth_for_scoring=truth,reading1=values[0],reading2=values[1],reading3=values[2],median=median,threshold=1,alarm=alarm,incomplete=degraded,scored_correct=None if alarm is None else alarm==(truth>=1)))
assert sr[4]['alarm']is False and sr[4]['scored_correct']is False
assert reading_only([None,.5,1.5])==(1.,True,True)
assert reading_only([None,None,None])==(None,None,True)
assert reading_only([.5,.5,1.5])==reading_only([1.5,.5,.5])
for r in wr:assert abs(r['recovered_kg_day']+r['uncollected_kg_day']+r['reject_kg_day']-100)<1e-12
for r in cr:assert abs(r['retained_constituent_kg']+r['residual_constituent_kg']-r['incoming_constituent_kg'])<1e-12
assert water(0,1,1,0)['energy_kWh_day']==0 and water(100,1,1,.02)['makeup_kg_day']==0 and water(100,0,1,.02)['makeup_kg_day']==100
assert constituent(0,1,1)['retained_constituent_kg']==0
bad=0
for fn,args in [(water,(-1,1,1,1)),(water,(100,1.1,1,1)),(water,(100,1,-.1,1)),(water,(100,1,1,-1)),(constituent,(1,1.1,.5)),(constituent,(1,.5,-1))]:
 try:fn(*args)
 except ValueError:bad+=1
 else:raise AssertionError(args)
def save(n,r):
 with(P/n).open('w',newline='',encoding='utf8')as f:w=csv.DictWriter(f,fieldnames=list(r[0]));w.writeheader();w.writerows(r)
for n,r in [('water-ledger.csv',wr),('constituent-ledger.csv',cr),('sensor-cases.csv',sr)]:save(n,r)
plt.rcParams.update({'font.size':10,'axes.spines.top':False,'axes.spines.right':False,'savefig.dpi':180})
fig,ax=plt.subplots(1,2,figsize=(11,4.5),layout='constrained')
for c,dash in [(0.5,':'),(.8,'--'),(1,'-')]:
 r=[r for r in wr if r['collection_fraction']==c];ax[0].plot([a['recovery_fraction']for a in r],[a['makeup_kg_day']for a in r],ls=dash,label=f'Collection {c:.0%}')
ax[0].set(xlabel='Recovery from collected stream (fraction)',ylabel='Required makeup (kg/day)',title='A  Collection and treatment are separate');ax[0].legend()
c=np.linspace(0,1,101);ax[1].plot(c,.02*100*c);ax[1].set(xlabel='Collected fraction',ylabel='Assigned treatment energy (kWh/day)',title='B  0.02 kWh per processed kg');fig.suptitle('Hypothetical steady throughput 100 kg/day; no water-quality model')
for e in ['png','svg','pdf']:fig.savefig(P/f'water-loop.{e}')
plt.close(fig)
fig,ax=plt.subplots(figsize=(8,4.8),layout='constrained');y=np.linspace(0,1,101);ax.plot(y,30*y,label='Retained constituent',ls='-');ax.plot(y,30*(1-y),label='Residual constituent',ls='--');ax.axhline(30,color='gray',ls=':',label='Incoming constituent');ax.set(xlabel='Retention fraction',ylabel='Tracked constituent mass (kg)',title='100 kg feed with 30% tracked constituent\nConservation alone establishes neither edible product nor fuel');ax.legend();ax.grid(alpha=.2)
for e in ['png','svg','pdf']:fig.savefig(P/f'constituent-balance.{e}')
plt.close(fig)
(P/'validation.json').write_text(json.dumps(dict(all_pass=True,water_rows=441,constituent_rows=441,sensor_cases=11,conservation=True,endpoints=True,readings_only_shared_miss_preserved=True,missing_median_policy=True,invalid_guards=bad),indent=2)+'\n')
(P/'study.json').write_text(json.dumps(dict(runtime=dict(python=sys.version,numpy=np.__version__,matplotlib=matplotlib.__version__),generator_sha256=hashlib.sha256(Path(__file__).read_bytes()).hexdigest(),synthetic=True,water=dict(Q_kg_day=100,energy_kWh_per_processed_kg=.02,collection_and_recovery_step=.05),constituent=dict(feed_kg=100,composition_and_retention_step=.05),sensors=dict(threshold=1,rule='Median of valid finite readings; >=1 alarm; fewer than2 readings =>unknown; incomplete if fewer than3',truth='Scoring only; never passed to algorithm',limits='No fault identity inferred, no safe physical action, no probability estimates')),indent=2)+'\n');print('PASS893rows')