oklahoma_fracking_earthquakes.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![bubble_map](data/oklahoma_bubbles.jpg)\n",
    "# Oklahoma: Variations in earthquake frequency and seismic activity\n",
    "\n",
    "The frequency of earthquakes and seismicity can change to do events other than natural causes (tectonic plate movement and stress release). For example, the state of Oklahoma began  natural gas fracking in 2010. The fracking itself, however, does not seem to change seismic behavior. It is the high pressure injection of waste-water deep into the earth that may cause the seismic activity (i.e. there may be a correletion between the injection wells as earthquakes.\n",
    "\n",
    "We will perform a quantitative analysis of seismic activity from 1980 to mid-2017 to investigate the effect that increased injection of saline wastewater due to oil mining in Oklahoma has had on the seismicity of the region.\n",
    "\n",
    "*image courtesy of <a href=\"http://www.behindenergy.com/confirmed-oklahoma-earthquakes-caused-by-fracking/?lang=en\">behindenergy.com</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import the packages\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import dc_stat_think as dcst\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## EDA\n",
    "Make a plot of all Oklahoma earthquakes between 1980 to the first half of 2017"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = r'/statistical_seismology/data'\n",
    "# Import the seimology data for Oklahoma: oklahoma\n",
    "oklahoma = pd.read_csv('data/oklahoma_earthquakes_1950-2017.csv')\n",
    "#print(oklahoma.head())\n",
    "from datetime import datetime as dt\n",
    "oklahoma['time'] = pd.to_datetime(oklahoma['time'])\n",
    "time_mags = oklahoma[(oklahoma['time'] >= '1980-01-01 01:00:00.0') & (oklahoma['mag'])]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mags: (8924,)\n",
      "time: (8924,)\n"
     ]
    }
   ],
   "source": [
    "# Create arrays of the times of earthquaks: time\n",
    "time = pd.read_csv('data/oklahoma_decimal_dates.csv', index_col=0)\n",
    "time = np.array(time['0'])\n",
    "\n",
    "# Create an array of the magnitudes of the earthquakes: mags\n",
    "mags = pd.read_csv('data/oklahoma_magnitudes.csv', index_col=0)\n",
    "mags = np.array(mags['0'])\n",
    "print('mags:',mags.shape)\n",
    "print('time:', time.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAK4CAYAAACMHlu3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3VuMHeld9/tfPXVcxz642257PB7Hk52EyYte2AQ2w0HsICFx84oIEkGEkgskLriAC5ASIqEoQogQoQjuOFxwkVwEoQRFvGgrNwmvFGXPBva7gTckmRxm4szBY7vb7u51qvNT+2J1L3e3V9ttu5erV/f3I41mdXWtp/5VT3Wt/vmp6sepqqoSAAAAAAA1MHUXAAAAAAA4uwilAAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUxqu7gF3r6/2ZtLu01NTm5mgmbWO26Lv5RL/NL/puftF384l+m1/03fyi7+q1utqZuvzUj5R6nlt3CXhM9N18ot/mF303v+i7+US/zS/6bn7RdyfTqQ+lAAAAAICTi1AKAAAAAKjNTJ8p/au/+it95StfUZ7n+uAHP6gPfOADs9wcAAAAAGDOzCyU/vM//7P+7d/+TZ/73OcUx7H+5m/+ZlabAgAAAADMKaeqqmoWDX/605+W4zj67ne/q8FgoI985CP64R/+4UPXL4qSB48BAAAA4IyZ2Ujp5uambty4ob/8y7/UG2+8od/6rd/Sl770JTmOc8j6s/nTzKurnZlNN4PZou/mE/02v+i7+UXfzSf6bX7Rd/OLvqvXYVPCzCyULi4u6tq1awqCQNeuXVMYhrp7967OnTs3q00CAAAAAObMzP767o/92I/pq1/9qqqq0q1btxTHsRYXF2e1OQAAAADAHJrZSOl73/te/eu//qve//73q6oqffzjH5fr8swoAAAAAOCemU4J85GPfGSWzQMAAAAA5tzMbt8FAAAAAOBhCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAAAAAFAbQikAAAAAoDaEUgAAAACogbWV4rSQtVXdpdTKq7sAAAAAADhrrK10/WZfWV4q8F1dXevIGKfusmrBSCkAAAAAPGVpXirLS0lSlpdKd16fRYRSAAAAAHjKQt9V4LuSpMB3Fe68Pou4fRcAAAAAnjJjHF1d6yjNS4W+e2Zv3ZUIpQAAAABQC2McNUIiGbfvAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAMApYm2lOC1kbVV3KcCRMCkOAAAAcEpYW+n6zb6yvFTgu7q61pExTt1lAQ/ESCkAAABwSqR5qSwvJUlZXirdeQ2cZIRSAAAA4JQIfVeB70qSAt9VuPMaOMm4fRcAAAA4JYxxdHWtozQvFfout+5iLhBKAQAAgFPEGEeNkF/zMT+4fRcAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAAAAAFAbQikAAAAAoDaEUgAAAABAbQilAAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAAAAAFAbQikAAAAAoDaEUgAAAABAbQilAAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAAAAAFAbQikAAAAAHANrK8VpIWurukuZK17dBQAAAADAvLO20vWbfWV5qcB3dXWtI2OcusuaC4yUAgAAAMATSvNSWV5KkrK8VLrzGg9HKAUAAACAJxT6rgLflSQFvqtw5zUejtt3AQAAAOAJGePo6lpHaV4q9F1u3X0EhFIAAAAAOAbGOGqERKxHxe27AAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKjNTP801Pve9z51Oh1J0uXLl/XJT35ylpsDAAAAAMyZmYXSNE0lSZ/97GdntQkAAAAAODWsrc7kPKczu3335ZdfVhzH+o3f+A19+MMf1r//+7/PalMAAAAAMNesrXT9Zl8/uNnX9Zt9WVvVXdJT41RVNZO9/fa3v63/+I//0Ac+8AFdv35dv/mbv6kvfelL8rzpg7NFUcrz3FmUAgAAAAAn2ijJ9cob25Ovn7+8oGbk11jR0zOz23ff9ra36bnnnpPjOHrb296mxcVFra+v6+LFi1PX39wczaSO1dWO1tf7M2kbs0XfzSf6bX7Rd/OLvptP9Nv8ou/m10nuO2srjYapsrxU4Lvqb8ca9pO6yzpWq6udqctndvvu5z//ef3Jn/yJJOnWrVsaDAZaXV2d1eYAAAAAYG4Z4+jqWkfPrXV0da1zpp4pndlI6fvf/3597GMf0wc/+EE5jqM//uM/PvTWXQAAAAA464xx1AjPXmaa2R4HQaBPf/rTs2oeAAAAAHAKzOz2XQAAAAAAHoZQCgAAAACoDaEUAAAAAFAbQikAAAAAoDaEUgAAAABAbQilAAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAACcIkVhdbeXqChs3aUAR+LVXQAAAACA41EUVi9986bipFAj8vTiC2vyPMahcLJxhgIAAACnRG+UKU4KSVKcFOqNsporAh6OUAoAAACcEt1moEY0vhmyEXnqNoOaKwIejtt3AQAAgFPC84xefGFNvVGmbjPg1l3MBUIpAAAAcIp4ntFyN6q7DODI+KcTAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA1TwgAAAADAI7C2UpwWkqTQd5WXVqHvyhin5srmE6EUAAAAAI7I2kqv3ujp9dt9lVWlwHO1shApCjxdXesQTB8Dt+8CAAAAwBGlealhkisvrNK0VH+YKS+ssrxUmpd1lzeXGCkFAAAAgCMKfVetyJfvGRnXUeC58j2jwHcV+m7d5c0lQikAAAAAHJExjq5d6uriuaYknik9DoRSAAAAAHgExjhqNfzJ157HU5FPgqMHAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAHCmWFspTgtZW9VdCgBgj5N2fT5p9ZxmXt0FAADwtFhb6frNvrK8VOC7urrWkTFO3WUBwJl30q7PJ62e046RUgDAmZHmpbK8lCRleal05zUAoF4n7fp80uo57QilAIAzI/RdBb4rSQp8V+HOawBAvU7a9fmk1XPacfsuAODMMMbR1bWO0rxU6LvcigUAJ8RJuz6ftHpOO0IpAOBMMcZRI+TjDwBOmpN2fT5p9Zxm3L4LAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAE+RtZXitJC1Vd2lzERd+3fc253W3t5lu6+zrNTdXqKisE9U42k/Lx6EiXcAAACAp8TaStdv9pXlpQLf1dW1joxx6i7r2NS1f8e93WntSZos81wjOVKSFPrOm1tqhp5aDV8vvrAmz5s+7vegGk/7efEwjJQCAAAAT0mal8ryUpKU5aXSndenRV37d9zbndbe3mXDJNcwzjXKCg3jXGVZKU4K9UbZY9V42s+LhyGUAgAAAE9J6LsKfFeSFPiuwp3Xp0Vd+3fc253W3t5lrchXq+GrGYxHSF3XUSPy1G0Gj1XjaT8vHobbdwEAAICnxBhHV9c6SvNSoe+euls069q/497uYe3tXSaNRzjf8eyiBkmubjM49Nbdh9V42s+LhyGUAgAAAE+RMY4a4en9Nbyu/Tvu7U5r7+Cy3dfLwdFGNh9U42k/Lx6E23cBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKgNoRQAAAAAUBtCKQAAAACgNoRSAAAAAHPH2kpxWqgorOK0kLXVsbY7rb1p23zQ+jiaszkRDgAAAIC5ZW2l6zf7StJCd/qJznUjRYGnq2sdGeM8cbtZXirw3X3tTdtm4LmSIxWFvW99HB0jpQAAAADmSpqXyvJSeWkVJ4XywirLS6V5eSztSrqvvWnbHCa5hnE+dX0cHSOlAAAAAOZK6LsKfFfWVmpEnnzPKPBdhb57LO3ujpTubW/qNg+MlD7p9s8qQikAAACAuWKMo6trHaV5qXe4i8pLq9B3n/jW2b3tHmzvsG1Kmro+jo5QCgAAAGDuGOOoEY7jjOcd31OJe9s96jYPWx9HwzOlAAAAAIDaEEoBAAAAALUhlAIAAAAAakMoBQAAAADUhlAKAAAAAKgNoRQAAAAAUBtCKQBAklTaSnFayNqq7lIAAHPMzuDzxNpKwzjXMM6P1O5hNWSF1cZWrP4wm3zvSeud9v7D6i2mbP9xzeI414UJdQAAsrbSK29s6dbtvgLf1dW1DhOAAwAembWVrt/sK8vLY/s8sbbSqzd6ev12X5WkKxc6unaxe2i7h9VQFFb/4/99TS9//45c19GPvn1Vb7vY1Wu3B49d77RtSZpar7WVvvaNt3RjfTjZ/vPPLDzW8ZnFca4TI6UAAKV5qTQrJUlZXirNy5orAgDMozQvleXH+3mS5qWGSa68sCoKq2GcP7Ddw2rojTJtDzIVpVWaldocpOqNsieqd9q2Dqu3N8o0GOb7tv+4x2cWx7lOhFIAgELfVRi4kqTAdxX6bs0VAQDmUei7Cvzj/TwJfVetyJfvGXmeUavhP7Ddw2roNgMttAN5rlEYuFpqh+o2gyeqd9q2Dqu32wzUbvn7tv+4x2cWx7lO3L4LAJAxjp6/vKjIjD/o5vkWIABAfYxxdHWtozQvj+3zxBhH1y51dfFcU5LUCL0HtntYDZ5n9H++54qePdccB8eG/8T1Hvb+afUa4+in331RW4N03/YfxyyOc50IpQAASZJrHDVCPhYAAE/GzODzxBhHrYb/xDUEntHKYuNI6z7Jtg6r15uy/ePc7rzi9l0AAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1IZQCgAAAACoDaEUAPDErK0Up4WsreouBQCeqnm9/s1D3XtrnFbvcezDbhtFYae2P4xzDeP8vhqKwmpjK1Z/mB1a35Pu8+N8f16djoltAAC1sbbS9Zt9ZXmpwHd1da0z95N4A8BRzOv1bx7q3luj5xrJkYrCTuqV9MT7sLuNJC10p5/oXDdSFHiT9l+90dPrt/uqJF1ebcsYR0VhZRxHN+4O9dbGUK7r6L9eW5HrmX31Pc7xfFi/zEO/PS5GSgEATyTNS2V5KUnK8lLpzmsAOO3m9fo3D3XvrXGYjEcrpXv1Hsc+7LaRl1ZxUigv7L72h0muvLAqCqutQTqpYWuYaquXqiit0qzU+nZ8X31Pus/T2pmHfntchFIAwBMJfVeB70qSAt9VuPMaAE67eb3+zUPde2tsRb5aDV/SvXqPYx922/Bdo0bkyffMvvZbkS/fM/I8o8V2OKlhsRVqsRvKc43CwNXqQuO++p50n6e1Mw/99ri4fRcA8ESMcXR1raM0LxX67qm5lQgAHmZer3/zUPfBGiXdV++T7sPebbzDXVRe2n1tXbvU1cVzTUlSI/T21XDtUldbg3QcXncC6ZMez4f1yzz02+MilAIAnpgxzuQDGwDOknm9/s1D3QdrPFjvcezD3jY8z9z3vd3AebAGYxytLDamfu+46nmc788rbt8FAAAAANSGUAoAAAAAqA2hFAAAAABQG0IpAAAAAKA2hFIAAAAAQG1mGkrv3Lmjn/u5n9Mrr7wyy80AAAAAAObUzEJpnuf6+Mc/riiKZrUJAAAAAAdYWylOC1lb1V3KsZq2X4ctG8a5hnEua6t9X2eF3bf+wXUlqSis7vYSFYXd1/7Bbe1+Xexp81GO/VHaOytmNsnNpz71Kf3ar/2a/vqv/3pWmwAAAACwh7WVrt/sK8tLBb6rq2sdGePUXdYTm7ZfkqYue/VGT6/f7quSdHm1LUl6c32gsqq0fKOvyJWiwNOV821dv9mfrHvlQkdXVtv655dvKU4KhaGrS8stlbaS5xrJGQfWwHd15Xxbr90eKEkL3eknOteNFHjuvnUedOwP7s+09qLAOzX99zAzCaV///d/r+XlZf3sz/7skUPp0lJTnufOohytrnZm0i5mj76bT/Tb/KLv5hd9N5/ot/l1UvtulORq9jM1d77uLDTUjPxaazoO0/ZL0tRl3p2RwiiQJFU7gS6MAmV5qd4w1fLFrqLAk9/w5YXeZF0v8GQ9I9fz1G57yvJS1jha7DaUZIUkqd0exye/4avZCmV8V1txoWYrnNS6u86Djv3B/ZnWXhR4p6b/HmYmofQLX/iCHMfRSy+9pG9961v66Ec/qr/4i7/Q6urqoe/Z3BzNohStrna0vt6fSduYLfpuPtFv84u+m1/03Xyi3+bXSe47ayuNhulkBK6/HWvYT+ou64lN2y9JU5cVaaE0yVRJcjrjwJkm2XikdLGl0TCVzUvlLX/fukUWyhRWZVFMRkqNrbS1NZqMlG6NMgW+q7zlazRMlaSFyqLQaJhORkp313nQsT+4P9Pas3l5avpv12H/mONUVTXTm5U/9KEP6ROf+ISef/75B643qx/sk3zRwIPRd/OJfptf9N38ou/mE/02v05631lbKc1Lhb57qm79nLZfhy2L0/HIZiMcj8Htfn3p0qJu3tyerH9wXWMcFYVVb5Sp2wxkjDNpX9K+be1u23eN8tJOXedR9mdae6ep/6TDQ+nMnikFAAAA8PQZ40zC2Gkybb8OW9Zq7L/ldffrwDP71p+2rucZLXfv/bHWvesffO/u155npq7zKPtzWHtnwczP1s9+9rOz3gQAAAAAYE6drQgOAAAAADhRCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqc/r+LBeAU+O0/kl7zC/OSQDHbRbXldN2rdo7bUvou/umXxnGudK81GI7lOeZfdO57P0LtkVhtTVI5TcCFYVVXlq5jqPeKJPvGnk7f5XX2kpbg1Sh76rV8Kcev4PHN8tKbfRirXQbCgL30Br2vk/SZPqXOC327cNZRCgFcCJZW+n6zf5kUumra51T8cGK+cU5CeC4zeK6ctquVdZWevVGT6/f7qusKgWeq5WFSIHnqiitvv7qhvKy0uXVtn7inef1r9+5rTgp1Ig8vfjC2iSofu0bb+nG+lCd9i11Ql9LnVDffXNLcVpolBa6dnFBl1fburk50lsbQ7muox99+6qef2Zh3/E7eHwvLTf1f/3LDzSMc7Uavn7xPVf0P7+3fl8Ne9/nuUZypCwrdWt7pNt3Y5V2vA8/9e61MxlMz94eA5gLaV4qy0tJUpaXSndeA3XhnARw3GZxXTlt16o0LzVMcuWFVZqW6g8z5YXVMMm1sR0ryUqVpVV/mOnG3aHiZDyiGieFeqNMktQbZRoMcxWlVX9U6O4gVW+UqTfMlKSlkrRUPx63t9VLVZRWaVZqc5Ded/wOHt8bd4caxrmk8ajtD273p9aw933DJNcwzpWXVtu9THFSTPZhd/2zhlAK4EQKfVfBzu0tge9ObnUB6sI5CeC4zeK6ctquVaHvqhX58j2jMHTVaQXyPaNW5GtloaEocOW6Rp1WoEvLLTWi8Y2gjchTtxlIkrrNQO2WL8816jQ9LbdDdZuBuq1AUegqCl11GuP2FruhPNcoDFwttcP7jt/B43tpuaVWw5cktRq+njvfmVrD3ve1Il+thi/fNVroBmpE3mQfdtc/a5yqqqq6i5Ck9fX+TNpdXe3MrG3MFn03n46z307bMzEnHT9zD3dSz0n6bj7Rb/PrpH/WndRr1eM6zmdKL15cUNxPeKa0JqurnanLeaYUwIlljKNGyGUKJwfnJIDjNovrymm7VhnjTEYjJe0Lbp1WoL0xx/OMlrvRfW14ntHKYkOL7VB5nE3aWAka921rZbFx3/sPrrP3+AaBq0sr7YfWcPB9u6873v59OItOdxQHAAAAAJxohFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGpDKAUAAAAA1Ob0/K1oAAAwN07bPIo4W8qdeTNP6vl71J+vvfN/NkKvtn05WO+0+Tx35xR90Pyhku6bI3T362bgqTfKlJTS5t2BPM9oqRPJ84ysre6b73RvXb5rlJdWvmuU5qWk8Xypu69910xqC31XgyS/b47So+z/3u3szsW6dz+n1XlarqWEUgAA8FRZW+n6zb6yvFTgu7q61pnrX6Zwtlhb6ZU3tnTrdv9Enr9H/fmyttKrN3p6/XZflaQrFzq6drH71PflYL1Xzrf12u2BsryU5xrJkZKk0Lff2FSSlfI9ox99+6qef2bhvlqLwuqlb95UnBRqRJ5+/B3n9a/fua3+KNOb6wOlWantUaayqLS62NALV5f14rvXdP1mX//xyrrystLl1bZ+6t1rMsbR9Zt9JWmhO/1ES+1Qd/qJisKqkuS5RmVpVVaVeoNsHJyN5BlX7aavVsPXiy+sPTSY7u7/3u1sDlKd60aKAm/Sf9ZW+t4b2/vq/MkfuqA3Noan4lpKKAUAAE9VmpfKdkYYsrxUmpf7JpQHTrI0L5VmJ/f8PerPV5qXGia58sJK0mQE7mnvy8F6e6Ns8vUwySfr9IaZHI3D2eYgnVprb5QpTsYjv3FS6MbdoeKkUJZZbQ8zZUWlOClVWqthWujuINVGL9bWMFWy06f9YabeKFMj9JTlpfLSKk7Go+KDYT7Z1qDI5XtGeWG11U/luUajspScQo3QU5wU6o0yLXejI+3/3u3ESaG8aWWce/2X5uV9dW704lNzLZ3PqgEAwNwKfVeB707+dX/3Fj1gHoS+q3Cc407k+XvUn6/Qd9WKfPmeUSWp1fBr2ZeD9XabgbYG42DainzJkVzHUbcVTEZKl9rh1Fq7zUCNyJuMlF5abunNjaEKa7XQCpRmparKqiyMWqGn5XaolW5Dg1GhKHCVl5U6rUDdZiBjHAW+K2srNSJPzdBTu+VPRkpbTV9laeUHRos2VJqXikIjz7hyXUeNyFO3GRx5//duJ83H+7m3/0Lf1WIr3FfnSrehJLOn4lrqVFVV1V2EJK2v92fS7upqZ2ZtY7bou/lEv80v+m5+zWPfnZbnoJ7EPPYbxpbPtfXmja0Te/7yTOnhz5R2F5o8U1qj1dXO1OWMlAIAgKfOGGdubzMD3BN+/h7158sYR62G/xQqengde+s9+PXu65Wg8dC2PM/su2V279dR5Gl1taPowICiMY46rUAH49LeOnYD5t6guff13tqWg0cbsXzYdh5U52m5ljIlDAAAAACgNoRSAAAAAEBtCKUAAAAAgNoQSgEAAAAAtSGUAgAAAABqQygFAAAAANRm/v9+MHBKzdu8UwAAzKtH/cwtd+b3nDav5kn4zD44V+fTNm2e0aMem6Kw2hqkU+cj3btfxjhT5+zcO++qJMVpoSQptDlMdXmlfd+27vYSOY600ApVVpVcx9HdXqI4zRV4roLA1VJnPKXMbl2N0FN/mOlOP9Gl5ZbSvNT6dqxznUjGONoepmqEnpY6keK0UD/OdH6xqWDKVDFZVur21kit0FcQjOc5fdhcpXuP7UmYZ/Y4EEqBE8jaStdv9pXlpQLf1dW1zlxfaAAAOKke9TPX2kqvvLGlW7f7CnxXV8639drtwYn5zC4Kq5e+eVNxUqgReXrxhbWnGkz3Hk/PNZIzrukox6YorL72jbd0Y30o13X0o29f1fPPLMgYZ99+haGr892G/vP6HeVlpcurbf3kD13Qa7cHev12X5Wky6vjAPrKW1v6ny+vyzWOzi829dsf/LHJtr76nzf07etbsrI6v9DU/3Z5Qd95Y1Ov3Ohps5/IM64urbT0rueWZBxHN++OZIy02An1n6/eVZoVClxHlYz6cSbXGYfDUVqqGfm6cr6p7UGh3FpdPNfSf3vx6r5gmmWl/vtL39dbG0NlRaV3X11SI/K13Am1OUh1rhspCrzJcTt4bG1V6c31gSpJVy50dO1id25/XySUAidQmpfK8lKSlOWl0rw8FRMjAwBw0jzqZ26al0qze+v3RtmJ+szujTLFyXj0LE4K9UaZlrvRU9v+3uM5THJJUui7Rzo2vVGmwTBXUVoVpbQ5SCfv2btfg2GusrBKdvqhP8y00Ys1THLlhZU0HtWUpM3tTElaKgpc9UeZ3rzd01LDV2+UabuXKStKldZqc5jqTi/RxnaiJCmVF5UKFRomuW7djeV5ZlxXZjWMC8VJoaqStka5KkeqKim2pZLcyjiOkrTQW3dGKitHgWcmNV7aM1q70YvVG2bKS6s4LXWnl2ihrBQFruKkUN60Ms6943bw2KZ5Odnf3VHjef2wsxcSAAAgAElEQVR9kWdKgRMo9F0FO7dlBL47uUUDAAAcr0f9zA19V2Fwb/1uMzhRn9ndZqBGNA4mjchTtxk81e3vPZ6tyFer4Us62rHpNgO1W7481ygMXC21w8l79u5Xu+Xr4nJLUeDKdY06rUAr3YZakS/fM/I8o8V2qMV2qKWFQFHoynUddZqBnjnfnbS30A0UeOP+XGqFOteNtLIQKYpc+Z6jKPDUinxdWG5odSGS5xo1I0/PnG+pEXlyXWmx4+/cPuyoHXpaagfyPaMo9HTxXFOdaHwL8m6Ne610G+q2AvmuUSP0dK4bqdMK1Aw9NSJPvmf2HbeDx3axHU72t9Xwaz/3noRTVVVVdxGStL7en0m7q6udmbWN2TrrfXfSnk85qrPeb/OMvptf9N18ot9Ojkf9zF0+19abN7Z4pvQQJ/mZ0mefXZr83PFM6dO3utqZupxQihOLvptP9Nv8ou/mF303n+i3+UXfzS/6rl6HhVJu3wUAAAAA1IZQCgAAAACoDaEUAAAAAFAbQikAAAAAoDaEUgAAAABAbQilAAAAAIDaeHUXAACnzUmbrw4AcDx2r+9ZYRWnxYm9zmdZqY1erJVuQ8Y4T33O0oPzpB7lc3HaOgeX7Z3DtBF6Gsb5fXOAWlupP8wmyyXpxsZAZWnlGaO3NhMVSaYgcNWMfLmOoxt3h1pqhYoiT6Hvahjn2uwn8l0jzzOKgntzgBrjKPRdpXkpa6vJ/m4PU7mOo8JaBd64XmMcecaoH2dqhJ4WWqF6o0xVVcnI0eYw1VIrlOM66jYD5aWVJPmu0SDJ1Y58lVW175jsHtt25E/Wn5c5Sh+EUAoAx8jaStdv9pXlpQLf1dW1ztx/UAAA7l3fk7TQd2705DtSFHgn7jqfZaX+8f+5rmGcqxG6OtdtKC+sGpGnF19Ym3kwLQqrl755U3FSqBF5+j/edUFvbAwf+Lk47bNT0r5ll1daeulbN3VjfShjpOVupFdv9JQXVhfPtfTfXrwqzzP6zmtb+urX31SaW51fjGSt1beub2p7lMlWlXzXleM4Or/Y0DMrTd3aTJRkhayk97xjVV7g6ruvbeqtjaGKUlpo+Vpoh1pohzLG0XI3ku8Z5UWpu71EZWn15sZQvUGqUVrKdR25zji4RoGrtLCqqnFwXGz7KkqpN0o0GJUyjlTJ0TsuL6jZ8LXYCmWdSv3huO+SrNQ7Ly+qEfm6utaRtZVe+uZNDeNcwyTXQjOUcR1dudDRtYvdE3UePipCKQAcozQvleWlJCnLS6V5qUbIpRYA5t3u9T0vrYaZVSswMs7Ju85v9GIN41yStNlPVVXSQitUnBTqjTItd6OZbr83yhQnhSQpTgpt9OKHfi5O++zcfb37/41erMEwV1FaFZlVmo0UJ4WMcdQfZtroxVrqRNroxYrT8fvWN2NlRam0qFQUVkVZqfIdqarUj3Pd2oq1PczlGkdZYXV7OxmPbI4KZUWlvCg1SBzJGFWO1PA9Db18UvcoKTSKcw1HhZLCKstLmdKRcaSikvKiUmZLecZIlZSXVr7najSyitNcrmdUWWlzkCorrFzXyKmk/jCT64Tj8JkWcl2jNC8Vp4XipFBZVuoNMxlj1NoZMT5p5+Gjmt/KAeAECn1Xge9O/mU39N26SwIAHIPd67u1lVoNI9/RibzOr3QbajV8DeNcS51w30hptxnMfPvdZqBG5E1GSle6DSWZfeDn4mGfnXuXrXQbard89UaZAt9ouRspzgrlhVWnFWil25DnGa10G2qErtLcanVnpHSrnyrxjIxbyXeNHMdRp+HrwmJDqhwlWSHjujq/EMkLXG0PEw1GjozjqR15Wmj6WmiNR0pbTX8yUpoWnkLfaJQVKm0p67tTRkrNfSOl1hay8icjpUvtUM2Gr07kyzqVyqpSEBi15KsVepNj4rtGjciTjXN1W4E6kS/jOmo1/BN3Hj4qp6qqqu4iJGl9vT+TdldXOzNrG7NF380n+m1+nyml7+YXfTef6Lf5s3t9X1tb0M2b2yf2Os8zpYc/Uxo0Q54prdHqamfqckIpTiz6bj7Rb/OLvptf9N18ot/mF303v+i7eh0WSpkSBgAAAABQG0IpAAAAAKA2hFIAAAAAQG0IpQAAAACA2hBKAQAAAAC1IZQCAAAAAGrj1V3AWTCvcxYCAB4f1/6z5TT192nal8fxoP3f/V5pH39GxYPtz+J4H5wndFb21m5tNdmmMY7itJB0bw5Na6upy9K8lO8apXk5+Z4kpXkp1xnPseq7432I00KOI0WBJ88zch1HG71YnUagTiuQJA3jXHFaTNrM81KDJFe3GSgKPPWzUkk/1TDN1Qg9tUJf12/2lBallluhCmvVTzL14lyu6yhwXXXDQEFg1Ix8BZ6r7WGqO/2RXnurr7dfWtJiK9Qrt7e11m1qUObyHVfvfHZJnme00Yu13I5UVpWKwj50LlZjHPmumey355n7jmFRWOWl1WI7lDHOqfh5JZTOmLWVrt/sK8tLBb6rq2uduT5hAAAPx7X/bDlN/X2a9uVxPGj/934vsdJC+Ogh4GD7V8639drtwbEe76KweumbNxUnhRqRpxdfWJtJMN27L8ZxdHNzpCQtFIau1habunFnqErSlQsdXb3Q0fWbfb1+u79v2Wu3B0rSQuu9WEVh5TiOLq+2ZYyjJCn07Tc2FaeFBnGmykr9uJBUaaEV6tm1tm7cHirJSgWBq5/9LxflGEf/65V1rW8lGqW58rzUnV6qqpICz1Ez8tWIPG32M/meUeQbJVmuu71MaWa1k32Vlvv31XelVuir3XRlHKO72yMNst3vvinXkcqdf6fwXMkzRlcvtPXc2oLSolBWWD2z2tZrt3oqS8kYaaEVSqq0Pcxk7XjZYjvSUjdUf5gpyQqN0kLXLi7o6sXu5Bh+/+a2vv9WT+0o0KXVlp4511Jpq7n/eX2kULq9va2FhYVZ1XIqpXmpbOdffrK8VJqXk38BAgCcTlz7z5bT1N+naV8ex4P2f+/30qxUavTIx+Zg+71RduzHuzfKFCfjEck4KdQbZVruRk/U5jR792VrmKo/HAe9wTDXuhMrL6yk8chlb5RpmOT3LcvyUnlpNRjmkiTfM9oapAp9V2leqjfMVJaVBnGhrKhU7GzPmFwbm4m2h5lc48imld66O1TguxomheKs0HBU7KvRWqmwkq2kJM1VWk9lbjXIxnXZSirLg3s5VpRSWpRSKkmlsnz/98s9A+dlKTmqtNFPFYUjdVq++qNM65sjjZJCdnwIVElyKmm401eS5LqZXMfR9jCVJCVpqX6c7zuGw1GhJC0VuKW2eqmaoaduM5j7n9cj/bPJt771Lf3iL/6ifumXfkm3bt3SL/zCL+gb3/jGrGs7FULfVeC7kqTAdxXuvAYAnF5c+8+W09Tfp2lfHseD9n/v98Lg8Y7Nwfa7zeDYj3e3GagRjYNJIxoHllnYuy+LrXBy+2y75Wt1oSHfG9962mr46jYDtSL/vmWB78p3jdotX1HoyvOMFtuhWg1fzcBTtxUoCl21G56WWr6i0FMUuuo0fK0sRVpoBfJcoyj0dHG5pZWFhlqRp0bgqdX01Gl64214rqLAU7fpq9XwFIW+osBVq+lpoTmuyzhS4I7/O8hzpdBz1Q5ddSNfgb//++6ewUnXlVzjaKUT6sJSU65x1GkGWl1qqhl541uGG54WW6EW2oFajXvLOo1A3U6gxU6oKHQn+7r3GLaa42MQ+K4Wu6GW2qGk+f95daqqeuhN8b/+67+uP/zDP9Tv/d7v6Ytf/KK+9rWv6c/+7M/0+c9//tgKWV/vH1tbe62udmbW9lGd9WczHtdJ6Ds8OvptftF3x+tpXvvpu/o9Tn+f1H4767+3HOWZ0mcuLeruncGxtM8zpU/3mdKgFfBMaY1WVztTlx9pfDeOYz3//POTr3/6p39an/rUp46nsjPAGGduh9IBAI+Ha//Zcpr6+zTty+N40P7vfs99gl/+D7Y/i+PteWYmt+wetLd2Y5x922w1/PvWnbZs9/0Hw/Pu8pWgMVm20Anvq+FytD/kdFr3Auo0u/8YtLpn2X95+8qh609zbqmha1rUj//QvWXPPTv9EcdLK+1HanvX3v3eNe0YSo9+G/lJdKR/OllcXNTLL78sxxn/AP7DP/wDz5YCAAAAAJ7YkWL1Jz7xCX30ox/Vd7/7Xb3nPe/Rc889pz/90z+ddW0AAAAAgFPuSKH0ypUr+tznPqfRaCRrrdrtxxuGBgAAAABgrweG0g996EOTW3an+cxnPnPsBQEAAAAAzo4HhtLf/u3fliT93d/9naIo0vve9z55nqd//Md/VJqmT6VAAAAAAMDp9cBQ+hM/8ROSpE996lP6whe+MFn+Iz/yI/rlX/7l2VYGAAAAADj1jvRMaZqm+v73v6+3ve1tkqRvf/vbKopipoUBAAAc5qzPpYmj2zsXZl7a+86Zg/NnSnqkc+tg+9O2c9Tz9TjOa2sr9YeZ7mzHcl1HF8+15XnmsfYp9F1JeuDxO+w9Rzmm0+YpDX1333ayrNSNjYHyvJQxjqLAm8znGXj31pWk/jDT5iCR6xhdWG4qCFwlSaEbd4dqeq62s0zdN3vq92OtdCIFgavBINP1Oz2tdVtqNjw1I1+jQa7/9YN1XVppy/ONItdTpUqN0FfkuXpza6jQOHrjbl9L7Ybyymo4zLTcGb/O0lIXF1vqtAIlWaG8sFpoh2pGvja3E72+2deFTlMXV9sKfVcbvVitcDzVy/YwVSP0tNSJNIxzrW/HOteJ7uvD3flMB0muduSrrKq5vh4eKZT+/u//vj70oQ/pwoULqqpKd+7c0ac//elZ1wYAAHAfaytdv9lXlpcKfFdX1zpz+4sYZmv3XEnSQnf6ic51I0WBNzlnrK306o2eXr/dVyXp8mpbxjgqCnukc+tg+0vtUJuDdN92JB3pfD2O89raSt95bUv/499f0w9uDRUGnl54blH/9flVVdIj7VOWl/JcIzlSlpVTj99h77FVpTfXBw88pnuP3XovVlHYcY2eq5WF8XYuLTf13196VV//3oa2Rrlc4yj0jSLfG88X6xo1Qk/GSLa0eu3WQIOkkGuklYWGrl5o6eU3euoPEsWZVO3U6xopCoyagaM7/VJlJTmSuk1Pvim0Mbh3PHbnzwx8I89YlZWjNKtkH9IXRlKrYWQryZGjZuSqE7l6806svJQ81+iFKwsKQ09pZpVkhVRJw7RQM/L17GpDb67H6sWZPONodbGhSo6MkRbbkZa6oQajXI3QVZKVeuflRTUif26vh0cKpT/zMz+jr3zlK/rOd74jx3H0zne+U543/5O0AgCA+ZPmpbKdUZUsL5Xm5amYPB7Hb/dcyUurOCmUN62Mc++cSfNSwyRXXowjxtYgVei7Cn33SOfWwfZD371vO5KOdL4ex3md5qU2erG2R4WKspLJS61vxVrfjrWy0HikfZKkYZJPlk87foe9J83Lhx7TvcduMLy3nTQttdAKZJxSN+4OtbGdKCkqFYVV6YwDcFk5Kq1kjDMJlElWKs0LFaWVtY56o1yv30mUpIVKq30hsrRSmltZK5U7SbXabeNAntt9X1FalaVkq4cH0t33xYmVMZLjOIpTKS1KZfl4W4Wsbt5N1G75CjxXg1Ghwo5bTtJCb26M1I8LFUWltCpl+qkawfiYu24m13E0THK5TqhhnGuYFnJ3Rpzn8Xp4pIo/9rGPTV3+yU9+8liLAQAAeJjQdxXs/IIb7PyyC0yze65YW6kRefI9s++cCX1XrciX7xlVkhbb4b5RvYedWwfbb+6ErYPbOcr5ehzndei7Wuk2tND0tNV35PuuVhcbWl1oTEZKj7pPWV6qFfmTkdJpx++w9zRCT8M4f+Ax3Xvs2i1/30jp7nYuLTe1shDprfWBEs/sGSl17xspbYVGw5GnvByPlHabvp49F2kYZ8rzQkb7R0pDfzxSmub3RkqjwJVvCu3J4pORUs999JHSRnRvpLQRjkdK0/TeSOnacjQZKW03vclIaRR6embl3khpZDyd64STkdJOI1C3E8i4joLAqCVfrdCb6+vhkULp7h88kqSiKPTlL39Z165dm1lRAAAAhzHG0dW1Ds+U4qH2nivvcBfveybSGEfXLnV18VxT0qM/Uzqt/WnPXh7lfD2O89oYR++4sqiL55qP/UzpwTp2j8e04/eg9zyz0pJ0+DHdd+yeXTz0mdJf+unn9ePvvHB8z5S2mzxTegI5VVVVD19tv6qq9MEPflB/+7d/e2yFrK/3j62tvVZXOzNrG7NF380n+m1+0Xfzi76bT/Tb/KLv5hd9V6/V1c7U5Wbq0od45ZVXdPv27ScqCAAAAACAI92++653vUuOMx4KrqpKy8vL+t3f/d2ZFgYAAAAAOP2OFEpffvnl+5ZlWXbsxQAAAAAAzpYj3b77q7/6q/u+ttbqV37lV2ZSEAAAAADg7HjgSOmHP/xh/cu//Iuk8S28kzd5nn7+539+tpUBAAAAAE69B4bSz3zmM5KkP/qjP9If/MEfPJWCAAAAAABnxwND6T/90z/pve99r9797nfri1/84n3ff9/73jezwgAAwNlkbcUcpDhWRzmnHuW8K22lOC0e+Ry1ttIwzpXmpRbboTzPHNv5vredLCv1xsZAS61QUeSpEY7n9ty7jrR/7tCDdUxbd9ocrNNqKHOr1zcGeu58R1HkKU4LSeM5SOO0mOy/MY7SvJRTaTJXp+cZGeOoEXqyttLdXqI8L7XVT7Tej3Vpqa1LK21tjVJVeaXvb/TUjXwNi1yjUaF205dTGV1cbGorztT0XPXzXKFx5TjSs5nVxq2ebvSG8l2jOCs0HGY6v9BSI/A0KnI5cpSXVpLkONLG9khFUmppsalOM5BxHBnHkes42ohH2rjdV1kaXXtmSRcWm3p9ayDPdVSWlYzjaLERKgpc3R0m6se5ltuRrl1cUF5avXmrr41BrLVuU6m1Gma5BsNc73pmSUHkyXeNyqpStxloGOe600+0ttjUIMmVl6VWF5oqq0pFYZWXdnJezZsHhtKvf/3reu973zu5hfcgQikAADhO1la6frOvLC8V+K6urnUIpngiRzmnHuW8s7bSK29s6dbt/iOdo9ZW+t4b2/qPV9aVl5Uur7b1kz90QW9sDJ/4fN9bvy0r/d/feEu3t0aylfQT7zyva5cXdfVCR6/dHijLS3mukRypKKwC39WV8+3J9w5+vbtulpW60090rhspCrz7at2tYWuQ6Ev//JqqSmpGrn7hf7+i9V6ssqpkjKP1zVilrXRppaVnzrUUp4X+5eVbipNCWWH1zEpLq0tNPbPS0o27Q33r1Tv6/lub2h5YyZE8V7p4rq1G6Oh7b/ZVlFJ14Hg4klwjBZ6U5uNgWVVS6Bv5vpQkVml5/3F0dH9b0+xGPnvwG//fW3IdyVb32nEk+a5kHCktxl8HgaOr5zsKfUffem1b2ZRaQt/RO59ZUOk4urDYkO85euP2SEmey5aVwsCTVaVznUiXL3T0+u2+2lGgy+fb+ql3r81dMH1gKP2d3/kdSdInP/nJp1IMAAA429K8VJaPf0PL8lJpXqoRHmmyAGCqo5xTj3LepXmpNHv0czTNS20NUyU77+0PM2304mM53/fWf3NzqO1hprIcL9/op7oQ5+qNssk6wySXNB65zPJy3/cOfr27riTFSaG8aWWc+2vdreH25kijpFAUuBqMCl2/1VPgu8oLq1FSKE4Kua6jrV6qZugpy0v1R5mslUZJoa1hqk4z0MZ2rDubiQZpoTS3KiU5lZQX0tYgUZJ5KqcEUmm8rLBSlUtlNX5fJSkvrUqrqSFw931HcV8Y3aM80Mh4u/deV5KKotLtfizPMSoPaSwvKq1vp2pEnvpJrjwtNEhyqZIGaa6orBR4ru4OMoXBSElaKnBL9YeZeqNMy93oiHtzMhzprP/qV7+qP//zP9f29raq6t6R/vKXvzyzwgAAwNkT+q6CnV+UA9+d3DoIPK6jnFOPct6FvqtwJ0g8yjka+q4WW+PbOPOyUqcVaKXb+P/Zu7MYS/K70PPf2Lez5l57dfVKg1kHL8gIW+O5SFzJXPDI0mgkroSw8MOIJyRA+IEnWyBeWCReR2C4aNAgS2OEZozV92oGNYu43nC33d3uruqqysrK9Wyx/yNiHk7mqcyszKqTVXk6z8n6faRWdUae+Mcv/r+IyPM//4jzI8nKpz7e98e/0g5oBj3SXOHqBgt1h8CzaPg2ncFwsBm41oGZ0v2/O/zz3muzrMBzTSxTPzLWvRiW2j6+a1JVEHgG15cbbPRidEPD80zKnYqirGg1HNo1hzhV1H2bOFHgmsM+ckwWmh5ZUbK+E+JYOkn6YKa0VXPxHI2tXkJ1gplSyxjOlFJNaKYUxpopNU2NpbqHY2l0+inFEbFYpsZi06HQNOquhVWzyYvhTGlDs0YzpXM1m8W2T7Kbu3pg0/DtMfZiumjV/lHmMX7+53+e3/7t3+bFF19E0x5M01+6dOnUAtnY6J9aW/stLtYn1raYLMndbJK8zS7J3ew6b7l7Vp4pPW95m2an/Uzp3HyNu6sdeab0mBim+pnSK3PyTOkZWlysH7l8rJnSdrvNJz/5yVMNSAghhBDiKHtvSIU4LeMcUyc57ownPEZ1XaMe2NQPLTuN431/O65r8sLl1iNfAxz4/8O/O+61jxrwjNZx4Idqc6PlgWeN/r9uHtz/vXYvuw8PVnRdY2nOB+DS8sHf12rD2cArlxrHxnPtiGWLi3Xqts5zPNw/p+W4mK7RfGhZs+6M3a5tG7Sbw9ty9/b/vBjrDPipn/opvvSlL/GzP/uzOM6Djvvpn/7piQUmhBBCCCGEEOL8G2tQ+u1vfxuAN954Y7RM07RRHVMhhBBCCCGEEOJJjDUo/Yu/+ItJxyGEEEIIIYQQ4hk01qD0V37lVw78rGkaruty48YNPv/5z9NsPnx/tBBCCCGEEEII8ThjDUqff/55TNPkM5/5DABf/epXWVtbY3l5md/93d/lT//0TycapBBCCCGEEEKI82msQem3vvUt/vZv/3b08yuvvMJnPvMZ/vAP/5CvfOUrEwtOzIZn5av7hRDiMLn+CXHQtJ4Tk4rrUSVWTmvdk8ZelhX3N0K+f3eHVy63WVoIHlkGZ9wSNUctV6qkF2U0fBvT1Ec/u6ZBJ0pp+Q6JKrB1nfVezMW5ANs2Ru0oVbLeifAskzQvhrVVbRPXHv5sGBqWpnNzo8/ldsBaJ2IrignDjJVWjcA2udnpUeUVzbpDripsWyfLSmxbJ0kLPNtkIXDJypKaa7HZzfjX79+hZtvUfZu8KLiz1sVzbbzAghI810SrNMqq5N56n6VWjecuNNnqxdzZ7tPwHVRZYVoavUGGKhTRICPwHGoNm14vo91wAY227+A6BioreX+nz0ojwHWGOTdNHdc0WetG2LqOpkG74TLf9Oj2U25t9LmyUMO2DMJ0WEomUQW+bdKLMsIwY30QsxC4mJZOq+ZSD+ypOvfGNdagNM9z3n77bV588UUA3n77bcqyJEkS8jyfaIBiupVlxc21/qjY8vWV+kyeCEIIcVJy/RPioGk9JyYV1/52TUMHbThIG2cb46570tjLsuKbb2/w5//390jzEs82+M//4RU+9OLCkYPdd+50+dYPNsiLisuLNT76Q8vc2Qwf2t5RcZRlxetvrBEnCs81+emXlvjXt9bphSlv3+nS8m06UcaNlTr/fnOHwLUIPJOP/8iFYfHPCv795hb3Ngds91PKoiRVJa5l4LsWaBqOWbG6laJR0gszcgXFE+TKNsC1DYqiIMqgeoI2LA1UdfJ1DQ1sHZLiwbqeBYZpYBmQq4pclRQFuI7JXMPmxYt13rjVI0xzDA2W52q4jk5ZwguXmtzbDgmjlLfv9lBFNSyd0/K5uBDwiR+/xItXWlNx7p3EWIPSL3zhC3zuc59jfn6esizp9Xr8wR/8AX/yJ3/CL/7iL046RjHF0rwgy4eXhywvSPNCassJIZ4Jcv0T4qBpPScmFdf+dsNkOEnjWMZY2xh33ZPGnuYFb9/pkKty9PP373Z46Xr7ofXSvKATpiTZsP1+mLHZi4/c3lFxxKkiThQAcaJY3Q6JE0UUF4SxwjQ0wlixthMRJQrbMuhHGavbERfnAzZ7MZ1+SqpKwkRRVBWUFUVVkRcVtmUQJRCnOaalkSkox03OIVkBuipR+ZMNSAHyJ1yxqCA+NJJOc7CqkrKoyBVU1XDAmxclYVrw7v2IOMupKkhVQSdMaGkuWV6w3U/ohhm9gUIVFUUBRVkxSDL6sc1GN+bqSn0qzr2TGCvaj3zkI/zDP/wDb731Frqu8/zzz2NZFj/5kz+Jps3WKFycLscysHcvorZljG5BEUKI806uf0IcNK3nxKTi2t9u4FoHZjsft41x1z1p7I5l8OLlFq9/d400L3Esg5cvtY5cz7EMWoGDaxvkRUU9sFloeCRZ+dD2jorDMnQ81xzNlF6cC7i7GZIXBYFnUnMsVFGx0vbZ7KZYhk7gmVyc8wFYaHi06hFxmhO45rEzpXFqoVFimzzVTKlt6hhaQTEFM6WOBYapYxmg68OZUgqGfeQY3Fj2eSNRhGm+mycX19axTZ25ukuaF5haxVZPA4YzpTXXpu5ZLDa9qTn3TkKrquqxfXvz5k2+/OUvE0URVVVRliV37tzhL//yL08tkI2N/qm1td/iYn1ibYuhST2nIbmbTZK32SW5O7lpeX5OcjebzmPepuWcOOy049rLnTxTOnvPlHq+K8+UnqHFxfqRy8calP7yL/8yn/jEJ3jttdf4pV/6Jb72ta/x/PPP83u/93unFqAMSsVhkrvZJHmbXZK72SW5m02St9kluZtdkruzddygdOwvOvqN3/gNlFK8+uqrfPaznx2VhxFCCCGEEEIIIZ6UPs6LPM8jyzKuX7/Od7/7XVzXnXRcQgghhNP4GwsAACAASURBVBBCCCGeAWMNSj/96U/z+c9/nk984hN8+ctf5td+7ddYWVmZdGxCCCGEEEIIIc65sW7f/YVf+AXKsuSv/uqv+PCHP8x3vvMdPv7xj086NiGEEEIIIYQQ59xYg9LPfe5zvPzyy1y8eJELFy5w4cKFScclhBBCCCGEEOIZMHZV1S9+8YuTjEMIIYQQQgghxDNorEHppz71Kf7mb/6Gj370oxjGg2KsFy9enFhgQsyaaa3NdhrO875NkvSbEEKcnsddU09SX/Ok9upv+rZJlCma7WD0uyRRrG6HXJwLcN0Hb62P2m5ZVsSpoiwrdF3Dc4avT/MCy9DJi3L0+qeJu9NJ+NbNTV6+2DqyTun+OAB0XcMydAZJTsMf1rl8XF/HqQLAc8zH1lYdty/KsqIzSLEMfdQvWgXvrXdZXe9hmQbthkfTc8jLglSVWMaw1ud8w2WrG7PWCwlci6wo6YcpKleUCuoNl0TlhIOETidjoeljewZhmOJZFl7dIcsUDd+lMirW1nqYmsHSUo2yAtcx6CcZO5shnmsz1/boxhl5rGg0HdbW+sOasxpoaNiORa1mEUaK+bpLnCmWah49lZOrElPXaTg2WaG434m4MtfAtQy+v7bDSxfarCwEbPZiDE0b9WGmCmxzOBZTZUnDd4iinDfvbnOpVePicm3q65QeZ6xBaRRFfPGLX6Tdbo+WaZrG17/+9YkFJsQsKcuKm2t9srzAtgyur9Rn8oJwlPO8b5Mk/SaEEKfncdfU435/GtdipUpef2ONfpRxbyvi0nzA7a2IH77SQqmSv/z6W4RxTuBZ/K//40u4rnnkdgHeXe1x636PrV7CfMPjylINXdfIsoKtfsJ8w8W1Ta4u1Xh/ffBEcXc6CV/6L/9GlCgsU+c//4dX+NCLCwcGg3txbHZjQKNVdxhEOYFr4rkmF+cCirI6tq/fXe1xe71PBVxdrnPjQuNA+/v3/ah9OaovLi0ErO1ErG4M6EUpVBq5Krh5b4edsHpoPzWgAnQNHBNMHaIMiodfelxPjfvCD4wOaDrovMdzF+pUQC/KsXTIi4qyqtAYjsNsy8A2Ku5uxWQKTF3jR56b53/66Su8eKU1c+85xhqUvvbaa7z++utSCkaIY6R5QZYXAGR5QZoXo08/Z9153rdJkn4TQojT87hr6nG/P41rcS/KiBNFlpWEcU6cF4SxohdldAYpYZwDEMY5q9shNy42j9wuQJjkJGlBmhUkqaIzSIeza0CcKHK/RNcKelH2xHF/59YWSVZQAZkq+f7dDi9dbx+Yld2LI0oUGho6GmGS41oGgzBnx05p+PaxfR0mw9m+vf3e/5rD+37UvhzVF5vdmE4vJc0LBrGiKCFXJUl69Chzb2lZQV5AlkMxVg9NrxLQK1AVbHQTfMcizQtSDapi2N9VVYGmUwL9uCDLh32hyoqNbsxGN+bqSn3m3nOMVRLm0qVLdLvdSccixMxyLAN794+KbRmjPzDnwXnet0mSfhNCiNPzuGvqcb8/jWtxw7fxXBPb1gk8C88yCDyThm9zcS4g8CwAAs/i4lxw7HYdyyBwLVzHwLENXMekVXMIPAvL0PFcE8vUsS2Dhm8/cdwfujaPaxtogG3qvHypdWD9/XH4uzOjjbpNPbAxDI1aYNGuOcdue299y9QxzWGfHG5/f+xH7ctRfbHQ9Gg1HBzLoOaZNH2Lpm/iOkfP+O0t1TWwDPAdMGZrcvAhOqBpw1nfxaaL7w77qmYbuLaJbQ3/9WwD1zKYr1nY1t46GotNj8WmN5PvObSqqh47yf2rv/qrfPvb3+bFF1/EsqzR8j//8z8/tUA2Nvqn1tZ+i4v1ibUtJmvWcneenx88yb7NWt4madaOCcnd7JLczSbJ28lM0zOlz12bp7sTAvJMqTxTKs+UnsTiYv3I5WMNSv/lX/7lyOUf/vCHny6qfWRQKg6T3M0mydvsktzNLsndbJK8zS7J3eyS3J2t4walY91sfJqDTyGEEEIIIYQQYs9Yz5QKIYQQQgghhBCTIINSIYQQQgghhBBnRgalQgghhBBCCCHOjAxKhRBCCCGEEEKcmYlVVS2Kgi984Qu89957GIbBl770Ja5evTqpzQkhhBBCCCGEmEETG5S+9tprAPz1X/81//zP/8yXvvQl/uzP/mxSm5uoaak1OC1xCHFW5Bx4tEn0z/42Aen/KSDnwWw7fE497jWnkePH1ZQUj3ZUPc1id5mhaXQGKZoGdc8myhQN38Y0x78ZcS/fh+uUniQ2pUryoqTh2xRVRZGX3Frvs9T08H0LxzJGbQNjXdePu/4Dxx5Pe+tURcVaJ+LiXIBtG/TDjH6cMVdzSVRxoI/KsiKMc+JUoVWwE6a0A4eiqsjzgs1+zE6YEqUpg17OxcU6lQEqLbnWrvPt2/fRK4Mbl1pspQnhICVLcvr9hDAtWWzarG4PUFlF4FrUmj4122Zu3mOrM2B7KyIrYC4w2Y4UizWbQarQdJPFBZ+iKKAC0GjUXDIK+jsJlqOzvR0TOBZJGvGdd3v4vsuP3mhybydmPvAI6jZxErO2nfOTLyxTq1nc3uyhozPX8Lg0XyMKc/717TUuLNQo9IokSgkcl5cvt8nKEk0DWze41424NBeg6RpFUbLdi7m3E3KhVaMWWMw1vKmvUfooExuUfupTn+ITn/gEAKurqywsLExqUxNVlhU31/pkeYFtGVxfqZ9JsqclDiHOipwDjzaJ/tnfpmnooA0LyEv/nx05D2bb4fwtHFGv77RzXJYV7672uL3epwKuLte5caEhx82Y9vrv1v0eW72E+YbHlaUa27Fiba3H925vs95JqKoSxzS5vFgj8C0+9urKWAPTvXwnqWKrnzDfcHFtc6y878X23lqX9+71CBwL1zG5uhzw//zLbeK0AA1+5odXqPk2C00X2zRG1/JHXdePu/6bhk5ZVdzdGDx0PO2t0wtT/us376KhEXgmH355iX/+3v3dgT28dKVNPRj2ka5rvHOnyzffWef+dsTaTowOFGVFzTW4vxPRT6pDe752TI/ceWx/QwoMxnjdE4oTbm0luz/0DvzqG+92D/ysMRyI5aMl9w/83tKhVXNQpSLNhx+IaMBy22OrE9KJKqrddhYaDs9davDJH7/Mi1daM3l+T2xQCmCaJr/1W7/F1772Nf74j//4ka9tt31M8/hPDZ/GcUVaxxElOX4/w9/9ud708F3rdAKbwTg+aE+TO3F2JpG3Z/UcGNdp9c/+3O1vM8mGn4rXauZTtS+ezqPyLNfL6Xc4f2mmHsrbaV/roiTH3IpwXBsA0zbl/D2Bvf4zTBPdMDAsg0rXGEQZlmsS5xVoGkUBcV5gWAaGaWJ5Fgstf6z2/X6Gbhl0YoUfOLhj5mgvtgqdCo1K08jLiu1BRpqXw4FkUdKNcvzAwQ+c0bq1mvnI6/px1/8kU6SZOvJ42ltna5CSFxWBZ5KqknvdmBINNI1EFVQ6oz7yXYvqXo8CnVSVpKrE1DXyooQE8uLwgPT8qNg/IH1YXkKSF6BV5EWFrWmosqQXKbJyOCDdaydRJUlWklbVzJ7fEx2UAvz+7/8+v/mbv8lnP/tZ/u7v/g7fP/oE3dmJJrL9xcU6Gxv9J16/LCuiMB19YtnvxoT95PErnrJpieOD9LS5E2djUnl7Fs+BkziN/jmcu/1t7n1S3oky6f8zdFye5Xo5Gw7nz7HNh/J22te6sqxQqSJNMipAZY6cvyew13+FUpRFQZEXaGVFzbcZ9BI8S6NfVRg6OKZBkRcUliKPczbyx5+Te/lOdrcRhSllXoyVo73YNEo0KrSqwtI15mo2jqUTp8Nrd9O30KvhdvZmSjtR9sjr+nHXf9PQ0arqyONpbx1LA8vQyLOCwDO50PR4/14PqgrX1NFKKNSwj/qpQisqDEocU8cxdXRAR6fmGoSxRpKfz4HpwzOlB1k6uJaBKissAzQNLEOn4ZvkmUbEg5lS19RxbR1H06b+/D7uA1StqqqJZPorX/kK9+/f59d//dcZDAZ8+tOf5u///u9xHOfI10/qj+lp/KGelud3piWOD4q8yZpNk8zbs3YOnNTT9s9RuZNnSqfPUXmW6+Xs2J+/5eXGkXmTZ0qny1HPlC4s1rm72pFnSuWZUnmm9IQ+8EFpFEX8zu/8Dpubmyil+NznPsenPvWpY18/zYNScTYkd7NJ8ja7JHezS3I3myRvs0tyN7skd2fruEHpxG7f9X2fP/qjP5pU80IIIYQQQgghzoHx7y0QQgghhBBCCCFOmQxKhRBCCCGEEEKcGRmUCiGEEEIIIYQ4MzIoFUIIIYQQQghxZmRQKoQQQgghhBDizEzs23eFOKnTrMs2jfUspzEm8WQkl48nfSSeBeepju9R5+y0ncfTFs+eceLKsoLNXsxczaXSeOi149SSPe41h7dflhVbOzG3Nwe8dKlFrWYfqIUaxjn9OGOp5WPbxti5H2fZ3s+GpjFIcmquRZoXpPmwNmlelKN6r441rF96byvEsQzmGu6odqlSJf04w7NMdvoJ9/oDOtsxuqERxzmB57A0H1ABG90IlSp81+HKfI2kKqnbFmuDkNt3ewSeSRRHvPXeNnlp8CMvLGIZGrquk6uM99f6WIaB79o4pkZeabhmwZu3u9Qdm+W5Ou2WS2VU3F3dYaeX0aoFrCy5vHt3wMsX2jjzFrdvdanbJqpUbA8ymq6NaZs4tkamdJ6/3MKydLaihN5OQrvlkamSvFJsb8QszdW4uFAjVjmb3YSlhk/ds7FtnZprE6Y5eVrwfqePYxjUXJtm4OAZBt+5vcXLF1r4NZtWzTlRndxpIYNSMRXKsuLmWp8sL7Atg+srR9cwetK2zvqP1zTGJJ6M5PLxpI/Es2D/cW4aOmjDN9KzeMwf9zd4ms7jab2ujBNXlhV89Z9uMogy8qLiw68sE3jW6LVlWfHuao/b630q4OpynRsXGg8NWo96DRzM09WlGm+8t83/8d/eIcsKWjWH/+0/fYjNQUqSKu51Qt653SVXJRfmA/7jR66xuh09NvfjLLu6VOP99QFRnPPW3Q6ubRAlOUpBUZU4lkHDt+mEKe26S6VVfPP766xtx2jAypzP5cUaFRWrWxFZpuiGGTu9DFWNnxPHhN2x+xEK7vzb2jG/C49YlvLmWvrw4o0uvNcF4I07R613hG/ce8wLNh5a4lo6gWtgmTplWbDRVex1hWWAa0BSDH/WNZ0Pv7LI81fa/MwPr8zcwFQGpWIqpHlBlg/Pqmz3E7XTbMtzzvZQn8aYxJORXD6e9JF4Fuw/zsMkB4azX7N4zB/3N3iazuNpva6ME9dmLyaMc4qyIoxzumGKZeqj16Z5QZjk5KoEIIzzh9o57jV72937txdl3FrvkWYPjs037+zQDBzyomR7JyVOFLqu0Q8zVrdDsrw8EP/hNsdd1osysrwgyhRhnKOjsdNPQdMwdY04UZQVpFlBGOX0ooxBqCiKiqqq6AwyfDeBCqJYoVRJkuYnGpACZMcOSGdPpkqMbPhvRcn+rlAFJCWoCnQNSq1kbSdmaS6gF2XMNdwzi/tJzNYQWpxbjmVg797+ZFvG6Faos27rtExjTOLJSC4fT/pIPAv2H+eBaxF4FjCbx/xR5+y0ncfTFs+eceJaaHgEnoWhawSeRTNwDrzWsQwC18IydUxTJ/Csh9o57jWHt9/wba4tNXBsA43hsflDl9vYloFl6My1HTx3eOtvPbC5OBeMlftxljV8G9sy8G2TwLOwbZ123aHuDpe36g6twMaxDQLfYmXRpxaYGIaGZeq0ajZzdZd2w8H3TFxHx3UszBNOiNtn/1nFqbFNHdc2CFwTfzene0wDXGv4r66DqeustD3qgU3Dt88s5ielVVV1ws8fJmNjoz+RdhcX6xNrW5yuw88lPE3upvG5k2mMaRKehXPuvObyNHN3XvtoWj0L5900etpnSqcpb/JM6cnsz508UyrPlMozpeNbXDz6ET0ZlIqpJbmbTZK32SW5m12Su9kkeZtdkrvZJbk7W8cNSqd3GC2EEEIIIYQQ4tyTQakQQgghhBBCiDMjg1IhhBBCCCGEEGdGBqVCCCGEEEIIIc6MDEqFEEIIIYQQQpwZGZQKIYQQQgghhDgz56i8rBBnYxZqu4nJmeZcP21s07xvQsyKpzmPpvUcVKqkF2U0fHuq6iE+aX8dV2vzqHaUKukMUgzHIoxzPGf4VjpOFVlW0I8zmoFDPbAPrHu4zw7XHN3fRpjmLDQ8iqoa1fTcq0s6zj7utZ0kip0wZaXlk6jhOp5jkubFaLv727QMnbwoH9kP4/SVUiXr2xFhmnNlt/zH++t9PMvAdU2agUNRVViGTpwqwjhH08A2DdK8QNPA1HW6YTpcrht8a3WDe3c6VKokzgragYvjmpimTlEUhEnF0kKAaeroFqzfD7ENWN/s8v5mgqNDXoDrWLx4rc1ONyTNK8gTbq6X1Gsac4GFQqNZ86i7OmvbIXGSYzkOK02X7X7ETj/FNHQsy8DSNcoKqrIgLTRcvWAjqghMwLR5bsFnoAooSuK8wHdsFtsBmlGyvh0TOA7LizVw4c77PeYCB93WqAqN+TmfsgLdhE43wTFNfN/C0DW2+wmlKnAti8tzNTIqLrQDXNc8cKzMGhmUCvEUyrLi5lqfLC+wLYPrK8OL7+Fls3hxEI93VP6nJddPG9s075sQs+JpzqNpPQeVKnn9jTXiROG5Jh97dWUqBqZP2l+H17u6VOP99cGR7ShV8o/fvcfd9QEFcKHtc3V5+Hf/vbUu//6DLVRR0a67/NyPXeTFKy10XXuozz7yyjLvrw+4vd6nAi4v1kZtfPfdbSxTx7NNrizX6YYp8w2Payt1ri/Xj41t//68u9rj7dUd/u17G2hUUGlcXa5j2zpLLZ+yqtC04bK9NpNUsdVPmG+4uLZ5ZD8Aj+2rsqx47dt3+Kdv30eVJRfaLpqmcWczJMlKri3XWJkLePlKi81BwtpmyO31kKoq0bXdQV5VkmQFSVaQK8XOoDgic+ERyzYef6CEOe9trz+0uNeruNvLdn9KD/025s3VeH8v7/53tGEUGTe3skO/SYDeyWN+BB0IfJOaY/Hqc3OszAdcW6lz40JjKq4XJyGDUiGeQpoXZLufOGZ5Mfr08fCyvU9BxflyVP6nJddPG9s075sQs+JpzqNpPQd7UUacDGf44kTRizLmGu4ZR/Xk/XV4vV6UHdtOL8oYhPlwnaKiH+d0BsMBTK+fEyXDWb4oydnoxlxdqeM55kN9ttmLCZOcXA0HNvvbiNOCqoIsz/C7MUVZkezOJj4qtv37EyY5O92MJC0wdY1UKbpRRlCabHQTAsfEMvUDbeZFSZwocr9E147uh73/f1RfxaliYysZvX69mwEVuarIVUknzHAck26U0e1ldAcZaVZQFCXF7qxwVVaEeU5VVCh11IBU7CmBLCuINdjqJrRqDmGcT8314iTO/qMtIWaYYxnYlgGAbRk4u/8dXibOp2nO9dPGNs37JsSseJrzaFrPwYZv47nDN7uea9Lw7TOOaOhJ++vweg3fPradhm9TCywcy8B3TeqeRavm0Ko5NOoWvmtgmTq+a7HY9EbrHu6zhYZH4FpYpo5p6gfa8BwD29JpBjbzTQ/HNnAdk8CzHhnb/v0JXIt208Z1DCxLo+ZaNH0b3zVZbLq4joFp6gfatAwdzx0OVo/rh3H6quHbLM67OJaBYWgsNW2WWy6WqWGZOq3Aph04NH2bZsOmWbNxbAPPNWgEFr5r4nsGddcaxjUlx/200gHbNvBsi/mmOzpWpuV6cRJaVVXVWQcBsLHRn0i7i4v1ibUtJmtWcifPlB40K3k7LdOc65PGdjh307xv4qBn7bybJY86jx6Xt2k9B+WZ0pTl5Qa9TiTPlMozpfJM6Qkt7h4Th8mgVEwtyd1skrzNLsnd7JLczSbJ2+yS3M0uyd3ZOm5QOj0fbwkhhBBCCCGEeObIoFQIIYQQQgghxJmRQakQQgghhBBCiDMjg1IhhBBCCCGEEGdGBqVCCCGEEEIIIc6MDEqFEEIIIYQQQpwZ86wDEGdjWmufTdqzut9PS/pt6Gn6QfpQCDHLJnkNO622j6u3edJtHVVT9FHrZFnBeicicKxhnczdWpFlWbHZiUlzRavmUg9syrKiM0hxLAPHMuhFGZahY5o6jmUQxjn9OGOh4VFpHNgX4EBt0706oXGqUKokL0paNWdUOzZJFKvbIRfnhjUsT1LjdP82DvftcXVOD9cq7QzS0b4ZmsbadoimwWJzWDfVNQ06UUrNtljrRDiGDjo0fGfU7/0w4/Z2j9t3O0T9kNvbMU3XoRbYFGWFaxkkqsA2TSqV8v5WhKvrqEIR5dCuudRs+O7tBFXCUtvkpWtNNnYSlho+SRrxr2+FJAxn6pbrOjeW6yR5wp2tFM8E0zRp+A62pbG+PcByfF64FPDN729wdwAWsFIH3bJY9DXeupfhegYfe3UJw7Pobfe5eW+HXgwNz6LhG/TikiuLTZYu1Lhz+y7/8lbGUkvn6mId09BxbZMg8HjlUpv1JMHSdUwN3ry5gaMbaLqGbZsENYc4zSiyiivLTZ672KRZd2byfYbUKX0GlWXFzbU+WV5gWwbXV+pTefCedu5mZb+nzUn77byec09z/MzKsXdec/cskNzNplnJ2ySvYafV9l47SarY6ifMN1xc2zzQ3jjbUqrk9TfWiBOF55p85JVl7myGD62zl7ssK/i/Xn+Pe5shaV5yaSFgse1zaSHgztaAb3xvgzDNuTAX8PEfucBGP+HeZoimVRiaTqoKolRxfaWBbmj84E6XLC/wHYv/4ZUlelHGfMPFNg3KquLuxoAKuLpc5/pynZtrfd5b6/LevR411+byUo2f+eEVlCr5y6+/RRjnBJ7F//LJF1nrxI/c97KseHe1x+31/mgbNy40AEZ9u9GLUapE07RRDO+vDw60W5YV//jde9xdHzBIci4v1biz1md1K6KsKlo1hxcvNfnBvR4N1+TdtT5VVTFIFM3ApuHbXJz32RkkvPP+Nt3kiQ6tc0MDxhms6Rq8cKnBf/r4DV662p7K9xlwfJ1SmSl9BqV5Qbb7KVeWF6R5geec/0PhWd3vpyX9NvQ0/SB9KISYZZO8hp1W23vt5EVJnChyv0TXDrY3zrZ6UUacDGcK40Sx2Ysfuc5mL6YXZuRFSZQoOmFK3bfZ7Masb8UkWYFSFf045/ZmnyQtUUVJkilUWWHpOkla0BlkKDWMvaKiH2VsdWNUUZH7JbkqSfOCXJUAhHFOL8oIk5wwUiRpgW0U9MOMXpTRGaSEcT567a31PpqmPXLf07wgTPID20j37XtelAzCYZuWqY9iONw/caoYhMN1k7Rgcydhq5+Rq5KirOiGGRvdmDBWlGVJlCrQhh8IRIlCQ2PTjOnFBWl+4kPh3Bl39rCsYLuXsrodce1CY+beZ8gzpc8gxzKwd28DsXdvH3kWPKv7/bSk34aeph+kD4UQs2yS17DTanuvHcvQ8VwTy9Qfam+cbTV8G88dvpn3XJOFhvfIdRYaHo3AxjJ0fNekFTi4jslC02Np3sO1DUxTo+5ZXFmo02o4mIZO3bdYqLu4joHrGLRqNkvzHp5r7v7eZr7pjfYlcC1aNQfLHN4OG3gWDd8mcC0C38R1hvtW35tpnAsIPAuAwLO4tlR/7L47lkHgWge2sXeb8V7f1gIL1zEOxHC43YZvUwuG67qOwULbZb5u7+ZEpxnYLDY9As+k5dv4jolnDdv0XXPYN02Pds3EsZ7ocDhXxp3v1DWYazhcnPNn8n2G3L77jJqF59smkbtZ2O9pdJJ+O8/n3Hl/pvQ85+68k9zNplnKmzxTenCd/bmTZ0rlmVJ5pnR8x92+K4NSMbUkd7NJ8ja7JHezS3I3myRvs0tyN7skd2fruEGp3L4rhBBCCCGEEOLMyKBUCCGEEEIIIcSZkUGpEEIIIYQQQogzI4NSIYQQQgghhBBnRgalQgghhBBCCCHOjAxKhRBCCCGEEEKcGfOsAxCzbX9tKmDq6zA+jQ+yzuQs1LQUx5P8iVnzQR+zco5Mn8flZFI5m8Zj4biYDr/n2VsWp4qyrNB1DcvQGST5qMbpXg3TumfjOSa9KMOxDALPOlDTc7MTs9WJ2OjHLNQ86jWbC/M1TFM/NpYwzolTheeYo5qnjmXgOSb9MGOjG7PY9A7Urdxfh9SxDNK8GMV+VN3RsqwO1Gw93Bd7dVQP14bdq5G61PBIVIGhaaR5gWFo1D2bQZLjWAaWobM9SKjZFuu9mMAy6cYZnmWQqYJunNIZZGhmxZ3bO2x3IvqDlFwzWKpZ5BV4tkGUKpr1gLmmxd31LjudhH4ftooHebUZ1hQNgYYFizW4vQPZ7u8NoG5C4OtcbHuoqqQoKoqiIM4KirzkTh8q4HIboj5sDruSFQ96CTQCGAxgADRMMHWwTTBNjZZvYTo2vq54ez1hselyfbnFdneDf3q3YKUJV+YCtrshuWbz3JzLZpSiodEIfC6u1MlVTpYUaJqG61nolk6pSpqBRzNwWGx5LLb9A7maFTIoFU+sLCturvXJ8gLT0EEbXlhty+D6Sn1q/richv37Oun9+yC3JU6f5E/Mmg/6mJVzZPo8LieTytk0HgvHxXR4+cJinbKseHe1x637PbZ6Cc2awyDKCVyTwLP4qRcW+ft/vcXaVoRhaMzVHfKiwjJ1fuKFRZ6/1KQsK/7bd+7yT9+6y837EUUFhg4X5nx+5MYcP/b8IhU8FMs7d7p88511tnoZrZqFqetkRYlpaLTqDt/5wRb9KKfp2/zCR6/x0tU2AO+u9ri93qeoKkxDJ1cFO/2U+YbHtZU615frvL8+IMsLdE1jbSciSRWea/KxV1cwTX3UF0mq2OontGsOO4OU+YaLa5ustDz+y2tv0w9TwkRxeSFgoxtTVhqBa2CbBr5roesVYVygaxXvX4AutgAAIABJREFUrQ2wDNjqpViGTpIrlKrIiuMyVXC3c/iX6SNzm/FgANrLobdzuEXoKOj0Su72wke29f6hddfi4b/R4MGyntq3YSru9PZHABthwhuraw/a6MJaNxyt9P7Wg9dCAm9tPzIm29RYafv83I9f5Gd/7NLMDUxnK1oxVdK8IMuHF4QwyQnjHIAsL0jzY68iM2n/vk56/z7IbYnTJ/kTs+aDPmblHJk+j8vJpHI2jcfCcTE9tDxTpHlBmOQkaUGaFfT6Gf0woygq4kRxa71PL8woyopBnLPRicnzkjQr2BmkpHlBL8rY3knpxQVFNYyhLIfvqzY6MRvd+MhYOrsDvlwVdAcZm/2EPC+JEsXd9ZAoVhRFRZQqVrcj0t31wyQnVyVpWtDppYSRIs0KklQRxjm9KBvtZydM6YfDgVGcKHpRdqAv8qIkThRRqogTRa5Ksrzg1nqfMM5RZUUYK7b6KWFcEKeKMCnYHmREmaIzyOiEGYMoJ0oVYVoMY8sL8vxRA1JxFFVUhKni3k48ytUskZlS8cQcy8C2DLK8IHCtAzOl+29tOQ/27+uk9++D3JY4fZI/MWs+6GNWzpHp87icTCpn03gsHBfTQ8vt4S2zgWvhOgZOatCo2eiGhmFoeK7JtaU6b9zaJkoUNc86MFParjmj21fn2g4Nz2C7B0UFug6Baw1vxWx6o5nS/bG0AofANUmykuYRM6Xb/YQiyvEdk4tz/mjdwLWwTB3d0Ah8i1wVZEWB6wxndxu+TWcwHJi2Aodkd8DquSYN3z7QF2VZ4bkmvmMOb+U1dWzLYKXlEXgWZVgSeCbzdYeyLA/OlNomumtgGcOZUt8xsQxIUh3L0Kmo0B45UyoOMw2NwDG50PZGuZolWlVV1VkHAbCx0Z9Iu4uL9Ym1LSb7TOm05U6eKR3PtOXtLMxq/iR3s+tpcyfPlJ6NaTrn5JnSB8Z5pnR5ucHGRl+eKZVnSgF5pvQkFhfrRy6XQamYWpK72SR5m12Su9kluZtNkrfZJbmbXZK7s3XcoHR6h9FCCCGEEEIIIc49GZQKIYQQQgghhDgzMigVQgghhBBCCHFmZFAqhBBCCCGEEOLMyKBUCCGEEEIIIcSZkTqlQkypafya/JN4kvhnfZ/FB+fwsSLHjjiPpu24PsuyMEqVR5YmOattFbslWZQqR+0pVbLVT2gHDpqh0ao5AHQG6YFSLbfX+1QVXFmosR2lBJZJP81pBw6VBpahj8qnNAOHoqpGy2BYymWvDEsY5+z0EyxDx7YNmoFDXpQoVRKniqqqME2ddt1F17WHysHs///9pWEeVQomThVKleRFSavmYJo6SpV0BimWoaPrGlGUc3trwJWFGr5njdbTNPBda1QSpywruv2UuxsDsqKgKCt8y8QwNJbnAixD580729xa2+Cbb27gGAaerdGLFb5tU/cN7t6PKQydn7jRRgFRmLDeiel1SzYP1Tk1GM7IacCNCzZVlvH9rQe/t3f/U4Bng2tAnEPNgyiGnhquWwDl7jp1Cy41NW5tVsS7yy40da62XaJS48qcz/1un41ORmAZlKVidavC8TQ+dK1JN465t5USOBbzTZckzUjzilbDp+mbDJICx9QpNY2aZxLFGZmCpfmAVKW8d7uDrps8t9Lm1RuLPHe5iW2ffc3fk5KSMGJqPcu5K8uKm2v9UZHu6yv1qXhDMo7FxTr37/dOHP8s7/N5MSvn3OFj5epSjffXB8/0sTMruRMHPSpv03ZNnFQ847SrVMnrb6wRJwrPNfnYqytPNDA9rW2VZcVOrPjvb97j3dUuRQFVVbDVy8hVQVVpvHylxaXlGloFa9sRug6tusM3vn+fOxsxZVFg2ybtmsXOICdwLQBW2sO6nlU1HIAstXxevtJiO0xRqqQCbNNgru5wvxvx1q0d7m2GFKXGhXmP5bZPPbC5udalO8hJ84LFpscr19pcnA9Y3QopqgrT0CmKYXumoZOrgp1+ynzD49JCwHonJkkP9kFZVry72uO9tS7v3etRc20uL9X48MtL/PP373N3fUA/zvAsne+8t0OuhgPnH70xzyDN6UfDQenVpTo/+dIiz11o8Nb7Hb76+g+4udYnyYcDPkOHmmez2LKIU8WdzfTEuX5WWYbGz/3oRf7nT74wtQPT40rCyEypEFMozQuy3U8ws7wgzQs8Z3ZO1yeJf9b3WXxwDh8rvSiTY0ecO9N2TZxUPOO024sy4mQ4wxcnil6UMddwz2xbaV7QjzLCSBElirKEOCsIkxyAoqjoRhnuTkKlgSpKVFYSxoreQFGUFWUJSVbQT3RyVZJkw23uhBllUaHpGoausROmdKOMQZg/2H5a4NoG2zsp/UiRqQpVFPTjHMNMyIuKQayIE0VelISpYqObYJnDbeWqZKByrN3B9kAN206zgiRVbHZj+mGGZeoH+iDNh/sYRookLbCNgn6YsbodMgiHA+AoUez0C9Js2M9pWbC6E6GhkaQKXdfoxRk7g5T5KGNtJ6QfFhS7M5oVUJSQqYLtAWTpg/0Wj6eKitubAzZ7MRcXamcdzonIM6VCTCHHMrCt4SdctmXgWNP5addxniT+Wd9n8cE5fKw0fFuOHXHuTNs1cVLxjNNuw7fx3OHg0XNNGr59pttyLIO6bxP4Jr5r4jkm7ZpJ4FrYpo5rmzR9m/m2y2LTxTR0fNfk0lJAo2Zi6BqGDq5tUHcNrN11XNukHQzb9RwTzzFoBw5N36YWWLiOgeMY1AMb3zGZazvUfRPb1HAsk7pnsVB3mWs41DwTzzVxbIPAMVlsuiw2PSxTx3EMWg1n1F6r4RD4w9e6jslC06Me2A/1gWMZBK5F4Ju4zrAv64HNxbmAWmDhWAa+a3JxzsOxh7dHO5bBxbZP3TdxneE2Gp5Nu+bQ8G1W2gH1wMDYTcXeTKltGszVLJq1J8v1s8o0NK4s1FhoeGcdyonJ7btiaj3ruZu2Z4nGtZc3eaZ09szSOSfPlB40S7kTDzwub9N2XMszpQ/Mzde4fWdHnimVZ0rlmdITOu72XRmUiqkluZtNkrfZJbmbXZK72SR5m12Su9kluTtbxw1K5fZdIYQQQgghhBBnRgalQgghhBBCCCHOjAxKhRBCCCGEEEKcGRmUCiGEEEIIIc6FvS9VKsup+NocMSYp5CaEEEIIIYSYeWVZcXOtT5YX2JbB9ZX6VHxztXg8mSkVQgghhBBCzLw0L8h2S91keTEqeyOmn8yUCvEUpq2G3FFmIUYhhJhGp1UbU4xvEn+z9rcJPFWd0r1bQ90oG9X13FsexjlhnJPnBVGmWGr5NOvOge0oVbLdS9A0RjVIDU2jF2U4lkHgWWPF+Kj9tAydvCgxNI1BklNzLYqqOtDeOP28t68AnmNOrOb4Xt/1w4w4zTFNnQvzNUxzWJs1jRVv3t3hucU6nmeRFyUN3x7VY82LElvXeftehyRXDEKFaxtYtoG1Dq5tYtsG7bpLWVa8fWubO6sdLMNA18E0dXzbohtnzO/Wc93qxmR5QbPm4O/239424UFt1ywr6McZdc/Gtoe1aNfWBnz9m7dYagW8cLmF7RgstwN0XWNtM2SzHzNfc6n5NnGq2OnF3OuG2KaBZRrUHZvltk9RVexV7tyrNQscqAe7v/6ra5uYpn6iXE0TGZQK8YRm4RaRWYhRCCGmkVIlr7+xRpwoPNfkY6+uyMB0wibxN2t/m6ahgzbM7XHtPyrvZVnx7mqP2+t9avUubd/ixoUGAO/c6fKNt+/z7mqPrV6Kaegstjz+40ev8dLVNrquoVTJ//vvq3z/ZoeSkqWmz4uXm7yz2iXJCixT58duLGCY+iNjfNR+Jqliq5/Q9G3eWe3i2gZJVvDy5Raea3F9ZVgj8nH9vH9fK+Dqcp0bFxqPjeWkOSzLinfudPnvb63x5q0OYaKoeRavXm/zY88v0o1T/s/XfkCSF1DBy5dbNOsOjm1Q8yzeX+/jGBrfu7XDdj8jL8HUwHd0HGc4zDENg4sLAS9dbbG1E/Ht97bZ6SVogGka2JaGhj4c1BngWDrrnQRVVviOwULTZ2XOx3VMWoEDOpiGTpzlvPHuNrkqsS2TDz03h+Pr/O9ffWu0f7apcW25zvMXG9RqFv/fN+7RT3IsQ2ep5dILE+53cvaeftUA24RW4NCs2SR5ga4Nj6WXrrbQNY17W+FoINwNU/qRAiqagcPzl5pcv9AYK1fTRgalQjyho24R8ZzpOqVmIUYhhJhGvSgjToazRHGi6EUZcw33jKM63ybxN2t/m2GSA8NZruPaf1Te07wgTHJyVZLnJWGcj24P7YQpvSgniofbq6qKQZyzuh1x7UIDzzHpRRndXkamCoqyZCdM2eol9MIMjeGs10Y3pu7bj4zxUfuZFyVxotCAMM7R0YYzuKnCMPRRvI/r5/37ym5b48Ry0hymeUEnTNkZZMRZQa5K0rxgoxOz0Y1Z34lIsmF/qqJko5dg2wZxqsjykiQtCFVBlObshkpRQapKSgoMHVRRESY59zYjNrsRqSooSihLKCmpKqgoQddAVQxSyFVJUVbEQDfMCHyLJCswDB3b0BmonChVRMlebDmdMOXercGB/ctVRT/OWO/GbPZiElVQlhVZqdgeZKS5Yv/XMVWAKiBMFegaZVWhayVhqri/HY9mj6NEkauSKClIdmezdT2nvztbP4vv92YrWiGmiGMZ2Lt/NGzLGN0WNE1mIUYhhJhGDd/Gc83RjFnDt886pHNvEn+z9rcZuNaBmdKj2n9U3h3LIHAtLFPHsnQCzxq10QocGr6F7xnEuYFp6NQ8i4tz/ug1Dd+m2bCxtw1KNNqBw3zDZWeQjmZKF5vegZnScftgbz/LshrFHXgWtq0TYBE45oH2HtfP+/e1ggP7Om5/jxO/Yxm0Aod2zWbNHsbvWAaLLY/Fpodt68PZ3rzA0gwWGy6eY45mSgdJhuObrO9YJNnuTKkOrmXg2AYVw5nSwLW4sOBjG7DZTYj0HHP31t3RTKlpPJgpVSVaWeE5Bs3Apu5auI5J3bVAh8C3cDIdf8sYzZS2AoflRY9v/mBntH+WqVH3bJaaHrWaxZ37IbkqsQyDuZpNLywZRAdnSk0DAsek6VujmdLAMVme89A1DaUKfNek7tkYerq7ZkXds6h71ti5mjZatXez8hnb2OhPpN3FxfrE2haTNQu5m4XnNT/oGGchb+JokrvZJbmbjEk/Uyp5e9isPFM6v1Aj6ifyTOkj4jiLZ0ozVRDYFu3asN8zVWCbxoFnSiNVyjOlZ2hxsX7kchmUiqkluZtNkrfZJbmbXZK72SR5m12Su+k0zgcQkruzddygVG7fFUIIIYQQQsy0k36plZgu8jVyQgghhBBCiJl2+Eutwnj4xVZSr3Q2yKBUCCGEEEIIMdP2vmQJIHCt0fO58kWPs0Fu3xVCCCGEEELMNF3XuL5SP9GXWonpIYNSIYQQQgghxMzTde1Afc5Zq9X5LJPbd4UQQgghhBBCnBn5+EA8sx5XR+sk9bmetJbXkziuftos1Ew97yQHQgjxZKbp+qlUSWeQjuqGPqqO517MZVk99Lf5cN3QvX3ba98y9GPrSo7bH4973eNq7R63/km2H8Y5aV7QqjkAo+3puvbQ/gPHvl9KEsXNtR6OZdAOHG5vDpirO7TrLmvbIaoouThfIytLaq5FGOfc34moexa+a5HmBVVVYZo6dc9mu5fQjzOWWj6BZzFIclzT4Na9LmmYYZo6WVaw0Y2xNI0wV1xoBwDc2RpgahqqKpkLXFzXxLUfztPezGxZVmx2YsI4I0kVYaa4slBH07VRTHt1aff3Q1lWZFnB7Y0+SVbw/MUm7YY7qjPbi7Ijj5O9erd5XlBS0Q5cNnsxmgYX5mvY9uw9QyuDUvFM2v+14Ud9VXhZVry72uP2ep8KuLpc58aFxrGD13Ff+7SUKnn9jTXiROG5Jh97dQXT1B+7P2LyJAdCCPFkpun6qVTJP373HqsbIYah8RMvLPL8peZD8RT7YtY1jbWdiCR98LdZ1zVurvVJUsVWP2G+4f7/7d15cGRnfe//z9l679Y2Gmnk8cx4Yow9BhyMiW3AhhCDWcwyBsLqkJDCCbHLoSpsMWa7uCC/SqDqFksFuPwIWyAu5+LwI5eAgbCEGLgxweANbOwZzy7NaOm9z/r7Q6O2pGlpWiO1jlr9fv0zmtbT5zznfPvoPN8+z6JUwtb2LVnd9cBRHRovq1z3tHtbn3ZtKyxoN7R7Pk5Xbqk2w+nev5L9P3xwRvf8dkJeEGl0KCNLs4loMmlpbDArzw+bx5+wLYVRpEMT5VPaS/W6ry98+wE9cqQk3/Mkw5RpWjKNUAPZpKarnoIwVCGb1O/+zhZVG672HynreKkuU9JAPqkwklw/1FA+KcuUjk7VVfcCDeWT2jVaUDpl6aGDM0qnE6pUXI30J/XboyUVy3VVG7NJc9IyJNNQsezKCyTbkpIJS9uGsurPJdWXS8o4eSoMw9CWvrTO2pLVwRNl3X3/MR06XlTDlQxzdoKlbYMZeUGkLYWUtg5kdMHOASUTtsIo0oHxko5NlvXA/klNTLtSJA32JfXcp56trQNpPXRoWrWGr2rDX/A5CcNIP7r3sO7/7Qkdnqwqm7TkB5GCaLZOF+4a0MueubvrElOSUvSk+dOGz00VPn/cQcMLVKl78vxQkprfArYam7CSsqtVrLqq1We/YazVfRWrrgYLqdMeDzqPGADAmdlIfz+LVVfliic/COUH0lS50bI+Dddv1nm60lCp4sqxzea9OZ205XqBvCBUre7Ly4QyjUDHizWVK7PthHojUKnmndJuaPd8nK7cUm2G071/JfufrjRUd2fLnpiqy3ZMZZO2yhVPU4nZp81zx+/5oRpe0LK9dHiyoqmyqyCI5AWRPN9XMiFJkSZKDUVhpEhSqeppstxQsVxXserJDyJFYaipsivbNBVGkWZqvnw/VM31FYZSseLp2HRNg/mUylVfhmGoWvd1dCpUtT5b1vcjmUaoMAgVyVQQSqGkIJRcL1Sx4skwTUWG5Jizib1tmaqnfB2fqWn8RE3VeiDfl/xQMiPJjQJNVRpyTEului+n0lCl4cs/+RS63gg0WWyoXPcVBJEMQ6rUfB2YKCmVtFSszJ6PxZ+TWsPXTNFVueHL9WYrWg982aYlxzJ1oljX8WJNY1tybX/uNwJaTOhJc9OGz30LuHiq8KRjKZty5NimIknZtLPkdOIrKbtahUxC6ZTd/NazkEm0dTzoPGIAAGdmI/39LGQSymUdFauuLMvQQC7Zsj7JhN2sc382qbobNJ+UznVdTZzs1ptO2XJsUwnH0pZCWrmso5mKJS8ZKp92Tmk3tHs+TlduqTbD6d6/kv33Z5NKJSx5QaShgVTzSWku62ggl5Tnh48fv20pnbRVqXmntJfGBrMayCU0VWrICQ05ttXySWk+42gwl1TKMVSrB3KDQKZlaSCXUBjNJsrZpKlE1lIYhqp7gQpZRyP9aaVTlnIZW6mUoyCINNKfnE0SfV92EMqxTSUtSzIN+b6rMJAsU0o4pgpZR30ZR33ZhU9KU0lbW/rSqvuBDhwryrZnE9m5J6UD2aS8IFI+ZWsgm1Q2aSuZsJVO2ipVXQ0Wkjo2ZatWD6VIyqZtnT2cVyGTUCGbmO3iq2jB58SxTPUVEsolbRUdU9mkpWRgNp+UDhVS2lJIt/V530iMKIqiuCshSRMTpY5sd3g437Fto7M6HTvGlHZGL19zGyUGZ6qXY9ftiF13Im6P20h/P9sZUzo8nNexY0XGlG6gMaW1hq+DJ8oKg0iFTFKFjKNKwztlTKkSFmNKYzQ8nG/5OkkpNixi152IW/cidt2L2HUn4ta9iN3GU2v42n/08ZjsHM237HZM7OK1VFLKkjAAAAAAutpct2NJsXcDx8oxphQAAABAVzNNQ7tG8xumGzhWhqQUAAAAQNebG+OJ7kP3XQAAAABAbEhKAQAAAACxISkFAAAAsCnMLdMXhhtigRG0iU7XwElLrSeG3tLO2mwbfa1YoJtw3WwMvRqH+X/P59bVXOk5mH/uJJ2yjcXndr3OdbvtmpXUp52yc2uX1hq+ko7VXGNTWnhu5s59LuXIC0JVq56OTlU1OpBRImHJC2bXIz0+XdNkua5tA1mlUrNrdSYdSw0vUBhGC9YKnS43NFVqKAgjmYahwUJS/bmkfD/U4cmKRvszKlddzZQa8oJQ/bmkTNPQTKmhE6W6xgazKlVc/fbojHYM5RSaUjbpKIii5nl0LFNBFKmQSTTXHZ2La6Xm6dhERQ+NT6k/ldL24ZxyJ8vlUo6qrn9K28H3Qz12uKgHj5xQ3knqvB0DypxcH3fxWq7zE+75v6/XfT02XlI+7Wh4INNyTdqNjqQU0OxFvu9oSfWGrxOluoYKKaUStnaN5nvq5tzr5j4Hrhco4Vgt4+/7oe66/6hqdV/plK3L94zKts223gtgIa6bjaFX4zD/73kyaWlsMKsgjFZ0DuafO9syJWN2u3PbkLTg3O7YmtNj4+WOn+t22zUriX07ZcMw0sMHZ/SLh8d1fKYhyzL0O2N92jGSbyaiCcfS9i1Z/fTBY6rUPFXqnpKOqbvuOybv5HncvjWvQsZRudrQ/vGK6m6gTMLWnnMGtXUwI8c25fmBpkoNDRXS2jaY0eETFdUavmqurydu79dvDk4rk3SUSJnaf6ikSt1TGEV66p5t+s2jx5VLJTQ2nNVQLqk77z6guuvLUqjjRU+VmqvIMHXWUEaRYWgg66juh0o5lhp+qJH+tNIpR+ef3a90ytGOrTk9crion953QD/85YSCkw9p+9K2xoZzGh1My/UjnbUlq2zGabYdfD/Ut36+T1///j554ex7BrKOLj5vq0a3ZLVzNK/d2wrNpPeRw0XtP1bUiWJdQ4W0do7mNTaQ0Ze/+2sdnKjIMg1d8eRtuuKis7ouMSUpBTT7zZ3rBfKCULW6Ly8TyjQCNbyAWdx6yNznQJJcr3X8i1VXtbovSarVfRWrrgYLqbbeC2AhrpuNoVfjMP/vebniaSrRUCGTWNE5mH/uKnVP0ux6mXPbkLTg3Bar7rqc63bbNSuJfTtlG16g6UpDlbqvhufL8AyVap6myw0lT64d6nqBjhdrqtV9BUGkYsWV54equ4EURfIDX1OlhkxDGp9x1XB9haFUcz2dmKkrm3Ie358bqFb39OsD05opN+Q4s4nYiWJdXhDJsgxNnKhruuLKMqVy1dfh8ZLqjUAJK9B0saFSpaFaI1AYStNVT9WGr0izXy7M1FzZpinDkMIgkutHCoNQpbonP4hUafiyLFPFqqvpSkOHTjSaCakk1VxfMxVX6ZQl349U81Iy60az7VCsutp/qCw/fPw91Yav8Zmq+vNJVWpe8zw3vECVuqd6I1DDDVRv+KrUPO33SpqpuAqCSEEQ6chUrbn9btJdKTTQIXMLLjuWqXTKlmObLLzcg9pZeLuQSSidmr0Jp1O2CplE2+8FsBDXzcbQq3GY//c8l3U0kEtKWtk5mH/usilH2bSzYBuLz20hk1iXc91uu2YlsW+nbNKx1J9NKpuylXRspVO28mlH/bnkgnOzpZBWOmXLsgwVsgmNbckolbBkmYZSCVsD+aRyaUdb+xJKJma7qKYTjob6UspmHPUXkspmbCUTs92DM+nZfYVhpJHBtJ567rC2D+dkGoaGh1LqzyZkGoZyGVtjW/NKJWePpb+Q1K6RgtJJS5YlDeUdZZK2DEm2baovnVAm5ag/4yiTtlU4uZ98ylE+m1A2aTfj2p9N6qyhpKx5D4/TCVt92YT6MgnlMwmlHWtB26GQSWjnWTnNf6iZSdra2pdRKmkrm3aa5znpWMqmHKWSlpIJq/n7nVvz6ssmZFmGEo6pbQPp5va7iRFF0YYYBTwxUerIdoeH8x3bNjprvWPHmNK10e3XXC+PKe322PWybo5dt183q7GR4tarcTjTMaXzY8eY0lPLrOeYUkk6Nl1Tve7LDyPt2TmgRMJaENv5Y0qHhnM6erTImNKYDA/nW75OUooNi9h1J+LWvYhd9yJ23Ym4dS9it7GsJLEmdvFaKind/AMGAAAAAGxac08N0b027rNdAAAAAMCmR1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiwzRVQJvanW68VbkzWZNsJfurNXxJOmU9q3a206tr0wEAsBLz77eD4cpWVFxqzdC1uAevZn3ttXAm67wvXtu11vDluoEqDU+DuZSq7uwap0nH0vHpmmYqDeVSjlIpW33ZpLwglKTmjLtz66Kmk7bSSbtZD2nh2qj1uq/fPDapZCQZlrGgjGOZqjV8VWqeoiiSbZvKJh1NluvKpxNKJ22V696Cc+i6gcanq0o7s/s0DEmhdHCyrCiKZMjQ1r60+gupZr0Wn6d63df9+0/oyGRZhXRSQ4W0do/1KZWym2ufNryguabqXF0bXrDgHCzXFuwGJKVAG8Iw0r6jJbleoIRjaddofskEb3E5SW2990z398jhog6MlxRJ2jGS1+5thQXJ8HLbaXc/AAD0ssX321ooDabba/zP3WvrDV8nSnUNFVJKJWzt2JrTY+PlVd2D27mP+36ou+4/qlrdVzpl6/I9o2uWmC51bMsdy/w625apMIr06NEZ3ffIpGzLUBBEGu7PyLYlhZEeODCtYsWVbVnaOZLTyEBGA/mUTMvQ9uGcwjDSL387oRNFVwP5hLYNZjU8kFbCtiRj9vgTjqXR/rS+8u8PyQtC1Ru+nvu725VJO5Ixm1wem6nq6PGKHjtWVsMLNJhLyPUjJRxLtm1oqJBSPpNQNu3o8j2jCsNI/99dj+rwRFlTZVdJ21QQBjo+XVe5FsiPJMeS8pmEnnruFm3fmtdQX0pT5UbzPI32p/W5b/5Kdz803Tw/2ZSlJ+0a0nXPf6IOHq/ont9OyAsijW3J6qyhrDw/1ESxJt8PZRiz50CZ0hIOAAAgAElEQVSSDk2UW7YFuwVJKdCGhhfIPfmNlOsFanhBy/WwWpWb+/l07z3T/VXqnjx/9hvDuW/T5sqebjvt7gcAgF62+H5brrrK2u2tjTl3r/WCULW6Ly8TyjQCFavuqu/B7dzHi1VXtfrsU7Ra3Vex6mqwkFrRfk63/8XHttyxzK9zpT7bbimWPNUagRzbVL3hK51yZZmGSlVXdTdQEEQKAl/TFVeWbcp2LGWTtqbLDbleoErdl+cHKlc9TSYaKuQSKtc8JWxTqYQt1wu0f7ykSs2TZZuq1HwdL9U1NC9vmym6mim7qjcCuX6g6aovzwuUzyYkb/bJeCbpNM9h3fVVrLhq+KGqdV+ebSrwQ9Xd2YRUkoJQqruBxqdrymQSSqfsBedp/3hJx6bcBefH9UIdL9a1f7ykUs1T3Z09V9PFhjJJW0nHUrniSZIc29R0uSFJS7YFu0V31RaISdKxlHCs5jeRc9092i3XznvPdH/ZlCPHNhVJyqadBWVPt5129wMAQC9bfL/NZRJt3zPn7rVhGCmdsuXYphKOpUImoemyu6p7cDv38cK8ZCidslXIJFa8n9Ptf/GxLXcs8+ucTTlKJ21N5etKJy3ZliHHcpRLJWTbUj5la6binnyqaqk/m9CWfEr5lCPTMtSfSyoMIx0+XlbdDZXLOOrPJjRZrCuSlLBn9zX3VDKbduQFobJpW1vyKWVSjz8p7SskVHM9zVRcGabUn7Hl+taCJ6WWZTTPYS7lqJBNqFx1lUnZzSelddeSH8wmppYppRKWtvanNZhLKpO0Z7vfnjxPo/1pjQwkdPB4tXl+Eo6pLYWUdm7N6+DxilIJS14Qqb+Q1EAuKc8Plcs6zSel/bmkpNlktFVbsFsYURStrFN8h0xMlDqy3eHhfMe2jc7aaLFjTGl7Nlrc0D5i172IXXcibt1l/v327O0DmjxRXtF7GVPaus7S2o4plaRHjhTl2KYUSaNDmeZ4zHrdVy2KGFMao+HhfMvXSUqxYRG77kTcuhex617ErjsRt+5F7DYu3w/10KFpRZFajm9dj9gxieTSlkpK6b4LAAAAoOuFYaTHxsuKQimKpOG+tRk7u9I6MInkypGUAgAAAOh68ydRGp+uyg9D5TOJdU0MmUTyzKxtp3IAAAAAiMHcJEpeECrS7Oy081dDWM86SGISyRUgbQcAAADQ9UzT0K7RvGoNX9m001yjdD0Tw7k6MKZ0ZUhKAQAAAGwKpmkom3a0e1shtsTQNNtbwxaP69jZ8jxPN998sw4dOiTXdfWWt7xFf/AHf9Cp3QEAAADocfNnviUx7B4di9TXv/519ff362//9m81NTWlvXv3kpRi3XVyfdBOW496bJRj3SwWn0/OLwBsbCv9Oz1XPgijU15rtY2V/K6T94y13HY79V6qzFJrmc6tARuGkUxzdv3QuXKSWta91T5qDV9HJqvNbrtzExwtXhd1uTrPrQE6V5e5dT/nr1M7Vz/LMJprl86tIVopufrpQ0eUSthK2JZGB7LaOVqQbZvNYwzDSF4Qqj+X1PR0Xd+/51HV65HOPXtQZ4/klTu5FurcfpZa33bxeqXd3NboWFL6ghe8QFdffXXz/5bFIF+srzOZknujTOO9HvXYKMe6WSw+nzu25vTYeJnzCwAb1Ervg/PL10OpLznbtl1qG8ttfz3vGWt5v2+n3ovPyVyZesPXiVJdQ4XUgvVDwzDSI4eL2n+sqBPFugbyKSUdS1v6UkrYlmSoZZLZah+lqqsTxbpGBzPNCY6SjrWg7MBQ7pTzMVfnesPXRLEm1ws0VWpoqJDWztG8do3kte9oSQfGSwqiSAnbUn82oYcPzyiTtJVO2RobzGqiVNWX/s+Dqs/mrjIkZVO2Lr9wVL977hYdnqxoYrqqmYqnfDqhXM7Uv/zgMc19xfG9X41ruJDQOWN92n1Wn1IJW4P5pKbKjeZ5m38+J4o1+X4owzC0YySv3dsKXdvW6FhSms1mJUnlclk33XST3vrWty5bfmAgI9vuTOK61CKt2PhWE7tq3VOm5Cpz8v/5vrQyKWfN39MJ61GPTu6jF6+5xefTSTvKZJOxf5ZWqhdjt1kQu+5E3OKz0vvg/PINN1B+a06SltzGcttfz3vGWt7v26m3pJZlTMfSdM1XJptUKmE361Gte7JPVGXZtkzLUmQYCg1DmWyyud9czl5Q96XqkUon5IZSJptUXy6ps8b61XD9BWWL5fqSdTYdS1GxociIZFqWLMeSnbDlpB3ZSVvJVEKuFyg0JFmmQhlKpROKIik0DR2brMvzHz9fkSQ/CDVRrKnsB7JsW15oKIgk0za1/0hFjz9zn1Vu+Kq4vvwwUmgYcpK2rHrQPG/zz2dUbMhyTCXs2Xp2S1ujlY52tD5y5IhuuOEGve51r9NLXvKSZctOTVU7Uofh4bwmJkod2TY6a7WxC8NI1Uqj+U1YaaamSqm+5u/phPWoR6f20avX3OLz6WWdDfFZWolejd1mQOy6E3GL10rvg/PLj2zNqzRTk6Qlt7Hc9tfznrGW9/t26r34nMyVqTd8Bb6vaqWh0Aua9QjDSP7J34VBICOKZEaz+5l7UjpddRfUfbl6DGQcDecSSjqmDh2elmOZC8oWdg3qt/unWta53vBlKJQRhQqDQIEXyHd9eTVPfsNXo+42n5QqCGUqUr3mKp2yZYaRRgZTcmwpmPek1LZMDRfSytmWin5djhnJMqTQD7VzW1a/+M2JBYlpLmkrm7Blm4bMKJK36LzNP5+GQvleqLofyHeTXdPWaMWIomhxgr4mjh8/ruuuu07vfe97dfnll5+2fKf+KPMHv3utRewYU7r+++jla67bx5T2cuy6HbHrTsQtfmc6pvSssX5NniifdhuMKV3fMaXzx40u7t47t72RkYKOHSsypjQm656U3nrrrfrmN7+p3bt3N1/7zGc+o1Qq1bI8SSkWI3bdibh1L2LXvYhddyJu3YvYbWy1hq/9Rx+Pz87RfHMmXmIXr6WS0o51373lllt0yy23dGrzAAAAAHCKpGMp4VjNJ6VzT0+X0m09mzYjFu8BAAAAsGmYpqFdo/m2Ek1WI9gYzLgrAAAAAABraf540OU0vEDuyXGZc8vIYP2RlAIAAADoSXNdfSW11dUXnUH3XQAAAAA9aSVdfdE5JKUAAAAAetZcV1/Eh7MPrEKrdbekpdfTmlvfaqVrSZ3prHC+H6pYdVU4ud4VlsfsewCA1eA+srTTtYPqdV+PjZeUTzsaHsh0rN0ShJEqNe+UdUjn6rhUu26pNVYXq9d9/eKRcT3y2KTy6ZSe/DvD2j6abx7P4s/I0aNlff+XjymXdpTNJHXhziFtGcws2Mdcey6XchRE0ab8fJGUAmdobra2esPXiVJdQ4WUErYlGbN/PObP4BaGkR45XNSB8ZIiSTtG8tq9rdD2Qt1nMiuc74e66/6jqtV9pVO2Lt8zSmK6DGbfAwCsBveRpZ2uHVSv+/rinQ/q4ERFlmnoiidv0xUXndV2u2V+oie1fjgwV+6hx6b0iweP6kSxrqFCWjtHZ+siacl2nesGzddSCXvJ2Nbrvj799f/WLx55fB3Uf7v7oF56+Tn6g6efLdM0FnxGUjL0ni/8X0UnyxqSBvMH9EfPf6Iu/J0hmabRbM9Vap6qDV9P3N6vdMrZdJ8vklLgDM3N1uYFoWp1X14mlOeHkmYHzc/N4JZO2mp4gSp1r/n7Ss1r/q7d/UhasM3TKVZd1eqz30jW6r6KVVeDhdSZHu6md6bnGQAAifvIck7XDjo8WdFMxVUQRAqCSEemam23W+Z/GWBb5pIPBxpeoDCMVKq6qjcCNdxA9YbfrIukJdt1kpqvmcbSsT08WdHhSXfRsYd65FhRT6+6SiftBZ+R//vQRDMhlaRIUqXh6eEjMzp3R7/SSbvZnguC2Se8lYYvyzI33eeLxybAGZqbrc2xTKVTthzbVDblKJt2JC2cwS3pWMqmHDm2Kds2lU07bc/udqazwhUyCaVTs3+s0ilbhUxipYfYU5h9DwCwGtxHlna6dtDYYFZ92YQsy1DCMbVtIN12u2X+lwGVuqdKzZP0+BcDc0nr/qMlHTlRVSblKJW0lExYSiXtZl3m4mcZhmzLlG0azXbd/LbecrEdG8xqbHBhvZOOqd0jBRUyiVM+I1fs2ab5zzoNSdmko3O39TX3MdeesyxD2bSjbNLelJ8vI4qi6PTFOm9ionT6QmdgeDjfsW2js7ohdowpPVU3xG0pvT4WqJtj1+uIXXcibt1rqdj1+n1kOZ0aU3q6J6UNL9D+o4/H6ql7RnX8eLnlmFLfD/XQoWkFQSTLNPWE7X0yTYMxpWtoeDjf8vXN88wXiMH82drm//Fs1Z3CNI3mU9TV7GclbNuky+4KMPseAGA1uI8s7XTtoFTK1nk7Bs5ou/OXdJEWPhyYezo5N44zPa9X22JeEMqQIdsymv9P23bLtt5yx3HZnjFdtmdsyfrO/4yMjub0mtE9y26zF9pzXDUAAAAAutbiRG/+z4uTVmuZJ4yLE9jN1kV2IyMpBQAAALBptfsEe3ECu1wXWbpqry2SUgAAAACbXhhGqtYfX6O0lXYS2Lklbip1T9mUo91j7S3zh6WRlAIAAADYFJZ6gjk3IVKm5Kpaaaxqnc9aw9eB8ZI8P5Rjm9o2lDnjeUMwi6QUAAAAQNebPxPv/HVKpceXjslobdaRjRb9i9UhKUVPWqtxAHEsudIK4xo2htUs+wMAm1W33KPWq55xnI/NEoP5v5dOXYJv/pqlixPPla4ju3jZv/lLwqSTtnaM5FWpecqmnQXJ7fHjVf3wV/sVuKEGB3J6+vmjKhSSS+5nrrwZWXrS7mGNDGeVTtoLlp+Zf9xhGG2ItudaIylFz1nuW7SV8P1Qd91/VLW6r3TK1uV7RmP547BWx4PVmRtfcmC8pEjSjpG8dm9jjAmA3tYt96j1qmcc52OzxOB065G2Wv5lfuI5N4lRvi+t0kzttJMY7TtaUr3h60SproFcUlPlhoYKKaUStnaNzt7jFyfFx49XdfP/+on85paO6bv/fUjveu3TWiami8t/8+6DesaFo9q1rU/DA2mlErZ2bM3psfGyXC+QaRg6OlVVvRFv27MTNsdRACvQ6lu0M1GsuqrVZ/+M1Oq+ilV3zeq4Emt1PFidhheoUvfk+aF8P1Sl5hELAD2vW+5R61XPOM7HZonB/N9X6p4qNe+UsnOJ587RfMvk2zQNZVLOaZPyuX15Qaha3Ve14atW9+X5YXN/cxMizd/Wfz08Pi8hnVWsenrg4FTL/Swu7wfSgYmKJsuN5r6KVbd53NOVhkqV2fZmnG3PTiApRc9ZafeNpRQyCaVTs50N0ilbhUxizeq4Emt1PFidpGMpm3Lk2KZs21Q27RALAD2vW+5R61XPOM7HZonB/N9nU05zYqFWT0RXO4Rmbl+OZSqdspVJ2kqnbDm2uew5vOTcrad0Qy1kHF2wfaCt8rYlnT2c1WAu2dxXIZNoHnd/Nql8dra9GWfbsxOMKIo2xPjciYlSR7Y7PJzv2LbRWZ2MHWNKO6eXr7luH1Pay7HrdsSuO/VK3DbSPWo5K6nnamLHmNKlrXZMaTvbOF3slhpLOn9M6XLnkDGlyxsezrd8nTGl6EntLqJ8OrZtarCQWoMarc5aHQ9WxzQNpoQHgEW65R61XvWM43xslhgs/v3isqsdP7vU++eSv3aSwC1bMrr29y9oe59LlZ+/r/nHbZrGhmh7rrXuS68BAAAAYJHVjp/tlvG3mxFJKQAAAICut9rxs90y/nYz2vjP8QEAAADgNOZm3z3T8bOrfT/OHE9KAQAAAGwKq519dy1m752b+DAMN8R8sl2BJ6UAAAAANoW4Zxpe7WRLvYqkFAAAAEDX2wgJYavJkrph5uO4cYaAk1byzdpS62TN/3mpNbOWKrMW9Y3720EAAHoN997TW69z1OmEsJ3jmJmq64c/36cwlHZt69eO4dyy69pPTtb0H/cfkBEaOntrn3ZuK6gvnzztOq2b7bNGUgpoZd+szS9rW6ZkSL4fLvh5/jaWKr+ab/Ba1VdS7N8OAgDQSzbCk7mNbj3P0dzsuXP7WsvZc9s5jvHxit79Dz/T3EjS/3p4Sg8fntGF5wzJCyKlU7Yu3zPaTEwnJ2t6z/97l2r+bHnLOKBzz+rTy551js7bMXDKA47N/FljoiNAK1uXan7ZSt1Tpead8vP8bSxVfjXrX7WqL2trAQCwvrj3nt56nqO52XN3jubXPGlr5zh++ptjmj+1URhJh05UdWSqJkmq1X0Vq27z9//9yIQa/uPlg0iaLDV0eLJ6yvY3+2eNpBTQytalml82m3KUTTun/Dx/G0uVX803eK3qy9paAACsL+69p7dZzlE7x3HpeSOanwZbpnTWUEbbBtKSpHTKViGTaP7+qbuHNb93sWVIg/mkxgYzp2x/s5zHpRhRFG2IuYonJkod2e7wcL5j20ZnrXfsGFO6Nrjmuhex617ErjsRt+41P3Yb4d670a3XOWqni+vgUE6HDk+fUV1OdxxhGOnAoaJ+8sABWYatc88e1Pm7BmXbJmNKTxoezrd8nTGlwElz61KdSdmlfm6n/JlqVd+VHAMAAFg97r2nt17n6HQTHYVhpN8enNax8dIZjctc7jiaCXEY6ulP2nHKtgcLqZbvGxxM66XPOm/V++92dN8FAAAA0PVO18W14QVquGszLtP3Q00W6/L9sLntzTzms9M2Z6oNAAAAoKfMTXS0VBfXpGMpOZtDKuFYcixTtYZ/2u6wi7vN+n6ou+4/qlrdb86o28mZf3sBSSkAAACATWG5Lq6maeh3tvcrZUqOZeqx8fJpl1hpNU61WHVVq89Omzs3o+5gIbVsQozl0X0XAAAAQE+wTiatXhC21d22VbfcQiahdGo28Z0/o+5cQkxCunI8KQUAAADQU9rtbtuqnGkaunzP6JIz6mLlSEoBAAAAdL0wjFRrzHarbfXEMgwjVeuewjA67fjTOUuV46no2iIpRc9b7ZpPvh+29U1ZXGtLbYY1rQAA2CzW4768ke/9S9Vtpa+32u4jh4s6MF5SJGnHSF67txUWrOW+72hJmZKraqXRHEPazhIrpmko6VgL1px/5HBRlbqnbMrR7rHH91MsNvTv9+3T/v0z2rF1UM+5eLv6+1svByNJ5bKrn/z2iCozni49b1SJjN2TT19JStHT2llkeTmtZl9r9Udktfs5U3HtFwAAnGo97ssb+d6/VN1W+norDS9Qpe7JO7lES6XmLVindG5saEat1zBdSb2H+1I6MF6S54dybFPbhjLKph0Viw2953/9WCV39n2/2FfWj+47qvdcd0nLxLRcdvU3X/6JDk/NPt395s8O6Oqnna0tQ5kl25SbVe8cKdDCateUajX7Wif2c6ZYMwsAgI1jPe7LG/nev1TdVvp6K0nHUjblyLFN2bapbNpZME70dGuYrrTe0cnfRfPKPXBwSoubguWaq1/tP9Fyuw8fmVGx5jf/7/mRHjteXrZNuVnxpBQ9bbVrSs3Nvjb3pHRu9rW13s+ZYs0sAAA2jvW4L2/ke/9SdVvp662YpqHdYwVtG8pIOnVM6dzY0HxfWqWZ2oqeHi+uR38uqR0jeVVqnrJpp/nE9YLtA8ok1HxSKkm5dEJP3jnUcrvnbutTIW2rfPIBh2Mb2rElt2ybcrMyoiiKTl+s8yYmSh3Z7vBwvmPbRmetV+wYU7q2uOa6F7HrXsSuOxG37rWa2DGmtDNjStt1prFbXI+l6sWY0uUND+dbvs6TUvS8dge5L8W2TQ0Wlv5js1b7OVNx7RcAAJxqPe7LG/nev1TdVvr6emu3HoVCUi+7/InS5e1tN5dL6KqLdq6ydt0v/ggDAAAAQIfMf6q5Vttb6WRSG/np9UZAUgoAAABgU1qcQA4M5VRr+KtKDltNfLTcU9SNPCPyRrE5OysDAAAA6HnzE8h6w9eDj57Q/qMl7TtaUhie2dQ6K53JdyPPiLxR8KQUAAAAwKaUdCzZtqlKzZNtmfJPdqOd6057JuNV52bybbc77kaeEXmjICkFAAAAsHmdfCBqGYYmpmo6PllVOmXrPKv/jDe5kgmYVprE9iKSUgAAAACbUsML5Aehko6lhhdooJCSFWXk2Ka8IFy3pVc2yizCGxVnBlikE7OjhWGkWmN2YeTFizkDAABgbc215xzLbHadzaYcZTMJFYv1Zjfa07X7lvt9q/ZdGEY6NlHRzx46pFol0sVP2KrdO/qXTH7LZVc/e+SoqkVXe3Zs0Y6xwqZdo3Q5JKXAPJ2YHS0MIz1yuKgD4yVFknaM5LV7W4HEFAAAoAMWt+d2bM3JC0I5lqmZxslJhqLZco+Nl5ds9y3XLmzVvts1ktcvf3tcf/+1e+Wf3MZ37zmsvc84R8+/bOcpyWa57Or/+cef6tCkJ0n6P85BvfQZu3TV03f0XGLaW0cLnEYnZkdreIEqdU+eH8r3Q1VqHrOuAQAAdMji9pwXhEonbXlBKM+f7crrB6GKVXfZdt9y7cJW7bti1dVDB6ebCakkBaH04MFpFavuKfV8+MiMirXHS7t+qIcPz7Qsu9mRlALzrHSK73a3mU05cmxTtm0qm3aYdQ0AAKBDlmrPJR1LycTjrxcyiWXbfcu1C1u17wqZhJ6wvX9BV1TLlM7f3q9CJnFKPc/d1qdC+vHSCdvUuWN9LctudkYURWe2QM8am5godWS7w8P5jm0bnRVX7BhTujpcc92L2HUvYtediFv3InYb31LtucGhnA4dnm6+zpjS9TU8nG/5OmNKgUU6MTuaaRrKpp013SYAAABaW6o9Zy16/XTtvuV+36p9Z5qGto3k9LKRJ7ZVz1wuoec+ZUdbZTezzZuGAwAAAAA2PJJSAAAAAEBsSEoBAAAAALEhKQUAAAAAxIakFAAAAAAQG5JSAAAAAEBsWBIG6DJrvY5qJ9ZlBQAA8eL+Hg/XDfTgvhN66MBxpRNJXfzEEW3dkl0yBuWyq98cntZQJqmy62vn1rwsx+y5uJGUAl0kDCPtO1qS6wVKOJZ2jeZX9QdrrbcHAADix/09Hq4b6Lbv3a/v/WKi+dq3fn5Ib3rBBXryE7acEoNy2dX//N/36PhURTUv0lAhpWzK0Ysu26lCLtlTcaP7LtBFGl4g1wskSa4XqHHy542yPQAAED/u7/E4Xqzp0aO1Ba/VXV+/PjTdMgYPH5lRpeYrlBSEkWqup0rd07Gpas/FjaQU6CJJx1LCsSRJCcdS8uTPG2V7AAAgftzf47GlkNY5o+kFr6UStp54Vn/LGJy7rU/ZtC1TkmUaSiccZVOORgYyPRc3uu8CXcQ0De0aza/ZGJG13h4AAIgf9/d4JBKW/vC5e/SU3e2NKc3lEvrLay9iTKlISoGuY5qG0sm1u3TXensAACB+3N/jkUhYesp5W/WU87a2VT6XS+jiNstuZnTfBQAAAADEhqQUAAAAABAbklIAAAAAQGxISgEAAAAAsSEpBQAAAADEhqQUAAAAABAb5okGFgnDiHW9AAAAsGJhGGniRFX3HTyhjGPr7OG8Roaysu3WzwKrVU+PHi1qdCCjvnxSXhD2ZBuUpBSYJwwj7TtakusFSjiWdo3me+6PAgAAAFYuDCP98uHj+vy/3adiNZRhSMP9KV118dl69lPPOiUxrVY9feob92qy6CrhGHrGnm3aNpxVKmH3XBuU7rvAPA0vkOsFkiTXC9Q4+TMAAACwnIYX6JEjM3K9UJGkMJIqNU+HT1RUrLqnlN8/XlK56iuKItXqgQ6fqMjzw55sg5KUAvMkHUsJx5IkJRxLyZM/AwAAAMtJOpZ2b+tTwjFlSDINKZt2NDaUVSGTOKX8zq155TK2DMNQOmVpbCgrxzZ7sg1K911gHtM0tGs0z5hSAAAArIhpGnrKuVv0rtc+va0xpZmMoz+75kmMKRVJKXAK0zSUTnJpAAAAYGVM09DIcFYjw9m2ymcyji7cPdT8/1ITIm12vXnUAAAAAIANgaQUAAAAABAbklIAAAAAQGxISgEAAAAAsSEpBQAAAADEhqQUAAAAABAb1r0AAAAAgDXg+6GOHq9ovFhVXzqps0fySiSsJcuXy65+c3haZ2/Jaag/3XPrk84hKQUAAACAVfL9UN+954C+/7ODKtZ9pRxLv3f+iPZeubtlYlouu/qf//sezZRdJRKWXvv7T9AF5wz2ZGJK910AAAAAWKVi1dWxiZrqfqAwjOT5gY5NVXW8WGtZ/uEjM6rUfEWSGm6gR48V1fCC9a30BkFSCgAAAACrVMgkNDKcVsq2ZJqGHNvSyEBGWwrpluXP3danbNqWISmZsHTOSEFJZ+muvpsZ3XcBAAAAYJVs29QfXHS2LjxrsK0xpblcQn957UWMKRVJKQAAAACsCds2tX00r+2j+bbK53IJXXze1g7XauOj+y4AAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2LAkDAAAAACsAd8PdXi8rEPTZW3JpbVztLDkOqWSVK/7emy8pHza0fBARrbdm88MSUoBAAAAYJV8P9Sdv3hM37nrMVUagVIJW5ftGdW1V+5umZjW676+eOeDOjhRkWUauuLJ23TFRWf1ZGLae0cMAAAAAGusWHV15FhVrh8qiiTXC3RksqzjxVrL8ocnK5qpuAqCSK4X6shUTcWqu8613hhISgEAAABglQqZhLaNZJSwTRmGlHAsbRvMaUsh3bL82GBWfdmELMtQwjG1bSCtQiaxzrXeGOi+CwAAAACrZNumnve7O3Th2GBbY0pTKVvXPe98xpSKpBQAAAAA1oRtm9oxVtCOsUJb5VMpW+ftGOhwrTa+3kzFAQAAAAAbAkkpAAAAACA2JKUAAAAAgNiQlAIAAAAAYkNSCgAAAACITUeT0nvuuUfXXXddJ3cBAAAAABuG64uZWwIAAA+kSURBVAY6fLws1w1OWzYMI1Vqnio1T2EYrUPtNqaOLQnzmc98Rl//+teVTrdeLBYAAAAANhPXDfSNn+xTpeYpm3Z0zWW7llynNAwjPXK4qAPjJUWSdozktXtbQaZprG+lN4COJaU7duzQxz72Mb3jHe9oq/zAQEa23TpgqzU8nO/IdtF5xK47EbfuRey6F7HrTsStexG77tXJ2B04VlQoQ+l0QqGk0DaX3F+17sk+UVUylZAk2Qlb+b60MimnY/XbqDqWlF599dU6ePBg2+WnpqodqcfwcF4TE6WObBudRey6E3HrXsSuexG77kTcuhex616djp3phzIVNZ+Umn645P7CMJLf8NWou4ok+W5SpZmaKqV6x+oXt6US9I4lpQAAAADQSxIJS9dctkvHizVtKaSX7LorSaZpaPdYQduGMpKkdNLuya67EkkpAAAAAKyZRMLS2JZcW2VN01A23XvddRdjSRgAAAAAQGw6mpRu375dt912Wyd3AQAAAADoYjwpBQAAAADEhqQUAAAAABAbklIAAAAAQGxISgEAAAAAsSEpBTapMIxUa/gKwyjuqgAAAPSMMIxUqXmq1Ly22mG02VinFNiUwjDSvqMluV6ghGNp12i+ZxdjBgAAWC9hGOmRw0UdGC8pkrRjJK/d2wpLtsNos83iSSmwCTW8QK4XSJJcL1Dj5M8AAADonIYXqFL35PmhfD9UpeYt2w6jzTaLpBTYhJKOpYRjSZISjqXkyZ8BAADQOUnHUjblyLFN2bapbNpZth1Gm20W3XeBTcg0De0azavhBUo6Vk92AwEAAFhvpmlo91hB24YykqR00l62HUabbRZJKbBJmaahdJJLHAAAYD2ZpqFs2llR+V5vs9F9FwAAAAAQG5JSAAAAAEBsSEoBAAAAALEhKQUAAAAAxIakFAAAAAAQG5JSAAAAAEBsSEoBAAAAIAZhGKnW8BWGUdxViVVvL4gDAAAAADEIw0j7jpbkeoESjqVdo3mZphF3tWLBk1IAAAAAWGcNL5DrBZIk1wvUOPlzLyIpBQAAAIB1lnQsJRxLkpRwLCVP/tyL6L4LAAAAAOvMNA3tGs2r4QVKOlbPdt2VSEoBAAAAIBamaSidJCWj+y4AAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAAACA2BhRFEVxVwIAAAAA0Jt4UgoAAAAAiA1JKQAAAAAgNiSlAAAAAIDYkJQCAAAAAGJDUgoAAAAAiA1JKQAAAAAgNiSlAAAAAIDY2HFXYDXuuece/d3f/Z2++MUv6r777tP73vc+JRIJXXDBBXr3u98t0zT14Q9/WHfffbdM09Q73/lOPe1pT9Pk5KTe9ra3qV6va+vWrfrwhz+sdDod9+H0lDON3fT0tK6++mqdd955kqSrrrpKb3zjG2M+mt7RTtxuvfVW/fznP1c2m9Xb3vY2XXTRRdq/f7/e9a53yTAMPeEJT9D73vc+mSbfia2nM43dfffdpz//8z/Xrl27JEmvfe1r9aIXvSjeg+kRnufp5ptv1qFDh+S6rt7ylrfo3HPPbXktffzjH9f3v/992batm2++WU95ylO47mK02thx3cVjJXGTpP379+uGG27QN77xDUmifRmj1caO9uUGEHWpT3/609E111wTvepVr4qiKIr27t0b3X333VEURdFHP/rR6I477ogeeOCB6FWvelUUhmH06KOPRnv37o2iKIo++MEPRv/8z/8cRVEUfepTn4o+97nPxXIMvWo1sfvxj38c/Y//8T9iq3svaydu3/ve96I3velNURAE0YkTJ5px+7M/+7PoJz/5SRRFUfSe97wn+va3vx3PQfSo1cTutttuiz772c/GVvdedvvtt0e33nprFEVRNDk5GT372c9ueS3de++90XXXXReFYRgdOnQouvbaa6Mo4rqL02pjx3UXj3bjFkVR9LWvfS3au3dv9IxnPKP5ftqX8Vlt7Ghfxq9rvzLdsWOHPvaxjzX/f+zYMV188cWSpIsvvlh33323tm7dqlQqJdd1VS6XZduzD4bvvvtuXXHFFZKkK6+8Uv/5n/+5/gfQw1YTu3vvvVf33Xef3vCGN+imm27S+Ph4LMfQi9qJ28MPP6wrrrhCpmlqcHBQlmVpYmJC9913n37v935PEtdcHFYTu3vvvVff//739frXv14333yzyuVyXIfRc17wghfoL//yL5v/tyyr5bV0991361nPepYMw9DY2JiCINDk5CTXXYxWGzuuu3i0GzdJ6uvr05e+9KUF76d9GZ/Vxo72Zfy6Nim9+uqrm4mKJJ199tn62c9+Jkn693//d9VqNdm2LdM09cIXvlB/8id/oje96U2SpHK5rHw+L0nKZrMqlUrrfwA9bDWx2717t2666SZ96Utf0lVXXaVbb701lmPoRe3E7YILLtCPfvQjeZ6nAwcO6OGHH1atVlMURTIMQxLXXBxWE7unPOUpesc73qEvf/nLOvvss/WJT3wirsPoOdlsVrlcTuVyWTfddJPe+ta3tryWyuWycrncgveVSiWuuxitNnZcd/FoN26S9Pu///vKZDIL3k/7Mj6rjR3ty/h1bVK62Ic+9CF96lOf0vXXX6+hoSENDAzojjvu0JYtW3TnnXfqu9/9rj7+8Y/r2LFjyuVyqlQqkqRKpaJCoRBz7XvbSmJ32WWX6dJLL5UkPe95z9P9998fc+17V6u4PetZz9Ill1yiN77xjfrc5z6nCy+8UP39/QvGsXHNxW8lsXve856nJz3pSZK45uJw5MgR/dEf/ZFe9rKX6SUveUnLa2n+PW3u9Xw+z3UXs9XEjusuPu3EbSm0L+O1mtjRvozfpklKf/CDH+hDH/qQPv3pT2t6elrPfOYzVSgUlMlkZFmWstmsEomEKpWKLr74Yv3gBz+QJP3whz/U0572tJhr39tWErtbbrlF3/rWtyRJd911ly688MKYa9+7WsXt0Ucf1dDQkP7xH/9Rb37zm2UYhgqFgvbs2aOf/vSnkmavuUsuuSTm2ve2lcTuT//0T/XLX/5SEtfcejt+/Lje9KY36e1vf7te+cpXSlLLa+niiy/Wf/zHfygMQx0+fFhhGGpwcJDrLkarjR3XXTzajdtSaF/GZ7Wxo30Zv66efXe+nTt36vrrr1c6ndall16qZz/72QqCQD//+c/1mte8RkEQ6CUveYl2796tt7zlLXrnO9+p2267TQMDA/rIRz4Sd/V72kpi91d/9Ve6+eab9ZWvfEXpdJruFTFqFbdGo6Ef/ehHuv3225VMJvXe975XkvTOd75T73nPe/TRj35Uu3fv1tVXXx1z7XvbSmL3/ve/Xx/84AflOI62bNmiD37wgzHXvnf8/d//vYrFoj75yU/qk5/8pCTp3e9+t2699dYF15JlWbrkkkv06le/WmEYct1tAKuNHdddPNqN21JoX8ZntbGjfRk/I4qiKO5KAAAAAAB606bpvgsAAAAA6D4kpQAAAACA2JCUAgAAAABiQ1IKAAAAAIgNSSkAAAAAIDYkpQAAzFMqlXTDDTdIko4dO6Y3v/nNa76PIAh04403qlarrfm2pdmF4m+88UYFQdCR7QMAsJZISgEAmGdmZkYPPPCAJGlkZESf+cxn1nwfX/nKV/SsZz1L6XR6zbctSdlsVpdffrm++tWvdmT7AACsJZJSAADmufXWWzU+Pq4bbrhBBw8e1HOf+1xJ0rve9S594AMf0Ktf/Wq98IUv1J133qkbb7xRV111lf7mb/5G0uwT0A9/+MPau3evXvrSl+of/uEfTtl+FEX64he/qBe/+MUql8u69NJLVS6XJUkHDx7Ui170IknSHXfcob179+plL3uZbr75ZjUaDUnSl770Jb3qVa/SNddco7179+qRRx6RJD33uc/VW9/6Vl199dU6ceKErrnmGn3hC18Qy5EDADY6klIAAOa55ZZbtHXrVn3iE5845Xfj4+P6p3/6J11//fX667/+a33gAx/QHXfcodtuu02lUkm33XabJOlrX/uabr/9dn33u9/Vf/3Xfy3YxoMPPqh8Pq98Pq9cLqfnPOc5+rd/+zdJs4noy1/+cj300EO67bbb9NWvflX/8i//oqGhIX32s59VuVzWd77zHX3xi1/UN77xDT3nOc/Rl7/85ea2r7zySn3rW9/S0NCQ+vr6lMlk9Otf/7qDZwsAgNWz464AAADd4sorr5QkjY2N6QlPeIKGhoYkSf39/ZqZmdFdd92lBx54QD/5yU8kSdVqVb/+9a91ySWXNLexb98+jY6ONv//ile8Qh/72Mf0yle+Ut/4xjf0+c9/Xnfeeaf279+vP/zDP5QkeZ6nPXv2KJfL6SMf+Yj+9V//Vfv27dOPfvQjXXDBBc1tXXTRRQvqOzY2pn379un888/vzAkBAGANkJQCANAmx3GaP9v2qbfQIAj09re/Xc9//vMlSZOTk8pmswvKGIax4L1Pf/rTNT4+rm9/+9vavn27RkZGFASBXvjCF+qWW26RNDtxURAEOnLkiK677jq94Q1v0JVXXqktW7Y0x79KUjKZXLAvy7JkmnSKAgBsbNypAACYx7Zt+b5/Ru+97LLLdNttt8nzPFUqFb3uda/TL37xiwVldu7cqUOHDjX/bxiGXv7yl+vWW2/VtddeK0m69NJLdeedd+rEiROKokjvf//79fnPf16/+tWvtHPnTv3xH/+xnvzkJ+s73/nOsjPsHjp0SDt27DijYwEAYL2QlAIAMM/Q0JDGxsZ03XXXrfi9r3nNa7Rr1y7t3btXr3jFK3Tttdfq0ksvXVDm/PPP19TUlEqlUvO1F7/4xarVarrqqquaZW688Ua98Y1v1Itf/GKFYajrr79ez3zmMxWGoV70ohdp7969Ouecc3Tw4MGWdSkWiyqXy3TdBQBseEbEtHwAAKyrL3zhCzJNU294wxsUhqG+8pWv6NFHH212110Ln//852Xbtl7/+tev2TYBAOgEnpQCALDOXvva1+rHP/6xarWabrzxRt1+++36i7/4izXbfqVS0V133aVXv/rVa7ZNAAA6hSelAAAAAIDY8KQUAAAAABAbklIAAAAAQGxISgEAAAAAsSEpBQAAAADEhqQUAAAAABCb/x+WBaveV1Y75AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x864 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.style.use('seaborn')\n",
    "plt.figure(figsize=(16, 12))\n",
    "plt.plot(time, mags, marker='.', linestyle='none', alpha=0.25)\n",
    "plt.xlabel('time (year)')\n",
    "plt.ylabel('magnitude')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot indicates a distinct increase in the frequency and magnitueds of earthquakes after 2010.\n",
    "\n",
    "For further analysis, we need to compare the estimate and inter-earthquake times pre-2010 and post-2010 for earthquakes  of magnitude 3 . We also need to include the 95% confidence intervals of the means."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1980 through 2009\n",
      "mean time gap: 204.61 days\n",
      "95% conf int [139.93, 277.13] days\n",
      "\n",
      "2010 through 2017\n",
      "mean time gap: 1.12 days\n",
      "95% conf int [0.97, 1.29] days\n"
     ]
    }
   ],
   "source": [
    "# Create arrays of the pre-2010 and post-2010 time gaps between earthquakes: dt_pre, dt_post\n",
    "dt_pre = pd.read_csv('data/pre2010_timegaps.csv', index_col=0)\n",
    "dt_pre = np.array(dt_pre['0'])\n",
    "\n",
    "dt_post = pd.read_csv('data/post2010_timegaps.csv', index_col=0)\n",
    "dt_post = np.array(dt_post['0'])\n",
    "\n",
    "# Compute the mean inter-earthquake times: mean_dt_pre, mean_dt_post\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 means: bs_reps_pre, bs_reps_post\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 intervales: conf_int_pre, conf_int_post\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(\"\"\"\n",
    "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 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": [
    "There is almost a 200-fold increase in earthquakes after 2010\n",
    "\n",
    "## Hypothesis testing\n",
    "\n",
    "There is emperical evidence of a massive increase in earthquake frequency after wastewater injection began in Oklahoma.\n",
    "\n",
    "H$_{0}$: Assume that the inter-earthquake times have the same mean pre and post-2010. What is the p-value associated with this hypothesis.\n",
    "\n",
    "Test statistic: the mean of pre-2010 time gaps minus the mean of post-2010 time gaps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p = 0.00\n"
     ]
    }
   ],
   "source": [
    "# Compute teh observed test statistic: mean_dt_diff\n",
    "mean_dt_diff = mean_dt_pre - mean_dt_post\n",
    "\n",
    "# Shift the post-2010 mean to have the same mean as the pre-2010 data: dt_post_shift\n",
    "dt_post_shift = dt_post - mean_dt_post + mean_dt_pre\n",
    "\n",
    "# Compute 10,000 bootstrap replicates of each mean of dt_pre\n",
    "# and dt_post_shift: bs_reps_pre, bs_reps_post\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",
    "# Compute replicates difference of means: bs_reps\n",
    "bs_reps = bs_reps_pre - bs_reps_post\n",
    "\n",
    "# Compute the p-value\n",
    "p_val = np.sum(bs_reps >= mean_dt_diff) / 10000\n",
    "print('p = {0:.2f}'.format(p_val))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "0 out of 10,000 samples had a test statistic greater than what was observed.\n",
    "\n",
    "# Earthquake magnitudes in Oklahoma\n",
    "\n",
    "Do the post-2010 earthquakes follow the Gutenberg-Richter Law. We need to study the magnitudes of the pre and post fracking quakes\n",
    "\n",
    "## EDA: Comparing magnitueds before and after 2010\n",
    "\n",
    "We need to make an ECDF of earthquake magnitudes for both periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe8AAAFXCAYAAACLEMbVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XdgU1XfB/Bv0jRNm+42lLbQUipT9lBQCioiiFt4AIEqoggIooKCqFCQIdOFICAP4IO+Ci6mgiAoQ2QXQShgKWXF0l2a0TRN3j9KA6FtOtLbm5t+P/+Yk5Pc/O4x5Zdzzz3nyKxWqxVEREQkGXKxAyAiIqKqYfImIiKSGCZvIiIiiWHyJiIikhgmbyIiIolh8iYiIpIYhdgBVFZ6+nWxQyglKMgH2dl6scOQNLah89iGzmMbOo9t6Lzb21Cj8Sv3tex5O0Gh8BA7BMljGzqPbeg8tqHz2IbOq0obMnkTERFJDJM3ERGRxDB5ExERSQyTNxERkcQweRMREUkMkzcREZHEMHkTERFJDJN3Dfj775MYO/YlW/nMmSSMGPEsXn75RXz44TxYLBYAwKJFH2DEiGfx0kvD8NdfiQCAnJwcvP76GLz88ouYOnUyjEaj3bELCgqwadN6AMB//7sM69d/J/j5jB37ElJTL5Rbn5+fj4kTX8fYsS9h5MjncfLkXwCAkydPYMSI5zB69HCsXLkcAGCxWDB//myMHPk8xo59CZcvXwIAHDz4J55/fjCeeeYZrF69QvBzIiJyJ4Im7+PHjyM+Pr7U8zt37kS/fv0wcOBArFu3TsgQSjGazEi+mgujyVwjx/vqqy8wd+4MmEwm23Pz5s3CuHETsGTJCqjVvti+fSvOnTuLEyf+wvLlX2DKlPfw0UcLAACrV3+OXr36YMmSFWjSpBk2bPje7vhZWZm25O0q1q79Cp06dcanny7HO+8k4IMP5gIAFix4H9OmzcKSJf/FqVMnceZMEvbs+Q0mkwnLlq3CqFGv4NNPP4TFYsGcOTMwc+Y8fP3117h4MRXHjyeKfFZERNIh2PKon3/+OTZu3Ahvb2+75wsLC/H+++/ju+++g7e3N5555hncf//90Gg0QoViYzSZMeOLw9Bm6hEe4oMpz3WCSulcE0RGNsCsWfMxY8ZU23Pp6dfQunVbAEDr1m2xd+/vuPvue6BSqWAymaDT6aBQFH/uX38lIj7+eQBAly73YPnyxRg4cIjtWP/730pcuJCCVas+BwDs2bMbu3b9itzcXLz44ih069Yd/fo9iujoRoiOjsGAAc9gzpwZMJvNkMlkePXVN9CkSVM8/nhvbNy4DQCQkDAZTzzRD3fe2QozZiQgMzMd9eqFITHxGDZs2AoAWLlyObKzs2AwGDBt2ixERjawxTRgwGAolZ4AALO5CEqlF3S6fBQWmmyvu+uurjhy5CAyMzNw991dAQCtWrVGUtJp5ObmwM/P3/ba1q3b4q+/EtG2bTun/l8QuSujyYwrGTpEhqqd/jerzOObC6DVpSFcHQaVwqta7z+XmQGV2a/M9zt7fFdQ0TnU9jkK1vOOiorCokWLSj2fnJyMqKgoBAQEQKlUomPHjjh8+LBQYdi5kqGDNrN43Vhtph5XMnROH/O++3raEnGJiIhIHDt2BACwb98eGI0GeHh4QCaTY8iQ/njttZfxzDNDAQA6nQ6+vr4AAB8fH+Tn59sd69lnh6NRoxg8//wIAIBGo8HHH3+GcePG2y6hX7uWhoSEmXj11QlYvPgj9O8/EIsXf45XX52AOXNmlBv7hg0/IiIiAp99thLDh49EdnaWre6ee7rhk0+WokuXe/Dbb7/avc/Pzw9eXipkZmZgxowpGDlyDHQ6HXx81LbXlJyLTqeDWu1re14ul8PPzx8FBUakpl5AUVER9u/fB6PRULkGJ3Ixey7tx5u/T8OeS/vLrE/KPIcFhz9FUua5ah3DaDJj+pq9eP+n9Zi+Zm+1rhpe02dgY/LPuKbPKH18cwHmHvoYC458irmHPobRXFClYxvNBZh3eBHe2TEP8w4vKvX+kvoFRz4ts14KKjoHMc5RsJ537969cfny5VLP5+fnw8/v5mLrarW6VMIqS1CQj9Nr5/r6e6NBPV9cvpaPBvV80bZ5fXh7OdcEGo0fCgrU8PT0sC0iP3/+XMyaNQvfffd/aN26Na5fv469e3cgIiIMa9ashk6nw+DBg9Gjxz0ICPCHt7ccISF+yMyUISQkyG4x+luPrVZ7ISqqHTQaP9xxRxQsFjM0Gj8EBQXhjjsaAgAuX76Inj27IyDADxpNJ2RkXING4we5XGY7rlKpQGCgD9LSLqN79+7QaPyg0bRGcHAwNBo/KJUKdO3aCRqNH6KjI5GRkVFqgfwzZ85gwoTxmDhxInr06IH8/HyYTEbb6+TyItSvHwq53AKFwmJ7XiYDwsODsHDhAixcOA/+/v5o3rwJ6tcPc7gIPznGtqu+n5J+xbrfN2FAq8fQt3nPUvWrDq/F1uTf0Cf2PjzfaaBd3fZzu/HNuR8BAN+c+xF+fir0atLdVn9CexqLjhdfNVt0/HNM6T4OrcNbVOkYh8+lIjf6Z3jJrci1nEKavg06RUZX+vz+vX4N03fOAwBsS92FT/pOR32/erb6v9Ou4pqhOKlfM2QgT56JhppmlT7+ucwMpOmvAQDS9NdgVF5Hw5DQStdLQW2eY2X/lmt9VzFfX1/odDd7vDqdzi6Zl6emdqt5e2gH2+Wn/DwDKv7ZUD6Nxg/p6deRlaVDYWGRbeezLVu24Y033kFoqAYffjgPXbrcg9zcXMjlnsjK0qOoyAKZzAOXL19DixatsHnzNvTt+xh++mk7mjdvZbeDWna2ASZTIdLTr0OnK4BKVYD09OvIztbDZDLfeK3M9p4GDaKwc+cedOvWA+fOnUFgYDDS06/DZCpEamoaPD09kZR0Bjk5ekRERGHfvgNo2/ZuXLlyGdnZ2Tdea0Z2th7p6deRn2+ETldgF1NKynm8886bmD79fTRp0tRWJ5N5IDHxNCIiIrFz5294/vmX4O8fgl9++RWdO8fh5MkTaNQoFunp17F16w7MmrUQDRtqMGLEKNx3X2+X3DlOCkq+h1S2xUdX4FTOWbQMbIoxHV60q9uZugffJ28CAKw+/h2u5xfggeg4W/260+vxu/YPAMDPyb9BrzdhQIsnbfX/S/zB7nj/S/wB7QLb28pfJtrfr/Jl4nq8oWhg99xXxzeUKt96jAt5ZyGTWwEAMrkVF/LOIjo9uHInD2BL8m/25b9/w+OxD9vK2Tn2/7Zm5+iRLq/890ll9kOYTz2k6a8hzKceVCb772NF9VJQW+d4+9+yo0Re68k7NjYWqampyMnJgY+PDw4fPowXXnih1j5fpVQgNiJA0M9o0CAKb7zxKlQqFdq374iuXbuhqKgIJ04cx6hRw1FUVISHHuqDqKhGeO65FzBz5jRs2vQjAgICkZAwy+5YQUFBKCw0Y8mST+DlVfE4ypgxr2Hu3Jn4+usvYTabMXnyFADAgAHPYOTIYYiIiET9+uEAgEcffQKzZk3HmDEjUL9+fSiVykqd37Jln8JkMuHjj4tvuvP19cWcOR/gjTcmY/r0d2GxWNC58924885WsFha4tChAxg1ajisVivefjsBAKDR1MPo0S/A19cH99/fC40bx1aydYnsOUrOJXUAcCrnLBYfXWH3mvXJP9m9fn3yT3bJuyRx31q+NXl7W9QoQIFd+VaPxvS29bxLyrd7vHFvW8+7pHyrdmEt8X3yBhRZi+Ah80C7sJaljuFIl/DO2Ja6y658q2j/htB4hyLdkAGNdyii/RtW6fgqhRcmdnoFRuV1qEylx7xL6qU85l3ROYhxjjKr1WoV6uCXL1/G+PHjsW7dOmzatAl6vR4DBw7Ezp07sXjxYlitVvTr1w9Dhgyp8Fiu+EtN6j2eEyeOw2Aw4K67uuDSpYuYMOEVrFu3oeI31iCpt6ErcPc2/CVlF366sAN9Gz2Ih2Lut6u7NTkDKJXAx+ycWOp4ix+YZ3s8budkFKHIVvaABz554H1befnhdUjMPQyZDLBagXYBnfBSpwG2+qT0ZCz6axkgA2AFXmkzEs019j9EkzLPYXPKNjwa0xvNQ5qUeY57Lu3HxvPb8Hjj3ohr2LVUfW5BHv7OTMKdIc0R4OVf5jEcuabPwJ/aQ+gS3hn1fEpfzq2Jm63c/XtYG6rS8xY0edckV/xSSP3LmpmZgWnT3oHZXAiz2YwXXhiFLl3uqdUYpN6GrsCd2/CXlF3YkPKzrfxEzMN2Cbyi5FxR/Z9XDmPNmZvTVeObDUCXyE62stFkxoQfP4c5KAWK7BgsfGpEqbu9k9KTsSl5Bx6LfbBU4q5L3Pl7WFtc+rI5uY6QkFAsWrRM7DCojnPUM701cZeUb+99O/J0o6fwfcqPtp5zv5in7OpLEvXmC1vxaKM+dokbKB5mW/jUCIfTtJprYut00iZxcIU1IhJNUuY5LDr+OVLyLmLR8c8dTqcqS7ug1kDJtUPrjfItGqlawni2JcxGTxjPtkQjVenx4i6RnbDsybmlEneJkvtkhJhfTVRdTN5EJJolx1c6LFfkgbDHYbgSCrMZMFwJxQNhj9vVR4aqEa5ojsK/eiJc0RyRoepyjkQkLfwpSUSCcnQz1K03i5VVDlc0wNXCy7bL3hGe9tOsIkPVCC/oDu3R4lUTb0/OKqUCU57rJOjqZERi4DeZiARjNBdg4p5ptmlO8+KmVelu5rEdRuLNLfMhD8qDJdsfYx8ZaVdfmeRcG9NDiWobL5s7obCwEDNmTMHLL7+IESOexd69vwMALl++hNGjX8DLL7+IBQvet+0qVlIXH39zqom77ypW4vad1xISJmPs2JcQHx+P/v0fQ0LCZEHOhcQ1a/8CFFmLe9NF1iLM2r+gSu/PzDOiKOUeFB7tg6KUe5CZZyz1Go5JU11U55K30VyAlNyLNbL27MaNG+HvH4glS1ZgwYJP8MEHxVNQirf+HI0lS1bAarViz57ipL516xYkJLyN3Nxc2zHcfVcxoOyd16ZPfx+ffrocn376KXx9/fDKKxNEORcSVlZhrsNyEEJQMlnVai0u3yoyVI3wEB8AKPOyOFFdVaeSd00vHt+nTx+MGDHKVvbwKP7lf+ZMEtq37wigeKeww4cPAgD8/Pzx6af2PdK//kq07bp162tLlLWr2KuvjsawYYOxd+9uAEC/fo9i/Pix+PjjhdBqr+LVV0djzJgRGDv2JZw7V7yAxeOP31y1KSFhMo4ePYyCAiPefXcSRo8efmOnsT6216xcuRzjxo3CiBHP4coV+zXqBwwYjCeeeBpA2buKyWQy265iwM2d18qyaNEi9O8/AKGh0lrrmGrGkNgRMGZ4w2wGjBneGBI7wq6+5LL4O892rJFdAIncRZ1K3lpdmt3i8VpdmlPHU6vV8PFRQ6/X4d13J2HEiNEAAKvVCplMBgDw8VFDpyteQf3ee+NKbZHq7ruKAWXvvAYA2dlZ2L9/Px5++LFyYyT3FhPuj7DcPig82gdhuX0QE1569TBeFicqrU79NYSrw+wWjw9Xhzl9zLS0f/H222/iqaf646GHinuucvnN30R6/c3kXBa1Wg29Xg8vLxX0en2Fm7Q0a1a8I1FISKhtfDwgIBABAYEAgAsXLqBt2w4AgCZNmuHatdI/UEouU6ampuDuu4tXVIuOboTAwKAyPicEmZmZpY6RnPwPEhLexpgxr6J9+47Q6fJhMNzc4ECv18PX1/G57Nr1Kx599FF4eDi3WxyJZ9YfC3DVeA0Rqnp45543qvx+lVKBhGGdeTc4URXVqZ53yeLxb3Qci4mdXnF68fiMjAyMHz8Wo0e/gkcffcL2fJMmzXD0aPEe5X/++Qfatm1f3iHQunVb7N+/z/baNm3a2dXLZHJYrZZbyqWPceuPhUaNGuGvv44BAM6dO4Pg4OIxRLPZDL1ej8LCQqSkJAMAGjeOtd1sduXKZeTm5tzyOWV80A0pKecxZcokJCTMRNeu9wIA1GpfKBSeuHLlMqxWKw4e3O/wvAHg8OGD6N69u8PXkOsqSdwAcNV4DbP+qNrNaCXYsyaqujr316JSeCEmIKpGjrV06VJcv34dq1evwOrVKwAACxd+grFjX8O8ebOwbNliREc3wn33ld4juIS77yrmyMWLqWjYsCEKhN+3ngRQkrjLKwOAGn7Q4bpdmYicx41JnCD1hfi5q5h7EKsNK9r0AwDScvMwbe9MyL0BiwGY1u1dhAVUfVcsofF76Dy2ofO4MQlVSkREJKZNewerVi2H2WzG+PGTxA6J3ExWjhkFfz8ImXc+rAZfZLUyI4zrpRA5jcm7DuOuYiQ0U2ERYFHAqgu8WSYip9WpG9aIqHYpPT0clomoetjzJqJyLT22Eieyk9A6qDlGtR9e5ffHhPsjLNgbaVkGhAV7lzmPm4iqjsmbiMpUkrgB4ER2EpYeW1nlBM553ETC4F8SEZWpJHGXV64s7upFVPM45k1E1dI9+CG7TUW6Bz8kbkBEdQiTNxFVS+874mD6pznMBR4w/dMcve+IEzskojqDl82JqFq0mTpYshvBkt3IVg70dW7JYSKqHPa8iahabp+zzTncRLWHyZuIqoVzuInEw+RNRNVSMocbAOdwE9UyjnkTUbVwDjeRePjXRlSHLT66AqdyzqJlYFOM6fBild/POdxE4uBlc6I6qiRxA8CpnLNYfHSFyBERUWUxeRPVUSWJu7wyEbkuJm8iKlN8swEOy0QkHo55E1GZukR2AgBsOL8VTzTuYysTkfiYvImoXF0iOzFpE7kgXjYncmNGkxnJV3NhNJnFDoWIahB73kRuymgyY/rqQ0jLMiAs2BsJwzpzLjaRm2DPm8hNpWjzkJZlAACkZRmQos0TOSIiqilM3kRuihuHELkvJm8iN8WNQ4jcF5M3kZsKD1FDHvIPPNtvgzzkH4SHqMUOiYhqCJM3kZvadHY7lI3/gcLTCmXjf7Dp7HaxQyKiGsLkTeSm/szZDZms+LFMVlwmIvfA5E1UR3UMaeOwTESui8mbqI4a3naoLWF3DGmD4W2HihwREVUWV2wgqsOGtx2K4WIHQURVxp43ERGRxDB5ExERSQyTN5GEGU1mnEnN4sYjRHUMx7yJJIobjxDVXYL1vC0WC6ZOnYqBAwciPj4eqampdvX//e9/8fTTT6Nfv37Yvp2LRxBVFTceIaq7BPuZvmPHDphMJqxduxaJiYmYM2cOPvvsMwBAXl4e1qxZg19++QUGgwFPPvkkevXqJVQoREREbkWwnveRI0cQFxcHAGjXrh1Onjxpq/P29kZERAQMBgMMBgNkJctAEVGlhYeoIb/xpyOXgWuXE9UhgvW88/Pz4evrayt7eHjAbDZDoSj+yPDwcDzyyCMoKirCyJEjKzxeUJAPFArX2xVJo/ETOwTJYxtWjzbXCIu1+LHFChiKrGhSQVuyrcvHtnEe29B5lW1DwZK3r68vdDqdrWyxWGyJe/fu3bh27Rp+/fVXAMALL7yADh06oE2b8pdnzM7WCxVqtWk0fkhPvy52GJLGNqy+a+nXgej98AzJRWFmAK6lt0Z6gMrhe9jWZeP30HlsQ+fd3oaOErlgl807dOiA3buLN0JITExE06ZNbXUBAQFQqVRQKpXw8vKCn58f8vJ4sw1RVWy68g1U9XKhUACqernYdOUbsUMioloiWM+7V69e2LdvHwYNGgSr1YrZs2dj1apViIqKQs+ePfHHH39gwIABkMvl6NChA+69916hQiFyS9qiS3a7hmmLLokbEBHVGsGSt1wux3vvvWf3XGxsrO3xuHHjMG7cOKE+nqjOG9lqGJadXG1XJiL3wBXWiNxUm3otMbLVMIT7hGFkq2FoU6+l2CERUQ3hckxEbqxNvZZM2kRuiD1vIiIiiWHyJiIikhgmbyIiIolh8iYiIpIYJm8iIiKJYfImIiKSGCZvIiIiiWHyJnJhRpMZyVdzYTSZxQ6FiFwIF2khclFGkxnTVx9CWpYBYcHeSBjWGSol/2SJiD1vIpeVos1DWpYBAJCWZUCKljvvEVExJm8iF2UqLHJYJqK6i8mbyEUpPT0clomo7mLyJnJR4SFqyIMuwLPtdsiDLiA8RC12SETkIpi8iVzUtn/2QHlHEhReRVDekYRt/+wROyQichFM3kQuanfWL5DJih/LZMVlIiKAyZtIsloHNXdYJiL3xeRNJFGj2g+3JezWQc0xqv1wkSMiotrCFR+IJGxU++HQaPyQnn5d7FCIqBax501ERCQxTN5EREQSw+RNJCJuPEJE1cExbyKRcOMRIqou9ryJRMKNR4ioupi8iYiIJIbJm0gk4SFqyG+soCaXgWuXE1GlMXkTiUSbqYPFWvzYYi0uExFVBpM3kUi4XzcRVReTN5FIuF83EVUXkzeRSDjmTUTVxeRNJBKOeRNRdTF5E4mEY95EVF1M3kQi4Zg3EVUX12IkEkl4iBrymN/gEWREUbYK4SH3ih0SEUkEe95EIpl/8EMoQ41QKABlqBHzD34odkhEJBFM3kQiyUEWZDfuNpfJisu3qu8V6rBMRHUXkzeRi5py70Rbwq7vFYop904UOSIichUc8yZyYUzYRFQW9ryJiIgkhsmbSEBGkxnJV3NhNJnFDoWI3AgvmxMJxGgyY8YXh6HN1CM8xAdTnusElZJ/ckTkPPa8iQRyJUMHbaYeAKDN1ONKBpc/JaKaweRNJJAQfxU8buw84iGXIcRfJXJEROQumLyJBKLN1KHoxs4jRRYrNx4hohrD5E0kEG48QkRCEezuGYvFgmnTpuHMmTNQKpWYOXMmoqOjbfW///47Fi9eDABo2bIlEhISICtZborIDXDjESISimA97x07dsBkMmHt2rWYMGEC5syZY6vLz8/H/PnzsXTpUqxbtw6RkZHIzs4WKhQiUYSHqHFjyBtyWXGZiKgmCJa8jxw5gri4OABAu3btcPLkSVvdsWPH0LRpU8ydOxeDBw9GaGgogoODhQqFSBTaTB1uDHnDYgXHvImoxgh22Tw/Px++vr62soeHB8xmMxQKBbKzs3HgwAGsX78ePj4+GDJkCNq1a4eYmJhyjxcU5AOFwvUuO2o0fmKHIHnu2oYXbpsa5q32qvBcq9sW7tqGtYlt6Dy2ofMq24aCJW9fX1/odDf/8bJYLFAoij8uMDAQrVu3hkajAQB06tQJp0+fdpi8s7P1QoVabRqNH9LTr4sdhqS5cxsadAWlyhWda3Xawp3bsLawDZ3HNnTe7W3oKJELdtm8Q4cO2L17NwAgMTERTZs2tdW1atUKZ8+eRVZWFsxmM44fP4477rhDqFCIRMExbyISimA97169emHfvn0YNGgQrFYrZs+ejVWrViEqKgo9e/bEhAkT8OKLLwIA+vTpY5fcidxBWWPegb5e4gZFRG5BsOQtl8vx3nvv2T0XGxtre/zII4/gkUceEerjiUTHed5EJBQu0kIkEM7zJiKhMHkTCYRj3kQkFCZvIoFwnjcRCYXJm0ggHPMmIqEweRNVk9FkRvLVXBhN5jLrOeZNREIR7G5zIndmNJkxffUhpGUZEBbsjYRhnaFS2v85xYT7IyzY2/aamHB/kaIlInfDnjdRNaRo85CWZQAApGUZkKLNK/UalVKBBh1Pw/uurWjQ8XSp5E5EVF1M3kTVUJnx7KXHVuJUbhIA4FRuEpYeW1krsRGR+2PyJqqGyoxnn8hOclgmIqouJm+iaqiJOdyDmjzlsExEVB6HyXvHjh21FQeRpNTEHO64hl0xqMlT8PHwwaAmTyGuYdcajpKI3JXD5L148WLb4/HjxwseDJFU1NQc7riGXTG/xzQmbiKqEofJ22q12h6npKQIHgyRVHAONxGJyWHylslktRUHkaRw3XIiEpPDiaeFhYXQarWwWCy2x7f2xiMiIgQPkMgVca9uIhKTw+St1+sxdOhQW8IeMmSIrU4mk+HXX38VNjoiF8V1y4lITA6T986dO2srDiJJ4Zg3EYmpwvUac3JysHnzZpw/fx5eXl5o0qQJHn74YXh7e9dGfEQuqWTM22LlmDcR1T6HN6z9/fff6NOnD3bv3g1vb294eHjg559/Rq9evXDmzJnaipHI5XCvbiISk8Oe9wcffIC5c+eiR48eds/v3LkTc+fOxcqVXKuZ6iaOeRORmBz2vP/9999SiRsAHnjgAWRlZQkWFJErcLRfN8e8iUhMDnveSqWy3DrOASd3VtF+3dyrm4jEVKl53rfO7b61jshdlbVfd4voYFu9SqlAwrDOuJKhQ2Somnt1E1GtqtI871ux503urDJj2iqlArERAbUVEhGRDed5E5WBY9pE5Moq3M97//79SE5OtpXXrFmD/fv3CxoUkdi4djkRuTKHyfunn37C1KlTYTAYbM+FhIRgypQp2LZtm+DBEYmF87iJyJU5TN4rVqzAl19+iVatWtme69u3L1avXo1ly5YJHhyRWDiPm4hcWYX7eYeFhZV6vkGDBrBYLIIFRSQ2jnkTkSurMHnrdKUvF+bn53OqGLk1jnkTkStzmLyfeOIJvP7667h06ZLtuX///RdvvPEG+vTpI3hwRGLhmDcRuTKHU8Wef/55ZGdn47HHHoOnpyeUSiUMBgOGDh2KsWPH1laMRLWOY95E5MocJu//+7//w/jx4zFq1CicOnUKPj4+iI2NhZeXF2bNmoV33nmntuIkqlUc8yYiV+bwsvm3334LAPDx8cGsWbPQsmVLeHl5AQAOHz4sfHREIuGYNxG5sgpvWCvrMZG745g3EbmyCldYK8G1zKku4Zg3Ebkyh8mbCZvqKo55E5Erc3jD2rlz59CzZ08AQFpamu2x1WpFenq68NERCchoMpe7pWd4iBpyOWCxAHI5x7yJyLU4TN5cv5zcldFkxvTVh5CWZUBYsDcShnW2S+CZeUaULCJosRSXA3297I6x9NhKnMhOQuug5hjVfnhthk9EdZzD5B0ZGVlbcRDVqhRtHtKyijfcScsyIEWbhxbRwbb6yFA1wkN8oM3UIzzEB5Gh9j3vksQNACeyk7D02EomcCKqNQ6TN1FdpVIqMOW5TuVeVi9J3OWViYiEVOm7zYncSWXmcauUCsRGBJRK3JW8ZGshAAAZtElEQVQ1qMlTDstERNXFnjfVSWXN4759TNtZcQ27AgA2nt+Gxxv3tpWJiJzF5E11Um3N445r2JVJm4hqHC+bExERSQyTN9VJXISFiKRMsORtsVgwdepUDBw4EPHx8UhNTS3zNS+++CK+/vprocIgKhM3HiEiKRMsee/YsQMmkwlr167FhAkTMGfOnFKv+eijj5CbmytUCETl4sYjRCRlgiXvI0eOIC4uDgDQrl07nDx50q5+69atkMlk6N69u1AhEJWLG48QkZQJdrd5fn4+fH19bWUPDw+YzWYoFAqcPXsWmzdvxieffILFixdX6nhBQT5QKFxvXFKj8RM7BMkTow3r5Rjty6F+Tsch5neB30PnsQ2dxzZ0XmXbULDk7evrC53u5qVIi8UChaL449avX4+0tDQ899xzuHLlCjw9PREZGemwF56drRcq1GrTaPyQnn5d7DAkTaw29FbIIJcVXzKXy4rLzsYh1neB30PnsQ2dxzZ03u1t6CiRC5a8O3TogF27dqFv375ITExE06ZNbXUTJ060PV60aBFCQ0N5+ZxqnKNdw2pjkRYiIqEIlrx79eqFffv2YdCgQbBarZg9ezZWrVqFqKgo29aiREIxmsyY8cVh28YiU57rVO1lTomIXI1g/5rJ5XK89957ds/FxsaWet0rr7wiVAhUh13J0EGbWTzUos3U40qGDrERAbb6mHB/hAV727YEjQn3FytUIqIqY1eE3FKIvwoechmKLFZ4yGUI8VfZ1auUCiQM61zuZXUiIlfGFdbILWkzdSi6MahdZLGWOY/b2V3DiIjEwuRNbonzuInInTF5k1vi2uVE5M6YvMktce1yInJnTN7klrh2ORG5MyZvcksc8yYid8bkTW6JY95E5M6YvMktccybiNwZkze5JY55E5E7Y/Imt8QxbyJyZ0ze5JY45k1E7ozJm9wSx7yJyJ0xeZNb4pg3EbkzJm9ySxzzJiJ3xuRNkmU0mZF8NRdGk7lUHce8icidcS9EkiSjyYzpqw8hLcuAsGBvJAzrbLe1Z0y4P8KCvW31MeH+IkZLRFSzmLxJklK0eUjLMgAA0rIMSNHmoUV0sK1epVQgYVhnXMnQITJUzT27icit8F80kqTKjGmrlArERgTUVkhERLWGY94kSRzTJqK6jMmbJInzuImoLmPyJkniPG4iqsuYvEmSOI+biOoyJm+SpNoY884tyMMfVw8ityCvxo9NROQM3m1OklQy5m2xCjPmnVuQh7f3zbSVZ9/7LgK8OFeciFwDe94kSUKPec/580OHZSIiMTF5kyQJPeadV6RzWG7gXd9hmYhISEzeJEliz/Oe3HW8LWE38K6PyV3H1+rnE1HdxjFvkiShx7wrgwmbiMTCnje5LEe7hnGeNxHVZex5k0symsyY8cVhaDP1CA/xwZTnOnFzESKiG9jzJpd0JUMHbaYeAKDN1ONKhn3PumTLTwDc8pOI6hx2Zcglhfir4CGXochihYdchhB/lV09t/wkorqMPW9ySdpMHYpuDGoXWaxljmmXbPnJxE1EdQ2TN7kkrl1ORFQ+Jm9ySWLP4yYicmVM3uSSuF83EVH5mLzJJXEeNxFR+Zi8ySVxzJuIqHxM3uSSOOZNRFQ+Jm9ySRzzJiIqH5M3uSSOeRMRlY/Jm1wSx7yJiMrH5E2iMZrMOJOaVeauYRzzJiIqH9eVJFEYTWZMX30IaVkGhAV7I2FYZ7tlTks2Himp58YjREQ3MXmTKFK0eUjLMgAA0rIMSNHmoUV0sK2eG48QEZWP/yKSKCozpl2y8QgREdkTLHlbLBZMmzYNZ86cgVKpxMyZMxEdHW2rX716NbZs2QIA6NGjB8aOHStUKOSCOKZNRFR9gt2wtmPHDphMJqxduxYTJkzAnDlzbHWXLl3Cxo0b8c0332Dt2rXYu3cvkpKShAqFXBDncRMRVZ9gyfvIkSOIi4sDALRr1w4nT5601dWvXx8rVqyAh4cH5HI5zGYzvLy8hAqFXBDncRMRVZ9gl83z8/Ph6+trK3t4eMBsNkOhUMDT0xPBwcGwWq2YN28eWrZsiZiYGIfHCwrygULhepdWNRo/sUOQpAsZ9snaW+3l8m3pyvG5cmxSwTZ0HtvQeZVtQ8GSt6+vL3S6m/9AWywWKBQ3P66goABvv/021Go1EhISKjxedrZekDidodH4IT39uthhSFJujr5U2dXb0lXj4/fQeWxD57ENnXd7GzpK5IJdNu/QoQN2794NAEhMTETTpk1tdVarFS+//DKaNWuG9957Dx4ertejJmHxhjUiouoTrOfdq1cv7Nu3D4MGDYLVasXs2bOxatUqREVFwWKx4ODBgzCZTNizZw8AYPz48Wjfvr1Q4ZAIjCZzufO0w0PUkMsBiwWQy3nDGhFRVQiWvOVyOd577z2752JjY22PT5w4IdRHkwuoaAW1zDwjLJbixxZLcTnQlzctEhFVBtc2J0GUtYLarSJD1QgP8QEAhIf4IDKUPW8iosriCmskCpVSgSnPdYLebIWPQsblT4mIqoA9bxJEZRZhUSkVaBYdzMRNRFRFTN4kCC7CQkQkHCZvEkRlNh4hIqLqYfImQXAeNxGRcJi8SRDceISISDhM3iQIjnkTEQmHyZsEwTFvIiLhMHmTIDjmTUQkHCZvqjajyYzkq7kwmsyl6krWLge4djkRUU3j6hhULVy7nIhIPOx5U7W4w9rlPyVvx7idk/FT8nYYzQVIyb0Io7lA7LCIiCrEnjcJomTt8vK2BBXbT8nbsSV1OwBgS+p2/HblD+jMOoT51MPETq+IHB0RkWOu9S8qSUbJPG6L1fHa5bERASJEV7GSxF1CZy6eypamvwatLk2MkIiIKo2Xzala3HUed5hPPYSrw9BV08nu+dvLRERiYs+bqsUd53G/0XEswtVhUCm8MLT1AOAEsD/9MLpqOhWXiYhcBJM3VYs7zuOOCYiyKw9tPQBDwaRNRK6Hl82pWqSwdnlZd5DnFuThj6sHRYyKiMh57HlTtZQ15u1K87iN5gJM3DMNRdYieMg8MC9uGgqKCjDlj/dRZJX+JX4iqtvY86ZqcfUx71n7F9iSdJG1CLP2L8DfmUlM3ETkFpi8qUyOlj4FXH/MO6swt1T5zpDm8JCVHWeUT0RthEVEVCN42ZxKqWjpUwCICfdHWLC37TUx4f4iRVt5AV7+mHHPZPydmYQ7Q5pj6bGVuKi/iiifCEzq8prY4RERVRqTN5VS1tKnLaKD7V6jUiqQMKyzy66gVp4AL3/cE3EXADBhE5FkSeNfXKpVlR3PFnsFtQ8PLcE/1y+ggU84Xu/0MlSK4hvmuD45Ebk7jnlTKa4+ng3cTNwAcFmvxVt734PRXACjuQDzDi8SNzgiIoGx502lVGbdcrGVJO4ShZZC25rkafprpV4vh6w2wiIiqhXseVMpUl23PFwdhnB1GMJ86pWqm9T5VREiIiISBnveVIqrzOFefHQFTuWcRcvAphjT4cUKX18y5j2x0yvQ6tJgtVqw9+qfeKBhdzTw41QwInIfTN5UiiuMeZckbgA4lXMWi4+usCXwim5IUym8bOuUNw5sJGicRERi4GVzKsUV1i0vSdxllbnfNhHVdex511FGk7ncOdq1tW75NX0Gtib+gnqK+mgZ0gwAcC77PNL06Q7fF64Og7/cD3mW67bnhrcYUuPxERG5KibvOshoMmPGF4ehzdQjPMQHU57rVOuLrFzTZ2D6n/Ns5RCvIAAyZBZkVfhelcILCd0mYtfFPfhDewhPNu6LjuFtBYyWiMi1MHnXQVcydNBm6gEA2kw9rmTo7BZbqY2lT+cf/MSunFmQXaX3qxReeLjxg3i48YM1GRYRkSQweddBIf4qeMhlKLJY4SGXIcRfZVdfU0ufJmWew+aUbXg0pjeahzSB0VyA1LxLAAC9xVg6Lq/gcnvezQPuqFYMRETuiMm7DtJm6lB0Y1C7yGItc0zb2aVPkzLPYdHxzwEAi45/jpGthuGH5M1IN2SU+563737dNuZdX63Brxd+x9m882gecAde6fhStWMhInI3TN51UG3M415yfKVd+fOT/4MFFofvUSm80FrTAq3RAgDQKrRFjcdFROQOmLzroJqax73oyHIk5f6DSO/6uCuiE+qrNWjoF4l/dddQBPsfBBZYoPEOLbfnzcviRESVx+TtphxNBQsPUUMuBywWQC6v3jzuksQNAFcM/+LH5M0AADnk5faw3+r8qm3MO9q/IT4//gWScv/hZXEioipi8nZDRpMZ01cfst0tnjCss10Cz8wzwnIjv1osxeXbx7xLknNT/8boc+OO7mj/higoKsDfmUm2xH278hJ3kDwQKoUXmgXf7GG/0vElaDR+SE+/XuZ7iIiobEzebihFm4e0LAMAIC3LgBRtHlpEB9vqI0PVCA/xsc3zjgxVw2gugFaXBrWnDz49sgKZhcV3fZ/NO4+zicsBFN8NnmPKRZG1/DHy8nre73abUJOnSERUpzF510EqpQKhbRORk38eZg81rurr44tT3yDD6HiBlIoWUBnd5nnbmLfOpEdWQQ7qqzW4I7CxbdMQIiJyHpO3RFV1TPvWOdY/Jf+Cf/IvAADyi3RYeHRJpT7TUc97cufXbDt3BXjV/KIuRER0E5O3BFVqTDt6LzyD8lGY7YvT1xrg+0tfwlBUemGUygjzCsXAFk/bjXnfGdIc10352HlpN7fcJCKqZUzeEpSizUOaxxF4driItEtRSNE2sxvT/i71C3iH5gMyQBGajy8vrKjS8R+OfhAHrxxFpjkLjdQN8Obd42x1KoUX7om4C0BxD/vZloNq5qSIiKjSmLxd1C/nf8dPF3agb6MH8VDjHnZ1e679ClWji5DJAI9GF7Hn2q9oEf0fW/1Fw2Xgxpaetv86EOoZjDEdX8Sf2kPoEt4Z9XxC8WjsQzV4NkREVJME28/bYrFg6tSpGDhwIOLj45GammpXv27dOjz99NMYMGAAdu3aJVQY1WI0mZF8NRdGk7nax7h8/Sr+d+obXL5+tcz6r//+HmN2TsTXf39fqu6X879jQ8oWFKIAG1K24Jfzv9vVH79+CLIbSVkmKy5Xlr+nH2bf+y4aeNcHADTwro/pcW+hnk8oHo99GPV8Qit9LCIiEodgPe8dO3bAZDJh7dq1SExMxJw5c/DZZ58BANLT07FmzRp8//33KCgowODBg3HvvfdCqVQKFY6dkmlR4eqwUndBG01mTPvyN2R6nkNIYRNMG3pfmRtznM+5gBWnf8cD4T3QOLCRXd3l61fx/qGPAAAH/j1qdzMXUJy496YdAADbf5+5s5+tftOFn+16zpsu/Fyq9+1It7C7bccFgI4hbRDoE4BG/tFoGdIMKoUXJncdX+njERGRaxGs533kyBHExcUBANq1a4eTJ0/a6v766y+0b98eSqUSfn5+iIqKQlJSklCh2DGaCzDv8CIsOPIp5h1eBKO5wK7+5NVLyIv+BcrIFORF/4KTVy+VOsb5nAtYeHQJjmn/xsKjS3A+54Jd/f/+XuuwfGtiLat8+zzp28vtgu90WH7mzn7oFnY3gOJEPrztUDzd5DF0CGvDKVtERG5AsJ53fn4+fH19bWUPDw+YzWYoFArk5+fDz8/PVqdWq5Gfn+/weEFBPlAoqrcG963OZWYgTX8NAJCmvwaj8joahty8VJx+KdnuknS6RzI0mjZ2x1hx2v4y9k7t77i7SWtb+apea1d/Va+FRnPzfHvFdMP2lL125Vvrh7Xtj9XHv7Mr31r/dq+xWLh7GQ5oE3F3eDtM6D6y1HmOu28YxmFYue3gSm49N6oetqHz2IbOYxs6r7JtKFjy9vX1hU6ns5UtFgsUCkWZdTqdzi6ZlyU7W18jcanMfgjzqYc0/TWE+dSDymS/PGfHkHbYfO4Xu/Lty3c+EN4Dx7R/25VvfU3f6F7Ykrrdrnxr/ZMxj8OgL8TetAPoFnY3nox53K6+c8hduB5bgM0pv+DRmIfQOeSuUjE822Iwnm0xGAAkvbwol0d1HtvQeWxD57ENnXd7GzpK5IIl7w4dOmDXrl3o27cvEhMT0bRpU1tdmzZt8NFHH6GgoAAmkwnJycl29UJSKbwwsdMr5Y551/MJRUKXiXZ3Xt+ucWAjTOjwMnZqyx7z7hvbCwCwNXUn+kQ/YCvf6pk7+9mNc9/ugeg4PBAdV40zJCIidyezWq1WIQ5ssVgwbdo0nD17FlarFbNnz8bu3bsRFRWFnj17Yt26dVi7di2sVitGjhyJ3r17OzyeK/6i4y9N57ENncc2dB7b0HlsQ+dVpectWPKuaa74peCX1XlsQ+exDZ3HNnQe29B5VUnegt1tTkRERMJg8iYiIpIYJm8iIiKJYfImIiKSGCZvIiIiiWHyJiIikhgmbyIiIolh8iYiIpIYJm8iIiKJkcwKa0RERFSMPW8iIiKJYfImIiKSGCZvIiIiiWHyJiIikhgmbyIiIolh8iYiIpIYJu9qOn78OOLj48UOQ5IKCwvx5ptvYvDgwejfvz9+/fVXsUOSpKKiIkyePBmDBg3CkCFDcPHiRbFDkqTMzEz06NEDycnJYociWU8++STi4+MRHx+PyZMnix2OJC1btgwDBw7E008/jW+//bbC1ytqISa38/nnn2Pjxo3w9vYWOxRJ2rhxIwIDAzF//nxkZ2fjqaeeQs+ePcUOS3J27doFAPjmm29w4MABvP/++/jss89EjkpaCgsLMXXqVKhUKrFDkayCggIAwJo1a0SORLoOHDiAY8eO4euvv4bBYMDKlSsrfA973tUQFRWFRYsWiR2GZPXp0wevvvqqrezh4SFiNNL14IMPYsaMGQCAq1evIjQ0VOSIpGfu3LkYNGgQ6tWrJ3YokpWUlASDwYDhw4fj2WefRWJiotghSc7evXvRtGlTjBkzBqNGjcJ9991X4XvY866G3r174/Lly2KHIVlqtRoAkJ+fj3HjxuG1114TOSLpUigUmDRpErZv345PPvlE7HAk5YcffkBwcDDi4uKwfPlyscORLJVKhRdeeAH/+c9/cOHCBYwYMQJbt26FQsH0UlnZ2dm4evUqli5disuXL2P06NHYunUrZDJZue9hz5tEodVq8eyzz+KJJ57AY489JnY4kjZ37lxs27YNU6ZMgV6vFzscyfj+++/xxx9/ID4+HqdPn8akSZOQnp4udliSExMTg8cffxwymQwxMTEIDAxkO1ZRYGAgunXrBqVSicaNG8PLywtZWVkO38PkTbUuIyMDw4cPx5tvvon+/fuLHY5krV+/HsuWLQMAeHt7QyaTcQiiCr766it8+eWXWLNmDVq0aIG5c+dCo9GIHZbkfPfdd5gzZw4AIC0tDfn5+WzHKurYsSP27NkDq9WKtLQ0GAwGBAYGOnwPr2tQrVu6dCny8vKwZMkSLFmyBEDxTYC8aahqHnroIUyePBlDhgyB2WzG22+/DS8vL7HDojqmf//+mDx5Mp555hnIZDLMnj2bl8yr6P7778ehQ4fQv39/WK1WTJ06tcIf4txVjIiISGJ42ZyIiEhimLyJiIgkhsmbiIhIYpi8iYiIJIbJm4iISGKYvImoUj7++GPbJjJV3ZTnwIED3MiHqAZxMh4RVcqt69EfPHhQxEiIiMmbSMIOHDiApUuXwtPTE5cvX8YDDzwAHx8f7NixAwCwfPlybN26FRs2bIDBYICnpycWLlyIxo0b48CBA5g5cyY8PDzQrl07JCcnY82aNYiPj0fr1q1x5MgRZGVl4d1330WPHj3w1ltv4a677sKpU6cAAP/5z3/w7bffolmzZjhz5gyA4vXCDx48iDlz5mDv3r14//334eXlhZiYGFvMqampmDZtGnJycqBSqTBlyhS0bNmy9huPSMJ42ZxI4o4fP47p06fj+++/x1dffYXg4GD88MMPaNasGbZs2YIdO3ZgzZo12Lx5M+677z589dVXKCwsxMSJEzF//nysX7++1IpYhYWFWLt2LSZPnoyPP/7Yru7dd98FAId7DptMJrz11lv45JNP8MMPP9itnjdp0iS8+eab+PHHHzFjxgy8/vrrNdgaRHUDe95EEte0aVOEh4cDAIKCgtC1a1cAQEREBPLy8rBw4UJs2bIFFy5cwJ49e9CiRQucPXsWISEhaN68OYDiJS5nzZplO2ZcXBwAoEmTJsjJyalyTGfOnEG9evUQGxsLAHjqqafw8ccfQ6fT4eTJk5g8ebLttXq9HtnZ2QgKCqpeAxDVQUzeRBLn6elpV751TWStVouBAwdi6NCh6N69O0JDQ3H69Gl4eHjAYrGUe8ySNdIdbUl4K6vVCplMBrPZbHvfrSsvl8RksVigVCqxYcMGW92///5b4SYMRGSPl82J3NiJEycQHR2NYcOGoXXr1tixYweKiorQuHFj5OXl2caqN23aVKXjenh42BJ1UFAQzp07B6vVip07dwIAmjVrhoyMDCQlJQEAtmzZAgDw8/NDo0aNbMl73759GDJkSI2cK1Fdwp43kRvr1q0bkpKS0LdvX1itVnTu3Bnnzp2DUqnEvHnzMGnSJMjlcsTExFRpV7eePXviiSeewA8//IAJEyZg1KhRCA0NRceOHZGdnQ1PT0988MEHePPNN6FQKOxuSJs/fz6mTZuGFStWwNPTEx9++GGle/hEVIy7ihHVQRaLBQsWLMDYsWPh4+ODVatWIS0tDW+99ZbYoRFRJbDnTVQHyeVyBAYGon///vD09ERkZKTdDWtE5NrY8yYiIpIY3rBGREQkMUzeREREEsPkTUREJDFM3kRERBLD5E1ERCQxTN5EREQS8/+MQcqyXCTvIgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x396 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Slice the magnitudes for before and after 2010: mags_pre, mags_post\n",
    "mags_pre = mags[time < 2010]\n",
    "mags_post = mags[time >= 2010]\n",
    "\n",
    "# Generate the ECDFs\n",
    "_ = plt.plot(*dcst.ecdf(mags_pre), marker='.', linestyle='none')\n",
    "_ = plt.plot(*dcst.ecdf(mags_post), marker='.', linestyle='none')\n",
    "\n",
    "# Label the axes and show the plot\n",
    "_ = plt.xlabel('magnitude')\n",
    "_ = plt.ylabel('ECDF')\n",
    "plt.legend(('1980 through 2009', '2010 through 2017'), loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both curves appear to follow the Gutenberg-Richter law but with different completeness thresholds. The difference could be attributed to improvements in sensing capabilities in more recent years.\n",
    "\n",
    "We can safely use a completeness threshold of mt=3 and compute the b-values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function that computes the b-value: b_value\n",
    "def b_value(mags, mt, perc, n_reps):\n",
    "    \"\"\"Compute the b-value and optionally its confidence interval\"\"\"\n",
    "    \n",
    "    # Extract the magnitudes above the completeness threshold: m\n",
    "    m = mags[mags >= mt]\n",
    "    \n",
    "    # Compute the b-value: b\n",
    "    b = np.mean((m) - mt) * np.log(10)\n",
    "    \n",
    "    # Draw the bootstrap replicates: m_bs_reps\n",
    "    if n_reps is None:\n",
    "        return b\n",
    "    else:\n",
    "        m_bs_reps = dcst.draw_bs_reps(m, np.mean, size=n_reps)\n",
    "        \n",
    "        # Compute the b-value from the replicates: b_bs_reps\n",
    "        b_bs_reps = (m_bs_reps - mt) * np.log(10)\n",
    "        \n",
    "        # Compute the confidence interval: conf_int\n",
    "        conf_int = np.percentile(b_bs_reps, perc)\n",
    "        \n",
    "        return b, conf_int"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "1980 through 2009\n",
      "b-value: 0.74\n",
      "95% conf int: [0.54, 0.97]\n",
      "\n",
      "2010 through mid-2017\n",
      "b-value: 0.62\n",
      "95% conf_int: [0.60, 0.65]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "mt=3\n",
    "# Use the b_value() function to compute b-value and the\n",
    "# Confidence interval for pre-2010, b_pre, conf_int_pre\n",
    "b_pre, conf_int_pre = b_value(mags_pre, mt, perc=[2.5, 97.5], n_reps=10000)\n",
    "\n",
    "# Compute the b-value and confidence interval for post-2010: b_post, conf_int_post\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 eartherquakes. The confidence intervals, however, do overlap.\n",
    "\n",
    "## Hypothesis testing\n",
    "We can perform a hypothesis test to see if we might get these results if the b-values were actually the same.\n",
    "\n",
    "For the hypothesis test, we should only include earthquakes with a magnitude $\\geq$ the completeness threshold. below the threshold, we are not comparing earthquakes before and after 2010. We want to compare <i>observed</i> earthquake and we do not have a complete dataset below the completeness threshold.\n",
    "\n",
    "A permutation test is appropriate because our hypothesis assumes both sets of magnitueds are Exponentially distributed and we are only using data for earthquakes above the completeness theshold. We are also using the G-R Law (our model) to deal with missing data.\n",
    "\n",
    "Test statistic: the difference between the mean post-2010 magnitude and the pre-2010 magnitude is a reasonable test stat because the mean value is directly proportional to the difference in b-value.\n",
    "\n",
    "We do not need to worry about the difference in the number of earthquakes before and after 2010 because we are only interested in summary statistics about their magnitudes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p= 0.1072\n"
     ]
    }
   ],
   "source": [
    "# Slice out the magnitudes above the completeness threshold: mags_pre, mags_post\n",
    "mags_pre = mags_pre[mags_pre >= mt]\n",
    "mags_post = mags_post[mags_post >= mt]\n",
    "\n",
    "# Compute the observed difference in mean magnitudes: diff_obs\n",
    "diff_obs = np.mean(mags_post) - np.mean(mags_pre)\n",
    "\n",
    "# Generate 10,000 permutation replicateds, use dcst.diff_of_means, as the kwarg: perm_reps\n",
    "perm_reps = dcst.draw_perm_reps(mags_post, mags_pre, dcst.diff_of_means, size=10000)\n",
    "\n",
    "# Compute and print the p-value: p_val\n",
    "p_val = np.sum(perm_reps <= diff_obs) / 10000\n",
    "print('p=', p_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A p-value that close to 0.1 suggests that the observed magnitudes are indicitive of there being no change in b-value after waste-water injection began.\n",
    "\n",
    "## Conclusions from the analysis of earthquakes:\n",
    "From 2010 to mid-2017 Oklahoma had twice as many magnitude 3+ earthquakes than the entire state of California.\n",
    "\n",
    "The seismicity, as measured by the b-value, is comparable before and after water injection. The earthquakes became moe frequent after fracking but the magnitudes did not change much. Earthquakes are over 100 times more frequent in Oklahoma after widespread waste-water injection began.\n",
    "\n",
    "Oklahoma has a b-value smaller than the <a href=\"https://github.com/btab2273/statistical_seismology/blob/master/statistical_seismology_parkfield_ca.ipynb\">Parkfield region</a>. This means a randomly selected earthquake  above magnitude 3 in Oklahoma more likely than not has a smaller magnitude than an earthquake above magnitude 3 randomly selected from the Parkfield region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save the arrays as .csv files\n",
    "times_doc = pd.DataFrame(time)\n",
    "mags_doc = pd.DataFrame(mags)\n",
    "mags_arr = np.array(mags_doc)\n",
    "times_doc.to_csv('oklahoma_decimal_dates.csv')\n",
    "mags_doc.to_csv('oklahoma_magnitudes.csv')"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}