{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "34d16d4a",
   "metadata": {},
   "source": [
    "# HPLC 分離最佳化教學筆記本\n",
    "\n",
    "本筆記本以互動、可修改的方式，帶你重現以下論文所提出的 HPLC（高效液相層析）\n",
    "分離最佳化方法：\n",
    "\n",
    "> Zisi, Ch., Pappa-Louisi, A., & Nikitas, P. (2020). Separation optimization in\n",
    "> HPLC analysis implemented in R programming language. *Journal of\n",
    "> Chromatography A*, 1617, 460823.\n",
    "\n",
    "**版權聲明**：論文正文與其補充材料（含 Data.xlsx 等檔案）受版權保護，**未隨附於\n",
    "此**。本筆記本僅內嵌一小段可歸因的事實性數據（9 個溶質在 4 種有機修飾劑比例下\n",
    "的滯留時間表）供教學重現使用，不構成對原著作權材料的重製。\n",
    "\n",
    "## 這個筆記本教你什麼？\n",
    "\n",
    "1. 什麼是滯留因子 k，怎麼從滯留時間 tR 算出來\n",
    "2. 怎麼用簡單的線性迴歸模型描述「有機修飾劑比例 f」與「滯留因子 k」的關係\n",
    "3. 怎麼用這個模型模擬層析圖（chromatogram）\n",
    "4. 怎麼定義並計算「解析度」（Resolution, Rs），以及如何掃描找出最佳分離條件\n",
    "\n",
    "## 需要的套件\n",
    "\n",
    "本筆記本只使用下列 Python 套件（皆為資料科學常見套件，不需要 plotly 或\n",
    "seaborn）：\n",
    "\n",
    "- `numpy`\n",
    "- `pandas`\n",
    "- `scipy`（只用 `scipy.stats.linregress` 做線性迴歸）\n",
    "- `matplotlib`\n",
    "\n",
    "建議初學者：本筆記本設計成**由上而下逐格執行**（Run All 或依序按 Shift+Enter）。\n",
    "每一步都會印出或畫出中間結果，讓你確認自己看得懂每一步在做什麼。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ca3f173f",
   "metadata": {},
   "source": [
    "## 資料\n",
    "\n",
    "我們使用論文補充材料 `Data.xlsx` 中 `i-ret.fit` 工作表內的一小段實驗數據：\n",
    "\n",
    "- 管柱：Kinetex 2.6 um XB-C18，150 x 4.6 mm\n",
    "- 移動相：乙腈（acetonitrile）／水，pH 5.7\n",
    "- 管柱死時間（void time）：t0 = 1.4 分鐘\n",
    "- 9 個溶質（A1 ~ A9），在 4 種有機修飾劑比例 f（0.40、0.45、0.50、0.60）下量測\n",
    "  滯留時間 tR（分鐘）\n",
    "\n",
    "下面直接把這張表格用 `pandas.DataFrame` 打進程式碼裡（寬格式：每一列是一個 f\n",
    "值，每一欄是一個溶質的 tR）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed5317ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "# 管柱死時間（分鐘）\n",
    "T0 = 1.4\n",
    "\n",
    "# 峰形參數（高斯峰模型；下一步會用到）：\n",
    "#   峰高     h = max(0.014, h0 + h1 * tR)\n",
    "#   峰寬參數 s = s0 + s1 * tR\n",
    "SHAPE = dict(h0=0.138, h1=-0.0038, s0=0.013, s1=0.0105)\n",
    "\n",
    "# 內嵌的滯留時間資料（寬格式）：欄位 f, A1..A9，數值為 tR（分鐘）\n",
    "retention_wide = pd.DataFrame(\n",
    "    {\n",
    "        \"f\":  [0.40,   0.45,   0.50,  0.60],\n",
    "        \"A1\": [3.735,  3.000,  2.537, 2.013],\n",
    "        \"A2\": [4.671,  3.650,  3.010, 2.265],\n",
    "        \"A3\": [6.646,  5.200,  4.208, 3.000],\n",
    "        \"A4\": [7.662,  5.600,  4.349, 2.940],\n",
    "        \"A5\": [8.498,  5.900,  4.430, 2.910],\n",
    "        \"A6\": [10.760, 6.850,  4.842, 3.061],\n",
    "        \"A7\": [11.090, 8.000,  6.051, 3.866],\n",
    "        \"A8\": [15.990, 10.380, 7.239, 4.198],\n",
    "        \"A9\": [19.000, 12.690, 8.921, 5.085],\n",
    "    }\n",
    ")\n",
    "\n",
    "retention_wide.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "545217e3",
   "metadata": {},
   "source": [
    "## 第一步：算 k 值\n",
    "\n",
    "滯留因子（retention factor）k 描述一個溶質「相對於流動相多花了多久時間」才\n",
    "從管柱中沖提出來，定義為：\n",
    "\n",
    "$$ k = \\frac{t_R - t_0}{t_0} $$\n",
    "\n",
    "其中：\n",
    "- \\(t_R\\) 是該溶質的滯留時間（實驗量測值）\n",
    "- \\(t_0\\) 是管柱死時間（在本例中 t0 = 1.4 分鐘）\n",
    "\n",
    "k 值不受管柱長度、流速等因素影響，是層析學裡描述「滯留強弱」最基本的量。\n",
    "\n",
    "下面的程式碼會把剛剛的寬格式資料轉成「長格式」（每一列是一個 solute + f 的\n",
    "組合），再算出每一列的 k 值。\n",
    "\n",
    "> **練習提示**：改成 t0 = 1.2 再跑一次，k 會怎麼變？（提示：t0 變小 -> 分母變小\n",
    "> 且 tR - t0 變大，k 應該會全面上升。你可以複製下面的 cell，把 `T0` 換成\n",
    "> `1.2` 試試看。）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5f12445d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 把寬格式 (f, A1..A9) 轉成長格式 (solute, f, tR)\n",
    "retention_long = retention_wide.melt(id_vars=\"f\", var_name=\"solute\", value_name=\"tR\")\n",
    "\n",
    "# 依照 A1..A9 的順序排列 solute（melt 預設會照欄位順序，這裡明確設定以確保穩定）\n",
    "solute_order = list(retention_wide.columns[1:])\n",
    "retention_long[\"solute\"] = pd.Categorical(retention_long[\"solute\"], categories=solute_order, ordered=True)\n",
    "retention_long = retention_long.sort_values([\"solute\", \"f\"]).reset_index(drop=True)\n",
    "\n",
    "# 計算 k = (tR - t0) / t0，主流程使用 t0 = 1.4\n",
    "retention_long[\"k\"] = (retention_long[\"tR\"] - T0) / T0\n",
    "\n",
    "retention_long.head(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e67dac97",
   "metadata": {},
   "source": [
    "## 第二步：擬合模型 1\n",
    "\n",
    "論文中的「模型 1」假設 ln k 與有機修飾劑比例 f 呈線性關係：\n",
    "\n",
    "$$ \\ln k = c_0 - c_1 \\cdot f $$\n",
    "\n",
    "對每一個溶質，我們用它在 4 個 f 值下量到的 (f, ln k) 資料點做一次簡單線性迴歸\n",
    "（OLS），求出截距與斜率：\n",
    "\n",
    "- 迴歸式寫成 `ln k = intercept + slope * f`\n",
    "- 對照模型 1：`c0 = intercept`，`c1 = -slope`（因為模型裡是「減」c1*f）\n",
    "\n",
    "我們用 `scipy.stats.linregress` 對 9 個溶質各自迴歸一次，並整理成一張參數表。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5a36706",
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import stats\n",
    "\n",
    "rows = []\n",
    "for solute, grp in retention_long.groupby(\"solute\", sort=False, observed=True):\n",
    "    f = grp[\"f\"].to_numpy(dtype=float)\n",
    "    lnk = np.log(grp[\"k\"].to_numpy(dtype=float))\n",
    "    n = len(f)\n",
    "\n",
    "    result = stats.linregress(f, lnk)\n",
    "    c0 = result.intercept\n",
    "    c1 = -result.slope\n",
    "\n",
    "    fitted = result.intercept + result.slope * f\n",
    "    ssr = np.sum((lnk - fitted) ** 2)\n",
    "    sigma = np.sqrt(ssr / (n - 2))\n",
    "\n",
    "    rows.append(dict(solute=solute, c0=c0, c1=c1, r_squared=result.rvalue ** 2, sigma=sigma))\n",
    "\n",
    "params = pd.DataFrame(rows)\n",
    "\n",
    "display_params = params.copy()\n",
    "for col in [\"c0\", \"c1\", \"r_squared\", \"sigma\"]:\n",
    "    display_params[col] = display_params[col].round(4)\n",
    "print(display_params.to_string(index=False))\n",
    "\n",
    "print(f\"\\nc0 範圍：{params['c0'].min():.4f} ~ {params['c0'].max():.4f}\")\n",
    "print(f\"c1 範圍：{params['c1'].min():.4f} ~ {params['c1'].max():.4f}\")\n",
    "print(\"論文 i-optim 工作表公布值：c0 3.1404-5.6030，c1 5.9113-8.5366\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cca518c6",
   "metadata": {},
   "source": [
    "## 第三步：畫 ln k 對 f\n",
    "\n",
    "有了每個溶質的 c0、c1，我們可以把「實測的 (f, ln k) 點」和「擬合出來的直線」\n",
    "畫在同一張圖上，直觀檢查模型 1 是否適合描述這些資料（點應該幾乎都落在直線上）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffc42ca8",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# 圖表中文字型：優先使用 Windows 內建的正黑體／微軟雅黑，若系統沒有這些字型則\n",
    "# 靜默退回預設字型（僅影響圖片上文字是否顯示為方框，不影響任何計算結果）。\n",
    "plt.rcParams[\"font.sans-serif\"] = [\"Microsoft JhengHei\", \"Microsoft YaHei\", \"SimHei\", \"DejaVu Sans\"]\n",
    "plt.rcParams[\"axes.unicode_minus\"] = False\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 4.5))\n",
    "\n",
    "solutes = list(params[\"solute\"])\n",
    "colors = plt.cm.tab10(np.linspace(0, 1, len(solutes)))\n",
    "f_line = np.linspace(retention_long[\"f\"].min() - 0.02, retention_long[\"f\"].max() + 0.02, 100)\n",
    "\n",
    "for solute, color in zip(solutes, colors):\n",
    "    grp = retention_long[retention_long[\"solute\"] == solute]\n",
    "    lnk = np.log(grp[\"k\"].to_numpy(dtype=float))\n",
    "    ax.scatter(grp[\"f\"], lnk, color=color, s=25, zorder=3)\n",
    "\n",
    "    row = params.loc[params[\"solute\"] == solute].iloc[0]\n",
    "    ax.plot(f_line, row[\"c0\"] - row[\"c1\"] * f_line, color=color, linewidth=1.3, label=str(solute))\n",
    "\n",
    "ax.set_xlabel(\"有機修飾劑比例 f\")\n",
    "ax.set_ylabel(\"ln k\")\n",
    "ax.set_title(\"模型 1：ln k = c0 - c1 * f\")\n",
    "ax.legend(fontsize=8, ncol=3, loc=\"upper right\")\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86311e29",
   "metadata": {},
   "source": [
    "## 第四步：模擬層析圖\n",
    "\n",
    "有了模型 1 的參數，在給定 f 之下，我們可以：\n",
    "\n",
    "1. 用 `k = exp(c0 - c1 * f)` 算出每個溶質的 k 值\n",
    "2. 用 `tR = t0 * (1 + k)` 算出滯留時間\n",
    "3. 用峰形參數算出每個峰的高度與寬度：\n",
    "   - 峰高 `h = max(0.014, h0 + h1 * tR)`（設下限避免晚洗出的峰高變負值）\n",
    "   - 峰寬參數 `s = s0 + s1 * tR`，基底峰寬 `w = 4*s/sqrt(2)`\n",
    "4. 把每個溶質的訊號視為一個高斯峰，加總起來就是模擬層析圖：\n",
    "\n",
    "$$ \\text{signal}(t) = \\sum_j h_j \\cdot \\exp\\left(-\\left(\\frac{t - t_{R,j}}{s_j}\\right)^2\\right) $$\n",
    "\n",
    "下面這個 cell 是**互動式**的：把 `f_choice` 改成別的數值（例如 0.40、0.50、\n",
    "0.60），再重新執行這個 cell，觀察層析圖如何變化（峰的位置、間距、是否重疊）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d25dbc50",
   "metadata": {},
   "outputs": [],
   "source": [
    "f_choice = 0.45  # <-- 學生可以改這個值，改完重新執行本 cell 看看層析圖怎麼變\n",
    "\n",
    "# 由模型 1 參數預測在 f_choice 之下的滯留時間與峰形\n",
    "peaks = params[[\"solute\", \"c0\", \"c1\"]].copy()\n",
    "peaks[\"f\"] = f_choice\n",
    "peaks[\"k\"] = np.exp(peaks[\"c0\"] - peaks[\"c1\"] * f_choice)\n",
    "peaks[\"tR\"] = T0 * (1 + peaks[\"k\"])\n",
    "peaks[\"s\"] = SHAPE[\"s0\"] + SHAPE[\"s1\"] * peaks[\"tR\"]\n",
    "peaks[\"w\"] = 4 * peaks[\"s\"] / np.sqrt(2)\n",
    "peaks[\"h\"] = np.maximum(0.014, SHAPE[\"h0\"] + SHAPE[\"h1\"] * peaks[\"tR\"])\n",
    "peaks = peaks.drop(columns=[\"c0\", \"c1\"]).sort_values(\"tR\").reset_index(drop=True)\n",
    "\n",
    "print(peaks.round(4))\n",
    "\n",
    "# 模擬層析圖：時間軸從 0 到最大 tR 的 1.15 倍\n",
    "dt = 0.002\n",
    "t_max_sim = float(peaks[\"tR\"].max()) * 1.15\n",
    "t = np.arange(0, t_max_sim + dt, dt)\n",
    "\n",
    "signal = np.zeros_like(t)\n",
    "for _, row in peaks.iterrows():\n",
    "    signal += row[\"h\"] * np.exp(-(((t - row[\"tR\"]) / row[\"s\"]) ** 2))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6.5, 4.5))\n",
    "ax.plot(t, signal, color=\"#0d7377\", linewidth=1.0)\n",
    "for _, row in peaks.iterrows():\n",
    "    ax.annotate(str(row[\"solute\"]), (row[\"tR\"], row[\"h\"]), textcoords=\"offset points\",\n",
    "                xytext=(0, 4), ha=\"center\", fontsize=8, color=\"grey\")\n",
    "ax.set_xlabel(\"時間 t / 分鐘\")\n",
    "ax.set_ylabel(\"訊號\")\n",
    "ax.set_title(f\"模擬層析圖（f = {f_choice}）\")\n",
    "ax.set_ylim(top=peaks[\"h\"].max() * 1.25)\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b033b81f",
   "metadata": {},
   "source": [
    "## 第五步：解析度與最佳化\n",
    "\n",
    "「解析度」（Resolution, Rs）描述兩個相鄰的峰分離得好不好，定義為：\n",
    "\n",
    "$$ R_s = \\frac{2 (t_{R,2} - t_{R,1})}{w_1 + w_2} $$\n",
    "\n",
    "Rs 越大代表兩個峰分得越開；一般認為 Rs >= 1.5 算是「基線分離」。對一個混合物\n",
    "而言，我們關心的是**所有相鄰峰對中最差（最小）的 Rs**，因為那就是整個分離方法\n",
    "的瓶頸。\n",
    "\n",
    "**最佳化的想法**：掃描一系列的 f 值（例如 0.30 到 0.60，每次增加 0.005），對每\n",
    "個 f：\n",
    "\n",
    "1. 用模型 1 算出所有溶質的 tR 與峰形\n",
    "2. 算出所有相鄰峰對的 Rs，取其中最小值\n",
    "3. 檢查最晚出峰的時間 `tR_max` 是否在時間預算 `t_max` 之內（可行）\n",
    "\n",
    "在所有「可行」（`tR_max <= t_max`）的 f 之中，選出「最小 Rs 最大」的那個 f，\n",
    "就是在時間限制之內最好的分離條件。\n",
    "\n",
    "下面的程式碼實作這個掃描（用簡單的迴圈，重點在於清楚易懂，不追求效能），並畫\n",
    "出「最小 Rs 對 f」與「tR_max 對 f」的關係圖。你會看到在 f ≈ 0.35 附近，出現一\n",
    "個解析度驟降的「共沖提」（co-elution）低點 —— 這是因為在那附近兩個峰的出峰\n",
    "順序正好交換，中間有一瞬間幾乎完全重疊。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7fd0d90",
   "metadata": {},
   "outputs": [],
   "source": [
    "def predict_isocratic(params, f, t0=T0):\n",
    "    '''由模型 1 參數，預測恆溶劑（isocratic）條件下每個溶質的滯留時間與峰形。'''\n",
    "    out = params[[\"solute\", \"c0\", \"c1\"]].copy()\n",
    "    out[\"f\"] = f\n",
    "    out[\"k\"] = np.exp(out[\"c0\"] - out[\"c1\"] * f)\n",
    "    out[\"tR\"] = t0 * (1 + out[\"k\"])\n",
    "    out[\"s\"] = SHAPE[\"s0\"] + SHAPE[\"s1\"] * out[\"tR\"]\n",
    "    out[\"w\"] = 4 * out[\"s\"] / np.sqrt(2)\n",
    "    out[\"h\"] = np.maximum(0.014, SHAPE[\"h0\"] + SHAPE[\"h1\"] * out[\"tR\"])\n",
    "    out = out.drop(columns=[\"c0\", \"c1\"]).sort_values(\"tR\").reset_index(drop=True)\n",
    "    return out\n",
    "\n",
    "\n",
    "def min_resolution(peaks_df):\n",
    "    '''相鄰峰對中解析度最差（最小 Rs）的一對。'''\n",
    "    p = peaks_df.sort_values(\"tR\").reset_index(drop=True)\n",
    "    solute = p[\"solute\"].astype(str).to_numpy()\n",
    "    tR = p[\"tR\"].to_numpy(dtype=float)\n",
    "    w = p[\"w\"].to_numpy(dtype=float)\n",
    "    if len(p) < 2:\n",
    "        return dict(Rs=np.nan, pair=None)\n",
    "    Rs_values = 2 * (tR[1:] - tR[:-1]) / (w[1:] + w[:-1])\n",
    "    i = int(np.argmin(Rs_values))\n",
    "    pair = f\"{solute[i]}/{solute[i + 1]}\"\n",
    "    return dict(Rs=float(Rs_values[i]), pair=pair)\n",
    "\n",
    "\n",
    "def optimise_isocratic(params, f_range=(0.30, 0.60), df=0.005, t_max=20, t0=T0):\n",
    "    '''掃描 f，找出時間預算內解析度最好的條件。'''\n",
    "    f_grid = np.arange(f_range[0], f_range[1] + df / 2, df)\n",
    "    rows = []\n",
    "    for f in f_grid:\n",
    "        pk = predict_isocratic(params, float(f), t0)\n",
    "        worst = min_resolution(pk)\n",
    "        rows.append(dict(f=float(f), Rs=worst[\"Rs\"], pair=worst[\"pair\"], tR_max=float(pk[\"tR\"].max())))\n",
    "    scan = pd.DataFrame(rows)\n",
    "    scan[\"feasible\"] = scan[\"tR_max\"] <= t_max\n",
    "\n",
    "    feasible = scan[scan[\"feasible\"]]\n",
    "    best = feasible.loc[feasible[\"Rs\"].idxmax()] if not feasible.empty else None\n",
    "    return scan, best\n",
    "\n",
    "\n",
    "t_max = 20\n",
    "scan, best = optimise_isocratic(params, t_max=t_max)\n",
    "\n",
    "print(f\"== 最佳化結果 (t_max = {t_max}) ==\")\n",
    "print(f\" 最佳 f = {best['f']:.3f}   最小 Rs = {best['Rs']:.4f}   \"\n",
    "      f\"tR_max = {best['tR_max']:.2f} 分鐘   最難分開的峰對 = {best['pair']}\")\n",
    "\n",
    "# 畫出「最小 Rs 對 f」與「tR_max 對 f」（雙 y 軸）\n",
    "fig, ax1 = plt.subplots(figsize=(6.5, 4.5))\n",
    "ax1.axhline(1.5, linestyle=\"--\", color=\"#2e7d4f\", linewidth=1, label=\"Rs = 1.5 參考線\")\n",
    "ax1.plot(scan[\"f\"], scan[\"Rs\"], color=\"#0d7377\", linewidth=1.5, label=\"最小 Rs\")\n",
    "ax1.scatter([best[\"f\"]], [best[\"Rs\"]], color=\"#2e7d4f\", s=60, zorder=5, label=\"最佳 f\")\n",
    "ax1.set_xlabel(\"有機修飾劑比例 f\")\n",
    "ax1.set_ylabel(\"最小 Rs\", color=\"#0d7377\")\n",
    "\n",
    "ax2 = ax1.twinx()\n",
    "ax2.plot(scan[\"f\"], scan[\"tR_max\"], color=\"#c8632b\", linewidth=1.5, label=\"tR_max\")\n",
    "ax2.axhline(t_max, linestyle=\":\", color=\"#c8632b\", linewidth=1)\n",
    "ax2.set_ylabel(\"tR_max / 分鐘\", color=\"#c8632b\")\n",
    "\n",
    "ax1.set_title(f\"掃描最佳化：最佳 f = {best['f']:.3f}, Rs = {best['Rs']:.3f}\")\n",
    "fig.tight_layout()\n",
    "plt.show()\n",
    "\n",
    "# 共沖提檢查：f = 0.350 附近的最小 Rs 值\n",
    "worst_350 = min_resolution(predict_isocratic(params, 0.350))\n",
    "print(f\"\\n== f = 0.350 的共沖提檢查 ==\")\n",
    "print(f\" 最小 Rs = {worst_350['Rs']:.4f}（{worst_350['pair']}）—— 這就是圖上 f≈0.35 附近的驟降點\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6087619b",
   "metadata": {},
   "source": [
    "## 練習題\n",
    "\n",
    "以下三題請你動手修改、重新執行程式碼來回答。可以複製上面的 cell 內容，或直接\n",
    "呼叫上面已經定義好的函式（`optimise_isocratic`、`predict_isocratic`、\n",
    "`min_resolution` 等）。"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5330c0e9",
   "metadata": {},
   "source": [
    "### (a) 把 t_max 從 20 改成 12，看看新的最佳 f 和 Rs 是多少\n",
    "\n",
    "**提示**：重新呼叫 `optimise_isocratic(params, t_max=12)`，比較新的 `best['f']`\n",
    "和 `best['Rs']` 跟 t_max=20 時的結果有什麼不同。時間預算變嚴格時，最佳 f\n",
    "通常會需要往哪個方向調整？"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4a89932",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 你的程式碼："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e32bc079",
   "metadata": {},
   "source": [
    "### (b) 從混合物中移除 A6，看看最佳條件會不會改變\n",
    "\n",
    "**提示**：先用 `params_no_A6 = params[params['solute'] != 'A6'].reset_index(drop=True)`\n",
    "拿掉 A6 那一列，再用 `optimise_isocratic(params_no_A6, t_max=20)` 重新最佳化，\n",
    "比較最佳 f、Rs、以及最難分開的峰對跟原本的結果有何不同。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f1715e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 你的程式碼："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25a9aac5",
   "metadata": {},
   "source": [
    "### (c) 改變峰寬參數 s1，觀察 Rs 如何變化\n",
    "\n",
    "**提示**：`SHAPE` 字典裡的 `s1` 決定了峰寬隨滯留時間增加的速度\n",
    "（`s = s0 + s1 * tR`）。你可以複製一份 `SHAPE`（例如\n",
    "`shape2 = dict(SHAPE); shape2['s1'] = 0.02`），修改 `predict_isocratic` 或直接\n",
    "在 cell 裡重新計算 `s`、`w`，再看看同一個 f 之下的最小 Rs 會如何隨 s1 變化\n",
    "（s1 越大，峰越寬，Rs 應該會怎麼變？）。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "924cca51",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 你的程式碼："
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e605149c",
   "metadata": {},
   "source": [
    "## 延伸閱讀\n",
    "\n",
    "本筆記本示範的是簡化過、方便教學的版本。如果你想看完整、包含梯度沖提\n",
    "（gradient elution）數值解、更多繪圖與正確性檢查的「生產版」實作，可以參考同一\n",
    "個資料夾中的：\n",
    "\n",
    "- `rchromoptim_modern.py` —— 完整的 Python 版本（向量化、模組化的函式庫）\n",
    "- `rchromoptim_modern.R` —— 功能相同的 R 語言版本\n",
    "\n",
    "這兩個檔案實作了與本筆記本相同的核心邏輯（模型 1 擬合、恆溶劑與梯度沖提滯留時\n",
    "間預測、解析度計算、最佳化掃描），並額外提供梯度沖提模擬、更完整的繪圖函式，\n",
    "適合想深入研究或應用在自己資料上的同學參考。"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}