XGBOOST_training.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T08:04:42.237208Z",
     "iopub.status.busy": "2023-09-15T08:04:42.237054Z",
     "iopub.status.idle": "2023-09-15T08:04:45.316072Z",
     "shell.execute_reply": "2023-09-15T08:04:45.315299Z",
     "shell.execute_reply.started": "2023-09-15T08:04:42.237208Z"
    }
   },
   "outputs": [],
   "source": [
    "# import pandas for data wrangling\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "#xgboost and sklearn for learning and validation\n",
    "import xgboost as xgb\n",
    "from sklearn.metrics import accuracy_score\n",
    "from xgboost import XGBRegressor\n",
    "from xgboost import plot_importance\n",
    "from sklearn.feature_selection import SelectFromModel\n",
    "from sklearn.metrics import r2_score, mean_squared_error\n",
    "from sklearn.model_selection import StratifiedKFold, cross_val_score\n",
    "\n",
    "# import packages for hyperparameters tuning\n",
    "from hyperopt import tpe, STATUS_OK, Trials, hp, fmin, STATUS_OK, space_eval\n",
    "from hyperopt.pyll.base import scope\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import random"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T08:04:45.438534Z",
     "iopub.status.busy": "2023-09-15T08:04:45.438064Z",
     "iopub.status.idle": "2023-09-15T08:04:45.441803Z",
     "shell.execute_reply": "2023-09-15T08:04:45.441123Z",
     "shell.execute_reply.started": "2023-09-15T08:04:45.438516Z"
    }
   },
   "outputs": [],
   "source": [
    "# Set the seed\n",
    "seed_value = 42\n",
    "np.random.seed(seed_value)\n",
    "random.seed(seed_value)\n",
    "\n",
    "out='results/xgboost/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T08:04:45.442964Z",
     "iopub.status.busy": "2023-09-15T08:04:45.442646Z",
     "iopub.status.idle": "2023-09-15T08:04:45.718451Z",
     "shell.execute_reply": "2023-09-15T08:04:45.717642Z",
     "shell.execute_reply.started": "2023-09-15T08:04:45.442948Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>al</th>\n",
       "      <th>CAS</th>\n",
       "      <th>smiles</th>\n",
       "      <th>ABC</th>\n",
       "      <th>ABCGG</th>\n",
       "      <th>SpAbs_A</th>\n",
       "      <th>SpMax_A</th>\n",
       "      <th>SpDiam_A</th>\n",
       "      <th>SpAD_A</th>\n",
       "      <th>SpMAD_A</th>\n",
       "      <th>...</th>\n",
       "      <th>AMW</th>\n",
       "      <th>WPath</th>\n",
       "      <th>WPol</th>\n",
       "      <th>Zagreb1</th>\n",
       "      <th>Zagreb2</th>\n",
       "      <th>mZagreb1</th>\n",
       "      <th>mZagreb2</th>\n",
       "      <th>δD</th>\n",
       "      <th>δP</th>\n",
       "      <th>δH</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1,1,1,2-Tetrachloroethane</td>\n",
       "      <td>b'630-20-6'</td>\n",
       "      <td>ClCC(Cl)(Cl)Cl</td>\n",
       "      <td>4.012290</td>\n",
       "      <td>4.284816</td>\n",
       "      <td>5.818626</td>\n",
       "      <td>2.074313</td>\n",
       "      <td>4.148627</td>\n",
       "      <td>5.818626</td>\n",
       "      <td>0.969771</td>\n",
       "      <td>...</td>\n",
       "      <td>20.736383</td>\n",
       "      <td>28</td>\n",
       "      <td>3</td>\n",
       "      <td>24.0</td>\n",
       "      <td>22.0</td>\n",
       "      <td>4.312500</td>\n",
       "      <td>1.375000</td>\n",
       "      <td>18.0</td>\n",
       "      <td>4.4</td>\n",
       "      <td>4.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1,1,1-Trichloroethane</td>\n",
       "      <td>b'71-55-6'</td>\n",
       "      <td>CC(Cl)(Cl)Cl</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>...</td>\n",
       "      <td>16.491254</td>\n",
       "      <td>16</td>\n",
       "      <td>0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>4.062500</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>16.8</td>\n",
       "      <td>4.3</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1,1,1-Trifluoroethane</td>\n",
       "      <td>b'420-46-2'</td>\n",
       "      <td>CC(F)(F)F</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.800000</td>\n",
       "      <td>...</td>\n",
       "      <td>10.502336</td>\n",
       "      <td>16</td>\n",
       "      <td>0</td>\n",
       "      <td>20.0</td>\n",
       "      <td>16.0</td>\n",
       "      <td>4.062500</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>14.6</td>\n",
       "      <td>10.0</td>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1,1,2,2-Tetrabromoethane</td>\n",
       "      <td>b'79-27-6'</td>\n",
       "      <td>BrC(Br)C(Br)Br</td>\n",
       "      <td>3.932653</td>\n",
       "      <td>4.244375</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>42.711125</td>\n",
       "      <td>29</td>\n",
       "      <td>4</td>\n",
       "      <td>22.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>4.222222</td>\n",
       "      <td>1.444444</td>\n",
       "      <td>21.0</td>\n",
       "      <td>7.0</td>\n",
       "      <td>8.2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1,1,2,2-Tetrachloroethane</td>\n",
       "      <td>b'79-34-5'</td>\n",
       "      <td>ClC(Cl)C(Cl)Cl</td>\n",
       "      <td>3.932653</td>\n",
       "      <td>4.244375</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>2.000000</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>6.000000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>...</td>\n",
       "      <td>20.736383</td>\n",
       "      <td>29</td>\n",
       "      <td>4</td>\n",
       "      <td>22.0</td>\n",
       "      <td>21.0</td>\n",
       "      <td>4.222222</td>\n",
       "      <td>1.444444</td>\n",
       "      <td>18.8</td>\n",
       "      <td>5.1</td>\n",
       "      <td>5.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1187</th>\n",
       "      <td>Quinine</td>\n",
       "      <td>b'130-95-0'</td>\n",
       "      <td>[H][C@@]1([C@@H](C2=CC=NC3=CC=C(C=C23)OC)O)C[C...</td>\n",
       "      <td>18.958632</td>\n",
       "      <td>14.851657</td>\n",
       "      <td>32.595631</td>\n",
       "      <td>2.494948</td>\n",
       "      <td>4.989896</td>\n",
       "      <td>32.595631</td>\n",
       "      <td>1.358151</td>\n",
       "      <td>...</td>\n",
       "      <td>6.753829</td>\n",
       "      <td>1286</td>\n",
       "      <td>42</td>\n",
       "      <td>132.0</td>\n",
       "      <td>161.0</td>\n",
       "      <td>7.000000</td>\n",
       "      <td>5.333333</td>\n",
       "      <td>19.0</td>\n",
       "      <td>6.6</td>\n",
       "      <td>11.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1188</th>\n",
       "      <td>Sulfur Dioxide</td>\n",
       "      <td>b'9/5/7446'</td>\n",
       "      <td>O=S=O</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>0.942809</td>\n",
       "      <td>...</td>\n",
       "      <td>21.320633</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>2.250000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>15.8</td>\n",
       "      <td>8.4</td>\n",
       "      <td>10.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1189</th>\n",
       "      <td>Thionyl Chloride</td>\n",
       "      <td>b'9/7/7719'</td>\n",
       "      <td>O=S(Cl)Cl</td>\n",
       "      <td>2.449490</td>\n",
       "      <td>2.449490</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>1.732051</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>3.464102</td>\n",
       "      <td>0.866025</td>\n",
       "      <td>...</td>\n",
       "      <td>29.476173</td>\n",
       "      <td>9</td>\n",
       "      <td>0</td>\n",
       "      <td>12.0</td>\n",
       "      <td>9.0</td>\n",
       "      <td>3.111111</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>16.9</td>\n",
       "      <td>6.4</td>\n",
       "      <td>6.1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1190</th>\n",
       "      <td>Triethylene Glycol Monooleyl Ether</td>\n",
       "      <td>b'5274-66-8'</td>\n",
       "      <td>COCCOCCOCCO</td>\n",
       "      <td>7.071068</td>\n",
       "      <td>6.765664</td>\n",
       "      <td>13.191508</td>\n",
       "      <td>1.931852</td>\n",
       "      <td>3.863703</td>\n",
       "      <td>13.191508</td>\n",
       "      <td>1.199228</td>\n",
       "      <td>...</td>\n",
       "      <td>6.077958</td>\n",
       "      <td>220</td>\n",
       "      <td>8</td>\n",
       "      <td>38.0</td>\n",
       "      <td>36.0</td>\n",
       "      <td>4.250000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>16.0</td>\n",
       "      <td>3.1</td>\n",
       "      <td>8.4</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1191</th>\n",
       "      <td>Vinyl Silane</td>\n",
       "      <td>b'7291-09-0'</td>\n",
       "      <td>C=C[SiH3]</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>1.414214</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>2.828427</td>\n",
       "      <td>0.942809</td>\n",
       "      <td>...</td>\n",
       "      <td>6.447097</td>\n",
       "      <td>4</td>\n",
       "      <td>0</td>\n",
       "      <td>6.0</td>\n",
       "      <td>4.0</td>\n",
       "      <td>2.250000</td>\n",
       "      <td>1.000000</td>\n",
       "      <td>15.5</td>\n",
       "      <td>2.6</td>\n",
       "      <td>4.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>1192 rows × 875 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                      al           CAS  \\\n",
       "0              1,1,1,2-Tetrachloroethane   b'630-20-6'   \n",
       "1                  1,1,1-Trichloroethane    b'71-55-6'   \n",
       "2                  1,1,1-Trifluoroethane   b'420-46-2'   \n",
       "3               1,1,2,2-Tetrabromoethane    b'79-27-6'   \n",
       "4              1,1,2,2-Tetrachloroethane    b'79-34-5'   \n",
       "...                                  ...           ...   \n",
       "1187                             Quinine   b'130-95-0'   \n",
       "1188                      Sulfur Dioxide   b'9/5/7446'   \n",
       "1189                    Thionyl Chloride   b'9/7/7719'   \n",
       "1190  Triethylene Glycol Monooleyl Ether  b'5274-66-8'   \n",
       "1191                        Vinyl Silane  b'7291-09-0'   \n",
       "\n",
       "                                                 smiles        ABC      ABCGG  \\\n",
       "0                                        ClCC(Cl)(Cl)Cl   4.012290   4.284816   \n",
       "1                                          CC(Cl)(Cl)Cl   3.464102   3.464102   \n",
       "2                                             CC(F)(F)F   3.464102   3.464102   \n",
       "3                                        BrC(Br)C(Br)Br   3.932653   4.244375   \n",
       "4                                        ClC(Cl)C(Cl)Cl   3.932653   4.244375   \n",
       "...                                                 ...        ...        ...   \n",
       "1187  [H][C@@]1([C@@H](C2=CC=NC3=CC=C(C=C23)OC)O)C[C...  18.958632  14.851657   \n",
       "1188                                              O=S=O   1.414214   1.414214   \n",
       "1189                                          O=S(Cl)Cl   2.449490   2.449490   \n",
       "1190                                        COCCOCCOCCO   7.071068   6.765664   \n",
       "1191                                          C=C[SiH3]   1.414214   1.414214   \n",
       "\n",
       "        SpAbs_A   SpMax_A  SpDiam_A     SpAD_A   SpMAD_A  ...        AMW  \\\n",
       "0      5.818626  2.074313  4.148627   5.818626  0.969771  ...  20.736383   \n",
       "1      4.000000  2.000000  4.000000   4.000000  0.800000  ...  16.491254   \n",
       "2      4.000000  2.000000  4.000000   4.000000  0.800000  ...  10.502336   \n",
       "3      6.000000  2.000000  4.000000   6.000000  1.000000  ...  42.711125   \n",
       "4      6.000000  2.000000  4.000000   6.000000  1.000000  ...  20.736383   \n",
       "...         ...       ...       ...        ...       ...  ...        ...   \n",
       "1187  32.595631  2.494948  4.989896  32.595631  1.358151  ...   6.753829   \n",
       "1188   2.828427  1.414214  2.828427   2.828427  0.942809  ...  21.320633   \n",
       "1189   3.464102  1.732051  3.464102   3.464102  0.866025  ...  29.476173   \n",
       "1190  13.191508  1.931852  3.863703  13.191508  1.199228  ...   6.077958   \n",
       "1191   2.828427  1.414214  2.828427   2.828427  0.942809  ...   6.447097   \n",
       "\n",
       "      WPath  WPol  Zagreb1  Zagreb2  mZagreb1  mZagreb2    δD    δP    δH  \n",
       "0        28     3     24.0     22.0  4.312500  1.375000  18.0   4.4   4.2  \n",
       "1        16     0     20.0     16.0  4.062500  1.000000  16.8   4.3   2.0  \n",
       "2        16     0     20.0     16.0  4.062500  1.000000  14.6  10.0   0.0  \n",
       "3        29     4     22.0     21.0  4.222222  1.444444  21.0   7.0   8.2  \n",
       "4        29     4     22.0     21.0  4.222222  1.444444  18.8   5.1   5.3  \n",
       "...     ...   ...      ...      ...       ...       ...   ...   ...   ...  \n",
       "1187   1286    42    132.0    161.0  7.000000  5.333333  19.0   6.6  11.0  \n",
       "1188      4     0      6.0      4.0  2.250000  1.000000  15.8   8.4  10.0  \n",
       "1189      9     0     12.0      9.0  3.111111  1.000000  16.9   6.4   6.1  \n",
       "1190    220     8     38.0     36.0  4.250000  3.000000  16.0   3.1   8.4  \n",
       "1191      4     0      6.0      4.0  2.250000  1.000000  15.5   2.6   4.0  \n",
       "\n",
       "[1192 rows x 875 columns]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset = pd.read_csv('data/HSP_descriptors_NO_ERRS_ZEROS.csv',index_col=[0] )\n",
    "dataset\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T08:04:45.730732Z",
     "iopub.status.busy": "2023-09-15T08:04:45.730522Z",
     "iopub.status.idle": "2023-09-15T08:04:45.738251Z",
     "shell.execute_reply": "2023-09-15T08:04:45.737576Z",
     "shell.execute_reply.started": "2023-09-15T08:04:45.730715Z"
    }
   },
   "outputs": [],
   "source": [
    "X = dataset.iloc[:,3:-3]\n",
    "y = dataset.iloc[:,-3:]\n",
    "\n",
    "y_D = y['δD']\n",
    "y_P = y['δP']\n",
    "y_H = y['δH']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train, test, val split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T08:04:45.739537Z",
     "iopub.status.busy": "2023-09-15T08:04:45.739264Z",
     "iopub.status.idle": "2023-09-15T08:04:45.749447Z",
     "shell.execute_reply": "2023-09-15T08:04:45.748823Z",
     "shell.execute_reply.started": "2023-09-15T08:04:45.739490Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "X_train_val, X_test, y_train_val, y_test = train_test_split(X, y, test_size=0.2, random_state=seed_value)\n",
    "\n",
    "y_train_val_D,y_train_val_P,y_train_val_H = y_train_val['δD'],y_train_val['δP'],y_train_val['δH']\n",
    "y_test_D,y_test_P,y_test_H = y_test['δD'],y_test['δP'],y_test['δH']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(953, 869)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_train_val.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(239, 869)"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Additional functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T11:39:25.035063Z",
     "iopub.status.busy": "2023-09-15T11:39:25.034337Z",
     "iopub.status.idle": "2023-09-15T11:39:25.038950Z",
     "shell.execute_reply": "2023-09-15T11:39:25.038534Z",
     "shell.execute_reply.started": "2023-09-15T11:39:25.035042Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "/notebooks/NEW_FILTERED/70/\n"
     ]
    }
   ],
   "source": [
    "# number of features(descriptors) to keep\n",
    "n_features=10\n",
    "out_f=out+str(n_features)+'/' #output folder: /results/xgboost/{n_features}/\n",
    "print(out_f)\n",
    "\n",
    "#function to extract the most important features from the xgboost model\n",
    "def to_keep(X,model,n_features=n_features):\n",
    "    # select the best features for the current training set\n",
    "    fi=pd.DataFrame()\n",
    "    fi['columns']=X.columns\n",
    "    fi['importances'] = model.best_estimator_.feature_importances_\n",
    "    fi.sort_values(by='importances',ascending=False,inplace=True)\n",
    "\n",
    "    feets_to_keep = fi[:n_features]['columns']\n",
    "\n",
    "    return feets_to_keep\n",
    "\n",
    "#function to retrieve best model and its statistics\n",
    "def get_best(best,X,y,X_train,y_train,X_test,y_test):\n",
    "\n",
    "    best['max_depth']=int(best['max_depth'])\n",
    "    best_model = XGBRegressor(**best)\n",
    "    best_model.fit(X_train,y_train)\n",
    "\n",
    "    predictions = best_model.predict(X)\n",
    "    test_predictions = best_model.predict(X_test)\n",
    "\n",
    "    results = { \"best\": best,\n",
    "                \"RMSE all\" :mean_squared_error(y,predictions,squared=False),\n",
    "               \"R2 all\":r2_score(y,predictions),\n",
    "               \"RMSE test\": mean_squared_error(y_test,test_predictions,squared=False),\n",
    "               \"R2 test\": r2_score(y_test,test_predictions) }\n",
    "\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-09-15T11:39:25.148893Z",
     "iopub.status.busy": "2023-09-15T11:39:25.148154Z"
    }
   },
   "outputs": [],
   "source": [
    "from sklearn.model_selection import RandomizedSearchCV, KFold\n",
    "\n",
    "# Search space for the weaker XGBoost model used for feature selection\n",
    "\n",
    "model1_params = {\n",
    "    # Maximum depth of the tree, increasing it increases the model complexity.\n",
    "    \"max_depth\": [4,5,6],\n",
    "    # Percentage of columns to be randomly samples for each tree.\n",
    "    \"colsample_bytree\": [0.3,0.5,0.8],\n",
    "    # reg_alpha provides l1 regularization to the weight, higher values result in more conservative models\n",
    "    \"reg_alpha\": [0.1, 1, 10, 100],\n",
    "    # reg_lambda provides l2 regularization to the weight, higher values result in more conservative models\n",
    "    \"reg_lambda\": [0.1, 1, 10, 100]}\n",
    "\n",
    "nsplits1=5\n",
    "# define the outer cross-validation loop\n",
    "outer_cv = KFold(n_splits=nsplits1, shuffle=True)\n",
    "\n",
    "nsplits2=3\n",
    "# define the inner cross-validation loop\n",
    "inner_cv = KFold(n_splits=nsplits2, shuffle=True)\n",
    "\n",
    "# define the weaker models for each HSP parameter\n",
    "fmd=XGBRegressor(tree_method='gpu_hist',seed=seed_value)\n",
    "fmp=XGBRegressor(tree_method='gpu_hist',seed=seed_value)\n",
    "fmh=XGBRegressor(tree_method='gpu_hist',seed=seed_value)\n",
    "\n",
    "# Define the search space for weaker models (inner cross validation)\n",
    "gsd = RandomizedSearchCV(fmd, param_distributions=model1_params, cv=inner_cv, scoring='neg_mean_squared_error')\n",
    "gsp = RandomizedSearchCV(fmp, param_distributions=model1_params, cv=inner_cv, scoring='neg_mean_squared_error')\n",
    "gsh = RandomizedSearchCV(fmh, param_distributions=model1_params, cv=inner_cv, scoring='neg_mean_squared_error')\n",
    "\n",
    "# \"stronger\" model search space\n",
    "\n",
    "space={'max_depth': scope.int(hp.quniform(\"max_depth\", 3, 20, 1)),\n",
    "        'gamma': hp.uniform ('gamma', 0,20),\n",
    "        'reg_alpha' : hp.quniform('reg_alpha', 0,180,1),\n",
    "        'reg_lambda' : hp.uniform('reg_lambda', 0,1),\n",
    "        'colsample_bytree' : hp.uniform('colsample_bytree', 0.5,1), #fraction of columns selected for each tree\n",
    "        'min_child_weight' : hp.quniform('min_child_weight', 0, 20, 1),\n",
    "        'n_estimators': 1000,\n",
    "        'seed':seed_value\n",
    "    }       # 'seed': 0\n",
    "\n",
    "outer_scores = {\"D\":[],\"P\":[],\"H\":[]}\n",
    "\n",
    "max_evals=500 # maximum number of evaluations for the hyperopt search\n",
    "\n",
    "i=1\n",
    "\n",
    "###### TRAINING LOOP #######\n",
    "for train_ix, val_ix in outer_cv.split(X_train_val):\n",
    "    \n",
    "    print(f\" CV:{i}/{nsplits1}\")\n",
    "\n",
    "    # split data into train and val sets\n",
    "    X_train, X_val = X_train_val.iloc[train_ix], X_train_val.iloc[val_ix]\n",
    "    y_train, y_val = y_train_val.iloc[train_ix], y_train_val.iloc[val_ix]\n",
    "\n",
    "    y_train_D = y_train['δD']\n",
    "    y_train_P = y_train['δP']\n",
    "    y_train_H = y_train['δH']\n",
    "\n",
    "    y_val_D = y_val['δD']\n",
    "    y_val_P = y_val['δP']\n",
    "    y_val_H = y_val['δH']\n",
    "\n",
    "    # fit the feature selection grid search on the current training set\n",
    "    gsd.fit(X_train, y_train_D)\n",
    "    gsp.fit(X_train, y_train_P)\n",
    "    gsh.fit(X_train, y_train_H)\n",
    "\n",
    "    # get the most important n_features\n",
    "    D_to_keep = to_keep(X_train,gsd)\n",
    "    P_to_keep = to_keep(X_train,gsp)\n",
    "    H_to_keep = to_keep(X_train,gsh)\n",
    "    \n",
    "    # Use the features for training stronger models\n",
    "    X_D = X_train_val[D_to_keep]\n",
    "    X_P = X_train_val[P_to_keep]\n",
    "    X_H = X_train_val[H_to_keep]\n",
    "\n",
    "    X_train_D,X_train_P,X_train_H = X_train[D_to_keep],X_train[P_to_keep],X_train[H_to_keep]\n",
    "    X_val_D,X_val_P,X_val_H = X_val[D_to_keep], X_val[P_to_keep], X_val[H_to_keep] \n",
    "\n",
    "    ### DEFINE THE HYPEROPT OBJECTIVES FOR EACH PARAMETER\n",
    "    \n",
    "    def hyperparameter_tuning_D(space):\n",
    "\n",
    "        model = xgb.XGBRegressor(**space, tree_method='gpu_hist', eval_metric=\"rmse\",\n",
    "                early_stopping_rounds=100)\n",
    "\n",
    "        #Define evaluation datasets.\n",
    "        evaluation = [(X_train_D, y_train_D), (X_val_D, y_val_D)]\n",
    "\n",
    "        #Fit the model. Define evaluation sets, early_stopping_rounds, and eval_metric.\n",
    "        model.fit(X_train_D, y_train_D,\n",
    "                eval_set=evaluation,verbose=False)\n",
    "\n",
    "        #Obtain prediction and rmse score.\n",
    "        pred = model.predict(X_val_D)\n",
    "        rmse = mean_squared_error(y_val_D, pred, squared=False)\n",
    "        r2 = r2_score(y_val_D,pred)\n",
    "        # print (\"RMSE:\", rmse, 'R2: ', r2)\n",
    "\n",
    "        #Specify what the loss is for each model.\n",
    "        return {'loss':rmse, 'status': STATUS_OK, 'model': model}\n",
    "\n",
    "    def hyperparameter_tuning_P(space):\n",
    "\n",
    "        model = xgb.XGBRegressor(**space, tree_method='gpu_hist', eval_metric=\"rmse\",\n",
    "                early_stopping_rounds=100)\n",
    "\n",
    "        #Define evaluation datasets.\n",
    "        evaluation = [(X_train_P, y_train_P), (X_val_P, y_val_P)]\n",
    "\n",
    "        #Fit the model. Define evaluation sets, early_stopping_rounds, and eval_metric.\n",
    "        model.fit(X_train_P, y_train_P,\n",
    "                eval_set=evaluation,verbose=False)\n",
    "\n",
    "        #Obtain prediction and rmse score.\n",
    "        pred = model.predict(X_val_P)\n",
    "        rmse = mean_squared_error(y_val_P, pred, squared=False)\n",
    "        r2 = r2_score(y_val_P,pred)\n",
    "        #print (\"RMSE:\", rmse, 'R2: ', r2)\n",
    "\n",
    "        #Specify what the loss is for each model.\n",
    "        return {'loss':rmse, 'status': STATUS_OK, 'model': model}\n",
    "\n",
    "    def hyperparameter_tuning_H(space):\n",
    "\n",
    "        model = xgb.XGBRegressor(**space, tree_method='gpu_hist', eval_metric=\"rmse\",\n",
    "                early_stopping_rounds=100)\n",
    "\n",
    "        #Define evaluation datasets.\n",
    "        evaluation = [(X_train_H, y_train_H), (X_val_H, y_val_H)]\n",
    "\n",
    "        #Fit the model. Define evaluation sets, early_stopping_rounds, and eval_metric.\n",
    "        model.fit(X_train_H, y_train_H,\n",
    "                eval_set=evaluation,verbose=False)\n",
    "\n",
    "        #Obtain prediction and rmse score.\n",
    "        pred = model.predict(X_val_H)\n",
    "        rmse = mean_squared_error(y_val_H, pred, squared=False)\n",
    "        r2 = r2_score(y_val_H,pred)\n",
    "        #print (\"RMSE:\", rmse, 'R2: ', r2)\n",
    "\n",
    "        #Specify what the loss is for each model.\n",
    "        return {'loss':rmse, 'status': STATUS_OK, 'model': model}\n",
    "\n",
    "    #### TRAINING + HYPERPARAMETER OPTIMIZATION\n",
    "    \n",
    "    trials = Trials()\n",
    "    best_D = fmin(fn=hyperparameter_tuning_D,\n",
    "                space=space,\n",
    "                algo=tpe.suggest,\n",
    "                max_evals=max_evals,\n",
    "                trials=trials)\n",
    "    trials = Trials()\n",
    "    best_P = fmin(fn=hyperparameter_tuning_P,\n",
    "                space=space,\n",
    "                algo=tpe.suggest,\n",
    "                max_evals=max_evals,\n",
    "                trials=trials)\n",
    "    trials = Trials()\n",
    "    best_H = fmin(fn=hyperparameter_tuning_H,\n",
    "                space=space,\n",
    "                algo=tpe.suggest,\n",
    "                max_evals=max_evals,\n",
    "                trials=trials)\n",
    "\n",
    "    r_D={\"res\":get_best(best_D,X_D,y_train_val_D,X_train_D,y_train_D,X_val_D,y_val_D),\"features\":D_to_keep}  \n",
    "    r_P={\"res\":get_best(best_P,X_P,y_train_val_P,X_train_P,y_train_P,X_val_P,y_val_P),\"features\":P_to_keep}\n",
    "    r_H={\"res\":get_best(best_H,X_H,y_train_val_H,X_train_H,y_train_H,X_val_H,y_val_H),\"features\":H_to_keep}\n",
    "\n",
    "    outer_scores[\"D\"].append(r_D)\n",
    "    outer_scores[\"P\"].append(r_P)\n",
    "    outer_scores[\"H\"].append(r_H)\n",
    "\n",
    "\n",
    "    i+=1\n",
    "# print the cross-validation scores\n",
    "print('Nested cross-validation scores:', outer_scores)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Make dataframes with performance statistics (optional)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r={\"D\":{\"RMSE all\":[],\"RMSE val\":[],\"R2 all\":[],\"R2 val\":[]},\n",
    "  \"P\":{\"RMSE all\":[],\"RMSE val\":[],\"R2 all\":[],\"R2 val\":[]},\n",
    "  \"H\":{\"RMSE all\":[],\"RMSE val\":[],\"R2 all\":[],\"R2 val\":[]}}\n",
    "\n",
    "for par in r.keys():\n",
    "    for i in outer_scores[par]:\n",
    "        r[par][\"RMSE all\"].append(i['res']['RMSE all'])\n",
    "        r[par][\"RMSE val\"].append(i['res']['RMSE test'])\n",
    "        r[par][\"R2 all\"].append(i['res']['R2 all'])\n",
    "        r[par][\"R2 val\"].append(i['res']['R2 test'])\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cross validation results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cv_df = pd.DataFrame(\n",
    "   {\n",
    "      \"<RMSE all>\": [np.mean(r[\"D\"][\"RMSE all\"]),np.mean(r[\"P\"][\"RMSE all\"]),np.mean(r[\"H\"][\"RMSE all\"])],\n",
    "      \"std(RMSE all)\": [np.std(r[\"D\"][\"RMSE all\"]),np.std(r[\"P\"][\"RMSE all\"]),np.std(r[\"H\"][\"RMSE all\"])],\n",
    "      \"<R2 all>\": [np.mean(r[\"D\"][\"R2 all\"]),np.mean(r[\"P\"][\"R2 all\"]),np.mean(r[\"H\"][\"R2 all\"])],\n",
    "       \"std(R2 all)\": [np.std(r[\"D\"][\"R2 all\"]),np.std(r[\"P\"][\"R2 all\"]),np.std(r[\"H\"][\"R2 all\"])],\n",
    "       \"<RMSE val>\": [np.mean(r[\"D\"][\"RMSE val\"]),np.mean(r[\"P\"][\"RMSE val\"]),np.mean(r[\"H\"][\"RMSE val\"])],\n",
    "      \"std(RMSE val)\": [np.std(r[\"D\"][\"RMSE val\"]),np.std(r[\"P\"][\"RMSE val\"]),np.std(r[\"H\"][\"RMSE val\"])],\n",
    "      \"<R2 val>\": [np.mean(r[\"D\"][\"R2 val\"]),np.mean(r[\"P\"][\"R2 val\"]),np.mean(r[\"H\"][\"R2 val\"])],\n",
    "       \"std(R2 val)\":[np.std(r[\"D\"][\"R2 val\"]),np.std(r[\"P\"][\"R2 val\"]),np.std(r[\"H\"][\"R2 val\"])]\n",
    "   },\n",
    "   index=[\"D\", \"P\", \"H\"]\n",
    ")\n",
    "\n",
    "cv_df.to_csv(out_f+f\"CV_{nsplits1}fold_results.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Final training and validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_TRAIN, X_VAL, y_TRAIN, y_VAL = train_test_split(X_train_val, y_train_val, test_size=0.2, random_state=seed_value)\n",
    "\n",
    "y_TRAIN_D,y_TRAIN_P,y_TRAIN_H = y_TRAIN['δD'],y_TRAIN['δP'],y_TRAIN['δH']\n",
    "y_VAL_D,y_VAL_P,y_VAL_H = y_VAL['δD'],y_VAL['δP'],y_VAL['δH']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find the best model and features for the HSP parameter :\n",
    "\n",
    "best_D=outer_scores[\"D\"][np.argmax(r[\"D\"][\"R2 val\"])]['res']['best']\n",
    "D_features=outer_scores[\"D\"][np.argmax(r[\"D\"][\"R2 val\"])]['features']\n",
    "\n",
    "#use those features\n",
    "X_TRAIN_D, X_VAL_D = X_TRAIN[D_features], X_VAL[D_features]\n",
    "\n",
    "# define the model\n",
    "best_D_model = XGBRegressor(**best_D,tree_method='gpu_hist', eval_metric=\"rmse\", early_stopping_rounds=100 )\n",
    "\n",
    "evaluation = [(X_TRAIN_D, y_TRAIN_D), (X_VAL_D, y_VAL_D)]\n",
    "#MODEL FIT\n",
    "best_D_model.fit(X_TRAIN_D, y_TRAIN_D,eval_set=evaluation,verbose=False)\n",
    "\n",
    "\n",
    "# Evaluate the performance and print out the results:\n",
    "\n",
    "predictions_D = best_D_model.predict(X[D_features])\n",
    "test_predictions_D = best_D_model.predict(X_test[D_features])\n",
    "train_val_predictions_D = best_D_model.predict(X_train_val[D_features])\n",
    "\n",
    "best_D_RMSE_all = mean_squared_error(y_D,predictions_D,squared=False)\n",
    "best_D_r2_all = r2_score(y_D,predictions_D)\n",
    "\n",
    "best_D_RMSE_test = mean_squared_error(y_test_D,test_predictions_D,squared=False)\n",
    "best_D_r2_test = r2_score(y_test_D,test_predictions_D)\n",
    "\n",
    "best_D_RMSE_train_val = mean_squared_error(y_train_val_D,train_val_predictions_D,squared=False)\n",
    "best_D_r2_train_val = r2_score(y_train_val_D,train_val_predictions_D)\n",
    "\n",
    "print(\"D:\")\n",
    "print(\"RMSE all:\", best_D_RMSE_all)\n",
    "print(\"R2 all:\", best_D_r2_all)\n",
    "\n",
    "print(\"RMSE test:\", best_D_RMSE_test)\n",
    "print(\"R2 test:\", best_D_r2_test)\n",
    "\n",
    "print(\"RMSE train-val:\", best_D_RMSE_train_val)\n",
    "print(\"R2 train-val:\", best_D_r2_train_val)\n",
    "\n",
    "\n",
    "# Do the same for P and H parameters\n",
    "\n",
    "best_P=outer_scores[\"P\"][np.argmax(r[\"P\"][\"R2 val\"])]['res']['best']\n",
    "P_features=outer_scores[\"P\"][np.argmax(r[\"P\"][\"R2 val\"])]['features']\n",
    "\n",
    "X_TRAIN_P, X_VAL_P = X_TRAIN[P_features], X_VAL[P_features]\n",
    "\n",
    "best_P_model = XGBRegressor(**best_P,tree_method='gpu_hist', eval_metric=\"rmse\", early_stopping_rounds=100 )\n",
    "\n",
    "evaluation = [(X_TRAIN_P, y_TRAIN_P), (X_VAL_P, y_VAL_P)]\n",
    "#MODEL FIT\n",
    "best_P_model.fit(X_TRAIN_P, y_TRAIN_P,eval_set=evaluation,verbose=False)\n",
    "\n",
    "\n",
    "predictions_P = best_P_model.predict(X[P_features])\n",
    "test_predictions_P = best_P_model.predict(X_test[P_features])\n",
    "train_val_predictions_P = best_P_model.predict(X_train_val[P_features])\n",
    "\n",
    "best_P_RMSE_all = mean_squared_error(y_P,predictions_P,squared=False)\n",
    "best_P_r2_all = r2_score(y_P,predictions_P)\n",
    "\n",
    "best_P_RMSE_test = mean_squared_error(y_test_P,test_predictions_P,squared=False)\n",
    "best_P_r2_test = r2_score(y_test_P,test_predictions_P)\n",
    "\n",
    "best_P_RMSE_train_val = mean_squared_error(y_train_val_P,train_val_predictions_P,squared=False)\n",
    "best_P_r2_train_val = r2_score(y_train_val_P,train_val_predictions_P)\n",
    "\n",
    "print(\"P:\")\n",
    "print(\"RMSE all:\", best_P_RMSE_all)\n",
    "print(\"R2 all:\", best_P_r2_all)\n",
    "\n",
    "print(\"RMSE test:\", best_P_RMSE_test)\n",
    "print(\"R2 test:\", best_P_r2_test)\n",
    "\n",
    "print(\"RMSE train-val:\", best_P_RMSE_train_val)\n",
    "print(\"R2 train-val:\", best_P_r2_train_val)\n",
    "\n",
    "\n",
    "\n",
    "best_H=outer_scores[\"H\"][np.argmax(r[\"H\"][\"R2 val\"])]['res']['best']\n",
    "H_features=outer_scores[\"H\"][np.argmax(r[\"H\"][\"R2 val\"])]['features']\n",
    "\n",
    "X_TRAIN_H, X_VAL_H = X_TRAIN[H_features], X_VAL[H_features]\n",
    "best_H_model = XGBRegressor(**best_H,tree_method='gpu_hist', eval_metric=\"rmse\", early_stopping_rounds=100)\n",
    "\n",
    "\n",
    "evaluation = [(X_TRAIN_H, y_TRAIN_H), (X_VAL_H, y_VAL_H)]\n",
    "#MODEL FIT\n",
    "best_H_model.fit(X_TRAIN_H, y_TRAIN_H,eval_set=evaluation,verbose=False)\n",
    "\n",
    "\n",
    "predictions_H = best_H_model.predict(X[H_features])\n",
    "test_predictions_H = best_H_model.predict(X_test[H_features])\n",
    "train_val_predictions_H = best_H_model.predict(X_train_val[H_features])\n",
    "\n",
    "best_H_RMSE_all = mean_squared_error(y_H,predictions_H,squared=False)\n",
    "best_H_r2_all = r2_score(y_H,predictions_H)\n",
    "\n",
    "best_H_RMSE_test = mean_squared_error(y_test_H,test_predictions_H,squared=False)\n",
    "best_H_r2_test = r2_score(y_test_H,test_predictions_H)\n",
    "\n",
    "best_H_RMSE_train_val = mean_squared_error(y_train_val_H,train_val_predictions_H,squared=False)\n",
    "best_H_r2_train_val = r2_score(y_train_val_H,train_val_predictions_H)\n",
    "\n",
    "print(\"D:\")\n",
    "print(\"RMSE all:\", best_H_RMSE_all)\n",
    "print(\"R2 all:\", best_H_r2_all)\n",
    "\n",
    "print(\"RMSE test:\", best_H_RMSE_test)\n",
    "print(\"R2 test:\", best_H_r2_test)\n",
    "\n",
    "print(\"RMSE train-val:\", best_H_RMSE_train_val)\n",
    "print(\"R2 train-val:\", best_H_r2_train_val)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save the models:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "best_P_model.save_model(out_f+f\"MODEL_P_XGBOOST.json\")\n",
    "best_D_model.save_model(out_f+f\"MODEL_D_XGBOOST.json\")\n",
    "best_H_model.save_model(out_f+f\"MODEL_H_XGBOOST.json\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optional dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "best_cv_df = pd.DataFrame(\n",
    "   {\n",
    "      \"RMSE all\": [best_D_RMSE_all,best_P_RMSE_all,best_H_RMSE_all],\n",
    "      \n",
    "      \"R2 all\": [best_D_r2_all,best_P_r2_all,best_H_r2_all],\n",
    "\n",
    "      \"RMSE train-val\": [best_D_RMSE_train_val,best_P_RMSE_train_val,best_H_RMSE_train_val],\n",
    "\n",
    "      \"R2 train-val\": [best_D_r2_train_val,best_P_r2_train_val,best_H_r2_train_val],\n",
    "       \n",
    "      \"RMSE test\": [best_D_RMSE_test,best_P_RMSE_test,best_H_RMSE_test],\n",
    "      \n",
    "      \"R2 test\": [best_D_r2_test,best_P_r2_test,best_H_r2_test]},\n",
    "       \n",
    "   index=[\"D\", \"P\", \"H\"] )\n",
    "\n",
    "best_cv_df.to_csv(out_f+\"best_CV_results.csv\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = {'D':{'true':y_train_val_D.values,'pred':train_val_predictions_D},\n",
    "           'P':{'true':y_train_val_P.values,'pred':train_val_predictions_P},\n",
    "           'H':{'true':y_train_val_H.values,'pred':train_val_predictions_H}\n",
    "           }\n",
    "y_test = {'D':{'true':y_test_D.values,'pred':test_predictions_D},\n",
    "           'P':{'true':y_test_P.values,'pred':test_predictions_P},\n",
    "           'H':{'true':y_test_H.values,'pred':test_predictions_H}\n",
    "           }\n",
    "y_train_df=pd.DataFrame.from_dict({(i,j): y_train[i][j] \n",
    "                           for i in y_train.keys() \n",
    "                           for j in y_train[i].keys()},\n",
    "                       orient='index').T\n",
    "y_train_df.to_csv(out_f+\"y_train_XGBOOST_xgboost_nfeats{n_features}.csv\")\n",
    "\n",
    "y_test_df=pd.DataFrame.from_dict({(i,j): y_test[i][j] \n",
    "                           for i in y_test.keys() \n",
    "                           for j in y_test[i].keys()},\n",
    "                       orient='index').T\n",
    "y_test_df.to_csv(out_f+f\"y_test_XGBOOST_xgboost_nfeats{n_features}.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Feature importances:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "D_IMP=pd.DataFrame({'Name': D_to_keep,'Importance':best_D_model.feature_importances_})\n",
    "P_IMP=pd.DataFrame({'Name': P_to_keep,'Importance':best_P_model.feature_importances_})\n",
    "H_IMP=pd.DataFrame({'Name': H_to_keep,'Importance':best_H_model.feature_importances_})\n",
    "\n",
    "D_IMP.to_csv(out_f+f'D_XGBOOST_NONONE_xgboost_feature_importances_nfeats{n_features}.csv')\n",
    "P_IMP.to_csv(out_f+f'P_XGBOOST_NONONE_xgboost_feature_importances_nfeats{n_features}.csv')\n",
    "H_IMP.to_csv(out_f+f'H_XGBOOST_NONONE_xgboost_feature_importances_nfeats{n_features}.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "y_test_means_stds_df = pd.DataFrame(\n",
    "   {\n",
    "      \"<TRUE>\": [ np.mean(y_test['D']['true']),np.mean(y_test['P']['true']),np.mean(y_test['H']['true']) ],\n",
    "      \"std(True)\": [np.std(y_test['D']['true']),np.std(y_test['P']['true']),np.std(y_test['H']['true'])],\n",
    "      \"<PREDICTED>\": [ np.mean(y_test['D']['pred']),np.mean(y_test['P']['pred']),np.mean(y_test['H']['pred']) ],\n",
    "       \"std(PREDICTED)\": [np.std(y_test['D']['pred']),np.std(y_test['P']['pred']),np.std(y_test['H']['pred'])]\n",
    "       },\n",
    "   index=[\"D\", \"P\", \"H\"]\n",
    ")\n",
    "\n",
    "y_test_means_stds_df.to_csv(out_f+'y_test_means_stds.csv')\n",
    "\n",
    "y_train_means_stds_df = pd.DataFrame(\n",
    "   {\n",
    "      \"<TRUE>\": [ np.mean(y_train['D']['true']),np.mean(y_train['P']['true']),np.mean(y_train['H']['true']) ],\n",
    "      \"std(True)\": [np.std(y_train['D']['true']),np.std(y_train['P']['true']),np.std(y_train['H']['true'])],\n",
    "      \"<PREDICTED>\": [ np.mean(y_train['D']['pred']),np.mean(y_train['P']['pred']),np.mean(y_train['H']['pred']) ],\n",
    "       \"std(PREDICTED)\": [np.std(y_train['D']['pred']),np.std(y_train['P']['pred']),np.std(y_train['H']['pred'])]\n",
    "       },\n",
    "   index=[\"D\", \"P\", \"H\"]\n",
    ")\n",
    "\n",
    "y_train_means_stds_df.to_csv(out_f+'y_train_means_stds.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Optional figure:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ynp=pd.DataFrame.to_numpy(y)\n",
    "predictions = [predictions_D,predictions_P,predictions_H]\n",
    "sns.set_style('whitegrid')\n",
    "sns.set_context(\"paper\")\n",
    "fig,axs = plt.subplots(nrows=1,ncols=3,figsize=(20,5))\n",
    "fig.suptitle('Hansen solubility parameters XGBOOST predictions and true values')\n",
    "colors = ['b','y','g']\n",
    "for i in range(3):\n",
    "    sns.scatterplot(x=ynp[:,i],y=predictions[i],ax=axs[i],color=colors[i])\n",
    "    sns.lineplot(x=ynp[:,i],y=ynp[:,i],ax=axs[i],color='r',alpha=0.8, label=f'{y.columns[i]}-{y.columns[i]}')\n",
    "    axs[i].set_xlabel(y.columns[i])\n",
    "    axs[i].set_ylabel(f'{y.columns[i]} XGBOOST prediction')\n",
    "\n",
    "#fig.savefig(out_f+f'XGBOOST_xgboost_NONONE_1_{n_features}.png')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  },
  "vscode": {
   "interpreter": {
    "hash": "85dfa1d5ff96de0fad6e949b2d0fc07a90dc20502f07b7bbd90aaffa3f3c36c4"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}