feature_selection.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-21T10:08:06.340700Z",
     "iopub.status.busy": "2023-05-21T10:08:06.340287Z",
     "iopub.status.idle": "2023-05-21T10:08:06.947033Z",
     "shell.execute_reply": "2023-05-21T10:08:06.946264Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd \n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "plt.style.use(\"default\")\n",
    "plt.rc(\"text\", usetex=True)\n",
    "plt.rc(\"font\", family=\"cm\")\n",
    "plt.rcParams[\"grid.color\"] = (0.5, 0.5, 0.5, 0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-21T10:08:06.950870Z",
     "iopub.status.busy": "2023-05-21T10:08:06.950383Z",
     "iopub.status.idle": "2023-05-21T10:09:58.859289Z",
     "shell.execute_reply": "2023-05-21T10:09:58.858307Z"
    }
   },
   "outputs": [],
   "source": [
    "X_dataframe = pd.read_csv(\"/mnt/ferracci/features_dataframe_new.csv.gz\")\n",
    "X = np.load(\"/mnt/ferracci/features_new.npz\", allow_pickle=True)['a']\n",
    "y = np.array(pd.read_csv(\"/mnt/ferracci/targets_dataframe_new.csv.gz\")[\"Qedep\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-21T10:09:58.863407Z",
     "iopub.status.busy": "2023-05-21T10:09:58.863159Z",
     "iopub.status.idle": "2023-05-21T10:15:54.609628Z",
     "shell.execute_reply": "2023-05-21T10:15:54.608907Z"
    }
   },
   "outputs": [],
   "source": [
    "import xgboost as xgb \n",
    "\n",
    "dm_train = xgb.DMatrix(X, label=y)\n",
    "\n",
    "params = {\"max_depth\": 8, \n",
    "          \"min_child_weight\": 1, \n",
    "          \"eta\": 0.08, \n",
    "          \"subsample\": 0.8, \n",
    "          \"colsample_bytree\": 1, \n",
    "          \"objective\": \"reg:squarederror\", \n",
    "          \"eval_metric\": \"mape\",\n",
    "          \"tree_method\": \"gpu_hist\",\n",
    "}\n",
    "\n",
    "xgb_cv = xgb.cv(params, dtrain=dm_train, num_boost_round=1000, nfold=5, verbose_eval=False,\n",
    "                early_stopping_rounds=5, as_pandas=True)\n",
    "\n",
    "mape_all_features = xgb_cv[\"test-mape-mean\"].min()\n",
    "std_mape_all_features = xgb_cv.loc[xgb_cv[\"test-mape-mean\"].idxmin(), \"test-mape-std\"]\n",
    "print(f\"({100 * mape_all_features:.4f} +/- {100 * std_mape_all_features:.4f})%\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-21T10:15:54.645735Z",
     "iopub.status.busy": "2023-05-21T10:15:54.645258Z",
     "iopub.status.idle": "2023-05-21T10:15:54.654431Z",
     "shell.execute_reply": "2023-05-21T10:15:54.653938Z"
    }
   },
   "outputs": [],
   "source": [
    "# since the features are highly correlated, not all of them are necessary\n",
    "def feature_selection(features_dataframe, mape_all_features, std_mape_all_features):\n",
    "    \"\"\"\n",
    "    Adds one by one the feature that maximizes stepwise perfomance gain, until the mape is comparable (inside one\n",
    "    standard deviation) of the mape obtained using all features.\n",
    "\n",
    "    Parameters:\n",
    "        features_dataframe (pd.DataFrame): pandas dataframe with columns corresponding to the computed features\n",
    "        mape_all_features (float): mape obtained using all features\n",
    "        std_mape_all_features (float): standard deviation on the mape obtained using all features\n",
    "\n",
    "    Returns:\n",
    "        selected_features_names (list): list of strings representing the names of selected features\n",
    "        trailing_best_mape (list): best mape obtained at each iteration of feature selection\n",
    "        trailing_std_best_mape (list): standard deviation of best mape for each iteration of feature selection\n",
    "    \"\"\"\n",
    "    selected_features = np.empty((len(features_dataframe), 0)) # empty column vector\n",
    "    selected_features_names = []\n",
    "    trailing_best_mape = []\n",
    "    trailing_std_best_mape = []\n",
    "\n",
    "    # first iteration: we need to find the feature that provides the best performance on its own\n",
    "    best_mape = 1 \n",
    "    for name, feature in features_dataframe.items():\n",
    "        feature = np.array(feature).reshape(-1, 1)\n",
    "        # add feature to the stack\n",
    "        selected_features = np.column_stack((selected_features, feature))\n",
    "        dm_train = xgb.DMatrix(selected_features, label=y)\n",
    "        xgb_cv = xgb.cv(params, dtrain=dm_train, num_boost_round=1000, nfold=5, verbose_eval=False,\n",
    "                    early_stopping_rounds=5, as_pandas=True)\n",
    "        mape = xgb_cv[\"test-mape-mean\"].min()\n",
    "        if (mape < best_mape):\n",
    "            best_mape = mape\n",
    "            std_mape = xgb_cv.loc[xgb_cv[\"test-mape-mean\"].idxmin(), \"test-mape-std\"]\n",
    "            best_feature = feature\n",
    "            best_feature_name = name\n",
    "        # pop back feature after loop\n",
    "        selected_features = np.delete(selected_features, -1, axis=1)\n",
    "\n",
    "    selected_features = np.column_stack((selected_features, best_feature))\n",
    "    selected_features_names.append(best_feature_name)\n",
    "    trailing_best_mape.append(best_mape)\n",
    "    trailing_std_best_mape.append(std_mape)\n",
    "    \n",
    "    # print first selected feature and trailing mape:\n",
    "    print(f\"Adding feature: {selected_features_names[0]} | \" +\n",
    "          f\"MAPE: ({100 * trailing_best_mape[0]:.2f} +/- {100 * trailing_std_best_mape[0]:.2f})%\") \n",
    "\n",
    "    # successive iterations: we need to find the feature that provides the best stepwise performance gain\n",
    "    while np.abs(trailing_best_mape[-1] - mape_all_features) > std_mape_all_features:\n",
    "        epsilon = 1e-6 # tolerance on performace gain\n",
    "        for name, feature in features_dataframe.drop(selected_features_names, axis=1).items():\n",
    "            feature = np.array(feature).reshape(-1, 1)\n",
    "            # add feature to the stack\n",
    "            selected_features = np.column_stack((selected_features, feature))\n",
    "            dm_train = xgb.DMatrix(selected_features, label=y)\n",
    "            xgb_cv = xgb.cv(params, dtrain=dm_train, num_boost_round=1000, nfold=5, verbose_eval=False,\n",
    "                        early_stopping_rounds=5, as_pandas=True)\n",
    "            mape = xgb_cv[\"test-mape-mean\"].min()\n",
    "            if (trailing_best_mape[-1] - mape) > epsilon:\n",
    "                epsilon = trailing_best_mape[-1] - mape # updated highest performange gain\n",
    "                best_mape = mape\n",
    "                std_mape = xgb_cv.loc[xgb_cv[\"test-mape-mean\"].idxmin(), \"test-mape-std\"]\n",
    "                best_feature = feature\n",
    "                best_feature_name = name\n",
    "            # pop back feature after loop\n",
    "            selected_features = np.delete(selected_features, -1, axis=1)\n",
    "\n",
    "        selected_features = np.column_stack((selected_features, best_feature))\n",
    "        selected_features_names.append(best_feature_name)\n",
    "        trailing_best_mape.append(best_mape)\n",
    "        trailing_std_best_mape.append(std_mape)\n",
    "\n",
    "        # print selected feature and trailing mape:\n",
    "        print(f\"Adding feature: {selected_features_names[-1]} | \" +\n",
    "              f\"MAPE: ({100 * trailing_best_mape[-1]:.4f} +/- {100 * trailing_std_best_mape[-1]:.4f})%\") \n",
    "\n",
    "    return selected_features_names, trailing_best_mape, trailing_std_best_mape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-21T10:15:54.657252Z",
     "iopub.status.busy": "2023-05-21T10:15:54.656869Z",
     "iopub.status.idle": "2023-05-22T16:09:34.729679Z",
     "shell.execute_reply": "2023-05-22T16:09:34.728038Z"
    }
   },
   "outputs": [],
   "source": [
    "selected_features_names, trailing_best_mape, trailing_std_best_mape = feature_selection(X_dataframe,\n",
    "                                                                                        mape_all_features,\n",
    "                                                                                        std_mape_all_features)\n",
    "\n",
    "with open('/home/ferracci/new_dataset/features_list.txt', 'w') as f:\n",
    "    f.write(str(selected_features_names))\n",
    "    f.write('\\n')\n",
    "    f.write(str(trailing_best_mape))\n",
    "    f.write('\\n')\n",
    "    f.write(str(trailing_std_best_mape))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2023-05-22T16:09:34.736728Z",
     "iopub.status.busy": "2023-05-22T16:09:34.735930Z",
     "iopub.status.idle": "2023-05-22T16:09:36.982078Z",
     "shell.execute_reply": "2023-05-22T16:09:36.980699Z"
    }
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(11,4), width_ratios=(1,2), dpi=150)\n",
    "\n",
    "with open('/home/ferracci/new_dataset/features_list.txt', 'r') as f:\n",
    "    file_content = f.read()\n",
    "\n",
    "features_list = file_content.split('\\n')\n",
    "selected_features_names = eval(features_list[0])\n",
    "trailing_best_mape = eval(features_list[1])\n",
    "trailing_std_best_mape = eval(features_list[2])\n",
    "\n",
    "ax = axes[0]\n",
    "ax.errorbar([0, 1, 2, 3], np.array(trailing_best_mape[:4])*100, yerr=np.array(trailing_std_best_mape[:4])*100,\n",
    "            marker=\"o\", markersize=3, capsize=5., color=\"k\")\n",
    "ax.axhline(mape_all_features*100, linewidth=1., linestyle=\"--\", dashes=(5, 5))\n",
    "ax.fill_between(x=range(-1, len(selected_features_names)+1), y1=(mape_all_features-std_mape_all_features)*100, \n",
    "                y2=(mape_all_features+std_mape_all_features)*100, alpha=0.3)\n",
    "\n",
    "ax.set_xlim([-0.2, 3.2])\n",
    "\n",
    "ax.set_xticks([0, 1, 2, 3])\n",
    "ax.tick_params(axis=\"both\", which=\"major\", labelsize=12)\n",
    "ax.tick_params(axis=\"both\", which='minor', labelsize=12)\n",
    "\n",
    "ax.set_xticklabels([\"AccumCharge\", \"$pe_{std}$\", \"$\\\\rho_{cc}$\", \"$ht_{95\\%-90\\%}$\"], \n",
    "                   rotation=-45, ha=\"left\", rotation_mode=\"anchor\")\n",
    "\n",
    "ax.grid()\n",
    "\n",
    "ax = axes[1]\n",
    "ax.errorbar(selected_features_names[3:], np.array(trailing_best_mape[3:])*100, yerr=np.array(trailing_std_best_mape[3:])*100,\n",
    "            marker=\"o\", markersize=3, capsize=3., color=\"k\")\n",
    "ax.axhline(mape_all_features*100, linewidth=1., linestyle=\"--\", dashes=(5, 5))\n",
    "ax.fill_between(x=range(-1, len(selected_features_names[2:])+1), y1=(mape_all_features-std_mape_all_features)*100, \n",
    "                y2=(mape_all_features+std_mape_all_features)*100, alpha=0.3)\n",
    "\n",
    "ax.set_xlim([-0.2, len(selected_features_names[3:])-0.8])\n",
    "\n",
    "ax.tick_params(axis=\"both\", which=\"major\", labelsize=12)\n",
    "ax.tick_params(axis=\"both\", which='minor', labelsize=12)\n",
    "\n",
    "ax.set_xticklabels([\"$ht_{95\\%-90\\%}$\", \"nPMTs\", \"$z_{cc}$\", \"$ht_{entropy}$\", \"$ht_{5\\%-2\\%}$\", \"$ht_{kurtosis}$\", \"$pe_{mean}$\", \n",
    "                    \"$ht_{10\\%-5\\%}$\", \"$pe_{15\\%}$\", \"$ht_{80\\%-75\\%}$\"], \n",
    "                   rotation=-45, ha=\"left\", rotation_mode=\"anchor\")\n",
    "\n",
    "ax.grid()\n",
    "\n",
    "fig.supxlabel(\"Added feature\", fontsize=15, y=-0.15)\n",
    "fig.supylabel(\"MAPE, \\%\", fontsize=15, x=0.07)\n",
    "\n",
    "fig.savefig(\"/home/ferracci/new_dataset/images/BDT_feature_selection.png\", dpi=300, bbox_inches=\"tight\", pad_inches=0.2);"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "deep_learning",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}