"""Statistical validation suite — faithful Python port of the Rust ml-validation/statistical.rs (Bailey &
López de Prado). The deterministic parts (PPF/CDF/DSR/PBO) match the Rust bit-for-bit; the permutation
test matches the algorithm (its RNG differs, so its p-value is tested by property, not bit-equality).

This is the platform's honesty backbone — a richer Deflated Sharpe (skew/kurtosis-adjusted standard
error) plus Probability of Backtest Overfitting (CSCV) and a Monte-Carlo permutation test.
"""
from __future__ import annotations

import itertools
import math
from dataclasses import dataclass

import numpy as np

_EULER = 0.5772156649015329


def normal_cdf(x: float) -> float:
    """Standard-normal CDF (Abramowitz & Stegun 26.2.17), the exact Rust approximation."""
    if x <= -8.0:
        return 0.0
    if x >= 8.0:
        return 1.0
    b1, b2, b3, b4, b5 = 0.31938153, -0.356563782, 1.781477937, -1.821255978, 1.330274429
    ax = abs(x)
    t = 1.0 / (1.0 + 0.2316419 * ax)
    poly = b1 * t + b2 * t**2 + b3 * t**3 + b4 * t**4 + b5 * t**5
    pdf = math.exp(-0.5 * ax * ax) / math.sqrt(2.0 * math.pi)
    cdf_right = 1.0 - pdf * poly
    return cdf_right if x >= 0.0 else 1.0 - cdf_right


def normal_ppf(p: float) -> float:
    """Inverse standard-normal CDF — Peter Acklam's rational approximation (the exact Rust port)."""
    if p <= 0.0:
        return -8.0
    if p >= 1.0:
        return 8.0
    if abs(p - 0.5) < 1e-15:
        return 0.0
    a = (-3.969683028665376e1, 2.209460984245205e2, -2.759285104469687e2,
         1.38357751867269e2, -3.066479806614716e1, 2.506628277459239)
    b = (-5.447609879822406e1, 1.615858368580409e2, -1.556989798598866e2,
         6.680131188771972e1, -1.328068155288572e1)
    c = (-7.784894002430293e-3, -3.223964580411365e-1, -2.400758277161838,
         -2.549732539343734, 4.374664141464968, 2.938163982698783)
    d = (7.784695709041462e-3, 3.224671290700398e-1, 2.445134137142996, 3.754408661907416)
    p_low = 0.02425
    if p < p_low:
        q = math.sqrt(-2.0 * math.log(p))
        num = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]
        den = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0
        return num / den
    if p <= 1.0 - p_low:
        q = p - 0.5
        r = q * q
        num = (((((a[0]*r + a[1])*r + a[2])*r + a[3])*r + a[4])*r + a[5])*q
        den = ((((b[0]*r + b[1])*r + b[2])*r + b[3])*r + b[4])*r + 1.0
        return num / den
    q = math.sqrt(-2.0 * math.log(1.0 - p))
    num = ((((c[0]*q + c[1])*q + c[2])*q + c[3])*q + c[4])*q + c[5]
    den = (((d[0]*q + d[1])*q + d[2])*q + d[3])*q + 1.0
    return -(num / den)


def sharpe_ratio(returns: np.ndarray) -> float:
    """Per-period Sharpe = mean / sample-std (ddof=1)."""
    r = np.asarray(returns, float)
    if r.size < 2:
        return 0.0
    std = float(r.std(ddof=1))
    return float(r.mean()) / std if std >= 1e-15 else 0.0


def skewness(returns: np.ndarray) -> float:
    """Adjusted Fisher-Pearson skewness."""
    r = np.asarray(returns, float)
    n = r.size
    if n < 3:
        return 0.0
    d = r - r.mean()
    m2 = float((d**2).mean())
    if m2 < 1e-30:
        return 0.0
    skew_raw = float((d**3).mean()) / m2**1.5
    return skew_raw * math.sqrt(n * (n - 1.0)) / (n - 2.0)


def excess_kurtosis(returns: np.ndarray) -> float:
    """Bias-corrected excess kurtosis."""
    r = np.asarray(returns, float)
    n = r.size
    if n < 4:
        return 0.0
    d = r - r.mean()
    m2 = float((d**2).mean())
    if m2 < 1e-30:
        return 0.0
    kurt_raw = float((d**4).mean()) / (m2 * m2)
    denom = (n - 2.0) * (n - 3.0)
    if abs(denom) < 1e-30:
        return 0.0
    return (n - 1.0) * ((n + 1.0) * (kurt_raw - 3.0) + 6.0) / denom


@dataclass(frozen=True, slots=True)
class DsrResult:
    observed_sharpe: float
    expected_max_sharpe: float
    sharpe_std_error: float
    deflated_sharpe: float          # the z-statistic
    pvalue: float                   # 1 - Φ(z); DSR "probability" = Φ(z) = 1 - pvalue


def deflated_sharpe_ratio(observed_sharpe: float, num_trials: int, sharpe_variance: float,
                          skew: float, kurt: float, num_observations: int) -> DsrResult:
    """Deflated Sharpe (Bailey & LdP 2014) with the skew/kurtosis-adjusted standard error."""
    n = float(max(num_trials, 1))
    obs = float(max(num_observations, 2))
    sr_std = math.sqrt(abs(sharpe_variance))
    term1 = (1.0 - _EULER) * normal_ppf(1.0 - 1.0 / n) if n > 1.0 else 0.0
    n_e = n * math.e
    term2 = _EULER * normal_ppf(1.0 - 1.0 / n_e) if n_e > 1.0 else 0.0
    expected_max = sr_std * (term1 + term2)

    se_inner = max(1.0 - skew * observed_sharpe + kurt / 4.0 * observed_sharpe**2, 1e-10)
    se = math.sqrt(se_inner / (obs - 1.0))
    z = (observed_sharpe - expected_max) / se if se > 1e-15 else 0.0
    return DsrResult(observed_sharpe, expected_max, se, z, 1.0 - normal_cdf(z))


@dataclass(frozen=True, slots=True)
class PermutationResult:
    observed_sharpe: float
    pvalue: float
    num_permutations: int
    null_mean: float
    null_std: float


def permutation_test(daily_returns: np.ndarray, num_permutations: int, seed: int) -> PermutationResult:
    """Monte-Carlo permutation test: shuffle returns, fraction of null Sharpes ≥ observed."""
    r = np.asarray(daily_returns, float)
    observed = sharpe_ratio(r)
    if r.size < 2 or num_permutations == 0:
        return PermutationResult(observed, 1.0, 0, 0.0, 0.0)
    rng = np.random.default_rng(seed)
    nulls = np.empty(num_permutations)
    work = r.copy()
    for i in range(num_permutations):
        rng.shuffle(work)
        nulls[i] = sharpe_ratio(work)
    pvalue = float(np.count_nonzero(nulls >= observed)) / num_permutations
    null_std = float(nulls.std(ddof=1)) if num_permutations > 1 else 0.0
    return PermutationResult(observed, pvalue, num_permutations, float(nulls.mean()), null_std)


def binomial_coefficient(n: int, k: int) -> int:
    return math.comb(n, k) if 0 <= k <= n else 0


@dataclass(frozen=True, slots=True)
class PboResult:
    pbo: float
    num_combinations: int
    logit_distribution: list[float]


def probability_of_backtest_overfitting(per_fold_sharpes: list[float]) -> PboResult:
    """PBO via CSCV (Bailey et al. 2017): over all C(N, N/2) IS/OOS splits, the fraction where the
    in-sample mean Sharpe beats the out-of-sample mean. ~0.5 = no overfit, >>0.5 = temporal overfit."""
    n = len(per_fold_sharpes)
    if n < 4 or n % 2 != 0:
        return PboResult(0.5, 0, [])
    half = n // 2
    total = binomial_coefficient(n, half)
    if total <= 10_000:
        combos = list(itertools.combinations(range(n), half))
    else:
        rng = np.random.default_rng(42)
        seen, combos = set(), []
        while len(combos) < 10_000:
            combo = tuple(sorted(rng.choice(n, size=half, replace=False).tolist()))
            if combo not in seen:
                seen.add(combo)
                combos.append(combo)
    s = per_fold_sharpes
    logits, overfit = [], 0
    for combo in combos:
        is_set = set(combo)
        is_mean = sum(s[i] for i in combo) / half
        oos_mean = sum(s[i] for i in range(n) if i not in is_set) / (n - half)
        degradation = is_mean - oos_mean
        denom = abs(is_mean) + abs(oos_mean)
        ratio = min(max(abs(oos_mean) / denom, 0.01), 0.99) if denom > 1e-15 else 0.5
        signed = min(ratio, 0.49) if degradation > 0.0 else max(ratio, 0.51)
        logits.append(math.log(signed / (1.0 - signed)))
        if is_mean > oos_mean:
            overfit += 1
    pbo = overfit / len(logits) if logits else 0.5
    return PboResult(pbo, len(combos), logits)
