code/temporal_significance_test.py

# -*- coding: utf-8 -*-
"""

=============================================================================

KIỂM ĐỊNH THỐNG KÊ Ý NGHĨA VỚI MOVING-BLOCK BOOTSTRAP (M-BB)

BẢN CẢI TIẾN: HỖ TRỢ ĐA KHỐI (48 BƯỚC / 24H & 336 BƯỚC / 7 NGÀY)

VÀ XỬ LÝ CHẶT CHẼ RANH GIỚI CÁC PHÂN ĐOẠN LIÊN TỤC (CONTIGUOUS SEGMENTS)

=============================================================================

Đặc điểm:

1. Xác định 16 phân đoạn liên tục trong tập đánh giá 14,385 mẫu.

2. Khi lấy mẫu khối (Block Sampling), chỉ lấy các block nằm trọn trong cùng một

   phân đoạn liên tục (không lấy block bắc cầu qua khoảng đứt gãy thời gian).

3. Đánh giá đa quy mô:

   - Block = 48 bước (1 ngày - chu kỳ phụ tải ngày).

   - Block = 336 bước (7 ngày - chu kỳ tuần).

4. Tính toán đầy đủ:

   - ΔRMSE và ΔMAE

   - Khoảng tin cậy 95% CI

   - Tỷ lệ mô hình cây tốt hơn Persistence

=============================================================================

"""

import os
import sys
import json
from pathlib import Path

if sys.platform.startswith("win"):
    try:
        sys.stdout.reconfigure(encoding="utf-8")
        sys.stderr.reconfigure(encoding="utf-8")
    except Exception:
        pass

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.metrics import mean_absolute_error, mean_squared_error
from sklearn.impute import SimpleImputer
import xgboost as xgb
from sklearn.ensemble import HistGradientBoostingRegressor

sns.set_theme(style="whitegrid", palette="muted")
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False


def build_contiguous_blocks(segment_ids, block_length):
    """

    Xác định tất cả các điểm bắt đầu hợp lệ sao cho toàn bộ block dài block_length

    nằm trọn trong cùng một phân đoạn thời gian liên tục.

    """
    valid_starts = []
    n = len(segment_ids)
    for i in range(n - block_length + 1):
        if segment_ids[i] == segment_ids[i + block_length - 1]:
            valid_starts.append(i)
    return np.array(valid_starts)


def contiguous_moving_block_bootstrap(

    y_true, y_pred_cand, y_pred_base, segment_ids,

    block_length=48, n_boot=1000, seed=42

):
    np.random.seed(seed)
    n = len(y_true)
    valid_starts = build_contiguous_blocks(segment_ids, block_length)

    if len(valid_starts) == 0:
        raise ValueError(f"Không tìm thấy phân đoạn liên tục nào đủ dài {block_length} bước.")

    num_blocks = int(np.ceil(n / block_length))

    diff_rmse_list = []
    diff_mae_list = []

    rmse_base_orig = np.sqrt(mean_squared_error(y_true, y_pred_base))
    rmse_cand_orig = np.sqrt(mean_squared_error(y_true, y_pred_cand))
    mae_base_orig = mean_absolute_error(y_true, y_pred_base)
    mae_cand_orig = mean_absolute_error(y_true, y_pred_cand)

    for b in range(n_boot):
        chosen_starts = np.random.choice(valid_starts, size=num_blocks, replace=True)
        boot_indices = np.concatenate([np.arange(s, s + block_length) for s in chosen_starts])[:n]

        yt = y_true[boot_indices]
        yp_c = y_pred_cand[boot_indices]
        yp_b = y_pred_base[boot_indices]

        rmse_c = np.sqrt(np.mean((yt - yp_c) ** 2))
        rmse_b = np.sqrt(np.mean((yt - yp_b) ** 2))
        diff_rmse = rmse_b - rmse_c  # > 0 nghĩa là Candidate tốt hơn (giảm RMSE)

        mae_c = np.mean(np.abs(yt - yp_c))
        mae_b = np.mean(np.abs(yt - yp_b))
        diff_mae = mae_b - mae_c

        diff_rmse_list.append(diff_rmse)
        diff_mae_list.append(diff_mae)

    diff_rmse_arr = np.array(diff_rmse_list)
    diff_mae_arr = np.array(diff_mae_list)

    ci_rmse_low, ci_rmse_high = np.percentile(diff_rmse_arr, [2.5, 97.5])
    ci_mae_low, ci_mae_high = np.percentile(diff_mae_arr, [2.5, 97.5])

    return {
        "block_length": block_length,
        "valid_blocks_count": len(valid_starts),
        "rmse_base": rmse_base_orig,
        "rmse_cand": rmse_cand_orig,
        "delta_rmse_mean": np.mean(diff_rmse_arr),
        "ci_rmse_95": (ci_rmse_low, ci_rmse_high),
        "rmse_better_pct": np.mean(diff_rmse_arr > 0) * 100,
        "mae_base": mae_base_orig,
        "mae_cand": mae_cand_orig,
        "delta_mae_mean": np.mean(diff_mae_arr),
        "ci_mae_95": (ci_mae_low, ci_mae_high),
        "mae_better_pct": np.mean(diff_mae_arr > 0) * 100,
        "diff_rmse_samples": diff_rmse_arr,
        "diff_mae_samples": diff_mae_arr
    }


def main():
    base_dir = Path("output_preprocessing")
    data_dir = base_dir / "data"
    reports_dir = base_dir / "reports"
    figures_dir = base_dir / "figures"

    print("=" * 80)
    print("CHẠY KIỂM ĐỊNH CONTIGUOUS MOVING-BLOCK BOOTSTRAP (BLOCK 48 & BLOCK 336)")
    print("=" * 80)

    # 1. Đọc dữ liệu
    test_df = pd.read_csv(data_dir / "test.csv")
    train_df = pd.read_csv(data_dir / "train.csv")
    val_df = pd.read_csv(data_dir / "validation.csv")

    test_df["TIME"] = pd.to_datetime(test_df["TIME"])
    test_df = test_df.sort_values("TIME").reset_index(drop=True)

    candidate_features = [
        "UA", "UB", "UC", "IA", "IB", "IC", "CosF", "P_kW", "Q_kVAr",
        "U_avg", "U_max", "U_min", "U_unbalance_pct",
        "I_avg", "I_max", "I_min", "I_unbalance_pct",
        "S_kVA", "lag_1", "lag_2", "lag_3", "lag_48", "lag_336",
        "rolling_mean_3", "rolling_mean_6", "rolling_mean_48",
        "hour", "dayofweek", "dayofmonth", "month", "is_weekend",
        "sin_hour", "cos_hour"
    ]
    feature_cols = [c for c in candidate_features if c in train_df.columns]
    target_col = "target_P_kW"

    # Lọc mask chung 14,385 mẫu
    pred_persistence = test_df["P_kW"].values
    pred_s48 = test_df["lag_48"].values
    pred_s336 = test_df["lag_336"].values
    pred_ma3 = test_df["rolling_mean_3"].values
    y_test = test_df[target_col].values

    mask_common = (
        np.isfinite(y_test) &
        np.isfinite(pred_persistence) &
        np.isfinite(pred_s48) &
        np.isfinite(pred_s336) &
        np.isfinite(pred_ma3)
    )

    eval_df = test_df[mask_common].copy().reset_index(drop=True)
    y_eval = y_test[mask_common]
    p_eval = pred_persistence[mask_common]

    # Tính toán chính xác các phân đoạn thời gian liên tục (Contiguous Segments)
    diffs = eval_df["TIME"].diff()
    segment_ids = (diffs != pd.Timedelta(minutes=30)).cumsum().values
    n_segments = len(np.unique(segment_ids))
    print(f"-> Tập đánh giá chuẩn: {len(y_eval):,} mẫu, gồm {n_segments} phân đoạn thời gian liên tục.")

    # Huấn luyện XGBoost & HistGBDT
    imputer = SimpleImputer(strategy="median")
    X_train_imp = imputer.fit_transform(train_df[feature_cols])
    y_train = train_df[target_col].values
    X_val_imp = imputer.transform(val_df[feature_cols])
    y_val = val_df[target_col].values
    X_test_imp = imputer.transform(test_df[feature_cols])

    print("\nHuấn luyện XGBoost Regressor...")
    model_xgb = xgb.XGBRegressor(
        n_estimators=300, learning_rate=0.05, max_depth=6, subsample=0.8,
        colsample_bytree=0.8, random_state=42, n_jobs=-1, early_stopping_rounds=30
    )
    model_xgb.fit(X_train_imp, y_train, eval_set=[(X_val_imp, y_val)], verbose=False)
    pred_xgb = model_xgb.predict(X_test_imp)[mask_common]

    print("Huấn luyện HistGradientBoosting Regressor...")
    model_hgb = HistGradientBoostingRegressor(
        max_iter=300, learning_rate=0.05, max_depth=6, random_state=42, early_stopping=True
    )
    model_hgb.fit(X_train_imp, y_train)
    pred_hgb = model_hgb.predict(X_test_imp)[mask_common]

    # Kiểm định với Block 48 (24h) và Block 336 (7 ngày)
    configs = [
        {"name": "24h (48 bước)", "block": 48},
        {"name": "7 ngày (336 bước)", "block": 336}
    ]

    all_results = {}

    for cfg in configs:
        b_len = cfg["block"]
        b_name = cfg["name"]
        print(f"\n" + "-" * 70)
        print(f"TIẾN HÀNH CONTIGUOUS M-BB CHO KHUNG {b_name}...")
        print("-" * 70)

        res_xgb = contiguous_moving_block_bootstrap(
            y_eval, pred_xgb, p_eval, segment_ids, block_length=b_len, n_boot=1000
        )
        res_hgb = contiguous_moving_block_bootstrap(
            y_eval, pred_hgb, p_eval, segment_ids, block_length=b_len, n_boot=1000
        )

        print(f"XGBoost vs Persistence (Block {b_len}):")
        print(f"  ΔRMSE trung bình: {res_xgb['delta_rmse_mean']:+.2f} kW | 95% CI: [{res_xgb['ci_rmse_95'][0]:.2f}; {res_xgb['ci_rmse_95'][1]:.2f}] kW")
        print(f"  Tỷ lệ thắng RMSE: {res_xgb['rmse_better_pct']:.1f}%")
        print(f"  ΔMAE trung bình:  {res_xgb['delta_mae_mean']:+.2f} kW | 95% CI: [{res_xgb['ci_mae_95'][0]:.2f}; {res_xgb['ci_mae_95'][1]:.2f}] kW")

        print(f"HistGBDT vs Persistence (Block {b_len}):")
        print(f"  ΔRMSE trung bình: {res_hgb['delta_rmse_mean']:+.2f} kW | 95% CI: [{res_hgb['ci_rmse_95'][0]:.2f}; {res_hgb['ci_rmse_95'][1]:.2f}] kW")
        print(f"  Tỷ lệ thắng RMSE: {res_hgb['rmse_better_pct']:.1f}%")
        print(f"  ΔMAE trung bình:  {res_hgb['delta_mae_mean']:+.2f} kW | 95% CI: [{res_hgb['ci_mae_95'][0]:.2f}; {res_hgb['ci_mae_95'][1]:.2f}] kW")

        all_results[f"block_{b_len}"] = {
            "block_length_steps": b_len,
            "block_duration": b_name,
            "valid_blocks_count": res_xgb["valid_blocks_count"],
            "XGBoost_vs_Persistence": {
                "delta_rmse_mean_kW": round(float(res_xgb['delta_rmse_mean']), 2),
                "ci_rmse_95_kW": [round(float(x), 2) for x in res_xgb['ci_rmse_95']],
                "rmse_better_pct": round(float(res_xgb['rmse_better_pct']), 1),
                "delta_mae_mean_kW": round(float(res_xgb['delta_mae_mean']), 2),
                "ci_mae_95_kW": [round(float(x), 2) for x in res_xgb['ci_mae_95']],
                "is_significant_at_005": bool(res_xgb['ci_rmse_95'][0] > 0)
            },
            "HistGBDT_vs_Persistence": {
                "delta_rmse_mean_kW": round(float(res_hgb['delta_rmse_mean']), 2),
                "ci_rmse_95_kW": [round(float(x), 2) for x in res_hgb['ci_rmse_95']],
                "rmse_better_pct": round(float(res_hgb['rmse_better_pct']), 1),
                "delta_mae_mean_kW": round(float(res_hgb['delta_mae_mean']), 2),
                "ci_mae_95_kW": [round(float(x), 2) for x in res_hgb['ci_mae_95']],
                "is_significant_at_005": bool(res_hgb['ci_rmse_95'][0] > 0)
            }
        }

    # Lưu JSON báo cáo tổng hợp
    full_report = {
        "method": "Contiguous Moving-Block Bootstrap (Zero Cross-Gap Leakage, B=1000)",
        "evaluated_samples": len(y_eval),
        "contiguous_segments_count": n_segments,
        "results": all_results
    }
    with open(reports_dir / "moving_block_bootstrap_summary.json", "w", encoding="utf-8") as f:
        json.dump(full_report, f, indent=4, ensure_ascii=False)

    # Vẽ biểu đồ so sánh 2x2: (Block 48 & Block 336) cho XGBoost và HistGBDT
    fig, axes = plt.subplots(2, 2, figsize=(15, 10))

    # Hàng 1: Block 48
    r48_xgb = contiguous_moving_block_bootstrap(y_eval, pred_xgb, p_eval, segment_ids, block_length=48, n_boot=1000)
    r48_hgb = contiguous_moving_block_bootstrap(y_eval, pred_hgb, p_eval, segment_ids, block_length=48, n_boot=1000)

    sns.histplot(r48_xgb["diff_rmse_samples"], bins=40, kde=True, ax=axes[0, 0], color="#ff7f0e", alpha=0.6)
    axes[0, 0].axvline(0, color="black", ls="--", lw=1.5, label="Không chênh lệch (0 kW)")
    axes[0, 0].axvline(r48_xgb['ci_rmse_95'][0], color="red", ls=":", lw=1.5, label=f"2.5%: {r48_xgb['ci_rmse_95'][0]:.1f}")
    axes[0, 0].axvline(r48_xgb['ci_rmse_95'][1], color="red", ls=":", lw=1.5, label=f"97.5%: {r48_xgb['ci_rmse_95'][1]:.1f}")
    axes[0, 0].set_title(f"A. XGBoost vs Persistence [Block 48 (24h)]\nΔRMSE: {r48_xgb['delta_rmse_mean']:+.1f} kW | Thắng: {r48_xgb['rmse_better_pct']:.1f}%", fontweight="bold")
    axes[0, 0].set_xlabel("ΔRMSE (Persistence - XGBoost) [kW]")
    axes[0, 0].legend()

    sns.histplot(r48_hgb["diff_rmse_samples"], bins=40, kde=True, ax=axes[0, 1], color="#2ca02c", alpha=0.6)
    axes[0, 1].axvline(0, color="black", ls="--", lw=1.5, label="Không chênh lệch (0 kW)")
    axes[0, 1].axvline(r48_hgb['ci_rmse_95'][0], color="red", ls=":", lw=1.5, label=f"2.5%: {r48_hgb['ci_rmse_95'][0]:.1f}")
    axes[0, 1].axvline(r48_hgb['ci_rmse_95'][1], color="red", ls=":", lw=1.5, label=f"97.5%: {r48_hgb['ci_rmse_95'][1]:.1f}")
    axes[0, 1].set_title(f"B. HistGBDT vs Persistence [Block 48 (24h)]\nΔRMSE: {r48_hgb['delta_rmse_mean']:+.1f} kW | Thắng: {r48_hgb['rmse_better_pct']:.1f}%", fontweight="bold")
    axes[0, 1].set_xlabel("ΔRMSE (Persistence - HistGBDT) [kW]")
    axes[0, 1].legend()

    # Hàng 2: Block 336
    r336_xgb = contiguous_moving_block_bootstrap(y_eval, pred_xgb, p_eval, segment_ids, block_length=336, n_boot=1000)
    r336_hgb = contiguous_moving_block_bootstrap(y_eval, pred_hgb, p_eval, segment_ids, block_length=336, n_boot=1000)

    sns.histplot(r336_xgb["diff_rmse_samples"], bins=40, kde=True, ax=axes[1, 0], color="#ff7f0e", alpha=0.6)
    axes[1, 0].axvline(0, color="black", ls="--", lw=1.5, label="Không chênh lệch (0 kW)")
    axes[1, 0].axvline(r336_xgb['ci_rmse_95'][0], color="red", ls=":", lw=1.5, label=f"2.5%: {r336_xgb['ci_rmse_95'][0]:.1f}")
    axes[1, 0].axvline(r336_xgb['ci_rmse_95'][1], color="red", ls=":", lw=1.5, label=f"97.5%: {r336_xgb['ci_rmse_95'][1]:.1f}")
    axes[1, 0].set_title(f"C. XGBoost vs Persistence [Block 336 (7 ngày)]\nΔRMSE: {r336_xgb['delta_rmse_mean']:+.1f} kW | Thắng: {r336_xgb['rmse_better_pct']:.1f}%", fontweight="bold")
    axes[1, 0].set_xlabel("ΔRMSE (Persistence - XGBoost) [kW]")
    axes[1, 0].legend()

    sns.histplot(r336_hgb["diff_rmse_samples"], bins=40, kde=True, ax=axes[1, 1], color="#2ca02c", alpha=0.6)
    axes[1, 1].axvline(0, color="black", ls="--", lw=1.5, label="Không chênh lệch (0 kW)")
    axes[1, 1].axvline(r336_hgb['ci_rmse_95'][0], color="red", ls=":", lw=1.5, label=f"2.5%: {r336_hgb['ci_rmse_95'][0]:.1f}")
    axes[1, 1].axvline(r336_hgb['ci_rmse_95'][1], color="red", ls=":", lw=1.5, label=f"97.5%: {r336_hgb['ci_rmse_95'][1]:.1f}")
    axes[1, 1].set_title(f"D. HistGBDT vs Persistence [Block 336 (7 ngày)]\nΔRMSE: {r336_hgb['delta_rmse_mean']:+.1f} kW | Thắng: {r336_hgb['rmse_better_pct']:.1f}%", fontweight="bold")
    axes[1, 1].set_xlabel("ΔRMSE (Persistence - HistGBDT) [kW]")
    axes[1, 1].legend()

    fig.suptitle("13. Kiểm định ý nghĩa thống kê Đa chu kỳ (Block 24h & Block 7 ngày) qua Contiguous Moving-Block Bootstrap", fontsize=13, fontweight="bold")
    fig.tight_layout()
    fig.savefig(figures_dir / "13_moving_block_bootstrap_rmse.png", dpi=200)
    plt.close(fig)

    print("\n" + "=" * 80)
    print("✅ ĐÃ HOÀN THÀNH TOÀN DIỆN CONTIGUOUS MOVING-BLOCK BOOTSTRAP!")
    print(f"📊 Báo cáo JSON: {reports_dir / 'moving_block_bootstrap_summary.json'}")
    print(f"📈 Biểu đồ 2x2: {figures_dir / '13_moving_block_bootstrap_rmse.png'}")
    print("=" * 80)


if __name__ == "__main__":
    main()