"""Finite analytical witnesses for articles 16--20; no training or external datasets."""

import itertools

import json

import math

from pathlib import Path

from statistics import NormalDist

import numpy as np



def run():

    topics=json.loads((Path(__file__).resolve().parent.parent/'topics.json').read_text(encoding='utf-8'))

    panels=[]; results=[]; points=[]; count=0

    def ck(condition,message):

        nonlocal count

        count+=1

        if not bool(condition): raise AssertionError(message)

    def panel(n,k,title,xlabel,ylabel,x,curves,caption):

        slug=topics[n-1][0]

        p=dict(article=n,number=k,slug=slug,title=title,xlabel=xlabel,ylabel=ylabel,

               series=[dict(label=label,x=np.asarray(x,dtype=float).tolist(),y=np.asarray(y,dtype=float).tolist()) for label,y in curves],

               caption=caption,xscale='linear',yscale='linear',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)

        panels.append(p)

    def result(n,label,start,witness,scope):

        results.append(dict(article=n,slug=topics[n-1][0],check=label,assertions=count-start,witness=witness,scope=scope,status='PASS'))



    # 16: exact linear extrema; deliberately unsuccessful first-order search.

    start=count

    w=np.full(4,0.5); corners=np.array(list(itertools.product([-1.,1.],repeat=4)))

    minima=[]

    for eps in np.linspace(0,.5,51):

        scores=1+eps*corners@w

        expected=1-2*eps

        ck(abs(float(scores.min())-expected)<1e-12,'linear box minimum')

        ck(np.all(np.abs(eps*corners)<=eps+1e-14),'feasible box corners')

        minima.append(float(scores.min()))

    ck(np.all(np.diff(minima)<=1e-14),'nested threat sets')

    delta=0.

    for _ in range(20):

        delta=float(np.clip(delta+0.1*2*delta,-1,1))

        ck(delta==0.,'stationary ascent remains at zero')

    ck(delta**2==0 and (-1.)**2==1,'failed search has a feasible better point')

    normal=NormalDist()

    p1=normal.cdf(1.)

    for sigma in [.25,.5,1.]:

        ck(abs(sigma*normal.inv_cdf(p1)-sigma)<1e-12,'normal quantile reference')

    eps=np.linspace(0,.5,101)

    panel(16,1,'Coordinate budgets accumulate with dimension','Infinity-norm budget epsilon','Worst signed linear margin',eps,

          [(f'd = {d}',1-eps*math.sqrt(d)) for d in [4,16,64]],

          'Constructed unit-Euclidean-norm weights, clean signed margin 1, and no input-domain clipping. Negative margins cross the linear decision boundary; these are different finite-dimensional models, not measurements of neural networks.')

    ps=np.unique(np.r_[np.linspace(.5,.99,100),p1])

    z=np.array([normal.inv_cdf(float(p)) for p in ps])

    panel(16,2,'Smoothing radius at supplied class probabilities','Leading-class probability p_A','Binary smoothing radius',ps,

          [(f'sigma = {s:g}',s*z) for s in [.25,.5,1.]],

          'Analytical binary Gaussian-smoothing formula R = sigma Phi^{-1}(p_A), with p_B = 1-p_A. Probabilities are supplied parameters, not finite-sample estimates; changing sigma need not preserve a real classifier probability.')

    points.extend([[16,1,0,.25,.5],[16,2,0,p1,.25]])

    result(16,'Linear corner extrema, stationary-gradient counterexample, and smoothing quantile',start,

           dict(dimensions=4,corners=16,epsilon_values=51,margin_at_epsilon_point25=.5,gradient_steps=20,stationary_loss=0.,endpoint_loss=1.,p_A_at_normal_quantile_1=p1),

           'All 16 corners of a 4D linear box at 51 radii; 20 ascent iterations on delta^2; supplied Gaussian probabilities. No trained network or attack benchmark.')



    # 17: finite outcome population and exact conditional decision costs.

    start=count

    y=np.r_[1.,np.zeros(9),np.ones(9),0.]

    q=np.r_[np.full(10,.1),np.full(10,.9)]

    constant=np.full(20,.5)

    ck(abs(y.mean()-.5)<1e-14,'aggregate probability')

    ck(abs(np.mean((constant-y)**2)-.25)<1e-14,'constant Brier')

    ck(abs(np.mean((q-y)**2)-.09)<1e-14,'informative Brier')

    for sl,truth in [(slice(0,10),.1),(slice(10,20),.9)]:

        ck(abs(y[sl].mean()-truth)<1e-14,'constructed subgroup')

        ck(abs(abs(constant[sl].mean()-y[sl].mean())-.4)<1e-14,'subgroup error')

    yp=np.array([0.,1.,0.,1.]); pp=np.array([.25,.25,.75,.75])

    fine=.5*abs(yp[:2].mean()-pp[:2].mean())+.5*abs(yp[2:].mean()-pp[2:].mean())

    coarse=abs(yp.mean()-pp.mean())

    ck(abs(fine-.25)<1e-14 and coarse==0.,'coarse bins hide event probability error')

    logits=np.array([-2.,.5,3.])

    for t in [.1,.25,1.,4.,10.]:

        probs=np.exp(logits/t-np.max(logits/t));probs/=probs.sum()

        ck(np.array_equal(np.argsort(probs),np.argsort(logits)),'positive temperature ordering')

    alpha=np.linspace(0,1,101)

    brier=[]

    for a in alpha:

        pred=.5+a*(q-.5)

        measured=float(np.mean((pred-y)**2))

        ck(abs(measured-(.09+.16*(1-a)**2))<1e-12,'interpolation Brier')

        brier.append(measured)

    panel(17,1,'Calibration alone does not reward information','Interpolation alpha toward true group probability','Expected Brier loss',alpha,

          [('Interpolated forecast',brier),('Conditional irreducible loss',np.full_like(alpha,.09))],

          'Two equally frequent Bernoulli groups have probabilities 0.1 and 0.9. Forecasts are 0.5 + alpha(q-0.5). The constant forecast at alpha=0 is aggregate calibrated but uninformative; all values are assigned population quantities.')

    costs=np.linspace(0,10,201); pooled=[]; informed=[]

    for c in costs:

        threshold=c/(c+1)

        actions=q>threshold

        cost=np.where(actions,c*(1-y),y).mean()

        pool_action=.5>threshold

        pool_cost=float(np.where(pool_action,c*(1-y),y).mean())

        expected=.5*(min(.9*c,.1)+min(.1*c,.9))

        ck(abs(cost-expected)<1e-12,'conditional decision cost')

        ck(abs(pool_cost-min(.5*c,.5))<1e-12,'pooled decision cost')

        informed.append(float(cost));pooled.append(pool_cost)

    panel(17,2,'Equal calibration can support unequal decisions','False-positive cost; false-negative cost = 1','Expected decision cost',costs,

          [('Group-informed decisions',informed),('Constant-probability decisions',pooled)],

          'Exact binary Bayes costs in the same two-group population. Each action minimizes expected cost for its available information; the constant and informed forecasts are both calibrated relative to their own information.')

    points.extend([[17,1,0,1.,.09],[17,2,0,1.,.1]])

    result(17,'Calibrated-but-uninformative population, bin merging, and cost-optimal actions',start,

           dict(cases=20,group_positive_counts=[1,9],constant_Brier=.25,informed_Brier=.09,coarse_probability_ECE=0.,fine_probability_ECE=.25,cost_at_equal_error_costs={'pooled':.5,'informed':.1}),

           'Assigned finite outcomes and population expectations; binary event-probability bins, not multiclass confidence ECE. No validation-set fitting or real-user outcomes.')



    # 18: weighting identity, clipping bias, inverse conditioning, indistinguishable worlds.

    start=count

    p=np.array([.9,.1]);qdist=np.array([.5,.5]);loss=np.array([0.,1.])

    weights=qdist/p

    weighted=float(np.sum(p*weights*loss)); target=float(qdist@loss)

    ck(abs(weighted-target)<1e-14 and target==.5,'change of measure')

    clipped=float(np.sum(p*np.minimum(weights,2)*loss))

    ck(abs(clipped-.2)<1e-14,'clipping introduces explicit bias')

    ck(abs(1/np.sum(p*weights**2)-.36)<1e-14,'weight second moment')

    a=np.linspace(.01,.99,99)

    curves=[]

    for tau in [.25,.5,.75]:

        eff=1/(tau*tau/a+(1-tau)**2/(1-a))

        for ai,ei in zip(a,eff):

            wp=np.array([tau/ai,(1-tau)/(1-ai)])

            ck(abs(float(ei)-1/float(np.dot([ai,1-ai],wp**2)))<1e-12,'ESS weight moments')

        curves.append((f'Target mass {tau:g}',eff))

    panel(18,1,'Poor overlap reduces the weight-only sample fraction','Source probability of region','Asymptotic n_eff / n',a,curves,

          'Binary source mass a and target mass tau give weights tau/a and (1-tau)/(1-a). Curves show 1/E_P[w^2], the population limit of the weight-only diagnostic, not a universal variance or finite-sample guarantee.')

    cvals=np.linspace(.51,1,99);gains=[]

    for c in cvals:

        C=np.array([[c,1-c],[1-c,c]])

        prior=np.array([.7,.3])

        mu=C@prior

        recovered=np.linalg.solve(C,mu)

        ck(np.allclose(recovered,prior,atol=1e-12,rtol=0),'label prior recovery')

        shifted=mu+np.array([-.001,.001])

        shifted_prior=np.linalg.solve(C,shifted)

        gain=(shifted_prior[1]-prior[1])/.001

        ck(abs(gain-1/(2*c-1))<1e-9,'binary inversion sensitivity')

        gains.append(float(1/(2*c-1)))

    C=np.array([[.8,.2],[.2,.8]])

    ck(abs(float((C@np.array([.7,.3]))[1])-.38)<1e-14,'known target prediction fraction')

    singular=np.full((2,2),.5)

    ck(np.linalg.matrix_rank(singular)==1,'uninformative predictions are singular')

    panel(18,2,'An invertible label-shift operator can amplify noise','Symmetric classifier accuracy c','Prior-error / prediction-frequency-error',cvals,

          [('Exact local amplification',gains)],

          'Binary conditional confusion matrix with diagonal c and off-diagonal 1-c. Recovery slope is 1/(2c-1); stable class conditionals and prior support are assumed. The divergence near c=0.5 reflects lost information.')

    x=np.array([0,1]);ya=x.copy();yb=1-x

    ck(np.mean(x!=ya)==0. and np.mean(x!=yb)==1.,'same inputs opposite risks')

    train=np.array([[0,0],[1,1]]); labels=np.array([0,1])

    ck(np.array_equal(train[:,0],labels) and np.array_equal(train[:,1],labels),'two invariant training rules')

    targetx=np.array([[0,1],[1,0]])

    ck(np.mean(targetx[:,0]!=labels)==0. and np.mean(targetx[:,1]!=labels)==1.,'ambiguous stable mechanism')

    points.extend([[18,1,1,.5,1.],[18,2,0,.75,2.]])

    result(18,'Change of measure, clipping bias, prior inversion, and nonidentifiability witnesses',start,

           dict(target_risk=.5,source_risk=.1,clipped_weighted_risk=.2,weight_ESS_fraction=.36,target_prior_positive=.3,observed_positive_fraction=.38,indistinguishable_world_risks=[0.,1.]),

           'Finite known distributions and fixed conditional operators; no learned density ratios, target labels, or neural adaptation benchmark.')



    # 19: finite KL optimum and preference cycle.

    start=count

    margin=np.linspace(-4,4,161)

    panel(19,1,'The same pair margin has beta-dependent loss','Policy-to-reference log-ratio margin','Preferred-pair logistic loss',margin,

          [(f'beta = {b:g}',np.logaddexp(0,-b*margin)) for b in [.5,1.,2.]],

          'Declared DPO pair loss log(1+exp(-beta*m)); m is the difference of sequence log policy-to-reference ratios. These curves are mathematical loss values, not human agreement measurements.')

    betas=np.linspace(.25,4,151);gap=math.log(3)

    policy=1/(1+np.exp(-gap/betas))

    panel(19,2,'Reference regularization selects a response mixture','KL regularization beta','Optimal probability of preferred response',betas,

          [('Equal reference; reward gap log(3)',policy)],

          'Exact two-response KL-regularized optimum with positive equal reference probabilities. Reward difference log(3) is assigned; its scale and the reference are part of this constructed objective.')

    grid=np.linspace(.0001,.9999,10001)

    for b in [.25,.5,1.,2.,4.]:

        opt=1/(1+math.exp(-gap/b))

        def objective(pv):

            return pv*gap-b*(pv*np.log(pv/.5)+(1-pv)*np.log((1-pv)/.5))

        ck(float(objective(opt))>=float(np.max(objective(grid)))-1e-12,'concave policy optimum vs grid')

        ck(abs(float(np.logaddexp(0,0))-math.log(2))<1e-14,'zero margin logistic loss')

    ck(abs(float(np.logaddexp(0,-math.log(3)))-math.log(4/3))<1e-14,'positive nonzero margin exact anchor')
    ck(abs(float(np.logaddexp(0,math.log(3)))-math.log(4))<1e-14,'negative nonzero margin exact anchor')
    ref=np.array([.2,.3,.5]);rewards=np.array([-1.,.2,1.3])

    for b in [.25,1.,4.]:

        raw=ref*np.exp(rewards/b);base=raw/raw.sum()

        for offset in [-3.,0.,2.,5.]:

            raw2=ref*np.exp((rewards+offset)/b);shift=raw2/raw2.sum()

            ck(np.allclose(shift,base,atol=1e-12,rtol=0),'additive reward invariance')

            ck(abs(float(shift.sum())-1)<1e-14,'policy normalization')

    rankings=[('A','B','C'),('B','C','A'),('C','A','B')]

    probs=[]

    for first,second in [('A','B'),('B','C'),('C','A')]:

        prob=sum(r.index(first)<r.index(second) for r in rankings)/3

        ck(abs(prob-2/3)<1e-14,'mixture preference cycle')

        probs.append(prob)

    cycle=float(sum(math.log(p/(1-p)) for p in probs))

    ck(abs(cycle-3*math.log(2))<1e-12 and cycle>0,'scalar log-odds incompatibility')

    ck(abs(1/(1+math.exp(-gap))-.75)<1e-14,'two-response reference point')

    points.extend([[19,1,0,0.,math.log(2)],[19,2,0,1.,.75]])

    result(19,'Finite KL optimum, reward-offset invariance, and cyclic population preferences',start,

           dict(policy_grid_points=10001,beta_values_checked=[.25,.5,1.,2.,4.],optimal_preferred_probability_beta1=.75,cycle_pair_probabilities=probs,cycle_logodds_sum=cycle,scalar_score_required_sum=0.),

           'Assigned finite rewards and three deterministic ranking groups. Grid comparison is a numerical consistency check, not neural convergence or measured human alignment.')



    # 20: compensated rotations, complete 2^3 code enumeration, causal-use witness.

    start=count

    theta=np.linspace(0,math.pi/2,101)

    panel(20,1,'Coordinate contributions change while output does not','Hidden-basis rotation angle (radians)','Contribution to fixed scalar output',theta,

          [('First rotated coordinate',np.cos(theta)**2),('Second rotated coordinate',np.sin(theta)**2),('Total output',np.ones_like(theta))],

          'Hidden input (1,0), readout (1,0), and a matched orthogonal rotation of both vectors. Contributions are cosine squared and sine squared. This compensated linear-interface symmetry need not survive an intervening elementwise nonlinearity.')

    hidden=[np.array([1.,0.]),np.array([0.,1.]),np.array([2.,-3.]),np.array([-.5,.7])]

    v=np.array([1.,-.4])

    for angle in theta:

        R=np.array([[math.cos(angle),-math.sin(angle)],[math.sin(angle),math.cos(angle)]])

        for h in hidden:

            ck(abs(float((R@v)@(R@h)-v@h))<1e-12,'compensated rotation invariance')

    W=np.array([[1.,-.5,-.5],[0.,math.sqrt(3)/2,-math.sqrt(3)/2]])

    patterns=np.array(list(itertools.product([0.,1.],repeat=3)))

    errors_linear=[];errors_relu=[]

    for feat in patterns:

        decoded=W.T@(W@feat);rect=np.maximum(decoded,0)

        k=int(feat.sum())

        er=float(np.sum((rect-feat)**2));el=float(np.sum((decoded-feat)**2))

        expected={0:0.,1:0.,2:.5,3:3.}[k]

        ck(abs(er-expected)<1e-12,'all binary rectified reconstruction cases')

        errors_linear.append(el);errors_relu.append(er)

    ck(np.allclose(W@np.ones(3),W@np.zeros(3),atol=1e-14),'dense sparse encoding collision')

    probabilities=np.linspace(0,1,101);linear=[];rectified=[]

    for prob in probabilities:

        masses=np.array([prob**int(x.sum())*(1-prob)**(3-int(x.sum())) for x in patterns])

        er=float(masses@np.array(errors_relu));el=float(masses@np.array(errors_linear))

        ck(abs(float(masses.sum())-1)<1e-12,'Bernoulli pattern normalization')

        ck(abs(er-(1.5*prob**2+1.5*prob**3))<1e-12,'rectified expected error')

        ck(abs(el-(1.5*prob+1.5*prob**2))<1e-12,'linear expected error')

        rectified.append(er);linear.append(el)

    panel(20,2,'Sparse success can coexist with dense information loss','Independent feature activation probability p','Expected total squared reconstruction error',probabilities,

          [('Rectified decoder',rectified),('Linear decoder',linear)],

          'All eight binary inputs of a fixed three-feature code in two dimensions are enumerated. Directions are 120 degrees apart. Independent Bernoulli(p) activations are assumed; the code is constructed, not trained.')

    observational=np.array([[-1.,-1.],[0.,0.],[1.,1.]])

    ck(np.array_equal(observational[:,0],observational[:,1]),'both coordinates decode observed signal')

    for s,t in observational:

        for d in [-2.,.5,3.]:

            original=s;unused_changed=s;used_changed=s+d

            ck(unused_changed==original,'unused correlated coordinate has no causal effect')

            ck(abs((used_changed-original)-d)<1e-14,'used coordinate changes readout')

    points.extend([[20,1,0,math.pi/4,.5],[20,2,0,.5,.5625]])

    result(20,'Basis invariance, eight-pattern superposition errors, and unused decodable signal',start,

           dict(rotation_angles=101,hidden_inputs=4,binary_patterns=8,rectified_errors_by_active_count={'0':0.,'1':0.,'2':.5,'3':3.},expected_rectified_error_p_half=.5625,dense_and_empty_encode_same=True),

           'Fixed linear interfaces and assigned binary features; no large-model feature extraction, sparse-autoencoder training, or completeness claim.')



    ck(len(panels)==10 and len(results)==5 and len(points)==10,'module artifact counts')

    return dict(panels=panels,results=results,points=points)



if __name__=='__main__':

    out=run()

    print(json.dumps(dict(status='PASS',panels=len(out['panels']),groups=len(out['results']),

                         assertions=sum(r['assertions'] for r in out['results']),points=len(out['points'])),indent=2))