cossatot/evenly_truncated_mfd_debug.ipynb

## evenly_truncated_mfd_debug.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Diagnosing rate problems with evenly-distributed MFDs\n",
    "\n",
    "\n",
    "Based on manual inspection and quantitative testing of the XML fault source files produced through OpenQuake code, we have found that in many cases the bin rates for the MFDs are too high.  This leads to very high total annual moment release rates that are above the moment accumulation rates (determined by the slip rate, fault area and rigidity).\n",
    "\n",
    "This notebook focuses on one aspect of the problem that is probably dominant over the other possibilities (inaccurate fault area calculations, etc.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "import openquake.mbt as oqmbt\n",
    "from openquake.mbt.tools.faults import (rates_for_double_truncated_mfd, \n",
    "                                        _get_rate_above_m_low, _get_cumul_rate_truncated)\n",
    "from openquake.mbt.tools.mfd import EvenlyDiscretizedMFD, mo_to_mag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The major function that does these calculations is called `rates_for_double_truncated_mfd` and it is in the `oq-mbtk/openquake/mbt/tools/faults.py` source file.  It calls two helper functions, `_get_rate_above_m_low` and `_get_cumul_rate_truncated`.  These functions are copied below, so the source code is visible:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```python\n",
    "def _get_rate_above_m_low(seismic_moment, m_low, m_upp, b_gr):\n",
    "    \"\"\"\n",
    "    :parameter seismic_moment:\n",
    "        Seismic moment in Nm\n",
    "    :parameter m_low:\n",
    "        Lower magnitude threshold\n",
    "    :parameter m_upp:\n",
    "        Upper magnitude threshold\n",
    "    :parameter b_gr:\n",
    "        b value of the Gutenberg-Richter relationship\n",
    "    \"\"\"\n",
    "    a_m = 9.1\n",
    "    b_m = 1.5\n",
    "    beta = b_gr * numpy.log(10.)\n",
    "    x = (-seismic_moment*(b_m*numpy.log(10.) - beta) /\n",
    "         (beta*(10**(a_m + b_m*m_low) -\n",
    "          10**(a_m + b_m*m_upp)*numpy.exp(beta*(m_low - m_upp)))))\n",
    "    rate_m_low = x * (1-numpy.exp(-beta*(m_upp-m_low)))\n",
    "    return rate_m_low\n",
    "\n",
    "\n",
    "def _get_cumul_rate_truncated(m, m_low, m_upp, rate_gt_m_low, b_gr):\n",
    "    \"\"\"\n",
    "    This is basically equation 9 of Youngs and Coppersmith (1985)\n",
    "    \"\"\"\n",
    "    beta = b_gr * numpy.log(10.)\n",
    "    nmr1 = numpy.exp(-beta*(m-m_low))\n",
    "    nmr2 = numpy.exp(-beta*(m_upp-m_low))\n",
    "    den1 = 1-numpy.exp(-beta*(m_upp-m_low))\n",
    "    rate = rate_gt_m_low * (nmr1 - nmr2) / den1\n",
    "    return rate\n",
    "\n",
    "\n",
    "def rates_for_double_truncated_mfd(area, slip_rate, m_low, m_upp, b_gr,\n",
    "                                   bin_width=0.1, rigidity=32e9):\n",
    "    \"\"\"\n",
    "    :parameter area:\n",
    "        Area of the fault surface\n",
    "        float [km2]\n",
    "    :parameter slip_rate:\n",
    "        Slip-rate\n",
    "        float [mm/tr]\n",
    "    :parameter m_low:\n",
    "        Lower magnitude\n",
    "        float\n",
    "    :parameter m_upp:\n",
    "        Upper magnitude\n",
    "    :parameter b_gr:\n",
    "        b-value of Gutenber-Richter relationship\n",
    "    :parameter bin_width:\n",
    "        Bin width\n",
    "    :parameter rigidity:\n",
    "        Rigidity [Pa]\n",
    "    :return:\n",
    "        A list containing the rates per bin starting from m_low\n",
    "    \"\"\"\n",
    "    #\n",
    "    # Compute moment\n",
    "    slip_m = slip_rate * 1e-3  # [mm/yr] -> [m/yr]\n",
    "    area_m2 = area * 1e6\n",
    "    moment_from_slip = (rigidity * area_m2 * slip_m)\n",
    "    #\n",
    "    # Compute total rate\n",
    "    rate_above = _get_rate_above_m_low(moment_from_slip, m_low, m_upp, b_gr)\n",
    "    #\n",
    "    # Compute rate per bin\n",
    "    rrr = []\n",
    "    mma = []\n",
    "    for mmm in numpy.arange(m_low, m_upp, bin_width):\n",
    "        rte = (_get_cumul_rate_truncated(mmm, m_low, m_upp, rate_above, b_gr) -\n",
    "               _get_cumul_rate_truncated(mmm+bin_width, m_low,\n",
    "                                         m_upp, rate_above, b_gr))\n",
    "        mma.append(mmm+bin_width/2)\n",
    "        rrr.append(rte)\n",
    "return rrr\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we need to do some comparisions. We will write a few helper functions to let us calculate the annual moment accumulation and release from the fault data that we have."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_moment_from_mfd(mfd, threshold=-1):\n",
    "    \"\"\"\n",
    "    This computes the total scalar seismic moment released per year by a\n",
    "    source\n",
    "    :parameter mfd:\n",
    "        An instance of openquake.hazardlib.mfd\n",
    "    :param threshold:\n",
    "        Lower threshold magnitude\n",
    "    :returns:\n",
    "        A float corresponding to the rate of scalar moment released\n",
    "    \"\"\"\n",
    "    \n",
    "    # Note that this functions uses 9.1 instead of 9.05 to calculate\n",
    "    # moment from magnitude, which is consistent with the rest of\n",
    "    # the functions in the mbt.tools.mfd library and with the\n",
    "    # calculation of moment from fault area used here.\n",
    "    \n",
    "    if isinstance(mfd, (EvenlyDiscretizedMFD)):\n",
    "        occ_list = mfd.get_annual_occurrence_rates()\n",
    "        mo_tot = 0.0\n",
    "        for occ in occ_list:\n",
    "            if occ[0] > threshold:\n",
    "                mo_tot += occ[1] * 10.**(1.5*occ[0] + 9.1)\n",
    "    else:\n",
    "        raise ValueError('Unrecognised MFD type: %s' % type(mfd))\n",
    "    return mo_tot\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "def calc_moment_release_from_bin_rates(min_bin_midpt, bin_width, \n",
    "                                       occurrence_rates):\n",
    "    \n",
    "    EDMFD = EvenlyDiscretizedMFD(min_bin_midpt, bin_width, occurrence_rates)\n",
    "    \n",
    "    return get_moment_from_mfd(EDMFD)\n",
    "\n",
    "\n",
    "def get_bin_midpts(M_low, M_upp, bin_width):\n",
    "    \n",
    "    '''\n",
    "    Copied from `rates_for_double_truncated_mfd`\n",
    "    '''\n",
    "    mma = []\n",
    "    for mmm in numpy.arange(M_low, M_upp, bin_width):\n",
    "        mma.append(mmm + bin_width/2)\n",
    "    \n",
    "    return mma\n",
    "\n",
    "\n",
    "def calc_moment_release_from_fault_params(area, slip_rate, M_low, M_upp,\n",
    "                                          b_gr, bin_width=0.1, rigidity=32e9):\n",
    "    \n",
    "    bin_rates = rates_for_double_truncated_mfd(area, slip_rate, M_low, M_upp,\n",
    "                                              b_gr, bin_width, rigidity)\n",
    "    \n",
    "    bin_mids = get_bin_midpts(M_low, M_upp, bin_width)\n",
    "    first_bin_mid = bin_mids[0]\n",
    "    \n",
    "    return calc_moment_release_from_bin_rates(first_bin_mid, bin_width, bin_rates)\n",
    "\n",
    "\n",
    "def annual_mo_from_fault_params(area, slip_rate, rigidity=32e9):\n",
    "    slip_m = slip_rate * 1e-3  # [mm/yr] -> [m/yr]\n",
    "    area_m2 = area * 1e6\n",
    "    moment_from_slip = (rigidity * area_m2 * slip_m)\n",
    "    \n",
    "    return moment_from_slip\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Testing the sensitivity of annual moment release to `M_max`\n",
    "\n",
    "Let's do a test where we loop through `M_max` values, starting very close to `M_min` and moving up to over 3 orders of magnitude higher. The step size is much smaller than the `bin_width` to avoid aliasing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAI2CAYAAACBqkT0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3XmYXGWV+PHv6U4ngRDCqgOihF0CWQlhcSKgrCphExX0ZwIKg4AiqCCDCDqig8qijAyDjixREQVBEFTQMbK5kEBA9kUCRJAlQGcn2/n9Ubc7lU51dwWqUunu7+d56knXXc+9XQ2nT5/7vpGZSJIkSaq9pkYHIEmSJPVWJtuSJElSnZhsS5IkSXVisi1JkiTVicm2JEmSVCcm25IkSVKdmGxLkiRJdWKyLalPiYjLI+JrjY6jr4mIQyLi2YiYGxGjV+N5PxoRt6yu85Wd910R8XhxvQev7vNLWnOYbEt9SETMiIhFEbFRh+X3RkRGxNDGRFaZifHKImLPiJjZ6DjegG8DJ2bmOpl5bz1OEBFDi89xv7ZlmfnjzNy3HufrxleB/yqu9/qOK3vaz6KkN85kW+p7ngKOaHsTEcOBtRsXTu2UJ1k9UZT01v8ubw482OggVqNqrrfX/ixKWq63/kddUucmAx8vez8RuLJ8g4gYEhFXRsRLEfF0RHypLQmMiEkRcWdEXBARr0XE3yNi92L5sxHxYkRMLDvWgIj4dkQ8ExEvRMQlEbFWsW7PiJgZEZ8r9ns+Io4q1h0LfBQ4tfhT/I2VLqaoAp4QEY8DjxfL3hkRt0bEKxHxaER8qLObEREfiIjpxbXcFREjytZ9MSKejIg5EfFQRBxStm7riPhjRLRGxMsRcXXZulU5/5SIOCci7gTmA1tGxFER8XBx3r9HxL8V2w4Cfg1sWtyTuRGxaUQ0lcU6KyJ+FhEbdHK+9SPiV8X39tXi683K1k8qzjknIp6KiI92cpxxEfGn4r49HxH/FRH9K2w3ICLmAs3AfRHxZLE8I2Lrsu3a/4rR1eeiWL9WRJxXfDZbI+KO4jN1W7HJa8W92a24njvK9t09Iu4u9rs7Inbv8L34j+LzPScibokOlecO13ZMRDxRfJ9viIhNi+VPAlsCNxZxDOjkEN3+LHZx7ssj4uKI+HVxjjsj4l8i4sLi+/pIlLXrdPNZ/u+IuLbs/bkR8fuIiGpikdSNzPTly1cfeQEzgL2BR4HtKSVAMylV4RIYWmx3JfBLYDAwFHgM+ESxbhKwBDiq2P9rwDPA94ABwL7AHGCdYvsLgBuADYrj3Qh8o1i3Z3GsrwItwPsoJZzrF+svB77WzTUlcGtx/LWAQcCzRXz9gNHAy8Cwjscs1r0I7FJcy8TiHg0o1h8ObEqpMPFhYB6wSbHuKuCMYt1A4F+L5V2ev0L8U4r7t0OxfQvwfmArIIA9insypuyezexwjJOAPwObFd+D/wGu6uR8GwKHUaqgDgZ+DlxfFvtsYLvi/SbADp0cZydg1yLmocDDwGe7+T5t3cX78u9Ld5+L7xX37W3F92334rqHFsftV3bcScAdxdcbAK8C/6+I+4ji/YZl34sngW0pfZamAP/ZyfW8p/i+jinOfRFwW8eftTf7s9jF/pcX59+J0ufv/yhVyj/O8p/LP5Rt39VneW1KP+OTgPHFcTdr9H+vfPnqLS8r21Lf1FZR24dSkvSPthUR0Qx8BDg9M+dk5gzgPEoJSpunMvOyzFwKXA28HfhqZr6embcAi4Cti8rYscDJmflKZs4Bvl4cv83iYt/FmXkzMBfYbhWv5xvF8RcAHwBmFPEtyVJ/8LWUko2OjgX+JzP/kplLM/MK4HVKSSSZ+fPMfC4zl2Xm1ZQq5+PK4t4c2DQzF2ZmW/V0Vc7f5vLMfLDYfnFm3pSZT2bJH4FbKCVBnTkOOCMzZ2bm68DZwAejQltNZs7KzGszc37x/TiHUkLfZhmwY0SslZnPZ2bFVojMnJaZfy5inkEpwd+j0rZvUMXPRZT+wnI0cFJm/qP4vt1VXHd33g88npmTi7ivAh4BDizb5rLMfKz4LP0MGNXJsT4K/DAz7ynOfTqwW6x6r3WnP4tVuK74PiwErgMWZuaVZT+X7ZXtrj7LmTmf0s/3+cCPgE9nZk98LkBaI5lsS33TZOBISpWsjn+23ohSNfHpsmVPU6oitnmh7OsFAJnZcdk6wMaUqmbTinaD14DfFMvbzMrMJWXv5xf7ropny77eHNil7XzFOT8K/EuF/TYHPtdh27dTqgASER+P5S0mrwE7Uro/AKdSqjz/NSIejIij38D5K8VPRBwQEX8u2hNeo1TZ7bSdoTjndWXnexhYCry144YRsXZE/E/RgjGbUuvFehHRnJnzKFU9jwOej4ibIuKdlU4YEdsWLSj/LI7z9W5iXFWdfS42olTJffINHHNTVvxcw8qf7X9WOGe3x8rMucCsDseqRlc/i93p+DNX6WcQ6PazTGb+Bfg7pc/0z1YxDkldMNmW+qDMfJrSn5zfB/yiw+qXWV61bfMOVq3iVn6sBZRaEdYrXkMys9pkOt/Ads8Cfyw733pZGhHiUxX2exY4p8O2a2fmVRGxOfB94ERKbQbrAQ9QSkbIzH9m5jGZuSnwb8DFRQ/yqpx/pfiL/t5rKY3e8dbivDe3nbeTe/IscECHcw7MzErfs89R+svBLpm5LvDutlMX1/XbzNyHUgvJI8U9qOS/i/XbFMf597IYqzGfFR8G7OqXkXIvAwsptdl01N3n5TlW/FzDG/9sr3CsKPXTb7iqx+rmZ7EmuvssF9ucQKkd5jlKv0hKqhGTbanv+gTwnqKa2a74E/TPgHMiYnDxP+pTKP15eZVk5jJK/5O/ICLeAhARb4uI/ao8xAuUHjRbFb8Cto2I/xcRLcVr54jYvsK23weOi4hdomRQRLw/IgZT6l9O4KUi7qMoVQMp3h8eyx8sfLXYdtkqnr+S/pSSnpeAJRFxAKU++DYvABtGxJCyZZdQ+n5tXsS2cUQc1MnxB1P6Bei1KD1EeVbZNb01Ig4qEsfXKbVuLOviOLOBuUX1u6tfJiqZDhwZEc0RsT9VtqAUn6kfAudH6eHQ5ig9CNl2z5bR+WfmZkrfmyMjol9EfBgYRul7tqquAo6KiFHFub8O/KVoqVlVFX8Wa6i7z/K2lHq8P0apneTUiOisfUbSKjLZlvqooid4aierP03pAaq/A3cAP6GU4LwRpwFPAH8u2g1+R/U92f8LDCv+9L3SWMWVFH3I+1LqC3+OUlvAuZQS2I7bTgWOAf6LUsL8BKU/55OZD1HqVf8TpQR3OHBn2e47A3+J0kgbN1DqIf77qpy/i/g/Q+kXnlcptRjcULb+EUqJ3t+L+7Ip8J1im1siYg6lhyV36eQUF1J6+O/lYrvflK1rovSL1XPAK5QS4M6S6M8Xsc2h9EvL1Z1s15mTKPVKt7XZVPX9LTv334C7izjPBZqK3uNzgDuLe7Nr+U6ZOYtST/3nKLV8nAp8IDNfXsXYyczfAWdS+ivE85Qq7R/pcqfOj9XVz+Kb1tVnuejr/xFwbmbel5mPU/orxeTofBQVSasgMqv9K60kSZKkVWFlW5IkSaoTk21JktZAxSg3cyu8Kk40JGnNZBuJJEmSVCdWtiVJkqQ6MdmW1CtExDci4rMNOvf4iHi0EeeuhYiYEhGf7GTdO4rWheY6nfvXETGxHseup4j4a0Ts0Og4JK35TLYl9XgRsTGlKa//p3i/Z0RkRFzXYbuRxfIptTx/Zt6emas6xXy3itkeL46IlyOiNSJu62LbKRGxsKyvtybJf2Y+U0zKs7QWx6tw/AMy84p6HLujiBgQET+MiNnFzJendLHtjhHx2+LeV+q3/Dbw1fpFK6m3MNmW1BtMAm7OzAVly14CdouIDcuWTQQeW52BvUmXAhsA2xf/ntzN9icWifE69Uj+e4GzgW0ozfy4F6XJW/bvZNvFlMY6/0Qn628A9oqIame+lNRHmWxL6g0OAP7YYdkiShOlfASgaIP4MPDjag5YVMCPj4jHI2JORPxHRGwVEXcVldGfRUT/Yts9I2Jm2b4zIuLzEXF/UZG+OiIGrsoFFbMyTgCOzcyXMnNpZk5blWOsoq2K1ojZEfHLYnZJImJocS/6Fe+nFPfizuK+3BIRG3VzLQMj4kcRMauYbObuiHhr2fE+WXx9X4dRNzIi9izW7Vrc+9eK7fZ8A9c4EfiPzHw1Mx+mNBnPpEobZuajmfm/wIOdrF8ITAOqnQ1VUh9lsi2pNxgOVGqbuJJSewmUkqIHKM2OWK39gJ2AXSnNNngppSmt305puusjutj3Q8D+wBbACIqkruiBfq2L15HF/uOAp4GvFK0Mf4uIw7qJ9xvFtne+gWT048DRwCbAEuC7XWx7JHAU8BZK08t/vptjTwSGULpvGwLHUZoyfgWZObKtMk9pJstHgXsi4m3ATZSmFN+gON+1RfsQRatNZ/fz/mKb9Ytru6/slPcBb6bv+mFg5JvYX1IfYLItqTdYj9K04SvIzLuADSJiO0rJ5JWreNxvZubszHyQUqJ+SzEleyvwa2B0F/t+NzOfy8xXgBuBUUVMz2Tmel28flLsvxmlhL4V2BQ4EbgiIrbv5HynAVsCb6P0S8GNEbHVKlzr5Mx8IDPnUZqG/ENdPBR5WWY+VrTt/Kzt2rqwmFKSvXVbhT4zZ3e2cUT8K6XEekKx3ccotQndnJnLMvNWYCrwPoDMPL6L+zmiOOw6xb+tZadqBQZ3E3tX5lD67ElSp0y2JfUGr9J50jSZUqK6F3BdJ9t05oWyrxdUeL8Onftn2dfzu9m2kgWUktSvZeaizPwj8Adg30obZ+ZfMnNOZr5ePHB4J0UyWqVny75+GmgBOmsPWdVrmwz8FvhpRDwXEd+MiJZKG0bE2ykl8BMzs62/fnPg8PKKNfCvlCrV1Zpb/Ltu2bJ1qfBL2ioYDLz2JvaX1AeYbEvqDe4Htu1k3WTgeEqV0fmrL6TKYvlQep292mYHvL/C7qsyC1kCsQrbv73s63dQSvRfXoX9Ow8kc3FmfiUzhwG7Ax9geXtPu4hYi1Kf/YWZ+euyVc9SqryXV6wHZeZ/Fvtd0sX9fLCI4VXgeVZs+xhJJz3ZVdqeFdtSJGklJtuSeoObgT0qrcjMp4p1Z6zWiDpRNpReZ6+2BzhvA54BTo+IfhHxLkrV+d92PGZErBcR+xUPIvYrEvZ3A78p1rc95Di0i9A+FhHDImJtSkPaXVOr4f4iYq+IGF60pcymlMgvq7DpD4FHMvObHZb/CDiwuMbm4jr3jIjNADLzuC7uZ3lP9pXAlyJi/eIB1GOAyzuJOYqHWtsegh0YEQPK1g+k1M9/6xu4JZL6EJNtSb3BlcD7isroSjLzjsxclQcjGy4zFwMHUWoFaaU0csbHM/MRgIj494hoq/62UOpxfolSNfrTwMFlbRhvp9Qa8o8uTjmZUuL5T2Ag8JkaXs6/ANdQSrQfpjRyzOQK230EOKRDZXp8Zj5L6V78O6VrfBb4Aqv+/7CzgCcp3Ys/At/KzLZfSNr+4vCOYtvNKbXytFW+F7DiQ7gHAlN62udK0uoXmavyV0lJWjNFxNeBFzPzwkbHsqaJiC8BL2Xm/zQ6lt4iIv4CfCIzH2h0LJLWbCbbkiRJUp3YRiKpT4qI8Z09VNfo2HqiiPhoVw8oSlJfZWVbkiRJqhMr25IkSVKd9Gt0ALW00UYb5dChQxsdhiRJknq5adOmvZyZG3e3Xa9KtocOHcrUqVMbHYYkSZJ6uYh4uprtekUbSUQcGBGXtra2NjoUSZIkqV2vSLYz88bMPHbIkCGNDkWSJElq1yuSbUmSJGlN1Kt6tiVJ0uq3ePFiZs6cycKFCxsdilRzAwcOZLPNNqOlpeUN7W+yLUmS3pSZM2cyePBghg4dSkQ0OhypZjKTWbNmMXPmTLbYYos3dAzbSCRJ0puycOFCNtxwQxNt9ToRwYYbbvim/mpjsi1Jkt40E231Vm/2s22yLUmSeryZM2dy0EEHsc0227DVVltx0kknsWjRovb1d9xxB+PGjeOd73wn2223HRdffHGnxxo6dCjjx49fYdmoUaPYcccdAZgyZQpDhgxh1KhRjBgxgr333psXX3yRyy67jFGjRjFq1Cj69+/P8OHDGTVqFF/84hfrc9E1dv311/PQQw+9qWNMmTKFD3zgAzWKqPamTJnCXXfdtVrPabItSZJ6tMzk0EMP5eCDD+bxxx/nscceY+7cuZxxxhkA/POf/+TII4/kkksu4ZFHHuHOO+/kf//3f7nuuus6PeacOXN49tlnAXj44YdXWj9+/HimT5/O/fffz84778z3vvc9jjrqKKZPn8706dPZdNNN+cMf/sD06dP5z//8z/pceI1Vm2wvWbJkNUTzxnUVn8m2JEnSKvq///s/Bg4cyFFHHQVAc3MzF1xwAT/84Q+ZP38+3/ve95g0aRJjxowBYKONNuKb3/wm3/rWtzo95oc+9CGuvvpqAK666iqOOOKIittlJnPmzGH99devOt7LL7+cgw8+mH322YehQ4fyX//1X5x//vmMHj2aXXfdlVdeeQWA6dOns+uuuzJixAgOOeQQXn31VQD23HNPTj75ZMaOHcv222/P3XffzaGHHso222zDl770pfbz/OhHP2LcuHGMGjWKf/u3f2Pp0qUArLPOOpxxxhmMHDmSXXfdlRdeeIG77rqLG264gS984QuMGjWKJ598coWYJ02axHHHHccuu+zCqaeeyrx58zj66KMZN24co0eP5pe//OVK19nZNjNmzGD8+PGMGTOGMWPGtCe/zz//PO9+97vb/4pw++23A3DLLbew2267MWbMGA4//HDmzp270rn23HNPPvvZzzJ27Fi+853vcOONN7LLLrswevRo9t57b1544QVmzJjBJZdcwgUXXMCoUaO4/fbbeemllzjssMPYeeed2Xnnnbnzzjur/j5Wy9FIJElSzXzlxgd56LnZNT3msE3X5awDd+h0/YMPPshOO+20wrJ1112Xd7zjHTzxxBM8+OCDTJw4cYX1Y8eO7bKKe9hhh3HUUUfx+c9/nhtvvJEf//jHTJ48uX397bffzqhRo5g1axaDBg3i61//+ipd0wMPPMC9997LwoUL2XrrrTn33HO59957Ofnkk7nyyiv57Gc/y8c//nEuuugi9thjD7785S/zla98hQsvvBCA/v37M3XqVL7zne9w0EEHMW3aNDbYYAO22morTj75ZF588UWuvvpq7rzzTlpaWjj++OP58Y9/zMc//nHmzZvHrrvuyjnnnMOpp57K97//fb70pS8xYcIEPvCBD/DBD36wYswzZ87krrvuorm5mX//93/nPe95Dz/84Q957bXXGDduHHvvvfcK259zzjkVt3nLW97CrbfeysCBA3n88cc54ogjmDp1Kj/5yU/Yb7/9OOOMM1i6dCnz58/n5Zdf5mtf+xq/+93vGDRoEOeeey7nn38+X/7yl1eKb9GiRUydOhWAV199lT//+c9EBD/4wQ/45je/yXnnncdxxx3HOuusw+c//3kAjjzySE4++WT+9V//lWeeeYb99tuv4l8y3gyTbUmSpA423HBD1l9/fX7605+y/fbbs/baa6+wfvz48fzqV78C4Nxzz+XUU0/lkksuqfr4e+21F4MHD2bw4MEMGTKEAw88EIDhw4dz//3309raymuvvcYee+wBwMSJEzn88MPb958wYUL79jvssAObbLIJAFtuuSXPPvssd9xxB9OmTWPnnXcGYMGCBbzlLW8BSol6W1/1TjvtxK233lpVzIcffjjNzc1Aqdp8ww038O1vfxsojUjzzDPPrLB9Z9tsuummnHjiiUyfPp3m5mYee+wxAHbeeWeOPvpoFi9ezMEHH8yoUaP44x//yEMPPcS73vUuoJRQ77bbbhXj+/CHP9z+9cyZM/nwhz/M888/z6JFizodtu93v/vdCr90zZ49m7lz57LOOutUdU+qYbItSZJqpqsKdL0MGzaMa665ZoVls2fP5plnnmHrrbdm2LBhTJs2jYMOOqh9/bRp0xg7dixLly5tr4pPmDCBr371q+3bfPjDH+aEE07g8ssv7/L8EyZM4LDDDlulmAcMGND+dVNTU/v7pqamqnqiy7fveKwlS5aQmUycOJFvfOMbK+3b0tLSPsJGc3Nz1T3YgwYNav86M7n22mvZbrvtVtjmhRde6Habs88+m7e+9a3cd999LFu2jIEDBwLw7ne/m9tuu42bbrqJSZMmccopp7D++uuzzz77cNVVV61SfJ/+9Kc55ZRTmDBhAlOmTOHss8+uuM+yZcv485//3B5DPdizLUmSerT3vve9zJ8/nyuvvBKApUuX8rnPfY5Jkyax9tprtyfM06dPB2DWrFmcccYZnHnmmTQ3N7c/1FieaAMccsghnHrqqey3335dnv+OO+5gq622quk1DRkyhPXXX7+9b3ny5MntVe5qvPe97+Waa67hxRdfBOCVV17h6aef7nKfwYMHM2fOnKqOv99++3HRRReRmQDce++9VW/T2trKJptsQlNTE5MnT27vJX/66ad561vfyjHHHMMnP/lJ7rnnHnbddVfuvPNOnnjiCaDUB95WCe9Ka2srb3vb2wC44oorOr3Gfffdl4suuqj9fdtnpJZMtiVJUo8WEVx33XX8/Oc/Z5tttmHbbbdl4MCB7X3Um2yyCT/60Y849thj2W677dh00035zGc+023yOnjwYE477TT69++/0rq2nu2RI0cyefJkzjvvvJpf1xVXXMEXvvAFRowYwfTp0yv2KXdm2LBhfO1rX2PfffdlxIgR7LPPPjz//PNd7vORj3yEb33rW4wePXqlByQ7OvPMM1m8eDEjRoxghx124Mwzz6x6m+OPP54rrriCkSNH8sgjj7RXpKdMmcLIkSMZPXo0V199NSeddBIbb7wxl19+OUcccQQjRoxgt91245FHHun2+s8++2wOP/xwdtppJzbaaKP25QceeCDXXXdd+wOS3/3ud5k6dSojRoxg2LBhq9QKVK1o+22jNxg7dmy2NcZLkqTV4+GHH2b77bdvdBhVu/jii/nv//5vbrvttlUaRUR9V6XPeERMy8yx3e1rZVuSJPUpxx9/PH/7299MtLVamGxLkiRJdWKyLUmSJNWJyXYNfPyHf+XLv3yg0WFIkiRpDeM42zXw/GsLGNS/udFhSJIkaQ1jZbsGIqAXDeoiSZKkGjHZroGmCBKzbUmSGiUi+NjHPtb+fsmSJWy88cbt05KvbtOnT+fmm29uyLnfrD333JPuhlK+8MILmT9/fvv7973vfbz22mv1Dq1dT7q/Jts1ssxcW5Kkhhk0aBAPPPAACxYsAODWW29tn0GwEXpSMvhGdEy2b775ZtZbb72anqOraeR70v012a6BiLCNRJKkBnvf+97HTTfdBMBVV13FEUcc0b7ulVde4eCDD2bEiBHsuuuu3H///UBppsGJEycyfvx4Nt98c37xi19w6qmnMnz4cPbff38WL14MwLRp09hjjz3Yaaed2G+//dpnY9xzzz057bTTGDduHNtuuy233347ixYt4stf/jJXX301o0aN4uqrr14hzhkzZjB+/HjGjBnDmDFjuOuuu9rXnXvuuQwfPpyRI0fyxS9+EYAnnniCvffem5EjRzJmzBiefPJJpkyZskLV/sQTT+Tyyy8HYOjQoZx++umMGjWKsWPHcs8997Dffvux1VZbtc+Q2NX+5T71qU8xduxYdthhB8466ywAvvvd7/Lcc8+x1157sddee7Wf8+WXXwbg/PPPZ8cdd2THHXfkwgsvbL/m7bffnmOOOYYddtiBfffdt/0Xo3KTJk3iuOOOY5ddduHUU0/lr3/9K7vtthujR49m991359FHH614f+fNm8fRRx/NuHHjGD16NL/85S87/6CsbpnZa1477bRTNsIBF96Wn7j8rw05tyRJjfbQQw81OoQcNGhQ3nfffXnYYYflggULcuTIkfmHP/wh3//+92dm5oknnphnn312Zmb+/ve/z5EjR2Zm5llnnZXvete7ctGiRTl9+vRca6218uabb87MzIMPPjivu+66XLRoUe6222754osvZmbmT3/60zzqqKMyM3OPPfbIU045JTMzb7rppnzve9+bmZmXXXZZnnDCCRVjnTdvXi5YsCAzMx977LFsy19uvvnm3G233XLevHmZmTlr1qzMzBw3blz+4he/yMzMBQsW5Lx581a4tszME044IS+77LLMzNx8883z4osvzszMz372szl8+PCcPXt2vvjii/mWt7wlM7PL/ffYY4+8++67V4hhyZIluccee+R9993Xfo6XXnqpff+291OnTs0dd9wx586dm3PmzMlhw4blPffck0899VQ2Nzfnvffem5mZhx9+eE6ePHmlezNx4sR8//vfn0uWLMnMzNbW1ly8eHFmZt5666156KGHVry/p59+evvxXn311dxmm21y7ty5Fe//G1HpMw5MzSryU0cjqYGmJh+QlCSp3eWXw4wZtTve0KEwaVK3m40YMYIZM2Zw1VVX8b73vW+FdXfccQfXXnstAO95z3uYNWsWs2fPBuCAAw6gpaWF4cOHs3TpUvbff38Ahg8fzowZM3j00Ud54IEH2GeffQBYunQpm2yySfuxDz30UAB22mknZlRx3YsXL+bEE09k+vTpNDc389hjjwHwu9/9jqOOOoq1114bgA022IA5c+bwj3/8g0MOOQSAgQMHdnt8gAkTJrRfw9y5cxk8eDCDBw9mwIABq9Rb/bOf/YxLL72UJUuW8Pzzz/PQQw8xYsSITre/4447OOSQQxg0aBBQuje33347EyZMYIsttmDUqFFA1/fq8MMPp7m5NMpba2srEydO5PHHHyci2v/S0NEtt9zCDTfcwLe//W0AFi5cyDPPPLPSFOuNYLJdA0GwzGxbkqSSKhLjepkwYQKf//znmTJlCrNmzapqnwEDBgDQ1NRES0sLEdH+fsmSJWQmO+ywA3/605+63L+5ubnLPuM2F1xwAW9961u57777WLZsWdUJdLl+/fqxbNmy9vcLFy7s9Jravi6/pu72B3jqqaf49re/zd13383666/PpEkB7jvFAAAgAElEQVSTKm5XrfI4mpubK7aRAO2JOsCZZ57JXnvtxXXXXceMGTPYc889K+6TmVx77bVst912bzi+erFnuwYicCwSSZLWAEcffTRnnXUWw4cPX2H5+PHj+fGPfwyU+pU32mgj1l133aqOud122/HSSy+1J9uLFy/mwQcf7HKfwYMHM2fOnIrrWltb2WSTTWhqamLy5MksXboUgH322YfLLrus/cHDV155hcGDB7PZZptx/fXXA/D6668zf/58Nt98cx566CFef/11XnvtNX7/+99XdS1tqtl/9uzZDBo0iCFDhvDCCy/w61//utvrGz9+PNdffz3z589n3rx5XHfddYwfP36VYivX2tra/qBreU95x/Pvt99+XHTRRWRR/Lz33nvf8DlrzWS7BnxAUpKkNcNmm23GZz7zmZWWn3322UybNo0RI0bwxS9+kSuuuKLqY/bv359rrrmG0047jZEjRzJq1KgVHmqsZK+99uKhhx6q+IDk8ccfzxVXXMHIkSN55JFH2iu5+++/PxMmTGDs2LGMGjWqvSVi8uTJfPe732XEiBHsvvvu/POf/+Ttb387H/rQh9hxxx350Ic+xOjRo6u+HqCq/UeOHMno0aN55zvfyZFHHsm73vWu9nXHHnss+++/f/sDkm3GjBnDpEmTGDduHLvssguf/OQnVzm2cqeeeiqnn346o0ePXuGvBh3v75lnnsnixYsZMWIEO+ywA2eeeeYbPmetRfaiLHHs2LHZ3biQ9XDw9+5k8MB+TP7ELqv93JIkNdrDDz+8RvTGSvVS6TMeEdMyc2x3+1rZroGitUuSJElagcl2DQSORiJJkqSVmWzXgNO1S5IkqRKT7RqIgLLRcyRJ6nN60zNgUrk3+9k22a6BwMq2JKnvGjhwILNmzTLhVq+TmcyaNesNjYXexkltaiDCnm1JUt+12WabMXPmTF566aVGhyLV3MCBA9lss83e8P4m2zVgG4kkqS9raWlhiy22aHQY0hqpV7SRRMSBEXFpa2trY85vG4kkSZIq6BXJdmbemJnHDhkypCHnb2qyjUSSJEkr6xXJdqMFwTKzbUmSJHVgsl0DEdhEIkmSpJWYbNeIhW1JkiR1ZLJdA6UZJCVJkqQVmWzXQGmcbdNtSZIkrchkuwYC20gkSZK0MpPtGii1kZhtS5IkaUUm2zXgDJKSJEmqxGS7JnxAUpIkSSsz2a6BJh+QlCRJUgUm2zVQGo2k0VFIkiRpTWOyXQOBD0hKkiRpZSbbNWBlW5IkSZWYbNeAM0hKkiSpEpPtWghYZmlbkiRJHZhs10AAlrYlSZLUkcl2DdhGIkmSpEpMtmsgbCORJElSBSbbNRA4GokkSZJWZrJdAxGOsy1JkqSVmWzXQAQsW9boKCRJkrSmMdmugSiNRyJJkiStwGS7BkozSNpGIkmSpBWZbNdAUzjMtiRJklZmsl0DQTj0nyRJklZisl0DpTaSRkchSZKkNY3Jdg2EM0hKkiSpApPtGvABSUmSJFVisl0DziApSZKkSky2ayAcjUSSJEkVmGzXQFOEbSSSJElaicl2DQSwzFxbkiRJHZhs10BY2ZYkSVIFJts1YM+2JEmSKjHZroEgHI1EkiRJKzHZrgHH2ZYkSVIlJts10GQbiSRJkiow2a6BiGCZlW1JkiR1YLJdA84gKUmSpEpMtmvBNhJJkiRVYLJdA02O/SdJkqQKTLZroDSDpNm2JEmSVmSyXQMWtiVJklSJyXYNNDlduyRJkiow2a6BUhtJo6OQJEnSmsZkuxYiGh2BJEmS1kAm2zXQVOTatpJIkiSpnMl2DQSlbNtWEkmSJJUz2a6BsLItSZKkCky2a6CtY9tUW5IkSeVMtmugqWjatrAtSZKkcibbNeQskpIkSSpnsl0DjvwnSZKkSky2a6ApbCORJEnSyky2a6CtsG0biSRJksqZbNdA+9B/jQ1DkiRJaxiT7Rpom9TGcbYlSZJUzmS7BqxsS5IkqRKT7RqItgcklzU4EEmSJK1RTLZrYPkMkta2JUmStJzJdg00tbWRmGtLkiSpTL9GB9AmIrYEzgCGZOYHi2XjgY9SinNYZu7ewBA71dZG4tB/kiRJKlfXynZE/DAiXoyIBzos3z8iHo2IJyLiiwCZ+ffM/ET5dpl5e2YeB/wKuKKesb4ZPiApSZKkSurdRnI5sH/5gohoBr4HHAAMA46IiGHdHOdI4Cf1CLAWwhkkJUmSVEFdk+3MvA14pcPiccATRSV7EfBT4KDOjhER7wBaM3NOJ+uPjYipETH1pZdeqlXoq6T9AUmzbUmSJJVpxAOSbwOeLXs/E3hbRGwYEZcAoyPi9LL1nwAu6+xgmXlpZo7NzLEbb7xxfSLuhm0kkiRJqmSNeUAyM2cBx1VYflYDwlkly2eQbHAgkiRJWqM0orL9D+DtZe83K5b1WO1D/1nbliRJUplGJNt3A9tExBYR0R/4CHBDA+KombY2kmXm2pIkSSpT76H/rgL+BGwXETMj4hOZuQQ4Efgt8DDws8x8sJ5x1NvyNhKzbUmSJC1X157tzDyik+U3AzfX89yrUziDpCRJkipwuvYacJxtSZIkVWKyXQPt42z7gKQkSZLKmGzXQFNxF61sS5IkqVyvSLYj4sCIuLS1tbUx5y9q28vMtiVJklSmVyTbmXljZh47ZMiQhpzfGSQlSZJUSa9IttcUFrYlSZJUzmS7Bppi+SOSkiRJUhuT7RpwBklJkiRVYrJdA8tnkGxwIJIkSVqjmGzXQFP7A5Jm25IkSVrOZLsG2ttIljU2DkmSJK1ZTLZromgjsbItSZKkMibbNdDeRmKuLUmSpDIm2zUQ4QOSkiRJWpnJdg00tQ/9Z7YtSZKk5XpFsh0RB0bEpa2trQ05f1ORbS812ZYkSVKZXpFsZ+aNmXnskCFDGnL+5qKNZJmz2kiSJKlMr0i2G625qGyba0uSJKmcyXYNtI2zvdRsW5IkSWVMtmugvY3Enm1JkiSVMdmugbY2EivbkiRJKmeyXQNNTVa2JUmStDKT7Rposo1EkiRJFZhs10Bbz/bSZQ0ORJIkSWsUk+0aaCruoj3bkiRJKmeyXQPN9mxLkiSpApPtGnDoP0mSJFXSK5LtiDgwIi5tbW1t1PkB20gkSZK0ol6RbGfmjZl57JAhQxpyfttIJEmSVEmvSLYbzdFIJEmSVInJdg20jUayzDYSSZIklTHZrgHbSCRJklSJyXYNtM0gudRkW5IkSWVMtmugfbp220gkSZJUxmS7BtraSBz6T5IkSeVMtmtg+aQ2DQ5EkiRJaxST7RqIttFI7NmWJElSmX5drYyI3YCPAeOBTYAFwAPATcCPMrMxUzauYZqdQVKSJEkVdFrZjohfA58EfgvsTynZHgZ8CRgI/DIiJqyOINd07T3bVrYlSZJUpqvK9v/LzJc7LJsL3FO8zouIjeoWWQ/iaCSSJEmqpNPKdma+HBHNEfGHrrapT1g9y/JJbRociCRJktYoXT4gmZlLgWURMWQ1xdMjFbm2PduSJElaQZcPSBbmAn+LiFuBeW0LM/MzdYtqFUXEgcCBW2+9daPOT4SjkUiSJGlF1STbvyhea6zMvBG4cezYscc0KobmCCvbkiRJWkE1yfYrwE2ZuazewfRkTU3haCSSJElaQTWT2nwYeDwivhkR76x3QD1VcwTm2pIkSSrXbbKdmR8DRgNPApdHxJ8i4tiIGFz36HqQpvABSUmSJK2oqunaM3M2cA3wU0qT2xwC3BMRn65jbD1KU5M925IkSVpRt8l2REyIiOuAKUALMC4zDwBGAp+rb3g9R3NTOBqJJEmSVlDNA5KHARdk5m3lCzNzfkR8oj5h9TzNYbItSZKkFXWbbGfmxC7W/b624fRcEcFSx2uRJElSmap6ttW95iZYZs+2JEmSyphs10hzOM62JEmSVtRpsh0Rl0bEIQ7xV52mprCyLUmSpBV0Vdn+X0ojjtwcEb+PiNMiYuRqiqvHcTQSSZIkddTpA5KZ+RfgL8DZEbEhsC/wuYgYDtwL/CYzf7Z6wlzzNUWw1FxbkiRJZaoZ+o/MnAVcVbyIiJ2A/esYV4/TFD4gKUmSpBVVlWx3lJnTgGk1jqVHa3YGSUmSJHXgaCQ10uRoJJIkSeqgVyTbEXFgRFza2trasBiam4I02ZYkSVKZbttIIqIZeD8wtHz7zDy/fmGtmsy8Ebhx7NixxzQqhqawjUSSJEkrqqZn+0ZgIfA3wAnJO9HU5GgkkiRJWlE1yfZmmTmi7pH0cM2ORiJJkqQOqunZ/nVE7Fv3SHq4fk1NLFlm4V+SJEnLVVPZ/jNwXUQ0AYuBADIz161rZD1Mc1OYbEuSJGkF1STb5wO7AX9Lh9voVL/mYOESb48kSZKWq6aN5FngARPtrvVzUhtJkiR1UE1l++/AlIj4NfB628I1aei/NUFzUxNLHI5EkiRJZapJtp8qXv2Llyqwsi1JkqSOuk22M/MrqyOQnq65OVjsA5KSJEkq0yuma18TtFjZliRJUgcm2zViz7YkSZI6MtmuEXu2JUmS1FGnPdsR8eUu9svM/I86xNNjNTcHS0y2JUmSVKarByTnVVi2NvBJYEPAZLtMP2eQlCRJUgedJtuZeV7b1xExGDgJOBr4KXBeZ/v1Vf2amlhqz7YkSZLKdDn0X0RsAJwCfBS4AhiTma+ujsB6mn62kUiSJKmDrnq2vwUcClwKDM/Muastqh6o2QckJUmS1EFXo5F8DtgU+BLwXETMLl5zImL26gmv57BnW5IkSR111bPtsICroLkpWJawbFnS1BSNDkeSJElrgE4T6ohYp7udq9mmr2hpLt1K+7YlSZLUpqvq9S8j4ryIeHdEDGpbGBFbRsQnIuK3wP71D7F7EXFgRFza2trasBiai2q2fduSJElq02mynZnvBX4P/BvwYES0RsQs4EfAvwATM/Oa1RNm1zLzxsw8dsiQIQ2LoV+RbNu3LUmSpDZdDv2XmTcDN6+mWHo0K9uSJEnqyIcga6Stsr3YiW0kSZJUMNmukeam0q20si1JkqQ2Jts10q/Znm1JkiStqKpkOyL+NSKOKr7eOCK2qG9YPU8/e7YlSZLUQbfJdkScBZwGnF4saqE0IonKNLePRmKyLUmSpJJqKtuHABOAeQCZ+RwwuJ5B9UT9ip7tJT4gKUmSpEI1yfaizEwgAconuNFyzY6zLUmSpA6qSbZ/FhH/A6wXEccAvwN+UN+wep6WZnu2JUmStKIuJ7UByMxvR8Q+wGxgO+DLmXlr3SPrYezZliRJUkfdJtsRcW5mngbcWmGZCv0cZ1uSJEkdVNNGsk+FZQfUOpCerrl9Bkl7tiVJklTSaWU7Ij4FHA9sGRH3l60aDNxZ78B6mv79ijYSRyORJElSoas2kp8Avwa+AXyxbPmczHylrlH1QC3NpT8SWNmWJElSm06T7cxsBVqBIwAi4i3AQGCdiFgnM59ZPSH2DCbbkiRJ6qiaGSQPjIjHgaeAPwIzKFW8VaYt2V5kG4kkSZIK1Twg+TVgV+CxzNwCeC/w57pG1QP1b6tsL7GyLUmSpJJqku3FmTkLaIqIpsz8AzC2znH1OC39HI1EkiRJK+p2nG3gtYhYB7gN+HFEvAjMq29YPY8925IkSeqomsr2QcB84GTgN8CTwIH1DKonsmdbkiRJHXVZ2Y6IZuBXmbkXsAy4YrVE1QO1NLeNs21lW5IkSSVdVrYzcymwLCKGrKZ4eizbSCRJktRRNT3bc4G/RcStlPVqZ+Zn6hZVD9SvmK7dNhJJkiS1qSbZ/kXxUhcigv7NTVa2JUmS1K7bZDsz7dOuUktzOM62JEmS2lUzGskar5jl8tLW1taGxtHSz8q2JEmSlusVyXZm3piZxw4Z0tjnOFuam+zZliRJUruqk+2IWLuegfQG9mxLkiSpXLfJdkTsHhEPAY8U70dGxMV1j6wHamkOk21JkiS1q6ayfQGwHzALIDPvA95dz6B6qhYr25IkSSpTVRtJZj7bYdHSOsTS47U0N7FoiT3bkiRJKqlmnO1nI2J3ICOiBTgJeLi+YfVMjkYiSZKkctVUto8DTgDeBvwDGFW8Vwf97dmWJElSmWomtXkZ+OhqiKXHs2dbkiRJ5aoZjeSbEbFuRLRExO8j4qWI+NjqCK6ncZxtSZIklaumjWTfzJwNfACYAWwNfKGeQfVULc1NTtcuSZKkdtUk222tJu8Hfp6ZjZ0TfQ3Wv58925IkSVqumtFIfhURjwALgE9FxMbAwvqG1TO1NDexZJltJJIkSSrptrKdmV8EdgfGZuZiYB5wUL0D64lK42xb2ZYkSVJJNZVtgE2BvSNiYNmyK+sQT4/maCSSJEkq122yHRFnAXsCw4CbgQOAOzDZXonjbEuSJKlcNQ9IfhB4L/DPzDwKGAkMqWtUPVSpsm3PtiRJkkqqSbYXZOYyYElErAu8CLy9vmH1TC39mlhkZVuSJEmFanq2p0bEesD3gWnAXOBPdY2qh2rr2c5MIqLR4UiSJKnBqpmu/fjiy0si4jfAupl5f33D6pn6NweZsHRZ0q/ZZFuSJKmvq2a69oiIj0XElzNzBvBaRIyrf2g9T0tz6XbaSiJJkiSormf7YmA34Iji/Rzge3WLqAfr3690Oxcv8SFJSZIkVdezvUtmjomIewEy89WI6F/nuHqkAf2aAXh9yVKgpbHBSJIkqeGqqWwvjohmIAGK6drtk6hgQFHZft1ZJCVJkkR1yfZ3geuAt0TEOZQmtPl6XaPqoQa0lG7nwsVLGxyJJEmS1gTVjEby44iYRmlimwAOzsyH6x5ZD7S8jcTKtiRJkrpItiNig7K3LwJXla/LzFfqGVhPtLyNxMq2JEmSuq5sT6PUp10+YHTb+wS2rGNcPVJ7sr3YyrYkSZK6SLYzc4vVGUhvMKDFNhJJkiQttyqT2pxZvH+Hk9pUZhuJJEmSyq3KpDZHFu+d1KYTDv0nSZKkck5qU0PtbST2bEuSJAkntakp20gkSZJUzkltasg2EkmSJJVzUpsaGli0kTiDpCRJkqCbZLtoH3kwM98JPLJ6Quq5+jUFTWFlW5IkSSVdtpFk5lLg0Yh4x2qKp0eLCAb0azbZliRJElDdaCTrAw9GxF+BeW0LM3NC3aLqwQa0NPG6bSSSJEmiumT7zLpH0YsM6NdkZVuSJElAdQ9I/nF1BPJmRMSBwIFbb711o0OxjUSSJEntqhn6b42XmTdm5rFDhgxpdChFZds2EkmSJPWSZHtNUurZtrItSZKkKpLtiDipmmUqsY1EkiRJbaqpbE+ssGxSjePoNWwjkSRJUptOH5CMiCOAI4EtIuKGslWDgVfqHVhPNaBfE3MWLml0GJIkSVoDdDUayV3A88BGwHlly+cA99czqJ6s1EZiZVuSJEldJNuZ+TTwNLDb6gun5xvY0sRCH5CUJEkS1T0geWhEPB4RrRExOyLmRMTs1RFcT7RW/2YWOIOkJEmSqG4GyW8CB2bmw/UOpjdYq6UfCxaZbEuSJKm60UheMNGu3lr9m1iweCmZ2ehQJEmS1GDVVLanRsTVwPXA620LM/MXdYuqB1u7fz+WLksWLV3GgH7NjQ5HkiRJDVRNsr0uMB/Yt2xZAibbFQxsKSXYCxeZbEuSJPV13SbbmXnU6gikt1i7fynBnr94CUNoaXA0kiRJaqRqRiPZNiJ+HxEPFO9HRMSX6h9az7RWUdn2IUlJkiRV84Dk94HTgcUAmXk/8JF6BtWTrdVW2TbZliRJ6vOqSbbXzsy/dljmfOSdaK9sO9a2JElSn1dNsv1yRGxF6aFIIuKDlKZxVwVtPdu2kUiSJKma0UhOAC4F3hkR/wCeAj5W16h6sLbRSGwjkSRJUjWjkfwd2DsiBgFNmTmn/mH1XG2V7YW2kUiSJPV53SbbEbEe8HFgKNAvIgDIzM/UNbIeygckJUmS1KaaNpKbgT8DfwOW1Tecnm/tltIt9QFJSZIkVZNsD8zMU+oeSS8xsH/pmdMFixywRZIkqa+rZjSSyRFxTERsEhEbtL3qHlkP1b+5ieamsLItSZKkqirbi4BvAWdQDP9X/LtlvYLqySKCtVqa7dmWJElSVcn254CtM/PlegfTW6zVv9nRSCRJklRVG8kTwPx6B9KbWNmWJEkSVFfZngdMj4g/AK+3LXTov86t3b/ZGSQlSZJUVbJ9ffFSlQa2NPuApCRJkqqaQfKKiOgPbFssejQzF9c3rJ7NyrYkSZKguhkk9wSuAGYAAbw9IiZm5m31Da3nWqulmdfm+/uIJElSX1dNG8l5wL6Z+ShARGwLXAXsVM/AejJHI5EkSRJUNxpJS1uiDZCZjwEt9Qup53M0EkmSJEF1le2pEfED4EfF+48CU+sXUs83aEA/5jlduyRJUp9XTbL9KeAEoG2ov9uBi+sWUS+wzoB+zHt9CZlJRDQ6HEmSJDVINaORvA6cX7xUhUED+rEsYcHipazdv5rfZyRJktQbdduzHREfiIh7I+KViJgdEXMiYvbqCK6nWmdAMwBzX7eVRJIkqS+rpux6IXAo8LfMzDrH0ysMGlC6rfNeXwqDGxyMJEmSGqaa0UieBR4w0a7eOu3JtpVtSZKkvqyayvapwM0R8Ufg9baFmWkPdyfaku05C022JUmS+rJqku1zgLnAQKB/fcPpHQZZ2ZYkSRLVJdubZuaOdY+kF1lnYJFsO9a2JElSn1ZNz/bNEbFv3SPpRWwjkSRJElSXbH8K+E1ELHDov+rYRiJJkiSoblIbB69bRWu3NBNhsi1JktTXVVPZ1ipqagoG9e/HHJNtSZKkPs1ku04GDWi2si1JktTHmWzXyaAB/UozSEqSJKnP6rRnOyI26GrHzHyl9uH0HoMH2EYiSZLU13X1gOQ0IIGosC6BLesSUS9RqmybbEuSJPVlnSbbmbnF6gyktxk0oB+vzJvf6DAkSZLUQNXMIElErA9sQ2nKdgAy87Z6BdUbDB7Qj7lWtiVJkvq0bpPtiPgkcBKwGTAd2BX4E/Ce+obWsw0y2ZYkSerzqhmN5CRgZ+DpzNwLGA28VteoeoF1BvZj7sIlZGajQ5EkSVKDVJNsL8zMhQARMSAzHwG2q29YPd+6A1tYsixZsNjh/yRJkvqqanq2Z0bEesD1wK0R8SrwdH3D6vmGrNUCQOuCxazdv6rWeEmSJPUy3WaBmXlI8eXZEfEHYAjwm7pG1Qusu1bp1s5esIRNhjQ4GEmSJDVENQ9IvqPs7VPFv/8CPFOXiHqJtsr27IWLGxyJJEmSGqWa/oabWD65zUBgC+BRYIc6xtXjrTuwaCOZb7ItSZLUV1XTRjK8/H1EjAGOr1tEvcS6VrYlSZL6vGpGI1lBZt4D7FKHWHqV8gckJUmS1DdV07N9StnbJmAM8FzdIuolBg9c/oCkJEmS+qZqerYHl329hFIP97W1DiQitgTOAIZk5geLZU3AfwDrAlMz84pan7deWpqbGNS/2cq2JElSH1ZNz/ZX3ujBI+KHwAeAFzNzx7Ll+wPfAZqBH2Tmf2bm34FPRMQ1ZYc4iNI08bOAmW80jkZZd60We7YlSZL6sG57tiNi24i4NCJuiYj/a3tVefzLgf07HK8Z+B5wADAMOCIihnWy/3bAXZl5CvCpKs+5xlh3YAuzrWxLkiT1WdW0kfwcuAT4AbBKc49n5m0RMbTD4nHAE0Ulm4j4KaUK9kMVDjETWFR83ePmPR+yVottJJIkSX1YNcn2ksz87xqe823As2XvZwK7RMSGwDnA6Ig4PTO/AfwCuCgixgO3VTpYRBwLHAvwjne8o9ImDbPuWv34x2sLGx2GJEmSGqSaZPvGiDgeuA54vW1hZr5Sy0AycxZwXIdl84FPdLPfpcClAGPHjs1axvRmrbtWCw8/P6fRYUiSJKlBqkm2Jxb/fqFsWQJbvsFz/gN4e9n7zYplvY4925IkSX1bNaORbFHjc94NbBMRW1BKsj8CHFnjc6wRhqzVwpzXl7B0WdLcFI0OR5IkSatZNZVtImJ3YOj/b+++w6Os8jaOf89kUkglISH0XkKXLr1ZQBd776BiV6you6+r26xrWQuKKIJdEUVQQKX33gWkhd5CCTWkzHn/mCSEkDKTzCQxuT/XtdeayTDnPMmTmfs5z++ck/P51toxHvy7L4HeQKwxZifwd2vtR8aYB4EpuJf++9hau9b7rpd9WVu2H0tJo3JoUCn3RkRERERKmic7SH4KNARWcGZFEAsUGrattTfm8/jPwM+ed/PPKTJzF8nkUwrbIiIiIhWRJyPbHYDm1toyNfnwzyAqc2RbW7aLiIiIVEyFbmoDrAGq+bsj5VFWGYnW2hYRERGpmDwZ2Y4FfjfGLOLspf8u81uvyonozNKRI6dSC3mmiIiIiJRHnoTt5/3difIqOsw9sn34hMK2iIiISEXkydJ/M3N+bYzpDtwIzMz7X5Q8Y8xAYGCjRo1KuytnqVzJPbJ9+KTKSEREREQqIk9qtjHGtDXGvGqMSQT+Cazza6+8ZK2dYK0dEhUVVdpdOUuQ00FEsJNDGtkWERERqZDyHdk2xjTBPYJ9I5AEfA0Ya22fEupbuRAdFsThkwrbIiIiIhVRQWUk64HZwF+stZsAjDGPlkivypHosCCNbIuIiIhUUAWVkVwF7AGmG2M+NMb0A7TnuJdiQgM1si0iIiJSQeUbtq21P1hrbwASgOnAUKCqMWa4Meaikurgn110WBCHT2iCpIiIiEhFVOgESWvtCWvtF9bagUAtYDkwzO89KydiQlWzLSIiIlJRebQaSRZr7WFr7QhrbT9/dai8iQ4L4mRqBilpGaXdFREREREpYV6FbfFeTFjWWtsa3RYRERGpaBS2/Sxry3atSCIiIiJS8Shs+1n2yLYmSYqIiIhUOOUibBtjBhpjRiQnJ5d2V84RExYIwCGVkYiIiIhUOOUibJfV7drhTBnJYZWRiIiIiFQ45SJsl2VRlQIxRr+FVUcAACAASURBVBMkRURERCoihW0/cwY4iAwJ1ARJERERkQpIYbsEVAkP4qDCtoiIiEiFo7BdAuLCgzlw7HRpd0NERERESpjCdgmIjQgmSWFbREREpMJR2C4BceHBHDiusC0iIiJS0Shsl4C4iGCOpaSTkpZR2l0RERERkRKksF0C4sKDAVS3LSIiIlLBKGyXgLgId9hOUimJiIiISIWisF0CYjWyLSIiIlIhKWyXgKyRbU2SFBEREalYykXYNsYMNMaMSE5OLu2u5KlKeBAASce0sY2IiIhIRVIuwra1doK1dkhUVFRpdyVPgQEOokMDOXA8pbS7IiIiIiIlqFyE7T+DuAjtIikiIiJS0Shsl5C4iGCSjquMRERERKQiUdguIbHhGtkWERERqWgUtktInMK2iIiISIWjsF1C4iKCOZWWwYnT6aXdFREREREpIQrbJUQb24iIiIhUPArbJaRqpDts7zuq5f9EREREKgqF7RJSPSoEgL0K2yIiIiIVhsJ2CYmPzAzbyQrbIiIiIhWFwnYJiQgJJDzYqZFtERERkQpEYbsEVYsK0ci2iIiISAWisF2CqkeFsEdhW0RERKTCUNguQfGRGtkWERERqUjKRdg2xgw0xoxITk4u7a4UqHpUCAeOnyY9w1XaXRERERGRElAuwra1doK1dkhUVFRpd6VA1aJCyHBZko6nlnZXRERERKQElIuw/WeRtdb2nuRTpdwTERERESkJCtslKGutbe0iKSIiIlIxKGyXoOpRlQC0IomIiIhIBaGwXYKiQwMJcjq0IomIiIhIBaGwXYKMMVprW0RERKQCUdguYVprW0RERKTiUNguYTUrV2LXEa1GIiIiIlIRKGyXsFrRldh7NEUb24iIiIhUAArbJaxWdCUyXFZ12yIiIiIVgMJ2CasVHQqgUhIRERGRCkBhu4TVrOxea3vnYYVtERERkfJOYbuEVa8cgjGw8/DJ0u6KiIiIiPiZwnYJC3YGEB8RopFtERERkQpAYbsU1IquxC6FbREREZFyT2G7FNSKrsTOIyojERERESnvykXYNsYMNMaMSE5OLu2ueKRmdCX2HNFa2yIiIiLlXbkI29baCdbaIVFRUaXdFY/Uig4l3WXZd+x0aXdFRERERPyoXITtP5ta0ZnL/x3yrpRk5Owt3DFqkT+6JCIiIiJ+oLBdCrI2tvF2RZJ1e44xY8MBdmtDHBEREZE/BYXtUlCzciUcBrZ5ObJtsQDM3njAH90SERERER9T2C4FQU4HNaMrse3gCe/+oTtrM2tjku87JSIiIiI+p7BdSupVCSMxybuwnZm1mbspiQyXLfC5IiIiIlL6FLZLSd0qoSQe9LKMxLoD9pGTaaze9edY5lBERESkIlPYLiX1qoSRfCqNIydTPf43FogMcWIMzP5DddsiIiIiZZ3CdimpVyUMwOvR7eiwIFrWiGK26rZFREREyjyF7VJSL9a9/J83kyStBQP0aBzLsu2HOZaS5qfeiYiIiIgvKGyXklrRoRgDiUmej2xbwBhDj8ZxpLss8zcf9F8HRURERKTYFLZLSUhgADWiKpHo1ci2xQDt60YTGhTADNVti4iIiJRpCtulyL0iiRdhG8C41+nu0TiW6ev3Z69QIiIiIiJlj8J2KaoXG8Y2byZIZtZsA/RrFs+e5BR+33PUL30TERERkeJT2C5F9aqEcuhEKsmnPJvoaLEY447bfZpWxRiYum6/P7soIiIiIsWgsF2K6mYu/+fpiiQ2x8h2XEQwbWpVZup6hW0RERGRskphuxQ1jHOH7c0Hjnv0fGvBmDNf90uoysodR9h/LMUf3RMRERGRYlLYLkV1q4ThdBg27fcwbGMxnEnb/ZrFAzBjvVYlERERESmLFLZLUWCAgzpVQtm834sykhwj282qR1AjKoTf1u3zUw9FREREpDgUtktZo7hwNnlaRpLra2MMfZtVZfbGJFLSMnzfOREREREplnIRto0xA40xI5KTk0u7K15rVDWcxKQTpGW4Cn2ue2TbnPVYv2bxnErLYN7mJH91UURERESKqFyEbWvtBGvtkKioqNLuitcaVQ0n3WU9XG/bYnI90rVhFSKCnUxavdcf3RMRERGRYigXYfvPrFHVcACPJknmrtkGCHYG0K9ZVX5dt8+j0XERERERKTkK26WsYZw7bHuy/J/l3LAN0L9ldY6cTGPhlkM+7p2IiIiIFIfCdikLC3ZSPSqEzR6NbJ+99F+WXk3iqBQYwKQ1e/zRRREREREpIoXtMqBRVc9WJMlvZLtSUAB9EuKYsnYfGa7ca5aIiIiISGlR2C4DGsaFs3n/cawtOCjn3K49t/4tq5N0/DRLtx32fQdFREREpEgUtsuAxvHhnEjNYNeRUwU+z0LeQ9tA34SqBDkdKiURERERKUMUtsuAhGoRAGzYe6zA57lrtvMWHuykZ+NYpqzZW+gIeWEyXJb//LyOnYc9WY5QRERERPKjsF0GNIl3h+31hYRtyHdgG4ABLauzOzmFZduLV0qy/dBJRszawos/ry/W64iIiIhUdArbZUBESCC1YyoVGrYLG7C+qEU8wU4HP67YXaz+uDIb+mn1Htbu/vPtyikiIiJSVihslxFN4yNZv+dogc+xeewgmVNESCAXNItn4qo9pBdjg5ucof6NXzcW+XVEREREKjqF7TKiWfUItiSdICUtI9/nuHeQLChuw2Xn1eDgiVTmbj5Y5L5k1XwnVIvgt3X7WLnjSJFfS0RERKQiU9guIxKqRZLhsgVu217Q0n9ZejeNIyLEyfgVu4rcl6yB7UHd6hEdGsjrv/5R5NcSERERqcgUtsuIhOqeTZIsZGCbYGcAl7SszpQ1ezmVmv8oeUGyarYjQgK5p1dDZv5xgCWJ2gpeRERExFsK22VEvSphBDsdbNibf922u2a7sLFtuPy8GpxIzWDq+n1F6osrs9zbALd1qUtseBD//UWj2yIiIiLeUtguIwIchibxEQWObFtL4XUkQOcGVagaEcz4Iq5KYjMLSYwxhAY5ua93I+ZvOcjcTUlFej0RERGRikphuwxJqBbBuj0FhG08ytoEOAwD29Rgxob9JJ9M87ofWauRODIbu7lzHWpEhfDSpPW4XMXbMEdERESkIlHYLkMSqkeSdPw0B46dzvsJtvCa7SxXnFeTtAzLhFXej25nhe2slU9CAgN4/KKmrN6VzMTV2g5eRERExFMK22VIs8xJkuvyWW/b05ptgJY1I2kaH8G3S3d63Y+sCZKOHE1d0bYmzapH8uqU9ZxOL9rEy7ys2nmEX9bu9dnriYiIiJQlCttlSIsaUQCs3pX3ro3Wi5FtYwzXdqjFyh1H+GNf4dvAn9VO9muceSzAYXh6QAI7Dp3iswXbvXq9goyYtYV7P1vKoq1a7URERETKH4XtMiSqUiB1q4SyJr+wjedhG9yj0U6H4dslO7zqR9bIdu4NdHo2jqV7o1jenraR5FPe14LnJcNlcVkY+tXyItWXi4iIiJRlCttlTMuaUQWMbHteRgIQGx5M34SqfL98F2lebN+eXbOd63Fj3KPbR06m8f7MzR6/XkFc1hIZ4mT/sdM8PW5V9u6VIiIiIuWBwnYZ06pmFDsPn+LwidRzvuftyDbAtR1qk3Q8lRkbDnj8b2x2zfa5jbWsGcWVbWvy8Zyt7D5yyrvO5MFloVZ0KE9c3JRJa/by1WLvRuFFREREyjKF7TKmVU133faa3eeObhdl0Ld30zhiw4P4xotSkrxqtnN6/KImWAuvTdngfYdyt2UtDgcM6dGA7o1ieWHCWjbt967GXERERKSsUtguY1oWMEnSPbLt3dB2YICDq9rVYvr6/SQdz2dJwVyy1tLOa2Qb3CPRd/aoz7jlu1i+/bBX/TmnLetux+EwvH5dG0KDnDz4xXJS0ny34gm4j+mbxTt8VmsuIiIi4gmF7TImKjSQOjH5TJK01ouK7TOubV+LdJdl3DLPlgHMHtku4DkP9GlE1Yhgnp/we7E2uslw2ewLiKqRIbx2bWvW7z3GS5PWF/k18/LH/mM89d0q7vxkMadSfRvkRURERPKjsF0GtcpnkmRRarYBGsdH0L5uNF8s3O5RMM5vNZKcwoOdDOufwModRxi3fJf3ncrRVs71vPsmxDOoWz0+mZfIFB+uv52W7j6mJdsO8+AXy0j3YsKoiIiISFGVi7BtjBlojBmRnJz3Kh5/Ni1rRrHj0CmOnDx7kqS1nm3Xnpdbzq9D4sGTzNt8sPAn59quPT9Xtq1Jm9qVeXnyeo6fTi9Sv6w9t1zl6QEJtK4VxRPfrmT7wZNFet3csi4gLmwez9T1+3lm3GqtfCIiIiJ+Vy7CtrV2grV2SFRUVGl3xSeyJknmHt22WK9rtrMMaFmd6NBAPluwrdDnunJt154fh8Pw/MDmHDh2mnembSpSv3KPbAMEOwN496Z2OIzhvs+X+qR+OyMzWN/cuQ6P9GvMt0t38vLk4k/wzMvEVbt5b8YmhXkREREpH2G7vGlVKwpjYOWOI2c9XpyR7ZDAAK7rUJtf1+1jb3JKgc+1nLtde37a1onmqnbupQATk0543S+XzfsConZMKK9f14a1u4/yj4m/e/26ueVcznDoBY25uXMd3p+5mZGztxT7tXP7dslOXpm8gafGrlK5ioiISAWnsF0GRVUKpFFcOMu25xG2i5q2gRs71SHDZfm6kLWsz4xse/a6T/dPIDDA8K+fvA/F7tVI8v5ev2bx3Ne7IV8s3M73yz2b3JmfrMwb4DAYY/jH5S25pFU1/vXTumK/dm4uawlyOvh26U7u/3yZz1dWyUmTPUVERMo2he0yqm2dyizffvisUoTMaYtFfs16sWH0aBzLl4u2Fzji6skEyZyqRobwUL/G/LZuP7/+vs+rPllr811iEODxC5vQuX4Mz45bwx/7ir7+doYr65jcXwc4DG9cfx5dG1bhyW9XMX3D/iK/dl5tta4ZxfMDm/PL7/sY/MniIte0F2TzgeO0+PtkXvx5XfbxiYiISNmisF1GtasTzeGTaSTmmCBorS3WyDbALefXZe/RFKauLyBc5rNde0Hu7F6fJvHhPP/jWk6meh4sXXlMkMzJGeDg7RvbEhbs5L7PlnKiyBMx3QcVkKOtYGcAH9zanoTqEdz76VIWbPFg8qgHMlwWh8NwR7f6vHF9GxZuPcTNHy7gUB67ghZH0rHTuCx8MGsLd45e7Nc1xA8eP+3xOu0iIiJyhsJ2GdW2TjQAy7advWlMMbM2/RKqUi0ypMCJkq4CtmvPT2CAg39d0YpdR07xv6meT5Z0eXABUTUyhLdvbMvWpBMM+25VkSYeZk2QdOSqWYkICWT0oE7UiQnlzk8WF3uTHnAfU1aov7JtLT64pT3r9x7jug/msye5+FvcZ8k6pqvb1WLOxiSufHcumw8c99nr5/TEtyvp8+oMJq7a7ZfXFxERKa8UtsuoxlXDiQh2snzHmfBX3JptcI8U39S5DrM3JuW7Lbr1smY7S6f6MVzbvhYjZ29hw17PSj4KG9nO0qVhFZ68OIGJq/YwYpb3kxpd2csZnttWlfBgPrurM7ERwdz+8SJ+333U69fPKcNlCcgR6i9oHs+YwZ3Yl5zCNcPns8VHgdiVWQl0fcfafHH3+SSfSuOKd+YyvaC7FkV0LCWdY6fTefCL5TwzbrXfasX3H0vhlpELmbxmj19eX0REpKQpbJdRDofhvDqVWbbtzCRJi8UUe2wbbupchyCng1FzE/P8flFGtrM8c0kzwkOc/N8PazwagbZ5LP2Xn3t7NeDS1tV5efJ6Zv1xwKt+ndmCPu/vx0eG8PldnQkPdnLrRwvZtL/ogTjDnjuC3rlBFb4ccj4paRlc+/78vHcI9VJ6ZtoOcBg61Y/hx4e6U6dKKINHL+b9mZt9uvRgusvSpUEV7u3VkC8Xbefyd+cUq4Y+P2t3H2XOpiTu/WwZQ79aTvJJ/5XGfLpgW7HWiBcREfGEwnYZ1rZ2ZdbvPZpdp+yLkW2A2PBgrjivBt8t23nOxjlwZrv2oogJC+KZAQksSjzE2KWFr/LhKmSCZE7GGF69pjVN4iN46MvlbDvo+VKDWRcQAQUk+1rRoXx2V2eMMdw8ckGRN9Sx1hKQRzMta0bx7b1dCAkM4PoP5nt9wZBb7mOqWbkSY+/tyqWtqvPSpPUM/XqFz1ZCyXBZQgIdPD0ggdGDO3HweCqXvTOHrxZt92moz8hwv9bANjWYuGoPF74x0y8j9QBfLdrO8Bmbuej1mUxd593EXm8kJp1g4Ntz+HzhNr9PZD2aklbkeQ0iIuIfCttlWNu60bgsrNrpHgUt6nbteRncvT4paS6+XHTuMoC2GCPbANe2r037utH85+d1HC5kUqDL5fmqJwChQU5G3NoBY2DIGM8nTGa4PDumBnHhfHZXJ06nu7hp5IIi1VjnLiPJ/frj7u9K7ZhQBn+ymHHLir7sYPZyhjmOqVJQAG/f2JYnL27Kjyt3c+37vqkTdx+T++2iV5M4Jj3Sg/Z1o3l63Goe+nI5x1J8MwKdnvl7urdXA354oBuVQwMZ9Mliho1d5bM2smS4LAnVIggPcXLn6CU88Pky9h8teA36oli7+yirdyXz1+/XcOV7c1mRa/18Xxo8ajHdXp7GJ3O3kubHNd4PHj/Ni5PW5VuK5ksnU9M5na4lLkXkz0thuwxrW7syAMsyJ+35cgQxoVok3RpVYfS8xHM+lItas53F4TD8+8qWHEtJ55+FrL2d1w6ShalTJZR3bmzHxv3HeHLsSo9+LgXVbOeWUC2SMYM7ceRkGjd/uJADx7xbhSPDVfBOn/GRIXxzbxc61Y/hsW9W8u70ou02mZGjjCQnYwwP9GnEh7d2YGvSCQa+PZdFWw95/fpnt2UJyPFuUTUyhDGDO/PkxU2ZtGYvl/5vDqt2Fj9EZl0UOR0OWtaMYsJD3bmvd0O+XbqD/m/OZu6mpGK3kSXdZWkYF87Eh3rw5MVN+XXdPvq9PpMvFm7PLjvyTTvu39NjFzZhb3IKV743l6e/W+Xz1WkADhw/zcnTGTw/4Xf6vznLp0ta5jRjwwE+mLmFi96YxTPjVvvlIiXL9R8soPvL0xkzP5HUdP9dQGzaf4z7P1/K9A37/b776+LEQz4pJSuMy2V9fpEqIt5T2C7DKocG0SAu7KwVSXxRs53lzu712Xs0hUlr9p71uDfBND8J1SK5r3dDxi3bVeAHvvVwgmRu3RvH8syAZvy8ei/vzdhc6POz69A9PONb16rMqEEd2ZPsnrDnzbJ3OVcjyU9kSCCfDOrEZW1q8OqUDTw3fq3XJQY5N+rJywXN4/n+/q5EhDi56cMFjJq7tcghIsNanLl+eAEOd6j/esj5pGe4uHr4PEbO3lKsoJqRqzQm2BnAsP4JjL2vK8FOBzePXMhz49d4tbxkvm1l3oEIcjp4oE8jpgztScsaUTz7/WquHzHfZ6O2Wb/Xy8+rwdTHe3FX9/qMXbqTPq/N4NMFvi0tSc+w/KVNdUbe1gGXhUGjFnP7x4vY6OP6+qwL9Mva1GDs0h30enUGr//6h1/q3/ckp3D0VBrPjV9Lv9dn8P3ynX4px5m3+SA/r97LoFGLuXr4POZuSvJb6L7vs6X85e053PnJYp9cpOZn/MpdtH7hFx74fBnr9hRv4ndh7h6zhAe/WFbsCeaFWbs7mWe/X+3XO0RZfly5m3mb/XceZDlyMpX5mw+WyH4Ju4+c8utGa1mstX69OM7djr9/R8WlsF3Gdawbw+LEQ7hcmZuo+y5r07tJVRrEhvHxnK1nPe7Ndu0FebBvIxpVDeev41bn+yHssvmXXBTmrh71ufy8Grz2ywamrS+45ja7vtmLYN+xXgwf3d6BbYdOcPOHCznoYeAuqIwkpyCngzevP497ejbg0wXbuP/zpV69CZ4Jpvk/p3F8BOMf7EbvplV5YcLvDP16RZGCatba4XnpUC+Gnx/pQZ+mVfnXT+u4c/Rij39W57bjfnN25mqrXZ1ofnq4B4O61WPM/G0MeGs2SxKLN1qf7nLhzFFcXz82jC/u7syr17Rm4/7jDHhrNm/8+kexSxjSM+vQnQEOIkIC+eulzZn0SA+aV4/k/35Yw+XvzmHptuIvOQnuYwp0OLigeTxThvbkb5c2Y9n2w/R/azbPjV/js9H0rHKfZy9txm+P9aJfs6r8b+pGer86nU/nn3u3rHhtubiuQ21GDepIRHAgj369kkvems2vv+/z6QdsWubv6dlLEtiTnMLNIxdyw4gFxb4rlJeTqRk0iQ9n6fbDXPbOXAZ/spiVfgiPu4+kYC3M2LCfAW/N5u4xS1i90z8j6lPX7WPiqj1c8r/ZDP5ksc/O6dymrN3HFwu3c8W7c7npwwXM2eifMJya7uKRr5Zz04cLueLdufy8eo/fwvDHcxO58cMF9PvvDD5bsM1vYfhUagZ9/zuDbi9N463fNvrl7lqWd6dvotXzU3hm3OpiLThQmP3HTtPkb5P4YtF2v7XhCwrbZVznBjEcTUlnw75jYH2atXE4DIO61WPFjiMs3XbmA8Xb7drzE+wM4OWrW7PnaAqvTF6f53MyirFRjzGGl65qTYsakTz85QrW781/ROXMDpLeNda1USwf396RbYdOcJOHgduVx2ok+XE4DM9c0ozn/uLebfLmkQsLrXPPbseVFbYL/jOODAlkxK3teeKiJvy4cjdXvTePxCTPJ5eC++eXOwDnVDk0iA9ubc8/Lm/B3E0HufjN2cwoQglDVjDN62KlUlAAfx/Ygi/vPp8Ml+XaD+bz759+L/IHU0bGucdkjOHaDrX57bFeXNqqOm9N3ciAt2azsBgbHqVnl8acaatxfARf3N2Zt29sy4Fjp7l6+Dye/HZlsTcOSs+w2RcQQU4Hd/VowMwn+3Bz5zp8vnA7vV+dzkdzthZ7xClrB1qnw0HdKmG8c1M7fnigGw3iwvm/8Wu5+I1ZTF6z1ychKCPzmPo0rcrEh7rz9o1tSc1wcfeYJVw9fB7zN/tmM6qsY7rl/LpMf6I3fx/YnM0HTnDdB/O59aOFPh1JTc+w9E2IZ/ZTfXjy4qYs236Yy9/1fejO+nuaM6wvj17QhEVbDzHwnTncMWqRT8Owy2VxWffd0scvbMLy7Ye5evg8bhgx3+dhOC3DhdNh+Oslzdi0/zi3fLSQy97xfRhOy3BhLXRvFEvyqTTu/3wZF2SWmfk6DB9LSSMowEFUpUD+9sMaur00jf9N3ejxZ4GnTqamk5LmIiQwgDd++4MuL07lr9+v9tlytDntPHyKDJflu2U7ueD1mQz+ZLFf7hJkvZcFFjTqVAaU7d4JnerHALBwy8HMCZK+jNtwdftaVA4N5P2ZZ9auzvpj8EVb7etGc0dX92hkXiNERS0jyVIpKICRt3UkLDiAOz9Zkm99ddbfd1FG0bs2iuWj2zuSePAEN48sPHBnuPJejaQgg7vX550b27F6ZzLXvD+PHYcKXwklK8R5MlrvcBge7NuYTwZ1Yu/RFAa+M8erFTg8Ga03xnBbl3qMf7AbMWGB3DFqMS9MWOvdaL0r/7CdpUvDKkwe2pMbO9Xhw9lb+cvbc4q0GVFajkmfucWGB/PmDW0ZPbgTqekurh+xgKe/W1WkpQgLqq0f2KYG0x7vzT29GvD98l30fW0Go+clZgc/b6XncVEUExbEPy5vyaRHetCmdmX+OdFdzz11XdFHhrMvIHKc6OfVrszXQ87no9s74HAY7v1sKde8P/+sC/miSHO5sj9IHQ73z+yXR3vy4lWt2H0khRs/XMCtHy0s9ojtmYsiByGBAQzqVp/ZT/Xh2UsSWLv7KFe8O5c7P1lc7Fpray2pGS6CAgwRIYE80KcRc4b1PSt0Dxq1yCfhPt3lwhiIDgvikQsaM2dYH57q35RVO5O5evg8bvpwAfM3Hyx2AErLPMdjwoJ4qF9j5gzry98ubcbWpBPc8tFCrnhvHr+s3euTuRDpGS6CnA7u7tmA2cP68NJVrTiW4g7DF74+k68Xb/dJ+ULWhUrfhKpMfbw3793cjogQJ89+v5ruL0/n3embfLZrb1qGi4gQJz880I2vhpxPm9qVef3XP+j60jT+Pn6NR58HnrXjPqYH+zbit8d6clW7mny7dCf9Xp/JkDFLWJJ4yGdhOC3DEh8Zwryn+zL0gsas3HGEmz5cyF/ensMPy3f57O5X1t9tkMK2FEet6FBqVq7Eosw/At9GbffqHnd0rcevv+/LruvMniDpozaevLgptaIrMey7VecEr6JMkMytWlQII2/ryKETqQz5dEme4e7MaiRFa6Nbo1g+vqMjW5Pcgbug228FlVwU5NLW1RlzZyf2HzvNle/NK3SEK3vtcC/+ins1iWPCg93dO2aOXsLrv/7h0QdgusvlcQlOs+qR/Phgd+7oWo9RcxO5/J25Bd51OLudc0eB8xIe7OQ/V7Zi9OBOnDidztXD5/Hvn373arOdwkbrwf3z+uXRntzTs0H2h9KElbu9+kDK+oALzOcXFRbs5JkBzZg8tCeta1Xm7z+uZeA7c4tUJpOe4cKZz4dOk/gIxgzuxKg7OoKBO0cv4daPFnn8u8kpv2MyxtCvWTyTH+nBS1e1Ysehk1w9fD73fLqkyLubpudxByIwwMGNneow48nePHtJAqt3JTPwnTk88PmyIreT9eEfmOMColJQAEN6NmRW5gj04sRD/OXtOdz76VKPN+7KLXsScI7fU3iw86zQvXzHEa7wQehOzXCdNeIXERLI/b0bMWdYH/52aTM27j/OjR8u4LrMpUiLGrSyzoeswBMW7OSuHg2Y9VQf/nNlKw6fSGXIp0sZ8NZsxq/YVeSLyay2so4p2BnADZ3qMPXx3rx7UzsqBQUw7LvV9HxlOiNnbynWMpipOc6HAIfhklbVGf9AN764qzPNa0Ty6pQNdH1xKv/+6Xf2JhdvgnBauvvujTGG8xtU4eM7OvLLoz25tHV1vli0Q315ygAAHnZJREFUnV6vTueBL5YVu77/zDnuoFHVCF68qjVzh/XloT6NWJR4iGven89Vw+cxyQd3CdIyXAQGGGLDgxl6QRPmPt2XF69qRUpaBkO/XkGPl6fz/szNxb5gyXlMZVnZ7p0A7tHtRVsP4fLROtu53d6lHpUCAxg+0z3R8EzNtm8aCw1y8tJVrdmadII3f9t41ve8WWe7IK1qRfHG9eexfPsRnhx77pbuxdmoJ0vOwH3ThwvyDdyeTJDMz/kNqjDuvq6EBDq4fsT8AndSzD2Z0FO1Y0L57r6uXNO+Fv+bupHBoxfnud76WW25IMCL4fqQwACev6wFowZ15OCJVC57Zy4fz9laaLD3ZGQ7p6wwfEPmKPeAt2Z5XPKRnuHyqJ3QICfPXNKM8Q90o3pUCA99uZzbRy32eJ337GMq5OfXqGo4n97ZifdubseRk6lc8/58HvtmhVer4eQ1sp2TMYY+CVWZMrQnfx/YnNW7krnkrdn89fvVXtXZZ9fW53NMzgAHN2SG4ccvbMKcjUlc9MYs/vbDaq+Ox1rrPqZ8PkhDAs+E4Yf7NmL6hv1c9MYsho1dxe4j3i15mZ7hvnuT1x297DD8dF8e6deYOZuS6P/WLB76crnX4T6vuwLntDOsL0/1b8qKzNB9x6hFRbp7k55hCczjfAgNcofh2U/14R+Xt2DX4VPc9vEirnhvXpHueGSXFeU6pmBnADd1rsO0x3vx5vXnYbE88tUK+r0+ky8XbS/SfIisEJdTgMNwaevqTHyoO2MGd6JebCj/+mkdXV+axuu//lGk2uS8Qpwxhq6NYhkzuBM/PdydC5rH8/HcRHq8Mo0nv11Z5EnVOe/eZGkSH8Fr17Zh9lN9ubtnA2ZtOMBl78zlxhELirxiTl4XlHERwTx2UVPmPd2Xf17egkMnUrnv82X0eW0GY+YnFnkyelquC72QwABu7FSHXx/txag7OtIgLoyXJq2n64tTeWHC2iKP3p8pI/FDOPIhhe0/gc71Y0g6nsr2Qyd9PrIN7luMN3aqw48rdrPz8MnsbcB9Gey7N47l+g61+XD2lrNGarxdZ7sg/VtWY1j/BCas3M1bU88N9eB5LXV+umWWlBQ0wl2cSZ/gruf9/v5u7hVdPl/GiFl57waZ7mUwzSkkMIBXr2nNv65oydxNSQx8Zw5rd+d/ezzD5Sp0FDgvfZpWZfLQHvRoFMs/Jv7OHZ8sLnCZuJy38j0VERLIf65sxRd3d8Zl4foRC3hu/JpCV8bwZGQ7p5Y1o/jhgW78fWBzlm07zIVvzOKt3zYWGhjS8pn0mRdj3CNoUx/vxf29GzJh5W76vjaDkbO3eHRr3B1MC28nMMDBoG71mflkb27rUo+vFu+g96szGDFrs0cBKGsks7BjCg1y8lC/xsx8yl03/tWiHfR+dTpv/bbRo1HHrPMhr8CYU2RIII9d1JRZT/Xhti51+X75Lnq/NoN/Tvzd44uINA/O8ciQQB69sAmzn+rDvb0a8tvv+7jw9Zk8/s1KjzfByg48BZzj4cFO7u/diNmZoXvljiNc+d48r0N3QXc6wP0+cFuXesx4sg8vXtWKQydOc+foJVz6vzlMWr3H47KP1OywnXdbzgAHV7StyeRHevLBre2JqhTIM+NW0+uVGXw0Z6tXgS53iMvJGEPPJnF8NaQL4+7vSqf6Mfxv6ka6vTSNFyas9eoCLKuMJL+2WtSI4q0b2jLjid7c1KkOE1bt5oLXZ3H3mCVe18OnZdh8yyCqRYXwzIBmzHumL89eksDWpBMMGrWY/m/OZuzSnV6VzKQVcEyhQU5u7VKPaY/35v1b2lElPIjnxq+l60vTeG3KBvYf8270Pr/fk8PhvuD/4u7zmfhQdy5qUY1P52+j16vTuf/zpdlLHXvTTn7HVJaU7d4JcKZuG3xfs53lrh71ARg5e2v2DpK+GtnO8uylzYiPCOaxb87sbOjNdu2euLdXA65tX4s3f9vI+BW7sh/P+swo6ohzTt0bxzLy9g5sOXCcmz5ccM6ktgwfXEDERQTz1ZDzuaRldf7z83qe/X7NOTVu2RMki9iWMYZbzq/L1/d0IS3dctV78/h6cd47Qma4in4HIjY8mJG3d+CfV7Rk0daD9M9cSSIv2fXNRRil6NowlslDezC4W30+XbCNi9+YxeyN+e/Sme6yXrcT4DAM6lafqY/34sLm8bzx2x8MKGT97wwPg2lOoUFOnuqfwJShPWlfL5p//bSO/m/NKnDSqbX2rM2HPFE5NIjnL2vBlKE96FAvmv/8vJ6L3pjFpNV7Chw5S3e58h0FzktseDD/uLwlvzzak55N4njjtz/o/doMvli4vcCSgpwruXjazt8HtmDaE724rE0NRs3dSs9XpvPGr38Uut50eo7yhMJEhwUxrH8Cs4f1YXC3+kxctZu+/53BM+NWsauQQJd9oeLBuZczdA/rn5Adum//eJFHoSTN5dkxBTndZTnTHu/Na9e2ISUtg/s+X0b/t2YxfsWuQksK0rPLSAo+JofDcHGLaox/oBuf3ukegf7nxN/p/vJ03pm20aOSgpyTgAvSrk40H97WgV8f7cmAVu5A1/OV6Tzx7UqPVsdIzWe0PrfaMaG8cHlL5g5z3/VYnHiIq4fP49r33XcJPLlgSUt3FdpOREhg9l2c/17bBoAnvl1Jz1emM2LWZo/WU/ckmAY4DP1bVuf7+7vx3X1d6Fw/hndnbKL7S9MZNnaVx8uIpmVYAp0Fn3sta7rvSM8Z1pchPRsyZ2MSV703j6uHz2PyGs9KWQq6gChLynbvBHAvSRYbHgz4djWSnGpUrsQVbWvy1eLtRV62rTBRlQJ59do2bDlwgpczVydxFXOCZG7GGP59ZSs61Y/hybGrskcYiluznVuPxnF8fId70uT1H8xnX47RWvfIdvHbCAl07wZ5f++GfLloO4M/WczRHG+o3pZc5KddnWgmPtydDvWiGfbdap74dtU5I03ejgLnZozh1vPrMvGh7lSLDOHuMUv46/erz6mxzspdRW0rNMjJcwObM/beLgQHOrj1o0U8NXZlnh/ixTmm+MgQ3r2pHWMGdyLDWm4euZBHvlqe5+hPce5ANIgL55NBnfj4jg5YC3eMWsxdoxfnuZqMp6PAeWlUNYJRgzoxenAngp0O7vt8Gdd9MD/feuG86qg90SAunOG3tOe7+7pSNyaUZ79fzcVvzmLK2rxXLsm6K+DtLeJa0aG8dm0bpgztSY/Gcbw1dSO9XnXfIchvwq57FNi7dmLDg/nbX5ozK3Pk/rulu+jz6gyeG7/mrPeE3O24j8nzN4nwYCf39W7InMzQvWrnEa7yIHSnpZ9bclGQwAAH17Svxa+P9eKtG84D4JGvVnDh6zMZu3RnvpPasi+KPLzQM8bQo7F7BPq7+7pwXu3KvPbLH3R/aRqvTF5f4Ko8uevQC9M4PoLXrzuPGU/25pbz6zJx1W4ufGMm93y6pMB5MVmbUXk68a5KeDCPXtiEeU/35e8Dm7P7SAp3jl5C/7dmFToCnZ5HGUl+gpwOrm5fi8lDezBqUEfqx4bxn5/X0/XFabz487oC68dT8ygjKUj7ujF8cGsHpj3em+s61uKHFbu48I1ZDP5kcaETa9MyJwF7olpUCE8PSGD+M/14fmBz9h9L4d7P3KUsn8zdWuBdsLxKY8oihe0/AWMMnbNGt/14Pt3bqwEpaS5GzU0Eil9ykZdujWKzJ87N25Tkrtn28VkY5HTwwS3tqRHlDnWJSSfObEHvw2Pq1iiWMYM7szc5hes+mM/Ow+7byO7VSHzTjsNheKp/Aq9c05r5mw9y9XtnVirxVdgGd2gYM7gzj/RrzLjlO7ni3blnjf5kFLM0JkujqhF8/0BX7unZgC8WbefSt2eftYpEfit3eKt93Rh+frgH9/VuyHfLdnHRGzP5LcdoenYtcDFPvp5N4pgytCcP92vMpNV76fffmYyZn3jWiEx6ZnlCce529E2IZ/LQHu4PpM0HueiNWbw0af1ZpTJ5TbzzVq8mcfz8cA/+c2Urtiad4Ip35/Lwl8uzz+0saV6MAuelfd1ovr23CyNubY8F7vl0KVcPn3fOikXpRbgrkFPj+Ajev7U94x/oRvPqkfzrp3X0fW0GXy3afk5wTM0o+vkQHxnCC5e3ZPqTvbm6fS2+WLidnq9M558Tfz8nOKZlXRQVIRyE5QjdTw9wTwy96r153PZx3kv5eVpWlFuAw3D5ee6yj/dvaUdIYABPfLuSvv+dwZeLzl3tw9NR4Ly0rxvDx3d05KeHu9OzaRzDZ26m+8vTeP7HvMs+3HXo3v+eakWH8vxlLbInBM7ffJDLC1irOy29aH9PoUFOBnWrz4wne/Pm9efhMIYnvl1Jr1fzn7SZmpH/vIT8GONeDvPLIefz44Pd6NU0jg9nb6HHK9N4/JuVeU7gLaw0Jj/1Y8P41xWtmP9MPx67sAkrdxzhxg8XcNk7c/lx5e48704VVO6Tn7BgJ3d0q8+MJ/ow/OZ2xIYH8fyE3+ny4lRenrw+zwuJ7AuIQkbRS1vZ7p1k69zAHbZ3HvJu0o83GlWN4KLm8ezNHJHxV64f1j+BBrFhPDl2FafSMvxSGhMdFsSoQZ2w1nLHqEUcOO6urfZ1aUyn+jF8dldnDp9I5foPFrDt4AlcRVyNpCDXdajNmMGd2HfUvd338u2HizxBMj8BDsOjFzZh9KBOJB1P5bJ35mSX4ni6UY8ngp0BPHNJMz6/szMnT2dw5XtzeWfaRtIzXF4tZ1iYkED37pPf39+V6NAg7hqzhEe+Ws6hE6k5toX3TTuPXdiEyUN70KZWZZ4bv5Yr35ubfRGR7qOfXbAzgHt7NWT6E70Z2KYG78/cTN/XZjBu2U5cLpsdHot7TM4ABzd1rsOMJ/vwYJ9GTFm7l77/nclLk9Zn31nJcHk2ubQgxhgualGNX4b25KWrWrHryCmu+2A+gz9ZnL1CSnohtcCealO7Mp/d1Zkv7upMXGQIT49bzYWvz2T8il3Zt/nT85h4562alSvx4lWtmPa4+3c0au5Werw8nZcmrc+e35FzjfKiCgt2cm+vhsx+qg9PD0hgzS73Un65Q3dRAk9OjsySgp8e7s5Ht3cgJizYXWv96nRGz0vMvkvg7ShwXlrUiOLdm9rx22O9GNi6Bp8tcNfxDhu7iq057uSkFeEORE5VwjMnBD7Tj79e0ozNB/Jeq7uod1WyBGbWqU96xD0CXbfKmUmbr03ZcNZFWFq656PAeWldqzLv3NSOGU/04aZOdfhp9W4ufnMWg0YtOmsEurj1zTFhQTzcr3H2yiInUtN5+Mvl2XeNcl78pxbjgjzAYRjQqjrj7u/Gd/d1pXvjWD6YuZker0zjsa9XnLVLaXqulXDKqoDnn3++tPvgMyNGjHh+yJAhpd0NvwgNCuDTBdvYdeQUQy9o4rd26sSE8tXiHQDc3aMB4SFOn7cRGOCgda0oPpqzldPpLtrUjqJvQrzP24kODaJjvRhGz0vMXs3l/t6NCAkM8Gk71aMq0aNxHN8s2cF3y3Zx/HQ6bWtH07NJnE/bqR0TyoXNq/HT6j2MnufeZWzXkVM81LdxscNITnWrhHH5eTVZnHiIj+cmknT8NCt3JtOlQRW6NIz1WTu1Y0K5tn1tdh4+xah5iczZlITLwh/7jvHohb47x+MjQ7iuQ20CHIbPF2znmyU7iI8MYcraffRoEkvHejGFv4gHosOCuLJtTRrEhfHTqr2MmruVQydSOZ6Swab9x3mgTyOftBMW7OTiFtXo2SSOpdsOM2b+NmZtPEDdKmH8sHwXfROq0rZOdLHbCXI66Noolqvb1SLp2Gk+XbCNrxfvoFJQAHuOpLDt0AmG9GxY7HYcDkPLmlHc0rku4cFOxq/Yxai5iew4dIraMaF8u2QnF7eIp2XNqGK3VTsmlBs61qZlzSgWJx7iswXbmbxmL/GRwWzYe4x9x1IY3K1+sduJCg3kohbVGNimBknHT/PFou18Nn8bJ1PTiY8M4ZslO7m0dXUSqkUWq50gp4MO9WK45fy6RFUKZNKavYyel8jSbYepExPKksTDHD2Vzq1d6harHWMMDeLCuaFjbdrXjWb93qN8vnA7Xy/ZgcO4a/+/XbqTK9rWpGFceLHaigkL4qIW1biqXU3S0l2MXbaTUXO3sunACerHhjFzYxLpGZYbOtUpVjtBTgft60ZzW5e61KpciflbDvL5wu1MXLmbkEAHlQKdjFu+i+s61KZ2TGiR2zHGUD82jGva16ZXkzj2HT3Nl4u3M3peIvuOptAwLpxJa/Zml4cUR1RoIH0SqnJz57qEBTv5Ze0+Pl2wjekb9hMR4sTg3oL+ps51qB5VqcjtOAMctKoZxa3n16VVzSg27j/Ol4t28OmCbSSfSqNx1Qh+WL6LqEqBXHZezWIdU43Klbi0dQ2ubFuLDJdl/MrdfDIvkcWJh4gJDeJUWgY/r97LbV3rZpfblqQXXnhhz/PPPz+isOeZsr6fvDc6dOhglyxZUtrd8AtrLfWf+Zlujarw+V3n+7Wtek//BMDCZ/sRHxnit3b++8sG3p62idu61OUfl7f0WzuT1+zhvs+XYS2sev4iIkMC/dLOhr3HuHnkQpKOn+aeng145pJmfmnn4PHT3PvZUhYnukewNv17gE/Ddpa0DBevTtnAiFnuDY+GXtDYbxd641fs4v9+WMPRlHQCAwwb/32JX9pZt+coT41dxerMjUmevSTBJ4Ext+RTafz3lw18umAb1rrnK6z8+0U+b8flsoxdtpNXJp8ZKfvn5S24tUs9n7e1emcy//rpdxZmlnpUiwxhwbP9fN7OkZOpDJ+xmVGZm/u4LPz32jbFDiK5uVyWCat28+ZvG9madIIAh6FuTCjTnujt03YANu47xptTN/LTqj0EBhjSMizDb27HgFbVfdrOidPpfLZgGyNmbeHgiVSCAhw0jg/np4d7+LQday0Lthzi7WkbmZdjB89RgzrSp2lVn7Z14NhpPpqzlc8WbOP46XSCne7Bmm/v7erTdjJclslr9vLejE2s3X2USoEBnErL4Jt7upy1SIEvbD5wnJGzt/Dd0l3Z9dqdG1RhzOBOPm0nJS2D75bt5MNZW0g8eJKwoABOpGYw8aHuPrl4zWnFjiN8OHsLk1bvwWEMDoehX0JVht/S3qftJJ9K48tF2/lkbiJ7j6YQEezk2Ol0pj7eq9gXekVhjFlqre1Q2PM0sv0nkTXB7Op2tf0+6/YvratjgEtaVffb6icAHerGsHHfcXo3rUrTahF+a6dR1QiqhAVx4Phpbj2/rl9q0cFd99yvWVXmbEzi4hbVaFHDt29mWUKDnFx+Xg0OHD3N8dPpDO5W3y+/pwCHexJT8+qRzNmUxGVtatCsevFG4vKTUC2SK9rWZOP+44QGBXBz5+KNxOUnLiKY6zrUIjzEyYrtR7iuY20aVfX9G3RIYAB9EqrSN6Eqa3cfJTo0kBuLORKXF2MMLWpEcWOn2mS4LKt3JnNrl3rUqxLm87biI0O4pn0tWtSIZO3uo1SPqsS1HWr7vJ2QwAB6NI7j6na1OJqSxro9R7mzewNqRhd9JC4vxhgSqkVyS+c61IyuxJpdydSLC+fKtsUbictLlfBgLm1Vnf4tq3Hg2Gk2HzjOkJ4NfD6YkTXSfWuXukSHBrF611ESqkUwsE0Nn7ZjjKF2TChXt69Fj8axbNx/nL3JKTzQpxFVfDy6GBbspHvjWG7uXJfQoADW7E6mTa3K9G9ZzaftOIyhSXwEN3WqQ/u60ew8fIrdmXcOoyr5doAmJiyIC5rFc33H2jgcDn7ffZSO9WPo18y3d3idAQ5a16rMrV3q0bx6BIkHT3Lw+Gke6tuY8GDf3rWuFhXCpa2qc1U79wj0uj1H6dYw1ud3eEMCA+hQL4bbutSjYVw4mw8c51hKGg/3bezzu9aeqJgj2/Xr2yW3317a3RApMdZav14QlXQ78ueQ9blRns69kjwmX86BKIgrc9fhcnVMLosxJXNM7snN/q8Fdi+t6v9jyloe1B93QnMrj8eUF/PCCx6NbJevsF2Oy0hEREREpOzwtIykbE/fFBERERH5E1PYFhERERHxE4VtERERERE/UdgWEREREfEThW0RERERET9R2BYRERER8ROFbRERERERP1HYFhERERHxE4VtERERERE/UdgWEREREfGTchG2jTEDjTEjkpOTS7srIiIiIiLZykXYttZOsNYOiYqKKu2uiIiIiIhkKxdhW0RERESkLFLYFhERERHxE4VtERERERE/UdgWEREREfEThW0RERERET9R2BYRERER8ROFbRERERERP1HYFhERERHxE2OtLe0++Iwx5gCwrZSajwWSSqltKXt0PkhOOh8kN50TkpPOhz+nutbauMKeVK7Cdmkyxiyx1nYo7X5I2aDzQXLS+SC56ZyQnHQ+lG8qIxERERER8ROFbRERERERP1HY9p0Rpd0BKVN0PkhOOh8kN50TkpPOh3JMNdsiIiIiIn6ikW0RERERET9R2PaCMaayMWasMWa9MWadMaZLru/3NsYkG2NWZP7vudLqq/ifMaZpjt/1CmPMUWPM0FzPMcaY/xljNhljVhlj2pVWf8W/PDwf9B5RwRhjHjXGrDXGrDHGfGmMCcn1/WBjzNeZ7xELjTH1SqenUhI8OB/uMMYcyPEecVdp9VV8x1naHfiTeQuYbK29xhgTBITm8ZzZ1tq/lHC/pBRYazcA5wEYYwKAXcD3uZ42AGic+b/OwPDM/5dyxsPzAfQeUWEYY2oCDwPNrbWnjDHfADcAn+R42p3AYWttI2PMDcDLwPUl3lnxOw/PB4CvrbUPlnT/xH80su0hY0wU0BP4CMBam2qtPVK6vZIypB+w2Vqbe1Oly4Ex1m0BUNkYU73kuyclLL/zQSoeJ1DJGOPEPUCzO9f3LwdGZ/73WKCfMcaUYP+kZBV2Pkg5pLDtufrAAWCUMWa5MWakMSYsj+d1McasNMZMMsa0KOE+Sum5Afgyj8drAjtyfL0z8zEp3/I7H0DvERWGtXYX8BqwHdgDJFtrf8n1tOz3CGttOpAMVCnJfkrJ8PB8ALg6s+xwrDGmdol2UvxCYdtzTqAdMNxa2xY4ATyd6znLcG/d2QZ4G/ihZLsopSGzpOgy4NvS7ouUvkLOB71HVCDGmGjcI9f1gRpAmDHmltLtlZQWD8+HCUA9a21r4FfO3PWQPzGFbc/tBHZaaxdmfj0Wd/jOZq09aq09nvnfPwOBxpjYku2mlIIBwDJr7b48vrcLyDkyUSvzMSm/8j0f9B5R4VwAbLXWHrDWpgHjgK65npP9HpFZWhAFHCzRXkpJKfR8sNYetNaezvxyJNC+hPsofqCw7SFr7V5ghzGmaeZD/YDfcz7HGFMtq9bOGNMJ989Xb5rl343kXzLwI3Bb5qok5+O+bbin5LompSDf80HvERXOduB8Y0xo5u+9H7Au13N+BG7P/O9rgGlWG2CUV4WeD7nm9FyW+/vy56TVSLzzEPB55m3iLcAgY8y9ANba93G/Ud5njEkHTgE36E2zfMus278QuCfHYznPiZ+BS4BNwElgUCl0U0qIB+eD3iMqEGvtQmPMWNzlQ+nAcmCEMeYfwBJr7Y+4J91/aozZBBzCXe8v5ZCH58PDxpjLMr9/CLijtPorvqMdJEVERERE/ERlJCIiIiIifqKwLSIiIiLiJwrbIiIiIiJ+orAtIiIiIuInCtsiIiIiIn6isC0iIiIi4icK2yIif3LGGGuM+SzH105jzAFjzMTS7JeIiChsi4iUByeAlsaYSplfX4h7G3ARESllCtsiIuXDz8Clmf+d75bxWYwxzxtjRhtjZhtjthljrjLGvGKMWW2MmWyMCcx83nPGmMXGmDXGmBHGzZn5WO/M57xojPm3Pw9OROTPSmFbRKR8+Aq4wRgTArQGFnrwbxoCfYHLgM+A6dbaVri3ks8K7u9Yaztaa1sClYC/WGvTcW8jPdwYcwHQH3jBlwcjIlJeKGyLiJQD1tpVQD3co9o/e/jPJllr04DVQAAwOfPx1ZmvBdDHGLPQGLMadzBvkdneWuBTYCIw2Fqb6oPDEBEpd5yl3QEREfGZH4HXgN5AFQ+efxrAWusyxqRZa23m4y7AmTlK/h7QwVq7wxjzPBCS49+3Ao4AVX3TfRGR8kcj2yIi5cfHwAvW2tU+er2sYJ1kjAkHrsn6hjHmKiAG6Am8bYyp7KM2RUTKFYVtEZFywlq701r7Px++3hHgQ2ANMAVYDGCMiQVeAu6y1v4BvAO85at2RUTKE3PmrqGIiIiIiPiSRrZFRERERPxEEyRFRMoxY8wg4JFcD8+11j5QGv0REaloVEYiIiIiIuInKiMREREREfEThW0RERERET9R2BYRERER8ROFbRERERERP1HYFhERERHxk/8H0Cg5o/aw/hUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110020898>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "_M_min = 6.5\n",
    "_M_max = 8.8\n",
    "_M_max_step = 0.001\n",
    "_bin_size = 0.1\n",
    "\n",
    "# fault params from an example in NE Asia\n",
    "fault_area = 317.11\n",
    "fault_slip_rate = 0.225\n",
    "fault_b_gr = 1.\n",
    "\n",
    "# test max ranges\n",
    "muppets = numpy.arange(_M_min, _M_max, _M_max_step)[1:]\n",
    "\n",
    "annual_mos = [calc_moment_release_from_fault_params(fault_area, fault_slip_rate,\n",
    "                                                    _M_min, M_upp, fault_b_gr,\n",
    "                                                    bin_width=_bin_size)\n",
    "              for M_upp in muppets]\n",
    "\n",
    "\n",
    "annual_mo_accumulation_rate = annual_mo_from_fault_params(fault_area, fault_slip_rate)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(12,9))\n",
    "plt.plot(muppets, annual_mos, label='OQ-MBT moment release rate')\n",
    "plt.axhline(annual_mo_accumulation_rate, color='r', lw=0.5,\n",
    "            label='Moment accumulation rate')\n",
    "\n",
    "plt.gca().set_yscale('log')\n",
    "\n",
    "plt.xlabel('M_max')\n",
    "plt.ylabel('Annual moment release rate (N m / yr)')\n",
    "plt.title('Moment release rate as a function of M_max\\n(M_min={0}, bin_size={1})'.format(\n",
    "                                                                        _M_min, _bin_size))\n",
    "plt.legend(loc='upper right')\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the annual moment release rate is only equal to the accumulation rate on rare occurrences. The problem is most severe for small faults, with `M_max` values very close to `M_min`. \n",
    "\n",
    "Please note this is a log scale.  For very small faults, with `M_max` <0.05 magnitude units above `M_min`, the annual moment release rate is 1-2 *orders of magnitude* higher than the accumulation rate.  This causes bullseyes around very small and otherwise forgettable faults in the hazard models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A correction\n",
    "\n",
    "My proposed correction is to round the `M_max` up to the next multiple of `bin_size` above `M_min`.  For an evenly-distributed MFD, it doesn't make sense to have the upper end of the range not at the end of a bin.  This creates big problems, as seen above--The bin midpoints need to actually be in the center of the bins, and evenly-spaced, otherwise we need to use a different MFD altogether.\n",
    "\n",
    "The most basic objection here is that we are 'forcing' earthquakes on a fault bigger than the `M_max`. This is true, but keep in mind 2 points:\n",
    "   1. There is no empirical or theoretical value for maximum earthquake magnitude on a fault.\n",
    "   2. We use rupture area-magnitude scaling relationships to infer `M_max` based on a median full-fault rupture. There is a *huge* amount of scatter in this data. The scatter is about an order of magnitude, while the bin sizes are typically 0.1 orders of magnitude. We can round up and still be *well within* the uncertainty in the scaling regressions.\n",
    "   3. Rounding up lets us __a)__ calculate correct rates (which is really important) and __b)__ accommodate slightly-larger but still realistic values into the hazard model, which is a good, conservative practice.\n",
    "   \n",
    "   \n",
    "Another objection is that the bin size in a model may be very coarse (say, 0.5 M), and this can lead to forcing much larger earthquakes.  Well, garbage in equals garbage out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This correction requires three things:\n",
    "  1. A function to round the magnitude up.\n",
    "  2. A function to generate the bins correctly. There is a known bug in `numpy.arange`, which is listed in the documenation for that function: it does not consistently add or omit the final value if the value is the min plus the multiple of a step.\n",
    "  3. A change in the `rates_from_double_truncated_mfd` function to use the new functions.\n",
    "  \n",
    "All of this is illustrated below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def round_M_max(M_max, M_min, bin_size, tol=0.0001):\n",
    "    \"\"\"\n",
    "    Rounds `M_max` up to `M_min` plus an integer multiple of `bin_size`.\n",
    "    \n",
    "    :param M_max:\n",
    "        Initial value for maximum earthquake magnitude.\n",
    "        \n",
    "    :type M_max:\n",
    "        float\n",
    "        \n",
    "    :param M_min:\n",
    "        Minimum earthquake magnitude.\n",
    "    \n",
    "    :param bin_size:\n",
    "        Size (or width) of the bins in an evenly-discretized,\n",
    "        truncated Gutenberg-Richter distribution.\n",
    "        \n",
    "    :type bin_size:\n",
    "        float\n",
    "        \n",
    "    :param tol:\n",
    "        Tolerance for deciding whether `M_max` falls on a bin edge.\n",
    "        Should be larger than the floating-point precision but much\n",
    "        smaller than the bin_size.\n",
    "        \n",
    "    :type tol:\n",
    "        float\n",
    "        \n",
    "    :returns:\n",
    "        New (rounded-up) M_max.\n",
    "        \n",
    "    :rtype:\n",
    "        float\n",
    "    \"\"\"\n",
    "    \n",
    "    M_diff = M_max - M_min\n",
    "    if M_diff > 0:\n",
    "        n_bins = M_diff // bin_size\n",
    "        bin_remainder = M_diff % bin_size\n",
    "    \n",
    "        if bin_remainder < tol:\n",
    "            n_bins = n_bins\n",
    "        else:\n",
    "            n_bins = n_bins + 1\n",
    "    \n",
    "        new_M_max = M_min + n_bins * bin_size\n",
    "        \n",
    "    else:\n",
    "        new_M_max = M_max\n",
    "    \n",
    "    return new_M_max\n",
    "\n",
    "\n",
    "def _make_range(vmin, vmax, step, tol=0.0001):\n",
    "    \n",
    "    \"\"\"\n",
    "    Makes a list of equally-spaced values that consistently\n",
    "    omits the final value, unlike numpy.arange\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    num_list = [vmin]\n",
    "    \n",
    "    while num_list[-1] < (vmax - step - tol):\n",
    "        num_list.append(num_list[-1] + step)\n",
    "        \n",
    "    return num_list\n",
    "\n",
    "\n",
    "def rates_for_double_truncated_mfd_roundup(area, slip_rate, m_low, m_upp, b_gr,\n",
    "                                           bin_width=0.1, rigidity=32e9):\n",
    "    \"\"\"\n",
    "    :parameter area:\n",
    "        Area of the fault surface\n",
    "        float [km2]\n",
    "    :parameter slip_rate:\n",
    "        Slip-rate\n",
    "        float [mm/tr]\n",
    "    :parameter m_low:\n",
    "        Lower magnitude\n",
    "        float\n",
    "    :parameter m_upp:\n",
    "        Upper magnitude\n",
    "    :parameter b_gr:\n",
    "        b-value of Gutenber-Richter relationship\n",
    "    :parameter bin_width:\n",
    "        Bin width\n",
    "    :parameter rigidity:\n",
    "        Rigidity [Pa]\n",
    "    :return:\n",
    "        A list containing the rates per bin starting from m_low\n",
    "    \"\"\"\n",
    "    #\n",
    "    # Compute moment\n",
    "    slip_m = slip_rate * 1e-3  # [mm/yr] -> [m/yr]\n",
    "    area_m2 = area * 1e6\n",
    "    moment_from_slip = (rigidity * area_m2 * slip_m)\n",
    "    \n",
    "    # Round up the maximum magnitude\n",
    "    m_upp = round_M_max(m_upp, m_low, bin_width)\n",
    "    \n",
    "    # Compute total rate\n",
    "    rate_above = _get_rate_above_m_low(moment_from_slip, m_low, m_upp, b_gr)\n",
    "    #\n",
    "    # Compute rate per bin\n",
    "    rrr = []\n",
    "    mma = []\n",
    "    for mmm in _make_range(m_low, m_upp, bin_width):\n",
    "        rte = (_get_cumul_rate_truncated(mmm, m_low, m_upp, rate_above, b_gr) -\n",
    "               _get_cumul_rate_truncated(mmm+bin_width, m_low,\n",
    "                                         m_upp, rate_above, b_gr))\n",
    "        mma.append(mmm+bin_width/2)\n",
    "        rrr.append(rte)\n",
    "    return rrr\n",
    "\n",
    "\n",
    "\n",
    "\"\"\"\n",
    "The function below is just a comparison helper function\n",
    "\"\"\"\n",
    "\n",
    "def calc_moment_release_from_fault_params_roundup(area, slip_rate, M_low, M_upp,\n",
    "                                                  b_gr, bin_width=0.1, rigidity=32e9):\n",
    "    \n",
    "    bin_rates = rates_for_double_truncated_mfd_roundup(area, slip_rate, M_low, M_upp,\n",
    "                                                       b_gr, bin_width, rigidity)\n",
    "    \n",
    "    bin_mids = get_bin_midpts(M_low, M_upp, bin_width)\n",
    "    first_bin_mid = bin_mids[0]\n",
    "    \n",
    "    return calc_moment_release_from_bin_rates(first_bin_mid, bin_width, bin_rates)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtsAAAI2CAYAAACBqkT0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xl8VOXZ//HvlSGAYIyIS0EtKCqyJQFCBG0EVBYXELSoqD8BW6lb3aooVZT20VqtW7Vaa58qiJZitahUbF0eKYtaAcUNcasoVJRNAgEiWa7fH+dkmIRJMsEZJgmf9+s1L2fOdl/nzIxcuec6923uLgAAAADJl5HuAAAAAICmimQbAAAASBGSbQAAACBFSLYBAACAFCHZBgAAAFKEZBsAAABIEZJtAAAAIEVItgHsVsxsipndnO44djdmNtLMVphZsZn13IXtnmNmL+yq9mLaPcbMPg7Pd8Subh9Aw0GyDexGzGy5mW0zs32rLX/LzNzMOqYnsvhIjHdkZgPMbGW649gJd0i61N33dPe3UtGAmXUMP8fNKpe5++PuPjgV7dXhl5J+F57v09VXNrbvIoCdR7IN7H4+kzS68oWZ9ZDUKn3hJE9sktUYWaCp/n+5g6T30x3ELpTI+TbZ7yKA7Zrq/9QB1GyapPNiXo+R9GjsBmaWbWaPmtkaM/vczG6oTALNbKyZLTCzu81sg5n9x8yODpevMLPVZjYm5lgtzOwOM/vCzL42swfNbI9w3QAzW2lmPwv3W2Vm48J14yWdI2lC+FP8rHgnE/YCXmJmH0v6OFx2pJm9aGbrzexDMzujpothZqeY2ZLwXF41s5yYddeZ2admtsnMlprZyJh1h5nZv8ysyMzWmtmMmHX1aX+Omd1iZgskbZF0qJmNM7MPwnb/Y2Y/CbdtLel5Se3Da1JsZu3NLCMm1nVm9oSZ7VNDe23M7O/he/tN+PygmPVjwzY3mdlnZnZODccpMLPXwuu2ysx+Z2bN42zXwsyKJUUkvW1mn4bL3cwOi9ku+itGbZ+LcP0eZnZn+NksMrP54WdqbrjJhvDa9AvPZ37Mvkeb2cJwv4VmdnS19+J/ws/3JjN7war1PFc7twvM7JPwfX7WzNqHyz+VdKikWWEcLWo4RJ3fxVranmJmD5jZ82EbC8zse2Z2T/i+LrOYcp06Psu/N7OnYl7fZmYvm5klEguAOrg7Dx48dpOHpOWSTpD0oaQuChKglQp64VxSx3C7RyU9IylLUkdJH0n6UbhurKQySePC/W+W9IWk+yW1kDRY0iZJe4bb3y3pWUn7hMebJenWcN2A8Fi/lJQp6SQFCWebcP0USTfXcU4u6cXw+HtIai1pRRhfM0k9Ja2V1LX6McN1qyUdFZ7LmPAatQjXj5LUXkHHxJmSNktqF66bLun6cF1LST8Il9fafpz454TXr1u4faakkyV1kmSS+ofXpFfMNVtZ7RiXS3pd0kHhe/AHSdNraK+tpNMV9KBmSfqrpKdjYt8oqXP4up2kbjUcp7ekvmHMHSV9IOmKOt6nw2p5Hfu+1PW5uD+8bgeG79vR4Xl3DI/bLOa4YyXND5/vI+kbSf8vjHt0+LptzHvxqaQjFHyW5kj6dQ3nc1z4vvYK275P0tzq37Xv+l2sZf8pYfu9FXz+/k9BT/l52v69fCVm+9o+y60UfMfHSioMj3tQuv9/xYNHU3nQsw3snip71AYpSJL+W7nCzCKSzpI00d03uftySXcqSFAqfebuj7h7uaQZkg6W9Et3/9bdX5C0TdJhYc/YeElXuvt6d98k6Vfh8SuVhvuWuvtsScWSOtfzfG4Nj79V0imSlofxlXlQH/yUgmSjuvGS/uDu/3b3cnefKulbBUmk3P2v7v6lu1e4+wwFPecFMXF3kNTe3UvcvbL3tD7tV5ri7u+H25e6+3Pu/qkH/iXpBQVJUE0ulHS9u690928lTZb0Q4tTVuPu69z9KXffEr4ftyhI6CtVSOpuZnu4+yp3j1sK4e6L3f31MOblChL8/vG23UlxPxcW/MJyvqTL3f2/4fv2anjedTlZ0sfuPi2Me7qkZZKGxWzziLt/FH6WnpCUV8OxzpH0sLu/GbY9UVI/q3+tdY3fxQTMDN+HEkkzJZW4+6Mx38toz3Ztn2V336Lg+32XpMck/dTdG+N9AUCDRLIN7J6mSTpbQU9W9Z+t91XQm/h5zLLPFfQiVvo65vlWSXL36sv2lLSfgl6zxWG5wQZJ/wiXV1rn7mUxr7eE+9bHipjnHSQdVdle2OY5kr4XZ78Okn5WbduDFfQAyszOs+0lJhskdVdwfSRpgoKe5zfM7H0zO38n2o8Xv8zsRDN7PSxP2KCgZ7fGcoawzZkx7X0gqVzSAdU3NLNWZvaHsARjo4LSi73NLOLumxX0el4oaZWZPWdmR8Zr0MyOCEtQvgqP86s6Yqyvmj4X+yroyf10J47ZXlU/19KOn+2v4rRZ57HcvVjSumrHSkRt38W6VP/OxfsOSqrzsyx3/7ek/yj4TD9RzzgA1IJkG9gNufvnCn5yPknS36qtXqvtvbaVvq/69bjFHmurglKEvcNHtrsnmkz7Tmy3QtK/Ytrb24MRIS6Ks98KSbdU27aVu083sw6S/ijpUgVlBntLek9BMiJ3/8rdL3D39pJ+IumBsAa5Pu3vEH9Y3/uUgtE7DgjbnV3Zbg3XZIWkE6u12dLd471nP1Pwy8FR7r6XpGMrmw7P65/uPkhBCcmy8BrE8/tw/eHhcX4eE2MitqjqzYC1/TESa62kEgVlNtXV9Xn5UlU/19LOf7arHMuCevq29T1WHd/FpKjrsxxuc4mCcpgvFfwhCSBJSLaB3dePJB0X9mZGhT9BPyHpFjPLCv+hvkrBz8v14u4VCv6Rv9vM9pckMzvQzIYkeIivFdxoVh9/l3SEmf0/M8sMH33MrEucbf8o6UIzO8oCrc3sZDPLUlC/7JLWhHGPU9AbqPD1KNt+Y+E34bYV9Ww/nuYKkp41ksrM7EQFdfCVvpbU1syyY5Y9qOD96hDGtp+ZnVrD8bMU/AG0wYKbKG+KOacDzOzUMHH8VkHpRkUtx9koqTjs/a7tj4l4lkg628wiZjZUCZaghJ+phyXdZcHNoRELboSsvGYVqvkzM1vBe3O2mTUzszMldVXwntXXdEnjzCwvbPtXkv4dltTUV9zvYhLV9Vk+QkGN97kKykkmmFlN5TMA6olkG9hNhTXBi2pY/VMFN1D9R9J8SX9WkODsjGslfSLp9bDc4CUlXpP9J0ldw5++dxirOJ6wDnmwgrrwLxWUBdymIIGtvu0iSRdI+p2ChPkTBT/ny92XKqhVf01BgttD0oKY3ftI+rcFI208q6CG+D/1ab+W+C9T8AfPNwpKDJ6NWb9MQaL3n/C6tJf023CbF8xsk4KbJY+qoYl7FNz8tzbc7h8x6zIU/GH1paT1ChLgmpLoq8PYNin4o2VGDdvV5HIFtdKVZTYJvb8xbb8raWEY522SMsLa41skLQivTd/Yndx9nYKa+p8pKPmYIOkUd19bz9jl7i9JmqTgV4hVCnraz6p1p5qPVdt38Tur7bMc1vU/Juk2d3/b3T9W8CvFNKt5FBUA9WDuif5KCwAAAKA+6NkGAAAAUoRkGwCABigc5aY4ziPuREMAGibKSAAAAIAUoWcbAAAASBGSbQBNgpndamZXpKntQjP7MB1tJ4OZzTGzH9ew7vth6UIkRW0/b2ZjUnHsVDKzN8ysW7rjANDwkWwDaPTMbD8FU17/IXw9wMzczGZW2y43XD4nme27+zx3r+8U83UKZ3t8wMzWmlmRmc2tZds5ZlYSU9eblOTf3b8IJ+UpT8bx4hz/RHefmopjV2dmLczsYTPbGM58eVUt23Y3s3+G1z5eveUdkn6ZumgBNBUk2wCagrGSZrv71phlayT1M7O2McvGSPpoVwb2HT0kaR9JXcL/XlnH9peGifGeqUj+m4DJkg5XMPPjQAWTtwytYdtSBWOd/6iG9c9KGmhmic58CWA3RbINoCk4UdK/qi3bpmCilLMkKSyDOFPS44kcMOwBv9jMPjazTWb2P2bWycxeDXtGnzCz5uG2A8xsZcy+y83sajN7J+yRnmFmLetzQuGsjMMljXf3Ne5e7u6L63OMeuoUlkZsNLNnwtklZWYdw2vRLHw9J7wWC8Lr8oKZ7VvHubQ0s8fMbF042cxCMzsg5ng/Dp+/XW3UDTezAeG6vuG13xBuN2AnznGMpP9x92/c/QMFk/GMjbehu3/o7n+S9H4N60skLZaU6GyoAHZTJNsAmoIekuKVTTyqoLxECpKi9xTMjpioIZJ6S+qrYLbBhxRMaX2wgumuR9ey7xmShko6RFKOwqQurIHeUMvj7HD/AkmfS/pFWMrwrpmdXke8t4bbLtiJZPQ8SedLaiepTNK9tWx7tqRxkvZXML381XUce4ykbAXXra2kCxVMGV+Fu+dW9swrmMnyQ0lvmtmBkp5TMKX4PmF7T4XlQwpLbWq6nu+E27QJz+3tmCbflvRd6q4/kJT7HfYHsBsg2QbQFOytYNrwKtz9VUn7mFlnBcnko/U87u3uvtHd31eQqL8QTsleJOl5ST1r2fded//S3ddLmiUpL4zpC3ffu5bHn8P9D1KQ0BdJai/pUklTzaxLDe1dK+lQSQcq+KNglpl1qse5TnP399x9s4JpyM+o5abIR9z9o7Bs54nKc6tFqYIk+7DKHnp331jTxmb2AwWJ9fBwu3MVlAnNdvcKd39R0iJJJ0mSu19cy/XMCQ+7Z/jfopimiiRl1RF7bTYp+OwBQI1ItgE0Bd+o5qRpmoJEdaCkmTVsU5OvY55vjfN6T9Xsq5jnW+rYNp6tCpLUm919m7v/S9IrkgbH29jd/+3um9z92/CGwwUKk9EErYh5/rmkTEk1lYfU99ymSfqnpL+Y2ZdmdruZZcbb0MwOVpDAj3H3yvr6DpJGxfZYS/qBgp7qRBWH/90rZtleivNHWj1kSdrwHfYHsBsg2QbQFLwj6Yga1k2TdLGCntEtuy6k+Gz7UHo1PSpnB3wnzu71mYXMJVk9tj845vn3FST6a+uxf82BuJe6+y/cvaukoyWdou3lPVFmtoeCOvt73P35mFUrFPS8x/ZYt3b3X4f7PVjL9Xw/jOEbSatUtewjVzXUZCeoi6qWpQDADki2ATQFsyX1j7fC3T8L112/SyOqQcxQejU9Km/gnCvpC0kTzayZmR2joHf+n9WPaWZ7m9mQ8EbEZmHCfqykf4TrK29y7FhLaOeaWVcza6VgSLsnkzXcn5kNNLMeYVnKRgWJfEWcTR+WtMzdb6+2/DFJw8JzjITnOcDMDpIkd7+wlusZW5P9qKQbzKxNeAPqBZKm1BCzhTe1Vt4E29LMWsSsb6mgnv/FnbgkAHYjJNsAmoJHJZ0U9ozuwN3nu3t9boxMO3cvlXSqglKQIgUjZ5zn7sskycx+bmaVvb+ZCmqc1yjojf6ppBExZRgHKygN+W8tTU5TkHh+JamlpMuSeDrfk/SkgkT7AwUjx0yLs91ZkkZW65kudPcVCq7FzxWc4wpJ16j+/4bdJOlTBdfiX5J+4+6Vf5BU/uLw/XDbDgpKeSp7vreq6k24wyTNaWyfKwC7nrnX51dJAGiYzOxXkla7+z3pjqWhMbMbJK1x9z+kO5amwsz+LelH7v5eumMB0LCRbAMAAAApQhkJgN2SmRXWdFNdumNrjMzsnNpuUASA3RU92wAAAECK0LMNAAAApEizdAeQTPvuu6937Ngx3WEAAACgiVu8ePFad9+vru2aVLLdsWNHLVq0KN1hAAAAoIkzs88T2a5JlJGY2TAze6ioqCjdoQAAAABRTSLZdvdZ7j4+Ozs73aEAAAAAUU0i2QYAAAAaoiZVsw0AAJKvtLRUK1euVElJSbpDAXa5li1b6qCDDlJmZuZO7U+yDQAAarVy5UplZWWpY8eOMrN0hwPsMu6udevWaeXKlTrkkEN26hiUkQAAgFqVlJSobdu2JNrY7ZiZ2rZt+51+1SHZBgAAdSLRxu7qu372SbYBAECDt3LlSp166qk6/PDD1alTJ11++eXatm1bdP38+fNVUFCgI488Up07d9YDDzxQ47E6duyowsLCKsvy8vLUvXt3SdKcOXOUnZ2tvLw85eTk6IQTTtDq1av1yCOPKC8vT3l5eWrevLl69OihvLw8XXfddak56SR7+umntXTp0u90jDlz5uiUU05JUkTJN2fOHL366qvpDqMKkm0AANCgubtOO+00jRgxQh9//LE++ugjFRcX6/rrr5ckffXVVzr77LP14IMPatmyZVqwYIH+9Kc/aebMmTUec9OmTVqxYoUk6YMPPthhfWFhoZYsWaJ33nlHffr00f33369x48ZpyZIlWrJkidq3b69XXnlFS5Ys0a9//evUnHiSJZpsl5WV7YJodl5t8ZFsAwAA1NP//d//qWXLlho3bpwkKRKJ6O6779bDDz+sLVu26P7779fYsWPVq1cvSdK+++6r22+/Xb/5zW9qPOYZZ5yhGTNmSJKmT5+u0aNHx93O3bVp0ya1adMm4XinTJmiESNGaNCgQerYsaN+97vf6a677lLPnj3Vt29frV+/XpK0ZMkS9e3bVzk5ORo5cqS++eYbSdKAAQN05ZVXKj8/X126dNHChQt12mmn6fDDD9cNN9wQbeexxx5TQUGB8vLy9JOf/ETl5eWSpD333FPXX3+9cnNz1bdvX3399dd69dVX9eyzz+qaa65RXl6ePv300yoxjx07VhdeeKGOOuooTZgwQZs3b9b555+vgoIC9ezZU88888wO51nTNsuXL1dhYaF69eqlXr16RZPfVatW6dhjj43+ijBv3jxJ0gsvvKB+/fqpV69eGjVqlIqLi3doa8CAAbriiiuUn5+v3/72t5o1a5aOOuoo9ezZUyeccIK+/vprLV++XA8++KDuvvtu5eXlad68eVqzZo1OP/109enTR3369NGCBQsSfh+ThdFIAABAwn4x630t/XJjUo/Ztf1eumlYtxrXv//+++rdu3eVZXvttZe+//3v65NPPtH777+vMWPGVFmfn59fay/u6aefrnHjxunqq6/WrFmz9Pjjj2vatGnR9fPmzVNeXp7WrVun1q1b61e/+lW9zum9997TW2+9pZKSEh122GG67bbb9NZbb+nKK6/Uo48+qiuuuELnnXee7rvvPvXv31833nijfvGLX+iee+6RJDVv3lyLFi3Sb3/7W5166qlavHix9tlnH3Xq1ElXXnmlVq9erRkzZmjBggXKzMzUxRdfrMcff1znnXeeNm/erL59++qWW27RhAkT9Mc//lE33HCDhg8frlNOOUU//OEP48a8cuVKvfrqq4pEIvr5z3+u4447Tg8//LA2bNiggoICnXDCCVW2v+WWW+Jus//+++vFF19Uy5Yt9fHHH2v06NFatGiR/vznP2vIkCG6/vrrVV5eri1btmjt2rW6+eab9dJLL6l169a67bbbdNddd+nGG2/cIb5t27Zp0aJFkqRvvvlGr7/+usxM//u//6vbb79dd955py688ELtueeeuvrqqyVJZ599tq688kr94Ac/0BdffKEhQ4bE/SUjlUi2AQDAbqdt27Zq06aN/vKXv6hLly5q1apVlfWFhYX6+9//Lkm67bbbNGHCBD344IMJH3/gwIHKyspSVlaWsrOzNWzYMElSjx499M4776ioqEgbNmxQ//79JUljxozRqFGjovsPHz48un23bt3Url07SdKhhx6qFStWaP78+Vq8eLH69OkjSdq6dav2339/SUGiXllX3bt3b7344osJxTxq1ChFIhFJQW/zs88+qzvuuENSMCLNF198UWX7mrZp3769Lr30Ui1ZskSRSEQfffSRJKlPnz46//zzVVpaqhEjRigvL0//+te/tHTpUh1zzDGSgoS6X79+ceM788wzo89XrlypM888U6tWrdK2bdtqHJbvpZdeqvJH18aNG1VcXKw999wzoWuSDCTbAAAgYbX1QKdK165d9eSTT1ZZtnHjRn3xxRc67LDD1LVrVy1evFinnnpqdP3ixYuVn5+v8vLyaK/48OHD9ctf/jK6zZlnnqlLLrlEU6ZMqbX94cOH6/TTT69XzC1atIg+z8jIiL7OyMhIqCY6dvvqxyorK5O7a8yYMbr11lt32DczMzM6gkYkEkm4Brt169bR5+6up556Sp07d66yzddff13nNpMnT9YBBxygt99+WxUVFWrZsqUk6dhjj9XcuXP13HPPaezYsbrqqqvUpk0bDRo0SNOnT69XfD/96U911VVXafjw4ZozZ44mT54cd5+Kigq9/vrr0RjSgZptAADQoB1//PHasmWLHn30UUlSeXm5fvazn2ns2LFq1apVNGFesmSJJGndunW6/vrrNWnSJEUikehNjbGJtiSNHDlSEyZM0JAhQ2ptf/78+erUqVNSzyk7O1tt2rSJ1i1PmzYt2sudiOOPP15PPvmkVq9eLUlav369Pv/881r3ycrK0qZNmxI6/pAhQ3TffffJ3SVJb731VsLbFBUVqV27dsrIyNC0adOiteSff/65DjjgAF1wwQX68Y9/rDfffFN9+/bVggUL9Mknn0gK6sAre8JrU1RUpAMPPFCSNHXq1BrPcfDgwbrvvvuirys/I7sSyTYAAGjQzEwzZ87UX//6Vx1++OE64ogj1LJly2gddbt27fTYY49p/Pjx6ty5s9q3b6/LLruszuQ1KytL1157rZo3b77Dusqa7dzcXE2bNk133nln0s9r6tSpuuaaa5STk6MlS5bErVOuSdeuXXXzzTdr8ODBysnJ0aBBg7Rq1apa9znrrLP0m9/8Rj179tzhBsnqJk2apNLSUuXk5Khbt26aNGlSwttcfPHFmjp1qnJzc7Vs2bJoj/ScOXOUm5urnj17asaMGbr88su13377acqUKRo9erRycnLUr18/LVu2rM7znzx5skaNGqXevXtr3333jS4fNmyYZs6cGb1B8t5779WiRYuUk5Ojrl271qsUKFms8q+RpiA/P98rC+cBAEByfPDBB+rSpUu6w0jYAw88oN///veaO3duvUYRAWoS7ztgZovdPb+ufenZBgAATcrFF1+sd999l0QbDQLJNgAAAJAiJNsAAABAipBsJ8F5D7+hG595L91hAAAAoIFhnO0kWLVhq1o3j6Q7DAAAADQw9GwngZnUhAZ1AQAAQJKQbCdBhplcZNsAAKTKV199pbPOOkudOnVS7969ddJJJyU0+UmyLFmyRLNnz673fgMGDFBjGZZ4w4YNeuCBB77zcRr6OVeOz76rkGwnSQW5NgAAKeHuGjlypAYMGKBPP/1Uixcv1q233lpl6vDaVJ+u3N1VUVFRrxh2NtluTBJNtnfm+u1KdU1PT7LdCJkZZSQAAKTIK6+8oszMTF144YXRZbm5uSosLJS765prrlH37t3Vo0cPzZgxQ1IwW2FhYaGGDx+url27avny5ercubPOO+88de/eXStWrNALL7ygfv36qVevXho1apSKi4slSQsXLtTRRx+t3NxcFRQUqKioSDfeeKNmzJihvLw8zZgxQ5s3b9b555+vgoIC9ezZU88884wkaevWrTrrrLPUpUsXjRw5Ulu3bo17Th07dtTEiROVl5en/Px8vfnmmxoyZIg6deoUneWwtnPr37+/Tj31VB166KG67rrr9Pjjj6ugoEA9evSIzg65Zs0anX766erTp4/69OmjBQsWSApmXzz//PM1YMAAHXroobr33nslSdddd50+/fRT5eXl6ZprrqkSb32uX6yatvnlL3+pPn36qHv37ho/fnx0yvd7771XXbt2VU5Ojs466yxJqvFax6r+fkvSiBEj1Lt3b3Xr1k0PPfRQ9By3bt2qvLw8nXPOOZKkxx57TAUFBcrLy9NPfvKT6PTySePuTebRu3dvT4cT75nrP5ryRlraBgAg1ZYuXbr9xexr3R8+KbmP2dfW2v5vf/tbv+KKK+Kue/LJJ/2EE07wsrIy/+qrr/zggw/2L7/80l955RVv1aqV/+c//3F3988++8zNzF977TV3d1+zZo0XFhZ6cXGxu7v/+te/9l/84hf+7bff+iGHHOJvvBH8u15UVOSlpaX+yCOP+CWXXBJtd+LEiT5t2jR3d//mm2/88MMP9+LiYr/zzjt93Lhx7u7+9ttveyQS8YULF+4Qd4cOHfyBBx5wd/crrrjCe/To4Rs3bvTVq1f7/vvvX+e5ZWdn+5dffuklJSXevn17v/HGG93d/Z577vHLL7/c3d1Hjx7t8+bNc3f3zz//3I888kh3d7/pppu8X79+XlJS4mvWrPF99tnHt23b5p999pl369Yt7nVO9Pq5u/fv398XLlxY6zbr1q2LHvvcc8/1Z5991t3d27Vr5yUlJdHrWtu1jlX9/Y5tY8uWLd6tWzdfu3atu7u3bt06us3SpUv9lFNO8W3btrm7+0UXXeRTp07d4fyrfAdCkhZ5Avkpo5EkQUYGN0gCAJAO8+fP1+jRoxWJRHTAAQeof//+Wrhwofbaay8VFBTokEMOiW7boUMH9e3bV5L0+uuva+nSpTrmmGMkSdu2bVO/fv304Ycfql27durTp48kaa+99orb7gsvvKBnn31Wd9xxhySppKREX3zxhebOnavLLrtMkpSTk6OcnJwaYx8+fLgkqUePHiouLlZWVpaysrLUokULbdiwodZz69Onj9q1aydJ6tSpkwYPHhw91iuvvCJJeumll7R06dJoexs3boz2LJ988slq0aKFWrRoof333z+hkpxErl+s2rZ55ZVXdPvtt2vLli1av369unXrpmHDhiknJ0fnnHOORowYoREjRtR6ratPn179/b733ns1c+ZMSdKKFSv08ccfq23btlX2efnll7V48eLo+71161btv//+dV6L+iDZTgKTqYJsGwCwOzjx17u8yW7duunJJ5+s936tW7eu8bW7a9CgQZo+fXqVbd59992Eju3ueuqpp9S5c+d6x1WpRYsWkqSMjIzo88rXddUdV98+9liV+1ZUVOj1119Xy5Yta90/EonU2Z6U2PWLVdM2JSUluvjii7Vo0SIdfPDBmjx5skpKSiS9/iNNAAAgAElEQVRJzz33nObOnatZs2bplltu0bvvvpvwtY6Nb86cOXrppZf02muvqVWrVhowYEC0jeoxjhkzRrfeemud57+zqNlOAjMxFgkAACly3HHH6dtvv43W3UrSO++8o3nz5qmwsFAzZsxQeXm51qxZo7lz56qgoKDOY/bt21cLFizQJ598IimoC/7oo4/UuXNnrVq1SgsXLpQkbdq0SWVlZcrKytKmTZui+w8ZMkT33XdftNb4rbfekiQde+yx+vOf/yxJeu+99/TOO+/s9Hnv7LlVGjx4sO67777o6yVLltS6ffVzrE1N1y+RbSqT3n333VfFxcXRP6QqKiq0YsUKDRw4ULfddpuKiopUXFxc47WuTVFRkdq0aaNWrVpp2bJlev3116PrMjMzVVpaKkk6/vjj9eSTT2r16tWSpPXr1+vzzz9P6BokimQ7CbhBEgCA1DEzzZw5Uy+99JI6deqkbt26aeLEifre976nkSNHKicnR7m5uTruuON0++2363vf+16dx9xvv/00ZcoUjR49Wjk5OerXr5+WLVum5s2ba8aMGfrpT3+q3NxcDRo0SCUlJRo4cKCWLl0avUFy0qRJKi0tVU5Ojrp166ZJkyZJki666CIVFxerS5cuuvHGG9W7d++dPu+dPbdK9957rxYtWqScnBx17do1euNlTdq2batjjjlG3bt33+EGyepqun6JbLP33nvrggsuUPfu3TVkyJBoCUd5ebnOPfdc9ejRQz179tRll12mvffeu8ZrXZuhQ4eqrKxMXbp00XXXXRctf5Gk8ePHR8tVunbtqptvvlmDBw9WTk6OBg0apFWrVtV5/Powb0JZYn5+vqdjXMcR9y9QVstmmvajo3Z52wAApNoHH3ywQ30ssDuJ9x0ws8Xunl/XvvRsJ4FZuiMAAABAQ0SynQQmRiMBAADAjki2k4Dp2gEAABAPyXYSmEkNeNZSAAAApAnJdhKY6NkGAADAjki2k8CMmm0AAADsiGQ7CUi2AQBILTPTueeeG31dVlam/fbbT6ecckpa4lmyZIlmz56dlra/qwEDBqiuoZLvuecebdmyJfr6pJNO0oYNG1IdWlRjvr7VNYlk28yGmdlDRUVF6WmfMhIAAFKqdevWeu+997R161ZJ0osvvqgDDzwwbfE0pWQwnurJ9uzZs7X33nsntY3apohvSte3SSTb7j7L3cdnZ2enpf2MDHq2AQBItZNOOknPPfecJGn69OkaPXp0dN369es1YsQI5eTkqG/fvtFp0idPnqwxY8aosLBQHTp00N/+9jdNmDBBPXr00NChQ6PTdi9evFj9+/dX7969NWTIkOgsggMGDNC1116rgoICHXHEEZo3b562bdumG2+8UTNmzIjOKBlr+fLlKiwsVK9evdSrVy+9+uqr0XW33XabevToodzcXF133XWSpE8++UQnnHCCcnNz1atXL3366aeaM2dOlV77Sy+9VFOmTJEkdezYURMnTlReXp7y8/P15ptvasiQIerUqVN0lsja9o910UUXKT8/X926ddNNN90kKZh58ssvv9TAgQM1cODAaJtr166VJN11113q3r27unfvrnvuuSd6zl26dNEFF1ygbt26afDgwdE/jGKNHTtWF154oY466ihNmDBBb7zxhvr166eePXvq6KOP1ocffhj3+m7evFnnn3++CgoK1LNnTz3zzDM1f1AaGndvMo/evXt7Opzzx9d95P3z09I2AACptnTp0nSH4K1bt/a3337bTz/9dN+6davn5ub6K6+84ieffLK7u1966aU+efJkd3d/+eWXPTc3193db7rpJj/mmGN827ZtvmTJEt9jjz189uzZ7u4+YsQInzlzpm/bts379evnq1evdnf3v/zlLz5u3Dh3d+/fv79fddVV7u7+3HPP+fHHH+/u7o888ohfcsklcWPdvHmzb9261d3dP/roI6/MT2bPnu39+vXzzZs3u7v7unXr3N29oKDA//a3v7m7+9atW33z5s1Vzs3d/ZJLLvFHHnnE3d07dOjgDzzwgLu7X3HFFd6jRw/fuHGjr1692vfff39391r379+/vy9cuLBKDGVlZd6/f39/++23o22sWbMmun/l60WLFnn37t29uLjYN23a5F27dvU333zTP/vsM49EIv7WW2+5u/uoUaN82rRpO1ybMWPG+Mknn+xlZWXu7l5UVOSlpaXu7v7iiy/6aaedFvf6Tpw4MXq8b775xg8//HAvLi6Oe/1TId53QNIiTyA/bZbuZL8pMBNFJACA3ceUKdLy5ck7XseO0tixdW6Wk5Oj5cuXa/r06TrppJOqrJs/f76eeuopSdJxxx2ndevWaePGjZKkE088UZmZmerRo4fKy8s1dOhQSVKPHj20fPlyffjhh3rvvfc0aNAgSVJ5ebnatWsXPfZpp50mSerdu7eWJ3DepaWluvTSS7VkyRJFIhF99NFHkqSXXnpJ48aNU6tWrSRJ++yzjzZt2qT//ve/GjlypCSpZcuWdR5fkoYPHx49h+LiYmVlZSkrK0stWrSoV231E088oYceekhlZWVatWqVli5dqpycnBq3nz9/vkaOHKnWrVtLCq7NvHnzNHz4cB1yyCHKy8uTVPu1GjVqlCKRiCSpqKhIY8aM0ccffywzi/7SUN0LL7ygZ599VnfccYckqaSkRF988cUOU6g3RCTbSUIZCQBgt5FAYpwqw4cP19VXX605c+Zo3bp1Ce3TokULSVJGRoYyMzNlZtHXZWVlcnd169ZNr732Wq37RyKRWuuMK91999064IAD9Pbbb6uioiLhBDpWs2bNVBEziUdJSUmN51T5PPac6tpfkj777DPdcccdWrhwodq0aaOxY8fG3S5RsXFEIpG4ZSSSoom6JE2aNEkDBw7UzJkztXz5cg0YMCDuPu6up556Sp07d97p+NKlSdRsp1swgyQAAEi1888/XzfddJN69OhRZXlhYaEef/xxSUG98r777qu99toroWN27txZa9asiSbbpaWlev/992vdJysrS5s2bYq7rqioSO3atVNGRoamTZum8vJySdKgQYP0yCOPRG88XL9+vbKysnTQQQfp6aefliR9++232rJlizp06KClS5fq22+/1YYNG/Tyyy8ndC6VEtl/48aNat26tbKzs/X111/r+eefr/P8CgsL9fTTT2vLli3avHmzZs6cqcLCwnrFFquoqCh6o2tsTXn19ocMGaL77rtPHvZuvvXWWzvd5q5Gsp0EwdB/pNsAAKTaQQcdpMsuu2yH5ZMnT9bixYuVk5Oj6667TlOnTk34mM2bN9eTTz6pa6+9Vrm5ucrLy6tyU2M8AwcO1NKlS+PeIHnxxRdr6tSpys3N1bJly6I9uUOHDtXw4cOVn5+vvLy8aEnEtGnTdO+99yonJ0dHH320vvrqKx188ME644wz1L17d51xxhnq2bNnwucjKaH9c3Nz1bNnTx155JE6++yzdcwxx0TXjR8/XkOHDo3eIFmpV69eGjt2rAoKCnTUUUfpxz/+cb1jizVhwgRNnDhRPXv2rPKrQfXrO2nSJJWWlionJ0fdunXTpEmTdrrNXc2aUpKYn5/vdY0bmQrjHnlDa4u3adZPf7DL2wYAINU++OCDRlEbC6RKvO+AmS129/y69qVnOwmCMpKm80cLAAAAkoNkOwnMpJh7EAAAAABJJNtJwg2SAAAA2BHJdhJkcIMkAKCJ49857K6+62efZDsJgtFI0h0FAACp0bJlS61bt46EG7sdd9e6det2aqz0SkxqkwQmbpAEADRdBx10kFauXKk1a9akOxRgl2vZsqUOOuignd6fZDsJ6NkGADRlmZmZOuSQQ9IdBtAoUUaSBMwgCQAAgHhItpPBpAq6tgEAAFANyXYSmCS6tgEAAFAdyXYSUEYCAACAeEi2k8AoIwEAAEAcJNtJYGI0EgAAAOyIZDsJzBhnGwAAADsi2U4CM6miIt1RAAAAoKEh2U4CC8YjAQAAAKog2U6CYAZJykgAAABQFcl2EmQYw2wDAABgRyTbSWAyhv4DAADADki2kyAoI0l3FAAAAGhoSLaTwJhBEgAAAHGQbCcBN0gCAAAgHpLtJGAGSQAAAMRDsp0ExmgkAAAAiINkOwkyzCgjAQAAwA5ItpPAJFWQawMAAKAaku0kMHq2AQAAEAfJdhJQsw0AAIB4SLaTwGSMRgIAAIAdkGwnAeNsAwAAIB6S7STIoIwEAAAAcZBsJ4GZqYKebQAAAFRDsp0EzCAJAACAeEi2k4EyEgAAAMRBsp0EGYz9BwAAgDhItpMgmEGSbBsAAABVkWwnAR3bAAAAiIdkOwkymK4dAAAAcZBsJ0FQRpLuKAAAANDQkGwng1m6IwAAAEADRLKdBBlhrk0pCQAAAGKRbCeBKci2KSUBAABALJLtJDB6tgEAABAHyXYSVFZsk2oDAAAgFsl2EmSERdt0bAMAACAWyXYSMYskAAAAYpFsJwEj/wEAACAeku0kyDDKSAAAALAjku0kqOzYpowEAAAAsUi2kyA69F96wwAAAEADQ7KdBJWT2jDONgAAAGKRbCcBPdsAAACIh2Q7CazyBsmKNAcCAACABoVkOwm2zyBJ3zYAAAC2I9lOgozKMhJybQAAAMRolu4AKpnZoZKul5Tt7j8MlxVKOkdBnF3d/eg0hlijyjIShv4DAABArJT2bJvZw2a22szeq7Z8qJl9aGafmNl1kuTu/3H3H8Vu5+7z3P1CSX+XNDWVsX4X3CAJAACAeFJdRjJF0tDYBWYWkXS/pBMldZU02sy61nGcsyX9ORUBJoMxgyQAAADiSGmy7e5zJa2vtrhA0idhT/Y2SX+RdGpNxzCz70sqcvdNNawfb2aLzGzRmjVrkhV6vURvkCTbBgAAQIx03CB5oKQVMa9XSjrQzNqa2YOSeprZxJj1P5L0SE0Hc/eH3D3f3fP322+/1ERcB8pIAAAAEE+DuUHS3ddJujDO8pvSEE69bJ9BMs2BAAAAoEFJR8/2fyUdHPP6oHBZoxUd+o++bQAAAMRIR7K9UNLhZnaImTWXdJakZ9MQR9JUlpFUkGsDAAAgRqqH/psu6TVJnc1spZn9yN3LJF0q6Z+SPpD0hLu/n8o4Um17GQnZNgAAALZLac22u4+uYflsSbNT2fauZMwgCQAAgDiYrj0JGGcbAAAA8ZBsJ0F0nG1ukAQAAEAMku0kyAivIj3bAAAAiNUkkm0zG2ZmDxUVFaWn/bBvu4JsGwAAADGaRLLt7rPcfXx2dnZa2mcGSQAAAMTTJJLthoKObQAAAMQi2U6CDNt+iyQAAABQiWQ7CZhBEgAAAPGQbCfB9hkk0xwIAAAAGhSS7STIiN4gSbYNAACA7Ui2kyBaRlKR3jgAAADQsJBsJ0VYRkLPNgAAAGKQbCdBtIyEXBsAAAAxSLaTwIwbJAEAALAjku0kyIgO/Ue2DQAAgO2aRLJtZsPM7KGioqK0tJ8RZtvlJNsAAACI0SSSbXef5e7js7Oz09J+JCwjqWBWGwAAAMRoEsl2ukXCnm1ybQAAAMQi2U6CynG2y8m2AQAAEINkOwmiZSTUbAMAACAGyXYSVJaR0LMNAACAWCTbSZCRQc82AAAAdkSynQQZlJEAAAAgDpLtJKis2S6vSHMgAAAAaFBItpMgI7yK1GwDAAAgFsl2EkSo2QYAAEAcJNtJwNB/AAAAiKdJJNtmNszMHioqKkpX+5IoIwEAAEBVTSLZdvdZ7j4+Ozs7Le1TRgIAAIB4mkSynW6MRgIAAIB4SLaToHI0kgrKSAAAABCDZDsJKCMBAABAPCTbSVA5g2Q5yTYAAABikGwnQXS6dspIAAAAEINkOwkqy0gY+g8AAACxSLaTYPukNmkOBAAAAA0KyXYSWOVoJNRsAwAAIEaz2laaWT9J50oqlNRO0lZJ70l6TtJj7p6eKRsbmAgzSAIAACCOGnu2zex5ST+W9E9JQxUk210l3SCppaRnzGz4rgiyoYvWbNOzDQAAgBi19Wz/P3dfW21ZsaQ3w8edZrZvyiJrRBiNBAAAAPHU2LPt7mvNLGJmr9S2TWrCaly2T2qT5kAAAADQoNR6g6S7l0uqMLPsXRRPoxTm2tRsAwAAoIpab5AMFUt618xelLS5cqG7X5ayqOrJzIZJGnbYYYelq32ZMRoJAAAAqkok2f5b+Giw3H2WpFn5+fkXpCuGiBk92wAAAKgikWR7vaTn3L0i1cE0ZhkZxmgkAAAAqCKRSW3OlPSxmd1uZkemOqDGKmImcm0AAADEqjPZdvdzJfWU9KmkKWb2mpmNN7OslEfXiGQYN0gCAACgqoSma3f3jZKelPQXBZPbjJT0ppn9NIWxNSoZGdRsAwAAoKo6k20zG25mMyXNkZQpqcDdT5SUK+lnqQ2v8YhkGKORAAAAoIpEbpA8XdLd7j43dqG7bzGzH6UmrMYnYiTbAAAAqKrOZNvdx9Sy7uXkhtN4mZnKGa8FAAAAMRKq2UbdIhlSBTXbAAAAiEGynSQRY5xtAAAAVFVjsm1mD5nZSIb4S0xGhtGzDQAAgCpq69n+k4IRR2ab2ctmdq2Z5e6iuBodRiMBAABAdTXeIOnu/5b0b0mTzaytpMGSfmZmPSS9Jekf7v7Ergmz4cswUzm5NgAAAGIkMvSf3H2dpOnhQ2bWW9LQFMbV6GQYN0gCAACgqoSS7ercfbGkxUmOpVGLMIMkAAAAqmE0kiTJYDQSAAAAVNMkkm0zG2ZmDxUVFaUthkiGyUm2AQAAEKPOMhIzi0g6WVLH2O3d/a7UhVU/7j5L0qz8/PwL0hVDhlFGAgAAgKoSqdmeJalE0ruSmJC8BhkZjEYCAACAqhJJtg9y95yUR9LIRRiNBAAAANUkUrP9vJkNTnkkjVyzjAyVVdDxDwAAgO0S6dl+XdJMM8uQVCrJJLm775XSyBqZSIaRbAMAAKCKRJLtuyT1k/SuM9xGjZpFTCVlXB4AAABsl0gZyQpJ75Fo164Zk9oAAACgmkR6tv8jaY6ZPS/p28qFDWnov4YgkpGhMoYjAQAAQIxEku3Pwkfz8IE46NkGAABAdXUm2+7+i10RSGMXiZhKuUESAAAAMZrEdO0NQSY92wAAAKiGZDtJqNkGAABAdSTbSULNNgAAAKqrsWbbzG6sZT939/9JQTyNViRiKiPZBgAAQIzabpDcHGdZK0k/ltRWEsl2jGbMIAkAAIBqaky23f3OyudmliXpcknnS/qLpDtr2m931SwjQ+XUbAMAACBGrUP/mdk+kq6SdI6kqZJ6ufs3uyKwxqYZZSQAAACopraa7d9IOk3SQ5J6uHvxLouqEYpwgyQAAACqqW00kp9Jai/pBklfmtnG8LHJzDbumvAaD2q2AQAAUF1tNdsMC1gPkQxThUsVFa6MDEt3OAAAAGgAakyozWzPunZOZJvdRWYkuJTUbQMAAKBSbb3Xz5jZnWZ2rJm1rlxoZoea2Y/M7J+ShqY+xLqZ2TAze6ioqChtMUTC3mzqtgEAAFCpxmTb3Y+X9LKkn0h638yKzGydpMckfU/SGHd/cteEWTt3n+Xu47Ozs9MWQ7Mw2aZuGwAAAJVqHfrP3WdLmr2LYmnU6NkGAABAddwEmSSVPdulTGwDAACAEMl2kkQygktJzzYAAAAqkWwnSbMINdsAAACoKqFk28x+YGbjwuf7mdkhqQ2r8WlGzTYAAACqqTPZNrObJF0raWK4KFPBiCSIEYmORkKyDQAAgEAiPdsjJQ2XtFmS3P1LSVmpDKoxahbWbJdxgyQAAABCiSTb29zdJbkkxU5wg+0ijLMNAACAahJJtp8wsz9I2tvMLpD0kqT/TW1YjU9mhJptAAAAVFXrpDaS5O53mNkgSRsldZZ0o7u/mPLIGhlqtgEAAFBdncm2md3m7tdKejHOMoSaMc42AAAAqkmkjGRQnGUnJjuQxi4SnUGSmm0AAAAEauzZNrOLJF0s6VAzeydmVZakBakOrLFp3iwsI2E0EgAAAIRqKyP5s6TnJd0q6bqY5ZvcfX1Ko2qEMiPBjwT0bAMAAKBSjcm2uxdJKpI0WpLMbH9JLSXtaWZ7uvsXuybExoFkGwAAANUlMoPkMDP7WNJnkv4labmCHm/EqEy2t1FGAgAAgFAiN0jeLKmvpI/c/RBJx0t6PaVRNULNK3u2y+jZBgAAQCCRZLvU3ddJyjCzDHd/RVJ+iuNqdDKbMRoJAAAAqqpznG1JG8xsT0lzJT1uZqslbU5tWI0PNdsAAACoLpGe7VMlbZF0paR/SPpU0rBUBtUYUbMNAACA6mrt2TaziKS/u/tASRWSpu6SqBqhzEjlONv0bAMAACBQa8+2u5dLqjCz7F0UT6NFGQkAAACqS6Rmu1jSu2b2omJqtd39spRF1Qg1C6drp4wEAAAAlRJJtv8WPlALM1PzSAY92wAAAIiqM9l2d+q0E5QZMcbZBgAAQFQio5E0eOEslw8VFRWlNY7MZvRsAwAAYLsmkWy7+yx3H5+dnd77ODMjGdRsAwAAICrhZNvMWqUykKaAmm0AAADEqjPZNrOjzWyppGXh61wzeyDlkTVCmREj2QYAAEBUIj3bd0saImmdJLn725KOTWVQjVUmPdsAAACIkVAZibuvqLaoPAWxNHqZkQxtK6NmGwAAAIFExtleYWZHS3Izy5R0uaQPUhtW48RoJAAAAIiVSM/2hZIukXSgpP9Kygtfo5rm1GwDAAAgRiKT2qyVdM4uiKXRo2YbAAAAsRIZjeR2M9vLzDLN7GUzW2Nm5+6K4BobxtkGAABArETKSAa7+0ZJp0haLukwSdekMqjGKjOSwXTtAAAAiEok2a4sNTlZ0l/dPb1zojdgzZtRsw0AAIDtEhmN5O9mtkzSVkkXmdl+kkpSG1bjlBnJUFkFZSQAAAAI1Nmz7e7XSTpaUr67l0raLOnUVAfWGAXjbNOzDQAAgEAiPduS1F7SCWbWMmbZoymIp1FjNBIAAADEqjPZNrObJA2Q1FXSbEknSpovku0dMM42AAAAYiVyg+QPJR0v6St3HycpV1J2SqNqpIKebWq2AQAAEEgk2d7q7hWSysxsL0mrJR2c2rAap8xmGdpGzzYAAABCidRsLzKzvSX9UdJiScWSXktpVI1UZc22u8vM0h0OAAAA0iyR6dovDp8+aGb/kLSXu7+T2rAap+YRk7tUXuFqFiHZBgAA2N0lMl27mdm5Znajuy+XtMHMClIfWuOTGQkuJ6UkAAAAkBKr2X5AUj9Jo8PXmyTdn7KIGrHmzYLLWVrGTZIAAABIrGb7KHfvZWZvSZK7f2NmzVMcV6PUollEkvRtWbmkzPQGAwAAgLRLpGe71MwiklySwunaqZOIo0XYs/0ts0gCAABAiSXb90qaKWl/M7tFwYQ2v0ppVI1Ui8zgcpaUlqc5EgAAADQEiYxG8riZLVYwsY1JGuHuH6Q8skZoexkJPdsAAACoJdk2s31iXq6WND12nbuvT2VgjdH2MhJ6tgEAAFB7z/ZiBXXasQNGV752SYemMK5GKZpsl9KzDQAAgFqSbXc/ZFcG0hS0yKSMBAAAANvVZ1KbSeHr7zOpTXyUkQAAACBWfSa1OTt8zaQ2NWDoPwAAAMRiUpskipaRULMNAAAAMalNUlFGAgAAgFhMapNElJEAAAAgFpPaJFHLsIyEGSQBAAAg1ZFsh+Uj77v7kZKW7ZqQGq9mGaYMo2cbAAAAgVrLSNy9XNKHZvb9XRRPo2ZmatEsQrINAAAASYmNRtJG0vtm9oakzZUL3X14yqJqxFpkZuhbykgAAACgxJLtSSmPoglp0SyDnm0AAABISuwGyX/tikC+CzMbJmnYYYcdlu5QKCMBAABAVCJD/zV47j7L3cdnZ2enO5SwZ5syEgAAADSRZLshCWq26dkGAABAAsm2mV2eyDIEKCMBAABApUR6tsfEWTY2yXE0GZSRAAAAoFKNN0ia2WhJZ0s6xMyejVmVJWl9qgNrrFo0y9CmkrJ0hwEAAIAGoLbRSF6VtErSvpLujFm+SdI7qQyqMQvKSOjZBgAAQC3Jtrt/LulzSf12XTiNX8vMDJVwgyQAAACU2A2Sp5nZx2ZWZGYbzWyTmW3cFcE1Rns0j2grM0gCAABAic0gebukYe7+QaqDaQr2yGymrdtItgEAAJDYaCRfk2gnbo/mGdpaWi53T3coAAAASLNEerYXmdkMSU9L+rZyobv/LWVRNWKtmjdTeYVrW3mFWjSLpDscAAAApFEiyfZekrZIGhyzzCWRbMfRMjNIsEu2kWwDAADs7upMtt193K4IpKlo1TxIsLeUlilbmWmOBgAAAOmUyGgkR5jZy2b2Xvg6x8xuSH1ojdMeYc82N0kCAAAgkRsk/yhpoqRSSXL3dySdlcqgGrM9Knu2SbYBAAB2e4kk263c/Y1qy5iPvAbRnm3G2gYAANjtJZJsrzWzTgpuipSZ/VDBNO6Io7JmmzISAAAA/H/27js8qmr7//h7Tya9kRB67x2kCNKbBVTsFSuoXLtY0Xu/P6/eZr2Wa0cUwa6IIgjY6L1XAem9hRJqSJn9+2MyIYSUmWQmE5LP63l8JJNh9j7JYWadddZe25tuJA8Aw4GmxpidwGbg1oDO6hzm6UaiMhIRERER8aYbySbgQmNMNOCw1h4N/LTOXZ7MdqrKSERERETKvUKDbWNMBeB2oC7gNMYAYK19OKAzO0dpgaSIiIiIeHhTRjIRmAesBFyBnc65LyrU/SPVAkkRERER8SbYjrDWPhbwmZQREWHuNacn09SwRURERKS886YbyafGmHuMMdWMMYme/wI+s3NUWIiDEIdRZltEREREvMpspwGvAH8jq/1f1v/rB2pS5zJjDJGhIarZFhERERGvgu3HgYbW2uRAT6asiAwLUTcSEREREfGqjGQDcCLQEylLlNkWEREREfAus30cWGaMmQqc8jyo1n/5iwoL0Q6SIiIiIuJVsB+pTTMAACAASURBVP1D1n/ipYjQEC2QFBERERGvdpAcZYwJAxpnPbTOWpse2Gmd25TZFhERERHwbgfJXsAoYAtggFrGmDustTMCO7VzV2RoCIdP6HpEREREpLzzpozkv8DF1tp1AMaYxsCXQPtATuxcpm4kIiIiIgLedSMJ9QTaANbaP4HQwE3p3KduJCIiIiIC3mW2FxljRgCfZX19C7AocFM690WHOzmu7dpFREREyj1vgu37gAcAT6u/mcC7AZtRGRAT7uT4qQystRhjgj0dEREREQkSb7qRnAJey/pPvBAd7sRl4WR6JlFh3lzPiIiIiEhZVGjNtjHmcmPMUmPMQWPMEWPMUWPMkZKY3LkqJjwEgGOnVEoiIiIiUp55k3Z9A7gGWGmttQGeT5kQHe7+sR4/lQmxQZ6MiIiIiASNN91ItgOrFGh7LyY72FZmW0RERKQ88yaz/RQw0RgzHTjledBaqxrufHiC7aOpCrZFREREyjNvgu1/A8eACCAssNMpG6KV2RYRERERvAu2q1trWwZ8JmVITERWsK1e2yIiIiLlmjc12xONMRcHfCZliMpIRERERAS8C7bvAyYbY06q9Z93VEYiIiIiIuDdpjZqXuejqNAQjFGwLSIiIlLeeZPZFh85HIboMCdHFWyLiIiIlGsKtgMkOjxEmW0RERGRck7BdoBEhzvdO0iKiIiISLmVb822MSaxoL9orT3o/+mUHbHhKiMRERERKe8KWiC5GLCAyeN7FqgfkBmVEe7MtoJtERERkfIs32DbWluvJCdS1kSHOzl4/ESwpyEiIiIiQeTNDpIYYxKARri3bAfAWjsjUJMqC2LDnRxTZltERESkXCs02DbG3A08AtQElgEXAHOBPoGd2rktWsG2iIiISLnnTTeSR4Dzga3W2t5AW+BwQGdVBsREODmWmoG1NthTEREREZEg8SbYTrXWpgIYY8KttWuBJoGd1rkvLiKUDJflZLra/4mIiIiUV97UbO8wxlQAfgB+NcYcArYGdlrnvvjIUABSTqYTFeZVabyIiIiIlDGFRoHW2quz/vicMWYqEA9MDuisyoC4SPeP9sjJDKrFB3kyIiIiIhIU3iyQrJ3jy81Z/68KbAvIjMoIT2b7SGp6kGciIiIiIsHiTX3DT5ze3CYCqAesA1oEcF7nvLiIrDKSEwq2RURERMorb8pIWuX82hjTDrg/YDMqI+KU2RYREREp97zpRnIGa+0SoFMA5lKm5FwgKSIiIiLlkzc124/l+NIBtAN2BWxGZURsxOkFkiIiIiJSPnlTsx2b488ZuGu4v/P3RIwx9YG/AfHW2uuyHnMA/wTigEXW2lH+HjdQQkMcRIeFKLMtIiIiUo55U7P9fFFf3BjzMXA5sM9a2zLH4/2AN4EQYIS19kVr7SbgLmPMmBwvcSXubeIPADuKOo9giYsMVc22iIiISDlWaM22MaaxMWa4MeYXY8wUz39evv4nQL9crxcCvAP0B5oDNxtjmufz95sAc6y1jwH3eTlmqREXEcoRZbZFREREyi1vyki+Bd4HRgA+7T1urZ1hjKmb6+GOwIasTDbGmK9wZ7D/yOMldgBpWX8+5/Y9j48MVRmJiIiISDnmTbCdYa19z49j1gC25/h6B9DJGFMR+DfQ1hjzjLX2BWAs8JYxpjswI68XM8YMAYYA1K5dO6+nBE1cpJOdh1ODPQ0RERERCRJvgu3xxpj7ge+BU54HrbUH/TkRa+0B4N5cj50A7irk7w0HhgN06NDB+nNOxRUXGcqa3UeDPQ0RERERCRJvgu07sv7/ZI7HLFC/iGPuBGrl+Lpm1mNljmq2RURERMo3b7qR1PPzmAuBRsaYeriD7JuAgX4eo1SIjwzl6KkMMl2WEIcJ9nREREREpIR5k9nGGNMFqJvz+dba0V78vS+BXkCSMWYH8Hdr7UfGmAeBn3G3/vvYWrva96mXfp4t24+mplMhKizIsxERERGRkubNDpKfAg2AZZzuCGKBQoNta+3N+Tw+EZjo/TTPTXFZu0imnFSwLSIiIlIeeZPZ7gA0t9aWqsWH54L4rMy2tmwXERERKZ8K3dQGWAVUDfREyiJPGYl6bYuIiIiUT95ktpOAP4wxCziz9d8VAZtVGZGQVTpy+GRaIc8UERERkbLIm2D7uUBPoqxKiHZntg8dV7AtIiIiUh550/pves6vjTHdgJuB6Xn/jZJnjBkADGjYsGGwp3KGCpHuzPahEyojERERESmPvKnZxhjT1hjzijFmC/BPYE1AZ+Uja+14a+2Q+Pj4YE/lDGFOB7HhTg4qsy0iIiJSLuWb2TbGNMadwb4ZSAa+Boy1tncJza1MSIgO49AJBdsiIiIi5VFBZSRrgZnA5dbaDQDGmEdLZFZlSEJ0mDLbIiIiIuVUQWUk1wC7ganGmA+NMX0B7Tnuo8SoUGW2RURERMqpfINta+0P1tqbgKbAVGAoUNkY854x5uKSmuC5LiE6jEPHtUBSREREpDwqdIGktfa4tfYLa+0AoCawFBgW8JmVEYlRqtkWERERKa+86kbiYa09ZK0dbq3tG6gJlTUJ0WGcSMskNT0z2FMRERERkRLmU7AtvkuM9vTaVnZbREREpLxRsB1gni3b1ZFEREREpPxRsB1g2ZltLZIUERERKXfKRLBtjBlgjBmekpIS7KmcJTE6FICDKiMRERERKXfKRLBdWrdrh9NlJIdURiIiIiJS7pSJYLs0i48MxRgtkBQREREpjxRsB5gzxEFcRKgWSIqIiIiUQwq2S0DFmDAOKNgWERERKXcUbJeASjHh7D96KtjTEBEREZESpmC7BCTFhpOsYFtERESk3FGwXQIqxYSz/5iCbREREZHyRsF2CagUG87R1AxS0zODPRURERERKUEKtktApZhwANVti4iIiJQzCrZLQKVYd7CdrFISERERkXJFwXYJSFJmW0RERKRcUrBdAjyZbS2SFBERESlfykSwbYwZYIwZnpKSEuyp5KliTBgAyUe1sY2IiIhIeVImgm1r7Xhr7ZD4+PhgTyVPoSEOEqJC2X8sNdhTEREREZESVCaC7XNBpVjtIikiIiJS3ijYLiGVYsNJPqYyEhEREZHyRMF2CUmKUWZbREREpLxRsF1CKinYFhERESl3FGyXkEqx4ZxMz+T4qYxgT0VERERESoiC7RKijW1EREREyh8F2yWkcpw72N57RO3/RERERMoLBdslpFp8BAB7FGyLiIiIlBsKtktIlbisYDtFwbaIiIhIeaFgu4TERoQSE+5UZltERESkHFGwXYKqxkcosy0iIiJSjijYLkHV4iPYrWBbREREpNxQsF2CqsQpsy0iIiJSnpSJYNsYM8AYMzwlJSXYUylQtfgI9h87RUamK9hTEREREZESUCaCbWvteGvtkPj4+GBPpUBV4yPIdFmSj6UFeyoiIiIiUgLKRLB9rvD02t6dcjLIMxERERGRkqBguwR5em1rF0kRERGR8kHBdgmqFh8JoI4kIiIiIuWEgu0SlBAVSpjToY4kIiIiIuWEgu0SZIxRr20RERGRckTBdglTr20RERGR8kPBdgmrUSGSnYfVjURERESkPFCwXcJqJkSy50iqNrYRERERKQcUbJewmgmRZLqs6rZFREREygEF2yWsZkIUgEpJRERERMoBBdslrEYFd6/tHYcUbIuIiIiUdQq2S1i1ChEYAzsOnQj2VEREREQkwBRsl7BwZwhVYiOU2RYREREpBxRsB0HNhEh2KtgWERERKfMUbAdBzYRIdhxWGYmIiIhIWVcmgm1jzABjzPCUlJRgT8UrNRIi2X1YvbZFREREyroyEWxba8dba4fEx8cHeypeqZkQRYbLsvfoqWBPRUREREQCqEwE2+eamglZ7f8O+lZKMmLmJu4cuSAQUxIRERGRAFCwHQSejW187UiyZvdRpq3bzy5tiCMiIiJyTlCwHQQ1KkTiMLDVx8y2xQIwc/3+QExLRERERPxMwXYQhDkd1EiIZOuB4779RXeszYz1yf6flIiIiIj4nYLtIKlbMZotyb4F21mxNrM3JJPpsgU+V0RERESCT8F2kNSpGMWWAz6WkVh3gH34RDord54bbQ5FREREyjMF20FSt2I0KSfTOXwizeu/Y4G4CCfGwMw/VbctIiIiUtop2A6SuhWjAXzObidEh9GyejwzVbctIiIiUuop2A6Suknu9n++LJK0FgzQvVESS7Yd4mhqeoBmJyIiIiL+oGA7SGomRGEMbEn2PrNtAWMM3RtVIsNlmbvxQOAmKCIiIiLFpmA7SCJCQ6geH8kWnzLbFgO0r5NAVFgI01S3LSIiIlKqKdgOIndHEh+CbQDj7tPdvVESU9fuy+5QIiIiIiKlj4LtIKqbFM1WXxZIZtVsA/RtVoXdKan8sftIQOYmIiIiIsWnYDuI6laM4uDxNFJOerfQ0WIxxh1u925SGWPg9zX7AjlFERERESkGBdtBVCer/Z+3HUlsjsx2pdhw2tSswO9rFWyLiIiIlFYKtoOoQSV3sL1x/zGvnm8tGHP6675NK7N8+2H2HU0NxPREREREpJgUbAdRnYrROB2GDfu8DLaxGE5H232bVQFg2lp1JREREREpjRRsB1FoiIPaFaPYuM+HMpIcme1m1WKpHh/Bb2v2BmiGIiIiIlIcCraDrGGlGDZ4W0aS62tjDH2aVWbm+mRS0zP9PzkRERERKZYyEWwbYwYYY4anpKQEeyo+a1g5hi3Jx0nPdBX6XHdm25zxWN9mVTiZnsmcjcmBmqKIiIiIFFGZCLatteOttUPi4+ODPRWfNawcQ4bLetlv22JyPdKlQUViw51MWrknENMTERERkWIoE8H2uaxh5RgArxZJ5q7ZBgh3htC3WWV+XbPXq+y4iIiIiJQcBdtB1qCSO9j2pv2f5exgG6Bfy2ocPpHO/E0H/Tw7ERERESkOBdtBFh3upFp8BBu9ymyf2frPo2fjSkSGhjBp1e5ATFFEREREikjBdinQsLJ3HUnyy2xHhoXQu2klfl69l0xX7p4lIiIiIhIsCrZLgQaVYti47xjWFhwo59yuPbd+LauRfOwUi7ce8v8ERURERKRIFGyXAo2qxHA8LZOdh08W+DwLeae2gT5NKxPmdKiURERERKQUUbBdCjStGgvAuj1HC3yeu2Y7bzHhTno0SuLnVXsKzZAXJtNl+c/ENew45E07QhERERHJj4LtUqBxFXewvbaQYBvyTWwD0L9lNXalpLJkW/FKSbYdPMHwGZt4YeLaYr2OiIiISHmnYLsUiI0IpVZiZKHBdmEJ64tbVCHc6eDHZbuKNR9X1kA/rdzN6l3n3q6cIiIiIqWFgu1SokmVONbuPlLgc2weO0jmFBsRyoXNqjBhxW4yirHBTc6g/vVf1xf5dURERETKOwXbpUSzarFsSj5Oanpmvs9x7yBZULgNV5xXnQPH05i98UCR5+Kp+W5aNZbf1uxl+fbDRX4tERERkfJMwXYp0bRqHJkuW+C27QW1/vPo1aQSsRFOxi3bWeS5eBLbg7rWJSEqlNd+/bPIryUiIiJSninYLiWaVvNukWQhiW3CnSFc2rIaP6/aw8m0/LPkBfHUbMdGhPKXng2Y/ud+Fm3RVvAiIiIivlKwXUrUrRhNuNPBuj351227a7YLy23DledV53haJr+v3Vukubiyyr0NcHvnOiTFhPHfX5TdFhEREfGVgu1SIsRhaFwltsDMtrUUXkcCdKpfkcqx4YwrYlcSm1VIYowhKszJfb0aMnfTAWZvSC7S64mIiIiUVwq2S5GmVWNZs7uAYBuvYm1CHIYBbaozbd0+Uk6k+zwPTzcSR9Zgt3SqTfX4CF6ctBaXq3gb5oiIiIiUJwq2S5Gm1eJIPnaK/UdP5f0EW3jNtsdV59UgPdMyfoXv2W1PsO3pfBIRGsLjFzdh5c4UJqzUdvAiIiIi3lKwXYo0y1okuSafftve1mwDtKwRR5MqsXy7eIfP8/AskHTkGOqqtjVoVi2OV35ey6mMoi28zMuKHYf5ZfUev72eiIiISGmiYLsUaVE9HoCVO/PetdH6kNk2xnB9h5os336YP/cWvg38GeNkv8bpx0Ichqf7N2X7wZN8Nm+bT69XkOEzNnHvZ4tZsFndTkRERKTsUbBdisRHhlKnYhSr8gu28T7YBnc22ukwfLtou0/z8GS2c2+g06NREt0aJvHWlPWknPS9FjwvmS6Ly8LQr5YWqb5cREREpDRTsF3KtKwRX0Bm2/syEoCkmHD6NK3M90t3ku7D9u3ZNdu5HjfGnd0+fCKd96dv9Pr1CuKylrgIJ/uOnuLpsSuyd68UERERKQsUbJcyrWrEs+PQSQ4dTzvre75mtgGu71CL5GNpTFu33+u/Y7Nrts8erGWNeK5uW4OPZ21m1+GTvk0mDy4LNROieOKSJkxatYevFvqWhRcREREpzRRslzKtarjrtlftOju7XZSkb68mlUiKCeMbH0pJ8qrZzunxixtjLbz68zrfJ5R7LGtxOGBI9/p0a5jE8+NXs2GfbzXmIiIiIqWVgu1SpmUBiyTdmW3fUtuhIQ6uaVeTqWv3kXwsn5aCuXh6aeeV2QZ3Jvqu7vUYu3QnS7cd8mk+Z41l3eM4HIbXbmhDVJiTB79YSmq6/zqegPuYvlm43W+15iIiIiLeULBdysRHhVI7MZ9Fktb6ULF92vXta5Lhsoxd4l0bwOzMdgHPeaB3QyrHhvPc+D+KtdFNpstmX0BUjovg1etbs3bPUV6ctLbIr5mXP/cd5anvVnDXJws5mebfQF5EREQkPwq2S6FW+SySLErNNkCjKrG0r5PAF/O3eRUY59eNJKeYcCfD+jVl+fbDjF260/dJ5RgrZz/vPk2rMKhrXT6Zs4Wf/dh/Oz3DfUyLth7iwS+WkOHDglERERGRoioTwbYxZoAxZnhKSt5dPM41LWvEs/3gSQ6fOHORpLXebdeel1svqM2WAyeYs/FA4U/OtV17fq5uW4M2tSrw0uS1HDuVUaR5WXt2ucrT/ZvSumY8T3y7nG0HThTpdXPzXEBc1LwKv6/dxzNjV6rziYiIiARcmQi2rbXjrbVD4uPjgz0Vv/Asksyd3bZYn2u2Pfq3rEZCVCifzdta6HNdubZrz4/DYXhuQHP2Hz3F21M2FGleuTPbAOHOEN4Z2A6HMdz3+WK/1G9nZgXWt3SqzSN9G/Ht4h28NLn4CzzzMmHFLt6dtkHBvIiIiJSNYLusaVUzHmNg+fbDZzxenMx2RGgIN3Soxa9r9rInJbXA51rO3q49P21rJ3BNO3crwC3Jx32el8vmfQFRKzGK125ow+pdR/jHhD98ft3ccrYzHHphI27pVJv3p29kxMxNxX7t3L5dtIOXJ6/jqTErVK4iIiJSzinYLoXiI0NpWCmGJdvyCLaLGm0DN3esTabL8nUhvaxPZ7a9e92n+zUlNMTwr598D4rd3Ujy/l7fZlW4r1cDvpi/je+Xere4Mz+emDfEYTDG8I8rW3Jpq6r866c1xX7t3FzWEuZ08O3iHdz/+RK/d1bJSYs9RURESjcF26VU29oVWLrt0BmlCFnLFov8mnWTouneKIkvF2wrMOPqzQLJnCrHRfBQ30b8tmYfv/6x16c5WWvzbTEI8PhFjelUL5G/jl3Fn3uL3n870+U5JvfXIQ7D6zeeR5cGFXny2xVMXbevyK+d11ita8Tz3IDm/PLHXgZ/srDINe0F2bj/GC3+PpkXJq7JPj4REREpXRRsl1Ltaidw6EQ6W3IsELTWFiuzDXDrBXXYcySV39cWEFzms117Qe7qVo/GVWJ47sfVnEjzPrB05bFAMidniIO3bm5LdLiT+z5bzPEiL8R0H1RIjrHCnSF8cFt7mlaL5d5PFzNvkxeLR72Q6bI4HIY7u9bj9RvbMH/zQW75cB4H89gVtDiSj57CZeGDGZu4a9TCgPYQP3DslNd92kVEROQ0BdulVNvaCQAs2XrmpjHFjLXp27QyVeMiClwo6Spgu/b8hIY4+NdVrdh5+CT/+937xZIuLy4gKsdF8NbNbdmcfJxh360o0sJDzwJJR66aldiIUEYN6kjtxCju+mRhsTfpAfcxeYL6q9vW5INb27N2z1Fu+GAuu1OKv8W9h+eYrm1Xk1nrk7n6ndls3H/Mb6+f0xPfLqf3K9OYsGJXQF5fRESkrFKwXUo1qhxDbLiTpdtPB3/FrdkGd6Z4YKfazFyfnO+26NbHmm2PjvUSub59TUbM3MS6Pd6VfBSW2fbo3KAiT17SlAkrdjN8hu+LGl3Z7QzPHqtiTDif3d2JpNhw7vh4AX/sOuLz6+eU6bKE5AjqL2xehdGDO7I3JZXr3pvLJj8FxK6sSqAbz6/FF/dcQMrJdK56ezZTC7prUURHUzM4eiqDB79YyjNjVwasVnzf0VRuHTGfyat2B+T1RURESpqC7VLK4TCcV7sCS7aeXiRpsZhi57ZhYKfahDkdjJy9Jc/vFyWz7fHMpc2IiXDy/35Y5VUG2ubR+i8/9/asz2Wtq/HS5LXM+HO/T/M6vQV93t+vEhfB53d3IibcyW0fzWfDvqIHxJn27Ax6p/oV+XLIBaSmZ3L9+3Pz3iHURxlZ0XaIw9CxXiI/PtSN2hWjGDxqIe9P3+jX1oMZLkvn+hW5t2cDvlywjSvfmVWsGvr8rN51hFkbkrn3syUM/WopKScCVxrz6bytxeoRLyIi4g0F26VY21oVWLvnSHadsj8y2wBJMeFcdV51vluy46yNc+D0du1FkRgdxjP9m7Jgy0HGLC68y4erkAWSORljeOW61jSuEstDXy5l6wHvWw16LiBCCojsayZE8dndnTDGcMuIeUXeUMdaS0gew7SsEc+393YmIjSEGz+Y6/MFQ265j6lGhUjG3NuFy1pV48VJaxn69TK/dULJdFkiQh083b8powZ35MCxNK54exZfLdjm16A+M9P9WgPaVGfCit1c9Pr0gGTqAb5asI33pm3k4tem8/sa3xb2+mJL8nEGvDWLz+dvDfhC1iOp6UVe1yAiIoGhYLsUa1snAZeFFTvcWdCibteel8Hd6pGa7uLLBWe3AbTFyGwDXN++Fu3rJPCfiWs4VMiiQJfL+64nAFFhTobf1gFjYMho7xdMZrq8O6b6lWL47O6OnMpwMXDEvCLVWOcuI8n9+mPv70KtxCgGf7KQsUuK3nYwu51hjmOKDAvhrZvb8uQlTfhx+S6uf98/deLuY3K/XfRsXIlJj3SnfZ0Enh67koe+XMrRVP9koDOyfk/39qzPDw90pUJUKIM+WciwMSv8NoZHpsvStGosMRFO7hq1iAc+X8K+IwX3oC+K1buOsHJnCn/7fhVXvzubZbn65/vT4JEL6frSFD6ZvZn0APZ4P3DsFC9MWpNvKZo/nUjL4FSGWlyKyLlLwXYp1rZWBQCWZC3a82cGsWnVOLo2rMioOVvO+lAuas22h8Nh+PfVLTmamsE/C+m9ndcOkoWpXTGKt29ux/p9R3lyzHKvfi4F1Wzn1rRqHKMHd+TwiXRu+XA++4/61oUj01XwTp9V4iL45t7OdKyXyGPfLOedqUXbbTIzRxlJTsYYHujdkA9v68Dm5OMMeGs2CzYf9Pn1zxzLEpLj3aJyXASjB3fiyUuaMGnVHi773yxW7Ch+EOm5KHI6HLSsEc/4h7pxX68GfLt4O/3emMnsDcnFHsMjw2VpUCmGCQ9158lLmvDrmr30fW06X8zfll125J9x3L+nxy5qzJ6UVK5+dzZPf7fC791pAPYfO8WJU5k8N/4P+r0xw68tLXOatm4/H0zfxMWvz+CZsSsDcpHiceMH8+j20lRGz91CWkbgLiA27DvK/Z8vZuq6fQHf/XXhloN+KSUrjMtl/X6RKiK+U7BdilWICqN+pegzOpL4o2bb465u9dhzJJVJq/ac8bgvgWl+mlaN475eDRi7ZGeBH/jWywWSuXVrlMQz/ZsxceUe3p22sdDnZ9ehe3nGt65ZgZGDzmd3invBni9t73J2I8lPXEQonwzqyBVtqvPKz+t4dtxqn0sMcm7Uk5cLm1fh+/u7EBvhZOCH8xg5e3ORg4hMa3Hm+uGFONxB/ddDLiAj08W1781hxMxNxQpUM3OVxoQ7QxjWrylj7utCuNPBLSPm8+y4VT61l8x3rKw7EGFOBw/0bsjPQ3vQsno8f/1+JTcOn+u3rK3n93rledX5/fGe3N2tHmMW76D3q9P4dJ5/S0syMi2Xt6nGiNs74LIwaORC7vh4Aev9XF/vuUC/ok11xizeTs9XpvHar38GpP59d0oqR06m8+y41fR9bRrfL90RkHKcORsPMHHlHgaNXMi1781h9obkgAXd9322mMvfmsVdnyz0y0VqfsYt30nr53/hgc+XsGZ38RZ+F+ae0Yt48IslxV5gXpjVu1L46/crA3qHyOPH5buYszFw54HH4RNpzN14oET2S9h1+GRAN1rzsNYG9OI49ziB/h0Vl4LtUu78Ooks3HIQlytrE3X/xdr0alyZ+knRfDxr8xmP+7Jde0Ee7NOQhpVj+NvYlfl+CLts/iUXhbm7ez2uPK86r/6yjilrC665za5v9iGwP79uIh/d0YGtB49zy4fzOeBlwF1QGUlOYU4Hb9x4Hn/pUZ9P523l/s8X+/QmeDowzf85jarEMu7BrvRqUpnnx//B0K+XFSlQ9fQOz0uHuolMfKQ7vZtU5l8/reGuUQu9/lmdPY77zdmZa6x2tRP46eHuDOpal9Fzt9L/zZks2lK8bH2Gy4UzR3F9vaRovrinE69c15r1+47R/82ZvP7rn8UuYcjIqkN3hjiIjQjlb5c1Z9Ij3WleLY7/98MqrnxnFou3Fr/lJLiPKdTh4MLmVfh5aA/+77JmLNl2iH5vzuTZcav8lk33lPv89bJm/PZYT/o2q8z/fl9Pr1em8uncs++WFW8sFzd0qMXIQecTGx7Ko18v59I3Z/LrH3v9+gGbnvV7D4WrygAAIABJREFU+uulTdmdksotI+Zz0/B5xb4rlJcTaZk0rhLD4m2HuOLt2Qz+ZCHLAxA87jqcirUwbd0++r85k3tGL2LljsBk1H9fs5cJK3Zz6f9mMviThX47p3P7efVevpi/javemc3AD+cxa31gguG0DBePfLWUgR/O56p3ZjNx5e6ABcMfz97CzR/Oo+9/p/HZvK0BC4ZPpmXS57/T6PriFN78bX1A7q55vDN1A62e+5lnxq4sVsOBwuw7eorG/zeJLxZsC9gY/qBgu5TrVD+RI6kZrNt7FKxfY20cDsOgrnVZtv0wi7ee/kDxdbv2/IQ7Q3jp2tbsPpLKy5PX5vmczGJs1GOM4cVrWtOiehwPf7mMtXvyz6ic3kHSt8G6NEzi4zvOZ+vB4wz0MuB25dGNJD8Oh+GZS5vx7OXu3SZvGTG/0Dr37HFcnmC74H/GcRGhDL+tPU9c3Jgfl+/imnfnsCXZ+8Wl4P755Q6Ac6oQFcYHt7XnH1e2YPaGA1zyxkymFaGEwROY5nWxEhkWwt8HtODLey4g02W5/oO5/PunP4r8wZSZefYxGWO4vkMtfnusJ5e1qsabv6+n/5szmV+MDY8ysktjTo/VqEosX9zTibdubsv+o6e49r05PPnt8mJvHJSRabMvIMKcDu7uXp/pT/bmlk61+Xz+Nnq9MpWPZm0udsbJswOt0+GgTsVo3h7Yjh8e6Er9SjH8v3GrueT1GUxetccvQVBm1jH1blKZCQ91462b25KW6eKe0Yu49r05zN3on82oPMd06wV1mPpEL/4+oDkb9x/nhg/mcttH8/2aSc3ItPRpWoWZT/XmyUuasGTbIa58x/9Bt+ff06xhfXj0wsYs2HyQAW/P4s6RC/waDLtcFpd13y19/KLGLN12iGvfm8NNw+f6PRhOz3ThdBj+dmkzNuw7xq0fzeeKt/0fDKdnurAWujVMIuVkOvd/voQLs8rM/B0MH01NJyzEQXxkKP/3wyq6vjiF//2+3uvPAm+dSMsgNd1FRGgIr//2J51f+J2/fb/Sb+1oc9px6CSZLst3S3Zw4WvTGfzJwoDcJfC8l4UWlHUqBUr37ISO9RIBmL/pQNYCSX+G23Bt+5pUiArl/emne1d7/jH4Y6z2dRK4s4s7G5lXhqioZSQekWEhjLj9fKLDQ7jrk0X51ld7/n0XJYvepWESH91xPlsOHOeWEYUH3JmuvLuRFGRwt3q8fXM7Vu5I4br357D9YOGdUDxBnDfZeofD8GCfRnwyqCN7jqQy4O1ZPnXg8CZbb4zh9s51GfdgVxKjQ7lz5EKeH7/at2y9K/9g26Nzg4pMHtqDmzvW5sOZm7n8rVlF2owoPceiz9ySYsJ546a2jBrckbQMFzcOn8fT360oUivCgmrrB7SpzpTHe/GXnvX5fulO+rw6jVFztmQHfr7KyOOiKDE6jH9c2ZJJj3SnTa0K/HOCu5779zVFzwxnX0DkONHPq1WBr4dcwEd3dMDhMNz72WKue3/uGRfyRZHucmV/kDoc7p/ZL4/24IVrWrHrcCo3fziP2z6aX+yM7emLIgcRoSEM6lqPmU/15q+XNmX1riNc9c5s7vpkYbFrra21pGW6CAsxxEaE8kDvhswa1ueMoHvQyAV+Ce4zXC6MgYToMB65sBGzhvXmqX5NWLEjhWvfm8PAD+cxd+OBYgdA6VnneGJ0GA/1bcSsYX34v8uasTn5OLd+NJ+r3p3DL6v3+GUtREamizCng3t61GfmsN68eE0rjqa6g+GLXpvO1wu3+aV8wXOh0qdpZX5/vBfv3tKO2Agnf/1+Jd1emso7Uzf4bdfe9EwXsRFOfnigK18NuYA2tSrw2q9/0uXFKfx93CqvPg+8G8d9TA/2achvj/XgmnY1+HbxDvq+Np0hoxexaMtBvwXD6ZmWKnERzHm6D0MvbMTy7YcZ+OF8Ln9rFj8s3em3u1+ef7dhCralOGomRFGjQiQLsv4R+DfUdnf3uLNLXX79Y292XWf2Akk/jfHkJU2omRDJsO9WnBV4FWWBZG5V4yMYcfv5HDyexpBPF+UZ3J3uRlK0Mbo2TOLjO89nc7I74C7o9ltBJRcFuax1NUbf1ZF9R09x9btzCs1wZfcO9+Ffcc/GlRj/YDf3jpmjFvHar3969QGY4XJ5XYLTrFocPz7YjTu71GXk7C1c+fbsAu86nDnO2VngvMSEO/nP1a0YNbgjx09lcO17c/j3T3/4tNlOYdl6cP+8fnm0B3/pUT/7Q2n88l0+fSB5PuBC8/lFRYc7eaZ/MyYP7UHrmhX4+4+rGfD27CKVyWRkunDm86HTuEosowd3ZOSd54OBu0Yt4raPFnj9u8kpv2MyxtC3WRUmP9KdF69pxfaDJ7j2vbn85dNFRd7dNCOPOxChIQ5u7libaU/24q+XNmXlzhQGvD2LBz5fUuRxPB/+oTkuICLDQhjSowEzsjLQC7cc5PK3ZnHvp4u93rgrt+xFwDl+TzHhzjOC7qXbD3OVH4LutEzXGRm/2IhQ7u/VkFnDevN/lzVj/b5j3PzhPG7IakVa1EDLcz54Ap7ocCd3d6/PjKd685+rW3HoeBpDPl1M/zdnMm7ZziJfTHrG8hxTuDOEmzrW5vfHe/HOwHZEhoUw7LuV9Hh5KiNmbipWG8y0HOdDiMNwaatqjHugK1/c3Ynm1eN45ed1dHnhd/790x/sSSneAuH0DPfdG2MMF9SvyMd3ns8vj/bgstbV+GLBNnq+MpUHvlhS7Pr+0+e4g4aVY3nhmtbMHtaHh3o3ZMGWg1z3/lyueW8Ok/xwlyA900VoiCEpJpyhFzZm9tN9eOGaVqSmZzL062V0f2kq70/fWOwLlpzHVJqV7tkJ4M5uL9h8EJef+mzndkfnukSGhvDedPdCw9M12/4ZLCrMyYvXtGZz8nHe+G39Gd/zpc92QVrVjOf1G89j6bbDPDnm7C3di7NRj0fOgHvgh/PyDbi9WSCZnwvqV2TsfV2ICHVw4/C5Be6kmHsxobdqJUbx3X1duK59Tf73+3oGj1qYZ7/1M8ZyQYgP6fqI0BCeu6IFIwedz4HjaVzx9mw+nrW50MDem8x2Tp5g+KasLHf/N2d4XfKRkenyapyoMCfPXNqMcQ90pVp8BA99uZQ7Ri70us979jEV8vNrWDmGT+/qyLu3tOPwiTSue38uj32zzKduOHlltnMyxtC7aWV+HtqDvw9ozsqdKVz65kz+9v1Kn+rss2vr8zkmZ4iDm7KC4ccvasys9clc/PoM/u+HlT4dj7XWfUz5fJBGhJ4Ohh/u05Cp6/Zx8eszGDZmBbsO+9byMiPTffcmrzt62cHw0314pG8jZm1Ipt+bM3joy6U+B/d53RU4a5xhfXiqXxOWZQXdd45cUKS7NxmZltA8zoeoMHcwPPOp3vzjyhbsPHSS2z9ewFXvzinSHY/ssqJcxxTuDGFgp9pMebwnb9x4HhbLI18to+9r0/lywbYirYfwBHE5hTgMl7WuxoSHujF6cEfqJkXxr5/W0OXFKbz2659Fqk3OK4gzxtClYRKjB3fkp4e7cWHzKnw8ewvdX57Ck98uL/Ki6px3bzwaV4nl1evbMPOpPtzToz4z1u3nirdnc/PweUXumJPXBWWl2HAeu7gJc57uwz+vbMHB42nc9/kSer86jdFztxR5MXp6rgu9iNAQbu5Ym18f7cnIO8+nfqVoXpy0li4v/M7z41cXOXt/uowkAMGRHynYPgd0qpdI8rE0th084ffMNrhvMd7csTY/LtvFjkMnsrcB92dg361REjd2qMWHMzedkanxtc92Qfq1rMqwfk0Zv3wXb/5+dlAP3tdS56drVklJQRnu4iz6BHc97/f3d3V3dPl8CcNn5L0bZIaPgWlOEaEhvHJda/51VUtmb0hmwNuzWL0r/9vjmS5XoVngvPRuUpnJQ7vTvWES/5jwB3d+srDANnE5b+V7KzYilP9c3Yov7umEy8KNw+fx7LhVhXbG8CaznVPLGvH88EBX/j6gOUu2HuKi12fw5m/rCw0Y0vNZ9JkXY9wZtN8f78n9vRowfvku+rw6jREzN3l1a9wdmBY+TmiIg0Fd6zH9yV7c3rkuXy3cTq9XpjF8xkavAiBPJrOwY4oKc/JQ30ZMf8pdN/7Vgu30emUqb/623quso+d8yCtgzCkuIpTHLm7CjKd6c3vnOny/dCe9Xp3GPyf84fVFRLoX53hcRCiPXtSYmU/15t6eDfjtj71c9Np0Hv9mudebYGUHPAWc4zHhTu7v1ZCZWUH38u2HufrdOT4H3QXd6QD3+8Dtnesy7cnevHBNKw4eP8VdoxZx2f9mMWnlbq/LPtKyg+28x3KGOLiqbQ0mP9KDD25rT3xkKM+MXUnPl6fx0azNPgV0uYO4nIwx9Ghcia+GdGbs/V3oWC+R//2+nq4vTuH58at9ugDzlJHkN1aL6vG8eVNbpj3Ri4EdazN+xS4ufG0G94xe5HM9fHqmzbcMomp8BM/0b8acZ/rw10ubsjn5OINGLqTfGzMZs3iHTyUz6QUcU1SYk9s612XK4714/9Z2VIwJ49lxq+ny4hRe/Xkd+476lr3P7/fkcLgv+L+45wImPNSNi1tU5dO5W+n5ylTu/3xxdqtjX8bJ75hKk9I9OwFO122D/2u2Pe7uXg+AETM3Z+8g6a/MtsdfL2tGldhwHvvm9M6GvmzX7o17e9bn+vY1eeO39YxbtjP7cc9nRlEzzjl1a5TEiDs6sGn/MQZ+OO+sRW2ZfriAqBQbzldDLuDSltX4z8S1/PX7VWfVuGUvkCziWMYYbr2gDl//pTPpGZZr3p3D1wvz3hEy01X0OxBJMeGMuKMD/7yqJQs2H6BfVieJvGTXNxchS9GlQRKTh3ZncNd6fDpvK5e8PoOZ6/PfpTPDZX0eJ8RhGNS1Hr8/3pOLmlfh9d/+pH8h/b8zvQxMc4oKc/JUv6b8PLQH7esm8K+f1tDvzRkFLjq11p6x+ZA3KkSF8dwVLfh5aHc61E3gPxPXcvHrM5i0cneBmbMMlyvfLHBekmLC+ceVLfnl0R70aFyJ13/7k16vTuOL+dsKLCnI2cnF23H+PqAFU57oyRVtqjNy9mZ6vDyV13/9s9B+0xk5yhMKkxAdxrB+TZk5rDeDu9Zjwopd9PnvNJ4Zu4KdhQR02RcqXpx7OYPuYf2aZgfdd3y8wKugJN3l3TGFOd1lOVMe78Wr17chNT2T+z5fQr83ZzBu2c5CSwoysstICj4mh8NwSYuqjHugK5/e5c5A/3PCH3R7aSpvT1nvVUlBzkXABWlXO4EPb+/Ar4/2oH8rd0DX4+WpPPHtcq+6Y6Tlk63PrVZiFM9f2ZLZw9x3PRZuOci1783h+vfddwm8uWBJz3AVOk5sRGj2XZz/Xt8GgCe+XU6Pl6cyfMZGr/qpexOYhjgM/VpW4/v7u/LdfZ3pVC+Rd6ZtoNuLUxk2ZoXXbUTTMy2hzoLPvZY13HekZw3rw5AeDZi1Pplr3p3Dte/NYfIq70pZCrqAKE1K9+wEcLckS4oJB/zbjSSn6hUiuaptDb5auK3IbdsKEx8ZyivXt2HT/uO8lNWdxFXMBZK5GWP499Wt6FgvkSfHrMjOMBS3Zju37o0q8fGd7kWTN34wl705srXuzHbxx4gIde8GeX+vBny5YBuDP1nIkRxvqL6WXOSnXe0EJjzcjQ51Exj23Uqe+HbFWZkmX7PAuRljuO2COkx4qBtV4yK4Z/Qi/vb9yrNqrD1xV1HHigpz8uyA5oy5tzPhoQ5u+2gBT41ZnueHeHGOqUpcBO8MbMfowR3JtJZbRsznka+W5pn9Kc4diPqVYvhkUEc+vrMD1sKdIxdy96iFeXaT8TYLnJeGlWMZOagjowZ3JNzp4L7Pl3DDB3PzrRfOq47aG/UrxfDere357r4u1EmM4q/fr+SSN2bw8+q8O5d47gr4eou4ZkIUr17fhp+H9qB7o0q8+ft6er7ivkOQ34JddxbYt3GSYsL5v8ubMyMrc//d4p30fmUaz45bdcZ7Qu5x3Mfk/ZtETLiT+3o1YFZW0L1ix2Gu8SLoTs84u+SiIKEhDq5rX5NfH+vJmzedB8AjXy3jotemM2bxjnwXtWVfFHl5oWeMoXsjdwb6u/s6c16tCrz6y590e3EKL09eW2BXntx16IVpVCWW1244j2lP9uLWC+owYcUuLnp9On/5dFGB62I8m1F5u/CuYkw4j17UmDlP9+HvA5qz63Aqd41aRL83ZxSagc7Io4wkP2FOB9e2r8nkod0ZOeh86iVF85+Ja+nywhRemLimwPrxtDzKSArSvk4iH9zWgSmP9+KG82vyw7KdXPT6DAZ/srDQhbXpWYuAvVE1PoKn+zdl7jN9eW5Ac/YdTeXez9ylLJ/M3lzgXbC8SmNKIwXb5wBjDJ082e0Ank/39qxParqLkbO3AMUvuchL14ZJ2Qvn5mxIdtds+/ksDHM6+ODW9lSPdwd1W5KPn96C3o/H1LVhEqMHd2JPSio3fDCXHYfct5Hd3Uj8M47DYXiqX1Nevq41czce4Np3T3cq8VewDe6gYfTgTjzStxFjl+7gqndmn5H9ySxmaYxHw8qxfP9AF/7Soz5fLNjGZW/NPKOLRH6dO3zVvk4iEx/uzn29GvDdkp1c/Pp0fsuRTc+uBS7mydejcSV+HtqDh/s2YtLKPfT973RGz91yRkYmI6s8oTh3O/o0rcLkod3dH0gbD3Dx6zN4cdLaM0pl8lp456uejSsx8eHu/OfqVmxOPs5V78zm4S+XZp/bHuk+ZIHz0r5OAt/e25nht7XHAn/5dDHXvjfnrI5FGUW4K5BToyqxvH9be8Y90JXm1eL4109r6PPqNL5asO2swDEts+jnQ5W4CJ6/siVTn+zFte1r8sX8bfR4eSr/nPDHWYFjuueiqAjBQXSOoPvp/u6Fode8O4fbP867lZ+3ZUW5hTgMV57nLvt4/9Z2RISG8MS3y+nz32l8ueDsbh/eZoHz0r5OIh/feT4/PdyNHk0q8d70jXR7aQrP/Zh32Ye7Dt3331PNhCieu6JF9oLAuRsPcGUBvbrTM4r27ykqzMmgrvWY9mQv3rjxPBzG8MS3y+n5Sv6LNtMy81+XkB9j3O0wvxxyAT8+2JWeTSrx4cxNdH95Co9/szzPBbyFlcbkp15SNP+6qhVzn+nLYxc1Zvn2w9z84TyueHs2Py7flefdqYLKffITHe7kzq71mPZEb967pR1JMWE8N/4POr/wOy9NXpvnhUT2BUQhWfRgK92zk2yd6ruD7R0HfVv044uGlWO5uHkV9mRlZAIV1w/r15T6SdE8OWYFJ9MzA1IakxAdxshBHbHWcufIBew/5q6t9ndpTMd6iXx2dycOHU/jxg/msfXAcVxF7EZSkBs61GL04I7sPeLe7nvptkNFXiCZnxCH4dGLGjNqUEeSj6VxxduzsktxvN2oxxvhzhCeubQZn9/ViROnMrn63dm8PWU9GZkun9oZFiYi1L375Pf3dyEhKoy7Ry/ika+WcvB4Wo5t4f0zzmMXNWby0O60qVmBZ8et5up3Z2dfRGT46WcX7gzh3p4NmPpELwa0qc770zfS59VpjF2yA5fLZgePxT0mZ4iDgZ1qM+3J3jzYuyE/r95Dn/9O58VJa7PvrGS6vFtcWhBjDBe3qMovQ3vw4jWt2Hn4JDd8MJfBnyzM7pCSUUgtsLfa1KrAZ3d34ou7O1EpLoKnx67kotemM27Zzuzb/Bl5LLzzVY0KkbxwTSumPO7+HY2cvZnuL03lxUlrs9d35OxRXlTR4U7u7dmAmU/15un+TVm1093KL3fQXZSAJydHVknBTw9346M7OpAYHe6utX5lKqPmbMm+S+BrFjgvLarH887Advz2WE8GtK7OZ/PcdbzDxqxgc447OelFuAORU8WYrAWBz/Tlb5c2Y+P+vHt1F/WuikdoVp36pEfcGeg6FU8v2nz153VnXISlZ3ifBc5L65oVeHtgO6Y90ZuBHWvz08pdXPLGDAaNXHBGBrq49c2J0WE83LdRdmeR42kZPPzl0uy7Rjkv/tOKcUEe4jD0b1WNsfd35bv7utCtURIfTN9I95en8NjXy87YpTQjVyec0irkueeeC/Yc/Gb48OHPDRkyJNjTCIiosBA+nbeVnYdPMvTCxgEbp3ZiFF8t3A7APd3rExPh9PsYoSEOWteM56NZmzmV4aJNrXj6NK3i93ESosI4v24io+Zsye7mcn+vhkSEhvh1nGrxkXRvVIlvFm3nuyU7OXYqg7a1EujRuJJfx6mVGMVFzavy08rdjJrj3mVs5+GTPNSnUbGDkZzqVIzmyvNqsHDLQT6evYXkY6dYviOFzvUr0rlBkt/GqZUYxfXta7Hj0ElGztnCrA3JuCz8ufcoj17kv3O8SlwEN3SoRYjD8Pm8bXyzaDtV4iL4efVeujdO4vy6iYW/iBcSosO4um0N6leK5qcVexg5ezMHj6dxLDWTDfuO8UDvhn4ZJzrcySUtqtKjcSUWbz3E6LlbmbF+P3UqRvPD0p30aVqZtrUTij1OmNNBl4ZJXNuuJslHT/HpvK18vXA7kWEh7D6cytaDxxnSo0Gxx3E4DC1rxHNrpzrEhDsZt2wnI2dvYfvBk9RKjOLbRTu4pEUVWtaIL/ZYtRKjuOn8WrSsEc/CLQf5bN42Jq/aQ5W4cNbtOcreo6kM7lqv2OPER4VycYuqDGhTneRjp/hiwTY+m7uVE2kZVImL4JtFO7isdTWaVo0r1jhhTgcd6iZy6wV1iI8MZdKqPYyas4XFWw9ROzGKRVsOceRkBrd1rlOscYwx1K8Uw03n16J9nQTW7jnC5/O38fWi7TiMu/b/28U7uKptDRpUiinWWInRYVzcoirXtKtBeoaLMUt2MHL2ZjbsP069pGimr08mI9NyU8faxRonzOmgfZ0Ebu9ch5oVIpm76QCfz9/GhOW7iAh1EBnqZOzSndzQoRa1EqOKPI4xhnpJ0VzXvhY9G1di75FTfLlwG6PmbGHvkVQaVIph0qo92eUhxREfFUrvppW5pVMdosOd/LJ6L5/O28rUdfuIjXBicG9BP7BTbarFRxZ5HGeIg1Y14rntgjq0qhHP+n3H+HLBdj6dt5WUk+k0qhzLD0t3Eh8ZyhXn1SjWMVWvEMllratzdduaZLos45bv4pM5W1i45SCJUWGcTM9k4so93N6lTna5bUl6/vnndz/33HPDC3ueKe37yfuiQ4cOdtGiRSU/8I5FkO6fpvP5sdYycMR8WlaP42+XNQ/oWDd/OA+Adwa2IzE6LGDjfLNoO98v3cnFzaswyA8fcPlZsPkAr2e1HBxxRweiw/x/AQGw7eAJ/v3TGo6kpnN56+rc0ql4Hwb5STmZzuu//uneVRT47K5Ofss655Thsny9cDsTVuwC4Np2NbmumB8G+Zm9IZmPZ2/mRFomIQ7DZ3d1Csg4Ww+c4IMZG7MzZQM71WZA6+p+H+fYqQy+XbSdX7JKV6LDnIy4o4Pfx3G5LNPXJ/PVgm3ZWedBXepycYuqfh9r0/7jfDZvK2uyss4JUWG8e0s7v49zLDWDcct3MXn1HjIzXVjgvp4N/H7x6nJZ5m46wJjFO9hzJBWHMVSJDee1G8/z6zgAOw6dYMySnczfdIAQhyHTZRnatxGd6lf06zgn0zP57Q/3tulHUtNxOhzUTHBn2/3JWssfu4/y/dIdrM6RZXyqX1Pa1qrg17EOn0hn4srd/LpmL6npmYSGOKifFM1zV7Tw6ziZLsvCLQcZt2wXWw4cJyzEQVqmi2cvb0GzarF+HWvn4ZNMXLmbGX8muzs9hThoVi2OZ/o39es4aRkuZqzfz4QVu9l7JJUIp4PUDBf/uboV9ZKi/TrWhn3HmLByNws2HcBhDMZhaFergl8TJ+B+b52ydh+TV+3h0Ik0IkNDOJmeyQtDrqFuXf8kNHxhjFlsrS30zV3Btj+80wn2570duYiIiIgEzomL/0tUl7tLfFxvg+3ApPiC5cABCEZZzJE24PLvVbaIN6y1AWsHGYxx5NzgSdKUpXOvJI/JX5t5eTOOoWSOyZ/rOgqiYyo6ay2Z1j9rVQrjbhdbcscUNX4N/PJcQMcqDmW2RURERER85G1mu3Qv3xQREREROYcp2BYRERERCRAF2yIiIiIiAaJgW0REREQkQBRsi4iIiIgEiIJtEREREZEAUbAtIiIiIhIgCrZFRERERAJEwbaIiIiISIAo2BYRERERCZAyEWwbYwYYY4anpKQEeyoiIiIiItnKRLBtrR1vrR0SHx8f7KmIiIiIiGQrE8G2iIiIiEhppGBbRERERCRAFGyLiIiIiASIgm0RERERkQBRsC0iIiIiEiAKtkVEREREAkTBtoiIiIhIgCjYFhEREREJEGOtDfYc/MYYsx/YGqThk4DkII0tpY/OB8lJ54PkpnNCctL5cG6qY62tVNiTylSwHUzGmEXW2g7BnoeUDjofJCedD5KbzgnJSedD2aYyEhERERGRAFGwLSIiIiISIAq2/Wd4sCcgpYrOB8lJ54PkpnNCctL5UIapZltEREREJECU2RYRERERCRAF2z4wxlQwxowxxqw1xqwxxnTO9f1expgUY8yyrP+eDdZcJfCMMU1y/K6XGWOOGGOG5nqOMcb8zxizwRizwhjTLljzlcDy8nzQe0Q5Y4x51Biz2hizyhjzpTEmItf3w40xX2e9R8w3xtQNzkylJHhxPtxpjNmf4z3i7mDNVfzHGewJnGPeBCZba68zxoQBUXk8Z6a19vISnpcEgbV2HXAegDEmBNgJfJ/raf2BRln/dQLey/q/lDFeng+g94hywxhTA3gYaG6tPWmM+Qa4Cfgkx9PuAg5ZaxsaY24CXgJuLPHJSsB5eT4AfG2tfbCk5yeBo8y2l4wx8UAP4CMAa22atfZwcGclpUhfYKO1NvemSlcCo63bPKCCMaZayU/wV2p7AAAEkklEQVRPSlh+54OUP04g0hjjxJ2g2ZXr+1cCo7L+PAboa4wxJTg/KVmFnQ9SBinY9l49YD8w0hiz1BgzwhgTncfzOhtjlhtjJhljWpTwHCV4bgK+zOPxGsD2HF/vyHpMyrb8zgfQe0S5Ya3dCbwKbAN2AynW2l9yPS37PcJamwGkABVLcp5SMrw8HwCuzSo7HGOMqVWik5SAULDtPSfQDnjPWtsWOA48nes5S3Bv3dkGeAv4oWSnKMGQVVJ0BfBtsOciwVfI+aD3iHLEGJOAO3NdD6gORBtjbg3urCRYvDwfxgN1rbWtgV85fddDzmEKtr23A9hhrZ2f9fUY3MF3NmvtEWvtsaw/TwRCjTFJJTtNCYL+wBJr7d48vrcTyJmZqJn1mJRd+Z4Peo8ody4ENltr91tr04GxQJdcz8l+j8gqLYgHDpToLKWkFHo+WGsPWGtPZX05AmhfwnOUAFCw7SVr7R5guzGmSdZDfYE/cj7HGFPVU2tnjOmI++erN82y72byLxn4Ebg9qyvJBbhvG+4uualJEOR7Pug9otzZBlxgjInK+r33Bdbkes6PwB1Zf74OmGK1AUZZVej5kGtNzxW5vy/nJnUj8c1DwOdZt4k3AYOMMfcCWGvfx/1GeZ8xJgM4CdykN82yLatu/yLgLzkey3lOTAQuBTYAJ4BBQZimlBAvzge9R5Qj1tr5xpgxuMuHMoClwHBjzD+ARdbaH3Evuv/UGLMBOIi73l/KIC/Ph4eNMVdkff8gcGew5iv+ox0kRUREREQCRGUkIiIiIiIBomBbRERERCRAFGyLiIiIiASIgm0RERERkQBRsC0iIiIiEiAKtkVEREREAkTBtojIOc4YY40xn+X42mmM2W+MmRDMeYmIiIJtEZGy4DjQ0hgTmfX1Rbi3ARcRkSBTsC0iUjZMBC7L+nO+W8Z7GGOeM8aMMsbMNMZsNcZcY4x52Riz0hgz2RgTmvW8Z40xC40xq4wxw42bM+uxXlnPecEY8+9AHpyIyLlKwbaISNnwFXCTMSYCaA3M9+LvNOD/t3e3vFUEYRiG7zdUQIIghaDxtA5sQwkC0SAIfwCCQ/QvtA6B6UeoIMHALyCEgkIhmjoOGCxIRC2U9EGcPYFUVcymyfS+5OzOl3vy7mwGbgP3gNfAxySLTK+SnwX37SQ3kywAF4CVJH+YXiO9U1V3gLvAesvNSFIvDNuS1IEkn4FrTKva707YbTfJITABzgHvh/bJMBbAclXtVdWEaTC/Psz3FXgFvAUeJfndYBuS1J25016AJKmZN8Az4BZw+QTv/wJIclRVh0kytB8Bc0OV/DlwI8n3qloDzv/XfxE4AK62Wb4k9cfKtiT14yWwnmTSaLxZsP5ZVReBB7MHVXUfmAeWgK2qutRoTknqimFbkjqR5EeSzYbjHQAvgC/AB2AfoKquAE+Bx0m+AdvARqt5Jakn9e+roSRJkqSWrGxLkiRJI/EHSUnqWFU9BFaPNX9K8uQ01iNJZ43HSCRJkqSReIxEkiRJGolhW5IkSRqJYVuSJEkaiWFbkiRJGolhW5IkSRrJX7wAIDjWtr/+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x110571cf8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "annual_mos = [calc_moment_release_from_fault_params(fault_area, fault_slip_rate,\n",
    "                                                    _M_min, M_upp, fault_b_gr,\n",
    "                                                    bin_width=_bin_size)\n",
    "              for M_upp in muppets]\n",
    "\n",
    "annual_mos_roundup = [calc_moment_release_from_fault_params_roundup(fault_area, \n",
    "                                                    fault_slip_rate,\n",
    "                                                    _M_min, M_upp, fault_b_gr,\n",
    "                                                    bin_width=_bin_size)\n",
    "              for M_upp in muppets]\n",
    "\n",
    "\n",
    "annual_mo_accumulation_rate = annual_mo_from_fault_params(fault_area, fault_slip_rate)\n",
    "\n",
    "\n",
    "plt.figure(figsize=(12,9))\n",
    "\n",
    "plt.plot(muppets, annual_mos, label='OQ-MBT moment release rate')\n",
    "plt.plot(muppets, annual_mos_roundup, label='Corrected moment release rate')\n",
    "\n",
    "\n",
    "plt.axhline(annual_mo_accumulation_rate, color='r', lw=0.5,\n",
    "            label='Moment accumulation rate')\n",
    "\n",
    "plt.gca().set_yscale('log')\n",
    "\n",
    "plt.xlabel('M_max')\n",
    "plt.ylabel('Annual moment release rate (N m / yr)')\n",
    "plt.title('Moment release rate as a function of M_max\\n(M_min={0}, bin_size={1})'.format(\n",
    "                                                                        _M_min, _bin_size))\n",
    "plt.legend(loc='upper right')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see from this graph that the corrected moment release rate equals the moment accumulation rate regardless of the initial `M_max`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "\n",
    "Rounding up `M_max` so that it falls on the edge of a bin solves these problems.  It is easy to implement (as is evident here)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}