code/power_load_preprocessing_pipeline.py

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

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

PIPELINE TIỀN XỬ LÝ DỮ LIỆU DỰ BÁO PHỤ TẢI ĐIỆN 30 PHÚT (BẢN CHUẨN XÁC THỜI GIAN)

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

Khắc phục triệt để 3 vấn đề lớn:

1. Loại bỏ hoàn toàn mốc thời gian bất thường sau 30/08/2025 (khống chế chặt dải 2020 - 2025).

2. Tách bạch rõ 3 mức độ thiếu:

   - Thiếu trong dữ liệu thô (raw missing).

   - Thiếu do tạo lưới 30 phút (grid gaps).

   - Thiếu còn lại sau khi nội suy có kiểm soát (<= 2 bước).

3. Chia Train/Validation/Test TRƯỚC KHI NỘI SUY (Chỉ nội suy độc lập trong từng tập,

   triệt tiêu hoàn toàn nguy cơ rò rỉ dữ liệu qua ranh giới tập dữ liệu - Zero Leakage).

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

"""

import os
import sys
import re
import json
import math
import argparse
from pathlib import Path
from datetime import datetime

# Đảm bảo in tiếng Việt chuẩn trên Windows console
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

# Cấu hình thẩm mỹ đồ thị
sns.set_theme(style="whitegrid", palette="muted")
plt.rcParams['font.sans-serif'] = 'DejaVu Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.titlesize'] = 14
plt.rcParams['axes.labelsize'] = 11

NUMERIC_COLS = [
    "UA", "UB", "UC", "IA", "IB", "IC", "CosF",
    "P_kW", "Q_kVAr", "P_MW", "Q_MVAr",
    "P_MW_raw", "P_kW_raw", "Q_MVAr_raw", "Q_kVAr_raw"
]
MIN_VALID_DATE = pd.to_datetime("2020-01-01 00:00:00")
MAX_VALID_DATE = pd.to_datetime("2025-08-31 23:59:59")


# =============================================================================
# 1. TÌM KIẾM VÀ ĐỌC FILE CSV HỢP LỆ
# =============================================================================
def find_csv_files(data_dir: Path) -> list[Path]:
    if not data_dir.exists():
        raise FileNotFoundError(f"Thư mục không tồn tại: {data_dir}")

    all_csvs = sorted(data_dir.rglob("*.csv"))
    exclude_keywords = ["report", "summary", "cleaned", "featured", "supervised", "train", "validation", "test", "merged_with_features"]
    valid_csvs = [p for p in all_csvs if not any(kw in p.name.lower() for kw in exclude_keywords)]
    return valid_csvs


# =============================================================================
# 2. CHUẨN HÓA THỜI GIAN VÀ DỮ LIỆU SỐ (CHẶT CHẼ, CHỐNG NHẦM NĂM)
# =============================================================================
def parse_time_robust(df: pd.DataFrame) -> pd.Series:
    """

    Parse thời gian chặt chẽ:

    - Bắt buộc chuỗi thời gian phải chứa năm 2020-2025.

    - Tuyệt đối không dùng pd.to_datetime tự động điền ngày hiện tại (2026) khi chuỗi chỉ có mm:ss.

    """
    time_candidates = ["TIME", "Time", "ThoiGian", "THOIGIAN", "Gio", "GIO"]
    time_col = next((c for c in time_candidates if c in df.columns), None)
    
    if not time_col:
        return pd.Series([pd.NaT] * len(df))

    s_time = df[time_col].astype(str).str.strip()

    # Kiểm tra xem chuỗi có chứa năm 4 chữ số (2020-2025) hay không
    has_year = s_time.str.contains(r"\b202[0-5]\b", regex=True)

    # Nếu chuỗi không có năm, kiểm tra xem có cột DATE/Ngay đi kèm không
    date_candidates = ["DATE", "Date", "Ngay", "NGAY", "Ngày"]
    date_col = next((c for c in date_candidates if c in df.columns), None)

    if date_col:
        s_date = pd.to_datetime(df[date_col], errors="coerce", format="mixed")
        s_time_fixed = s_time.str.replace(r"\s+", "", regex=True)
        s_time_fixed = s_time_fixed.where(s_time_fixed.str.count(":") == 2, "00:" + s_time_fixed)
        dt = s_date + pd.to_timedelta(s_time_fixed, errors="coerce")
        return dt

    # Nếu không có cột DATE và s_time cũng không chứa năm 202x -> coi như NaT
    # (Tránh trường hợp pandas tự điền năm 2026 hiện tại)
    parsed = pd.Series([pd.NaT] * len(df), index=df.index)
    if has_year.any():
        parsed.loc[has_year] = pd.to_datetime(s_time.loc[has_year], errors="coerce", format="mixed")

    return parsed


def clean_single_numeric_series(series: pd.Series) -> pd.Series:
    """

    Xử lý số thập phân tiếng Việt (dấu phẩy), loại bỏ dấu chấm phân tách hàng nghìn nếu có.

    """
    s = series.astype(str).str.strip().str.replace(",", ".", regex=False)
    # Nếu có nhiều hơn 1 dấu chấm (ví dụ: '79.760.600') -> định dạng ngăn cách nghìn
    s = s.apply(lambda x: ".".join(x.split(".")[:2]) if x.count(".") > 1 else x)
    return pd.to_numeric(s, errors="coerce")


def load_and_standardize_single_file(file_path: Path) -> pd.DataFrame:
    try:
        df = pd.read_csv(file_path, low_memory=False, encoding="utf-8")
    except UnicodeDecodeError:
        df = pd.read_csv(file_path, low_memory=False, encoding="latin1")

    df.columns = [str(c).strip() for c in df.columns]

    # Parse TIME chặt chẽ
    df["TIME"] = parse_time_robust(df)
    df = df.dropna(subset=["TIME"]).copy()

    # Khống chế thời gian thực tế: Chỉ nhận dữ liệu từ 2020-01-01 đến 2025-08-31
    df = df[(df["TIME"] >= MIN_VALID_DATE) & (df["TIME"] <= MAX_VALID_DATE)].copy()
    if len(df) == 0:
        return pd.DataFrame()

    raw_num_cols = ["UA", "UB", "UC", "IA", "IB", "IC", "CosF", "P_kW", "Q_kVAr"]
    for c in raw_num_cols:
        if c in df.columns:
            df[c] = clean_single_numeric_series(df[c])

    # CHUẨN HÓA ĐƠN VỊ VẬT LÝ VÀ PHÂN ĐỊNH NGOẠI LAI (PHYSICAL CONTROL):
    # - Lưu trữ nguyên bản giá trị thô: P_MW_raw, P_kW_raw
    # - Đánh dấu P_outlier = 1 nếu công suất vượt ngưỡng kiểm soát dữ liệu (20 MW)
    #   hoặc không tương thích với điện áp và dòng điện đo được đồng thời
    # - Giá trị bất thường chuyển sang NaN để xử lý bằng nội suy có kiểm soát (<= 2 bước).
    if "P_kW" in df.columns:
        df["P_MW_raw"] = df["P_kW"].copy()
        df["P_kW_raw"] = df["P_MW_raw"] * 1000.0

        # Kiểm tra ngưỡng kiểm soát dữ liệu 20 MW (trung bình lộ chỉ 8 MW, trần tải dòng <= 12 MW)
        is_outlier = df["P_MW_raw"] > 20.0
        df["P_outlier"] = is_outlier.astype(int)

        # Chuyển P_MW và P_kW của bản ghi ngoại lai thành NaN để nội suy an toàn (hoặc loại bỏ nếu kéo dài)
        df["P_MW"] = df["P_MW_raw"].where(~is_outlier, np.nan)
        df["P_kW"] = df["P_MW"] * 1000.0
    else:
        df["P_outlier"] = 0

    if "Q_kVAr" in df.columns:
        df["Q_MVAr_raw"] = df["Q_kVAr"].copy()
        df["Q_kVAr_raw"] = df["Q_MVAr_raw"] * 1000.0
        # Nếu P bị đánh dấu ngoại lai, Q cũng chuyển sang NaN
        df["Q_MVAr"] = df["Q_MVAr_raw"].where(~df["P_outlier"].astype(bool), np.nan)
        df["Q_kVAr"] = df["Q_MVAr"] * 1000.0

    if "MaDiemDo" in df.columns:
        df["MaDiemDo"] = pd.to_numeric(df["MaDiemDo"], errors="coerce")
        df = df.dropna(subset=["MaDiemDo"]).copy()
        df["MaDiemDo"] = df["MaDiemDo"].astype(int)

    if "TenDiemDo" in df.columns:
        df["TenDiemDo"] = df["TenDiemDo"].astype(str).str.strip()

    df["source_file"] = file_path.name
    return df


# =============================================================================
# 3. KIỂM TRA CHẤT LƯỢNG DỮ LIỆU THÔ (RAW QUALITY REPORT)
# =============================================================================
def audit_data_quality(df: pd.DataFrame, reports_dir: Path) -> dict:
    total_rows = len(df)
    unique_points = df["MaDiemDo"].nunique() if "MaDiemDo" in df.columns else 0
    t_min = df["TIME"].min()
    t_max = df["TIME"].max()

    dup_count = int(df.duplicated(subset=["TIME", "MaDiemDo"]).sum())
    missing_stats = df.isna().sum()
    missing_pct = (missing_stats / total_rows * 100).round(2)

    report_data = []
    for col in df.columns:
        report_data.append({
            "column_name": col,
            "data_type": str(df[col].dtype),
            "missing_count": int(missing_stats[col]),
            "missing_pct": float(missing_pct[col]),
            "unique_values": int(df[col].nunique())
        })
    df_report = pd.DataFrame(report_data)
    df_report.to_csv(reports_dir / "quality_report.csv", index=False, encoding="utf-8-sig")

    return {
        "total_rows": total_rows,
        "unique_points": unique_points,
        "time_min": str(t_min),
        "time_max": str(t_max),
        "duplicate_rows": dup_count,
        "missing_pct_P_kW": float(missing_pct.get("P_kW", 0.0))
    }


# =============================================================================
# 4. TẠO LƯỚI THỜI GIAN 30 PHÚT CHUẨN (KHÔNG NỘI SUY MÙ)
# =============================================================================
def build_regular_30m_grid(df_raw: pd.DataFrame) -> pd.DataFrame:
    """

    Quy chuẩn dữ liệu về lưới 30 phút đều đặn:

    - Làm tròn TIME về mốc 30 phút gần nhất (00:00, 00:30, ...)

    - Khử trùng lặp trên cặp (MaDiemDo, TIME_30m)

    - Reindex đầy đủ từ t_min đến t_max theo chu kỳ 30 phút (CHƯA NỘI SUY).

    """
    df = df_raw.copy()
    df["TIME_30m"] = df["TIME"].dt.round("30min")
    df_sorted = df.sort_values(by=["MaDiemDo", "TIME_30m", "TIME"]).reset_index(drop=True)
    df_dedup = df_sorted.drop_duplicates(subset=["MaDiemDo", "TIME_30m"], keep="first").copy()

    # Dải thời gian chung từ min đến max thực tế
    t_start = df_dedup["TIME_30m"].min()
    t_end = df_dedup["TIME_30m"].max()
    full_timeline = pd.date_range(start=t_start, end=t_end, freq="30min", name="TIME")

    grid_list = []
    for point_id, group in df_dedup.groupby("MaDiemDo"):
        if "TIME" in group.columns:
            group = group.drop(columns=["TIME"])
        group = group.set_index("TIME_30m").sort_index()

        group_reindexed = group.reindex(full_timeline)
        group_reindexed["MaDiemDo"] = point_id
        if "TenDiemDo" in group.columns:
            group_reindexed["TenDiemDo"] = group["TenDiemDo"].dropna().iloc[0] if len(group["TenDiemDo"].dropna()) > 0 else f"Point_{point_id}"
        if "MBA" in group.columns:
            group_reindexed["MBA"] = group["MBA"].dropna().iloc[0] if len(group["MBA"].dropna()) > 0 else "Unknown"

        # Đánh dấu dữ liệu thiếu trên lưới trước khi can thiệp
        group_reindexed["was_raw_missing"] = group_reindexed["P_kW"].isna().astype(int)
        if "P_outlier" in group_reindexed.columns:
            group_reindexed["P_outlier"] = group_reindexed["P_outlier"].fillna(0).astype(int)
        grid_list.append(group_reindexed.reset_index())

    df_grid = pd.concat(grid_list, ignore_index=True)
    return df_grid


# =============================================================================
# 5. CHIA TẬP TRAIN / VAL / TEST TRƯỚC KHI NỘI SUY (CHỐNG RÒ RỈ DỮ LIỆU)
# =============================================================================
def split_timeline_into_subsets(df_grid: pd.DataFrame, train_ratio=0.70, val_ratio=0.15):
    """

    Cắt tập dữ liệu theo mốc thời gian tuyệt đối trước khi bất kỳ thao tác nội suy nào diễn ra.

    """
    timeline_unique = pd.Series(df_grid["TIME"].unique()).sort_values().reset_index(drop=True)
    total_slots = len(timeline_unique)

    train_end_idx = int(total_slots * train_ratio)
    val_end_idx = int(total_slots * (train_ratio + val_ratio))

    t_train_cutoff = timeline_unique.iloc[train_end_idx - 1]
    t_val_start = timeline_unique.iloc[train_end_idx]
    t_val_cutoff = timeline_unique.iloc[val_end_idx - 1]
    t_test_start = timeline_unique.iloc[val_end_idx]

    train_mask = df_grid["TIME"] <= t_train_cutoff
    val_mask = (df_grid["TIME"] >= t_val_start) & (df_grid["TIME"] <= t_val_cutoff)
    test_mask = df_grid["TIME"] >= t_test_start

    train_subset = df_grid[train_mask].copy()
    val_subset = df_grid[val_mask].copy()
    test_subset = df_grid[test_mask].copy()

    split_meta = {
        "train_start": str(train_subset["TIME"].min()),
        "train_end": str(train_subset["TIME"].max()),
        "val_start": str(val_subset["TIME"].min()),
        "val_end": str(val_subset["TIME"].max()),
        "test_start": str(test_subset["TIME"].min()),
        "test_end": str(test_subset["TIME"].max()),
    }

    return train_subset, val_subset, test_subset, split_meta


# =============================================================================
# 6. NỘI SUY CÓ KIỂM SOÁT ĐỘC LẬP TRONG TỪNG TẬP (ISOLATED IMPUTATION)
# =============================================================================
def impute_isolated_subset(df_subset: pd.DataFrame, limit_steps: int = 2) -> pd.DataFrame:
    """

    Nội suy tuyến tính CHỈ trong phạm vi nội bộ của tập dữ liệu đó:

    - Không nhìn sang tương lai của tập kế tiếp.

    - Giới hạn tối đa limit_steps (2 bước = 1 tiếng).

    - Các khoảng trống dài hơn giữ nguyên NaN.

    """
    imputed_groups = []
    for point_id, group in df_subset.groupby("MaDiemDo"):
        group = group.set_index("TIME").sort_index().copy()
        present_nums = [c for c in NUMERIC_COLS if c in group.columns]

        was_missing = group["P_kW"].isna()
        # Nội suy có kiểm soát: method='time', limit=limit_steps, limit_direction='forward'
        imputed_vals = group[present_nums].interpolate(method="time", limit=limit_steps, limit_direction="forward")
        is_imputed = was_missing & imputed_vals["P_kW"].notna()

        group[present_nums] = imputed_vals
        group["is_imputed"] = is_imputed.astype(int)
        imputed_groups.append(group.reset_index())

    return pd.concat(imputed_groups, ignore_index=True)


# =============================================================================
# 7. TRÍCH XUẤT ĐẶC TRƯNG CHUỖI THỜI GIAN (FEATURE ENGINEERING)
# =============================================================================
def engineer_features(df_long: pd.DataFrame) -> pd.DataFrame:
    df = df_long.copy()
    eps = 1e-6

    # Đặc trưng 3 pha
    if set(["UA", "UB", "UC"]).issubset(df.columns):
        df["U_avg"] = df[["UA", "UB", "UC"]].mean(axis=1)
        df["U_max"] = df[["UA", "UB", "UC"]].max(axis=1)
        df["U_min"] = df[["UA", "UB", "UC"]].min(axis=1)
        df["U_unbalance_pct"] = (df["U_max"] - df["U_min"]) / (df["U_avg"].abs() + eps) * 100

    if set(["IA", "IB", "IC"]).issubset(df.columns):
        df["I_avg"] = df[["IA", "IB", "IC"]].mean(axis=1)
        df["I_max"] = df[["IA", "IB", "IC"]].max(axis=1)
        df["I_min"] = df[["IA", "IB", "IC"]].min(axis=1)
        df["I_unbalance_pct"] = (df["I_max"] - df["I_min"]) / (df["I_avg"].abs() + eps) * 100

    if "P_kW" in df.columns and "Q_kVAr" in df.columns:
        df["S_kVA"] = np.sqrt(df["P_kW"].fillna(0)**2 + df["Q_kVAr"].fillna(0)**2)
        df["S_MVA"] = df["S_kVA"] / 1000.0

    df["hour"] = df["TIME"].dt.hour + df["TIME"].dt.minute / 60.0
    df["dayofweek"] = df["TIME"].dt.dayofweek
    df["dayofmonth"] = df["TIME"].dt.day
    df["month"] = df["TIME"].dt.month
    df["is_weekend"] = (df["dayofweek"] >= 5).astype(int)
    df["sin_hour"] = np.sin(2 * np.pi * df["hour"] / 24.0)
    df["cos_hour"] = np.cos(2 * np.pi * df["hour"] / 24.0)

    # Lags & Rolling không rò rỉ dữ liệu (shift 1)
    feat_groups = []
    for point_id, grp in df.groupby("MaDiemDo"):
        grp = grp.sort_values("TIME").copy()
        for lag in [1, 2, 3, 48, 336]:
            grp[f"lag_{lag}"] = grp["P_kW"].shift(lag)
            if "P_MW" in grp.columns:
                grp[f"lag_mw_{lag}"] = grp["P_MW"].shift(lag)
        p_prev = grp["P_kW"].shift(1)
        grp["rolling_mean_3"] = p_prev.rolling(3, min_periods=1).mean()
        grp["rolling_mean_6"] = p_prev.rolling(6, min_periods=1).mean()
        grp["rolling_mean_48"] = p_prev.rolling(48, min_periods=1).mean()
        if "P_MW" in grp.columns:
            p_mw_prev = grp["P_MW"].shift(1)
            grp["rolling_mean_mw_3"] = p_mw_prev.rolling(3, min_periods=1).mean()
            grp["rolling_mean_mw_6"] = p_mw_prev.rolling(6, min_periods=1).mean()
            grp["rolling_mean_mw_48"] = p_mw_prev.rolling(48, min_periods=1).mean()
        feat_groups.append(grp)

    return pd.concat(feat_groups, ignore_index=True)


# =============================================================================
# 8. BẢNG TỔNG HỢP TRẠM VÀ PHÂN TÍCH RÕ NGUYÊN NHÂN THIẾU
# =============================================================================
def generate_station_summary(df_raw: pd.DataFrame, df_grid: pd.DataFrame, df_clean: pd.DataFrame, reports_dir: Path) -> pd.DataFrame:
    """

    Phân tách rạch ròi 3 mức độ thiếu:

    1. Thiếu trong file gốc (Raw Missing P_kW).

    2. Thiếu do tạo lưới 30m liên tục (Grid Gap Missing).

    3. Thiếu còn lại sau nội suy <= 2 bước (Remaining Missing).

    Đồng thời hiển thị công suất theo cả 2 đơn vị chuẩn: kW và MW.

    """
    summary = []
    total_slots = df_grid["TIME"].nunique()

    for pid in sorted(df_grid["MaDiemDo"].unique()):
        raw_sub = df_raw[df_raw["MaDiemDo"] == pid]
        grid_sub = df_grid[df_grid["MaDiemDo"] == pid]
        clean_sub = df_clean[df_clean["MaDiemDo"] == pid]

        name = grid_sub["TenDiemDo"].iloc[0] if "TenDiemDo" in grid_sub.columns else f"Point_{pid}"
        raw_rows = len(raw_sub)
        raw_nulls = raw_sub["P_kW"].isna().sum() if "P_kW" in raw_sub.columns else 0
        raw_null_pct = (raw_nulls / raw_rows * 100) if raw_rows > 0 else 0

        grid_missing = grid_sub["P_kW"].isna().sum()
        grid_missing_pct = (grid_missing / total_slots) * 100

        clean_missing = clean_sub["P_kW"].isna().sum()
        clean_missing_pct = (clean_missing / total_slots) * 100

        imputed_count = int(clean_sub["is_imputed"].sum())
        outlier_count = int(clean_sub["P_outlier"].sum()) if "P_outlier" in clean_sub.columns else 0
        p_mean_kw = clean_sub["P_kW"].mean()
        p_max_kw = clean_sub["P_kW"].max()
        p_mean_mw = clean_sub["P_MW"].mean() if "P_MW" in clean_sub.columns else (p_mean_kw / 1000.0 if pd.notna(p_mean_kw) else None)
        p_max_mw = clean_sub["P_MW"].max() if "P_MW" in clean_sub.columns else (p_max_kw / 1000.0 if pd.notna(p_max_kw) else None)
        p_raw_max_mw = raw_sub["P_MW_raw"].max() if "P_MW_raw" in raw_sub.columns else (raw_sub["P_kW"].max() if "P_kW" in raw_sub.columns else None)

        summary.append({
            "MaDiemDo": pid,
            "TenDiemDo": name,
            "raw_recorded_rows": raw_rows,
            "raw_null_pct": round(float(raw_null_pct), 2),
            "outliers_flagged_count": outlier_count,
            "total_30m_grid_slots": total_slots,
            "grid_missing_pct": round(float(grid_missing_pct), 2),
            "imputed_steps_count": imputed_count,
            "remaining_missing_pct": round(float(clean_missing_pct), 2),
            "P_mean_kW": round(float(p_mean_kw), 1) if pd.notna(p_mean_kw) else None,
            "P_max_kW": round(float(p_max_kw), 1) if pd.notna(p_max_kw) else None,
            "P_mean_MW": round(float(p_mean_mw), 3) if pd.notna(p_mean_mw) else None,
            "P_max_MW": round(float(p_max_mw), 3) if pd.notna(p_max_mw) else None,
            "P_raw_max_MW": round(float(p_raw_max_mw), 3) if pd.notna(p_raw_max_mw) else None
        })

    df_summary = pd.DataFrame(summary).sort_values(by="P_mean_kW", ascending=False).reset_index(drop=True)
    df_summary.to_csv(reports_dir / "station_summary.csv", index=False, encoding="utf-8-sig")
    return df_summary


# =============================================================================
# 9. TẠO BIẾN MỤC TIÊU CHO BÀI TOÁN HỌC CÓ GIÁM SÁT
# =============================================================================
def build_supervised_subset(df_featured_subset: pd.DataFrame, target_point: int, threshold_kw: float = 50.0):
    """

    Tạo biến mục tiêu:

    - target_P_kW: Công suất tại thời điểm t + 30 phút (kW)

    - target_P_MW: Công suất tại thời điểm t + 30 phút (MW)

    - target_delta_P_kW: Mức biến thiên Delta P = P[t+1] - P[t] (kW)

    - target_direction: Phân loại TĂNG / GIẢM / KHÔNG ĐỔI dựa trên ngưỡng biến thiên threshold_kw

      (Mặc định 50 kW ~ 0.05 MW, tương đương ~0.6% công suất trung bình trạm 8,000 kW)

    """
    sub = df_featured_subset[df_featured_subset["MaDiemDo"] == target_point].sort_values("TIME").copy()

    # Nhãn mục tiêu P[t+1] theo cả kW và MW
    sub["target_P_kW"] = sub["P_kW"].shift(-1)
    if "P_MW" in sub.columns:
        sub["target_P_MW"] = sub["P_MW"].shift(-1)
    else:
        sub["target_P_MW"] = sub["target_P_kW"] / 1000.0

    # Loại bỏ dòng không có target (dòng cuối cùng của tập)
    sub = sub.dropna(subset=["target_P_kW"]).copy()

    # Mức biến thiên và phân loại xu hướng
    sub["target_delta_P_kW"] = sub["target_P_kW"] - sub["P_kW"]
    sub["target_delta_P_MW"] = sub["target_P_MW"] - sub["P_MW"]
    sub["target_delta_P"] = sub["target_delta_P_kW"]  # Biến delta mặc định theo kW

    conditions = [
        sub["target_delta_P_kW"] > threshold_kw,
        sub["target_delta_P_kW"] < -threshold_kw
    ]
    sub["target_direction"] = np.select(conditions, ["TĂNG", "GIẢM"], default="KHÔNG ĐỔI")
    return sub


# =============================================================================
# 10. BỘ 8 BIỂU ĐỒ ĐÁNH GIÁ CHẤT LƯỢNG CAO (200 DPI - CHUẨN KHOA HỌC)
# =============================================================================
def generate_evaluation_figures(df_raw: pd.DataFrame, df_clean: pd.DataFrame, df_supervised: pd.DataFrame, target_point: int, figures_dir: Path):
    figures_dir.mkdir(parents=True, exist_ok=True)
    clean_target = df_clean[df_clean["MaDiemDo"] == target_point].set_index("TIME").sort_index()
    raw_target = df_raw[df_raw["MaDiemDo"] == target_point].set_index("TIME").sort_index()

    # --- 1. Toàn cảnh chuỗi phụ tải & Đối chiếu phục hồi nội suy ---
    fig1, (ax1a, ax1b) = plt.subplots(2, 1, figsize=(16, 9))
    sample_month = clean_target.iloc[:1440]  # 30 ngày mẫu
    ax1a.plot(sample_month.index, sample_month["P_kW"], color="#1f77b4", lw=1.2, label=f"Điểm đo {target_point} (30 ngày)")
    ax1a.set_title(f"1A. Toàn cảnh chuỗi phụ tải P (kW) sau làm sạch (Khung 30 ngày - Điểm đo {target_point})", fontsize=13, fontweight="bold")
    ax1a.set_ylabel("Công suất P (kW)")
    ax1a.legend(loc="upper right")

    t0 = clean_target.index.min()
    t1 = t0 + pd.Timedelta(days=7)
    r_slice = raw_target.loc[t0:t1]
    c_slice = clean_target.loc[t0:t1]
    ax1b.plot(r_slice.index, r_slice["P_kW"], "r.", markersize=8, alpha=0.6, label="Dữ liệu gốc (Bị thiếu / Lệch giây)")
    ax1b.plot(c_slice.index, c_slice["P_kW"], color="#2ca02c", lw=1.5, alpha=0.9, label="Đã quy chuẩn lưới 30m & Nội suy ngắn (<=2 bước)")
    ax1b.set_title("1B. Kiểm chứng phục hồi dữ liệu: So sánh Gốc vs Sau xử lý (Zoom 7 ngày)", fontsize=13, fontweight="bold")
    ax1b.set_ylabel("Công suất P (kW)")
    ax1b.legend(loc="upper right")
    fig1.tight_layout()
    fig1.savefig(figures_dir / "01_load_timeline_and_imputation.png", dpi=200)
    plt.close(fig1)

    # --- 2. Tỷ lệ khuyết dữ liệu gốc (%) ---
    fig2, ax2 = plt.subplots(figsize=(11, 4.5))
    missing_cols = [c for c in ["UA", "UB", "UC", "IA", "IB", "IC", "CosF", "P_kW", "Q_kVAr"] if c in df_raw.columns]
    missing_p = (df_raw[missing_cols].isna().mean() * 100).sort_values(ascending=True)
    bars2 = ax2.barh(missing_p.index, missing_p.values, color=sns.color_palette("Blues_r", len(missing_p)))
    ax2.set_title("2. Đánh giá tỷ lệ khuyết dữ liệu gốc theo từng biến số (%)", fontsize=13, fontweight="bold")
    for b in bars2:
        w = b.get_width()
        ax2.text(w + 0.05, b.get_y() + b.get_height()/2, f"{w:.2f}%", va="center", fontsize=10, fontweight="bold")
    fig2.tight_layout()
    fig2.savefig(figures_dir / "02_data_completeness_missing_matrix.png", dpi=200)
    plt.close(fig2)

    # --- 3. Phân phối xác suất & Boxplot ---
    fig3, axes3 = plt.subplots(1, 2, figsize=(15, 5.5))
    p_data = clean_target["P_kW"].dropna()
    p_mean_val = p_data.mean()
    p_med_val = p_data.median()
    sns.histplot(p_data, bins=60, kde=True, ax=axes3[0], color="#1f77b4", edgecolor="black", alpha=0.6)
    axes3[0].axvline(p_mean_val, color="red", linestyle="--", lw=2, label=f"Mean = {p_mean_val:,.1f} kW ({p_mean_val/1000:.2f} MW)")
    axes3[0].axvline(p_med_val, color="green", linestyle="-.", lw=2, label=f"Median = {p_med_val:,.1f} kW ({p_med_val/1000:.2f} MW)")
    axes3[0].set_title("Phân phối xác suất (Histogram & KDE) của P (kW)", fontsize=12, fontweight="bold")
    axes3[0].set_xlabel("Công suất P (kW)")
    axes3[0].legend()

    sns.boxplot(y=p_data, ax=axes3[1], color="#98df8a")
    axes3[1].set_title("Biểu đồ Boxplot: Khoảng tứ phân vị & Ngoại lai của P (kW)", fontsize=12, fontweight="bold")
    axes3[1].set_ylabel("Công suất P (kW)")
    fig3.suptitle("3. Kiểm định phân phối và biên độ công suất thực tế P (kW)", fontsize=15, fontweight="bold")
    fig3.tight_layout()
    fig3.savefig(figures_dir / "03_power_distribution_and_anomaly.png", dpi=200)
    plt.close(fig3)

    # --- 4. Ma trận tương quan đặc trưng ---
    fig4, ax4 = plt.subplots(figsize=(10, 7))
    corr_cols = [c for c in ["P_kW", "target_P_kW", "lag_1", "lag_2", "lag_48", "rolling_mean_3", "U_avg", "I_avg", "S_kVA", "CosF"] if c in df_supervised.columns]
    sns.heatmap(df_supervised[corr_cols].corr(), annot=True, fmt=".2f", cmap="vlag", ax=ax4)
    ax4.set_title("4. Ma trận tương quan đặc trưng với biến mục tiêu (target_P_kW)", fontsize=13, fontweight="bold")
    fig4.tight_layout()
    fig4.savefig(figures_dir / "04_correlation_heatmap_features.png", dpi=200)
    plt.close(fig4)

    # --- 5. Đánh giá biến mục tiêu (Delta P & Nhãn xu hướng) ---
    fig5, axes5 = plt.subplots(1, 2, figsize=(15, 5.5))
    sns.histplot(df_supervised["target_delta_P_kW"].dropna(), bins=70, kde=True, ax=axes5[0], color="#9467bd")
    axes5[0].axvline(0, color="red", linestyle="--", lw=1.5)
    axes5[0].set_title("Phân phối Delta P = P[t+1] - P[t] (kW)", fontsize=12, fontweight="bold")
    axes5[0].set_xlabel("Delta P (kW)")

    d_counts = df_supervised["target_direction"].value_counts()
    bars5 = axes5[1].bar(d_counts.index, d_counts.values, color=["#2ca02c", "#d62728", "#7f7f7f"][:len(d_counts)], edgecolor="black")
    axes5[1].set_title("Tỷ lệ nhãn xu hướng (Ngưỡng biến thiên |Delta P| >= 50 kW)", fontsize=12, fontweight="bold")
    for b in bars5:
        h = b.get_height()
        axes5[1].text(b.get_x() + b.get_width()/2, h + 10, f"{h:,}\n({h/len(df_supervised)*100:.1f}%)", ha="center", fontsize=10, fontweight="bold")
    fig5.suptitle("5. Kiểm định biến mục tiêu dự báo (Target Verification)", fontsize=15, fontweight="bold")
    fig5.tight_layout()
    fig5.savefig(figures_dir / "05_target_evaluation_and_direction.png", dpi=200)
    plt.close(fig5)

    # --- 6. Scatter P[t] vs P[t+1] ---
    fig6, ax6 = plt.subplots(figsize=(7, 7))
    samp = df_supervised.sample(n=min(3000, len(df_supervised)), random_state=42)
    ax6.scatter(samp["P_kW"], samp["target_P_kW"], alpha=0.3, color="#1f77b4", s=25)
    m_val = min(samp["P_kW"].min(), samp["target_P_kW"].min())
    mx_val = max(samp["P_kW"].max(), samp["target_P_kW"].max())
    ax6.plot([m_val, mx_val], [m_val, mx_val], "r--", lw=2, label="Đường y = x")
    r = df_supervised[["P_kW", "target_P_kW"]].corr().iloc[0, 1]
    ax6.set_title(f"6. Tương quan tuyến tính P[t] vs P[t+1] (r = {r:.4f})", fontsize=13, fontweight="bold")
    ax6.set_xlabel("Công suất hiện tại P[t] (kW)")
    ax6.set_ylabel("Công suất sau 30 phút P[t+1] (kW)")
    ax6.legend(loc="upper left")
    fig6.tight_layout()
    fig6.savefig(figures_dir / "06_scatter_lag_vs_target.png", dpi=200)
    plt.close(fig6)

    # --- 7. Tính mùa vụ Giờ & Thứ ---
    fig7, axes7 = plt.subplots(1, 2, figsize=(16, 5.5))
    h_m = clean_target.groupby(clean_target.index.hour)["P_kW"].mean()
    h_s = clean_target.groupby(clean_target.index.hour)["P_kW"].std()
    axes7[0].plot(h_m.index, h_m.values, marker="o", color="#ff7f0e", lw=2.5)
    axes7[0].fill_between(h_m.index, h_m.values - h_s.values, h_m.values + h_s.values, color="#ff7f0e", alpha=0.2)
    axes7[0].set_title("Phụ tải trung bình theo 24 giờ trong ngày (kW)", fontsize=12, fontweight="bold")
    axes7[0].set_xlabel("Giờ trong ngày (0 - 23)")
    axes7[0].set_ylabel("Công suất P trung bình (kW)")
    axes7[0].set_xticks(range(0, 24))

    dow_names = ["Thứ 2", "Thứ 3", "Thứ 4", "Thứ 5", "Thứ 6", "Thứ 7", "Chủ Nhật"]
    d_m = clean_target.groupby(clean_target.index.dayofweek)["P_kW"].mean()
    axes7[1].bar(dow_names[:len(d_m)], d_m.values, color=sns.color_palette("Set2", len(d_m)), edgecolor="black")
    axes7[1].set_title("Phụ tải trung bình theo ngày trong tuần (kW)", fontsize=12, fontweight="bold")
    axes7[1].set_ylabel("Công suất P trung bình (kW)")
    fig7.suptitle("7. Đánh giá tính mùa vụ phụ tải điện (Seasonality Profile)", fontsize=15, fontweight="bold")
    fig7.tight_layout()
    fig7.savefig(figures_dir / "07_daily_and_weekly_seasonality.png", dpi=200)
    plt.close(fig7)

    # --- 8. Timeline chia tập Train - Val - Test (Thực tế đến tháng 8/2025) ---
    fig8, ax8 = plt.subplots(figsize=(15, 4))
    df_sup_sorted = df_supervised.sort_values("TIME").reset_index(drop=True)
    n = len(df_sup_sorted)
    t_idx = int(n * 0.70)
    v_idx = int(n * 0.85)

    ax8.plot(df_sup_sorted.iloc[:t_idx]["TIME"], df_sup_sorted.iloc[:t_idx]["P_kW"], color="#1f77b4", label="Train (70%)", lw=1)
    ax8.plot(df_sup_sorted.iloc[t_idx:v_idx]["TIME"], df_sup_sorted.iloc[t_idx:v_idx]["P_kW"], color="#ff7f0e", label="Validation (15%)", lw=1)
    ax8.plot(df_sup_sorted.iloc[v_idx:]["TIME"], df_sup_sorted.iloc[v_idx:]["P_kW"], color="#2ca02c", label="Test (15%)", lw=1)
    ax8.set_title("8. Phân chia tập dữ liệu tuần tự theo thời gian (Thực tế: 03/2020 -> 08/2025)", fontsize=13, fontweight="bold")
    ax8.set_xlabel("Thời gian")
    ax8.set_ylabel("Công suất P (kW)")
    ax8.legend(loc="upper right")
    fig8.tight_layout()
    fig8.savefig(figures_dir / "08_data_split_timeline.png", dpi=200)
    plt.close(fig8)


# =============================================================================
# HÀM THỰC THI CHÍNH TOÀN BỘ PIPELINE
# =============================================================================
def run_pipeline(

    data_dir_str: str,

    output_dir_str: str = "output_preprocessing",

    target_point: int = None,

    threshold_kw: float = 50.0

):
    base_out = Path(output_dir_str)
    data_out = base_out / "data"
    reports_out = base_out / "reports"
    figures_out = base_out / "figures"

    for d in [data_out, reports_out, figures_out]:
        d.mkdir(parents=True, exist_ok=True)

    data_dir = Path(data_dir_str)

    print("=" * 80)
    print("PIPELINE TIỀN XỬ LÝ PHỤ TẢI ĐIỆN 30 PHÚT (BẢN SỬA LỖI THỜI GIAN & CHIA TẬP)")
    print(f"📁 Dữ liệu đầu vào:      {data_dir.resolve()}")
    print(f"💾 Thư mục dữ liệu sạch: {data_out.resolve()}")
    print(f"📑 Thư mục báo cáo CSV:  {reports_out.resolve()}")
    print(f"📊 Thư mục biểu đồ (DPI cao): {figures_out.resolve()}")
    print("=" * 80)

    # 1. Quét CSV & đọc chuẩn hóa
    csv_files = find_csv_files(data_dir)
    print(f"\n[1/7] Tìm thấy {len(csv_files)} file CSV hợp lệ. Đang chuẩn hóa...")
    dfs = []
    for f in csv_files:
        df_item = load_and_standardize_single_file(f)
        if len(df_item) > 0:
            dfs.append(df_item)

    df_raw = pd.concat(dfs, ignore_index=True)
    t_min_real = df_raw["TIME"].min()
    t_max_real = df_raw["TIME"].max()
    print(f"      -> Tổng số dòng dữ liệu thô hợp lệ: {len(df_raw):,}")
    print(f"      -> Dải thời gian thực tế: {t_min_real} đến {t_max_real}")

    # 2. Báo cáo chất lượng thô
    print("\n[2/7] Xuất báo cáo chất lượng: reports/quality_report.csv...")
    audit_summary = audit_data_quality(df_raw, reports_out)

    # 3. Tạo lưới 30 phút CHƯA NỘI SUY
    print("\n[3/7] Quy chuẩn lưới thời gian 30 phút đồng nhất cho từng điểm đo (CHƯA nội suy)...")
    df_grid = build_regular_30m_grid(df_raw)
    total_slots = df_grid["TIME"].nunique()
    print(f"      -> Tổng số mốc 30 phút chuẩn: {total_slots:,} mốc")

    # 4. CHIA TẬP TRAIN / VAL / TEST TRƯỚC KHI NỘI SUY (CHỐNG RÒ RỈ DỮ LIỆU)
    print("\n[4/7] Chia Train (70%) - Validation (15%) - Test (15%) THEO THỜI GIAN trước khi nội suy...")
    train_raw_grid, val_raw_grid, test_raw_grid, split_meta = split_timeline_into_subsets(df_grid, 0.70, 0.15)
    print(f"      -> Train subset: {len(train_raw_grid):,} dòng ({split_meta['train_start']} -> {split_meta['train_end']})")
    print(f"      -> Val subset:   {len(val_raw_grid):,} dòng ({split_meta['val_start']} -> {split_meta['val_end']})")
    print(f"      -> Test subset:  {len(test_raw_grid):,} dòng ({split_meta['test_start']} -> {split_meta['test_end']})")

    # 5. NỘI SUY ĐỘC LẬP TRONG TỪNG TẬP (ISOLATED IMPUTATION)
    print("\n[5/7] Nội suy kiểm soát (<= 2 bước) độc lập trong từng tập (Zero Data Leakage)...")
    train_clean = impute_isolated_subset(train_raw_grid, limit_steps=2)
    val_clean = impute_isolated_subset(val_raw_grid, limit_steps=2)
    test_clean = impute_isolated_subset(test_raw_grid, limit_steps=2)

    train_clean["split"] = "train"
    val_clean["split"] = "val"
    test_clean["split"] = "test"

    df_clean_long = pd.concat([train_clean, val_clean, test_clean], ignore_index=True)
    df_clean_long.to_csv(data_out / "cleaned_long.csv", index=False, encoding="utf-8-sig")

    df_wide = df_clean_long.pivot_table(index="TIME", columns="MaDiemDo", values=["P_kW", "P_MW", "Q_kVAr"], aggfunc="first")
    df_wide.columns = [f"{var}_{pt}" for (var, pt) in df_wide.columns]
    df_wide.reset_index().to_csv(data_out / "cleaned_wide.csv", index=False, encoding="utf-8-sig")

    # 6. Trích xuất đặc trưng & Báo cáo trạm
    print("\n[6/7] Trích xuất đặc trưng vật lý & chuỗi thời gian, xuất station_summary.csv...")
    df_featured = engineer_features(df_clean_long)
    df_featured.to_csv(data_out / "featured_data.csv", index=False, encoding="utf-8-sig")

    df_station_summary = generate_station_summary(df_raw, df_grid, df_clean_long, reports_out)
    if target_point is None:
        target_point = int(df_station_summary.iloc[0]["MaDiemDo"])
    print(f"      -> Trạm mục tiêu lớn nhất được chọn: MaDiemDo = {target_point} ({df_station_summary.iloc[0]['TenDiemDo']})")

    # Tạo tập học có giám sát và chia train / val / test
    feat_train = df_featured[df_featured["split"] == "train"]
    feat_val = df_featured[df_featured["split"] == "val"]
    feat_test = df_featured[df_featured["split"] == "test"]

    sup_train = build_supervised_subset(feat_train, target_point, threshold_kw)
    sup_val = build_supervised_subset(feat_val, target_point, threshold_kw)
    sup_test = build_supervised_subset(feat_test, target_point, threshold_kw)

    df_supervised = pd.concat([sup_train, sup_val, sup_test], ignore_index=True)
    df_supervised.to_csv(data_out / "supervised_data.csv", index=False, encoding="utf-8-sig")
    sup_train.to_csv(data_out / "train.csv", index=False, encoding="utf-8-sig")
    sup_val.to_csv(data_out / "validation.csv", index=False, encoding="utf-8-sig")
    sup_test.to_csv(data_out / "test.csv", index=False, encoding="utf-8-sig")

    print(f"      -> Supervised Data: Train {len(sup_train):,} | Val {len(sup_val):,} | Test {len(sup_test):,}")

    # 7. Sinh 8 biểu đồ đánh giá chuẩn khoa học vào figures/
    print("\n[7/7] Sinh 8 biểu đồ đánh giá chất lượng cao (DPI 200) vào figures/...")
    generate_evaluation_figures(df_raw, df_clean_long, df_supervised, target_point, figures_out)

    # Lưu metadata tổng kết
    summary_final = {
        "timestamp": datetime.now().isoformat(),
        "total_files_read": len(csv_files),
        "raw_rows": len(df_raw),
        "clean_rows": len(df_clean_long),
        "selected_target_point": target_point,
        "time_range": {
            "start": str(t_min_real),
            "end": str(t_max_real)
        },
        "split_info": {
            "train_samples": len(sup_train),
            "train_range": [split_meta["train_start"], split_meta["train_end"]],
            "val_samples": len(sup_val),
            "val_range": [split_meta["val_start"], split_meta["val_end"]],
            "test_samples": len(sup_test),
            "test_range": [split_meta["test_start"], split_meta["test_end"]],
        },
        "paths": {
            "clean_data_dir": str(data_out.resolve()),
            "reports_dir": str(reports_out.resolve()),
            "figures_dir": str(figures_out.resolve())
        }
    }
    with open(reports_out / "preprocessing_summary.json", "w", encoding="utf-8") as f:
        json.dump(summary_final, f, indent=4, ensure_ascii=False)

    print("\n" + "=" * 80)
    print("✅ ĐÃ HOÀN TẤT VÀ KHẮC PHỤC TRIỆT ĐỂ LỖI THỜI GIAN & CHIA TẬP KHÔNG RÒ RỈ!")
    print("=" * 80)
    return summary_final


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description="Pipeline tiền xử lý phụ tải điện 30 phút")
    parser.add_argument("--data_dir", type=str, default=None)
    parser.add_argument("--output_dir", type=str, default="output_preprocessing")
    parser.add_argument("--target_point", type=int, default=None)
    parser.add_argument("--threshold_kw", type=float, default=50.0)
    args = parser.parse_args()

    if args.data_dir:
        input_path = args.data_dir
    elif Path("/content/drive/MyDrive/Tam_TLU/datasets").exists():
        input_path = "/content/drive/MyDrive/Tam_TLU/datasets"
    elif Path("./datasets").exists():
        input_path = "./datasets"
    else:
        input_path = "."

    run_pipeline(
        data_dir_str=input_path,
        output_dir_str=args.output_dir,
        target_point=args.target_point,
        threshold_kw=args.threshold_kw
    )