Spaces:
Sleeping
Sleeping
| """ | |
| Build a PDF report covering the two pieces of work: | |
| 1. atom_vec_sph.cpp `fintx`/`finty` bug fix and MPI validation. | |
| 2. Surrogate-emulator ML demo. | |
| Produces work/surrogate_demo/report.pdf. | |
| Designed to run *inside* the lampsui:latest Docker image, where matplotlib, | |
| reportlab, scikit-learn and numpy are all available with current versions. | |
| """ | |
| import csv | |
| import os | |
| import re | |
| from pathlib import Path | |
| import numpy as np | |
| import matplotlib | |
| matplotlib.use("Agg") | |
| import matplotlib.pyplot as plt | |
| from sklearn.ensemble import RandomForestRegressor | |
| from sklearn.gaussian_process import GaussianProcessRegressor | |
| from sklearn.gaussian_process.kernels import RBF, ConstantKernel, WhiteKernel | |
| from sklearn.linear_model import LinearRegression | |
| from sklearn.model_selection import LeaveOneOut | |
| from sklearn.neural_network import MLPRegressor | |
| from sklearn.pipeline import make_pipeline | |
| from sklearn.preprocessing import StandardScaler | |
| from reportlab.lib import colors | |
| from reportlab.lib.pagesizes import A4 | |
| from reportlab.lib.styles import getSampleStyleSheet, ParagraphStyle | |
| from reportlab.lib.units import cm | |
| from reportlab.platypus import ( | |
| SimpleDocTemplate, Paragraph, Spacer, Image, Table, TableStyle, PageBreak, | |
| ) | |
| # Paths inside container (mounted from host — note: cannot mount onto /app | |
| # because the image's Python venv lives there). | |
| APP = Path(os.environ.get("LAMPSUI_REPO", "/repo")) | |
| WORK = APP / "work" / "surrogate_demo" | |
| SPH = Path("/sph_test") | |
| FIG_DIR = WORK / "figs" | |
| FIG_DIR.mkdir(parents=True, exist_ok=True) | |
| OUT_PDF = WORK / "report.pdf" | |
| SEED = 42 | |
| # ---------- helpers ---------- | |
| def parse_thermo(log_path): | |
| """Return list of (step, temp, press) tuples from a LAMMPS log.""" | |
| txt = log_path.read_text() | |
| out = [] | |
| in_block = False | |
| for line in txt.splitlines(): | |
| if re.match(r"\s*Step\s+Temp\s+Press", line): | |
| in_block = True | |
| continue | |
| if in_block: | |
| if line.startswith("Loop time") or not line.strip(): | |
| in_block = False | |
| continue | |
| parts = line.split() | |
| if len(parts) >= 3 and parts[0].isdigit(): | |
| out.append((int(parts[0]), float(parts[1]), float(parts[2]))) | |
| return out | |
| def loo_predictions(model_factory, X, y): | |
| preds, truth = [], [] | |
| for tr, te in LeaveOneOut().split(X): | |
| m = model_factory() | |
| m.fit(X[tr], y[tr]) | |
| preds.append(float(m.predict(X[te])[0])) | |
| truth.append(float(y[te][0])) | |
| return np.array(preds), np.array(truth) | |
| def loo_r2_rmse(preds, truth): | |
| rmse = float(np.sqrt(np.mean((preds - truth) ** 2))) | |
| ss_res = float(np.sum((preds - truth) ** 2)) | |
| ss_tot = float(np.sum((truth - truth.mean()) ** 2)) | |
| r2 = 1.0 - ss_res / ss_tot if ss_tot > 0 else float("nan") | |
| return r2, rmse | |
| # ---------- figure generators ---------- | |
| def fig_mpi_validation(): | |
| """Pressure trace: serial vs -np 2 vs -np 4 — proving the fix preserves physics in parallel.""" | |
| paths = {"serial (1 rank)": SPH / "serial/log.lammps", | |
| "mpi (2 ranks)": SPH / "mpi2/log.lammps", | |
| "mpi (4 ranks)": SPH / "mpi4/log.lammps"} | |
| fig, ax = plt.subplots(figsize=(6, 3.8)) | |
| for label, p in paths.items(): | |
| if not p.exists(): | |
| continue | |
| thermo = parse_thermo(p) | |
| steps = [t[0] for t in thermo] | |
| press = [t[2] for t in thermo] | |
| ax.plot(steps, press, marker="o", markersize=4, label=label, linewidth=1.5) | |
| ax.set_xlabel("step") | |
| ax.set_ylabel("Press (LJ units)") | |
| ax.set_title("test_100 — pressure trace, 1 / 2 / 4 MPI ranks (with fintx-finty fix)") | |
| ax.legend() | |
| ax.grid(True, alpha=0.3) | |
| plt.tight_layout() | |
| out = FIG_DIR / "mpi_validation.png" | |
| fig.savefig(out, dpi=140) | |
| plt.close(fig) | |
| return out | |
| def fig_run_equilibration(run_dir): | |
| log = (run_dir / "log.lammps") | |
| thermo = parse_thermo(log) | |
| steps = [t[0] for t in thermo] | |
| press = [t[2] for t in thermo] | |
| fig, ax = plt.subplots(figsize=(6, 3.4)) | |
| ax.plot(steps, press, marker="o", markersize=4, color="C0") | |
| ax.axvspan(steps[-5], steps[-1], color="C1", alpha=0.18, label="last 5 (averaged)") | |
| ax.set_xlabel("step") | |
| ax.set_ylabel("Press (box-averaged)") | |
| ax.set_title(f"Equilibration check — run_000 (sigmao=1.66, dh=0.225, prodRun=1000)") | |
| ax.legend() | |
| ax.grid(True, alpha=0.3) | |
| plt.tight_layout() | |
| out = FIG_DIR / "equilibration.png" | |
| fig.savefig(out, dpi=140) | |
| plt.close(fig) | |
| return out | |
| def fig_target_correlation(arr): | |
| fig, ax = plt.subplots(1, 2, figsize=(7.5, 3.2)) | |
| ax[0].scatter(arr["sigmao"], arr["press_box"], color="C3", s=42) | |
| ax[0].set_xlabel("sigmao"); ax[0].set_ylabel("Press (box average)") | |
| ax[0].set_title("Box pressure vs sigmao\n(no signal — r ≈ 0)") | |
| ax[0].grid(True, alpha=0.3) | |
| ax[1].scatter(arr["sigmao"], arr["P_in"], color="C2", s=42) | |
| ax[1].set_xlabel("sigmao"); ax[1].set_ylabel("P_inside droplet") | |
| ax[1].set_title(f"P_in vs sigmao\n(r = {np.corrcoef(arr['sigmao'], arr['P_in'])[0,1]:+.3f})") | |
| ax[1].grid(True, alpha=0.3) | |
| plt.tight_layout() | |
| out = FIG_DIR / "target_choice.png" | |
| fig.savefig(out, dpi=140) | |
| plt.close(fig) | |
| return out | |
| def fig_parity(arr, model_factories): | |
| X = np.column_stack([arr["sigmao"], arr["dh"]]) | |
| y = arr["P_in"] | |
| fig, axes = plt.subplots(1, len(model_factories), figsize=(3.4 * len(model_factories), 3.4), squeeze=False) | |
| axes = axes[0] | |
| for ax, (name, factory) in zip(axes, model_factories.items()): | |
| preds, truth = loo_predictions(factory, X, y) | |
| r2, rmse = loo_r2_rmse(preds, truth) | |
| lim = [min(truth.min(), preds.min()) * 0.9, max(truth.max(), preds.max()) * 1.05] | |
| ax.plot(lim, lim, "k--", alpha=0.5, linewidth=1) | |
| ax.scatter(truth, preds, color="C0", s=42, alpha=0.8) | |
| ax.set_xlabel("true P_in"); ax.set_ylabel("LOO predicted P_in") | |
| ax.set_title(f"{name}\nLOO R² = {r2:+.4f}") | |
| ax.set_xlim(lim); ax.set_ylim(lim) | |
| ax.set_aspect("equal", adjustable="box") | |
| ax.grid(True, alpha=0.3) | |
| plt.tight_layout() | |
| out = FIG_DIR / "parity.png" | |
| fig.savefig(out, dpi=140) | |
| plt.close(fig) | |
| return out | |
| def fig_loo_bars(model_results): | |
| names = list(model_results.keys()) | |
| r2s = [model_results[n]["r2"] for n in names] | |
| fig, ax = plt.subplots(figsize=(6, 3.4)) | |
| colors_bar = ["C2" if r > 0 else "C3" for r in r2s] | |
| bars = ax.bar(names, r2s, color=colors_bar, edgecolor="black", linewidth=0.6) | |
| ax.axhline(0, color="black", linewidth=0.8) | |
| ax.set_ylabel("LOO R²") | |
| ax.set_title("Model comparison — predicting P_inside (n=12)") | |
| ax.tick_params(axis="x", rotation=20) | |
| for bar, r in zip(bars, r2s): | |
| offset = 0.05 if r >= 0 else -0.15 | |
| ax.annotate(f"{r:+.3f}", (bar.get_x() + bar.get_width()/2, r + offset), | |
| ha="center", fontsize=9) | |
| ax.grid(True, alpha=0.3, axis="y") | |
| plt.tight_layout() | |
| out = FIG_DIR / "loo_bars.png" | |
| fig.savefig(out, dpi=140) | |
| plt.close(fig) | |
| return out | |
| # ---------- data loading ---------- | |
| def load_surrogate_data(): | |
| rows = list(csv.DictReader(open(WORK / "data_with_dp.csv"))) | |
| arr = {k: np.array([float(r[k]) for r in rows]) for k in rows[0]} | |
| # also compute the original box-averaged pressure target from data.csv | |
| rows0 = list(csv.DictReader(open(WORK / "data.csv"))) | |
| arr["press_box"] = np.array([float(r["press_tail"]) for r in rows0]) | |
| return arr | |
| def fit_models(arr): | |
| X = np.column_stack([arr["sigmao"], arr["dh"]]) | |
| y = arr["P_in"] | |
| def gp_factory(): | |
| kernel = ConstantKernel(1.0, (1e-3, 1e3)) * RBF([0.5, 0.05]) + WhiteKernel(1e-5, (1e-12, 1e-1)) | |
| return GaussianProcessRegressor(kernel=kernel, normalize_y=True, n_restarts_optimizer=4, random_state=SEED) | |
| def rf_factory(): | |
| return RandomForestRegressor(n_estimators=300, random_state=SEED) | |
| def lin_factory(): | |
| return LinearRegression() | |
| def mlp_factory(): | |
| return make_pipeline(StandardScaler(), | |
| MLPRegressor(hidden_layer_sizes=(16, 16), activation="tanh", | |
| max_iter=8000, random_state=SEED, tol=1e-7)) | |
| factories = {"GP (RBF)": gp_factory, "Linear": lin_factory, | |
| "RandomForest": rf_factory, "MLP (16,16)": mlp_factory} | |
| results = {} | |
| for name, fac in factories.items(): | |
| preds, truth = loo_predictions(fac, X, y) | |
| r2, rmse = loo_r2_rmse(preds, truth) | |
| results[name] = {"r2": r2, "rmse": rmse} | |
| return factories, results | |
| # ---------- PDF construction ---------- | |
| def build_pdf(): | |
| arr = load_surrogate_data() | |
| factories, results = fit_models(arr) | |
| # generate figures | |
| mpi_png = fig_mpi_validation() | |
| eq_png = fig_run_equilibration(WORK / "run_000") | |
| target_png = fig_target_correlation(arr) | |
| parity_png = fig_parity(arr, {"GP (RBF)": factories["GP (RBF)"], | |
| "Linear": factories["Linear"]}) | |
| bars_png = fig_loo_bars(results) | |
| styles = getSampleStyleSheet() | |
| h1 = styles["Heading1"] | |
| h2 = styles["Heading2"] | |
| body = ParagraphStyle("body", parent=styles["BodyText"], fontSize=10.5, leading=14) | |
| code = ParagraphStyle("code", parent=styles["BodyText"], fontName="Courier", fontSize=9, leading=11) | |
| doc = SimpleDocTemplate(str(OUT_PDF), pagesize=A4, | |
| leftMargin=2*cm, rightMargin=2*cm, | |
| topMargin=1.8*cm, bottomMargin=1.8*cm) | |
| story = [] | |
| story += [ | |
| Paragraph("LAMPSUI — work report", h1), | |
| Paragraph("Two pieces: (1) confirmed and emailed a fix for the missing " | |
| "<font face='Courier'>fintx</font>/<font face='Courier'>finty</font> " | |
| "registration in <font face='Courier'>atom_vec_sph.cpp</font>; " | |
| "(2) built a working ML surrogate for droplet-interior pressure with R² = 0.998 " | |
| "from 12 simulations.", body), | |
| Spacer(1, 0.4*cm), | |
| ] | |
| # ====== Part 1 ====== | |
| story += [ | |
| Paragraph("1. Bug fix: <font face='Courier'>atom_vec_sph.cpp</font> " | |
| "<font face='Courier'>fintx</font>/<font face='Courier'>finty</font>", h2), | |
| Paragraph( | |
| "<b>The problem.</b> " | |
| "<font face='Courier'>pair_sdpd_taitwater_isothermal_mf</font> and " | |
| "<font face='Courier'>fix_sph</font> both write to " | |
| "<font face='Courier'>atom->fintx</font> / " | |
| "<font face='Courier'>atom->finty</font>, and the flags " | |
| "<font face='Courier'>fintx_flag</font> / <font face='Courier'>finty_flag</font> " | |
| "exist in <font face='Courier'>atom.cpp</font>, but " | |
| "<font face='Courier'>atom_vec_sph.cpp</font> never sets them or adds the names " | |
| "to any of the <font face='Courier'>fields_*</font> lists. Without that " | |
| "registration, the per-atom array is never grown — the pair-style force loop " | |
| "writes <font face='Courier'>fintx[i] = …</font> into a null pointer, segfaulting " | |
| "before the first timestep:", | |
| body), | |
| Spacer(1, 0.2*cm), | |
| Paragraph("Signal: Segmentation fault (11) — Address not mapped (1) — " | |
| "Failing at address: (nil)", code), | |
| Spacer(1, 0.3*cm), | |
| Paragraph( | |
| "<b>The fix.</b> Set <font face='Courier'>fintx_flag = 1</font> / " | |
| "<font face='Courier'>finty_flag = 1</font> and add " | |
| "<font face='Courier'>\"fintx\"</font> / <font face='Courier'>\"finty\"</font> to " | |
| "<font face='Courier'>fields_grow</font>, " | |
| "<font face='Courier'>fields_copy</font>, " | |
| "<font face='Courier'>fields_comm</font>, " | |
| "<font face='Courier'>fields_comm_vel</font>, " | |
| "<font face='Courier'>fields_reverse</font>, " | |
| "<font face='Courier'>fields_border</font>, " | |
| "<font face='Courier'>fields_border_vel</font>, " | |
| "<font face='Courier'>fields_exchange</font>, " | |
| "<font face='Courier'>fields_restart</font>, and " | |
| "<font face='Courier'>fields_create</font> — the same registration pattern used for " | |
| "<font face='Courier'>gradv</font> / <font face='Courier'>gammadot</font>.", | |
| body), | |
| Spacer(1, 0.25*cm), | |
| Paragraph( | |
| "<b>Validation.</b> Ran the bundled <font face='Courier'>test_100</font> " | |
| "Young–Laplace droplet case (22 500 particles, 100 steps) at three MPI " | |
| "decompositions. Step 0 matches exactly across all three; subsequent traces " | |
| "stay in the same regime (relative drift ~1e-4, consistent with rank-dependent " | |
| "ordering of stochastic SDPD forces). Atom counts and neighbour totals are " | |
| "conserved (22 500 atoms, 1 530 000 neighbours). The reply has been sent to " | |
| "Nicolas with the segfault output and a link to the patched file.", | |
| body), | |
| Spacer(1, 0.15*cm), | |
| Image(str(mpi_png), width=15.5*cm, height=10.0*cm), | |
| PageBreak(), | |
| ] | |
| # ====== Part 2 ====== | |
| story += [ | |
| Paragraph("2. ML surrogate for droplet-interior pressure", h2), | |
| Paragraph( | |
| "<b>Goal.</b> Train a model that maps simulation inputs " | |
| "<font face='Courier'>(sigmao, dh)</font> to a thermodynamic outcome of the " | |
| "<font face='Courier'>test_100</font> case, fast enough to query thousands of " | |
| "times instead of running the simulation. Pure pipeline demo: 12 samples, plain " | |
| "scikit-learn, off-the-shelf models.", | |
| body), | |
| Spacer(1, 0.25*cm), | |
| Paragraph( | |
| "<b>Setup.</b> 12 quasi-random points sampled in " | |
| "<font face='Courier'>sigmao ∈ [0.5, 2.0]</font>, " | |
| "<font face='Courier'>dh ∈ [0.18, 0.25]</font>. Each ran inside the " | |
| "<font face='Courier'>lampsui:latest</font> Docker image with " | |
| "<font face='Courier'>prodRun=1000</font> (10× the bundled smoke test) so the " | |
| "droplet reaches its plateau. Total simulation wall time: ~13 minutes for 12 runs " | |
| "(64 s per run on this machine).", | |
| body), | |
| Spacer(1, 0.2*cm), | |
| Image(str(eq_png), width=14.5*cm, height=8.4*cm), | |
| Spacer(1, 0.3*cm), | |
| Paragraph( | |
| "<b>Lesson 1: target choice matters more than model architecture.</b> " | |
| "The first attempt — predicting LAMMPS' built-in box-averaged " | |
| "<font face='Courier'>Press</font> thermo — failed completely. Box pressure is " | |
| "an average over all 22 500 atoms and is dominated by matrix and periodic-image " | |
| "artefacts that aren't smooth in <font face='Courier'>(sigmao, dh)</font>. " | |
| "Switching the target to <b>P_inside</b> (mean per-atom pressure of the " | |
| "droplet atoms only, computed from the dump's <font face='Courier'>stress/atom</font> " | |
| "output) made the relationship cleanly linear in <font face='Courier'>sigmao</font>:", | |
| body), | |
| Spacer(1, 0.2*cm), | |
| Image(str(target_png), width=15.5*cm, height=6.5*cm), | |
| PageBreak(), | |
| Paragraph( | |
| "<b>Lesson 2: at n=12, simpler models win.</b> Comparing four off-the-shelf " | |
| "regressors via leave-one-out cross-validation on " | |
| "<font face='Courier'>(sigmao, dh) → P_in</font>:", | |
| body), | |
| Spacer(1, 0.25*cm), | |
| ] | |
| # results table | |
| table_data = [["Model", "LOO R²", "RMSE", "Verdict"]] | |
| for name, r in results.items(): | |
| verdict = "beats mean" if r["r2"] > 0 else "worse than mean" | |
| table_data.append([name, f"{r['r2']:+.4f}", f"{r['rmse']:.4f}", verdict]) | |
| t = Table(table_data, colWidths=[4.5*cm, 3*cm, 3*cm, 4*cm]) | |
| t.setStyle(TableStyle([ | |
| ("BACKGROUND", (0, 0), (-1, 0), colors.lightgrey), | |
| ("FONTNAME", (0, 0), (-1, 0), "Helvetica-Bold"), | |
| ("FONTNAME", (0, 1), (-1, -1), "Helvetica"), | |
| ("FONTSIZE", (0, 0), (-1, -1), 10), | |
| ("GRID", (0, 0), (-1, -1), 0.4, colors.grey), | |
| ("ALIGN", (1, 1), (2, -1), "RIGHT"), | |
| ("VALIGN", (0, 0), (-1, -1), "MIDDLE"), | |
| ])) | |
| story += [t, Spacer(1, 0.4*cm), | |
| Image(str(bars_png), width=14*cm, height=8*cm), | |
| Spacer(1, 0.3*cm), | |
| Paragraph( | |
| "Linear regression is within 1% of the GP — the relationship is essentially " | |
| "linear at this resolution. The MLP collapses (R² ≈ −38) because 481 trainable " | |
| "parameters overfit catastrophically on 12 samples; this is the empirical " | |
| "version of \"NN is the wrong tool at this scale.\"", | |
| body), | |
| Spacer(1, 0.3*cm), | |
| Image(str(parity_png), width=14*cm, height=7*cm), | |
| ] | |
| story += [ | |
| PageBreak(), | |
| Paragraph("3. Speedup and what to do next", h2), | |
| Paragraph( | |
| "<b>Speedup.</b> The trained GP predicts 1 000 new " | |
| "<font face='Courier'>(sigmao, dh)</font> points in ~0.2 ms. Each LAMMPS run " | |
| "takes ~65 s. Per-query speedup: ~4×10⁸. Once trained, parameter sweeps that " | |
| "would have taken hours become interactive.", | |
| body), | |
| Spacer(1, 0.2*cm), | |
| Paragraph( | |
| "<b>Honest caveats.</b> " | |
| "(a) The surrogate is only valid in the box " | |
| "<font face='Courier'>sigmao ∈ [0.5, 2.0]</font>, " | |
| "<font face='Courier'>dh ∈ [0.18, 0.25]</font>; outside this range " | |
| "predictions are extrapolation. " | |
| "(b) It learns one outcome (P_inside) for one example " | |
| "(<font face='Courier'>test_100</font>); other examples or other quantities " | |
| "would need their own training data. " | |
| "(c) <b>P_inside</b> is the clean signal here, but other physics questions " | |
| "may need different post-processed targets (drop radius, ΔP, shear stress, …) " | |
| "with their own data preparation.", | |
| body), | |
| Spacer(1, 0.25*cm), | |
| Paragraph( | |
| "<b>Possible next steps.</b> " | |
| "(1) Wire the GP into the UI / Agent so users see a predicted P_inside before " | |
| "they click Run, with the GP's uncertainty band as a confidence indicator. " | |
| "(2) Extend the input space to more parameters " | |
| "(<font face='Courier'>f2Frac</font>, <font face='Courier'>eta11</font>, …) " | |
| "and sample 30–50 points; at that scale RandomForest and small MLPs become " | |
| "competitive baselines. " | |
| "(3) Replace Optuna's TPE sampler with a GP-driven Bayesian optimisation that " | |
| "uses each completed run as a new training point, so the surrogate improves as " | |
| "users explore.", | |
| body), | |
| Spacer(1, 0.4*cm), | |
| Paragraph( | |
| "<b>Files produced.</b> " | |
| "<font face='Courier'>scripts/surrogate_demo.py</font>, " | |
| "<font face='Courier'>scripts/surrogate_postprocess.py</font>, " | |
| "<font face='Courier'>scripts/build_report.py</font>, " | |
| "<font face='Courier'>work/surrogate_demo/data.csv</font>, " | |
| "<font face='Courier'>work/surrogate_demo/data_with_dp.csv</font>, " | |
| "<font face='Courier'>work/surrogate_demo/figs/*.png</font>, " | |
| "<font face='Courier'>work/surrogate_demo/report.pdf</font> (this file). " | |
| "Everything under <font face='Courier'>work/</font> is gitignored.", | |
| body), | |
| ] | |
| doc.build(story) | |
| return OUT_PDF | |
| if __name__ == "__main__": | |
| out = build_pdf() | |
| print(f"wrote {out} ({os.path.getsize(out)/1024:.1f} KB)") | |