Northern-Thailand-Ping-River-Monitor/scripts/forecast_parity.py

#!/usr/bin/env python3
"""Compare ping-collector's forecast with the Python reference on the live database.

src/ml/predict.py is the reference implementation the Rust forecaster
(rust/forecast) was ported from. Run after any change to src/ml or
rust/forecast; scripts/install.sh runs it on every install. On the server, as
the service user from the app dir:

    sudo -u water-monitor .venv/bin/python scripts/forecast_parity.py

Both sides get the same moment and the same Open-Meteo rain series (fetched
once here and handed to the collector as a file), so a difference is a porting
bug, not timing. Exit 0 when every row matches: station, horizon, as_of,
source, model and every level exactly; probabilities within a few units in the
last place (numpy's and libm's exp may differ in the last bit).
"""

import argparse
import datetime
import json
import math
import os
import subprocess
import sys
import tempfile
from pathlib import Path

REPO = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(REPO))

ULPS = 4
EXACT = (
    "station_code",
    "horizon_hours",
    "as_of",
    "source",
    "model_version",
    "trained_at",
    "current_level",
    "predicted_max_level",
    "threshold_warning",
    "threshold_danger",
)


def ulps(a: float, b: float) -> float:
    if a == b:
        return 0
    return abs(a - b) / math.ulp(max(abs(a), abs(b)))


def compare(py: list, rs: list) -> list:
    problems = []
    if len(py) != len(rs):
        problems.append(f"row count: python {len(py)}, collector {len(rs)}")
    for p, r in zip(py, rs):
        where = f"{p.get('station_code')} {p.get('horizon_hours')}h"
        for k in EXACT:
            if p.get(k) != r.get(k):
                problems.append(
                    f"{where} {k}: python {p.get(k)!r}, collector {r.get(k)!r}"
                )
        for k in ("p_warning", "p_danger"):
            d = ulps(float(p[k]), float(r[k]))
            if d > ULPS:
                problems.append(
                    f"{where} {k}: python {p[k]!r}, collector {r[k]!r} ({d:.0f} ulp)"
                )
        ps, rs_ = p.get("stages"), r.get("stages")
        if (ps is None) != (rs_ is None) or (ps and len(ps) != len(rs_)):
            problems.append(f"{where} stages differ")
        elif ps:
            for a, b in zip(ps, rs_):
                if a["stage"] != b["stage"] or a["level"] != b["level"]:
                    problems.append(f"{where} stage {a['stage']} differs")
                elif ulps(a["p_exceed"], b["p_exceed"]) > ULPS:
                    problems.append(
                        f"{where} stage {a['stage']} p_exceed {a['p_exceed']} vs {b['p_exceed']}"
                    )
    return problems


def main() -> int:
    parser = argparse.ArgumentParser(description=__doc__.split("\n\n")[0])
    parser.add_argument("--collector", default=str(REPO / "bin" / "ping-collector"))
    parser.add_argument("--models-dir", default=str(REPO / "models"))
    parser.add_argument(
        "--db-url", default=None, help="default: the app's configuration"
    )
    parser.add_argument("--now", default=None, help="local time, default the clock")
    args = parser.parse_args()

    from src.ml import predict
    from src.ml import rain as rain_mod
    from src.ml.export import json_path_for

    models = Path(args.models_dir)
    stale = [
        p.name
        for p in models.glob("flood_*.joblib")
        if not json_path_for(p).exists()
        or json_path_for(p).stat().st_mtime < p.stat().st_mtime
    ]
    if stale:
        print(
            f"not exported for the collector: {', '.join(stale)}; run `python -m src.ml.export`"
        )
        return 2

    now = (
        datetime.datetime.fromisoformat(args.now)
        if args.now
        else datetime.datetime.now()
    )
    rain = rain_mod.serving_series()
    if rain is None:
        import pandas as pd

        rain = pd.Series(dtype=float)
    py = predict.get_latest_forecasts(
        db_url=args.db_url, models_dir=models, now=now, rain=rain
    )
    py = json.loads(json.dumps(py, default=float))

    with tempfile.TemporaryDirectory() as tmp:
        rain_file = Path(tmp) / "rain.json"
        rain_file.write_text(
            json.dumps(
                {
                    "timestamps": [t.strftime("%Y-%m-%dT%H:%M:%S") for t in rain.index],
                    "values": [None if math.isnan(v) else float(v) for v in rain],
                }
            )
        )
        env = dict(os.environ, MODELS_DIR=str(models), RUST_LOG="error")
        if args.db_url:
            env["POSTGRES_CONNECTION_STRING"] = args.db_url
        out = subprocess.run(
            [
                args.collector,
                "forecast",
                "--json",
                "--now",
                now.isoformat(),
                "--rain-file",
                str(rain_file),
            ],
            capture_output=True,
            text=True,
            env=env,
            cwd=REPO,
            check=False,
        )
    if out.returncode != 0:
        print(f"ping-collector forecast failed:\n{out.stderr}")
        return 2
    rs = json.loads(out.stdout)

    problems = compare(py, rs)
    models_used = sorted({r["model_version"] for r in py})
    n_model = sum(r["source"] == "model" for r in py)
    print(
        f"as of {now:%Y-%m-%d %H:%M:%S}: {len(py)} rows ({n_model} from models {models_used}), "
        f"rain {int(rain.notna().sum())} h"
    )
    if problems:
        print(f"{len(problems)} DIFFERENCES:")
        for p in problems[:40]:
            print("  " + p)
        return 1
    exact = sum(p[k] == r[k] for p, r in zip(py, rs) for k in ("p_warning", "p_danger"))
    print(
        f"identical: every level and field; probabilities {exact}/{2 * len(py)} bit-identical, "
        f"the rest within {ULPS} ulp"
    )
    return 0


if __name__ == "__main__":
    sys.exit(main())