"""Finite analytical examples for training and optimization. No learned-model benchmark."""
from pathlib import Path
import itertools,json,math
import numpy as np
ROOT=Path(__file__).resolve().parents[1]
TOPICS=json.loads((ROOT/'topics.json').read_text(encoding='utf-8'))
def run():
 panels=[];results=[];points=[];count=0
 def check(v,msg):
  nonlocal count
  assert bool(v),msg
  count+=1
 def series(label,x,y):return {'label':label,'x':np.asarray(x,dtype=float).tolist(),'y':np.asarray(y,dtype=float).tolist()}
 def panel(n,j,title,xlabel,ylabel,ss,caption,xscale='linear',yscale='linear'):
  slug=TOPICS[n-1][0];stem=f'figures/{slug}-{j}'
  panels.append(dict(article=n,number=j,slug=slug,title=title,xlabel=xlabel,ylabel=ylabel,series=ss,caption=caption,xscale=xscale,yscale=yscale,kind='constructed_analytic_example',image=stem+'.png',vector=stem+'.svg',data=stem+'.csv',alt=title+'. '+caption))
 def finish(n,start,label,witness,scope):
  results.append(dict(article=n,slug=TOPICS[n-1][0],check=label,assertions=count-start,witness=witness,scope=scope,status='PASS'))
 start=count;g=np.array([1.,3.,5.,7.]);vv=[]
 for b in range(1,5):
  means=[np.mean(z) for z in itertools.combinations(g,b)];v=float(np.var(means));vv.append(v)
  check(np.isclose(v,5*(4-b)/(b*3)),f'finite population B={b}')
 eta=.5;b=4;h=1;V=eta/(b*h*(2-eta*h));v=0.
 for _ in range(200):v=(1-eta*h)**2*v+eta**2/b
 check(np.isclose(v,1/12),'independent recurrence variance')
 check(np.isclose(V,1/12),'closed form variance')
 H=np.array([[2.,.3],[.3,1.]])
 S=np.array([[1.,.2],[.2,.7]]);eta=.2;A=np.eye(2)-eta*H
 C=np.linalg.solve(np.eye(4)-np.kron(A,A), (eta**2*S/4).reshape(4)).reshape((2,2));C2=np.zeros((2,2))
 for _ in range(400):C2=A@C2@A.T+eta**2*S/4
 check(np.allclose(C,C2,atol=1e-12),'Lyapunov iteration vs direct solution')
 check(np.min(np.linalg.eigvalsh(C))>0,'positive covariance')
 check(abs(1-2.1)**20>1,'unstable step grows')
 check(not np.isclose(.1/(2-.1),.4/(4*(2-.4))),'fixed eta/B not exact')
 x=np.arange(1,5)
 panel(1,1,'Finite sampling exhausts the noise','Batch size B (four gradients)','Variance of batch mean',[series('Without replacement',x,vv),series('With replacement',x,5/x)],'Exact population {1,3,5,7}; variance denominator N. Uniform distinct subsets versus independent draws; gradient units squared.')
 panels[-1]['xticks']=[1,2,3,4];panels[-1]['marker']='o'
 e=np.linspace(.02,1.9,95)
 panel(1,2,'Finite steps change stationary variance','Learning rate eta (h = 1)','Stationary parameter variance',[series('Exact discrete recursion',e,e/(4*(2-e))),series('Small-step approximation',e,e/8)],'Constructed quadratic with h=1, additive noise variance 1/B and B=4. Independent noise, constant step, eta<2; no neural-network measurement.')
 points.extend([[1,1,0,2,5/3],[1,2,0,.5,1/12]])
 finish(1,start,'Sampling and discrete covariance identities',{'finite_population_variances':vv,'scalar_stationary_variance':V,'matrix_covariance':C.tolist()},'Enumerates all nonempty subset sizes of a four-gradient population and verifies one stable 2D covariance. No posterior-calibration or generalization evidence.')
 start=count
 Z=np.array([[1.,0.],[0.,2.]])
 wh=np.array([1.,.5]);gm=np.min(Z@wh)/np.linalg.norm(wh)
 check(np.isclose(gm,2/np.sqrt(5)),'optimal constructed margin')
 for c in [.01,.1,1,10,100]:
  check(np.isclose(np.min(Z@(c*wh))/np.linalg.norm(c*wh),gm),'positive scale invariance')
  check(np.all(Z@(c*wh)>0),'zero classification error')
 w=np.array([.7,.4]);grad=-(Z.T@(1/(1+np.exp(Z@w))))/2
 AA=np.diag([2.,.5]);theta=np.linalg.solve(AA,w);gtheta=AA.T@grad;eta=.1
 check(np.allclose(AA@(theta-eta*gtheta),w-eta*AA@AA.T@grad),'coordinate chain rule')
 check(not np.allclose(w-eta*grad,AA@(theta-eta*gtheta)),'different induced update')
 # Independent finite-difference logistic gradient.
 loss=lambda u:float(np.logaddexp(0,-Z@u).mean())
 fd=np.array([(loss(w+1e-6*np.eye(2)[j])-loss(w-1e-6*np.eye(2)[j]))/2e-6 for j in range(2)])
 check(np.allclose(grad,fd,atol=1e-9),'logistic gradient')
 wt=np.zeros(2);trace=[]
 for t in range(10000):
  wt+=.5*(Z.T@(1/(1+np.exp(Z@wt))))/2
  if t in [0,9,99,999,9999]:trace.append([t+1,loss(wt),float(np.min(Z@wt)/np.linalg.norm(wt))])
 check(trace[-1][1]<trace[0][1],'finite run loss decreases')
 for score in np.linspace(-5,5,101):
  check((np.logaddexp(0,score)+np.logaddexp(0,-score))/2>=math.log(2)-1e-12,'contradictory pair optimum')
 c=np.linspace(.1,10,100)
 panel(2,1,'Confidence can rise without moving the boundary','Positive weight multiplier c','Mean logistic loss',[series('Two equal functional margins',c,np.logaddexp(0,-c))],'Signed vectors (1,0),(0,2), base w=(1,1/2). Positive rescaling preserves every classification and the normalized margin; loss alone changes.')
 theta=np.linspace(1,89,89);rad=theta*np.pi/180
 panel(2,2,'Zero error does not select a direction','Weight direction angle (degrees)','Minimum normalized training margin',[series('Unit direction (cos theta, sin theta)',theta,np.minimum(np.cos(rad),2*np.sin(rad)))],'Two signed constraints; all displayed directions classify both points correctly. Exact optimal angle arctan(1/2) is about 26.565 degrees; no test population specified.')
 points.extend([[2,1,0,1,math.log1p(math.exp(-1))],[2,2,0,45,1/math.sqrt(2)]])
 finish(2,start,'Margins, coordinate changes, and contradictory labels',{'maximum_margin':float(gm),'trajectory_step_loss_margin':trace,'direct_gradient':grad.tolist()},'Finite two-constraint geometry, explicit gradient steps and 101 contradictory-pair scores. Finite descent is not proof of the cited asymptotic theorems.')
 start=count
 x=np.array([1.,2.,4.]);d=3;P=np.eye(d)-np.ones((d,d))/d;cc=P@x;eps=.1;v=cc@cc/d;s=np.sqrt(v+eps)
 norm=lambda xx:(P@xx)/np.sqrt((P@xx)@(P@xx)/d+eps)
 J=P/s-np.outer(cc,cc)/(d*s**3)
 check(np.allclose(norm(x),norm(x+10)),'shift invariant')
 check(np.allclose(J@np.ones(d),0),'mean null direction')
 check(np.allclose(J@cc,eps*cc/s**3),'radial direction')
 fd=np.column_stack([(norm(x+1e-6*np.eye(d)[j])-norm(x-1e-6*np.eye(d)[j]))/2e-6 for j in range(d)])
 check(np.allclose(fd,J,atol=1e-9),'normalization Jacobian')
 for aa in [.01,.1,1,10]:
  check(np.allclose(norm(aa*x),aa*cc/np.sqrt(aa**2*v+eps)),'positive scaling expression')
 check(not np.allclose(norm(.1*x),norm(x)),'epsilon breaks scale invariance')
 bn=lambda b:-b/np.sqrt(b*b+4*.01)
 check(bn(2)<0 and bn(-2)>0,'companion reversal')
 scale=np.logspace(-3,2,101)
 panel(3,1,'A stabilizer limits scale invariance','Positive input scale a','Output norm / norm at a = 1',[series(f'epsilon = {ee}',scale,scale*np.sqrt(1+ee)/np.sqrt(scale**2+ee)) for ee in [.01,.1]],'Centered vector has variance v=1. Norm ratio is dimensionless. Mean shift cancels exactly; positive scale cancels only in the zero-epsilon nonconstant limit.',xscale='log')
 bb=np.linspace(-4,4,161)
 panel(3,2,'The same example changes with its companion','Companion scalar value b','Normalized value of fixed sample zero',[series(f'epsilon = {ee}',bb,-bb/np.sqrt(bb**2+4*ee)) for ee in [.01,.1]],'Two-example batch {0,b}, population variance, no learned gain or bias. Only the companion changes; sign reversal follows directly from the shared statistics.')
 points.extend([[3,1,0,1,1],[3,2,0,2,-2/math.sqrt(4.04)]])
 finish(3,start,'Normalization invariance and Jacobian checks',{'jacobian':J.tolist(),'max_finite_difference_error':float(np.max(np.abs(J-fd))),'companion_outputs':[float(bn(2)),float(bn(-2))]},'Three-dimensional normalization and two-sample batch normalization with stated population variance and epsilon; excludes network training and mixed-precision implementation.')
 start=count
 for L in [1,2,4,8,16,32,64,128]:
  for lam in [.5,1.,2.]:
   h0=1.7;hh=h0
   for _ in range(L):hh*=1-lam/L
   check(np.isclose(hh,(1-lam/L)**L*h0),'Euler recurrence')
   der=-h0*(1-lam/L)**(L-1)
   fun=lambda la:(1-la/L)**L*h0
   fd=(fun(lam+1e-6)-fun(lam-1e-6))/2e-6
   check(np.isclose(der,fd,atol=2e-9),'general h0 Euler derivative')
 check(abs(1-2.1)**10>1,'unstable discrete multiplier')
 check(np.isclose((1-1)**1,0),'contractive but noninvertible update')
 check(math.exp(-1)*math.sqrt(17)>1,'non-normal transient amplification')
 L=np.arange(1,129)
 panel(4,1,'Stability alone does not ensure accuracy','Equal steps L over T = 1','Absolute multiplier error',[series('Euler versus exp(-1)',L,np.abs((1-1/L)**L-math.exp(-1)))],'Scalar dh/dt=-h, fixed horizon T=1, fixed field. Error concerns the multiplier h(T)/h(0), not a learned task outcome.',xscale='log',yscale='log')
 t=np.linspace(0,5,101)
 panel(4,2,'Stable eigenvalues allow transient growth','Dimensionless time t','State norm / initial norm',[series(f'Off-diagonal coupling k = {k}',t,np.exp(-t)*np.sqrt(1+k*k*t*t)) for k in [0,2,4]],'A=[[-1,k],[0,-1]], h(0)=(0,1). Eigenvalues are -1 for every curve; non-normal coupling causes finite-horizon amplification.')
 points.extend([[4,1,0,1,math.exp(-1)],[4,2,2,1,math.sqrt(17)/math.e]])
 finish(4,start,'Discrete stability, sensitivity, and transient growth',{'general_initial_state':1.7,'non_normal_norm_t1_k4':math.sqrt(17)/math.e},'24 scalar Euler cases with recurrence and finite-difference sensitivities plus explicit instability, noninvertibility and non-normal examples. Not a residual-network benchmark.')
 start=count
 def optimum(a,b,al,be,C,k=1):
  N=(al*a/(be*b))**(1/(al+be))*(C/k)**(be/(al+be));return N,C/(k*N)
 for a,b,al,be in [(1,1,.5,.5),(2,3,.3,.7),(.5,2,.8,.2)]:
  for C in [100,10000,1e6]:
   N,D=optimum(a,b,al,be,C)
   check(np.isclose(al*a*N**(-al),be*b*D**(-be)),'weighted marginals')
   check(al*al*a*N**(-al)+be*be*b*D**(-be)>0,'positive log-size curvature')
   f=lambda z:a*z**(-al)+b*(C/z)**(-be)
   check(f(N)<f(N*.9) and f(N)<f(N*1.1),'local allocation comparison')
   grid=N*np.logspace(-3,3,10001)
   check(np.all(f(grid)>=f(N)-1e-12),'independent dense log-grid minimum')
 N,D=optimum(1,1,.5,.5,10000)
 check(np.isclose(N,100) and np.isclose(D,100),'symmetric optimum')
 cap=25;Nc=10000/cap;f=lambda z:z**(-.5)+(10000/z)**(-.5)
 check(f(Nc)<f(Nc*1.1),'active data cap at minimum feasible N')
 def inferN(q,C=10000.):
  lo=0.;hi=math.sqrt(C) if q==0 else min(math.sqrt(C),C/q)
  for _ in range(100):
   mid=(lo+hi)/2
   val=C-q*mid-(C*mid)**(2/3)
   if val>0:lo=mid
   else:hi=mid
  return (lo+hi)/2
 qs=np.linspace(0,100,101);Ns=np.array([inferN(q) for q in qs])
 for q,n in zip(qs,Ns):
  dd=10000/n-q;check(dd>0,'inference feasible')
  check(np.isclose(10000,np.sqrt(n)*dd**1.5,rtol=1e-10),'inference derivative zero')
  ff=lambda nn:nn**(-.5)+(10000/nn-q)**(-.5)
  check(ff(n)<=ff(n*.99) and ff(n)<=ff(n*1.01),'inference nearby allocation')
 check(np.all(np.diff(Ns)<0),'inference charge lowers optimum in chosen model')
 nn=np.logspace(0,3,121)
 panel(5,1,'Allocation changes loss at fixed compute','Normalized parameter count N','Constructed excess loss',[series(f'Product budget C = {C:g}',nn,1/np.sqrt(nn)+np.sqrt(nn/C)) for C in [100,10000,1e6]],'a=b=k=1, alpha=beta=1/2, D=C/N. Dimensionless resource units; no coefficients fitted to observed model training.',xscale='log',yscale='log')
 panel(5,2,'Inference accounting shifts the optimum','Normalized inference charge q','Optimal normalized N',[series('C total = 10,000',qs,Ns)],'Minimize N^(-1/2)+D^(-1/2) subject to N(D+q)=10,000 and D>0. qN is an illustrative lifetime charge; all units normalized, not serving measurements.')
 points.extend([[5,1,1,100,.2],[5,2,0,0,100]])
 finish(5,start,'Allocation derivatives, data cap, and inference budget',{'symmetric_optimum':[N,D],'data_cap_optimum':[Nc,cap],'inference_q100_optimum':[float(Ns[-1]),float(10000/Ns[-1]-100)]},'Nine power-law coefficient/budget cases and 101 inference charges with exact stationary conditions. Does not validate the loss or cost model empirically.')
 return {'panels':panels,'results':results,'points':points}
if __name__=='__main__':
 r=run();print(json.dumps({'groups':len(r['results']),'assertions':sum(x['assertions'] for x in r['results']),'panels':len(r['panels'])}))