import pandas as pd
import numpy as np


def add_returns(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy()
    df['return_1'] = df['Close'].pct_change(1)
    df['return_5'] = df['Close'].pct_change(5)
    df['return_20'] = df['Close'].pct_change(20)
    return df


def add_sma(df: pd.DataFrame, windows: list[int] | None = None) -> pd.DataFrame:
    if windows is None:
        windows = [5, 10, 20, 50]
    df = df.copy()
    for w in windows:
        df[f'sma_{w}'] = df['Close'].rolling(window=w).mean()
        df[f'close_to_sma_{w}'] = df['Close'] / df[f'sma_{w}'] - 1
    return df


def add_ema(df: pd.DataFrame, windows: list[int] | None = None) -> pd.DataFrame:
    if windows is None:
        windows = [12, 26]
    df = df.copy()
    for w in windows:
        df[f'ema_{w}'] = df['Close'].ewm(span=w, adjust=False).mean()
    if 12 in windows and 26 in windows:
        df['macd'] = df['ema_12'] - df['ema_26']
        df['macd_signal'] = df['macd'].ewm(span=9, adjust=False).mean()
        df['macd_hist'] = df['macd'] - df['macd_signal']
    return df


def add_rsi(df: pd.DataFrame, period: int = 14) -> pd.DataFrame:
    df = df.copy()
    delta = df['Close'].diff()
    gain = delta.where(delta > 0, 0.0)
    loss = (-delta).where(delta < 0, 0.0)
    avg_gain = gain.rolling(window=period, min_periods=1).mean()
    avg_loss = loss.rolling(window=period, min_periods=1).mean()
    rs = avg_gain / avg_loss.replace(0, np.nan)
    df['rsi'] = 100 - (100 / (1 + rs))
    df['rsi'] = df['rsi'].fillna(50)
    return df


def add_bollinger(df: pd.DataFrame, period: int = 20, num_std: float = 2.0) -> pd.DataFrame:
    df = df.copy()
    df['bb_mid'] = df['Close'].rolling(window=period).mean()
    bb_std = df['Close'].rolling(window=period).std()
    df['bb_upper'] = df['bb_mid'] + num_std * bb_std
    df['bb_lower'] = df['bb_mid'] - num_std * bb_std
    df['bb_width'] = (df['bb_upper'] - df['bb_lower']) / df['bb_mid']
    df['bb_position'] = (df['Close'] - df['bb_lower']) / (df['bb_upper'] - df['bb_lower'] + 1e-10)
    return df


def add_atr(df: pd.DataFrame, period: int = 14) -> pd.DataFrame:
    df = df.copy()
    high_low = df['High'] - df['Low']
    high_close = (df['High'] - df['Close'].shift()).abs()
    low_close = (df['Low'] - df['Close'].shift()).abs()
    tr = pd.concat([high_low, high_close, low_close], axis=1).max(axis=1)
    df['atr'] = tr.rolling(window=period, min_periods=1).mean()
    df['atr_pct'] = df['atr'] / df['Close']
    return df


def add_volume_features(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy()
    df['volume_sma_5'] = df['Volume'].rolling(window=5).mean()
    df['volume_ratio'] = df['Volume'] / df['volume_sma_5'].replace(0, np.nan)
    df['volume_ratio'] = df['volume_ratio'].fillna(1)
    return df


def add_interaction_features(df: pd.DataFrame) -> pd.DataFrame:
    df = df.copy()
    req = ['bb_width', 'atr_pct', 'rsi', 'bb_position', 'close_to_sma_50', 'return_20', 'volume_ratio']
    if all(c in df.columns for c in req):
        df['bb_width_x_atr_pct'] = df['bb_width'] * df['atr_pct']
        df['rsi_x_bb_pos'] = (df['rsi'] - 50) * (df['bb_position'] - 0.5)
        df['cts50_x_ret20'] = df['close_to_sma_50'] * df['return_20']
        df['bb_width_x_rsi'] = df['bb_width'] * (df['rsi'] - 50)
        df['atr_pct_x_vol_ratio'] = df['atr_pct'] * df['volume_ratio']
    return df


INTERACTION_COLS = [
    'bb_width_x_atr_pct',
    'rsi_x_bb_pos',
    'cts50_x_ret20',
    'bb_width_x_rsi',
    'atr_pct_x_vol_ratio',
]


def add_target(df: pd.DataFrame, horizon: int = 1) -> pd.DataFrame:
    df = df.copy()
    df['target'] = df['Close'].shift(-horizon) / df['Close'] - 1
    df['target_direction'] = (df['target'] > 0).astype(int)
    return df


FEATURE_COLS = [
    'return_1', 'return_5', 'return_20',
    'close_to_sma_5', 'close_to_sma_10', 'close_to_sma_20', 'close_to_sma_50',
    'rsi',
    'bb_width', 'bb_position',
    'macd', 'macd_signal', 'macd_hist',
    'atr_pct',
    'volume_ratio',
    'kst_norm', 'efficiency_ratio', 'choppiness',
]


def add_adx(df: pd.DataFrame, period: int = 14) -> pd.DataFrame:
    """
    Average Directional Index — сила и направление тренда.
    Использует Wilders Smoothing (SMMA) по оригинальному методу Уайлдера (1978).
    Добавляет колонки: adx, pdi, mdi.
    """
    df = df.copy()
    high, low, close = df['High'].values, df['Low'].values, df['Close'].values
    n = len(df)

    # True Range
    tr = np.zeros(n)
    for i in range(1, n):
        hl = high[i] - low[i]
        hc = abs(high[i] - close[i-1])
        lc = abs(low[i] - close[i-1])
        tr[i] = max(hl, hc, lc)

    # ATR через Wilders Smoothing
    atr_ = pd.Series(tr).ewm(alpha=1.0/period, adjust=False).mean().values

    # Directional Movement
    up_move = np.zeros(n)
    down_move = np.zeros(n)
    for i in range(1, n):
        up_move[i] = max(high[i] - high[i-1], 0)
        down_move[i] = max(low[i-1] - low[i], 0)

    # +DM и -DM через Wilders Smoothing
    plus_dm = pd.Series(up_move).ewm(alpha=1.0/period, adjust=False).mean().values
    minus_dm = pd.Series(down_move).ewm(alpha=1.0/period, adjust=False).mean().values

    # +DI и -DI
    atr_safe = np.where(atr_ > 1e-10, atr_, 1e-10)
    df['pdi'] = 100 * plus_dm / atr_safe
    df['mdi'] = 100 * minus_dm / atr_safe

    # DX и ADX (через Wilders Smoothing)
    dx = 100 * np.abs(df['pdi'] - df['mdi']) / np.where(df['pdi'] + df['mdi'] > 1e-10, df['pdi'] + df['mdi'], 1e-10)
    df['adx'] = pd.Series(dx).ewm(alpha=1.0/period, adjust=False).mean().values
    return df


def add_realized_volatility(df: pd.DataFrame, periods: list[int] = [5, 10, 20]) -> pd.DataFrame:
    """Реализованная волатильность (std of returns) за несколько окон."""
    df = df.copy()
    returns = df['Close'].pct_change().values
    for p in periods:
        df[f'rv_{p}'] = pd.Series(returns).rolling(p).std().values
    return df


def add_parkinson_volatility(df: pd.DataFrame, period: int = 20) -> pd.DataFrame:
    """Parkinson volatility estimator (HL-based).
    Записывается в parkinson_vol (period=20) или parkinson_vol_10 (period=10)."""
    df = df.copy()
    high, low = df['High'].values, df['Low'].values
    n = len(df)
    hl_ratio = np.zeros(n)
    for i in range(1, n):
        if high[i] > low[i] and high[i] > 0:
            hl_ratio[i] = np.log(high[i] / low[i]) ** 2
    col = 'parkinson_vol_10' if period == 10 else 'parkinson_vol'
    df[col] = np.sqrt(
        pd.Series(hl_ratio).rolling(period).sum() / (4 * np.log(2) * period)
    ).values
    return df


def add_garman_klass_volatility(df: pd.DataFrame, period: int = 20) -> pd.DataFrame:
    """Garman-Klass volatility estimator.
    Записывается в garman_klass_vol (period=20) или garman_klass_vol_10 (period=10)."""
    df = df.copy()
    high, low, open_, close = df['High'].values, df['Low'].values, df['Open'].values, df['Close'].values
    n = len(df)
    gk = np.zeros(n)
    for i in range(1, n):
        hl = np.log(high[i] / low[i]) ** 2
        co = np.log(close[i] / open_[i]) ** 2
        gk[i] = 0.5 * hl - (2 * np.log(2) - 1) * co
    col = 'garman_klass_vol_10' if period == 10 else 'garman_klass_vol'
    df[col] = np.sqrt(pd.Series(gk).rolling(period).mean()).values
    return df


def add_yang_zhang_volatility(df: pd.DataFrame, period: int = 20) -> pd.DataFrame:
    """Yang-Zhang volatility estimator (overnight + Rogers-Satchell).
    Записывается в yang_zhang_vol (period=20) или yang_zhang_vol_10 (period=10)."""
    df = df.copy()
    high, low, open_, close = df['High'].values, df['Low'].values, df['Open'].values, df['Close'].values
    n = len(df)
    rs = np.zeros(n)
    overnight = np.zeros(n)
    for i in range(1, n):
        hl = np.log(high[i] / low[i])
        co = np.log(close[i] / open_[i])
        rs[i] = hl ** 2 + co ** 2
        oc = np.log(open_[i] / close[i-1])
        overnight[i] = oc ** 2
    col = 'yang_zhang_vol_10' if period == 10 else 'yang_zhang_vol'
    df[col] = np.sqrt(
        pd.Series(rs).rolling(period).mean() + pd.Series(overnight).rolling(period).mean()
    ).values
    return df


# ─────────────────────────────────────────────────────────────────────
# v13 features (audit 2026-08-03)
# ─────────────────────────────────────────────────────────────────────

def _rs_hurst(returns: np.ndarray) -> float:
    """Hurst exponent via Rescaled Range (R/S) analysis.

    Классический метод Hurst (1951):
      Для каждого lag L: делим series на chunk'и размера L, для каждого
      считаем R = max(cumdev) - min(cumdev), где cumdev — кумулятивная
      сумма отклонений от среднего chunk'a; S = std chunk'a.
      R/S = mean(R/S). Hurst = slope log(R/S) vs log(L).

    H < 0.5 → mean-reverting, H > 0.5 → trending, H ≈ 0.5 → random walk.
    На малых выборках возвращает 0.5 (random walk fallback).
    """
    n = len(returns)
    if n < 30:
        return 0.5

    # Lag values — геометрическая прогрессия для best log-log fit
    lags = [2, 4, 8, 16, 32, 64]
    lags = [L for L in lags if L < n // 2]
    if len(lags) < 3:
        return 0.5

    rs_values = []
    for L in lags:
        # Число полных chunk'ов
        n_chunks = n // L
        if n_chunks < 1:
            continue
        rs_chunk = []
        for k in range(n_chunks):
            chunk = returns[k * L: (k + 1) * L]
            mean = chunk.mean()
            dev = chunk - mean
            cumdev = np.cumsum(dev)
            R = cumdev.max() - cumdev.min()
            S = chunk.std()
            if S > 1e-12:
                rs_chunk.append(R / S)
        if rs_chunk:
            rs_values.append((L, np.mean(rs_chunk)))

    if len(rs_values) < 3:
        return 0.5

    log_L = np.log([L for L, _ in rs_values])
    log_RS = np.log([rs for _, rs in rs_values])
    try:
        poly = np.polyfit(log_L, log_RS, 1)
        return float(np.clip(poly[0], 0.0, 1.0))
    except Exception:
        return 0.5


def add_hurst(df: pd.DataFrame, window: int = 100) -> pd.DataFrame:
    """Hurst exponent over rolling window.

    Цель (audit 2026-08-03): модель должна различать trending vs
    mean-reverting режимы. SBER в up-тренде давал 12/12 SHORT SL —
    модель не учитывала, что H>0.5 = trending, SHORT mean-revert
    стратегия обречена.
    """
    df = df.copy()
    returns = np.log(df['Close'] / df['Close'].shift(1)).fillna(0).values
    n = len(returns)

    hurst_vals = np.full(n, 0.5)
    for i in range(window, n):
        hurst_vals[i] = _rs_hurst(returns[i - window:i])

    df['hurst_100'] = hurst_vals
    # Категориальные флаги: 1 = trending (H>0.55), 1 = mean-reverting (H<0.45)
    # Нейтральная зона [0.45, 0.55] — оба флага = 0 (random walk).
    df['is_trending'] = (df['hurst_100'] > 0.55).astype(float)
    df['is_mean_reverting'] = (df['hurst_100'] < 0.45).astype(float)
    # Confidence = |H - 0.5| × 2 (0 = random, 1 = strong trend/MR)
    df['hurst_confidence'] = (df['hurst_100'] - 0.5).abs() * 2.0
    return df


HURST_COLS = ['hurst_100', 'is_trending', 'is_mean_reverting', 'hurst_confidence']


def add_vwap_features(df: pd.DataFrame, window: int = 24) -> pd.DataFrame:
    """VWAP (Volume Weighted Average Price) distance features.

    Цель: использовать объём-on-price информацию как фильтр.
    `vwap_dist_pct` > 0 — цена выше VWAP (бычий pressure).
    `above_vwap` — бинарный флаг.
    """
    df = df.copy()
    typical_price = (df['High'] + df['Low'] + df['Close']) / 3.0
    vol = df['Volume'].replace(0, np.nan)
    vp = typical_price * vol

    vol_sum = vol.rolling(window, min_periods=1).sum()
    vp_sum = vp.rolling(window, min_periods=1).sum()
    vwap = (vp_sum / vol_sum).bfill()

    df['vwap_24'] = vwap
    df['vwap_dist_pct'] = ((df['Close'] - vwap) / vwap * 100).fillna(0)
    df['above_vwap'] = (df['Close'] > vwap).astype(float)
    return df


VWAP_COLS = ['vwap_24', 'vwap_dist_pct', 'above_vwap']


def add_temporal_features(df: pd.DataFrame) -> pd.DataFrame:
    """Добавляет временные признаки: час, день недели, месяц, циклические кодировки."""
    df = df.copy()

    # Парсим дату и время из колонок Date (YYYY.MM.DD) и Time (HH:MM)
    # Если есть timestamp, используем его
    if 'timestamp' in df.columns:
        dt = pd.to_datetime(df['timestamp'], unit='s', utc=True)
        # Смещение на MSK (UTC+3) для MOEX
        try:
            dt_msk = dt.dt.tz_convert('Europe/Moscow')
        except (AttributeError, TypeError):
            # Fallback: timestamp уже в MSK или tz_convert недоступен
            dt_msk = dt
    else:
        dt = pd.to_datetime(df['Date'] + ' ' + df['Time'], format='%Y.%m.%d %H:%M')
        dt_msk = dt

    df['hour'] = dt_msk.dt.hour.astype(float)
    df['day_of_week'] = dt_msk.dt.dayofweek.astype(float)  # 0=Monday
    df['month'] = dt_msk.dt.month.astype(float)
    df['day_of_month'] = dt_msk.dt.day.astype(float)
    df['quarter'] = dt_msk.dt.quarter.astype(float)
    df['is_weekend'] = (dt_msk.dt.dayofweek >= 5).astype(float)

    # Циклическое кодирование часа (для учёта цикличности времени суток)
    hour_rad = 2 * np.pi * df['hour'] / 24
    df['hour_sin'] = np.sin(hour_rad)
    df['hour_cos'] = np.cos(hour_rad)

    # Циклическое кодирование дня недели
    dow_rad = 2 * np.pi * df['day_of_week'] / 7
    df['dow_sin'] = np.sin(dow_rad)
    df['dow_cos'] = np.cos(dow_rad)

    # Циклическое кодирование месяца
    month_rad = 2 * np.pi * df['month'] / 12
    df['month_sin'] = np.sin(month_rad)
    df['month_cos'] = np.cos(month_rad)

    # Флаги торговых сессий MOEX (MSK)
    # Основная сессия: 10:00-18:45
    # Вечерняя сессия: 19:00-23:50
    df['is_main_session'] = ((df['hour'] >= 10) & (df['hour'] < 19)).astype(float)
    df['is_evening_session'] = ((df['hour'] >= 19) & (df['hour'] < 24)).astype(float)
    df['is_market_open'] = ((df['is_main_session'] + df['is_evening_session']) >= 0.5).astype(float)

    # Время от открытия сессии (в часах)
    # Основная сессия: 10:00 → session_hour=0, 18:00 → session_hour=8
    # Вечерняя сессия: 19:00 → session_hour=0, 23:00 → session_hour=4
    df['session_hour'] = np.where(
        df['is_main_session'] == 1, df['hour'] - 10,
        np.where(df['is_evening_session'] == 1, df['hour'] - 19, -1)
    )

    return df


TEMPORAL_COLS = [
    'hour', 'day_of_week', 'month', 'day_of_month', 'quarter', 'is_weekend',
    'hour_sin', 'hour_cos', 'dow_sin', 'dow_cos', 'month_sin', 'month_cos',
    'is_main_session', 'is_evening_session', 'is_market_open', 'session_hour',
]


def compute_kst_array(df: pd.DataFrame) -> np.ndarray:
    """KST (Know Sure Thing) oscillator — взвешенная сумма сглаженных ROC.

    ROC periods: 10, 15, 20, 30
    SMA (smoothing) periods: 10, 10, 10, 15
    Weights: 1, 2, 3, 4
    KST = RCMA1 + 2*RCMA2 + 3*RCMA3 + 4*RCMA4

    Returns np.ndarray normalized to ~[-1, +1].
    """
    close = df['Close'].values
    n = len(close)

    roc10 = np.zeros(n)
    roc15 = np.zeros(n)
    roc20 = np.zeros(n)
    roc30 = np.zeros(n)
    for i in range(30, n):
        if close[i - 10] > 1e-10:
            roc10[i] = (close[i] / close[i - 10] - 1) * 100
        if close[i - 15] > 1e-10:
            roc15[i] = (close[i] / close[i - 15] - 1) * 100
        if close[i - 20] > 1e-10:
            roc20[i] = (close[i] / close[i - 20] - 1) * 100
        if close[i - 30] > 1e-10:
            roc30[i] = (close[i] / close[i - 30] - 1) * 100

    s_roc10 = pd.Series(roc10).rolling(10, min_periods=1).mean().values
    s_roc15 = pd.Series(roc15).rolling(10, min_periods=1).mean().values
    s_roc20 = pd.Series(roc20).rolling(10, min_periods=1).mean().values
    s_roc30 = pd.Series(roc30).rolling(15, min_periods=1).mean().values

    kst = s_roc10 + 2 * s_roc15 + 3 * s_roc20 + 4 * s_roc30
    return kst / 100.0


def add_kst(df: pd.DataFrame) -> pd.DataFrame:
    """KST (Know Sure Thing) oscillator — добавляет kst_norm в DataFrame."""
    df = df.copy()
    df['kst_norm'] = compute_kst_array(df)
    return df


def compute_efficiency_ratio_array(df: pd.DataFrame, period: int = 10) -> np.ndarray:
    """Kaufman Efficiency Ratio. Возвращает np.ndarray [0, 1]."""
    close = df['Close'].values
    n = len(close)
    er = np.zeros(n)
    for i in range(period, n):
        total_move = abs(close[i] - close[i - period])
        changes = np.abs(np.diff(close[i - period:i + 1]))
        noise = np.sum(changes)
        er[i] = total_move / noise if noise > 1e-10 else 0.0
    return np.clip(er, 0, 1)


def add_efficiency_ratio(df: pd.DataFrame, period: int = 10) -> pd.DataFrame:
    """Kaufman Efficiency Ratio — добавляет efficiency_ratio в DataFrame."""
    df = df.copy()
    df['efficiency_ratio'] = compute_efficiency_ratio_array(df, period)
    return df


def compute_choppiness_array(df: pd.DataFrame, period: int = 14) -> np.ndarray:
    """Choppiness Index. Возвращает np.ndarray [0, 1].
    
    v12.7: Vectorized — uses pandas rolling sum instead of nested loop.
    O(n) instead of O(n × period).
    """
    high = df['High'].values
    low = df['Low'].values
    close = df['Close'].values
    n = len(high)

    # True Range (vectorized)
    close_back = np.zeros(n); close_back[1:] = close[:-1]
    tr = np.maximum.reduce([
        high - low,
        np.abs(high - close_back),
        np.abs(low - close_back),
    ])

    # TR rolling sum (vectorized)
    tr_sum = pd.Series(tr).rolling(period, min_periods=1).sum().values

    # HH/LL rolling (vectorized)
    high_series = pd.Series(high)
    low_series = pd.Series(low)
    hh = high_series.rolling(period, min_periods=1).max().values
    ll = low_series.rolling(period, min_periods=1).min().values

    # Choppiness formula
    with np.errstate(divide='ignore', invalid='ignore'):
        valid_mask = ((hh - ll) > 1e-10) & (tr_sum > 1e-10)
        chop = np.where(
            valid_mask,
            100 * np.log10(tr_sum / (hh - ll)) / np.log10(period),
            50
        )

    return np.clip(chop / 100, 0, 1)


def add_choppiness(df: pd.DataFrame, period: int = 14) -> pd.DataFrame:
    """Choppiness Index — добавляет choppiness в DataFrame."""
    df = df.copy()
    df['choppiness'] = compute_choppiness_array(df, period)
    return df


def engineer_features(df: pd.DataFrame, skip_temporal: bool = False) -> pd.DataFrame:
    """
    Добавляет все технические индикаторы.

    Parameters
    ----------
    skip_temporal : bool
        Если True, пропускает временные признаки (hour, day_of_week и т.д.).
        Рекомендуется для D1/W1, где час всегда 0, день недели не информативен.
    """
    df = df.copy()
    df = add_returns(df)
    df = add_sma(df)
    df = add_ema(df)
    df = add_rsi(df)
    df = add_bollinger(df)
    df = add_atr(df)
    df = add_volume_features(df)
    df = add_interaction_features(df)
    df = add_kst(df)
    df = add_efficiency_ratio(df)
    df = add_choppiness(df)
    # ADX (Wilders) — сила и направление тренда
    df = add_adx(df)
    # Волатильностные метрики (Realized, Park/GK/YZ)
    df = add_realized_volatility(df)
    df = add_parkinson_volatility(df)
    df = add_garman_klass_volatility(df)
    df = add_yang_zhang_volatility(df)
    # Дополнительные периоды для vol-экспертов (period=10)
    df = add_parkinson_volatility(df, period=10)   # → parkinson_vol_10
    df = add_garman_klass_volatility(df, period=10)  # → garman_klass_vol_10
    df = add_yang_zhang_volatility(df, period=10)  # → yang_zhang_vol_10
    # SMA 200 / EMA 200 для режимного фильтра MoE (не входят в FEATURE_COLS)
    df['sma_200'] = df['Close'].rolling(window=200).mean()
    df['ema_200'] = df['Close'].ewm(span=200, adjust=False).mean()
    # v13 features (audit 2026-08-03): Hurst exponent + VWAP distance
    df = add_hurst(df)
    df = add_vwap_features(df)
    if not skip_temporal:
        df = add_temporal_features(df)
    return df