File size: 19,996 Bytes
b871dba
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
"""Shared data, metric, configuration, and serialization utilities."""

from __future__ import annotations

import json
import os
import random
from pathlib import Path
from typing import Any, Dict, Iterable, Mapping, Optional, Sequence, Tuple

import numpy as np
import torch
import yaml
from torch import Tensor, nn
from torch.utils.data import DataLoader, Dataset


PAPER_SAMPLE_COUNT = 2595
AVAILABLE_SAMPLE_COUNT = 2215


def load_config(config_path: Path) -> Dict[str, Any]:
    """Load and minimally validate the experiment YAML."""
    config_path = Path(config_path).expanduser().resolve()
    if not config_path.is_file():
        raise FileNotFoundError(f"Configuration file not found: {config_path}")
    with config_path.open("r", encoding="utf-8") as handle:
        config = yaml.safe_load(handle)
    if not isinstance(config, dict):
        raise ValueError(f"Configuration must be a mapping: {config_path}")
    required_sections = ("experiment", "paths", "data", "model", "training", "evaluation")
    missing = [section for section in required_sections if section not in config]
    if missing:
        raise ValueError(f"Configuration is missing sections: {missing}")

    data_config = config["data"]
    if data_config.get("coordinate_normalization") != "none":
        raise ValueError("Paper fidelity requires coordinate_normalization: none")
    if data_config.get("target_normalization") != "train_minmax":
        raise ValueError("Paper fidelity requires target_normalization: train_minmax")
    if list(data_config.get("target_names", [])) != ["u", "v", "p"]:
        raise ValueError("Model target order must be [u, v, p]")
    if config["training"].get("scheduler") != "none":
        raise ValueError("The paper does not specify a learning-rate scheduler")
    return config


def resolve_path(project_root: Path, configured_path: str) -> Path:
    """Resolve an absolute path or a path relative to the project root."""
    path = Path(configured_path).expanduser()
    return path.resolve() if path.is_absolute() else (project_root / path).resolve()


def configured_paths(config: Mapping[str, Any], project_root: Path) -> Dict[str, Path]:
    """Resolve all data and output paths without consulting environment variables."""
    paths = config["paths"]
    data_dir = resolve_path(project_root, str(paths["data_dir"]))
    return {
        "data": data_dir / str(paths["data_file"]),
        "train_indices": data_dir / str(paths["train_indices"]),
        "validation_indices": data_dir / str(paths["validation_indices"]),
        "test_indices": data_dir / str(paths["test_indices"]),
        "checkpoint": resolve_path(project_root, str(paths["checkpoint"])),
        "results_dir": resolve_path(project_root, str(paths["results_dir"])),
    }


def set_deterministic_seed(seed: int) -> None:
    """Seed every RNG used by this project."""
    random.seed(seed)
    np.random.seed(seed)
    torch.manual_seed(seed)
    if torch.cuda.is_available():
        torch.cuda.manual_seed_all(seed)
    if hasattr(torch.backends, "cudnn"):
        torch.backends.cudnn.deterministic = True
        torch.backends.cudnn.benchmark = False
    try:
        torch.use_deterministic_algorithms(True, warn_only=True)
    except (AttributeError, TypeError):
        pass


def choose_device(requested: str) -> torch.device:
    """Resolve ``auto``, CPU, or an explicit CUDA device."""
    requested = requested.strip().lower()
    if requested == "auto":
        return torch.device("cuda" if torch.cuda.is_available() else "cpu")
    device = torch.device(requested)
    if device.type == "cuda" and not torch.cuda.is_available():
        raise RuntimeError(f"CUDA device requested but CUDA is unavailable: {requested}")
    return device


def _load_index_array(path: Path) -> np.ndarray:
    if not path.is_file():
        raise FileNotFoundError(f"Split index file not found: {path}")
    indices = np.load(path, allow_pickle=False)
    if indices.ndim != 1 or not np.issubdtype(indices.dtype, np.integer):
        raise ValueError(f"Split indices must be a one-dimensional integer array: {path}")
    return np.asarray(indices, dtype=np.int64)


def load_data_and_splits(
    config: Mapping[str, Any], project_root: Path
) -> Tuple[np.ndarray, Dict[str, np.ndarray], Dict[str, Path]]:
    """Load the CFD array read-only and validate the fixed supplied splits."""
    paths = configured_paths(config, project_root)
    if not paths["data"].is_file():
        raise FileNotFoundError(f"CFD data file not found: {paths['data']}")
    data = np.load(paths["data"], mmap_mode="r", allow_pickle=False)
    expected_points = int(config["data"]["num_points"])
    expected_channels = len(config["data"]["source_channels"])
    if data.ndim != 3 or data.shape[1:] != (expected_points, expected_channels):
        raise ValueError(
            f"Expected CFD array [cases,{expected_points},{expected_channels}], "
            f"got {data.shape}"
        )
    if data.dtype != np.float32:
        raise ValueError(f"Expected float32 CFD data, got {data.dtype}")
    if not np.isfinite(data).all():
        raise ValueError("CFD data contains NaN or infinite values")

    splits = {
        "train": _load_index_array(paths["train_indices"]),
        "validation": _load_index_array(paths["validation_indices"]),
        "test": _load_index_array(paths["test_indices"]),
    }
    total_cases = int(data.shape[0])
    for name, indices in splits.items():
        if indices.size == 0:
            raise ValueError(f"Split {name} is empty")
        if indices.min() < 0 or indices.max() >= total_cases:
            raise ValueError(f"Split {name} contains out-of-range indices")
        if np.unique(indices).size != indices.size:
            raise ValueError(f"Split {name} contains duplicate indices")
    split_sets = {name: set(values.tolist()) for name, values in splits.items()}
    if split_sets["train"] & split_sets["validation"]:
        raise ValueError("Training and validation splits overlap")
    if split_sets["train"] & split_sets["test"]:
        raise ValueError("Training and test splits overlap")
    if split_sets["validation"] & split_sets["test"]:
        raise ValueError("Validation and test splits overlap")
    covered = split_sets["train"] | split_sets["validation"] | split_sets["test"]
    if covered != set(range(total_cases)):
        raise ValueError("Fixed splits do not cover the dataset exactly")
    return data, splits, paths


def fit_target_minmax(
    data: np.ndarray, train_indices: np.ndarray, target_indices: Sequence[int]
) -> Tuple[np.ndarray, np.ndarray]:
    """Fit per-variable extrema using the training split and no other cases."""
    training_cases = np.asarray(data[train_indices], dtype=np.float32)
    training_targets = np.take(training_cases, target_indices, axis=-1)
    target_min = training_targets.min(axis=(0, 1)).astype(np.float32)
    target_max = training_targets.max(axis=(0, 1)).astype(np.float32)
    if not np.isfinite(target_min).all() or not np.isfinite(target_max).all():
        raise ValueError("Target normalization extrema are not finite")
    if np.any(target_max <= target_min):
        raise ValueError(
            f"Every target must have positive min-max span: min={target_min}, max={target_max}"
        )
    return target_min, target_max


class PointCFDDataset(Dataset):
    """A fixed case-index view of the supplied PointCFD NumPy array."""

    def __init__(
        self,
        data: np.ndarray,
        case_indices: np.ndarray,
        input_indices: Sequence[int],
        target_indices: Sequence[int],
        target_min: np.ndarray,
        target_max: np.ndarray,
    ) -> None:
        self.data = data
        self.case_indices = np.asarray(case_indices, dtype=np.int64)
        self.input_indices = tuple(int(index) for index in input_indices)
        self.target_indices = tuple(int(index) for index in target_indices)
        self.target_min = np.asarray(target_min, dtype=np.float32)
        self.target_max = np.asarray(target_max, dtype=np.float32)
        if self.target_min.shape != (len(self.target_indices),):
            raise ValueError("target_min has the wrong shape")
        if self.target_max.shape != self.target_min.shape:
            raise ValueError("target_max has the wrong shape")
        self.target_span = self.target_max - self.target_min
        if np.any(self.target_span <= 0):
            raise ValueError("Target min-max span must be positive")

    def __len__(self) -> int:
        return int(self.case_indices.size)

    def __getitem__(self, item: int) -> Tuple[Tensor, Tensor, Tensor]:
        case_index = int(self.case_indices[item])
        sample = self.data[case_index]
        coordinates = np.ascontiguousarray(sample[:, self.input_indices], dtype=np.float32)
        targets = np.ascontiguousarray(sample[:, self.target_indices], dtype=np.float32)
        normalized = np.ascontiguousarray(
            (targets - self.target_min) / self.target_span, dtype=np.float32
        )
        return (
            torch.from_numpy(coordinates),
            torch.from_numpy(normalized),
            torch.tensor(case_index, dtype=torch.int64),
        )


def make_loader(
    dataset: Dataset,
    batch_size: int,
    shuffle: bool,
    num_workers: int,
    seed: int,
    pin_memory: bool,
) -> DataLoader:
    if batch_size <= 0:
        raise ValueError("batch_size must be positive")
    if num_workers < 0:
        raise ValueError("num_workers cannot be negative")
    generator = torch.Generator()
    generator.manual_seed(seed)

    def seed_worker(worker_id: int) -> None:
        worker_seed = (seed + worker_id) % (2**32)
        np.random.seed(worker_seed)
        random.seed(worker_seed)

    return DataLoader(
        dataset,
        batch_size=batch_size,
        shuffle=shuffle,
        num_workers=num_workers,
        pin_memory=pin_memory,
        drop_last=False,
        generator=generator,
        worker_init_fn=seed_worker if num_workers else None,
    )


def inverse_target_minmax(
    normalized: np.ndarray, target_min: np.ndarray, target_max: np.ndarray
) -> np.ndarray:
    return normalized * (target_max - target_min) + target_min


def compute_field_metrics(
    predictions_normalized: np.ndarray,
    targets_normalized: np.ndarray,
    target_min: np.ndarray,
    target_max: np.ndarray,
    target_names: Sequence[str],
    relative_l2_epsilon: float,
) -> Tuple[Dict[str, Any], np.ndarray, np.ndarray]:
    """Compute paper-style field error metrics after inverse normalization."""
    predictions_normalized = np.asarray(predictions_normalized, dtype=np.float32)
    targets_normalized = np.asarray(targets_normalized, dtype=np.float32)
    if predictions_normalized.shape != targets_normalized.shape:
        raise ValueError("Prediction and target shapes do not match")
    if predictions_normalized.ndim != 3:
        raise ValueError("Evaluation arrays must be [cases, points, variables]")
    if predictions_normalized.shape[-1] != len(target_names):
        raise ValueError("Target names do not match the evaluated variables")
    if relative_l2_epsilon <= 0:
        raise ValueError("relative_l2_epsilon must be positive")

    normalized_error = predictions_normalized - targets_normalized
    normalized_mse = float(np.mean(np.square(normalized_error), dtype=np.float64))
    predictions = inverse_target_minmax(predictions_normalized, target_min, target_max)
    targets = inverse_target_minmax(targets_normalized, target_min, target_max)
    error = predictions - targets

    rmse: Dict[str, float] = {}
    relative_l2: Dict[str, Dict[str, float]] = {}
    for channel, name in enumerate(target_names):
        channel_error = error[:, :, channel].astype(np.float64)
        channel_target = targets[:, :, channel].astype(np.float64)
        rmse[str(name)] = float(np.sqrt(np.mean(np.square(channel_error))))
        numerator = np.linalg.norm(channel_error, axis=1)
        denominator = np.maximum(
            np.linalg.norm(channel_target, axis=1), relative_l2_epsilon
        )
        per_case = numerator / denominator
        relative_l2[str(name)] = {
            "mean": float(np.mean(per_case)),
            "max": float(np.max(per_case)),
            "min": float(np.min(per_case)),
        }
    metrics = {
        "normalized_mse": normalized_mse,
        "rmse": rmse,
        "relative_l2": relative_l2,
    }
    return metrics, predictions.astype(np.float32), targets.astype(np.float32)


def evaluate_model(
    model: nn.Module,
    loader: DataLoader,
    device: torch.device,
    target_min: np.ndarray,
    target_max: np.ndarray,
    target_names: Sequence[str],
    relative_l2_epsilon: float,
) -> Tuple[Dict[str, Any], Dict[str, np.ndarray]]:
    """Evaluate a complete split and retain arrays needed for reporting."""
    model.eval()
    coordinates_batches = []
    prediction_batches = []
    target_batches = []
    index_batches = []
    with torch.no_grad():
        for coordinates, targets, case_indices in loader:
            predictions = model(coordinates.to(device, non_blocking=True))
            coordinates_batches.append(coordinates.numpy())
            prediction_batches.append(predictions.cpu().numpy())
            target_batches.append(targets.numpy())
            index_batches.append(case_indices.numpy())
    if not prediction_batches:
        raise ValueError("Evaluation loader produced no batches")
    coordinates_array = np.concatenate(coordinates_batches, axis=0)
    predictions_normalized = np.concatenate(prediction_batches, axis=0)
    targets_normalized = np.concatenate(target_batches, axis=0)
    case_indices_array = np.concatenate(index_batches, axis=0)
    metrics, predictions, targets = compute_field_metrics(
        predictions_normalized,
        targets_normalized,
        target_min,
        target_max,
        target_names,
        relative_l2_epsilon,
    )
    arrays = {
        "coordinates": coordinates_array.astype(np.float32),
        "predictions": predictions,
        "targets": targets,
        "predictions_normalized": predictions_normalized.astype(np.float32),
        "targets_normalized": targets_normalized.astype(np.float32),
        "case_indices": case_indices_array.astype(np.int64),
    }
    return metrics, arrays


def _json_safe(value: Any) -> Any:
    if isinstance(value, Mapping):
        return {str(key): _json_safe(item) for key, item in value.items()}
    if isinstance(value, (list, tuple)):
        return [_json_safe(item) for item in value]
    if isinstance(value, Path):
        return str(value)
    if isinstance(value, np.ndarray):
        return value.tolist()
    if isinstance(value, np.generic):
        return value.item()
    if isinstance(value, torch.device):
        return str(value)
    if isinstance(value, float) and not np.isfinite(value):
        raise ValueError("JSON output cannot contain NaN or Infinity")
    return value


def write_json(path: Path, payload: Mapping[str, Any]) -> None:
    """Atomically write strict, human-readable JSON."""
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    temporary = path.with_name(path.name + ".tmp")
    with temporary.open("w", encoding="utf-8") as handle:
        json.dump(_json_safe(payload), handle, indent=2, sort_keys=True, allow_nan=False)
        handle.write("\n")
    os.replace(temporary, path)


def append_jsonl(path: Path, payload: Mapping[str, Any]) -> None:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    with path.open("a", encoding="utf-8") as handle:
        handle.write(json.dumps(_json_safe(payload), sort_keys=True, allow_nan=False) + "\n")
        handle.flush()


def write_npz(path: Path, **arrays: np.ndarray) -> None:
    """Atomically write a compressed NumPy archive."""
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    temporary = path.with_name(path.name + ".tmp.npz")
    np.savez_compressed(temporary, **arrays)
    os.replace(temporary, path)


def save_checkpoint(path: Path, payload: Mapping[str, Any]) -> None:
    path = Path(path)
    path.parent.mkdir(parents=True, exist_ok=True)
    temporary = path.with_name(path.name + ".tmp")
    torch.save(dict(payload), temporary)
    os.replace(temporary, path)


def load_checkpoint(path: Path, device: torch.device) -> Dict[str, Any]:
    path = Path(path).expanduser().resolve()
    if not path.is_file():
        raise FileNotFoundError(f"Checkpoint not found: {path}")
    try:
        checkpoint = torch.load(path, map_location=device, weights_only=False)
    except TypeError:
        checkpoint = torch.load(path, map_location=device)
    if not isinstance(checkpoint, dict) or "model_state_dict" not in checkpoint:
        raise ValueError(f"Unsupported checkpoint payload: {path}")
    return checkpoint


def rng_state() -> Dict[str, Any]:
    state: Dict[str, Any] = {
        "python": random.getstate(),
        "numpy": np.random.get_state(),
        "torch": torch.get_rng_state(),
    }
    if torch.cuda.is_available():
        state["cuda"] = torch.cuda.get_rng_state_all()
    return state


def restore_rng_state(state: Optional[Mapping[str, Any]]) -> None:
    if not state:
        return
    random.setstate(state["python"])
    np.random.set_state(state["numpy"])
    torch.set_rng_state(state["torch"])
    if torch.cuda.is_available() and "cuda" in state:
        torch.cuda.set_rng_state_all(state["cuda"])


def checkpoint_metadata_matches(
    checkpoint: Mapping[str, Any],
    config: Mapping[str, Any],
    target_min: np.ndarray,
    target_max: np.ndarray,
) -> None:
    """Reject silent channel or normalization changes on resume."""
    expected = {
        "source_channels": list(config["data"]["source_channels"]),
        "input_names": list(config["data"]["input_names"]),
        "target_names": list(config["data"]["target_names"]),
    }
    for key, value in expected.items():
        if list(checkpoint.get(key, [])) != value:
            raise ValueError(f"Checkpoint {key} metadata does not match the configuration")
    if not np.allclose(np.asarray(checkpoint.get("target_min")), target_min, rtol=0, atol=0):
        raise ValueError("Checkpoint target_min differs from the current training split")
    if not np.allclose(np.asarray(checkpoint.get("target_max")), target_max, rtol=0, atol=0):
        raise ValueError("Checkpoint target_max differs from the current training split")


def selected_indices(indices: np.ndarray, maximum_cases: Optional[int]) -> np.ndarray:
    if maximum_cases is None:
        return indices
    if maximum_cases <= 0:
        raise ValueError("maximum_cases must be positive")
    return indices[: min(maximum_cases, indices.size)]


def prepare_datasets(
    config: Mapping[str, Any],
    project_root: Path,
    smoke_test: bool = False,
) -> Tuple[Dict[str, PointCFDDataset], np.ndarray, np.ndarray, Dict[str, Path], Dict[str, int]]:
    """Build all split datasets under one normalization contract."""
    data, splits, paths = load_data_and_splits(config, project_root)
    if smoke_test:
        active_splits = {
            "train": selected_indices(splits["train"], 4),
            "validation": selected_indices(splits["validation"], 2),
            "test": selected_indices(splits["test"], 2),
        }
    else:
        active_splits = splits
    target_indices = tuple(int(i) for i in config["data"]["target_indices"])
    target_min, target_max = fit_target_minmax(data, active_splits["train"], target_indices)
    datasets = {
        name: PointCFDDataset(
            data,
            indices,
            config["data"]["input_indices"],
            target_indices,
            target_min,
            target_max,
        )
        for name, indices in active_splits.items()
    }
    counts = {name: len(dataset) for name, dataset in datasets.items()}
    return datasets, target_min, target_max, paths, counts