import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
from scipy.signal import butter, filtfilt, find_peaks, savgol_filter, medfilt, welch
from scipy.stats import kurtosis, skew
from typing import Tuple, Dict, Optional
import struct
import warnings
import os
from matplotlib.backends.backend_pdf import PdfPages
from datetime import datetime

class EMCtoECGConverter:
    def __init__(self, sample_rate: float = 250.0):
        """
        Inicializa o conversor EMC para ECG com configurações otimizadas.
        """
        self.sample_rate = sample_rate
        self.qrs_width = int(0.1 * sample_rate)  # Largura típica do QRS em amostras
        self.min_bpm = 30
        self.max_bpm = 200
        self.min_rr = 60 / self.max_bpm
        self.max_rr = 60 / self.min_bpm

    def read_emc_file(self, filename: str) -> np.ndarray:
        """
        Lê arquivo EMC com validação extensiva dos dados.
        """
        if not os.path.exists(filename):
            raise FileNotFoundError(f"Arquivo não encontrado: {filename}")

        try:
            formats_to_try = [
                ('float32', 4),
                ('int16', 2),
                ('int32', 4),
                ('float64', 8)
            ]

            best_data = None
            best_quality = -1

            for dtype, _ in formats_to_try:
                try:
                    data = np.fromfile(filename, dtype=dtype)
                    if len(data) > self.sample_rate * 5:  # Mínimo 5 segundos
                        quality = self.evaluate_signal_quality(data)
                        if quality > best_quality:
                            best_quality = quality
                            best_data = data
                except:
                    continue

            if best_data is not None and best_quality > 0.5:  # Qualidade mínima aceitável
                return best_data
            raise ValueError("Dados não atingem qualidade mínima necessária")

        except Exception as e:
            raise ValueError(f"Erro ao ler arquivo EMC: {str(e)}")

    def evaluate_signal_quality(self, data: np.ndarray) -> float:
        """
        Avalia a qualidade do sinal com múltiplos critérios.
        """
        try:
            # 1. Verificações básicas
            if len(data) < self.sample_rate * 5:
                return 0
            if not np.all(np.isfinite(data)):
                return 0
            if np.std(data) < 1e-6:
                return 0

            # 2. Análise espectral
            freqs, psd = welch(data, fs=self.sample_rate)
            power_cardiac = np.sum(psd[(freqs >= 0.5) & (freqs <= 40)])
            power_total = np.sum(psd)
            spectral_ratio = power_cardiac / power_total if power_total > 0 else 0

            # 3. Análise estatística
            kurt = abs(kurtosis(data))
            skewness = abs(skew(data))

            # 4. Detecção de descontinuidades
            diff = np.diff(data)
            discontinuities = np.sum(np.abs(diff) > 5 * np.std(diff))
            continuity_score = 1 - (discontinuities / len(data))

            # 5. SNR estimado
            noise_estimate = np.std(diff) / np.sqrt(2)
            signal_power = np.std(data)
            snr = 20 * np.log10(signal_power/noise_estimate) if noise_estimate > 0 else 0
            snr_score = min(max(snr / 40.0, 0), 1)

            # Combinação dos critérios
            quality_score = (
                0.3 * spectral_ratio +
                0.2 * (1 / (1 + kurt/10)) +
                0.2 * (1 / (1 + skewness/5)) +
                0.15 * continuity_score +
                0.15 * snr_score
            )

            return quality_score

        except:
            return 0

    def advanced_preprocess_signal(self, data: np.ndarray) -> np.ndarray:
        """
        Pré-processamento avançado com múltiplas etapas de limpeza.
        """
        if not isinstance(data, np.ndarray) or len(data) == 0:
            raise ValueError("Dados inválidos")

        # 1. Remoção de valores inválidos e interpolação
        invalid_mask = ~np.isfinite(data)
        if np.any(invalid_mask):
            x = np.arange(len(data))
            data = np.interp(x, x[~invalid_mask], data[~invalid_mask])

        # 2. Remoção de tendência usando detrend
        data = signal.detrend(data)

        # 3. Filtragem adaptativa para remover ruído de baixa frequência
        window_size = int(2 * self.sample_rate)  # 2 segundos
        if window_size % 2 == 0:
            window_size += 1
        baseline = signal.savgol_filter(data, window_size, 3)
        data = data - baseline

        # 4. Remoção de outliers com método Hampel
        window = int(0.2 * self.sample_rate)
        k = 3  # Número de desvios padrão

        rolling_median = np.array([np.median(data[max(0, i-window):min(len(data), i+window)])
                                 for i in range(len(data))])
        rolling_std = np.array([np.std(data[max(0, i-window):min(len(data), i+window)])
                              for i in range(len(data))])

        outlier_mask = np.abs(data - rolling_median) > (k * rolling_std)
        data[outlier_mask] = rolling_median[outlier_mask]

        # 5. Suavização adaptativa
        window_size = int(0.05 * self.sample_rate)
        if window_size % 2 == 0:
            window_size += 1

        # Aplica filtro mediano para remover ruído impulsivo
        data = medfilt(data, window_size)

        # Aplica Savitzky-Golay para preservar morfologia
        data = savgol_filter(data, window_size, 3)

        # 6. Normalização por segmentos com sobreposição
        segment_size = self.sample_rate  # 1 segundo
        overlap = segment_size // 2
        normalized_data = np.zeros_like(data)
        window = np.hanning(segment_size)

        for i in range(0, len(data)-segment_size, overlap):
            segment = data[i:i+segment_size]
            # Remove DC offset
            segment = segment - np.mean(segment)
            # Normalização
            if np.std(segment) > 1e-6:
                segment = segment / np.std(segment)

            # Aplica janela de sobreposição suave
            if i == 0:
                normalized_data[i:i+segment_size] = segment
            else:
                # Mistura suave com sobreposição
                normalized_data[i:i+segment_size] = (
                    normalized_data[i:i+segment_size] * (1 - window) +
                    segment * window
                )

        # 7. Filtragem final para suavização
        b, a = butter(4, [0.5/(self.sample_rate/2), 40.0/(self.sample_rate/2)], btype='band')
        normalized_data = filtfilt(b, a, normalized_data)

        return normalized_data

    def enhanced_bandpass_filter(self, data: np.ndarray) -> np.ndarray:
        """
        Filtragem passa-banda aprimorada com múltiplas etapas.
        """
        nyquist = 0.5 * self.sample_rate

        # 1. Remoção de ruído de alta frequência
        b_high, a_high = butter(4, 40.0/nyquist, btype='low')
        data = filtfilt(b_high, a_high, data)

        # 2. Remoção de baseline wander
        b_low, a_low = butter(4, 0.5/nyquist, btype='high')
        data = filtfilt(b_low, a_low, data)

        # 3. Filtros notch para frequências específicas
        for freq in [50.0, 60.0]:  # Frequências de rede comuns
            b_notch, a_notch = butter(4, [(freq-0.5)/nyquist, (freq+0.5)/nyquist], btype='bandstop')
            data = filtfilt(b_notch, a_notch, data)

        # 4. Filtro passa-banda específico para ECG
        b_ecg, a_ecg = butter(4, [5.0/nyquist, 30.0/nyquist], btype='band')
        data = filtfilt(b_ecg, a_ecg, data)

        return data

    def robust_qrs_detection(self, data: np.ndarray) -> np.ndarray:
        """
        Detecção robusta de complexos QRS com ajustes otimizados.
        """
        # 1. Preparação do sinal
        # Aplica um filtro passa-banda específico para QRS
        nyquist = 0.5 * self.sample_rate
        b_qrs, a_qrs = butter(4, [5.0/nyquist, 15.0/nyquist], btype='band')
        filtered_qrs = filtfilt(b_qrs, a_qrs, data)

        # Derivada para realçar QRS
        diff = np.gradient(filtered_qrs)
        squared = diff ** 2

        # Integração móvel com janela menor
        window_size = int(0.08 * self.sample_rate)  # Reduzido para 80ms
        window = np.ones(window_size)/window_size
        integrated = np.convolve(squared, window, mode='same')

        # 2. Detecção inicial de picos com limiar mais sensível
        threshold = 0.3 * np.mean(integrated)  # Reduzido para 0.3
        min_distance = int(0.2 * self.sample_rate)  # 200ms entre picos

        peaks, properties = find_peaks(integrated,
                                     height=threshold,
                                     distance=min_distance,
                                     prominence=0.1)  # Adiciona critério de proeminência

        # 3. Validação dos picos com critérios ajustados
        valid_peaks = []
        prev_peak = -np.inf
        qrs_width = int(0.12 * self.sample_rate)  # Aumentado para 120ms

        for peak in peaks:
            if peak < qrs_width//2 or peak > len(data)-qrs_width//2:
                continue

            # Janela em torno do pico
            window_start = peak - qrs_width//2
            window_end = peak + qrs_width//2
            window = data[window_start:window_end]

            # Critérios de validação ajustados
            amplitude = np.abs(np.max(window) - np.min(window))
            slope = np.max(np.abs(np.diff(window)))

            # Calcula largura do QRS
            above_half = window > (np.max(window) + np.min(window))/2
            width = np.sum(above_half) / self.sample_rate * 1000

            # Verificações com limiares mais flexíveis
            amplitude_ok = amplitude > 0.3 * np.std(data)  # Reduzido para 0.3
            slope_ok = slope > 0.2 * np.std(np.diff(data))  # Reduzido para 0.2
            width_ok = 40 < width < 160  # Range mais amplo

            # Intervalo R-R com tolerância maior
            if len(valid_peaks) > 0:
                rr_interval = (peak - prev_peak) / self.sample_rate
                rr_ok = 0.25 < rr_interval < 2.5  # Range mais amplo
            else:
                rr_ok = True

            if amplitude_ok and slope_ok and width_ok and rr_ok:
                valid_peaks.append(peak)
                prev_peak = peak

        # 4. Refinamento final com janela maior
        refined_peaks = []
        search_window = int(0.05 * self.sample_rate)  # 50ms

        for peak in valid_peaks:
            window_start = max(0, peak - search_window)
            window_end = min(len(data), peak + search_window)
            window = data[window_start:window_end]
            max_pos = np.argmax(np.abs(window))
            refined_peak = window_start + max_pos
            refined_peaks.append(refined_peak)

        # 5. Verificação final de intervalos RR
        if len(refined_peaks) > 2:
            rr_intervals = np.diff(refined_peaks) / self.sample_rate
            mean_rr = np.mean(rr_intervals)
            std_rr = np.std(rr_intervals)

            # Remove picos que geram intervalos RR muito discrepantes
            valid_rr = []
            for i in range(len(refined_peaks)-1):
                rr = (refined_peaks[i+1] - refined_peaks[i]) / self.sample_rate
                if abs(rr - mean_rr) <= 3 * std_rr:  # Regra dos 3 sigma
                    if not valid_rr or refined_peaks[i] > valid_rr[-1]:
                        valid_rr.append(refined_peaks[i])

            if valid_rr:
                valid_rr.append(refined_peaks[-1])
            refined_peaks = valid_rr

        return np.array(refined_peaks)

    def calculate_signal_quality_metrics(self, data: np.ndarray, qrs_peaks: np.ndarray) -> Dict:
        """
        Calcula métricas de qualidade do sinal processado.
        """
        metrics = {}

        try:
            # 1. SNR
            signal_windows = []
            noise_windows = []

            for peak in qrs_peaks:
                if peak < self.qrs_width or peak >= len(data) - self.qrs_width:
                    continue
                # Janela do sinal (QRS)
                signal_windows.append(data[peak-self.qrs_width//2:peak+self.qrs_width//2])
                # Janelas de ruído (entre QRS)
                if len(signal_windows) > 1:
                    prev_peak = qrs_peaks[qrs_peaks < peak][-1]
                    mid_point = prev_peak + (peak - prev_peak)//2
                    noise_windows.append(data[mid_point-self.qrs_width//2:mid_point+self.qrs_width//2])

            if signal_windows and noise_windows:
                signal_power = np.mean([np.var(w) for w in signal_windows])
                noise_power = np.mean([np.var(w) for w in noise_windows])
                metrics['snr'] = 10 * np.log10(signal_power/noise_power) if noise_power > 0 else float('inf')
            else:
                metrics['snr'] = 0

            # 2. Variabilidade RR
            if len(qrs_peaks) > 1:
                rr_intervals = np.diff(qrs_peaks) / self.sample_rate
                metrics['rr_mean'] = np.mean(rr_intervals)
                metrics['rr_std'] = np.std(rr_intervals)
                metrics['rr_cv'] = metrics['rr_std'] / metrics['rr_mean'] * 100
            else:
                metrics['rr_mean'] = 0
                metrics['rr_std'] = 0
                metrics['rr_cv'] = 0

            # 3. Qualidade morfológica - Corrigido para garantir kernel ímpar
            kernel_size = self.sample_rate
            if kernel_size % 2 == 0:
                kernel_size += 1
            metrics['baseline_variation'] = np.std(medfilt(data, kernel_size))

            # 4. Frequência cardíaca
            if len(qrs_peaks) > 1:
                duration = len(data) / self.sample_rate
                metrics['heart_rate'] = len(qrs_peaks) * 60 / duration
            else:
                metrics['heart_rate'] = 0

        except Exception as e:
            print(f"\nAviso: Erro no cálculo de métricas: {str(e)}")
            # Valores padrão em caso de erro
            metrics = {
                'snr': 0,
                'rr_mean': 0,
                'rr_std': 0,
                'rr_cv': 0,
                'baseline_variation': 0,
                'heart_rate': 0
            }

        return metrics

    def convert_to_ecg(self, data: np.ndarray) -> Dict[str, np.ndarray]:
        """
        Converte sinal EMC para ECG com purificação máxima.
        """
        # 1. Pré-processamento avançado
        processed_data = self.advanced_preprocess_signal(data)

        # 2. Filtragem aprimorada
        filtered_data = self.enhanced_bandpass_filter(processed_data)

        # 3. Detecção robusta de QRS
        qrs_peaks = self.robust_qrs_detection(filtered_data)

        # 4. Análise de qualidade
        quality_metrics = self.calculate_signal_quality_metrics(filtered_data, qrs_peaks)

        return {
            'ecg_signal': filtered_data,
            'qrs_peaks': qrs_peaks,
            'quality_metrics': quality_metrics
        }

    def plot_five_windows(self, emc_data: np.ndarray, ecg_data: Dict[str, np.ndarray],
                         pdf_pages: PdfPages) -> None:
        """
        Plota 5 janelas de 10 segundos do sinal ECG purificado.
        """
        # Configuração da figura
        fig = plt.figure(figsize=(15, 25))
        fig.suptitle('Análise de ECG Purificado - 5 Janelas de 10 Segundos', fontsize=16)

        window_duration = 10.0  # segundos
        samples_per_window = int(window_duration * self.sample_rate)

        for i in range(5):
            # Calcula índices para a janela atual
            start_idx = i * samples_per_window
            end_idx = start_idx + samples_per_window

            if start_idx >= len(emc_data):
                break

            # Ajusta end_idx se necessário
            end_idx = min(end_idx, len(emc_data))

            # Cria subplot
            ax = plt.subplot(5, 1, i+1)

            # Tempo para este segmento
            time = np.arange(start_idx, end_idx) / self.sample_rate

            # Plota o sinal ECG
            ax.plot(time, ecg_data['ecg_signal'][start_idx:end_idx], 'b-',
                   label='ECG Purificado', linewidth=1)

            # Marca os picos QRS nesta janela
            mask = (ecg_data['qrs_peaks'] >= start_idx) & (ecg_data['qrs_peaks'] < end_idx)
            if np.any(mask):
                qrs_times = ecg_data['qrs_peaks'][mask] / self.sample_rate
                qrs_amplitudes = ecg_data['ecg_signal'][ecg_data['qrs_peaks'][mask]]
                ax.plot(qrs_times, qrs_amplitudes, 'ro',
                       label='Complexos QRS', markersize=8)

                # Adiciona intervalos R-R e qualidade
                for j in range(len(qrs_times)-1):
                    rr_interval = qrs_times[j+1] - qrs_times[j]
                    midpoint_time = (qrs_times[j] + qrs_times[j+1]) / 2
                    midpoint_amp = np.max(ecg_data['ecg_signal'][start_idx:end_idx]) * 1.1

                    # Calcula qualidade do intervalo
                    segment = ecg_data['ecg_signal'][int(qrs_times[j]*self.sample_rate):
                                                    int(qrs_times[j+1]*self.sample_rate)]
                    quality = self.evaluate_signal_quality(segment)

                    # Cor baseada na qualidade
                    color = 'green' if quality > 0.8 else 'orange' if quality > 0.6 else 'red'

                    ax.annotate(f'{rr_interval:.2f}s (Q:{quality:.2f})',
                              xy=(midpoint_time, midpoint_amp),
                              ha='center', va='bottom',
                              color=color)

            # Adiciona grade de referência
            ax.grid(True, which='both', alpha=0.3)
            ax.grid(True, which='major', alpha=0.5)

            # Configurações do gráfico
            ax.set_title(f'Janela {i+1}: {i*10:.1f}-{(i+1)*10:.1f} segundos', fontsize=12)
            ax.set_xlabel('Tempo (s)', fontsize=10)
            ax.set_ylabel('Amplitude Normalizada', fontsize=10)

            # Adiciona indicadores de qualidade
            if 'quality_metrics' in ecg_data:
                quality_text = f"SNR: {ecg_data['quality_metrics'].get('snr', 'N/A'):.1f}dB"
                ax.text(0.02, 0.98, quality_text,
                       transform=ax.transAxes,
                       verticalalignment='top',
                       fontsize=8)

            if i == 0:  # Legenda apenas no primeiro gráfico
                ax.legend(loc='upper right')

            # Define limites do eixo x para cada janela
            ax.set_xlim(i*10, (i+1)*10)

        plt.tight_layout()
        pdf_pages.savefig(fig)
        plt.close()

    def plot_quality_summary(self, ecg_data: Dict[str, np.ndarray], pdf_pages: PdfPages) -> None:
        """
        Plota resumo das métricas de qualidade do sinal.
        """
        fig = plt.figure(figsize=(10, 8))
        plt.axis('off')

        if 'quality_metrics' in ecg_data:
            metrics = ecg_data['quality_metrics']
            summary_text = (
                "Métricas de Qualidade do Sinal\n\n"
                f"SNR: {metrics.get('snr', 'N/A'):.1f} dB\n"
                f"Frequência Cardíaca Média: {metrics.get('heart_rate', 'N/A'):.1f} bpm\n"
                f"Intervalo R-R Médio: {metrics.get('rr_mean', 'N/A'):.3f} s\n"
                f"Variabilidade R-R (CV): {metrics.get('rr_cv', 'N/A'):.1f}%\n"
                f"Variação da Linha de Base: {metrics.get('baseline_variation', 'N/A'):.3f}\n"
            )
            plt.text(0.1, 0.9, summary_text, fontsize=12, va='top')

        pdf_pages.savefig(fig)
        plt.close()

    def process_file(self, filename: str) -> None:
        """
        Processa arquivo EMC e salva visualização em PDF com análise de qualidade.
        """
        try:
            # Cria nome do arquivo PDF
            base_filename = os.path.basename(os.path.splitext(filename)[0])
            pdf_filename = os.path.join('/kaggle/working/', f'{base_filename}_analise_purificada.pdf')

            print("\n=== Iniciando Processamento do ECG ===")
            print(f"Arquivo de entrada: {filename}")
            print("Configurando parâmetros iniciais...")
            print(f"- Taxa de amostragem: {self.sample_rate} Hz")
            print(f"- Largura QRS esperada: {self.qrs_width} amostras")
            print(f"- Range de BPM aceito: {self.min_bpm}-{self.max_bpm} bpm")

            # Garante que o diretório existe
            os.makedirs('/kaggle/working/', exist_ok=True)
            print(f"\nDiretório de saída configurado: /kaggle/working/")

            with PdfPages(pdf_filename) as pdf:
                # Leitura do arquivo
                print("\n[1/7] Lendo arquivo EMC...")
                print("- Testando diferentes formatos de dados")
                emc_data = self.read_emc_file(filename)
                print(f"✓ Dados lidos com sucesso")
                print(f"- Total de amostras: {len(emc_data):,}")
                print(f"- Duração do sinal: {len(emc_data)/self.sample_rate:.1f} segundos")

                # Pré-processamento
                print("\n[2/7] Iniciando pré-processamento do sinal...")
                print("- Removendo valores inválidos")
                print("- Aplicando detrend")
                print("- Removendo ruído de baixa frequência")
                print("- Removendo outliers")
                print("- Aplicando normalização")
                processed_data = self.advanced_preprocess_signal(emc_data)
                print("✓ Pré-processamento concluído")

                # Filtragem
                print("\n[3/7] Aplicando filtros...")
                print("- Filtro passa-baixa (40 Hz)")
                print("- Filtro passa-alta (0.5 Hz)")
                print("- Filtros notch (50/60 Hz)")
                print("- Filtro passa-banda ECG (5-30 Hz)")
                filtered_data = self.enhanced_bandpass_filter(processed_data)
                print("✓ Filtragem concluída")

                # Detecção QRS
                print("\n[4/7] Detectando complexos QRS...")
                print("- Calculando derivada do sinal")
                print("- Aplicando limiar adaptativo")
                print("- Validando picos candidatos")
                print("- Refinando posições dos picos")
                qrs_peaks = self.robust_qrs_detection(filtered_data)
                print(f"✓ Detecção concluída: {len(qrs_peaks)} complexos QRS encontrados")

                # Análise de qualidade
                print("\n[5/7] Calculando métricas de qualidade...")
                quality_metrics = self.calculate_signal_quality_metrics(filtered_data, qrs_peaks)
                print("- SNR calculado")
                print("- Variabilidade R-R analisada")
                print("- Qualidade morfológica avaliada")
                print(f"✓ SNR médio: {quality_metrics.get('snr', 'N/A'):.1f} dB")
                print(f"✓ Frequência cardíaca média: {quality_metrics.get('heart_rate', 'N/A'):.1f} bpm")

                # Resultados
                ecg_data = {
                    'ecg_signal': filtered_data,
                    'qrs_peaks': qrs_peaks,
                    'quality_metrics': quality_metrics
                }

                # Gerando PDF
                print("\n[6/7] Gerando relatório PDF...")
                print("- Criando página de informações")
                fig_info = plt.figure(figsize=(8.5, 11))
                plt.axis('off')

                info_text = (
                    f"Análise de ECG Purificado\n\n"
                    f"Arquivo: {os.path.basename(filename)}\n"
                    f"Data: {datetime.now().strftime('%d/%m/%Y %H:%M:%S')}\n"
                    f"Dispositivo: KT-01181\n"
                    f"Taxa de amostragem: {self.sample_rate} Hz\n"
                    f"Total de amostras: {len(emc_data):,}\n"
                    f"Duração total: {len(emc_data)/self.sample_rate:.1f} segundos\n"
                    f"Complexos QRS detectados: {len(qrs_peaks)}\n"
                    f"SNR médio: {quality_metrics.get('snr', 'N/A'):.1f} dB\n"
                    f"Frequência cardíaca média: {quality_metrics.get('heart_rate', 'N/A'):.1f} bpm\n"
                    f"Intervalo R-R médio: {quality_metrics.get('rr_mean', 'N/A'):.3f} s\n"
                    f"Variabilidade R-R: {quality_metrics.get('rr_cv', 'N/A'):.1f}%\n"
                    f"Qualidade global do sinal: {self.evaluate_signal_quality(filtered_data):.2%}\n"
                )
                plt.text(0.1, 0.9, info_text, fontsize=12, va='top')
                pdf.savefig(fig_info)
                plt.close()

                print("- Gerando visualizações das janelas temporais")
                self.plot_five_windows(emc_data, ecg_data, pdf)

                print("- Adicionando resumo de qualidade")
                self.plot_quality_summary(ecg_data, pdf)

                print("\n[7/7] Finalização...")
                print(f"✓ Processamento concluído com sucesso!")
                print(f"✓ Relatório salvo em: {pdf_filename}")
                print("\n=== Processamento Finalizado ===")

        except Exception as e:
            print("\n❌ ERRO DURANTE O PROCESSAMENTO")
            print(f"Descrição do erro: {str(e)}")
            print("\nDetalhes técnicos do erro:")
            import traceback
            traceback.print_exc()

# Exemplo de uso
if __name__ == "__main__":
    try:
        # Instancia o conversor
        converter = EMCtoECGConverter(sample_rate=250)

        # Processa o arquivo
        filename = "/kaggle/input/1728646157-8nqhdzjoei477swp-cardios/1728646157-8nqHdzjoei477swP.xcm"
        converter.process_file(filename)

    except Exception as e:
        print(f"Erro: {str(e)}")