目錄
程式碼品質問題
問題 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. 使用相對路徑和 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") # 檔案不存在會崩潰
為什麼這是問題：
- 容錯性差：一個錯誤會導致整個批次處理失敗
- 除錯困難：錯誤訊息不明確，難以定位問題
- 資料損失：處理中途失敗可能導致已處理資料遺失
- 使用者體驗差：無法友善地告知使用者問題所在
改良方案：
# 改進的 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 中，難以重用和測試。
為什麼這是問題：
- 重用性差：相同功能需要複製貼上
- 測試困難：Notebook 難以進行單元測試
- 版本控制不友善：Git 難以追蹤 Notebook 的變更
- 效能問題：每次都要執行所有 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 是兩個完全獨立的專案，無法互相利用。
為什麼這是問題：
- 資料重複：兩個專案可能處理相同的分子
- 功能重複：分子處理邏輯重複實現
- 無法協同：無法結合兩個模型的優勢
- 維護成本高：需要分別維護兩套系統
改良方案：
# 統一的架構設計
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 或命令列。
為什麼這是問題：
- 整合困難：無法整合到其他應用程式
- 批次處理不便：難以自動化大量預測
- 遠端呼叫不可能：無法提供網路服務
- 使用者體驗差：需要了解內部實現細節
改良方案：
# 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：資料驗證不足
問題描述： 輸入資料缺乏驗證，可能導致下游錯誤。
為什麼這是問題：
- 垃圾進垃圾出：錯誤的輸入產生錯誤的輸出
- 難以除錯：錯誤在後續步驟才顯現
- 資料品質無法保證：無法追蹤資料問題的來源
改良方案：
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
為什麼這是問題：
- 模型偏差：模型會偏向大組，小組預測不準
- 評估不公平：大組主導整體指標
- 過擬合風險：小組容易過擬合
改良方案：
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：缺乏資料版本控制
問題描述： 無法追蹤資料變更歷史，難以重現實驗。
為什麼這是問題：
- 不可重現：無法確定使用的是哪個版本的資料
- 實驗混亂：不同版本的資料產生不同結果
- 協作困難：多人使用不同版本的資料
改良方案：
# 使用 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)
問題：
- 搜索空間太大：200 trials 可能不足以覆蓋
- 層數固定：只搜索 2-4 層，可能錯過更深/更淺的網絡
- 缺乏 Dropout：沒有正則化超參數
- 缺乏學習率衰減：可能影響訓練穩定性
改良方案：
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 模型但未使用集成預測。
為什麼這是問題：
- 浪費資源：訓練了 10 個模型但只用一個
- 預測不穩定：單一模型容易受到訓練隨機性影響
- 錯失改進機會：集成通常能提升 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：缺乏實驗追蹤
問題描述： 無法系統化追蹤實驗參數、結果和比較。
為什麼這是問題：
- 無法比較：不清楚哪個實驗最好
- 無法重現：不知道確切的參數設置
- 浪費時間：可能重複相同實驗
改良方案：
# 使用 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 總共需要數天。
為什麼這是問題：
- 開發迭代慢：測試新想法需要等很久
- 資源浪費：可能早停但設置了太多 epochs
- 實驗次數受限：無法進行大量實驗
改良方案：
# 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 計算流程（耗時數小時）。
為什麼這是問題：
- 使用者體驗差：無法即時預測
- 計算成本高：需要 Gaussian 軟體和 HPC
- 限制應用場景：無法用於快速篩選
改良方案：
# 方案 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（高）- 立即實施
- 配置管理
- 使用
config.yaml替換硬編碼路徑 -
實施 pathlib 統一路徑處理 2. 錯誤處理
-
為所有外部 API 呼叫（PubChem、Gaussian）添加 try-except
-
實施日誌系統 3. 模組化重構
-
將 Notebook 中的核心邏輯提取為 Python 模組
- 建立可重用的 API
優先級 2（中）- 短期改進
- 集成學習
- 實施 10-fold 模型的集成預測
-
提供不確定性估計 2. 資料驗證
-
使用 Pydantic 驗證輸入
-
實施資料品質檢查 3. 實驗追蹤
-
整合 MLflow
- 記錄所有超參數和結果
優先級 3（低）- 長期優化
- 快速預測模型
- 訓練 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