"""Reproduce 40 declared mathematical figures. No network, product benchmark or empirical dataset."""
from pathlib import Path
import json,csv,math,textwrap
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
ROOT=Path(__file__).resolve().parents[1]
TOPICS=json.loads((ROOT/"topics.json").read_text(encoding="utf-8"))
PANELS=[]
def series(label,x,y):return {"label":label,"x":[float(v) for v in x],"y":[float(v) for v in y]}
def panel(n,j,title,xlabel,ylabel,ss,caption,xscale="linear",yscale="linear",**extra):
 slug=TOPICS[n-1][0]
 PANELS.append(dict(article=n,slug=slug,number=j,title=title,xlabel=xlabel,ylabel=ylabel,series=ss,caption=caption,
  xscale=xscale,yscale=yscale,kind="constructed_analytic_example",image=f"figures/{slug}-{j}.png",
  vector=f"figures/{slug}-{j}.svg",data=f"figures/{slug}-{j}.csv",
  alt=title+". "+caption,**extra))
def bin_tail(n,q,p):return sum(math.comb(n,k)*p**k*(1-p)**(n-k) for k in range(q,n+1))
def entropy(p):
 return 0.0 if p in (0,1) else -p*math.log2(p)-(1-p)*math.log2(1-p)
def create_panels():
 p=np.linspace(.5,1,101)
 panel(1,1,"Majorities under independent availability","Per-replica availability","Quorum availability",
  [series(f"{n} replicas",p,[bin_tail(n,n//2+1,x) for x in p]) for n in [3,5,7]],
  "Independent identical replica availability; majority acknowledgements only. Persistent voting rules, partitions and common-cause failures are outside this probability model.")
 t=np.linspace(0,60,121)
 panel(1,2,"An ideal quorum collector","Elapsed time (ms)","Probability quorum has replied",
  [series(f"{n} replicas",t,[bin_tail(n,n//2+1,1-math.exp(-x/10)) for x in t]) for n in [3,5,7]],
  "Independent exponential acknowledgement times with mean 10 ms. No leader election, queueing, retransmission or correlated delay; this is not measured Raft latency.")
 p=np.linspace(0,1,101)
 panel(2,1,"Simultaneous withdrawals can violate an invariant","Per-reviewer withdrawal probability","Probability all withdraw",
  [series(f"{n} reviewers",p,p**n) for n in [2,4,8]],
  "All n transactions read the same initially valid snapshot, independently withdraw with probability p and update distinct rows. At least two reviewers; no serializable guard.")
 a=np.linspace(0,.95,96)
 panel(2,2,"Retries consume an explicit budget","Independent abort probability","Expected attempts per success",[series("Geometric retry model",a,1/(1-a))],
  "Independent identical abort probability and unlimited retries. Contention feedback, deadlines and permanent errors are excluded.")
 p=np.linspace(0,.9,91)
 panel(3,1,"Lost acknowledgements amplify unguarded effects","Acknowledgement-loss probability","Expected extra effects",[series("Effect repeated on every retry",p,p/(1-p))],
  "Every attempt performs the effect before an independently lost acknowledgement; retries continue until acknowledged. Durable deduplication would change this model.")
 h=np.linspace(0,24,97)
 panel(3,2,"Receipt retention has a storage cost","Retention horizon (hours)","Receipt payload (GB, decimal)",
  [series(f"{r:,} operations/s",h,h*3600*r*128/1e9) for r in [100,1000,10000]],
  "128 payload bytes per completed operation, constant rate and complete retention. Indexes, replication, allocator overhead and garbage-collection delay are excluded.")
 t=np.linspace(0,40,101)
 panel(4,1,"Propagation under a Poisson contact model","Time since update (minutes)","Probability update is not propagated",
  [series(f"Contact rate {r}/min",t,np.exp(-r*t)) for r in [.1,.5,1]],
  "One successful contact suffices and contacts form a Poisson process. This is a single-link survival model, not a bound for an arbitrary CRDT network.")
 u=np.linspace(0,2_000_000,101)
 panel(4,2,"Metadata growth before compaction","Retained operations (millions)","Metadata payload (MB, decimal)",
  [series(f"{m} bytes/operation",u/1e6,u*m/1e6) for m in [16,32,64]],
  "Linear accounting with declared per-operation metadata. No compaction, compression or implementation-specific overhead; not measured CRDT memory.")
 b=np.arange(1,33)
 panel(5,1,"Bit width expands the raw state space","State bits","Possible bit patterns",[series("2 to the power of b",b,2.0**b)],
  "All bit patterns counted, without reachability reduction or symmetry. The count is not a model-checker runtime prediction.",yscale="log")
 N=np.logspace(0,6,121)
 panel(5,2,"Random testing can miss one rare bad state","Independent uniform samples","Probability bad state is never sampled",
  [series(f"{b}-bit state",N,np.exp(N*np.log1p(-2.0**(-b)))) for b in [8,12,16]],
  "Exactly one bad state among all equally likely b-bit patterns, sampling with replacement. This model excludes systematic tests and biased sampling.",xscale="log")
 panel(6,1,"Reachable outcomes depend on the memory model","Model","Distinct load-result pairs",
  [series("Finite store-buffering model",[0,1,2],[3,4,3])],
  "Two threads each store then load the other location. Sequential consistency allows three pairs; the explicit store-buffer model allows four; flushing before loads restores three. Counts, not frequencies.",
  plot="bar",ticklabels=["Sequential\nconsistency","Store\nbuffers","Buffers +\nflush-before-load"])
 m=np.arange(1,16)
 panel(6,2,"Interleavings grow before state reduction","Events in each of two threads","Order-preserving interleavings",[series("Central binomial coefficient",m,[math.comb(2*int(x),int(x)) for x in m])],
  "Two fixed per-thread orders of m distinct events, with no cross-thread constraints or equivalent-state reduction.",yscale="log")
 x=np.arange(256)
 panel(7,1,"One boundary value defeats an optimization","Unsigned 8-bit input x","Value after addition",
  [series("Mathematical x + 1",x,x+1),series("8-bit modular x + 1",x,(x+1)%256)],
  "Unsigned 8-bit modular arithmetic. At x = 255 the increment wraps to zero, so replacing (x + 1 > x) by true is invalid under these semantics.")
 k=np.arange(1,21);a=10.0**k
 panel(7,2,"Evaluation order changes floating-point results","Exponent k in a = 10^k","Expression result",
  [series("(a - a) + 1",k,(a-a)+1),series("a + (-a + 1)",k,a+(-a+1))],
  "IEEE-style binary64 calculations in the supplied NumPy environment. These constructed expressions are not a benchmark of LLVM or CompCert.")
 g=np.array([8,16,32,64,128,256])
 panel(8,1,"One tag bit is not the complete memory cost","Tagged granule (bytes)","Raw tag storage (%)",[series("One bit per granule",g,100/(8*g))],
  "Raw bit accounting only; physical storage granularity, ECC, capability width and metadata caches are excluded. Not a measured CHERI overhead.",xscale="log")
 r=np.linspace(0,1000,101)
 panel(8,2,"Delayed reuse retains memory","Freeing rate (MB/s, decimal)","Quarantine payload (MB, decimal)",
  [series(f"Delay {ms} ms",r,r*ms/1000) for ms in [1,5,10]],
  "Steady freeing rate times a fixed quarantine delay. Reallocation, fragmentation and bursty frees are excluded; this is not a Cornucopia benchmark.")
 p=np.linspace(0,.5,101)
 panel(9,1,"A noisy bit can still carry information","Independent bit-flip probability","Mutual information (bits)",[series("Uniform bit through binary symmetric channel",p,[1-entropy(x) for x in p])],
  "Uniform secret bit and independent symmetric errors. This hypothetical channel is not a measured processor side channel.")
 n=np.arange(1,32,2)
 panel(9,2,"Repeated observations can reduce decision error","Odd number of independent observations","Majority decision error",
  [series(f"Bit-flip probability {p}",n,[bin_tail(int(v),int(v)//2+1,p) for v in n]) for p in [.1,.25,.4]],
  "Independent repeated observations of the same bit and majority decoding. Correlated noise and nonstationary secrets invalidate this calculation.",yscale="log")
 k=np.arange(1,1001)
 panel(10,1,"Repeated Gaussian releases accumulate a bound","Number of releases","Sufficient epsilon bound",
  [series(f"Noise sigma = {s}",k,k/(2*s*s)+2*np.sqrt(k/(2*s*s)*math.log(1e6))) for s in [4,8,16]],
  "Sensitivity 1, delta = 10^-6, no subsampling. RDP composition followed by continuous-order optimization gives a sufficient bound, not optimal analytical calibration.",xscale="log")
 r=np.arange(1,1001)
 panel(10,2,"Averaging reduces independent release noise","Repeated releases of one statistic","Standard deviation of the average",[series("Independent noise sigma = 8",r,8/np.sqrt(r))],
  "The same fixed statistic is released repeatedly with independent Gaussian noise. Accuracy improves while privacy cost accumulates; correlated or reused noise differs.",xscale="log",yscale="log")
 b=np.arange(8,129)
 panel(11,1,"Polynomial identity testing needs a field budget","Approximate field-size bits b","Root-probability upper bound",
  [series(f"Polynomial degree {d:,}",b,np.minimum(1,d/2.0**b)) for d in [1,256,65536]],
  "Degree/field-size union bound for a nonzero polynomial and uniform challenge. The horizontal variable uses idealized cardinality 2^b; it is not a security estimate for a full proof system.",yscale="log")
 N=2.0**np.arange(4,21)
 panel(11,2,"Succinct verification does not erase prover work","Problem size N","Declared abstract operation count",
  [series("N",N,N),series("N log2 N",N,N*np.log2(N)),series("log2 N",N,np.log2(N))],
  "Illustrative growth functions with coefficient one. Curves are not measured or asserted complexities of Groth16, PLONK or STARK implementations.",xscale="log",yscale="log")
 w=2.0**np.arange(4,17)
 panel(12,1,"Width sets a Count-Min error scale","Counter width per row","Epsilon in width >= e/epsilon",[series("e / width",w,math.e/w)],
  "Standard insertion-only Count-Min sizing relation. A probability guarantee additionally needs appropriate independent hashing and depth; the toy appendix hashes do not establish that guarantee.",xscale="log",yscale="log")
 m=2.0**np.arange(4,19)
 panel(12,2,"HyperLogLog trades registers for precision","Number of registers","Asymptotic relative standard error (%)",[series("104 / sqrt(m)",m,104/np.sqrt(m))],
  "Classical approximate error scale under suitable hashing and asymptotic conditions. Small-cardinality corrections, bias and implementation choices are excluded.",xscale="log",yscale="log")
 x=np.linspace(0,1,257)
 panel(13,1,"A stale rank model can miss its search window","Query coordinate x","Continuous rank target",
  [series("Original uniform rank: 1024 x",x,1024*x),series("Shifted rank: 1024 sqrt(x)",x,1024*np.sqrt(x))],
  "Declared continuous population targets for N = 1025. Their maximum gap is 256 at x = 0.25; finite lower-bound ranks are checked separately in the appendix.")
 e=np.arange(1,257)
 panel(13,2,"A bounded residual narrows a comparison budget","Certified absolute rank error (slots)","Illustrative comparison budget",
  [series("ceil(log2(2 epsilon + 1))",e,np.ceil(np.log2(2*e+1))),series("ceil(log2(1025))",e,np.full(len(e),11))],
  "Approximate binary-search comparison budgets, not timed implementations. A window is sound only if its residual bound or fallback validation is valid.")
 c=np.arange(1,501)
 ss=[]
 for alpha in [.7,1.1,1.5]:
  probs=np.arange(1,1001,dtype=float)**(-alpha);probs/=probs.sum()
  ss.append(series(f"Zipf exponent {alpha}",c,np.cumsum(probs)[:500]))
 panel(14,1,"An ideal static cache under a known popularity law","Cache capacity (equal-size objects)","Oracle static hit probability",ss,
  "1,000 equal-size objects with independent stationary Zipf requests and known popularity. Top-C storage is an oracle reference, not a measured LRU, ARC, TinyLFU or LeCaR hit rate.")
 t=np.linspace(0,500,101)
 panel(14,2,"Forgetting controls adaptation speed","Time since last observation (arbitrary units)","Decayed historical frequency",
  [series(f"Half-life {H}",t,64*2.0**(-t/H)) for H in [10,50,200]],
  "Constructed exponential decay from count 64 with no new observations. A shorter half-life forgets more quickly; the curves do not implement a cited cache policy.")
 n=np.unique(np.round(np.logspace(0,3,121)).astype(int))
 panel(15,1,"Fan-out amplifies a fixed branch tail","Required branches","Probability at least one branch exceeds threshold",
  [series(f"Branch tail probability {p}",n,1-(1-p)**n) for p in [.001,.01,.05]],
  "Independent identical branches; all must complete and every branch uses the same fixed latency threshold. Correlation and quorum completion are excluded.",xscale="log")
 ss=[]
 for h in [0,.1,.25]:
  rho=np.linspace(0,.98/(1+h),121)
  ss.append(series(f"Full-work duplicate fraction {h}",rho,1/(1-(1+h)*rho)))
 panel(15,2,"Extra work erodes queue capacity","Original utilization rho","Mean response / mean service time",ss,
  "M/M/1 reference model with full additional work on fraction h of requests and fixed service rate. Curves stop below instability; successful early cancellation is excluded.",yscale="log")
 p=np.logspace(-4,-1,121)
 panel(16,1,"An MDS threshold under independent unavailability","Per-fragment unavailability","Probability reconstruction threshold is missed",
  [series(f"n={n}, k=6; storage {n/6:.2f}x",p,[bin_tail(n,n-6+1,x) for x in p]) for n in [8,10,12]],
  "Instantaneous independent identical unavailability and ideal MDS fragments. This is not annual data-loss probability; correlated failures, keys and repair dynamics are excluded.",xscale="log",yscale="log")
 d=np.arange(6,14)
 panel(16,2,"More helpers lower an information-flow bound","Repair helpers d (k = 6, n = 14)","Download / one-fragment size",
  [series("Minimum-storage lower bound",d,d/(d-6+1)),series("Download k full fragments",d,np.full(len(d),6))],
  "Regenerating-code minimum-storage bound with alpha = M/k. Exact-repair achievability, subpacketization and disk-read overhead are additional requirements; not a Hitchhiker measurement.")
 q=np.linspace(0,2000,101)
 panel(17,1,"Buffer bytes translate into delay","Queued payload (kB, decimal)","Fluid queue delay (microseconds)",
  [series(f"Link {C} Gbit/s",q,q*1000*8/(C*1e9)*1e6) for C in [10,40,100]],
  "One fixed-capacity fluid queue; propagation, scheduler effects and packet overhead are excluded. Decimal kB and Gbit/s are used.")
 share=np.linspace(0,1,101)
 panel(17,2,"Full utilization can conceal unequal allocation","First flow share of total rate","Jain equality index",
  [series("Two flows, fixed positive total",share,1/(2*(share**2+(1-share)**2)))],
  "Equal weights and two nonnegative rates summing to one. This describes rate equality, not weighted fairness or deadline satisfaction.")
 tau=np.linspace(0,300,121)
 panel(18,1,"Idle retention changes a simple cold probability","Idle timeout (seconds)","Probability an arrival finds a cold environment",
  [series(f"Mean gap {g} s",tau,np.exp(-tau/g)) for g in [60,10,1]],
  "Poisson arrivals, one environment, negligible service time, timeout reset at every arrival and no other eviction. Concurrency and scale-out are excluded.")
 n=np.arange(1,101)
 panel(18,2,"A function graph amplifies cold components","Required components","Probability at least one component is cold",
  [series(f"Component cold probability {p}",n,1-(1-p)**n) for p in [.01,.05,.2]],
  "Independent identical cold events and all-component completion. Shared placement, dependency caches and correlated eviction are excluded.")
 T=np.arange(128,16385,128)
 panel(19,1,"Key-value state grows with retained tokens","Context per request (tokens)","KV state only (GiB)",
  [series(f"{batch} request"+("s" if batch!=1 else "")+"; 8 KV heads",T,2*32*8*128*2*T*batch/2**30) for batch in [1,8,16]],
  "32 layers, 8 key-value heads, head dimension 128 and 2 bytes/scalar. Excludes weights, activations, workspace, metadata and padding; no named GPU capacity is implied.")
 B=np.array([4,8,16,32,64])
 panel(19,2,"Page tails retain a bounded allocation waste","Page size (token slots)","Mean unused slots per request",
  [series("Paged allocation",B,(B-1)/2),series("Reserve 512 slots per request",B,np.full(len(B),255.5))],
  "Uniform integer lengths 1 through 512. Page sizes divide 512. Counts include only tail slots, not metadata, external fragmentation or sharing; not a vLLM benchmark.")
 J=np.logspace(0,4,121)
 panel(20,1,"More inner iterations cannot remove run variation","Iterations per independent run","Standard error (arbitrary latency units)",
  [series(f"{R} independent runs",J,np.sqrt(4/R+100/(R*J))) for R in [5,20]],
  "Constructed independent run effects with standard deviation 2 and within-run noise with standard deviation 10. This is known-model uncertainty, not a fitted confidence interval.",xscale="log")
 w=np.linspace(0,1,101)
 panel(20,2,"Workload weights change the aggregate comparison","Hard-workload fraction w","Mean latency (arbitrary units)",
  [series("System A: easy 1, hard 10",w,1+9*w),series("System B: easy 2, hard 12",w,2+10*w)],
  "Constructed workloads: A is faster in each class. Comparing A at w=0.9 (9.1) with B at w=0.1 (3.0) reverses the aggregate because the workloads differ.")
def render():
 out=ROOT/"figures";out.mkdir(exist_ok=True)
 plt.rcParams.update({"font.family":"DejaVu Sans","font.size":10,"axes.spines.top":False,"axes.spines.right":False,
  "axes.edgecolor":"#a8ada5","axes.labelcolor":"#263b3b","xtick.color":"#465957","ytick.color":"#465957",
  "svg.hashsalt":"auturnal-cs-20-v1","figure.facecolor":"#f7f6f0","axes.facecolor":"#f7f6f0","savefig.facecolor":"#f7f6f0"})
 colors=["#205e67","#ad683c","#657fa3","#85566f"];styles=["-","--","-.",":"]
 for p in PANELS:
  fig,ax=plt.subplots(figsize=(7.2,4.8),dpi=180)
  fig.subplots_adjust(left=.13,right=.97,top=.78,bottom=.18)
  fig.text(.13,.94,f"{p['article']:02}.{p['number']}  /  DECLARED MODEL",fontsize=8.5,color="#8b7144",weight="bold")
  fig.text(.13,.875,"\n".join(textwrap.wrap(p["title"],52)),fontsize=13,color="#233d3d",weight="bold",va="center")
  for i,s in enumerate(p["series"]):
   if p.get("plot")=="bar":ax.bar(s["x"],s["y"],width=.55,color=colors[0]);
   else:ax.plot(s["x"],s["y"],label=s["label"],lw=2,color=colors[i%4],ls=styles[i%4])
  ax.set(xlabel=p["xlabel"],ylabel=p["ylabel"],xscale=p["xscale"],yscale=p["yscale"])
  if p.get("plot")=="bar":ax.set_xticks(p["series"][0]["x"],p["ticklabels"]);ax.set_yticks(range(5))
  ax.grid(axis="y",alpha=.19);ax.set_axisbelow(True)
  if p.get("plot")!="bar":ax.legend(frameon=True,facecolor="#f7f6f0",edgecolor="#dedfd7",fontsize=8,loc="best")
  fig.text(.13,.035,"CONSTRUCTED EXAMPLE  /  See article for assumptions and downloadable data.",fontsize=7,color="#66736d")
  fig.savefig(ROOT/p["image"],dpi=180,metadata={"Software":"Auturnal declared-model generator"})
  fig.savefig(ROOT/p["vector"],metadata={"Date":None,"Creator":"Auturnal declared-model generator"})
  plt.close(fig)
  with (ROOT/p["data"]).open("w",encoding="utf-8",newline="") as f:
   w=csv.writer(f);w.writerow(["series","x","y"])
   for s in p["series"]:
    w.writerows((s["label"],format(x,".17g"),format(y,".17g")) for x,y in zip(s["x"],s["y"]))
 (out/"panels.json").write_text(json.dumps(PANELS,ensure_ascii=False,indent=2,allow_nan=False)+"\n",encoding="utf-8")
if __name__=="__main__":
 create_panels();assert len(PANELS)==40
 render();print(json.dumps({"figures":len(PANELS),"numeric_points":sum(len(s["x"]) for p in PANELS for s in p["series"])}))