"""
When will this invoice be paid? — a payment-date model for B2B receivables, with a cash forecast.

Data: the HighRadius B2B invoice dataset (50,000 real, anonymised invoices from one company's receivables,
posted Dec 2018 – May 2020): posting/baseline/due dates, payment-terms code, amount, currency, business
unit, masked customer, and the clearing date for the 40,000 that were paid; 10,099 are still open.

Pipeline
  1. load + clean (mixed date encodings, masked customer names → cust_number, exact-duplicate rows dropped)
  2. features known on the posting day; customer history uses ONLY invoices cleared before that day (no lookahead)
  3. time-based split: train on invoices posted before 2019-12-01, test on Dec 2019 – Feb 2020
  4. models: quantile gradient boosting for days-late (p10 / median / p90) + a classifier for "> 7 days late";
     baseline = the customer's own historic median lateness
  5. evaluation: MAE vs baseline, by segment; p10–p90 coverage; classifier AUC / calibration
  6. the WATCH: score the 10,099 open invoices → weekly cash forecast with Monte-Carlo bands vs the due-date schedule
  7. export: trees → model.json; curves, tables, forecast → evaluation.json
"""
import json, math
import numpy as np
import pandas as pd
from sklearn.ensemble import GradientBoostingRegressor, GradientBoostingClassifier
from sklearn.metrics import mean_absolute_error, roc_auc_score, average_precision_score, brier_score_loss, roc_curve
from sklearn.calibration import calibration_curve

SEED = 42
rng = np.random.default_rng(SEED)
SPLIT = pd.Timestamp("2019-12-01")
LATE7 = 7

# ---------------------------------------------------------------- 1. load + clean
raw = pd.read_csv("data/h2h_invoices.csv")
n_raw = len(raw)
raw = raw.drop_duplicates()
def parse_date(s):
    s = s.astype(str).str.replace(r"\.0$", "", regex=True)
    return pd.to_datetime(s, format="%Y%m%d", errors="coerce").fillna(pd.to_datetime(s, errors="coerce"))
for k in ["clear_date", "due_in_date", "posting_date", "baseline_create_date"]:
    raw[k] = parse_date(raw[k])
raw["cust_number"] = raw["cust_number"].astype(str).str.lstrip("0")
raw["amount"] = raw["total_open_amount"].astype(float)
raw["is_open"] = raw["isOpen"].astype(int)
raw["days_late"] = (raw["clear_date"] - raw["due_in_date"]).dt.days
raw = raw.sort_values(["posting_date", "doc_id"]).reset_index(drop=True)
hygiene = {"rows_raw": int(n_raw), "rows_after_dedupe": int(len(raw)), "masked_name_variants": int(raw["name_customer"].nunique()), "customers": int(raw["cust_number"].nunique()),
           "open": int(raw["is_open"].sum()), "cleared": int((raw["is_open"] == 0).sum()), "area_business_all_null": bool(raw["area_business"].isna().all()),
           "baseline_equals_posting_share": float((raw["baseline_create_date"] == raw["posting_date"]).mean()),
           "open_amount_total": float(raw.loc[raw["is_open"] == 1, "amount"].sum())}
print(hygiene)

# ---------------------------------------------------------------- 2. features (history strictly before the posting day)
GLOBAL_MED = float(raw.loc[(raw["is_open"] == 0) & (raw["posting_date"] < SPLIT), "days_late"].median())
hist_cols = ["n_prior", "prior_median_late", "prior_mean_late", "prior_share_late", "prior_share_late7", "last_late", "last3_mean_late", "prior_std_late",
             "days_since_last_clear", "open_count_at_posting", "open_amount_at_posting", "prior_amount_median"]
H = np.full((len(raw), len(hist_cols)), np.nan)
for cust, g in raw.groupby("cust_number", sort=False):
    idx = g.index.values; post = g["posting_date"].values; clear = g["clear_date"].values; dl = g["days_late"].values; amt = g["amount"].values
    for j, i in enumerate(idx):
        t = post[j]
        done = (clear < t) & ~np.isnat(clear)
        n = int(done.sum())
        openm = ~done & (post < t)                      # posted earlier, not yet cleared on this posting day
        row = [n, GLOBAL_MED, GLOBAL_MED, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, int(openm.sum()), float(amt[openm].sum()), np.nan]
        if n:
            d = dl[done]; c = clear[done]
            order = np.argsort(c)
            row[1] = float(np.median(d)); row[2] = float(d.mean()); row[3] = float((d > 0).mean()); row[4] = float((d > LATE7).mean())
            row[5] = float(d[order[-1]]); row[6] = float(d[order[-3:]].mean()); row[7] = float(d.std()) if n > 1 else 0.0
            row[8] = float((t - c.max()) / np.timedelta64(1, "D")); row[11] = float(np.median(amt[done]))
        H[i] = row
hist = pd.DataFrame(H, columns=hist_cols, index=raw.index)

def build_features(d, h):
    f = pd.DataFrame(index=d.index)
    f["log_amount"] = np.log1p(d["amount"])
    f["term_days"] = (d["due_in_date"] - d["baseline_create_date"]).dt.days.clip(0, 120)
    f["baseline_offset"] = (d["baseline_create_date"] - d["posting_date"]).dt.days.clip(-60, 60)
    for b in ["U001", "CA02", "U013"]:
        f["bu_" + b] = (d["business_code"] == b).astype(int)
    f["due_dow"] = d["due_in_date"].dt.dayofweek
    f["due_dom"] = d["due_in_date"].dt.day
    f["due_month"] = d["due_in_date"].dt.month
    f["post_dow"] = d["posting_date"].dt.dayofweek
    for c in hist_cols:
        f[c] = h[c]
    return f

X_all = build_features(raw, hist).fillna(-1)                      # -1 = "no history yet"; n_prior == 0 is the indicator
FEATURES = list(X_all.columns)
cleared = raw["is_open"] == 0
tr = cleared & (raw["posting_date"] < SPLIT)
te = cleared & (raw["posting_date"] >= SPLIT)
X_tr, X_te = X_all[tr], X_all[te]; y_tr, y_te = raw.loc[tr, "days_late"].values, raw.loc[te, "days_late"].values
print("train", tr.sum(), "test", te.sum(), "open", (~cleared).sum())

# ---------------------------------------------------------------- 3/4. models
def gbr(loss, alpha=None):
    kw = dict(n_estimators=220, learning_rate=0.05, max_depth=3, subsample=0.8, min_samples_leaf=40, random_state=SEED, loss=loss)
    if alpha is not None: kw["alpha"] = alpha
    return GradientBoostingRegressor(**kw)
models = {"median": gbr("absolute_error"), "q10": gbr("quantile", 0.10), "q90": gbr("quantile", 0.90)}
for k, m in models.items(): m.fit(X_tr, y_tr)
clf = GradientBoostingClassifier(n_estimators=220, learning_rate=0.05, max_depth=3, subsample=0.8, min_samples_leaf=40, random_state=SEED)
clf.fit(X_tr, (y_tr > LATE7).astype(int))

pred = {k: m.predict(X_te) for k, m in models.items()}
p_late7 = clf.predict_proba(X_te)[:, 1]
y7 = (y_te > LATE7).astype(int)
# baseline: the customer's own prior median (fallback: global median from the training period)
base = X_te["prior_median_late"].fillna(GLOBAL_MED).values
base_term = pd.Series(y_tr).groupby(X_tr["term_days"].values).median()
base2 = X_te["term_days"].map(base_term).fillna(GLOBAL_MED).values
def mae(a, b): return float(mean_absolute_error(a, b))
res = {
    "n_train": int(tr.sum()), "n_test": int(te.sum()), "n_open": int((~cleared).sum()), "split": str(SPLIT.date()),
    "late_share_test": float((y_te > 0).mean()), "late7_share_test": float(y7.mean()),
    "mae": {"customer_median_baseline": mae(y_te, base), "terms_median_baseline": mae(y_te, base2), "zero_baseline": mae(y_te, np.zeros_like(y_te)), "gbm_median": mae(y_te, pred["median"])},
    "coverage_p10_p90": float(((y_te >= pred["q10"]) & (y_te <= pred["q90"])).mean()),
    "band_width_median_days": float(np.median(pred["q90"] - pred["q10"])),
    "late7": {"auc": float(roc_auc_score(y7, p_late7)), "pr_auc": float(average_precision_score(y7, p_late7)), "brier": float(brier_score_loss(y7, p_late7)),
              "auc_baseline_prior_share": float(roc_auc_score(y7, X_te["prior_share_late7"].fillna(0).values))},
}
print(json.dumps(res, indent=1))

# MAE by segment
seg = pd.DataFrame({"y": y_te, "p": pred["median"], "b": base, "n_prior": X_te["n_prior"].values, "amount": raw.loc[te, "amount"].values, "term": X_te["term_days"].values})
def seg_table(col, bins, labels):
    out = []
    for lab, m in zip(labels, [(seg[col] >= lo) & (seg[col] < hi) for lo, hi in bins]):
        if m.sum(): out.append({"segment": lab, "n": int(m.sum()), "mae_gbm": mae(seg.y[m], seg.p[m]), "mae_baseline": mae(seg.y[m], seg.b[m]), "late7_share": float((seg.y[m] > LATE7).mean())})
    return out
mae_by_history = seg_table("n_prior", [(0, 1), (1, 6), (6, 21), (21, 10 ** 9)], ["no history", "1–5 prior", "6–20 prior", "21+ prior"])
mae_by_amount = seg_table("amount", [(0, 5000), (5000, 20000), (20000, 60000), (60000, 10 ** 12)], ["< 5k", "5–20k", "20–60k", "60k+"])

# error distribution: actual vs predicted (sampled points for a scatter) + residual histogram
samp = rng.choice(len(y_te), 600, replace=False)
scatter = [{"y": int(y_te[i]), "p": round(float(pred["median"][i]), 1), "lo": round(float(pred["q10"][i]), 1), "hi": round(float(pred["q90"][i]), 1)} for i in samp]
resid = pred["median"] - y_te
resid_hist = [{"bin": f"{lo}", "n": int(((resid >= lo) & (resid < hi)).sum())} for lo, hi in [(-100, -14), (-14, -7), (-7, -3), (-3, -1), (-1, 1), (1, 3), (3, 7), (7, 14), (14, 100)]]

# classifier curves
fpr, tpr, _ = roc_curve(y7, p_late7)
frac, mean_pred = calibration_curve(y7, p_late7, n_bins=8, strategy="quantile")
dec = pd.qcut(pd.Series(p_late7).rank(method="first"), 10, labels=False)
late7_by_decile = [float(y7[dec == d].mean()) for d in range(9, -1, -1)]

# ---------------------------------------------------------------- 5. descriptive tables (cleared invoices, all periods)
cl = raw[cleared]
dl_hist = [{"bin": lab, "n": int(m.sum()), "amount": float(cl.loc[m, "amount"].sum())} for lab, m in [
    ("early (< −7)", cl.days_late < -7), ("−7 … −1", cl.days_late.between(-7, -1)), ("on the day", cl.days_late == 0), ("1–3", cl.days_late.between(1, 3)),
    ("4–7", cl.days_late.between(4, 7)), ("8–14", cl.days_late.between(8, 14)), ("15–30", cl.days_late.between(15, 30)), ("31–60", cl.days_late.between(31, 60)), ("60+", cl.days_late > 60)]]
tt = cl.assign(term=(cl.due_in_date - cl.baseline_create_date).dt.days)
by_terms = [{"terms": f"net {int(k)}", "n": int(len(g)), "median_late": float(g.days_late.median()), "late7_share": float((g.days_late > LATE7).mean())} for k, g in tt.groupby("term") if len(g) >= 300]
by_bu = [{"bu": k, "n": int(len(g)), "median_late": float(g.days_late.median()), "late7_share": float((g.days_late > LATE7).mean())} for k, g in cl.groupby("business_code") if len(g) >= 100]
top = cl.groupby("cust_number").agg(n=("days_late", "size"), amount=("amount", "sum"), median_late=("days_late", "median"), late7=("days_late", lambda s: float((s > LATE7).mean()))).sort_values("amount", ascending=False).head(15)
top_customers = [{"rank": i + 1, "n": int(r.n), "amount": float(r.amount), "median_late": float(r.median_late), "late7_share": float(r.late7)} for i, r in enumerate(top.itertuples())]
hb = hist[cleared]["prior_share_late7"]
history_predicts = [{"band": lab, "n": int(m.sum()), "late7_share": float((cl.days_late[m] > LATE7).mean())} for lab, m in [
    ("no history", hb.isna()), ("0%", hb == 0), ("0–10%", (hb > 0) & (hb <= .1)), ("10–25%", (hb > .1) & (hb <= .25)), ("25%+", hb > .25)]]
by_month = [{"month": str(k), "n": int(len(g)), "late7_share": float((g.days_late > LATE7).mean()), "median_late": float(g.days_late.median())} for k, g in cl.groupby(cl.posting_date.dt.to_period("M"))]

# ---------------------------------------------------------------- 6. the WATCH: forecast the open book
op = raw[~cleared].copy(); X_op = X_all[~cleared]
for k, m in models.items(): op["pred_" + k] = m.predict(X_op)
op["p_late7"] = clf.predict_proba(X_op)[:, 1]
def sample_days(lo, md, hi, u):
    """piecewise-linear quantile function through (0.1, lo), (0.5, md), (0.9, hi), extended linearly in the tails"""
    lo = np.minimum(lo, md); hi = np.maximum(hi, md)
    out = np.where(u < 0.1, lo - (0.1 - u) / 0.4 * (md - lo), np.where(u < 0.5, lo + (u - 0.1) / 0.4 * (md - lo), np.where(u < 0.9, md + (u - 0.5) / 0.4 * (hi - md), hi + (u - 0.9) / 0.4 * (hi - md))))
    return out
SIMS = 400
week0 = pd.Timestamp("2020-02-24")                       # Monday of the week the open book starts
WEEKS = 26
due_week = ((op["due_in_date"] - week0).dt.days // 7).clip(0, WEEKS - 1).values
amt = op["amount"].values
due_sched = np.bincount(due_week, weights=amt, minlength=WEEKS)
sims = np.zeros((SIMS, WEEKS))
lo, md, hi = op["pred_q10"].values, op["pred_median"].values, op["pred_q90"].values
for s in range(SIMS):
    days = np.round(sample_days(lo, md, hi, rng.uniform(size=len(op))))
    wk = ((op["due_in_date"] + pd.to_timedelta(days, unit="D") - week0).dt.days // 7).clip(0, WEEKS - 1).values
    sims[s] = np.bincount(wk, weights=amt, minlength=WEEKS)
p50_week = op["due_in_date"] + pd.to_timedelta(np.round(md), unit="D")
pred_sched = np.bincount(((p50_week - week0).dt.days // 7).clip(0, WEEKS - 1).values, weights=amt, minlength=WEEKS)
forecast = {"week0": str(week0.date()), "weeks": WEEKS, "due": [float(v) for v in due_sched], "pred_median_date": [float(v) for v in pred_sched],
            "p10": [float(v) for v in np.percentile(sims, 10, axis=0)], "p50": [float(v) for v in np.percentile(sims, 50, axis=0)], "p90": [float(v) for v in np.percentile(sims, 90, axis=0)],
            "cum_due": [float(v) for v in np.cumsum(due_sched)], "cum_p10": [float(v) for v in np.percentile(np.cumsum(sims, axis=1), 10, axis=0)],
            "cum_p50": [float(v) for v in np.percentile(np.cumsum(sims, axis=1), 50, axis=0)], "cum_p90": [float(v) for v in np.percentile(np.cumsum(sims, axis=1), 90, axis=0)],
            "open_total": float(amt.sum()), "expected_late7_amount": float((amt * op["p_late7"].values).sum()), "expected_late7_count": float(op["p_late7"].sum()),
            "sims": SIMS}
# predicted ageing of the open book: amount by predicted lateness bucket (median)
ageing = [{"bucket": lab, "amount": float(amt[m].sum()), "n": int(m.sum())} for lab, m in [("early / on time", md <= 0), ("1–7 days", (md > 0) & (md <= 7)), ("8–30 days", (md > 7) & (md <= 30)), ("30+ days", md > 30)]]
# top open invoices to chase: highest amount × P(late > 7)
chase = op.assign(ev=op["amount"] * op["p_late7"]).sort_values("ev", ascending=False).head(12)
chase_list = [{"amount": float(r.amount), "due": str(r.due_in_date.date()), "p_late7": float(r.p_late7), "pred_days": float(r.pred_median), "lo": float(r.pred_q10), "hi": float(r.pred_q90), "bu": r.business_code, "n_prior": int(X_op.loc[i, "n_prior"])} for i, r in chase.iterrows()]

imp_reg = sorted(zip(FEATURES, models["median"].feature_importances_.tolist()), key=lambda t: -t[1])
imp_clf = sorted(zip(FEATURES, clf.feature_importances_.tolist()), key=lambda t: -t[1])

# ---------------------------------------------------------------- 7. export
def tree_to_json(t):
    return {"left": t.children_left.tolist(), "right": t.children_right.tolist(), "feat": t.feature.tolist(),
            "thr": [float(v) for v in t.threshold.tolist()], "val": [round(float(v[0][0]), 6) for v in t.value.tolist()]}
def export_reg(m): return {"kind": "regressor", "init": float(m.init_.constant_[0][0]), "learning_rate": m.learning_rate, "trees": [tree_to_json(e[0].tree_) for e in m.estimators_]}
model_json = {"features": FEATURES, "global_median_late": GLOBAL_MED, "late_threshold_days": LATE7,
              "models": {k: export_reg(m) for k, m in models.items()},
              "note": "regressors: days_late = init + lr * sum(leaf); classifier: P(late > 7) = sigmoid(init_raw + lr * sum(leaf)). History features are -1 when the customer has no cleared invoice yet (n_prior = 0)."}
model_json["models"]["late7"] = {"kind": "classifier", "init_raw": float(np.log((y_tr > LATE7).mean() / (1 - (y_tr > LATE7).mean()))), "learning_rate": clf.learning_rate, "trees": [tree_to_json(e[0].tree_) for e in clf.estimators_]}
def predict_from_json(mj, row):
    raw_ = mj["init"] if mj["kind"] == "regressor" else mj["init_raw"]
    for t in mj["trees"]:
        node = 0
        while t["left"][node] != -1:
            node = t["left"][node] if row[t["feat"][node]] <= t["thr"][node] else t["right"][node]
        raw_ += mj["learning_rate"] * t["val"][node]
    return raw_ if mj["kind"] == "regressor" else 1 / (1 + math.exp(-raw_))
rows = X_te.values[:300]
for k in ["median", "q10", "q90"]:
    dev = float(np.abs(np.array([predict_from_json(model_json["models"][k], r) for r in rows]) - pred[k][:300]).max()); print(k, "export dev", dev); assert dev < 1e-6
dev = float(np.abs(np.array([predict_from_json(model_json["models"]["late7"], r) for r in rows]) - p_late7[:300]).max()); print("late7 export dev", dev); assert dev < 1e-6

# presets: three real test invoices (features only) + a reference sample for ranking
presets = {}
for name, cond in [("clean_repeat", (X_te["n_prior"] >= 20) & (X_te["prior_share_late7"] == 0)), ("slipping", (X_te["n_prior"] >= 5) & (X_te["last3_mean_late"] > 7)), ("first_invoice", X_te["n_prior"] == 0)]:
    i = X_te[cond].index[0]
    presets[name] = {f: float(X_te.loc[i, f]) for f in FEATURES}
    presets[name].update({"pred_median": float(pred["median"][list(X_te.index).index(i)]), "p_late7": float(p_late7[list(X_te.index).index(i)])})
reference = [{"p": round(float(pred["median"][i]), 1), "a": float(raw.loc[te, "amount"].values[i]), "y": int(y_te[i]), "pl": round(float(p_late7[i]), 4)} for i in rng.choice(len(y_te), 1000, replace=False)]

evaluation = {"results": res, "hygiene": hygiene, "mae_by_history": mae_by_history, "mae_by_amount": mae_by_amount, "scatter": scatter, "resid_hist": resid_hist,
              "late7_roc": {"fpr": [round(float(v), 4) for v in fpr[::8]], "tpr": [round(float(v), 4) for v in tpr[::8]]},
              "late7_calibration": {"mean_pred": [round(float(v), 4) for v in mean_pred], "frac_pos": [round(float(v), 4) for v in frac]}, "late7_by_decile": late7_by_decile,
              "tables": {"days_late_hist": dl_hist, "by_terms": by_terms, "by_bu": by_bu, "top_customers": top_customers, "history_predicts": history_predicts, "by_month": by_month},
              "forecast": forecast, "ageing": ageing, "chase_list": chase_list, "importance_regressor": imp_reg, "importance_classifier": imp_clf, "presets": presets, "reference": reference}
json.dump(model_json, open("model.json", "w")); json.dump(evaluation, open("evaluation.json", "w"))
print("wrote model.json, evaluation.json")
print("mae by history:", mae_by_history); print("history predicts:", history_predicts); print("ageing:", ageing)
print("forecast: open", round(forecast["open_total"] / 1e6, 1), "M; expected >7d late", round(forecast["expected_late7_amount"] / 1e6, 1), "M")
print("importance reg:", imp_reg[:8]); print("importance clf:", imp_clf[:8])