"""Métricas de §6: Sharpe, DSR, PBO (CSCV) y las métricas descriptivas.

Todas se calculan con código propio a partir de métodos publicados (Bailey y López de Prado,
2014; Bailey, Borwein, López de Prado y Zhu, 2017). Solo usa la biblioteca estándar y numpy.
"""

from __future__ import annotations

import math
from collections.abc import Sequence
from dataclasses import dataclass
from datetime import date
from itertools import combinations
from statistics import NormalDist

import numpy as np

from estudio01.engine import Run

TRADING_DAYS = 252  # protocolo §6
EULER_GAMMA = 0.5772156649  # protocolo §6: constante de Euler-Mascheroni
N_BLOCKS = 16  # protocolo §6
DAYS_PER_YEAR = 365.25  # convencion del implementador (ver run_metrics)


def total_return(returns: Sequence[float]) -> float:
    """Retorno total: producto de (1 + R_t) menos 1, multiplicado en orden."""
    return math.prod(1.0 + value for value in returns) - 1.0


def is_constant(values: Sequence[float]) -> bool:
    """True si la desviación estándar es exactamente 0: todos los valores son idénticos.

    Se detecta así y no comparando la desviación calculada con 0, porque el ruido de punto flotante
    da una desviación diminuta (pero no nula) a una serie constante distinta de cero.
    """
    array = np.asarray(values, dtype=np.float64)
    return array.size < 2 or bool(array.max() == array.min())


def sharpe_per_bar(returns: Sequence[float]) -> float:
    """Media entre desviación estándar muestral (divisor n - 1); 0 si la desviación es 0.

    La desviación es 0 exactamente cuando todos los valores son idénticos (ver `is_constant`).
    """
    values = np.asarray(returns, dtype=np.float64)
    if is_constant(values):
        return 0.0
    return float(np.mean(values)) / float(np.std(values, ddof=1))


def sharpe_annualized(returns: Sequence[float]) -> float:
    """Sharpe por barra multiplicado por la raíz de 252 (tasa libre de riesgo 0)."""
    return sharpe_per_bar(returns) * math.sqrt(TRADING_DAYS)


def skew_kurtosis(returns: Sequence[float]) -> tuple[float, float]:
    """Asimetría m3 / m2^(3/2) y curtosis de Pearson m4 / m2^2, con momentos sin corrección."""
    values = np.asarray(returns, dtype=np.float64)
    if is_constant(values):
        raise ValueError("los retornos no tienen varianza: asimetría y curtosis no están definidas")
    centered = values - np.mean(values)
    m2 = float(np.mean(centered**2))
    m3 = float(np.mean(centered**3))
    m4 = float(np.mean(centered**4))
    return m3 / m2**1.5, m4 / m2**2


def sample_variance(values: Sequence[float]) -> float:
    """Varianza muestral (divisor N - 1) de los Sharpe por barra de las N configuraciones."""
    return float(np.var(np.asarray(values, dtype=np.float64), ddof=1))


def expected_max_sharpe(variance: float, n_trials: int) -> float:
    """SR0 = sqrt(V) * [(1 - g) * inv(1 - 1/N) + g * inv(1 - 1/(N * e))], g = Euler-Mascheroni."""
    normal = NormalDist()
    return math.sqrt(variance) * (
        (1.0 - EULER_GAMMA) * normal.inv_cdf(1.0 - 1.0 / n_trials)
        + EULER_GAMMA * normal.inv_cdf(1.0 - 1.0 / (n_trials * math.e))
    )


def probabilistic_sharpe(
    sr: float, reference: float, n_returns: int, skew: float, kurt: float
) -> float:
    """Phi((SR - ref) * sqrt(T - 1) / sqrt(1 - skew * SR + (kurt - 1) / 4 * SR^2))."""
    radicand = 1.0 - skew * sr + (kurt - 1.0) / 4.0 * sr * sr
    if radicand <= 0.0:
        raise ValueError(f"el radicando del DSR no es positivo ({radicand})")
    z = (sr - reference) * math.sqrt(n_returns - 1) / math.sqrt(radicand)
    return NormalDist().cdf(z)


@dataclass(frozen=True)
class Dsr:
    sr_hat: float
    skew: float
    kurt: float
    variance: float
    sr0: float
    dsr: float
    psr0: float
    n_returns: int
    n_trials: int
    passes: bool  # DSR >= 0.95 evaluado sobre el valor de 64 bits, sin redondear


DSR_THRESHOLD = 0.95  # protocolo §7


def deflated_sharpe(
    sr_hat: float,
    n_returns: int,
    skew: float,
    kurt: float,
    variance: float,
    n_trials: int,
) -> Dsr:
    """DSR de la variante 1 con su SR0, el PSR contra cero y el veredicto de §7."""
    sr0 = expected_max_sharpe(variance, n_trials)
    dsr = probabilistic_sharpe(sr_hat, sr0, n_returns, skew, kurt)
    psr0 = probabilistic_sharpe(sr_hat, 0.0, n_returns, skew, kurt)
    return Dsr(
        sr_hat, skew, kurt, variance, sr0, dsr, psr0, n_returns, n_trials, dsr >= DSR_THRESHOLD
    )


def average_rank(values: Sequence[float], index: int) -> float:
    """Rango de values[index] con 1 = menor; los empates reciben el rango promedio."""
    target = values[index]
    below = sum(1 for value in values if value < target)
    equal = sum(1 for value in values if value == target)
    return below + (equal + 1) / 2.0


@dataclass(frozen=True)
class Pbo:
    pbo: float
    n_combinations: int
    n_returns: int  # T
    n_used: int  # retornos que quedan tras descartar los T mod 16 más antiguos
    n_dropped: int
    block_length: int
    omega_counts: tuple[tuple[float, int], ...]  # (omega, número de combinaciones), ascendente


def sharpe_columns(selected: np.ndarray) -> list[float]:
    """Sharpe por barra de cada columna (divisor n - 1); 0 si la columna es constante."""
    mean = selected.mean(axis=0)
    std = selected.std(axis=0, ddof=1)
    constant = selected.max(axis=0) == selected.min(axis=0)  # desviación exactamente 0
    return [
        0.0 if flat else float(m) / float(s) for flat, m, s in zip(constant, mean, std, strict=True)
    ]


def combination_omega(blocks: np.ndarray, chosen: Sequence[int]) -> float:
    """omega = rango fuera de muestra / (N + 1) de la configuración ganadora dentro de muestra.

    `blocks` tiene forma (bloques, retornos por bloque, N); `chosen` son los bloques dentro de
    muestra y el resto es fuera de muestra. Gana dentro de muestra la configuración de mayor
    Sharpe; un empate lo gana la de menor número de variante (la primera columna).
    """
    n_blocks, _, n_configs = blocks.shape
    rest = sorted(set(range(n_blocks)).difference(chosen))
    in_sample = sharpe_columns(blocks[list(chosen)].reshape(-1, n_configs))
    out_sample = sharpe_columns(blocks[rest].reshape(-1, n_configs))
    winner = int(np.argmax(in_sample))  # argmax devuelve la primera de las empatadas
    return average_rank(out_sample, winner) / (n_configs + 1)


def pbo_cscv(matrix: np.ndarray, n_blocks: int = N_BLOCKS) -> Pbo:
    """PBO por validación cruzada combinatoria simétrica sobre una matriz T x N de retornos.

    Se descartan los (T mod 16) retornos más antiguos y el resto se divide en 16 bloques
    contiguos iguales. En cada una de las C(16, 8) combinaciones, 8 bloques son dentro de muestra
    y 8 fuera (ver `combination_omega`). lambda = ln(omega / (1 - omega)) y PBO = proporción de
    combinaciones con lambda <= 0.
    """
    returns = np.asarray(matrix, dtype=np.float64)
    if returns.ndim != 2:
        raise ValueError("la matriz de retornos debe ser T x N")
    n_returns, n_configs = returns.shape
    dropped = n_returns % n_blocks
    used = returns[dropped:]
    block_length = used.shape[0] // n_blocks
    if block_length < 1:
        raise ValueError("no hay retornos suficientes para 16 bloques")
    blocks = used.reshape(n_blocks, block_length, n_configs)

    counts: dict[float, int] = {}
    non_positive = 0
    total = 0
    for chosen in combinations(range(n_blocks), n_blocks // 2):
        omega = combination_omega(blocks, chosen)
        lam = math.log(omega / (1.0 - omega))
        non_positive += 1 if lam <= 0.0 else 0
        counts[omega] = counts.get(omega, 0) + 1
        total += 1
    return Pbo(
        pbo=non_positive / total,
        n_combinations=total,
        n_returns=n_returns,
        n_used=used.shape[0],
        n_dropped=dropped,
        block_length=block_length,
        omega_counts=tuple(sorted(counts.items())),
    )


def equity_curve(returns: Sequence[float]) -> list[float]:
    """Curva de capital que arranca en 1 y compone los retornos."""
    curve = [1.0]
    for value in returns:
        curve.append(curve[-1] * (1.0 + value))
    return curve


def max_drawdown(curve: Sequence[float]) -> float:
    """Mayor caída desde un máximo previo, como fracción positiva del máximo (0 si no hay)."""
    peak, worst = curve[0], 0.0
    for value in curve:
        peak = max(peak, value)
        worst = max(worst, (peak - value) / peak)
    return worst


@dataclass(frozen=True)
class RunMetrics:
    n_returns: int
    total_return: float
    sharpe_per_bar: float
    sharpe_annualized: float
    max_drawdown: float
    executions: int
    fraction_in_position: float  # entre 0 y 1; run_study la escribe como porcentaje (0 a 100)
    fixes_per_year: float
    zero_deviation: bool  # la desviación estándar fue 0 y el Sharpe se reporta como 0 (§6)


def run_metrics(run: Run, start: date | None = None, end: date | None = None) -> RunMetrics:
    """Métricas de una corrida, en toda la ventana o solo en las barras con start <= fecha <= end.

    Con start y end es el tramo final: la posición se hereda de la corrida completa, el retorno
    total es el producto de (1 + R_t) menos 1, el Sharpe usa raíz de 252 como en la ventana y el
    drawdown máximo se calcula sobre la curva reiniciada en 1 al inicio del tramo. El drawdown de
    toda la ventana incluye también el capital inicial 1. Una ejecución es una barra con al menos
    un lado, incluidos la entrada inicial y el cierre. «Fijos por año» es T dividido entre los
    años (de 365.25 días) que van de la primera a la última barra considerada.
    """
    picks = [
        i
        for i, day in enumerate(run.dates)
        if (start is None or day >= start) and (end is None or day <= end)
    ]
    if not picks:
        return RunMetrics(0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 0.0, True)
    returns = [run.ret[i] for i in picks]
    span_days = (run.dates[picks[-1]] - run.dates[picks[0]]).days
    return RunMetrics(
        n_returns=len(picks),
        total_return=total_return(returns),
        sharpe_per_bar=sharpe_per_bar(returns),
        sharpe_annualized=sharpe_annualized(returns),
        max_drawdown=max_drawdown(equity_curve(returns)),
        executions=sum(1 for i in picks if run.sides[i] > 0),
        fraction_in_position=sum(1 for i in picks if run.d[i] != 0) / len(picks),
        fixes_per_year=0.0 if span_days == 0 else len(picks) / (span_days / DAYS_PER_YEAR),
        zero_deviation=is_constant(returns),
    )