OneScience-Group/MassConservingCNN

1"""Compute paper-aligned metrics and plot input, target, and prediction."""2 3import json4from pathlib import Path5 6import matplotlib7matplotlib.use("Agg")8import matplotlib.pyplot as plt9import numpy as np10import yaml11 12 13ROOT = Path(__file__).resolve().parents[1]14 15 16def metrics(candidate_n, candidate_p, target_n, target_p):17    rmse = np.sqrt(np.mean((candidate_n - target_n) ** 2, axis=(0, 2)))18    sample_variable_rmse = np.sqrt(np.mean((candidate_n - target_n) ** 2, axis=2))19    mass_h = np.mean(np.abs(candidate_p[:, 1].sum(1) - target_p[:, 1].sum(1)) / candidate_p.shape[2])20    mass_r = np.mean(np.abs(candidate_p[:, 2].sum(1) - target_p[:, 2].sum(1)) / candidate_p.shape[2])21    h_bias = np.mean(candidate_p[:, 1] - target_p[:, 1])22    return {"J": float(sample_variable_rmse.mean()), "rmse_u": float(rmse[0]),23            "rmse_h": float(rmse[1]), "rmse_r": float(rmse[2]),24            "mass_error_h_per_point": float(mass_h), "mass_error_r_per_point": float(mass_r),25            "h_bias": float(h_bias)}26 27 28def main():29    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())30    data = np.load(ROOT / config["paths"]["inference"])31    if str(data["format_version"]) != config["data"]["format_version"]:32        raise ValueError("prediction format mismatch")33    input_n, target_n, prediction_n = data["inputs"][:, :3], data["targets"], data["predictions"]34    input_p, target_p, prediction_p = data["xa"], data["targets_physical"], data["predictions_physical"]35    expected = (len(target_n), 3, 250)36    if any(array.shape != expected for array in (input_n, target_n, prediction_n, input_p, target_p, prediction_p)):37        raise ValueError("evaluation arrays must have shape [B,3,250]")38    baseline = metrics(input_n, input_p, target_n, target_p)39    prediction = metrics(prediction_n, prediction_p, target_n, target_p)40    improvement = {key: float(100 * (baseline[key] - prediction[key]) / baseline[key])41                   for key in ("J", "rmse_u", "rmse_h", "rmse_r", "mass_error_h_per_point", "mass_error_r_per_point")42                   if baseline[key] != 0}43    if baseline["h_bias"] != 0:44        improvement["absolute_h_bias"] = float(45            100 * (abs(baseline["h_bias"]) - abs(prediction["h_bias"])) / abs(baseline["h_bias"])46        )47    report = {"samples": len(target_n), "baseline_input": baseline, "prediction": prediction,48              "relative_improvement_percent": improvement,49              "note": "Structured synthetic engineering validation; not paper performance."}50    values = list(baseline.values()) + list(prediction.values()) + list(improvement.values())51    if not np.isfinite(values).all():52        raise FloatingPointError("non-finite evaluation metric")53    output = ROOT / config["paths"]["evaluation_dir"]54    output.mkdir(parents=True, exist_ok=True)55    (output / "metrics.json").write_text(json.dumps(report, indent=2) + "\n")56    x = np.arange(250) * float(config["data"]["domain_km"]) / 25057    figure, axes = plt.subplots(3, 1, figsize=(11, 8), sharex=True)58    for index, (axis, variable) in enumerate(zip(axes, ("u", "h", "r"))):59        axis.plot(x, input_p[0, index], color="steelblue", label="input X^a", linewidth=1.4)60        axis.plot(x, target_p[0, index], color="black", label="QPEns target", linewidth=1.5)61        axis.plot(x, prediction_p[0, index], color="firebrick", label="CNN prediction", linewidth=1.3)62        if variable == "r":63            axis.fill_between(x, 0, data["radar"][0, 0] * max(target_p[0, 2].max(), 1e-6), color="gold", alpha=0.2, label="radar mask")64        axis.set_ylabel(variable); axis.grid(alpha=0.2)65    axes[0].legend(ncol=3); axes[-1].set_xlabel("distance (km)")66    figure.suptitle("MassConservingCNN structured synthetic validation")67    figure.tight_layout(); figure.savefig(output / "input_target_prediction.png", dpi=150); plt.close(figure)68    print(f"evaluation={output.relative_to(ROOT)} J={prediction['J']:.6f} h_mass={prediction['mass_error_h_per_point']:.6f}")69 70 71if __name__ == "__main__":72    main()73