"""Reproduce the 80 analytical figures. All inputs are constructed, not observations.
Run: python analysis/reproduce.py
Dependencies: numpy>=2.0, matplotlib>=3.8. No network access is used.
"""
from pathlib import Path
import csv, hashlib, json, math, textwrap
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt

ROOT=Path(__file__).resolve().parents[1]
L=lambda a,b,n=201:np.linspace(a,b,n)
G=lambda a,b,n=201:np.geomspace(a,b,n)
PANELS=[]
def add(i,n,title,xlabel,ylabel,x,ys,caption,xscale="linear",yscale="linear"):
    slug=TOPICS[i-1][0]
    PANELS.append(dict(article=i,slug=slug,number=n,title=title,xlabel=xlabel,ylabel=ylabel,
        xscale=xscale,yscale=yscale,caption=caption,kind="constructed_analytic_example",
        series=[dict(label=label,x=np.asarray(x,dtype=float).tolist(),y=np.asarray(y,dtype=float).tolist()) for label,y in ys]))

def models():
    x=L(1,10)
    add(1,1,"A lever rule changes the phase fraction","Mean concentration (arbitrary units)","Dense-phase volume fraction",(x),
        [("Dilute = 1; dense = 10",(x-1)/9)],
        "Two coexisting concentrations are held at 1 and 10 arbitrary units. Mass balance gives the dense-phase fraction (mean minus 1)/9. This equilibrium construction does not establish phase separation in a cell.")
    t=L(0,20)
    add(1,2,"Recovery kinetics do not identify a mechanism","Time (arbitrary units)","Normalized fluorescence",t,
        [("Mobile 1.0; rate 0.2",1-np.exp(-.2*t)),("Mobile 0.6; rate 0.2",.6*(1-np.exp(-.2*t))),("Mobile 1.0; rate 1.0",1-np.exp(-t))],
        "Exponential recovery with a mobile fraction of 0.6 or 1 and a rate of 0.2 or 1 per time unit. The curves isolate mobility and recovery rate; they are not fitted FRAP measurements.")
    p=L(0,1)
    add(2,1,"A changing mixture can change the clock input","Fraction of cell type A","Bulk methylation fraction",p,[("Type A = 0.2; type B = 0.8",.2*p+.8*(1-p))],
        "Cell-type methylation is fixed at 0.2 and 0.8. A change in composition changes the bulk assay without any within-cell methylation change; no clock coefficients are applied.")
    x=L(-15,15)
    add(2,2,"An extreme baseline is not an intervention","Baseline deviation from mean (score units)","Expected repeat deviation (score units)",x,
        [("Repeat correlation 0.8",.8*x),("Perfect repeat correlation",x)],
        "The conditional expectation of a repeated centered score is rho times the baseline deviation under a joint Gaussian model with equal marginal variance. The model has rho = 0.8 and no treatment effect.")
    t=L(0,40)
    add(3,1,"Persistence depends on removal as well as input","Time (arbitrary units)","Senescent burden (arbitrary units)",t,
        [("Removal rate 0.1",(1-np.exp(-.1*t))/.1),("Removal rate 0.2",(1-np.exp(-.2*t))/.2)],
        "Constant input I = 1 burden unit per time unit; initial burden is zero. The linear clearance model approaches I/k. It omits tissue-specific benefits and damage.")
    s=L(0,30)
    add(3,2,"A benefit-cost balance need not be monotonic","Burden (arbitrary units)","Model contribution (arbitrary units)",s,
        [("Saturating benefit",4*s/(2+s)),("Accumulating cost",.2*s),("Benefit minus cost",4*s/(2+s)-.2*s)],
        "Illustrative benefit 4S/(2+S) and cost 0.2S. Coefficients are chosen to show competing effects, not to identify a therapeutic target.")
    h=L(0,1)
    add(4,1,"The effective bottleneck controls drift","Initial heteroplasmy","One-generation sampling variance",h,
        [("Effective copies = 20",h*(1-h)/20),("Effective copies = 100",h*(1-h)/100)],
        "Binomial sampling variance h(1-h)/Ne is plotted for two effective copy numbers. Molecular correlations and selection are excluded.")
    add(4,2,"Selection changes the expected fraction","Initial heteroplasmy","Fraction after selection",h,
        [("Relative advantage = -0.2",h*.8/(1-.2*h)),("Neutral",h),("Relative advantage = 0.2",h*1.2/(1+.2*h))],
        "A deterministic two-type weighting uses h' = h(1+s)/(1+sh). This expected shift is distinct from sampling variance and from sequencing error.")
    n=L(2,1000,999)
    add(5,1,"Collision risk grows with pair count","Number of recorded cells","Expected colliding pairs",n,[("10,000 equiprobable labels",n*(n-1)/20000)],
        "Independent labels are drawn uniformly from 10,000 states. The plotted expectation counts colliding pairs; it is not the probability of at least one collision.")
    t=L(0,20)
    add(5,2,"Recording windows and retention compete","Time or number of serial steps","Remaining fraction",t,
        [("Unedited: rate 0.2 per time unit",np.exp(-.2*t)),("Retained: 10% loss per step",.9**t)],
        "The first curve treats time continuously; the second is a continuous interpolation of a discrete retention model. A common horizontal scale is illustrative and does not equate time units with biological divisions.")
    p=G(.0001,.2)
    ppv=lambda p,se,sp:p*se/(p*se+(1-p)*(1-sp))
    add(6,1,"Prevalence changes the meaning of a positive","Disease prevalence","Positive predictive value",p,
        [("Sensitivity 0.8; specificity 0.99",ppv(p,.8,.99)),("Sensitivity 0.8; specificity 0.999",ppv(p,.8,.999))],
        "Bayes calculation for a binary test with constant sensitivity and specificity. Population spectrum, correlated tests and follow-up outcomes are not modeled.",xscale="log")
    sp=L(.98,.9999)
    add(6,2,"False positives can dominate screening","Specificity","Expected people per 100,000 screened",sp,
        [("False positives",100000*.995*(1-sp)),("True positives: prevalence 0.005",np.full_like(sp,400))],
        "Prevalence is 0.005 and sensitivity is 0.8. Counts are mathematical expectations, so noninteger values are permitted; no screening program was simulated or evaluated.")
    p=L(0,1)
    add(7,1,"Purity changes the apparent allele fraction","Tumor purity","Variant allele fraction",p,
        [("Clonal; one mutant copy; diploid",p/2),("Half the tumor cells carry variant",p/4)],
        "Tumor and normal cells are diploid; there is one mutant copy per carrier. Copy-number changes, allelic imbalance and sequencing error are excluded.")
    k=np.arange(1,16)
    add(7,2,"Several samples can still miss a subclone","Number of independent regions sampled","Probability of missing every region",k,[("20% occupied regions",.8**k)],
        "Each sampled region independently contains the subclone with probability 0.2. Spatial dependence can make this simple miss-probability calculation optimistic.")
    stages=np.arange(1,5)
    add(8,1,"Conditional losses multiply across delivery","Delivery stage (1 to 4)","Fraction of initial administered material",stages,
        [("Stage efficiencies 0.8, 0.5, 0.1, 0.5",np.cumprod([.8,.5,.1,.5])),("Third-stage efficiency doubled",np.cumprod([.8,.5,.2,.5]))],
        "Efficiencies are conditional on reaching the previous stage; they are selected illustrative values. Final yields are 0.02 and 0.04, not measured mRNA delivery rates.")
    t=L(0,30)
    add(8,2,"RNA and protein follow different time courses","Time (arbitrary units)","Normalized molecular abundance",t,
        [("mRNA: decay rate 0.5",np.exp(-.5*t)),("Protein: production 1; decay 0.1",(np.exp(-.5*t)-np.exp(-.1*t))/(.1-.5))],
        "Initial mRNA is 1 and initial protein is 0. Linear production and first-order loss are assumed. Protein abundance uses the same arbitrary reference amount; it need not remain below 1.")
    t=L(0,30);s=.99*np.exp(-.5*t);r=.01*np.exp(.1*t)
    add(9,1,"A response can conceal a growing escape clone","Time (arbitrary units)","Normalized tumor burden",t,
        [("Sensitive cells",s),("Resistant cells",r),("Total burden",s+r)],
        "The sensitive population starts at 0.99 with net rate -0.5, and resistant cells at 0.01 with net rate +0.1. The model omits immune interactions, capacity limits and treatment changes.",yscale="log")
    add(9,2,"Composition changes before total regrowth","Time (arbitrary units)","Resistant fraction of remaining burden",t,[("Resistant / total",r/(s+r))],
        "This fraction is derived from the same two-population model as Figure 1. It is not a response curve from a clinical trial.")
    l=L(0,5)
    add(10,1,"An unchanged genome need not imply no transmission","Expected substitutions along a lineage","Probability of zero substitutions",l,[("Poisson mutation model",np.exp(-l))],
        "A Poisson count with mean lambda has zero-count probability exp(-lambda). Identical sampled genomes alone do not establish direction or a direct transmission event.")
    q=L(0,1)
    add(10,2,"Sampling gaps compound across a chain","Independent sampling probability per case","Probability all five cases are sampled",q,[("Five-case chain",q**5)],
        "All five cases must be sampled independently. Selection correlated with symptoms, place or time invalidates this simple expression.")
    n=np.arange(1,51)
    add(11,1,"Repeat visits help only under stated timing assumptions","Number of independent overpasses","Probability of at least one detection",n,[("Source active 10% of the time",1-.9**n)],
        "Each overpass independently intersects an emitting episode with probability 0.1 and detects it perfectly. Weather, detection thresholds and temporal correlation are excluded.")
    x=L(0,1)
    add(11,2,"A second view can break a source ambiguity","Source 1 contribution (normalized units)","Predicted observation (normalized units)",x,
        [("Observation 1: x1 + x2",np.ones_like(x)),("Observation 2: x1 + 0.2 x2",x+.2*(1-x))],
        "The first observation fixes x1+x2=1 but leaves the partition unidentified. An independent second sensitivity pattern changes along that ambiguity; no real transport operator is fitted.")
    dh=L(-5,5)
    add(12,1,"Storage conversion depends on specific yield","Water-table change (m)","Storage change (million cubic metres)",dh,
        [("Specific yield 0.05",100e6*.05*dh/1e6),("Specific yield 0.10",100e6*.1*dh/1e6),("Specific yield 0.20",100e6*.2*dh/1e6)],
        "Unconfined area is fixed at 100 square kilometres. Uniform water-table change and constant specific yield are assumed; confined-aquifer storage uses a different coefficient.")
    t=L(0,10);draw=np.where(t<=5,t,10-t);maxdraw=np.maximum.accumulate(draw);comp=.1*np.maximum(maxdraw-2,0)
    add(12,2,"Recovered water levels can coexist with lasting compaction","Time (arbitrary units)","Displacement (m)",t,
        [("Drawdown",draw),("Irreversible compaction",comp)],
        "A constructed drawdown-recovery cycle peaks at 5 m. Permanent compaction is 0.1 times historic drawdown beyond a 2 m threshold. This phenomenological memory rule is not a calibrated aquifer law.")
    h=L(.5,1.5)
    add(13,1,"A nonlinear flux amplifies thickness changes","Thickness / reference thickness","Flux / reference flux",h,
        [("Exponent 2",h**2),("Exponent 4",h**4),("Exponent 6",h**6)],
        "Dimensionless power laws isolate sensitivity to thickness. They are not forecasts for a named glacier and omit evolving buttressing and bed geometry.")
    add(13,2,"The same thickness can support different fluxes","Thickness / reference thickness","Normalized flux",h,
        [("Flux multiplier 0.5",.5*h**4),("Flux multiplier 0.75",.75*h**4),("Flux multiplier 1",h**4)],
        "The exponent is held at 4 while a dimensionless multiplier changes. This is a sensitivity factor, not a physically calibrated measure of ice-shelf buttressing.")
    mu=L(0,2)
    exceed=np.array([.5*math.erfc((1.2815515655446004-v)/math.sqrt(2)) for v in mu])
    add(14,1,"A fixed threshold detects a shifting mean","Mean warming / baseline standard deviation","Daily exceedance probability",mu,
        [("Fixed baseline 90th percentile",exceed),("Threshold tracks the same distribution",np.full_like(mu,.1))],
        "Daily anomalies are Gaussian with unchanged unit variance. A moving threshold follows the mean exactly. Neither curve alone specifies a multi-day marine heatwave event.")
    r=L(.1,.95)
    add(14,2,"Persistence controls runs of hot days","Probability next day stays above threshold","Probability five specified days all exceed",r,
        [("Initial exceedance probability 0.1",.1*r**4),("Independent days",np.full_like(r,1e-5))],
        "A two-state Markov construction starts with exceedance probability 0.1. Conditional persistence controls the next four days; the curve is not a count of distinct heatwave episodes.")
    lam=L(.15,1.5)
    add(15,1,"Increasing variance has more than one cause","Recovery rate (per arbitrary time unit)","Stationary variance (arbitrary units squared)",lam,
        [("Noise amplitude 1",1/(2*lam)),("Noise amplitude 2",4/(2*lam))],
        "A stationary Ornstein-Uhlenbeck process has variance sigma squared/(2 lambda). Changes in forcing amplitude can mimic a variance signal of weaker recovery.")
    add(15,2,"Autocorrelation depends on the sampling interval","Recovery rate (per arbitrary time unit)","Lag-one correlation",lam,
        [("Sampling interval 1",np.exp(-lam)),("Sampling interval 2",np.exp(-2*lam))],
        "Autocorrelation is exp(-lambda times the sampling interval). Stationarity and the one-dimensional linear model are assumptions; no AMOC tipping date is estimated.")
    kt=L(0,5)
    add(16,1,"Equilibrium capacity is not a cycle yield","Rate constant times exposure duration","Fraction of equilibrium uptake",kt,
        [("First-order equilibration",1-np.exp(-kt))],
        "The uptake fraction follows 1-exp(-kt). A short cycle can leave capacity unused even when equilibrium uptake is high; desorption and water collection are not included.")
    cyc=np.arange(1,7)
    add(16,2,"Delivered water requires collection efficiency","Completed cycles per day","Delivered water (kg per day)",cyc,
        [("Collection efficiency 0.7",1*.2*cyc*.7),("Ideal collection",1*.2*cyc)],
        "Sorbent mass is 1 kg and usable uptake swing is 0.2 kg water per kg sorbent per cycle. The same full uptake swing is assumed on every cycle; cycle count is not freely achievable.")
    dt=L(0,15)
    add(17,1,"A cold surface and a useful cooling load differ","Temperature below ambient (K)","Available cooling power (W per square metre)",dt,
        [("Net ambient radiative power 80",80-8*dt),("20 W per square metre solar penalty",60-8*dt)],
        "A linear loss coefficient of 8 W per square metre per K opposes a fixed ambient net cooling term. Negative values indicate heating in this simplified model.")
    al=L(0,.1)
    add(17,2,"Small absorptance changes have a large solar cost","Additional solar absorptance","Additional heat input (W per square metre)",al,
        [("Irradiance 500 W per square metre",500*al),("Irradiance 1,000 W per square metre",1000*al)],
        "Additional solar heat is incident irradiance times added absorptance. This excludes wavelength dependence and atmospheric-window exchange.")
    eta=L(.5,1)
    add(18,1,"First-cycle loss consumes a finite inventory","First-cycle inventory retention","Capacity (normalized units)",eta,
        [("Remaining inventory",100*eta),("Irreversibly consumed inventory",100*(1-eta))],
        "The starting sodium inventory is 100 normalized capacity units. Retention accounts only for the first cycle and does not predict later aging.")
    f=L(.2,1)
    add(18,2,"Inactive mass changes the usable metric","Cell mass / total pack mass","Pack specific energy (Wh per kg)",f,
        [("Cell specific energy 150 Wh per kg",150*f)],
        "All cells are assigned an illustrative 150 Wh per kg. Packaging reduces the gravimetric result by the cell-mass fraction; voltage, thermal design and usable state-of-charge limits are omitted.")
    z=L(0,5)
    add(19,1,"High zT approaches an ideal thermal limit","Effective average zT","Ideal conversion efficiency (%)",z,
        [("Hot 600 K; cold 300 K",50*(np.sqrt(1+z)-1)/(np.sqrt(1+z)+.5))],
        "The constant-property ideal formula gives 10.8% at zT=1 and 20.0% at zT=3. Contact losses and the power needed for heat rejection are omitted.")
    r=G(.03,30)
    add(19,2,"Load matching maximizes this power model","Load resistance / internal resistance","Power / maximum power",r,
        [("Fixed thermoelectric voltage",4*r/(1+r)**2)],
        "A fixed source voltage and internal resistance give P/Pmax=4r/(1+r)^2. Thermal feedback is absent, so this curve does not determine an efficiency optimum.",xscale="log")
    f=L(0,.1)
    add(20,1,"A small defect path can erase selectivity","Defective area fraction","Effective permeance selectivity",f,
        [("Selective path 100:1; defect path 100:100",100/((1-f)+100*f))],
        "Parallel paths contribute by area. The selective path has A/B permeances 100/1 and the nonselective path 100/100 in common arbitrary units; no membrane data are fitted.")
    thick=G(.1,10)
    add(20,2,"Thinning raises ideal permeance, not selectivity","Selective-layer thickness / reference","Permeance / reference",thick,[("Constant permeability",1/thick)],
        "With constant permeability and no support resistance or defects, permeance is inversely proportional to thickness. Mechanical stability and process recovery are outside the model.",xscale="log",yscale="log")

    tau=G(1,10000)
    add(21,1,"Averaging suppresses white noise under an ideal model","Averaging time / reference time","Relative instability / reference",tau,
        [("White-frequency noise",1/np.sqrt(tau))],
        "The curve is normalized tau to the power -1/2. It does not include a systematic floor or establish a clock's absolute accuracy.",xscale="log",yscale="log")
    rho=L(-.9,.99)
    add(21,2,"Correlations determine comparison uncertainty","Correlation of two equal-uncertainty errors","Uncertainty of difference / single-clock uncertainty",rho,
        [("Equal individual uncertainties",np.sqrt(2-2*rho))],
        "For two estimates with uncertainty u and correlation rho, the difference has uncertainty u times sqrt(2-2rho). Shared errors can cancel in a difference while remaining in absolute values.")
    t=L(0,3)
    add(22,1,"Longer interrogation amplifies the ideal signal","Interrogation time / reference time","Relative acceleration scale factor",t,[("Fixed effective wave vector",t**2)],
        "The ideal scale factor is proportional to T squared. The normalized curve assumes fixed effective wave vector and neglects pulse-duration and trajectory corrections.")
    add(22,2,"Lost contrast can reverse the time advantage","Interrogation time / contrast-decay time","Normalized information-amplitude proxy",t,
        [("Contrast times scale factor",t**2*np.exp(-t*t))],
        "Contrast is assigned exp[-(T/Tc)^2]. The plotted amplitude proxy T squared times contrast peaks at T/Tc=1; atom loss, cycle rate and detection noise are not modeled.")
    u=L(8,24)
    add(23,1,"The spin scale falls as interaction grows","Interaction U / tunneling t","Exchange J / tunneling t",u,[("Strong-coupling approximation",4/u)],
        "Second-order strong-coupling theory gives J=4t squared/U. The range U/t from 8 to 24 emphasizes the asymptotic regime; higher-order corrections are omitted.")
    e=L(0,.3)
    add(23,2,"Independent readout errors suppress correlations","Probability of flipping each spin label","Observed / true two-spin correlation",e,
        [("Independent symmetric label flips",(1-2*e)**2)],
        "Each of the two binary labels is independently flipped with probability e. Correlated errors and detection selection require a different correction.")
    r=G(.001,100)
    add(24,1,"A conducting path can dominate a parallel measurement","Resistance of path 1 (ohm)","Measured parallel resistance (ohm)",r,
        [("Path 2 fixed at 100 ohm",100*r/(100+r))],
        "A two-resistor network illustrates current-path sensitivity. It does not assign a superconducting volume fraction or represent a measured nickelate sample.",xscale="log",yscale="log")
    x=L(0,1)
    add(24,2,"An uncontrolled variable can reshape an apparent dome","Composition coordinate (normalized)","Transition-scale proxy / reference",x,
        [(f"Illustrative width {w}",np.maximum(0,1-((x-.5)/w)**2)) for w in [.1,.2,.3]],
        "Three clipped parabolas share center 0.5 and peak 1 but differ in width. These are schematic functions, not transition-temperature data or a phase diagram.")
    w=L(0,4)
    def lor(g):
        return sum(g/math.pi/((w-c)**2+g*g) for c in [1,2,3])
    add(25,1,"Broadening can turn discrete peaks into a continuum","Excitation energy (arbitrary units)","Spectral density (inverse energy units)",w,
        [("Lorentzian width 0.1",lor(.1)),("Lorentzian width 0.6",lor(.6))],
        "Three unit-area Lorentzian components are centered at 1, 2 and 3. Only the width changes. A broad spectrum alone cannot identify fractional excitations.")
    bg=L(-.3,.3)
    add(25,2,"A reference value is not a channel assignment","Added background / thermal conductance quantum","Total / thermal conductance quantum",bg,
        [("Reference 0.5 plus background",.5+bg),("Reference only",np.full_like(bg,.5))],
        "A reference of one-half quantum is combined with an arbitrary additive background. The calculation is an identifiability example, not evidence of a quantized or topological channel.")
    mu=L(2,30);kb=1.380649e-23;mh=1.6735575e-27
    add(26,1,"Composition and gravity enter a spectral scale","Mean molecular weight (hydrogen-mass units)","Atmospheric scale height (km)",mu,
        [("Temperature 1,000 K; gravity 10 m/s squared",kb*1000/(mu*mh*10)/1000)],
        "An isothermal ideal-gas scale height is calculated with constant gravity. Cloud altitude, reference radius and opacity still influence an actual transmission spectrum.")
    wave=L(1,5)
    add(26,2,"A small stellar contrast can mimic a planetary feature","Wavelength (micrometres)","Apparent transit depth (%)",wave,
        [("Featureless planet; varying stellar factor",100*.01*(1+.02*np.sin(2*np.pi*(wave-1)/4))),("Planet alone",np.full_like(wave,1.0))],
        "The planet's true depth is 0.01 and a purely constructed stellar multiplier varies sinusoidally by 2%. The sinusoid is not a stellar-atmosphere model or observed spectrum.")
    dm=L(-.3,.3)
    add(27,1,"A distance calibration shift propagates to H0","Distance-modulus shift (mag)","Value / original value",dm,
        [("Inferred distance",10**(dm/5)),("Inferred H0 at fixed redshift",10**(-dm/5))],
        "A modulus offset changes distance by 10^(offset/5). At fixed low redshift the inferred expansion rate changes inversely; the plot does not resolve the observed Hubble tension.")
    rho=L(-.8,.9)
    add(27,2,"Tension depends on the full covariance","Correlation between two estimates","Difference / standard uncertainty",rho,
        [("Difference 6; individual uncertainties 1",6/np.sqrt(2-2*rho))],
        "Both illustrative standard uncertainties equal 1 in the same units, and the mean difference is 6. This is a covariance sensitivity calculation, not a significance estimate for published data.")
    expo=L(1,10)
    add(28,1,"A background-free limit scales with exposure","Exposure / reference exposure","90% upper rate limit (per reference exposure)",expo,
        [("Zero observed events; known zero background",-math.log(.1)/expo)],
        "The one-sided Poisson upper mean for zero counts is -ln(0.1). Known zero background and fixed efficiency are required; real searches use additional information and nuisance models.")
    eth=L(0,8)
    add(28,2,"Threshold choices select different signal fractions","Recoil threshold (keV)","Surviving signal fraction",eth,
        [("Illustrative exponential scale 2 keV",np.exp(-eth/2))],
        "The signal spectrum is assigned a simple exponential with scale 2 keV and unit efficiency above threshold. No dark-matter mass, cross section or exclusion contour is inferred.")
    z=L(0,3)
    add(29,1,"Host dispersion is diluted by redshift","Source redshift","Observed host dispersion (pc per cubic centimetre)",z,
        [("Rest-frame host contribution 100",100/(1+z))],
        "A fixed rest-frame host contribution of 100 pc per cubic centimetre is divided by 1+z. Real hosts differ, and selection can correlate their contribution with redshift.")
    n=L(1,400,400)
    add(29,2,"Independent scatter averages slowly","Number of independent bursts","Standard error (pc per cubic centimetre)",n,
        [("Per-burst scatter 200",200/np.sqrt(n))],
        "A common standard deviation of 200 and independent sampling produce the square-root law. Shared foreground errors and selection biases do not disappear this way.")
    r=L(0,.99)
    add(30,1,"Foreground flexibility can absorb a signal","Absolute correlation with a nuisance template","Fraction of signal norm squared retained",r,
        [("Orthogonal residual after one-template fit",1-r*r)],
        "For two normalized templates in a Euclidean inner product, projection removes the squared correlation. An experiment's actual noise-weighted templates must be evaluated separately.")
    frac=G(1e-6,1e-3)
    add(30,2,"A tiny fractional residual can be astrophysically large","Fractional foreground calibration error","Residual temperature (K)",frac,
        [("Foreground temperature 1,000 K",1000*frac)],
        "A fixed 1,000 K foreground multiplied by the fractional error gives residual amplitude. Spectral structure and coupling to the antenna are excluded.",xscale="log",yscale="log")
    r=G(.1,10)
    add(31,1,"Selection changes the apparent state population","State-2 / state-1 selection probability","Observed fraction of state 2",r,
        [("True populations are equal",r/(1+r))],
        "Two conformations have equal true populations. Only their relative selection probability changes, illustrating a bias from picking or reconstruction.",xscale="log")
    add(31,2,"A population bias becomes an apparent energy gap","State-2 / state-1 selection probability","Apparent free-energy difference / kBT",r,
        [("Equal true free energies",-np.log(r))],
        "Applying an equilibrium Boltzmann relation to selected fractions creates an apparent gap -ln(r). This does not justify equilibrium free-energy inference from non-equilibrium or selected particles.",xscale="log")
    comp=np.arange(1,5)
    add(32,1,"Identical acquired data can hide distinct images","Image component index","Component amplitude (arbitrary units)",comp,
        [("Candidate x",np.array([1,2,0,0])),("Candidate x plus null-space component",np.array([1,2,1,-1]))],
        "The illustrative acquisition retains only components 1 and 2. Both candidates therefore produce the same measurements although components 3 and 4 differ. Lines only connect component markers.")
    n=G(100,1000000)
    add(32,2,"A small global error can hide a missing feature","Number of image pixels","Mean squared error",n,
        [("One unit-amplitude pixel is deleted",1/n)],
        "All other pixels are exact and a single pixel of amplitude 1 is replaced by zero. This deliberately simple construction separates average error from preservation of a small feature.",xscale="log",yscale="log")
    p=G(.001,.5)
    add(33,1,"A transported score meets a new prevalence","Condition prevalence","Positive predictive value",p,
        [("Sensitivity = specificity = 0.95",ppv(p,.95,.95))],
        "Sensitivity and specificity are held fixed to isolate prevalence. Real site shifts may change both, making recalibration alone insufficient.",xscale="log")
    m=L(1,200,200)
    add(33,2,"Many tiles are not many independent patients","Tiles per patient","Approximate effective tile-sample size",m,
        [("100 patients; within-patient correlation 0.2",100*m/(1+(m-1)*.2)),("100 patients; correlation 0.8",100*m/(1+(m-1)*.8))],
        "The design-effect approximation assumes equal cluster sizes and exchangeable correlations. This effective size does not license splitting a patient's tiles between training and testing.")
    x=L(0,2*np.pi,1501)
    add(34,1,"A small energy error can oscillate rapidly","Coordinate (arbitrary length units)","Energy perturbation (arbitrary energy units)",x,
        [("Frequency 1",.01*np.sin(x)),("Frequency 10",.01*np.sin(10*x))],
        "Both displayed energy perturbations have amplitude 0.01. Frequencies are expressed per coordinate unit; these are analytic functions, not errors measured from a trained potential.")
    k=G(1,100)
    add(34,2,"Equal energy error can hide unequal force error","Spatial frequency (per arbitrary length unit)","Peak force-error amplitude (energy per length unit)",k,
        [("Energy-error amplitude fixed at 0.01",.01*k)],
        "Differentiating an energy perturbation 0.01 sin(kx) gives force-error amplitude 0.01k. Increasing frequency from 1 to 100 raises the peak force error from 0.01 to 1 without changing peak energy error.",xscale="log",yscale="log")
    p=G(.001,.5)
    post=lambda p,lr:lr*p/(1-p+lr*p)
    add(35,1,"The same likelihood ratio implies different posteriors","Prior probability","Posterior probability",p,
        [("Likelihood ratio 10",post(p,10)),("No update",p)],
        "A calibrated likelihood ratio of 10 is assumed for illustration. A variant-effect score is not automatically a clinical likelihood ratio.",xscale="log")
    add(35,2,"Correlated predictions cannot be counted twice","Prior probability","Posterior probability",p,
        [("One valid likelihood ratio of 10",post(p,10)),("Incorrect duplicate multiplication: 100",post(p,100))],
        "The second curve deliberately demonstrates double counting if two perfectly redundant reports are treated as independent. It is not a recommended evidence-combination rule.",xscale="log")
    l=L(0,15)
    add(36,1,"Mean coverage leaves a dropout tail","Mean independent reads per strand","Probability of zero reads",l,
        [("Poisson read sampling",np.exp(-l))],
        "A Poisson model gives dropout exp(-lambda). Uneven strand abundance, synthesis bias and correlated errors can produce a heavier tail.",yscale="log")
    n=G(1,10000000)
    add(36,2,"Small strand failure rates compound across an archive","Number of required strands","Probability every strand is recovered",n,
        [("Independent residual failure probability 0.000001",np.exp(n*np.log1p(-1e-6)))],
        "The calculation assumes every strand is required and fails independently after decoding. Erasure codes and correlated failures change the archive-level recovery model.",xscale="log")
    r=G(1,1000)
    add(37,1,"Compute-only savings meet a system-level ceiling","Improvement factor of accelerated operations","End-to-end energy improvement factor",r,
        [("Accelerated share of original energy 0.8",1/(.2+.8/r)),("Accelerated share 0.2",1/(.8+.2/r))],
        "Amdahl-style energy accounting assumes all unaccelerated energy stays fixed. Converters, data movement and control must be included when estimating the accelerated share.",xscale="log")
    d=G(.001,1)
    add(37,2,"Quantization noise follows a declared error model","Quantization step (normalized units)","Quantization-error variance",d,
        [("Uniform independent error within one step",d*d/12)],
        "Variance is step squared/12 for a uniform error on minus half a step to plus half a step. Clipping, correlated error and device drift are excluded.",xscale="log",yscale="log")
    extra=L(0,2)
    add(38,1,"An active perturbation separates identical contacts","Known additional force (N)","Contact deformation (m)",extra,
        [("Initial force 1 N; stiffness 10 N/m",(1+extra)/10),("Initial force 2 N; stiffness 20 N/m",(2+extra)/20)],
        "Two linear contacts have the same initial deformation of 0.1 m. A known force increment produces different slopes, resolving this idealized force-stiffness ambiguity.")
    mu=L(.1,1)
    add(38,2,"Grip requirements depend on uncertain friction","Coulomb friction coefficient","Minimum normal force (N)",mu,
        [("Tangential load 2 N",2/mu)],
        "The static inequality uses tangential load 2 N and one resultant normal force. Dynamics, compliance, multiple fingers and object damage limits are outside the model.")
    dt=G(.005,.5)
    add(39,1,"Finer sampling can amplify derivative noise","Sample interval (arbitrary time units)","Noise standard deviation of forward derivative",dt,
        [("Independent observation noise SD 0.01",np.sqrt(2)*.01/dt)],
        "A forward difference of two independent noisy observations has variance 2 sigma squared/dt squared. Correlated noise and smoothed derivatives require different formulas.",xscale="log",yscale="log")
    r=L(0,.995)
    add(39,2,"Correlated candidate terms destabilize coefficients","Absolute correlation between two unit-norm columns","Condition number of the Gram matrix",r,
        [("Two-column normalized library",(1+r)/(1-r))],
        "The Gram matrix has diagonal entries 1 and off-diagonal correlation r. This is the Gram condition number, not the square-root condition number of the design matrix.",yscale="log")
    temp=L(1,600)
    add(40,1,"The ideal erasure scale grows with temperature","Reservoir temperature (K)","Minimum mean erasure work (zeptojoules)",temp,
        [("One unbiased bit; reversible limit",kb*temp*math.log(2)/1e-21)],
        "The reversible isothermal bound kBT ln(2) applies to resetting one unbiased bit with a specified logical and thermodynamic cycle. It is not an energy requirement for every logic gate.")
    eps=L(0,.5)
    h=np.zeros_like(eps);mask=eps>0
    h[mask]=-eps[mask]*np.log(eps[mask])-(1-eps[mask])*np.log1p(-eps[mask])
    add(40,2,"Permitted error changes the information erased","Final probability of the undesired bit state","Work bound / full-erasure bound",eps,
        [("Symmetric memory; same-state endpoint energies",(math.log(2)-h)/math.log(2))],
        "For an initially unbiased bit and equal logical-state endpoint energies, reducing entropy to binary entropy H(epsilon) costs at least kBT[ln(2)-H]. Finite-time dissipation and control overhead are omitted.")

def render():
    out=ROOT/"figures";out.mkdir(exist_ok=True)
    plt.rcParams.update({"font.family":"DejaVu Sans","font.size":10,"axes.labelsize":10,
        "axes.titlesize":14,"axes.titleweight":"bold","axes.spines.top":False,"axes.spines.right":False,
        "axes.edgecolor":"#8d9994","axes.labelcolor":"#233b3a","xtick.color":"#4f615d",
        "ytick.color":"#4f615d","text.color":"#233b3a","figure.facecolor":"#faf9f4",
        "axes.facecolor":"#faf9f4","savefig.facecolor":"#faf9f4","svg.hashsalt":"auturnal-english-frontiers-40-v1"})
    colors=["#205e67","#ad683c","#657fa3","#85566f"]
    for p in PANELS:
        fig,ax=plt.subplots(figsize=(7.2,4.8),layout="constrained")
        fig.set_constrained_layout_pads(w_pad=.12,h_pad=.14)
        for k,s in enumerate(p["series"]):
            x=np.asarray(s["x"]);y=np.asarray(s["y"])
            assert len(x)==len(y) and np.isfinite(x).all() and np.isfinite(y).all()
            ax.plot(x,y,color=colors[k%4],lw=2.4,label=s["label"],linestyle=["-","--","-.",":"][k%4],
                marker="o" if len(x)<=15 else None,markersize=4)
        ax.set(xlabel=p["xlabel"],ylabel=p["ylabel"],xscale=p["xscale"],yscale=p["yscale"])
        ax.set_title("\n".join(textwrap.wrap(p["title"],48)),loc="left",pad=15)
        ax.grid(alpha=.2,color="#82938d",lw=.6)
        handles,labels=ax.get_legend_handles_labels()
        ax.legend(handles,["\n".join(textwrap.wrap(v,46)) for v in labels],loc="best",frameon=True,
            facecolor="#faf9f4",edgecolor="#e0e4de",fontsize=8.5)
        stem=p["slug"]+f"-{p['number']:02}"
        fig.savefig(out/(stem+".png"),dpi=180,metadata={"Software":"Auturnal analytical figures; constructed examples"})
        fig.savefig(out/(stem+".svg"),metadata={"Date":None,"Creator":"Auturnal reproducible analytical figures"})
        plt.close(fig)
        with (out/(stem+".csv")).open("w",encoding="utf-8",newline="") as f:
            w=csv.writer(f);w.writerow(["series","x","y"])
            for s in p["series"]:
                for x,y in zip(s["x"],s["y"]):w.writerow([s["label"],format(x,".15g"),format(y,".15g")])
        p["image"]="figures/"+stem+".png";p["vector"]="figures/"+stem+".svg";p["data"]="figures/"+stem+".csv"
        p["alt"]=p["title"]+". "+p["caption"]
    (out/"panels.json").write_text(json.dumps(PANELS,indent=2,ensure_ascii=False)+"\n",encoding="utf-8")
    print(json.dumps({"figures":len(PANELS),"articles":len(set(p["slug"] for p in PANELS))}))

if __name__=="__main__":
    TOPICS=json.loads((ROOT/"topics.json").read_text(encoding="utf-8"))
    models()
    assert len(PANELS)==80
    render()