10_RF_bootstrap_class_weighted_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bootstrap class weighted random forest model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a similar model as the weighted random forest. The difference is that the class weight is computed for each bootstrap sample instead of over the entire training dataset"
   ]
  },
  {
   "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_bootstrap_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: 85.3min\n",
      "[Parallel(n_jobs=-1)]: Done  70 out of  70 | elapsed: 280.7min 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_subsample',\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_subsample'))\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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEHCAYAAACgHI2PAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlc0lEQVR4nO3de3TV5Z3v8fd379xDIJAEEAIJSCI3FSEiKlSsRbB6dJyprZfWnt4Yq/Z0ujrt2LOOp7XrdI0z2hkdR8dxrNU6VtuOWql1AG1F8IIYFFGIYOQaEQgIAZJAbt/zx97gJiZkB3b47ez9ea3FSn63vb8PXX749fk9v+cxd0dERFJXKOgCRESkbynoRURSnIJeRCTFKehFRFKcgl5EJMUp6EVEUlxGPCeZ2TzgbiAMPOjut3dxzmzgLiAT2OXuF0T3bwL2A+1Am7tX9fR9xcXFXl5eHk9pIiICrFy5cpe7l3R1rMegN7MwcC8wB6gD3jCzBe6+NuacQuA+YJ67bzGzoZ0+5kJ33xVvweXl5VRXV8d7uohI2jOzzd0di6frZjpQ6+4b3L0FeAK4otM51wJPufsWAHffebzFiohIYsUT9COBrTHbddF9sSqBwWa2xMxWmtn1McccWBzdP//EyhURkd6Kp4/eutjXed6EDGAacBGQC7xmZsvdfT1wvrtvi3bnPG9m77n70k99SeQfgfkAo0eP7k0bRETkGOIJ+jpgVMx2KbCti3N2uXsj0GhmS4EzgfXuvg0i3Tlm9jSRrqBPBb27PwA8AFBVVaUJeERSXGtrK3V1dRw8eDDoUvqVnJwcSktLyczMjPuaeIL+DaDCzMYAHwJXE+mTj/UM8K9mlgFkAecA/2xm+UDI3fdHf78Y+Gnc1YlIyqqrq6OgoIDy8nLMuuo4kM7cnd27d1NXV8eYMWPivq7HoHf3NjO7GVhEZHjlQ+6+xsxuiB6/391rzGwhsBroIDIE810zGws8Hf0fMQP4tbsv7HXrRCTlHDx4UCHfS2ZGUVER9fX1vbournH07v4c8Fynffd32r4DuKPTvg1EunBERD5FId97x/N3ljJvxh5qa+ffX/qAZe/37l86EZEg3XXXXTQ1NfXpd6RM0GeGQjywdANPrqwLuhQRkbgp6HshFDLOH1fMy7W70apZIhKPTZs2MX78eL75zW8yefJkrrvuOl544QXOP/98KioqWLFiBY2NjXz961/n7LPP5qyzzuKZZ545cu2sWbOYOnUqU6dO5dVXXwVgyZIlzJ49my984QuMHz+e6667rttM+pd/+Re2bdvGhRdeyIUXXgjA4sWLOffcc5k6dSpXXXUVBw4cOPGGunvS/Zk2bZofj9+8scXL/u5Zr/mo4biuF5GTZ+3atUGX4Bs3bvRwOOyrV6/29vZ2nzp1qn/ta1/zjo4O//3vf+9XXHGF/+hHP/JHH33U3d337NnjFRUVfuDAAW9sbPTm5mZ3d1+/fr0fzq0XX3zRBw4c6Fu3bvX29nafMWOGL1u2rNsaysrKvL6+3t3d6+vrfdasWX7gwAF3d7/99tv9tttu+9Q1Xf3dAdXeTabG9TC2v5g5rhiAl9/fxfjhAwOuRkTiddsf1rB2276EfubEEQP58f+Y1ON5Y8aM4fTTTwdg0qRJXHTRRZgZp59+Ops2baKuro4FCxZw5513ApHRQlu2bGHEiBHcfPPNrFq1inA4zPr164985vTp0yktLQVgypQpbNq0iZkzZ/ZYy/Lly1m7di3nn38+AC0tLZx77rm9bntnKRX0IwpzGVuSz7L3d/HNWWODLkdE+oHs7Owjv4dCoSPboVCItrY2wuEwTz75JKeddtpR1/3kJz9h2LBhvP3223R0dJCTk9PlZ4bDYdra2uKqxd2ZM2cOjz/++Ik06VNSKugBZo0r5rfVdRxqayc7Ixx0OSISh3juvIMyd+5c7rnnHu655x7MjLfeeouzzjqLhoYGSktLCYVCPPLII7S3tx/X5xcUFLB//36Ki4uZMWMGN910E7W1tYwbN46mpibq6uqorKw8oTakzMPYw2ZWlNDc2s6bm/cGXYqIpIBbb72V1tZWzjjjDCZPnsytt94KwI033sgjjzzCjBkzWL9+Pfn5+cf1+fPnz+eSSy7hwgsvpKSkhIcffphrrrmGM844gxkzZvDee++dcBvMk3CESlVVlR/vfPT7DrZy1k+f54YLxvKDueMTXJmIJEpNTQ0TJkwIuox+qau/OzNb6d0s7JRyd/QDczKZMqqQl9+Pe50TEZGUlnJBD5HRN6s/bKChqTXoUkREALjyyiuZMmXKUX8WLVp0Ur475R7GAsyqKObuP73Pqx/s4pLTTwm6HBERnn766cC+OyXv6M8cVciA7AyW1ar7RkQkJYM+Mxxixtgh6qcXSXLJOBgk2R3P31lKBj1E+um3fNzElt19O1mQiByfnJwcdu/W3FS94dGFR2JfzopHSvbRQ2Q8PcDLtbu4tkhr0Iokm9LSUurq6nq9iEa6O7yUYG+kbNCfWpLP8IE5vFxbz7XnKOhFkk1mZmavlsOT45eyXTdmxsyKYl6p3U17h/6voYikr5QNeogMs2xobuXdDxuCLkVEJDApHfTnnRqdtljDLEUkjaV00JcUZDN+eIGGWYpIWkvpoIdI983KzXtobjm+KURFRPq7lA/6mRUltLR38PrG3UGXIiISiJQP+unlQ8gKh3hF/fQikqZSPuhzs8JMKxvMMvXTi0iaSvmgB5hZUcx72/dTv/9Q0KWIiJx0cQW9mc0zs3VmVmtmt3RzzmwzW2Vma8zspU7Hwmb2lpk9m4iie2tWRWSYpbpvRCQd9Rj0ZhYG7gUuASYC15jZxE7nFAL3AZe7+yTgqk4f812gJhEFH49JIwZRmJep8fQikpbiuaOfDtS6+wZ3bwGeAK7odM61wFPuvgXA3XcePmBmpcClwIOJKbn3wiHjvFOLePn9XZopT0TSTjxBPxLYGrNdF90XqxIYbGZLzGylmV0fc+wu4IdAx7G+xMzmm1m1mVX3xWx2M8eVsH3fQT6oP5DwzxYRSWbxBL11sa/zbXEGMI3Inftc4FYzqzSzy4Cd7r6ypy9x9wfcvcrdq0pKSuIoq3cO99Nr9I2IpJt4gr4OGBWzXQps6+Kche7e6O67gKXAmcD5wOVmtolIl89nzew/T7jq4zBqSB5lRXl6ICsiaSeeoH8DqDCzMWaWBVwNLOh0zjPALDPLMLM84Bygxt1/5O6l7l4eve7P7v7lBNbfK+ePK2b5ho9pbT9mL5KISErpMejdvQ24GVhEZOTMb919jZndYGY3RM+pARYCq4EVwIPu/m7flX18Zo0r5sChNlZt3Rt0KSIiJ01cK0y5+3PAc5323d9p+w7gjmN8xhJgSa8rTKDzTi0mZJF++rPLhwRZiojISZMWb8YeNigvk9NLC9VPLyJpJa2CHmDmuCJWbd3LvoOtQZciInJSpGHQl9De4Sz/QNMWi0h6SLugn1pWSG5mWNMhiEjaSLugz84Ic87YIQp6EUkbaRf0ADPHFbOhvpFte5uDLkVEpM+lZ9BHp0PQouEikg7SMuhPG1ZASUE2y9R9IyJpIC2D3syYOa6YV2t30dGhaYtFJLWlZdBDZN6b3Y0t1GzfF3QpIiJ9Km2DfuY49dOLSHpI26AfPiiHiqEDNMxSRFJe2gY9REbfrNj4MU0tbUGXIiLSZ9I66OdMGMahtg6WrtddvYikrrQO+uljhjAoN5PFa7YHXYqISJ9J66DPCIe4aMJQXqjZoVWnRCRlpXXQA8ydNJx9B9tYsfHjoEsREekTaR/0n6koISczxCJ134hIikr7oM/NCvOZihIWr9mht2RFJCWlfdBDpPtm+76DvPNhQ9CliIgknIIeuGjCUMIhU/eNiKQkBT1QmJfFOWOGKOhFJCUp6KPmThrOB/WN1O48EHQpIiIJpaCPmjNxGACL1+quXkRSi4I+akRhLmeUDmLRmh1BlyIiklBxBb2ZzTOzdWZWa2a3dHPObDNbZWZrzOyl6L4cM1thZm9H99+WyOITbe6k4by9dS/bGw4GXYqISML0GPRmFgbuBS4BJgLXmNnETucUAvcBl7v7JOCq6KFDwGfd/UxgCjDPzGYkrPoEmzsp0n3zvLpvRCSFxHNHPx2odfcN7t4CPAFc0emca4Gn3H0LgLvvjP50dz/8dDMz+idp30oaN7SAsSX56r4RkZQST9CPBLbGbNdF98WqBAab2RIzW2lm1x8+YGZhM1sF7ASed/fXT7DmPnXxxOEs37CbhqbWoEsREUmIeILeutjX+a48A5gGXArMBW41s0oAd2939ylAKTDdzCZ3+SVm882s2syq6+vr460/4eZOGkZbh/PndbqrF5HUEE/Q1wGjYrZLgW1dnLPQ3RvdfRewFDgz9gR33wssAeZ19SXu/oC7V7l7VUlJSXzV94EzSwsZNjCbRe8q6EUkNcQT9G8AFWY2xsyygKuBBZ3OeQaYZWYZZpYHnAPUmFlJ9EEtZpYLfA54L2HV94FQyJgzcRgvra/nYGt70OWIiJywHoPe3duAm4FFQA3wW3dfY2Y3mNkN0XNqgIXAamAF8KC7vwucArxoZquJ/IPxvLs/2zdNSZy5k4bT3NrOsve1xKCI9H8Z8Zzk7s8Bz3Xad3+n7TuAOzrtWw2cdYI1nnQzxhZRkJPBojXbj7wxKyLSX+nN2C5khkNcNH4of6rZQZuWGBSRfk5B3425k4azp6mVNzbtCboUEZEToqDvxgWnlZCdoSUGRaT/U9B3Iy8rg1kVxTy/dgfuSfsyr4hIjxT0x3DxpOF8uLeZNdv2BV2KiMhxU9Afw+cmDCNkqPtGRPo1Bf0xDMnPYvqYISzWJGci0o8p6Htw8cThrNuxn027GoMuRUTkuCjoe3BxdI56dd+ISH+loO9B6eA8Jo8cyOK16r4Rkf5JQR+HiycO580te9i5X0sMikj/o6CPw9xJw3GHF9buDLoUEZFeU9DHoXLYAMqL8tRPLyL9koI+DmbGxZOG8+oHu9h3UEsMikj/oqCP09xJw2htd5asC26ZQxGR46Ggj9NZowZTPCBb3Tci0u8o6ON0eInBJe/t1BKDItKvKOh7Ye6kYTS2tLN0vbpvRKT/UND3wnmnFjOyMJd/fbFWUxeLSL+hoO+FrIwQ35tTyeq6Bp57R331ItI/KOh76cqzRlI5bAB3Ll5Hq9aTFZF+QEHfS+GQ8YO549m4q5HfVdcFXY6ISI8U9MfhcxOGUlU2mLteWE9zi0bgiEhyU9AfBzPj7y4Zz879h/jlqxuDLkdE5JgU9Mfp7PIhXDR+KP+25AP2NrUEXY6ISLfiCnozm2dm68ys1sxu6eac2Wa2yszWmNlL0X2jzOxFM6uJ7v9uIosP2g/mncaBQ23820sfBF2KiEi3egx6MwsD9wKXABOBa8xsYqdzCoH7gMvdfRJwVfRQG/B9d58AzABu6nxtfzZ++ECunDKSh1/ZxEcNzUGXIyLSpXju6KcDte6+wd1bgCeAKzqdcy3wlLtvAXD3ndGfH7n7m9Hf9wM1wMhEFZ8MvjenEne4+4X3gy5FRKRL8QT9SGBrzHYdnw7rSmCwmS0xs5Vmdn3nDzGzcuAs4PXjrDUpjRqSx3UzRvPb6q3U7jwQdDkiIp8ST9BbF/s6v/+fAUwDLgXmAreaWeWRDzAbADwJ/I277+vyS8zmm1m1mVXX1/evuWRuvnAcuZlh7ly0LuhSREQ+JZ6grwNGxWyXAtu6OGehuze6+y5gKXAmgJllEgn5x9z9qe6+xN0fcPcqd68qKSnpTRsCVzQgm299ZiwL12xn1da9QZcjInKUeIL+DaDCzMaYWRZwNbCg0znPALPMLMPM8oBzgBozM+AXQI27/1MiC08235w1lqL8LP7hv9/ThGciklR6DHp3bwNuBhYReZj6W3dfY2Y3mNkN0XNqgIXAamAF8KC7vwucD3wF+Gx06OUqM/t8H7UlUAOyM/jOZ8fx2obdLH1/V9DliIgcYcl491lVVeXV1dVBl9FrLW0dfPbnSxiYk8mz35lJKNTV4w0RkcQzs5XuXtXVMb0Zm0BZGSG+f3Elaz/ax7PvfBR0OSIigII+4a44cyTjhxfw88XraGnTNMYiEjwFfYKFQsbfzRvP5t1N/OaNLUGXIyKioO8Ls08rYfqYIdz9p1oaD7UFXY6IpDkFfR8wi9zV7zpwiF++ommMRSRYCvo+Mq1sMHMmDuPfX9rAnkZNYywiwVHQ96Efzj2NxpY27vlzbdCliEgaU9D3oYphBXzp7NH88tWNvFKrl6hEJBgK+j72fy6dwKklA/juE2+xc9/BoMsRkTSkoO9j+dkZ/Nt1U2k81M53Hn+LtnaNrReRk0tBfxJUDCvgZ1dO5vWNH/PPL6wPuhwRSTMK+pPkL6eWcvXZo7j3xQ94cd3OoMsRkTSioD+JfnL5JMYPL+B7v1nFtr1aY1ZETg4F/UmUkxnmvuum0tbu3PTrNzUXjoicFAr6k2xsyQBu/6vTeWvLXv5x4XtBlyMiaUBBH4DLzhjBV88t48GXN7J4zfagyxGRFKegD8j/vnQCZ5QO4vu/e5stu5uCLkdEUpiCPiDZGWHuvXYqBtz06zc51NYedEkikqIU9AEaNSSPn39xCu982MDP/lgTdDkikqIU9AGbM3EY8z8zll+9tpk/vL0t6HJEJAUp6JPAD+aexrSywdzy5Go21B8IuhwRSTEK+iSQGQ5xzzVnkZUR4sbH3uRgq/rrRSRxFPRJYkRhLv/8pSm8t30/P35mTdDliEgKUdAnkdmnDeXmC8fxm+qt/PC/3qapRevNisiJywi6ADna9+ZU4jj3LfmA6s17uOeas5g0YlDQZYlIP6Y7+iQTDhk/mDuex75xDgcOtnHlva/yy1c24u5BlyYi/VRcQW9m88xsnZnVmtkt3Zwz28xWmdkaM3spZv9DZrbTzN5NVNHp4Lxxxfz3d2cxs6KY2/6wlm/9qpqPtci4iByHHoPezMLAvcAlwETgGjOb2OmcQuA+4HJ3nwRcFXP4YWBegupNK0UDsvnFV6v4v5dNZOn6XVxy91Je+2B30GWJSD8Tzx39dKDW3Te4ewvwBHBFp3OuBZ5y9y0A7n5kZQ13Xwp8nKB6046Z8fWZY3jqxvPIz8rg2geX8/PF67QkoYjELZ6gHwlsjdmui+6LVQkMNrMlZrbSzK5PVIESMXnkIP7wnZn81dRS7vlzLV96YDl1ezQZmoj0LJ6gty72dX4ymAFMAy4F5gK3mlllbwoxs/lmVm1m1fX19b25NG3kZ2dw51VncvfVU1i3fT+fv3sZ//3OR0GXJSJJLp6grwNGxWyXAp0nZakDFrp7o7vvApYCZ/amEHd/wN2r3L2qpKSkN5emnSumjOSP/2smY4rz+fZjb/Kjp96huUVv04pI1+IJ+jeACjMbY2ZZwNXAgk7nPAPMMrMMM8sDzgE0HWMfKivK53c3nMdfXzCWx1dsoer/Pc+Nj63k6bfq2Nuk0Tki8okeX5hy9zYzuxlYBISBh9x9jZndED1+v7vXmNlCYDXQATzo7u8CmNnjwGyg2MzqgB+7+y/6pjnpJSsjxI8umcCcCcN48s0PeaFmB8+9s51wyKgqG8ycicOYM3EYZUX5QZcqIgGyZHwRp6qqyqurq4Muo9/p6HBWf9jAC2t38ELNDt7bvh+AiqED+Fw09KeUFhIKdfXYRUT6MzNb6e5VXR5T0KeuLbubeKEmEvqvb/yY9g6neEA2F40fytzJw7igcihhhb5ISlDQCw1NrSxZv5Pn1+7gpXX17D/UxqghufzP88bwxapSCnIygy5RRE6Agl6O0tLWwQs1O3jo5Y1Ub97DgOwMvlg1iq+dX86oIXlBlycix0FBL916e+tefvnKRp5d/REd7syZOIxvzBzL2eWDMVO3jkh/oaCXHm1vOMijyzfx2Otb2NvUyuSRA/nGzDFcevoIsjI0yalIslPQS9yaW9p5+q0PeeiVjdTuPEBJQTbXzyjjuhllDMnPCro8EemGgl56zd1Z+v4ufvHyRpauryc7I8TlZ47g+nPLOb1UC6GIJJtjBb1WmJIumRkXVJZwQWUJ7+/Yzy9f3cTv3/qQ362s48xRhXxlRhmXnXEKOZnhoEsVkR7ojl7itu9gK0+trOPR5Zv5oL6RwXmZfLFqFNedU8boIo3WEQmSum4kodyd1z7YzaPLN7N47Q463JldWcJXzi3TS1giAVHQS5/Z3nCQX6/YwuMrtlC//xCjhuRy3TllfLFqlB7eipxECnrpc63tHSxes4NfvbaJ1zd+TFZGiMvOOIVvzBzDpBF6eCvS1xT0clKt37GfR1/bzJNv1tHU0s55pxbxrVljuaCyRBOqifQRBb0EoqGplV+v2MLDr25kx75DnFqSzzdmjuUvp47UaB2RBFPQS6Ba2jp47p2P+I9lG1izbR9D8rP48owyvjKjjJKC7KDLE0kJCnpJCu7O8g0f84uXN/BCzU6yMkJcOWUk35g1hsphBUGXJ9Kv6YUpSQpmxrmnFnHuqUV8UH+Ah17eyH+trOM31Vu5oLKEb80ay8yK4qDLFEk5uqOXQH3c2MJjyzfzyGub2XXgEL/6+nQ+U6nF4UV661h39JqWUAI1JD+L71xUwZIfzAZg1da9gdYjkooU9JIUBmRncMqgHDbtagy6FJGUo6CXpFFelM+m3Qp6kURT0EvSKC/OY9PupqDLEEk5CnpJGuVF+Xzc2EJDc2vQpYikFAW9JI3y4nwANqv7RiShFPSSNMqLIkG/UQ9kRRIqrqA3s3lmts7Mas3slm7OmW1mq8xsjZm91JtrRQDKoouXbNqlfnqRROrxzVgzCwP3AnOAOuANM1vg7mtjzikE7gPmufsWMxsa77Uih+VkhhkxKEcjb0QSLJ47+ulArbtvcPcW4Angik7nXAs85e5bANx9Zy+uFTmivFhDLEUSLZ6gHwlsjdmui+6LVQkMNrMlZrbSzK7vxbUiR5QV5eulKZEEi2dSs65Wiug8QU4GMA24CMgFXjOz5XFeG/kSs/nAfIDRo0fHUZakojHFeexpaqWhqZVBeZlBlyOSEuK5o68DRsVslwLbujhnobs3uvsuYClwZpzXAuDuD7h7lbtXlZRoUqt0dXjkjbpvRBInnqB/A6gwszFmlgVcDSzodM4zwCwzyzCzPOAcoCbOa0WOODyWXkEvkjg9dt24e5uZ3QwsAsLAQ+6+xsxuiB6/391rzGwhsBroAB5093cBurq2j9oiKWD0kDzMNJZeJJHiWnjE3Z8Dnuu07/5O23cAd8RzrUh3IkMsc/VAViSB9GasJB1NbiaSWAp6STplmq5YJKEU9JJ0xhTls7eplb1NLUGXIpISFPSSdD4ZeaPuG5FEUNBL0ik/MrmZum9EEkFBL0lnlIZYiiSUgl6SzpEhlnogK5IQCnpJShpiKZI4CnpJSuWaxVIkYRT0kpTGFOfT0NzKnkYNsRQ5UQp6SUqaxVIkcRT0kpTKi6NDLBX0IidMQS9JadSQPEIGG7VQuMgJU9BLUsrOCDOiULNYiiSCgl6SVnlRPpvVdSNywuKaj14kCOXFeSxYtQ13x6yr5YdFkldrewcNza00NLeyt6mVfc2t7G1uoaGplb3R/Q1N0ePR7ZzMEM9+Z1bCa1HQS9IqL8pn38E29jS1MiQ/K+hyJA11dDj7D7WxLyawGw4HdmxQxwT2vubIzKuNLe3H/OyC7AwG5mZSmJfJoNxMKocNYGhBTp+0Q0EvSWtMzPqxCno5EQdb2z8J46aWTqF89F11Q6fjHd7952ZlhCjMjQR1YV4mIwtzmHjKwCPbh38OzM2MOS+LgTkZZIRPXs+5gl6SVtnhsfS7Gpk6enDA1UjQ2to72Hewrfugjv5siN5txwZ4S1tHt58bMo4K4kF5WYwuyj8qwDsH9eH9OZnhk/g3cPwU9JK0RkeHWGrkTf/i7hxq6+BgazvNre00t0R+Hmxtp7mlI7KvtZ2D0f2Hz4k9v7Gl7agA39fcyv5Dbcf83vyscORuORrKp5YMOHJXHdtFUpibdVSAF2RnEAql9jMgBb0krayMECMH52pyswRqa48N2o5ug/aTYP4kjI/e7ug+qFvb8WN0d3QnOyNEblaY3MwwedHQHj4wh9OGFTDoSEhnMigvEtax4T0wJ5OsDA0i7I6CXpJaeZqtH+vuNEf7k/c2tdLY0tZ98B4VtN3cQXfabm3vfQKHQ0ZeZpicaAjnHvk9RPGALHKzwuRkfnLsqO3oNTkxv0f2h446JzsjTDjF76qDpKCXpFZelM/vV33Y74ZYujsHDh3d/bC3ueWoB4J7Dw+ziznWU39yLDM+Cd5OQTowN5NhA7OPCtLcTwV116Hb+ZzMk/jQUPqGgl6SWnlxPvsPtvFxYwtFA7JP+vd3dDj7D7axp6mFvbEjNpo+Ce/D46L3Rs85vN1+jOEaeVnhaDdEVqf+5CwK8z558FeQk9llEOdkhsnOCPWrf/wkOAp6SWpjjkxu1nRCQd8WfXll75HRGC0xYR0ZUre3+dPbDc2tx+xvLsjOiPQZR/uNTynMpTDadzz4yOiMmPCO9ilnZ/SP0RqSGuIKejObB9wNhIEH3f32TsdnA88AG6O7nnL3n0aPfRf4FmDAf7j7XYkoXNJD7BDLaWWDOdTWfvQd9VF30V13hzQ0HXvEhhkxD/oid9jlRXlHbRfmZX7qjntgbqa6NaRf6DHozSwM3AvMAeqAN8xsgbuv7XTqMne/rNO1k4mE/HSgBVhoZn909/cTUr2kvFGDI0Msf7xgDbc+8y5Nx3jbMByymFEZmQwtyKFyaGTExuDoXfWRO+zcT+7CC3JSf3idpLd47uinA7XuvgHAzJ4ArgA6B31XJgDL3b0peu1LwJXAPx5fuZJusjJC/O3c09hY3xi9q/5kDHRhbtZRbx8OyM5Qn7VIF+IJ+pHA1pjtOuCcLs4718zeBrYBf+vua4B3gZ+ZWRHQDHweqD6xkiXd3Dh7XNAliPRr8QR9V7dInR9PvQmUufsBM/s88Hugwt1rzOwfgOeBA8DbQJedpWY2H5gPMHr06PiqFxGRHsXzJKkOGBWzXUrkrv0Id9/n7geivz8HZJpZcXT7F+4+1d0/A3wMdNk/7+4PuHuVu1eVlJQcR1NERKQr8QT9G0CFmY0xsyzgamBB7AlmNtyinaNmNj36ubuj20OjP0cDfwk8nrjyRUSkJz123bh7m5ndDCwiMrzyIXdfY2Y3RI/fD3wB+LaZtRHpi7/a/cjo4yejffStwE3uvqcvGiIiIl0zP57Zh/pYVVWVV1frma2ISLzMbKW7V3V1TG97iIikOAW9iEiKU9CLiKS4pOyjN7N6YHMPpxUDu05COclG7U4vand6OZF2l7l7l2PTkzLo42Fm1d09eEhland6UbvTS1+1W103IiIpTkEvIpLi+nPQPxB0AQFRu9OL2p1e+qTd/baPXkRE4tOf7+hFRCQO/S7ozWyema0zs1ozuyXoevqKmT1kZjvN7N2YfUPM7Hkzez/6c3CQNfYFMxtlZi+aWY2ZrYkuRZnybTezHDNbYWZvR9t9W3R/Srf7MDMLm9lbZvZsdDtd2r3JzN4xs1VmVh3dl/C296ugj1nW8BJgInCNmU0Mtqo+8zAwr9O+W4A/uXsF8KfodqppA77v7hOAGcBN0f+NU73th4DPuvuZwBRgnpnNIPXbfdh3gZqY7XRpN8CF7j4lZlhlwtver4KemGUN3b0FOLysYcpx96VE5u+PdQXwSPT3R4C/OJk1nQzu/pG7vxn9fT+R//hHkuJt94gD0c3M6B8nxdsNYGalwKXAgzG7U77dx5Dwtve3oO9qWcORAdUShGHu/hFEAhEYGnA9fcrMyoGzgNdJg7ZHuy9WATuB5909LdoN3AX8EOiI2ZcO7YbIP+aLzWxldJU96IO2x7OUYDKJZ1lDSQFmNgB4Evgbd9+XDot+u3s7MMXMCoGnzWxywCX1OTO7DNjp7ivNbHbA5QThfHffFl2g6Xkze68vvqS/3dH3uKxhitthZqcARH/uDLiePmFmmURC/jF3fyq6Oy3aDuDue4ElRJ7RpHq7zwcuN7NNRLpiP2tm/0nqtxsAd98W/bkTeJpI93TC297fgr7HZQ1T3ALgq9Hfvwo8E2AtfSK6JOUvgBp3/6eYQynddjMrid7JY2a5wOeA90jxdrv7j9y91N3Lifz3/Gd3/zIp3m4AM8s3s4LDvwMXA+/SB23vdy9MmdnnifTpHV7W8GfBVtQ3zOxxYDaR2ex2AD8Gfg/8FhgNbAGucvfOD2z7NTObCSwD3uGTPtv/TaSfPmXbbmZnEHnwFiZyA/Zbd/9pdBnOlG13rGjXzd+6+2Xp0G4zG0vkLh4i3ei/dvef9UXb+13Qi4hI7/S3rhsREeklBb2ISIpT0IuIpDgFvYhIilPQi4ikOAW9pDUzKzSzG4OuQ6QvKegl3RUCnwr66EypIilBQS/p7nbg1Oh84G9E58L/NfBOdJKxO6L7V5vZXx++yMx+ELP/8Nzx+Wb2x+ic8u+a2ZeCapRIrP42qZlIot0CTHb3KdE3M/8Y3d4YnU2wwd3PNrNs4BUzWwxURP9MJzLR3gIz+wxQAmxz90sBzGzQyW+OyKcp6EWOtsLdN0Z/vxg4w8y+EN0eRCTgL47+eSu6f0B0/zLgTjP7B+BZd1928soW6Z6CXuRojTG/G/Add18Ue4KZzQX+3t3/vfPFZjYN+Dzw92a22N1/2qfVisRBffSS7vYDBd0cWwR8OzptMmZWGZ1lcBHw9eic+ZjZSDMbamYjgCZ3/0/gTmBq35cv0jPd0Utac/fdZvaKRRZhbyYyU+hhDwLlwJvR6ZPrgb9w98VmNgF4LbogygHgy8A44A4z6wBagW+fvJaIdE+zV4qIpDh13YiIpDgFvYhIilPQi4ikOAW9iEiKU9CLiKQ4Bb2ISIpT0IuIpDgFvYhIivv/09vWcq1Fo6sAAAAASUVORK5CYII=\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.869729</td>\n",
       "      <td>0.648184</td>\n",
       "      <td>0.006621</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>0.913657</td>\n",
       "      <td>0.634728</td>\n",
       "      <td>0.007819</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>0.936296</td>\n",
       "      <td>0.623262</td>\n",
       "      <td>0.004499</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>0.951303</td>\n",
       "      <td>0.616497</td>\n",
       "      <td>0.003070</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>0.961392</td>\n",
       "      <td>0.611727</td>\n",
       "      <td>0.003258</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>11</td>\n",
       "      <td>0.968767</td>\n",
       "      <td>0.609092</td>\n",
       "      <td>0.002965</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>13</td>\n",
       "      <td>0.974395</td>\n",
       "      <td>0.607987</td>\n",
       "      <td>0.004230</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>15</td>\n",
       "      <td>0.978791</td>\n",
       "      <td>0.606400</td>\n",
       "      <td>0.003416</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>17</td>\n",
       "      <td>0.982301</td>\n",
       "      <td>0.604899</td>\n",
       "      <td>0.003877</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19</td>\n",
       "      <td>0.984962</td>\n",
       "      <td>0.603284</td>\n",
       "      <td>0.003800</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  trees   mean_tr   mean_te    std_te\n",
       "0     1  0.869729  0.648184  0.006621\n",
       "1     3  0.913657  0.634728  0.007819\n",
       "2     5  0.936296  0.623262  0.004499\n",
       "3     7  0.951303  0.616497  0.003070\n",
       "4     9  0.961392  0.611727  0.003258\n",
       "5    11  0.968767  0.609092  0.002965\n",
       "6    13  0.974395  0.607987  0.004230\n",
       "7    15  0.978791  0.606400  0.003416\n",
       "8    17  0.982301  0.604899  0.003877\n",
       "9    19  0.984962  0.603284  0.003800"
      ]
     },
     "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.53%\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:  1910\n",
      "True negative:  734058\n",
      "False positive:  3103\n",
      "False negative:  4158\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/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAs5klEQVR4nO3de7hVVb3/8feH+x3ZoIggXpL0kOWN8NLRo1KAHUvr0RNqSR47mpnWSY8/rU6aRkc7ZqWmHVMSTROyi5YpEmZmIQpqGhiBV5Cb3EEE9uX7+2OOBYvt2muvDXvDZO/P63nms9Ycc4wxx94bvmusMcccUxGBmZnlT7ud3QAzMyvNAdrMLKccoM3McsoB2swspxygzcxyygHazCynHKBtu0jqKuk3klZL+vl21HOWpEebs207i6RjJc3Z2e2wXZ88D7ptkHQm8BXgIGAt8DwwLiKe3M56PwNcBBwTETXb2868kxTAkIiYt7PbYq2fe9BtgKSvAN8Hvg30BwYDtwCnNEP1+wD/aAvBuRKSOuzsNlgrEhHeWvEG9AbWAaeXydOZLIAvTNv3gc7p2PHAAuASYCmwCDgnHfsmsAmoTuc4F7gK+GlR3fsCAXRI+58FXiHrxb8KnFWU/mRRuWOAZ4DV6fWYomOPA9cAf071PAr0a+BnK7T/sqL2nwp8FPgHsAL4alH+4cA0YFXKezPQKR17Iv0sb6ef91NF9f8/YDFwdyEtlXlPOsfhaX8vYBlw/M7+t+Et/5t70K3f0UAX4Fdl8nwNOAo4FDiELEh9vej4nmSBfiBZEP6hpD4RcSVZr3xiRPSIiDvKNURSd+BG4KSI6EkWhJ8vka8KeCjl7QvcADwkqW9RtjOBc4A9gE7ApWVOvSfZ72Ag8A3gx8CngSOAY4FvSNo/5a0F/hPoR/a7GwF8ASAijkt5Dkk/78Si+qvIvk2cV3ziiHiZLHjfI6kb8BPgzoh4vEx7zQAPcbQFfYFlUX4I4izg6ohYGhFvkfWMP1N0vDodr46I35H1Hg/cxvbUAQdL6hoRiyJiVok8/wrMjYi7I6ImIn4G/B34WFGen0TEPyLiHWAS2YdLQ6rJxturgfvIgu8PImJtOv8s4AMAETEzIp5K530N+D/gXyr4ma6MiI2pPVuJiB8Dc4HpwACyD0SzRjlAt37LgX6NjI3uBbxetP96SttcR70Avx7o0dSGRMTbZMMCnwcWSXpI0kEVtKfQpoFF+4ub0J7lEVGb3hcC6JKi4+8Uykt6r6TfSlosaQ3ZN4R+ZeoGeCsiNjSS58fAwcBNEbGxkbxmgAN0WzAN2EA27tqQhWRfzwsGp7Rt8TbQrWh/z+KDETE5Ij5C1pP8O1ngaqw9hTa9uY1taopbydo1JCJ6AV8F1EiZslOhJPUgG9e/A7gqDeGYNcoBupWLiNVk464/lHSqpG6SOko6SdJ3UrafAV+XtLukfin/T7fxlM8Dx0kaLKk3cEXhgKT+kj6exqI3kg2V1Jao43fAeyWdKamDpE8BQ4HfbmObmqInsAZYl3r3F9Q7vgTY/12lyvsBMDMiPkc2tv6j7W6ltQkO0G1ARNxANgf668BbwHzgi8CvU5ZvATOAF4AXgWdT2racawowMdU1k62Dajuy2SALyWY2/AvpAly9OpYDJ6e8y8lmYJwcEcu2pU1NdCnZBci1ZL37ifWOXwVMkLRK0r81VpmkU4DRZMM6kP0dDpd0VrO12Fot36hiZpZT7kGbmeWUA7SZWU45QJuZ5ZQDtJlZTrXphV36VbWPfffuuLObYU3wjxe6NZ7JcmUtK5dFxO7bWn7UCd1j+YpSszHfbeYLGydHxOhtPVfetOkAve/eHXl68t47uxnWBKP2OnRnN8Ga6Pdxf/27Qptk+Ypanp48uKK87QfMbeyuz11Kmw7QZpZ/AdRRt7ObsVM4QJtZrgVBdVQ2xNHaOECbWe65B21mlkNBUNtG73h2gDaz3Ksrv2Bgq+V50GaWawHUEhVtjZF0oKTni7Y1kr4sqUrSFElz02ufojJXSJonaY6kUUXpR0h6MR27UZJSemdJE1P6dEn7FpUZm84xV9LYxtrrAG1muVdHVLQ1JiLmRMShEXEo2SPP1pM9Du5yYGpEDAGmpn0kDQXGAO8jW5XwFkntU3W3kj3ibEjaCvOvzwVWRsQBwPeA61JdVcCVwJFkj5W7sviDoBQHaDPLtQCqIyrammgE8HJEvE72hPsJKX0CWx5wcQpwX3qc2avAPGC4pAFAr4iYFtmSoHfVK1Oo635gROpdjwKmRMSKiFgJTGFLUC/JY9BmlmtR4fBF0k/SjKL92yLitgbyjiF7WAVA/4hYBBARiyTtkdIHAk8VlVmQ0qrT+/rphTLzU101klaTPRt0c3qJMiU5QJtZvgXUVt45XhYRwxrLJKkT8HGKnvjTUNbSLWowfVvLlOQhDjPLtexOwsq2JjgJeDYiCg8PXpKGLUivS1P6AqB4PYhBZE8EWpDe10/fqkx6WHNvsicINVRXgxygzSznRG2FWxOcwZbhDYAHgcKsirHAA0XpY9LMjP3ILgY+nYZD1ko6Ko0vn12vTKGu04DH0jj1ZGCkpD7p4uDIlNYgD3GYWa5lFwmbFHzLktQN+AhwflHytcAkSecCbwCnA0TELEmTgNlADXBhxOb7zi8A7gS6Ag+nDbKnt98taR5Zz3lMqmuFpGuAZ1K+qyNiRbm2OkCbWa5l86CbL0BHxHqyi3bFacvJZnWUyj8OGFcifQZwcIn0DaQAX+LYeGB8pW11gDaz3Ktrxh70rsQB2sxyrbl70LsSB2gzy7VA1LbR+QwO0GaWex7iMDPLoUBsivaNZ2yFHKDNLNeyG1U8xGFmlku+SGhmlkMRojbcgzYzy6U696DNzPInu0jYNkNV2/ypzWyX4YuEZmY5Vut50GZm+eM7Cc3McqzOszjMzPInWyzJAdrMLHcCUe1bvc3M8icC36hiZpZP8o0qZmZ5FLgHbWaWW75IaGaWQ4Ha7IL9bfNjycx2GQFUR4eKtkpI2k3S/ZL+LuklSUdLqpI0RdLc9NqnKP8VkuZJmiNpVFH6EZJeTMdulKSU3lnSxJQ+XdK+RWXGpnPMlTS2sbY6QJtZzonaCrcK/QB4JCIOAg4BXgIuB6ZGxBBgatpH0lBgDPA+YDRwi6TCnL9bgfOAIWkbndLPBVZGxAHA94DrUl1VwJXAkcBw4MriD4JSHKDNLNeC7E7CSrbGSOoFHAfcARARmyJiFXAKMCFlmwCcmt6fAtwXERsj4lVgHjBc0gCgV0RMi4gA7qpXplDX/cCI1LseBUyJiBURsRKYwpagXpIDtJnlXhN60P0kzSjazqtX1f7AW8BPJD0n6XZJ3YH+EbEIIL3ukfIPBOYXlV+Q0gam9/XTtyoTETXAaqBvmboa5IuEZpZrEWrKWhzLImJYmeMdgMOBiyJiuqQfkIYzGlBq3CTKpG9rmZLcgzazXMsuEravaKvAAmBBRExP+/eTBewladiC9Lq0KP/eReUHAQtT+qAS6VuVkdQB6A2sKFNXgxygzSznsmcSVrI1JiIWA/MlHZiSRgCzgQeBwqyKscAD6f2DwJg0M2M/souBT6dhkLWSjkrjy2fXK1Oo6zTgsTROPRkYKalPujg4MqU1yEMcZpZr2UXCZp0HfRFwj6ROwCvAOWSd1UmSzgXeAE4HiIhZkiaRBfEa4MKIqE31XADcCXQFHk4bZBcg75Y0j6znPCbVtULSNcAzKd/VEbGiXEMdoM0s95rzTsKIeB4oNU49ooH844BxJdJnAAeXSN9ACvAljo0HxlfaVgdoM8u1tnwnoQO0meWeHxprZpZDEVBd5wBtZpY72RCHA7SZWS41YZ2NVsUBOkfmz+vMtz+/7+b9xW904jP/tZi1K9szbXJvJNitXzWXfv8N+u5Zsznf0gUd+Y/jD+LTlyzm9AveAmDuC125/suD2bihHcNPXMMF17yJBI9OrOL2a/ai757VAHz8nLc46axsps/t1wxg+tReRJ04/Li1m8tY03XsXMd3fzmPjp2C9h2CPz20G3dfvyfHnryKz1yymL2HbOTijw5h7gvdNpf51BeXMPqMFdTWiVu/vhcz/9gLgHH3vELVHtW07xD8bXoPbv7qQOrq2s4fpgWm2e0yWixASwrghoi4JO1fCvSIiKsqLH8VsC4irm/COX8HnBkRqySti4geTW/5zrP3ARu59fdzAKithbMOfx8fOmkVPXrXMvayxQD8+vZ+/PR7e/Kl67YsA/CjqwbywRPXblXXjZcP4kvfmc8/HbGer396f2b8oefmPMd9fCVf/PabW+Wf9Uw3Zj3TnR9Nzc5/yalDeGFaDw45Zl2L/bytWfVGcdnp72HD+va07xDc8Ot5PPNYT177exeu/ty+XHzdgq3yDx6ygeNPWcV5JxxIVf9qrp34Cuf+c0/q6sS48/dh/br2QPDfP36dYz+2ij8+UHYRtFam7Q5xtORPvRH4pKR+LXiOrUTER9PKVLu85//UkwH7bKT/oGq696zbnL7hnXZb9Wr/8nBvBgzexD7v3bA5bfmSDqxf256hw9YjwYdPW8FfHuld9nwSbNrYjppNonqjqKkWfXavbvafq+0QG9Zntx536Bi07xhEwPx5XVjwcpd35T561Goef2A3qje1Y8n8zix8rRMHHrYeIAVnaN8BOnSKRlZvaJ3q0nMJG9tam5YM0DXAbcB/1j8gaR9JUyW9kF4HN1DHUEmPS3pF0sVF5X8taaakWcWrVUl6bUd+ILSkxx/YjeNPXbV5/yfX7slZRwzlsV/24ez/WgTAhvXtmHTLHnz6ksVblV2+uCP9BmwJrv32qmbZ4o6b9//8u934/IgDueY/9mXpm1n60GHrOeSYdZxx2MGccdjBHHH8GgYP2diCP2Hr165dcMuUOUx8YRbPPdGDOc91bzBvvwHVvLWw0+b9ZYs6bR6GAhh378tMfGEW76xrx59+u1tLNjt3slkc7SvaWpuW/t7wQ+AsSfW7bzcDd0XEB4B7gBsbKH8Q2RqqhcWtC1Hm3yPiCLK7gS6W1LfSBkk6r7AU4VvLaxsvsBNUbxJPPdqb4z62anPaOZcv5p6Zsznxkyt5cPzuANz1v3vyif94i67d67YqHyV6WIW+xVEfWc2E6bP50dQ5HHbsWq7/cvbZ+OarnZg/rzP3zJzFvc/O4q9/7smLTzUcUKxxdXXiCx85kLOOGMqBh65nnwPfaThzI+ucfe3M93DGYUPp2Ck49J/b1rBT4UaVSrbWpkUDdESsIVvI+uJ6h44G7k3v7wb+uYEqHkoLZS8jW12qf0q/WNJfgafIVoca0oQ23RYRwyJi2O598/mJ+8xjPTng/evps3vNu46d8ImVPPm77PPu7891445v7cXZw4fyq9t3576b+vPA+H70G1DNskVbeszLFnbc3BvrVVVLp87Z//yTzlq++SLVXx7uzUGHr6dr9zq6dq9j2AlreGmmA3RzeHtNe/46rQcfPGFtg3mWLezI7ntt2rzfb8Amli/puFWe6o3tmPZoL44etbrF2ppXHuJoOd8newRMuf/tDY2qFX/HrgU6SDoe+DBwdEQcAjwHvHtQbxf2+K/7bDW88eYrW776PjW5N3sfkP1abvj1PO56ejZ3PT2bT3zuLcZctIRT/n0ZffvX0K1HHS/N7EYE/P7+qs3/qZcv2XJd+KlHezN4SDZ2vfvAal6Y1oPaGqiphhef6rH5mDVd76oauvfKvqF16lLH4ceuY/68hv+ZPvVob44/ZRUdO9XRf++NDNxvE3Oe60aXbrVU7ZF9uLZrHwwfsaZsPa1RYRZHW+xBt/g0u7SC0ySyIF1YJOQvZCs83Q2cBTzZhCp7kz3va72kg4CjmrO9O9uG9eLZP/XkS9/Z8uCFO769Fwte7ky7drDHwE3vmgFQykXXzuf6Lw9m04Z2DDthzeYZHA/csTvTHu1F+w7Qc7caLvneGwAce/Iq/vrnHpx/4kFIMOyENRw1ck3L/JBtQFX/ai79wRu0awft2sETv+nN9N/34pjRq/nCt96kd98arrn7VV6e1YWvnfkeXv9HF574zW7c9vgcamu1eSpdl251XHXnq9l0vfbB83/uwW/vqnhEr9Voq7M4FKUGLJuj4qJpbpL6A68C34mIq9JTbscD/cgeP3NORLxRr/xVFE2zk/Q34GRgEfBrskfFzAF2B66KiMclvQYMi4hllUyzG3ZIl3h68t7lsljOjNrr0J3dBGui38f9Mxt5yklZfQ7aI04cf1pFeX/5oVu361x502I96OLgGBFLgG5F+68BJzZS/qp6+8XL+p3UQJl9S53fzHZtrXH4ohK+k9DMcs13EpqZ5ZgDtJlZDnnBfjOzHGuNc5wr4QBtZrkWATVesN/MLJ/a6hBH2/xYMrNdRnOvxZEWVXtR0vOSZqS0KklTJM1Nr32K8l8haZ6kOZJGFaUfkeqZJ+lGKVtnUlJnSRNT+vR030ehzNh0jrmSxjbWVgdoM8u9CFW0NcEJEXFo0U0tlwNTI2IIMDXtI2ko2V3P7wNGA7dIKizicytwHtlaQEPSccjuml4ZEQcA3wOuS3VVAVcCR7JlAbiyC3s7QJtZ7u2AxZJOASak9xOAU4vS70uLtr0KzAOGSxoA9IqIaZHdjn1XvTKFuu4HRqTe9ShgSkSsiIiVwBS2BPWSHKDNLNcimrRYUr/CcsJpO69UlcCjaU35wvH+EbEoO18sAvZI6QOB+UVlF6S0gel9/fStykREDbAa6Fumrgb5IqGZ5ZyorXwWx7IK1uL4UEQslLQHMEXS38ue/N2iTPq2linJPWgzy73mHIOOiIXpdSnwK7Lx4CVp2IL0ujRlX0C25nzBIGBhSh9UIn2rMpI6kK3AuaJMXQ1ygDazXGvO9aAldZfUs/AeGAn8DXgQKMyqGAs8kN4/CIxJMzP2I7sY+HQaBlkr6ag0vnx2vTKFuk4DHkvj1JOBkZL6pIuDI1NagzzEYWb5FqUf47aN+gO/SjPiOgD3RsQjkp4BJkk6F3gDOB0gImal9exnkz1n9cKIKDwr7wLgTqAr8HDaAO4A7pY0j6znPCbVtULSNcAzKd/VEbGiXGMdoM0s95rrVu+IeAU4pET6cmBEA2XGAeNKpM8ADi6RvoEU4EscG8+WB5c0ygHazHItmnaRsFVxgDaz3GuhBz/lngO0meVeE+8SbDUcoM0s1yIcoM3McqutrmbnAG1muecxaDOzHApEnWdxmJnlUxvtQDtAm1nO+SKhmVmOtdEutAO0meWee9D1SLqJMp9bEXFxi7TIzKxIAHV1DtD1zdhhrTAza0gA7kFvLSImFO9L6h4Rb7d8k8zMttZW50E3OrlQ0tGSZgMvpf1DJN3S4i0zMyuICrdWppLZ398nexrtcoCI+CtwXAu2ycysSGWPu2qNFxIrmsUREfPTEwgKahvKa2bW7Fph77gSlQTo+ZKOAUJSJ+Bi0nCHmVmLC4g2OoujkiGOzwMXAgOBN4FD076Z2Q6iCrfWpdEedEQsA87aAW0xMyutjQ5xVDKLY39Jv5H0lqSlkh6QtP+OaJyZGdDsszgktZf0nKTfpv0qSVMkzU2vfYryXiFpnqQ5kkYVpR8h6cV07EalC3WSOkuamNKnS9q3qMzYdI65ksY21s5KhjjuBSYBA4C9gJ8DP6vw92Bmtn0KN6pUslXuS2x9Le1yYGpEDAGmpn0kDQXGAO8DRgO3SGqfytwKnAcMSdvolH4usDIiDgC+B1yX6qoCrgSOBIYDVxZ/EJRSSYBWRNwdETVp+ylt9guHme0M2WOvGt8qIWkQ8K/A7UXJpwCFm/MmAKcWpd8XERsj4lVgHjBc0gCgV0RMi4gA7qpXplDX/cCI1LseBUyJiBURsRKYwpagXlK5tTiq0ts/SLocuI8sMH8KeKhcpWZmzap5Z3F8H7gM6FmU1j8iFgFExCJJe6T0gcBTRfkWpLTq9L5+eqHM/FRXjaTVQN/i9BJlSip3kXAmWUAu/GbOLzoWwDXlKjYzay6q/Dt7P0nF6wjdFhG3ba5HOhlYGhEzJR1fyalLpEWZ9G0tU1K5tTj2K1fQzGyHaNoFwGURMazM8Q8BH5f0UaAL0EvST4Elkgak3vMAYGnKvwDYu6j8IGBhSh9UIr24zAJJHYDewIqUfny9Mo+X+2EqetCXpIMl/ZukswtbJeXMzLZfhRcIK7hIGBFXRMSgiNiX7OLfYxHxaeBBoDCrYizwQHr/IDAmzczYj+xi4NNpOGStpKPS+PLZ9coU6jotnSOAycBISX3SxcGRKa1Bjc6DlnQlWdQfCvwOOAl4kmxQ3Mys5bX8tIRrgUmSzgXeAE4HiIhZkiYBs4Ea4MKIKCx1cQFwJ9AVeDhtAHcAd0uaR9ZzHpPqWiHpGuCZlO/qiFhRrlGV3Op9GnAI8FxEnCOpP1tf/TQza1l1zV9lRDxOGmKIiOXAiAbyjQPGlUifARxcIn0DKcCXODYeGF9pGysJ0O9ERJ2kGkm9yMZmfKOKme0YXrC/rBmSdgN+TDazYx3wdEs2ysysWBNmcbQqlazF8YX09keSHiGbnP1CyzbLzKyIA/TWJB1e7lhEPNsyTTIzMyjfg/5umWMBnNjMbdnh/vFCN0btdejOboaZNcJDHPVExAk7siFmZiUFzX2r9y6jokdemZntVO5Bm5nlk4c4zMzyqo0G6EqeqCJJn5b0jbQ/WNLwlm+amVnSzE9U2VVUsljSLcDRwBlpfy3wwxZrkZlZEUXlW2tTyRDHkRFxuKTnACJipaROLdwuM7MtPIujQdXpGVwBIGl3WmTpEjOz0lpj77gSlQxx3Aj8CthD0jiypUa/3aKtMjMr1kbHoCtZi+MeSTPJluITcGpEvNRIMTOz5tFKx5crUcmC/YOB9cBvitMi4o2WbJiZ2WYO0A16iC0PPOwC7AfMAd7Xgu0yM9tMbfSqVyVDHO8v3k+r3J3fQHYzM2smTb6TMCKelfTBlmiMmVlJHuIoTdJXinbbAYcDb7VYi8zMivkiYVk9i97XkI1J/6JlmmNmVkIbDdBl50GnG1R6RMQ30zYuIu5JT601M9sxmmketKQukp6W9FdJsyR9M6VXSZoiaW567VNU5gpJ8yTNkTSqKP0ISS+mYzdKUkrvLGliSp8uad+iMmPTOeZKGttYexsM0JI6REQt2ZCGmdlOIbJZHJVsFdgInBgRhwCHAqMlHQVcDkyNiCHA1LSPpKHAGLJZa6OBW1LHFeBW4DxgSNpGp/RzgZURcQDwPeC6VFcVcCVwJDAcuLL4g6CUcj3owpO7n5f0oKTPSPpkYavoV2Fmtr2acbGkyKxLux3TFsApwISUPgE4Nb0/BbgvIjZGxKvAPGC4pAFkD9CeFhEB3FWvTKGu+4ERqXc9CpgSESsiYiUwhS1BvaRKxqCrgOVkzyAszIcO4JcVlDUz236Vj0H3kzSjaP+2iLitOEPqAc8EDgB+GBHTJfWPiEUAEbFI0h4p+0DgqaLiC1JadXpfP71QZn6qq0bSaqBvcXqJMiWVC9B7pBkcf2NLYC5oo0P2ZrZTVB5xlkXEsLJVZUO3h0raDfiVpIPLZC+1jF79eFi/ldtSpqRyQxztgR5p61n0vrCZme0QLbEedESsAh4nG2ZYkoYtSK9LU7YFwN5FxQYBC1P6oBLpW5WR1AHoDawoU1eDyvWgF0XE1eUKm5ntEM30nT0tl1wdEaskdQU+THYR70FgLHBten0gFXkQuFfSDcBeZBcDn46IWklr0wXG6cDZwE1FZcYC04DTgMciIiRNBr5ddGFwJHBFufaWC9Btc4VsM8uXaNa1OAYAE9I4dDtgUkT8VtI0YJKkc4E3gNMBImKWpEnAbLL7QC5MQyQAFwB3Al2Bh9MGcAdwt6R5ZD3nMamuFZKuAZ5J+a6OiBXlGlsuQI+o/Gc2M2tBzdSDjogXgMNKpC+ngZgXEeOAcSXSZwDvGr9O94mc3kBd44Hxlba3wQDdWGQ3M9tRfKu3mVleOUCbmeVQK32cVSUcoM0s14SHOMzMcssB2swsrxygzcxyygHazCyH/EQVM7Mcc4A2M8unZrzVe5fiAG1muechDjOzPPKNKmZmOeYAbWaWP76T0Mwsx1TXNiO0A7SZ5ZvHoM3M8stDHGZmeeUAbWaWT+5Bm5nllQO0mVkONe9TvXcp7XZ2A8zMyinMg65ka7QuaW9Jf5D0kqRZkr6U0qskTZE0N732KSpzhaR5kuZIGlWUfoSkF9OxGyUppXeWNDGlT5e0b1GZsekccyWNbay9DtBmln8RlW2NqwEuiYh/Ao4CLpQ0FLgcmBoRQ4CpaZ90bAzwPmA0cIuk9qmuW4HzgCFpG53SzwVWRsQBwPeA61JdVcCVwJHAcODK4g+CUhygzSz3mqsHHRGLIuLZ9H4t8BIwEDgFmJCyTQBOTe9PAe6LiI0R8SowDxguaQDQKyKmRUQAd9UrU6jrfmBE6l2PAqZExIqIWAlMYUtQL8lj0LuQdu2Cmx75B8sXdeQbY/fn2JNX8ZlLFrP3kI1c/NEhzH2hGwD9B23ix3/8Owte6QzA32d258bLBwFw/KkrGXPRUiJgxZKOXHfRYNas8D+DlvCVG97gyA+vZdWyDpx/4oEA7D/0HS66dgFdu9exZEEnrrtwMOvXtadnnxr++7bXeO+h7zBlUh9++LVBm+s54P3rufT78+ncpY6nH+vFrf+9F9kX/zaiaTeq9JM0o2j/toi4rVTGNPRwGDAd6B8RiyAL4pL2SNkGAk8VFVuQ0qrT+/rphTLzU101klYDfYvTS5QpKbc9aEnr6u1/VtLNjZT5uKTCV5OrJF3akm3c0U793DLmz+2yef+1v3fh6s/ty4tPdX9X3kWvd+YLHzmQL3zkwM3BuV374IKrF3LZ6e/hgg8fyCsvdeHj5yzbYe1vax6dWMXXztpvq7QvXz+f8d8ewOdHHMifH+7FaRcsBWDTBjHhf/fkx1cPeFc9F1+7gB9cNohzPnQQA/fbyLAT1u6Q9ueJ6irbgGURMaxoayg49wB+AXw5ItaUO3WJtCiTvq1lSsptgN4WEfFgRFy7s9vREvoN2MTwEWt4+N6qzWnz53VhwctdypTamtLVli5d64Cge486li/u2PyNNQD+Nr0Ha1du/e1k0Hs2bv5Afe6Jnvzzv64GYOM77Zn1dA82bdz6v2TVHtV061nHSzO7A+L39/fhmNGrd0j786QJAbrxuqSOZMH5noj4ZUpekoYtSK9LU/oCYO+i4oOAhSl9UIn0rcpI6gD0BlaUqatBu2SAlvSxdHX0OUm/l9Q/pTfay95Vff6bC7n9WwOIusq+2u45eBM/fHQO//uLeRw8PPsyUlsjbrp8ED96bA73Pjebwe/dwOSfVTVSkzWn1+d04ehRWYft2JNXs/te1WXz992zmmWLtnyILlvYkX57li/T6gTNdpEwjQXfAbwUETcUHXoQKMyqGAs8UJQ+Js3M2I/sYuDTaThkraSjUp1n1ytTqOs04LE0Tj0ZGCmpT7o4ODKlNSjPAbqrpOcLG3B10bEngaMi4jDgPuCySiuVdJ6kGZJmVLOxeVvcQo788BpWLevAvBe7VZR/xdIOfPqD/8SFIw/k/67ai8tveYNuPWpp3yE4+ezlXDjyvZx52FBefakLn7poaeMVWrO54St787HPLuPmR/5B1x611Gwq/4GrEoejLY0/J811kRD4EPAZ4MSi+PJR4FrgI5LmAh9J+0TELGASMBt4BLgwImpTXRcAt5NdOHwZeDil3wH0lTQP+AppRkhErACuAZ5J29UprUF5vjr0TkQcWtiR9FlgWNodBExMX0U6Aa9WWmkak7oNoJeqdon7k4Z+8G2OGrmGD46YTafOQbeetVx20+t856J9Suav3tSO6k3ZZ++8F7ux8LVODNx/4+b/7Itezy4e/vHB3fjUFx2gd6T587rw1TPeA8DA/Tdy5Ihyw5+wbFFH+g3Y0mPut1c1yxfn+b9tC2mm/6kR8SQNX2Ed0UCZccC4EukzgINLpG8ATm+grvHA+Erbm+cedDk3ATdHxPuB84HKB2J3QT/5nwF8ethQxh45lP+5YB/++mSPBoMzQO+qGtq1y/5F7zl4IwP328jiNzqxbHFHBr93A72ragA4/Li1W110tJbXu28WbKXgzC8t4bd39y2bf8XSjqxf146DDn8bCD582kqmTe69A1qaH815o8quZlf9KO4NvJneN3o3Tmt1zOjVfOFbb9K7bw3X3P0qL8/qwtfOfA/vP2odZ//XYmprRG2duPHyQaxdlf2p77mhP9f/ah411WLpm524/st7N3IW21aX3/I6Hzh6Hb2ravjpjNnc/d3+dO1Wx8c+m82c+fPDvXn0vi3XACZMn033HnV06BQcPWoNXz1jf96Y24WbLh/Epd+fT6cudcz4Q0+eeaznzvqRdo6INrtgv6Kyu292OEnrIqJH0f5ngWER8UVJp5DdofMm2RzFD0bE8fXyXAWsi4jrGzpHL1XFkSr5rcbMmsnv4/6ZETGs8Zyl9dxtUBx23Jcqyvun31y2XefKm9z2oIuDc9q/E7gzvX+ALVdMG8pzVcu20Mx2lNY4fFGJ3AZoMzMgu0DYRoc4HKDNLP/aZnx2gDaz/PMQh5lZTrXVWRwO0GaWb01bza5VcYA2s1zLblRpmxHaAdrM8q+NPpPQAdrMcs89aDOzPPIYtJlZXrXdtTgcoM0s/zzEYWaWQ1H546xaGwdoM8s/96DNzHKqbcZnB2gzyz/Vtc0xDgdoM8u3wDeqmJnlkYg2e6PKrvrQWDNrSyIq2xohabykpZL+VpRWJWmKpLnptU/RsSskzZM0R9KoovQjJL2Yjt0oSSm9s6SJKX26pH2LyoxN55grqaJnqTpAm1n+NVOAJnsk3uh6aZcDUyNiCDA17SNpKDAGeF8qc4uk9qnMrcB5wJC0Feo8F1gZEQeQPTf1ulRXFXAlcCQwHLiy+IOgIQ7QZpZvhTHoSrbGqop4AlhRL/kUYEJ6PwE4tSj9vojYGBGvAvOA4ZIGAL0iYlpkT92+q16ZQl33AyNS73oUMCUiVkTESmAK7/6geBePQZtZ7jVhFkc/STOK9m+LiNsaKdM/IhYBRMQiSXuk9IHAU0X5FqS06vS+fnqhzPxUV42k1UDf4vQSZRrkAG1mOVfx8AXAsogY1kwnVunGNJi+rWUa5CEOM8u3oDnHoEtZkoYtSK9LU/oCYO+ifIOAhSl9UIn0rcpI6gD0JhtSaaiushygzSz/mmkMugEPAoVZFWOBB4rSx6SZGfuRXQx8Og2HrJV0VBpfPrtemUJdpwGPpXHqycBISX3SxcGRKa0sD3GYWe411zxoST8Djicbq15ANrPiWmCSpHOBN4DTASJilqRJwGygBrgwImpTVReQzQjpCjycNoA7gLslzSPrOY9Jda2QdA3wTMp3dUTUv1j5Lg7QZpZ/zRSgI+KMBg6NaCD/OGBcifQZwMEl0jeQAnyJY+OB8RU3FgdoM8u7CKhtm/d6O0CbWf610Vu9HaDNLP8coM3McigAP5PQzCyPAsJj0GZm+RP4IqGZWW55DNrMLKccoM3M8mi71tnYpTlAm1m+BeCHxpqZ5ZR70GZmeeRbvc3M8ikgPA/azCynfCehmVlOeQzazCyHIjyLw8wst9yDNjPLoyBqaxvP1go5QJtZvnm5UTOzHPM0OzOz/Akg3IM2M8uh8IL9Zma51VYvEira6PQVAElvAa/v7Ha0gH7Asp3dCGuS1vw32ycidt/WwpIeIfv9VGJZRIze1nPlTZsO0K2VpBkRMWxnt8Mq57+ZldJuZzfAzMxKc4A2M8spB+jW6bad3QBrMv/N7F08Bm1mllPuQZuZ5ZQDtJlZTjlA55CkkPTdov1LJV3VhPJXSbq0ief8naTd0vt1TSlr5dX/fUr6rKSbGynzcUmXp/dN/nta6+AAnU8bgU9KqnRy/naLiI9GxKoddT4rLyIejIhrd3Y7bOdygM6nGrKr+v9Z/4CkfSRNlfRCeh3cQB1DJT0u6RVJFxeV/7WkmZJmSTqvKP21HfmBYBlJH5M0XdJzkn4vqX9Kb7SXba2fA3R+/RA4S1Lveuk3A3dFxAeAe4AbGyh/EDAKGA5cKaljSv/3iDgCGAZcLKlv8zfd6ukq6fnCBlxddOxJ4KiIOAy4D7hsZzTQ8smLJeVURKyRdBdwMfBO0aGjgU+m93cD32mgiociYiOwUdJSoD+wgCwofyLl2RsYAixv7vbbVt6JiEMLO5I+S/YBCTAImChpANAJeHWHt85yyz3ofPs+cC7QvUyehiaybyx6Xwt0kHQ88GHg6Ig4BHgO6LLdrbTtcRNwc0S8Hzgf/z2siAN0jkXECmASWZAu+AswJr0/i+wrcqV6AysjYr2kg4CjmqWhtj16A2+m92N3ZkMsfxyg8++7bL3U4sXAOZJeAD4DfKkJdT1C1pN+AbgGeKrZWmnb6irg55L+ROtdbtS2kW/1NjPLKfegzcxyygHazCynHKDNzHLKAdrMLKccoM3McsoB2sqSVJtuUf6bpJ9L6rYddd0p6bT0/nZJQ8vkPV7SMdtwjpJrilSy1khTV/HzKnPW0hygrTHvRMShEXEwsAn4fPFBSe23pdKI+FxEzC6T5XigyQHarDVxgLam+BNwQOrd/kHSvcCLktpL+l9Jz6RV9s4HUOZmSbMlPQTsUagorbQ3LL0fLelZSX9NK/TtS/ZB8J+p936spN0l/SKd4xlJH0pl+0p6NK0G93+AGvshGlrRLx37bmrLVEm7p7T3SHoklflTugvTrMV5sSSriKQOwElkdyNCtkrewRHxagpyqyPig5I6A3+W9ChwGHAg8H6yxZpmA+Pr1bs78GPguFRXVUSskPQjYF1EXJ/y3Qt8LyKeTEusTgb+CbgSeDIirpb0r8BWAbcB/57O0RV4RtIvImI52Zonz0bEJZK+ker+ItnSr5+PiLmSjgRuAU7chl+jWZM4QFtjuqYlMiHrQd9BNvTwdEQUVl4bCXygML5Mtr7EEOA44GcRUQsslPRYifqPAp4o1JXWHynlw2RrXBf2e0nqmc7xyVT2IUkrK/iZGlrRrw6YmNJ/CvxSUo/08/686NydKziH2XZzgLbGbLVUJkAKVG8XJwEXRcTkevk+SsOr7RWXrWS9gXZkq/AVL71aaEvF6xXUW9FvvaTHaXgFuUjnXVX/d2C2I3gM2prDZOCCwkMBJL1XUnfgCWBMGqMeAJxQouw04F8k7ZfKVqX0tUDPonyPkg03kPIdmt4+QbaqH5JOAvo00tZyK/q1AwrfAs4kGzpZA7wq6fR0Dkk6pJFzmDULB2hrDreTjS8/K+lvwP+RfTv7FTAXeBG4Ffhj/YIR8RbZuPEvJf2VLUMMvwE+UbhISLaK37B0EXI2W2aTfBM4TtKzZEMtbzTS1nIr+r0NvE/STLIx5sKTT84Czk3tmwWcUsHvxGy7eTU7M7Occg/azCynHKDNzHLKAdrMLKccoM3McsoB2swspxygzcxyygHazCyn/j/Jq1URSVE1tAAAAABJRU5ErkJggg==\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): 31.47659854976928\n",
      "False Alarm Rate (FAR): 0.42093925207654775\n",
      "Positive Predictive Value (PPV): 0.3810093756233792\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": [
    "There is a minor improvement respect to the class weighted model but anyway the performance is very poor"
   ]
  },
  {
   "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": [
    "features = ['HRV_norm', 'HRV_norm_text', 'IR_108', 'IR_108_text', 'WV_062-IR_108', 'WV_062-IR_108_text']\n",
    "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/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfI0lEQVR4nO3df5zVVb3v8dcbVBBTNKkOgkoZmhaJOvmDfnmvHn+W6M2jeCqlTPrt0cRf6elwzRLSksysyLxmeckfqdHVxOqolSI5KAxImSCoII9QrDkSSgqf+8d3TWzn7Jm1Z/bM7L1n3s/HYx6zv2t91/qutTfO2/Xda88oIjAzM7OODar1AMzMzOqdw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmbWYyR9UdK1tR6HWU+TP2dpVh8krQTeBGwqKd4zIp6tss9PRMSvqhtd45E0DXhrRHyk1mOxxueVpVl9+WBEvK7kq9tB2RMkbVXL63dXo47b6pfD0qzOSRou6QeS1khaLelSSYNT3R6S/lPSOknPS7pR0o6p7kfAbsDPJa2XdJ6kQyWtatf/SkmHp8fTJN0q6ceS/guY3Nn1y4x1mqQfp8djJIWkj0l6RtJfJH1K0rsktUj6q6SrS9pOlvSApG9JapX0R0mHldTvImmOpBckLZN0Rrvrlo77U8AXgZPT3Bel8z4m6Q+SXpT0pKRPlvRxqKRVks6RtDbN92Ml9dtK+rqkp9L4fidp21R3sKQH05wWSTq0Gy+11TGHpVn9+yHwKvBWYD/gCOATqU7AZcAuwN7ArsA0gIj4KPA0W1arX6vwehOBW4EdgRsz16/EQcBY4GRgJnARcDjwduAkSe9vd+6TwAjgP4DbJL0+1c0GVqW5ngh8tTRM2437B8BXgZvS3PdN56wFPgDsAHwMuFLS/iV9/BMwHBgFnA58W9JOqe4K4ABgAvB64Dxgs6RRwJ3Apal8KvBTSW/ownNkdc5haVZf7kirk79KukPSm4CjgbMi4m8RsRa4EpgEEBHLIuKXEbExIp4DvgG8v+PuKzIvIu6IiM0UodLh9Sv05Yh4OSLuAf4GzI6ItRGxGvgtRQC3WQvMjIhXIuIm4HHgWEm7Au8Bzk99LQSuBT5abtwR8VK5gUTEnRGxPAr3A/cA7y055RXgknT9u4D1wF6SBgEfB/4tIlZHxKaIeDAiNgIfAe6KiLvStX8JNAPHdOE5sjrn+/pm9eX40s04kg4EtgbWSGorHgQ8k+rfCFxF8QN/+1T3lyrH8EzJ4907u36F/lzy+KUyx68rOV4dr911+BTFSnIX4IWIeLFdXVMH4y5L0tEUK9Y9KeYxDFhccsq6iHi15HhDGt8IYCiwvEy3uwP/IumDJWVbA/fmxmONw2FpVt+eATYCI9r9EG9zGRDAOyNinaTjgatL6ttvd/8bRUAAkN57bH+7sLRN7vo9bZQklQTmbsAc4Fng9ZK2LwnM3YDVJW3bz/U1x5KGAD8FTgV+FhGvSLqD4lZ2zvPAy8AewKJ2dc8AP4qIM/5bK+s3fBvWrI5FxBqKW4Vfl7SDpEFpU0/brdbtKW4V/jW9d3Zuuy7+DLyl5PhPwFBJx0raGrgYGFLF9XvaG4EzJW0t6V8o3oe9KyKeAR4ELpM0VNI7Kd5TvLGTvv4MjEm3UAG2oZjrc8CraZV5RCWDSrekrwO+kTYaDZZ0SArgHwMflHRkKh+aNguN7vr0rV45LM3q36kUP+iXUtxivRUYmer+N7A/0EqxyeS2dm0vAy5O74FOjYhW4DMU7/etplhprqJznV2/p82n2Az0PPAV4MSIWJfqTgHGUKwybwf+I70/2JFb0vd1kh5JK9IzgZsp5vGvFKvWSk2luGX7MPACMAMYlIJ8IsXu2+coVprn4p+v/Yp/KYGZ1QVJkyl+gcJ7aj0Ws/b8fz5mZmYZDkszM7MM34Y1MzPL8MrSzMwsw5+z7IdGjBgRY8aMqfUwzMwayoIFC56PiLK/ptBh2Q+NGTOG5ubmWg/DzKyhSHqqozrfhjUzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluFfStAPLV7dypgL7qz1MBreyunH1noIZlYnvLI0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDKyYSnpSklnlRzPlXRtyfHXJYWkvdq1mynpvE76vVDSMkmPSzqypHwbSbMk/UnSHyV9KJV/QdJSSS2Sfi1p9w76nSZpanp8vaQVkhZKWiTpsMxcP5fGFJJGlJRL0lWprkXS/iV1Z0t6TNISSbMlDe2k/8mSdulsDJnxfbG7bc3MrPsqWVk+CEwAkDQIGAG8vaR+AnA3MKmtIJ13InBTuQ4l7ZPOfztwFHCNpMGp+iJgbUTsCewD3J/KHwWaIuKdwK3A1yoYO8C5ETEeOAv4bubcB4DDgafalR8NjE1fU4DvpHmMAs5M43oHMJiS56GMyUC3wxJwWJqZ1UAlYfkAKSwpwm0J8KKknSQNAfYGzue1IfE+YGVEtA+dNhOBn0TExohYASwDDkx1HwcuA4iIzRHxfHp8b0RsSOc8BIyuZIIl5gGjOjshIh6NiJUdjPeGKDwE7ChpZKrbCthW0lbAMODZcn1LOhFoAm5MK91tJR0g6X5JC9KKfaSk4Wm1vVdqN1vSGZKmp+sslHRjmf6nSGqW1LxpQ2tlz4iZmVUkG5YR8SzwqqTdKEJzHjAfOITih39LRLQAmyXtm5pNAmZ30u0o4JmS41XAKEk7puMvS3pE0i2S3lSm/enAL3Jjb+co4I4utmlTdrwRsRq4AngaWAO0RsQ95TqIiFuBZuDDaaX7KvAt4MSIOAC4DvhKRLQCnwOulzQJ2Ckivh8RFwAvRcT4iPhwmf5nRURTRDQNHja8m9M0M7NyKt3g07a6bAvLeSXHD6ZzZgOT0gprInBLJ/2pTFlQrNJGAw9ExP7pOle8pqH0EYqQvrzCsV8u6Ungx8BXK2xT0Xgl7UQx1zdT3F7dLo2vEnsB7wB+KWkhcDFptRwRvwQWA98GPtHNMZuZWQ+pNCzb3rccR3Eb9iGKleUEiiCFIixPonjPryUi1nbS3ypg15Lj0RS3L9cBG4DbU/ktQOlmmsMp3tM8LiI2prKvpFuTCzu41rnAWynC6IcVzLUr4z0cWBERz0XEK8BtbLllnSPgsbRSHB8R4yLiCPjHe757Ay8Br+/mmM3MrId0ZWX5AeCFiNgUES8AO1IE5jyAiFhOEXbT6fwWLMAcilXoEElvptg48/uICODnwKHpvMOApQCS9gO+RxGU/wjiiLioLXA6ulhEbAa+CQwq3XnbBXOAU9Ou2IMpbreuobj9erCkYZKUxvuHTvp5Edg+PX4ceIOkQ9L8tpbUtnHq7NTPKcB1krZO5a+UPDYzsz5SaVguptgF+1C7sta2DTjJbOBtbFkZlhURjwE3UwTh3cBnI2JTqj4fmCapBfgocE4qvxx4HXBLWknOqXDsbdcM4FKgs4+znClpFcXKsaXkIzJ3AU9SbET6PvCZ1Od8ip25j1A8H4OAWZ0M43rgu2kVPJhix/AMSYuAhcAESXtS3Ho9JyJ+C/yGYlVM6rul3AYfMzPrPSoyxPqTISPHxsjTZtZ6GA1v5fRjaz0EM+tDkhZERFO5Ov8GHzMzs4yterPz9P7gjHbFKyLihN68bo6k2yl2sJY6PyLm9lD/3wbe3a74mxHxf3qifzMz61u9GpYpfHokgHpSb4d1RHy2N/s3M7O+5duwZmZmGQ5LMzOzDIelmZlZhsPSzMwso1c3+FhtjBs1nGZ/RtDMrMd4ZWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhneDdsPLV7dypgL7qz1MKyB+S+umL2WV5ZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWUa/DEtJ69P3MZJekrRQ0lJJN0jaupN2O0u6V9J6SVe3qztA0mJJyyRdJUmpfLfU5lFJLZKO6aT/MZL+tYp5HSppQnfbm5lZ9/TLsGxneUSMB8YBo4GTOjn3ZeDfgall6r4DTAHGpq+jUvnFwM0RsR8wCbimk/7HAN0OS+BQwGFpZtbHBkJYAhARm4DfA6M6OedvEfE7itD8B0kjgR0iYl5EBHADcHxbM2CH9Hg48Gwnw5gOvDetdM+WNFjS5ZIeTqvST6brfUHSdenxOElLJO0DfAo4O7V/b7sxTpHULKl504bWip4TMzOrzID5E12ShgIHAf/WjeajgFUlx6vYErrTgHskfR7YDji8k34uAKZGxAfSmKYArRHxLklDgAck3QPMBO6TdAJwEfDJiFgq6bvA+oi4on3HETELmAUwZOTY6MYczcysAwNhZbmHpIXAOuDpiGjpRh8qU9YWSKcA10fEaOAY4EeSKn1ejwBOTeObD+wMjI2IzcBk4EfA/RHxQDfGbGZmPWQgrCyXR8T4dCv1PknHRcScLvaxiuL9zjaj2XK79XTS+5cRMS+tYEcAayvoV8DnI2JumbqxwHpgly6O1czMethAWFkCEBFrKG6DXtjNti9KOjjtgj0V+Fmqfho4DEDS3sBQ4LkOunoR2L7keC7w6bYdupL2lLSdpOHAN4H3ATtLOrGD9mZm1gcGTFgmdwDD2m+OKSVpJfANYLKkVWljDcCngWuBZcBy4Bep/BzgDEmLgNnA5LQJqJwW4FVJiySdnfpbCjwiaQnwPYrV/pXANRHxJ4qV63RJbwR+DpxQboOPmZn1HnX8c90a1ZCRY2PkaTNrPQxrYCunH1vrIZj1OUkLIqKpXN1AW1mamZl12UDY4PPfSDoSmNGueEVEnNBD/Y+j2MlaamNEHNQT/ZuZWd8akGGZdp+W24HaU/0vBsb3Vv9mZta3fBvWzMwsw2FpZmaW4bA0MzPLcFiamZllDMgNPv3duFHDafbn5MzMeoxXlmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluHdsP3Q4tWtjLngzloPwwYo/8US64+8sjQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMuoqLCWtb3c8WdLV6fE0SaslLZS0VNIpqXyFpL3atZsp6by+G3llJI2XdEwV7Y+XtE9PjsnMzPLqKiwrcGVEjAcmAt+TtDXwE2BS2wmSBgEnAjf1xgAkDa6i+Xig22EJHA84LM3M+lijhSUAEfEEsAHYCZhNSVgC7wNWRsRT5dqm1eptku6W9ISkr5XUnSJpsaQlkmaUlK+XdImk+cAh6XiGpAWSfiXpQEn3SXpS0nEdXHcb4BLg5LQ6PlnSdpKuk/SwpEclTUznXiXpS+nxkZJ+I2kCcBxweWq/R7v+p0hqltS8aUNrl59TMzPrWL39ia5tJS0sOX49MKf9SZL2B56IiLXAWkmbJe0bEYsognN25jrjgf2AjcDjkr4FbAJmAAcAfwHukXR8RNwBbAcsiYi2ANsOuC8izpd0O3Ap8M8Uq74flhtzRPw9BWBTRHwu9fNV4D8j4uOSdgR+L+lXwAXAw5J+C1wFHBMRyyXNAf5fRNxapv9ZwCyAISPHRmb+ZmbWBfUWli+l26xAsQoEmkrqz5Z0BvAW4KiS8tnAJEmPUdyi/VLmOr+OiNZ0jaXA7sDOFAH4XCq/kWKVegdFkP60pP3fgbvT48XAxoh4RdJiYEyFcwU4AjhO0tR0PBTYLSL+kOb5G+DsiFjehT7NzKyH1VtY5lwZEVdI+l/ADZL2iIiXKcLyHuB+oCWtODuzseTxJornQZ2c/3JEbCo5fiUi2lZvm9v6i4jNkrrynAr4UEQ8XqZuHLAO2KUL/ZmZWS9o1PcsbwOagdPS8XKKYJlO/hZsR+YD75c0Im3iOYUifHvSi8D2Jcdzgc9LEoCk/dL33YFzKG4VHy3poA7am5lZH2jIsEwuAb6Qdr9CEZJvA27vTmcRsQa4ELgXWAQ8EhE/64mBlrgX2Kdtgw/wZWBroEXSEuDLKTh/AEyNiGeB04FrJQ2l2Pl7btoMtEcH1zAzsx6mLXcTrb8YMnJsjDxtZq2HYQPUyunH1noIZt0iaUFENJWra+SVpZmZWZ9otA0+FZN0JMVHQUqtiIgT+vO1zcys5/XbsIyIuRQbaAbUtc3MrOf5NqyZmVmGw9LMzCzDYWlmZpbhsDQzM8votxt8BrJxo4bT7M+6mZn1GK8szczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw7th+6HFq1sZc8GdtR6GWV3yX0Wx7vDK0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFh2g6T16fsYSS9JWihpqaQbJG3dSbudJd0rab2kq9vVHSBpsaRlkq6SpFS+W2rzqKQWScf07uzMzKw9h2X1lkfEeGAcMBo4qZNzXwb+HZhapu47wBRgbPo6KpVfDNwcEfsBk4BrembYZmZWKYdlD4mITcDvgVGdnPO3iPgdRWj+g6SRwA4RMS8iArgBOL6tGbBDejwceLZc35KmSGqW1LxpQ2tVczEzs9dyWPYQSUOBg4C7u9F8FLCq5HgVW0J3GvARSauAu4DPl+sgImZFRFNENA0eNrwbQzAzs444LKu3h6SFwDrg6Yho6UYfKlMW6fspwPURMRo4BviRJL9uZmZ9yD90q9f2nuVbgYMlHdeNPlZRvN/ZZjRbbreeDtwMEBHzgKHAiG6P1szMusxh2UMiYg1wAXBhN9u+KOngtAv2VOBnqfpp4DAASXtThOVzPTJoMzOriMOyZ90BDJP03o5OkLQS+AYwWdIqSfukqk8D1wLLgOXAL1L5OcAZkhYBs4HJaROQmZn1ka1qPYBGFBGvS99XAu8oKQ9g30zbMR2UN5f2VVK+FHh390drZmbV8srSzMwswyvLXiDpSGBGu+IVEXFCLcZjZmbVcVj2goiYC8yt9TjMzKxn+DasmZlZhsPSzMwsw2FpZmaW4bA0MzPL8AaffmjcqOE0Tz+21sMwM+s3vLI0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGPzrSDy1e3cqYC+6s9TDMrBes9MfCasIrSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZRs3DUtKVks4qOZ4r6dqS469LCkl7tWs3U9J5nfR7oaRlkh6XdGRJ+TaSZkn6k6Q/SvpQKv+CpKWSWiT9WtLuHfQ7TdLU9Ph6SSskLZS0SNJhmbl+Lo0pJI0oKZekq1Jdi6T9S+rOlvSYpCWSZksa2tk1zMys59U8LIEHgQkAkgYBI4C3l9RPAO4GJrUVpPNOBG4q16GkfdL5bweOAq6RNDhVXwSsjYg9gX2A+1P5o0BTRLwTuBX4WoXjPzcixgNnAd/NnPsAcDjwVLvyo4Gx6WsK8J00j1HAmWlc7wAGU/I8mJlZ36iHsHyAFJYU4bYEeFHSTpKGAHsD5/PakHgfsDIi2odOm4nATyJiY0SsAJYBB6a6jwOXAUTE5oh4Pj2+NyI2pHMeAkZ3cR7zgFGdnRARj0bEyg7Ge0MUHgJ2lDQy1W0FbCtpK2AY8Gy5viVNkdQsqXnThtYuDt3MzDpT87CMiGeBVyXtRhGa84D5wCFAE9ASES3AZkn7pmaTgNmddDsKeKbkeBUwStKO6fjLkh6RdIukN5Vpfzrwiy5O5Sjgji62aVN2vBGxGrgCeBpYA7RGxD3lOoiIWRHRFBFNg4cN7+YwzMysnJqHZdK2umwLy3klxw+mc2YDk9IKayJwSyf9qUxZUKzSRgMPRMT+6TpXvKah9BGKkL68wrFfLulJ4MfAVytsU9F4Je1EMdc3A7sA26XxmZlZH6qXsGx733IcxW3YhyhWlhMoghSKsDyJ4j2/lohY20l/q4BdS45HU9y+XAdsAG5P5bcApZtpDqd4T/O4iNiYyr6SNvAs7OBa5wJvBS4GfljBXLsy3sOBFRHxXES8AtzGllvWZmbWR+olLB8APgC8EBGbIuIFYEeKwJwHEBHLKcJuOp3fggWYQ7EKHSLpzRQbZ34fEQH8HDg0nXcYsBRA0n7A9yiC8h9BHBEXRcT4tImnrIjYDHwTGFS687YL5gCnpl2xB1Pcbl1Dcfv1YEnDJCmN9w/d6N/MzKpQL2G5mGIX7EPtylrbNuAks4G3sWVlWFZEPAbcTBGEdwOfjYhNqfp8YJqkFuCjwDmp/HLgdcAtaSU5pysTSEF8KdDZx1nOlLSKYuXYUvIRmbuAJyk2In0f+Ezqcz7FztxHKJ6PQcCsrozLzMyqp+JnvPUnQ0aOjZGnzaz1MMysF6ycfmyth9BvSVoQEU3l6uplZWlmZla3tqr1AKqR3h+c0a54RUScUIvxtJF0O8UO1lLnR8TcWozHzMyq09BhmcKn7gKo1mFtZmY9y7dhzczMMhyWZmZmGQ5LMzOzDIelmZlZRkNv8LHyxo0aTrM/i2Vm1mO8sjQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYY/OtIPLV7dypgL7qz1MMzM+lRv/vkyryzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkDKiwlrW93PFnS1enxNEmrJS2UtFTSKal8haS92rWbKem8vhu5mZnV0oAKywpcGRHjgYnA9yRtDfwEmNR2gqRBwInATb0xAEmDe6NfMzPrPodlGRHxBLAB2AmYTUlYAu8DVkbEU+XaptXqbZLulvSEpK+V1J0iabGkJZJmlJSvl3SJpPnAIel4hqQFkn4l6UBJ90l6UtJxHVx3iqRmSc2bNrT2wLNgZmZtBlpYbptusy6UtBC4pNxJkvYHnoiItRHRAmyWtG+qnkQRoJ0ZD5wMjANOlrSrpF2AGcD/TPXvknR8On87YElEHBQRv0vH90XEAcCLwKXAPwMndDTmiJgVEU0R0TR42PDc82BmZl0w0P5E10vpNitQrAKBppL6syWdAbwFOKqkfDYwSdJjFLdov5S5zq8jojVdYymwO7AzRQA+l8pvpFil3gFsAn5a0v7vwN3p8WJgY0S8ImkxMKbCuZqZWQ8ZaCvLnCsjYi+KVeENkoam8tnAScDhQEtErM30s7Hk8SaK/ylRJ+e/HBGbSo5fiYhIjze39RcRmxl4/4NjZlZzDssyIuI2oBk4LR0vB9YB08nfgu3IfOD9kkakTTynAPf3wHDNzKyXOSw7dgnwhbT7FYqQfBtwe3c6i4g1wIXAvcAi4JGI+FlPDNTMzHqXttzts/5iyMixMfK0mbUehplZn1o5/diq2ktaEBFN5eq8sjQzM8vwZpFuknQkxUdBSq2IiBNqMR4zM+s9Dstuioi5wNxaj8PMzHqfb8OamZllOCzNzMwyHJZmZmYZDkszM7MMb/Dph8aNGk5zlZ83MjOzLbyyNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLUETUegzWwyS9CDxe63H0gRHA87UeRC8bCHMEz7M/aeQ57h4RbyhX4V931z89HhFNtR5Eb5PU3N/nORDmCJ5nf9Jf5+jbsGZmZhkOSzMzswyHZf80q9YD6CMDYZ4DYY7gefYn/XKO3uBjZmaW4ZWlmZlZhsPSzMwsw2HZYCQdJelxScskXVCmXpKuSvUtkvavtG29qHKOKyUtlrRQUnPfjrxrKpjn2yTNk7RR0tSutK0nVc6zIV7PCub44fRvtUXSg5L2rbRtPalyng3xWnYoIvzVIF/AYGA58BZgG2ARsE+7c44BfgEIOBiYX2nbeviqZo6pbiUwotbz6KF5vhF4F/AVYGpX2tbLVzXzbJTXs8I5TgB2So+PbrT/LqudZ6O8lp19eWXZWA4ElkXEkxHxd+AnwMR250wEbojCQ8COkkZW2LYeVDPHRpKdZ0SsjYiHgVe62raOVDPPRlHJHB+MiL+kw4eA0ZW2rSPVzLPhOSwbyyjgmZLjVamsknMqaVsPqpkjQAD3SFogaUqvjbJ61bwejfJaQvVjbYTXs6tzPJ3izkh32tZSNfOExngtO+Rfd9dYVKas/Wd/Ojqnkrb1oJo5Arw7Ip6V9Ebgl5L+GBG/6dER9oxqXo9GeS2h+rE2wutZ8Rwl/Q+KEHlPV9vWgWrmCY3xWnbIK8vGsgrYteR4NPBshedU0rYeVDNHIqLt+1rgdopbR/WomtejUV5LqHKsDfJ6VjRHSe8ErgUmRsS6rrStE9XMs1Feyw45LBvLw8BYSW+WtA0wCZjT7pw5wKlpx+jBQGtErKmwbT3o9hwlbSdpewBJ2wFHAEv6cvBdUM3r0SivJVQx1gZ6PbNzlLQbcBvw0Yj4U1fa1pFuz7OBXssO+TZsA4mIVyV9DphLsTPtuoh4TNKnUv13gbsodosuAzYAH+usbQ2m0alq5gi8CbhdEhT/tv9vRNzdx1OoSCXzlPRPQDOwA7BZ0lkUuw//qxFeS6hunhR/6qnuX88K/81+CdgZuCbN59WIaGqU/y6hunnSQP9tdsS/7s7MzCzDt2HNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8v4/0kFA2ZHMr61AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.barh(np.arange(6), f_importance_ordered, xerr=std[inds], align=\"center\")\n",
    "plt.yticks(np.arange(6), 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
}