import numpy as np
import pandas as pd
from sklearn.cluster import KMeans
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import LeaveOneOut, cross_val_score
from scipy.stats import welch_t_test, ttest_ind
from statsmodels.stats.multitest import multipletests

class MCIDrivingDetector:
    """
    MCI驾驶行为检测器
    
    论文: PMC Sensors, 2025
    数据: 51名驾驶员（10 MCI, 41 normal）
    特征: GPS+IMU+OBD-II遥测数据
    方法: 聚类+RF排序+L1逻辑回归+LOOCV
    
    最强标志物: 油门位置变异性（Cohen's d=0.86）
    """
    
    # 论文提取的驾驶行为特征
    FEATURES = [
        # 速度相关
        'mean_speed', 'speed_var', 'max_speed', 'speed_percentile_85',
        # 油门相关（最强预测器）
        'throttle_mean', 'throttle_var', 'throttle_jerk',
        # 刹车相关
        'brake_count', 'brake_force_mean', 'hard_brake_count',
        # 转向相关
        'steering_var', 'lane_deviation',
        # 行程模式
        'night_trip_ratio', 'trip_distance_mean', 'trip_count_per_day',
        # 位置模式
        'location_entropy', 'radius_of_gyration',
        # 时间模式
        'driving_time_window', 'weekend_ratio'
    ]
    
    def __init__(self):
        self.clustering = KMeans(n_clusters=3, random_state=42)
        self.rf = RandomForestClassifier(n_estimators=200, random_state=42)
        self.lr = LogisticRegression(
            penalty='l1', solver='liblinear', C=0.1, random_state=42
        )
        
    def feature_analysis(self, X: np.ndarray, y: np.ndarray) -> dict:
        """
        完整特征分析管道
        
        Args:
            X: 特征矩阵, shape=(n_participants, n_features)
            y: 标签, 0=Normal, 1=MCI
        """
        results = {}
        
        # 1. 聚类分析
        cluster_labels = self.clustering.fit_predict(X)
        results['clusters'] = cluster_labels
        
        # 2. Random Forest特征排序
        self.rf.fit(X, y)
        importance = self.rf.feature_importances_
        results['feature_importance'] = dict(zip(self.FEATURES, importance))
        
        # 3. Welch t-test + BH校正
        p_values = []
        effect_sizes = []
        for i, feat in enumerate(self.FEATURES):
            normal_vals = X[y == 0, i]
            mci_vals = X[y == 1, i]
            
            _, p_val = ttest_ind(normal_vals, mci_vals, equal_var=False)
            
            # Cohen's d
            pooled_std = np.sqrt(
                ((len(normal_vals)-1)*np.var(normal_vals, ddof=1) +
                 (len(mci_vals)-1)*np.var(mci_vals, ddof=1)) /
                (len(normal_vals) + len(mci_vals) - 2)
            )
            d = (np.mean(mci_vals) - np.mean(normal_vals)) / (pooled_std + 1e-10)
            
            p_values.append(p_val)
            effect_sizes.append(d)
        
        # Benjamini-Hochberg FDR校正
        _, corrected_p, _, _ = multipletests(p_values, method='fdr_bh')
        
        results['effect_sizes'] = dict(zip(self.FEATURES, effect_sizes))
        results['p_values'] = dict(zip(self.FEATURES, p_values))
        results['corrected_p'] = dict(zip(self.FEATURES, corrected_p))
        
        # 4. L1逻辑回归 + LOOCV
        loo = LeaveOneOut()
        scores = cross_val_score(self.lr, X, y, cv=loo, scoring='roc_auc')
        results['loocv_auc'] = scores.mean()
        results['loocv_auc_std'] = scores.std()
        
        return results
    
    def cold_start_analysis(self, X: np.ndarray, y: np.ndarray, 
                            max_trips: int = 200) -> dict:
        """
        冷启动分析：最少需要多少次行程
        
        论文发现: ~50次行程≈4个月自然驾驶
        """
        trip_counts = [5, 10, 20, 30, 40, 50, 75, 100, 150, 200]
        auc_by_trips = {}
        
        for n in trip_counts:
            # 模拟前n次行程的累积数据
            # 实际应按时间顺序取前n次
            subset_X = X[:n]  # 简化
            subset_y = y[:n] if len(y) >= n else y
            
            if len(np.unique(subset_y)) < 2:
                auc_by_trips[n] = None
                continue
            
            loo = LeaveOneOut()
            scores = cross_val_score(self.lr, subset_X, subset_y, 
                                     cv=min(loo.get_n_splits(subset_X), 5),
                                     scoring='roc_auc')
            auc_by_trips[n] = scores.mean()
        
        return auc_by_trips


# 实际测试
if __name__ == "__main__":
    np.random.seed(42)
    
    # 模拟51名参与者数据
    n_normal, n_mci = 41, 10
    n_features = 15
    
    # 正常组
    X_normal = np.random.normal(0, 1, (n_normal, n_features))
    # MCI组：油门变异性和均值偏移
    X_mci = np.random.normal(0, 1, (n_mci, n_features))
    X_mci[:, 5] += 0.86  # throttle_var (Cohen's d=0.86)
    X_mci[:, 4] += 0.86  # throttle_mean
    X_mci[:, 11] += 0.69  # night_trip_ratio (d=0.69)
    
    X = np.vstack([X_normal, X_mci])
    y = np.array([0]*n_normal + [1]*n_mci)
    
    detector = MCIDrivingDetector()
    results = detector.feature_analysis(X, y)
    
    print("=== 特征重要性 Top 5 ===")
    importance = sorted(results['feature_importance'].items(), 
                       key=lambda x: -x[1])
    for feat, imp in importance[:5]:
        d = results['effect_sizes'][feat]
        print(f"  {feat}: importance={imp:.4f}, Cohen's d={d:.2f}")
    
    print(f"\n=== LOOCV AUC ===")
    print(f"  AUC = {results['loocv_auc']:.3f} ± {results['loocv_auc_std']:.3f}")
    
    print(f"\n=== 论文结果 ===")
    print(f"  AUC = 0.698 (95% CI: 0.493-0.872)")
    print(f"  Sensitivity = 0.80")
    print(f"  Specificity = 0.66")
    print(f"  最强标志物: 油门位置变异性 (d=0.86)")
    print(f"  冷启动: ~50次行程≈4个月")
    print(f"  性别混淆: 排除后AUC降至0.598")

class ThrottleVariabilityMonitor:
    """
    油门变异性监测器
    
    论文发现: 油门变异性是MCI最强行为标志物
    IMS应用: 作为认知状态的辅助指标
    
    注意: 油门变异性也受路况、交通、驾驶风格影响
    必须建立个体基线并排除环境因素
    """
    
    def __init__(self, baseline_window_days=30):
        self.baseline_days = baseline_window_days
        self.throttle_history = []
        
    def compute_throttle_variability(self, obd_data: np.ndarray) -> float:
        """
        计算油门位置变异性
        
        Args:
            obd_data: OBD-II油门位置序列 (%)
            
        Returns:
            variability: 变异性指标
        """
        # 标准差
        std = np.std(obd_data)
        
        # 变异系数（CV）
        mean = np.mean(obd_data)
        cv = std / (mean + 1e-10)
        
        # 频域熵（操作频率的复杂性）
        from scipy.signal import welch
        freqs, psd = welch(obd_data, fs=10)  # 10Hz采样
        psd_norm = psd / (np.sum(psd) + 1e-10)
        entropy = -np.sum(psd_norm * np.log(psd_norm + 1e-10))
        
        return {
            'std': std,
            'cv': cv,
            'spectral_entropy': entropy,
            'jerk_mean': np.mean(np.abs(np.diff(obd_data)))
        }
    
    def assess_cognitive_risk(self, current_variability: dict, 
                               baseline: dict) -> dict:
        """
        评估认知风险
        
        策略: 与个体基线比较，而非绝对阈值
        """
        # 标准化偏差
        std_z = (current_variability['std'] - baseline['std']) / (baseline['std'] + 1e-10)
        cv_z = (current_variability['cv'] - baseline['cv']) / (baseline['cv'] + 1e-10)
        entropy_z = (current_variability['spectral_entropy'] - baseline['spectral_entropy']) / \
                    (baseline['spectral_entropy'] + 1e-10)
        
        # 综合风险
        risk_score = 0.4 * std_z + 0.3 * cv_z + 0.3 * entropy_z
        
        if risk_score > 1.5:
            return {'risk': 'high', 'score': risk_score, 'recommendation': '建议认知评估'}
        elif risk_score > 0.8:
            return {'risk': 'moderate', 'score': risk_score, 'recommendation': '持续监测'}
        else:
            return {'risk': 'low', 'score': risk_score, 'recommendation': '正常'}