| |
| """Publication figures for the HFX hierarchical fracture study. Run after the study stages; re-run after the holdouts.""" |
| import os, json, sys, math |
| import numpy as np, pandas as pd |
| import matplotlib; matplotlib.use('Agg') |
| import matplotlib.pyplot as plt |
| from matplotlib import cm, colors, gridspec |
| from matplotlib.lines import Line2D |
| from hfx_analysis import load_log, agg, load_run, savefig, ROOT, FIG, PROC, MODE_COL |
| from render_structure import render |
| import matplotlib.image as mpimg |
|
|
| df = load_log(); df = df[df.status == 'ok'].copy() |
| LS = os.path.join(PROC, 'load_sharing.csv') |
| if os.path.exists(LS): df = df.set_index('runId').join(pd.read_csv(LS).set_index('runId')).reset_index() |
| df['f0'] = df['levelMaterial'].apply(lambda l: l[0] if isinstance(l, list) and len(l) else np.nan) |
| df['Wf_uJ'] = df['Wf'] * 1e6; df['Fpeak_m'] = df['Fpeak'] * 1e3 |
| A = agg(df) |
| A['f0'] = df.groupby('designId')['f0'].mean() |
| A['dml'] = df.groupby('designId')['p_disorderMinLevel'].first() |
| A['fmode'] = df.groupby('designId')['fractureMode'].agg(lambda s: s.mode().iloc[0]) |
| main = df[df.stage != 'holdout'] |
| print('runs', len(df), 'designs', df.designId.nunique(), 'seeds per design: min', df.groupby('designId').size().min(), 'median', df.groupby('designId').size().median()) |
|
|
| def errbar(ax, x, y, ystd, **kw): |
| ax.errorbar(x, y, yerr=ystd, capsize=2, **kw) |
|
|
| def design_rows(ids): |
| return A.loc[[i for i in ids if i in A.index]] |
|
|
| |
| def fig1(): |
| fig = plt.figure(figsize=(12, 6.4)) |
| gs = gridspec.GridSpec(2, 5, height_ratios=[1, 1.15], hspace=0.05, wspace=0.12) |
| refs = ['ref1_nested_rings_wire_mesh', 'ref2_organic_thread_net', 'ref3_polygonal_cells_inner_venation', 'ref4_triangulated_lattice_defects', 'ref5_leaf_venation_hierarchy'] |
| caps = ['(1) nested rings:\nclosed loops, slender wires', '(2) thread net:\norganic disorder, redundancy', '(3) cells with inner veins:\nnested subdivision', '(4) triangulated lattice:\norder, braces, defects', '(5) leaf venation:\nthickness hierarchy, T-junctions'] |
| for i, (r, c) in enumerate(zip(refs, caps)): |
| ax = fig.add_subplot(gs[0, i]); ax.imshow(mpimg.imread(os.path.join(ROOT, 'reference_images', r + '_512.jpg'))); ax.set_xticks([]); ax.set_yticks([]); ax.set_title(c, fontsize=8.5) |
| specs = [('B2_ref_ord', 'ordered nested (δ = 0)'), ('E4_L3_dis0.35_rt0.6', 'level-selective (finest level organic)'), ('B2_ref', 'disordered (δ = 0.35)'), ('R_dis1', 'organic (δ = 1)'), ('E5_ord_red1', 'ordered + full bracing')] |
| for i, (d, t) in enumerate(specs): |
| ax = fig.add_subplot(gs[1, i]) |
| row = df[df.designId == d].iloc[0]; res, geom, snaps = load_run(row['dir']) |
| render(geom, ax, color_by='level', title=t, frame=True) |
| ax.set_xlabel('D=%d b=%d r_t=%.2g ρ=%.2f' % (row.p_depth, row.p_branching, row.p_thicknessRatio, row.s_relDensityDesign), fontsize=8) |
| handles = [Line2D([0], [0], color=c, lw=3, label='level %d' % k) for k, c in enumerate(['#d62728', '#ff7f0e', '#bcbd22', '#2ca02c'])] |
| fig.legend(handles=handles, loc='lower center', ncol=4, frameon=False, bbox_to_anchor=(0.5, 0.01)) |
| fig.suptitle('Reference images (top) and members of the generated hierarchical cellular-network family coloured by hierarchy level (bottom; line width ∝ member thickness)', fontsize=10, y=0.99) |
| savefig(fig, 'fig1_references_and_family') |
|
|
| |
| def fig2(): |
| fig, axes = plt.subplots(2, 4, figsize=(12, 5.6)); axes = axes.ravel() |
| pars = [('p_depth', 'hierarchy depth D'), ('p_branching', 'branching ratio b'), ('p_disorder', 'geometric disorder δ'), ('p_thicknessRatio', 'thickness ratio r_t'), |
| ('p_interfaceStrengthRatio', 'interface strength s_int'), ('p_redundancy', 'redundancy (braces)'), ('p_strengthCV', 'strength CV'), ('f0', 'level-0 material fraction f0')] |
| stage_col = {'baseline': '#7f7f7f', 'reconnaissance': '#1f77b4', 'hypothesis-tests': '#2ca02c', 'refinement': '#ff7f0e', 'holdout': '#d62728'} |
| for ax, (k, lab) in zip(axes, pars): |
| for st, col in stage_col.items(): |
| sub = df[df.stage == st] |
| if len(sub): ax.scatter(sub[k], sub['Wf_uJ'], s=12, alpha=0.7, color=col, label=st, edgecolor='none') |
| ax.set_xlabel(lab); ax.set_ylabel('work to failure W_f [10⁻⁶ E_s W² b_z]'); ax.set_yscale('log') |
| axes[0].legend(frameon=False, fontsize=7) |
| fig.suptitle('Coverage of the design space: every simulated run (n = %d runs, %d designs) by study stage; vertical axis = work to failure' % (len(df), df.designId.nunique()), fontsize=10) |
| fig.tight_layout(); savefig(fig, 'fig2_design_space_coverage') |
|
|
| |
| def pareto_front(x, y): |
| idx = np.argsort(-x); front = []; best = -np.inf |
| for i in idx: |
| if y[i] > best: front.append(i); best = y[i] |
| return np.array(front) |
|
|
| def fig3(): |
| fig, axes = plt.subplots(1, 3, figsize=(13, 4.3)) |
| B = A[A.stage != 'holdout'] if 'stage' in A else A |
| |
| ax = axes[0] |
| sc = ax.scatter(main['K0'], main['Wf_uJ'], c=main['abruptness'], cmap='coolwarm', vmin=0.15, vmax=1.0, s=10, alpha=0.5, edgecolor='none') |
| ax.scatter(B['K0_mean'], B['Wf_mean'] * 1e6, c=B['abruptness_mean'], cmap='coolwarm', vmin=0.15, vmax=1.0, s=34, edgecolor='k', linewidth=0.4) |
| fr = pareto_front(B['K0_mean'].values, B['Wf_mean'].values * 1e6) |
| o = np.argsort(B['K0_mean'].values[fr]); ax.plot(B['K0_mean'].values[fr][o], B['Wf_mean'].values[fr][o] * 1e6, 'k--', lw=1, label='Pareto front (max K0, max W_f)') |
| for i in fr: ax.annotate(B.index[i].replace('_', ' '), (B['K0_mean'].values[i], B['Wf_mean'].values[i] * 1e6), fontsize=6, xytext=(3, 3), textcoords='offset points') |
| ax.set_xlabel('initial stiffness K_0 [E_s b_z]'); ax.set_ylabel('work to failure W_f [10⁻⁶ E_s W² b_z]'); ax.set_yscale('log'); ax.legend(frameon=False, fontsize=7) |
| cb = fig.colorbar(sc, ax=ax, fraction=0.05, pad=0.02); cb.set_label('abruptness (max drop / F_peak)') |
| ax.set_title('(a) stiffness – work to failure', fontsize=9) |
| |
| ax = axes[1] |
| ax.scatter(main['abruptness'], main['Wf_uJ'], c=main['K0'], cmap='viridis', s=10, alpha=0.5, edgecolor='none') |
| sc2 = ax.scatter(B['abruptness_mean'], B['Wf_mean'] * 1e6, c=B['K0_mean'], cmap='viridis', s=34, edgecolor='k', linewidth=0.4) |
| fr = pareto_front(-B['abruptness_mean'].values, B['Wf_mean'].values * 1e6) |
| o = np.argsort(B['abruptness_mean'].values[fr]); ax.plot(B['abruptness_mean'].values[fr][o], B['Wf_mean'].values[fr][o] * 1e6, 'k--', lw=1, label='Pareto front (min abruptness, max W_f)') |
| for i in fr: ax.annotate(B.index[i].replace('_', ' '), (B['abruptness_mean'].values[i], B['Wf_mean'].values[i] * 1e6), fontsize=6, xytext=(3, 3), textcoords='offset points') |
| ax.set_xlabel('abruptness A (largest force drop / F_peak)'); ax.set_ylabel('work to failure W_f [10⁻⁶ E_s W² b_z]'); ax.set_yscale('log'); ax.legend(frameon=False, fontsize=7) |
| cb = fig.colorbar(sc2, ax=ax, fraction=0.05, pad=0.02); cb.set_label('K_0 [E_s b_z]') |
| ax.set_title('(b) damage tolerance – work to failure', fontsize=9) |
| |
| ax = axes[2] |
| Ue = 0.5 * main['Fpeak'] * main['dPeak'] * 1e6 |
| ax.scatter(Ue, main['Wf_uJ'], c=[MODE_COL.get(m, '#333') for m in main['fractureMode']], s=12, alpha=0.7, edgecolor='none') |
| lim = [min(Ue.min(), main['Wf_uJ'].min()) * 0.8, max(Ue.max(), main['Wf_uJ'].max()) * 1.3] |
| ax.plot(lim, lim, 'k-', lw=0.8, label='W_f = ½ F_peak d_peak'); ax.plot(lim, [3 * l for l in lim], 'k:', lw=0.8, label='3 × ') |
| ax.set_xscale('log'); ax.set_yscale('log'); ax.set_xlabel('elastic energy at peak ½ F_peak d_peak [10⁻⁶ E_s W² b_z]'); ax.set_ylabel('work to failure W_f [10⁻⁶ E_s W² b_z]') |
| ax.legend(handles=[Line2D([0], [0], marker='o', color='w', markerfacecolor=c, label=m, markersize=6) for m, c in MODE_COL.items() if m != 'no-failure'] + [Line2D([0], [0], color='k', lw=0.8, label='W_f = ½F_peak d_peak'), Line2D([0], [0], color='k', ls=':', lw=0.8, label='3×')], frameon=False, fontsize=7) |
| ax.set_title('(c) energy budget: post-peak contribution', fontsize=9) |
| fig.tight_layout(); savefig(fig, 'fig3_stiffness_toughness_pareto') |
|
|
| |
| def fig4(): |
| sel = [('B0_flat16_ord', 'flat ordered grid (D=0)', '#1f77b4'), ('B2_ref_ord', 'ordered nested hierarchy', '#9467bd'), ('E4_L3_dis0.35_rt0.6', 'level-selective (finest level organic)', '#2ca02c'), ('B2_ref', 'disordered hierarchy (δ=0.35)', '#d62728'), ('E5_ord_red1', 'ordered + full bracing', '#ff7f0e'), ('R_frame', 'disordered + side stringers', '#8c564b'), ('R_rt0.3', 'disordered, strong thickness hierarchy r_t=0.3', '#17becf')] |
| fig, axes = plt.subplots(1, 2, figsize=(12, 4.4)) |
| for d, lab, col in sel: |
| rows = df[df.designId == d] |
| for k, (_, r) in enumerate(rows.iterrows()): |
| res = json.load(open(os.path.join(ROOT, r['dir'], 'results.json'))) |
| dd = [s['d'] for s in res['steps']]; FF = [s['F'] * 1e3 for s in res['steps']] |
| for ax in axes: ax.plot(dd, FF, color=col, lw=1.3 if k == 0 else 0.6, alpha=1 if k == 0 else 0.45, label=lab + ' (n=%d seeds)' % len(rows) if k == 0 else None) |
| axes[0].set_xlim(0, 0.03); axes[0].set_ylim(0, None); axes[1].set_xlim(0, 0.03); axes[1].set_ylim(1e-3, 3); axes[1].set_yscale('log') |
| for ax in axes: ax.set_xlabel('applied displacement d [W]'); ax.set_ylabel('reaction force F [10⁻³ E_s W b_z]') |
| axes[0].legend(frameon=False, fontsize=7); axes[0].set_title('(a) linear scale', fontsize=9); axes[1].set_title('(b) logarithmic force scale (post-peak tails)', fontsize=9) |
| fig.tight_layout(); savefig(fig, 'fig4_load_displacement_curves') |
|
|
| |
| def fig5(): |
| sel = [('B2_ref_ord', 'ordered nested hierarchy: abrupt, localised'), ('B2_ref', 'disordered hierarchy: progressive, distributed'), ('E4_L3_dis0.35_rt0.6', 'level-selective: intermediate')] |
| keys = [('peak', 'at peak force (before the drop)'), ('f1', 'at 50 % of the final-failure displacement'), ('final', 'at final failure')] |
| fig, axes = plt.subplots(3, 3, figsize=(11, 11.2)) |
| for i, (d, lab) in enumerate(sel): |
| row = df[(df.designId == d) & (df.seed == 1)].iloc[0]; res, geom, snaps = load_run(row['dir']); m = res['metrics'] |
| for j, (k, kl) in enumerate(keys): |
| s = snaps[k]; ax = axes[i, j] |
| render(geom, ax, u=np.array(s['u']), scale=1.0, color_by='util', state=np.array(s['state']), util=np.array(s['util']), show_failed='ghost', frame=True) |
| nf = int((np.array(s['state']) == 2).sum()) |
| ax.set_title('%s\nd = %.4f, F = %.2e, failed %d' % (kl, s['d'], s['F'], nf), fontsize=7.5) |
| axes[i, 0].set_ylabel('%s\nA = %.2f, LI = %.2f, pre-peak failures %d' % (lab, m['abruptness'], m['localizationIndex'], m['prePeakFailures']), fontsize=8) |
| sm = plt.cm.ScalarMappable(cmap='viridis', norm=plt.Normalize(0, 1)); sm.set_array([]) |
| cb = fig.colorbar(sm, ax=axes, fraction=0.02, pad=0.02); cb.set_label('member utilisation ε_max/ε_c (failed members: dashed grey; deformation × 1)') |
| savefig(fig, 'fig5_fracture_snapshots') |
|
|
| |
| def sweep(ax, ids, xs, ykey, yscale=1, label=None, color=None, marker='o', ls='-'): |
| rows = design_rows(ids) |
| x = [xs[i] for i in rows.index]; y = rows[ykey + '_mean'].values * yscale; e = rows[ykey + '_std'].fillna(0).values * yscale |
| errbar(ax, x, y, e, fmt=marker + ls, color=color, label=label, ms=4, lw=1.1) |
| |
| for i in rows.index: |
| sub = df[df.designId == i]; ax.scatter([xs[i]] * len(sub), sub[ykey] * yscale, s=6, color=color, alpha=0.35, edgecolor='none') |
|
|
| def fig6(): |
| fig, axes = plt.subplots(4, 5, figsize=(15, 11.5)) |
| metrics = [('K0', 1, 'K_0 [E_s b_z]'), ('Fpeak', 1e3, 'F_peak [10⁻³ E_s W b_z]'), ('Wf', 1e6, 'W_f [10⁻⁶ E_s W² b_z]'), ('abruptness', 1, 'abruptness A')] |
| |
| dis_ids = {'R_dis0': 0, 'E1_dis0.05': 0.05, 'E1_dis0.1': 0.1, 'R_dis0.15': 0.15, 'E1_dis0.2': 0.2, 'B2_ref': 0.35, 'E1_dis0.5': 0.5, 'R_dis0.7': 0.7, 'R_dis1': 1.0} |
| cv_ids = {'R_cv0': 0, 'B2_ref': 0.05, 'R_cv0.2': 0.2, 'F6_dis_cv0.4': 0.4} |
| ocv_ids = {'B2_ref_ord': 0.05, 'E7_ord_cv0.1': 0.1, 'E7_ord_cv0.2': 0.2, 'E7_ord_cv0.4': 0.4} |
| for r, (mk, sc, lab) in enumerate(metrics): |
| ax = axes[r, 0]; sweep(ax, list(dis_ids), dis_ids, mk, sc, 'geometric disorder δ (r_t=0.6)', '#d62728') |
| ax.set_xlabel('geometric disorder δ'); ax.set_ylabel(lab) |
| if r == 0: ax.set_title('(a) geometric disorder', fontsize=9) |
| ax2 = ax.twiny(); sweep(ax2, list(ocv_ids), ocv_ids, mk, sc, 'strength CV, ordered geometry', '#1f77b4', 's', '--'); ax2.set_xlabel('strength CV (ordered geometry, δ=0)', color='#1f77b4', fontsize=8); ax2.tick_params(axis='x', colors='#1f77b4', labelsize=7) |
| if r == 0: ax.legend(frameon=False, fontsize=6.5, loc='center right'); ax2.legend(frameon=False, fontsize=6.5, loc='upper right') |
| |
| for r, (mk, sc, lab) in enumerate(metrics): |
| ax = axes[r, 1] |
| sub = A[(A.p_disorder == 0.35) & (A.dml == 0) & (A.p_redundancy == 0) & (A.p_interfaceStrengthRatio == 0.8) & (A.p_strengthCV == 0.05) & (A.p_relDensity == 0.2) & (A.p_layout == 'grid') & (~A.index.str.contains('notch|con|frame|H2|flat|noif|movie|_s3|movie|n0'))] |
| for D, col, mkr in [(0, '#7f7f7f', 'x'), (1, '#9467bd', 'v'), (2, '#17becf', '^'), (3, '#d62728', 'o'), (6, '#2ca02c', 's')]: |
| s2 = sub[sub.p_depth == D] |
| if not len(s2): continue |
| ax.errorbar(s2['f0'], s2[mk + '_mean'] * sc, yerr=s2[mk + '_std'].fillna(0) * sc, fmt=mkr, color=col, ms=4, capsize=2, lw=0.8, ls='none', label='D = %d' % D) |
| ax.set_xlabel('level-0 (coarsest) material fraction f_0'); ax.set_ylabel(lab) |
| if r == 0: ax.set_title('(b) thickness hierarchy (δ=0.35): collapse on f_0', fontsize=9); ax.legend(frameon=False, fontsize=6.5) |
| |
| sint_d = {'R_sint0.3': 0.3, 'R_sint0.5': 0.5, 'B2_ref': 0.8, 'R_sint1': 1.0, 'R_sint1.5': 1.5} |
| sint_o = {'E3_ord_sint0.3': 0.3, 'E3_ord_sint0.5': 0.5, 'B2_ref_ord': 0.8, 'E3_ord_sint1': 1.0, 'E3_ord_sint1.5': 1.5} |
| sint_l = {'F4_L3_sint0.5': 0.5, 'E4_L3_dis0.35_rt0.6': 0.8, 'F4_L3_sint1.5': 1.5} |
| for r, (mk, sc, lab) in enumerate(metrics): |
| ax = axes[r, 2]; sweep(ax, list(sint_o), sint_o, mk, sc, 'ordered (δ=0)', '#1f77b4', 's'); sweep(ax, list(sint_d), sint_d, mk, sc, 'disordered (δ=0.35)', '#d62728'); sweep(ax, list(sint_l), sint_l, mk, sc, 'level-selective L3', '#2ca02c', '^') |
| ax.set_xlabel('inter-level (interface) strength ratio s_int'); ax.set_ylabel(lab) |
| if r == 0: ax.set_title('(c) inter-level strength', fontsize=9); ax.legend(frameon=False, fontsize=6.5) |
| if mk in ('K0', 'Fpeak', 'Wf'): ax.set_yscale('log') |
| |
| red_o = {'B2_ref_ord': 0, 'E5_ord_red0.5': 0.5, 'E5_ord_red1': 1.0} |
| red_d = {'B2_ref': 0, 'E5_dis035_red0.25': 0.25, 'R_red0.5': 0.5, 'E5_dis035_red0.75': 0.75, 'R_red1': 1.0} |
| red_1 = {'R_dis1': 0, 'E5_dis1_red0.5': 0.5, 'E5_dis1_red1': 1.0} |
| red_l = {'E4_L3_dis0.35_rt0.6': 0, 'F3_L3_red0.5': 0.5, 'F3_L3_red1': 1.0} |
| for r, (mk, sc, lab) in enumerate(metrics): |
| ax = axes[r, 3]; sweep(ax, list(red_o), red_o, mk, sc, 'ordered (δ=0)', '#1f77b4', 's'); sweep(ax, list(red_d), red_d, mk, sc, 'δ=0.35', '#d62728'); sweep(ax, list(red_1), red_1, mk, sc, 'δ=1', '#8c564b', 'v'); sweep(ax, list(red_l), red_l, mk, sc, 'level-selective L3', '#2ca02c', '^') |
| ax.set_xlabel('redundancy (fraction of leaf-cell braces)'); ax.set_ylabel(lab) |
| if r == 0: ax.set_title('(d) redundancy (bracing)', fontsize=9); ax.legend(frameon=False, fontsize=6.5) |
| if mk in ('K0', 'Fpeak', 'Wf'): ax.set_yscale('log') |
| |
| depC = {'R_depthC1': 1, 'R_depthC2': 2, 'B2_ref': 3, 'R_depthC6': 6} |
| lsel = {'B2_ref': 0, 'E4_L1_dis0.35_rt0.6': 1, 'E4_L2_dis0.35_rt0.6': 2, 'E4_L3_dis0.35_rt0.6': 3, 'B2_ref_ord': 4} |
| lsel1 = {'E4_L1_dis0.35_rt1': 1, 'E4_L2_dis0.35_rt1': 2, 'E4_L3_dis0.35_rt1': 3} |
| lselD6 = {'F7_D6_L2': 2, 'F7_D6_L4': 4, 'F7_D6_L5': 5, 'F7_D6_L6': 6} |
| for r, (mk, sc, lab) in enumerate(metrics): |
| ax = axes[r, 4]; sweep(ax, list(lsel), lsel, mk, sc, 'levels ≥ k disordered (D=3, r_t=0.6); k=4 ≡ ordered', '#2ca02c', '^'); sweep(ax, list(lsel1), lsel1, mk, sc, 'same, r_t=1', '#bcbd22', 'v'); sweep(ax, list(lselD6), lselD6, mk, sc, 'D=6 (fixed contrast), levels ≥ k', '#17becf', 'd') |
| ax.set_xlabel('first disordered level k'); ax.set_ylabel(lab) |
| if r == 0: ax.set_title('(e) level-selective disorder', fontsize=9); ax.legend(frameon=False, fontsize=6.5) |
| if mk in ('K0', 'Fpeak', 'Wf'): ax.set_yscale('log') |
| fig.suptitle('Factor effects at fixed material (ρ = 0.2, 256-cell family). Markers: design means ± s.d. over seeds; small dots: individual seeds.', fontsize=10, y=0.995) |
| fig.tight_layout(); savefig(fig, 'fig6_factor_effects') |
| |
| fig, axes = plt.subplots(1, 4, figsize=(12, 3)) |
| for ax, (mk, sc, lab) in zip(axes, metrics): |
| sweep(ax, list(depC), depC, mk, sc, 'fixed contrast t_D/t_0 = 0.216', '#d62728'); dep = {'R_depth0': 0, 'R_depth1': 1, 'R_depth2': 2, 'B2_ref': 3, 'R_depth6': 6}; sweep(ax, list(dep), dep, mk, sc, 'fixed per-level ratio r_t = 0.6', '#1f77b4', 's') |
| ax.set_xlabel('hierarchy depth D (256 leaf cells)'); ax.set_ylabel(lab) |
| axes[0].legend(frameon=False, fontsize=7); fig.suptitle('Hierarchy depth at fixed finest resolution and fixed material (δ = 0.35)', fontsize=10) |
| fig.tight_layout(); savefig(fig, 'fig6b_depth_effect') |
|
|
| |
| def fig8(): |
| if 'ls_frac_ge0.8' not in df: print('no load-sharing data'); return |
| fig, axes = plt.subplots(1, 3, figsize=(13, 4.2)) |
| B = A.join(df.groupby('designId')[['ls_frac_ge0.8', 'lp_max_share', 'lp_top3_share', 'ls_frac_ge0.5']].mean()) |
| B = B[B.stage != 'holdout'] |
| ax = axes[0] |
| sc = ax.scatter(B['ls_frac_ge0.8'] + 1e-3, B['lp_max_share'], c=B['abruptness_mean'], cmap='coolwarm', vmin=0.15, vmax=1, s=30, edgecolor='k', linewidth=0.4) |
| ax.set_xscale('log'); ax.set_xlabel('fraction of members within 20 % of failure at the first event, φ_0.8 (+10⁻³)'); ax.set_ylabel('largest single-member share of the reaction force, s_max') |
| ax.axvline(0.05 + 1e-3, color='k', ls=':', lw=0.8); ax.axhline(0.30, color='k', ls=':', lw=0.8) |
| ax.text(0.06, 0.9, 'route (i): uniform\nload sharing', transform=ax.transAxes, fontsize=7.5); ax.text(0.06, 0.55, 'route (ii): dominant\ncoarse load path', transform=ax.transAxes, fontsize=7.5); ax.text(0.06, 0.08, 'progressive region', transform=ax.transAxes, fontsize=7.5) |
| cb = fig.colorbar(sc, ax=ax, fraction=0.05, pad=0.02); cb.set_label('abruptness A (design mean)'); ax.set_title('(a) two routes to abrupt failure', fontsize=9) |
| ax = axes[1] |
| for lab, mask, col, mk in [('ordered (φ_0.8 ≥ 0.05)', B['ls_frac_ge0.8'] >= 0.05, '#1f77b4', 's'), ('disordered / hierarchical', B['ls_frac_ge0.8'] < 0.05, '#d62728', 'o')]: |
| s2 = B[mask]; ax.errorbar(s2['lp_max_share'], s2['abruptness_mean'], yerr=s2['abruptness_std'].fillna(0), fmt=mk, color=col, ms=4, capsize=2, lw=0.7, ls='none', label=lab) |
| ax.set_xlabel('largest single-member force share s_max (initial state)'); ax.set_ylabel('abruptness A'); ax.legend(frameon=False, fontsize=7); ax.set_title('(b) abruptness vs dominant load path', fontsize=9) |
| ax = axes[2] |
| ax.errorbar(B['ls_frac_ge0.8'] + 1e-3, B['abruptness_mean'], yerr=B['abruptness_std'].fillna(0), fmt='o', color='#7f7f7f', ms=3, capsize=1.5, lw=0.6, ls='none') |
| ax.set_xscale('log'); ax.set_xlabel('φ_0.8 (+10⁻³)'); ax.set_ylabel('abruptness A'); ax.set_title('(c) abruptness vs load-sharing uniformity', fontsize=9) |
| fig.tight_layout(); savefig(fig, 'fig8_mechanism_two_routes') |
|
|
| |
| def fig7(): |
| P = json.load(open(os.path.join(ROOT, 'study', 'predictions.json')))['predictions'] |
| hold = df[df.stage == 'holdout'] |
| fig, axes = plt.subplots(1, 4, figsize=(14, 3.7)) |
| for ax, (mk, sc, lab, log) in zip(axes, [('K0', 1, 'K_0 [E_s b_z]', True), ('Fpeak', 1e3, 'F_peak [10⁻³ E_s W b_z]', True), ('Wf', 1e6, 'W_f [10⁻⁶ E_s W² b_z]', True), ('abruptness', 1, 'abruptness A', False)]): |
| xs, ys, xe, ye, ok = [], [], [], [], [] |
| for i, (k, v) in enumerate(P.items()): |
| pt, lo, hi = [q * sc for q in v[mk]] |
| obs = hold[hold.designId == k][mk] * sc |
| xs.append(pt); ys.append(obs.mean()); xe.append([pt - lo, hi - pt]); ye.append(obs.std() if len(obs) > 1 else 0); ok.append(lo <= obs.mean() <= hi) |
| ax.annotate(k.split('_')[0], (pt, obs.mean()), fontsize=6, xytext=(3, 2), textcoords='offset points') |
| ax.scatter([pt] * len(obs), obs, s=6, color='#7f7f7f', alpha=0.5, zorder=1) |
| xe = np.array(xe).T |
| ax.errorbar(xs, ys, xerr=xe, yerr=ye, fmt='o', ms=4, capsize=2, lw=0.8, color='#1f77b4', ecolor='#9ecae1', zorder=2) |
| for x, y, o in zip(xs, ys, ok): |
| if not o: ax.scatter([x], [y], s=60, facecolor='none', edgecolor='#d62728', linewidth=1.2, zorder=3) |
| lim = [min(min(xs), min(ys)) * 0.8, max(max(xs), max(ys)) * 1.2] |
| ax.plot(lim, lim, 'k--', lw=0.8) |
| if log: ax.set_xscale('log'); ax.set_yscale('log') |
| ax.set_xlabel('predicted (bar = stated range)'); ax.set_ylabel('observed (mean of 3 seeds ± s.d.)') |
| ax.set_title('%s: %d/%d within range' % (lab, sum(ok), len(ok)), fontsize=9) |
| fig.suptitle('Holdout designs (12 × 3 unseen seeds): predictions recorded before the runs vs observations; red circles = misses; grey dots = seeds', fontsize=9.5, y=1.02) |
| fig.tight_layout(); savefig(fig, 'fig7_holdout_predictions') |
|
|
| |
| def fig9(): |
| ids = {'R_rho0.1': 0.1, 'B2_ref': 0.2, 'R_rho0.3': 0.3} |
| rows = design_rows(list(ids)); rho = np.array([ids[i] for i in rows.index]) |
| fig, axes = plt.subplots(1, 3, figsize=(10, 3.2)) |
| for ax, (mk, sc, lab) in zip(axes, [('K0', 1, 'K_0 [E_s b_z]'), ('Fpeak', 1e3, 'F_peak [10⁻³ E_s W b_z]'), ('Wf', 1e6, 'W_f [10⁻⁶ E_s W² b_z]')]): |
| y = rows[mk + '_mean'].values * sc; e = rows[mk + '_std'].fillna(0).values * sc |
| ax.errorbar(rho, y, yerr=e, fmt='o', color='#d62728', capsize=2, ms=4) |
| p = np.polyfit(np.log(rho), np.log(y), 1); xx = np.linspace(0.08, 0.35, 20); ax.plot(xx, np.exp(p[1]) * xx ** p[0], 'k--', lw=0.8, label='∝ ρ^%.2f' % p[0]) |
| ax.set_xscale('log'); ax.set_yscale('log'); ax.set_xlabel('relative density ρ'); ax.set_ylabel(lab); ax.legend(frameon=False, fontsize=8) |
| fig.suptitle('Density scaling of the disordered reference (δ = 0.35): bending-dominated cellular-solid exponents', fontsize=9.5) |
| fig.tight_layout(); savefig(fig, 'fig9_density_scaling') |
|
|
| if __name__ == '__main__': |
| which = sys.argv[1:] or ['1', '2', '3', '4', '5', '6', '7', '8', '9'] |
| for w in which: |
| globals()['fig' + w]() |
| A.to_csv(os.path.join(PROC, 'design_summary.csv')) |
|
|
