"""Fetch one country's inputs and print suiren 1's forecast through the end of tomorrow, UTC. Demand from Energy-Charts, weather from Open-Meteo, holidays from Nager.Date: the sources the shipped weights trained on. Open-Meteo's free API is for non-commercial use. See GUIDE.md. python -m eval.forecast --zone de --seeds 0,1,2 > de.csv python -m eval.forecast --zone de --json > de.json """ from __future__ import annotations import argparse import datetime as dt import http.client import json import logging import sys import time import urllib.error import urllib.parse import urllib.request from collections.abc import Callable, Sequence from pathlib import Path from typing import NamedTuple, TypedDict import numpy as np from eval.suiren import Suiren log = logging.getLogger("forecast") ROOT = Path(__file__).resolve().parents[1] CONTEXT = 2048 LEVELS = (0.05, 0.5, 0.95) COLUMNS = ("hour_utc", "p05_mw", "p50_mw", "p95_mw") MODEL = "suiren-1" STAMP = "%Y-%m-%dT%H:%M:%SZ" HOUR = np.timedelta64(1, "h") LOAD_URL = "https://api.energy-charts.info/public_power" WEATHER_URL = "https://historical-forecast-api.open-meteo.com/v1/forecast" HOLIDAY_URL = "https://date.nager.at/api/v3/PublicHolidays/{year}/{country}" WEATHER = ("temperature_2m", "wind_speed_10m", "shortwave_radiation", "cloud_cover") # gb is left out: its target is Elexon's demand outturn, which Energy-Charts does not serve. ZONES = ("at ba be bg ch cz de dk ee es fi fr ge gr hr hu ie it lt lu lv md nl no pl pt ro rs se " "si sk xk").split() # Nager.Date does not serve Kosovo, so training read every Kosovo day as a working day. NO_HOLIDAY_TABLE = ("xk",) # Energy-Charts answers 429 with Retry-After after a burst of about a dozen requests, and a dropped # connection is transport noise; ATTEMPTS is shared by both. ATTEMPTS = 4 RETRY_AFTER_DEFAULT = 15 RETRY_AFTER_CAP = 60 SERVER_ERRORS = (500, 502, 503, 504) # The model forecasts a day ahead, so a feed this late is a dead feed. MAX_LAG_HOURS = 48 sleep = time.sleep # a name of its own so a test can swap it class Row(TypedDict): hour_utc: str # "2026-10-26T08:00Z" p05_mw: int | None # None when the band is NaN p50_mw: int | None p95_mw: int | None class Forecast(TypedDict): model: str zone: str issued_at: str # the time this zone's run started first_hour_utc: str unit: str seeds: list[int] missing_history_hours: int rows: list[Row] class Latest(TypedDict): model: str run_issued_at: str # the last run, success or not zones: dict[str, Forecast] # last good record per zone failed: dict[str, dict] # this run only: {zone: {"error": str, "at": str}} zone_info: dict[str, dict] # {zone: {"country": str, "note": str}} def retry_wait(failure: urllib.error.HTTPError, attempt: int) -> int | None: """Seconds to wait before the next try, or None when the status is not worth a retry.""" if failure.code == 429: try: seconds = int(failure.headers.get("Retry-After", RETRY_AFTER_DEFAULT)) except (TypeError, ValueError): # Retry-After may also be an HTTP date seconds = RETRY_AFTER_DEFAULT return min(max(seconds, 0), RETRY_AFTER_CAP) return 2 ** attempt if failure.code in SERVER_ERRORS else None def get_json(url: str, params: dict | None = None) -> dict | list: query = f"?{urllib.parse.urlencode(params)}" if params else "" request = urllib.request.Request(url + query, headers={"User-Agent": "suiren-1-forecast"}) for attempt in range(ATTEMPTS): try: with urllib.request.urlopen(request, timeout=120) as answer: return json.load(answer) except urllib.error.HTTPError as failure: wait, reason = retry_wait(failure, attempt), f"HTTP {failure.code}" if wait is None or attempt == ATTEMPTS - 1: raise except (urllib.error.URLError, TimeoutError, ConnectionError, http.client.IncompleteRead) as failure: wait, reason = 2 ** attempt, type(failure).__name__ if attempt == ATTEMPTS - 1: raise log.info("%s from %s; waiting %d s (try %d of %d)", reason, urllib.parse.urlsplit(url).netloc, wait, attempt + 1, ATTEMPTS) sleep(wait) def as_floats(values: list) -> np.ndarray: return np.array([np.nan if value is None else value for value in values], dtype=float) def hourly_load(payload: dict) -> tuple[np.ndarray, np.ndarray]: """Sub-hourly points average to the hour; an hour missing any point is a gap, as in training.""" seconds = np.asarray(payload["unix_seconds"], dtype=np.int64) load = next(item["data"] for item in payload["production_types"] if item["name"] == "Load") values = as_floats(load) per_hour = 3600 // int(np.median(np.diff(seconds))) hours, index = np.unique((seconds // 3600).astype("datetime64[h]"), return_inverse=True) present = ~np.isnan(values) counts = np.bincount(index, weights=present, minlength=hours.size) sums = np.bincount(index, weights=np.where(present, values, 0.0), minlength=hours.size) return hours, np.where(counts == per_hour, sums / np.maximum(counts, 1), np.nan) def on_grid(hours: np.ndarray, values: np.ndarray, grid: np.ndarray) -> np.ndarray: out = np.full(grid.size, np.nan) position = ((hours - grid[0]) // HOUR).astype(int) keep = (position >= 0) & (position < grid.size) out[position[keep]] = values[keep] return out def calendar(grid: np.ndarray, holidays: np.ndarray) -> np.ndarray: """Columns 0 to 3 of the covariates. Hours are UTC, so the holiday flag follows the UTC date.""" days = grid.astype("datetime64[D]") return np.column_stack([ grid.astype(int) % 24, (days.astype(int) + 3) % 7, grid.astype("datetime64[M]").astype(int) % 12 + 1, np.isin(days, holidays).astype(float), ]).astype(float) def fetch_load(zone: str, first: np.datetime64, last: np.datetime64) -> tuple[np.ndarray, np.ndarray]: payload = get_json(LOAD_URL, {"country": zone, "start": str(first.astype("datetime64[D]")), "end": str(last.astype("datetime64[D]"))}) return hourly_load(payload) def fetch_weather(point: dict, grid: np.ndarray) -> np.ndarray: """The forecast issued the day before each hour, as in training; for tomorrow, today's.""" columns = [f"{name}_previous_day1" for name in WEATHER] payload = get_json(WEATHER_URL, { "latitude": point["latitude"], "longitude": point["longitude"], "start_date": str(grid[0].astype("datetime64[D]")), "end_date": str(grid[-1].astype("datetime64[D]")), "hourly": ",".join(columns), "timezone": "UTC"}) hourly = payload["hourly"] hours = np.array(hourly["time"], dtype="datetime64[h]") return np.column_stack([on_grid(hours, as_floats(hourly[name]), grid) for name in columns]) def fetch_holidays(zone: str, grid: np.ndarray) -> np.ndarray: """Only days typed Public and flagged global count, as in training.""" if zone in NO_HOLIDAY_TABLE: return np.array([], dtype="datetime64[D]") years = range(int(str(grid[0])[:4]), int(str(grid[-1])[:4]) + 1) days = {item["date"] for year in years for item in get_json(HOLIDAY_URL.format(year=year, country=zone.upper())) if item.get("global") and "Public" in item.get("types", ())} return np.array(sorted(days), dtype="datetime64[D]") def weather_point(zone: str) -> dict: shipped = ROOT / "weather_points.json" path = shipped if shipped.exists() else ROOT / "etl" / "weather_points.json" return next(point for point in json.loads(path.read_text()) if point["zone"] == zone) class Fetchers(NamedTuple): load: Callable[[str, np.datetime64, np.datetime64], tuple[np.ndarray, np.ndarray]] weather: Callable[[dict, np.ndarray], np.ndarray] holidays: Callable[[str, np.ndarray], np.ndarray] FETCH = Fetchers(fetch_load, fetch_weather, fetch_holidays) def inputs(zone: str, now: np.datetime64, hours: int | None, fetch: Fetchers = FETCH) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """The grid, the 2,048-hour history and the covariates, starting at the first unpublished hour.""" seen, load = fetch.load(zone, now - (CONTEXT + 72) * HOUR, now) origin = seen[~np.isnan(load)][-1] + HOUR lag = int((now - origin) // HOUR) if lag > MAX_LAG_HOURS: raise ValueError(f"{zone}: the last published load hour is {origin - HOUR} UTC, and the first missing " f"hour is {lag} hours old (limit {MAX_LAG_HOURS}). The model forecasts a day ahead, " "so it refuses older data.") end = (now.astype("datetime64[D]") + 2).astype("datetime64[h]") horizon = hours or int((end - origin) // HOUR) grid = np.arange(origin - CONTEXT * HOUR, origin + horizon * HOUR, HOUR) history = on_grid(seen, load, grid)[:CONTEXT] log.info("%s: forecast from %s UTC for %d hours; %d of %d history hours missing", zone, origin, horizon, int(np.isnan(history).sum()), CONTEXT) covariates = np.column_stack([calendar(grid, fetch.holidays(zone, grid)), fetch.weather(weather_point(zone), grid)]) return grid[CONTEXT:], history, covariates def load_models(seeds: Sequence[int], device: str | None, weights: Path) -> dict[int, Suiren]: return {seed: Suiren.from_checkpoint(weights / f"v31_final_stepC_seed{seed}.pt", device=device) for seed in seeds} def megawatts(value: float) -> int | None: return None if np.isnan(value) else int(round(float(value))) def forecast(zone: str, models: dict[int, Suiren], now: dt.datetime | None = None, hours: int | None = None, fetch: Fetchers = FETCH) -> Forecast: """`now` must be timezone-aware. The bands are the mean over the models.""" if now is not None and now.tzinfo is None: raise ValueError("now must be timezone-aware: a naive time would be read as local") moment = (now or dt.datetime.now(dt.timezone.utc)).astimezone(dt.timezone.utc) start = np.datetime64(moment.replace(tzinfo=None), "h") target, history, covariates = inputs(zone, start, hours, fetch) runs = [model.predict(history, horizon=target.size, levels=LEVELS, covariates=covariates, zone=zone) for model in models.values()] bands = {level: np.mean([run[level] for run in runs], axis=0) for level in LEVELS} rows = [Row(hour_utc=np.datetime_as_string(hour, unit="m") + "Z", **{name: megawatts(bands[level][index]) for level, name in zip(LEVELS, COLUMNS[1:])}) for index, hour in enumerate(target)] return Forecast(model=MODEL, zone=zone, issued_at=moment.strftime(STAMP), first_hour_utc=rows[0]["hour_utc"], unit="MW", seeds=list(models), missing_history_hours=int(np.isnan(history).sum()), rows=rows) def main(argv: list[str] | None = None) -> int: parser = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) parser.add_argument("--zone", required=True, choices=ZONES) parser.add_argument("--seeds", default="0", help="weight files to average, e.g. 0,1,2") parser.add_argument("--hours", type=int, default=None, help="hours to forecast; default runs to the end of tomorrow, UTC") parser.add_argument("--device", default=None, help="cpu or cuda; default uses a GPU when torch finds one") parser.add_argument("--weights", type=Path, default=ROOT / "weights", help="folder holding the .pt files") parser.add_argument("--json", action="store_true", help="print one JSON record instead of CSV") args = parser.parse_args(argv) logging.basicConfig(level=logging.INFO, format="%(message)s", stream=sys.stderr) models = load_models([int(seed) for seed in args.seeds.split(",")], args.device, args.weights) # Passed by name so a test can swap FETCH; a default argument would bind the original. record = forecast(args.zone, models, hours=args.hours, fetch=FETCH) if args.json: sys.stdout.write(json.dumps(record, allow_nan=False) + "\n") return 0 sys.stdout.write(",".join(COLUMNS) + "\n") for row in record["rows"]: bands = ("nan" if row[name] is None else f"{row[name]:.0f}" for name in COLUMNS[1:]) sys.stdout.write(f"{row['hour_utc'][:13]},{','.join(bands)}\n") return 0 if __name__ == "__main__": sys.exit(main())