"""Reproduce all computations used in the curved-fold paper.""" from __future__ import annotations import argparse, csv, json, os, platform, sys, time from dataclasses import replace from pathlib import Path import numpy as np import scipy ROOT=Path(__file__).resolve().parents[1] os.environ.setdefault('MPLCONFIGDIR',str(ROOT/'.cache'/'matplotlib')) import matplotlib matplotlib.use('Agg') matplotlib.rcParams['pdf.fonttype']=42 import matplotlib.pyplot as plt from fold_model import (Parameters,mellin_coefficient,mellin_quadrature, optimized_jump_bound,power_generator,power_drift_constant, near_jump_bound,critical_barrier_constants,critical_barrier_generator) from simulate_fold import simulate,summarize from path_data import write_paths def write_csv(path,rows): with path.open('w',newline='') as fp: writer=csv.DictWriter(fp,fieldnames=list(rows[0]));writer.writeheader();writer.writerows(rows) def main(): ap=argparse.ArgumentParser();ap.add_argument('--skip-simulation',action='store_true') ap.add_argument('--n',type=int,default=20000) ap.add_argument('--output-root',type=Path,default=ROOT.parent, help='Write data, tables and figures below this directory.') args=ap.parse_args() if args.n<1:ap.error('--n must be positive') start=time.time();data=args.output_root/'anc'/'data';paper=args.output_root data.mkdir(parents=True,exist_ok=True) for d in ['tables','figures']:(paper/d).mkdir(parents=True,exist_ok=True) par=Parameters();rows=[] for q in [.2,.6,.85]: for p in [-(1-q)/2,.1,.2]: direct=mellin_quadrature(p,q);formula=mellin_coefficient(p,q) rows.append(dict(q=q,p=p,quadrature=direct,gamma_formula=formula, absolute_error=abs(direct-formula))) write_csv(data/'mellin_quadrature.csv',rows) maxerr=max(r['absolute_error'] for r in rows) checks=[] for d,p in [(1.2,.2),(3.,-.1)]: pa=replace(par,delta0=d);h=np.geomspace(1e-10,1e5,701)[:,None] v=np.linspace(-20,20,63)[None,:];apower=power_drift_constant(p,pa) normalized=power_generator(h,v,p,pa)/h**(p-1) checks.append(dict(check='power',dimension=d,q=pa.q,p=p, tested_points=int(normalized.size),maximum=float(normalized.max()), threshold=-apower)) assert normalized.max()<=-apower+1e-10 for q in [.2,.6,.85]: pa=replace(par,alpha=2*q,delta0=2.2);L,K=critical_barrier_constants(pa) h=np.geomspace(L*1e-10,L,300)[:,None];v=np.linspace(-20,20,51)[None,:] val=critical_barrier_generator(h,v,pa)*h**q checks.append(dict(check='critical',dimension=2.2,q=q,p=0, tested_points=int(val.size),maximum=float(val.max()),threshold=0.)) assert val.max()<=1e-9 write_csv(data/'generator_checks.csv',checks) bounds=[] for h in [1e-2,1e-4,1e-6,1e-8]: for kappa in [.01,.0625,.25,1.]: bound,L=optimized_jump_bound(h,.2,replace(par,kappa=kappa)) bounds.append(dict(h0=h,kappa=kappa,jump_upper=bound,contact_lower=1-bound,optimal_collar=L)) write_csv(data/'curvature_bounds.csv',bounds) tex=[] for r in bounds: if r['h0'] in [1e-4,1e-8]: tex.append(f"$10^{{{int(round(np.log10(r['h0'])))}}}$ & {r['kappa']:g} & {r['jump_upper']:.5f} & {r['contact_lower']:.5f} & {r['optimal_collar']:.5g} "+r"\\") (paper/'tables/curvature.tex').write_text('\n'.join(tex)+'\n') fig,ax=plt.subplots(figsize=(6.2,3.8)) ks=np.geomspace(1e-4,1,100) for h0 in [1e-2,1e-4,1e-6]: vals=[optimized_jump_bound(h0,.2,replace(par,kappa=float(k)))[0] for k in ks] ax.loglog(ks,vals,label=fr'$h_0=10^{{{int(round(np.log10(h0)))}}}$') ax.set_xlabel(r'Quadratic curvature parameter $\kappa$') ax.set_ylabel('Upper bound on eventual jump probability') ax.legend();ax.grid(True,which='both',alpha=.25);fig.tight_layout() fig.savefig(paper/'figures/curvature_bound.pdf');plt.close(fig) fig,ax=plt.subplots(figsize=(6.2,3.8));eps=np.geomspace(1e-5,.1,100) for d,p,label in [(1.2,.2,'Subcritical, $p=0.2$'),(3.,-.1,'Supercritical, $p=-0.1$')]: vals=near_jump_bound(.05,eps,p,replace(par,delta0=d)) ax.loglog(eps,vals,label=label) ax.set_xlabel(r'Pre-terminal radius $\epsilon$');ax.set_ylabel('Upper bound on near-fold jump termination') ax.legend();ax.grid(True,which='both',alpha=.25);fig.tight_layout() fig.savefig(paper/'figures/near_terminal_bound.pdf');plt.close(fig) configs=[('sub_c08_d128',1.2,.25,.08,1/128), ('sub_c04_d128',1.2,.25,.04,1/128), ('sub_c04_d256',1.2,.25,.04,1/256), ('sub_c02_d256',1.2,.25,.02,1/256), ('critical_c04_d256',2.2,.25,.04,1/256), ('super_c08_d128',3.,.25,.08,1/128), ('super_c04_d256',3.,.25,.04,1/256), ('super_c02_d256',3.,.25,.02,1/256), ('sub_flat_c02_d256',1.2,.0625,.02,1/256), ('sub_curved_c02_d256',1.2,1.,.02,1/256)] summary=[] for i,(name,d,kappa,cutoff,step) in enumerate(configs): fp=(ROOT/'data' if args.skip_simulation else data)/(name+'.json') if not args.skip_simulation: pa=replace(par,delta0=d,kappa=kappa) paths,meta=simulate(args.n,pa,seed=20261020+1009*i,cutoff=cutoff,max_step=step) write_paths(data/(name+'.csv'),paths) meta['summary']=summarize(paths);fp.write_text(json.dumps(meta,indent=2)) elif not fp.exists(): continue else:meta=json.loads(fp.read_text()) sm=meta['summary'] row=dict(name=name,dimension=d,kappa=kappa,cutoff=cutoff,max_step=step,n=meta['n']) for cause in ['contact','jump','survival']: for key in ['count','estimate','standard_error']:row[cause+'_'+key]=sm[cause][key] row['omitted_normal_mean_upper']=meta['omitted_normal_mean_upper'] summary.append(row);print(name,sm,flush=True) if summary: write_csv(data/'simulation_summary.csv',summary) tex=[] for r in summary: tex.append(f"{r['dimension']:g} & {r['kappa']:g} & {r['cutoff']:g} & 1/{round(1/r['max_step'])} & {r['contact_estimate']:.5f} & {r['jump_estimate']:.5f} & {r['survival_estimate']:.5f} "+r"\\") (paper/'tables/simulations.tex').write_text('\n'.join(tex)+'\n') meta=dict(python=sys.version,numpy=np.__version__,scipy=scipy.__version__,platform=platform.platform(), max_mellin_absolute_error=maxerr,generator_points=sum(r['tested_points'] for r in checks), total_simulated_paths=sum(r['n'] for r in summary),elapsed_seconds=time.time()-start) (data/'reproduction.json').write_text(json.dumps(meta,indent=2));print(meta,flush=True) if __name__=='__main__':main()