yzt15806542928 commited on
Commit
80cf062
·
verified ·
1 Parent(s): 434e30c

Upload folder using huggingface_hub

Browse files
README.md ADDED
@@ -0,0 +1,153 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ ---
2
+ frameworks: PyTorch
3
+ language:
4
+ - en
5
+ license: apache-2.0
6
+ tags:
7
+ - OneScience
8
+ - Earth Science
9
+ - Weather Forecast
10
+ - Masked Autoencoder
11
+ - ERA5
12
+ - W-MAE
13
+ tasks: []
14
+ datasets:
15
+ - OneScience/ERA5
16
+ ---
17
+ <p align="center">
18
+ <strong>
19
+ <span style="font-size: 30px;">W-MAE</span>
20
+ </strong>
21
+ </p>
22
+
23
+ # Model Introduction
24
+
25
+ W-MAE (Weather Masked AutoEncoder) is a pretraining model for multivariable weather forecasting.
26
+
27
+ Paper: W-MAE: Pre-trained weather model with masked autoencoder for multi-variable weather forecasting
28
+
29
+ https://arxiv.org/abs/2304.08754
30
+
31
+ # Model Description
32
+
33
+ W-MAE first learns spatial relationships among weather variables through masked reconstruction, then learns temporal dependencies by fine-tuning on a forecasting task.
34
+
35
+ # Use Cases
36
+
37
+ | Scenario | Description |
38
+ | :---: | :--- |
39
+ | Masked weather-field pretraining | Train the W-MAE reconstruction model with ERA5 HDF5 data that follows this project's protocol. |
40
+ | Local quick validation | Use synthetic HDF5 data to check loading, training, inference, and visualization of inference results. |
41
+ | ModelScope / OneCode execution | Download the standalone model package, configure data, and run the scripts directly. |
42
+ | Multi-GPU training | Launch PyTorch DistributedDataParallel with `torchrun`. |
43
+
44
+ # Usage Guide
45
+
46
+ ## 1. OneCode Usage
47
+
48
+ Experience intelligent one-click AI4S programming through the OneCode online environment:
49
+
50
+ [Click to Experience Intelligent One-Click AI4S Programming](https://web-2069360198568017922-iaaj.ksai.scnet.cn:58043/home)
51
+
52
+ ## 2. Manual Installation and Usage
53
+
54
+ **Hardware Requirements**
55
+
56
+ - A GPU or DCU is recommended.
57
+ - CPU can be used for import and small-scale connectivity verification; full training and inference will be slow.
58
+ - DCU users must install DTK in advance. DTK 25.04.2 or above, or the OneScience recommended version matching your cluster, is recommended.
59
+
60
+ ### Download the Model Package
61
+
62
+ ```bash
63
+ hf download OneScience-Group/W-MAE --local-dir ./W-MAE
64
+ cd W-MAE
65
+ ```
66
+
67
+ ### Install the Runtime Environment
68
+
69
+ **DCU Environment**
70
+
71
+ ```bash
72
+ # Please activate DTK and CONDA first
73
+ conda create -n onescience311 python=3.11 -y
74
+ conda activate onescience311
75
+ # uv installation is supported
76
+ pip install onescience[earth-dcu] -i http://mirrors.onescience.ai:3141/pypi/simple/ --trusted-host mirrors.onescience.ai
77
+ ```
78
+
79
+ **GPU Environment**
80
+ ```bash
81
+ # Please activate CONDA first
82
+ conda create -n onescience311 python=3.11 -y libstdcxx-ng=12 libgcc-ng=12 gcc_linux-64=12 gxx_linux-64=12
83
+ conda activate onescience311
84
+ # uv installation is supported
85
+ pip install onescience[earth-gpu] -i http://mirrors.onescience.ai:3141/pypi/simple/ --trusted-host mirrors.onescience.ai
86
+ ```
87
+
88
+ ### Training Data Introduction
89
+
90
+ The OneScience community provides an ERA5 data slice for training. Download it to the directory configured by `data.dataset_dir` in `conf/config.yaml`:
91
+
92
+ ```bash
93
+ hf download --repo-type dataset OneScience-Group/ERA5 --local-dir ./data/era5
94
+ ```
95
+
96
+ The W-MAE adapter expects yearly HDF5 files under `data/era5/data/`. Each file must contain a `fields` dataset, an ordered list of 20 channel names, six-hour time steps, and normalization statistics. Verify the physical ERA5 variable order before scientific training.
97
+
98
+ ### Generate Synthetic Data
99
+
100
+ When real ERA5 data is unavailable, generate protocol-compatible files for pipeline checks:
101
+
102
+ ```bash
103
+ python scripts/fake_data.py
104
+ ```
105
+
106
+ The synthetic files use placeholder channel names and must not be used for scientific evaluation.
107
+
108
+ ### Training
109
+
110
+ Single GPU:
111
+
112
+ ```bash
113
+ python scripts/train.py
114
+ ```
115
+
116
+ Multi-GPU:
117
+
118
+ ```bash
119
+ torchrun --nproc_per_node=8 scripts/train.py
120
+ ```
121
+
122
+ Training starts from random initialization and saves `data/checkpoint/model_bak.pth` by default. A compatible checkpoint can be supplied explicitly when continuing training.
123
+
124
+ ### Training Weights
125
+
126
+ This repository provides a `weight/` directory for W-MAE checkpoints. The weight files will be uploaded soon and are expected to be available in the near future.
127
+
128
+ ### Inference
129
+
130
+ Inference reads `data/checkpoint/model_bak.pth` by default and writes compressed reconstruction samples to `outputs/inference/`:
131
+
132
+ ```bash
133
+ python scripts/inference.py
134
+ ```
135
+
136
+ ### Evaluation and Visualization
137
+
138
+ ```bash
139
+ python scripts/result.py
140
+ ```
141
+
142
+ The result script validates reconstruction files and writes diagnostic figures under `outputs/inference/diagnostics/`.
143
+
144
+ # Official OneScience Resources
145
+
146
+ | Platform | OneScience Main Repository | Skills Repository |
147
+ | --- | --- | --- |
148
+ | Gitee | https://gitee.com/onescience-ai/onescience | https://gitee.com/onescience-ai/oneskills |
149
+ | GitHub | https://github.com/onescience-ai/OneScience | https://github.com/onescience-ai/oneskills |
150
+
151
+ # Citation and License
152
+
153
+ - This repository is an independent reproduction of the original W-MAE paper.
conf/config.yaml ADDED
@@ -0,0 +1,89 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ project:
2
+ name: W-MAE
3
+ output_dir: ./outputs
4
+
5
+ data:
6
+ dataset_dir: ./data/era5
7
+ train_years: [2000, 2001]
8
+ val_years: [2002]
9
+ test_years: [2003]
10
+ # These names match the generated test HDF5 files. Replace the complete list
11
+ # with the verified channel order when using real ERA5 data.
12
+ variables:
13
+ - synthetic_channel_00
14
+ - synthetic_channel_01
15
+ - synthetic_channel_02
16
+ - synthetic_channel_03
17
+ - synthetic_channel_04
18
+ - synthetic_channel_05
19
+ - synthetic_channel_06
20
+ - synthetic_channel_07
21
+ - synthetic_channel_08
22
+ - synthetic_channel_09
23
+ - synthetic_channel_10
24
+ - synthetic_channel_11
25
+ - synthetic_channel_12
26
+ - synthetic_channel_13
27
+ - synthetic_channel_14
28
+ - synthetic_channel_15
29
+ - synthetic_channel_16
30
+ - synthetic_channel_17
31
+ - synthetic_channel_18
32
+ - synthetic_channel_19
33
+ expected_channels: 20
34
+ input_steps: 1
35
+ output_steps: 1
36
+ time_step_hours: 6
37
+ normalize: true
38
+ source_size: [721, 1440]
39
+ model_size: [720, 1440]
40
+ latitude_crop: drop_last
41
+ batch_size: 4
42
+ num_workers: 0
43
+
44
+ fake_data:
45
+ years: [2000, 2001, 2002, 2003]
46
+ time_steps: 40
47
+ variables:
48
+ - synthetic_channel_00
49
+ - synthetic_channel_01
50
+ - synthetic_channel_02
51
+ - synthetic_channel_03
52
+ - synthetic_channel_04
53
+ - synthetic_channel_05
54
+ - synthetic_channel_06
55
+ - synthetic_channel_07
56
+ - synthetic_channel_08
57
+ - synthetic_channel_09
58
+ - synthetic_channel_10
59
+ - synthetic_channel_11
60
+ - synthetic_channel_12
61
+ - synthetic_channel_13
62
+ - synthetic_channel_14
63
+ - synthetic_channel_15
64
+ - synthetic_channel_16
65
+ - synthetic_channel_17
66
+ - synthetic_channel_18
67
+ - synthetic_channel_19
68
+
69
+ model:
70
+ image_size: [720, 1440]
71
+ patch_size: [8, 8]
72
+ input_channels: 20
73
+ embed_dim: 768
74
+ encoder_depth: 12
75
+ decoder_embed_dim: 512
76
+ decoder_depth: 6
77
+ mask_ratio: 0.75
78
+ norm_pixel_loss: false
79
+
80
+ inference:
81
+ output_dir: ./outputs/inference
82
+ split: test
83
+ samples: 1
84
+ mask_ratio: 0.75
85
+
86
+ validation:
87
+ expected_channels: 20
88
+ mask_ratio: 0.75
89
+ diagnostic_channel: 0
config.json ADDED
@@ -0,0 +1,85 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "model_name": "W-MAE",
3
+ "model_type": "w_mae",
4
+ "architectures": [
5
+ "MaskedAutoencoderAFNO"
6
+ ],
7
+ "framework": "PyTorch",
8
+ "domain": "atmosphere",
9
+ "task": "masked-weather-pretraining",
10
+ "implementation": {
11
+ "entry_point": "model/w_mae.py",
12
+ "scope": "AFNO-based masked autoencoder for multi-variable weather fields; this metadata describes the local PyTorch implementation and its ERA5 adapter"
13
+ },
14
+ "architecture": {
15
+ "family": "masked autoencoder with Adaptive Fourier Neural Operator (AFNO) blocks",
16
+ "input_format": "BCHW",
17
+ "input_grid_shape": [
18
+ 720,
19
+ 1440
20
+ ],
21
+ "patch_size": [
22
+ 8,
23
+ 8
24
+ ],
25
+ "input_channels": 20,
26
+ "encoder": {
27
+ "embed_dim": 768,
28
+ "depth": 12,
29
+ "mlp_ratio": 4.0,
30
+ "afno_num_blocks": 8
31
+ },
32
+ "decoder": {
33
+ "embed_dim": 512,
34
+ "depth": 6,
35
+ "mlp_ratio": 4.0
36
+ },
37
+ "mask_ratio": 0.75,
38
+ "supported_mask_ratios": [
39
+ 0.0,
40
+ 0.75
41
+ ],
42
+ "normalization": "LayerNorm",
43
+ "activation": "GELU"
44
+ },
45
+ "data": {
46
+ "dataset": "ERA5",
47
+ "source_grid_shape": [
48
+ 721,
49
+ 1440
50
+ ],
51
+ "preprocessing": "drop the last latitude row to obtain the model grid",
52
+ "time_step_hours": 6,
53
+ "input_steps": 1,
54
+ "output_steps": 1,
55
+ "channel_names": [
56
+ "synthetic_channel_00",
57
+ "synthetic_channel_01",
58
+ "synthetic_channel_02",
59
+ "synthetic_channel_03",
60
+ "synthetic_channel_04",
61
+ "synthetic_channel_05",
62
+ "synthetic_channel_06",
63
+ "synthetic_channel_07",
64
+ "synthetic_channel_08",
65
+ "synthetic_channel_09",
66
+ "synthetic_channel_10",
67
+ "synthetic_channel_11",
68
+ "synthetic_channel_12",
69
+ "synthetic_channel_13",
70
+ "synthetic_channel_14",
71
+ "synthetic_channel_15",
72
+ "synthetic_channel_16",
73
+ "synthetic_channel_17",
74
+ "synthetic_channel_18",
75
+ "synthetic_channel_19"
76
+ ],
77
+ "channel_names_note": "The repository configuration intentionally uses synthetic names; physical ERA5 variable names must be supplied and verified before scientific training."
78
+ },
79
+ "configuration_sources": [
80
+ "conf/config.yaml",
81
+ "model/w_mae.py",
82
+ "model/afno.py",
83
+ "scripts/era5_adapter.py"
84
+ ]
85
+ }
configuration.json ADDED
@@ -0,0 +1,12 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ {
2
+ "framework": "PyTorch",
3
+ "task": "masked_weather_pretraining",
4
+ "model": "W-MAE",
5
+ "input_format": "BCHW",
6
+ "protocol": "era5_20_channel_720x1440",
7
+ "default_config": "conf/config.yaml",
8
+ "train": "scripts/train.py",
9
+ "inference": "scripts/inference.py",
10
+ "evaluation": "scripts/result.py",
11
+ "visualization": "scripts/result.py"
12
+ }
model/afno.py ADDED
@@ -0,0 +1,151 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import torch
4
+ from torch import nn
5
+ from torch.nn import functional as F
6
+
7
+
8
+ class Mlp(nn.Module):
9
+ def __init__(
10
+ self,
11
+ in_features: int,
12
+ hidden_features: int | None = None,
13
+ out_features: int | None = None,
14
+ act_layer: type[nn.Module] = nn.GELU,
15
+ drop: float = 0.0,
16
+ ) -> None:
17
+ super().__init__()
18
+ hidden_features = hidden_features or in_features
19
+ out_features = out_features or in_features
20
+ self.fc1 = nn.Linear(in_features, hidden_features)
21
+ self.act = act_layer()
22
+ self.fc2 = nn.Linear(hidden_features, out_features)
23
+ self.drop = nn.Dropout(drop)
24
+
25
+ def forward(self, x: torch.Tensor) -> torch.Tensor:
26
+ x = self.fc1(x)
27
+ x = self.act(x)
28
+ x = self.drop(x)
29
+ x = self.fc2(x)
30
+ return self.drop(x)
31
+
32
+
33
+ class AFNO2D(nn.Module):
34
+ """Adaptive Fourier mixing over a channel-last 2D patch grid."""
35
+
36
+ def __init__(
37
+ self,
38
+ hidden_size: int,
39
+ num_blocks: int = 8,
40
+ sparsity_threshold: float = 0.01,
41
+ hard_thresholding_fraction: float = 1.0,
42
+ hidden_size_factor: int = 1,
43
+ ) -> None:
44
+ super().__init__()
45
+ if hidden_size % num_blocks != 0:
46
+ raise ValueError(f"hidden_size={hidden_size} must be divisible by num_blocks={num_blocks}.")
47
+ if not 0.0 < hard_thresholding_fraction <= 1.0:
48
+ raise ValueError("hard_thresholding_fraction must be in (0, 1].")
49
+
50
+ self.hidden_size = hidden_size
51
+ self.sparsity_threshold = sparsity_threshold
52
+ self.num_blocks = num_blocks
53
+ self.block_size = hidden_size // num_blocks
54
+ self.hard_thresholding_fraction = hard_thresholding_fraction
55
+ self.hidden_size_factor = hidden_size_factor
56
+ scale = 0.02
57
+
58
+ self.w1 = nn.Parameter(
59
+ scale * torch.randn(2, num_blocks, self.block_size, self.block_size * hidden_size_factor)
60
+ )
61
+ self.b1 = nn.Parameter(scale * torch.randn(2, num_blocks, self.block_size * hidden_size_factor))
62
+ self.w2 = nn.Parameter(
63
+ scale * torch.randn(2, num_blocks, self.block_size * hidden_size_factor, self.block_size)
64
+ )
65
+ self.b2 = nn.Parameter(scale * torch.randn(2, num_blocks, self.block_size))
66
+
67
+ def forward(self, x: torch.Tensor) -> torch.Tensor:
68
+ if x.ndim != 4 or x.shape[-1] != self.hidden_size:
69
+ raise ValueError(f"AFNO2D expects [B,H,W,{self.hidden_size}], got {tuple(x.shape)}.")
70
+
71
+ bias = x
72
+ dtype = x.dtype
73
+ x = torch.fft.rfft2(x.float(), dim=(1, 2), norm="ortho")
74
+ batch, height, freq_width, _ = x.shape
75
+ x = x.reshape(batch, height, freq_width, self.num_blocks, self.block_size)
76
+
77
+ hidden_block = self.block_size * self.hidden_size_factor
78
+ o1_real = torch.zeros(
79
+ batch, height, freq_width, self.num_blocks, hidden_block, device=x.device, dtype=x.real.dtype
80
+ )
81
+ o1_imag = torch.zeros_like(o1_real)
82
+ o2_real = torch.zeros_like(x.real)
83
+ o2_imag = torch.zeros_like(x.real)
84
+
85
+ total_modes = height // 2 + 1
86
+ kept_modes = max(1, int(total_modes * self.hard_thresholding_fraction))
87
+ row_slice = slice(total_modes - kept_modes, total_modes + kept_modes)
88
+ col_slice = slice(0, min(kept_modes, freq_width))
89
+ selected = x[:, row_slice, col_slice]
90
+
91
+ o1_real[:, row_slice, col_slice] = F.relu(
92
+ torch.einsum("...bi,bio->...bo", selected.real, self.w1[0])
93
+ - torch.einsum("...bi,bio->...bo", selected.imag, self.w1[1])
94
+ + self.b1[0]
95
+ )
96
+ o1_imag[:, row_slice, col_slice] = F.relu(
97
+ torch.einsum("...bi,bio->...bo", selected.imag, self.w1[0])
98
+ + torch.einsum("...bi,bio->...bo", selected.real, self.w1[1])
99
+ + self.b1[1]
100
+ )
101
+
102
+ hidden_real = o1_real[:, row_slice, col_slice]
103
+ hidden_imag = o1_imag[:, row_slice, col_slice]
104
+ o2_real[:, row_slice, col_slice] = (
105
+ torch.einsum("...bi,bio->...bo", hidden_real, self.w2[0])
106
+ - torch.einsum("...bi,bio->...bo", hidden_imag, self.w2[1])
107
+ + self.b2[0]
108
+ )
109
+ o2_imag[:, row_slice, col_slice] = (
110
+ torch.einsum("...bi,bio->...bo", hidden_imag, self.w2[0])
111
+ + torch.einsum("...bi,bio->...bo", hidden_real, self.w2[1])
112
+ + self.b2[1]
113
+ )
114
+
115
+ x = torch.stack((o2_real, o2_imag), dim=-1)
116
+ x = F.softshrink(x, lambd=self.sparsity_threshold)
117
+ x = torch.view_as_complex(x)
118
+ x = x.reshape(batch, height, freq_width, self.hidden_size)
119
+ x = torch.fft.irfft2(x, s=bias.shape[1:3], dim=(1, 2), norm="ortho")
120
+ return x.to(dtype=dtype) + bias
121
+
122
+
123
+ class Block(nn.Module):
124
+ def __init__(
125
+ self,
126
+ dim: int,
127
+ mlp_ratio: float = 4.0,
128
+ drop: float = 0.0,
129
+ act_layer: type[nn.Module] = nn.GELU,
130
+ norm_layer: type[nn.Module] = nn.LayerNorm,
131
+ double_skip: bool = True,
132
+ num_blocks: int = 8,
133
+ sparsity_threshold: float = 0.01,
134
+ hard_thresholding_fraction: float = 1.0,
135
+ ) -> None:
136
+ super().__init__()
137
+ self.norm1 = norm_layer(dim)
138
+ self.filter = AFNO2D(dim, num_blocks, sparsity_threshold, hard_thresholding_fraction)
139
+ self.drop_path = nn.Identity()
140
+ self.norm2 = norm_layer(dim)
141
+ self.mlp = Mlp(dim, int(dim * mlp_ratio), act_layer=act_layer, drop=drop)
142
+ self.double_skip = double_skip
143
+
144
+ def forward(self, x: torch.Tensor) -> torch.Tensor:
145
+ residual = x
146
+ x = self.filter(self.norm1(x))
147
+ if self.double_skip:
148
+ x = x + residual
149
+ residual = x
150
+ x = self.drop_path(self.mlp(self.norm2(x)))
151
+ return x + residual
model/check_checkpoint.py ADDED
@@ -0,0 +1,117 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import argparse
4
+ import json
5
+ from pathlib import Path
6
+ import re
7
+ import subprocess
8
+ import sys
9
+ import zipfile
10
+
11
+
12
+ PROJECT_ROOT = Path(__file__).resolve().parents[1]
13
+ if str(PROJECT_ROOT) not in sys.path:
14
+ sys.path.insert(0, str(PROJECT_ROOT))
15
+
16
+
17
+ REQUIRED_KEYS = {
18
+ "cls_token",
19
+ "pos_embed",
20
+ "mask_token",
21
+ "decoder_pos_embed",
22
+ "patch_embed.proj.weight",
23
+ "decoder_embed.weight",
24
+ "decoder_norm.weight",
25
+ "decoder_pred.weight",
26
+ }
27
+
28
+ PARAMETER_KEY_PATTERN = re.compile(
29
+ r"(?:decoder_blocks|blocks)\.\d+\.[A-Za-z0-9_.]*?(?:weight|bias|w1|b1|w2|b2)"
30
+ r"|(?:patch_embed\.proj|decoder_embed|decoder_norm|decoder_pred)\.(?:weight|bias)"
31
+ r"|(?:cls_token|pos_embed|mask_token|decoder_pos_embed)"
32
+ )
33
+
34
+
35
+ def inspect_archive(path: Path) -> dict:
36
+ if not zipfile.is_zipfile(path):
37
+ raise ValueError(f"{path} is not a zip-based PyTorch checkpoint.")
38
+ with zipfile.ZipFile(path) as archive:
39
+ names = archive.namelist()
40
+ data_pickle = next((name for name in names if name.endswith("/data.pkl")), None)
41
+ if data_pickle is None:
42
+ raise ValueError("Checkpoint archive does not contain data.pkl.")
43
+ payload = archive.read(data_pickle)
44
+
45
+ strings = subprocess.run(
46
+ ["strings"], input=payload, stdout=subprocess.PIPE, stderr=subprocess.PIPE, check=True
47
+ ).stdout.decode("utf-8", errors="ignore").splitlines()
48
+ parameter_keys = {
49
+ match.group(0)
50
+ for line in strings
51
+ for match in [PARAMETER_KEY_PATTERN.search(line)]
52
+ if match
53
+ }
54
+ encoder_indices = sorted(
55
+ {int(match.group(1)) for key in parameter_keys if (match := re.match(r"blocks\.(\d+)\.", key))}
56
+ )
57
+ decoder_indices = sorted(
58
+ {int(match.group(1)) for key in parameter_keys if (match := re.match(r"decoder_blocks\.(\d+)\.", key))}
59
+ )
60
+ missing = sorted(REQUIRED_KEYS - parameter_keys)
61
+ return {
62
+ "path": str(path),
63
+ "archive_entries": len(names),
64
+ "parameter_key_count": len(parameter_keys),
65
+ "encoder_blocks": len(encoder_indices),
66
+ "decoder_blocks": len(decoder_indices),
67
+ "missing_required_keys": missing,
68
+ "status": "compatible_structure" if not missing and len(encoder_indices) == 12 and len(decoder_indices) == 6 else "mismatch",
69
+ }
70
+
71
+
72
+ def inspect_with_torch(path: Path) -> dict:
73
+ try:
74
+ import torch
75
+ from model.w_mae import w_mae_base
76
+ except (ImportError, OSError) as error:
77
+ return {"status": "unavailable", "reason": str(error)}
78
+
79
+ try:
80
+ checkpoint = torch.load(path, map_location="cpu")
81
+ except Exception as error:
82
+ return {"status": "unavailable", "reason": f"checkpoint deserialization failed: {error}"}
83
+ state_dict = checkpoint.get("model", checkpoint)
84
+ model = w_mae_base()
85
+ model_state = model.state_dict()
86
+ missing = sorted(set(model_state) - set(state_dict))
87
+ unexpected = sorted(set(state_dict) - set(model_state))
88
+ shape_mismatches = {
89
+ key: {"model": list(model_state[key].shape), "checkpoint": list(state_dict[key].shape)}
90
+ for key in model_state.keys() & state_dict.keys()
91
+ if tuple(model_state[key].shape) != tuple(state_dict[key].shape)
92
+ }
93
+ return {
94
+ "status": "compatible" if not missing and not unexpected and not shape_mismatches else "mismatch",
95
+ "missing_keys": missing,
96
+ "unexpected_keys": unexpected,
97
+ "shape_mismatches": shape_mismatches,
98
+ }
99
+
100
+
101
+ def main() -> None:
102
+ parser = argparse.ArgumentParser(description="Inspect W-MAE checkpoint structure and model compatibility.")
103
+ parser.add_argument("checkpoint", type=Path)
104
+ parser.add_argument("--torch-check", action="store_true")
105
+ args = parser.parse_args()
106
+
107
+ report = {"archive": inspect_archive(args.checkpoint)}
108
+ if args.torch_check:
109
+ report["torch"] = inspect_with_torch(args.checkpoint)
110
+ print(json.dumps(report, indent=2))
111
+ archive_ok = report["archive"]["status"] == "compatible_structure"
112
+ torch_ok = not args.torch_check or report["torch"]["status"] == "compatible"
113
+ raise SystemExit(0 if archive_ok and torch_ok else 1)
114
+
115
+
116
+ if __name__ == "__main__":
117
+ main()
model/checkpoint.py ADDED
@@ -0,0 +1,83 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ from contextlib import contextmanager
4
+ from pathlib import Path
5
+ import sys
6
+ from typing import Any, Iterator, Mapping
7
+
8
+ import torch
9
+ from torch import nn
10
+
11
+
12
+ @contextmanager
13
+ def _temporary_module_path(source_root: str | Path | None) -> Iterator[None]:
14
+ if source_root is None:
15
+ yield
16
+ return
17
+
18
+ path = str(Path(source_root).resolve())
19
+ sys.path.insert(0, path)
20
+ try:
21
+ yield
22
+ finally:
23
+ try:
24
+ sys.path.remove(path)
25
+ except ValueError:
26
+ pass
27
+
28
+
29
+ def _extract_state_dict(checkpoint: Any) -> Mapping[str, torch.Tensor]:
30
+ if isinstance(checkpoint, Mapping):
31
+ for key in ("model", "model_state", "state_dict"):
32
+ candidate = checkpoint.get(key)
33
+ if isinstance(candidate, Mapping):
34
+ return candidate
35
+ if checkpoint and all(isinstance(key, str) for key in checkpoint):
36
+ return checkpoint
37
+ raise ValueError("Checkpoint does not contain a model state dictionary.")
38
+
39
+
40
+ def load_checkpoint(
41
+ model: nn.Module,
42
+ checkpoint_path: str | Path,
43
+ *,
44
+ strict: bool = True,
45
+ map_location: str | torch.device = "cpu",
46
+ source_root: str | Path | None = None,
47
+ ) -> dict[str, Any]:
48
+ """Load a trusted W-MAE checkpoint and return loading metadata.
49
+
50
+ ``source_root`` is needed only for legacy checkpoints whose pickle payload
51
+ references modules from the original repository, such as ``utils.YParams``.
52
+ """
53
+ checkpoint_path = Path(checkpoint_path)
54
+ if not checkpoint_path.is_file():
55
+ raise FileNotFoundError(f"Checkpoint does not exist: {checkpoint_path}")
56
+
57
+ with _temporary_module_path(source_root):
58
+ checkpoint = torch.load(checkpoint_path, map_location=map_location, weights_only=False)
59
+
60
+ state_dict = _extract_state_dict(checkpoint)
61
+ model_state = model.state_dict()
62
+ missing = sorted(set(model_state) - set(state_dict))
63
+ unexpected = sorted(set(state_dict) - set(model_state))
64
+ shape_mismatches = {
65
+ key: {"model": list(model_state[key].shape), "checkpoint": list(state_dict[key].shape)}
66
+ for key in model_state.keys() & state_dict.keys()
67
+ if tuple(model_state[key].shape) != tuple(state_dict[key].shape)
68
+ }
69
+ if strict and (missing or unexpected or shape_mismatches):
70
+ raise RuntimeError(
71
+ "Checkpoint is incompatible with the model: "
72
+ f"missing={missing}, unexpected={unexpected}, shape_mismatches={shape_mismatches}"
73
+ )
74
+
75
+ incompatible = model.load_state_dict(state_dict, strict=strict)
76
+ return {
77
+ "checkpoint_path": str(checkpoint_path),
78
+ "strict": strict,
79
+ "missing_keys": list(incompatible.missing_keys),
80
+ "unexpected_keys": list(incompatible.unexpected_keys),
81
+ "shape_mismatches": shape_mismatches,
82
+ "checkpoint_keys": len(state_dict),
83
+ }
model/w_mae.py ADDED
@@ -0,0 +1,260 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ from functools import partial
4
+ from typing import NamedTuple
5
+
6
+ import numpy as np
7
+ import torch
8
+ from torch import nn
9
+
10
+ from .afno import Block
11
+
12
+
13
+ class WMAEOutput(NamedTuple):
14
+ loss: torch.Tensor
15
+ prediction: torch.Tensor
16
+ mask: torch.Tensor
17
+
18
+
19
+ def _sincos_1d(embed_dim: int, positions: np.ndarray) -> np.ndarray:
20
+ if embed_dim % 2 != 0:
21
+ raise ValueError("The 1D sine-cosine embedding dimension must be even.")
22
+ omega = np.arange(embed_dim // 2, dtype=np.float64)
23
+ omega = 1.0 / (10000 ** (omega / (embed_dim / 2.0)))
24
+ values = np.einsum("m,d->md", positions.reshape(-1), omega)
25
+ return np.concatenate((np.sin(values), np.cos(values)), axis=1)
26
+
27
+
28
+ def build_2d_sincos_position_embedding(
29
+ embed_dim: int, grid_size: tuple[int, int], include_cls_token: bool
30
+ ) -> torch.Tensor:
31
+ if embed_dim % 4 != 0:
32
+ raise ValueError("The 2D sine-cosine embedding dimension must be divisible by four.")
33
+ grid_h = np.arange(grid_size[0], dtype=np.float32)
34
+ grid_w = np.arange(grid_size[1], dtype=np.float32)
35
+ grid = np.meshgrid(grid_w, grid_h)
36
+ embedding = np.concatenate(
37
+ (_sincos_1d(embed_dim // 2, grid[0]), _sincos_1d(embed_dim // 2, grid[1])), axis=1
38
+ )
39
+ if include_cls_token:
40
+ embedding = np.concatenate((np.zeros((1, embed_dim)), embedding), axis=0)
41
+ return torch.from_numpy(embedding).float().unsqueeze(0)
42
+
43
+
44
+ class PatchEmbed(nn.Module):
45
+ def __init__(
46
+ self,
47
+ img_size: tuple[int, int],
48
+ patch_size: tuple[int, int],
49
+ in_chans: int,
50
+ embed_dim: int,
51
+ ) -> None:
52
+ super().__init__()
53
+ if img_size[0] % patch_size[0] or img_size[1] % patch_size[1]:
54
+ raise ValueError(f"img_size={img_size} must be divisible by patch_size={patch_size}.")
55
+ self.img_size = img_size
56
+ self.patch_size = patch_size
57
+ self.grid_size = (img_size[0] // patch_size[0], img_size[1] // patch_size[1])
58
+ self.num_patches = self.grid_size[0] * self.grid_size[1]
59
+ self.proj = nn.Conv2d(in_chans, embed_dim, kernel_size=patch_size, stride=patch_size)
60
+
61
+ def forward(self, x: torch.Tensor) -> torch.Tensor:
62
+ if x.ndim != 4 or tuple(x.shape[-2:]) != self.img_size:
63
+ raise ValueError(f"PatchEmbed expects [B,C,{self.img_size[0]},{self.img_size[1]}], got {tuple(x.shape)}.")
64
+ return self.proj(x).flatten(2).transpose(1, 2)
65
+
66
+
67
+ class MaskedAutoencoderAFNO(nn.Module):
68
+ """W-MAE pretraining model reconstructed from the official AFNO source."""
69
+
70
+ def __init__(
71
+ self,
72
+ img_size: tuple[int, int] = (720, 1440),
73
+ patch_size: tuple[int, int] = (8, 8),
74
+ in_chans: int = 20,
75
+ embed_dim: int = 768,
76
+ depth: int = 12,
77
+ decoder_embed_dim: int = 512,
78
+ decoder_depth: int = 6,
79
+ mlp_ratio: float = 4.0,
80
+ norm_layer: type[nn.Module] = nn.LayerNorm,
81
+ norm_pix_loss: bool = False,
82
+ num_blocks: int = 8,
83
+ sparsity_threshold: float = 0.01,
84
+ hard_thresholding_fraction: float = 1.0,
85
+ ) -> None:
86
+ super().__init__()
87
+ self.img_size = tuple(img_size)
88
+ self.patch_size = tuple(patch_size)
89
+ self.in_chans = in_chans
90
+ self.embed_dim = embed_dim
91
+ self.decoder_embed_dim = decoder_embed_dim
92
+ self.norm_pix_loss = norm_pix_loss
93
+
94
+ self.patch_embed = PatchEmbed(self.img_size, self.patch_size, in_chans, embed_dim)
95
+ num_patches = self.patch_embed.num_patches
96
+ # The official AFNO path keeps this checkpoint key but does not prepend
97
+ # a class token in forward_encoder, so it must not enter DDP reduction.
98
+ self.cls_token = nn.Parameter(torch.zeros(1, 1, embed_dim), requires_grad=False)
99
+ self.pos_embed = nn.Parameter(torch.zeros(1, num_patches + 1, embed_dim), requires_grad=False)
100
+ block_args = dict(
101
+ mlp_ratio=mlp_ratio,
102
+ norm_layer=norm_layer,
103
+ num_blocks=num_blocks,
104
+ sparsity_threshold=sparsity_threshold,
105
+ hard_thresholding_fraction=hard_thresholding_fraction,
106
+ )
107
+ self.blocks = nn.ModuleList([Block(dim=embed_dim, **block_args) for _ in range(depth)])
108
+ self.norm = norm_layer(embed_dim)
109
+
110
+ self.decoder_embed = nn.Linear(embed_dim, decoder_embed_dim)
111
+ self.mask_token = nn.Parameter(torch.zeros(1, 1, decoder_embed_dim))
112
+ self.decoder_pos_embed = nn.Parameter(
113
+ torch.zeros(1, num_patches, decoder_embed_dim), requires_grad=False
114
+ )
115
+ decoder_args = dict(block_args)
116
+ decoder_args["norm_layer"] = norm_layer
117
+ self.decoder_blocks = nn.ModuleList(
118
+ [Block(dim=decoder_embed_dim, **decoder_args) for _ in range(decoder_depth)]
119
+ )
120
+ self.decoder_norm = norm_layer(decoder_embed_dim)
121
+ patch_area = self.patch_size[0] * self.patch_size[1]
122
+ self.decoder_pred = nn.Linear(decoder_embed_dim, in_chans * patch_area)
123
+ self.initialize_weights()
124
+
125
+ def initialize_weights(self) -> None:
126
+ grid_size = self.patch_embed.grid_size
127
+ self.pos_embed.data.copy_(build_2d_sincos_position_embedding(self.embed_dim, grid_size, True))
128
+ self.decoder_pos_embed.data.copy_(
129
+ build_2d_sincos_position_embedding(self.decoder_embed_dim, grid_size, False)
130
+ )
131
+ nn.init.xavier_uniform_(self.patch_embed.proj.weight.data.flatten(1))
132
+ nn.init.normal_(self.cls_token, std=0.02)
133
+ nn.init.normal_(self.mask_token, std=0.02)
134
+ self.apply(self._init_weights)
135
+
136
+ @staticmethod
137
+ def _init_weights(module: nn.Module) -> None:
138
+ if isinstance(module, nn.Linear):
139
+ nn.init.xavier_uniform_(module.weight)
140
+ if module.bias is not None:
141
+ nn.init.zeros_(module.bias)
142
+ elif isinstance(module, nn.LayerNorm):
143
+ nn.init.zeros_(module.bias)
144
+ nn.init.ones_(module.weight)
145
+
146
+ def patchify(self, images: torch.Tensor) -> torch.Tensor:
147
+ if images.ndim != 4 or images.shape[1] != self.in_chans or tuple(images.shape[-2:]) != self.img_size:
148
+ raise ValueError(
149
+ f"patchify expects [B,{self.in_chans},{self.img_size[0]},{self.img_size[1]}], "
150
+ f"got {tuple(images.shape)}."
151
+ )
152
+ ph, pw = self.patch_size
153
+ gh, gw = self.patch_embed.grid_size
154
+ x = images.reshape(images.shape[0], self.in_chans, gh, ph, gw, pw)
155
+ x = torch.einsum("nchpwq->nhwpqc", x)
156
+ return x.reshape(images.shape[0], gh * gw, ph * pw * self.in_chans)
157
+
158
+ def unpatchify(self, patches: torch.Tensor) -> torch.Tensor:
159
+ ph, pw = self.patch_size
160
+ gh, gw = self.patch_embed.grid_size
161
+ expected_dim = ph * pw * self.in_chans
162
+ if patches.ndim != 3 or patches.shape[1:] != (gh * gw, expected_dim):
163
+ raise ValueError(f"unpatchify expects [B,{gh * gw},{expected_dim}], got {tuple(patches.shape)}.")
164
+ x = patches.reshape(patches.shape[0], gh, gw, ph, pw, self.in_chans)
165
+ x = torch.einsum("nhwpqc->nchpwq", x)
166
+ return x.reshape(patches.shape[0], self.in_chans, gh * ph, gw * pw)
167
+
168
+ @staticmethod
169
+ def random_masking(
170
+ x: torch.Tensor, mask_ratio: float
171
+ ) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
172
+ if mask_ratio not in {0.0, 0.75}:
173
+ raise ValueError("Official W-MAE AFNO grid reshaping supports mask_ratio 0.0 or 0.75 only.")
174
+ batch, length, channels = x.shape
175
+ len_keep = int(length * (1.0 - mask_ratio))
176
+ noise = torch.rand(batch, length, device=x.device)
177
+ ids_shuffle = torch.argsort(noise, dim=1)
178
+ ids_restore = torch.argsort(ids_shuffle, dim=1)
179
+ ids_keep = ids_shuffle[:, :len_keep]
180
+ x_masked = torch.gather(x, 1, ids_keep.unsqueeze(-1).expand(-1, -1, channels))
181
+ mask = torch.ones(batch, length, device=x.device)
182
+ mask[:, :len_keep] = 0
183
+ return x_masked, torch.gather(mask, 1, ids_restore), ids_restore
184
+
185
+ def forward_encoder(
186
+ self, images: torch.Tensor, mask_ratio: float
187
+ ) -> tuple[torch.Tensor, torch.Tensor, torch.Tensor]:
188
+ x = self.patch_embed(images) + self.pos_embed[:, 1:]
189
+ x, mask, ids_restore = self.random_masking(x, mask_ratio)
190
+ grid_h, grid_w = self.patch_embed.grid_size
191
+ divisor = 2 if mask_ratio == 0.75 else 1
192
+ expected_tokens = (grid_h // divisor) * (grid_w // divisor)
193
+ if x.shape[1] != expected_tokens:
194
+ raise ValueError("The selected mask ratio does not form the rectangular AFNO grid expected by W-MAE.")
195
+ x = x.reshape(x.shape[0], grid_h // divisor, grid_w // divisor, self.embed_dim)
196
+ for block in self.blocks:
197
+ x = block(x)
198
+ x = self.norm(x)
199
+ return x.flatten(1, 2), mask, ids_restore
200
+
201
+ def forward_decoder(self, latent: torch.Tensor, ids_restore: torch.Tensor) -> torch.Tensor:
202
+ x = self.decoder_embed(latent)
203
+ missing_tokens = ids_restore.shape[1] - x.shape[1]
204
+ if missing_tokens < 0:
205
+ raise ValueError("Latent token count exceeds the decoder target token count.")
206
+ mask_tokens = self.mask_token.expand(x.shape[0], missing_tokens, -1)
207
+ x = torch.cat((x, mask_tokens), dim=1)
208
+ x = torch.gather(x, 1, ids_restore.unsqueeze(-1).expand(-1, -1, x.shape[-1]))
209
+ x = x + self.decoder_pos_embed
210
+ grid_h, grid_w = self.patch_embed.grid_size
211
+ x = x.reshape(x.shape[0], grid_h, grid_w, self.decoder_embed_dim)
212
+ for block in self.decoder_blocks:
213
+ x = block(x)
214
+ x = self.decoder_pred(self.decoder_norm(x))
215
+ return x.flatten(1, 2)
216
+
217
+ def forward_loss(
218
+ self, images: torch.Tensor, prediction: torch.Tensor, mask: torch.Tensor, mask_ratio: float
219
+ ) -> torch.Tensor:
220
+ target = self.patchify(images)
221
+ if self.norm_pix_loss:
222
+ mean = target.mean(dim=-1, keepdim=True)
223
+ variance = target.var(dim=-1, keepdim=True)
224
+ target = (target - mean) / torch.sqrt(variance + 1e-6)
225
+ loss = (prediction - target).pow(2).mean(dim=-1)
226
+ if mask_ratio == 0.0:
227
+ return loss.mean()
228
+ masked_count = mask.sum()
229
+ if masked_count.item() == 0:
230
+ raise ValueError("Masked reconstruction loss requires at least one masked patch.")
231
+ return (loss * mask).sum() / masked_count
232
+
233
+ def forward(self, images: torch.Tensor, mask_ratio: float = 0.75) -> WMAEOutput:
234
+ latent, mask, ids_restore = self.forward_encoder(images, mask_ratio)
235
+ prediction = self.forward_decoder(latent, ids_restore)
236
+ loss = self.forward_loss(images, prediction, mask, mask_ratio)
237
+ return WMAEOutput(loss, prediction, mask)
238
+
239
+
240
+ def w_mae_base(
241
+ embed_dim: int = 768,
242
+ depth: int = 12,
243
+ decoder_embed_dim: int = 512,
244
+ decoder_depth: int = 6,
245
+ mlp_ratio: float = 4.0,
246
+ norm_layer: type[nn.Module] = partial(nn.LayerNorm, eps=1e-6),
247
+ **kwargs: object,
248
+ ) -> MaskedAutoencoderAFNO:
249
+ return MaskedAutoencoderAFNO(
250
+ embed_dim=embed_dim,
251
+ depth=depth,
252
+ decoder_embed_dim=decoder_embed_dim,
253
+ decoder_depth=decoder_depth,
254
+ mlp_ratio=mlp_ratio,
255
+ norm_layer=norm_layer,
256
+ **kwargs,
257
+ )
258
+
259
+
260
+ mae_vit_base_patch16 = w_mae_base
scripts/era5_adapter.py ADDED
@@ -0,0 +1,135 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ from pathlib import Path
4
+ from typing import Any, Sequence
5
+
6
+ import h5py
7
+
8
+
9
+ EXPECTED_CHANNELS = 20
10
+ SOURCE_SIZE = (721, 1440)
11
+ MODEL_SIZE = (720, 1440)
12
+
13
+
14
+ def _decode_variables(values: Sequence[Any]) -> list[str]:
15
+ return [value.decode() if isinstance(value, bytes) else str(value) for value in values]
16
+
17
+
18
+ def validate_variables(variables: Sequence[str], expected_channels: int = EXPECTED_CHANNELS) -> list[str]:
19
+ names = list(variables)
20
+ if len(names) != expected_channels:
21
+ raise ValueError(
22
+ f"W-MAE requires exactly {expected_channels} explicitly ordered variables; "
23
+ f"received {len(names)}. The official repository does not publish a complete "
24
+ "channel-name mapping, so no default mapping is assumed."
25
+ )
26
+ if len(set(names)) != len(names):
27
+ raise ValueError("W-MAE variable names must be unique and explicitly ordered.")
28
+ return names
29
+
30
+
31
+ def inspect_era5_contract(
32
+ dataset_dir: str | Path,
33
+ years: Sequence[int],
34
+ variables: Sequence[str],
35
+ expected_channels: int = EXPECTED_CHANNELS,
36
+ ) -> dict[str, Any]:
37
+ """Validate the HDF5 metadata before importing the torch-based datapipe."""
38
+ dataset_dir = Path(dataset_dir)
39
+ names = validate_variables(variables, expected_channels)
40
+ missing_years = [year for year in years if not (dataset_dir / "data" / f"{year}.h5").is_file()]
41
+ if missing_years:
42
+ raise ValueError(f"ERA5 year files are missing: {missing_years}")
43
+
44
+ first_file = dataset_dir / "data" / f"{years[0]}.h5"
45
+ with h5py.File(first_file, "r") as handle:
46
+ if "fields" not in handle:
47
+ raise ValueError(f"{first_file} does not contain a 'fields' dataset.")
48
+ fields = handle["fields"]
49
+ if fields.ndim != 4:
50
+ raise ValueError(f"fields must have shape [T,C,H,W], got {fields.shape}.")
51
+ if tuple(fields.shape[-2:]) != SOURCE_SIZE:
52
+ raise ValueError(f"W-MAE expects source grid {SOURCE_SIZE}, got {fields.shape[-2:]}.")
53
+ if "variables" not in fields.attrs or "time_step" not in fields.attrs:
54
+ raise ValueError("fields must define 'variables' and 'time_step' attributes.")
55
+ available = _decode_variables(fields.attrs["variables"])
56
+ fields_shape = tuple(fields.shape)
57
+ time_step = int(fields.attrs["time_step"])
58
+ missing_variables = [name for name in names if name not in available]
59
+ if missing_variables:
60
+ raise ValueError(f"Configured variables are absent from HDF5 metadata: {missing_variables}")
61
+ if "global_means" not in handle or "global_stds" not in handle:
62
+ stats_dir = dataset_dir / "stats"
63
+ if not (stats_dir / "global_means.npy").is_file() or not (stats_dir / "global_stds.npy").is_file():
64
+ raise ValueError("Normalization statistics are missing from HDF5 and dataset_dir/stats.")
65
+
66
+ return {
67
+ "file": str(first_file),
68
+ "fields_shape": fields_shape,
69
+ "time_step": time_step,
70
+ "selected_variables": names,
71
+ "channel_indices": [available.index(name) for name in names],
72
+ "crop": "fields[..., :720, :]",
73
+ "model_size": MODEL_SIZE,
74
+ }
75
+
76
+
77
+ def _crop_last_latitude(value: Any) -> Any:
78
+ if tuple(value.shape[-2:]) == MODEL_SIZE:
79
+ return value
80
+ if tuple(value.shape[-2:]) != SOURCE_SIZE:
81
+ raise ValueError(f"Expected trailing spatial shape {SOURCE_SIZE}, got {tuple(value.shape[-2:])}.")
82
+ return value[..., : MODEL_SIZE[0], :]
83
+
84
+
85
+ class WMAEERA5Dataset:
86
+ """Thin W-MAE adapter around OneScience ERA5Dataset.
87
+
88
+ The official W-MAE loader removes the final latitude row from a
89
+ 721x1440 ERA5 field. This wrapper preserves OneScience loading and
90
+ normalization while applying that exact spatial convention.
91
+ """
92
+
93
+ def __init__(
94
+ self,
95
+ dataset_dir: str | Path,
96
+ years: Sequence[int],
97
+ variables: Sequence[str],
98
+ task: str = "forecast",
99
+ input_steps: int = 1,
100
+ output_steps: int = 1,
101
+ normalize: bool = True,
102
+ ) -> None:
103
+ if task not in {"pretrain", "forecast"}:
104
+ raise ValueError("task must be either 'pretrain' or 'forecast'.")
105
+ names = validate_variables(variables)
106
+ inspect_era5_contract(dataset_dir, years, names)
107
+
108
+ try:
109
+ from onescience.datapipes.climate.era5 import ERA5Dataset
110
+ except (ImportError, OSError) as error:
111
+ raise RuntimeError(
112
+ "OneScience ERA5Dataset could not be imported. Verify the active "
113
+ "OneScience/PyTorch runtime before constructing WMAEERA5Dataset."
114
+ ) from error
115
+
116
+ self.task = task
117
+ self.dataset = ERA5Dataset(
118
+ dataset_dir=str(dataset_dir),
119
+ used_years=list(years),
120
+ used_variables=names,
121
+ input_steps=input_steps,
122
+ output_steps=output_steps,
123
+ normalize=normalize,
124
+ )
125
+
126
+ def __len__(self) -> int:
127
+ return len(self.dataset)
128
+
129
+ def __getitem__(self, index: int) -> tuple[Any, Any, Any, int, list[str]]:
130
+ invar, outvar, cos_zenith, step_idx, time_index = self.dataset[index]
131
+ invar = _crop_last_latitude(invar)
132
+ outvar = _crop_last_latitude(outvar)
133
+ cos_zenith = _crop_last_latitude(cos_zenith)
134
+ target = invar if self.task == "pretrain" else outvar
135
+ return invar, target, cos_zenith, step_idx, time_index
scripts/fake_data.py ADDED
@@ -0,0 +1,84 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import argparse
4
+ from pathlib import Path
5
+ import sys
6
+
7
+ import h5py
8
+ import numpy as np
9
+ import yaml
10
+
11
+
12
+ PROJECT_ROOT = Path(__file__).resolve().parents[1]
13
+ sys.path.insert(0, str(PROJECT_ROOT))
14
+ sys.path.insert(0, str(PROJECT_ROOT / "scripts"))
15
+
16
+ from era5_adapter import inspect_era5_contract
17
+
18
+
19
+ def load_config(path: Path) -> dict:
20
+ with path.open("r", encoding="utf-8") as handle:
21
+ return yaml.safe_load(handle)
22
+
23
+
24
+ def generate_year_file(
25
+ output_path: Path,
26
+ variables: list[str],
27
+ time_steps: int,
28
+ height: int,
29
+ width: int,
30
+ time_step_hours: int,
31
+ ) -> None:
32
+ output_path.parent.mkdir(parents=True, exist_ok=True)
33
+ channels = len(variables)
34
+ means = np.zeros((1, channels, 1, 1), dtype=np.float32)
35
+ stds = np.ones((1, channels, 1, 1), dtype=np.float32)
36
+
37
+ with h5py.File(output_path, "w") as handle:
38
+ fields = handle.create_dataset(
39
+ "fields",
40
+ shape=(time_steps, channels, height, width),
41
+ dtype="float32",
42
+ chunks=(1, 1, height, width),
43
+ fillvalue=0.0,
44
+ )
45
+ fields.attrs["variables"] = variables
46
+ fields.attrs["time_step"] = time_step_hours
47
+ handle.create_dataset("global_means", data=means)
48
+ handle.create_dataset("global_stds", data=stds)
49
+
50
+
51
+ def main() -> None:
52
+ parser = argparse.ArgumentParser(description="Generate lightweight ERA5-style HDF5 files for W-MAE.")
53
+ parser.add_argument("--config", type=Path, default=PROJECT_ROOT / "conf" / "config.yaml")
54
+ parser.add_argument("--output-dir", type=Path, default=None)
55
+ args = parser.parse_args()
56
+
57
+ config = load_config(args.config)
58
+ data_config = config["data"]
59
+ fake_config = config["fake_data"]
60
+ output_dir = args.output_dir or PROJECT_ROOT / data_config["dataset_dir"]
61
+ variables = list(fake_config["variables"])
62
+ years = list(fake_config["years"])
63
+ source_height, source_width = data_config["source_size"]
64
+
65
+ for year in years:
66
+ path = output_dir / "data" / f"{year}.h5"
67
+ generate_year_file(
68
+ output_path=path,
69
+ variables=variables,
70
+ time_steps=int(fake_config["time_steps"]),
71
+ height=int(source_height),
72
+ width=int(source_width),
73
+ time_step_hours=int(data_config["time_step_hours"]),
74
+ )
75
+ print(f"created {path} ({path.stat().st_size / 1024:.1f} KiB physical size)")
76
+
77
+ report = inspect_era5_contract(output_dir, years, variables)
78
+ print(f"validated fields shape: {report['fields_shape']}")
79
+ print(f"W-MAE spatial adapter: {report['crop']} -> {report['model_size']}")
80
+ print("Synthetic channel names are test-only and must not be used for real training.")
81
+
82
+
83
+ if __name__ == "__main__":
84
+ main()
scripts/inference.py ADDED
@@ -0,0 +1,177 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import argparse
4
+ from pathlib import Path
5
+ import sys
6
+ from typing import Any
7
+
8
+ import numpy as np
9
+ import torch
10
+ from torch.utils.data import DataLoader
11
+
12
+
13
+ PROJECT_ROOT = Path(__file__).resolve().parents[1]
14
+ if str(PROJECT_ROOT) not in sys.path:
15
+ sys.path.insert(0, str(PROJECT_ROOT))
16
+ SCRIPTS_DIR = PROJECT_ROOT / "scripts"
17
+ if str(SCRIPTS_DIR) not in sys.path:
18
+ sys.path.insert(0, str(SCRIPTS_DIR))
19
+
20
+ from model.checkpoint import load_checkpoint
21
+ from era5_adapter import WMAEERA5Dataset
22
+ from train import build_model, device_summary, load_config, resolve_device, resolve_project_path, validate_config
23
+
24
+
25
+ def build_inference_dataset(config: dict[str, Any], config_path: Path) -> WMAEERA5Dataset:
26
+ data = config["data"]
27
+ if int(data["input_steps"]) != 1 or int(data["output_steps"]) != 1:
28
+ raise ValueError("W-MAE inference currently requires input_steps=1 and output_steps=1.")
29
+ if not data["variables"]:
30
+ raise ValueError(
31
+ "data.variables is empty. Supply the verified, explicitly ordered 20-channel list before inference."
32
+ )
33
+ return WMAEERA5Dataset(
34
+ dataset_dir=resolve_project_path(data["dataset_dir"], config_path),
35
+ years=data["test_years"],
36
+ variables=data["variables"],
37
+ task="pretrain",
38
+ input_steps=1,
39
+ output_steps=1,
40
+ normalize=bool(data["normalize"]),
41
+ )
42
+
43
+
44
+ def sample_time_index(time_index: Any, batch_index: int, batch_size: int) -> Any:
45
+ """Extract one sample from DataLoader's sequence-major default collation."""
46
+ if isinstance(time_index, (list, tuple)):
47
+ if time_index and isinstance(time_index[0], (list, tuple)):
48
+ return time_index[0][batch_index]
49
+ if len(time_index) == batch_size:
50
+ return time_index[batch_index]
51
+ return time_index
52
+
53
+
54
+ def run_inference(
55
+ model: torch.nn.Module,
56
+ loader: DataLoader,
57
+ output_dir: Path,
58
+ mask_ratio: float,
59
+ device: torch.device,
60
+ max_samples: int | None = None,
61
+ log_interval: int = 1,
62
+ ) -> int:
63
+ output_dir.mkdir(parents=True, exist_ok=True)
64
+ model.eval()
65
+ written = 0
66
+ with torch.no_grad():
67
+ total_batches = len(loader)
68
+ for batch_index, (inputs, _, _, step_idx, time_index) in enumerate(loader, start=1):
69
+ if inputs.ndim != 4:
70
+ raise ValueError(f"W-MAE inference expects a 4D batch, got {tuple(inputs.shape)}.")
71
+ inputs = inputs.to(device)
72
+ result = model(inputs, mask_ratio=mask_ratio)
73
+ reconstruction = model.unpatchify(result.prediction)
74
+ error = reconstruction - inputs
75
+ mask = result.mask.reshape(result.mask.shape[0], *model.patch_embed.grid_size)
76
+ for sample_index in range(inputs.shape[0]):
77
+ if max_samples is not None and written >= max_samples:
78
+ print(
79
+ f"inference batch={batch_index}/{total_batches} samples_written={written}",
80
+ flush=True,
81
+ )
82
+ return written
83
+ sample_time = sample_time_index(time_index, sample_index, inputs.shape[0])
84
+ if isinstance(sample_time, (list, tuple)):
85
+ sample_time = sample_time[0]
86
+ np.savez_compressed(
87
+ output_dir / f"sample_{written:06d}.npz",
88
+ input=inputs[sample_index].cpu().numpy(),
89
+ reconstruction=reconstruction[sample_index].cpu().numpy(),
90
+ error=error[sample_index].cpu().numpy(),
91
+ mask=mask[sample_index].cpu().numpy(),
92
+ step_idx=np.asarray(step_idx[sample_index].item()),
93
+ time_index=np.asarray(sample_time),
94
+ )
95
+ written += 1
96
+ if batch_index == 1 or batch_index % log_interval == 0 or batch_index == total_batches:
97
+ print(
98
+ f"inference batch={batch_index}/{total_batches} samples_written={written}",
99
+ flush=True,
100
+ )
101
+ return written
102
+
103
+
104
+ def main() -> None:
105
+ parser = argparse.ArgumentParser(description="Run W-MAE reconstruction inference on ERA5 samples.")
106
+ parser.add_argument("--config", type=Path, default=PROJECT_ROOT / "conf" / "config.yaml")
107
+ parser.add_argument("--checkpoint", type=Path, default="./data/checkpoint/model_bak.pth")
108
+ parser.add_argument("--checkpoint-source-root", type=Path, default=None)
109
+ parser.add_argument("--non-strict-checkpoint", action="store_true")
110
+ parser.add_argument("--output-dir", type=Path, default=None)
111
+ parser.add_argument(
112
+ "--device",
113
+ default="auto",
114
+ help="Inference device: auto (default), cuda (AMD DCU through HIP), or cpu.",
115
+ )
116
+ parser.add_argument("--mask-ratio", type=float, default=None)
117
+ parser.add_argument("--max-samples", type=int, default=None)
118
+ parser.add_argument("--log-interval", type=int, default=1)
119
+ args = parser.parse_args()
120
+ if args.max_samples is not None and args.max_samples <= 0:
121
+ raise ValueError("--max-samples must be positive when provided.")
122
+ if args.log_interval <= 0:
123
+ raise ValueError("--log-interval must be positive.")
124
+
125
+ config_path = args.config.resolve()
126
+ print(
127
+ f"inference started: config={config_path} checkpoint={args.checkpoint} device={args.device}",
128
+ flush=True,
129
+ )
130
+ config = load_config(config_path)
131
+ validate_config(config)
132
+ if args.checkpoint is None:
133
+ raise ValueError("Inference requires an explicit --checkpoint path.")
134
+
135
+ device = resolve_device(args.device)
136
+ print(f"runtime: {device_summary(device)}", flush=True)
137
+ print("building model", flush=True)
138
+ model = build_model(config).to(device)
139
+ print(f"model ready: parameters_device={next(model.parameters()).device}", flush=True)
140
+ source_root = (
141
+ resolve_project_path(args.checkpoint_source_root, config_path)
142
+ if args.checkpoint_source_root
143
+ else None
144
+ )
145
+ report = load_checkpoint(
146
+ model,
147
+ resolve_project_path(args.checkpoint, config_path),
148
+ strict=not args.non_strict_checkpoint,
149
+ map_location=device,
150
+ source_root=source_root,
151
+ )
152
+ print(f"loaded checkpoint: {report}", flush=True)
153
+
154
+ data_config = config["data"]
155
+ print("building inference dataset", flush=True)
156
+ dataset = build_inference_dataset(config, config_path)
157
+ loader = DataLoader(
158
+ dataset,
159
+ batch_size=int(data_config["batch_size"]),
160
+ shuffle=False,
161
+ num_workers=int(data_config["num_workers"]),
162
+ )
163
+ print(f"dataset ready: samples={len(dataset)} batches={len(loader)}", flush=True)
164
+ output_dir = resolve_project_path(args.output_dir, config_path) if args.output_dir else (
165
+ resolve_project_path(config["project"]["output_dir"], config_path) / "inference"
166
+ )
167
+ mask_ratio = float(config["model"]["mask_ratio"] if args.mask_ratio is None else args.mask_ratio)
168
+ if mask_ratio not in {0.0, 0.75}:
169
+ raise ValueError("W-MAE supports mask_ratio 0.0 or 0.75 only.")
170
+ print(f"inference running: mask_ratio={mask_ratio} output_dir={output_dir}", flush=True)
171
+ count = run_inference(model, loader, output_dir, mask_ratio, device, args.max_samples, args.log_interval)
172
+ print(f"wrote {count} reconstruction samples to {output_dir}", flush=True)
173
+ print("inference completed successfully", flush=True)
174
+
175
+
176
+ if __name__ == "__main__":
177
+ main()
scripts/result.py ADDED
@@ -0,0 +1,180 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import argparse
4
+ import json
5
+ from pathlib import Path
6
+ from typing import Any
7
+
8
+ import numpy as np
9
+
10
+
11
+ PROJECT_ROOT = Path(__file__).resolve().parents[1]
12
+ DEFAULT_RESULTS_DIR = PROJECT_ROOT / "outputs" / "inference"
13
+
14
+
15
+ def _load_array(path: Path, name: str) -> np.ndarray:
16
+ value = np.load(path)[name]
17
+ if not np.all(np.isfinite(value)):
18
+ raise ValueError(f"{path.name}:{name} contains NaN or Inf values.")
19
+ return value
20
+
21
+
22
+ def validate_sample(path: Path, expected_channels: int = 20) -> dict[str, Any]:
23
+ with np.load(path) as result:
24
+ required = {"input", "reconstruction", "error", "mask", "step_idx", "time_index"}
25
+ missing = required - set(result.files)
26
+ if missing:
27
+ raise ValueError(f"{path} is missing fields: {sorted(missing)}")
28
+ inputs = _load_array(path, "input")
29
+ reconstruction = _load_array(path, "reconstruction")
30
+ error = _load_array(path, "error")
31
+ mask = _load_array(path, "mask")
32
+
33
+ expected_shape = (expected_channels, 720, 1440)
34
+ if inputs.shape != expected_shape:
35
+ raise ValueError(f"input shape must be {expected_shape}, got {inputs.shape}.")
36
+ if reconstruction.shape != inputs.shape or error.shape != inputs.shape:
37
+ raise ValueError("reconstruction and error must have the same shape as input.")
38
+ if mask.shape != (90, 180):
39
+ raise ValueError(f"mask must be a patch grid with shape (90, 180), got {mask.shape}.")
40
+ if not np.array_equal(mask, mask.astype(bool)):
41
+ raise ValueError("mask must contain only binary values.")
42
+ if not np.allclose(error, reconstruction - inputs):
43
+ raise ValueError("error does not equal reconstruction - input.")
44
+
45
+ spatial_mask = np.repeat(np.repeat(mask.astype(bool), 8, axis=0), 8, axis=1)
46
+ visible_error = np.abs(error[:, ~spatial_mask])
47
+ masked_error = np.abs(error[:, spatial_mask])
48
+ mask_fraction = float(mask.mean())
49
+ return {
50
+ "file": str(path),
51
+ "shape": list(inputs.shape),
52
+ "mask_fraction": mask_fraction,
53
+ "masked_pixels": int(spatial_mask.sum()),
54
+ "visible_mae": np.mean(visible_error, axis=1).tolist() if visible_error.size else [0.0] * expected_channels,
55
+ "masked_mae": np.mean(masked_error, axis=1).tolist() if masked_error.size else [0.0] * expected_channels,
56
+ "mae": np.mean(np.abs(error), axis=(1, 2)).tolist(),
57
+ "rmse": np.sqrt(np.mean(error**2, axis=(1, 2))).tolist(),
58
+ }
59
+
60
+
61
+ def save_diagnostic_png(path: Path, output_path: Path, channel: int) -> None:
62
+ import matplotlib
63
+
64
+ matplotlib.use("Agg")
65
+ import matplotlib.pyplot as plt
66
+
67
+ with np.load(path) as result:
68
+ input_field = result["input"][channel]
69
+ reconstruction = result["reconstruction"][channel]
70
+ error = result["error"][channel]
71
+ mask = result["mask"].astype(bool)
72
+ time_index = result["time_index"].item()
73
+
74
+ input_low, input_high = np.nanpercentile(input_field, [1, 99])
75
+ reconstruction_low, reconstruction_high = np.nanpercentile(reconstruction, [1, 99])
76
+ field_low = min(input_low, reconstruction_low)
77
+ field_high = max(input_high, reconstruction_high)
78
+ if not np.isfinite(field_low) or not np.isfinite(field_high) or field_high <= field_low:
79
+ field_low, field_high = float(np.nanmin(input_field)), float(np.nanmax(input_field))
80
+ error_limit = float(np.nanpercentile(np.abs(error), 99))
81
+ if not np.isfinite(error_limit) or error_limit <= 0:
82
+ error_limit = max(float(np.nanmax(np.abs(error))), 1.0)
83
+
84
+ latitude = np.linspace(90, -90, input_field.shape[0])
85
+ longitude = np.linspace(0, 360, input_field.shape[1], endpoint=False)
86
+ extent = [longitude[0], longitude[-1], latitude[-1], latitude[0]]
87
+ masked_pixels = np.repeat(np.repeat(mask, 8, axis=0), 8, axis=1)
88
+ masked_overlay = np.ma.masked_where(~masked_pixels, masked_pixels)
89
+
90
+ figure, axes = plt.subplots(2, 3, figsize=(18, 8), constrained_layout=True)
91
+ field_kwargs = {"cmap": "viridis", "vmin": field_low, "vmax": field_high, "extent": extent, "aspect": "auto"}
92
+ error_kwargs = {"cmap": "RdBu_r", "vmin": -error_limit, "vmax": error_limit, "extent": extent, "aspect": "auto"}
93
+ panels = [
94
+ (axes[0, 0], input_field, "Input (normalized)", field_kwargs),
95
+ (axes[0, 1], reconstruction, "Reconstruction (normalized)", field_kwargs),
96
+ (axes[0, 2], error, "Signed error", error_kwargs),
97
+ (axes[1, 0], np.abs(error), "Absolute error", {**error_kwargs, "cmap": "magma", "vmin": 0}),
98
+ ]
99
+ for axis, field, title, kwargs in panels:
100
+ image = axis.imshow(field, **kwargs)
101
+ axis.set_title(title)
102
+ axis.set_xlabel("Longitude (degrees east)")
103
+ axis.set_ylabel("Latitude (degrees)")
104
+ axis.contour(masked_pixels, levels=[0.5], colors="white", linewidths=0.35, extent=extent)
105
+ figure.colorbar(image, ax=axis, shrink=0.8)
106
+
107
+ mask_image = axes[1, 1].imshow(mask.astype(float), cmap="gray_r", vmin=0, vmax=1, extent=extent, aspect="auto")
108
+ axes[1, 1].set_title("Patch mask (white = masked)")
109
+ axes[1, 1].set_xlabel("Longitude (degrees east)")
110
+ axes[1, 1].set_ylabel("Latitude (degrees)")
111
+ figure.colorbar(mask_image, ax=axes[1, 1], ticks=[0, 1], shrink=0.8)
112
+
113
+ axes[1, 2].hist(error.ravel(), bins=80, color="#315f8c", alpha=0.85)
114
+ axes[1, 2].axvline(0, color="black", linewidth=0.8)
115
+ axes[1, 2].set_title(f"Error distribution\nMAE={np.mean(np.abs(error)):.4f}, RMSE={np.sqrt(np.mean(error**2)):.4f}")
116
+ axes[1, 2].set_xlabel("Signed error (normalized)")
117
+ axes[1, 2].set_ylabel("Pixel count")
118
+
119
+ figure.suptitle(f"W-MAE reconstruction diagnostic | channel={channel} | time={time_index}")
120
+ output_path.parent.mkdir(parents=True, exist_ok=True)
121
+ figure.savefig(output_path, dpi=150)
122
+ plt.close(figure)
123
+
124
+
125
+ def main() -> None:
126
+ parser = argparse.ArgumentParser(description="Validate W-MAE inference outputs and make diagnostics.")
127
+ parser.add_argument(
128
+ "results_dir",
129
+ type=Path,
130
+ nargs="?",
131
+ default=DEFAULT_RESULTS_DIR,
132
+ help="Inference output directory (default: outputs/inference).",
133
+ )
134
+ parser.add_argument("--expected-channels", type=int, default=20)
135
+ parser.add_argument("--mask-ratio", type=float, default=None)
136
+ parser.add_argument(
137
+ "--channel",
138
+ type=int,
139
+ default=0,
140
+ help="Channel for the diagnostic image (default: 0). Use --channel -1 to disable it.",
141
+ )
142
+ parser.add_argument("--diagnostic-dir", type=Path, default=None)
143
+ parser.add_argument(
144
+ "--visualization-samples",
145
+ type=int,
146
+ default=1,
147
+ help="Number of samples to visualize; validation still covers all samples (default: 1).",
148
+ )
149
+ args = parser.parse_args()
150
+ if args.visualization_samples <= 0:
151
+ raise ValueError("--visualization-samples must be positive.")
152
+ files = sorted(args.results_dir.glob("sample_*.npz"))
153
+ if not files:
154
+ raise FileNotFoundError(f"No sample_*.npz files found in {args.results_dir}.")
155
+ reports = [validate_sample(path, args.expected_channels) for path in files]
156
+ if args.mask_ratio is not None:
157
+ for report in reports:
158
+ if not np.isclose(report["mask_fraction"], args.mask_ratio):
159
+ raise ValueError(
160
+ f"{report['file']} mask fraction {report['mask_fraction']} "
161
+ f"does not match expected {args.mask_ratio}."
162
+ )
163
+ if args.channel is not None and args.channel >= 0:
164
+ if not 0 <= args.channel < args.expected_channels:
165
+ raise ValueError("--channel is outside the configured channel range.")
166
+ diagnostic_dir = args.diagnostic_dir or args.results_dir / "diagnostics"
167
+ for path in files[:args.visualization_samples]:
168
+ save_diagnostic_png(path, diagnostic_dir / f"{path.stem}_channel_{args.channel:02d}.png", args.channel)
169
+ summary = {
170
+ "samples": len(reports),
171
+ "mean_mask_fraction": float(np.mean([report["mask_fraction"] for report in reports])),
172
+ "mean_mae": np.mean([report["mae"] for report in reports], axis=0).tolist(),
173
+ "mean_rmse": np.mean([report["rmse"] for report in reports], axis=0).tolist(),
174
+ "reports": reports,
175
+ }
176
+ print(json.dumps(summary, indent=2), flush=True)
177
+
178
+
179
+ if __name__ == "__main__":
180
+ main()
scripts/train.py ADDED
@@ -0,0 +1,384 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ from __future__ import annotations
2
+
3
+ import argparse
4
+ import os
5
+ from pathlib import Path
6
+ import sys
7
+ import traceback
8
+ from typing import Any
9
+
10
+ import torch
11
+ import torch.distributed as dist
12
+ from torch.nn.parallel import DistributedDataParallel
13
+ from torch.utils.data import DataLoader
14
+ from torch.utils.data.distributed import DistributedSampler
15
+ import yaml
16
+
17
+
18
+ PROJECT_ROOT = Path(__file__).resolve().parents[1]
19
+ if str(PROJECT_ROOT) not in sys.path:
20
+ sys.path.insert(0, str(PROJECT_ROOT))
21
+ SCRIPTS_DIR = PROJECT_ROOT / "scripts"
22
+ if str(SCRIPTS_DIR) not in sys.path:
23
+ sys.path.insert(0, str(SCRIPTS_DIR))
24
+
25
+ from model.checkpoint import load_checkpoint
26
+ from model.w_mae import w_mae_base
27
+ from era5_adapter import WMAEERA5Dataset
28
+
29
+
30
+ def load_config(path: Path) -> dict[str, Any]:
31
+ with path.open("r", encoding="utf-8") as handle:
32
+ config = yaml.safe_load(handle)
33
+ if not isinstance(config, dict):
34
+ raise ValueError(f"Configuration must contain a mapping: {path}")
35
+ return config
36
+
37
+
38
+ def resolve_project_path(value: str | Path, config_path: Path) -> Path:
39
+ path = Path(value)
40
+ return path if path.is_absolute() else (config_path.parent.parent / path).resolve()
41
+
42
+
43
+ def validate_config(config: dict[str, Any]) -> None:
44
+ data = config["data"]
45
+ model = config["model"]
46
+ if int(data["expected_channels"]) != int(model["input_channels"]):
47
+ raise ValueError("data.expected_channels must equal model.input_channels.")
48
+ if tuple(data["model_size"]) != tuple(model["image_size"]):
49
+ raise ValueError("data.model_size must equal model.image_size.")
50
+ if tuple(model["patch_size"])[0] <= 0 or tuple(model["patch_size"])[1] <= 0:
51
+ raise ValueError("model.patch_size values must be positive.")
52
+ if float(model["mask_ratio"]) not in {0.0, 0.75}:
53
+ raise ValueError("W-MAE supports mask_ratio 0.0 or 0.75 only.")
54
+
55
+
56
+ def build_model(config: dict[str, Any]) -> torch.nn.Module:
57
+ model_config = config["model"]
58
+ return w_mae_base(
59
+ img_size=tuple(model_config["image_size"]),
60
+ patch_size=tuple(model_config["patch_size"]),
61
+ in_chans=int(model_config["input_channels"]),
62
+ embed_dim=int(model_config["embed_dim"]),
63
+ depth=int(model_config["encoder_depth"]),
64
+ decoder_embed_dim=int(model_config["decoder_embed_dim"]),
65
+ decoder_depth=int(model_config["decoder_depth"]),
66
+ mlp_ratio=float(model_config.get("mlp_ratio", 4.0)),
67
+ norm_pix_loss=bool(model_config["norm_pixel_loss"]),
68
+ )
69
+
70
+
71
+ def build_dataset(config: dict[str, Any], config_path: Path, split: str) -> WMAEERA5Dataset:
72
+ data = config["data"]
73
+ years = data[f"{split}_years"]
74
+ if not data["variables"]:
75
+ raise ValueError(
76
+ "data.variables is empty. Supply the verified, explicitly ordered 20-channel list before training."
77
+ )
78
+ return WMAEERA5Dataset(
79
+ dataset_dir=resolve_project_path(data["dataset_dir"], config_path),
80
+ years=years,
81
+ variables=data["variables"],
82
+ task="pretrain",
83
+ input_steps=int(data["input_steps"]),
84
+ output_steps=int(data["output_steps"]),
85
+ normalize=bool(data["normalize"]),
86
+ )
87
+
88
+
89
+ def validate_architecture_reference(
90
+ model: torch.nn.Module,
91
+ reference_path: Path,
92
+ config_path: Path,
93
+ source_root: Path | None,
94
+ ) -> None:
95
+ reference_model = build_model(load_config(config_path))
96
+ load_checkpoint(reference_model, reference_path, source_root=source_root)
97
+ model_state = model.state_dict()
98
+ reference_state = reference_model.state_dict()
99
+ if set(model_state) != set(reference_state):
100
+ raise RuntimeError("Model parameter names do not match the official architecture reference.")
101
+ mismatches = {
102
+ key: (tuple(model_state[key].shape), tuple(reference_state[key].shape))
103
+ for key in model_state
104
+ if tuple(model_state[key].shape) != tuple(reference_state[key].shape)
105
+ }
106
+ if mismatches:
107
+ raise RuntimeError(f"Model parameter shapes do not match the official architecture: {mismatches}")
108
+
109
+
110
+ def save_checkpoint(
111
+ path: Path,
112
+ model: torch.nn.Module,
113
+ optimizer: torch.optim.Optimizer,
114
+ epoch: int,
115
+ loss: float,
116
+ config: dict[str, Any],
117
+ ) -> None:
118
+ path.parent.mkdir(parents=True, exist_ok=True)
119
+ torch.save(
120
+ {
121
+ "model": model.state_dict(),
122
+ "optimizer": optimizer.state_dict(),
123
+ "epoch": epoch,
124
+ "loss": loss,
125
+ "config": config,
126
+ },
127
+ path,
128
+ )
129
+
130
+
131
+ def resolve_device(requested: str) -> torch.device:
132
+ if requested == "auto":
133
+ requested = "cuda" if torch.cuda.is_available() else "cpu"
134
+ device = torch.device(requested)
135
+ if device.type == "cuda" and not torch.cuda.is_available():
136
+ raise RuntimeError("The requested cuda device is unavailable. Use --device cpu or check the DCU runtime.")
137
+ return device
138
+
139
+
140
+ def device_summary(device: torch.device) -> str:
141
+ if device.type != "cuda":
142
+ return f"device={device}"
143
+ backend = f"HIP {torch.version.hip}" if torch.version.hip else f"CUDA {torch.version.cuda}"
144
+ return f"device={device} name={torch.cuda.get_device_name(device)} backend={backend}"
145
+
146
+
147
+ def memory_summary(device: torch.device) -> str:
148
+ if device.type != "cuda":
149
+ return "memory=unavailable(cpu)"
150
+ allocated = torch.cuda.memory_allocated(device) / (1024**3)
151
+ reserved = torch.cuda.memory_reserved(device) / (1024**3)
152
+ peak = torch.cuda.max_memory_allocated(device) / (1024**3)
153
+ return f"memory_allocated={allocated:.2f}GB memory_reserved={reserved:.2f}GB peak_allocated={peak:.2f}GB"
154
+
155
+
156
+ def setup_distributed(requested_device: str) -> tuple[torch.device, int, int, int]:
157
+ world_size = int(os.environ.get("WORLD_SIZE", "1"))
158
+ rank = int(os.environ.get("RANK", "0"))
159
+ local_rank = int(os.environ.get("LOCAL_RANK", "0"))
160
+ if world_size == 1:
161
+ return resolve_device(requested_device), rank, world_size, local_rank
162
+
163
+ if requested_device == "cpu":
164
+ device = torch.device("cpu")
165
+ backend = "gloo"
166
+ else:
167
+ if not torch.cuda.is_available():
168
+ raise RuntimeError("torchrun requires one accessible DCU/GPU per local rank.")
169
+ if local_rank >= torch.cuda.device_count():
170
+ raise RuntimeError(
171
+ f"LOCAL_RANK={local_rank} exceeds the {torch.cuda.device_count()} visible accelerator devices."
172
+ )
173
+ torch.cuda.set_device(local_rank)
174
+ device = torch.device("cuda", local_rank)
175
+ backend = "nccl"
176
+ dist.init_process_group(backend=backend, init_method="env://")
177
+ return device, rank, world_size, local_rank
178
+
179
+
180
+ def rank_print(message: str, rank: int) -> None:
181
+ if rank == 0:
182
+ print(message, flush=True)
183
+
184
+
185
+ def all_ranks_finite(value: torch.Tensor, world_size: int) -> bool:
186
+ finite = torch.isfinite(value).to(dtype=torch.int32)
187
+ if world_size > 1:
188
+ dist.all_reduce(finite, op=dist.ReduceOp.MIN)
189
+ return bool(finite.item())
190
+
191
+
192
+ def unwrap_model(model: torch.nn.Module) -> torch.nn.Module:
193
+ return model.module if isinstance(model, DistributedDataParallel) else model
194
+
195
+
196
+ def train_one_epoch(
197
+ model: torch.nn.Module,
198
+ loader: DataLoader,
199
+ optimizer: torch.optim.Optimizer,
200
+ mask_ratio: float,
201
+ device: torch.device,
202
+ epoch: int,
203
+ log_interval: int,
204
+ rank: int,
205
+ world_size: int,
206
+ ) -> float:
207
+ model.train()
208
+ total_loss = 0.0
209
+ batches = 0
210
+ total_batches = len(loader)
211
+ if device.type == "cuda":
212
+ torch.cuda.reset_peak_memory_stats(device)
213
+ for batch_index, (inputs, _, _, _, _) in enumerate(loader, start=1):
214
+ inputs = inputs.to(device)
215
+ output = model(inputs, mask_ratio=mask_ratio)
216
+ if not all_ranks_finite(output.loss.detach(), world_size):
217
+ raise FloatingPointError("Training loss is NaN or Inf on at least one rank.")
218
+ optimizer.zero_grad(set_to_none=True)
219
+ output.loss.backward()
220
+ optimizer.step()
221
+ total_loss += float(output.loss.detach().cpu())
222
+ batches += 1
223
+ if batch_index == 1 or batch_index % log_interval == 0 or batch_index == total_batches:
224
+ display_loss = output.loss.detach()
225
+ if world_size > 1:
226
+ dist.all_reduce(display_loss, op=dist.ReduceOp.SUM)
227
+ display_loss /= world_size
228
+ rank_print(
229
+ f"epoch={epoch} batch={batch_index}/{total_batches} loss={float(display_loss):.6f} "
230
+ f"{memory_summary(device)}",
231
+ rank,
232
+ )
233
+ if batches == 0:
234
+ raise ValueError("Training DataLoader produced no batches.")
235
+ totals = torch.tensor([total_loss, batches], dtype=torch.float64, device=device)
236
+ if world_size > 1:
237
+ dist.all_reduce(totals, op=dist.ReduceOp.SUM)
238
+ return float((totals[0] / totals[1]).cpu())
239
+
240
+
241
+ def main() -> None:
242
+ parser = argparse.ArgumentParser(description="Train W-MAE using the configured ERA5 adapter.")
243
+ parser.add_argument("--config", type=Path, default=PROJECT_ROOT / "conf" / "config.yaml")
244
+ parser.add_argument("--epochs", type=int, default=1)
245
+ parser.add_argument("--learning-rate", type=float, default=1e-4)
246
+ parser.add_argument(
247
+ "--device",
248
+ default="auto",
249
+ help="Training device: auto (default), cuda (AMD DCU through HIP), or cpu.",
250
+ )
251
+ parser.add_argument("--checkpoint", type=Path, default=None)
252
+ parser.add_argument("--checkpoint-source-root", type=Path, default=None)
253
+ parser.add_argument("--non-strict-checkpoint", action="store_true")
254
+ parser.add_argument(
255
+ "--output-checkpoint",
256
+ type=Path,
257
+ default=PROJECT_ROOT / "data" / "checkpoint" / "model_bak.pth",
258
+ )
259
+ parser.add_argument("--architecture-reference-checkpoint", type=Path, default=None)
260
+ parser.add_argument("--log-interval", type=int, default=1)
261
+ args = parser.parse_args()
262
+ if args.epochs <= 0:
263
+ raise ValueError("--epochs must be positive.")
264
+ if args.log_interval <= 0:
265
+ raise ValueError("--log-interval must be positive.")
266
+
267
+ config_path = args.config.resolve()
268
+ device, rank, world_size, local_rank = setup_distributed(args.device)
269
+ rank_print(
270
+ f"training started: config={config_path} epochs={args.epochs} device={args.device}", rank
271
+ )
272
+ rank_print(
273
+ f"distributed: enabled={world_size > 1} world_size={world_size} "
274
+ f"rank={rank} local_rank={local_rank}",
275
+ rank,
276
+ )
277
+ config = load_config(config_path)
278
+ validate_config(config)
279
+ data_config = config["data"]
280
+ print(f"rank={rank} local_rank={local_rank} runtime: {device_summary(device)}", flush=True)
281
+ rank_print("building model", rank)
282
+ model = build_model(config).to(device)
283
+ print(
284
+ f"rank={rank} local_rank={local_rank} model_parameters_device="
285
+ f"{next(model.parameters()).device}",
286
+ flush=True,
287
+ )
288
+
289
+ if args.checkpoint is not None:
290
+ checkpoint_path = resolve_project_path(args.checkpoint, config_path)
291
+ source_root = (
292
+ resolve_project_path(args.checkpoint_source_root, config_path)
293
+ if args.checkpoint_source_root
294
+ else None
295
+ )
296
+ report = load_checkpoint(
297
+ model,
298
+ checkpoint_path,
299
+ strict=not args.non_strict_checkpoint,
300
+ map_location=device,
301
+ source_root=source_root,
302
+ )
303
+ rank_print(f"loaded checkpoint: {report}", rank)
304
+ else:
305
+ rank_print("training from randomly initialized weights", rank)
306
+
307
+ if world_size > 1:
308
+ if device.type == "cuda":
309
+ model = DistributedDataParallel(model, device_ids=[local_rank], output_device=local_rank)
310
+ else:
311
+ model = DistributedDataParallel(model)
312
+ rank_print("building dataset", rank)
313
+ dataset = build_dataset(config, config_path, "train")
314
+ sampler = DistributedSampler(dataset, shuffle=True) if world_size > 1 else None
315
+ loader = DataLoader(
316
+ dataset,
317
+ batch_size=int(data_config["batch_size"]),
318
+ shuffle=sampler is None,
319
+ sampler=sampler,
320
+ num_workers=int(data_config["num_workers"]),
321
+ )
322
+ rank_print(
323
+ f"dataset ready: samples={len(dataset)} batches_per_rank={len(loader)} "
324
+ f"batch_size_per_rank={data_config['batch_size']} "
325
+ f"global_batch_size={int(data_config['batch_size']) * world_size}",
326
+ rank,
327
+ )
328
+ optimizer = torch.optim.AdamW(model.parameters(), lr=args.learning_rate)
329
+ for epoch in range(args.epochs):
330
+ if sampler is not None:
331
+ sampler.set_epoch(epoch)
332
+ rank_print(f"epoch={epoch + 1}/{args.epochs} started", rank)
333
+ loss = train_one_epoch(
334
+ model,
335
+ loader,
336
+ optimizer,
337
+ float(config["model"]["mask_ratio"]),
338
+ device,
339
+ epoch + 1,
340
+ args.log_interval,
341
+ rank,
342
+ world_size,
343
+ )
344
+ rank_print(f"epoch={epoch + 1}/{args.epochs} finished loss={loss:.6f}", rank)
345
+
346
+ if world_size > 1:
347
+ dist.barrier()
348
+
349
+ if rank != 0:
350
+ return
351
+
352
+ checkpoint_model = unwrap_model(model)
353
+
354
+ reference_path = (
355
+ resolve_project_path(args.architecture_reference_checkpoint, config_path)
356
+ if args.architecture_reference_checkpoint
357
+ else None
358
+ )
359
+ reference_source_root = (
360
+ resolve_project_path(args.checkpoint_source_root, config_path)
361
+ if args.checkpoint_source_root
362
+ else None
363
+ )
364
+ if reference_path is not None:
365
+ print(f"checking architecture against {reference_path}", flush=True)
366
+ validate_architecture_reference(checkpoint_model, reference_path, config_path, reference_source_root)
367
+ print("architecture reference check passed", flush=True)
368
+ else:
369
+ print("architecture reference check skipped: no reference checkpoint supplied", flush=True)
370
+ output_path = resolve_project_path(args.output_checkpoint, config_path)
371
+ save_checkpoint(output_path, checkpoint_model, optimizer, args.epochs, loss, config)
372
+ print(f"checkpoint saved: {output_path}", flush=True)
373
+ print("training completed successfully", flush=True)
374
+
375
+
376
+ if __name__ == "__main__":
377
+ try:
378
+ main()
379
+ except Exception:
380
+ traceback.print_exc()
381
+ raise
382
+ finally:
383
+ if dist.is_available() and dist.is_initialized():
384
+ dist.destroy_process_group()
weight/.gitkeep ADDED
@@ -0,0 +1 @@
 
 
1
+ # Reserved for model weights.