#!/usr/bin/env python3
"""P(Trendumkehr) — kalibriertes Modell nach dem Vorbild von P(break).

User-Frage 2026-07-31: „können wir die Trendumkehr-Wahrscheinlichkeit berechnen?"

WARUM DIESER AUFBAU: Das einzige belegt-prädiktive Modell im Projekt ist P(break)
(AUC 0,71 out-of-sample, Kalibrierung trifft). Es gelang, weil das Ereignis
**mechanisch scharf** definiert war — ein Barriere-Rennen (Level ± 0,5×ATR binnen
12 Bars), nicht „der Kurs bricht aus". Genau daran scheitern Umkehr-Modelle sonst.

EREIGNIS-DEFINITION (dieselbe Bauart):
  Trendrichtung d = sign(EMA12 − EMA50) mit Totband.
  Ab Bar i läuft ein Rennen über die nächsten _W Bars:
      erreicht der Kurs zuerst  −_REV×ATR GEGEN d   → UMKEHR (y=1)
      erreicht er zuerst        +_REV×ATR MIT   d   → Fortsetzung (y=0)
      keins von beidem binnen _W                    → Fortsetzung (y=0, wie im
                                                      P(break)-Training beim Timeout)
  Symmetrisch → die Basisrate sagt direkt, ob der Trend überhaupt persistiert.

MERKMALE (alle kausal aus Bars <= i, ATR-normiert):
  stretch    Überdehnung in Trendrichtung   (C − EMA50)/ATR × d
  angle      Regressionswinkel × d          (dreht der Winkel schon?)
  mom3/mom6  kurzes/mittleres Momentum × d
  strength   Trendstärke |EMA12−EMA50|/ATR
  age        Bars seit dem letzten EMA-Kreuz (gedeckelt/log-skaliert)
  gap        Abstand zum nächsten M30-Gegenlevel in ATR
  volr       ATR relativ zum eigenen Mittel (Vola-Regime)

METHODE: Logistische Regression (Gradientenabstieg auf standardisierten Merkmalen),
Fit auf H1, **AUC + Kalibrierung gemessen auf H2**. Ein Modell wird nur dann
empfohlen, wenn es out-of-sample trennt UND kalibriert ist — sonst ist es eine
weitere unbelegte Zahl (die Lehre aus `analyze_verdict_calibration.py`, wo conf_pct
in H2 invertiert war).

⚠ Erwartungsmanagement: Antizyklisch am Extrem ist im Projekt 6× als Münzwurf
gemessen, die Bounce-Anzeige hat ~50 % Richtungstreffer. Ein AUC nahe 0,5 wäre also
die BESTÄTIGUNG des bisherigen Befunds, kein Rückschlag.
"""
import sys
from bisect import bisect_left
from math import exp, log

import MetaTrader5 as mt5

_W = 24              # Renn-Fenster in M5-Bars (2 h)
_REV = 1.0           # Barriere in ×ATR (symmetrisch)
_ATRMIN = 0.06
_DEAD = 0.15         # Totband für die Trendrichtung (wie _TREND_DEADBAND)
_ANGLE_LR = 20
_M30 = 6
_PIV_K = 3
_PIV_LOOKBACK = 50
_FEATS = ("stretch", "angle", "mom3", "mom6", "strength", "age", "gap", "volr")


def _ema(v, p):
    k = 2.0 / (p + 1); o = []; e = v[0]
    for i, x in enumerate(v):
        e = x if i == 0 else x * k + e * (1 - k)
        o.append(e)
    return o


def _atr_series(H, L, C, p=14):
    t = [0.0]
    for i in range(1, len(C)):
        t.append(max(H[i] - L[i], abs(H[i] - C[i - 1]), abs(L[i] - C[i - 1])))
    out = []
    for i in range(len(C)):
        if not i:
            out.append(None); continue
        w = t[max(1, i - p + 1):i + 1]
        out.append(sum(w) / max(1, len(w)))
    return out


def _slope(ys):
    """Steigung einer LS-Regression über ys (x = 0..n-1)."""
    n = len(ys)
    if n < 3:
        return 0.0
    mx = (n - 1) / 2.0
    my = sum(ys) / n
    num = sum((i - mx) * (y - my) for i, y in enumerate(ys))
    den = sum((i - mx) ** 2 for i in range(n))
    return num / den if den else 0.0


def build_levels(H, L):
    from collections import deque
    nm = len(H) // _M30
    mh = [max(H[i * _M30:(i + 1) * _M30]) for i in range(nm)]
    ml = [min(L[i * _M30:(i + 1) * _M30]) for i in range(nm)]
    born = [[] for _ in range(nm)]
    for p in range(_PIV_K, nm - _PIV_K):
        if mh[p] == max(mh[p - _PIV_K:p + _PIV_K + 1]) and p + _PIV_K < nm:
            born[p + _PIV_K].append(mh[p])
        if ml[p] == min(ml[p - _PIV_K:p + _PIV_K + 1]) and p + _PIV_K < nm:
            born[p + _PIV_K].append(ml[p])
    lv = []; win = deque(); cur = []
    for m in range(nm):
        win.append(born[m]); cur.extend(born[m])
        while len(win) > _PIV_LOOKBACK:
            for x in win.popleft():
                try: cur.remove(x)
                except ValueError: pass
        lv.append(sorted(cur))
    return lv


def samples(O, H, L, C, AT, EF, ES, LV, a, b, step=2):
    """→ [(feature-dict, y)] — kausal: Merkmale aus <=i, Ausgang aus >i."""
    out = []
    atr_mean = None
    age = 0; last_d = 0
    # `age` muss lückenlos mitlaufen → über den ganzen Bereich iterieren,
    # aber nur jeden `step`-ten Bar als Stichprobe nehmen.
    for i in range(max(a, _ANGLE_LR + 60), min(b, len(C) - _W - 1)):
        atr = AT[i]
        if not atr:
            continue
        atr = max(atr, _ATRMIN)
        diff = EF[i] - ES[i]
        d = 1 if diff > _DEAD * atr else -1 if diff < -_DEAD * atr else 0
        if d != last_d and d != 0:
            age = 0; last_d = d
        else:
            age += 1
        if d == 0:
            continue                      # kein Trend → keine „Umkehr" definierbar
        if (i - a) % step:
            continue
        # ── Ausgang: Barriere-Rennen
        up = C[i] + _REV * atr
        dn = C[i] - _REV * atr
        hit_rev = hit_con = None
        for j in range(i + 1, min(i + 1 + _W, len(C))):
            if d > 0:
                if L[j] <= dn: hit_rev = j; break
                if H[j] >= up: hit_con = j; break
            else:
                if H[j] >= up: hit_rev = j; break
                if L[j] <= dn: hit_con = j; break
        y = 1 if hit_rev is not None else 0
        # ── Merkmale
        w = C[max(0, i - _ANGLE_LR):i + 1]
        ang = _slope(w) / atr * _ANGLE_LR          # normierte Steigung
        atr_w = [x for x in AT[max(1, i - 200):i + 1] if x]
        atr_mean = (sum(atr_w) / len(atr_w)) if atr_w else atr
        lv = LV[min(i // _M30, len(LV) - 1)]
        gap = 3.0
        if lv:
            k = bisect_left(lv, C[i])
            # Gegenlevel = das, auf das der Trend zuläuft
            tgt = None
            if d > 0 and k < len(lv):
                tgt = lv[k]
            elif d < 0 and k > 0:
                tgt = lv[k - 1]
            if tgt is not None:
                gap = min(3.0, abs(tgt - C[i]) / atr)
        f = {
            "stretch": (C[i] - ES[i]) / atr * d,
            "angle": ang * d,
            "mom3": (C[i] - C[i - 3]) / atr * d,
            "mom6": (C[i] - C[i - 6]) / atr * d,
            "strength": abs(diff) / atr,
            "age": log(1 + min(age, 300)),
            "gap": gap,
            "volr": atr / atr_mean if atr_mean else 1.0,
        }
        out.append((f, y))
    return out


def fit(rows, iters=400, lr=0.35):
    """Logistische Regression auf standardisierten Merkmalen → (mu, sd, w)."""
    n = len(rows)
    mu = {}; sd = {}
    for k in _FEATS:
        xs = [r[0][k] for r in rows]
        m = sum(xs) / n
        v = sum((x - m) ** 2 for x in xs) / n
        mu[k] = m; sd[k] = (v ** 0.5) or 1.0
    X = [[(r[0][k] - mu[k]) / sd[k] for k in _FEATS] for r in rows]
    Y = [r[1] for r in rows]
    w = [0.0] * (len(_FEATS) + 1)          # w[0] = Bias
    for _ in range(iters):
        g = [0.0] * len(w)
        for xi, yi in zip(X, Y):
            z = w[0] + sum(wk * xk for wk, xk in zip(w[1:], xi))
            p = 1.0 / (1.0 + exp(-max(-30, min(30, z))))
            e = p - yi
            g[0] += e
            for t in range(len(xi)):
                g[t + 1] += e * xi[t]
        for t in range(len(w)):
            w[t] -= lr * g[t] / n
    return mu, sd, w


def predict(f, mu, sd, w):
    z = w[0] + sum(w[t + 1] * ((f[k] - mu[k]) / sd[k]) for t, k in enumerate(_FEATS))
    return 1.0 / (1.0 + exp(-max(-30, min(30, z))))


def auc(ps, ys):
    pairs = sorted(zip(ps, ys))
    pos = sum(ys); neg = len(ys) - pos
    if not pos or not neg:
        return 0.5
    rank = 0.0; i = 0
    while i < len(pairs):
        j = i
        while j + 1 < len(pairs) and pairs[j + 1][0] == pairs[i][0]:
            j += 1
        r = (i + j) / 2.0 + 1
        for k in range(i, j + 1):
            if pairs[k][1] == 1:
                rank += r
        i = j + 1
    return (rank - pos * (pos + 1) / 2) / (pos * neg)


def main():
    global _REV, _W
    n = int(sys.argv[1]) if len(sys.argv) > 1 else 80000
    # Ereignis-Definition variabel — ein Null-Ergebnis ist nur dann belastbar, wenn
    # es nicht an EINER willkürlichen Schwelle hängt.
    if len(sys.argv) > 2:
        _REV = float(sys.argv[2])
    if len(sys.argv) > 3:
        _W = int(sys.argv[3])
    mt5.initialize(); sym = None
    for c in ("SpotCrude", "USOIL", "WTI", "XTIUSD"):
        if mt5.symbol_info(c):
            sym = c; break
    bars = mt5.copy_rates_from_pos(sym, mt5.TIMEFRAME_M5, 0, n)
    mt5.shutdown()
    O = [float(b["open"]) for b in bars]; H = [float(b["high"]) for b in bars]
    L = [float(b["low"]) for b in bars]; C = [float(b["close"]) for b in bars]
    AT = _atr_series(H, L, C); EF = _ema(C, 12); ES = _ema(C, 50)
    LV = build_levels(H, L)
    N = len(C); mid = N // 2

    print("=" * 96)
    print(f"  P(TRENDUMKEHR) — {sym} M5 ({N} Bars)")
    print(f"  Ereignis: erreicht der Kurs zuerst −{_REV}×ATR GEGEN den Trend (Umkehr)")
    print(f"            oder +{_REV}×ATR MIT dem Trend, binnen {_W} Bars = {_W*5} min?")
    print("  Fit auf H1, AUC + Kalibrierung gemessen auf H2 (out-of-sample).")
    print("=" * 96)

    tr = samples(O, H, L, C, AT, EF, ES, LV, 100, mid)
    te = samples(O, H, L, C, AT, EF, ES, LV, mid, N)
    if not tr or not te:
        print("zu wenig Daten"); return
    br_tr = sum(y for _, y in tr) / len(tr)
    br_te = sum(y for _, y in te) / len(te)
    print(f"\n  Stichprobe   H1 n={len(tr)}  Basisrate Umkehr {100*br_tr:.1f}%")
    print(f"               H2 n={len(te)}  Basisrate Umkehr {100*br_te:.1f}%")
    print("  (50 % = der Trend persistiert gar nicht; <50 % = Trendfolge trägt)")

    mu, sd, w = fit(tr)
    print("\n  Gewichte (standardisiert, + = erhöht P(Umkehr)):")
    print(f"    {'bias':<10}{w[0]:+.3f}")
    for t, k in enumerate(_FEATS):
        print(f"    {k:<10}{w[t+1]:+.3f}")

    ps = [predict(f, mu, sd, w) for f, _ in te]
    ys = [y for _, y in te]
    a = auc(ps, ys)
    print(f"\n  AUC out-of-sample (H2): {a:.3f}   "
          f"(0,50 = Münzwurf · P(break) erreichte 0,71)")

    print("\n  Kalibrierung H2 — trifft die vorhergesagte Rate die echte?")
    print(f"    {'Bucket':<14}{'n':>7}{'vorhergesagt':>14}{'real':>9}")
    edges = [0.0, 0.35, 0.45, 0.50, 0.55, 0.65, 1.01]
    for lo, hi in zip(edges, edges[1:]):
        sel = [(p, y) for p, y in zip(ps, ys) if lo <= p < hi]
        if len(sel) < 50:
            continue
        pm = sum(p for p, _ in sel) / len(sel)
        rm = sum(y for _, y in sel) / len(sel)
        print(f"    {lo:.2f}–{hi:<9.2f}{len(sel):>7}{100*pm:>13.1f}%{100*rm:>8.1f}%")

    # In-sample-Kontrolle: trennt es wenigstens auf den Trainingsdaten?
    ps_tr = [predict(f, mu, sd, w) for f, _ in tr]
    print(f"\n  (Kontrolle) AUC in-sample H1: {auc(ps_tr, [y for _, y in tr]):.3f}")

    print("\n  Einzelmerkmale — AUC allein (H2), zeigt wer überhaupt trennt:")
    for k in _FEATS:
        xs = [f[k] for f, _ in te]
        ak = auc(xs, ys)
        print(f"    {k:<10}{ak:.3f}" + ("   <-- traegt" if abs(ak - 0.5) >= 0.03 else ""))


if __name__ == "__main__":
    main()