import numpy as np, pandas as pd, json
from itertools import product
pd.set_option("display.width", 200)

CAUSES = ["picking_error", "transport_damage", "wrong_address", "stale_stock", "supplier_failure"]
SYMPTOMS = ["damaged_package", "wrong_product", "not_delivered", "delay"]
P_WAREHOUSE = {"zaragoza": 0.60, "getafe": 0.40}
P_CAUSE = {"zaragoza": [0.35, 0.30, 0.15, 0.08, 0.12], "getafe": [0.25, 0.20, 0.12, 0.33, 0.10]}     # CPT from 06-03
P_SYMPTOM = {"damaged_package": [0.10, 0.80, 0.02, 0.02, 0.15],                                      # P(yes | cause)
             "wrong_product":   [0.70, 0.02, 0.01, 0.10, 0.30],
             "not_delivered":   [0.05, 0.10, 0.85, 0.40, 0.30],
             "delay":           [0.10, 0.30, 0.60, 0.90, 0.70]}
MEAN_DAYS = {"picking_error": 3, "transport_damage": 4, "wrong_address": 8, "stale_stock": 9, "supplier_failure": 7}
P_CARRIER = {"zaragoza": {"TransNova": 0.5, "RapidShip": 0.3, "LogiSur": 0.2},
             "getafe":   {"RapidShip": 0.5, "LogiSur": 0.3, "TransNova": 0.2}}

def generate_incidents(n=2000, seed=42, labelled_frac=0.5, p_cause=None):
    """Synthetic incidents.csv consistent with the network from 06-03. 'cause' is the truth; 'recorded_cause' only in labelled_frac."""
    rng = np.random.default_rng(seed); p_cause = p_cause or P_CAUSE            # p_cause: to simulate a different month (drift)
    warehouse = rng.choice(list(P_WAREHOUSE), size=n, p=list(P_WAREHOUSE.values()))
    cause = np.array([rng.choice(CAUSES, p=p_cause[w]) for w in warehouse])
    idx = np.array([CAUSES.index(c) for c in cause])
    df = pd.DataFrame({"incident_id": [f"I{7000 + i}" for i in range(n)], "warehouse": warehouse})
    for s in SYMPTOMS:
        df[s] = (rng.random(n) < np.array(P_SYMPTOM[s])[idx]).astype(int)
    df["days_since_shipping"] = np.clip(rng.poisson([MEAN_DAYS[c] for c in cause]) + 1, 1, 30)
    df["carrier"] = [rng.choice(list(P_CARRIER[w]), p=list(P_CARRIER[w].values())) for w in warehouse]
    df["order_value"] = np.round(rng.gamma(2.0, 60.0, n) + 5, 2)
    df["cause"] = cause
    df["recorded_cause"] = np.where(rng.random(n) < labelled_frac, cause, None)
    return df

inc = generate_incidents(); print(inc.shape, "| labelled:", inc["recorded_cause"].notna().sum())
print(inc.drop(columns=["cause"]).head(3).to_string())
print(inc.groupby("warehouse")["cause"].value_counts(normalize=True).unstack()[CAUSES].round(2).to_string())
print(inc.groupby("cause")[SYMPTOMS].mean().loc[CAUSES].round(2))
print("mean days per cause:", inc.groupby("cause")["days_since_shipping"].mean().round(1).to_dict())
labelled = inc[inc["recorded_cause"].notna()].copy(); y = labelled["recorded_cause"]

(2000, 11) | labelled: 1003
  incident_id warehouse  damaged_package  wrong_product  not_delivered  delay  days_since_shipping    carrier  order_value    recorded_cause
0       I7000    getafe                0              0              0      1                    9  RapidShip        23.26       stale_stock
1       I7001  zaragoza                1              0              0      0                    5  TransNova       100.02  transport_damage
2       I7002    getafe                0              0              1      1                    7  RapidShip       119.75               NaN
cause      picking_error  transport_damage  wrong_address  stale_stock  supplier_failure
warehouse
getafe              0.26              0.19           0.12         0.32              0.11
zaragoza            0.37              0.30           0.14         0.08              0.11
                  damaged_package  wrong_product  not_delivered  delay
cause
picking_error                0.08           0.70           0.06   0.10
transport_damage             0.83           0.01           0.11   0.32
wrong_address                0.02           0.02           0.83   0.58
stale_stock                  0.03           0.11           0.40   0.90
supplier_failure             0.13           0.25           0.30   0.68
mean days per cause: {'picking_error': 4.0, 'stale_stock': 10.2, 'supplier_failure': 8.0, 'transport_damage': 5.1, 'wrong_address': 9.0}

from sklearn.metrics import accuracy_score, f1_score, log_loss, classification_report

def binary_cpt(lst): return {(c,): {1: p, 0: round(1 - p, 4)} for c, p in zip(CAUSES, lst)}
NETWORK = {"warehouse": ([], list(P_WAREHOUSE), {(): P_WAREHOUSE}),
           "cause": (["warehouse"], CAUSES, {(w,): dict(zip(CAUSES, P_CAUSE[w])) for w in P_WAREHOUSE})}
for s in SYMPTOMS: NETWORK[s] = (["cause"], [1, 0], binary_cpt(P_SYMPTOM[s]))
ORDER = ["warehouse", "cause"] + SYMPTOMS

def joint_prob(assignment):
    p = 1.0
    for node in ORDER:
        parents, _, cpt = NETWORK[node]; p *= cpt[tuple(assignment[q] for q in parents)][assignment[node]]
    return p

def infer(query, evidence):
    hidden = [n for n in ORDER if n != query and n not in evidence]; res = {}
    for v in NETWORK[query][1]:
        res[v] = sum(joint_prob(dict(evidence, **{query: v}, **dict(zip(hidden, comb))))
                     for comb in product(*[NETWORK[n][1] for n in hidden]))
    z = sum(res.values()); return {v: p / z for v, p in res.items()}

def bayes_diagnosis(df):
    cache, rows = {}, []
    for _, r in df.iterrows():
        key = (r["warehouse"],) + tuple(int(r[s]) for s in SYMPTOMS)
        if key not in cache:
            cache[key] = infer("cause", dict({"warehouse": r["warehouse"]}, **{s: int(r[s]) for s in SYMPTOMS}))
        rows.append(cache[key])
    return pd.DataFrame(rows, index=df.index)[CAUSES]

CLASSES = sorted(CAUSES)                                   # the alphabetical order scikit-learn uses
P_bayes = bayes_diagnosis(labelled); pred_bayes = P_bayes.idxmax(axis=1)
print("Majority (picking_error):", round((y == "picking_error").mean(), 3))
print("Bayesian network: accuracy", round(accuracy_score(y, pred_bayes), 3), "| macro F1", round(f1_score(y, pred_bayes, average="macro"), 3),
      "| log-loss", round(log_loss(y, P_bayes[CLASSES].values, labels=CLASSES), 3))
print(classification_report(y, pred_bayes, digits=2))
maxp = P_bayes.max(axis=1)
print(f"Mean max P {maxp.mean():.3f} | rows with P ≥ 0.6: {(maxp >= 0.6).mean():.1%}, accuracy {accuracy_score(y[maxp >= 0.6], pred_bayes[maxp >= 0.6]):.3f} | P < 0.6: accuracy {accuracy_score(y[maxp < 0.6], pred_bayes[maxp < 0.6]):.3f}")

Majority (picking_error): 0.314
Bayesian network: accuracy 0.719 | macro F1 0.606 | log-loss 0.789
                  precision    recall  f1-score   support
   picking_error       0.78      0.90      0.83       315
     stale_stock       0.54      0.77      0.64       179
supplier_failure       0.50      0.05      0.09       118
transport_damage       0.88      0.82      0.85       257
   wrong_address       0.62      0.62      0.62       134
        accuracy                           0.72      1003
       macro avg       0.66      0.63      0.61      1003
Mean max P 0.723 | rows with P ≥ 0.6: 67.3%, accuracy 0.822 | P < 0.6: accuracy 0.506

from sklearn.model_selection import StratifiedKFold, cross_val_predict
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import OneHotEncoder, StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, HistGradientBoostingClassifier

CATS, NUMS, BINS = ["warehouse", "carrier"], ["days_since_shipping", "order_value"], SYMPTOMS
COLS = CATS + NUMS + BINS
def build_model(classifier=None, cats=CATS, nums=NUMS, bins=BINS):
    parts = ([("cat", OneHotEncoder(handle_unknown="ignore"), cats)] if cats else []) + \
            ([("num", StandardScaler(), nums)] if nums else []) + ([("bin", "passthrough", bins)] if bins else [])
    return Pipeline([("prep", ColumnTransformer(parts)), ("model", classifier if classifier is not None else LogisticRegression(max_iter=2000))])

cv = StratifiedKFold(5, shuffle=True, random_state=42)
def evaluate_cv(pipe, Xd, name):
    P = cross_val_predict(pipe, Xd, y, cv=cv, method="predict_proba"); pred = np.array(CLASSES)[P.argmax(1)]
    print(f"{name:34s} accuracy {accuracy_score(y, pred):.3f}  macro F1 {f1_score(y, pred, average='macro'):.3f}  log-loss {log_loss(y, P, labels=CLASSES):.3f}")
    return P, pred
models = {"Multiclass logistic": LogisticRegression(max_iter=2000), "Tree (depth 5)": DecisionTreeClassifier(max_depth=5, random_state=42),
          "Forest (300, leaf 5)": RandomForestClassifier(n_estimators=300, min_samples_leaf=5, random_state=42),
          "HGB (depth 3)": HistGradientBoostingClassifier(max_depth=3, learning_rate=0.05, max_iter=200, random_state=42)}
results = {n: evaluate_cv(build_model(m), labelled[COLS], n) for n, m in models.items()}
evaluate_cv(build_model(cats=["warehouse"], nums=[], bins=BINS), labelled[["warehouse"] + BINS], "Logistic ONLY warehouse + symptoms")
evaluate_cv(build_model(cats=[], nums=["days_since_shipping"], bins=[]), labelled[["days_since_shipping"]], "Logistic ONLY days")
P_log, pred_log = results["Multiclass logistic"]
print(pd.crosstab(y, pred_log, rownames=["actual"], colnames=["logistic"]).to_string())
accs = [(accuracy_score(y.iloc[te], pred_bayes.iloc[te]), accuracy_score(y.iloc[te], pred_log[te])) for _, te in cv.split(labelled[COLS], y)]
print("Per fold (bayes, logistic):", [(round(a, 3), round(b, 3)) for a, b in accs], "| mean advantage", round(np.mean([b - a for a, b in accs]), 3))
for n_lab in (100, 250, 500, 1003):
    sub = labelled.sample(n_lab, random_state=1) if n_lab < len(labelled) else labelled
    P = cross_val_predict(build_model(), sub[COLS], sub["recorded_cause"], cv=cv, method="predict_proba")
    print(f"n labelled={n_lab:4d}: logistic {accuracy_score(sub['recorded_cause'], np.array(CLASSES)[P.argmax(1)]):.3f}  Bayesian network {accuracy_score(sub['recorded_cause'], bayes_diagnosis(sub).idxmax(axis=1)):.3f}")

Multiclass logistic                accuracy 0.738  macro F1 0.658  log-loss 0.690
Tree (depth 5)                     accuracy 0.690  macro F1 0.598  log-loss 1.796
Forest (300, leaf 5)               accuracy 0.735  macro F1 0.646  log-loss 0.759
HGB (depth 3)                      accuracy 0.724  macro F1 0.655  log-loss 0.802
Logistic ONLY warehouse + symptoms accuracy 0.709  macro F1 0.612  log-loss 0.804
Logistic ONLY days                 accuracy 0.454  macro F1 0.292  log-loss 1.206
logistic          picking_error  stale_stock  supplier_failure  transport_damage  wrong_address
actual
picking_error               291            3                 9                10              2
stale_stock                   5          127                20                 3             24
supplier_failure             20           33                26                15             24
transport_damage             26            5                12               211              3
wrong_address                 4           28                15                 2             85
Per fold (bayes, logistic): [(0.677, 0.716), (0.731, 0.721), (0.751, 0.761), (0.73, 0.77), (0.705, 0.72)] | mean advantage 0.019
n labelled= 100: logistic 0.780  Bayesian network 0.790
n labelled= 250: logistic 0.752  Bayesian network 0.744
n labelled= 500: logistic 0.764  Bayesian network 0.746
n labelled=1003: logistic 0.738  Bayesian network 0.719

AUTO_THRESHOLD, HIGH_VALUE, MAX_DAYS = 0.60, 300.0, 10
def decide(row, probs):
    cause, p = CLASSES[int(np.argmax(probs))], float(np.max(probs)); actions, trace = [], []
    if row["order_value"] > HIGH_VALUE:
        queue = "human_review"; trace.append(f"R0: value {row['order_value']:.0f} € > {HIGH_VALUE:.0f} € → a person always reviews")
    elif p < AUTO_THRESHOLD:
        queue = "human_review"; trace.append(f"R1: P({cause})={p:.2f} < {AUTO_THRESHOLD} → grey zone")
    else:
        queue = "automatic"; trace.append(f"R2: P({cause})={p:.2f} ≥ {AUTO_THRESHOLD} → cause accepted as a fact")
    if cause == "transport_damage" and p >= AUTO_THRESHOLD:
        actions.append(f"open_trace:{row['carrier']}"); trace.append("R3: transport damage → trace with the carrier")
    if cause == "picking_error" and p >= AUTO_THRESHOLD:
        actions.append(f"reship_from:{row['warehouse']}"); trace.append("R4: picking error → immediate reshipment and notice to the warehouse")
    if row["not_delivered"] == 1 and row["days_since_shipping"] > MAX_DAYS:
        actions.append("carrier_claim"); trace.append(f"R5: not delivered and {row['days_since_shipping']} days > {MAX_DAYS} → claim, whatever the cause")
    if cause == "stale_stock" and p >= AUTO_THRESHOLD:
        actions.append("notify_purchasing"); trace.append("R6: stale stock → notice to purchasing and alternative for the customer")
    return {"queue": queue, "proposed_cause": cause, "confidence": round(p, 3), "actions": actions, "trace": trace}

decisions = pd.DataFrame([decide(r, P_log[i]) for i, (_, r) in enumerate(labelled.iterrows())], index=labelled.index)
labelled2 = pd.concat([labelled, decisions], axis=1)
for queue, g in labelled2.groupby("queue"):
    print(f"{queue:16s} n={len(g):4d} ({len(g) / len(labelled2):.1%})  accuracy of the proposed cause {accuracy_score(g['recorded_cause'], g['proposed_cause']):.3f}")
print("actions triggered:", pd.Series([a.split(':')[0] for l in labelled2["actions"] for a in l]).value_counts().to_dict())
r = labelled2.iloc[6]; print(r["incident_id"], "| actual:", r["recorded_cause"], "→", r["queue"], r["proposed_cause"], r["confidence"], r["actions"]); [print("   ", t) for t in r["trace"]]
for t in (0.5, 0.6, 0.7, 0.8):
    auto = labelled2[(P_log.max(1) >= t) & (labelled2["order_value"] <= HIGH_VALUE)]
    print(f"threshold {t}: automatic {len(auto) / len(labelled2):.1%}, accuracy {accuracy_score(auto['recorded_cause'], auto['proposed_cause']):.3f}")
model = build_model().fit(labelled[COLS], y)                                     # final model with the 1,003 labelled rows
unlabelled = inc[inc["recorded_cause"].isna()].copy(); P_unl = model.predict_proba(unlabelled[COLS])
dec_unl = pd.DataFrame([decide(r, P_unl[i]) for i, (_, r) in enumerate(unlabelled.iterrows())], index=unlabelled.index)
auto = dec_unl["queue"] == "automatic"
print("Unlabelled:", dec_unl["queue"].value_counts().to_dict(), "| true (hidden) accuracy automatic:", round(accuracy_score(unlabelled.loc[auto, "cause"], dec_unl.loc[auto, "proposed_cause"]), 3),
      "| human:", round(accuracy_score(unlabelled.loc[~auto, "cause"], dec_unl.loc[~auto, "proposed_cause"]), 3))

automatic        n= 649 (64.7%)  accuracy of the proposed cause 0.858
human_review     n= 354 (35.3%)  accuracy of the proposed cause 0.517
actions triggered: {'reship_from': 290, 'open_trace': 223, 'notify_purchasing': 104, 'carrier_claim': 67}
I7010 | actual: picking_error → human_review picking_error 0.808 ['reship_from:zaragoza']
    R0: value 349 € > 300 € → a person always reviews
    R4: picking error → immediate reshipment and notice to the warehouse
threshold 0.5: automatic 76.5%, accuracy 0.821
threshold 0.6: automatic 64.7%, accuracy 0.858
threshold 0.7: automatic 53.8%, accuracy 0.906
threshold 0.8: automatic 45.8%, accuracy 0.919
Unlabelled: {'automatic': 697, 'human_review': 300} | true (hidden) accuracy automatic: 0.852 | human: 0.52

import joblib
card = {"name": "NovaMarket incident cause diagnosis (case 9)", "version": "0.1.0-prototype", "date": "2026-08-18",
        "model": "Multiclass logistic regression: warehouse, carrier, 4 symptoms, days since shipping, order value",
        "data": {"incidents": int(len(inc)), "labelled": int(len(labelled)), "classes": CLASSES},
        "metrics_cv5": {"accuracy": 0.738, "f1_macro": 0.658, "baseline_bayesian_network": 0.719, "baseline_majority": 0.314},
        "queues": {"automatic_threshold": AUTO_THRESHOLD, "high_value_review": HIGH_VALUE, "automatic_fraction": 0.647, "automatic_accuracy": 0.858},
        "intended_use": "Propose cause and actions to the agent; decide automatically only in the automatic queue and with reversible actions",
        "out_of_scope_use": "Attributing responsibility to a carrier or employee without human review; using the cause to penalise",
        "limitations": ["supplier_failure is poorly detected (recall 0.22)", "only 1,003 labels, recorded in half of the cases",
                        "synthetic data consistent with the CPTs of 06-03: validate with real data before the pilot"]}
reference = {"symptom_rates": labelled[SYMPTOMS].mean().round(3).to_dict(), "automatic_frac": 0.647,
             "proposed_cause_dist": labelled2["proposed_cause"].value_counts(normalize=True).round(3).to_dict(),
             "days_edges": [0, 2, 4, 6, 8, 10, 14, None], "days_freq": None}
PSI_WARNING, PSI_ALERT, RATE_DIFF = 0.10, 0.25, 0.08
def psi(ref, new, eps=1e-4):
    r, n = np.clip(np.asarray(ref, float), eps, None), np.clip(np.asarray(new, float), eps, None)
    return float(np.sum((n - r) * np.log(n / r)))
def frequencies(values, edges):
    b = [-np.inf if x is None else x for x in edges]; b[-1] = np.inf
    c, _ = np.histogram(values, bins=b); return (c / c.sum()).tolist()
reference["days_freq"] = frequencies(labelled["days_since_shipping"], reference["days_edges"])
json.dump(card, open("model_card_incidents.json", "w"), indent=2, ensure_ascii=False)
json.dump(reference, open("drift_reference_incidents.json", "w"), indent=2); joblib.dump(model, "incidents_model.joblib")

def check_drift(batch, reference, model):
    P = model.predict_proba(batch[COLS]); dec = pd.DataFrame([decide(r, P[i]) for i, (_, r) in enumerate(batch.iterrows())])
    warnings, v = [], psi(reference["days_freq"], frequencies(batch["days_since_shipping"], reference["days_edges"]))
    if v > PSI_WARNING: warnings.append(f"PSI days = {v:.3f}")
    for s, t in reference["symptom_rates"].items():
        if abs(batch[s].mean() - t) > RATE_DIFF: warnings.append(f"rate {s}: {batch[s].mean():.3f} versus {t:.3f}")
    fa = (dec["queue"] == "automatic").mean()
    if abs(fa - reference["automatic_frac"]) > RATE_DIFF: warnings.append(f"automatic fraction {fa:.3f} versus {reference['automatic_frac']:.3f}")
    dc = dec["proposed_cause"].value_counts(normalize=True)
    for c, t in reference["proposed_cause_dist"].items():
        if abs(dc.get(c, 0) - t) > RATE_DIFF: warnings.append(f"proposed cause {c}: {dc.get(c, 0):.3f} versus {t:.3f}")
    level = "ALERT" if v > PSI_ALERT or len(warnings) >= 2 else ("WARNING" if warnings else "OK")
    return level, warnings
print("Normal month:", check_drift(generate_incidents(400, seed=7), reference, model))
p_transport = {"zaragoza": [0.20, 0.55, 0.10, 0.06, 0.09], "getafe": [0.15, 0.45, 0.08, 0.25, 0.07]}   # transport damage shoots up
print("Month with a transport problem:", check_drift(generate_incidents(400, seed=8, p_cause=p_transport), reference, model))

Normal month: ('OK', [])
Month with a transport problem: ('ALERT', ['rate damaged_package: 0.525 versus 0.259', 'rate wrong_product: 0.150 versus 0.282', 'rate not_delivered: 0.177 versus 0.260', 'automatic fraction 0.782 versus 0.647', 'proposed cause picking_error: 0.233 versus 0.345', 'proposed cause transport_damage: 0.497 versus 0.240'])