"""
⚡ ULTRON Systems // Stationary Feature Engineering Engine for OHLCV Bars
Transforms non-stationary price levels into normalized, stationary statistical distributions:
- Multi-scale log returns
- Intra-bar candle morphology (body, wicks, range)
- Parkinson & Garman-Klass micro-volatility
- Normalized volume dynamics
- Multi-horizon forward return & directional targets
"""
import numpy as np
import pandas as pd
from typing import Tuple, List

def compute_stationary_features(df: pd.DataFrame, pred_horizons: List[int] = [1, 3, 5, 10, 15]) -> Tuple[np.ndarray, np.ndarray]:
    """
    Computes 8 stationary features and multi-horizon target returns from raw OHLCV.
    Input df columns: ['timestamp', 'open', 'high', 'low', 'close', 'volume', 'ticker']
    Returns:
        features: np.ndarray [N_bars, 8], float32
        targets:  np.ndarray [N_bars, len(pred_horizons)], float32 (forward log returns)
    """
    # Sort strictly by timestamp
    df = df.sort_values("timestamp").reset_index(drop=True)
    
    close = df["close"].values.astype(np.float64)
    open_p = df["open"].values.astype(np.float64)
    high = df["high"].values.astype(np.float64)
    low = df["low"].values.astype(np.float64)
    volume = df["volume"].values.astype(np.float64)
    
    eps = 1e-7
    n_bars = len(close)
    if n_bars < 60:
        return np.empty((0, 8), dtype=np.float32), np.empty((0, len(pred_horizons)), dtype=np.float32)
        
    # 1. Log Returns: ln(C_t / C_{t-1})
    close_prev = np.roll(close, 1)
    close_prev[0] = close[0]
    log_ret = np.log(np.maximum(close, eps) / np.maximum(close_prev, eps))
    log_ret[0] = 0.0
    
    # 2. Inter-bar Gap: (O_t - C_{t-1}) / C_{t-1}
    co_gap = (open_p - close_prev) / np.maximum(close_prev, eps)
    co_gap[0] = 0.0
    
    # 3. High-Low Spread / Range: (H_t - L_t) / C_t
    hl_range = high - low
    hl_spread = hl_range / np.maximum(close, eps)
    
    # 4. Normalized Candle Body: (C_t - O_t) / (H_t - L_t + eps) -> range in [-1, 1]
    candle_body = (close - open_p) / (hl_range + eps)
    
    # 5. Upper Wick: (H_t - max(O_t, C_t)) / (H_t - L_t + eps) -> range in [0, 1]
    max_oc = np.maximum(open_p, close)
    upper_wick = (high - max_oc) / (hl_range + eps)
    
    # 6. Lower Wick: (min(O_t, C_t) - L_t) / (H_t - L_t + eps) -> range in [0, 1]
    min_oc = np.minimum(open_p, close)
    lower_wick = (min_oc - low) / (hl_range + eps)
    
    # 7. Parkinson Micro-Volatility proxy: sqrt((ln(H/L))^2 / (4 * ln(2)))
    hl_ratio = np.maximum(high, eps) / np.maximum(low, eps)
    parkinson_vol = np.sqrt(np.maximum(0.0, np.log(hl_ratio)**2 / (4.0 * np.log(2.0))))
    
    # 8. Normalized Volume Z-Score relative to rolling window
    log_vol = np.log1p(np.maximum(volume, 0.0))
    vol_series = pd.Series(log_vol)
    vol_mean = vol_series.rolling(window=30, min_periods=1).mean().values
    vol_std = vol_series.rolling(window=30, min_periods=1).std().fillna(1.0).values
    vol_zscore = (log_vol - vol_mean) / np.maximum(vol_std, 1e-4)
    
    # Stack features: shape [N, 8]
    feats = np.column_stack([
        log_ret,        # Feature 0: 1-bar log return
        co_gap,         # Feature 1: Overnight / inter-bar gap
        hl_spread,      # Feature 2: Relative bar range
        candle_body,    # Feature 3: Directional body strength [-1, 1]
        upper_wick,     # Feature 4: Overhead rejection wick [0, 1]
        lower_wick,     # Feature 5: Underlying absorption wick [0, 1]
        parkinson_vol,  # Feature 6: Intraday continuous volatility
        vol_zscore      # Feature 7: Relative volume burst
    ])
    
    # Clip extreme values for numerical stability (prevent gradient explosions)
    feats = np.clip(feats, -10.0, 10.0).astype(np.float32)
    
    # Compute multi-horizon forward returns: ln(C_{t + h} / C_t)
    targets_list = []
    max_horizon = max(pred_horizons)
    
    for h in pred_horizons:
        fwd_close = np.roll(close, -h)
        fwd_ret = np.log(np.maximum(fwd_close, eps) / np.maximum(close, eps))
        # Mask out boundary wrap-around
        fwd_ret[-h:] = np.nan
        targets_list.append(fwd_ret)
        
    targets = np.column_stack(targets_list).astype(np.float32)
    
    # Remove final `max_horizon` rows where future is unknown
    valid_len = n_bars - max_horizon
    if valid_len <= 0:
        return np.empty((0, 8), dtype=np.float32), np.empty((0, len(pred_horizons)), dtype=np.float32)
        
    return feats[:valid_len], targets[:valid_len]