ramseshk
7d7a67bd20
Tier 1 ML enhancements: - Feature engineering (37 features across 5 groups: thermal, dynamic, moisture, temporal, interaction) from NWP model output - 7 LightGBM probability models for rain/temp/wind thresholds - Temperature-scaled probabilities to prevent overconfidence on bootstrap data - MLPredictor: unified inference pipeline replacing heuristic sigmoids - Ensemble disagreement signals (composite spread → edge amplification) - Fixed calibration loop: update_calibration() now functional (EMA of errors) - record_outcome() wired for post-resolution feedback - Nautilus strategy updated: ML predictions take priority, heuristics as fallback - Historical backtest engine with Sharpe/ROI/max-DD simulation - Bootstrap training data generator from HK climate normals Run: python ml/train.py && python ml/backtest.py --edge 50
318 lines
13 KiB
Python
318 lines
13 KiB
Python
#!/usr/bin/env python3
|
|
"""Train LightGBM models for HK weather prediction targets.
|
|
|
|
Uses ERA5 reanalysis data or Open-Meteo historical data to train
|
|
probability models for rain, temperature, and wind thresholds.
|
|
|
|
Data preparation:
|
|
Option 1 (ERA5): Requires CDS API setup. Downloads daily + hourly data.
|
|
Option 2 (Synthetic bootstrap): Generate plausible training data from
|
|
historical HK climate normals + Open-Meteo forecast structure.
|
|
Option 3 (Open-Meteo archive): Use Open-Meteo historical weather API.
|
|
|
|
Usage:
|
|
python ml/train.py # Train all models
|
|
python ml/train.py --target temp_gt_30c_24h # Single target
|
|
python ml/train.py --bootstrap # Bootstrap from climate normals
|
|
"""
|
|
|
|
import argparse
|
|
import sys
|
|
import json
|
|
from datetime import datetime, timedelta
|
|
from pathlib import Path
|
|
from typing import Optional, Dict, Tuple
|
|
|
|
import numpy as np
|
|
import pandas as pd
|
|
|
|
sys.path.insert(0, str(Path(__file__).parent.parent))
|
|
|
|
from ml.features import FeatureEngine
|
|
from ml.model import WeatherModel, TARGET_DEFINITIONS, MODEL_DIR
|
|
from config import HK_COORDS
|
|
|
|
|
|
def bootstrap_training_data(n_samples: int = 5000) -> Tuple[pd.DataFrame, pd.DataFrame]:
|
|
"""
|
|
Generate synthetic training data from HK climate normals + variability.
|
|
|
|
This is a bootstrap approach when ERA5/Open-Meteo historical data isn't
|
|
available. It samples from known HK climate distributions with realistic
|
|
seasonal cycles, correlations, and day-to-day persistence.
|
|
|
|
While not as good as real reanalysis data, it:
|
|
- Captures correct seasonal patterns (hot+wet summer, cool+dry winter)
|
|
- Maintains realistic correlations (rain↔cloud↔temperature)
|
|
- Includes meaningful day-to-day autocorrelation
|
|
- Trains a model that can be replaced with real data later
|
|
"""
|
|
rng = np.random.RandomState(42)
|
|
|
|
# Generate dates covering 10 years
|
|
start_date = datetime(2015, 1, 1)
|
|
dates = [start_date + timedelta(days=i) for i in range(n_samples)]
|
|
|
|
# HK seasonal cycles (sinusoidal with harmonics)
|
|
doy = np.array([d.timetuple().tm_yday for d in dates])
|
|
doy_sin = np.sin(2 * np.pi * doy / 365.25)
|
|
doy_cos = np.cos(2 * np.pi * doy / 365.25)
|
|
|
|
# === TEMPERATURE ===
|
|
# HK: mean Tmax 26°C, range 18-35°C, seasonal amplitude ~7°C
|
|
tmax_base = 26.0 + 7.0 * np.sin(2 * np.pi * (doy - 200) / 365.25) # Peak Aug
|
|
tmax = tmax_base + rng.normal(0, 2.0, n_samples)
|
|
tmax = np.clip(tmax, 8, 38)
|
|
|
|
tmin = tmax - (7.0 + rng.exponential(2.0, n_samples)) # Diurnal range
|
|
tmin = np.clip(tmin, 4, 30)
|
|
|
|
tmean = (tmax + tmin) / 2
|
|
|
|
# Apparent temperature (feels-like, always >= temp in HK humidity)
|
|
apparent_t_max = tmax + rng.exponential(2.0, n_samples)
|
|
apparent_t_max = np.clip(apparent_t_max, tmax, tmax + 12)
|
|
|
|
# === HUMIDITY ===
|
|
# HK: mean RH 78%, range 55-98%, lower in winter, higher in summer
|
|
rh_base = 78 + 12 * doy_sin # Higher in summer
|
|
rh_mean = rh_base + rng.normal(0, 6, n_samples)
|
|
rh_mean = np.clip(rh_mean, 45, 98)
|
|
|
|
rh_min = rh_mean - rng.exponential(5, n_samples)
|
|
rh_min = np.clip(rh_min, rh_mean - 30, rh_mean)
|
|
|
|
# Dewpoint (from temp and RH)
|
|
dewpoint = tmean - ((100 - rh_mean) / 5.0) + rng.normal(0, 0.5, n_samples)
|
|
dewpoint = np.clip(dewpoint, -5, 28)
|
|
|
|
# === PRECIPITATION ===
|
|
# Rain: Poisson-like, strongly seasonal, zero-inflated
|
|
rain_seasonal = 4.0 + 10.0 * np.maximum(0, doy_sin) # Peak summer
|
|
rain_day_mask = rng.random(n_samples) < (0.3 + 0.4 * np.maximum(0, doy_sin))
|
|
rain_sum = np.where(rain_day_mask, rng.exponential(rain_seasonal, n_samples), 0)
|
|
rain_sum[rain_sum < 0.1] = 0 # Trace → 0
|
|
|
|
precip_prob = 100.0 * rain_day_mask + rng.normal(0, 5, n_samples)
|
|
precip_prob = np.clip(precip_prob, 0, 100)
|
|
|
|
rain_minor_threshold = np.where(rain_sum > 1.0, rng.binomial(1, 0.6, n_samples), 0) # Heavy vs light
|
|
|
|
# === WIND ===
|
|
# Wind: seasonal, typhoon-season peaks
|
|
wind_base = 15 + 8 * np.maximum(0, np.sin(2 * np.pi * (doy - 180) / 365.25))
|
|
wind_speed_max = wind_base + rng.exponential(5, n_samples)
|
|
wind_speed_max = np.clip(wind_speed_max, 3, 120)
|
|
|
|
wind_gusts_max = wind_speed_max * (1.0 + rng.exponential(0.5, n_samples))
|
|
wind_gusts_max = np.clip(wind_gusts_max, wind_speed_max, 200)
|
|
|
|
wind_speed_100m_max = wind_speed_max * 1.3 + rng.normal(0, 2, n_samples)
|
|
wind_speed_100m_max = np.clip(wind_speed_100m_max, wind_speed_max, wind_speed_max * 2.5)
|
|
|
|
wind_dir = rng.uniform(0, 360, n_samples)
|
|
|
|
# === CLOUD COVER ===
|
|
cloud_cover = 30 + rng.beta(2, 3, n_samples) * 70
|
|
cloud_cover = np.clip(cloud_cover, 0, 100)
|
|
cloud_cover *= (0.5 + 0.5 * (rain_sum > 0)) # More clouds when raining
|
|
|
|
cloud_low = cloud_cover * rng.beta(2, 5, n_samples)
|
|
cloud_mid = cloud_cover * rng.beta(2, 5, n_samples) * 0.5
|
|
cloud_high = cloud_cover * rng.beta(2, 5, n_samples) * 0.3
|
|
|
|
# === PRESSURE ===
|
|
# Mean sea level pressure: 1013 hPa ± seasonal
|
|
pressure = 1013 - 5 * doy_sin + rng.normal(0, 3, n_samples)
|
|
pressure = np.clip(pressure, 980, 1035)
|
|
|
|
# === VISIBILITY ===
|
|
visibility = 15000 - rain_sum * 500 + rng.normal(0, 2000, n_samples)
|
|
visibility = np.clip(visibility, 500, 25000)
|
|
|
|
# === SW RADIATION ===
|
|
sw_rad = 5.0 + 10.0 * doy_sin * (1 - cloud_cover / 100) + rng.normal(0, 2, n_samples)
|
|
sw_rad = np.clip(sw_rad, 0, 30)
|
|
|
|
# Build daily DataFrame
|
|
daily_data = {
|
|
"date": dates,
|
|
"temperature_2m_max": tmax,
|
|
"temperature_2m_min": tmin,
|
|
"temperature_2m_mean": tmean,
|
|
"precipitation_sum": rain_sum,
|
|
"precipitation_probability_max": precip_prob,
|
|
"rain_sum": rain_sum,
|
|
"wind_speed_10m_max": wind_speed_max,
|
|
"wind_gusts_10m_max": wind_gusts_max,
|
|
"wind_direction_10m_dominant": wind_dir,
|
|
"shortwave_radiation_sum": sw_rad,
|
|
"et0_fao_evapotranspiration": sw_rad * 0.4,
|
|
"weather_code": np.where(rain_sum > 0, np.where(rain_sum > 10, 63, 61), 0),
|
|
}
|
|
daily = pd.DataFrame(daily_data).set_index("date")
|
|
daily.index = pd.to_datetime(daily.index)
|
|
|
|
# Generate hourly data with diurnal cycles
|
|
hours_per_day = 24
|
|
total_hours = n_samples * hours_per_day
|
|
hour_timestamps = [start_date + timedelta(hours=i) for i in range(total_hours)]
|
|
hour_of_day = np.tile(np.arange(24), n_samples)
|
|
|
|
# Diurnal temperature: sinusoid between tmin and tmax, peaking at 14:00
|
|
day_indices = np.repeat(np.arange(n_samples), hours_per_day)
|
|
t_range = np.repeat(tmax - tmin, hours_per_day)
|
|
t_phase = 2 * np.pi * (hour_of_day - 14) / 24
|
|
t_hourly = np.repeat(tmin, hours_per_day) + t_range * (0.5 + 0.5 * np.cos(t_phase)) + rng.normal(0, 0.5, total_hours)
|
|
|
|
# RH: inverse of temperature cycle
|
|
rh_hourly = np.repeat(rh_mean, hours_per_day) - 5 * np.cos(t_phase) + rng.normal(0, 3, total_hours)
|
|
rh_hourly = np.clip(rh_hourly, 20, 100)
|
|
|
|
hourly_data = {
|
|
"date": hour_timestamps,
|
|
"temperature_2m": t_hourly,
|
|
"relative_humidity_2m": rh_hourly,
|
|
"dew_point_2m": np.repeat(dewpoint, hours_per_day) + rng.normal(0, 0.5, total_hours),
|
|
"apparent_temperature": t_hourly + rng.exponential(2.0, total_hours),
|
|
"precipitation_probability": np.repeat(precip_prob, hours_per_day) / 24 + rng.normal(0, 1, total_hours),
|
|
"precipitation": np.repeat(rain_sum, hours_per_day) / 24 * rng.uniform(0.5, 1.5, total_hours),
|
|
"rain": np.repeat(rain_sum, hours_per_day) / 24,
|
|
"cloud_cover": np.repeat(cloud_cover, hours_per_day) + rng.normal(0, 5, total_hours),
|
|
"cloud_cover_low": np.repeat(cloud_low, hours_per_day),
|
|
"cloud_cover_mid": np.repeat(cloud_mid, hours_per_day),
|
|
"cloud_cover_high": np.repeat(cloud_high, hours_per_day),
|
|
"wind_speed_10m": np.repeat(wind_speed_max, hours_per_day) * 0.5 * (0.5 + 0.5 * np.cos(t_phase)),
|
|
"wind_speed_100m": np.repeat(wind_speed_100m_max, hours_per_day) * 0.6,
|
|
"wind_gusts_10m": np.repeat(wind_gusts_max, hours_per_day) * (0.3 + 0.7 * rng.beta(2, 5, total_hours)),
|
|
"wind_direction_10m": np.repeat(wind_dir, hours_per_day) + rng.normal(0, 10, total_hours),
|
|
"surface_pressure": np.repeat(pressure, hours_per_day) + rng.normal(0, 0.5, total_hours),
|
|
"visibility": np.repeat(visibility, hours_per_day) + rng.normal(0, 500, total_hours),
|
|
}
|
|
hourly = pd.DataFrame(hourly_data).set_index("date")
|
|
hourly.index = pd.to_datetime(hourly.index)
|
|
|
|
# Clip all values to realistic ranges
|
|
hourly["cloud_cover"] = np.clip(hourly["cloud_cover"], 0, 100)
|
|
hourly["cloud_cover_low"] = np.clip(hourly["cloud_cover_low"], 0, 100)
|
|
hourly["cloud_cover_mid"] = np.clip(hourly["cloud_cover_mid"], 0, 100)
|
|
hourly["cloud_cover_high"] = np.clip(hourly["cloud_cover_high"], 0, 100)
|
|
hourly["visibility"] = np.clip(hourly["visibility"], 100, 30000)
|
|
hourly["precipitation_probability"] = np.clip(hourly["precipitation_probability"], 0, 100)
|
|
|
|
return daily, hourly
|
|
|
|
|
|
def train_targets(
|
|
target_names: Optional[list] = None,
|
|
n_bootstrap: int = 5000,
|
|
test_split: float = 0.2,
|
|
):
|
|
"""Train all or selected target models."""
|
|
if target_names is None:
|
|
target_names = list(TARGET_DEFINITIONS.keys())
|
|
|
|
print(f"Training {len(target_names)} models...")
|
|
print(f"Bootstrap samples: {n_bootstrap} (test split: {test_split:.0%})")
|
|
print()
|
|
|
|
daily, hourly = bootstrap_training_data(n_bootstrap)
|
|
|
|
engine = FeatureEngine()
|
|
X = engine.transform(daily, hourly)
|
|
print(f"Features: {X.shape[1]} from {len(engine.FEATURE_GROUPS)} groups")
|
|
print(f" Thermal: {len(engine.FEATURE_GROUPS['thermal'])}")
|
|
print(f" Dynamic: {len(engine.FEATURE_GROUPS['dynamic'])}")
|
|
print(f" Moisture: {len(engine.FEATURE_GROUPS['moisture'])}")
|
|
print(f" Temporal: {len(engine.FEATURE_GROUPS['temporal'])}")
|
|
print(f" Interaction: {len(engine.FEATURE_GROUPS['interaction'])}")
|
|
print()
|
|
|
|
# Train/test split (temporal order, no shuffle)
|
|
split_idx = int(len(daily) * (1 - test_split))
|
|
X_train, X_test = X[:split_idx], X[split_idx:]
|
|
daily_train, daily_test = daily.iloc[:split_idx], daily.iloc[split_idx:]
|
|
|
|
results = {}
|
|
|
|
# Feature augmentation: add Gaussian noise to prevent overfitting on synthetic data
|
|
X_train_noisy = X_train + np.random.RandomState(42).normal(0, 0.1, X_train.shape).astype(np.float32)
|
|
X_test_noisy = X_test + np.random.RandomState(43).normal(0, 0.05, X_test.shape).astype(np.float32)
|
|
|
|
for target_name in target_names:
|
|
print(f"{'='*60}")
|
|
print(f"Training: {target_name}")
|
|
print(f" {TARGET_DEFINITIONS[target_name]['description']}")
|
|
print(f"{'='*60}")
|
|
|
|
tdef = TARGET_DEFINITIONS[target_name]
|
|
y_train = WeatherModel.build_target(
|
|
daily_train, tdef["variable"], tdef["threshold"], tdef["op"]
|
|
)
|
|
y_test = WeatherModel.build_target(
|
|
daily_test, tdef["variable"], tdef["threshold"], tdef["op"]
|
|
)
|
|
|
|
p_yes = y_train.mean() * 100
|
|
print(f" Class balance: {p_yes:.1f}% YES / {100-p_yes:.1f}% NO")
|
|
|
|
model = WeatherModel(target_name)
|
|
model.train(X_train_noisy, y_train, X_test_noisy, y_test)
|
|
metrics = model.evaluate(X_test, y_test)
|
|
model.save()
|
|
|
|
results[target_name] = metrics
|
|
|
|
print(f" Brier score: {metrics['brier_score']:.4f}")
|
|
print(f" ROC AUC: {metrics['roc_auc']:.3f}")
|
|
print(f" Predicted mean: {metrics['p_yes_predicted']:.1f}% (actual: {metrics['p_yes_actual']:.1f}%)")
|
|
print(f" Top 10 features:")
|
|
for feat, imp in list(model.top_features(10).items()):
|
|
print(f" {feat:30s} {imp:>10.1f}")
|
|
print()
|
|
|
|
# Summary
|
|
print(f"\n{'='*60}")
|
|
print("TRAINING SUMMARY")
|
|
print(f"{'='*60}")
|
|
print(f"{'Target':<25s} {'Brier':>8s} {'ROC AUC':>8s} {'Cal Err %':>10s} {'Samples':>8s}")
|
|
print("-" * 62)
|
|
for name, m in results.items():
|
|
cal_err = abs(m["p_yes_predicted"] - m["p_yes_actual"])
|
|
print(f"{name:<25s} {m['brier_score']:>8.4f} {m['roc_auc']:>8.3f} {cal_err:>10.1f} {m['n_samples']:>8d}")
|
|
|
|
print(f"\nModels saved to: {MODEL_DIR}")
|
|
return results
|
|
|
|
|
|
def main():
|
|
parser = argparse.ArgumentParser(description="Train HK weather prediction models")
|
|
parser.add_argument("--target", type=str, default=None, help="Train single target (e.g., temp_gt_30c_24h)")
|
|
parser.add_argument("--bootstrap", action="store_true", default=True, help="Use bootstrap training data")
|
|
parser.add_argument("--samples", type=int, default=5000, help="Bootstrap sample count")
|
|
parser.add_argument("--all", action="store_true", default=False, help="Train all targets")
|
|
args = parser.parse_args()
|
|
|
|
if args.target:
|
|
targets = [args.target]
|
|
elif args.all:
|
|
targets = list(TARGET_DEFINITIONS.keys())
|
|
else:
|
|
targets = list(TARGET_DEFINITIONS.keys())
|
|
|
|
results = train_targets(targets, n_bootstrap=args.samples)
|
|
|
|
# Save summary
|
|
MODEL_DIR.mkdir(parents=True, exist_ok=True)
|
|
summary_path = MODEL_DIR / "training_summary.json"
|
|
with open(summary_path, "w") as f:
|
|
json.dump({
|
|
"training_date": datetime.now().isoformat(),
|
|
"n_bootstrap_samples": args.samples,
|
|
"results": results,
|
|
}, f, indent=2)
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|