Unexecuted reproducibility supplement

def field_layer(count,cut,rtol,name): lo,hi = .08,.24 omega = np.linspace(lo,hi,count) width,mid = hi-lo,(hi+lo)/2 h = np.column_stack([np.ones(count)/np.sqrt(width), np.sqrt(12)*(omega-mid)/width**1.5]) bound = radial_tail(cut,0,lo) jwrite(name+"_radial_tail",bound) runs = [radial(om,0,cut,rtol) for om in omega] tau = np.array([r[0] for r in runs]) reflection = np.array([r[1] for r in runs]) S = np.array([r[2] for r in runs]) flux = np.array([r[3] for r in runs]) unitary = np.array([r[4] for r in runs]) weights = simpson_weights(omega) q = np.abs(tau)**2/np.expm1(4*np.pi*omega) K = h.conj().T@((weights*q)[:,None]*h) energy = h.conj().T@((weights*omega)[:,None]*h) save(name,omega=omega,tau=tau,reflection=reflection,S=S, K=K,C=K.T,H=energy,packets=h,weights=weights, flux_residual=flux,unitarity_residual=unitary, ODE_function_evaluations=np.array([r[5] for r in runs])) jwrite(name+"_basis","packet_order":["constant","linear"], "frequency_grid_order":"ascending", "S_rows":["outgoing horizon","outgoing infinity"], "S_columns":["incoming horizon","incoming infinity"], "K_indices":"one-particle ket, bra","C":"K transpose", "radial_tail_artifact":name+"_radial_tail.json") need(max(flux)<1e-5 and max(unitary)<1e-5,"radial flux residual") check(name+"_packet_gram",h.conj().T@(weights[:,None]*h),I,1e-10) check(name+"_energy",energy, np.array([[mid,width/np.sqrt(12)],[width/np.sqrt(12),mid]]),1e-10) return K,energy,"id":name,"grid":count,"cut":cut,"rtol":rtol, "radial_tail":bound,"max_flux_residual":float(max(flux)), "max_unitarity_residual":float(max(unitary)), "certified_total_amplitude_error":False

def field_pilot(): K17,H17,r17 = field_layer(17,80,1e-10,"field_grid17") K33,H33,r33 = field_layer(33,80,1e-10,"field_grid33") K65,H65,r65 = field_layer(65,80,1e-10,"field_grid65") Kb,Hb,rb = field_layer(33,40,1e-10,"field_boundary40") Kt,Ht,rt = field_layer(33,80,1e-7,"field_coarse_tolerance") rho,basis,tail,p,Cn = gaussian_cutoff(K65,2,"field_finite_state") a = 1-1/float(np.linalg.det(I+K65).real) mu = float(np.trace(K65).real) DATA["field"] = {"dataset_id":"field_grid65", "parameters":"r_h":1,"beta":4*np.pi,"ell":0, "band":[.08,.24],"J":2,"N":2, "runs":[r17,r33,r65,rb,rt], "diagnostics_not_certificates": "quadrature_17_33":float(np.linalg.norm(K17-K33,2)), "quadrature_33_65":float(np.linalg.norm(K33-K65,2)), "boundary_40_80":float(np.linalg.norm(Kb-K33,2)), "ODE_tolerance":float(np.linalg.norm(Kt-K33,2)), "mean_occupation":mu,"selected_energy":float(np.trace(H65@K65).real), "normal_mode_means":np.linalg.eigvalsh(K65), "Gaussian_tail_estimate":tail,"Markov_tail":min(1,mu/3), "discarded_mean":mu-float(np.dot(np.arange(3),p)), "offdiagonal_magnitude":float(abs(K65[0,1])), "a":a,"benchmark_uses_only":"a, not full K or finite state"} dump() return a