{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.0"
  }
 },
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "# Gas Permeability Prediction from ATR-FTIR Spectroscopy\n### CatBoost / GBR / SVR with Leverage-Based Applicability Domain Filtering\n\nThis notebook trains a CatBoost model to predict gas permeability (log10-transformed, Barrer)  \nfrom ATR-FTIR spectral features combined with one-hot encoded polymer-family metadata.\n\n**Pipeline overview**\n1. Load permeability data and polymer metadata\n2. Load and interpolate ATR-FTIR spectra onto a common wavenumber grid\n3. Compress spectra to PCA scores (retaining 95% variance)\n4. Train CatBoost with 5-fold cross-validation (fixed hyperparameters)\n5. Apply Williams plot (leverage method) to identify out-of-domain records\n6. Retrain on in-domain data and evaluate on both scales (log and original)\n7. Save trained model bundle for inference\n\n**Required files**\n- `Dataset.csv` - polymer permeability values and metadata\n- `graphs.csv` - spectrum-to-polymer mapping (new_key, Item)\n- `FTIR_Spectra/` - folder of individual ATR-FTIR CSV files (no header; columns: wavenumber, absorbance)\n"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 1. Dependencies"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# Install CatBoost if not already present (Colab)\n!pip install catboost -q\n\nimport os\nimport re\nimport json\nimport joblib\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom datetime import datetime\nfrom scipy.interpolate import interp1d\nfrom scipy.signal import savgol_filter\nfrom sklearn.decomposition import PCA\nfrom sklearn.model_selection import KFold, GridSearchCV\nfrom sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.svm import SVR\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.ensemble import GradientBoostingRegressor\nfrom catboost import CatBoostRegressor"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 2. Configuration\n\nSet all run parameters here. No changes are needed elsewhere in the notebook.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# \u2500\u2500 Target gas \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n# Options: 'CO2', 'He', 'CH4', 'H2', 'N2', 'O2'\nGAS_TARGET = 'CO2'\n\n# \u2500\u2500 Model family \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n# 'catboost' : CatBoostRegressor (recommended; best performance on this dataset)\n# 'gbr'      : sklearn GradientBoostingRegressor\n# 'svr'      : Support Vector Regression (with StandardScaler pre-processing)\nMODEL_FAMILY = 'catboost'\n\n# \u2500\u2500 Grid search toggle \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n# True  : run GridSearchCV at both training steps (slower, finds best params)\n# False : use fixed hyperparameters from FIXED_PARAMS below (faster)\nUSE_GRID_SEARCH = False\n\n# \u2500\u2500 Fixed hyperparameters (used when USE_GRID_SEARCH = False) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n# These were selected from grid search on the CO2 training set.\n# If switching gas or model family, consider running a grid search first.\nFIXED_PARAMS = {\n    'catboost': {\n        'depth'        : 6,\n        'iterations'   : 400,\n        'l2_leaf_reg'  : 5,\n        'learning_rate': 0.05,\n    },\n    'gbr': {\n        'n_estimators'    : 200,\n        'learning_rate'   : 0.05,\n        'max_depth'       : 4,\n        'subsample'       : 0.9,\n        'min_samples_split': 8,\n        'min_samples_leaf' : 3,\n        'max_features'    : 'sqrt',\n        'loss'            : 'squared_error',\n    },\n    'svr': {\n        # SVR params are prefixed with 'svr__' because the model is wrapped\n        # in a Pipeline with a StandardScaler step.\n        'svr__kernel' : 'rbf',\n        'svr__C'      : 10.0,\n        'svr__gamma'  : 'scale',\n        'svr__epsilon': 0.05,\n    },\n}\n\n# \u2500\u2500 FTIR feature representation \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n# 'pca'    : PCA scores retaining 95% cumulative variance (recommended)\n# 'smooth' : Savitzky-Golay smoothed spectrum + first derivative (concatenated)\n# 'raw'    : Interpolated absorbance values with no dimensionality reduction\nFTIR_MODE = 'pca'\n\n# \u2500\u2500 Data paths \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nMETA_PATH   = 'Dataset.csv'       # polymer metadata and permeability\nGRAPHS_PATH = 'graphs.csv'        # spectrum-to-polymer mapping (new_key, Item)\nFTIR_DIR    = '/content/Data'     # folder of individual FTIR CSV files\n\n# \u2500\u2500 Cross-validation \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nN_FOLDS         = 5\nCV_RANDOM_STATE = 40\n\n# \u2500\u2500 PCA variance threshold \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nPCA_VARIANCE_THRESHOLD = 0.95\n\n# \u2500\u2500 Williams plot thresholds \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nSTD_RESID_THRESHOLD = 2.0   # |standardised residual| > this flags a record\n# Leverage threshold h* = 2p/n is computed automatically"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 3. Model Definition and Parameter Grids\n\nParameter grids are used when `USE_GRID_SEARCH = True`. Each model has a standard grid  \nsuitable for this dataset size (~800 spectra, ~70-90 PCA features).\n\n**CatBoost** is the recommended choice. SVR requires StandardScaler (handled automatically  \nvia a Pipeline). GBR is sklearn's gradient boosting and requires no special wrapping.\n\nWhen `USE_GRID_SEARCH = False`, the fixed values in `FIXED_PARAMS` are used instead  \nand no grid search is performed.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# \u2500\u2500 Parameter grids (used when USE_GRID_SEARCH = True) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\n\nparam_grids = {\n    'catboost': {\n        'depth'        : [4, 6],\n        'learning_rate': [0.05, 0.1],\n        'iterations'   : [400, 600],\n        'l2_leaf_reg'  : [3, 5, 10],\n    },\n    'gbr': {\n        'n_estimators'    : [200, 500],\n        'learning_rate'   : [0.05, 0.1],\n        'max_depth'       : [3, 4, 6],\n        'subsample'       : [0.9],\n        'min_samples_split': [8, 10],\n        'min_samples_leaf' : [3, 4],\n        'max_features'    : ['sqrt'],\n        'loss'            : ['squared_error'],\n    },\n    'svr': {\n        # Prefixed with 'svr__' to match the Pipeline step name\n        'svr__kernel' : ['rbf'],\n        'svr__C'      : [0.1, 1.0, 10.0, 50.0],\n        'svr__gamma'  : ['scale', 0.01, 0.05],\n        'svr__epsilon': [0.01, 0.05, 0.1],\n    },\n}\n\n\ndef build_model(model_family):\n    \"\"\"\n    Return a freshly instantiated model for the given family.\n\n    SVR is wrapped in a Pipeline with StandardScaler so that feature scaling\n    is applied consistently within each CV fold.\n\n    Parameters\n    ----------\n    model_family : str  One of 'catboost', 'gbr', 'svr'\n\n    Returns\n    -------\n    model : estimator instance (unfitted)\n    \"\"\"\n    if model_family == 'catboost':\n        return CatBoostRegressor(\n            loss_function='RMSE',\n            verbose=False,\n            random_seed=42,\n        )\n    elif model_family == 'gbr':\n        return GradientBoostingRegressor(random_state=42)\n    elif model_family == 'svr':\n        return Pipeline([\n            ('scaler', StandardScaler()),\n            ('svr', SVR()),\n        ])\n    else:\n        raise ValueError(f\"Unknown MODEL_FAMILY '{model_family}'. \"\n                         f\"Choose 'catboost', 'gbr', or 'svr'.\")\n\n\ndef get_params(model_family, use_grid_search):\n    \"\"\"\n    Return either the parameter grid (for grid search) or a single-value\n    grid wrapping the fixed params (for fixed-param CV).\n\n    Returns a dict suitable for passing to GridSearchCV in both cases,\n    which keeps the training logic uniform regardless of USE_GRID_SEARCH.\n    \"\"\"\n    if use_grid_search:\n        return param_grids[model_family]\n    else:\n        # Wrap each fixed value in a list so GridSearchCV accepts it\n        # but only evaluates one combination (no real search)\n        return {k: [v] for k, v in FIXED_PARAMS[model_family].items()}\n\n\n# Validate config\nbuild_model(MODEL_FAMILY)   # will raise immediately if MODEL_FAMILY is invalid\nprint(f\"Model family : {MODEL_FAMILY}\")\nprint(f\"Grid search  : {USE_GRID_SEARCH}\")\nif USE_GRID_SEARCH:\n    print(f\"Parameter grid ({MODEL_FAMILY}):\")\n    for k, v in param_grids[MODEL_FAMILY].items():\n        print(f\"  {k}: {v}\")\nelse:\n    print(f\"Fixed parameters ({MODEL_FAMILY}):\")\n    for k, v in FIXED_PARAMS[MODEL_FAMILY].items():\n        print(f\"  {k}: {v}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 4. Load Polymer Metadata and Permeability Data\n\n`Dataset.csv` contains permeability values (Barrer) for six gases and one-hot encoded  \npolymer-family metadata. `graphs.csv` maps each FTIR spectrum ID (`new_key`) to its  \npolymer ID (`Item`).\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# Load polymer metadata (permeability + one-hot encoded features)\nmeta = pd.read_csv(META_PATH)\nmeta['Item'] = meta['Item'].astype(int)\n\n# Load spectrum-to-polymer mapping\ngraphs_df = pd.read_csv(GRAPHS_PATH, header=None, names=['new_key', 'Item'])\ngraphs_df['Item']    = graphs_df['Item'].astype(int)\ngraphs_df['new_key'] = graphs_df['new_key'].astype(int)\n\n# Inner join: keep only spectra that have a matching polymer record\nmerged = graphs_df.merge(meta, on='Item', how='left')\n\n# One-hot encoded polymer-family columns used as supplementary model features\nohe_columns = [\n    'Category_Amorphous', 'Category_Semi-crystalline',\n    'Chemical_Family_Acrylic', 'Chemical_Family_Fluoropolymer',\n    'Chemical_Family_Phenolic', 'Chemical_Family_Polyacetal',\n    'Chemical_Family_Polyamide', 'Chemical_Family_Polyanhydride',\n    'Chemical_Family_Polycarbonate', 'Chemical_Family_Polydiene',\n    'Chemical_Family_Polyester', 'Chemical_Family_Polyether',\n    'Chemical_Family_Polyketone', 'Chemical_Family_Polyolefin',\n    'Chemical_Family_Polystyrene', 'Chemical_Family_Polysulfone',\n    'Chemical_Family_Polyurethane', 'Chemical_Family_Polyvinyl',\n]\n\nmissing = [c for c in ohe_columns if c not in merged.columns]\nif missing:\n    print(f\"Warning: OHE columns not found in metadata: {missing}\")\n\n# Polymer-family feature matrix keyed by new_key\npolymer_features_df = merged.set_index('new_key')[ohe_columns]\n\n# Drop rows with no permeability value for the selected gas\ngas_col = 'He (Barrer)' if GAS_TARGET == 'He' else GAS_TARGET\nmerged = merged.dropna(subset=[gas_col])\n\n# Target series keyed by new_key\ny_series = merged.set_index('new_key')[gas_col]\n\nprint(f\"Polymers with {GAS_TARGET} permeability data: {merged['Item'].nunique()}\")\nprint(f\"Spectra with {GAS_TARGET} permeability data : {len(y_series)}\")\nprint(f\"\\n{GAS_TARGET} permeability summary (Barrer):\")\nprint(y_series.describe().round(3))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 5. Load and Interpolate ATR-FTIR Spectra\n\nEach spectrum is stored as an individual CSV file (no header; columns: wavenumber in cm-1,  \nabsorbance in arbitrary units). All spectra are interpolated onto a common wavenumber grid  \nspanning the union of all recorded ranges.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "def extract_number(filename):\n    \"\"\"Extract the leading integer from a filename (used as spectrum ID).\"\"\"\n    match = re.search(r'\\d+', filename)\n    return int(match.group()) if match else float('inf')\n\n\n# Load all FTIR CSV files from the data directory\nfile_names = sorted(os.listdir(FTIR_DIR), key=extract_number)\ncsv_data = {}\n\nfor fname in file_names:\n    if fname.endswith('.csv'):\n        key = fname.rsplit('.', 1)[0]\n        fpath = os.path.join(FTIR_DIR, fname)\n        df = pd.read_csv(fpath, names=['wavelength', 'absorbance'], header=None)\n        csv_data[key] = df\n\nprint(f\"FTIR spectra loaded: {len(csv_data)}\")\n\n# Determine common wavenumber grid from the union of all spectral ranges\nmin_wl = min(df.iloc[:, 0].min() for df in csv_data.values())\nmax_wl = max(df.iloc[:, 0].max() for df in csv_data.values())\nmax_pts = max(len(df) for df in csv_data.values())\ncommon_wavelengths = np.linspace(min_wl, max_wl, max_pts)\n\nprint(f\"Common wavenumber grid: {min_wl:.1f} to {max_wl:.1f} cm-1 ({max_pts} points)\")\n\n# Build nearest-neighbour interpolation functions for each spectrum\ninterp_func_dict = {\n    key: interp1d(df.iloc[:, 0], df.iloc[:, 1], kind='nearest', fill_value='extrapolate')\n    for key, df in csv_data.items()\n    if len(df) <= max_pts\n}\n\n# Assemble interpolated spectra into a matrix (rows = spectra, index = new_key)\nintensity_list, key_list = [], []\nfor key in interp_func_dict:\n    intensity_list.append(interp_func_dict[key](common_wavelengths))\n    key_list.append(int(extract_number(key)))\n\ncombined_df = pd.DataFrame(intensity_list, index=key_list)\ncombined_df.index.name = 'new_key'\ncombined_df = combined_df.dropna()\n\nprint(f\"Interpolated FTIR matrix shape: {combined_df.shape}  (spectra x wavenumber points)\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 6. Feature Engineering\n\nFTIR spectra are compressed using one of three representations controlled by `FTIR_MODE`:\n\n- **pca**: Principal Component Analysis retaining 95% of cumulative variance (default)\n- **smooth**: Savitzky-Golay smoothed spectrum concatenated with its first derivative\n- **raw**: Interpolated absorbance values without dimensionality reduction\n\nThe resulting spectral features are then concatenated with the one-hot encoded  \npolymer-family metadata to form the final feature matrix.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "if FTIR_MODE == 'pca':\n    # Fit PCA on all spectra; retain components up to the variance threshold\n    pca_full = PCA().fit(combined_df)\n    cumulative_variance = np.cumsum(pca_full.explained_variance_ratio_)\n    n_components = int(np.argmax(cumulative_variance >= PCA_VARIANCE_THRESHOLD) + 1)\n\n    pca_final = PCA(n_components=n_components).fit(combined_df)\n    X_spec = pd.DataFrame(\n        pca_final.transform(combined_df),\n        index=combined_df.index\n    )\n    print(f\"[PCA] Components retained: {n_components} ({PCA_VARIANCE_THRESHOLD:.0%} variance)\")\n    print(f\"X_spec shape: {X_spec.shape}\")\n\nelif FTIR_MODE == 'smooth':\n    # Downsample by factor 2 for speed, then apply Savitzky-Golay filter\n    STEP = 2\n    X_ds = combined_df.iloc[:, ::STEP].copy()\n    X_mat = X_ds.to_numpy(dtype=np.float32)\n\n    X_smooth_mat = savgol_filter(X_mat, window_length=11, polyorder=3, axis=1)\n    X_deriv_mat  = savgol_filter(X_mat, window_length=11, polyorder=3, deriv=1, axis=1)\n\n    wl = X_ds.columns.astype(str)\n    X_smooth = pd.DataFrame(X_smooth_mat, index=X_ds.index, columns=wl)\n    X_deriv  = pd.DataFrame(X_deriv_mat,  index=X_ds.index, columns=[c + '_d1' for c in wl])\n    X_spec   = pd.concat([X_smooth, X_deriv], axis=1)\n    pca_final = None\n    print(f\"[Smooth] X_spec shape: {X_spec.shape}\")\n\nelif FTIR_MODE == 'raw':\n    X_spec = combined_df.copy()\n    pca_final = None\n    print(f\"[Raw] X_spec shape: {X_spec.shape}\")\n\nelse:\n    raise ValueError(f\"Unknown FTIR_MODE: '{FTIR_MODE}'. Choose 'pca', 'smooth', or 'raw'.\")\n\nX_spec.index = combined_df.index\nX_spec.columns = X_spec.columns.astype(str)\n\n# Align spectral features with polymer-family metadata and permeability targets\ncommon_keys = X_spec.index.intersection(y_series.index)\nX_ftir   = X_spec.loc[common_keys]\nX_meta   = polymer_features_df.loc[common_keys]\nX_all    = pd.concat([X_ftir, X_meta], axis=1)\ny_all    = y_series.loc[common_keys]\ny_all_log = np.log10(y_all)\n\nprint(f\"\\nFinal feature matrix shape : {X_all.shape}\")\nprint(f\"Target variable            : log10({GAS_TARGET}) permeability\")\nprint(f\"Samples                    : {len(y_all_log)}\")\nprint(f\"\\nlog10({GAS_TARGET}) summary:\")\nprint(y_all_log.describe().round(3))"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 7. Initial Cross-Validation (Full Dataset)\n\n5-fold cross-validation on the full dataset before leverage filtering. When\n`USE_GRID_SEARCH = True`, a grid search runs inside each outer fold (3-fold inner CV)\nto find the best hyperparameters. When `USE_GRID_SEARCH = False`, the fixed values\nfrom `FIXED_PARAMS` are used directly.\n\nThe mean RMSE from this step is used as the sigma estimate in the Williams plot\nstandardised residual calculation.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "def run_cv(X, y, model_family, use_grid_search, n_folds=5, random_state=40):\n    \"\"\"\n    Run KFold cross-validation with GridSearchCV inside each fold.\n\n    When use_grid_search=False the parameter grid contains only the fixed values,\n    so GridSearchCV evaluates a single combination -- equivalent to plain CV\n    but keeping the code path identical in both modes.\n\n    Parameters\n    ----------\n    X              : pd.DataFrame  Feature matrix\n    y              : pd.Series     Log10-transformed target\n    model_family   : str           'catboost', 'gbr', or 'svr'\n    use_grid_search: bool          Whether to search over multiple param combos\n    n_folds        : int           Number of CV folds\n    random_state   : int           KFold random seed\n\n    Returns\n    -------\n    metrics_df  : pd.DataFrame  Per-fold R2, RMSE, MAE and best params\n    y_pred_oof  : pd.Series     Out-of-fold predictions (same index as y)\n    best_params : dict          Best parameters from the final fold (indicative)\n    \"\"\"\n    kf = KFold(n_splits=n_folds, shuffle=True, random_state=random_state)\n    param_grid = get_params(model_family, use_grid_search)\n    fold_metrics = []\n    y_pred_oof = pd.Series(np.nan, index=y.index, dtype=float)\n\n    for fold, (train_idx, test_idx) in enumerate(kf.split(X), 1):\n        X_tr, X_te = X.iloc[train_idx], X.iloc[test_idx]\n        y_tr, y_te = y.iloc[train_idx], y.iloc[test_idx]\n\n        gs = GridSearchCV(\n            build_model(model_family),\n            param_grid,\n            cv=3,               # inner CV for grid search\n            scoring='r2',\n            refit=True,\n            n_jobs=1,\n        )\n        gs.fit(X_tr, y_tr.squeeze())\n        y_pred = gs.predict(X_te)\n        y_pred_oof.iloc[test_idx] = y_pred\n\n        metrics = {\n            'Fold'       : fold,\n            'R2'         : r2_score(y_te, y_pred),\n            'RMSE'       : np.sqrt(mean_squared_error(y_te, y_pred)),\n            'MAE'        : mean_absolute_error(y_te, y_pred),\n            'best_params': gs.best_params_,\n        }\n        fold_metrics.append(metrics)\n        print(f\"  Fold {fold}: R2={metrics['R2']:.3f}  \"\n              f\"RMSE={metrics['RMSE']:.3f}  MAE={metrics['MAE']:.3f}  \"\n              f\"params={gs.best_params_}\")\n\n    metrics_df = pd.DataFrame(fold_metrics).set_index('Fold')\n    print(\"\\n  Mean across folds:\")\n    print(metrics_df[['R2', 'RMSE', 'MAE']].mean().round(3).to_string())\n\n    best_params = fold_metrics[-1]['best_params']\n    return metrics_df, y_pred_oof, best_params\n\n\n# \u2500\u2500 Run initial CV on full dataset (pre-leverage filtering) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nprint(f\"Initial CV: {GAS_TARGET}, model={MODEL_FAMILY}, \"\n      f\"grid_search={USE_GRID_SEARCH}, n={len(y_all_log)}\")\nprint(\"-\" * 60)\ncv_metrics_initial, y_pred_oof_initial, best_params_initial = run_cv(\n    X_all, y_all_log, MODEL_FAMILY, USE_GRID_SEARCH, N_FOLDS, CV_RANDOM_STATE\n)\n\n# Train a full-dataset model to generate predictions for the Williams plot.\n# Uses best params from CV if grid search was on, otherwise fixed params.\nparams_for_williams = best_params_initial if USE_GRID_SEARCH else FIXED_PARAMS[MODEL_FAMILY]\n\nif MODEL_FAMILY == 'catboost':\n    model_initial = CatBoostRegressor(loss_function='RMSE', verbose=False,\n                                      random_seed=42, **params_for_williams)\nelif MODEL_FAMILY == 'gbr':\n    model_initial = GradientBoostingRegressor(random_state=42, **params_for_williams)\nelif MODEL_FAMILY == 'svr':\n    model_initial = Pipeline([('scaler', StandardScaler()), ('svr', SVR(**{\n        k.replace('svr__', ''): v for k, v in params_for_williams.items()\n    }))])\n\nmodel_initial.fit(X_all, y_all_log)\ny_pred_full = pd.Series(model_initial.predict(X_all), index=X_all.index)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 8. Applicability Domain: Williams Plot (Leverage Method)\n\nThe Williams plot identifies records that are either:\n- **High-leverage**: structurally unusual relative to the training set (leverage h_ii > h*)\n- **High-residual**: poorly predicted regardless of structural novelty (|standardised residual| > 2)\n\nRecords outside both thresholds are excluded before final model training.  \nLeverage is computed from the Hat matrix H = X(X'X)^-1 X', where X includes an  \nintercept column. The Moore-Penrose pseudoinverse is used in place of the standard  \ninverse to handle rank-deficient matrices (common when PCA retains many components).\n\n**Threshold**: h* = 2p/n, where p = number of features + 1 (intercept), n = number of samples.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# Build Hat matrix using pseudoinverse (robust to rank-deficient PCA feature sets)\nX_np = X_all.to_numpy()\nX_int = np.hstack([np.ones((X_np.shape[0], 1)), X_np])   # add intercept column\nH = X_int @ np.linalg.pinv(X_int.T @ X_int) @ X_int.T   # Hat matrix\n\nleverage_values = np.diag(H)\n\n# Residuals from the full-dataset model predictions\nresiduals = y_all_log.to_numpy() - y_pred_full.to_numpy()\n\n# Standardised residuals: e_i / (RMSE * sqrt(1 - h_ii))\nsigma_hat = cv_metrics_initial['RMSE'].mean()\nstd_residuals = residuals / (sigma_hat * np.sqrt(np.clip(1 - leverage_values, 1e-10, None)))\n\n# Leverage threshold h* = 2p/n\np = X_int.shape[1]\nn = X_int.shape[0]\nh_star = 2 * p / n\n\nleverage_df = pd.DataFrame({\n    'y_true'      : y_all_log.values,\n    'y_pred'      : y_pred_full.values,\n    'leverage'    : leverage_values,\n    'std_residual': std_residuals,\n}, index=X_all.index)\n\n# Flag out-of-domain records\nflagged_mask = (\n    (np.abs(leverage_df['std_residual']) > STD_RESID_THRESHOLD) |\n    (leverage_df['leverage'] > h_star)\n)\n\nprint(f\"Leverage threshold h*   : {h_star:.4f}  (p={p}, n={n})\")\nprint(f\"Std. residual threshold : +/-{STD_RESID_THRESHOLD}\")\nprint(f\"Records flagged         : {flagged_mask.sum()} of {n}\")\nprint(f\"Records retained        : {(~flagged_mask).sum()}\")\n\n# Williams plot\nfig, ax = plt.subplots(figsize=(8, 6))\nax.scatter(\n    leverage_df.loc[~flagged_mask, 'leverage'],\n    leverage_df.loc[~flagged_mask, 'std_residual'],\n    s=30, alpha=0.6, label='Within domain'\n)\nax.scatter(\n    leverage_df.loc[flagged_mask, 'leverage'],\n    leverage_df.loc[flagged_mask, 'std_residual'],\n    s=60, marker='o', facecolors='none', edgecolors='purple', linewidths=1.2,\n    label='Out of domain (flagged)'\n)\nax.axhline( STD_RESID_THRESHOLD, color='red',   linestyle='--', linewidth=1, label=f'|std. residual| = {STD_RESID_THRESHOLD}')\nax.axhline(-STD_RESID_THRESHOLD, color='red',   linestyle='--', linewidth=1)\nax.axvline(h_star,               color='green', linestyle='--', linewidth=1, label=f'h* = {h_star:.3f}')\nax.set_xlabel('Leverage (h_ii)', fontsize=12)\nax.set_ylabel('Standardised Residual', fontsize=12)\nax.set_title(f'Williams Plot - {GAS_TARGET} Permeability', fontsize=13, fontweight='bold')\nax.legend()\nax.grid(True, linestyle=':', alpha=0.5)\nplt.tight_layout()\nplt.savefig(f'williams_plot_{GAS_TARGET}.png', dpi=600, bbox_inches='tight')\nplt.show()"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 9. Final Cross-Validation (Within-Domain Records)\n\nThe model is retrained and evaluated on the leverage-filtered dataset using the same\n`run_cv` function and settings as Section 7. Results reported here represent\nwithin-applicability-domain performance and are the numbers cited in the paper.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# Remove flagged records\nflagged_indices = leverage_df.index[flagged_mask]\nX_clean = X_all.drop(index=flagged_indices)\ny_clean = y_all_log.drop(index=flagged_indices)\n\nprint(f\"Dataset size before leverage filtering : {len(X_all)}\")\nprint(f\"Records removed                        : {len(flagged_indices)}\")\nprint(f\"Dataset size after leverage filtering  : {len(X_clean)}\")\nprint()\n\nprint(f\"Final CV: {GAS_TARGET}, model={MODEL_FAMILY}, \"\n      f\"grid_search={USE_GRID_SEARCH}, n={len(y_clean)}\")\nprint(\"-\" * 60)\ncv_metrics_clean, _, best_params_clean = run_cv(\n    X_clean, y_clean, MODEL_FAMILY, USE_GRID_SEARCH, N_FOLDS, CV_RANDOM_STATE\n)"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 10. Train Final Model on Full Within-Domain Dataset\n\nThe final model is trained on all within-domain records using the best parameters\nfound during grid search (if `USE_GRID_SEARCH = True`) or the fixed parameters\n(if `USE_GRID_SEARCH = False`). In-sample metrics are reported alongside the\nCV estimates from Section 9.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# Use grid-search best params if available, otherwise fall back to fixed params\nparams_final = best_params_clean if USE_GRID_SEARCH else FIXED_PARAMS[MODEL_FAMILY]\n\nif MODEL_FAMILY == 'catboost':\n    final_model = CatBoostRegressor(loss_function='RMSE', verbose=False,\n                                    random_seed=42, **params_final)\nelif MODEL_FAMILY == 'gbr':\n    final_model = GradientBoostingRegressor(random_state=42, **params_final)\nelif MODEL_FAMILY == 'svr':\n    final_model = Pipeline([('scaler', StandardScaler()), ('svr', SVR(**{\n        k.replace('svr__', ''): v for k, v in params_final.items()\n    }))])\n\nfinal_model.fit(X_clean, y_clean)\n\nprint(f\"Final model trained: {MODEL_FAMILY}\")\nprint(f\"Parameters used    : {params_final}\")\nprint(f\"Training samples   : {len(X_clean)}\")\n\n# In-sample predictions (log scale)\ny_pred_log = pd.Series(final_model.predict(X_clean), index=X_clean.index)\n\n# Metrics on log scale\nr2_log   = r2_score(y_clean, y_pred_log)\nrmse_log = np.sqrt(mean_squared_error(y_clean, y_pred_log))\nmae_log  = mean_absolute_error(y_clean, y_pred_log)\n\n# Metrics on original (Barrer) scale\ny_true_orig = 10 ** y_clean\ny_pred_orig = 10 ** y_pred_log\nr2_orig    = r2_score(y_true_orig, y_pred_orig)\nrmse_orig  = np.sqrt(mean_squared_error(y_true_orig, y_pred_orig))\nmae_orig   = mean_absolute_error(y_true_orig, y_pred_orig)\nmed_factor = float(np.median(10 ** np.abs(y_pred_log - y_clean)))\n\nprint(f\"\\n===== In-sample metrics: {GAS_TARGET} permeability =====\")\nprint(f\"  Log scale  : R2={r2_log:.3f}  RMSE={rmse_log:.3f}  MAE={mae_log:.3f}\")\nprint(f\"  Orig scale : R2={r2_orig:.3f}  RMSE={rmse_orig:.3f}  MAE={mae_orig:.3f}\")\nprint(f\"  Median error factor : {med_factor:.3f}\")\nprint()\nprint(\"Note: use CV metrics from the Final CV section for held-out performance estimates.\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 11. Parity and Residual Plots\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "fig, axes = plt.subplots(1, 2, figsize=(13, 6))\n\n# \u2500\u2500 Parity plot (log scale) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nax = axes[0]\nax.scatter(y_clean, y_pred_log, s=20, alpha=0.5)\nlims = [min(y_clean.min(), y_pred_log.min()) - 0.2,\n        max(y_clean.max(), y_pred_log.max()) + 0.2]\nax.plot(lims, lims, 'k--', linewidth=1, label='Parity')\nax.set_xlim(lims); ax.set_ylim(lims)\nax.set_xlabel(f'Actual log10({GAS_TARGET}) [Barrer]', fontsize=12)\nax.set_ylabel(f'Predicted log10({GAS_TARGET}) [Barrer]', fontsize=12)\nax.set_title('Parity Plot (log scale)', fontsize=13, fontweight='bold')\nax.text(0.05, 0.92, f'R2 = {r2_log:.3f}', transform=ax.transAxes, fontsize=11)\nax.legend(); ax.grid(True, linestyle=':', alpha=0.5)\n\n# \u2500\u2500 Residual plot (log scale) \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nax = axes[1]\nresiduals_clean = y_clean - y_pred_log\nax.scatter(y_pred_log, residuals_clean, s=20, alpha=0.5)\nax.axhline(0, color='k', linestyle='--', linewidth=1)\nax.set_xlabel(f'Predicted log10({GAS_TARGET}) [Barrer]', fontsize=12)\nax.set_ylabel('Residual (actual - predicted)', fontsize=12)\nax.set_title('Residual Plot (log scale)', fontsize=13, fontweight='bold')\nax.grid(True, linestyle=':', alpha=0.5)\n\nplt.tight_layout()\nplt.savefig(f'parity_residual_{GAS_TARGET}.png', dpi=600, bbox_inches='tight')\nplt.show()"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 12. Save Model Bundle\n\nThe saved bundle contains everything needed to make predictions on new spectra  \nwithout re-running this notebook.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "bundle = {\n    'final_model'              : final_model,\n    'model_family'             : MODEL_FAMILY,\n    'params_used'              : params_final,\n    'X_columns'                : list(X_clean.columns.astype(str)),\n    'FTIR_MODE'                : FTIR_MODE,\n    'pca_final'                : pca_final if FTIR_MODE == 'pca' else None,\n    'original_ftir_wavenumbers': common_wavelengths,\n    'gas_target'               : GAS_TARGET,\n    'ohe_columns'              : ohe_columns,\n}\n\nfname = (f\"model_{GAS_TARGET}_{MODEL_FAMILY}_{FTIR_MODE}_\"\n         f\"{datetime.now().strftime('%Y%m%d')}.joblib\")\njoblib.dump(bundle, fname)\nprint(f\"Model bundle saved: {fname}\")\nprint(f\"  Model family : {MODEL_FAMILY}\")\nprint(f\"  Parameters   : {params_final}\")"
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": "## 13. Predict Permeability for a New Spectrum\n\nProvide the path to a new ATR-FTIR CSV file (no header; columns: wavenumber, absorbance)  \nand set the polymer-family metadata flags. The cell loads the saved model bundle  \nand returns a permeability prediction in Barrer.\n"
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": "# \u2500\u2500 Inputs \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nMODEL_PATH = fname                       # path to the saved .joblib bundle\nNEW_FTIR_PATH = 'new_polymer.csv'        # replace with your FTIR file path\n\n# Set the relevant one-hot encoded metadata flags to 1; leave others at 0.\n# Example for a semi-crystalline polyamide:\n#   METADATA_FLAGS = {'Category_Semi-crystalline': 1, 'Chemical_Family_Polyamide': 1}\nMETADATA_FLAGS = {}\n\n# \u2500\u2500 Load bundle \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nbundle = joblib.load(MODEL_PATH)\nmodel_infer    = bundle['final_model']\nx_columns      = [str(c) for c in bundle['X_columns']]\npca_infer      = bundle['pca_final']\nwavenumbers    = np.asarray(bundle['original_ftir_wavenumbers'], dtype=float)\nftir_mode_infer = bundle['FTIR_MODE']\n\n# \u2500\u2500 Load and interpolate new spectrum \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nftir_new = pd.read_csv(NEW_FTIR_PATH, names=['wavenumber', 'absorbance'], header=None)\nftir_new = ftir_new.apply(pd.to_numeric, errors='coerce').dropna()\nftir_interp_fn = interp1d(\n    ftir_new['wavenumber'], ftir_new['absorbance'],\n    kind='nearest', fill_value='extrapolate'\n)\nftir_interp = ftir_interp_fn(wavenumbers).reshape(1, -1)\n\n# \u2500\u2500 Apply FTIR feature transform \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nif ftir_mode_infer == 'pca':\n    X_ftir_new = pd.DataFrame(\n        pca_infer.transform(ftir_interp),\n        columns=[str(i) for i in range(pca_infer.n_components_)]\n    )\nelse:\n    X_ftir_new = pd.DataFrame(ftir_interp, columns=combined_df.columns.astype(str))\n\n# \u2500\u2500 Build full input row \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\nX_new = pd.DataFrame(np.zeros((1, len(x_columns))), columns=x_columns, dtype=float)\nfor col in X_ftir_new.columns:\n    if col in X_new.columns:\n        X_new.loc[0, col] = X_ftir_new.loc[0, col]\nfor feat, val in METADATA_FLAGS.items():\n    if feat not in X_new.columns:\n        raise ValueError(f\"Metadata feature '{feat}' not found in trained model columns.\")\n    X_new.loc[0, feat] = val\n\n# \u2500\u2500 Predict \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\npred_log10    = float(model_infer.predict(X_new)[0])\npred_original = 10 ** pred_log10\n\nprint(f\"Predicted log10({bundle['gas_target']}) permeability : {pred_log10:.4f}\")\nprint(f\"Predicted {bundle['gas_target']} permeability (Barrer) : {pred_original:.4f}\")"
  }
 ]
}