AWS-S3-Dask-XGBoost.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Overview\" data-toc-modified-id=\"Overview-1\">Overview</a></span></li><li><span><a href=\"#Load-libraries\" data-toc-modified-id=\"Load-libraries-2\">Load libraries</a></span></li><li><span><a href=\"#Load-Data\" data-toc-modified-id=\"Load-Data-3\">Load Data</a></span></li><li><span><a href=\"#Train-the-XGB-model\" data-toc-modified-id=\"Train-the-XGB-model-4\">Train the XGB model</a></span><ul class=\"toc-item\"><li><span><a href=\"#Determine-feature-importance\" data-toc-modified-id=\"Determine-feature-importance-4.1\">Determine feature importance</a></span></li></ul></li><li><span><a href=\"#Test-model-performance\" data-toc-modified-id=\"Test-model-performance-5\">Test model performance</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Overview"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Author**:\n",
    "**Zhonghua Zheng** (zzheng25@illinois.edu)   \n",
    "GitHub: **zzheng93**   \n",
    "\n",
    "**Motivation**:\n",
    "This notebook shows the workflow of \n",
    "* Loading Climate Data ([CESM-LE from AWS S3](https://medium.com/pangeo/cesm-lens-on-aws-4e2a996397a1))\n",
    "* Scalable [Dask-XGBoost](https://examples.dask.org/machine-learning/xgboost.html) *regression* implementation (the official documentation of [Dask-XGBoost](https://examples.dask.org/machine-learning/xgboost.html) only provides the *classfication* implementation.\n",
    "\n",
    "**Task**:\n",
    "How to predict **maximum reference height temperature (TREFHTMX)** based on related features?\n",
    "\n",
    "* Response (Y):   \n",
    "\"TREFHTMX\": Maximum reference height temperature over output period  \n",
    "\n",
    "* Features (X):   \n",
    "\"PRECT\": Total (convective and large-scale) precipitation rate (liq + ice)  \n",
    "\"WSPDSRFAV\": Horizontal total wind speed average at the surface  \n",
    "\"TS\": Surface temperature (radiative)  \n",
    "\"TREFHT\": Reference height temperature    \n",
    "\n",
    "**Prerequisite**  \n",
    "How to create environment for this workflow?  \n",
    "Recommendation: \n",
    "```bash\n",
    "git clone http://github.com/dask/dask-tutorial\n",
    "cd dask-tutorial\n",
    "conda env create -f binder/environment.yml \n",
    "conda activate dask-tutorial\n",
    "conda install xarray cftime=1.0.3.4\n",
    "conda install -c conda-forge dask-ml dask-xgboost\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we load the necessary libraries, connect to cluster, and set up anonymous access to s3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import s3fs\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import dask\n",
    "import dask.array as da\n",
    "import dask.dataframe as dd\n",
    "import dask_xgboost\n",
    "import matplotlib.pyplot as plt\n",
    "from dask.distributed import Client\n",
    "from dask_ml.model_selection import train_test_split\n",
    "from sklearn import metrics\n",
    "import xgboost\n",
    "\n",
    "client = Client() # connect to cluster\n",
    "s3 = s3fs.S3FileSystem(anon=True) # anonymous access to s3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style=\"border: 2px solid white;\">\n",
       "<tr>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Client</h3>\n",
       "<ul style=\"text-align: left; list-style: none; margin: 0; padding: 0;\">\n",
       "  <li><b>Scheduler: </b>tcp://127.0.0.1:34195</li>\n",
       "  <li><b>Dashboard: </b><a href='http://127.0.0.1:8787/status' target='_blank'>http://127.0.0.1:8787/status</a>\n",
       "</ul>\n",
       "</td>\n",
       "<td style=\"vertical-align: top; border: 0px solid white\">\n",
       "<h3 style=\"text-align: left;\">Cluster</h3>\n",
       "<ul style=\"text-align: left; list-style:none; margin: 0; padding: 0;\">\n",
       "  <li><b>Workers: </b>4</li>\n",
       "  <li><b>Cores: </b>16</li>\n",
       "  <li><b>Memory: </b>32.66 GB</li>\n",
       "</ul>\n",
       "</td>\n",
       "</tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Client: 'tcp://127.0.0.1:34195' processes=4 threads=16, memory=32.66 GB>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "client"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load Data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we load the CESM-LE RCP8.5 data from [AWS S3](https://medium.com/pangeo/cesm-lens-on-aws-4e2a996397a1).   \n",
    "The original format of the data is [Zarr](https://zarr.readthedocs.io/en/stable/), but Dask-XGBoost only accepts Dask.array or Dask.dataframe collections.\n",
    "\n",
    "So we use [xarray](http://xarray.pydata.org/en/stable/) to load the data, and convert to [Dask Dataframe](https://docs.dask.org/en/latest/dataframe.html).  \n",
    "\n",
    "As an illustration, here I only use \"member_id=1\" and \"time=\"2010-01-01\"\" for this tutorial.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.55 s, sys: 308 ms, total: 6.85 s\n",
      "Wall time: 1min 34s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "zarr_ls = [xr.open_zarr(s3fs.S3Map(root=\"ncar-cesm-lens/atm/daily/cesmLE-RCP85-\"+var+\".zarr\", s3=s3))[var]\n",
    "           .sel(member_id=1,time=\"2010-01-01\") \n",
    "           for var in [\"TREFHT\",\"TREFHTMX\",\"PRECT\",\"TS\",\"WSPDSRFAV\"]]  # read data from s3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "dddf=xr.merge(zarr_ls).to_dask_dataframe()  # merge arrays and convert to dask dataframe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "X=dddf[[\"TREFHT\",\"PRECT\",\"TS\",\"WSPDSRFAV\"]]  # features for Dask-XGBoost\n",
    "Y=dddf[\"TREFHTMX\"]  # response for Dask-XGBoost"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.15)  # split into training and testing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train the XGB model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This parameters are available at [Learning Task Parameters](https://xgboost.readthedocs.io/en/latest/parameter.html).\n",
    "* **reg:squarederror**: regression with squared loss.\n",
    "* **eta (learning_rate)**: Step size shrinkage used in update to prevents overfitting. After each boosting step, we can directly get the weights of new features, and eta shrinks the feature weights to make the boosting process more conservative.  \n",
    "* **max_depth**: Maximum depth of a tree. Increasing this value will make the model more complex and more likely to overfit. 0 is only accepted in lossguided growing policy when tree_method is set as hist and it indicates no limit on depth. Beware that XGBoost aggressively consumes memory when training a deep tree. \n",
    "* **num_boost_round**:  [Number of boosting](iterationshttps://xgboost.readthedocs.io/en/latest/python/python_api.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 306 ms, sys: 57.8 ms, total: 364 ms\n",
      "Wall time: 13 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "params = {'objective': 'reg:squarederror',\n",
    "          'max_depth': 6, 'eta': 0.01}\n",
    "\n",
    "bst = dask_xgboost.train(client, params, X_train, y_train, num_boost_round=500)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Determine feature importance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, importance provides a score that indicates how useful or valuable each feature was in the construction of the boosted decision trees within the model. The more an attribute is used to make key decisions with decision trees, the higher its relative importance (reference: [Feature Importance and Feature Selection With XGBoost in Python](https://machinelearningmastery.com/feature-importance-and-feature-selection-with-xgboost-in-python/))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEWCAYAAAAHC8LZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3dfZxWdZ3/8dcbEMQAiXBYRkRCWVEBR8Gb1GzCUMxQLGMle4iakSu6628rpNVfUdZqmGnlritlhneorCne5U0i4ZqKWCNihsqNiaJAgtzfDZ/945zB4zg31+DMXJyL9/PxOI8553vuPt9r4HrPubnOpYjAzMwsb9oUuwAzM7Md4QAzM7NccoCZmVkuOcDMzCyXHGBmZpZLDjAzM8slB5jtEiR9WtL8YtdRF0mVkpY0MP8YSa9KWitpZGvW9lHtzK+75Z8DzHZqkhZL2pC+edcM1xWwXkjav2Y6Ip6MiANaqMbfSPphS2w79QPguojoFBH3fpQNpa/n55qprka15OveVI39oWD5067YBZgVYERE/L7YRRTRvsBLxS4CQFK7iNha7DqaSpLf60qQj8AstyTtL+kPkt6TtELSnWn7rHSRF9Ijtn+q/dd3eiTybUlzJa2TdKOkHpJ+J2mNpN9L+nhm+WmS3k73NUvSwWn7WOBMYHy6r/vT9nJJd0taLmmRpH/JbKtjetS2UtJfgMMb6OMCoC9wf7r9DpL2TOtdKulNST+U1DZdfj9JMyT9PX1NbpPUNZ13C9A7s63xdR2VZI/SJE2U9D+SbpW0GjhbUhtJEyQtSPdzl6Ru9dS/w6+7pD7pkfRYSW+l/f1mZlsdJF2bznsrHe+Q3a+kSyS9DUwFfgeUZ47kyyUdIelpSavS7V8nqX1mHyHpfCWncFdK+k9Jysz/uqSX09r/Iumwxn7/1owiwoOHnXYAFgOfq2feVOBSkj/EdgeOzcwLYP/MdCWwpNZ2nwF6AHsDy4A/AYcCHYAZwPcyy58LdE7nXQtUZeb9BvhhZroN8DzwXaA9SQAtBE5M518JPAl0A/YB5mVra+w1AO4FbgA+BpQBs4FvpPP2B4alde4FzAKubWBblbX3nV0GmAhsAUam/eoIXJy+dr3S/dwATK2n9h1+3YE+6e9xatrXgcDyTG0/SLdVlvb1j8Dlmf1uBX6cbrdjPX0dDBxFcjaqD/AycHGtf0cPAF1Jwn85MDyd92XgTZI/QJS+9vs29vv30IzvD8UuwIOHhob0DW8tsCozfD2ddzMwGehVx3qFBNiZmem7gesz0xcB99ZTU9d0+3um07/hgwF2JPC3Wut8B7gpHV9Y8yaYTo+t/cZax2tQ86bdA9gEdMzMHw08Uc+6I4E/17Wtul6XOvY3EZhVa/7LwPGZ6Z4kIdeujv3v8OvO+wHWPzN/EnBjOr4A+Hxm3onA4sx+NwO7N9TXOuq9GLin1r+j7B9GdwET0vFHgH+tYxsN/v49NN/g88KWByOj7mtg44HLgdmSVgJXR8Svm7DddzLjG+qY7gSQnp77Eclf3HsB29JlugPv1bHdfUlOVa3KtLUlOeoCKAfeyMx7vQk17wvsBizNnMlqU7M9SWXAz4FPkxwxtgFWNmH7dXmj1vS+wD2StmXaqknC9c0CtlfQ617P/l8nORKD5HV8vda88sz08ojY2FAhkv4R+CkwBNiD5Ejs+VqLvZ0ZX5+pbx+SEK2tsd+/NRNfA7Pcioi3I+LrEVEOfAP4L2XuPGxGXwFOBT4H7ElyZADJaSNI/krPegNYFBFdM0PniPh8On8pyZtfjd5NqOUNkiOw7pltd4mIg9P5V6T1DIqILsBXM3XWVes6kjfupENJWO9Va5m6+ndSrf7tHhGFhNeOqP1avZWOv0USFnXNgw/XXddXb1wP/BXol75e/84HX6+GvAHsV097Q79/ayYOMMstSV+W1CudXEnyBlWdTr9Dcu2hOXQmCY2/k7zZ/0et+bX3NRtYnd5A0FFSW0kDJNXcrHEX8B1JH0/rv6jQQiJiKfAocLWkLukNFftJ+kym1rXAKkl7A99upNZXgN0lnSxpN+AykmtGDflv4EeS9gWQtJekUwvtww74/5L2UHLjzDnAnWn7VOCydP/dSa453drAdt4BPiFpz0xbZ2A1sFZSf+Cfm1DXr4BvSRqsxP7pa9LY79+aiQPM8qDmrrma4Z60/XDgWUlrgftIrkcsSudNBKakd5eN+oj7v5nk9NSbwF9IbhzIuhE4KN3XvRFRDYwAKoBFwAqSN7uaN87vp9tbRBJGtzSxnrNIbg74C0lw/w/JdaiabR9GcmrzQeC3tda9guRNf5Wkb0XEe8AFaX1vkhyRNfZZqZ+RvN6PSlpD8noc2cQ+NMUfgNeAx4GfRMSjafsPgTnAXOBFkptB6v08XkT8lST0Fqb9Lwe+RXKEvQb4Je+HY6MiYhrJqeXb0/XvBboV8Pu3ZqL0AqOZ2U5FUh+SANgtcvjZM2t5PgIzM7NccoCZmVku+RSimZnlko/AzMwsl/xB5lbUtWvX2H//lviY0s5p3bp1fOxjHyt2Ga3G/S1t7m/xPP/88ysiovbnEx1gralHjx7MmTOn2GW0mpkzZ1JZWVnsMlqN+1va3N/ikVTn02p8CtHMzHLJAWZmZrnkADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxxgZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55AAzM7NccoCZmVkuOcDMzCyXHGBmZpZLDjAzM8slB5iZmeWSA8zMzHLJAWZmZrnkADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxQRxa5hl9G77/7RZtTPil1Gq/nmwK1c/WK7YpfRatzf0rar93fxlScXrRZJz0fEkNrtPgIzM7OC9enTh4EDB1JRUcGQIUmmvPvuuwwbNox+/foxbNgwVq5cCcBtt91GRUXF9qFNmzZUVVUBMHz4cA455BAOPvhgzj//fKqrq5tcSy4DTNInJFWlw9uS3sxMR/pznqT7JXVN1+kjaUNmuSpJZ6XzFkt6MdN+dLr8vFr7nSjpW5L+M13uL7W2eXoxXg8zs9b0xBNPUFVVxZw5cwC48sorOf7443n11Vc5/vjjufLKKwE488wzqaqqoqqqiltuuYU+ffpQUVEBwF133cULL7zAvHnzWL58OdOmTWtyHbk8Ho6IvwMVkIQKsDYifpJOr42ImnlTgHHAj9JVF9TMq8NnI2JFzYSkPg3sf1xmmQca2KaZWcmbPn06M2fOBGDMmDFUVlby4x//+APLTJ06ldGjR2+f7tKlCwBbt25l8+bNSGryfnN5BNYETwN7F7sIM7NSIYkTTjiBwYMHM3nyZADeeecdevbsCUDPnj1ZtmzZh9a78847PxBgACeeeCJlZWV07tyZ009v+gmsXB6BFUJSW+B44MZM836SqjLTF0XEk+n4E5KqgU0RcWQ9y/8D8JMm1jEWGAvw8U/sRZemrGxmtpN56qmnKC8vZ9myZQwbNoz+/fs3us6zzz7LHnvswYABAz7Q/sgjj7Bx40bOPPNMZsyYwbBhw5pUSykegXVMQ+fvQDfgscy8BRFRkRmezMz7bNp2ZH3LA//d1GIiYnJEDImIIZ26OL7MLN/Ky8sBKCsr47TTTmP27Nn06NGDpUuXArB06VLKyso+sM4dd9zxoaOvGrvvvjunnHIK06dPb3ItpRhgG9Kw2RdoT3INzMzMPqJ169axZs2a7eOPPvooAwYM4JRTTmHKlCkATJkyhVNPPXX7Otu2bWPatGmcccYZ29vWrl27PfC2bt3KQw89VNCRXG0lewoxIt6T9C/AdEnXF7seM7O8e+eddzjttNOAJHi+8pWvMHz4cA4//HBGjRrFjTfeSO/evT9wR+GsWbPo1asXffv23d62bt06TjnlFDZt2kR1dTVDhw7l/PPPb3I9uf8gcz13IXbKzL8fuAt4EngZmJ9Z/dcR8XNJi4EhddyF+EBEDMi01d7Xh5ZpyAEHHBDz589vfMESMXPmTCorK4tdRqtxf0ub+1s89X2QOfdHYBExsdZ0p1rTIzKTHevZRp862hYDA2q11d7Xh5YxM7PWUYrXwMzMbBfgADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxxgZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55AAzM7NccoCZmVkuOcDMzCyXHGBmZpZLDjAzM8slB5iZmeWSA8zMzHLJAWZmZrnkADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmudSu2AXsSjZsqabPhAeLXUar+ebArZzt/pYs97f4Fl95MgDV1dUMGTKEvffemwceeACAX/ziF1x33XW0a9eOk08+mUmTJvHYY48xYcIENm/eTPv27bnqqqsYOnQoAJdeeik333wzK1euZO3atUXrU1O0WIBJugZ4PSKuTacfAd6IiPPS6auBN4E+wFAggI3AqIhYJGkxsAbYBrwDnBURb2faAdoCvwUuj4hNktoA1zayvQBWptt7Pa2lGngxU/7IiFiczvsZcDqwT0Rsk9QH+F+gd0Rsy/S3ChgbEbM/8otnZtYEP/vZzzjwwANZvXo1AE888QTTp09n7ty5dOjQgWXLlgHQvXt37r//fsrLy5k3bx4nnngib775JgAjRozgwgsvpF+/fkXrR1O15CnEPwJHA6TB0h04ODP/aKAzUA4MioiBwGnAqswyn42IQ4A5wL/Xah8IHAH0BSan7f9UwPYGATOByzLtGyKiIjMsztR9GvAGcBxAOu8N4NM1K0vqD3R2eJlZa1uyZAkPPvgg55133va266+/ngkTJtChQwcAysrKADj00EMpLy8H4OCDD2bjxo1s2rQJgKOOOoqePXu2cvUfTUsG2FOkAUYSXPOANZI+LqkDcCCwAVhacyQTEUsiYmUd25oF7F+7MSLWAucDIyV1A3oWuL2ngb0L6MNn07qvB0Zn2qcCZ2Smz0jbzMxa1cUXX8ykSZNo0+b9t/NXXnmFJ598kiOPPJLPfOYzPPfccx9a7+677+bQQw/dHnJ51GIBFhFvAVsl9SYJsqeBZ4FPAUOAucDtwAhJVZKulnRoPZv7Ah88xZfdz2pgEdAPuKvA7Q0H7s1Md0zXqZJ0T6Z9NEkw3QN8QdJuaftdJKFZcwr2n4A76tqRpLGS5kiaszY9vDczaw4PPPAAZWVlDB48+APtW7duZeXKlTzzzDNcddVVjBo1iojYPv+ll17ikksu4YYbbmjtkptVS9/EUXMUdjTwU5KjnqOB94A/RsQSSQeQXLMaCjwu6csR8Xi6/hPp9am5fPCUX22C5IirgO31AJZRxynED2xQag98Hvh/EbFG0rPACcCD6bW4l4DjJb0DbImIeXUVFhGTSU9x9u67f9S1jJnZjnjqqae47777eOihh9i4cSOrV6/mq1/9Kr169eKLX/wikjjiiCNo06YNK1asYK+99mLJkiWcdtpp3Hzzzey3337F7sJH0tK30ddcBxtIciruGZIjsKNJwo2I2BQRv4uIbwP/AYzMrP/Z9JrUWRGxijpI6kxyI8grhWwP2Bd4CfhBI7UPB/YEXkxvADmWuk8j+vShmRXFFVdcwZIlS1i8eDF33HEHQ4cO5dZbb2XkyJHMmDEDSE4nbt68me7du7Nq1SpOPvlkrrjiCo455pgiV//RNTnA0mtYgwpc/CmS03/vRkR1RLwLdCUJsaclHSapPN1uG2AQ8HoTaukE/Bdwb0SsLGR7EbEBuBg4K71uVp/RwHkR0Sci+gCfBE6QtEc6/26SI7R6Tx+amRXDueeey8KFCxkwYABnnHEGU6ZMQRLXXXcdr732GpdffjkVFRVUVFRsv0Nx/Pjx9OrVi/Xr19OrVy9+85vfFLcTBSjoFKKkmcAp6fJVwHJJf4iIf2tk1RdJ7j68vVZbp4hYIWkI8Mv0pg6A2cB1BZT0hCSRBPA9wOVpe1kh24uIpZKmAuMy626XhtSJwDcy66yT9L/ACODOiFgl6RmgR0QsKqBmM7MWU1lZSWVlJQDt27fn1ltv/dAyl112GZddVvfVmEmTJjFp0qTt0zNnzmyJMpuVshf26l1I+nNEHCrpPJLPQ31P0tz0lnQr0AEHHBDz588vdhmtZubMmdv/Q+0K3N/S5v4Wj6TnI2JI7fZCTyG2k9QTGAU80KyVmZmZ7YBCA+wHwCPAgoh4TlJf4NWWK8vMzKxhBV0Di4hpwLTM9ELgSy1VlJmZWWMKOgKT9I+SHpc0L50eJKmhz2WZmZm1qEJPIf4S+A6wBSAi5vLBRymZmZm1qkIDbI86HlS7tbmLMTMzK1ShAbZC0n4kX0WCpNOBpS1WlZmZWSMKfRbiOJLn+fWX9CbJw3PPbLGqzMzMGtFogKWPZBoSEZ+T9DGgTUSsaWw9MzOzltToKcT0u7UuTMfXObzMzGxnUOg1sMckfUvSPpK61QwtWpmZmVkDCr0Gdm76c1ymLYC+zVuOmZlZYQp9EscnW7oQMzOzpij061TOqqs9Im5u3nLMzMwKU+gpxMMz47sDxwN/AhxgZmZWFIWeQrwoOy1pT+CWFqnIzMysAIXehVjbeqBfcxZiZmbWFIVeA7uf9DFSJKF3EJmvVzEzM2tthV4D+0lmfCvwekQsaYF6zMzMClLoKcTPR8Qf0uGpiFgi6cctWpmZmVkDCg2wYXW0ndSchZiZmTVFg6cQJf0zcAHQV9LczKzOwFMtWZiZmVlDGrsGdjvwO+AKYEKmfU1EvNtiVZmZmTWiwQCLiPeA94DRAJLKSD7I3ElSp4j4W8uXaGZm9mEFXQOTNELSqyRfZPkHYDHJkZmZmVlRFHoTxw+Bo4BX0gf7Ho+vgZmZWREVGmBbIuLvQBtJbSLiCaCiBesyMzNrUKEfZF4lqRPwJHCbpGUkH2g2MzMrikKPwE4lef7hxcDDwAJgREsVZWZm1phCn0a/TtK+QL+ImCJpD6Bty5ZmZmZWv0LvQvw68D/ADWnT3sC9LVWUmZlZYwo9hTgOOAZYDRARrwJlLVWUmZlZYwoNsE0RsblmQlI73v96FTMzs1ZX6F2If5D070BHScNIno94f8uVVZo2bKmmz4QHi11Gq/nmwK2c7f6WrFLs7+IrTy52CdYEhR6BTQCWAy8C3wAeAi5rqaLMzIph48aNHHHEERxyyCEcfPDBfO973wOgqqqKCy64gIqKCoYMGcLs2bMB2LJlC2PGjGHgwIEceOCBXHHFFdu3demll7LPPvvQqVOnovRlV9DY0+h7R8TfImIb8Mt02GVI+gTweDr5D0A1SZAD3AOMStu2Ad+IiGdbvUgzazYdOnRgxowZdOrUiS1btnDsscdy0kkn8d3vfpcxY8ZwySWX8NBDDzF+/HhmzpzJtGnT2LRpEy+++CLr16/noIMOYvTo0fTp04cRI0Zw4YUX0q9fv2J3q2Q1dgrxXuAwAEl3R8SXWr6knUf69JEKAEkTgbUR8RNJnwJ+ChwWEZskdQfaF69SM2sOkrYfMW3ZsoUtW7YgCUmsW7cOgPfee4/y8vLty69bt46tW7eyYcMG2rdvT5cuXQA46qijitOJXUhjAabMeN+WLCRnegIrImITQESsKHI9ZtZMqqurGTx4MK+99hrjxo3jyCOP5Nprr6WyspKbbrqJbdu28cc//hGA008/nenTp9OzZ0/Wr1/PNddcQ7du3Yrcg11HY9fAop7xXd2jwD6SXpH0X5I+U9+CksZKmiNpztrVq1uxRDPbEW3btqWqqoolS5Ywe/Zs5s2bx/XXX88FF1zAG2+8wTXXXMPXvvY1AGbPnk3btm156623WLRoEVdffTULFy4scg92HY0F2CGSVktaAwxKx1dLWiNpl303joi1wGBgLMk1sTslnV3PspMjYkhEDOmUnlows51f165dqays5OGHH2bKlCkcd9xxAHz5y1/efhPH7bffzvDhw9ltt90oKyvjmGOOYc6cOcUse5fSYIBFRNuI6BIRnSOiXTpeM71LvxtHRHVEzIyI7wEXArvU9UGzUrR8+XJWrVoFwIYNG/j9739P//79KS8v54UXXgBgxowZ22/M6N27NzNmzCAiWLduHc888wz9+/cvWv27mkI/B2YZkg4AtqVPJIHkRo/Xi1iSmTWDpUuXMmbMGKqrq9m2bRujRo3iC1/4Al27duXcc8/lpptuYvfdd2fy5MkAjBs3jnPOOYcBAwYQEZxzzjkMGjQIgPHjx3P77bezfv16evXqxXnnncfEiROL2LvSowhf2ipErbsQBwO/ALqSfK3Ma8DYxm7mOOCAA2L+/PktXuvOYubMmVRWVha7jFbj/pY297d4JD0fEUNqt/sIrEARMTEz/jxwdPGqMTOzQp/EYWZmtlNxgJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxxgZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55AAzM7NccoCZmVkuOcDMzCyXHGBmZpZLDjAzM8slB5iZmeWSA8zMzHLJAWZmZrnkADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAtbJzzz2XsrIyBgwYUOxSzMxyrSQDTFK1pCpJ8yRNk7RHHe33S+qatveRtCGdVzOclc7rJOkGSQskvSRplqQjM8u9LenNzHT7hmo7++yzefjhh1v+RTAzK3Htil1AC9kQERUAkm4Dzgd+Wqt9CjAO+FG6zoKaebX8ClgE9IuIbZL6AgdmtjMRWBsRPymksOOOO47FixfvcMfMzCxRqgGW9SQwqI72p+tp307SfsCRwJkRsQ0gIhYCC5u7SDMza5qSDjBJ7YCTgIdrtbcFjgduzDTvJ6kqM30R8HGgKiKqP0INY4GxAD169NjRzZiZWS0leQ0M6JiG0Rzgb7wfVDXtfwe6AY9l1lkQERWZ4cnmKCQiJkfEkIgYsueeezbHJs3MjNINsA2ZILooIjZn24F9gfYk18Aa8hJwiKRSfZ3MzHJrl3xjjoj3gH8BviVptwaWW0ByFPd9SQKQ1E/SqTu679GjR/OpT32K+fPn06tXL2688cbGVzIzsw8p6WtgDYmIP0t6ATiD5EaP2tfAfh0RPwfOA64GXpO0nuT047d3dL9Tp079CFWbmVmNkgywiOhUSHtEjMhMdqxnndXA1xvY18QdKNHMzD6iXfIUopmZ5Z8DzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxxgZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55AAzM7NccoCZmVkuOcDMzCyXHGBmZpZLDjAzM8slB5iZmeWSA8zMzHLJAWZmZrnkADMzs1xygJmZWS45wMzMLJccYGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLjnAzMwslxxgZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55AAzM7NccoCZmVkuOcDMzCyXFBHFrmGXIWkNML/YdbSi7sCKYhfRitzf0ub+Fs++EbFX7cZ2xahkFzY/IoYUu4jWImmO+1u63N/Slof++hSimZnlkgPMzMxyyQHWuiYXu4BW5v6WNve3tO30/fVNHGZmlks+AjMzs1xygJmZWS45wFqBpOGS5kt6TdKEYtezoyT9WtIySfMybd0kPSbp1fTnxzPzvpP2eb6kEzPtgyW9mM77uSS1dl8KIWkfSU9IelnSS5L+NW0vyT5L2l3SbEkvpP39ftpekv0FkNRW0p8lPZBOl2xfASQtTmutkjQnbctvnyPCQwsOQFtgAdAXaA+8ABxU7Lp2sC/HAYcB8zJtk4AJ6fgE4Mfp+EFpXzsAn0xfg7bpvNnApwABvwNOKnbf6ulvT+CwdLwz8Erar5Lsc1pbp3R8N+BZ4KhS7W9a578BtwMPlPq/57TWxUD3Wm257bOPwFreEcBrEbEwIjYDdwCnFrmmHRIRs4B3azWfCkxJx6cAIzPtd0TEpohYBLwGHCGpJ9AlIp6O5H/CzZl1dioRsTQi/pSOrwFeBvamRPscibXp5G7pEJRofyX1Ak4GfpVpLsm+NiK3fXaAtby9gTcy00vStlLRIyKWQvKGD5Sl7fX1e+90vHb7Tk1SH+BQkqOSku1zekqtClgGPBYRpdzfa4HxwLZMW6n2tUYAj0p6XtLYtC23ffajpFpeXeeGd4XPLtTX79y9HpI6AXcDF0fE6gZO9+e+zxFRDVRI6grcI2lAA4vntr+SvgAsi4jnJVUWskodbbnoay3HRMRbksqAxyT9tYFld/o++wis5S0B9slM9wLeKlItLeGd9JQC6c9laXt9/V6Sjtdu3ylJ2o0kvG6LiN+mzSXdZ4CIWAXMBIZTmv09BjhF0mKS0/pDJd1KafZ1u4h4K/25DLiH5BJHbvvsAGt5zwH9JH1SUnvgDOC+ItfUnO4DxqTjY4DpmfYzJHWQ9EmgHzA7PUWxRtJR6Z1LZ2XW2amk9d0IvBwRP83MKsk+S9orPfJCUkfgc8BfKcH+RsR3IqJXRPQh+T85IyK+Sgn2tYakj0nqXDMOnADMI899LsadI7vaAHye5A62BcClxa7nI/RjKrAU2ELyV9jXgE8AjwOvpj+7ZZa/NO3zfDJ3KQFDSP7jLACuI30izM42AMeSnBqZC1Slw+dLtc/AIODPaX/nAd9N20uyv5laK3n/LsSS7SvJndAvpMNLNe9Fee6zHyVlZma55FOIZmaWSw4wMzPLJQeYmZnlkgPMzMxyyQFmZma55CdxmOWcpGrgxUzTyIhYXKRyzFqNb6M3yzlJayOiUyvur11EbG2t/ZnVx6cQzUqcpJ6SZqXfATVP0qfT9uGS/qTk+78eT9u6SbpX0lxJz0galLZPlDRZ0qPAzelDf6+S9I2xHbMAAAFrSURBVFy67DeK2EXbRfkUoln+dUyfIA+wKCJOqzX/K8AjEfEjSW2BPSTtBfwSOC4iFknqli77feDPETFS0lCSr8qoSOcNBo6NiA3pk8zfi4jDJXUAnpL0aCRfu2HWKhxgZvm3ISIqGpj/HPDr9MHE90ZEVfoE9lk1gRMRNd/zdizwpbRthqRPSNoznXdfRGxIx08ABkk6PZ3ek+RZeQ4wazUOMLMSFxGzJB1H8uWNt0i6ClhF3V+B0dBXZayrtdxFEfFIsxZr1gS+BmZW4iTtS/LdV78kebr+YcDTwGfSp4yTOYU4CzgzbasEVkTE6jo2+wjwz+lRHZL+MX3CuVmr8RGYWemrBL4taQuwFjgrIpan17F+K6kNyXdADQMmAjdJmgus5/2v2ajtV0Af4E/pV2osp0hfK2+7Lt9Gb2ZmueRTiGZmlksOMDMzyyUHmJmZ5ZIDzMzMcskBZmZmueQAMzOzXHKAmZlZLv0foUJK2/EFsL8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "ax = xgboost.plot_importance(bst)\n",
    "ax.grid(False, axis=\"y\")\n",
    "ax.set_title('Estimated feature importance')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Test model performance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we test our XGBoost emulator performance using the testing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_hat = dask_xgboost.predict(client, bst, X_test).persist()\n",
    "y_test, y_hat = dask.compute(y_test, y_hat) # get the predicted data and true data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU0AAAE9CAYAAACP0jAFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nOy9eXhT5533/bm9SF4kg7EVY3aDA45jIIAD8SRMS5oFOkygaZ6EbGWaJszThDaetgN0UubqlPadkM6ToS1t5yVp3qHNQmgWkvI0tDChHdJxSVjCZoQRi4GgGFvYaLEtWdb9/iGdgyzLZqltjPl9rsuXpXPuc3SUwJff/VuV1hpBEATh4ki50g8gCIJwNSGiKQiCcAmIaAqCIFwCIpqCIAiXgIimIAjCJSCiKQiCcAmkXekH+EvIz8/XY8aMudKPIQjCAGPnzp0NWmtHsnNXtWiOGTOGHTt2XOnHEARhgKGUqu3qnGzPBUEQLgERTUEQhEtARFMQBOESENEUBEG4BEQ0BUEQLgERTUEQhEtARFMQBOESENEUBEG4BEQ0BUEQLgERTUEQrmrqfcE+/TwRTUEQrlrqfUGWb9jfp8IpoikIwlWLw25lxfwyHHZrn32miKYgCFc1fSmYIKIpCIJwSYhoCoIwoOlpf6eIpiAIAxYjUOR0e3vsniKagiAMWDz+IJV3XM/KTc4eszhFNAVBGJA43V4WvvQhjYFQj95Xaa179IZ9SXl5uZZxF4IgJKPeF8TjD5Jni0bXLyXKrpTaqbUuT3ZOLE1BEAYE8dtvw5cJsOSNPXj8PRcMEtEUBOGqJ7EyyEh6z7NZCYUjrNhYLT5NQRAEg2SVQQ67FYfdyvK5pT36WSKagiAMCJL5LOt9UZ+mJa3npE5EUxCEAUm8X/O5+yb3WLmliKYgCP2arnyRF/JROuxWKu+4HkACQYIgXBvU+4IseWMPVa6GTscv1BLO6faycpOTJ365g0de3N5jVUEimoIg9GsCwTBfe203Tre3U3S8qy13vS/Iqi2HuX/aSLItqfzjXRMoKczpkecR0RQEoV9S7wvisFv56cPT+MmDU8izWTulFQEdxNTA2Jqv33mSQKidH/7+kFiagiAMXBK332u2HQXoZF063V4eevHPVK7b3Uk4SwpzeO6+yTx770R+8uCUHrM003rkLoIgCD3Mivll5uvGQKhDSaRBns3KuHwby+eWJt2qe/xR8bVlpFFcYO+RCLpYmoIg9CsMK9NV52PJG3t4v7qOA6e9LHtzb1KL0paRZoppYinlio3VHK33c+vYPEk5EgRhYGL4I1dvddHUHOK1j05w47Acls4uSbrWyMFMVkq5fG4pS+4u4YebazpF4C8X2Z4LgtDvMGrGtYYTjc0svauE4gJ7t5U9RkTdwIigr5hfRk5mOhXF+T3ybGJpCoLQL7FlpLF0dglL7yrhh78/hKvOx9LZJR222UYeZ/y23HhviKjHH+SdPad7rGGHWJqCIPQ7HHYrS2eXsGJjNaFwhBFDsnh+cw22jLROJZH+1nC391q5ydlJbP8SxNIUBKFfYmzRw5EID5aPRClYNHNsB/Hz+IMc9wTYcews0NHHmXivnkJEUxCEfkGy7XNbe4SaOj//z6aDHHSfY/VWV4d1JYU5fO+eMr77mwNmkntie7iebNYBIpqCIPQD6n3BpOlE4XZNsK2dL1eMoWToIB6ZMbrTdeVFQ1j72HTybNZO/k3AjKz3FCKagiBccVx1PmrqfGY3IjNtKMfKD+ZP5E9HPDQ2h3hmwz5TXI0gUOW63d1uvy+mucelIKIpCMIVxen2sraqlh8viNaXJ/bBzMlMx1nn41h9gC9OHc6qBVPMruxLZ5eYaUjG+8St+IWae1wqIpqCIFwRDGtx5SYnCytGk5ttMcVy1ngHDrs1mmb01h6us1koGJTB2qpaXHU+8x5Gfblxv1VbDie1KHvSpykpR4Ig9DmGNVl5x/UEgmFWb3VhSUth0cyxvFJ1nB9vdeE+18KfjngIBNtpDTXTHoEF00cmTVJfvmE/K+aXdWlRJgaI/hLE0hQEoce5mK7qxrTI9NQUls8tZensEp7fXMPqP7iYf9MwXvnwBDeNHES7hlAE8u1Wfr3zVKdyyPjtd1eCKT5NQRD6LRcrUobAWdJSyLNZKSnM4bFbi8i2pnG0PkBORhrrd57Cmhpdn5Gewop7ypJamt1ZkUYtu1iagiD0Sy4l8OKwW1k0cywQFdt39pzm2S9Moi2iafCHKMyxYremA1B7toU3dp3q1Ez4QuLsdHu79HVeDiKagiD0OBdr1TndXr6+bjeV63az49hZKu+4Hm9LGyfPBhg6KINzLWEamtvM9beOy+sggE63t0MKUjz1vqA5J2hhxegeszQlECQIwhWjpDCHX31lBntONPH067sZmZvJ8bPNWNNSQGvc3iApQATIz0pnxtg8Hq4YA0QF85vrP+a4J4CrzsfqrS4zHclIlg+E2jl1tpklb+7lrSdvlSbEgiBcfSQ2CgbYWlPPk58ZR3pKCjoCI3OzuHH4ICAqmAANzW18Ze1H7Dh2lifWfsSyN/dS7fYRCEU48Mk5qt3eDqN629ojHDjVREMgxMnGFtb80dUjzy+iKQjCJWEI3YV8hMnOxweJjIqeZW/tY97kYfxyey21jQEKcqykpyre2n260/XNbRGe3XSQ/Z+cI8+Wjo4d9wfDqLh1DruV780ro2zEINIVZKSl8MVpIy/3K3eg10RTKZWhlPpQKbVHKXVAKfUvseNDlFKblVKHY79z4675tlLKpZQ6pJS6u7eeTRCEy8MQPafb222EPDGCbohkfCTbYbdy/7SRVJ8+xyeNzXxpxmha2qJb8kdnjKE4P7PL57Cmp/LhsUZTwAKhMK8+cUuH4WklhTms+dLN/ObrM9nw1K09NlhNaa0vvOpybqyUArK11n6lVDrwAfA0cC9wVmv9rFJqGZCrtV6qlCoFXgOmA8OALcB4rXV7V59RXl6ud+zY0SvPLwhCcgzxc7q93QqRsc7wLwIsnlXM2qraDtH1H2w8wH9W1TKhwEbNpz5CEfh8WQG/3V+X9L6VtxfzpyMeAL5x53jed9bxzNwbe/Q7KqV2aq3Lk53rNUtTR/HH3qbHfjQwD1gbO74WmB97PQ9Yp7UOaq2PAS6iAioIQj/CEMILpfEkBl1C4Qirt7pMS7PeF6TK1cCbuz9h/HU25k8eTiTmwOxKMDPTUpg9sZCfPTKNnz0yjYri/B4XzAvRq9FzpVQqsBMoBn6qtd6ulCrQWrsBtNZupdR1seXDgT/HXX4qdkwQhH5GfC5mvEWZKJTGsVULpuDxR+vM82xWc3uvFIzKzWLuxEJ+9P5huuvBbrPAG1/tuW325dKrgSCtdbvW+iZgBDBdKVXWzXKV5Fgn34FSapFSaodSakd9fX1PPaogCJdI/ATIZD5OI9DjdHtx2KMVP0tnl+DxR0frVn/qJRSO0NLWzv/Z7MQX7NITB8AP75t6xQUT+ih6rrVuAv4AzAbqlFKFALHfZ2LLTgHx4a0RQKfwmdZ6jda6XGtd7nA4evW5BUHoHsPiLCnM6VQF5PEHCYUjrNhYbXZVX7GxmpWbnNw8OpeWYDufHe+g1uOnpfsxPwAMzkrvxW9y8fRm9NyhlBoce50J3AE4gXeBhbFlC4F3Yq/fBRYopaxKqSLgeuDD3no+QRB6hsStuREpX7XlMItnFRMKR1i5yYnHH8SSlsLs0qG8+D9H0cArH54g2Vw0BeRlnxfJAruV4gJ733yhC9CblmYhsFUptRf4CNistd4IPAvcqZQ6DNwZe4/W+gCwHqgGNgFPdRc5FwShb0lWpghQ5WowSxmNbbrHH2TF/DJysy3mQDRje/7u3tNE2qO+t3p/KOln3V6SjydwvnzS2xrskLh+JenN6PlerfUUrfUkrXWZ1vp7seMerfXntNbXx36fjbvmB1rrcVrrCVrr93rr2QRBuDQSfZbG+ypXA0+9tov9n5zDVedj1ZbDLKwYzaoth03fZU2dj3/fUmM22rhn0jDa2yNdfta9U4Zx39RRHY6pflSHI7XngiB0i7HVjk8VMt6XFObwg3kTeXl7LbnZFtOvWVxgx2G3snxuKY2BEKu3uvjnd/ZzuM5HS1uEYIJmKsBhs3DGH+KDww0cqvOTrqBNw4jBGRQOyuzRMbx/Cb2W3N4XSHK7IPQN5qCzOL9lfPR81ngHW2vqqbzjelPcXHU+1lbVMm/yMD5pbOa28Q427XOz6v3ONeD52RZG5GaSZ0vnv5wNWFKijYchmsz+cMWYHh1ZcSG6S24XS1MQhItiyRt7zHk88cGfWeMd/NvmQ3zvnjJWbnISCkcIhSMc9wT4+5ljWfzaLto1DPuf45w+19rhntmWVDLTUmhqCREMt5NtHczD00fy9u5PSEmBghxrnwvmhRDRFAThghhBGFedjzXbjrJo5ljWVtWysGI0K3/vZFRuFuVFQygvGtIpYLPAE+CVD092Esy0FHjitiL+UFNPhj/ID++bTFNzG8ve3ktzW4RsayrP3jupXwkmiGgKghCjq8qeKlcDa6tqWTRzLM9vrgHg+c015oyfcfk2vnHneACz6geg2JHNr3edoqm5c05RYY4VfzDMH2rqOXzGTyDUzvajHjYfrOM/Hp5GU3MbRY7sfpHMnoiIpiAIHaZDrtzk5Ln7JptNOSpf/5hVD9xEU3Mbxz0BvnXnBP5t8yHz2raIZvXWqJ8yFI6wYn4ZH9TU8/33nJ0+Z+TgDP56vIPjnmYWzyqmuMCOq87Hd39zgPf2uzl0JsAJTzMPzBjV6dr+goimIAhmZU8iJYU5rH0s2jfn6+t28717ypgzqZDJowaTZ7Pi8Qc54Qnw04em0tTcxjPv7ANg9daaTvfKSIXT51q5rdhBTd0xM8IO8KMFUwD45vqPub20oLe+Zo/Qf5KfBEG4ohg9Lg0r04iY59ms5NmsjMnLpsiRDUBjIMQjL1QBcENsC/3y9loKczIAaGrtmFNUmGMlPT2NLEt0tORxTwCPP1o59OTLO81GHt/5m9J+58NMRERTEK5hkrV2i++BWeVqYPmG/bjqfAD88zv7eW+vmyd+tYNDZwJ8UFPP4lnFrNl2lKbmEEfqAyx7c2+H7jvWVEVzKMx35tzAmkfLmTOpkB8vmEJJYQ4ef5AjDX4WzRyLxx/kqdd2dZo22d8Q0RSEa5Tu5pO76nzRSY6/O8S8ycNYs+0obe0RDn3qZc22Iyz/fCkTCmwMz81izbaj5Gal420N0xqO4PzU26E9WTiiybdZeWPXKdZsO4rT7WVtVS1Ot5c8m5Ubhuacryu/CtLGJbldEK4xEptrAEmj5eWjBvPzPx6hpDCH5XNLAXh63W7OeIOMzsviwZtH8eymg6Qq1WHMLmBW8zw8fSR7Tp4j05KKJS2F5XNLKSnMMWeRG35U4/Mv1A2+r7gindsFQehfGN2HEq1Lo4YcoqL1tdd2M2/yMLbW1LP6oamm0DUGQpwNhHjqs+M40djMS386SktbeyfBVIA9K53PleRz+lwrmZZUVswvY/ncUrNaKL6VXLxg9wfBvBAimoIwgIkfbLZ8w36ATh3XF1aMpvL1j6lyNZBnszJhaHSrfPBTLyc8zazc5DS31N+6cwI3Dh/Ek389jhNnA7SGO+9UJ48YxNK7Stj/iY/7p43EkhaVmRUbq81uSNB5HMbVgoimIAxQ4q3KxCmQRoCn3hekojifVQ/cxNqqWiA6/Gz9zpMUDsrk2U0HOeNtpTEQonzUYL7/22oe+cV2/m3zoU6Ng1OA78wpIduaxu2lBax9bLoZbTcIhbvubnS1IHmagjBASZzjYyStu+p8fPXVnfz8oWnmuYrifDMYs3qry/RhLntrHzVnzrHoVzsIBNuxpCnaNbTHWZg5VoU/qMnJTOP3B+swwiSNgRBrth01hdKIsl/tiKUpCAOYxC2wxx/k+c01+Fuj/b0TfZwef5CDn3ppDESbAz9YPpL0lFS+cmsRKSmKUJLteCisyc1K51+/MIl/uCNaTrnj2Fm+9tpubh2bh1LRrfmabUdZOrvkqt2WG0j0XBCuEYwUH4imFFUU5wNR4TRqxu+fNpI1247QEgpT5w0SAcLtETLTUjsFfAy+M6eErbF8zdVbXRxwn+N6hx1fSxt1viDF19nMOvWrRTC7i56LaArCNYAxGRIwK36M40aH9UAwzHFPgKH2DA6d8fPsFybiPtfCe/vdHGsImP0tDcYMyaD2bCuFgzJY/jelrN95kqWzS8zzKzZWm7XoV0NUPB7ppykI1yhGEMhht3L/tJEUObI75ESu2Fhtrl06u4R/31LDXTcU8P33nPzx8Bl+u78OiNaNJ7J0dinf/c0BTp9r5dlNB8lMS+1gTRopRleLdXmxiE9TEAYYxoAzw7qs9wWpcjXw9XW7TP+l0+1l5SYngWCY5XNLWT63lOICO96WEMNzs0gDUzABWuNGHFpSYOUXJlLkyKZwcCYAd5UWcOiMnx3HzprPYLSIG2iIaArCAMKoGX/oxT+z49hZQuEIHn+Q4gI7EwqiCeUef3TGz/3TRnLybDPH6gOs3OTklarjOOsCfPvtvXQ3hlxrGJWXxfIN+5lbNhQFTB01hJ89NJU5kwr76JteOUQ0BWEAYQwzG5dv4+XttebMcVedj9PnWmgMhKJJ5ndcT3nREEbmZfPy9loam9v4zd7TKOBcYgJmDKMJR5uG7Uc9fHyqieG5WUweMYjyoiEdBDO+W9JAQ3yagjDAKCnM4WePTAOiKURG38vReVkAHHSfTynKtqTyyIzRrN56mCMNzQBkptEpcT2eJ24bw8MVY3h9x0lzxEUycRyIggliaQrCgMRht5rbcFedj+Ub9nOkwU9TcxvF19miCewb9hMKR1j1XzWcazmfTtSVYBpiUeyw80rVcdzeIK9UHR+w4tgVYmkKwlVMsrk+Rj/KlZucLJo51mzrpnW0UfA/3DGe3GwLKzc5uXVsXtKxFAYpQFF+FqF2TYqC2rMt5GSm88CMUQRCYSrvmtAXX7NfIZamIFylGNU8Tre3Q1T8kRe38+239hIKRygusPPcfZN59ouTGJdvY/GsYrPG/P5pI7u8d3rsdwS4ZWweL3ypnBG5WZQNs1NeNIR6X5CTja1Je3EOdCS5XRCuYpJZmu/tdZuJ5kZSudPt5aEX/8xPH4zO8lm99TAH3D7GDMng+NnoaF2bReEPaUYMzsBmSaP2rJ+cLCt13iArvzDRnAsUnxg/ULfm0k9TEPohvWGl1fuCnQTToDAng+/+5gD/9PZeqt3R8RWGYAL89fjryEyP+kPz7Faev38qhTkZpCp4bcdJVm5ydphpPlAF80KIaArCFaC7UROJ65K9jr+H0+01k9gN8mxWs+lwvS9aJnm0wYfrjJ9ZJdclnSoxdWQumZZ0ltxdwvK5pZQXDeGFhTfz8ldm8MKXylk6u4RVWw73+xk+vY2IpiBcAeLbtiUjvsu6Ud2TKLJdjd016r+XvLGHynW78fiD3DNpGME2SE+BbYcbzLXZ1mj2ZeXtxew40cRPH5xqBomMayuK83HYrZQU5lB5x/Ws2nL4mvRlGohoCsIVojvBNLqsL6wYzaothwG6FNlVWw6bQvnkyztZsbEajz/Iopljow0zNlbzy+21RIjOLUtR52dFBoJRm3PC0Gi1UG62hVVbDrNo5lggGoGP7/4eP6biWkVSjgShnxFvQa6tqjU7rne31mG34nR7OdLgZ1RuFis2VtPga6U9oln98DTW/PcRDpz2EmyHs4HzVuK9U4YxY0we63eeZHBWutnzsqTw/IRII9C0fMP+a14wQaLngtCvuZQItdHmLc9m5f3qOpa+vQ+IToQ8fMbHh8ebgGj/yz8cPsOek+dY82g5a7Ydxd8axpKWQiAY5oWFNwOdLeGBHC1PRKLngnCVYlh5Fwq+GI06GgMhHHYrt5cWYI21c3vlw5PsP32uw3qn289Xbi2iojif5+6bzM8emcbiWcWcaGw2K4gSP/NaEcwLIaIpCFeAi41AO91eKtft5tFfbO/2Gledj90nz/LVV3fy3l43qzYfItjecY0iWmt+4/BB5NvS+dkfj1DlajD7bRYX2LlhaHRbLgGfrhGfpiD0AJeydXW6vSx86UPWPjY9aUfz+IT1VVsOm0POjFEViWsBvv9/q2kOaQZnpvGN9bvIsqZ3WFeYk8m37irh5e21sXtl8PTnJpgjLyBqSa5aMMUUUfFfJkcsTUH4C7nYnEuDksKcbgUzcexuSWEOeTZrp88w1nr8QRr80a5Fp8+10hKG3Mzz9tBtxUPMdm/xEyHLi4Z0+vx4kRTBTI4EggShB/hLgiSJ18ZbmvERa0MwjbWvbz8BgK+1rUPTjfxsCxnpKZxqilb7VN5eTOGgTFb+zskNhTnmGArjXtdSgOdikUCQIPQyf4lgdmWlJkteNwI0r28/wdK397H07X38cLMTSyz1Mh0YMSSL0XnZ3DtlGMX5WcyeWMi7e0/jsFtZPKvYtHATm30IF4eIpiD0IYnilOg7TCZi8dv1WeMd/PM7+zvcIxiGUGzD2Aa0tUd4ZMZoGpvbWP3wNPJsVgKhdo6c8fPvW2o6bP3zbOK7vFQkECQIfURXCeLxXYNWbTnMvMnDzGPx1uZ7e90si+VeTh4xqMO90+D8XB+t4wI+0fu88KVyXHU+igvs5pbcGHz23H2Te+PrDlhENAWhj7hQRNpht7KwYjSVr39MkSO7Q6Do717aTnu7NhttFNgzUGC+DwM2awopKYrv/E2pWc0DmEIdf8yY4WO8Fi4eCQQJwhUkWRDG6fZ26IO550QTS9/ex5ghmRw/2wJAlkXRHDr/d3f0kExG5GbxyIzRFDmyO/S9jO/kDgzYgWc9SXeBILE0BeEy6ImIs7FdX1gxukO+pCGY7+11808b9nGdzQLArcX5HP/wJACtcYKZZUnh2XsnAbB6q8sU3cWziikusJtiuXR2yTXTRLg3kUCQIFwC8S3bLjbi3FUlj8NuZd7kYVS+/nGnNe/tdbPsrb0E285POfvU22K+jgCWFJhQYCM9NfrXeM22oyyfW8rqh6ayeFYxX311Jx5/kKWzS8wmHPGCKVHzy0NEUxAukviWbYm+ya7Ex6j+SSacRpf1VQ/c1MF/6XR7WfbWXlrbwqSnpdIeiQBw8HTHe6yYN5EfLZjC+Ovs5GZbzONrth2lqbmNQLCdY/UBlr25t0OLN7iwf1Xoml4TTaXUSKXUVqXUQaXUAaXU07HjNyml/qyU+lgptUMpNT3umm8rpVxKqUNKqbt769kE4VIxtrKG0CQKZldWW3fVPwbFBXbTgjXEta09TLAd5k4sxNUQtTBPe0PmNePys7i9tIA8mxVLWgp5tmhgx4iWlxcN4Zdfns7grHSq3d6kQ9REMC+P3vRphoFvaq13KaXswE6l1GbgOeBftNbvKaU+H3v/WaVUKbAAuBEYBmxRSo3XWrd39QGC0BckSxWK9wdeyGrrSjAdditLZ5fg8UfTf/ytYY7W+1lydwmjcm04z/jJsqRis0CsShJLKiy5q4Qbhw8yU4csaSnm/eB8oMd4/6MHplBeNET6YfYQvWZpaq3dWutdsdc+4CAwnGiWhPGnaBBwOvZ6HrBOax3UWh8DXMB0BOEKkyiKXY2e6Iqutu5GruSKjdUsmjmWb9w5nqGDM/mXjfvJyYraMy98cNwUzLQU+NEDU7ltvIPK1z82OxQtnV3SZc240+3lnT3Rv2IimD1Dn/g0lVJjgCnAdqAS+KFS6iTwb8C3Y8uGAyfjLjsVOyYIV5xEUepKgBIFMtnQs/j7LJ1dQigc4fnNNazc5MTd1EJzmzYbBgPE2mKy7O4S5kwqpKQwh1UP3MSabUdxur1dtnAzzhmd30Uwe4ZeF02llA14E6jUWnuBrwL/oLUeCfwD8AtjaZLLOyWRKqUWxXyhO+rr63vrsQWhW7oSzHiBNCzJUDiSdG29L0hjIERbe4TG5hA1dT7ONrd1WvvtOSX8/KGpPP6ZceYxI1E9z2ZNOg7DqC4yuiQJPUev5mkqpdKJCuYrWuu3YocXAk/HXv8aeDH2+hQQ760ewfmtu4nWeg2wBqLJ7b3w2IJwUVwoz9GwJONzI43rlryxh0AwzKFPvQSC7UQ0OGwWmttCHe7x8PSR7DjRxI3DB3W6t1HRs2rL4aSlmbId7x16M3quiFqRB7XWz8edOg18Jvb6duBw7PW7wAKllFUpVQRcD3zYW88nCH8JXfk14/2LhrWXiCF435tXhtaasIaUFKjznxfMe6cMw25N5binuctcTmPLvbBidFJxFMHsHXrT0rwVeBTYp5T6OHbsn4AngB8ppdKAVmARgNb6gFJqPVBNNPL+lETOhf5KMkvOEMn4tKTENU6317Q836+uwxuMbt0zLSn4Ws9v40flZvHr//1XQDT6nliLblDlaqDy9Y8vmNYk9By9Jppa6w9I7qcEmNbFNT8AftBbzyQIvUm8SBpb90TBfOjFPzMqN4sHbx6Fq95nnmtv7+j3LByUSZ7NypI39vDcfZO77PK+tqq2U3K80LtIRZAgXIBkkWnDL5lsy2z0xDTOGb/zbFZyM9Op/uQcy97exydN58sim9vg82UFvPb4DFZ+YSIPzBh1wecyRDq+bl3ofUQ0BaEb4n2XieIZCkfM8sQqV4O5fsXGahZWjGbVlsNUuRp44P+twun24qrzcaShmZCGOWUFBMPnvU8FNgv/Mm8iFcX5pmAmy8FMJH4MhtA3iGgKQjfENwGOD/wYoyOWzi7BVefj7/7zI6pcDXj8QWrqfORmW1gxv4ym5jbOtYbZc6KJf99yyLzvvlPnuOuGQvP9EzPHAh2tWsNH2p0oSuONvkdawwlCAslSiRLn9TjdXr6+bjfjC+wsn1vKf/7dzeY2+VdfmWFeY/DL7bWcbmw2359rbWNUXhY3jRzM3LKhzJs6giVv7AE6lkFeKG1IUov6HrE0BSGORMst/r3Hf/51SWEOP14wheVzS1m15bDZdAPgg5p6s7PRrhNngWiT4KmjB5uf87eThrG2qpald0/gtvEOMw0psUHwxYihCGbfIpamIMSRaLnFW5grNzlNH2O9L8iabUd57r7JrH2zAKUAACAASURBVJhfhqvOx9qqWspHDeb77zn5zpxo/8qqo1Ff52/312GzRgsiC3OsVN45AY8/6v88+KmXVx+/RSLgVwliaQpCAsm25q46X6dyyFA4gscftUArX/+YWeMd/N/9bgDsGenU+4IcOeM3198wNIefPzSVd782E4fdSklsBvkNQ3PMlm5C/0dEU7imuJyAieG/DATDHZr5GtHzPJuVVQ/cxNaaehaURyPfv9xeS+W63aSnpZr3ubu0gDmTCjuJ8qoFU2SLfRUhoilcM1wo0vzeXnfS4yWFOfzqKzN49ouTzHpvV50PS1qK2dy3uMDOivllrN5aA8Bht5flc0u5eUxu9B4FWR0absD5ru6uOh/C1YOIpnDN0F2k+b29bp56bRevbz+R9No8m5VVWw7jqvNRuW43f//yDj4918IzG/bxxC938MgLVbjqfHzaFLNCY61kHr9tHOkK8myZncTaaPG2tqpWUoauIi5aNJVSsn8Qrnq6SgafM6mQf50/ka019Unn+Tjs0RZsa6tquWfSMCIaPM0hlt5dwokGP4fOBHhs7UdkZ0T/SqXFCohXb3Vx4/BBLJ5VnFSsK4rzJWXoKuOCoqmUmq6U2kesG5FSarJS6ie9/mSC0At0t0V/YMYoKu+4vtMQMuM6o3flu3tPE9Gab999A+5zLZxtiU6MbGmLEIzVkGdYUmkMnO9atGbb0S6tSRHMq4uLsTR/DMwFPABa6z3ArN58KEH4S+lOoOKT1BNJjGJXuRrMOvPKdbtpDIS4Z9IwmkMRfrL1MJsOdGz5+vczi7FZFD/84mTWVtWyeFYxudmWC5ZDClcPFyOaKVrr2oRj0rJN6LdcTGlhV+fjm/tWuRp49KXtuOp8LJo5lrpzLTz5yi7c56KNNk41teKsO1/lo4D/rDrOD++bwjt7TlN5x/UUF9i77FIkXJ1cjGiejI3Z1UqpVKVUJVDTy88lCJfNhUoLLzTjx1XnY/mG/TQ1t6E1nPA0s/J3h3A1NONtaWPV+66k99VAXraF8qIhrJhfRp7Nas5JFwYOF1MR9FWiW/RRQB2wJXZMEPotXQlid1tkp9vLsrf2cbzBz88fPt/ydd2OEzxYPopwe5j9p/1dXm9R57f3xudIkGfgcUFLU2t9Rmu9QGudH/tZoLVu6IuHE4TLoav+l0aPy/jXieeP1fuIaMjNtlBcYKdkqJ22sGbZ2/u6FMzPlxUAsPDWMfz04WmXXDsuXF1cTPT8BaXUmsSfvng4QbhUuvJnGilDxswe43XiumfvncyEAjt5tmjp5HU5Gcy/aRiWVJLy84em8rNHyvnOnBKemXujiOQ1wMVsz7fEvc4AvkDH+eSC0G/ozl9ZUpjT4Vz8a48/yNF6P4Oz0rGkpeCq8/Hkq7tYeMtofrzVRTAh9KkAaxoUObIBOlX7CAOXi9mevx73sxa4Fyjt/UcThMuju27mxrnlG/bj8Z9fk2ezMjIvm9xsC6FwhAOfnKOpuY01HxzB3xo211lSofL2Yp6+vZhgGI7VB3r9+wj9i8spoywCRvf0gwhCT5USGjN6uvJtGlv1+CR2jz/IsXofe040UVPn4ydbXaSnQns7xPc20hH4qLaRhyvG8LOHpjJnUmGnzxAGNhfj02xUSp2N/TQBm4mO4hWEHqOnxjYYIyIq77i+0xY9/jPybNYOrd4aAyECoQi/2n6cEbmZfPmvxhCOgI7oDvdo1/DIjOiccRHMa5NufZpKKQVMBj6JHYporXU3lwjCZXG5YxsS04gMK/JCn1HvC3YQzeICO0NzrITDGre3lVc/PEFEg0qBdAVtESh2ZGOzplFeNOTyvqQwIOjW0owJ5Nta6/bYjwim0GtcjmAueWNPB+vU6fayYmM1D73456SNNzz+oDm64rgnYPo136+u41RTK8c9fjLSUznjj9aNByMwNCeDtBTFN++cwAsLb5YI+TXOxfg0P1RKTe31JxGEyyA+SFPvC7Jyk5PFs4qTdkN3ur08+ovtPPnyTvJsVn71lRlmeaOrPtrTsh2FJ06EM9IU2ZY0frJgCuVFQ0Qwha5FUyllbN1vIyqch5RSu5RSu5VSu/rm8QSha1x1vg7WokFxbEJkosDl2az8fWxUbiLNoWhOkc2SShiYPmYwClBKcTg2skJG5QrQvU/zQ2AqML+PnkUQgAuXOxpr1lbV8r17ykxr0Wi24arzsWbbUZbOLjHP1fuCPLH2I/acOkemVfFK1XFe+fCEOW73rd3R1OMx+Vk0nvRS5w0yvsBGZno0q73IkZ00uCRce3S3PVcAWusjyX766PmEa4yLjaIbAZ/1O092WOvxB/naa7tpag516ovZ1BxCA+E2zY/ed5FtTSPPZiXPZmX44GiS+meuv47MNFg2+waGDc7k2S9O4tkvTmLlJmfSPpvCtUd3lqZDKfWNrk5qrZ/vhecRrnEuJYqebIJjSWEOLz8+wzwfX/HjbQ1TYLfS1BJEaaj1NOOq81FRnM8375zAV1/dReGgTCaOyKW8aEgHH6bRLk4sTaE7SzMVsAH2Ln4EoVe4GGEytvDP3Te5w/p6X5DGQIgVG6s7rC8pzOEH8ycSCIUJt8OwQRnYMlLJzbYA8EljtC/mL7fXYuSIJKYyiWAK0L2l6dZaf6/PnkQQLhJjC79ifhke/3n/Z70vSOW63ew64SEjLb3DOYDBWemE29uxpiuyrWkMzkoHolH1H//hMFlWxeLPFvPy9sSe24Jwngv6NAWhr7kYf6YhmI/+YruZj+mwW7ln0jBa2uCpz47rEASq9wX5zoZ9tIYhLzuDU43NHPc0s2JjNY2BEGkpKfyfL0bTiixpMqRV6Jru/nR8rs+eQhBidFc3nkiezcp1dmsH3+aovCxSFPiDYZxur3k/V52PU7Et+MThOaSlpgCaxbOKqSjO59XHb2HOpMKkW35BiKfL7bnW+mxfPoggON3eLuvG4zG25/MmD+PQGT+uOl+H9RENq9538cIHxygtzOEbd46nqbnNbO/2h5p6vvG58fz+YB3FBVH3fPwMHxFMoTtkHyL0C+IbbVxoCJmxPS9yZJOZlmIGcwysqcQS06E5FGbNtqMUObKZUGAjhWinot9V1yFFwcLlIKIp9AsMIbzYqY0ef7DDXHGICu93391PsB0emj6SH35xMoOzLGaS+19fn08EaAlHuLu0QHyXwmVxMZ3bBaFPiC+HvNAAtEd/sZ0hWekEgu2meHr8QWrPRsfrvrHjJDtqGxmUmU6ezYrT7eWZuTfSHGpnz6kmbhvvYGtNfe9+IWFAIv/UCv0Cp9vLQy/+mSfWfsRTr+zsNhBUUpjD9+4pI9OS1iHH44Oaelpj7d6CETh8xk8oHMFV52PhSx/y+vYTnD7XSlt7dF8ulqZwOcifGqHPSNaqzSDPZqUwJ4O29giH63ydmnDEU+8L8s6e03x2vIOIhqbmNup9QW4cPqjDumGDMnj2i5OoKM7nu397I89vqWHy8EHU1Pk5Vh+QKLlwWcj2XOgTnG4vC1/6kLWPTU/qt/T4g5xqbGHEkCzG5GcnLZGMZ2HFaFZvdZGZphiclc7yDfspKbCZ53OsKXiDbeb7OZMKGZyVTkVxPhOG5kjXdeGyEUtT6BNKCnO6FEyDtvYIaE1aStd/LI3Gw2u2HWVMXhYtYU1Tcxsr5pcxY2yeaQV4gxHa2zH9nUZXpHpfUART+IsQ0RR6HcM/mUwwjXN7TjTR3BYhEAx362s0ks9nlw7ltQ+j7dyMuvGVvztEOG7tE7cVmUJ5ueM0BCEREU2hV+mu1Vv8uZzMaB14OBKt0rmQuG2tqWdWST4Qrf5x2K08e+9EbhoxiHunDEMBM8bmdRBKEUyhJxDRFHqV7iw8oyemw26lvGgI19ksfNLUynffvXAZ5bzJw3Cd6ThzvKQwh6WzS/jDoXrsGWnmZwhCTyKBIKFXiO++3pVwGVVAK+aXsePYWTyBEMMHZ/BpbPBZsuvqfUGefHknNWd8hNqim/F9p5vM8xXF+bz6xC0cqw+wtqqW4gK7CKfQo4ilKfQ4ybbkySzH+HG7q7cepl3D9KIhaDrWNyZ2ZlcK2sLttMQcmGPzbR3W59msvLPntIynEHoFEU2hx0nckhsiWuVq6LDOsDQ9/iBfuGk4AO/tP004EumwxhBgY/0/3DGeonwbn4v5NLv6/IstyRSES6HXRFMpNVIptVUpdVApdUAp9XTcua/FplseUEo9F3f820opV+zc3b31bELvk9j1fGHFaCpf/7hDgrthaa7achh7RjQQlJuVQUpcmU+8NWqIIUSrfSqK8kkBbi8p6PbzBaEn6U2fZhj4ptZ6l1LKDuxUSm0GCoB5wCStdVApdR2AUqoUWADcCAwDtiilxmut23vxGYU+orjAbuZpxvs7SwpzWDG/jGVvfgzAGV8rGZbz/5Ybs8wBs4LHYbey9svTyc228NPcLCqKk1ucgtAb9JqlqbV2a613xV77gIPAcOCrwLNa62Ds3JnYJfOAdVrroNb6GOACpvfW8wm9j7GlNrbYxrElb+zp5Kf8H2PrHoERgzLNiiCH3crS2SWdSh5zsy08+ovtvLy9ViZECn1Kn/g0lVJjgCnAdmA8MFMptV0p9Uel1M2xZcOBk3GXnYodE65CnG4vS97Yw5I39gCY2/BkNeUlhTl8928nApCWqqj1NJvrDD9msmt+9ZUZPDJjtGzFhT6l11OOlFI24E2gUmvtVUqlAbnALcDNwHql1FiSzyTq1CZWKbUIWAQwatSoXntu4fIxhG7p7BJzjK7hmzRyKQ2hc7q9lBTm4GuN1om3xjoQfVBTT0lhTrd5no2BEF9/fTe/jNWUC0Jf0KuWplIqnahgvqK1fit2+BTwlo7yIRAB8mPHR8ZdPgI4nXhPrfUarXW51rrc4XD05uMLl4nHH2TF/DJTMOMj38ZIiypXg9nEo8rVwMb9nwKQmZ7ClJE5PP6Zceb9urIkc7Mt2CypnTq3C0Jv0pvRcwX8AjiotX4+7tQG4PbYmvGABWgA3gUWKKWsSqki4Hrgw956PqF3cLq9PPLidnYcO8vyDfvNwWYef9C0NBdWjObJV3YCsPax6VQU5/PZ62OWoo6w+6SXVb8/dMHPKinM4bVFFZJaJPQpvbk9vxV4FNinlPo4duyfgJeAl5RS+4EQsFBrrYEDSqn1QDXRyPtTEjnvv8RHwOOP5dmsjBiSxfqdJ80xE5V3XG9GwBfNHAuAL9Zx3dhW26zRP4pGwrrx/kKIYAp9TW9Gzz/QWiut9SSt9U2xn99qrUNa60e01mVa66la6/fjrvmB1nqc1nqC1vq93no24fIxttqJEXAjQu7xB0lPUaY/E6LC9tx9k7l/2kgqX/+YpuY2MiwdXdh1vtYO74fnZvX+lxGEy0AqgoSLJl4YQ+FIh3PxSejHPQEaA6FOpZTrd55k1QM3AZBKCqu3uswqoWJHdJTuhOuyqby9WHpeCv0WEU2hSxLzH+OFMbHnpdPtZeUmJ3k2Kz9eMIWK4vwOUe/4VKN/fnc/+TYLj8wYzdfX7abK1cCabUcAuHH4IH65vbbb0RiCcCUR0RSSkljzbRxbucnJio3VzC4d2qG2fOUmJ/7WMB5/kO//9iBOt7dDWpGRglRRnM+37pzAibPRxsFj8rLJzbYQjFmuO2sbaW+PJHkiQegfiGgKSYmv8zbE06jOaWoO8cy75xtwGMdtGWls2ufmwGkv3/r1xx1SjYzIOcDtpQXcNDKXwVnpWNJSaAyEaGqJ5mlOG53LWIftgjOCBOFKIaIpdLkVNpLS47fZeTYr6akpZKaldMiPzLNFhdNZ5+dzJfk88/lS8x7zJg8za86NYyvml7F6qwuI1qU/fmsRAG/vPo3rjL/baZSCcCUR0bwGifdVOt1eHv3F9m59iPGpRR5/kGxrGuMc53tYGlv5PJuVkgIb/+Vs4CtrP+LJl3fy3l43T7++m/f2uqlct9u0PhsDISxpKeZoiwlDc0iJBdTzsi1iaQr9FhHNa4BkqUHGsTyblfEF9m5FKtGnuXhWMempKazYWG1u2xdWRGvAK++awBO3jWGsIxtLWgpFjmxKC3MYnJXOwU+9uOp8VK7bzVOv7eL+aSNZWxUN+qzfeZICu5VRQzL5pKkFV52vd/+jCMJlIqI5wEnMqUzcbjvsVlYtmNLtSIrE1KGm5jaaQ2Ez7cjp9vK113bjdHt5b6+b/66p5/S5VhbPKibPZjW38TcMzaG4wM7yuaXcMDSH8qIhZrnlrWPzcHuDZKSl8OMFU6WWXOi3iGheg1xKV6D44WcOu5VFM8fy7Q17OVTnMzus59msjHVkc6w+wNfW7ebQmQBP/vU4iguiuZdLZ5ewtqrW3Irn2awsn1tqPsfyDfsZnpvFjYV2vntPmeRoCv0aEc0BjjEnPFnJo/G7qxG7cD5dyDhfXGBn8WeK+dVjM1jzpZvNddnWNIoc2fxkwRQKbBZuG+/o0BpuYcVo1mw7araMW7nJ2WFr/93fHOBLt4wx55QLQn9FplFeAxidhuLzKpe8sccUUyO1KJF6X5AVG6tNC9Hp9rL4lZ24GppZ+YWJ0a32hv2smF/G/dNGsmrLYUbmZlDnD/FmrPYcYOUmJ4FgGK2jVulz9002nwuiEyRXPXATa7Yd7dA2ThD6I2JpXgNcyJoEkp73+INUnz7H6q0us1vRcU8zBTlW3th1CledjxXzy/D4g3z3NwdYWDGaAnsGAP9dU8+KjdVAtElHempKhyqiRGE0tvISNRf6OyraYOjqpLy8XO/YseNKP8ZVQWJXosT3RjPg+ONOt5dvrd/DU7OKKS8agscf5OnXdnHoTICs9FTSUmH93/8VJYU5ZsrSsrf28fHJJr4zp4StNfU0NYcYnGVh+dxSUxCXb9hvJrsnWsBiZQr9AaXUTq11ebJzYmlegyRalIbf0rAmjVzKZW/u5eCnXpa9vZfKdbvJs1n5X9OifaLb29vRcc32j9UHWLXlMDePHgyAPSOdWeMd1JzxEwiGzYbEcH70RfzngUyQFK4OxKd5DWCIouG7NIIzhm9x5San2fvSSEeq90WT2Fc/OJXBWenm9vnXO0+SoiA1NYVRudEBaFWuBirXf8w/3jmePx31MC4/i1F5WVS+/jHfv6eM20sLzHsaPlDjc7oaZSEI/RWxNAcYiVZkfO23kTb03H2TOwRcQuEIjYFQh+1xfPs3o4Gwq87HkYYANmsqSkU7HXn8QSqK8/nHO8eztaaeRTPHMirWhGPtY9N5YMaoDjmhhkjGHxOEqwkRzQFEVwGfFfPLOnU4j08jag6FWfzqLirXRRPUjah5TZ2Pb6zfxZ5T59hx7CzFBXYGZ1rwtrbT1hYhFI6wcpOTKlcDP/3jEapPnwOigZ9VWw4nDeqISApXOyKaA4jE9CGn20vlut14/OfnjxtCaVieAIOzLCy5u4RHZow2x+WuWjCFf5pzAxnp6Tx+2xjmTCrEVecjEIrOo0hPTyXLksbS2SUUF9gpLcxh2ewbWLPtKKu3ujrcXxAGEuLTHIAs37CfhRWjWb3VRfXpcyx7cy/ZsZk7gWDYfG2UT94/bSTPbNgHCn764FTT//jajpNcZ7fwxq5PuL2kgOd+56SlLYI1Bf7hc9dz23iHacEun1vKqi2HWTRzLGu2HZXUIWHAIpbmAMMoe1xbVcvyuaX87OFp5GZHU34Wzyrmk6YW7pk0zFxf7wvy8vZafvrQVF59/Baz5ttV52PvySYCoXYisXLJf713EgU2C8EI/Pv7NfzzO+ddAUYQqaI4P2kFkiAMFEQ0ByCGgJUU5pCbbeG5+ybzQU09FcX5/HjBFN7de9oM8nj8QQ5+6qWpua3DPQ58co4IcKqplcKcTNZsO8qeE03U+UMAhMOaxBRfCe4I1wKyPR9AxEe/HfZoKtCTr+xi6uhB/Jcz2mXdnpEOnJ/xk2ezkpuZzj+9vReVonj18VsA+EmsQTDAfdNGMG/qCADKth+nqbmNT31BHru1SARSuOYQS3OAkBg5r/cFeX5zDeda2/gvZwPDBmVgz0hn6dv7KC20m/5MV52P2rPN3DdtBD+YN5GSwhxKCnP48l+NMe/9650ngagQ/9v/uonRedn85IEp0o1IuCYRS/Mqo6tSw/gWbhDddq+YX8amfW5W/8HF8r8pZXBWOikKXvqf49xeUgBAbraFNAW/+OA46zNOmYnsfzjcYN77sVvHdhh3YUlLobxoSB98W0Hof4ileRXRXeMNo6u60+01R1h8c/3HfFTbyE8WTGXOpEJysy3cMNTOxGGDosGhDfs5Vh8g1A4aKMzJMOf2jM3PMu89Ki/L/EyPPyiBHuGaRkTzKsKwHrsSLH9r2Ows9L17ynB7WwkEo3mVVa4GVmysJsuSxrNfnERJYQ6zxjvwtrQRAb4wZRhP3zHe9HVmpqea9z3haWb5huj0yYde/LMMPROuaUQ0rxKcbi8LX/qwW8GypKUQCIZZuclJedEQvjRjNG0RzZOv7uJ/v7KTpuYQ3pZo9Pu9vW6Wvr0P97mW6Pt9p3nmnX0smjkWgP+Obc8/V5LP7aUFrJhfRm62BXX1NsUShB5BRPMqoaQwh7WPTe8yadxhj46QyLamsWjmWHYcO8uP3ndR721FA1+uGAPAoTMBlr25F4AUBftONwHQGoald0Wrezz+IFmWqKVZUZTP8g37zWd45YlbOpVkCsK1hIhmP+JCYx7ybFaWb9jPe3vdna6p9wXJs1kJhNpZs+0oAJkWRVNLNP9ywtAc/mpcHgCfHe8AIKLBmhYVx+tsFiaPGsySN/awYmM16amKVAW/q67rEGASwRSudSR63k+Ib5vWlc/SYbcyb/Iwvv76bn4Zi3IveWMPS2eXsGrLYeZNHsaps838/cwyXvrTMWbfWMhbu08DsPJ3ToJt7ea9youGcNOIQYy/zs5vqSPTmtZhFAVEq4KKC+wS9BGEOEQ0+wnJeksmSy8qcmRjs6SaY3FD4Qh5tmi60cpNTn7y4BQAPj7ZxAF3E4poZPzE2WYy06NNg2eMjVqcD948ilF5WfC+i8GZ6Xj8wQ6WpIilIHRGtuf9hGTjJxLTi+p9UVH72cPTTHEzouN5NiuhcITiAju52RasadAc0hhxm1RgZG4mAE+9uosvv/Qhy97exwlPMzePzjWtVZkEKQjdI5ZmPyBxa57YOBjOd1+fNd7Bu3tPU1xg5/3qOg6cPmduo0PhCK46H89vrkGTApzfjrdpOFTXDIA1LYV/u38ye0408cCMUWZnddmKC8KFEdHsByRuzRPfxwvmd97djyUVdhw7y3O/czKhIMeMeDvrfKzc5CQ9NYVQnP8SMLfpAHXeII2BEFtr6k3BND5XEITuke15P6ErX2a81Tl51GCe+sw4rOlpDM5Kp6Qwh2f+5gYgOtjM1xqmpa0dh91CW0I+pQa+M6eEKSNzmDY6l+ICu8znEYTLQCzNfkaVq4G1VbXmiFvD6nTV+Vj86i7a0fzr/ElUFOeTm21h2Zt7SU9NoTEQ9UWePBugtsHf6b7DcqzMmxrtVmQMVhPBFIRLRyzNK0Cy4WcQFczK1z9m3uRhnYIya6tq+f78iRTlZTM4K9re7Vh9gAOnz+FracPtC5KioLlN09pxZw7A8NxoLbkxWE0EUxAuDxHNPqTeF+wUFTeCQE63l7VVtXz3b29kzqTCpD5OAFd9gCd+9RFVrgbW7zzJ0rtLWPXgFEqH5jA0J6PL/6H3TR0hvktB6AFENHuYrlJ26n1Blryxh5WbnB2i4kZLN4gOO1u/82TS/ExXnY9lb+0l3N5OcyhCU3MbgWCY//jvaEf1W8flcfpcK6kJ/0enjxmM3ZrK5FGDe/7LCsI1iPg0e5DuqnqMbbHxOvGao/V+Vj801RxDEY/T7eX5zTW0E2Gcw0ZaSrSfZZEjm19VHeeZd/Zhi9WKt0UgTUE4Fgja94mXkqF2GXQmCD2EiGYPkqyqJ/F8PIZVaklLYWRetnnc4+8YPV+5yQnADQWDeOzWItOn+a31e9jv9pIONMU5MsMa7Bkp+FojBMMRFs0cJ1tyQeghlE6cjnUVUV5ernfs2HGlH+OyMHIvV8wvw+OPWpuWtBRC4Qi2jLQOwZp6XxCPP5pb+dVXd9LWHuHFR29m9VYXs8Y7uHH4IJa/s59AMIzbe949kKrg27NLePwz467U1xSEqxKl1E6tdXmyc+LT7EW6ipI73V5WbnKysGI0DruVPJsVpaKzw1fML2Pp7JIOlqGrzseqLYfJzbZQYLMSCEY44YlW9/x650lOeJr55p0T8MVKKp+4bQy3FQ+hXcPaP9dKaaQg9CAimr1E4mgKI2puCKa/NcyabUdNK/JofYBj9QGWb9jPyk3ODtd99ZWdLKwYDcCcsvPDzGaNd3DoTIClb+9j9dbDDLVbqby9mGfm3sjLj1fw8PSRvPXkrbI1F4QeRHyavUS8f9PYihsJ6/Ht1xx2Kx5/kH+8awIvb6/lSIOfpXedtzQbAyECoXYOfHKOn/33EUblZpFtiUbDN+0731fzr8bl8eIHx3nhg6PMnlhISWEOP7h3Up9/b0EY6Iho9iKJzTeMCHZ8TXl8pc/PH5rGCU8zz2+pYVReFhXF+RQX2Bl/nY3fH6xjVG4WcycWsrWmnjyblQlDz7dxmzpqCM/MyeCnfzhyRb6rIFwr9Nr2XCk1Uim1VSl1UCl1QCn1dML5bymltFIqP+7Yt5VSLqXUIaXU3b31bH2JYXEaXdfjq38eevHPrN7q4vvzJzL+umhLt3f3nuaeyYWs2XaU17efMO/jbQnR0tbOD95zcvPoXBx2K95YV/YU4KU/HePG4YMoHZYj6UWC0Iv0pqUZBr6ptd6llLIDO5VSm7XW1UqpkcCdgKkKSqlSYAFwIzAM2KKUGq+1TlIU2L9Ilowef8z4Hd/6bfVWF6Nys1g8q5jcbAvZ1jQaAyF2nzrL/xzxMHxwBlsP1VN5roUDbh8pQEGO6keNuwAAEL9JREFUFQ2set/FhKE5bKr+lMIcK/5QmMduLWJtVS3L55aKD1MQepFeszS11m6t9a7Yax9wEBgeO/3vwBLOdysDmAes01oHtdbHABcwvbeer6dINou8q/nkiWLWFtGs3upi+Yb9LJ0dHWpWVjiY8ddl88P7JlNynY1AKBoRz8lIJScjrcO/cvdPG0mDP8joIVmUFw1hxfwymeEjCL1Mn0TPlVJjgCnAdqXUPcAnWus9CcuGAyfj3p/ivMj2W5IltCfrhxkvoA67lcWzijnd2MI9k4ZxtN7PnhNNOOxWahv85NszaGpu40wgSHMoamg3tbZzsrEFrWBcflQk50wq5JePzeD/+/IMHHarWJiC0Af0eiBIKWUD3gQqiW7ZnwHuSrY0ybFOmfdKqUXAIoBRo0b13IP2Ek63lwfXVDEmP5s1X7rZFLbcbAvhSIRReVnk2608884+Vm6q5mxLOw1+D03NIQpzMqi8cwItbe1sjgWCJgy1s+ivz1f4VBTnd/fxgiD0ML1qaSql0okK5ita67eAcUARsEcpdRwYAexSSg0lalmOjLt8BHA68Z5a6zVa63KtdbnD4ejNx78oLrQ9bwyEiGg43tCMx39+TWMgGthpam5DaUhPAWv6+X/DDrp9nGhsxlXn41hDgFG5WQSCYd7afZqnX9slCeuCcIXoNUtTKaWAXwAHtdbPA2it9wHXxa05DpRrrRuUUu8CryqlnicaCLoe+LC3nq+n6G57DrB6q4tn751EkSPb9Dc63V6amtuwZaThbWnDecZPqoK2WLOOdsCSAuPybTQ1t7H/Ey/jC2z8x6PlfFBTz7y4Nm+CIPQtvWlp3go8CtyulPo49vP5rhZrrQ8A64FqYBPwVH+NnF8owGMc8/iD7D99jpf+dMw87nR7eeiFP7Ps7b0s/kwxvtZo2tCCm0fyyC2jzXWhCMydWEh50RDGD7UzOMtCns3K45+R5huCcCXpzej5B1prpbWepLW+Kfbz24Q1Y7TWDXHvf6C1Hqe1nqC1fq+3nu1SuZjIeLK1jYEQWmta2sKs2FhtjuD90i2jiUQ0z/9XDb+rruPzZQUc9zTzUW0jIwZnmNe/sfMUHn+Q6+xWSSUShH6C1J5fgGQi2VX7t/i19b4ga7YdZdSQLNJTU2hrj269q1wN/OyPR7CmpdIcinCqsZnf7q/jXEuIxbOKmTRiEAA51hTqYp/53H2TJZVIEPoJIpoXIN5naYjixax12K3cP20kWZY00lJSMDrwFRfYWXLXBBoCIQDOxoJD6amK1Vtd/OmIB4BhgzMZnZdFnk1SiQShPyGieREYgnmhJsPxON1elr+zn1A4wn1TR2DLSMNV52PJG3uo87Wa67Izog2FQ+3R1nD/8fA0bhoxCHtGtEpIEIT+hYjmRRC/7e5qW564rqQwh588OAVLWgorf+fk1rF5PL+5hjPeVl744DgAgzPT/v/27j04qvKM4/j3CRBCSCAgMSYImpIiA8jdW6tWbMeitXipdbR11GGm9VLx0vFaqkNLOxXohXHUttp2SqtVW2rF6nhDxUsHEQEDUWMIAgpkIKJCEAjEvP3jvJssSy57IGd3w/4+Mzs5e/Zsnidvdp5933N5D+OHBsPx40oKGFHaj1MqBvHgFSdw/2UTmXfJePUyRTKMimYSEnuY7R0YShzKV5QU8vPzRlNWlMd9i2t5t247l598LL38afw7djex+P2Pyc2Bl6q3Ul23oyWervARyUwqmkmKL5jxB4YSC2Xs560LKrnx0ZUA9O+Ty+QRR7K3yVG3fTexc9ibgemTK5g/7SRGlfXX7EQi3YCKZpLaKpIx8QeJYr3F26aMAOD1mno+/GQXj6/czKCCXOa/sQEXd8Fo/c5GTqkYpKG4SDehopmEtnqXiWL3L5/9bDW3LgjmIllXv5NfPFPNR5/uBuCo/nnc/70J/PaiCRxXUkCOwbljytr9nSKSeVQ0k9DRPs3E57dNGcGci8byek09m/2dIfvlBfck/9boUipKCllYuZmZ3x7FmMH9qSgpTNFfISJdQUUzSbEh+JLaj9u8Ydozq+q47E9LuWthFbVbGni5pp5cPwzf19TMMQP7tFwzPuv80Qzom8uAvrlp/ItE5GCoaCapum4H1z60nOsfXcmkoUUtRXTWU+8yaWgRCys3c/XpX2LfF8088NoHTB1TRo/gFEx2NzkK8/Y/53LeojUH3KpXRDKfimYS6hsamf1sNbk9c7jqtC8x5/n3WVL7MbVbGqjatJ27n6tm0tAi7n91LbVbG7h44hAeWfYhM791PMMG5TP7guMpLmy9pjzW29SlkSLdjy45SUJxYe+W2+7WbmlgZGk/BvTNZd6iNUyfXMHTVXWcN+FoBg/I544nVrHp011UbtzOnn0fsG3XXsYOLeLMkSUHHHEXke5HRTOEbTsbmf7ISgr8UPu8sWX85InVHDMwH4Dy4r6MKu3fUkCL8ntx78u1un5c5DCi4XkS6hsa+cH8ZRxR0JsLxpexftsurnt4OfctXsMXrpnde5t46d0tzH62mjvPHQkEBVR3hxQ5/KinmYS31n1C5cbtPPBKLcs2fAZAU7Oj7rM9TDulnHteqmXGwtWMHlzEp5/vZfaz1UBw+pH2W4ocXlQ0k1CU3wsH/Pn19Qz2kwR/smsfsy8cw6TygSyuqeeq04dRXty35ai4huQihycVzTAMcnKMHOD6yRWcPaaU+oZGBvTNZVL5wFBTx4lI96R9mkl4Z9N2AIryejIgP5e8Xjm8XFPfMkHHnIvGthRKFUyRw5uKZhJOHV5MXg9jx54mzhhezJ59zUwdU6ZCKZKFNDxP0vFHFzHtq+XBE4MFKzYecO6liBz+1NPsROxqoF17mygv7ssfX11Ls4OLdO9xkaykotmJ2A3SarbuZF3951x1+jBygH59eqU7NRFJAxXNJEwqH8io0n78c/lHAOT37sFDSze0e+9zETl8aZ9mkh684gS27QxmNerZw7jspGM0PBfJQuppdiI2a/u2ncEdJuddMp7rvlbBwsrN6mmKZCEVzU7E38aivqGR2i0NzH2hhsnDi9XTFMlCKppJiL9LZEVJIfMuHtdycruIZBcVzSTEz6d55xNVTCofqMslRbKUDgQlKVYgZ50/er/nIpJd1NM8CLcuqNTQXCRLqWiKiISgohlScWFv3UVSJIupaIZU39DIvEVrNDwXyVIqmiFpomGR7KaieRBUMEWyl4qmiEgIKpoiIiGoaIqIhKCiKSISgoqmiEgIKpoiIiGoaIqIhKCiKSISgoqmiEgIKpoiIiGYcy7dORw0M6sHNkQcZhDwccQxMj2HdMfPhBzSHT8Tckh3/FTmcIxzrritF7p10UwFM3vLOTcpm3NId/xMyCHd8TMhh3THz5QcNDwXEQlBRVNEJAQVzc49kO4ESH8O6Y4P6c8h3fEh/TmkOz5kQA7apykiEoJ6miIiIWR90TSzIWb2spm9Z2bvmNkNfv0sM1tlZm+b2fNmVhb3njvMrNbM3jezb0YRP+71m83MmdmgKOJ3lIOZzTSzTb4N3jazc6LIoaM2MLPpPsY7ZjYnivgd5WBmj8X9/evN7O0ocugg/jgze8PHf8vMTowific5jDWzJWa22sz+a2b9osjBzPLM7E0zq/Txf+bXDzSzF8xsjf85IIr4SXPOZfUDKAUm+OVCoAYYCfSL2+Z64A9+eSRQCfQGyoG1QI+uju+fDwGeIzgXdVAU8Ttpg5nAzW1sn5I2ACYDi4De/rUjU90GCdv8BrgrxW3wPHC2X38OsDgNn4NlwNf8+mnArIjawIACv9wLWAqcDMwBbvfrbwdmR9UGyTyyvqfpnKtzzq3wyw3Ae8Bg59yOuM36ArGdv+cBjzrnGp1z64Ba4EQOUnvx/cu/A26Ni93l8ZPIoS2paoNrgLudc43+ta1RxO8kBwDMzICLgUeiyKGD+A6I9ez6A5ujiN9JDscBr/rNXgC+E0UOLrDTP+3lH87Hme/XzwfOjyJ+srK+aMYzs2OB8QTfcJjZL83sI+D7wF1+s8HAR3Fv20jHBeag4pvZVGCTc64yYbPI4ifm4Fdd53dT/CVuWJSSNgCGA6eZ2VIze8XMTog6fhs5xJwGbHHOrYk6h4T4NwJz/efw18AdUcdvI4cqYKp/6bsEI6BIcjCzHn4XyFbgBefcUqDEOVcHQWEHjowqfjJUND0zKwD+DdwY62U652Y454YADwPXxTZt4+2HfApCfHygCZhBa6Heb9Mo4ifm4Nvg98AwYBxQRzA8jSyHNuL3BAYQDNFuAf7pe3ypbIOYS2ntZRJVDm3Evwa4yX8ObwL+HGX8dnKYBvzIzJYTDNv3RpWDc+4L59w44GjgRDMb3VGqXR0/GSqagJn1IviQPOyce7yNTf5B65BkI63ftBD8czcf8I5Diz+MYB9NpZmt9zFWmNlRUcRvJwecc1v8h7gZeJDWoU8q2iAW53E/bHsTaCa49jhlbeDX9wQuBB6L2zxVbXAFEFv+FxH+D9rLwTlX7Zw7yzk3keCLY22UOfiYnwGLgSnAFjMr9fmVEvRCI43fWXJZ/SD4tvobMC9h/ZfjlqcDC/zyKPbf+fwBh77z+4D4Cdusp/VAUJfG76QNSuOWbyLYf5SyNgCuBn7ul4cTDMUslW3gX5sCvJKwLlVt8B5whl/+OrA8DZ+D2AG4HP/6tIjaoBgo8st9gNeAc4G57H8gaE5UbZBUnlEHyPQHcCpBl34V8LZ/nEPwbVvl1/+X4OBQ7D0zCL5t38cf2ezq+AnbrMcXza6O30kb/B1Y7dc/yf5FNPI2AHKBh/z/YQVwZqrbwL/2V+DqNt6TijY4FVjui8NSYGIaPgc3EBxJrwHuxl8UE0EbjAFW+vhVtJ6pcATwIrDG/xwYVRsk89AVQSIiIWifpohICCqaIiIhqGiKiISgoikiEoKKpohICCqaklHM7As/o0+Vmf3LzPIP4XedYWZP+eWpZnZ7B9sWmdm1BxFjppndfLA5SvejoimZZrdzbpxzbjTB5XpXx79ogdCfW+fck865uzvYpAgIXTQl+6hoSiZ7Dagws2P9HI/3E5zkPsTMzvJzPK7wPdICADObYmbVZvY6waWP+PVXmtm9frnEzP7j522sNLOvEJy0Pcz3cuf67W4xs2V+wpKfxf2uGX7+xkUEMwBJFlHRlIzkr/c+m+CKJAiK09+cc+OBz4GfAt9wzk0A3gJ+bGZ5BNfIf5tgVqKj2vn19xBcFjkWmAC8Q3B53lrfy73FzM4Cvkxwrfc4YKKZnW5mE4FLCGYAuhA4oc0Ictjqme4ERBL0sdbZ0V8jmNWnDNjgnHvDrz+ZYALa/wWTHpELLAFGAOucn77NzB4CfthGjDOByyGYVQfYHj8buHeWf6z0zwsIimgh8B/n3C4f48lD+mul21HRlEyz2wVTg7XwhfHz+FUEcy1emrDdOLpuajADfuWc+2NCjBu7MIZ0QxqeS3f0BvBVM6sAMLN8MxsOVAPlZjbMb3dpO+9/kWCeytikt/2ABoJeZMxzwLS4faWDzexIghnMLzCzPmZWSLArQLKIiqZ0O865euBK4BEzW0VQREc45/YQDMef9geCNrTzK24AJpvZaoIZhEY557YRDPerzGyuc+55gnlUl/jtFgCFLrgdxGMEMwD9m2AXgmQRzXIkIhKCepoiIiGoaIqIhKCiKSISgoqmiEgIKpoiIiGoaIqIhKCiKSISgoqmiEgI/wchTtqQ9eySNQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "rmse: 2.5904884\n",
      "r2_score: 0.98544579149015\n"
     ]
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.scatter(y_hat,y_test,s=0.1)\n",
    "ax.plot([y_test.min(), y_test.min()], [y_test.max(), y_test.max()],c=\"black\")\n",
    "ax.set(\n",
    "    xlabel=\"Predicted\",\n",
    "    ylabel=\"True\",\n",
    ")\n",
    "plt.show()\n",
    "print(\"rmse:\",np.sqrt(metrics.mean_squared_error(y_test, y_hat)))\n",
    "print(\"r2_score:\",metrics.r2_score(y_test,y_hat))"
   ]
  }
 ],
 "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.4"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}