legacy/code/run_experiments.py

"""Reproduce every reported AUREOLE v1 numerical experiment, without downloads.

Run from the package root: python code/run_experiments.py
No GPU, trained network, path tracer, DLSS SDK, or internet access is used.
"""
from pathlib import Path
import csv, json, platform, sys, time
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from scipy.stats import t as student_t
from aureole_core import (observe, query_value, future_metric, risk,
    transform_coding, sqrt_psd, batch_covariance, closed_observation_basis,
    area_mixture)

ROOT = Path(__file__).resolve().parents[1]
OUT, FIG = ROOT / "results", ROOT / "figures"
OUT.mkdir(exist_ok=True); FIG.mkdir(exist_ok=True)
plt.rcParams.update({"font.family":"DejaVu Sans", "font.size":10,
    "axes.spines.top":False,"axes.spines.right":False,"figure.dpi":140})
COLORS = {"frame":"#979BA8", "screen":"#E2A34A", "world":"#116E83", "versioned":"#614DA1"}


def savefig(name):
    plt.savefig(FIG / (name + ".png"), bbox_inches="tight")
    plt.savefig(FIG / (name + ".pdf"), bbox_inches="tight")
    plt.close()


def stats(x):
    x = np.asarray(x, float)
    mean = float(x.mean())
    se = float(x.std(ddof=1) / np.sqrt(len(x))) if len(x)>1 else 0.
    half = float(student_t.ppf(.975, len(x)-1)*se) if len(x)>1 else 0.
    return {"n":len(x), "mean":mean, "ci95":[mean-half,mean+half], "std":float(x.std(ddof=1)) if len(x)>1 else 0.}


def csv_out(name, rows):
    with (OUT / name).open("w",newline="") as f:
        writer=csv.DictWriter(f,fieldnames=list(rows[0])); writer.writeheader(); writer.writerows(rows)


def run_atlas(changed, nseeds=24):
    """Controlled 2.5D surface atlas; synthetic Gaussian noise, known identities.

    64x160 canonical cells, 64x48 orthographic viewport. A revisited region is
    absent for frames 24..523 (exactly 500 frames). Observations and noise are
    identical across baselines. Lighting is known and changes on return.
    screen retains only cells visible in the immediately preceding frame, with
    perfect reprojection. world retains all seen cells. versioned additionally
    receives an authoritative material-generation change event on return.
    """
    height,width,view=64,160,48
    yy,xx=np.mgrid[:height,:width]
    base=.48+.16*np.sin(xx*.23)*np.cos(yy*.19)+.16*((xx//5+yy//7)%2-.5)
    base=np.clip(base,.08,.92)
    new=base.copy(); new[:,:80]=1-base[:,:80]
    prior_mean,prior_var,noise=.5,.09,.0144
    p_sample=.125
    origins=[(t//3)%9 for t in range(24)]+[104+(t%8) for t in range(500)]+[(t//3)%9 for t in range(24)]
    methods=["frame","screen","world","versioned"]
    raw=[]; trace={m:[] for m in methods}; seed_images={}
    for seed in range(nseeds):
        rng=np.random.default_rng(1729+seed)
        means={m:np.full((height,width),prior_mean) for m in methods}
        variances={m:np.full((height,width),prior_var) for m in methods}
        previous=np.zeros((height,width),bool)
        local={m:[] for m in methods}
        for t,origin in enumerate(origins):
            vis=np.zeros((height,width),bool); vis[:,origin:origin+view]=True
            truth=new if changed and t>=524 else base
            lighting=.72 if t<524 else 1.18
            mask=(rng.random((height,view))<p_sample)
            target=lighting*truth[:,origin:origin+view]
            samples=target+rng.normal(0,np.sqrt(noise),target.shape)
            for method in methods:
                if method=="frame":
                    means[method].fill(prior_mean); variances[method].fill(prior_var)
                if method=="screen":
                    fresh=vis & ~previous
                    means[method][fresh]=prior_mean; variances[method][fresh]=prior_var
                if method=="versioned" and changed and t==524:
                    means[method][:,:80]=prior_mean; variances[method][:,:80]=prior_var
                m=means[method][:,origin:origin+view]
                p=variances[method][:,origin:origin+view]
                k=p*lighting/(noise+lighting*lighting*p)
                m[mask]+=k[mask]*(samples[mask]-lighting*m[mask])
                p[mask]*=(1-k[mask]*lighting)
                mse=float(np.mean((lighting*m-target)**2))
                if t>=524:
                    local[method].append(mse)
                    raw.append({"changed":changed,"seed":seed,"return_frame":t-524,"method":method,"mse":mse})
                if seed==0 and t==524:
                    seed_images[method]=lighting*m.copy()
            if seed==0 and t==524: seed_images["reference"]=target.copy()
            previous=vis
        for m in methods: trace[m].append(local[m])
    csv_out("atlas_"+("changed" if changed else "static")+".csv",raw)
    summary={m:{"first_return":stats(np.array(trace[m])[:,0]),"return_24_mean":stats(np.array(trace[m]).mean(1))} for m in methods}
    fig,ax=plt.subplots(figsize=(6.6,3.6))
    for m in methods:
        a=np.array(trace[m]); ax.plot(a.mean(0),label=m,color=COLORS[m],lw=2)
    ax.set(xlabel="Frames since return after 500 absent frames",ylabel="Mean squared radiance error",title="Unannounced to stale cache; signaled to versioned cache" if changed else "Static material; changed known illumination")
    ax.legend(ncol=2,frameon=False); savefig("atlas_"+("changed" if changed else "static"))
    if not changed:
        fig,axes=plt.subplots(1,4,figsize=(8.8,3.2))
        for ax,name in zip(axes,["reference","frame","screen","world"]):
            ax.imshow(seed_images[name],cmap="gray",vmin=0,vmax=1.2); ax.set_title(name); ax.axis("off")
        fig.suptitle("Controlled atlas: first revisit, seed 1729 (not a path-traced game)")
        savefig("atlas_views")
    return summary


def active_allocation(nseeds=2000):
    """Exact matched-prior Bayes experiment. Equal scalar sample costs."""
    d,budget=12,36
    prior=np.diag(np.linspace(.6,1.4,d))
    # Current image uses first mode; a declared future output/task set uses 4.
    future=np.diag([24.,12.,6.,3.,.15,.15,.15,.15,.15,.15,.15,.15])
    immediate=np.diag([1.]+[0.]*(d-1))
    hs=np.eye(d); noise=.25
    names=["uniform","entropy","current_only","future_loss"]
    curves={}; schedules={}; final={}
    for name in names:
        p=prior.copy(); curve=[risk(p,future)]; schedule=[]
        for step in range(budget):
            if name=="uniform": j=step%d
            elif name=="entropy": j=int(np.argmax(np.diag(p)))
            else:
                w=immediate if name=="current_only" else future
                j=int(np.argmax([query_value(p,w,h,noise) for h in hs]))
            schedule.append(j); _,p=observe(np.zeros(d),p,hs[j],0.,noise)
            curve.append(risk(p,future))
        curves[name]=curve; schedules[name]=schedule; final[name]=p
    rng=np.random.default_rng(8431)
    truths=rng.multivariate_normal(np.zeros(d),prior,size=nseeds)
    common_noise=rng.normal(0,np.sqrt(noise),(nseeds,budget,d))
    summary={}; raw=[]
    for name in names:
        m=np.zeros((nseeds,d)); p=prior.copy()
        for step,j in enumerate(schedules[name]):
            k=p[:,j]/(noise+p[j,j]); observation=truths[:,j]+common_noise[:,step,j]
            m+=(observation-m[:,j])[:,None]*k[None,:]
            _,p=observe(np.zeros(d),p,hs[j],0.,noise)
        err=truths-m; loss=np.einsum("ni,ij,nj->n",err,future,err)
        summary[name]={"expected_loss":risk(final[name],future),"monte_carlo":stats(loss),"allocation":np.bincount(schedules[name],minlength=d).tolist()}
        raw.extend({"seed_draw":j,"method":name,"loss":float(v)} for j,v in enumerate(loss))
    csv_out("active_samples.csv",raw)
    csv_out("active_curves.csv",[{"queries":i,**{n:curves[n][i] for n in names}} for i in range(budget+1)])
    plt.figure(figsize=(6.6,3.6))
    for n in names: plt.plot(curves[n],label=n.replace("_"," "),lw=2)
    plt.yscale("log"); plt.xlabel("Equal-cost scalar measurements"); plt.ylabel("Expected weighted future error"); plt.legend(frameon=False,ncol=2)
    plt.title("Known linear model and declared future query distribution")
    savefig("active_sampling")
    return summary


def compression_experiment():
    rng=np.random.default_rng(923)
    d=10; b=rng.normal(size=(d,d)); source=b@b.T/d+.1*np.eye(d)
    q,_=np.linalg.qr(rng.normal(size=(d,d))); w=q@np.diag(np.geomspace(30,.01,d))@q.T
    root=sqrt_psd(source); errors=[]
    gaussian=rng.normal(size=(100000,d))
    summary={}
    for rank in range(d+1):
        u,vals,cov=transform_coding(source,w,rank)
        e=gaussian@(np.eye(d)-u@u.T)@root
        realized=np.einsum("ni,ij,nj->n",e,w,e).mean()
        random=[]
        for _ in range(100):
            v,_=np.linalg.qr(rng.normal(size=(d,d)))
            random.append(risk(root@(np.eye(d)-v[:,:rank]@v[:,:rank].T)@root,w))
        errors.append({"rank":rank,"optimal":risk(cov,w),"eigen_tail":float(vals[rank:].sum()),"monte_carlo":float(realized),"random_mean":float(np.mean(random))})
    csv_out("compression.csv",errors)
    summary={"rank4":errors[4],"note":"Accessible Gaussian source innovation; not a guarantee about compressing an unknown scene."}
    plt.figure(figsize=(6.6,3.6))
    plt.plot([x["rank"] for x in errors],[x["optimal"] for x in errors],label="Task-weighted transform",lw=2)
    plt.plot([x["rank"] for x in errors],[x["random_mean"] for x in errors],label="Random rank-matched transform",lw=2)
    plt.xlabel("Retained exact linear coordinates");plt.ylabel("Weighted distortion");plt.legend(frameon=False);savefig("compression")
    return summary


def limits_experiment():
    p=np.eye(2); w=np.diag([1.,0.]); h1=np.array([1.,1.]); h2=np.array([0.,1.]); r=.1
    _,p1=observe(np.zeros(2),p,h1,0.,r)
    # Same task marginal, different nuisance uncertainty, therefore different update.
    _,known=observe(np.zeros(2),np.diag([1.,0.]),h1,0.,r)
    _,unknown=observe(np.zeros(2),p,h1,0.,r)
    # Chronoscopic teacher reveals an independent sign. Student has no evidence.
    rng=np.random.default_rng(42); bits=rng.choice([-1.,1.],size=200000)
    full=closed_observation_basis([np.eye(2)],np.vstack([np.array([1.,0.]),h1]))
    return {"query2_gain_alone":query_value(p,w,h2,r),"query2_gain_after_query1":query_value(p1,w,h2,r),
      "task_variance_if_nuisance_known":float(known[0,0]),"task_variance_if_nuisance_discarded":float(unknown[0,0]),
      "image_only_rank":1,"query_closed_rank":int(full.shape[1]),
      "future_teacher_bit_mse":0.,"causal_optimal_bit_mse":float(np.mean(bits**2)),
      "same_ray_naively_counted_20_times_variance":1/(1+20/r),"correct_one_ray_variance":1/(1+1/r)}


def weak_signals_experiment():
    # A correlated auxiliary observation must be conditioned on the existing signal.
    p=np.array([[1.,.45],[.45,1.2]]); w=np.diag([1.,2.])
    h=np.eye(2); noise=np.array([[.3,.24],[.24,.5]])
    joint=batch_covariance(p,h,noise)
    naive=batch_covariance(p,h,np.diag(np.diag(noise)))
    rng=np.random.default_rng(94); n=150000
    x=rng.multivariate_normal([0,0],p,n); e=rng.multivariate_normal([0,0],noise,n); y=x+e
    k_correct=p@np.linalg.inv(p+noise); k_naive=p@np.linalg.inv(p+np.diag(np.diag(noise)))
    ec=x-y@k_correct.T; en=x-y@k_naive.T
    return {"correct_predicted_risk":risk(joint,w),"correct_measured_risk":float(np.einsum('ni,ij,nj->n',ec,w,ec).mean()),
      "independent_assumption_predicted_risk":risk(naive,w),"independent_assumption_measured_risk":float(np.einsum('ni,ij,nj->n',en,w,en).mean())}


def main():
    started=time.perf_counter()
    report={"release":"1.0.0","date":"2026-09-19","scope":"CPU synthetic proof-of-mechanism, no trained neural renderer", "experiments":{}}
    report["experiments"]["E1_static_revisit"]=run_atlas(False)
    report["experiments"]["E2_hidden_material_change"]=run_atlas(True)
    report["experiments"]["E3_active_sampling"]=active_allocation()
    report["experiments"]["E4_transform_coding"]=compression_experiment()
    report["experiments"]["E5_counterexamples"]=limits_experiment()
    report["experiments"]["E6_correlated_weak_signals"]=weak_signals_experiment()
    report["environment"]={"python":sys.version,"numpy":np.__version__,"platform":platform.platform(),"elapsed_seconds":time.perf_counter()-started}
    (OUT/"experiment_report.json").write_text(json.dumps(report,indent=2))
    print(json.dumps(report,indent=2))


if __name__=="__main__": main()