import json, math, random, subprocess
from fractions import Fraction
from pathlib import Path
import scipy
from scipy.special import gammaincc
from scipy.stats import poisson
root = Path(__file__).resolve().parent
rng = random.Random(638141)
cases=[]
# Probe transition regions throughout the UI domain, including its boundaries.
for h in [1, 5, 10]:
    for mu in [.005, .016, .05]:
        for cv in [.05, .1, .5, 1]:
            r = h/mu
            for z in [-12,-6,-3,-1,0,1,3,6,12]:
                n = max(1, round(r + z*cv*math.sqrt(r)))
                cases.append(('cdf',n,h,mu,cv))
for _ in range(10000):
    h = rng.uniform(1,10)
    mu = rng.uniform(.005,.05)
    cv = rng.uniform(.05,1)
    r = h/mu
    n = max(1,round(r+rng.uniform(-8,8)*cv*math.sqrt(r)))
    cases.append(('cdf',n,h,mu,cv))
# Exhaustive integer-count monotonicity along selected scenarios.
monotonic_spans=[]
for h,mu,cv in [(1,.05,.05),(1,.05,1),(5,.016,.5),(10,.005,.05),(10,.005,1)]:
    start=len(cases)
    limit=math.ceil(h/mu+12*cv*math.sqrt(h/mu)+50)
    cases.extend(('cdf',n,h,mu,cv) for n in range(limit+1))
    monotonic_spans.append((start,len(cases),[h,mu,cv]))
quantile_specs=[(p,h,mu,cv) for h,mu,cv in [(5,.016,.5),(1,.05,.05),(1,.05,1),(10,.005,.05),(10,.005,1)] for p in [.1,.5,.9]]
quantile_start=len(cases)
cases.extend(('quantile',*args) for args in quantile_specs)
outputs=json.loads(subprocess.check_output(['node',str(root/'validate-model.cjs')],input=json.dumps(cases).encode()))
maxerr=(0,None)
for args,out in zip(cases[:quantile_start],outputs[:quantile_start]):
    _,n,h,mu,cv=args
    expected=0 if n<1 else float(gammaincc(n/cv**2,(h/mu)/cv**2))
    error=abs(expected-out)
    if error>maxerr[0]: maxerr=(error,{'args':args,'javascript':out,'scipy':expected})
    assert math.isfinite(out) and 0<=out<=1
    assert error<1e-7, (args,out,expected)
monotonic=[]
for start,end,scenario in monotonic_spans:
    sequence=outputs[start:end]
    drops=[(i,sequence[i-1],sequence[i]) for i in range(1,len(sequence)) if sequence[i]<sequence[i-1]]
    assert not drops, (scenario,drops)
    monotonic.append({'h_mu_cv':scenario,'counts_tested':len(sequence),'violations':len(drops)})
quantiles=[]
for args,actual in zip(quantile_specs,outputs[quantile_start:]):
    p,h,mu,cv=args
    lo=0; hi=max(1,math.ceil(h/mu))
    def f(n): return 0 if n<1 else float(gammaincc(n/cv**2,(h/mu)/cv**2))
    while f(hi)<p: hi*=2
    while hi-lo>1:
        mid=(hi+lo)//2
        if f(mid)>=p: hi=mid
        else: lo=mid
    assert actual==hi, (args,actual,hi)
    quantiles.append({'p_h_mu_cv':args,'count':actual,'cdf_previous':f(actual-1),'cdf_count':f(actual)})
poisson_cases=[('cdf',n,1,1/r,1) for r in [20,50,312.5,2000] for n in range(max(1,int(r-5*math.sqrt(r))),int(r+5*math.sqrt(r))+1)]
poisson_outputs=json.loads(subprocess.check_output(['node',str(root/'validate-model.cjs')],input=json.dumps(poisson_cases).encode()))
poisson_maxerr=max(abs(out-poisson.cdf(n-1,h/mu)) for (_,n,h,mu,_),out in zip(poisson_cases,poisson_outputs))
assert poisson_maxerr<1e-7
deterministic_cases=[]
deterministic_expected=[]
scenarios=0
for h_tenths in range(10,101):
    for mu_thousandths in range(5,51):
        scenarios+=1
        h=float(Fraction(h_tenths,10))
        mu=float(Fraction(mu_thousandths,1000))
        expected=math.ceil(Fraction(h_tenths*100,mu_thousandths))
        deterministic_cases.extend([('equalCount',h,mu),('cdf',expected-1,h,mu,0),('cdf',expected,h,mu,0),('quantile',.1,h,mu,0),('quantile',.5,h,mu,0),('quantile',.9,h,mu,0)])
        deterministic_expected.extend([expected,0,1,expected,expected,expected])
deterministic_outputs=json.loads(subprocess.check_output(['node',str(root/'validate-model.cjs')],input=json.dumps(deterministic_cases).encode()))
assert deterministic_outputs==deterministic_expected
result={'reference':'SciPy '+scipy.__version__,'cdf_checks':quantile_start,'maximum_absolute_cdf_error':maxerr[0],'worst_case':maxerr[1],'monotonicity':monotonic,'quantiles':quantiles,'poisson_identity_checks':len(poisson_cases),'maximum_absolute_poisson_error':poisson_maxerr,'deterministic_exact_rational_checks':len(deterministic_cases),'deterministic_scenarios':scenarios}
(root/'model-validation.json').write_text(json.dumps(result,indent=2)+'\n')
print(json.dumps(result,indent=2))