""" Impact of white noise amplitude on Prony series fit scatter. Supplementary data figure. This code is first generated by a Claude Opus model and modified and vetted by Joshua A. Dijksman. Model: G(t) = E0 + E1*exp(-t/tau1) + E2*exp(-t/tau2) + E3*exp(-t/tau3) fitted as: sigma(t) = C_inf + A1*exp(-t/tau1) + A2*exp(-t/tau2) + A3*exp(-t/tau3) Parameters match experiment: tau = [20, 500, 1500] s, ramp_time = 20 s, hold = 3000 s (~ 1.3 x tau3) For each noise level: draw N_mc realisations of Gaussian noise on the hold-phase stress, refit all 7 parameters, collect the distribution. """ import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from scipy.optimize import curve_fit, nnls import warnings warnings.filterwarnings('ignore') # ============================================================================= # Parameters — edit these to match your experiment # ============================================================================= TAU = np.array([20., 500., 1500.]) # relaxation times [s] E0 = 50. # permanent spring modulus [Pa] E = np.array([300., 150., 100.]) # Maxwell arm moduli [Pa] ETA = E * TAU # viscosities [Pa·s] EPS = 0.05 # applied strain [-] RAMP_TIME = 20. # ramp duration [s] HOLD_TIME = 3000. # hold duration [s] # Noise levels to sweep: std as fraction of peak hold-phase stress NOISE_FRACS = np.array([0.001, 0.002, 0.005, 0.01, 0.02, 0.05]) N_MC = 400 # Monte Carlo draws per noise level RNG_SEED = 42 # ============================================================================= # Simulation (forward Euler, minimal) # ============================================================================= def simulate_hold(n_pts=4000): """Return (t_hold, sigma_hold) — time and stress on the hold phase only.""" t_total = RAMP_TIME + HOLD_TIME t = np.linspace(0, t_total, n_pts) dt = t[1] - t[0] strain = np.where(t <= RAMP_TIME, EPS * t / RAMP_TIME, EPS) ev = np.zeros(3) sigma = np.zeros(n_pts) for i in range(n_pts): ee = strain[i] - ev smx = E * ee sigma[i] = E0 * strain[i] + smx.sum() if i < n_pts - 1: ev = ev + (smx / ETA) * dt mask = t >= RAMP_TIME return t[mask] - RAMP_TIME, sigma[mask] # ============================================================================= # Fit model and helpers # ============================================================================= def model(t, Ci, A1, A2, A3, t1, t2, t3): return Ci + A1*np.exp(-t/t1) + A2*np.exp(-t/t2) + A3*np.exp(-t/t3) def initial_guess(t, s): """NNLS warm-start for [C_inf, A1, A2, A3] with tau fixed, then perturb tau.""" Phi = np.column_stack([np.ones_like(t), np.exp(-t/TAU[0]), np.exp(-t/TAU[1]), np.exp(-t/TAU[2])]) c, _ = nnls(Phi, s) p0 = [c[0]*0.9, c[1]*0.9, c[2]*0.9, c[3]*0.9, TAU[0]*1.2, TAU[1]*1.2, TAU[2]*1.2] return p0 # ============================================================================= # Main computation # ============================================================================= t_hold, s_clean = simulate_hold() # Log-spaced subsample — captures all timescales, fits run fast pos = t_hold > 0 t_pos = t_hold[pos]; s_pos = s_clean[pos] n_sub = min(350, len(t_pos)) sub_idx = np.unique(np.round( np.logspace(0, np.log10(len(t_pos)-1), n_sub)).astype(int)) t_sub = t_pos[sub_idx]; s_sub = s_pos[sub_idx] peak = s_clean[0] p0_base = initial_guess(t_sub, s_sub) bounds = ([0,0,0,0, 0.5, 5., 50.], [np.inf, np.inf, np.inf, np.inf, TAU[0]*10, TAU[1]*10, TAU[2]*10]) fit_kw = dict(maxfev=2000, ftol=1e-7, xtol=1e-7) # True parameter values (ideal step) ref = np.array([E0*EPS, E[0]*EPS, E[1]*EPS, E[2]*EPS, TAU[0], TAU[1], TAU[2]]) # p_store[noise_level, mc_draw, param] p_store = np.full((len(NOISE_FRACS), N_MC, 7), np.nan) rng = np.random.default_rng(RNG_SEED) for ni, nf in enumerate(NOISE_FRACS): sigma_n = nf * peak p0 = p0_base.copy() ok = 0 for attempt in range(N_MC * 6): noise = rng.normal(0, sigma_n, len(t_sub)) try: pf, _ = curve_fit(model, t_sub, s_sub + noise, p0=p0, bounds=bounds, **fit_kw) p0 = list(np.maximum(pf, [0.01]*4 + [0.5]*3)) p_store[ni, ok] = pf ok += 1 if ok >= N_MC: break except RuntimeError: pass print(f"noise={nf*100:.1f}% fits={ok}/{N_MC} " f"CV: C_inf={100*np.nanstd(p_store[ni,:ok,0])/ref[0]:.0f}% " f"tau1={100*np.nanstd(p_store[ni,:ok,4])/ref[4]:.0f}% " f"tau2={100*np.nanstd(p_store[ni,:ok,5])/ref[5]:.0f}% " f"tau3={100*np.nanstd(p_store[ni,:ok,6])/ref[6]:.0f}%") # ============================================================================= # Figure # ============================================================================= PNAMES = ['C_inf\n(=E₀ε)', 'A₁\n(=E₁ε)', 'A₂\n(=E₂ε)', 'A₃\n(=E₃ε)', 'τ₁ (20 s)', 'τ₂ (500 s)', 'τ₃ (1500 s)'] COLORS = ['#222222', '#2176ae', '#c1121f', '#2d9e3a', '#7bb3d4', '#e07a80', '#7ed1a0'] noise_pct = NOISE_FRACS * 100 fig = plt.figure(figsize=(14, 9)) fig.suptitle( f'Effect of white-noise amplitude on Prony fit scatter\n' f'τ = {list(TAU.astype(int))} s, ramp = {RAMP_TIME:.0f} s, ' f'hold = {HOLD_TIME:.0f} s ({HOLD_TIME/TAU[2]:.1f}×τ₃), ' f'N = {N_MC} MC draws per noise level', fontsize=10, y=0.99, ) gs = gridspec.GridSpec(2, 2, figure=fig, hspace=0.44, wspace=0.32) # ── (1) CV vs noise level ──────────────────────────────────────────────────── ax1 = fig.add_subplot(gs[0, 0]) for k in range(7): cv = np.array([100 * np.nanstd(p_store[ni, :, k]) / ref[k] for ni in range(len(NOISE_FRACS))]) ax1.plot(noise_pct, cv, color=COLORS[k], lw=2, marker='o', ms=5, ls='-' if k < 4 else '--', label=PNAMES[k]) ax1.axhline(10, color='#999', ls=':', lw=1.2) ax1.axhline(20, color='#999', ls=':', lw=1.2) ax1.text(noise_pct[-1]*0.88, 10.6, '10%', color='#888', fontsize=7, ha='right') ax1.text(noise_pct[-1]*0.88, 20.6, '20%', color='#888', fontsize=7, ha='right') ax1.set_xscale('log') ax1.set_yscale('log') ax1.set_ylim(bottom=0) ax1.set(xlabel='Noise std [% of peak hold stress]', ylabel='Ratio of true / fitted value [%]')#, # title='(1) Parameter scatter vs noise level\n' # 'Solid = amplitudes, dashed = timescales') #ax1.legend(fontsize=7, ncol=2, loc='upper left') #ax1.grid(True, which='both', alpha=0.28) # ── (2) Histograms at 3 representative noise levels ───────────────────────── ax2 = fig.add_subplot(gs[0, 1]) noise_show = [1, 2, 4] # indices into NOISE_FRACS (0.5%, 1%, 2%) hatch_list = ['', '///', 'xxx'] # Show only tau3 (most diagnostic) and A1 (ramp-dominated) for k, (ki, pname, col) in enumerate([(6, 'τ₃', '#7ed1a0'), (4, 'τ₁', '#7bb3d4'), (0, 'C_inf', '#222222')]): for j, (ni, hatch) in enumerate(zip(noise_show, hatch_list)): data = p_store[ni, :, ki] data = data[~np.isnan(data)] label = f'{pname} ({NOISE_FRACS[ni]*100:.1f}% noise)' if j == 0 else None # normalise by true value so all params on same axis ax2.hist(data / ref[ki], bins=30, density=True, alpha=0.45, color=col, edgecolor='none', hatch=hatch, label=label if j==0 else f' {NOISE_FRACS[ni]*100:.1f}% noise') ax2.axvline(1.0, color='k', lw=1.5, ls='--', label='True value') ax2.set(xlabel='Fitted / true value', ylabel='Density', title='(2) Distribution of fitted τ₃, τ₁, C_inf\n' 'at three noise levels (normalised by true value)') ax2.legend(fontsize=6.5, loc='upper right') ax2.grid(True, axis='y', alpha=0.28) # ── (3) Scatter plot: tau2 vs tau3 at two noise levels ────────────────────── ax3 = fig.add_subplot(gs[1, 0]) scatter_ni = [1, 3] # 0.2% and 1% noise scatter_c = ['#2176ae', '#c1121f'] scatter_a = [0.5, 0.4] for ni, sc, sa in zip(scatter_ni, scatter_c, scatter_a): d2 = p_store[ni, :, 5]; d3 = p_store[ni, :, 6] ok = ~(np.isnan(d2) | np.isnan(d3)) ax3.scatter(d2[ok], d3[ok], s=8, color=sc, alpha=sa, linewidths=0, label=f'{NOISE_FRACS[ni]*100:.1f}% noise (N={ok.sum()})') ax3.axvline(TAU[1], color='k', ls='--', lw=1.2) ax3.axhline(TAU[2], color='k', ls='--', lw=1.2) ax3.plot(TAU[1], TAU[2], 'k+', ms=12, mew=2, label='True values') ax3.set(xlabel='Fitted τ₂ [s]', ylabel='Fitted τ₃ [s]', title='(3) Joint scatter of τ₂ and τ₃\n' 'Correlated — fitter trades off slow timescales against each other') ax3.legend(fontsize=7, loc='upper right') ax3.grid(True, alpha=0.28) ax3.set_yscale('log') ax3.set_xscale('log') # ── (4) Stress curve + noise realisations + fit band ──────────────────────── ax4 = fig.add_subplot(gs[1, 1]) # show at 1% noise (index 3) ni_demo = 3 sigma_n = NOISE_FRACS[ni_demo] * peak rng2 = np.random.default_rng(0) # fit curves from stored parameters fits_demo = np.array([ model(t_pos, *p_store[ni_demo, i]) for i in range(N_MC) if not np.any(np.isnan(p_store[ni_demo, i])) ]) blo = np.percentile(fits_demo, 2.5, axis=0) bhi = np.percentile(fits_demo, 97.5, axis=0) bmed = np.median(fits_demo, axis=0) # a few noisy realisations for visual context for _ in range(8): nd = rng2.normal(0, sigma_n, len(t_sub)) ax4.semilogx(t_sub, (s_sub + nd), '.', color='#bbbbbb', ms=2, alpha=0.6) ax4.fill_between(t_pos, blo, bhi, color='#a8c8f0', alpha=0.6, label='95% fit band') ax4.semilogx(t_pos, bmed, 'b-', lw=2, label='MC median fit') ax4.semilogx(t_pos, s_pos, 'k-', lw=2, label='True (no noise)') ax4.axhline(E0*EPS, color='steelblue', ls=':', lw=1.3, label=f'True G(∞) = {E0*EPS:.2f} Pa') ax4.axvline(HOLD_TIME, color='gray', ls='--', lw=1.2) ax4.text(HOLD_TIME*0.88, bhi[-1], f'End of\nexperiment\n({HOLD_TIME:.0f} s)', color='gray', fontsize=7, ha='right', va='top') for tk, lbl in zip(TAU, ['τ₁','τ₂','τ₃']): ax4.axvline(tk, color='silver', ls='-', lw=0.7, alpha=0.6) ax4.text(tk*1.06, s_pos[0]*0.97, lbl, color='gray', fontsize=7, va='top') ax4.set(xlabel='Time since ramp end [s] (log scale)', ylabel='Stress [Pa]', title=f'(4) Noisy data + fit band at {NOISE_FRACS[ni_demo]*100:.0f}% noise\n' f'Band widens at long times — slow arm poorly constrained') ax4.legend(fontsize=7, loc='upper right') ax4.grid(True, which='both', alpha=0.28) plt.savefig('noise_impact.png', dpi=150, bbox_inches='tight') print('\nFigure saved: noise_impact.png') plt.show() # second simple figure fig = plt.figure(figsize=(14, 9)) gs = gridspec.GridSpec(1, 2, figure=fig, hspace=0.44, wspace=0.32) # ── (1) CV vs noise level ──────────────────────────────────────────────────── ax1 = fig.add_subplot(gs[0, 0]) for k in range(7): cv = np.array([100 * np.nanstd(p_store[ni, :, k]) / ref[k] for ni in range(len(NOISE_FRACS))]) ax1.plot(noise_pct, cv, color=COLORS[k], lw=2, marker='o', ms=5, ls='-' if k < 4 else '--', label=PNAMES[k]) #ax1.axhline(10, color='#999', ls=':', lw=1.2) #ax1.axhline(20, color='#999', ls=':', lw=1.2) #ax1.text(noise_pct[-1]*0.88, 10.6, '10%', color='#888', fontsize=7, ha='right') #ax1.text(noise_pct[-1]*0.88, 20.6, '20%', color='#888', fontsize=7, ha='right') ax1.set_xscale('log') ax1.set_yscale('log') ax1.set_ylim(bottom=0) ax1.set(xlabel='Noise std [% of peak hold stress]', ylabel='Ratio of true / fitted value [%]') ax1.xaxis.label.set_size(20) ax1.yaxis.label.set_size(20) ax1.tick_params(axis='both', labelsize=20) # ── (3) Scatter plot: tau2 vs tau3 at two noise levels ────────────────────── ax3 = fig.add_subplot(gs[0, 1]) scatter_ni = [1, 3] # 0.2% and 1% noise scatter_c = ['#2176ae', '#c1121f'] scatter_a = [0.5, 0.4] for ni, sc, sa in zip(scatter_ni, scatter_c, scatter_a): d2 = p_store[ni, :, 5]; d3 = p_store[ni, :, 6] ok = ~(np.isnan(d2) | np.isnan(d3)) ax3.scatter(d2[ok], d3[ok], s=8, color=sc, alpha=sa, linewidths=0, label=f'{NOISE_FRACS[ni]*100:.1f}% noise (N={ok.sum()})') ax3.axvline(TAU[1], color='k', ls='--', lw=1.2) ax3.axhline(TAU[2], color='k', ls='--', lw=1.2) ax3.plot(TAU[1], TAU[2], 'k+', ms=12, mew=2, label='True values') ax3.set(xlabel='Fitted τ₂ [s]', ylabel='Fitted τ₃ [s]') ax3.xaxis.label.set_size(20) ax3.yaxis.label.set_size(20) ax3.tick_params(axis='both', labelsize=20) ax3.grid(True, alpha=0.28) ax3.set_yscale('log') ax3.set_xscale('log') ax1.text(0.02, 0.98, '(a)', transform=ax1.transAxes, fontsize=20, va='top', ha='left') ax3.text(0.02, 0.98, '(b)', transform=ax3.transAxes, fontsize=20, va='top', ha='left') plt.savefig('noise_impact_simple.png', dpi=150, bbox_inches='tight') print('\nFigure saved: noise_impact_simple.png') plt.show()