#!/usr/bin/env python3 """Create the final FAST-DeepONet Darcy figures and reproduction report.""" from __future__ import annotations import argparse import json import os import sys from pathlib import Path import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt import numpy as np import yaml PROJECT_ROOT = Path(__file__).resolve().parents[1] if str(PROJECT_ROOT) not in sys.path: sys.path.insert(0, str(PROJECT_ROOT)) def load_config(path: Path) -> dict[str, object]: """Load visualization configuration without importing the Torch runtime.""" with path.resolve().open("r", encoding="utf-8") as handle: config = yaml.safe_load(handle) if not isinstance(config, dict): raise ValueError(f"Config must be a mapping: {path}") return config def atomic_save_figure(figure: plt.Figure, path: Path, dpi: int = 180) -> None: temporary = path.with_name(path.stem + ".tmp" + path.suffix) figure.savefig(temporary, dpi=dpi, bbox_inches="tight") os.replace(temporary, path) def plot_training(history: list[dict[str, float]], results_dir: Path) -> None: epochs = np.asarray([row["epoch"] for row in history]) train_mse = np.asarray([row["train_mse"] for row in history]) tune_mse = np.asarray([row["tune_mse"] for row in history]) train_total = np.asarray([row["train_total_loss"] for row in history]) figure, axis = plt.subplots(figsize=(7.5, 4.8)) axis.semilogy(epochs, train_total, label="train objective", linewidth=1.4) axis.semilogy(epochs, train_mse, label="train MSE", linewidth=1.4) axis.semilogy(epochs, tune_mse, label="tune MSE", linewidth=1.6) axis.set_xlabel("Epoch") axis.set_ylabel("Loss") axis.set_title("FAST-DeepONet Darcy training") axis.grid(True, which="both", alpha=0.25) axis.legend() atomic_save_figure(figure, results_dir / "training_curves.png") plt.close(figure) def plot_sample_metrics(metrics: dict[str, object], results_dir: Path) -> None: values = np.asarray([sample["relative_l2"] for sample in metrics["samples"]], dtype=float) figure, axis = plt.subplots(figsize=(7.5, 4.5)) bars = axis.bar(np.arange(values.size), values, color="#4C78A8") axis.axhline(float(metrics["mean_relative_l2"]), color="#E45756", linestyle="--", label="mean") axis.set_xticks(np.arange(values.size), [str(index) for index in range(values.size)]) axis.set_xlabel("Test sample index") axis.set_ylabel("Relative L2") axis.set_title("Per-sample test error") axis.grid(True, axis="y", alpha=0.25) axis.legend() for bar, value in zip(bars, values): axis.text(bar.get_x() + bar.get_width() / 2, value, f"{value:.3f}", ha="center", va="bottom", fontsize=8) atomic_save_figure(figure, results_dir / "test_relative_l2.png") plt.close(figure) def plot_field_comparison(prediction_path: Path, grid_size: int, results_dir: Path) -> int: with np.load(prediction_path, allow_pickle=False) as data: errors = np.asarray(data["relative_l2"]) sample_index = int(np.argsort(errors)[errors.size // 2]) log_permeability = np.asarray(data["inputs"][sample_index]).reshape(grid_size, grid_size) target = np.asarray(data["targets"][sample_index]).reshape(grid_size, grid_size) prediction = np.asarray(data["predictions"][sample_index]).reshape(grid_size, grid_size) mask = np.asarray(data["masks"][sample_index]).reshape(grid_size, grid_size) filename = str(data["filenames"][sample_index]) absolute_error = np.abs(target - prediction) figure, axes = plt.subplots(1, 4, figsize=(16, 3.8), constrained_layout=True) panels = [ (log_permeability, "log permeability", "coolwarm"), (target, "pressure target", "coolwarm"), (prediction, "pressure prediction", "coolwarm"), (absolute_error, "absolute error", "magma"), ] pressure_min = min(float(target.min()), float(prediction.min())) pressure_max = max(float(target.max()), float(prediction.max())) for panel_index, (axis, (field, title, cmap)) in enumerate(zip(axes, panels)): kwargs = {"origin": "lower", "extent": (0, 1, 0, 1), "cmap": cmap} if panel_index in (1, 2): kwargs.update(vmin=pressure_min, vmax=pressure_max) image = axis.imshow(field.T, **kwargs) if np.any(mask > 0.5) and np.any(mask <= 0.5): axis.contour( np.linspace(0, 1, grid_size), np.linspace(0, 1, grid_size), mask.T, levels=[0.5], colors="white" if panel_index != 0 else "black", linewidths=0.5, alpha=0.7, ) axis.set_title(title) axis.set_xlabel("x") if panel_index == 0: axis.set_ylabel("y") figure.colorbar(image, ax=axis, shrink=0.82) figure.suptitle(f"Test sample {sample_index}: {filename} (relative L2={errors[sample_index]:.5f})") atomic_save_figure(figure, results_dir / "darcy_field_comparison.png") plt.close(figure) return sample_index def write_report( path: Path, training_summary: dict[str, object], metrics: dict[str, object], representative_index: int, ) -> None: lines = [ "# FAST-DeepONet Darcy reproduction report", "", "## Formal training", "", f"- Status: `{training_summary['status']}`", f"- Completed epochs: {training_summary['completed_epochs']} / {training_summary['requested_epochs']}", f"- Best epoch: {training_summary['best_epoch']}", f"- Best tune MSE: {float(training_summary['best_tune_mse']):.10e}", f"- Effective factor rank: {training_summary['effective_factor_rank']}", f"- Device: `{training_summary['device_name']}`", f"- SLURM job: `{training_summary['slurm_job_id']}`", f"- Checkpoint: `{training_summary['checkpoint']}`", "", "## Test relative L2", "", ] for sample in metrics["samples"]: lines.append( f"- {int(sample['sample_index']):02d} `{sample['filename']}`: {float(sample['relative_l2']):.10e}" ) lines.extend( [ "", f"Mean relative L2: **{float(metrics['mean_relative_l2']):.10e}**", "", "## Artifacts and interpretation", "", f"- The four-panel field figure uses representative sample index {representative_index}.", "- `training_curves.png` shows the actual 300-epoch objective/train MSE/tune MSE history.", "- `test_relative_l2.png` shows all ten fixed test fields.", "- Calibration data were not used because this task requests point prediction, not UQ calibration.", "- This run uses 10 test fields and one model seed. Table 2 reports 100 independent test fields and summarizes five fits, so its 0.0925±0.0019 value is not directly comparable to this fixed-data run.", "", ] ) temporary = path.with_suffix(path.suffix + ".tmp") temporary.write_text("\n".join(lines), encoding="utf-8") os.replace(temporary, path) def run(config: dict[str, object]) -> None: results_dir = Path(config["paths"]["results_dir"]) history_path = results_dir / "training_history.json" summary_path = results_dir / "training_summary.json" metrics_path = results_dir / "test_metrics.json" predictions_path = results_dir / "test_predictions.npz" for path in (history_path, summary_path, metrics_path, predictions_path): if not path.is_file(): raise FileNotFoundError(f"Required result is missing: {path}") history = json.loads(history_path.read_text(encoding="utf-8")) training_summary = json.loads(summary_path.read_text(encoding="utf-8")) metrics = json.loads(metrics_path.read_text(encoding="utf-8")) if len(history) != 300 or int(training_summary["completed_epochs"]) != 300: raise ValueError("Visualization requires the completed 300-epoch formal run") if len(metrics["samples"]) != 10: raise ValueError("Visualization requires all ten test metrics") plot_training(history, results_dir) plot_sample_metrics(metrics, results_dir) representative_index = plot_field_comparison( predictions_path, int(config["data"]["grid_size"]), results_dir ) write_report(results_dir / "reproduction_report.md", training_summary, metrics, representative_index) print( f"visualization_complete mean_relative_l2={float(metrics['mean_relative_l2']):.10e} " f"results_dir={results_dir}", flush=True, ) def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--config", type=Path, default=PROJECT_ROOT / "config" / "config.yaml") return parser.parse_args() def main() -> None: args = parse_args() run(load_config(args.config)) if __name__ == "__main__": main()