目錄

  1. 程式碼品質問題
  2. 架構設計問題
  3. 資料處理問題
  4. 模型設計問題
  5. 訓練流程問題
  6. 可擴展性問題
  7. 綜合改良建議

程式碼品質問題

問題 1.1：硬編碼路徑

問題描述： 多處使用硬編碼的絕對路徑，嚴重影響可移植性。

位置：

  • grouping/code/grouping_COSMO_f.ipynb: python train_exp_final = pd.read_csv(r"C:\Users\Hung\Desktop\2025_COSMO\grouping\data\train_exp_datapoint.csv") train_group_all = pd.read_csv(r"C:\Users\Hung\Desktop\2025_COSMO\grouping\data\train_group_all_V2.csv")
  • mol_to_cosmo/name_to_mol.py: python data = pd.read_csv("name.csv") mol_list = GetFileName(".\\mol\\", ".mol")

為什麼這是問題：

  1. 不可移植：程式碼無法在其他電腦或作業系統上執行
  2. 協作困難：多人協作時每個人都要修改路徑
  3. 維護困難：路徑變更時需要修改多處程式碼
  4. 錯誤風險：容易因路徑錯誤導致執行失敗

改良方案：

# 1. 使用相對路徑和 pathlib
from pathlib import Path
import os

# 定義專案根目錄
PROJECT_ROOT = Path(__file__).parent.parent
DATA_DIR = PROJECT_ROOT / "data"
MODEL_DIR = PROJECT_ROOT / "models"

# 使用配置檔
import yaml

class Config:
    def __init__(self, config_path='config.yaml'):
        with open(config_path, 'r') as f:
            self.config = yaml.safe_load(f)

    @property
    def data_dir(self):
        return Path(self.config['paths']['data_dir'])

    @property
    def model_dir(self):
        return Path(self.config['paths']['model_dir'])

# config.yaml
"""
paths:
  data_dir: ./data
  model_dir: ./models
  output_dir: ./results

  # 子目錄
  s_profiles: ${paths.data_dir}/s-profiles-all
  train_data: ${paths.data_dir}/train_exp_data_F.csv
"""

# 使用範例
config = Config()
train_data = pd.read_csv(config.data_dir / 'train_exp_data_F.csv')

問題 1.2：缺乏錯誤處理

問題描述： 許多關鍵操作缺乏 try-except 錯誤處理。

位置：

  • name_to_mol.py: python sm = c[0].isomeric_smiles # 如果查詢失敗會崩潰
  • process_gjf.py: python mol = open(f"mol/{j}", "r") # 檔案不存在會崩潰

為什麼這是問題：

  1. 容錯性差：一個錯誤會導致整個批次處理失敗
  2. 除錯困難：錯誤訊息不明確，難以定位問題
  3. 資料損失：處理中途失敗可能導致已處理資料遺失
  4. 使用者體驗差：無法友善地告知使用者問題所在

改良方案：

# 改進的 name_to_mol.py
import logging
from typing import Optional, List
from dataclasses import dataclass

# 設置日誌
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s - %(levelname)s - %(message)s',
    handlers=[
        logging.FileHandler('processing.log'),
        logging.StreamHandler()
    ]
)

@dataclass
class ProcessingResult:
    cid: str
    success: bool
    error_message: Optional[str] = None
    smiles: Optional[str] = None

def fetch_smiles_safe(cid: str, max_retries: int = 3) -> Optional[str]:
    """
    安全地從 PubChem 獲取 SMILES，帶重試機制
    """
    for attempt in range(max_retries):
        try:
            c = get_compounds(str(cid), 'cid')
            if not c:
                logging.warning(f"CID {cid}: No compound found")
                return None

            smiles = c[0].isomeric_smiles
            logging.info(f"CID {cid}: Successfully retrieved SMILES")
            return smiles

        except Exception as e:
            logging.warning(f"CID {cid}: Attempt {attempt+1}/{max_retries} failed: {e}")
            if attempt < max_retries - 1:
                time.sleep(2 ** attempt)  # 指數退避
            else:
                logging.error(f"CID {cid}: All attempts failed")
                return None

def process_molecule(cid: str, output_dir: Path) -> ProcessingResult:
    """
    處理單個分子，包含完整錯誤處理
    """
    try:
        # 1. 獲取 SMILES
        smiles = fetch_smiles_safe(cid)
        if smiles is None:
            return ProcessingResult(cid, False, "Failed to fetch SMILES")

        # 2. 轉換為 3D 結構
        try:
            mol = MolFromSmiles(smiles)
            if mol is None:
                return ProcessingResult(cid, False, "Invalid SMILES string", smiles)

            mh = AddHs(mol)
            EmbedMolecule(mh, useBasicKnowledge=True, useExpTorsionAnglePrefs=False)

        except Exception as e:
            return ProcessingResult(cid, False, f"3D generation failed: {e}", smiles)

        # 3. 幾何優化
        try:
            if not UFFHasAllMoleculeParams(mh):
                return ProcessingResult(cid, False, "UFF parameters incomplete", smiles)

            UFFOptimizeMolecule(mh)

        except Exception as e:
            return ProcessingResult(cid, False, f"Optimization failed: {e}", smiles)

        # 4. 儲存 MOL 檔
        try:
            output_path = output_dir / f"{cid}.mol"
            MolToMolFile(mh, str(output_path))

        except Exception as e:
            return ProcessingResult(cid, False, f"File save failed: {e}", smiles)

        return ProcessingResult(cid, True, None, smiles)

    except Exception as e:
        logging.exception(f"CID {cid}: Unexpected error")
        return ProcessingResult(cid, False, f"Unexpected error: {e}")

def batch_process_molecules(cid_list: List[str], output_dir: Path) -> dict:
    """
    批次處理分子，支援進度保存和續傳
    """
    output_dir.mkdir(parents=True, exist_ok=True)

    # 檢查已處理的分子
    existing = {f.stem for f in output_dir.glob("*.mol")}
    to_process = [cid for cid in cid_list if str(cid) not in existing]

    logging.info(f"Total: {len(cid_list)}, Already done: {len(existing)}, To process: {len(to_process)}")

    results = {
        'success': [],
        'failed': []
    }

    # 使用 tqdm 顯示進度
    from tqdm import tqdm

    for cid in tqdm(to_process, desc="Processing molecules"):
        result = process_molecule(str(cid), output_dir)

        if result.success:
            results['success'].append(result)
        else:
            results['failed'].append(result)
            logging.error(f"CID {cid} failed: {result.error_message}")

    # 儲存失敗清單
    if results['failed']:
        with open('failed_molecules.csv', 'w') as f:
            f.write("CID,SMILES,Error\n")
            for r in results['failed']:
                f.write(f"{r.cid},{r.smiles or 'N/A'},{r.error_message}\n")

    logging.info(f"Processing complete: {len(results['success'])} success, {len(results['failed'])} failed")

    return results

問題 1.3：缺乏程式碼模組化

問題描述： 許多功能混雜在 Jupyter Notebook 中，難以重用和測試。

為什麼這是問題：

  1. 重用性差：相同功能需要複製貼上
  2. 測試困難：Notebook 難以進行單元測試
  3. 版本控制不友善：Git 難以追蹤 Notebook 的變更
  4. 效能問題：每次都要執行所有 cell

改良方案：

# 建議的專案結構
AI_for_Chemistry/
├── src/
│   ├── __init__.py
│   ├── data/
│   │   ├── __init__.py
│   │   ├── pubchem.py          # PubChem 查詢
│   │   ├── mol_generator.py    # 3D 分子生成
│   │   ├── gaussian_writer.py  # GJF 檔案生成
│   │   └── cosmo_parser.py     # COSMO 檔案解析
│   ├── features/
│   │   ├── __init__.py
│   │   ├── sigma_profile.py    # σ-profile 提取
│   │   └── feature_engineer.py # 特徵工程
│   ├── models/
│   │   ├── __init__.py
│   │   ├── chemprop_wrapper.py # ChemProp 封裝
│   │   └── mlp_vle.py          # MLP VLE 模型
│   ├── training/
│   │   ├── __init__.py
│   │   ├── trainer.py          # 訓練邏輯
│   │   └── optimizer.py        # 超參數優化
│   └── utils/
│       ├── __init__.py
│       ├── config.py           # 配置管理
│       ├── logging.py          # 日誌工具
│       └── metrics.py          # 評估指標
├── notebooks/
│   ├── 01_data_preparation.ipynb
│   ├── 02_feature_extraction.ipynb
│   ├── 03_model_training.ipynb
│   └── 04_analysis.ipynb
├── scripts/
│   ├── prepare_data.py
│   ├── train_model.py
│   └── predict.py
├── tests/
│   ├── test_data.py
│   ├── test_features.py
│   └── test_models.py
├── config/
│   ├── config.yaml
│   └── hyperparams.yaml
├── requirements.txt
├── setup.py
└── README.md

範例模組：

# src/data/mol_generator.py
from pathlib import Path
from typing import Optional, Tuple
from rdkit import Chem
from rdkit.Chem import AllChem
import logging

class MoleculeGenerator:
    """3D 分子結構生成器"""

    def __init__(self, output_dir: Path, force_field: str = 'UFF'):
        self.output_dir = Path(output_dir)
        self.output_dir.mkdir(parents=True, exist_ok=True)
        self.force_field = force_field
        self.logger = logging.getLogger(__name__)

    def smiles_to_3d(self, smiles: str, optimize: bool = True) -> Optional[Chem.Mol]:
        """
        將 SMILES 轉換為 3D 結構

        Args:
            smiles: SMILES 字串
            optimize: 是否進行幾何優化

        Returns:
            3D 分子物件，失敗則返回 None
        """
        try:
            mol = Chem.MolFromSmiles(smiles)
            if mol is None:
                self.logger.error(f"Invalid SMILES: {smiles}")
                return None

            mol = Chem.AddHs(mol)

            # 嵌入 3D 座標
            result = AllChem.EmbedMolecule(
                mol,
                useBasicKnowledge=True,
                useExpTorsionAnglePrefs=False,
                randomSeed=42
            )

            if result != 0:
                self.logger.warning(f"Embedding failed for: {smiles}")
                return None

            # 幾何優化
            if optimize:
                if not AllChem.UFFHasAllMoleculeParams(mol):
                    self.logger.warning(f"UFF parameters incomplete: {smiles}")
                    return mol

                AllChem.UFFOptimizeMolecule(mol)

            return mol

        except Exception as e:
            self.logger.exception(f"Error processing {smiles}: {e}")
            return None

    def save_mol(self, mol: Chem.Mol, filename: str) -> bool:
        """儲存分子為 MOL 檔案"""
        try:
            output_path = self.output_dir / f"{filename}.mol"
            Chem.MolToMolFile(mol, str(output_path))
            return True
        except Exception as e:
            self.logger.error(f"Failed to save {filename}: {e}")
            return False

    def batch_process(self, smiles_dict: dict) -> Tuple[list, list]:
        """
        批次處理 SMILES

        Args:
            smiles_dict: {id: smiles}

        Returns:
            (success_ids, failed_ids)
        """
        success = []
        failed = []

        for mol_id, smiles in smiles_dict.items():
            mol = self.smiles_to_3d(smiles)
            if mol is not None and self.save_mol(mol, str(mol_id)):
                success.append(mol_id)
            else:
                failed.append(mol_id)

        return success, failed

# 使用範例
if __name__ == "__main__":
    generator = MoleculeGenerator(output_dir="output/mol")

    smiles_dict = {
        "ethanol": "CCO",
        "benzene": "c1ccccc1"
    }

    success, failed = generator.batch_process(smiles_dict)
    print(f"Success: {len(success)}, Failed: {len(failed)}")

架構設計問題

問題 2.1：兩個獨立專案缺乏整合

問題描述： ChemProp_VP 和 2025_COSMO 是兩個完全獨立的專案，無法互相利用。

為什麼這是問題：

  1. 資料重複：兩個專案可能處理相同的分子
  2. 功能重複：分子處理邏輯重複實現
  3. 無法協同：無法結合兩個模型的優勢
  4. 維護成本高：需要分別維護兩套系統

改良方案：

# 統一的架構設計
from abc import ABC, abstractmethod
from typing import Dict, Any
import torch
import tensorflow as tf

class BasePredictor(ABC):
    """預測器基類"""

    @abstractmethod
    def prepare_input(self, molecules: list, conditions: dict) -> Any:
        """準備模型輸入"""
        pass

    @abstractmethod
    def predict(self, input_data: Any) -> Dict[str, float]:
        """執行預測"""
        pass

    @abstractmethod
    def load_model(self, model_path: str):
        """載入模型"""
        pass

class ChemPropVPPredictor(BasePredictor):
    """ChemProp 蒸氣壓預測器"""

    def __init__(self, checkpoint_dir: str):
        self.model = None
        self.load_model(checkpoint_dir)

    def prepare_input(self, molecules: list, conditions: dict):
        """
        Args:
            molecules: [SMILES 列表]
            conditions: {'temperature': T}
        """
        from chemprop.data import MoleculeDataset

        data = []
        for smiles in molecules:
            data.append({
                'smiles': smiles,
                'temperature': conditions['temperature']
            })

        return MoleculeDataset(data)

    def predict(self, input_data):
        with torch.no_grad():
            preds = self.model(input_data)

        return {
            'vapor_pressure': preds.cpu().numpy()
        }

    def load_model(self, model_path: str):
        from chemprop.train import load_model
        self.model = load_model(model_path)

class COSMOVLEPredictor(BasePredictor):
    """COSMO VLE 預測器"""

    def __init__(self, model_dir: str, fold: int = 1):
        self.model = None
        self.scaler = None
        self.load_model(f"{model_dir}/MLP_Vali_COSMO_{fold}.h5")
        self.load_scaler(f"{model_dir}/pipelineMLP_COSMO.joblib")

    def prepare_input(self, molecules: list, conditions: dict):
        """
        Args:
            molecules: [{'CID_A': id, 'CID_B': id}]
            conditions: {'temperature': T, 'xA': xA}
        """
        from src.features.sigma_profile import load_sigma_profile

        features = []
        for pair in molecules:
            # 載入 σ-profiles
            profile_A = load_sigma_profile(pair['CID_A'])
            profile_B = load_sigma_profile(pair['CID_B'])

            # 組合特徵
            feature = np.concatenate([
                profile_A['sigma'],
                profile_B['sigma'],
                [profile_A['area'], profile_B['area']],
                [profile_A['volume'], profile_B['volume']],
                [conditions['temperature'], conditions['xA']]
            ])

            features.append(feature)

        return np.array(features)

    def predict(self, input_data):
        # 正規化
        X_scaled = self.scaler.transform(input_data)

        # 預測
        y_pred = self.model.predict(X_scaled)

        # 反正規化
        yA = self.inverse_transform(y_pred[:, 0], 'yA')
        P = self.inverse_transform(y_pred[:, 1], 'P')

        return {
            'vapor_composition_A': yA,
            'total_pressure': P
        }

    def load_model(self, model_path: str):
        from tensorflow.keras.models import load_model
        self.model = load_model(model_path, compile=False)

    def load_scaler(self, scaler_path: str):
        from joblib import load
        self.scaler = load(scaler_path)

class EnsemblePredictor:
    """集成預測器 - 結合多個模型"""

    def __init__(self):
        self.predictors = {}

    def add_predictor(self, name: str, predictor: BasePredictor):
        self.predictors[name] = predictor

    def predict_vapor_pressure(self, smiles: str, temperature: float) -> dict:
        """預測單一組分蒸氣壓"""
        predictor = self.predictors['chemprop_vp']

        input_data = predictor.prepare_input(
            [smiles],
            {'temperature': temperature}
        )

        return predictor.predict(input_data)

    def predict_vle(self, cid_a: str, cid_b: str, temperature: float, xa: float) -> dict:
        """預測二元混合物 VLE"""
        predictor = self.predictors['cosmo_vle']

        input_data = predictor.prepare_input(
            [{'CID_A': cid_a, 'CID_B': cid_b}],
            {'temperature': temperature, 'xA': xa}
        )

        return predictor.predict(input_data)

    def predict_full_system(self, cid_a: str, cid_b: str, temperature: float) -> dict:
        """
        完整系統預測 - 結合純組分蒸氣壓和混合物 VLE
        """
        # 1. 預測純組分蒸氣壓
        vp_a = self.predict_vapor_pressure(self.get_smiles(cid_a), temperature)
        vp_b = self.predict_vapor_pressure(self.get_smiles(cid_b), temperature)

        # 2. 預測 VLE（多個組成點）
        xa_range = np.linspace(0, 1, 11)
        vle_results = []

        for xa in xa_range:
            result = self.predict_vle(cid_a, cid_b, temperature, xa)
            vle_results.append({
                'xA': xa,
                'yA': result['vapor_composition_A'][0],
                'P': result['total_pressure'][0]
            })

        return {
            'pure_vapor_pressure': {
                'A': vp_a['vapor_pressure'][0],
                'B': vp_b['vapor_pressure'][0]
            },
            'vle_data': vle_results,
            'ideal_vle': self.calculate_ideal_vle(
                vp_a['vapor_pressure'][0],
                vp_b['vapor_pressure'][0],
                xa_range
            )
        }

    @staticmethod
    def calculate_ideal_vle(vp_a: float, vp_b: float, xa_range: np.ndarray) -> list:
        """計算理想溶液 VLE（Raoult's Law）"""
        results = []
        for xa in xa_range:
            xb = 1 - xa
            P_ideal = xa * vp_a + xb * vp_b
            yA_ideal = (xa * vp_a) / P_ideal if P_ideal > 0 else 0

            results.append({
                'xA': xa,
                'yA_ideal': yA_ideal,
                'P_ideal': P_ideal
            })
        return results

# 使用範例
ensemble = EnsemblePredictor()
ensemble.add_predictor('chemprop_vp', ChemPropVPPredictor('models/chemprop'))
ensemble.add_predictor('cosmo_vle', COSMOVLEPredictor('models/cosmo'))

# 預測完整系統
results = ensemble.predict_full_system(
    cid_a='702',    # Ethanol
    cid_b='1031',   # Benzene
    temperature=298.15
)

# 繪圖
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# xy 圖
xa = [d['xA'] for d in results['vle_data']]
ya = [d['yA'] for d in results['vle_data']]
ya_ideal = [d['yA_ideal'] for d in results['ideal_vle']]

ax1.plot(xa, ya, 'o-', label='Predicted')
ax1.plot(xa, ya_ideal, '--', label='Ideal (Raoult)')
ax1.plot([0, 1], [0, 1], 'k:', label='Diagonal')
ax1.set_xlabel('xA (liquid)')
ax1.set_ylabel('yA (vapor)')
ax1.legend()
ax1.grid(True)

# Px 圖
p = [d['P'] for d in results['vle_data']]
p_ideal = [d['P_ideal'] for d in results['ideal_vle']]

ax2.plot(xa, p, 'o-', label='Predicted')
ax2.plot(xa, p_ideal, '--', label='Ideal (Raoult)')
ax2.axhline(results['pure_vapor_pressure']['A'], color='r', linestyle=':', label='P_A*')
ax2.axhline(results['pure_vapor_pressure']['B'], color='b', linestyle=':', label='P_B*')
ax2.set_xlabel('xA')
ax2.set_ylabel('P (Pa)')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.savefig('vle_prediction.png')

問題 2.2：缺乏 API 介面

問題描述： 無法透過程式化方式呼叫模型，只能透過 Notebook 或命令列。

為什麼這是問題：

  1. 整合困難：無法整合到其他應用程式
  2. 批次處理不便：難以自動化大量預測
  3. 遠端呼叫不可能：無法提供網路服務
  4. 使用者體驗差：需要了解內部實現細節

改良方案：

# REST API 實現 (使用 FastAPI)
from fastapi import FastAPI, HTTPException
from pydantic import BaseModel, Field
from typing import List, Optional
import uvicorn

app = FastAPI(
    title="Chemistry AI API",
    description="AI models for chemical property prediction",
    version="1.0.0"
)

# 請求模型
class VaporPressureRequest(BaseModel):
    smiles: str = Field(..., example="CCO")
    temperature: float = Field(..., ge=200, le=600, example=298.15)

class VLERequest(BaseModel):
    cid_a: str = Field(..., example="702")
    cid_b: str = Field(..., example="1031")
    temperature: float = Field(..., ge=200, le=600, example=298.15)
    xa_values: Optional[List[float]] = Field(None, example=[0.1, 0.3, 0.5, 0.7, 0.9])

# 回應模型
class VaporPressureResponse(BaseModel):
    smiles: str
    temperature: float
    vapor_pressure: float
    unit: str = "Pa"

class VLEResponse(BaseModel):
    cid_a: str
    cid_b: str
    temperature: float
    data: List[dict]

# 全域預測器
predictor = None

@app.on_event("startup")
async def startup_event():
    """啟動時載入模型"""
    global predictor
    predictor = EnsemblePredictor()
    predictor.add_predictor('chemprop_vp', ChemPropVPPredictor('models/chemprop'))
    predictor.add_predictor('cosmo_vle', COSMOVLEPredictor('models/cosmo'))
    print("Models loaded successfully")

@app.get("/")
async def root():
    return {"message": "Chemistry AI API", "version": "1.0.0"}

@app.post("/predict/vapor-pressure", response_model=VaporPressureResponse)
async def predict_vapor_pressure(request: VaporPressureRequest):
    """
    預測單一組分蒸氣壓

    - **smiles**: SMILES 字串
    - **temperature**: 溫度 (K)
    """
    try:
        result = predictor.predict_vapor_pressure(
            request.smiles,
            request.temperature
        )

        return VaporPressureResponse(
            smiles=request.smiles,
            temperature=request.temperature,
            vapor_pressure=float(result['vapor_pressure'][0])
        )

    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))

@app.post("/predict/vle", response_model=VLEResponse)
async def predict_vle(request: VLERequest):
    """
    預測二元混合物汽液平衡

    - **cid_a**: 組分 A 的 PubChem CID
    - **cid_b**: 組分 B 的 PubChem CID
    - **temperature**: 溫度 (K)
    - **xa_values**: 液相組成列表（可選，預設為 [0.1, 0.3, 0.5, 0.7, 0.9]）
    """
    try:
        xa_values = request.xa_values or [0.1, 0.3, 0.5, 0.7, 0.9]

        vle_data = []
        for xa in xa_values:
            result = predictor.predict_vle(
                request.cid_a,
                request.cid_b,
                request.temperature,
                xa
            )

            vle_data.append({
                'xA': xa,
                'yA': float(result['vapor_composition_A'][0]),
                'P': float(result['total_pressure'][0])
            })

        return VLEResponse(
            cid_a=request.cid_a,
            cid_b=request.cid_b,
            temperature=request.temperature,
            data=vle_data
        )

    except Exception as e:
        raise HTTPException(status_code=500, detail=str(e))

@app.get("/health")
async def health_check():
    """健康檢查"""
    return {"status": "healthy", "models_loaded": predictor is not None}

# 運行伺服器
if __name__ == "__main__":
    uvicorn.run(app, host="0.0.0.0", port=8000)

客戶端使用範例：

import requests
import json

# API 端點
BASE_URL = "http://localhost:8000"

# 1. 預測蒸氣壓
vp_request = {
    "smiles": "CCO",
    "temperature": 298.15
}

response = requests.post(f"{BASE_URL}/predict/vapor-pressure", json=vp_request)
print(json.dumps(response.json(), indent=2))

# 2. 預測 VLE
vle_request = {
    "cid_a": "702",
    "cid_b": "1031",
    "temperature": 298.15,
    "xa_values": [0.1, 0.3, 0.5, 0.7, 0.9]
}

response = requests.post(f"{BASE_URL}/predict/vle", json=vle_request)
print(json.dumps(response.json(), indent=2))

資料處理問題

問題 3.1：資料驗證不足

問題描述： 輸入資料缺乏驗證，可能導致下游錯誤。

為什麼這是問題：

  1. 垃圾進垃圾出：錯誤的輸入產生錯誤的輸出
  2. 難以除錯：錯誤在後續步驟才顯現
  3. 資料品質無法保證：無法追蹤資料問題的來源

改良方案：

from pydantic import BaseModel, validator, Field
from rdkit import Chem
from typing import Optional
import re

class MoleculeData(BaseModel):
    """分子資料模型，包含自動驗證"""

    cid: str = Field(..., description="PubChem CID")
    smiles: Optional[str] = Field(None, description="SMILES string")
    temperature: Optional[float] = Field(None, ge=200, le=600, description="Temperature in K")

    @validator('cid')
    def validate_cid(cls, v):
        """驗證 CID 格式"""
        if not v.isdigit():
            raise ValueError(f"CID must be numeric, got: {v}")
        if int(v) <= 0:
            raise ValueError(f"CID must be positive, got: {v}")
        return v

    @validator('smiles')
    def validate_smiles(cls, v):
        """驗證 SMILES 有效性"""
        if v is None:
            return v

        mol = Chem.MolFromSmiles(v)
        if mol is None:
            raise ValueError(f"Invalid SMILES: {v}")

        # 檢查是否為合理的有機分子
        num_atoms = mol.GetNumAtoms()
        if num_atoms < 2:
            raise ValueError(f"Molecule too small (< 2 atoms): {v}")
        if num_atoms > 100:
            raise ValueError(f"Molecule too large (> 100 atoms): {v}")

        return v

    @validator('temperature')
    def validate_temperature(cls, v):
        """驗證溫度合理性"""
        if v is not None:
            if v < 200:
                raise ValueError(f"Temperature too low (< 200 K): {v}")
            if v > 600:
                raise ValueError(f"Temperature too high (> 600 K): {v}")
        return v

class SigmaProfile(BaseModel):
    """σ-profile 資料驗證"""

    sigma: List[float] = Field(..., min_items=51, max_items=51)
    area: float = Field(..., gt=0)
    volume: float = Field(..., gt=0)

    @validator('sigma')
    def validate_sigma(cls, v):
        """驗證 σ-profile 合理性"""
        # 檢查是否為正規化
        total = sum(v)
        if not (0.95 <= total <= 1.05):
            raise ValueError(f"σ-profile not normalized: sum={total}")

        # 檢查是否有負值
        if any(x < 0 for x in v):
            raise ValueError("σ-profile contains negative values")

        return v

    @validator('area')
    def validate_area(cls, v):
        """驗證表面積合理性"""
        if v < 50 or v > 1000:  # Ų
            raise ValueError(f"Unusual surface area: {v} Ų")
        return v

def validate_vle_data(df: pd.DataFrame) -> Tuple[pd.DataFrame, List[str]]:
    """
    驗證 VLE 實驗資料

    Returns:
        (cleaned_df, errors)
    """
    errors = []

    # 必要欄位檢查
    required_cols = ['comp.A', 'comp.B', 'T', 'xA', 'yA', 'P']
    missing_cols = set(required_cols) - set(df.columns)
    if missing_cols:
        raise ValueError(f"Missing required columns: {missing_cols}")

    # 移除 NaN
    initial_rows = len(df)
    df = df.dropna(subset=required_cols)
    if len(df) < initial_rows:
        errors.append(f"Removed {initial_rows - len(df)} rows with NaN values")

    # 驗證數值範圍
    invalid_rows = []

    # 溫度範圍
    temp_mask = (df['T'] < 200) | (df['T'] > 600)
    if temp_mask.any():
        invalid_rows.extend(df[temp_mask].index.tolist())
        errors.append(f"Invalid temperature in rows: {df[temp_mask].index.tolist()}")

    # 組成範圍 [0, 1]
    xa_mask = (df['xA'] < 0) | (df['xA'] > 1)
    ya_mask = (df['yA'] < 0) | (df['yA'] > 1)

    if xa_mask.any():
        invalid_rows.extend(df[xa_mask].index.tolist())
        errors.append(f"Invalid xA in rows: {df[xa_mask].index.tolist()}")

    if ya_mask.any():
        invalid_rows.extend(df[ya_mask].index.tolist())
        errors.append(f"Invalid yA in rows: {df[ya_mask].index.tolist()}")

    # 壓力範圍
    p_mask = (df['P'] <= 0) | (df['P'] > 1e7)  # Pa
    if p_mask.any():
        invalid_rows.extend(df[p_mask].index.tolist())
        errors.append(f"Invalid pressure in rows: {df[p_mask].index.tolist()}")

    # 熱力學一致性檢查：yA 應該在 xA 兩側（輕組分）
    # 簡化檢查：對於大多數系統，yA != xA（除了理想溶液在特定點）

    # 移除所有無效行
    invalid_rows = list(set(invalid_rows))
    if invalid_rows:
        df = df.drop(invalid_rows)
        errors.append(f"Removed {len(invalid_rows)} invalid rows total")

    return df, errors

# 使用範例
try:
    mol = MoleculeData(
        cid="702",
        smiles="CCO",
        temperature=298.15
    )
    print(f"Valid molecule: {mol}")
except ValueError as e:
    print(f"Validation error: {e}")

# 驗證 VLE 資料
df = pd.read_csv('vle_data.csv')
df_clean, errors = validate_vle_data(df)

if errors:
    print("Validation errors:")
    for error in errors:
        print(f"  - {error}")

print(f"Clean data shape: {df_clean.shape}")

問題 3.2：資料不平衡

問題描述： 分組資料嚴重不平衡，Non-HB_Non-HB 組遠大於其他組。

位置：

Non-HB_Non-HB: 222,285 bytes (CSV)
HB-A_HB-A: 58,017 bytes
OH_OH: 57,008 bytes
...
COOH_COOH: 5,072 bytes

為什麼這是問題：

  1. 模型偏差：模型會偏向大組，小組預測不準
  2. 評估不公平：大組主導整體指標
  3. 過擬合風險：小組容易過擬合

改良方案：

from sklearn.utils import resample
from imblearn.over_sampling import SMOTE
from imblearn.under_sampling import RandomUnderSampler
from imblearn.pipeline import Pipeline as ImbPipeline

class BalancedGroupTrainer:
    """處理不平衡分組資料的訓練器"""

    def __init__(self, strategy='hybrid'):
        """
        Args:
            strategy: 'oversample', 'undersample', 'hybrid', 'weighted'
        """
        self.strategy = strategy

    def get_group_statistics(self, df: pd.DataFrame, group_col: str = 'sheet_name'):
        """分析分組統計"""
        stats = df.groupby(group_col).size().sort_values(ascending=False)

        print("Group Statistics:")
        print("-" * 50)
        for group, count in stats.items():
            percentage = count / len(df) * 100
            print(f"{group:30s}: {count:6d} ({percentage:5.2f}%)")

        print("-" * 50)
        print(f"Total: {len(df)}")
        print(f"Max/Min ratio: {stats.max() / stats.min():.2f}x")

        return stats

    def balance_by_oversample(self, X, y, groups, target_ratio=0.5):
        """
        過採樣（SMOTE）

        Args:
            target_ratio: 最小組相對於最大組的目標比例
        """
        print("Applying SMOTE oversampling...")

        # 計算每組的目標樣本數
        group_counts = pd.Series(groups).value_counts()
        max_count = group_counts.max()
        target_count = int(max_count * target_ratio)

        X_resampled = []
        y_resampled = []
        groups_resampled = []

        for group_name in group_counts.index:
            mask = np.array(groups) == group_name
            X_group = X[mask]
            y_group = y[mask]

            current_count = len(X_group)

            if current_count < target_count:
                # 需要過採樣
                # SMOTE 需要至少 k+1 個樣本（預設 k=5）
                if current_count >= 6:
                    smote = SMOTE(random_state=42, k_neighbors=min(5, current_count-1))
                    X_group_res, y_group_res = smote.fit_resample(X_group, y_group)

                    # 如果還不夠，隨機複製
                    if len(X_group_res) < target_count:
                        additional = target_count - len(X_group_res)
                        indices = np.random.choice(len(X_group_res), additional, replace=True)
                        X_group_res = np.vstack([X_group_res, X_group_res[indices]])
                        y_group_res = np.vstack([y_group_res, y_group_res[indices]])
                else:
                    # 樣本太少，只能隨機複製
                    X_group_res, y_group_res = resample(
                        X_group, y_group,
                        n_samples=target_count,
                        replace=True,
                        random_state=42
                    )
            else:
                X_group_res = X_group
                y_group_res = y_group

            X_resampled.append(X_group_res)
            y_resampled.append(y_group_res)
            groups_resampled.extend([group_name] * len(X_group_res))

        X_balanced = np.vstack(X_resampled)
        y_balanced = np.vstack(y_resampled)

        print(f"Original: {len(X)} samples")
        print(f"Balanced: {len(X_balanced)} samples")

        return X_balanced, y_balanced, groups_resampled

    def balance_by_undersample(self, X, y, groups, target_ratio=2.0):
        """
        欠採樣

        Args:
            target_ratio: 最大組相對於最小組的目標比例
        """
        print("Applying undersampling...")

        group_counts = pd.Series(groups).value_counts()
        min_count = group_counts.min()
        target_count = int(min_count * target_ratio)

        X_resampled = []
        y_resampled = []
        groups_resampled = []

        for group_name in group_counts.index:
            mask = np.array(groups) == group_name
            X_group = X[mask]
            y_group = y[mask]

            if len(X_group) > target_count:
                X_group, y_group = resample(
                    X_group, y_group,
                    n_samples=target_count,
                    replace=False,
                    random_state=42
                )

            X_resampled.append(X_group)
            y_resampled.append(y_group)
            groups_resampled.extend([group_name] * len(X_group))

        X_balanced = np.vstack(X_resampled)
        y_balanced = np.vstack(y_resampled)

        print(f"Original: {len(X)} samples")
        print(f"Balanced: {len(X_balanced)} samples")

        return X_balanced, y_balanced, groups_resampled

    def get_class_weights(self, groups):
        """計算類別權重（用於加權損失）"""
        group_counts = pd.Series(groups).value_counts()
        total = len(groups)

        weights = {}
        for group, count in group_counts.items():
            # 反比權重
            weights[group] = total / (len(group_counts) * count)

        # 正規化
        max_weight = max(weights.values())
        weights = {k: v/max_weight for k, v in weights.items()}

        return weights

    def train_with_strategy(self, X, y, groups, model):
        """根據策略訓練模型"""

        if self.strategy == 'oversample':
            X_train, y_train, _ = self.balance_by_oversample(X, y, groups)
            model.fit(X_train, y_train)

        elif self.strategy == 'undersample':
            X_train, y_train, _ = self.balance_by_undersample(X, y, groups)
            model.fit(X_train, y_train)

        elif self.strategy == 'hybrid':
            # 先過採樣小組，再欠採樣大組
            X_over, y_over, groups_over = self.balance_by_oversample(X, y, groups, target_ratio=0.3)
            X_train, y_train, _ = self.balance_by_undersample(X_over, y_over, groups_over, target_ratio=3.0)
            model.fit(X_train, y_train)

        elif self.strategy == 'weighted':
            # 使用樣本權重
            weights = self.get_class_weights(groups)
            sample_weights = np.array([weights[g] for g in groups])

            # Keras 模型需要在 fit 中傳入 sample_weight
            if hasattr(model, 'fit'):
                model.fit(X, y, sample_weight=sample_weights)

        else:
            raise ValueError(f"Unknown strategy: {self.strategy}")

        return model

# 使用範例
df = pd.read_csv('train_group_all_V2.csv')

# 分析分組統計
trainer = BalancedGroupTrainer(strategy='hybrid')
trainer.get_group_statistics(df, 'sheet_name')

# 提取特徵和標籤
X = df.iloc[:, :-3].values  # 特徵
y = df.iloc[:, -2:].values  # yA, P
groups = df['sheet_name'].values

# 平衡並訓練
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import Dense

model = Sequential([
    Dense(128, activation='relu', input_dim=X.shape[1]),
    Dense(128, activation='relu'),
    Dense(2)
])
model.compile(optimizer='adam', loss='mse')

model = trainer.train_with_strategy(X, y, groups, model)

問題 3.3：缺乏資料版本控制

問題描述： 無法追蹤資料變更歷史，難以重現實驗。

為什麼這是問題：

  1. 不可重現：無法確定使用的是哪個版本的資料
  2. 實驗混亂：不同版本的資料產生不同結果
  3. 協作困難：多人使用不同版本的資料

改良方案：

# 使用 DVC (Data Version Control)
"""
安裝：
pip install dvc[all]

初始化：
git init
dvc init

追蹤資料：
dvc add data/train_exp_data_F.csv
dvc add data/s-profiles-all/
git add data/train_exp_data_F.csv.dvc data/.gitignore
git commit -m "Add training data"

配置遠端儲存（例如 S3）：
dvc remote add -d storage s3://my-bucket/dvc-storage
dvc push

其他人拉取資料：
git clone <repo>
dvc pull
"""

# 資料版本管理器
import hashlib
import json
from pathlib import Path
from datetime import datetime

class DataVersionManager:
    """簡單的資料版本管理"""

    def __init__(self, data_dir: Path, version_file: Path = None):
        self.data_dir = Path(data_dir)
        self.version_file = version_file or (self.data_dir / 'data_versions.json')
        self.versions = self.load_versions()

    def load_versions(self) -> dict:
        if self.version_file.exists():
            with open(self.version_file, 'r') as f:
                return json.load(f)
        return {}

    def save_versions(self):
        with open(self.version_file, 'w') as f:
            json.dump(self.versions, f, indent=2)

    def compute_hash(self, file_path: Path) -> str:
        """計算檔案的 MD5 hash"""
        md5 = hashlib.md5()
        with open(file_path, 'rb') as f:
            for chunk in iter(lambda: f.read(4096), b""):
                md5.update(chunk)
        return md5.hexdigest()

    def register_file(self, file_path: Path, description: str = ""):
        """註冊新版本的檔案"""
        file_path = Path(file_path)
        rel_path = file_path.relative_to(self.data_dir)

        file_hash = self.compute_hash(file_path)
        file_size = file_path.stat().st_size

        key = str(rel_path)
        if key not in self.versions:
            self.versions[key] = {'history': []}

        # 檢查是否已存在相同 hash
        if self.versions[key]['history'] and self.versions[key]['history'][-1]['hash'] == file_hash:
            print(f"File {key} unchanged (hash matches)")
            return

        version_info = {
            'version': len(self.versions[key]['history']) + 1,
            'hash': file_hash,
            'size': file_size,
            'timestamp': datetime.now().isoformat(),
            'description': description
        }

        self.versions[key]['history'].append(version_info)
        self.versions[key]['current'] = version_info

        self.save_versions()
        print(f"Registered {key} version {version_info['version']}")

    def verify_file(self, file_path: Path) -> bool:
        """驗證檔案是否與記錄的 hash 匹配"""
        file_path = Path(file_path)
        rel_path = str(file_path.relative_to(self.data_dir))

        if rel_path not in self.versions:
            print(f"Warning: {rel_path} not registered")
            return False

        current_hash = self.compute_hash(file_path)
        expected_hash = self.versions[rel_path]['current']['hash']

        if current_hash != expected_hash:
            print(f"ERROR: {rel_path} hash mismatch!")
            print(f"  Expected: {expected_hash}")
            print(f"  Got: {current_hash}")
            return False

        print(f"✓ {rel_path} verified")
        return True

    def get_manifest(self) -> dict:
        """生成當前資料清單"""
        manifest = {
            'timestamp': datetime.now().isoformat(),
            'files': {}
        }

        for file_key, info in self.versions.items():
            manifest['files'][file_key] = info['current']

        return manifest

    def save_manifest(self, output_path: Path):
        """儲存清單檔案"""
        manifest = self.get_manifest()
        with open(output_path, 'w') as f:
            json.dump(manifest, f, indent=2)
        print(f"Manifest saved to {output_path}")

# 使用範例
data_dir = Path("data")
dvm = DataVersionManager(data_dir)

# 註冊資料檔案
dvm.register_file(
    data_dir / "train_exp_data_F.csv",
    description="Training data with COSMO features, version after preprocessing"
)

dvm.register_file(
    data_dir / "train_group_all_V2.csv",
    description="Grouped training data by hydrogen bonding type"
)

# 驗證所有檔案
for file_path in data_dir.glob("*.csv"):
    dvm.verify_file(file_path)

# 生成實驗清單
dvm.save_manifest(Path("manifests") / f"manifest_{datetime.now().strftime('%Y%m%d')}.json")

模型設計問題

問題 4.1：超參數搜索空間可能不夠優化

問題描述： Optuna 搜索空間可能過大或不夠精細，導致搜索效率低。

當前設置（MLP_ln_V1_COSMO.ipynb）：

n_layers = trial.suggest_int('n_layers', 2, 4)
units = trial.suggest_int(f'units_layer_{i}', 64, 512)
learning_rate = trial.suggest_float('learning_rate', 1e-5, 1e-2, log=True)
batch_size = trial.suggest_categorical('batch_size', [32, 64, 128, 256])
epochs = trial.suggest_int('epochs', 500, 2000)

問題：

  1. 搜索空間太大：200 trials 可能不足以覆蓋
  2. 層數固定：只搜索 2-4 層，可能錯過更深/更淺的網絡
  3. 缺乏 Dropout：沒有正則化超參數
  4. 缺乏學習率衰減：可能影響訓練穩定性

改良方案：

from optuna.samplers import TPESampler
from optuna.pruners import MedianPruner
import optuna

def create_advanced_model(trial, input_dim=108, output_dim=2):
    """
    更進階的模型搜索空間
    """
    # 網絡深度
    n_layers = trial.suggest_int('n_layers', 1, 5)

    # 隱藏層架構策略
    architecture_type = trial.suggest_categorical('architecture', [
        'constant',      # 所有層相同大小
        'decreasing',    # 遞減
        'increasing',    # 遞增
        'bottleneck'     # 中間小兩端大
    ])

    # 基礎單元數
    base_units = trial.suggest_int('base_units', 64, 512, step=32)

    # 根據架構類型建立層
    units = []
    if architecture_type == 'constant':
        units = [base_units] * n_layers

    elif architecture_type == 'decreasing':
        decay_rate = trial.suggest_float('decay_rate', 0.5, 0.9)
        for i in range(n_layers):
            units.append(int(base_units * (decay_rate ** i)))

    elif architecture_type == 'increasing':
        growth_rate = trial.suggest_float('growth_rate', 1.1, 2.0)
        for i in range(n_layers):
            units.append(int(base_units * (growth_rate ** i)))

    elif architecture_type == 'bottleneck':
        bottleneck_ratio = trial.suggest_float('bottleneck_ratio', 0.2, 0.5)
        mid = n_layers // 2
        for i in range(n_layers):
            if i == mid:
                units.append(int(base_units * bottleneck_ratio))
            else:
                units.append(base_units)

    # 正則化
    dropout_rate = trial.suggest_float('dropout', 0.0, 0.5)
    use_batch_norm = trial.suggest_categorical('batch_norm', [True, False])
    l2_reg = trial.suggest_float('l2_reg', 1e-6, 1e-3, log=True)

    # 激活函數
    activation = trial.suggest_categorical('activation', ['relu', 'elu', 'selu', 'tanh'])

    # 建立模型
    from tensorflow.keras.layers import Dense, Dropout, BatchNormalization
    from tensorflow.keras.models import Sequential
    from tensorflow.keras.regularizers import l2

    model = Sequential()

    # 輸入層
    model.add(Dense(
        units[0],
        activation=activation,
        input_dim=input_dim,
        kernel_regularizer=l2(l2_reg)
    ))

    if use_batch_norm:
        model.add(BatchNormalization())

    if dropout_rate > 0:
        model.add(Dropout(dropout_rate))

    # 隱藏層
    for i in range(1, n_layers):
        model.add(Dense(
            units[i],
            activation=activation,
            kernel_regularizer=l2(l2_reg)
        ))

        if use_batch_norm:
            model.add(BatchNormalization())

        if dropout_rate > 0:
            model.add(Dropout(dropout_rate))

    # 輸出層
    model.add(Dense(output_dim))

    # 優化器配置
    optimizer_name = trial.suggest_categorical('optimizer', ['adam', 'adamw', 'rmsprop'])
    learning_rate = trial.suggest_float('learning_rate', 1e-5, 1e-2, log=True)

    if optimizer_name == 'adam':
        optimizer = tf.keras.optimizers.Adam(learning_rate)
    elif optimizer_name == 'adamw':
        optimizer = tf.keras.optimizers.AdamW(learning_rate)
    elif optimizer_name == 'rmsprop':
        optimizer = tf.keras.optimizers.RMSprop(learning_rate)

    model.compile(optimizer=optimizer, loss='mse', metrics=['mae'])

    return model

def objective(trial, X_train, y_train, X_val, y_val):
    """
    優化目標函數，加入早停和學習率衰減
    """
    # 建立模型
    model = create_advanced_model(trial, input_dim=X_train.shape[1])

    # 訓練配置
    batch_size = trial.suggest_categorical('batch_size', [32, 64, 128, 256])
    max_epochs = 2000

    # Callbacks
    callbacks = []

    # 1. Early Stopping
    early_stop = tf.keras.callbacks.EarlyStopping(
        monitor='val_loss',
        patience=100,
        restore_best_weights=True
    )
    callbacks.append(early_stop)

    # 2. Learning Rate Scheduling
    use_lr_schedule = trial.suggest_categorical('lr_schedule', [True, False])
    if use_lr_schedule:
        schedule_type = trial.suggest_categorical('schedule_type', [
            'step',
            'exponential',
            'cosine'
        ])

        if schedule_type == 'step':
            step_decay = tf.keras.callbacks.LearningRateScheduler(
                lambda epoch: trial.params['learning_rate'] * (0.95 ** (epoch // 50))
            )
            callbacks.append(step_decay)

        elif schedule_type == 'exponential':
            exp_decay = tf.keras.callbacks.LearningRateScheduler(
                lambda epoch: trial.params['learning_rate'] * np.exp(-0.01 * epoch)
            )
            callbacks.append(exp_decay)

        elif schedule_type == 'cosine':
            cosine_decay = tf.keras.optimizers.schedules.CosineDecay(
                trial.params['learning_rate'],
                decay_steps=max_epochs
            )
            # 需要重新編譯模型
            model.optimizer.learning_rate = cosine_decay

    # 3. Optuna Pruning Callback
    pruning_callback = optuna.integration.TFKerasPruningCallback(trial, 'val_loss')
    callbacks.append(pruning_callback)

    # 訓練
    history = model.fit(
        X_train, y_train,
        validation_data=(X_val, y_val),
        epochs=max_epochs,
        batch_size=batch_size,
        callbacks=callbacks,
        verbose=0
    )

    # 評估
    y_pred = model.predict(X_val)

    # 自定義評估指標（AAD-y + AARD-P）
    # 需要反正規化
    y_val_denorm = denormalize(y_val, scaler)
    y_pred_denorm = denormalize(y_pred, scaler)

    aad_y = np.mean(np.abs(y_val_denorm[:, 0] - y_pred_denorm[:, 0]))
    aard_p = np.mean(np.abs((y_val_denorm[:, 1] - y_pred_denorm[:, 1]) / y_val_denorm[:, 1])) * 100

    # 組合指標
    score = aad_y * 100 + aard_p  # AAD-y(%) + AARD-P(%)

    return score

# 執行優化
study = optuna.create_study(
    direction='minimize',
    sampler=TPESampler(seed=42),
    pruner=MedianPruner(n_startup_trials=10, n_warmup_steps=5)
)

study.optimize(
    lambda trial: objective(trial, X_train, y_train, X_val, y_val),
    n_trials=300,  # 增加 trials
    timeout=86400,  # 24 小時
    show_progress_bar=True
)

# 分析結果
print("Best trial:")
print(f"  Value: {study.best_value:.4f}")
print("  Params:")
for key, value in study.best_params.items():
    print(f"    {key}: {value}")

# 視覺化
import optuna.visualization as vis

# 優化歷史
fig = vis.plot_optimization_history(study)
fig.write_html("optimization_history.html")

# 參數重要性
fig = vis.plot_param_importances(study)
fig.write_html("param_importances.html")

# 平行座標圖
fig = vis.plot_parallel_coordinate(study)
fig.write_html("parallel_coordinate.html")

問題 4.2：缺乏集成學習（Ensemble）

問題描述： COSMO-MLP 訓練了 10-fold 模型但未使用集成預測。

為什麼這是問題：

  1. 浪費資源：訓練了 10 個模型但只用一個
  2. 預測不穩定：單一模型容易受到訓練隨機性影響
  3. 錯失改進機會：集成通常能提升 1-3% 準確度

改良方案：

from typing import List
import numpy as np

class EnsembleVLEPredictor:
    """集成 VLE 預測器"""

    def __init__(self, model_dir: str, n_folds: int = 10):
        self.models = []
        self.scalers = []
        self.n_folds = n_folds

        # 載入所有 fold 模型
        for fold in range(1, n_folds + 1):
            model = load_model(f"{model_dir}/MLP_Vali_COSMO_{fold}.h5", compile=False)
            scaler = load(f"{model_dir}/pipelineMLP_COSMO_fold{fold}.joblib")

            self.models.append(model)
            self.scalers.append(scaler)

        print(f"Loaded {len(self.models)} models")

    def predict_mean(self, X: np.ndarray) -> np.ndarray:
        """平均集成"""
        predictions = []

        for model, scaler in zip(self.models, self.scalers):
            X_scaled = scaler.transform(X)
            y_pred = model.predict(X_scaled, verbose=0)
            predictions.append(y_pred)

        return np.mean(predictions, axis=0)

    def predict_median(self, X: np.ndarray) -> np.ndarray:
        """中位數集成（對異常值更魯棒）"""
        predictions = []

        for model, scaler in zip(self.models, self.scalers):
            X_scaled = scaler.transform(X)
            y_pred = model.predict(X_scaled, verbose=0)
            predictions.append(y_pred)

        return np.median(predictions, axis=0)

    def predict_weighted(self, X: np.ndarray, weights: List[float] = None) -> np.ndarray:
        """加權集成"""
        if weights is None:
            # 根據驗證性能計算權重
            weights = self.calculate_weights()

        weights = np.array(weights)
        weights = weights / weights.sum()  # 正規化

        predictions = []

        for model, scaler in zip(self.models, self.scalers):
            X_scaled = scaler.transform(X)
            y_pred = model.predict(X_scaled, verbose=0)
            predictions.append(y_pred)

        predictions = np.array(predictions)
        return np.average(predictions, axis=0, weights=weights)

    def predict_with_uncertainty(self, X: np.ndarray) -> tuple:
        """預測並返回不確定性估計"""
        predictions = []

        for model, scaler in zip(self.models, self.scalers):
            X_scaled = scaler.transform(X)
            y_pred = model.predict(X_scaled, verbose=0)
            predictions.append(y_pred)

        predictions = np.array(predictions)

        mean = np.mean(predictions, axis=0)
        std = np.std(predictions, axis=0)

        return mean, std

    def calculate_weights(self) -> np.ndarray:
        """
        根據驗證集性能計算每個模型的權重
        需要事先載入各 fold 的驗證結果
        """
        # 這裡簡化為等權重，實際應根據驗證性能
        # 例如：weights = 1 / val_losses
        return np.ones(self.n_folds)

    def predict_bootstrap(self, X: np.ndarray, n_bootstrap: int = 100) -> tuple:
        """
        Bootstrap 集成
        隨機抽樣模型進行預測，提供更穩健的不確定性估計
        """
        bootstrap_predictions = []

        for _ in range(n_bootstrap):
            # 隨機選擇模型子集
            indices = np.random.choice(self.n_folds, size=self.n_folds, replace=True)

            fold_predictions = []
            for idx in indices:
                model = self.models[idx]
                scaler = self.scalers[idx]

                X_scaled = scaler.transform(X)
                y_pred = model.predict(X_scaled, verbose=0)
                fold_predictions.append(y_pred)

            bootstrap_predictions.append(np.mean(fold_predictions, axis=0))

        bootstrap_predictions = np.array(bootstrap_predictions)

        mean = np.mean(bootstrap_predictions, axis=0)
        std = np.std(bootstrap_predictions, axis=0)

        # 95% 信賴區間
        lower = np.percentile(bootstrap_predictions, 2.5, axis=0)
        upper = np.percentile(bootstrap_predictions, 97.5, axis=0)

        return mean, std, lower, upper

# 使用範例
ensemble = EnsembleVLEPredictor('results_V1/model/h5', n_folds=10)

# 1. 基本預測（平均）
y_pred_mean = ensemble.predict_mean(X_test)

# 2. 預測與不確定性
y_pred, y_std = ensemble.predict_with_uncertainty(X_test)

print(f"Predicted yA: {y_pred[0, 0]:.4f} ± {y_std[0, 0]:.4f}")
print(f"Predicted P: {y_pred[0, 1]:.2f} ± {y_std[0, 1]:.2f} Pa")

# 3. Bootstrap 集成（最穩健）
y_pred, y_std, y_lower, y_upper = ensemble.predict_bootstrap(X_test, n_bootstrap=100)

print(f"Predicted yA: {y_pred[0, 0]:.4f} [{y_lower[0, 0]:.4f}, {y_upper[0, 0]:.4f}]")
print(f"Predicted P: {y_pred[0, 1]:.2f} [{y_lower[0, 1]:.2f}, {y_upper[0, 1]:.2f}] Pa")

# 4. 視覺化不確定性
import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# yA
axes[0].errorbar(
    range(len(y_pred)),
    y_pred[:, 0],
    yerr=y_std[:, 0],
    fmt='o',
    capsize=5,
    label='Predicted ± std'
)
axes[0].plot(y_test[:, 0], 's', label='Actual', alpha=0.5)
axes[0].set_xlabel('Sample')
axes[0].set_ylabel('yA')
axes[0].legend()
axes[0].grid(True)

# P
axes[1].errorbar(
    range(len(y_pred)),
    y_pred[:, 1],
    yerr=y_std[:, 1],
    fmt='o',
    capsize=5,
    label='Predicted ± std'
)
axes[1].plot(y_test[:, 1], 's', label='Actual', alpha=0.5)
axes[1].set_xlabel('Sample')
axes[1].set_ylabel('P (Pa)')
axes[1].legend()
axes[1].grid(True)

plt.tight_layout()
plt.savefig('ensemble_uncertainty.png', dpi=300)

訓練流程問題

問題 5.1：缺乏實驗追蹤

問題描述： 無法系統化追蹤實驗參數、結果和比較。

為什麼這是問題：

  1. 無法比較：不清楚哪個實驗最好
  2. 無法重現：不知道確切的參數設置
  3. 浪費時間：可能重複相同實驗

改良方案：

# 使用 MLflow 進行實驗追蹤
import mlflow
import mlflow.tensorflow

# 設置追蹤 URI
mlflow.set_tracking_uri("file:./mlruns")
mlflow.set_experiment("COSMO_VLE_Prediction")

def train_with_mlflow(X_train, y_train, X_val, y_val, params):
    """
    使用 MLflow 追蹤訓練過程
    """
    with mlflow.start_run(run_name=f"fold_{params['fold']}"):
        # 記錄參數
        mlflow.log_params(params)

        # 記錄資料集資訊
        mlflow.log_param("train_size", len(X_train))
        mlflow.log_param("val_size", len(X_val))
        mlflow.log_param("n_features", X_train.shape[1])

        # 建立模型
        model = build_model(params)

        # 自定義 Callback 記錄指標
        class MLflowCallback(tf.keras.callbacks.Callback):
            def on_epoch_end(self, epoch, logs=None):
                mlflow.log_metrics({
                    "train_loss": logs['loss'],
                    "val_loss": logs['val_loss'],
                    "train_mae": logs['mae'],
                    "val_mae": logs['val_mae']
                }, step=epoch)

        # 訓練
        history = model.fit(
            X_train, y_train,
            validation_data=(X_val, y_val),
            epochs=params['epochs'],
            batch_size=params['batch_size'],
            callbacks=[
                tf.keras.callbacks.EarlyStopping(patience=100),
                MLflowCallback()
            ],
            verbose=1
        )

        # 評估
        y_pred = model.predict(X_val)

        # 計算指標
        aad_y = np.mean(np.abs(y_val[:, 0] - y_pred[:, 0]))
        aard_p = np.mean(np.abs((y_val[:, 1] - y_pred[:, 1]) / y_val[:, 1])) * 100

        # 記錄最終指標
        mlflow.log_metrics({
            "final_val_loss": min(history.history['val_loss']),
            "AAD_y": aad_y,
            "AARD_P": aard_p,
            "combined_score": aad_y * 100 + aard_p
        })

        # 儲存模型
        mlflow.tensorflow.log_model(model, "model")

        # 儲存圖表
        fig, axes = plt.subplots(1, 2, figsize=(12, 5))

        axes[0].plot(history.history['loss'], label='Train')
        axes[0].plot(history.history['val_loss'], label='Val')
        axes[0].set_title('Loss')
        axes[0].legend()

        axes[1].scatter(y_val[:, 0], y_pred[:, 0], alpha=0.5)
        axes[1].plot([0, 1], [0, 1], 'r--')
        axes[1].set_title('yA: Predicted vs Actual')

        plt.tight_layout()
        mlflow.log_figure(fig, "training_results.png")
        plt.close()

        # 記錄配置檔
        with open('model_config.json', 'w') as f:
            json.dump(params, f, indent=2)
        mlflow.log_artifact('model_config.json')

        return model, history

# 使用範例
kf = KFold(n_splits=10, shuffle=True, random_state=42)

for fold, (train_idx, val_idx) in enumerate(kf.split(X)):
    params = {
        'fold': fold + 1,
        'n_layers': 2,
        'units_layer_1': 129,
        'units_layer_2': 236,
        'learning_rate': 0.000252,
        'batch_size': 64,
        'epochs': 1500
    }

    X_train, X_val = X[train_idx], X[val_idx]
    y_train, y_val = y[train_idx], y[val_idx]

    model, history = train_with_mlflow(X_train, y_train, X_val, y_val, params)

# 比較實驗
from mlflow import search_runs

# 獲取所有實驗
runs = search_runs(experiment_names=["COSMO_VLE_Prediction"])

# 排序找出最佳
best_runs = runs.sort_values('metrics.combined_score').head(5)

print("Top 5 runs:")
print(best_runs[['run_id', 'metrics.AAD_y', 'metrics.AARD_P', 'metrics.combined_score']])

# 載入最佳模型
best_run_id = best_runs.iloc[0]['run_id']
best_model = mlflow.tensorflow.load_model(f"runs:/{best_run_id}/model")

問題 5.2：訓練時間過長

問題描述： 每個 fold 訓練 4-12 小時，10-fold 總共需要數天。

為什麼這是問題：

  1. 開發迭代慢：測試新想法需要等很久
  2. 資源浪費：可能早停但設置了太多 epochs
  3. 實驗次數受限：無法進行大量實驗

改良方案：

# 1. 使用學習率預熱和早停
def create_model_with_lr_warmup(params):
    """帶學習率預熱的模型"""
    model = build_model(params)

    # 自定義學習率調度
    class WarmUpSchedule(tf.keras.optimizers.schedules.LearningRateSchedule):
        def __init__(self, initial_learning_rate, warmup_steps):
            super().__init__()
            self.initial_learning_rate = initial_learning_rate
            self.warmup_steps = warmup_steps

        def __call__(self, step):
            arg1 = tf.math.rsqrt(step + 1)
            arg2 = step * (self.warmup_steps ** -1.5)
            return self.initial_learning_rate * tf.math.minimum(arg1, arg2)

    lr_schedule = WarmUpSchedule(params['learning_rate'], warmup_steps=1000)

    optimizer = tf.keras.optimizers.Adam(lr_schedule)
    model.compile(optimizer=optimizer, loss='mse', metrics=['mae'])

    return model

# 2. 使用更智能的早停策略
class AdaptiveEarlyStopping(tf.keras.callbacks.Callback):
    """自適應早停 - 根據改善速度動態調整耐心"""

    def __init__(self, monitor='val_loss', min_delta=0.0001, 
                 initial_patience=50, max_patience=200):
        super().__init__()
        self.monitor = monitor
        self.min_delta = min_delta
        self.initial_patience = initial_patience
        self.max_patience = max_patience
        self.patience = initial_patience
        self.wait = 0
        self.best = np.inf
        self.best_weights = None

    def on_epoch_end(self, epoch, logs=None):
        current = logs.get(self.monitor)

        if current < self.best - self.min_delta:
            # 顯著改善
            self.best = current
            self.wait = 0
            self.best_weights = self.model.get_weights()

            # 重置耐心
            self.patience = self.initial_patience
        else:
            # 沒有改善
            self.wait += 1

            # 增加耐心（可能在平台期）
            if self.wait > self.patience * 0.8:
                self.patience = min(self.patience * 1.2, self.max_patience)

            if self.wait >= self.patience:
                print(f"\nEarly stopping at epoch {epoch+1}")
                self.model.stop_training = True
                self.model.set_weights(self.best_weights)

# 3. 使用混合精度訓練（加速 2x）
from tensorflow.keras import mixed_precision

# 啟用混合精度
policy = mixed_precision.Policy('mixed_float16')
mixed_precision.set_global_policy(policy)

# 建立模型時輸出層保持 float32
model = Sequential([
    Dense(128, activation='relu', input_dim=108),
    Dense(128, activation='relu'),
    Dense(2, dtype='float32')  # 輸出層用 float32
])

# 4. 使用資料預處理管道（減少 I/O 時間）
def create_tf_dataset(X, y, batch_size=64, shuffle=True):
    """建立高效的 TensorFlow Dataset"""
    dataset = tf.data.Dataset.from_tensor_slices((X, y))

    if shuffle:
        dataset = dataset.shuffle(buffer_size=1000)

    dataset = dataset.batch(batch_size)
    dataset = dataset.prefetch(tf.data.AUTOTUNE)

    return dataset

# 使用範例
train_dataset = create_tf_dataset(X_train, y_train, batch_size=64)
val_dataset = create_tf_dataset(X_val, y_val, batch_size=64, shuffle=False)

# 5. 平行訓練多個 fold（如果有多張 GPU）
import multiprocessing as mp

def train_single_fold(fold_idx, X_train, y_train, X_val, y_val, gpu_id):
    """在指定 GPU 上訓練單個 fold"""
    # 設置使用的 GPU
    os.environ['CUDA_VISIBLE_DEVICES'] = str(gpu_id)

    # 訓練邏輯
    model = create_model_with_lr_warmup(params)

    history = model.fit(
        X_train, y_train,
        validation_data=(X_val, y_val),
        epochs=2000,
        batch_size=64,
        callbacks=[
            AdaptiveEarlyStopping(initial_patience=50, max_patience=200)
        ],
        verbose=1
    )

    model.save(f'models/fold_{fold_idx+1}.h5')

    return fold_idx, history

# 平行訓練（假設有 4 張 GPU）
from concurrent.futures import ProcessPoolExecutor

n_gpus = 4
kf = KFold(n_splits=10, shuffle=True, random_state=42)

with ProcessPoolExecutor(max_workers=n_gpus) as executor:
    futures = []

    for fold_idx, (train_idx, val_idx) in enumerate(kf.split(X)):
        gpu_id = fold_idx % n_gpus

        X_train, X_val = X[train_idx], X[val_idx]
        y_train, y_val = y[train_idx], y[val_idx]

        future = executor.submit(
            train_single_fold,
            fold_idx, X_train, y_train, X_val, y_val, gpu_id
        )
        futures.append(future)

    # 等待所有 fold 完成
    results = [f.result() for f in futures]

print(f"Training complete: {len(results)} folds trained")

可擴展性問題

問題 6.1：無法處理新分子

問題描述： 如果要預測新分子，需要完整跑一遍 COSMO 計算流程（耗時數小時）。

為什麼這是問題：

  1. 使用者體驗差：無法即時預測
  2. 計算成本高：需要 Gaussian 軟體和 HPC
  3. 限制應用場景：無法用於快速篩選

改良方案：

# 方案 1：預訓練 σ-profile 生成模型
from tensorflow.keras.layers import Input, Dense, LSTM, Bidirectional
from tensorflow.keras.models import Model

class SigmaProfileGenerator:
    """
    從 SMILES 直接生成 σ-profile 的神經網絡
    跳過 DFT 計算
    """

    def __init__(self):
        self.model = self.build_model()
        self.char_to_idx = self.build_vocabulary()

    def build_vocabulary(self):
        """SMILES 字符映射"""
        chars = ['C', 'N', 'O', 'S', 'F', 'Cl', 'Br', 'I', 
                '(', ')', '[', ']', '=', '#', '@', '+', '-', 
                '1', '2', '3', '4', '5', '6', '7', '8', '9', 'c', 'n', 'o', 's']
        return {c: i+1 for i, c in enumerate(chars)}  # 0 留給 padding

    def build_model(self):
        """建立 SMILES -> σ-profile 模型"""
        # 編碼器（SMILES 序列）
        smiles_input = Input(shape=(None,), name='smiles')
        x = tf.keras.layers.Embedding(len(self.char_to_idx) + 1, 128)(smiles_input)
        x = Bidirectional(LSTM(256, return_sequences=True))(x)
        x = Bidirectional(LSTM(256))(x)

        # 分子級特徵預測
        molecular_features = Dense(128, activation='relu')(x)

        # σ-profile 預測（51 bins）
        sigma_profile = Dense(51, activation='softmax', name='sigma_profile')(molecular_features)

        # 面積和體積預測
        area = Dense(1, activation='relu', name='area')(molecular_features)
        volume = Dense(1, activation='relu', name='volume')(molecular_features)

        model = Model(
            inputs=smiles_input,
            outputs=[sigma_profile, area, volume]
        )

        model.compile(
            optimizer='adam',
            loss={
                'sigma_profile': 'mse',
                'area': 'mse',
                'volume': 'mse'
            },
            loss_weights={
                'sigma_profile': 1.0,
                'area': 0.1,
                'volume': 0.1
            }
        )

        return model

    def smiles_to_sequence(self, smiles: str, max_length: int = 100) -> np.ndarray:
        """將 SMILES 轉換為序列"""
        seq = [self.char_to_idx.get(c, 0) for c in smiles]

        # Padding
        if len(seq) < max_length:
            seq = seq + [0] * (max_length - len(seq))
        else:
            seq = seq[:max_length]

        return np.array([seq])

    def predict_sigma_profile(self, smiles: str) -> dict:
        """從 SMILES 預測 σ-profile"""
        seq = self.smiles_to_sequence(smiles)

        sigma, area, volume = self.model.predict(seq, verbose=0)

        return {
            'sigma': sigma[0],
            'area': area[0, 0],
            'volume': volume[0, 0]
        }

    def train_from_database(self, smiles_list: List[str], 
                           sigma_profiles: List[np.ndarray],
                           areas: List[float],
                           volumes: List[float]):
        """
        從現有 COSMO 資料庫訓練模型

        Args:
            smiles_list: SMILES 列表
            sigma_profiles: 對應的 σ-profile
            areas: 對應的表面積
            volumes: 對應的體積
        """
        # 準備輸入
        X = np.array([self.smiles_to_sequence(s) for s in smiles_list]).squeeze(1)

        # 準備輸出
        y = {
            'sigma_profile': np.array(sigma_profiles),
            'area': np.array(areas).reshape(-1, 1),
            'volume': np.array(volumes).reshape(-1, 1)
        }

        # 訓練
        history = self.model.fit(
            X, y,
            epochs=100,
            batch_size=32,
            validation_split=0.2,
            callbacks=[
                tf.keras.callbacks.EarlyStopping(patience=10),
                tf.keras.callbacks.ReduceLROnPlateau(patience=5)
            ]
        )

        return history

# 使用範例
generator = SigmaProfileGenerator()

# 1. 從現有資料庫訓練
# 載入所有 COSMO 資料
cosmo_db = load_cosmo_database()  # 假設此函數存在

generator.train_from_database(
    smiles_list=cosmo_db['smiles'],
    sigma_profiles=cosmo_db['sigma_profiles'],
    areas=cosmo_db['areas'],
    volumes=cosmo_db['volumes']
)

# 2. 儲存模型
generator.model.save('sigma_profile_generator.h5')

# 3. 快速預測新分子
new_smiles = "CCO"
sigma_info = generator.predict_sigma_profile(new_smiles)

print(f"σ-profile: {sigma_info['sigma']}")
print(f"Area: {sigma_info['area']:.2f} Ų")
print(f"Volume: {sigma_info['volume']:.2f} ų")

# 方案 2：使用遷移學習（Graph Neural Network）
from torch_geometric.nn import GCNConv, global_mean_pool
import torch
import torch.nn as nn

class GNN_SigmaProfile(nn.Module):
    """
    使用 GNN 從分子圖直接預測 σ-profile
    """

    def __init__(self, num_node_features=9, hidden_channels=64):
        super().__init__()
        self.conv1 = GCNConv(num_node_features, hidden_channels)
        self.conv2 = GCNConv(hidden_channels, hidden_channels)
        self.conv3 = GCNConv(hidden_channels, hidden_channels)

        # σ-profile 預測
        self.sigma_head = nn.Sequential(
            nn.Linear(hidden_channels, 128),
            nn.ReLU(),
            nn.Linear(128, 51),
            nn.Softmax(dim=1)
        )

        # 分子性質預測
        self.property_head = nn.Sequential(
            nn.Linear(hidden_channels, 64),
            nn.ReLU(),
            nn.Linear(64, 2)  # area, volume
        )

    def forward(self, x, edge_index, batch):
        # GNN 編碼
        x = self.conv1(x, edge_index)
        x = F.relu(x)
        x = self.conv2(x, edge_index)
        x = F.relu(x)
        x = self.conv3(x, edge_index)

        # 圖級聚合
        x = global_mean_pool(x, batch)

        # 預測
        sigma = self.sigma_head(x)
        properties = self.property_head(x)

        return sigma, properties

# 使用 GNN 進行預測
model = GNN_SigmaProfile()
# ... 訓練過程 ...

# 快速預測
from rdkit import Chem
from torch_geometric.data import Data

def smiles_to_graph(smiles):
    """將 SMILES 轉換為圖"""
    mol = Chem.MolFromSmiles(smiles)

    # 節點特徵（原子）
    node_features = []
    for atom in mol.GetAtoms():
        features = [
            atom.GetAtomicNum(),
            atom.GetDegree(),
            atom.GetFormalCharge(),
            atom.GetHybridization(),
            atom.GetIsAromatic(),
            # ... 更多特徵
        ]
        node_features.append(features)

    # 邊（鍵）
    edge_index = []
    for bond in mol.GetBonds():
        edge_index.append([bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()])
        edge_index.append([bond.GetEndAtomIdx(), bond.GetBeginAtomIdx()])

    return Data(
        x=torch.tensor(node_features, dtype=torch.float),
        edge_index=torch.tensor(edge_index, dtype=torch.long).t()
    )

# 預測新分子
graph = smiles_to_graph("CCO")
sigma, properties = model(graph.x, graph.edge_index, torch.zeros(len(graph.x), dtype=torch.long))

print(f"Predicted σ-profile: {sigma.detach().numpy()}")
print(f"Predicted area: {properties[0, 0].item():.2f}")
print(f"Predicted volume: {properties[0, 1].item():.2f}")

綜合改良建議

優先級 1（高）- 立即實施

  1. 配置管理
  • 使用 config.yaml 替換硬編碼路徑

  • 實施 pathlib 統一路徑處理 2. 錯誤處理

  • 為所有外部 API 呼叫（PubChem、Gaussian）添加 try-except

  • 實施日誌系統 3. 模組化重構

  • 將 Notebook 中的核心邏輯提取為 Python 模組

  • 建立可重用的 API

優先級 2（中）- 短期改進

  1. 集成學習
  • 實施 10-fold 模型的集成預測

  • 提供不確定性估計 2. 資料驗證

  • 使用 Pydantic 驗證輸入

  • 實施資料品質檢查 3. 實驗追蹤

  • 整合 MLflow

  • 記錄所有超參數和結果

優先級 3（低）- 長期優化

  1. 快速預測模型
  • 訓練 SMILES → σ-profile 生成模型

  • 跳過耗時的 DFT 計算 2. Web API

  • 建立 FastAPI 服務

  • 提供 RESTful 介面 3. 進階模型

  • 嘗試 Transformer 架構

  • 探索 Graph Neural Networks

實施路線圖

第 1 週：
- 建立 config.yaml
- 實施日誌系統
- 添加錯誤處理

第 2-3 週：
- 模組化重構
- 建立 src/ 目錄結構
- 撰寫單元測試

第 4 週：
- 實施集成預測
- 整合 MLflow

第 5-6 週：
- 建立 FastAPI 服務
- 編寫 API 文檔

第 7-8 週：
- 訓練快速預測模型
- 優化推理速度

長期：
- 持續優化模型架構
- 擴展到更多化學系統

文件版本: 1.0 更新日期: 2026-02-12