08_RF_model_balanced_accuracy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Random Forest model with balanced accuracy optimization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again a classical random forest approach but this time we try to optimize for the balanced accuracy, i.e. the average of the recall for each class in order to compensate for the imbalance between classes "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import pickle\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.model_selection import ShuffleSplit\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import plot_confusion_matrix\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some variables"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Due to limitations in the computation power we do not do a grid search as exhaustive as would be desirable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = '/data/ml_course/05_Capstone_project/rf_data/'\n",
    "model_name = 'RF_model_balanced_accuracy.pickle'\n",
    "test_size = 0.1 # 10% of points used for testing\n",
    "val_size = 0.1  # 10% of points (after removing test points) used for validation\n",
    "n_splits = 5  # number of splits in folding\n",
    "n_estimators = np.r_[np.arange(1, 21, 2), np.arange(20, 60, 10)]  # number of trees\n",
    "\n",
    "features = ['HRV_norm', 'HRV_norm_text', 'IR_108', 'IR_108_text', 'WV_062-IR_108', 'WV_062-IR_108_text']\n",
    "nfeatures = len(features)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "202009\r"
     ]
    }
   ],
   "source": [
    "years = ['2018', '2019', '2020']\n",
    "months = ['04', '05', '06', '07', '08', '09']\n",
    "\n",
    "for year in years:\n",
    "    for month in months:\n",
    "        print(year+month, end=\"\\r\", flush=True)\n",
    "        fname = datapath+year+month+'_data.npz'\n",
    "        if not os.path.isfile(fname):\n",
    "            continue\n",
    "        \n",
    "        X = None\n",
    "        with np.load(fname, allow_pickle=False) as npz_file:\n",
    "            # Load the arrays\n",
    "            if X is None:\n",
    "                X = npz_file['features']\n",
    "                y = npz_file['targets']\n",
    "            else:\n",
    "                X = np.concatenate((X, npz_file['features']), axis=0)\n",
    "                y = np.concatenate((y, npz_file['targets']), axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix shape: (7435724, 6)\n",
      "target matrix shape: (7435724,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix shape:', X.shape)\n",
    "print('target matrix shape:', y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove NaNs\n",
    "# NaNs can appear in texture fields.\n",
    "# They could safely be considered 0s but we remove them\n",
    "ind_row, ind_col = np.where(np.isnan(X))\n",
    "X = np.delete(X, ind_row, axis=0)\n",
    "y = np.delete(y, ind_row, axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix shape: (7432290, 6)\n",
      "target matrix shape: (7432290,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix shape:', X.shape)\n",
    "print('target matrix shape:', y.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1], dtype=int8)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Group together 0 (No POH was computed) and 1 (POH below 90%)\n",
    "np.unique(y)\n",
    "y[y == 1] = 0\n",
    "y[y == 2] = 1\n",
    "np.unique(y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "npixels: 7432290\n",
      "npixels with hail: 61171\n",
      "% of pixels with hail: 0.8230437725115678\n"
     ]
    }
   ],
   "source": [
    "npixels = y.size\n",
    "npixels_hail = y[y==1].size\n",
    "\n",
    "print('npixels:', npixels)\n",
    "print('npixels with hail:', npixels_hail)\n",
    "print('% of pixels with hail:', npixels_hail/npixels*100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The dataset is highly imbalanced toward non-hail pixels!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Split data into test and training"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split data\n",
    "X_tr, X_te, y_tr, y_te = train_test_split(X, y, test_size=test_size, random_state=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "features matrix for training shape: (6689061, 6)\n",
      "target matrix for training shape: (6689061,)\n"
     ]
    }
   ],
   "source": [
    "print('features matrix for training shape:', X_tr.shape)\n",
    "print('target matrix for training shape:', y_tr.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prepare data folds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ss = ShuffleSplit(n_splits=n_splits, test_size=val_size, random_state=0)\n",
    "ss.get_n_splits()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitting 5 folds for each of 14 candidates, totalling 70 fits\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=-1)]: Using backend LokyBackend with 4 concurrent workers.\n",
      "[Parallel(n_jobs=-1)]: Done  42 tasks      | elapsed: 102.8min\n",
      "[Parallel(n_jobs=-1)]: Done  70 out of  70 | elapsed: 323.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(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": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create pipeline\n",
    "# Weighted random forest: modifies the weight by penalizing more heavily the under-represented class\n",
    "# balanced: modifies the weight according to the overall distribution\n",
    "# balanced_subsample: modies the weight according to the distribution in the bootstrapped subsamples\n",
    "pipe_RF = Pipeline([\n",
    "    ('scaler', None),\n",
    "    ('RF', RandomForestClassifier(random_state=0, max_depth=None, class_weight=None))\n",
    "])\n",
    "\n",
    "# Create cross-validation object\n",
    "grid_RF = {'RF__n_estimators': n_estimators}\n",
    "\n",
    "grid_RF_cv = GridSearchCV(pipe_RF, grid_RF, 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": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect results in a DataFrame\n",
    "RF_results = pd.DataFrame({\n",
    "    'trees': grid_RF_cv.cv_results_['param_RF__n_estimators'],\n",
    "    'mean_tr': grid_RF_cv.cv_results_['mean_train_score'],\n",
    "    'mean_te': grid_RF_cv.cv_results_['mean_test_score'],\n",
    "    'std_te': grid_RF_cv.cv_results_['std_test_score']\n",
    "})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkxElEQVR4nO3de3yU5Z338c9vJodJICEJBIRMaMCCisghhojrCbVWbK1WW7tSrW3dLqtdq32e3e3j7j72sE8P60vX1Xa1lqW09mDdtupKW4uHrecTB0GUg0oBIYAQDuGQI0l+zx8ziUMMMsAM92Tm+3698srch5n5Xfjyy8V1X/d1m7sjIiLZKxR0ASIikl4KehGRLKegFxHJcgp6EZEsp6AXEclyeUEX0J9hw4Z5TU1N0GWIiAwYS5Ys2e7ulf0dy8igr6mpYfHixUGXISIyYJjZOwc7pqEbEZEsp6AXEclyCnoRkSyXkWP0IpL99u/fT0NDA21tbUGXMqBEIhGi0Sj5+flJv0dBLyKBaGhooKSkhJqaGsws6HIGBHdnx44dNDQ0MGbMmKTfp6EbEQlEW1sbQ4cOVcgfBjNj6NChh/2vIAW9iARGIX/4juTPLGuCvr2zi3uf+TPPvd0YdCkiIhklqaA3s5lm9qaZrTGzmw9yzgwzW2ZmK8zsmfi+E+L7en72mNlXU1h/r4JwiDnPruW/l25Ox8eLiKTFnXfeSUtLS1q/45BBb2Zh4G7gImACMMvMJvQ5pwy4B7jE3U8GrgBw9zfdfYq7TwFOBVqAh1PZgIQaqPtQOYvW70zHx4uIpEVGBD1QD6xx97Xu3gE8AFza55zPAg+5+wYAd9/Wz+ecD/zZ3Q96m+7Rqh9TwYadLWzdo+laInJo69ev58QTT+RLX/oSEydO5KqrruLJJ5/kjDPOYNy4cSxcuJDm5mauvfZapk2bxtSpU3nkkUd633vWWWdRW1tLbW0tL774IgBPP/00M2bM4NOf/jQnnngiV111FQd7kt/3v/99Nm/ezLnnnsu5554LwOOPP87pp59ObW0tV1xxBfv27TvqdiYzvbIK2Jiw3QCc1uec8UC+mT0NlAB3ufvP+pxzJfCrg32Jmc0GZgOMHj06ibLeb1pNBQCL1u/k4kmjjugzROTY+9bvVrBy856UfuaEUaV84xMnH/K8NWvW8Jvf/IY5c+Ywbdo07r//fp5//nnmz5/Pd7/7XSZMmMB5553HvHnzaGpqor6+no985CMMHz6cJ554gkgkwttvv82sWbN61+haunQpK1asYNSoUZxxxhm88MILnHnmme/77htvvJE77riDp556imHDhrF9+3a+/e1v8+STTzJo0CBuvfVW7rjjDr7+9a8f1Z9FMkHf3yXevn895REbmjkfKAJeMrOX3f0tADMrAC4B/vFgX+Luc4A5AHV1dUf0INsJo0opyg+zaJ2CXkSSM2bMGE455RQATj75ZM4//3zMjFNOOYX169fT0NDA/Pnzuf3224HYtNANGzYwatQobrjhBpYtW0Y4HOatt97q/cz6+nqi0SgAU6ZMYf369f0GfV8vv/wyK1eu5IwzzgCgo6OD008//ajbmEzQNwDVCdtRoO8VzwZgu7s3A81m9iwwGehp+UXAq+6+9Sjr/UD54RC1Hypj4fpd6fwaEUmxZHre6VJYWNj7OhQK9W6HQiE6OzsJh8M8+OCDnHDCCQe875vf/CYjRozgtddeo7u7m0gk0u9nhsNhOjs7k6rF3bngggv41a8OOvhxRJIZo18EjDOzMfGe+ZXA/D7nPAKcZWZ5ZlZMbGhnVcLxWXzAsE0qTaupYPW7e9jduv9YfJ2IZLkLL7yQH/zgB73j7EuXLgVg9+7djBw5klAoxM9//nO6urqO6PNLSkrYu3cvANOnT+eFF15gzZo1ALS0tBzwL4Ujdcigd/dO4AbgMWLh/Wt3X2Fm15nZdfFzVgELgOXAQmCuu78BEA/+C4CHjrraJEyrqcAdXt2gXr2IHL1bbrmF/fv3M2nSJCZOnMgtt9wCwJe//GXuu+8+pk+fzltvvcWgQYOO6PNnz57NRRddxLnnnktlZSU//elPmTVrFpMmTWL69OmsXr36qNtgB7saHKS6ujo/0gePtHR0MumbjzP77LF8beaJKa5MRFJl1apVnHTSSUGXMSD192dnZkvcva6/87PmztgexQV5nFw1RPPpRUTisnL1yvqacu578R3a9ncRyQ8HXY6ICJdddhnr1q07YN+tt97KhRdemPbvzsqgr6up4D+fW8frm3b3zq0XEQnSww+nZVGApGTd0A28d+PUwnUavhHJZJl4jTDTHcmfWVYGfcWgAj48fDCLNU4vkrEikQg7duxQ2B+GngePJM7ZT0ZWDt0ATKsp5/fLt9DV7YRDWvNaJNNEo1EaGhpobNTS4oej51GChyOLg76CXy3cyJvv7mXCqNKgyxGRPvLz8w/rcXhy5LJy6AYOXOBMRCSXZW3QR8uLGDkkoqAXkZyXtUFvZtTVVLBo/U5d7BGRnJa1QQ+xG6e27mln487WoEsREQlMVgf9tDEapxcRyeqgHz+8hNJInoJeRHJaVgd9KBQbp1+ooBeRHJbVQQ+xaZZrG5vZvq896FJERAKR9UFfP6YcgMV6vKCI5KisD/qJVUMoyAtpnF5EclbWB31hXpgp1WUKehHJWVkf9AD1NRWs2LyH5vbknsQuIpJNciLo62rK6ep2lm5oCroUEZFjLieC/tQPlRMyNM1SRHJSTgR9SSSfk0aWskhPnBKRHJQTQQ+x+fRLN+5if1d30KWIiBxTORX0bfu7eWPT7qBLERE5pnIn6OM3TmmapYjkmpwJ+uElEWqGFrNId8iKSI7JmaAHqKupYPH6nXR360EkIpI7ciro62sq2NWynz837gu6FBGRYyapoDezmWb2ppmtMbObD3LODDNbZmYrzOyZhP1lZvZbM1ttZqvM7PRUFX+43nsQiYZvRCR3HDLozSwM3A1cBEwAZpnZhD7nlAH3AJe4+8nAFQmH7wIWuPuJwGRgVWpKP3w1Q4sZNrhQF2RFJKck06OvB9a4+1p37wAeAC7tc85ngYfcfQOAu28DMLNS4Gzgx/H9He7elKLaD5uZMa2mnIW6cUpEckgyQV8FbEzYbojvSzQeKDezp81siZldE98/FmgEfmJmS81srpkN6u9LzGy2mS02s8WNjY2H2YzkTaupYFNTK5ub9MBwEckNyQS99bOv77SVPOBU4OPAhcAtZjY+vr8W+KG7TwWagX7H+N19jrvXuXtdZWVlsvUftno9MFxEckwyQd8AVCdsR4HN/ZyzwN2b3X078Cyx8fgGoMHdX4mf91tiwR+YE48rYVBBWEEvIjkjmaBfBIwzszFmVgBcCczvc84jwFlmlmdmxcBpwCp3fxfYaGYnxM87H1iZotqPSF44RO2Hylm0TjNvRCQ3HDLo3b0TuAF4jNiMmV+7+wozu87MroufswpYACwHFgJz3f2N+Ed8BfilmS0HpgDfTXkrDlN9TQVvbt3L7pb9QZciIpJ2ecmc5O6PAo/22Xdvn+3bgNv6ee8yoO7IS0y9nvn0i9/ZyfknjQi4GhGR9MqpO2N7TKkuIz9sehCJiOSEnAz6SH6YU6qG6EEkIpITcjLoITZ88/qm3bTt7wq6FBGRtMrdoP9QBfu7nGUbm4IuRUQkrXI26Otq4g8i0fCNiGS5nA36suICThhRwqJ3NJ9eRLJbzgY9xB4vuGT9To3Ti0hWy+mgn3nySJo7unhy1dagSxERSZucDvrTjx/KcaURHnp1U9CliIikTU4HfThkfHJqFc+81Ujj3vagyxERSYucDnqAy2ur6Op25r/Wd0FOEZHskPNBP35ECadUDeGhVxuCLkVEJC1yPugh1qtfsXkPq9/dE3QpIiIpp6AHPjF5FHkh42FdlBWRLKSgB4YNLmTGCZU8vHQTXd19n5IoIjKwKejjLq+Nsm1vOy+s2R50KSIiKaWgjzvvxOGURvJ0UVZEso6CPi6SH+biyaNYsOJd9rV3Bl2OiEjKKOgTfKq2irb93fzx9S1BlyIikjIK+gS1o8v50NBiLYkgIllFQZ/AzLh8apSX1u6gYVdL0OWIiKSEgr6Py6ZWAfDIMi2JICLZQUHfx+ihxdTXVPDgqw24a069iAx8Cvp+XF5bxdrGZl5r2B10KSIiR01B34+PTRpJQV5Ic+pFJCso6PtRGsnnoxNGMP+1zXR0dgddjojIUVHQH8SnaqM0teznqTe3BV2KiMhRUdAfxFnjhjFscIGGb0RkwEsq6M1sppm9aWZrzOzmg5wzw8yWmdkKM3smYf96M3s9fmxxqgpPt7xwiEunVPGn1dvY1dwRdDkiIkfskEFvZmHgbuAiYAIwy8wm9DmnDLgHuMTdTwau6PMx57r7FHevS0nVx8jltVXs73J+v1xz6kVk4EqmR18PrHH3te7eATwAXNrnnM8CD7n7BgB3z4qB7QkjSznxuBIe1JIIIjKAJRP0VcDGhO2G+L5E44FyM3vazJaY2TUJxxx4PL5/9sG+xMxmm9liM1vc2NiYbP1pZWZcXlvFso1N/LlxX9DliIgckWSC3vrZ1/eW0TzgVODjwIXALWY2Pn7sDHevJTb087dmdnZ/X+Luc9y9zt3rKisrk6v+GLh0ShUhQ48ZFJEBK5mgbwCqE7ajQN9B6wZggbs3u/t24FlgMoC7b47/3gY8TGwoaMAYURrhzHGxxwx26zGDIjIAJRP0i4BxZjbGzAqAK4H5fc55BDjLzPLMrBg4DVhlZoPMrATAzAYBHwXeSF35x8anaqvY1NTKK+t2Bl2KiMhhyzvUCe7eaWY3AI8BYWCeu68ws+vix+9191VmtgBYDnQDc939DTMbCzxsZj3fdb+7L0hXY9LloxOOY1BBmIdebeD044cGXY6IyGGxTFyhsa6uzhcvzqwp9//wm9d49PUtLP6/F1BUEA66HBGRA5jZkoNNYdedsUm6vDZKc0cXj698N+hSREQOi4I+SaeNqaCqrEhz6kVkwFHQJykUMi6bWsXzbzeydU9b0OWIiCRNQX8YLq+twoFb/7haT58SkQFDQX8YxlYO5sbzxvHQ0k388pUNQZcjIpIUBf1huun8ccw4oZJv/W4FSzfsCrocEZFDUtAfplDIuPMvpzCiNMKXf/kqO/a1B12SiMgHUtAfgbLiAu69+lR2NHdw0wPL6NLSCCKSwRT0R2hi1RC+felEnl+znTueeDPockREDkpBfxQ+M62aWfXV3P3Un3l8hW6kEpHMpKA/St/4xMmcUjWEv/v1a6zb3hx0OSIi76OgP0qR/DA/vLqWcNi4/hdLaOnoDLokEZEDKOhTIFpezF1XTuXNrXv554ff0M1UIpJRFPQpcs74Sv7XR8bz8NJN/Pzld4IuR0Skl4I+hW4498Ocd+Jw/t/vV/KqbqYSkQyhoE+hUMj4989MYeSQIr78i1fZrpupRCQDKOhTbEhxPj+8upZdLR185f6ldHZ1B12SiOQ4BX0anDxqCN+57BReWruD2x9/K+hyRCTHKejT5NOnRrnqtNHc+8yfWfDGlqDLEZEcpqBPo69/YgKTq8v46n8tY+G6nUGXIyI5SkGfRoV5YX78+TpGlRVx7U8X8drGpqBLEpEcpKBPs2GDC7n/S9MpH5TPNfMWsmrLnqBLEpEco6A/Bo4bEuH+L02nKD/M5378Cmsb9wVdkojkEAX9MVJdUcwv//o0AK6a+wobd7YEXJGI5AoF/TF0fOVgfv5Xp9HS0cVVc19h6562oEsSkRygoD/GThpZyn3X1rNjXztXzX1FjyIUkbRT0AdgSnUZ874wjYZdLXzuxwvZ3bo/6JJEJIsp6ANy2tih/Ohzdby9bS9f+MlCmtu1jr2IpEdSQW9mM83sTTNbY2Y3H+ScGWa2zMxWmNkzfY6FzWypmf0+FUVni3PGV/KDWbUsb9jNl+5bTNv+rqBLEpEsdMigN7MwcDdwETABmGVmE/qcUwbcA1zi7icDV/T5mJuAVakoONvMnHgc/3bFZF5et4Prf7GEjk4tgiYiqZVMj74eWOPua929A3gAuLTPOZ8FHnL3DQDuvq3ngJlFgY8Dc1NTcvb55NQqvnvZKTz1ZiNf/S+teCkiqZVM0FcBGxO2G+L7Eo0Hys3saTNbYmbXJBy7E/ga8IHpZWazzWyxmS1ubGxMoqzsMqt+NLdcPIFHX3+Xr/12Od3dehyhiKRGXhLnWD/7+qZQHnAqcD5QBLxkZi8T+wtgm7svMbMZH/Ql7j4HmANQV1eXkyn3V2eOoaW9k3974i3Kigu45eKTMOvvj19EJHnJBH0DUJ2wHQU293POdndvBprN7FlgMlALXGJmHwMiQKmZ/cLdrz760rPTDed9mJ0tHcx7YR0jSgv5m3OOD7okERngkhm6WQSMM7MxZlYAXAnM73POI8BZZpZnZsXAacAqd/9Hd4+6e038fX9SyH8wM+OWj0/g4kkj+d4fV/PQqw1BlyQiA9whe/Tu3mlmNwCPAWFgnruvMLPr4sfvdfdVZrYAWE5sLH6uu7+RzsKzWShk/NtnJrOzuYOv/XY5FYMKmHHC8KDLEpEBytwzbzi8rq7OFy9eHHQZgdvbtp/P/Ohl3tnRzK/+ejqTq8uCLklEMpSZLXH3uv6O6c7YDFYSyee+L06jYlABX/zpItZtbw66JBEZgBT0GW54aYSfXVsPwDXzXmHbXq14KSKHR0E/AIytHMy8L0xj+94OvjBvEXvbtAiaiCRPQT9ATKku456ra3lr616u+8US2ju1Lo6IJEdBP4Cce8Jwbv3UJF5Ys4O//43unhWR5CRzw5RkkE+dGmXb3nZuXbCaysGFuntWRA5JQT8AXXfOWLbuadPdsyKSFAX9AGRmfP3iCTTua+d7f1xNZUkhl9dGgy5LRDKUgn6ACoWMOz4zmZ37YnfPhkPGJZNHaRhHRN5HF2MHsMK8MD+65lQmjCrlpgeW8akfvsgra3cEXZaIZBgF/QBXGsnnwev/gu9edgqbmlr5yzkv84WfLGTF5t1BlyYiGUJr3WSR1o4u7ntpPfc8tYY9bZ1cMnkU//uC8dQMGxR0aSKSZh+01o2CPgvtbtnPj579M/NeWEdnl3NlfTU3njeO4aWRoEsTkTRR0OeobXva+P6f3uaBhRvJD4e49swaZp99PEOK8oMuTURSTEGf49Zvb+aOJ95i/mubGVKUz5dnHM/n/6KGSH446NJEJEW0THGOqxk2iO/PmsofbjyT2tFlfO+Pqzn39qd56s1tQZcmIseAgj6HnDxqCD/5Yj3/NXs6JZE8vviTRdz84HKthimS5RT0Oei0sUP53VfO5LpzjufXizcy887neHHN9qDLEpE0UdDnqMK8MDdfdCK/vf4vKMwL8dm5r/D1R96gpaMz6NJEJMUU9DmudnQ5f7jxLK49Yww/e+kdLrrrORav3xl0WSKSQgp6oaggzNc/MYEHZk+n250rfvQS3/nDStr26+EmItlAQS+9po8dyh9vOptZ9aP5z+fWcfEPnue1jU1BlyUiR0lBLwcYXJjHdy87hZ9dW09zeyeX//BFbn/sTTo6u4MuTUSOkIJe+nX2+EoWfPVsPjmliv94ag2X/MfzPLJsE1v3tAVdmogcJt0ZK4f0xMqt/PPDr7NtbzsAY4cN4rSxQ5k+toLpY4cyQmvoiAROSyDIUevs6mbllj28vHYHL6/dyaJ1O9nbHpuKqeAXCZ6CXlKuq9tZubkn+HewsJ/gn1pdxnFDIhw3JMKI0gilkTw9AUskTRT0knYfFPw9ivLD8dAv5LjSCCOGRBhR8t5fBNHyIv1rQOQIHXXQm9lM4C4gDMx193/t55wZwJ1APrDd3c8xswjwLFBI7Pm0v3X3bxzq+xT0A19Xt7NpVyvv7mnj3T1tbN3d9r7X2/a009F14Gyen11bz9njKwOqWmTg+qCgP+TDwc0sDNwNXAA0AIvMbL67r0w4pwy4B5jp7hvMbHj8UDtwnrvvM7N84Hkz+6O7v3x0TZJMFw4Zo4cWM3po8UHPcXd2NnewdU87G3e18Dc/X8LKLXsU9CIpdsigB+qBNe6+FsDMHgAuBVYmnPNZ4CF33wDg7tvivx3YFz8nP/6TeWNFEggzY+jgQoYOLmTCqFLKivNp2NUSdFkiWSeZefRVwMaE7Yb4vkTjgXIze9rMlpjZNT0HzCxsZsuAbcAT7v5Kf19iZrPNbLGZLW5sbDysRkh2iJYXsXFna9BliGSdZIK+v2kSfXvlecCpwMeBC4FbzGw8gLt3ufsUIArUm9nE/r7E3ee4e52711VW6p/uuShaVqwevUgaJBP0DUB1wnYU2NzPOQvcvdndtxO7ADs58QR3bwKeBmYeabGS3aorimjY1UomzgQTGciSCfpFwDgzG2NmBcCVwPw+5zwCnGVmeWZWDJwGrDKzyviFWsysCPgIsDpl1UtWiZYX097ZzfZ9HUGXIpJVDnkx1t07zewG4DFi0yvnufsKM7sufvxed19lZguA5UA3sSmYb5jZJOC++MydEPBrd/992lojA1q0vAiAjbtaqCwpDLgakeyRzKwb3P1R4NE+++7ts30bcFuffcuBqUdZo+SI6orYVMyGXa3Uji4PuBqR7KHVKyVjVJXFevS6ICuSWgp6yRiDCvOoGFRAwy5NsRRJJQW9ZJTYXHr16EVSSUEvGaW6vJhN6tGLpJSCXjJKtLyIhqZWurs1l14kVRT0klGi5UV0dHbTuK896FJEsoaCXjJKtHeKpcbpRVJFQS8Zpbq8Z4qlxulFUkVBLxmlquy9m6ZEJDUU9JJRigrCDBtcoCmWIimkoJeMEy0vVo9eJIUU9JJxouVFuhgrkkIKesk40fJiNmkuvUjKKOgl40TLi9jf5Wzd2xZ0KSJZQUEvGSdxuWIROXoKesk40XItVyySSgp6yTg969Jv3KkevUgqKOgl40TywwwvKVSPXiRFFPSSkWJTLNWjF0kFBb1kJN00JZI6CnrJSNHyIjY3tdKlufQiR01BLxmpuqKYzm7n3T2aSy9ytBT0kpF6p1hqcTORo6agl4wULddNUyKpkhd0ASL9GVUWwQw2aoqlpIC709ntdMV/Ont/d9PdDZ3d3b37u7uP7NwDt7v7/YzE7W53OrveO7fbneKCPG65eELK26+gl4xUmBdmRElEPfos19rRxaamFhp2tdKwq5VNTa1s2tXKnrb9sUDsioejJwRoVzwke0Kzz3ZXQjD3vM6ka/rhkBE2Ixwy8kJGOPze9rDBhWn5TgW9ZCwtVzzw7W7dT8OuFjYlhPimpvdCfWdzxwHn54WMkWURyosLYiEY/ykIhXtf9wbk+7ZDhEOQFwr1f67FQrX3XINwOHTA8bzwe+8J9W6H+mwfeLz3dSjU+/mxbXv/dsgws2P+3yGpoDezmcBdQBiY6+7/2s85M4A7gXxgu7ufY2bVwM+A44BuYI6735WSyiXrRcuLWPzOrqDLkINwd7bv60gI8FjPPDHU97Z3HvCewrwQVeVFVJUVMbGqlGh5MVVlRb37RpRGCIeOfRBmu0MGvZmFgbuBC4AGYJGZzXf3lQnnlAH3ADPdfYOZDY8f6gT+zt1fNbMSYImZPZH4XpGDqa4o5nfLt9DZ1U1eWPMGjrWubmfrnrbe0G7Y1XJAb3zTrlbaO7sPeE9JYR5V5UVEy4s4bUxFPMCLe/cNHVQQSI821yXTo68H1rj7WgAzewC4FEgM688CD7n7BgB33xb/vQXYEn+918xWAVV93ivSr2h5EV3dzpbdbb1LF0vqtHd2saUpIcj79Mzf3d1GZ5/B7aGDCqgqL+KEESWcd8LweIC/1ysfUpQfUGvkgyQT9FXAxoTtBuC0PueMB/LN7GmgBLjL3X+WeIKZ1QBTgVeOtFjJLYlTLBX0h6+lo7NPgMd74/Ge+ba97XhCjpvBiJIIVeVF1I4uJ1r+3pBKtLyIUWVFFBfost5AlMx/tf7+ndX3GnYecCpwPlAEvGRmL7v7WwBmNhh4EPiqu+/p90vMZgOzAUaPHp1c9ZLVem6a2rirhdMZGnA1mcXd2dPaSUNTS3xY5cCLnR90oTNaVsxZ4yp7e+HR+O+RQ4ooyNMQWTZKJugbgOqE7SiwuZ9ztrt7M9BsZs8Ck4G3zCyfWMj/0t0fOtiXuPscYA5AXV1dBk2GkqCMHFJEyHLzpqnubmd7c/v7Zqskhvq+Phc6I/mheHgXM7FqCNH4uHhPoA8v0YXOXJVM0C8CxpnZGGATcCWxMflEjwD/YWZ5QAGxoZ1/t9hVlx8Dq9z9jtSVLbmgIC/EcaWRrJxi2drRxebdrWxuiv1samrrfb25qZXNu9vo6HuhM5JHVVkR1RVFTB9bERsbTwhyXeiUgzlk0Lt7p5ndADxGbHrlPHdfYWbXxY/f6+6rzGwBsJzYNMq57v6GmZ0JfA543cyWxT/yn9z90XQ0RrLPQFyuuKc3vjkhvDf1hnhbv8MqPePjo8oiTKwawoUnH8eosti4uC50ytFK6spKPJgf7bPv3j7btwG39dn3PP2P8YskJVpexMtrdwRdxgH6641v2tXTE29lS1MbHV0H9saLC8JUxYN7YtUQqsoiBwT5iNKIxsclbXQJXTJatKKYd5dtoqOz+5gEYd/eeM/YeE+Ib25qe19vPGQwojQW3JOiZcycGImF+pD3gry0KE/DKhIYBb1ktGh5Ed0O7+5uY/TQo59iGVtbpfWA8fDe8fGD9MYHFYSpik8vnBQti/fMI71BftyQCPm6oUsymIJeMlrvuvS7Wg4Z9N3dzvZ97fEgbztwbPwweuPR+JBKz09pRL1xGdgU9JLRquM3TW3c1UJLR+cBFzgTe+ObmlrZsruV/V0HzswdXJjX2wOfHC3rHUoZFd83olS9ccl+CnrJaCOHxOZ+3/LfK/g/D75+wLGQwXHx3viU6jI+dsrI3oucPUMtpRHNVBFR0EtGywuH+PuPnsCmppY+vfEiRpQUarEzkSQo6CXjXT/j+KBLEBnQ1B0SEclyCnoRkSynoBcRyXIKehGRLKegFxHJcgp6EZEsp6AXEclyCnoRkSxn7pn31D4zawTeOcRpw4Dtx6CcTKN25xa1O7ccTbs/5O6V/R3IyKBPhpktdve6oOs41tTu3KJ255Z0tVtDNyIiWU5BLyKS5QZy0M8JuoCAqN25Re3OLWlp94AdoxcRkeQM5B69iIgkQUEvIpLlBlzQm9lMM3vTzNaY2c1B15MuZjbPzLaZ2RsJ+yrM7Akzezv+uzzIGtPBzKrN7CkzW2VmK8zspvj+rG67mUXMbKGZvRZv97fi+7O63T3MLGxmS83s9/HtXGn3ejN73cyWmdni+L6Ut31ABb2ZhYG7gYuACcAsM5sQbFVp81NgZp99NwP/4+7jgP+Jb2ebTuDv3P0kYDrwt/H/xtne9nbgPHefDEwBZprZdLK/3T1uAlYlbOdKuwHOdfcpCfPnU972ARX0QD2wxt3XunsH8ABwacA1pYW7Pwvs7LP7UuC++Ov7gE8ey5qOBXff4u6vxl/vJfY/fxVZ3naP2RffzI//OFnebgAziwIfB+Ym7M76dn+AlLd9oAV9FbAxYbshvi9XjHD3LRALRGB4wPWklZnVAFOBV8iBtseHL5YB24An3D0n2g3cCXwN6E7Ylwvththf5o+b2RIzmx3fl/K2D7SHg1s/+zQ/NAuZ2WDgQeCr7r7HrL//9NnF3buAKWZWBjxsZhMDLintzOxiYJu7LzGzGQGXE4Qz3H2zmQ0HnjCz1en4koHWo28AqhO2o8DmgGoJwlYzGwkQ/70t4HrSwszyiYX8L939ofjunGg7gLs3AU8Tu0aT7e0+A7jEzNYTG4o9z8x+Qfa3GwB33xz/vQ14mNjwdMrbPtCCfhEwzszGmFkBcCUwP+CajqX5wOfjrz8PPBJgLWlhsa77j4FV7n5HwqGsbruZVcZ78phZEfARYDVZ3m53/0d3j7p7DbH/n//k7leT5e0GMLNBZlbS8xr4KPAGaWj7gLsz1sw+RmxMLwzMc/fvBFtRepjZr4AZxJYt3Qp8A/hv4NfAaGADcIW7971gO6CZ2ZnAc8DrvDdm+0/Exumztu1mNonYhbcwsQ7Yr939X8xsKFnc7kTxoZu/d/eLc6HdZjaWWC8eYsPo97v7d9LR9gEX9CIicngG2tCNiIgcJgW9iEiWU9CLiGQ5Bb2ISJZT0IuIZDkFveQ0Myszsy8HXYdIOinoJdeVAe8L+vhKqSJZQUEvue5fgePj64Eviq+Ffz/wenyRsdvi+5eb2d/0vMnM/iFhf8/a8YPM7A/xNeXfMLO/DKpRIokG2qJmIql2MzDR3afE78z8Q3x7XXw1wd3uPs3MCoEXzOxxYFz8p57YQnvzzexsoBLY7O4fBzCzIce+OSLvp6AXOdBCd18Xf/1RYJKZfTq+PYRYwH80/rM0vn9wfP9zwO1mdivwe3d/7tiVLXJwCnqRAzUnvDbgK+7+WOIJZnYh8D13/1HfN5vZqcDHgO+Z2ePu/i9prVYkCRqjl1y3Fyg5yLHHgOvjyyZjZuPjqww+BlwbXzMfM6sys+FmNgpocfdfALcDtekvX+TQ1KOXnObuO8zsBYs9hL2V2EqhPeYCNcCr8eWTG4FPuvvjZnYS8FL8gSj7gKuBDwO3mVk3sB+4/ti1ROTgtHqliEiW09CNiEiWU9CLiGQ5Bb2ISJZT0IuIZDkFvYhIllPQi4hkOQW9iEiW+/+E5sxv5nWcgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot test curve\n",
    "RF_results.plot(x='trees', y='mean_te')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>trees</th>\n",
       "      <th>mean_tr</th>\n",
       "      <th>mean_te</th>\n",
       "      <th>std_te</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0.877805</td>\n",
       "      <td>0.670746</td>\n",
       "      <td>0.006704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>3</td>\n",
       "      <td>0.921436</td>\n",
       "      <td>0.656556</td>\n",
       "      <td>0.003472</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>5</td>\n",
       "      <td>0.943111</td>\n",
       "      <td>0.647394</td>\n",
       "      <td>0.004272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>7</td>\n",
       "      <td>0.957388</td>\n",
       "      <td>0.643077</td>\n",
       "      <td>0.003147</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>0.967030</td>\n",
       "      <td>0.639547</td>\n",
       "      <td>0.002806</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>11</td>\n",
       "      <td>0.973599</td>\n",
       "      <td>0.635807</td>\n",
       "      <td>0.002325</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>13</td>\n",
       "      <td>0.978553</td>\n",
       "      <td>0.634020</td>\n",
       "      <td>0.003119</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>15</td>\n",
       "      <td>0.982266</td>\n",
       "      <td>0.631742</td>\n",
       "      <td>0.002071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>17</td>\n",
       "      <td>0.985342</td>\n",
       "      <td>0.630901</td>\n",
       "      <td>0.002694</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19</td>\n",
       "      <td>0.987744</td>\n",
       "      <td>0.630547</td>\n",
       "      <td>0.002910</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  trees   mean_tr   mean_te    std_te\n",
       "0     1  0.877805  0.670746  0.006704\n",
       "1     3  0.921436  0.656556  0.003472\n",
       "2     5  0.943111  0.647394  0.004272\n",
       "3     7  0.957388  0.643077  0.003147\n",
       "4     9  0.967030  0.639547  0.002806\n",
       "5    11  0.973599  0.635807  0.002325\n",
       "6    13  0.978553  0.634020  0.003119\n",
       "7    15  0.982266  0.631742  0.002071\n",
       "8    17  0.985342  0.630901  0.002694\n",
       "9    19  0.987744  0.630547  0.002910"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Ten best combinations according to the mean test score\n",
    "RF_results.sort_values(by='mean_te', ascending=False).head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 67.82%\n"
     ]
    }
   ],
   "source": [
    "# Report test score\n",
    "acc_RF = 100*grid_RF_cv.score(X_te, y_te)\n",
    "print('Test accuracy: {:.2f}%'.format(acc_RF))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = grid_RF_cv.predict(X_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn, fp, fn, tp = confusion_matrix(y_te, y_pred).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True positive:  2195\n",
      "True negative:  733177\n",
      "False positive:  3984\n",
      "False negative:  3873\n",
      "Positive pixels:  6068\n",
      "Negative pixels:  737161\n"
     ]
    }
   ],
   "source": [
    "print('True positive: ', tp)\n",
    "print('True negative: ', tn)\n",
    "print('False positive: ', fp)\n",
    "print('False negative: ', fn)\n",
    "print('Positive pixels: ', y_te[y_te == 1].size)\n",
    "print('Negative pixels: ', y_te[y_te == 0].size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEWCAYAAABLzQ1kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAtQklEQVR4nO3de5yVVd3//9ebGc7nAUEEzxKGlqTmIctUSrQyra9+b9I7qdu+HrLsLsuvdt8F6U0/7VtqZlqU5Dkx81TmAQ/carcnPKSiEpMYoIjCAIIchJnP749rbdiMe/bsgRnmYub9fDyux97XutZa19qMfvba61rXuhQRmJlZ/nRp7waYmVlpDtBmZjnlAG1mllMO0GZmOeUAbWaWUw7QZmY55QBtW0RST0l/krRc0h+2oJ6TJN3Xmm1rL5I+IWl2e7fDtn3yPOjOQdKJwHeAPYEVwHPA5Ih4dAvr/TLwTeBjEbF+S9uZd5ICGBkRte3dFuv43IPuBCR9B7gU+DEwFNgJuAI4thWq3xn4e2cIzpWQVN3ebbAOJCK8deAN6A+sBE4ok6c7WQB/I22XAt3TscOABcDZwFvAQuCr6diPgPeAdekcpwCTgOuL6t4FCKA67X8FeJWsFz8XOKko/dGich8DngKWp9ePFR2bAVwA/DXVcx8wuInPVmj/OUXtPw74DPB3oA74flH+A4DHgGUp7+VAt3Ts4fRZ3k2f91+K6v+/wJvAdYW0VGb3dI590/4OwGLgsPb+b8Nb/jf3oDu+g4EewG1l8vwHcBAwBtiHLEj9Z9Hx7ckC/XCyIPxLSQMjYiJZr3xaRPSJiKvKNURSb+Ay4OiI6EsWhJ8rka8GuCvlHQRcDNwlaVBRthOBrwJDgG7Ad8ucenuyf4PhwA+B3wD/CuwHfAL4oaTdUt564NvAYLJ/u7HA1wEi4tCUZ5/0eacV1V9D9mvi1OITR8Q/yIL3DZJ6Ab8Dro6IGWXaawZ4iKMzGAQsjvJDECcB50fEWxHxNlnP+MtFx9el4+si4i9kvcdRm9meBmBvST0jYmFEzCqR57PAnIi4LiLWR8TvgVeAY4ry/C4i/h4Rq4Gbyb5cmrKObLx9HXATWfD9eUSsSOefBXwYICKejojH03lfA34NfLKCzzQxItam9mwiIn4DzAGeAIaRfSGaNcsBuuNbAgxuZmx0B+CfRfv/TGkb6mgU4FcBfVrakIh4l2xY4HRgoaS7JO1ZQXsKbRpetP9mC9qzJCLq0/tCAF1UdHx1obykD0j6s6Q3Jb1D9gthcJm6Ad6OiDXN5PkNsDfwi4hY20xeM8ABujN4DFhDNu7alDfIfp4X7JTSNse7QK+i/e2LD0bEvRHxabKe5Ctkgau59hTa9PpmtqklriRr18iI6Ad8H1AzZcpOhZLUh2xc/ypgUhrCMWuWA3QHFxHLycZdfynpOEm9JHWVdLSkn6Rsvwf+U9J2kgan/Ndv5imfAw6VtJOk/sB5hQOShkr6fBqLXks2VFJfoo6/AB+QdKKkakn/AowG/ryZbWqJvsA7wMrUuz+j0fFFwG7vK1Xez4GnI+JrZGPrv9riVlqn4ADdCUTExWRzoP8TeBuYD3wDuD1l+S9gJvA88ALwTErbnHNNB6alup5m06DahWw2yBtkMxs+SboA16iOJcDnUt4lZDMwPhcRizenTS30XbILkCvIevfTGh2fBFwjaZmk/91cZZKOBY4iG9aB7O+wr6STWq3F1mH5RhUzs5xyD9rMLKccoM3McsoB2swspxygzcxyqlMv7DK4pip22bFrezfDWuDvz/dqPpPlygqWLo6I7Ta3/LjDe8eSulKzMd/v6efX3hsRR23uufKmUwfoXXbsypP37tjezbAWGLfDmPZugrXQ/XFL47tCW2RJXT1P3rtTRXmrhs1p7q7PbUqnDtBmln8BNNDQ3s1oFw7QZpZrQbAuKhvi6GgcoM0s99yDNjPLoSCo76R3PDtAm1nuNZRfMLDD8jxoM8u1AOqJirbmSBol6bmi7R1J/y6pRtJ0SXPS68CiMudJqpU0W9K4ovT9JL2Qjl0mSSm9u6RpKf0JSbsUlZmQzjFH0oTm2usAbWa510BUtDUnImZHxJiIGEP2yLNVZI+DOxd4ICJGAg+kfSSNBsYDe5GtSniFpKpU3ZVkjzgbmbbC/OtTgKURsQdwCXBRqqsGmAgcSPZYuYnFXwSlOECbWa4FsC6ioq2FxgL/iIh/kj3h/pqUfg0bH3BxLHBTepzZXKAWOEDSMKBfRDwW2ZKg1zYqU6jrFmBs6l2PA6ZHRF1ELAWmszGol+QxaDPLtahw+CIZLGlm0f6UiJjSRN7xZA+rABgaEQsBImKhpCEpfTjweFGZBSltXXrfOL1QZn6qa72k5WTPBt2QXqJMSQ7QZpZvAfWVd44XR8T+zWWS1A34PEVP/Gkqa+kWNZm+uWVK8hCHmeVadidhZVsLHA08ExGFhwcvSsMWpNe3UvoCoHg9iBFkTwRakN43Tt+kTHpYc3+yJwg1VVeTHKDNLOdEfYVbC3yJjcMbAHcChVkVE4A7itLHp5kZu5JdDHwyDYeskHRQGl8+uVGZQl3HAw+mcep7gSMlDUwXB49MaU3yEIeZ5Vp2kbBFwbcsSb2ATwOnFSVfCNws6RRgHnACQETMknQz8BKwHjgzYsN952cAVwM9gbvTBtnT26+TVEvWcx6f6qqTdAHwVMp3fkTUlWurA7SZ5Vo2D7r1AnRErCK7aFectoRsVkep/JOBySXSZwJ7l0hfQwrwJY5NBaZW2lYHaDPLvYZW7EFvSxygzSzXWrsHvS1xgDazXAtEfSedz+AAbWa55yEOM7McCsR7UdV8xg7IAdrMci27UcVDHGZmueSLhGZmORQh6sM9aDOzXGpwD9rMLH+yi4SdM1R1zk9tZtsMXyQ0M8uxes+DNjPLH99JaGaWYw2exWFmlj/ZYkkO0GZmuROIdb7V28wsfyLwjSpmZvkk36hiZpZHgXvQZma55YuEZmY5FKjTLtjfOb+WzGybEcC6qK5oq4SkAZJukfSKpJclHSypRtJ0SXPS68Ci/OdJqpU0W9K4ovT9JL2Qjl0mSSm9u6RpKf0JSbsUlZmQzjFH0oTm2uoAbWY5J+or3Cr0c+CeiNgT2Ad4GTgXeCAiRgIPpH0kjQbGA3sBRwFXSCrM+bsSOBUYmbajUvopwNKI2AO4BLgo1VUDTAQOBA4AJhZ/EZTiAG1muRZkdxJWsjVHUj/gUOAqgIh4LyKWAccC16Rs1wDHpffHAjdFxNqImAvUAgdIGgb0i4jHIiKAaxuVKdR1CzA29a7HAdMjoi4ilgLT2RjUS3KANrPca0EPerCkmUXbqY2q2g14G/idpGcl/VZSb2BoRCwESK9DUv7hwPyi8gtS2vD0vnH6JmUiYj2wHBhUpq4m+SKhmeVahFqyFsfiiNi/zPFqYF/gmxHxhKSfk4YzmlBq3CTKpG9umZLcgzazXMsuElZVtFVgAbAgIp5I+7eQBexFadiC9PpWUf4di8qPAN5I6SNKpG9SRlI10B+oK1NXkxygzSznsmcSVrI1JyLeBOZLGpWSxgIvAXcChVkVE4A70vs7gfFpZsauZBcDn0zDICskHZTGl09uVKZQ1/HAg2mc+l7gSEkD08XBI1NakzzEYWa5ll0kbNV50N8EbpDUDXgV+CpZZ/VmSacA84ATACJilqSbyYL4euDMiKhP9ZwBXA30BO5OG2QXIK+TVEvWcx6f6qqTdAHwVMp3fkTUlWuoA7SZ5V5r3kkYEc8BpcapxzaRfzIwuUT6TGDvEulrSAG+xLGpwNRK2+oAbWa51pnvJHSANrPc80NjzcxyKALWNThAm5nlTjbE4QBtZpZLLVhno0NxgM6R+bXd+fHpu2zYf3NeN778vTdZsbSKx+7tjwQDBq/ju5fOY9D263nl2V78/HvZvPcAvnz2mxxy9HIAfnfh9tz/hxpWLq/ijtoXNtT5q4k78Le/9gVg7RqxbHFXbn3lBZ77ax9+PXHjXafz/9Gd71/xTz6W6rOW6dq9gZ/dWkvXbkFVdfDIXQO47qfbs9vo1XzzwgX07N3AogXduOjMnVi1soqq6uDbP53PHh9aTVV1cP8fBjLt8qGb1Dnp6rkM2+k9TjtiVBNn7ZjaYJrdNqPNArSkAC6OiLPT/neBPhExqcLyk4CVEfHTFpzzL8CJEbFM0sqI6NPylrefHfdYy5X3zwagvh5O2ncvDjl6GX361zPhnDcBuP23g7n+ku351kUL2GXUai6/ZzZV1bBkUTVnfGoUB316OVXVcNCn3+HzX13Mvx3ywU3OcfqPNt64dMdVg6l9sScAYw5ZueHc7yyt4quHfJB9P/nO1vjYHdK6teKcE3Znzaos+F58ey1PPdiXr//X6/zm/B144fE+HDl+Ccef8RbX/r9hHHrMMrp2D04fO4ruPRuYMuMVZtw+kEULugFwyNHLWPNu5/yZTyce4mjLT70W+KKkwW14jk1ExGfSylTbvOce6cuwndcydMQ6evdt2JC+ZnUXlDoTPXoFVekrdt3ajekAH9xvFYOGri97joduH8hhxy19X/qjdw3go4e/Q49eZZcJsLLEmlXZrcfVXYOqrkEEjNh9LS883huAZx/uy8c/m/1CiYAevRroUhV069HA+vfEqpXZ/549etXzxdPe5sZLh5Y+VSfQkJ5L2NzW0bRlgF4PTAG+3fiApJ0lPSDp+fS6UxN1jJY0Q9Krks4qKn+7pKclzSperUrSa1vzC6EtzbhjAIcdt2zD/u8u3J6T9hvNg7cO5OTvLdyQ/sozvfg/h43itCNGcdZFCzYE7OYsWtCVRfO7MebjK5s9t22eLl2CK6bPZtrzs3j24T7MfrY3/5zdg4PHZb9MPvG55Wy3wzoAHvnzANas6sLvn5vF9U+9zC2/GsKKZdkfc8I5b/LHXw1h7erO2YvMZnFUVbR1NG39F/8lcJKk/o3SLweujYgPAzcAlzVRfk+yNVQLi1t3Ten/FhH7kd0NdJakQZU2SNKphaUI315S33yBdrDuPfH4ff059JhlG9K+eu6b3PD0SxzxxaXcOXW7Del77ruK38yYzS/u/js3/WII762prBcx4/aBfPyzy6hq9N/0kkXVvPZyT/Y/zMMbW6qhQXz906M4ab/RjBqzip1Hrebi7+zIMV9ZzOX3/J2efepZ/1729xr1kVU01MOJH9mLkw/ck/91+ttsv9NadttrNTvs+h7/c0/j/4U6j8KNKpVsHU2bBuiIeIdsIeuzGh06GLgxvb8O+HgTVdyVFspeTLa6VOE33lmS/gY8TrY61MgWtGlKROwfEftvNyif37hPPdiXPT60ioHbvX+I4vAvLOXRv7z/f9adRq6lR68GXpvdo6Jz/PcdA0oObzz8pwF87OhlVHctUcg2y7vvVPG3x/rw0cNXML+2B9//0u5846gPMOP2gSz8ZzbGfPgXljLzob7UrxfLl3Tlpad68YF9VjN6v3cZ+aFVXPPES/zs9lqG77aWn9xS286faOvzEEfbuZTsETC9y+RparBzbdH7eqBa0mHAp4CDI2If4Fmgsqi0jZhx+8BNhhhef7XbhveP39ufHffI/lnenNeN+hTDFy3oyoJ/9GDoiPearX9+bXdWLq9m9P6rmj23bZ7+Nevp3S/7hdatRwP7fmIl82t70H9QNqQhBSd+axF/vi778ff264XhpqB7z3r23HcV82u78+drB3Pivnsx4cDRnH3cHrz+anfOOX6P9vpY7aIwi6Mz9qDbfJpdWsHpZrIgXVgk5H/IVni6DjgJeLQFVfYne97XKkl7Age1Znvb25pV4plH+vKtn2x88MJVP96BBf/oTpcuMGT4e5x1UfYghxef7M20y3elujob7/zmjxfQf1AWFH57wTAeun0ga1d34aT9RnPUl+r48nezmSAzbh/IJ49duslFRYA353fj7Te68uGD3z8ubS1TM3Qd3/35PLp0gS5d4OE/9eeJ+/tx3Clvc8xXFgPw17v7c99NNQDc+btBnH3JfKY8NBsE902rYe7LPdvzI+RKZ53FoWyZ0jaouGiam6ShwFzgJxExKT3ldiowmOzxM1+NiHmNyk+iaJqdpBeBzwELgdvJHhUzG9gOmBQRMyS9BuwfEYsrmWa3/z494sl7dyyXxXJm3A5j2rsJ1kL3xy1PN/OUk7IG7jkkjph6fEV5bz3kyi06V960WQ+6ODhGxCKgV9H+a8ARzZSf1Gi/eFm/o5sos0up85vZtq0jDl9UwncSmlmu+U5CM7Mcc4A2M8shL9hvZpZjHXGOcyUcoM0s1yJgvRfsNzPLp846xNE5v5bMbJvR2mtxpEXVXpD0nKSZKa1G0nRJc9LrwKL850mqlTRb0rii9P1SPbWSLpOyW78kdZc0LaU/ke77KJSZkM4xR9KE5trqAG1muRehirYWODwixhTd1HIu8EBEjAQeSPtIGk121/NewFHAFZIKi/hcCZxKthbQyHQcsruml0bEHsAlwEWprhpgInAgGxeA2/BFUIoDtJnl3lZYLOlY4Jr0/hrguKL0m9KibXOBWuAAScOAfhHxWGS3Y1/bqEyhrluAsal3PQ6YHhF1EbEUmM7GoF6SA7SZ5VpEixZLGlxYTjhtp5aqErgvrSlfOD40IhZm54uFwJCUPhyYX1R2QUobnt43Tt+kTESsB5YDg8rU1SRfJDSznBP1lc/iWFzBWhyHRMQbkoYA0yW9Uvbk7xdl0je3TEnuQZtZ7rXmGHREvJFe3wJuIxsPXpSGLUivb6XsC8jWnC8YAbyR0keUSN+kjKRqshU468rU1SQHaDPLtdZcD1pSb0l9C++BI4EXgTuBwqyKCcAd6f2dwPg0M2NXsouBT6ZhkBWSDkrjyyc3KlOo63jgwTROfS9wpKSB6eLgkSmtSR7iMLN8i2wcupUMBW5LM+KqgRsj4h5JTwE3SzoFmAecABARs9J69i+RPWf1zIgoPCvvDOBqoCdwd9oArgKuk1RL1nMen+qqk3QB8FTKd35E1JVrrAO0meVea93qHRGvAvuUSF8CjG2izGRgcon0mcDeJdLXkAJ8iWNT2fjgkmY5QJtZrkXLLhJ2KA7QZpZ7bfTgp9xzgDaz3GvhXYIdhgO0meVahAO0mVluddbV7BygzSz3PAZtZpZDgWjwLA4zs3zqpB1oB2gzyzlfJDQzy7FO2oV2gDaz3HMPuhFJv6DM91ZEnNUmLTIzKxJAQ4MDdGMzt1orzMyaEoB70JuKiGuK9yX1joh3275JZmab6qzzoJudXCjpYEkvAS+n/X0kXdHmLTMzK4gKtw6mktnfl5I9jXYJQET8DTi0DdtkZlakssdddcQLiRXN4oiI+ekJBAX1TeU1M2t1HbB3XIlKAvR8SR8DQlI34CzScIeZWZsLiE46i6OSIY7TgTOB4cDrwJi0b2a2lajCrWNptgcdEYuBk7ZCW8zMSuukQxyVzOLYTdKfJL0t6S1Jd0jabWs0zswMaPVZHJKqJD0r6c9pv0bSdElz0uvAorznSaqVNFvSuKL0/SS9kI5dpnShTlJ3SdNS+hOSdikqMyGdY46kCc21s5IhjhuBm4FhwA7AH4DfV/jvYGa2ZQo3qlSyVe5bbHot7VzggYgYCTyQ9pE0GhgP7AUcBVwhqSqVuRI4FRiZtqNS+inA0ojYA7gEuCjVVQNMBA4EDgAmFn8RlFJJgFZEXBcR69N2PZ32B4eZtYfssVfNb5WQNAL4LPDbouRjgcLNedcAxxWl3xQRayNiLlALHCBpGNAvIh6LiACubVSmUNctwNjUux4HTI+IuohYCkxnY1AvqdxaHDXp7UOSzgVuIgvM/wLcVa5SM7NWVfksjsGSipepmBIRUxrluRQ4B+hblDY0IhYCRMRCSUNS+nDg8aJ8C1LauvS+cXqhzPxU13pJy4FBxeklypRU7iLh02QBufAvc1rRsQAuKFexmVlrUeW/2RdHxP5N1iN9DngrIp6WdFglpy6RFmXSN7dMSeXW4ti1XEEzs62idW/jPgT4vKTPAD2AfpKuBxZJGpZ6z8OAt1L+BcCOReVHAG+k9BEl0ovLLJBUDfQH6lL6YY3KzCjX2Ioe9CVpb0n/W9LJha2ScmZmW67CC4QVXCSMiPMiYkRE7EJ28e/BiPhX4E6gMKtiAnBHen8nMD7NzNiV7GLgk2k4ZIWkg9L48smNyhTqOj6dI4B7gSMlDUwXB49MaU1qdh60pIlkUX808BfgaOBRskFxM7O21/bTEi4EbpZ0CjAPOAEgImZJuhl4CVgPnBkRhaUuzgCuBnoCd6cN4CrgOkm1ZD3n8amuOkkXAE+lfOdHRF25RlVyq/fxwD7AsxHxVUlD2fTqp5lZ22po/SojYgZpiCEilgBjm8g3GZhcIn0msHeJ9DWkAF/i2FRgaqVtrCRAr46IBknrJfUjG5vxjSpmtnV4wf6yZkoaAPyGbGbHSuDJtmyUmVmxFszi6FAqWYvj6+ntryTdQzY5+/m2bZaZWREH6E1J2rfcsYh4pm2aZGZmUL4H/bMyxwI4opXbstX9/flejNthTHs3w8ya4SGORiLi8K3ZEDOzkoKW3OrdoVT0yCszs3blHrSZWT55iMPMLK86aYCu5IkqkvSvkn6Y9neSdEDbN83MLGnlJ6psKypZLOkK4GDgS2l/BfDLNmuRmVkRReVbR1PJEMeBEbGvpGcBImKppG5t3C4zs408i6NJ69IzuAJA0na0ydIlZmaldcTecSUqGeK4DLgNGCJpMtlSoz9u01aZmRXrpGPQlazFcYOkp8mW4hNwXES83EwxM7PW0UHHlytRyYL9OwGrgD8Vp0XEvLZsmJnZBg7QTbqLjQ887AHsCswG9mrDdpmZbaBOetWrkiGODxXvp1XuTmsiu5mZtZIW30kYEc9I+mhbNMbMrCQPcZQm6TtFu12AfYG326xFZmbFfJGwrL5F79eTjUn/sW2aY2ZWQicN0GXnQacbVPpExI/SNjkibkhPrTUz2zpaaR60pB6SnpT0N0mzJP0opddImi5pTnodWFTmPEm1kmZLGleUvp+kF9KxyyQppXeXNC2lPyFpl6IyE9I55kia0Fx7mwzQkqojop5sSMPMrF2IbBZHJVsF1gJHRMQ+wBjgKEkHAecCD0TESOCBtI+k0cB4sllrRwFXpI4rwJXAqcDItB2V0k8BlkbEHsAlwEWprhpgInAgcAAwsfiLoJRyPejCk7ufk3SnpC9L+mJhq+ifwsxsS7XiYkmRWZl2u6YtgGOBa1L6NcBx6f2xwE0RsTYi5gK1wAGShpE9QPuxiAjg2kZlCnXdAoxNvetxwPSIqIuIpcB0Ngb1kioZg64BlpA9g7AwHzqAWysoa2a25Sofgx4saWbR/pSImFKcIfWAnwb2AH4ZEU9IGhoRCwEiYqGkISn7cODxouILUtq69L5xeqHM/FTXeknLgUHF6SXKlFQuQA9JMzheZGNgLuikQ/Zm1i4qjziLI2L/slVlQ7djJA0AbpO0d5nspZbRaxwPG7dyc8qUVG6Iowrok7a+Re8Lm5nZVtEW60FHxDJgBtkww6I0bEF6fStlWwDsWFRsBPBGSh9RIn2TMpKqgf5AXZm6mlSuB70wIs4vV9jMbKtopd/sabnkdRGxTFJP4FNkF/HuBCYAF6bXO1KRO4EbJV0M7EB2MfDJiKiXtCJdYHwCOBn4RVGZCcBjwPHAgxERku4Fflx0YfBI4Lxy7S0XoDvnCtlmli/RqmtxDAOuSePQXYCbI+LPkh4DbpZ0CjAPOAEgImZJuhl4iew+kDPTEAnAGcDVQE/g7rQBXAVcJ6mWrOc8PtVVJ+kC4KmU7/yIqCvX2HIBemzln9nMrA21Ug86Ip4HPlIifQlNxLyImAxMLpE+E3jf+HW6T+SEJuqaCkyttL1NBujmIruZ2dbiW73NzPLKAdrMLIc66OOsKuEAbWa5JjzEYWaWWw7QZmZ55QBtZpZTDtBmZjnkJ6qYmeWYA7SZWT614q3e2xQHaDPLPQ9xmJnlkW9UMTPLMQdoM7P88Z2EZmY5pobOGaEdoM0s3zwGbWaWXx7iMDPLKwdoM7N8cg/azCyvHKDNzHKodZ/qvU3p0t4NMDMrpzAPupKt2bqkHSU9JOllSbMkfSul10iaLmlOeh1YVOY8SbWSZksaV5S+n6QX0rHLJCmld5c0LaU/IWmXojIT0jnmSJrQXHsdoM0s/yIq25q3Hjg7Ij4IHAScKWk0cC7wQESMBB5I+6Rj44G9gKOAKyRVpbquBE4FRqbtqJR+CrA0IvYALgEuSnXVABOBA4EDgInFXwSlOECbWe61Vg86IhZGxDPp/QrgZWA4cCxwTcp2DXBcen8scFNErI2IuUAtcICkYUC/iHgsIgK4tlGZQl23AGNT73ocMD0i6iJiKTCdjUG9JI9BbwO6dm/gZ7fW0rVbUFUdPHLXAK776fbsttdqzrpwAd16NFC/Xlx+3ghmP9eLw7+wlBO+/taG8rt+cA1njvsAr87qyeQbXqVmyDqqqoMXn+jD5d8fTkOD2vHTdUzb7fAe3/v5PAYOWU80wF+uH8TtV23HJz63jC+f/SY7jlzLWZ8ZyZznewFQ3bWBb/1kASM/vJpogCt/OJznH+sDwE9uqaVm6HreW5P9nc4bvxvLl3Rtt8+21bXsRpXBkmYW7U+JiCmlMqahh48ATwBDI2IhZEFc0pCUbTjweFGxBSltXXrfOL1QZn6qa72k5cCg4vQSZUrKbYCWtDIi+hTtfwXYPyK+UabM54HREXGhpEnAyoj4aZs3to2tWyvOOWF31qyqoqo6uPj2Wp56sC8nf+9Nrr94KDMf6sdHj3iHU/7zDc45fg8eum0gD92W/XLaZc/VTPrda7w6qycAk0/bmVUrq4DgB7/5J584Zhn/fUfZX1m2GerXiynn70DtC73o2buey+/5O8883JfXXunB+V/bhbMuWrBJ/qNPqgPg9LGj6D9oHZNvmMs3jx5JRBaULzpzpw3BvDNqwUXCxRGxf7P1SX2APwL/HhHvpOHjkllLpEWZ9M0tU1KHGuKIiDsj4sL2bkfrE2tWZcNe1V2Dqq6xYcitd996AHr3q6du0ft7VYcft4wZtw/YsJ8FZ6iqhupu0WmnL7W1ure6UvtCFlBXv1vF/NoeDB62jvm1PVjwjx7vy7/TB9bw7CN9AVi+pCsrl1fxgX1Wb9U255kaKtsqqkvqShacb4iIW1PyojRsQXot/ARdAOxYVHwE8EZKH1EifZMykqqB/kBdmbqatE0GaEnHpKujz0q6X9LQlP4VSZe3d/vaQpcuwRXTZzPt+Vk8+3AfZj/bm1/9cDhf+8FCrp/5Ev/nB28w9cfD3lfu0M8v46GiAA0w+cZ/MO35Waxe2YVH/jzgfWWsdQ0d8R67772aV55pugf86qyeHDxuOV2qgqE7rmXkh1ex3Q7vbTh+9iXzuWL6bE7890V0um/VoNUuEqax4KuAlyPi4qJDdwKFWRUTgDuK0senmRm7kl0MfDINh6yQdFCq8+RGZQp1HQ88mMap7wWOlDQwXRw8MqU1KbdDHEBPSc8V7deQfXCAR4GDIiIkfQ04Bzi7kkolnUp25ZUebDs/GRsaxNc/PYre/eqZeNVcdh61ms/8ax2/nrgDj/5lAIces4zvXDyfc/9l9w1lRn3kXdau7sI/Z/fcpK7/OHF3unZv4NzL5zHm4yt55uG+W/vjdBo9etXzg9++xq9+uMOGXy+l3HtTDTuNXMPl9/ydtxZ046WZvamvT8Mb39iZJW92pWfvrK5PHd+V+2+p2VofIRda8U7CQ4AvAy8UxZfvAxcCN0s6BZgHnAAQEbMk3Qy8RDYD5MyIqE/lzgCuBnoCd6cNsi+A6yTVkvWcx6e66iRdADyV8p0fEXXlGpvnAL06IsYUdgpj0Gl3BDAt/RTpBsyttNJ0wWAKQD/VbHNdkXffqeJvj/Xho4ev4NMn1HHlD3YA4OE/9efffzp/k7yHHbvp8EaxdWu78Nh9/Th43HIH6DZSVR384Lev8eCtA/nr3QPK5m2oF7+etPF60SV3zuH1V7sDsOTNbOhq9btVPHTbQEZ9ZFWnC9Ct9aMhIh6l9FgwwNgmykwGJpdInwnsXSJ9DSnAlzg2FZhaaXu3ySEO4BfA5RHxIeA04P2Deh1I/5r19O6XfWl369HAvp9YyfzaHixZ1JUPH/wuAGM+vpI35nbfUEYKPvG55cy4Y8CGtB696qkZsg6ALlXBAWPfYX5th/6na0fBd342n/lzenDrlO2azd29ZwPde2Z/430PXUH9ejFvTg+6VAX9atYDWcA/8FPv8Nornetv1po3qmxr8tyDLqc/8Hp63+zdONu6mqHr+O7P59GlC3TpkvWWn7i/Hyvf6cIZ579BVVXw3touXPq9jdcsPnTQuyxe2JU3520M2j16NTDp6rnZdL2q4Lm/9uHP1w5qj4/U4e11wLt86oSlvPpSD66YPhuA3/1/w+jaLfj6f71O/0HrueC6ufxjVg/+48TdGTBoPZN//yrRkPWYf/LNnQDo2q2BH9/4KlXV2d/smUf6cvcNnexvFtFpF+xXVHb3zVZXbpqdpGPJ7tB5nWyO4kcj4rBGeSbRzDS7fqqJA1XyV42ZtZL745anK5n61pS+A0bERw79VkV5H/nTOVt0rrzJbQ+6ODin/avJBuSJiDvYeMW0qTyT2raFZra1dMThi0rkNkCbmQHZBcJOOsThAG1m+dc547MDtJnln4c4zMxyqrPO4nCANrN8a9lqdh2KA7SZ5Vp2o0rnjNAO0GaWf530mYQO0GaWe+5Bm5nlkcegzczyqvOuxeEAbWb55yEOM7McihY9k7BDcYA2s/xzD9rMLKc6Z3x2gDaz/FND5xzjcIA2s3wLfKOKmVkeiei0N6psqw+NNbPOJKKyrRmSpkp6S9KLRWk1kqZLmpNeBxYdO09SraTZksYVpe8n6YV07DJJSundJU1L6U9I2qWozIR0jjmSKnqWqgO0meVfKwVoskfiHdUo7VzggYgYCTyQ9pE0GhgP7JXKXCGpKpW5EjgVGJm2Qp2nAEsjYg+y56ZelOqqASYCBwIHABOLvwia4gBtZvlWGIOuZGuuqoiHgbpGyccC16T31wDHFaXfFBFrI2IuUAscIGkY0C8iHovsqdvXNipTqOsWYGzqXY8DpkdEXUQsBabz/i+K9/EYtJnlXgtmcQyWNLNof0pETGmmzNCIWAgQEQslDUnpw4HHi/ItSGnr0vvG6YUy81Nd6yUtBwYVp5co0yQHaDPLuYqHLwAWR8T+rXRilW5Mk+mbW6ZJHuIws3wLWnMMupRFadiC9PpWSl8A7FiUbwTwRkofUSJ9kzKSqoH+ZEMqTdVVlgO0meVfK41BN+FOoDCrYgJwR1H6+DQzY1eyi4FPpuGQFZIOSuPLJzcqU6jreODBNE59L3CkpIHp4uCRKa0sD3GYWe611jxoSb8HDiMbq15ANrPiQuBmSacA84ATACJilqSbgZeA9cCZEVGfqjqDbEZIT+DutAFcBVwnqZas5zw+1VUn6QLgqZTv/IhofLHyfRygzSz/WilAR8SXmjg0ton8k4HJJdJnAnuXSF9DCvAljk0FplbcWBygzSzvIqC+c97r7QBtZvnXSW/1doA2s/xzgDYzy6EA/ExCM7M8CgiPQZuZ5U/gi4RmZrnlMWgzs5xygDYzy6MtWmdjm+YAbWb5FoAfGmtmllPuQZuZ5ZFv9TYzy6eA8DxoM7Oc8p2EZmY55TFoM7McivAsDjOz3HIP2swsj4Kor28+WwfkAG1m+eblRs3McszT7MzM8ieAcA/azCyHwgv2m5nlVme9SKjopNNXACS9DfyzvdvRBgYDi9u7EdYiHflvtnNEbLe5hSXdQ/bvU4nFEXHU5p4rbzp1gO6oJM2MiP3bux1WOf/NrJQu7d0AMzMrzQHazCynHKA7pint3QBrMf/N7H08Bm1mllPuQZuZ5ZQDtJlZTjlA55CkkPSzov3vSprUgvKTJH23hef8i6QB6f3KlpS18hr/e0r6iqTLmynzeUnnpvct/ntax+AAnU9rgS9KqnRy/haLiM9ExLKtdT4rLyLujIgL27sd1r4coPNpPdlV/W83PiBpZ0kPSHo+ve7URB2jJc2Q9Kqks4rK3y7paUmzJJ1alP7a1vxCsIykYyQ9IelZSfdLGprSm+1lW8fnAJ1fvwROktS/UfrlwLUR8WHgBuCyJsrvCYwDDgAmSuqa0v8tIvYD9gfOkjSo9ZtujfSU9FxhA84vOvYocFBEfAS4CTinPRpo+eTFknIqIt6RdC1wFrC66NDBwBfT++uAnzRRxV0RsRZYK+ktYCiwgCwofyHl2REYCSxp7fbbJlZHxJjCjqSvkH1BAowApkkaBnQD5m711lluuQedb5cCpwC9y+RpaiL72qL39UC1pMOATwEHR8Q+wLNAjy1upW2JXwCXR8SHgNPw38OKOEDnWETUATeTBemC/wHGp/cnkf1ErlR/YGlErJK0J3BQqzTUtkR/4PX0fkJ7NsTyxwE6/37GpkstngV8VdLzwJeBb7WgrnvIetLPAxcAj7daK21zTQL+IOkROu5yo7aZfKu3mVlOuQdtZpZTDtBmZjnlAG1mllMO0GZmOeUAbWaWUw7QVpak+nSL8ouS/iCp1xbUdbWk49P730oaXSbvYZI+thnnKLmmSCVrjbR0FT+vMmdtzQHamrM6IsZExN7Ae8DpxQclVW1OpRHxtYh4qUyWw4AWB2izjsQB2lriEWCP1Lt9SNKNwAuSqiT9P0lPpVX2TgNQ5nJJL0m6CxhSqCittLd/en+UpGck/S2t0LcL2RfBt1Pv/ROStpP0x3SOpyQdksoOknRfWg3u14Ca+xBNreiXjv0steUBSdultN0l3ZPKPJLuwjRrc14sySoiqRo4muxuRMhWyds7IuamILc8Ij4qqTvwV0n3AR8BRgEfIlus6SVgaqN6twN+Axya6qqJiDpJvwJWRsRPU74bgUsi4tG0xOq9wAeBicCjEXG+pM8CmwTcJvxbOkdP4ClJf4yIJWRrnjwTEWdL+mGq+xtkS7+eHhFzJB0IXAEcsRn/jGYt4gBtzemZlsiErAd9FdnQw5MRUVh57Ujgw4XxZbL1JUYChwK/j4h64A1JD5ao/yDg4UJdaf2RUj5FtsZ1Yb+fpL7pHF9MZe+StLSCz9TUin4NwLSUfj1wq6Q+6fP+oejc3Ss4h9kWc4C25myyVCZAClTvFicB34yIexvl+wxNr7ZXXLaS9Qa6kK3CV7z0aqEtFa9X0GhFv1WSZtD0CnKRzrus8b+B2dbgMWhrDfcCZxQeCiDpA5J6Aw8D49MY9TDg8BJlHwM+KWnXVLYmpa8A+hblu49suIGUb0x6+zDZqn5IOhoY2Exby63o1wUo/Ao4kWzo5B1grqQT0jkkaZ9mzmHWKhygrTX8lmx8+RlJLwK/Jvt1dhswB3gBuBL478YFI+JtsnHjWyX9jY1DDH8CvlC4SEi2it/+6SLkS2ycTfIj4FBJz5ANtcxrpq3lVvR7F9hL0tNkY8yFJ5+cBJyS2jcLOLaCfxOzLebV7MzMcso9aDOznHKANjPLKQdoM7OccoA2M8spB2gzs5xygDYzyykHaDOznPr/AeJue08993yEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "disp = plot_confusion_matrix(grid_RF_cv, X_te, y_te, labels=[0, 1], display_labels=['No hail', 'Hail'], normalize=None)\n",
    "disp.ax_.set_title('Confusion matrix')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of Detection (POD): 36.17336849044166\n",
      "False Alarm Rate (FAR): 0.5404518144611558\n",
      "Positive Predictive Value (PPV): 0.35523547499595404\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": [
    "When using the balanced accuracy as a target score we do improve significantly the probability of detection without penalizing too much the false alarm rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(model_name, 'wb') as bfile: \n",
    "    s = pickle.dump(grid_RF_cv.best_estimator_, bfile)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(model_name, 'rb') as bfile:\n",
    "    rf_model_best = pickle.load(bfile)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get features importance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_importance = rf_model_best['RF'].feature_importances_\n",
    "std = np.std([tree.feature_importances_ for tree in rf_model_best['RF'].estimators_], axis=0)\n",
    "inds = np.argsort(f_importance)[::-1]\n",
    "f_importance_ordered = f_importance[inds]\n",
    "features_ordered = np.array(features)[inds]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAEICAYAAAAwft9dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgaElEQVR4nO3de7xVZb3v8c8XVPCKF6qNoFKGpkWikrdunqPba4me3Iq7Usqku2miYrr35pglqCWZ3cg8XvJFXvJCR1OrrVaK5EJhgZgKggpyRLEIQ0nhd/4Yz9LhbK75rMtcay7W+r5fr/Vac1yeZ/zGYKz15RnzmaCIwMzMzFrXr9EFmJmZ9XQOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmlndSPqGpCsaXYdZvcmfszTrGSQtAd4BrCut3iUinutkn5+LiN92rroNj6RJwLsj4lONrsU2fB5ZmvUsH4+ILUpfHQ7KepC0USOP31Ebat3WczkszXo4SYMk/UzScknLJF0gqX/atrOk/5a0UtKLkq6TtHXadi2wI/ArSS9LOkvSgZKWVvS/RNLB6fUkSTdJ+rmkvwHjah2/Sq2TJP08vR4uKSR9RtKzkv4i6QuSPiCpWdJfJV1eajtO0v2Svi9plaQ/SzqotH17STMkvSRpoaRTKo5brvsLwDeA49O5z037fUbSY5JWS3pK0udLfRwoaamkMyStSOf7mdL2TSV9R9LTqb4/Sto0bdtP0gPpnOZKOrADf9TWgzkszXq+q4HXgXcDewKHAJ9L2wRcCGwP7AbsAEwCiIhPA8/w5mj1ojYebwxwE7A1cF3m+G2xLzACOB6YCpwLHAy8FzhO0kcr9n0KGAz8F3CzpG3TtunA0nSuxwLfLodpRd0/A74NXJ/OfY+0zwrgY8BWwGeASyXtVerjX4BBwFDgZOAHkrZJ2y4B9gYOALYFzgLWSxoK3A5ckNZPAH4p6W3tuEbWwzkszXqWW9Po5K+SbpX0DuBw4LSI+HtErAAuBcYCRMTCiPhNRKyNiBeA7wIfbb37NpkZEbdGxHqKUGn1+G30zYh4NSLuBv4OTI+IFRGxDPgDRQC3WAFMjYjXIuJ64HHgSEk7AB8Czk59zQGuAD5dre6IeKVaIRFxe0QsisJ9wN3Ah0u7vAacn45/B/AysKukfsBnga9FxLKIWBcRD0TEWuBTwB0RcUc69m+AJuCIdlwj6+H8XN+sZzm6PBlH0j7AxsBySS2r+wHPpu1vBy6j+IW/Zdr2l07W8Gzp9U61jt9Gz5dev1JleYvS8rJ466zDpylGktsDL0XE6opto1upuypJh1OMWHehOI/NgHmlXVZGxOul5TWpvsHAQGBRlW53Av5N0sdL6zYG7snVYxsOh6VZz/YssBYYXPFLvMWFQADvj4iVko4GLi9tr5zu/neKgAAgvfdY+biw3CZ3/HobKkmlwNwRmAE8B2wractSYO4ILCu1rTzXtyxLGgD8EjgRuC0iXpN0K8Wj7JwXgVeBnYG5FdueBa6NiFP+qZX1Gn4Ma9aDRcRyikeF35G0laR+aVJPy6PWLSkeFf41vXd2ZkUXzwPvKi0/AQyUdKSkjYHzgAGdOH69vR04VdLGkv6N4n3YOyLiWeAB4EJJAyW9n+I9xetq9PU8MDw9QgXYhOJcXwBeT6PMQ9pSVHokfSXw3TTRqL+k/VMA/xz4uKRD0/qBabLQsPafvvVUDkuznu9Eil/0Cygesd4EDEnb/jewF7CKYpLJzRVtLwTOS++BToiIVcCXKN7vW0Yx0lxKbbWOX2+zKCYDvQh8Czg2IlambScAwylGmbcA/5XeH2zNjen7SkkPpxHpqcANFOfx7xSj1raaQPHI9iHgJWAK0C8F+RiK2bcvUIw0z8S/X3sV/6MEZtYjSBpH8Q8ofKjRtZhV8t98zMzMMhyWZmZmGX4Ma2ZmluGRpZmZWYY/Z9kLDR48OIYPH97oMszMNiizZ89+MSKq/jOFDsteaPjw4TQ1NTW6DDOzDYqkp1vb5sewZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwy/I8S9ELzlq1i+MTbG12GWZ+wZPKRjS7BuoFHlmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZRo8KS0kvVyyPk3R5ej1J0jJJcyQtkHRCWr9Y0q4V7aZKOqv7Km8bSaMkHdGJ9kdL2r2eNZmZWV6PCss2uDQiRgFjgJ9I2hj4BTC2ZQdJ/YBjgeu7ogBJ/TvRfBTQ4bAEjgYclmZm3WxDC0sAIuJJYA2wDTCdUlgCHwGWRMTT1dqm0erNku6U9KSki0rbTpA0T9J8SVNK61+WdL6kWcD+aXmKpNmSfitpH0n3SnpK0lGtHHcT4Hzg+DQ6Pl7S5pKulPSQpEckjUn7XibpP9PrQyX9XtIBwFHAxan9zhX9j5fUJKlp3ZpV7b6mZmbWup72X3RtKmlOaXlbYEblTpL2Ap6MiBXACknrJe0REXMpgnN65jijgD2BtcDjkr4PrAOmAHsDfwHulnR0RNwKbA7Mj4iWANscuDcizpZ0C3AB8K8Uo76rq9UcEf9IATg6Ir6S+vk28N8R8VlJWwN/kvRbYCLwkKQ/AJcBR0TEIkkzgP8bETdV6X8aMA1gwJARkTl/MzNrh54Wlq+kx6xAMQoERpe2ny7pFOBdwGGl9dOBsZIepXhE+5+Z4/wuIlalYywAdgK2owjAF9L66yhGqbdSBOkvS+3/AdyZXs8D1kbEa5LmAcPbeK4AhwBHSZqQlgcCO0bEY+k8fw+cHhGL2tGnmZnVWU8Ly5xLI+ISSf8LuEbSzhHxKkVY3g3cBzSnEWcta0uv11FcB9XY/9WIWFdafi0iWkZv61v6i4j1ktpzTQV8IiIer7JtJLAS2L4d/ZmZWRfYUN+zvBloAk5Ky4sogmUy+UewrZkFfFTS4DSJ5wSK8K2n1cCWpeW7gK9KEoCkPdP3nYAzKB4VHy5p31bam5lZN9ggwzI5H/h6mv0KRUi+B7ilI51FxHLgHOAeYC7wcETcVo9CS+4Bdm+Z4AN8E9gYaJY0H/hmCs6fARMi4jngZOAKSQMpZv6emSYD7dzKMczMrM705tNE6y0GDBkRQ06a2ugyzPqEJZOPbHQJVieSZkfE6GrbNuSRpZmZWbfY0Cb4tJmkQyk+ClK2OCKO6c3HNjOz+uu1YRkRd1FMoOlTxzYzs/rzY1gzM7MMh6WZmVmGw9LMzCzDYWlmZpbRayf49GUjhw6iyZ/9MjOrG48szczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw7Nhe6F5y1YxfOLtjS7DrM/y/0TS+3hkaWZmluGwNDMzy3BYmpmZZTgszczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllZMNS0qWSTist3yXpitLydySFpF0r2k2VdFaNfs+RtFDS45IOLa3fRNI0SU9I+rOkT6T1X5e0QFKzpN9J2qmVfidJmpBeXyVpsaQ5kuZKOihzrl9JNYWkwaX1knRZ2tYsaa/SttMlPSppvqTpkgbW6H+cpO1r1ZCp7xsdbWtmZh3XlpHlA8ABAJL6AYOB95a2HwDcCYxtWZH2Oxa4vlqHknZP+78XOAz4oaT+afO5wIqI2AXYHbgvrX8EGB0R7wduAi5qQ+0AZ0bEKOA04MeZfe8HDgaerlh/ODAifY0HfpTOYyhwaqrrfUB/StehinFAh8MScFiamTVAW8LyflJYUoTbfGC1pG0kDQB2A87mrSHxEWBJRFSGTosxwC8iYm1ELAYWAvukbZ8FLgSIiPUR8WJ6fU9ErEn7PAgMa8sJlswEhtbaISIeiYglrdR7TRQeBLaWNCRt2wjYVNJGwGbAc9X6lnQsMBq4Lo10N5W0t6T7JM1OI/Yhkgal0fauqd10SadImpyOM0fSdVX6Hy+pSVLTujWr2nZFzMysTbJhGRHPAa9L2pEiNGcCs4D9KX75N0dEM7Be0h6p2Vhgeo1uhwLPlpaXAkMlbZ2WvynpYUk3SnpHlfYnA7/O1V7hMODWdrZpUbXeiFgGXAI8AywHVkXE3dU6iIibgCbgk2mk+zrwfeDYiNgbuBL4VkSsAr4CXCVpLLBNRPw0IiYCr0TEqIj4ZJX+p0XE6IgY3X+zQR08TTMzq6atE3xaRpctYTmztPxA2mc6MDaNsMYAN9boT1XWBcUobRhwf0TslY5zyVsaSp+iCOmL21j7xZKeAn4OfLuNbdpUr6RtKM71nRSPVzdP9bXFrsD7gN9ImgOcRxotR8RvgHnAD4DPdbBmMzOrk7aGZcv7liMpHsM+SDGyPIAiSKEIy+Mo3vNrjogVNfpbCuxQWh5G8fhyJbAGuCWtvxEoT6Y5mOI9zaMiYm1a9630aHJOK8c6E3g3RRhd3YZzbU+9BwOLI+KFiHgNuJk3H1nnCHg0jRRHRcTIiDgE3njPdzfgFWDbDtZsZmZ10p6R5ceAlyJiXUS8BGxNEZgzASJiEUXYTab2I1iAGRSj0AGS3kkxceZPERHAr4AD034HAQsAJO0J/IQiKN8I4og4tyVwWjtYRKwHvgf0K8+8bYcZwIlpVux+FI9bl1M8ft1P0maSlOp9rEY/q4Et0+vHgbdJ2j+d38aSWiZOnZ76OQG4UtLGaf1rpddmZtZN2hqW8yhmwT5YsW5VywScZDrwHt4cGVYVEY8CN1AE4Z3AlyNiXdp8NjBJUjPwaeCMtP5iYAvgxjSSnNHG2luOGcAFQK2Ps5wqaSnFyLG59BGZO4CnKCYi/RT4UupzFsXM3Icprkc/YFqNMq4CfpxGwf0pZgxPkTQXmAMcIGkXikevZ0TEH4DfU4yKSX03V5vgY2ZmXUdFhlhvMmDIiBhy0tRGl2HWZy2ZfGSjS7AOkDQ7IkZX2+Z/wcfMzCxjo67sPL0/OKVi9eKIOKYrj5sj6RaKGaxlZ0fEXXXq/wfABytWfy8i/k89+jczs+7VpWGZwqcuAVRPXR3WEfHlruzfzMy6lx/DmpmZZTgszczMMhyWZmZmGQ5LMzOzjC6d4GONMXLoIJr8OS8zs7rxyNLMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDI8G7YXmrdsFcMn3t7oMsz6LP+vI72PR5ZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWUafCktJL1csj5N0eXo9SdIySXMkLZB0Qlq/WNKuFe2mSjqr+yo3M7NG6lNh2QaXRsQoYAzwE0kbA78AxrbsIKkfcCxwfVcUIKl/V/RrZmYd57CsIiKeBNYA2wDTKYUl8BFgSUQ8Xa1tGq3eLOlOSU9Kuqi07QRJ8yTNlzSltP5lSedLmgXsn5anSJot6beS9pF0r6SnJB3VynHHS2qS1LRuzao6XAUzM2vR18Jy0/SYdY6kOcD51XaStBfwZESsiIhmYL2kPdLmsRQBWsso4HhgJHC8pB0kbQ9MAf5n2v4BSUen/TcH5kfEvhHxx7R8b0TsDawGLgD+FTimtZojYlpEjI6I0f03G5S7DmZm1g597f+zfCU9ZgWKUSAwurT9dEmnAO8CDiutnw6MlfQoxSPa/8wc53cRsSodYwGwE7AdRQC+kNZfRzFKvRVYB/yy1P4fwJ3p9TxgbUS8JmkeMLyN52pmZnXS10aWOZdGxK4Uo8JrJA1M66cDxwEHA80RsSLTz9rS63UUfylRjf1fjYh1peXXIiLS6/Ut/UXEevreX3DMzBrOYVlFRNwMNAEnpeVFwEpgMvlHsK2ZBXxU0uA0iecE4L46lGtmZl3MYdm684Gvp9mvUITke4BbOtJZRCwHzgHuAeYCD0fEbfUo1MzMupbefNpnvcWAISNiyElTG12GWZ+1ZPKRjS7BOkDS7IgYXW2bR5ZmZmYZnizSQZIOpfgoSNniiDimEfWYmVnXcVh2UETcBdzV6DrMzKzr+TGsmZlZhsPSzMwsw2FpZmaW4bA0MzPL8ASfXmjk0EE0+XNeZmZ145GlmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZlleDZsLzRv2SqGT7y90WWY9Wn+n0d6F48szczMMhyWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7OMXhmWkl5O34dLekXSHEkLJF0jaeMa7baTdI+klyVdXrFtb0nzJC2UdJkkpfU7pjaPSGqWdESN/odL+vdOnNeBkg7oaHszM+uYXhmWFRZFxChgJDAMOK7Gvq8C/wFMqLLtR8B4YET6OiytPw+4ISL2BMYCP6zR/3Cgw2EJHAg4LM3MullfCEsAImId8CdgaI19/h4Rf6QIzTdIGgJsFREzIyKAa4CjW5oBW6XXg4DnapQxGfhwGumeLqm/pIslPZRGpZ9Px/u6pCvT65GS5kvaHfgCcHpq/+GKGsdLapLUtG7NqjZdEzMza5s+8190SRoI7At8rQPNhwJLS8tLeTN0JwF3S/oqsDlwcI1+JgITIuJjqabxwKqI+ICkAcD9ku4GpgL3SjoGOBf4fEQskPRj4OWIuKSy44iYBkwDGDBkRHTgHM3MrBV9YWS5s6Q5wErgmYho7kAfqrKuJZBOAK6KiGHAEcC1ktp6XQ8BTkz1zQK2A0ZExHpgHHAtcF9E3N+Bms3MrE76wshyUUSMSo9S75V0VETMaGcfSyne72wxjDcft55Mev8yImamEexgYEUb+hXw1Yi4q8q2EcDLwPbtrNXMzOqsL4wsAYiI5RSPQc/pYNvVkvZLs2BPBG5Lm58BDgKQtBswEHihla5WA1uWlu8CvtgyQ1fSLpI2lzQI+B7wEWA7Sce20t7MzLpBnwnL5FZgs8rJMWWSlgDfBcZJWpom1gB8EbgCWAgsAn6d1p8BnCJpLjAdGJcmAVXTDLwuaa6k01N/C4CHJc0HfkIx2r8U+GFEPEExcp0s6e3Ar4Bjqk3wMTOzrqPWf6/bhmrAkBEx5KSpjS7DrE9bMvnIRpdg7SRpdkSMrratr40szczM2q0vTPD5J5IOBaZUrF4cEcfUqf+RFDNZy9ZGxL716N/MzLpXnwzLNPu02gzUevU/DxjVVf2bmVn38mNYMzOzDIelmZlZhsPSzMwsw2FpZmaW0Scn+PR2I4cOosmf8TIzqxuPLM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmluGPjvRC85atYvjE2xtdhplZh/W0/+LMI0szM7MMh6WZmVmGw9LMzCzDYWlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbDsgMkvZy+D5f0iqQ5khZIukbSxjXabSfpHkkvS7q8YtvekuZJWijpMklK63dMbR6R1CzpiK49OzMzq+Sw7LxFETEKGAkMA46rse+rwH8AE6ps+xEwHhiRvg5L688DboiIPYGxwA/rU7aZmbWVw7JOImId8CdgaI19/h4Rf6QIzTdIGgJsFREzIyKAa4CjW5oBW6XXg4DnqvUtabykJklN69as6tS5mJnZWzks60TSQGBf4M4ONB8KLC0tL+XN0J0EfErSUuAO4KvVOoiIaRExOiJG999sUAdKMDOz1jgsO29nSXOAlcAzEdHcgT5UZV2k7ycAV0XEMOAI4FpJ/nMzM+tG/qXbeS3vWb4b2E/SUR3oYynF+50thvHm49aTgRsAImImMBAY3OFqzcys3RyWdRIRy4GJwDkdbLta0n5pFuyJwG1p8zPAQQCSdqMIyxfqUrSZmbWJw7K+bgU2k/Th1naQtAT4LjBO0lJJu6dNXwSuABYCi4Bfp/VnAKdImgtMB8alSUBmZtZNNmp0ARuiiNgifV8CvK+0PoA9Mm2Ht7K+qdxXaf0C4IMdr9bMzDrLI0szM7MMjyy7gKRDgSkVqxdHxDGNqMfMzDrHYdkFIuIu4K5G12FmZvXhx7BmZmYZDkszM7MMh6WZmVmGw9LMzCzDE3x6oZFDB9E0+chGl2Fm1mt4ZGlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswx/dKQXmrdsFcMn3t7oMszMutWSLvzInEeWZmZmGQ5LMzOzDIelmZlZhsPSzMwsw2FpZmaW4bA0MzPLcFiamZllOCzNzMwyHJZmZmYZDkszM7OMhoelpEslnVZavkvSFaXl70gKSbtWtJsq6awa/Z4jaaGkxyUdWlq/iaRpkp6Q9GdJn0jrvy5pgaRmSb+TtFMr/U6SNCG9vkrSYklzJM2VdFDmXL+SagpJg0vrJemytK1Z0l6lbadLelTSfEnTJQ2sdQwzM6u/hocl8ABwAICkfsBg4L2l7QcAdwJjW1ak/Y4Frq/WoaTd0/7vBQ4Dfiipf9p8LrAiInYBdgfuS+sfAUZHxPuBm4CL2lj/mRExCjgN+HFm3/uBg4GnK9YfDoxIX+OBH6XzGAqcmup6H9Cf0nUwM7Pu0RPC8n5SWFKE23xgtaRtJA0AdgPO5q0h8RFgSURUhk6LMcAvImJtRCwGFgL7pG2fBS4EiIj1EfFien1PRKxJ+zwIDGvnecwEhtbaISIeiYglrdR7TRQeBLaWNCRt2wjYVNJGwGbAc9X6ljReUpOkpnVrVrWzdDMzq6XhYRkRzwGvS9qRIjRnArOA/YHRQHNENAPrJe2Rmo0FptfodijwbGl5KTBU0tZp+ZuSHpZ0o6R3VGl/MvDrdp7KYcCt7WzTomq9EbEMuAR4BlgOrIqIu6t1EBHTImJ0RIzuv9mgDpZhZmbVNDwsk5bRZUtYziwtP5D2mQ6MTSOsMcCNNfpTlXVBMUobBtwfEXul41zylobSpyhC+uI21n6xpKeAnwPfbmObNtUraRuKc30nsD2wearPzMy6UU8Jy5b3LUdSPIZ9kGJkeQBFkEIRlsdRvOfXHBEravS3FNihtDyM4vHlSmANcEtafyNQnkxzMMV7mkdFxNq07ltpAs+cVo51JvBu4Dzg6jaca3vqPRhYHBEvRMRrwM28+cjazMy6SU8Jy/uBjwEvRcS6iHgJ2JoiMGcCRMQiirCbTO1HsAAzKEahAyS9k2LizJ8iIoBfAQem/Q4CFgBI2hP4CUVQvhHEEXFuRIxKk3iqioj1wPeAfuWZt+0wAzgxzYrdj+Jx63KKx6/7SdpMklK9j3WgfzMz64SeEpbzKGbBPlixblXLBJxkOvAe3hwZVhURjwI3UAThncCXI2Jd2nw2MElSM/Bp4Iy0/mJgC+DGNJKc0Z4TSEF8AVDr4yynSlpKMXJsLn1E5g7gKYqJSD8FvpT6nEUxM/dhiuvRD5jWnrrMzKzzVPyOt95kwJARMeSkqY0uw8ysWy2ZfGSn2kuaHRGjq23rKSNLMzOzHmujRhfQGen9wSkVqxdHxDGNqKeFpFsoZrCWnR0RdzWiHjMz65wNOixT+PS4AGp0WJuZWX35MayZmVmGw9LMzCzDYWlmZpbhsDQzM8vYoCf4WHUjhw6iqZOfNzIzszd5ZGlmZpbhsDQzM8twWJqZmWU4LM3MzDIclmZmZhkOSzMzswyHpZmZWYbD0szMLMNhaWZmlqGIaHQNVmeSVgOPN7qOKgYDLza6iCpcV/u4rvZxXe3TyLp2ioi3Vdvgf+6ud3o8IkY3uohKkppcV9u5rvZxXe3jutrHj2HNzMwyHJZmZmYZDsveaVqjC2iF62of19U+rqt9XFc7eIKPmZlZhkeWZmZmGQ5LMzOzDIdlDyfpMEmPS1ooaWKV7ZJ0WdreLGmvXFtJ20r6jaQn0/dtuqsuSTtIukfSY5IelfS1UptJkpZJmpO+juiuutK2JZLmpWM3ldZ3+np1pjZJu5auyRxJf5N0WtrWHdfsPZJmSloraUJb2nbTPVa1rh5wj9W6Xl12j3XiejX6/vpkut+bJT0gaY9c23r9TLZLRPirh34B/YFFwLuATYC5wO4V+xwB/BoQsB8wK9cWuAiYmF5PBKZ0Y11DgL3S6y2BJ0p1TQImNOJ6pW1LgMFV+u3U9apHbRX9/D+KD0931zV7O/AB4FvlY/WAe6y1uhp9j1Wtqyvvsc7W1eD76wBgm/T6cLrhd1hHvjyy7Nn2ARZGxFMR8Q/gF8CYin3GANdE4UFga0lDMm3HAFen11cDR3dXXRGxPCIeBoiI1cBjwNB2Hr/udWX67ez1qmdtBwGLIuLpDtTQoboiYkVEPAS81o62XX6PtVZXo++xGterloZdrwqNuL8eiIi/pMUHgWFtaFuPn8l2cVj2bEOBZ0vLS/nnH/rW9qnV9h0RsRyKXywUf+PsrrreIGk4sCcwq7T6K+lxzJUdeLTS2boCuFvSbEnjS/t09nrVo7YWY4HpFeu6+pp1pG133GNZDbrHaumqe6wu14vG318nUzxdybWtx89kuzgsezZVWVf5WZ/W9mlL247qTF3FRmkL4JfAaRHxt7T6R8DOwChgOfCdbq7rgxGxF8WjoC9L+kg7j9+VtSFpE+Ao4MbS9u64Zl3Rtsv7buA9VktX3WP1uF4Nvb8k/Q+KsDy7vW27g8OyZ1sK7FBaHgY818Z9arV9vuXxXvq+ohvrQtLGFL/ErouIm1t2iIjnI2JdRKwHfkrxGKbb6oqIlu8rgFtKx+/s9ep0bcnhwMMR8XzLim66Zh1p2x33WKsafI+1qgvvsU7VlTTs/pL0fuAKYExErGxD23r8TLaLw7JnewgYIemd6W99Y4EZFfvMAE5UYT9gVXosUavtDOCk9Pok4LbuqkuSgJ8Bj0XEd8sNKt6fOwaY3411bS5py1TH5sAhpeN39np1qrbS9hOoeETWTdesI2274x6rqgfcY63V1ZX3WGf+HFs05P6StCNwM/DpiHiijW3r8TPZPl05e8hfnf+imCH5BMWssHPTui8AX0ivBfwgbZ8HjK7VNq3fDvgd8GT6vm131QV8iOJRSjMwJ30dkbZdm/ZtpvhhGNKNdb2LYrbdXODRel+vOvxZbgasBAZV9Nkd1+xfKP6W/zfgr+n1Vj3gHqtaVw+4x1qrq0vvsU7+OTby/roC+Evpz6qpVtt6/ky258v/3J2ZmVmGH8OamZllOCzNzMwyHJZmZmYZDkszM7MMh6WZmVmGw9LMzCzDYWlmZpbx/wEwKKI2HtFGWgAAAABJRU5ErkJggg==\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
}