diff --git a/.dockerignore b/.dockerignore index 92d233b..4a4dbd4 100644 --- a/.dockerignore +++ b/.dockerignore @@ -28,19 +28,6 @@ venv/ .venv/ .env -# contrail_ml and contrail_flights are excluded entirely — ML training/inference -# and the OpenSky real-traffic loader are separate concerns with their own heavy -# deps (scipy/sklearn/xgboost/mlflow; pyopensky/pandas/pyarrow). The server image -# has one job: run the CP-SAT gRPC service. Both run offline in their own envs. -contrail_ml/ -contrail_flights/ - -# contrail_ml artifacts -mlruns/ -data/processed/ -*.joblib -*.npz - # Editor / OS noise and local notes .DS_Store *.ipynb_checkpoints/ diff --git a/.env.example b/.env.example deleted file mode 100644 index 711a7fe..0000000 --- a/.env.example +++ /dev/null @@ -1,21 +0,0 @@ -# contrail_ml configuration (12-factor). Copy to `.env` and fill in as needed. -# Every MLConfig field can be overridden via a CONTRAIL_ML_ variable. -# NO secrets belong in the repo — .env is gitignored. - -# Folder of IAGOS per-flight NetCDF files (the ground-truth labels). -CONTRAIL_ML_IAGOS_DIR= - -# MLflow tracking backend. Local file store by default; point at a server URI -# (e.g. http://mlflow.internal:5000) to share runs with a team. -CONTRAIL_ML_MLFLOW_TRACKING_URI=file:./mlruns - -# Geographic region (degrees) and cruise pressure levels (hPa, comma-separated). -CONTRAIL_ML_LON_MIN=-30 -CONTRAIL_ML_LON_MAX=30 -CONTRAIL_ML_LAT_MIN=30 -CONTRAIL_ML_LAT_MAX=70 -CONTRAIL_ML_PRESSURE_LEVELS_HPA=150,175,200,225,250,300 - -# Geographic anchor: local sim (x=0, y=0) -> (lat, lon). -CONTRAIL_ML_ORIGIN_LAT=43.0 -CONTRAIL_ML_ORIGIN_LON=-5.0 diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a14ce3f..d1edabd 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -27,41 +27,11 @@ jobs: run: ruff check . - name: Mypy (type-check) - run: mypy contrail_env service contrail_flights + run: mypy contrail_env service - name: Pytest run: pytest -q - ml: - # The contrail_ml ISSR model + MLOps lifecycle. Runs the HERMETIC ML tests - # (synthetic fallback only — no network, no credentials). Separate job so a - # heavy [ml] install never slows the core lint-type-test gate, and the ML - # deps never touch the runtime Docker image. - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - - uses: actions/setup-python@v5 - with: - python-version: "3.11" - - - name: Install package with ML + dev extras - run: | - python -m pip install --upgrade pip - pip install -e ".[ml,dev]" - - - name: Generate gRPC stubs - run: bash scripts/gen_proto.sh - - - name: Ruff (lint contrail_ml) - run: ruff check contrail_ml - - - name: Mypy (type-check contrail_ml) - run: mypy contrail_ml - - - name: Pytest (hermetic ML tests) - run: pytest tests/test_ml_*.py -q - docker-build: runs-on: ubuntu-latest steps: @@ -88,8 +58,8 @@ jobs: docker-publish: runs-on: ubuntu-latest # Only publish on pushes to main — not on pull requests. - # The full suite (core + ML) AND the image smoke-test must pass first. - needs: [lint-type-test, ml, docker-build] + # The full suite AND the image smoke-test must pass first. + needs: [lint-type-test, docker-build] if: github.event_name == 'push' && github.ref == 'refs/heads/main' permissions: contents: read diff --git a/.gitignore b/.gitignore index 5e019f3..2df127c 100644 --- a/.gitignore +++ b/.gitignore @@ -17,10 +17,5 @@ build/ # Generated gRPC stubs (regenerated by scripts/gen_proto.sh) service/generated/ -# contrail_ml — MLflow runs, versioned datasets, trained artifacts (all rebuilt) -mlruns/ -data/processed/ -*.joblib -*.npz -*.dvc -.dvc/ +# Benchmark output +results.csv diff --git a/Dockerfile b/Dockerfile index c21d730..e9c8f30 100644 --- a/Dockerfile +++ b/Dockerfile @@ -32,12 +32,10 @@ COPY pyproject.toml README.md ./ COPY contrail_env/ ./contrail_env/ COPY service/ ./service/ COPY scripts/ ./scripts/ -# contrail_ml and the quantum solvers (pasqal_analog, xanadu_gbs) are -# intentionally excluded. This image has one job: run the CP-SAT gRPC server. -# ML training/inference and quantum sampling are separate concerns that carry -# heavy deps (scipy/sklearn/xgboost/pulser/strawberryfields) and run in their -# own environments. Including their source here would ship dead code that can -# never execute (deps absent) and would crash at request time if called. +# Only contrail_env + service ship here — the image's one job is the CP-SAT gRPC +# server. The quantum solver modules (pasqal_analog, xanadu_gbs) ride along in +# contrail_env and run on their built-in fallbacks; the optional [quantum] SDKs +# (pulser/strawberryfields) are not installed, keeping the runtime image lean. RUN pip install . # The gRPC stubs are gitignored — generate them from solver.proto at build time. diff --git a/README.md b/README.md index bdf780c..478c820 100644 --- a/README.md +++ b/README.md @@ -84,7 +84,7 @@ The dashboard has six tabs: live CP-SAT convergence (over ZMQ), the conflict-gra topology, QUBO matrix statistics (size, sparsity, penalty constants), the chosen-option trade-offs, the quantum benchmark (CP-SAT vs Pasqal vs Xanadu over N seeds, with live convergence curves for the BO loop and the GBS sampler), and a -geographic map — the predicted ISSR risk as a marker overlay on a real Plotly +geographic map — the ISSR risk as a marker overlay on a real Plotly `geo` basemap (country borders / coastlines, drawn with SVG and bundled offline vectors, so it needs no WebGL or network) with the chosen vs context routes on top. @@ -94,51 +94,24 @@ The benchmark also runs headless: python -m contrail_env.benchmark --flights 4 --seeds 5 --csv results.csv ``` -## ISSR model — real weather instead of synthetic blobs +## ISSR field -By default the airspace's contrail zones are synthetic Gaussian blobs. The -`contrail_ml` package replaces them with a trained model that predicts -ice-supersaturated regions (ISSRs) from real weather, bias-corrected against -in-situ IAGOS humidity, and wraps it in a full MLOps lifecycle (versioned data → -train → calibrate → register with MLflow → serve → monitor). - -The model plugs in through the **same `ISSRField` interface** the synthetic -field uses, so `World` and the QUBO assembly are untouched — you flip a switch: - -```python -from contrail_env import default_european_world -world = default_european_world(issr_source="ml", issr_kwargs={...}) # vs "synthetic" -``` - -Over gRPC the `ScenarioConfig` gained an `issr_source` field (default -`"synthetic"`, so existing clients are unaffected). Install the extra and try the -whole serving seam offline: - -``` -pip install -e ".[ml]" -python -m contrail_ml serve-check # ML ISSR field -> CP-SAT solve -python -m contrail_ml train --synthetic --no-mlflow # CV, fit, calibrate, baseline table -``` - -The science, the model, and how to read the model-vs-baselines table are in -[docs/ML.md](docs/ML.md); data sources (IAGOS / ARCO-ERA5 / GFS) in -[docs/DATA.md](docs/DATA.md). The hermetic tests use a guarded synthetic -fallback — no network, no credentials. +The airspace's contrail zones (ice-supersaturated regions) are synthetic +Gaussian blobs — a controllable obstacle field for the optimizer to route +around. The field is consumed through a small `ISSRField` interface +(`rhi_excess`, `is_inside`, `mask_grid`), so the *source* of the field is +pluggable without touching `World` or the QUBO assembly. ## Layout ``` -contrail_env/ synthetic environment, QUBO assembly, CP-SAT solver, - quantum pipelines (pasqal_analog, xanadu_gbs, quantum_common, - bayes_opt) and the benchmark protocol (benchmark.py) -contrail_ml/ ISSR model (features, RHiCorrector, calibration), MLflow - registry, serving (MLIssrField), monitoring, and the data - loaders (IAGOS/ERA5/GFS) — the [ml] extra +contrail_env/ synthetic environment, ISSR field, geo anchor, QUBO assembly, + CP-SAT solver, quantum pipelines (pasqal_analog, xanadu_gbs, + quantum_common, bayes_opt) and the benchmark protocol (benchmark.py) service/ gRPC service, ZMQ progress streaming, client gui/ PyQt6 dashboard tests/ environment build, CP-SAT vs brute-force, quantum solvers vs - brute-force, benchmark round-trip, gRPC round-trip, and the - hermetic contrail_ml suite (test_ml_*.py) + brute-force, benchmark round-trip, gRPC round-trip, GUI map panel ``` ## Development @@ -154,11 +127,6 @@ CI runs the same checks on every push and pull request. ## Roadmap -- **Done:** real ISSR model (`contrail_ml`) replacing the synthetic field, with - the train → calibrate → register → serve → monitor MLOps lifecycle (see - [docs/ML.md](docs/ML.md)). Next: run it end-to-end on real IAGOS + ERA5 (the - loaders are written; they need portal registration) and report the honest - real-data comparison table. - Run the Pasqal pipeline on real Pulser hardware: needs `[quantum]` extras plus a conflict graph that embeds as a valid unit-disk register (auto-detected; the built-in simulator is the fallback). diff --git a/contrail_flights/geo.py b/contrail_env/geo.py similarity index 55% rename from contrail_flights/geo.py rename to contrail_env/geo.py index 02a6fc9..f997e78 100644 --- a/contrail_flights/geo.py +++ b/contrail_env/geo.py @@ -1,20 +1,13 @@ """ -geo.py — local sim frame <-> real (lon, lat), for the real-flight loader. +geo.py — Where the synthetic world sits on Earth (local sim frame <-> lon/lat). -This is a DELIBERATE, independent mirror of -`contrail_ml.features.GeoAnchor`: the brief requires `contrail_flights` to be -a standalone sibling of `contrail_ml` (no import between the two), so the -transform is duplicated here rather than shared. The two MUST stay in sync — -same formula, same default origin — so the predicted ISSR field (placed via -the contrail_ml anchor) and the real flights (placed via this anchor) agree on -where things are on the map. `tests/test_geo_transform_roundtrip.py` checks the -inverse, and a cross-check test guards against drift from the contrail_ml copy. +The optimizer works entirely in the local 1500x800 km Cartesian box; geography +only matters for the GUI Map tab, which needs to place that box on a real +basemap. This module owns that one mapping so the map has a home for it without +reaching into any heavier package. lat = origin_lat + y_km / KM_PER_DEG_LAT lon = origin_lon + x_km / (KM_PER_DEG_LAT * cos(lat)) - -Only numpy is needed, so this module imports cleanly in the lean install (the -[flights] extra is only for the actual OpenSky pull, not the geometry). """ from __future__ import annotations @@ -23,20 +16,26 @@ import numpy as np -# Must match contrail_ml.config.MLConfig defaults so flights and the predicted -# field share one coordinate system. -DEFAULT_ORIGIN_LAT = 43.0 -DEFAULT_ORIGIN_LON = -5.0 -KM_PER_DEG_LAT = 111.0 + +def _maybe_scalar(arr: np.ndarray, like: np.ndarray | float) -> np.ndarray | float: + """Return a Python float when the input was scalar, else the array.""" + if np.isscalar(like) or (isinstance(like, np.ndarray) and like.ndim == 0): + return float(arr) + return arr @dataclass(frozen=True) class GeoAnchor: - """Maps the local sim frame (x_km east, y_km north) to/from (lon, lat).""" + """Maps the local sim frame (x_km east, y_km north) to geography. - origin_lat: float = DEFAULT_ORIGIN_LAT - origin_lon: float = DEFAULT_ORIGIN_LON - km_per_deg_lat: float = KM_PER_DEG_LAT + The sim is a flat Cartesian box; the map lives on (lon, lat). One + small-angle anchor ties them together so the risk overlay and the routes + render over the same place. + """ + + origin_lat: float = 43.0 + origin_lon: float = -5.0 + km_per_deg_lat: float = 111.0 def local_to_geo( self, x_km: np.ndarray | float, y_km: np.ndarray | float @@ -59,8 +58,6 @@ def geo_to_local( return _maybe_scalar(x, lon_deg), _maybe_scalar(y, lat_deg) -def _maybe_scalar(arr: np.ndarray, like: np.ndarray | float) -> np.ndarray | float: - """Return a Python float when the input was scalar, else the array.""" - if np.isscalar(like) or (isinstance(like, np.ndarray) and like.ndim == 0): - return float(arr) - return arr +# The canonical anchor: places the 1500x800 km box over south-west -> central +# Europe (roughly Madrid to Frankfurt), matching the default world geometry. +EUROPEAN_ANCHOR = GeoAnchor(origin_lat=43.0, origin_lon=-5.0) diff --git a/contrail_env/world.py b/contrail_env/world.py index 23ae6c0..47595d2 100644 --- a/contrail_env/world.py +++ b/contrail_env/world.py @@ -178,14 +178,8 @@ def default_european_world( - Uniform 3x3 sector grid with capacity 3 issr_source: - "synthetic" (default) — the random Gaussian-blob ISSR field, exactly - as before, so all existing callers/tests are unaffected. - "ml" — a trained model's predicted field (contrail_ml.MLIssrField), - which honours the SAME ISSRField interface, so World/qubo need no - changes. contrail_ml is imported LAZILY here, so the core install - (no [ml] extra) still imports contrail_env fine as long as the ml - source is not requested. - `issr_kwargs` is forwarded to the chosen ISSR-field factory. + "synthetic" (default) — the random Gaussian-blob ISSR field. + `issr_kwargs` is forwarded to the ISSR-field factory. """ from .synthetic_issr import random_issr_field @@ -205,14 +199,9 @@ def default_european_world( seed=seed, **kw, ) - elif issr_source == "ml": - # Lazy import: only the "ml" branch pulls in the [ml] extra. - from contrail_ml.issr_field import ml_issr_field - - issr = ml_issr_field(**kw) else: raise ValueError( - f"unknown issr_source {issr_source!r} (use 'synthetic' or 'ml')" + f"unknown issr_source {issr_source!r} (use 'synthetic')" ) from .airspace import uniform_sector_grid diff --git a/contrail_flights/__init__.py b/contrail_flights/__init__.py deleted file mode 100644 index b40a3d5..0000000 --- a/contrail_flights/__init__.py +++ /dev/null @@ -1,46 +0,0 @@ -""" -contrail_flights — real historical traffic (OpenSky) as a Flight source. - -An alternative to the synthetic `build_random_flights`: pull actual flown -ADS-B tracks for a European region + time window from the OpenSky Network and -reduce each to a `contrail_env.Flight`, so the existing option enumeration, -QUBO builder, all three solvers, and the GUI map consume them unchanged. - -DESIGN INVARIANTS -================= -1. Real data only — never fabricate traffic. Missing deps/credentials/network, - or an empty result, raise `OpenSkyUnavailableError` (see opensky_client). -2. Independent of contrail_ml: this package is a sibling, not a dependency. It - carries its own geo transform (geo.py), kept in sync with the model's anchor. -3. Lean import surface: only numpy + contrail_env at import time. pyopensky / - pandas / pyarrow (the [flights] extra) are imported lazily where used. -4. Default stays synthetic: the service only builds real flights when explicitly - asked (flight_source="real"). - -HONEST NAMING -============= -OpenSky gives real flown HISTORICAL tracks, not published schedules. This is -"demonstrated on real historical European traffic" — the same evaluation method -used in real contrail-avoidance trials — NOT "schedule optimization". -""" - -from __future__ import annotations - -from .config import DEFAULT_CONFIG, FlightsConfig -from .geo import GeoAnchor - -__all__ = ["FlightsConfig", "DEFAULT_CONFIG", "GeoAnchor", "flights_for", "reduce_to_flight"] - - -def __getattr__(name: str): - # Lazy re-export so the common helpers are importable from the package root - # without importing the heavier orchestration modules at package load. - if name == "flights_for": - from .build_dataset import flights_for - - return flights_for - if name == "reduce_to_flight": - from .reduce_to_flight import reduce_to_flight - - return reduce_to_flight - raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/contrail_flights/build_dataset.py b/contrail_flights/build_dataset.py deleted file mode 100644 index 2f2ee83..0000000 --- a/contrail_flights/build_dataset.py +++ /dev/null @@ -1,100 +0,0 @@ -""" -build_dataset.py — orchestrate: OpenSky -> cleaned tracks -> list[Flight]. - -`flights_for` is the one entry point the service/CLI call: it pulls (or loads -from cache) the tracks for the configured region+window, cleans them, and -reduces each to a `contrail_env.Flight`, all anchored to one shared snapshot -window so the flights genuinely compete for the same airspace. - -Real European traffic naturally overlaps (it follows a small set of airways), -so — unlike the synthetic generator — no artificial "corridor" trick is needed. -""" - -from __future__ import annotations - -import logging - -from contrail_env import Flight - -from .config import FlightsConfig -from .reduce_to_flight import reduce_to_flight -from .tracks import Track, clean_tracks - -log = logging.getLogger(__name__) - - -def flights_for( - config: FlightsConfig, - *, - tracks: list[Track] | None = None, - use_cache: bool = True, - max_flights: int | None = None, -) -> list[Flight]: - """Build a list of real `Flight`s for `config`. - - Parameters - ---------- - tracks : pre-supplied cleaned tracks (used by tests/offline runs). When - None, tracks come from the cache or a live OpenSky pull. - use_cache : read/write the on-disk parquet cache around a live pull. - max_flights : cap the number of flights returned (largest tracks first). - - Raises OpenSkyUnavailableError (from the client) if a live pull is needed - but unavailable — never silently returns synthetic or empty data. - """ - raw = tracks if tracks is not None else _load_or_pull(config, use_cache=use_cache) - cleaned = clean_tracks(raw, config.min_track_points, config.min_track_km) - if not cleaned: - raise ValueError( - "no tracks survived cleaning — every track was too short/sparse for " - f"the thresholds (min_points={config.min_track_points}, " - f"min_track_km={config.min_track_km})." - ) - - # Largest (longest) tracks first, then optionally cap. - cleaned.sort(key=lambda t: t.ground_distance_km(), reverse=True) - if max_flights is not None: - cleaned = cleaned[:max_flights] - - # Shared reference epoch = earliest in-window time, so the first flight - # departs at t=0 and the rest are offset relative to it. - reference_epoch_s = min(float(t.time_s[0]) for t in cleaned) - - flights: list[Flight] = [] - n_fallback = 0 - for i, tr in enumerate(cleaned): - try: - reduced = reduce_to_flight( - tr, config, - name=f"R{i + 1}", - reference_epoch_s=reference_epoch_s, - ) - except ValueError as exc: - log.warning("skipping track %s: %s", tr.icao24, exc) - continue - if not reduced.used_observed_baseline: - n_fallback += 1 - flights.append(reduced.flight) - - if n_fallback: - log.info("%d/%d flights used the synthetic baseline fallback " - "(observed altitude too sparse)", n_fallback, len(flights)) - if not flights: - raise ValueError("no usable flights after reduction (all tracks skipped).") - return flights - - -def _load_or_pull(config: FlightsConfig, *, use_cache: bool) -> list[Track]: - """Return raw tracks from cache, else pull from OpenSky and cache them.""" - from .cache import TrackCache - from .opensky_client import OpenSkyClient - - cache = TrackCache(config) - if use_cache and cache.exists(): - log.info("loading tracks from cache %s", cache.path) - return cache.load() - - tracks = OpenSkyClient(config).fetch_tracks() - if use_cache: - cache.save(tracks) - return tracks diff --git a/contrail_flights/cache.py b/contrail_flights/cache.py deleted file mode 100644 index 397fbf0..0000000 --- a/contrail_flights/cache.py +++ /dev/null @@ -1,81 +0,0 @@ -""" -cache.py — local parquet cache for pulled OpenSky tracks. - -OpenSky bulk-history queries are rate-limited and credit-metered, so we cache -the cleaned tracks on disk keyed by (bbox, time window). Re-running the same -scenario then hits the cache instead of the network. - -pandas/pyarrow are in the [flights] extra and imported lazily, so this module -imports in the lean install; only `load`/`save` touch the heavy deps. -""" - -from __future__ import annotations - -import hashlib -import os - -import numpy as np - -from .config import FlightsConfig -from .tracks import Track - -_COLUMNS = ("icao24", "callsign", "time_s", "lat", "lon", "baro_altitude_m") - - -def _key(config: FlightsConfig) -> str: - """Stable filename for a config's bbox + time window.""" - raw = f"{config.bbox}|{config.start_time}|{config.end_time}" - return hashlib.sha1(raw.encode()).hexdigest()[:16] - - -class TrackCache: - """Reads/writes lists of `Track` as a single long-format parquet file.""" - - def __init__(self, config: FlightsConfig) -> None: - self._config = config - self._path = os.path.join(config.cache_dir, f"tracks_{_key(config)}.parquet") - - @property - def path(self) -> str: - return self._path - - def exists(self) -> bool: - return os.path.exists(self._path) - - def load(self) -> list[Track]: - """Load cached tracks; raises FileNotFoundError if absent.""" - import pandas as pd - - df = pd.read_parquet(self._path) - tracks: list[Track] = [] - for (icao24, callsign), grp in df.groupby(["icao24", "callsign"], sort=False): - tracks.append(Track( - icao24=str(icao24), - callsign=str(callsign), - time_s=grp["time_s"].to_numpy(dtype=float), - lat=grp["lat"].to_numpy(dtype=float), - lon=grp["lon"].to_numpy(dtype=float), - baro_altitude_m=grp["baro_altitude_m"].to_numpy(dtype=float), - )) - return tracks - - def save(self, tracks: list[Track]) -> None: - """Write tracks to the cache (creating the directory if needed).""" - import pandas as pd - - os.makedirs(self._config.cache_dir, exist_ok=True) - frames = [] - for tr in tracks: - n = tr.n_points - frames.append(pd.DataFrame({ - "icao24": np.repeat(tr.icao24, n), - "callsign": np.repeat(tr.callsign, n), - "time_s": tr.time_s, - "lat": tr.lat, - "lon": tr.lon, - "baro_altitude_m": tr.baro_altitude_m, - })) - out = pd.concat(frames, ignore_index=True) if frames else pd.DataFrame( - columns=list(_COLUMNS) - ) - out.to_parquet(self._path, index=False) diff --git a/contrail_flights/config.py b/contrail_flights/config.py deleted file mode 100644 index 3a36488..0000000 --- a/contrail_flights/config.py +++ /dev/null @@ -1,109 +0,0 @@ -""" -config.py — Typed configuration for the real-flight (OpenSky) loader. - -One dataclass (`FlightsConfig`) carries the planning box, the time window, the -cleaning thresholds, the snapshot-window normalization, the cache directory, -and the OpenSky credentials (read from the environment — never hardcoded). The -geographic box and anchor default to MATCH `contrail_ml.config.MLConfig`, so a -real flight's local (x_km, y_km) lands in the same frame the predicted ISSR -field uses. -""" - -from __future__ import annotations - -import os -from dataclasses import dataclass, field, replace -from typing import Any - -from .geo import DEFAULT_ORIGIN_LAT, DEFAULT_ORIGIN_LON, GeoAnchor - - -@dataclass(frozen=True) -class FlightsConfig: - """All configuration for pulling, cleaning, and reducing real tracks.""" - - # ---- Planning box (deg). Matches the default sim grid extent (0..1500 km - # east, 0..800 km north) under the shared anchor below. ---------------- - lon_min: float = -5.0 - lon_max: float = 15.0 - lat_min: float = 43.0 - lat_max: float = 50.0 - - # ---- Geographic anchor: local (x=0, y=0) -> (origin_lon, origin_lat). - # MUST match contrail_ml.config.MLConfig so flights and field agree. ---- - origin_lat: float = DEFAULT_ORIGIN_LAT - origin_lon: float = DEFAULT_ORIGIN_LON - - # ---- Time window (ISO 8601) for the historical pull -------------------- - start_time: str = "2023-06-01T12:00:00Z" - end_time: str = "2023-06-01T13:00:00Z" - - # ---- Track cleaning ---------------------------------------------------- - min_track_points: int = 10 # drop sparse tracks - min_track_km: float = 200.0 # drop tracks that barely cross the box - - # ---- Departure normalization (mirrors build_random_flights) ----------- - # Real entry times span the whole query window; we compress them into a - # short snapshot window so flights compete for the same airspace. - snapshot_window_s: float = 300.0 - - # ---- Cache + credentials (read from env, never hardcoded) ------------- - cache_dir: str = "data/flights_cache" - opensky_username: str | None = None - opensky_password: str | None = None - - extra: dict[str, Any] = field(default_factory=dict) - - # ------------------------------------------------------------------ # - - @property - def anchor(self) -> GeoAnchor: - """The geographic anchor placing the local frame on (lon, lat).""" - return GeoAnchor(origin_lat=self.origin_lat, origin_lon=self.origin_lon) - - @property - def bbox(self) -> tuple[float, float, float, float]: - """(lon_min, lon_max, lat_min, lat_max) for the OpenSky query.""" - return (self.lon_min, self.lon_max, self.lat_min, self.lat_max) - - def in_bbox(self, lon: float, lat: float) -> bool: - """Is this (lon, lat) inside the planning box?""" - return (self.lon_min <= lon <= self.lon_max - and self.lat_min <= lat <= self.lat_max) - - def with_overrides(self, **kwargs: Any) -> FlightsConfig: - """Return a copy with the given fields replaced.""" - return replace(self, **kwargs) - - @classmethod - def from_env(cls, prefix: str = "CONTRAIL_FLIGHTS_") -> FlightsConfig: - """Build a config, overriding any field from `${PREFIX}FIELD` env vars. - - OpenSky credentials also fall back to the conventional OPENSKY_USERNAME - / OPENSKY_PASSWORD names so they line up with pyopensky's own config. - """ - base = cls() - overrides: dict[str, Any] = {} - for f in base.__dataclass_fields__.values(): # type: ignore[attr-defined] - if f.name == "extra": - continue - env_key = f"{prefix}{f.name.upper()}" - if env_key in os.environ: - overrides[f.name] = _coerce(getattr(base, f.name), os.environ[env_key]) - overrides.setdefault("opensky_username", os.environ.get("OPENSKY_USERNAME")) - overrides.setdefault("opensky_password", os.environ.get("OPENSKY_PASSWORD")) - return replace(base, **overrides) - - -def _coerce(reference: Any, raw: str) -> Any: - """Coerce an env-var string to the type of the dataclass default value.""" - if isinstance(reference, bool): - return raw.strip().lower() in ("1", "true", "yes", "on") - if isinstance(reference, int) and not isinstance(reference, bool): - return int(raw) - if isinstance(reference, float): - return float(raw) - return raw - - -DEFAULT_CONFIG = FlightsConfig() diff --git a/contrail_flights/opensky_client.py b/contrail_flights/opensky_client.py deleted file mode 100644 index cac45ef..0000000 --- a/contrail_flights/opensky_client.py +++ /dev/null @@ -1,111 +0,0 @@ -""" -opensky_client.py — thin wrapper around the OpenSky Network historical API. - -ACCESS NOTES (see docs/DATA.md) -=============================== -* The public OpenSky REST API works unauthenticated for light, RECENT queries - (a rolling ~1-2 hour window) but is heavily rate-limited. -* Bulk HISTORICAL queries (arbitrary past dates, whole regions) require an - OpenSky research/Trino account — apply via https://opensky-network.org. Put - the credentials in OPENSKY_USERNAME / OPENSKY_PASSWORD (or a pyopensky config - file); never hardcode them. - -HONESTY GUARANTEE -================= -This client NEVER fabricates traffic. If the `pyopensky` dependency is missing, -credentials are absent, or the network/query fails, it raises -`OpenSkyUnavailableError` with an actionable message. An empty result is also -surfaced as an error (a real region/time should have traffic) rather than -silently returning zero flights. -""" - -from __future__ import annotations - -from .config import FlightsConfig -from .tracks import Track, make_track - - -class OpenSkyUnavailableError(RuntimeError): - """Raised when real OpenSky data cannot be obtained (deps/creds/network).""" - - -class OpenSkyClient: - """Fetches and assembles raw `Track`s for a region + time window. - - The heavy `pyopensky` import is lazy and lives inside `_require_api`, so - importing this module is cheap and the lean install is unaffected until a - real pull is actually requested. - """ - - def __init__(self, config: FlightsConfig) -> None: - self._config = config - - def _require_api(self): - """Import pyopensky and build an authenticated handle, or fail loudly.""" - try: - from pyopensky.trino import Trino - except ImportError as exc: # dependency not installed - raise OpenSkyUnavailableError( - "the 'flights' extra is not installed — real OpenSky access needs " - "pyopensky. Install with: pip install -e '.[flights]'" - ) from exc - - cfg = self._config - if not (cfg.opensky_username and cfg.opensky_password): - raise OpenSkyUnavailableError( - "OpenSky credentials missing — set OPENSKY_USERNAME / " - "OPENSKY_PASSWORD (research/Trino access required for historical " - "bulk queries; see docs/DATA.md)." - ) - try: - return Trino() - except Exception as exc: # network / auth / config failure - raise OpenSkyUnavailableError( - f"could not connect to OpenSky Trino: {exc}" - ) from exc - - def fetch_tracks(self) -> list[Track]: - """Pull state vectors for the configured bbox+window and group them - into per-aircraft `Track`s. Raises OpenSkyUnavailableError on any - failure (including an empty result).""" - api = self._require_api() - cfg = self._config - try: - df = api.history( - start=cfg.start_time, - stop=cfg.end_time, - bounds=cfg.bbox, - ) - except Exception as exc: - raise OpenSkyUnavailableError(f"OpenSky history query failed: {exc}") from exc - - if df is None or len(df) == 0: - raise OpenSkyUnavailableError( - f"OpenSky returned no traffic for {cfg.bbox} in " - f"[{cfg.start_time}, {cfg.end_time}] — widen the window/box or " - f"check access (it should not be empty for a real region/time)." - ) - return tracks_from_state_dataframe(df) - - -def tracks_from_state_dataframe(df) -> list[Track]: - """Group an OpenSky state-vector dataframe into per-icao24 `Track`s. - - Expects the pyopensky column names (icao24, callsign, time, latitude, - longitude, baroaltitude). Kept separate from the network call so it can be - unit-tested with a small in-memory frame. - """ - tracks: list[Track] = [] - for icao24, grp in df.groupby("icao24"): - callsign = "" - if "callsign" in grp.columns and len(grp["callsign"]): - callsign = str(grp["callsign"].iloc[0] or "").strip() - tracks.append(make_track( - icao24=str(icao24), - callsign=callsign, - time_s=grp["time"].to_numpy(), - lat=grp["latitude"].to_numpy(), - lon=grp["longitude"].to_numpy(), - baro_altitude_m=grp["baroaltitude"].to_numpy(), - )) - return tracks diff --git a/contrail_flights/reduce_to_flight.py b/contrail_flights/reduce_to_flight.py deleted file mode 100644 index 3ae06c6..0000000 --- a/contrail_flights/reduce_to_flight.py +++ /dev/null @@ -1,176 +0,0 @@ -""" -reduce_to_flight.py — one cleaned ADS-B track -> one contrail_env.Flight. - -This is the seam that lets REAL historical traffic flow through the exact same -optimizer the synthetic generator feeds: the output is an ordinary -`contrail_env.Flight`, so option enumeration, the QUBO builder, all three -solvers, and the GUI map need no changes. - -Honesty notes -============= -* The baseline is PREFERABLY the real observed altitude profile (resampled - from the track's barometric altitude), so "baseline" means the actually-flown - plan. Only when the altitude data is too sparse/noisy do we fall back to the - synthetic fuel-optimal `build_baseline_profile`, and we say so (a flag on the - returned info, logged by build_dataset). -* Aircraft type -> performance profile is an explicit, documented lookup. We do - NOT invent a new performance model per type; unmapped types use a320_like. - -Pure numpy + contrail_env (core) — no pandas/network — so it imports lean. -""" - -from __future__ import annotations - -from dataclasses import dataclass - -import numpy as np - -from contrail_env import ( - Aircraft, - AltitudeProfile, - AltitudeSegment, - Flight, - a320_like, - build_baseline_profile, - m_to_fl, -) - -from .config import FlightsConfig -from .tracks import Track - -# Explicit aircraft-type -> Aircraft factory table. Only the a320_like linear -# performance model exists today, so every mapped type points at it; the table -# is the place to add real BADA-class profiles later. Keys are ICAO type -# designators (as OpenSky's aircraft DB reports them). -AIRCRAFT_TYPE_TO_FACTORY = { - "A319": a320_like, - "A320": a320_like, - "A321": a320_like, - "B737": a320_like, - "B738": a320_like, -} - -# FL snap grid (RVSM 2000-ft steps) and the physical envelope we clamp to. -_FL_STEP = 20 -_FL_MIN, _FL_MAX = 280, 420 - - -@dataclass(frozen=True) -class ReducedFlight: - """A Flight plus provenance: did we use the real profile or the fallback?""" - - flight: Flight - used_observed_baseline: bool - - -def aircraft_for(tail: str, aircraft_type: str | None) -> Aircraft: - """Map an ICAO aircraft type to an Aircraft, defaulting to a320_like.""" - factory = AIRCRAFT_TYPE_TO_FACTORY.get((aircraft_type or "").upper(), a320_like) - return factory(tail) - - -def reduce_to_flight( - track: Track, - config: FlightsConfig, - *, - name: str | None = None, - aircraft_type: str | None = None, - reference_epoch_s: float = 0.0, -) -> ReducedFlight: - """Convert one cleaned `Track` into a `contrail_env.Flight`. - - reference_epoch_s anchors departure times: departure_s = entry_time - - reference_epoch_s (clamped at 0). build_dataset passes the earliest entry - time across the batch so flights share one snapshot window. - - Raises ValueError if the track never enters the planning box (the caller - should have filtered, but we fail loudly rather than fabricate a position). - """ - anchor = config.anchor - - inside = ( - (track.lon >= config.lon_min) & (track.lon <= config.lon_max) - & (track.lat >= config.lat_min) & (track.lat <= config.lat_max) - ) - idx = np.flatnonzero(inside) - if idx.size < 2: - raise ValueError( - f"track {track.icao24!r} has fewer than 2 points inside the planning " - f"box {config.bbox}; cannot derive origin/destination" - ) - - i0, i1 = int(idx[0]), int(idx[-1]) - - # origin/destination: first/last in-box position -> local km. - ox, oy = anchor.geo_to_local(float(track.lon[i0]), float(track.lat[i0])) - dx, dy = anchor.geo_to_local(float(track.lon[i1]), float(track.lat[i1])) - - departure_s = max(0.0, float(track.time_s[i0]) - reference_epoch_s) - - # In-box slice, times relative to entry (t=0 at entry). - tt = track.time_s[i0:i1 + 1] - float(track.time_s[i0]) - alt = track.baro_altitude_m[i0:i1 + 1] - total_duration_s = float(tt[-1]) - - ac = aircraft_for(name or track.callsign or track.icao24, aircraft_type) - - profile = _observed_profile(tt, alt) - used_observed = profile is not None - if profile is None: - # Too sparse/noisy to trust the observed altitudes — fall back to the - # synthetic fuel-optimal climb. (build_dataset logs when this happens.) - duration = total_duration_s if total_duration_s > 0 else 3600.0 - profile = build_baseline_profile(ac, total_duration_s=duration) - - flight = Flight( - name=name or (track.callsign.strip() or track.icao24), - origin_km=(float(ox), float(oy)), - destination_km=(float(dx), float(dy)), - departure_s=departure_s, - aircraft=ac, - baseline=profile, - ) - return ReducedFlight(flight=flight, used_observed_baseline=used_observed) - - -def _observed_profile(tt: np.ndarray, alt_m: np.ndarray) -> AltitudeProfile | None: - """Resample observed barometric altitude into contiguous AltitudeSegments. - - Returns None when the data can't support a trustworthy profile (zero - duration, or every altitude sample missing), so the caller can fall back. - """ - if tt.size < 2 or tt[-1] <= tt[0]: - return None - finite = np.isfinite(alt_m) - if not finite.any(): - return None - - # Snap each finite altitude to the RVSM FL grid; carry the last good FL - # across any NaN gaps (forward-fill) so a few missing reports don't break - # the profile. - fls = np.empty(alt_m.shape, dtype=int) - last = None - for k in range(alt_m.size): - if finite[k]: - last = _snap_fl(float(alt_m[k])) - fls[k] = last if last is not None else _snap_fl(float(alt_m[finite][0])) - - # Group consecutive equal-FL runs into segments, boundaries at the time of - # each FL change. Segments are contiguous: each end == next start. - segments: list[AltitudeSegment] = [] - seg_start_t = float(tt[0]) - cur_fl = int(fls[0]) - for k in range(1, fls.size): - if int(fls[k]) != cur_fl: - segments.append(AltitudeSegment(cur_fl, seg_start_t, float(tt[k]))) - seg_start_t = float(tt[k]) - cur_fl = int(fls[k]) - segments.append(AltitudeSegment(cur_fl, seg_start_t, float(tt[-1]))) - - return AltitudeProfile(segments=tuple(segments)) - - -def _snap_fl(altitude_m: float) -> int: - """Snap an altitude (m) to the nearest in-envelope RVSM flight level.""" - fl = int(round(m_to_fl(altitude_m) / _FL_STEP) * _FL_STEP) - return max(_FL_MIN, min(_FL_MAX, fl)) diff --git a/contrail_flights/tracks.py b/contrail_flights/tracks.py deleted file mode 100644 index bba68e4..0000000 --- a/contrail_flights/tracks.py +++ /dev/null @@ -1,107 +0,0 @@ -""" -tracks.py — Raw ADS-B tracks and their cleaning. - -A `Track` is one aircraft's observed path: parallel arrays of (time, lat, lon, -baro_altitude) plus identity (icao24, callsign). The OpenSky client produces -these; `clean_tracks` drops the unusable ones (too few points, too short, or -flat noise) and normalizes each (sorted by time, de-duplicated timestamps). - -Pure numpy — no pandas, no network — so it imports in the lean install and the -hermetic tests can build fixtures directly. -""" - -from __future__ import annotations - -from dataclasses import dataclass - -import numpy as np - -# Mean Earth radius (km) for the haversine ground-distance estimate. -_EARTH_R_KM = 6371.0 - - -@dataclass(frozen=True) -class Track: - """One cleaned aircraft track. Arrays are parallel and time-ascending. - - time_s: UTC seconds (epoch). lat/lon: degrees. baro_altitude_m: metres - (may contain NaN where the report lacked altitude). - """ - - icao24: str - callsign: str - time_s: np.ndarray - lat: np.ndarray - lon: np.ndarray - baro_altitude_m: np.ndarray - - @property - def n_points(self) -> int: - return int(self.time_s.size) - - def ground_distance_km(self) -> float: - """Great-circle path length summed over consecutive points (km).""" - if self.n_points < 2: - return 0.0 - return float(np.sum(_haversine_km( - self.lat[:-1], self.lon[:-1], self.lat[1:], self.lon[1:] - ))) - - -def make_track( - icao24: str, - callsign: str, - time_s, - lat, - lon, - baro_altitude_m, -) -> Track: - """Build a Track from raw sequences, sorted by time with duplicate - timestamps removed (keeping the first). No filtering of short/sparse - tracks here — that is `clean_tracks`' job.""" - t = np.asarray(time_s, dtype=float) - order = np.argsort(t, kind="stable") - t = t[order] - lat_a = np.asarray(lat, dtype=float)[order] - lon_a = np.asarray(lon, dtype=float)[order] - alt_a = np.asarray(baro_altitude_m, dtype=float)[order] - - # Drop duplicate timestamps (ADS-B often repeats the last state vector). - keep = np.concatenate(([True], np.diff(t) > 0)) - return Track( - icao24=icao24, - callsign=callsign, - time_s=t[keep], - lat=lat_a[keep], - lon=lon_a[keep], - baro_altitude_m=alt_a[keep], - ) - - -def clean_tracks( - tracks: list[Track], - min_points: int = 10, - min_track_km: float = 200.0, -) -> list[Track]: - """Keep only tracks with enough points AND enough ground distance. - - A loud, explicit filter — never silently fabricate or pad a thin track. - """ - out: list[Track] = [] - for tr in tracks: - if tr.n_points < min_points: - continue - if tr.ground_distance_km() < min_track_km: - continue - out.append(tr) - return out - - -def _haversine_km(lat1, lon1, lat2, lon2) -> np.ndarray: - """Great-circle distance (km) between two arrays of (lat, lon) degrees.""" - p1 = np.radians(np.asarray(lat1, dtype=float)) - p2 = np.radians(np.asarray(lat2, dtype=float)) - dphi = p2 - p1 - dlam = np.radians(np.asarray(lon2, dtype=float) - np.asarray(lon1, dtype=float)) - a = np.sin(dphi / 2.0) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dlam / 2.0) ** 2 - return 2.0 * _EARTH_R_KM * np.arcsin(np.sqrt(np.clip(a, 0.0, 1.0))) diff --git a/contrail_ml/__init__.py b/contrail_ml/__init__.py deleted file mode 100644 index dfe7470..0000000 --- a/contrail_ml/__init__.py +++ /dev/null @@ -1,44 +0,0 @@ -""" -contrail_ml — Observation-grounded ISSR model + MLOps lifecycle. - -This package replaces the synthetic ISSR generator -(`contrail_env.synthetic_issr`) with a machine-learning model that predicts -ice-supersaturated regions (ISSRs) from numerical-weather-prediction (NWP) -fields, bias-corrected against in-situ aircraft humidity (IAGOS) and served -into the existing optimizer through the SAME duck-typed `ISSRField` interface. - -DESIGN INVARIANTS (see CONTRAIL_ML_IMPLEMENTATION_PLAN.md §0) -============================================================ -1. Never fabricate observational data. The only synthetic data is the - explicitly-named `synthetic_fallback` used SOLELY for tests/CI. -2. No training/serving skew: feature computation lives in ONE module - (`contrail_ml.features`), shared by the table builder and the predictor. -3. Don't claim to beat physics until measured: `evaluate` compares against - raw-ERA5 and quantile-mapping baselines before any "improvement". -4. Keep the deployed solver image lean: all ML deps live in the `[ml]` extra. -5. Don't break the existing system: the synthetic path stays the default. - -LIGHT IMPORT SURFACE -==================== -Importing `contrail_ml` pulls in only `features`/`config`, which depend on -numpy (+ optional pandas). The heavy estimators (xgboost, scikit-learn, -mlflow) and the geo stack (pycontrails, xarray) are imported lazily inside -the modules that need them, so `import contrail_ml` is cheap and the world -hook in `contrail_env` can import it without dragging in the full extra. -""" - -from __future__ import annotations - -from . import config, features - -__all__ = ["config", "features", "MLIssrField", "ml_issr_field"] - - -def __getattr__(name: str): - # Lazy re-export so `from contrail_ml import MLIssrField` works without - # forcing scipy/sklearn to import at package import time. - if name in ("MLIssrField", "ml_issr_field"): - from . import issr_field - - return getattr(issr_field, name) - raise AttributeError(f"module {__name__!r} has no attribute {name!r}") diff --git a/contrail_ml/__main__.py b/contrail_ml/__main__.py deleted file mode 100644 index d406d8c..0000000 --- a/contrail_ml/__main__.py +++ /dev/null @@ -1,10 +0,0 @@ -"""Entry point so `python -m contrail_ml ` works.""" - -from __future__ import annotations - -import sys - -from .cli import main - -if __name__ == "__main__": - sys.exit(main()) diff --git a/contrail_ml/baselines.py b/contrail_ml/baselines.py deleted file mode 100644 index 5e75283..0000000 --- a/contrail_ml/baselines.py +++ /dev/null @@ -1,135 +0,0 @@ -""" -baselines.py — The honest competition (plan §0.3, §6). - -Before claiming the ML model "beats physics", it has to beat the cheap things. -These three baselines all expose the same `predict_rhi` / `predict_proba` -surface as the model, so `evaluate.py` can score them side by side: - - RawERA5Baseline : use the NWP RHi as-is (the uncorrected forecast). - MultiplicativeFactorBaseline : one global factor f minimising ||f*rhi - truth||. - QuantileMappingBaseline : bivariate (T, RHi) empirical quantile mapping — the - standard statistical bias-correction the ML must beat. - -A baseline's classification score is its corrected RHi passed through a fixed -logistic centred at the ISSR threshold, so ROC/PR/Brier are all well defined. -""" - -from __future__ import annotations - -from dataclasses import dataclass, field - -import numpy as np -import pandas as pd - -# Logistic scale (in % RHi) used to turn a corrected-RHi baseline into a -# probability. Deliberately fixed (uncalibrated) — the point of the baselines -# is to be simple; calibration is the model's advantage. -_LOGISTIC_SCALE = 8.0 - - -def _logistic_from_rhi(rhi: np.ndarray, threshold: float) -> np.ndarray: - return 1.0 / (1.0 + np.exp(-(np.asarray(rhi, dtype=float) - threshold) / _LOGISTIC_SCALE)) - - -@dataclass -class RawERA5Baseline: - """The uncorrected NWP forecast: corrected RHi == raw RHi.""" - - issr_threshold: float = 100.0 - name: str = "raw_era5" - - def fit(self, df: pd.DataFrame) -> RawERA5Baseline: - return self - - def predict_rhi(self, df: pd.DataFrame) -> np.ndarray: - return df["rhi"].to_numpy(dtype=float) - - def predict_proba(self, df: pd.DataFrame) -> np.ndarray: - return _logistic_from_rhi(self.predict_rhi(df), self.issr_threshold) - - -@dataclass -class MultiplicativeFactorBaseline: - """Scale the NWP RHi by a single least-squares factor f.""" - - issr_threshold: float = 100.0 - name: str = "x_factor" - factor: float = 1.0 - - def fit(self, df: pd.DataFrame) -> MultiplicativeFactorBaseline: - rhi = df["rhi"].to_numpy(dtype=float) - truth = df["rhi_iagos"].to_numpy(dtype=float) - denom = float(np.sum(rhi * rhi)) - self.factor = float(np.sum(rhi * truth) / denom) if denom > 0 else 1.0 - return self - - def predict_rhi(self, df: pd.DataFrame) -> np.ndarray: - return self.factor * df["rhi"].to_numpy(dtype=float) - - def predict_proba(self, df: pd.DataFrame) -> np.ndarray: - return _logistic_from_rhi(self.predict_rhi(df), self.issr_threshold) - - -@dataclass -class QuantileMappingBaseline: - """Bivariate (T, RHi) empirical quantile mapping. - - Stratify by temperature bin; within each bin learn the monotone map from - the NWP-RHi quantiles to the IAGOS-RHi quantiles, then apply it at predict - time. This is the conventional statistical bias-correction the ML version - has to outperform to justify itself. - """ - - issr_threshold: float = 100.0 - n_temp_bins: int = 8 - n_quantiles: int = 50 - name: str = "quantile_map" - _t_edges: np.ndarray = field(default_factory=lambda: np.array([])) - _src_q: list = field(default_factory=list) - _dst_q: list = field(default_factory=list) - - def fit(self, df: pd.DataFrame) -> QuantileMappingBaseline: - t = df["T"].to_numpy(dtype=float) - rhi = df["rhi"].to_numpy(dtype=float) - truth = df["rhi_iagos"].to_numpy(dtype=float) - - self._t_edges = np.quantile(t, np.linspace(0, 1, self.n_temp_bins + 1)) - self._t_edges[0] -= 1e-6 - self._t_edges[-1] += 1e-6 - qs = np.linspace(0, 1, self.n_quantiles) - self._src_q, self._dst_q = [], [] - for lo, hi in zip(self._t_edges[:-1], self._t_edges[1:], strict=False): - m = (t >= lo) & (t < hi) - if m.sum() < 5: - # too few points: identity map for this bin - self._src_q.append(None) - self._dst_q.append(None) - continue - self._src_q.append(np.quantile(rhi[m], qs)) - self._dst_q.append(np.quantile(truth[m], qs)) - return self - - def predict_rhi(self, df: pd.DataFrame) -> np.ndarray: - t = df["T"].to_numpy(dtype=float) - rhi = df["rhi"].to_numpy(dtype=float) - out = rhi.copy() - bin_idx = np.clip(np.digitize(t, self._t_edges) - 1, 0, len(self._src_q) - 1) - for b in range(len(self._src_q)): - if self._src_q[b] is None: - continue - sel = bin_idx == b - if sel.any(): - out[sel] = np.interp(rhi[sel], self._src_q[b], self._dst_q[b]) - return out - - def predict_proba(self, df: pd.DataFrame) -> np.ndarray: - return _logistic_from_rhi(self.predict_rhi(df), self.issr_threshold) - - -def default_baselines(issr_threshold: float = 100.0) -> list: - """The standard baseline set used by the comparison table.""" - return [ - RawERA5Baseline(issr_threshold=issr_threshold), - MultiplicativeFactorBaseline(issr_threshold=issr_threshold), - QuantileMappingBaseline(issr_threshold=issr_threshold), - ] diff --git a/contrail_ml/calibrate.py b/contrail_ml/calibrate.py deleted file mode 100644 index 816fc7f..0000000 --- a/contrail_ml/calibrate.py +++ /dev/null @@ -1,138 +0,0 @@ -""" -calibrate.py — Make the probabilities mean what they say. - -A classifier that outputs 0.7 should be right ~70 % of the time at that score. -Raw tree-ensemble probabilities usually aren't; we fix that with a monotone -post-hoc map fitted on a held-out fold: - - isotonic : non-parametric, monotone step function (default; flexible). - platt : a 1-D logistic squashing (smoother, needs less data). - -We also provide a SPLIT-CONFORMAL interval on corrected RHi: a distribution- -free band rhi_hat ± q such that the true RHi lands inside ~(1 - alpha) of the -time, calibrated on held-out residuals. And `expected_calibration_error` (ECE) -to quantify how honest the probabilities are — the headline calibration metric. -""" - -from __future__ import annotations - -from dataclasses import dataclass -from typing import Any - -import numpy as np - - -@dataclass -class ProbabilityCalibrator: - """Post-hoc calibration map for P(ISSR). Fit on a held-out fold, then apply - to every future probability.""" - - method: str = "isotonic" - _model: Any = None - _fitted: bool = False - - def fit(self, p_raw: np.ndarray, y_true: np.ndarray) -> ProbabilityCalibrator: - p = np.asarray(p_raw, dtype=float).reshape(-1) - y = np.asarray(y_true, dtype=int).reshape(-1) - - if len(np.unique(y)) < 2: - # Can't calibrate against a single class; act as identity. - self._model = None - self._fitted = True - return self - - if self.method == "isotonic": - from sklearn.isotonic import IsotonicRegression - - iso = IsotonicRegression(out_of_bounds="clip", y_min=0.0, y_max=1.0) - iso.fit(p, y) - self._model = iso - elif self.method == "platt": - from sklearn.linear_model import LogisticRegression - - lr = LogisticRegression(C=1e6, solver="lbfgs") - lr.fit(p.reshape(-1, 1), y) - self._model = lr - else: - raise ValueError(f"unknown calibration method {self.method!r}") - self._fitted = True - return self - - def transform(self, p_raw: np.ndarray) -> np.ndarray: - if not self._fitted: - raise RuntimeError("ProbabilityCalibrator.transform before fit().") - p = np.asarray(p_raw, dtype=float).reshape(-1) - if self._model is None: - return np.clip(p, 0.0, 1.0) - if self.method == "isotonic": - return np.clip(self._model.predict(p), 0.0, 1.0) - return self._model.predict_proba(p.reshape(-1, 1))[:, 1] - - # convenience - def fit_transform(self, p_raw, y_true) -> np.ndarray: - return self.fit(p_raw, y_true).transform(p_raw) - - -def expected_calibration_error( - p: np.ndarray, y_true: np.ndarray, n_bins: int = 10 -) -> float: - """Expected Calibration Error: mean |confidence - accuracy| over equal-width - probability bins, weighted by bin population. Lower is better; the plan's - target is ECE < 0.05.""" - p = np.asarray(p, dtype=float).reshape(-1) - y = np.asarray(y_true, dtype=float).reshape(-1) - edges = np.linspace(0.0, 1.0, n_bins + 1) - ece = 0.0 - n = len(p) - for lo, hi in zip(edges[:-1], edges[1:], strict=False): - in_bin = (p >= lo) & (p < hi) if hi < 1.0 else (p >= lo) & (p <= hi) - if not in_bin.any(): - continue - conf = p[in_bin].mean() - acc = y[in_bin].mean() - ece += (in_bin.sum() / n) * abs(conf - acc) - return float(ece) - - -def reliability_curve( - p: np.ndarray, y_true: np.ndarray, n_bins: int = 10 -) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Return (bin_confidence, bin_accuracy, bin_count) for a reliability - diagram. Empty bins are dropped.""" - p = np.asarray(p, dtype=float).reshape(-1) - y = np.asarray(y_true, dtype=float).reshape(-1) - edges = np.linspace(0.0, 1.0, n_bins + 1) - conf, acc, cnt = [], [], [] - for lo, hi in zip(edges[:-1], edges[1:], strict=False): - in_bin = (p >= lo) & (p < hi) if hi < 1.0 else (p >= lo) & (p <= hi) - if not in_bin.any(): - continue - conf.append(p[in_bin].mean()) - acc.append(y[in_bin].mean()) - cnt.append(int(in_bin.sum())) - return np.array(conf), np.array(acc), np.array(cnt) - - -@dataclass -class ConformalRHi: - """Split-conformal interval half-width for corrected RHi. - - Calibrate on held-out absolute residuals |rhi_true - rhi_hat|; the (1-alpha) - empirical quantile is a distribution-free half-width q such that the true - value lies within rhi_hat ± q about (1-alpha) of the time. - """ - - alpha: float = 0.1 - half_width: float = float("nan") - - def fit(self, rhi_true: np.ndarray, rhi_hat: np.ndarray) -> ConformalRHi: - resid = np.abs(np.asarray(rhi_true, dtype=float) - np.asarray(rhi_hat, dtype=float)) - n = len(resid) - # finite-sample-adjusted quantile level - level = min(1.0, np.ceil((n + 1) * (1 - self.alpha)) / n) - self.half_width = float(np.quantile(resid, level)) - return self - - def interval(self, rhi_hat: np.ndarray) -> tuple[np.ndarray, np.ndarray]: - r = np.asarray(rhi_hat, dtype=float) - return r - self.half_width, r + self.half_width diff --git a/contrail_ml/cli.py b/contrail_ml/cli.py deleted file mode 100644 index 98f6e80..0000000 --- a/contrail_ml/cli.py +++ /dev/null @@ -1,264 +0,0 @@ -""" -cli.py — `python -m contrail_ml ` entry points (plan §9). - -Commands -======== - build-dataset : IAGOS + ERA5 -> versioned parquet (+ dataset card). - --synthetic writes a fallback table for offline work. - train : CV, fit, calibrate, evaluate, log to MLflow, register Staging. - evaluate : recompute the model-vs-baselines table for a dataset. - predict : build an ISSR field for a region/time and save it (NPZ). - serve-check : build an MLIssrField, drop it into the world via issr_source - ="ml", run the existing CP-SAT solve — a smoke test of the - whole serving seam. - monitor : rolling skill + feature drift over a window; print a retrain - recommendation. - -`--synthetic` keeps train/predict/serve-check hermetic (no network/credentials), -behind the same loud guards as the rest of the package. -""" - -from __future__ import annotations - -import argparse -import json -import sys -import warnings -from typing import Any - -import numpy as np -import pandas as pd - -from .config import MLConfig - - -def _cfg(args: argparse.Namespace) -> MLConfig: - cfg = MLConfig.from_env() - if getattr(args, "mlflow_uri", None): - cfg = cfg.with_overrides(mlflow_tracking_uri=args.mlflow_uri) - return cfg - - -# --------------------------------------------------------------------------- # -# build-dataset -# --------------------------------------------------------------------------- # -def cmd_build_dataset(args: argparse.Namespace) -> int: - cfg = _cfg(args) - if args.synthetic: - from .data.synthetic_fallback import make_synthetic_training_table - - df = make_synthetic_training_table(n=args.n, seed=args.seed, allow_synthetic=True) - import os - - os.makedirs(f"{cfg.data_dir}/processed", exist_ok=True) - path = f"{cfg.data_dir}/processed/synthetic_{args.seed}.parquet" - df.to_parquet(path, index=False) - print(f"[synthetic] wrote {len(df)} rows -> {path}") - return 0 - - from .data.build_dataset import build_and_write - - parquet, card = build_and_write(cfg, sampling_weight=args.sampling_weight) - print(f"wrote {parquet}\ncard {card}") - return 0 - - -# --------------------------------------------------------------------------- # -# train -# --------------------------------------------------------------------------- # -def cmd_train(args: argparse.Namespace) -> int: - cfg = _cfg(args) - df = _load_dataset(cfg, args) - from .train import run_training - - res = run_training(cfg, df, log_mlflow=not args.no_mlflow, do_cv=not args.no_cv) - print("dataset_hash:", res.dataset_hash) - if res.cv_metrics: - print("cv:", {k: round(v, 4) for k, v in res.cv_metrics.items()}) - print("holdout:", {k: round(v, 4) for k, v in res.holdout_metrics.items()}) - print("\nmodel vs baselines:") - print(res.comparison.round(4).to_string()) - if res.mlflow: - print("\nmlflow:", res.mlflow) - return 0 - - -# --------------------------------------------------------------------------- # -# evaluate -# --------------------------------------------------------------------------- # -def cmd_evaluate(args: argparse.Namespace) -> int: - cfg = _cfg(args) - df = _load_dataset(cfg, args) - from .baselines import default_baselines - from .evaluate import comparison_table - from .registry import load_model - - model = load_model(cfg, stage=args.stage) - baselines = [b.fit(df) for b in default_baselines(cfg.issr_rhi_threshold)] - table = comparison_table(df, model, baselines, p_threshold=cfg.issr_p_threshold) - print(table.round(4).to_string()) - return 0 - - -# --------------------------------------------------------------------------- # -# predict -# --------------------------------------------------------------------------- # -def cmd_predict(args: argparse.Namespace) -> int: - cfg = _cfg(args) - from .issr_field import ml_issr_field - - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - field = ml_issr_field( - cfg, - met_source="synthetic" if args.synthetic else args.met_source, - time=args.time, - allow_synthetic=args.synthetic, - grid_res_deg=args.grid_res, - ) - out = args.out or "issr_field.npz" - np.savez( - out, - lon=field.lon_axis, lat=field.lat_axis, pressure=field.pressure_axis, - rhi_excess=field.rhi_excess_cube, p_issr=field.p_issr_cube, - ) - frac = float((field.p_issr_cube >= field.p_threshold).mean()) - print(f"wrote {out} cube={field.rhi_excess_cube.shape} " - f"ISSR fraction={frac:.3f}") - return 0 - - -# --------------------------------------------------------------------------- # -# serve-check -# --------------------------------------------------------------------------- # -def cmd_serve_check(args: argparse.Namespace) -> int: - cfg = _cfg(args) - from contrail_env import ( - build_and_evaluate_flight, - build_capacity_buckets, - build_conflict_graph, - build_random_flights, - default_european_world, - solve_cpsat, - ) - - issr_kwargs: dict[str, Any] = dict( - config=cfg, met_source="synthetic" if args.synthetic else args.met_source, - allow_synthetic=args.synthetic, grid_res_deg=args.grid_res, - ) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - world = default_european_world(seed=args.seed, issr_source="ml", - issr_kwargs=issr_kwargs) - flights = build_random_flights(n_flights=args.n_flights, world=world, - seed=args.seed, corridor_frac=0.05, - snapshot_window_s=(0.0, 300.0)) - evals = [] - for f in flights: - evals.extend(build_and_evaluate_flight(f, world)) - conflicts = build_conflict_graph(evals, world) - buckets = build_capacity_buckets(evals, world) - result = solve_cpsat(evals, conflicts, buckets, time_limit_s=10.0) - - print(f"ISSR source : ml ({issr_kwargs['met_source']})") - print(f"status : {result.status}") - print(f"objective : {result.objective:.1f}") - print(f"conflicts : {len(conflicts)}") - print("chosen options:") - for i in result.chosen_eval_indices: - ev = evals[i] - print(f" {ev.flight_name}: option {ev.option_index} " - f"fuel={ev.fuel_kg:.0f}kg contrail={ev.contrail_cells} " - f"disrupt={ev.disruption_FLmin:.1f}") - return 0 - - -# --------------------------------------------------------------------------- # -# monitor -# --------------------------------------------------------------------------- # -def cmd_monitor(args: argparse.Namespace) -> int: - cfg = _cfg(args) - from .monitor import monitor - from .registry import load_model - - reference = pd.read_parquet(args.reference) - current = pd.read_parquet(args.current) - model = load_model(cfg, stage=args.stage) - report = monitor(model, reference, current, cfg) - print(json.dumps(report.to_dict(), indent=2, default=float)) - return 0 - - -# --------------------------------------------------------------------------- # -def _load_dataset(cfg: MLConfig, args: argparse.Namespace) -> pd.DataFrame: - if getattr(args, "synthetic", False): - from .data.synthetic_fallback import make_synthetic_training_table - - return make_synthetic_training_table(n=args.n, seed=args.seed, allow_synthetic=True) - if not getattr(args, "data", None): - raise SystemExit("provide --data or --synthetic") - return pd.read_parquet(args.data) - - -def build_parser() -> argparse.ArgumentParser: - p = argparse.ArgumentParser(prog="python -m contrail_ml") - p.add_argument("--mlflow-uri", default=None, help="override MLflow tracking URI") - sub = p.add_subparsers(dest="command", required=True) - - bd = sub.add_parser("build-dataset", help="build the training table") - bd.add_argument("--synthetic", action="store_true") - bd.add_argument("--n", type=int, default=5000) - bd.add_argument("--seed", type=int, default=0) - bd.add_argument("--sampling-weight", action="store_true") - bd.set_defaults(func=cmd_build_dataset) - - tr = sub.add_parser("train", help="train + calibrate + register") - tr.add_argument("--data", default=None, help="training parquet path") - tr.add_argument("--synthetic", action="store_true") - tr.add_argument("--n", type=int, default=8000) - tr.add_argument("--seed", type=int, default=0) - tr.add_argument("--no-mlflow", action="store_true") - tr.add_argument("--no-cv", action="store_true") - tr.set_defaults(func=cmd_train) - - ev = sub.add_parser("evaluate", help="model-vs-baselines table") - ev.add_argument("--data", default=None) - ev.add_argument("--synthetic", action="store_true") - ev.add_argument("--n", type=int, default=8000) - ev.add_argument("--seed", type=int, default=0) - ev.add_argument("--stage", default="Staging") - ev.set_defaults(func=cmd_evaluate) - - pr = sub.add_parser("predict", help="build + save an ISSR field") - pr.add_argument("--synthetic", action="store_true") - pr.add_argument("--met-source", default="gfs", choices=["gfs", "era5"]) - pr.add_argument("--time", default=None) - pr.add_argument("--grid-res", type=float, default=1.0) - pr.add_argument("--out", default=None) - pr.set_defaults(func=cmd_predict) - - sc = sub.add_parser("serve-check", help="solve with the ML ISSR field") - sc.add_argument("--synthetic", action="store_true", default=True) - sc.add_argument("--real", dest="synthetic", action="store_false") - sc.add_argument("--met-source", default="gfs", choices=["gfs", "era5"]) - sc.add_argument("--grid-res", type=float, default=2.0) - sc.add_argument("--seed", type=int, default=1) - sc.add_argument("--n-flights", type=int, default=4) - sc.set_defaults(func=cmd_serve_check) - - mo = sub.add_parser("monitor", help="rolling skill + drift") - mo.add_argument("--reference", required=True, help="training reference parquet") - mo.add_argument("--current", required=True, help="new labelled window parquet") - mo.add_argument("--stage", default="Staging") - mo.set_defaults(func=cmd_monitor) - - return p - - -def main(argv: list[str] | None = None) -> int: - args = build_parser().parse_args(argv) - return int(args.func(args)) - - -if __name__ == "__main__": - sys.exit(main()) diff --git a/contrail_ml/config.py b/contrail_ml/config.py deleted file mode 100644 index 373eb1e..0000000 --- a/contrail_ml/config.py +++ /dev/null @@ -1,145 +0,0 @@ -""" -config.py — Typed configuration for the contrail_ml pipeline. - -One dataclass (`MLConfig`) carries every knob: the geographic region and -pressure levels to pull, the temporal train/test split, model hyperparameters, -the MLflow tracking URI, and the geographic anchor that maps the local sim -frame to real (lon, lat). It loads from environment variables (12-factor) or a -YAML file, and hardcodes NO secrets — credentials/paths come from the -environment. - -Why a single config object? - Both the training-table builder and the live predictor read the SAME - region/levels/anchor from here, which (together with the shared - `features.py`) is what guarantees train/serve consistency. -""" - -from __future__ import annotations - -import os -from dataclasses import dataclass, field, replace -from typing import Any - -from .features import GeoAnchor - -# The local sim domain, copied from contrail_env.synthetic_issr.ISSRField.domain -# (x_km, y_km, z_m). Kept here so the ML grid and the optimizer agree on extent. -SIM_DOMAIN: tuple[float, float, float, float, float, float] = ( - 0.0, 1500.0, 0.0, 800.0, 9000.0, 12500.0 -) - - -@dataclass(frozen=True) -class MLConfig: - """All configuration for building, training, serving and monitoring.""" - - # ---- Geographic region (degrees) and pressure levels (hPa) ---------- - # Europe / North-Atlantic box covering the cruise region. - lon_min: float = -30.0 - lon_max: float = 30.0 - lat_min: float = 30.0 - lat_max: float = 70.0 - pressure_levels_hpa: tuple[int, ...] = (150, 175, 200, 225, 250, 300) - - # ---- Temporal split (inclusive year bounds) ------------------------- - train_years: tuple[int, int] = (2011, 2019) - test_years: tuple[int, int] = (2020, 2022) - - # ---- Geographic anchor: local (x=0, y=0) -> (lat, lon), + ref time --- - # Anchors the 1500x800 km sim box over south-west -> central Europe. - origin_lat: float = 43.0 - origin_lon: float = -5.0 - ref_time: str = "2019-07-01T12:00:00" # ISO; used for serving snapshots - - # ---- Model hyperparameters ------------------------------------------ - regime_split_rhi: float = 85.0 # % RHi boundary between dry/humid sub-models - n_ensemble: int = 10 # bootstrap re-fits for predictive spread - issr_rhi_threshold: float = 100.0 # RHi % at/above which a cell is ISSR - issr_p_threshold: float = 0.5 # calibrated P(ISSR) cut for is_inside() - xgb_max_depth: int = 6 - xgb_n_estimators: int = 300 - xgb_learning_rate: float = 0.05 - mlp_hidden: tuple[int, ...] = (64, 32) - random_state: int = 42 - - # ---- Paths / MLflow -------------------------------------------------- - data_dir: str = "data" - mlflow_tracking_uri: str = "file:./mlruns" - registered_model_name: str = "contrail-issr-rhi-corrector" - - # ---- Credentials / data paths (read from env, never hardcoded) ------ - iagos_dir: str | None = None # local dir of IAGOS NetCDF files - - # Free-form extra knobs without widening the schema. - extra: dict[str, Any] = field(default_factory=dict) - - # ------------------------------------------------------------------ # - - @property - def anchor(self) -> GeoAnchor: - """The geographic anchor used by BOTH the model grid and the GUI.""" - return GeoAnchor(origin_lat=self.origin_lat, origin_lon=self.origin_lon) - - @property - def bbox(self) -> tuple[float, float, float, float]: - """(lon_min, lon_max, lat_min, lat_max) for the weather loaders.""" - return (self.lon_min, self.lon_max, self.lat_min, self.lat_max) - - def with_overrides(self, **kwargs: Any) -> MLConfig: - """Return a copy with the given fields replaced.""" - return replace(self, **kwargs) - - # ------------------------------------------------------------------ # - # Loaders - # ------------------------------------------------------------------ # - - @classmethod - def from_env(cls, prefix: str = "CONTRAIL_ML_") -> MLConfig: - """Build a config, overriding any field from `${PREFIX}FIELD` env vars. - - Example: CONTRAIL_ML_MLFLOW_TRACKING_URI, CONTRAIL_ML_IAGOS_DIR. - Only fields present in the environment are overridden; types are - coerced from the dataclass defaults. - """ - base = cls() - overrides: dict[str, Any] = {} - for f in base.__dataclass_fields__.values(): # type: ignore[attr-defined] - if f.name == "extra": - continue - env_key = f"{prefix}{f.name.upper()}" - if env_key in os.environ: - overrides[f.name] = _coerce(getattr(base, f.name), os.environ[env_key]) - return replace(base, **overrides) - - @classmethod - def from_yaml(cls, path: str) -> MLConfig: - """Load a config from a YAML file (only known fields are applied).""" - import yaml # lazy: PyYAML is a transitive dep, not always present - - with open(path, encoding="utf-8") as fh: - data = yaml.safe_load(fh) or {} - known = {f.name for f in cls.__dataclass_fields__.values()} # type: ignore[attr-defined] - clean = {k: v for k, v in data.items() if k in known} - return cls(**clean) - - -def _coerce(reference: Any, raw: str) -> Any: - """Coerce an env-var string to the type of the dataclass default value.""" - if isinstance(reference, bool): - return raw.strip().lower() in ("1", "true", "yes", "on") - if isinstance(reference, int) and not isinstance(reference, bool): - return int(raw) - if isinstance(reference, float): - return float(raw) - if isinstance(reference, tuple): - parts = [p.strip() for p in raw.split(",") if p.strip()] - # Preserve element type from the first reference element if available. - if reference and isinstance(reference[0], int): - return tuple(int(p) for p in parts) - if reference and isinstance(reference[0], float): - return tuple(float(p) for p in parts) - return tuple(parts) - return raw - - -DEFAULT_CONFIG = MLConfig() diff --git a/contrail_ml/data/__init__.py b/contrail_ml/data/__init__.py deleted file mode 100644 index 620977f..0000000 --- a/contrail_ml/data/__init__.py +++ /dev/null @@ -1,6 +0,0 @@ -"""contrail_ml.data — data loaders, collocation, and the training-table builder. - -Real loaders (iagos/era5/gfs) need credentials/network and are NOT exercised in -CI. The hermetic test fixture is `synthetic_fallback`, which is loudly guarded -so it can never be mistaken for real observations. -""" diff --git a/contrail_ml/data/build_dataset.py b/contrail_ml/data/build_dataset.py deleted file mode 100644 index f80cc89..0000000 --- a/contrail_ml/data/build_dataset.py +++ /dev/null @@ -1,119 +0,0 @@ -""" -build_dataset.py — Orchestrate IAGOS + ERA5 -> one versioned training table. - -Pipeline: load IAGOS waypoints (truth) -> collocate ERA5 onto them (inputs) -> -compute the shared features -> attach labels (rhi_iagos, y_issr, residual delta) --> enforce the schema -> write a content-addressed parquet plus a dataset_card -describing coverage and class balance. If DVC is available the parquet is -`dvc add`-ed; otherwise the content hash is the version. - -The training table is the single source of truth every model run records (by -hash), so any reported metric is reproducible from a known dataset. -""" - -from __future__ import annotations - -import hashlib -import json -import os -from typing import TYPE_CHECKING - -import pandas as pd - -from ..features import add_derived_features -from .schema import ( - LABEL_DELTA, - LABEL_ISSR, - LABEL_RHI, - TABLE_COLUMNS, - issr_label, - validate_table, -) - -if TYPE_CHECKING: - from ..config import MLConfig - - -def build_training_table(cfg: MLConfig, *, sampling_weight: bool = False) -> pd.DataFrame: - """Build the labelled training table from real IAGOS + ERA5 data.""" - from .collocate import add_sampling_weight, collocate_era5 - from .iagos import load_iagos_waypoints - - iagos = load_iagos_waypoints(cfg) - collocated = collocate_era5(iagos, cfg) - df = add_derived_features(collocated, pressure_col="pressure_hpa") - - # Labels: IAGOS measured RHi is the truth; the residual is what the model - # learns; y_issr is the binary supervised target. - df[LABEL_RHI] = collocated["rhi_iagos"].to_numpy() - df[LABEL_DELTA] = df[LABEL_RHI].to_numpy() - df["rhi"].to_numpy() - df[LABEL_ISSR] = issr_label(df[LABEL_RHI].to_numpy(), cfg.issr_rhi_threshold) - - df["qc_flag"] = collocated.get("qc_flag", 0) - df["source"] = "iagos+era5" - df = add_sampling_weight(df, enable=sampling_weight) - - table = df.loc[:, list(TABLE_COLUMNS)].reset_index(drop=True) - validate_table(table) - return table - - -def content_hash(df: pd.DataFrame) -> str: - h = hashlib.sha256(pd.util.hash_pandas_object(df, index=False).values.tobytes()) - return h.hexdigest()[:16] - - -def write_versioned(df: pd.DataFrame, cfg: MLConfig) -> tuple[str, str]: - """Write the table to a content-addressed parquet + a dataset_card.json. - - Returns (parquet_path, card_path). - """ - out_dir = os.path.join(cfg.data_dir, "processed") - os.makedirs(out_dir, exist_ok=True) - digest = content_hash(df) - parquet_path = os.path.join(out_dir, f"issr_training_{digest}.parquet") - df.to_parquet(parquet_path, index=False) - - card = _dataset_card(df, cfg, digest) - card_path = os.path.join(out_dir, f"issr_training_{digest}.card.json") - with open(card_path, "w", encoding="utf-8") as fh: - json.dump(card, fh, indent=2, default=str) - - _maybe_dvc_add(parquet_path) - return parquet_path, card_path - - -def _dataset_card(df: pd.DataFrame, cfg: MLConfig, digest: str) -> dict: - t = pd.to_datetime(df["time"]) - return { - "hash": digest, - "n_rows": int(len(df)), - "time_start": str(t.min()), - "time_end": str(t.max()), - "lon_range": [float(df["lon"].min()), float(df["lon"].max())], - "lat_range": [float(df["lat"].min()), float(df["lat"].max())], - "pressure_levels_hpa": list(cfg.pressure_levels_hpa), - "issr_rate": float(df[LABEL_ISSR].mean()), - "rhi_bias_raw": float((df["rhi"] - df[LABEL_RHI]).mean()), - "source": "iagos+era5", - "schema": list(TABLE_COLUMNS), - } - - -def _maybe_dvc_add(path: str) -> None: - """Best-effort `dvc add` for data versioning; silent if DVC isn't set up.""" - import shutil - import subprocess - - if shutil.which("dvc") is None: - return - try: - subprocess.run(["dvc", "add", path], check=False, - stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL) - except Exception: - pass - - -def build_and_write(cfg: MLConfig, *, sampling_weight: bool = False) -> tuple[str, str]: - df = build_training_table(cfg, sampling_weight=sampling_weight) - return write_versioned(df, cfg) diff --git a/contrail_ml/data/collocate.py b/contrail_ml/data/collocate.py deleted file mode 100644 index bd310d0..0000000 --- a/contrail_ml/data/collocate.py +++ /dev/null @@ -1,78 +0,0 @@ -""" -collocate.py — Put ERA5 onto the IAGOS waypoints (plan §4). - -For each IAGOS measurement we need the NWP fields AT THAT point in space-time, so -the model input row lines up with the truth label. This is a 4-D (time, level, -lat, lon) linear interpolation of the ERA5 dataset onto the waypoint coordinates -via xarray `.interp`. - -Known caveat (documented, not hidden): IAGOS aircraft avoid deep convection, so -the sample is biased away from the most intense humidity. We surface this as an -optional `sample_weight` column and call it out in docs/ML.md rather than -pretending the sample is unbiased. -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING - -import numpy as np -import pandas as pd - -from ..features import RAW_NWP_COLUMNS - -if TYPE_CHECKING: - from ..config import MLConfig - - -def collocate_era5(iagos_df: pd.DataFrame, cfg: MLConfig) -> pd.DataFrame: - """Interpolate ERA5 raw NWP fields onto IAGOS waypoints. - - Returns the IAGOS frame with the RAW_NWP_COLUMNS added (the model inputs), - so the caller can compute features and attach labels. Heavy imports are - deferred; fails loudly via era5.open_era5 if the geo stack is missing. - """ - import xarray as xr - - from .era5 import ERA5_VARIABLES, open_era5 - - if iagos_df.empty: - raise ValueError("collocate_era5 got an empty IAGOS frame") - - times = pd.to_datetime(iagos_df["time"]) - t0, t1 = times.min(), times.max() - met = open_era5(cfg, (np.datetime64(t0), np.datetime64(t1))) - ds = getattr(met, "data", met) - - pts = { - "time": xr.DataArray(times.values, dims="points"), - "level": xr.DataArray(iagos_df["pressure_hpa"].to_numpy(), dims="points"), - "latitude": xr.DataArray(iagos_df["lat"].to_numpy(), dims="points"), - "longitude": xr.DataArray(iagos_df["lon"].to_numpy(), dims="points"), - } - interp = ds.interp(**pts) - - out = iagos_df.copy().reset_index(drop=True) - for name in RAW_NWP_COLUMNS: - cf = ERA5_VARIABLES[name] - out[name] = np.asarray(interp[cf].values, dtype=float) - - # The IAGOS T/q are the truth; the ERA5 T/q are the model inputs. Keep the - # ERA5 versions under the canonical feature names, and the measured RHi as - # the label (already on iagos_df as rhi_iagos). - return out.dropna(subset=list(RAW_NWP_COLUMNS)).reset_index(drop=True) - - -def add_sampling_weight(df: pd.DataFrame, enable: bool = False) -> pd.DataFrame: - """Optionally add a `sample_weight` correcting (crudely) for IAGOS avoiding - the most humid/convective air. Off by default; documented in ML.md.""" - out = df.copy() - if not enable: - out["sample_weight"] = 1.0 - return out - # Up-weight the under-sampled humid tail so the loss isn't dominated by the - # over-represented dry regime. - rhi = out["rhi"].to_numpy() if "rhi" in out.columns else out["rhi_iagos"].to_numpy() - w = 1.0 + np.clip(rhi - 90.0, 0.0, None) / 30.0 - out["sample_weight"] = w - return out diff --git a/contrail_ml/data/era5.py b/contrail_ml/data/era5.py deleted file mode 100644 index ccd5387..0000000 --- a/contrail_ml/data/era5.py +++ /dev/null @@ -1,119 +0,0 @@ -""" -era5.py — Pull ERA5 reanalysis over the region/levels/time (training input). - -ERA5 is the best reconstruction of PAST weather (reanalysis): it assimilates -real observations, so it is the right source for the model INPUTS at training -time. We use the ARCO-ERA5 zarr mirror on a public Google Cloud bucket via -pycontrails, which needs NO Copernicus CDS credentials. - -NOT hermetic: this needs the `[ml]` geo stack (pycontrails/xarray/zarr/gcsfs) -and network. Every entry point fails loudly with an actionable message if a -dependency or the data is missing — it never fabricates fields (plan §0.1). The -heavy imports are deferred into the functions so importing this module is cheap. - -VERIFY AGAINST DOCS before a real run: the pycontrails ARCO-ERA5 loader class -and the exact CF variable names occasionally change between releases. -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING, Any - -import numpy as np - -from ..features import RAW_NWP_COLUMNS - -if TYPE_CHECKING: - from ..config import MLConfig - from ..predict import MetCube - -# Map our raw feature names -> ERA5/CF standard names used by pycontrails. -# (Vertical velocity in ERA5 is omega = Lagrangian tendency of pressure, Pa/s.) -ERA5_VARIABLES: dict[str, str] = { - "T": "air_temperature", - "specific_humidity": "specific_humidity", - "u": "eastward_wind", - "v": "northward_wind", - "w": "lagrangian_tendency_of_air_pressure", - "cloud_cover": "fraction_of_cloud_cover", - "cloud_ice": "specific_cloud_ice_water_content", -} - - -def _require_pycontrails() -> Any: - try: - import pycontrails # noqa: F401 - except ImportError as exc: # pragma: no cover - exercised only with extras - raise ImportError( - "ERA5 access needs the [ml] extra (pycontrails, xarray, zarr, gcsfs). " - "Install with: pip install -e \".[ml]\". See docs/DATA.md." - ) from exc - from pycontrails.datalib.ecmwf import arco_era5 - - return arco_era5 - - -def open_era5( - cfg: MLConfig, - time: Any, - *, - variables: tuple[str, ...] = RAW_NWP_COLUMNS, -): - """Open an ARCO-ERA5 MetDataset for `time` over the configured region/levels. - - `time` is a single timestamp or a (start, end) pair. Returns a pycontrails - MetDataset (xarray-backed). Variable selection is restricted to the - ERA5∩GFS intersection (RAW_NWP_COLUMNS) so nothing is missing at serve time. - """ - arco_era5 = _require_pycontrails() - cf_vars = [ERA5_VARIABLES[v] for v in variables] - # NOTE: verify the exact constructor/kwargs against the installed - # pycontrails version before a real pull. - met = arco_era5.open_arco_era5( - time=time, - variables=cf_vars, - pressure_levels=list(cfg.pressure_levels_hpa), - ) - return _subset_region(met, cfg) - - -def _subset_region(met: Any, cfg: MLConfig) -> Any: - """Crop a MetDataset/xarray to the configured lon/lat box.""" - ds = getattr(met, "data", met) - lon_min, lon_max, lat_min, lat_max = cfg.bbox - return ds.sel(longitude=slice(lon_min, lon_max), latitude=slice(lat_max, lat_min)) - - -def load_met_cube( - cfg: MLConfig, - time: str | None = None, - *, - grid_res_deg: float = 1.0, -) -> MetCube: - """Load a single-time ERA5 snapshot as a serving MetCube (retrospective - replay alternative to GFS).""" - - t = time or cfg.ref_time - ds = _subset_region(open_era5(cfg, t), cfg) - return _dataset_to_metcube(ds, cfg, grid_res_deg, t) - - -def _dataset_to_metcube(ds: Any, cfg: MLConfig, grid_res_deg: float, t: Any) -> MetCube: - """Regrid an xarray dataset onto an ascending (lon, lat, pressure) grid and - pack it into a MetCube.""" - from ..predict import MetCube - - lon_axis = np.arange(cfg.lon_min, cfg.lon_max + 1e-6, grid_res_deg) - lat_axis = np.arange(cfg.lat_min, cfg.lat_max + 1e-6, grid_res_deg) - pressure_axis = np.array(sorted(cfg.pressure_levels_hpa), dtype=float) - - interp = ds.interp(longitude=lon_axis, latitude=lat_axis, level=pressure_axis) - fields: dict[str, np.ndarray] = {} - for name in RAW_NWP_COLUMNS: - cf = ERA5_VARIABLES[name] - da = interp[cf] - # squeeze any singleton time dim and order axes (lon, lat, level) - arr = np.asarray(da.transpose("longitude", "latitude", "level").values, dtype=float) - fields[name] = arr - return MetCube(lon_axis=lon_axis, lat_axis=lat_axis, pressure_axis=pressure_axis, - fields=fields, time=t) diff --git a/contrail_ml/data/gfs.py b/contrail_ml/data/gfs.py deleted file mode 100644 index c3f387c..0000000 --- a/contrail_ml/data/gfs.py +++ /dev/null @@ -1,90 +0,0 @@ -""" -gfs.py — Pull a GFS forecast for the serving features (plan §7). - -At serve time you need to know where ISSR WILL be when the aircraft arrives, not -where it is now — so the operational input is a FORECAST. GFS is NOAA's global -forecast model: free, public-domain, refreshed every 6 hours, delivered as GRIB -(hence cfgrib in the [ml] extra). pycontrails ships a `GFSForecast` loader. - -Same honesty contract as era5.py: lazy heavy imports, loud actionable failure -if deps/data are missing, never fabricate. Features are restricted to the -ERA5∩GFS intersection so the model never sees a missing column. - -VERIFY AGAINST DOCS: GFS variable names / the GFSForecast signature before a -real pull. -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING, Any - -from ..features import RAW_NWP_COLUMNS - -if TYPE_CHECKING: - from ..config import MLConfig - from ..predict import MetCube - -# GFS exposes the same physical fields under its own names. w is the GFS -# vertical velocity (Pa/s), matching ERA5's omega so the feature is consistent. -GFS_VARIABLES: dict[str, str] = { - "T": "air_temperature", - "specific_humidity": "specific_humidity", - "u": "eastward_wind", - "v": "northward_wind", - "w": "lagrangian_tendency_of_air_pressure", - "cloud_cover": "fraction_of_cloud_cover", - "cloud_ice": "specific_cloud_ice_water_content", -} - - -def _require_gfs() -> Any: - try: - from pycontrails.datalib.gfs import GFSForecast - except ImportError as exc: # pragma: no cover - needs extras - raise ImportError( - "GFS access needs the [ml] extra (pycontrails, cfgrib). Install with: " - "pip install -e \".[ml]\". See docs/DATA.md." - ) from exc - return GFSForecast - - -def open_gfs(cfg: MLConfig, time: Any, *, variables: tuple[str, ...] = RAW_NWP_COLUMNS): - """Open a GFSForecast MetDataset for `time` over the configured region/levels.""" - GFSForecast = _require_gfs() - cf_vars = [GFS_VARIABLES[v] for v in variables] - # NOTE: verify constructor kwargs against the installed pycontrails version. - gfs = GFSForecast( - time=time, - variables=cf_vars, - pressure_levels=list(cfg.pressure_levels_hpa), - ) - met = gfs.open_metdataset() - return getattr(met, "data", met) - - -def load_forecast_cube( - cfg: MLConfig, - time: str | None = None, - *, - grid_res_deg: float = 1.0, -) -> MetCube: - """Load a single forecast valid-time as a serving MetCube.""" - from .era5 import _dataset_to_metcube # shared regrid/pack logic - - t = time or cfg.ref_time - ds = open_gfs(cfg, t) - ds = ds.sel( - longitude=slice(cfg.lon_min, cfg.lon_max), - latitude=slice(cfg.lat_max, cfg.lat_min), - ) - # GFS CF names match the ERA5 mapping values here, so the shared packer works. - cube = _dataset_to_metcube(ds, cfg, grid_res_deg, t) - return cube - - -def _gfs_names_match_era5() -> bool: - """Guard used by the packer: the two mappings must share CF names so the - shared regrid code can read either dataset.""" - from .era5 import ERA5_VARIABLES - - return all(GFS_VARIABLES[k] == ERA5_VARIABLES[k] for k in RAW_NWP_COLUMNS) diff --git a/contrail_ml/data/iagos.py b/contrail_ml/data/iagos.py deleted file mode 100644 index 86f5fb4..0000000 --- a/contrail_ml/data/iagos.py +++ /dev/null @@ -1,113 +0,0 @@ -""" -iagos.py — Load IAGOS in-situ aircraft humidity as ground truth (plan §4). - -IAGOS (In-service Aircraft for a Global Observing System) equips commercial -aircraft with research-grade humidity sensors. Its measured RHi near the -tropopause is the closest thing to truth for ISSR — and the label this whole -project trains against. IAGOS distributes one NetCDF per flight; access needs a -(free) registration on the IAGOS portal (see docs/DATA.md). The data path comes -from config/env, never hardcoded, and no credentials live in the repo. - -We RECOMPUTE RHi from the measured T, humidity and pressure with the SAME -thermodynamics as the features (contrail_ml.features), rather than trusting a -pre-computed RHi column, so the label and the features share one definition (no -label skew). The function fails loudly if the data directory is absent — it -never invents measurements. - -VERIFY AGAINST DOCS: IAGOS variable names differ between products (IAGOS-CORE, -MOZAIC, CARIBIC). The mapping below is a documented default; override via config. -""" - -from __future__ import annotations - -import glob -import os -from typing import TYPE_CHECKING - -import numpy as np -import pandas as pd - -from ..features import EPSILON, rhi_percent - -if TYPE_CHECKING: - from ..config import MLConfig - -# Default IAGOS NetCDF variable names -> our canonical columns. Override per -# product through cfg.extra["iagos_variables"]. -IAGOS_VARIABLES: dict[str, str] = { - "time": "UTC_time", - "lat": "lat", - "lon": "lon", - "pressure_pa": "air_press_AC", # static pressure, Pa - "temperature_k": "air_temp_AC", # air temperature, K - "h2o_vmr_ppmv": "H2O_gas_P", # water-vapour volume mixing ratio, ppmv - "qc": "H2O_gas_P_val", # validation/quality flag (product-specific) -} - -# Validation flags that we keep (product-specific; documented in DATA.md). -GOOD_QC_VALUES: tuple[int, ...] = (0, 1) - - -def load_iagos_waypoints(cfg: MLConfig) -> pd.DataFrame: - """Load every IAGOS flight NetCDF under cfg.iagos_dir into one waypoint - dataframe with columns: lat, lon, pressure_hpa, time, T, specific_humidity, - rhi_iagos, qc_flag. - - Raises a clear error if the directory is missing/empty — never fabricates. - """ - data_dir = cfg.iagos_dir or os.environ.get("CONTRAIL_ML_IAGOS_DIR") - if not data_dir: - raise FileNotFoundError( - "IAGOS directory not set. Set cfg.iagos_dir or CONTRAIL_ML_IAGOS_DIR " - "to a folder of IAGOS NetCDF files. See docs/DATA.md for registration." - ) - if not os.path.isdir(data_dir): - raise FileNotFoundError(f"IAGOS directory does not exist: {data_dir}") - files = sorted(glob.glob(os.path.join(data_dir, "*.nc"))) - if not files: - raise FileNotFoundError(f"no IAGOS *.nc files found in {data_dir}") - - var = {**IAGOS_VARIABLES, **(cfg.extra.get("iagos_variables", {}))} - frames = [_load_one(f, var, cfg) for f in files] - df = pd.concat(frames, ignore_index=True) - return df.dropna(subset=["T", "specific_humidity", "pressure_hpa"]).reset_index(drop=True) - - -def _load_one(path: str, var: dict[str, str], cfg: MLConfig) -> pd.DataFrame: - import xarray as xr # lazy: part of the [ml] extra - - ds = xr.open_dataset(path) - try: - pressure_pa = np.asarray(ds[var["pressure_pa"]].values, dtype=float) - temperature = np.asarray(ds[var["temperature_k"]].values, dtype=float) - vmr_ppmv = np.asarray(ds[var["h2o_vmr_ppmv"]].values, dtype=float) - lat = np.asarray(ds[var["lat"]].values, dtype=float) - lon = np.asarray(ds[var["lon"]].values, dtype=float) - time = pd.to_datetime(np.asarray(ds[var["time"]].values)) - qc = (np.asarray(ds[var["qc"]].values) if var["qc"] in ds else - np.zeros_like(temperature, dtype=int)) - finally: - ds.close() - - pressure_hpa = pressure_pa / 100.0 - # ppmv volume mixing ratio -> specific humidity (kg/kg). - vmr = vmr_ppmv * 1e-6 - mmr = EPSILON * vmr # mass mixing ratio - q = mmr / (1.0 + mmr) - - rhi = np.asarray(rhi_percent(temperature, q, pressure_hpa), dtype=float) - - df = pd.DataFrame( - { - "lat": lat, "lon": lon, "pressure_hpa": pressure_hpa, "time": time, - "T": temperature, "specific_humidity": q, "rhi_iagos": rhi, - "qc_flag": np.asarray(qc).astype(int), - } - ) - # Quality screening + region/level window. - df = df[df["qc_flag"].isin(GOOD_QC_VALUES)] - lon_min, lon_max, lat_min, lat_max = cfg.bbox - df = df[(df["lon"].between(lon_min, lon_max)) & (df["lat"].between(lat_min, lat_max))] - p_lo, p_hi = min(cfg.pressure_levels_hpa), max(cfg.pressure_levels_hpa) - df = df[df["pressure_hpa"].between(p_lo - 25, p_hi + 25)] - return df.reset_index(drop=True) diff --git a/contrail_ml/data/schema.py b/contrail_ml/data/schema.py deleted file mode 100644 index 267b1d1..0000000 --- a/contrail_ml/data/schema.py +++ /dev/null @@ -1,60 +0,0 @@ -""" -schema.py — Column contract for the training table. - -One row = one IAGOS waypoint, with ERA5 fields collocated onto it. Keeping the -column names, dtypes and groupings in one place lets every stage -(build_dataset, model, evaluate, monitor) agree without guessing. - -Groups -====== -keys : where/when the row is (for grouped CV and provenance). -features : the model inputs (defined in features.FEATURES — single source). -labels : the supervised targets derived from the IAGOS measurement. -meta : quality flags, source tag, optional sampling weight. -""" - -from __future__ import annotations - -from ..features import FEATURES - -# ---- key columns ---------------------------------------------------------- -KEY_COLUMNS: tuple[str, ...] = ("lat", "lon", "pressure_hpa", "time") - -# ---- label columns -------------------------------------------------------- -# rhi_iagos : measured RHi (%) — the ground truth, recomputed from IAGOS T/q/p -# with the SAME thermodynamics as the features (no label skew). -# y_issr : (rhi_iagos >= 100) as int — the classification target. -# delta : rhi_iagos - rhi(nwp) — the residual the corrector predicts. -LABEL_RHI = "rhi_iagos" -LABEL_ISSR = "y_issr" -LABEL_DELTA = "delta" -LABEL_COLUMNS: tuple[str, ...] = (LABEL_RHI, LABEL_ISSR, LABEL_DELTA) - -# ---- meta columns --------------------------------------------------------- -META_COLUMNS: tuple[str, ...] = ("qc_flag", "source", "sample_weight") - -# Full ordered table schema. -TABLE_COLUMNS: tuple[str, ...] = ( - *KEY_COLUMNS, - *FEATURES, - *LABEL_COLUMNS, - *META_COLUMNS, -) - - -def issr_label(rhi_iagos, rhi_threshold: float = 100.0): - """y_issr = (measured RHi >= threshold). Works on scalars or arrays.""" - import numpy as np - - return (np.asarray(rhi_iagos, dtype=float) >= rhi_threshold).astype(int) - - -def validate_table(df) -> None: - """Raise if the dataframe is missing any required column. A loud schema - check beats a model that silently trains on the wrong thing.""" - missing = [c for c in TABLE_COLUMNS if c not in df.columns] - if missing: - raise ValueError( - f"training table is missing required columns {missing}; " - f"expected schema = {list(TABLE_COLUMNS)}" - ) diff --git a/contrail_ml/data/synthetic_fallback.py b/contrail_ml/data/synthetic_fallback.py deleted file mode 100644 index d9df6b3..0000000 --- a/contrail_ml/data/synthetic_fallback.py +++ /dev/null @@ -1,158 +0,0 @@ -""" -synthetic_fallback.py — Hermetic, clearly-labelled fake training data. - -PURPOSE AND THE HARD RULE -========================= -The whole project's credibility rests on never passing synthetic data off as -real observations (plan §0.1). This module is the *single* exception: a tiny -labelled table generated from a KNOWN function, used SOLELY by unit tests / CI -so the model, calibration, and integration code can be exercised with no -network and no credentials. - -Two guards enforce that it can't leak into real work: - 1. It refuses to run unless the caller passes `allow_synthetic=True`. - 2. It emits a loud warning every time it is used, and stamps every row with - source == "synthetic_fallback". - -WHAT IT MODELS -============== -A physically-flavoured bias-correction problem. Each row has a "true" -(IAGOS-like) relative humidity over ice and a DRY-BIASED NWP humidity — the -real, documented failure mode of weather models near the tropopause. The bias -is small in dry air and grows in the humid regime (RHi > ~85%), and depends on -cloud-ice / vertical-velocity too, so: - * the residual `delta = rhi_iagos - rhi(nwp)` is genuinely learnable, - * the regime split at 85% RHi separates two different bias behaviours, - * noise keeps P(ISSR) genuinely uncertain near the threshold (so calibration - has something to do). -""" - -from __future__ import annotations - -import warnings - -import numpy as np -import pandas as pd - -from ..features import ( - EPSILON, - add_derived_features, - e_si_pa, -) -from .schema import LABEL_DELTA, LABEL_ISSR, LABEL_RHI, TABLE_COLUMNS, issr_label - -SOURCE_TAG = "synthetic_fallback" - - -def make_synthetic_training_table( - n: int = 4000, - seed: int = 0, - *, - allow_synthetic: bool = False, - rhi_threshold: float = 100.0, -) -> pd.DataFrame: - """Generate a labelled training table that matches the real schema. - - Parameters - ---------- - n, seed : size and reproducibility. - allow_synthetic : MUST be True. The guard exists so production code can - never call this by accident. - rhi_threshold : RHi % used to derive the binary ISSR label. - - Returns a DataFrame with exactly `schema.TABLE_COLUMNS`. - """ - if not allow_synthetic: - raise RuntimeError( - "synthetic_fallback refused: pass allow_synthetic=True. This data " - "is for tests/CI only and must never stand in for real IAGOS/ERA5." - ) - warnings.warn( - "USING SYNTHETIC FALLBACK DATA — not real observations. " - "Valid for tests/CI only; never report metrics from this as science.", - stacklevel=2, - ) - - rng = np.random.default_rng(seed) - - # ---- keys & geometry -------------------------------------------------- - pressure_hpa = rng.choice([150, 175, 200, 225, 250, 300], size=n).astype(float) - lon = rng.uniform(-30.0, 30.0, size=n) - lat = rng.uniform(30.0, 70.0, size=n) - # Spread timestamps across several years and times of day. - base = np.datetime64("2015-01-01T00:00:00") - offsets = rng.integers(0, 6 * 365 * 24 * 60, size=n) # minutes over ~6 yr - time = base + offsets.astype("timedelta64[m]") - - # ---- temperature: colder higher up, with weather noise ---------------- - # ~225 K near 225 hPa, falling ~0.05 K/hPa toward lower pressure. - T = 225.0 - 0.05 * (pressure_hpa - 225.0) + rng.normal(0.0, 4.0, size=n) - T = np.clip(T, 200.0, 245.0) - - # ---- pick a TRUE RHi spanning dry..supersaturated --------------------- - # Mixture: mostly sub-saturated, a meaningful humid/ISSR tail. - true_rhi = rng.uniform(20.0, 150.0, size=n) - - # Auxiliary fields, partly correlated with humidity (physical flavour). - cloud_ice = np.clip( - rng.normal(0.0, 1.0, size=n) * 1e-6 + 8e-6 * np.clip(true_rhi - 90.0, 0.0, None) / 60.0, - 0.0, None, - ) - w = rng.normal(0.0, 0.1, size=n) - 0.05 * np.clip(true_rhi - 100.0, 0.0, None) / 50.0 - cloud_cover = np.clip(0.15 + 0.6 * np.clip(true_rhi - 80.0, 0.0, None) / 70.0 - + rng.normal(0.0, 0.1, size=n), 0.0, 1.0) - u = rng.normal(10.0, 12.0, size=n) - v = rng.normal(0.0, 12.0, size=n) - - # ---- the DRY BIAS: model underestimates humidity, worse when humid ---- - humid_excess = np.clip(true_rhi - 85.0, 0.0, None) - bias = ( - 2.0 # small baseline dry bias - + 0.28 * humid_excess # grows in the humid regime - + 1.5e6 * cloud_ice # ice-laden air biased more - - 12.0 * np.clip(-w, 0.0, None) # strong updraughts biased more - + rng.normal(0.0, 4.0, size=n) # irreducible noise - ) - nwp_rhi = np.clip(true_rhi - bias, 1.0, 200.0) - - # Back out the NWP specific humidity consistent with nwp_rhi (so that the - # `rhi` feature recomputed by features.py reproduces nwp_rhi exactly). - q_nwp = _q_from_rhi(nwp_rhi, T, pressure_hpa) - - raw = pd.DataFrame( - { - "lat": lat, - "lon": lon, - "pressure_hpa": pressure_hpa, - "time": time, - "T": T, - "specific_humidity": q_nwp, - "u": u, - "v": v, - "w": w, - "cloud_cover": cloud_cover, - "cloud_ice": cloud_ice, - } - ) - - df = add_derived_features(raw, pressure_col="pressure_hpa") - - # ---- labels (truth) + residual + meta -------------------------------- - df[LABEL_RHI] = true_rhi - df[LABEL_DELTA] = true_rhi - df["rhi"].to_numpy() # = bias the model learns - df[LABEL_ISSR] = issr_label(true_rhi, rhi_threshold) - df["qc_flag"] = "ok" - df["source"] = SOURCE_TAG - df["sample_weight"] = 1.0 - - return df.loc[:, list(TABLE_COLUMNS)].reset_index(drop=True) - - -def _q_from_rhi(rhi_pct: np.ndarray, temperature_k: np.ndarray, - pressure_hpa: np.ndarray) -> np.ndarray: - """Invert RHi -> specific humidity q (kg/kg), the inverse of - features.vapor_pressure_pa composed with rhi_percent.""" - esi = np.asarray(e_si_pa(temperature_k), dtype=float) - e = (rhi_pct / 100.0) * esi # vapour pressure, Pa - p_pa = pressure_hpa * 100.0 - return EPSILON * e / (p_pa - (1.0 - EPSILON) * e) diff --git a/contrail_ml/evaluate.py b/contrail_ml/evaluate.py deleted file mode 100644 index be1dc06..0000000 --- a/contrail_ml/evaluate.py +++ /dev/null @@ -1,174 +0,0 @@ -""" -evaluate.py — Metrics, baseline comparison, and diagnostic plots. - -This is where the project earns the right to claim anything (plan §0.3, §6). -Nothing here trains; it scores predictions against the IAGOS truth on the -held-out set and tabulates the model next to the raw-ERA5 and quantile-mapping -baselines. The headline numbers: - - * RHi dry-bias reduction (mean error vs IAGOS) and RMSE/MAE, - * ISSR skill: ETS, F1, ROC-AUC, PR-AUC, - * probability honesty: Brier score and Expected Calibration Error. - -Report whatever the numbers are — modest gains, reported rigorously with -calibrated probabilities, are themselves the contribution. -""" - -from __future__ import annotations - -from typing import TYPE_CHECKING, Any - -import numpy as np -import pandas as pd - -from .calibrate import expected_calibration_error, reliability_curve - -if TYPE_CHECKING: - from matplotlib.figure import Figure - - -def equitable_threat_score(y_true: np.ndarray, y_pred: np.ndarray) -> float: - """ETS (Gilbert skill score) for a binary forecast. Chance-corrected hit - rate; 0 = no skill, 1 = perfect. The plan's bar is beating raw ERA5's - ~0.2-0.4.""" - yt = np.asarray(y_true, dtype=int) - yp = np.asarray(y_pred, dtype=int) - hits = int(np.sum((yp == 1) & (yt == 1))) - false_alarms = int(np.sum((yp == 1) & (yt == 0))) - misses = int(np.sum((yp == 0) & (yt == 1))) - n = len(yt) - pred_pos = hits + false_alarms - obs_pos = hits + misses - hits_random = pred_pos * obs_pos / n if n else 0.0 - denom = hits + false_alarms + misses - hits_random - return float((hits - hits_random) / denom) if denom != 0 else 0.0 - - -def score_predictions( - rhi_true: np.ndarray, - y_true: np.ndarray, - rhi_hat: np.ndarray, - p_issr: np.ndarray, - *, - p_threshold: float = 0.5, -) -> dict[str, float]: - """All metrics for one predictor, as a flat dict.""" - rhi_true = np.asarray(rhi_true, dtype=float) - y_true = np.asarray(y_true, dtype=int) - rhi_hat = np.asarray(rhi_hat, dtype=float) - p_issr = np.asarray(p_issr, dtype=float) - y_pred = (p_issr >= p_threshold).astype(int) - - err = rhi_hat - rhi_true - out: dict[str, float] = { - "rhi_rmse": float(np.sqrt(np.mean(err ** 2))), - "rhi_mae": float(np.mean(np.abs(err))), - "rhi_bias": float(np.mean(err)), - "ets": equitable_threat_score(y_true, y_pred), - } - - # sklearn metrics; guard the single-class case (AUC undefined). - from sklearn.metrics import ( - average_precision_score, - brier_score_loss, - f1_score, - precision_score, - recall_score, - roc_auc_score, - ) - - out["f1"] = float(f1_score(y_true, y_pred, zero_division=0)) - out["precision"] = float(precision_score(y_true, y_pred, zero_division=0)) - out["recall"] = float(recall_score(y_true, y_pred, zero_division=0)) - if len(np.unique(y_true)) == 2: - out["roc_auc"] = float(roc_auc_score(y_true, p_issr)) - out["pr_auc"] = float(average_precision_score(y_true, p_issr)) - out["brier"] = float(brier_score_loss(y_true, np.clip(p_issr, 0, 1))) - else: - out["roc_auc"] = float("nan") - out["pr_auc"] = float("nan") - out["brier"] = float("nan") - out["ece"] = expected_calibration_error(np.clip(p_issr, 0, 1), y_true) - return out - - -def comparison_table( - holdout: pd.DataFrame, - model: Any, - baselines: list, - *, - p_threshold: float = 0.5, -) -> pd.DataFrame: - """One row per predictor (model + each baseline), one column per metric. - - `model` is an RHiCorrector-like object (.predict -> rhi_hat, rhi_std, - p_issr); each baseline exposes predict_rhi/predict_proba. The baselines must - already be fitted on the training split. - """ - rhi_true = holdout["rhi_iagos"].to_numpy(dtype=float) - y_true = holdout["y_issr"].to_numpy(dtype=int) - - rows: dict[str, dict[str, float]] = {} - - rhi_hat, _std, p_issr = model.predict(holdout) - rows["ml_corrector"] = score_predictions( - rhi_true, y_true, rhi_hat, p_issr, p_threshold=p_threshold - ) - for b in baselines: - rows[b.name] = score_predictions( - rhi_true, y_true, b.predict_rhi(holdout), b.predict_proba(holdout), - p_threshold=p_threshold, - ) - - table = pd.DataFrame(rows).T - table.index.name = "predictor" - return table - - -# =========================================================================== -# Diagnostic plots (logged to MLflow as artifacts) -# =========================================================================== - - -def plot_reliability(p_issr: np.ndarray, y_true: np.ndarray) -> Figure: - """Reliability diagram (calibration curve) with the y=x ideal.""" - import matplotlib - - matplotlib.use("Agg") - import matplotlib.pyplot as plt - - conf, acc, _cnt = reliability_curve(np.clip(p_issr, 0, 1), y_true) - fig, ax = plt.subplots(figsize=(4.5, 4.5)) - ax.plot([0, 1], [0, 1], "k--", lw=1, label="perfect") - ax.plot(conf, acc, "o-", label="model") - ax.set_xlabel("predicted P(ISSR)") - ax.set_ylabel("observed ISSR frequency") - ax.set_title("Reliability") - ax.legend() - fig.tight_layout() - return fig - - -def plot_rhi_scatter(rhi_true: np.ndarray, rhi_hat: np.ndarray, - rhi_raw: np.ndarray) -> Figure: - """Corrected vs raw RHi against IAGOS truth — visualises the bias removal.""" - import matplotlib - - matplotlib.use("Agg") - import matplotlib.pyplot as plt - - fig, ax = plt.subplots(figsize=(5, 5)) - lim = (0.0, max(160.0, float(np.max(rhi_true)) + 5)) - ax.plot(lim, lim, "k--", lw=1) - ax.scatter(rhi_true, rhi_raw, s=4, alpha=0.3, label="raw NWP") - ax.scatter(rhi_true, rhi_hat, s=4, alpha=0.3, label="corrected") - ax.axvline(100, color="grey", lw=0.7) - ax.axhline(100, color="grey", lw=0.7) - ax.set_xlabel("IAGOS RHi (%)") - ax.set_ylabel("predicted RHi (%)") - ax.set_title("RHi: raw vs corrected") - ax.set_xlim(lim) - ax.set_ylim(lim) - ax.legend() - fig.tight_layout() - return fig diff --git a/contrail_ml/features.py b/contrail_ml/features.py deleted file mode 100644 index 3f20113..0000000 --- a/contrail_ml/features.py +++ /dev/null @@ -1,312 +0,0 @@ -""" -features.py — SHARED feature engineering (training AND serving). - -This is the anti-skew guarantee (plan §0.2): every feature the model ever sees -is computed here, once, by code used identically by the training-table builder -(`data/build_dataset.py`) and the live predictor (`predict.py`). If a feature -were computed one way at train time and another at serve time, the model would -silently degrade in production. Centralizing it makes that class of bug -impossible. - -Contents -======== -1. Thermodynamics — the physics that turns (T, q, p) into relative humidity - over ice (RHi). Murphy & Koop (2005) ice-saturation vapour pressure. -2. Standard-atmosphere altitude <-> pressure (ISA barometric formula). -3. Cyclical time encodings (solar time, day-of-year) as sin/cos pairs. -4. The geographic anchor mapping the local sim frame <-> real (lon, lat). -5. `FEATURES` — the exact, ordered list of model-input columns, restricted - to variables present in BOTH ERA5 and GFS so nothing is missing at serve - time. - -All thermodynamic functions accept floats or numpy arrays and return the same -shape, so they vectorize over a whole grid or a whole training table. -""" - -from __future__ import annotations - -from dataclasses import dataclass -from typing import TYPE_CHECKING - -import numpy as np - -if TYPE_CHECKING: # pandas is in the [ml] extra; keep import-time deps light. - import pandas as pd - -# =========================================================================== -# 1. THERMODYNAMICS -# =========================================================================== - -# Dry-air / water-vapour molecular-mass ratio (Rd/Rv). -EPSILON = 0.621981 - - -def e_si_pa(temperature_k: np.ndarray | float) -> np.ndarray | float: - """Saturation vapour pressure OVER ICE, in Pa. Murphy & Koop (2005). - - e_si = exp(9.550426 - 5723.265/T + 3.53068*ln(T) - 0.00728332*T) - - Valid for T in roughly [110, 273.16] K, which covers the entire upper - troposphere where ISSRs live. T in kelvin. - """ - t = np.asarray(temperature_k, dtype=float) - out = np.exp( - 9.550426 - - 5723.265 / t - + 3.53068 * np.log(t) - - 0.00728332 * t - ) - return _maybe_scalar(out, temperature_k) - - -def vapor_pressure_pa( - specific_humidity: np.ndarray | float, - pressure_pa: np.ndarray | float, -) -> np.ndarray | float: - """Water-vapour partial pressure (Pa) from specific humidity q and total - pressure p: - - e = q * p / (EPSILON + (1 - EPSILON) * q) - - q is dimensionless (kg water / kg moist air); p and e share units (Pa). - """ - q = np.asarray(specific_humidity, dtype=float) - p = np.asarray(pressure_pa, dtype=float) - e = q * p / (EPSILON + (1.0 - EPSILON) * q) - return _maybe_scalar(e, specific_humidity) - - -def rhi_percent( - temperature_k: np.ndarray | float, - specific_humidity: np.ndarray | float, - pressure_hpa: np.ndarray | float, -) -> np.ndarray | float: - """Relative humidity over ice, in PERCENT, from (T[K], q, p[hPa]). - - RHi = 100 * e / e_si(T) - - RHi >= 100 % is the supersaturation that lets a contrail persist. - """ - p_pa = np.asarray(pressure_hpa, dtype=float) * 100.0 - e = np.asarray(vapor_pressure_pa(specific_humidity, p_pa), dtype=float) - esi = np.asarray(e_si_pa(temperature_k), dtype=float) - rhi = 100.0 * e / esi - return _maybe_scalar(rhi, temperature_k) - - -def rhi_excess(rhi_pct: np.ndarray | float) -> np.ndarray | float: - """ISSRField semantics: max(0, RHi/100 - 1). - - 0 outside supersaturation; grows with how far above ice-saturation the - air is. This is exactly the quantity the synthetic ISSR blobs return, so - the ML field is a drop-in. - """ - r = np.asarray(rhi_pct, dtype=float) - out = np.clip(r / 100.0 - 1.0, 0.0, None) - return _maybe_scalar(out, rhi_pct) - - -# =========================================================================== -# 2. STANDARD ATMOSPHERE — altitude <-> pressure -# =========================================================================== - -_T0 = 288.15 # sea-level standard temperature (K) -_L = 0.0065 # tropospheric lapse rate (K/m) -_P0 = 1013.25 # sea-level standard pressure (hPa) -_G = 9.80665 # gravity (m/s^2) -_R = 287.05 # specific gas constant for dry air (J/kg/K) -_H_TROP = 11_000.0 # tropopause altitude (m) -_EXP = _G / (_R * _L) # ~5.2559 -_P_TROP = _P0 * (1.0 - _L * _H_TROP / _T0) ** _EXP # ~226.32 hPa -_T_TROP = _T0 - _L * _H_TROP # 216.65 K - - -def altitude_to_pressure_hpa(altitude_m: np.ndarray | float) -> np.ndarray | float: - """ISA pressure (hPa) at a geometric altitude (m). Piecewise: lapse-rate - troposphere below 11 km, isothermal stratosphere above.""" - h = np.asarray(altitude_m, dtype=float) - trop = _P0 * np.clip(1.0 - _L * h / _T0, 1e-9, None) ** _EXP - strat = _P_TROP * np.exp(-_G * (h - _H_TROP) / (_R * _T_TROP)) - out = np.where(h <= _H_TROP, trop, strat) - return _maybe_scalar(out, altitude_m) - - -def pressure_to_altitude_m(pressure_hpa: np.ndarray | float) -> np.ndarray | float: - """Inverse of altitude_to_pressure_hpa (ISA).""" - p = np.asarray(pressure_hpa, dtype=float) - trop = (_T0 / _L) * (1.0 - (p / _P0) ** (1.0 / _EXP)) - strat = _H_TROP - (_R * _T_TROP / _G) * np.log(np.clip(p / _P_TROP, 1e-9, None)) - out = np.where(p >= _P_TROP, trop, strat) - return _maybe_scalar(out, pressure_hpa) - - -# =========================================================================== -# 3. CYCLICAL TIME ENCODINGS -# =========================================================================== -# A model must not learn "11pm and midnight are 23 hours apart". Encoding a -# periodic quantity as a (sin, cos) pair puts the two ends of the cycle next -# to each other in feature space. - - -def solar_time_encoding( - utc_hour: np.ndarray | float, - lon_deg: np.ndarray | float, -) -> tuple[np.ndarray | float, np.ndarray | float]: - """(sin, cos) of LOCAL solar time. Local solar hour ~= UTC + lon/15.""" - h = np.asarray(utc_hour, dtype=float) - lon = np.asarray(lon_deg, dtype=float) - local = np.mod(h + lon / 15.0, 24.0) - ang = 2.0 * np.pi * local / 24.0 - return _maybe_scalar(np.sin(ang), utc_hour), _maybe_scalar(np.cos(ang), utc_hour) - - -def day_of_year_encoding( - day_of_year: np.ndarray | float, -) -> tuple[np.ndarray | float, np.ndarray | float]: - """(sin, cos) of the day-of-year (seasonal cycle).""" - d = np.asarray(day_of_year, dtype=float) - ang = 2.0 * np.pi * d / 365.25 - return _maybe_scalar(np.sin(ang), day_of_year), _maybe_scalar(np.cos(ang), day_of_year) - - -# =========================================================================== -# 4. GEOGRAPHIC ANCHOR — local sim frame <-> real (lon, lat) -# =========================================================================== - - -@dataclass(frozen=True) -class GeoAnchor: - """Maps the local sim frame (x_km east, y_km north) to geography. - - The sim is a flat 1500x800 km Cartesian box; real weather lives on - (lon, lat). One small-angle anchor ties them together so the model, the - predicted field, and the GUI map all agree on where things are. - - lat = origin_lat + y_km / KM_PER_DEG_LAT - lon = origin_lon + x_km / (KM_PER_DEG_LAT * cos(lat)) - """ - - origin_lat: float - origin_lon: float - km_per_deg_lat: float = 111.0 - - def local_to_geo( - self, x_km: np.ndarray | float, y_km: np.ndarray | float - ) -> tuple[np.ndarray | float, np.ndarray | float]: - """(x_km, y_km) -> (lon, lat) in degrees.""" - x = np.asarray(x_km, dtype=float) - y = np.asarray(y_km, dtype=float) - lat = self.origin_lat + y / self.km_per_deg_lat - lon = self.origin_lon + x / (self.km_per_deg_lat * np.cos(np.radians(lat))) - return _maybe_scalar(lon, x_km), _maybe_scalar(lat, y_km) - - def geo_to_local( - self, lon_deg: np.ndarray | float, lat_deg: np.ndarray | float - ) -> tuple[np.ndarray | float, np.ndarray | float]: - """(lon, lat) -> (x_km, y_km). Inverse of local_to_geo.""" - lon = np.asarray(lon_deg, dtype=float) - lat = np.asarray(lat_deg, dtype=float) - y = (lat - self.origin_lat) * self.km_per_deg_lat - x = (lon - self.origin_lon) * self.km_per_deg_lat * np.cos(np.radians(lat)) - return _maybe_scalar(x, lon_deg), _maybe_scalar(y, lat_deg) - - -# =========================================================================== -# 5. THE FEATURE LIST + the shared table builder -# =========================================================================== - -# Raw NWP columns expected on an input row (present in BOTH ERA5 and GFS). -RAW_NWP_COLUMNS: tuple[str, ...] = ( - "T", # temperature, K - "specific_humidity", # kg/kg - "u", "v", "w", # wind components (m/s; w = vertical velocity Pa/s) - "cloud_cover", # fraction [0,1] - "cloud_ice", # specific cloud ice water content, kg/kg -) - -# The exact, ordered model-input columns. `rhi` is the NWP relative humidity -# over ice (the quantity the model corrects) — computed from THIS row's T/q/p, -# whether that row came from ERA5 (train) or GFS (serve). One name, no skew. -FEATURES: tuple[str, ...] = ( - "T", - "specific_humidity", - "rhi", - "u", - "v", - "w", - "cloud_cover", - "cloud_ice", - "altitude_m", - "sin_solar", - "cos_solar", - "sin_doy", - "cos_doy", -) - - -def add_derived_features(df: pd.DataFrame, pressure_col: str = "pressure_hpa") -> pd.DataFrame: - """Add every derived feature column to `df` IN PLACE-SAFE fashion. - - Requires raw columns: T, specific_humidity, , u, v, w, - cloud_cover, cloud_ice, lon, and a `time` column (datetime-like or a - numeric UTC-hour fallback). Returns a NEW frame with the FEATURES columns - present, so the same call site serves both the training table and the - live grid. - """ - out = df.copy() - - out["rhi"] = rhi_percent(out["T"].to_numpy(), out["specific_humidity"].to_numpy(), - out[pressure_col].to_numpy()) - out["rhi_excess"] = rhi_excess(out["rhi"].to_numpy()) - out["altitude_m"] = pressure_to_altitude_m(out[pressure_col].to_numpy()) - - utc_hour, doy = _time_parts(out) - sin_solar, cos_solar = solar_time_encoding(utc_hour, out["lon"].to_numpy()) - sin_doy, cos_doy = day_of_year_encoding(doy) - out["sin_solar"] = sin_solar - out["cos_solar"] = cos_solar - out["sin_doy"] = sin_doy - out["cos_doy"] = cos_doy - return out - - -def feature_matrix(df: pd.DataFrame) -> pd.DataFrame: - """Select exactly the FEATURES columns, in order. Raises if any is absent - (a loud failure beats silently training on the wrong columns).""" - missing = [c for c in FEATURES if c not in df.columns] - if missing: - raise KeyError( - f"feature_matrix: missing columns {missing}; " - f"did you call add_derived_features() first?" - ) - return df.loc[:, list(FEATURES)] - - -# =========================================================================== -# Internal helpers -# =========================================================================== - - -def _maybe_scalar(arr: np.ndarray, like: np.ndarray | float) -> np.ndarray | float: - """Return a Python float when the input was scalar, else the array.""" - if np.isscalar(like) or (isinstance(like, np.ndarray) and like.ndim == 0): - return float(arr) - return arr - - -def _time_parts(df: pd.DataFrame) -> tuple[np.ndarray, np.ndarray]: - """Extract (utc_hour, day_of_year) arrays from a `time` column. - - Accepts a real datetime column or, as a fallback for synthetic tests, - explicit `utc_hour`/`day_of_year` numeric columns. - """ - import pandas as pd - - if "time" in df.columns: - t = pd.to_datetime(df["time"]) - utc_hour = (t.dt.hour + t.dt.minute / 60.0).to_numpy() - doy = t.dt.dayofyear.to_numpy().astype(float) - return utc_hour, doy - if "utc_hour" in df.columns and "day_of_year" in df.columns: - return df["utc_hour"].to_numpy(), df["day_of_year"].to_numpy() - raise KeyError("add_derived_features needs a 'time' column (or utc_hour + day_of_year).") diff --git a/contrail_ml/issr_field.py b/contrail_ml/issr_field.py deleted file mode 100644 index 6a54e08..0000000 --- a/contrail_ml/issr_field.py +++ /dev/null @@ -1,181 +0,0 @@ -""" -issr_field.py — MLIssrField: the model, dressed up as an ISSRField. - -THE SEAM (plan §1) -================== -Everything in the optimizer that touches ISSR goes through the duck-typed -`ISSRField` interface from `contrail_env.synthetic_issr`: - - rhi_excess(x_km, y_km, z_m) -> float - is_inside(x_km, y_km, z_m) -> bool - rhi_excess_grid(x, y, z: ndarray) -> ndarray - mask_grid(x, y, z: ndarray) -> ndarray - threshold: float - -`World.is_issr_cell` and `qubo.build_conflict_graph` consume ONLY that -interface. So the trained model integrates by implementing the same contract — -no change to `World` or `qubo.py`. This file is that implementation. - -The model lives on a real geographic grid `(lon, lat, pressure)`. When the -optimizer asks about a local sim point `(x_km, y_km, z_m)`, we map it through -the shared `GeoAnchor` (x/y -> lon/lat) and the standard atmosphere -(z_m -> pressure), then interpolate the predicted field. The local<->geo -transform is the SAME one the features and the GUI use, so everyone agrees on -where things are. -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any - -import numpy as np - -from .features import GeoAnchor, altitude_to_pressure_hpa - -if TYPE_CHECKING: - from .config import MLConfig - -# Local sim box (x_km, y_km, z_m), matching ISSRField.domain. -_SIM_DOMAIN = (0.0, 1500.0, 0.0, 800.0, 9000.0, 12500.0) - - -@dataclass -class MLIssrField: - """A predicted ISSR-risk field exposing the `ISSRField` interface. - - Holds the corrected-RHi-excess and calibrated P(ISSR) on an ascending - `(lon, lat, pressure)` grid and interpolates them on query. - - Parameters - ---------- - lon_axis, lat_axis, pressure_axis : 1-D ascending grid axes. - rhi_excess_cube, p_issr_cube : arrays of shape (n_lon, n_lat, n_pressure). - anchor : the geographic anchor (must match the one used for features/GUI). - threshold : RHi-excess cut for the pure-RHi mode (parity with ISSRField). - p_threshold : calibrated-probability cut for the default risk-aware mode. - mode : "prob" (is_inside = P(ISSR) >= p_threshold) or "rhi" - (is_inside = rhi_excess > threshold). - """ - - lon_axis: np.ndarray - lat_axis: np.ndarray - pressure_axis: np.ndarray - rhi_excess_cube: np.ndarray - p_issr_cube: np.ndarray - anchor: GeoAnchor - threshold: float = 0.3 - p_threshold: float = 0.5 - mode: str = "prob" - domain: tuple[float, float, float, float, float, float] = _SIM_DOMAIN - source: str = "ml" - - _rhi_interp: Any = field(default=None, repr=False) - _p_interp: Any = field(default=None, repr=False) - - def __post_init__(self) -> None: - from scipy.interpolate import RegularGridInterpolator - - axes = (np.asarray(self.lon_axis, dtype=float), - np.asarray(self.lat_axis, dtype=float), - np.asarray(self.pressure_axis, dtype=float)) - self._rhi_interp = RegularGridInterpolator( - axes, np.asarray(self.rhi_excess_cube, dtype=float), - bounds_error=False, fill_value=0.0, - ) - self._p_interp = RegularGridInterpolator( - axes, np.asarray(self.p_issr_cube, dtype=float), - bounds_error=False, fill_value=0.0, - ) - - # ------------------------------------------------------------------ # - # coordinate mapping - # ------------------------------------------------------------------ # - def _to_query_points(self, x_km, y_km, z_m) -> tuple[np.ndarray, tuple[int, ...]]: - """Map local (x_km, y_km, z_m) to interpolator points (lon, lat, p).""" - x = np.asarray(x_km, dtype=float) - y = np.asarray(y_km, dtype=float) - z = np.asarray(z_m, dtype=float) - shape = np.broadcast(x, y, z).shape - xb, yb, zb = np.broadcast_arrays(x, y, z) - lon, lat = self.anchor.local_to_geo(xb.ravel(), yb.ravel()) - pressure = altitude_to_pressure_hpa(zb.ravel()) - pts = np.column_stack([np.asarray(lon), np.asarray(lat), np.asarray(pressure)]) - return pts, shape - - # ------------------------------------------------------------------ # - # ISSRField interface — pointwise - # ------------------------------------------------------------------ # - def rhi_excess(self, x_km: float, y_km: float, z_m: float) -> float: - pts, _ = self._to_query_points(x_km, y_km, z_m) - return float(self._rhi_interp(pts)[0]) - - def p_issr(self, x_km: float, y_km: float, z_m: float) -> float: - pts, _ = self._to_query_points(x_km, y_km, z_m) - return float(self._p_interp(pts)[0]) - - def is_inside(self, x_km: float, y_km: float, z_m: float) -> bool: - if self.mode == "rhi": - return self.rhi_excess(x_km, y_km, z_m) > self.threshold - return self.p_issr(x_km, y_km, z_m) >= self.p_threshold - - # ------------------------------------------------------------------ # - # ISSRField interface — vectorized - # ------------------------------------------------------------------ # - def rhi_excess_grid(self, x_km: np.ndarray, y_km: np.ndarray, - z_m: np.ndarray) -> np.ndarray: - pts, shape = self._to_query_points(x_km, y_km, z_m) - return self._rhi_interp(pts).reshape(shape) - - def p_issr_grid(self, x_km: np.ndarray, y_km: np.ndarray, - z_m: np.ndarray) -> np.ndarray: - pts, shape = self._to_query_points(x_km, y_km, z_m) - return self._p_interp(pts).reshape(shape) - - def mask_grid(self, x_km: np.ndarray, y_km: np.ndarray, - z_m: np.ndarray) -> np.ndarray: - if self.mode == "rhi": - return self.rhi_excess_grid(x_km, y_km, z_m) > self.threshold - return self.p_issr_grid(x_km, y_km, z_m) >= self.p_threshold - - -# =========================================================================== -# Factory — build an MLIssrField for a region/time (plan §7) -# =========================================================================== - - -def ml_issr_field( - config: MLConfig | None = None, - met_source: str = "gfs", - time: str | None = None, - *, - model: Any = None, - allow_synthetic: bool = False, - grid_res_deg: float = 1.0, - seed: int = 0, -) -> MLIssrField: - """Build an `MLIssrField` for the configured region. - - met_source: - "gfs" / "era5" — operational: load a registered model and real weather - (where ISSR *will* be, from forecast, or replay from reanalysis). - "synthetic" — hermetic test/demo path: a clearly-labelled fake met - cube and a model trained on the synthetic fallback. Requires - allow_synthetic=True (the guard keeps fake weather out of real runs). - - Delegates the actual grid construction to `predict.py` so the data and - model plumbing lives in one place. - """ - from . import predict # lazy: avoids a heavy import at module load - from .config import DEFAULT_CONFIG - - cfg = config or DEFAULT_CONFIG - if met_source == "synthetic": - return predict.synthetic_ml_issr_field( - cfg, model=model, allow_synthetic=allow_synthetic, - grid_res_deg=grid_res_deg, seed=seed, - ) - return predict.real_ml_issr_field( - cfg, met_source=met_source, time=time, model=model, - grid_res_deg=grid_res_deg, - ) diff --git a/contrail_ml/model.py b/contrail_ml/model.py deleted file mode 100644 index 4715785..0000000 --- a/contrail_ml/model.py +++ /dev/null @@ -1,285 +0,0 @@ -""" -model.py — RHiCorrector: a regime-split, uncertainty-aware bias corrector. - -THE IDEA (plan §5) -================== -Weather models are dry-biased near the tropopause, and the bias behaves -differently in dry vs humid air. So we split on the NWP relative-humidity-over- -ice at 85 % and fit two sub-models: - - dry regime (rhi < 85): XGBoost regressor — fast, tabular. - humid regime (rhi >= 85): MLP regressor — where the ISSR signal and - the nonlinearity live. - -We predict the RESIDUAL delta = rhi_iagos - rhi_nwp (the correction), not the -absolute RHi — easier and more stable. Corrected RHi = rhi_nwp + delta_hat. - -A separate XGBoost CLASSIFIER head gives P(ISSR), which `calibrate.py` then -calibrates. A bootstrap ENSEMBLE of the regressor gives a predictive spread -(rhi_std), so downstream code knows where the model is unsure. - -PHYSICS-INFORMED, THE LEGITIMATE WAY -==================================== -The thermodynamics enter as engineered features (`rhi`, computed from T/q/p by -the shared features module), the regime split is a physical boundary, and -outputs are clamped to thermodynamic plausibility. There is no PDE-residual -loss — there is no governing PDE for this threshold problem, so a PINN would be -cargo-culting. - -The class exposes a scikit-learn-style surface: fit / predict / save / load -plus predict_proba_issr, so it slots into the training and serving code without -special-casing. -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any - -import numpy as np - -from .features import FEATURES, feature_matrix - -if TYPE_CHECKING: - import pandas as pd - - from .config import MLConfig - -# RHi is clamped to this range on output: 0 % (bone dry) up to a generous -# supersaturation ceiling. Real upper-tropospheric RHi rarely exceeds ~160 %. -_RHI_MIN, _RHI_MAX = 0.0, 180.0 - - -@dataclass -class RHiCorrector: - """Regime-split residual corrector with a calibrated-ready ISSR head. - - Hyperparameters default to small, CPU-friendly values; `from_config` - pulls the production values from an `MLConfig`. - """ - - regime_split_rhi: float = 85.0 - issr_rhi_threshold: float = 100.0 - n_ensemble: int = 10 - xgb_max_depth: int = 6 - xgb_n_estimators: int = 300 - xgb_learning_rate: float = 0.05 - mlp_hidden: tuple[int, ...] = (64, 32) - random_state: int = 42 - min_regime_samples: int = 50 # below this, fall back to one combined model - - # learned state (populated by fit) - _ensemble: list[dict[str, Any]] = field(default_factory=list, repr=False) - _classifier: Any = field(default=None, repr=False) - _fitted: bool = field(default=False, repr=False) - - # ------------------------------------------------------------------ # - @classmethod - def from_config(cls, cfg: MLConfig) -> RHiCorrector: - return cls( - regime_split_rhi=cfg.regime_split_rhi, - issr_rhi_threshold=cfg.issr_rhi_threshold, - n_ensemble=cfg.n_ensemble, - xgb_max_depth=cfg.xgb_max_depth, - xgb_n_estimators=cfg.xgb_n_estimators, - xgb_learning_rate=cfg.xgb_learning_rate, - mlp_hidden=cfg.mlp_hidden, - random_state=cfg.random_state, - ) - - # ------------------------------------------------------------------ # - # Fit - # ------------------------------------------------------------------ # - def fit(self, X: pd.DataFrame, y, sample_weight=None) -> RHiCorrector: - """Fit the corrector. - - X : DataFrame containing the FEATURES columns (must include `rhi`). - y : the residual target delta = rhi_iagos - rhi_nwp (per row). - sample_weight : optional per-row weights. - - The ISSR classification label is derived internally as - (rhi_nwp + delta >= issr_rhi_threshold), so the caller passes only the - residual — one target, sklearn-style. - """ - Xmat = feature_matrix(X).to_numpy(dtype=float) - rhi_nwp = X["rhi"].to_numpy(dtype=float) - delta = np.asarray(y, dtype=float) - w = None if sample_weight is None else np.asarray(sample_weight, dtype=float) - - # Derived classification label from the truth = nwp + residual. - y_issr = (rhi_nwp + delta >= self.issr_rhi_threshold).astype(int) - - rng = np.random.default_rng(self.random_state) - n = len(Xmat) - - # Bootstrap ensemble of the regime-split regressor. - self._ensemble = [] - for b in range(self.n_ensemble): - if b == 0: - idx = np.arange(n) # first member = full-data fit (stable mean) - else: - idx = rng.integers(0, n, size=n) - member = self._fit_regime_regressor( - Xmat[idx], rhi_nwp[idx], delta[idx], - None if w is None else w[idx], seed=self.random_state + b, - ) - self._ensemble.append(member) - - # Single classification head (calibrated later, externally). - self._classifier = self._fit_classifier(Xmat, y_issr, w) - self._fitted = True - return self - - def _fit_regime_regressor(self, Xmat, rhi_nwp, delta, w, seed) -> dict[str, Any]: - """Fit dry + humid residual sub-models, with a combined fallback when a - regime is too sparse to fit on its own.""" - from sklearn.neural_network import MLPRegressor - from xgboost import XGBRegressor - - humid = rhi_nwp >= self.regime_split_rhi - dry = ~humid - - def _xgb() -> Any: - return XGBRegressor( - max_depth=self.xgb_max_depth, - n_estimators=self.xgb_n_estimators, - learning_rate=self.xgb_learning_rate, - subsample=0.9, - colsample_bytree=0.9, - random_state=seed, - n_jobs=1, - verbosity=0, - ) - - def _mlp() -> Any: - return MLPRegressor( - hidden_layer_sizes=tuple(self.mlp_hidden), - max_iter=500, - random_state=seed, - early_stopping=False, - ) - - if dry.sum() < self.min_regime_samples or humid.sum() < self.min_regime_samples: - # Too few samples in some regime: one combined XGBoost model. - combined = _xgb() - _xgb_fit(combined, Xmat, delta, w) - return {"mode": "combined", "model": combined} - - dry_model = _xgb() - _xgb_fit(dry_model, Xmat[dry], delta[dry], None if w is None else w[dry]) - humid_model = _mlp() - # MLPRegressor has no sample_weight; weights apply to XGB only. - humid_model.fit(Xmat[humid], delta[humid]) - return {"mode": "split", "dry": dry_model, "humid": humid_model} - - def _fit_classifier(self, Xmat, y_issr, w) -> Any: - from xgboost import XGBClassifier - - clf = XGBClassifier( - max_depth=self.xgb_max_depth, - n_estimators=self.xgb_n_estimators, - learning_rate=self.xgb_learning_rate, - subsample=0.9, - colsample_bytree=0.9, - random_state=self.random_state, - n_jobs=1, - verbosity=0, - eval_metric="logloss", - ) - # Degenerate label (all one class) — fall back to a constant predictor. - if len(np.unique(y_issr)) < 2: - return _ConstantClassifier(float(y_issr.mean())) - if w is None: - clf.fit(Xmat, y_issr) - else: - clf.fit(Xmat, y_issr, sample_weight=w) - return clf - - # ------------------------------------------------------------------ # - # Predict - # ------------------------------------------------------------------ # - def _member_delta(self, member: dict[str, Any], Xmat, rhi_nwp) -> np.ndarray: - if member["mode"] == "combined": - return member["model"].predict(Xmat) - humid = rhi_nwp >= self.regime_split_rhi - out = np.empty(len(Xmat), dtype=float) - if (~humid).any(): - out[~humid] = member["dry"].predict(Xmat[~humid]) - if humid.any(): - out[humid] = member["humid"].predict(Xmat[humid]) - return out - - def predict(self, X: pd.DataFrame) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - """Return (rhi_hat, rhi_std, p_issr). - - rhi_hat : corrected RHi (%), ensemble mean, clamped to plausibility. - rhi_std : ensemble spread on corrected RHi (predictive uncertainty). - p_issr : raw P(ISSR) from the classifier head (calibrate externally). - """ - self._check_fitted() - Xmat = feature_matrix(X).to_numpy(dtype=float) - rhi_nwp = X["rhi"].to_numpy(dtype=float) - - members = np.stack( - [rhi_nwp + self._member_delta(m, Xmat, rhi_nwp) for m in self._ensemble] - ) - members = np.clip(members, _RHI_MIN, _RHI_MAX) - rhi_hat = members.mean(axis=0) - rhi_std = members.std(axis=0) - p_issr = self._classifier.predict_proba(Xmat)[:, 1] - return rhi_hat, rhi_std, p_issr - - def predict_corrected_rhi(self, X: pd.DataFrame) -> np.ndarray: - return self.predict(X)[0] - - def predict_proba_issr(self, X: pd.DataFrame) -> np.ndarray: - """Raw P(ISSR) in [0, 1] (one column).""" - self._check_fitted() - Xmat = feature_matrix(X).to_numpy(dtype=float) - return self._classifier.predict_proba(Xmat)[:, 1] - - # ------------------------------------------------------------------ # - # Persistence - # ------------------------------------------------------------------ # - def save(self, path: str) -> None: - import joblib - - joblib.dump(self, path) - - @staticmethod - def load(path: str) -> RHiCorrector: - import joblib - - obj = joblib.load(path) - if not isinstance(obj, RHiCorrector): - raise TypeError(f"{path} did not contain an RHiCorrector") - return obj - - # ------------------------------------------------------------------ # - def _check_fitted(self) -> None: - if not self._fitted: - raise RuntimeError("RHiCorrector.predict called before fit().") - - @property - def feature_names(self) -> tuple[str, ...]: - return FEATURES - - -def _xgb_fit(model, X, y, w) -> None: - if w is None: - model.fit(X, y) - else: - model.fit(X, y, sample_weight=w) - - -@dataclass -class _ConstantClassifier: - """Fallback when the training labels are single-class. Mimics the slice of - the scikit-learn API the corrector uses (`predict_proba[:, 1]`).""" - - p: float - - def predict_proba(self, X) -> np.ndarray: - n = len(X) - col1 = np.full(n, self.p) - return np.column_stack([1.0 - col1, col1]) diff --git a/contrail_ml/monitor.py b/contrail_ml/monitor.py deleted file mode 100644 index 4b0927b..0000000 --- a/contrail_ml/monitor.py +++ /dev/null @@ -1,112 +0,0 @@ -""" -monitor.py — Production monitoring with a real ground-truth feedback loop. - -This is not a static dashboard (plan §8). As new IAGOS flights arrive, we score -the model that was actually serving against the fresh truth (rolling ETS / -PR-AUC / ECE) AND watch the input distribution drift away from the training -reference (PSI per feature). When skill drops below a floor or drift breaches a -bound, we emit a structured `retrain_recommended` signal — the trigger an -orchestrator would act on. -""" - -from __future__ import annotations - -from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any - -import numpy as np -import pandas as pd - -from .evaluate import score_predictions -from .features import FEATURES - -if TYPE_CHECKING: - from .config import MLConfig - - -def population_stability_index( - expected: np.ndarray, actual: np.ndarray, n_bins: int = 10 -) -> float: - """PSI between a reference and a current sample of one feature. - - Rule of thumb: < 0.1 stable, 0.1-0.25 moderate shift, > 0.25 significant. - Bin edges come from the reference quantiles. - """ - expected = np.asarray(expected, dtype=float) - actual = np.asarray(actual, dtype=float) - edges = np.quantile(expected, np.linspace(0, 1, n_bins + 1)) - edges[0], edges[-1] = -np.inf, np.inf - e_hist, _ = np.histogram(expected, bins=edges) - a_hist, _ = np.histogram(actual, bins=edges) - e_frac = np.clip(e_hist / max(1, e_hist.sum()), 1e-6, None) - a_frac = np.clip(a_hist / max(1, a_hist.sum()), 1e-6, None) - return float(np.sum((a_frac - e_frac) * np.log(a_frac / e_frac))) - - -def feature_drift( - reference: pd.DataFrame, current: pd.DataFrame, features: tuple[str, ...] = FEATURES -) -> dict[str, float]: - """PSI per feature between the training reference and a current window.""" - return { - f: population_stability_index(reference[f].to_numpy(), current[f].to_numpy()) - for f in features - if f in reference.columns and f in current.columns - } - - -@dataclass -class MonitorReport: - n_new: int - skill: dict[str, float] - drift: dict[str, float] - retrain_recommended: bool - reasons: list[str] = field(default_factory=list) - - def to_dict(self) -> dict[str, Any]: - return { - "n_new": self.n_new, - "skill": self.skill, - "max_drift_psi": max(self.drift.values()) if self.drift else 0.0, - "drift": self.drift, - "retrain_recommended": self.retrain_recommended, - "reasons": self.reasons, - } - - -def monitor( - model: Any, - reference: pd.DataFrame, - current: pd.DataFrame, - cfg: MLConfig, - *, - min_ets: float = 0.2, - max_psi: float = 0.25, -) -> MonitorReport: - """Score the serving model on fresh labelled data and check feature drift. - - `current` must carry the labels (rhi_iagos, y_issr) so skill is measured - against truth, plus the FEATURES so the model can predict and drift can be - computed. - """ - reasons: list[str] = [] - - rhi_hat, _std, p = model.predict(current) - skill = score_predictions( - current["rhi_iagos"].to_numpy(), current["y_issr"].to_numpy(), rhi_hat, p, - p_threshold=cfg.issr_p_threshold, - ) - drift = feature_drift(reference, current) - - if skill["ets"] < min_ets: - reasons.append(f"ETS {skill['ets']:.3f} < floor {min_ets}") - worst = max(drift.items(), key=lambda kv: kv[1], default=(None, 0.0)) - if worst[1] > max_psi: - reasons.append(f"feature '{worst[0]}' PSI {worst[1]:.3f} > bound {max_psi}") - - return MonitorReport( - n_new=int(len(current)), - skill=skill, - drift=drift, - retrain_recommended=bool(reasons), - reasons=reasons, - ) diff --git a/contrail_ml/predict.py b/contrail_ml/predict.py deleted file mode 100644 index 130e5d7..0000000 --- a/contrail_ml/predict.py +++ /dev/null @@ -1,293 +0,0 @@ -""" -predict.py — Turn a trained model + a weather snapshot into an MLIssrField. - -Serving runs the model over a `(lon, lat, pressure)` grid and packs the -corrected-RHi-excess and calibrated P(ISSR) cubes into an `MLIssrField` -(issr_field.py), which the optimizer consumes through the ISSRField interface. - -By default we serve from a GFS FORECAST — operationally you want to know where -ISSR *will* be when the aircraft arrives, not where it is now. ERA5 replay is -available for retrospective studies. Every feature is computed by the shared -`features.add_derived_features`, restricted to the ERA5∩GFS variable set, so no -feature is ever missing at serve time. - -A guarded `synthetic` path builds a clearly-labelled fake met cube and a model -trained on the synthetic fallback, so `serve-check` and the integration tests -exercise the entire seam with no network and no credentials. -""" - -from __future__ import annotations - -import warnings -from dataclasses import dataclass -from typing import TYPE_CHECKING, Any - -import numpy as np -import pandas as pd - -from .features import RAW_NWP_COLUMNS, add_derived_features, rhi_excess -from .issr_field import MLIssrField - -if TYPE_CHECKING: - from .config import MLConfig - - -@dataclass -class MetCube: - """A weather snapshot on an ascending `(lon, lat, pressure)` grid. - - `fields` maps each raw NWP variable name (features.RAW_NWP_COLUMNS) to an - array of shape (n_lon, n_lat, n_pressure). `time` is a single valid time - applied to every cell. - """ - - lon_axis: np.ndarray - lat_axis: np.ndarray - pressure_axis: np.ndarray - fields: dict[str, np.ndarray] - time: Any # np.datetime64 / str / pandas Timestamp - - def shape(self) -> tuple[int, int, int]: - return (len(self.lon_axis), len(self.lat_axis), len(self.pressure_axis)) - - -# =========================================================================== -# Core: met cube + model -> MLIssrField -# =========================================================================== - - -def field_from_metcube( - model: Any, - cube: MetCube, - cfg: MLConfig, - *, - mode: str = "prob", -) -> MLIssrField: - """Run `model` over every cell of `cube` and build an MLIssrField.""" - nlon, nlat, npr = cube.shape() - lon3, lat3, p3 = np.meshgrid(cube.lon_axis, cube.lat_axis, cube.pressure_axis, - indexing="ij") - rows = { - "lon": lon3.ravel(), - "lat": lat3.ravel(), - "pressure_hpa": p3.ravel(), - "time": pd.to_datetime(cube.time), - } - for name in RAW_NWP_COLUMNS: - if name not in cube.fields: - raise KeyError(f"met cube missing required field {name!r}") - rows[name] = np.asarray(cube.fields[name], dtype=float).ravel() - - df = add_derived_features(pd.DataFrame(rows), pressure_col="pressure_hpa") - rhi_hat, _rhi_std, p_issr = model.predict(df) - - rhi_exc_cube = np.asarray(rhi_excess(rhi_hat), dtype=float).reshape(nlon, nlat, npr) - p_issr_cube = np.asarray(p_issr, dtype=float).reshape(nlon, nlat, npr) - - return MLIssrField( - lon_axis=np.asarray(cube.lon_axis, dtype=float), - lat_axis=np.asarray(cube.lat_axis, dtype=float), - pressure_axis=np.asarray(cube.pressure_axis, dtype=float), - rhi_excess_cube=rhi_exc_cube, - p_issr_cube=p_issr_cube, - anchor=cfg.anchor, - threshold=0.3, # RHi-excess cut for the pure-RHi mode (parity w/ ISSRField) - p_threshold=cfg.issr_p_threshold, - mode=mode, - ) - - -# =========================================================================== -# Real serving path (GFS forecast / ERA5 replay) -# =========================================================================== - - -def real_ml_issr_field( - cfg: MLConfig, - met_source: str = "gfs", - time: str | None = None, - *, - model: Any = None, - grid_res_deg: float = 1.0, - mode: str = "prob", -) -> MLIssrField: - """Load a registered model + real weather and build the field. - - Not hermetic: needs the `[ml]` geo stack and network/credentials. Kept thin - so the offline tests don't depend on it. - """ - if model is None: - from .registry import load_model - - model = load_model(cfg) - - if met_source == "gfs": - from .data.gfs import load_forecast_cube - - cube = load_forecast_cube(cfg, time=time, grid_res_deg=grid_res_deg) - elif met_source == "era5": - from .data.era5 import load_met_cube - - cube = load_met_cube(cfg, time=time, grid_res_deg=grid_res_deg) - else: - raise ValueError(f"unknown met_source {met_source!r} (use gfs/era5/synthetic)") - - return field_from_metcube(model, cube, cfg, mode=mode) - - -# =========================================================================== -# Hermetic synthetic path (tests / serve-check demo) — clearly labelled -# =========================================================================== - - -def synthetic_ml_issr_field( - cfg: MLConfig, - *, - model: Any = None, - allow_synthetic: bool = False, - grid_res_deg: float = 2.0, - seed: int = 0, - mode: str = "prob", -) -> MLIssrField: - """Build an MLIssrField from a fake met cube + fallback-trained model. - - For tests and offline demos ONLY; guarded so fake weather can't slip into a - real run. - """ - if not allow_synthetic: - raise RuntimeError( - "synthetic_ml_issr_field refused: pass allow_synthetic=True. The " - "synthetic met cube is for tests/serve-check demos only." - ) - if model is None: - model = _quick_fallback_model(cfg, seed=seed) - cube = synthetic_met_cube(cfg, grid_res_deg=grid_res_deg, seed=seed, - allow_synthetic=True) - return field_from_metcube(model, cube, cfg, mode=mode) - - -def _quick_fallback_model(cfg: MLConfig, seed: int = 0) -> Any: - """Train a small RHiCorrector on the synthetic fallback (tests/demo).""" - from .data.schema import LABEL_DELTA - from .data.synthetic_fallback import make_synthetic_training_table - from .model import RHiCorrector - - df = make_synthetic_training_table(n=3000, seed=seed, allow_synthetic=True) - m = RHiCorrector( - regime_split_rhi=cfg.regime_split_rhi, - issr_rhi_threshold=cfg.issr_rhi_threshold, - n_ensemble=3, - xgb_n_estimators=80, - min_regime_samples=30, - random_state=cfg.random_state, - ) - m.fit(df, df[LABEL_DELTA].to_numpy()) - return m - - -def synthetic_met_cube( - cfg: MLConfig, - grid_res_deg: float = 2.0, - seed: int = 0, - *, - allow_synthetic: bool = False, -) -> MetCube: - """Fabricate a physically-flavoured fake weather cube (tests/demo ONLY). - - Builds a few smooth humid 'blobs' in (lon, lat, pressure) whose NWP RHi - peaks just below saturation; the model's dry-bias correction then pushes - the blob cores past 100 % — a nice illustration of recovering hidden ISSR - from a dry-biased forecast. Auxiliary fields (cloud ice, vertical velocity, - cloud cover) are correlated with humidity the same way the training - fallback is, so the learned relationship transfers. - """ - if not allow_synthetic: - raise RuntimeError("synthetic_met_cube refused: pass allow_synthetic=True.") - warnings.warn( - "USING SYNTHETIC MET CUBE — fabricated weather, not a real forecast. " - "Valid for tests/serve-check demos only.", - stacklevel=2, - ) - from .config import SIM_DOMAIN - from .features import EPSILON, e_si_pa - - rng = np.random.default_rng(seed) - lon_axis = np.arange(cfg.lon_min, cfg.lon_max + 1e-6, grid_res_deg) - lat_axis = np.arange(cfg.lat_min, cfg.lat_max + 1e-6, grid_res_deg) - pressure_axis = np.array(sorted(cfg.pressure_levels_hpa), dtype=float) - lon3, lat3, p3 = np.meshgrid(lon_axis, lat_axis, pressure_axis, indexing="ij") - - # Temperature: colder aloft, smooth latitude gradient. - T = 225.0 - 0.05 * (p3 - 225.0) - 0.15 * (lat3 - 50.0) - - # Place the humid blobs INSIDE the geographic footprint of the local sim - # box (so the optimizer actually sees ISSR to route around), not scattered - # across the whole continental bbox. The footprint comes from the SAME - # anchor the field/GUI use. - x0, x1, y0, y1, _z0, _z1 = SIM_DOMAIN - corners = [(x0, y0), (x1, y0), (x0, y1), (x1, y1)] - geo = [cfg.anchor.local_to_geo(cx, cy) for cx, cy in corners] - box_lon_lo = min(g[0] for g in geo) - box_lon_hi = max(g[0] for g in geo) - box_lat_lo = min(g[1] for g in geo) - box_lat_hi = max(g[1] for g in geo) - - # Cruise pressures: only the levels the FL340-400 band actually maps onto - # (~187-250 hPa) carry blob cores, so the aircraft fly through the ISSR. - cruise_levels = [p for p in pressure_axis if 180.0 <= p <= 255.0] or [225.0] - # Keep cores comfortably inside the box footprint (margin > blob sigma). - lon_lo, lon_hi = box_lon_lo + 2.0, box_lon_hi - 2.0 - lat_lo, lat_hi = box_lat_lo + 1.5, box_lat_hi - 1.5 - - # A handful of Gaussian humid blobs -> a "true" RHi field. - true_rhi = np.full(lon3.shape, 35.0) - for _ in range(6): - c_lon = rng.uniform(lon_lo, lon_hi) - c_lat = rng.uniform(lat_lo, lat_hi) - c_p = float(rng.choice(cruise_levels)) - amp = rng.uniform(70.0, 100.0) - true_rhi += amp * np.exp( - -(((lon3 - c_lon) / 6.0) ** 2 - + ((lat3 - c_lat) / 4.0) ** 2 - + ((p3 - c_p) / 30.0) ** 2) - ) - true_rhi = np.clip(true_rhi, 5.0, 150.0) - - humid = np.clip(true_rhi - 90.0, 0.0, None) - cloud_ice = np.clip(8e-6 * humid / 60.0, 0.0, None) - w = -0.05 * np.clip(true_rhi - 100.0, 0.0, None) / 50.0 - cloud_cover = np.clip(0.15 + 0.6 * np.clip(true_rhi - 80.0, 0.0, None) / 70.0, 0.0, 1.0) - u = np.full(lon3.shape, 12.0) - v = np.zeros(lon3.shape) - - # Apply the SAME dry bias the training fallback uses, so the NWP cube is - # dry-biased and the model has the correction to apply. - bias = 2.0 + 0.28 * np.clip(true_rhi - 85.0, 0.0, None) + 1.5e6 * cloud_ice - nwp_rhi = np.clip(true_rhi - bias, 1.0, 200.0) - - esi = np.asarray(e_si_pa(T), dtype=float) - e = (nwp_rhi / 100.0) * esi - q_nwp = EPSILON * e / (p3 * 100.0 - (1.0 - EPSILON) * e) - - fields = { - "T": T, "specific_humidity": q_nwp, "u": u, "v": v, "w": w, - "cloud_cover": cloud_cover, "cloud_ice": cloud_ice, - } - return MetCube(lon_axis=lon_axis, lat_axis=lat_axis, pressure_axis=pressure_axis, - fields=fields, time=cfg.ref_time) - - -# =========================================================================== -# Tabular prediction (used by the monitor) -# =========================================================================== - - -def predict_table(model: Any, df: pd.DataFrame) -> pd.DataFrame: - """Append model outputs to a feature dataframe (corrected RHi, std, p_issr).""" - feat = add_derived_features(df) if "rhi" not in df.columns else df - rhi_hat, rhi_std, p_issr = model.predict(feat) - out = feat.copy() - out["rhi_hat"] = rhi_hat - out["rhi_std"] = rhi_std - out["p_issr"] = p_issr - return out diff --git a/contrail_ml/registry.py b/contrail_ml/registry.py deleted file mode 100644 index 15efe7d..0000000 --- a/contrail_ml/registry.py +++ /dev/null @@ -1,255 +0,0 @@ -""" -registry.py — MLflow experiment tracking + model registry (plan §8). - -A served artifact is the corrector PLUS its probability calibrator bundled as a -`ServedModel`, so serving applies the calibrated P(ISSR), not the raw head. - -`log_run` records everything needed to reproduce a model — params, metrics, the -comparison table, diagnostic plots, the git SHA, the dataset hash, and library -versions — logs the bundle as an artifact, registers it under -`cfg.registered_model_name`, and moves the new version to Staging. - -`load_model` pulls the latest registered version back for serving. The registry -defaults to a local file backend (`mlruns/`) and is configurable to a remote -tracking URI, so the same code runs on a laptop or against a team server. -""" - -from __future__ import annotations - -import subprocess -from dataclasses import dataclass -from typing import TYPE_CHECKING, Any - -import numpy as np -import pandas as pd - -if TYPE_CHECKING: - from .calibrate import ProbabilityCalibrator - from .config import MLConfig - from .model import RHiCorrector - - -@dataclass -class ServedModel: - """Serving bundle: the corrector + the probability calibrator. - - `.predict(df)` mirrors RHiCorrector.predict but returns the CALIBRATED - P(ISSR), so the whole serving path (predict.field_from_metcube) is unchanged - whether it gets a raw corrector or a calibrated bundle. - """ - - corrector: RHiCorrector - calibrator: ProbabilityCalibrator | None = None - conformal_half_width: float = float("nan") - - def predict(self, df) -> tuple[np.ndarray, np.ndarray, np.ndarray]: - rhi_hat, rhi_std, p_raw = self.corrector.predict(df) - p = self.calibrator.transform(p_raw) if self.calibrator is not None else p_raw - return rhi_hat, rhi_std, np.asarray(p, dtype=float) - - def predict_proba_issr(self, df) -> np.ndarray: - return self.predict(df)[2] - - def save(self, path: str) -> None: - import joblib - - joblib.dump(self, path) - - @staticmethod - def load(path: str) -> ServedModel: - import joblib - - obj = joblib.load(path) - if not isinstance(obj, ServedModel): - raise TypeError(f"{path} did not contain a ServedModel") - return obj - - -# =========================================================================== -# MLflow helpers -# =========================================================================== - -_ARTIFACT_DIR = "model" -_ARTIFACT_NAME = "served_model.joblib" - - -def set_tracking(cfg: MLConfig) -> None: - import os - - import mlflow - - # The local file backend (the spec's default) is in maintenance mode on - # newer MLflow and raises unless this opt-out is set. Honour it for file - # URIs; a remote/db tracking URI is unaffected. - if cfg.mlflow_tracking_uri.startswith(("file:", "./", "../", "mlruns")): - os.environ.setdefault("MLFLOW_ALLOW_FILE_STORE", "true") - mlflow.set_tracking_uri(cfg.mlflow_tracking_uri) - mlflow.set_experiment(cfg.registered_model_name) - - -def _git_sha() -> str: - try: - return subprocess.check_output( - ["git", "rev-parse", "HEAD"], stderr=subprocess.DEVNULL - ).decode().strip() - except Exception: - return "unknown" - - -def _lib_versions() -> dict[str, str]: - import importlib - - out = {} - for p in ("numpy", "pandas", "scipy", "sklearn", "xgboost", "mlflow"): - try: - out[p] = importlib.import_module(p).__version__ # type: ignore[attr-defined] - except Exception: - out[p] = "absent" - return out - - -def log_run( - cfg: MLConfig, - served_model: ServedModel, - *, - params: dict[str, Any], - metrics: dict[str, float], - comparison: pd.DataFrame | None = None, - figures: dict[str, Any] | None = None, - dataset_hash: str = "unknown", - register: bool = True, -) -> dict[str, Any]: - """Log a training run to MLflow and (optionally) register + stage it. - - Returns {run_id, model_version (or None)}. - """ - import json - import os - import tempfile - - import mlflow - - set_tracking(cfg) - result: dict[str, Any] = {"run_id": None, "model_version": None} - - with mlflow.start_run() as run: - result["run_id"] = run.info.run_id - mlflow.log_params({k: _short(v) for k, v in params.items()}) - mlflow.set_tags( - {"git_sha": _git_sha(), "dataset_hash": dataset_hash, **_lib_versions()} - ) - mlflow.log_metrics({k: float(v) for k, v in metrics.items() - if v is not None and np.isfinite(v)}) - - # Log the bundle as a proper pyfunc MODEL (so it is registrable on - # modern MLflow, which requires a logged model — a raw artifact path is - # no longer enough). - info = mlflow.pyfunc.log_model( - artifact_path=_ARTIFACT_DIR, - python_model=_served_pyfunc(served_model), - ) - - with tempfile.TemporaryDirectory() as td: - if comparison is not None: - cpath = os.path.join(td, "comparison_table.csv") - comparison.to_csv(cpath) - mlflow.log_artifact(cpath) - try: # markdown is nicer but needs `tabulate`; CSV is the source of truth - mlflow.log_text(comparison.to_markdown(), "comparison_table.md") - except ImportError: - pass - - meta = {"dataset_hash": dataset_hash, "git_sha": _git_sha(), - "config": _config_dict(cfg)} - jpath = os.path.join(td, "run_metadata.json") - with open(jpath, "w", encoding="utf-8") as fh: - json.dump(meta, fh, indent=2, default=str) - mlflow.log_artifact(jpath) - - for name, fig in (figures or {}).items(): - mlflow.log_figure(fig, f"plots/{name}.png") - - if register: - mv = mlflow.register_model(info.model_uri, cfg.registered_model_name) - result["model_version"] = mv.version - _transition(cfg, mv.version, "Staging") - - return result - - -def _served_pyfunc(served: ServedModel): - """Wrap a ServedModel as an mlflow PythonModel (predict -> DataFrame). - - Defined inside a factory so `import mlflow` stays lazy; cloudpickle captures - the served bundle (the fitted estimators) with the model. - """ - import mlflow - - class _ServedPyfunc(mlflow.pyfunc.PythonModel): # type: ignore[misc,name-defined] - def predict(self, context, model_input, params=None): # noqa: ARG002 - rhi_hat, rhi_std, p = served.predict(model_input) - return pd.DataFrame({"rhi_hat": rhi_hat, "rhi_std": rhi_std, "p_issr": p}) - - return _ServedPyfunc() - - -def _transition(cfg: MLConfig, version: str, stage: str) -> None: - from mlflow.tracking import MlflowClient - - MlflowClient().transition_model_version_stage( - name=cfg.registered_model_name, version=version, stage=stage, - archive_existing_versions=False, - ) - - -def promote(cfg: MLConfig, version: str, stage: str = "Production") -> None: - """Promote a registered version to a stage (Staging/Production/Archived).""" - set_tracking(cfg) - _transition(cfg, version, stage) - - -@dataclass -class _PyfuncAdapter: - """Adapts a loaded pyfunc model back to the `.predict -> (rhi_hat, rhi_std, - p_issr)` tuple interface the serving code (predict.field_from_metcube) and - the monitor expect.""" - - pyfunc: Any - - def predict(self, df) -> tuple: - out = self.pyfunc.predict(df) - return (out["rhi_hat"].to_numpy(), out["rhi_std"].to_numpy(), - out["p_issr"].to_numpy()) - - def predict_proba_issr(self, df): - return self.predict(df)[2] - - -def load_model(cfg: MLConfig, stage: str = "Staging") -> _PyfuncAdapter: - """Load the latest registered model for a stage (default Staging) and adapt - it to the native serving interface.""" - import mlflow - from mlflow.tracking import MlflowClient - - set_tracking(cfg) - client = MlflowClient() - versions = client.get_latest_versions(cfg.registered_model_name, stages=[stage]) - if not versions: - raise RuntimeError( - f"no registered model '{cfg.registered_model_name}' in stage {stage!r}. " - f"Train one first: python -m contrail_ml train" - ) - mv = versions[0] - loaded = mlflow.pyfunc.load_model(f"models:/{cfg.registered_model_name}/{mv.version}") - return _PyfuncAdapter(loaded) - - -def _short(v: Any) -> Any: - s = str(v) - return s[:250] - - -def _config_dict(cfg: MLConfig) -> dict[str, Any]: - return {f.name: getattr(cfg, f.name) - for f in cfg.__dataclass_fields__.values() # type: ignore[attr-defined] - if f.name != "extra"} diff --git a/contrail_ml/train.py b/contrail_ml/train.py deleted file mode 100644 index 65fd063..0000000 --- a/contrail_ml/train.py +++ /dev/null @@ -1,208 +0,0 @@ -""" -train.py — Cross-validate, fit, calibrate, evaluate, and register (plan §6). - -Pipeline -======== -1. Temporal split: train on older years, hold out recent years untouched, so - the headline metrics are measured on data the model never saw. -2. Grouped CV: GroupKFold over (year, ~5x5deg spatial tile) so neither time nor - place leaks across folds — the optimistic-CV trap for spatial data. -3. Fit the corrector on a fit-subset; calibrate P(ISSR) and fit the conformal - RHi interval on a disjoint calibration-subset (no calibration leakage). -4. Evaluate the calibrated model AND every baseline on the temporal holdout; - build the comparison table. -5. Log params, metrics, the table, and diagnostic plots to MLflow; register the - model and move it to Staging. - -`run_training(..., log_mlflow=False)` runs everything except the MLflow side, so -the training logic is unit-testable with no tracking server. -""" - -from __future__ import annotations - -import hashlib -from dataclasses import dataclass, field -from typing import TYPE_CHECKING, Any - -import numpy as np -import pandas as pd - -from .baselines import default_baselines -from .calibrate import ConformalRHi, ProbabilityCalibrator -from .data.schema import LABEL_DELTA, LABEL_ISSR, LABEL_RHI -from .evaluate import comparison_table, plot_reliability, plot_rhi_scatter, score_predictions -from .model import RHiCorrector -from .registry import ServedModel - -if TYPE_CHECKING: - from .config import MLConfig - - -@dataclass -class TrainResult: - served: ServedModel - conformal: ConformalRHi - comparison: pd.DataFrame - cv_metrics: dict[str, float] - holdout_metrics: dict[str, float] - dataset_hash: str - mlflow: dict[str, Any] = field(default_factory=dict) - - -# =========================================================================== -# splits & grouping -# =========================================================================== - - -def _year(df: pd.DataFrame) -> np.ndarray: - return pd.to_datetime(df["time"]).dt.year.to_numpy() - - -def temporal_split(df: pd.DataFrame, cfg: MLConfig) -> tuple[pd.DataFrame, pd.DataFrame]: - """Older years -> train, recent years -> held-out test. Falls back to a - random 75/25 split if the year ranges don't select anything.""" - yr = _year(df) - tr = (yr >= cfg.train_years[0]) & (yr <= cfg.train_years[1]) - te = (yr >= cfg.test_years[0]) & (yr <= cfg.test_years[1]) - if tr.sum() < 20 or te.sum() < 20: - rng = np.random.default_rng(cfg.random_state) - mask = rng.random(len(df)) < 0.75 - return df.loc[mask].reset_index(drop=True), df.loc[~mask].reset_index(drop=True) - return df.loc[tr].reset_index(drop=True), df.loc[te].reset_index(drop=True) - - -def cv_groups(df: pd.DataFrame) -> np.ndarray: - """Group label = year + ~5x5deg spatial tile (no time/space leakage).""" - yr = _year(df) - tile_lat = np.floor(df["lat"].to_numpy() / 5.0).astype(int) - tile_lon = np.floor(df["lon"].to_numpy() / 5.0).astype(int) - return np.array([f"{a}_{b}_{c}" for a, b, c in zip(yr, tile_lat, tile_lon, strict=False)]) - - -def dataset_hash(df: pd.DataFrame) -> str: - """Content hash of the training table for provenance.""" - h = hashlib.sha256(pd.util.hash_pandas_object(df, index=True).values.tobytes()) - return h.hexdigest()[:16] - - -# =========================================================================== -# CV & fit -# =========================================================================== - - -def cross_validate(cfg: MLConfig, train_df: pd.DataFrame, n_splits: int = 4) -> dict[str, float]: - """Grouped CV; returns mean RHi MAE and ISSR ETS across folds. Uses a - single-member corrector (no ensemble) for speed — CV measures signal, not - predictive spread.""" - from sklearn.model_selection import GroupKFold - - groups = cv_groups(train_df) - n_groups = len(np.unique(groups)) - n_splits = max(2, min(n_splits, n_groups)) - gkf = GroupKFold(n_splits=n_splits) - - maes, etss = [], [] - for tr_idx, va_idx in gkf.split(train_df, groups=groups): - tr, va = train_df.iloc[tr_idx], train_df.iloc[va_idx] - m = _new_corrector(cfg, n_ensemble=1) - m.fit(tr, tr[LABEL_DELTA].to_numpy()) - rhi_hat, _s, p = m.predict(va) - s = score_predictions(va[LABEL_RHI].to_numpy(), va[LABEL_ISSR].to_numpy(), - rhi_hat, p) - maes.append(s["rhi_mae"]) - etss.append(s["ets"]) - return {"rhi_mae": float(np.mean(maes)), "ets": float(np.mean(etss)), - "n_splits": float(n_splits)} - - -def fit_and_calibrate( - cfg: MLConfig, train_df: pd.DataFrame, calib_frac: float = 0.25 -) -> tuple[ServedModel, ConformalRHi]: - """Fit the corrector on a fit-subset, then calibrate P(ISSR) and the - conformal RHi interval on a disjoint, group-held-out calibration subset.""" - from sklearn.model_selection import GroupShuffleSplit - - groups = cv_groups(train_df) - gss = GroupShuffleSplit(n_splits=1, test_size=calib_frac, random_state=cfg.random_state) - fit_idx, cal_idx = next(gss.split(train_df, groups=groups)) - fit_df, cal_df = train_df.iloc[fit_idx], train_df.iloc[cal_idx] - - corrector = _new_corrector(cfg) - corrector.fit(fit_df, fit_df[LABEL_DELTA].to_numpy()) - - rhi_hat, _s, p_raw = corrector.predict(cal_df) - calibrator = ProbabilityCalibrator(method="isotonic").fit( - p_raw, cal_df[LABEL_ISSR].to_numpy() - ) - conformal = ConformalRHi(alpha=0.1).fit(cal_df[LABEL_RHI].to_numpy(), rhi_hat) - - served = ServedModel(corrector=corrector, calibrator=calibrator, - conformal_half_width=conformal.half_width) - return served, conformal - - -def _new_corrector(cfg: MLConfig, n_ensemble: int | None = None) -> RHiCorrector: - m = RHiCorrector.from_config(cfg) - if n_ensemble is not None: - m.n_ensemble = n_ensemble - return m - - -# =========================================================================== -# top-level -# =========================================================================== - - -def run_training( - cfg: MLConfig, - df: pd.DataFrame, - *, - log_mlflow: bool = True, - do_cv: bool = True, -) -> TrainResult: - """Full training run. Returns a TrainResult; logs to MLflow unless disabled.""" - from .data.schema import validate_table - - validate_table(df) - dh = dataset_hash(df) - train_df, test_df = temporal_split(df, cfg) - - cv_metrics = cross_validate(cfg, train_df) if do_cv else {} - served, conformal = fit_and_calibrate(cfg, train_df) - - baselines = [b.fit(train_df) for b in default_baselines(cfg.issr_rhi_threshold)] - comparison = comparison_table(test_df, served, baselines, - p_threshold=cfg.issr_p_threshold) - - # Headline metrics = the model's holdout row. - rhi_hat, _std, p = served.predict(test_df) - holdout = score_predictions(test_df[LABEL_RHI].to_numpy(), - test_df[LABEL_ISSR].to_numpy(), rhi_hat, p, - p_threshold=cfg.issr_p_threshold) - - result = TrainResult(served=served, conformal=conformal, comparison=comparison, - cv_metrics=cv_metrics, holdout_metrics=holdout, dataset_hash=dh) - - if log_mlflow: - from .registry import log_run - - figures = { - "reliability": plot_reliability(p, test_df[LABEL_ISSR].to_numpy()), - "rhi_scatter": plot_rhi_scatter(test_df[LABEL_RHI].to_numpy(), rhi_hat, - test_df["rhi"].to_numpy()), - } - params = { - "regime_split_rhi": cfg.regime_split_rhi, - "n_ensemble": cfg.n_ensemble, - "xgb_n_estimators": cfg.xgb_n_estimators, - "xgb_max_depth": cfg.xgb_max_depth, - "issr_p_threshold": cfg.issr_p_threshold, - "n_train_rows": len(train_df), - "n_test_rows": len(test_df), - } - metrics = {**{f"cv_{k}": v for k, v in cv_metrics.items()}, - **{f"holdout_{k}": v for k, v in holdout.items()}} - result.mlflow = log_run(cfg, served, params=params, metrics=metrics, - comparison=comparison, figures=figures, - dataset_hash=dh) - return result diff --git a/docs/DATA.md b/docs/DATA.md deleted file mode 100644 index ab0095d..0000000 --- a/docs/DATA.md +++ /dev/null @@ -1,136 +0,0 @@ -# Data sources for `contrail_ml` - -The ISSR model trains on three free-to-academic data sources. None of them are -committed to the repo, and **no credentials live in the codebase** — paths and -keys come from the environment (see `.env.example`). The hermetic tests use the -guarded `synthetic_fallback` and need none of this. - -| Source | Role | Cost | Access | -|--------|------|------|--------| -| **IAGOS** | ground-truth in-situ humidity (labels) | free, registration | portal sign-up | -| **ARCO-ERA5** | reanalysis weather (training inputs) | free | public GCS bucket, no creds | -| **GFS** | forecast weather (serving inputs) | free, public domain | NOAA, via pycontrails | -| ECMWF HRES | best forecast (optional, not used) | paid/licensed | — | - -## IAGOS (the labels) - -IAGOS (In-service Aircraft for a Global Observing System, -) equips airliners with research humidity sensors. Its -measured relative-humidity-over-ice (RHi) near the tropopause is the truth this -project trains against. - -1. Register (free) on the IAGOS data portal and download the per-flight NetCDF - files for the years/region you want. -2. Put them in a folder and point the loader at it: - ``` - export CONTRAIL_ML_IAGOS_DIR=/path/to/iagos/netcdf - ``` -3. `contrail_ml/data/iagos.py` recomputes RHi from the measured T, water-vapour - and pressure with the SAME thermodynamics as the features (no label skew). - -**Variable names differ between IAGOS products** (IAGOS-CORE, MOZAIC, CARIBIC). -The defaults in `IAGOS_VARIABLES` target IAGOS-CORE; override per product via -`cfg.extra["iagos_variables"]`. Verify the names in your NetCDF before a real run. - -## ARCO-ERA5 (the training inputs) - -ERA5 is ECMWF's reanalysis — the best reconstruction of *past* weather. We read -the **ARCO-ERA5** zarr mirror on a public Google Cloud bucket through -pycontrails, which needs **no Copernicus CDS account**. `era5.py` pulls the -fields, `collocate.py` interpolates them onto the IAGOS waypoints. - -Requires the `[ml]` extra (`pycontrails`, `xarray`, `zarr`, `gcsfs`). - -## GFS (the serving inputs) - -At serve time you need where ISSR *will* be, so the operational input is NOAA's -**GFS forecast** — free, public-domain, refreshed every 6 h, delivered as GRIB -(hence `cfgrib`). `gfs.py` loads it via pycontrails' `GFSForecast`. - -## Config keys - -All configuration is `contrail_ml.config.MLConfig`, overridable from -`CONTRAIL_ML_*` environment variables (see `.env.example`) or a YAML file. Key -ones: `iagos_dir`, `mlflow_tracking_uri`, `pressure_levels_hpa`, the lon/lat -box, `train_years`/`test_years`, and the geographic anchor (`origin_lat`, -`origin_lon`). - -## Reproducibility - -`build-dataset` writes a **content-addressed** parquet plus a `dataset_card.json` -(row count, time/region coverage, class balance, dry-bias). Every training run -records that dataset hash, so any reported number is traceable to an exact table. -If DVC is installed the parquet is `dvc add`-ed automatically. - ---- - -# Real flights for `contrail_flights` (OpenSky) - -`contrail_flights` is an **independent** sibling package (it shares no imports -with `contrail_ml`). It turns real flown **historical** ADS-B tracks into the -same `contrail_env.Flight` objects the synthetic generator produces, so the -optimizer, solvers, and GUI map consume them unchanged. It is opt-in via -`flight_source="real"`; the default stays synthetic. Requires the `[flights]` -extra (`pyopensky`, `pandas`, `pyarrow`); the hermetic tests need none of it. - -> **Honest naming.** OpenSky gives real flown *historical* traffic (actual ADS-B -> tracks), **not** published schedules or planned flights. The accurate claim is -> "demonstrated on real historical European traffic" — the same evaluation method -> used in real contrail-avoidance trials (Google / American Airlines) — not -> "schedule optimization". - -## OpenSky access - -| Path | Window | Auth | Use | -|------|--------|------|-----| -| Public REST API | recent rolling ~1–2 h | none (rate-limited) | small/recent pulls | -| Research / Trino | arbitrary history, whole regions | account required | bulk historical | - -1. **Recent/light:** the unauthenticated REST API covers only a short rolling - window and is heavily rate-limited — fine for a quick demo, not for history. -2. **Bulk historical:** apply for an OpenSky **research/Trino account** at - (free for academic use). Put credentials in the - environment — never in the repo: - ``` - export OPENSKY_USERNAME=... - export OPENSKY_PASSWORD=... - ``` -3. Pulled tracks are cached as parquet (keyed by bbox + time window) under - `cache_dir`, since bulk queries are rate-limited/credit-metered. - -**No fabrication.** If `pyopensky` is missing, credentials are absent, or the -query fails/returns empty, `OpenSkyClient` raises `OpenSkyUnavailableError` with -an actionable message — it never returns invented or zero traffic. - -## Aircraft type → performance profile - -The optimizer needs an `Aircraft` performance model per flight. OpenSky's free -aircraft-database CSV maps `icao24` → ICAO type designator; we map that type to -the **nearest existing** profile via an explicit table -(`reduce_to_flight.AIRCRAFT_TYPE_TO_FACTORY`): - -| ICAO type | Profile | Note | -|-----------|---------|------| -| A319 / A320 / A321 | `a320_like` | linear narrowbody surrogate | -| B737 / B738 | `a320_like` | same class | -| *(anything else)* | `a320_like` (default) | documented approximation | - -**Limitation:** we do **not** synthesize a new performance model per type — every -mapped type currently points at the one `a320_like` linear surrogate. The table -is the single place to add real BADA-class profiles later. - -## Baseline = the real flown profile - -For real flights the `baseline` is, by preference, the **observed** altitude -profile: the track's barometric altitude is snapped to the RVSM flight-level grid -and resampled into `AltitudeSegment`s. Only when altitude data is too sparse/noisy -do we fall back to the synthetic fuel-optimal `build_baseline_profile`, and -`build_dataset` logs how many flights hit that fallback. - -## Config keys - -`contrail_flights.config.FlightsConfig`, overridable from `CONTRAIL_FLIGHTS_*` -env vars (plus `OPENSKY_USERNAME`/`OPENSKY_PASSWORD`). Key ones: the lon/lat -`bbox`, `start_time`/`end_time`, `min_track_points`/`min_track_km`, -`snapshot_window_s`, `cache_dir`, and the anchor (`origin_lat`, `origin_lon`, -which **must** match `MLConfig` so flights and the predicted field align). diff --git a/docs/ML.md b/docs/ML.md deleted file mode 100644 index 35cef8c..0000000 --- a/docs/ML.md +++ /dev/null @@ -1,105 +0,0 @@ -# The ISSR model (`contrail_ml`) - -This replaces the synthetic Gaussian-blob ISSR field with a machine-learning -model that predicts ice-supersaturated regions (ISSRs) from real weather, and -wraps it in a full MLOps lifecycle. The optimizer is unchanged — the model -plugs in through the same `ISSRField` interface. - -## Why this is a real problem - -A contrail persists (and warms) only where the air is **ice-supersaturated**: -relative humidity over ice `RHi >= 100 %`. The trouble is that weather models -carry a well-documented **dry bias in RHi near the tropopause** — exactly where -ISSRs form — so raw ERA5 has only a weak ISSR skill (equitable threat score -~0.2–0.4 against in-situ IAGOS measurements). The standard fix is to -**bias-correct RHi against observations**. This project builds a calibrated ML -version of that correction. - -## The pipeline - -``` -IAGOS truth ─┐ - ├─ collocate ─► training table ─► train ─► register ─► serve ─► monitor -ARCO-ERA5 ──┘ (parquet) (MLflow) (MLIssrField) -``` - -1. **Features** (`features.py`) — one shared module computes RHi (Murphy–Koop - ice-saturation), altitude↔pressure, cyclical time encodings, and the - local↔geo transform. Used identically at train and serve time, so there is - **no training/serving skew**. -2. **Model** (`model.py`) — `RHiCorrector` is **regime-split**: an XGBoost - regressor in dry air, a small MLP in the humid regime (split at 85 % RHi, - where the ISSR signal lives). It predicts the **residual** - `delta = RHi_IAGOS − RHi_NWP` (the correction), so corrected RHi = - `RHi_NWP + delta_hat`. A bootstrap ensemble gives a predictive spread; an - XGBoost head gives `P(ISSR)`. -3. **Calibration** (`calibrate.py`) — isotonic regression makes `P(ISSR)` honest - (target ECE < 0.05); split-conformal gives a distribution-free RHi interval. -4. **Serving** (`predict.py`, `issr_field.py`) — the model runs over a - `(lon, lat, pressure)` grid; `MLIssrField` wraps the result and answers the - optimizer's `is_inside` / `rhi_excess` queries by interpolation. Default - serving source is the **GFS forecast** (operational); ERA5 replay is - available for retrospective studies. - -### Physics-informed, the legitimate way - -The thermodynamics enter as **engineered features** (`rhi`, `e_si(T)`), the -regime split is a **physical boundary**, and outputs are clamped to -thermodynamic plausibility. There is no PDE-residual loss — there is no -governing PDE for this threshold problem, so a PINN would be cargo-culting. - -## Reading the comparison table - -`train`/`evaluate` print one row per predictor on the temporal holdout: - -| column | meaning | better | -|--------|---------|--------| -| `rhi_mae`, `rhi_rmse` | corrected-RHi error vs IAGOS | lower | -| `rhi_bias` | mean error — **the dry-bias number** | nearer 0 | -| `ets` | ISSR equitable threat score | higher | -| `f1`, `roc_auc`, `pr_auc` | ISSR detection skill | higher | -| `brier`, `ece` | probability honesty / calibration | lower | - -The model is compared against three baselines it must beat to justify itself: -`raw_era5` (uncorrected), `x_factor` (one global scale), and `quantile_map` -(bivariate T/RHi quantile mapping — the standard statistical correction). The -headline claim is only ever "corrected RHi reduces the dry bias and lifts ISSR -ETS over raw ERA5, with calibrated probabilities" — **reported honestly whatever -the numbers are**. - -## Known limitations (stated, not hidden) - -- **IAGOS sampling bias** — aircraft avoid deep convection, so the truth set - under-samples the most intense humidity. Surfaced as an optional - `sample_weight`; documented, not pretended away. -- **Reanalysis→forecast gap** — we train on ERA5 reanalysis but serve from GFS - forecast; the domain shift is a real, measured effect (watch it with - `monitor`). -- **No guarantee of beating quantile mapping** — on real data the gains may be - modest. The rigorous benchmark + calibration is itself the contribution. -- The bundled `serve-check`/tests use a **synthetic** met cube (clearly - labelled, guarded) so the seam is exercisable offline. It is never a stand-in - for real weather. - -## MLOps lifecycle - -- **Data versioning** — content-addressed parquet + dataset card (+ DVC if set up). -- **Tracking + registry** — MLflow (`mlruns/` by default): params, metrics, the - comparison table, reliability/scatter plots, git SHA, dataset hash, library - versions. Model registered as `contrail-issr-rhi-corrector`, moved to Staging. -- **Monitoring** (`monitor.py`) — as new IAGOS arrives, rolling ETS/PR-AUC/ECE - against fresh truth + per-feature PSI drift, emitting a `retrain_recommended` - signal. A real ground-truth feedback loop, not a static dashboard. - -## Commands - -``` -python -m contrail_ml build-dataset # IAGOS+ERA5 -> versioned parquet -python -m contrail_ml build-dataset --synthetic # offline fallback table -python -m contrail_ml train --data # CV, fit, calibrate, register -python -m contrail_ml train --synthetic --no-mlflow # hermetic dry run -python -m contrail_ml evaluate --data # model-vs-baselines table -python -m contrail_ml predict --synthetic # build + save an ISSR field -python -m contrail_ml serve-check # full seam: ML field -> CP-SAT -python -m contrail_ml monitor --reference a.parquet --current b.parquet -``` diff --git a/gui/app.py b/gui/app.py index ff88bc7..b1b871b 100644 --- a/gui/app.py +++ b/gui/app.py @@ -86,6 +86,7 @@ from contrail_env import assemble_qubo, fl_to_m, waypoints_for from contrail_env.benchmark import SOLVER_NAMES, run_benchmark +from contrail_env.geo import EUROPEAN_ANCHOR from contrail_env.pasqal_analog import MAX_STATEVECTOR_QUBITS from service.client import DEFAULT_SERVER_ADDRESS, SolverClient from service.generated import solver_pb2 @@ -110,6 +111,11 @@ _MAP_GRID_NX, _MAP_GRID_NY = 100, 68 _MAP_MARKER_SIZE = 9 _MAP_RISK_FLOOR_FRAC = 0.04 +# Animation: number of frames swept across the planning window, and the +# per-frame dwell (ms) when playing. The route "reveals" (its trail grows) and +# the aircraft marker glides along it as the time cursor advances. +_MAP_ANIM_FRAMES = 48 +_MAP_FRAME_MS = 90 # Bundled Natural Earth vectors for the geo basemap (see _geo_assets_script). _GEO_TOPOJSON_NAME = "world_50m" _GEO_ASSET_PATH = Path(__file__).resolve().parent / "assets" / f"{_GEO_TOPOJSON_NAME}.json" @@ -141,6 +147,24 @@ def _geo_assets_script() -> str: ) +def _autoplay_script() -> str: + """Auto-run the route-reveal animation once the plot is live. + + Plotly doesn't autoplay frames, so we poll for the graph div + library and + then kick off Plotly.animate from the first frame. Because _render_map + rewrites this page on every solve, the reveal replays each time a solution + lands. No-op for a static (frame-less) figure. + """ + return ( + "" + ) + + def _mean_best_cost(instances, solver: str) -> str: """Mean best cost of one solver across the benchmark instances.""" costs = [ @@ -158,14 +182,67 @@ def _mean_best_cost(instances, solver: str) -> str: @dataclass(frozen=True) class RouteLine: - """One flight's ground track in geographic coordinates for the map.""" + """One flight's ground track in geographic coordinates for the map. + + `t_s` is the world time (s) at each vertex; it drives the animation (the + trail reveals for vertices with t <= cursor, the aircraft marker rides the + interpolated position at the cursor). + """ name: str lon: np.ndarray lat: np.ndarray + t_s: np.ndarray chosen: bool +def _route_state_at(r: RouteLine, t_cursor: float): + """A route's revealed trail + aircraft-head position at time `t_cursor`. + + The trail is every vertex flown so far (t <= cursor) plus the interpolated + current point, so the line ends exactly under the moving marker. np.interp + clamps outside the flight window, so before departure the head sits at the + origin and after arrival at the destination. + """ + t = np.asarray(r.t_s, dtype=float) + lon = np.asarray(r.lon, dtype=float) + lat = np.asarray(r.lat, dtype=float) + h_lon = float(np.interp(t_cursor, t, lon)) + h_lat = float(np.interp(t_cursor, t, lat)) + flown = t <= t_cursor + trail_lon = np.append(lon[flown], h_lon) + trail_lat = np.append(lat[flown], h_lat) + return trail_lon, trail_lat, h_lon, h_lat + + +def _trail_trace(go, r: RouteLine, trail_lon, trail_lat): + """The growing ground-track line (bright green + thick if chosen).""" + return go.Scattergeo( + lon=trail_lon, lat=trail_lat, mode="lines", + line=dict( + width=4 if r.chosen else 1.4, + color="#28dc5a" if r.chosen else "rgba(170,172,184,0.5)", + ), + name=f"{r.name} (chosen)" if r.chosen else r.name, + hovertemplate=f"{r.name}", + ) + + +def _head_trace(go, r: RouteLine, h_lon, h_lat): + """The aircraft marker riding the front of the trail.""" + return go.Scattergeo( + lon=[h_lon], lat=[h_lat], mode="markers", + marker=dict( + size=14 if r.chosen else 8, + symbol="triangle-up", + color="#39ff7a" if r.chosen else "rgba(210,212,224,0.85)", + line=dict(width=1.2 if r.chosen else 0.6, color="#0e0e10"), + ), + name=r.name, showlegend=False, + hovertemplate=f"{r.name}", + ) + + def build_map_figure( *, source: str, @@ -173,17 +250,25 @@ def build_map_figure( lat: np.ndarray, risk: np.ndarray, routes: list[RouteLine], + animate: bool = True, ): - """Build the geographic Map figure: ISSR risk overlay + flight routes. + """Build the geographic Map figure: ISSR risk overlay + animated flight routes. Pure function of plain arrays so it carries no Qt/web dependency and can be unit-tested directly. `lon`/`lat`/`risk` are the (same-shape) sampled risk - field; `routes` are per-flight ground tracks already in lon/lat. + field; `routes` are per-flight ground tracks in lon/lat with per-vertex + times (`RouteLine.t_s`). + + When `animate` and the routes span a time window, the figure carries Plotly + `frames` (one per time step) plus a play/pause control and a slider: each + route's trail grows and an aircraft marker glides along it as the cursor + sweeps the planning window. The risk overlay is static (trace 0); each route + contributes a trail trace then a head trace, so frames update indices 1..2N + and leave the overlay untouched. Rendered on Plotly's `geo` subplot (SVG / Natural Earth vectors) rather than a WebGL tile map, so it draws real country borders + coastlines and displays - even where WebGL is blocked (locked-down / headless Chromium). Returns a - Plotly Figure framed on the data's bounding box. + even where WebGL is blocked (locked-down / headless Chromium). """ import plotly.graph_objects as go # lazy: keep module import lean @@ -193,23 +278,26 @@ def build_map_figure( fig = go.Figure() - # ISSR risk as an Inferno marker cloud over the basemap. `geo` has no native - # density trace, so we draw the grid samples as markers and skip near-zero - # "clear sky" cells, keeping the map readable where there's no risk. + # --- trace 0: ISSR risk as a graded marker cloud over the basemap. `geo` + # has no native density trace, so we draw the grid samples as markers, skip + # near-zero "clear sky" cells, and scale each marker by its risk so hotter + # air reads as bigger + brighter (a softer, denser look than flat dots). rmax = float(np.nanmax(risk_f)) if risk_f.size else 0.0 keep = risk_f > (_MAP_RISK_FLOOR_FRAC * rmax) if rmax > 0 else np.zeros_like(risk_f, bool) + kr = risk_f[keep] + sizes = _MAP_MARKER_SIZE * (0.45 + 1.1 * (kr / rmax)) if rmax > 0 else _MAP_MARKER_SIZE fig.add_trace( go.Scattergeo( lon=lon_f[keep], lat=lat_f[keep], mode="markers", marker=dict( - size=_MAP_MARKER_SIZE, - color=risk_f[keep], + size=sizes, + color=kr, colorscale="Inferno", cmin=0.0, cmax=rmax if rmax > 0 else 1.0, - opacity=0.55, + opacity=0.6, line=dict(width=0), colorbar=dict(title="ISSR
risk", thickness=12, len=0.6, x=0.99), ), @@ -218,20 +306,68 @@ def build_map_figure( ) ) - # Routes: muted gray for context/baseline, bright green for the chosen plan. + # Decide whether to animate: need routes with a non-degenerate time span. + times = None + if animate and routes: + t0 = min(float(np.asarray(r.t_s, dtype=float)[0]) for r in routes) + t1 = max(float(np.asarray(r.t_s, dtype=float)[-1]) for r in routes) + if t1 > t0: + times = np.linspace(t0, t1, _MAP_ANIM_FRAMES) + + # --- traces 1..2N: per route, a trail line then a head marker. The base + # figure shows the first frame (t0) when animating, else the full route. for r in routes: - fig.add_trace( - go.Scattergeo( - lon=np.asarray(r.lon, dtype=float), - lat=np.asarray(r.lat, dtype=float), - mode="lines", - line=dict( - width=3 if r.chosen else 1.5, - color="#28dc5a" if r.chosen else "rgba(200,200,210,0.6)", - ), - name=f"{r.name} (chosen)" if r.chosen else r.name, - hovertemplate=f"{r.name}", - ) + if times is not None: + tl_lon, tl_lat, h_lon, h_lat = _route_state_at(r, float(times[0])) + else: + tl_lon = np.asarray(r.lon, dtype=float) + tl_lat = np.asarray(r.lat, dtype=float) + h_lon, h_lat = float(tl_lon[-1]), float(tl_lat[-1]) + fig.add_trace(_trail_trace(go, r, tl_lon, tl_lat)) + fig.add_trace(_head_trace(go, r, h_lon, h_lat)) + + # --- frames + play/pause + slider (features: route reveal AND time-step) --- + if times is not None: + anim_indices = list(range(1, 1 + 2 * len(routes))) + frames = [] + for k, tc in enumerate(times): + data = [] + for r in routes: + tl_lon, tl_lat, h_lon, h_lat = _route_state_at(r, float(tc)) + data.append(_trail_trace(go, r, tl_lon, tl_lat)) + data.append(_head_trace(go, r, h_lon, h_lat)) + frames.append(go.Frame(name=str(k), data=data, traces=anim_indices)) + fig.frames = frames + + play = dict( + label="▶ Play", method="animate", + args=[None, dict(frame=dict(duration=_MAP_FRAME_MS, redraw=True), + fromcurrent=True, transition=dict(duration=0), + mode="immediate")], + ) + pause = dict( + label="⏸ Pause", method="animate", + args=[[None], dict(frame=dict(duration=0, redraw=False), + mode="immediate", transition=dict(duration=0))], + ) + fig.update_layout( + updatemenus=[dict( + type="buttons", direction="left", showactive=False, + x=0.01, y=0.02, xanchor="left", yanchor="bottom", + bgcolor="rgba(20,20,24,0.6)", font=dict(color="#dddde2"), + buttons=[play, pause], + )], + sliders=[dict( + active=0, x=0.12, len=0.8, y=0.02, yanchor="bottom", + bgcolor="rgba(20,20,24,0.6)", font=dict(color="#bdbdc6", size=10), + currentvalue=dict(prefix="t = ", suffix=" min", font=dict(color="#dddde2")), + steps=[dict( + method="animate", label=f"{tc / 60:.0f}", + args=[[str(k)], dict(mode="immediate", + frame=dict(duration=0, redraw=True), + transition=dict(duration=0))], + ) for k, tc in enumerate(times)], + )], ) # Frame the view on the data box with a small margin; dark land/ocean makes @@ -703,18 +839,13 @@ def _build_map_tab(self) -> QWidget: return self.map_view def _active_anchor(self, field): - """The GeoAnchor placing this ISSR field on the map. + """The GeoAnchor placing the synthetic ISSR field on the map. - The MLIssrField carries its own anchor; the synthetic field has none, - so we place it with the SAME canonical anchor the model/features use, - so synthetic and ml render over the same European box. + The field carries no anchor of its own (geography only matters for the + map), so we place it with the canonical European anchor that matches the + default world geometry. """ - anchor = getattr(field, "anchor", None) - if anchor is not None: - return anchor - from contrail_ml.config import DEFAULT_CONFIG - - return DEFAULT_CONFIG.anchor + return getattr(field, "anchor", None) or EUROPEAN_ANCHOR def _map_routes(self, chosen_by_flight: dict[str, int] | None, anchor) -> list[RouteLine]: """Per-flight ground tracks in lon/lat (chosen profile, else baseline).""" @@ -731,12 +862,14 @@ def _map_routes(self, chosen_by_flight: dict[str, int] | None, anchor) -> list[R wps = waypoints_for(flight, profile) wx = np.array([w[0] for w in wps], dtype=float) wy = np.array([w[1] for w in wps], dtype=float) + wt = np.array([w[3] for w in wps], dtype=float) # world time per vertex lon, lat = anchor.local_to_geo(wx, wy) routes.append( RouteLine( name=flight.name, lon=np.asarray(lon, dtype=float), lat=np.asarray(lat, dtype=float), + t_s=wt, chosen=is_chosen, ) ) @@ -745,8 +878,7 @@ def _map_routes(self, chosen_by_flight: dict[str, int] | None, anchor) -> list[R def _render_map(self, chosen_by_flight: dict[str, int] | None) -> None: """Render the ISSR-risk overlay + flight routes onto the geographic map. - Works for BOTH the synthetic ISSRField and the MLIssrField, since they - share the rhi_excess_grid interface. On a 2-D map the altitude dimension + Uses the ISSRField's rhi_excess_grid interface. On a 2-D map the altitude dimension collapses, so each flight is one ground track; once solved, chosen routes are drawn green over the muted context. Builds a Plotly figure and loads its self-contained HTML into the embedded web view. @@ -779,6 +911,8 @@ def _render_map(self, chosen_by_flight: dict[str, int] | None) -> None: assets = _geo_assets_script() if assets: html = html.replace("", assets + "", 1) + # Auto-run the reveal once the plot is live (replays on every solve). + html = html.replace("", _autoplay_script() + "", 1) with open(self._map_html_path, "w", encoding="utf-8") as fh: fh.write(html) if self.map_view is not None: diff --git a/pyproject.toml b/pyproject.toml index e94045b..1d7fa07 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -41,39 +41,6 @@ dev = [ "mypy>=1.8", "ruff>=0.4", ] -# Observation-grounded ISSR model + full MLOps lifecycle. Install with: -# pip install -e ".[ml]" -# Deliberately NOT in the core deps: the deployed solver image (Dockerfile -# runtime stage) must stay lean and free of the heavy scientific/geo stack. -# The geo loaders (pycontrails/xarray/cfgrib/zarr/gcsfs/netCDF4/dask) are only -# needed to BUILD a dataset from real weather; training/serving/tests on a -# prepared parquet need only the tabular stack (pandas/scipy/sklearn/xgboost). -ml = [ - "pycontrails>=0.60", - "xarray>=2024.0", - "netCDF4", - "cfgrib", # GFS is GRIB - "zarr", - "gcsfs", # ARCO-ERA5 is zarr on a public GCS bucket - "dask", - "pandas>=2.0", - "pyarrow", - "scipy", - "scikit-learn>=1.4", - "xgboost>=2.0", - "mlflow>=2.12", - "matplotlib", # reliability / diagnostic plots logged to MLflow - "tabulate", # markdown rendering of the comparison table for MLflow -] -# Real historical-traffic loader (OpenSky). The package's geometry + reduction -# run on numpy alone; this extra is only needed to actually PULL and CACHE real -# tracks (pyopensky + the tabular stack). Like [ml], it is NOT in the core deps, -# so the deployed solver image stays lean. -flights = [ - "pyopensky>=2.0", - "pandas>=2.0", - "pyarrow", -] # Real quantum SDK backends. The pipelines run WITHOUT these (built-in # statevector / exact-hafnian samplers); installing them switches # pasqal_analog to Pulser's QuTiP emulator and xanadu_gbs to Strawberry @@ -94,10 +61,7 @@ Repository = "https://github.com/jaewonyun1234/Flight_path_optimization_Contrail # the `service` directory so the generated subpackage imports fine once built. # --------------------------------------------------------------------------- [tool.setuptools.packages.find] -# Auto-discover packages present in the build context. Locally this finds -# contrail_env, service, contrail_ml, and contrail_ml.data. Inside the Docker -# build context, contrail_ml is excluded via .dockerignore so only -# contrail_env and service are found — pip install . succeeds either way. +# Auto-discover packages present in the build context: contrail_env and service. where = ["."] # --------------------------------------------------------------------------- @@ -109,6 +73,7 @@ where = ["."] [tool.ruff] line-length = 100 extend-exclude = [ + "playground", # scratch / sandbox notebooks — not linted "service/generated", "contrail_env/__init__.py", "contrail_env/demo.py", @@ -126,12 +91,6 @@ extend-exclude = [ select = ["E", "F", "I", "W", "UP", "B"] ignore = ["E501"] -# The ML test modules call pytest.importorskip(...) before importing contrail_ml -# so the core lint-type-test job (no [ml] extra) skips them instead of erroring. -# That guard must run before the imports it protects, which is exactly E402. -[tool.ruff.lint.per-file-ignores] -"tests/test_ml_*.py" = ["E402"] - # --------------------------------------------------------------------------- # Mypy — fully type-check the new code. The pre-existing modules cannot be # edited, so their (best-effort) type errors are not enforced; the new diff --git a/service/proto/solver.proto b/service/proto/solver.proto index 5ac06a1..0eb55ab 100644 --- a/service/proto/solver.proto +++ b/service/proto/solver.proto @@ -16,17 +16,14 @@ message ScenarioConfig { double time_limit_s = 9; // CP-SAT time limit string progress_topic = 10; // ZMQ topic to publish progress on double issr_threshold = 11; // RHi-excess above which a cell forms a contrail - // ML ISSR source (plan §7). Default "" == "synthetic", so existing clients - // are unaffected. "ml" builds the world's ISSR from a trained model. - string issr_source = 12; // "synthetic" (default) or "ml" - string issr_time = 13; // ISO valid-time for the ML forecast (optional) - double issr_p_threshold = 14; // calibrated P(ISSR) cut for is_inside (0 -> default) - // Real historical-traffic source (contrail_flights / OpenSky). Default - // "" == "synthetic", so existing clients are unaffected. "real" pulls actual - // flown ADS-B tracks for the window below and reduces them to flights. - string flight_source = 15; // "synthetic" (default) or "real" - string flight_start_time = 16; // ISO start of the historical pull window - string flight_end_time = 17; // ISO end of the historical pull window + // ISSR source. Only "synthetic" (default; "" == "synthetic") is supported. + string issr_source = 12; // "synthetic" (default) + string issr_time = 13; // reserved (unused) + double issr_p_threshold = 14; // reserved (unused) + // Flight source. Only "synthetic" (default; "" == "synthetic") is supported. + string flight_source = 15; // "synthetic" (default) + string flight_start_time = 16; // reserved (unused) + string flight_end_time = 17; // reserved (unused) } message FlightChoice { diff --git a/service/scenario.py b/service/scenario.py index 34859bf..420d99b 100644 --- a/service/scenario.py +++ b/service/scenario.py @@ -59,46 +59,18 @@ def build_world_and_flights(cfg: solver_pb2.ScenarioConfig) -> tuple[World, list * its own baseline cruise flight level. """ source = cfg.issr_source or "synthetic" - if source == "synthetic": - world = default_european_world(seed=cfg.seed, n_issr_blobs=cfg.n_issr_blobs) - # Threshold drives what counts as a contrail (cell RHi-excess > threshold), - # so it affects the solve, not just the picture. 0 means "use the default". - if cfg.issr_threshold > 0: - world.issr.threshold = float(cfg.issr_threshold) - elif source == "ml": - try: - from contrail_ml.config import MLConfig - except ImportError as exc: - raise ValueError( - "issr_source='ml' is not available in this deployment — " - "the server image ships only the CP-SAT solver. " - "Install contrail_ml locally with: pip install -e '.[ml]'" - ) from exc - - ml_cfg = MLConfig() - if cfg.issr_p_threshold > 0: - ml_cfg = ml_cfg.with_overrides(issr_p_threshold=cfg.issr_p_threshold) - world = default_european_world( - seed=cfg.seed, - issr_source="ml", - issr_kwargs=dict( - config=ml_cfg, - met_source="gfs", # operational forecast - time=cfg.issr_time or None, - ), - ) - else: - raise ValueError(f"unknown issr_source {source!r} (use 'synthetic' or 'ml')") + if source != "synthetic": + raise ValueError(f"unknown issr_source {source!r} (use 'synthetic')") + world = default_european_world(seed=cfg.seed, n_issr_blobs=cfg.n_issr_blobs) + # Threshold drives what counts as a contrail (cell RHi-excess > threshold), + # so it affects the solve, not just the picture. 0 means "use the default". + if cfg.issr_threshold > 0: + world.issr.threshold = float(cfg.issr_threshold) - # Flight source: synthetic random generator (default) or real historical - # OpenSky traffic. The real path keeps its own origin/destination/baseline, - # so it returns before the synthetic chord-geometry randomization below. flight_source = cfg.flight_source or "synthetic" - if flight_source == "real": - return world, _build_real_flights(cfg, world) if flight_source != "synthetic": raise ValueError( - f"unknown flight_source {flight_source!r} (use 'synthetic' or 'real')" + f"unknown flight_source {flight_source!r} (use 'synthetic')" ) flights = build_random_flights( @@ -149,49 +121,6 @@ def _clamp(v: float, lo: float, hi: float) -> float: return world, flights -def _build_real_flights(cfg: solver_pb2.ScenarioConfig, world: World) -> list[Flight]: - """Build real historical flights via contrail_flights (OpenSky). - - contrail_flights is imported lazily (and only here), so the core server - image without the [flights] extra is unaffected unless a client explicitly - asks for flight_source="real". Missing deps / credentials / network surface - as a clear ValueError (mapped to a gRPC error by the Solve handler), never a - silent crash or fabricated traffic. - """ - try: - from contrail_flights.build_dataset import flights_for - from contrail_flights.config import FlightsConfig - from contrail_flights.opensky_client import OpenSkyUnavailableError - except ImportError as exc: - raise ValueError( - "flight_source='real' is not available in this deployment — the " - "server image ships only the synthetic flight generator. Install " - "contrail_flights locally with: pip install -e '.[flights]'" - ) from exc - - # Derive the geo planning box from the world grid via the shared anchor, so - # the real flights land in the same local frame the optimizer/grid use. - fcfg = FlightsConfig() - g = world.grid - lon0, lat0 = fcfg.anchor.local_to_geo(g.x_min_km, g.y_min_km) - lon1, lat1 = fcfg.anchor.local_to_geo(g.x_max_km, g.y_max_km) - overrides: dict[str, object] = dict( - lon_min=min(lon0, lon1), lon_max=max(lon0, lon1), - lat_min=min(lat0, lat1), lat_max=max(lat0, lat1), - snapshot_window_s=cfg.snapshot_window_s or 300.0, - ) - if cfg.flight_start_time: - overrides["start_time"] = cfg.flight_start_time - if cfg.flight_end_time: - overrides["end_time"] = cfg.flight_end_time - fcfg = fcfg.with_overrides(**overrides) - - try: - return flights_for(fcfg, max_flights=cfg.n_flights or None) - except OpenSkyUnavailableError as exc: - raise ValueError(str(exc)) from exc - - def build_scenario_full( cfg: solver_pb2.ScenarioConfig, ) -> tuple[World, list[Flight], list[EvaluatedOption], list[ConflictEdge], list[CapacityBucket]]: diff --git a/tests/test_flight_reduction.py b/tests/test_flight_reduction.py deleted file mode 100644 index 39dcf9f..0000000 --- a/tests/test_flight_reduction.py +++ /dev/null @@ -1,77 +0,0 @@ -"""reduce_to_flight: a fake ADS-B track -> a valid contrail_env.Flight. - -Hermetic — builds Track fixtures by hand, no network/credentials. Checks the -core conversion: origin/destination inside the planning box, eastbound geometry, -the observed climb preserved as the baseline, and the loud failures (no in-box -points; all-NaN altitude falls back instead of fabricating). -""" - -import numpy as np -import pytest - -from contrail_env import fl_to_m -from contrail_flights.config import FlightsConfig -from contrail_flights.reduce_to_flight import aircraft_for, reduce_to_flight -from contrail_flights.tracks import make_track - - -def _climbing_track(n: int = 20): - """West->east crossing of the box at ~46.5N, climbing FL340 -> FL360.""" - t = np.arange(n) * 60.0 - lon = np.linspace(-4.0, 13.0, n) - lat = np.full(n, 46.5) - alt = np.where(np.arange(n) < n // 2, fl_to_m(340), fl_to_m(360)) - return make_track("abc123", "TEST123 ", t, lat, lon, alt) - - -def test_reduce_produces_valid_eastbound_flight(): - cfg = FlightsConfig() - reduced = reduce_to_flight(_climbing_track(), cfg, name="R1") - f = reduced.flight - - assert reduced.used_observed_baseline - assert 0.0 <= f.origin_km[0] <= 1500.0 - assert 0.0 <= f.origin_km[1] <= 800.0 - assert f.destination_km[0] > f.origin_km[0] # eastbound - assert f.departure_s == 0.0 # first (and only) track -> reference epoch - - -def test_observed_climb_becomes_the_baseline(): - cfg = FlightsConfig() - reduced = reduce_to_flight(_climbing_track(), cfg) - fls = [s.fl for s in reduced.flight.baseline.segments] - assert fls == [340, 360] - segs = reduced.flight.baseline.segments - assert segs[0].t_start_s == 0.0 - # Contiguous: each segment's end is the next one's start. - assert segs[0].t_end_s == segs[1].t_start_s - - -def test_track_with_no_in_box_points_raises(): - cfg = FlightsConfig() - n = 12 - t = np.arange(n) * 60.0 - lon = np.linspace(-4.0, 13.0, n) - lat = np.full(n, 10.0) # far south of the box -> never inside - alt = np.full(n, fl_to_m(360)) - with pytest.raises(ValueError): - reduce_to_flight(make_track("x", "Y", t, lat, lon, alt), cfg) - - -def test_all_nan_altitude_falls_back_not_fabricates(): - cfg = FlightsConfig() - n = 12 - t = np.arange(n) * 60.0 - lon = np.linspace(-4.0, 13.0, n) - lat = np.full(n, 46.5) - alt = np.full(n, np.nan) - reduced = reduce_to_flight(make_track("x", "Y", t, lat, lon, alt), cfg) - assert not reduced.used_observed_baseline - assert reduced.flight.baseline.n_segments >= 1 - - -def test_aircraft_type_mapping_defaults_to_a320(): - a = aircraft_for("UNMAPPED", "Z999") - assert a.tail == "UNMAPPED" - # A mapped type and an unmapped one both resolve to a valid Aircraft. - assert aircraft_for("X", "A320").initial_mass_kg > 0 diff --git a/tests/test_geo_transform_roundtrip.py b/tests/test_geo_transform_roundtrip.py index 87424d8..05f48e6 100644 --- a/tests/test_geo_transform_roundtrip.py +++ b/tests/test_geo_transform_roundtrip.py @@ -1,13 +1,13 @@ -"""The flight loader's geo transform inverts cleanly and matches the model's. +"""The map's geo transform inverts cleanly: geo_to_local(local_to_geo(x,y)) == (x,y). -geo_to_local(local_to_geo(x, y)) == (x, y), and the (deliberately duplicated) -contrail_flights anchor agrees with contrail_ml's so the map and the real -flights share one coordinate system. +The Map tab places the local sim box on a lon/lat basemap via contrail_env.geo; +this guards that the forward/inverse pair is consistent so routes and the risk +overlay land where they should. """ import numpy as np -from contrail_flights.geo import GeoAnchor +from contrail_env.geo import EUROPEAN_ANCHOR, GeoAnchor def test_roundtrip_vectorized(): @@ -29,16 +29,8 @@ def test_roundtrip_scalar_returns_floats(): assert abs(x - 750.0) < 1e-6 and abs(y - 400.0) < 1e-6 -def test_matches_contrail_ml_anchor(): - # The independent copy MUST agree with the model's anchor (same origin), - # or the predicted field and the real flights would land in different places. - from contrail_ml.config import MLConfig - from contrail_ml.features import GeoAnchor as MLAnchor - - cfg = MLConfig() - flights_anchor = GeoAnchor() - ml_anchor = MLAnchor(origin_lat=cfg.origin_lat, origin_lon=cfg.origin_lon) - lon_f, lat_f = flights_anchor.local_to_geo(500.0, 300.0) - lon_m, lat_m = ml_anchor.local_to_geo(500.0, 300.0) - assert abs(lon_f - lon_m) < 1e-9 - assert abs(lat_f - lat_m) < 1e-9 +def test_canonical_anchor_sits_over_europe(): + # The default box anchor should place (x=0, y=0) at south-west Europe. + lon, lat = EUROPEAN_ANCHOR.local_to_geo(0.0, 0.0) + assert abs(lon - (-5.0)) < 1e-9 + assert abs(lat - 43.0) < 1e-9 diff --git a/tests/test_gui_map_panel.py b/tests/test_gui_map_panel.py index c1251a2..5e61203 100644 --- a/tests/test_gui_map_panel.py +++ b/tests/test_gui_map_panel.py @@ -26,7 +26,18 @@ from PyQt6.QtWidgets import QApplication # noqa: E402 from contrail_env import build_random_flights, default_european_world # noqa: E402 -from gui.app import MainWindow, RouteLine, _geo_assets_script, build_map_figure # noqa: E402 +from gui.app import ( # noqa: E402 + _MAP_ANIM_FRAMES, + MainWindow, + RouteLine, + _geo_assets_script, + build_map_figure, +) + + +def _route(name, lon, lat, t_s, chosen): + return RouteLine(name, np.asarray(lon, float), np.asarray(lat, float), + np.asarray(t_s, float), chosen) @pytest.fixture(scope="module") @@ -40,32 +51,54 @@ def qapp(): # Pure figure builder — no Qt, no web view. # # --------------------------------------------------------------------------- # -def test_build_map_figure_has_risk_and_route_traces(): +def test_build_map_figure_animates_routes_with_frames(): lon, lat = np.meshgrid(np.linspace(-2, 18, 20), np.linspace(43, 50, 15), indexing="ij") risk = np.abs(np.sin(lon) * np.cos(lat)) + t = np.array([0.0, 400.0, 800.0, 1200.0]) routes = [ - RouteLine("AB123", np.array([0.0, 5.0, 10.0]), np.array([45.0, 46.0, 47.0]), chosen=True), - RouteLine("CD456", np.array([1.0, 6.0]), np.array([44.0, 48.0]), chosen=False), + _route("AB123", [0, 5, 8, 10], [45, 46, 46.5, 47], t, chosen=True), + _route("CD456", [1, 3, 5, 6], [44, 45, 47, 48], t, chosen=False), ] fig = build_map_figure(source="synthetic", lon=lon, lat=lat, risk=risk, routes=routes) - # SVG `geo` traces (no WebGL): one marker overlay for risk + one line per route. - assert all(t.type == "scattergeo" for t in fig.data) - modes = [t.mode for t in fig.data] - assert modes.count("markers") == 1 # the ISSR risk overlay - assert modes.count("lines") == 2 # one trace per route - # The figure must serialise to a self-contained page (this is what the - # web view loads); inlined plotly.js makes it well over the setHtml limit. + # SVG `geo` traces (no WebGL). trace 0 = risk overlay; then per route a + # trail line + a head marker, so markers = risk + 2 heads, lines = 2 trails. + assert all(tr.type == "scattergeo" for tr in fig.data) + assert fig.data[0].mode == "markers" + modes = [tr.mode for tr in fig.data] + assert modes.count("lines") == 2 + assert modes.count("markers") == 3 + # Animation machinery: one frame per time step, a play/pause control + slider. + assert len(fig.frames) == _MAP_ANIM_FRAMES + assert fig.layout.updatemenus and fig.layout.sliders + # Frames update only the 2N route traces, never the risk overlay (trace 0). + assert list(fig.frames[0].traces) == [1, 2, 3, 4] + # Self-contained page (what the web view loads); inlined plotly.js is large. html = fig.to_html(include_plotlyjs="inline", full_html=True) assert " 2_000_000 +def test_build_map_figure_static_when_not_animated(): + lon, lat = np.meshgrid(np.linspace(-2, 18, 10), np.linspace(43, 50, 8), indexing="ij") + risk = np.abs(np.sin(lon)) + t = np.array([0.0, 600.0, 1200.0]) + routes = [_route("AB123", [0, 5, 10], [45, 46, 47], t, chosen=True)] + fig = build_map_figure(source="synthetic", lon=lon, lat=lat, risk=risk, + routes=routes, animate=False) + # No frames; the route is drawn in full (one trail line + one head marker). + assert not fig.frames + modes = [tr.mode for tr in fig.data] + assert modes.count("lines") == 1 + assert modes.count("markers") == 2 # risk overlay + head marker + + def test_build_map_figure_handles_empty_routes(): lon, lat = np.meshgrid(np.linspace(0, 10, 8), np.linspace(44, 49, 6), indexing="ij") fig = build_map_figure(source="synthetic", lon=lon, lat=lat, risk=np.zeros_like(lon), routes=[]) - # Only the (empty, all-clear-sky) risk overlay; no routes. + # Only the (empty, all-clear-sky) risk overlay; no routes, no frames. assert [t.mode for t in fig.data] == ["markers"] + assert not fig.frames def test_geo_basemap_is_bundled_offline(): @@ -98,33 +131,17 @@ def test_render_map_synthetic_writes_html(qapp): assert os.path.getsize(win._map_html_path) > 2_000_000 -def test_render_map_ml_writes_html(qapp): - # The ML field needs the [ml] extra (scipy interpolation + the model). - pytest.importorskip("scipy") - pytest.importorskip("sklearn") - import warnings - - from contrail_ml.config import MLConfig - from contrail_ml.issr_field import MLIssrField - - cfg = MLConfig() - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - world = default_european_world( - seed=1, - issr_source="ml", - issr_kwargs=dict(config=cfg, met_source="synthetic", - allow_synthetic=True, grid_res_deg=3.0, seed=1), - ) - assert isinstance(world.issr, MLIssrField) +def test_render_map_uses_canonical_anchor(qapp): + # The synthetic field carries no anchor of its own, so the map falls back to + # the canonical European anchor. + from contrail_env.geo import EUROPEAN_ANCHOR + + world = default_european_world(seed=1) flights = build_random_flights(n_flights=2, world=world, seed=1, corridor_frac=0.05, snapshot_window_s=(0.0, 300.0)) - win = MainWindow() win._world = world win._flights = flights - # Renders the ML risk field + routes without raising; the active anchor - # comes from the field itself (MLIssrField.anchor), not the synthetic default. win._render_map(None) - assert win._active_anchor(world.issr) is world.issr.anchor + assert win._active_anchor(world.issr) is EUROPEAN_ANCHOR assert os.path.exists(win._map_html_path) diff --git a/tests/test_ml_calibrate.py b/tests/test_ml_calibrate.py deleted file mode 100644 index 671dfb8..0000000 --- a/tests/test_ml_calibrate.py +++ /dev/null @@ -1,48 +0,0 @@ -"""Calibration: isotonic lowers ECE; conformal interval covers ~(1-alpha).""" - -import numpy as np -import pytest - -# Needs the [ml] extra (scikit-learn). The core lint-type-test CI job installs -# only [dev], so skip there; the dedicated `ml` job installs [ml] and runs these. -pytest.importorskip("sklearn") - -from contrail_ml.calibrate import ( - ConformalRHi, - ProbabilityCalibrator, - expected_calibration_error, -) - - -def test_isotonic_lowers_ece_on_heldout(): - rng = np.random.default_rng(0) - n = 4000 - # Latent score -> true label; observed prob is a miscalibrated (squashed) - # version of the latent, so calibration has real work to do. - latent = rng.normal(size=n) - y = (latent + rng.normal(scale=0.5, size=n) > 0).astype(int) - p_true = 1.0 / (1.0 + np.exp(-latent)) - p_miscal = p_true ** 2 # systematically under-confident - - half = n // 2 - cal = ProbabilityCalibrator(method="isotonic").fit(p_miscal[:half], y[:half]) - p_cal = cal.transform(p_miscal[half:]) - - ece_before = expected_calibration_error(p_miscal[half:], y[half:]) - ece_after = expected_calibration_error(p_cal, y[half:]) - assert ece_after < ece_before - assert np.all((p_cal >= 0.0) & (p_cal <= 1.0)) - - -def test_conformal_interval_covers_target(): - rng = np.random.default_rng(1) - n = 5000 - rhi_hat = rng.uniform(40, 140, size=n) - rhi_true = rhi_hat + rng.normal(scale=6.0, size=n) - - half = n // 2 - conf = ConformalRHi(alpha=0.1).fit(rhi_true[:half], rhi_hat[:half]) - lo, hi = conf.interval(rhi_hat[half:]) - covered = np.mean((rhi_true[half:] >= lo) & (rhi_true[half:] <= hi)) - assert conf.half_width > 0 - assert covered >= 0.85 # ~90% target, allow sampling slack diff --git a/tests/test_ml_features.py b/tests/test_ml_features.py deleted file mode 100644 index 7b95462..0000000 --- a/tests/test_ml_features.py +++ /dev/null @@ -1,67 +0,0 @@ -"""Shared feature engineering: thermodynamics, atmosphere, encodings, geo. - -These are the anti-skew keystone (one definition for train and serve), so they -get known-value and round-trip checks. -""" - -import numpy as np -import pytest - -from contrail_ml.features import ( - GeoAnchor, - altitude_to_pressure_hpa, - day_of_year_encoding, - e_si_pa, - pressure_to_altitude_m, - rhi_excess, - rhi_percent, - solar_time_encoding, -) - - -def test_e_si_monotonic_and_positive(): - temps = np.array([210.0, 220.0, 230.0, 240.0, 250.0]) - esi = np.asarray(e_si_pa(temps)) - assert np.all(esi > 0) - assert np.all(np.diff(esi) > 0) # warmer air holds more vapour over ice - - -def test_rhi_hits_100_at_ice_saturation(): - # If vapour pressure equals the ice-saturation pressure, RHi must be 100%. - T, p_hpa = 225.0, 225.0 - esi = float(e_si_pa(T)) # Pa - e = esi # exactly saturated - p_pa = p_hpa * 100.0 - eps = 0.621981 - q = eps * e / (p_pa - (1.0 - eps) * e) # invert vapour-pressure formula - assert abs(float(rhi_percent(T, q, p_hpa)) - 100.0) < 1e-6 - - -def test_rhi_excess_semantics(): - assert float(rhi_excess(130.0)) == pytest.approx(0.3) - assert float(rhi_excess(80.0)) == 0.0 - assert float(rhi_excess(100.0)) == 0.0 - - -def test_altitude_pressure_roundtrip(): - for h in [9000.0, 10363.0, 10973.0, 11582.0, 12192.0]: - p = float(altitude_to_pressure_hpa(h)) - h2 = float(pressure_to_altitude_m(p)) - assert abs(h - h2) < 1.0 # within a metre - - -def test_cyclical_encodings_unit_circle(): - s, c = solar_time_encoding(np.array([0.0, 6.0, 18.0]), np.array([0.0, 0.0, 0.0])) - assert np.allclose(np.asarray(s) ** 2 + np.asarray(c) ** 2, 1.0) - sd, cd = day_of_year_encoding(np.array([1.0, 90.0, 180.0, 365.0])) - assert np.allclose(np.asarray(sd) ** 2 + np.asarray(cd) ** 2, 1.0) - - -def test_geo_anchor_roundtrip(): - anchor = GeoAnchor(origin_lat=43.0, origin_lon=-5.0) - xs = np.array([0.0, 750.0, 1500.0]) - ys = np.array([0.0, 400.0, 800.0]) - lon, lat = anchor.local_to_geo(xs, ys) - x2, y2 = anchor.geo_to_local(lon, lat) - assert np.allclose(np.asarray(x2), xs, atol=1e-6) - assert np.allclose(np.asarray(y2), ys, atol=1e-6) diff --git a/tests/test_ml_guards.py b/tests/test_ml_guards.py deleted file mode 100644 index 727db69..0000000 --- a/tests/test_ml_guards.py +++ /dev/null @@ -1,44 +0,0 @@ -"""Honesty guards (plan §0.1): synthetic data is opt-in and loud; real loaders -fail clearly when data is absent rather than fabricating it.""" - -import warnings - -import pytest - -from contrail_ml.config import MLConfig -from contrail_ml.data.synthetic_fallback import make_synthetic_training_table - - -def test_synthetic_fallback_refuses_without_flag(): - with pytest.raises(RuntimeError, match="allow_synthetic"): - make_synthetic_training_table(n=100) - - -def test_synthetic_fallback_warns_loudly_when_allowed(): - with pytest.warns(UserWarning, match="SYNTHETIC"): - df = make_synthetic_training_table(n=200, allow_synthetic=True) - assert (df["source"] == "synthetic_fallback").all() - - -def test_synthetic_met_cube_refuses_without_flag(): - from contrail_ml import predict - - with pytest.raises(RuntimeError, match="allow_synthetic"): - predict.synthetic_met_cube(MLConfig()) - - -def test_synthetic_ml_issr_field_refuses_without_flag(): - from contrail_ml.issr_field import ml_issr_field - - with pytest.raises(RuntimeError, match="allow_synthetic"): - ml_issr_field(MLConfig(), met_source="synthetic") - - -def test_iagos_loader_raises_when_dir_missing(): - from contrail_ml.data.iagos import load_iagos_waypoints - - cfg = MLConfig(iagos_dir=None) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - with pytest.raises(FileNotFoundError): - load_iagos_waypoints(cfg) diff --git a/tests/test_ml_issr_field.py b/tests/test_ml_issr_field.py deleted file mode 100644 index ef16682..0000000 --- a/tests/test_ml_issr_field.py +++ /dev/null @@ -1,71 +0,0 @@ -"""MLIssrField conforms to the ISSRField interface, with correct geo mapping.""" - -import numpy as np -import pytest - -# Needs the [ml] extra (scipy interpolation). The core lint-type-test CI job -# installs only [dev], so skip there; the dedicated `ml` job runs these fully. -pytest.importorskip("scipy") - -from contrail_env import fl_to_m -from contrail_ml.features import GeoAnchor, altitude_to_pressure_hpa -from contrail_ml.issr_field import MLIssrField - - -def _field_with_blob_at(local_x, local_y, fl): - """Build an MLIssrField whose ISSR blob is centred on the geo point that - the given local (x_km, y_km, FL) maps to.""" - anchor = GeoAnchor(origin_lat=43.0, origin_lon=-5.0) - lon0, lat0 = anchor.local_to_geo(local_x, local_y) - p0 = float(altitude_to_pressure_hpa(fl_to_m(fl))) - - lon_axis = np.arange(-30.0, 30.01, 2.0) - lat_axis = np.arange(30.0, 70.01, 2.0) - pressure_axis = np.array([150, 175, 200, 225, 250, 300], dtype=float) - lon3, lat3, p3 = np.meshgrid(lon_axis, lat_axis, pressure_axis, indexing="ij") - - blob = np.exp(-(((lon3 - lon0) / 4.0) ** 2 - + ((lat3 - lat0) / 3.0) ** 2 - + ((p3 - p0) / 25.0) ** 2)) - return MLIssrField( - lon_axis=lon_axis, lat_axis=lat_axis, pressure_axis=pressure_axis, - rhi_excess_cube=0.4 * blob, p_issr_cube=blob, anchor=anchor, - p_threshold=0.5, mode="prob", - ) - - -def test_implements_issrfield_interface(): - field = _field_with_blob_at(750.0, 400.0, 360) - for attr in ("rhi_excess", "is_inside", "rhi_excess_grid", "mask_grid", "threshold"): - assert hasattr(field, attr) - - -def test_point_inside_known_blob(): - field = _field_with_blob_at(750.0, 400.0, 360) - z = fl_to_m(360) - # The blob centre maps back to this local point -> inside. - assert field.is_inside(750.0, 400.0, z) is True - # A point near the opposite corner of the box -> outside the blob. - assert field.is_inside(50.0, 50.0, z) is False - - -def test_grid_query_shapes(): - field = _field_with_blob_at(750.0, 400.0, 360) - xs = np.linspace(50, 1450, 20) - ys = np.linspace(50, 750, 8) - xx, yy = np.meshgrid(xs, ys, indexing="ij") - zz = np.full_like(xx, fl_to_m(360)) - mask = field.mask_grid(xx, yy, zz) - exc = field.rhi_excess_grid(xx, yy, zz) - assert mask.shape == xx.shape - assert exc.shape == xx.shape - assert mask.any() # the blob shows up somewhere on the slice - - -def test_rhi_mode_uses_threshold(): - field = _field_with_blob_at(750.0, 400.0, 360) - field.mode = "rhi" - field.threshold = 0.3 - z = fl_to_m(360) - # rhi_excess peaks at 0.4 at the centre (> 0.3) -> inside in rhi mode. - assert field.is_inside(750.0, 400.0, z) is True diff --git a/tests/test_ml_model.py b/tests/test_ml_model.py deleted file mode 100644 index e26790e..0000000 --- a/tests/test_ml_model.py +++ /dev/null @@ -1,73 +0,0 @@ -"""RHiCorrector: fits, predicts, routes regimes, and reduces the dry bias. - -All hermetic — trains on the guarded synthetic fallback, never on real data. -""" - -import warnings - -import numpy as np -import pytest - -# Needs the [ml] extra (pandas/scikit-learn/xgboost). The core lint-type-test CI -# job installs only [dev], so skip there; the dedicated `ml` job runs these fully. -pytest.importorskip("pandas") -pytest.importorskip("sklearn") -pytest.importorskip("xgboost") - -from contrail_ml.data.schema import LABEL_DELTA -from contrail_ml.data.synthetic_fallback import make_synthetic_training_table -from contrail_ml.model import RHiCorrector - - -@pytest.fixture(scope="module") -def table(): - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - return make_synthetic_training_table(n=4000, seed=7, allow_synthetic=True) - - -def _small_model(): - return RHiCorrector(n_ensemble=3, xgb_n_estimators=60, min_regime_samples=30, - regime_split_rhi=85.0, random_state=0) - - -def test_fit_predict_shapes_and_ranges(table): - train, test = table.iloc[:3000], table.iloc[3000:] - m = _small_model().fit(train, train[LABEL_DELTA].to_numpy()) - rhi_hat, rhi_std, p_issr = m.predict(test) - - n = len(test) - assert rhi_hat.shape == (n,) and rhi_std.shape == (n,) and p_issr.shape == (n,) - assert np.all(p_issr >= 0.0) and np.all(p_issr <= 1.0) - assert np.all(rhi_std >= 0.0) and rhi_std.max() > 0.0 # ensemble has spread - - -def test_model_reduces_dry_bias(table): - train, test = table.iloc[:3000], table.iloc[3000:] - m = _small_model().fit(train, train[LABEL_DELTA].to_numpy()) - rhi_hat, _s, _p = m.predict(test) - - truth = test["rhi_iagos"].to_numpy() - raw_mae = np.abs(test["rhi"].to_numpy() - truth).mean() - model_mae = np.abs(rhi_hat - truth).mean() - assert model_mae < raw_mae # correcting helps - - -def test_predict_before_fit_raises(table): - with pytest.raises(RuntimeError): - _small_model().predict(table) - - -def test_predict_proba_issr_in_unit_interval(table): - m = _small_model().fit(table, table[LABEL_DELTA].to_numpy()) - p = m.predict_proba_issr(table) - assert np.all((p >= 0.0) & (p <= 1.0)) - - -def test_combined_fallback_when_regime_sparse(table): - # A very high min_regime_samples forces the combined-model branch. - m = RHiCorrector(n_ensemble=2, xgb_n_estimators=40, min_regime_samples=10_000) - m.fit(table, table[LABEL_DELTA].to_numpy()) - rhi_hat, _s, _p = m.predict(table) - assert rhi_hat.shape[0] == len(table) - assert all(member["mode"] == "combined" for member in m._ensemble) diff --git a/tests/test_ml_train.py b/tests/test_ml_train.py deleted file mode 100644 index 2c9b240..0000000 --- a/tests/test_ml_train.py +++ /dev/null @@ -1,45 +0,0 @@ -"""Training pipeline (no MLflow): CV runs, model beats raw ERA5, table is built. - -Exercises temporal split, grouped CV, fit+calibrate, and the model-vs-baselines -comparison — the core of the MLOps story — without a tracking server. -""" - -import warnings - -import pytest - -# Needs the [ml] extra (pandas/scikit-learn/xgboost). The core lint-type-test CI -# job installs only [dev], so skip there; the dedicated `ml` job runs these fully. -pytest.importorskip("pandas") -pytest.importorskip("sklearn") -pytest.importorskip("xgboost") - -from contrail_ml.config import MLConfig -from contrail_ml.data.synthetic_fallback import make_synthetic_training_table -from contrail_ml.train import run_training - - -def test_run_training_builds_comparison_and_beats_raw(): - cfg = MLConfig(n_ensemble=3, xgb_n_estimators=60) - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - df = make_synthetic_training_table(n=5000, seed=11, allow_synthetic=True) - - res = run_training(cfg, df, log_mlflow=False, do_cv=True) - - # CV produced metrics. - assert "rhi_mae" in res.cv_metrics and res.cv_metrics["rhi_mae"] > 0 - - # The comparison table has the model row and every baseline. - table = res.comparison - for name in ("ml_corrector", "raw_era5", "x_factor", "quantile_map"): - assert name in table.index - - ml = table.loc["ml_corrector"] - raw = table.loc["raw_era5"] - # Bias removed and ISSR skill not worse than the uncorrected forecast. - assert abs(ml["rhi_bias"]) < abs(raw["rhi_bias"]) - assert ml["ets"] >= raw["ets"] - assert ml["rhi_mae"] < raw["rhi_mae"] - # Calibrated probabilities are honest. - assert ml["ece"] < 0.05 diff --git a/tests/test_ml_world_integration.py b/tests/test_ml_world_integration.py deleted file mode 100644 index 9474831..0000000 --- a/tests/test_ml_world_integration.py +++ /dev/null @@ -1,68 +0,0 @@ -"""The seam: World(issr_source="ml") -> qubo -> CP-SAT, with no qubo/World edits. - -Mirrors tests/test_solver_cpsat.py but with the ML ISSR field swapped in via the -duck-typed interface. The whole point of the project: the trained model drops -into the existing optimizer untouched. -""" - -import warnings - -import pytest - -# Needs the [ml] extra (pandas/scipy/scikit-learn/xgboost). The core lint-type-test -# CI job installs only [dev], so skip there; the dedicated `ml` job runs these fully. -pytest.importorskip("pandas") -pytest.importorskip("scipy") -pytest.importorskip("sklearn") -pytest.importorskip("xgboost") - -from contrail_env import ( - build_and_evaluate_flight, - build_capacity_buckets, - build_conflict_graph, - build_random_flights, - default_european_world, - solve_cpsat, -) -from contrail_ml.config import MLConfig -from contrail_ml.issr_field import MLIssrField - - -def _ml_world(): - cfg = MLConfig() - with warnings.catch_warnings(): - warnings.simplefilter("ignore") - return default_european_world( - seed=1, - issr_source="ml", - issr_kwargs=dict(config=cfg, met_source="synthetic", - allow_synthetic=True, grid_res_deg=3.0, seed=1), - ) - - -def test_ml_world_builds_with_ml_field(): - world = _ml_world() - assert isinstance(world.issr, MLIssrField) - - -def test_ml_world_solves_feasibly(): - world = _ml_world() - flights = build_random_flights(n_flights=4, world=world, seed=1, - corridor_frac=0.05, snapshot_window_s=(0.0, 300.0)) - evals = [] - for f in flights: - evals.extend(build_and_evaluate_flight(f, world)) - conflicts = build_conflict_graph(evals, world) - buckets = build_capacity_buckets(evals, world) - - result = solve_cpsat(evals, conflicts, buckets, time_limit_s=10.0) - selected = set(result.chosen_eval_indices) - - # One option per flight. - assert len({evals[i].flight_name for i in selected}) == 4 - # No conflict edge fully chosen. - for e in conflicts: - assert not (e.i in selected and e.j in selected) - # No capacity bucket exceeded. - for b in buckets: - assert sum(1 for m in b.members if m in selected) <= b.capacity diff --git a/tests/test_opensky_guard.py b/tests/test_opensky_guard.py deleted file mode 100644 index 32bb341..0000000 --- a/tests/test_opensky_guard.py +++ /dev/null @@ -1,39 +0,0 @@ -"""The OpenSky client fails loudly — it never returns empty/fabricated data. - -Whether or not pyopensky is installed, fetching without credentials must raise -OpenSkyUnavailableError (missing dependency OR missing credentials), so a real -pull can never silently degrade into zero or invented flights. -""" - -import pytest - -from contrail_flights.config import FlightsConfig -from contrail_flights.opensky_client import ( - OpenSkyClient, - OpenSkyUnavailableError, - tracks_from_state_dataframe, -) - - -def test_fetch_raises_without_deps_or_credentials(): - cfg = FlightsConfig(opensky_username=None, opensky_password=None) - client = OpenSkyClient(cfg) - with pytest.raises(OpenSkyUnavailableError): - client.fetch_tracks() - - -def test_dataframe_grouping_is_unit_testable(): - # The pure parsing step works on an in-memory frame (no network), so the - # grouping logic is covered even though the live query can't be in CI. - pd = pytest.importorskip("pandas") - df = pd.DataFrame({ - "icao24": ["a", "a", "b"], - "callsign": ["AAA", "AAA", "BBB"], - "time": [0.0, 60.0, 0.0], - "latitude": [46.0, 46.1, 47.0], - "longitude": [0.0, 0.5, 1.0], - "baroaltitude": [10000.0, 10500.0, 11000.0], - }) - tracks = tracks_from_state_dataframe(df) - assert {t.icao24 for t in tracks} == {"a", "b"} - assert next(t for t in tracks if t.icao24 == "a").n_points == 2 diff --git a/tests/test_real_flight_world_integration.py b/tests/test_real_flight_world_integration.py deleted file mode 100644 index 618e44e..0000000 --- a/tests/test_real_flight_world_integration.py +++ /dev/null @@ -1,67 +0,0 @@ -"""The seam: real (fixture) flights -> qubo -> CP-SAT, with no qubo/World edits. - -Mirrors test_ml_world_integration.py but swaps in real-flight objects built from -hand-made ADS-B tracks via contrail_flights. The point of the feature: real -historical traffic drops into the existing optimizer untouched. Hermetic — no -network/credentials (tracks are passed directly to flights_for). -""" - -import numpy as np - -from contrail_env import ( - build_and_evaluate_flight, - build_capacity_buckets, - build_conflict_graph, - default_european_world, - fl_to_m, - solve_cpsat, -) -from contrail_flights.build_dataset import flights_for -from contrail_flights.config import FlightsConfig -from contrail_flights.tracks import make_track - - -def _fixture_tracks(n_flights: int = 3): - """A handful of overlapping west->east crossings at adjacent latitudes.""" - tracks = [] - for i in range(n_flights): - n = 24 - t = np.arange(n) * 60.0 - lon = np.linspace(-4.5, 13.0, n) # stays inside 0..1500 km east - lat = np.full(n, 45.5 + i * 0.6) - alt = np.full(n, fl_to_m(360)) - tracks.append(make_track(f"ac{i}", f"F{i}", t, lat, lon, alt)) - return tracks - - -def test_real_flights_build(): - cfg = FlightsConfig() - flights = flights_for(cfg, tracks=_fixture_tracks(3), use_cache=False) - assert len(flights) == 3 - for f in flights: - assert 0.0 <= f.origin_km[0] <= 1500.0 - assert f.baseline.n_segments >= 1 - - -def test_real_flights_solve_feasibly(): - cfg = FlightsConfig() - flights = flights_for(cfg, tracks=_fixture_tracks(3), use_cache=False) - world = default_european_world(seed=1) - - evals = [] - for f in flights: - evals.extend(build_and_evaluate_flight(f, world)) - conflicts = build_conflict_graph(evals, world) - buckets = build_capacity_buckets(evals, world) - - result = solve_cpsat(evals, conflicts, buckets, time_limit_s=10.0) - selected = set(result.chosen_eval_indices) - - # One option per flight. - assert len({evals[i].flight_name for i in selected}) == 3 - # No conflict edge fully chosen. - for e in conflicts: - assert not (e.i in selected and e.j in selected) - # No capacity bucket exceeded. - for b in buckets: - assert sum(1 for m in b.members if m in selected) <= b.capacity