Upload folder using huggingface_hub

.gitattributes

@@ -1,35 +1,3 @@
-*.7z filter=lfs diff=lfs merge=lfs -text
-*.arrow filter=lfs diff=lfs merge=lfs -text
-*.bin filter=lfs diff=lfs merge=lfs -text
-*.bz2 filter=lfs diff=lfs merge=lfs -text
-*.ckpt filter=lfs diff=lfs merge=lfs -text
-*.ftz filter=lfs diff=lfs merge=lfs -text
-*.gz filter=lfs diff=lfs merge=lfs -text
-*.h5 filter=lfs diff=lfs merge=lfs -text
-*.joblib filter=lfs diff=lfs merge=lfs -text
-*.lfs.* filter=lfs diff=lfs merge=lfs -text
-*.mlmodel filter=lfs diff=lfs merge=lfs -text
-*.model filter=lfs diff=lfs merge=lfs -text
-*.msgpack filter=lfs diff=lfs merge=lfs -text
-*.npy filter=lfs diff=lfs merge=lfs -text
-*.npz filter=lfs diff=lfs merge=lfs -text
-*.onnx filter=lfs diff=lfs merge=lfs -text
-*.ot filter=lfs diff=lfs merge=lfs -text
-*.parquet filter=lfs diff=lfs merge=lfs -text
-*.pb filter=lfs diff=lfs merge=lfs -text
-*.pickle filter=lfs diff=lfs merge=lfs -text
-*.pkl filter=lfs diff=lfs merge=lfs -text
 *.pt filter=lfs diff=lfs merge=lfs -text
-*.pth filter=lfs diff=lfs merge=lfs -text
-*.rar filter=lfs diff=lfs merge=lfs -text
-*.safetensors filter=lfs diff=lfs merge=lfs -text
-saved_model/**/* filter=lfs diff=lfs merge=lfs -text
-*.tar.* filter=lfs diff=lfs merge=lfs -text
-*.tar filter=lfs diff=lfs merge=lfs -text
-*.tflite filter=lfs diff=lfs merge=lfs -text
-*.tgz filter=lfs diff=lfs merge=lfs -text
-*.wasm filter=lfs diff=lfs merge=lfs -text
-*.xz filter=lfs diff=lfs merge=lfs -text
-*.zip filter=lfs diff=lfs merge=lfs -text
-*.zst filter=lfs diff=lfs merge=lfs -text
-*tfevents* filter=lfs diff=lfs merge=lfs -text
+*.npz filter=lfs diff=lfs merge=lfs -text
+*.npy filter=lfs diff=lfs merge=lfs -text

README.md

@@ -0,0 +1,149 @@
+---
+license: apache-2.0
+language:
+- en
+tags:
+- OneScience
+- Earth Science
+- Ensemble Forecast Postprocessing
+- Probabilistic Neural Network
+frameworks: PyTorch
+---
+
+<p align="center">
+  <strong><span style="font-size: 30px;">PPNN</span></strong>
+</p>
+
+# Model Introduction
+
+PPNN converts ensemble weather forecasts and station information into calibrated probability distributions for 2 m temperature, supporting station-level ensemble postprocessing, uncertainty representation, and probabilistic weather forecasting.
+
+Paper: Neural Networks for Postprocessing Ensemble Weather Forecasts  
+https://doi.org/10.1175/MWR-D-18-0187.1
+
+# Model Description
+
+The method was proposed by research teams at Ludwig-Maximilians-Universität München, Karlsruhe Institute of Technology, and the Heidelberg Institute for Theoretical Studies. The paper trains on ECMWF TIGGE ensemble forecasts and Deutscher Wetterdienst station observations of 2 m temperature. It is intended for station-temperature probabilistic postprocessing and ensemble calibration at a fixed 48-hour lead.
+
+# Use Cases
+
+| Use Case | Description |
+| :---: | :--- |
+| Ensemble forecast postprocessing | Generate Gaussian temperature forecasts from ensemble moments and station information. |
+| Probabilistic training | Train mean and scale parameters with Gaussian CRPS. |
+| ModelScope/OneCode execution | Validate feature construction, training, inference, calibration evaluation, and visualization in ModelScope or OneCode environments. |
+| Multi-GPU training | Train independent network replicas through `torchrun`. |
+
+# Usage Instructions
+
+## 1.OneCode
+
+Experience intelligent, one-click AI4S programming through the OneCode online environment:
+
+[Try intelligent, one-click AI4S programming](https://web-2069360198568017922-iaaj.ksai.scnet.cn:58043/home)
+
+## 2. Download and Installation
+
+```bash
+hf download OneScience-Group/PPNN --local-dir ./PPNN
+cd PPNN
+```
+
+### Environment Dependencies
+
+**Hardware Requirements**
+
+- A GPU or DCU is recommended.
+- A CPU can be used for connectivity validation with the default small-sample configuration.
+- DCU users must install DTK first. DTK 25.04.2 or later, or the OneScience-recommended version matching the cluster, is recommended.
+
+**DCU Environment**
+
+```bash
+# Activate DTK and Conda first
+conda create -n onescience311 python=3.11 -y
+conda activate onescience311
+pip install onescience[earth-dcu] -i http://mirrors.onescience.ai:3141/pypi/simple/ --trusted-host mirrors.onescience.ai
+```
+
+**GPU Environment**
+
+```bash
+# Activate Conda first
+conda create -n onescience311 python=3.11 -y libstdcxx-ng=12 libgcc-ng=12 gcc_linux-64=12 gxx_linux-64=12
+conda activate onescience311
+pip install onescience[earth-gpu] -i http://mirrors.onescience.ai:3141/pypi/simple/ --trusted-host mirrors.onescience.ai
+```
+
+### Training Data
+
+Training data consist of ECMWF TIGGE ensemble forecasts and DWD station observations of 2 m temperature. Each sample contains 50 ensemble members, 18 forecast variables, and station location and elevation information, with same-valid-time station temperature as the target. Forecasts are initialized daily at 00 UTC with a fixed 48-hour lead, and temperature is measured in degrees Celsius. This repository uses a small number of synthetic samples to validate training, inference, and probabilistic evaluation; they do not represent the official TIGGE or DWD distribution, training scale, or formal paper performance.
+
+```bash
+python scripts/fake_data.py
+```
+
+### Training
+
+For single-device training, use:
+
+```bash
+python scripts/train.py
+```
+
+For multi-GPU training, use:
+
+```bash
+torchrun --nproc_per_node=8 --nnodes=1 --rdzv_id=1000 --rdzv_backend=c10d --max_restarts=0 --master_addr="localhost" --master_port=29500 scripts/train.py
+```
+
+The default engineering configuration reduces the hidden width from the paper's long-training value of 512 to 32, the independent network count from 10 to 2, and the training epochs and sample count. It does not reduce the 50 members, 18 variables, metadata for 537 stations, two-dimensional station embedding, or 48-hour lead. A formal experiment requires the real TIGGE and DWD data, the paper-scale network ensemble, and the complete training period.
+
+Training outputs are saved to:
+
+```text
+result/checkpoints/ppnn.pt
+result/training/metrics.json
+```
+
+### Trained Weights
+
+No weights are bundled under `weight/`, and no directly loadable pretrained checkpoint released by the paper authors was identified. Locally trained engineering checkpoints are written to `result/checkpoints/ppnn.pt` and must not be represented as official pretrained weights.
+
+### Inference
+
+```bash
+python scripts/inference.py
+```
+
+Inference loads the trained checkpoint and generates Gaussian temperature forecasts at the fixed 48-hour lead from 50-member ensemble statistics and station information. Complete numerical results are saved to:
+
+```text
+result/output/predictions.npz
+```
+
+### Evaluation and Visualization
+
+```bash
+python scripts/result.py
+```
+
+Evaluation follows the paper's probabilistic protocol for Gaussian CRPS, raw-ensemble CRPS, CRPSS, PIT calibration, and spread-error diagnostics, saving overall and station-sample aggregate results; the task has one fixed 48-hour lead and no classes, so no per-lead or per-class metrics are saved. It also generates PIT-calibration and probabilistic-error comparisons. Synthetic-data results validate the engineering workflow only and do not represent formal paper performance.
+
+```text
+result/evaluation/metrics.json
+result/evaluation/ppnn_evaluation.png
+```
+
+# Official OneScience Information
+
+| Platform | OneScience Main Repository | Skills Repository |
+| --- | --- | --- |
+| Gitee | https://gitee.com/onescience-ai/onescience | https://gitee.com/onescience-ai/oneskills |
+| GitHub | https://github.com/onescience-ai/OneScience | https://github.com/onescience-ai/oneskills |
+
+# Citation and License
+
+This repository is an independent engineering reproduction of the public PPNN specifications.
+
+Use of this repository's code, official model weights, and data remains subject to the licenses and terms of their respective projects.

conf/config.yaml

@@ -0,0 +1,34 @@
+seed: 42
+data:
+  root: data
+  file: data/ppnn_virtual.npz
+  format_version: ppnn_nn_aux_emb_v1
+  dates: [2020-01-15, 2020-07-15, 2021-01-15, 2021-07-15]
+  active_station_count: 12
+  station_count: 537
+  member_count: 50
+  variable_count: 18
+  lead_hours: 48
+  target_variable: t2m
+model:
+  name: PPNN
+  variant: NN-aux-emb
+  hidden_size: 32
+  paper_hidden_size: 512
+  embedding_dim: 2
+  continuous_features: 40
+  input_features: 42
+  sigma_epsilon: 0.0001
+train:
+  replicas: 2
+  paper_replicas: 10
+  epochs: 12
+  batch_size: 64
+  learning_rate: 0.003
+runtime:
+  device: auto
+paths:
+  checkpoint: result/checkpoints/ppnn.pt
+  training_metrics: result/training/metrics.json
+  inference: result/output/predictions.npz
+  evaluation_dir: result/evaluation

config.json

@@ -0,0 +1,71 @@
+{
+  "model_name": "PPNN",
+  "model_type": "ppnn",
+  "architectures": [
+    "PPNN"
+  ],
+  "framework": "PyTorch",
+  "domain": "weather forecasting",
+  "task": "probabilistic ensemble forecast postprocessing",
+  "implementation": {
+    "entry_point": "model/ppnn.py",
+    "scope": "Core-method and full-data-dimension reduced-scale end-to-end engineering reproduction",
+    "train_script": "scripts/train.py",
+    "inference_script": "scripts/inference.py",
+    "evaluation_script": "scripts/result.py",
+    "synthetic_data_script": "scripts/fake_data.py"
+  },
+  "architecture": {
+    "variant": "NN-aux-emb",
+    "continuous_features": 40,
+    "station_embedding_dim": 2,
+    "network_input_features": 42,
+    "engineering_hidden_size": 32,
+    "paper_long_training_hidden_size": 512,
+    "activation": "ReLU",
+    "outputs": [
+      "mu",
+      "raw_sigma"
+    ],
+    "scale_transform": "abs(raw_sigma) + epsilon"
+  },
+  "data": {
+    "datasets": [
+      "ECMWF TIGGE ensemble forecasts",
+      "DWD station observations"
+    ],
+    "protocol": "Daily 00 UTC initialization at a fixed 48-hour lead; one sample is a valid date-station pair",
+    "format_version": "ppnn_nn_aux_emb_v1",
+    "raw_ensemble_shape": [
+      "N",
+      50,
+      18
+    ],
+    "continuous_input_shape": [
+      "N",
+      40
+    ],
+    "output_shape": [
+      "N",
+      2
+    ],
+    "ensemble_members": 50,
+    "forecast_variables": 18,
+    "station_count": 537,
+    "lead_time_hours": 48,
+    "target": "2 m temperature",
+    "target_unit": "degree Celsius",
+    "ensemble_statistics": [
+      "mean",
+      "sample standard deviation (ddof=1)"
+    ]
+  },
+  "configuration_sources": [
+    "conf/config.yaml",
+    "model/ppnn.py",
+    "scripts/fake_data.py",
+    "scripts/train.py",
+    "scripts/inference.py",
+    "scripts/result.py"
+  ]
+}

configuration.json

@@ -0,0 +1,12 @@
+{
+  "framework": "PyTorch",
+  "task": "ensemble_weather_forecast_postprocessing",
+  "model": "PPNN",
+  "input_format": "ensemble [N,50,18], auxiliary [N,4], station_index [N]",
+  "protocol": "NN-aux-emb Gaussian postprocessing at a fixed 48-hour lead with 537-station metadata",
+  "default_config": "conf/config.yaml",
+  "training": "scripts/train.py",
+  "inference": "scripts/inference.py",
+  "evaluation": "scripts/result.py",
+  "visualization": "scripts/result.py"
+}

model/ppnn.py

@@ -0,0 +1,74 @@
+"""PPNN NN-aux-emb model and proper scoring rules."""
+
+from __future__ import annotations
+
+import math
+
+import torch
+from torch import nn
+
+
+FORMAT_VERSION = "ppnn_nn_aux_emb_v1"
+MEMBER_COUNT = 50
+VARIABLE_COUNT = 18
+STATION_COUNT = 537
+LEAD_HOURS = 48
+AUX_FEATURE_COUNT = 4
+CONTINUOUS_FEATURE_COUNT = VARIABLE_COUNT * 2 + AUX_FEATURE_COUNT
+EMBEDDING_DIM = 2
+INPUT_FEATURE_COUNT = CONTINUOUS_FEATURE_COUNT + EMBEDDING_DIM
+
+
+class PPNN(nn.Module):
+    """NN-aux-emb: ensemble moments, auxiliary predictors, and station embedding."""
+
+    def __init__(self, hidden_size: int = 32, station_count: int = STATION_COUNT, eps: float = 1e-4):
+        super().__init__()
+        if station_count != STATION_COUNT:
+            raise ValueError(f"station_count must be the paper value {STATION_COUNT}")
+        self.station_count = station_count
+        self.eps = float(eps)
+        self.station_embedding = nn.Embedding(station_count, EMBEDDING_DIM)
+        self.hidden = nn.Linear(INPUT_FEATURE_COUNT, hidden_size)
+        self.output = nn.Linear(hidden_size, 2)
+
+    def forward(self, continuous: torch.Tensor, station_index: torch.Tensor):
+        if continuous.ndim != 2 or continuous.shape[-1] != CONTINUOUS_FEATURE_COUNT:
+            raise ValueError(f"continuous features must have shape [batch,{CONTINUOUS_FEATURE_COUNT}]")
+        if station_index.ndim != 1 or station_index.shape[0] != continuous.shape[0]:
+            raise ValueError("station_index must have shape [batch]")
+        if station_index.numel() and (station_index.min() < 0 or station_index.max() >= self.station_count):
+            raise ValueError("station index outside [0, 536]")
+        features = torch.cat((continuous, self.station_embedding(station_index)), dim=-1)
+        mu, raw_sigma = self.output(torch.relu(self.hidden(features))).unbind(dim=-1)
+        return mu, raw_sigma.abs() + self.eps
+
+
+def ensemble_features(ensemble: torch.Tensor, auxiliary: torch.Tensor) -> torch.Tensor:
+    """Reduce [batch, 50, 18] to 36 moments and append four auxiliary features."""
+    if ensemble.ndim != 3 or ensemble.shape[1:] != (MEMBER_COUNT, VARIABLE_COUNT):
+        raise ValueError(f"ensemble must have shape [batch,{MEMBER_COUNT},{VARIABLE_COUNT}]")
+    if auxiliary.ndim != 2 or auxiliary.shape != (ensemble.shape[0], AUX_FEATURE_COUNT):
+        raise ValueError(f"auxiliary must have shape [batch,{AUX_FEATURE_COUNT}]")
+    moments = torch.cat((ensemble.mean(dim=1), ensemble.std(dim=1, correction=1)), dim=-1)
+    return torch.cat((moments, auxiliary), dim=-1)
+
+
+def gaussian_crps(mu: torch.Tensor, sigma: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
+    """Closed-form CRPS for a Gaussian predictive distribution."""
+    if torch.any(sigma <= 0):
+        raise ValueError("sigma must be strictly positive")
+    z = (target - mu) / sigma
+    pdf = torch.exp(-0.5 * z.square()) / math.sqrt(2.0 * math.pi)
+    cdf = 0.5 * (1.0 + torch.erf(z / math.sqrt(2.0)))
+    return sigma * (z * (2.0 * cdf - 1.0) + 2.0 * pdf - 1.0 / math.sqrt(math.pi))
+
+
+def ensemble_crps(ensemble: torch.Tensor, target: torch.Tensor) -> torch.Tensor:
+    """Empirical ensemble CRPS without constructing a member-pair matrix."""
+    if ensemble.shape[-1] != MEMBER_COUNT:
+        raise ValueError(f"last ensemble dimension must be {MEMBER_COUNT}")
+    sorted_members = ensemble.sort(dim=-1).values
+    weights = (2 * torch.arange(1, MEMBER_COUNT + 1, device=ensemble.device) - MEMBER_COUNT - 1).to(ensemble.dtype)
+    pair_term = (sorted_members * weights).sum(dim=-1) / (MEMBER_COUNT * MEMBER_COUNT)
+    return (ensemble - target.unsqueeze(-1)).abs().mean(dim=-1) - pair_term

scripts/fake_data.py

@@ -0,0 +1,58 @@
+"""Generate a small, structured PPNN dataset while preserving paper dimensions."""
+
+from pathlib import Path
+
+import numpy as np
+import yaml
+
+ROOT = Path(__file__).resolve().parents[1]
+
+
+def main():
+    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())
+    data = config["data"]
+    if (data["station_count"], data["member_count"], data["variable_count"], data["lead_hours"]) != (537, 50, 18, 48):
+        raise ValueError("paper station/member/variable/lead dimensions cannot be reduced")
+    rng = np.random.default_rng(int(config["seed"]))
+    station_index_all = np.arange(537)
+    latitude = (47.0 + 8.0 * ((station_index_all * 0.61803398875) % 1)).astype(np.float32)
+    longitude = (5.5 + 9.5 * ((station_index_all * 0.41421356237) % 1)).astype(np.float32)
+    elevation = (5.0 + 1450.0 * ((station_index_all * 0.2718281828) % 1)).astype(np.float32)
+    station_id = np.asarray([f"DWD-{i:04d}" for i in station_index_all])
+    active = np.linspace(0, 536, int(data["active_station_count"]), dtype=np.int64)
+    dates = np.asarray(data["dates"], dtype="datetime64[D]")
+    date_id, station_id_index = np.meshgrid(np.arange(len(dates)), active, indexing="ij")
+    date_id, station_id_index = date_id.ravel(), station_id_index.ravel()
+    n = len(date_id)
+    day_of_year = (dates[date_id] - dates[date_id].astype("datetime64[Y]")).astype(int) + 1
+    seasonal = np.sin(2 * np.pi * (day_of_year - 172) / 365.25).astype(np.float32)
+    lat = latitude[station_id_index]
+    lon = longitude[station_id_index]
+    elev = elevation[station_id_index]
+    spatial = (-0.42 * (lat - 51.0) + 0.13 * (lon - 10.0) - 0.0062 * elev).astype(np.float32)
+    common_weather = rng.normal(0, 2.2, len(dates)).astype(np.float32)[date_id]
+    target = (11.0 + 9.0 * seasonal + spatial + common_weather + 0.18 * np.sin(station_id_index * 0.17)).astype(np.float32)
+    variable_offset = np.linspace(-4.0, 5.0, 18, dtype=np.float32)
+    variable_scale = np.linspace(0.72, 1.25, 18, dtype=np.float32)
+    centre = target[:, None] * variable_scale[None] + variable_offset[None] + 0.3 * common_weather[:, None]
+    spread = (0.45 + 0.025 * np.abs(common_weather) + 0.00045 * elev + 0.08 * (1.0 + seasonal)).astype(np.float32)
+    member_axis = np.linspace(-2.1, 2.1, 50, dtype=np.float32)
+    wave = np.sin((station_id_index[:, None, None] + 1) * 0.13 + np.arange(18)[None, None] * 0.31 + member_axis[None, :, None] * 1.7)
+    ensemble = centre[:, None, :] + spread[:, None, None] * variable_scale[None, None, :] * (member_axis[None, :, None] + 0.22 * wave)
+    ensemble[:, :, 0] += (0.35 + 0.14 * seasonal)[:, None]
+    auxiliary = np.stack((latitude[station_id_index], longitude[station_id_index], elevation[station_id_index], day_of_year / 366.0), axis=1).astype(np.float32)
+    output = ROOT / data["file"]
+    output.parent.mkdir(parents=True, exist_ok=True)
+    np.savez_compressed(
+        output, format_version=np.asarray(data["format_version"]), ensemble=ensemble.astype(np.float32),
+        auxiliary=auxiliary, target=target, station_index=station_id_index.astype(np.int64), date_index=date_id.astype(np.int64),
+        dates=dates.astype("U10"), lead_hours=np.asarray(data["lead_hours"], np.int64), variables=np.asarray(["t2m"] + [f"var_{i:02d}" for i in range(1, 18)]),
+        station_id=station_id, station_latitude=latitude, station_longitude=longitude, station_elevation_m=elevation,
+    )
+    if ensemble.shape != (n, 50, 18) or len(station_id) != 537 or not np.isfinite(ensemble).all():
+        raise RuntimeError("generated data failed shape or finite-value validation")
+    print(f"saved={output.relative_to(ROOT)} samples={n} ensemble={ensemble.shape} metadata_stations=537 lead=48h")
+
+
+if __name__ == "__main__":
+    main()

scripts/inference.py

@@ -0,0 +1,52 @@
+"""Run all PPNN replicas and save the complete verification payload."""
+
+import sys
+from pathlib import Path
+
+import numpy as np
+import torch
+import yaml
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+from model.ppnn import FORMAT_VERSION, PPNN, ensemble_features
+
+
+def main():
+    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())
+    data = np.load(ROOT / config["data"]["file"])
+    checkpoint = torch.load(ROOT / config["paths"]["checkpoint"], map_location="cpu", weights_only=False)
+    if str(data["format_version"]) != FORMAT_VERSION or checkpoint.get("format_version") != FORMAT_VERSION:
+        raise ValueError("data/checkpoint format version mismatch")
+    if data["ensemble"].shape[1:] != (50, 18) or len(data["station_id"]) != 537 or int(data["lead_hours"]) != 48:
+        raise ValueError("inference requires 50 members, 18 variables, 537 stations, and 48h lead")
+    scaling = checkpoint["scaling"]
+    continuous = ensemble_features(torch.from_numpy(data["ensemble"]), torch.from_numpy(data["auxiliary"]))
+    continuous = (continuous - scaling["feature_mean"]) / scaling["feature_std"]
+    station = torch.from_numpy(data["station_index"])
+    mus, sigmas = [], []
+    for state in checkpoint["model"]:
+        model = PPNN(int(checkpoint["model_config"]["hidden_size"]), eps=float(checkpoint["model_config"]["sigma_epsilon"]))
+        model.load_state_dict(state)
+        model.eval()
+        with torch.no_grad():
+            mu_scaled, sigma_scaled = model(continuous, station)
+        mus.append(mu_scaled * scaling["target_std"] + scaling["target_mean"])
+        sigmas.append(sigma_scaled * scaling["target_std"])
+    # The paper averages Gaussian parameters from random networks; this is not a mixture distribution.
+    mu = torch.stack(mus).mean(0).numpy()
+    sigma = torch.stack(sigmas).mean(0).numpy()
+    target = data["target"]
+    raw_t2m = data["ensemble"][:, :, 0]
+    if not all(np.isfinite(x).all() for x in (mu, sigma, target, raw_t2m)) or not np.all(sigma > 0):
+        raise RuntimeError("inference output is non-finite or has non-positive sigma")
+    output = ROOT / config["paths"]["inference"]
+    output.parent.mkdir(parents=True, exist_ok=True)
+    np.savez_compressed(output, mu=mu, sigma=sigma, targets=target, raw_ensemble_t2m=raw_t2m,
+                        station_index=data["station_index"], station_id=data["station_id"], date_index=data["date_index"],
+                        dates=data["dates"], lead_hours=data["lead_hours"], format_version=np.asarray(FORMAT_VERSION))
+    print(f"saved={output.relative_to(ROOT)} samples={len(mu)} replicas={len(mus)} raw_ensemble={raw_t2m.shape}")
+
+
+if __name__ == "__main__":
+    main()

scripts/result.py

@@ -0,0 +1,68 @@
+"""Evaluate Gaussian and raw-ensemble forecasts and create calibration plots."""
+
+import json
+import math
+import sys
+from pathlib import Path
+
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+import torch
+import yaml
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+from model.ppnn import FORMAT_VERSION, ensemble_crps, gaussian_crps
+
+
+def main():
+    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())
+    data = np.load(ROOT / config["paths"]["inference"])
+    if str(data["format_version"]) != FORMAT_VERSION or data["raw_ensemble_t2m"].shape != (len(data["mu"]), 50):
+        raise ValueError("prediction version/shape mismatch")
+    mu = torch.from_numpy(data["mu"])
+    sigma = torch.from_numpy(data["sigma"])
+    target = torch.from_numpy(data["targets"])
+    raw = torch.from_numpy(data["raw_ensemble_t2m"])
+    model_scores = gaussian_crps(mu, sigma, target)
+    raw_scores = ensemble_crps(raw, target)
+    model_mean, raw_mean = float(model_scores.mean()), float(raw_scores.mean())
+    pit = 0.5 * (1.0 + torch.erf((target - mu) / sigma / math.sqrt(2.0))).numpy()
+    bins = np.linspace(0.0, 1.0, 11)
+    counts, _ = np.histogram(pit, bins=bins)
+    observed = np.asarray([(pit <= edge).mean() for edge in bins[1:-1]])
+    nominal = bins[1:-1]
+    calibration_mae = float(np.mean(np.abs(observed - nominal)))
+    rmse = float(torch.sqrt(torch.mean((mu - target).square())))
+    spread_error_ratio = float(sigma.mean()) / rmse
+    metrics = {
+        "note": "structured virtual-data engineering validation, not paper performance",
+        "samples": len(mu), "station_metadata_count": len(data["station_id"]), "lead_hours": int(data["lead_hours"]),
+        "mean_gaussian_crps": model_mean, "mean_raw_ensemble_crps": raw_mean,
+        "crpss": 1.0 - model_mean / raw_mean, "pit_histogram_counts": counts.tolist(),
+        "pit_bin_edges": bins.tolist(), "calibration_mae": calibration_mae,
+        "rmse_mu": rmse, "mean_sigma": float(sigma.mean()), "spread_error_ratio": spread_error_ratio,
+    }
+    numeric = [model_mean, raw_mean, metrics["crpss"], calibration_mae, rmse, metrics["mean_sigma"], spread_error_ratio]
+    if not np.isfinite(numeric).all() or raw_mean <= 0 or rmse <= 0:
+        raise RuntimeError("evaluation metrics are not finite and valid")
+    output = ROOT / config["paths"]["evaluation_dir"]
+    output.mkdir(parents=True, exist_ok=True)
+    (output / "metrics.json").write_text(json.dumps(metrics, indent=2) + "\n")
+    fig, axes = plt.subplots(1, 3, figsize=(13, 4), constrained_layout=True)
+    axes[0].bar((bins[:-1] + bins[1:]) / 2, counts / counts.sum(), width=0.085)
+    axes[0].axhline(0.1, color="black", linestyle="--", linewidth=1)
+    axes[0].set(xlabel="PIT", ylabel="relative frequency", title="PIT histogram")
+    axes[1].plot(nominal, observed, marker="o"); axes[1].plot([0, 1], [0, 1], "k--")
+    axes[1].set(xlabel="nominal CDF", ylabel="observed frequency", title="Calibration")
+    axes[2].scatter(sigma.numpy(), np.abs(mu.numpy() - target.numpy()), c=data["date_index"], cmap="viridis", s=22)
+    axes[2].set(xlabel="predictive sigma", ylabel="absolute error", title=f"Spread/error ratio = {spread_error_ratio:.2f}")
+    fig.savefig(output / "ppnn_evaluation.png", dpi=150)
+    plt.close(fig)
+    print(f"gaussian_crps={model_mean:.6f} raw_crps={raw_mean:.6f} crpss={metrics['crpss']:.6f} spread_error_ratio={spread_error_ratio:.6f}")
+
+
+if __name__ == "__main__":
+    main()

scripts/train.py

@@ -0,0 +1,104 @@
+"""Train independent PPNN replicas with Gaussian CRPS and optional rank allocation."""
+
+import json
+import os
+import sys
+from pathlib import Path
+
+import numpy as np
+import torch
+import yaml
+from torch.utils.data import DataLoader, TensorDataset
+
+ROOT = Path(__file__).resolve().parents[1]
+sys.path.insert(0, str(ROOT))
+from model.ppnn import FORMAT_VERSION, PPNN, ensemble_features, gaussian_crps
+
+
+def load_data(config):
+    data = np.load(ROOT / config["data"]["file"])
+    expected = (len(data["target"]), 50, 18)
+    if str(data["format_version"]) != FORMAT_VERSION or data["ensemble"].shape != expected:
+        raise ValueError(f"dataset version/shape mismatch: expected {FORMAT_VERSION} and {expected}")
+    if len(data["station_id"]) != 537 or int(data["lead_hours"]) != 48:
+        raise ValueError("dataset must retain 537 stations and 48h lead")
+    arrays = [data["ensemble"], data["auxiliary"], data["target"]]
+    if not all(np.isfinite(x).all() for x in arrays):
+        raise ValueError("dataset contains non-finite values")
+    return data
+
+
+def main():
+    config = yaml.safe_load((ROOT / "conf/config.yaml").read_text())
+    data = load_data(config)
+    distributed = int(os.environ.get("WORLD_SIZE", "1")) > 1
+    if distributed:
+        torch.distributed.init_process_group("gloo")
+    rank = torch.distributed.get_rank() if distributed else 0
+    world_size = torch.distributed.get_world_size() if distributed else 1
+    device_name = config["runtime"]["device"]
+    local_rank = int(os.environ.get("LOCAL_RANK", "0"))
+    enough_accelerators = torch.cuda.is_available() and torch.cuda.device_count() >= world_size
+    device = torch.device("cuda", local_rank) if device_name == "auto" and enough_accelerators else torch.device("cpu")
+    ensemble = torch.from_numpy(data["ensemble"])
+    auxiliary = torch.from_numpy(data["auxiliary"])
+    station = torch.from_numpy(data["station_index"])
+    target = torch.from_numpy(data["target"])
+    continuous_raw = ensemble_features(ensemble, auxiliary)
+    feature_mean = continuous_raw.mean(0)
+    feature_std = continuous_raw.std(0, correction=1).clamp_min(1e-6)
+    target_mean, target_std = target.mean(), target.std(correction=1).clamp_min(1e-6)
+    continuous = (continuous_raw - feature_mean) / feature_std
+    target_scaled = (target - target_mean) / target_std
+    dataset = TensorDataset(continuous, station, target_scaled)
+    local = []
+    replicas = int(config["train"]["replicas"])
+    for replica in range(replicas):
+        if replica % world_size != rank:
+            continue
+        seed = int(config["seed"]) + replica
+        torch.manual_seed(seed)
+        generator = torch.Generator().manual_seed(seed)
+        loader = DataLoader(dataset, batch_size=int(config["train"]["batch_size"]), shuffle=True, generator=generator)
+        model = PPNN(int(config["model"]["hidden_size"]), eps=float(config["model"]["sigma_epsilon"])).to(device)
+        optimizer = torch.optim.Adam(model.parameters(), lr=float(config["train"]["learning_rate"]))
+        history = []
+        for epoch in range(int(config["train"]["epochs"])):
+            total = 0.0
+            for x, s, y in loader:
+                mu, sigma = model(x.to(device), s.to(device))
+                loss = gaussian_crps(mu, sigma, y.to(device)).mean()
+                if not torch.isfinite(loss):
+                    raise RuntimeError("non-finite Gaussian CRPS loss")
+                optimizer.zero_grad(set_to_none=True)
+                loss.backward()
+                optimizer.step()
+                total += float(loss.detach()) * len(y)
+            history.append(total / len(dataset))
+        local.append({"replica": replica, "state_dict": {k: v.cpu() for k, v in model.state_dict().items()}, "loss": history})
+        print(f"rank={rank} replica={replica} final_scaled_crps={history[-1]:.6f}")
+    if distributed:
+        gathered = [None] * world_size
+        torch.distributed.all_gather_object(gathered, local)
+        trained = [item for rank_items in gathered for item in rank_items]
+    else:
+        trained = local
+    if rank == 0:
+        trained.sort(key=lambda x: x["replica"])
+        if [x["replica"] for x in trained] != list(range(replicas)):
+            raise RuntimeError("distributed ranks did not return every replica")
+        checkpoint = ROOT / config["paths"]["checkpoint"]
+        metrics = ROOT / config["paths"]["training_metrics"]
+        checkpoint.parent.mkdir(parents=True, exist_ok=True)
+        metrics.parent.mkdir(parents=True, exist_ok=True)
+        scaling = {"feature_mean": feature_mean, "feature_std": feature_std, "target_mean": target_mean, "target_std": target_std}
+        states = [x["state_dict"] for x in trained]
+        torch.save({"format_version": FORMAT_VERSION, "model": states,
+                    "model_config": config["model"], "config": config, "scaling": scaling}, checkpoint)
+        metrics.write_text(json.dumps({"replicas": [{"id": x["replica"], "loss": x["loss"]} for x in trained], "world_size": world_size}, indent=2) + "\n")
+    if distributed:
+        torch.distributed.destroy_process_group()
+
+
+if __name__ == "__main__":
+    main()