Boosting.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "accompanied-joining",
   "metadata": {},
   "source": [
    "## Boosting with applications in astronomy\n",
    "\n",
    "This is a short introduction to Boosting algorithms with an application in astronomy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "demanding-harvest",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Basic packages and plotting\n",
    "import numpy as np \n",
    "import scipy.io as sio\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from matplotlib import rcParams\n",
    "import pandas as pd\n",
    "#import seaborn as sns\n",
    "rcParams['font.family'] = 'serif'\n",
    "\n",
    "# Adjust rc parameters to make plots pretty\n",
    "def plot_pretty(dpi=200, fontsize=9):\n",
    "    \n",
    "    import matplotlib.pyplot as plt\n",
    "\n",
    "    plt.rc(\"savefig\", dpi=dpi)       # dpi resolution of saved image files\n",
    "    plt.rc('text', usetex=True)      # use LaTeX to process labels\n",
    "    plt.rc('font', size=fontsize)    # fontsize\n",
    "    plt.rc('xtick', direction='in')  # make axes ticks point inward\n",
    "    plt.rc('ytick', direction='in')\n",
    "    plt.rc('xtick.major', pad=10) \n",
    "    plt.rc('xtick.minor', pad=5)\n",
    "    plt.rc('ytick.major', pad=10) \n",
    "    plt.rc('ytick.minor', pad=5)\n",
    "    plt.rc('lines', dotted_pattern = [0.5, 1.1]) # fix dotted lines\n",
    "\n",
    "    return\n",
    "\n",
    "plot_pretty()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "conditional-retailer",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Skikit-learn and XGBoost\n",
    "from sklearn.ensemble import AdaBoostClassifier #Adaptive Boosting\n",
    "from sklearn.ensemble import GradientBoostingClassifier #Gradient Boosting\n",
    "#from xgboost import XGBClassifier #Extreme Gradient Boosting\n",
    "from sklearn.linear_model import LogisticRegression # Linear regression for comparison\n",
    "\n",
    "# Metrics\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import precision_score\n",
    "from sklearn.metrics import recall_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "\n",
    "# Preprocessing\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import train_test_split"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "human-cricket",
   "metadata": {},
   "source": [
    "### Astronomy-related Example\n",
    "\n",
    "Import a catalog of SDSS objects (Stars, Galaxies, QSOs).\n",
    "\n",
    "Specifically, a feature matrix (`X_feat`) that contains magnitudes $u,g,r,i,z$, colors (all the possible differences of magnitudes), and ratios of magnitudes.\n",
    "\n",
    "We also \n",
    "\n",
    "- Stars $\\to y = 0$\n",
    "\n",
    "- Galaxies $\\to y = 1$\n",
    "\n",
    "- Quasars $\\to y = 2$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "painful-yield",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import feature matrix\n",
    "X_feat = np.load(\"Feature_matrix.npy\")\n",
    "# Import labels\n",
    "y_lab = np.load(\"Labels.npy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "minute-reminder",
   "metadata": {},
   "source": [
    "Now, keep only stars and galaxies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "statewide-arbitration",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_feat_bn = X_feat[y_lab!=2] #For binary classification, Stars and Galaxies only\n",
    "y_bn = y_lab[y_lab!=2] # Labels for binary classification"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "advisory-outreach",
   "metadata": {},
   "source": [
    "**Split into Training-Validation-Test sets**\n",
    "\n",
    "Using `sklearn`'s `train_test_split` this can be done in two steps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "emerging-association",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train_bn, X_valtest_bn, y_train_bn, y_valtest_bn = train_test_split(\n",
    "    X_feat_bn, y_bn, test_size=0.3, random_state=42)\n",
    "X_val_bn, X_test_bn, y_val_bn, y_test_bn = train_test_split(\n",
    "    X_valtest_bn, y_valtest_bn, test_size=0.5, random_state=42)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "communist-target",
   "metadata": {},
   "source": [
    "**Rescale features**: We rescale all the features, in such a way that all have mean 0, and standard deviation 1 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "charitable-tokyo",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaler = StandardScaler()\n",
    "\n",
    "scaler.fit(X_train_bn)\n",
    "X_train_bn = scaler.transform(X_train_bn) # Train, rescaled\n",
    "X_val_bn = scaler.transform(X_val_bn) # Validation, rescaled\n",
    "X_test_bn = scaler.transform(X_test_bn) # Test, rescaled"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "separate-dinner",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "14000\n"
     ]
    }
   ],
   "source": [
    "print(len(X_train_bn))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "collect-marble",
   "metadata": {},
   "source": [
    "### **Baseline: Logistic Regression**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "controversial-modem",
   "metadata": {},
   "outputs": [],
   "source": [
    "LR = LogisticRegression(penalty='none', max_iter=10000)\n",
    "LR.fit(X_train_bn, y_train_bn)\n",
    "y_pred_LR = LR.predict(X_test_bn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "declared-albuquerque",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.8156666666666667\n",
      "Precision: 0.8235294117647058\n",
      "Recall: 0.8133159268929504\n"
     ]
    }
   ],
   "source": [
    "print('Accuracy:',accuracy_score(y_test_bn,y_pred_LR))\n",
    "print('Precision:',precision_score(y_test_bn,y_pred_LR))\n",
    "print('Recall:',recall_score(y_test_bn,y_pred_LR))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ambient-fashion",
   "metadata": {},
   "source": [
    "### AdaBoost\n",
    "\n",
    "Adaptive Boosting\n",
    "\n",
    "![Adaboost](Adaboost.png)\n",
    "(not mine)\n",
    "\n",
    "![Boost_example](Boost_example.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "medical-least",
   "metadata": {},
   "outputs": [],
   "source": [
    "ADA = AdaBoostClassifier(random_state=0)\n",
    "ADA.fit(X_train_bn, y_train_bn)\n",
    "y_pred_ADA = ADA.predict(X_test_bn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "collectible-ideal",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.8923333333333333\n",
      "Precision: 0.8847867600254615\n",
      "Recall: 0.9073107049608355\n"
     ]
    }
   ],
   "source": [
    "print('Accuracy:',accuracy_score(y_test_bn,y_pred_ADA))\n",
    "print('Precision:',precision_score(y_test_bn,y_pred_ADA))\n",
    "print('Recall:',recall_score(y_test_bn,y_pred_ADA))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "third-investor",
   "metadata": {},
   "source": [
    "### Gradient Boosting (classifier)\n",
    "\n",
    "![Boost_1](Boost_1.png)\n",
    "![Boost_2](Boost_2.png)\n",
    "![Boost_3](Boost_3.png)\n",
    "![Boost_4](Boost_4.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "essential-montana",
   "metadata": {},
   "outputs": [],
   "source": [
    "GRAD = GradientBoostingClassifier(n_estimators=100, learning_rate=1., max_depth=5, random_state=0)\n",
    "GRAD.fit(X_train_bn, y_train_bn)\n",
    "y_pred_GRAD = GRAD.predict(X_test_bn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "elementary-principal",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy: 0.916\n",
      "Precision: 0.91343669250646\n",
      "Recall: 0.922976501305483\n"
     ]
    }
   ],
   "source": [
    "print('Accuracy:',accuracy_score(y_test_bn,y_pred_GRAD))\n",
    "print('Precision:',precision_score(y_test_bn,y_pred_GRAD))\n",
    "print('Recall:',recall_score(y_test_bn,y_pred_GRAD))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "modular-updating",
   "metadata": {},
   "source": [
    "### XGBoost (Extreme Gradient Boosting)\n",
    "\n",
    "![XGBoost](XGboost.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "dying-address",
   "metadata": {},
   "outputs": [
    {
     "ename": "XGBoostError",
     "evalue": "\nXGBoost Library (libxgboost.dylib) could not be loaded.\nLikely causes:\n  * OpenMP runtime is not installed\n    - vcomp140.dll or libgomp-1.dll for Windows\n    - libomp.dylib for Mac OSX\n    - libgomp.so for Linux and other UNIX-like OSes\n    Mac OSX users: Run `brew install libomp` to install OpenMP runtime.\n\n  * You are running 32-bit Python on a 64-bit OS\n\nError message(s): ['dlopen(/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/lib/libxgboost.dylib, 6): Library not loaded: /usr/local/opt/libomp/lib/libomp.dylib\\n  Referenced from: /Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/lib/libxgboost.dylib\\n  Reason: image not found']\n",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mXGBoostError\u001b[0m                              Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-14-b1ab54ea0f05>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0mxgboost\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mXGBClassifier\u001b[0m \u001b[0;31m#Extreme Gradient Boosting\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/__init__.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      7\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      8\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 9\u001b[0;31m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m\u001b[0mcore\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mDMatrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDeviceQuantileDMatrix\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mBooster\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mDataIter\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     10\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m\u001b[0mtraining\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mtrain\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcv\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;32mfrom\u001b[0m \u001b[0;34m.\u001b[0m \u001b[0;32mimport\u001b[0m \u001b[0mrabit\u001b[0m  \u001b[0;31m# noqa\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/core.py\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m    201\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m \u001b[0;31m# load the XGBoost library globally\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 203\u001b[0;31m \u001b[0m_LIB\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_load_lib\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    204\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    205\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/core.py\u001b[0m in \u001b[0;36m_load_lib\u001b[0;34m()\u001b[0m\n\u001b[1;32m    179\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mlib_success\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    180\u001b[0m         \u001b[0mlibname\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbasename\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlib_paths\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 181\u001b[0;31m         raise XGBoostError(\n\u001b[0m\u001b[1;32m    182\u001b[0m             f\"\"\"\n\u001b[1;32m    183\u001b[0m \u001b[0mXGBoost\u001b[0m \u001b[0mLibrary\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mlibname\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0mcould\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mbe\u001b[0m \u001b[0mloaded\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mXGBoostError\u001b[0m: \nXGBoost Library (libxgboost.dylib) could not be loaded.\nLikely causes:\n  * OpenMP runtime is not installed\n    - vcomp140.dll or libgomp-1.dll for Windows\n    - libomp.dylib for Mac OSX\n    - libgomp.so for Linux and other UNIX-like OSes\n    Mac OSX users: Run `brew install libomp` to install OpenMP runtime.\n\n  * You are running 32-bit Python on a 64-bit OS\n\nError message(s): ['dlopen(/Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/lib/libxgboost.dylib, 6): Library not loaded: /usr/local/opt/libomp/lib/libomp.dylib\\n  Referenced from: /Library/Frameworks/Python.framework/Versions/3.9/lib/python3.9/site-packages/xgboost/lib/libxgboost.dylib\\n  Reason: image not found']\n"
     ]
    }
   ],
   "source": [
    "from xgboost import XGBClassifier #Extreme Gradient Boosting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "sorted-unknown",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'XGBClassifier' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-15-888cc4ce996d>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mXGB\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mXGBClassifier\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmax_depth\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlearning_rate\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1.0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn_estimators\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m100\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrandom_state\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0mXGB\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train_bn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my_train_bn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0my_pred_XGB\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mXGB\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpredict\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_test_bn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'XGBClassifier' is not defined"
     ]
    }
   ],
   "source": [
    "XGB = XGBClassifier(max_depth=1, learning_rate=1.0, n_estimators=100, random_state=0)\n",
    "XGB.fit(X_train_bn, y_train_bn)\n",
    "y_pred_XGB = XGB.predict(X_test_bn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tutorial-quantum",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "equipped-riverside",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "pleasant-equivalent",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}