HyperparameterTuning.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "60afa7d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mp\n",
    "\n",
    "from sklearn.model_selection import train_test_split \n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.ensemble import GradientBoostingRegressor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "a2487c21",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Getting the well (feature) and streamflow (target) data \n",
    "ibound = np.loadtxt('base/ibound_lyr1.txt')            \n",
    "all_well_loc = np.asarray(np.where(ibound == 1)).T\n",
    "\n",
    "#Grab the pickles of ReachQ (streamflow) and well heads\n",
    "ReachQdf = pd.read_pickle('outputs/SS_reachQ1000_test11.pkl')\n",
    "well_head_df= pd.read_pickle('outputs/heads/WellHead_SS1000test11.pkl')\n",
    "\n",
    "#Calculating the difference in heads from the wells and in the flows from with the well and without\n",
    "NoPumpingReachQ = pd.read_pickle('outputs/SS_reachQ1000_test11NoPumping.pkl')\n",
    "No_pumping_well_head_df = pd.read_pickle('outputs/heads/WellHead_SS1000test11NoPumping.pkl')\n",
    "members = np.loadtxt(\"mf_notebooks/ModflowSfrWellOutput1000test11NoPumping.txt\")\n",
    "\n",
    "diff_head_df = No_pumping_well_head_df - well_head_df.loc[members.astype('int')]\n",
    "diff_reachq_df = NoPumpingReachQ - ReachQdf.loc[members.astype('int')]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "f4d65f3b",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Setting up the head inputs - dropping out the downstream most CH boundary and the river cells \n",
    "river_well_loc = np.where(all_well_loc[:,0] == 25)\n",
    "CHBound_loc = np.where(all_well_loc[:,1] > 47)\n",
    "drop_locs= np.concatenate((river_well_loc, CHBound_loc), axis = 1)\n",
    "\n",
    "#currently not dropping any\n",
    "drop_locs = [[]]\n",
    "\n",
    "well_head_noRiv = well_head_df.drop(drop_locs[0], axis = 1)\n",
    "streamflow = ReachQdf.loc[well_head_noRiv.index]\n",
    "\n",
    "\n",
    "well_head_noRiv = diff_head_df.drop(drop_locs[0], axis = 1)\n",
    "streamflow = diff_reachq_df.loc[well_head_noRiv.index]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f06bdff",
   "metadata": {},
   "source": [
    "### Decision Tree Hyperparameter Tuning \n",
    "\n",
    "Prints out the best hyperparameter combo from the search as a dictionary - which we then use in the other notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "21662f61",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'criterion': ['squared_error', 'friedman_mse'], 'max_depth': [4, 5, 6, 7, 8], 'min_samples_split': [2, 5, 10, 15], 'min_samples_leaf': [1, 2, 5, 10], 'max_leaf_nodes': [25, 39, 53, 67, None], 'min_impurity_decrease': [5, 15, 25, 35, 45]}\n",
      "Number of possible parameter combinations:  4000\n",
      "Fitting 3 folds for each of 100 candidates, totalling 300 fits\n",
      "Best DT pars:  {'min_samples_split': 2, 'min_samples_leaf': 2, 'min_impurity_decrease': 5, 'max_leaf_nodes': None, 'max_depth': 8, 'criterion': 'squared_error'}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAHFCAYAAAAT5Oa6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB3hUlEQVR4nO3dd3RU1drH8e+kF5IhvVBDEYl0pAREeq+velERFLEgqBgRBfUq4PWCYFcs2IiCAhakSi4oiCBEEIhSlKvSISGkh/Ry3j9yMxLSJiGQ9vuslbWYM/ucs+cQmce9n/1sk2EYBiIiIiJ1mE1Vd0BERESkqikgEhERkTpPAZGIiIjUeQqIREREpM5TQCQiIiJ1ngIiERERqfMUEImIiEidp4BIRERE6jwFRCIiIlLnKSASuYrCwsIwmUyWHycnJ/z9/enbty/z588nJibmit7/+PHjmEwmwsLCynXexIkTadq06RXpU2n3vPhZlfQzceLEq9qvS13aH3d3d3r06MHy5curtF9XU0V/r0SqE5O27hC5esLCwrj77rtZsmQJ1157LdnZ2cTExLBjxw6WLFmCra0tK1euZMCAAVfk/pmZmezfv5/mzZvj4+Nj9Xl//fUXycnJdOzY8Yr0q6R7nj9/3vJ63759PPjgg8ybN4++fftajvv4+NC8efOr1q9LmUwmbrnlFh577DEMw+DYsWPMmzePQ4cO8emnnzJu3Lgq69vVUtHfK5HqRAGRyFVUEBDt2bOH66+/vtB7J0+e5IYbbiAxMZE//vgDPz+/Kupl9fT999/Tt29fvvjiC2655ZYS26Wnp+Pk5ITJZLoq/TKZTDz44IMsWrTIcuzEiRM0bdqUG2+8kW3btl2VfhRIS0vDxcXlqt5TpDbQlJlINdG4cWNefvllUlJSWLx4caH3fv75Z0aNGoWnpydOTk507NiRzz//vMg1zpw5w/3330+jRo1wcHAgMDCQW265hXPnzgHFT22cP3/eco6joyM+Pj707NmTb7/91tKmuCmzjIwMnnzySYKCgnBwcKBBgwY8+OCDJCYmFmrXtGlTRowYQXh4OJ06dcLZ2Zlrr72Wjz766PIeGH9PQW7atIlJkybh4+ODi4sLmZmZAKxcuZKQkBBcXV2pV68egwcPZv/+/UWuY+3ztVaTJk3w8fGxPPcCycnJzJgxo9AzCw0NJTU1tVC7xMRE7rnnHjw9PalXrx7Dhw/n6NGjmEwm5syZY2k3Z84cTCYT+/bt45ZbbsHDw8MyWmYYBm+//TYdOnTA2dkZDw8PbrnlFo4ePVroXvv372fEiBH4+vri6OhIYGAgw4cP5/Tp05Y2X3zxBd26dcNsNuPi4kKzZs2YNGmS5f2Spsx27NhB//79cXNzw8XFhR49erBhw4ZCbQr+Drdu3cqUKVPw9vbGy8uLm266ibNnz5b72YtUlAIikWpk2LBh2Nra8sMPP1iObd26lZ49e5KYmMi7777LmjVr6NChA7feemuhL6AzZ87QpUsXvv76a6ZPn87GjRt57bXXMJvNJCQklHjPCRMmsHr1ap599lk2bdrEBx98wIABA4iLiyvxHMMwGDNmDC+99BITJkxgw4YNTJ8+nY8//ph+/fpZApICv/zyC4899hiPPvooa9asoV27dtxzzz2FPuflmDRpEvb29ixdupQvv/wSe3t75s2bx+23305wcDCff/45S5cuJSUlhV69enH48GHLudY+3/JISkoiPj6ea665xnIsLS2N3r178/HHHzNt2jQ2btzIzJkzCQsLY9SoURQM1ufl5TFy5Eg+++wzZs6cyddff023bt0YMmRIife76aabaNGiBV988QXvvvsuAJMnTyY0NJQBAwawevVq3n77bQ4dOkSPHj0sgVpqaioDBw7k3LlzvPXWW2zevJnXXnuNxo0bk5KSAsCuXbu49dZbadasGStWrGDDhg08++yz5OTklPoMtm3bRr9+/UhKSuLDDz9k+fLluLm5MXLkSFauXFmk/b333ou9vT2fffYZCxcu5Pvvv2f8+PHle/Ail8MQkatmyZIlBmDs2bOnxDZ+fn5G69atLa+vvfZao2PHjkZ2dnahdiNGjDACAgKM3NxcwzAMY9KkSYa9vb1x+PDhEq997NgxAzCWLFliOVavXj0jNDS01H7fddddRpMmTSyvw8PDDcBYuHBhoXYrV640AOO9996zHGvSpInh5ORknDhxwnIsPT3d8PT0NCZPnlzqfS+2detWAzC++OILy7GC53nnnXcWanvy5EnDzs7OePjhhwsdT0lJMfz9/Y2xY8dajln7fEsCGFOnTjWys7ONrKws47///a8xatQow83Nzfj5558t7ebPn2/Y2NgU+bv/8ssvDcD45ptvDMMwjA0bNhiA8c477xRqN3/+fAMwZs+ebTk2e/ZsAzCeffbZQm137dplAMbLL79c6PipU6cMZ2dn44knnjAMwzB+/vlnAzBWr15d4ud76aWXDMBITEwssU1xv1fdu3c3fH19jZSUFMuxnJwco02bNkbDhg2NvLw8wzD+/jucOnVqoWsuXLjQAIyoqKgS7ytSmTRCJFLNGBel9f3555/8/vvv3HHHHQDk5ORYfoYNG0ZUVBRHjhwBYOPGjfTt25fWrVuX635du3YlLCyM559/noiICLKzs8s8Z8uWLQBFVnj94x//wNXVle+++67Q8Q4dOtC4cWPLaycnJ6655hpOnDhRrr6W5Oabby70+j//+Q85OTnceeedhZ6Zk5MTvXv35vvvvwfK93xL8/bbb2Nvb4+DgwPXXHMNGzduZPny5XTu3NnSZv369bRp04YOHToUus/gwYMxmUyWPhXkHI0dO7bQPW6//XarP//69esxmUyMHz++0L38/f1p37695V4tWrTAw8ODmTNn8u677xYaOSvQpUsXS38+//xzzpw5U+bzSE1N5aeffuKWW26hXr16luO2trZMmDCB06dPF3muo0aNKvS6Xbt2AJX2OyJSFgVEItVIamoqcXFxBAYGAlimNmbMmIG9vX2hn6lTpwIQGxsL5OcCNWzYsNz3XLlyJXfddRcffPABISEheHp6cueddxIdHV3iOXFxcdjZ2RVZUWQymfD39y8y3ebl5VXkGo6OjqSnp5e7v8UJCAgo9LrguXXp0qXIc1u5cqXlmZXn+ZZm7Nix7Nmzh507d7J48WLc3Ny47bbb+OOPPwr16ddffy1yHzc3NwzDsNyn4Nl6enoWukdpSfbFfX7DMPDz8ytyv4iICMu9zGYz27Zto0OHDjz11FNcd911BAYGMnv2bEtgfOONN7J69WpLgNmwYUPatGlTalmBhIQEDMMo0i/A8rtd1u+Io6MjQKX9joiUxa6qOyAif9uwYQO5ubn06dMHAG9vbwCefPJJbrrppmLPadWqFZC//PziRFhreXt789prr/Haa69x8uRJ1q5dy6xZs4iJiSE8PLzYc7y8vMjJyeH8+fOFgiLDMIiOjraMKlwtl64oK3huX375JU2aNCnxvPI839L4+PhYVg2GhITQunVrevfuzaOPPsr69est93J2di4xmbygLwXPNj4+vlBQVFqAWtznN5lMbN++3RJYXOziY23btmXFihUYhsGvv/5KWFgYzz33HM7OzsyaNQuA0aNHM3r0aDIzM4mIiGD+/PmMGzeOpk2bEhISUuT6Hh4e2NjYEBUVVeS9gkTpgs8rUl0oIBKpJk6ePMmMGTMwm81MnjwZyP8ybtmyJb/88gvz5s0r9fyhQ4eydOlSjhw5YtWXeHEaN27MQw89xHfffcePP/5YYrv+/fuzcOFCli1bxqOPPmo5/tVXX5Gamkr//v0rdP/KMnjwYOzs7Pjrr7+KTCddrDzPtzx69erFnXfeyccff8yuXbsICQlhxIgRzJs3Dy8vL4KCgko8t3fv3ixcuJCVK1cyZcoUy/EVK1ZYff8RI0bwwgsvcObMmSJTbyUxmUy0b9+eV199lbCwMPbt21ekjaOjI71796Z+/fr85z//Yf/+/cUGRK6urnTr1o1Vq1bx0ksv4ezsDOQnjC9btoyGDRsWSjgXqQ4UEIlUgYMHD1ryOmJiYti+fbulMOPXX39daNRl8eLFDB06lMGDBzNx4kQaNGhAfHw8v/32G/v27eOLL74A4LnnnmPjxo3ceOONPPXUU7Rt25bExETCw8OZPn061157bZF+JCUl0bdvX8aNG8e1116Lm5sbe/bsITw8vMQRE4CBAwcyePBgZs6cSXJyMj179uTXX39l9uzZdOzYkQkTJlT+QyuHpk2b8txzz/H0009z9OhRhgwZgoeHB+fOnWP37t24uroyd+5cwPrnW17/+te/WLlyJc888wzffvstoaGhfPXVV9x44408+uijtGvXjry8PE6ePMmmTZt47LHHLKvJevbsyWOPPUZycjKdO3dm165dfPLJJwDY2JSd6dCzZ0/uv/9+7r77bn7++WduvPFGXF1diYqKYseOHbRt25YpU6awfv163n77bcaMGUOzZs0wDINVq1aRmJjIwIEDAXj22Wc5ffo0/fv3p2HDhiQmJvL6669jb29P7969S+zD/PnzGThwIH379mXGjBk4ODjw9ttvc/DgQZYvX37V6kSJWK0KE7pF6pyCFTUFPw4ODoavr6/Ru3dvY968eUZMTEyx5/3yyy/G2LFjDV9fX8Pe3t7w9/c3+vXrZ7z77ruF2p06dcqYNGmS4e/vb9jb2xuBgYHG2LFjjXPnzhmGUXQ1UEZGhvHAAw8Y7dq1M9zd3Q1nZ2ejVatWxuzZs43U1FTLdS9dZWYY+SvFZs6caTRp0sSwt7c3AgICjClTphgJCQmF2jVp0sQYPnx4kc/Uu3dvo3fv3lY/u9JWmZW0am/16tVG3759DXd3d8PR0dFo0qSJccsttxjffvttoXbWPt/iAMaDDz5Y7HuPP/64ARjbtm0zDMMwLly4YPzzn/80WrVqZTg4OBhms9lo27at8eijjxrR0dGW8+Lj4427777bqF+/vuHi4mIMHDjQiIiIMADj9ddft7QrWGV2/vz5Yu//0UcfGd26dTNcXV0NZ2dno3nz5sadd95pWf32+++/G7fffrvRvHlzw9nZ2TCbzUbXrl2NsLAwyzXWr19vDB061GjQoIHl93XYsGHG9u3bLW2KW2VmGIaxfft2o1+/fpb7d+/e3Vi3bl2hNiX9HRb8fW/durWEJy9SuVSpWkSkBvjss8+44447+PHHH+nRo0dVd0ek1lFAJCJSzSxfvpwzZ87Qtm1bbGxsiIiI4MUXX6Rjx45XfSsQkbpCOUQiItWMm5sbK1as4Pnnnyc1NZWAgAAmTpzI888/X9VdE6m1NEIkIiIidZ4KM4qIiEidp4BIRERE6jwFRCIiIlLnKanaSnl5eZw9exY3NzcVFBMREakhDMMgJSWFwMDAUgubKiCy0tmzZ2nUqFFVd0NEREQq4NSpU6VugK2AyEpubm5A/gN1d3ev4t6IiIiINZKTk2nUqJHle7wkCoisVDBN5u7uroBIRESkhikr3UVJ1SIiIlLnKSASERGROk8BkYiIiNR5CohERESkzlNAJCIiInWeAiIRERGp8xQQiYiISJ2ngEhERETqPAVEIiIiUucpIBIREZE6TwGRiIiI1HkKiERERKTOU0AkIiIiVSojI4NNmzZVaR8UEImIiEiV+eOPP+jRowfDhg1jx44dVdYPBUQiIiJSJVasWEGnTp3Yv38/Hh4eZGRkVFlfFBCJiIjIVZWens7kyZO5/fbbuXDhAr169SIyMpIBAwZUWZ8UEImIiMhV8/vvv9OtWzfee+89TCYT//znP9myZQsNGjSo0n7ZVendRUREpM745JNPmDJlCmlpafj6+rJs2TIGDhxY1d0CFBCJiIhICXLzDCKOxrHrrzgMw8Dd2Y7kjBxMmAhp7kX3Zl7Y2pjy2/0Vx66jsXDJewCpqak89NBDhIWFAdCnb1+m//tN0pzrs+uvOLoGeVraVhWTYRhGlfaghkhOTsZsNpOUlIS7u3tVd0dEROSKCj8YxaxVB0hMyy6xTX0Xe27p2JBlu0+QkZ1X5L0XbmpLI1M8Y8eO5fDhw9jY2DBo/EPEtxzOuQt/X9fT1Z7/69CAAcH+lR4cWfv9rYDISgqIRESkrgg/GMUDy/Zd1jUMwyD1wGaSt7xHdmYGtq4eeI96HKfG7Uo9L8DsxOyRwQxpE3BZ9y9g7fe3kqpFRETEIjfPYM7aw5d1jbysdOLWv0zcxjfIzszAqWlHAu5+s8xgCCA6KYMpy/YRfjDqsvpQXsohEhEREYvdx+KJTq54PaCsmKOcX7OAnPgzYLKhfq/xuHe/BZPJujEYAzABc9cdZmCw/1XLLVJAJCIiIhYxKRULhgzD4MIv4cR/+x7kZmNbzwvv0U/g1PC68l8LiErKYPexeEKae1WoP+WlgEhEREQsfN2cyn1OXmYaceFvkvb7dgCcm12P1/BHsXUxX1ZfKhqcVUSV5hDNmTMHk8lU6Mff39/yvmEYzJkzh8DAQJydnenTpw+HDh0qdI3MzEwefvhhvL29cXV1ZdSoUZw+fbpQm4SEBCZMmIDZbMZsNjNhwgQSExOvxkcUERGpUboGeeLvbn1QlBn9J1Fhj+QHQza21O8zCZ9bnr3sYAgqFpxVVJUnVV933XVERUVZfg4cOGB5b+HChbzyyissWrSIPXv24O/vz8CBA0lJSbG0CQ0N5euvv2bFihXs2LGDCxcuMGLECHJzcy1txo0bR2RkJOHh4YSHhxMZGcmECROu6ucUERGpCWxtTMwZFVxmO8MwSN67juhlM8hJjMLW3Qf/cS9g7naT1flCJTGRv9qsa5DnZV2nXPesymX3c+bMYfXq1URGRhZ5zzAMAgMDCQ0NZebMmUD+aJCfnx8LFixg8uTJJCUl4ePjw9KlS7n11lsBOHv2LI0aNeKbb75h8ODB/PbbbwQHBxMREUG3bt0AiIiIICQkhN9//51WrVpZ1VctuxcRkbqktDpEeRkXiNv4Bmn/3QmAc4tueA0LxdbZzdLG0c6GiT2a8N4Px4D8vCBrFKRQvzO+U6Usvbf2+7vKc4j++OMPAgMDcXR0pFu3bsybN49mzZpx7NgxoqOjGTRokKWto6MjvXv3ZufOnUyePJm9e/eSnZ1dqE1gYCBt2rRh586dDB48mF27dmE2my3BEED37t0xm83s3LnT6oBIRESkLhnSJoCBwf5FKlUf/GUfn857jLRzZ7C1s8Pc+27qdR6FyZQfypiA4W39ef32TtjamOjY2IO56w4TlfR3PpC/uyO3d21MUno2qyPPEp+a9fd7lVyHyFpVGhB169aNTz75hGuuuYZz587x/PPP06NHDw4dOkR0dDQAfn5+hc7x8/PjxIkTAERHR+Pg4ICHh0eRNgXnR0dH4+vrW+Tevr6+ljbFyczMJDMz0/I6OTm5Yh9SRESkhrK1MdGzhTc9W3hjGAavvfYai2fOJDs7m6CgIFauXEn7jp1Zuus4J+LTaOLpwoSQpjjY/T1lVhBY7T4WT0xKBr5uToWqUT89PLjE966mKg2Ihg4davlz27ZtCQkJoXnz5nz88cd0794dwBJxFjAMo8ixS13aprj2ZV1n/vz5zJ0716rPISIiUlvl5hl8u/8vnpk+lT0/bAbg5ptvZvF77/PfhDw2HowiONDMxJ5BJQYytjamEpfPl/be1VTlU2YXc3V1pW3btvzxxx+MGTMGyB/hCQj4e9gsJibGMmrk7+9PVlYWCQkJhUaJYmJi6NGjh6XNuXPnitzr/PnzRUafLvbkk08yffp0y+vk5GQaNWp0WZ9PRETkasnNM6waebm0XecmHuw9kUBMSgbHY1N574uNHFn+b3JTzoOtHU2HT6HT3ZMZsXhfoWmwyt5y42qrVgFRZmYmv/32G7169SIoKAh/f382b95Mx44dAcjKymLbtm0sWLAAgM6dO2Nvb8/mzZsZO3YsAFFRURw8eJCFCxcCEBISQlJSErt376Zr164A/PTTTyQlJVmCpuI4Ojri6Oh4JT+uiIjIFRF+MKpI3k5xAUtx7WxMkGeAYeSRvHsVids+ASMPO48AfEbPwvBrznvbjxe5Z8GWG5WVDH21VekqsxkzZjBy5EgaN25MTEwMzz//PNu2bePAgQM0adKEBQsWMH/+fJYsWULLli2ZN28e33//PUeOHMHNLT+TfcqUKaxfv56wsDA8PT2ZMWMGcXFx7N27F1tbWyB/au7s2bMsXrwYgPvvv58mTZqwbt06q/uqVWYiIlIThB+MYsqyfUVWdV26equkdgC5aUnEbniFjKN7AXBpfSNegx/CxtGl1HubyE+K3jGzX5XkARWnRqwyO336NLfffjuxsbH4+PjQvXt3IiIiaNKkCQBPPPEE6enpTJ06lYSEBLp168amTZsswRDAq6++ip2dHWPHjiU9PZ3+/fsTFhZmCYYAPv30U6ZNm2ZZjTZq1CgWLVp0dT+siIjIFZabZzB33eFig5yL9wjrd61fie0yTh0kdu1Cci/EY7JzwKP//dRrP7jM/N2Ce1ztLTcqS5WOENUkGiESEZHqbtdfcdz+fkSZ7Z4Z3pp/bfit0DHDyCN51xck7vg0f4rMsyE+o2fi4BtU7n68emsHRrUPrBarx2rECJGIiIhUHmv3/vrhj9hCr3NTE4hd/woZx/cD4HpdXzwHTcXGwblC/Xhm9QHmrjtUqKhjdU+6rvKtO0RERKRyWLv317b/nrf8Of3EL0QtmUbG8f2Y7BzxGvoIXsOnVzgYAriQmVukwnVB0nX4wagKX/dK0giRiIhILdE1yJMAs1OhVWPFMQEYuST8uIKkH1cABvZejfEeMwsH78ZXpG8X5zANDPavNknXBTRCJCIiUkvY2piYPbLsjVmzL8QTteKfJP24HDBwbTsQ/7teuWLBUIGLk66rG40QiYiI1CJD2gRwT8+mfPjj8WLfTz+2j9j1L5OXloSDkzONRj5CTrMbrmofrc11upo0QiQiIlKL5OYZBNYvmv9j5OWS8MMnxHw+m7y0JOx9mhK2ZgtHVszjmeGtr2ofrc11upo0QiQiIlJLFFd5GiAnOZbYdS+SefoQAPU6DCH4/x5m7IBu2NqYmNgziA92HCsz9+hiJsD0v6rW5TnH35y/BL+60QiRiIhIDZabZ7DrrzieW3eIB5btKxLUpP+1h6iwaWSePoTJwRnvUU/gNfghxlzf1JLYXJB7ZG2ac0G7+3oF5QdG5Thn9sjgapdQDRohEhERqbFKGhECMHJzSPzhE5J3rwLAwa853qNnYu8RCMDaX6J4YkhrS3AypE0A74zvVOL1LuZ/UU2hjo09ipxT38UeoNDSe/9qXodIlaqtpErVIiJSnZS2F1lOcgyxaxaSefZ3ANw6jcCj7z2Y7OwLtVt+X/ciW2zk5hnsPhZPdFI68alZeNZzxLeeI5gg9kJmsVWnC865uCo1oErVIiIicuWUtmdZ2h8/EffNq+RlXMDk6IrX0Gm4tupZ7HWKW+1la2Mq9z5kJZ1Tk/YzU0AkIiJSRYobWbFmFCXiaFyRaS0jN5uE78NI+XkNAA4BLfEeNRP7+v4lXqc6rvaqKgqIREREqkBx+T/W7PcVfjCKWV8dKHQsOzGa2LULyIr6AwC360fj0WciJlv74i5RrVd7VRUFRCIiIldZSfk/Bft9vTO+U7FBUXHnpR75kbiNb2BkpmLjVA+vYY/i0rJbifeu7qu9qooCIhERkauotPyf0vb7uvQ8IyeLhK0fkrJvAwCOgdfiPeoJ7My+pd6/uq/2qioKiERERK6i3cfiS13WXrDfV9iPx/B2c7TkFl18XnbCWWLXLCDr3F8AuHe7mfq9JmCyLflr/aG+zenZwqfKVntVdwqIREREriJr9/H614bfLH92dbSlhY8rAKmHtxH3n0UYWenYOLvjPfxRnJt3KfN6Lf3catSqr6tNAZGIiMhVVJGVXamZuew/dp6E797nwi/hADg2vA7vkY9j5+59xe5blyggEhERuYq6BnkSYHYiOimj2Dyi4mTHneb8mhfIPn8cMOEeMpb6N4zDZGNb5rlaUWYd7WUmIiJyFRXsGwbW7QF24dBWoj4OJfv8cWxc6uM79jk8bpxgdTAEWlFmDQVEIiIiV1nBvmH+5pKnsfKyM4j95jXi1r+MkZ2BY+N2BNz9Bs5BHUs8x9O1cN0hf7NTiUv4pTBNmYmIiFxBJVWjHtImgIHB/pb3YlMyLYnUWedPELtmAdlxJwET5p63Y+5xa5mjQs+MuA5/d6cq3z+sJlJAJCIiUskKgqBvD0fzdeQZ4lP/3vX94mrUF+8Blptn8P72o/y5Yz3xm9/FyMnE1tUD75GP49SknVX39Xd30kqyClJAJCIiUomK25LjYiVVo05PS8V117vEbfwKAKemHfEeMR1bVw+r7hugxOnLooBIRESkkpS0JcfFiqtG/euvvzJ27FiOHDmCjY0Nvn3vxKHzTZhM1qf6PjNcidOXQ0nVIiIiZcjNM9j1VxxrIs+w6684cvOKhjylbclxqYJq1Lv+jGXmv1/h+i5dOXLkCA0aNOD777/n9KaPmD7wWlwcyl5JVsDD1cH6DyRFaIRIRESkFGXtSl+QL/Tjn+dL3ZLjUnmZaYy8eSyJB78HwKlZZ/xvf4pUjxbY2ph4ZEBLHurXglc2H+GtrX+VeT1rK2BL8RQQiYiIlKCsXenvvzGItb9ElSsQAsg69xfn17xATkIUmGyo3/su3Lv+H3E5NoXyi2xtTNzQwseqgEiVqC+PAiIREZFilLUrPcDiH46V65qGYXBh/wbit3wAuTnYuvvgM+oJHBu0tlz30vyisipbqxJ15VAOkYiISDHK2pW+vPIyLhC7ej7xm9+F3BycW3QjYOIblmCoQEF+0e5j8UDpla1VibryKCASEREpRmXm5GRG/ZeosEdI++9OTDZ2ePS7D5+b/omts5tV9y+psrUqUVceTZmJiIgUozJycgzDIOXntSR8vwTycvD0b8iCtz7i+d1Z5b7/pZWtVYm6cikgEhERKUZFdqW/WG56CnHfvEb6nz8B4HJND1yHTiOwRRsCjhyuUE7QxZWtpXJpykxERKQY5d2V/mKZZ34jKmxafjBka4fnwAfwHvMktk71+NeGwzwzvHWx11VOUNVRQCQiIlICa3alv5hh5JH001dEfzaL3OTz2NUPwH/8S7h1GoHJZLIkTHu4OionqJrRlJmIiEgpCnJ3wn48ZtmNvji5aUnEbXiV9KM/A+BybS+8hjyMjaNLkbYxKRmM7tBAOUHViAIiERGptQqqSF9uwGFrY8LbzbHE9zNOHSR27YvkXogDW3s8B0ymXvvBmEzF36sgYVo5QdWHAiIREamVytpyo7y86xUNiAwjj+SIL0ncvgyMPOw8G9L81qfJdG+kIoo1jHKIRESk1inYcuPSwooFW26EH4wq/0UviXByUxOJ+Xw2iT98AkYertf1JeCuV5n8f/0AJUzXNAqIRESkRrt0J/qsnLwyt9yYu+5wsTvWlyY2NdPy54wTvxK15GEyju/HZOeI19BH8Bo+HRsHZ5r5uCphugbSlJmIiNRYxU2Lebo6EJ9acuHDi7fGKE/+jq+bE0ZeLkk7V5K0cwUYedh7NcZ79EwcfJoUahfS3EsJ0zWMAiIREamRStqJvrRg6GLl3ZqjkVMmiV/NIfnofgBc2w7Ac8AD2DgUHglK+N/9lTBds2jKTEREapzSdqK31sVbY1w67XbpdNrmzZvp3KkjyUf3Y7J3wmv4dLyHhRYJhgD+taH803FS9TRCJCIiNc7l7kRvY4KE/+UElbYabcC1PsyZM4d58+ZhGAbNWwWT3uth7L0alXjtikzHSdVTQCQiIjXO5e5En2fAg5/t5/7Tibz3w7EiI03RSRnc9/Ym6v/0Dgf35e9FNnnyZPpNnMETq49c8f7J1aeASEREahxrd6I3UWS1vIUBvL+9aDAEkPbXz8RueIXT6cm4ubnx3nvvcdttt7Hrr7hK7Z9UHwqIRESkxilrJ3oT+avN4spIsL401cfIzSFx+1KSf/oKAAe/5nzw6WeM7d/V6vuq8GLNpKRqERGpcUrbib7g9egOgeW6Zk5yDOc+m2UJhtw6Dcd//Is4ejUo131VeLFmUkAkIiI1Ukk70RcUQBwY7G/1tdL++ImoJdPIPPs7JkdXvMc8iefAKZjsHIpMf5V1XxVerJk0ZSYiIjVWwU70xRVAzM0zSp3eAjDlZRO/NYzkn9cA4BDQEu9RM7Gv71/q9Fdp95WaSQGRiIjUSJfuZD+iXf4U2cXHnh56LQ+tiCxyrgnITozGtPV1kv97AAC360fj0WciJlt7q6a/VHixdlFAJCIiNU5xtYPqu9gDkJiWbTlW0oCNw6k9nF/3KqkpydRzr0/g6OlkBnayvO//vzpEmv6qOxQQiYhIjVLSlh0XB0IFiqwiy8kmYeuHpOxbD0D37t1ZsWIFDRs11vRXHVdtkqrnz5+PyWQiNDTUcswwDObMmUNgYCDOzs706dOHQ4cOFTovMzOThx9+GG9vb1xdXRk1ahSnT58u1CYhIYEJEyZgNpsxm81MmDCBxMTEq/CpRESkMl3Olh3ZCWeJXjbDEgwNvv1+pr++nLM59QAIae7F6A4NCGnupWCoDqoWAdGePXt47733aNeuXaHjCxcu5JVXXmHRokXs2bMHf39/Bg4cSEpKiqVNaGgoX3/9NStWrGDHjh1cuHCBESNGkJuba2kzbtw4IiMjCQ8PJzw8nMjISCZMmHDVPp+IiFSOim7ZkfrbdqLCHiHr3F/YOLvje8tsfm88isdXHeL29yO4YcEWwg9GAWXvaya1k8kwjCr9m75w4QKdOnXi7bff5vnnn6dDhw689tprGIZBYGAgoaGhzJw5E8gfDfLz82PBggVMnjyZpKQkfHx8WLp0KbfeeisAZ8+epVGjRnzzzTcMHjyY3377jeDgYCIiIujWrRsAERERhISE8Pvvv9OqVSur+pmcnIzZbCYpKQl3d/cr8zBERKRUayLP8EgxSdIlycvOJGHL+1yIDAfAsWEw3iOfwM7du0hbE3D/jUGs/SWq2H3NlE9UM1n7/V3lI0QPPvggw4cPZ8CAAYWOHzt2jOjoaAYNGmQ55ujoSO/evdm5cycAe/fuJTs7u1CbwMBA2rRpY2mza9cuzGazJRiC/Dljs9lsaSMiIjVDebbEyI47TfTSx/4XDJlwDxmL3+3ziw2GIH8rj8U/HCsyAhWdlMGUZfssI0hSO1VpUvWKFSvYt28fe/bsKfJedHQ0AH5+foWO+/n5ceLECUsbBwcHPDw8irQpOD86OhpfX98i1/f19bW0KU5mZiaZmZmW18nJyVZ+KhERuVLK2jqjwIVDW4n/z1sY2RnYuJjxHvEYzkGdSjmjZAb5o0dz1x1mYLC/8otqqSobITp16hSPPPIIy5Ytw8mp5IjfZCr8i2cYRpFjl7q0TXHty7rO/PnzLUnYZrOZRo0alXpPERG58krbOgMgLzuD2G9eJ279yxjZGTg2bkvA3W9WOBgqYABRSRnsPhZ/WdeR6qvKAqK9e/cSExND586dsbOzw87Ojm3btvHGG29gZ2dnGRm6dBQnJibG8p6/vz9ZWVkkJCSU2ubcuXNF7n/+/Pkio08Xe/LJJ0lKSrL8nDp16rI+r4iIVI6Sts7Iij1J9CfTST2wGTBh7nk7frc+j129yttoNSal/AndUjNU2ZRZ//79OXDgQKFjd999N9deey0zZ86kWbNm+Pv7s3nzZjp27AhAVlYW27ZtY8GCBQB07twZe3t7Nm/ezNixYwGIiori4MGDLFy4EICQkBCSkpLYvXs3Xbvm71b8008/kZSURI8ePUrsn6OjI46OjpX+uUVExDq5eQYRR+PY9VccYBDSzJvu/1sSP6RNAP2u9ePjncfZfSyO9V8u51z42xg5mdi6euA1cgbOTdpXep+syWG6tIK2ahrVDFUWELm5udGmTZtCx1xdXfHy8rIcDw0NZd68ebRs2ZKWLVsyb948XFxcGDduHABms5l77rmHxx57DC8vLzw9PZkxYwZt27a1JGm3bt2aIUOGcN9997F48WIA7r//fkaMGGH1CjMREbm6wg9GMWvVgULFFhdt/Yv6Lva8cFNbAGatOkB8YjLxm98h9eAWAJyadsR7xHRsXT2KvW5Flbav2aX9vrSCtlap1QzVulL1E088QXp6OlOnTiUhIYFu3bqxadMm3NzcLG1effVV7OzsGDt2LOnp6fTv35+wsDBsbW0tbT799FOmTZtmWY02atQoFi1adNU/j4iIlC38YBQPLNtX7HuJadmW97LOH+f86hfIiT8NJhvq33AH7iH/wGSq3GwQa/Y1K+h3cRW0C1apvTO+k4KiaqzK6xDVFKpDJCJScdZOI+XmGfR8YQvRySXn6hiGwYVf/kPCd+9h5GRhW88T71FP4NSoTYnnlIenqwPxqVmW19aM8OTmGdywYEuJRSMLRph2zOyn6bOrzNrv72o9QiQiIjVfeaaRdh+LLzUYystMI+4/i0j77QcAnJp1xnv4dGxdzJfdz4KgZdvjfdl7IqFcOUBlVdC+eJVaSHOvy+6rVD4FRCIicsWUdxqptFVcWef+4vyaF8hJiMqfIut9J+5db6qUKbKLp8Uc7GzKHbRYu/pMq9SqLwVEIiJyRZS2EWtJxQ6LW8VlGAYX9n9D/Jb3ITcHWzef/Cmyhq0rra/+l5n4bG0F7fJU2parSwGRiIhcERWZRuoa5Im/u5Nl2iwvM5W4jW+QduRHAJxbdMVr2KPYOruVdFmrPT2sNb7ujpWyNL6sCtrWrlKTqlPle5mJiEjtVJFpJFsbE3NG5Veizoz6L1FLpuUHQzZ2ePS7F5+bnqmUYAjA192R0R0aEPK/2kaXo7QK2tauUpOqpYBIRESuiIpOIw2+zp8h7OPcp0+Qk3QOW7Mf/ncswL3LGDxcHXh3fCceHdDyqvXPWiVV0PY3O2nJfQ2gKTMREbkiKjKNFB8fz6RJk1izZg0AHW8cTNvbZuLhUZ+bOzakR0tvbG1MDAz2Z/nuU6WuSCvJlZy+GtImgIHB/qpUXQNphEhERK6I8k4jRURE0LFjR9asWYOdvQNNRz5EXPeH2HYijdWRZ3li1a9sPhxtufacUcGYSrm2tfetbLY2JkKae1XadJxcHQqIRETkirFmGikvL4+XXnqJXr16cfLkSQIaNcVn3EKM4CGYTH8HEwVL9cMPRpV57XfHd+JdTV9JOahStZVUqVpEpOJKqlQdGxvLxIkT2bBhAwBjx97K0WvHcT7TttjrFFfxubQq2NpoVaz9/lZAZCUFRCIil+/iACXqSCT/nvEAp0+fxtHRkddff522/W9m3Ac/lXmd5fd1V8VnsYq27hARkWqlYAuPs4lpJEd8SeL2ZWDk0aBpMzasXkX79u1ZE3nGqmup4rNUNgVEIiJyxRVs4ZGTmkjs+pfJOL4fANfr+mI7aCpRtr60x/ql8MdjU69gb6UuUlK1iIhcUQVbeKSf/JWosGlkHN+Pyc4Rr6HT8Bo+HRsHZ+auO0xunmGpVF2W5btPkpunjA+pPAqIRETkitr153l++2YJ51b8k9wL8dh7Ncb/zleo124QJpOp0BYetjYmbu/auMxrRidnsvtY/JXvvNQZmjITEZFCLl2Z1bmJB3uOxbPraCyQX2One7Oi9XWKW9F1PuYcd/5jDEm/5idKu7YdgOeAB7BxKDoKVJAX1NTbxap+Ko9IKpMCIhERsShIfL54U1aTCS5ej7xo65/Ud7HnhZvaWur5FHeey/lDnFy1kNTEOEz2TngOmkq9Nv1KvHdB/pB2jpeqoIBIRESAvxOfL83MKa44S2JaNg8s28e74zsBFDrPyMslacdnnNj1OWBg79MUn9EzsfdqVOx9L91KQzvHS1VQQCQiIpbE5/KmKc9ddxjDMCzn5aTEErvuJTJPHQSgXvshePS/Dxt7xxKvYVB4K42CLT+mLNuH6X/vFyhu6w0VX5TKoIBIRETYfSy+0HSXtS4+J/3oXmLXv0xeejImB2e8Bj+Ea3DvMq8xqWfTIltpFGzLcek0nL/Zidkjg0udqgu4pI2INRQQiYjIZSUoG7k5JG5fRvJPXwJg79ssf4rMs4FV5w8M9i/2eFk7x5c0xVew55n2LJPyUEAkIiIVTlDOST5P7NqFZJ75DQC3TsPx6HsPJjuHMs+1JheoYOf4S5U2xWf879pz1x1mYLC/ps/EKgqIRESkzETm4qT9+RMJ37xGTnoKJgcXvIZOw/XaG6w6t7hcoPIoa4rv4tpG2vNMrKHCjCIiYklkhr+DlZIYudnEb/mA81/9i5z0FFpe157Au9+wOhiC/JGhy5nSsnaKT7WKxFoKiEREBPg7kdnfXHj6zHRRhJSdGE30pzNJ2bMagNDQUA7u280HDw3D09Xeqvs8M7w1O2b2u6z8HtUqksqmKTMREbEoLpG5oFL1e0uX8/myp8lKTaF+/fqEhYUxevRoy3n9rvWj+/zviE/NKvbaBTlDE3sGXXZej2oVSWVTQCQiUkfl5hlEHI3jxz9jOZuYToP6zvRo4U33Zl6WvJvcPIMffjvLM0/N4se1ywDo1r07K1esoEmTJoWu52Bnw7z/a8OUZfuAsusHXY7y1ioSKYvJMIqrQSqXSk5Oxmw2k5SUhLu7e1V3R0TksoQfjGLWqgMkpmUXec/V0ZYXb26HjY2JRz/4D0dX/pus6D8BcO96Ew0G3M3kvq14qF+LYgOOq1kbSHWIpCzWfn8rILKSAiIRqS3CD0bxwP9GcUqT+tt24sLfwMhKx8bZHa/hj+LSvIvl/Uv3M7vY1awerUrVUhoFRJVMAZGIVHfWBAa5eQY9X9hCdHIpS9Zzsoj/7n0uRG4EwLFhMN4jn8DO3bvY9u+qAKJUY9Z+fyuHSESkFrB26mj3sfhSg6Hs+DOcX/MC2THHAHDv/g/q9xqPyca2xHNUAFFqAwVEIiI10MWjQcdjU3n12z+KtCluC4vS6vJcOLSV+P+8hZGdgY2LGe/h03Fu1rnMvqgAotQGCohERGqY4kaDilPcFhbF1eXJy84g4dv3uPDrJgAcG7fFe8QM7NysD3BUAFFqOgVEIiI1SEkbmpbk0i0sugZ54u/uZJk2y4o9SeyaF8iOPQmYMPe4DXPP20qdIiuOCiBKTadK1SIiNURpG5qWpWAEx9bGxJxR+Vt0XDjwLdGfPEp27ElsXT3wve156ve6o1zBkIn8XCUVQJSaTiNEIiI1RFkbmpbG183JkneUmJyC774POLF5NQBOTTrgPfIxbF09ynVNFUCU2kQBkYhIDVHRPB0PF3sSUrO4YcEWTvz5O+dXv0BO/Gkw2TDsrmkE9rmN7X/Gk5GdV67r+qsAotQiCohERGqI47FpFTovIS2bKZ/u5cKvm0j4djFGTha29TzxGfUEh/3a8HDnxrw7oSu7j8Xz45+xLNr6Z5nXfGZ460rZk0ykulBAJCJSA+TmGSzffbJC5+ZlphH3n7dI+20bAE5BnfEeMR1bFzPw9yq0gqTrr/adLnPTVAVDUtsoqVpEpAYoq6BiSbLOHSXq49D8YMhkQ/3eE/H9x2xLMHTxKjT4e9NU+DtHqIByhqQ2U0AkIlIDlDd/yDAMUvZ/Q9TSx8hJOIutmw9+4xZg7n4LJlPRf/ovvv6QNgG8M74T/ubCS+n9zU6FijyK1CaaMhMRqQHKU+cnLzOVuI1vknZkBwDOLbriNSwUW+eS93G69PpD2gQwMNhfm6ZKnaGASESkBuga5El9F3sS07JLbZcZ9QexaxeQkxgNNrZ49J6IW5cxmEzFBzIFOUHF1RGytTFpOw6pMxQQiYjUAoZhkLJ3LQlbl0BeDrZmP3xGPYFjYKvSz0M5QSKggEhEpEbYfSy+xNGh3IwLxH3zGul/RADgfE0I3kMfwcapXpnXre9iX6n9FKmplFQtIlIDlJRUnXn2CFFLpuUHQ7Z2eAyYjM+Yp6wKhgCS0rKZsmwf4QejKrO7IjWORohERGqAS5OeDSOPlD2rSdj2MeTlYlc/AO/RM3H0b1Gu6xrk5xEV1CLS1JnUVQqIREQqUcF+YZW9MqtzEw88Xe2JT80mNz2ZuA2vkv7XHgBcru2F15CHsHF0rdC1L65FpCRqqasUEImIVJLwg1HMXXe40AasAZWw31fBdeNTs8k4fZjYtQvJTYkFW3s8+99HvQ5DS1xFVh4V3StNpDZQQCQiUgnCD0YxZdm+IttdRCdlMGXZvhILGl48ouRdzxEMiE3NtIwu/edgNFM/24dh5JH801ck/rAUjDzsPBvgM3omDr7NKu0zlKfWkUhto4BIROQy5eYZzF13uNi9v0rL0SluROliZmc7kjNyyE1NJHbDK2Qc2weAa3AfPAdNxcbRpVL6X1otIpG6QgGRiMhl2n0svsSgBorP0SlpROliSek5ZJw8QOy6F8m9EI/JzhGPAZOp125gpUyRgfYnEymggEhE5DJZm3tT0K60EaUCRl4uSbs+J+nH5WDkYe/VCO/RM3HwaXr5Hb6IfyXkOInUBlYHRL/++qvVF23Xrp1V7d555x3eeecdjh8/DsB1113Hs88+y9ChQ4H8yqtz587lvffeIyEhgW7duvHWW29x3XXXWa6RmZnJjBkzWL58Oenp6fTv35+3336bhg0bWtokJCQwbdo01q5dC8CoUaN48803qV+/vtWfSUSkJNbm3hS0K2tEKfdCArHrXyLjxC8AuLYZgOfAB7BxqHiOz9PDWuPr7oi3qyOYIPZCpvYnE7mI1QFRhw4dMJlMGIZR5lBtbm6uVdds2LAhL7zwAi1a5NfN+Pjjjxk9ejT79+/nuuuuY+HChbzyyiuEhYVxzTXX8PzzzzNw4ECOHDmCm5sbAKGhoaxbt44VK1bg5eXFY489xogRI9i7dy+2trYAjBs3jtOnTxMeHg7A/fffz4QJE1i3bp21H19EpERdgzwJMDsRnZRR4qiPp6s9nZt4kJtnsOOP8yVeK/14JLHrXyIvNRGTvSOeg6ZSr03/CvetID9o0g1BCnxESmEyDKO0UVuLEydOWP68f/9+ZsyYweOPP05ISAgAu3bt4uWXX2bhwoWMGTOmwh3y9PTkxRdfZNKkSQQGBhIaGsrMmTOB/NEgPz8/FixYwOTJk0lKSsLHx4elS5dy6623AnD27FkaNWrEN998w+DBg/ntt98IDg4mIiKCbt26ARAREUFISAi///47rVqVvs9PgeTkZMxmM0lJSbi7l7xjtIjUTQU5QUCJQVF9F3uycvJIyyr6P41GXi5JPy4naedKwMDeuwk+Y2Zh79Wown0qCH9KWuEmUhdY+/1t9QhRkyZNLH/+xz/+wRtvvMGwYcMsx9q1a0ejRo145plnKhQQ5ebm8sUXX5CamkpISAjHjh0jOjqaQYMGWdo4OjrSu3dvdu7cyeTJk9m7dy/Z2dmF2gQGBtKmTRt27tzJ4MGD2bVrF2az2RIMAXTv3h2z2czOnTutDohEREozpE0A74zvVOqqsZL2IstJiSN23YtknjoIQL32g/Hofz829o4l3s/GBPfcEMT6X6NKvJ/yg0SsV6Gk6gMHDhAUFFTkeFBQEIcPHy73tUJCQsjIyKBevXp8/fXXBAcHs3PnTgD8/PwKtffz87OMVkVHR+Pg4ICHh0eRNtHR0ZY2vr6+Re7r6+traVOczMxMMjMzLa+Tk5PL9blEpO4ZGOyPq4MdD3y6l9RM61IH0o/uJXb9y+SlJ2NycMZr8EO4Bvcu87xFt3dkWLtAZg1t/XcdI+UHiVRYhQKi1q1b8/zzz/Phhx/i5JSf5JeZmcnzzz9P69aty3WtVq1aERkZSWJiIl999RV33XUX27Zts7x/ab6SNTlMl7Yprn1Z15k/fz5z58619mOISB1XVk2hSxl5uSRuX0pyxJcA2Ps2w2f0TOw9G5R6nperA6M7BOLh6khunoGtjUnbbYhUggoFRO+++y4jR46kUaNGtG/fHoBffvkFk8nE+vXry3UtBwcHS1L19ddfz549e3j99dcteUPR0dEEBPw93BsTE2MZNfL39ycrK4uEhIRCo0QxMTH06NHD0ubcuXNF7nv+/Pkio08Xe/LJJ5k+fbrldXJyMo0aVXwuX0RqL2tqCl0sJ/k8sWtfJPNM/oh6vY7D8ex3DyY7h2LbT+3TjMzsPL6OPENcahYf/Xicj348XinbgohIPpuKnNS1a1eOHTvGv//9b9q1a0fbtm2ZN28ex44do2vXrpfVIcMwyMzMJCgoCH9/fzZv3mx5Lysri23btlmCnc6dO2Nvb1+oTVRUFAcPHrS0CQkJISkpid27d1va/PTTTyQlJVnaFMfR0RF3d/dCPyIil7KmptDF0v7cTdSSaWSeOYzJwQXv0bPwGjSlxGAI4JNdJ/nwx+PEpxbOQSrYFiT8YNRlfAIRgcsozOji4sL9999/WTd/6qmnGDp0KI0aNSIlJYUVK1bw/fffEx4ejslkIjQ0lHnz5tGyZUtatmzJvHnzcHFxYdy4cQCYzWbuueceHnvsMby8vPD09GTGjBm0bduWAQMGAPnTe0OGDOG+++5j8eLFQP6y+xEjRiihWkQuW1k1hQoYudkkbvuE5D1fA+Dg3xLv0TOxr+9f5rkXMnOKvyYlbwsiIuVT4YBo6dKlLF68mKNHj7Jr1y6aNGnCq6++SrNmzRg9erRV1zh37hwTJkwgKioKs9lMu3btCA8PZ+DAgQA88cQTpKenM3XqVEthxk2bNllqEAG8+uqr2NnZMXbsWEthxrCwMEsNIoBPP/2UadOmWVajjRo1ikWLFlX0o4uIWFhTpTon6Rzn1ywkK+oIAG6dR+HR525MdvaXff/itgURkfKzug7Rxd555x2effZZQkNDef755zl06BDNmjUjLCyMjz/+mK1bt16JvlYp1SESkeLs+iuO29+PKPH9tP/uIu6b18jLTMXG0RWvYaG4XBNSpJ2TnQ0ZOXkV7sfrt3VgdIfSE7JF6iJrv78rlEP05ptv8v777/P0009jZ/f3INP111/PgQMHKnJJEZEaqXMTDzxdi+b/GDnZxH+7mPNf/5u8zFQcAloRcPebxQZDAI72Ffrn2MLa7UNEpHgVmjI7duwYHTt2LHLc0dGR1NTUy+6UiEhlyM0zLDV6rkRdnoKl9vGpWYWOZydEEbt2AVnRfwLg3vUm6t94Jybbkv/JTUovPk+oLAVbc3QN8qzQ+SKSr0IBUVBQEJGRkYWqVwNs3LiR4ODgSumYiMjlKK4uUGUuUy9pqX3q7zuI2/gGRlYaNk5ueI2YjkvzLlZds76zPUnp2VavWCsI7WaPDFZCtchlqlBA9Pjjj/Pggw+SkZGBYRjs3r2b5cuXM3/+fD744IPK7qOISLmUFKwULFO/3L29iltqb+RkEb/lAy7s/wYAxwbBeI96HDt3H6uvO7FHU17/7g9MFN4PreB1fRf7Qtt/aGsOkcpToYDo7rvvJicnhyeeeIK0tDTGjRtHgwYNeP3117ntttsqu48iIlYrrS5QZS1Tv3SpfXb8Gc6veYHsmGMAuHf/B/V7jcdkY1vSJYrVJciz2P3QCgKfgcH+V3QKUKQuq/Cy+/vuu4/77ruP2NhY8vLyit0vTETkaiurLlBlLFO/eKl96uHvifvPWxhZ6di4mPEePh3nZp0rdN3YC5mM7tCg1MBHS+tFrowKBUT9+vVj1apV1K9fH29vb8vx5ORkxowZw5YtWyqtgyIi5WFNXaDytCuOr5sTedkZJHz7Hhd+3QSAY+O2eI+YgZ1bxQOWgpVi2p9M5OqrUED0/fffk5WVVeR4RkYG27dvv+xOiYhUlLXLzy9nmbpbxjliP51B+rnjgAlzj9sw97yt3FNkBbRSTKTqlSsg+vXXXy1/Pnz4MNHR0ZbXubm5hIeH06CBCoOJSNXpGuRJgNmJ6KSMYvOILjf4+Pjjj5k6dSrpaWnYuNbHZ8QMnJp2qHB/tVJMpHooV0DUoUMHTCYTJpOJfv36FXnf2dmZN998s9I6JyJSXrY2JmaPDGbKsn3FrtaCv4MPa+sU5eYZbD14klmPhbL329UA9Ovfn3v++TJv7ootlLNkY4K8EtbNuzrYYm9no5ViItVQuQKiY8eOYRgGzZo1Y/fu3fj4/L2c1MHBAV9f30J7iImIVIUhbQJKXa01pE2A1XWKwg9G8cg76zi28nmy406ByQbzDeOIvWEcnt6+7JjZrlBQ1bmJB3tPJBCdnEFsSiYJaVnYmCCkmTfd/5cXpJViItVPhfYyq4u0l5lIzVPSCFBJdYoKwpKCOkUbD5xl/BMvkPDtYoycLGzreeI98nGcGre1nPPuZdY0EpEry9rv7wolVc+fPx8/Pz8mTZpU6PhHH33E+fPnmTlzZkUuKyJSqYpbrWVtnaIuDVyYMGEC8b/kr5p1CuqM94jp2LqYC50zZ+2hy6ppJCLVQ4V2E1y8eDHXXnttkePXXXcd77777mV3SkTkSrGmTtGJ/x7muvYdiftlC5hsqN97Ir7/mF0kGAKITs5k97H4K9hjEbkaKjRCFB0dTUBA0SFiHx8foqKiLrtTIiIVVVaidGn1hwzD4ELkRuK/ex9ys7F188Z71BM4NSx9j8bLqWkkItVDhQKiRo0a8eOPPxIUFFTo+I8//khgYGCldExEpLysSZQuqf5QXmYqceGLSPs9v5Zaux79iL/+Xmydy84ZvJyaRiJSPVRoyuzee+8lNDSUJUuWcOLECU6cOMFHH33Eo48+yn333VfZfRQRKVNBovSl02EFG7qGH8wfve7cxINL030yo/8kKuyR/GDIxhbPfvewa8t/aOBX9pZE/u6OKqgoUgtUaIToiSeeID4+nqlTp1oqVjs5OTFz5kyefPLJSu2giEhZyrOh694TCZY6QYZhkLJvPQlbP4TcHGzdffEZPRPHwFb8cjqJOaOCeWDZvlLvPWfUdUqoFqkFKhQQmUwmFixYwDPPPMNvv/2Gs7MzLVu2xNHRsbL7JyJSJms3dA378Rg//HEegNyMC8RtfJ30/+4CwLlld7yGhWLrVA/Izwsa3aEB747vxKxVBwoVUwSo72LPCze11ZJ7kVqiwrvdA9SrV48uXbpUVl9ERCrE2qTmf234DYDMs0c4v3YhuUnnwNYOj7734NZpBCbT3yM9BXlBQ9oEMDDYn4ijcez6Kw4wLEUWNTIkUntYHRDddNNNhIWF4e7uzk033VRq21WrVl12x0RErGVtUrNhGKTsWU3CtjDIy8Wuvj/eo2fh6N+iUDtPV/tCeUG2NiZ6tvCmZwvvyuy2iFQjVgdEZrPZ8n9PZnPRWhwiIlUlITUTkwlKq7ufm55M3IZXSf9rDwAurW7Aa+jD2Di6Fmn7fx0aaPRHpI7R1h1W0tYdItVT+MGoMhOfM04fJnbti+SmnAdbezz730e9DkMLTZFdbPl93YtUuBaRmumKbt0hInI1lVRsMTfPYM7aQyWeZxh5JP/0FYk/LAUjDzuPQHxGz8LBr1mx7U3kbwCrZfQidY/VAVHHjh1L/L+pS+3bV/r/rYmIWKu0YotmZweikzOLPS83LYnY9a+QcWwvAC7BvfEa9CA2ji7Fti/41232yGBNl4nUQVYHRGPGjLH8OSMjg7fffpvg4GBCQkIAiIiI4NChQ0ydOrXSOykidVNJu9IXFFvse61PsedlnDpI7NqF5F6Ix2TngMeAB6jXbmCp/1Pnf0lFaxGpWyqUQ3TvvfcSEBDAv/71r0LHZ8+ezalTp/joo48qrYPVhXKIRK6u3DyDGxZsKbW+kI0JS5FFACMvl6SIL0ja8RkYedh7NcJ79EwcfJoWOddEfn2iST2bMjDYv8ieZyJSO1zRHKIvvviCn3/+ucjx8ePHc/3119fKgEhErq6yii1C4WAoNzWB2HUvkXHiFwBc2/THc+AUbByKX5KvESERuViFAiJnZ2d27NhBy5YtCx3fsWMHTk7a5FBELl95dpBPPx5J7PqXyEtNxGTviOfAqdRr279Iu5Ht/BkQ7F8oMVtEBCoYEIWGhjJlyhT27t1L9+7dgfwcoo8++ohnn322UjsoInWTNcUWjbxckn5cQdLOFYCBvXcTfEbPwt67UaF2Nia4r1cQTw4LvkK9FZGarkIB0axZs2jWrBmvv/46n332GQCtW7cmLCyMsWPHVmoHRaRu6hrkSX0X+yJ7iBXISYkjdv1LZJ48AEC9doPwGHA/NvaFA6knh17L3T2DcLCzueJ9FpGaq8J1iMaOHavgR0SqRPrRvcRueIW8tCRMDs54DX4Q1+A+xbZt17C+giERKVOFA6LExES+/PJLjh49yowZM/D09GTfvn34+fnRoEGDyuyjiNRBu4/FFxkdMvJySdy+lOSILwGw9w3KnyLzLPnfnPLkIolI3VWhgOjXX39lwIABmM1mjh8/zr333ounpydff/01J06c4JNPPqnsfopIDVdStemS3rs0kMlJPk/s2hfJPHMYgHodh+PZ7x5Mdg6l3tfajV9FpG6rUEA0ffp0Jk6cyMKFC3Fzc7McHzp0KOPGjau0zolI7VBatWmg2Pdu69LY8jrtrz3ErX+FvIwUTA4ueA2dhuu1N5R5X393R23DISJWqVBAtGfPHhYvXlzkeIMGDYiOjr7sTolI7VFatemSNmWNTsrg1W//i5GbQ+IPn5C8exUADv4t8B41E3sP62oHzRl1nZbWi4hVKhQQOTk5kZycXOT4kSNH8PEpvpS+iNQ9uXkGc9cdLhIMAcUeu/i9nKQYzq9ZQFbUEQDcOo/Eo88kTHb2Zd63vos9L9zUVkUXRcRqFQqIRo8ezXPPPcfnn38OgMlk4uTJk8yaNYubb765UjsoIjWXNdWmi5P2313EffMaeZmp2Di64jUsFJdrQso8r2dzL6b2aUH35l4aGRKRcqnQWtSXXnqJ8+fP4+vrS3p6Or1796ZFixa4ubnx73//u7L7KCI1VHlXeBm52cR/+x7nv/43eZmpOAS0IuDuN6wKhgCa+9ajZ0tvBUMiUm4VGiFyd3dnx44dbNmyhX379pGXl0enTp0YMGBAZfdPRGqo3DyD2JRMq9tnJ0YTu2YBWdF/AODe5f+o3/tOTLZlT5EVaOLpUu5+iohABQKinJwcnJyciIyMpF+/fvTr1+9K9EtEarDiVpWVJvX3HcRtfAMjKw0bJze8hj+KS4uu5bqnjQkmhDStQG9FRCoQENnZ2dGkSRNyc3OvRH9EpIYraVVZcYycLOK3fMiF/RsAcGwQjPeox7FzL7w4w8XBlrSs0v/Nua+XtucQkYqr0L8e//znP3nyySeJj4+v7P6ISA1W2qqyS2XHnyFq6QxLMOTe/R/43T6vSDAE8P6d1/PpPd3o2MjMpdlBNiaYfKM2bhWRy1OhHKI33niDP//8k8DAQJo0aYKrq2uh9/ftK762iIjUbtauKks9vI24/yzCyErHxtkd7xGP4dysc5F2JsDf7ET3Zvmrxnq2vIGsnDyW7jrOifg0mni6MCGkqUaGROSyVSggGjNmDCaTCcOw5v8DRaSuiE5KL/X9vOxMEr57jwu//AcAx0Zt8B75OHZuXkXaFowEzR4ZXGjVmIOdDff0alZpfRYRgXIGRGlpaTz++OOsXr2a7Oxs+vfvz5tvvom3t/eV6p+I1BDhB6P414bfSnw/O+4U51e/QHbsCUwmE4F97sD2+n9gsrEF8qe+8i76fyz//23toeKKInI1lCsgmj17NmFhYdxxxx04Ozvz2WefMWXKFL744osr1T8RqQHKSqS+cPA74je9jZGdiX09D9Z9tZIBAwYU2tC1cxMP9p5IKHbzVxGRK61cAdGqVav48MMPue222wC444476NmzJ7m5udja2l6RDopI9VZaInVeVgbxm98l9eC3ADg1ac+HYR8zuE97AEKaF54qu/S1iMjVUq6A6NSpU/Tq1cvyumvXrtjZ2XH27FkaNWpU6Z0TkeqvpETqrPMniF2zgOy4k2CywdzzdswhY/H09q2CXoqIlK5cAVFubi4ODg6FL2BnR05OTqV2SkRqjku35zAMgwu/bibh28UYOZnY1vPEe+TjODVuiwmYu+4wA4P9NR0mItVKuQIiwzCYOHEijo6OlmMZGRk88MADhZber1q1qvJ6KCLVUm6ewe5j8fxxLsVyLC8zjfhNb5N6+HsAnII64T18Orau9YH8XeyjkjLYfSxe02MiUq2UKyC66667ihwbP358pXVGRKqvggAoJiWDY+dT+XjXcRLSsi3vZ8Uc5fyaBeTEnwGTDfVvnIB7t5sxmYrWCCrvpq8iIldauQKiJUuWXKl+iEg1VtreZIZhcOGXcOK/fQ9ys7F188Z71OM4NbyuxOv5ujldye6KiJRblZZ3nT9/Pl26dMHNzQ1fX1/GjBnDkSNHCrUxDIM5c+YQGBiIs7Mzffr04dChQ4XaZGZm8vDDD+Pt7Y2rqyujRo3i9OnThdokJCQwYcIEzGYzZrOZCRMmkJiYeKU/okiNV7CkvrhgKC8zjdi1C4n/z1uQm41z8y4E3P1GqcFQgDl/Sb2ISHVSpQHRtm3bePDBB4mIiGDz5s3k5OQwaNAgUlNTLW0WLlzIK6+8wqJFi9izZw/+/v4MHDiQlJS/8xZCQ0P5+uuvWbFiBTt27ODChQuMGDGi0Aa048aNIzIykvDwcMLDw4mMjGTChAlX9fOK1DSlLanPjP6TqLBHSPt9O9jY4tF3Ej43P4uts3up17y08rSISHVgMqrR/hvnz5/H19eXbdu2ceONN2IYBoGBgYSGhjJz5kwgfzTIz8+PBQsWMHnyZJKSkvDx8WHp0qXceuutAJYyAN988w2DBw/mt99+Izg4mIiICLp16wZAREQEISEh/P7777Rq1arMviUnJ2M2m0lKSsLdvfR/8EVqi11/xXH7+xGFjhmGQcq+9SRs/RByc7B198Vn1BM4Nri21GvVd7HnhZvaqvK0iFxV1n5/V6sdEZOSkgDw9MwfTj927BjR0dEMGjTI0sbR0ZHevXuzc+dOAPbu3Ut2dnahNoGBgbRp08bSZteuXZjNZkswBNC9e3fMZrOljYgUdWnyc17GBWJXzyfh28WQm4Nzy+4E3P1GmcHQ0DZ+7P3nQAVDIlJtVWhz1yvBMAymT5/ODTfcQJs2bQCIjo4GwM/Pr1BbPz8/Tpw4YWnj4OCAh4dHkTYF50dHR+PrW7QYnK+vr6XNpTIzM8nMzLS8Tk5OruAnE6m5Lk5+zjx7hPNrF5KbdA5s7PDoOwm3ziMxmcqe/rozJEjTZCJSrVWbgOihhx7i119/ZceOHUXeu/QfXMMwyvxH+NI2xbUv7Trz589n7ty51nRdpNbqGuSJv7sj//12BQnbwiAvF7v6/niPmoljQMsyzzeRv0mrkqhFpLqrFlNmDz/8MGvXrmXr1q00bNjQctzf3x+gyChOTEyMZdTI39+frKwsEhISSm1z7ty5Ivc9f/58kdGnAk8++SRJSUmWn1OnTlX8A4rUUEmJCbB5YX6+UF4uLq16EjDxdauDIVAStYjUDFUaEBmGwUMPPcSqVavYsmULQUFBhd4PCgrC39+fzZs3W45lZWWxbds2evToAUDnzp2xt7cv1CYqKoqDBw9a2oSEhJCUlMTu3bstbX766SeSkpIsbS7l6OiIu7t7oR+RumTnzp106NCBn77fjMnOHs9BU/EePQsbR9eyTyZ/ZOid8Z2UNyQiNUKVTpk9+OCDfPbZZ6xZswY3NzfLSJDZbMbZ2RmTyURoaCjz5s2jZcuWtGzZknnz5uHi4sK4ceMsbe+55x4ee+wxvLy88PT0ZMaMGbRt25YBAwYA0Lp1a4YMGcJ9993H4sWLAbj//vsZMWKEVSvMRGqzrJw8lu46zon4NJp4unBHt8a88dorPPXUU+Tm5tKoaTNy+oTi4Nes1Os8Paw1wQHuxKZm4uuWP02mkSERqSmqNCB65513AOjTp0+h40uWLGHixIkAPPHEE6SnpzN16lQSEhLo1q0bmzZtws3NzdL+1Vdfxc7OjrFjx5Kenk7//v0JCwvD1tbW0ubTTz9l2rRpltVoo0aNYtGiRVf2A4pUc/O/Ocz724+R97/iG7lpSTw8cSzpR/cC+fW7Rk59llnr/izzWlFJ6dx3Y+lBk4hIdVWt6hBVZ6pDJLXN/G8Os/iHY5bXGacOErt2IbkX4jHZOTBm6tN89dozRByNL1KLqDherg7sfnqARoVEpFqx9vu72qwyE5Er6+LNWT1dHHjvf8GQkZdLUsQXJO34DIw87Dwb4jNmFpEuTcnONega5Imnqz3xqdmlXj8uNUu72ItIjaWASKQOKGlz1tzUBGLXvUzGiUgAXNv0x3PgFGwcnMgzYOmu49zTqxn/16EBH/54vMz7aBd7EampFBCJ1HIFm7NeOjeefuIX4ta9RG5qAiZ7RzwHTqVe2/6F2pyITwNgQLC/VQGRdrEXkZpKAZFILZaVk8dTXx8sFAwZebkk/biCpJ0rAAN778Z4j56Fg3fjIuc38XQB8gs0BpidiE7KKHajVxVgFJGarloUZhSRyhd+MIru878lPjXLciwnJY5zK/9J0s7lgEG9doPwv/OVYoMhAF/3/BEfWxsTs0cGA38XXCygAowiUhsoIBKphQqmyS5OhE4/to+osGlknjyAycEZ75Ez8Bo6DRv7kqe5pi3fT/jBKACGtAngnfGd8DcXbq8CjCJSG2jKTKSWyc0zmLvusGVqy8jLJXHHpyTv+gIwsPcNwmf0LOw9G1h1vbnrDjMw2B9bGxND2gQwMNjfslpNBRhFpLZQQCRSgxUspT+bkEbk6UTyJ7AMy2qynORYYtctJPP0YQDqdRyGZ797Mdk5WHV9A4hKyiDsx2NM7Jm/Y72tjUlL60Wk1lFhRiupMKNUNyUtpS+Q/tceYje8Sl56MiYHZ7yGTMO1da8K3y/A7MTskcGaGhORGsXa72/lEInUQAU5QsUFQ0ZuDglbPyLmy7nkpSfj4N+CgIlvXFYwBPkjRVOW7bPkFImI1CaaMhOpYS7NEbpYTlIM59cuIOvsEQDcOo/Eo88kTHb2JV7PxgSGQbHXu5RB4ZwiEZHaQiNEIjXM7mPxxY4Mpf0RQVTYNLLOHsHG0RWf/3sKzwGTSwyGTP/7ua9XkFXBUIGopAx2H4uvUN9FRKorjRCJ1DCXbo9h5GaTsHUJKXvXAuAQcA0+o2diZ/Yr9Tr+F+UEZecafGRFJeqS+iAiUtMpIBKpYY7Hpln+nJ0YTeyaBWRF/wGAe5f/o37vOzHZ/j0qNKF7Y65v6om3qyOYIPZCZpHl8gOD/csVEGmLDhGpbRQQidQguXkGy3efBCD1yI/EffM6RlYaNk5ueA0PxaVFt0LtbUzwzIjrcLArfXa8rK05LhagLTpEpBZSQCRSTRXUGLq4AOLuY/FExSeTsPVDUvZtAMCxQWu8Rz2OnbtvkWvc1yuozGAI/t6aY8qyfaW2M6EtOkSkdlJAJFINFVdjKMDsRDfPTKKXPU7Wub8AcO9+C/VvGI/Jtuh/yiPb+fPksGCr71mwNUdJtY1Uh0hEajMFRCLVTEGNoUunrv7c9R9++s8ijKx0bJzd8R7xGM7NOpd4nQHB/uW+98Vbc0QnZxB/IRNPVwf8zc7aokNEajUFRCLVSHE1hvKyM0n47n0u/BIOgGOjNniPnIGdm3ep16po4rO25hCRukgBkUg1cmmNoey4U5xfs4Ds88cBE+aQsZhvGIfJxrbU6yjxWUSkfBQQiVSR4pKmL67vc+HgFuI3vY2RnYGNa328R8zAuWkHAPpd68OW388Xe10lPouIlJ8CIpEqUFLS9PVN6pOXlUH8t++SeuBbAJyatMN7xOPY1vOwtL2vV3Nu6dSIf645SHxqVqFrKPFZRKT8FBCJXGUlJU1HJWXw1Xc/EbtmAdlxJ8Fkg7nn7ZhDxlqmyEzkV5guSHAe3Ma/yCiTRoZERMpPAZHIVVTSxqyGYZB6YDPxmxdj5GRiW88T75EzcGrcrsg1nhneWkGQiEglU0AkchUVtzFrXlY68ZveJvXQVgCcmnbEe8Rj2LrWL3L+8HYB/GvDb0Wm2jRNJiJyebTbvchVtPlwdKHXWTHHiPo4ND8YMtlQv/dd+I6dW2wwBLD+16giAVV0UgZTlu0j/GDUleq2iEitpxEikUpU3MoxWxsTuXkGi7b8YdlA1TAMLvwSTvy370FuNrZu3niPehynhteV+54G+blFc9cdZmCwv6bPREQqQAGRSCUpaeXYqPYBrImMIjo5/3heZhpx4W+S9vt2AJybd8FrWCi2LuYK39sgPyl797F4FVUUEakABUQilaC0lWOLfzhmeZ0Z/SexaxeQkxAFNrbUv/Eu3LuOwWQqefbaBGXuQF/g4jpGIiJiPeUQiVymklaOXcwwDJL3riN62QxyEqKwdffFf9wCzN1uKhQM9b/WhwBz4S03/M1OPDqgpVV9qeh2HSIidZ1GiEQuU3Erxy6Wl3GBuI1vkPbfnQA4t+yeP0XmVK9I23t7NadrkGeRPCSAFXtOEZ2UUWzgdXF9IhERKT8FRCKXqbRpqsyzRzi/diG5SefAxg6Pvnfj1nkUJlPhxOdLCy4Wlwc0e2QwU5btKzKFZrrofSVUi4hUjKbMRC5TcdNUhmGQvGc10Z/OJDfpHHZmP/zHL8T9+tFFgqECZQU0Q9oE8M74TvgXM6X2zvhOqkMkInIZNEIkUg65eQYRR+PY9VccYBDSzJsuQZ4EmJ0s01m56SnEffMa6X/+BIBLq554DZ2GjaNrsdf0d3dkzqjrrApohrQJYGCwtusQEalsCohErBR+MIpZqw6QmJZtObZo61/Ud7GnS1MPopIyyDzzG+fXLCQ35TzY2uPZ/z7qdRha4qjQowOu4aF+LcoV0JQ0pSYiIhWngEjkf0oqqgj5wdADy/YVe15iWjabDkWTvHsVids+ASMPO49AfEbPwsGvGQODfTl4JlnbbYiIVGMKiEQouaji7JHBDAz2Z87awyWem5uWROyGV8g4uhcAl9a98Rr8IDaOLgAcPJPMtsf7svdEgqa5RESqKQVEUueVVFSxYI+w0AHXWKpMXyrj1EFi175I7oU4THYOePS/n3rtBxeaIotKymDviQRNc4mIVGMKiKROK62oYsEeYR/9eLToe0Yeybu+IHHHp/lTZJ4N8RkzCwefpsXeRxWkRUSqNwVEUqeVVVTRAJLScwody01NIHb9K2Qc3w+Aa5t+eA6cgo2Dc4nXUQVpEZHqTQGR1GnlHbnJOPErseteJDc1AZO9I54Dp1Cv7YAS26uCtIhIzaCASOo0a0dujLxcknauJGnnCjDysPdujPfoWTh4Ny7xHFWQFhGpORQQSZ3W9ZKiisXJuRBP7LqXyDz5KwD12g3CY8D92Nj/HUy5Otri5mhHdHKm5Zi/ltaLiNQYCoikTrO1MZW4RxhA+rH9xK5/mby0REz2TngOfpB61/Utcp2X/9FeFaRFRGowBURS5xXsEXZxHSIjL5eMn1YQ88MKwMDeNwif0bOw92xQ5PxJPZtaRoG0tF5EpGZSQCR1VkFl6ujkDGJTMhjTPpCo5AxS4s6x9d3ZxBzML7RYr8NQPPrdi429Y7HXGRjsfzW7LSIiV4ACIqmTiqtMDZD+1x5iN7xKXnoyJgdnvIdOw+XaXsVeQyvIRERqDwVEUucUV5nayM0h8YdPSN69CgAHv+Z4j56JvUdgsdfQCjIRkdpFAZHUCbl5BhFH49j5ZyxhO48XCoZykmOIXbOQzLO/A+DWeSQefSZhsrO3tLExQd5FJ2kFmYhI7aKASGq98INRzFp1gMS07CLvpf3xE3HfvEpexgVMjq54D30El1Y9irTLM+CZ4a3xdnPUCjIRkVpIAZHUauEHo3hg2b4ix43cbBK+DyPl5zUAOAS0xHvUTOzrl5wg7e3myOgORVeZiYhIzaeASGqt3DyDOWsPFzmenRhN7NoFZEX9AYBblzF49L4Lk619kbYX035kIiK1lwIiqbUKltRfLO3ITmI3vo6RmYqNUz28hj+KS4tuZV4rQKvJRERqNQVEUmtdvHGrkZNFwtaPSNm3HgDHBq3xHvU4du6+ZV7HhFaTiYjUdjZVefMffviBkSNHEhgYiMlkYvXq1YXeNwyDOXPmEBgYiLOzM3369OHQoUOF2mRmZvLwww/j7e2Nq6sro0aN4vTp04XaJCQkMGHCBMxmM2azmQkTJpCYmHiFP51UtYIpruyEs0Qve9wSDLl3uwW/2+dbFQwFmJ14Z3wnrSYTEanlqjQgSk1NpX379ixatKjY9xcuXMgrr7zCokWL2LNnD/7+/gwcOJCUlBRLm9DQUL7++mtWrFjBjh07uHDhAiNGjCA3N9fSZty4cURGRhIeHk54eDiRkZFMmDDhin8+qVpdgzyxP76LqLBHyDr3FzbO7vjeMgePPhMx2ZY8OFrP0ZaX/9Ge5fd1Z8fMfgqGRETqAJNhGCVt8n1VmUwmvv76a8aMGQPkjw4FBgYSGhrKzJkzgfzRID8/PxYsWMDkyZNJSkrCx8eHpUuXcuuttwJw9uxZGjVqxDfffMPgwYP57bffCA4OJiIigm7d8nNFIiIiCAkJ4ffff6dVq1ZW9S85ORmz2UxSUhLu7u6V/wCkUqWnp/Poo4+yePFiABwbXpc/RebmXeI5BRNiGhESEak9rP3+rtIRotIcO3aM6OhoBg0aZDnm6OhI79692blzJwB79+4lOzu7UJvAwEDatGljabNr1y7MZrMlGALo3r07ZrPZ0qY4mZmZJCcnF/qRmuHIkSN0796dxYsXYzKZuP3+UFpNerFIMHRpRpC/psdEROqsaptUHR0dDYCfn1+h435+fpw4ccLSxsHBAQ8PjyJtCs6Pjo7G17doroivr6+lTXHmz5/P3LlzL+szyNW3bNkyHnjgAVJTU/H19WXZsmUMHDjQUql6119xgEFIM2+6BHmy90QCMSkZKrYoIlLHVduAqIDJVPgLyjCMIscudWmb4tqXdZ0nn3yS6dOnW14nJyfTqFEja7stV1laWhoPP/wwH330EQB9+/bl008/JSAgf7TH1sZEzxbe9GxReJQopLnXVe+riIhUP9V2yszfP79i8KWjODExMZZRI39/f7KyskhISCi1zblz54pc//z580VGny7m6OiIu7t7oR+png4fPkyXLl346KOPsLGxYe7cuWzevNkSDImIiJSl2gZEQUFB+Pv7s3nzZsuxrKwstm3bRo8e+XtNde7cGXt7+0JtoqKiOHjwoKVNSEgISUlJ7N6929Lmp59+IikpydJGaibDMFiyZAnXX389hw8fJiAggO+++45nn30WW1vbqu6eiIjUIFU6ZXbhwgX+/PNPy+tjx44RGRmJp6cnjRs3JjQ0lHnz5tGyZUtatmzJvHnzcHFxYdy4cQCYzWbuueceHnvsMby8vPD09GTGjBm0bduWAQMGANC6dWuGDBnCfffdZ1lxdP/99zNixAirV5hJ9XPhwgWmTp3K0qVLARg0aBBLly4tNl9MRESkTEYV2rp1qwEU+bnrrrsMwzCMvLw8Y/bs2Ya/v7/h6Oho3HjjjcaBAwcKXSM9Pd146KGHDE9PT8PZ2dkYMWKEcfLkyUJt4uLijDvuuMNwc3Mz3NzcjDvuuMNISEgoV1+TkpIMwEhKSrqcjyyV4JdffjFatWplAIatra0xb948Izc3t6q7JSIi1ZC139/Vpg5Rdac6RFXPMAzef/99HnnkETIyMmjQoAErVqzghhtuqOquiYhINWXt93e1X2UmAvm/0JMnT2bFihUADBs2jI8//hhv75ILLYqIiFir2iZVixTYv38/nTt3ZsWKFdjZ2fHiiy+ybt06BUMiIlJpNEIk1ZZhGLz99ttMnz6drKwsGjduzIoVKwgJCanqromISC2jgEiqpcTERO677z6+/PJLAEaNGsWSJUvw9PSs4p6JiEhtpCkzqXb27NlDp06d+PLLL7G3t+e1115j9erVCoZEROSK0QiRVBuGYfD666/zxBNPkJ2dTVBQECtXrqRLly5V3TUREanlFBBJtRAfH8+kSZNYs2YNADfffDMffPAB9evXr9qOiYhInaApM6lyERERdOzYkTVr1uDg4MBbb73FF198oWBIRESuGgVEUmXy8vJ48cUX6dWrFydPnqRFixZEREQwdepUTCZTVXdPRETqEE2ZSZWIjY3lrrvu4ptvvgHgtttuY/HixaoCLiIiVUIBkVx127dv5/bbb+fMmTM4OTnxxhtvcO+992pUSEREqoymzOSqycvLY968efTt25czZ87QqlUrfvrpJ+677z4FQyIiUqU0QiRXRUxMDBMmTGDTpk0ATJgwgbfffpt69epVcc9EREQUEMlV8P333zNu3DiioqJwdnbmrbfeYuLEiRoVEhGRakNTZnLF5ObmMnfuXPr3709UVBTBwcH8/PPP3H333QqGRESkWtEIkVwR0dHR3HHHHWzZsgWASZMm8eabb+Li4lLFPRMRESlKAZFUum+//ZY77riDmJgYXF1deffddxk/fnxVd0tERKREmjKTSpOTk8M///lPBg0aRExMDO3atWPv3r0KhkREpNrTCJFUijNnzjBu3Dh++OEHACZPnsyrr76Ks7NzFfdMRESkbAqI5LJt3LiRO++8k9jYWNzc3Hj//fe59dZbq7pbIiIiVtOUmVRYdnY2M2fOZNiwYcTGxtKpUyf27dunYEhERGocjRBJhZw8eZLbb7+dnTt3AvDQQw/x0ksv4ejoWMU9ExERKT8FRFJu69at46677iIhIQGz2cyHH37IzTffXNXdEhERqTBNmYnVsrKyeOyxxxg1ahQJCQl06dKF/fv3KxgSEZEaTwGRWOXYsWP06tWLV155BYBHH32UHTt2EBQUVMU9ExERuXyaMpMyrVq1ikmTJpGUlISHhwdhYWGMGjWqqrslIiJSaTRCJCXKzMzk4Ycf5uabbyYpKYmQkBAiIyMVDImISK2jgEiK9eeff9KjRw8WLVoEwBNPPMG2bdto3LhxFfdMRESk8mnKTIr4/PPPuffee0lJScHLy4tPPvmEYcOGVXW3RERErhiNEIlFeno6U6ZM4dZbbyUlJYVevXoRGRmpYEhERGo9BUQCwJEjR+jevTvvvvsuJpOJp59+mi1bttCwYcOq7pqIiMgVpykz4dNPP2Xy5Mmkpqbi6+vLsmXLGDhwYFV3S0RE5KrRCFEdlpaWxr333sv48eNJTU2lb9++REZGKhgSEZE6RwFRHXX48GG6du3Khx9+iMlkYs6cOWzevJmAgICq7pqIiMhVpymzOigsLIwHH3yQtLQ0/P39+eyzz+jbt29Vd0tERKTKaISoDrlw4QJ33XUXd999N2lpaQwcOJBffvlFwZCIiNR5CojqiAMHDtClSxc++eQTbGxs+Pe//014eDi+vr5V3TUREZEqpymzWs4wDD744AOmTZtGRkYGDRo0YPny5fTq1auquyYiIlJtKCCqxVJSUpg8eTLLly8HYOjQoXzyySd4e3tXcc9ERESqF02Z1VL79++nU6dOLF++HFtbWxYuXMj69esVDImIiBRDI0S1jGEYvPPOOzz66KNkZWXRuHFjVqxYQUhISFV3TUREpNpSQFSLJCUlce+99/Lll18CMGrUKJYsWYKnp2cV90xERKR605RZLfHzzz/TsWNHvvzyS+zt7Xn11VdZvXq1giEREREraISohjMMgzfeeIPHH3+c7OxsgoKCWLlyJV26dKnqromIiNQYCohqsISEBCZNmsTq1asBuPnmm/nggw+oX79+lfZLRESkptGUWQ0VERFBx44dWb16NQ4ODixatIgvvvhCwZCIiEgFKCCqYfLy8njppZfo1asXJ06coHnz5uzatYsHH3wQk8lU1d0TERGpkTRlVoPExcVx1113sWHDBgBuvfVW3nvvPdzd3au4ZyIiIjWbRohqiB07dtChQwc2bNiAo6MjixcvZvny5QqGREREKoEComouLy+P+fPn06dPH06fPk2rVq3YvXs3999/v6bIREREKommzKqxmJgYJkyYwKZNmwAYP34877zzDvXq1avinomIiNQuCoiqqe+//55x48YRFRWFs7Mzb731FhMnTtSokIiIyBWgKbNqJjc3l+eee47+/fsTFRVFcHAwe/bs4e6771YwJCIicoXUqYDo7bffJigoCCcnJzp37sz27durukuFREdHM2jQIGbPnk1eXh6TJk1iz549XHfddVXdNRERkVqtzgREK1euJDQ0lKeffpr9+/fTq1cvhg4dysmTJ6u6awB8++23dOjQgS1btuDq6sonn3zChx9+iIuLS1V3TUREpNYzGYZhVHUnroZu3brRqVMn3nnnHcux1q1bM2bMGObPn1/m+cnJyZjNZpKSkip1qXtOTg5z587l3//+N4Zh0LZtWz7//HOuvfbaSruHiIhIXWXt93edSKrOyspi7969zJo1q9DxQYMGsXPnzmLPyczMJDMz0/I6OTm50vuVnZ3NoEGD+P777wGYPHkyr776Ks7OzpV+LxERESlZnZgyi42NJTc3Fz8/v0LH/fz8iI6OLvac+fPnYzabLT+NGjWq9H7Z29vTuXNn3NzcWL58Oe+++66CIRERkSpQJwKiApeu0jIMo8SVW08++SRJSUmWn1OnTl2RPs2bN4/IyEhuu+22K3J9ERERKVudmDLz9vbG1ta2yGhQTExMkVGjAo6Ojjg6Ol7xvjk4ONCsWbMrfh8REREpWZ0YIXJwcKBz585s3ry50PHNmzfTo0ePKuqViIiIVBd1YoQIYPr06UyYMIHrr7+ekJAQ3nvvPU6ePMkDDzxQ1V0TERGRKlZnAqJbb72VuLg4nnvuOaKiomjTpg3ffPMNTZo0qequiYiISBWrM3WILteVqkMkIiIiV4613991IodIREREpDQKiERERKTOU0AkIiIidZ4CIhEREanzFBCJiIhInaeASEREROo8BUQiIiJS5ykgEhERkTpPAZGIiIjUeXVm647LVVDQOzk5uYp7IiIiItYq+N4ua2MOBURWSklJAaBRo0ZV3BMREREpr5SUFMxmc4nvay8zK+Xl5XH27Fnc3NwwmUyVdt3k5GQaNWrEqVOntEdaJdJzvTL0XK8MPdcrQ8/1yqhpz9UwDFJSUggMDMTGpuRMIY0QWcnGxoaGDRteseu7u7vXiF+smkbP9crQc70y9FyvDD3XK6MmPdfSRoYKKKlaRERE6jwFRCIiIlLnKSCqYo6OjsyePRtHR8eq7kqtoud6Zei5Xhl6rleGnuuVUVufq5KqRUREpM7TCJGIiIjUeQqIREREpM5TQCQiIiJ1ngIiERERqfMUEFWxt99+m6CgIJycnOjcuTPbt2+v6i5VGz/88AMjR44kMDAQk8nE6tWrC71vGAZz5swhMDAQZ2dn+vTpw6FDhwq1yczM5OGHH8bb2xtXV1dGjRrF6dOnC7VJSEhgwoQJmM1mzGYzEyZMIDEx8Qp/uqozf/58unTpgpubG76+vowZM4YjR44UaqNnW37vvPMO7dq1sxSrCwkJYePGjZb39Uwv3/z58zGZTISGhlqO6blWzJw5czCZTIV+/P39Le/XyedqSJVZsWKFYW9vb7z//vvG4cOHjUceecRwdXU1Tpw4UdVdqxa++eYb4+mnnza++uorAzC+/vrrQu+/8MILhpubm/HVV18ZBw4cMG699VYjICDASE5OtrR54IEHjAYNGhibN2829u3bZ/Tt29do3769kZOTY2kzZMgQo02bNsbOnTuNnTt3Gm3atDFGjBhxtT7mVTd48GBjyZIlxsGDB43IyEhj+PDhRuPGjY0LFy5Y2ujZlt/atWuNDRs2GEeOHDGOHDliPPXUU4a9vb1x8OBBwzD0TC/X7t27jaZNmxrt2rUzHnnkEctxPdeKmT17tnHdddcZUVFRlp+YmBjL+3XxuSogqkJdu3Y1HnjggULHrr32WmPWrFlV1KPq69KAKC8vz/D39zdeeOEFy7GMjAzDbDYb7777rmEYhpGYmGjY29sbK1assLQ5c+aMYWNjY4SHhxuGYRiHDx82ACMiIsLSZteuXQZg/P7771f4U1UPMTExBmBs27bNMAw928rk4eFhfPDBB3qmlyklJcVo2bKlsXnzZqN3796WgEjPteJmz55ttG/fvtj36upz1ZRZFcnKymLv3r0MGjSo0PFBgwaxc+fOKupVzXHs2DGio6MLPT9HR0d69+5teX579+4lOzu7UJvAwEDatGljabNr1y7MZjPdunWztOnevTtms7nO/D0kJSUB4OnpCejZVobc3FxWrFhBamoqISEheqaX6cEHH2T48OEMGDCg0HE918vzxx9/EBgYSFBQELfddhtHjx4F6u5z1eauVSQ2Npbc3Fz8/PwKHffz8yM6OrqKelVzFDyj4p7fiRMnLG0cHBzw8PAo0qbg/OjoaHx9fYtc39fXt078PRiGwfTp07nhhhto06YNoGd7OQ4cOEBISAgZGRnUq1ePr7/+muDgYMs//nqm5bdixQr27dvHnj17iryn39WK69atG5988gnXXHMN586d4/nnn6dHjx4cOnSozj5XBURVzGQyFXptGEaRY1Kyijy/S9sU176u/D089NBD/Prrr+zYsaPIe3q25deqVSsiIyNJTEzkq6++4q677mLbtm2W9/VMy+fUqVM88sgjbNq0CScnpxLb6bmW39ChQy1/btu2LSEhITRv3pyPP/6Y7t27A3XvuWrKrIp4e3tja2tbJEqOiYkpEpVLUQWrIUp7fv7+/mRlZZGQkFBqm3PnzhW5/vnz52v938PDDz/M2rVr2bp1Kw0bNrQc17OtOAcHB1q0aMH111/P/Pnzad++Pa+//rqeaQXt3buXmJgYOnfujJ2dHXZ2dmzbto033ngDOzs7y2fWc718rq6utG3blj/++KPO/r4qIKoiDg4OdO7cmc2bNxc6vnnzZnr06FFFvao5goKC8Pf3L/T8srKy2LZtm+X5de7cGXt7+0JtoqKiOHjwoKVNSEgISUlJ7N6929Lmp59+Iikpqdb+PRiGwUMPPcSqVavYsmULQUFBhd7Xs608hmGQmZmpZ1pB/fv358CBA0RGRlp+rr/+eu644w4iIyNp1qyZnmslyczM5LfffiMgIKDu/r5e5SRuuUjBsvsPP/zQOHz4sBEaGmq4uroax48fr+quVQspKSnG/v37jf379xuA8corrxj79++3lCV44YUXDLPZbKxatco4cOCAcfvttxe7LLRhw4bGt99+a+zbt8/o169fsctC27VrZ+zatcvYtWuX0bZt22q7LLQyTJkyxTCbzcb3339faMltWlqapY2ebfk9+eSTxg8//GAcO3bM+PXXX42nnnrKsLGxMTZt2mQYhp5pZbl4lZlh6LlW1GOPPWZ8//33xtGjR42IiAhjxIgRhpubm+X7py4+VwVEVeytt94ymjRpYjg4OBidOnWyLH0Ww9i6dasBFPm56667DMPIXxo6e/Zsw9/f33B0dDRuvPFG48CBA4WukZ6ebjz00EOGp6en4ezsbIwYMcI4efJkoTZxcXHGHXfcYbi5uRlubm7GHXfcYSQkJFylT3n1FfdMAWPJkiWWNnq25Tdp0iTLf8s+Pj5G//79LcGQYeiZVpZLAyI914opqCtkb29vBAYGGjfddJNx6NAhy/t18bmaDMMwqmZsSkRERKR6UA6RiIiI1HkKiERERKTOU0AkIiIidZ4CIhEREanzFBCJiIhInaeASEREROo8BUQiIiJS5ykgEhGpRCaTidWrV1d1N0SknBQQiUiNtXPnTmxtbRkyZEi5zmvatCmvvfbalemUiNRICohEpMb66KOPePjhh9mxYwcnT56s6u6ISA2mgEhEaqTU1FQ+//xzpkyZwogRIwgLCyv0/tq1a7n++utxcnLC29ubm266CYA+ffpw4sQJHn30UUwmEyaTCYA5c+bQoUOHQtd47bXXaNq0qeX1nj17GDhwIN7e3pjNZnr37s2+ffuu5McUkatEAZGI1EgrV66kVatWtGrVivHjx7NkyRIKtmbcsGEDN910E8OHD2f//v189913XH/99QCsWrWKhg0b8txzzxEVFUVUVJTV90xJSeGuu+5i+/btRERE0LJlS4YNG0ZKSsoV+YwicvXYVXUHREQq4sMPP2T8+PEADBkyhAsXLvDdd98xYMAA/v3vf3Pbbbcxd+5cS/v27dsD4Onpia2tLW5ubvj7+5frnv369Sv0evHixXh4eLBt2zZGjBhxmZ9IRKqSRohEpMY5cuQIu3fv5rbbbgPAzs6OW2+9lY8++giAyMhI+vfvX+n3jYmJ4YEHHuCaa67BbDZjNpu5cOGC8pdEagGNEIlIjfPhhx+Sk5NDgwYNLMcMw8De3p6EhAScnZ3LfU0bGxvLlFuB7OzsQq8nTpzI+fPnee2112jSpAmOjo6EhISQlZVVsQ8iItWGRohEpEbJycnhk08+4eWXXyYyMtLy88svv9CkSRM+/fRT2rVrx3fffVfiNRwcHMjNzS10zMfHh+jo6EJBUWRkZKE227dvZ9q0aQwbNozrrrsOR0dHYmNjK/XziUjV0AiRiNQo69evJyEhgXvuuQez2VzovVtuuYUPP/yQV199lf79+9O8eXNuu+02cnJy2LhxI0888QSQX4fohx9+4LbbbsPR0RFvb2/69OnD+fPnWbhwIbfccgvh4eFs3LgRd3d3y/VbtGjB0qVLuf7660lOTubxxx+v0GiUiFQ/GiESkRrlww8/ZMCAAUWCIYCbb76ZyMhI3N3d+eKLL1i7di0dOnSgX79+/PTTT5Z2zz33HMePH6d58+b4+PgA0Lp1a95++23eeust2rdvz+7du5kxY0ah63/00UckJCTQsWNHJkyYwLRp0/D19b2yH1hErgqTcemkuYiIiEgdoxEiERERqfMUEImIiEidp4BIRERE6jwFRCIiIlLnKSASERGROk8BkYiIiNR5CohERESkzlNAJCIiInWeAiIRERGp8xQQiYiISJ2ngEhERETqPAVEIiIiUuf9P+O50UXL8cM2AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Finding the hyperparameters\n",
    "\n",
    "#set this so we can maintain the train test split in the actual run of the model\n",
    "random_state = 29\n",
    "\n",
    "#reach number to use \n",
    "rn = 32\n",
    "\n",
    "X_train, X_test, y_train, y_test = train_test_split(well_head_noRiv, streamflow[rn], test_size=0.2, random_state = random_state)\n",
    "\n",
    "#split quality function\n",
    "split_criterion = ['squared_error', 'friedman_mse']\n",
    "\n",
    "#Number of features to consider at every split\n",
    "max_features = [1.0]\n",
    "\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [int(x) for x in range(4,9)] #14 was the max depth, 3 or greater is recommened, 8 seems high enough \n",
    "#max_depth.append(None)\n",
    "\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [2, 5, 10, 15]\n",
    "\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [1, 2, 5, 10]\n",
    "\n",
    "#maximum \n",
    "maxSqrt = int(np.max(y_train)**(1/2))\n",
    "max_leaf_nodes = [int(x) for x in range(25,maxSqrt,int(maxSqrt/5))]\n",
    "max_leaf_nodes.append(None)\n",
    "\n",
    "#min weight fraction leaf - need to look into this one \n",
    "#min_weight_fraction_leaf = [0.0]\n",
    "\n",
    "min_imp_decrease = [int(x) for x in range(5,55,10)]\n",
    "#min_imp_decrease.append(0)\n",
    "#cost complexcity prunning \n",
    "ccp_alpha = [0.0]\n",
    "\n",
    "splitter = ['best']\n",
    "# Create the random grid\n",
    "random_grid = {#'n_estimators': n_estimators,\n",
    "                'criterion' : split_criterion,\n",
    "               #'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf,\n",
    "              # 'min_weight_fraction_leaf': min_weight_fraction_leaf,\n",
    "                'max_leaf_nodes': max_leaf_nodes,\n",
    "                'min_impurity_decrease':min_imp_decrease}\n",
    "                #'splitter': splitter}\n",
    "print(random_grid)\n",
    "\n",
    "#looking at how many possible parameter combinations there are \n",
    "num_par_combos = 1\n",
    "\n",
    "for i in random_grid.keys():\n",
    "    num_par_combos *= len(random_grid[i])\n",
    "print(\"Number of possible parameter combinations: \", num_par_combos )  \n",
    "\n",
    "# Use the random grid to search for best hyperparameters\n",
    "# First create the base model to tune\n",
    "dtr = DecisionTreeRegressor()\n",
    "\n",
    "# Random search of parameters, using 3 fold cross validation, \n",
    "# search across 100 different combinations, and use all available cores\n",
    "dtr_random = RandomizedSearchCV(estimator = dtr, scoring = 'neg_root_mean_squared_error',param_distributions = random_grid, n_iter = 100, cv = 3, verbose=2, random_state=29, n_jobs = -1)\n",
    "\n",
    "# Fit the random search model\n",
    "dtr_random.fit(X_train, y_train)\n",
    "\n",
    "best_pars_dt = dtr_random.best_params_\n",
    "print(\"Best DT pars: \", best_pars_dt) # This is the most important part\n",
    "\n",
    "best_random = dtr_random.best_estimator_\n",
    "y_pred = best_random.predict(X_test)\n",
    "#well head Tree NSE\n",
    "\n",
    "#label='Prediction to Test Values for Hyperparameter Tuned Tree'\n",
    "plt.scatter(y_test, y_pred)  \n",
    "outmax = np.max(y_test)\n",
    "plt.plot([0, outmax], [0, outmax], color='k', linestyle='-')\n",
    "plt.xlabel(\"Actual\")\n",
    "plt.ylabel(\"Predicted\")\n",
    "plt.title(\"Decision Tree Regression\")\n",
    "# plt.xlim(60,100)\n",
    "# plt.ylim(60,100)\n",
    "#plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# _, ax = plt.subplots(figsize=(25, 25))\n",
    "# _ = plot_tree(best_random, ax=ax)#, feature_names=col_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fb771112",
   "metadata": {},
   "source": [
    "### Gradient Boosting Hyperparameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "3a408202",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['alpha', 'ccp_alpha', 'criterion', 'init', 'learning_rate', 'loss', 'max_depth', 'max_features', 'max_leaf_nodes', 'min_impurity_decrease', 'min_samples_leaf', 'min_samples_split', 'min_weight_fraction_leaf', 'n_estimators', 'n_iter_no_change', 'random_state', 'subsample', 'tol', 'validation_fraction', 'verbose', 'warm_start'])\n",
      "{'loss': ['squared_error', 'absolute_error', 'huber', 'quantile'], 'subsample': [0.5, 0.6, 0.7, 0.8, 0.9], 'n_estimators': [25, 75, 125, 175, 225], 'criterion': ['squared_error', 'friedman_mse'], 'max_depth': [4, 5, 6, 7, 8], 'min_samples_split': [2, 5, 10, 15], 'min_samples_leaf': [1, 2, 5, 10], 'learning_rate': [0.01, 0.1, 0.2, 0.5], 'max_leaf_nodes': [50, None], 'min_impurity_decrease': [5, 15, 25, 35, 45]}\n",
      "Number of possible parameter combinations:  640000\n",
      "Fitting 3 folds for each of 100 candidates, totalling 300 fits\n",
      "{'subsample': 0.7, 'n_estimators': 75, 'min_samples_split': 5, 'min_samples_leaf': 5, 'min_impurity_decrease': 5, 'max_leaf_nodes': None, 'max_depth': 5, 'loss': 'squared_error', 'learning_rate': 0.2, 'criterion': 'friedman_mse'}\n"
     ]
    }
   ],
   "source": [
    "new_zone_q = ReachQdf.loc[well_head_noRiv.index]\n",
    "rn =25\n",
    "print(GradientBoostingRegressor().get_params().keys())\n",
    "#choosing one well and finding the hyperparameters that best work for it \n",
    "rand_reach = np.random.randint(low=1,high=48)\n",
    "#print(rand_well)\n",
    "X_train, X_test, y_train, y_test = train_test_split(well_head_noRiv, streamflow[rn], test_size=0.2,random_state=29)\n",
    "loss = ['squared_error', 'absolute_error', 'huber', 'quantile']\n",
    "n_estimators = [int(x) for x in range(25,275,50)] #this was found by looking at oob predictions - actually could have been lower than this\n",
    "subsample = [0.5,0.6,0.7,0.8,0.9]             \n",
    "#split quality function\n",
    "split_criterion = ['squared_error', 'friedman_mse']\n",
    "splitter = ['best']\n",
    "#Number of features to consider at every split\n",
    "max_features = [1.0, 'sqrt','log2']\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [int(x) for x in range(4,9)] #14 was the max depth, 3 is recommened, 10 seems high enough \n",
    "#max_depth.append(None)\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [2, 5, 10, 15]\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [1, 2, 5, 10]\n",
    "#maximum \n",
    "maxSqrt = int(np.max(y_train)**(1/2))\n",
    "max_leaf_nodes = [int(x) for x in range(50,maxSqrt,int(maxSqrt/5))]\n",
    "max_leaf_nodes.append(None)\n",
    "#min weight fraction leaf - need to look into this one \n",
    "#min_weight_fraction_leaf = [0.0]\n",
    "learning_rate = [0.01,0.1,0.2,0.5]\n",
    "min_imp_decrease = [int(x) for x in range(5,55,10)]\n",
    "#min_imp_decrease.append(0)\n",
    "#cost complexcity prunning \n",
    "ccp_alpha = [0.0]\n",
    "\n",
    "# Create the random grid\n",
    "random_grid = {#'n_estimators': n_estimators,\n",
    "                'loss' : loss,\n",
    "                'subsample':subsample,\n",
    "                'n_estimators':n_estimators,\n",
    "                'criterion' : split_criterion,\n",
    "              # 'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf,\n",
    "              # 'min_weight_fraction_leaf': min_weight_fraction_leaf,\n",
    "                'learning_rate': learning_rate,\n",
    "                'max_leaf_nodes': max_leaf_nodes,\n",
    "                'min_impurity_decrease':min_imp_decrease}\n",
    "print(random_grid)\n",
    "#looking at how many possible parameter combinations there are \n",
    "num_par_combos = 1\n",
    "\n",
    "for i in random_grid.keys():\n",
    "    num_par_combos *= len(random_grid[i])\n",
    "print(\"Number of possible parameter combinations: \", num_par_combos )  \n",
    "\n",
    "# Use the random grid to search for best hyperparameters\n",
    "# First create the base model to tune\n",
    "xgr = GradientBoostingRegressor()\n",
    "# Random search of parameters, using 3 fold cross validation, \n",
    "# search across 100 different combinations, and use all available cores\n",
    "xgr_random = RandomizedSearchCV(estimator = xgr, scoring = 'neg_root_mean_squared_error',param_distributions = random_grid, n_iter = 100, cv = 3, verbose=2, random_state=42, n_jobs = -1)\n",
    "# Fit the random search model\n",
    "xgr_random.fit(X_train, y_train)\n",
    "best_pars_gr = xgr_random.best_params_\n",
    "print(best_pars_gr)\n",
    "\n",
    "#from sklearn.tree import plot_tree\n",
    "\n",
    "\n",
    "best_random_gr = xgr_random.best_estimator_\n",
    "\n",
    "\n",
    "#_, ax = plt.subplots(figsize=(25, 25))\n",
    "#_ = plot_tree(best_random, ax=ax)#, feature_names=col_names)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e39ce142",
   "metadata": {},
   "source": [
    "### Random Forest Hyperparameter tuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a6617042",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "dict_keys(['bootstrap', 'ccp_alpha', 'criterion', 'max_depth', 'max_features', 'max_leaf_nodes', 'max_samples', 'min_impurity_decrease', 'min_samples_leaf', 'min_samples_split', 'min_weight_fraction_leaf', 'n_estimators', 'n_jobs', 'oob_score', 'random_state', 'verbose', 'warm_start'])\n",
      "{'max_samples': [0.1, 0.3, 0.5, 0.7, 0.9], 'n_estimators': [25, 75, 125, 175, 225], 'criterion': ['squared_error', 'friedman_mse'], 'max_depth': [4, 5, 6, 7, 8], 'min_samples_split': [2, 5, 10, 15], 'min_samples_leaf': [1, 2, 5, 10], 'max_leaf_nodes': [50, None], 'min_impurity_decrease': [5, 15, 25, 35, 45]}\n",
      "Number of possible parameter combinations:  40000\n",
      "Fitting 3 folds for each of 100 candidates, totalling 300 fits\n",
      "{'n_estimators': 25, 'min_samples_split': 2, 'min_samples_leaf': 1, 'min_impurity_decrease': 45, 'max_samples': 0.9, 'max_leaf_nodes': None, 'max_depth': 7, 'criterion': 'friedman_mse'}\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "rn = 25\n",
    "random_state = 29\n",
    "new_zone_q = ReachQdf.loc[well_head_noRiv.index]\n",
    "\n",
    "print(RandomForestRegressor().get_params().keys())\n",
    "#choosing one well and finding the hyperparameters that best work for it \n",
    "rand_reach = np.random.randint(low=1,high=48)\n",
    "#print(rand_well)\n",
    "X_train, X_test, y_train, y_test = train_test_split(well_head_noRiv, streamflow[rn], test_size=0.2, random_state = random_state)\n",
    "n_estimators = [int(x) for x in range(25,275,50)]\n",
    "max_sample = [0.1,0.3,0.5,0.7,0.9]             \n",
    "#split quality function\n",
    "split_criterion = ['squared_error', 'friedman_mse']\n",
    "splitter = ['best']\n",
    "bootstrap = ['True']\n",
    "#Number of features to consider at every split\n",
    "max_features = [1.0]\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [int(x) for x in range(4,9)] #14 was the max depth, 3 is recommened, 10 seems high enough \n",
    "#max_depth.append(None)\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [2, 5, 10,15]\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [1, 2, 5, 10]\n",
    "#maximum \n",
    "maxSqrt = int(np.max(y_train)**(1/2))\n",
    "max_leaf_nodes = [int(x) for x in range(50,maxSqrt,int(maxSqrt/5))]\n",
    "max_leaf_nodes.append(None)\n",
    "\n",
    "#min weight fraction leaf - need to look into this one \n",
    "#min_weight_fraction_leaf = [0.0]\n",
    "\n",
    "min_imp_decrease = [int(x) for x in range(5,55,10)]\n",
    "#min_imp_decrease.append(0)\n",
    "#cost complexcity prunning \n",
    "ccp_alpha = [0.0]\n",
    "\n",
    "# Create the random grid\n",
    "random_grid = {#'n_estimators': n_estimators,\n",
    "                #'bootstrap' : bootstrap,\n",
    "                'max_samples':max_sample,\n",
    "                'n_estimators':n_estimators,\n",
    "                'criterion' : split_criterion,\n",
    "               #'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf,\n",
    "              # 'min_weight_fraction_leaf': min_weight_fraction_leaf,\n",
    "                'max_leaf_nodes': max_leaf_nodes,\n",
    "                'min_impurity_decrease':min_imp_decrease}\n",
    "print(random_grid)\n",
    "\n",
    "#looking at how many possible parameter combinations there are \n",
    "num_par_combos = 1\n",
    "\n",
    "for i in random_grid.keys():\n",
    "    num_par_combos *= len(random_grid[i])\n",
    "    \n",
    "print(\"Number of possible parameter combinations: \", num_par_combos ) \n",
    "\n",
    "# Use the random grid to search for best hyperparameters\n",
    "# First create the base model to tune\n",
    "rfr = RandomForestRegressor()\n",
    "# Random search of parameters, using 3 fold cross validation, \n",
    "# search across 100 different combinations, and use all available cores\n",
    "rfr_random = RandomizedSearchCV(estimator = rfr, scoring = 'neg_root_mean_squared_error',param_distributions = random_grid, n_iter = 100, cv = 3, verbose=2, random_state=42, n_jobs = -1)\n",
    "# Fit the random search model\n",
    "rfr_random.fit(X_train, y_train)\n",
    "best_pars_rf = rfr_random.best_params_\n",
    "print(best_pars_rf)\n",
    "\n",
    "#from sklearn.tree import plot_tree\n",
    "\n",
    "\n",
    "best_random_rf = rfr_random.best_estimator_\n",
    "y_pred = best_random_rf.predict(X_test)\n",
    "\n",
    "#_, ax = plt.subplots(figsize=(25, 25))\n",
    "#_ = plot_tree(best_random, ax=ax)#, feature_names=col_names)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6989a647",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}