import pandas as pd
import numpy as np
from scipy.stats import chi2_contingency
import statsmodels.api as sm

class MaritimeFatigueAnalyzer:
    """
    海事疲劳盲区分析器
    
    数据源: Paris/Tokyo/Caribbean MoU PSC检查记录
    时间范围: 2021-2025
    伤亡数: 1,022起
    
    分析方法: 逻辑回归，调整缺陷数/滞留状态/船型/区域
    """
    
    def __init__(self):
        # PSC检查数据结构
        self.inspection_data = pd.DataFrame({
            'vessel_imo': [],         # IMO编号
            'inspection_region': [],   # Paris/Tokyo/Caribbean
            'ship_type': [],           # 船型
            'deficiency_count': [],    # 缺陷总数
            'detained': [],            # 是否滞留 (0/1)
            'bnwas_deficiency': [],     # BNWAS缺陷 (0/1)
            'stcw_citation': [],       # STCW Ch.VIII引用 (0/1)
            'mlc_citation': [],        # MLC Reg 2.3引用 (0/1)
            'casualty_12m': [],        # 12月内伤亡 (0/1)
        })
        
    def logistic_regression_analysis(self, df: pd.DataFrame) -> dict:
        """
        逻辑回归分析：BNWAS缺陷对12月内伤亡概率的影响
        
        论文核心方法: 调整缺陷数/滞留/船型/区域后的逻辑回归
        """
        # 自变量
        X = df[['bnwas_deficiency', 'deficiency_count', 'detained', 
                'ship_type', 'inspection_region']]
        X = pd.get_dummies(X, columns=['ship_type', 'inspection_region'], 
                          drop_first=True)
        X = sm.add_constant(X)
        
        # 因变量
        y = df['casualty_12m']
        
        # 逻辑回归
        model = sm.Logit(y, X)
        result = model.fit()
        
        # 提取BNWAS的OR和CI
        bnwas_idx = list(X.columns).index('bnwas_deficiency')
        coef = result.params[bnwas_idx]
        se = result.bse[bnwas_idx]
        
        or_point = np.exp(coef)
        or_lower = np.exp(coef - 1.96 * se)
        or_upper = np.exp(coef + 1.96 * se)
        
        return {
            'odds_ratio': or_point,
            'ci_lower': or_lower,
            'ci_upper': or_upper,
            'p_value': result.pvalues[bnwas_idx],
            'interpretation': f"BNWAS缺陷→12月内伤亡 OR={or_point:.2f} "
                            f"(95% CI: {or_lower:.2f}-{or_upper:.2f})"
        }
    
    def regional_equity_analysis(self, df: pd.DataFrame) -> dict:
        """
        区域执法公平性分析
        
        论文发现: Caribbean MoU引用率仅为Paris的45%
        """
        regions = ['Paris', 'Tokyo', 'Caribbean']
        results = {}
        
        for region in regions:
            region_data = df[df['inspection_region'] == region]
            total_inspections = len(region_data)
            fatigue_citations = region_data['bnwas_deficiency'].sum() + \
                               region_data['stcw_citation'].sum()
            
            rate_per_1000 = fatigue_citations / total_inspections * 1000
            results[region] = {
                'total_inspections': total_inspections,
                'fatigue_citations': int(fatigue_citations),
                'rate_per_1000': rate_per_1000,
                'enforcement_gap': 'low' if rate_per_1000 < 4 else 'adequate'
            }
        
        return results


# 核心统计结果
if __name__ == "__main__":
    print("=== 核心发现 ===")
    print()
    print("1. BNWAS保护效应:")
    print("   OR = 0.26 (95% CI: 0.08-0.81)")
    print("   → 被查出BNWAS缺陷的船，后续12月伤亡概率降低74%")
    print("   → 这不是缺陷安全，而是检查触发的监管链产生了保护")
    print()
    print("2. BNWAS缺陷基线风险:")
    print("   OR = 3.88 (95% CI: 1.23-12.20)")
    print("   → 有BNWAS缺陷的船，基线伤亡风险是正常船的3.88倍")
    print("   → 这些船最需要检查保护，但最可能逃过检查")
    print()
    print("3. STCW vs MLC引用比:")
    print("   STCW Ch.VIII: 44次/5年 (需要主观判断)")
    print("   MLC Reg 2.3: 1,414次/5年 (只需文件审查)")
    print("   比率: 1:32 → 反映检查难度差异")
    print()
    print("4. 区域执法差距:")
    print("   Paris MoU:  6.1/1000次检查")
    print("   Tokyo MoU:  7.1/1000次检查")
    print("   Caribbean:  2.8/1000次检查 → 不足Paris一半")

class FatigueMeasurementFramework:
    """
    疲劳测量框架：从代理指标到真实指标
    
    海事教训: 休息小时记录(代理指标) → 实际疲劳检测率仅1.4%
    IMS方案: 必须直接测量疲劳状态，不能用驾驶时长替代
    """
    
    # 代理指标 vs 真实指标
    PROXY_METRICS = {
        'driving_hours': '驾驶时长',
        'break_frequency': '休息频率',
        'time_of_day': '时段',
        'days_since_last_sleep': '距上次睡眠'
    }
    
    DIRECT_METRICS = {
        'perclos': 'PERCLOS闭眼比例(60s窗)',
        'blink_duration': '眨眼时长(ms)',
        'hrv_rmssd': '心率变异性RMSSD(ms)',
        'breathing_rate': '呼吸率(bpm)',
        'reaction_time': '反应时间(ms)',
        'gaze_entropy': '注视熵(散乱度)'
    }
    
    # IMS必须用真实指标，不能用代理指标
    REQUIRED_METRICS = ['perclos', 'hrv_rmssd', 'breathing_rate']
    
    def evaluate_fatigue_system(self, system_metrics: dict) -> dict:
        """评估DMS系统是否直接测量疲劳"""
        has_direct = any(m in system_metrics for m in self.DIRECT_METRICS)
        only_proxy = all(m in self.PROXY_METRICS for m in system_metrics)
        
        if only_proxy:
            return {
                'rating': 'insufficient',
                'warning': '海事教训: 仅代理指标无法检测真实疲劳',
                'recommendation': '增加PERCLOS或HRV等直接测量'
            }
        elif has_direct:
            return {
                'rating': 'adequate',
                'direct_metrics': [m for m in system_metrics if m in self.DIRECT_METRICS]
            }
        else:
            return {'rating': 'unknown', 'system_metrics': system_metrics}