09_RF_class_weighted_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Class weighted random forest model with balanced accuracy grid search"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is another attempt to try to compensate the fact that the data is highly imbalanced. A weighted random forest will try to compensate the imbalance by changing the weight given to each class when computing the impurity score at each split point. In our case we give a weight that is the inverse of the frequency of occurrence for each class in the training dataset. We again use balanced accuracy as the scoring in the grid search"
   ]
  },
  {
   "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": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# suppress anoying iypthon warnings. Not ideal since we suppress also potentially relevant warnings\n",
    "#import warnings\n",
    "#warnings.filterwarnings('ignore')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = '/data/ml_course/05_Capstone_project/rf_data/'\n",
    "model_name = 'RF_class_weighted.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": 4,
   "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": 5,
   "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": 6,
   "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": 7,
   "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": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1], dtype=int8)"
      ]
     },
     "execution_count": 8,
     "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": 9,
   "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": [
    "## Split data into test and training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "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": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 12,
     "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": 13,
   "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: 76.5min\n",
      "[Parallel(n_jobs=-1)]: Done  70 out of  70 | elapsed: 266.6min 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(class_weight='balanced',\n",
       "                                                               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, scoring='balanced_accuracy',\n",
       "             verbose=True)"
      ]
     },
     "execution_count": 13,
     "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='balanced'))\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, scoring='balanced_accuracy', 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": 14,
   "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": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlPElEQVR4nO3de3xU9Z3/8ddnJjcSQoAQLkmAoILINUBEWKQF6wWrK3RXd0Wt219/LbWurd3t77c/97drW/fhtmvtw7W2WsuqrbVV1/68QNWCl4q31mqQO1GMiBBRSbiTkPvn98dM4jAEM4EJM5l5Px+PPDLne86Z+Xzx4XtOvuec7zF3R0REUlcg0QWIiEjvUtCLiKQ4Bb2ISIpT0IuIpDgFvYhIistIdAFdGTJkiJeVlSW6DBGRPmP16tV17l7U1bqkDPqysjIqKysTXYaISJ9hZu8fa52GbkREUlxMQW9mC8zsbTOrNrMbjrHNPDNba2abzOzFiPZtZrYhvE6H6SIiJ1m3QzdmFgTuBM4DaoA3zGy5u2+O2GYgcBewwN23m9nQqLeZ7+518StbRERiFcsY/Uyg2t23ApjZw8BCYHPENlcAj7n7dgB33xXvQkUktbS0tFBTU0NjY2OiS+lTcnJyKC0tJTMzM+Z9Ygn6EmBHxHINcFbUNuOATDNbBeQDP3b3X4XXOfCMmTnwc3df2tWHmNkSYAnAqFGjYu6AiPRNNTU15OfnU1ZWhpklupw+wd3ZvXs3NTU1jBkzJub9Yhmj7+q/QPRMaBnADOAi4ALgRjMbF143x92nAxcCf29mn+nqQ9x9qbtXuHtFUVGXVwiJSAppbGyksLBQId8DZkZhYWGP/wqKJehrgJERy6XAzi62WeHu9eGx+JeAqQDuvjP8exfwOKGhIBERhfxxOJ5/s1iC/g1grJmNMbMs4HJgedQ2y4C5ZpZhZrmEhnaqzCzPzPLDxeUB5wMbe1xlDJpb2/nZqnd5aUttb7y9iEif1W3Qu3srcB2wEqgCHnH3TWZ2jZldE96mClgBrAdeB+5x943AMOAVM1sXbn/K3Vf0Rkcyg8bSl97lqfUf9sbbi4j0ittvv52GhoZe/YyY7ox196eBp6Pa7o5avhW4NaptK+EhnN5mZkwuHcj6D/afjI8TEYmL22+/nauuuorc3Nxe+4yUujN2SkkBWz4+SGNLW6JLEZE+YNu2bYwfP56vfOUrTJo0iSuvvJLnnnuOOXPmMHbsWF5//XXq6+v58pe/zJlnnsm0adNYtmxZ575z585l+vTpTJ8+nT/+8Y8ArFq1innz5nHppZcyfvx4rrzySo71JL877riDnTt3Mn/+fObPnw/AM888w+zZs5k+fTqXXXYZhw4dOuF+JuVcN8drcmkBbe3O5g8PMH3UoESXIyIxuul3m9i880Bc33NC8QC++5cTu92uurqa3/72tyxdupQzzzyTBx98kFdeeYXly5fz/e9/nwkTJnDOOedw3333sW/fPmbOnMm5557L0KFDefbZZ8nJyeGdd95h8eLFnXN0rVmzhk2bNlFcXMycOXN49dVXOfvss4/67G9+85vcdtttvPDCCwwZMoS6ujpuvvlmnnvuOfLy8rjlllu47bbb+M53vnNC/xYpFfRTSgsA2PjBfgW9iMRkzJgxTJ48GYCJEyfyuc99LjQUPHky27Zto6amhuXLl/OjH/0ICF0Wun37doqLi7nuuutYu3YtwWCQLVu2dL7nzJkzKS0tBaC8vJxt27Z1GfTRXnvtNTZv3sycOXMAaG5uZvbs2Sfcx5QK+uEDchjSP5v1NRqnF+lLYjny7i3Z2dmdrwOBQOdyIBCgtbWVYDDIo48+yumnn37Eft/73vcYNmwY69ato729nZycnC7fMxgM0traGlMt7s55553HQw89dCJdOkpKjdGbGVNKC9igoBeROLngggv4yU9+0jnOvmbNGgD279/PiBEjCAQCPPDAA7S1Hd+5wfz8fA4ePAjArFmzePXVV6murgagoaHhiL8UjldKBT3A5JIC3tl1kIbm2L5BRUQ+zY033khLSwtTpkxh0qRJ3HjjjQBce+213H///cyaNYstW7aQl5d3XO+/ZMkSLrzwQubPn09RURG//OUvWbx4MVOmTGHWrFm89dZbJ9wHO9bZ4ESqqKjw433wyHObP+Yrv6rk/10zm4qywXGuTETipaqqijPOOCPRZfRJXf3bmdlqd6/oavvUO6IPn5DVOL2ISEhKnYwFGDYgh2EDstmgG6dEJIl84Qtf4L333jui7ZZbbuGCCy7o9c9OuaAHmFwykPU1+xJdhohIp8cffzxhn51yQzcQup5+a109BxtbEl2KiHyKZDxHmOyO598sJYN+cmkB7rApznfaiUj85OTksHv3boV9D3Q8eCTymv1YpOjQTeiE7Iaa/cw6pTDB1YhIV0pLS6mpqaG2VlOL90THowR7IiWDfkj/bEoG9tNMliJJLDMzs0ePw5Pjl5JDNxA6qt+ooBcRSeGgLy3gvbp69h/WCVkRSW8pG/QdM1lu0lG9iKS5lA36ScXhO2QV9CKS5lI26AflZTFycD/NZCkiaS9lgx5gSslA1n+wL9FliIgkVEoH/eTSAnbsOcze+uZElyIikjApHfRTOm6c0ji9iKSxlA76iQp6EZHUDvqCfpmMGZKnmSxFJK2ldNBD6A5ZXXkjIuks5YN+SmkBO/c3UneoKdGliIgkRExBb2YLzOxtM6s2sxuOsc08M1trZpvM7MWodUEzW2NmT8aj6J6YrHF6EUlz3Qa9mQWBO4ELgQnAYjObELXNQOAu4BJ3nwhcFvU21wNV8Si4pyaWFGCGhm9EJG3FckQ/E6h2963u3gw8DCyM2uYK4DF33w7g7rs6VphZKXARcE98Su6Z/tkZnDIkTw8LF5G0FUvQlwA7IpZrwm2RxgGDzGyVma02s6sj1t0O/BPQ/mkfYmZLzKzSzCrj/SCCKaUD2aA7ZEUkTcUS9NZFW/SzvzKAGYSO3C8AbjSzcWZ2MbDL3Vd39yHuvtTdK9y9oqioKIayYje5pICPDzTx8YHGuL6viEhfEEvQ1wAjI5ZLgZ1dbLPC3evdvQ54CZgKzAEuMbNthIZ8zjGzX59w1T3UMWWxxulFJB3FEvRvAGPNbIyZZQGXA8ujtlkGzDWzDDPLBc4Cqtz9n9291N3Lwvv9wd2vimP9MZlQPICAacpiEUlP3T4z1t1bzew6YCUQBO5z901mdk14/d3uXmVmK4D1hMbi73H3jb1ZeE/kZmUwdmg+G3SHrIikoZgeDu7uTwNPR7XdHbV8K3Drp7zHKmBVjyuMk8mlBax6exfujllXpx1ERFJTyt8Z22FKaQF1h5r5cL9OyIpIekmboO+4Q1bX04tIukmboD9jxACCAdP19CKSdtIm6HMyg4wbls+GDw4kuhQRkZMqbYIeQk+c2lCzD/fo+71ERFJXWgX95NIC9ja0ULP3cKJLERE5adIq6DvvkNWNUyKSRtIq6E8fnk9m0HTljYiklbQK+uyMIOOHD9CVNyKSVtIq6CE0Tr++Zr9OyIpI2ki7oJ9SUsDBxlbe392Q6FJERE6KtAv6yeETsprJUkTSRdoF/bhh+WRlBDSTpYikjbQL+sxggDNGDNCVNyKSNtIu6CE0Tr9p5wHa23VCVkRSX1oG/eTSAg41tfLe7vpElyIi0uvSMuj1DFkRSSdpGfSnFfUnJzOgcXoRSQtpGfQZwQATiwt0h6yIpIW0DHoIPXFq4wcHaNMJWRFJcWkb9FNKCzjc0sa7tYcSXYqISK9K66AHPUNWRFJf2gb9mCH9ycsK6g5ZEUl5aRv0wYAxsbhAc96ISMpL26AHmDoydIfsoabWRJciItJr0jroL5g4nObWdlZu/CjRpYiI9JqYgt7MFpjZ22ZWbWY3HGObeWa21sw2mdmL4bYcM3vdzNaF22+KZ/EnasboQZQO6scTaz9IdCkiIr2m26A3syBwJ3AhMAFYbGYTorYZCNwFXOLuE4HLwquagHPcfSpQDiwws1lxq/4EmRkLy4t5tbqOXQcbE12OiEiviOWIfiZQ7e5b3b0ZeBhYGLXNFcBj7r4dwN13hX+7u3dcqJ4Z/kmqO5QWlZfQ7vDkug8TXYqISK+IJehLgB0RyzXhtkjjgEFmtsrMVpvZ1R0rzCxoZmuBXcCz7v7nrj7EzJaYWaWZVdbW1vaoEydi7LB8JhYP0PCNiKSsWILeumiLPirPAGYAFwEXADea2TgAd29z93KgFJhpZpO6+hB3X+ruFe5eUVRUFGv9cbGovIT1NfvZqrtkRSQFxRL0NcDIiOVSYGcX26xw93p3rwNeAqZGbuDu+4BVwILjLba3/OXUYszgibXR3RIR6ftiCfo3gLFmNsbMsoDLgeVR2ywD5ppZhpnlAmcBVWZWFD5Ri5n1A84F3opb9XEyvCCH2acUsmztB7gn1SkEEZET1m3Qu3srcB2wEqgCHnH3TWZ2jZldE96mClgBrAdeB+5x943ACOAFM1tP6AvjWXd/sne6cmIWlZfw/u4G1u7Yl+hSRETiypLxCLaiosIrKytP6mceaGyh4ubnWHzmSG5a2OVpBBGRpGVmq929oqt1aX1nbKQBOZmce8ZQnlz/IS1t7YkuR0QkbhT0ERaWl7C7vplXqusSXYqISNwo6CPMO72IATkZLFuja+pFJHUo6CNkZwS5aMoIntn8MQ3NmtFSRFKDgj7KovISGprbeHbzx4kuRUQkLhT0Uc4sG0xxQQ6Pa/hGRFKEgj5KIGBcUl7Cy+/UUXeoKdHliIicMAV9FxZNK6at3XlqvWa0FJG+T0HfhfHDBzB+eL5mtBSRlKCgP4ZF00pYs30f7++uT3QpIiInREF/DJd0zGi5RjNaikjfpqA/huKB/ZhZNlgzWopIn6eg/xSLppWwta6eDR/sT3QpIiLHTUH/KT4/aQRZwYCGb0SkT1PQf4qC3Ezmjy/id+t30tau4RsR6ZsU9N1YVF5C7cEm/viuZrQUkb5JQd+N+eOHkp+ToSkRRKTPUtB3IyczyIWThrNy40ccbm5LdDkiIj2moI/BovIS6pvbeK5KM1qKSN+joI/BWacUMnxADss0JYKI9EEK+hgEA8Yl5cWseruWPfXNiS5HRKRHFPQxWlheTGu789QGzWgpIn2Lgj5GE0YMYOzQ/nqerIj0OQr6GJkZi6aVUPn+XnbsaUh0OSIiMVPQ98DC8mKCAeOm322mXXfKikgfoaDvgdJBufzrRWfwXNXH3P78O4kuR0QkJjEFvZktMLO3zazazG44xjbzzGytmW0ysxfDbSPN7AUzqwq3Xx/P4hPhS39RxqUzSrnj+XdYsVEnZkUk+WV0t4GZBYE7gfOAGuANM1vu7psjthkI3AUscPftZjY0vKoV+La7v2lm+cBqM3s2ct++xsy4edEkqncd4h8fWUfZkDzGDx+Q6LJERI4pliP6mUC1u29192bgYWBh1DZXAI+5+3YAd98V/v2hu78Zfn0QqAJK4lV8ouRkBvn5F2fQPzuDJb9azb4GXVsvIskrlqAvAXZELNdwdFiPAwaZ2SozW21mV0e/iZmVAdOAP3f1IWa2xMwqzayytrY2puITadiAHO7+4gw+2t/IdQ+uobWtPdEliYh0KZagty7aoi85yQBmABcBFwA3mtm4zjcw6w88CnzL3Q909SHuvtTdK9y9oqioKKbiE236qEHcvGgSr1TX8YPfv5XockREutTtGD2hI/iREculQPQjl2qAOnevB+rN7CVgKrDFzDIJhfxv3P2xONScVP7mzJFs/vAA977yHhNGDOCvZ5QmuiQRkSPEckT/BjDWzMaYWRZwObA8aptlwFwzyzCzXOAsoMrMDLgXqHL32+JZeDL5l4vOYPYphfzz4xtYt2NfossRETlCt0Hv7q3AdcBKQidTH3H3TWZ2jZldE96mClgBrAdeB+5x943AHOCLwDnhSy/Xmtnne6kvCZMZDHDnldMp6p/N1x5Yza6DjYkuSUSkk7kn3x2eFRUVXllZmegyemzTzv1c+rM/MaF4AA9+9SyyM4KJLklE0oSZrXb3iq7W6c7YOJpYXMCtl01h9ft7+d7yTSTjl6iIpJ9YTsZKD1w8pZjNOw9w16p3mVBcwBdnjU50SSKS5nRE3wu+ff7pzD+9iJuWb+LPW3cnuhwRSXMK+l4QDBg/XjyNUYW5XPubN/lg3+FElyQiaUxB30sG5GTyX1dX0Nzazlfvr6S+qTXRJYlImlLQ96JTi/rzkyum8fbHB/nGQ5omQUQSQ0Hfy+adPpSbLpnIH97axU2/26wrcUTkpNNVNyfBVbNGs2NPAz9/aSujC3P5ytxTEl2SiKQRBf1J8n8WjGfH3gb+/ekqSgf1Y8GkEYkuSUTShIZuTpJAwLjtb8opHzmQ6x9ey5rtexNdkoikCQX9SZSTGeS/rq5g6IBsvnJ/JTv2NCS6JBFJAwr6k2xI/2x+8aWZtLY7X/rF6+xvaEl0SSKS4hT0CXDa0P78/Isz2L6nga/9upLmVl12KSK9R0GfILNOKeSHl07hta17uOHR9brsUkR6ja66SaAvTCtl++7D/OdzWxhVmMu3zh3X/U4iIj2koE+wb37uNLbvaeD2595h5KBcPYpQROJOQZ9gZsYP/moyH+4/zA2Prad4YD9mn1qY6LJEJIVojD4JZGUE+NlVMxhdmMfXHqiketfBRJckIilEQZ8kCvpl8osvnUlWRoAv/eIN3qurT3RJIpIiFPRJZOTgXO75uzM52NjKxXe8zKOra3Q1joicMAV9kikfOZDfXz+XiSUFfPu36/jWf6/lYKNuqhKR46egT0LFA/vx0Fdn8Y/njePJ9R/y+Tte5k3NjSMix0lBn6SCAeObnxvLI1+bRXs7XHb3n7jzhWra2jWUIyI9o6BPcjNGD+bp6+dy4aTh3Lryba685zU+2t+Y6LJEpA9R0PcBBf0y+cniafzw0imsr9nPgh+/xDObPkp0WSLSRyjo+wgz428qRvLkN86mdFA/ljywmn99YgONLW2JLk1EklxMQW9mC8zsbTOrNrMbjrHNPDNba2abzOzFiPb7zGyXmW2MV9Hp7JSi/jz69b/gq3PH8OvXtnPJT1/h7Y90g5WIHFu3QW9mQeBO4EJgArDYzCZEbTMQuAu4xN0nApdFrP4lsCBO9QqQnRHkXy6awP1fnsme+hb+8qevcMfz77C3vjnRpYlIEorliH4mUO3uW929GXgYWBi1zRXAY+6+HcDdd3WscPeXgD1xqlcifHZcESu+NZd544q47dktzP6P5/nXJzawtfZQoksTkSQSy6RmJcCOiOUa4KyobcYBmWa2CsgHfuzuv+pJIWa2BFgCMGrUqJ7smtaG9M9m6dUVvPXRAe575T0eeaOG3/x5O58bP5T/efYpzDplMGaW6DJFJIFiOaLvKiWiL+bOAGYAFwEXADeaWY8mV3f3pe5e4e4VRUVFPdlVgPHDB/DDS6fy6g3n8I1zxvLm9n0s/q/XuPgnr/D4mho9xUokjcUS9DXAyIjlUmBnF9uscPd6d68DXgKmxqdE6Ymi/Gz+8bxx/PGGc/jBX02msaWNf/jvdcz94R+4a1U1+xo0ji+SbmIJ+jeAsWY2xsyygMuB5VHbLAPmmlmGmeUSGtqpim+p0hM5mUEWzxzFs//wWX7xP85k7NB8frjibWb/4A98Z9lGtml2TJG00e0Yvbu3mtl1wEogCNzn7pvM7Jrw+rvdvcrMVgDrgXbgHnffCGBmDwHzgCFmVgN8193v7Z3uSLRAwJh/+lDmnz6Uqg9D4/gPv76DB157nwsnDefrnz2NyaUFiS5TRHqRJeM0uBUVFV5ZWZnoMlLWroON/PLVbTzwp/c52NTK3LFD+Pq8U5l9SqFO3Ir0UWa22t0rulynoE9fBxpb+M1r27n3lfeoO9TE1JEDuXbeqZx3xjACAQW+SF+ioJdP1djSxqNv1vDzF7eyfU8Dpw3tzzWfPZWF5cVkBjVLhkhfoKCXmLS2tfP0xo/42ap3qfrwAMUFOXz1M6fwt2eOJDdLz5EXSWYKeukRd2fVllp+9sK7vL5tD4NyM/nSX4zhirNGUZSfnejyRKQLCno5bpXb9vCzVe/y/Fu7CAaMOacNYVF5MedPHE7/bB3liyQLBb2csOpdh3h8TQ3L1u6kZu9hcjIDnDdhOIvKi/nMuCKN5YskmIJe4sbdWf3+Xp5Y+wFPrf+QvQ0tDMrN5KIpI1hYXsKMUYN0xY5IAijopVe0tLXz0pZalq3dyTObP6KxpZ2Sgf1YWF7MomkljBuWn+gSRdKGgl56XX1TK89s/ogn1uzkleo62tqdU4vymDZqEFNLC5hSOpDxI/LJzggmulSRlKSgl5Oq9mATT63fyYtballfs5/d4QeiZAUDjB+Rz9TSgUwpLWDqyIGcWtSfoIZ6RE6Ygl4Sxt35YN9h1tfsZ13NPtbt2MfGDw5wqKkVgLysIBNLCphaWsC0UYM4b8IwndgVOQ4Kekkq7e3O1rpDrNsRDv+a/VTtPEBzWzt3LJ7GJVOLE12iSJ/zaUGvC6HlpAsEjNOG5nPa0Hz+ekYpAA3NrUz87kqqd+kxiCLxpr+RJSnkZmVQXNCP7bs1T75IvCnoJWmUDcll2+6GRJchknIU9JI0Rhfm8b6O6EXiTkEvSWP04Fz2NrSw/3BLoksRSSkKekkaowvzANiu4RuRuFLQS9IoG5ILwDYN34jElYJeksaowaGg1zi9SHwp6CVp5GZlMGxAtq68EYkzBb0kldGDdeWNSLwp6CWpjC7M5X0d0YvElYJekkrZkDx2HWyiobk10aWIpAwFvSSV0YUdJ2R1VC8SLwp6SSqjB4eupdc4vUj8xBT0ZrbAzN42s2ozu+EY28wzs7VmtsnMXuzJviIdRumIXiTuup2m2MyCwJ3AeUAN8IaZLXf3zRHbDATuAha4+3YzGxrrviKRCvplMjgvS5dYisRRLEf0M4Fqd9/q7s3Aw8DCqG2uAB5z9+0A7r6rB/uKHCF05Y2GbkTiJZagLwF2RCzXhNsijQMGmdkqM1ttZlf3YF+RI5QV5mnoRiSOYgn6rp7cHP38wQxgBnARcAFwo5mNi3Hf0IeYLTGzSjOrrK2tjaEsSVWjBueyc/9hmlrbEl2KSEqIJehrgJERy6XAzi62WeHu9e5eB7wETI1xXwDcfam7V7h7RVFRUaz1SwoqG5KLO+zYczjRpYikhFiC/g1grJmNMbMs4HJgedQ2y4C5ZpZhZrnAWUBVjPuKHKFjumKN04vER7dX3bh7q5ldB6wEgsB97r7JzK4Jr7/b3avMbAWwHmgH7nH3jQBd7dtLfZEUURYOel15IxIf3QY9gLs/DTwd1XZ31PKtwK2x7CvyaQblZpKfnaEHhYvEie6MlaRjZozWg8JF4kZBL0lJDwoXiR8FvSSlssJcavYeprWtPdGliPR5CnpJSqML82htd3bua0x0KSJ9noJektLowXpQuEi8KOglKZUN0bX0IvGioJekNDQ/m5zMgK68EYkDBb0kJTPT5GYicaKgl6Q1arCmKxaJBwW9JK2yIXm8v6eB9vYuJzwVkRgp6CVpjS7Mpbm1nY8O6BJLkROhoJekVdY5i6XG6UVOREyTmokkwujOB4XXM/vUwgRXI72hvqmV2oNN1B1qovZgE7WHmqg72MThljayMgJkBYNkZQTIDBrZGYFQW1R7VkYgtC7c1vGTGTSyI9qCga6eg5QeFPSStEYU9CMzaLrEso9pbGmj7lATdYeaQ+Ed/jkizMOvG5qPfopYwCA7I0hzWzttcTw/EwwYWcGOL4fgJ18cwcARXxpZGUGygoHO9Z3tEV8a2eH9Ot4rK+K9sjMCZAYDR7R17BPdnhk0zHr/C0hBL0krGDBG6sqbpNDa1s7u+ubOoD4qvCNeH2hs7fI9BuVmMqR/NkX52ZSPHEhR/2yG5GdTFG7rWDc4L6vz6Lut3Wlpa6eptZ3m1naa28K/o5cj29vaaGl1mo7ato2WNqe5Nfr92mhube9ct/9wS3ifNprb2mlp9aj3j+/8S5FfBsMG5PD76+fG9f1BQS9JrqwwT0f0vaS93dnb0PzJkfehxnBgNx8V3nsamvEuDq77Z2dQFA7r8cMHcPZpWaHl/CPDuzAvm6yMnp8SDAaMYCBITmYwDj2OD/dPgr/jy6Hji6TjC6Szva3tqC+WlqgvpqbO92onJ6N3+qmgl6Q2ujCXP2/djbuflD9x+zp350Bj6zGPtjuOxjuGVroaGsnOCHQG9ajBuUwfPajzqLsjvIeGf/fLSp4APlnMjOyMINm9FMq9QUEvSW304Fzqm9uoO9RMUX52ostJmIbmo8O7NnIMPHwSs/ZQE82tRw8tZASs8+h6aH42E4sHdB6JRw6fFOVn0z87Q1+qKUZBL0ltdMTkZskU9B1/vje1ttPU0k5jS+jP9qbWts62ztet7TS1fPL6iG1b2o+9X0sb+w+3UHuwifouTlqaQWFeVmeAnzok76ghk44wL+iXSSCNrzpJdwp6SWqRDwqvKBt8xLrWtujg7CI8o4K0c9sjgrjrwG3sbO86sE9Ux5UY2eErQLIzP3mdkxmgIDeLUYV54aPurCOOuov6h05aZgR1K4x0T0EvSa1kYD+CAePmpzbzo5VvHxG6J3rpXeja7HDIZgTIzox4nRGkf3YGhXldBHFEIHe1X2h9ZGgfHeRZwYCOsOWkUdBLUsvKCPDt88ex5aODXYbokUfBR7cda5+sYEBHw5I2FPSS9K6dd1qiSxDp03RIIyKS4hT0IiIpTkEvIpLiFPQiIikupqA3swVm9raZVZvZDV2sn2dm+81sbfjnOxHrrjezjWa2ycy+FcfaRUQkBt1edWNmQeBO4DygBnjDzJa7++aoTV9294uj9p0EfBWYCTQDK8zsKXd/Jy7Vi4hIt2I5op8JVLv7VndvBh4GFsb4/mcAr7l7g7u3Ai8CXzi+UkVE5HjEEvQlwI6I5ZpwW7TZZrbOzH5vZhPDbRuBz5hZoZnlAp8HRnb1IWa2xMwqzayytra2B10QEZFPE8sNU13dpx197/mbwGh3P2RmnweeAMa6e5WZ3QI8CxwC1gFdPpXA3ZcCSwHMrNbM3u+mriFAXQz1pxr1O72o3+nlRPo9+lgrYgn6Go48Ci8FdkZu4O4HIl4/bWZ3mdkQd69z93uBewHM7Pvh9/tU7l7U3TZmVunuFTHUn1LU7/SifqeX3up3LEM3bwBjzWyMmWUBlwPLo4obbuEJrM1sZvh9d4eXh4Z/jwL+CngofuWLiEh3uj2id/dWM7sOWAkEgfvcfZOZXRNefzdwKfB1M2sFDgOXu3c+eOxRMysEWoC/d/e9vdERERHpWkyTmrn708DTUW13R7z+KfDTY+wb/yfdhiztpfdNdup3elG/00uv9Nu8qyf+iohIytAUCCIiKU5BLyKS4vpc0Hc3706qMLP7zGyXmW2MaBtsZs+a2Tvh34MSWWNvMLORZvaCmVWF50e6Ptye0n03sxwzez180+EmM7sp3J7S/e5gZkEzW2NmT4aX06Xf28xsQ3iOsMpwW9z73qeCPmLenQuBCcBiM5uQ2Kp6zS+BBVFtNwDPu/tY4PnwcqppBb7t7mcAs4C/D/83TvW+NwHnuPtUoBxYYGazSP1+d7geqIpYTpd+A8x39/KI6+fj3vc+FfSc2Lw7fYq7vwTsiWpeCNwffn0/sOhk1nQyuPuH7v5m+PVBQv/zl5DiffeQQ+HFzPCPk+L9BjCzUuAi4J6I5pTv96eIe9/7WtDHOu9Oqhrm7h9CKBCBoQmup1eZWRkwDfgzadD38PDFWmAX8Ky7p0W/gduBfwLaI9rSod8Q+jJ/xsxWm9mScFvc+97XHg4ey7w7kgLMrD/wKPAtdz8QvvE6pbl7G1BuZgOBx8PTfKc0M7sY2OXuq81sXoLLSYQ57r4zPIPAs2b2Vm98SF87ou923p0U97GZjQAI/96V4Hp6hZllEgr537j7Y+HmtOg7gLvvA1YROkeT6v2eA1xiZtsIDcWeY2a/JvX7DYC77wz/3gU8Tmh4Ou5972tB3+28OyluOfB34dd/ByxLYC29Ijxn0r1AlbvfFrEqpftuZkXhI3nMrB9wLvAWKd5vd/9ndy919zJC/z//wd2vIsX7DWBmeWaW3/EaOJ/Q1O5x73ufuzM2PA3y7Xwy786/J7ai3mFmDwHzCE1b+jHwXULTPz8CjAK2A5e5e/QJ2z7NzM4GXgY28MmY7f8lNE6fsn03symETrwFCR2APeLu/xaeJypl+x0pPHTzv9z94nTot5mdQugoHkLD6A+6+7/3Rt/7XNCLiEjP9LWhGxER6SEFvYhIilPQi4ikOAW9iEiKU9CLiKQ4Bb2kNTMbaGbXJroOkd6koJd0NxA4KujDM6WKpAQFvaS7/wBODc8H/kZ4LvwHgQ3hScZuDbevN7OvdexkZv87or1j7vg8M3sqPKf8RjP720R1SiRSX5vUTCTebgAmuXt5+M7Mp8LL74VnE9zv7meaWTbwqpk9A4wN/8wkNNHecjP7DFAE7HT3iwDMrODkd0fkaAp6kSO97u7vhV+fD0wxs0vDywWEAv788M+acHv/cPvLwI/M7BbgSXd/+eSVLXJsCnqRI9VHvDbgG+6+MnIDM7sA+IG7/zx6ZzObAXwe+IGZPePu/9ar1YrEQGP0ku4OAvnHWLcS+Hp42mTMbFx4lsGVwJfDc+ZjZiVmNtTMioEGd/818CNgeu+XL9I9HdFLWnP33Wb2qoUewn6Y0EyhHe4ByoA3w9Mn1wKL3P0ZMzsD+FP4gSiHgKuA04BbzawdaAG+fvJ6InJsmr1SRCTFaehGRCTFKehFRFKcgl5EJMUp6EVEUpyCXkQkxSnoRURSnIJeRCTF/X9wR0j5DG8SFAAAAABJRU5ErkJggg==\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": 16,
   "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>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.869553</td>\n",
       "      <td>0.648764</td>\n",
       "      <td>0.006034</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>0.913380</td>\n",
       "      <td>0.631325</td>\n",
       "      <td>0.004110</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>0.935883</td>\n",
       "      <td>0.621396</td>\n",
       "      <td>0.002962</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>0.950982</td>\n",
       "      <td>0.617027</td>\n",
       "      <td>0.003669</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>0.961451</td>\n",
       "      <td>0.612936</td>\n",
       "      <td>0.003046</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>11</td>\n",
       "      <td>0.969070</td>\n",
       "      <td>0.610380</td>\n",
       "      <td>0.002130</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>13</td>\n",
       "      <td>0.974540</td>\n",
       "      <td>0.608806</td>\n",
       "      <td>0.002136</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>15</td>\n",
       "      <td>0.979006</td>\n",
       "      <td>0.607142</td>\n",
       "      <td>0.002585</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>17</td>\n",
       "      <td>0.982240</td>\n",
       "      <td>0.605886</td>\n",
       "      <td>0.001904</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19</td>\n",
       "      <td>0.985004</td>\n",
       "      <td>0.605151</td>\n",
       "      <td>0.001974</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  trees   mean_tr   mean_te    std_te\n",
       "0     1  0.869553  0.648764  0.006034\n",
       "1     3  0.913380  0.631325  0.004110\n",
       "2     5  0.935883  0.621396  0.002962\n",
       "3     7  0.950982  0.617027  0.003669\n",
       "4     9  0.961451  0.612936  0.003046\n",
       "5    11  0.969070  0.610380  0.002130\n",
       "6    13  0.974540  0.608806  0.002136\n",
       "7    15  0.979006  0.607142  0.002585\n",
       "8    17  0.982240  0.605886  0.001904\n",
       "9    19  0.985004  0.605151  0.001974"
      ]
     },
     "execution_count": 16,
     "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": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 65.24%\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": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = grid_RF_cv.predict(X_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn, fp, fn, tp = confusion_matrix(y_te, y_pred).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True positive:  1876\n",
      "True negative:  733874\n",
      "False positive:  3287\n",
      "False negative:  4192\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": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEWCAYAAABLzQ1kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtW0lEQVR4nO3deZxdRZ338c833dn3hYQsIDhEMCBrhKAjAhECjgr6wEwUITpxEERwQx54VIIwYcBRQVZFCatKIgpEEUIMZACHLewkEBJZkpCwJJ2VLKS7f88fp25y09zuvp10p0+6v+/X67zuPXWq6tTtwO/WrVOnjiICMzPLnw6t3QAzMyvNAdrMLKccoM3McsoB2swspxygzcxyygHazCynHKBtm0jqKunPklZK+sM21HOSpPuas22tRdInJM1t7XbYjk+eB90+SPoS8F1gL2A18AwwMSIe3sZ6TwbOBD4WEdXb2s68kxTA8IiY39ptsbbPPeh2QNJ3gcuBi4FBwK7ANcBxzVD9B4CX20NwLoekytZug7UhEeGtDW9Ab2ANcGIDeTqTBfDFabsc6JyOHQ4sAr4HvA0sAb6ajv0YeA/YmM4xHrgAuLWo7t2AACrT/leAV8h68a8CJxWlP1xU7mPAE8DK9PqxomMzgYuAv6d67gMG1PPZCu0/p6j9xwOfBl4GqoD/V5T/YOARYEXKexXQKR17MH2Wd9Pn/bei+v8v8CZwSyEtlfmndI4D0/4QYClweGv/t+Et/5t70G3foUAX4I4G8vwAGAXsD+xHFqR+WHR8Z7JAP5QsCF8tqW9ETCDrlU+OiB4RcX1DDZHUHbgCODYiepIF4WdK5OsH3J3y9gd+DtwtqX9Rti8BXwUGAp2Asxs49c5kf4OhwPnAr4EvAwcBnwDOl/TBlLcG+A4wgOxvNxr4BkBEHJby7Jc+7+Si+vuR/Zo4tfjEEfEPsuD9W0ndgBuAGyNiZgPtNQM8xNEe9AeWRsNDECcBF0bE2xHxDlnP+OSi4xvT8Y0R8Vey3uOeW9meWmAfSV0jYklEzC6R51+AeRFxS0RUR8TvgZeAzxbluSEiXo6IdcAUsi+X+mwkG2/fCNxGFnx/ERGr0/lnA/sCRMSTEfFoOu9rwK+AT5bxmSZExIbUni1ExK+BecBjwGCyL0SzRjlAt33LgAGNjI0OAV4v2n89pW2qo06AXwv0aGpDIuJdsmGB04Alku6WtFcZ7Sm0aWjR/ptNaM+yiKhJ7wsB9K2i4+sK5SV9SNJfJL0paRXZL4QBDdQN8E5ErG8kz6+BfYArI2JDI3nNAAfo9uARYD3ZuGt9FpP9PC/YNaVtjXeBbkX7OxcfjIhpEXEUWU/yJbLA1Vh7Cm16Yyvb1BTXkrVreET0Av4foEbKNDgVSlIPsnH964EL0hCOWaMcoNu4iFhJNu56taTjJXWT1FHSsZJ+krL9HvihpJ0kDUj5b93KUz4DHCZpV0m9gfMKByQNkvS5NBa9gWyopKZEHX8FPiTpS5IqJf0bMAL4y1a2qSl6AquANal3f3qd428BH3xfqYb9AngyIr5GNrb+y21upbULDtDtQET8nGwO9A+Bd4CFwDeBO1OW/wRmAc8BzwNPpbStOdd0YHKq60m2DKodyGaDLCab2fBJ0gW4OnUsAz6T8i4jm4HxmYhYujVtaqKzyS5Aribr3U+uc/wC4CZJKyT9a2OVSToOOIZsWAeyf4cDJZ3UbC22Nss3qpiZ5ZR70GZmOeUAbWaWUw7QZmY55QBtZpZT7XphlwH9KmK3XTq2djOsCV5+rlvjmSxXVrN8aUTstLXlxxzRPZZVlZqN+X5PPrdhWkQcs7Xnypt2HaB326Ujj0/bpbWbYU0wZsj+rd0Ea6K/xe117wptkmVVNTw+bdey8lYMntfYXZ87lHYdoM0s/wKopba1m9EqHKDNLNeCYGOUN8TR1jhAm1nuuQdtZpZDQVDTTu94doA2s9yrbXjBwDbL86DNLNcCqCHK2hojaU9JzxRtqyR9W1I/SdMlzUuvfYvKnCdpvqS5ksYUpR8k6fl07ApJSumdJU1O6Y9J2q2ozLh0jnmSxjXWXgdoM8u9WqKsrTERMTci9o+I/ckeebaW7HFw5wIzImI4MCPtI2kEMBbYm2xVwmskVaTqriV7xNnwtBXmX48HlkfEHsBlwKWprn7ABOAQssfKTSj+IijFAdrMci2AjRFlbU00GvhHRLxO9oT7m1L6TWx+wMVxwG3pcWavAvOBgyUNBnpFxCORLQl6c50yhbpuB0an3vUYYHpEVEXEcmA6m4N6SR6DNrNcizKHL5IBkmYV7V8XEdfVk3cs2cMqAAZFxBKAiFgiaWBKHwo8WlRmUUrbmN7XTS+UWZjqqpa0kuzZoJvSS5QpyQHazPItoKb8zvHSiBjZWCZJnYDPUfTEn/qylm5RvelbW6YkD3GYWa5ldxKWtzXBscBTEVF4ePBbadiC9Pp2Sl8EFK8HMYzsiUCL0vu66VuUSQ9r7k32BKH66qqXA7SZ5ZyoKXNrgi+yeXgDYCpQmFUxDrirKH1smpmxO9nFwMfTcMhqSaPS+PIpdcoU6joBuD+NU08DjpbUN10cPDql1ctDHGaWa9lFwiYF3wZJ6gYcBXy9KPkSYIqk8cAC4ESAiJgtaQowB6gGzojYdN/56cCNQFfgnrRB9vT2WyTNJ+s5j011VUm6CHgi5bswIqoaaqsDtJnlWjYPuvkCdESsJbtoV5y2jGxWR6n8E4GJJdJnAfuUSF9PCvAljk0CJpXbVgdoM8u92mbsQe9IHKDNLNeauwe9I3GANrNcC0RNO53P4ABtZrnnIQ4zsxwKxHtR0XjGNsgB2sxyLbtRxUMcZma55IuEZmY5FCFqwj1oM7NcqnUP2swsf7KLhO0zVLXPT21mOwxfJDQzy7Eaz4M2M8sf30loZpZjtZ7FYWaWP9liSQ7QZma5E4iNvtXbzCx/IvCNKmZm+STfqGJmlkeBe9BmZrnli4RmZjkUqN0u2N8+v5bMbIcRwMaoLGsrh6Q+km6X9JKkFyUdKqmfpOmS5qXXvkX5z5M0X9JcSWOK0g+S9Hw6doUkpfTOkian9Mck7VZUZlw6xzxJ4xprqwO0meWcqClzK9MvgHsjYi9gP+BF4FxgRkQMB2akfSSNAMYCewPHANdIKsz5uxY4FRietmNS+nhgeUTsAVwGXJrq6gdMAA4BDgYmFH8RlOIAbWa5FmR3EpazNUZSL+Aw4HqAiHgvIlYAxwE3pWw3Acen98cBt0XEhoh4FZgPHCxpMNArIh6JiABurlOmUNftwOjUux4DTI+IqohYDkxnc1AvyQHazHKvCT3oAZJmFW2n1qnqg8A7wA2Snpb0G0ndgUERsQQgvQ5M+YcCC4vKL0ppQ9P7uulblImIamAl0L+Buurli4RmlmsRaspaHEsjYmQDxyuBA4EzI+IxSb8gDWfUo9S4STSQvrVlSnIP2sxyLbtIWFHWVoZFwKKIeCzt304WsN9Kwxak17eL8u9SVH4YsDilDyuRvkUZSZVAb6Cqgbrq5QBtZjmXPZOwnK0xEfEmsFDSnilpNDAHmAoUZlWMA+5K76cCY9PMjN3JLgY+noZBVksalcaXT6lTplDXCcD9aZx6GnC0pL7p4uDRKa1eHuIws1zLLhI26zzoM4HfSuoEvAJ8layzOkXSeGABcCJARMyWNIUsiFcDZ0RETarndOBGoCtwT9oguwB5i6T5ZD3nsamuKkkXAU+kfBdGRFVDDXWANrPca847CSPiGaDUOPXoevJPBCaWSJ8F7FMifT0pwJc4NgmYVG5bHaDNLNfa852EDtBmlnt+aKyZWQ5FwMZaB2gzs9zJhjgcoM3McqkJ62y0KQ7QObJwfmcuPm23TftvLujEyd9/k9XLK3hkWm8k6DNgI2dfvoD+O1fz0tPd+MX3s3nvAZz8vTf5+LErAXjgjj7cduUgJOg3aCP/98rX6d2/hl9OGMKzf+8JwIb1YsXSjvzppec3nfPd1R34j0/uxceOWck3L35ju332tqZj51p+9qf5dOwUVFQGD93dh1t+ujNf+9FiRh21io3viSWvd+Jn39mVd1dVUFEZfOenC9njI+uoqAz+9oe+TL5qEF271/CzO+dvqnfA4I3c/8e+/HJCg3cItyktMM1uh9FiAVpSAD+PiO+l/bOBHhFxQZnlLwDWRMRPm3DOvwJfiogVktZERI+mt7z17LLHBq7921wAamrgpAP35uPHrqBH7xrGnfMmAHf+ZgC3XrYz37p0EbvtuY6r7p1LRSUse6uS0z+1J6OOygL0tecP5dczX6J3/xp+c9Fgpt6wEyef/San/XjzjUt3XT+A+S903aINN/9kMB8Z9e52+sRt18YN4pwT/4n1a7Pg+/M75/PE/T156sGeTLp4MLU1YvwPFjP2zLe4fuIQDvvsCjp2Dk4bvSedu9Zy3cyXmHlnX95a1IlvHLXnpnqvuvdlHv5r71b8ZK2h/Q5xtOSn3gB8QdKAFjzHFiLi02llqh3eMw/1ZPAHNjBo2Ea696zdlL5+XQeUOhNdugUV6St244bN6RFAiPXrOhAB766poP/OG993jgfu7Mvhxy/ftD/vua4sf6eSgz65uqU+Vjsi1q/Nbj2u7BhUdAwi4Kn/6UltTfYP9eKT3RkwOPt3iYAu3WrpUBF06lJL9Xti7Zot//ccsvsG+gyo5oXHum/fj5IDtem5hI1tbU1LBuhq4DrgO3UPSPqApBmSnkuvu9ZTxwhJMyW9IumsovJ3SnpS0uzi1aokvbY9vxBa0sy7+nD48Ss27d9wyc6cdNAI7v9TX075/pJN6S891Y3/OHxPvn7knpx16SIqKqGyI5x5yUJOO3IvvnTA3ix4uQtjvrhsi/rfWtSRtxZ2Yv9/XgNAbS1c9+OhfO1HDS4NYE3QoUNwzfS5TH5uNk8/2IO5T28ZWMd8sYon7u8FwEN/6cP6tR34/TOzufWJF7n9lwNZvWLLH7hHHL+c/5nah9Jr7rRd2SyOirK2tqalfzdcDZwkqe5vsquAmyNiX+C3wBX1lN+LbA3VwuLWHVP6v0fEQWR3A50lqX+5DZJ0amEpwneW1TReoBVsfE88el9vDvvsik1pXz33TX775ByO/MJypk7aaVP6Xgeu5dcz53LlPS9z25UDeW+9qN4If7l5AFffN5ffPT2b3T+8jslXDtriHDPv7Ms//8sKKtJ/03++cQAfPXIVA4e+v6dtW6e2VnzjqD056aAR7Ln/Wj6w57pNx7541lvUVMP9f+oDwJ4HrKW2Br50wN6ccshe/J/T3mHnXTdsUd8nj1vBA3f02Y6fIB8KN6qUs7U1LRqgI2IV2ULWZ9U5dCjwu/T+FuCf66ni7rRQ9lKy1aUKUeYsSc8Cj5KtDjW8CW26LiJGRsTInfrn8xv3ift7ssdH1tJ3p+r3HTvi88tLjkHuOnwDXbrV8trcLvxjdjauPGS395Dgk59bwZxZW/be/ueuPlsMb7z4ZDem3jCAUw4ewa8vHMKM2/tx/cTBzfzJ2qd3V1Xw7CM9+OgR2dDRp06s4uBPreLSb36AQm/4iM8vZ9YDPampFiuXdWTOE9340H6bA/oHR6yjoiKY/3y31vgIra69DnFsj1kclwNPATc0kKe+NVGLuxA1QKWkw4FPAYdGxFpJM4Eu29zKHJl5Z98thjfeeKUTQz/4HgCPTuvNLntkf5Y3F3RipyHvUVGZDVks+kcXBg17j+qNYsHLXVixrII+/Wt46sGe7DJ8/ab6Fs7vzJqVlYwYuXZT2rlXL9j0/r7J/Xj52a6M/8HmoRRrmt79qqmuFu+uqqBTl1oO/MQaplw9kJGHr+Jfz3ib739hDzas29w/eueNbLhpxh/70rlrLXsduJY7fr35l9Lhxy9n5l0NPh2pzfIsjhaUVnCaQvacrsIiIf9LtsLTLcBJwMNNqLI32fO+1kraCxjVnO1tbevXiqce6sm3frL5wQvXXzyERf/oTIcOMHDoe5x1afYghxce787kq3ansjIb7zzz4kX07p8N25z03Tc5+/PDqewYDBz6HmdfvjkAz7yzL588bvmmi4rW/PoN2sjZv1hAhw7QoQM8+OfePPa3Xtzw9xfp2Dn4r8n/AOClJ7tzxbnDmHpDf7532UKue2AuKPuSfPXFzTNsDvvsSn508u6t9XFaXXudxaFsmdIWqLhompukQcCrwE8i4oL0lNtJwACyx898NSIW1Cl/AUXT7CS9AHwGWALcSfaomLnATsAFETFT0mvAyIhYWs40u5H7dYnHp+3SUBbLmTFD9m/tJlgT/S1uf7KRp5w0qO9eA+PISSeUlfdPH792m86VNy3Wgy4OjhHxFtCtaP814MhGyl9QZ794Wb9j6ymzW6nzm9mOzUMcZmY55DFoM7Mcc4A2M8shL9hvZpZjbXGOczkcoM0s1yKg2gv2m5nlU3sd4mifX0tmtsNo7rU40qJqz0t6RtKslNZP0nRJ89Jr36L850maL2mupDFF6QeleuZLukLKbv2S1FnS5JT+WLrvo1BmXDrHPEnjGmurA7SZ5V6Eytqa4IiI2L/oppZzgRkRMRyYkfaRNILsrue9gWOAayQVFvG5FjiVbC2g4ek4ZHdNL4+IPYDLgEtTXf2ACcAhbF4ArsH79x2gzSz3tsNiSccBN6X3NwHHF6XflhZtexWYDxwsaTDQKyIeiex27JvrlCnUdTswOvWuxwDTI6IqIpYD09kc1EtygDazXIugKUMcAwrLCaft1FJVAvelNeULxwdFxJLsfLEEGJjShwILi8ouSmlD0/u66VuUiYhqYCXQv4G66uWLhGaWc6Km/FkcS8tYi+PjEbFY0kBguqSXGjz5+0UD6VtbpiT3oM0s95pzDDoiFqfXt4E7yMaD30rDFqTXt1P2RWRrzhcMAxan9GEl0rcoI6mSbAXOqgbqqpcDtJnlWmEtjuaYxSGpu6SehffA0cALwFSgMKtiHHBXej8VGJtmZuxOdjHw8TQMslrSqDS+fEqdMoW6TgDuT+PU04CjJfVNFwePTmn18hCHmeVbpAchN49BwB1pRlwl8LuIuFfSE8AUSeOBBcCJABExO61nP4fsOatnREThWXmnAzcCXYF70gZwPXCLpPlkPeexqa4qSRcBT6R8F0ZEVUONdYA2s9xrrlu9I+IVYL8S6cuA0fWUmQhMLJE+C9inRPp6UoAvcWwSmx9c0igHaDPLtWjaRcI2xQHazHKvhR78lHsO0GaWe028S7DNcIA2s1yLcIA2M8ut9rqanQO0meWex6DNzHIoELWexWFmlk/ttAPtAG1mOeeLhGZmOdZOu9AO0GaWe+5B1yHpShr43oqIs1qkRWZmRQKorXWArmvWdmuFmVl9AnAPeksRcVPxvqTuEfFuyzfJzGxL7XUedKOTCyUdKmkO8GLa30/SNS3eMjOzgihza2PKmf19OdnTaJcBRMSzwGEt2CYzsyLlPe6qLV5ILGsWR0QsTE8gKKipL6+ZWbNrg73jcpQToBdK+hgQkjoBZ5GGO8zMWlxAtNNZHOUMcZwGnAEMBd4A9k/7Zmbbicrc2pZGe9ARsRQ4aTu0xcystHY6xFHOLI4PSvqzpHckvS3pLkkf3B6NMzMDmn0Wh6QKSU9L+kva7ydpuqR56bVvUd7zJM2XNFfSmKL0gyQ9n45doXShTlJnSZNT+mOSdisqMy6dY56kcY21s5whjt8BU4DBwBDgD8Dvy/w7mJltm8KNKuVs5fsWW15LOxeYERHDgRlpH0kjgLHA3sAxwDWSKlKZa4FTgeFpOyaljweWR8QewGXApamufsAE4BDgYGBC8RdBKeUEaEXELRFRnbZbabc/OMysNWSPvWp8K4ekYcC/AL8pSj4OKNycdxNwfFH6bRGxISJeBeYDB0saDPSKiEciIoCb65Qp1HU7MDr1rscA0yOiKiKWA9PZHNRLamgtjn7p7QOSzgVuIwvM/wbc3VClZmbNqnlncVwOnAP0LEobFBFLACJiiaSBKX0o8GhRvkUpbWN6Xze9UGZhqqta0kqgf3F6iTIlNXSR8EmygFz4y3y96FgAFzVUsZlZc1H5v9kHSCpeR+i6iLhuUz3SZ4C3I+JJSYeXc+oSadFA+taWKamhtTh2b6igmdl20bQLgEsjYmQDxz8OfE7Sp4EuQC9JtwJvSRqces+DgbdT/kXALkXlhwGLU/qwEunFZRZJqgR6A1Up/fA6ZWY29GHKetCXpH0k/aukUwpbOeXMzLZdmRcIy7hIGBHnRcSwiNiN7OLf/RHxZWAqUJhVMQ64K72fCoxNMzN2J7sY+HgaDlktaVQaXz6lTplCXSekcwQwDThaUt90cfDolFavRudBS5pAFvVHAH8FjgUeJhsUNzNreS0/LeESYIqk8cAC4ESAiJgtaQowB6gGzoiIwlIXpwM3Al2Be9IGcD1wi6T5ZD3nsamuKkkXAU+kfBdGRFVDjSrnVu8TgP2ApyPiq5IGseXVTzOzllXb/FVGxEzSEENELANG15NvIjCxRPosYJ8S6etJAb7EsUnApHLbWE6AXhcRtZKqJfUiG5vxjSpmtn14wf4GzZLUB/g12cyONcDjLdkoM7NiTZjF0aaUsxbHN9LbX0q6l2xy9nMt2ywzsyIO0FuSdGBDxyLiqZZpkpmZQcM96J81cCyAI5u5Ldvdy891Y8yQ/Vu7GWbWCA9x1BERR2zPhpiZlRQ0963eO4yyHnllZtaq3IM2M8snD3GYmeVVOw3Q5TxRRZK+LOn8tL+rpINbvmlmZkkzP1FlR1HOYknXAIcCX0z7q4GrW6xFZmZFFOVvbU05QxyHRMSBkp4GiIjlkjq1cLvMzDbzLI56bUzP4AoASTvRIkuXmJmV1hZ7x+UoZ4jjCuAOYKCkiWRLjV7coq0yMyvWTsegy1mL47eSniRbik/A8RHxYiPFzMyaRxsdXy5HOQv27wqsBf5cnBYRC1qyYWZmmzhA1+tuNj/wsAuwOzAX2LsF22Vmtona6VWvcoY4PlK8n1a5+3o92c3MrJk0+U7CiHhK0kdbojFmZiV5iKM0Sd8t2u0AHAi802ItMjMr5ouEDepZ9L6abEz6jy3THDOzEtppgG5wHnS6QaVHRPw4bRMj4rfpqbVmZttHM82DltRF0uOSnpU0W9KPU3o/SdMlzUuvfYvKnCdpvqS5ksYUpR8k6fl07ApJSumdJU1O6Y9J2q2ozLh0jnmSxjXW3noDtKTKiKghG9IwM2sVIpvFUc5Whg3AkRGxH7A/cIykUcC5wIyIGA7MSPtIGgGMJZu1dgxwTeq4AlwLnAoMT9sxKX08sDwi9gAuAy5NdfUDJgCHAAcDE4q/CEppqAddeHL3M5KmSjpZ0hcKW1l/CjOzbdWMiyVFZk3a7Zi2AI4DbkrpNwHHp/fHAbdFxIaIeBWYDxwsaTDZA7QfiYgAbq5TplDX7cDo1LseA0yPiKqIWA5MZ3NQL6mcMeh+wDKyZxAW5kMH8KcyypqZbbvyx6AHSJpVtH9dRFxXnCH1gJ8E9gCujojHJA2KiCUAEbFE0sCUfSjwaFHxRSltY3pfN71QZmGqq1rSSqB/cXqJMiU1FKAHphkcL7A5MBe00yF7M2sV5UecpRExssGqsqHb/SX1Ae6QtE8D2Usto1c3HtZt5daUKamhIY4KoEfaeha9L2xmZttFS6wHHRErgJlkwwxvpWEL0uvbKdsiYJeiYsOAxSl9WIn0LcpIqgR6A1UN1FWvhnrQSyLiwoYKm5ltF830mz0tl7wxIlZI6gp8iuwi3lRgHHBJer0rFZkK/E7Sz4EhZBcDH4+IGkmr0wXGx4BTgCuLyowDHgFOAO6PiJA0Dbi46MLg0cB5DbW3oQDdPlfINrN8iWZdi2MwcFMah+4ATImIv0h6BJgiaTywADgRICJmS5oCzCG7D+SMNEQCcDpwI9AVuCdtANcDt0iaT9ZzHpvqqpJ0EfBEyndhRFQ11NiGAvTo8j+zmVkLaqYedEQ8BxxQIn0Z9cS8iJgITCyRPgt43/h1uk/kxHrqmgRMKre99QboxiK7mdn24lu9zczyygHazCyH2ujjrMrhAG1muSY8xGFmllsO0GZmeeUAbWaWUw7QZmY55CeqmJnlmAO0mVk+NeOt3jsUB2gzyz0PcZiZ5ZFvVDEzyzEHaDOz/PGdhGZmOaba9hmhHaDNLN88Bm1mll8e4jAzyysHaDOzfHIP2swsrxygzcxyqHmf6r1D6dDaDTAza0hhHnQ5W6N1SbtIekDSi5JmS/pWSu8nabqkeem1b1GZ8yTNlzRX0pii9IMkPZ+OXSFJKb2zpMkp/TFJuxWVGZfOMU/SuMba6wBtZvkXUd7WuGrgexHxYWAUcIakEcC5wIyIGA7MSPukY2OBvYFjgGskVaS6rgVOBYan7ZiUPh5YHhF7AJcBl6a6+gETgEOAg4EJxV8EpThAm1nuNVcPOiKWRMRT6f1q4EVgKHAccFPKdhNwfHp/HHBbRGyIiFeB+cDBkgYDvSLikYgI4OY6ZQp13Q6MTr3rMcD0iKiKiOXAdDYH9ZI8Br0D6dAhuPLel1m2pCPnj/sgn/jMCk7+3pvsMnwDZ316OPOe6wZAZcdavvWTRQzfdx1RC9eeP5TnHulB5661/OBXrzFkt/eorYFHp/di0sVDWvlTtV3f/fkCDvnUalYsreTrR+4JwAf3XsdZlyyiU5daaqrFVecNY+4z3Tji88s58Rtvbyq7+4fXc8aYD/HK7K5UdqzljIlvsO+ha4gQN16yMw//tU8rfapW0LQbVQZImlW0f11EXFcqYxp6OAB4DBgUEUsgC+KSBqZsQ4FHi4otSmkb0/u66YUyC1Nd1ZJWAv2L00uUKSm3AVrSmojoUbT/FWBkRHyzgTKfA0ZExCWSLgDWRMRPW7yx28nxX1vKwnld6NajBoDXXurChV/bjbMuXbRFvmNPqgLgtNF70rv/Rib+9lXOPHY4AH/85UCe/d8eVHas5dIprzDyiFXMeqDX9v0g7cR9k/sx9YYBfP8Xm/+f/NoPF3Przwcx64FefPTIVYz/4WLOOWEPHrijLw/ckf3a3W2vdVxww2u8MrsrAF/81tusWFrJ+E98GCno2bemVT5Pa2rCRcKlETGy0fqkHsAfgW9HxKo0fFwya4m0aCB9a8uU1KaGOCJiakRc0trtaAkDBr/HwaNXcc/v+m1KWzi/C4v+0eV9eXf90HqefqgnACuXdWTNygo+tN86NqzrwLP/m33nVW/swLznu7LT4I3b5wO0Qy881oPVy7fsA0VA955ZgO3eq4aqtzq+r9wRx69g5p19Nu2PGVvFbVcOTOXFqqrc9qtajGrL28qqS+pIFpx/GxF/SslvpWEL0mvh58wiYJei4sOAxSl9WIn0LcpIqgR6A1UN1FWvHTJAS/psujr6tKS/SRqU0r8i6arWbl9LOO3Hi/nNfw4mauv9pt/kldldOXTMSjpUBIN22cDwfdey05D3tsjTvVcNo45axdMP96inFmsJvzx/KF/70RJunTWH//jRYiZdPPh9eQ773AoeSAG6e68smI87502umvYyP/jVa/QZ0M6+VINmu0iYxoKvB16MiJ8XHZoKFGZVjAPuKkofm2Zm7E52MfDxNByyWtKoVOcpdcoU6joBuD+NU08DjpbUN10cPDql1SvPAbqrpGcKG3Bh0bGHgVERcQBwG3BOuZVKOlXSLEmzNrKheVvcQg751CpWLK1k/vPdyso/7bZ+LF3SkavufZnTL1zMnFndqanZHNg7VATnXfM6d10/gDcXdG6pZlsJnxm3jF9NGMKXR47gVxcM5bs/X7jF8T0PeJcN6zrw+txseKOiMthpyEbmPNGdb475EC8+2Z3/OH9JazS9VTXXRULg48DJwJFF8eXTwCXAUZLmAUelfSJiNjAFmAPcC5wREYUxptOB35BdOPwHcE9Kvx7oL2k+8F3SjJCIqAIuAp5I24UprV55/q20LiL2L+wUxqDT7jBgcvop0gl4tdxK0wWD6wB6qd8OcX/SiI++y6ijV/HR0XPo1Dno1rOGc658nZ+c+YGS+WtrxK8u2Hzt4bKp83jjlc2B+Nv/vZA3Xu3MHb/ZqcXbbls66sQqrv1RdmH2wT/35ts/3TJAH37clsMbq6oqWL+2A3+/pzcAD/2lN8d8cdl2a29uNNP/qRHxMKXHggFG11NmIjCxRPosYJ8S6euBE+upaxIwqdz25rkH3ZArgasi4iPA14H3D8S2ITf812C+PHIE4w4ZwX+d/gGefbhHvcEZoHPXWjp3zb7kDzxsNTXVYsG87E807pwldO9Zyy/P9+yN1rDsrY7se+i7AOz/z2tY/OrmL04p+MRnVjLzrj5FJcSj03ux78fWbCrz+stt+j/392nOG1V2NHnuQTekN/BGet/o3Tht1ceOWck3/vMNevev5qJbXuUfs7vwgy/9E336VzPx968QtbDszY785MxdgexC45e+/TYL5nXm6vteBmDqDQO493f9W/NjtFnnXvM6+x66ht79qrl11hxu+dkgLv/+ME6/cDEVFcF7Gzpw+fc3X2f6yKh3Wbqk4/uGna7/z8Gcc+UCTvvxYlYuq+Rn392l7qnatoh2u2C/ory7b7a7hqbZSTqO7A6dN8jmKH40Ig6vk+cCGplm10v94hCV/FVjZs3kb3H7k+VMfatPzz7D4oDDvlVW3of+fM42nStvctuDLg7Oaf9G4Mb0/i42XzGtL88FLdtCM9te2uLwRTlyG6DNzIDsAmE7HeJwgDaz/Guf8dkB2szyz0McZmY51V5ncThAm1m+NW01uzbFAdrMci27UaV9RmgHaDPLv3b6TEIHaDPLPfegzczyyGPQZmZ51X7X4nCANrP88xCHmVkORZOeSdimOECbWf65B21mllPtMz47QJtZ/qm2fY5xOECbWb4F7fZGlR31mYRm1k6IQFHe1mhd0iRJb0t6oSitn6Tpkual175Fx86TNF/SXEljitIPkvR8OnaFJKX0zpImp/THJO1WVGZcOsc8SWU9qs8B2szyL6K8rXE3AsfUSTsXmBERw4EZaR9JI4CxwN6pzDWSKlKZa4FTgeFpK9Q5HlgeEXuQPZbv0lRXP2ACcAhwMDCh+IugPg7QZpZ/zRSgI+JBoKpO8nHATen9TcDxRem3RcSGiHgVmA8cLGkw0CsiHonsoa431ylTqOt2YHTqXY8BpkdEVUQsB6bz/i+K9/EYtJnlW8uPQQ+KiCUAEbFE0sCUPpTsodQFi1LaxvS+bnqhzMJUV7WklUD/4vQSZerlAG1mudeEWRwDJM0q2r8uIq7b2tOWSIsG0re2TL0coM0s58oeXwZYGhEjm3iCtyQNTr3nwcDbKX0RsEtRvmHA4pQ+rER6cZlFkiqB3mRDKouAw+uUmdlYwzwGbWb5FjTnRcJSpgKFWRXjgLuK0semmRm7k10MfDwNh6yWNCqNL59Sp0yhrhOA+9M49TTgaEl908XBo1Nag9yDNrP8a6YxaEm/J+vJDpC0iGxmxSXAFEnjgQXAiQARMVvSFGAOUA2cERE1qarTyWaEdAXuSRvA9cAtkuaT9ZzHprqqJF0EPJHyXRgRdS9Wvo8DtJnlXnMt2B8RX6zn0Oh68k8EJpZInwXsUyJ9PSnAlzg2CZhUdmNxgDazHYEXSzIzy6EIqGmf93o7QJtZ/rkHbWaWUw7QZmY5FICfSWhmlkcB4TFoM7P8CXyR0MwstzwGbWaWUw7QZmZ5tE3rbOzQHKDNLN8C8ENjzcxyyj1oM7M88q3eZmb5FBCeB21mllO+k9DMLKc8Bm1mlkMRnsVhZpZb7kGbmeVREDU1jWdrgxygzSzfvNyomVmOeZqdmVn+BBDuQZuZ5VB4wX4zs9xqrxcJFe10+gqApHeA11u7HS1gALC0tRthTdKW/80+EBE7bW1hSfeS/X3KsTQijtnac+VNuw7QbZWkWRExsrXbYeXzv5mV0qG1G2BmZqU5QJuZ5ZQDdNt0XWs3wJrM/2b2Ph6DNjPLKfegzcxyygHazCynHKBzSFJI+lnR/tmSLmhC+Qsknd3Ec/5VUp/0fk1TylrD6v49JX1F0lWNlPmcpHPT+yb/e1rb4ACdTxuAL0gqd3L+NouIT0fEiu11PmtYREyNiEtaux3Wuhyg86ma7Kr+d+oekPQBSTMkPZded62njhGSZkp6RdJZReXvlPSkpNmSTi1Kf217fiFYRtJnJT0m6WlJf5M0KKU32su2ts8BOr+uBk6S1LtO+lXAzRGxL/Bb4Ip6yu8FjAEOBiZI6pjS/z0iDgJGAmdJ6t/8Tbc6ukp6prABFxYdexgYFREHALcB57RGAy2fvFhSTkXEKkk3A2cB64oOHQp8Ib2/BfhJPVXcHREbgA2S3gYGAYvIgvLnU55dgOHAsuZuv21hXUTsX9iR9BWyL0iAYcBkSYOBTsCr2711llvuQefb5cB4oHsDeeqbyL6h6H0NUCnpcOBTwKERsR/wNNBlm1tp2+JK4KqI+AjwdfzvYUUcoHMsIqqAKWRBuuB/gbHp/UlkP5HL1RtYHhFrJe0FjGqWhtq26A28kd6Pa82GWP44QOffz9hyqcWzgK9Keg44GfhWE+q6l6wn/RxwEfBos7XSttYFwB8kPUTbXW7UtpJv9TYzyyn3oM3McsoB2swspxygzcxyygHazCynHKDNzHLKAdoaJKkm3aL8gqQ/SOq2DXXdKOmE9P43kkY0kPdwSR/binOUXFOknLVGmrqKn1eZs5bmAG2NWRcR+0fEPsB7wGnFByVVbE2lEfG1iJjTQJbDgSYHaLO2xAHamuIhYI/Uu31A0u+A5yVVSPpvSU+kVfa+DqDMVZLmSLobGFioKK20NzK9P0bSU5KeTSv07Ub2RfCd1Hv/hKSdJP0xneMJSR9PZftLui+tBvcrQI19iPpW9EvHfpbaMkPSTintnyTdm8o8lO7CNGtxXizJyiKpEjiW7G5EyFbJ2yciXk1BbmVEfFRSZ+Dvku4DDgD2BD5CtljTHGBSnXp3An4NHJbq6hcRVZJ+CayJiJ+mfL8DLouIh9MSq9OADwMTgIcj4kJJ/wJsEXDr8e/pHF2BJyT9MSKWka158lREfE/S+anub5It/XpaRMyTdAhwDXDkVvwZzZrEAdoa0zUtkQlZD/p6sqGHxyOisPLa0cC+hfFlsvUlhgOHAb+PiBpgsaT7S9Q/CniwUFdaf6SUT5GtcV3Y7yWpZzrHF1LZuyUtL+Mz1beiXy0wOaXfCvxJUo/0ef9QdO7OZZzDbJs5QFtjtlgqEyAFqneLk4AzI2JanXyfpv7V9orLlrPeQAeyVfiKl14ttKXs9QrqrOi3VtJM6l9BLtJ5V9T9G5htDx6DtuYwDTi98FAASR+S1B14EBibxqgHA0eUKPsI8ElJu6ey/VL6aqBnUb77yIYbSPn2T28fJFvVD0nHAn0baWtDK/p1AAq/Ar5ENnSyCnhV0onpHJK0XyPnMGsWDtDWHH5DNr78lKQXgF+R/Tq7A5gHPA9cC/xP3YIR8Q7ZuPGfJD3L5iGGPwOfL1wkJFvFb2S6CDmHzbNJfgwcJukpsqGWBY20taEV/d4F9pb0JNkYc+HJJycB41P7ZgPHlfE3MdtmXs3OzCyn3IM2M8spB2gzs5xygDYzyykHaDOznHKANjPLKQdoM7OccoA2M8up/w+JjHI14v2xWAAAAABJRU5ErkJggg==\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": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of Detection (POD): 30.91628213579433\n",
      "False Alarm Rate (FAR): 0.44589987804563724\n",
      "Positive Predictive Value (PPV): 0.36335463877590546\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": [
    "The POD is actually worst than with a classical RF model optimized for balanced accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "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": 24,
   "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": 25,
   "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": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAEICAYAAAAwft9dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfE0lEQVR4nO3df5zVVb3v8dcbVBBTNKkOgjplaFok6uQP+uW9evxZojeP4qmUMum3RxMV09PhmiWkJZlZkXnN8pI/UqOridVRK0VyUBiQMkFQQR6hWHMklBA+94/vmtjO2TNrz8ye2Xtm3s/HYx6zv2t91/qutTfO2/Xda88oIjAzM7P2Dar1AMzMzOqdw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZWNZK+KOm6Wo/DrNrkz1ma1QdJK4E3AZtLiveOiOe62ecnIuJX3Rtd3yNpGvDWiPhIrcdifZ9Xlmb15YMR8bqSry4HZTVI2qaW1++qvjpuq18OS7M6J2m4pB9IWiNptaTLJA1OdXtJ+k9J6yS9IOkmSTunuh8BewA/l7Re0gWSDpe0qk3/KyUdmR5Pk3SbpB9L+i9gUkfXLzPWaZJ+nB43SApJH5P0rKS/SPqUpHdJapb0V0nXlLSdJOlBSd+S1CLpj5KOKKnfTdIcSS9KWibprDbXLR33p4AvAqemuS9K531M0h8kvSTpKUmfLOnjcEmrJJ0naW2a78dK6reX9HVJT6fx/U7S9qnuUEkPpTktknR4F15qq2MOS7P690PgVeCtwAHAUcAnUp2Ay4HdgH2B3YFpABHxUeAZtq5Wv1bh9SYAtwE7Azdlrl+JQ4AxwKnATOBi4Ejg7cApkt7f5tyngBHAfwC3S3p9qpsNrEpzPRn4ammYthn3D4CvAjenue+fzlkLfADYCfgYcJWkA0v6+CdgODAKOBP4tqRdUt2VwEHAeOD1wAXAFkmjgLuAy1L5FOCnkt7QiefI6pzD0qy+3JlWJ3+VdKekNwHHAudExN8iYi1wFTARICKWRcQvI2JjRDwPfAN4f/vdV2ReRNwZEVsoQqXd61foyxHxSkTcC/wNmB0RayNiNfBbigButRaYGRGbIuJm4AngeEm7A+8BLkx9LQSuAz5abtwR8XK5gUTEXRGxPAoPAPcC7y05ZRNwabr+3cB6YB9Jg4CPA/8WEasjYnNEPBQRG4GPAHdHxN3p2r8EmoDjOvEcWZ3zfX2z+nJi6WYcSQcD2wJrJLUWDwKeTfVvBK6m+IG/Y6r7SzfH8GzJ4z07un6F/lzy+OUyx68rOV4dr911+DTFSnI34MWIeKlNXWM74y5L0rEUK9a9KeYxDFhccsq6iHi15HhDGt8IYCiwvEy3ewL/IumDJWXbAvflxmN9h8PSrL49C2wERrT5Id7qciCAd0bEOkknAteU1Lfd7v43ioAAIL332PZ2YWmb3PWrbZQklQTmHsAc4Dng9ZJ2LAnMPYDVJW3bzvU1x5KGAD8FTgd+FhGbJN1JcSs75wXgFWAvYFGbumeBH0XEWf+tlfUbvg1rVsciYg3FrcKvS9pJ0qC0qaf1VuuOFLcK/5reOzu/TRd/Bt5ScvwnYKik4yVtC1wCDOnG9avtjcDZkraV9C8U78PeHRHPAg8Bl0saKumdFO8p3tRBX38GGtItVIDtKOb6PPBqWmUeVcmg0i3p64FvpI1GgyUdlgL4x8AHJR2dyoemzUKjOz99q1cOS7P6dzrFD/qlFLdYbwNGprr/DRwItFBsMrm9TdvLgUvSe6BTIqIF+AzF+32rKVaaq+hYR9evtvkUm4FeAL4CnBwR61LdaUADxSrzDuA/0vuD7bk1fV8n6dG0Ij0buIViHv9KsWqt1BSKW7aPAC8CM4BBKcgnUOy+fZ5ipXk+/vnar/iXEphZXZA0ieIXKLyn1mMxa8v/52NmZpbhsDQzM8vwbVgzM7MMryzNzMwy/DnLfmjEiBHR0NBQ62GYmfUpCxYseCEiyv6aQodlP9TQ0EBTU1Oth2Fm1qdIerq9Ot+GNTMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4V9K0A8tXt1Cw9S7aj2MAWvl9ONrPQQzqzKvLM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzjGxYSrpK0jklx3MlXVdy/HVJIWmfNu1mSrqgg34vkrRM0hOSji4p307SLEl/kvRHSR9K5V+QtFRSs6RfS9qznX6nSZqSHt8gaYWkhZIWSToiM9fPpTGFpBEl5ZJ0daprlnRgSd25kh6XtETSbElDO+h/kqTdOhpDZnxf7GpbMzPrukpWlg8B4wEkDQJGAG8vqR8P3ANMbC1I550M3FyuQ0n7pfPfDhwDXCtpcKq+GFgbEXsD+wEPpPLHgMaIeCdwG/C1CsYOcH5EjAPOAb6bOfdB4Ejg6TblxwJj0tdk4DtpHqOAs9O43gEMpuR5KGMS0OWwBByWZmY1UElYPkgKS4pwWwK8JGkXSUOAfYELeW1IvA9YGRFtQ6fVBOAnEbExIlYAy4CDU93HgcsBImJLRLyQHt8XERvSOQ8DoyuZYIl5wKiOToiIxyJiZTvjvTEKDwM7SxqZ6rYBtpe0DTAMeK5c35JOBhqBm9JKd3tJB0l6QNKCtGIfKWl4Wm3vk9rNlnSWpOnpOgsl3VSm/8mSmiQ1bd7QUtkzYmZmFcmGZUQ8B7wqaQ+K0JwHzAcOo/jh3xwRzcAWSfunZhOB2R10Owp4tuR4FTBK0s7p+MuSHpV0q6Q3lWl/JvCL3NjbOAa4s5NtWpUdb0SsBq4EngHWAC0RcW+5DiLiNqAJ+HBa6b4KfAs4OSIOAq4HvhIRLcDngBskTQR2iYjvR8RU4OWIGBcRHy7T/6yIaIyIxsHDhndxmmZmVk6lG3xaV5etYTmv5PihdM5sYGJaYU0Abu2gP5UpC4pV2mjgwYg4MF3nytc0lD5CEdJXVDj2KyQ9BfwY+GqFbSoar6RdKOb6Zorbqzuk8VViH+AdwC8lLQQuIa2WI+KXwGLg28AnujhmMzOrkkrDsvV9y7EUt2EfplhZjqcIUijC8hSK9/yaI2JtB/2tAnYvOR5NcftyHbABuCOV3wqUbqY5kuI9zRMiYmMq+0q6NbmwnWudD7yVIox+WMFcOzPeI4EVEfF8RGwCbmfrLescAY+nleK4iBgbEUfBP97z3Rd4GXh9F8dsZmZV0pmV5QeAFyNic0S8COxMEZjzACJiOUXYTafjW7AAcyhWoUMkvZli48zvIyKAnwOHp/OOAJYCSDoA+B5FUP4jiCPi4tbAae9iEbEF+CYwqHTnbSfMAU5Pu2IPpbjduobi9uuhkoZJUhrvHzro5yVgx/T4CeANkg5L89tWUuvGqXNTP6cB10vaNpVvKnlsZma9pNKwXEyxC/bhNmUtrRtwktnA29i6MiwrIh4HbqEIwnuAz0bE5lR9ITBNUjPwUeC8VH4F8Drg1rSSnFPh2FuvGcBlQEcfZzlb0iqKlWNzyUdk7gaeotiI9H3gM6nP+RQ7cx+leD4GAbM6GMYNwHfTKngwxY7hGZIWAQuB8ZL2prj1el5E/Bb4DcWqmNR3c7kNPmZm1nNUZIj1J0NGjomRZ8ys9TAGrJXTj6/1EMysCyQtiIjGcnX+DT5mZmYZ2/Rk5+n9wRltildExEk9ed0cSXdQ7GAtdWFEzK1S/98G3t2m+JsR8X+q0b+ZmfWuHg3LFD5VCaBq6umwjojP9mT/ZmbWu3wb1szMLMNhaWZmluGwNDMzy3BYmpmZZfToBh+rjbGjhtPkz/qZmVWNV5ZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbh3bD90OLVLTRMvavWw7Au8F8sMatPXlmamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhn9MiwlrU/fGyS9LGmhpKWSbpS0bQftdpV0n6T1kq5pU3eQpMWSlkm6WpJS+R6pzWOSmiUd10H/DZL+tRvzOlzS+K62NzOzrumXYdnG8ogYB4wFRgOndHDuK8C/A1PK1H0HmAyMSV/HpPJLgFsi4gBgInBtB/03AF0OS+BwwGFpZtbLBkJYAhARm4HfA6M6OOdvEfE7itD8B0kjgZ0iYl5EBHAjcGJrM2Cn9Hg48FwHw5gOvDetdM+VNFjSFZIeSavST6brfUHS9enxWElLJO0HfAo4N7V/b5sxTpbUJKlp84aWip4TMzOrzID5E12ShgKHAP/WheajgFUlx6vYGrrTgHslfR7YATiyg36mAlMi4gNpTJOBloh4l6QhwIOS7gVmAvdLOgm4GPhkRCyV9F1gfURc2bbjiJgFzAIYMnJMdGGOZmbWjoGwstxL0kJgHfBMRDR3oQ+VKWsNpNOAGyJiNHAc8CNJlT6vRwGnp/HNB3YFxkTEFmAS8CPggYh4sAtjNjOzKhkIK8vlETEu3Uq9X9IJETGnk32soni/s9Vott5uPZP0/mVEzEsr2BHA2gr6FfD5iJhbpm4MsB7YrZNjNTOzKhsIK0sAImINxW3Qi7rY9iVJh6ZdsKcDP0vVzwBHAEjaFxgKPN9OVy8BO5YczwU+3bpDV9LeknaQNBz4JvA+YFdJJ7fT3szMesGACcvkTmBY280xpSStBL4BTJK0Km2sAfg0cB2wDFgO/CKVnwecJWkRMBuYlDYBldMMvCppkaRzU39LgUclLQG+R7Havwq4NiL+RLFynS7pjcDPgZPKbfAxM7Oeo/Z/rltfNWTkmBh5xsxaD8O6YOX042s9BLMBS9KCiGgsVzfQVpZmZmadNhA2+Pw3ko4GZrQpXhERJ1Wp/7EUO1lLbYyIQ6rRv5mZ9a4BGZZp92m5HajV6n8xMK6n+jczs97l27BmZmYZDkszM7MMh6WZmVmGw9LMzCxjQG7w6e/GjhpOkz+vZ2ZWNV5ZmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhnfD9kOLV7fQMPWuWg/DrCz/ZRXri7yyNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZdIGl9+t4g6WVJCyUtlXSjpG07aLerpPskrZd0TZu6gyQtlrRM0tWSlMr3SG0ek9Qs6bienZ2ZmbXlsOy+5RExDhgLjAZO6eDcV4B/B6aUqfsOMBkYk76OSeWXALdExAHARODa6gzbzMwq5bCskojYDPweGNXBOX+LiN9RhOY/SBoJ7BQR8yIigBuBE1ubATulx8OB58r1LWmypCZJTZs3tHRrLmZm9loOyyqRNBQ4BLinC81HAatKjlexNXSnAR+RtAq4G/h8uQ4iYlZENEZE4+Bhw7swBDMza4/Dsvv2krQQWAc8ExHNXehDZcoifT8NuCEiRgPHAT+S5NfNzKwX+Ydu97W+Z/lW4FBJJ3Shj1UU73e2Gs3W261nArcARMQ8YCgwosujNTOzTnNYVklErAGmAhd1se1Lkg5Nu2BPB36Wqp8BjgCQtC9FWD5flUGbmVlFHJbVdScwTNJ72ztB0krgG8AkSask7ZeqPg1cBywDlgO/SOXnAWdJWgTMBialTUBmZtZLtqn1APqiiHhd+r4SeEdJeQD7Z9o2tFPeVNpXSflS4N1dH62ZmXWXV5ZmZmYZXln2AElHAzPaFK+IiJNqMR4zM+seh2UPiIi5wNxaj8PMzKrDt2HNzMwyHJZmZmYZDkszM7MMh6WZmVmGN/j0Q2NHDadp+vG1HoaZWb/hlaWZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWV4N2w/tHh1Cw1T76r1MMx61UrvALce5JWlmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpZRV2EpaX2b40mSrkmPp0laLWmhpKWSTkvlKyTt06bdTEkX9N7IKyNpnKTjutH+REn7VXNMZmaWV1dhWYGrImIcMAH4nqRtgZ8AE1tPkDQIOBm4uScGIGlwN5qPA7oclsCJgMPSzKyX9bWwBCAingQ2ALsAsykJS+B9wMqIeLpc27RavV3SPZKelPS1krrTJC2WtETSjJLy9ZIulTQfOCwdz5C0QNKvJB0s6X5JT0k6oZ3rbgdcCpyaVsenStpB0vWSHpH0mKQJ6dyrJX0pPT5a0m8kjQdOAK5I7fdq0/9kSU2SmjZvaOn0c2pmZu2rtz/Rtb2khSXHrwfmtD1J0oHAkxGxFlgraYuk/SNiEUVwzs5cZxxwALAReELSt4DNwAzgIOAvwL2SToyIO4EdgCUR0RpgOwD3R8SFku4ALgP+mWLV98NyY46Iv6cAbIyIz6V+vgr8Z0R8XNLOwO8l/QqYCjwi6bfA1cBxEbFc0hzg/0XEbWX6nwXMAhgyckxk5m9mZp1Qb2H5crrNChSrQKCxpP5cSWcBbwGOKSmfDUyU9DjFLdovZa7z64hoSddYCuwJ7EoRgM+n8psoVql3UgTpT0va/x24Jz1eDGyMiE2SFgMNFc4V4CjgBElT0vFQYI+I+EOa52+AcyNieSf6NDOzKqu3sMy5KiKulPS/gBsl7RURr1CE5b3AA0BzWnF2ZGPJ480Uz4M6OP+ViNhccrwpIlpXb1ta+4uILZI685wK+FBEPFGmbiywDtitE/2ZmVkP6KvvWd4ONAFnpOPlFMEynfwt2PbMB94vaUTaxHMaRfhW00vAjiXHc4HPSxKApAPS9z2B8yhuFR8r6ZB22puZWS/ok2GZXAp8Ie1+hSIk3wbc0ZXOImINcBFwH7AIeDQiflaNgZa4D9ivdYMP8GVgW6BZ0hLgyyk4fwBMiYjngDOB6yQNpdj5e37aDLRXO9cwM7Mq09a7idZfDBk5JkaeMbPWwzDrVSunH1/rIVgfJ2lBRDSWq+vLK0szM7Ne0dc2+FRM0tEUHwUptSIiTurP1zYzs+rrt2EZEXMpNtAMqGubmVn1+TasmZlZhsPSzMwsw2FpZmaW4bA0MzPL6LcbfAaysaOG0+TPnJmZVY1XlmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy/BHR/qhxatbaJh6V62HYWZWsXr/E2teWZqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMmoelpKuknROyfFcSdeVHH9dUkjap027mZIu6KDfiyQtk/SEpKNLyreTNEvSnyT9UdKHUvkXJC2V1Czp15L2bKffaZKmpMc3SFohaaGkRZKOyMz1c2lMIWlESbkkXZ3qmiUdWFJ3rqTHJS2RNFvS0I6uYWZm1VfzsAQeAsYDSBoEjADeXlI/HrgHmNhakM47Gbi5XIeS9kvnvx04BrhW0uBUfTGwNiL2BvYDHkjljwGNEfFO4DbgaxWO//yIGAecA3w3c+6DwJHA023KjwXGpK/JwHfSPEYBZ6dxvQMYTMnzYGZmvaMewvJBUlhShNsS4CVJu0gaAuwLXMhrQ+J9wMqIaBs6rSYAP4mIjRGxAlgGHJzqPg5cDhARWyLihfT4vojYkM55GBjdyXnMA0Z1dEJEPBYRK9sZ741ReBjYWdLIVLcNsL2kbYBhwHPl+pY0WVKTpKbNG1o6OXQzM+tIzcMyIp4DXpW0B0VozgPmA4cBjUBzRDQDWyTtn5pNBGZ30O0o4NmS41XAKEk7p+MvS3pU0q2S3lSm/ZnALzo5lWOAOzvZplXZ8UbEauBK4BlgDdASEfeW6yAiZkVEY0Q0Dh42vIvDMDOzcmoelknr6rI1LOeVHD+UzpkNTEwrrAnArR30pzJlQbFKGw08GBEHputc+ZqG0kcoQvqKCsd+haSngB8DX62wTUXjlbQLxVzfDOwG7JDGZ2ZmvahewrL1fcuxFLdhH6ZYWY6nCFIowvIUivf8miNibQf9rQJ2LzkeTXH7ch2wAbgjld8KlG6mOZLiPc0TImJjKvtK2sCzsJ1rnQ+8FbgE+GEFc+3MeI8EVkTE8xGxCbidrbeszcysl9RLWD4IfAB4MSI2R8SLwM4UgTkPICKWU4TddDq+BQswh2IVOkTSmyk2zvw+IgL4OXB4Ou8IYCmApAOA71EE5T+COCIujohxaRNPWRGxBfgmMKh0520nzAFOT7tiD6W43bqG4vbroZKGSVIa7x+60L+ZmXVDvYTlYopdsA+3KWtp3YCTzAbextaVYVkR8ThwC0UQ3gN8NiI2p+oLgWmSmoGPAuel8iuA1wG3ppXknM5MIAXxZUBHH2c5W9IqipVjc8lHZO4GnqLYiPR94DOpz/kUO3MfpXg+BgGzOjMuMzPrPhU/460/GTJyTIw8Y2ath2FmVrGV04+v9RCQtCAiGsvV1cvK0szMrG5tU+sBdEd6f3BGm+IVEXFSLcbTStIdFDtYS10YEXNrMR4zM+uePh2WKXzqLoBqHdZmZlZdvg1rZmaW4bA0MzPLcFiamZllOCzNzMwy+vQGHytv7KjhNNXBZ5bMzPoLryzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhj470Q4tXt9Aw9a5aD8PMrFf15J/58srSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpYxoMJS0vo2x5MkXZMeT5O0WtJCSUslnZbKV0jap027mZIu6L2Rm5lZLQ2osKzAVRExDpgAfE/StsBPgImtJ0gaBJwM3NwTA5A0uCf6NTOzrnNYlhERTwIbgF2A2ZSEJfA+YGVEPF2ubVqt3i7pHklPSvpaSd1pkhZLWiJpRkn5ekmXSpoPHJaOZ0haIOlXkg6WdL+kpySd0M51J0tqktS0eUNLFZ4FMzNrNdDCcvt0m3WhpIXApeVOknQg8GRErI2IZmCLpP1T9USKAO3IOOBUYCxwqqTdJe0GzAD+Z6p/l6QT0/k7AEsi4pCI+F06vj8iDgJeAi4D/hk4qb0xR8SsiGiMiMbBw4bnngczM+uEgfYnul5Ot1mBYhUINJbUnyvpLOAtwDEl5bOBiZIep7hF+6XMdX4dES3pGkuBPYFdKQLw+VR+E8Uq9U5gM/DTkvZ/B+5JjxcDGyNik6TFQEOFczUzsyoZaCvLnKsiYh+KVeGNkoam8tnAKcCRQHNErM30s7Hk8WaK/ylRB+e/EhGbS443RUSkx1ta+4uILQy8/8ExM6s5h2UZEXE70ASckY6XA+uA6eRvwbZnPvB+SSPSJp7TgAeqMFwzM+thDsv2XQp8Ie1+hSIk3wbc0ZXOImINcBFwH7AIeDQiflaNgZqZWc/S1rt91l8MGTkmRp4xs9bDMDPrVSunH9+t9pIWRERjuTqvLM3MzDK8WaSLJB1N8VGQUisi4qRajMfMzHqOw7KLImIuMLfW4zAzs57n27BmZmYZDkszM7MMh6WZmVmGw9LMzCzDG3z6obGjhtPUzc8bmZnZVl5ZmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllKCJqPQarMkkvAU/Uehy9ZATwQq0H0YsG0nw91/6pnue6Z0S8oVyFf91d//RERDTWehC9QVLTQJkrDKz5eq79U1+dq2/DmpmZZTgszczMMhyW/dOsWg+gFw2kucLAmq/n2j/1ybl6g4+ZmVmGV5ZmZmYZDkszM7MMh2UfI+kYSU9IWiZpapl6Sbo61TdLOrDStvWmm3NdKWmxpIWSmnp35J1XwVzfJmmepI2SpnSmbb3p5lz71OsKFc33w+nfb7OkhyTtX2nbetPNudb3axsR/uojX8BgYDnwFmA7YBGwX5tzjgN+AQg4FJhfadt6+urOXFPdSmBEredRxbm+EXgX8BVgSmfa1tNXd+ba117XTsx3PLBLenxsP/9vtuxc+8Jr65Vl33IwsCwinoqIvwM/ASa0OWcCcGMUHgZ2ljSywrb1pDtz7Wuyc42ItRHxCLCps23rTHfm2hdVMt+HIuIv6fBhYHSlbetMd+Za9xyWfcso4NmS41WprJJzKmlbT7ozV4AA7pW0QNLkHhtldXTntemPr2tH+tLrCp2f75kUd0u60rbWujNXqPPX1r/urm9RmbK2n/1p75xK2taT7swV4N0R8ZykNwK/lPTHiPhNVUdYPd15bfrj69qRvvS6QifmK+l/UATIezrbtk50Z65Q56+tV5Z9yypg95Lj0cBzFZ5TSdt60p25EhGt39cCd1DcIqpX3Xlt+uPr2q4+9rpChfOV9E7gOmBCRKzrTNs60p251v1r67DsWx4Bxkh6s6TtgInAnDbnzAFOTztFDwVaImJNhW3rSZfnKmkHSTsCSNoBOApY0puD76TuvDb98XUtqw++rlDBfCXtAdwOfDQi/tSZtnWmy3PtC6+tb8P2IRHxqqTPAXMpdp5dHxGPS/pUqv8ucDfFLtFlwAbgYx21rcE0KtKduQJvAu6QBMW/8f8bEff08hQqVslcJf0T0ATsBGyRdA7FTsP/6m+va3tzpfjTTn3mdYWK/x1/CdgVuDbN7dWIaOyn/82WnSt94L9Z/7o7MzOzDN+GNTMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMws4/8DOUgDZgLH1+sAAAAASUVORK5CYII=\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
}