"""Offline meta-labeling pipeline over the EA's exported candidate datasets.

Input: Common\\Files\\Warrior_EA\\MetaExport\\<SYM>_<PERIOD>.f32 (+ .meta.csv sidecar),
written by CSignalMETA::ExportMetaDataset - one row per resolved, labeled candidate:
    [barTime int64][family int32][pattern int32][side int32][won int32][width float32s]
byte-equivalent to what the EA's pass 2 trains on (same window builder, descriptor, labels).

Discipline mirrors the EA exactly:
  - chronological split per symbol: oldest 55% train | purge | 15% calibration | purge | 30% test,
    purges = horizonBars * period (triple-barrier labels look that far forward -> leakage band);
  - operating point fitted on the CALIBRATION slice only (coverage x (p - BE), the EA's
    FitDirConfThreshold objective, with the EA's 25%-coverage floor);
  - the TEST slice is touched once, by the fitted threshold.

Usage:
  py -3 research/meta_pool.py stats                    # per-file inventory + base rates
  py -3 research/meta_pool.py eval XAUUSD_16388        # single-symbol train+eval
  py -3 research/meta_pool.py pool                     # pooled train, per-symbol eval
"""
import csv
import glob
import os
import sys

import numpy as np

EXPORT_DIR = r"C:\Users\admin\AppData\Roaming\MetaQuotes\Terminal\Common\Files\Warrior_EA\MetaExport"
PERIOD_SECONDS = {16385: 3600, 16386: 7200, 16387: 10800, 16388: 14400, 16390: 21600,
                  16392: 28800, 16396: 43200, 16408: 86400}
MIN_COVERAGE_FRACTION = 0.25   # the EA's MIN_COVERAGE_FRACTION_OF_BASE_RATE against a 100% base
EDGE_MIN_SIGMAS = 2.0          # the EA's deployability floor
SEED = 20260813


def load(tag):
    meta_path = os.path.join(EXPORT_DIR, tag + ".meta.csv")
    # MQL5's FileOpen(FILE_CSV) without FILE_ANSI writes UTF-16LE
    with open(meta_path, encoding="utf-16") as f:
        r = list(csv.DictReader(f))[0]
    width = int(r["width"])
    dt = np.dtype([("t", "<i8"), ("fam", "<i4"), ("pat", "<i4"), ("side", "<i4"),
                   ("won", "<i4"), ("x", "<f4", (width,))])
    # memory-map: 6.5 GB of exports on a box that also hosts the training charts. The
    # export walks candidate ids (sweep order, oldest first), so the sort below is a
    # no-op in practice and the mmap is used as-is; only slices ever materialize.
    raw = np.memmap(os.path.join(EXPORT_DIR, tag + ".f32"), dtype=dt, mode="r")
    if np.any(np.diff(raw["t"]) < 0):
        order = np.argsort(raw["t"], kind="stable")
        raw = np.asarray(raw)[order]
    meta = {k: (float(v) if k not in ("symbol",) else v) for k, v in r.items()}
    meta["period_sec"] = PERIOD_SECONDS[int(float(r["period"]))]
    return raw, meta


def chrono_split(rows, meta):
    """Return boolean masks (train, calib, test) with horizon purges between slices."""
    t = rows["t"]
    purge = int(float(meta["horizonBars"])) * meta["period_sec"]
    n = len(t)
    i_train_end = int(n * 0.55)
    i_calib_end = int(n * 0.70)
    t_train_hi = t[i_train_end - 1]
    t_calib_lo = t_train_hi + purge
    t_calib_hi = t[min(i_calib_end, n - 1)]
    t_test_lo = t_calib_hi + purge
    train = t <= t_train_hi
    calib = (t >= t_calib_lo) & (t <= t_calib_hi)
    test = t >= t_test_lo
    return train, calib, test


class MLP:
    """width -> h1 -> h2 -> 2 softmax, input standardization + Adam + CE. Mirrors the EA's
    head in spirit (dense+norm+bounded scale); small enough to train in seconds."""

    def __init__(self, width, h1=32, h2=16, seed=SEED):
        rng = np.random.default_rng(seed)
        s1 = np.sqrt(2.0 / width)
        s2 = np.sqrt(2.0 / h1)
        s3 = np.sqrt(2.0 / h2)
        self.W1 = rng.normal(0, s1, (width, h1)); self.b1 = np.zeros(h1)
        self.W2 = rng.normal(0, s2, (h1, h2));    self.b2 = np.zeros(h2)
        self.W3 = rng.normal(0, s3, (h2, 2));     self.b3 = np.zeros(2)
        self.mu = None
        self.sd = None

    def _fwd(self, X):
        Z1 = X @ self.W1 + self.b1
        A1 = np.maximum(Z1, 0.01 * Z1)
        Z2 = A1 @ self.W2 + self.b2
        A2 = np.maximum(Z2, 0.01 * Z2)
        Z3 = A2 @ self.W3 + self.b3
        Z3 = Z3 - Z3.max(axis=1, keepdims=True)
        P = np.exp(Z3)
        P /= P.sum(axis=1, keepdims=True)
        return Z1, A1, Z2, A2, P

    def fit(self, X, y, epochs=30, batch=256, lr=1e-3, l2=3e-3, verbose=False):
        # X stays float32 (the pooled set is ~4 GB); standardization and the float64
        # cast happen per mini-batch, where they cost ~2 MB instead of 2x the corpus.
        self.mu = X.mean(axis=0, dtype=np.float64)
        self.sd = X.std(axis=0, dtype=np.float64) + 1e-8
        Y = np.stack([y, 1 - y], axis=1).astype(np.float64)           # slot0 = P(win)
        rng = np.random.default_rng(SEED)
        params = [self.W1, self.b1, self.W2, self.b2, self.W3, self.b3]
        m = [np.zeros_like(p) for p in params]
        v = [np.zeros_like(p) for p in params]
        b1m, b2m, eps, step = 0.9, 0.999, 1e-8, 0
        for ep in range(epochs):
            idx = rng.permutation(len(X))
            for lo in range(0, len(idx), batch):
                sel = idx[lo:lo + batch]
                Xb = (X[sel].astype(np.float64) - self.mu) / self.sd
                Yb = Y[sel]
                Z1, A1, Z2, A2, P = self._fwd(Xb)
                dZ3 = (P - Yb) / len(sel)
                dW3 = A2.T @ dZ3 + l2 * self.W3; db3 = dZ3.sum(0)
                dA2 = dZ3 @ self.W3.T
                dZ2 = dA2 * np.where(Z2 > 0, 1.0, 0.01)
                dW2 = A1.T @ dZ2 + l2 * self.W2; db2 = dZ2.sum(0)
                dA1 = dZ2 @ self.W2.T
                dZ1 = dA1 * np.where(Z1 > 0, 1.0, 0.01)
                dW1 = Xb.T @ dZ1 + l2 * self.W1; db1 = dZ1.sum(0)
                grads = [dW1, db1, dW2, db2, dW3, db3]
                step += 1
                for p, g, mi, vi in zip(params, grads, m, v):
                    mi *= b1m; mi += (1 - b1m) * g
                    vi *= b2m; vi += (1 - b2m) * g * g
                    mh = mi / (1 - b1m ** step)
                    vh = vi / (1 - b2m ** step)
                    p -= lr * mh / (np.sqrt(vh) + eps)
            if verbose and ep % 10 == 9:
                sub = rng.choice(len(X), min(len(X), 200_000), replace=False)
                p = self.p_win(X[sub])
                ce = -np.mean(np.log(np.clip(np.where(y[sub] == 1, p, 1 - p), 1e-9, 1)))
                print(f"  epoch {ep+1}: train CE {ce:.4f} (200k subsample)", flush=True)
        return self

    def p_win(self, X, chunk=65536):
        out = np.empty(len(X))
        for lo in range(0, len(X), chunk):
            Xs = (X[lo:lo + chunk].astype(np.float64) - self.mu) / self.sd
            out[lo:lo + chunk] = self._fwd(Xs)[4][:, 0]
        return out


def fit_threshold(p, won, be_pct):
    """The EA's operating point: maximize coverage x (precision - BE) over the calib slice,
    subject to the 25% coverage floor; ties to the lower threshold."""
    be = be_pct / 100.0
    best = (-np.inf, 0.0)
    for thr in np.linspace(0.0, 0.98, 50):
        sel = p >= thr
        cov = sel.mean()
        if cov < MIN_COVERAGE_FRACTION or sel.sum() == 0:
            continue
        prec = won[sel].mean()
        score = cov * (prec - be)
        if score > best[0]:
            best = (score, thr)
    return best[1]


def evaluate(tag, p_test, rows_test, meta, label=""):
    won = rows_test["won"].astype(np.float64)
    be = float(meta["breakEvenPct"])
    base = 100.0 * won.mean()
    out = [f"{tag}{label}: test n={len(won)}, base win {base:.1f}% vs BE {be:.1f}% (gap {base-be:+.1f}pp)"]
    for thr_label, thr in [("argmax(0.5)", 0.5), ("fitted", None)]:
        t = thr if thr is not None else evaluate.fitted_thr
        sel = p_test >= t
        if sel.sum() == 0:
            out.append(f"  {thr_label:<12} thr {t:.2f}: no trades")
            continue
        cov = 100.0 * sel.mean()
        prec = 100.0 * won[sel].mean()
        n = int(sel.sum())
        se = 100.0 * np.sqrt((base / 100) * (1 - base / 100) / n)
        score = cov * (prec - be) / 100.0
        verdict = "DEPLOYABLE" if (cov >= 25.0 and prec - base > EDGE_MIN_SIGMAS * se
                                   and prec > be) else ""
        out.append(f"  {thr_label:<12} thr {t:.2f}: trades {n} ({cov:.1f}%) win {prec:.1f}% | "
                   f"cov x (p-BE) = {score:+.2f} | skill vs base {prec-base:+.1f}pp "
                   f"(2 sigma = {EDGE_MIN_SIGMAS*se:.1f}pp) {verdict}")
    print("\n".join(out))


def tags_on_disk():
    return sorted(os.path.basename(f)[:-4] for f in glob.glob(os.path.join(EXPORT_DIR, "*.f32")))


def cmd_stats():
    for tag in tags_on_disk():
        rows, meta = load(tag)
        t0 = np.datetime64(int(rows["t"].min()), "s")
        t1 = np.datetime64(int(rows["t"].max()), "s")
        won = rows["won"].astype(float)
        long_mask = rows["side"] > 0
        print(f"{tag}: {len(rows)} rows {t0}..{t1} | base {100*won.mean():.1f}% "
              f"(long {100*won[long_mask].mean():.1f}% short {100*won[~long_mask].mean():.1f}%) "
              f"vs BE {meta['breakEvenPct']:.1f}% | width {int(meta['width'])} "
              f"geom {meta['slMult']:.2f}/{meta['tpMult']:.2f}")


def _standardize(xb, mu, sd, chunk=65536):
    out = np.empty(xb.shape, dtype=np.float32)
    for lo in range(0, len(xb), chunk):
        out[lo:lo + chunk] = (xb[lo:lo + chunk] - mu) / sd
    return out


def run_eval(train_sets, eval_tags, label, per_symbol_norm=False, h1=32, h2=16, l2=3e-3):
    """per_symbol_norm: standardize each symbol's features by ITS OWN train-slice mu/sd
    before pooling - the cross-sectional literature's standard fix for pooled models
    collapsing to the prior when instruments' feature scales differ (the EA's own
    batch-norm does the per-chart equivalent). Fitted on train only - no leakage."""
    norms = {}
    Xs, ys = [], []
    for rows, meta, mask in train_sets:
        xb = rows["x"][mask]                # float32 - cast happens per batch in fit()
        if per_symbol_norm:
            mu = xb.mean(axis=0, dtype=np.float64)
            sd = xb.std(axis=0, dtype=np.float64) + 1e-8
            norms[id(rows)] = (mu, sd)
            xb = _standardize(xb, mu, sd)
        Xs.append(xb)
        ys.append(rows["won"][mask])
    X = np.concatenate(Xs)
    y = np.concatenate(ys)
    del Xs
    print(f"[{label}] training on {len(X)} rows from {len(train_sets)} file(s)...")
    net = MLP(X.shape[1], h1=h1, h2=h2).fit(X, y, l2=l2, verbose=True)
    for tag, rows, meta, calib, test in eval_tags:
        def prep(sl):
            xb = rows["x"][sl]
            if per_symbol_norm:
                mu, sd = norms[id(rows)]
                xb = _standardize(xb, mu, sd)
            return xb
        p_cal = net.p_win(prep(calib))
        evaluate.fitted_thr = fit_threshold(p_cal, rows["won"][calib].astype(float), float(meta["breakEvenPct"]))
        p_test = net.p_win(prep(test))
        evaluate(tag, p_test, rows[test], meta, label=f" [{label}]")
    return net, norms


def cmd_curve(data):
    """Dose-response readout for the pooled-norm model: precision vs threshold on a FIXED
    ladder, calib and test side by side. No fitting, no selection - the question is whether
    the ranking ever crosses each symbol's break-even at any coverage, or flattens. Same
    seed as pool2 -> byte-identical model; this adds no new training attempt."""
    train_sets = [(r, m, tr) for (r, m, tr, _, _) in data.values()]
    eval_tags = [(t, r, m, c, te) for t, (r, m, _, c, te) in data.items()]
    net, norms = run_eval(train_sets, eval_tags, "pooled-norm", per_symbol_norm=True,
                          h1=64, h2=32, l2=1e-3)
    for tag, (rows, meta, _, calib, test) in data.items():
        be = float(meta["breakEvenPct"])
        mu, sd = norms[id(rows)]
        p_cal = net.p_win(_standardize(rows["x"][calib], mu, sd))
        p_tst = net.p_win(_standardize(rows["x"][test], mu, sd))
        w_cal = rows["won"][calib].astype(float)
        w_tst = rows["won"][test].astype(float)
        print(f"\n{tag} threshold curve (BE {be:.1f}%):")
        print(f"  {'thr':>5} | {'cal cov%':>8} {'cal win%':>8} | {'test cov%':>9} {'test n':>7} {'test win%':>9} {'vs BE':>6}")
        for thr in np.arange(0.50, 0.96, 0.05):
            sc, st = p_cal >= thr, p_tst >= thr
            if st.sum() < 200:
                break
            print(f"  {thr:5.2f} | {100*sc.mean():8.1f} {100*w_cal[sc].mean() if sc.sum() else 0:8.1f} | "
                  f"{100*st.mean():9.1f} {int(st.sum()):7d} {100*w_tst[st].mean():9.1f} "
                  f"{100*w_tst[st].mean()-be:+6.1f}")


def main():
    cmd = sys.argv[1] if len(sys.argv) > 1 else "stats"
    if cmd == "stats":
        cmd_stats()
        return
    if cmd == "eval":
        tag = sys.argv[2]
        rows, meta = load(tag)
        train, calib, test = chrono_split(rows, meta)
        print(f"{tag}: split train {train.sum()} / calib {calib.sum()} / test {test.sum()}")
        run_eval([(rows, meta, train)], [(tag, rows, meta, calib, test)], "self")
        return
    if cmd in ("pool", "pool2", "curve"):
        data = {}
        for tag in tags_on_disk():
            rows, meta = load(tag)
            if int(meta["width"]) != 1103:
                print(f"{tag}: width {int(meta['width'])} != 1103 - skipped (pool needs one input shape)")
                continue
            data[tag] = (rows, meta, *chrono_split(rows, meta))
        train_sets = [(r, m, tr) for (r, m, tr, _, _) in data.values()]
        eval_tags = [(t, r, m, c, te) for t, (r, m, _, c, te) in data.items()]
        if cmd == "curve":
            cmd_curve(data)
            return
        if cmd == "pool2":
            # attempt #2 on the same test slices (the naive pooled model underfit to the
            # prior: train CE pinned at base-rate entropy while every self model fit).
            # Per-symbol standardization + capacity 64/32 + l2 1e-3. Verdicts must be
            # read family-wise across BOTH attempts.
            run_eval(train_sets, eval_tags, "pooled-norm", per_symbol_norm=True,
                     h1=64, h2=32, l2=1e-3)
            return
        run_eval(train_sets, eval_tags, "pooled")
        for t, (r, m, tr, c, te) in data.items():
            run_eval([(r, m, tr)], [(t, r, m, c, te)], f"self:{t}")
        return
    print(__doc__)


if __name__ == "__main__":
    main()