07_RF_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Forest model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is our first attempt with the Random Forest model. It is a classical Random Forest model. In the grid search we try to optimize the accuracy "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import plot_confusion_matrix\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to limitations in the computation power we do not do a grid search as exhaustive as would be desirable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = '/data/ml_course/05_Capstone_project/rf_data/'\n",
    "model_name = 'RF_model.pickle'\n",
    "test_size = 0.1 # 10% of points used for testing\n",
    "val_size = 0.1  # 10% of points (after removing test points) used for validation\n",
    "n_splits = 5  # number of splits in folding\n",
    "n_estimators = np.r_[np.arange(1, 21, 2), np.arange(20, 60, 10)]  # number of trees\n",
    "\n",
    "features = ['HRV_norm', 'HRV_norm_text', 'IR_108', 'IR_108_text', 'WV_062-IR_108', 'WV_062-IR_108_text']\n",
    "nfeatures = len(features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "202009\r"
     ]
    }
   ],
   "source": [
    "years = ['2018', '2019', '2020']\n",
    "months = ['04', '05', '06', '07', '08', '09']\n",
    "\n",
    "for year in years:\n",
    "    for month in months:\n",
    "        print(year+month, end=\"\\r\", flush=True)\n",
    "        fname = datapath+year+month+'_data.npz'\n",
    "        if not os.path.isfile(fname):\n",
    "            continue\n",
    "        \n",
    "        X = None\n",
    "        with np.load(fname, allow_pickle=False) as npz_file:\n",
    "            # Load the arrays\n",
    "            if X is None:\n",
    "                X = npz_file['features']\n",
    "                y = npz_file['targets']\n",
    "            else:\n",
    "                X = np.concatenate((X, npz_file['features']), axis=0)\n",
    "                y = np.concatenate((y, npz_file['targets']), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix shape: (7435724, 6)\n",
      "target matrix shape: (7435724,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix shape:', X.shape)\n",
    "print('target matrix shape:', y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove NaNs\n",
    "# NaNs can appear in texture fields.\n",
    "# They could safely be considered 0s but we remove them\n",
    "ind_row, ind_col = np.where(np.isnan(X))\n",
    "X = np.delete(X, ind_row, axis=0)\n",
    "y = np.delete(y, ind_row, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix shape: (7432290, 6)\n",
      "target matrix shape: (7432290,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix shape:', X.shape)\n",
    "print('target matrix shape:', y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1], dtype=int8)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Group together 0 (No POH was computed) and 1 (POH below 90%)\n",
    "np.unique(y)\n",
    "y[y == 1] = 0\n",
    "y[y == 2] = 1\n",
    "np.unique(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "npixels: 7432290\n",
      "npixels with hail: 61171\n",
      "% of pixels with hail: 0.8230437725115678\n"
     ]
    }
   ],
   "source": [
    "npixels = y.size\n",
    "npixels_hail = y[y==1].size\n",
    "\n",
    "print('npixels:', npixels)\n",
    "print('npixels with hail:', npixels_hail)\n",
    "print('% of pixels with hail:', npixels_hail/npixels*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is highly imbalanced toward non-hail pixels!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Split data into test and training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split data\n",
    "X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=test_size, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix for training shape: (6689061, 6)\n",
      "target matrix for training shape: (6689061,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix for training shape:', X_tr.shape)\n",
    "print('target matrix for training shape:', y_tr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare data folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss = ShuffleSplit(n_splits=n_splits, test_size=val_size, random_state=0)\n",
    "ss.get_n_splits()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 5 folds for each of 14 candidates, totalling 70 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.\n",
      "[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed: 77.8min\n",
      "[Parallel(n_jobs=-1)]: Done  70 out of  70 | elapsed: 256.1min finished\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=ShuffleSplit(n_splits=5, random_state=0, test_size=0.1, train_size=None),\n",
       "             estimator=Pipeline(steps=[('scaler', None),\n",
       "                                       ('RF',\n",
       "                                        RandomForestClassifier(random_state=0))]),\n",
       "             n_jobs=-1,\n",
       "             param_grid={'RF__n_estimators': array([ 1,  3,  5,  7,  9, 11, 13, 15, 17, 19, 20, 30, 40, 50])},\n",
       "             return_train_score=True, verbose=True)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create pipeline\n",
    "# Weighted random forest: modifies the weight by penalizing more heavily the under-represented class\n",
    "# balanced: modifies the weight according to the overall distribution\n",
    "# balanced_subsample: modies the weight according to the distribution in the bootstrapped subsamples\n",
    "pipe_RF = Pipeline([\n",
    "    ('scaler', None),\n",
    "    ('RF', RandomForestClassifier(random_state=0, max_depth=None, class_weight=None))\n",
    "])\n",
    "\n",
    "# Create cross-validation object\n",
    "grid_RF = {'RF__n_estimators': n_estimators}\n",
    "\n",
    "grid_RF_cv = GridSearchCV(pipe_RF, grid_RF, cv=ss, return_train_score=True, verbose=True, n_jobs=-1)\n",
    "\n",
    "# Fit estimator\n",
    "grid_RF_cv.fit(X_tr, y_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect results in a DataFrame\n",
    "RF_results = pd.DataFrame({\n",
    "    'trees': grid_RF_cv.cv_results_['param_RF__n_estimators'],\n",
    "    'mean_tr': grid_RF_cv.cv_results_['mean_train_score'],\n",
    "    'mean_te': grid_RF_cv.cv_results_['mean_test_score'],\n",
    "    'std_te': grid_RF_cv.cv_results_['std_test_score']\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEHCAYAAACncpHfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAiQ0lEQVR4nO3de5hcVZ3u8e+vqm+5d8htIJ2QiAESQkBo2yjIDHKQgAwR1BFE8UEwh5Gb58zoYXAQxms4cBxA8GE4gILcdMQcM4JA5EQRjgFCSCA3IQ+JoScggXR3Qld31+13/ti7OpWiQ1enu9N0rffzPP1U7WutFXG9e+3L2ubuiIhIeBJDXQARERkaCgARkUApAEREAqUAEBEJlAJARCRQCgARkUBVlbOSmS0AbgSSwO3uvrhk+XjgTuAQoBP4kruvjZddDnwZMOB/u/sNJdv+I3AdMMnd33y3ckycONFnzJhRTpFFRCT23HPPvenuk0rn9xoAZpYEbgFOBpqBZ81sqbuvL1rtSmC1u59pZofH659kZnOJGv8mIA08YmYPufvL8b6nxfvdWk4lZsyYwcqVK8tZVUREYmb2557ml3MKqAnY5O6vuHsaeABYWLLOHOBxAHffCMwwsynAbGCFu6fcPQv8HjizaLt/Bb4O6Gk0EZH9rJwAmAq8WjTdHM8rtgY4C8DMmoCDgQZgLXCCmU0ws5HAacC0eL0zgP909zXv9uNmtsjMVprZyu3bt5dRXBERKUc5AWA9zCs9Yl8MjDez1cClwPNA1t03ANcCy4BHiIIiG4fBN4Bv9vbj7n6buze6e+OkSe84hSUiIvuonIvAzcRH7bEGYFvxCu6+EzgfwMwM2Bz/4e53AHfEy74X7+8QYCawJlqdBmCVmTW5++v9qI+IiJSpnAB4FphlZjOB/wTOBj5XvIKZ1QOp+BrBhcATcShgZpPd/Q0zm050mujD7t4CTC7afgvQ2NtdQCIiMnB6DQB3z5rZJcCjRLeB3unu68zsonj5rUQXe+82sxywHrigaBcPmtkEIANcHDf+IiIyxGw4DQfd2Njoug1URKRvzOw5d28snV/Wg2Cy/3Rlc7R1ZNjZkaGtIxt/ZrrnVVclGD+ymvEjaxg/qib6HFlN/cgakomerteLiPRMATAEOjM5fv3Cazy27nV2tKfZ2bm7ke/M5Pdpn2Ywtq6aA0bVUD+ymgNG1lA/soaxI6qoTiaoShhVCSOZSFCVNJLxdFXCSMbLkwmjOmnd+yn8ja6tIr5YLyIVRAGwH216Yxf3Pr2VB59rZmdnlmkHjKChfiTvmziacSOqGTeymrF1VYwbUc3YEdXRvPhv7IhqxtZVk83n2dGepjWVYUd7mpZUmpb2NC2pDK2pNDviz9d3drLx9V3s7MiQzTu5vJPJ59mXM341ycQegVD4mzAq6oWMqYv+Myrs2+O7hLunu+dH8u50ZfN0ZXJ0ZnJ0ZHJ0ZvJ07vGZozMbTXekc5w696+49KRZ/fjXF5FSCoBB1pXN8cja17n36a08s3kH1UljwdwDOfdD0/nQzAP6fGRdQ4KRNVU0jN+38uTzvkcg5HJF07k82bzT1pFhR3sXb70dBcxb7Wl2vJ1mR3uaHak0r7ak2PF2ml1d2X0rRImEwYjqJHXxX211grqqJCNqktRVJ8jmnf+17CWmjh/BWcc0DMhviogCYNBsebOd+5/Zyr8/18yO9jQHTxjJFacezqePbWDi6NohK1ciYdTE1wpGkOzXvrqyOVpTGd7uynY/LVgItN3T8Wc8pzBdW52IGvyqJNVJe9cgzObynHv701y55EXmHDSWw/9qbL/KLSIR3QU0gDK5PMvW/4X7nt7Kk5veJJkwPj5nCp/70HSOO2QiCV2k3Wdv7OrkEzc9yejaKpZechxj6qqHukgiw4buAhpkv3iumWsf2cj2XV1MrR/BP5x8KH/3wWlMGVs31EWrCJPH1HHzOR/gc7c/zdd/8QI/OvcYXZiuMO5OOpcnnc2TyTnpbPS9MC+dy5PJ5cnno6tMeXfco8+8R9sXTxeWe9F0Pj7gzbuTzxetQ/E679xXz/suzNs9nS/Zlzt7Ke/ubT2ejr7vOZ33aFt35+IT388RB40b0H9zBcAAeGX721z5yxc5YupYrv3Ukfz1oZN1S+Yg+ND7JvD1Uw7j+7/ZyB1PbubCj75vqIs0LOXzcUNbaFizUcOazubpKvpeaHAL83c3zLl4mRfN33NfXT3sO13yWbxNYX/DmVl06jNhRsIMM3r8TFh0qrT4M2HRSVIzI5Eo2ke834QZqXRuwMusAOgnd+fqpeuorUrwb184lsljdMQ/mBad8D6e+3MLi3+zkaOn1dM444ChLlKv3KML7YWGtCub2/09kyedy9GViRrNrkzJ8myhAc69Y/vdjWnhyDm31yPn4gY3mx/YhrYmmaA6adRUJbr/qpMJapIJagvfqxKMrqvq/l6bLFqvaJvaqkTR/pK79xnvvzqZKGlMCw1poTHd3aBCvE6ipIG1PRvp0n3xjn3v2Xh3b0fJ9DDskSoA+uk3a1/nDy+/ydV/O0eN/35gZlz/d0dxxg+f5OL7VvHQZR9914vqhaPdqIGNG9oeG9Xd83pdXtRYR/vO7WX73fsYiDa3OmlRo1qdpKak4Sw0qiOqk4ytq4qXRRfYdzeqiT0a6e59lC4r3Xdh/+9osKPyDMeGTyIKgH5o78ryrf9Yz+wDx/KF+QcPdXGCMbaumh+deyxn/ugpFt78FAeMqtmjsS5ufAfitIIZ1FYlqI2PSGurCg1isvv7mLqq7unaqgS11Yk9GuvSecX7KOyzdL+F3yxsp5sIZKApAPrhpsdf5vWdndxy7geoSpbzagUZKHMOGsvNnzuGu/+4haqE7dmQVieoSUYN594a2tqSBry4oS1tqKsS736bqshwpQDYRy//ZRd3PLmZzxzbwLEHv/fPQ1eik+dM4eQ5U4a6GCLDlg5b94G7c9Wv1jKqtoorTj18qIsjIrJPFAD7YOmabax4ZQdfO+UwJgzhU70iIv2hAOijnZ0ZvvPQBuY1jOOcpulDXRwRkX2mawB99K/LXuLNt7u4/bxGPewlIsOaegB9sH7bTu76f1s4p2k6R02rH+riiIj0iwKgTPm8881fraV+ZA1fP+WwoS6OiEi/KQDK9OCqZlb+uYUrFhxO/ciaoS6OiEi/KQDK0JbKsPg3Gzlmej2fPlYvJBGRyqCLwGW47rGNtKTS3H1Bkx7HF5GKoR5AL15sbuPep7dy3odnDPhY3CIiQ0kB8C7yeeeff7WWCaNq+e8fP3SoiyMiMqAUAO/igWdfZc2rrXzjE4czVq8gFJEKowDYix3taf7noxv50MwD+OTRU4e6OCIiA04BsBfX/mYjuzqzfPuTczUUsIhUJAVAD7a1dvCzla9y/kdmcOiUMUNdHBGRQaEA6MFrbR0AHD9r4hCXRERk8CgAetDSngFgvJ74FZEKpgDoQWtHFAD1I3Xnj4hULgVAD1pTaQCN+SMiFU0B0IPWVIaEwZhajZQhIpVLAdCDllSa+pE1GvdHRCpaWQFgZgvM7E9mtsnMruhh+XgzW2JmL5jZM2Y2t2jZ5Wa21szWmdlXi+Z/O15/tZk9ZmYHDUiNBkBrR0bn/0Wk4vUaAGaWBG4BTgXmAOeY2ZyS1a4EVrv7POA84MZ427nAl4Em4CjgdDObFW9znbvPc/ejgV8D3+x/dQZGaypN/QgFgIhUtnJ6AE3AJnd/xd3TwAPAwpJ15gCPA7j7RmCGmU0BZgMr3D3l7lng98CZ8Xo7i7YfBXi/ajKAWtozugVURCpeOQEwFXi1aLo5nldsDXAWgJk1AQcDDcBa4AQzm2BmI4HTgGmFjczsu2b2KnAue+kBmNkiM1tpZiu3b99eXq36qa0jwzidAhKRCldOAPR0JbT0aH0xMN7MVgOXAs8DWXffAFwLLAMeIQqKbPdO3L/h7tOAe4FLevpxd7/N3RvdvXHSpEllFLf/WlJp9QBEpOKVEwDNFB21Ex3Zbytewd13uvv58fn884BJwOZ42R3ufoy7nwDsAF7u4TfuAz7V9+IPvK5sjlQ6p2sAIlLxygmAZ4FZZjbTzGqAs4GlxSuYWX28DOBC4InCOX4zmxx/Tic6TXR/PD2raBdnABv7U5GB0paKnwIepR6AiFS2Xp90cvesmV0CPAokgTvdfZ2ZXRQvv5XoYu/dZpYD1gMXFO3iQTObAGSAi929JZ6/2MwOA/LAn4GLBqpS/dE9DIR6ACJS4cp61NXdHwYeLpl3a9H3PwKzSreLl310L/PfE6d8SrW0R8NA6BqAiFQ6PQlcoiWlgeBEJAwKgBJtHYWB4BQAIlLZFAAlCj0AnQISkUqnACjRmspQnTRG1iSHuigiIoNKAVCiNR4JVC+CF5FKpwAo0ZrK6BZQEQmCAqCEhoEQkVAoAEq06V0AIhIIBUCJ6G1gCgARqXwKgCLuTktK7wIQkTAoAIp0ZvKks3m9C0BEgqAAKNKS0jhAIhIOBUCR1pRGAhWRcCgAirSmCuMAqQcgIpVPAVCk8C6A8aPUAxCRyqcAKFK4BlA/Qj0AEal8CoAirXoXgIgERAFQpDWVpq46QV21RgIVkcqnACiih8BEJCQKgCKtqQzjdAuoiARCAVCkVSOBikhAFABFWjsyugVURIKhACjSmkozTreAikggFAAxd6c1lWG8bgEVkUAoAGJvd2XJ5l3PAIhIMBQAsd0PgekUkIiEQQEQ00igIhIaBUCs+10Ao9QDEJEwKABi3SOB6hqAiARCARArvAtAt4GKSCgUADGNBCoioVEAxFpSaUbXVlGd1D+JiIRBrV2sNZXR0b+IBEUBEGtNpRUAIhKUsgLAzBaY2Z/MbJOZXdHD8vFmtsTMXjCzZ8xsbtGyy81srZmtM7OvFs2/zsw2xtssMbP6gajQvtK7AEQkNL0GgJklgVuAU4E5wDlmNqdktSuB1e4+DzgPuDHedi7wZaAJOAo43cxmxdssA+bG27wE/FP/q7Pv2jr0LgARCUs5PYAmYJO7v+LuaeABYGHJOnOAxwHcfSMww8ymALOBFe6ecvcs8HvgzHi9x+J5ACuAhn7Xph9a9C4AEQlMOQEwFXi1aLo5nldsDXAWgJk1AQcTNehrgRPMbIKZjQROA6b18BtfAn7T04+b2SIzW2lmK7dv315Gcfsun3faOjQSqIiEpZwAsB7mecn0YmC8ma0GLgWeB7LuvgG4luh0zyNEQZEt3tDMvhHPu7enH3f329y90d0bJ02aVEZx+25nZwZ3GKcegIgEpKqMdZrZ86i9AdhWvIK77wTOBzAzAzbHf7j7HcAd8bLvxfsjnv4icDpwkruXhsp+05LSMBAiEp5yegDPArPMbKaZ1QBnA0uLVzCz+ngZwIXAE3EoYGaT48/pRKeJ7o+nFwD/AzjD3VMDUZl9VRgGQreBikhIeu0BuHvWzC4BHgWSwJ3uvs7MLoqX30p0sfduM8sB64ELinbxoJlNADLAxe7eEs+/GagFlkWdBla4+0UDVK8+0bsARCRE5ZwCwt0fBh4umXdr0fc/ArNKt4uXfXQv899ffjEHV2tH3APQbaAiEhA9CQy0tBeuAagHICLhUAAQvQvADMaqByAiAVEAEF0EHltXTTLR0x2vIiKVSQFAYRwgHf2LSFgUAEQ9AD0EJiKhUQAQ3QaqHoCIhEYBQHQbqG4BFZHQKACA1vaMHgITkeAEHwCZXJ5dXVk9AyAiwQk+ANo6CsNA6BSQiIQl+ADQQHAiEioFgAaCE5FABR8AeheAiIQq+ADoPgU0Qj0AEQmLAqBwCmiUegAiEhYFQEeaZMIYU1vWqxFERCpG8AHQkspQP6Ka+K1kIiLBCD4A2lIZ3QIqIkEKPgBaUmndAioiQVIAaCRQEQlU8AHQlkozTreAikiAgg8A9QBEJFRBB0BnJkdHJsf4UeoBiEh4gg6Awkig4/QyGBEJUNAB0No9DpB6ACISnqADoEVDQYtIwIIOAL0LQERCFngA6F0AIhKuoANA7wIQkZAFHQCtHWlqkglGVCeHuigiIvtd2AHQHg0Ep5FARSREYQdAR1q3gIpIsIIOgJZUhnE6/y8igSorAMxsgZn9ycw2mdkVPSwfb2ZLzOwFM3vGzOYWLbvczNaa2Toz+2rR/M/E8/Jm1jggtemj1lRaF4BFJFi9BoCZJYFbgFOBOcA5ZjanZLUrgdXuPg84D7gx3nYu8GWgCTgKON3MZsXbrAXOAp4YgHrsk9ZURi+DF5FgldMDaAI2ufsr7p4GHgAWlqwzB3gcwN03AjPMbAowG1jh7il3zwK/B86M19vg7n8aoHr0mbtHAaCXwYtIoMoJgKnAq0XTzfG8YmuIjuYxsybgYKCB6Cj/BDObYGYjgdOAaf0t9EDoyORI5/LqAYhIsKrKWKeneyS9ZHoxcKOZrQZeBJ4Hsu6+wcyuBZYBbxMFRbYvBTSzRcAigOnTp/dl03elh8BEJHTl9ACa2fOovQHYVryCu+909/Pd/WiiawCTgM3xsjvc/Rh3PwHYAbzclwK6+23u3ujujZMmTerLpu9q9zhA6gGISJjKCYBngVlmNtPMaoCzgaXFK5hZfbwM4ELgCXffGS+bHH9OJzpNdP9AFb4/do8DpB6AiISp11NA7p41s0uAR4EkcKe7rzOzi+LltxJd7L3bzHLAeuCCol08aGYTgAxwsbu3AJjZmcAPiXoLD5nZanc/ZQDr9q70LgARCV051wBw94eBh0vm3Vr0/Y/ArNLt4mUf3cv8JcCSsks6wPQuABEJXbBPAheuAeh1kCISqoADIMOI6iR1GglURAIVbAC0pDK6BVREghZsALR1pHULqIgELdgAaElldAFYRIIWbABEI4GqByAi4Qo4APQuABEJW5AB4O60dugisIiELcgA2NWVJZd3jQQqIkELMgBa2zUOkIhImAHQET0FrIvAIhKyIAOgRSOBioiEGQB6F4CISLABoB6AiEiQAdA9FLRGAhWRgAUZAK2pDGNqq6hKBll9EREg2ABIUz9KR/8iErYwA6Ajo1tARSR4QQZASyqjN4GJSPCCDIA2jQQqIhJmAOhdACIiAQZALu/s7MzoITARCV5wAbCzI4O7ngEQEQkuAAoPgY3XbaAiErjgAqC1Ix4GQu8CEJHAhRcA3QPBqQcgImELMACiHoBuAxWR0AUXAHoXgIhIJLgAaEulMYOxdQoAEQlbcAFQGAYikbChLoqIyJAKMAA0DISICAQYAG0dGghORAQCDICoB6AAEBEJLgBaU3oXgIgIlBkAZrbAzP5kZpvM7Ioelo83syVm9oKZPWNmc4uWXW5ma81snZl9tWj+AWa2zMxejj/HD0iNetGayjBOPQARkd4DwMySwC3AqcAc4Bwzm1Oy2pXAanefB5wH3BhvOxf4MtAEHAWcbmaz4m2uAB5391nA4/H0oMrk8rzdlVUPQESE8noATcAmd3/F3dPAA8DCknXmEDXiuPtGYIaZTQFmAyvcPeXuWeD3wJnxNguBu+LvdwGf7E9FytGqh8BERLqVEwBTgVeLppvjecXWAGcBmFkTcDDQAKwFTjCzCWY2EjgNmBZvM8XdXwOIPyf39ONmtsjMVprZyu3bt5dXq73YPQ6QegAiIuUEQE9PTHnJ9GJgvJmtBi4Fngey7r4BuBZYBjxCFBTZvhTQ3W9z90Z3b5w0aVJfNn2H3SOBqgcgIlJVxjrN7D5qh+jIflvxCu6+EzgfwMwM2Bz/4e53AHfEy74X7w/gL2Z2oLu/ZmYHAm/0ox5laWmP3wWgHoCISFk9gGeBWWY208xqgLOBpcUrmFl9vAzgQuCJOBQws8nx53Si00T3x+stBb4Yf/8i8Kv+VKQc3T0AXQMQEem9B+DuWTO7BHgUSAJ3uvs6M7soXn4r0cXeu80sB6wHLijaxYNmNgHIABe7e0s8fzHwczO7ANgKfGagKrU3eheAiMhu5ZwCwt0fBh4umXdr0fc/ArNKt4uXfXQv898CTiq7pAOgNZWhKmGMri2r2iIiFS2oJ4FbUhnqR1YTXaYQEQlbUAHQ1pHWLaAiIrGgAqClPaNbQEVEYmEFQEo9ABGRgqACoK0jo6GgRURiQQVA1ANQAIiIQEAB0JnJ0ZnJ6xSQiEgsmADQSKAiInsKJwA6NA6QiEixYAKgpV0jgYqIFAsmAPQuABGRPYUTAPFIoONHqQcgIgIBBUBLoQcwQj0AEREIKADaUhlqqxKMqEkOdVFERN4TghkXWQ+BiQwfmUyG5uZmOjs7h7oow0pdXR0NDQ1UV5fX1gUTAK2pjG4BFRkmmpubGTNmDDNmzNDw7WVyd9566y2am5uZOXNmWdsEcwqoNZVhnG4BFRkWOjs7mTBhghr/PjAzJkyY0KdeUzAB0JJKqwcgMoyo8e+7vv6bBRMArR0ZXQMQESkSRAC4O616F4CIDCM33HADqVRqUH8jiABIpXNkcq53AYjIsLE/AiCIu4C6HwJTAIgMO//yH+tYv23ngO5zzkFjufpvj9jr8i1btrBgwQKOP/54VqxYwVFHHcX555/P1VdfzRtvvMG9997LEUccwaWXXsqLL75INpvlmmuuYeHChWzZsoUvfOELtLe3A3DzzTfzkY98hN/97ndcc801TJw4kbVr13Lsscdyzz339Hje/qabbmLbtm2ceOKJTJw4keXLl/PYY49x9dVX09XVxSGHHMKPf/xjRo8e3a9/hyB6ALuHgtYpIBEpz6ZNm7j88st54YUX2LhxI/fddx9PPvkk119/Pd/73vf47ne/y8c+9jGeffZZli9fzte+9jXa29uZPHkyy5YtY9WqVfzsZz/jsssu697n888/zw033MD69et55ZVXeOqpp3r87csuu4yDDjqI5cuXs3z5ct58802+853v8Nvf/pZVq1bR2NjID37wg37XMYgeQHcA6DZQkWHn3Y7UB9PMmTM58sgjATjiiCM46aSTMDOOPPJItmzZQnNzM0uXLuX6668HoltXt27dykEHHcQll1zC6tWrSSaTvPTSS937bGpqoqGhAYCjjz6aLVu2cPzxx/dalhUrVrB+/XqOO+44ANLpNB/+8If7XccwAqDwLoBR6gGISHlqa2u7vycSie7pRCJBNpslmUzy4IMPcthhh+2x3TXXXMOUKVNYs2YN+Xyeurq6HveZTCbJZrNllcXdOfnkk7n//vv7U6V3COIUUIt6ACIywE455RR++MMf4u5AdHoHoK2tjQMPPJBEIsFPf/pTcrncPu1/zJgx7Nq1C4D58+fz1FNPsWnTJgBSqdQePYt9FUQAtLbrXQAiMrCuuuoqMpkM8+bNY+7cuVx11VUAfOUrX+Guu+5i/vz5vPTSS4waNWqf9r9o0SJOPfVUTjzxRCZNmsRPfvITzjnnHObNm8f8+fPZuHFjv+tghfQaDhobG33lypV93u7bv17PA89sZd23FgxCqURkoG3YsIHZs2cPdTGGpZ7+7czsOXdvLF03iB7AoVNG84l5Bw51MURE3lOCuAj82Q9O57MfnD7UxRAReYczzzyTzZs37zHv2muv5ZRTThn03w4iAERE3quWLFkyZL8dxCkgERl+htP1yfeKvv6bKQBE5D2nrq6Ot956SyHQB4UXwhQ/d9AbnQISkfechoYGmpub2b59+1AXZVgpvBKyXGUFgJktAG4EksDt7r64ZPl44E7gEKAT+JK7r42X/TfgQsCBF4Hz3b3TzI4CbgVGA1uAc919YEd8EpFhqbq6uuzXGsq+6/UUkJklgVuAU4E5wDlmNqdktSuB1e4+DziPKCwws6nAZUCju88lCpCz421uB65w9yOBJcDX+l8dEREpVznXAJqATe7+irungQeAhSXrzAEeB3D3jcAMM5sSL6sCRphZFTAS2BbPPwx4Iv6+DPjUPtdCRET6rJwAmAq8WjTdHM8rtgY4C8DMmoCDgQZ3/0/gemAr8BrQ5u6PxdusBc6Iv38GmNbTj5vZIjNbaWYrdT5QRGTglHMNoKe3DJdeml8M3Ghmq4nO8z8PZONrAwuBmUAr8O9m9nl3vwf4EnCTmX0TWAqke/pxd78NuA3AzLab2Z97Ke9E4M0y6lVpVO+wqN7h6U/dD+5pZjkB0MyeR+cN7D6NA0B88fZ8AIteb7M5/jsF2Ozu2+NlvwQ+AtwTnyr6eDz/UOATvRXE3Sf1to6ZrexpzItKp3qHRfUOz2DUvZxTQM8Cs8xsppnVEF3EXVpSsPp4GUR3/DwRh8JWYL6ZjYyD4SRgQ7zN5PgzAfwz0R1BIiKyn/QaAO6eBS4BHiVqvH/u7uvM7CIzuyhebTawzsw2Et0tdHm87dPAL4BVRKeGEsSnc4juJnoJ2EjUo/jxgNVKRER6NayGgy6HmS2KrxsERfUOi+odnsGoe8UFgIiIlEdjAYmIBEoBICISqIoJADNbYGZ/MrNNZnbFUJdnMJnZnWb2hpmtLZp3gJktM7OX48/xQ1nGwWBm08xsuZltMLN1ZnZ5PL+i625mdWb2jJmtiev9L/H8iq43REPRmNnzZvbreLri6wxgZlvM7EUzW21mK+N5A173igiAMscrqiQ/AUpfcHwF8Li7zyIalqMSQzAL/IO7zwbmAxfH/ztXet27gI+5+1HA0cACM5tP5dcbojsKNxRNh1DnghPd/eiie/8HvO4VEQCUN15RxXD3J4AdJbMXAnfF3+8CPrk/y7Q/uPtr7r4q/r6LqGGYSoXX3SNvx5PV8Z9T4fU2swaiB0RvL5pd0XXuxYDXvVICoJzxiirdFHd/DaKGEpg8xOUZVGY2A/gA8DQB1D0+FbIaeANYFj9jU+n1vgH4OpAvmlfpdS5w4DEze87MFsXzBrzulfJCmHLGK5IKYWajgQeBr7r7zugh88rm7jngaDOrB5aY2dwhLtKgMrPTgTfc/Tkz+5shLs5QOM7dt8UjJiyLH7IdcJXSA+h1vKIA/MXMDgSIP98Y4vIMCjOrJmr873X3X8azg6g7gLu3Ar8jugZUyfU+DjjDzLYQndL9mJndQ2XXuZu7b4s/3yB6X0oTg1D3SgmAXscrCsBS4Ivx9y8CvxrCsgyKeDypO4AN7v6DokUVXXczmxQf+WNmI4D/QjSESsXW293/yd0b3H0G0f+f/6+7f54KrnOBmY0yszGF70SDZq5lEOpeMU8Cm9lpROcMk8Cd7v7doS3R4DGz+4G/IRoe9i/A1cD/AX4OTCcahO8z7l56oXhYM7PjgT8QjStVOC98JdF1gIqtu5nNI7rolyQ6aPu5u3/LzCZQwfUuiE8B/aO7nx5Cnc3sfURH/RCdpr/P3b87GHWvmAAQEZG+qZRTQCIi0kcKABGRQCkAREQCpQAQEQmUAkBEJFAKAJEexO+5/spQl0NkMCkARHpWD7wjAOKRZ0UqggJApGeLgUPi8difjd9DcB/wYjww23Xx/BfM7L8WNjKzrxXNL4zbP8rMHorH819rZp8dqkqJFKuUweBEBtoVwFx3Pzp+EvWheHpzPDpjm7t/0MxqgafM7DFgVvzXRDRA4VIzOwGYBGxz908AmNm4/V8dkXdSAIiU5xl33xx//zgwz8w+HU+PI2r4Px7/PR/PHx3P/wNwvZldC/za3f+w/4otsncKAJHytBd9N+BSd3+0eAUzOwX4vrv/W+nGZnYscBrwfTN7zN2/NailFSmDrgGI9GwXMGYvyx4F/j4emhozOzQetfFR4Evx+wows6lmNtnMDgJS7n4PcD1wzOAXX6R36gGI9MDd3zKzp8xsLdBBNOpqwe3ADGBVPET1duCT7v6Ymc0G/hi/pOZt4PPA+4HrzCwPZIC/3381Edk7jQYqIhIonQISEQmUAkBEJFAKABGRQCkAREQCpQAQEQmUAkBEJFAKABGRQP1/VYxo5o1sFvEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot test curve\n",
    "RF_results.plot(x='trees', y='mean_te')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "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>trees</th>\n",
       "      <th>mean_tr</th>\n",
       "      <th>mean_te</th>\n",
       "      <th>std_te</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>0.999448</td>\n",
       "      <td>0.993840</td>\n",
       "      <td>0.000116</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>0.999283</td>\n",
       "      <td>0.993821</td>\n",
       "      <td>0.000148</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>11</td>\n",
       "      <td>0.999560</td>\n",
       "      <td>0.993818</td>\n",
       "      <td>0.000135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>13</td>\n",
       "      <td>0.999643</td>\n",
       "      <td>0.993813</td>\n",
       "      <td>0.000144</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19</td>\n",
       "      <td>0.999797</td>\n",
       "      <td>0.993801</td>\n",
       "      <td>0.000141</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>17</td>\n",
       "      <td>0.999757</td>\n",
       "      <td>0.993790</td>\n",
       "      <td>0.000121</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>15</td>\n",
       "      <td>0.999705</td>\n",
       "      <td>0.993790</td>\n",
       "      <td>0.000124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>0.999023</td>\n",
       "      <td>0.993704</td>\n",
       "      <td>0.000139</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>50</td>\n",
       "      <td>0.999966</td>\n",
       "      <td>0.993704</td>\n",
       "      <td>0.000124</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>40</td>\n",
       "      <td>0.999937</td>\n",
       "      <td>0.993694</td>\n",
       "      <td>0.000120</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   trees   mean_tr   mean_te    std_te\n",
       "4      9  0.999448  0.993840  0.000116\n",
       "3      7  0.999283  0.993821  0.000148\n",
       "5     11  0.999560  0.993818  0.000135\n",
       "6     13  0.999643  0.993813  0.000144\n",
       "9     19  0.999797  0.993801  0.000141\n",
       "8     17  0.999757  0.993790  0.000121\n",
       "7     15  0.999705  0.993790  0.000124\n",
       "2      5  0.999023  0.993704  0.000139\n",
       "13    50  0.999966  0.993704  0.000124\n",
       "12    40  0.999937  0.993694  0.000120"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Ten best combinations according to the mean test score\n",
    "RF_results.sort_values(by='mean_te', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 99.40%\n"
     ]
    }
   ],
   "source": [
    "# Report test score\n",
    "acc_RF = 100*grid_RF_cv.score(X_te, y_te)\n",
    "print('Test accuracy: {:.2f}%'.format(acc_RF))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = grid_RF_cv.predict(X_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn, fp, fn, tp = confusion_matrix(y_te, y_pred).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True positive:  1771\n",
      "True negative:  736967\n",
      "False positive:  194\n",
      "False negative:  4297\n",
      "Positive pixels:  6068\n",
      "Negative pixels:  737161\n"
     ]
    }
   ],
   "source": [
    "print('True positive: ', tp)\n",
    "print('True negative: ', tn)\n",
    "print('False positive: ', fp)\n",
    "print('False negative: ', fn)\n",
    "print('Positive pixels: ', y_te[y_te == 1].size)\n",
    "print('Negative pixels: ', y_te[y_te == 0].size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEWCAYAAABLzQ1kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsKElEQVR4nO3de7xVVb338c+XO8odFLkpXjBD8x6XfDQUD2A3zEfPoTyJHTqmmXq6Hu2pID101HPUUtOyxHsJWSplioSZWd7wkhcMIU1AUITNVQVh83v+mGPBYrP22mvD3uzJ3t/36zVfa80xxxhzrI3+1lhjjjmmIgIzM8ufVk3dADMzK80B2swspxygzcxyygHazCynHKDNzHLKAdrMLKccoG2HSOoo6TeSVkn65Q7Uc7qkBxuybU1F0rGS5jZ1O2zXJ8+DbhkkfRb4KnAQsAZ4DpgcEY/uYL2fA84DPhIRG3e0nXknKYBBETG/qdtizZ970C2ApK8CPwC+D/QG9gauA8Y2QPX7AK+0hOBcCUltmroN1oxEhLdmvAFdgbXAaWXytCcL4IvT9gOgfTo2AlgEfA1YCiwBPp+OfQ94H9iQzjEBmATcXlT3QCCANmn/TOBVsl78a8DpRemPFpX7CPAUsCq9fqTo2MPAJcCfUz0PAr1q+WyF9n+zqP0nAx8DXgGqgG8V5R8CPAasTHmvBdqlY4+kz/JO+rz/UlT/fwJvArcV0lKZ/dM5jkz7fYFlwIim/m/DW/4396Cbv+FAB+DuMnn+HzAMOBw4jCxIfbvo+F5kgb4fWRD+kaTuETGRrFc+NSI6RcSN5RoiaXfgauCkiOhMFoSfK5GvB3BfytsTuBK4T1LPomyfBT4P7Am0A75e5tR7kf0N+gHfBX4K/CtwFHAs8F1J+6W81cBXgF5kf7uRwJcAIuK4lOew9HmnFtXfg+zXxFnFJ46Iv5MF7zsk7QbcBNwcEQ+Xaa8Z4CGOlqAnsCzKD0GcDlwcEUsj4m2ynvHnio5vSMc3RMTvyHqPH9jO9mwCDpHUMSKWRMRLJfJ8HJgXEbdFxMaI+AXwN+CTRXluiohXIuI9YBrZl0ttNpCNt28A7iQLvj+MiDXp/C8BhwJExNMR8Xg67z+AnwAfreAzTYyI9ak9W4mInwLzgCeAPmRfiGZ1coBu/pYDveoYG+0LvF60/3pK21xHjQD/LtCpvg2JiHfIhgXOBpZIuk/SQRW0p9CmfkX7b9ajPcsjojq9LwTQt4qOv1coL+lASb+V9Kak1WS/EHqVqRvg7YhYV0eenwKHANdExPo68poBDtAtwWPAOrJx19osJvt5XrB3Stse7wC7Fe3vVXwwImZExD+R9ST/Rha46mpPoU1vbGeb6uN6snYNioguwLcA1VGm7FQoSZ3IxvVvBCalIRyzOjlAN3MRsYps3PVHkk6WtJuktpJOknR5yvYL4NuS9pDUK+W/fTtP+RxwnKS9JXUFLiockNRb0qfSWPR6sqGS6hJ1/A44UNJnJbWR9C/AYOC329mm+ugMrAbWpt79OTWOvwXst02p8n4IPB0RXyAbW//xDrfSWgQH6BYgIq4kmwP9beBtYCHwZeCelOW/gNnA88ALwDMpbXvONROYmup6mq2Daiuy2SCLyWY2fJR0Aa5GHcuBT6S8y8lmYHwiIpZtT5vq6etkFyDXkPXup9Y4Pgm4RdJKSf9cV2WSxgJjyIZ1IPt3OFLS6Q3WYmu2fKOKmVlOuQdtZpZTDtBmZjnlAG1mllMO0GZmOdWiF3bp1aN1DBzQtqmbYfXwyvO71Z3JcmUNK5ZFxB7bW3708bvH8qpSszG39fTz62dExJjtPVfetOgAPXBAW56cMaCpm2H1MLrv4U3dBKun38ddNe8KrZflVdU8OWPvivK27jOvrrs+dyktOkCbWf4FsIlNTd2MJuEAbWa5FgQborIhjubGAdrMcs89aDOzHAqC6hZ6x7MDtJnl3qbyCwY2Ww7QZpZrAVS30ADtG1XMLPc2ERVtdZH0AUnPFW2rJf2HpB6SZkqal167F5W5SNJ8SXMljS5KP0rSC+nY1ZKU0ttLmprSn5A0sKjM+HSOeZLG19VeB2gzy7UANkRUtNVZV8TciDg8Ig4neyblu2TP67wQmBURg4BZaR9Jg4FxwMFky8ZeJ6l1qu56smdQDkpb4QaZCcCKiDgAuAq4LNXVA5gIDCV77ufE4i+CUhygzSzXgqC6wq2eRgJ/j4jXgbHALSn9FrY8gWgscGd63uRrwHxgiKQ+QJeIeCyyNZtvrVGmUNddwMjUux4NzIyIqohYAcxkS1AvyWPQZpZvAdWVx95ekmYX7d8QETfUkncc2dOEAHpHxBKAiFgiac+U3g94vKjMopS2Ib2vmV4oszDVtVHSKrKHN29OL1GmJAdoM8u17E7Cii2LiKPryiSpHfApih7JVlvWWppUW/r2linJQxxmlnOiusKtHk4CnomIwtPd30rDFqTXpSl9EVC8YE9/ske2LUrva6ZvVUZSG6Ar2SPeaqurVg7QZpZr2UVCVbTVw2fYMrwBMB0ozKoYD9xblD4uzczYl+xi4JNpOGSNpGFpfPmMGmUKdZ0KPJTGqWcAoyR1TxcHR6W0WnmIw8xyLZsHXa/gW5ak3YB/Ar5YlHwpME3SBGABcBpARLwkaRowB9gInBuxeWGQc4CbgY7A/WkDuBG4TdJ8sp7zuFRXlaRLgKdSvosjoqpcWx2gzSz3NtWvd1xWRLxLdtGuOG052ayOUvknA5NLpM8GDimRvo4U4EscmwJMqbStDtBmlmsN3YPelThAm1muBaK6hV4uc4A2s9xryCGOXYkDtJnlWiDej9Z1Z2yGHKDNLNeyG1U8xGFmlku+SGhmlkMRojrcgzYzy6VN7kGbmeVPdpGwZYaqlvmpzWyX4YuEZmY5Vu150GZm+eM7Cc3McmyTZ3GYmeVPtliSA7SZWe4EYoNv9TYzy58IfKOKmVk+yTeqmJnlUeAetJlZbvkioZlZDgXygv1mZnkUwIYWuhZHy/zdYGa7EFFd4VZRbVI3SXdJ+puklyUNl9RD0kxJ89Jr96L8F0maL2mupNFF6UdJeiEdu1qSUnp7SVNT+hOSBhaVGZ/OMU/S+Lra6gBtZrkWZHcSVrJV6IfAAxFxEHAY8DJwITArIgYBs9I+kgYD44CDgTHAdZIKk7KvB84CBqVtTEqfAKyIiAOAq4DLUl09gInAUGAIMLH4i6AUB2gzy72G6kFL6gIcB9wIEBHvR8RKYCxwS8p2C3Byej8WuDMi1kfEa8B8YIikPkCXiHgsIgK4tUaZQl13ASNT73o0MDMiqiJiBTCTLUG9JAdoM8u1CNWnB91L0uyi7awa1e0HvA3cJOlZST+TtDvQOyKWZOeLJcCeKX8/YGFR+UUprV96XzN9qzIRsRFYBfQsU1etWubIu5ntMrKLhBXf6r0sIo4uc7wNcCRwXkQ8IemHpOGMWpTqlkeZ9O0tU5J70GaWc9kzCSvZKrAIWBQRT6T9u8gC9ltp2IL0urQo/4Ci8v2BxSm9f4n0rcpIagN0BarK1FUrB2gzy7XsIqEq2uqsK+JNYKGkD6SkkcAcYDpQmFUxHrg3vZ8OjEszM/Yluxj4ZBoGWSNpWBpfPqNGmUJdpwIPpXHqGcAoSd3TxcFRKa1WHuIws9xr4DsJzwPukNQOeBX4PFlndZqkCcAC4DSAiHhJ0jSyIL4RODciqlM95wA3Ax2B+9MG2QXI2yTNJ+s5j0t1VUm6BHgq5bs4IqrKNdQB2sxyraHvJIyI54BS49Qja8k/GZhcIn02cEiJ9HWkAF/i2BRgSqVtdYA2s9zzQ2PNzHIoAjZscoA2M8udbIjDAdrMLJcqXWejuXGAzpGF89vz/bMHbt5/c0E7PveNN1mzojWPzeiKBN16beDrP1hAz702AvDqnA5c/Z8DeGdNK1q1gmt+9wrtOgQP39uNO6/uTXU1DB25mi98Z8nmev84vRu3X7EXKNhv8Douuu51nvtzJ34ycctNTQv/3p5vXfc6Hzlp1U77/M3ZV69cwNAT17ByWRu+eEI2w2u/we9x3qWL6Lj7Jt5a1I7Lzt2bd9duuSFjj37v89OH53L7Fb2568d71lZ1s1eYZtcSNVqAlhTAlRHxtbT/daBTREyqsPwkYG1E/G89zvk74LMRsVLS2ojoVP+WN50BB6zn+t/PBaC6Gk4/8mCOOWklnbpWM/6bbwJwz896cftVe3HBZYuo3giXn7cP37j6dfY/eB2rq1rTum2wuqo1P7ukL9fOmEu3ntX8zwV78+yfOnHEsWt549V2TL1mT668dx6du1Wzcln2n8Dhx6zdfO7VK1rz+WM+yJEfXd00f4hm6MGpPZh+Uy++8cMtd/r+x/8u5KcX9+WFxzsxatxyTj1nKbf+T5/Nx8+etJinHurcFM3NmZY7xNGYn3o9cIqkXo14jq1ExMfSwie7vOf+1Jk++6ynd/8N7N550+b0de+1Qqkz8fQfO7PvB99j/4PXAdClRzWtW8OSBe3ot996uvXMpmsecewaHv1dNwDuv6MnnzxzGZ27Zce69dq4zbkfva8bHz5+NR12K3sXqtXDi090Ys2KrftD/fdfzwuP7w7As4905v98fMuvleFjVrFkQTtef6XDTm1nXm1KzyWsa2tuGjNAbwRuAL5S84CkfSTNkvR8et27ljoGS3pY0quSzi8qf4+kpyW9VLwYiqR/7MwvhMb08L3dGHHyys37N126F6cfNZiHft2dM76RDVcserUDEnzrM/tx7qgDmfaj7Gdw34Hvs+jv7XlzYTuqN8JfHujK22+03VzmjVfb85VPHcAFnxjEU3/YtodW89zWOF6f24Hho7NfKcd+YhV79N0AQPuO1fzzl5Zy+xW9m7J5uZHN4mhd0dbcNPbvhh8Bp0vqWiP9WuDWiDgUuAO4upbyB5Et0VdYO7VtSv+3iDiKbLL5+ZJ6VtogSWcVVrp6e3l13QWawIb3xeMPduW4T67cnPb5C9/kjqfncMIpK5g+ZQ8AqjfCi0/uzn9e+zpX3DOPvzzQlWf/1InO3ao5778X8f2z9+Frnx5E7wHv07pN1huuroY3XmvP//xqPhdd9zo/+PoA1q7a8h/28rfa8I+XO3L0CA9vNLYrvzqAT565jGsfeIWOnarZ+H7WAzzjG29x90/3YN27zS/gbI/CjSoNcav3rqZRLxJGxGpJtwLnA+8VHRoOnJLe3wZcXksV90XEemC9pKVAb7IFR86X9OmUZwDZ/fHLK2zTDWQ9e44+rEMuf8M/9VBnDvjQu3TfY9vhh+M/vYLvfG4/zvjGm+zRZwOHDn+Hrmko48MnrGb+Cx054ti1DBu1mmGjsiD7u9t70rpV9lF79dnAB498lzZtYa+936f//ut547V2fODw7J/nkd904yMnraRN221ObQ1s4fwOfOsz+wPQb7/1DB2Z/XsddMS7/J+Pr2TCtxfTqUs1sUm8v74V029qFj8Ot0tzHL6oxM4Yef8B2RMGdi+Tp7ZAub7ofTXQRtII4ERgeEQcBjwLNKuBuofv6b7VEMMbr7bb/P7xGV0ZcED2ZzlqxBpem9OBde+K6o3w/GOd2PvA7Fjh4t+ala35zc29GPPZ7Jb/j4xZxV//kl07XbW8NYv+3p4+e79f67mt8XTtmQ1pSMFnL3iL396W/RD82qcPYPzQwYwfOpi7f7YHd16zZ4sOzg25WNKuptGn2aUFQqaRBenCPeh/IVtA5DbgdODRelTZlexxMu9KOggY1pDtbWrr3hXP/KkzF1y+5Wr/jd/vy6K/t6dVK9iz3/ucf1m2TnjnbtWc8sW3Oe9jByLBkBNWM/TErBd2/Xf68eqcjgCc/pU36b9/FriPHrGGZ/7YmX//6EG0ah38+3cW06VH1gN/c2E73l7clkOHr92ZH7lFuPC61zl0+Fq69tjI7bPncNsVvem42yY+eeYyAP58f1cevLNHE7cyv1rqLA5lq+A1QsVF09wk9QZeAy6PiEnpIYpTgF5kTzf4fEQsqFF+EkXT7CS9CHwCWALcQ/YkgrnAHsCkiHhY0j+AoyNiWSXT7I4+rEM8OWNAuSyWM6P7Ht7UTbB6+n3c9XQdi+iX1f2gPeOEKadWlPfXx1y/Q+fKm0brQRcHx4h4C9itaP8fwAl1lJ9UY7941aiTaikzsNT5zWzX1hyHLyrhOwnNLNd8J6GZWY45QJuZ5VBDL9i/K3GANrPca6nzoB2gzSzXImCjF+w3M8unljrE0TK/lsxsl9HQa3GkRdVekPScpNkprYekmZLmpdfuRfkvkjRf0lxJo4vSj0r1zJd0tZStMympvaSpKf2JdN9Hocz4dI55ksbX1VYHaDPLvQhVtNXD8RFxeNFNLRcCsyJiEDAr7SNpMNldzwcDY4DrJBVWsboeOItsLaBB6Thkd02viIgDgKuAy1JdPYCJwFC2LAC3+YugFAdoM8u9nbAe9FjglvT+FuDkovQ7I2J9RLwGzAeGSOoDdImIxyK7HfvWGmUKdd0FjEy969HAzIioiogVwEy2BPWSHKDNLNci6rVYUq/CcsJpO6tUlcCDaU35wvHeEbEkO18sAQrPGOsHLCwquyil9Uvva6ZvVSYiNgKrgJ5l6qqVLxKaWc6J6spncSyrYC2OYyJisaQ9gZmS/lb25NuKMunbW6Yk96DNLPcacgw6Ihan16XA3WTjwW+lYQvS69KUfRHZmvMF/YHFKb1/ifStykhqQ7YCZ1WZumrlAG1mudaQ60FL2l1S58J7YBTwIjAdKMyqGA/cm95PB8almRn7kl0MfDINg6yRNCyNL59Ro0yhrlOBh9I49QxglKTu6eLgqJRWKw9xmFm+RTYO3UB6A3enGXFtgJ9HxAOSngKmSZoALABOA4iIl9J69nPInrN6bkQUnpV3DnAz0BG4P20ANwK3SZpP1nMel+qqknQJ8FTKd3FEVJVrrAO0meVeQ93qHRGvAoeVSF8OjKylzGRgcon02cAhJdLXkQJ8iWNT2PLgkjo5QJtZrkX9LhI2Kw7QZpZ7jfTgp9xzgDaz3KvnXYLNhgO0meVahAO0mVlutdTV7BygzSz3PAZtZpZDgdjkWRxmZvnUQjvQDtBmlnO+SGhmlmMttAvtAG1muecedA2SrqHM91ZEnN8oLTIzKxLApk0O0DXN3mmtMDOrTQDuQW8tIm4p3pe0e0S80/hNMjPbWkudB13n5EJJwyXNAV5O+4dJuq7RW2ZmVhAVbs1MJbO/f0D2NNrlABHxV+C4RmyTmVmRyh531RwvJFY0iyMiFqYnEBRU15bXzKzBNcPecSUqCdALJX0ECEntgPNJwx1mZo0uIFroLI5KhjjOBs4F+gFvAIenfTOznUQVbs1LnT3oiFgGnL4T2mJmVloLHeKoZBbHfpJ+I+ltSUsl3Stpv53RODMzwLM4yvg5MA3oA/QFfgn8ojEbZWa2WeFGlUq2CklqLelZSb9N+z0kzZQ0L712L8p7kaT5kuZKGl2UfpSkF9Kxq5VmUkhqL2lqSn9C0sCiMuPTOeZJGl9XOysJ0IqI2yJiY9pup1l+V5lZXmWPvap7q4cL2Hqyw4XArIgYBMxK+0gaDIwDDgbGANdJap3KXA+cBQxK25iUPgFYEREHAFcBl6W6egATgaHAEGBi8RdBKbUG6PSN0gP4g6QLJQ2UtI+kbwL3VfY3MDNrAJtU2VYBSf2BjwM/K0oeCxTunr4FOLko/c6IWB8RrwHzgSGS+gBdIuKxiAjg1hplCnXdBYxMvevRwMyIqIqIFcBMtgT1kspdJHyarKdc+NRfLDoWwCXlKjYzayiqvHfcS1LxOkI3RMQNNfL8APgm0LkorXdELAGIiCWS9kzp/YDHi/ItSmkb0vua6YUyC1NdGyWtAnoWp5coU1K5tTj2LVfQzGynqN8FwGURcXRtByV9AlgaEU9LGlFBfaW65VEmfXvLlFTRnYSSDgEGAx021xpxayVlzcx2TP0uANbhGOBTkj5GFs+6SLodeEtSn9R77gMsTfkXAQOKyvcHFqf0/iXSi8ssktQG6ApUpfQRNco8XK6xlUyzmwhck7bjgcuBT9VVzsyswTTQNLuIuCgi+kfEQLKLfw9FxL8C04HCrIrxwL3p/XRgXJqZsS/ZxcAn03DIGknD0vjyGTXKFOo6NZ0jgBnAKEnd08XBUSmtVpX0oE8FDgOejYjPS+rN1oPrZmaNa1Ojn+FSYJqkCcAC4DSAiHhJ0jRgDrARODciCmsRnQPcDHQE7k8bwI3AbZLmk/Wcx6W6qiRdAjyV8l0cEVXlGlVJgH4vIjZJ2iipC1nX3zeqmNnO0UgL9kfEw6QhhohYDoysJd9kYHKJ9NnAISXS15ECfIljU4AplbaxkgA9W1I34KdkMzvWAk9WegIzsx1Vj1kczUola3F8Kb39saQHyOb+Pd+4zTIzK+IAvTVJR5Y7FhHPNE6TzMwMyvegryhzLIATGrgtO90rz+/G6H5HNHUzrF5aaFeqhfMQRw0RcfzObIiZWUlBxbdxNzcV3ahiZtak3IM2M8snD3GYmeVVCw3QldzqLUn/Kum7aX9vSUMav2lmZomfqFKr64DhwGfS/hrgR43WIjOzIorKt+amkiGOoRFxpKRnASJihaR2jdwuM7MtPIujVhvSI14CQNIe7IylS8zMkubYO65EJUMcVwN3A3tKmgw8Cny/UVtlZlashY5BV7IWxx2SniZb6UnAyRHxch3FzMwaRjMdX65EnQFa0t7Au8BvitMiYkFjNszMbDMH6Frdx5bnaXUA9gXmkj2G3Mys0amFXvWqZIjjQ8X7aZW7L9aS3czMGki97ySMiGckfbgxGmNmVpKHOEqT9NWi3VbAkcDbjdYiM7NivkhYVuei9xvJxqR/1TjNMTMroYUG6LLzoNMNKp0i4ntpmxwRd6SHIpqZ7RwNNA9aUgdJT0r6q6SXJH0vpfeQNFPSvPTavajMRZLmS5oraXRR+lGSXkjHrpaklN5e0tSU/oSkgUVlxqdzzJM0vq721hqgJbVJjxev9dFXZmaNTWSzOCrZKrAeOCEiDgMOB8ZIGgZcCMyKiEHArLSPpMHAOLJZa2OA61LHFeB64CxgUNrGpPQJwIqIOAC4Crgs1dUDmAgMBYYAE4u/CEop14MuPLn7OUnTJX1O0imFraI/hZnZjmrAxZIiszbttk1bAGOBW1L6LcDJ6f1Y4M6IWB8RrwHzgSGS+pA9QPuxiAjg1hplCnXdBYxMvevRwMyIqIqIFcBMtgT1kioZg+4BLCd7BmFhPnQAv66grJnZjqt8DLqXpNlF+zdExA3FGVIP+GngAOBHEfGEpN4RsQQgIpZI2jNl7wc8XlR8UUrbkN7XTC+UWZjq2ihpFdCzOL1EmZLKBeg90wyOF9kSmAta6JC9mTWJyiPOsog4umxV2dDt4ZK6AXdLOqRM9lLL6NWMhzVbuT1lSio3xNEa6JS2zkXvC5uZ2U7RGOtBR8RK4GGyYYa30rAF6XVpyrYIGFBUrD+wOKX3L5G+VRlJbYCuQFWZumpVrge9JCIuLlfYzGynaKDf7Gm55A0RsVJSR+BEsot404HxwKXp9d5UZDrwc0lXAn3JLgY+GRHVktakC4xPAGcA1xSVGQ88BpwKPBQRIWkG8P2iC4OjgIvKtbdcgG6ZK2SbWb5Eg67F0Qe4JY1DtwKmRcRvJT0GTJM0AVgAnAYQES9JmgbMIbsP5Nw0RAJwDnAz0BG4P20ANwK3SZpP1nMel+qqknQJ8FTKd3FEVJVrbLkAPbLyz2xm1ogaqAcdEc8DR5RIX04tMS8iJgOTS6TPBrYZv073iZxWS11TgCmVtrfWAF1XZDcz21l8q7eZWV45QJuZ5VAzfZxVJRygzSzXhIc4zMxyywHazCyvHKDNzHLKAdrMLIf8RBUzsxxzgDYzy6cGvNV7l+IAbWa55yEOM7M88o0qZmY55gBtZpY/vpPQzCzHtKllRmgHaDPLN49Bm5nll4c4zMzyygHazCyf3IM2M8srB2gzsxxq2Kd671JaNXUDzMzKKcyDrmSrsy5pgKQ/SHpZ0kuSLkjpPSTNlDQvvXYvKnORpPmS5koaXZR+lKQX0rGrJSmlt5c0NaU/IWlgUZnx6RzzJI2vq70O0GaWfxGVbXXbCHwtIj4IDAPOlTQYuBCYFRGDgFlpn3RsHHAwMAa4TlLrVNf1wFnAoLSNSekTgBURcQBwFXBZqqsHMBEYCgwBJhZ/EZTiAG1muddQPeiIWBIRz6T3a4CXgX7AWOCWlO0W4OT0fixwZ0Ssj4jXgPnAEEl9gC4R8VhEBHBrjTKFuu4CRqbe9WhgZkRURcQKYCZbgnpJHoPehbRqFVxz/yssf7Mt3x2/H1/49hsM+6fVbHhfLHm9PVd8dQDvrG5Dm7abuOCyRQw69F0i4Prv9uP5xzrTcfdqrrh73ub6evXZwEO/7s6PJ/Zvwk/VfH31ygUMPXENK5e14YsnfACAb/34H/Tffz0Au3ep5p3VrfnSP32A4z+9gtO+tHRz2X0/uI5zRx/Iqy915Mz/XMKJp62gU9dqTh70oSb5LE2qfjeq9JI0u2j/hoi4oVTGNPRwBPAE0DsilkAWxCXtmbL1Ax4vKrYopW1I72umF8osTHVtlLQK6FmcXqJMSbkN0JLWRkSnov0zgaMj4stlynwKGBwRl0qaBKyNiP9t9MbuJCd/4W0WzmvPbp2zKybPPNKZKf/dl03VYsK3FjPuy0u58ft9OemzywE4+8SD6NpzA5Nvf5XzPnYg773Tmi+NOmhzfdfeP5dHf9etKT5Ki/Dg1B5Mv6kX3/jhlv8nv3/2wM3vz/ruYt5Zk/2I/cPd3fnD3dmv3YEHvcekm/7Bqy91BODxmV2YflMvpvz5bzuv8TlTj4uEyyLi6DrrkzoBvwL+IyJWp+HjkllLpEWZ9O0tU1KzGuKIiOkRcWlTt6Mx9OrzPkNGrub+X/TcnPbMI13YVJ39m7/8zG706rMBgL0PXM+zj2bfbauWt2Xt6tYceNi7W9XXd9/1dOu1kRef2H0nfYKW58UnOrFmRW19oOC4T63kD/dsOwR5/Mkrefiebpv3//bM7lQtbds4jdxFaFNlW0V1SW3JgvMdEfHrlPxWGrYgvRZ+ziwCBhQV7w8sTun9S6RvVUZSG6ArUFWmrlrtkgFa0ifT1dFnJf1eUu+Ufqaka5u6fY3h7O+9wc/+qy9Ry3+Eo8dV8dQfOgPw6pwODB+9ilatg94D1jPoQ++yR98NW+U/fuwK/ji9G6W/1K2xHTL0HVa83YbFr7Xf5lgWuLvt/EblVdBgFwnTWPCNwMsRcWXRoelAYVbFeODeovRxaWbGvmQXA59MwyFrJA1LdZ5Ro0yhrlOBh9I49QxglKTu6eLgqJRWq9wOcQAdJT1XtN+D7IMDPAoMi4iQ9AXgm8DXKqlU0llkV17pwG4N19pGNPTEVaxc1ob5L+zGocPXbHP8M+e/SfVG8dCvs97YjDt7sveg9Vx7/1yWLmrHnNm7U71x60D80bEruPz8fXZK+21bNXvJBR844h3Wv9eK1+d23PmNyrEGvJPwGOBzwAtF8eVbwKXANEkTgAXAaQAR8ZKkacAcshkg50ZEdSp3DnAz0BG4P22QfQHcJmk+Wc95XKqrStIlwFMp38URUVWusXkO0O9FxOGFncIYdNrtD0xNP0XaAa9VWmm6YHADQBf12CXuTxp89DsMG7WaD5/wEu3aB7t1ruabV7/O5efvw4mnVTHkxNVc+M8HUOgNb6oWP5m05drDVfe+whtFPbX9Br9H6zYw/4Vd4wuquWnVOjjmY6v48phB2xwbMbZ04G7xGuj/1Ih4lNp/No6spcxkYHKJ9NnAISXS15ECfIljU4AplbY3zwG6nGuAKyNiuqQRwKQmbU0ju+nSvtx0aV8ADh2+hlPPfpvLz9+Ho0es5p+/9Bbf+L+DWL9uy2hV+w6bQMH691pz5LFrqN4oFszrsPn4iLErHASa0JHHrmHh/PYsW9Juq3QpOPYTq/j6Kfs3UcvyyQv273q6Am+k93XejdNcnftfi2jbPvjvO+cD2cWkqy8cQLdeG5j881eJTbD8zbbbDGUc98mVfOdz+zVFk1uUC697nUOHr6Vrj43cPnsOt13Rmxm/6MlHa+klf2jYOyxb0pY3F2w9Lj3h24s5/uSVtO+4idtnz+GBX/Tg9iv22kmfIgciWuyC/YrK7r7Z6cpNs5M0luwOnTfI5ih+OCJG1MgziTqm2XVRjxja6sTG/BjW0HL636vV7vdx19OVTH2rTedu/eOI4y6oKO+ffvPNHTpX3uS2B10cnNP+zWQD8kTEvWy5YlpbnkmN20Iz21k8xGFmlkcBtNAhDgdoM8u/lhmfHaDNLP88xGFmllMtdRaHA7SZ5Vv9VrNrVhygzSzXshtVWmaEdoA2s/xroc8kdIA2s9xzD9rMLI88Bm1mllctdy0OB2gzyz8PcZiZ5VDU65mEzYoDtJnln3vQZmY51TLjswO0meWfNrXMMQ4HaDPLt6DF3qjSqu4sZmZNRwSKyrY665KmSFoq6cWitB6SZkqal167Fx27SNJ8SXMljS5KP0rSC+nY1ZKU0ttLmprSn5A0sKjM+HSOeZIqelSfA7SZ5V9EZVvdbgbG1Ei7EJgVEYOAWWkfSYOBccDBqcx1klqnMtcDZwGD0laocwKwIiIOIHss32Wprh7ARGAoMASYWPxFUBsHaDPLvwYK0BHxCFBVI3kscEt6fwtwclH6nRGxPiJeA+YDQyT1AbpExGORPdT11hplCnXdBYxMvevRwMyIqIqIFcBMtv2i2IbHoM0s3+o3Bt1L0uyi/Rsi4oY6yvSOiCUAEbFE0p4pvR/ZQ6kLFqW0Del9zfRCmYWpro2SVgE9i9NLlKmVA7SZ5V49ZnEsa8CneqtEWpRJ394ytfIQh5nlXIXDG9t/M8tbadiC9Lo0pS8CBhTl6w8sTun9S6RvVUZSG6Ar2ZBKbXWV5QBtZvkWNHaAng4UZlWMB+4tSh+XZmbsS3Yx8Mk0HLJG0rA0vnxGjTKFuk4FHkrj1DOAUZK6p4uDo1JaWR7iMLP8a6B50JJ+AYwgG6teRDaz4lJgmqQJwALgNICIeEnSNGAOsBE4NyKqU1XnkM0I6QjcnzaAG4HbJM0n6zmPS3VVSboEeCrluzgial6s3IYDtJnlXkMt2B8Rn6nl0Mha8k8GJpdInw0cUiJ9HSnAlzg2BZhScWNxgDazXYEXSzIzy6EIqG6Z93o7QJtZ/rkHbWaWUw7QZmY5FICfSWhmlkcB4TFoM7P8CXyR0MwstzwGbWaWUw7QZmZ5tEPrbOzSHKDNLN8C8ENjzcxyyj1oM7M88q3eZmb5FBCeB21mllO+k9DMLKc8Bm1mlkMRnsVhZpZb7kGbmeVRENXVdWdrhhygzSzfvNyomVmOeZqdmVn+BBDuQZuZ5VB4wX4zs9xqqRcJFS10+gqApLeB15u6HY2gF7CsqRth9dKc/832iYg9trewpAfI/j6VWBYRY7b3XHnTogN0cyVpdkQc3dTtsMr538xKadXUDTAzs9IcoM3McsoBunm6oakbYPXmfzPbhsegzcxyyj1oM7OccoA2M8spB+gckhSSrija/7qkSfUoP0nS1+t5zt9J6pber61PWSuv5t9T0pmSrq2jzKckXZje1/vf05oHB+h8Wg+cIqnSyfk7LCI+FhErd9b5rLyImB4RlzZ1O6xpOUDn00ayq/pfqXlA0j6SZkl6Pr3uXUsdgyU9LOlVSecXlb9H0tOSXpJ0VlH6P3bmF4JlJH1S0hOSnpX0e0m9U3qdvWxr/hyg8+tHwOmSutZIvxa4NSIOBe4Arq6l/EHAaGAIMFFS25T+bxFxFHA0cL6kng3fdKuho6TnChtwcdGxR4FhEXEEcCfwzaZooOWTF0vKqYhYLelW4HzgvaJDw4FT0vvbgMtrqeK+iFgPrJe0FOgNLCILyp9OeQYAg4DlDd1+28p7EXF4YUfSmWRfkAD9gamS+gDtgNd2eusst9yDzrcfABOA3cvkqW0i+/qi99VAG0kjgBOB4RFxGPAs0GGHW2k74hrg2oj4EPBF/O9hRRygcywiqoBpZEG64C/AuPT+dLKfyJXqCqyIiHclHQQMa5CG2o7oCryR3o9vyoZY/jhA598VbL3U4vnA5yU9D3wOuKAedT1A1pN+HrgEeLzBWmnbaxLwS0l/ovkuN2rbybd6m5nllHvQZmY55QBtZpZTDtBmZjnlAG1mllMO0GZmOeUAbWVJqk63KL8o6ZeSdtuBum6WdGp6/zNJg8vkHSHpI9txjpJrilSy1kh9V/HzKnPW2BygrS7vRcThEXEI8D5wdvFBSa23p9KI+EJEzCmTZQRQ7wBt1pw4QFt9/Ak4IPVu/yDp58ALklpL+h9JT6VV9r4IoMy1kuZIug/Ys1BRWmnv6PR+jKRnJP01rdA3kOyL4Cup936spD0k/Sqd4ylJx6SyPSU9mFaD+wmguj5EbSv6pWNXpLbMkrRHSttf0gOpzJ/SXZhmjc6LJVlFJLUBTiK7GxGyVfIOiYjXUpBbFREfltQe+LOkB4EjgA8AHyJbrGkOMKVGvXsAPwWOS3X1iIgqST8G1kbE/6Z8PweuiohH0xKrM4APAhOBRyPiYkkfB7YKuLX4t3SOjsBTkn4VEcvJ1jx5JiK+Jum7qe4vky39enZEzJM0FLgOOGE7/oxm9eIAbXXpmJbIhKwHfSPZ0MOTEVFYeW0UcGhhfJlsfYlBwHHALyKiGlgs6aES9Q8DHinUldYfKeVEsjWuC/tdJHVO5zgllb1P0ooKPlNtK/ptAqam9NuBX0vqlD7vL4vO3b6Cc5jtMAdoq8tWS2UCpED1TnEScF5EzKiR72PUvtpecdlK1htoRbYKX/HSq4W2VLxeQY0V/d6V9DC1ryAX6bwra/4NzHYGj0FbQ5gBnFN4KICkAyXtDjwCjEtj1H2A40uUfQz4qKR9U9keKX0N0Lko34Nkww2kfIent4+QreqHpJOA7nW0tdyKfq2Awq+Az5INnawGXpN0WjqHJB1WxznMGoQDtDWEn5GNLz8j6UXgJ2S/zu4G5gEvANcDf6xZMCLeJhs3/rWkv7JliOE3wKcLFwnJVvE7Ol2EnMOW2STfA46T9AzZUMuCOtpabkW/d4CDJT1NNsZcePLJ6cCE1L6XgLEV/E3MdphXszMzyyn3oM3McsoB2swspxygzcxyygHazCynHKDNzHLKAdrMLKccoM3Mcur/A7ENDRgLmVt0AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "disp = plot_confusion_matrix(grid_RF_cv, X_te, y_te, labels=[0, 1], display_labels=['No hail', 'Hail'], normalize=None)\n",
    "disp.ax_.set_title('Confusion matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of Detection (POD): 29.185893210283453\n",
      "False Alarm Rate (FAR): 0.026317181728279168\n",
      "Positive Predictive Value (PPV): 0.9012722646310433\n"
     ]
    }
   ],
   "source": [
    "pod = 100*tp/(tp+fn)\n",
    "far = 100*fp/(fp+tn)\n",
    "ppv = tp/(tp+fp)\n",
    "print('Probability of Detection (POD):', pod)\n",
    "print('False Alarm Rate (FAR):', far)\n",
    "print('Positive Predictive Value (PPV):', ppv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Basically the model has a hard time detecting hail. Hence the low POD and FAR."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(model_name, 'wb') as bfile: \n",
    "    s = pickle.dump(grid_RF_cv.best_estimator_, bfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(model_name, 'rb') as bfile:\n",
    "    rf_model_best = pickle.load(bfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_importance = rf_model_best['RF'].feature_importances_\n",
    "std = np.std([tree.feature_importances_ for tree in rf_model_best['RF'].estimators_], axis=0)\n",
    "inds = np.argsort(f_importance)[::-1]\n",
    "f_importance_ordered = f_importance[inds]\n",
    "features_ordered = np.array(features)[inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAEICAYAAAAwft9dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhb0lEQVR4nO3deZhcVb3u8e8bxAQQwhCVpiO0YkBRJGAUxKPGB44MKoErynAVoghOOCAgUTyeiFMiKBHxcozoZRDiAQSMghD00A5A0A4knRBkCAmQMASCtmAgYvI7f+zVUFSqe1VXV3VVOu/nefrp2sNa+7d3qvvN2rWqWhGBmZmZ9W1EswswMzNrdQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiaWd1I+pKk85tdh1m9ye+zNGsNkpYBLwfWlqzeJSIeGmSfH42I3wyuug2PpKnAqyPig82uxTZ8HlmatZb3RsRLSr5qDsp6kPSiZh6/Vhtq3da6HJZmLU7SaEk/lvSwpBWSvi5pk7RtZ0n/I2mVpMclXSJp67TtYmBH4JeSnpL0BUkTJS0v63+ZpP3T46mSrpD0U0l/Byb3d/wKtU6V9NP0uENSSPqwpAcl/VXSxyW9SVK3pL9JOrek7WRJN0n6vqQeSX+RtF/J9h0kzZb0hKR7JR1fdtzSuj8OfAk4Ip37grTfhyXdKelJSfdJ+lhJHxMlLZd0sqSV6Xw/XLJ9M0nfkXR/qu+PkjZL2/aRdHM6pwWSJtbwT20tzGFp1vouBP4FvBrYE3gX8NG0TcC3gB2A1wKvAKYCRMSHgAd4frT67SqPNwm4AtgauCRz/GrsDYwDjgBmAKcD+wOvAz4g6R1l+94HjAH+E7hS0rZp2yxgeTrXw4FvloZpWd0/Br4J/Hc69z3SPiuB9wBbAR8Gzpa0V0kf2wOjgXbgOOAHkrZJ284C3gjsC2wLfAFYJ6kduAb4elp/CvBzSS8dwDWyFuewNGstV6fRyd8kXS3p5cBBwOci4h8RsRI4GzgSICLujYgbImJNRDwGfBd4R9/dV+WWiLg6ItZRhEqfx6/S1yLimYiYA/wDmBURKyNiBfAHigDutRKYERHPRsR/A3cB75b0CuDfgNNSX/OB84EPVao7Ip6uVEhEXBMRS6LwO2AO8LaSXZ4FzkjHvxZ4CthV0gjgI8BnI2JFRKyNiJsjYg3wQeDaiLg2HfsGoAs4eADXyFqc7+ubtZZDSyfjSHozsCnwsKTe1SOAB9P2lwHnUPzC3zJt++sga3iw5PFO/R2/So+WPH66wvJLSpZXxAtnHd5PMZLcAXgiIp4s2zahj7orknQQxYh1F4rz2BxYWLLLqoj4V8ny6lTfGGAUsKRCtzsB75f03pJ1mwI35uqxDYfD0qy1PQisAcaU/RLv9S0ggDdExCpJhwLnlmwvn+7+D4qAACC99lh+u7C0Te749dYuSSWBuSMwG3gI2FbSliWBuSOwoqRt+bm+YFnSSODnwDHALyLiWUlXU9zKznkceAbYGVhQtu1B4OKIOH69VjZs+DasWQuLiIcpbhV+R9JWkkakST29t1q3pLhV+Lf02tmpZV08CryqZPluYJSkd0vaFPgyMHIQx6+3lwGfkbSppPdTvA57bUQ8CNwMfEvSKElvoHhN8ZJ++noU6Ei3UAFeTHGujwH/SqPMd1VTVLol/RPgu2mi0SaS3pIC+KfAeyUdkNaPSpOFxg789K1VOSzNWt8xFL/oF1PcYr0CaEvbvgrsBfRQTDK5sqztt4Avp9dAT4mIHuCTFK/3raAYaS6nf/0dv95upZgM9DjwDeDwiFiVth0FdFCMMq8C/jO9PtiXy9P3VZJuSyPSzwCXUZzH0RSj1mqdQnHL9s/AE8B0YEQK8kkUs28foxhpnop/vw4r/lACM2sJkiZTfIDCvzW7FrNy/p+PmZlZhsPSzMwsw7dhzczMMjyyNDMzy/D7LIehMWPGREdHR7PLMDPboMybN+/xiKj4MYUOy2Goo6ODrq6uZpdhZrZBkXR/X9t8G9bMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYY/lGAYWriih44p1zS7DLONxrJp7252CdZgHlmamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhnZsJR0tqTPlSxfL+n8kuXvSApJu5a1myHpC/30+0VJ90q6S9IBJetfLGmmpLsl/UXS+9L6z0taLKlb0m8l7dRHv1MlnZIeXyBpqaT5khZI2i9zriemmkLSmJL1knRO2tYtaa+SbSdJukPSIkmzJI3qp//Jknbor4ZMfV+qta2ZmdWumg9Svxl4PzBD0ghgDLBVyfZ9geuAI4GvAqT9DgfeWqlDSbul/V8H7AD8RtIuEbEWOB1YGRG7pH62Tc1uByZExGpJnwC+DRxRRf2nRsQVkt4JzATG9bPvTcCvgM6y9QelduOAvYHzgL0ltQOfAXaLiKclXZbO64I++p8MLAIeqqLuSr4EfLPGtmZWJ49cOuUFyxPnnrnePp2dnUNUjQ2Fam7D3kQRiFCE2yLgSUnbSBoJvBY4jSIker0dWBYR9/fR5yTgZxGxJiKWAvcCb07bPgJ8CyAi1kXE4+nxjRGxOu0zFxhbzQmWuAVo72+HiLg9Ipb1Ue9FUZgLbC2pLW17EbCZpBcBm9NHEEo6HJgAXJJGuptJeqOk30mal0bsbZJGp9H2rqndLEnHS5qWjjNf0iUV+j9BUpekrrWre6q7ImZmVpXsyDIiHpL0L0k7UoRmb+i8BegBuiOiW9I6SXtExAKK4JzVT7ftFIHXaznQLmnrtPw1SROBJcCJEfFoWfvjgF/nai9zIHD1ANv0agceLFleDrRHRJeks4AHgKeBORExp1IHaXR7InBKarcp8H1gUkQ8JukI4BsR8ZG03wWSvgdsExE/guI2cUSM76P/mRQjZ0a2jYsaz9PMqrD90dNesNzpP9E17FU7wad3dNkblreULN+c9pkFHJlGWJOAy/vpTxXWBUV4jwVuioi90nHOekFD6YMUI7T173tUdqak+4CfUvstzIr1StqG4lxfSXE7eYtUXzV2BV4P3CBpPvBl0mg5Im4AFgI/AD5aY81mZlYn1YblzRTBuDvFbdi5FCPLfSmCFIqw/ACwP8Voc2U//S0HXlGyPJbi9uUqYDVwVVp/OVA6mWZ/itc0D4mINWndN9Ktyfl9HOtU4NUUYXRhFec6kHr3B5ZGxGMR8SxwJc/fss4RcEdEjE9fu0fEu+C513xfSzFa3ba/TszMrPEGMrJ8D/BERKyNiCeArSkC8xaAiFhCEXbT6P8WLMBsilHoSEmvpJg486eICOCXwMS0337AYgBJewI/pAjK54I4Ik7vDZy+DhYR64DvASNKZ94OwGzgmDQrdh+gJyIeprj9uo+kzSUp1XtnP/08CWyZHt8FvFTSW9L5bSrpdWnbSamfo4CfpFu2AM+WPDYzsyFSbVgupJgFO7dsXU/vBJxkFvAanh8ZVhQRdwCXUQThdcCn0kxYKCYLTZXUDXwIODmtPxN4CXB5GknOrrL23mMG8HWgv7ezfEbScoqRY7eef4vMtcB9FBORfgR8MvV5K3AFcBvF9RhBet2wDxcA/5VGwZtQzBieLmkBMB/YV9IuFLdeT46IPwC/pxgVk/rurjTBx8zMGkdFhthwMrJtXLQdO6PZZZhtNJZ5gs+wIGleREyotM2f4GNmZpZRzYcS1Cy9Pji9bPXSiDiskcfNkXQVxQzWUqdFxPV16v8HrP+BDN+LiP9fj/7NzGxoNTQsU/jUJYDqqdFhHRGfamT/ZmY2tHwb1szMLMNhaWZmluGwNDMzy3BYmpmZZTR0go81x+7to+ny+77MzOrGI0szM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8vwbNhhaOGKHjqmXNPsMsws8V8l2fB5ZGlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZWxUYSnpqbLlyZLOTY+nSlohab6kxZKOSuuXStq1rN0MSV8YusrNzKyZ/EHqL3R2RJwlaRwwT9IVwM+AI4GvAkgaARwOvLURBUjaJCLWNqJvMxtaj1w6BYCJc898bl1nZ2eTqrHB2KhGltWKiHuA1cA2wCyKsOz1dmBZRNxfqW0arV4p6TpJ90j6dsm2oyQtlLRI0vSS9U9JOkPSrcBb0vJ0SfMk/UbSmyV1SrpP0iF9HPcESV2Sutau7qnDVTAzs14b28hyM0nzS5a3BWaX7yRpL+CeiFgJrJS0TtIeEbGAIjhnZY4zHtgTWAPcJen7wFpgOvBG4K/AHEmHRsTVwBbAooj4Sjr+FkBnRJwm6Srg68C/A7sBF1aqOSJmAjMBRraNiyquhZk12PZHTwOg03+ia4O3sYXl0xExvndB0mRgQsn2kyQdD7wKOLBk/SzgSEl3AJOAr2SO89uI6EnHWAzsBGxHEYCPpfWXUIxSr6YI0p+XtP8ncF16vBBYExHPSloIdFR5rmZmVie+DftCZ0fErsARwEWSRqX1s4APAPsD3WnE2Z81JY/XUvynRP3s/0zZ65TPRkTv6HBdb38RsY6N7z84ZmZN57CsICKuBLqAY9PyEmAVMI38Ldi+3Aq8Q9IYSZsARwG/q0O5ZmbWYA7Lvp0BfD7NfoUiJF8DXFVLZxHxMPBF4EZgAXBbRPyiHoWamVlj6fm7fTZcjGwbF23Hzmh2GWaWLPMEnw2CpHkRMaHSNo8szczMMjxZpEaSDqB4K0ippRFxWDPqMTOzxnFY1igirgeub3YdZmbWeL4Na2ZmluGwNDMzy3BYmpmZZTgszczMMjzBZxjavX00XX5fl5lZ3XhkaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGZ4NOwwtXNFDx5Rrml2GmZXwXx7ZsHlkaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZlltFRYSnqqbHmypHPT46mSVkiaL2mxpKPS+qWSdi1rN0PSF4au8upIGi/p4EG0P1TSbvWsyczM8ja0D1I/OyLOkjQOmCfpCuBnwJHAVwEkjQAOB97aiAIkbRIRa2tsPh6YAFxbY/tDgV8Bi2tsb2ZD4JFLp6y3buLcM9db19nZOQTVWD201MiyWhFxD7Aa2AaYRRGWvd4OLIuI+yu1TaPVKyVdJ+keSd8u2XaUpIWSFkmaXrL+KUlnSLoVeEtani5pnqTfSHqzpE5J90k6pI/jvhg4AzgijY6PkLSFpJ9I+rOk2yVNSvueI+kr6fEBkn4vaV/gEODM1H7nsv5PkNQlqWvt6p4BX1MzM+tbq40sN5M0v2R5W2B2+U6S9gLuiYiVwEpJ6yTtERELKIJzVuY444E9gTXAXZK+D6wFpgNvBP4KzJF0aERcDWwBLIqI3gDbAuiMiNMkXQV8Hfh3YDfgwko1R8Q/UwBOiIgTUz/fBP4nIj4iaWvgT5J+A0wB/izpD8A5wMERsUTSbOBXEXFFhf5nAjMBRraNi8z5m1kDbX/0tPXWdfpPdG3QWi0sn46I8b0LkiZT3LbsdZKk44FXAQeWrJ8FHCnpDmAS8JXMcX4bET3pGIuBnYDtKALwsbT+EopR6tUUQfrzkvb/BK5LjxcCayLiWUkLgY4qzxXgXcAhkk5Jy6OAHSPiznSevwdOioglA+jTzMzqrNXCMqf3Ncv/A1wkaeeIeIYiLOcAvwO604izP2tKHq+luA7qZ/9nyl6nfDYiekdv63r7i4h1kgZyTQW8LyLuqrBtd2AVsMMA+jMzswbYUF+zvBLoAo5Ny0sogmUa+VuwfbkVeIekMZI2AY6iCN96ehLYsmT5euDTkgQgac/0fSfgZIpbxQdJ2ruP9mZmNgQ2yLBMzgA+n2a/QhGSrwGuqqWziHgY+CJwI7AAuC0iflGPQkvcCOzWO8EH+BqwKdAtaRHwtRScPwZOiYiHgOOA8yWNopj5e2qaDLRzH8cwM7M60/N3E224GNk2LtqOndHsMsysxDJP8Gl5kuZFxIRK2zbkkaWZmdmQ2NAm+FRN0gEUbwUptTQiDhvOxzYzs/obtmEZEddTTKDZqI5tZmb159uwZmZmGQ5LMzOzDIelmZlZhsPSzMwsY9hO8NmY7d4+mi6/p8vMrG48sjQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswzPhh2GFq7ooWPKNc0uw8zK+C+PbLg8sjQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMoZlWEp6Kn3vkPS0pPmSFku6SNKm/bTbTtKNkp6SdG7ZtjdKWijpXknnSFJav2Nqc7ukbkkH99N/h6SjB3FeEyXtW2t7MzOrzcbwQepLImK8pE2AG4APAJf0se8zwH8Ar09fpc4DTgDmAtcCBwK/Br4MXBYR50naLW3r6KP/DuBo4NIaz2Ui8BRwc43tzWyIPHLplPXWTZx7ZsV9Ozs7G1yNDdawHFlWEhFrgT8B7f3s84+I+CNFaD5HUhuwVUTcEhEBXAQc2tsM2Co9Hg081E8Z04C3pZHuSZI2kXSmpD+nUenH0vE+L+kn6fHukhalIP44cFJq/7ayGk+Q1CWpa+3qnqquiZmZVWdjGFkCIGkUsDfw2RqatwPLS5aX83zoTgXmSPo0sAWwfz/9TAFOiYj3pJpOAHoi4k2SRgI3SZoDzAA6JR0GnA58LCIWS/ov4KmIOKu844iYCcwEGNk2Lmo4RzOro+2Pnrbeuk7/ia4N1sYwstxZ0nxgFfBARHTX0IcqrOsNpKOACyJiLHAwcLGkaq/ru4BjUn23AtsB4yJiHTAZuBj4XUTcVEPNZmZWJxvDyLL3Ncs2itHaIRExe4B9LAfGliyP5fnbrcdRvH5JRNySRrBjgJVV9Cvg0xFxfYVt4yhen9xhgLWamVmdbQwjSwAi4mGK26BfrLHtk5L2SbNgjwF+kTY/AOwHIOm1wCjgsT66ehLYsmT5euATvTN0Je0iaQtJo4HvAW8HtpN0eB/tzcxsCGw0YZlcDWxePjmmlKRlwHeByZKWp4k1AJ8AzgfuBZZQzIQFOBk4XtICYBYwOU0CqqQb+JekBZJOSv0tBm6TtAj4IcVo/2zg/0XE3RQj12mSXgb8Ejis0gQfMzNrHPX9e902VCPbxkXbsTOaXYaZlVnmCT4tTdK8iJhQadvGNrI0MzMbsI1hgs96JB0ATC9bvTQiDqtT/7tTzGQttSYi9q5H/2ZmNrQ2yrBMs08rzUCtV/8LgfGN6t/MzIaWb8OamZllOCzNzMwyHJZmZmYZDkszM7OMjXKCz3C3e/touvx+LjOzuvHI0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGX7ryDC0cEUPHVOuaXYZZmYN0Yw/deaRpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FZA0lPpe8dkp6WNF/SYkkXSdq0n3bbSbpR0lOSzi3b9kZJCyXdK+kcSUrrd0xtbpfULengxp6dmZmV82fDDt6SiBgvaRPgBuADwCV97PsM8B/A69NXqfOAE4C5wLXAgcCvgS8Dl0XEeZJ2S9s66n0SZmat6JFLp6y3buLcM1+w3NnZ2fA6PLKsk4hYC/wJaO9nn39ExB8pQvM5ktqArSLilogI4CLg0N5mwFbp8WjgoUp9SzpBUpekrrWrewZ1LmZm9kIeWdaJpFHA3sBna2jeDiwvWV7O86E7FZgj6dPAFsD+lTqIiJnATICRbeOihhrMzFrO9kdPW29dp//qyAZpZ0nzgVXAAxHRXUMfqrCuN/COAi6IiLHAwcDFkvzvZmY2hPxLd/CWRMR44NXAPpIOqaGP5cDYkuWxPH+79TjgMoCIuAUYBYypuVozMxswh2WdRMTDwBTgizW2fVLSPmkW7DHAL9LmB4D9ACS9liIsH6tL0WZmVhWHZX1dDWwu6W197SBpGfBdYLKk5WmGK8AngPOBe4ElFDNhAU4Gjpe0AJgFTE6TgMzMbIh4gk8NIuIl6fsySt4CkkJsj0zbjj7Wd7H+20mIiMXAW2uv1szMBssjSzMzswyPLBtA0gHA9LLVSyPisGbUY2Zmg+OwbICIuB64vtl1mJlZffg2rJmZWYbD0szMLMNhaWZmluGwNDMzy/AEn2Fo9/bRdDXhg4bNzIYrjyzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbht44MQwtX9NAx5Zpml2Fm1hDLmvDWOI8szczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZmZmYZTQ9LSWdL+lzJ8vWSzi9Z/o6kkLRrWbsZkr7QT79flHSvpLskHVCy/sWSZkq6W9JfJL0vrf+8pMWSuiX9VtJOffQ7VdIp6fEFkpZKmi9pgaT9Mud6YqopJI0pWS9J56Rt3ZL2Ktl2kqQ7JC2SNEvSqP6OYWZm9dcKnw17M/B+YIakEcAYYKuS7fsC1wFHAl8FSPsdDry1UoeSdkv7vw7YAfiNpF0iYi1wOrAyInZJ/Wybmt0OTIiI1ZI+AXwbOKKK+k+NiCskvROYCYzrZ9+bgF8BnWXrD0rtxgF7A+cBe0tqBz4D7BYRT0u6LJ3XBVXUZWa2wXrk0il9bps498yK6zs7OxtUTQuMLCkCZN/0+HXAIuBJSdtIGgm8FjiNIiR6vR1YFhH399HnJOBnEbEmIpYC9wJvTts+AnwLICLWRcTj6fGNEbE67TMXGDvA87gFaO9vh4i4PSKW9VHvRVGYC2wtqS1texGwmaQXAZsDD1XqW9IJkrokda1d3TPA0s3MrD9NH1lGxEOS/iVpR4rQ7A2dtwA9QHdEdEtaJ2mPiFhAEZyz+um2nSLwei0H2iVtnZa/JmkisAQ4MSIeLWt/HPDrAZ7KgcDVA2zTqx14sGR5OdAeEV2SzgIeAJ4G5kTEnEodRMRMipEtI9vGRY11mJm1hO2Pntbnts6N+K+O9I4ue8PylpLlm9M+s4Aj0whrEnB5P/2pwrqg+M/BWOCmiNgrHeesFzSUPghMACqP89d3pqT7gJ8C36yyTVX1StqG4lxfSXE7eYtUn5mZDaFWCcubKYJxd4rbsHMpRpb7UgQpFGH5AWB/itHmyn76Ww68omR5LMXty1XAauCqtP5yoHQyzf4Ur2keEhFr0rpvpAk88/s41qnAq4EvAxdWca4DqXd/YGlEPBYRzwJX8vwtazMzGyKtEpY3Ae8BnoiItRHxBLA1RWDeAhARSyjCbhr934IFmE0xCh0p6ZUUE2f+FBEB/BKYmPbbD1gMIGlP4IcUQflcEEfE6RExPiLG93WwiFgHfA8YUTrzdgBmA8ekWbH7AD0R8TDF7dd9JG0uSaneO2vo38zMBqFVwnIhxSzYuWXrenon4CSzgNfw/Miwooi4A7iMIgivAz6VZsJCMVloqqRu4EPAyWn9mcBLgMvTSHL2QE4gBfHXgf7ezvIZScspRo7dJW+RuRa4j2Ii0o+AT6Y+bwWuAG6juB4jSK9LmpnZ0FHxO96Gk5Ft46Lt2BnNLsPMrCGWNWiCj6R5ETGh0rZWGVmamZm1rKa/dWQw0uuD08tWL42Iw5pRTy9JV1HMYC11WkRc34x6zMxscDbosEzh03IB1OywNjOz+vJtWDMzswyHpZmZWYbD0szMLMNhaWZmlrFBT/CxynZvH01XEz5o2MxsuPLI0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsQxHR7BqsziQ9CdzV7Dr6MAZ4vNlFVOC6Bq5Va2vVuqB1a3NdhZ0i4qWVNvjj7oanuyJiQrOLqERSVyvW5roGrlVra9W6oHVrc115vg1rZmaW4bA0MzPLcFgOTzObXUA/WrU21zVwrVpbq9YFrVub68rwBB8zM7MMjyzNzMwyHJZmZmYZDssWJ+lASXdJulfSlArbJemctL1b0l65tpK2lXSDpHvS922GsjZJr5B0o6Q7Jd0h6bMlbaZKWiFpfvo6eKjqStuWSVqYjt1Vsr7Z12zXkmsyX9LfJX0ubRuKa/YaSbdIWiPplGra1uOa1VpXo59jg6ktbWvY82wQ16yhz7Eqa/u/6XnfLelmSXvk2tbrZzMrIvzVol/AJsAS4FXAi4EFwG5l+xwM/BoQsA9wa64t8G1gSno8BZg+xLW1AXulx1sCd5fUNhU4pRnXLG1bBoyp0G9Tr1mFfh6heAP1UF2zlwFvAr5ReqxGPs8GWVfDnmODra2Rz7PB1tWo59gAatsX2CY9Pogh+n1WzZdHlq3tzcC9EXFfRPwT+BkwqWyfScBFUZgLbC2pLdN2EnBhenwhcOhQ1hYRD0fEbQAR8SRwJ9BeQw11rSvTb1OvWdk++wFLIuL+Gmqoqa6IWBkRfwaeHUDbwV6zmutq8HNsULVlNO2alan3c6za2m6OiL+mxbnA2Cra1uNnM8th2dragQdLlpez/g98X/v01/blEfEwFL9UKP6nOZS1PUdSB7AncGvJ6hPTbZif1HBLZbB1BTBH0jxJJ5Ts0zLXDDgSmFW2rtHXrJa2g71mg6nrOQ14jtWjtkY9z+pyzaj/c6yW2o6juMuSa1uPn80sh2VrU4V15e/16WufatoOxmBqKzZKLwF+DnwuIv6eVp8H7AyMBx4GvjPEdb01IvaiuAX0KUlvH+DxG1kbkl4MHAJcXrJ9KK5ZI9o2vO8GPcfqUVujnmf1uGaNeI4NqDZJ76QIy9MG2rZRHJatbTnwipLlscBDVe7TX9tHe2/tpe8rh7g2JG1K8Uvskoi4sneHiHg0ItZGxDrgRxS3X4asrojo/b4SuKrk+E2/ZslBwG0R8WjviiG6ZrW0Hew1G0xdjXyODbq2Bj7PBlVX0ojnWNW1SXoDcD4wKSJWVdG2Hj+bWQ7L1vZnYJykV6b/7R0JzC7bZzZwjAr7AD3pVkR/bWcDx6bHxwK/GMraJAn4MXBnRHy3tEHZ63OHAYuGsK4tJG2Z6tgCeFfJ8Zt6zUq2H0XZ7bEhuma1tB3sNau5rgY/xwZbWyOfZ4P5t+zViOdYVbVJ2hG4EvhQRNxdZdt6/GzmNWLWkL/q90UxO/Juiplgp6d1Hwc+nh4L+EHavhCY0F/btH474LfAPen7tkNZG/BvFLdQuoH56evgtO3itG83xQ9B2xDW9SqKWXYLgDta6ZqlbZsDq4DRZX0OxTXbnuJ/938H/pYeb9Xo51mtdTX6OTbI2hr6PBvkv2XDnmNV1nY+8NeSf7Ou/trW82cz9+WPuzMzM8vwbVgzM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDL+F72d8pGdwVV+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.barh(np.arange(nfeatures), f_importance_ordered, xerr=std[inds], align=\"center\")\n",
    "plt.yticks(np.arange(nfeatures), features_ordered)\n",
    "plt.title('Feature importance')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}