#!/usr/bin/env python3 """ R1 - Closed-form and Monte-Carlo numbers used in the refined hypotheses (refined/R1..R30). Nothing here is an experimental observation. Sections A OD600 stray-light compression and the bias in fitted growth rate B Jensen bias when a mean response is inverted through a curved calibration C Spurious ratio-denominator correlation (Pearson) D Weibull area scaling of dielectric breakdown strength E Unstirred-water-layer series resistance in permeability assays F DSC thermal-lag shift of an onset temperature G Error-bar deflation under a forgotten linear trend H Inverse-variance weighting bias when variances are estimated (Monte Carlo) I Two-boundary TDTR sensitivity weight Run: python3 investigations/R1_refinement_numbers.py """ import numpy as np, json, os, csv rng = np.random.default_rng(11) HERE = os.path.dirname(os.path.abspath(__file__)); ROOT = os.path.dirname(HERE) OUT = os.path.join(ROOT, 'data'); os.makedirs(OUT, exist_ok=True) res = {} # ---------- A. OD600 stray light ---------- def od_meas(od_true, stray): # I/I0 = 10^-od ; measured includes an additive stray fraction of I0 return -np.log10(10.0**(-od_true) + stray) def mu_bias(od_lo, od_hi, stray, n=41): """apparent slope of ln(OD) vs true OD, in the exponential phase, relative to the true slope (OD proportional to cell number, so ln OD = mu t + const).""" od = np.linspace(od_lo, od_hi, n) return np.polyfit(od, np.log(od_meas(od, stray)/od), 1)[0] A = {} for stray in [1e-4, 1e-3, 1e-2]: for (lo, hi) in [(0.05, 0.3), (0.1, 0.6), (0.2, 1.0), (0.3, 1.4)]: A[f'stray={stray:g}_window={lo}-{hi}'] = round(float(mu_bias(lo, hi, stray)), 4) res['A_OD_bias_fraction_per_OD'] = A # high-OD regime, where stray light actually bites A2 = {} for stray in [1e-3, 1e-2]: for (lo, hi) in [(1.0, 2.0), (1.5, 2.5), (2.0, 3.0)]: A2[f'stray={stray:g}_window={lo}-{hi}'] = round(float(np.exp(mu_bias(lo,hi,stray)*(hi-lo))-1)*100, 1) res['A2_suppression_pct_high_OD'] = A2 res['A_example_suppression_pct'] = { 'stray=1e-3, window 0.2-1.0': round(100*float(np.exp(mu_bias(0.2,1.0,1e-3)*0.8)-1), 1), 'stray=1e-2, window 0.2-1.0': round(100*float(np.exp(mu_bias(0.2,1.0,1e-2)*0.8)-1), 1)} # ---------- B. Jensen bias ---------- # calibration y = a x + b x^2 ; x has mean m and variance s^2 ; we invert B = {} for b_over_a in [0.05, 0.2, 0.5]: # curvature scale, 1/unit for cv in [0.1, 0.3, 0.6]: # x = m(1+eps), eps ~ N(0,cv^2), m=1 m, s = 1.0, cv bias = b_over_a*s**2/(1+2*b_over_a*m) # first-order inverse bias B[f'b/a={b_over_a}_cv={cv}'] = round(float(bias), 4) res['B_Jensen_bias_in_x_units'] = B # ---------- C. spurious ratio-denominator correlation ---------- C = {} for cvx, cvy in [(0.1,0.1),(0.3,0.3),(0.1,0.3),(0.3,0.1),(0.5,0.2)]: r = -cvy/np.sqrt(cvx**2+cvy**2+cvx**2*cvy**2) C[f'cvX={cvx}_cvY={cvy}'] = round(float(r),3) res['C_spurious_corr'] = C # ---------- D. Weibull area scaling ---------- D = {} for m in [5, 10, 20]: for ratio in [10, 100]: D[f'm={m}_A_ratio={ratio}'] = round(float(ratio**(-1.0/m)), 3) res['D_scale_factor'] = D # ---------- E. UWL ---------- # P_app^-1 = P_true^-1 + P_uwl^-1 ; P_uwl = D/h, D~5e-6 cm2/s, h in um D_ = 5e-6 E = {} for ptrue in [5e-6, 5e-5, 5e-4]: # cm/s row = {} for h_um in [200, 500, 1000, 2000]: puwl = D_/(h_um*1e-4) papp = 1/(1/ptrue + 1/puwl) row[f'h={h_um}um'] = round(float(papp/ptrue), 3) E[f'Ptrue={ptrue:.0e}cm/s'] = row res['E_Papp_over_Ptrue'] = E # ---------- F. DSC lag ---------- # onset shift = tau * beta, tau = m*cp*Rth the effective pan-sample time constant F = {} for tau in [5, 10, 20]: for beta in [5, 10, 20]: F[f'tau={tau}s_beta={beta}K/min'] = round(float(tau*beta/60.0), 2) res['F_DSC_onset_shift_K'] = F # ---------- G. forgotten trend ---------- G = {} # y = a t + noise, t in [0,1]; naive SE of mean uses residual sd/sqrt(N) for a_over_sigma in [0.5, 1, 2, 4]: # var of mean of a linear ramp = a^2/12 ; residual var ~ sigma^2 true_se = np.sqrt(a_over_sigma**2/12 + 1.0) naive_se = 1.0 G[f'trend/sigma={a_over_sigma}'] = round(float(true_se/naive_se), 2) res['G_SE_underestimate_factor'] = G # ---------- H. inverse-variance weighting with estimated variances ---------- H = {} N = 40000 for k in [2, 5, 10, 30]: est = [] for _ in range(N): s2 = rng.chisquare(k, size=4)/k # 4 measurements, k dof each x = rng.normal(0, 1, size=4) w = 1/s2 est.append((w*x).sum()/w.sum()) est = np.array(est) H[f'dof={k}'] = dict(mean=round(float(est.mean()),3), sd=round(float(est.std()),3), quoted_se=round(float(1/np.sqrt(4)),3), quoted_over_true=round(float((1/np.sqrt(4))/est.std()),3)) res['H_weighting_bias'] = H # ---------- I. two-boundary TDTR sensitivity (thermal quadrupole) ---------- # 1D frequency-domain multilayer: surface temperature response to a modulated # flux, with two interface conductances G1 (transducer/film) and G2 (film/substrate). def quadrupole(k, C, d, w): gam = np.sqrt(1j*w*C/k) ch, sh = np.cosh(gam*d), np.sinh(gam*d) return np.array([[ch, -sh/(k*gam)], [-k*gam*sh, ch]]) def resistance(R): return np.array([[1.0, -R], [0.0, 1.0]]) def surf_Z(f, G1, G2, kf=1.0, Cf=1.6e6, df=100e-9, kt=150.0, Ct=2.4e6, dt=80e-9, ks=150.0, Cs=1.6e6): w = 2*np.pi*f M = quadrupole(kt, Ct, dt, w) # transducer M = M @ resistance(1.0/G1) # transducer/film interface M = M @ quadrupole(kf, Cf, df, w) # film M = M @ resistance(1.0/G2) # film/substrate interface gam_s = np.sqrt(1j*w*Cs/ks) Zsub = 1.0/(ks*gam_s) return (M[0,0]*Zsub + M[0,1])/(M[1,0]*Zsub + M[1,1]) I = {} G1, G2 = 1.0e8, 5.0e7 for f in [1e5, 3e5, 1e6, 3e6, 1e7, 3e7]: h = 1e-4 Z0 = surf_Z(f, G1, G2) s1 = abs((np.log(surf_Z(f, G1*(1+h), G2))-np.log(surf_Z(f, G1*(1-h), G2)))/(2*h)) s2 = abs((np.log(surf_Z(f, G1, G2*(1+h)))-np.log(surf_Z(f, G1, G2*(1-h))))/(2*h)) # normalised sensitivity matrix condition number for the pair J = np.array([[s1, s2]]) cond = float(np.linalg.cond(J @ J.T)) I[f'f={f:.0e}Hz'] = dict(S_G1=round(float(s1),3), S_G2=round(float(s2),3), ratio=round(float(s2/s1),3)) res['I_TDTR_sensitivity'] = I json.dump(res, open(os.path.join(OUT,'r1_refinement_numbers.json'),'w'), indent=1) print(json.dumps(res, indent=1))