"""Exact finite and analytic checks for articles 11--15.
Run from any directory with Python 3 + NumPy. No network, training, or files written.
Curves are declared mathematical constructions, not observed model benchmarks.
"""
import itertools
import math
import numpy as np

TOPICS={
11:"diffusion-score-error",12:"flow-matching-transport",
13:"variational-posterior-collapse",14:"low-rank-adaptation-capacity",
15:"sparse-expert-routing"}

def run():
    panels=[]; results=[]; points=[]; count=0
    def ck(ok,message):
        nonlocal count
        count+=1
        if not bool(ok): raise AssertionError(message)
    def close(a,b,atol=1e-10):
        return bool(np.allclose(a,b,rtol=1e-10,atol=atol))
    def series(label,x,y):
        return {"label":label,"x":[float(a) for a in x],"y":[float(a) for a in y]}
    def panel(n,k,title,xlabel,ylabel,seq,caption,xscale="linear",yscale="linear"):
        slug=TOPICS[n]
        panels.append({"article":n,"number":k,"slug":slug,"title":title,"xlabel":xlabel,"ylabel":ylabel,
            "series":seq,"caption":caption,"xscale":xscale,"yscale":yscale,
            "kind":"constructed_analytic_example","image":f"figures/{slug}-{k}.png",
            "vector":f"figures/{slug}-{k}.svg","data":f"figures/{slug}-{k}.csv",
            "alt":title+". "+caption})
    def result(n,start,label,witness,scope):
        results.append({"article":n,"slug":TOPICS[n],"check":label,"assertions":count-start,
            "witness":witness,"scope":scope,"status":"PASS"})

    # 11: additive-noise conversion and reverse-clock stationary OU moments.
    start=count
    sigma=np.logspace(-2,1,61);delta=.1;s_ref=1. # one fixed data unit
    panel(11,1,"One noise perturbation, two normalized errors","Noise ratio sigma / s_ref","Dimensionless normalized error",
        [series("Normalized score error",sigma/s_ref,s_ref*delta/sigma),series("Normalized clean-estimate error",sigma/s_ref,delta*sigma/s_ref)],
        "Constructed additive-corruption calculation: dimensionless noise-prediction error delta=0.1 and fixed s_ref=one data unit. Plot s_ref*abs(score error)=s_ref*delta/sigma and abs(clean-estimate error)/s_ref=delta*sigma/s_ref. Both axes are dimensionless; these are distinct normalized errors, not interchangeable quality metrics or image-model measurements.",
        "log","log")
    points.append([11,1,0,1.,.1])
    tau=np.linspace(0,4,81)
    panel(11,2,"A converged process can retain score bias","Reverse time tau","Mean displacement",
        [series(f"Score bias b={b:g}",tau,2*b*(-np.expm1(-tau/2))) for b in [0,.1,.25]],
        "Exact Gaussian reverse-clock process dX=(-X/2+b)dt+dW, initialized N(0,1). Mean=2b(1-exp(-t/2)); continuous variance remains one. Curves are analytical, not sampler measurements.")
    points.append([11,2,2,2.,.5*(1-math.exp(-1))])
    for s in [.1,.5,1.,2.]:
        var=1+s*s
        for y in [-2.,-.3,0.,.7,2.]:
            h=1e-5
            logp=lambda z: -.5*math.log(2*math.pi*var)-z*z/(2*var)
            numerical=(logp(y+h)-logp(y-h))/(2*h)
            oracle=-y/var
            ck(close(numerical,oracle,1e-8),"Gaussian score derivative")
            noise=s*y/var
            ck(close(-noise/s,oracle),"conditional mean noise conversion")
    errors=[]
    for b in [0.,.1,.25,-.2]:
        errs=[]
        for N in [4,8,16,32,64,128]:
            h=2/N;r=1-h/2;m=0.;v=1.
            for _ in range(N):m=r*m+h*b;v=r*r*v+h
            formula=2*b*(1-r**N)
            vv=r**(2*N)+h*(1-r**(2*N))/(1-r*r)
            ck(close(m,formula),"OU Euler mean recurrence")
            ck(close(v,vv),"OU Euler variance recurrence")
            ck(v>1,"Euler finite-step variance is not exactly stationary")
            errs.append(abs(m-2*b*(1-math.exp(-1))))
        ck(all(errs[i+1]<=errs[i]+1e-14 for i in range(len(errs)-1)),"OU refinement converges to biased mean")
        errors.append({"b":b,"mean_error_at_N128":errs[-1]})
    ck(.5*(1-math.exp(-1))>.3,"nonzero limiting mean with bias .25")
    result(11,start,"Noise-to-score identity and exact biased OU mean/variance recurrences",
        {"T":2,"step_counts":[4,8,16,32,64,128],"biases":[0,.1,.25,-.2],
         "limiting_mean_b025":.5*(1-math.exp(-1)),"refinement":errors},
        "Scalar Gaussian fields and listed Euler discretizations only; no learned predictor or image sampler is tested.")

    # 12: covariance algebra, regular Gaussian marginal field, and Euler/Heun.
    start=count
    t=np.linspace(0,1,101);v=(1-t)**2+t*t
    panel(12,1,"Same Gaussian endpoints, different interior variances","Interpolation time t","Marginal variance",
        [series("Independent coupling",t,v),series("Identity coupling",t,np.ones_like(t)),
         series("Anticorrelated coupling",t,(1-2*t)**2)],
        "Constructed standard-normal endpoint couplings. The anticorrelated case collapses at t=0.5 and is not a regular invertible ODE flow through that time. Independent coupling has midpoint variance 0.5.")
    points.append([12,1,0,.5,.5])
    def a(t):return (2*t-1)/((1-t)**2+t*t)
    def solve(N,method):
        h=1/N;x=1.
        for i in range(N):
            u=i*h;k1=a(u)*x
            if method=="Euler":x+=h*k1
            else:x+=h*(k1+a(u+h)*(x+h*k1))/2
        return x
    Ns=list(range(1,65));ey=[abs(solve(n,"Euler")-1) for n in Ns];hy=[abs(solve(n,"Heun")-1) for n in Ns]
    panel(12,2,"Exact velocity does not remove integration error","Uniform step count N","Absolute endpoint scale error",
        [series("Euler (N evaluations)",Ns,ey),series("Heun (2N evaluations)",Ns,hy)],
        "Independent-Gaussian marginal ODE x'=a(t)x, x(0)=1. Exact x(1)=1. Error is computed from deterministic solver trajectories; no network is trained. Heun uses twice as many field evaluations per step.","log","linear")
    points.append([12,2,0,1.,1.])
    for tt in np.linspace(.01,.99,51):
        weights=np.array([1-tt,tt]);disp=np.array([-1.,1.])
        cov=float(weights@disp);variance=float(weights@weights)
        ck(close(cov/variance,a(tt)),"Gaussian conditional velocity coefficient")
        h=1e-6
        exact=lambda u:math.sqrt((1-u)**2+u*u)
        derivative=(exact(tt+h)-exact(tt-h))/(2*h)
        ck(close(derivative,a(tt)*exact(tt),1e-8),"analytic trajectory satisfies ODE")
    ck(close(solve(1,"Euler"),0),"one-step Euler collapse")
    ck(close(solve(1,"Heun"),.5),"one-step Heun scale")
    ck(ey[-1]<ey[0] and hy[-1]<hy[0],"endpoint refinement improvement")
    ck(close(v[50],.5) and close((1-2*t[50])**2,0),"midpoint variance witnesses")
    result(12,start,"Conditional covariance and independently known Gaussian flow endpoints",
        {"independent_midpoint_variance":.5,"anticorrelated_midpoint_variance":0,
         "Euler_N1_scale":solve(1,"Euler"),"Heun_N1_scale":solve(1,"Heun"),
         "Euler_N64_error":ey[-1],"Heun_N64_error":hy[-1]},
        "Gaussian coupling identities and two explicitly coded scalar solvers; the singular anticorrelated path is only a counterexample, not a valid regular flow.")

    # 13: finite information decomposition and scalar constrained Gaussian MLE.
    start=count
    def entropy_term(p):return 0. if p==0 else p*math.log(p)
    def binary_info(e):return math.log(2)+entropy_term(e)+entropy_term(1-e)
    e=np.linspace(0,.5,101)
    biased=.9*math.log(1.8)+.1*math.log(.2)
    panel(13,1,"Positive coding rate need not encode the input","Binary-channel flip probability epsilon","Information or rate (nats)",
        [series("Informative channel: I = rate",e,[binary_info(float(x)) for x in e]),
         series("Input-independent biased encoder: rate",e,np.full_like(e,biased)),
         series("Input-independent biased encoder: I",e,np.zeros_like(e))],
        "Exact finite construction: uniform binary X and prior. Informative channel flips X with probability epsilon and has aggregate-prior KL zero. Separate input-independent encoder emits Z=1 with probability 0.9; its positive rate is entirely aggregate mismatch.")
    points.append([13,1,0,0.,math.log(2)])
    ss=np.linspace(.1,6,60)
    panel(13,2,"Observation noise sets a latent-loading threshold","Fixed observation-noise variance s","Optimal loading energy a squared",
        [series(f"Data variance lambda={lam}",ss,np.maximum(lam-ss,0)) for lam in [1,3,5]],
        "Constructed scalar Gaussian model z~N(0,1), x|z~N(az,s). Expected marginal likelihood for centered data variance lambda is maximized at a^2=max(lambda-s,0). This is not a trained deep VAE result.")
    points.append([13,2,1,3.,0.])
    def information(q):
        px=np.array([.5,.5]);prior=np.array([.5,.5]);joint=px[:,None]*q;qz=joint.sum(axis=0)
        rate=0.;mi=0.;agg=0.
        for x in range(2):
            for z in range(2):
                if q[x,z]>0:
                    rate+=joint[x,z]*math.log(q[x,z]/prior[z])
                    mi+=joint[x,z]*math.log(q[x,z]/qz[z])
        for z in range(2):
            if qz[z]>0:agg+=qz[z]*math.log(qz[z]/prior[z])
        return rate,mi,agg
    for eps in e:
        rate,mi,agg=information(np.array([[1-eps,eps],[eps,1-eps]]))
        ck(close(rate,mi+agg),"finite rate decomposition")
        ck(close(mi,binary_info(float(eps))),"finite mutual information versus entropy")
        ck(close(agg,0),"aggregate prior matched")
    br,bi,ba=information(np.array([[.1,.9],[.1,.9]]))
    ck(close(bi,0) and br>0 and close(br,ba),"positive rate without input information")
    for lam in [.5,1.,3.,5.]:
        for s in [.1,.5,1.,3.,6.]:
            best=max(lam-s,0)
            f=lambda energy:.5*(math.log(energy+s)+lam/(energy+s))
            # Every candidate on this fixed grid must be no better than the analytic optimum.
            candidates=np.linspace(0,8,161)
            ck(all(f(best)<=f(float(c))+1e-12 for c in candidates),"constrained Gaussian optimum against independent grid")
            if best>0:ck(close(1/(best+s)-lam/(best+s)**2,0),"interior likelihood derivative")
            else:ck((1/s-lam/s**2)>=-1e-12,"boundary likelihood derivative")
    result(13,start,"Finite binary KL decomposition and constrained Gaussian loading optimum",
        {"binary_channel_grid":101,"independent_encoder_rate_nats":br,"independent_encoder_information_nats":bi,
         "lambda3_noise3_loading_energy":0,"Gaussian_parameter_pairs":20,"candidate_energies_per_pair":161},
        "Exact finite information sums and scalar Gaussian expected likelihood. No deep encoder, decoder, or mutual-information estimator is trained.")

    # 14: SVD residuals, gauge equivalence, and covariance-weighted error.
    start=count
    ranks=np.arange(9);flat=np.ones(8);decay=.5**np.arange(8)
    tail=lambda spectrum: [float(np.sum(spectrum[r:]**2)/np.sum(spectrum**2)) for r in ranks]
    panel(14,1,"Rank efficiency depends on the declared spectrum","Retained rank r","Fraction of squared spectral energy discarded",
        [series("Flat singular values",ranks,tail(flat)),series("Geometric singular values",ranks,tail(decay))],
        "Constructed 8x8 target matrices. Flat spectrum is all ones; geometric spectrum is (1,1/2,...,1/128). Curves are optimal Frobenius residual energy fractions, not measured adaptation quality.")
    points.append([14,1,0,4.,.5])
    rho=np.arange(1,101,dtype=float)
    panel(14,2,"Weight-space accuracy can miss important inputs","Input variance ratio rho","Expected squared output error",
        [series("Unweighted rank-one choice",rho,rho),series("Input-weighted rank-one optimum",rho,np.minimum(16,rho))],
        "Constructed Delta W=diag(4,1), centered input covariance diag(1,rho). Unweighted truncation keeps the first direction and incurs error rho; the output-optimal rank-one residual is min(16,rho).","log","linear")
    points.append([14,2,1,100.,16.])
    rng=np.random.default_rng(1401)
    U=np.linalg.qr(rng.normal(size=(8,8)))[0];V=np.linalg.qr(rng.normal(size=(8,8)))[0]
    for spectrum in [flat,decay]:
        M=U@np.diag(spectrum)@V.T;u,s,vh=np.linalg.svd(M)
        for r in ranks:
            mr=(u[:,:r]*s[:r])@vh[:r,:]
            ck(close(np.linalg.norm(M-mr,"fro")**2,np.sum(spectrum[r:]**2)),"SVD residual energy")
            ck(np.linalg.matrix_rank(mr,tol=1e-10)<=r,"rank truncation bound")
    B=np.array([[1.,2.],[3.,4.],[0.,1.]])
    A=np.array([[1.,-1.,2.],[2.,0.,1.]])
    Q=np.array([[2.,1.],[0.,3.]])
    ck(close((B@Q)@(np.linalg.inv(Q)@A),B@A),"factor gauge identity")
    D=np.diag([4.,1.]);weightblind=np.diag([4.,0.])
    for rr in [1.,4.,16.,25.,100.]:
        S=np.diag([1.,rr]);L=np.diag([1.,math.sqrt(rr)])
        u,s,vh=np.linalg.svd(D@L);M=((u[:,:1]*s[:1])@vh[:1,:])@np.linalg.inv(L)
        residual=D-M
        error=float(np.trace(residual@S@residual.T))
        blind=float(np.trace((D-weightblind)@S@(D-weightblind).T))
        ck(close(error,min(16,rr)),"weighted SVD optimum")
        ck(close(blind,rr),"unweighted update functional error")
    result(14,start,"Rotated SVD residuals, factor gauge invariance, and weighted linear-output error",
        {"matrix_dimension":8,"ranks":list(map(int,ranks)),"rho100_unweighted_error":100,
         "rho100_weighted_error":16,"factor_gauge_matrix":[[2,1],[0,3]]},
        "Constructed matrices and a centered linear-input criterion only. No claim about trained adapters, optimizer behavior, or accelerator savings.")

    # 15: full enumeration of all 2^8 uniform top-one routing assignments.
    start=count
    assignments=list(itertools.product([0,1],repeat=8))
    counts=np.array([[a.count(0),a.count(1)] for a in assignments])
    for pair in counts:ck(int(pair.sum())==8,"routing assignment conservation")
    capacities=np.arange(1,9)
    iid=[]
    for C in capacities:
        exact=float(np.maximum(counts-C,0).sum(axis=1).mean()/8)
        combinatorial=sum(math.comb(8,j)*(max(j-C,0)+max(8-j-C,0)) for j in range(9))/(8*256)
        ck(close(exact,combinatorial),"enumeration versus binomial overflow")
        iid.append(exact)
    panel(15,1,"Uniform routing still overflows finite capacity","Capacity C per expert (tokens)","Expected dropped-token fraction",
        [series("Independent uniform routing",capacities,iid),
         series("Perfectly balanced groups",capacities,np.maximum(4-capacities,0)/4),
         series("All tokens choose one expert",capacities,np.maximum(8-capacities,0)/8)],
        "Exact 8-token, 2-expert, top-one drop-overflow model. Independent routing averages all 256 assignments; balanced and concentrated curves are deterministic load cases. No neural model or hardware is measured.")
    points.append([15,1,0,4.,35/256])
    mixing=np.linspace(0,1,101)
    def mixed(C):
        return (1-mixing)*iid[C-1]+mixing*max(8-C,0)/8
    panel(15,2,"Global balance does not reveal group correlation","Probability rho of fully correlated group routing","Expected dropped-token fraction",
        [series("Capacity C=4",mixing,mixed(4)),series("Capacity C=6",mixing,mixed(6))],
        "Mixture of independent fair assignments and a fair single expert for the entire group. Every rho preserves expected expert shares of one half. Curves differ only in within-group dependence and stated capacity.")
    points.append([15,2,0,1.,.5])
    ck(close(iid[3],35/256),"capacity-four rational expectation")
    ck(close(counts.mean(axis=0),[4,4]),"independent uniform expert marginals")
    ck(close(np.array([[8,0],[0,8]]).mean(axis=0),[4,4]),"correlated expert marginals agree")
    for C in [4,6]:
        for rr in [0.,.25,.5,.75,1.]:
            # Build complete finite mixture probabilities and sum overflow independently.
            fullcounts=np.vstack([counts,[[8,0],[0,8]]])
            probabilities=np.array([(1-rr)/256]*256+[rr/2,rr/2])
            expected=float(probabilities@(np.maximum(fullcounts-C,0).sum(axis=1)/8))
            ck(close(expected,(1-rr)*iid[C-1]+rr*max(8-C,0)/8),"finite mixture overflow")
            ck(close(probabilities@fullcounts,[4,4]),"mixture preserves marginal balance")
    result(15,start,"Exhaustive top-one routing and matched-marginal correlated counterexample",
        {"assignments_enumerated":256,"tokens":8,"experts":2,"capacities":list(map(int,capacities)),
         "iid_drop_at_capacity4":iid[3],"correlated_drop_at_capacity4":.5,
         "expected_loads_both":[4,4]},
        "Eight-token, two-expert top-one routing with a drop policy; integer counting and expectations only, not a trained MoE or distributed-system benchmark.")
    assert len(panels)==10 and len(results)==5 and len(points)==10
    return {"panels":panels,"results":results,"points":points}

if __name__=="__main__":
    import json
    out=run()
    print(json.dumps({"status":"PASS","panels":len(out["panels"]),"results":out["results"],"points":out["points"]},indent=2))