scripts/result.py

"""Summarize DKE growth, compensation, correlations, maps and spectra."""

import json
from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import yaml


ROOT = Path(__file__).resolve().parents[1]


def weighted_map_correlation_matrix(maps, latitudes):
    weights = np.broadcast_to(np.cos(np.deg2rad(latitudes))[:, None], maps.shape[1:]).reshape(-1)
    weights = weights / weights.sum()
    flattened = maps.reshape(maps.shape[0], -1).astype(np.float64)
    centered = flattened - np.sum(flattened * weights, axis=1, keepdims=True)
    covariance = (centered * weights) @ centered.T
    scale = np.sqrt(np.maximum(np.diag(covariance), np.finfo(float).tiny))
    return covariance / np.outer(scale, scale)


def main():
    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())
    data = np.load(ROOT / config["paths"]["inference"])
    dke, maps, spectra = data["global_dke"], data["dke_maps"], data["spectra"]
    if dke.shape != (5, 73) or maps.shape != (5, 73, 721, 1440) or spectra.shape != (5, 73, 720):
        raise ValueError("inference product dimensions are incomplete")
    growth = dke / np.maximum(dke[:, :1], np.finfo(float).tiny)
    compensation = dke / np.maximum(dke[:1], np.finfo(float).tiny)
    per_step = dke[:, 1:] / np.maximum(dke[:, :-1], np.finfo(float).tiny)
    map_index = int(config["evaluation"]["map_hour"] // config["data"]["step_hours"])
    corr = weighted_map_correlation_matrix(maps[:, map_index], data["latitudes_degrees"])
    report = {
        "training_required": False,
        "experiments": data["experiments"].tolist(),
        "logical_field_shape": data["logical_field_shape"].tolist(),
        "logical_spectral_shape": data["logical_spectral_shape"].tolist(),
        "global_dke_0h": dke[:, 0].tolist(),
        "global_dke_72h": dke[:, -1].tolist(),
        "growth_factor_72h": growth[:, -1].tolist(),
        "scaling_compensation_ratio_72h": compensation[:, -1].tolist(),
        "coslat_weighted_72h_spatial_correlation_matrix": corr.tolist(),
        "per_step_growth_factor": per_step.tolist(),
        "map_storage": "all 73 hourly times on the complete 721 x 1440 grid",
        "spectrum_storage": "all 73 times x T719 total wavenumbers",
    }
    output = ROOT / config["paths"]["evaluation_dir"]
    output.mkdir(parents=True, exist_ok=True)
    (output / "metrics.json").write_text(json.dumps(report, indent=2) + "\n")
    names, hours = data["experiments"], data["times_hours"]
    figure, axes = plt.subplots(1, 3, figsize=(16, 4.5))
    for index, name in enumerate(names):
        axes[0].semilogy(hours, dke[index], label=str(name))
    axes[0].set(xlabel="Lead time (h)", ylabel="Global DKE (m2 s-2)", title="Ensemble DKE growth")
    axes[0].legend(fontsize=7)
    image = axes[1].imshow(maps[-1, map_index], extent=[0, 360, -90, 90], origin="lower", cmap="magma", aspect="auto")
    axes[1].set(xlabel="Longitude", ylabel="Latitude", title=f"{names[-1]} DKE at 72 h")
    figure.colorbar(image, ax=axes[1], label="m2 s-2")
    wave = data["total_wavenumber"]
    for index, name in enumerate(names):
        axes[2].loglog(wave[1:], spectra[index, -1, 1:], label=str(name))
    axes[2].set(xlabel="Total wavenumber", ylabel="Spectral DKE", title="T719 spectrum at 72 h")
    axes[2].legend(fontsize=7)
    figure.tight_layout()
    figure.savefig(output / "dke_diagnostics.png", dpi=160)
    plt.close(figure)
    print(f"evaluation={output.relative_to(ROOT)}")


if __name__ == "__main__":
    main()