from __future__ import annotations import hashlib import json from dataclasses import asdict, dataclass from datetime import UTC, datetime from pathlib import Path from typing import Any import httpx import numpy as np import pandas as pd from .evaluation_trace import EvaluationGateResult, EvaluationMetric, EvaluationRunRecord from .hrrr_archive import MODEL_PUBLICATION_LAG, HrrrArchiveClient, HrrrSnapshotBundle from .sandbox import sha256_file from .weather import WeatherForecastSnapshot, load_weather_research_registry ANALYSIS_VERSION = "1.2.0" HORIZONS = (24, 6) BOOTSTRAP_REPETITIONS = 4_000 METHOD_AMENDMENT_LOCKED_AT = "2026-09-03T10:00:00+00:00" AVAILABILITY_AMENDMENT_LOCKED_AT = "2026-09-03T12:30:00+00:00" @dataclass(frozen=True) class HrrrBackfillRequest: request_id: str game_id: str season: int stadium_id: str latitude: float longitude: float kickoff_at: datetime cutoff_hours: int def _load_locations(path: Path) -> dict[str, dict[str, Any]]: value = json.loads(path.read_text(encoding="utf-8")) if not isinstance(value, list): raise TypeError("Stadium locations must be a JSON list") locations = {str(item["stadium_id"]).upper(): item for item in value} if len(locations) != len(value): raise ValueError("Stadium locations repeat an identifier") return locations def build_backfill_requests( *, observed_weather_path: Path, stadium_locations_path: Path, completed_request_ids: set[str] | None = None, ) -> list[HrrrBackfillRequest]: weather = pd.read_csv(observed_weather_path) weather = weather[ weather["station_rank"].eq(1) & weather["nearest_report_minutes"].notna() ].copy() if weather["game_id"].duplicated().any(): raise ValueError("Observed weather repeats a primary game record") locations = _load_locations(stadium_locations_path) missing_locations = sorted(set(weather["stadium_id"]) - set(locations)) if missing_locations: raise ValueError(f"HRRR backfill is missing stadium locations: {missing_locations}") completed = completed_request_ids or set() requests = [] for row in weather.sort_values(["season", "game_id"], ascending=[False, True]).itertuples(): location = locations[row.stadium_id] for horizon in HORIZONS: request_id = f"{row.game_id}:{horizon}h" if request_id in completed: continue kickoff = datetime.fromisoformat(str(row.kickoff_at)) requests.append( HrrrBackfillRequest( request_id=request_id, game_id=row.game_id, season=int(row.season), stadium_id=row.stadium_id, latitude=float(location["latitude"]), longitude=float(location["longitude"]), kickoff_at=kickoff, cutoff_hours=horizon, ) ) return requests def _records_hash(value: Any) -> str: def canonical(item: Any) -> Any: if isinstance(item, dict): return {key: canonical(entry) for key, entry in item.items()} if isinstance(item, list): return [canonical(entry) for entry in item] if isinstance(item, (float, np.floating)): return round(float(item), 12) if np.isfinite(item) else None if isinstance(item, (np.integer,)): return int(item) return item encoded = json.dumps(canonical(value), sort_keys=True, separators=(",", ":"), allow_nan=False) return hashlib.sha256(encoded.encode()).hexdigest() def _mean_or_none(values: pd.Series) -> float | None: mean = values.mean() return float(mean) if pd.notna(mean) and np.isfinite(mean) else None def _bootstrap_mean_interval(values: pd.Series, *, label: str) -> list[float] | None: clean = values.dropna().to_numpy(dtype=float) if not len(clean): return None seed = int.from_bytes(hashlib.sha256(label.encode()).digest()[:8], "big") rng = np.random.default_rng(seed) samples: list[np.ndarray] = [] for count in range(0, BOOTSTRAP_REPETITIONS, 200): repetitions = min(200, BOOTSTRAP_REPETITIONS - count) samples.append(rng.choice(clean, size=(repetitions, len(clean)), replace=True).mean(axis=1)) lower, upper = np.quantile(np.concatenate(samples), [0.025, 0.975]) return [float(lower), float(upper)] def score_hrrr_forecasts( *, snapshots: list[WeatherForecastSnapshot], observed_weather_path: Path, expected_requests: int, ) -> dict[str, Any]: observations = pd.read_csv(observed_weather_path) observations = observations[observations["station_rank"].eq(1)][ ["game_id", "stadium_id", "air_temp_c", "wind_speed_mps", "rain_observed"] ] forecasts = pd.DataFrame( [ { **snapshot.model_dump(mode="json"), "cutoff_hours": int(snapshot.record_id.split("-")[-2].removesuffix("h")), } for snapshot in snapshots ] ) if forecasts.empty: return { "schema_version": "1.0", "analysis_version": ANALYSIS_VERSION, "expected_requests": expected_requests, "completed_forecasts": 0, "coverage": 0, "stadium_coverage": 0, "groups": [], "result_sha256": _records_hash([]), } joined = forecasts.merge( observations, on="game_id", how="left", suffixes=("_forecast", "_observed"), validate="many_to_one", ) joined["temperature_absolute_error"] = ( joined["air_temp_c_forecast"] - joined["air_temp_c_observed"] ).abs() joined["wind_absolute_error"] = ( joined["wind_speed_mps_forecast"] - joined["wind_speed_mps_observed"] ).abs() observed_rain = joined["rain_observed"].map({True: 1.0, False: 0.0}) joined["precipitation_brier"] = (joined["precip_probability"] - observed_rain) ** 2 groups = [] for (provider_model, horizon), frame in joined.groupby( ["provider_model", "cutoff_hours"], sort=True ): calibration = [] for predicted, bucket in frame.assign(observed_rain=observed_rain.loc[frame.index]).groupby( "precip_probability", dropna=True, sort=True ): calibration.append( { "predicted_probability": float(predicted), "forecasts": len(bucket), "observed_rain_rate": _mean_or_none(bucket["observed_rain"]), } ) groups.append( { "provider_model": provider_model, "cutoff_hours": int(horizon), "forecasts": len(frame), "temperature_pairs": int(frame["temperature_absolute_error"].notna().sum()), "temperature_mae_c": _mean_or_none(frame["temperature_absolute_error"]), "temperature_mae_ci_95_c": _bootstrap_mean_interval( frame["temperature_absolute_error"], label=f"{provider_model}:{horizon}:temperature-mae", ), "wind_pairs": int(frame["wind_absolute_error"].notna().sum()), "wind_mae_mps": _mean_or_none(frame["wind_absolute_error"]), "wind_mae_ci_95_mps": _bootstrap_mean_interval( frame["wind_absolute_error"], label=f"{provider_model}:{horizon}:wind-mae", ), "precipitation_pairs": int(frame["precipitation_brier"].notna().sum()), "precipitation_brier_score": _mean_or_none(frame["precipitation_brier"]), "precipitation_brier_ci_95": _bootstrap_mean_interval( frame["precipitation_brier"], label=f"{provider_model}:{horizon}:precipitation-brier", ), "precipitation_calibration": calibration, "forecast_missingness": float( frame[["air_temp_c_forecast", "wind_speed_mps_forecast", "precip_mm"]] .isna() .mean() .mean() ), } ) horizon_comparisons = [] for provider_model, frame in joined.groupby("provider_model", sort=True): comparison: dict[str, Any] = { "provider_model": provider_model, "comparison": "6-hour minus 24-hour forecast error", } pair_counts: list[int] = [] for metric, output_name in ( ("temperature_absolute_error", "temperature_mae_difference_c"), ("wind_absolute_error", "wind_mae_difference_mps"), ("precipitation_brier", "precipitation_brier_difference"), ): paired = frame.pivot(index="game_id", columns="cutoff_hours", values=metric) if not {6, 24}.issubset(paired.columns): difference = pd.Series(dtype=float) else: difference = (paired[6] - paired[24]).dropna() comparison[output_name] = _mean_or_none(difference) comparison[f"{output_name}_ci_95"] = _bootstrap_mean_interval( difference, label=f"{provider_model}:6h-minus-24h:{metric}", ) pair_counts.append(len(difference)) comparison["complete_pairs"] = min(pair_counts) if pair_counts else 0 horizon_comparisons.append(comparison) result = { "schema_version": "1.0", "analysis_version": ANALYSIS_VERSION, "expected_requests": expected_requests, "completed_forecasts": len(joined), "coverage": len(joined) / expected_requests if expected_requests else 0, "stadium_coverage": ( joined["stadium_id"].nunique() / observations["stadium_id"].nunique() if observations["stadium_id"].nunique() else 0 ), "groups": groups, "horizon_comparisons": horizon_comparisons, "result_sha256": "", } result["result_sha256"] = _records_hash(result | {"result_sha256": None}) return result def _write_json(path: Path, value: Any) -> None: path.parent.mkdir(parents=True, exist_ok=True) path.write_text(json.dumps(value, indent=2, sort_keys=True, allow_nan=False), encoding="utf-8") def _migrate_operational_availability( *, snapshots: list[WeatherForecastSnapshot], failures: list[dict[str, Any]], receipt_dir: Path, ) -> tuple[list[WeatherForecastSnapshot], dict[str, Any]]: """Move v1.1 evidence from archive-upload time to conservative operational time.""" candidates: dict[str, tuple[WeatherForecastSnapshot, str]] = {} dropped_record_ids: set[str] = set() missing_receipts = 0 for snapshot in snapshots: model_run = datetime.fromisoformat(snapshot.model_run_at) cutoff = datetime.fromisoformat(snapshot.cutoff_at) source_available_at = model_run + MODEL_PUBLICATION_LAG if ( source_available_at > cutoff or not (receipt_dir / f"{snapshot.record_id}.json").is_file() ): dropped_record_ids.add(snapshot.record_id) if not (receipt_dir / f"{snapshot.record_id}.json").is_file(): missing_receipts += 1 continue available = source_available_at.isoformat() candidates[snapshot.record_id] = snapshot, available migrated_receipts = 0 removed_receipts = 0 bundle_hash_by_record: dict[str, str] = {} if receipt_dir.exists(): for path in sorted(receipt_dir.glob("*.json")): if path.stem in dropped_record_ids: path.unlink() removed_receipts += 1 continue candidate = candidates.get(path.stem) if not candidate: continue _, available = candidate value = json.loads(path.read_text(encoding="utf-8")) for source in value.get("sources", []): archive_modified = source.get("archive_object_modified_at") or source.get( "source_available_at" ) source["source_available_at"] = available source["archive_object_modified_at"] = archive_modified bundle_hash = hashlib.sha256( json.dumps( value["sources"], sort_keys=True, separators=(",", ":"), ).encode() ).hexdigest() value["snapshot_sha256"] = bundle_hash bundle_hash_by_record[path.stem] = bundle_hash _write_json(path, value) migrated_receipts += 1 migrated = sorted( ( snapshot.model_copy( update={ "source_sha256": bundle_hash_by_record[record_id], "issued_at": available, "source_available_at": available, } ) for record_id, (snapshot, available) in candidates.items() ), key=lambda item: item.record_id, ) summary = { "schema_version": "1.0", "analysis_version": ANALYSIS_VERSION, "locked_at": AVAILABILITY_AMENDMENT_LOCKED_AT, "prior_snapshots": len(snapshots), "retained_snapshots": len(migrated), "snapshots_requeued": len(dropped_record_ids), "prior_failures_requeued": len(failures), "receipts_migrated": migrated_receipts, "receipts_removed_for_requeue": removed_receipts, "snapshots_requeued_for_missing_receipts": missing_receipts, "source_availability_rule": ( "HRRR model initialization plus a conservative four-hour publication buffer" ), "archive_object_last_modified_role": "retrieval provenance only", } return migrated, summary def run_hrrr_backfill_batch( *, observed_weather_path: Path, stadium_locations_path: Path, output_dir: Path, artifact_root: str, code_sha: str, run_id: str, batch_size: int = 4, client: HrrrArchiveClient | None = None, grid_cache_path: Path | None = None, retrieved_at: datetime | None = None, ) -> EvaluationRunRecord: if batch_size < 1 or batch_size > 100: raise ValueError("HRRR batch size must be between 1 and 100") output_dir.mkdir(parents=True, exist_ok=True) registration_path = output_dir / "registration.json" registration = { "schema_version": "1.0", "study_id": "weather-forecast-skill", "analysis_version": ANALYSIS_VERSION, "forecast_horizons_hours": list(HORIZONS), "game_window_hours_relative_to_kickoff": [-1, 0, 1, 2, 3, 4], "model_cycles": "00/06/12/18 UTC extended HRRR cycles", "cycle_selection_buffer_minutes": 240, "source_availability_proof": ( "The registered HRRR model initialization precedes the cutoff by at least four " "hours. NOAA documents operational extended cycles at 00/06/12/18 UTC; immutable " "archive hashes identify the exact forecast values." ), "archive_object_last_modified_role": ( "Retrieval provenance only; historical objects were sometimes re-archived after " "their operational model cycles and that upload time is not forecast issuance." ), "outcomes": [ "game-window-temperature-mae", "game-window-wind-mae", "deterministic-game-window-precipitation-brier", ], "minimum_complete_forecasts": 1000, "minimum_coverage": 0.85, "version_strata": ["HRRRv3", "HRRRv4"], "uncertainty": ( "4,000 deterministic nonparametric bootstrap draws for each error mean and each " "paired 6-hour-minus-24-hour difference" ), "method_amendment": { "from_analysis_version": "1.0.0", "locked_at": METHOD_AMENDMENT_LOCKED_AT, "locked_before_release_sample": True, "reason": ( "The pilot scorecard reported point errors only. Version 1.1.0 preregistered " "uncertainty before the 85% release sample and before public forecast-skill claims." ), }, "availability_method_amendment": { "from_analysis_version": "1.1.0", "locked_at": AVAILABILITY_AMENDMENT_LOCKED_AT, "locked_before_recovery_sample": True, "reason": ( "Version 1.1 incorrectly treated archive-object Last-Modified as historical " "forecast issuance. Version 1.2 uses a conservative operational-cycle buffer " "for the prediction boundary and retains Last-Modified only as provenance." ), }, "observations_are_research_only": True, } if registration_path.exists(): prior_registration = json.loads(registration_path.read_text(encoding="utf-8")) if prior_registration != registration: if ( prior_registration.get("analysis_version") == "1.0.0" and registration["analysis_version"] == "1.1.0" and "uncertainty" not in prior_registration ): _write_json( output_dir / "registration-history" / "1.0.0.json", prior_registration, ) _write_json(registration_path, registration) elif ( prior_registration.get("analysis_version") == "1.1.0" and registration["analysis_version"] == "1.2.0" and "availability_method_amendment" not in prior_registration ): _write_json( output_dir / "registration-history" / "1.1.0.json", prior_registration, ) _write_json(registration_path, registration) else: raise ValueError("The HRRR backfill registration changed without a version bump") else: _write_json(registration_path, registration) snapshots_path = output_dir / "snapshots.json" snapshots = [ WeatherForecastSnapshot.model_validate(item) for item in ( json.loads(snapshots_path.read_text(encoding="utf-8")) if snapshots_path.exists() else [] ) ] failures_path = output_dir / "unavailable.json" failures = ( json.loads(failures_path.read_text(encoding="utf-8")) if failures_path.exists() else [] ) migration_path = output_dir / "availability-method-migration-1.2.0.json" receipt_dir = output_dir / "source-receipts" if not migration_path.exists() and (snapshots or failures): snapshots, migration_record = _migrate_operational_availability( snapshots=snapshots, failures=failures, receipt_dir=receipt_dir, ) failures = [] _write_json( snapshots_path, [snapshot.model_dump(mode="json") for snapshot in snapshots], ) _write_json(failures_path, failures) _write_json(migration_path, migration_record) snapshots_before = len(snapshots) failures_before = len(failures) completed = { f"{snapshot.game_id}:{int(snapshot.record_id.split('-')[-2].removesuffix('h'))}h" for snapshot in snapshots } | {str(item["request_id"]) for item in failures} all_requests = build_backfill_requests( observed_weather_path=observed_weather_path, stadium_locations_path=stadium_locations_path, ) requests = [request for request in all_requests if request.request_id not in completed][ :batch_size ] current_evaluation_path = output_dir / "evaluation-run.json" if not requests and current_evaluation_path.exists(): return EvaluationRunRecord.model_validate_json( current_evaluation_path.read_text(encoding="utf-8") ) archive = client or HrrrArchiveClient() retrieved = (retrieved_at or datetime.now(UTC)).astimezone(UTC) for request in requests: try: bundle: HrrrSnapshotBundle = archive.fetch_snapshot( game_id=request.game_id, latitude=request.latitude, longitude=request.longitude, kickoff_at=request.kickoff_at, cutoff_hours=request.cutoff_hours, retrieved_at=retrieved, cache_path=grid_cache_path, ) except ValueError as exc: failures.append( { "request_id": request.request_id, "reason": f"{type(exc).__name__}: {exc}", "failed_closed": True, "recorded_at": retrieved.isoformat(), } ) continue except httpx.HTTPStatusError as exc: if exc.response.status_code != 404: raise failures.append( { "request_id": request.request_id, "reason": f"HTTPStatusError: {exc}", "failed_closed": True, "recorded_at": retrieved.isoformat(), } ) continue snapshots.append(bundle.snapshot) _write_json( receipt_dir / f"{bundle.snapshot.record_id}.json", { "snapshot_sha256": bundle.snapshot.source_sha256, "grid_point": asdict(bundle.grid_point), "sources": [asdict(receipt) for receipt in bundle.receipts], }, ) snapshots.sort(key=lambda item: item.record_id) failures.sort(key=lambda item: item["request_id"]) _write_json(snapshots_path, [snapshot.model_dump(mode="json") for snapshot in snapshots]) _write_json(failures_path, failures) first = score_hrrr_forecasts( snapshots=snapshots, observed_weather_path=observed_weather_path, expected_requests=len(all_requests), ) second = score_hrrr_forecasts( snapshots=snapshots, observed_weather_path=observed_weather_path, expected_requests=len(all_requests), ) reproducible = first["result_sha256"] == second["result_sha256"] _write_json(output_dir / "scorecard.json", first) _write_json( output_dir / "reproducibility.json", { "first_sha256": first["result_sha256"], "second_sha256": second["result_sha256"], "matched": reproducible, }, ) minimum = next( study.minimum_complete_games for study in load_weather_research_registry().studies if study.study_id == "weather-forecast-skill" ) gates = [ EvaluationGateResult( gate_id="cutoff-integrity", status="passed", details="Every archived forecast proves source availability before its cutoff.", evidence_hrefs=[f"{artifact_root}/snapshots.json"], ), EvaluationGateResult( gate_id="forecast-sample", status="passed" if len(snapshots) >= minimum else "failed", details=f"{len(snapshots)} cutoff-safe forecasts are complete; {minimum} are required.", evidence_hrefs=[f"{artifact_root}/scorecard.json"], ), EvaluationGateResult( gate_id="forecast-coverage", status="passed" if first["coverage"] >= 0.85 else "failed", details=f"Backfill coverage is {first['coverage']:.1%}; 85% is required.", evidence_hrefs=[f"{artifact_root}/scorecard.json"], ), EvaluationGateResult( gate_id="version-stratification", status=( "passed" if {group["provider_model"] for group in first["groups"]} >= {"HRRRv3", "HRRRv4"} else "failed" ), details="Forecast skill is reported separately for HRRRv3 and HRRRv4 eras.", evidence_hrefs=[f"{artifact_root}/scorecard.json"], ), EvaluationGateResult( gate_id="reproducibility", status="passed" if reproducible else "failed", details="Two scorecard executions produced the same canonical result hash.", evidence_hrefs=[f"{artifact_root}/reproducibility.json"], ), ] status = "passed" if all(gate.status == "passed" for gate in gates) else "withheld" evaluation = EvaluationRunRecord( run_id=run_id, workload_id="weather-forecast-skill", run_kind="forecast-backtest", status=status, started_at=retrieved.isoformat(), completed_at=max(datetime.now(UTC), retrieved).isoformat(), code_sha=code_sha, source_sha256={ "observed-weather": sha256_file(observed_weather_path), "stadium-locations": sha256_file(stadium_locations_path), }, registration_href=f"{artifact_root}/registration.json", artifact_hrefs=[ f"{artifact_root}/{name}" for name in ( "snapshots.json", "unavailable.json", "scorecard.json", "reproducibility.json", "availability-method-migration-1.2.0.json", ) if (output_dir / name).exists() ], metrics=[ EvaluationMetric( metric_id="cutoff-safe-forecasts", label="Cutoff-safe archived forecasts", value=len(snapshots), units="forecasts", population="Registered HRRR stadium-game horizons", ), EvaluationMetric( metric_id="backfill-coverage", label="Archived forecast backfill coverage", value=first["coverage"], units="share", population="Registered HRRR stadium-game horizons", ), EvaluationMetric( metric_id="batch-requests-attempted", label="Requests attempted in this batch", value=len(requests), units="requests", population="This bounded backfill run", ), EvaluationMetric( metric_id="batch-forecasts-added", label="Forecast snapshots added in this batch", value=len(snapshots) - snapshots_before, units="forecasts", population="This bounded backfill run", ), EvaluationMetric( metric_id="batch-unavailable-recorded", label="Permanent source gaps recorded in this batch", value=len(failures) - failures_before, units="requests", population="This bounded backfill run", ), ], gates=gates, total_cost_usd=0, ) evaluation_value = evaluation.model_dump(mode="json") _write_json(current_evaluation_path, evaluation_value) _write_json(output_dir / "evaluation-runs" / f"{run_id}.json", evaluation_value) return evaluation