Blame for scripts/analysis/check_forecast_lead5_anomaly.py

"""
Direct check of the raw previous_runs_forecast.csv archive, independent
of the model/evaluation pipeline -- built specifically because
evaluate_flood_detection.py's --dump-fp-details showed something
systematic and graph-wide at the 5-day horizon (185/185 real test
examples produced at least one real false positive there, often
several stations simultaneously, including stations on different
rivers with no surface connection). --check-forecast-fp-correlation
already ruled out forecast AVAILABILITY as the driver; this checks
whether the real forecast VALUES themselves are anomalous specifically
at lead time 5, compared to the other real lead times in the same file.

Usage:
    python -m scripts.analysis.check_forecast_lead5_anomaly --data-root datasets
"""
import argparse
from pathlib import Path

import numpy as np
import pandas as pd

FORECAST_LEAD_TIMES = [1, 2, 3, 5, 7]  # matches dynamic_features.py's real, confirmed coverage


def summarize_lead_time_column(df: pd.DataFrame, lead_time: int) -> dict:
    """
    Real summary statistics for one real precipitation_previous_dayN
    column -- real coverage, real central tendency, real spread, and
    explicit counts of two real, concrete anomaly types (negative
    precipitation, which is physically impossible, and implausibly
    extreme values) rather than just eyeballing a mean/max pair.
    """
    col = f"precipitation_previous_day{lead_time}"
    if col not in df.columns:
        return {"lead_time": lead_time, "column_present": False}

    values = df[col].dropna()
    n_negative = int((values < 0).sum())
    # A real, extreme daily precipitation threshold -- the wettest days
    # ever recorded anywhere are on the order of a few hundred mm; used
    # here only to flag values worth a second look, not as a hard
    # physical cutoff.
    n_extreme = int((values > 300).sum())

    return {
        "lead_time": lead_time,
        "column_present": True,
        "n_real": len(values),
        "mean": float(values.mean()) if len(values) else float("nan"),
        "median": float(values.median()) if len(values) else float("nan"),
        "std": float(values.std()) if len(values) else float("nan"),
        "min": float(values.min()) if len(values) else float("nan"),
        "max": float(values.max()) if len(values) else float("nan"),
        "n_negative": n_negative,
        "n_extreme_gt_300mm": n_extreme,
    }


def main() -> None:
    parser = argparse.ArgumentParser(description="Check the real ECMWF forecast archive for a lead-time-5-specific anomaly")
    parser.add_argument("--data-root", type=Path, default=Path("datasets"))
    parser.add_argument("--forecast-path", type=Path, default=None,
                         help="Defaults to <data-root>/previous_runs_forecast.csv")
    args = parser.parse_args()

    forecast_path = args.forecast_path or (args.data_root / "previous_runs_forecast.csv")
    if not forecast_path.exists():
        print(f"No real forecast archive found at {forecast_path}.")
        return

    df = pd.read_csv(forecast_path)
    print(f"Loaded {len(df)} real row(s) from {forecast_path}")
    print()

    rows = []
    for lead_time in FORECAST_LEAD_TIMES:
        summary = summarize_lead_time_column(df, lead_time)
        rows.append(summary)

    result_df = pd.DataFrame(rows)
    print("=" * 70)
    print("Real forecast precipitation, by lead time (mm)")
    print("=" * 70)
    print(result_df.to_string(index=False))

    present = result_df[result_df["column_present"]]
    if len(present) < 2:
        print("\nFewer than 2 real lead-time columns present -- nothing to compare.")
        return

    if 5 in present["lead_time"].values:
        lead5 = present[present["lead_time"] == 5].iloc[0]
        others = present[present["lead_time"] != 5]
        others_mean_of_means = others["mean"].mean()
        others_mean_of_maxes = others["max"].mean()
        print(f"\nLead time 5 mean ({lead5['mean']:.2f}mm) vs. real average of other lead times' means "
              f"({others_mean_of_means:.2f}mm)")
        print(f"Lead time 5 max ({lead5['max']:.2f}mm) vs. real average of other lead times' maxes "
              f"({others_mean_of_maxes:.2f}mm)")
        print(f"Lead time 5 negative values: {int(lead5['n_negative'])}; extreme (>300mm) values: "
              f"{int(lead5['n_extreme_gt_300mm'])}")

        mean_ratio = lead5["mean"] / others_mean_of_means if others_mean_of_means > 0 else float("inf")
        max_ratio = lead5["max"] / others_mean_of_maxes if others_mean_of_maxes > 0 else float("inf")
        if mean_ratio > 2.0 or max_ratio > 2.0 or lead5["n_negative"] > 0:
            print(f"\nLead time 5 shows a real, meaningful anomaly relative to the other real lead times "
                  f"(mean ratio {mean_ratio:.2f}x, max ratio {max_ratio:.2f}x) -- real, direct evidence "
                  f"worth investigating further as the actual cause of the day-5 false-positive pattern.")
        else:
            print(f"\nLead time 5 does NOT show a meaningful anomaly relative to the other real lead "
                  f"times (mean ratio {mean_ratio:.2f}x, max ratio {max_ratio:.2f}x) -- the raw forecast "
                  f"values themselves aren't the explanation; the real cause likely lies in how this "
                  f"lead time's data is being consumed downstream (e.g. the anchor/target-date mapping, "
                  f"or the model's own horizon-conditioning for h=5 specifically), not the raw archive.")
    else:
        print("\nLead time 5 column not present in this real file -- can't run the real comparison.")


if __name__ == "__main__":
    main()