diagnostic_prediction.py

"""
PINN-VLE 預測診斷腳本
=====================
用途：分析預測結果不佳的根本原因
- 診斷 1：Sigma-profile 檔案缺失分析
- 診斷 2：按二元系統分組的誤差分析
- 診斷 3：訓練集 vs 預測集的分布比較
- 診斷 4：特徵範圍與外推風險檢查

使用方式：
    python diagnostic_prediction.py
    或在 Jupyter Notebook 中逐段執行
"""

import os
import csv
import numpy as np
import pandas as pd
from pathlib import Path
from collections import Counter

# ============================================================
# 設定路徑（請根據實際環境修改）
# ============================================================
ROOT = Path(__file__).resolve().parent.parent
DATA_DIR = ROOT / "data"
SIGMA_DIR = DATA_DIR / "s-profile-area"
PREDICTION_CSV = DATA_DIR / "grouping_alldata_V3.csv"
TRAINING_GROUP_CSV = DATA_DIR / "train_group_all_V2.csv"
COMPOUND_MAP_CSV = DATA_DIR / "COSMO_CID_SMILES.csv"
OUTPUT_DIR = ROOT / "code" / "diagnostic_output"
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

SIGMA_PROFILE_DIM = 51


def load_prediction_data() -> pd.DataFrame:
    df = pd.read_csv(PREDICTION_CSV)
    df.columns = df.columns.str.strip()
    for col in ["comp.A", "comp.B"]:
        df[col] = df[col].astype(int)
    return df


def load_training_groups() -> pd.DataFrame:
    df = pd.read_csv(TRAINING_GROUP_CSV)
    df.columns = df.columns.str.strip()
    for col in ["comp.A", "comp.B"]:
        df[col] = df[col].astype(int)
    return df


def load_compound_map() -> pd.DataFrame:
    df = pd.read_csv(COMPOUND_MAP_CSV)
    df.columns = df.columns.str.strip()
    return df


# ============================================================
# 診斷 1：Sigma-Profile 檔案缺失分析
# ============================================================
def diagnose_missing_sigma_profiles(df_pred: pd.DataFrame) -> pd.DataFrame:
    print("=" * 60)
    print("診斷 1：Sigma-Profile 檔案缺失分析")
    print("=" * 60)

    all_compounds = set(df_pred["comp.A"].unique()) | set(df_pred["comp.B"].unique())
    print(f"預測資料集中共有 {len(all_compounds)} 個不重複化合物")

    missing = []
    existing = []
    for cid in sorted(all_compounds):
        fpath = SIGMA_DIR / f"{cid}-opt-b3lyp.csv"
        if fpath.exists():
            existing.append(cid)
        else:
            missing.append(cid)

    print(f"有 sigma-profile 的化合物：{len(existing)}")
    print(f"缺少 sigma-profile 的化合物：{len(missing)}")

    if missing:
        print(f"\n缺失的化合物 ID 清單：")
        for cid in missing:
            print(f"  - {cid}")

    affected_a = df_pred["comp.A"].isin(missing)
    affected_b = df_pred["comp.B"].isin(missing)
    affected_any = affected_a | affected_b
    affected_both = affected_a & affected_b

    total = len(df_pred)
    n_any = affected_any.sum()
    n_both = affected_both.sum()
    n_a_only = (affected_a & ~affected_b).sum()
    n_b_only = (affected_b & ~affected_a).sum()
    n_clean = total - n_any

    print(f"\n對預測資料的影響：")
    print(f"  總資料筆數：           {total}")
    print(f"  兩者皆有 sigma-profile：{n_clean} ({100*n_clean/total:.1f}%)")
    print(f"  至少一方缺失：         {n_any} ({100*n_any/total:.1f}%)")
    print(f"    僅 comp.A 缺失：     {n_a_only}")
    print(f"    僅 comp.B 缺失：     {n_b_only}")
    print(f"    兩者皆缺失：         {n_both}")

    df_pred["missing_sigma_A"] = affected_a
    df_pred["missing_sigma_B"] = affected_b
    df_pred["missing_sigma_any"] = affected_any

    if missing:
        compound_map = load_compound_map()
        missing_info = compound_map[compound_map["COSMO"].isin(missing)][["COSMO", "IUPAC", "SMILES"]]
        if not missing_info.empty:
            print(f"\n缺失化合物的詳細資訊：")
            print(missing_info.to_string(index=False))
        out_path = OUTPUT_DIR / "missing_sigma_profiles.csv"
        pd.DataFrame({"compound_id": missing}).to_csv(out_path, index=False)
        print(f"\n缺失清單已儲存至：{out_path}")

    print()
    return df_pred


# ============================================================
# 診斷 2：按二元系統分組的誤差分析
# ============================================================
def diagnose_error_by_system(df_pred: pd.DataFrame) -> pd.DataFrame:
    print("=" * 60)
    print("診斷 2：按二元系統分組的誤差分析")
    print("=" * 60)

    df = df_pred.copy()
    df["system"] = df["comp.A"].astype(str) + "-" + df["comp.B"].astype(str)

    has_exp = df["P_exp"].notna() & df["yA_exp"].notna()
    df_valid = df[has_exp].copy()

    P_exp = df_valid["P_exp"].values.astype(float)
    y_exp = df_valid["yA_exp"].values.astype(float)

    P_kPa = np.where(P_exp > 500, P_exp / 1000.0, P_exp)
    df_valid["P_exp_kPa"] = P_kPa

    print(f"\n壓力單位統計：")
    n_pa = (P_exp > 500).sum()
    n_kpa = (P_exp <= 500).sum()
    print(f"  被判定為 Pa（P > 500）的資料筆數：{n_pa} ({100*n_pa/len(P_exp):.1f}%)")
    print(f"  被判定為 kPa（P <= 500）的資料筆數：{n_kpa} ({100*n_kpa/len(P_exp):.1f}%)")

    ambiguous = (P_exp > 100) & (P_exp <= 1000)
    n_ambiguous = ambiguous.sum()
    if n_ambiguous > 0:
        print(f"\n  ⚠ 壓力值在 100-1000 之間（單位模糊區間）的資料：{n_ambiguous} 筆")
        print(f"    這些資料的壓力單位判斷可能有誤！")

    system_stats = df_valid.groupby("system").agg(
        count=("P_exp_kPa", "size"),
        T_min=("T", "min"),
        T_max=("T", "max"),
        xA_min=("xA", "min"),
        xA_max=("xA", "max"),
        P_min=("P_exp_kPa", "min"),
        P_max=("P_exp_kPa", "max"),
        missing_sigma=("missing_sigma_any", "any"),
    ).reset_index()

    system_stats = system_stats.sort_values("count", ascending=False)

    print(f"\n資料量前 20 的二元系統：")
    top20 = system_stats.head(20)
    print(top20[["system", "count", "T_min", "T_max", "P_min", "P_max", "missing_sigma"]].to_string(index=False))

    missing_systems = system_stats[system_stats["missing_sigma"]]
    if not missing_systems.empty:
        print(f"\n⚠ 含有缺失 sigma-profile 的系統數量：{len(missing_systems)}")
        affected_points = missing_systems["count"].sum()
        print(f"  受影響的資料筆數：{affected_points} ({100*affected_points/len(df_valid):.1f}%)")

    out_path = OUTPUT_DIR / "system_statistics.csv"
    system_stats.to_csv(out_path, index=False)
    print(f"\n系統統計已儲存至：{out_path}")

    print()
    return df_valid


# ============================================================
# 診斷 3：訓練集 vs 預測集的分布比較
# ============================================================
def diagnose_distribution_comparison(df_pred: pd.DataFrame):
    print("=" * 60)
    print("診斷 3：訓練集 vs 預測集的分布比較")
    print("=" * 60)

    if not TRAINING_GROUP_CSV.exists():
        print("⚠ 找不到訓練分組檔案，跳過此診斷。")
        print()
        return

    df_train = load_training_groups()

    train_compounds_A = set(df_train["comp.A"].unique())
    train_compounds_B = set(df_train["comp.B"].unique())
    train_compounds = train_compounds_A | train_compounds_B

    pred_compounds_A = set(df_pred["comp.A"].unique())
    pred_compounds_B = set(df_pred["comp.B"].unique())
    pred_compounds = pred_compounds_A | pred_compounds_B

    common = train_compounds & pred_compounds
    pred_only = pred_compounds - train_compounds
    train_only = train_compounds - pred_compounds

    print(f"\n化合物覆蓋分析：")
    print(f"  訓練集化合物數：{len(train_compounds)}")
    print(f"  預測集化合物數：{len(pred_compounds)}")
    print(f"  兩者共有：{len(common)}")
    print(f"  僅在預測集中（新化合物）：{len(pred_only)}")
    print(f"  僅在訓練集中：{len(train_only)}")

    if pred_only:
        print(f"\n  預測集中的新化合物（訓練時未見過）：")
        for cid in sorted(pred_only):
            sigma_exists = (SIGMA_DIR / f"{cid}-opt-b3lyp.csv").exists()
            status = "有" if sigma_exists else "缺"
            print(f"    - {cid} (sigma-profile: {status})")

    affected_new = df_pred["comp.A"].isin(pred_only) | df_pred["comp.B"].isin(pred_only)
    n_new = affected_new.sum()
    print(f"\n  涉及新化合物的預測資料筆數：{n_new} ({100*n_new/len(df_pred):.1f}%)")

    train_systems = set(
        df_train.apply(lambda r: f"{r['comp.A']}-{r['comp.B']}", axis=1)
    )
    pred_systems = set(
        df_pred.apply(lambda r: f"{r['comp.A']}-{r['comp.B']}", axis=1)
    )

    common_sys = train_systems & pred_systems
    pred_only_sys = pred_systems - train_systems
    print(f"\n二元系統覆蓋分析：")
    print(f"  訓練集系統數：{len(train_systems)}")
    print(f"  預測集系統數：{len(pred_systems)}")
    print(f"  兩者共有：{len(common_sys)}")
    print(f"  僅在預測集中（新系統）：{len(pred_only_sys)}")

    if "T" in df_train.columns:
        print(f"\n溫度分布比較：")
        print(f"  訓練集 T 範圍：{df_train['T'].min():.2f} ~ {df_train['T'].max():.2f} K")
        print(f"  預測集 T 範圍：{df_pred['T'].min():.2f} ~ {df_pred['T'].max():.2f} K")

        pred_T_below = (df_pred["T"] < df_train["T"].min()).sum()
        pred_T_above = (df_pred["T"] > df_train["T"].max()).sum()
        if pred_T_below > 0 or pred_T_above > 0:
            print(f"  ⚠ 預測集溫度超出訓練範圍：低於下限 {pred_T_below} 筆，高於上限 {pred_T_above} 筆")

    if pred_only_sys:
        new_sys_df = df_pred[
            df_pred.apply(lambda r: f"{r['comp.A']}-{r['comp.B']}", axis=1).isin(pred_only_sys)
        ]
        n_new_sys_points = len(new_sys_df)
        print(f"\n  涉及新系統的預測資料筆數：{n_new_sys_points} ({100*n_new_sys_points/len(df_pred):.1f}%)")

        out_path = OUTPUT_DIR / "new_systems_not_in_training.csv"
        pd.DataFrame({"system": sorted(pred_only_sys)}).to_csv(out_path, index=False)
        print(f"  新系統清單已儲存至：{out_path}")

    print()


# ============================================================
# 診斷 4：特徵範圍與外推風險檢查
# ============================================================
def diagnose_feature_ranges(df_pred: pd.DataFrame):
    print("=" * 60)
    print("診斷 4：特徵範圍與外推風險檢查")
    print("=" * 60)

    print("\n物理條件範圍：")
    for col in ["T", "xA", "yA_exp", "P_exp"]:
        if col in df_pred.columns:
            vals = pd.to_numeric(df_pred[col], errors="coerce").dropna()
            print(f"  {col:>8s}: min={vals.min():12.4f}  max={vals.max():12.4f}  mean={vals.mean():12.4f}  std={vals.std():12.4f}")

    print(f"\n壓力分布區間：")
    P_vals = pd.to_numeric(df_pred["P_exp"], errors="coerce").dropna()
    bins = [0, 1, 10, 50, 100, 500, 1000, 5000, 10000, 100000, 1e9]
    labels = ["0-1", "1-10", "10-50", "50-100", "100-500", "500-1k", "1k-5k", "5k-10k", "10k-100k", ">100k"]
    P_binned = pd.cut(P_vals, bins=bins, labels=labels, right=True)
    dist = P_binned.value_counts().sort_index()
    for label, count in dist.items():
        if count > 0:
            bar = "#" * min(count // 50, 50)
            print(f"  {label:>10s}: {count:>6d} {bar}")

    print(f"\n  ⚠ 壓力 > 500 的資料被自動判定為 Pa 並除以 1000")
    print(f"    但壓力在 500-1000 之間的資料，可能本身就是 kPa！")
    print(f"    這會導致這些資料的壓力被錯誤地縮小 1000 倍。")

    if "VpA" in df_pred.columns and "VpB" in df_pred.columns:
        print(f"\n蒸氣壓範圍：")
        VpA = pd.to_numeric(df_pred["VpA"], errors="coerce").dropna()
        VpB = pd.to_numeric(df_pred["VpB"], errors="coerce").dropna()
        print(f"  VpA: min={VpA.min():.4f}  max={VpA.max():.4f}")
        print(f"  VpB: min={VpB.min():.4f}  max={VpB.max():.4f}")

        n_vpa_zero = (VpA <= 0).sum()
        n_vpb_zero = (VpB <= 0).sum()
        if n_vpa_zero > 0 or n_vpb_zero > 0:
            print(f"  ⚠ VpA <= 0 的資料：{n_vpa_zero} 筆（ln 轉換會出問題）")
            print(f"  ⚠ VpB <= 0 的資料：{n_vpb_zero} 筆（ln 轉換會出問題）")

    if "gammaA" in df_pred.columns and "gammaB" in df_pred.columns:
        print(f"\n活度係數範圍：")
        gA = pd.to_numeric(df_pred["gammaA"], errors="coerce").dropna()
        gB = pd.to_numeric(df_pred["gammaB"], errors="coerce").dropna()
        print(f"  gammaA: min={gA.min():.6f}  max={gA.max():.4f}  mean={gA.mean():.4f}")
        print(f"  gammaB: min={gB.min():.6f}  max={gB.max():.4f}  mean={gB.mean():.4f}")

        n_extreme_gA = ((gA > 100) | (gA < 0.01)).sum()
        n_extreme_gB = ((gB > 100) | (gB < 0.01)).sum()
        if n_extreme_gA > 0 or n_extreme_gB > 0:
            print(f"  ⚠ gammaA 極端值（<0.01 或 >100）：{n_extreme_gA} 筆")
            print(f"  ⚠ gammaB 極端值（<0.01 或 >100）：{n_extreme_gB} 筆")
            print(f"    這些系統高度非理想，Margules 模型可能不適用")

    print(f"\n組成 xA 分布：")
    xA = df_pred["xA"].values.astype(float)
    near_zero = (xA < 0.01).sum()
    near_one = (xA > 0.99).sum()
    print(f"  xA < 0.01（接近純 B）：{near_zero} 筆")
    print(f"  xA > 0.99（接近純 A）：{near_one} 筆")
    if near_zero + near_one > 0:
        print(f"  ⚠ 極端組成區間的預測通常較不準確")

    print()


# ============================================================
# 診斷 5：Sigma-Profile 品質檢查
# ============================================================
def diagnose_sigma_profile_quality(df_pred: pd.DataFrame):
    print("=" * 60)
    print("診斷 5：Sigma-Profile 品質檢查")
    print("=" * 60)

    all_compounds = sorted(set(df_pred["comp.A"].unique()) | set(df_pred["comp.B"].unique()))

    zero_profiles = []
    sparse_profiles = []
    normal_profiles = []
    missing_profiles = []

    for cid in all_compounds:
        fpath = SIGMA_DIR / f"{cid}-opt-b3lyp.csv"
        if not fpath.exists():
            missing_profiles.append(cid)
            continue

        try:
            with open(fpath, "r", encoding="utf-8") as f:
                vals = [float(row[1]) for row in csv.reader(f)][:SIGMA_PROFILE_DIM]
        except (IndexError, ValueError):
            missing_profiles.append(cid)
            continue

        arr = np.array(vals)
        if np.allclose(arr, 0):
            zero_profiles.append(cid)
        elif np.count_nonzero(arr) < 10:
            sparse_profiles.append(cid)
        else:
            normal_profiles.append(cid)

    print(f"\n  正常的 sigma-profile：{len(normal_profiles)}")
    print(f"  稀疏的 sigma-profile（非零值 < 10 個）：{len(sparse_profiles)}")
    print(f"  全零的 sigma-profile：{len(zero_profiles)}")
    print(f"  檔案缺失或無法讀取：{len(missing_profiles)}")

    if zero_profiles:
        print(f"\n  全零 sigma-profile 的化合物：{zero_profiles}")
    if sparse_profiles:
        print(f"  稀疏 sigma-profile 的化合物：{sparse_profiles}")

    print()


# ============================================================
# 診斷 6：分組誤差分析（按化學分類）
# ============================================================
def diagnose_error_by_chemical_group(df_pred: pd.DataFrame):
    print("=" * 60)
    print("診斷 6：按化學分類的資料分布分析")
    print("=" * 60)

    group_files = list((ROOT / "grouping_COSMO").glob("*.csv"))
    if not group_files:
        print("⚠ 找不到分組 CSV 檔案，跳過此診斷。")
        print()
        return

    group_data = []
    for gf in sorted(group_files):
        gdf = pd.read_csv(gf)
        gdf.columns = gdf.columns.str.strip()
        group_name = gf.stem
        group_data.append({
            "group": group_name,
            "count": len(gdf),
            "systems": gdf.apply(lambda r: f"{r['comp.A']}-{r['comp.B']}", axis=1).nunique(),
            "T_min": gdf["T"].min(),
            "T_max": gdf["T"].max(),
        })

    df_groups = pd.DataFrame(group_data).sort_values("count", ascending=False)
    total = df_groups["count"].sum()

    print(f"\n化學分類群組統計（共 {total} 筆）：")
    print(f"{'群組':<20s} {'資料筆數':>8s} {'佔比':>8s} {'系統數':>6s} {'T範圍':>20s}")
    print("-" * 65)
    for _, row in df_groups.iterrows():
        pct = 100 * row["count"] / total
        t_range = f"{row['T_min']:.1f}-{row['T_max']:.1f} K"
        print(f"{row['group']:<20s} {row['count']:>8d} {pct:>7.1f}% {row['systems']:>6d} {t_range:>20s}")

    out_path = OUTPUT_DIR / "chemical_group_statistics.csv"
    df_groups.to_csv(out_path, index=False)
    print(f"\n群組統計已儲存至：{out_path}")
    print()


# ============================================================
# 總結報告
# ============================================================
def generate_summary(df_pred: pd.DataFrame):
    print("=" * 60)
    print("診斷總結報告")
    print("=" * 60)

    total = len(df_pred)

    if "missing_sigma_any" in df_pred.columns:
        n_missing = df_pred["missing_sigma_any"].sum()
    else:
        all_compounds = set(df_pred["comp.A"].unique()) | set(df_pred["comp.B"].unique())
        missing = [c for c in all_compounds if not (SIGMA_DIR / f"{c}-opt-b3lyp.csv").exists()]
        n_missing = (df_pred["comp.A"].isin(missing) | df_pred["comp.B"].isin(missing)).sum()

    P_vals = pd.to_numeric(df_pred["P_exp"], errors="coerce").dropna()
    n_ambiguous_P = ((P_vals > 100) & (P_vals <= 1000)).sum()

    if "gammaA" in df_pred.columns:
        gA = pd.to_numeric(df_pred["gammaA"], errors="coerce").dropna()
        n_extreme_gamma = ((gA > 100) | (gA < 0.01)).sum()
    else:
        n_extreme_gamma = 0

    print(f"""
┌──────────────────────────────────────────────────────────┐
│  預測資料總筆數：{total:>8d}                              │
│                                                          │
│  問題類別                    受影響筆數     佔比          │
│  ─────────────────────────────────────────────────        │
│  缺少 sigma-profile          {n_missing:>8d}     {100*n_missing/total:>5.1f}%        │
│  壓力單位模糊（100-1000）    {n_ambiguous_P:>8d}     {100*n_ambiguous_P/total:>5.1f}%        │
│  活度係數極端值              {n_extreme_gamma:>8d}     {100*n_extreme_gamma/total:>5.1f}%        │
└──────────────────────────────────────────────────────────┘

建議的行動優先順序：
  1. 排除缺少 sigma-profile 的資料，重新計算誤差指標
  2. 檢查壓力單位模糊區間的資料，確認單位是否正確
  3. 對活度係數極端值的系統，考慮使用 NRTL/Wilson 替代 Margules
  4. 分析不在訓練集中的新系統，評估外推風險
""")

    out_path = OUTPUT_DIR / "diagnostic_summary.txt"
    with open(out_path, "w", encoding="utf-8") as f:
        f.write(f"預測診斷總結\n")
        f.write(f"{'='*40}\n")
        f.write(f"總資料筆數：{total}\n")
        f.write(f"缺少 sigma-profile：{n_missing} ({100*n_missing/total:.1f}%)\n")
        f.write(f"壓力單位模糊：{n_ambiguous_P} ({100*n_ambiguous_P/total:.1f}%)\n")
        f.write(f"活度係數極端值：{n_extreme_gamma} ({100*n_extreme_gamma/total:.1f}%)\n")
    print(f"總結報告已儲存至：{out_path}")


# ============================================================
# 主程式
# ============================================================
def main():
    print("\n" + "█" * 60)
    print("  PINN-VLE 預測品質診斷工具")
    print("█" * 60 + "\n")

    print(f"資料路徑：{DATA_DIR}")
    print(f"Sigma-profile 資料夾：{SIGMA_DIR}")
    print(f"預測輸入檔：{PREDICTION_CSV}")
    print()

    if not PREDICTION_CSV.exists():
        print(f"錯誤：找不到預測輸入檔 {PREDICTION_CSV}")
        return

    df_pred = load_prediction_data()
    print(f"成功載入預測資料：{len(df_pred)} 筆\n")

    df_pred = diagnose_missing_sigma_profiles(df_pred)
    diagnose_error_by_system(df_pred)
    diagnose_distribution_comparison(df_pred)
    diagnose_feature_ranges(df_pred)
    diagnose_sigma_profile_quality(df_pred)
    diagnose_error_by_chemical_group(df_pred)
    generate_summary(df_pred)

    print("\n診斷完成！所有輸出檔案位於：")
    print(f"  {OUTPUT_DIR}")


if __name__ == "__main__":
    main()