#!/usr/bin/env python3 """Tier-2 item 15 -- extend the epistemic design with the excluded inputs. Why this run matters -------------------- The six-input box covers chromatin modulus scale, shell Young modulus, shell shrinkage coupling, pretension, natural volume ratio and probe stress. Two things that plainly belong in an epistemic design are missing, and one of them is a config field that could have been swept at no extra cost: * `shell_thickness_um` = 0.12 um, classified in Table S1 as "model input / assumption" with no literature provenance. The reported endpoint is N = sigma * t_current, so it is DIRECTLY proportional to the thickness that was never varied, while inputs with far weaker leverage were. * partial core-shell interface coupling. Real chromatin-lamina attachment is partial and discrete (LAD contacts); the model ties core to shell perfectly. The design instead varies "shell shrinkage coupling", which is the more arbitrary of the two: it is an eigenstrain exponent fraction (nonlinear_solver.py:1849, `alpha * cfg.stimulus.shell_eigenstrain_fraction`), i.e. how much of the core's logarithmic natural shrinkage the shell adopts -- not a mechanical attachment at all. How much does the missing thickness matter? Take the published nominal endpoints (control 0.060352, hyperosmotic 0.037475 mN/m) and the fact that the prestress contributes S0_t * t0 = pretension = 0.05 mN/m to the reported resultant independently of t0 (prestress_pa = 1000 * pretension / t0, nonlinear_solver.py: 1043, and the resultant multiplies by the current thickness, lines 1562-1578, so the two factors cancel). The elastic part is what scales with thickness: control(t) = 0.05 + 0.010352 * (t / 0.12) difference(t) = -0.022877 * (t / 0.12) percent(t) = 100 * difference(t) / control(t) t = 0.06 um -> about -20.7 % t = 0.12 um -> -37.9 % (the published headline) t = 0.25 um -> about -66.6 % That swing is as large as the strongest input in the published design (shell_young_pa spans -17.0 % to -74.6 %; pretension spans -65.9 % to -20.5 %), so the headline percentage rests on an unswept assumption with top-rank leverage. The estimate above is first order -- it holds the elastic solution fixed, which is not exact because the prestress magnitude and the chromatin inner_axes both depend on t0 -- and this script prints it next to the solved value so the approximation is testable rather than asserted. Why BOTH absolute and percent are reported ------------------------------------------ N = sigma * t is linear in thickness, so it is tempting to assume thickness cancels from a relative change. It does not. The reported resultant uses the CURRENT thickness, t0 * lambda_n (nonlinear_solver.py:1562-1563), and the two conditions sit at different thickness stretches, so t0 is weighted differently in the two endpoints. On top of that, t0 re-enters the mechanics twice: the shell prestress is inversely proportional to it (lines 1043, 1416-1419, 1526-1529) and the chromatin map is sampled on inner_axes = outer_axes - t0 (lines 2247-2249). The absolute difference is very nearly proportional to t0; the percent change is a hyperbolic function of it, because the t0-invariant prestress offset sits in the denominator. Reporting only one of the two hides half the sensitivity. What this script does NOT do ---------------------------- It does not regenerate the mesh. `shell_thickness_um` is a config field and the solver reads it for the prestress, the chromatin sampling and the reported resultant, but the core/shell partition is baked into the gmsh mesh through cell tags, and nothing in the solver checks the config against it (verified: no mesh-vs-thickness consistency test exists). So this is a sweep of the REPORTED and PRESTRESSED thickness on fixed geometry, which is the honest description and is exactly what the Table S1 entry feeds. The fully consistent experiment regenerates `ellipsoid_core_shell.msh` per thickness so the shell's element geometry moves too; the mesh generator is not shipped in the package, so that variant cannot be scripted from here. The dry-run prints this caveat. Core-shell interface coupling: NO SUCH KNOB EXISTS -------------------------------------------------- Read before believing any script that claims to sweep it. In `nuclear_envelope_fem/nonlinear_solver.py` (2298 lines): * `_build_material_regions` (line 2270) builds core and shell as ELEMENT SUBSETS of one global basis: `basis.with_elements(core_elements)` and `basis.with_elements(shell_elements)`. Both regions share the same displacement dofs, so the nodes on the core-shell interface are literally the same unknowns. Continuity there is exact, unconditional, and not parameterised by anything. * `_assemble_system` (line 1024) loops over volume regions only. The single surface assembly in the file is `_assemble_outer_nominal_shear_load` (line 2028) on `tags.OUTER_NE_FACET` -- the probe. There is no interface facet integral. * `_region_forms` (line 1064) contains the Neo-Hookean residual/tangent and the prestress term; no cohesive, penalty or spring contribution. * `NonlinearSolveOptions` (lines 49-50, a frozen dataclass) has 31 fields and not one of them concerns the interface. The only cell tags consumed anywhere are `tags.CORE_VOLUME`, `tags.SHELL_VOLUME` (line 589) and `tags.OUTER_NE_FACET` (line 608); the word "interface" occurs once, in a mesh-tag validation message (line 605). No cohesive/tie/contact/attachment/LAD term exists. So the honest answer is: partial interface coupling cannot be turned on from a config or an option. Introducing it is a code change, minimally: 1. mesh: duplicate the nodes on the core-shell interface so core and shell have independent displacements there, and emit an interface facet tag. This is a change to the (unshipped) gmsh generator plus `nuclear_envelope_fem/tags.py`. 2. options: add `interface_coupling_fraction` (1.0 = today's perfect tie) and `interface_stiffness_pa_per_um` to `NonlinearSolveOptions` (line 50). 3. solver: add `_assemble_interface_coupling(...)` producing the residual and tangent of a traction t = k * f * [[u]] on an interface `FacetBasis`, and call it from `_assemble_system` (line 1024) beside the volume loop; the tangent contribution is mandatory or Newton loses quadratic convergence. 4. warm start: `_project_p1_nodal_warm_start` (line 1131) asserts one displacement per mesh vertex; duplicated interface nodes break that assumption and every stored `.npy` warm start becomes stale. 5. post-processing: the outer-surface proxy is unaffected, but `_deformed_geometry_summary` and the compression-area diagnostics should be re-checked for the now-discontinuous field. Step 1 alone changes the mesh, which invalidates every published equilibrium. That is a real result for the manuscript: the model cannot represent partial chromatin-lamina attachment, and the sentence about LADs should say so instead of implying `shell_coupling` covers it. Usage: python3 15_extended_parameter_box.py --project-root /path/to/project \\ --out-root /path/to/out/extended --mode oat --dry-run python3 15_extended_parameter_box.py --project-root /path/to/project \\ --out-root /path/to/out/extended --mode lhs --points 64 --dry-run """ from __future__ import annotations import argparse import copy import csv import json import math import sys from pathlib import Path import numpy as np sys.path.insert(0, str(Path(__file__).resolve().parents[1] / "tier1_scripts")) from _common import ( # noqa: E402 add_common_arguments, import_runner, patch_options, ) PARAMETERS = ( "chromatin_modulus_scale", "shell_young_pa", "shell_coupling", "pretension_mn_per_m", "natural_volume_ratio", "probe_shear_stress_pa", ) THICKNESS = "shell_thickness_um" EXTENDED_PARAMETERS = PARAMETERS + (THICKNESS,) CHANGE_COLUMN = "paired_mean_resultant_change_percent" DIFFERENCE_COLUMN = "paired_mean_resultant_difference_mn_per_m" CONTROL_COLUMN = "control_mean_resultant_mn_per_m" # Published design ranges (reynolds_finan_epistemic_uq.json) and the nominal # anchor (reynolds_finan_boundary_stress.json), reproduced so --dry-run works # without the project tree. Validated against those files at run time. RANGES = { "chromatin_modulus_scale": (0.5, 2.0, "log"), "shell_young_pa": (1000.0, 4000.0, "log"), "shell_coupling": (0.0, 0.7, "linear"), "pretension_mn_per_m": (0.025, 0.1, "log"), "natural_volume_ratio": (0.65, 0.85, "linear"), "probe_shear_stress_pa": (250.0, 500.0, "linear"), } NOMINAL = { "chromatin_modulus_scale": 1.0, "shell_young_pa": 2000.0, "shell_coupling": 0.35, "pretension_mn_per_m": 0.05, "natural_volume_ratio": 0.7255172413793103, "probe_shear_stress_pa": 375.0, } NOMINAL_THICKNESS_UM = 0.12 PUBLISHED_ANCHOR = { "control_mean_resultant_mn_per_m": 0.06035217853263504, "hyperosmotic_mean_resultant_mn_per_m": 0.03747527607923647, "paired_mean_resultant_difference_mn_per_m": -0.022876902453398566, "paired_mean_resultant_change_percent": -37.905677988124374, "discretization": "coarse_h1p25_P1", } PUBLISHED_DESIGN_NAMES = ( "epistemic_uq_design.csv", "FigureS5_UQ_design_source_data.csv", ) def main() -> int: parser = argparse.ArgumentParser(description=__doc__) add_common_arguments(parser) # Table S4/S5 rows are coarse_h1p25 P1. Merging a refined thickness sweep # into them would confound thickness with the mesh bias that item 13 measures. parser.set_defaults(primary_mesh="coarse_h1p25") parser.add_argument( "--mode", choices=("oat", "lhs"), default="oat", help="oat: thickness ladder around the nominal anchor. lhs: fresh 7-D design.", ) parser.add_argument("--thickness-min", type=float, default=0.06) parser.add_argument("--thickness-max", type=float, default=0.25) parser.add_argument( "--thickness-um", type=float, nargs="*", default=None, help="Explicit oat thickness levels (um), replacing the default ladder.", ) parser.add_argument( "--nominal-thickness-um", type=float, default=NOMINAL_THICKNESS_UM, help="Published thickness; the oat anchor and the reference for ratios.", ) parser.add_argument("--points", type=int, default=64, help="lhs sample count.") parser.add_argument( "--seed", type=int, default=20260810, help=( "lhs seed. The published seed reproduces the published six columns " "exactly and only adds the thickness coordinate; see --verify-design." ), ) parser.add_argument( "--verify-design", type=Path, default=None, help=( "Published design CSV to check the lhs draw against. Default: search " f"the tree for {', '.join(PUBLISHED_DESIGN_NAMES)}." ), ) parser.add_argument( "--element-order", choices=("p1", "p2"), default="p1", help="p1 matches the published screening discretization.", ) parser.add_argument("--warm-start-root", type=Path, default=None) parser.add_argument("--cuda-device-index", type=int, default=0) args = parser.parse_args() if not 0.0 < args.thickness_min < args.thickness_max: raise SystemExit("--thickness-min must be positive and below --thickness-max") if args.points < 2: raise SystemExit("--points must be at least 2") if args.mode == "oat": design = oat_design(args) else: design = lhs_design( args.points, args.seed, thickness_range=(args.thickness_min, args.thickness_max), ) degree = 2 if args.element_order == "p2" else 1 discretization = f"{args.primary_mesh}_P{degree}" print(f"Extended parameter box ({args.mode}) -- planned runs") print(f" project root : {args.project_root}") print(f" output root : {args.out_root}") print(f" discretization: {discretization}") print( f" new input : {THICKNESS} " f"{args.thickness_min:g}-{args.thickness_max:g} um (log), " f"published value {args.nominal_thickness_um:g}" ) print(f" design points : {len(design)} pairs = {2 * len(design)} equilibria") _print_design(design) if args.mode == "lhs": _report_design_verification(design, args) _print_first_order_prediction(design, args) _print_scope_notes() if args.dry_run: print("\n--dry-run: nothing solved.") return 0 uq = _import_uq(args.project_root) if tuple(uq.PARAMETERS) != PARAMETERS: raise SystemExit( f"Upstream PARAMETERS order changed ({uq.PARAMETERS}); " "this script needs updating." ) if degree == 2: # Upstream _options couples quadrature order and load steps to the degree. patch_options(uq, element_degree=2, quadrature_order=3, initial_load_steps=16) base = json.loads(uq.BASE_CONFIG.read_text(encoding="utf-8")) _check_thickness_field(base) _preflight_thickness_propagation(uq, base, design[0][THICKNESS]) control_nodal, hyper_nodal, nominal_volume_ratio = _load_warm_start( _resolve_warm_start_root(args, uq) ) args.out_root.mkdir(parents=True, exist_ok=True) design_path = args.out_root / f"extended_box_{args.mode}_design.csv" result_path = args.out_root / f"extended_box_{args.mode}_results.csv" failure_path = args.out_root / f"extended_box_{args.mode}_failures.csv" _write_csv(design_path, design) element_label = f"vector_P{degree}_tetrahedron" rows: list[dict] = [] failures: list[dict] = [] for point in design: thickness = float(point[THICKNESS]) print(f"\n=== {point['sample']} (t = {thickness:g} um) ===", flush=True) point_base = copy.deepcopy(base) point_base["geometry"][THICKNESS] = thickness try: row = uq._run_sample( args.out_root, point_base, dict(point), resume=args.resume, sparse_backend=args.sparse_backend, cuda_device_index=args.cuda_device_index, rigid_constraint_linear_solver=args.rigid_constraint_linear_solver, nominal_control_nodal=control_nodal, nominal_hyperosmotic_nodal=hyper_nodal, nominal_volume_ratio=nominal_volume_ratio, mesh_name=args.primary_mesh, ) except Exception as exc: # noqa: BLE001 -- never silently drop a point print(f" FAILED: {type(exc).__name__}: {exc}") failures.append( { "sample": point["sample"], THICKNESS: thickness, "error_type": type(exc).__name__, "error_message": str(exc), "eligible_for_summary": False, } ) _write_csv(failure_path, failures) continue # _run_sample hard-codes the P1 label and the coarse-screen interpretation. row["element"] = element_label row["mesh"] = args.primary_mesh rows.append(row) _write_csv(result_path, rows) print( f" difference {float(row[DIFFERENCE_COLUMN]):+.6f} mN/m, " f"change {float(row[CHANGE_COLUMN]):+.2f} %", flush=True, ) if not rows: raise SystemExit("No design point converged; nothing to summarise.") rows = _add_thickness_diagnostics(rows, args) _write_csv(result_path, rows) summary = summarize(rows, failures, design, discretization, args) summary_path = args.out_root / f"extended_box_{args.mode}_summary.json" summary_path.write_text( json.dumps(summary, indent=2, ensure_ascii=False) + "\n", encoding="utf-8" ) _print_summary(summary, failures) print(f"\nWrote {design_path}") print(f"Wrote {result_path}") if failures: print(f"Wrote {failure_path}") print(f"Wrote {summary_path}") return 0 def oat_design(args: argparse.Namespace) -> list[dict]: """Anchor plus a thickness ladder, all other inputs at the nominal point.""" anchor = float(args.nominal_thickness_um) if args.thickness_um: levels = [float(value) for value in args.thickness_um] else: # Low bound, high bound and the log-midpoint of anchor and high bound; # the log-midpoint of the full range coincides with the anchor, so it # would carry no information. levels = [ args.thickness_min, math.sqrt(anchor * args.thickness_max), args.thickness_max, ] design = [_oat_point(anchor, args)] for level in levels: if level <= 0.0: raise SystemExit("--thickness-um values must be positive") if abs(math.log(level / anchor)) < 1.0e-3: continue design.append(_oat_point(level, args)) if len({point["sample"] for point in design}) != len(design): raise SystemExit("Duplicate thickness levels requested") return design def _oat_point(thickness: float, args: argparse.Namespace) -> dict: anchor = float(args.nominal_thickness_um) if abs(math.log(thickness / anchor)) < 1.0e-3: design_class, signature = "nominal_anchor", "nominal" elif thickness <= args.thickness_min: design_class, signature = "one_at_a_time_boundary", f"{THICKNESS}:low" elif thickness >= args.thickness_max: design_class, signature = "one_at_a_time_boundary", f"{THICKNESS}:high" else: design_class, signature = "one_at_a_time_interior", f"{THICKNESS}:interior" return { "sample": f"ext_t{thickness:.4f}".replace(".", "p"), "design_class": design_class, "boundary_signature": signature, "boundary_dimension_count": 0 if design_class == "nominal_anchor" else 1, **{name: float(NOMINAL[name]) for name in PARAMETERS}, THICKNESS: float(thickness), } def lhs_design( points: int, seed: int, *, thickness_range: tuple[float, float] ) -> list[dict]: """Fixed-seed LHS over the seven inputs. Reimplements `run_reynolds_finan_epistemic_uq.latin_hypercube_design` so the design can be printed without the project tree (that module needs scipy and the FEM package). The generator draws one permutation and one uniform vector per column in order, so appending thickness as the seventh column leaves the first six columns bit-identical to the published 64-point design at the same seed and count -- verified by --verify-design. """ spec = dict(RANGES) spec[THICKNESS] = (thickness_range[0], thickness_range[1], "log") rng = np.random.default_rng(int(seed)) unit = np.empty((points, len(EXTENDED_PARAMETERS)), dtype=float) for column in range(len(EXTENDED_PARAMETERS)): unit[:, column] = (rng.permutation(points) + rng.random(points)) / points design = [] for index in range(points): row = { "sample": f"ext_lhs_{index + 1:03d}", "design_class": "extended_latin_hypercube", "boundary_signature": "", "boundary_dimension_count": len(EXTENDED_PARAMETERS), "published_sample_id": f"lhs_{index + 1:03d}", } for column, parameter in enumerate(EXTENDED_PARAMETERS): row[parameter] = _transform_unit_interval(unit[index, column], spec[parameter]) design.append(row) return design def _transform_unit_interval(value: float, spec: tuple[float, float, str]) -> float: lower, upper, scale = spec if scale == "linear": return lower + value * (upper - lower) if scale == "log": return math.exp(math.log(lower) + value * (math.log(upper) - math.log(lower))) raise ValueError(f"Unsupported parameter scale: {scale}") def summarize( rows: list[dict], failures: list[dict], design: list[dict], discretization: str, args: argparse.Namespace, ) -> dict: changes = np.asarray([float(row[CHANGE_COLUMN]) for row in rows]) differences = np.asarray([float(row[DIFFERENCE_COLUMN]) for row in rows]) thickness = np.asarray([float(row[THICKNESS]) for row in rows]) anchor = _anchor_row(rows, args) summary = { "analysis": "reynolds_finan_extended_parameter_box", "mode": args.mode, "discretization": discretization, "planned_pairs": len(design), "converged_pairs": len(rows), "failed_pairs": len(failures), "failed_samples": [failure["sample"] for failure in failures], "new_input": { THICKNESS: { "min": args.thickness_min, "max": args.thickness_max, "scale": "log", "published_value": args.nominal_thickness_um, "basis": ( "Table S1 lists 0.12 um as a model input with no literature " "provenance; the reported resultant is proportional to it." ), } }, "excluded_input_requiring_code_change": { "name": "partial_core_shell_interface_coupling", "status": "not_implementable_from_configuration", "evidence": ( "nonlinear_solver.py:2270 _build_material_regions ties core and " "shell through a shared basis; nonlinear_solver.py:1024 " "_assemble_system assembles volume regions only; the sole facet " "assembly is _assemble_outer_nominal_shear_load at line 2028; " "NonlinearSolveOptions (line 50) has no interface field." ), "minimal_change": ( "duplicate interface nodes and tag the interface facets in the " "mesh generator, add interface_coupling_fraction and " "interface_stiffness_pa_per_um to NonlinearSolveOptions, add " "_assemble_interface_coupling with its tangent and call it from " "_assemble_system, then regenerate every warm start because " "_project_p1_nodal_warm_start (line 1131) assumes one " "displacement per mesh vertex." ), }, "paired_mean_resultant_change_percent": { "minimum": float(np.min(changes)), "median": float(np.median(changes)), "maximum": float(np.max(changes)), }, "paired_mean_resultant_difference_mn_per_m": { "minimum": float(np.min(differences)), "median": float(np.median(differences)), "maximum": float(np.max(differences)), }, "all_pairs_unload": bool(np.all(changes < 0.0)), "thickness_um_span": [float(np.min(thickness)), float(np.max(thickness))], "geometry_caveat": ( "shell_thickness_um was varied in the case configuration only. The " "mesh, and therefore the shell's element geometry, is unchanged; the " "sweep covers the prestress magnitude, the chromatin inner_axes and " "the reported resultant scale." ), "interpretation": ( "Deterministic coverage of an input the published design excluded. " "Not a probability distribution, and it does not promote local field " "or compression endpoints." ), } if anchor is not None: summary["anchor"] = { "sample": anchor["sample"], THICKNESS: float(anchor[THICKNESS]), CHANGE_COLUMN: float(anchor[CHANGE_COLUMN]), DIFFERENCE_COLUMN: float(anchor[DIFFERENCE_COLUMN]), "published_change_percent": PUBLISHED_ANCHOR[CHANGE_COLUMN], "published_discretization": PUBLISHED_ANCHOR["discretization"], "reproduces_published_anchor_within_0p01_pp": bool( abs(float(anchor[CHANGE_COLUMN]) - PUBLISHED_ANCHOR[CHANGE_COLUMN]) < 0.01 ), "note": ( "Only comparable when this run used coarse_h1p25 P1; a mismatch " "there means the mesh differs, not that the thickness plumbing " "misfired." ), } residuals = [ row["linear_in_thickness_residual_percent"] for row in rows if row.get("linear_in_thickness_residual_percent") is not None ] if residuals: summary["absolute_difference_linear_in_thickness"] = { "maximum_absolute_residual_percent": float( np.max(np.abs(np.asarray(residuals))) ), "claim": ( "difference(t) / difference(anchor) equals t / t_anchor to " "within this residual, i.e. the absolute endpoint is linear " "in thickness while the percent change is not." ), } if args.mode == "lhs" and len(rows) >= 3: summary["descriptive_rank_correlation"] = [ { "parameter": parameter, "spearman_rho": _spearman_rho( np.asarray([float(row[parameter]) for row in rows]), changes ), "interpretation": ( "descriptive_parameter_ranking_not_biological_inference; " "no p value is computed here (numpy only)" ), } for parameter in EXTENDED_PARAMETERS ] summary["descriptive_rank_correlation"].sort( key=lambda entry: abs(entry["spearman_rho"]), reverse=True ) for rank, entry in enumerate(summary["descriptive_rank_correlation"], start=1): entry["rank"] = rank return summary def _add_thickness_diagnostics(rows: list[dict], args: argparse.Namespace) -> list[dict]: """Attach anchor-relative and first-order-prediction columns.""" anchor = _anchor_row(rows, args) for row in rows: thickness = float(row[THICKNESS]) difference = float(row[DIFFERENCE_COLUMN]) row["thickness_ratio_to_anchor"] = thickness / float(args.nominal_thickness_um) row["absolute_difference_mn_per_m"] = difference if anchor is None or row is anchor: row["difference_ratio_to_anchor"] = 1.0 if anchor is not None else None row["linear_in_thickness_residual_percent"] = 0.0 if anchor is not None else None row["first_order_predicted_change_percent"] = None row["change_minus_first_order_prediction_percentage_points"] = None continue anchor_difference = float(anchor[DIFFERENCE_COLUMN]) ratio = difference / anchor_difference if anchor_difference else float("nan") row["difference_ratio_to_anchor"] = ratio row["linear_in_thickness_residual_percent"] = 100.0 * ( ratio / row["thickness_ratio_to_anchor"] - 1.0 ) predicted = _first_order_change_percent( float(anchor[CONTROL_COLUMN]), anchor_difference, float(row["pretension_mn_per_m"]), row["thickness_ratio_to_anchor"], ) row["first_order_predicted_change_percent"] = predicted row["change_minus_first_order_prediction_percentage_points"] = ( float(row[CHANGE_COLUMN]) - predicted if predicted is not None else None ) return rows def _first_order_change_percent( anchor_control: float, anchor_difference: float, pretension: float, ratio: float ) -> float | None: """Percent change if only the elastic part scaled with thickness. The prestress contributes S0_t * t0 = pretension to the reported resultant whatever the thickness, so it stays in the denominator while everything else scales. """ elastic_control = anchor_control - pretension control = pretension + elastic_control * ratio if abs(control) < 1.0e-12: return None return 100.0 * anchor_difference * ratio / control def _anchor_row(rows: list[dict], args: argparse.Namespace) -> dict | None: anchor = float(args.nominal_thickness_um) for row in rows: if abs(math.log(float(row[THICKNESS]) / anchor)) < 1.0e-3 and all( abs(float(row[name]) - NOMINAL[name]) <= 1.0e-9 * max(1.0, abs(NOMINAL[name])) for name in PARAMETERS ): return row return None def _spearman_rho(values: np.ndarray, outcome: np.ndarray) -> float: """Spearman rho with average ranks; numpy only, no p value (no scipy here).""" ranked = [_average_ranks(values), _average_ranks(outcome)] first, second = (array - array.mean() for array in ranked) denominator = math.sqrt(float(np.sum(first**2)) * float(np.sum(second**2))) if denominator == 0.0: return float("nan") return float(np.sum(first * second) / denominator) def _average_ranks(values: np.ndarray) -> np.ndarray: order = np.argsort(values, kind="mergesort") ranks = np.empty(len(values), dtype=float) ranks[order] = np.arange(1, len(values) + 1, dtype=float) unique, inverse, counts = np.unique(values, return_inverse=True, return_counts=True) if len(unique) != len(values): sums = np.zeros(len(unique), dtype=float) np.add.at(sums, inverse, ranks) ranks = (sums / counts)[inverse] return ranks def _print_design(design: list[dict]) -> None: columns = EXTENDED_PARAMETERS header = ( f"\n {'sample':<18}{'design_class':<26}" + "".join(f"{name[:11]:>13}" for name in columns) ) print(header) print(" " + "-" * (len(header) - 3)) for point in design: line = f" {point['sample']:<18}{point['design_class']:<26}" line += "".join(f"{float(point[name]):>13.5g}" for name in columns) print(line) def _print_first_order_prediction(design: list[dict], args: argparse.Namespace) -> None: print( "\n First-order prediction from the published anchor " f"({PUBLISHED_ANCHOR['discretization']}, control " f"{PUBLISHED_ANCHOR[CONTROL_COLUMN]:.6f} mN/m, difference " f"{PUBLISHED_ANCHOR[DIFFERENCE_COLUMN]:.6f} mN/m):" ) print(f" {'thickness um':>14}{'difference':>14}{'change %':>12}") levels = sorted({round(float(point[THICKNESS]), 6) for point in design}) if len(levels) > 8: levels = [levels[0], levels[len(levels) // 2], levels[-1]] for thickness in levels: ratio = thickness / float(args.nominal_thickness_um) predicted = _first_order_change_percent( PUBLISHED_ANCHOR[CONTROL_COLUMN], PUBLISHED_ANCHOR[DIFFERENCE_COLUMN], NOMINAL["pretension_mn_per_m"], ratio, ) difference = PUBLISHED_ANCHOR[DIFFERENCE_COLUMN] * ratio predicted_text = "n/a" if predicted is None else f"{predicted:.2f}" print(f" {thickness:>14.4f}{difference:>14.6f}{predicted_text:>12}") if len(levels) < len({round(float(point[THICKNESS]), 6) for point in design}): print(" (min / median / max of the sampled thicknesses)") print( " Holds the elastic solution fixed, so it is falsifiable, not a result: " "the solve moves it through the prestress and the chromatin inner_axes." ) def _report_design_verification(design: list[dict], args: argparse.Namespace) -> None: path = args.verify_design or _find_file(args.project_root, PUBLISHED_DESIGN_NAMES) if path is None: print( "\n Published design CSV not found; skipping the reproduction check. " "Pass --verify-design to enable it." ) return with Path(path).open(newline="", encoding="utf-8") as handle: published = list(csv.DictReader(handle)) if len(published) != len(design): print( f"\n Published design has {len(published)} rows against " f"{len(design)} planned; reproduction check skipped " "(use --points to match)." ) return deviation = max( abs(float(design[index][name]) - float(published[index][name])) for index in range(len(design)) for name in PARAMETERS ) print(f"\n Reproduction check against {path}") print(f" maximum deviation over the six published columns: {deviation:.3e}") if deviation == 0.0: print( " The extended design reproduces the published points exactly and " "only adds the thickness coordinate, so the new rows are directly " "comparable to Table S4 row by row." ) else: print( " NOT a superset of the published design -- seed or count differs. " "Report the extended design as a separate sample, not as an " "augmentation of Table S4." ) def _print_scope_notes() -> None: print("\n Scope notes") print( " * shell_thickness_um is swept in the case config only. The mesh is not " "regenerated, so the shell's element geometry is unchanged; the sweep " "moves the prestress magnitude (nonlinear_solver.py:1043), the chromatin " "inner_axes (line 2247) and the reported resultant scale (line 1562)." ) print( " * partial core-shell interface coupling is NOT implemented here " "because no such knob exists: core and shell share displacement dofs " "through one basis (_build_material_regions, nonlinear_solver.py:2270), " "_assemble_system (line 1024) has no interface facet integral, and " "NonlinearSolveOptions (line 50) has no interface field. Adding it needs " "duplicated interface nodes in the mesh, two new option fields, a new " "_assemble_interface_coupling residual/tangent pair called from " "_assemble_system, and regenerated warm starts " "(_project_p1_nodal_warm_start, line 1131). That invalidates every " "published equilibrium, so it is a manuscript statement, not a run." ) def _print_summary(summary: dict, failures: list[dict]) -> None: change = summary["paired_mean_resultant_change_percent"] difference = summary["paired_mean_resultant_difference_mn_per_m"] print("\n=== Extended box summary ===") print( f" change % : min {change['minimum']:+.2f}, " f"median {change['median']:+.2f}, max {change['maximum']:+.2f}" ) print( f" difference : min {difference['minimum']:+.6f}, " f"median {difference['median']:+.6f}, max {difference['maximum']:+.6f} mN/m" ) linear = summary.get("absolute_difference_linear_in_thickness") if linear: print( " linearity : absolute difference tracks t / t_anchor to within " f"{linear['maximum_absolute_residual_percent']:.2f} %" ) anchor = summary.get("anchor") if anchor: print( f" anchor check : {anchor[CHANGE_COLUMN]:+.4f} % against published " f"{anchor['published_change_percent']:+.4f} % " f"({'match' if anchor['reproduces_published_anchor_within_0p01_pp'] else 'DIFFERS'})" ) if failures: print("\n=== Points that did NOT converge ===") for failure in failures: print( f" {failure['sample']:<18}{failure['error_type']}: " f"{failure['error_message']}" ) def _import_uq(project_root: Path): try: import_runner(_source_root(project_root)) except ImportError as exc: raise SystemExit( f"The runner under {project_root} could not be imported: {exc}. " "The submission package is incomplete -- it ships no src/, no meshes " "and no chromatin .npz. Point --project-root at the real working tree." ) from exc try: import run_reynolds_finan_epistemic_uq as uq # noqa: PLC0415 except ImportError as exc: raise SystemExit( f"Could not import run_reynolds_finan_epistemic_uq: {exc}. " "It needs src/nuclear_envelope_fem, run_extreme_finite_deformation, " "run_vaziri_finan_two_group and scipy on the real project tree." ) from exc return uq def _source_root(project_root: Path) -> Path: for candidate in ( project_root / "scripts", project_root / "05_Source_Data", project_root, ): if (candidate / "run_reynolds_finan_heterogeneous_two_group.py").exists(): return candidate return project_root def _check_thickness_field(base: dict) -> None: if THICKNESS not in base.get("geometry", {}): raise SystemExit( "The base case config has no geometry.shell_thickness_um, so this " "sweep would silently do nothing. The solver reads that field at " "nonlinear_solver.py:1043 (prestress), 2247 (chromatin inner_axes) " "and 1562 (reported resultant); find where it now lives and update " "this script rather than running it." ) def _preflight_thickness_propagation(uq, base: dict, thickness: float) -> None: """Fail loudly if the template builder drops the thickness override.""" probe = copy.deepcopy(base) probe["geometry"][THICKNESS] = thickness try: control, hyper = uq.build_template_case_configs( probe, reference_osmolarity_mosm=380.0, hyperosmotic_osmolarity_mosm=580.0, solid_fraction=0.204, core_young_pa=800.0, shell_young_pa=float(NOMINAL["shell_young_pa"]), shell_coupling=float(NOMINAL["shell_coupling"]), ) except Exception as exc: # noqa: BLE001 print(f" note: could not pre-verify thickness propagation ({exc}).") return for label, config in (("Control", control), ("Hyperosmotic", hyper)): value = config.get("geometry", {}).get(THICKNESS) if value is None or abs(float(value) - thickness) > 1.0e-12: raise SystemExit( f"build_template_case_configs dropped the thickness override in " f"the {label} config ({value!r} instead of {thickness}). The " "sweep would be a no-op; fix the propagation before running." ) print(f" preflight: thickness override reaches both case configs ({thickness} um).") def _resolve_warm_start_root(args: argparse.Namespace, uq) -> Path: if args.warm_start_root is not None: return args.warm_start_root if args.primary_mesh == "coarse_h1p25": return Path(uq.NOMINAL_COARSE_P1_ROOT) return ( args.project_root / "results" / "reynolds_finan_probe_optimized_20260807" / "model_runs" / args.primary_mesh / "p1" ) def _load_warm_start(root: Path) -> tuple[np.ndarray, np.ndarray, float]: control = root / "control" / "nonlinear_nodal_displacement.npy" hyper = root / "hyperosmotic" / "nonlinear_nodal_displacement.npy" config = root / "hyperosmotic" / "case_config.json" missing = [str(path) for path in (control, hyper, config) if not path.exists()] if missing: raise SystemExit( "Missing nominal warm start on the target mesh:\n " + "\n ".join(missing) + "\nPass --warm-start-root pointing at the nominal run for this mesh. " "The array must be vertex-valued on the SAME mesh as the target run." ) volume_ratio = float( json.loads(config.read_text(encoding="utf-8"))["stimulus"]["volume_ratio"] ) return np.load(control), np.load(hyper), volume_ratio def _find_file(root: Path, names: tuple[str, ...]): for name in names: for candidate in (root / name, root / "05_Source_Data" / name): if candidate.exists(): return candidate matches = sorted(root.rglob(name)) if matches: return matches[0] return None def _write_csv(path: Path, rows: list[dict]) -> None: if not rows: return path.parent.mkdir(parents=True, exist_ok=True) fieldnames: list[str] = [] for row in rows: for key in row: if key not in fieldnames: fieldnames.append(key) with path.open("w", newline="", encoding="utf-8") as handle: writer = csv.DictWriter(handle, fieldnames=fieldnames) writer.writeheader() writer.writerows(rows) if __name__ == "__main__": raise SystemExit(main())