{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Earthquakes and oil mining in Oklahoma\n",
    "> Of course, earthquakes have a big impact on society, and recently are connected to human activity. In this final chapter, you'll investigate the effect that increased injection of saline wastewater due to oil mining in Oklahoma has had on the seismicity of the region. This is the Summary of lecture \"Case Studies in Statistical Thinking\", via datacamp.\n",
    "\n",
    "- toc: true \n",
    "- badges: true\n",
    "- comments: true\n",
    "- author: Chanseok Kang\n",
    "- categories: [Python, Datacamp, Statistics]\n",
    "- image: images/compare_mag_okh.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import dc_stat_think as dcst\n",
    "\n",
    "plt.rcParams['figure.figsize'] = (10, 5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variations in earthquake frequency and seismicity\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDA: Plotting earthquakes over time\n",
    "Make a plot where the y-axis is the magnitude and the x-axis is the time of all earthquakes in Oklahoma between 1980 and the first half of 2017. Each dot in the plot represents a single earthquake. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>depth</th>\n",
       "      <th>mag</th>\n",
       "      <th>magType</th>\n",
       "      <th>nst</th>\n",
       "      <th>gap</th>\n",
       "      <th>dmin</th>\n",
       "      <th>rms</th>\n",
       "      <th>...</th>\n",
       "      <th>depthError</th>\n",
       "      <th>magError</th>\n",
       "      <th>magNst</th>\n",
       "      <th>status</th>\n",
       "      <th>locationSource</th>\n",
       "      <th>magSource</th>\n",
       "      <th>loc_name</th>\n",
       "      <th>loc_admin1</th>\n",
       "      <th>loc_admin2</th>\n",
       "      <th>loc_cc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1974-12-16 02:30:21.400</td>\n",
       "      <td>35.330</td>\n",
       "      <td>-97.480</td>\n",
       "      <td>10.0</td>\n",
       "      <td>2.6</td>\n",
       "      <td>ml</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>m</td>\n",
       "      <td>tul</td>\n",
       "      <td>Moore</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Cleveland County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1975-09-13 01:25:02.800</td>\n",
       "      <td>34.139</td>\n",
       "      <td>-97.369</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3.4</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Wilson</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Carter County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1975-10-12 02:58:11.200</td>\n",
       "      <td>34.816</td>\n",
       "      <td>-97.406</td>\n",
       "      <td>20.0</td>\n",
       "      <td>3.2</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Maysville</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Garvin County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1975-11-29 14:29:40.900</td>\n",
       "      <td>34.521</td>\n",
       "      <td>-97.347</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3.5</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>slm</td>\n",
       "      <td>Wynnewood</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Garvin County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1976-04-16 18:59:44.200</td>\n",
       "      <td>36.107</td>\n",
       "      <td>-99.875</td>\n",
       "      <td>5.0</td>\n",
       "      <td>3.4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Arnett</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Ellis County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                      time  latitude  longitude  depth  mag magType  nst  gap  \\\n",
       "0  1974-12-16 02:30:21.400    35.330    -97.480   10.0  2.6      ml  NaN  NaN   \n",
       "1  1975-09-13 01:25:02.800    34.139    -97.369    5.0  3.4      lg  NaN  NaN   \n",
       "2  1975-10-12 02:58:11.200    34.816    -97.406   20.0  3.2      lg  NaN  NaN   \n",
       "3  1975-11-29 14:29:40.900    34.521    -97.347    5.0  3.5      lg  NaN  NaN   \n",
       "4  1976-04-16 18:59:44.200    36.107    -99.875    5.0  3.4     NaN  NaN  NaN   \n",
       "\n",
       "   dmin  rms  ... depthError magError magNst    status locationSource  \\\n",
       "0   NaN  NaN  ...        NaN      NaN    NaN  reviewed              m   \n",
       "1   NaN  NaN  ...        NaN      NaN    NaN  reviewed             us   \n",
       "2   NaN  NaN  ...        NaN      NaN    NaN  reviewed             us   \n",
       "3   NaN  NaN  ...        NaN      NaN    NaN  reviewed             us   \n",
       "4   NaN  NaN  ...        NaN      NaN    NaN  reviewed             us   \n",
       "\n",
       "   magSource   loc_name  loc_admin1        loc_admin2 loc_cc  \n",
       "0        tul      Moore    Oklahoma  Cleveland County     US  \n",
       "1        tul     Wilson    Oklahoma     Carter County     US  \n",
       "2        tul  Maysville    Oklahoma     Garvin County     US  \n",
       "3        slm  Wynnewood    Oklahoma     Garvin County     US  \n",
       "4        tul     Arnett    Oklahoma      Ellis County     US  \n",
       "\n",
       "[5 rows x 26 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_csv('./dataset/oklahoma_earthquakes_1950-2017.csv', skiprows=2)\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from datetime import datetime as dt\n",
    "import time\n",
    "# Resource: https://stackoverflow.com/questions/6451655/how-to-convert-python-datetime-dates-to-decimal-float-years\n",
    "\n",
    "def toYearFraction(date):\n",
    "    def sinceEpoch(date): # returns seconds since epoch\n",
    "        return time.mktime(date.timetuple())\n",
    "    s = sinceEpoch\n",
    "\n",
    "    year = date.year\n",
    "    startOfThisYear = dt(year=year, month=1, day=1)\n",
    "    startOfNextYear = dt(year=year+1, month=1, day=1)\n",
    "\n",
    "    yearElapsed = s(date) - s(startOfThisYear)\n",
    "    yearDuration = s(startOfNextYear) - s(startOfThisYear)\n",
    "    fraction = yearElapsed/yearDuration\n",
    "\n",
    "    return date.year + fraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "times = pd.to_datetime(df['time']).apply(toYearFraction).to_numpy()\n",
    "mags = df['mag'].to_numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlcAAAE9CAYAAAAmvEclAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzda4xk6X3f9+/znOdc6lRV33t2dmZ3Z5ZcSqZIiVxqSYe0Y8eSYluxJMiRgSCIYEQOohdOFAWOkQRxECUvHMCB7VyAWIhiIBAiB0aiWIGhxLIAJzItkZS1a0okpRW55F5md2Z7pq9VXZdzec7z5MXp6unp6Z6+TFV31fT/Awxm+1LnPOdW8+w5v/r/lfceIYQQQggxHvqyByCEEEII8SyRyZUQQgghxBjJ5EoIIYQQYoxkciWEEEIIMUYyuRJCCCGEGCOZXAkhhBBCjJGZ5MKVUgvA3wU+CXjgL3nvv3zc76+srPjbt29PckhCCCGEEGPxxhtvbHjvVw9/f6KTK+C/B37Ne/8XlFIRkD7pl2/fvs3rr78+4SEJIYQQQjw9pdR7R31/YpMrpdQc8CeAfxvAe18AxaTWJ4QQQggxDSaZufoIsA78L0qpryql/q5SqjnB9QkhhBBCXLpJTq4M8Bng5733rwJ94D89/EtKqZ9WSr2ulHp9fX19gsMRQgghhJi8SU6uPgA+8N7/9t7Xv0w92XqE9/4XvPevee9fW119LBMmhBBCCDFTJja58t6vAe8rpb5771s/CPzBpNYnhBBCCDENJv1pwZ8B/t7eJwXfBn5qwusTQgghhLhUE51cee9/F3htkusQQgghhJgmUqFdCCGEEGKMJv1YUAghhBDiqRXWUVaOMNBEZrrvDcnkSgghhBBTrbCOu9sDPKCAm4vpVE+wpndkQgghhBBAWTk80IwNfu/raSaTKyGEEEJMtTDQKKCfW9Te19NMHgsKIYQQYqpFRnNzMZXMlRBCCCEmb5aC3k8jMrOzfTK5EkIIIWbUrAW9rwo5AkIIIcSMmrWg91UhkyshhBBiRs1a0PuqkMeCQgghxIyataD3VSGTKyGEEGKGzVLQ+6qQoyGEEEIIMUYyuRJCCCGEGCOZXAkhhBBCjJFMroQQQgghxkgmV0IIIYQQYySTKyGEEEKIMZLJlRBCCCHEGMnkSgghhBBijGRyJYQQQggxRjK5EkIIIYQYI5lcCSGEEEKMkUyuhBBCCCHGSCZXQgghhBBjJJMrIYQQQogxksmVEEIIIcQYyeRKCCGEEGKMZHIlhBBCCDFGMrkSQgghhBgjmVwJIYQQQoyRuewBCCGEEEIAFNZRVo4w0ERmdu//yORKCCGEEJeusI672wM8oICbi+nMTrBmc9RCCCGEeKaUlcMDzdjg976eVRO9c6WUehfYBSrAeu9fm+T6hBBCCDGbwkCjgH5uUXtfz6qLeCz4p7z3GxewHiGEEELMqMhobi6mJ2auZiGXJZkrIYQQQkyFyDx5wjQruaxJj8gDv66UekMp9dNH/YJS6qeVUq8rpV5fX1+f8HCEEEIIMatmJZc16cnVH/Pefwb4YeDfU0r9icO/4L3/Be/9a97711ZXVyc8HCGEEELMqlnJZU30saD3/t7e3w+UUr8CfA744iTXKYQQQohn02lzWZdtYqNSSjWVUu3RfwN/GvjGpNYnhBBCiGdfZDTN2EztxAome+fqOeBXlFKj9fxv3vtfm+D6hBBCCCEu3cQmV977t4FPTWr5QgghhBDTaHrvqQkhhBBCzCCpcyWEEELMsF5mGRSWNDK0EvlnfRrIURBCCCFmVC+zfPW9LRwejeLVW0sywZoC8lhQCCGEmFGDwuLwLLcSHJ5BYS97SAKZXAkhhBAzK40MGsVmL0OjSCO5azUN5CgIIYQQM6qVGF69tTTTmavCOvq53S8MOu01rE5j9o6CEEIIIfa1ktmcVEE9sXp3o8cHWwM2ejnL7ZgXF1Nur7RmeoI1uyMXQgghxEwrK0dZeYzRRCYg1Jqy8lPbkPm0ZHIlhBBCiEsRBpowUFjrKGxF6RxhoKa2IfNpzeZ9RCGEEELMvMhobq+0WG0nkrkSQgghhBiHyGgiE132MMZqtqeGQgghhBBTRiZXQgghhBBjJJMrIYQQQogxksyVEEKcoLBuP2w760FbIc5j3NfAs35NyeRKCCGeoLCOu9sDPKCAm4vpM/mPgRDHGfc1cBWuqWdra4QQYszKyuGBZmzwe18LcZWM+xq4CteUTK6EEOIJwkCjgH5uUXtfC3GVjPsauArXlDwWFEKIJ4iM5uZi+kznQ8R0uaw80nHrPe01cJZxL6R1Xasw0JSVo7COQVE3b55vRI/0Sjy83FnIa8nkSgghTlAXOZzON3HxbLmsPNJJ6z3pGjjtuA/+nt17HOiBO+s91vsFgYbVdsLnXl6hlZjHlrvaTljfzaY+rzV9IxJCCCGuqMvKIz3tek/7+oO/Vzdo9kRGM7AVgYb5NMZWMCjskcsdFHYm8loyuRJCCCGmxGXlkZ52vad9/cHfqxs0KwrrSE1A5aAzyDEBpJE5crlpZGYiryWPBYUQQogpcVkZv/Oud5R/8v5hjuqkxssHf6+wjs6w4HtuLlC6elkrzWQ/c3XUuCIz/RlImVwJIYQQU+SyMn5nXe8oD1VUjgedjGvzCVGgacZHTy0O56fCQHNvZ8BaJ8NWFSjFzcWUnWFBM3k4QTs8rlnIQE736IQQQggxlUZ5qMhoHJ4o0KfOW43yU2XlaUQBxgTYihOXMStkciWEEEKIMxvloQrr0CiKyp06bzXKT4WBYlhUWFthAk5cxqyQx4JCCCGEOLODeaibCylK8cQc1FH5qdsrLVbbCbA3+TphGbNCJldCiKk3C0UDhZh1Z7nOepllUFjSyDxW8LP+JGC9jO1eQScrSENDMzGEhzJZdX4qemT9h5dfWkcnK5hPIhZb0Ynj955Ln6TJ5EoIMdWuQpNXIS7bWa6zXmb56ntbODwaxau3lo4s+NmMDP/0mw8onKXTr/j8R1dYaIZHLvvga4eF5X43I9CKwbBko18ShQqN5oe+5/qRE6zjwvWX9X4h71BCiKl2FZq8CnHZznKdDQqLw7PcSnD4Ywt+bvRzHI6VVgPnKyrcscs++Nph4bAVLLcSdgtLVlZcn09x1HewnjT+04brJ00mV0KIqXYVmrwKcdnOcp2lkUGj2OxlaNSxBT9XmjEazUZviFYBAfrYZR98bSPSmAA2exntyJCEAWudARrNfHL0Y8GzhusnTR4LCiGm2mmKG0omS8ySizhfz7qOszRnVgq++/ocg9KShgal6u8ftYwf+p7rdLKCUGtQ1H9TP1rsDOvvh4FmUNYTsshorrVb3FpuHZu5Om7bRsVJTxOunzSZXAkhpt6TigZKJkvMkos4X8+7jtM2Zx7lmhZbEe9vDh7LNx1cxmIropkY3t3osdbJAM9cI+RBN2N7UJLbiqoCHUBvaPnYtTlur6bcXmnRSpJHlnPctgFT9x4w8bUrpQKl1FeVUr866XUJIa4eyWSJWXIR5+uk1nE416RRpy4eOioWmkSG3awiKz0LaQx4BrmlEQQYo0D5vYbOpytEWi97+t4DLmJq97PAmxewHiHEFSSZLDFLLuJ8ndQ6DueaRhOs0xQPHRULzQpLOwlIQsXOIAcUaWwYVhXWevBqr6Hz6QqR1suevveAiT4WVEq9APw54K8Df2WS6xJCXE2X1ehWiPOYxPl6OIM0qWvi4HKvtRJK57g5n9bP4vaM6lONntGFWhOFmoVGhFKAhzQ2rLYSBqWtM1dG161wbD3e+cbRpRYONogebdvhcU3Le8CkM1f/HfAfA+0Jr0cIcYXNQiNXIUbGeb4el6+a1DUxWuaDboYH7N4jOLPXhPl+N6O0Fd95MODWcpNBYbm10uS9jT6NWHNnY7j//Veea9PL6uzWeiejcI7OoGS5HfPiYp27iow+Muu13SvqrNfg6KzXZZvYSJRSPwI88N6/ccLv/bRS6nWl1Ovr6+uTGo4QQgjxzLmMvNHBddb5KP9IfapGHOF8RWgUDkduHQ5HZIJHvl85//DRonN7nxYMCLV+JHd13qzXZZrkNO+PAT+mlHoX+PvADyilfunwL3nvf8F7/5r3/rXV1dUJDkcIIYR4tlxG3ujgOut8lHqkPtUwL9AqoLQejSY2Go2msNUj3w+0ejhR0vVEqbAVpXOP5K7Om/W6TMp7P/mVKPWvAH/Ve/8jT/q91157zb/++usTH48QV4nUgJos2b9iEsbR52+StnsFG/2MNKzX18lK0jCon036esLTyYr9HNWgsCg0jUhTOk8aBqy0E8JA73/irzMs6A5LGmHAfBox34iITP3zonRs9DMGecVKK2Y+jS69lhWAUuoN7/1rh78vda6EeIZJDajJkv0rJuEs51VhHeu7df6pn1siM/lzsJdZfvfOFg96Of3C0umXtJsha1tDPv3iEs1Y00pC3nqwy1pnyEavINCauUaA8oprcw2SUPGJGwvcWmmyMyjYGRZ8+a1Nolg/spxrcw2MUby91uWt9T5aeZZaCT/8yRtPbOJ82S7kXcB7/xsn3bUSQozfNNZ/eZbI/hWTcJbz6jLOwUFhKZxjIY2pnCMrHXORwXlHFNX5qZ1hgVGKyBi880QajAoYFBVxqIhCzbCs77h5wDlwvnpsOcPSEhnNbmHRyrMy18BWHNtjcFrI/2IJ8QybxvovzxLZv2ISznJeXcY5mEaGSGt2BjmB1iShpltYtNIURZ2fWmhEWO8prEVpReHA+oo0CshLT1E6GqEhjQwK0Bq0Ch5bTiM0FNbRjgzOKza6Q0zAsT0Gp8WFZK5OSzJXQoyfZIImS/avmISznFeXcQ4e7A0IPFKzKo0MkdHsDIr9mlelc4S6/tnovxea0X6phbJy9DN75HJG9a06w3p5K81kah4JSuZKiCtq2uq/PGtmcf/KhPDynHbfj+O8GlfQ/WABT6XA+/rRYFk9nCRdMwmFdfTzep0bPcuHnQFF4Vida/D8fIPS1csZBdX7uWV99+FySufq8YaGNDaP1etqJWZ/m7Z7BVFYV2fvZ3a/sXMzMVNxbsvkSgghrhAJ4V+eSez745bZyyxffW9rv2zBq7eWzjXBOqqA54dbAz7czegOS5aaMZ96YZFbK03u7Qz4YGvAu1t9/uBul/e2dikrz825Bp94YZF2bIhCxWIaca3d4EF3yB+u7WKdp7SObmZR2mO04rO3VvjY9dZ+IVFgf5uKqmKrV/LxG/NYV/EH97oEus5tffz5eebS8NLPbbmihBDiCpEQ/uWZxL4/bpmDwuLwLLcSHL5+PPcUyz9YwHNgK8DTTiKMfhhMLyuPMZqq8pTWEgUBaRSCgs1+jkcxn8ZkpWc3L0ErlFL1oz/rqFxFKzZoFVBRPdbAebRNc40Yh0Mp2M0qbAXX51NKV7GbF1NxbsvkSgghrhAJ4V+eSez745aZRgaNYrOXoVGk0fkeVB1VwDM1AaDYzQqsexhMDwOFtY4gUITGUFQVg6IED8vNGIWnM8hJQkU7DsF5vPcU1hEaTaADernF+YqA4LEGzqNt6g5zNBrvoZ0EmADWOgNCHdCOo6k4tyXQLsQESKZFTDM5PyfrSft3lBkaNTQexzHY7hX7maODQe+zZq5GmSmo74SNxjXanp1+QScrmU/qu1EbvWw/pD4qHjr6+52NLm9vDJhPDLdX2oSBYlBU5IUlLyuCIKC5N6ZOLyOvPN45dooKBbywkLLSjkljsx9gf9DJ+M7GLsorllsRaWyYb0QMcstGP2c+CfebQKeh2Q/MT5IE2oW4IJJpEdNuFkP4s+JJ1/+o4Ocov3RtPiEK9FO9R/Qyyzfu7uw9shs8kq1qJacPshfW8e5Gj7VOBniuzzf2806RqUPjX/r2Bg6Hc/DKtTa7WckHOwM2dws8YLTildU2vTzjy+/soPHEoaF0ngfdgntbA97v9NkdWrxSLDZCXlxo8p3NXQINDzpDgkCTW5hLAhabCZ94fp4Xlhp89vYSv/b1NTZ6Q3pDx8dvzvHytSbX5+oq74FWfGute2zz54smV5cQYyaZFiGuridd/4fzS+NoPDzObFVZeRpRQBKZx/JOnazA4fazTTvDHGM0ePDeExlNoBTWOzb7FRrPjcUmeVlxf6eo61kZsKUnNoaGCSgqT99aFIqGMVRegQtIAg0qoKwczURjK3hva0BZWVbmGnhdkduKUGuGRT3uyOgnNn++aDK5EmLMJNMixNX1pOv/cH5pHI2Hx5mtCgPFsKjICvtY3mk+idDo/WzTQiPGWgcKlFIU1lF5j1Ga5WaAQ3Fvu08cBjy3ENWV2C2YUJFby9BWRIGiaQwez9BaAuVBV2SVA1/tlVlwmABuLaWEgWGjO0S5gNgElM7RiOpxF9Y9sfnzRZPHgkKMWWTq2/xXJdMi+R0hHjrp+m/GhqhyXLu5MJbMVSsxfGS1xVo3Yy4x+3euzlp2ITKa1VZSN0PWmvn00SKdzcTwhY+t0BnW2a4w0Gz0M+YahuFyhVaKxWZUfzJwN+C1Fx1r3ZIk0jSCgI8uNbi91MC7Re71Coqi4vpCzAsLLYbFMpuDkn5u2emXDIqKhTSkEQdEgWapYVjrDLkxZ1hqtXlhIaWZGBT1mEdjeWEpJdy7mxYavV9P6zLI5EqICbgqmRbJlwnxuKOu/ydlmp7Gdq/gS9/eYGBL3t+oGx7Pp+bMda1G2a2DNaSaccDNxRRgv9bV5m5Bbh1/eK/LwFa8+6AHSrPYNDTDkM1ezh/e3+aDzR7b/QqnoBkGzKcRNxabbPUzKge5tczFIS8tt1lqRvXrHuyQDUq6paNpNHnlmU9DKlsRhAFUMJeGvLDYIo0DVudj5u6GrKQRu0U9MXt5uUUS1xPcwroLaWR9FHkXFEKcm+TLhDidkzJN5zXKQi004v2Gx+fJXh1VQ2p0TR/OirkKSlcRKQUKIg1xYNjpF2S2wmiDI0QFdf4Jrcgqh9FQOQUoYhMSBAGZLckqS+UduIDAhKAUKtB4VfcbrJSmqiA0hiAIGJQlpXMspTHD3NHJC1pxhDF15stWjCXP9jRkciWEODfJlwlxOidlms5rlIXaGeb7DY/Pk706qobU6Jo+nBXTAYQ6oPAePBQO8sqy0IxITIB1Fk2Jr+r8E86TBBrrINAe8OS2pKoqEhOSBIZAadAVlS3Be3zlUJ665pV3BAGU1lJVFWkYEmrN1iCnEWvm44heXmBtnfkyAWPJsz0NeSx4jKNyJFc1W3JVt1uc7Krly85CrpuHLnpfTMu+PzyOGwspzdg80rT4uNeM+vidtA2LrYgvvLLCWjfjCx9ZZT6NzlTXqqwcRekonePFpZRBWfHR1YDQaEJdN1UeFBaloBUZ4uWUzrDg+lxEP3fc+lgDWznQmpV2xLC03FyMud/J6QxKBpllaT5ioRGy1StZbhh2ByV9C8tpwsurDbZ6BXle8FIrZtgwhAG8fG2R1bmEQCsCrVB4cutpRgHPtRIAjAm4vdIkDDR3tvuoFcX1+YTSOYZFxaJkrqbLUTkS4EpmSyRTI05yVfJlZyHXzUMXvS+mZd8fHsdqO2F9N6s/zYZjoRkd+5qz1MHqZZa313v1o8Dc8vxCeuqJ1d3tAf3C8ubdLo1Yc2djyK3lJt2sYLkVYzRkZcXO0OJx5LZ+JPit+13e2+oRmwDvPe0oQgfgnceg+dZWh91eSS+vQCvSwLOTAQ4Gj44evrF+5Pjm39/l+24uoJUiDQ09W6LRNIwmt45r7QaNWPO9vQU6ecnu0KJ8XVdrUFjCAJZaCT/8yRuPFFa9KFfzaj/BUTmSq5otuarbLcTTkOvmoYveF9Oy7w+PY1DYE8d1njpY561zNVqXVgqHIzIBzleERlG6ClAYE9DPHYGGZhwxzB25dWhdZ6EaJqAoofC+Lo1QeXJv8VVAFBlUUJdpIDRUDqL4dPtOA7aCvPRYC0EA+HoiFyhFVUEcKbQKGNqSYe5ohAFJUue+PLAy18BWdSbtMpxqcqVqP6mU+i/2vn5JKfW5yQ7t8hyVI7mq2ZKrut1CPA25bh666H0xLfv+8DjSyJw4rvPUwTpvnavRupz3aDSFrdAqoLSeUAeAx9qKZqypHPTzgkasiY3GuToLNbQVUQiRUnVRz0ARK4MKKorC4qu6wCilJdBQ5Kfbdw4wAcShwhioKkBVKKWovCcIIC88zlc0TEgj1gzLiiyrc18K2OgOMUGdSbsMp+otqJT6eert/QHv/ceVUovAr3vvPzvOwUxTb0HJXD10VbdbiKch181Dkrmqx3GacZ01cwVn7yF4eF2jzNXoGWaodf0sk3r9g8JSVo40NJTO0RkWdAYFtoK5Rt1XsLSeNArqnoO7GZvdnMrVy0/iiICKe72C3WFJb1AxKCyrczE3FmI+3Mm5381pGIV1kCYBn35hgecXm4SBejg2ABRhUH9KMQ1NXY/L792h8tS9Bst6f4x6Ek7S0/YW/KPe+88opb4K4L3fVkpdznRwwg6e/M14tiJpk3pDkUyNEGc369fNON9PLnpfXOT6jpoMjULgo+bMR42rsI6dQUFZuf1il0f92zNqpvyk4zDI7X5D4xeXUlbnkiPrbI2aMoeB3h9rMzaEwaFJ1t5YR3fMBgeerDVjwyC3rG0P6OSW5TRmoRnh8bxzv0PfOlbSCFtVvLeTkRcViy2Lc56tXkEcaZSG3Dry0rI1CNjNKoqiojeosM4R9TVGK+73MmIT0MsslfO8tNRgdb7Bdq+kcp5Aa5IowASKMAhIo2B/kpVGpi4oeklOO3solVL1fUJAKbVKfSfrmXJSw81pCEkeZ9rHJ4SYHfJ+cjpHBdBt5bjfzfDeP1aM8+C/J2/d7/L7d7t4HItpxM2FJo04eGR/n+Y4POhk/K9fepc317bJSs8nX1jkz37iOt9zc/6R9Y0KmNqqAqXqZQV6P2g/Cra3U8PuwPLxG/NERjEoLN9a61E6i0JzvR3zf/zzO/z+/R26A0tsApbnYsqyZKtvwXsKC85B/+QHY0+wTgBEQAE0DAQaXlhJKQqLq+o8V6MRkgSalWaMCRTX51JK7/jISpsXlhI+9/LKmavVj8Npr5b/AfgV4JpS6q8Dvwn81xMb1SU5TcPNyw5JHmfaxyeEmB3yfnI6RwXQh4Xbe1z2eDHOg68bFo7YaObTmKz0DEv72P4+zXHY6OcMioJ2I6YRaYqqYjcvHlvfqICpMcEjRTZHQftRsD2NzP64y8qzm1XERtOMI2wFW4OS3dISBYYoDKgUVBaKyqMCRRjVjwk99QRDHRjreabneu9FSVzX3ioLhyJEaY3WGu0VeI8O6myW144ARRIqbMW5G1k/rVNN57z3f08p9Qbwg9T76se9929OdGSX4DQNNy87JHmcaR+fEGJ2yPvJ6RwVQG9EGpNxZDHOg69rRHVJgcxaFtOIRmge29+nOQ4rzZg0inhvq09WeqIgoB1Hj61vVMDUVhUmqMcaBZo0qtc7CrYPCrs/7sgo2knAhzuO0llMoFlKQ9qh4U5lKcqK2AQEBiKv6GWe0lrY2+bDU8HzTNH3nlKS7X1qMYw0RVHiHTilcMqD0riqDsErp7E4stJjAs7dyPppPTHQrpRaetKLvfdb4xzMNATan5QzmJaQ5HGmfXxCiOkkH+A5v1GW6WBeqrB16BsPaWzqZs2H9qv37AfFR5mrfm7377SkoSEK9YnB9sI6vvbBFr/3/g6tKOB7X1jk+YX0kVwVsL/ssnKPNDZuJWZ/G76z3uX+TsZyK2YuDcEr5hshH3YHbHQLVtox82nIuw+6/L9/uMFaZ5fFJGa+mfDBZo/7wx5GBywnMc4rdssMW9VlHnaznMHQYUKF1p5eF/CQpKADKArY7UIGLMXQXtA0o5h2ZCidJwkDnm83ySys93poDTfmWnzsxhyxVvStw2jF9XadNwuN4eWVJjeX0oke//MG2t9g//MDvARs7/33AnAHeHnM47x0TwpCTntAddrHJ4SYPsfleuT95PR2BsX+/htNpApbP9IrBsV+QP2kXO+9nQG/f7e7n2/6vhcWH8trHVRYx++9v8U//Oo97ncHNMKQ3dzx/bc8vazk2nyy/1jOA/e2BhTO0RmULLdjXlxMub3SAuDND3f45X/+Ad08p3Ce5+caXJtPyHNHNyvJS4utYLkV852NDv/i3Q5FBeVjoyqpp0jHOXRD54jyDIMcuO+AIQAB9V0vTYfqwO9FdHnl3g6hVpTWU3lYbIaspAkfudbizlZ/OouIeu9f9t5/BPjHwI9671e898vAjwD/4CIGKIQQYnIkX/V0zlJ0+qRc7yiHNco3HZXXOrzunUF9p2uhmWCMopvVn6QbZcDqBtG+nvC5um5WZAJCPfpZPd7NXonXsNxqYCtFXjrmkojdvCSvHIutBIejV5T0co9WcFE58YB6snLUc7aicOTF3icHTUBVQe4d7cRMfxFR4LPe+/9n9IX3/h8Bf3IyQxJCCHFRJF/1dM5SdPqkXO8oh9XPC0zAkXmtw+teSOsZzk4/w1rPXBISaLWfAasbRCsK64h0HWIvbEXpRj+rx7vcClEONntDTOCJQ003K2jHIXGg2e5laDStKKQVK5yH7IKy4hX1nSt1xM+iSBNHiso5MlsRBBArzW5mL7WI6GnnnRtKqf8c+CXqyeNPApsTG9Ulk6zB9JBjIWbJLJ6vV6359riPUWTqcgajIp6jZR61T0f7elRv6vBYPvbc3P6yQq3r5snB8Rm4yGg+9eIS1jl+5+0tGlHAFz66zK2VFt5DZ1jU4wrrmk/XWgmoh5+gm29E+/mwG/MpP/rqde5sDYkDRcMELLQSFtKQtz7c4UGvZC41LMQRt1caLIQx39raQKGJAsd2z5H36xY3UQxJw7DSSMhtxWY/Y3fXUxYQBpC2oCghzyAOIUrqrzcH9QPFRWBhAeI4JNQV1huuNSPaYcJmNqByFa0oZjFt8MrzbZYaITt5xWa3oLQVC62Q28tNPvXi4qU8EoTTT67+TeDnqMsxAHxx73vPHKnvMj3kWIhZMsvn61XJV03iGBXW7Tdk7ueWyJycWdsZFMc2Z15II/q5pVdYHqzXP1e79evM3p2vg+Ne72b8g9+5y4IEc4EAACAASURBVDc+3Eaj+c6DHj/zg99NGhm+9NYGhbN0+hWf/+gKC82Q1Xaynwd7b7PH3e0B672cD7czNnoFO4OMe90hS42E5TTCOs8H3R6b3YJAKcJYs90Zcn842ppDjyyHe392LNB7fOMroHPg65zHclfbwPYOPEx0FbxJAfSIqOtdJWHOQiPjy+9s8fxizCAvWdvOKV0dkP/EjXnWewU/0YimL3M14r3f8t7/rPf+1b0/PzvuTwpOC8kfTA85FmKWyPk6/SZxjM66zJOaMx/181Fu6qh1rHWzvTpX9Sf/OlnFWjejkxU4HCutBs5XVLjHGkgPi7pkQSuOKJ2jrOoegJqANNYMS8t2VhCbEBNovNZ4q7DUOajLYAwYBWUFURjiFdhKUdq68ntsFIFSlBV0s/LSMlenunOllPr/OCJL5r3/gbGP6JJJ/mB6yLEQs0TO1+k3iWN01mUeVRsrCvRjuayDPw+DOm101DquzyWkUcTuZh+N5uWVYO97Bo1mozdEq4AAvd9Aup9b+rmlEWmSULHey+vHkMFezSsqBrljOY0Ig4APuj1s5QiUQsUaA498au8iWQte12H6oixRHkzgCY3DOyhcXVA0DGAuCac+c/VXD/x3AvwEcDllTyfsquUPptm0HYvLytPMQo5nFsY4adN2vorHnfcYndSA+ajM1WG9zLLRy8DX1dGHpWW5FbHYiEhjU/f22+tJqBS0IkPrWl0mYb5RTxCOWsfNpZQfe/V55tOQuSTgh7/vJqtzCf3c8okX5tgdWtov1pmraG+8o9xXKzYspFG9LVozKC2bvZxhblFKc20uBgXvrHXZyR8WXbi3k/HFb91lba0kiqDRBu/AllDk9d9aQ7XXAmc0WXDU7WwS6rtPqPp1o5/t8PAh4wKwMFffqVIK4igkCiCrYCGJeK7ZpN2ISCJNbAxaQ2dg+XC7R2w0f/SVFf7VT9yc7syV9/6NQ9/6LaXUP53AeKbCVckfzIJpORaXlaeZhRzPLIzxokzL+SqOd9ZjdPj8HvXiO+rrg5mrg3qZ5ctvPeCbD3qUZUVReUITEAfw8mqLW8stjFHcWe+x3i/wOCoH1+caJKFmuRUTBhoT6MfW8dZal1/8rXfZGOQEaBZbMf284kF3yFsPeiRxQLdfstxKSELFJ24scGulyfpuxlonAzzX5xusLCQUu47FZkxpPYutiPudjMI5NjPLg07Bei/nna0O33inx2C0cUdkpoBjb22NIll4jq6tsGcH2Oke/E75yH8H9ElN/Riwcr7+ZKWHzEGaaN7Zzrg23+BPNq9fyjV5qjUqpZYO/FlRSv0Z4PqExybE1LisPM0s5HhmYYxCnNfh8/tgZumor486/weFZWAr5pKIMAro5ZYogFYa0s8sw9ISGc3AVgQamnFEXtbTtySqs1HHZa7e2xqQO8f1+ZQgUHzYydjNS9AKE2gircitIw4VUVjfMasrtde9BpPIUFZ+fztGWS+N2q+LhQdb1f0GizI4VxubcVOA8+A0aKPQwd5/B9CMQ6zz3NkaXtr70WkfCx6s1G6Bd4B/50kvUEol1J8qjPfW88ve+587/1CFuDyXlaeZhRzPLIxRiPM6fH4fzCwd9fVR538aGVIT8H42pCwrWrGhqKA3KFldjWmEdQua1AT084p+XhCHCvBkhd27c6WOXMetpZRYa9Y6AwI0z88ntOOQYW6xlaNwitho8tKjcDRCQxoZwqBgWFSAZ6ER7m/HKOvl8ES6LjyKqj+p6D1EYXWuBszjVjebBu3Yv3Olfd2LsJ+XzCchLy01Lu396Im9Bfd/SanEe58d+l7svT/qZuDo5wpoeu97SqkQ+E3gZ733XznuNdPQW1CI40jm6nizMEYhzuukzNVpzv+Dmas0MgzKusHxyl4vvMP9BtPQ7FfNHLXPOW4db611+eZal9V2widuLjzWpzDUui4aqjULzWh/zKN6Wwdb9ozGoVRdxHRUwb2s6jzYILfc2+rzxbc2ubfTZTGJeXGlTaAV1nnu7+Rs9QfkVclGv0+eQRgr0jhEK4i0RitFYkIUisJVBAqiIKRfFmwNh1QVLKQxK0mDdprQSgIUGuscznlCo7m+kPB8O0HhGVpHqBWB1uwMCxyKT9yY53tfXJz4+9F5ewuOfAn4zKHvffmI7+3z9axtVOQi3Ptz8kxOnMk4/lEb1z+Mz/o/sJeVpxnXeid5fJ52jGcZ2+F/AJ6V8+1Zv36Ocp5JysHXHp4cTMrh8/ukr48a66CwhIEmDc1e2QNHaR0buxkr7YTWXi+Z0d8Ht++odYzqa3WGBYPc0o4NjTCg3NuHUAfhDzZvHoXmy8pRlI7Suf2A/Ki34WYvx2gFqm4evdJKKCvHna0eg6FlmFve3Rmws5uzWxYoPMWap9pLpn+422N9s6QsofSQRDCfQD8r6HehLCFNYWnREgWG0lXYqqKohmxuVmxbaAF5e8hGMCRNI5Ziw+p8m1AH3O/0GNqCrd4c7oV5nktjAO51M3b6lqW5kJvtBsPSsdMvuDafPM2hP7cnTq6UUteBm0BDKfUqD6vPzwEntppWSgXUjxRfAf5H7/1vP91wxUHjCBKPK4wsoebpNs3H5yxjG/3ucQUYZ9U0H59JOSkoftJ58O5G75FA9u2V1lTus8I63rrf3W/InFuPq+DdzR79wvLiYsonby7whVdWH5lYPWn7Cuv4g3s7fPGbG7y3scs317s0TEi7YfjsrSU+eq2N0YBS3FxMH2ne/KCT0WoYvr3WY7EdEWvNJ28u8O31Lv/3733IZnfIWi9ntdVgqR1yY77B5m7B1+5usrk75P52Rf+RLcw5sljonk4B9x8cuq/Sh/f6JUe1fYa9MPvu6IsCKOBhhH5PBl9/QCuAqhq1eK7PneVGwK2VJp97eZmf+mMfuZQJ1kln4p8B/ibwAvC3gb+19+evAP/ZSQv33lfe+0/vvf5zSqlPHv4dpdRPK6VeV0q9vr6+ftbxX2njCBKPK4wsoebpNs3H5yxjO6kA46ya5uMzKScFxU86Dw4Hsqd1nx1uyDzMHbl1RKauKxWbgKyo9h/hjV7zpO0rK8duVqEVhEbjnWYuCTBKs9ErAYUxAbbisebNjnpy53DMJyEOTycr2BlYtFI00wiPwhhFHBi2ehmD0mJMCCpETdn81fs6xA51mQcA5z3J3gcHNvrHppcm6ol3rrz3vwj8olLqJ7z3/+d5V+K931FK/QbwZ4FvHPrZLwC/AHXm6rzruIrGESQeVxhZQs3TbZqPz1nGdlIBxlk1zcdnUk4Kip90HoSBeiSQPa377GBD5tJZGrHGVaNHopbcViRRQBqZR17zpO0LA007CXAeSutQ2tHNKtoNxUqrTuBYW2GCR4uQjq4bHYBG08lKYq2ZTyIWUoPznv6gftRnrSevLDfmG/jdAmtL8OV+XappoVR95wrq+1sK0EqRFfUHB1aa8eWM60mBdqXUT3rvf0kp9R9xdIX2v/2E164C5d7EqgH8OvA3vPe/etxrJNB+dpK5Ov16xz2+WcvITPNxHmfmataOy8isjvskT8pGjTtzdfD1cHwA/Dzb8DTLKqxjZ1DsB9VL59joZQzyivlGyPMLaZ2/6mekoWF1rn6Mtb6bMSgsK81kvxhmL7P7jZ07WcFmL+fOeocPdgquL8R89vYqg9KyO7Q8N5cwn+7dz/FQOvfwc/97X49C7oV1/OHaDu8+qB/xJSag1YgItOfe9oCtoaVynu7A8vX313nzTp8d6nxQulfcM+akh4TnkwKxgsLXy28Aiw147lrKy0vz7PYL3u/toFXIUmxohBEfuz7Hv/7aLW6vtsY8mkedN9De3Pv7PKN7nvquV0D9+PF/f9LESpzPOMLO4wpMX0bg+7RZlXFnWmYxIzOO4zOp7T7L2J70u7N4XEaexQKkJ2WHzhoMP6j+3YfVtw8ee7v3CO2oRsfn2YanPacio7k29zD3M5oYJqFBAYPc8k/eXGOzn6O94gsfW+XWcpO3H/RweNZ2Ml69tQTAV9/b2nu057k216Cfl/yjb6xT4ijf9uwMK+5uD1lohgRe8/lXVmg3DA86GYutiO1eHfIe5bBMUH/CblBYvvZBh6/f2aHcq281n8S8+WCHrKjIywqtYbeb8+GBJ20DYLB352jck6pH1nHg9s4usDuEO+8N+J33Dmax6oGkAfzOBx2ajYh/6196eT/LdpFOeiz4P+39/V+ddcHe+68Br55zXEKcysHcRj+vP8J81BvfaX9v3Ot91kz7dk/7+K6ag9mhOkflJ3ZMDh77jd36X//59OnPg0mcU4eXudHPKV3FajuhOyzYzQs6WZ2HWm4lbPay/UzW6Ht3t/sMS0s3qyi947l2g3s7A+53c6rK8/xcg7tbfbYGBcvtaL8w6CirOCiq/X20sZuzm1W4yhOFGl0phnlF7ku8UzSikLyEsqrIp/Ry0jxsnRNHGo3ivc0Bg8JO3+RqZO8R378L3D74Gu/9X5rMsIQ4ndNmVcadabmKGRmY/u2e9vFdNReZjTp47J/U6PhpljupZs8rzZhQB6zvZmivaMcR80mEZsBmL0Oj9jNZGsVmLyPSmkZomEsCQqW53xviHDw3F3N3e8iH3SFhYFhKo0cKgx7VDDoMFO0kQAdqv0RDEECsQpT2DAsLVIQBxFOWuRo5OKy8cASJ4tZy+kiW7SKdtojol4B/Rl1WYb9j0NOE3I8imStxHheRuToqy/Gs1Vk6rUnWmTrpGJ3mGD6r2aVZdTgbBU/OQo37Og0DvV9nKo3Mue5inDUbdtw1st0r9nNVAHe2+0Q64COrrbqW1HafUGtWWjFpbOgMCj7YHmC0YrkVP2zgXFrmk4hmUjd8/vb9Ll99f5vKeW4vNUnjAK30fuZqVAQUeKQ4aRjo/e+X1vHN+zt86/4uWeFZSA23l1sU1vLWep+sqOpPdWYV3/hgg99+r4c9tN0N6k/sdc68h8+mASzGYEJIGprKOoYDSFPNc/Mxc3GDj9+Y58c+89LUZq5GUu/9fzLmMQkxFqfNaZw30zKpLMesGm3vuLNNJ2VbTpt9eRazS7PsYDZqXMf4yet6NMPVy+x+TkmjePXW0pknWAeXe9ptOFyLrRmZ/VxVnjsG1rI5KAmB774xx6deWGQ3K1nrDgFFM9K8db/P5qBgs5uxOpfw0WttXlhscGulyc6woJkYvIc76z3e+M4W73cGxFrz8nNzfP7lZfqF5dpcwnonwykI9yZpoybQo/ezonJ85TsbvPHOFt950KFXeBabEd+12ma+GbHWGfL2gw79vGK3X9E95p7MflPmCRsCw1HD6N6Be1a549vbQ+aTIbulp/DwFz8/nXWuRn5VKfWvTXQkQkypg/mIUb2Yq1SP6CiTqMt00jKvYi2oZ81lHONBYfdzSg7/SD2p8zjtNhyuxXYwV2W9ozPMaUcB7TRkd1iwM8wxRhMoTaAVhfUMy4o01JhAoTU458mK6pH6boPCspOX6EDRiAxBEFC6isJZbMV+A+ZQP94EevR+ppViNyux3qO1IQg0cRgwKC3bg5xAeVAGAoW/nKdsZ+IdaK0YFpdX5+q0k6ufpZ5gDZVSXaXUrlKqO8mBCTEtDmc5jmugepVcRA7l8DIlTzX7LuMYp5HZzykdzC6d12m34WAttsO5KqM0842Y3aJid1DSbkQsNGKsdVTeUTlPZBSNMGBQOmzlcXsThiQK9pcZBpo0MizEIa6qs1FVVRHqgEgbTMB+A+bS1U2gG5Hefw8bvZ8572knIUYpnLNUlSMvK9LQsJjGVF6Bt1B51NPNTS+E2puINqIprXN10SRzdTVMQybmrGOYVP2c05iG/XXUWODx/TCOekCTyuOI6XCRx/ioPnonPRI8bRPmzrB4rBHy4Wujn1kGpSUNDc3E0M8snawgDQ0bvYyv391hoRHx4lLKbmYxWtGIAgZFhaL+dF6nX6IUzKcxS62INDKEpu5ROCjrmldpaHh7o8vb632W04iPXGuz0kpIY4NSPLLe+TSqGzDvvS4M9H79q+9sdPn9e112egVGa/7IjTZxoHjjTodhXrIylzAsKn7rWx/yxt2LeAD4ZPPUd4hK6ieEAbCSwPJyzGdeWOHHv/8lPvXS0kTH8FSZK6XUUQ2aO8B73vsZmMeKaTENdYjOM4ajshwXYRr215PGMgooj2usJ+WlJE81+y7qGJ/nfDzNa0YNjg/W7rqxkD7WF9F7+OZal6Kq2OqVfPzGPM044NZyi/Vuxv/1L+7SzQqyzNZ5qDBgkJe8uNQkt57dYUHlPFor5hohi2nMSismDjXLrYi8sLy7OUQrjzEBvawk1Ir3N4dcX0jpF3a/OOg317r7tbEWmzFb/ZyNXs58Wn968+ZiyrCwvPnhLm/e7fL7d3eIgoDfevs+nX5BL3NY51hshfR2c6ZgXgU8HpwvgQ8y+OBuzjfv3eV+r+Cv/blPcnPpxFbIY3faM/jvAF8B/ue9P18B/j7wLaXUn57Q2MQzaBpyM9MwhtOaprFKJkrMkvOcj6d5zVF9/47qizjKes01YhwOpdj/2Vo3w7qKG0spmXN0soqVVoLzil5mCbQCpdEKjFYYrXHOYZ0n0AqFoptVaOVZmWvQ6VuGheX6QorzjmFZPTaO5VZC4Ry7eYkxmsgEKNR+/8HdrGKYO6JAAYo0CchLzW7hSJOQKAwYFo7+jFzWXsP6bsFaN7uU9Z92cvUu8Kr3/vu9998PfJq6R+APAf/NhMYmnkHTkJuZhjGc1jSNVTJRYpac53w8zWsO1u7Kijq3lEbmsdeNsl7dYY5G12UZ9n52fS7B6IB7WwMSrZlPAjZ6GVp5Womhch68w3mwzmOdQ2uN0YrKeTyeuSTAecVGd8h809CIDGs7A7TSNMLgsXGMamO14xBrHYWt8HhMUH9asJ0ENGJNUXnAM8gq4tDRjjSDrKQoKxqRpjkjl7VysNqOuD538Z8UhNPXufpd7/2nj/reUT87L8lcXQ3TkJuZhjGc1jSNVTJRYpac53w8bS21J/U1HL3uYB/AKNSP/Ozu1oD3twcsNyPCQLPWzZhLHmaiSvvwFtFokhTuvbae4Gk6e21rVpoJZeXY6OfMJyELe8s8PI40qsc6qjQ/Ws6oFlc/s3zYHbC5m5GXnmtzCR7Pmx/uEOqAV55rMygtv/HmA/7x197jbhcS6r5/5d5UoseBYpjUd3DaHF37KqCui3XUU8ZlQAWwsbeweepcVbb3mpcWNS8upcRBTLcYMrAlGk0aBiynTT7zkSV+6OM3Jv5I8GnrXH1TKfXz1I8CAf4N6keCMfVjzpkgb/zT4bSZikkWqzzNGM57voz7PJumnNFRYzm8vdMy1mfZpD84MKsuar8c7mv48Hv6se8pZR4JuY/Wg4Lrc0k9YfKw3IrqZsqVYz6pJ0evv/uAr72/y2LL8C9/7DluNNP9Za3vZmz0MkKt+bAzoJuVzCUR82mE9/XPS+v2w+/N2FCUD4upHmznMyom2hkWbPUKqsqThHXLojQy3FxMebAz5Bt3t+lmlt7QEkV1VfRdoO/rBsSRgurQ/RrH0ROrhHqidFx8axOgqidnwd56RtPNAvj2tuPudo/5qEdrDtI4wmhPFBjw0B1YNnrZpeSt4PR3rhrAXwb+OPWdzd+kzmFl1AVGx9KvcZJ3rqYpGCxOdlwhvos6buc9X67aeXbVtncaPO0+f1aP2eHtWm0nj4XMT+rgMKnm7geLDw8Ly93tAR92M+5s9FltN7i7M6AZhax3BlxfSFmei7i3mfEPf/ddOhkYDZ+4OcfP/OB38dx8g86w4De/tcmdzV2GZUW/rAi1Zr4R8fmPLtOKDev9nDsbfW4uNklCzSvX2ry30WexHRF4WGzGdLMSW1WgFEbDr/zOXR5kQ+5tDlloxbSigFBr7u70We/mdHMeq8o+rULgo9ca/I2/8OmJfmLwuDtXpzpzvPdD7/3f8t7/ee/9j3vv/6b3fuC9d+OaWE2ahG1ny3GF+C7quJ33fLlq59lV295p8LT7/Fk9Zoe366iQ+Vle/7T75bjiw8PCkZWehjGYQOO8o7SOyCicrlu6OAd3t3tUChqxwhhFNy/ZGVrKyrMzsHjvaSYhzkNeOlpxSBxq1nsF/czuLz8MwFaQW4fDMZ+E+8H2RhRgTICtoDusC4+mJkQrhVEa6zTd0uG9Rht96kdd0yAIoF84vv3gcqYop5pcKaU+ppT6ZaXUHyil3h79mfTgxknCtrPluEJ8F3Xcznu+XLXz7Kpt7zR42n3+rB6zw9t1VMj8LK8fZ3P3g8WHG5EmCRVDa7GVQytNaDSF9WgHtgSt4eZii8DDMPdY65mLQxYahjBQLKQGpRT9rEQriENNLy/JS8dqq+45OFp+WYEJIDYajaaTlfvB9mFRYW2FCWCuURceHdgS5z3WO4x2zIUapRzOupm5awVQVdCMNK9cm2xvweOc9rHgbwI/B/y3wI8CP7X32p8b52AmHWh/VnMGz6pJZq7Osv7zZK4Oh10vY3wXdb7P+nV1mQViz+tZzlyNq3HzKGR+lmvxuN8/GAg/rgjpUePe7hV0suKRJstF6biz3eP99T79oiQOQ+YbBq0Uw7xkrZfTzx2hUaztZKxtD7g2F/PHv+sat1ZadYV5Bd9a2+GtB32aUcByI2QnKzFBwAuLKXONcK8MQ90qZ64RgqqbM6eRYb4R0ckKNns5jTAABYO84p0HO9zZKkB75pKIFxcbNOOAr7y9yVv3OngcoPnmh9vcmfJnVt/7XMpf/lMf5Yc//dJE1/O0gfaG9/6fKKWU9/494L9USv0z6gnXzJCw7Wy57OP1NOvfGRT4vb8nkWmZdAPcs7js4/Q0ZrUp99Pu82k9ZuNu3AxnvxYP/35h3YmNn48ad2Ed37i7s/e6Aa/eqnM/X/n2A7747Q3eut9hd2hpJiHthuGlhRZvbXZ4/0GPnaHDAQ2jWEgj1voFvRKem+uw0ooZlJbffmeDnUGB84rr7ZhQawaVxaD5yGqb2GiWWjG9rCSNDIPCstyOWW1FtOKQr97ZIS9LNvolaWj4YHOX72z0qJwjq+Ajy02eazfo5SVfvbPB5rCu0zUrD5G/fn/Az3/xO3z85hK3Vy/+7tVpz9pMKaWBt5RS/75S6s8D1yY4LiFm1kVkWqSg53hIU+7pMsnc03lzV6dp/Hza1w0Ky25hCQONDg1Ka+JAY61nUFUUBaggIDJgFKA9KI3z4LzfLyTay+qmzK0oJNQBvczh8MQ6QGuFo75rZYK6z15uHVpBKzZkpWeznxMGilYakZcVXjkK71Hq4acgg0AzLC3r/QICTazrieMsedAt+GDncsrJn3Zy9R8CKfAfAN8P/CTwFyc1KCFm2UVkWqSg53hIU+7pMsnc03lzV6dp/Hza16WRoR3VjwddafHOkVcOYxRpEBBF4KuKwoL1gFPg64mRVmq/kGgrqZsy94qS0lW0Eo1GkbsK53xd7ymqg+re1Xkr56GXW5JQsdyMKStPb1AQhwHKayKl8N5R2AKAqnI0QsNqM4LKkdftB2fKtbmIFxYal7Lu02auXgP+GnCL+hOOAN57/33jHIwUERXT6mkaPZ/08e9x5UvGueyrZBYzV9PqNNmkk4z7vD3PtXs4d/Wk7XpSc+gHnYyNfs5KM+bafEIvq2svdQYFW72cO1t9ukPLUivij1yfxwSKb7y/xd2dnNDo/7+9Ow+S5DzvO/998qq7uvqc6bkHwAAECJAEOZYokrIp2xKPXZuSTW6I2pVpUWHuIa3ttaQNWt6w9I/s3ZCs8Fqr3Q1aJkWGZcmWpZVJmaZMHTRvCEOCJC4CmAtz9vRdd1Ve7/6RVY3qnuprpmq6uuf5RCAwlVWV+Wbn5DtvZ/7yecmlHYpph9lihulimrGMt1ZIdLHW4uZyA8dOPldtBjTaIfm0y3QxvZbNcq0kMN+9cjbWqYXVLUCadR2COKbc8JkvN2n4ETnPIZ/xODyWxrUt/vPz1/nmxVU8Fx6cLfLUKws8dbHM3kwss96kgJvsKmJDqwVuGt56Zoqf/iuPcuZwcajbv9vM1W8BPwc8y/655arUQAxioudBrXc32xjVXM2o2atJuQ+aWivcNpu0E4P+e3sn69uYu8qn+w8W+9Xjq7dDPCe7bsLkpWobYO2170esNNr8+YVlLi3XKaY9Xj3e5PsemGSpEVNpBVTbIbPFLC0/ZiKbwbEtpovptX2ZKaZ5aKbI5cUalxdqfHeuigHGcy7tyHByKofX+U53nxarPovVNofHMpyayq+tyw+TB4cKGe+2vuj6coOXbta5Vm1AZ+7D+ZrP1JiHAzx2pEgYGb56folaZ3SQseFI0cGyhForYrkeIzYUHGgbKHgOfhzz0HSBk5M50p7NYi1gqV7n+UtVKj0/49mcRdqxmSsH6wqOWsB4xuLYRJYHJnLYts1CrcFKPWAil2Ysk6HcDPHDeE/O6Z1uccEY82ljzCVjzKvd/4baMqVGxLDyS5qLUgfJTrJJ+8Fuzsut6vFt/Hks1ttrrxthxFItAIGMl8wLWGkFXC83sQRyKRcQxDbYIiCmkwm8PVsZRAYsQUTIeMlkzK3OU4Jm7TO3Tzbdu66t9nmu0qIZ+JSyHlnP5mYtqbU1kUthOxYLtYCqH2Lbr12tMQYihIofgy14NtgCsZ3U8MqkHBBJpswRYbXRwhaDiL025Ut3XXU/JuW5tw1WBLAswRhhsRliJMa2k5BZNpXcBl1u+HvWr+50cPULIvIbIvJBEfkb3f+G2jKlRsSw8kuai1IHyU6ySfvBbs7Lrerxbfx5TOVSa6+zjs1k3gUDTT+k6YcU0y5HxzLEBurtADCYSIiMASOdTODt2UrXFogNxhiafjIZc9qz17Wl32TTvevaap8PF9NkXI/Vhk/Dj5jNJ7W2luttojBmOu9S8Byi6LXq7SaGKIwpOBZEBj9KpsWxoqSGV7MdgjG4AhhDKZsmMoIx0VruScUzsAAAIABJREFUqLuunGfR9oPbbpkZII4NIoapjIMYiyhKQmaNdpJTm8h6e9av7jRz9a+B1wHP89ptQWOM+fAgG6OZK3Wv7bQujh/GrHZ+CxrLeLu63bFx0tSN29vvGZ+DsA9qMAaRudprd1Iba7N6fPPlFjcrTcbSLkc65RnKTR/XSvqBi4tVyjWffMbj6HiGsay3LkvVm7Nybeu2Pqm73SCKWay1qDQDJnMppgppRFjXpt5tl3LeWv+TDJpCuvkE17LIppL5+boZsgvzFb50foFWO2K2lGax0uLVxQYnp7I8dqTIjUqLC7fqXF+u045iMq4wO5HngaksK/WI5VqbSs0nlUlu8c2XW7iW4XXHJjg5kWMs69L2I1YbbS4uN3nm8gLldsAbjpR46+lJyn6SZ3vpRo2LSyuUUi4PzpQopD1OT+d4/GiJm5UGL81VKTcCTGw4NZ3jbQ/OMDOWHurfl7vNXL3RGPPEgNuk1J7qNxfZjdUGc+UWYNblEvww5qWblV3nSXpzKHFkmClmyKTsdbmG/TwgOajz1Kk7s1k2ab/Y+Pc5l9p+XzY7h2utkOevrzJfayPA69sRJ6dy+GHMiu/z4vUK4wWPXNqllPW4WW5xaaHKtXKLajPEETgTGWbHsqw212fAgHXtLGU85lZbye3YdsRUIQmi96vh5hNTynn4Ycwrtyp8++oqC9UWcQylrEfKtRjLOJSbITNjaZoNn2dvVHlhbplL800sE1HxY4ppj29dXeaPvutg4mR/XctQ9w2ubXFhsc6zV9OkPOHWUo1KAGEQUwuTiZgN8PXLK5yeLnBoLMNE1uPl+TJzy3XmKhGWA5XGMsuNkDcdn+TacoOnLi3QCmJu0KbWhgdnisRiyLo2T19e4dkbS1y81cCxbQ6X0sxX2nzg7EnG895tx2fYdtoLfl1EHhtqS5S6x/rNRbZZLuFO8yS93/PjmGYQHqiMlebG1EEyyL/PDT/Ej2NK2RSem9SM6s53aIn0necPy6LZjsm4NvmsS8uPKLf829q0sZ3lln9b/7RdDbcgimn6MY4IGS+5RRnGBtsSojiZj3As7bLaDqi0AlzbxbYF49gApLwkH+X7Bs9xiLGIsRAbStkUlpXMe+iITcMIti04XnLry7XBs5Nbhe3QIJK0qdaKEdvFscG1wLIsGn5EtRVQaSaTcOfSLrEt1H1DJmVhi83NSoO6H+LaLpYtuK6F5wjL9ZByy7+bvxJ3bKeDq3cA3xKRl0TkOyLyrIh8Z5gNU2rY+s1Ftlku4U7zJL3f8yyLjOscqIyV5sbUQTLIv89Zz8GzLFYbbfwgqRnVne8wNqbvPH/EMZmURTOIqDUC0p7NWNq7rU0b2zmW9m7rn7ar4ebaFhnPIjSGpp+E6x1LiGKDbSXzEZZbAaWUSzHtEkQBUWSQMAKg7Sf5KM8TWn5AbGJSjuBiEcYRrhVjWUJoIrJiiCJD6CeDjiACP0pC7ilHMCZpUz5tYaKAMIIghjiOyXo2hbRLMZNO9qcVYEWGnCc02zGRiZgtZsl5DkEUEEeGIIjxQ8NEzmEsfe+vWsHOM1cn+y0f9BODmrlS99pu5iK70zzJVpmrg0AzV+ogGeTf51or7Jtz2lgXy3OstX4Hk1yJwsBUIU0+7fRt08Zl/fqn7Wq49WauuvWwoDOo7MlcNdohFxarlKttxBKCIKQVw+xYhnzaYW6lwVLNx3UsWlHEeCbFVCGFY0tSyJSYG0s12jFkHeFmLclFPTid53Aps/YAQBDFXLhV4VatBQiHChkenMkzlvUIopiX51a5stzkUCHNickcSDKwnC6mWa37a22MjGGqmOHhQ8Wh3xK8q8yVll1QB1U3L+GHMSv15PLxZiHWO82TbPzefhuA7KQTV/vPQRwUD2Kf7jYH2T0/XMvCcy2mC+l16+v9s2sspDOnTHL7LgmoT3cyU93P9mtTv/psQZT88tat7dR9v/sz6WbIen9OM506WH4Ys1BtUW76ZF2HsayHa1mUmz6L1Ra+H5HPuKQ9h2qzDe2YtOswlU8zlvVotEOurzbJp5zX7gIEERlPyKY8xjIelVZItdGiGMY4lkXKsag2QwoZN6lgn3JoBCGXlus0WhGTOW+t/WMZj4cPlXh1sc4zr67QDELe+cgs4/kkP4bAiYkcjVyK6yuNOz5+g6I9o7rv+WHM5cVa3yD7/W5jwLeU8XomoxUeP1paC9tqoH3/OIgPIozCPnUfYPGjiOVawKNHxsil7HVt6Vd4tDtnn4G1ZZ5t7XgfNuvDgNt+Jpste+HGKl98aZGlWhOxhCeOlKgHIZVWyJdfuoUfxdTaIbmMzXLZJ+U45LMOD04WePBQgZdulJkuZgBDLuVyfr5KFEesNAJOThS4tFzBs4UXr60Sk9x+zKVdDpXSTGQ8/sLpKTwH/t03LnNxrk0MFFLwppMTfP9Dh8in4PPP3eLLF5aJYsilhe/OVfnJdzzEQq3F89crLDVafO3CIr4f49gWf/Hh6ZEPtCt1YG1XYO9+tl1wtl/YVo2+g/ggwijsU/cBlmImRUyclEPY0JZ+hUe7gfN+xUh3YrM+rN/PZLNl1VaEJVDMeFhi0wwD2oEhDGJihHTKBWMRtMGIkM+4xEYo+z6ODYjp5MUCqm0fxzFk0g5iCSEhUQSWWEQCrmvjuEJgYhyxELFoBBE3V1r4PqTcJPBujFBvhbguLNUjbtUDXBvSHhiE6yvJ9EJNP07mT4yTgH0h55L17H0RaFfqwNquwN79bLvgbL+wrRp9B/FBhFHYp+4DLJVmGwsrqTO1oS39Co92A+f9ipHuxGZ9WL+fyWbLCmmb2ECl6RObiIzjknIFx7WwMLTaAUiMmwIxhlozwBLDmOcRRuBYNmnboZRNcWIyg8RCsxViYoODg21DbGJsA0EQEQYGVyxCE2NMTNa1mR1P43nQDpLAu4ghl3YIApjM2RzKuQQRtHwQDEfH00zlUmQ8i3YYY1lJwL5aD2j40egH2u8VDbSrvbIxyA79J/Dd6SS/mwVQN4bl+wXqRy0Hs13mahTbrLZ3EI/bKOzTxsxVv4dlXNvqW+RzYzHP3u/3s7E/6obTs65zW4B+szA8JAOt7pQ93e+PZT2CMGax3uLV+Qo3qj62BVPZNNdX61xZanFqMs2TpyYJY0OzHSJiMVNMvVYMdalBzQ+ZLmbIpWzmVptcLzdZrgYYDON5j4msy8nJAsWMw3ylyQs3Kjx1foHldpvHDo/x3ieOcWQ8y1jG46nzc/zrpy5TaYX8xYdneP/ZU0mBVcuiEYQ02hENP+DKYg3HdnjLqQlOTeeHerzvtoioUgdaEv5MfsPZLLvRu7y3KN/GfEe/7wO3ZSKOlLIsVFvripj2vh6VHMzG4Gy/gP4otFPtzkE8bqOwT/0efNku19kNk3dzWON5j5Wav2X2ql8R5NWG33cbm4Xh+7XroZnia+9Va1yer/G1iyusNHyaYUQpleLZa0sYsfj6BcPFpQalTJqaH2BbcLSU41AhTSnr8NTlFWIxWDdqnDmc48ZKi1cXqtyoNSEWPNvi9EyBKDZcXmzywtwyL16psBIkt9VeXWxRaUa85w3HCMM2v/LZF5nrzN68VLtOpRXxyJExqvWQExMZ6kFI3Y+4vtKklHNpX0iC8Jq5UmoEbJbd2K4o31bf75eJ6BYU7C1iuteZEaXU4O104uRu5spCts1e7aYI8p20q/tecq1NSLkWnu1Q99u0Y8NEPoXYsFwP8KMQ104Ket5YqXNpucLFpRrtOGI6l8KyYm5VfMIowvVsxNhIp/SDI4b5WnutUGlsd2+lJm0stwJWmwHXVwOaBjwryWMFUcxCrcV4xiOIQvwoxharc2sVJrIpwoiDl7kSkeMi8mci8qKIPC8if29Y21JqkDbLbmxXlG+r7/fLRHQLCvYWMd3rzIhSavB2OnFyN3PVHWBtlb3aTRHkO2lX973kGpyhHcT4UUjOS5GyhOVaGxPBRM7Fsx1aQcRcuc5CrU2lHlJvJWUlFupt4tjiUNHDsW0CP8JIhIljgjAmNMJMPrVWqFRCiIB2BFEMhbRLKeNytOSSEfDjJI/l2hbT+TQrTR/XdvBsi8jEnQEqLDfaODYHL3MlIrPArDHmmyJSAL4B/LAx5oXNvqOZK3UvbJbL2Emeaqef2Wzy560yV8CuJowdBTvNoCl1v9iqf+lO/p51nS3zWN2r3d3Xm/UH/fKQG4uWbtdOP4hpBOHahM3dLFhvG8pNn+vLdVp+zHQxQ8P3ubLS4FAhzUMzRYI4ZqnW5tJ8jXIrpJBxyDgOJ6YyhJGhmPaYHctQbvlcW65TawRExuBYQimf5th4lnLL55UbZa5XWvzpi3MsrdQ4NFHgQ+84zcmpAkFo+MwzV/jct6/huDbvfcNR3vXEUUAYy7i4jpVcvVv7mdjMjmVGu4jonTDG3ARudv5cFZEXgaPApoMrpYZtJ3mq7vJ+k7b2K9y31fp719Gb69q4vjuZMHav7TSDptT9YrtaW/V2cvvq4nxtLU+1MWvZzU711sHyGv1zV739UTe3ZWBtcubt2rmx1la3Zt2N5QaxgGsJk/kUQRRzo9xmud5mvt7mcDHD2VNTydW2zu3JaiMgMIZ2GOL4wnjGI+06rLR9bFtYqLWS0Hw74sJSg3TKpt4MedBYvHyrQhjBYq3FVy/O8czVJFh1qVal8icv830PzXDhZpkvXCx39iDiE1+9QtOPOTM7xqOzY6zW2/hxTLkRMFlIcXw8S24Pixzfk95PRE4BTwJP9XnvIyJyTkTOLSws3IvmqPvYTvJUd5N3utP1jEKNnt3aaQZNqfvFVudxv/pW/bKW3de7rXm1mz5ku1pbfhzjWkLac2j6SQ0sR4RC2kMQWn50e30ugWOlDI/OjnF0PMtMMU0u5aytv+FHLNd9AmNwbAvPEsLY4DpCK4w6k0bbVOvJ3TQPsIGVZkgYGlbb0YZ9gKV6Mg9hFBv8OLmF6jk2rtVt2971Q0MfXIlIHvg94O8bYyob3zfGfMwYc9YYc3Z6enrYzVH3uZ3kqe4m73Sn6xmFGj27tdMMmlL3i63O4371rfplLbuvd1vzajd9yHa1tjzLIogNLT8k4yU1sEJjqLZ8DIa0Z2/6ndVGgCVQaQXU2yEWQt0PWam1afghS7U27SDCj5PbgkFoSDt2Z9LoiLF850ocSfZqPOviOEIpZW/YB5jMuTg22JbgWckg1A8jgrjbtr3rh4Za50pEXOAPgT8yxvzqdp/XzJW6F3aSubqbW1p3up5RqNGzW5q5Umq9rc7j7nu99a22qne38XN3s+3NPtu7DWBtWfeqTzei0FtDq5vN2vidctNnvtJmqpCi4YeMZ71k4mc/ZKXhk0s5LFbb5FJ2EsC3rLWaXkEcr2Wm/stLN3nhWpUHZnK864mjuI5FEMacuzTPF15cIOUKP/j6Wd58enotv9abE9sqpzZo9zxzJSIC/CvgxZ0MrJS6VzarhbPZ8t5OaGOH0y+A3pujSq7q9F/vxmKco1CjZ6PtOuvtMmhK3W+2Oo83vtfbR9yez9z9ubSb723VD3bP+656K+RmuYEYYWzCu62OV/c73eKoDf+1q2dBpzhqtdkN7CdPSmddh+6kir3950rNx7aEMEzmMgQ4UsoSRDHpqxaWWGQ8i9nxHEfHs9RaIauNnYX476Vhpr3eDvw48KyIfKuz7OeNMZ8d4jaVGqje4OfGkCfAUs1ns6KAWwVbuxO8dh+5fvLkxG0d1l4bhYlwlTqoRvX82lhYNOPZfOvqKt+9UcECHjs2xg+/6fi6p/B6vxNGEcWMR8Z1uL7a4Nmrq8Qk09UcL+WotgNepcFivcXhYoa0a631n/VWyK987kV+69y1tXV/+cIiP/NDj3J+vswv/afza8u/dH6Z/+P9bySMYL7WRoDXHylx5nBhJH6Ow3xa8MusjUuV2p96g59+HJNx7bWQJ0DGszvh0SQ82T2pe8Ol3aklek/47gSvk/k0S7UWDT8cucHVdvuglLpzo3p+9RYWNcBqI6DWjBjLuCDQaAeUW/66wdX679iAgCRB81o7Yjzv4lhCzQ9YrPuMZVwWK20mc2lKOYdyM+DGaoNGO+LKSg2LJBBugJVGwPXVFhfmG9BZHgN1P+SFm1VOTeUoZVM0g4BmMDo/x9HqzZUaMb3Bz25gEz+87cpVKeNuW0i0V3eC1+4EyFlv9E7F/RiyV2q/GNXzq7ewKBhKWYd8xubaSoAFHJ/I3laYs/udSisgimCmkAIMt1abNP2QxnLIdCHF5HiKlUZAEMY4tkUUx6zWWizWA8qtgKAdcjifIya52gVJoP1oKU3TT6YR6y7PeQ6PzRbWXbnKuM7I/BxHr0dX6h7aSaaolPEot3xef7SE2/lMNx9xpNS/6KfnJHVpNpsw1XMsnjw5sS5zNWo224demxVMVWq/GNaDJDvpW7bqI+60PXf7fc+xODWVZ7qQBpLz+qHpIhdOVIgiOD2VZzzv3RaIn86nafoRBsh6NqWMBwgnp/Ms13xmS2mm8mkm8h6VZsihYorJfIqmH2OwOVzMsFRr8d+94wRz5TLfuVrjgek8v/gjT/DokRLf88AUF+ZqfOG5OcaKaf6Hdz7A48cmCMKYQ/UWGJjq/NI7CkavR1fqHtlJ5qHWCnnu+uqm2aiNRUF79QuvbtxePp0e7E4N2FYB2e0mo1Vq1A0r97TT9e6kj9jtE8eD2J/bCh6nHSZzaQzJ7bhaK1w30fTMWJooMtiWRS7VKdPgWGQ7txYPj6U5UkracmQsS6VVodoKuVlukUvZlJshtgMmjPnkuWt86XINgGfmavzOU5f5R3/tDXz+2Sv8/jNz+MDcUot/89RVDMlVNQOd4qEtjo9nR6If0l5Q3bd2UnSvNxsVk0y2PMzt7SfbTUar1Kgb1jm5V8WE79X+9Ct0GmO4WW5wZbnOfLmFaydX5g4V0+tmwrhZblJp+tRaIZadFCY9c6jAQ9N5zhwucmk5GVh5ncT2y7dqNPyQb16tYoCClwxcluoB9TCkFZiRKh7apYMrdd/aSeZhkNmoUc1Y3KntJqNVatQN65zcq2LC92p/+hU6tRCmCimmC2nG8x4iyeCrNy4QRDG2LRQyHoGJabZigjhmLONypJRlqpDmkZkCAH6nBOfDh/JkPYc3Hy8gQNVPcleTOZec45B2ZaSKh3bpbUF13+rNPPTWsNqYnXpktrhWPA9Yq0vT/TNsnzfq5hOmC+mR+K1qEPplM/b6UrxSu7GTXOGw1tsvG3W37RnU/mxsW7/1ek7y+mgpiwj4Qcxz11epEGAhnJq8fb2ubeHZFpM5j5xbYiqfwnWstf40n3b44Ped4ptX5rm5GvHI4Qw//rYHAfiRv/AAz1wp82fP3WBqPMdH3vkgZ09P49rJhM3dPnpUal3p4Erd17on4VaTOXdzBeN5j5U5f22S0yCKN61z1WuzCY5XG/7I1La5U/0mo1ZqPxlW8d7t8oqbZaPutj13+/3N2tavYPD67YSM51KIsBZy326QBq/1vasNnzCK+Ud/8E1eXkrmEfzG9Sa/8cWX+a/fdJxyrcHvPn2DJnC1WefX/vQV/tl/U+DkVJ6FariWB9tqwup7af/26koNyHaTOXdzBRaybpLTph/vKG+kExwrpXqNcv7yTttmDKzU2yzWWqzU2/hBMki7VWklv6SGr90Z6F7l7m5LJClr88KNMgsVH0gGdgZYqCWTM3czVxkHXJL6V68uN0b2Z6mDK3Xf224y526uoDvA6k5YmvGsHeWNdIJjpVSvUc5f3mnbRGBmLM2JyRwzY2mC+LVBjx/FrDb8tQEWJP2qH8ZUmj5PX1zmwkKFpUabYiZJsnevnM2OeTg2a5mrZggBSf2rkxPZkf1Z6m1Bdd/zHIvpQnqt5lS//MN4xqMRhMzk0+TSSaE6P4wpZbef06rfpfCtcl5KqYNtq2zUXk/83m1bN0/ab939tuXaFmFoqDbbZFyHrOew2vC5udpkud4GkgHQ0fGkGGg3crHaCICY6UKGvB/y4bc9zL/8wne5WYbHj2b5228/w4nJHPm0w8XFBp/9zhUmC3n+x3c+xJnDRYCh5Obulg6u1H3PD2MWqq3knn07xHPW5x/8MOalucq6WlfA2nd84m3v8/fLQYzivGJKqXujX59wp3WqhlGva7XhJ7flqi0gyYkKMF1Ir/V9GzOq85UmfhzjWQFHShkAwsh0Cos6nfItyS+UfhRTaYWU6y2ev1FhPN8iCuBbVxa4WE7a8PT1Bv/i8y/yT9//Zr700hy/++dXWGkZrq9W+fiXLnF6usip6fzQcnN3Y7Rao9Qe2O6efb9aV6Naj0YptX/tVX2srda3MSfarXG1cVsNP6lbdXQ8h2ULi/UWQWSYyHu4lkWlGazdtnNtiygyNNoBmbTH0fEMD04WmCi4XF1prWvLyws1blaaPHutjI8hl7FwXVhuBlxbbd7Vfg6TXrlS9707mQdwVOvRKKX2r72qj7XV+lw7yUD11riqt8PbttXbT8aRYbUesFBrEUQxM4U00/nUuvjEickcADU/IAyh3A4AOD2Z5fKqv9aWo8U0bT9mKu/gAOVmjAATGZdjnatjo0iMMXvdhjVnz541586d2+tmqPvQdnmFWiu8bR7AvZ4DTCl18Ox15qrf+oBtM1fwWj8ZhDEvzZVZqPtcX27w0EyR180WbitX44dJ0P36SpNU56ns8azHJ774Ct+6tsrjR4q878lTHJnIUG+HXLi1yrnLZcazDu9+4iinpvN3vZ93S0S+YYw5u3G5XrlSitfyD92JiDd2Gvl0/8mVu0+8wO4LaI5iTkAptbfutF8YdH/Sr67VdtvKp5N+cKHaohnGhHFMyrZIuxYNP2K14VPKeusGaFnPoZR1qfshbT8mO+bwvQ9MUm7G2GKxWG8ylnMRSB4hBDKedVezZdwLo906pe6h3YRCddJipZRar7fwshUbmn5EaGCukmSjMp5NvR3eForPeQ7n56uIJXzyK+f5l1+4RL0zkPry+UX+yd98E8u1Jj//e9+hFhgs4NvXynz0PY8zM5beq93dkg6ulOroDXHW20lofbPBUu+kxUmo02z5eaWUOui6fWgp6zE7kaOUT1HIuFRbIWnHppT1qLfDdaH4ejukEYR4jsVkPs1/WKgSGUhbEBmo+BE3Ky0uLdSJMIzlXOqtgJsVn8V6WwdXSo263YRCeyctBkMp42ooXSl1X+vtQ7OeTdazcTrzCcLmofixtIdFg6VaiwenC/ypLFLvPPA4lnJ4YCpHzhVshHI9wAJmix5TudQe7en2dHClVMduJj3dzaTFuwmaashdKXUnRqHv2Kxgcr9QfHfi5+7rJ09O0PBD3nB0nOVae22C5r/7Vx/msaMlHjta4h+8q8HnnrtJMe3yI28+PjLzCPajgyuleuwmFLqTSYt3m+PSwqJKqd0apb5juyD8Zp/rhuE/+ZVX+DdP3yAErjXrfOqrl3j4cIkgivnid5c4v9AgimGleQGM8L0PTo1kPzl6LVLqANlNcT8tLKqUuhMHpe8IopgXb9aSkLuTXP2ZqwbMVVpcW21SDUJyKY+s59AKY25WWyO7rzq4UmqIdpvj0sKiSqndOih9h2tbPDqbT6YiCyEEDhdcDhfTHCtlKLgO9bZPww9JOxazhfTI7qveFlSqYxiZhd3muEZxAlKl1Gg7KH2H51h86O1nuLHc5Msv3+KBmSI/+57Xc3Qimez5x952khOv5Ml5Fm8/c4hHZosju686uFKK4WYWdp/jGs3OQik1ug5K3/GVl+f47LNztKKIxcsrfO7Z68yWsixWW3z6mZuEcUTTj3ny5AQL1RaeM5rZ1NFrkVJ74KBkFpRSaj977kaVCMNUIQMiXFpq0vBDrq02iYk5NpUnIma1GY50X62DK6U4OJkFpZTazx4/UsBGWKw2wRhOT2bIeg7HShksLK4t1rCxKGWcke6r9bagUhyczIJSSu1nb3/4MB96e5mnLy3zxuMlfuL7z3Tmds3z3jcc4rkbVV53KMebT02OdF+tgyulOg5KZkEppfYjP4z54+eu8YffuYXBsPLKEj/waIU3npjglbkK//bPrxER8/y1Mqemi5w5XNzrJm9K/yVRSiml1J4LopgrK21E4PhEjjiC8/M1AF5dbhAR88BMkYiYV5cbe9zarengSimllFJ7zrUtToynMAauLtexbHhoJg/AyYksNhYX5yvYWJzslGcYVXpbUCmllFJ7znMs/urjxxgvpLm62OCR2SJvPDEBwJnDRX7m3a/j1eUGJyeyI31LEPTKlVJKKaVGhOdYZF2HwMSsNHxqrXDtPde2SLv2yD4h2EuvXCmllFJqJHz7yjL/+A+epdoOEQM//vZTvP8tJ1mstvj1P3mFmBgLi5/6K2c4NZ3f6+ZuamjDPxH5uIjMi8hzw9qGUkoppQ6O8/M1wthwaCyD4wg3ltcXET01UyAm5tpqc6+buqVhXlv7TeDdQ1y/UkoppQ6Qh2byOJZwq9wkDA1HJtYXEb08X8XC4lgps9dN3dLQbgsaY74oIqeGtX6llFJKHSxvPDHBz77rEb51tcypyQw/8OiRtSKiP/H9p7m4WOeBqdxI3xIEDbQrpZRSakSs1Hzmyj4zxRStAIIwmTuw1gpZqvkU0g5LtfVB91G054MrEfmIiJwTkXMLCwt73RyllFJK7ZFyyycm5vBYlpiYcssHoOGHxBgm82liDA1fB1dbMsZ8zBhz1hhzdnp6eq+bo5RSSqk9Mpb2sLCYKzewsBhLewBkPQcLYanWwkLIeqNd7GC0W6fUAeWHsU4SrZRSG4znPf7SIzPcrDQZS7vk0skwJZ92ePxoiXLLZyztkU+P9vBlmKUYfhv4GvCIiFwTkZ8c1raU2k/8MOb6SoNblRbXVxr4nUyBUkrd7/wwZrXpU24EXFmuc3mxhh/Ga8tjA6tNf+T7zWE+LfjBYa1bqf0siGIMkEuXK6aJAAAK+ElEQVQ51NshQRTr1SullCLpH4PIkPFsDBBEhiBKBlL7qd8c7etqSh1Arm0hQL0dIp3XSimlkv7QtYWmHwGGUsZd6yP3U7+pgyul7jHPsTg6ntXMlVJKbeA5Fqem8kwX0kByparbR04X0jT8kKznjHy/qYMrpfaA5+igSiml+kn6R2/dMj+MWai2MCRXrzwnO9J96Oi2TCmllFKK9VlV03k9ynRwpZRSSqmRtt+yqnpbUCmllFIjbb9lVXVwpZRSSqmRt5+yqvujlUoppZRS+4QOrpRSSimlBkgHV0oppZRSA6SDK6WUUkqpAdLBlVJKKaXUAOngSimllFJqgHRwpZRSSik1QDq4UkoppZQaIB1cKaWUUkoNkA6ulFJKKaUGSAdXSimllFIDpIMrpZRSSqkB0sGVUkoppdQA6eBKKaWUUmqAdHCllFJKKTVAOrhSSimllBogHVwppZRSSg2QDq6UUkoppQZIB1dKKaWUUgOkgyullFJKqQHSwZVSSiml1ADp4EoppZRSaoB0cKWUUkopNUA6uFJKKaWUGiAdXCmllFJKDZAOrpRSSimlBkgHV0oppZRSAzTUwZWIvFtEXhKR8yLy0WFuSymllFJqFAxtcCUiNvDrwHuAx4APishjw9qeUkoptZ/4YUy9HeKH8a7eu5PPqXvLGeK6vwc4b4y5CCAivwO8D3hhiNtUSimlRp4fxlxfaWAAAY6OZ/Eca9v3droOtbeGeRSOAld7Xl/rLFtHRD4iIudE5NzCwsIQm6OUUkqNhiCKMUAu5WA6r3fy3k7XofbWMAdX0meZuW2BMR8zxpw1xpydnp4eYnOUUkqp0eDaFgLU2yHSeb2T93a6DrW3hnlb8BpwvOf1MeDGELenlFJK7QueY3F0PEsQxbi2te523lbv7XQdam8Nc3D1NHBGRE4D14EfBX5siNtTSiml9g3P2XrgtJPB0k4/p+6toQ2ujDGhiPw08EeADXzcGPP8sLanlFJKKTUKhnnlCmPMZ4HPDnMbSimllFKjRK8lKqWUUkoNkA6ulFJKKaUGSAdXSimllFIDpIMrpZRSSqkB0sGVUkoppdQA6eBKKaWUUmqAdHCllFJKKTVAYsxt0/3tGRFZAF7d63YcMFPA4l43Qg2FHtuDS4/twaXH9mA5aYy5bWLkkRpcqcETkXPGmLN73Q41eHpsDy49tgeXHtv7g94WVEoppZQaIB1cKaWUUkoNkA6uDr6P7XUD1NDosT249NgeXHps7wOauVJKKaWUGiC9cqWUUkopNUA6uNqHROTjIjIvIs/1LHujiHxNRJ4Vkc+ISLGz3BWRT3aWvygi/7DnO+8WkZdE5LyIfHQv9kWtt8tj64nIJzrLvy0i7+z5zls6y8+LyL8QEdmD3VEdInJcRP6scw4+LyJ/r7N8QkQ+LyKvdP4/3lkuneN2XkS+IyJv7lnXhzqff0VEPrRX+6QSd3BsX9c5n9si8rMb1qV98gGhg6v96TeBd29Y9hvAR40xTwD/H/BzneUfAFKd5W8B/nsROSUiNvDrwHuAx4APishj96Lxaku/yc6P7d8B6Cz/QeCfiUj3nP5/gI8AZzr/bVynurdC4GeMMY8CbwV+qnO+fRT4E2PMGeBPOq8hOS+7x+4jJMcTEZkAfgH4XuB7gF/o/qOt9sxuj+0y8HeBX+ldifbJB4sOrvYhY8wXSU7QXo8AX+z8+fPA3+x+HMiJiANkAB+okHTM540xF40xPvA7wPuG3Xa1tV0e28dIOm2MMfPAKnBWRGaBojHmayYJVX4K+OFht11tzhhz0xjzzc6fq8CLwFGSc+6TnY99kteO0/uAT5nE14FS57i+C/i8MWbZGLNC8vdBB857aLfH1hgzb4x5Ggg2rEr75ANEB1cHx3PAX+/8+QPA8c6f/z1QB24CV4BfMcYsk5z8V3u+f62zTI2ezY7tt4H3iYgjIqdJrkweJzmO13q+r8d2hIjIKeBJ4CngkDHmJiT/SAMznY9tdn7qeTvCdnhsN6PH9gDRwdXB8WGSy9HfAAokV6gg+W0oAo4Ap4GfEZEHgH4ZHH10dDRtdmw/TtIBnwP+OfBVklsUemxHlIjkgd8D/r4xprLVR/ssM1ssV3tsF8d201X0WabHdp9y9roBajCMMd8FfghARB4G/qvOWz8GfM4YEwDzIvIV4CzJb0jHe1ZxDLhx71qsdmqzY2uMCYH/pfs5Efkq8AqwQnI8u/TYjgARcUn+8f0tY8zvdxbfEpFZY8zNzm2/+c7ya/Q/P68B79yw/AvDbLfa3i6P7WY2O+ZqH9IrVweEiMx0/m8B/xvw/3beugL85c7TRzmSwOV3gaeBMyJyWkQ84EeBT9/7lqvtbHZsRSTbOaaIyA8CoTHmhc4tiKqIvLXzlODfAv7D3rReQfL0H/CvgBeNMb/a89ange4Tfx/iteP0aeBvdc7btwLlznH9I+CHRGS8E2T/oc4ytUfu4NhuRvvkA0SLiO5DIvLbJL+9TgG3SJ4eygM/1fnI7wP/0BhjOpeqP0ESfhbgE8aYX+6s570kt5Ns4OPGmF+6l/uhbrfLY3uK5B/WGLgO/KQx5tXOes6SPHmYAf4T8D8bPdn3jIi8A/gS8CzJ8QL4eZJszr8DTpD8IvQBY8xy5x/s/4skrN4AfsIYc66zrg93vgvwS8aYT9yzHVG3uYNje5jkVn6x8/ka8JgxpqJ98sGhgyullFJKqQHS24JKKaWUUgOkgyullFJKqQHSwZVSSiml1ADp4EoppZRSaoB0cKWUUkopNUA6uFJK3XMiUhKR/6nn9RER+fdD2tYPi8g/Hsa6e7bxhIj85jC3oZTaP7QUg1LqnuvU6PpDY8zj92BbXwX+ujFmcUjrd4wxoYj8MfBhY8yVYWxHKbV/6JUrpdRe+N+BB0XkWyLyyyJySkSeAxCRvy0ifyAinxGRSyLy0yLyD0TkGRH5uohMdD73oIh8TkS+ISJfEpHXbdxIZ7qgtjFmUUQKnfW5nfeKInJZRNzN1iUif01Enups+49F5FBn+S+KyMdE5D8Dn+ps7jMkVbWVUvc5HVwppfbCR4ELxpg3GWN+rs/7j5PMi/k9wC8BDWPMk8DXSKbzAfgYSeX5twA/C/zffdbzduCbAMaYKsk8fN15N38U+L3OvJubrevLwFs72/4d4H/tWfdbgPcZY36s8/oc8P07/gkopQ4snbhZKTWK/qwzGKqKSJnkqhAkU4y8oTOt09uA301migEg1Wc9s8BCz+vfIBkg/QHwE8Df2WZdx4B/25l41wMu9azr08aYZs/reeDIbndUKXXw6OBKKTWK2j1/jntexyT9lgWsGmPetM16msBY94Ux5iudW5B/CbCNMc+JSHGLdf0a8KvGmE+LyDuBX+x5r77hs+nO9pRS9zm9LaiU2gtVoHCnXzbGVIBLIvIBAEm8sc9HXwQe2rDsU8Bvk0xovt26xkgmxQb40DbNehh4brf7opQ6eHRwpZS654wxS8BXROQ5EfnlO1zNfwv8pIh8G3geeF+fz3wReFJ67vcBvwWMkwywtlvXL5LcLvwSsN3Thj8A/Mfd7oRS6uDRUgxKqQNNRP5P4DPGmD/uvH4/SRD9xwe4jRTwX4B3GGPCQa1XKbU/aeZKKXXQ/RPgewFE5NeA9wDvHfA2TgAf1YGVUgr0ypVSSiml1EBp5koppZRSaoB0cKWUUkopNUA6uFJKKaWUGiAdXCmllFJKDZAOrpRSSimlBkgHV0oppZRSA/T/Azo9h9GFXZOMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot time vs. magnitude\n",
    "_ = plt.plot(times, mags, marker='.', linestyle='none', alpha=0.1)\n",
    "\n",
    "# Label axes\n",
    "_ = plt.xlabel('time (year)')\n",
    "_ = plt.ylabel('magnitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Estimates of the mean interearthquake times\n",
    "The graphical EDA in the last exercise shows an obvious change in earthquake frequency around 2010. To compare, compute the mean time between earthquakes of magnitude 3 and larger from 1980 through 2009 and also from 2010 through mid-2017. Also include 95% confidence intervals of the mean. The variables dt_pre and dt_post respectively contain the time gap between all earthquakes of magnitude at least 3 from pre-2010 and post-2010 in units of days.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/chanseok/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:2: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time</th>\n",
       "      <th>latitude</th>\n",
       "      <th>longitude</th>\n",
       "      <th>depth</th>\n",
       "      <th>mag</th>\n",
       "      <th>magType</th>\n",
       "      <th>nst</th>\n",
       "      <th>gap</th>\n",
       "      <th>dmin</th>\n",
       "      <th>rms</th>\n",
       "      <th>...</th>\n",
       "      <th>depthError</th>\n",
       "      <th>magError</th>\n",
       "      <th>magNst</th>\n",
       "      <th>status</th>\n",
       "      <th>locationSource</th>\n",
       "      <th>magSource</th>\n",
       "      <th>loc_name</th>\n",
       "      <th>loc_admin1</th>\n",
       "      <th>loc_admin2</th>\n",
       "      <th>loc_cc</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1975-09-13 01:25:02.800</td>\n",
       "      <td>34.1390</td>\n",
       "      <td>-97.3690</td>\n",
       "      <td>5.000</td>\n",
       "      <td>3.4</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Wilson</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Carter County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1975-10-12 02:58:11.200</td>\n",
       "      <td>34.8160</td>\n",
       "      <td>-97.4060</td>\n",
       "      <td>20.000</td>\n",
       "      <td>3.2</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Maysville</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Garvin County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1975-11-29 14:29:40.900</td>\n",
       "      <td>34.5210</td>\n",
       "      <td>-97.3470</td>\n",
       "      <td>5.000</td>\n",
       "      <td>3.5</td>\n",
       "      <td>lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>slm</td>\n",
       "      <td>Wynnewood</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Garvin County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>1976-04-16 18:59:44.200</td>\n",
       "      <td>36.1070</td>\n",
       "      <td>-99.8750</td>\n",
       "      <td>5.000</td>\n",
       "      <td>3.4</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Arnett</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Ellis County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>1976-04-19 04:42:42.200</td>\n",
       "      <td>36.1340</td>\n",
       "      <td>-99.8410</td>\n",
       "      <td>5.000</td>\n",
       "      <td>3.5</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>...</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>tul</td>\n",
       "      <td>Arnett</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Ellis County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8972</th>\n",
       "      <td>2017-07-17 03:29:22.050</td>\n",
       "      <td>36.2482</td>\n",
       "      <td>-98.4378</td>\n",
       "      <td>1.990</td>\n",
       "      <td>3.2</td>\n",
       "      <td>mb_lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>24.0</td>\n",
       "      <td>0.245</td>\n",
       "      <td>0.30</td>\n",
       "      <td>...</td>\n",
       "      <td>6.8</td>\n",
       "      <td>0.051</td>\n",
       "      <td>100.0</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>us</td>\n",
       "      <td>Fairview</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Major County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8973</th>\n",
       "      <td>2017-07-17 03:51:18.700</td>\n",
       "      <td>35.8628</td>\n",
       "      <td>-96.6819</td>\n",
       "      <td>6.806</td>\n",
       "      <td>3.2</td>\n",
       "      <td>ml</td>\n",
       "      <td>NaN</td>\n",
       "      <td>42.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.37</td>\n",
       "      <td>...</td>\n",
       "      <td>3.2</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>tul</td>\n",
       "      <td>tul</td>\n",
       "      <td>Stroud</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Lincoln County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8974</th>\n",
       "      <td>2017-07-17 21:06:34.200</td>\n",
       "      <td>36.4454</td>\n",
       "      <td>-97.0626</td>\n",
       "      <td>5.966</td>\n",
       "      <td>3.3</td>\n",
       "      <td>ml</td>\n",
       "      <td>NaN</td>\n",
       "      <td>33.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.34</td>\n",
       "      <td>...</td>\n",
       "      <td>4.1</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>tul</td>\n",
       "      <td>tul</td>\n",
       "      <td>McCord</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Osage County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8975</th>\n",
       "      <td>2017-07-18 02:32:30.360</td>\n",
       "      <td>36.6221</td>\n",
       "      <td>-98.4302</td>\n",
       "      <td>5.000</td>\n",
       "      <td>3.6</td>\n",
       "      <td>mb_lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>30.0</td>\n",
       "      <td>0.240</td>\n",
       "      <td>0.56</td>\n",
       "      <td>...</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.054</td>\n",
       "      <td>90.0</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>us</td>\n",
       "      <td>Cherokee</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Alfalfa County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8977</th>\n",
       "      <td>2017-07-18 21:59:52.270</td>\n",
       "      <td>35.8749</td>\n",
       "      <td>-96.6953</td>\n",
       "      <td>4.950</td>\n",
       "      <td>3.0</td>\n",
       "      <td>mb_lg</td>\n",
       "      <td>NaN</td>\n",
       "      <td>50.0</td>\n",
       "      <td>0.089</td>\n",
       "      <td>0.30</td>\n",
       "      <td>...</td>\n",
       "      <td>3.3</td>\n",
       "      <td>0.050</td>\n",
       "      <td>106.0</td>\n",
       "      <td>reviewed</td>\n",
       "      <td>us</td>\n",
       "      <td>us</td>\n",
       "      <td>Cushing</td>\n",
       "      <td>Oklahoma</td>\n",
       "      <td>Payne County</td>\n",
       "      <td>US</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2508 rows × 26 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                        time  latitude  longitude   depth  mag magType  nst  \\\n",
       "1    1975-09-13 01:25:02.800   34.1390   -97.3690   5.000  3.4      lg  NaN   \n",
       "2    1975-10-12 02:58:11.200   34.8160   -97.4060  20.000  3.2      lg  NaN   \n",
       "3    1975-11-29 14:29:40.900   34.5210   -97.3470   5.000  3.5      lg  NaN   \n",
       "4    1976-04-16 18:59:44.200   36.1070   -99.8750   5.000  3.4     NaN  NaN   \n",
       "5    1976-04-19 04:42:42.200   36.1340   -99.8410   5.000  3.5     NaN  NaN   \n",
       "...                      ...       ...        ...     ...  ...     ...  ...   \n",
       "8972 2017-07-17 03:29:22.050   36.2482   -98.4378   1.990  3.2   mb_lg  NaN   \n",
       "8973 2017-07-17 03:51:18.700   35.8628   -96.6819   6.806  3.2      ml  NaN   \n",
       "8974 2017-07-17 21:06:34.200   36.4454   -97.0626   5.966  3.3      ml  NaN   \n",
       "8975 2017-07-18 02:32:30.360   36.6221   -98.4302   5.000  3.6   mb_lg  NaN   \n",
       "8977 2017-07-18 21:59:52.270   35.8749   -96.6953   4.950  3.0   mb_lg  NaN   \n",
       "\n",
       "       gap   dmin   rms  ... depthError magError magNst    status  \\\n",
       "1      NaN    NaN   NaN  ...        NaN      NaN    NaN  reviewed   \n",
       "2      NaN    NaN   NaN  ...        NaN      NaN    NaN  reviewed   \n",
       "3      NaN    NaN   NaN  ...        NaN      NaN    NaN  reviewed   \n",
       "4      NaN    NaN   NaN  ...        NaN      NaN    NaN  reviewed   \n",
       "5      NaN    NaN   NaN  ...        NaN      NaN    NaN  reviewed   \n",
       "...    ...    ...   ...  ...        ...      ...    ...       ...   \n",
       "8972  24.0  0.245  0.30  ...        6.8    0.051  100.0  reviewed   \n",
       "8973  42.0    NaN  0.37  ...        3.2      NaN    NaN  reviewed   \n",
       "8974  33.0    NaN  0.34  ...        4.1      NaN    NaN  reviewed   \n",
       "8975  30.0  0.240  0.56  ...        2.0    0.054   90.0  reviewed   \n",
       "8977  50.0  0.089  0.30  ...        3.3    0.050  106.0  reviewed   \n",
       "\n",
       "     locationSource  magSource   loc_name  loc_admin1      loc_admin2 loc_cc  \n",
       "1                us        tul     Wilson    Oklahoma   Carter County     US  \n",
       "2                us        tul  Maysville    Oklahoma   Garvin County     US  \n",
       "3                us        slm  Wynnewood    Oklahoma   Garvin County     US  \n",
       "4                us        tul     Arnett    Oklahoma    Ellis County     US  \n",
       "5                us        tul     Arnett    Oklahoma    Ellis County     US  \n",
       "...             ...        ...        ...         ...             ...    ...  \n",
       "8972             us         us   Fairview    Oklahoma    Major County     US  \n",
       "8973            tul        tul     Stroud    Oklahoma  Lincoln County     US  \n",
       "8974            tul        tul     McCord    Oklahoma    Osage County     US  \n",
       "8975             us         us   Cherokee    Oklahoma  Alfalfa County     US  \n",
       "8977             us         us    Cushing    Oklahoma    Payne County     US  \n",
       "\n",
       "[2508 rows x 26 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df_over3 = df[df['mag'] >= 3]\n",
    "df_over3['time'] = pd.to_datetime(df_over3['time']).copy()\n",
    "df_over3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt_pre_df = df_over3[(df_over3['time'].dt.year < 2010) & (df_over3['time'].dt.year >= 1980)]\n",
    "dt_post_df = df_over3[df_over3['time'].dt.year >= 2010]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "dt_pre = dt_pre_df.time.diff().dt.total_seconds().apply(lambda x: x / 86400).to_numpy()[1:]\n",
    "dt_post = dt_post_df.time.diff().dt.total_seconds().apply(lambda x: x / 86400).to_numpy()[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2.51464274e+02, 2.95448234e+02, 6.40863624e+02, 1.08647905e+03,\n",
       "       3.17398881e+02, 5.90184145e+02, 4.83233923e+02, 6.97192766e+01,\n",
       "       6.93095571e+02, 2.84084049e+01, 4.69432503e+02, 2.64515749e+02,\n",
       "       7.65691760e+01, 5.69709848e+01, 1.05820879e+02, 7.22962820e+02,\n",
       "       2.33607648e+02, 7.01886896e+01, 1.14955992e+02, 3.60235141e+02,\n",
       "       8.36699482e+02, 1.11743014e+02, 1.41681284e+02, 5.96914548e+02,\n",
       "       1.67977123e+02, 1.50232531e+02, 3.27134280e+02, 2.14277463e+01,\n",
       "       1.84143676e+02, 2.32951451e+02, 3.79080845e+02, 1.42725851e+02,\n",
       "       8.97876019e+01, 5.96111736e+00, 1.89721845e+01, 2.77162708e+00,\n",
       "       1.13697293e+01, 9.83503366e+01, 1.19944992e+01, 4.82748252e+00,\n",
       "       2.03827727e+01, 3.62473927e+01, 7.40873495e-01, 5.60745802e+01,\n",
       "       2.23031436e+01, 1.99971991e+00, 6.99775810e-01, 1.19004729e+01,\n",
       "       8.67183507e+00, 4.74282720e+00, 8.78014236e-01, 6.85203079e+00])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt_pre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1.21727083, 0.00597188, 1.92131238, ..., 0.7189294 , 0.22634444,\n",
       "       0.81067025])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dt_post"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1980 through 2009\n",
      "mean time gap: 204.61 days\n",
      "95% conf int: [139.84, 276.32] days\n",
      "\n",
      "2010 through mid-2017\n",
      "mean time gap: 1.12 days\n",
      "95% conf int: [0.97, 1.29] days\n"
     ]
    }
   ],
   "source": [
    "# Compute mean interearthquake time\n",
    "mean_dt_pre = np.mean(dt_pre)\n",
    "mean_dt_post = np.mean(dt_post)\n",
    "\n",
    "# Draw 10,000 bootstrap replicates of the mean\n",
    "bs_reps_pre = dcst.draw_bs_reps(dt_pre, np.mean, size=10000)\n",
    "bs_reps_post = dcst.draw_bs_reps(dt_post, np.mean, size=10000)\n",
    "\n",
    "# Compute the confidence interval\n",
    "conf_int_pre = np.percentile(bs_reps_pre, [2.5, 97.5])\n",
    "conf_int_post = np.percentile(bs_reps_post, [2.5, 97.5])\n",
    "\n",
    "# Print the results\n",
    "print(\"\"\"1980 through 2009\n",
    "mean time gap: {0:.2f} days\n",
    "95% conf int: [{1:.2f}, {2:.2f}] days\"\"\".format(mean_dt_pre, *conf_int_pre))\n",
    "\n",
    "print(\"\"\"\n",
    "2010 through mid-2017\n",
    "mean time gap: {0:.2f} days\n",
    "95% conf int: [{1:.2f}, {2:.2f}] days\"\"\".format(mean_dt_post, *conf_int_post))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hypothesis test: did earthquake frequency change?\n",
    "Obviously, there was a massive increase in earthquake frequency once wastewater injection began. Nonetheless, you will still do a hypothesis test for practice. You will not test the hypothesis that the interearthquake times have the same distribution before and after 2010, since wastewater injection may affect the distribution. Instead, you will assume that they have the same mean. So, compute the p-value associated with the hypothesis that the pre- and post-2010 interearthquake times have the same mean, using the mean of pre-2010 time gaps minus the mean of post-2010 time gaps as your test statistic.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 0.0\n"
     ]
    }
   ],
   "source": [
    "# Compute the observed test statistic\n",
    "mean_dt_diff = mean_dt_pre - mean_dt_post\n",
    "\n",
    "# Shift the post-2010 data to have the same mean as the pre-2010 data\n",
    "dt_post_shift = dt_post - mean_dt_post + mean_dt_pre\n",
    "\n",
    "# Compute 10,000 bootstrap replicates from arrays\n",
    "bs_reps_pre = dcst.draw_bs_reps(dt_pre, np.mean, size=10000)\n",
    "bs_reps_post = dcst.draw_bs_reps(dt_post_shift, np.mean, size=10000)\n",
    "\n",
    "# Get replicates of difference of means\n",
    "bs_reps = bs_reps_pre - bs_reps_post\n",
    "\n",
    "# Compute and print the p-value\n",
    "p_val = np.sum(bs_reps >= mean_dt_diff) / 10000\n",
    "print('p =', p_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In 10,000 samples, not one had a test statistic greater than was was observed. The p-value is, predictably based on what we have done so far, is tiiiiiny!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Earthquake magnitudes in Oklahoma\n",
    "- Magnitudes in Oklahoma\n",
    "    - Verify that the Gutenberg-Richter Law holds before and after 2010\n",
    "    - Compute b-values\n",
    "    - Perform hypothesis test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### EDA: Comparing magnitudes before and after 2010\n",
    "Make an ECDF of earthquake magnitudes from 1980 through 2009. On the same plot, show an ECDF of magnitudes of earthquakes from 2010 through mid-2017."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmEAAAE9CAYAAABDUbVaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df3xU1Z3/8dcnkwRQEQOIFQMJFBCQKJCAIJVVKVW7ClJoC2UVFGS12q/a1ta1rQar1l3rVtuy361KC1p+VNEW2lqVFq1gSzFBlAL+4AukiVhFjAgohsmc7x93EiaTSWYCubmTyfv5ePC4c+beufOZEcM75557jjnnEBEREZG2lRV0ASIiIiIdkUKYiIiISAAUwkREREQCoBAmIiIiEgCFMBEREZEAKISJiIiIBCA76AJaqmfPnq6wsDDoMkRERESSKi8vf885d3Kife0uhBUWFlJWVhZ0GSIiIiJJmVlFU/t0OVJEREQkAAphIiIiIgFQCBMREREJQLsbE5bI4cOHqaqq4tChQ0GXIhmsc+fO5Ofnk5OTE3QpIiKSATIihFVVVdG1a1cKCwsxs6DLkQzknGPv3r1UVVXRr1+/oMsREZEMkBGXIw8dOkSPHj0UwMQ3ZkaPHj3U2yoiIq0mI0IYoAAmvtPfMRERaU2+hTAz+7mZvWtmf29iv5nZj81su5m9amYj/aqlLVx11VX06tWLYcOGNXj+lVdeYezYsRQVFXHppZfy4YcfAt44tlmzZlFUVMSQIUP4wQ9+UP+ap59+mtNPP50BAwZwzz33JHy/RYsWsXv37vp2YWEh7733ng+f7Ijzzjsv6RxtmzZtYuzYsZxxxhmceeaZ/OpXv6rft3PnTs4++2wGDhzIl7/8ZWpqagD45JNP+PKXv8yAAQM4++yz2bVrFwA1NTVceeWVFBUVcdZZZ/H888/79dFERETanJ89YYuAi5rZfzEwMPpnHvB/fazFd7Nnz+bpp59u9PzcuXO555572Lx5M1OmTOHee+8F4PHHH+eTTz5h8+bNlJeX87Of/Yxdu3ZRW1vLddddxx/+8Ae2bt3KsmXL2Lp1a6PzxoewdHHcccfxyCOPsGXLFp5++mluvPFGPvjgAwC+/e1vc9NNN/Hmm2+Sl5fHwoULAVi4cCF5eXls376dm266iW9/+9sAPPTQQwBs3ryZ1atX841vfINIJBLMBxMRiSqvqGbBc9spr6j2943KFsGjU7ztsajcAGvv87ZHs789SuUzpcHn9i2EOedeAN5v5pDJwCPOsx44ycxO9asev40fP57u3bs3ev71119n/PjxAEycOJEnnngC8C5tHTx4kHA4zMcff0xubi4nnngiGzZsYMCAAfTv35/c3FymT5/OypUrG5xzxYoVlJWVMXPmTIYPH87HH38MwE9+8hNGjhxJUVERr732GgDvv/8+l112GWeeeSZjxozh1VdfBaC0tJQf/vCH9eccNmxYfQ/U97//fQYPHszEiROZMWNGg+Mef/xxRo8ezaBBg1i7dm2jzzto0CAGDhwIQO/evenVqxd79uzBOceaNWuYNm0aALNmzeI3v/kNACtXrmTWrFkATJs2jT/96U8459i6dSsTJkwAoFevXpx00klaLUFEmnZHTyjt5m2b8uAFcEcPb5vIPYXeOe4pTLi7vKKa3zx8J0VrZvObh+88tiC2+nb48QhvG69sEfzuBvh/a7zt0Qaxyg2weBKsucvbxgeOZPvbo1Q+U5p87iDHhJ0GVMa0q6LPNWJm88yszMzK9uzZ0ypv3la/yQwbNoxVq1YBXoCprPQ+8rRp0zj++OM59dRT6du3L9/85jfp3r07b731Fn369Kl/fX5+Pm+99VaDc06bNo2SkhKWLFnCpk2b6NKlCwA9e/Zk48aNXHvttfXB6fbbb2fEiBG8+uqr3H333VxxxRXN1ltWVsYTTzzByy+/zJNPPtko9ITDYTZs2MD999/P/Pnzmz3Xhg0bqKmp4dOf/jR79+7lpJNOIjs7u9Hniv3M2dnZdOvWjb1793LWWWexcuVKwuEwO3fupLy8vP77E5EMUtrtyJ8mj8mLHpOXeP8dPXGRwzjARQ4nDmIPXgC7yyES9rbxQeyeQjgU/TfhUHXCILZv3UPckfUw52Zt5o6sh9m37qFUPmFjq2+HF++H93d42/ggtm1l8+1U7VoLtTXgar3trrUt298epfKZ0uRzBxnCEo1ydokOdM496Jwrcc6VnHxywjUwW6S8opqZD6/nvmdfZ+bD630NYj//+c9ZsGABxcXF7N+/n9zcXMALKKFQiN27d7Nz507uu+8+duzYgXONv4JUB4R/4QtfAKC4uLi+V2vdunVcfvnlAFxwwQXs3buXffv2NXmOdevWMXnyZLp06ULXrl259NJLk75HIm+//TaXX345v/jFL8jKymr2czW176qrriI/P5+SkhJuvPFGzjnnnPoQJyJpIpVLOs2FrPjnEh6ThyPiBSwiCYNYJHIYXPQfFhdtx/vnK823D1U33waKD74AQN2P5bp2i21b1Xx7yOTm26kqPBdCuWAhb1t4bsv2t0epfKY0+dxB/otWBfSJaecDbTLIaf2OvdSEI0QcHA5HWL9jL8UFTfx2dYwGDx7Ms88+C8Abb7zB73//ewCWLl3KRRddRE5ODr169WLcuHGUlZXRp0+fBr09VVVV9O7dO6X36tSpEwChUIhwOAw0HXCys7MbjK+qm3oh0fHJ3iPehx9+yL/+679y5513MmbMGMDrpfvggw8Ih8NkZ2c3+Fz5+flUVlaSn59POBxm3759dO/eHTPjRz/6Uf15zznnnPpLnSKSBio3wMKJR9pzVkOf0Q2PKe3W4LdrK+0GpU3/IphIhAjmvODjHEQs0qgHIQKE4tqNehk+dZbXAxbbjtU5r2Hw6tz434VuI6fhdq/1PpN57aMyZJLXAxbbjlUy29tuW+kFsLp2S/UZDbNWeT09hec2/u+TbH97lMpnSpPPHWRP2CrgiuhdkmOAfc65t9vijcf070FudhYhg5zsLMb07+Hbe7377rsARCIR7rzzTq655hoA+vbty5o1a3DOcfDgQdavX8/gwYMZNWoUb775Jjt37qSmpobly5czadKkRuft2rUr+/fvT/r+48ePZ8mSJQA8//zz9OzZkxNPPJHCwkI2btwIwMaNG9m5cycAn/nMZ/jtb3/LoUOHOHDgQH1oTFVNTQ1Tpkzhiiuu4Itf/GL982bG+eefz4oVKwBYvHgxkyd7v9lNmjSJxYsXA954twsuuAAz46OPPuLgwYMArF69muzsbIYOHdqiekTkGCS7TBgbwBK18cJQgx6qoygjQl2vecN2rKpz7m5wTF27gXlroHcxZGV723lrGu6/ZdeR4NU5z2vHK5mNXfIA9ukLsEseOPpwNHE+jLsRuvf3thMTDO8omQ2X//ro36NOn9Fw7jeaDhrJ9rdHqXymNPjcvvWEmdky4Dygp5lVAbcDOQDOuf8FngI+D2wHPgKu9KuWeMUFeSyZO4b1O/Yypn+PVukFmzFjBs8//zzvvfce+fn5zJ8/nzlz5rBs2TIWLFgAeJfyrrzS+5jXXXcdV155JcOGDcM5x5VXXsmZZ54JwE9/+lMuvPBCamtrueqqqzjjjDMavd/s2bO55ppr6NKlC3/961+brKu0tLT+3Mcdd1x92Jk6dSqPPPIIw4cPZ9SoUQwaNAiAUaNGMWnSJM466ywKCgooKSmhW7dmxmnEeeyxx3jhhRfYu3cvixYtArw7OYcPH85//ud/Mn36dL773e8yYsQI5syZA8CcOXO4/PLLGTBgAN27d2f58uWAF2AvvPBCsrKyOO2003j00UdTrkNEjlGiy4QJerAcXsCq28ar66E6EqDifvu/5AFv4HlsO84rV+6k6OF+ZJsj7IzNc3dSHHdM4eeuYxdg21bhhkyi8HPXJf5c8cErXqLgFa9k9rEHI/CCV6LwJR2GJbv8lG5KSkpc/GDxbdu2MWTIkIAqyjwHDhzghBNO4KOPPmL8+PE8+OCDjBzZrqdxazX6uyYZIzZkJbo8GL2M2CBgxR0XKe3W4DKhM8iKO+bWX2/mjvLPkJUFkQjcVryOu6cUNXyvskVJL7uVV1S36i/OIm3FzMqdcyWJ9mmUszQyb948tm7dyqFDh5g1a5YCmEimSaGXKwJYTAqLWOPxK++ccAaf2r+lvpfrnRPOIH6eoakj8xlatozDNY6ckLFsZH7jelLoWSouyFP4koyjECaNLF26NOgSRCRg8ZcXE10z2f3F3/PPhy/gDKtgiysg8sXfNwphxQV5LJs3Vr1YIgkohImIZKIklxuTjeV6N7eAT31SUd/L9W5uQcKAVT53DQ8lCVjqxRJJLGMW8BYRkagkc2/V3bFIdJvojsXdM//MG7W9qXXGG7W92T3zzwnfqrggj+vOH6CQJXIU1BMmIpKBmuvpSmVOreKCPMqv/gs/02VEEd+oJ0xEJMMk6+la13Omt8s1bMdTL5eIvxTCWkFlZSXnn38+Q4YM4YwzzuCBB47Mc/P+++8zceJEBg4cyMSJE6mu9mZjfu211xg7diydOnVqsEA2wNNPP83pp5/OgAEDuOeeexK+56JFi9i9+8gCA4WFhbz33ns+fLqmzZ49u37y1WO1atWqJj/rCSeckPD5//7v/2bo0KGceeaZTJgwgYqKivp9ixcvZuDAgQwcOLB+bjSA73znO/Tp06fROW+66SaGDx/O8OHDGTRoECeddFIrfCoRn9w7yLvEeO+go3p510vv4me1l7LTncLPai+l66V3tXKBIpIKhbBWkJ2dzX333ce2bdtYv349CxYsYOvWrQDcc889TJgwgTfffJMJEybUB43u3bvz4x//mG9+85sNzlVbW8t1113HH/7wB7Zu3cqyZcvqzxUrPoSloqllhtLBpEmTuOWWW1r0mhEjRlBWVsarr77KtGnT+Na3vgV4wXf+/Pn87W9/Y8OGDcyfP78+/F566aVs2NB4jbsf/ehHbNq0iU2bNvG1r32tfo1MkbRz7yA4+I73+OA7RxXEigvyGHX1T/jD+U8x6uqfqKdLJCAKYa3g1FNPrZ9Lq2vXrgwZMoS33noLgJUrVzJr1iwAZs2axW9+8xsAevXqxahRo8jJyWlwrg0bNjBgwAD69+9Pbm4u06dPZ+XKlQ2OWbFiBWVlZcycOZPhw4fz8ccfA/CTn/yEkSNHUlRUxGuvvQZ4M+bPmzePz33uc1xxxRUcOnSIK6+8kqKiIkaMGMFzzz0HeKHu+uuvr3+PSy65hOeffx6AhQsXMmjQIM477zyuvvrqBse98MILnHPOOfTv3z9hr9iuXbsYPHgwc+fOZdiwYcycOZM//vGPjBs3joEDB9YHotj337lzJ2PHjmXUqFF873vfa/J7P//88znuuOMAGDNmDFVVVQA888wzTJw4ke7du5OXl8fEiRN5+umn64879dT4e7waWrZsGTNmzGj2GJHA1AWwptop0qVGkeB13BBWuQHW3udtW9GuXbt4+eWXOfvsswF455136v/RP/XUU+vXkmzKW2+9RZ8+R9Y1z8/Prw90daZNm0ZJSQlLlixh06ZNdOnSBfAWyd64cSPXXnttg0uc5eXlrFy5kqVLl9YvobR582aWLVvGrFmz6hfvTmT37t18//vfZ/369axevbo+3NV5++23WbduHb/73e+a7Mnavn07N9xwA6+++iqvvfYaS5cuZd26dfzwhz/k7rsbr+92ww03cO211/LSSy/xqU99qtnvq87ChQu5+OKLgdS+w6ZUVFSwc+dOLrjggpSOF0lHdWPAYpcKEpH00zFDWOUGWDwJ1tzlbVspiB04cICpU6dy//33c+KJJx7VORItI2WWaBafxuouoRUXF7Nr16765ydNmlQf1NatW8fll18OwODBgykoKOCNN95o8pwbNmzgX/7lX+jevTs5OTkNFuUGuOyyy8jKymLo0KG8807i38j79etHUVERWVlZnHHGGUyYMAEzo6ioqEGddV588cX6nqi6Wpvzy1/+krKyMm6++Wbg2L7D5cuXM23aNEKhUPKDRdLUYxdvprbWC2G1tV5bRNJPxwxhu9ZCbQ24Wm+7a+0xn/Lw4cNMnTqVmTNnNhhPdMopp/D2228DXq9Rr169mj1Pfn4+lZWV9e2qqip69+6dUg2dOnUCIBQKNRj/dfzxx9c/bmqt0OzsbCKRI78v1/WOJVtbtO49mzs29pisrKz6dlZWVpPj1BKFpu985zv1g+fr/PGPf+Suu+5i1apV9ec9lu9w+fLluhQp7d6W3fsYEF5K/5qlDAgvZcvuBGtDikjgOmYIKzwXQrlgIW9beO4xnc45x5w5cxgyZAhf//rXG+ybNGlS/d15ixcvZvLkyc2ea9SoUbz55pvs3LmTmpoali9fzqRJkxod17VrV/bv39/iWsePH8+SJUsAeOONN/jHP/7B6aefTmFhIZs2bSISiVBZWVk/Vmv06NH8+c9/prq6mnA4zBNPPNHi92ypcePGsXz5coD6WgHuuuuu+sHzAC+//DL//u//zqpVqxqE2wsvvJBnn32W6upqqqurefbZZ7nwwguTvu/rr79OdXU1Y8eObeVPJNK29uz/pNm2iKSHjjlZa5/RMGuV1wNWeK7XPgYvvvgijz76KEVFRfW9NHfffTef//znueWWW/jSl77EwoUL6du3L48//jgA//znPykpKeHDDz8kKyuL+++/n61bt3LiiSfy05/+lAsvvJDa2lquuuoqzjjjjEbvOXv2bK655hq6dOnCX//615Rr/epXv8o111xDUVER2dnZLFq0iE6dOjFu3Lj6y4bDhg2rv9HgtNNO49Zbb+Xss8+md+/eDB06lG7duiV5l2PzwAMP8JWvfIUHHniAqVOnNnnczTffzIEDB+ovkfbt25dVq1bRvXt3vve97zFq1CgAbrvtNrp37w7At771LZYuXcpHH31Efn4+c+fOpbS0FPAG5E+fPj3lS5civkiy3BAkX3JIRNoHS3a5Kd2UlJS4srKyBs9t27aNIUOGBFRR5jtw4AAnnHAC4XCYKVOmcNVVVzFlypSgywqE/q6Jr+KXG4JGQSxS2g1zYOaN+XIGWXHH3PrrzSz92z/q2185uy93TynypWQRaZ6ZlTvnShLt65iXI6VFSktLGT58OMOGDaNfv35cdtllQZckIs2YOjKf3JBhQG7ImDoyP+iSRCSBjnk5UlokfkZ/EQlOmBA51NZPPxEmRG7cMcUFeSybN5b1WvdRJK2pJ0xEpB1ZcfEmDtWGcA4O1YZYcfGmhMdpMlaR9JcxPWHOOQ2oFl+1t/GTkpm27N7HreFH69szNf2ESLuVET1hnTt3Zu/evfpHUnzjnGPv3r107tw56FKkg9P0EyKZIyN6wvLz86mqqmLPnj1BlyIZrHPnzuTna4CziIi0jowIYTk5OfTr1y/oMkRERERSlhEhTESkXUhhIlYR6TgyYkyYiEjai5+INdHErCno2bVTs20RaT8UwkRE0oiL28Yb1rtbs20RaT8UwkRE0kSYEDhvOSJctB2n+qOa+vUis6JtEWmfFMJERNLE7Wf9mZraLJyDmtosbj/rz42OGdO/B51ysggZ5OZkMaZ/jwAqFZHWoIH5IiJtxAEWs423Z/8nnB7+ZX37cwnmACsuyGPJ3DFakkgkAyiEiYi0gQhgMSksYo0vRaQ66L64IE/hSyQD6HKkiEiamDoyn9yQYUBuyJg6UpMDi2Qy9YSJiKSJ4oI8ls0bq0uNIh2EQpiISBrRpUaRjkMhTESktdw7CA6+A8efAje/EXQ1IpLmNCZMRKQ13DsId/AdHOAOvuMFMhGRZiiEiYi0gsjBd8BFp55w0baISDMUwkRE2kAkunWuYVtEOi6FMBGRNpBduo9aA2dQa15bRDo2DcwXEWkjdcFLv/2KCOhngYiIiEggFMJEREREAuBrCDOzi8zsdTPbbma3JNjf18yeM7OXzexVM/u8n/WIiIiIpAvfQpiZhYAFwMXAUGCGmQ2NO+y7wGPOuRHAdOB//KpHREREJJ342RM2GtjunNvhnKsBlgOT445xwInRx92A3T7WIyIiIpI2/Lw78jSgMqZdBZwdd0wp8KyZfQ04Hvisj/WIiIiIpA0/e8IswXMurj0DWOScywc+DzxqZo1qMrN5ZlZmZmV79uzxoVQRERGRtuVnT1gV0CemnU/jy41zgIsAnHN/NbPOQE/g3diDnHMPAg8ClJSUxAc5ERHfhUu7kYU3070mWhWR1uBnT9hLwEAz62dmuXgD71fFHfMPYAKAmQ0BOgPq6hKRtBIu7UbIgTkIOa8tInKsfAthzrkwcD3wDLAN7y7ILWZ2h5lNih72DeBqM3sFWAbMds6pp0tE0krdD0qzhu1Gx9QNwjBNwigiyfm6bJFz7ingqbjnbot5vBUY52cNIiJtYs5qbOFEIJrF5qwOtBwRSX9aO1JEpDX0Ge0Fr11rofBcry0i0gyFMBGR1tJntMKXiKRMwxZEREREAqAQJiIiIhIAhTARERGRACiEiYiIiARAIUxEREQkAAphIiIiIgFQCBMREREJgEKYiIiISAAUwkREREQCoBAmIiIiEgCFMBGRR6bAnad4WxGRNqIQJiId2yNTcDvW4MKHcDvWKIiJSJtRCBORDi2yYw04MAAXbYuItAGFMBGRJCLRrXMN2yIix0IhTEQkiezSfdQaOINa89oiIscqO+gCRETag7rgpd9cRaS16OeJiIiISAAUwkREREQCoBAmIiIiEgCFMBEREZEAKISJiIiIBEAhTERERCQACmEiIiIiAVAIE5GMVl5RzYLntlNeUR10KSIiDWiyVhHJWOUV1cx4aD2HwxFysrNYdvUYigvygi5LRARQT5iIZLAnN1ZRE47ggJpwhCc3VjU6Jguiq3d7W/1QFJG2op83IpKx9uz/pNl2HYvbioi0BYUwERERkQBoTJiIZLTt2V8hKwsiEfgqfwq6HBGReuoJE5GM9T87JhAKgRmEQl5bRCRdKISJSMaq+wFn1rAtIpIO9DNJRDq2UKfm2yIiPlEIE5GO7XvvHgleoU5eW0SkDWhgvoiIgpeIBEA9YSIiIiIBUAgTERERCYBCmIiIiEgAFMJEREREAuBrCDOzi8zsdTPbbma3NHHMl8xsq5ltMbOlftYjIiIiki58uzvSzELAAmAiUAW8ZGarnHNbY44ZCPwHMM45V21mvfyqR0QyU3lFNet37GVM/x4UF+QFXY6ISMr8nKJiNLDdObcDwMyWA5OBrTHHXA0scM5VAzjndJ+4iKSsvKKaGQ+t53A4Qk52FsuuHqMgJiLthp+XI08DKmPaVdHnYg0CBpnZi2a23swu8rEeEckwT26soiYcwQE14QhPbqwKuiQRkZT52RNmCZ5zCd5/IHAekA+sNbNhzrkPGpzIbB4wD6Bv376tX6mItEt79n/SbDsLcHU/iUx3IolIevHzZ1IV0CemnQ/sTnDMSufcYefcTuB1vFDWgHPuQedciXOu5OSTT/atYBHJPBa3FRFJF36GsJeAgWbWz8xygenAqrhjfgOcD2BmPfEuT+7wsSYRERGRtOBbCHPOhYHrgWeAbcBjzrktZnaHmU2KHvYMsNfMtgLPATc75/b6VZOIZJaeXTs12xYRSWe+LuDtnHsKeCruudtiHjvg69E/IiItMqx3t2bbIiLpTONURaTd2rJ7X7NtEZF0phAmIu1WsrsjRUTSmUKYiIiISAAUwkREREQCoBAmIiIiEgBf744UEfFTz66d2J79FbKyIBKB27quC7okEZGUqSdMRNqtO175DKEQmEEo5LVFRNoLhTARSVvlFdUseG475RXVCffX/QAza9gWEWkPmr0caWbZ0ZnvRUTaVHlFNTMeWs/hcISc7CyWXT2G4oK8lp2kcx4cqm7YFhFJE8l+cdxQ98DMfuJzLSIi9Z7cWEVNOIIDasIRntxY1fKT3LLrSPDqnOe1RUTSRLKB+RbzeJyfhYiIxEplItYswNX9lLImfqtU8BKRNJWsJ8y1SRUiIkfJ4rYiIu1Fsp6wwWb2Kt7Pt09HHxNtO+fcmb5WJyIdVs+unZpti4i0d8lC2JA2qUJEJM6w3t2abYuItHfNhjDnXAWAmZ0EDIw+/YZzbp/fhYlIx7Zl975m2yIi7V2yKSpygQeBy4CdeJchC8zs18A1zrka/0sUkY4olYH5IiLtWbKB+d8FcoA+zrkRzrnhQF+88PY9v4sTERERyVTJQtgXgKudc/vrnog+/iowxc/CRERERDJZshAWcc59FP+kc+4Amr5CRERE5KgluzvSmVkeiafgifhQj4iIiEiHkCyEdQPKSRzC1BMmIr7RPGEikumSTVFR2EZ1iIg0oHnCRCTTNTsmzMwuNLNpCZ7/iplN9K8sEcl05RXVLHhuO+UV1Qn3a54wEcl0yS5HzgcuTfD8GuDXwOpWr0hEMl55RTUzHlrP4XCEnOwsll09huKCvAbHxI930PgHEck0ye6OPM45tyf+SefcP4Hj/SlJRDLdkxurqAlHcEBNOMKTG6saHTN1ZD65IcOA3JAxdWR+m9cpIuKnZD1hnc0s2zkXjn3SzHKALv6VJSKZLJXZ8IsL8lg2byzrd+xlTP8ejXrKRETau2Qh7EngITO73jl3EMDMjgd+HN0nIuKb4oI8hS8RyVipLFv0DlBhZuVmVg7sAvZE94mItJimnxARST5FRRi4xczmAwOiT293zn3se2UikrE0/YSISPIpKr4FEA1dg51zm+sCmJnd3Qb1iUgG0vQTIiLJL0dOj3n8H3H7LmrlWkSkg0hlYL6ISKZLFsKsiceJ2iIiIiKSomQhzDXxOFFbRERERFKUbIqKs8zsQ7xery7Rx0TbnX2tTERERCSDJbs7MtRWhYhIx5HyFBWlMXdNlmrwvohklmSXI0VEWl1KU1SUdmu+LSLSzimEiUib0xQVIiIKYSLig/KKahY8t53yiuqE+1ttiorcrs23RUTSWLKB+SIiLVJeUc2Mh9ZzOBwhJzuLZVePabT+Y6stW3RrFdydDzX7vQB2a9XRli0i0uYUwkSkVT25sYqacASAmnCEJzdWNQphU0fms6KsksO1jpyQMaGdhoIAABRiSURBVHVk/tG/oYKXiLRTCmEi0qpSudRYXJDHsnljWb9jL2P692gU0kREOgJfx4SZ2UVm9rqZbTezW5o5bpqZOTMr8bMeEfFfqpcaiwvyuO78AQpgItJh+RbCzCwELAAuBoYCM8xsaILjugL/B/ibX7WISNtJafoJERHxtSdsNLDdObfDOVcDLAcmJzju+8B/AYd8rEVE2oimnxARSY2fIew0oDKmXRV9rp6ZjQD6OOd+52MdItKGWm36CRGRDOdnCLMEz9Uv+m1mWcCPgG8kPZHZPDMrM7OyPXv2tGKJIiIiIsHwM4RVAX1i2vnA7ph2V2AY8LyZ7QLGAKsSDc53zj3onCtxzpWcfPLJPpYsIiIi0jb8DGEvAQPNrJ+Z5QLTgVV1O51z+5xzPZ1zhc65QmA9MMk5V+ZjTSIiIiJpwbcQ5pwLA9cDzwDbgMecc1vM7A4zm+TX+4qIiIi0B75O1uqcewp4Ku6525o49jw/axGRttFqSxKJiGQ4LeAtIi3W3ALdmidMRCQ1WrZIRFok2QLd1R/VYHi3QmdF2yIi0ph6wkSkReoW6HYcWaA71pj+PeiUk0XIIDcnizH9ewRTqIhImlNPmIi0SLLJWIsL8lgyd4wW5xYRSUIhTERaXXFBnsKXiEgSuhwpIiIiEgCFMBFpEU1BISLSOhTCRKRFNAWFiEjrUAgTkRbZsntfs20REUmNQpiItEiyuyNFRCQ1CmEiIiIiAVAIExEREQmAQpiIiIhIABTCRKRFNEWFiEjrUAgTkRbRFBUiIq1DIUxEGiivqGbBc9spr6hOuF9TVIiItA6tHSki9corqpnx0HoOhyPkZGex7OoxjdaAdHGviW+LiEhq1BMmIvWe3FhFTTiCA2rCEZ7cWNXomKkj88kNGQbkhoypI/PbvE4RkUygnjARqZfKRKzFBXksmzeW9Tv2MqZ/j0Y9ZSIikhqFMBFpseKCPIUvEZFjpBAmIvXadPqJ0pi7Kks1uF9EOh6NCRORem02/URpt+bbIiIdgEKYiNTT9BMiIm1HIUxE6qUyML/NZB/XfFtEpJ1TCBOR9PTdt48Er+zjvLaISAbRwHwRSV8KXiKSwdQTJiIiIhIAhTARqdemU1SIiHRwCmEiHUiyxbnbbIoKERHRmDCRjiKVxbmrP6rB8Bblzoq2RUTEH+oJE+kgUlmce0z/HnTKySJkkJuTxZj+Pdq+UBGRDkI9YSIdRKqLcy+ZO0aLc4uItAGFMBFpQItzi4i0DV2OFBEREQmAQpiIiIhIABTCRDoIzQEmIpJeFMJEOgjNASYikl4UwkQ6iC279zXbFhGRtqUQJtJBpDJFhYiItB2FMJEMkWxJIhERSS++zhNmZhcBDwAh4GHn3D1x+78OzAXCwB7gKudchZ81iWSiVJYk0sB8EZH04ltPmJmFgAXAxcBQYIaZDY077GWgxDl3JrAC+C+/6hHJZKksSTR1ZD65IcOA3JAxdWR+m9cpIiJH+NkTNhrY7pzbAWBmy4HJwNa6A5xzz8Ucvx74Nx/rEclYLkkbvJnwl80bqyWJRETShJ8h7DSgMqZdBZzdzPFzgD/4WI9Ixkp1+gktSSQikj78DGGW4LlEv6BjZv8GlAD/0sT+ecA8gL59+7ZWfSIZQ9NPiIi0P37eHVkF9Ilp5wO74w8ys88C3wEmOecS3jPvnHvQOVfinCs5+eSTfSlWpD3T9BMiIu2PnyHsJWCgmfUzs1xgOrAq9gAzGwH8DC+AvetjLSIiIiJpxbcQ5pwLA9cDzwDbgMecc1vM7A4zmxQ97F7gBOBxM9tkZquaOJ2IiIhIRvF1njDn3FPAU3HP3Rbz+LN+vr+IiIhIutKM+SLtRHMz4msiVhGR9sfXnjARaR3JZsSfOjKfFWWVHK515GgiVhGRdkEhTKQdqJsRH47MiB8bwjQRq4hI+6MQJtIOpDIFhSZiFRFpXzQmTERERCQACmEiIiIiAVAIExEREQmAQphIO6ApKEREMo9CmEg7MKx3t2bbIiLS/iiEibQDW3bva7YtIiLtj6aoEEkD5RXVzc7xlcoUFWmlcgPsWguF53rtusd9Rgdbl4hIGlEIEwlYstnwoZ2NCavcAAsnHmmHOkEkDKFcmLVKQUxEJEqXI0UCVjcbvuPIbPjxpo7MJzdkGJCb7ssSxQYwgNpPwNVCbY3XIyYiIoB6wkQCl+ps+O16WSILeT1hdZcnRUREIUykvWjXyxJd8J2GY8JK90FpzB2epbrRQEQ6HoUwkYC1q/FeR+vcbzR+TsFLRDo4hTCRgLXLOcDuHQQH34HjT4Gb3/CeK1sE21YGWpaISHuiECYSsHY3B1hdAANve+8gOP9W+N0NwdYlItLO6O5IkYC1uznA6gJYbFs9YCIiLaYQJuKz8opqFjy3nfKK6qBL8c+QyUFXICLS7uhypIiPMm4i1qaUzPa221Z6gSz20qQG4IuIJKQQJuKjuolY4chErPEhbOrIfFaUVXK41pGTDhOxlp4EOO/xnNUNZ7iv3ND060pmHwljdVsREWmSQpiIj1ySNqTZRKyxAQy82e/rgljlBlg8KbDSREQyjUKYiI9SnX4ifSZiTRATd631Qtiutd7SQyIi0io0MF/ER+1u+olE6pYaKjzXW3qoEWvTckREMoV6wkR8lHbTTxzNUkF1Y8L6jIZZq7wesT/dEd1pUPpBq5YoItJRKISJdBSl3Rq3Y4NYc4Pu6/QZ7f1JtAyRiIi0iC5Hiohn19qgKxAR6VDUEyaSKe7oCZHDYFlw1TNH7mh88X54b3vy19eN/Yp1/CmtX6eIiADqCRM5Jslmw2+ziVjrAhiAi3hTS5Qtgl9cDK/9Ht57Pfk5+oz2pqPIOd5rxy7OLSIirU49YSJHKZXZ8NtsIta6ABZr20qIhFt2nj6j4Tu7W6cmERFplkKYyFFKZTb8Vp2INf7OxsoN3jiuLj0SHz9kMux8oZkgpqklRESCpBAmcpRSnX6iVSZiTXRnY3YXCH8CRBK/pmQ2nDL0yJiwngPhtd9Fd2pqCRGRoCmEiaSD2JDVtTd0L4TPzvfaryxN/JraGpoMYHX6jIbpTbxeREQCpRAm0ozyimr/13SM7+Xav9v78/OLvDsdE433Am/2+qZ6wnoXt3qZIiLSuhTCRJpQXlHNzIfXUxOOkJudxZK5DQfep3znY2zIKrkSMDhrBryz1Rs83xRX6/1pSt3s9V16wMd74ZVfwfv/Dz51Fsxbk8InFBGRICmEiTRh/Y691IQjRJw38H79jr0NQljdnY83soSLQy/ROTIFyobA+v+BQ/sguxN8UNHwpGW/8Lbli5sPWAAWaronrFvfI7PX19Es9iIi7YpCmEgT8o7LJeK8xxHntWMVF+SxNWcGoWjbNv8vbE7x5MkCWME5DceE7XnDuxsyctgLYDel+kYiIpKuFMKkw0o23mvL7n38NftaTsnax4eRznz4fD94ORvCNfV3Gh71/0AWajqI9b8Arvj1kXZsb5eIiGQMhTDpkFKZaPXmv0+iW8hb4Lpb6BDdPt4GH0d3pjIDfVNKrmo4JmzIZNi6Ev7xF+h7TsMAJiIiGUshTDJSsl6uJzdWsZXpZOVCJAK3b1zX6LgTa98HwAycS3Fq03E3wut/ODIm7FNFMXNz4U2yWqfPaG8uLziyFRGRDsOcc/6d3Owi4AEgBDzsnLsnbn8n4BGgGNgLfNk5t6u5c5aUlLiysjJ/CpZj0ibTOQDceSqEP4Ls4+C7byes46yHC8kybyzXK3N3NaonXNqNUMxf/VqD7NiABFDaDYcXvuq2ScWfQ0REOjQzK3fOlSTa59sC3mYWAhYAFwNDgRlmNjTusDlAtXNuAPAj4D/9qidlq2+HH4/wtnJE5QZYe5+3TaC8opryh7/G55/7POUPf63JBa15ZArceYq3bUrpSd60DqUnNd5XF8DA2955aqNDzvpFIaEsrwcrlOW149X9xTdr2I5ncVs+VQQ9T4fBl3iLXTeoWwFMRERS5+flyNHAdufcDgAzWw5MBrbGHDMZKI0+XgH81MzM+dk915zVt3tLvMCR7cT5gZSSVio3wOJJ3gztoVxvfqq4weJu9e1cnfVbAK6231K+ujvM/XHD8zwyBXZE56/ascZrx49/Kj0Jr98Jb1t6UsPldeoCWFNtGgYs5xIHrCzAxSSshCFszmpYOLFhO36QvIKXiIgcJd96woDTgMqYdlX0uYTHOOfCwD6g0WrEZjbPzMrMrGzPnj0+lQtsW9V8u6PatdYLYK7W2+5a2+iQYfv+DBzpWaprN/CPvzTfBo4EsKbayWUBWPSVTQWsSx5o2Mt1yQONj+kz2gteE25LHMBERESOgZ8hLNEQmvh/UVM5Bufcg865Eudcycknn9wqxSU0ZFLz7Y6q8FyvB8xC3rbw3EaHdD7zsgbBp/OZlzU+T99zmm8TfXFz7f4XNN8GKN2HRV9p0XYjJbO94PXpC7xtUwPj+4z2JkFVABMRkVbm28B8MxsLlDrnLoy2/wPAOfeDmGOeiR7zVzPLBv4JnNzc5UjfB+avvt3rARsySZciY1Vu8HrACs9tOpCk8t09MiX5VAz1lySt4aXIlpxDREQkDTQ3MN/PEJYNvAFMAN4CXgK+4pzbEnPMdUCRc+4aM5sOfME596Xmzqu7I0VERKS9aC6E+TYw3zkXNrPrgWfwpqj4uXNui5ndAZQ551YBC4FHzWw78D4w3a96RERERNKJr5O1OueeAp6Ke+62mMeHgC/6WYOIiIhIOvJzYL6IiIiINEEhTERERCQACmEiIiIiAVAIExEREQmAQpiIiIhIABTCRERERAKgECYiIiISAN9mzPeLme0BKnx+m57Aez6/R0em79c/+m79pe/XP/pu/aXv1z/JvtsC51zCha/bXQhrC2ZW1tQSA3Ls9P36R9+tv/T9+kffrb/0/frnWL5bXY4UERERCYBCmIiIiEgAFMISezDoAjKcvl//6Lv1l75f/+i79Ze+X/8c9XerMWEiIiIiAVBPmIiIiEgAFMLimNlFZva6mW03s1uCrieTmNnPzexdM/t70LVkGjPrY2bPmdk2M9tiZjcEXVOmMLPOZrbBzF6Jfrfzg64pE5lZyMxeNrPfBV1LJjGzXWa22cw2mVlZ0PVkGjM7ycxWmNlr0Z+/Y1v0el2OPMLMQsAbwESgCngJmOGc2xpoYRnCzMYDB4BHnHPDgq4nk5jZqcCpzrmNZtYVKAcu09/dY2dmBhzvnDtgZjnAOuAG59z6gEvLKGb2daAEONE5d0nQ9WQKM9sFlDjnNEeYD8xsMbDWOfewmeUCxznnPkj19eoJa2g0sN05t8M5VwMsByYHXFPGcM69ALwfdB2ZyDn3tnNuY/TxfmAbcFqwVWUG5zkQbeZE/+i311ZkZvnAvwIPB12LSKrM7ERgPLAQwDlX05IABgph8U4DKmPaVegfMmlnzKwQGAH8LdhKMkf0Utkm4F1gtXNO323ruh/4FhAJupAM5IBnzazczOYFXUyG6Q/sAX4RvZT+sJkd35ITKIQ1ZAme02+80m6Y2QnAE8CNzrkPg64nUzjnap1zw4F8YLSZ6XJ6KzGzS4B3nXPlQdeSocY550YCFwPXRYeFSOvIBkYC/9c5NwI4CLRoLLlCWENVQJ+Ydj6wO6BaRFokOl7pCWCJc+7JoOvJRNFLDc8DFwVcSiYZB0yKjl1aDlxgZr8MtqTM4ZzbHd2+C/wab9iNtI4qoCqmZ3wFXihLmUJYQy8BA82sX3SA3XRgVcA1iSQVHTy+ENjmnPvvoOvJJGZ2spmdFH3cBfgs8FqwVWUO59x/OOfynXOFeD9z1zjn/i3gsjKCmR0fvVGH6GWyzwG6O72VOOf+CVSa2enRpyYALboZKrvVq2rHnHNhM7seeAYIAT93zm0JuKyMYWbLgPOAnmZWBdzunFsYbFUZYxxwObA5OnYJ4Fbn3FMB1pQpTgUWR++ezgIec85pGgVpD04Bfu39jkY2sNQ593SwJWWcrwFLoh03O4ArW/JiTVEhIiIiEgBdjhQREREJgEKYiIiISAAUwkREREQCoBAmIiIiEgCFMBEREZEAKISJiMQws2vM7Iro49lm1vsozrHLzHq2fnUikkk0T5iISAzn3P/GNGfjTW6plTNEpNWpJ0xE0pqZFZrZa9HFcf9uZkvM7LNm9qKZvWlmo6N//hJdRPcvdTNYm9lxZvaYmb1qZr8ys7+ZWUl03wEzu8vMXjGz9WZ2SvT5UjP7pplNA0rwJmLcZGZdYnu4zKzEzJ6PPu5hZs9G3/9nxKxDa2b/ZmYbouf4WXTSVxERhTARaRcGAA8AZwKDga8AnwG+CdyKt4zQ+OgiurcBd0df91Wg2jl3JvB9oDjmnMcD651zZwEvAFfHvqFzbgVQBsx0zg13zn3cTH23A+ui778K6AtgZkOAL+MtojwcqAVmHtU3ICIZR5cjRaQ92Omc2wxgZluAPznnnJltBgqBbnhLCw0EHJATfd1n8MIbzrm/m9mrMeesAeqWHyoHJh5DfeOBL0Tf5/dmVh19fgJe8HspunRMF+DdY3gfEckgCmEi0h58EvM4EtOO4P0c+z7wnHNuipkVAs9H9xtNO+yOrNtWS2o/D8McuYLQOW5fojXgDFjsnPuPFM4tIh2MLkeKSCboBrwVfTw75vl1wJcAzGwoUNTC8+4Husa0d3HkkubUmOdfIHqZ0cwuBvKiz/8JmGZmvaL7uptZQQtrEJEMpRAmIpngv4AfmNmLQOzA9/8BTo5ehvw28CqwrwXnXQT8b93AfGA+8ICZrcXrPaszHxhvZhuBzwH/AHDObQW+CzwbrWE1cOpRfD4RyUB2pDdeRCSzRO9EzHHOHTKzT+P1TA1yztUEXJqIiMaEiUhGOw54zsxy8MZnXasAJiLpQj1hIiIiIgHQmDARERGRACiEiYiIiARAIUxEREQkAAphIiIiIgFQCBMREREJgEKYiIiISAD+P0gqiSB+aCnFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Get magnitudes before and after 2010\n",
    "mags_pre = mags[times < 2010]\n",
    "mags_post = mags[times >= 2010]\n",
    "\n",
    "# Generate ECDFs\n",
    "_ = plt.plot(*dcst.ecdf(mags_pre), marker='.', linestyle='none')\n",
    "_ = plt.plot(*dcst.ecdf(mags_post), marker='.', linestyle='none')\n",
    "\n",
    "# Label axes\n",
    "_ = plt.xlabel('magnitude')\n",
    "_ = plt.ylabel('ECDF')\n",
    "plt.legend(('1980 though 2009', '2010 through mid-2017'), loc='upper left');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quantification of the b-values\n",
    "Based on the plot you generated in the previous exercise, you can safely use a completeness threshold of mt = 3. Using this threshold, compute b-values for the period between 1980 and 2009 and for 2010 through mid-2017."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "def b_value(mags, mt, perc=[2.5, 97.5], n_reps=None):\n",
    "    \"\"\"Compute the b-value and optionally its confidence interval.\"\"\"\n",
    "    # Extract magnitudes above completeness threshold: m\n",
    "    m = mags[mags >= mt]\n",
    "    \n",
    "    # Compute b-value: b\n",
    "    b = (np.mean(m) - mt) * np.log(10)\n",
    "    \n",
    "    # Draw bootstrap replicates\n",
    "    if n_reps is None:\n",
    "        return b\n",
    "    else:\n",
    "        m_bs_reps = dcst.draw_bs_reps(m, np.mean, n_reps)\n",
    "        \n",
    "        # Compute b-value from replicates: b_bs_reps\n",
    "        b_bs_reps = (m_bs_reps - mt) * np.log(10)\n",
    "        \n",
    "        # Compute confidence interval: conf_int\n",
    "        conf_int = np.percentile(b_bs_reps, perc)\n",
    "        \n",
    "        return b, conf_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "mt = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1980 through 2009\n",
      "b-value: 0.75\n",
      "95% conf int: [0.56, 0.94]\n",
      "\n",
      "2010 through mid-2017\n",
      "b-value: 0.62\n",
      "95% conf int: [0.60, 0.65]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Compute b-value and confidence interval for pre-2010\n",
    "b_pre, conf_int_pre = b_value(mags_pre, mt, perc=[2.5, 97.5], n_reps=10000)\n",
    "\n",
    "# Compute b-value and confidence interval for post-2010\n",
    "b_post, conf_int_post = b_value(mags_post, mt, perc=[2.5, 97.5], n_reps=10000)\n",
    "\n",
    "# Report the results\n",
    "print(\"\"\"\n",
    "1980 through 2009\n",
    "b-value: {0:.2f}\n",
    "95% conf int: [{1:.2f}, {2:.2f}]\n",
    "\n",
    "2010 through mid-2017\n",
    "b-value: {3:.2f}\n",
    "95% conf int: [{4:.2f}, {5:.2f}]\n",
    "\"\"\".format(b_pre, *conf_int_pre, b_post, *conf_int_post))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The confidence interval for the b-value for recent earthquakes is tighter than for earlier ones because there are many more recent ones. Still, the confidence intervals overlap, and we can perform a hypothesis test to see if we might get these results if the b-values are actually the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Hypothesis test: are the b-values different?\n",
    "Perform the hypothesis test sketched out on the previous exercise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 0.0817\n"
     ]
    }
   ],
   "source": [
    "# Only magnitudes above completeness threshold\n",
    "mags_pre = mags_pre[mags_pre >= mt]\n",
    "mags_post = mags_post[mags_post >= mt]\n",
    "\n",
    "# Observed difference in mean magnitudes: diff_obs\n",
    "diff_obs = np.mean(mags_post) - np.mean(mags_pre)\n",
    "\n",
    "# Generate permutation replicates: perm_reps\n",
    "perm_reps = dcst.draw_perm_reps(mags_post, mags_pre, dcst.diff_of_means, size=10000)\n",
    "\n",
    "# Compute and print p-value\n",
    "p_val = np.sum(perm_reps < diff_obs) / 10000\n",
    "print('p =', p_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A p-value around 0.1 suggests that the observed magnitudes are commensurate with there being no change in b-value after wastewater injection began."
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}