#!/usr/bin/env python3
"""How good does a forecast have to be before the target scores are reachable?

Run this BEFORE fitting anything. Two independent facts come out of it:

1. The target set is self-contradictory at FBI = 1. With a + b = a + c the false
   alarms equal the misses, so FAR = 1 - POD exactly, and {POD > 0.6, FAR < 0.2,
   FBI ~ 1} silently means POD > 0.8. Any procedure that pins FBI to 1 locks
   itself onto FAR = 1 - POD and can never reach FAR < 0.2.

2. Given the base rate, there is a minimum AUC below which NO threshold on ANY
   score can deliver the required (POD, FAR) pair. If the achievable AUC is far
   below it, the answer is a better information source, not a better model class.
   The bound here assumes binormal ROC with equal variances, which is optimistic:
   a real classifier needs at least this much.
"""

from __future__ import annotations

import sys
import os

import numpy as np
from scipy.stats import norm

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
import cfs_common as C
import build_obs as B

POD_T, POFD_T, FAR_T, CSI_T = 0.60, 0.20, 0.20, 0.40


def far_to_pofd(pod, far, s):
    """POFD implied by (POD, FAR) at base rate s:  b = a(1-FAR)^-1 FAR."""
    a = pod * s
    b = a * far / (1.0 - far)
    return b / (1.0 - s)


def auc_needed(pod, pofd):
    """Binormal equal-variance AUC of the ROC passing through (POFD, POD)."""
    d = norm.ppf(pod) - norm.ppf(pofd)
    return float(norm.cdf(d / np.sqrt(2.0)))


def weekly_rates(tag="train"):
    """Observed EHDE base rate per week window, over the working init sample."""
    ehd, z = B.ehd_of(tag)
    row = {int(o): r for r, o in enumerate(z["ordinal"])}
    inits = C.both_inits()
    out = {}
    for wname, idx in C.WEEKS:
        hits = tot = 0
        for s in inits:
            base = C.to_date(s).toordinal()
            rows = [row.get(base + int(k), -1) for k in idx]
            if min(rows) < 0:
                continue
            hits += int(ehd[rows].any(axis=0).sum())
            tot += ehd.shape[1]
        out[wname] = hits / tot
    return out


def main():
    print("=" * 74)
    print("1. IS THE TARGET SET SELF-CONSISTENT?")
    print("=" * 74)
    print("At FBI = 1: b = c, so FAR = 1 - POD exactly.")
    print("  POD > %.2f  =>  FAR > %.2f, which violates FAR < %.2f."
          % (POD_T, 1 - POD_T, FAR_T))
    print("  FAR < %.2f at FBI = 1 requires POD > %.2f, not %.2f."
          % (FAR_T, 1 - FAR_T, POD_T))
    print("=> Do not pin FBI to 1. Sweep the alert rate r instead and report the")
    print("   under-forecasting (FBI < 1) that buying precision costs.")

    print()
    print("=" * 74)
    print("2. OBSERVED EHDE BASE RATE PER WEEK (541 points, %d inits)"
          % len(C.both_inits()))
    print("=" * 74)
    rates = weekly_rates("train")
    for w in C.WEEK_NAMES:
        print("  %s  %.4f" % (w, rates[w]))

    print()
    print("=" * 74)
    print("3. AUC REQUIRED FOR POD > %.2f WITH FAR < %.2f" % (POD_T, FAR_T))
    print("=" * 74)
    print("%-5s %9s %9s %9s %9s" % ("week", "base", "POFD_req", "AUC_req", "AUC_req*"))
    print("      %s" % "(AUC_req* is for the FBI=1 reading, POD>0.8)")
    for w in C.WEEK_NAMES:
        s = rates[w]
        pofd = far_to_pofd(POD_T, FAR_T, s)
        a1 = auc_needed(POD_T, min(pofd, 0.999))
        pofd2 = far_to_pofd(0.80, FAR_T, s)
        a2 = auc_needed(0.80, min(pofd2, 0.999))
        print("%-5s %9.4f %9.4f %9.4f %9.4f" % (w, s, pofd, a1, a2))

    print()
    print("Read this against the AUC the forecasts actually deliver. If the gap is")
    print("large at long lead, stop tuning thresholds and model families: the")
    print("information is not in the forecast. Release probabilities there instead.")


if __name__ == "__main__":
    main()