"""Utility condivise per l'analisi PaPP v2.

I dati di input sono i due CSV prodotti dallo script MQL5 `PaPP_CrossExport.mq5`:
  - PaPP_crosses_<SYM>_D1.csv  (un record per incrocio)
  - PaPP_bars_<SYM>_D1.csv     (un record per giorno D1 = baseline)

Per default vengono cercati in ../data. Si possono passare percorsi diversi via argv.
"""
import os, sys, numpy as np, pandas as pd

HERE = os.path.dirname(os.path.abspath(__file__))
DATA = os.path.normpath(os.path.join(HERE, "..", "data"))
HZ = [1, 3, 5, 10, 20]


def load(path):
    df = pd.read_csv(path)
    df["time"] = pd.to_datetime(df["time"], format="%Y.%m.%d", errors="coerce")
    if df["time"].isna().all():
        df["time"] = pd.to_datetime(df["time"], errors="coerce")
    return df.dropna(subset=["time"]).sort_values("time").reset_index(drop=True)


def load_pair(sym="EURUSD", year_min=1999, data_dir=None):
    d = data_dir or DATA
    cr = load(os.path.join(d, f"PaPP_crosses_{sym}_D1.csv"))
    ba = load(os.path.join(d, f"PaPP_bars_{sym}_D1.csv"))
    if year_min:
        cr = cr[cr.time.dt.year >= year_min].reset_index(drop=True)
        ba = ba[ba.time.dt.year >= year_min].reset_index(drop=True)
    return cr, ba


def add_regime(df, edges):
    """Regime = trend(sopra/sotto MA365) x tercile(cluster) x tercile(velocita')."""
    tr = (df["trend"] > 0).astype(int)
    cl = np.digitize(df["cluster_pct"], edges["cl"])
    ve = np.digitize(df["vel_med"], edges["ve"])
    out = df.copy()
    out["regime"] = tr.astype(str) + "_" + cl.astype(str) + "_" + ve.astype(str)
    return out


def regime_edges(ba):
    return {
        "cl": np.quantile(ba["cluster_pct"], [1/3, 2/3]),
        "ve": np.quantile(ba["vel_med"], [1/3, 2/3]),
    }


def block_bootstrap_p(x, nb=2000, bs=20, seed=42):
    """p-value bilaterale per media != 0 con block bootstrap (gestisce overlap)."""
    rng = np.random.default_rng(seed)
    x = np.asarray(x, float)
    n = len(x)
    if n < 10:
        return np.nan
    nblocks = int(np.ceil(n / bs))
    idx = np.arange(n)
    means = np.empty(nb)
    for k in range(nb):
        starts = rng.integers(0, n, size=nblocks)
        pick = np.concatenate([
            idx[s:s+bs] if s+bs <= n else np.concatenate([idx[s:], idx[:s+bs-n]])
            for s in starts])[:n]
        means[k] = x[pick].mean()
    return 2 * min((means <= 0).mean(), (means >= 0).mean())


def benjamini_hochberg(p):
    """q-value BH per array di p (NaN ammessi)."""
    p = np.asarray(p, float)
    out = np.full_like(p, np.nan)
    idx = np.where(~np.isnan(p))[0]
    if len(idx) == 0:
        return out
    m = len(idx)
    order = idx[np.argsort(p[idx])]
    q = p[order] * m / np.arange(1, m + 1)
    q = np.minimum.accumulate(q[::-1])[::-1]
    out[order] = np.clip(q, 0, 1)
    return out