from __future__ import print_function
import numpy as np
import itertools as itr
import scipy.stats as ss
import scipy.special as spec
import statsmodels.distributions.empirical_distribution as smd
import matplotlib.pyplot as plt
import collections as cls
import pandas as pd
import joblib as job
import psutil as ps
import util.metrics as perf


PBO = cls.namedtuple(
    "PBO",
    [
        "pbo",
        "prob_oos_loss",
        "linear_model",
        "stochastic",
        "Cs",
        "J",
        "J_bar",
        "R",
        "R_bar",
        "R_rank",
        "R_bar_rank",
        "rn",
        "rn_bar",
        "w_bar",
        "logits",
        "R_n_star",
        "R_bar_n_star",
    ],
)


PBOCore = cls.namedtuple(
    "PBOCore",
    [
        "J",
        "J_bar",
        "R",
        "R_bar",
        "R_rank",
        "R_bar_rank",
        "rn",
        "rn_bar",
        "w_bar",
        "logits",
    ],
)


def pbo(M, S, metric_func, threshold=0.0, n_jobs=1, verbose=False):
    """
    Based on http://papers.ssrn.com/sol3/papers.cfm?abstract_id=2326253

    Features:
    * training and test sets are of equal size, providing comparable accuracy
    to both IS and OOS Sharpe ratios.
    * CSCV is symmetric, decline in performance can only result from
    overfitting, not arbitrary discrepancies between the training and test
    sets.
    * CSCV respects the time-dependence and other season-dependent features
    present in the data.
    * Results are deterministic, can be replicated.
    * Dispersion in the distribution of logits conveys relevant info regarding
    the robustness of the strategy selection process.
    * Model-free, non-parametric. Logits distribution resembles the cumulative
    Normal distribution if w_bar are close to uniform distribution (i.e. the
    backtest appears to be information-less). Therefore, for good backtesting,
    the distribution of logits will be centered in a significantly positive
    value, and its tail will marginally cover the region of negative logit
    values.

    Limitations:
    * CSCV is symmetric, for some strategies, K-fold CV might be better.
    * Not suitable for time series with strong auto-correlation, especially
    when S is large.
    * Assumes all the sample statistics carry the same weight.
    * Entirely possible that all the N strategy configs have high but similar
    Sharpe ratios. Therefore, PBO may appear high, however, 'overfitting' here
    is among many 'skilful' strategies.

    Parameters:

    M:
        returns data, numpy or dataframe format.
    S:
        chuncks to devided M into, must be even number. Paper suggests setting
        S = 16. See paper for details of choice of S.
    metric_func:
        evaluation function for returns data
    threshold:
        used as prob. of OOS Loss calculation cutoff. For Sharpe ratio,
        this should be 0 to indicate probabilty of loss.
    n_jobs:
        if greater than 1 then enable parallel mode
    hist:
        Default False, whether to plot histogram for rank of logits.
        Some problems exist when S >= 10. Need to look at why numpy /
        matplotlib does it.

    Returns:
    PBO result in namedtuple, instance of PBO.
    """
    if S % 2 == 1:
        raise ValueError(
            "S must be an even integer, {:.1f} was given".format(S)
        )

    n_jobs = int(n_jobs)
    if n_jobs < 0:
        n_jobs = max(1, ps.cpu_count(logical=False))

    if isinstance(M, pd.DataFrame):
        # conver to numpy values
        if verbose:
            print("Convert from DataFrame to numpy array.")
        M = M.values

    # Paper suggests T should be 2x the no. of observations used by investor
    # to choose a model config, due to the fact that CSCV compares combinations
    # of T/2 observations with their complements.
    T, N = M.shape
    residual = T % S
    if residual != 0:
        M = M[residual:]
        T, N = M.shape

    sub_T = T // S

    if verbose:
        print("Total sample size: {:,d}, chunck size: {:,d}".format(T, sub_T))

    # generate subsets, each of length sub_T
    Ms = []
    Ms_values = []
    for i in range(S):
        start, end = i * sub_T, (i + 1) * sub_T
        Ms.append((i, M[start:end, :]))
        Ms_values.append(M[start:end, :])
    Ms_values = np.array(Ms_values)

    if verbose:
        print("No. of Chuncks: {:,d}".format(len(Ms)))

    # generate combinations
    Cs = [x for x in itr.combinations(Ms, S // 2)]
    if verbose:
        print("No. of combinations = {:,d}".format(len(Cs)))

    # Ms_index used to find J_bar (complementary OOS part)
    Ms_index = set([x for x in range(len(Ms))])

    # create J and J_bar
    if n_jobs < 2:
        J = []
        J_bar = []

        for i in range(len(Cs)):
            # make sure chucks are concatenated in their original order
            order = [x for x, _ in Cs[i]]
            sort_ind = np.argsort(order)

            Cs_values = np.array([v for _, v in Cs[i]])
            # if verbose:
            #     print('Cs index = {}, '.format(order), end='')
            joined = np.concatenate(Cs_values[sort_ind, :])
            J.append(joined)

            # find Cs_bar
            Cs_bar_index = list(sorted(Ms_index - set(order)))
            # if verbose:
            # print('Cs_bar_index = {}'.format(Cs_bar_index))
            J_bar.append(np.concatenate(Ms_values[Cs_bar_index, :]))

        # compute matrices for J and J_bar, e.g. Sharpe ratio
        R = [metric_func(j) for j in J]
        R_bar = [metric_func(j) for j in J_bar]

        # compute ranks of metrics
        R_rank = [ss.rankdata(x) for x in R]
        R_bar_rank = [ss.rankdata(x) for x in R_bar]

        # find highest metric, rn contains the index position of max value
        # in each set of R (IS)
        rn = [np.argmax(r) for r in R_rank]
        # use above index to find R_bar (OOS) in same index position
        # i.e. the same config / setting
        rn_bar = [R_bar_rank[i][rn[i]] for i in range(len(R_bar_rank))]

        # formula in paper used N+1 as the denominator for w_bar.
        w_bar = [float(r) / N for r in rn_bar]
        # logit(.5) gives 0 so if w_bar value is equal to median logits is 0
        logits = [spec.logit(w) for w in w_bar]
    else:
        # use joblib for parallel calc
        # print('Run in parallel mode.')
        cores = job.Parallel(n_jobs=n_jobs)(
            job.delayed(pbo_core_calc)(
                Cs_x, Ms, Ms_values, Ms_index, metric_func, verbose
            )
            for Cs_x in Cs
        )

        J = [c.J for c in cores]
        J_bar = [c.J_bar for c in cores]
        R = [c.R for c in cores]
        R_bar = [c.R_bar for c in cores]
        R_rank = [c.R_rank for c in cores]
        R_bar_rank = [c.R_bar_rank for c in cores]
        rn = [c.rn for c in cores]
        rn_bar = [c.rn_bar for c in cores]
        w_bar = [c.w_bar for c in cores]
        logits = [c.logits for c in cores]

    # prob of overfitting
    phi = np.array([1.0 if lam <= 0 else 0.0 for lam in logits]) / len(Cs)
    pbo_test = np.sum(phi)

    # performance degradation
    R_n_star = np.array([R[i][rn[i]] for i in range(len(R))])
    R_bar_n_star = np.array([R_bar[i][rn[i]] for i in range(len(R_bar))])
    lm = ss.linregress(x=R_n_star, y=R_bar_n_star)

    prob_oos_loss = np.sum(
        [1.0 if r < threshold else 0.0 for r in R_bar_n_star]
    ) / len(R_bar_n_star)

    # Stochastic dominance
    y = np.linspace(
        min(R_bar_n_star), max(R_bar_n_star), endpoint=True, num=1000
    )
    R_bar_n_star_cdf = smd.ECDF(R_bar_n_star)
    optimized = R_bar_n_star_cdf(y)

    R_bar_cdf = smd.ECDF(np.concatenate(R_bar))
    non_optimized = R_bar_cdf(y)

    dom_df = pd.DataFrame(
        dict(optimized_IS=optimized, non_optimized_OOS=non_optimized)
    )
    dom_df.index = y
    # visually, non_optimized curve above optimized curve indicates good
    # backtest with low overfitting.
    dom_df["SD2"] = dom_df.non_optimized_OOS - dom_df.optimized_IS

    result = PBO(
        pbo_test,
        prob_oos_loss,
        lm,
        dom_df,
        Cs,
        J,
        J_bar,
        R,
        R_bar,
        R_rank,
        R_bar_rank,
        rn,
        rn_bar,
        w_bar,
        logits,
        R_n_star,
        R_bar_n_star,
    )
    return result


def pbo_core_calc(Cs, Ms, Ms_values, Ms_index, metric_func, verbose=False):
    # make sure chucks are concatenated in their original order
    order = [x for x, _ in Cs]
    sort_ind = np.argsort(order)

    Cs_values = np.array([v for _, v in Cs])
    if verbose:
        print("Cs index = {}, ".format(order), end="")
    J_x = np.concatenate(Cs_values[sort_ind, :])

    # find Cs_bar
    Cs_bar_index = list(sorted(Ms_index - set(order)))
    if verbose:
        print("Cs_bar_index = {}".format(Cs_bar_index))
    J_bar_x = np.concatenate(Ms_values[Cs_bar_index, :])

    R_x = metric_func(J_x)
    R_bar_x = metric_func(J_bar_x)

    R_rank_x = ss.rankdata(R_x)
    R_bar_rank_x = ss.rankdata(R_bar_x)

    rn_x = np.argmax(R_rank_x)
    rn_bar_x = R_bar_rank_x[rn_x]

    w_bar_x = float(rn_bar_x) / len(R_bar_rank_x)
    logit_x = spec.logit(w_bar_x)

    core = PBOCore(
        J_x,
        J_bar_x,
        R_x,
        R_bar_x,
        R_rank_x,
        R_bar_rank_x,
        rn_x,
        rn_bar_x,
        w_bar_x,
        logit_x,
    )

    return core

def psr_from_returns(returns, risk_free=0, target_sharpe=0):
    """
    PSR from return series.

    Parameters:
        returns:
            return series
        risk_free:
            risk free or benchmark rate for sharpe ratio calculation,
            default 0.
        target_sharpe:
            minimum sharpe ratio

    Returns:
        PSR probabilities.
    """
    T = len(returns)
    sharpe = perf.sharpe_iid(returns, bench=risk_free, factor=1)
    skew = returns.skew()
    kurtosis = returns.kurtosis() + 3

    return psr(
        sharpe=sharpe,
        T=T,
        skew=skew,
        kurtosis=kurtosis,
        target_sharpe=target_sharpe,
    )


def psr(sharpe, T, skew, kurtosis, target_sharpe=0):
    """
    Probabilistic Sharpe Ratio.

    Parameters:
        sharpe:
            observed sharpe ratio, in same frequency as T.
        T:
            no. of observations, should match return / sharpe sampling period.
        skew:
            sharpe ratio skew
        kurtosis:
            sharpe ratio kurtosis
        target_sharpe:
            target sharpe ratio

    Returns:
        Cumulative probabilities for observed sharpe ratios under standard
        Normal distribution.
    """
    value = (
        (sharpe - target_sharpe)
        * np.sqrt(T - 1)
        / np.sqrt(1.0 - skew * sharpe + sharpe ** 2 * (kurtosis - 1) / 4.0)
    )
    # print(value)
    psr = ss.norm.cdf(value, 0, 1)
    return psr


def dsr(test_sharpe, sharpe_std, N, T, skew, kurtosis):
    """
    Deflated Sharpe Ratio statistic. DSR = PSR(SR_0).
    See paper for definition of SR_0. http://ssrn.com/abstract=2460551

    Parameters:
        test_sharpe :
            reported sharpe, to be tested.
        sharpe_std :
            standard deviation of sharpe ratios from N trials / configurations
        N :
            number of backtest configurations
        T :
            number of observations
        skew :
            skew of returns
        kurtosis :
            kurtosis of returns

    Returns:
        DSR statistic
    """
    # sharpe_std = np.std(sharpe_n, ddof=1)
    target_sharpe = sharpe_std * expected_max(N)

    dsr_stat = psr(test_sharpe, T, skew, kurtosis, target_sharpe)

    return dsr_stat


def dsr_from_returns(test_sharpe, returns_df, risk_free=0):
    """
    Calculate DSR based on a set of given returns_df.

    Parameters:
        test_sharpe :
            Reported sharpe, to be tested.
        returns_df :
            Log return series
        risk_free :
            Risk free return, default 0.
    Returns:
        DSR statistic
    """
    T, N = returns_df.shape
    sharpe = perf.sharpe_iid(returns_df, bench=risk_free, factor=1)
    sharpe_std = np.std(sharpe, ddof=1)
    skew = returns_df.skew()
    kurtosis = returns_df.kurtosis() + 3

    out = dsr(
        test_sharpe,
        sharpe_std=sharpe_std,
        N=N,
        T=T,
        skew=skew,
        kurtosis=kurtosis,
    )

    return out


def minTRL(sharpe, skew, kurtosis, target_sharpe=0, prob=0.95):
    """
    Minimum Track Record Length.

    Parameters:
        sharpe :
            observed sharpe ratio, in same frequency as observations.
        skew :
            sharpe ratio skew
        kurtosis :
            sharpe ratio kurtosis
        target_sharpe :
            target sharpe ratio
        prob :
            minimum probability for estimating TRL.

    Returns:
        minTRL, in terms of number of observations.
    """
    min_track = (
        1
        + (1 - skew * sharpe + sharpe ** 2 * (kurtosis - 1) / 4.0)
        * (ss.norm.ppf(prob) / (sharpe - target_sharpe)) ** 2
    )
    return min_track


def expected_max(N):
    """
    Expected maximum of IID random variance X_n ~ Z, n = 1,...,N,
    where Z is the CDF of the standard Normal distribution,
    E[MAX_n] = E[max{x_n}]. Computed for a large N.

    """
    if N < 5:
        raise AssertionError("Condition N >> 1 not satisfied.")
    return (1 - np.euler_gamma) * ss.norm.ppf(
        1 - 1.0 / N
    ) + np.euler_gamma * ss.norm.ppf(1 - np.exp(-1) / N)


def minBTL(N, sharpe_IS):
    """
    Minimum backtest length. minBTL should be considered a necessary,
    non-sufficient condition to avoid overfitting. See PBO for a more precise
    measure of backtest overfitting.

    Paramters:
        N :
            number of backtest configurations
        sharpe_IS :
            In-Sample observed Sharpe ratio

    Returns:
        btl :
            minimum back test length
        upper_bound :
            upper bound for minBTL
    """
    exp_max = expected_max(N)

    btl = (exp_max / sharpe_IS) ** 2

    upper_bound = 2 * np.log(N) / sharpe_IS ** 2

    return (btl, upper_bound)
