lab2/src/eda_analysis.py - master - dakone22/ppa

lab2/src/eda_analysis.py
380 строк15 KB
dakone22
created: lab2, lab3, lab4
04 июн 2026, 14:19
04 июн 2026, 14:1903322e8
"""
Скрипт для разведочного анализа данных (EDA) и визуализации.
Версия: исправленная и дополненная (Granger Causality, Persistent Homology fix).
"""

import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')  # Use non-interactive backend to avoid Tkinter issues
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.decomposition import PCA
from sklearn.manifold import TSNE
import umap
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import RandomForestClassifier
import shap
import os
import warnings
warnings.filterwarnings('ignore')

# Настройка стиля
plt.style.use('seaborn-v0_8-darkgrid')
sns.set_palette("husl")

# Параметры
import os
BASE_DIR = os.path.dirname(os.path.abspath(__file__))
DATA_DIR = os.path.join(BASE_DIR, "..", "data", "processed")
DATASET_PATH = os.path.join(DATA_DIR, "balanced_dataset.parquet")  # Сбалансированный датасет (включает Label 1)
OUTPUT_DIR = os.path.join(BASE_DIR, "..", "plots")
os.makedirs(OUTPUT_DIR, exist_ok=True)

def load_data():
    """Загрузка и подготовка данных."""
    print("Загрузка данных...")
    df = pd.read_parquet(DATASET_PATH)
    print(f"Форма данных: {df.shape}")
    print(f"Признаки: {list(df.columns)}")
    return df

def plot_label_distribution(df):
    """Визуализация распределения классов."""
    print("Построение распределения классов...")
    plt.figure(figsize=(10, 6))
    ax = df['label'].value_counts(normalize=True).sort_index().plot(kind='bar', color='steelblue')
    plt.title('Распределение классов (Label Distribution)', fontsize=16)
    plt.xlabel('Label', fontsize=12)
    plt.ylabel('Доля от общего числа', fontsize=12)
    plt.xticks(rotation=0)

    # Добавляем проценты
    total = len(df)
    for p in ax.patches:
        percentage = f'{100 * p.get_height():.1f}%'
        x = p.get_x() + p.get_width() / 2
        y = p.get_height()
        ax.annotate(percentage, (x, y), ha='center', va='bottom', fontsize=10)

    plt.tight_layout()
    plt.savefig(f"{OUTPUT_DIR}/label_distribution.png", dpi=300)
    plt.close()
    print("  > Сохранено: label_distribution.png")

def plot_correlation_heatmap(df):
    """Тепловая карта корреляции признаков."""
    print("Построение тепловой карты корреляции...")
    # Выбираем только числовые признаки (исключая ID)
    numeric_cols = df.select_dtypes(include=[np.number]).columns.tolist()
    if 'label' in numeric_cols:
        numeric_cols.remove('label')

    # Ограничиваем количество признаков для читаемости (топ-20 по дисперсии)
    variances = df[numeric_cols].var()
    top_features = variances.nlargest(20).index.tolist()

    plt.figure(figsize=(16, 14))
    corr = df[top_features].corr()
    sns.heatmap(corr, annot=True, fmt=".2f", cmap='coolwarm', center=0,
                square=True, linewidths=.5, cbar_kws={"shrink": .8})
    plt.title('Тепловая карта корреляции признаков (Top 20)', fontsize=16)
    plt.xticks(rotation=45, ha='right')
    plt.yticks(rotation=0)
    plt.tight_layout()
    plt.savefig(f"{OUTPUT_DIR}/correlation_heatmap.png", dpi=300)
    plt.close()
    print("  > Сохранено: correlation_heatmap.png")

def plot_pca_3d(df):
    """3D визуализация PCA."""
    print("Построение PCA 3D...")
    # Сэмплируем для скорости
    sample_df = df.sample(min(50000, len(df)), random_state=42)

    numeric_cols = sample_df.select_dtypes(include=[np.number]).columns.tolist()
    if 'label' in numeric_cols:
        numeric_cols.remove('label')

    X = sample_df[numeric_cols]
    y = sample_df['label']

    # Обработка NaN и масштабирование
    X = X.fillna(0)
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)

    pca = PCA(n_components=3, random_state=42)
    X_pca = pca.fit_transform(X_scaled)

    fig = plt.figure(figsize=(12, 10))
    ax = fig.add_subplot(111, projection='3d')

    # Цвета для классов
    colors = {0: 'red', 1: 'blue', 2: 'green', 3: 'orange'}
    for label in sorted(y.unique()):
        idx = y == label
        ax.scatter(X_pca[idx, 0], X_pca[idx, 1], X_pca[idx, 2],
                   c=colors.get(label, 'gray'), label=f'Label {label}',
                   alpha=0.5, s=10)

    ax.set_xlabel('PC1')
    ax.set_ylabel('PC2')
    ax.set_zlabel('PC3')
    ax.set_title('PCA 3D Projection')
    ax.legend()
    plt.savefig(f"{OUTPUT_DIR}/pca_3d.png", dpi=300)
    plt.close()
    print("  > Сохранено: pca_3d.png")
    print(f"  > Explained variance: {pca.explained_variance_ratio_}")

def plot_tsne(df):
    """t-SNE визуализация."""
    print("Построение t-SNE...")
    sample_df = df.sample(min(8000, len(df)), random_state=42)

    numeric_cols = sample_df.select_dtypes(include=[np.number]).columns.tolist()
    if 'label' in numeric_cols:
        numeric_cols.remove('label')

    X = sample_df[numeric_cols].fillna(0)
    y = sample_df['label']

    tsne = TSNE(n_components=2, random_state=42, perplexity=30, max_iter=1000)
    X_tsne = tsne.fit_transform(X)

    plt.figure(figsize=(12, 8))
    colors = {0: 'red', 1: 'blue', 2: 'green', 3: 'orange'}
    for label in sorted(y.unique()):
        idx = y == label
        plt.scatter(X_tsne[idx, 0], X_tsne[idx, 1],
                   c=colors.get(label, 'gray'), label=f'Label {label}',
                   alpha=0.6, s=15)

    plt.title('t-SNE Projection', fontsize=16)
    plt.xlabel('t-SNE 1')
    plt.ylabel('t-SNE 2')
    plt.legend()
    plt.tight_layout()
    plt.savefig(f"{OUTPUT_DIR}/tsne.png", dpi=300)
    plt.close()
    print("  > Сохранено: tsne.png")

def plot_umap(df):
    """UMAP визуализация."""
    print("Построение UMAP...")
    sample_df = df.sample(min(8000, len(df)), random_state=42)

    numeric_cols = sample_df.select_dtypes(include=[np.number]).columns.tolist()
    if 'label' in numeric_cols:
        numeric_cols.remove('label')

    X = sample_df[numeric_cols].fillna(0)
    y = sample_df['label']

    reducer = umap.UMAP(random_state=42, n_neighbors=15, min_dist=0.1, n_jobs=1)
    X_umap = reducer.fit_transform(X)

    plt.figure(figsize=(12, 8))
    colors = {0: 'red', 1: 'blue', 2: 'green', 3: 'orange'}
    for label in sorted(y.unique()):
        idx = y == label
        plt.scatter(X_umap[idx, 0], X_umap[idx, 1],
                   c=colors.get(label, 'gray'), label=f'Label {label}',
                   alpha=0.6, s=15)

    plt.title('UMAP Projection', fontsize=16)
    plt.xlabel('UMAP 1')
    plt.ylabel('UMAP 2')
    plt.legend()
    plt.tight_layout()
    plt.savefig(f"{OUTPUT_DIR}/umap.png", dpi=300)
    plt.close()
    print("  > Сохранено: umap.png")

def plot_shap_summary(df):
    """SHAP values для интерпретации модели."""
    print("Построение SHAP summary plot...")

    # Используем сэмпл для скорости
    sample_df = df.sample(min(5000, len(df)), random_state=42)

    numeric_cols = sample_df.select_dtypes(include=[np.number]).columns.tolist()
    if 'label' in numeric_cols:
        numeric_cols.remove('label')

    X = sample_df[numeric_cols].fillna(0)
    y = sample_df['label']

    # Обучаем модель (для SHAP лучше использовать классификатор, например RF)
    model = RandomForestClassifier(n_estimators=50, random_state=42, n_jobs=-1)
    model.fit(X, y)

    explainer = shap.TreeExplainer(model)
    shap_values = explainer.shap_values(X)

    # Для многоклассовой классификации берем SHAP values для класса 0 (или усредняем)
    if isinstance(shap_values, list):
        shap_values_avg = np.abs(shap_values).mean(axis=0)  # Усредняем по классам
    else:
        shap_values_avg = np.abs(shap_values).mean(axis=0)

    plt.figure(figsize=(12, 8))
    shap.summary_plot(shap_values, X, show=False, plot_type="bar")
    plt.tight_layout()
    plt.savefig(f"{OUTPUT_DIR}/shap_summary.png", dpi=300, bbox_inches='tight')
    plt.close()
    print("  > Сохранено: shap_summary.png")

def plot_persistent_homology():
    """Persistent Homology (заглушка или попытка установки)."""
    print("Построение Persistent Homology...")
    plt.figure(figsize=(10, 6))

    try:
        # Пытаемся использовать ripser или gudhi
        from ripser import ripser
        from persim import plot_diagrams

        # Генерируем тестовые данные
        data = np.random.rand(100, 2)
        diagrams = ripser(data)['dgms']

        plt.figure(figsize=(12, 5))
        plot_diagrams(diagrams, show=False)
        plt.suptitle('Persistent Homology (Ripser)', fontsize=16)
        plt.savefig(f"{OUTPUT_DIR}/persistent_homology.png", dpi=300)
        plt.close()
        print("  > Сохранено: persistent_homology.png (используя ripser)")

    except ImportError:
        # Заглушка
        plt.text(0.5, 0.5, 'Persistent Homology\nNot available on Windows\n(requires ripser/gudhi)',
                 ha='center', va='center', fontsize=14)
        plt.title('Persistent Homology (Not Available)', fontsize=16)
        plt.axis('off')
        plt.savefig(f"{OUTPUT_DIR}/persistent_homology.png", dpi=300)
        plt.close()
        print("  > Сохранено: persistent_homology.png (заглушка)")

def granger_causality_test(df, messages_df):
    """
    Проведение Granger Causality test для временных рядов сообщений.
    Проверяет, влияет ли активность в одном канале на активность в другом.
    """
    print("\nПроведение Granger Causality test...")

    try:
        from statsmodels.tsa.stattools import grangercausalitytests
        import statsmodels.api as sm

        # Подготовка временных рядов: агрегируем количество сообщений по часам для каждого канала
        if messages_df is None:
            print("  > Нет данных сообщений, пропускаем Granger test")
            return

        # Преобразуем timestamp в datetime
        messages_df['datetime'] = pd.to_datetime(messages_df['timestamp'], unit='s')
        messages_df['hour'] = messages_df['datetime'].dt.floor('h')

        # Создаем временной ряд: количество постов в час для каждого канала
        hourly_counts = messages_df.groupby(['hour', 'channel']).size().reset_index(name='count')
        pivot_df = hourly_counts.pivot(index='hour', columns='channel', values='count').fillna(0)

        # Выбираем два самых активных канала для теста
        top_channels = pivot_df.sum().nlargest(2).index.tolist()
        if len(top_channels) < 2:
            print("  > Недостаточно каналов для теста")
            return

        # Подготовка данных для теста
        data_for_test = pivot_df[top_channels].values

        # Granger test
        max_lag = 5
        print(f"  > Тестируем каналы: {top_channels[0]} -> {top_channels[1]}")
        test_result = grangercausalitytests(data_for_test, maxlag=max_lag, verbose=False)

        # Извлечение p-value для лага 1
        p_values = []
        for lag in range(1, max_lag + 1):
            p_val = test_result[lag][0]['ssr_ftest'][1]
            p_values.append(p_val)

        # Визуализация
        plt.figure(figsize=(10, 6))
        plt.bar(range(1, max_lag + 1), p_values, color='steelblue')
        plt.axhline(y=0.05, color='red', linestyle='--', label='alpha=0.05')
        plt.xlabel('Lag')
        plt.ylabel('p-value')
        plt.title(f'Granger Causality: {top_channels[0]} -> {top_channels[1]}', fontsize=14)
        plt.xticks(range(1, max_lag + 1))
        plt.legend()
        plt.tight_layout()
        plt.savefig(f"{OUTPUT_DIR}/granger_causality.png", dpi=300)
        plt.close()
        print(f"  > Сохранено: granger_causality.png")
        print(f"  > P-values for lags 1-{max_lag}: {[f'{p:.4f}' for p in p_values]}")

    except ImportError:
        print("  > statsmodels не установлен. Установите: pip install statsmodels")
        # Создаем заглушку
        plt.figure(figsize=(10, 6))
        plt.text(0.5, 0.5, 'Granger Causality Test\nRequires statsmodels\ntry: pip install statsmodels',
                 ha='center', va='center', fontsize=14)
        plt.title('Granger Causality (Not Available)', fontsize=16)
        plt.axis('off')
        plt.savefig(f"{OUTPUT_DIR}/granger_causality.png", dpi=300)
        plt.close()
        print("  > Сохранено: granger_causality.png (заглушка)")
    except Exception as e:
        print(f"  > Ошибка при проведении теста: {e}")

def main():
    print("="*50)
    print("НАЧАЛО EDA АНАЛИЗА")
    print("="*50)

    # Загрузка данных пар
    df = load_data()

    # Загрузка данных сообщений для Granger test
    messages_df = None
    try:
        messages_path = os.path.join(DATA_DIR, "messages_with_features.parquet")
        messages_df = pd.read_parquet(messages_path)
        print(f"Загружено {len(messages_df)} сообщений для Granger test")
    except Exception as e:
        print(f"Не удалось загрузить сообщения: {e}")

    # 1. Распределение классов
    plot_label_distribution(df)

    # 2. Корреляция
    plot_correlation_heatmap(df)

    # 3. PCA 3D
    plot_pca_3d(df)

    # 4. t-SNE
    plot_tsne(df)

    # 5. UMAP
    plot_umap(df)

    # 6. SHAP
    plot_shap_summary(df)

    # 7. Persistent Homology
    plot_persistent_homology()

    # 8. Granger Causality Test
    granger_causality_test(df, messages_df)

    print("\n" + "="*50)
    print("EDA АНАЛИЗ ЗАВЕРШЕН")
    print(f"Все графики сохранены в папку: {OUTPUT_DIR}/")
    print("="*50)

if __name__ == "__main__":
    main()