{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <p style=\"text-align: right;\"> &#9989; Put your name here</p>\n",
    "#### <p style=\"text-align: right;\"> &#9989; Put your group member names here</p>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Day 04 In-Class Assignment: Analyzing Covid Data\n",
    "\n",
    "This notebook is an modified version of one authored by [Brian O'Shea](https://web.pa.msu.edu/people/osheabr/) at CMSE in the early part of the COVID-19 epidemic.\n",
    "\n",
    "## Goals for today's In-class assignment \n",
    "\n",
    "1. Read JSON data into a Pandas DataFrame\n",
    "2. Perform calculations with the `datetime` class\n",
    "3. Read COVID data from the web into a Pandas DataFrame\n",
    "4. Calculate the Doubling Time for the virus\n",
    "5. Represent that data in plots and discuss meaningful representations of data\n",
    "\n",
    "## Information on the data we will use\n",
    "\n",
    "This notebook uses Pandas to retrieve JSON data files created by the [Johns Hopkins Coronavirus Resource Center](https://coronavirus.jhu.edu) (and stored at [The COVID Tracking Project](https://covidtracking.com/)) to make real-time tracking plots of the growth of the COVID-19 pandemic in the US and Michigan, as well as estimate the doubling time of the number of positive cases.  These data files are updated daily, so re-running the notebook will get the most current data.  It's easy to change the state that is plotted by changing the `my_state` variable.  There's also quite a bit more information in the data files!\n",
    "\n",
    "**Important note:**  The data file tracks the number of people who have _tested positive_ of COVID in the US and Michigan, not an estimate of the total number of positive cases in the US and Michigan."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The imports "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import datetime"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## 1. Working with JSON data\n",
    "\n",
    "It is a bit of a story to explain what JSON (pronounced \"J\" \"son\") data is but it is worth our time. To understand JSON you have to know a little bit about JavaScript. If you read the <a href=\"https://en.wikipedia.org/wiki/JavaScript\"> wikipedia page </a> you will read:\n",
    "\n",
    "> Alongside HTML and CSS, JavaScript is one of the core technologies of the World Wide Web.[8] JavaScript enables interactive web pages and is an essential part of web applications. The vast majority of websites use it for client-side page behavior,[9] and all major web browsers have a dedicated JavaScript engine to execute it. \n",
    "\n",
    "So the fact that modern browsers have a dynamic character (pages change based on all kinds of inputs) is largely due to an underlying programming language called JavaScript. \n",
    "\n",
    "One problem is that it is likely the case that you want your web server (not on your machine) to dynamically generate content, to be rendered using JavaScript. How to get the dynamic content from the server to your local machine? The exchange of information has to use text to be universal. How do you translate existing JavaScript objects into text and back? You guessed it, that is what JSON is.\n",
    "\n",
    "From <a href=\"https://en.wikipedia.org/wiki/JSON\"> wikipedia JSON page</a>:\n",
    "> JavaScript Object Notation (JSON) is an open standard file format, and data interchange format, that uses human-readable text to store and transmit data objects consisting of attribute–value pairs and array data types (or any other serializable value).\n",
    "\n",
    "> JSON is a language-independent data format. It was derived from JavaScript, but many modern programming languages include code to generate and parse JSON-format data. The official Internet media type for JSON is application/json. JSON filenames use the extension .json. \n",
    "\n",
    "The term `serializable` is worth noting. By this we mean that actual object instances, with all their variable settings and methods, can be stored into a format and then re-created \"as they were\" in another location. That is, the `state` of a system can be stored and then re-created at a later date. Note that JSON was originally set up for JavaScript, but as a \"state saver\" can be used in all kinds of place. Like Python!\n",
    "\n",
    "Bottom line: JSON is a way to exchange information across the internet. We are going to use that!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.1 DataFrame to JSON and back\n",
    "\n",
    "The `.read_json` and `.to_json` allow us to read and write a DataFrame to JSON. Let's do that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_dict = { \"tens\": [10, 11, 12, 13, 14, 15],\n",
    "            \"twentys\" : [20, 21, 22, 23, 24, 25],\n",
    "            \"thirtys\" : [30, 31, 32, 33, 34, 35]\n",
    "          }"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Do the following:\n",
    "* turn that dictionary into a DataFrame\n",
    "* print the DataFrame\n",
    "* write out the DataFrame as a JSON file (look up `to_json`)\n",
    "* take a look at that file with your editor\n",
    "* read that JSON file back into a DataFrame with a different name (look up `read_json`)\n",
    "* print out the new DataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1.2 Read JSON over the internet\n",
    "\n",
    "But here is the really cool stuff. You don't have to have a JSON file on your local machine to read it in. You can actually read it from a source on the internet and use it locally. Remember, JSON is an exchange file format so that makes perfect sense!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     sepalLength  sepalWidth  petalLength  petalWidth    species\n",
      "0            5.1         3.5          1.4         0.2     setosa\n",
      "1            4.9         3.0          1.4         0.2     setosa\n",
      "2            4.7         3.2          1.3         0.2     setosa\n",
      "3            4.6         3.1          1.5         0.2     setosa\n",
      "4            5.0         3.6          1.4         0.2     setosa\n",
      "..           ...         ...          ...         ...        ...\n",
      "145          6.7         3.0          5.2         2.3  virginica\n",
      "146          6.3         2.5          5.0         1.9  virginica\n",
      "147          6.5         3.0          5.2         2.0  virginica\n",
      "148          6.2         3.4          5.4         2.3  virginica\n",
      "149          5.9         3.0          5.1         1.8  virginica\n",
      "\n",
      "[150 rows x 5 columns]\n"
     ]
    }
   ],
   "source": [
    "iris_df = pd.read_json(\"https://raw.githubusercontent.com/domoritz/maps/master/data/iris.json\")\n",
    "print(iris_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ah, it's our favorite iris dataset but read in from a file on the internet. Open that url once in your browser. Should look pretty understandable at this point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-----\n",
    "## 2. Analyzing COVID data\n",
    "\n",
    "You must be living under a rock to not know about tracking the COVID-19 virus outbreak. Let's see if we can do some data analysis of a time-series, JSON file. \n",
    "\n",
    "Gets the data from [covidtracking.com](https://covidtracking.com) - specifically, the page [https://covidtracking.com/api](https://covidtracking.com/api).  This page has historical data for the US as a whole, for states, and for counties.  The JSON files have a ton of data in addition to what we're showing. The two URLs we are most interested in are:\n",
    "* <a href=\"https://covidtracking.com/api/us/daily\"> https://covidtracking.com/api/us/daily </a> which is the US daily report as a JSON file\n",
    "* <a href=\"https://covidtracking.com/api/v1/states/daily.json\">https://covidtracking.com/api/v1/states/daily.json </a> which is the data for all 50 states.\n",
    "\n",
    "Before you do anything else, you should open those two links in your browser and look at the data. JSON is pretty readable and you need to read a little about what's there so we can do some more work."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Do the following:\n",
    "* read the US daily into a DataFrame\n",
    "* read the state data into a DataFrame\n",
    "   * make a DataFrame that contains just the state of Michigan from the state data\n",
    "* print the `.head` of the US and MI DataFrames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.1 Convert data to DateTime\n",
    "\n",
    "There's a date field in the data but it is not in pandas datetime format. That would be helpful since this is a time series. \n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Let's do what we did before:\n",
    "* Let's make a Series from the `date` field using `to_datetime`. \n",
    "Look at the format of the date of the JSON file in your browser and write the appropriate format specifier to get it into datetime format.\n",
    "* assign the Series to the DataFrame under the name \"datetime\"\n",
    "* make a new DataFrame with all the positive cases in the US along with the datatime (just the two columns)\n",
    "* print the head of the 2-column, positive US cases\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do the same for the Michigan State Data. \n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Make a two column DataFrame with datetime and positives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "#your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.2 Plot cumulative positive tests for US and MI\n",
    "\n",
    "Again, this data is large and thus we need to represent in a way that makes it understandable to us.\n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Using the `.plot` method, make two plots: one for the MI positive cases and one for US positive cases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x7fe088a845e0>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEiCAYAAADnMZWTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xUVf7/8dcHEggIhI5A6FU6EhCW3RV7/6EuCpYFBcW26hbX1dWt6lrWtayuha+F4koRG+piWVzFgkCQ3lskAYRA6JCQ8vn9MTc6YAiQNjPJ+/l4zGNmzr3n3s9Nmc+ce+49x9wdERGRI6kS6QBERCS6KVGIiEiRlChERKRIShQiIlIkJQoRESmSEoWIiBRJiUKkBMxsr5m1LWL5UjMbVI4hiZQ6JQqJemaWamYHzazhYeULzMzNrHXwfqyZ3V+esbl7LXdfd6T9u3tXd/+kPGMSKW1KFBIr1gNXFLwxs+5AjciFI1J5KFFIrJgADA97PwIYf6yVzWyQmaWb2e/NbFvQSrkqbHmimY03swwz+8bM7jWzKsGy9mb2qZntCupODqvnwfLRwFXAncHpqHeC5almdqaZNTOzA2ZWP6xu72B78cH7kWa23Mx2mNkHZtaqiOP5sZl9aWY7zSzNzK4Jyi8ws/lmtjso/3NYnQQze8XMtgf15ppZk7Djf9HMNpvZRjO738yqHu34pXJQopBY8RVQx8xOCj7AhgKvHOc2TgQaAs0JJZoxZtYpWPYUkAi0BU4llJSuDZbdB3wI1AOSgnUP4e5jgH8DjwSnoy46bPkmYBbws7DiK4Gp7p5jZhcDvwcuBRoBnwETCzsIM2sJTA/iaAT0AhYEi/cFsdcFLgBuCrZNcMyJQAugAXAjcCBYNg7IBdoDvYGzgeuO9filYlOikFhS0Ko4C1gBbCzGNv7g7tnu/inwHnB5WOK52933uHsq8A/g50GdHKAV0Mzds9z982LG/yrB6TMzM2BYUAZwA/Cguy9391zgb0CvI7QqrgL+6+4T3T3H3be7+wIAd//E3Re7e767LyKUbE4NO44GQHt3z3P3ee6+O2hVnAf80t33uftW4PEgvtI8folRShQSSyYQ+hZ+Dcdx2inMDnffF/b+G6AZoVZGteB9+LLmwes7AQPmBFcxjSzGvgGmAgPMrBnwU8AJtRwg9EH8ZHBKaCeQGeyzeSHbaQGsLWwHZnaKmf0vOIW2i1CroeAigAnAB8AkM9tkZo8Ep71aAfHA5rD9Pw80LuXjlxilRCExw92/IdSpfT7wRjE2Uc/MTgh73xLYBGzj+2/N4cs2Bvv91t2vd/dmhL75P2Nm7QsL8Sjx7yR0CudyQglvon8/fHMacIO71w171HD3LwvZVBrQ7gi7eRWYBrRw90TgOUIf8gStj7+4exfgR8CFhFpoaUA20DBs33XcvetxHr9UUEoUEmtGAacf1jI4Hn8xs2pm9hNCH5SvuXseMAV4wMxqB6d7fk3QB2Jml5lZUlB/B6GEkFfItrcQ6uMoyquEPpx/xvennSD0gX63mXUN9ploZpcdYRv/Bs40s8vNLM7MGphZr2BZbSDT3bPMrB+hhESwzdPMrHtwqm03oeSY5+6bCSWwf5hZHTOrYmbtzOzU4zx+qaCUKCSmuPtad08pZvVvCX3QbSL0YXuju68Ilt1KqCN4HfA5oQ/xl4JlfYHZZraX0Lf12919fSHbfxHoEpy+eesIMUwDOgBb3H1h2HG9CTxM6LTQbmAJoX6DH3D3DYRaVb8hdIpqAdAzWHwz8Fcz2wP8kVACLHAiodNfu4HlwKd8f0HAcEKn35YFP6OpQNPjPH6poEwTF0llYKG7o19x96SjrSsih1KLQkREiqREISIiRdKpJxERKZJaFCIiUiQlChERKVJcpAMobQ0bNvTWrVtHOgwRkZgyb968be7eqLBlFS5RtG7dmpSU4l5mLyJSOZnZN0daplNPIiJSJCUKEREpkhKFiIgUqcL1URQmJyeH9PR0srKyIh1K1EtISCApKYn4+PhIhyIiUeKoicLMXiI0yuZWd+8WlP0duAg4SGhc/GuDIZQxs7sJjfCZB9zm7h8E5X2AsYTmOf4PoYHF3MyqE5pboA+wHRgaTByDmY0A7g1Cud/dxxXnINPT06lduzatW7cmNF+MFMbd2b59O+np6bRp0ybS4YhIlDiWU09jgXMPK/sI6ObuPYBVwN0AZtaF0KxYXYM6zxTMuws8C4wmNHJmh7BtjiI0oUx7QrNqPRxsqz7wJ+AUoB/wJzOrd/yHCFlZWTRo0EBJ4ijMjAYNGqjlJSKHOGqicPeZhIYyDi/7MJiuEUJzGReMyDkYmBRMNbkeWAP0M7OmQB13nxVM1DIeuDisTkFLYSpwRjBN5DnAR+6e6e47CCWnwxPWMVOSODb6OYlUPnn5RQ/lVBqd2SMJTfQOoWkb08KWpQdlzYPXh5cfUidIPrsIzet7pG39gJmNNrMUM0vJyMgo0cFEs+eee47x40MzgI4dO5ZNmzZ9t+y6665j2bJlkQpNRGJUVk4eN70yr8h1StSZbWb3ALmEJoGBYMrFw3gR5cWtc2ih+xhgDEBycnKFHeXwxhtv/O712LFj6datG82aNQPghRdeiFRYIhKjUrft487XFzFnfWaR6xW7RRF0NF8IXBU27286oYnfCyQRmk0sne9PT4WXH1LHzOKAREKnuo60rZiUmppK586dGTFiBD169GDIkCHs37+fGTNm0Lt3b7p3787IkSPJzs4G4K677qJLly706NGDO+64A4A///nPPProo0ydOpWUlBSuuuoqevXqxYEDBxg0aBApKSk8++yz3Hnnnd/td+zYsdx6660AvPLKK/Tr149evXpxww03kJen2SxFKqODufnc8dpCTv/HJyxM28mTw3oVuX6xWhRmdi7wO+BUd98ftmga8KqZPQY0I9RpPcfd88xsj5n1B2YTmnbxqbA6I4BZwBDg4+BqqA+Av4V1YJ9N0GleEn95ZynLNu0u6WYO0aVZHf50Udejrrdy5UpefPFFBg4cyMiRI3nsscd4/vnnmTFjBh07dmT48OE8++yzDB8+nDfffJMVK1ZgZuzcufOQ7QwZMoSnn36aRx99lOTk5B8sGzBgAI888ggAkydP5p577mH58uVMnjyZL774gvj4eG6++Wb+/e9/M3z48NL7QYhI1HN37n5jMa9/nc51P27D6FPb0rh2QpF1jtqiMLOJhD7EO5lZupmNAp4mNIn7R2a2wMyeCwJYSmiO3mXA+8AtwcT1ADcBLxDq4F7L9/0aLwINzGwNoQnt7wq2lQncB8wNHn8NymJWixYtGDhwIABXX301M2bMoE2bNnTs2BGAESNGMHPmTOrUqUNCQgLXXXcdb7zxBjVr1jzmfTRq1Ii2bdvy1VdfsX37dlauXMnAgQOZMWMG8+bNo2/fvvTq1YsZM2awbt26MjlOEYk+n6zcyvXjUxj6/Fe8/nU6vzqzI/de2OWoSQKOoUXh7lcUUvxiEes/ADxQSHkK0K2Q8izgsiNs6yW+n+C+VBzLN/+ycqxXFMXFxTFnzhxmzJjBpEmTePrpp/n444+PeT9Dhw5lypQpdO7cmUsuuQQzw90ZMWIEDz74YHHDF5EY5O6M/TKV+95dRuPaCTSsXY2bB7XjtjPaH/M2NIRHOdqwYQOzZs0CYOLEiZx55pmkpqayZs0aACZMmMCpp57K3r172bVrF+effz5PPPEECxYs+MG2ateuzZ49ewrdz6WXXspbb73FxIkTGTp0KABnnHEGU6dOZevWrQBkZmbyzTdHHCxSRCoAd+fxj1bxl3eWceZJTfj4jlN599afcOe5nY/rUvhKMYRHtDjppJMYN24cN9xwAx06dODJJ5+kf//+XHbZZeTm5tK3b19uvPFGMjMzGTx4MFlZWaFf9OOP/2Bb11xzDTfeeCM1atT4LvkUqFevHl26dGHZsmX069cPgC5dunD//fdz9tlnk5+fT3x8PP/6179o1apVuRy7iJSv/Hzn0Q9X8swnaxma3IIHL+1OlSrFu0+qws2ZnZyc7IfPR7F8+XJOOumkCEUUkpqayoUXXsiSJUsiGsexiIafl4gUX8aebH7z2kJmrspgWN8W/O2SoycJM5vn7smFLVOLQkSkgsjPdyanpPHQ9BVk5eTxwCXduLJfyxKPuKBEUU5at24dE60JEYlN2/dm88vJC/hs9TZOaVOfBy7pRvvGtUtl20oUIiIxLCsnj9fmpfPUjNXsPJBTaq2IcJUmUbi7Brw7BhWtz0qkItux7yCXPPMFqdv307tlXcZe248uzeqU+n4qRaJISEhg+/btGmr8KArmo0hIOPoNOCISWfn5zq+nLGDTzixevrYvgzo2KrPPt0qRKJKSkkhPT6cijyxbWgpmuBOR6PTpqgzGf5nKjv0H+XrDTv46uCundWpcpvusFIkiPj5eM7aJSEzLz3ee+ngNT8xYRdM6CTSqk8BNg9rx8/5lfy9UpUgUIiKxbNf+HH41ZQEfr9jKpb2b88Al3alRrerRK5YSJQoRkSi24tvd3DBhHpt2HuC+wV25un+rcu9rVaIQEYlS0xdv5tdTFlI7IY5Jo/vTp1X9iMShRCEiEmXcnednruOh6Ss4uWVdnru6D43rRO5qRCUKEZEokpuXzx+nLeXV2Ru4sEdTHr2sJwnx5dcfURglChGRKLFz/0Fun7SAT1dlcNOgdvz27E7FHvG1NClRiIhEgTnrM7l90ny27c3mwUu7c0W/lpEO6TtKFCIiEbZ6yx6GvzSbE+sk8MZNA+melBjpkA6hRCEiEkEHDuZxy6tfU6t6HFNuGBDRTusjUaIQESln32zfx9Mfr+E/izdzMC+fnDxn/Mh+UZkkQIlCRKTcZOXk8fTHa3ju07VUrWIM7tWMeidUo3eLevy0Y6NIh3dEShQiIuXgy7XbuOfNJazfto+fnZzE787tFLUtiMMpUYiIlKEFaTt5/tO1TF/yLS3r1+SVUafw4w4NIx3WcVGiEBEpA5+uyuCZ/61h9vpMaifEcevp7bl5UPtyHcyvtFQ52gpm9pKZbTWzJWFl9c3sIzNbHTzXC1t2t5mtMbOVZnZOWHkfM1scLPunBaNamVl1M5sclM82s9ZhdUYE+1htZiNK66BFRMrSC5+tY8RLc0jL3M+9F5zErLvP4Ddnd4rJJAHHkCiAscC5h5XdBcxw9w7AjOA9ZtYFGAZ0Deo8Y2YFP5lngdFAh+BRsM1RwA53bw88DjwcbKs+8CfgFKAf8KfwhCQiEm3cnec+Xcv97y3ngu5N+eS3p3HdT9pSq3psn7w5aqJw95lA5mHFg4FxwetxwMVh5ZPcPdvd1wNrgH5m1hSo4+6zPDQp8/jD6hRsaypwRtDaOAf4yN0z3X0H8BE/TFgiIlHhYG4+d72+mIemr+CCHk15YlgvqsUdy3fx6FfcNNfE3TcDuPtmMyuYh6858FXYeulBWU7w+vDygjppwbZyzWwX0CC8vJA6hzCz0YRaK7RsGT23vYtI5ZCTl8+Nr8zj4xVb+cVp7fn1WR2jYoym0lLa6a6wn4wXUV7cOocWuo9x92R3T27UKHqvRRaRiic/37nr9cV8vGIr91/cjTvOiY6B/EpTcVsUW8ysadCaaApsDcrTgRZh6yUBm4LypELKw+ukm1kckEjoVFc6MOiwOp8UM14RkVLl7kydl85zn65lbcY+fnlmB64uh/mrI6G4LYppQMFVSCOAt8PKhwVXMrUh1Gk9JzhNtcfM+gf9D8MPq1OwrSHAx0E/xgfA2WZWL+jEPjsoExGJqIO5+fzmtYX8duoiEuKr8uSwXtx+RodIh1VmjtqiMLOJhL7ZNzSzdEJXIj0ETDGzUcAG4DIAd19qZlOAZUAucIu75wWbuonQFVQ1gOnBA+BFYIKZrSHUkhgWbCvTzO4D5gbr/dXdD+9UFxEpV1+t286D/1nOwvRd/Pqsjtx6evtyn8O6vFnoy3vFkZyc7CkpKZEOQ0QqmCUbd/HIByuZuSqDJnWq84cLu3Bhj2aRDqvUmNk8d08ubFlsX9wrIlLG9h/M5aHpKxg/6xvq1YznnvNP4ucDWkV8etLypEQhInIESzbu4taJ80ndvo+RA9vwq7M6UDshPtJhlTslChGRw7g742d9wwPvLadBrWq8el1/BrRrEOmwIkaJQkQkzK4DOdz1+iKmL/mW0zs35tHLelL/hGqRDiuilChERAKfrc7gzqmLyNiTze/P78x1P25b4W6eKw4lChGp9HLy8vn7BysZM3Md7RqdwOs3/YieLepGOqyooUQhIpXa0k27uPetJczfsJOf92/FPRecVKmuaDoWShQiUum4O/PTdvLKrG94c8FG6taI56krenNRz4pzX0RpUqIQkUpl5/6D3Dl1ER8u28IJ1aoycmAbbju9A4k1K99lr8dKiUJEKo056zO5fdJ8tu3N5s5zOzF8QOuYn1SoPOgnJCIVXl6+86//reGJ/66iRf2avHHTQLonJUY6rJihRCEiFdr2vdn8cvICPlu9jcG9mnH/xd0q5d3VJaFEISIVVvqO/Vz1wmw278riwUu7M6xviwo/0mtZUKIQkQppzda9DH9xNnuzc5k0uj8nt6wX6ZBilhKFiFQo7s6b8zdy71tLqBFflVev70+35uqPKAklChGpMD5fvY1HP1zJgrSd9Gtdnyev6EXTxBqRDivmKVGISIWwKH0nI16eQ9PEBB64pBtDk1sQV7W4sz1LOCUKEYl52bl53PHaQhrVqs57t/2ExBq6qqk0KVGISEwq6It4e8Emtu3NZtWWvbx8bV8liTKgRCEiMeeb7fu4580lfL5mG20ankC9mvH89pxOnNapcaRDq5CUKEQkZuTk5fPi5+t54r+riKtShfsGd+WqU1ppzogypkQhIjFhYdpO7npjMcs37+asLk346+CuuqKpnChRiEhU252Vw2MfrmLcrFQa167Oc1f34dxuJ0Y6rEpFiUJEotLurBxe/jyVFz9fx57sXH7evxV3nNOJOhqnqdyVKFGY2a+A6wAHFgPXAjWByUBrIBW43N13BOvfDYwC8oDb3P2DoLwPMBaoAfwHuN3d3cyqA+OBPsB2YKi7p5YkZhGJbjl5+YyZuY7nP13L7qxczu7ShNvO6KC7qyOo2InCzJoDtwFd3P2AmU0BhgFdgBnu/pCZ3QXcBfzOzLoEy7sCzYD/mllHd88DngVGA18RShTnAtMJJZUd7t7ezIYBDwNDixuziES31G37+OXkBSxI28mZJzXhl2cqQUSDkt62GAfUMLM4Qi2JTcBgYFywfBxwcfB6MDDJ3bPdfT2wBuhnZk2BOu4+y92dUAsivE7BtqYCZ5iGfhSpcNydKSlpnP/Pz1iXsZenr+zNCyOSlSSiRLFbFO6+0cweBTYAB4AP3f1DM2vi7puDdTabWcGFzc0JtRgKpAdlOcHrw8sL6qQF28o1s11AA2BbeCxmNppQi4SWLVsW95BEJAJy8vL5w1tLmDQ3jf5t6/PY5b1oVldXM0WTYrcozKweoW/8bQidSjrBzK4uqkohZV5EeVF1Di1wH+Puye6e3KhRo6IDF5Gosf9gLteNS2HS3DR+cVp7/n1dfyWJKFSSzuwzgfXungFgZm8APwK2mFnToDXRFNgarJ8OtAirn0ToVFV68Prw8vA66cHprUQgswQxi0iU2H8wl5Fj5zJnfSYPXtqdK/rpbEC0KkkfxQagv5nVDPoNzgCWA9OAEcE6I4C3g9fTgGFmVt3M2gAdgDnBaao9ZtY/2M7ww+oUbGsI8HHQjyEiMWxtxl6u/L/ZzFmfyeNDeylJRLmS9FHMNrOpwNdALjAfGAPUAqaY2ShCyeSyYP2lwZVRy4L1bwmueAK4ie8vj50ePABeBCaY2RpCLYlhxY1XRCJv5bd7mDhnA6/O3kBCfBX+deXJnNe9aaTDkqOwivYFPTk52VNSUiIdhogEcvPyeXP+Rl6ds4H5G3YSX9W4qGcz7jqvM41rJ0Q6PAmY2Tx3Ty5sme7MFpEyszBtJ/e8tZglG3fTrtEJ3HvBSVzSuzkNalWPdGhyHJQoRKTU7TqQw2MfrmT8V9/QqFZ1nr6yNxd0b4pug4pNShQiUmp2Z+XwfzPXMfbLVPZl5zJiQGt+c3ZHamt8ppimRCEiJebuvLVgIw+8t4Jte7M5r9uJ3Hp6B7o0qxPp0KQUKFGISLG5O0s37ea+d5cxe30mPVvU5eVr+tI9SUNvVCRKFCJyXAqSw/Qlm5m+5FvWZewjsUY8f7ukO8P6ttBscxWQEoWIHJM9WTlMSUln8twNrNqyl6pVjP5t63PtwDZc2L0p9U6oFukQpYwoUYhIkfZl5zJm5jpe/mI9u7Ny6d2yLg9c0o3zujWlvpJDpaBEISKF2rB9P+8s2sTYL1PJ2JPNOV2bcPOg9vRsUTfSoUk5U6IQke/k5TtzUzMZM3MdH68IjefZv219xvy8D71b1otwdBIpShQilVzGnmxmrsrg01UZfLY6gx37c6hXM55fndmRIclJNNew35WeEoVIJbRm6x7eW/QtHy3/liUbdwPQsFY1TuvcmFM7NuKsLk2oWU0fDxKivwSRSmJ3Vg4fLt3C5LkbmJu6AzM4uWU9fntOJ07t2IguTevo0lYplBKFSAWUn+98k7mfRek7WZi2i8UbQ88H8/Jp1aAmvz+/Mxf3ak7jOhq9VY5OiUKkgnB3vlqXyYSvUvls9Tb2ZOUCkBBfha7NEhnxo1ac170pvVvU1eB8clyUKERiWFrmfr5Ys40v1m5n1tptbNt7kLo147mwRzN6tUikR1JdOjSuRVzVkkxmKZWdEoVIlMvLd9Zl7GX5t3vYn51L5v6DrPp2D/M27CAt8wAAjWpXZ2D7hpzasRHnd29KQnzVCEctFYkShUgUcXfSMg+wMH0ni9J3sih9F0s27mLfwbxD1muamEC35omMGtiGge0b0r5xLZ1OkjKjRCFSDnbuP8jijbv4dlcWuw7ksHN/Dtv3ZZOWeYDt+w6SnZPHgZw89mblsic71LdQLa4KXZrWYUifJLon1aVb8zrUSYinVkIcdTS/g5QjJQqRUuburM3Yy6x1mcxdn8mi9J2kbt9/yDpVDOqfUI3m9WrSvG4NalSrSkJcFWpWq0rHE2vTM6kuHZvUplqc+hYk8pQoREoolBj2MWvddr5at53Z67azbe9BAJrUqU7vFvW4vG8LeibVpWX9miTWjKdWtTjdsyAxQ4lC5Di5O0s27uZ/K7eyKH0nC9J2sW1vNgAn1kngJx0a0b9tffq3bUDL+jXVdyAxT4lC5Bi4O3NTd/DOwk38d/kWNu/KwgzaN6rFTzs0pF+b+gxop8QgFZMShUgRtu3N5o2v05k0N411GftIiK/CTzs04tdndeSMk5poPgapFEqUKMysLvAC0A1wYCSwEpgMtAZSgcvdfUew/t3AKCAPuM3dPwjK+wBjgRrAf4Db3d3NrDowHugDbAeGuntqSWIWKYy7s3rrXpZv3s2Kb/ewMnhs3Bm6T6FPq3o8MqQdF/ZoqsHypNIp6V/8k8D77j7EzKoBNYHfAzPc/SEzuwu4C/idmXUBhgFdgWbAf82so7vnAc8Co4GvCCWKc4HphJLKDndvb2bDgIeBoSWMWeQ7WTl5vLNwEy99kcryzaFRVOOqGO0a1aJPq3pceUpLzurShI5Nakc4UpHIKXaiMLM6wE+BawDc/SBw0MwGA4OC1cYBnwC/AwYDk9w9G1hvZmuAfmaWCtRx91nBdscDFxNKFIOBPwfbmgo8bWbm7l7cuEUA9h/M5aXP1zP2y1S27T1Ixya1uO/ibvRtXY+2DWvpslSRMCVpUbQFMoCXzawnMA+4HWji7psB3H2zmTUO1m9OqMVQID0oywleH15eUCct2Faume0CGgDbwgMxs9GEWiS0bNmyBIcklcG0hZu4791lZOzJZlCnRlz347YMbN9AndAiR1CSRBEHnAzc6u6zzexJQqeZjqSw/0IvoryoOocWuI8BxgAkJyertSGFytx3kAfeW87rX6fTMymRZ686meTW9SMdlkjUK0miSAfS3X128H4qoUSxxcyaBq2JpsDWsPVbhNVPAjYF5UmFlIfXSTezOCARyCxBzFLJHDiYx7xvdvDpqq28OnsDB3LyuO309tx2RgeNqCpyjIqdKNz9WzNLM7NO7r4SOANYFjxGAA8Fz28HVaYBr5rZY4Q6szsAc9w9z8z2mFl/YDYwHHgqrM4IYBYwBPhY/RNyLL5cs42Xv0zls9UZZOXkU8XgvO5Nuf2MDuqYFjlOJb3q6Vbg38EVT+uAa4EqwBQzGwVsAC4DcPelZjaFUCLJBW4JrngCuInvL4+dHjwAXgQmBB3fmYSumhI5orTM/dz/3jI+WLqFxrWrMzS5BYM6Nya5VT1qayA9kWKxivYFPTk52VNSUiIdhpSznLx8nvp4Dc99upaqZtxyWjuu+0lbzcsgcozMbJ67Jxe2THcOSczbujuLX7w6nzmpmfy/ns24+/zONE2sEemwRCoMJQqJadm5eVz94mzSMg/w5LBeDO7V/OiVROS4KFFITHvmf2tZtWUvL1/Tl9M6Nz56BRE5bro+UGLW8s27eeaTNVzcq5mShEgZUqKQmLQ4fRdXvTCbujWr8ceLukY6HJEKTaeeJGZk5eQxY/lWZqzYwvtLvqVezWq8ct0pGupbpIwpUUjU23Ugh3/OWM3kuWnszc6lXs14zul6Ir87tzMnJiZEOjyRCk+JQqKWu/PG1xt5cPpyMvcd5OJezRmSnMQpbRpQVfNNi5QbJQqJSss27eZP05YwN3UHvVvWZey1/ejWPDHSYYlUSkoUElV2Hcjh8Y9WMX5WKnVrVuORn/VgSJ8kqqgFIRIxShQSNd5ZuIm/vLOM7fuyufqUVvzm7I7UramOapFIU6KQiNux7yAPTV/B5JQ0eiYl8vI1femepNNMItFCiUIiIicvnxnLtzBpbhqfr95Gnju/OK09vzxT80SIRBslCil3s9dt546pC0nLPECzxASu+0lbBvdqxklN60Q6NBEphBKFlBt3558z1vDEjFW0rF+TMT/vwxknNdGlriJRTolCykV+vvPnd5YyftY3XNq7Ofdd3FrV6mEAABZPSURBVI0TquvPTyQW6D9VylxOXj53vLaQtxds4oaftuWu8zpjplaESKxQopAytS87l9smzmfGiq3ceW4nbh7UPtIhichxUqKQMrM2Yy83TpjH2oy9PHBJN646pVWkQxKRYlCikDLx5Zpt3PDKPOKrVmHCqFMY2L5hpEMSkWJSopBSlZWTxytffcPD76+gTcMTeOmaviTVqxnpsESkBJQopFRk5eQxac4GnvlkLVv3ZHNap0Y8Maw3iTXiIx2aiJSQEoWU2LJNuxk1bi6bd2XRr019nhzWmwHtGkQ6LBEpJUoUUiJLNu7i6hdnUyO+Kq9efwo/aqe+CJGKpsSD6phZVTObb2bvBu/rm9lHZrY6eK4Xtu7dZrbGzFaa2Tlh5X3MbHGw7J8WXGRvZtXNbHJQPtvMWpc0Xik9C9J2cuX/fcUJ1eKYPHqAkoRIBVUao6/dDiwPe38XMMPdOwAzgveYWRdgGNAVOBd4xsyqBnWeBUYDHYLHuUH5KGCHu7cHHgceLoV4pRTM+2YHP39hNok145l8Q39aNlCHtUhFVaJEYWZJwAXAC2HFg4FxwetxwMVh5ZPcPdvd1wNrgH5m1hSo4+6z3N2B8YfVKdjWVOAM0y29ETdnfSbDX5xNw9rVmXLDAF3VJFLBlbRF8QRwJ5AfVtbE3TcDBM+Ng/LmQFrYeulBWfPg9eHlh9Rx91xgF/CDXlIzG21mKWaWkpGRUcJDkqJ8uXYbI16aw4mJCUwa3Z+miTUiHZKIlLFiJwozuxDY6u7zjrVKIWVeRHlRdQ4tcB/j7snuntyoUaNjDEeO12erM7j25bm0qF+DSaMH0KROQqRDEpFyUJKrngYC/8/MzgcSgDpm9gqwxcyauvvm4LTS1mD9dKBFWP0kYFNQnlRIeXiddDOLAxKBzBLELMX0v5VbuWHCPNo2PIF/X3cKDWpVj3RIIlJOit2icPe73T3J3VsT6qT+2N2vBqYBI4LVRgBvB6+nAcOCK5naEOq0nhOcntpjZv2D/ofhh9Up2NaQYB8/aFFI2XF3/rN4MzeMn0eHxrWYeH1/JQmRSqYs7qN4CJhiZqOADcBlAO6+1MymAMuAXOAWd88L6twEjAVqANODB8CLwAQzW0OoJTGsDOKVI/h6ww4emr6COesz6ZmUyPiRp5BYU3dai1Q2VtG+oCcnJ3tKSkqkw4hp2bl5/OPDVbzw2Toa1qrOLae1Z2jfFiTEVz16ZRGJSWY2z92TC1umO7PlENv3ZnPDhHmkfLODK/q15Pfnd6Z2gloRIpWZEoV8Z/WWPYwcN5etu7N5+sreXNijWaRDEpEooEQhQOjS15tf+Zrq8VWZfMMAerWoG+mQRCRKKFEIE+ds4N63ltChcS1evKYvzevqJjoR+Z4SRSWWn+/8/cOVPPvJWk7t2Ih/XXUytarrT0JEDqVPhUoqKyePO15byLuLNnNFv5bcN7grcVVLY4xIEalolCgqoaycPEZPmMfMVRncdV5nbvhpWzTWoogciRJFJbN1Txa3T1zAV+u388jPenB53xZHryQilZoSRSXy1vyN/GnaUg7k5PGPy3py6clJR68kIpWeEkUlkJWTx5+nLWXS3DT6tKrHwz/rQfvGtSIdlojECCWKCu7AwTyuHTuHr9Zlcstp7fjVmR3VaS0ix0WJogLLzs3jhlfmMXt9Jk8M7cXFvZsfvZKIyGH01bKCysnL5xevzmfmqgwevrSHkoSIFJsSRQWUm5fPLycv4KNlW/jr4K66sklESkSJooLJz3funLqI9xZt5p7zT2L4gNaRDklEYpwSRQXi7tzz1hLemL+RX5/Vket/2jbSIYlIBaBEUUG4O399dxkT52zg5kHtuPX09pEOSUQqCCWKCuKpj9fw8hepjBzYht+e00lDcohIqVGiqAA+WraFxz5axSW9m/OHC09SkhCRUqVEEeMWpO3k15MX0K15HR68tLuShIiUOiWKGDZt4SaGPj+LxJrxPHd1HxLiq0Y6JBGpgHRndowa+8V6/vzOMvq1rs+zV59Mg1rVIx2SiFRQShQx6NlP1vLw+ys4p2sT/nlFb6rHqSUhImVHiSLGTJyzgYffX8FFPZvx+OU9NcCfiJS5Yn/KmFkLM/ufmS03s6VmdntQXt/MPjKz1cFzvbA6d5vZGjNbaWbnhJX3MbPFwbJ/WtAja2bVzWxyUD7bzFoX/1Bj3wdLv+WeNxczqFMjHlOSEJFyUpJPmlzgN+5+EtAfuMXMugB3ATPcvQMwI3hPsGwY0BU4F3jGzArOmTwLjAY6BI9zg/JRwA53bw88Djxcgnhj2ux127l14ny6J9XlmatOJl5JQkTKSbE/bdx9s7t/HbzeAywHmgODgXHBauOAi4PXg4FJ7p7t7uuBNUA/M2sK1HH3We7uwPjD6hRsaypwhlXC6z9Tt+3juvEpJNWrwcvX9KVmNZ0xFJHyUypfS4NTQr2B2UATd98MoWQCNA5Waw6khVVLD8qaB68PLz+kjrvnAruABoXsf7SZpZhZSkZGRmkcUtTIz3d+O3UhAOOu7Uf9E6pFOCIRqWxKnCjMrBbwOvBLd99d1KqFlHkR5UXVObTAfYy7J7t7cqNGjY4WckwZNyuVuak7+NNFXWlRv2akwxGRSqhEicLM4gkliX+7+xtB8ZbgdBLB89agPB0InxghCdgUlCcVUn5IHTOLAxKBzJLEHEuWbNzFw++v4PTOjfnZyZp4SEQioyRXPRnwIrDc3R8LWzQNGBG8HgG8HVY+LLiSqQ2hTus5wempPWbWP9jm8MPqFGxrCPBx0I9R4W3dk8X141OoX7MajwzpoaE5RCRiStIrOhD4ObDYzBYEZb8HHgKmmNkoYANwGYC7LzWzKcAyQldM3eLueUG9m4CxQA1gevCAUCKaYGZrCLUkhpUg3piRn+/cNnE+O/fnMPWmATTUXdciEkHFThTu/jmF9yEAnHGEOg8ADxRSngJ0K6Q8iyDRVCbjZ6Xy1bpMHvlZD7o2S4x0OCJSyeli/CizLmMvD72/gtM6NeKy5KSjVxARKWNKFFFkUfpOLn/+KxLiq/LgpeqXEJHooEQRJb5cs42hz39FQnwVpt44gBMTEyIdkogIoEEBo8KXa7cxctxcWtavySvXnULj2koSIhI9lCgibNWWPVw/LoUW9Wry6vX9dYWTiEQdnXqKoF0Hchg9PoUa1eKYMOoUJQkRiUpKFBFyMDefX7z6NRt3HuC5q09Wn4SIRC2deoqAgoH+Plu9jUeG9CC5df1IhyQickRqUZQzd+f+95bz9oJN3HluJy5PbnH0SiIiEaREUc6e/XQtL32xnmsHtuamU9tFOhwRkaNSoihHU1LSeOT9lQzu1Yw/XNBFN9SJSExQoign/122hbvfWMxPOjTk70N6UqWKkoSIxAYlinIwNzWTW179mm7N6vDc1X2oFqcfu4jEDn1ilbEV3+5m1Ni5NK9bg5eu6csJ1XWhmYjEFiWKMvSfxZsZ8uwsEuKrMm5kPxrohjoRiUH6elsGDubm87f/LGfsl6n0alGXp6/sTVI9zXctIrFJiaKUbd2dxegJ81iQtpORA9tw13md1SchIjFNiaIUrc3Yy4iX5pC57yDPXX0y53ZrGumQRERKTImilLy9YCP3vrWE6nFVmDx6AN2TNIWpiFQMShQllJWTx71vLWHqvHT6tKrHE0N70aK++iNEpOJQoiiBzH0HGTVuLvM37OS209tz2xkdiKuq/ggRqViUKErgj28vYenG3Txz1cmc3139ESJSMenrbzHNXJXBu4s2c/Np7ZQkRKRCU6IohkXpO/nD20to3aAmN2oEWBGp4GLi1JOZnQs8CVQFXnD3h8pz/3uzc5m7PpNZ67bz+eptLNu8mxOqVeXFa/qSEF+1PEMRESl3UZ8ozKwq8C/gLCAdmGtm09x9WWnvy93Ztvcgm3YeYMW3u1mQtouFaTtZuWUPeflOtapV6NWyLn+6qAtD+iRROyG+tEMQEYk6UZ8ogH7AGndfB2Bmk4DBQKGJYtvebJ7/dC25+U5OXj55+U5OnpObl/+DsgM5uew+kMuerBx2Hshh864sDubmf7etOglx9GxRl5tPakf/tg04uWU9alRTC0JEKpdYSBTNgbSw9+nAKeErmNloYDRAtRPb8+D0Fd8tq2IQV6UKcVWNuCpGfNUqVA2ea1SrSu2EOBJrVqNlgxM4t2sCTRMTaFq3Bu0b16JNgxM0b4SIVHqxkCgK+6T2Q964jwHGAPQ+uY9//pdziKtqxFepog96EZESioVEkQ60CHufBGw60spVq5jmfBARKUWxcHnsXKCDmbUxs2rAMGBahGMSEak0ov6rt7vnmtkvgA8IXR77krsvjXBYIiKVRtQnCgB3/w/wn0jHISJSGcXCqScREYkgJQoRESmSEoWIiBRJiUJERIpk7n70tWKIme0BVgKJwK7jqHq86xenTknXbwhsi2A8pb2Pwo6nvH+mZVHnaL+nkm6/vI/5WI4n2o+hMIcfV7QdQ3GPOZ5j//sL30cnd69d6FruXqEeQErwPOY46x3X+uWxj8PXLzi2SMVT2vso7HjK+2daRvso8vcUa8d8LMcT7cdwLL+naDuG4h7z8fz9HevvuSKfenqnjNcvj31E2/rlsY9oW7+4dcpy+9F4zNG2fnnsI9rWL7N9VMRTTynunhzpOMpCRTu2inY8BSracVW04ymg4zr2ehWxRTEm0gGUoYp2bBXteApUtOOqaMdTQMd1jPUqXItCRERKV0VsUYiISClSohARkSLFXKIwMzezCWHv48wsw8zejWRcpcXMLgmOsXOkYymJiv57AjCzvZGOoSwc7bjM7BMzi4lO4Iry/3Q4M7vHzJaa2SIzW2Bmpxy9VvHFXKIA9gHdzKxG8P4sYOPxbMDMonnU3CuAzwnNu3HMzCzaJvMu8e9JpBQU6/8pmpnZAOBC4GR37wGcyaHTRZe6WEwUANOBC4LXVwATCxaYWT8z+9LM5gfPnYLya8zsNTN7B/iw/EM+OjOrBQwERhH8YZvZIDObaWZvmtkyM3vOzKoEy/aa2V/NbDYwIHKRH1Fxfk+fmVmvsPW+MLMe5Rr1cQh+P++GvX/azK4JXqea2V/M7GszWxxL32qLOq5YUcT/05F+X+eb2Qoz+9zM/hnFrd+mwDZ3zwZw923uvsnM+pjZp2Y2z8w+MLOm8F0L8Ing/2yJmfU73h3GaqKYBAwzswSgBzA7bNkK4Kfu3hv4I/C3sGUDgBHufnq5RXp8Lgbed/dVQKaZnRyU9wN+A3QH2gGXBuUnAEvc/RR3/7zcoz264vyeXgCuATCzjkB1d19UbhGXvm3ufjLwLHBHpIOpZI70//QDwd/o88B57v5joFE5xVgcHwItzGyVmT1jZqeaWTzwFDDE3fsALwEPhNU5wd1/BNwcLDsuMZkogg+O1oS+pR4+oVEi8JqZLQEeB7qGLfvI3TPLJcjiuYLQhyvB8xXB6znuvs7d8wh9K/9xUJ4HvF6+IR67Yv6eXgMuDP7wRwJjyyXYsvNG8DyP0M9Cys+R/p8K0xlY5+7rg/cTi1g3otx9L9AHGA1kAJOBG4BuwEdmtgC4F0gKqzYxqDsTqGNmdY9nn9F8rv5opgGPAoOABmHl9wH/c/dLzKw18EnYsn3lFNtxM7MGwOmEzus7oWlfndAH7OE3uxS8zwqSRzQ7rt+Tu+83s4+AwcDlQLR3muZy6BeuhMOWZwfPecTW/9vRjiuqFfH/NI3Cj8vKN8KSCf7vPwE+MbPFwC3AUnc/0inoI32GHJOYbFEEXgL+6u6LDytP5PtO02vKNaKSGQKMd/dW7t7a3VsA6wm1HvqZWZugb2Iooc65WFGc39MLwD+BuVHeAgT4BuhiZtXNLBE4I9IBlZJYP64j/T9B4ce1AmgbfGmB0P9ZVDKzTmbWIayoF7AcaBR0dGNm8WYWfjZlaFD+Y2CXux/XqLSx9A3nEO6eDjxZyKJHgHFm9mvg4/KNqkSuAB46rOx14CZgVrCsOzATeLN8Qyu+4vye3H2eme0GXi6HEIsluHIu293TzGwKsAhYDcyPbGQlU4GO60j/T1cCPzgudz9gZjcD75vZNmBOOcZ6vGoBTwWnj3KBNYROQ40B/hkkwDjgCWBpUGeHmX0J1CF0Sve4aAiPKGdmg4A73P3CSMdSXsysGaFmdWd3z49wOIUys57A/7n7cV9BEs0q6nEdCzOr5e57zcyAfwGr3f3xSMdVUmb2CaHPkJTibiOWTz1JBWRmwwldHXVPFCeJGwl1Dt4b6VhKU0U9ruNwfdARvJTQqdHnIxxP1FCLQkREihT1LQoza2Fm/zOz5Ra6Zf32oLy+mX1kZquD53pB+VnBDSeLg+fTw7bVJyhfE9xQE1NXOoiIRELUJwpCnTW/cfeTgP7ALWbWBbgLmOHuHYAZwXsIzRV7kbt3B0YAE8K29SyhTp8OwePc8jkEEZHYFfWJwt03u/vXwes9hC4Da07oOvtxwWrjCN2FibvPd/dNQflSICG4FK4pUMfdZ3nofNv4gjoiInJkUZ8owgXXOPcm1NnZxN03QyiZAI0LqfIzYH4wJkpzID1sWXpQJiIiRYiZ+yiCAb5eB37p7ruP1r0Q3GzyMHB2QVEhq6knX0TkKGKiRRGM+/M68G93Lxg7Z0vY6IhNga1h6ycRuiltuLuvDYrTOXTskyRgEyIiUqSoTxTBlUkvAsvd/bGwRdMIdVYTPL8drF8XeA+4292/KFg5OD21x8z6B9scXlBHRESOLOrvowjGJvkMWAwU3ID1e0L9FFOAlsAG4DJ3zzSze4G7Cd2eX+Bsd99qoVm5xgI1CM2VcKtH+w9ARCTCoj5RiIhIZEX9qScREYksJQoRESmSEoWIiBRJiUJERIqkRCEiIkVSohA5Bmb2ZzO7o4jlFweDVR5tO4esZ2Z/NbMzSytOkbKgRCFSOi4GjpooDl/P3f/o7v8ts6hESoEShcgRmNk9ZrbSzP4LdArKrjezuWa20MxeN7OaZvYj4P8BfzezBWbWLni8H8yJ8pmZdT7CemPNbEiw7VQz+5uZzTKzFDM72cw+MLO1wexzBXH9NohhkZn9JQI/GqlkYmZQQJHyZGZ9gGGERiuOA74G5gFvuPv/BevcD4xy96fMbBrwrrtPDZbNAG5099VmdgrwjLufXsh6h+86zd0HmNnjhEYRGAgkEBoy/zkzO5vQXCr9CA10Oc3MfuruM8vshyGVnhKFSOF+Arzp7vsBgg94gG5BgqgL1AI+OLxiMNLxj4DXwhJB9WPcb8F+FgO1gjlY9phZVjCO2dnBY36wXi1CiUOJQsqMEoXIkRU2vs1Y4GJ3X2hm1wCDClmnCrDT3XsVY5/ZwXN+2OuC93GEWhEPuvvzxdi2SLGoj0KkcDOBS8yshpnVBi4KymsDm4Oh768KW39PsAx33w2sN7PLIDQCspn1PHy9YvoAGBm0WjCz5mZW2KRdIqVGiUKkEMH0u5OBBYTmQvksWPQHQiMXfwSsCKsyCfitmc03s3aEksgoM1tIqH9h8BHWO964PgReBWaZ2WJgKiVLPCJHpdFjRUSkSGpRiIhIkZQoRESkSEoUIiJSJCUKEREpkhKFiIgUSYlCRESKpEQhIiJFUqIQEZEi/X8m7OYmsFgcxgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEiCAYAAADQ05jiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3yV5f3/8dcHCBuZQYGwNyJECENxi7buURStihQVcdX+qraODuvX1jrqaJ1UBQUBB2CtdeMALSvsvQMElBVGGAkk+fz+OHdsjCE5gZyR5P18PM4jJ/d9n/v+nMPhnetc57qv29wdERGJX1ViXYCIiBRPQS0iEucU1CIicU5BLSIS5xTUIiJxTkEtIhLnFNRS7pnZi2b2+2LW329mL0ezJpGyZBpHLcUxMwc6uvvqAsseBDq4+7XB7/cDNwGJwC7gG3cfHINyMbMzgLHunhSL44tEglrUclTM7HrgOmCgu9cFUoApsa1KpGJRUMvR6gN87O5rANz9O3cfebiNzSzNzO4zs6VmttPMRplZzQLrbzKz1WaWYWbvmVnzYLmZ2VNmttXMdpvZQjPrHqwbbWYPm1kd4EOguZntDW7NzexBMxsbbPuRmd1eqKYFZnZ5cL+LmX0aHH+FmV1ZzHNpFNS/OXgu7wbLG5rZ+2a2LVj+vpklFXjcUDNba2aZZrbOzK4psG6YmS0LHvexmbUu6flLxaeglqM1AxhiZveYWYqZVQ3jMdcAPwHaA52A3wGY2VnAI8CVQDNgPTAheMy5wGnB9g2AwcCOgjt1933AecBmd68b3DYXOvY44Or8X8ysG9Aa+E8Q9J8G2zQNtnvezI4/zPMYA9QGjg+2fypYXgUYFey3FXAAeDY4Xh3g78B57l4POBmYH6y7FLgfuJxQN9I0YHy4z18qrogFtZm9Gvz1Xxzm9lcGrawlZjYuUnVJ2XL3scAdhIL3K2Crmd1bwsOedfeN7p4B/Jn/Bec1wKvuPtfds4H7gJPMrA1wCKgHdCH03coyd//2CEqeDCTnt1SDY04KjnchkObuo9w9x93nAhOBQYV3YmbNCP1RGOHuO939kLt/FbwmO9x9orvvd/fM4DmeXuDheUB3M6vl7t+6+5Jg+c3AI8FzywH+UqDWsnr+Ug5FskU9GvhpOBuaWUdC/ykHuPvxwK8iWJeUTi6QUGhZAqHgAMDd33D3gYRaeiOAh8zsJ8Xsc2OB++uB5sH95sHv+fvdS6jV2MLdPyfUKn0O2GJmI83smNI+mSA4/wNcFSy6CngjuN8a6Gdmu/JvhIL8uCJ21RLIcPedhVeYWW0ze8nM1pvZHmAq0MDMqgat/sGEXqdvzew/ZtalwPGfKXDsDMDK8vlL+RSxoHb3qYTeaN8zs/ZBH+EcM5tW4A16E/Bc/pve3bdGqi4ptQ1Am0LL2lIgUPMFrcq3gYVAcf2nLQvcbwXkd09sJhRWwPfdBI2BTcH+/+7uvQl1NXQC7ili3+EMYxoPXG1mJwG1gC+C5RuBr9y9QYFbXXe/pYh9bAQamVmDItbdBXQG+rn7MYS6LCAUurj7x+5+DqHuneXAPwvs8+ZCx6/l7v8txfOXCijafdQjgTuCN9vdwPPB8k5AJzP7xsxmmFlYLXGJijeB35lZkplVMbOBwEXAO/D9F2MXmFm9YP15hIJkZjH7vC3YXyNCfbJvBsvHAb8ws2Qzq0Hoo/9Md08zsz5m1s/MEoB9QBah1n5hW4DGZla/mON/QOgPwkPAm+6eFyx/n9D78DozSwhufcysa+EdBN0OHxLqw24YbJsfyPUI9UvvCp7jH/MfZ2bHmtnFwR+hbGBvgefxInBffp+4mdU3syuC++E+f6mI3D1iN0ItscXB/bqE3rzzC9yWBeveJ9R3mECotZYONIhkbbqF/W9YC3gcSAN2A3OBiwusvxz4BtgJ7AEWAUOL2V8aoW6upYTGXL8G1C6wfgSwhtCnsfeBpGD52YRa6nuB7YS6K+oG60YDDxfYx6uEukx2EepOeZDQ2OqCdbxCqPXdp9DyzoS6RrYF+/gcSD7Mc2kU1L8leP6TguXNgS+DWlcS6nt2oBqhVvRXwWu5K9iuW4F9Xhe8hnsItbBfLen561bxbxE94SX4Euh9d+8e9KetcPdmRWz3IjDD3UcHv08B7nX32RErTmLCzNKAG939s1jXIlJeRK3rw933AOsKfJQzM+sZrH4XODNY3oRQV8jaaNUmIhLPIjk8bzwwHehsZulmdgOhb9BvMLMFwBLgkmDzj4EdZraU0Bc797i7xoiKiKC5PkRE4p7OTBQRiXMKahGROFctEjtt0qSJt2nTJhK7FhGpkObMmbPd3ROLWheRoG7Tpg2pqamR2LWISIVkZj862zefuj5EROKcglpEJM4pqEVE4lxE+qiLcujQIdLT08nKyorWIculmjVrkpSUREJC4ZlFRaSyilpQp6enU69ePdq0aYOZReuw5Yq7s2PHDtLT02nbtm2syxGROBG1ro+srCwaN26skC6GmdG4cWN96hCRH4hqH7VCumR6jUQqn5Km8tCXiaXw4osv8vrrrwMwevRoNm/+33VTb7zxRpYuXRqr0kSknNqWmc0tY+cWu03U+qgrghEjRnx/f/To0XTv3p3mzUOX+3v55ZdjVZaIlEO79x/i+S9XM3bGeg7lqkUNQFpaGl26dOH666+nR48eDBo0iP379zNlyhROPPFETjjhBIYNG0Z2djYA9957L926daNHjx7cfffdADz44IM88cQTvPPOO6SmpnLNNdeQnJzMgQMHOOOMM0hNTeWFF17gN7/5zffHHT16NHfccQcAY8eOpW/fviQnJ3PzzTeTm6srKYlURv9esJmzn/ySf05by8Bux/LBnacWu31MWtR/+vcSlm7eU6b77Nb8GP540fHFbrNixQpeeeUVBgwYwLBhw3jyySd56aWXmDJlCp06dWLIkCG88MILDBkyhMmTJ7N8+XLMjF27dv1gP4MGDeLZZ5/liSeeICUl5UfrTjrpJB577DEA3nzzTR544AGWLVvGm2++yTfffENCQgK33norb7zxBkOGDCnT10FE4lfWoVzun7SISfM20bNlA14b1pfjmxd3ec+QStOiBmjZsiUDBgwA4Nprr2XKlCm0bduWTp06AXD99dczdepUjjnmGGrWrMmNN97IpEmTqF27dtjHSExMpF27dsyYMYMdO3awYsUKBgwYwJQpU5gzZw59+vQhOTmZKVOmsHatLmIjUlnsP5jDsNGzmTx/E78a2JGJI04KK6QhRi3qklq+kRLuiIpq1aoxa9YspkyZwoQJE3j22Wf5/PPPwz7O4MGDeeutt+jSpQuXXXYZZoa7c/311/PII48cafkiUk5l5+Ry0+upzFi7g79d0ZPLeyWV6vGVqkW9YcMGpk+fDsD48eMZOHAgaWlprF69GoAxY8Zw+umns3fvXnbv3s3555/P008/zfz583+0r3r16pGZmVnkcS6//HLeffddxo8fz+DBgwE4++yzeeedd9i6dSsAGRkZrF9/2MmyRKSCyM1z7n57Id+s3sHjg0of0hBmi9rMGgAvA90JXfZ+mLtPL/XRYqxr16689tpr3HzzzXTs2JFnnnmG/v37c8UVV5CTk0OfPn0YMWIEGRkZXHLJJWRlZeHuPPXUUz/a19ChQxkxYgS1atX6PvzzNWzYkG7durF06VL69u0LQLdu3Xj44Yc599xzycvLIyEhgeeee47WrVtH5bmLSPQdys3j128t4N8LNnPfeV34We/ShzSEec1EM3sNmObuL5tZdaC2u+863PYpKSleeD7qZcuW0bVr1yMqsiykpaVx4YUXsnjx4pjVEK5Yv1YicvTWbtvLb95ZSOr6ndx7XhdGnN6+2O3NbI67pxS1rsQWtZkdA5wGDAVw94PAwVJXLSJSSUxduY0RY+eQULUKTw9O5tITWxzV/sLp+mgHbANGmVlPYA5wp7vvO6ojR1mbNm3KRWtaRMqvvDxnzIz1PPyfpXRoWo9RQ/twXP2aR73fcL5MrAb0Al5w9xOBfcC9hTcys+Fmlmpmqdu2bTvqwkREypPd+w9x/ahZ/PG9JQzo0IQJw/uXSUhDeEGdDqS7+8zg93cIBfcPuPtId09x95TExCKvz1jixCOi10ikPErfuZ/LX/iGGWt38OfLujNqaB/q1yq7OeVLDGp3/w7YaGadg0VnA6WefahmzZrs2LFDQVSM/Pmoa9Ysm7/CIhJ5adv3ceWL09mWmc3YG/pxTb/WZT4LZrgnvNwBvBGM+FgL/KK0B0pKSiI9PR11ixQv/wovIhL/Vm3J5JqXZ5KT54wf3j/sMw1LK6ygdvf5QJHDRsKVkJCgq5aISIWxdPMerntlJlWqGG8O70/HY+tF7Fia5lREpJQWbNzFkFdnUbt6Vcbd1J+2TepE9HgKahGRUkhNy+AXo2bToE4C427sT8tG4U/adqQq1VwfIiJH498LNnPdK7NIrFeDt24+KSohDWpRi4iU6FBuHo98sJxXv1lH79YNeeHaXjStF73RWQpqEZFibN2TxW3j5jI7bSdDT27D/ed3pXq16HZGKKhFRA5j/sZd3PR6KnuzcnjmqmQuST66OTuOlIJaRKQIqWkZDB01m4Z1Ehhzw8l0Oe6YmNWioBYRKWTG2h0MGz2b446pybibym7OjiOlUR8iIgV8s3o7Q0fNonmDWmU6sdLRUItaRCTw1cptDH89lbZN6jD2xn40qVsj1iUBCmoREQCmLNvCLWPn0qFpXcbe2I9GdarHuqTvqetDRCq9f83fxIixc+h8XD3G3RRfIQ1qUYtIJebuPP/lGh7/eAV92zbin0NSynQe6bKioBaRSulQbh6/f3cxE2Zv5JLk5jw2qAc1qlWNdVlFUlCLSKWzNzuHW9+Yy9SV27j9zA7cdW6nMp/svywpqEWkUsnOyWXY6NnMWb+Tv15+Alf1bRXrkkqkoBaRSsPduX/SYmaty4jpKeGlpVEfIlJpvPjVWibOTefOszuWm5AGBbWIVBIfLf6Oxz5ezkU9m/OrgR1jXU6pKKhFpMJbvGk3/+/N+fRMasDjg3rE9ReHRVFQi0iF9t3uLG54bTYNaycwckhvaibE5xC84ujLRBGpsA4czP1+Pul3bjk5qldlKUsKahGpkPLynF+/NZ/Fm3fzz+tS6NosdvNJH62wgtrM0oBMIBfIcfeUSBYlInI03J1HP17Oh4u/43cXdGVgt2NjXdJRKU2L+kx33x6xSkREykBenvOnfy/htenr+Xm/VtxwSttYl3TU1PUhIhWG+/9C+qZT23LfeV3L3QiPooQ76sOBT8xsjpkNL2oDMxtuZqlmlrpt27ayq1BEJEwvTV3La9PXc8MpbXnggm5UqVL+QxrCD+oB7t4LOA+4zcxOK7yBu4909xR3T0lMTCzTIkVESvLuvE389cPlXNijGQ+c3zXW5ZSpsILa3TcHP7cCk4G+kSxKRKQ05qzP4J53FtC/XSP+dmXPCtOSzldiUJtZHTOrl38fOBdYHOnCRETCkbHvILe9MY9m9Wvx0nUpcTun9NEI58vEY4HJQYd8NWCcu38U0apERMJwMCeP28fNJWPfQSbdenJcXp2lLJQY1O6+FugZhVpERMLm7vx24kL+u2YHf7uiJ91b1I91SRGjuT5EpFx64pMVTJ63ibvO6cTPeifFupyIUlCLSLkzbuYGnvtiDVf1acntZ3WIdTkRp6AWkXLl8+Vb+N27izijcyIPX9q9QpzQUhIFtYiUG7PTMrjtjXl0a34Mz/28F9WqVo4IqxzPUkTKvQUbd/GLUbNpVr8mo4b2pU6NyjMDhoJaROLe0s17GPLqLBrWSeCNm/qRWK9GrEuKKgW1iMS1rZlZXD9qFrWrV2Xcjf1pVr9WrEuKusrz2UFEyp3cPOeX4+eRmXWId28bQMtGtWNdUkwoqEUkbj3/xWpmrM3gb1f0pMtx5fcKLUdLXR8iEpcWb9rNM1NWcXHP5hX+hJaSKKhFJO5k5+Ry11sLaFSnOg9dcnysy4k5dX2ISNx56tNVrNiSyaihfWhQu3qsy4k5tahFJK7MWZ/ByKlruLpvS87s0jTW5cQFBbWIxI39B3P49VsLaN6gFg9c0C3W5cQNdX2ISNz4v/eXsX7HfiYM70/dSnTmYUnUohaRuDB+1gbGz9rAiNPb079d41iXE1cU1CISc3PWZ/CHfy3m1I5NuOcnnWNdTtxRUItITG3Zk8WIsXNp3qAW/7j6RKpWsAvTlgV1AolIzGQdyuXmMXPYl53D2Bv6aSjeYSioRSRm/vCvxczfuIsXr+1F5+PqxbqcuKWuDxGJicnz0nkrNZ07zurAT7s3i3U5cU1BLSJRt2HHfn7/7hL6tmnErwZ2inU5cS/soDazqmY2z8zej2RBIlKxuTu/mbgAA54c3FNfHoahNC3qO4FlkSpERCqHt+ekM2NtBved35WkhpVzfunSCiuozSwJuAB4ObLliEhFtn1vNn/5YBl92jTkqj4tY11OuRFui/pp4DdA3uE2MLPhZpZqZqnbtm0rk+JEpGJ5+P2l7MvO4ZHLT6CKujzCVmJQm9mFwFZ3n1Pcdu4+0t1T3D0lMTGxzAoUkYrhixVbeXf+Zm49owMdmmooXmmE06IeAFxsZmnABOAsMxsb0apEpELZue8gv31nIZ2OrcutZ7aPdTnlTolB7e73uXuSu7cBrgI+d/drI16ZiFQI7s4D7y5i5/6DPDU4mRrVqsa6pHJH46hFJKImz9vEB4u+49fndOb45vVjXU65VKpTyN39S+DLiFQiIhXOxoz9/PFfS+jTpiHDT2sX63LKLbWoRSQi9mXncNPrqZjBk1cm68SWo6BJmUSkzLk7905axMotmYz6RV9aNtKJLUdDLWoRKXNvz0nn3ws2c9e5nTm9k4brHi0FtYiUqbXb9vLge0s4qV1jbjldQ/HKgoJaRMrMwZw87pwwn+rVqvDk4J46+7CMqI9aRMrM05+tZNGm3bx4bW+a1a8V63IqDLWoRaRMLN28h5emruWK3kn8tPtxsS6nQlFQi8hRy81z7p+8iAa1Enjggq6xLqfCUVCLyFF7Y+Z65m/cxe8v7KYL1EaAglpEjsp3u7N47KMVnNqxCZckN491ORWSglpEjlhennPvpIUcys3j4Uu7Y6ZRHpGgoBaRIzbqv2l8uWIbD1zQldaN68S6nApLQS0iR2Txpt08+uFyBnY9luv6t451ORWaglpESm3/wRzunDCPhnUSeGxQD3V5RJhOeBGRUnF3fv/uEtZu38fYG/rRqI5GeUSaWtQiUirjZm1g4tx0fnlWRwZ0aBLrcioFBbWIhG3Bxl386b2lnN4pkTvP7hjrcioNBbWIhCVj30FufWMuifVq8PTgZE24FEXqoxaRErk7d7+9gG2Z2Uy85WQaql86qtSiFpESvfL1Oj5fvpUHLujKCUm6QG20KahFpFgLNu7i0Y+W85Pjj2XISRovHQsKahE5rD1Zh7hj/Dya1qvJYz/rqfHSMVJiUJtZTTObZWYLzGyJmf0pGoWJSGy5O/dPWsSmXQf4+9XJ1K+dEOuSKq1wvkzMBs5y971mlgB8bWYfuvuMCNcmIjE0YfZG3l/4Lff8pDO9WzeKdTmVWolB7e4O7A1+TQhuHsmiRCS2VnyXyYPvLeHUjk10gdo4EFYftZlVNbP5wFbgU3efGdmyRCRW9mbncPu4udSrmcCTV2q8dDwIK6jdPdfdk4EkoK+ZdS+8jZkNN7NUM0vdtm1bWdcpIlFwMCePW8bOYe32fTxzVTKJ9WrEuiShlKM+3H0X8CXw0yLWjXT3FHdPSUxMLKPyRCRasnNyuXPCPKat2s5fLz9B83jEkXBGfSSaWYPgfi1gILA80oWJSPQcOJjLTa/P4cPF3/GHC7txRUrLWJckBYQz6qMZ8JqZVSUU7G+5+/uRLUtEoiUz6xA3jE5l9voMHv3ZCQzu0yrWJUkh4Yz6WAicGIVaRCTKDubkMfz1OczdsJN/XH0iF/bQxWnjkSZlEqmk3J17Jy5k+todPD04WSEdx3QKuUgl9dSnK5k0bxN3ndOJS09sEetypBgKapFK6K3ZG/n756sZnNKS28/qEOtypAQKapFKZtqqbdw/eRGndmzCw5d110RL5YCCWqQSWfbtHm4ZO5cOTevy/DW9SKiqCCgP9K8kUkls3nWAYaNnU6dGVUb9og/1amo2vPJCQS1SCWzdk8U1L89kb1YOrw7tQ7P6tWJdkpSChueJVHBp2/cxdNQstmZmM+aGfhzfXJfSKm8U1CIV2LwNO7nhtVTcnTE39KN364axLkmOgIJapIL6ZMl3/HJC6DJao3/Rh3aJdWNdkhwhBbVIBfTW7I3cO2khJ7SozytD+9CkrqYrLc8U1CIViLvz4ldrefSj5ZzasQkvXtubOjX037y807+gSAWRl+c88uEy/jltHRf3bM4TV/SkejUN7KoIFNQiFcDBnDzunbiQSfM2MfTkNvzhwm66hFYFoqAWKed27M3mljfmMmtdBnef24nbzuyg08IrGAW1SDk2deU27n57AbsOHOKZq5K5JFmz4FVECmqRcmjzrgP8+YNl/Gfht3RoWpfXhvWla7NjYl2WRIiCWqQcyc1zXp62lqc/W0WeO/9vYCduPr0dNROqxro0iSAFtUg5sX7HPu56awGp63cysOux/PGibrRsVDvWZUkUKKhF4tyh3DzGzljP4x+voGoV46nBPbk0uYW+MKxEFNQiccrd+WTpFv764XLWbd/HaZ0S+evlJ9C8gWa+q2wU1CJxxt35evV2nv5sFXPW76RD07q8OjSFMzs3VSu6kioxqM2sJfA6cByQB4x092ciXZhIZXPgYC7vzt/E6G/SWLElk2b1a/Lny7ozOKUl1XQllkotnBZ1DnCXu881s3rAHDP71N2XRrg2kUoh61Aub8zcwPNfrGbHvoN0bXYMjw3qwSXJzalRTaM5JIygdvdvgW+D+5lmtgxoASioRY5Cdk4uk+du4h+fr2bTrgOc3L4xd57dkb5tG6mLQ36gVH3UZtYGOBGYGYliRCqDPVmHGDdzA69+vY6tmdn0TKrPY4N6MKBDk1iXJnEq7KA2s7rAROBX7r6niPXDgeEArVq1KrMCRSqCvDxn5roM3k7dyAeLvyXrUB6ndGjC367sySkdmqgFLcUKK6jNLIFQSL/h7pOK2sbdRwIjAVJSUrzMKhQpx3bvP8TYmet5c/ZGNmTsp16NalzeK4mf921F9xa6dqGEJ5xRHwa8Aixz9ycjX5JI+bd9bzavfL2OMdPXszc7h5PbN+bX53TiJ8cfR63q+oJQSiecFvUA4DpgkZnND5bd7+4fRK4skfJp064DvDxtLeNnbSA7J48LezTntjPb0+U4TZgkRy6cUR9fA+pAEynG/I27eOXrdXyw6FsALk1uwa1ntqe9LigrZUBnJoocodw859Ol3/HytHWkrt9JvRrVGDagDdef3IakhposScqOglqklA7l5jFpbjrPfbGGDRn7SWpYi99f2I0rU5KoVzMh1uVJBaSgFgnTody80AkqX6xiY8YBeiTV577zenFOt2N1irdElIJapAR7s3OYPDedf05bx4aM/fRIqs+fLj5ekyRJ1CioRQ5j5ZZMxkxfz6S56ew7mEvPpPr88aIUzuqigJboUlCLFHAoN4+Pl3zHmOnrmbkug+rVqnBhj2YMOakNPZPqK6AlJhTUIoS6N0Z9vY7XZ6xnW2Y2LRvV4r7zunBFSksa1ake6/KkklNQS6WWdSiXMdPX88JXa8jYd5AzOidy/UltOK1TIlWrqPUs8UFBLZVSdk4ub6em84/PV7FlTzandmzC3ed2pmfLBrEuTeRHFNRSqWRmHeKNmRt45et1bMvMJqV1Q5656kT6t2sc69JEDktBLZXC1swsRn2TxtgZ68nMyuGUDk146spkBnRorC8IJe4pqKVC27zrAM99sZq356STk5vHed2bMeL09pyQpClGpfxQUEuFtDUzi+e/WMO4mRtwnEG9kxh+WnvaNqkT69JESk1BLRVKxr6DvPTVGl6bnsahXOeK3kncflYHTZIk5ZqCWiqETbsOMGb6esZMT2P/oVwuTW7BnWd3pI1a0FIBKKilXFuwcRcvfrWGj5d8B8B53Zvxq4Ed6XhsvRhXJlJ2FNRSLm3adYDHPlrOv+ZvpkHtBG46rR3X9W+tLg6pkBTUUq7sy87hxa/WMHLqWgBuP7MDI85oT90aeitLxaV3t5QLeXnOxLnpPP7xCrZmZnNxz+b89rwutGhQK9aliUScglriWl6e89GS73j6s5Ws3LKX5JYNePG63vRq1TDWpYlEjYJa4pK78+nSLTz12SqWfbuH9ol1+MfVJ3Jhj2Y6k1AqHQW1xJ3VW/fy24kLmbN+J22b1OHpwclc1LO5ZrOTSktBLXEjN8955eu1PPHJSmpXr8qjPzuBn/VK0vUIpdIrMajN7FXgQmCru3ePfElSGS1K383v/rWYBRt3cU63Y/nzZd1pWq9mrMsSiQvhtKhHA88Cr0e2FKmMdu8/xN8+XcGYGetpXKc6Tw9O5pLk5uqHFimgxKB296lm1ibypUhlknUol9enp/HcF2vIzDrEkP6t+fW5nalfKyHWpYnEnTLrozaz4cBwgFatWpXVbqWCyc1z3p23iSc/XcmmXQc4rVMiv/1pZ45vrmlHRQ6nzILa3UcCIwFSUlK8rPYrFce67fu4Y/xcFm/awwkt6vPYoB4M6NAk1mWJxD2N+pCoeH/hZu6duIhqVY1nrkrmoh7NqaLhdiJhUVBLRO3NzuEvHyxj3MwN9GrVgGd/3ovmOu1bpFTCGZ43HjgDaGJm6cAf3f2VSBcm5d/Uldu4b9IiNu8+wPDT2nHPTzqToDHRIqUWzqiPq6NRiFQcW/dk8dePljNp7ibaJdbhnREn0bt1o1iXJVJuqetDykx2Ti6vfp3Gs5+v4lCuc+sZ7fnl2R2pmVA11qWJlGsKajlq+RMo/fmDZazfsZ+BXY/ldxd01WWwRMqIglqOyppte3nwvSVMW7WdDk3r8vqwvpzWKTHWZYlUKApqOSK5ec6rX6/j8U9WULNaFf54UTeu7d9aXxaKRICCWkpt7ba93MiwPPEAAA2MSURBVPNOaBpSTaAkEnkKagnbwZw8XvtvGn/7dAXVq1bhqcE9uTS5hSZQEokwBbWUyN35bNlW/vLBMtZt38fZXZryl8tP4Nhj1IoWiQYFtRRryebd/Pk/y/jvmh20T6zDqKF9OKNzolrRIlGkoJYiHczJ46nPVvLSV2uoXyuBhy45nqv7ttKXhSIxoKCWH1m1JZNfvTmfJZv3MDilJfef35X6tTVPtEisKKjle3l5zuvT03jkw+XUrVGNkdf15tzjj4t1WSKVnoJaANiyJ4u7317AtFXbOatLUx79WQ8S69WIdVkigoK60nN3Js7dxP+9v5TsnFwevrQ71/RrpS8LReKIgroS25ixn/snL2Laqu2ktG7Io4N60D6xbqzLEpFCFNSV0L7sHJ77YjWvfL2OalWMhy45nmv7tdYVV0TilIK6kpm3YSf/7835pO3Yz6XJzfnNT7voiisicU5BXUkczMnjhS/X8PfPV3HcMTV5c3h/+rVrHOuyRCQMCupK4MsVW3no/aWs3baPy05swZ8uOZ5jampctEh5oaCuwJZ9u4fHP17B58u30rZJHV65PoWzux4b67JEpJQU1BXQiu8yeWbKSj5Y9B31alTj/vO7MPTktlSvptO/RcojBXUFcSg3j0+WbGHsjPVMX7uDujWq8cuzOnDDKe10+rdIOaegLuc27zrA+FkbmDB7I9sys2nRoBb3/KQz1/RrRYPa1WNdnoiUgbCC2sx+CjwDVAVedve/RrQqOayc3DwWb97D9DU7+O+a7XyzejsOnNm5Kdf2b8XpnZpSVeOhRSqUEoPazKoCzwHnAOnAbDN7z92XRrq4yszd2XMgh/Rd+1m9dS8rvstk6bd7SE3byd7sHAA6Nq3LiNPbc3XfVrRsVDvGFYtIpITTou4LrHb3tQBmNgG4BDhsUG/fm82rX6/7wTI/iiILcy/LvZUNd3Cc3DzIc8e90H138jw0Q11ecD83z8nOyWNfdk7odjCHfdm57MvOYcueLPYdzP1+/9WqGO0S63BxcnNOateY/u0aa9IkkUoinKBuAWws8Hs60K/wRmY2HBgOUP24Djz0vhrc+cygihlVzf53v0rofo1qValboyq1q1ejbo1qNKlbnVaNa3N650RaNKhF8wa16NC0Lm0a19GoDZFKKpygLqrD80dNWncfCYwEOLFXb//yD+eGt6cjVJaTu5XVrqqYUSUI46pVLPgdzUQnIkclnKBOB1oW+D0J2FzcA6pWMQ0JExEpI+F8lp4NdDSztmZWHbgKeC+yZYmISL4SW9TunmNmtwMfExqe96q7L4l4ZSIiAoQ5jtrdPwA+iHAtIiJSBA0jEBGJcwpqEZE4p6AWEYlzCmoRkThnkTgd28wygRWHWV0f2H2YdU2A7WFsF+7+IrVdfp3xXGPh7Qq+tpE+blnss3C98fAaFqe41zfSx47390NZvW/KSz4U3K6k90XB/XV293pFbuXBXBRleQNSi1k3MpzHFbdduPuL1Hb5dcZzjYW3O9J/kyPZriz2WbjeeHgNw3lPxHmNMXk/lNX7przkQ8HtSnpfhPtvEouuj39XkO1ieex43y6Wx47lcw5XvD8XvW/ibLtIdX2kuntKtB4XbeWlzoLKW82qN7LKW71Q8WsubttItahHRvlx0VZe6iyovNWseiOrvNULFb/mw24bkRa1iIiUHQ3PExGJcwpqEZE4F9WgNrNcM5tf4NammG2/NLOYfHFgZm5mYwr8Xs3MtpnZ+7GoJ1xmdllQe5dY13I45fW1zWdme2NdQ2mVVHMs/68VqiPu37+FmdkDZrbEzBYGmfajq1+VhWi3qA+4e3KBW1qUjx+ufUB3M6sV/H4OsKk0OzCzsGYmLGNXA18TmjM8bMEFjKPlqF9bqbCO6P0bK2Z2EnAh0MvdewAD+eFlC8tMzLs+zKy3mX1lZnPM7GMza1Zg9bVm9l8zW2xmfaNc2ofABcH9q4Hx+SvMrG9Q17zgZ+dg+VAze9vM/g18Es1izawuMAC4geCNbmZnmNlUM5tsZkvN7EUzqxKs22tmD5nZTOCkaNbKkb2208wsucB235hZj6hW/b9jn1HwE4CZPWtmQ4P7aWb2JzOba2aL4qV1WFzN8aCY9+/hXufzzWy5mX1tZn+P0SeyZsB2d88GcPft7r75cJkWfHJ5+kgyLdpBXatAt8dkM0sA/gEMcvfewKvAnwtsX8fdTwZuDdZF0wTgKjOrCfQAZhZYtxw4zd1PBP4A/KXAupOA6939rKhVGnIp8JG7rwQyzKxXsLwvcBdwAtAeuDxYXgdY7O793P3rKNd6JK/ty8BQADPrBNRw94VRq7h0trt7L+AF4O5YF1NOHO79+yPB++Yl4Dx3PwVIjFKNhX0CtDSzlWb2vJmdHqlMi/bH8wPuXrBV1B3oDnxqoQvAVgW+LbD9eAB3n2pmx5hZA3ffFY1C3X1h0Id+NT++aEJ94DUz60joQr8FLxD5qbtnRKPGQq4Gng7uTwh+/w8wy93XApjZeOAU4B0gF5gYgzqP9LV9G/i9md0DDANGR6XYIzMp+DmH//1hlOId7v1blC7AWndfF/w+Hhge2fJ+zN33mllv4FTgTOBN4GEikGmx6EctyIAl7n64j96FB3lHe9D3e8ATwBlA4wLL/w/4wt0vCwLnywLr9kWptu+ZWWPgLEJ9v07ozeGEQvBwr2GWu+dGr8ofKdVr6+77zexT4BLgSiCWX37l8MNPozULrc8OfuYS+/9j+UqqOWaKef++R9E1W3QrPLzg/9CXwJdmtgi4jQhkWqz7qFcAiUGnPGaWYGbHF1g/OFh+CrDb3cOdgausvAo85O6LCi2vz/++ABsa1YqKNgh43d1bu3sbd28JrCPUeu5roQsTVyH0eka7m+NwjuS1fRn4OzA7Rp9a8q0HuplZDTOrD5wdw1rCFc81H+79C0XXvBxoZ/8bNTY4qtUGzKxz8MkvXzKwjAhkWkz/2rv7QTMbBPw9+IeoRujjT/7Fc3ea2X+BYwh93I12fenAM0WseozQx/NfA59Ht6oiXQ38tdCyicAtwPRg3QnAVGBydEsr2pG8tu4+x8z2AKOiUOKPWGgkT7a7bzSzt4CFwCpgXizqCUc5qflw79+fAz+q2d0PmNmtwEdmth2YFcVaC6oL/MPMGhD6xLKaUBfMSMo403QKeQVmZmcAd7v7hbGupSyYWXNCHzO7uHteDI7fE/inu0d7BNIRK481h8PM6gZ9xAY8B6xy96diXVdxzOxLQv8fU0v72Fh3fYiExcyGEBod8kCMQnoEoS+CfhftYx+p8lhzKdxkZvMJtVTrExoFUmGpRS0iEuci2qI2s5Zm9oWZLbPQaZZ3BssbmdmnZrYq+NkwWH5OMEh8UfDzrAL76h0sXx0McI+bb35FRCIp0l0fOcBd7t4V6A/cZmbdgHuBKe7eEZgS/A6ha4td5O4nANcDYwrs6wVCHfUdg9tPI1y7iEhciGhQu/u37j43uJ9JaOhKC0JjYV8LNnuN0FlJuPs8d98cLF8C1AyG5jQDjnH36R7qq3k9/zEiIhVd1L5MDMY8nkjoC6Fj3f1bCIU50LSIh/wMmBecR98CSC+wLj1YJiJS4UVlHLWFJlyZCPzK3feU1L0cDBB/FDg3f1ERm+lbUBGpFCLeog4mKZkIvOHu+XMgbCkwo1QzYGuB7ZMInZQxxN3XBIvTgaQCu00CNiMiUglEetSHAa8Ay9z9yQKr3iP0ZSHBz38F2zcgNBHLfe7+Tf7GQfdIppn1D/Y5JP8xIiIVXUTHUQfns08DFgH5JyncT6if+i2gFbABuMLdM8zsd8B9hE4XzXeuu2+10BUoRgO1CM1nfIdrELiIVAI64UVEJM7pFHIRkTinoBYRiXMKahGROKegFhGJcwpqEZE4p6CWcsHMHjSzw17R28wuDSb8Kmk/P9jOzB4ys4FlVadIJCiopaK4FCgxqAtv5+5/cPfPIlaVSBlQUEvcMrMHzGyFmX0GdA6W3WRms81sgZlNNLPaZnYycDHwuJnNN7P2we2jYF7zaWbW5TDbjQ6u24mZpZnZX8xsupmlmlkvM/vYzNYEV0vJr+ueoIaFZvanGLw0UsnEy6XsRX7AzHoDVxGacbEaMBeYA0xy938G2zwM3ODu/zCz94D33f2dYN0UYIS7rzKzfsDz7n5WEdsVPvRGdz/JzJ4idCbsAKAmoWl3XzSzcwnNh96X0GRh75nZae4+NWIvhlR6CmqJV6cCk919P0AQsADdg4BuQOgq0B8XfmAwW+PJwNsFgrhGmMfNP84ioG4wj3qmmWUFc9GcG9zyr+Jdl1BwK6glYhTUEs+Kmt9gNHCpuy8ws6HAGUVsUwXY5e7JR3DM7OBnXoH7+b9XI9SKfsTdK/TFVCW+qI9a4tVU4DIzq2Vm9YCLguX1gG+D6XOvKbB9ZrAOd98DrDOzKyA0i6OZ9Sy83RH6GBgWtNoxsxZmVtSFL0TKjIJa4lJwCbc3gfmE5jOfFqz6PaHZFz8Flhd4yATgHjObZ2btCYX4DWa2gFD/8iWH2a60dX0CjAOmm9ki4B2OLvhFSqTZ80RE4pxa1CIicU5BLSIS5xTUIiJxTkEtIhLnFNQiInFOQS0iEucU1CIicU5BLSIS5/4/1bdrDIZ/3GAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "positive_MI_df.plot(x=\"datetime\", y=\"positive\", title=\"MI positive cases\")\n",
    "positive_us_df.plot(x=\"datetime\", y=\"positive\", title=\"US positive cases\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - What do you notice about the number of positive cases in the US and in Michigan? How could we characterize the rate at which infections are changing? In particular, what are ways that might be able to compare the two rates on similar scales?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2.3 Time Calculations\n",
    "\n",
    "One of the most common ways to compare these rates in epidemiology is [\"Doubling Time\"](https://en.wikipedia.org/wiki/Doubling_time). In the space below we are going to look at the Doubling Time of the virus. But to do that, we will have to make some calculations regarding time differences.\n",
    "\n",
    "You can subtract two `DateTime` elements which yields a `Timedelta`, a data type representing a time difference. Let's try that below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<class 'datetime.timedelta'>\n",
      "30 days, 0:00:00\n",
      "['__abs__', '__add__', '__bool__', '__class__', '__delattr__', '__dir__', '__divmod__', '__doc__', '__eq__', '__floordiv__', '__format__', '__ge__', '__getattribute__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__mod__', '__mul__', '__ne__', '__neg__', '__new__', '__pos__', '__radd__', '__rdivmod__', '__reduce__', '__reduce_ex__', '__repr__', '__rfloordiv__', '__rmod__', '__rmul__', '__rsub__', '__rtruediv__', '__setattr__', '__sizeof__', '__str__', '__sub__', '__subclasshook__', '__truediv__', 'days', 'max', 'microseconds', 'min', 'resolution', 'seconds', 'total_seconds']\n",
      "30\n"
     ]
    }
   ],
   "source": [
    "t1 = datetime.datetime(2020, 7, 25)\n",
    "t2 = datetime.datetime(2020, 6, 25)\n",
    "diff = t1 - t2\n",
    "print(type(diff))\n",
    "print(diff)\n",
    "print(dir(diff))\n",
    "print(diff.days)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By looking at the `dir` of the `Timedelta` instance named `diff` we see the various data members, one of which is `.days`. Calling `.days` yields an integer value we can do work with. You could have done `.seconds` or changed the `.resolution` but for the work we are doing, the values are recorded as days so that is what we will work with"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "------\n",
    "## 3. Creating the \"Doubling time\" plot\n",
    "\n",
    "Now we calculate an estimate of the [doubling time](https://en.wikipedia.org/wiki/Doubling_time) of positive COVID-19 tests (which gives us a sense of how quickly the virus is spreading).  The doubling time is the amount of time it takes the number of cases to increase by a factor of two.  We use the instantaneous approximation of the doubling time - for times $t_1$ and $t_2$, and number of cases at those times $n_1$ and $n_2$, we can approximately the doubling time $\\tau_d$ as:\n",
    "\n",
    "$\\tau_d \\simeq (t_2 - t_1) \\frac{ln(2)}{ln(n_2/n_1)}$\n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Given that equation, let's write a function called `d_time` which takes as arguments two rows from our time-series, assumed to have the column headers \"datetime\" and \"positive\". Return from that function the doubling time, in days, between those two time periods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here for d_time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# test your code, show that it works"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.2 Calculating doubling time \n",
    "\n",
    "We have the option of specifying the time delta (every day, every 5 days, every week) to see how the doubling time might change. \n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Let's write a function `d_plot` that takes two arguments.\n",
    "- first is a DataFrame, either the US or the MI that has the two columns, \"positive\" and \"datetime\"\n",
    "- time interval. Conveniently, each element in the DataFrame represents one day, so an interval of 5 means we are using DataFrame rows that are 5 indicies apart. \n",
    " \n",
    "`d_plot` should return an array of doubling times, from the first reported cases up to the latest cases.\n",
    "\n",
    "Here's a hint. The `zip` function will create a list of pairs. If you give it two lists, one a range of indicies from earliest to latest by the interval, and the other a range of values, it will return a list of pairs where the pairs are of (time, value). You do not _need_ to use `zip` to do this calculation, it's just relatively elegant."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.3 Plot doubling times\n",
    "\n",
    "A list of doubling times might not be the most useful way to represent our data. So let's instead make a plot to represent the data.\n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This -  Plot a side-by-side plot of:\n",
    "- doubling time\n",
    "- positives\n",
    "\n",
    "for a given data set (US or MI)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3.4 Comparing Intervals\n",
    "\n",
    "In the previous section, you picked an interval to use for your doubling time calculation. What happens if we shrink or lengthen that interval?\n",
    "\n",
    "<font size=8 color=\"#009600\">&#9998;</font> Do This - Make 3 doubling time plots for different choices of interval for the same data (US or MI). What happens to the representation of the data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font size=8 color=\"#009600\">&#9998;</font> Do This -  What do you notice about this plot? Do you see a problem with how we are representing the data? How might we fix this plot to make this comparison more meaningful?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------\n",
    "### Congratulations, you're done with your in-class assignment!\n",
    "\n",
    "Now, you just need to submit this assignment by uploading it to the course <a href=\"https://d2l.msu.edu/\">Desire2Learn</a> web page for today's in-class assignment submission folder. (Don't forget to add all of the appropriate names in the first cell)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
