"""From-scratch deflated Sharpe ratio (Bailey & Lopez de Prado 2014) and the
Probability of Backtest Overfitting via Combinatorial Symmetric CV (CSCV)."""
import numpy as np
from itertools import combinations
from scipy.stats import norm, skew, kurtosis

def prob_sharpe_ratio(sr, T, g3, g4, sr_star=0.0):
    """PSR: probability the true (per-period) Sharpe exceeds sr_star, adjusting
    for skew g3 and kurtosis g4 (non-normal returns)."""
    denom = np.sqrt(1 - g3*sr + (g4-1)/4*sr**2)
    return float(norm.cdf((sr - sr_star)*np.sqrt(T-1)/denom))

def expected_max_sharpe(sr_std, N):
    """Expected maximum of N independent Sharpe estimates under the null (skill=0)."""
    g = 0.5772156649015329                                  # Euler-Mascheroni
    return sr_std*((1-g)*norm.ppf(1-1.0/N) + g*norm.ppf(1-1.0/(N*np.e)))

def deflated_sharpe(returns_best, all_sr, T=None):
    """Deflated Sharpe = PSR of the selected strategy against the expected-max
    Sharpe implied by having tried len(all_sr) strategies."""
    x = np.asarray(returns_best); T = T or len(x)
    sr = x.mean()/x.std()
    sr0 = expected_max_sharpe(np.std(all_sr, ddof=1), len(all_sr))
    return prob_sharpe_ratio(sr, T, skew(x), kurtosis(x, fisher=False), sr0), sr, sr0

def pbo_cscv(M, S=8):
    """Probability of Backtest Overfitting via CSCV. M is (T, N) of per-period
    strategy returns. Returns (pbo, logits)."""
    T, N = M.shape
    blocks = np.array_split(np.arange(T), S)
    logits = []
    for tr in combinations(range(S), S//2):
        te = [b for b in range(S) if b not in tr]
        i_tr = np.concatenate([blocks[b] for b in tr]); i_te = np.concatenate([blocks[b] for b in te])
        sr_is = M[i_tr].mean(0)/(M[i_tr].std(0)+1e-12)
        n_star = int(np.argmax(sr_is))                     # best strategy in-sample
        sr_oos = M[i_te].mean(0)/(M[i_te].std(0)+1e-12)
        rank = (sr_oos.argsort().argsort()[n_star] + 1)/(N + 1)   # OOS relative rank of the IS-best
        rank = min(max(rank, 1e-6), 1-1e-6)
        logits.append(np.log(rank/(1-rank)))
    logits = np.array(logits)
    return float(np.mean(logits <= 0)), logits              # PBO = P(IS-best underperforms OOS median)