11_balanced_RF_model_balanced_accuracy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Balanced random forest model with balanced accuracy score"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use the balanced random forest model maximizing the balanced accuracy, i.e. the average of the recall of each class, in the grid search. The balanced Random Forest model is implemented in the [imbalanced-learn](https://imbalanced-learn.readthedocs.io/en/stable/index.html) software package."
   ]
  },
  {
   "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 imblearn.ensemble import BalancedRandomForestClassifier\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": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "datapath = '/data/ml_course/05_Capstone_project/rf_data/'\n",
    "model_name = 'balanced_RF_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": "code",
   "execution_count": 3,
   "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": [
    "## 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:  9.1min\n",
      "[Parallel(n_jobs=-1)]: Done  70 out of  70 | elapsed: 27.6min finished\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "GridSearchCV(cv=ShuffleSplit(n_splits=5, random_state=0, test_size=0.1, train_size=None),\n",
       "             estimator=Pipeline(steps=[('scaler', None),\n",
       "                                       ('RF',\n",
       "                                        BalancedRandomForestClassifier(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",
    "# Balanced random forest\n",
    "pipe_RF = Pipeline([\n",
    "    ('scaler', None),\n",
    "    ('RF', BalancedRandomForestClassifier(random_state=0, max_depth=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": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate model"
   ]
  },
  {
   "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": {
      "text/plain": [
       "dict_keys(['mean_fit_time', 'std_fit_time', 'mean_score_time', 'std_score_time', 'param_RF__n_estimators', 'params', 'split0_test_score', 'split1_test_score', 'split2_test_score', 'split3_test_score', 'split4_test_score', 'mean_test_score', 'std_test_score', 'rank_test_score', 'split0_train_score', 'split1_train_score', 'split2_train_score', 'split3_train_score', 'split4_train_score', 'mean_train_score', 'std_train_score'])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "grid_RF_cv.cv_results_.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEGCAYAAABrQF4qAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjeUlEQVR4nO3de3xV5Z3v8c8vd3IhIAlIiAoKclOBitSK7Ui9YW9o25nKHHtmmJnDaLGjczmjM3OmOuP01fHYztiKcxhPq3ZGq6etN1Ra7elo27HVghLkJjUCQthREjAJkOwkO/s3f+yVuAnBbCBhZ6/1fb9eeWWv296/B19+eXj2s55l7o6IiIRXXrYLEBGR4aWgFxEJOQW9iEjIKehFREJOQS8iEnIF2S5gIFVVVT558uRslyEikjNeffXVZnevHujYiAz6yZMns27dumyXISKSM8zs7aMd09CNiEjIKehFREJOQS8iEnIjcox+IN3d3TQ0NBCPx7NdSk4pKSmhtraWwsLCbJciIlmSM0Hf0NBARUUFkydPxsyyXU5OcHf27dtHQ0MDU6ZMyXY5IpIlOTN0E4/HGTdunEL+GJgZ48aN07+CRCIuZ4IeUMgfB/2ZiUjODN2ISDgkepIc7ExwIJ6gLd7NgXgi+Onu+92VSIIZBuSZYQYGqd9920ae0fe671hw3uHX2WHn5QWv6T2v771T2wTX5KW/92Hvk/b+A3x++mcf/hnv10xfLe9/fmF+HjMnjh7yP3MFfQ65++67Wb58OaWlpdkuRSKquyd5WCgPFNSpEO+mbYAAPxBP0N7Vk+1mjFhV5cWs+1+XDfn7KuhzyN133811112noJfj0pVIpoVuehgP0Kvu7A3yw4M63p0c9HNKCvOoKCmkoqSAipJCRpcUMLGyhIri9/elfr9/PH1feUkBxQX5QGpCgTsk3XHAHZzUvt7XSQ/O6z0eXOME1wXn4aTOJe09g+cuJdOucU+9J33nHf6ZyeCiw+tK+9372cnMaiatlsKC4RlNV9BnaOfOnSxevJiLL76Yl19+mTlz5rBs2TJuu+029u7dy8MPP8zs2bP58pe/zMaNG0kkEtx+++0sWbKEnTt38sUvfpFDhw4BsHLlSi666CJefPFFbr/9dqqqqti0aRPnn38+Dz300IDj6t/61reIxWIsWrSIqqoqXnjhBZ5//nluu+02Ojs7Oeuss3jggQcoLy8/2X80chJ0JnoGHOI4Wq/58F536nVnYvCQHlWYf0QI144Zddi+/mE9Om1feXEBRUMYVr1DJXnou6YTYSPxUYLz58/3/mvdbN26lZkzZwLwd09vZkusbUg/c1bNaG779OyjHt+5cydTp05l/fr1zJ49mwsuuIA5c+bwne98h9WrV/PAAw8wa9YsZs2axXXXXUdLSwsLFixg/fr1qXG5vDxKSkp48803Wbp0KevWrePFF19kyZIlbN68mZqaGhYuXMhdd93FxRdfPGANvWsAVVVV0dzczGc/+1l+9KMfUVZWxp133klnZydf+cpXjrgu/c9Osq8z0cP+Q13sO9jF/kOH/+w71MX+Q519r9s6UmHelUFIlxXlH9FbLi8peL/HXDxwWPcGdXlJAYX5OTU/Q9KY2avuPn+gY+rRH4MpU6Zw7rnnAjB79mwuvfRSzIxzzz2XnTt30tDQwOrVq/n6178OpKaE7tq1i5qaGm688Ubq6urIz8/nN7/5Td97LliwgNraWgDmzp3Lzp07jxr06V5++WW2bNnCwoULAejq6uIjH/nIUDdZMtDelTgstHvDet+hLvYfTN+X+jnYmRjwffIMTikr6vuZeepoKksL+/WaC44YAhkdBHp+nnq9MrCcDPoP6nkPp+Li4r7XeXl5fdt5eXkkEgny8/N57LHHmD59+mHX3X777UyYMIENGzaQTCYpKSkZ8D3z8/NJJAYOgf7cncsvv5xHHnnkRJok/bg7bfFEEMqdfQG+b6Be98Eu9rd3HXXcuig/ry+0x5UXcca40tTrsiJOKSvmlLLC4HdqX+WoQvIU1jIMcjLoR6orr7ySe+65h3vuuQczY/369cybN4/W1lZqa2vJy8vju9/9Lj09xzfroKKiggMHDlBVVcWFF17IihUrqK+vZ+rUqbS3t9PQ0MDZZ589xK3KbT1Jp6X98LDu7Wm/157W+w4C/b32Lrp7Bh7OLC3K7wvuqvJizp5Q0Rfa43p74uVFfa/Liwt0H4OMCAr6IfS3f/u33HzzzZx33nm4O5MnT+aZZ57hS1/6Ep/73Of4wQ9+wKJFiygrKzuu91++fDlXXXUVEydO5IUXXuDBBx9k6dKldHZ2AvAP//APoQ/6rkSS99r79a4Pdh4Z5MFPS3tXMIPiSBUlBX2hXDu2lDm1Yw4L6lRPu5hTyos4pbSIUUX5J7exIkMkJ7+MlWOTy39225sO8r9/vI033mlj36EuDsQHHtoyg7Gl6QF9+O9TylO97rGlqWGUsaVFQzo7RCTb9GWs5Jx4dw//8kI9q362neLCPH7r7Or3x7bLjwzyMaVF+jJS5CgU9CPQNddcw44dOw7bd+edd3LllVdmqaKT68Vte7lt9Wbe3tfO1XNr+OtPzmR8RcngF4rIgDIKejNbDHwTyAe+7e7/2O/4WOB+4CwgDvyBu2/K5Fo50hNPPJHtEk7I3rY4m2KtvNvWSXFBHsUF+ZQU5lFSmPqd2k5/nTq272AXdzyzhWc3NnJmVRkP/9GHWTi1KtvNEcl5gwa9meUD9wKXAw3AWjNb7e5b0k77a6DO3a8xsxnB+ZdmeG3G3F2zGI7RcH4Hk0w6u/a3sznWxuZYa/C7jeaDncf9nkUFefzZ5Wfzx791Zt9t8CJyYjLp0S8A6t19O4CZPQosAdLDehbwNQB3f8PMJpvZBODMDK7NSElJCfv27dOa9Meg98Ej6fP2j1d3T5L6vQcPC/WtsTYOBDf/FOQZU8eXc8n0ambXjGZ2TSW1Y0fRlUgST/QQ704S7+6hM5H6He/uobO791jqdY8718ybxBnjjm9WkogMLJOgnwTsTttuAD7c75wNwGeB/zSzBcAZQG2G1wJgZsuB5QCnn376Ecdra2tpaGigqakpg5KlV++jBI9VMulsaWzjpfpmXnprH2t37KejOzX/v6QwtZTq1fMm9YX6tAnllBSqBy4yEmUS9AN1n/uPB/wj8E0zqwM2AuuBRIbXpna63wfcB6nplf2PFxYW6nF4w8jd2bmvnZfqm/nlW8386q19vNfeDcDU8eX8zvxaPnTGWGbXjGZKVblmuIjkkEyCvgE4LW27Foiln+DubcAyAEuNq+wIfkoHu1ayZ29bnF++tY//rG/ml/XNxFpTjxycWFnCx2dM4OJp47jorComjNaMF5FclknQrwWmmdkUYA9wLfC76SeY2Rig3d27gD8Cfu7ubWY26LVycr3bFuf/rd3N0xtivLn3IABjSgv5yJnjuGFRFQvPGseUqjJ9DyISIoMGvbsnzOxG4DlSUyTvd/fNZnZ9cHwVMBP4NzPrIfVF6x9+0LXD0xQ5mmTS+c/6Zh5+5W3+/9a99CSdj5w5js+fX8vCqVXMmjhai2mJhFjOLIEgx675YCc/WNfAI7/exa797ZxSVsRvn1/L0gWnM7lKM1tEwkRLIESIu/PKjv08/Moufrypke4eZ8GUU/jzK85m8Tmnam66SAQp6EOipb2Lx17bw/deeZu3mg4xuqSA6y48g99dcDrTJlRkuzwRySIFfY7bvb+de1+o54n1e+hMJJl3+hju+vx5fOq8Gi2rKyKAgj5nvdsWZ+V/1PPo2l2YGZ8/v5b/9uHTmV1Tme3SRGSEUdDnmH0HO1n1s7f4t1+9TU/S+cIFp3Hjx6cysXJUtksTkRFKQZ8jWju6+b8/3879L+0g3t3DNfNquenSaZw+rjTbpYnICKegH+EOdSZ44KUd3Pfz7bTFE3zqvIncfNnZTB1fnu3SRCRHKOhHqHh3Dw+9/Db/8uJb7D/UxWUzJ/Bnl5/NrJrR2S5NRHKMgn6ESfQkeWTtbu756ZvsPdDJR6dV8edXTGfuaWOyXZqI5CgF/Qiys/kQf/r9OtbvauGCyWO5Z+k8PnzmuGyXJSI5TkE/Arg7D7+yi68+u5Wigjy+tXQenz5vohYWE5EhoaDPsr1tcf7ysdd5cVsTH51WxV2fn8OplVoWWESGjoI+i559vZG/eXIj8e4e7lgym+suPEO9eBEZcgr6LGjt6Oa2pzbxZF2MObWV/NMX5nJWtaZLisjwUNCfZC/VN/MXP9jA3gOd3HzZNFYsmkphfl62yxKREFPQnyTx7h7u/PEbPPDSTs6sLuPxGy5ijqZMishJoKA/CTY2tPKn36+jfu9Bfv+iydyyeIZWlhSRk0ZBP8z+/Vc7+bunt1BVXsy//+ECPjqtOtsliUjEKOiH0c9+08RXVm9m0fTx/PPvzKWytDDbJYlIBGX0LaCZLTazbWZWb2a3DnC80syeNrMNZrbZzJalHbvJzDYF+28ewtpHtN372/mTR9YzfUIFK393nkJeRLJm0KA3s3zgXuAqYBaw1Mxm9TttBbDF3ecAlwDfMLMiMzsH+B/AAmAO8CkzmzaE9Y9IHV09LP/3V3F3/vWL51NapH84iUj2ZNKjXwDUu/t2d+8CHgWW9DvHgQpL3e1TDuwHEsBM4GV3b3f3BPAz4Johq34Ecnf+6vHXeeOdNr65dB5njCvLdkkiEnGZBP0kYHfadkOwL91KUqEeAzYCN7l7EtgEfMzMxplZKfAJ4LSBPsTMlpvZOjNb19TUdIzNGDkeeGknT9bF+LPLzmbR9PHZLkdEJKOgH+iefO+3fSVQB9QAc4GVZjba3bcCdwI/AX4MbCDV0z/yDd3vc/f57j6/ujo3Z6a8vH0fX12zlctnTWDFoqnZLkdEBMgs6Bs4vBdeS6rnnm4Z8Lin1AM7gBkA7v4dd/+Qu3+M1JDOmyde9sjT2NrBjd97jTNOKeUbvzOHvDytWSMiI0MmQb8WmGZmU8ysCLgWWN3vnF3ApQBmNgGYDmwPtscHv08HPgs8MjSljxydiR5ueOg1Orp6uO+/n8/oEs2wEZGRY9DpIO6eMLMbgeeAfOB+d99sZtcHx1cBdwAPmtlGUkM9t7h7c/AWj5nZOKAbWOHu7w1HQ7Lp9tWbqdvdwqrrPsTU8RXZLkdE5DAZzftz9zXAmn77VqW9jgFXHOXaj55IgSPdI7/exSO/3s0Nl5zF4nMmZrscEZEjaNnEE7B+13vc9tRmPjqtir+4Ynq2yxERGZCC/jg1HejkhodeY/zoYr517Tzy9eWriIxQumXzOHT3JFnxvdd4r72Lx790EWPLirJdkojIUSnoj8PX1rzBr3fs55+/MIfZNZXZLkdE5ANp6OYYPVW3h/tf2sHvXzSZa+bVZrscEZFBKeiPQdOBTm557HUWTD6Fv/nkzGyXIyKSEQX9MXh6Q4x4d5KvXnOOnvMqIjlDaXUMnqrbw+ya0UyboJuiRCR3KOgztKP5EBsaWlkytybbpYiIHBMFfYaeqtuDGXxmTv8VmkVERjYFfQbcnafqYlw4ZRynVpZkuxwRkWOioM/A6w2t7Gg+xNXzNGwjIrlHQZ+Bp+piFOXnadEyEclJCvpB9CSdp1+PsWhGNZWjtM68iOQeBf0gfvlWM00HOrl6rr6EFZHcpKAfxJPrY1QUF7Bohh70LSK5SUH/AeLdPTy3+R0Wn3MqJYX52S5HROS4KOg/wE+37uVgZ4Kr52nYRkRyl4L+AzxZt4fxFcVceOa4bJciInLcMgp6M1tsZtvMrN7Mbh3geKWZPW1mG8xss5ktSzv2p8G+TWb2iJnlxB1HLe1dvLhtL5+eU6OnR4lIThs06M0sH7gXuAqYBSw1s1n9TlsBbHH3OcAlwDfMrMjMJgF/Asx393OAfODaIax/2Pxo0zt097hm24hIzsukR78AqHf37e7eBTwKLOl3jgMVZmZAObAfSATHCoBRZlYAlAKxIal8mD25fg9nVpdxzqTR2S5FROSEZBL0k4DdadsNwb50K4GZpEJ8I3CTuyfdfQ/wdWAX0Ai0uvvzA32ImS03s3Vmtq6pqekYmzG0Yi0dvLJjP0vmTCL1d5eISO7KJOgHSjrvt30lUAfUAHOBlWY22szGkur9TwmOlZnZdQN9iLvf5+7z3X1+dXV1huUPj9UbUv/o0JLEIhIGmQR9A3Ba2nYtRw6/LAMe95R6YAcwA7gM2OHuTe7eDTwOXHTiZQ+vp+pizD1tDJOryrJdiojICcsk6NcC08xsipkVkfoydXW/c3YBlwKY2QRgOrA92H+hmZUG4/eXAluHqvjh8Jt3D7C1sY2r1ZsXkZAoGOwEd0+Y2Y3Ac6Rmzdzv7pvN7Prg+CrgDuBBM9tIaqjnFndvBprN7IfAa6S+nF0P3Dc8TRkaT67fQ36e8cnzFPQiEg6DBj2Au68B1vTbtyrtdQy44ijX3gbcdgI1njS9DxhZOLWK6oribJcjIjIkdGdsmlfffo89LR0athGRUFHQp3mybg8lhXlcMfvUbJciIjJkFPSB7p4kz77eyGUzJ1BenNGIlohITlDQB37xZhPvtXdryQMRCR0FfeDJ9THGlBbysbOze7OWiMhQU9ADhzoT/GTLu3zi3IkUFeiPRETCRakGPL/lHTq6ezRsIyKhpKAnteTBpDGjmH/G2GyXIiIy5CIf9M0HO/nFm818Zm4NeXrAiIiEUOSD/tnXG+lJulaqFJHQinzQP1m3hxmnVjDjVD1gRETCKdJBv2tfO+t3tbBEX8KKSIhFOuif3dgIwKfnTMxyJSIiwyfSQf9W00FOHV1C7djSbJciIjJsIh30sZYOJo4pyXYZIiLDKtJB39gap6ZyVLbLEBEZVpENencn1tJBjXr0IhJykQ3699q76UwkmagevYiEXGSDPtbSAaAevYiEXkZBb2aLzWybmdWb2a0DHK80s6fNbIOZbTazZcH+6WZWl/bTZmY3D3EbjktjaxxAPXoRCb1BH6VkZvnAvcDlQAOw1sxWu/uWtNNWAFvc/dNmVg1sM7OH3X0bMDftffYATwxxG45LY2uqR69ZNyISdpn06BcA9e6+3d27gEeBJf3OcaDCzAwoB/YDiX7nXAq85e5vn2DNQyLWEqcw36gqK852KSIiwyqToJ8E7E7bbgj2pVsJzARiwEbgJndP9jvnWuCRo32ImS03s3Vmtq6pqSmDsk5MrKWDUytLtGKliIReJkE/UBJ6v+0rgTqghtRQzUoz61slzMyKgM8APzjah7j7fe4+393nV1cP/+P8Gls7NIdeRCIhk6BvAE5L264l1XNPtwx43FPqgR3AjLTjVwGvufu7J1LsUIq1xKkZo6AXkfDLJOjXAtPMbErQM78WWN3vnF2kxuAxswnAdGB72vGlfMCwzcnWk3TebYszsVJfxIpI+A0668bdE2Z2I/AckA/c7+6bzez64Pgq4A7gQTPbSGqo5xZ3bwYws1JSM3b+eJjacMyaD3aSSDoT1aMXkQgYNOgB3H0NsKbfvlVpr2PAFUe5th0YdwI1Drm+m6XUoxeRCIjknbG6WUpEoiSSQa/lD0QkSiIa9HFKi/KpHFWY7VJERIZdJIO+sbWDiZUlpG7kFREJt0gGfaxVc+hFJDoiGfSNLR2aQy8ikRG5oO9KJGk62KkZNyISGZEL+nfb4rhrxo2IREfkgr53aqV69CISFZEL+t6bpfRlrIhEReSCPtaqm6VEJFoiF/SNLXEqRxVSWpTRMj8iIjkvekHfqqmVIhItkQt6PXBERKImckGvHr2IRE2kgr6jq4f32rvVoxeRSIlU0GvGjYhEUaSCvrFFDxwRkeiJVND39egV9CISIRkFvZktNrNtZlZvZrcOcLzSzJ42sw1mttnMlqUdG2NmPzSzN8xsq5l9ZCgbcCx6e/QTKouzVYKIyEk3aNCbWT5wL3AVMAtYamaz+p22Atji7nOAS4BvmFlRcOybwI/dfQYwB9g6RLUfs8bWDqrKiykuyM9WCSIiJ10mPfoFQL27b3f3LuBRYEm/cxyosNQjm8qB/UDCzEYDHwO+A+DuXe7eMlTFH6vUA0f0RayIREsmQT8J2J223RDsS7cSmAnEgI3ATe6eBM4EmoAHzGy9mX3bzMoG+hAzW25m68xsXVNT07G2IyMxPXBERCIok6Af6MGq3m/7SqAOqAHmAiuD3nwB8CHg/7j7POAQcMQYP4C73+fu8919fnV1dWbVHwN3p7GlQ3PoRSRyMgn6BuC0tO1aUj33dMuAxz2lHtgBzAiubXD3V4Lzfkgq+E+6tniCQ109mnEjIpGTSdCvBaaZ2ZTgC9ZrgdX9ztkFXApgZhOA6cB2d38H2G1m04PzLgW2DEnlx6gxmFo5UWP0IhIxg67V6+4JM7sReA7IB+53981mdn1wfBVwB/CgmW0kNdRzi7s3B2/xZeDh4C+J7aR6/yedbpYSkajKaFF2d18DrOm3b1Xa6xhwxVGurQPmH3+JQ0PLH4hIVEXmztjGljj5ecb4CgW9iERLZII+1tLBhIpi8vMGmkQkIhJe0Qn6Vk2tFJFoikzQN7bGmaigF5EIikTQuzuNrXFqdFesiERQJIJ+36EuuhJJLX8gIpEUiaDvm0OvoRsRiaBIBL0eOCIiURaNoG/RzVIiEl2RCPrG1jjFBXmcUlY0+MkiIiETiaDvXYc+9VwUEZFoiUTQN7bGtZiZiERWNIK+pUPLE4tIZIU+6HuSzrsHOjXjRkQiK/RBv/dAnJ6kq0cvIpEV+qB/f2qlevQiEk0RCPrUXbEauhGRqAp90OtZsSISdaEP+lhLnPLiAkaXFGa7FBGRrMgo6M1ssZltM7N6M7t1gOOVZva0mW0ws81mtizt2E4z22hmdWa2biiLz0Rja4dWrRSRSBv04eBmlg/cC1wONABrzWy1u29JO20FsMXdP21m1cA2M3vY3buC44vcvXmoi8+EHjgiIlGXSY9+AVDv7tuD4H4UWNLvHAcqLLXGQDmwH0gMaaXHKdaiB46ISLRlEvSTgN1p2w3BvnQrgZlADNgI3OTuyeCYA8+b2atmtvxoH2Jmy81snZmta2pqyrgBH6Qz0UPzwU5NrRSRSMsk6AdaCcz7bV8J1AE1wFxgpZmNDo4tdPcPAVcBK8zsYwN9iLvf5+7z3X1+dXV1JrUP6p3W4IEj6tGLSIRlEvQNwGlp27Wkeu7plgGPe0o9sAOYAeDuseD3XuAJUkNBJ0XfHHr16EUkwjIJ+rXANDObYmZFwLXA6n7n7AIuBTCzCcB0YLuZlZlZRbC/DLgC2DRUxQ+mbw69evQiEmGDzrpx94SZ3Qg8B+QD97v7ZjO7Pji+CrgDeNDMNpIa6rnF3ZvN7EzgiWAd+ALge+7+42FqyxEa+4Zu1KMXkegaNOgB3H0NsKbfvlVpr2Okeuv9r9sOzDnBGo9brKWDsaWFjCrKz1YJIiJZF+o7Y/XAERGRkAd9rKVDX8SKSORFIOj1RayIRFtog/5QZ4K2eEJDNyISeaEN+t6plerRi0jUhTboe2+WUo9eRKIutEGvm6VERFJCG/SxljhmcKqCXkQiLrRB39jawfiKYgrzQ9tEEZGMhDYFYy26WUpEBMIc9K2aQy8iAiENenenUT16EREgpEHf2tFNR3ePZtyIiBDSoNcDR0RE3hfKoNccehGR94Uy6GPBA0cmqUcvIhLSoG/poDDfqCovznYpIiJZF8qgb2zpYMLoEvLyLNuliIhkXSiDPtYap0ZTK0VEgAyD3swWm9k2M6s3s1sHOF5pZk+b2QYz22xmy/odzzez9Wb2zFAV/kEaWzuYqJulRESADILezPKBe4GrgFnAUjOb1e+0FcAWd58DXAJ8w8yK0o7fBGwdkooHkUw67+hZsSIifTLp0S8A6t19u7t3AY8CS/qd40CFmRlQDuwHEgBmVgt8Evj2kFX9AZoPddLd41r+QEQkkEnQTwJ2p203BPvSrQRmAjFgI3CTuyeDY3cDfwkk+QBmttzM1pnZuqampgzKGlhj781S6tGLiACZBf1AU1e83/aVQB1QA8wFVprZaDP7FLDX3V8d7EPc/T53n+/u86urqzMoa2CxluBmKfXoRUSAzIK+ATgtbbuWVM893TLgcU+pB3YAM4CFwGfMbCepIZ+Pm9lDJ1z1B+i9WUo9ehGRlEyCfi0wzcymBF+wXgus7nfOLuBSADObAEwHtrv7X7l7rbtPDq77D3e/bsiqH0BjSwclhXmMKS0czo8REckZBYOd4O4JM7sReA7IB+53981mdn1wfBVwB/CgmW0kNdRzi7s3D2PdR9UYzKFPfS8sIiKDBj2Au68B1vTbtyrtdQy4YpD3eBF48ZgrPEYxzaEXETlM6O6M1QNHREQOF6qg7+5JsvdAXOvQi4ikCVXQv9sWJ+lQo3XoRUT6hCroG4OplRPVoxcR6ROqoO+9WUo9ehGR94Uq6NWjFxE5UriCvqWDipICyoszmjUqIhIJoQp6PXBERORIoQr6xtYOLU8sItJPqII+1hLX+LyISD+hCfpk0vmts6u5YPLYbJciIjKihOZby7w845+/MDfbZYiIjDih6dGLiMjAFPQiIiGnoBcRCTkFvYhIyCnoRURCTkEvIhJyCnoRkZBT0IuIhJy5e7ZrOIKZNQFvD3JaFdB8EsoZadTuaFG7o+VE2n2Gu1cPdGBEBn0mzGydu8/Pdh0nm9odLWp3tAxXuzV0IyIScgp6EZGQy+Wgvy/bBWSJ2h0tane0DEu7c3aMXkREMpPLPXoREcmAgl5EJORyLujNbLGZbTOzejO7Ndv1DBczu9/M9prZprR9p5jZT8zszeB36B6nZWanmdkLZrbVzDab2U3B/lC33cxKzOzXZrYhaPffBftD3e5eZpZvZuvN7JlgOyrt3mlmG82szszWBfuGvO05FfRmlg/cC1wFzAKWmtms7FY1bB4EFvfbdyvwU3efBvw02A6bBPDn7j4TuBBYEfw3DnvbO4GPu/scYC6w2MwuJPzt7nUTsDVtOyrtBljk7nPT5s8PedtzKuiBBUC9u2939y7gUWBJlmsaFu7+c2B/v91LgO8Gr78LXH0yazoZ3L3R3V8LXh8g9T//JELedk85GGwWBj9OyNsNYGa1wCeBb6ftDn27P8CQtz3Xgn4SsDttuyHYFxUT3L0RUoEIjM9yPcPKzCYD84BXiEDbg+GLOmAv8BN3j0S7gbuBvwSSafui0G5I/WX+vJm9ambLg31D3vZcezi4DbBP80NDyMzKgceAm929zWyg//Th4u49wFwzGwM8YWbnZLmkYWdmnwL2uvurZnZJlsvJhoXuHjOz8cBPzOyN4fiQXOvRNwCnpW3XArEs1ZIN75rZRIDg994s1zMszKyQVMg/7O6PB7sj0XYAd28BXiT1HU3Y270Q+IyZ7SQ1FPtxM3uI8LcbAHePBb/3Ak+QGp4e8rbnWtCvBaaZ2RQzKwKuBVZnuaaTaTXwe8Hr3wOeymItw8JSXffvAFvd/Z/SDoW67WZWHfTkMbNRwGXAG4S83e7+V+5e6+6TSf3//B/ufh0hbzeAmZWZWUXva+AKYBPD0PacuzPWzD5BakwvH7jf3b+a3YqGh5k9AlxCatnSd4HbgCeB7wOnA7uA33b3/l/Y5jQzuxj4BbCR98ds/5rUOH1o225m55H64i2fVAfs++7+92Y2jhC3O10wdPMX7v6pKLTbzM4k1YuH1DD699z9q8PR9pwLehEROTa5NnQjIiLHSEEvIhJyCnoRkZBT0IuIhJyCXkQk5BT0EmlmNsbMvpTtOkSGk4Jeom4McETQByulioSCgl6i7h+Bs4L1wNcGa+F/D9gYLDJ2V7D/dTP7496LzOx/pu3vXTu+zMyeDdaU32RmX8hWo0TS5dqiZiJD7VbgHHefG9yZ+WywvSNYTbDV3S8ws2LgJTN7HpgW/CwgtdDeajP7GFANxNz9kwBmVnnymyNyJAW9yOF+7e47gtdXAOeZ2eeD7UpSAX9F8LM+2F8e7P8F8HUzuxN4xt1/cfLKFjk6Bb3I4Q6lvTbgy+7+XPoJZnYl8DV3/9f+F5vZ+cAngK+Z2fPu/vfDWq1IBjRGL1F3AKg4yrHngBuCZZMxs7ODVQafA/4gWDMfM5tkZuPNrAZod/eHgK8DHxr+8kUGpx69RJq77zOzlyz1EPYOUiuF9vo2MBl4LVg+uQm42t2fN7OZwK+CB6IcBK4DpgJ3mVkS6AZuOHktETk6rV4pIhJyGroREQk5Bb2ISMgp6EVEQk5BLyIScgp6EZGQU9CLiIScgl5EJOT+C4yoQgpli4c1AAAAAElFTkSuQmCC\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": 17,
   "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>12</th>\n",
       "      <td>40</td>\n",
       "      <td>0.931180</td>\n",
       "      <td>0.904946</td>\n",
       "      <td>0.001739</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>50</td>\n",
       "      <td>0.930777</td>\n",
       "      <td>0.904835</td>\n",
       "      <td>0.001492</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>30</td>\n",
       "      <td>0.931807</td>\n",
       "      <td>0.903650</td>\n",
       "      <td>0.001831</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>20</td>\n",
       "      <td>0.932786</td>\n",
       "      <td>0.901838</td>\n",
       "      <td>0.002367</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>17</td>\n",
       "      <td>0.928165</td>\n",
       "      <td>0.900631</td>\n",
       "      <td>0.002120</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>19</td>\n",
       "      <td>0.928291</td>\n",
       "      <td>0.900477</td>\n",
       "      <td>0.002057</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>15</td>\n",
       "      <td>0.927913</td>\n",
       "      <td>0.899779</td>\n",
       "      <td>0.002007</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>13</td>\n",
       "      <td>0.927432</td>\n",
       "      <td>0.898618</td>\n",
       "      <td>0.002171</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>11</td>\n",
       "      <td>0.926978</td>\n",
       "      <td>0.896824</td>\n",
       "      <td>0.001588</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>9</td>\n",
       "      <td>0.926059</td>\n",
       "      <td>0.894534</td>\n",
       "      <td>0.002466</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   trees   mean_tr   mean_te    std_te\n",
       "12    40  0.931180  0.904946  0.001739\n",
       "13    50  0.930777  0.904835  0.001492\n",
       "11    30  0.931807  0.903650  0.001831\n",
       "10    20  0.932786  0.901838  0.002367\n",
       "8     17  0.928165  0.900631  0.002120\n",
       "9     19  0.928291  0.900477  0.002057\n",
       "7     15  0.927913  0.899779  0.002007\n",
       "6     13  0.927432  0.898618  0.002171\n",
       "5     11  0.926978  0.896824  0.001588\n",
       "4      9  0.926059  0.894534  0.002466"
      ]
     },
     "execution_count": 17,
     "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": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test accuracy: 90.77%\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": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = grid_RF_cv.predict(X_te)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "tn, fp, fn, tp = confusion_matrix(y_te, y_pred).ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True positive:  5782\n",
      "True negative:  635883\n",
      "False positive:  101278\n",
      "False negative:  286\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": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWgAAAEWCAYAAABLzQ1kAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAr3ElEQVR4nO3de5xXVb3/8deb4X6VuwheUDEz73K81Mk0S60sraMnypLKfqaZlpc8Wimmx47WMc1rmeL9Apop5ZVQTloq4iVRvICigNyEAeTmMJfP74+9Br4M35n5AjPMZub9fDz2Y/Zee62113cGPrNm7bXXVkRgZmb5066lG2BmZsU5QJuZ5ZQDtJlZTjlAm5nllAO0mVlOOUCbmeWUA7RtEkldJP1F0lJJ925CPcdLerwp29ZSJH1a0pst3Q7b8snzoNsGSd8EzgR2BZYBLwOXRMTTm1jvt4HTgE9GRNWmtjPvJAUwLCKmt3RbrPVzD7oNkHQmcCXwK2AgsB1wHXB0E1S/PfBWWwjOpZDUvqXbYK1IRHhrxRvQC1gOHNdAnk5kAXxO2q4EOqVzhwCzgbOABcBc4Lvp3C+B1UBlusaJwIXAHQV17wAE0D4dfwd4h6wXPwM4viD96YJynwSeB5amr58sODcRuBj4R6rncaBfPZ+ttv3nFLT/GOCLwFtAOfCzgvz7A88AS1Lea4CO6dzf02dZkT7v1wvq/y9gHnB7bVoqs1O6xr7peBtgIXBIS//b8Jb/zT3o1u8goDPw5wby/Bw4ENgb2IssSP2i4PzWZIF+MFkQvlZS74gYRdYrHxMR3SPipoYaIqkbcBXwhYjoQRaEXy6Srw/wUMrbF/gt8JCkvgXZvgl8FxgAdATObuDSW5N9DwYDFwB/BL4F7Ad8GrhA0o4pbzVwBtCP7Ht3GPBDgIg4OOXZK33eMQX19yH7a+KkwgtHxNtkwftOSV2Bm4FbImJiA+01AzzE0Rb0BRZGw0MQxwMXRcSCiPiArGf87YLzlel8ZUQ8TNZ7/NhGtqcG2F1Sl4iYGxGvFcnzJWBaRNweEVURcTfwBvDlgjw3R8RbEbEKGEv2y6U+lWTj7ZXAPWTB93cRsSxd/zVgT4CIeCEink3XfRf4A/CZEj7TqIioSO1ZR0T8EZgGPAcMIvuFaNYoB+jWbxHQr5Gx0W2A9wqO30tpa+qoE+BXAt03tCERsYJsWOBkYK6khyTtWkJ7ats0uOB43ga0Z1FEVKf92gA6v+D8qtryknaR9FdJ8yR9SPYXQr8G6gb4ICI+aiTPH4HdgasjoqKRvGaAA3Rb8AzwEdm4a33mkP15Xmu7lLYxVgBdC463LjwZEY9FxOfJepJvkAWuxtpT26b3N7JNG+J6snYNi4iewM8ANVKmwalQkrqTjevfBFyYhnDMGuUA3cpFxFKycddrJR0jqaukDpK+IOnXKdvdwC8k9ZfUL+W/YyMv+TJwsKTtJPUCzqs9IWmgpK+ksegKsqGS6iJ1PAzsIumbktpL+jqwG/DXjWzThugBfAgsT737U+qcnw/suF6phv0OeCEivk82tv77TW6ltQkO0G1ARPyWbA70L4APgFnAj4AHUpb/BiYDrwBTgBdT2sZcazwwJtX1AusG1XZks0HmkM1s+AzpBlydOhYBR6W8i8hmYBwVEQs3pk0b6GyyG5DLyHr3Y+qcvxC4VdISSf/ZWGWSjgaOJBvWgeznsK+k45usxdZq+UEVM7Occg/azCynHKDNzHLKAdrMLKccoM3McqpNL+zSr09Z7LBth5Zuhm2At6b3bTyT5cqylXMWRkT/jS1/xKHdYlF5sdmY63vhlYrHIuLIjb1W3rTpAL3Dth2Y9Ni2Ld0M2wBHHv3txjNZroyfNKruU6EbZFF5NZMe266kvGWDpjX21OcWpU0HaDPLvwBqqGnpZrQIB2gzy7UgqIzShjhaGwdoM8s996DNzHIoCKrb6BPPDtBmlns1DS8Y2Go5QJtZrgVQ7QBtZpZPbbUH7ScJzSzXAqiMKGkrhaStJN0n6Q1Jr0s6SFIfSeMlTUtfexfkP0/SdElvSjqiIH0/SVPSuaskKaV3kjQmpT8naYeCMiPTNaZJGtlYWx2gzSzXgqC6xK1EvwMejYhdyV6S/DpwLjAhIoYBE9IxknYDRgCfIFvX+zpJZame68leEjwsbbVPMJ4ILI6InYErgMtSXX2AUcABZC9mHlX4i6AYB2gzy7eA6hK3xkjqCRxM9voxImJ1RCwBjgZuTdluZe0r4o4G7kkvBJ4BTAf2lzQI6BkRz0S2qP5tdcrU1nUfcFjqXR8BjI+I8ohYDIxnbVAvygHazHIte5KwtK0EO5K9VehmSS9JujG9gm1gRMwFSF8HpPyDyd5AVGt2Shuc9uumr1MmvWx5KdC3gbrq5QBtZjknqkvcyN5gP7lgO6lOZe2BfYHrI2Ifspccn9vgxdcXDaRvbJmiPIvDzHItu0nY2IvV11gYEcMbOD8bmB0Rz6Xj+8gC9HxJgyJibhq+WFCQv3BFtSFk79ScnfbrpheWmS2pPdCL7B2cs4FD6pSZ2NCHcQ/azHItmwddcg+64boi5gGzJH0sJR0GTAXGAbWzKkYCD6b9ccCINDNjKNnNwElpGGSZpAPT+PIJdcrU1nUs8EQap34MOFxS73Rz8PCUVi/3oM0s92pK70GX4jTgTkkdgXeA75J1VsdKOhGYCRwHEBGvSRpLFsSrgFMj1qzcdApwC9AFeCRtkN2AvF3SdLKe84hUV7mki4HnU76LIqK8oYY6QJtZrtX2oJusvoiXgWLDIIfVk/8S4JIi6ZOB3Yukf0QK8EXOjQZGl9pWB2gzy7VAVLfR0VgHaDPLvSYe4thiOECbWa4FYnWUNZ6xFXKANrNcyx5U8RCHmVkuNeVNwi2JA7SZ5VqEqA73oM3McqnGPWgzs/zJbhK2zVDVNj+1mW0xfJPQzCzHqj0P2swsf/wkoZlZjtV4FoeZWf5kiyU5QJuZ5U4gKv2ot5lZ/kTgB1XMzPJJflDFzCyPAvegzcxyyzcJzcxyKJAX7Dczy6MAKr0Wh5lZHsnrQZuZ5VHgJwnNzHLLPWgzsxyKkHvQZmZ5lN0k9KPeZmY51HbfSdg2P7WZbTGym4QqaSuFpHclTZH0sqTJKa2PpPGSpqWvvQvynydpuqQ3JR1RkL5fqme6pKskKaV3kjQmpT8naYeCMiPTNaZJGtlYWx2gzSz3qmlX0rYBDo2IvSNieDo+F5gQEcOACekYSbsBI4BPAEcC10mqHW+5HjgJGJa2I1P6icDiiNgZuAK4LNXVBxgFHADsD4wq/EVQjAO0meVa7ZOETdWDrsfRwK1p/1bgmIL0eyKiIiJmANOB/SUNAnpGxDMREcBtdcrU1nUfcFjqXR8BjI+I8ohYDIxnbVAvygHazHKvhnYlbUA/SZMLtpOKVBfA45JeKDg/MCLmAqSvA1L6YGBWQdnZKW1w2q+bvk6ZiKgClgJ9G6irXr5JaGa5FgGVNSX3JRcWDFvU51MRMUfSAGC8pDcayFusWx4NpG9smaLcgzazXMuGONqVtJVUX8Sc9HUB8Gey8eD5adiC9HVByj4b2Lag+BBgTkofUiR9nTKS2gO9gPIG6qqXA7SZ5V51Wo+jsa0xkrpJ6lG7DxwOvAqMA2pnVYwEHkz744ARaWbGULKbgZPSMMgySQem8eUT6pSpretY4Ik0Tv0YcLik3unm4OEprV4e4siZ5UvLuOLsbXn3jc5IcOZvZ/L8Ez155rFeSLBVv0rOvnImfbeuYt6sjvy/z+zKkB0rANh1vxX8+LJsWOzJP2/FPVcPRII+Ayv5r6vfo1ffahbM7sBvfrIdK5aWUVMjvvezOex/2DLmz+7ARScOpaZaVFXB0d9byFEnLGrJb8UW5YzTnuGA4bNZsrQzJ5/+5U2u73OHvs03/vNVAO4euzt/e3InAM46/Z/ssft8VqzoCMDlVx3EOzP6bPL18qx2ml0TGQj8Oc2Iaw/cFRGPSnoeGCvpRGAmcBxARLwmaSwwFagCTo2I6lTXKcAtQBfgkbQB3ATcLmk6Wc95RKqrXNLFwPMp30URUd5QY5stQEsK4LcRcVY6PhvoHhEXllj+QmB5RPzvBlzzYeCbEbFE0vKI6L7hLW9Z118wmOGHfMj5f3yXytWiYlU7tv/YAkaeMw+AB27sxx1XbL0mEA/avoLr//bmOnVUV2X1/HHiG/TqW82NFw9i3M39+fbZ87jrdwM5+MtL+PLIRbz3VifO/9ZO3DZpKn0GVHHFuGl07BSsWtGOHxy6KwcdvpS+W1dt9u/Blmj8hB35y0O7cPZP/rlB5X79349z+VWfZP6Ctf9Uu3ev4PgRUzjtrC9AwNW/fYRnJw1h+YpOANx4y748/c/tm7T9+dZ0j3pHxDvAXkXSFwGH1VPmEuCSIumTgd2LpH9ECvBFzo0GRpfa3ubsQVcAX5P0PxGxsBmvs0ZEfHFzXKe5rFjWjinPduPsK2cC0KFj0KFj9Tp5PlrVDjXSmYgAQny0qh09o5oVy8vYZugqACRYuSybxrniwzL6DKxcc61alRWipqaJPlQb8erUgQwcsHydtEFbL+PUH0yiV88KKirKuPLaA5n9fq9G6xq+zxxeenlrli/PAvJLL2/N8H3nMPGpoc3S9i2B30nY9KqAG4AzgJ8XnpC0Pdlvkf7AB8B3I2JmkTp2kzQR2A64MiKuSuUfIBts7wz8LiJuSOnvAsM31y+EpjbvvU706lvF5WdsxzuvdWbYnqs45eL36dy1hpsv3Zq/3duHbj2r+fV909eWmdmRH35+F7r2qGHkf81ljwNW0L4DnHbpLE7+7K507lrDNkMr+NGvsh73t86ax8++sRPjbu7HRyvbcemYt9fUteD9Dlxwwo7MmdGJ758/x73nTXT6D5/j6uv3Z87cnnxsl4X86ORJnHv+5xst17fvKj5Y2G3N8cJFXenbd9Wa4+9862WO//oUXnpla26+dR8qq1r3OhXZLI7W/Rnr09xj0NcCr0j6dZ30a4DbIuJWSd8DrmLtJO9CuwKHAj2ANyVdHxGVwPfSeE4X4HlJf0p/ojQqzXs8CWC7wfkagq+uhulTunLqf7/Prvuu5PrzBzPmmgGMPGce3z032+65egDjRvfnhJ/Oo8+ASu54fio9+1Qz7ZUuXPjdodww8Q06da7hr7f149rH32TQ9qu59ueDGXP1QL75k/lMfKA3n//Pco49+QOmTu7Kr0/bnj88+Qbt2sGAwZX8fsKbLJrXngu/N5RPH7WE3v0dpDdG586V7LbrB/z8nKfWpHXokP019PnD3uaYo7KZXdsMWsZFFzxJVWU75i3ozsX/8xmy0cF1RUq6+fa9KV/chQ7tazj91Oc47j9e464xezb/B2pBfuVVM4mIDyXdBpwOrCo4dRDwtbR/O1A3gNd6KCIqgApJC8gG+GcDp0v6asqzLdmd1ZICdOpt3wAwfK/ODc5B3Nz6Daqk/6BKdt13JQD/ftQSxl4zYJ08h351Med/e0dO+Ok8OnYKOnbK/tMP23MV2+ywmvff6bTmP/M2O6wG4DNfWcKYawYC8OjdfbjkzncA2G34SlZXiA/L27NVv7WBuO/WVWy/y0e8+lw3Pn3U0mb9zK1VOwUrVnTg1DO+tN658RN2YvyE7KZfsTHohQu7sufu89cc9+u7kldezX5+5Yu7AlBZVcb4CTvyH8e83pwfIzfa6hDH5phmdyXZs+ndGshTX6CsKNivBtpLOgT4HHBQROwFvEQ21LHF6zOgin7brGbW9Gzs8eWnerDdsAref6fjmjzPPtaLbXfOvi1LFpVRnYao577XkfdndGTr7VbTb+tKZr7VmSWLsj8LX/x7D7Yd9hGQ9ZJffroHADOndWJ1RTt69a3igzkdqFiV/SdYtqSMqZO7MWSnwm+/bYiVqzoyb0F3Pv3J91JKMHSHxSWVnfzSNuy7z1y6d6uge7cK9t1nLpNf2gaAPr1XrqnvoANm8+7MrZq87XnT1IslbUma/W/8NBQxlixI1969/CfZ1JPbgeOBpzegyl5kC5GslLQrcGBTtrelnfrf73PZj7anqlJsvd1qzrpiJlecvS2z3+6UhiFWc3qawTHl2e7c9putKWsPZe2C0y+dTc/eWcQ+/sx5nP3VYbTvEAwYvHrNjceTRr3PlWdvy/1/7I+As6+YiZQF6z9etGP2rFPAsSd/wNCPf9RC34Utz7lnPcWeu8+nZ88Kbr/pfu64e08uu/zfOe2U5/jGf06hrH3wf09tz4x3G1wbB4Dlyztx15g9uOrybNbWnWP2XHPD8Jwz/0GvnhVIwTszenPV9Qc06+fKi7a6YL8imuev/MJpbpIGAjOAX0fEhWn5vdFAP+q5SVh3mp2kV4GjgLnAA2TPsL9JdqPxwoiYWHiTsJRpdsP36hyTHtu2oSyWM0ce/e2WboJtoPGTRr1QwuPX9eq964D47OhjS8p7/6eu36Rr5U2z9aALg2NEzAe6Fhy/C3y2kfIX1jkunG/4hXrK7FDs+ma2ZWuNwxelyNc0BjOzOpr4ScItigO0meWeA7SZWQ55HrSZWY611XnQDtBmlmsRUFX6gv2tigO0meWehzjMzHLIY9BmZjkWDtBmZvnkm4RmZjkU4TFoM7OcEtWexWFmlk8egzYzyyGvxWFmllex9pVfbY0DtJnlnmdxmJnlUPgmoZlZfrXVIY62+WvJzLYoESppK5WkMkkvSfprOu4jabykaelr74K850maLulNSUcUpO8naUo6d5UkpfROksak9OfSK/5qy4xM15gmaWRj7XSANrNci2j6AA38GHi94PhcYEJEDAMmpGMk7Ub2gutPAEcC10kqS2WuB04ChqXtyJR+ItmLrXcGrgAuS3X1AUYBBwD7A6MKfxEU4wBtZrlXEyppK4WkIcCXgBsLko8Gbk37twLHFKTfExEVETEDmA7sL2kQ0DMinonszdu31SlTW9d9wGGpd30EMD4iyiNiMTCetUG9KI9Bm1nubcAYdD9JkwuOb4iIG+rkuRI4B+hRkDYwIuZm14q5kgak9MHAswX5Zqe0yrRfN722zKxUV5WkpUDfwvQiZYpygDazXAtETemzOBZGxPD6Tko6ClgQES9IOqSE+op1y6OB9I0tU5SHOMws96LErQSfAr4i6V3gHuCzku4A5qdhC9LXBSn/bGDbgvJDgDkpfUiR9HXKSGoP9ALKG6irXg7QZpZvTXiTMCLOi4ghEbED2c2/JyLiW8A4oHZWxUjgwbQ/DhiRZmYMJbsZOCkNhyyTdGAaXz6hTpnauo5N1wjgMeBwSb3TzcHDU1q9PMRhZvnX/POgLwXGSjoRmAkcBxARr0kaC0wFqoBTI6I6lTkFuAXoAjySNoCbgNslTSfrOY9IdZVLuhh4PuW7KCLKG2qUA7SZ5V5zrGYXEROBiWl/EXBYPfkuAS4pkj4Z2L1I+kekAF/k3GhgdKltrDdAS7qaBn5vRcTppV7EzGxjBVBT47U46prcwDkzs80jAC83uq6IuLXwWFK3iFjR/E0yM1uX1+Koh6SDJE0lPRYpaS9J1zV7y8zMajXhPLstSSnT7K4ke0RxEUBE/As4uBnbZGZWoLQpdq3xtVglzeKIiFlpoaZa1fXlNTNrcq2wd1yKUgL0LEmfBEJSR+B01l0Fysys+QREG53FUcoQx8nAqWSLerwP7J2Ozcw2E5W4tS6N9qAjYiFw/GZoi5lZcW10iKOUWRw7SvqLpA8kLZD0oKQdN0fjzMwAz+JowF3AWGAQsA1wL3B3czbKzGyN2gdVStlamVICtCLi9oioStsdtMrfVWaWV9lrrxrfWpuG1uLok3aflHQu2dqpAXwdeGgztM3MLNNGZ3E0dJPwBdZ9C8APCs4FcHFzNcrMrJBaYe+4FA2txTF0czbEzKyoVnoDsBQlPUkoaXdgN6BzbVpE3NZcjTIzW6t13gAsRaMBWtIo4BCyAP0w8AXgabLXjJuZNb822oMuZRbHsWRvGpgXEd8F9gI6NWurzMwK1ZS4tTKlDHGsiogaSVWSepK97dYPqpjZ5uEF+xs0WdJWwB/JZnYsByY1Z6PMzAp5Fkc9IuKHaff3kh4FekbEK83bLDOzAg7Q65K0b0PnIuLF5mmSmZlBwz3oyxs4F8Bnm7gtm91br3TliG32bulm2IZoN7WlW2AtwEMcdUTEoZuzIWZmRQV+1NvMLLfaaA+6lHnQZmYtSlHa1mg9UmdJkyT9S9Jrkn6Z0vtIGi9pWvrau6DMeZKmS3pT0hEF6ftJmpLOXaX04lZJnSSNSenPSdqhoMzIdI1pkkY21l4HaDPLv6ZbsL8C+GxE7EX2+r4jJR0InAtMiIhhwIR0jKTdgBHAJ4AjgesklaW6rgdOAoal7ciUfiKwOCJ2Bq4ALkt19QFGAQcA+wOjCn8RFFPKG1Uk6VuSLkjH20nav4RvhJlZ02iiAB2Z5emwQ9oCOBq4NaXfChyT9o8G7omIioiYAUwH9pc0iGzK8TMREWRLXxSWqa3rPuCw1Ls+AhgfEeURsRgYz9qgXlQpPejrgIOAb6TjZcC1JZQzM9tkpQ5vpCGOfpImF2wnrVefVCbpZbKnosdHxHPAwIiYC5C+DkjZBwOzCorPTmmD037d9HXKREQVsBTo20Bd9SrlJuEBEbGvpJfSBRdL6lhCOTOzplH6LI6FETG8oQwRUQ3snZ6Q/nNarbM+xS4cDaRvbJmiSulBV6YxlwCQ1J9WuSyJmeVVU90kLBQRS4CJZMMM89OwBenrgpRtNrBtQbEhwJyUPqRI+jplJLUHegHlDdRVr1IC9FXAn4EBki4hW2r0VyWUMzNrGk00Bi2pf+o5I6kL8DngDWAcUDurYiTwYNofB4xIMzOGkt0MnJSGQZZJOjCNL59Qp0xtXccCT6Rx6seAwyX1TjcHD09p9SplLY47Jb1AtuSogGMi4vXGypmZNYmN6B03YBBwaxoVaAeMjYi/SnoGGCvpRGAmcBxARLwmaSwwFagCTk1DJACnALcAXYBH0gZwE3C7pOlkPecRqa5ySRcDz6d8F0VEeUONLWXB/u2AlcBfCtMiYmZjZc3MmkQTBei00Ns+RdIXkXVCi5W5BLikSPpkYL3x64j4iBTgi5wbDYwutb2l3CR8iLUD3J2BocCbZPMCzcyandroXa9Shjj2KDxOq9z9oJ7sZmbWRDZ4LY6IeFHSvzVHY8zMimqja3GUMgZ9ZsFhO2Bf4INma5GZWaGmvUm4RSmlB92jYL+KbEz6T83THDOzIhyg15emonSPiJ9upvaYma3PAXpdktpHRFVDr74yM2tuwrM4iplENt78sqRxwL3AitqTEXF/M7fNzMxj0I3oAywiewdh7XzoABygzWzzcIBez4A0g+NV1l+JqY1+u8ysRbTRiNNQgC4DurMRS+SZmTUlD3Gsb25EXLTZWmJmVh8H6PW0zfecm1m+hGdxFFN0ZSczs83OPeh1NbZOqZnZ5uIxaDOzvHKANjPLoRJfZ9UaOUCbWa4JD3GYmeWWA7SZWV45QJuZ5ZQDtJlZDnk1OzOzHHOANjPLJz/qbWaWU211iKNdSzfAzKxBsQFbIyRtK+lJSa9Lek3Sj1N6H0njJU1LX3sXlDlP0nRJb0o6oiB9P0lT0rmrJCmld5I0JqU/J2mHgjIj0zWmSRrZWHsdoM0s/5ooQANVwFkR8XHgQOBUSbsB5wITImIYMCEdk86NAD4BHAlcl16mDXA9cBIwLG1HpvQTgcURsTNwBXBZqqsPMAo4ANgfGFX4i6AYB2gzy7XaJwlL2RoTEXMj4sW0vwx4HRgMHA3cmrLdChyT9o8G7omIioiYAUwH9pc0COgZEc9ERAC31SlTW9d9wGGpd30EMD4iyiNiMTCetUG9KI9Bm1nuqabkQeh+kiYXHN8QETcUrTMbetgHeA4YGBFzIQvikgakbIOBZwuKzU5plWm/bnptmVmpripJS4G+helFyhTlAG1m+bZhiyUtjIjhjWWS1B34E/CTiPgwDR8XzVpPixp6FeDGlCnKQxxmlntNNcQBIKkDWXC+MyLuT8nz07AF6euClD4b2Lag+BBgTkofUiR9nTKS2gO9gPIG6qqXA7SZ5V/TzeIQcBPwekT8tuDUOKB2VsVI4MGC9BFpZsZQspuBk9JwyDJJB6Y6T6hTprauY4En0jj1Y8Dhknqnm4OHp7R6eYjDzHKvCedBfwr4NjBF0ssp7WfApcBYSScCM4HjACLiNUljgalkM0BOjYjqVO4U4BagC/BI2iD7BXC7pOlkPecRqa5ySRcDz6d8FzX25ioHaDPLvyYK0BHxNPW/ELvoe1gj4hLgkiLpk4Hdi6R/RArwRc6NBkaX2l4HaDPLN7/V28wsn/xGFTOzPIu2GaEdoM0s99yDti1G/21W89PfzaT3gCqiBh6+oy8P3NSfHT+xitMvnU3HzjVUV4lrzhvCmy93BWDox1dx+mWz6dajmpoacdoXh1FZ4VmWzenWZ15l1Yp21FSL6ipx2pd25WfXvcOQnSoA6NazmhUflvHDIz5OWfvgjN+8x857rKSsLPjbfX0Zc+3WdOpcw8//8A7bbF9BTbV49m+9GP0/DT581vr4rd75I2l5RHQvOP4OMDwiftRAma8Au0XEpZIuBJZHxP82e2M3s+oqccNF2zB9Sle6dKvmmkff4sW/9+D7v5jDHb8dyOQne/Jvn/2QE38xh3OO3Zl2ZcE5V8/kN6dvxztTu9CjdxXVlfU+OWVN6JzjduHDxWv/m/3qhzuu2T/p/NmsWJatu3PwUYvp0DE4+XO70alzDTc8OZWJD/ZmycIO/OkPA/nXP3vQvkMNl90zjeGHLmXyk702+2dpSb5J2ApExDiySeKtWvmCDpQv6ADAqhVlzJremX6DKomAbj2yKZrdelZTPj/Ls99nljHj9c68M7ULAMsWt6of+xYqOPjLiznn68Oyo4DOXatpVxZ07FxDVaVYubyMio/a8a9/9gCgqrId017tSv9BlS3Z8BbhAL0FkfRl4BdAR2ARcHxEzC+ll93aDByymp12X8UbL3bl9xcM5ld3v8P/u2AuUnDGV7L//EN2rCBCXHLX2/TqW83/PbgV9143oJGabZMF/OquaRDw0J39eeTOfmtO7X7AchZ/0IE5MzoD8NRDvTno8KXc/eIUOnep4fe/HMKyJev+9+zWs4oDP7eUB25qYz+7wDcJc6hLwZM+AH1Y2zt+GjgwIkLS94FzgLNKqVTSSWRruNKZrk3X2hbQuWs159/4Lr+/YBtWLi9j5Mh5/GHUNjz98FYc/OUlnPnbWZz79Z0oax/svv8KTvviMCpWtePSMW8z7ZUuvPx0j5b+CK3aGV/dhfL5HenVt5JL757OrOmdePW57Ht+6NGLmfjg2qWAP7b3Cmpq4Jv77UH3XlVcfv9bvPRUD+bN7ARAu7LgvGvf5cHRA9aktSVt9SZhnu8SrYqIvWs34IKCc0OAxyRNAX5Ktph2SSLihogYHhHDO7Dl/kMvax+cf+O7PHF/b/7xyFYAfP64cp5+OBub/PtferHL3isB+GBuB155phsflrenYlU7nn+iJzvvsaqlmt5mlM/vCMDSRR34x6O92DX9PNqVBZ/6whL+7y9rA/Shx5QzeWJPqqvE0kUdmPp8d3bZc+Wa8z+5bCbvz+jEn9ta77lW0y3Yv0XJc4BuyNXANRGxB/ADoHMLt2czC868fBazpnXm/hv6r0ldNL8Dex60AoC9/305c2Zkv4BemNiDobt9RKcuNbQrC/Y8aDkz32pj37LNrFOXarp0q16zv9/By3j3zex7vu+nP2TW251ZOLfjmvwfzOnI3p9cBgSdulSz674rmPV29vMb+dM5dOtZze9HDVnvOm1BUy7Yv6XJ8xBHQ3oB76f9Rt/r1dp8Yv8VfO64xbwztTPXjX8TgJv/ZxBX/nQIp1w0h7KyYHVFO678afYfevnS9tz/h/5c/fBbRIhJT/Rg0oSeLfkRWr3e/asYdeM7AJSVBU8+0JvJE7O/bj7zlcVMfGDdNx2Nu6U/Z/32PW6Y8DoIHh/blxmvd6XfoNV888fzmDmtE9c++saavI/e3Y82I2JDFuxvVbbUAH0hcK+k98nedjC0ZZuzeb02qTtHbLNX0XM/OnKXoulP3N+bJ+5v8PVn1oTmzezEKYd/vOi5y8/cYb20j1aWccnJO66XvnBuR44Ysm9TN2/L0zbjc34DdOEc6HR8C9nSfkTEg6xde7W+PBc2bwvNbHNpjcMXpchtgDYzA7Les4c4zMxyqm3GZwdoM8s/D3GYmeWUZ3GYmeVRK30IpRQO0GaWa9mDKm0zQjtAm1n+eTU7M7N8cg/azCyPPAZtZpZXXovDzCy/PMRhZpZD0XZfebWlrgdtZm1JRGlbIySNlrRA0qsFaX0kjZc0LX3tXXDuPEnTJb0p6YiC9P0kTUnnrpKklN5J0piU/pykHQrKjEzXmCappGWSHaDNLP+a7o0qtwBH1kk7F5gQEcOACekYSbsBI8je2HQkcJ2kslTmerJX5w1LW22dJwKLI2Jn4ArgslRXH2AUcACwPzCq8BdBfRygzSz3VFNT0taYiPg7UF4n+Wjg1rR/K3BMQfo9EVERETOA6cD+kgYBPSPimYgI4LY6ZWrrug84LPWujwDGR0R5RCwGxrP+L4r1eAzazPIt2JAHVfpJmlxwfENE3NBImYERMRcgIuZKqn3x42CyF4LUmp3SKtN+3fTaMrNSXVWSlgJ9C9OLlKmXA7SZ5ZqIDXlQZWFEDG+yS68vGkjf2DL18hCHmeVfE90krMf8NGxB+rogpc8Gti3INwSYk9KHFElfp4yk9mTvTy1voK4GOUCbWf41b4Aex9qXT49k7ev0xgEj0syMoWQ3Ayel4ZBlkg5M48sn1ClTW9exwBNpnPox4HBJvdPNwcNTWoM8xGFm+bZhY9ANknQ3cAjZWPVsspkVlwJjJZ0IzASOA4iI1ySNBaYCVcCpEVGdqjqFbEZIF+CRtAHcBNwuaTpZz3lEqqtc0sXA8ynfRRFR92blehygzSz3SpmhUYqI+EY9pw6rJ/8lwCVF0icDuxdJ/4gU4IucGw2MLrmxOECbWe5t0vDFFs0B2szyLXCANjPLrTa6FocDtJnlnhfsNzPLKwdoM7McioDqtjnG4QBtZvnnHrSZWU45QJuZ5VAAfiehmVkeBYTHoM3M8ifwTUIzs9zyGLSZWU45QJuZ5ZEXSzIzy6cAmmi50S2NA7SZ5Z970GZmeeRHvc3M8ikgPA/azCyn/CShmVlOeQzazCyHIjyLw8wst9yDNjPLoyCqq1u6ES3CAdrM8s3LjZqZ5Zin2ZmZ5U8A4R60mVkOhRfsNzPLrbZ6k1DRRqevAEj6AHivpdvRDPoBC1u6EbZBWvPPbPuI6L+xhSU9Svb9KcXCiDhyY6+VN206QLdWkiZHxPCWboeVzj8zK6ZdSzfAzMyKc4A2M8spB+jW6YaWboBtMP/MbD0egzYzyyn3oM3McsoB2swspxygc0hSSLq84PhsSRduQPkLJZ29gdd8WNJWaX/5hpS1htX9fkr6jqRrGinzFUnnpv0N/nla6+AAnU8VwNcklTo5f5NFxBcjYsnmup41LCLGRcSlLd0Oa1kO0PlURXZX/4y6JyRtL2mCpFfS1+3qqWM3SRMlvSPp9ILyD0h6QdJrkk4qSH93c/5CsIykL0t6TtJLkv4maWBKb7SXba2fA3R+XQscL6lXnfRrgNsiYk/gTuCqesrvChwB7A+MktQhpX8vIvYDhgOnS+rb9E23OrpIerl2Ay4qOPc0cGBE7APcA5zTEg20fPJiSTkVER9Kug04HVhVcOog4Gtp/3bg1/VU8VBEVAAVkhYAA4HZZEH5qynPtsAwYFFTt9/WsSoi9q49kPQdsl+QAEOAMZIGAR2BGZu9dZZb7kHn25XAiUC3BvLUN5G9omC/Gmgv6RDgc8BBEbEX8BLQeZNbaZviauCaiNgD+AH+eVgBB+gci4hyYCxZkK71T2BE2j+e7E/kUvUCFkfESkm7Agc2SUNtU/QC3k/7I1uyIZY/DtD5dznrLrV4OvBdSa8A3wZ+vAF1PUrWk34FuBh4tslaaRvrQuBeSU/RepcbtY3kR73NzHLKPWgzs5xygDYzyykHaDOznHKANjPLKQdoM7OccoC2BkmqTo8ovyrpXkldN6GuWyQdm/ZvlLRbA3kPkfTJjbhG0TVFSllrZENX8fMqc9bcHKCtMasiYu+I2B1YDZxceFJS2cZUGhHfj4ipDWQ5BNjgAG3WmjhA24Z4Ctg59W6flHQXMEVSmaTfSHo+rbL3AwBlrpE0VdJDwIDaitJKe8PT/pGSXpT0r7RC3w5kvwjOSL33T0vqL+lP6RrPS/pUKttX0uNpNbg/AGrsQ9S3ol86d3lqywRJ/VPaTpIeTWWeSk9hmjU7L5ZkJZHUHvgC2dOIkK2St3tEzEhBbmlE/JukTsA/JD0O7AN8DNiDbLGmqcDoOvX2B/4IHJzq6hMR5ZJ+DyyPiP9N+e4CroiIp9MSq48BHwdGAU9HxEWSvgSsE3Dr8b10jS7A85L+FBGLyNY8eTEizpJ0Qar7R2RLv54cEdMkHQBcB3x2I76NZhvEAdoa0yUtkQlZD/omsqGHSRFRu/La4cCetePLZOtLDAMOBu6OiGpgjqQnitR/IPD32rrS+iPFfI5sjeva456SeqRrfC2VfUjS4hI+U30r+tUAY1L6HcD9krqnz3tvwbU7lXANs03mAG2NWWepTIAUqFYUJgGnRcRjdfJ9kfpX2yssW8p6A+3IVuErXHq1ti0lr1dQZ0W/lZImUv8KcpGuu6Tu98Bsc/AYtDWFx4BTal8KIGkXSd2AvwMj0hj1IODQImWfAT4jaWgq2yelLwN6FOR7nGy4gZRv77T7d7JV/ZD0BaB3I21taEW/dkDtXwHfJBs6+RCYIem4dA1J2quRa5g1CQdoawo3ko0vvyjpVeAPZH+d/RmYBkwBrgf+r27BiPiAbNz4fkn/Yu0Qw1+Ar9beJCRbxW94ugk5lbWzSX4JHCzpRbKhlpmNtLWhFf1WAJ+Q9ALZGHPtm0+OB05M7XsNOLqE74nZJvNqdmZmOeUetJlZTjlAm5nllAO0mVlOOUCbmeWUA7SZWU45QJuZ5ZQDtJlZTv1/fcir6samoWMAAAAASUVORK5CYII=\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": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Probability of Detection (POD): 95.28675016479895\n",
      "False Alarm Rate (FAR): 13.738925417920916\n",
      "Positive Predictive Value (PPV): 0.05400709882308986\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 significant improvement of the results when using the balanced accuracy as target."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Save best model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "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": 25,
   "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": 26,
   "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": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcsAAAEICAYAAAAwft9dAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgTklEQVR4nO3de5xVdb3/8dcbNVBUvFDHcUinFE2MRCWvZfSQk5dK9Jd5+5WSpt3MjnnD9HRILSEsyezniayfWUpHTYjSBO04XQS0QWFA0gRBBTUUbY6EEsLn/LG+E5vdnll7ZvbM3jO8n4/HPGavy/e7PmuJ857vWt+9RxGBmZmZta1ftQswMzOrdQ5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LM6sYSV+RdHO16zCrNPl9lma1QdJy4F+ADQWr94mI57vY56cj4oGuVdf7SBoP7B0Rn6h2Ldb7eWRpVls+GhHbF3x1OigrQdLW1Tx+Z/XWuq12OSzNapykQZJ+KOkFSSslXSNpq7RtL0n/LWm1pJcl3SZpp7TtJ8AewC8lrZF0qaRRklYU9b9c0uj0erykuyT9VNL/AGPbO36JWsdL+ml63SApJH1K0nOSXpX0WUnvldQs6a+SbixoO1bSQ5K+K6lF0hOSji7YvrukGZJekbRE0rlFxy2s+7PAV4BT07kvSPt9StKfJL0m6WlJnynoY5SkFZIukrQqne+nCrZvK+lbkp5J9f1B0rZp22GSZqdzWiBpVCf+U1sNc1ia1b4fA28CewMHAh8CPp22CbgW2B3YD3g7MB4gIj4JPMum0eo3yzzeGOAuYCfgtpzjl+NQYChwKjAZuAIYDewPnCLpA0X7Pg0MBv4DuFvSLmnbVGBFOteTgW8UhmlR3T8EvgH8Vzr3A9I+q4CPADsCnwKul3RQQR+7AYOAeuAc4HuSdk7brgMOBo4AdgEuBTZKqgfuAa5J6y8Gfi7prR24RlbjHJZmtWV6Gp38VdJ0Sf8CHAf8W0T8LSJWAdcDpwFExJKIuD8i1kXES8C3gQ+03X1Z5kTE9IjYSBYqbR6/TFdHxBsRMQv4GzA1IlZFxErg92QB3GoVMDki1kfEfwFPAh+W9HbgfcBlqa/5wM3AJ0vVHRGvlyokIu6JiKWR+S0wC3h/wS7rgavS8e8F1gD7SuoHnA18KSJWRsSGiJgdEeuATwD3RsS96dj3A03A8R24RlbjfF/frLacWDgZR9IhwDbAC5JaV/cDnkvb3wbcQPYDf4e07dUu1vBcwes92zt+mf5S8Pr1EsvbFyyvjM1nHT5DNpLcHXglIl4r2jayjbpLknQc2Yh1H7Lz2A5YWLDL6oh4s2B5bapvMDAAWFqi2z2Bj0v6aMG6bYAH8+qx3sNhaVbbngPWAYOLfoi3uhYI4D0RsVrSicCNBduLp7v/jSwgAEjPHotvFxa2yTt+pdVLUkFg7gHMAJ4HdpG0Q0Fg7gGsLGhbfK6bLUvqD/wcOBP4RUSslzSd7FZ2npeBN4C9gAVF254DfhIR5/5TK+szfBvWrIZFxAtktwq/JWlHSf3SpJ7WW607kN0q/Gt6dnZJURd/Ad5ZsPxnYICkD0vaBrgS6N+F41fa24ALJG0j6eNkz2HvjYjngNnAtZIGSHoP2TPF29rp6y9AQ7qFCvAWsnN9CXgzjTI/VE5R6Zb0j4Bvp4lGW0k6PAXwT4GPSjomrR+QJgsN6fjpW61yWJrVvjPJftAvJrvFehdQl7Z9DTgIaCGbZHJ3UdtrgSvTM9CLI6IF+DzZ876VZCPNFbSvveNX2sNkk4FeBr4OnBwRq9O204EGslHmNOA/0vPBttyZvq+W9GgakV4A3EF2HmeQjVrLdTHZLds/Aq8AE4F+KcjHkM2+fYlspHkJ/vnap/hDCcysJkgaS/YBCu+rdi1mxfybj5mZWQ6HpZmZWQ7fhjUzM8vhkaWZmVkOv8+yDxo8eHA0NDRUuwwzs15l3rx5L0dEyY8pdFj2QQ0NDTQ1NVW7DDOzXkXSM21t821YMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHP5Qgj5o4coWGsbdU+0yat7yCR+udglm1kt4ZGlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpYjNywlXS/p3wqWZ0q6uWD5W5JC0r5F7SZLurSdfi+XtETSk5KOKVj/FklTJP1Z0hOSPpbWf1nSYknNkn4jac82+h0v6eL0+hZJyyTNl7RA0tE553p+qikkDS5YL0k3pG3Nkg4q2HahpMclLZI0VdKAdvofK2n39mrIqe8rnW1rZmadV84Hqc8GPg5MltQPGAzsWLD9COA+4DTgawBpv5OBI0t1KGlY2n9/YHfgAUn7RMQG4ApgVUTsk/rZJTV7DBgZEWslfQ74JnBqGfVfEhF3SfogMAUY2s6+DwG/AhqL1h+X2g0FDgVuAg6VVA9cAAyLiNcl3ZHO65Y2+h8LLAKeL6PuUr4CfKOTbbc4L94+rt3to+ZOanNbY2Njhasxs96snNuwD5EFImThtgh4TdLOkvoD+wGXkYVEq6OA5RHxTBt9jgF+FhHrImIZsAQ4JG07G7gWICI2RsTL6fWDEbE27TMXGFLOCRaYA9S3t0NEPBYRy9uo99bIzAV2klSXtm0NbCtpa2A72ghCSScDI4Hb0kh3W0kHS/qtpHlpxF4naVAabe+b2k2VdK6kCek48yXdVqL/8yQ1SWrasLalvCtiZmZlyR1ZRsTzkt6UtAdZaLaGzuFAC9AcEc2SNko6ICIWkAXn1Ha6rScLvFYrgHpJO6XlqyWNApYC50fEX4ranwP8Oq/2IscC0zvYplU98FzB8gqgPiKaJF0HPAu8DsyKiFmlOkij2/OBi1O7bYDvAmMi4iVJpwJfj4iz0363SPoOsHNE/ACy28QRMaKN/qeQjZzpXzc0OnmefcpuZ0xod3uj/0SXmZWp3Ak+raPL1rCcU7A8O+0zFTgtjbDGAHe2059KrAuy8B4CPBQRB6XjXLdZQ+kTZCO0tu+hbW6SpKeBn9L5W5gl65W0M9m5voPsdvLAVF859gXeDdwvaT5wJWm0HBH3AwuB7wGf7mTNZmZWIeWG5WyyYBxOdht2LtnI8giyIIUsLE8BRpONNle1098K4O0Fy0PIbl+uBtYC09L6O4HCyTSjyZ5pnhAR69K6r6dbk/PbONYlwN5kYfTjMs61I/WOBpZFxEsRsR64m023rPMIeDwiRqSv4RHxIfjHM9/9yEaru7TXiZmZdb+OjCw/ArwSERsi4hVgJ7LAnAMQEUvJwm4C7d+CBZhBNgrtL+kdZBNnHomIAH4JjEr7HQ0sBpB0IPB9sqD8RxBHxBWtgdPWwSJiI/AdoF/hzNsOmAGcmWbFHga0RMQLZLdfD5O0nSSlev/UTj+vATuk108Cb5V0eDq/bSTtn7ZdmPo5HfhRumULsL7gtZmZ9ZByw3Ih2SzYuUXrWlon4CRTgXexaWRYUkQ8DtxBFoT3AV9IM2Ehmyw0XlIz8EngorR+ErA9cGcaSc4os/bWYwZwDdDe21kukLSCbOTYrE1vkbkXeJpsItIPgM+nPh8G7gIeJbse/UjPDdtwC/CfaRS8FdmM4YmSFgDzgSMk7UN26/WiiPg98DuyUTGp7+ZSE3zMzKz7KMsQ60v61w2NurMmV7uMmrfcE3zMrICkeRExstQ2f4KPmZlZjnI+lKDT0vPBiUWrl0XESd153DySppHNYC10WUTMrFD/3+OfP5DhOxHx/yvRv5mZ9axuDcsUPhUJoErq7rCOiC90Z/9mZtazfBvWzMwsh8PSzMwsh8PSzMwsh8PSzMwsR7dO8LHqGF4/iCa/h9DMrGI8sjQzM8vhsDQzM8vhsDQzM8vhsDQzM8vhsDQzM8vh2bB90MKVLTSMu6faZVgR/5UTs97LI0szM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McDkszM7McfTIsJa1J3xskvS5pvqTFkm6VtE077XaV9KCkNZJuLNp2sKSFkpZIukGS0vo9UpvHJDVLOr6d/hskndGF8xol6YjOtjczs87ZEj5IfWlEjJC0FXA/cApwWxv7vgH8O/Du9FXoJuA8YC5wL3As8GvgSuCOiLhJ0rC0raGN/huAM4DbO3kuo4A1wOxOtrcKevH2cR3af9TcSWXv29jY2MFqzKw79cmRZSkRsQF4BKhvZ5+/RcQfyELzHyTVATtGxJyICOBW4MTWZsCO6fUg4Pl2ypgAvD+NdC+UtJWkSZL+mEaln0nH+7KkH6XXwyUtSkH8WeDC1P79RTWeJ6lJUtOGtS1lXRMzMyvPljCyBEDSAOBQ4EudaF4PrChYXsGm0B0PzJL0RWAgMLqdfsYBF0fER1JN5wEtEfFeSf2BhyTNAiYDjZJOAq4APhMRiyX9J7AmIq4r7jgipgBTAPrXDY1OnKN10G5nTOjQ/o3+E11mvdaWMLLcS9J8YDXwbEQ0d6IPlVjXGkinA7dExBDgeOAnksq9rh8Czkz1PQzsCgyNiI3AWOAnwG8j4qFO1GxmZhWyJYwsW59Z1pGN1k6IiBkd7GMFMKRgeQibbreeQ/b8koiYk0awg4FVZfQr4IsRMbPEtqFkzyd372CtZmZWYVvCyBKAiHiB7Dbo5Z1s+5qkw9Is2DOBX6TNzwJHA0jaDxgAvNRGV68BOxQszwQ+1zpDV9I+kgZKGgR8BzgK2FXSyW20NzOzHrDFhGUyHdiueHJMIUnLgW8DYyWtSBNrAD4H3AwsAZaSzYQFuAg4V9ICYCowNk0CKqUZeFPSAkkXpv4WA49KWgR8n2y0fz3w/yLiz2Qj1wmS3gb8Ejip1AQfMzPrPmr757r1Vv3rhkbdWZOrXYYVWe4JPmY1TdK8iBhZatuWNrI0MzPrsC1hgs8/kXQMMLFo9bKIOKlC/Q8nm8laaF1EHFqJ/s3MrGdtkWGZZp+WmoFaqf4XAiO6q38zM+tZvg1rZmaWw2FpZmaWw2FpZmaWw2FpZmaWY4uc4NPXDa8fRJPf02dmVjEeWZqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeXwbNg+aOHKFhrG3VPtMqyP8V9NsS2ZR5ZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5HJZmZmY5aiosJa0pWh4r6cb0eryklZLmS1os6fS0fpmkfYvaTZZ0ac9VXh5JIyQd34X2J0oaVsmazMwsX2/7IPXrI+I6SUOBeZLuAn4GnAZ8DUBSP+Bk4MjuKEDSVhGxoZPNRwAjgXs72f5E4FfA4k62ty3Ui7eP63Ifo+ZO6lL7xsbGLtdgVi01NbIsV0Q8BawFdgamkoVlq6OA5RHxTKm2abR6t6T7JD0l6ZsF206XtFDSIkkTC9avkXSVpIeBw9PyREnzJD0g6RBJjZKelnRCG8d9C3AVcGoaHZ8qaaCkH0n6o6THJI1J+94g6avp9TGSfifpCOAEYFJqv1dR/+dJapLUtGFtS4evqZmZta3WRpbbSppfsLwLMKN4J0kHAU9FxCpglaSNkg6IiAVkwTk15zgjgAOBdcCTkr4LbAAmAgcDrwKzJJ0YEdOBgcCiiGgNsIFAY0RcJmkacA3wr8Aw4Melao6Iv6cAHBkR56d+vgH8d0ScLWkn4BFJDwDjgD9K+j1wA3B8RCyVNAP4VUTcVaL/KcAUgP51QyPn/G0Ls9sZE7rcR6P/RJdtwWotLF+PiBGtC5LGkt22bHWhpHOBdwLHFqyfCpwm6XFgDPDVnOP8JiJa0jEWA3sCu5IF4Etp/W1ko9TpZEH684L2fwfuS68XAusiYr2khUBDmecK8CHgBEkXp+UBwB4R8ad0nr8DLoyIpR3o08zMKqzWwjJP6zPL/wPcKmmviHiDLCxnAb8FmtOIsz3rCl5vILsOamf/N4qeU66PiNbR28bW/iJio6SOXFMBH4uIJ0tsGw6sBnbvQH9mZtYNeuszy7uBJuCstLyULFgmkH8Lti0PAx+QNFjSVsDpZOFbSa8BOxQszwS+KEkAkg5M3/cELiK7VXycpEPbaG9mZj2gV4ZlchXw5TT7FbKQfBcwrTOdRcQLwOXAg8AC4NGI+EUlCi3wIDCsdYIPcDWwDdAsaRFwdQrOHwIXR8TzwDnAzZIGkM38vSRNBtqrjWOYmVmFadPdROsr+tcNjbqzJle7DOtjlnuCj/VxkuZFxMhS23rzyNLMzKxH9LYJPmWTdAzZW0EKLYuIk/rysc3MrPL6bFhGxEyyCTRb1LHNzKzyfBvWzMwsh8PSzMwsh8PSzMwsh8PSzMwsR5+d4LMlG14/iCa/J87MrGI8sjQzM8vhsDQzM8vhsDQzM8vhsDQzM8vhsDQzM8vh2bB90MKVLTSMu6faZZhtxn+1xHozjyzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxyOCzNzMxybFFhKWlN0fJYSTem1+MlrZQ0X9JiSaen9csk7VvUbrKkS3uucjMzqyZ/kPrmro+I6yQNBeZJugv4GXAa8DUASf2Ak4Eju6MASVtFxIbu6Ntq04u3j6t2CT1i1NxJ1S6h2zU2Nla7BOsmW9TIslwR8RSwFtgZmEoWlq2OApZHxDOl2qbR6t2S7pP0lKRvFmw7XdJCSYskTSxYv0bSVZIeBg5PyxMlzZP0gKRDJDVKelrSCW0c9zxJTZKaNqxtqcBVMDOzVlvayHJbSfMLlncBZhTvJOkg4KmIWAWskrRR0gERsYAsOKfmHGcEcCCwDnhS0neBDcBE4GDgVWCWpBMjYjowEFgUEV9Nxx8INEbEZZKmAdcA/woMA35cquaImAJMAehfNzTKuBZWI3Y7Y0K1S+gRjf4TXdaLbWlh+XpEjGhdkDQWGFmw/UJJ5wLvBI4tWD8VOE3S48AY4Ks5x/lNRLSkYywG9gR2JQvAl9L628hGqdPJgvTnBe3/DtyXXi8E1kXEekkLgYYyz9XMzCrEt2E3d31E7AucCtwqaUBaPxU4BRgNNKcRZ3vWFbzeQPZLidrZ/42i55TrI6J1dLixtb+I2MiW9wuOmVnVOSxLiIi7gSbgrLS8FFgNTCD/FmxbHgY+IGmwpK2A04HfVqBcMzPrZg7Ltl0FfDnNfoUsJN8FTOtMZxHxAnA58CCwAHg0In5RiULNzKx7adPdPusr+tcNjbqzJle7DLPNLPcEH6txkuZFxMhS2zyyNDMzy+HJIp0k6Riyt4IUWhYRJ1WjHjMz6z4Oy06KiJnAzGrXYWZm3c+3Yc3MzHI4LM3MzHI4LM3MzHI4LM3MzHJ4gk8fNLx+EE1+T5uZWcV4ZGlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDYWlmZpbDbx3pgxaubKFh3D3VLsOsQ/wnvKyWeWRpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2FpZmaWw2HZCZLWpO8Nkl6XNF/SYkm3StqmnXa7SnpQ0hpJNxZtO1jSQklLJN0gSWn9HqnNY5KaJR3fvWdnZmbF/NmwXbc0IkZI2gq4HzgFuK2Nfd8A/h14d/oqdBNwHjAXuBc4Fvg1cCVwR0TcJGlY2tZQ6ZOw6nrx9nHVLqHqRs2dVO0SqqqxsbHaJVg7PLKskIjYADwC1Lezz98i4g9kofkPkuqAHSNiTkQEcCtwYmszYMf0ehDwfKm+JZ0nqUlS04a1LV06FzMz25xHlhUiaQBwKPClTjSvB1YULK9gU+iOB2ZJ+iIwEBhdqoOImAJMAehfNzQ6UYNV0W5nTKh2CVXX6L86YjXMI8uu20vSfGA18GxENHeiD5VY1xp4pwO3RMQQ4HjgJ5L8383MrAf5h27XLY2IEcDewGGSTuhEHyuAIQXLQ9h0u/Uc4A6AiJgDDAAGd7paMzPrMIdlhUTEC8A44PJOtn1N0mFpFuyZwC/S5meBowEk7UcWli9VpGgzMyuLw7KypgPbSXp/WztIWg58GxgraUWa4QrwOeBmYAmwlGwmLMBFwLmSFgBTgbFpEpCZmfUQT/DphIjYPn1fTsFbQFKIHZDTtqGN9U3889tJiIjFwJGdr9bMzLrKI0szM7McHll2A0nHABOLVi+LiJOqUY+ZmXWNw7IbRMRMYGa16zAzs8rwbVgzM7McDkszM7McDkszM7McDkszM7McnuDTBw2vH0STP5TazKxiPLI0MzPL4bA0MzPL4bA0MzPL4bA0MzPL4bA0MzPL4bA0MzPL4beO9EELV7bQMO6eapdhZh2w3G/3qmkeWZqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeVwWJqZmeWoelhKul7SvxUsz5R0c8HytySFpH2L2k2WdGk7/V4uaYmkJyUdU7D+LZKmSPqzpCckfSyt/7KkxZKaJf1G0p5t9Dte0sXp9S2SlkmaL2mBpKNzzvX8VFNIGlywXpJuSNuaJR1UsO1CSY9LWiRpqqQB7R3DzMwqrxY+G3Y28HFgsqR+wGBgx4LtRwD3AacBXwNI+50MHFmqQ0nD0v77A7sDD0jaJyI2AFcAqyJin9TPLqnZY8DIiFgr6XPAN4FTy6j/koi4S9IHgSnA0Hb2fQj4FdBYtP641G4ocChwE3CopHrgAmBYRLwu6Y50XreUUZdZj3nx9nHVLqHXGzV3UrVL6PUaGxu7re+qjyzJAuSI9Hp/YBHwmqSdJfUH9gMuIwuJVkcByyPimTb6HAP8LCLWRcQyYAlwSNp2NnAtQERsjIiX0+sHI2Jt2mcuMKSD5zEHqG9vh4h4LCKWt1HvrZGZC+wkqS5t2xrYVtLWwHbA86X6lnSepCZJTRvWtnSwdDMza0/VR5YR8bykNyXtQRaaraFzONACNEdEs6SNkg6IiAVkwTm1nW7ryQKv1QqgXtJOaflqSaOApcD5EfGXovbnAL/u4KkcC0zvYJtW9cBzBcsrgPqIaJJ0HfAs8DowKyJmleogIqaQjWzpXzc0OlmHWafsdsaEapfQ6zX6r47UtFoYWcKm0WVrWM4pWJ6d9pkKnJZGWGOAO9vpTyXWBdkvB0OAhyLioHSc6zZrKH0CGAmUe09kkqSngZ8C3yizTVn1StqZ7FzfQXY7eWCqz8zMelCthOVssmAcTnYbdi7ZyPIIsiCFLCxPAUaTjTZXtdPfCuDtBctDyG5frgbWAtPS+juBwsk0o8meaZ4QEevSuq+nCTzz2zjWJcDewJXAj8s4147UOxpYFhEvRcR64G423bI2M7MeUith+RDwEeCViNgQEa8AO5EF5hyAiFhKFnYTaP8WLMAMslFof0nvIJs480hEBPBLYFTa72hgMYCkA4HvkwXlP4I4Iq6IiBERMaKtg0XERuA7QL/CmbcdMAM4M82KPQxoiYgXyG6/HiZpO0lK9f6pE/2bmVkX1EpYLiSbBTu3aF1L6wScZCrwLjaNDEuKiMeBO8iC8D7gC2kmLGSThcZLagY+CVyU1k8CtgfuTCPJGR05gRTE1wDtvZ3lAkkryEaOzQVvkbkXeJpsItIPgM+nPh8G7gIeJbse/UjPJc3MrOco+xlvfUn/uqFRd9bkapdhZh2w3BN8qk7SvIgYWWpbrYwszczMalbV3zrSFen54MSi1csi4qRq1NNK0jSyGayFLouImdWox8zMuqZXh2UKn5oLoGqHtZmZVZZvw5qZmeVwWJqZmeVwWJqZmeVwWJqZmeXo1RN8rLTh9YNo8nu2zMwqxiNLMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHA5LMzOzHIqIatdgFSbpNeDJatfRQYOBl6tdRAe55p7hmntOb6y7kjXvGRFvLbXBH3fXNz0ZESOrXURHSGpyzd3PNfeM3lgz9M66e6pm34Y1MzPL4bA0MzPL4bDsm6ZUu4BOcM09wzX3jN5YM/TOunukZk/wMTMzy+GRpZmZWQ6HpZmZWQ6HZS8j6VhJT0paImlcie2SdEPa3izpoHLb1mjNyyUtlDRfUlMN1fwuSXMkrZN0cUfa1mjNtXqd/2/6N9EsabakA8ptW6M11+p1HpPqnS+pSdL7ym1bozVX/jpHhL96yRewFbAUeCfwFmABMKxon+OBXwMCDgMeLrdtrdWcti0HBtfgdX4b8F7g68DFHWlbazXX+HU+Atg5vT6ul/x7LllzjV/n7dk0h+U9wBO94DqXrLm7rrNHlr3LIcCSiHg6Iv4O/AwYU7TPGODWyMwFdpJUV2bbWqu5WnJrjohVEfFHYH1H29ZgzdVSTs2zI+LVtDgXGFJu2xqsuVrKqXlNpJQBBgJRbtsarLlbOCx7l3rguYLlFWldOfuU07Y7dKVmyP4HmCVpnqTzuq3K8uvpzrZd0dXj9obrfA7ZHYjOtK2UrtQMNXydJZ0k6QngHuDsjrTtBl2pGbrhOvvj7noXlVhX/NtUW/uU07Y7dKVmgCMj4nlJbwPul/RERPyuohX+s65cq1q+zu2p6ess6YNkwdP6XKrmr3OJmqGGr3NETAOmSToKuBoYXW7bbtCVmqEbrrNHlr3LCuDtBctDgOfL3Kectt2hKzUTEa3fVwHTyG7PdLeuXKtavs5tquXrLOk9wM3AmIhY3ZG23aArNdf0dW6VQmUvSYM72raCulJz91zn7n5Q66+KPvTeGngaeAebHnrvX7TPh9l8sswj5batwZoHAjsUvJ4NHFsLNRfsO57NJ/jU7HVup+aavc7AHsAS4IjOnm8N1VzL13lvNk2WOQhYmf5/rOXr3FbN3XKdu/WE/dUt/4iOB/5MNlPsirTus8Bn02sB30vbFwIj22tbyzWTzYRbkL4er7GadyP77fd/gL+m1zvW+HUuWXONX+ebgVeB+emrqRf8ey5Zc41f58tSTfOBOcD7esF1Lllzd11nf9ydmZlZDj+zNDMzy+GwNDMzy+GwNDMzy+GwNDMzy+GwNDMzy+GwNDMzy+GwNDMzy/G/icsUxAEGy9MAAAAASUVORK5CYII=\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
}