"""
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 "
"fintx/finty "
"registration in atom_vec_sph.cpp; "
"(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: atom_vec_sph.cpp "
"fintx/finty", h2),
Paragraph(
"The problem. "
"pair_sdpd_taitwater_isothermal_mf and "
"fix_sph both write to "
"atom->fintx / "
"atom->finty, and the flags "
"fintx_flag / finty_flag "
"exist in atom.cpp, but "
"atom_vec_sph.cpp never sets them or adds the names "
"to any of the fields_* lists. Without that "
"registration, the per-atom array is never grown — the pair-style force loop "
"writes fintx[i] = … 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(
"The fix. Set fintx_flag = 1 / "
"finty_flag = 1 and add "
"\"fintx\" / \"finty\" to "
"fields_grow, "
"fields_copy, "
"fields_comm, "
"fields_comm_vel, "
"fields_reverse, "
"fields_border, "
"fields_border_vel, "
"fields_exchange, "
"fields_restart, and "
"fields_create — the same registration pattern used for "
"gradv / gammadot.",
body),
Spacer(1, 0.25*cm),
Paragraph(
"Validation. Ran the bundled test_100 "
"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(
"Goal. Train a model that maps simulation inputs "
"(sigmao, dh) to a thermodynamic outcome of the "
"test_100 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(
"Setup. 12 quasi-random points sampled in "
"sigmao ∈ [0.5, 2.0], "
"dh ∈ [0.18, 0.25]. Each ran inside the "
"lampsui:latest Docker image with "
"prodRun=1000 (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(
"Lesson 1: target choice matters more than model architecture. "
"The first attempt — predicting LAMMPS' built-in box-averaged "
"Press 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 (sigmao, dh). "
"Switching the target to P_inside (mean per-atom pressure of the "
"droplet atoms only, computed from the dump's stress/atom "
"output) made the relationship cleanly linear in sigmao:",
body),
Spacer(1, 0.2*cm),
Image(str(target_png), width=15.5*cm, height=6.5*cm),
PageBreak(),
Paragraph(
"Lesson 2: at n=12, simpler models win. Comparing four off-the-shelf "
"regressors via leave-one-out cross-validation on "
"(sigmao, dh) → P_in:",
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(
"Speedup. The trained GP predicts 1 000 new "
"(sigmao, dh) 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(
"Honest caveats. "
"(a) The surrogate is only valid in the box "
"sigmao ∈ [0.5, 2.0], "
"dh ∈ [0.18, 0.25]; outside this range "
"predictions are extrapolation. "
"(b) It learns one outcome (P_inside) for one example "
"(test_100); other examples or other quantities "
"would need their own training data. "
"(c) P_inside 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(
"Possible next steps. "
"(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 "
"(f2Frac, eta11, …) "
"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(
"Files produced. "
"scripts/surrogate_demo.py, "
"scripts/surrogate_postprocess.py, "
"scripts/build_report.py, "
"work/surrogate_demo/data.csv, "
"work/surrogate_demo/data_with_dp.csv, "
"work/surrogate_demo/figs/*.png, "
"work/surrogate_demo/report.pdf (this file). "
"Everything under work/ 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)")