"""Sequential measurements in a defined neutral-gel virtual laboratory. Run: python example.py Requires NumPy, SciPy, Matplotlib; all data are generated by explicit equations. The physical oracle generates labels. Selection sees queried labels only. """ from pathlib import Path import json, os os.environ.setdefault('OPENBLAS_NUM_THREADS','1') os.environ.setdefault('OMP_NUM_THREADS','1') import numpy as np from scipy.optimize import brentq from scipy.linalg import cho_factor, cho_solve, solve_triangular from scipy.special import logsumexp import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib.lines import Line2D OUT=Path('active_outputs');OUT.mkdir(exist_ok=True) NOISE=.04;AMP=1.5;PRIOR=3.5;JITTER=1e-10 LENGTHS=np.array([[.6,.15],[.6,.35],[.6,.7]]) SEEDS=list(range(1900,1912));BUDGETS=[8,12,16,20,24,28,32,36] STRATEGIES=['random','maximin','uncertainty'] axis=np.linspace(0,1,25) POOL=np.array([(a,b) for b in axis for a in axis]) def physical(x): """Return ln(volume swelling), polymer fraction and shear scale in kPa.""" x=np.atleast_2d(x);nu=2*np.exp(np.log(30)*x[:,0]);chi=.2+.35*x[:,1] phi=[] for n,c in zip(nu,chi): def equation(p):return np.log1p(-p)+p+c*p*p+n*18e-6*(p**(1/3)-p/2) phi.append(brentq(equation,1e-7,1-1e-10,xtol=1e-14)) phi=np.array(phi) return -np.log(phi),phi,nu*8.314462618*298*phi**(1/3)/1000 def kernel(a,b,length): d=(a[:,None,:]-b[None,:,:])/length return AMP**2*np.exp(-.5*np.sum(d*d,axis=2)) class MixtureGP: """Exact conditioning in each kernel, followed by Bayesian model averaging.""" def __init__(self,x,y): self.x=x.copy();self.y=y.copy();self.components=[];evidence=[] for length in LENGTHS: K=kernel(x,x,length)+(NOISE**2+JITTER)*np.eye(len(x)) factor=cho_factor(K,lower=True);alpha=cho_solve(factor,y-PRIOR) ll=-.5*((y-PRIOR)@alpha)-np.log(np.diag(factor[0])).sum()-.5*len(x)*np.log(2*np.pi) evidence.append(ll);self.components.append((length,factor,alpha)) self.log_evidence=np.array(evidence) self.weights=np.exp(self.log_evidence-logsumexp(self.log_evidence)) def predict(self,q): means=[];variances=[] for length,factor,alpha in self.components: cross=kernel(self.x,q,length) projected=solve_triangular(factor[0],cross,lower=True) means.append(PRIOR+cross.T@alpha) variances.append(np.maximum(0,AMP**2-np.sum(projected**2,axis=0))) means=np.array(means);variances=np.array(variances) mean=self.weights@means variance=self.weights@(variances+means**2)-mean**2 return mean,np.maximum(variance,0) def assess(model,x,truth): mean,var=model.predict(x) return float(np.sqrt(np.mean((mean-truth)**2))) def campaign(seed): rng=np.random.default_rng(seed) initial=rng.choice(len(POOL),8,replace=False).tolist() # Potential outcomes couple comparisons; the selector receives only queried values. pool_truth,_,_=physical(POOL);pool_y=pool_truth+rng.normal(0,NOISE,len(POOL)) cal_x=rng.random((40,2));test_x=rng.random((256,2));shift_x=rng.random((128,2)) shift_x[:,1]=1+shift_x[:,1]*(.1/.35) cal_truth=physical(cal_x)[0];test_truth=physical(test_x)[0];shift_truth=physical(shift_x)[0] cal_y=cal_truth+rng.normal(0,NOISE,40) test_y=test_truth+rng.normal(0,NOISE,256);shift_y=shift_truth+rng.normal(0,NOISE,128) random_order=rng.permutation(len(POOL)).tolist() runs={} for strategy in STRATEGIES: chosen=initial.copy();curve=[];acquired=[];weight_history=[] for n in range(8,37): model=MixtureGP(POOL[chosen],pool_y[chosen]) weight_history.append(model.weights.tolist()) if n in BUDGETS:curve.append(assess(model,test_x,test_truth)) if n==36:break if strategy=='random':idx=next(i for i in random_order if i not in chosen);score=None else: if strategy=='maximin':scores=np.min(np.sum((POOL[:,None,:]-POOL[chosen][None,:,:])**2,axis=2),axis=1) else:scores=model.predict(POOL)[1] scores[chosen]=-np.inf;idx=int(np.argmax(scores));score=float(scores[idx]) chosen.append(idx);acquired.append(dict(index=idx,score=score)) cm,cv=model.predict(cal_x);scores=abs(cal_y-cm)/np.sqrt(cv+NOISE**2) rank=int(np.ceil((len(scores)+1)*.9));qhat=float(np.sort(scores)[rank-1]) tests={} for key,xx,tt,yy in [('nominal',test_x,test_truth,test_y),('shifted',shift_x,shift_truth,shift_y)]: mean,var=model.predict(xx);half=qhat*np.sqrt(var+NOISE**2) tests[key]=dict(RMSE_logQ=float(np.sqrt(np.mean((mean-tt)**2))), coverage=float(np.mean(abs(yy-mean)<=half)),mean_interval_width_logQ=float(np.mean(2*half)), x=xx.tolist(),truth=tt.tolist(),observation=yy.tolist(),mean=mean.tolist(),variance=var.tolist(),half_width=half.tolist()) runs[strategy]=dict(chosen=chosen,observations=pool_y[chosen].tolist(),learning_RMSE=curve, acquisitions=acquired,kernel_weights=model.weights.tolist(),weight_history=weight_history, calibration_scores=scores.tolist(),calibration_rank=rank,qhat=qhat,tests=tests) if seed==SEEDS[0]: grid_mean,grid_var=model.predict(POOL) runs[strategy]['grid_mean']=grid_mean.tolist();runs[strategy]['grid_variance']=grid_var.tolist() slice_x=np.column_stack([np.full(120,.5),np.linspace(0,1+.1/.35,120)]) sm,sv=model.predict(slice_x) runs[strategy]['slice']=dict(x=slice_x.tolist(),truth=physical(slice_x)[0].tolist(),mean=sm.tolist(),half_width=(qhat*np.sqrt(sv+NOISE**2)).tolist()) return dict(seed=seed,initial=initial,cal_x=cal_x.tolist(),cal_y=cal_y.tolist(),runs=runs) campaigns=[] for seed in SEEDS: campaigns.append(campaign(seed));print('Completed paired campaign',seed,flush=True) summary={} for strategy in STRATEGIES: curves=np.array([c['runs'][strategy]['learning_RMSE'] for c in campaigns]) summary[strategy]=dict(learning_mean=curves.mean(0).tolist(),learning_min=curves.min(0).tolist(),learning_max=curves.max(0).tolist(), final_RMSE_each=curves[:,-1].tolist(),final_RMSE_mean=float(curves[:,-1].mean()),final_RMSE_sd=float(curves[:,-1].std(ddof=1))) for key in ['nominal','shifted']: summary[strategy][key]={m:float(np.mean([c['runs'][strategy]['tests'][key][m] for c in campaigns])) for m in ['coverage','RMSE_logQ','mean_interval_width_logQ']} truth,phi,G=physical(POOL) results=dict(task='Explicit neutral-gel virtual laboratory, not fitted experimental formulations',seeds=SEEDS,budgets=BUDGETS, model=dict(lengthscales=LENGTHS.tolist(),prior_mean=PRIOR,amplitude_sd=AMP,noise_sd=NOISE,jitter=JITTER), costs=dict(training_queries=36,independent_calibration=40,total_per_campaign=76,nominal_test_audit=256,shifted_test_audit=128), pool=POOL.tolist(),pool_truth=truth.tolist(),pool_phi=phi.tolist(),pool_G_kPa=G.tolist(),summary=summary,campaigns=campaigns) (OUT/'results.json').write_text(json.dumps(results,ensure_ascii=False,indent=2)+'\n',encoding='utf-8') (OUT/'samples.json').write_text(json.dumps(dict(pool=results['pool'],truth=results['pool_truth'],campaigns=campaigns),ensure_ascii=False,indent=2)+'\n',encoding='utf-8') for name,row in summary.items(): print(f"{name:12s} final latent RMSE={row['final_RMSE_mean']:.5f}; nominal coverage={row['nominal']['coverage']:.3f}; shifted coverage={row['shifted']['coverage']:.3f}") print('Each campaign: 36 training measurements + 40 independent calibration measurements') print('All 12 paired campaigns retained; independent audit: 256 nominal + 128 shifted objects') plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False, 'axes.labelcolor':'#30453f','text.color':'#243d35','axes.edgecolor':'#bac7c0','xtick.color':'#53675e', 'ytick.color':'#53675e','svg.hashsalt':'L19-active','figure.facecolor':'#fffef9','axes.facecolor':'#fffef9'}) COLORS=['#a79883','#3f8297','#267b69'];LABELS=['Random','Space filling','Uncertainty'] first=campaigns[0]['runs'] def decorate(ax,title,x='',y=''): ax.set_title(title,loc='left',fontweight='bold',pad=14);ax.set_xlabel(x);ax.set_ylabel(y) ax.grid(axis='y',color='#e5eae1',lw=.7);ax.set_axisbelow(True) def draw(i,fig,spec): if i==4: grid=spec.subgridspec(2,1,hspace=.04);axes=[fig.add_subplot(grid[j,0]) for j in range(2)] for j,(key,ylab) in enumerate([('coverage','Fraction covered'),('mean_interval_width_logQ','Mean full width (log Q)')]): ax=axes[j] for k,(name,c) in enumerate(zip(STRATEGIES,COLORS)): values=[summary[name][domain][key] for domain in ['nominal','shifted']] ax.bar(np.arange(2)+(k-1)*.23,values,.21,color=c,label=LABELS[k]) for x,y in zip(np.arange(2)+(k-1)*.23,values):ax.text(x,y+.015,f'{y:.3f}',ha='center',fontsize=8) ax.set_xticks([0,1],['Same parameter domain','New solvent range']) decorate(ax,'E Read coverage together with interval width' if j==0 else '',y=ylab) if j==0:ax.axhline(.9,color='#ba8452',ls='--',lw=1);ax.set_ylim(0,1.15);ax.legend(frameon=False,ncol=3,fontsize=8,loc='lower left') else:ax.set_ylim(0,1.4) return ax=fig.add_subplot(spec) if i==0: im=ax.imshow(truth.reshape(25,25),origin='lower',extent=[0,1,.2,.55],aspect='auto',cmap='YlGnBu') ids=first['uncertainty']['chosen'];xy=POOL[ids] ax.scatter(xy[8:,0],.2+.35*xy[8:,1],s=25,facecolors='#d19a5c',edgecolors='white',lw=.8,label='28 sequential choices') ax.scatter(xy[:8,0],.2+.35*xy[:8,1],s=60,marker='+',c='white',lw=1.8,label='8 shared initial points') ax.set_xticks([0,.5,1],['2',f'{2*np.sqrt(30):.2f}','60']) decorate(ax,'A Learn a swelling surface with 36 measurements','Dry-reference strand density (mol / m³; log scale)','Solvent interaction parameter χ') ax.grid(False);ax.legend(frameon=True,facecolor='#244d56',labelcolor='white',fontsize=8,loc='lower right') fig.colorbar(im,ax=ax,label='Physical ln Q',shrink=.8) elif i==1: for name,label,c in zip(STRATEGIES,LABELS,COLORS): row=summary[name];ax.plot(BUDGETS,row['learning_mean'],'o-',c=c,label=label,ms=4) ax.fill_between(BUDGETS,row['learning_min'],row['learning_max'],color=c,alpha=.09) ax.set_yscale('log');ax.set_xticks(BUDGETS) decorate(ax,'B Compare complete learning campaigns','Training measurements','Latent-function RMSE (log Q)') ax.legend(frameon=False);ax.text(.03,.04,'Line: mean of 12 campaigns · band: full range',transform=ax.transAxes,fontsize=9) elif i==2: vals=np.array([summary[s]['final_RMSE_each'] for s in STRATEGIES]).T for row in vals:ax.plot(range(3),row,color='#b5beb2',lw=.8,alpha=.7,zorder=1) for k,(label,c) in enumerate(zip(LABELS,COLORS)): ax.scatter(np.full(12,k),vals[:,k],s=26,c=c,alpha=.65,zorder=2) ax.plot([k-.17,k+.17],[vals[:,k].mean()]*2,c=c,lw=4,zorder=3) ax.text(k,vals[:,k].max()+.012,f'mean {vals[:,k].mean():.4f}',ha='center',fontsize=9) ax.set_xticks(range(3),LABELS);ax.set_xlim(-.45,2.45);ax.set_ylim(0,.19) decorate(ax,'C Keep every paired repeat at the final budget',y='RMSE after 36 measurements (log Q)') ax.text(.04,.92,'Each thin line connects the same campaign seed',transform=ax.transAxes,fontsize=9) elif i==3: t=first['uncertainty']['tests']['nominal'];obs=np.array(t['observation']);mean=np.array(t['mean']);half=np.array(t['half_width']) ax.errorbar(obs,mean,yerr=half,fmt='none',ecolor='#a6bb9d',alpha=.45,lw=.6) ax.scatter(obs,mean,s=15,c=COLORS[2],alpha=.7);ax.plot([1.5,5.6],[1.5,5.6],c='#b78351',ls='--',lw=1) ax.set(xlim=(1.5,5.6),ylim=(1.5,5.6));decorate(ax,'D Check new observations and their intervals','New observed ln Q','Predicted ln Q with calibrated interval') ax.text(.04,.92,f"Seed 1900 · 256 held-out observations\nCoverage {t['coverage']:.1%}",transform=ax.transAxes,fontsize=9) else: row=first['uncertainty'];s=row['slice'];x=np.array(s['x']);chi=.2+.35*x[:,1];mu=np.array(s['mean']);half=np.array(s['half_width']) ax.axvspan(.55,.65,color='#ece2cc',alpha=.6,label='Beyond sampled solvent range') ax.fill_between(chi,mu-half,mu+half,color='#90b29d',alpha=.45,label='Calibrated observation interval') ax.plot(chi,s['truth'],c='#354a43',lw=2,label='Physical oracle');ax.plot(chi,mu,c='#267b69',ls='--',lw=2,label='GP mixture mean') ax.axvline(.55,color='#b78351',ls=':',lw=1) decorate(ax,'F Follow one slice into a new solvent range','χ at ν = 10.95 mol / m³','ln Q') ax.legend(frameon=False,fontsize=8,loc='upper right') ax.text(.03,.04,'Final kernel weights (short / medium / long):\n'+ ' / '.join(f'{v:.3f}' for v in row['kernel_weights']),transform=ax.transAxes,fontsize=8) keys=['surface','learning','paired','intervals','coverage','slice'] for i,key in enumerate(keys): fig=plt.figure(figsize=(8.6,5.5),layout='constrained');draw(i,fig,fig.add_gridspec(1,1)[0]);fig.savefig(OUT/f'{key}.svg',metadata={'Date':None});plt.close(fig) fig=plt.figure(figsize=(18,16.5),layout='constrained');grid=fig.add_gridspec(3,2,hspace=.12,wspace=.1) for i in range(6):draw(i,fig,grid[i//2,i%2]) fig.savefig(OUT/'active_learning.svg',metadata={'Date':None});fig.savefig(OUT/'active_learning.png',dpi=160,metadata={'Software':'L19 reproducible example'});plt.close(fig) print('Full campaigns, calibration scores and six figures: active_outputs/')