Skip to content

Instantly share code, notes, and snippets.

@jdossgollin
Created June 14, 2019 18:14
Show Gist options
  • Save jdossgollin/c95f1e9f4a52d54746037477cf37249c to your computer and use it in GitHub Desktop.
Save jdossgollin/c95f1e9f4a52d54746037477cf37249c to your computer and use it in GitHub Desktop.
Do we over-use the GEV Distribution?
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In statistical analysis of hydroclimate data we're often interested in fitting distributions to observed variables such as rainfall and streamflow.\n",
"The General Extreme Value (GEV) distribution is commonly used because of its \"fat\" tails, but maybe that's not always necessary.\n",
"\n",
"In the following jupyter notebook I do a simple experiment.\n",
"First, I create several time series with random period and phases, which represent the indices of large-scale climate variables (such as ENSO, the PDO, etc).\n",
"Second, I draw fake streamflow data using a conditional log-normal model:\n",
"\n",
"$$\n",
"\\begin{align*}\n",
" \\mu &= \\mu_0 + X^T \\beta \\\\\n",
" \\sigma &= 0.1 \\mu \\\\\n",
" \\log Q &\\sim \\mathcal{N}(\\mu, \\sigma)\n",
"\\end{align*}\n",
"$$\n",
"\n",
"where $\\beta$ are the (also randomly chosen) dependences on the large-scale climate predictors and $X$ is the notation for these predictors.\n",
"\n",
"This is an incredibly over-simplified model and isn't meant to represent all the physics of streamflow -- just the fact that theres is variability present on many timescales (and here there's no trend!)\n",
"Interestingly, if you look at the resulting histogram of a long simulation from this model, you would think that a GEV fit was necessary to capture the \"fat\" tail of this model.\n",
"However, the real model doesn't have any GEV dependence -- the \"fat tail\" comes purely from the dependence on the climate predictors.\n",
"It's therefore worth asking ourselves whether we always need to use a GEV distribution (which comes with many problems including difficulty estimating the parameters) or whether we can instead use a simpler statistical distribution (such as the log-normal) plus the use of conditional estimation.\n",
"I won't get into the theory in this post, but it's certainly an interesting question!\n",
"\n",
"Here's the analysis:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from scipy.stats import genextreme as gev\n",
"from scipy.stats import norm\n",
"from scipy.integrate import trapz\n",
"import matplotlib.pyplot as plt\n",
"np.random.seed(523)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The question to answer here is whether we can trick ourselves into thinking that data is GEV, when in fact it's log-normal but conditional on several slowly varying predictors.\n",
"Here's a function that will give us sinusoidal functions using the model defined above."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def get_sinusoidal_signal(n, amplitude, period, phase, mu0=4.0, coeff_var=0.1):\n",
" time = np.linspace(0, n, n)\n",
" n_predictors = len(amplitude)\n",
" assert len(period) == n_predictors\n",
" assert len(phase) == n_predictors\n",
" predictors = np.array([np.sin(2 * np.pi * time / period[i] + phase[i]) for i in range(n_predictors)]).T\n",
" mu = np.repeat(mu0, time.size)\n",
" for i in range(n_predictors):\n",
" mu += amplitude[i] * predictors[:, i]\n",
" sigma = coeff_var * mu\n",
" sigma[np.where(sigma < 0.1)] = 0.1\n",
" observed = np.exp(np.random.normal(loc=mu, scale=sigma))\n",
" return time,observed"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's implement this with three signals"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"n_predictors = 3\n",
"time,observed = get_sinusoidal_signal(\n",
" n=100000,\n",
" period=[9, 13, 24],\n",
" amplitude=np.random.normal(loc=0, scale=0.5, size=n_predictors),\n",
" phase=np.random.uniform(low=0, high=2*np.pi, size=n_predictors),\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we will estimate a GEV fit, and a log-normal fit to the data.\n",
"(Note that because of how we fit the log-normal, we have to re-scale it)."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/usr/local/miniconda3/envs/jupyter-blog/lib/python3.7/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log\n",
" after removing the cwd from sys.path.\n"
]
}
],
"source": [
"muhat, sigmahat = norm.fit(np.log(observed))\n",
"shapehat, lochat, scalehat = gev.fit(observed)\n",
"sflo = np.linspace(0, 600, 1000) # we will make predictions here\n",
"norm_prob = norm.pdf(np.log(sflo), loc=muhat, scale=sigmahat)\n",
"norm_prob /= trapz(x=sflo, y=norm_prob)\n",
"gev_prob = gev.pdf(sflo, shapehat, scale=scalehat, loc=lochat)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can plot this fit, as well as the first 1000 years of the time series"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEyCAYAAADk0er5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd8FGX+wPHPNwVC7yJFARUsoICA4ClSLGAFC4KnIqd33NnO8/QO9PTkzrPzs3u2wwMVBcSCFeWAgEiR3pFeQiehJYGQ8vz+mNlkkszuzia72d3k+3698srulGef2Zmd+c7TRowxKKWUUkrFk4RoZ0AppZRSKlQawCillFIq7mgAo5RSSqm4owGMUkoppeKOBjBKKaWUijsawCillFIq7mgAEwNE5FER+U+08xEqEekpIr9U8GdeJCIbRCRTRAaWMY1RIvJhuPMWSSIyTETmRDsf4SIid4vIXns/NhIRIyJnRDtfquKIyGoR6R3tfMSTePmdiMip9m87MZKfowFMBbB3pO+vQESOOd7faox52hjz2wrKy10isk5EjtoXkG9EpE5Z0jLG/GiMOTPceQzin8Drxpjaxpgv3BawL/YrRSRbRPaIyJsiUr+C81mhAu1XERkrIv+Kdh59RCQZeBG4wt6P6dHOkwovEdkqIpeVmFYsCDfGtDfGpAZJp7V90U6KUFYrFRG5XERm2ueBdBFZJiIjRCTFnj9KRHJLXJMO2fPWicidLmk+ICKL/Hze1hLXs0wRaW6M2W7/tvPt5VJFJOzXOA1gKoC9I2sbY2oD24FrHdPGV1Q+RKQX8DRwizGmDnA2MKmMaUXrhNIKWO1vpog8BDwH/AWoB/Sw15kmItUqJIdU7PdT3v0ahX3ZFEghwH5UqiJUpsBIRAYBk4GPgFbGmEbAYKAlcIpj0YnOa5IxxndzNw4Y6pL07fY8f64tkd6u8m+NR8YY/avAP2ArcFmJaaOAD+3XrQED/AbYARwE/gB0A1YAh7BKIJzr3wmstZf9Huvgdfvsh4EvAuStOjAaK8jaC7wF1LDn9QbSgBHAHuAD3zTH+s2BT4H9wBbgj455FwCLgCN22i8GyMfvgI1ABvAl0NyevgkoAI4BmUD1EuvVtaffXGJ6bWAfcKfj+54MTASOAkuAjo7lRwA77Xm/AJfa0xOAkXY+0rGChIYl9ttd9vc3G5gK3FciL8uBG+zXZwHT7O38xZlvoJG97UeAn4EngTmh7ldgOJALnLC/m68cx+EI+5jKAZI87L95WMffbuB1oJpjvgHuATbY39uTwOn2Okfs76oa0A7IspfPBGY41j/Dfl0PeN/OxzbgMSDBnrcN6GK/vs1e7xz7/W/9fQ/6V7F/uJ/nhjmPYecy+Dk/2L8l37GSCVxo/w4fs4+FffaxUs+R7lB7XjrweInPGYX12//Q/qzfhuvY9vM9nA7MsPNyABgP1C/xHTxs/w4PY52TUhzz/2LnaRfWeb7wd1LicwTrevFQkP0yCvta4zKvJZCH4/qBdTN0AmjsdT/b01vbeU0CngLygeP2Pnw9UB5DOs6ifaBXtT8/P+zCg8qx49/Cuku9wt7xXwAnAS3sH20ve/mBWBf7s+2D5TFgrp/P7ol18f8HcBGlA4CXsS6aDYE6wFfAM/a83vbB/RxWoFMDRwCDdVJZDPwd60J1GrAZ6GfPnwfcbr+uDfTwk8e+9g/9fPtzXgNmB/r+HPP623lMcpk3DvjY8X3nAjcByVgnkC326zPtE0Fzx/443X79J2C+/UOvDrztSNO3394Hatnfz1DgJ0cezsE6SVa3l9mBFagm2dt7AGhvLzsB68RYC+iAFVD5C2CC7dexwL9cjsNlWHdmNTzsvy5YpVlJ9rauBf7kSM9gHTt1gfZYQdF0O516wBrgjhLfVVKJ9X0BzPvAFKxjsDWwHrjLMe8h+/U7WMHk3Y55D0b7N65/ZQpgXM8Pfo6VO7HOeafZy34GfGDPOwfrInmxfRyPxvqtOwOYXKzzZoJ97Ift2Hb5Hs4ALsf6zTfBurF5ucR38DPWzUND+7P/YM/rjxXMdcA6D3yE/wDmLHte6yD7ZRR+Ahh7/jTgMcf7Zwh801tqP7vtNyAV+G3Yj7NoH+hV7c/PD7vwoHLs+BaO+enAYMf7T30/MOA77JO7/T4ByMZ/KcyVWIHJIfuH/iKQiBXBZ2FfrO1lLwS22K97Y0XizruD3hQFMN2B7SU+6xHgv/br2VgXWNdI3rHOGOB5x/vaWCec1v6+P8eytwF7/Mx7Fpjm+L7nl/jOdmMFAmdgBYiXAckl0liLXRpjv29m58134jPAaY75dezvtJX9/ingPfv1YODHEum/DTxh749c4CzHvKfxE8AE2q/2vLG4BzB3Ot4H3H8un/cn4HPHewNc5Hi/GBjheP9/2CduAgQw9rbnYJeq2PN+D6Tar+8CvnTsj98CE+z324DzK+q3rH/+/+zjK9M+Hn1/2fgPYFzPD36OlenAPY73Zzp+h3/Hvqmw59XEOm85A5jZQfJe5mPbw/cyEFha4ju4zfH+eeAt+/V7wLOOee3wH8BcbM9znp8nOL53X3A4yv4+nPtlpmOd24Bf7NcJWCVg14ewn79w229EKIDRNjCxa6/j9TGX97Xt162AV0TkkN0YKwMrGGnhlqgx5jtjzLVY0f4ArLui32LdHdQEFjvSmmpP99lvjDnuJ7+tgOa+de31H8Vq7wDWhacdsE5EForINX7SaY51IfLlNxMrgHPdnhIOAI391Gs3s+f77HB8RgFW9VhzY8xGrBPYKGCfiEwQkeaObfzcsX1rsYpGm/pJ9yjwDTDEnjQEqwjZl1b3Et/XrcDJWN95kjMtHN+JmwD7NRBn+gH3n4i0E5Gv7UbRR7ACqsYl0vN6zAbSGOvO2bm92yja/7OAniJyMlawMxG4SERaY90NL/PwGapiDDTG1Pf9YVXD+OP1/AAlzhH26ySsY7U5xX+D2VjnDyfncR/RY1tETrLPITvttD90SXuP43W2I61i20Lgc4BvG5v5Jhhjhtjf+xKs34rPJOd+Mcb0ccz7DGgmIj2wblBrYp3DAnHu5zL1DC0rDWDi3w7g9yUOyBrGmLmBVjLGFBhjpmPVz3bAurgfw6rC8KVTz1gNjwtXC5KPLSXyUccYc5X9eRuMMbdgVYM9B0wWkVou6ezCupgCYC/TCKsKJZh5WHfvNzgn2mlciXXn5nOKY34CVrXQLjuvHxljLrbzYez8+rbxyhLbmGKMceat5Hf0MXCLiFyIVVw905HWrBJp1TbG3I3V9iOP4g3vTvWw/W771S1PbnkNuP+AN4F1QFtjTF2s4Ea85ClEB7Dupls5pp2Kvf/tADMb+CPWnfRRrAvAcKy7+4II5ElFWIDzg9uxW+wcgXV85GEFFbuxfssAiEgNrPNHsY8r8T6Sx/Yz9uedZ6d9Wwhp78b7OWAd1m/khgDLBGUHfJOxqr9vxyrdPFGeNH1JhyGNUjSAiX9vAY+ISHsAEalnt0YvRUQGiMgQEWkglguAXljVKQXAu8BLInKSvXwLEennMR8/A0fsLns1RCRRRDqISDc7rdtEpIn9OYfsdfJd0vkI+I2IdBKR6lh3QwuMMVuDZcAYcxirGPo1EekvIsn2nfknWCUsHzgW7yIiN9ilNX/CCnzmi8iZItLX/uzjWEGdL59vAU+JSCt7m5qIyIAg2foW62T7T6zW/74L7NdAOxG53c5nsoh0E5GzjdX18DNglIjUFJFzgDv8fUCg/Wovshervj6QgPsPqzrsCJApImcBdwdJr0zsbZ+E9T3Xsb/rP2PdufrMAu6z/4NVPO18r+JMgPPDfqyG+87j92PgQRFpIyK1sc4RE40xeVgX32tF5Fd2r8N/EDxgiOSxXQe7ikVEWmA1yvVqEjBMRM4RkZpY1cuujFVP8xDwhIj8znEuaEvxEmIvxmFVcd9I4N5HofByDgqZBjBxzhjzOdYdywS7iHIVVmmDm4NYPXw2YP1gPwReMEVduUdgNY6bb6f1P6z6ZS/5yAeuBTphNYg9APwHq1gfrAZpq0UkE3gFGOJWHWWXHjyO1c5nN1Yr/iEllwuQj+ex7qBG29u4AKt04VJjTI5j0SlYP9KDWHcaNxhjcrEa2z1r538P1h3ho/Y6r2A15vtBRI5iBQjdg+QnBysYuQwrOPNNP4rVQHsI1h3lHooaSIN1Qa5tTx8L/DfAxwTbr2OAc+yqIdexczzsv4eBX2P1wngXq+omUu7Haju0GZiD9b2955g/C+vCMNvPexV/XM8PdonAU8BP9vHbA+tY+ABrf2/ButG4H8AYs9p+PQHr/HEUq01bDv5F8tj+B1YD/cNYVTGfeV3RGPMdVseKGVjn5RlBlp8I3IxVyrMD6zc8Caux+yeORQdL8XFbMn03rbbZdn53GmMWes1vEK8AN4nIQRF5NUxpInYDG6WUUqpSsUtoDmFVD22Jdn5UeGkJjFJKqUpDRK61q15rYZXErsTqLaMqGQ1glFJKVSYDsKpldwFtsaqjtKqhEtIqJKWUUkrFHS2BUUoppVTciVgAY3dHXeb4OyIifxKRhiIyTUQ22P8b2MuLiLwqIhtFZIWInB+pvCmllFIqvlVIFZKIJGINstMduBfIMMY8KyIjgQbGmBEichVW97er7OVeMcYE7KLauHFj07p165Dzk5WVRa1abmOoxb/KvG1QubdPty1+lWX7Fi9efMAY0yT4krGlfv365owzzoh2NjyJp+NO82oJ6XcRjucRBPvDGu/iJ/v1L0Az+3Uzip678DZwi2OdwuX8/XXp0sWUxcyZM8u0XjyozNtmTOXePt22+FWW7QMWmQo4/4b7r127diFva7TE03GnebWE8ruoqDYwQ7BGTwRoaozZbQdPu7EGCgPrWSfO5z6k4e35N0oppZSqYiJehWQP57wL6xk7e0XkkLEeMOWbf9AY00BEvgGeMcbMsadPB/5qjFlcIr3hWM89oWnTpl0mTJgQcp4yMzOpXdvLc+XiT2XeNqjc26fbFr/Ksn19+vRZbIzpGqEshZXzvNukSZMukyZNinKOvImn407zagnpd+G1qKasf1h98n9wvNcqpAiqzNtmTOXePt22+KVVSLEpno47zasllN9FUnhjJ1e3UFR9BNazZO7Aet7MHVjPpPFNv09EJmA14j1s7KompZRSlVtubi5paWkcP17qEWllVq9ePdauXRu29CKpquU1JSWFli1bkpycXOY0IhrA2E/QvBz4vWPys8AkEbkL2A74npz8LVYPpI1ANvCbSOZNKaVU7EhLS6NOnTq0bt0akWAPkPbm6NGj1KlTJyxpRVpVyqsxhvT0dNLS0mjTpk2Z04loAGOsJ4k2KjEtHbjUZVmD1cVaKaVUFXP8+PGwBi8qdokIjRo1Yv/+/eVKR0fiVUopFRM0eKk6wrGvNYBRSimlsC6qDz30UOH70aNHM2rUqIDrfPnllzz77LPl/uzx48dz3333lTudklq3bs2NN95Y+H7y5MkMGzYs7J8TyKhRoxg9enTY09UARimllAKqV6/OZ599xoEDBzyvc9111zFy5MgI5qr8Fi1axOrVq8u0bl5eXphzEz4awCillFJAUlISw4cP56WXXio176uvvqJ79+507tyZyy67jL179wIwduxY7rvvPg4fPkzr1q0pKCgAIDs7m1NOOYXc3Fw2bdpE//796dKlCz179mTdunWe8/Txxx9z7rnn0qFDB0aMGFE4fcyYMbRr147evXvzu9/9LmDpzcMPP8zTTz9danpGRgYDBw7kvPPOo0ePHqxYsQKwSkyGDx/OFVdcwdChQxk7diwDBw7k2muvpU2bNrz99tu8+OKLdO7cmR49epCRkQHAu+++S7du3ejYsSM33ngj2dnZnrezLDSAUUoppWz33nsv48eP5/Dhw8WmX3zxxcyfP5+lS5cyZMgQnn/++WLz69WrR8eOHZk1axZgBTz9+vUjOTmZ4cOH89prr7F48WJGjx7NPffc4ykvu3btYsSIEcyYMYNly5axcOFCvvjiC3bt2sWTTz7J/PnzmTZtWtCA6Oabb2bJkiVs3Lix2PQnnniCzp07s2LFCp5++mmGDh1aOG/x4sVMmTKFjz76CIBVq1bx0Ucf8fPPP/Pkk09Ss2ZNli5dyoUXXsj7778PwA033MDChQtZvnw5Z599NmPGjPG0nWVVEePAKKWUUp7946vVrNl1pNzp5Ofnk5iYCMA5zevyxLXtg65Tt25dhg4dyquvvkqNGjUKp6elpTF48GB2797NiRMnXLv/Dh48mIkTJ9KnTx8mTJjAPffcQ2ZmJnPnzmXQoEGFy+Xk5HjK/8KFC+nduzdNmljPNrz11luZPXs2AL169aJhw4YADBo0iPXr1/tNJzExkb/85S8888wzXHnllYXT58yZw6effgpA3759SU9PLwzcrrvuumLb36dPH+rUqUOdOnWoW7cu1157LQDnnntuYcnNqlWreOyxxzh06BCZmZn069fP03aWlZbAKKWUUg5/+tOfGDNmDFlZWYXT7r//fu677z5WrlzJ22+/7Trg3nXXXcd3331HRkYGixcvpm/fvhQUFFC/fn2WLVtW+Od1EDjj51E//qbn5+fTqVMnOnXqxN///vdi826//XZmz57N9u3bA6bj6x1U8mnT1atXL3ydkJBQ+D4hIaGwncywYcN4/fXXWblyJU888URYByV0oyUwSimlYoqXkhIvyjrgWsOGDbn55psZM2YMd955JwCHDx+mRQvr+cLjxo1zXa927dpccMEFPPDAA1xzzTUkJiZSt25d2rRpwyeffMKgQYMwxrBixQo6duwYNB/du3fngQce4MCBAzRo0ICPP/6Y+++/n65du/Lggw9y8OBB6tSpw6effsq5555LYmIiy5Ytc00rOTmZBx98kGeffZa+ffsCcMkllzB+/Hgef/xxUlNTady4MXXr1g35+/I5evQozZo1Izc3l/Hjxxd+X5GiJTBKKaVUCQ899FCx3kijRo1i0KBB9OzZk8aNG/tdb/DgwXz44YcMHjy4cNr48eMZM2YMHTt2pH379kyZMsV13bFjx9KyZcvCv/z8fJ555hn69OlDx44dOf/88xkwYAAtWrTg0UcfpXv37lx22WWcc8451KtXL+g23XXXXcV6FY0aNYpFixZx3nnnMXLkSL+BmVdPPvkk3bt35/LLL+ess84qV1qeeH1oUiz+6cMcS6vM22ZM5d4+3bb4pQ9zLL81a9aEPc0jR46EPc1ICTWvR48eNcYYk5uba6655hrz2WefRSJbrsL1vbrt81B+F1oCo5RSSsWZUaNG0alTJzp06ECbNm0YOHBgtLNU4bQNjFJKKRVnIjGybbzREhillFJKxR0NYJRSSikVdzSAUUoppVTc0QBGKaWUUnFHAxillFIK2Lt3L7/+9a857bTT6NKlCxdeeCGff/45qamp1KtXr3CU206dOvG///2P3r178/333xdL4+WXX/b8rCNVPtoLSSmlVJVnjGHgwIHccccdhQ8w3LZtG19++SUNGjSgZ8+efP3118XW2bRpExMmTCj2zJ8JEybwwgsvVGjeqyotgVFKKVXlzZgxg2rVqvGHP/yhcFqrVq24//77/a5z00038fXXXxc+nHHr1q3s2rWLiy++OOL5VRrAKKWUUqxevZrzzz/f7/wff/yxWBXSpk2baNSoERdccAFTp04FrNKXwYMHFz4QUUWWViEppZSKLd+NhD0ry51Mjfw8SLQvcyefC1c+63nde++9lzlz5lCtWjVeeOEF1yokgFtuuYUJEyYwYMAAJkyYwHvvvVfufCtvtARGKaVUlde+fXuWLFlS+P6NN95g+vTp7N+/P+B6AwcOZPr06SxZsoRjx44FLMVR4aUlMEoppWJLCCUlgRw7epQ6dep4WrZv3748+uijvPnmm9x9990AZGdnB12vdu3a9O7dmzvvvJNbbrmlXPlVodESGKWUUlWeiPDFF18wa9Ys2rRpwwUXXMAdd9zBc889B5RuAzN58uTCdW+55RaWL1/OkCFDopX9KimiJTAiUh/4D9ABMMCdwC/ARKA1sBW42RhzUKxWT68AVwHZwDBjzBKXZJVSSqmwa9asGRMmTHCdd/jwYb/rXX/99RhjIpUt5UekS2BeAaYaY84COgJrgZHAdGNMW2C6/R7gSqCt/TcceDPCeVNKKaVUnIpYACMidYFLgDEAxpgTxphDwABgnL3YOGCg/XoA8L6xzAfqi0izSOVPKaWUUvErkiUwpwH7gf+KyFIR+Y+I1AKaGmN2A9j/T7KXbwHscKyfZk9TSikVIwqO+a9KUaoiRbINTBJwPnC/MWaBiLxCUXWRG7eRf0pVKorIcKwqJpo2bUpqamrIGcvMzCzTevGgMm8bVO7t022LX5V9+5zn3XObVY/IttarV4+jR4+GNc38/PywpxkpVTGvx48fL9exFMkAJg1IM8YssN9Pxgpg9opIM2PMbruKaJ9j+VMc67cEdpVM1BjzDvAOQNeuXU3v3r1DzlhqaiplWS8eVOZtg8q9fbpt8auyb5/zvHte85QynXeDWbt2recuz14dDaEbdbRVxbympKTQuXPnMq8fsSokY8weYIeInGlPuhRYA3wJ3GFPuwOYYr/+Ehgqlh7AYV9Vk1JKKaWUU6R7Id0PjBeRFUAn4GngWeByEdkAXG6/B/gW2AxsBN4F9HnkSikVYyrzU37S0tIYMGAAbdu25fTTT+eBBx7gxIkTjB07lvvuuy/a2Suldu3a0c5CVEV0HBhjzDKgq8usS12WNcC9kcyPUkqp+PDvZf8udxo5OTlUr14dgHs6Bb4nNsZwww03cPfddzNlyhTy8/MZPnw4f/vb32jfvn2581JSXl4eSUk6GH556Ei8SimlqrwZM2aQkpLCb37zGwASExN56aWXeO+998jOzmbHjh3079+fM888k3/84x8AZGVlcfXVV9OxY0c6dOjAxIkTAVi8eDG9evWiS5cu9OvXj927rdYQvXv35tFHH6VXr1489dRTtG7dmoKCAsB6bMEpp5xCbm4umzZton///nTp0oWePXuybt06ALZs2cKFF15It27dePzxxyv6K4o5Gv4ppZSq8lavXk2XLl2KTatbty6nnnoqeXl5/Pzzz6xatYqaNWvSrVs3rr76arZt20bz5s355ptvAGu03tzcXO6//36mTJlCkyZNmDhxIn/7298Kn1J96NAhZs2aBcCSJUuYNWsWffr04bvvvqNfv34kJyczfPhw3nrrLdq2bcuCBQu45557mDFjBg888AB33303Q4cO5Y033qjYLygGaQCjlFKqyjPGYD3Rxn365ZdfTqNGjQC44YYbmDNnDldddRUPP/wwI0aM4JprrqFnz56sWrWKVatWcfnllwNWl+NmzYrGZB08eHCx1xMnTqRPnz58+umn/PGPfyQzM5O5c+cyaNCgwuVycnIA+Omnn/j0008BuP322xkxYkT4v4g4ogGMUkqpKq99+/aFwYHPkSNH2LFjB4mJiaWCGxGhXbt2LF68mG+//ZZHHnmEK664guuvv5727dszb94818+pVatW4evrrruORx55hIyMDJYtW0bfvn3Jysqifv36LFu2zHV9tyCrqtI2MEoppaq8Sy+9lOzsbN5//33AKjl56KGHGDZsGDVr1mTatGlkZGRw7NgxvvjiCy666CJ27dpFzZo1ue2223j44YdZsmQJZ555Jvv37y8MYHJzc1m9erXrZ9auXZsLLriABx54gH79+pGYmEjdunVp06YNn3zyCWCVAC1fvhyAiy66qPBhk+PHj4/0VxLzNIBRSikVgsr51GUR4fPPP+eTTz6hbdu2tGvXjpSUFJ5++mkALr74Ym6//XY6derEjTfeSNeuXVm5ciUXXHABnTp14qmnnuKxxx6jWrVqTJ48mREjRtCxY0c6derE3Llz/X7u4MGD+fDDD7nxxhsLp40fP54xY8bQsWNH2rdvz5Qp1nBpr7zyCm+88QbdunUL+HTsqkKrkJRSSsWcYN2evQh1xNhTTjmFr776qtT0YcOGMWzYsFLT+/XrR79+/UpN79SpE7Nnzy413W3Y/JtuugljTLGh+du0acPUqVNLLdumTZtiVVMjRwZ6Ok/lpyUwSimllIo7GsAopZRSKu5oAKOUUkqpuKMBjFJKqZhgPVFGVQXh2NcawCillIq6lJQU0tPTNYipAowxpKenk5KSUq50tBeSUkqpqGvZsiVpaWns378/bGkeP3683BfJilLV8pqSkkLLli3LlYYGMEoppaIuOTmZNm3ahDXN1NRUOnfuHNY0I0XzGjqtQlJKKaVU3NEARimllFJxRwMYpZRSSsUdDWCUUkopFXc0gFFKKaVU3NEARimllFJxRwMYpZRSSsUdDWCUUkopFXc0gFFKKeWZRDsDStk0gFFKKaVU3IloACMiW0VkpYgsE5FF9rSGIjJNRDbY/xvY00VEXhWRjSKyQkTOj2TelFJKKRW/KqIEpo8xppMxpqv9fiQw3RjTFphuvwe4Emhr/w0H3qyAvCmllFIqDkWjCmkAMM5+PQ4Y6Jj+vrHMB+qLSLMo5E8ppZRSMS7SAYwBfhCRxSIy3J7W1BizG8D+f5I9vQWww7Fumj1NKaWUUqqYpGALiEhfYL4xJrsM6V9kjNklIicB00RkXaCPcplmXPIzHKuKiaZNm5KamhpypjIzM8u0XjyozNsGlXv7dNviV2XfPud5t2Oz5LjZ1njaL5rXMjDGBPwD3gfWA/OA54FrgQbB1nNJZxTwMPAL0Mye1gz4xX79NnCLY/nC5fz9denSxZTFzJkzy7RePKjM22ZM5d4+3bb4VZbtAxaZEM+jsfDXsVm1kLc1WuLpuNO8WkL5XQStQjLGDDXGtANuxKrWeQPYH2w9EaklInV8r4ErgFXAl8Ad9mJ3AFPs118CQ+3eSD2Aw8aualJKKaWUcvJShXQb0BM4FzgAvA786CHtpsDnIuL7nI+MMVNFZCEwSUTuArYDg+zlvwWuAjYC2cBvQtsUpZRSSlUVQQMY4GVgE/AWMNMYs9VLwsaYzUBHl+npwKUu0w1wr5e0lVJKKVW1ealCagzcCaQAT4nIzyLyQcRzppQk9fVFAAAgAElEQVRSKgaV6luhVFQEDWBEpC5wKtAKaA3UAwoimy2llFJKKf+8VCHNcfy9boxJi2yWlFJKKaUCCxrAGGPOA7B7FGnZoVJKKaWizksVUgcRWYrVBXqNPapuh8hnTSmllFLKnZdHCbwD/NkY08oYcyrwkD1NKaWUUioqvAQwtYwxM31vjDGpQK2I5UgppZRSKggvjXg3i8jjgK/r9G3AlshlSSmllFIqMC8lMHcCTYDPgM/t1zpKrlJKKaWixksvpIPAHysgL0oppWKcRDsDStn8BjAi8hUBuk0bY66LSI6UUkoppYIIVAIzusJyoZRSSikVgkABzN+NMZeKyHPGmBEVliOllFIxTMczVbEhUADTTER6AdeJyARKVH0aY5ZENGdKKaVijraBUbEiYAkMMBJoCbxYYp4B+kYqU0oppZRSgfgNYIwxk4HJIvK4MebJCsyTUkoppVRAXrpRPyki5wGtncsbYz6LYL6UUkoppfwKGsCIyHvAecBqoMCebLAGtlNKKaWUqnBeHiXQwxhzTsRzopRSSinlkZdHCcwTEQ1glFJKIdqNWsUILyUw47CCmD1ADlYvOmOMOS+iOVNKKaWU8sNLAPMecDuwkqI2MEoppZRSUeMlgNlujPky4jlRSimllPLISwCzTkQ+Ar7CqkICtBu1UkpVVcYYRHRMXhVdXhrx1sAKXK4ArrX/rvH6ASKSKCJLReRr+30bEVkgIhtEZKKIVLOnV7ffb7Tntw51Y5RSSkVegbbjVTHAy0B2vynnZzwArAXq2u+fA14yxkwQkbeAu4A37f8HjTFniMgQe7nB5fxspZRSYWT34kCfiqSiLWgJjIikiMi9IvJvEXnP9+clcRFpCVwN/Md+L1jPUJpsLzIOGGi/HmC/x55/qWgZpVJKxRijJTAqJnipQvoAOBnoB8zCerjjUY/pvwz8laLeS42AQ8aYPPt9GtDCft0C2AFgzz9sL6+UUiqGFBiNYFT0eWnEe4YxZpCIDDDGjLMb9H4fbCURuQbYZ4xZLCK9fZNdFjUe5jnTHQ4MB2jatCmpqakeNqG4zMzMMq0XDyrztkHl3j7dtvhV2bfPed7t0iyB2bNmUS3Jy/1vdMXTftG8loExJuAf8LP9fzbQAWgMbPaw3jNYJSxbgT1ANjAeOAAk2ctcCHxvv/4euNB+nWQvJ4E+o0uXLqYsZs6cWab14kFl3jZjKvf26bbFr7JsH7DIBDmPxuJfl2YJJut4TsjbGw3xdNxpXi2h/C68hNDviEgD4HHgS2AN8LyHwOgRY0xLY0xrYAgwwxhzKzATuMle7A5giv36S/s99vwZ9sYopZSKIQXaCEbFAC+9kP5jv5wFnBaGzxwBTBCRfwFLgTH29DHAByKyEcjACnqUUkrFGGN0UHYVfUEDGBFpCjwNNDfGXGk/2PFCY8yYIKsWMsakAqn2683ABS7LHAcGeU1TKaVUdGgJjIoFXqqQxmK1T2luv18P/ClSGVJKKRXbCgq0BEZFn5cAprExZhJ2V2hjdXHOj2iulFJKxay8fA1gVPR5CWCyRKQRdpdmEemBNUaLUkqpKihfS2BUDPAyDsyfsXoInS4iPwFNKOpFpJRSqorJK9BCeBV9AQMYEUkAUoBewJlYg839YozJrYC8KaWUikH5eVoCo6IvYABjjCkQkf8zxlwIrK6gPCmllIphedoLScUAL21gfhCRG/XBikoppQDy87UKSUWf1zYwtYA8ETlO4dPUTd2I5kwppVRMyi/IC76QUhHmZSTeOhWREaWUUvEhP09LYFT0Ba1CEpHpXqYppZSqGvLztQRGRZ/fEhgRSQFqAo3thzn62sDUpWhUXqWUUlVMgbaBUTEgUBXS77EeGdAcWExRAHMEeCPC+VJKKRWj8vN1JA0VfX4DGGPMK8ArInK/Mea1CsyTUkqpGKZtYFQs8NsGRkS6icjJvuBFRIaKyBQReVVEGlZcFpVSSsWSAm0Do2JAoEa8bwMnAETkEuBZ4H2s5yC9E/msKaWUikV5GsCoGBCoDUyiMSbDfj0YeMcY8ynwqYgsi3zWlFJKxaK8PG0Do6IvUAlMooj4ApxLgRmOeV4GwFNKKVUJ5eZqCYyKvkCByMfALBE5ABwDfgQQkTOwqpGUUkpVQbm5WgKjoi9QL6Sn7AHrmgE/GGN8T+9KAO6viMwppZSKPVoCo2JBsKdRz3eZtj5y2VFKKRXr8vI0gFHR5+Vp1EoppVShXG3Eq2KABjBKKaVCoiUwKhZoAKOUUiokWgKjYkGghzkeBYy/+caYuhHJkVJKqZimJTAqFgTqhVQHQET+CewBPsB6oOOtQJ1gCdtPs54NVLc/Z7Ix5gkRaQNMABoCS4DbjTEnRKQ61ki/XYB0YLAxZmvZN00ppVQk5GsAo2KAlyqkfsaYfxtjjhpjjhhj3gRu9LBeDtDXGNMR6AT0F5EewHPAS8aYtsBB4C57+buAg8aYM4CX7OWUUkrFmNzcE9HOglKeAph8EblVRBJFJEFEbgWCPorUWDLtt8n2nwH6ApPt6eOAgfbrAfZ77PmXioh43A6llFIV5MSJnGhnQSlPAcyvgZuBvfbfIHtaUHbQswzYB0wDNgGHjDG+8sc0oIX9ugWwA8Cefxho5G0zlFJKVZT8E8ejnQWlgj/TyG6HMqAsiRtj8oFOIlIf+Bw4220x+79baUupRsQiMhwYDtC0aVNSU1NDzldmZmaZ1osHlXnboHJvn25b/Krs2+c873ZplkBW5tG42N542i+a19AFDWBEpB3wJtDUGNNBRM4DrjPG/MvrhxhjDolIKtADqC8iSXYpS0tgl71YGnAKkGY/RLIekOGS1jvAOwBdu3Y1vXv39pqNQqmpqZRlvXhQmbcNKvf26bbFr8q+fcXOu80TTVKCiYvtjaf9onkNnZcqpHeBR4BcAGPMCmBIsJVEpIld8oKI1AAuA9YCM4Gb7MXuAKbYr7+032PPn+F4/pJSSqkYkZ+rbWBU9AUtgQFqGmN+LtGe1ksfumbAOBFJxAqUJhljvhaRNcAEEfkXsBQYYy8/BvhARDZilbwEDZKUUkpVPMk/QV5+AUmJOhaqih4vAcwBETkduz2KiNwE7A62kl1S09ll+mbgApfpx7EaCCullIphyeSRmZNH/ZrVop0VVYV5CWDuxar7PEtEdgJbsAazU0opVQVVI4+D2bkawKio8hLAbDPGXCYitYAEY8zRSGdKKaVU7KomeaRn5tCmca1oZ0VVYV4qMLeIyDtYPYgygy2slFKqMhOqkcuBTB2NV0WXlwDmTOB/WFVJW0TkdRG5OLLZUkopFYuMCMnkkZ6lPZFUdAUNYIwxx4wxk4wxN2A1yq0LzIp4zpRSSsUgoTq5ZGgJjIoyT33gRKSXiPwb6+nRKViPFlBKKVXFGIS6SSdIz9IARkWXl5F4twDLgEnAX4wxWRHPlVJKhUHqL/vo0qoBdVKSo52VSsNIAg2Sctl7RJ+HpKIrYAmMPQjdf40x1xtjPtbgRSkVL3YdOsaw/y7kwYnLo52VSsUg1E/KZeehY9HOiqriAgYw9sMY+1RQXpRSKmyyT+QDsHm/dp4MKxHqJuSw86AGMCq6vIwDM1dEXgcmAoUlMMaYJRHLlVJKqZhkSKBWgtUGJvtEHjWrebmMKBV+Xo68X9n//+mYZoC+4c+OUkqFi/0sWAm8lAqNEaGGsUpfdh06xhkn1YlyjlRVFTSAMcZoFZJSSinAKoGpVmA14N2RoQGMip6g3ahFpKmIjBGR7+z354jIXZHPmlJKqVhjREjOzwZg4z5tX6Six8s4MGOB74Hm9vv1wJ8ilSGlKtLCrRncNXYh+QUm2llRYWa0BikyJAE5kcVJtZLYsE8fjaeix0sA09gYMwkoADDG5AH5Ec2VUhXk7g+XMH3dPtIzdVj0ykpEQ5hwKiABMHRsIqzfqyUwKnq8BDBZItIIu0WciPQADkc0V0pVGL1NVyoURqzLRocG+Wzcl4kxWnqposNLAPNn4EvgdBH5CXgfuD+iuVLKj/m78mg98hsOZYd3GHPRCKbS0ctqZBhJBODsBvlk5uSxNT07yjlSVZWXXkhLRKQX1lOpBfjFGJMb8Zwp5eK7rdahtyPjGPVrVit3enrzWHlpG5jI8JXAnF0/H0hi2Y6DtGlcK7qZUlWSl15Ig4AaxpjVwEBgooicH/GcxYiZv+zjw/nbop0NZfM1tk1OCs9lyRe/JOhVTilPCuwApkX149SqlsjS7YeinCNVVXmpQnrcGHNURC4G+gHjgDcjm63Y8Zv/LuSxL1ZFOxvKlldg/U9K8PQg9aAK7Nt0beiplDfGvmwkHDvIeS3rs2T7wSjnSFVVXq4Cvh5HVwNvGmOmAOUvu1eqDPLtIpOkMBeZaEPEysfgC06jnJFKxkgi+ZIImXvo1qYha3YdCXubNKW88BLA7BSRt4GbgW9FpLrH9ZQKO18JTLj44hYdBkYp77JS6sCR3fRq15gCA3M2Hoh2llQV5CUQuRlrILv+xphDQEPgLxHNVRV2z/jFXPHSrGhnI2bl2RFHuOINU5ieRjBKeZVVox4c3UXHlvWpk5LE7PX7o50lVQV56YWULSJbgStFpD/wkzHmh4jnrIr6duWeaGchpuXbJTAFYaryMaVeqMqiwD5WtIt8+GWl1IUju0hKTKBn28ak/rKf/AJDoraGVxXISy+kv2M13G0ENAb+KyKPeVjvFBGZKSJrRWS1iDxgT28oItNEZIP9v4E9XUTkVRHZKCIrqlJPp8osMyePMXO2UBCmOhpfG5iwtVnRKqRKK1xBriotq0ZdOLIbjOHKDs3YdzSHhVszop0tVcV4qUK6BehmjHnCGPME0AO41cN6ecBDxpiz7XXuFZFzgJHAdGNMW2C6/R7gSqCt/TecKtTTKZg/T1rGkHfmuc7LyDrBgNfnkHYwNgeTevrbtTz59Rpm/rIvLOn52sCEOX7RKiQFWIHxU9+sYaM+4yegozXqQ24WZGdw6dknUSM5ka+W74p2tlQV4yWA2QqkON5XBzYFW8kYs9sYs8R+fRRYC7QABmCV6GD/H2i/HgC8byzzgfoi0szLRlR2ny3ZyfzN7nc3U5btZHnaYd6dvbmCc+XN4Wxr4LnsE+F5fFZ+hEpM9Ga98inqIu99nd2Hj/Puj1sYOubnCOWqcjhUu7H1In0jNaslcfk5Tfl25W6O5+pj8lTF8RvAiMhrIvIqkAOsFpGxIvJfYBUQ0hO8RKQ10BlYADQ1xuwGK8gBTrIXawHscKyWZk9TAfjqnPNj9QpsXzzCXZwftjYwdjpa3VD5lCfI1SrFwIoCmA0ADO52Cgezc/l25e4o5kpVNYEa8S6y/y8GPndMTw3lA0SkNvAp8CdjzJEAA4a5zSh1GhGR4VhVTDRt2pTU1JCyA0BmZmbI65Xlc8rD7fPcpm3cbpVwpO3cRWpqepm2LZL27zsOwJo1a6l3aEPY0l24cCF76yaWO528fOuOcd68+TSpGd3RAWJt34VTNLZt4yFr32ZlZXn+7IzjVh1lTk5OSPmtzPsOip93mzRpQsqRdhRIIjuWzmDL4ZYYY2hWS3jt+5U0PLIxyrktEk/7RfMaOr8BjDFmHICIpABnYAUTm4wxx70mLiLJWMHLeGPMZ/bkvSLSzBiz264i8jWOSANOcazeEihVqWqMeQd4B6Br166md+/eXrNTKDU1Fc/rTf0GwPvy5eX2eQHysGvBdlizkpNPbkbv3ueFtm0V4LPdS2HPLs4++2x6dw5DgZr9XZzfpSsdWtQrd3IJ06dCfj49evTglIY1y51eecTavgunaGxbnW0HYf5catWqRe/el3haZ/fhY5A6g2rVq4WU38q876D4eff0tqebfc0zOFyrIa1q5tDK3u67U7by9ymrqdX6PLq1bhjF3BaJp/2ieQ1doCqkJBF5HiuwGAd8COwQkeftwCQgsYpaxgBrjTEvOmZ9Cdxhv74DmOKYPtTujdQDOOyralL+Jdp7MD9Gy7x9BW7hbiQbvka8WoVUWZkyPCbC1+VaD4fgDtRrBruXF76/qUtLGteuxkvT1kcxV6oqCVRm/gLWoHVtjDFdjDGdgdOB+sBoD2lfBNwO9BWRZfbfVcCzwOUisgG43H4P8C2wGdgIvAvcU5YNqmjGGDJz8qL2+QkS221gfJeOcGcvXAGRL18x+vWF3cKtGZz2yDekZ+ZEOysRV5aYXh874N3++i3g8A7IskbhrVktibt7n8HcTenM25Qe5dypqiBQAHMN8Du7BxEAxpgjwN3AVcESNsbMMcaIMeY8Y0wn++9bY0y6MeZSY0xb+3+GvbwxxtxrjDndGHOuMWZRsM+IBRMX7qDDE9+zaX9I7ZrDxhfAxMIFeP7mdMbM2VJsmu/uN9wFRGHvhRTe5GLW27M2U2Bg8TZvD+D7esUuJi3cEXzBGFSesYKqyvFQHvvqtwTgqzlPFk67tfupnFw3hae/XRuzpcKq8ggUwBjjcgYwxuSjv+9C/1u7F4DN+7Oi8vmFvZBi4GQx5J35PPn1mmLTfDe0J8L8EKNwj8Trlt6KtEO8Oj20hsfHc/N55ru1ZJ8oKpWbumo3v/lv2brlrt97lCnLdpZp3UC8fnv3fbSUv366IuyfXxF8P4lQClUiVWJYGe2v35wChJMzthdOS0lO5NGrz2blzsOMX7AtirlTVUGgAGaNiAwtOVFEbgPWRS5LFc8Yw4q0QwGXWbPriJ91rf/RKnlOiJNu1Cfywjs+RNjawPieheSS3nWv/8SL09aHdCf/4fxtvD1rM2+mFg2V9IcPlzDzl7I9K+aKl2bzwIRlZVo3Vu07erxCxgspSwlMjP6KYlJucgr767eg5f7iw4Jde14zLj6jMS9M/YW9Rzz3+VAqZIECmHuxRs9NFZH/E5HRIjIL+CNWNVKlMXbuVq57/SfmbPD/RNVHPnO/Cy3LYFmREK6h+iMlNz+8+QvXowSKkvGfXk4IpUe+7Qx3iVO4RfNwveCp6dw+ZkHEP6csR4g25g5N2kln0PTgDjhedIMnIjw5sAN5BYY/T1oW8+cmFb/8BjDGmJ3GmO7AP7FG490O/NMYc4ExJvxl2lG0drf14ws4HH+QCCVaAYzv5BALVUhufPkK/0B23pddsDmdXYeOBVzGLXvV7C5eoZQWFPW6im3Rzt/Crd7a4JRHWW4uNH4JTVqTM0gwBbB1TrHpbRrXYtR15/DTxnTejtFRwlX88/I06hnAjArIS9QUVgMFONH5m2UK50cngolUgBAuvmyFf+h/7wkOfmc+NZITWftk/9Lp2P+d+Vuy/SAfzNtGcqJwIh+O5eZT3+NnFbWhKJ0/Y0xIXXrDta5TpALtz5emkV9gdaWNFQUeftel1yl6OpYKbnejVhxPTiFlzRQ4q3jfjpu7nsLs9QcY/cMvtG9el0vaNYlSLlVlFd2hR2OElyDE30mwoGjlqIjZti8lRLMEBqwgxE1hGxjHBevOsQv5fOlOjtvVQMdCeI5ToItleb6CvBgtYfN5cOJyHv5kefAFK1BZjrmq1q2+vAoSktjSrAOs+wZyi7d3ERGeu+k82p5Um3vHL2HDXn1ApgovDWBwnKwCXHwS/FyZCgfLCnOevMqPdgQVhO9rC3c9eNjawBSmV3pe9STr5+Ev+AmYrkt65Qni8sLchihawrXfvH1YxX1UVbax5blw4iis/67UvNrVkxgzrBvVkxMZ9t+F7AxSlatUKDSAoeju21+QAsHDg/IU7+flF3CgjAOL+QKYaDci9seXrXCXFIUjNWOMo4qrdIq+AOZ4rvcGuYUjubrMK08Ml1sQ3kbB0SphqMi2WoVtYEII7gsKS+SUV2lNzuBIjfrsnPlP/r3s36Xmt6hfg/8O68aR47n8+t357DmsPZNUeGgAA4Vnq7LEAOHoRv3k12vo+q//lWlE31ht+1JS+AeyK3+CEx0DtLmXwFgPiyxTI94YLYEJZ5xrjGHxtoyQSlWCVYVl5eQFHdLAq7Icc0VVSPHxu4oFRhJYddqFtDiwmUaH3Z/+cm7Lerx/5wWkZ57glnfnsyMjQIcJpTzSAAZHG5gynN19pTflKQGZunoPAJnHAwcwbtUwsdr7yMeXvXBXIXlNLtDnbjkQePDBar4qpBDawERKXn54SmDKc6yX9M3K3dz45jwmL07zvE6w4/W+j5Zw3es/heXxHDoOTMVZ26oruYnV6LJ+pt9lOp/agHF3diM9M4fr/z2XlWmHKzCHqjLSAAZv3S39ndiKSmDKfkXwrRvsDt2tGqawCqnMnx5Zvm0Kd0mR14tToKqrpMSib80tf8n2/DK1gXE5YsrzHeTGYKC61Q4At6a7B4LPTV3H6O9/KTYtWAnMku1W6UtuGMbRKV8vJBWKnGo1WX76RZyxcyXs8j/wYpdWDfn07l9RPSmBwe/MY9qavRWYS1XZaABDURASqA2Mvwumly7YwST4GroGC2DisATGrZtyWNL1mF6g7ycpoejwP5CZw08biw9kmJwYeiNeCfBsqvJ8B+EqgfEJx3W6IMjv5s3UTbw+c2PxdSrweC1TCYyvCinMeakKlrW9hOPJNeCHxwIeYG2b1uGze37F6U1q87v3F/Hsd+vCfnyrqkEDGLydrPwtU9RQsLStB7L4ZU/xroMFBYY3Uzdx5Hhu4bRAFz0n1wAmBu8YnRcOE6ESGK/pBbrjT3aUwNw5dhG3/scaHdY3tVpSGQayCzCvPO0qwjWScThL6gId+/5UZHfwsn1S7P2e4sWJ5BQWnNMPtv4ISz8MuGzTuil88ocL+XX3U3lr1iZu/c8Cdh/WHkoqNBrAEKYGey5n8d6jU+n38uxi02at389zU9cx6svVhdN8BQFuF2Vn3tyClVgcptsZaPk6z4S/G7XHvAS48CcmlD78jSmq/Ekux0i8bkL9CpztQLyWtE1ZtpPWI79h/1H3Xm0/hLHIvqiKxnsIE2g7MnPyOHzMCuzDEZiXJcDyZS8G7wviwurW3aDVRfD93+DQ9oDLpiQn8vT15/LS4I6sSDvMFS/OZtLCHdqAWnkWdCTeqsBLFZLfde3/XtvA5NgPNXQ22PV9rtu53TnNLQjw3dHGUjfqYnm2v9xwlxR5L4HxXzTtLIHxcSZb2AbmRCjdqC35BYbWI7/hr/3PdKRdOs97jxxn/9EcOrSoV2qes/Qu12MR+/gF1kVj475MmtSp7jnfZeHbnlB+N4GOg5emrS9arsBQUGAKH1ZaFkVDJLmnsfVAFjsOZtOzbdEIsXrtLCdJ4MO2PRi0cxGHx/bn856/Jz8xuXD2PZ3uKbXK9Z1bcv6pDfjr5BX89dMVfL1yN8/ccC4t6teoyJyrOKQlMA6hBAGLt2WQk5fvOpJrqJ9VFMAEbvjpdvfqC2qC3aDvPnzM7115uDnz7Hvl78LwvzV7y9Qbwbm9efkFbPPTkDRgI16Xi2OBMYVBiO8zytKI1xdwvPhD0UXZLSu9XpjJNa/NKT2D4gGW16qXxADHUrgVFAYw3tcJVCKW43hi+aOfreS0R78tc94geMlq79Gp3D7m52LTCseB0UimzI7Ubsz0Ljdz0qGd9Fr2OZjgwXerRrX4+Hc9+Md17Vm4JYNL/y+VV6dvqJCnlqv4pQEMjq7QAUpRjIGfNh4gL7+AjfsyufHNefzr67XF5pdV4QXTLUAJUoXkmxasiubCZ2bQ7an/lT2TIXBm03ch8Fd18Nv3F3Ht6+4X8CCfUvjq6W/X0euFVNcBsgI24k0sffg7F/etG8rFzFed4hpsuqQTaJC8ZEf+vDZyTEwo+vzWI7/h3vFLPK1XFoUllyFEMIFKxJymr9sXcn4KCkyxfVWuRwmEvKZy2trsHH4+6zLO2rGUX6361tMJMiFBuONXrZn250voe9ZJvDhtPZe9OIvvVu7WgFK50gCGonYagSzbcYhb/7OAF6etJyPrBGA9xdr3VN1QT5YLtmSweJu1rq80xu0u25msWz7zCktgIv8Dz8sv4JZ35jN304GAyxUUu4iUnhYOzq/Klx/ffnEKNACcWxVSgaMNTFkelOnbl4XBj588l3THez+XmuYsIfLaiNcXwAy10/tmpfvAYuFQlmZNgQJKtxsIr21/snLyOO3Rb3nD0evJF/M5d3OBHdi1HvmNazqhlqQq/xad2ZcVp11Ip00/0WPNVM93eS0b1OTft3bho991p1a1JO4ev4QBb/zErPX7NZBRxWgAQ9FJy8uFauO+zMIf0SI7AIHAv0239guHsnO58c25QPG75qL0SucpUCPesnRSOZ6bz94jxUstVu08zJupm1xLdNKzTjBvczq/G7coYLrFqpAqoBdSoK7sgS6A7o14S68bSg9PXxZyXYKfQCffWev3B0zXa8lFYoDSkOKfX/79Ee42MG7JeG3742vwPG7etsJpvmPY2ch4w77MgOno9TGMRJhz7jWsat2d8zfMpu+STyA/N/h6tl+d3phv/ngxz990HumZJ7jjvZ+5+e15zN10QAMZBWgAAxSdtLxcZA3+nnNTNHVF2iGmrtpT+P6oo8Gu20f47jyLBzC+dB2f4ToOjP95wdw5diHdn55ebNqr0zfw3NR1LNl+0M9akBVkZNrijXh9+Qs5ewEVq6YKMBpyoLYjbidBZxuY8pRu+ap8nKsu2X6ID+Zt9ZyG82Lv9VECgYKJcI8ZVJY2MIG2wy0Zr3n2lVY5Ax7f/nPmL1heC79yvT6GhyQwu+MAFtjVSYy7Fo7s8rx6UmICN3c9hZkP9+bJAe3Zlp7Nr99dwIA3fuLrFbt0/JgqTnsh4RxszUMAY9yDEOek617/qdi8n7ek079DM79pulUhueXJ7WJc2MunDBenuZvSS03z3cluPpBF19YNi83z+hlu7RAqYhwYtyqIQJ/rPtico8SrwBeEhJB3e2e67as/fLgYsIZUd+t1VCovjnOz15IIl2Y9hcI9BkuwgezchHqceg3cfMk6R/DNdykhCtblW6QDb0kAACAASURBVKuQIkCExWddyuHajbhi+dfw1sVw/TvQ9jLPSVRLSuD2C1szqOspfLokjXdnb+a+j5ZySsMa/Pbi0xjUtSU1q+nlrKrRPU7RBSo33/DYFyv5zUVtOL1Jbf/Lu5zk3C5ydVKSOHo8j6wc/yUW09fude2FVGAMiUixBvxuJ39f1UJ5uikbYwpP7A1qVgNgl8tj790+/3B2Lh/M38qwi9o48u5Mm3Lnz41x+QzXEpgAF0D3Xl+OdQur50LPe6AL9TWvzWHrs1cHTSNY8OomUBXSiQiN5htK771A2+EWXHh9CrfzN+zjK5Vcv7eoO3qg7wcc48B4+lQVio0tO3GgXnP6LfyIRuNvZO2pXZjb4SpyqtV0Xd6ty3VKciK3dm/FkG6nMm3NHt6evZknvlzN6O9/4YbzW/Dr7q048+Q6kd4UFSO0ComiE/G63Uf5cP72ID03jOvZzfU863LxLrnqXeMWFQ5k57zYupVcuF1wy1KFdDw3v1j3ROdJ33en73bhd7sof7l8J6N/WM9/ftxcOC3tYDbr9hwBioK9cFdZO4NI36tgVRBbDmTx7HfrHO1yXNJ16bYeSqGBbz94LTEJmFYZAphApSHOfeq2P/YdOU7rkd8w02MPoAKXEg6v63jltcTG9xtzBmm+dQ9mF7W7SAxWAqPdqCPqUJ2TmNzrXha3682ZO5Zyy/SXOHvrz0hBaN2lExOE/h2a8dndv2LyHy7k0rNP4uOfd9Dv5dnc9OZcPluSpl2wqwANYMDR68Q6+QUqZjbu8UvAO7ZgwUWCS9dbt3Y5gcaBOXo8j2e+XevpoX/nPzmN85+cVvje2UDUN+aJ2wXTbZpv+SzHqLHXvf4T/V/+0c6f/7yXhzNgLBqLpzTntt01diFvzdrEjgyrdCloCUy++8VsR0Y2rUd+w9crStfll6dKryRnGl7r+t3Gtrl9zAImL04rFlS5ZW/ZDutBiuMXbCs904XvewmpBCbE1uZeA0H3QSCLJu4/msMDE5aSdaL0U66dv08NWyIvPzGZBef045Pe93GkVkP6LPucm2e+Rqs960K+0xERurZuyMtDOjP/0Ut59KqzOJCZw58nLeeCp/7H2FU5LNqaoQFpJRWxAEZE3hORfSKyyjGtoYhME5EN9v8G9nQRkVdFZKOIrBCR8yOVr5IOZp1ghn3HWTiqrctyzeulAHY32yAny5KCVUFIYbuJopO1WyPeK1/5kb6jU4ut68vzL3uP8vbszaTuKH2CLin7RD7Zjoa4uXlFH+K7a3Ev7QmtPYkznYj2QnKZ5uPMs686IlCvM7eAseR2r95llS5NWeYWwFj/w9HexC2YCsZtTJYfNxzg4U+Wlwhg/AeoSS69s5yDzJXMXyi7NtAx5FaS47kNTJCA+7mp65iybJfrPitWQhogGC4oMHy/ek+xYF2VXXq9ZnzW8w9M7XYriQV5XD1/HDelvs7pO1cgHga/K6lhrWoMv+R0ZjzUm/G/7U7fs05i7u48bnprHpe8MJMXf/iFLQfcB7xU8SmSJTBjgf4lpo0Ephtj2gLT7fcAVwJt7b/hwJsRzFcxzt42vpOr8/xdv6Y1DHZKciLg64VU+vS2aqf/0WTdSlacfKftE44GiH+etIwHJiwtdeewucQPsOSFqCw1F7kuJTBen3zt+3h/7StC6eEVCuPyJliefY18C/Pkp1G0lGiIW3oxa0LawWOkHcwuPieMJTDOfe+1LUigKpJclypKJ9/2JrqMj/Psd+tKTStLcBqoO7h7LzL35X/eksFfJy8P2E3fuQ8CDUoY7PfpM3HRDn7/wWI+W5LmfyEVGhE2t+jAhEsfZGanG0jOP0G/hR9zy/9ehPlvwTH/vSH9SUgQLjqjMS8P6cwrfWryf4M60qphLV6buZE+o1MZ8MZPvD1rE9vTs4MnpmJaxAIYY8xsIKPE5AHAOPv1OGCgY/r7xjIfqC8i/rvthJFzmPiibpdFZ9KSA5IZ415c/dqMjaUnlkjDH99Nc44jgPlu1R6mLNsVtP1FqBfKkhdcKH6Xm2OPDOs1gPE55qdrdWFpR5h7Oxa76Nj/3T6jWAAjJZZ32Ryrl5kvCCmwlzNk5eRxz/jF7Dt6vPAit3b3ES5+bqbr54Wje2fxKiT37z49M6dYEB6okaozT24Xat/8ZJc01u89ys5Dx8h2VMGUZZDCUANZfyVZQ96Zx6RFaUUNrQNUrzrnByt1862SfSKf1iO/4VB20eCIOw9aVY/ONjUqPAoSElnbuhsTLn2Qqd1uJadaDZg6Av7vLPj8D7B9QZka0tVIEm7s0pIPf9udeSMv5ZErz6KgwPDMd+u45IWZXP3qj7w+YwMbg4wPpGJTRfdCamqM2Q1gjNktIifZ01sAOxzLpdnTwj6M6KvTN3Body697ffOC68vgHC2gSkocedmcL9zB/hq+S6a2VVNELg0w8l30XFrdBbshB9qD5lJi0rfPbqNneEawLh1Xba/Kn/PCwr3SLyCtQ+KX4z9X5ycF8CiZxx5u5gVlsAUGD5fupNvV+6hXo1q9Gzb2G/+CrvzlmVkQT9pWem5B0RD3pnPhn2ZbHnmKkTEcy+kgCUwLlVI2SfyuejZGZzXsh5f3nexPdV3rATbEsdnhDgOTLAqpPwCQ3KiezDqPF4DNch2G0DSZ3tGNvXtnnlJdsmUjj0SOUYS2NyiA5tbdOCek34Fi8fCyk9g+cfQ8DRofwN0uAFOOifkJ9ieXC+F3/c6nd/3Op0dGdlMXbWH71btZvQP6xn9w3raNa1N/w7N6N/+ZM5uViekp6yr6IiVbtRuR4rrmUtEhmNVM9G0aVNSU1ND+qAJc7NpmlJQuN7ybUV3U2m79wKQeeRI4fxcu+4/+5h195WRkcGyFUdc077/46XF3vsuCBs2/n975x0fR3Xu/e/ZKq16t6xiSe69G9vggoFQAwQCCSEE0klC8iZ575tAOiS5IeXmpkDCJSQkJCT0XEjoxQJjg3HBvcpdtmzJRb3tas/7x8zszuzOFsmSLVnn+/nsZ3ennjPt/OY5z3mePVT3aFmCtx6N7j9vb9G6n9ZviTbTr3rnnahp5jofPWaNpNvV3UV1dTU/ercDp4A0t+Ars8Ki6sD+/VHbW/nOu+w42cNjO7pJc2un4tDhw1RXW1MG7D4VFilGGfbs047f4aPRI1eqq6tp1IdjHzzawM8fe425I+wvudDxDkr+Z2MXHx7nYURadEPqFJKAFGzbuZvqLq0uHfq5WbNuHaf2OC3Lb2wIH29juXfefY9DGQ5274t+k1616h38fm16t19b9+ixY+zq1o5FXd0RtnQdsy07wJ492ht7c0sL8bC7bltbWy3Tt50IH+/122p4omU/hT7rMdldr3UpvvpGNR6n4Ghd7ISdq9eEIyhv276d3GbNahgISn77fhfZXu3cN9QfpbraarpvOKldo5tqm0JlPHxE21fNnj1x61VdXR2q20bT9R+5bG1tdNlXr1nL8d3OqOmGznjypTf5zfpOLhkVzngcvjbD1pOj9dr1efCQ+T1J48233ibdo9V9x0mrEF+/bh3/rA7iEnCqU9tpzb4DVFdb360iz925hvm5W1BQQEFdQYI1Tp/qxkbIuBbnvEspaHibwvoV5Kz4JWLFL2jzldJQcAHH8+fSml4Fwr4zId55GQuMnQQnq1JZd6yHtUfb+e3ru/nN67vJTRFML3AyvcDJpDwnHptu1f5mKF1Dg6WsZ1rAHBNCFOvWl2LAaPVqgTLTcqWAbbhGKeWDwIMAc+bMkUuXLu1VAbI3rkAG2jDW2/HmHtiuCYf0rByoP052dhZLly7UVnj1RSBIamoqtLeTk5PDxEnlsD75JHkVlVUsXToagLZNdbDBum55cSFbTxwlZ0QZ7NhjmXfeefPhTWs3hbnOjx5cC8fCDarH42XmeedT89Ir1uVf0nK/VFRUwJ7dlu3NnjOX7/3mbbp7IMXrAvwUFRWzdOk0y3Kpe0/A6nctZdgp9sDOHWTl5EK9NRz+Hcu7GFOYCU2N7DgZZMfJLl752nmMKzLFadDLZWzv7d3HWfvKatxpPr6wdDR1TR18ZG45bV0B2roDOF59HXpgVGUVs84rx+tykLLmLWhvZ+bMmUwszsTtdIQSIfq3HYN1WsPtSUmB9g5mzJrNlJIsduhltxzv+fPxrF0J/m79oRgkv6CQsZW5sG0rne5MtnY4gXBdzedjS3A37N6FN9UHrbEdBkPrvBTOyZOenm7ZlnN3A6zRchq9sM/PC/v8ofgxje3dbDvSjBCrkRJmzFtAYUYK1c1b4eB+231Omz4TdEE8fvwEls7Rbrndx1rY+MpboeXKSkaydOlUS/mEJwXosJT9+YaNUFvLqIpK2L0r6liYz211dTVLly6lddMR2PB+9LLAyrZtsH+ftcwzZjI3IqAigHj5eaQEb/FYGjo28ebRsMgxtvte5w7QxVVObh4cq6ekpBQOWI/PgoULyUv3AuDZcxzeWx2aN2/uXH7wG21E3TcvmwC7d1BcUsrSpZMs2zDqd65ifu6OHjtaNhTHT33RHzSY7jHyS2DiR0ntuoqqI1sYU7uJUQcep+LAY7R70zlQNJ6DReM5VDiWbnf4ha2go4CD2Qct27WLL3Odsc+WLpbvqOf1HcdYsfs4yw91keJ2cP7ofJZNLGTZhEKKs1IHorpD6hoaLGU90wLmOeBW4F79+1nT9DuEEI8B5wFNRldTf+NxOQiYXrzN5mPD/0OgRaRN97rC2Z5D8SF6P8IkUfdJulc7DSdao5MR2q168EQ7X338fW5dWGHbnTX97lcs/1tNoybs881IzXG5JzwcOl7ma7tp1TujH2itXYEok/ye+taQgLE4qfYE2XCokfQU7Vi8s/cE7+zVIgV/ZG45H7zvbfY2tJGqX7H+QJBpP3iF6WXZ4WB5QZj0vZc5f0wej35mvj7N5N+gW/7jpQgw+8CYlzOWfHdvpFtXxPGI2Ecs7vnXNr595cQE24q9jc8+sjaUSBS0YfSFGfF9YPwx4sBEjlxy24TzNY9UA3h45T6eXKd1R/ame9Bcpw2HGplRls1PX9rB76v38LnFVVHLJ+pCau7QbmaPK7rM5uv1te3au5KtE68lKmK8smsntz9i/Cj6Roc3na2V89laOZ+UrlbKj+1i1LGdVNZtY+LBdQSFg/qcUg7nV3EkrxJ/T++C2hVkeLlxbhk3zi2jK9DD6r0neUMXNEaG9InFmVw0oZBlEwuZXpqdMDiiYuAYMAEjhPgHsBTIF0LUAt9HEy5PCCE+DRwEbtAXfwG4AqgB2oFPDlS5PE4HnaaHqNlht1PvLqo91cGU77/Mjz80Jdx3rj+zJLLXjrNWh9PodY0ynLDJpmzXOPz+zRrWH2xk/cENLBmX2JT7uUfiJ18MBIP60NlgqJF7al0tv7hhOgA3PvAOI7NTuH52qWW9nqBMGOMmcvYXHl0fsiKYj8vPXtrBH1bs46fXT7Xdzt4GzZphnK3/elV74994qJHyXJ9leytrwikS7Eai2OUpCpfXxm8imLz/oLF+oob3Tyv3cdHEwrjL2O1z//E2HnnnADuOal1UDqEd4xOt3YzKDSYQMPF9YAzssnRHCrK7/7UtvK2IeRV3Ps+tC0bZbtt8XK69fyUPfWIOv6/WrCR24iLWKCQhBEgZEjC2mcVtR5lFb8sSUyhinvX60b6VgBkcdHrT2VU+i13lsxDBHopOHWLUsZ2UNOxhxu63mL2rmiB/oWHPSA7nV1GXV8mxnNLEG9bxupwsHlfA4nEFfP+Dk9jT0Mrr2+t5fUc9v39zD/ctryHH52bxuAKWji9g8diCkCVPcWYYMAEjpbwpxqyLbJaVwJcGqixmPC4HpsE+FovEKX3EwVE9Q7M5ZkQgGG70ehuMy/zwv+Pv70fNNxqTpo7kBIx5tFIyb7/mnEd2+YL8PdI2yV1dUwfFWam8t1+zOnxw+sjQvJe3HuXzf13HNTNGRq+YoPzheeHf2+o0v6KGltg+HBAriKDV0dqMuQE03rQNkRargYt03tNi/yR3zk/qIjSZsP2JhLDd/K8/sYH1BxtD/x1CEJSSG//nHWaVZzOrPCfm9szHIp6DcJo3+rEQb/izXRwVc1ZogIPNPWw41BhVp4Mnw6PibIMnJrjXmnQBs+tY9CiSeLnDzJjLHzn/g/e9HV5OPwaGpVYxeJAOJ0fzKjiaVwGAK9DNiJMHGLP/GNldW5les5JZu/Vu0nf+BiWztU/pHBgxDWKkMzAQQjCmMIMxhRl8fslomtr9vLm7geqd9by5s4FnNxxBCJhaksXScQUsGV/IjDJlnRloBosT7xnD7RRWAWOaZ0RorcxPY9/xNlpNWaTDcSTgD6aw+cmQrJWi0+bBaLeqWcBENgiJktHZdyHZv7kHpTXCrhHADeDR1Vq/cqLhh/Eaacvw1VDE3ribs0+kaXQhJRkLJBAaHh2/TOZp8Y7qr17bxcMr97Pm2xfz13e1htsc0ycWiSSRXVlcEd075kZ6/cFGZpTFFjDdAfuGOrJ7yK4LyVyfSDFnPmf3/HsbdnxvVSesWslPrrNa2MzXo925TWTtMASMHXb3nd0xt0TijXNSjGPd2OGnpdNPpz9IQYZ64x6MBFweagvH0tWzkIbiC3AFuilsrKXwVC0LpRdq18LWZ7SFhROKJsGI6TBiChRNhqIp4Iv2vTLI8rm5evpIrp4+kmBQsuVIE9U7NUFz3/IafvNGDdk+N4vGFrBUt+Koa6X/GXYCxuNy0Gl6YNs9r4yHqjnsuLmhMCeHS4ZEQ52NBiHZYdTmN8D+CJgW6JG2Q2f/uGIf180qCf03NxbbjmijUtJt3tYt205SwET6GsXCLi+msYpdYxewEzA9sfdlZ2kJSi3mSyx+9ZrmFP13Uwj+eA1rrH0dawvy4uY6Lp9aHLN8iY53vONnjfQsOdnWjdfliLIW2V1T5mszUmib9/nwyv3RZTLHs4nYtrkL185qlej6boxznO3uOztRY/GTinP8jHofberkkl++xdHmzqSScirOPgGXhyP5VRzJr2Kh4cTbcgyOrIfD67TP7pdhw9/CK2WM1AWNSdTkVoHLY9m2wyGYVprNtNJsvnLRWBrbu3lr93Gqd9bz1q4G/rVRs+RPLcli6Xitu2lGWY6yzvQDw07ArNpzgpZOyZbDTUwpybJ9QBoCwWKBMfLiIFk8rsDWaTUWiS0wuoCxCdduJ2rMYd17K2Bs33KDQWxeuvnTyn28uCXsS21OP3BcdzhO9IZsFzMjGJQ4HCKUwgEIhfi2ExDmbqV4tTUsaJH7MjCOlVHmRLmQDFq7ApayxqK3WjJy8e+t6qBrRdhHyG57dt071jLELoT5XEm0nFgl2akhXycD26BwpkntEfmEnll/OG6ZLCIy4nowOxDbdc34gzIUqynV46SxvTsUlwXgeGvsLke7S9Pu+HT3BLnrmc18bnFV3OurQx9Wf7ixIymBqhjkZBTB+Mu1j0FrPRzdDMe2wLGtcHQL7HkDgvo173BBTiXkj4OCcZA/XvudPxZSMgHI9nks1pmtR5qp3llP9a4G7l9ew2/fqCEvzcOyCYVcMqmIRWMLSPVEhwpQJGbYCZgWXZQ8vuYQ3/nfLcyvyotaxhANZqfasIUgcffAqDwfB0xhqhOJDGOu3QP8aFNn1DTzcr0dEWXnyxDokbb5bwDqTPuPbLggsa+HXUA3fzCI1+G0+AOd0o+13Vvzlfow1lgYouf7z22NmmdrgUkQlTXyzciu3nYkso5EIWFMYXqoG86wLnUFevC6nKE3NzPZqe6oaWbiC5jwvCY9muzhxo4oEfrr13fzqfMryfLZ76s9IupyPBEBERGFI65X85G2y7cU6Aky+fsvWQTU3z97Xmg9u5F7Bsk68W6qbeIf7x1kx9FmPnNB9EgoA6PeSrwMbX634XeJF0pzQ9UMqJqBoydATmsDeU115LQ2kNPSQPaRNWTtehGnKWdTa0ompzIKaEwvoCktj6a0PJrT82j25eDMc3PRAlg4y8H+oynUHPbx3KZOnlxXi8sZpGJEJzPTnMw8z09WgntcEWbYCRgDw1fBzopnF1U29NCWia0OPo8Lt1OEGoxEXUjGQ9Vuvwds8nWYH/RRYiqBnrGz6AR6gsTQLxbMFqlQWRI4NB5ujLaK2Am6QIS4MFOfwLE33vG184EJW2Cilw8GozM6J+u0mYzjrmVf0n4U1yOrDvDZxVW8uOVo1LzCBP3osfSsEFaHWGMUF9gL8uU767l2ZknUdIgddTkWfovzcOwupC6bcgR6ZFSdjKzZYC+eugNB/rXxiG3+qLomu+sxPDz6S3+PHd8pUrgphgdBp4sTWcWcyLJmt3EEe8hsO0lOSz3ZurDJaa1n3KH38QbMVmNBa2omzbqoaUrLpWlUHqcm5LGtYyRbj+ZQc9jHk4ddPPuj17hwQgHXzChh2YTCUA4+hT3DVsAY2MWPsHuQGkgk3QlGRnhcDktXTbJdSHb5hBptRiaZyxfZaCayx9g9hP1Badv9EonRWBhDd+32nwzxwuz3ZYTHsebYAsdigYkY4hzLAhOZFTnZBttOHJqPVSSf/staSrKjg2L9+IXtfNYmJgokPt6xRks5hIgpvO2mp7hjK9reNuRm4RTtAxP+fbKtG5/Hadm+naBdset4aLrdtfTIO/v50fPbbctiHmJvYOQ2SpSvy+7+lFKqkPPDlKDDSWNGAY0ZEaEspCSlu52sthNktp0gq+0EWW0nyWo7QcXR7fi6rAMf2r3pNGfl0JhezInUcby9L5Unt2fzB08R582YwUcumERlftoZrNnQYdgJmH9+cSEf+t2q0H87sRKvy+fwqQ6yfJ6Y8wG8TodtTp1YRAZOM2P3dmwus12jGQ+7xvhkgi4AA8PvxedxhYLj9UVwxDu+8cRjb5BSEgjK0Nu1eb+BYJBgUPKvDdFdNFJGB4OLlagyErtz4XI64nY52lmo4rGnoW+jvnqC0raLDexFkdcV+83vWHN0t2Y86lvCy/dECA7zv3UHTjG6II09DeEIxnZdnkaAw1i02FgK42H4WCVyILfrSmzuDCiTv8KKEHR60+j0pnEstzxqttvfSaYuaAxhk9l2ktL2fUxo28D5wW4wmpgN0Ph+GrUpI8gpHk1aYSVkl0FWmf5dDmn5vc4Lda4w7ATMjLJsy3+7bpF4HGnqxJfA18HrdljeuhMmZIzToNu9YXaZGsraU71rAO0ER7IPfMMCY27g7fwWEtEdCPK9Z7fYzuutIIvF/7y1l3tf3MHnl0RbMvw9kqpvvWC7np0PjF2AQTvshsH392Plhc3R3Upm+jIoze6acDsdMa05h05Gd2vG47JfhX2YIq09dz2z2fK/OCuVfcfbQvX43rP2oisevQ00ZwiYRPdShz9IRorLcr9Mv/sVNRJJ0Sv87hROZI/kRLY1hlZBXQE3XHo9tB6DpkPQeJDWY/vYt3s7TUf3UbxvO5WHVuLpiUhR4kqFrNIIYWP6nTESnOdmU39u1ioOQghmFjp5v15rKM1xTpIl0QPSEzGk5x/vHeIn102LsXT8YF12FonOuF1c8bGzwJj9IeJhlMXcJdYXi8mRpg4eiQh0ZvDP9+OPaEmWJ9dqSfvsnKDjZRMOStln0WF3bHvrZG3Q2ofrEvqW9bvNxrIQCAZjlr3O5pgmS6vdOHgTOWkeLdTBaQSLa+4MO9k6HSKhE71hIYo85h6X1XrW0R0gM8UdJfj7I5SBQgGAwwGZxdqnbB7pU2HmxZrIvvPpTby+4xj/fXUFH6oM6iLnUEjs0HQI6jZBuzUJL8IJmSOt4ib0Xa6JH/fA5HcaaIadgAEwRx1v7qUFBhKPQrLzq3lrVwNff2KD7fLxGjnbLqSIhjLFHfuBv2RcAW/uCg/53pugCyIZzM6RyQRsi6S3Jv6+YPhR2PlrxPMj6QlKW+fPZLCzHvW1cbvi1/FHXsUikb+VHea0AAZNHX7GfvtF2+XtRGGy/Gnlvrjz3U5xWuIFrBYltzOxgIkV/flDM0p4fG04e3V7dw8+m+Guo7/1AvdfFD+Sq0JxOhRkeHnwE3O45Jdv8vdNzXxo4UIojvFS3N0OTbXQdNAkcA5p0w6sgubDICPusbRCq7jJHmUVOvoQ8cHGsBcwfXnTTWiBsREwn/jTezGXjxem/en1tVHTIq0e6V4XnX6tm+OpXdYhnqfauynI8IYe0vttRjVF8sDHZ3H732KPxjDvvy8WhuYzMAzVsBLY+S3EswJ8+IF3GJGZEnN+PB5dfTCu024scnzukCOpwcGT7eSleTjZ3p10HibQusfKc300tnf3SZwb/OcLYSfYC8cXsNwU98huJE9/ETkCrC+Yr0/zvRELs8+NmciYOx3dPaFko5Ec71DpBRSnT6wh3t1+wZqdmew9nkV2Tgu/22D/MhyFC8jL1z7MBEAEe0jrbCajvZGM9lNkdOjf7Y1kHNhP+o4mXEHrs6PTnUqLL5uW1BxafNlkB8p4vvNpmtLy+diCu6KC+50phqWAsWsPPj6/nL+9e9BmDpxXmcvqfeEsxJFWh0jHQ3MX0s3nlYfC7sfidHIrAWSlukMOtpF0+ntwRzQKI7NSeOTT55GV6mbuj1+LWufCCb1PMtgb7nujJvR7flVuwgzPZpwCenO47EadJPJ7Op1kfYmcdu0ozEiJEjCg+d4sm1AYN4jeqjuXsfDeN0L/n99cR2GG11ZE9wbzyK5lE4ssAsYu71B/4XI6+H+XjufnL+/s8zaMnGYAuWmemPdGIiKtLS1dAYod9uK2sUt1Iyn6l2AQDjV4qTnsY8u+NLr8TsaXtbF4emPileMgHU5afTm0+nKoo9JmgSCpXW1k6qImveNUSOxktZ2g9HgNnsAqZhjN2hu/0rqhcqsgdzTkjdFSMxROhvTEyYZPh2EpYFq6rQ+bbJ+bkuzYJuDIsfiRVpt5lXkWAeM1DUFNJrhZX/0kDOL5oXQFglH5/FztRwAAIABJREFUczwuB2MK02P6ghgjUNK9rqQsVCMyU0IJMO0wBJ7HqYWt32lKxZAaJ87B5JGZNLb7LSN13A7oieNGEa87zSCR31Mip92q/DT2Hrd/a3c7BPHWvm5mCc9E+PkUZnotx8SMXbhxsyOpXZdGZX4amaluXt12LE5JksfndnL7ktE88KaWNbqv/jnJ4HYI5lbEzkGTDObrJTu172+GdvGFYom3fho8pxjG+HugtsHL4eNeahu0T5ffidMhGVPSztzxzYzM75sY7xXCQUdKBh0pGbajqJCSkkOpBNJ3kdV6nIszquDEHji5F7Y8BZ1N4WXTCrQ0DKVzoXwBlJ0H3vR+K+qwFDDGC7jX5aArEKSx3R837kXkvEi9EdnGeJzhRiVR6HfQho+eDnYWkdmjclh34FTMJI1gbRyN5Q32/eQK9h5v46L/ejPh/k+a3nhdDhElyD6zqIq7ntlMeoorlK05GTr8Pbid1vInGi2YkeKm0x9/WHhDksPGAb571SR+GJGgMFaEWrAXHGa+c9WkKAFTkB47OF1k/SP3Ybe/C8bk8+lFldy/vIb7l++JW55k6AnK07bomHng47PZdqSJ/AxvaJSREb3a5XSQnx5bdJTmpMYcLZSf7uV4axdHTAIm8zSGOCcbIiA3zcPsomH5KFX0gWAQGttcHG9009DkoUH/bmxxIdFepHMz/Iwr7WD0yHYqRnTicQ8iC58QdLszaMgt51huORcbuaUMWuu1NAz12+DYNji2GVb8EmQPuFJg4tVwyd2aY/FpMizvOqMLIsfnCVkO4sW9SBQNMbIRMTc68QTM+KKMmG/e1v07mFuRy4rdx23n33XFBEtYfghbflJczigfG8MvxByAqyIvzSJghBBU5CUXPMncZeK0ETDG4RmV54sSMPGObVuXNurD2EZQ2kdONpNMzJbILpkXvrKIF7fU8VtT15bBB6cXRwmYeF1o5kzOT39hIdf/Phxz6CfXTbUVAvGy1HZ094SiOk8vy2bjoUaLU6pdCgif14XP42LZhKKkBcxnF1XyhxX2DrbdPUG8MQRMYYY3YaTkSJaOL+CyKSN4cXM4z1a2z6MJGIegOCv2iAg7H5mfXT+N2RU5tHQGuPb+lRYLnMfVe5+aRWPzWbH7eNIhAlROPkUk/oCguc3JqVY3ja0uGltdnGp1c6rFRXObi6A0LhpJdnqAwuxu5uRCevkpSvK78HmHsEkvvVD7jL4wPK2rBQ6+C8t/DJufgK5m+Njjp72rYSlgjKGm5q6eWA9o0ERAJJ88vyKUeTcycmujyZ8h3Ru97neunMj8qjzWHzyVVJyL+26axbqDp2IKmNy06DdWo4nzuh1sO2jNpGwWGF+/ZBy/fHUXqR4Hs0flsNEUpt0szO65ZnLcshrDVY3urAkjMthxVBNnxu7GF2Xw/kFt+0Y3yNTSLNuQ+aU5qXznyok8+NZefXk3TR1+EtkBrpxabBk5kgyTRmbi8zhtBYxdF2C8dyFzcsLirLC/xD3XTOameeW2/jGFcZyGl+9swOdx4u/pYcnYfDYearRYBuxSQBjdStlxLEWRlOti1c7C0R0Ixkxh0BenW+NeM+vcNL3MLqcg1eOkJDvVNshf5L0GMLowndEF6dTUR78MCCFC2/rApCJeSaJb7Y4Lx7Cy5jjXzSpJ6L9m7EMxfAgGobXDSXO7i+Z2Jy2R320uOrqtz32PO0hOup+inG4mlLWTneGnIMtPXpYfj0u7EQrqCmgoHjgH+TOClNBxKjy0u/EQnKiBo5s0q4y/HZxemHpDv+xuWAqYgP7g9HnC1Y9n+rfrMrhwfCEPr9zPyKyU0EM12+emvbuHC8bmhxpRs4XBaLRHZKUwpSSLZzckF/PE5RRxHUMNK4WZScWZ1Bxr4c7LJvCRB9+1zDO/wRuNicvh4MnPL4jZONtZHe64cAz3Ldca/dKcVEveJiM6aVVBWkgwmh/0v/7oDHLTvIzMTuFnL0U7bL79zWUAPLVOG4WVl+ahqcOfsAvppvPKowTMt6+YyI9fsA8tb1ARI1S3rY+OlHxl2Rh+YyN4XA7By19dTHOn3+KfkpemCQC7LqGynNgWB7dTUFWQxpbDzaHuEPMwcFsLjCFgbLpPLp5YxGvboxvxScXaMMlTNl18QSm5flYpQgj++X6txTHa63ZSlpuaVCoKA+M68JnEveGnZdRn5Z3LeHzNQb75tDXQXV66J8r/yLiG7aydTiG4alox//PW3riCLjPFFRq1Na4og70/uTKmr1TkSDNlgTk3CPRAW6eTtg4nbZ1OWjuctOr/W03T2zqdSGk96V53kExfgAxfgJF53WSkBshM6yEn3U92RoBUT3BoBsuVElePH6+/A293B6ndraR2tVLQIKg8dYzUrlbY/ja0NYQ/3RF+Yt5MGDENZt8GI2fB2EsgNdt2d71lWAoY4+FTkOFlu27FlhFNt9kxtdSmgSnL9fHLG6czrzI3ZImZVprNI5+aB2hv369tr7dYdt656yLWHTjF4rH5AGyrC1tG5lXk8t7+k6FugkjiCZjJIzO5++rJXDypiPP1ESmZqS5W3XWRff1NT19DuEkpLdaDSAw/IPMb+hVTi7lveQ2XTCpiV0RX2PyqPH5xw3TKcn2hxJnmG3jZhCLAvsE0k62nbRiRlcLe421MznPynRsW8MiqAxahcvHEQl7bXm+xeoTKrjfoYwvT2V0fewTN1JIsNh9uskyze7sOSvjqxeOoPdXBkaYOyyiq7kCQ8SMyAGuU4jzdr8Nue3ZdJnNG5bD2wCn8PZKHb5vH6n0nbEWk3SkzhHmOTcqL+2+eyfjvvATA3VdPDqUXmDxSEzBZqW7aTN1wt8wfxcfOK8fhEHx4dinLI7rf/nTbXPLSPTy74Qjf/d/o6MrfuGy8rUAFWDqugJLsVAoyvKHQBmaLTlaEA+700izuv3kW8378umW6J46Ayfa5+cZlE/jUBZX87V374Img+c8YAsZ46YjlX7R63wmCvRw5qDg7SAld3cIqRIzfJqHS1umks9uuO1uSlhIkLaWH9NQeCrP9pOviJMMXINOnfXsHk48KhISHq6cbd6Abt/4d+h/owuPvxKt/PP4O/XcHXn9HaJ7H32HJuG0mKBx0eNIgq11LZ1AyW3PaNQfJyy6H1JwBS3UwLAWM8QK7eGw+b+lB3iJDsaSZ3g6r8qO9pouzUqicVQqEY7W8ZQoYN6cilzkVubxnGn6d7nWxZFx4WNkNs8tCb7NP3L6A5TvqmV2Rw9r9J5FSS/ZnYDfa5KIJhcyrzEUIwa0LK5IOmmZezLDetCQYWTJ7VA6gWTNe31HPU+tqKc1NZfs9l+F2Ci76pdXZNyglZbn6yC695bVrbJ0mi4Rd/JmiTM1yke1z89rXF7N381omj8zipx+ehhDw2BpNxPzu5tkEpbTtCpxVrqn9zy8ZzX88udEyz2xhWTA6L0rA2CHRxN4vPzKDWyPi+3x4dmnot3k4vV3SRmOUl2GtuqjcxesHtfPwxQtH86k/a+e/IMPLVdNG8rqN5cROEJXrx91OkHpdTq6aVqz5i8woCQmYFLeT3988iyklWSz62XJAs/788NopUWU2YySZG21jwXri8wvIT/fEFDBCCFbeqVnaPv9Xra7mEXNmi8lLX13EqNw0Um3uA8PvyGdjLfvGZRNwOgRFmSkW/6RIzJYdQ6ybxdR3rpzIo6sP8tCtc5h29yuADFkg+xLMUXF6+ANCt4Y4LFYTy6fDSXunk0AweiSN0yFDoiQ3I0BZYRfpKT2kpWrT0vTfad4e7pj1RZsS9BNSat0q3W28W/cGN4xYrP9v1QLSdbdpv/3Gb/0TWqYtvJy/zfS/lcRx2XVcKZCSpX9yIKsCUrJN0/RParYmUNIKeHvDTi646CrS7PqwzyDDUsAYDfiV04pDWWuNbo5Z5dmsP9hocUysKrA+nG8+r9zSNdRoE8PDIJ6D5rUzS/jq4+GAREb8FcM6YSDR3vgLM1OYV5HLB+97G9CsHOasxeY3xnX7ww65ty2sYPnO+lAXj1noGCKjKEHwtjGFGez60eV4XA4un1rM3VdPtrzxLhlXYEkPYPazMX7a+S+YG4nZo6KHzxqOqHsb2hhTmEGtSfDcdcXEkIAxO8ea44h85oJKJo/MCuWr2XDoVCjez4pvXGjpfvvGpeNDPjc/vGYyG2vtxYzZEmI+ltvuudTiDG4WFyExh2bdu352CQXpXn7wr23kpXvYf++VVFdXs6pO0uHvYURmtOAxd3ma+dn102ho7QrVebTpen3z/y3l9e313PPvbaFjdN/HZun1kNw4p5Qb5pQBcPnUYst2d/3o8qh9fe2ScdS3dFriwkB0N+vcihzmVeaG8ib5XHDXVVPwxhARxrVk7mIzi6XyXF9IvPzs+mlU5Kfx+JpDPL2+NiQAzeLn0c+cx6zyHIvgsXsJ+MwFlTz09j6Ltcc4b06HYHppFsdbu/nU+ZV8ZpF2rzn1+ZdPHcF9y2t6nZ1bYU8wCO1dTlo6wuIjJFAihEl3wO46kvi8mrXEl9JDToGfgh4vIr+VNK9VnKT0tksnGAyLCH+bVUxYBEUi0dFmXc7fjiE05gOsTlAOVwq4feBJB08aeHzad2qO6b8+L7Scvow7TV9G/6Rkad077t4H7gy46+wd8M4ww1LAGNZfc4NvNErluT7WH2wk1e0MxdowuiVumlfG5JFZfGyeVdHHi5diWBD6QsgRVkKqx8mnL6i05HmZWBw7vPNXLx4X+v2DqyfzAyZTcefzgCbcDBaMzuOhT8xh0bj8hOUxi4RIc/13r5rEZxdV8fzmOu59cYelmyrkA4Pm22CeZ34rtsvqO6s8m5vmlXPtjOghd1mpbtZ/9xJqIrqFvnThGK6ZMZJOf5AxhVbr2XeunMTf3j3IvIpci6gAawN4y4IKbtF/P/2FBWw70sx3dSdmc5ySK6cV83aN5lxtJzA+MKmIeZVWYfbut8Jde7edbw0kdceyMfz85Z2U5WoCxtyI23VlAtw4t4yVehnKclMt9RiVl8anLqjksikjovx5hBD87MPTbbdpzI9kRFYKD39yXuhaMpg8Mos/f3Iutz28BoAnb18IhIMCZnsFt8wfFXNfGXo9zfs0uuLA6kh/41xNcM2tyOHrHxhncWL/9hUT+elLO5hWmhVlrTG6I6+YOoKvXDSWju4eDuoCyy6MghCCZ++4IGp6ZqqLDn8P+frw9/7KoH4uYzi+NrW5aOnQHF5b9W6clnZtWmtHtG8JgNfdQ1pqkDRvD0U53aQV6xaSlLClxJcSxOftwekARzCgd5F0U3g0jdacOtyBbjyBLtxN3bhOdIW6VCzdK6HfXbh7/LgCXbgD3fD8D3Wh0QtcqfaCwpcfU1Bs33uIiVNnWUWHJ10XIvr652hSxr4yLI/GF6Z7WdWYHnKsBLhs6gj+tekI/3HpeBaOyWfOqByW6TFQhBDs/vHluBzC9qH+3B3ns+y/3rR10PR5XGSluqMa0mS49/ppfOnR9czSu29A6/LZf++VHG3qZISNv4fB1NKsmPN+cYO10bp4UlGMJbU3+GRyF7mdDspyfSGLitkCY4hDY0RI5Hr//vIFVBWk2R4/IQQ/uW5qzP3mpnmiBAJAaY59YMIUt7PX2YNnj8pl9qhcLppYRH1LV8hfBOCjc8uiMiqbefATc3q1ry8uHc3nFlfhdjp48JbZTBgR3tdIm24oA8PfJVZU53jrRrL6Wxf1KaHl0vGaBdFsVavIS+PWBaOY5IodTRgIheg3O86axa1dd5jd9fTZxVUWq6QZQ0h7nI7Qcd2tB6aLF0Yhkn98dj6vbDsWGpn1xaWjgfhZwocDUkJTm4uGJjcnm900tblCQ4ib210EgxHhJlxBMnw9ZKQGqCjoIj+1nSJ3M3nuVvKcbWSLVjJFOyk9nSFh4dFFhbulC/cpXWwY4kT/7ZTJW8T8Tg9+lwe/00PApf92eWlPyQhNm1w8z2q5CAkPn1VgGMLD7QNH8teTwbH2aiZOXtrr9YYzw1LATMxz8oXrz7NMy0xx89dPa9NunKM1fo98al7IShOv/7yqIJ0fXTuFORU5tvPXfPvimMNNv3vVJEbGECIzyrJDPgKRxBIvn57i4YpFs2OWFRIHWzMzKslYMJHb7rGxwNh1IQFMKQmLLY/LwZcvHNOrffYnz3xxoa3zK2giIFIICCGYWZ5tETWngxAiJOQ+MHmEZV6885aTplmv/P3gXJqoOxFg+X8stU3q+Pr/XWLpXnU4BHdfM4Xq6uq42zOOq2ERMbhqWjH/3lRnt0qvmTJSu84uNR1XI9puutfJ55dU0ZlEd1BVQTq3L9FeSAwxXF09/ARMV7fgQH0KtQ0pHD7u5XijG3+PA5Bk0k6p5yQTU09S6jvFyMxGipxN5NFEJm2kBdtJCWgOo972TjxNnTgS+GwEHC78Lm9IZPhdHrrdKbSlZEZM99Kt/05tKeBEQSfdLq8mVtz6t8tLwOkCkbgbZHJkoDbFoGFYChgzk0dmRnVBGCw2Odwm4uNxzOPxoph++gKbXBSnwaJSt60vCcC4onQ+saCiX/cXyXUzS3ljRz23Lxkdmja2SOsKmFqauJG387s4k8wqtxeh8fjnF88fgJLYk5+u5fbxRAjq/HQvhRlevnvVpDNSjsr8tJADr5nRBX0LE754rHavLRpr7cr87U0z+e1NM/u0zUgmjcxky92XWrrlPjSzhB11zXz9A+NtuzAV0fiD8NyqfPYd9jA6eJgZzs1cl3qYqsyjlMh6Cv0n8PboPoRd+ofwqJVOj49udwrtKZmcyiiky5NKl1v7dLtTQr+16Sn4XSn4XR6CfbBqFNQV0FDUkHhBxZBkUAkYIcRlwK8BJ/CQlPLegd7nv78c3cd9rvLK15YM+D6yfGFLlsGScQW8/n+X9LlxU4R5+5vLCEoZZc1yOx289+2Lz1KpTp+yXB81P748Km9XfweJixxFleJ2cvc1U2IsrbCjqT3A544+wqWetaRJLaRCoMdFizeHprRcdqZV0JKaTXtKOu0pGXR402n3ptPp8SVl8VAokmXQCBghhBO4H7gEqAXWCCGek1Jui7/mae93IDev0FHipX9IlNZiKBMpXhSDk9HiCNd6/TimXA9VF0LJLFzZFeQ4HPTefjmwVDdWc8OM/on6qhh8DBoBA8wDaqSUewGEEI8B1wADKmAUCoVCkTwCieOONZBVcraLohjmDKZXnhLAHAO+Vp+mUCgUikGCH5cSL4pBgZDxUuueQYQQNwCXSik/o/+/BZgnpfxyxHKfAz4HUFRUNPuxxx7r9b5aW1tJTz83uzTO5brBuV0/VbehS1/qd+GFF66TUvZunP1ZwvzcLS3Knf3Xx54+yyVKjqF03amyavTqvpBSDooPsAB42fT/LuCueOvMnj1b9oXly5f3ab2hwLlcNynP7fqpug1d+lI/YK0cBM/e3n7GjRvX67qeLYbSdafKqtGb+2IwdSGtAcYKISqFEB7go8BzZ7lMCoVCoVAoBiGDxolXShkQQtwBvIw2jPpPUsqtZ7lYCoVCoVAoBiGDRsAASClfAF442+VQKBQKhUIxuBlMXUgKhUKhUCgUSaEEjEKhUCgUiiGHEjAKhUKhUCiGHErAKBQKhUKhGHIoAaNQKBQKhWLIoQSMQqFQKBSKIcegSSXQF4QQDcCBPqyaDxzv5+IMFs7lusG5XT9Vt6FLX+o3SkpZMBCFGUiEEC3AzrNdjiQZStedKqtG0vfFkBYwfUUIsVYOkRwkveVcrhuc2/VTdRu6nOv1MzOU6qrKOjAMlrKqLiSFQqFQKBRDDiVgFAqFQqFQDDmGq4B58GwXYAA5l+sG53b9VN2GLud6/cwMpbqqsg4Mg6Ksw9IHRqFQKBQKxdBmuFpgFAqFQqFQDGGGlYARQlwmhNgphKgRQtx5tsvTF4QQZUKI5UKI7UKIrUKI/6NPzxVCvCqE2K1/5+jThRDiN3qdNwkhZp3dGiRGCOEUQrwvhPi3/r9SCLFar9vjQgiPPt2r/6/R51eczXInQgiRLYR4SgixQz9/C86x8/Y1/ZrcIoT4hxAiZaieOyHEn4QQ9UKILaZpvT5XQohb9eV3CyFuPRt1SUSi52K8cyWEuEufvlMIcWmy2xxkZd0vhNgshNgghFh7tssqhMgT2jO+VQhxX8Q6s/Wy1ujXnBjEZa3Wt7lB/xT2R1ktSCmHxQdwAnuAKsADbAQmne1y9aEexcAs/XcGsAuYBPwMuFOffifwU/33FcCLgADmA6vPdh2SqOPXgb8D/9b/PwF8VP/9APAF/fcXgQf03x8FHj/bZU9Qr78An9F/e4Dsc+W8ASXAPiDVdM5uG6rnDlgMzAK2mKb16lwBucBe/TtH/51ztusWUc+Ez8VY50p/7mwEvEClvh3nQD1rB6Ks+rz9QP4gOq5pwAXA7cB9Eeu8ByzQr7UXgcsHcVmrgTkDef0OJwvMPKBGSrlXStkNPAZcc5bL1GuklHVSyvX67xZgO1rjcQ1aA4n+fa3++xrgEanxLpAthCg+w8VOGiFEKXAl8JD+XwDLgKf0RSLrZtT5KeCi/noj6W+EEJlojeIfAaSU3VLKRs6R86bjAlKFEC7AB9QxRM+dlPIt4GTE5N6eq0uBV6WUJ6WUp4BXgcsGvvS9IpnnYqxzdQ3wmJSyS0q5D6jRtzdQz9qBKOtA0eeySinbpJRvA53mhfVrKlNK+Y7UFMIjhK/BQVXWM8VwEjAlwCHT/1p92pBFN+PNBFYDRVLKOtBEDmCY64ZavX8FfAMI6v/zgEYpZUD/by5/qG76/CZ9+cFIFdAAPCy07rGHhBBpnCPnTUp5GPgFcBBNuDQB6zg3zp1Bb8/VUDiHyZQx1rk60/UeiLICSOAVIcQ6IcTn+qGcp1vWeNusTbDNvjAQZTV4WO8++u5AvKAMJwFjd/CG7BAsIUQ68DTwVSllc7xFbaYNynoLIa4C6qWU68yTbRaVScwbbLjQuiR+L6WcCbShdUPEYijVDd0f5Bo08/xINNPy5TaLDsVzl4hYdRkKdUymjL2t30DVeyDKCnC+lHIW2vX6JSHE4r4XMWE5ervM6SyfLANRVoCbpZRTgUX655Y+lC0uw0nA1AJlpv+lwJGzVJbTQgjhRhMvj0opn9EnHzO6GPTven36UKr3+cDVQoj9aGbMZWgWmWy9WwKs5Q/VTZ+fRbTZf7BQC9RKKVfr/59CEzTnwnkDuBjYJ6VskFL6gWeAhZwb586gt+dqKJzDZMoY61yd6XoPRFmRUhrf9cA/6Z+updMpa7xtlibYZl8YiLIaVlnD1eHvDECX3XASMGuAsUIbFeFBc0R67iyXqdfoZrg/AtullL80zXoOMEY53Ao8a5r+CaExH2gyzOCDDSnlXVLKUillBdr5eUNKeTOwHPiwvlhk3Yw6f1hffrC94QIgpTwKHBJCjNcnXQRs4xw4bzoHgflCCJ9+jRr1G/LnzkRvz9XLwAeEEDm6heoD+rTBRDLPxVjn6jngo/oIlUpgLJqT6UA9a/u9rEKINCFEBoDepfsBYAunz+mU1Rb9mmoRQszX77FPEL4GB1VZhRAuIUS+/tsNXEX/HFcrA+khPNg+aKMFdqF5XH/7bJenj3W4AM10twnYoH+uQOuPfB3YrX/n6ssL4H69zpsZYK/wfqznUsKjkKrQHow1wJOAV5+eov+v0edXne1yJ6jTDGCtfu7+F21kyjlz3oC7gR1oD6q/oo34GJLnDvgHmi+PH+3t89N9OVfAp/Q61gCfPNv1ilHXqOcicA9wdaJzBXxbX28nphExA/Ws7e+y6tfnRv2zdRCVdT+ahaNVv/4m6dPn6PfXHuA+9GC0g62saF3I69CedVuBX6OP+urPj4rEq1AoFAqFYsgxnLqQFAqFQqFQnCMoAaNQKBQKhWLIoQSMQqFQKBSKIYcSMAqFQqFQKIYcSsAoFAqFQqEYcigBoxgS6FlPjaymR4UQh03/Vw3A/m4TQjQIIR7q721H7CdVr0O3ETdBoRhsCCF69Ot0ixDiSSGEr5frt/Zy+T8LIT5sM32OEOI3+u/bhJ4BWQhxuxDiE6bpI3u5PyNz8iahZYu/TwiRbZof9xkjhPhWgvkvCC0bfYUwZThPsmxLhRALTf9DdT0d9LJ0mJ6jG/Q4MP2Cvv2Pmf6Hzl1/4Uq8iEJx9pFSnkCLo4IQ4gdAq5TyFwO828ellHcM1MaFEE4pZQcwQ48+rFAMVjqklMb99yha9uFQIE09sJqQUgZjrN8vSCnXosVSipz+gOnvbWixUnobpfZmKeVavRH/CVqQuCX69hfGXRO+Bfxn5ETTcblC/58duUwSLEWLsbJKL8sDcZfuHXuM8zoAVAAfQ4vCG/PcnQ7KAqMY8hhvd/qbyptCiCeEELuEEPcKIW4WQrwnhNgshBitL1cghHhaCLFG/5yfxD5WCCFmmP6vFEJM0yN5/knfzvtCiGv0+RX6Ouv1z0JTGZcLIf6OFvRMoRhqrADG6Nf4diHE74D1QJkQ4ib9XtsihPipeSUhxH/p98LrQogCfdpn9Xtno35Pmi07F+v30C6h5Ukz7p9/RxZICPEDIcR/6FabOcCjukXhSiHEP03LXSKEeCZyfTNSy8j8DaBcCDFdX894xhQLId4yWaMWCSHuRcvCvkEI8WiM47JfhC2sLiHEX3Rrz1NGnc3L6NaKaqEl7L0d+Jq+/UVGXfXlZggh3tW39U+hRXw2LEo/1Z99u4QQi+Kf0uhjafq/Ra+TUa8/CCG2CiFeEUKk6suMEUK8pp/H9fqz9l5gkV7ur5nPnRAiVwjxv3q53xVCTDPt+096+fcKIb4Sr6xKwCjONaYD/weYipY8bJyUch7wEPBlfZlfA/8tpZwLXK/PS8RDaG92CCHGoUWU3YQW3fMNfVsXAj8XWkjyeuASqSWJ+whgNp3OQ4t2Oel0KqpQnGmElgfncsLiezzwiNQSlPp98uf4AAAFBElEQVSBn6LlMJsBzBVCXKsvlwas1++HN4Hv69OfkVLOlVJOB7ajRTw2qECzgFwJPCCESElUPinlU2hv+TfrloUXgImGYAI+CTycxHZ60KLzToiY9THgZX3b04ENUso70S1UUkt9YjkuUsoDEdsYDzwopZwGNANfjFOO/cADaM+rGVLKFRGLPAJ8U9/WZsLHFcClP/u+GjHdzGgR7j66P1Y5TIwF7pdSTgYa0Z6fAI/q06ej5UCrQ0tWu0Iv939HbOdu4H293N/S62EwAbgU7Tn5faGlIrBFCRjFucYaKWWdlLILLSz2K/r0zWgPRNASD94nhNiAluMjU+j5UOLwJHCVfjN9CvizPv0DwJ36tqrRQm6XA27gD0KIzfq6ZrHynpRyX59rqFCceVL1a3wtWt6rP+rTD0gp39V/zwWqpZbQM4DWqBmZnYPA4/rvv6GlRAGYoltZNgM3A5NN+3xCShmUUu4G9hItJhIitVDzfwU+LrTumwXAi0mubpeBeQ3wSaF1Y0+VWqJCO8zHJZJDUsqV+m/zsegVQogsIFtK+aY+6S+EjzdoCVVBC+lfEWMze3SBMUNK+aUkdrtPSrnBvF392VkipfwngJSyU0rZnmA7F6CdF6SUbwB5en0AnpdSdkkpj6O9CBbF2ojygVGca3SZfgdN/4OEr3cHsED3P0kKKWW7EOJV4BrgRjQzNWgPueullDvNy+sPuGNob2kOoNM0uy3Z/SoUg4SQD4yBEAKs17Jdgx8LI4fNn4FrpZQbhRC3ofl7RC4T63+yPAz8C+0efFIXV3ERQjjRrLjbLQWQ8i0hxGI0q9BfhRA/l1I+YrOJePd4rHoFCBsVElqbksB49vXQu7beXI7Ispifrz1AKr077wZ26xjHIXIfMcuuLDCK4cgrQMg5V5h8WxLwEFpX0BoppZFK/mXgy0J/mgshZurTs4A63anxFsDZHwVXKAYxq4ElQoh8XQDchNZdBFpbY4wq+hjwtv47A6jTLZs3Y+UGIYRD96eoQkvCmAwt+nYBkFIeQXPo/Q5hy2lM9LL8BM1Ssili3iigXkr5BzQr1Cx9lj9eV0cE5UKIBfrvmwgfi/3AbP339ablLfUxkFI2AadM/i23ED7ep8N+9HoJIWYBlfEWllI2A7VGd6HQMn77YpVb5y308y2EWAoc17fTK5SAUQxHvgLM0R3ItqE5ySVESrkOrc/a3If+Q7Tuok1CGx75Q33674BbhRDvAuNQVhfFOY6Usg64C1iO5j+yXkr5rD67DZgshFiH5iNzjz79u2jC51W0TOZmdqI1yC8Ct0spO0mOP6P5zGwwnEzRurMOSSm3xVnvUSHEJrQRTGlo1tZIlgIbhBDvo4mMX+vTH0R7BjyaRPm2oz0bNgG5wO/16XcDvxZCrECzPBj8C/iQ4cQbsa1b0fzuNqH5Hd3D6fM0kKt3GX4BLUt1Im4BvqKXYxUwAi0TdUB37P1axPI/QH8Gozn73tqXgqps1AqFDbo5e455GLXQYktUAxP6e7io0IZRz9H7fRUKRT8itHgx70sp/5hwYcWQQVlgFAp7OoDLhR7ITmiBo1ajjR7qN/Ei9EB2aFacAY2hoVAMR3SrzzQ0h1nFOYSywCgUCoVCoRhyKAuMQqFQKBSKIYcSMAqFQqFQKIYcSsAoFAqFQqEYcigBo1AoFAqFYsihBIxCoVAoFIohhxIwCoVCoVAohhz/H78xKLczBMubAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 648x324 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig,axes = plt.subplots(nrows=1, ncols=2, figsize=(9, 4.5), sharey=True, gridspec_kw = {'width_ratios':[2, 1]})\n",
"ax = axes[0]\n",
"ax.plot(time[0:1000], observed[0:1000])\n",
"ax.grid()\n",
"ax.set_title('Time Series of Observed Streamflow')\n",
"ax.set_xlabel('Time [year]')\n",
"ax.set_ylabel('Observed Streamflow')\n",
"\n",
"ax = axes[1]\n",
"ax.plot(norm_prob, sflo, label='Naive Log-Normal')\n",
"ax.plot(gev_prob, sflo, label='GEV')\n",
"ax.hist(observed, orientation=\"horizontal\", density=True, bins=30, alpha=0.5, label='Observed')\n",
"ax.grid()\n",
"ax.set_xlabel('Probability Distribution Function')\n",
"ax.set_title('Histogram and GEV Fit')\n",
"ax.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can see from the plot here that a GEV distribution fits far better than a log-normal distribution here, and might otherwise be inclined to model the data as GEV.\n",
"However, this isn't necessarily necessary and may in fact lead to confusion.\n",
"\n",
"While it may not seem like this is a huge deal, inappropriate use of the GEV can lead to poor estimates, particularly as we extrapolate, estimating the probability of floods for which we have fewer and fewer observations.\n",
"For example, let's calculate the 99th and 99.9th percentiles of the estimated GEV distribution, as well as the same percentiles of the observed distribution (we have 100000 samples, which should give us reasonable estimates).\n",
"These percentiles correspond to the 100- and 1000-year floods."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"q99 of GEV:242.85588789724653\n",
"q999 of GEV: 474.2608694127031\n",
"q99 of observations: 212.594974711907\n",
"q99.9 of observations: 316.5251547839271\n",
"ratio of q99 in GEV to observed: 1.1423406796249387\n",
"ratio of q999 in GEV to observed: 1.4983354790125694\n"
]
}
],
"source": [
"q99_gev = gev.isf(0.01, shapehat, scale=scalehat, loc=lochat)\n",
"q999_gev = gev.isf(0.001, shapehat, scale=scalehat, loc=lochat)\n",
"print(f\"q99 of GEV:{q99_gev}\")\n",
"print(f\"q999 of GEV: {q999_gev}\")\n",
"q99_obs = np.percentile(observed, 100*0.99)\n",
"q999_obs = np.percentile(observed, 100*0.999)\n",
"print(f\"q99 of observations: {q99_obs}\")\n",
"print(f\"q99.9 of observations: {q999_obs}\")\n",
"print(f\"ratio of q99 in GEV to observed: {q99_gev / q99_obs}\")\n",
"print(f\"ratio of q999 in GEV to observed: {q999_gev / q999_obs}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In other words, we can see that our GEV model is substantially over-estimating the 100- and 1000-year floods, _in this simple setup_.\n",
"\n",
"Of course, the real question is whether *real-world* data is better represented by a GEV distribution or a condional log-normal distribution.\n",
"The answer is probably that it depends on the specific data, and the purpose of this post is just to illustrate how using conditional information can be an alternative to using GEV models.\n",
"As always, when working with real data be sure to *check your model* when you make any assumption, including those suggested here!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.1"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment