Blame for example.py

"""Load the frozen firepredict model and compute a chance for one commune-day.

Run from this folder:
    uv run --with scikit-learn==1.9.1 --with joblib --with pandas python example.py

The model answers: what is the chance that at least one fire of 0.1 ha or more
starts in this Var commune on the target day? The forecast is issued on the
morning of the day before the target day and uses weather up to two days
before the target day.

The five inputs (one row per commune and target day):

  log_usual_rate  log(usual_rate + 0.01). usual_rate is the commune's average
                  number of days per year with at least one fire start of
                  0.1 ha or more (BDIFF, 2006-2025 without 2024). Take it from
                  usual_rates.csv by INSEE code.
  log_local_fwi   log(1 + FWI). FWI is the Fire Weather Index (no unit) of the
                  day two days before the target day, computed from the
                  commune's nearest Météo-France stations with daily values
                  (maximum temperature, minimum humidity, mean wind, daily
                  rain). It is an approximation, not the official FWI product.
  local_wb_30d    Rain minus potential evaporation (Météo-France ETP), in mm,
                  summed over the 30 days that end two days before the target
                  day, at the commune's nearest stations. Negative = dry.
  orange_period   1 if the target day is in February or March (brush burning
                  allowed under conditions, Var prefectoral orders), else 0.
  red_period      1 if the target day is in June to September (fire use
                  banned), else 0.

A missing input is replaced with the training median stored in
metadata.json ("fill_values"), as the frozen pipeline does.

To build these inputs from raw data (weather download, station matching,
FWI), use the code at https://github.com/guillaumegay13/firepredict
(firepredict.model_communes.build_panel and firepredict.freeze.predict).

This is a research preview, not an official warning tool. Fire danger
warnings in France come from Météo-France and the prefectures.
"""

import hashlib
import json
from pathlib import Path

import joblib
import numpy as np
import pandas as pd
import sklearn

HERE = Path(__file__).parent


def sha256(path: Path) -> str:
    return hashlib.sha256(path.read_bytes()).hexdigest()


# The pipeline was saved with scikit-learn 1.9.1; older versions fail to load it.
if sklearn.__version__ != "1.9.1":
    print(f"warning: scikit-learn {sklearn.__version__} found, the model was saved with 1.9.1")

metadata = json.loads((HERE / "metadata.json").read_text())

# joblib files are pickles: only load files whose fingerprint you trust.
assert sha256(HERE / "model.joblib") == metadata["model_sha256"], "model.joblib does not match metadata.json"
assert sha256(HERE / "forecasts_2026.csv") == metadata["forecasts_2026"]["sha256"], "forecasts_2026.csv changed"

model = joblib.load(HERE / "model.joblib")
rates = pd.read_csv(HERE / "usual_rates.csv", dtype={"insee": str}).set_index("insee").usual_rate


def inputs(insee: str, fwi: float, wb_30d_mm: float, target_month: int) -> pd.DataFrame:
    """One row of the five model inputs, in the order the model expects."""
    row = {
        "log_usual_rate": np.log(rates.get(insee, np.nan) + 0.01),
        "log_local_fwi": np.log1p(fwi),
        "local_wb_30d": wb_30d_mm,
        "orange_period": int(target_month in (2, 3)),
        "red_period": int(target_month in (6, 7, 8, 9)),
    }
    X = pd.DataFrame([row])[metadata["inputs"]]
    return X.fillna(pd.Series(metadata["fill_values"]))


examples = [
    # (description, INSEE code, FWI two days before, 30-day rain - evaporation in mm, month of target day)
    ("Toulon, dry July day", "83137", 35.0, -140.0, 7),
    ("Toulon, February day", "83137", 8.0, -20.0, 2),
    ("Toulon, wet November day", "83137", 2.0, 60.0, 11),
]

print(f"model {metadata['version']}, frozen {metadata['frozen_at']}, inputs {metadata['inputs']}")
print(f"usual rate of Toulon (83137): {rates['83137']:.3f} days with a fire start per year\n")
for label, insee, fwi, wb, month in examples:
    X = inputs(insee, fwi, wb, month)
    chance = model.predict_proba(X)[0, 1]
    # Same result by hand: standardise, weighted sum, logistic function.
    scaler, logit = model[0], model[-1]
    z = ((X.to_numpy()[0] - scaler.mean_) / scaler.scale_) @ logit.coef_[0] + logit.intercept_[0]
    assert np.isclose(chance, 1 / (1 + np.exp(-z)))
    print(f"{label:28s} {X.round(3).to_dict('records')[0]}")
    print(f"{'':28s} chance of a fire start: {chance:.4%}\n")

forecasts = pd.read_csv(HERE / "forecasts_2026.csv", dtype={"insee": str})
print(f"forecasts_2026.csv: {len(forecasts)} rows, {forecasts.target_day.min()} to {forecasts.target_day.max()}")