{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to predict catastrophic events ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At the end of that short but preciser tutorial, you'll understand what is a GEV (Generalized extreme value distribution) and how to fit it on your data in python. \n",
    "\n",
    "What we need to think about before speaking of extreme values statistics, is that in general statistics focus on average values and observable phenomenoms, for which we **do have datas**. Extreme values statistics is different from \"average\" statistics in this regard, indeed as we will see, we here try to model something for which **we have no data!** You will see soon what we mean by that. Before getting there I want to show you why extreme values modelling is an actual problem, of high interest for insurance companies and public services. Let's start with two examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <center> Montpellier, France, 2014\n",
    "**3 hours of rain = 50% of total annual average.** *Montpellier is not known to be a city with a high annual rain level, located in the dry location of Provence, therefore if we look at the average value we have a very common one. But it is not excluded that sometimes, rarely, very high values appear, as we can see:*\n",
    "\n",
    "<img align=\"center\" src=\"https://s1.lemde.fr/image2x/2014/10/07/0x0/4501632_5_f1d0_a-picture-taken-on-october-7-2014-shows-a_98f1f199a7769f5add31267568c48d85.jpg\" style=\"width: 700px;\"/>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <center> Saltina river in Brigue, Switzerland, 1993\n",
    "**Unprecedented measured flows and precipitations in the twentieth century**.\n",
    "*Same problem here, this rivers has reasonable flows in average, but sometimes the flow can be exceptional.*\n",
    "<img align=\"center\" src=\"http://www.planat.ch/uploads/pics/p0032_01.jpg\" style=\"width: 600px;\"/>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <center> Number of natural catastrophes since 1980 is growing \n",
    "**Those \"extreme\" events tends to become more and numerous those days, and make extreme value analysis \"popular\".**\n",
    "<img align=\"center\" src=\"https://cdn.static-economist.com/sites/default/files/images/2017/08/blogs/graphic-detail/20170902_woc210.png\" style=\"width: 600px;\"/>\n",
    "<br>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <center> Some could wonder what statistical model shall we use to predict a Nadal loss at Roland Garros \n",
    "\n",
    "<img align=\"center\" src=\"https://www.freebetreview.com.au/Images/blog/rafael-nadal-french-open-record-infographic.png\" style=\"width: 600px;\"/>\n",
    "<br>\n",
    "<center> (Almost) just kidding"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Question: What if we would like to quantify the probability of such events ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we said the particularity of the extreme value analysis is not to model an average phenomenom but an extreme phenomenom."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to understand what we mean by that, let's work with the Nidd river dataset here : https://www.stat.auckland.ac.nz/~wild/data/Rdatasets/csv/evir/nidd.thresh.csv, it contains $n=154$ trimestrial measures of the river's flows (in $m³/s$), over a period of 38.5 years (source: Natural Environment Research, 1975). \n",
    "\n",
    "Let's plot a histogram of the data.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAH59JREFUeJzt3Xm8HHWZ7/HPlyQQIEBYIwQ0AeKCRhnIIIviCaCIiAQvzEURgYsT2UFxrlEZFhFFr6A4OiIjXEAZAkQEFBe4kQPCIJIIGJLAECEQMLJJyIKAIc/9o35NiqZPd52lunNOfd+vV7+69np+VdX9VP1qU0RgZmbVtVanAzAzs85yIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4J4JekjRXUlen4+gkSQdJWiRpuaR/6MP4b0zjDuuh/5mSftxk/IWS9unF/C6U9K+9jbNMko6V9GRaDptKCknbdzouqyYngpxGfzCSjpR0e609It4eEd0tpjMu/bCHlxRqp30TOCEiRkXEPfU9U9nnSFor1+0rki4FiIjH0rivtCPYiDgmIs4uY9qSunu7YyBpBHA+8IG0HJ4tI7aqkPQpSQtSUv2VpK1y/X6Zutc+L0uaM0Dz/aKkKyXd0eB/4yZJHxiI+bSDE8EgtAYkmDcBc1sMsxVwaBtiaaqno44BmK7yia6XxgAjab0MrQVJ7wO+ChwIbAI8AlxZ6x8R+6VkOyoiRgH/BVwzEPOOiK9GxMeAC4C9cjGtD+wM3DoQ82kHJ4Jeyh81SNpF0ixJS9Nh/vlpsNvS95K0F7KbpLUknSbpUUlPSbpc0ka56X4y9XtW0r/WzedMSTMk/VjSUuDINO87JS2RtFjSdyWtnZteSDpO0kOSlkk6W9J2aZylkq7OD19XxoaxSlpH0nJgGHCfpD81WVTfAM5qlLTqj5gkjZd0a4rzZmCzuuEPzy2bL7VYP5dK+r6kX0haAUxO3b6S+s+X9OHc8MMlPSNpp9S+q6T/Ssv1vvzeftr7P0fSHcALwLZ18+5pe8gP82bgwdS6RNJvGgyzUVrmT6dyn1ZLOql959T8ibQcd0jtn5J0XdFYelh+Z0q6Jm1ry5Qd2b1Z0hfStrBIuT3dFOvFaRt8QtmR37DUbztJv0nr7RlJV0ganRt3oaTPSfqjpOclXSVpZJE4cw4AromIuRHxMnA2sKek7RqUbRzwXuBHPZS9tl0elcr5nKRjJP1jinGJpO/WjXMJcA5wfa7z3sAdEfFSX9dD20WEP+kDLAT2qet2JHB7o2GAO4HDU/MoYNfUPA4IYHhuvP8FLCD78xgFXAv8KPXbAVgOvAdYm6zq5e+5+ZyZ2qeQJe91yfY4dgWGp/nNB07JzS+AG4ANgbcDLwEz0/w3AuYBR/SwHHqMNTft7ZssxwAmALOBT6VuXwEubbR80nI8H1gH2BNYBvy4btnsmfqfD6ysX0+5eV8KPA/skZbVyNTtK6n/6cAVueH3Bx5IzWOBZ4EPpXHfn9o3T/27gcfS8hwOjKibd8PtoUGMjbaPV5cpcDnZH8sGadj/Bo7O9Ts1NV8E/Ak4NtfvM72JpUFsZwIvAvumMl5Otpf9JWAE8M/AI7nhrwN+AKwPbAH8Hvh06rd9WobrAJuT7SB9u+639Huyo8dNyLbhY1K/NwJLmnw+noY7D/j33DTHpmV5YIOynQ50Nyl7bb1cmLabD6RlcV0q21jgKeB9dePtBvws135hbhn0aT20+9PxANakT9owl9dtcC/QcyK4DTgL2KyHDSr/Q58JHJdrfwvZn/vwtIFemeu3HvAyr00Et7WI/RTgp7n2APbItc8GPp9rPy//o6ybVo+x5qbdKhFsT/aH+lj6I2iYCNIPfiWwfm78/2R1IjgdmJ7rt35+2TSY96XA5Q261RLB9mSJZr3UfgVwemr+PLmEl7r9mpQwyRLBl5uUu+H20GC4RttHbZkNI0vaO+T6fZr0BwYcDdyQmucDn6otH+BRYKfexNIgtjOBm3PtB5D9Joal9g1SrKPJqrheAtbNDf8x4JYepj0FuKfut/SJXPs3gAt7Ge/ewDPAO8l2kH4ArAI+1mDYBcCRBdbL2Fy3Z4H/mWv/CWmHC/guWTXTTOCDuWEeBbbpz3po98dVQ683JSJG1z7AcU2GPRp4M/CApLvzVQ4NbEW2gdQ8SvZHOCb1W1TrEREvkG2AeYvyLelw/eeS/qKsuuir1FWpAE/mmv/WoH1UH2ItLCJ+QZYIpjYZbCvguYhYUTe/fP/8slnB65dNvUU99YiIBWR/oAdIWg/4CFnigezcxyGpCmCJpCVkR2lbFpk2vdseerIZ2VFh/fIfm5pvBd4r6Q1kSeMqYI9U7bERcO8AxFK/nTwTq0/s/y19jyJbXiOAxbnl9QOyvWckbSFpeqoyWgr8mNdvo3/JNb9Az9tkQxExEziD7A/6UbLksgx4PD+cpPcAbwBmFJhsod9NRJwQEYdExN4R8as0n4nA0oiobScDsU2UzomgHyLiochOFm0BfB2YoexEUaNHuv6Z7IdTU9sTfhJYDGxd6yFpXWDT+tnVtX8feACYEBEbAl8E1PfSFI61t04jq1ZYr4f+i4GN03LLzy/ff5taS/rzrl829Vo9UvdKsj3XA4F5KTlA9if/o/yOQESsHxHnFpl2k+2hN54hO/qqX/5PpHksIPvDPInsKHEZ2Z/pVLIj11UDGEsri8iOCDbLLa8NI+Ltqf/XyJbXO9M2+gkKbqNafYlxT5/DasNGxPciYkJEbEGWEIYD99dN8gjg2ohY3r8it/Qh4MZcbO1YD/3mRNAP6WTd5unHtyR1fgV4muzwNH8y8UrgM8pOjI4i24O/KiJWku2lHCBpd2UncM+i9Q9mA2ApsFzSW4FjB6xgzWPtlcgutZ1D9kNs1P9RYBbZieW1057bAblBZgAflvSetGy+TP+32+lk9b/HsvpoALI91gMk7StpmKSRkrokbd1wKnWabA+FpT3vq4FzJG0g6U3AZ1NsNbcCJ7D6qpTuuvamsaSTtEf2Jq4eYl0M3AScJ2lDZRcZbKfsSh7IttHlZCfFxwL/0otp1y4x7ulzRSrLSEnvUOaNZOdNLoiI52rTSjtWh5BVEZZtf+AXuXn3e5toByeC/vkgMFfZlTQXAIdGxIupaucc4I50yLwrcAnZ1Qq3kZ18exE4ESAi5qbm6WR7wMvITkq91GTenwM+nob9D7IqgoHSY6x9dBrZycCefBx4N/BXssP8y2s90rI5nuwPezHwHHWH/b2V/sDuBHYnt9zS4fyBZEdXT5Pt8f4LxX8nDbeHPoR4IrACeBi4nazsl+T630r2J3tbD+09xpKS6abA7/oQVyOfJKvKmke2bmawuirtLGAnspP3N5JddDDQRpItn+VkJ57vBOpvHpySYrilhPm/StlVgG8ju0S1ZqC2iVIpndCwNUjaC19CVu3zSKfjsaEjHXEdn6orbABJ+ifg4Ij4p07H0ltOBGsISQeQXX0gsit63k12BYhXkNkgoOz+imURcWenY+mtTt+haqsdSFYdI7I680OdBMwGj4i4qdMx9JWPCMzMKs4ni83MKm5QVA1tttlmMW7cuE6H0RYrVqxg/fXXuMuM26LKZYdql99lL6fss2fPfiYiNm813KBIBOPGjWPWrFmdDqMturu76erq6nQYHVHlskO1y++yd5UybUmPth7KVUNmZpXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFTco7izuj3HTbmw9UAkWnrt/R+ZrZtZbPiIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIjAzqzgnAjOziis1EUj6jKS5ku6XdKWkkZLGS7pL0kOSrpK0dpkxmJlZc6UlAkljgZOASRHxDmAYcCjwdeBbETEBeA44uqwYzMystbKrhoYD60oaDqwHLAb2Amak/pcBU0qOwczMmlBElDdx6WTgHOBvwE3AycDvImL71H8b4JfpiKF+3KnAVIAxY8bsPH369D7FMOeJ5/sWfD9NHLtRn8Zbvnw5o0aNGuBoBocqlx2qXX6XvZyyT548eXZETGo13PBS5g5I2hg4EBgPLAGuAfZrMGjDTBQRFwEXAUyaNCm6urr6FMeR027s03j9tfCwrj6N193dTV/LOthVuexQ7fK77F0djaHMqqF9gEci4umI+DtwLbA7MDpVFQFsDfy5xBjMzKyFMhPBY8CuktaTJGBvYB5wC3BwGuYI4PoSYzAzsxZKSwQRcRfZSeE/AHPSvC4CPg98VtICYFPg4rJiMDOz1ko7RwAQEWcAZ9R1fhjYpcz5mplZcb6z2Mys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKziWiYCSetLWis1v1nSRySNKD80MzNrhyJHBLcBI9Mbx2YCRwGXlhmUmZm1T5FEoIh4Afgo8G8RcRCwQ7lhmZlZuxRKBJJ2Aw4Dam95KfVhdWZm1j5FEsEpwBeAn0bEXEnbkr1TwMzMhoCWe/YRcStwq6T1U/vDwEllB2ZmZu1R5Kqh3STNA+an9ndJ+vfSIzMzs7YoUjX0bWBf4FmAiLgP2LPMoMzMrH0K3VAWEYvqOr1SQixmZtYBRa7+WSRpdyAkrU12fmB+uWGZmVm7FDkiOAY4HhgLPA7smNrNzGwIaHpEIGkYcHhEHNameMzMrM2aHhFExCvAgW2KxczMOqDIOYI7JH0XuApYUesYEX8oLSozM2ubIolg9/T95Vy3APYa+HDMzKzditxZPLkdgZiZWWe0TASSTm/UPSK+3Ki7mZkNLkWqhlbkmkcCH8b3EZiZDRlFqobOy7dL+iZwQ2kRmZlZW/XlncXrAdsOdCBmZtYZRc4RzCG7SghgGLA5r72CyMzMBrEi5wg+nGteCTwZEStLisfMzNqsSNXQcOAvEfEoMAE4TtLocsMyM7N2KZIIfgK8Iml74GJgPPCfpUZlZmZtUyQRrEpVQR8Fvh0RnwG2LDcsMzNrlyKJ4O+SPgZ8Evh56jaivJDMzKydiiSCo4DdgHMi4hFJ44EfF5m4pNGSZkh6QNL89P7jTSTdLOmh9L1xfwpgZmb90zIRRMS8iDgpIq5M7Y9ExLkFp38B8KuIeCvwLrI7kqcBMyNiAjAztZuZWYe0TASSJqS9+nmSHq59Coy3IdlL7i8GiIiXI2IJ2fsNLkuDXQZM6Xv4ZmbWX4qI5gNItwNnAN8CDiCrKlJEnNFivB2Bi4B5ZEcDs4GTgSciYnRuuOci4nXVQ5KmAlMBxowZs/P06dN7UazV5jzxfJ/G66+JYzfq03jLly9n1KhRAxzN4FDlskO1y++yl1P2yZMnz46ISa2GK5IIZkfEzpLmRMTE1O23EfHeFuNNAn4H7BERd0m6AFgKnFgkEeRNmjQpZs2a1aosDY2bdmOfxuuvhefu36fxuru76erqGthgBokqlx2qXX6XvauUaaf/75aJoMjJ4hclrQU8JOkESQcBWxQY73Hg8Yi4K7XPAHYCnpS0ZQpyS+CpAtMyM7OSFEkEp5A9aO4kYGfgE8ARrUaKiL8AiyS9JXXam6ya6Ibc+EcA1/cyZjMzG0BFHkN9N4CkiIijejn9E4ErJK0NPEx2fmEt4GpJRwOPAYf0cppmZjaAijx9dDeyK39GAW+U9C7g0xFxXKtxI+JeoFH91N69DdTMzMpRpGro28C+wLMAEXEf2WWhZmY2BBR6MU1ELKrr9EoJsZiZWQcUeR/BIkm7A5Hq+k/C7yw2MxsyihwRHAMcD4wluyR0x9RuZmZDQJGrhp4BDmtDLGZm1gFFrhoaT3YZ6Lj88BHxkfLCMjOzdilyjuA6sstHfwasKjccMzNrtyKJ4MWI+E7pkZiZWUcUSQQXSDoDuAl4qdYxIv5QWlRmZtY2RRLBROBwYC9WVw1Fajczs0GuSCI4CNg2Il4uOxgzM2u/IvcR3AeMbjmUmZkNSkWOCMYAD0i6m9eeI/Dlo2ZmQ0CRRND0lZRmZja4Fbmz+NZ2BGJmZp1R6OmjZmY2dDkRmJlVXI+JQNLM9P319oVjZmbt1uwcwZaS3gd8RNJ0QPmevrPYzGxoaJYITgemAVsD59f1853FZmZDRI+JICJmADMk/WtEnN3GmMzMrI2KXD56tqSPsPqF9d0R8fNywzIzs3ZpedWQpK8BJwPz0ufk1M3MzIaAIncW7w/sGBGrACRdBtwDfKHMwMzMrD2K3keQf+jcRmUEYmZmnVHkiOBrwD2SbiG7hHRPfDRgZjZkFDlZfKWkbuAfyRLB5yPiL2UHZmZm7VHkiICIWAzcUHIsZmbWAX7WkJlZxTkRmJlVXNNEIGktSfe3KxgzM2u/pokg3Ttwn6Q3tikeMzNrsyIni7cE5kr6PbCi1tHvLDYzGxqKJIKzSo/CzMw6puXJ4vTO4oXAiNR8N1D4XQSShkm6R9LPU/t4SXdJekjSVZLW7mPsZmY2AIo8dO6fgRnAD1KnscB1vZjHycD8XPvXgW9FxATgOeDoXkzLzMwGWJHLR48H9gCWAkTEQ8AWRSYuaWuyh9b9MLWL7IU2M9IglwFTeheymZkNpCKJ4KWIeLnWImk42RvKivg28L+BVal9U2BJRKxM7Y+THWGYmVmHFDlZfKukLwLrSno/cBzws1YjSfow8FREzJbUVevcYNCGSUXSVGAqwJgxY+ju7i4Q6uudOnFl64FK0Nd4ly9f3udxB7sqlx2qXX6XvbujMSii+c69pLXI6vE/QPZH/mvgh9FixPTymsOBlcBIYEPgp8C+wBsiYqWk3YAzI2LfZtOaNGlSzJo1q1iJ6oybdmOfxuuvhefu36fxuru76erqGthgBokqlx2qXX6XvauUaUuaHRGTWg1X5Omjq9LLaO4i23t/sFUSSON9gfS46nRE8LmIOEzSNcDBwHTgCOD6VtMyM7PyFLlqaH/gT8B3gO8CCyTt1495fh74rKQFZOcMLu7HtMzMrJ+KnCM4D5gcEQsAJG0H3Aj8suhMIqIb6E7NDwO79DZQMzMrR5Grhp6qJYHkYeCpkuIxM7M26/GIQNJHU+NcSb8AriY7R3AI2d3FZmY2BDSrGjog1/wk8L7U/DSwcWkRmZlZW/WYCCLiqHYGYmZmndHyZLGk8cCJwLj88H4MtZnZ0FDkqqHryC7x/BmrHxVhZmZDRJFE8GJEfKf0SMzMrCOKJIILJJ0B3AS8VOsYEYXfSWBmZmuuIolgItkzg/ZiddVQpHYzMxvkiiSCg4Bt84+iNjOzoaPIncX3AaPLDsTMzDqjyBHBGOABSXfz2nMEvnzUzGwIKJIIzig9iiGor+9BOHXiSo7s5zsU+vouBDOrpiLvI7i1HYGYmVlnFLmzeBmrXye5NjACWBERG5YZmJmZtUeRI4IN8u2SpuD3CZiZDRlFrhp6jYi4Dt9DYGY2ZBSpGvpornUtYBKrq4rMzGyQK3LVUP69BCuBhcCBpURjZmZtV+Qcgd9LYGY2hDV7VeXpTcaLiDi7hHjMzKzNmh0RrGjQbX3gaGBTwInAzGwIaPaqyvNqzZI2AE4GjgKmA+f1NJ6ZmQ0uTc8RSNoE+CxwGHAZsFNEPNeOwMzMrD2anSP4P8BHgYuAiRGxvG1RmZlZ2zS7oexUYCvgNODPkpamzzJJS9sTnpmZla3ZOYJe33VsZmaDj//szcwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzonAzKziSksEkraRdIuk+ZLmSjo5dd9E0s2SHkrfG5cVg5mZtVbmEcFK4NSIeBuwK3C8pB2AacDMiJgAzEztZmbWIaUlgohYHBF/SM3LgPnAWLLXXF6WBrsMmFJWDGZm1poiyn8PvaRxwG3AO4DHImJ0rt9zEfG66iFJU4GpAGPGjNl5+vTpfZr3nCee79N4nTJmXXjyb/2bxsSxGw1MMG22fPlyRo0a1ekwOqbK5XfZyyn75MmTZ0fEpFbDlZ4IJI0CbgXOiYhrJS0pkgjyJk2aFLNmzerT/MdNu7FP43XKqRNXct6clq+SbmrhufsPUDTt1d3dTVdXV6fD6Jgql99l7ypl2pIKJYJSrxqSNAL4CXBFRFybOj8pacvUf0vgqTJjMDOz5sq8akjAxcD8iDg/1+sG4IjUfARwfVkxmJlZa/2rg2huD+BwYI6ke1O3LwLnAldLOhp4DDikxBjMzKyF0hJBRNwOqIfee5c1XzMz6x3fWWxmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxTgRmZhXnRGBmVnFOBGZmFVfay+utc8ZNu7HTIfTJqRNXcmQfY1947v4DHI1ZdfiIwMys4pwIzMwqzonAzKzinAjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwqzjeU2ZDQqZvofCObDQU+IjAzqzgnAjOzinMiMDOrOJ8jMBukqnhepIplbgcfEZiZVVxHEoGkD0p6UNICSdM6EYOZmWXaXjUkaRjwPeD9wOPA3ZJuiIh57Y7FrL8GsqqiP+9jaKcyqmfW9LKXWSXVrOztqpLqxBHBLsCCiHg4Il4GpgMHdiAOMzMDFBHtnaF0MPDBiPhUaj8ceHdEnFA33FRgamp9C/BgWwPtnM2AZzodRIdUuexQ7fK77OV4U0Rs3mqgTlw1pAbdXpeNIuIi4KLyw1mzSJoVEZM6HUcnVLnsUO3yu+ydLXsnqoYeB7bJtW8N/LkDcZiZGZ1JBHcDEySNl7Q2cChwQwfiMDMzOlA1FBErJZ0A/BoYBlwSEXPbHccarHLVYTlVLjtUu/wuewe1/WSxmZmtWXxnsZlZxTkRmJlVnBNBm0m6RNJTku7PddtE0s2SHkrfG6fukvSd9CiOP0raqXOR918PZT9T0hOS7k2fD+X6fSGV/UFJ+3Ym6oEhaRtJt0iaL2mupJNT9yG/7puUfcive0kjJf1e0n2p7Gel7uMl3ZXW+1XpwhkkrZPaF6T+49oSaET408YPsCewE3B/rts3gGmpeRrw9dT8IeCXZPde7Arc1en4Syj7mcDnGgy7A3AfsA4wHvgTMKzTZehH2bcEdkrNGwD/nco45Nd9k7IP+XWf1t+o1DwCuCutz6uBQ1P3C4FjU/NxwIWp+VDgqnbE6SOCNouI24C/1nU+ELgsNV8GTMl1vzwyvwNGS9qyPZEOvB7K3pMDgekR8VJEPAIsIHs8yaAUEYsj4g+peRkwHxhLBdZ9k7L3ZMis+7T+lqfWEekTwF7AjNS9fr3XtocZwN6SGt2EO6CcCNYMYyJiMWQ/GmCL1H0ssCg33OM0/wENViek6o9LalUjDOGyp8P9fyDbO6zUuq8rO1Rg3UsaJule4CngZrIjnCURsTINki/fq2VP/Z8HNi07RieCNVuhx3EMct8HtgN2BBYD56XuQ7LskkYBPwFOiYilzQZt0G1Ql79B2Sux7iPilYjYkewpCrsAb2s0WPruSNmdCNYMT9YO+9P3U6n7kH8cR0Q8mX4oq4D/YHUVwJAru6QRZH+EV0TEtalzJdZ9o7JXad0DRMQSoJvsHMFoSbUbevPle7Xsqf9GFK9O7TMngjXDDcARqfkI4Ppc90+mK0h2BZ6vVSMMFXX13gcBtSuKbgAOTVdRjAcmAL9vd3wDJdXzXgzMj4jzc72G/LrvqexVWPeSNpc0OjWvC+xDdo7kFuDgNFj9eq9tDwcDv4l05rhUnT6rXrUPcCXZYfDfybL/0WR1gDOBh9L3JrH6ioPvkdUpzgEmdTr+Esr+o1S2P5L9CLbMDf+lVPYHgf06HX8/y/4eskP8PwL3ps+HqrDum5R9yK974J3APamM9wOnp+7bkiW3BcA1wDqp+8jUviD137YdcfoRE2ZmFeeqITOzinMiMDOrOCcCM7OKcyIwM6s4JwIzs4pzIrAhQ9IruSdZ3itpnKQuST8fwHlsJWlG6yF7Pd190tNGr8tfXy/pV5IG7eMVbHBo+6sqzUr0t8hu5X/VQD/GNyL+zOobgQZyuv9PUhewFbASXr0BaZOIeGKg52eW5yMCq4z07P/r0kPOfifpnan7HEmj0128z0r6ZOr+I0n71E1jnNL7FCQdKenatNf+kKRv9DDfhZK+KulOSbMk7STp15L+JOmY2nARcRpwBbB96tRF9kgCJJ0raV6K/ZsDu2Ss6nxEYEPJuukpjwCPRMRBdf3PAu6JiCmS9gIuJ3vg2R3AHsCjwMPAe1O/XYFjW8xzR7Knab4EPCjp3yJiUYPhFkXEbpK+BVya5jcSmAtcKOnTwLuAN5PdZQuwH3CdpE3IHsHw1oiI2iMLzAaKE4ENJa+rGqrzHuB/AETEbyRtKmkj4LdkL815lOyJmFNTvfxfY/Wz5HsyMyKeB5A0D3gTr32Eck3tz30O2YtKlgHLJL0oaXRE/KDBOHsAnwNWAS8CP5R0IzBg5zzMwFVDVi09PeL3NrKjgPeSVcU8TXYe4LcFpvlSrvkVet65qg23qm6cVY3GkbQt2VHEy5E9l34Xsqd3TgF+VSAus8KcCKxKbgMOA0gnZp+JiKWpKmczYEJEPAzcTrYnXiQRlGU/0h9+eo7/RhHxC+AUsuooswHjqiGrkjOB/yvpj8ALrH7cL2RvzBqWmn8LfI0sIXTKB4ETU/MGwPWSRpId1XymY1HZkOSnj5qtYSStA9wREZM6HYtVgxOBmVnF+RyBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlV3P8HFmCh6woIA3kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "df = pd.read_csv('data/nidd.thresh.csv')\n",
    "df = df['x']\n",
    "df.hist();\n",
    "plt.xlabel('Flow in m³/s');\n",
    "plt.ylabel('Number of measures');\n",
    "plt.title('Histogram of Nidd river\\'s flows, mean=%d m³/s' % df.mean());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that most flows are under $80$ $m³/s$, a significant part is under $100$ $m³/s$, but we also see that even if the flows are in general low, it can also become $300$ $m³/s$ high.\n",
    "\n",
    "This is typically a problem where we cannot use the average value to represent anything, since most flows are under $80$ $m³/s$, and some are very high."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Warm-up questions: what is the probability to observe a flow superior to $160$ $m³/s$ ? Then $255$ $m³/s$ ?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 210,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADL9JREFUeJzt3W+IZfddx/H3x938aRt182cS1iQ4G1m0QbQJQ4hW+iApaBJxI6SwIrpIYEGtplaxWwVbnyWibRVKy9pEVglN6jaQYKoS0gTxgVsnf5p/a8yaxnSbNTvFJm19oI39+uCebYbtzM7dmTtz3e+8XzDcc849l/v7cZb3nHvu3LupKiRJZ77vmfYAJEmTYdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDWxdSOf7KKLLqrZ2dmNfEpJOuM99thjX62qmZX229Cgz87OMj8/v5FPKUlnvCT/Ps5+XnKRpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJjb0k6JrMbvvwak870u33zSV55Wk0+UZuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJamKsoCf5rSTPJnkmyaeTnJtkR5JDSV5Icm+Ss9d7sJKk5a0Y9CSXAr8JzFXVjwJbgN3AHcBHq2on8DXg1vUcqCTp1Ma95LIVeEuSrcBbgWPAdcDB4f4DwM2TH54kaVwrBr2qvgL8MfAyo5C/DjwGvFZVbwy7HQUuXerxSfYmmU8yv7CwMJlRS5K+yziXXM4HdgE7gB8A3gbcsMSutdTjq2p/Vc1V1dzMzMxaxipJOoVxLrm8G/hSVS1U1beA+4CfBLYNl2AALgNeWacxSpLGME7QXwauTfLWJAGuB54DHgFuGfbZA9y/PkOUJI1jnGvohxi9+fk48PTwmP3AB4D3JzkCXAjcuY7jlCStYOvKu0BVfQj40EmbXwSumfiIJEmr4idFJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNjBX0JNuSHEzyL0kOJ/mJJBckeSjJC8Pt+es9WEnS8sY9Q/9T4O+q6keAHwcOA/uAh6tqJ/DwsC5JmpIVg57k+4B3AXcCVNX/VNVrwC7gwLDbAeDm9RqkJGll45yhXwEsAH+R5Ikkn0ryNuCSqjoGMNxevI7jlCStYJygbwWuBj5RVVcB/8VpXF5JsjfJfJL5hYWFVQ5TkrSScYJ+FDhaVYeG9YOMAv9qku0Aw+3xpR5cVfuraq6q5mZmZiYxZknSElYMelX9B/DlJD88bLoeeA54ANgzbNsD3L8uI5QkjWXrmPv9BnB3krOBF4FfYfTL4DNJbgVeBt6zPkOUJI1jrKBX1ZPA3BJ3XT/Z4UiSVstPikpSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktTE2EFPsiXJE0n+ZljfkeRQkheS3Jvk7PUbpiRpJadzhn4bcHjR+h3AR6tqJ/A14NZJDkySdHrGCnqSy4CbgE8N6wGuAw4OuxwAbl6PAUqSxjPuGfrHgN8Fvj2sXwi8VlVvDOtHgUsnPDZJ0mlYMehJfhY4XlWPLd68xK61zOP3JplPMr+wsLDKYUqSVjLOGfo7gZ9L8hJwD6NLLR8DtiXZOuxzGfDKUg+uqv1VNVdVczMzMxMYsiRpKSsGvao+WFWXVdUssBv4fFX9IvAIcMuw2x7g/nUbpSRpRWv5O/QPAO9PcoTRNfU7JzMkSdJqbF15lzdV1aPAo8Pyi8A1kx+SJGk1/KSoJDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWritL6cazOa3ffg1J77pdtvmtpzSzrzeIYuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJlYMepLLkzyS5HCSZ5PcNmy/IMlDSV4Ybs9f/+FKkpYzzhn6G8BvV9XbgWuBX09yJbAPeLiqdgIPD+uSpClZMehVdayqHh+WvwEcBi4FdgEHht0OADev1yAlSSs7rWvoSWaBq4BDwCVVdQxG0QcuXuYxe5PMJ5lfWFhY22glScsaO+hJzgM+C7yvqr4+7uOqan9VzVXV3MzMzGrGKEkaw1hBT3IWo5jfXVX3DZtfTbJ9uH87cHx9hihJGsc4f+US4E7gcFV9ZNFdDwB7huU9wP2TH54kaVxbx9jnncAvAU8neXLY9nvA7cBnktwKvAy8Z32GKEkax4pBr6p/BLLM3ddPdjiSpNXyk6KS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpiRX/k2hNz+y+B6c9hA330u03TXsI0hnLM3RJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU34wSL9/5JM53mrpvO80gR5hi5JTRh0SWrCoEtSE15Dl6ZtM75vsBnnvAHWdIae5GeSPJ/kSJJ9kxqUJOn0rTroSbYAHwduAK4EfiHJlZMamCTp9Kzlkss1wJGqehEgyT3ALuC5SQxM2lDTugQwTc5542zQpZ61XHK5FPjyovWjwzZJ0hSs5Qx9qV913/VrKMleYO+w+s0kz6/hOc8kFwFfnfYgpmTVc29yzuix35yWn/vaXxn84Dg7rSXoR4HLF61fBrxy8k5VtR/Yv4bnOSMlma+quWmPYxo289xhc8/fuU937mu55PLPwM4kO5KcDewGHpjMsCRJp2vVZ+hV9UaS9wJ/D2wB7qqqZyc2MknSaVnTB4uq6nPA5yY0lm423WWmRTbz3GFzz9+5T1Gq+SenJGmz8LtcJKkJg75KSe5KcjzJM4u2XZDkoSQvDLfnD9uT5M+Gr0h4KsnV0xv52i0z9w8n+UqSJ4efGxfd98Fh7s8n+enpjHoyklye5JEkh5M8m+S2YXv7Y3+Kubc/9knOTfKFJF8c5v6Hw/YdSQ4Nx/3e4Q9ESHLOsH5kuH92QwZaVf6s4gd4F3A18MyibX8E7BuW9wF3DMs3An/L6M+srwUOTXv86zD3DwO/s8S+VwJfBM4BdgD/BmyZ9hzWMPftwNXD8vcC/zrMsf2xP8Xc2x/74fidNyyfBRwajudngN3D9k8Cvzos/xrwyWF5N3DvRozTM/RVqqp/AP7zpM27gAPD8gHg5kXb/7JG/gnYlmT7xox08paZ+3J2AfdU1X9X1ZeAI4y+NuKMVFXHqurxYfkbwGFGn5Buf+xPMffltDn2w/H75rB61vBTwHXAwWH7ycf9xL+Hg8D1yfp/74BBn6xLquoYjP7xAxcP2zfL1yS8d7iscNeJSw40nvvwMvoqRmdrm+rYnzR32ATHPsmWJE8Cx4GHGL3ieK2q3hh2WTy/78x9uP914ML1HqNB3xhjfU3CGe4TwA8B7wCOAX8ybG859yTnAZ8F3ldVXz/VrktsO6Pnv8TcN8Wxr6r/rap3MPpU/DXA25fabbidytwN+mS9euLl9HB7fNg+1tcknMmq6tXhH/y3gT/nzZfW7eae5CxGQbu7qu4bNm+KY7/U3DfTsQeoqteARxldQ9+W5MTneRbP7ztzH+7/fsa/TLlqBn2yHgD2DMt7gPsXbf/l4S8ergVeP/HyvIuTrgv/PHDiL2AeAHYP7/rvAHYCX9jo8U3KcB30TuBwVX1k0V3tj/1yc98Mxz7JTJJtw/JbgHczeg/hEeCWYbeTj/uJfw+3AJ+v4R3SdTXtd4/P1B/g04xeXn6L0W/jWxldI3sYeGG4vaDefIf844yuuT0NzE17/Osw978a5vYUo3/M2xft//vD3J8Hbpj2+Nc4959i9NL5KeDJ4efGzXDsTzH39sce+DHgiWGOzwB/MGy/gtEvqSPAXwPnDNvPHdaPDPdfsRHj9JOiktSEl1wkqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDXxf9x1MN6TAOLqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n, bins, patches = plt.hist(df)\n",
    "for i, x in enumerate((bins > 160)):\n",
    "    if x and i < len(bins)-1:\n",
    "        patches[i].set_fc('r')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The step for the last bin represent one flow, this way we even can count visually the population in all the red bins :\n",
    "\n",
    "$$P(X > 160) \\simeq nb(X_{i} > 160) / n = 11/154 = 1/14$$\n",
    "\n",
    "Using our empiric data we compute this way the probability to have a flow superior to $160$ $ m³/s$, which tells us that we have a flow of this type in average every $3.5$ years.\n",
    "\n",
    "Nothing unusual yet, right ?...\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 211,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADMlJREFUeJzt3W+IZfddx/H3x938aRt182cS1iQ4G1m0QbQJQ4hW+iApaBJxI6SwIrpIYEGtplaxWwVbnyWibRVKy9pEVglN6jaQYKoS0gTxgVsnf5p/a8yaxnSbNTvFJm19oI39+uCebYbtzM7dmTtz3e+8XzDce849l/v7cZb3nHvunLupKiRJZ77vmfYAJEmTYdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDWxdSNf7KKLLqrZ2dmNfElJOuM99thjX62qmZW229Cgz87OMj8/v5EvKUlnvCT/Ps52nnKRpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJjb0StG1mN334FRe96Xbb5rK60rS6fIIXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNTFW0JP8VpJnkzyT5NNJzk2yI8mhJC8kuTfJ2es9WEnS8lYMepJLgd8E5qrqR4EtwG7gDuCjVbUT+Bpw63oOVJJ0auOectkKvCXJVuCtwDHgOuDg8PgB4ObJD0+SNK4Vg15VXwH+GHiZUchfBx4DXquqN4bNjgKXLvX8JHuTzCeZX1hYmMyoJUnfZZxTLucDu4AdwA8AbwNuWGLTWur5VbW/quaqam5mZmYtY5UkncI4p1zeDXypqhaq6lvAfcBPAtuGUzAAlwGvrNMYJUljGCfoLwPXJnlrkgDXA88BjwC3DNvsAe5fnyFKksYxzjn0Q4w+/HwceHp4zn7gA8D7kxwBLgTuXMdxSpJWsHXlTaCqPgR86KTVLwLXTHxEkqRV8UpRSWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITYwU9ybYkB5P8S5LDSX4iyQVJHkrywnB7/noPVpK0vHGP0P8U+Luq+hHgx4HDwD7g4araCTw8LEuSpmTFoCf5PuBdwJ0AVfU/VfUasAs4MGx2ALh5vQYpSVrZOEfoVwALwF8keSLJp5K8Dbikqo4BDLcXr+M4JUkrGCfoW4GrgU9U1VXAf3Eap1eS7E0yn2R+YWFhlcOUJK1knKAfBY5W1aFh+SCjwL+aZDvAcHt8qSdX1f6qmququZmZmUmMWZK0hBWDXlX/AXw5yQ8Pq64HngMeAPYM6/YA96/LCCVJY9k65na/Adyd5GzgReBXGP0y+EySW4GXgfeszxAlSeMYK+hV9SQwt8RD1092OJKk1fJKUUlqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpoYO+hJtiR5IsnfDMs7khxK8kKSe5OcvX7DlCSt5HSO0G8DDi9avgP4aFXtBL4G3DrJgUmSTs9YQU9yGXAT8KlhOcB1wMFhkwPAzesxQEnSeMY9Qv8Y8LvAt4flC4HXquqNYfkocOmExyZJOg0rBj3JzwLHq+qxxauX2LSWef7eJPNJ5hcWFlY5TEnSSsY5Qn8n8HNJXgLuYXSq5WPAtiRbh20uA15Z6slVtb+q5qpqbmZmZgJDliQtZcWgV9UHq+qyqpoFdgOfr6pfBB4Bbhk22wPcv26jlCStaC1/h/4B4P1JjjA6p37nZIYkSVqNrStv8qaqehR4dLj/InDN5IckSVoNrxSVpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU2c1pdzbUaz+x6c2mu/dPtNU3ttSWcej9AlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUxIpBT3J5kkeSHE7ybJLbhvUXJHkoyQvD7fnrP1xJ0nLGOUJ/A/jtqno7cC3w60muBPYBD1fVTuDhYVmSNCUrBr2qjlXV48P9bwCHgUuBXcCBYbMDwM3rNUhJ0spO6xx6klngKuAQcElVHYNR9IGLl3nO3iTzSeYXFhbWNlpJ0rLGDnqS84DPAu+rqq+P+7yq2l9Vc1U1NzMzs5oxSpLGMFbQk5zFKOZ3V9V9w+pXk2wfHt8OHF+fIUqSxjHOX7kEuBM4XFUfWfTQA8Ce4f4e4P7JD0+SNK6tY2zzTuCXgKeTPDms+z3gduAzSW4FXgbesz5DlCSNY8WgV9U/Alnm4esnOxxJ0mp5pagkNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqYsX/JFrTM7vvwWkPYcO9dPtN0x6CdMbyCF2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhNeWKT/V6Z1MZUXNKkDj9AlqQmDLklNGHRJasJz6NKUbcrPDZLpvG7VdF53g6zpCD3JzyR5PsmRJPsmNShJ0ulbddCTbAE+DtwAXAn8QpIrJzUwSdLpWcspl2uAI1X1IkCSe4BdwHOTGJi0kTbjd89Pc84vTeuFm5/qWcspl0uBLy9aPjqskyRNwVqO0Jf6Vfddv4aS7AX2DovfTPL8Gl7zTHIR8NVpD2JKNvPcYXPPf6y5T+k4eb0tP/e1vzP4wXE2WkvQjwKXL1q+DHjl5I2qaj+wfw2vc0ZKMl9Vc9MexzRs5rnD5p6/c5/u3NdyyuWfgZ1JdiQ5G9gNPDCZYUmSTteqj9Cr6o0k7wX+HtgC3FVVz05sZJKk07KmC4uq6nPA5yY0lm423WmmRTbz3GFzz9+5T1Gq+ZVTkrRZ+F0uktSEQV+lJHclOZ7kmUXrLkjyUJIXhtvzh/VJ8mfDVyQ8leTq6Y187ZaZ+4eTfCXJk8PPjYse++Aw9+eT/PR0Rj0ZSS5P8kiSw0meTXLbsL79vj/F3Nvv+yTnJvlCki8Oc//DYf2OJIeG/X7v8AciJDlnWD4yPD67IQOtKn9W8QO8C7gaeGbRuj8C9g339wF3DPdvBP6W0Z/fXgscmvb412HuHwZ+Z4ltrwS+CJwD7AD+Ddgy7TmsYe7bgauH+98L/Oswx/b7/hRzb7/vh/133nD/LODQsD8/A+we1n8S+NXh/q8Bnxzu7wbu3YhxeoS+SlX1D8B/nrR6F3BguH8AuHnR+r+skX8CtiXZvjEjnbxl5r6cXcA9VfXfVfUl4Aijr404I1XVsap6fLj/DeAwoyuk2+/7U8x9OW32/bD/vjksnjX8FHAdcHBYf/J+P/Hv4SBwfbL+3ztg0Cfrkqo6BqN//MDFw/rN8jUJ7x1OK9x14pQDjec+vI2+itHR2qba9yfNHTbBvk+yJcmTwHHgIUbvOF6rqjeGTRbP7ztzHx5/Hbhwvcdo0DfGWF+TcIb7BPBDwDuAY8CfDOtbzj3JecBngfdV1ddPtekS687o+S8x902x76vqf6vqHYyuir8GePtSmw23U5m7QZ+sV0+8nR5ujw/rx/qahDNZVb06/IP/NvDnvPnWut3ck5zFKGh3V9V9w+pNse+Xmvtm2vcAVfUa8Cijc+jbkpy4nmfx/L4z9+Hx72f805SrZtAn6wFgz3B/D3D/ovW/PPzFw7XA6yfenndx0nnhnwdO/AXMA8Du4VP/HcBO4AsbPb5JGc6D3gkcrqqPLHqo/b5fbu6bYd8nmUmybbj/FuDdjD5DeAS4Zdjs5P1+4t/DLcDna/iEdF1N+9PjM/UH+DSjt5ffYvTb+FZG58geBl4Ybi+oNz8h/zijc25PA3PTHv86zP2vhrk9xegf8/ZF2//+MPfngRumPf41zv2nGL11fgp4cvi5cTPs+1PMvf2+B34MeGKY4zPAHwzrr2D0S+oI8NfAOcP6c4flI8PjV2zEOL1SVJKa8JSLJDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6Qm/g879zZUL19swgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n, bins, patches = plt.hist(df)\n",
    "for i, x in enumerate((bins > 255)):\n",
    "    if x and i < len(bins)-1:\n",
    "        patches[i].set_fc('r')\n",
    "        \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram:\n",
    "$$P(X > 255) \\simeq nb(X_{i} > 255) / n = 3/154 $$\n",
    "\n",
    "Here again we can answer (hardly) that we therefore have a flow superior to $255$ $m³/s$ in average every $12$,$8$ years.\n",
    "\n",
    "Ok enough played, let's come to the annoying question :"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**\"You have no way to answer\" question: Now, what is the probability to observe such a flow superior to $400$ $m³/s$ ?**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAD8CAYAAABXe05zAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAC4VJREFUeJzt3W2MpXdZx/HfZbc8CVqwqyHAOsUQIzEGmg1BMcQAidI1VhNe1ARFY7KJBgSjMUuMgu+KUSImRLICBoUAUkgkFh+IQIwvXGxLgZa1UmGFQqUQw9MbAbl8MffCuMzsnLZ79nDNfD7JyZxzn7tzrn/uyXfPueecaXV3AJjjOzY9AAD3j3ADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wzJF1fNOrr766t7a21vGtAQ6kW2+99XPdfXSVfdcS7q2trdxyyy3r+NYAB1JV/eeq+zpVAjCMcAMMI9wAwwg3wDDCDTCMcAMMI9wAwwg3wDDCDTDMWj45+WBsnbp5I4977sYTG3lcgPvLM26AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYJiVwl1Vv1FVd1bVHVX15qp62LoHA2B3+4a7qh6X5NeTHO/uH05yRZIb1j0YALtb9VTJkSQPr6ojSR6R5NPrGwmAi9k33N39qSR/mOQTSe5N8oXu/od1DwbA7o7st0NVPTrJ9UmuSfL5JG+rqud39xsv2O9kkpNJcuzYsTWMul5bp27e2GOfu/HExh4bmGeVUyXPSfLx7v5sd381yTuS/NiFO3X36e4+3t3Hjx49eqnnBGCxSrg/keTpVfWIqqokz05ydr1jAbCXVc5xn0lyU5Lbknx4+W9Or3kuAPaw7znuJOnulyV52ZpnAWAFPjkJMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMMxK4a6qq6rqpqr6t6o6W1U/uu7BANjdkRX3e1WSv+vu51XVQ5I8Yo0zAXAR+4a7qr4ryTOT/FKSdPdXknxlvWMBsJdVTpU8Mclnk/x5VX2gql5bVd+55rkA2MMqp0qOJLk2yYu6+0xVvSrJqSS/u3OnqjqZ5GSSHDt27FLPyRpsnbp5Y4997sYTG3tsmG6VZ9z3JLmnu88st2/Kdsj/n+4+3d3Hu/v40aNHL+WMAOywb7i7+7+SfLKqfnDZ9OwkH1nrVADsadV3lbwoyZuWd5R8LMkvr28kAC5mpXB39+1Jjq95FgBW4JOTAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMOsHO6quqKqPlBVf7POgQC4uPvzjPvFSc6uaxAAVrNSuKvq8UlOJHntescBYD+rPuP+4yS/neTra5wFgBUc2W+HqvrpJPd1961V9RMX2e9kkpNJcuzYsUs2IAfT1qmbN/K45248sZHHhUtplWfcz0jyM1V1Lslbkjyrqt544U7dfbq7j3f38aNHj17iMQE4b99wd/dLu/vx3b2V5IYk7+nu5699MgB25X3cAMPse457p+5+X5L3rWUSAFbiGTfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMPuGu6qeUFXvraqzVXVnVb34cgwGwO6OrLDP15L8ZnffVlWPSnJrVb27uz+y5tkA2MW+z7i7+97uvm25/qUkZ5M8bt2DAbC7+3WOu6q2kjw1yZl1DAPA/lY5VZIkqapHJnl7kpd09xd3uf9kkpNJcuzYsUs24GGwdermTY/Amm3yGJ+78cTGHpv1WOkZd1Vdme1ov6m737HbPt19uruPd/fxo0ePXsoZAdhhlXeVVJLXJTnb3a9c/0gAXMwqz7ifkeQXkjyrqm5fLteteS4A9rDvOe7u/uckdRlmAWAFPjkJMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwwj3ADDCDfAMMINMIxwAwyz7/+6DA6SrVM3b3qEy+4wrnlTzt144rI8jmfcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8Awwg0wjHADDCPcAMMIN8AwK4W7qn6qqu6qqrur6tS6hwJgb/uGu6quSPLqJM9N8uQkP19VT173YADsbpVn3E9Lcnd3f6y7v5LkLUmuX+9YAOxllXA/Lsknd9y+Z9kGwAYcWWGf2mVbf8tOVSeTnFxufrmq7nowgw1ydZLPbXqIDbH2w+swr3/PtdcrHtT3/f5Vd1wl3PckecKO249P8ukLd+ru00lOr/rAB0VV3dLdxzc9xyZY++Fce3K41//tsPZVTpX8a5InVdU1VfWQJDckeed6xwJgL/s+4+7ur1XVC5P8fZIrkry+u+9c+2QA7GqVUyXp7ncledeaZ5nq0J0e2sHaD6/DvP6Nr726v+X3jAB8G/ORd4BhhHsfVfX6qrqvqu7Yse0xVfXuqvro8vXRy/aqqj9Z/jTAh6rq2s1N/uDtsfaXV9Wnqur25XLdjvteuqz9rqr6yc1MfWlU1ROq6r1Vdbaq7qyqFy/bD/yxv8jaD/yxr6qHVdX7q+qDy9p/f9l+TVWdWY77W5c3aqSqHrrcvnu5f+uyDNrdLhe5JHlmkmuT3LFj2x8kObVcP5XkFcv165L8bbbf+/70JGc2Pf8a1v7yJL+1y75PTvLBJA9Nck2S/0hyxabX8CDW/tgk1y7XH5Xk35c1Hvhjf5G1H/hjvxy/Ry7Xr0xyZjmef5XkhmX7a5L86nL915K8Zrl+Q5K3Xo45PePeR3f/U5L/vmDz9UnesFx/Q5Kf3bH9L3rbvyS5qqoee3kmvfT2WPterk/ylu7+n+7+eJK7s/3nEkbq7nu7+7bl+peSnM32J4YP/LG/yNr3cmCO/XL8vrzcvHK5dJJnJblp2X7hcT//83BTkmdX1W4fWrykhPuB+b7uvjfZ/iFP8r3L9sPy5wFeuJwOeP35UwU5wGtfXv4+NdvPvg7Vsb9g7ckhOPZVdUVV3Z7kviTvzvYriM9399eWXXau7xtrX+7/QpLvWfeMwn1prfTnAYb70yQ/kOQpSe5N8kfL9gO59qp6ZJK3J3lJd3/xYrvusm30+ndZ+6E49t39v939lGx/SvxpSX5ot92WrxtZu3A/MJ85/zJ4+Xrfsn2lPw8wWXd/ZvnB/nqSP8s3XxIfuLVX1ZXZDtebuvsdy+ZDcex3W/thOvZJ0t2fT/K+bJ/jvqqqzn/uZef6vrH25f7vzuqnFx8w4X5g3pnkBcv1FyT56x3bf3F5h8HTk3zh/Mvqg+KC87Y/l+T8O07emeSG5bfs1yR5UpL3X+75LpXlPOXrkpzt7lfuuOvAH/u91n4Yjn1VHa2qq5brD0/ynGyf439vkuctu1143M//PDwvyXt6+U3lWm36t7jf7pckb872y8KvZvtf11/J9jmsf0zy0eXrY/qbv5F+dbbPiX04yfFNz7+Gtf/lsrYPZfuH9rE79v+dZe13JXnupud/kGv/8Wy/5P1QktuXy3WH4dhfZO0H/tgn+ZEkH1jWeEeS31u2PzHb/xjdneRtSR66bH/Ycvvu5f4nXo45fXISYBinSgCGEW6AYYQbYBjhBhhGuAGGEW6AYYQbYBjhBhjm/wCT0qwEUjGt3wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n, bins, patches = plt.hist(df)\n",
    "for i, x in enumerate((bins > 400)):\n",
    "    if x and i < len(bins)-1:\n",
    "        patches[i].set_fc('r')\n",
    "        \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histogram:\n",
    "$$P(X > 400) \\simeq nb(X_{i} > 400) / n = 0 $$\n",
    "\n",
    "The answer is not satisfying, because it's surely not \"impossible\" to observe a flow superior to $400$ $m³/s$.\n",
    "\n",
    "## Partial Conclusion : we cannot estimate the probability of such a flow with the histogram!\n",
    "\n",
    "And now we precisely get to what we talked about before, we want to compute the probability of an event for which we have absolutely no data, I hope this is clear now. In statistical modelling, usually, we train an algorithm with a dataset that is supposed to have the same distribution as the testing data, in other words we have datas for what we want to model. Extreme value theory comes into play here:\n",
    "\n",
    "- **\"Classic\" statistics doesn't allow to compute probabilities such that $P(X > x)$ when $x$ is beyond the maximum of our observations.**\n",
    "\n",
    "\n",
    "- **Hopefully Extreme value theory gives us tools to extrapolate beyond our datas.**\n",
    "\n",
    "A \"famous\" joke illustrates this problematic: during the night, someone is squatting at the foot of a lamppost, and a passerby asks him: \"What are you doing?\" \"I'm looking for my keys.\" \"Are you certain that you lost them around the lamppost?\", asks again the passerby. \"No, but, in fact, that's the only lighted place.\"\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extreme value statistics is kind of like this joke, we want to make statistics where there is no light, where there is no data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory: Extreme values theorem / Implementation: scipy.stats.genextreme"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The extreme values theorem can be seen as the analog of the central limit theorem.\n",
    "For a random variable $X$ and it's maximum over n observations: $X_{n,n}$. \n",
    "\n",
    "\n",
    "\n",
    "Under general conditions on the distribution of X, there is three parameters\n",
    "$a_{n}$, $b_{n}$ and $\\gamma$ such that :\n",
    "   $$\\lim_{n -> +\\infty} P(\\frac{X_{n,n} - a_{n}}{b_{n}} \\leq x) = H_{\\gamma}(x)$$\n",
    "   where:\n",
    "$$ H_{\\gamma}(x) = \\begin{cases}\n",
    "                \\exp(-\\exp(-x))            &\\text{for } \\gamma = 0\\\\\n",
    "                \\exp(-(1+\\gamma x)^{-1/\\gamma}_{+})   &\\text{for }\n",
    "                                                        \\gamma != 0\n",
    "              \\end{cases}\n",
    "$$\n",
    "   \n",
    "   with $y_{+} = max(0,y)$.\n",
    "   \n",
    "Useful vocabulary:\n",
    " - $H_{\\gamma}$ is the **generalized extreme value (GEV) distribution**.\n",
    " - $\\gamma$ is the shape parameter\n",
    " - $a_{n}$ and $b_{n}$ are scaling parameters.\n",
    " \n",
    "With different values of $\\gamma$ there is a type I, type II and type III GEV, also called the Gumbel ($\\gamma=0$), Fréchet ($\\gamma >0$) and Weibull ($\\gamma <0$) distributions.\n",
    "\n",
    "If you have a random variable $X$ and study the distribution of it's maximum over a number of samples $X_{n,n}$, this distribution must belong and must \"fall\" in the domain of attraction of one of these 3 types of GEV. \n",
    " \n",
    "(More here : https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution.)\n",
    "\n",
    "**We will use the genextreme class from scipy**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAHwCAYAAAAIDnN0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXlcVdX2wL/bWRNxKKc00RyS4XIBcZZBFEyQVDTLIXG20hzKpFeaDb6nZU7PyvcqxdTStDRTf2WmqKQ54zyQs4IKKDiBMuzfH+dy3r1wuRcUu5r7y+d+uOecffZeezjnrLv2OnsJKSUKhUKhUCgUCsdRwtECKBQKhUKhUDzqKIVMoVAoFAqFwsEohUyhUCgUCoXCwSiFTKFQKBQKhcLBKIVMoVAoFAqFwsEohUyhUCgUCoXCwSiFTKH4myGEiBFCDDZ97yOEWFfM+bsIIaQQolRx5vt3RwgRKYSIdWD5h4QQAcWUl8W4Mo2HhsWRtym/G0KIBsWVn0LxMKAUMoWiEAghXhBCbBdC3BRCXDZ9f0UIIRwtmy2klIullMF/ZZlCiNNCiHTTQzX3M6cQ5wUIIc7/FTIWFSFEOSFEqhCivZVjM4QQyx0hl6n8XAU5t60vCSFWCyE6mqeTUrpJKWMKmZdNZbs4x5X5Dwiz/CtKKU8WR/4KxcOCUsgUCjsIIV4HZgEfAzWBGsBwoA1Q5i+W5WGxSnUxPVRzPyOKI1NH1V9KmQEsBV7KI09J4EVggSPkykNlKWVFwBP4FVghhIgs7kIeojGoUDxUKIVMobCBEMIZeB94RUq5XEp5XWrslVL2kVLeNqUrK4SYJoQ4a7JQzBVClDcdCxBCnBdCvG6yriUKIQaYlVGYc8cLIS4C84UQVUwWkCQhxFXT9zoFyK9Pkwkh3sxjtcoUQkTn1lMI8ZVJtgtCiA9NygZCiJIm+ZKFECeB0Htoz8/NrUlCiKlCiN+EEI8B/wfUNpOvthBikhBiuRBikRDiGhAphCghhIgSQpwQQqQIIb4TQlQ15Zdr4RkghDhnap/hQghfIcR+k5VrTh6ZBgohjpjS/iKEqFeA+AuACCFEBbN9IWj30f8z5ZUr13UhxGEhRLcC2iGfJSqvpagIclkgpbwopZwFTAKmCiFKmPI7LYToYPreXAixSwhxzTTmpptO32z6n2rqg1amMfS70CyBV4BJwvr0a2chxEnTOPnYrNxJQohF1uouhJgMtAPmCDNLqjCbAjWNza9N4/2MEOIds7wjhRCxpvF5VQhxSgjxrFlZkSaZrpuO9SlMGyoUjkApZAqFbVoBZYEf7aSbCjQGjEBD4ElgotnxmoCzaf8g4FMhRJUinFsVqAcMRbtu55u2nwLSAbtTglLKj3ItVkBTIAn4znR4AZBlKt8LCAZylYMhQJhpfzOgh72ybPA6YDA9KNuhtUV/KeVN4FkgwcyqlmA65zlgOVAZWAy8BnQF/IHawFXg0zzltAAaAb2AmcDbQAfADXheCOEPIIToCvwD6A48AWwBvrUmuJRyK5BoSptLP+AbKWWWafsEmoLhDLwHLBJC1CpSCxVRLhv8AFQHmlg5NguYJaWsBDzN/8aBn+l/ZVMfbDNttwBOmvKbXEB53dDGhzdanw20J6CU8m20uo2wYUn9N1p7NkDr85eAAWbHWwDHgMeBj4CvhMZjwGzgWSmlE9AaiLMnk0LhMKSU6qM+6lPAB+gLXMyzbyuQiqYI+QECuAk8bZamFXDK9D3AlLaU2fHLQMtCnnsHKGdDRiNw1Ww7Bhhs+h4JxOZJXx7YDYw3bdcAbgPlzdK8CGw0fd8ADDc7FgxI8/rkyf80cMPURrmfIWbHmwNXgDPAi2b7A4DzefKaBGzOs+8IEGS2XQvIBEoBLibZnjQ7ngL0Mtv+Hhht+v5/wCCzYyWAW0C9Aur2DrDO9L2SKa2Xjb6JA57L2xdmcpqPCfN+K7Rc1vIy7S9n2t/GrF86mL5vRlMYH7eXl0nus3nSWYwr0zmdzLZfAX4z68NFBZVhXu88+TUESqKNTVezY8OAGDM5/jQ7VsF0bk3gMbSxF4HZ2FYf9XlQP8pCplDYJgV43HxqSUrZWkpZ2XSsBJoFowKw2zQllgr8bNqv5yP/Z0UB7eFasZDnJknNhwkAIUQFIcR/TNM319AerpVzpxgLwVfAMSnlVNN2PaA0kGgmw3/QrCGgWaHOmZ1/phBldJVSVjb7fJF7QEq5A83aIvifZcYW5/Js10Pzj8qV9QiQjaZY5nLJ7Hu6le2KZnnNMsvrikmuJwuQ5WsgUAjxJJql8E8p5d7cg0KIl4QQcWb5uaNZbopKUeWyRm7aK1aODUKzyh4VQuwUQoTZyStvH9hLcwZt3Nwrj6P5aZqPuTNYtsPF3C9SylumrxWlZnXthebvmSiEWCOEeKYYZFIo7gtKIVMobLMN7Rf6czbSJKM95N3MFBBnqU0N2qMw58o857yONg3VQmpTTrnTTHbf+BRCRJnOHWS2+xxaHR83k6GSlNLNdDwRqGuW/qlC1MuWDK+iTQMnAG+aHcpbz4L2n0ObhjJX+MpJKS/chTjngGF58iovtenJ/IJIeRZtiq0P2nTl12b1qgd8AYwAqpmU9oNY75ebpv/m/mg171auAuiGZok9ZqUe8VLKF9GU7qnActMUX2H7wBp5x0julPNNCq6nvbyT0ayf5v5zTwGF6msp5S9Syo5oVtSjaP2jUDyQKIVMobCBlDIVbWrnMyFEDyFERaE5lRvRpkSQUuag3ehnCCGqAwghnhRChBQi/7s51wlNiUsVmjP7u4Wpi8nZ+TU061W6mQyJwDrgEyFEJVP9ns71s0KzYr0mhKhj8nuLKkx5BcjQGPgQbSq4H/CmqS1Bs2JVE9qLFLaYC0zOdXIXQjwhhLClMNvL6y0hhJspL2chRE875yxAU7raoPm05ZKr0CSZ8hqAZiHLh5QyCU2p6Cu0lyYGovly3YtcmNLWEEKMQBsXb5nGWN40fYUQT5iOpZp2Z5tkz0Hz1yoq44T2wkldYBTaW6mgTdv6CSGeMvXtW3nOu1RQeVLKbLTxN1kI4WTq87HAImvpzTG1Q7hJ0byNNo2efRf1Uij+EpRCplDYQUr5EdpD4E00i8MltCm98Wj+ZJi+/wn8YZpGXI91Z2prFPXcmWh+YMnAH2hTnIWhF9pU6BHxvzcZ55qOvYQ2NXQYzUl+OZpVATSF8RdgH7AHzVncHj8Jyzc6V5imfRcBU6WU+6SU8WiO6wuFEGWllEfRHNdPmqbqCprymgWsAtYJIa6b2qBFIdvAAinlCjQL0RJT2x9Ee7nAFsuBKmg+UolmeR0GPkGzql4CPIDfbeQzBBiHNvXtxv/G0t3KlSqEuAkcADoDPaWU8wpI2wk4JIS4gdaeL0gpM0xTfpOB30190NJOmeb8iOabGAesQZsaR0r5K5pytt90fHWe82YBPUxvSc62ku9INCvbSSAW+AYoqF7mlECzJiegTdv6o/m2KRQPJELKwliiFQqFQqFQKBT3C2UhUygUCoVCoXAwSiFTKBQKhUKhcDBKIVMoFAqFQqFwMEohUygUCoVCoXAwSiFTKBQKhUKhcDCl7Cd5sHj88celi4uLo8VQKBQKhUKhsMvu3buTpZRP2Ev30ClkLi4u7Nq1y9FiKBQKhUKhUNhFCFGYcHNqylKhUCgUCoXC0SiFTKFQKBQKhcLBKIVMoVAoFAqFwsHcVx8yIUQntDhlJYEvpZRT8hyPBD5GC7ILMEdK+WVRy8nMzOT8+fNkZGTco8QKhUKhUDyalCtXjjp16lC6dGlHi/JIct8UMiFESeBToCNwHtgphFhlCsBrzlIp5Yh7Kev8+fM4OTnh4uKCEOJeslIoFAqF4pFDSklKSgrnz5+nfv36jhbnkeR+Tlk2B/6UUp6UUt4BlgDP3Y+CMjIyqFatmlLGFAqFQqG4C4QQVKtWTc00OZD7qZA9CZwz2z5v2peXCCHEfiHEciFEXWsZCSGGCiF2CSF2JSUlWS1MKWMKhUKhUNw96jnqWO6nQmatZ2We7Z8AFymlAVgPLLCWkZTyv1LKZlLKZk88YXdttYeCNWvWcODAAUeLoVAoFAqF4gHgfipk5wFzi1cdIME8gZQyRUp527T5BeBzH+W5r5QsWRKj0ah/Tp8+XWDan3/+mU2bNuHu7m4zz8jISJYvX15oGVJTU/nss88KPJ6eno6/vz/Z2dmFztMap06dokWLFjRq1IhevXpx584dq+n2799Pq1atcHNzw8PDQzeFd+jQgatXrxapzEmTJhEdHU1kZCQxMTEABAQEWCwSfPr0abttej84ffo033zzjd10CQkJ9OjRA4CYmBjCwsJs5hkQEEBMTAyRkZF28zZvi86dO5Oamlpg2pkzZ3Lr1q0Cjw8ePJjDhzVXz4oVK9ot25y4uDjWrl2rb69atYopU6bYOOPe6dGjBydPngTgxo0bDBs2jKeffho3Nzf8/PzYvn17sZVVnGPs2LFjBAQEYDQaadq0KUOHDgXyt2FhmTRpEtOmTbsrWQo7huF/9zp3d3d69uxpcywVhejoaBISEuwnLIADBw4U6lpRKB5U7qdCthNoJISoL4QoA7wArDJPIISoZbYZDhy5j/LcV8qXL09cXJz+yRveKSsrS//eqVMnPvroo2I3D9tTyObNm0f37t0pWbLkPZUzfvx4xowZQ3x8PFWqVOGrr77KlyYrK4u+ffsyd+5cDh06RExMjP7mTr9+/WzK+bBR2IdZ7dq1i6Rg3y1r166lcuXKBR63pZBlZ2fz5Zdf4urqeldl51UmwsPDiYqKuqu8CsOhQ4fIzs6mQYMGgKZMVq1alfj4eA4dOkR0dDTJycn3rfx74bXXXmPMmDHExcVx5MgRRo4cCdy9QnYvFEUhy73XHTx4kDJlyjB37txCl2Prx+C9KmQeHh6cP3+es2fP3nUeCoUjuW8KmZQyCxgB/IKmaH0npTwkhHhfCBFuSvaaEOKQEGIf8BoQeb/kcQTR0dH07NmTLl26EBwcDMDHH3+Mr68vBoOBd999V0/79ddfYzAY8PT0pF+/fvr+zZs307p1axo0aGDxMLeWT1RUFCdOnMBoNDJu3Lh88ixevJjnntPeq3B1dWXp0qX6sYEDB/Ldd9/ZrZOUkg0bNuiWnv79+7Ny5cp86datW6fXB6BatWq6IhgeHs63335rtyxzKlasSPny5XF2dqZMmTJ202dnZzNu3Di9jf7zn/8AmgUlKCgIb29vPDw8+PHHHwFNyTRXEidNmsQnn3xCv3799DQAffr0YdUqi98VREVFsWXLFoxGIzNmzOD06dO0a9cOb29vvL292bp1K1CwdWXTpk26ZdXLy4vr169TsmRJqlatSpkyZXB2ds53Tnp6Oi+88AIGg4FevXqRnp6uH3NxcSE5OZmbN28SGhqKp6cn7u7uLF26lNmzZ5OQkEBgYCCBgYF6206cOJEWLVqwbdu2fJbH119/HW9vb4KCgsj14TRPk5ycjIuLC3fu3GHixIksXboUo9HI0qVLiY6OZsQI7SXqM2fOEBQUhMFgICgoSH9wRkZG8tprr+Ub54mJifj5+enWmC1btuRrB/MxfeLECbZv386HH35IiRLara1BgwaEhoYC0LVrV3x8fHBzc+O///2vnoe5FXD58uW6leXSpUt069YNT09PPD099X7Mzs5myJAhuLm5ERwcTHp6OidOnMDb21vPJz4+Hh8f2wb/xMRE6tSpo297eHhYbcO8li93d3fdAj958mSaNGlChw4dOHbsmJ7mxIkTdOrUCR8fH9q1a8fRo0dttnXeMVxY2rVrx59//mm3fc3H1+7du/H398fHx4eQkBASExNZvnw5u3btok+fPhiNRtLT0/ntt9/w8vLCw8ODgQMHcvv27YLE0OnSpQtLliwptPwKxQOFlPKh+vj4+Mi8HD58WP8+So6S/sX8N0qOyldmXkqUKCE9PT2lp6en7Nq1q5RSyvnz58snn3xSpqSkSCml/OWXX+SQIUNkTk6OzM7OlqGhoXLTpk3y4MGDsnHjxjIpKUlKKfX0/fv3lz169JDZ2dny0KFD8umnn7aZz6lTp6Sbm5tV+W7fvi1r1KihbycmJkoPDw8ppZTXr1+XdevWlRkZGfLatWt6PfJ+Dh06JJOSknQ5pJTy7NmzVsucMWOG7Nu3rwwODpZeXl5y6tSpFscbNmwok5OT8503aNAguXPnTrvtLaWU/v7+snHjxrp8TZs21WX5z3/+Iz/44AMppZQZGRnSx8dHnjx5UmZmZsq0tDQppdTrkpOTI/fs2SP9/Pz0vJs2bSrPnDkjY2Ji5HPPPSellDI1NVW6uLjIzMxMCzk2btwoQ0ND9e2bN2/K9PR0KaWUx48fl7lj1rx/zM8JCwuTsbGxUkqtL/Lmb41PPvlEDhgwQEop5b59+2TJkiX1dqtXr55MSkqSy5cvl4MHD9bPSU1NtTieCyCXLl1q0a65eQFy0aJFUkop33vvPfnqq6/mS5OUlCTr1asnpdTGfG6avNthYWEyOjpaSinlV199pbdrQeN82rRp8sMPP5RSSpmVlSWvXbuWrx38/Pzk/v37pZRS/vjjj/q1Z43c6+rWrVvSzc1NH3+PPfaYnmbZsmWyf//+Ukopn3/+eTljxgy9/NTUVHnq1ClZsmRJuXfvXimllD179pQLFy6UUkoZEBCg73/rrbfk7NmzC5RFSinnzZsnK1WqJDt16iSnT58ur169arUN3333Xfnxxx/r225ubvLUqVNy165d0t3dXd68eVOmpaXJp59+Wk/Xvn17efz4cSmllH/88YcMDAyUUhbc1nnHsC1y2yszM1OGh4fLzz77zGb7mo+vO3fuyFatWsnLly9LKaVcsmSJPo7Nx1R6erqsU6eOPHbsmJRSyn79+ul9YYvY2FgZFhZWqHoorGP+PFUUD8AuWQj95qELLv6gkmvGz0vHjh2pWrUqoFmN1q1bh5eXF6BZa+Lj49m3bx89evTg8ccfB9DTg/ars0SJEri6unLp0iWb+Tz11FMFypecnGwxjVWzZk1q167NgQMH2L59O+Hh4ZQtW5ayZctarUcu1t5ytTb1mpWVRWxsLDt37qRChQoEBQXh4+NDUFAQANWrVychIYFq1apZnPfll0VbF3jx4sU0a9YM0CxQuX5Z69atY//+/boFIC0tjfj4eOrUqcM//vEPNm/eTIkSJbhw4QKXLl3Cy8uLy5cvk5CQQFJSElWqVOGpp57iqaee4tVXX+Xy5cv88MMPREREUKqU7csmMzOTESNGEBcXR8mSJTl+/LjN9G3atGHs2LH06dOH7t27W1hNCmLz5s289tprABgMBgwGQ740Hh4evPHGG4wfP56wsDDatWtnNa+SJUsSERFh9ViJEiXo1asXAH379qV79+52ZSuIbdu28cMPPwDatPWbb76pH7M2zn19fRk4cCCZmZl07doVo9GYL8/ExEQK+6LP7NmzWbFiBQDnzp0jPj4+3/gzZ8OGDXz99deA1kbOzs5cvXqV+vXr67L4+Pjo1qrBgwczf/58pk+fztKlS9mxY4dNeQYMGEBISAg///wzP/74I//5z3/Yt29foeoCsGXLFrp160aFChUAzfIM2v1g69at9OzZU09rbl2y1tZFIT09Xa9/u3btGDRoEFBw+5qPr2PHjnHw4EE6duwIaNbGWrVq5Svj2LFj1K9fn8aNGwOaJf7TTz9l9OjRNmXLva8oFA8jfzuFbCYzHS2CBY899pj+XUrJW2+9xbBhwyzSzJ49u0B/srJly1qcbysfWy8SlC9fPt/6Mj179uS7775j/fr1zJkzB4Dr168X+OD+5ptvaNq0KampqWRlZVGqVCnOnz9P7dq186WtU6cO/v7+upLZuXNn9uzZoytkGRkZlC9fvkB57xUpJf/+978JCQmx2B8dHU1SUhK7d++mdOnSuLi46O3So0cPli9fzsWLF3nhhRf0c/r168fixYtZsmQJ8+bNs1v2jBkzqFGjBvv27SMnJ4dy5crZTB8VFUVoaChr166lZcuWrF+/nmeeecZuOfZ8EBs3bszu3btZu3Ytb731FsHBwUycODFfunLlyhXarzC3zFKlSpGTkwNw1+sWmctvbZz7+fmxefNm1qxZQ79+/Rg3bhwvvfSSRR7m49rNzU1v89wpy1xiYmJYv34927Zto0KFCgQEBOjnmctRmLqYy1qyZEl9ujgiIoL33nuP9u3b4+PjY1PZy6V27doMHDiQgQMH4u7uzsGDB/OlMW/rvDJaGwM5OTlUrly5wB9W1tq6KFj78Wmrfc3Hl5QSNzc3tm3bZrOMu5EL7v99RaG4n6hYln8hISEhzJs3jxs3bgBw4cIFLl++TFBQEN999x0pKSkAXLly5a7ycXJy4vr161bPqVKlCtnZ2RY3865duzJv3jxu3bql+7s4OTlZvJxg/nF1dUUIQWBgoG55WrBgge7Dk1fG/fv3c+vWLbKysti0aZPuKC6l5OLFi/lefChOQkJC+Pzzz8nMzATg+PHj3Lx5k7S0NKpXr07p0qXZuHEjZ86c0c954YUXWLJkCcuXL9d95EDzu5k5U1P03dzc8pWVt93T0tKoVasWJUqUYOHChXbfaj1x4gQeHh6MHz+eZs2a6f4+tvDz82Px4sUAHDx4kP379+dLk5CQQIUKFejbty9vvPEGe/bssSqvLXJycvS+/uabb2jbti2g+ant3r0bwMK30VberVu31v17Fi9erOdVEGfOnKF69eoMGTKEQYMG6fKb07RpU92H6emnn6ZZs2a8++67+gM9Pj6eH3/8kbS0NKpUqUKFChU4evQof/zxh55HjRo1OHLkCDk5ObqFByAoKIjPP/8c0Cw5165dsylvuXLlCAkJ4eWXX2bAgAH6/rfeessi31x+/vlnfXxevHiRlJQUnnzyyXxt6OLiotd9z549nDp1CtDGwIoVK0hPT+f69ev89NNPAFSqVIn69euzbNkyQLve7Fne8pZ54cIF/cdTYbDVvuY0adKEpKQkXSHLzMzk0KFD+WR45plnOH36tN63CxcuxN/f364cx48fd8ib1gpFcaAUsr+Q4OBgevfuTatWrfDw8KBHjx5cv34dNzc33n77bfz9/fH09GTs2LF3lU+1atVo06YN7u7uVp36g4ODiY2N1berVauGq6sr/fv3L1I9pk6dyvTp02nYsCEpKSn6lMWqVat0C0yVKlUYO3Ysvr6+GI1GvL29defq3bt307JlS6tTf4MHD7ZwKL9bBg8ejKurK97e3ri7uzNs2DCysrLo06cPu3btolmzZixevNjCEuXm5sb169d58sknLaZRatSoQdOmTS0esuYYDAZKlSqFp6cnM2bM4JVXXmHBggW0bNmS48ePW1hJrTFz5kzc3d3x9PSkfPnyPPvss3br9/LLL3Pjxg0MBgMfffQRzZs3z5fmwIEDNG/eHKPRyOTJk3nnnXcAGDp0KM8++6zu1G+Lxx57jEOHDuHj48OGDRv0/n3jjTf4/PPPad26tcVbjIGBgRw+fFh3SDdn9uzZzJ8/H4PBwMKFC5k1a5bNsmNiYvQXHb7//ntGjRqVL01oaKi+DApoU94XL16kYcOGeHh4MGTIEGrXrk2nTp3IysrCYDAwYcIEWrZsqZ8zZcoUwsLCaN++vUW/z5o1i40bN+Lh4YGPj4+uONiiT58+CCH0l3hA64eaNWvmS7tu3Tq930NCQvj444+pWbNmvjaMiIjgypUrGI1GPv/8c30az9vbm169emE0GomIiLCwbC9evJivvvoKT09P3NzcLF5MsUbeMZyYmGh3at4cW+1rTpkyZVi+fDnjx4/H09MTo9GovywRGRnJ8OHDMRqNSCmZP38+PXv2xMPDgxIlSjB8+HAAJk6cqL9YY37PAdi4caN+n1EoHjbE3ZqGHUWzZs1k3gf2kSNHaNq0qYMkenjYu3cv06dPZ+HChfq+yMhIwsLCLCxC95tRo0YRHh5epF/gjuTWrVt4eHiwZ88eq288KhxHeno6gYGB/P777/e8nEtxMG3aNNLS0vjggw/0fSEhIfzyyy8OlKrozJkzh6eeekr3S3sYuH37Nv7+/sTGxhZJmVRYop6nxY8QYreUspm9dGrUPkJ4eXkRGBhIdna2Qx9e7u7uD40ytn79egYOHMjYsWOVMvYAUr58ed577z0uXLhg86WWv4Ju3bpx4sQJNmzYYLH/YVPGAH2pkoeJs2fPMmXKFKWMKR5alIVMoVAoFAoFoJ6n94PCWsiUD5lCoVAoFAqFg1EKmUKhUCgUCoWDUQqZQqFQKBQKhYNRCplCoVAoFAqFg1EKWTFx6dIlevfuTYMGDfDx8aFVq1ZWF4O8G8wDNBeW3ADTeZFS0r59e7uLXNrjypUrdOzYkUaNGtGxY0euXr0KwOrVqy2CphfE6dOnCQgIICYmRg/mXFA9O3fuTGpqqs388gbEziUuLo61a9daPcfWMXN27dqlhymy1xe59YmOjmbSpEl28zbvp9atW9tM+89//tPm8dx2KiiIuS1iYmL09aAA5s6dq4cNuh/kHYf38/qxhXnfFieRkZEWC+b+XSjovmJOVlYW//jHP2jUqBFGo1FfB+9BJiYmBmdnZ11eo9HI+vXrHS2WVVJTU/nss8/07YSEBLvLFpnfEw4cOKDfcxUPFkohKwaklHTt2hU/Pz9OnjzJ7t27WbJkCefPn3e0aPlYu3Ytnp6eVKpU6Z7ymTJlCkFBQcTHxxMUFMSUKVMAbaHOVatWcevWreIQF9BkNo/DWRSKQyFr1qwZs2fPvqvyi4K5QmSNghQyKSU5OTn31E55FbLhw4fnC1NUnJiPQ0deP39V3z4I2IsYUVy88847JCQkcODAAeLi4tiyZYsekeBBpl27dhaRSTp06OBokaySVyGrXbt2kZR/Dw8Pzp8/z9mzZ++HeIp7QClkxcCGDRsoU6aMvpI0QL169Rg5ciSQ36oSFhamry5esWJFxo8fj4+PDx06dGDHjh0EBATQoEEDfTVq0IL1durUiSZNmvDee+/p+xctWqSvxj5s2DC7N93FixfroY46duzI9OnT9WPvv/8+H330UaHq/OOPP+or/Pfv35+VK1cCWmxEOijIAAAgAElEQVS9gIAAVq9ebfP8kiVLUrVqVcqUKWN3fS/zX+UffPABzzzzDB07duTFF19k2rRperply5bRvHlzGjduzJYtW7hz5w4TJ05k6dKl+VaOt3Zsx44dtG7dGi8vL1q3bs2xY8cATVnJDVpuzrJly/SV1v38/AD0+pQvX56KFSvmOyclJYXg4GC8vLwYNmyYRcy+3PSJiYn4+flhNBpxd3dny5YtREVF6UGd+/Tpw+nTp2natCmvvPIK3t7enDt3zqKdsrKy6N+/PwaDgR49eugKsnmaXbt2ERAQwOnTp5k7dy4zZszAaDSyZcsWJk2apLdtXFwcLVu2xGAw0K1bN90aGhAQwPjx4y3aHODQoUP6mDQYDMTHx+drB/NxaO/6OX36NO3atcPb2xtvb29dcczbLyNGjCA6OhrQ4oO6urpiMBh44403Cuwv8zwK6v/o6Gi6d+9Op06daNSokUVQ9KIgpWTcuHG4u7vj4eGhj8ecnBxeeeUV3NzcCAsLo3PnzlYfsDExMQQEBNCjRw+eeeYZ+vTpo4+f3377DS8vLzw8PBg4cKAeTNzFxYX333+ftm3bsmzZMgICAhgzZgx+fn40bdqUnTt30r17dxo1aqRHcgAtrJqPjw9ubm7897//LXQdb926xRdffMG///1vPYark5OThbW4oLwLcy+Mjo6ma9eudOnShfr16zNnzhymT5+Ol5cXLVu21MPOffHFF/j6+uLp6UlERMRd/0DcuXMnBoOBjIwMbt68iZubGwcPHiQmJgY/Pz+6deuGq6srw4cP1+ONfvvtt3h4eODu7s748eMt6vf222/j6elJy5Yt9cDuSUlJRERE4Ovri6+vL7///jsAkyZNYuDAgXob5P5wiIqK4sSJExiNRsaNG2dh/SroWslLly5d9DBmigcIKeVD9fHx8ZF5OXz48P82Ro2S0t+/eD+jRuUr05xZs2bJ0aNHF3h8/vz58tVXX9W3Q0ND5caNG6WUUgJy7dq1Ukopu3btKjt27Cjv3Lkj4+LipKenp35+zZo1ZXJysrx165Z0c3OTO3fulIcPH5ZhYWHyzp07UkopX375ZblgwQIppZT16tWTSUlJ+WR56qmn5LVr16SUUt66dUs2atRI3r59W+bk5MiGDRvKhIQEKaWUbdu2lZ6envk+v/76q5RSSmdnZ4t8K1eurH9ftGiRHDFihJRSyh9//FFOmDDBZvsV1E655NZl586d0tPTU966dUteu3ZNNmzYUH788cdSSin9/f3l2LFjpZRSrlmzRgYFBdnM09qxtLQ0mZmZKaWU8tdff5Xdu3eXUkq5ceNGGRoamu8cd3d3ef78eSmllFevXi1UHUeOHCnfe+89KaWUq1evloDeT4899piUUspp06bJDz/8UEopZVZWlt5fucellPLUqVNSCCG3bduWr51OnTolARkbGyullHLAgAF6O5mPi507d0p/f38ppZTvvvuunibvtoeHh4yJiZFSSjlhwgQ5ynQ9FNTmI0aMkIsWLZJSSnn79m1569atfO1gPg7tXT83b96U6enpUkopjx8/LnPvAeb9IqWUr776qpw/f75MSUmRjRs3ljk5OVLK//WNtf4yz6Og/p8/f76sX7++TE1Nlenp6fKpp56SZ8+eLVBeKaXs37+/XLZsmcW+5cuXyw4dOsisrCx58eJFWbduXZmQkCCXLVsmn332WZmdnS0TExNl5cqV852bK2ulSpXkuXPnZHZ2tmzZsqXcsmWLTE9Pl3Xq1JHHjh2TUkrZr18/OWPGDCml1t9Tp07V8/D395dvvvmmlFLKmTNnylq1asmEhASZkZEhn3zySZmcnCyllDIlJUVKKfX7Te7+gu4ruezbt08ajUabbVNQ3oW9Fz799NPy2rVr8vLly7JSpUry888/l1JKOXr0aL3euXlKKeXbb78tZ8+ebVOm3LY1v9f9+eef+vmvv/66fOWVV+Q///lPPX3ZsmXliRMnZFZWluzQoYNctmyZvHDhgqxbt668fPmyzMzMlIGBgXLFihV6/VatWiWllHLcuHHygw8+kFJK+eKLL8otW7ZIKaU8c+aMfOaZZ6SU2jXYqlUrmZGRIZOSkmTVqlXlnTt35KlTp6Sbm5suu/l2QddK3nNiY2NlWFiY1baweJ4qigVglyyEfqOWNL4PvPrqq8TGxlKmTBl27txpM22ZMmXo1KkToJmSy5YtS+nSpfHw8OD06dN6uo4dO1KtWjUAunfvrocH2b17N76+voAWRqZ69eo2y7ty5QpOTk6Atsp5+/bt+eWXXyhfvjxNmjTRY/nlWjvuhurVq5OQkABAeHh4sYVfiY2N5bnnnqN8+fKA9ivPnO7duwPg4+Nj0XaFJS0tjf79+xMfH48Qwu40S5s2bYiMjOT555/Xy7bH5s2b+eGHHwBterdKlSr50vj6+jJw4EAyMzPp2rUrRqPRal716tUrMGZg3bp1adOmDQB9+/Zl9uzZuqWoKKSlpZGamqoHdu7fvz89e/bUj1tr81atWjF58mTOnz+vW1/yYj4O85L3+snMzGTEiBHExcVRsmRJjh8/blPmSpUqUa5cOQYPHkxoaKhuAbPXX7b6PygoSLfkurq6cubMGerWrWtTjrzExsby4osvUrJkSWrUqIG/vz87d+4kNjaWnj17UqJECT2WZUE0b96cOnXqAGA0Gjl9+jROTk7Ur19fj3HZv39/Pv30U0aPHg1Ar169LPLIvR49PDxwc3PTr/kGDRpw7tw5qlWrxuzZs3UfvnPnzhEfH6/ff4rC/PnzmTVrFikpKWzdupW6desWmHdh74WBgYE4OTnh5OSEs7Ozfh/w8PBg//79ABw8eJB33nmH1NRUbty4QUhIiF1Z27VrZ9WyP3HiRHx9fSlXrpzF9Hbz5s1p0KABAC+++CKxsbGULl2agIAAnnjiCUCLbbp582a6du1KmTJl9LHo4+PDr7/+CmjRQA4fPqzne+3aNT3IemhoKGXLlqVs2bJUr15dt6oVRGGvFfN7tOLB4e+nkM2c+ZcX6ebmxvfff69vf/rppyQnJ9OsmbYwb6lSpXRzNkBGRob+vXTp0gghAChRogRly5bVv2dlZenpctOYb0sp6d+/P//6178KLWuuLCVKaLPVPXv2JDo6mqysLIvg2e3atdNvCuZMmzaNDh06UKNGDRITE6lVqxaJiYkWimBGRoauNBUn0k5Uidy2K1mypEXbFZYJEyYQGBjIihUr9JcObDF37ly2b9/OmjVrMBqNxMXFFeqhlbcv8+Ln58fmzZtZs2YN/fr1Y9y4cVb9uWwFLbc2XsByLJqPw7vFWpv37t2bFi1asGbNGkJCQvjyyy9p3769xXnm49De9TNjxgxq1KjBvn37yMnJ0afCCrquSpUqxY4dO/jtt99YsmQJc+bMYcOGDVb7yxxb/Z9bz7x1LQoFjd+C9m/fvp1hw4YBmjtBpUqVrMph77rIO07M7zHm+eXec2JiYli/fj3btm2jQoUKBAQEFHqsNGzYkLNnz3L9+nWcnJwYMGAAAwYMwN3dnezsbJt5F/ZemFdma+kiIyNZuXIlnp6eREdHWwSgLypXrlzhxo0bZGZmkpGRobdnQffkgjCvn/kYysnJYdu2bVbvmUUddwVdK3m5X/doxb2hfMiKgfbt25ORkcHnn3+u7zP3WXBxcSEuLo6cnBzOnTvHjh07ilzGr7/+ypUrV0hPT2flypW0adOGoKAgli9fzuXLlwHtxnHmzBmb+TRp0oSTJ0/q2wEBAWzdupXY2FgLi9OWLVssHFzzOrqGh4ezYMECABYsWKD7AwEcP368yG/5FYa2bdvy008/kZGRwY0bN1izZo3dc5ycnKwqltaOpaWl8eSTTwLovki2OHHiBC1atOD999/n8ccf59y5c3bP8fPzY/HixQD83//9n+6PZc6ZM2eoXr06Q4YMYdCgQezZswfQbuiFdY4+e/Ys27ZtAzSflrZt2wLaWNy9ezeAhRJUUDs5OztTpUoV3WK6cOFC3VpWECdPnqRBgwa89tprhIeH61YLc8zHob3rJy0tjVq1alGiRAkWLlyo+0nWq1ePw4cPc/v2bdLS0vjtt98AuHHjBmlpaXTu3JmZM2fqipe9/ipq/wO89NJLhb6e/fz8WLp0KdnZ2SQlJbF582aaN29O27Zt+f7778nJyeHSpUu68tCiRQv9urNlZX7mmWc4ffo0f/75J1C4PrJFWloaVapUoUKFChw9epQ//vjDarqgoCAuXLhgsa9ChQoMGjSIESNG6IpWdnY2d+7cKVLe98r169epVasWmZmZ+vUGsGLFCt56660i5TV06FA++OAD+vTpY+ETtmPHDk6dOkVOTg5Lly6lbdu2tGjRgk2bNpGcnEx2djbffvut3b4IDg5mzpw5+nbeHwp5sXVPK+haycv9ukcr7g2lkBUDQghWrlzJpk2bqF+/Ps2bN6d///5MnToV0KZK6tevj4eHB2+88Qbe3t5FLqNt27b069cPo9FIREQEzZo1w9XVlQ8//JDg4GAMBgMdO3YkMTHRZj6hoaEWvxZLlixJ+/bt6datG2XKlCm0PFFRUfz66680atSIX3/9laioKP3Yxo0bCQ0NBWDVqlVMnDix0PlGR0dTp04d/WP+pp2vry/h4eF4enrSvXt3mjVrZveFgMDAQA4fPpzPqd/asTfffJO33nqLNm3aFOqNtHHjxunOu35+fnh6eto9591332Xz5s14e3uzbt06qwGxY2JiMBqNeHl58f333zNq1ChAezAYDAb69Oljt5ymTZuyYMECDAYDV65c4eWXX9bLHzVqFO3atbMIMN+lSxdWrFihO/Wbs2DBAsaNG4fBYCAuLs5ufy5duhR3d3eMRiNHjx61at0zH4f2rp9XXnmFBQsW0LJlS44fP65bKOrWrcvzzz+vt4mXlxegPYzDwsIwGAz4+/szY8YMwH5/FbX/Afbv369P+eVl2LBh+jhu1aoV3bp1w2Aw4OnpSfv27fnoo4+oWbMmERER1KlTB3d3d4YNG0aLFi2KFMi+XLlyzJ8/n549e+Lh4UGJEiUsXpAoKp06dSIrKwuDwcCECROsTovn5OTw559/UrVq1XzHJk+eTK1atXB3d8fLy4t27drRv39/ateuXai8i4MPPviAFi1a0LFjR5555hl9/4kTJwp8w3zLli0Wy14sX76cr7/+mlKlStG7d2+ioqLYuXOnHjy+VatWREVF4e7uTv369enWrRu1atXiX//6F4GBgXh6euLt7W3xY9Uas2fPZteuXRgMBlxdXZk7d67N9NWqVaNNmza4u7szbtw4i2MFXSt5Mb9HKx4gCuNo9iB97Dr1K2ySkJAgO3ToYLEvr0P3vXDx4kXZvn37YsnLGtevX5dSas6rPj4+cvfu3fetLMX9w9o4fNhIS0uTPXr0KJa8csd1cnKybNCggUxMTCyWfO8XBw4ckGPGjHG0GEWmT58+8vLly/ecT94XSh4mMjIyZIsWLfQXWPKinqfFD8qpX2GNWrVqMWTIEK5du3bPa5FZ4+zZs3zyySfFnm8uQ4cO5fDhw2RkZNC/f/+7sjYqHM/9Hod/BZUqVWLZsmXFkldYWBipqancuXOHCRMmULNmzWLJ937h7u5usWTOw8KiRYscLYLDOXv2LFOmTKFUKfX4f9AQ0o5D6INGs2bNZN4V2Y8cOULTpk0dJJFCoVAoFH8P1PO0+BFC7JZSNrOXTvmQKRQKhUKhUDgYZbNUKBSKhwAp4cYNuHoVbt2CnBxtX7lyUKkSODtDEd7LUSgUDxhKIVMoFIoHGCkhJQUuXIDMTBACHnsMSpfWjucqaUJA9epQqxYo9yCF4uFDXbYKhULxgHLrFpw5AzdvakpY3bqaJcxsxRKkhIwMuHwZLl3SlLd69cBKEAiFQvEAo3zIioExY8Yw0yxCQEhICIMHD9a3X3/9dbtvJLVu3RooOJC1LcwDQUdGRloNTAwwevRoNm/eXKS87XHq1ClatGhBo0aN6NWrl74AZF72799Pq1atcHNzw8PDQ180skOHDlYXR81LZGSkHlw5N4yKi4sLHh4e+rpBW7duJSEhgR49etjNr3PnzqSmppKamspnn31W+Aqb+Oc//1modLnlAFaDjZvj4uJi8d8W5n0+ceJE1q9fX2DalStXWoRmycvcuXP5+uuvAW2h4Lwvzdgib/sVtv3vhZkzZxYor3mgZVuY98vs2bNp2rRpodZ3Ky6io6N54oknMBqNuLm5WQSAz+XqVTh6FO7cARcXuHgxhqNHt1ooY6BZxsqX15QwV1fYunU1b775LklJf1l1FApFMaAUsmKgdevWbN26FdAWTExOTubQoUP68a1bt+pxBQsi9/z7xZUrV/jjjz/w8/Mr1nzHjx/PmDFjiI+Pp0qVKnz11Vf50mRlZdG3b1/mzp3LoUOHiImJobRpvqVfv353pRDlsnHjRn0189atW1O7du0CFVJz1q5dS+XKle+7QpZbzv3k/fff1yMoWMOWQpaVlcXw4cOtLt5aGPK2X2Hb/27Jyspi3rx59O7d+57yMe+Xzz77jLVr11qs6G5PhuKgV69exMXFcejQIcqUKaMvXCwlJCbCiROaouXqCo8/Dps2xdi9T1SoAEOHhvL776s4duwWdtaJVigUDxBKISsG2rRpo98oDx06hLu7O05OTly9epXbt29z5MgRfRXxjz/+GF9fXwwGA++++66eh7n15Nq1a3Tr1g1XV1eGDx+ux+szT7N8+XIiIyMLLePy5cv1wL0TJkzg+eef149t3ryZzp07F7neUko2bNigW0T69+/PypUr86Vbt26dvkI5aCtN564SHx4ezrfffmu3LGdnZ8qUKUPVqlUtVpjPi7mFJDo6mu7du9OpUycaNWrEm2++qadzcXEhOTmZqKgoTpw4gdFo1Fe9LqiPcomKiiI9PR2j0ahbVbp27YqPjw9ubm7897//zVeOOYmJifj5+WE0GnF3d9dXxs8NSJz7Py+TJ0+mSZMmdOjQgWPHjun7za2iUVFRuLq6YjAYeOONN9i6dSurVq1i3LhxGI1GTpw4QUBAAP/4xz/w9/dn1qxZFtY20NZqat26Ne7u7npYoLxp3N3dOX36dL72M2//jIwMBgwYgIeHB15eXmzcuNFmv2RnZxMZGYm7uzseHh76CvvmbNiwAW9v70KtoVSY/h8+fDgnT54kPDycGTNmcOXKFbp27YrBYKBly5Z62KdJkyYxdOhQgoODeemll4iOjqZr16506dKF+vXrM2fOHKZPn46XlxctW7bkypUrduXLJSsri5s3b+qB5g8eTKJPnwgGDvTlpZd82bHjd06fPs3cuXOZMWOGHk3hp59+okWLFnh5edGhQwc98HSpUoLg4AD27l3NhQuQZ/gpFIoHlL+dD9no0WAnFFiRMRptxyyvXbs2pUqV4uzZs2zdupVWrVpx4cIFtm3bhrOzMwaDgTJlyrBu3Tri4+PZsWMHUkrCw8PZvHlzPqvVjh07OHz4MPXq1aNTp0788MMP9zwN9Pvvv+t5fPDBB3Tp0oXt27fTokULvvrqKz2w+JgxY/QHpzkvvPCCRXgkgJSUFCpXrqw/HOvUqZMvth1ocdOEEISEhJCUlMQLL7ygPxyrVKnC7du3SUlJoVq1anTu3Jkvv/yS2rVrW+Qxa9YsAH744QeL/YGBgZQsWZKyZcuyffv2fGXHxcWxd+9eypYtS5MmTRg5ciR169bVj0+ZMoWDBw/q8eMK00dTpkxhzpw5FjHn5s2bR9WqVUlPT8fX15eIiIgCA41/8803hISE8Pbbb5Odna1PVe3cudPivzm7d+9myZIl7N27l6ysLLy9vfHx8bFIc+XKFVasWMHRo0cRQpCamkrlypUJDw8nLCzMYgylpqayadMmQFM2zLl58yZbt25l8+bNDBw4kIMHD1qth7X2y51OBi1IOMCBAwc4evQowcHBHD9+HLDeL5cvX+bChQt6eblTiub8/vvv+eptC3v9P3fuXH7++Wc2btzI448/zsiRI/Hy8mLlypVs2LCBl156Sa/b7t27iY2NpXz58kRHR3Pw4EH27t1LRkYGDRs2ZOrUqezdu5cxY8bw9ddfM3r0aJuyLV26lNjYWBITE2ncuDFdunTh0iUYN24UL788hh492nLu3FlCQkI4cuQIw4cPp2LFirzxxhsAXL16lT/++AMhBF9++SUfffSRviizr28ztm3bQnDw85w5o1nabMSiVygUDwB/O4XMUeRaybZu3crYsWO5cOECW7duxdnZWfcPW7duHevWrdOtZTdu3CA+Pj6fQta8eXMaNGgAwIsvvkhsbOw9K2SJiYkWlpc+ffrw3Xff0bRpUzZs2KBbdaxZJQrC2qLCQoh8+7KysoiNjWXnzp1UqFCBoKAgfHx8CAoKAqB69eokJCRQrVo11q5dW6R65T5ICyIoKEiPC+jq6sqZM2csHsh5KWwf5WX27NmsWLECgHPnzhEfH1+gQubr68vAgQPJzMyka9euGI1Gm3mDFmevW7duVKhQAcBqsOlKlSpRrlw5Bg8eTGhoqE1fxF69ehV47MUXXwS0YNjXrl2zqhgVhtjYWEaOHAloAbDr1aunK2TW+sXNzY2TJ08ycuRIQkNDCQ4OzpdnYmKixaKV1sab+b6i9n9sbKwedL19+/akpKSQlpYGaG1evnx5PW1gYCBOTk44OTnh7OxMly5dAPDw8LAaUD0vvXr1Ys6cOUgpefXVV3nvvY8JD49i5871XLx4mMmTtXTXrl2zGkz6/Pnz9OrVi8TERO7cuUP9+vX1Y9WrVycxMYEGDeDIEW36s2nT/72ZqVAoHjz+dgqZLUvW/STXj+zAgQO4u7tTt25dPvnkEypVqsTAgQMBTYF56623GDZsmM288j5kcrfN9+c6xReW8uXLW5wTFhbGpEmTaNSoEV27dqVs2bKAfQtZSEgIly5dolmzZnzxxRekpqaSlZVFqVKlOH/+fD7LFmiWM39/f11x6ty5M3v27NEVsoyMDIsHXXGSWy/QAqnb8/8pbB+ZExMTw/r169m2bRsVKlQgICDAZv/4+fmxefNm1qxZQ79+/Rg3blyhfLisKR/mlCpVih07dvDbb7+xZMkS5syZowdCzktBQYetlSOEoFSpUvrUORRu/NmKAmKtX6pUqcK+ffv45Zdf+PTTT/nuu++YN2+exXl5x3G1atUsXgq5cuWKhYJ+N/2fl9z2yNtm5nmXKFFC3y5RokSR/MyEEHTs2IVp0/5Nnz5RCJHDtm3b7F4TI0eOZOzYsYSHhxMTE2Nh6cy9pkqXhqef1l4OOHMGGjYstFgKheIvRvmQFRNt2rRh9erVuo9T1apVSU1NZdu2bbRq1QrQ3r6cN28eN27cAODChQtcvnw5X147duzg1KlT5OTksHTpUtq2bQtAjRo1OHLkCDk5Obo1prA0bdqUP//8U9+uWLEirq6uTJw4UVcYQbOQ5TrJm39ypyt/+eUX4uLi+PLLLxFCEBgYqPsvLViwgOeeey5f2SEhIezfv59bt26RlZXFpk2bcHV1BbQH4MWLFwv1ZuH9wMnJycL6UNg+Kl26NJmZmQCkpaVRpUoVKlSowNGjR/njjz9slnnmzBmqV6/OkCFDGDRoEHv27LErp5+fHytWrCA9PZ3r16/z008/5Utz48YN0tLS6Ny5MzNnztSn2vLW0R65zuWxsbE4Ozvj7OyMi4uLLueePXs4deqU3bz9/Px0R/njx49z9uxZmjRpUmC5ycnJ5OTkEBERwQcffGC1XfKO44CAABYtWqQrUgsWLCAwMLDQdbUlc0xMDI8//vg9xdqcM2cOc+bMsZkmKwvWro3lqaeepkEDCA4OtjinoH5MS0vjySefBLR6m3P8+HHdl++xx6B2bUhN1T4KheLBRClkxYSHhwfJycm0bNnSYp+zs7P+iz04OJjevXvTqlUrPDw86NGjh9WHWatWrYiKisLd3Z369evTrVs3QPPXCQsLo3379tSqVatI8oWGhhITE2OxLyIigho1aujTc3fD1KlTmT59Og0bNiQlJYVBgwYBsGrVKiZOnAhofmJjx47F19cXo9GIt7c3oaGhgOaX07JlS90PrXPnziQkJNy1PEWlWrVqtGnTBnd3d8aNG1foPho6dCgGg4E+ffrQqVMnsrKyMBgMTJgwwWIMWCMmJgaj0YiXlxfff/89o0aNsiunt7c3vXr1wmg0EhERQbt27fKluX79OmFhYRgMBvz9/fXp5xdeeIGPP/4YLy8vTpw4YbesKlWq0Lp1a4YPH66/NRsREcGVK1cwGo18/vnnNG7c2Gr7mfPKK6+QnZ2Nh4cHvXr1Ijo62sKqlJcLFy4QEBCA0WgkMjKSf/3rX/nSPPvssxZLtwwdOhQnJyc8PT3x9PTkxo0buo/V3TBp0iR27dqFwWAgKioqn6JTVI4ePVrg1PXSpUtNL3YYOHx4L//85wRKl9amv3NlcHV1Ze7cuQB06dKFFStW6E79kyZNomfPnrRr1y7ftP3GjRv1awygRg1tRf+zZyE7+56qpFAo7hMquPgjRNu2bVm9erX+un9MTAzTpk1j9erVDpNp1KhRhIeH69OXCoU9unXrxkcffUSjRo0cLYpdwsLC+OGHHyhTQEyjy5c1JaluXU1pKg4uXbpE7969+e233yz2X78Ox45pK/mbDGsKRT7U87T4UcHFFfn45JNPOHv2rKPFsMDd3V0pY4oiMWXKFBIfkgW2Vq9eXaAydvs2nD+vxaGsXr34yjx79qz+tqU5Tk5QrRpcvKit7K9QKB4slIVMoVAo/mKkhOPHtZBI7u5/XVDwO3fgwAFNMXOQ26biAUc9T4ufR85C9rAplgqF4tElKUmbQqxb969TxkAr64kntHiXt2//deUqHg7Uc9Sx/C0UsnLlypGSkqIGk0KheODJzIQLF7SpShtL6N03atbU/j8ks54pwUQAACAASURBVL6KvwgpJSkpKZQrV87Rojyy/C3WIatTpw7nz58nSUXTVSgUDzgpKXDjhrYUxdGjjpEhI0Mr+/p1KEQUKsUjQrly5ahTp46jxXhk+VtciqVLl7ZYpVqhUCgeRLZvh3bt4M03YepUx8lx/ry2YOygQWAWG16hUDiQv8WUpUKhUDzo5OTAyJHashPvvONYWerUgb59YcECMEWGUigUDkYpZAqFQvEXsGgR7NwJH32kLUHhaF55BW7dgoULHS2JQqGAv8myFwqFQvEgk5EBTZpobzju2AElHpCfws2ba/5shw6BnVCpCoXiLnnklr1QKBSKB5W5c7UV+adMeXCUMYCXX4YjR8AsGpVCoXAQD9CtQaFQKP5+XLsGkydDhw7a50GiVy+oUkU59isUDwJKIVMoFIr7yLRpkJwMVmKlO5wKFSAyEn74QQuppFAoHIdSyBQKheI+ceUKzJgBPXpAM7seJI5h6FDIyoIlSxwtiULxaKMUMoVCobhPzJypOc2/+66jJSmYZ54Bb2/45htHS6JQPNoohUyhUCjuA6mpMGsWdO+uBRB/kOndW1uSIz7e0ZIoFI8uSiFTKBSK+8C//6059E+Y4GhJ7NOrl7bsxbffOloSheLRRSlkCoVCUcxcu6b5joWHg9HoaGnsU6cO+Ptr05YP2dKUCsXfBqWQKRQKRTEzdy5cver4EElFoXdvOHYM9u51tCQKxaOJUsgUCsUDyVWuspnNfMmXHOawo8UpNHfuaL5jQUHg6+toaQpPRASULq2c+xUKR6EUMoVC8UBxm9u8xEtUpSr++DOEIbjhRic6sZGNjhbPLt98AwkJMG6coyUpGlWrQkgILFumpi0VCkegFDKFQvHAkEoqnejEQhYyhjGsZS2HOcyHfMh+9tOe9sxlrqPFLBAptYVgDQYIDna0NEXnuee0EE8HDjhaEoXi0aOUowVQKBQKgOtcpx3tOMYxFrGIPvTRj73N27zO6zzP87zMy2SQwWhGO1Ba6/zf/2mBur/++uEM1h0Wpv1ftUpTKhUKxV+HspApFIoHgg/5kIMc5Cd+slDGcilHOZaznAgiGMMYPudzB0hpm48/1t5YfOEFR0tyd9SsCS1awE8/OVoSheLRQylkCoXC4RzjGDOYwQAGEEJIgenKUIYlLCGUUMYwhgM8OHNrcXEQEwOvvaY5xz+sdOkCO3ZAYqKjJVEoHi2UQqZQKByKRDKa0ZSnPP/CfgTuUpRiPvOpTGV605t00v8CKe3z739rwboHD3a0JPdGeLj2f/Vqx8qhUDxqKIVMoVA4lJ/4iZ/5mUlMogY1CnXOEzxBNNEc5CBRRN1nCe2TnAyLF8NLL0GVKo6W5t5wdwcXFzVtqVD81SiFTKFQOJTJTKYxjRnBiCKd14lOjGQks5nNDnbcJ+kKxxdfwO3bMHKkQ8UoFoTQrGS//gq3bjlaGoXi0UEpZAqFwmEc5CA72MFwhlOaojteTWYyNajBWMYiccziWZmZ8Omn0KEDuLo6RIRip0sXyMiA335ztCQKxaODUsgUCoXDmM98SlOavvS9q/OdcOJDPuR3fud7vi9m6QrHihVw4QKMGuWQ4u8L7dpB+fKwfr2jJVEoHh2EfMiWZG7WrJnctWuXo8VQKBT3yB3uUIc6+OHHcpbfdT7ZZOOFFze48f/s3XmczWX/x/HXZQZD1mHse0QolaW6UxShRSplSFJR/ShuN3cL2ve9tN9JKzUYFEmyNLTchVIKibKEG2OXfcz398flWMfMWb7f7zkz5/38PTzObeac63Opfrxd1/X9XCxiEUkkuTjLvLVubZupLl0KCQm+lvZUhw7w11+2r5qIhM8Y84PjOM3yep9WyEQkKiYzmUwyuZmbIxongQSe53mWs5xXeMWl2QVn8WKYNQtuu61ghTGwW7CLFtnVPxHxngKZiETF27xNFarQjsjvGGpLW9rRjqd5ml34dxL9jTdsz7GbI8uUMaltW/uqc2Qi/lAgExHfrWMdn/EZN3IjiS7d4DaUoWSSyQhGuDJeXnbuhPfeg2uugZQUX0r66vTT7a9L58hE/KFAJiK+m8IUsskmlVTXxryAC2hJS57mafaxz7VxTyQtDbZtgz59PC8VFYUKQZs2NpDls6PGIvmSApmI+G4qU6lMZU7jNFfHHcpQVrOaD/jA1XFz8sYb0KgRtGzpeamoadvWXqG0eHG0ZyJS8CmQiYivDnCAaUyjHe0wGFfHbk97mtKUJ3mSLLJcHftI8+fDvHn2ML9x95cQUwLnyKZNi+48ROKBApmI+OpHfmQzm3O9RDxcBsNgBrOMZUzCu7t/3nkHihaF7t09KxETataEevUUyET8oEAmIr6aylQMhra09WT8TnSiJjV5iZc8GX/PHhg5Eq66CpKTPSkRUy66CL76Cg4ciPZMRAo2BTIR8dVUpnIWZ5GCN48mJpJIX/qSQQa/8Ivr40+cCFu2FMxWFzk5/3zYvh1+/TXaMxEp2DwNZMaYDsaYJcaYZcaYe3J53zXGGMcYk2cnWxHJv7aznf/yX0+2K4/Um94Uoxgv87LrY48YATVq2JWjeBB4aOGrr6I7D5GCzrNAZoxJAF4FLgEaAt2MMcddvWuMKQn0B773ai4iEhtmMpMDHHClGWxukknmeq5nJCPZzGbXxl21yp6nuvHGgteZ/0Rq1oTq1eHrr6M9E5GCzcsVshbAMsdx/nQcZx+QBnTK4X2PAE8Dezyci4jEgC/4ghKU4FzO9bxWP/qxm928xVuujfnee7Yn1403ujZkvtCypV0hUz8yEe94GciqAn8d8fPVB792iDHmTKC64zifejgPEYkRX/EV53M+RSjiea3TOI0LuID/8B+yyY54vOxs+3RlmzZQu7YLE8xHWraEtWthxYpoz0Sk4PIykOXUnefQ36+MMYWAF4BBeQ5kzK3GmHnGmHmZmZkuTlFE/LKDHSxkIWdztm81b+M2/uRPvuTLiMfKyIDly+PnMP+Rzj/fvmrbUsQ7Xgay1UD1I35eDVh7xM9LAo2BDGPMCuAcYGJOB/sdx3nTcZxmjuM0SymIl8aJxIEf+REHhxa08K3m1VxNMsm8yZsRj/X221C6tG13EW8aNYIyZXSwX8RLXgayuUA9Y0xtY0wRoCswMfBNx3G2OY5T3nGcWo7j1AK+A65wHGeeh3MSkSiZwxwAmtPct5pJJNGDHkxgApmEv7q+dSuMGwfXXQfFirk4wXyiUCE47zytkIl4ybNA5jhOFnAHMBVYDIxxHGehMeZhY8wVXtUVkdg0hznUoQ7lKe9r3Vu4hf3s533eD3uMtDTbELZXLxcnls+0bGnvtNy4MdozESmYjJPPHptp1qyZM2+eFtFE8pua1OQf/IOP+Mj32udxHpvYxGIWh3V/ZvPmsG8f/PRTwb67Mjdff23Pkk2YAFdeGe3ZiOQfxpgfHMfJs8+qOvWLiOfWsY5VrPL1/NiRbuEWlrCErwl9z23RInuR+E03xW8YA2jWDBITYc6caM9EpGBSIBMRz81lLkDUAtk1XMNJnMQHfBDyZ0eOtE1gu3XzYGL5SFISNGmiQCbiFQUyEfHcHOaQQAJncmZU6pegBFdzNWMYw54QelBnZ8OoUdCuHVSs6OEE84kWLWDuXPvPRUTcpUAmIp77nu85jdMoTvGozaEHPdjGNiYxKejPfP21vS7p+us9nFg+0ry5vWj899+jPRORgkeBTEQ8lU02c5kbte3KgIu4iCpUCWnb8oMP4KSToFNOl77FoRYH/xVq21LEfQpkIuKpZSxjK1ujHsgSSKA73ZnClKB6ku3ZA2PHwtVX21Am0KABlCihQCbiBQUyEfHUT/wEwFmcFeWZwA3cQBZZpJGW53snT4Zt27RdeaSEBPu0pQKZiPsUyETEUwtZSCEKcSqnRnsqNKYxZ3BGUE1iR46ESpXgoot8mFg+0qKF7ce2d2+0ZyJSsCiQiYinfuVX6lKXJJKiPRXArpLNYx6LWXzC92zebFfIunWzvbfksBYtYP9++PnnaM9EpGBRIBMRTy1kIY1oFO1pHNKNbhSiUK6H+8eOtaFD25XH08F+EW8okImIZ/awh2Usi6lAVolKtKc9oxhFNjk31Bo5Ek49Fc6MTtu0mFatmt3KVSATcZcCmYh4ZglLOMCBmApkYHuSrWIVs5l93PeWL7f9x3r0iO+rkk7EGLtKpkAm4i4FMhHxzEIWAvYwfSzpRCdKUjLHw/0ffmhfr7vO50nlI02b2uawf/8d7ZmIFBwKZCLimYUsJJFETuGUaE/lKMUpzjVcQzrp7GLXoa87jt2uvOACqFkzihOMcWecYf9ZLVgQ7ZmIFBwKZCLimYUspB71KEKRaE/lOD3owQ52MJGJh742fz789ht07x7FieUDgbN18+dHdx4iBYkCmYh45ld+jbnzYwEXcAFVqMJoRh/6WlqabXNxzTVRnFg+UK0aJCfbfmQi4g4FMhHxxC528Sd/xmwgSyCBLnThMz5jG9vIzobRo6F9exs25MSMsatkWiETcY8CmYh44jd+w8GJuQP9R+pKV/axj4/5mO++g1WroGvXaM8qfzjjDPj1V9uvTUQip0AmIp4IPGEZqytkAC1oQS1qkUYaaWmQlARXXBHtWeUPZ55pr09asiTaMxEpGBTIRMQTv/IrhSlMXepGeyonZDB0pStfHJhB2phsLrsMSpWK9qzyhzPOsK/athRxhwKZiHhiIQupT30KUzjaU8lVV7qSPaslmesLabsyBPXr2xVFHewXcYeuzRURTyxhCWcSo3cP/f47TJoEM2Zw+l9/ccPS/oxnB5cOOR9eS4bzzoMOHeDss3W7+AkkJsJpp2mFTMQtWiETEdftZz8rWEE96kV7Koc5Dnz6KTRvbpd3/v1vWLmS/XUa8DHX0uCUTzBnVLft5x9/HFq2tP0dHn8ctmyJ9uxj0pln2hUyx4n2TETyPwUyEXHdKlaRRRYnc3K0p2L99BOccw507AibN8NLL8GKFbBwIdNvG8v2vWWY91wab45pay9p3LgRxoyBJk1g6FCoUcMGs337ov0riSlnnGGz6qpV0Z6JSP6nQCYirvuDPwCif6DfceDFF+3W419/wVtv2Vb8/foduhspLQ3KlIEm7TaQRpr9XNmycO21MHWqDXPt2tlgdtZZ8O23UfwFxZZAx36dIxOJnAKZiLhuGcuAKAey3bvhqqvgX/+y58EWLIBevaBw4aPe8vHH0LkzXFfkGr7jO5az/OhxmjSBcePsmbMdO+D88+GZZ7RPhz1DZowCmYgbFMhExHXLWEYxilGZytGZwPbtcMklMHGiXSH7+GMoX/64t02ZYjNW167QhS4AjGFMzmNefrnthNq5M9x1F3TpYj8cx046CerUgYULoz0TkfxPgUxEXLeMZZzMyRiM/8U3b4Y2beCbb2DUKPjnP+0yTg7S0qBCBWjdGmpRi3M5l4/46MRjlyxp71d69lmYMAEuusieN4tjjRvbnCoikVEgExHX/cEf0dmu3LvXblMuWGADU7duJ3zr33/bhy6vvfZwZ4uudOVnfmYxi09cwxgYNAg++cQmkVatYM0al38h+UfjxrB0qf1HLyLhUyATEVdlkx2dQOY40Ls3zJ4N775rtxhz8emn9gxZly6Hv3Yt12IwjGZ03vUuuww+/9w+LHD++bB6dWTzz6caNYKsLNvaTUTCp0AmIq5awxr2stf/QPbwwzByJDz6aK4rYwHp6VCpku0BG1CZylzABYxmNA5BHNpv1QpmzLDblu3axeX2ZeODd8dr21IkMgpkIuKqwBOWvvYgmz4dHnwQbrgBhgzJ8+07d8Jnn9nz+QkJR38vlVR+4zd+JciE0by5fQLzzz/h0kvj7qB//fp2y1eBTCQyCmQi4irfe5BlZtogduqp8PrrJzzAf6QpU+x25TXXHP+9znSmEIWC27YMaNUKxo6FH3+0h9KyskL4BeRvRYrAKacokIlESoFMRFy1jGUUpjDVqe59McexvcU2bYKPPoLixYP6WHo6pKTYo1/HqkAFLuKi4LctAzp2tIFw6lTbFiOONG6s1hcikVIgExFXLWMZdahDAgl5vzlSb71ltwuffto2cA3C7t32QP/VVx+/XRmQSirLWMZ8Qrw5+5ZboH9/eOEFGDEitM/mY40a2R3bnTujPROR/EuBTERcFehB5rn16+1KVOvWNgQFaepUGxxy2q4MuJqrSSQxtG3LgOeeswf8+/Sx92LGgcaN7WLl4ly6hYhI7hTIRMQ1Do5/LS/+/W+brII8NxaQng7lytljXyeSTDIXczFjGBPatiXYE+5paVC1qu2psXlzaJ/Ph/SkpUjkFMhExDUb2MDf/O19IJsxw7a4uOceaNAg6I/t3WtvU7ryyqOutMxRKqmsYAVzCGOVq2xZ29F/7Vq48cYCf+/lySdD0aI6RyYSCQUyEXGNL5eK798PffvaFDB4cEgfnTbNdqXIbbsy4EqupAhFwtu2BGjRwl5CPmkSPP98eGPkEwkJ0LChVshEIqFAJiKu+ZM/AahDHe+KvPmmbQs/bBgUKxbSR8eOhTJl7BWUeSlNaTrQgTGMIZvs8Obavz906mR7oy1YEN4Y+USjRgpkIpFQIBMR16xiFQA1qOFNgR074KGH7EH+Sy8N6aP79tnrJzt1sr2zgpFKKmtYw7d8G/pcwZ5tGz7cbmFef32BvvCxUSN7e9T27dGeiUj+pEAmIq5ZyUoqUIFihLZyFbRnn7WNYJ9+OqSD/GCPnW3bFtx2ZUBHOpJEUvjblmAbno0YAb/8AvfdF/44MS5wlG/JkujOQyS/UiATEdesYpV3q2Pr1tmWEl262OuKQpSeDiVLwsUXB/+ZkpTkMi4jnXQOcCDkmodcdhncdpsNlF9/Hf44Max+ffuqQCYSHgUyEXHNSlZ6F8gef9xu+T32WMgf3b8fPv4YrrjCPg0YilRSWcc6ZjM75LpHefZZqFEDeveGPXsiGysGnXyyPdz/22/RnolI/qRAJiKucHBYxSpqUtP9wdets2exevaEuqE/wZmRYduBhbJdGXAZl3ESJ0W2bQlQooR9IGHJEnjkkcjGikFFithQpkAmEh4FMhFxxWY2s4td3qyQPf+8PZV/zz1hfTw93eah9u1D/2xxitORjoxjHFlEeGl4u3Zw003w1FMwP8RrmfKB+vW1ZSkSLgUyEXHFSlYCuL9CtmkTvPYadOsW1upYVhZMmACXXx5yl4xDUkllIxuZyczwBjjSc8/Zg/633AIHIjiXFoMaNIClSwvcL0vEFwpkIuIKz1peDBtmr0gaMiSsj3/1lX0wM5ztyoAOdKAUpSLftgTbAuP55+GHH+zl6AVIgwb2mN/KldGeiUj+o0AmIq7wZIVsxw546SW4+mrbCj4M6elQvDhcckn400giiU50Yjzj2ce+8AcK6NrVXqY5ZIhdASwgAk9a6hyZSOgUyETEFatYRTGKUY5y7g367ru2edhdd4X18QMHYPx420O2ePHIppJKKlvZyjSmRTYQ2B5qr7xif21hrvzFokAvMgUykdApkImIKwJPWBpCa9h6QtnZdrvy3HPh7LPDGuLbb+0DmpFsVwZczMWUpaw725YAjRtDv3726dF589wZM8rKlbM/dLBfJHQKZCLiCtd7kE2eDH/8AQMGhD1EejokJYV8y1KOilCEq7iKj/mYPbjUR+zBB6FCBbj9dhtAC4AGDbRCJhIOBTIRcYXrXfpffBGqV7fnx8KQnQ3jxkGHDrZDvxtSSWUHO/icz90ZsHRpeOYZmDMH3nnHnTGjTIFMJDwKZCISsT3sYT3r3TvQv2ABzJwJd9wBiYlhDfH997BmjTvblQEXcRHlKOfetiXYS8dbtrQ91jZvdm/cKKlfHzZsgC1boj0TkfxFgUxEIvYXfwEutrx4+WV7Cr9377CHGDvWdo+//HJ3pgSQSCKd6cwkJrGLXe4MGjjgv3lzgbh8XJeMi4RHgUxEIhboQebKCtn27fDRR7YRbHJyWEM4jj0/1q6d3RV0Uyqp7GQnk5ns3qBNmkCfPvCf/+T7/T49aSkSHgUyEYlYoAeZKytkH35oG8HeemvYQ8ydC3/95e52ZUArWlGRiu5uWwI88IBdFQzzeqhYUbs2FC6sFTKRUCmQiUjEVrEKg6Ea1SIbyHHsKlGTJtC8edjDpKfbo2dXXBHZdHKSQALXcA2TmcwOdrg3cEqKDWOffGKvF8inEhNtKFu2LNozEclfFMhEJGIrWUkVqlCYwpENNG8e/PQT3HabPVsVhsB2Zdu29pYiL6SSyh72MIlJ7g48YABUrQp33ml/IflU3boKZCKhUiATkYi51vLizTfttt1114U9xPz5sHy5N9uVAedxHlWp6v62ZfHi8Mgj9hHR9HR3x/ZRIJDl40wp4jsFMhGJmCuBbMcOe5i/a9eITuKnp0NCAnTqFNl0clOIQlzLtXzO52xlq7uD33ADnHYaDB4M+1y4NzMK6taFv/+27S9EJDgKZCISEQeHNayJ/PxYero9zN+rV/hzcWy7i4sugvLlI5tOXlJJZR/7+IRP3B04IQGeftreUvDGG+6O7ZO6de2rti1FgqdAJiIR2cpWdrObqlSNbKD334d69ezdlWH65RcbArzcrgw4m7OpSU33ty0B2reHNm3g4YftBeT5jAKZSOgUyEQkImtYA0AVqoQ/yIoVkJFht+vCPMwPdpGtUCG48srwpxIsg6ELXZjGNDaxyeXBjV0l27QJnnzS3bF9ULOmXehTIBMJngKZiERkLWsBIlshGznSvl5/fURzSU+HVq3sfd1+SCWVLLKYwAT3Bz/rLOjeHYYNg//9z/3xPVSkiA1lCmQiwVMgE5GIBFbIwg5kjmO3K1u3hlq1wp7HokWweLE/25UBZ3EWJ3OyN9uWAA89ZA/2P/GEN+N7qG5dWLo02rMQyT8UyEQkIoFAVpnK4Q3w3Xf2T+4bbohoHunpdqfvqqsiGiYkBkMqqcxkJhvw4JHCk0+Gm2+2zXJXrXJ/fA+p9YVIaBTIRCQia1lLOcqRRFJ4A7z/PhQrFvHSVno6tGwJlcPMheFKJZVsshnHOG8K3HuvfX3kEW/G90i9evZ5hE0uH68TKagUyEQkImtYE/525f79NkldcQWULBn2HJYssU9Y+rldGXAap9GABt5tW9aoAf/3f/DOO/nqUJaetBQJjQKZiERkLWvDf8Jy5kzYuBG6dYtoDoGm9ldfHdEwYQlsW85m9qEHHFw3eLA9Kf/QQ96M7wEFMpHQKJCJSEQiWiFLS7Nd+Tt0iGgO6em2fVm1CHvThiuVVBwc0vHouqNKlaBfPxg1yj69kA/Urm3P9CmQiQRHgUxEwpZFFutZH94K2d69MH68PYVftGjYc1i2zN5HHo3tyoBTOZXTOM27bUuAu+6CEiXggQe8q+GiokXtbqsCmUhwFMhEJGzrWU822eGtkH3+OWzfbu+ujMC4g2fpO3eOaJiIpZLKt3zLX/zlTYFy5eBf/7LLgfPne1PDZYEnLUUkbwpkIhK2iHqQjR5tL5y86KKI5pCeDs2b20ak0ZRKKgBjGONdkX/9C8qWhfvu866GixTIRIKnQCYiYQv72qRdu2DiRLusVbhw2PVXrIB586K7XRlQl7o0pSkf8qF3RcqUgX//GyZPhrlzvavjkrp1bduLrVujPROR2OdpIDPGdDDGLDHGLDPG3JPD9//PGPOLMeYnY8zXxpiGXs5HRNwV9rVJU6bAzp2QmhpR/cB2ZSwEMoDudOdHfuQ3fvOuyB132FWyfNCXrHZt+7p8eXTnIZIfeBbIjDEJwKvAJUBDoFsOgetDx3FOcxznDOBp4Hmv5iMi7lvDGhJJJIWU0D44bpzdrjz//Ijqp6fbKx/r1IloGNd0pSuFKMQoRnlXpFQpGDgQJk2CH3/0ro4LFMhEguflClkLYJnjOH86jrMPSAM6HfkGx3G2H/HTkwBdsiGSj6xhDZWpTKFQfivZuxc+/RSuvBISE8Ou/ddf9talWFkdA3t9VBvaMIpROF7+dtavn92+jPFVssDVpApkInnzMpBVhaMeN1p98GtHMcbcboz5A7tC1t/D+YiIy9ayNvTtyunTYceOiLu4jh9vX6P9dOWxutOd5Sznv/zXuyKlS8OAAfDxx/Dzz97ViVDZsnZBT4FMJG9eBjKTw9eO+yuj4zivOo5zMnA3cG+OAxlzqzFmnjFmXmZmpsvTFJFwrWFN6Af6x42zf0q3aRNR7fR0OP10OOWUiIZx3VVcRTGKebttCdC/v/3nGMOrZMbYbUsFMpG8eRnIVgPVj/h5Ncj1XpE04MqcvuE4zpuO4zRzHKdZSkqIZ1VExDMhr5BlZcEnn0DHjvYqoHDrroVvvomt7cqAUpTiCq5gNKPZz37vCpUta0PZuHH2Is8YVbu2fRpWRHLnZSCbC9QzxtQ2xhQBugITj3yDMabeET+9DFjq4XxExEU72ck2toUWyGbNgs2bI96unDABHCc2AxnYbctNbGIqU70t9K9/2e79jz7qbZ0IBAKZoxPCIrnyLJA5jpMF3AFMBRYDYxzHWWiMedgYc8XBt91hjFlojPkJGAj09Go+IuKusHqQjR8PxYq5cndlw4Zw6qkRDeOZ9rSnHOW837ZMTrYH/MeOjdk7LmvXtm3nNmyI9kxEYpunfcgcx/nMcZxTHMc52XGcxw5+7X7HcSYe/N//dBynkeM4ZziOc6HjOAu9nI+IuCfkHmSOYw+hd+gAxYuHXXf9epg9O3ZXxwCKUIQudOETPmEHO7wtNnCg/ecZo6tketJSJDjq1C8iYQl5hezHH+3hr06d8n5vLj7+GLKzYzuQgd223M1uJjDB20Lly8Ptt0NaGvz+u7e1wqBeZCLBUSATkbD8j/8BtvdWUCZOhEKF4NJLI6o7dqx9srJx44iG8dw/+Ae1qOX9tiXYVbKiReGpp7yvFSKtkIkER4FMRMKyjnUUoxilKBXcECEC0gAAIABJREFUByZOhH/8AyJ4UjozEzIy7OqYyamxTgwxGLrTnelMZx3rvC1WsSL07g3vvw+rVnlbK0QlSth/5XrSUiR3CmQiEpZ1rKMSlTA5thw8xl9/wU8/wRVX5P3eXHzyCRw4EPvblQHd6U422aSR5n2xO++0r88+632tEKkXmUjeFMhEJCyBQBaUSZPsa4SBLD3d3lt5xhkRDeObUzmVpjTlAz7wvliNGtCjBwwfbp98iCG1aimQieRFgUxEwhJSIJs40R78ql8/7HqbN8OMGflju/JIN3ADP/Ijv/Kr98XuucfeFfrii97XCkHt2nYn9cCBaM9EJHYpkIlIWIIOZNu3w8yZEa+OTZxoG/3nl+3KgG50I5FE3uM974udcgpcey28+ips3ep9vSDVrg3798OaNdGeiUjsUiATkZDtYx+b2BRcIJs2zf5p3LFjRDXHjoWaNaFZs4iG8V0KKVzO5XzAB2SR5X3BIUPs5e2vvOJ9rSCp9YVI3hTIRCRkG7Bt14MKZFOmQOnS9gnLMG3danNdftuuDOhJT9azni/4wvtiTZrAZZfZbcudO72vF4RAINOTliInpkAmIiFbjz00nmcgcxz4/HO4+GJITAy73ief2EW2a68Ne4ioupRLKUc5f7YtAYYOhU2b4M03/amXhxo1bJDWCpnIiSmQiUjIAn218gxkv/5qDw5dcklE9caOtX+ot2gR0TBRU4QidKc7H/MxW9jifcFzz4XWrW0LjL17va+Xh6JFoUoVBTKR3CiQiUjIgg5kU6bY1/btw661dSt88UX+3a4M6ElP9rGP0Yz2p+DQofaqqvff96deHmrUsO3oRCRnCmQiErJAIKtAhdzfOGUKnH46VA3yAvIcTJyYv7crA87kTBrTmHd515+CbdpA8+bw5JP28dQoq1Ej5i4REIkpCmQiErJ1rKMMZUgi6cRv2r4dvv7ale3K6tXh7LMjGibqDIae9OR7vmcJS3woaOwq2Z9/wpgx3tfLQ40asHq1vRheRI6nQCYiIQuqB9nMmXZlJoJAtm1bwdiuDOhOdxJI8O9wf8eO0KgRPP541JNQjRr2OFtmZlSnIRKzFMhEJGRBBbIpU6BkyYjaXUycCPv25f/tyoDKVKY97Xmf9zmAD23rCxWyfckWLjx8fVWUVK9uX7VtKZIzBTIRCVmegcxxbCBr2xYKFw67TkHZrjxST3qyhjXMZKY/Bbt0sReAPv64/fcSJTVq2FcFMpGcKZCJSMjyDGSLFtlH6iLcrpw61W5XFipAv1NdwRWUoYx/h/sTE+Guu2DOHPjyS39q5kCBTCR3Bei3ORHxw98H/y/XQBZod9GhQ9h1Jk0qWNuVAUkkcR3XMY5x/vQkA+jZEypXtqtkUZKcDMWLq/WFyIkokIlISILq0v/55/YweeDgUBjGjoVq1QrWdmVAb3qzl72MYpQ/BZOSYOBAmDED5s71p+YxjFHrC5HcBBXIjDHjjDGXGWMU4ETiXJ5NYf/+G776KqLtyu3bC+Z2ZcCZnMlZnMVwhuPg07mu226DsmXhiSf8qZcDBTKREwv2t7rXgeuApcaYJ40xDTyck4jEsDwD2cyZdq8xwu3KvXsL3nblkW7hFhawgB/4wZ+CJUvCHXfAhAmweLE/NY9RvboCmciJBBXIHMeZ7jhOd+AsYAUwzRjzrTHmJmNM+I9QiUi+k2cgmzIFTjoJWrYMu8bYsba5/znnhD1EzOtGN4pRjOEM969o//72INdTT/lX8wg1asD69TFxvaZIzAl6M8AYUw64EegNzAeGYQPaNE9mJiIxaR3rKEQhylP++G86jj0/1qaNvVE6DNu32yEK6nZlQGlK04UufMiH/M3f/hQtXx5uvRVGjYKVK/2peYTAk5arV/teWiTmBXuGbDzwFVAc6Og4zhWO44x2HKcfUMLLCYpIbFnHOlJIIYGE47+5ZAmsWBHR+bFPPy3425UBvenN3/zNWMb6V3TQIHvC/tln/at5kFpfiJxYsH//fMtxnIaO4zzhOM7/AIwxRQEcx2nm2exEJObk2oPsiy/sa/v2YY8f2K4899ywh8g3zuM8GtCAt3jLv6LVqkGPHvDWW7Bhg391USATyU2wgezRHL72XzcnIiL5Q66BbPp0OPlkqF07rLF37LBH0Dp3LtjblQEGQ2968y3fsohF/hW+6y67DDlsmH81sVkQFMhEcpLrb3nGmErGmKZAMWPMmcaYsw7+aI3dvhSROHPCQJaVBRkZ9rqkMMXTdmVAD3pQmML+rpLVr29T7yuv2CsRfJKUBBUqKJCJ5CSvv4O2B54FqgHPA88d/DEQGOLt1EQk1jg4rGc9Fal4/DfnzrVLXG3ahD3+2LFQpUpE95HnOxWoQCc68T7vsxcfHz8cPNg+QfH66/7VxG5bqlu/yPFyDWSO47znOM6FwI2O41x4xI8rHMcZ79McRSRGbGMb+9mfcyCbPt0eFr/wwrDG/vvv+NquPNIt3MImNvEJn/hX9Kyz7Fm/F16A3bt9K6vmsCI5y2vL8vqD/7OWMWbgsT98mJ+IxJAN2EPgFahw/DenT4czz7StFcIwcSLs2QNdukQyw/ypLW2pSU1/e5KBXSXbsAHeece3koFA5vh0QYFIfpHX30NPOvhaAiiZww8RiSMnDGQ7d8J//xvR+bG0NHvoO562KwMKUYje9GY60/md3/0rfMEF9nHWp5+G/ft9KVmjhv3PZfNmX8qJ5Bt5bVn+5+DrQzn98GeKIhIrThjIvvrK/oEe5vmxLVtsM9jU1PjbrgzoTW8SSeQN3vCvqDEwZIhtEpuW5kvJwJOWag4rcrRgG8M+bYwpZYwpbIyZYYzZeMR2pojEiRMGsunToUiRsK9LmjDB5rmuXSOdYf5ViUp0pjPv8i672OVf4csug9NOgyefhOxsz8tVrWpf16zxvJRIvhLs30XbOY6zHbgcWA2cAtzp2axEJCatZz0AKaQc/Y3p0+G88+w9iWFIS7Pty5o2jXSG+Vtf+rKFLYxmtH9FjYF77oFFi+yt7h5TIBPJWbCBLHCB+KXAR47jaPdfJA5tYAPJJFP40G8JQGYm/Pxz2NuVGzbAjBl2dcwYlyaaT53P+TSiEa/xmr+Fu3SBOnXg8cc9P21fubJ9VSATOVqwgWySMeY3oBkwwxiTAuzxbloiEos2sOH47cqZM+1rmAf609PtTlk8b1cGGAx96cs85jGXuf4VTkyEO++EOXPgyy89LVWkiG0Oq0AmcrSgApnjOPcA5wLNHMfZD+wEOnk5MRGJPTkGsunToXTpsPcb09KgUSNo3NiFCRYA13M9JSjh/yrZjTdCpUrwxBOel6paVYFM5FihPM90KpBqjLkBuAZo582URCRW5RjIZsywzWATE0Meb/Vq+4CmVscOK0Uprud60khjMz6eDklKgoEDbcCeN8/TUgpkIscL9inLD7BXKLUEmh/80czDeYlIDDoukP35JyxfHvb5sTFj7GtqqguTK0D60Ic97OFd3vW38P/9H5Qp4/kqWbVqCmQixwr2r7TNgIaOo97KIvFqP/vZzOajA9n06fY1zPNjaWl2p7NePRcmWICczum0pCWv8zoDGEChkDYzIlCyJNxxBzz6KCxeDKee6kmZqlVh0yZ7M0NSkiclRPKdYP+//FegkpcTEZHYtpGNAMcHsqpVoX79kMf74w97H7m2K3PWl74sYxnTme5v4f79oVgxeOopz0oEWl+sXetZCZF8J9hAVh5YZIyZaoyZGPjh5cREJLYEepAdulg8O9s+Ydm2bVj9KkYfbLUVj3dXBuNqrqYCFXiVV/0tnJICt94Ko0Z5dgu4epGJHC/YQPYgcCXwOPDcET9EJE4c16X/55/tvlOY58fS0mwv2Ro13JphwVKUotzCLUxiEn/yp7/FBw2yIfvZZz0ZXoFM5HjBtr2YBawACh/833OBHz2cl4jEmOMCWaD/WBiBbOFC+OUXbVfmpS99SSCBl3nZ38LVq8P118Pw4bZzr8sUyESOF+xTlrcA6cB/Dn6pKvCxV5MSkdhzXCDLyIBTToEqVUIea/Roe4n4Nde4OMECqApVSCWVEYxgO9v9LX733bB3Lwwb5vrQpUvbW7YUyEQOC3bL8nbgPLC/IziOsxSObUYkIgXZBjZQmMKUpjQcOACzZ0OrViGP4zh2u/LCC20fUsndAAawgx28wzv+Fq5fH66+Gl59Fba7GwaNUS8ykWMFG8j2Oo6zL/ATY0wioBYYInEk0IPMYOz5se3boXXrkMeZPx+WLtV2ZbCa0YyWtGQYwzjAAX+LDx4M27bB66+7PrQCmcjRgg1ks4wxQ4BixpiLgbHAJO+mJSKx5qimsBkZ9jWMFbK0NNvU/+qr3ZtbQTeAASxnOZP8/m23aVNo1w6efx5273Z1aAUykaMFG8juATKBX4DbgM+Ae72alIjEnqMC2axZULfu4dPZQcrOtufH2rWD5GQPJllAdaITNanJi7zof/HBg+3B/nfc3TKtWtX2IVO7cREr2Kcss7GH+Ps6jnON4zjD1bVfJL6sZ73tQRY4PxbGduV339nWVtquDE0iifSnP7OYxXzm+1u8VSs45xx45hnIynJt2KpVYd8+2LjRtSFF8rVcA5mxHjTGbAR+A5YYYzKNMff7Mz0RiQUOzuEVsgULYOvWsLcrixaFTp08mGQB14telKCE/6tkxsCQIbBihf0X6BK1vhA5Wl4rZAOwT1c2dxynnOM4ycDZwHnGmH95PjsRiQl/8zd72GMD2axZ9oshBrKsLLtdedllUKqUB5Ms4EpTmpu4iY/4iHWs87f4ZZdB48b20vHsbFeGVCATOVpegewGoJvjOMsDX3Ac50/g+oPfE5E4cFQPsowMqFPHNg8NwfTp9ijS9dd7MME40Z/+ZJHFa7zmb+FCheCee2DRIpjozq15CmQiR8srkBV2HOe4HX7HcTKBwt5MSURizaFAll0+7PNjo0ZBmTJw6aUuTy6O1KUuHenI67zOLnb5Wzw1FU4+GR591JWT+JUq2d1QBTIRK69Ati/M74lIARIIZDV/2Q5btoQcyHbuhAkT4Npr7RkyCd+/+Tcb2eh/o9jERPvE5Q8/wJQpEQ9XuDBUrAirV7swN5ECIK9A1sQYsz2HHzuA0/yYoIhEXyCQVZ71u/1CiOfHPvnEhjJtV0auJS05l3N5lmfJwr2nHoPSowfUrAmPPOLKKlmg9YWI5BHIHMdJcBynVA4/SjqOoy1LkTgRCGSlMn6E2rWhRo2QPj9ypD1y1rKlF7OLLwbDPdzDClYwhjH+Fi9SxJ4l++47mDEj4uEqVYL1612Yl0gBEGxjWBGJYxvYQKnsEiTM/ibk1bENG+CLL6B7d3s2XCJ3OZfTkIY8xVM4ft9id9NNdmnr4YcjHqpyZVjn8wOjIrFKvz2KSJ4yyeQfC8vApk0hnx8bPdr2ku3e3Zu5xaNCFOIu7mIBC/icz/0tXrQo3H03fPXV4RYoYQqskB3w+YpOkVikQCYiecokk7YZifYnIa6QjRoFTZrYNlbinm50oxrVeIqn/C/eu7c9kR/hKlmlSratmbr1iyiQiUgQMsnknFl77YHuWrWC/tzSpfD991od80IRijCIQcxiFt/xnb/FixWDu+6CmTPh22/DHqZSJfuqbUsRBTIRCUKms4HTZ4Xe7uLDD22vqW7dvJlXvOtNb8pSNjqrZLfdBuXL2ycuw1S5sn393/9cmpNIPqZAJiK5cnBIWZRJyY17QtqudBy7XXnhhVCtmocTjGMlKMEd3MHHfMxiFvtb/KSTYNAg+PxzmDMnrCG0QiZymAKZiORqG9s4L+Ngv6sQVsjmzrVbltqu9FY/+lGMYjzDM/4Xv/12SE623fvDoEAmcpgCmYjkKpNMWs2CndWTQzo/NmqUfSCvc2fv5iaQQgq96MUHfMBKVvpbvGRJGDAAJk2C+fND/njx4vaieQUyEQUyEclDprOBVrNga+sz7YGwIOzfDx99BB07QunSHk9QuIu7MBie5En/i/frZ/8lR7BKpjNkIgpkIpKH3b/Np+IG2N/q3KA/M2UKZGZCz54eTkwOqU51etGLEYxgFav8LV6mDPTvD+PHw6+/hvzxSpW0QiYCCmQikoekDNtSoUjrdkF/5r33oEIFaN/eq1nJsQYzGCA6q2QDBkCJEmGtkimQiVgKZCKSq7IZC/irGpSt0zSo92/aZI8Ude8OhXXjrW9qUIObuZkRjGA1q/0tnpwMd9wBY8bAokUhfbRyZW1ZioACmYjkxnGoNusPvmmVQDFTPKiPpKXZM2TarvTfYAbj4ERnlWzQINsK48EHQ/pYpUqwYwfs3OnNtETyCwUyETmxJUsotX4XP7UuE/RH3nvPXpXUpImH85Ic1aQmN3ETwxnu/ypZ+fJ263LsWFiwIOiPBVpfrF/v0bxE8gkFMhE5sYOXRy9rVTWoty9ebPuPaXUsegYzmGyyo9O9f+BA+8TlAw8E/ZFAt36dI5N4p0AmIieWkcGGKoXZUze4VvvvvQcJCXDddR7PS06oFrW4kRsZznDWsMbf4mXL2lD28cfwww9BfSSwQqZzZBLvPA1kxpgOxpglxphlxph7cvj+QGPMImPMAmPMDGNMTS/nIyIhcByYNYtvWyWSYirk+fYDB+CDD+CSS6BiRR/mJyc0hCEc4EB0VskGDLDBLMhVMnXrF7E8C2TGmATgVeASoCHQzRjT8Ji3zQeaOY5zOpAOPO3VfEQkREuXwv/+x7TW+0khJc+3T58Oa9dquzIW1KY2PenJm7zp/1myUqXgzjth8mT47rs8316+vF1VVSCTeOflClkLYJnjOH86jrMPSAM6HfkGx3G+dBxn18GffgfoCmKRWJGRAcC01llBBbL33rMLIx07ejwvCcq93Es22TzMw/4X79fPJq0gVskSEmzPOgUyiXdeBrKqwF9H/Hz1wa+dSC9giofzEZFQzJpFVqUUltYjz0C2bRtMmABdu9r7KyX6alGLPvThbd5mCUv8LV6iBNx9N3zxBXz9dZ5v1/VJIt4GspwuvXNyfKMx1wPNgGdO8P1bjTHzjDHzMjMzXZyiiOTIcSAjg62tzwADFcj9DNnYsbBnj7YrY81QhpJEEvdxn//F+/a1Sev++/N8q7r1i3gbyFYD1Y/4eTVg7bFvMsa0BYYCVziOszengRzHedNxnGaO4zRLScl760REIvTHH7B2LWtb1QXyXiF77z2oXx9atPBjchKsClRgEIMYy1jmMc/f4sWLw+DB8OWX9kcuKldWIBPxMpDNBeoZY2obY4oAXYGJR77BGHMm8B9sGNvg4VxEJBQHz48tbW1PGeQWyJYssbtSN90EJqd1cYmqQQyiHOUYwhD/i996K1StalfJnBw3SAC7QrZ+PWRn+zg3kRjjWSBzHCcLuAOYCiwGxjiOs9AY87Ax5oqDb3sGKAGMNcb8ZIyZeILhRMRPs2ZBxYosq58A5B7I3n7bHszWdmVsKkUphjKUaUxjBjP8LZ6UBEOH2sQ+bdoJ31apEmRl2XtQReKVcXL5W0ssatasmTNvns9L7yLxxHGgRg34xz/49+jqvMZr7GJXjm/dvx+qV4dzzrG9QCU27WEPp3AKlajE93yPyfGIr0f27oVTTrH7kv/9b47LqOnpcO218PPPcPrp/k1NxA/GmB8cx2mW1/vUqV9EjrZ8OaxeDa1bk0lmrqtjkyfbraZevXycn4QsiSQe4iHmMpcJTPC3eNGicN998P338NlnOb6lwsFnRvTMlsQzBTIROdrB82O0akUmmbk+YTlihF34uOQSf6Ym4etBD07lVIYylCyy/C3esyfUqQP33pvjQbFAINugk8QSxxTIRORoGRmQkgKnnprrCtmaNXbB48YbITHR1xlKGBJJ5DEe4zd+4y3e8rd44cLwyCPw008wZsxx3w48PK9AJvFMgUxEDjt4fyWtWoExuQay996zix033+zzHCVsV3IlF3AB93M/29jmb/GuXe0BsXvvtYcPj1C2rH0wRFuWEs8UyETksBUrYNUqaN0agA1syDGQZWfbpytbtYK6df2dooTPYHie59nIRh7jMX+LFyoETzxhe9yNGHHct8qX1wqZxDcFMhE5bNYs+9qqFTvZyW525xjIZs2yf6727u3z/CRiTWlKT3oyjGH8wR/+Fr/kEjj/fHjoIdh19JO7FSpohUzimwKZiByWkWGXKho2JBP7p2NOh/pHjIDSpaFzZ5/nJ654jMcoTGHu4i5/CxtjV8nWrYNhw476VkqKVsgkvimQichhgfNjhQodCmTHrpBt2QLjxkH37lCsWDQmKZGqQhXu5m7GM55ZzPK3+HnnQceO8NRTsHnzoS9rhUzinQKZiFgrV9ozZK1aAfb8GBwfyD780F4krt5j+dsgBlGNagxkINn4fGfRY4/B9u02lB2kFTKJdwpkImIFzo8dPNCf0wqZ49jtyjPOgLPO8nuC4qbiFOcpnuJHfuR93ve3+GmnwfXXw0sv2f4p2BWybdtg3z5/pyISKxTIRMTKyIDkZGjUCCDHM2Tz5sH8+TrMX1B0oxtnczaDGczf/O1v8YcfhgMH7AF/Dvci07alxCsFMhGxMjIOnR8Du2WZRBIncdKht7zxBpx0EvToEaU5iqsMhhd5kXWs42Ee9rd4rVrQt69dcl24UN36Je4pkImI7T22fPmh82PAoaawgYuot2yBjz6yh/lLlYrWRMVt53AON3MzL/ACC1nob/H77oOSJeHuu7VCJnFPgUxEjjs/Bhx3j+UHH8Du3fB//+fz3MRzT/IkJSnJ7dyOg+Nf4XLlYOhQmDyZCn/8F9AKmcQvBTIRsYGsbFl72PqgI69Nchy7XdmiBZx5ZrQmKV5JIYUneIJZzOJDPvS3eL9+ULMmFZ67G9AKmcQvBTIRsefHLrjg0PkxOPrapNmzYfFi6NMnSvMTz/WmNy1owSAGsZWt/hVOSoLHHqP0L19ROOGAVsgkbimQicS71avtPUhHbFfC0VuWb7wBZcpAly5RmJ/4IoEEXuM1NrCB+7nf3+LdumGaNiWFTDL/l+VvbZEYoUAmEu+OuL8yYCc72cUuUkhhwwbbmb9nTyhePEpzFF80pSl96curvMp85vtXuFAhePZZUg6sY8Oc5f7VFYkhCmQi8W7WLLv8dfrph750ZFPYt9+G/ft1mD9ePMqjlKc8fejjbwf/1q2pUN4h8/ctsHGjf3VFYoQCmUi8y8iA88+HhIRDXwoEsnLZKfznP3Y3s0GD6ExP/FWGMjzDM3zP97zJm77WTjm7DhsOlLNNY0XijAKZSDxbuxaWLs3x/BjA8i/qsWKFVsfiTQ960IY23MVd/MVfvtWtUK80mYWrwGuvwaJFvtUViQUKZCLxLIf+Y3D4YvHJr1enQgW46iqf5yVRZTC8yZsc4AB96ONbb7KUFNixvxi7TyoPAwfafisicUKBTCSeZWTYtvtNmhz15UwyYWUNZn5anF69oEiR6ExPoqcOdXiMx5jMZD7iI19qBq5PyvznozB1Knz2mS91RWKBAplIPAv0Hzvi/BjYQJbwWn+M0XZlPOtHP87hHPrT/9A2tpcOXZ90aU+oX9+uku3b53ldkVigQCYSr9auhd9/hwsvPP5bu7biDO/FVVcZatSIwtwkJiSQwAhGsIMd9Ke/5/UOXTC+pTA8/7z97/OVVzyvKxILFMhE4lVGhn095vwYwM+jGpO9pQz9vf8zWGJcQxpyL/eSRhoTmehpraMuGL/0UrjkEvvEpdr3SxxQIBOJV19+afuPHXN+zHHgj5cupeQZf9CyZZTmJjHlbu7mNE6jD308vVbp0ApZIH89/zzs3An33edZTZFYoUAmEq9OcH4sIwN2/1qXxv1nYkxUZiYxpghFeJu3Wcc6BjDAszolS9oHSA5dMN6gAdx+OwwfDj/95FldkVigQCYSj1avhmXLcjw/9tJLQPlMmndb5v+8JGY1oxlDGcp7vMd4xntSwxi7SnbUDuUDD0ByMvzzn2qDIQWaAplIPDrB+bHly2HiRAdufZPKSWV9n5bEtvu4j6Y05VZuZR3rPKmRknJMICtbFh5/HGbPho/8ab8hEg0KZCLxKCPD/kF3xP2VYBukGwP0eZ0KVIjK1CR2FaYwIxnJTnbSi16eNIxNSYFNm475Yq9e0Lw5DBoE27e7XlMkFiiQicSjL7+EVq2g0OHfAnbuhLfeggs7b4Fqa0ghJYoTlFjVgAY8zdN8xmee3HWZnJxDIEtIgFdfhfXr7RamSAGkQCYSb1atgj//PO782MiRsHUrtP/nbwAKZHJCt3M7F3MxAxnIUpa6Ona5crB5cw7faN4cbr0VXn4ZFixwtaZILFAgE4k3OZwfcxx7mL9pUyh/rv0DVluWciKFKMQ7vENRitKDHmSR5drYycmwZQscOJDDNx97zLZquf12HfCXAkeBTCTeZGTYZYjGjQ99aepUWLQI+veHTGNPVGuFTHJTlaq8zut8z/c8wiOujVuunM1a27ad4JtPPQVffw0ffOBaTZFYoEAmEm9yOD/27LNQpQp07WrvsSxKUUpQIoqTlPwglVR60pNHeISZzHRlzORk+3rcObKAm26Cs8+GO++0e+wiBYQCmUg8WbHC/jhiu/LHH2HGDNvmqUgRG8hSSMGgrrCSt1d4hfrUpzvdWc/6iMcrV86+5niODOxfJF57zXaPVQd/KUAUyETiSeD82BEH+p991nZIv+02+/MNbND5MQlaCUowmtFsZSs96EE22RGNl+cKGcBZZ0HfvjaYzZ0bUT2RWKFAJhJPMjKgfHlo2BCAlSthzBj78Frp0vYtgRUykWCdzukMYxjTmMZTPBXRWHmukAU89hhUrAi33AL790dUUyQWKJCJxAvHsefHWrc+dH7sxRdtI9h//vPw2xTIJBy3cAuppHIf9/E1X4c9TlArZGD/BvHKK/Dzz/Y/ZJF8ToFMJF6sWGF7kB08P7Zli72zuWuGC1YaAAAgAElEQVRXqF798Nu0ZSnhMBje5E1qUYuudGUjG8Map0wZ+5eEPFfIAK66Cjp1ss1ily8Pq55IrFAgE4kXx5wf+89/bHf+f//78Ft2Hfw/rZBJOEpRijGMIZNMutOdA+TUTCx3CQk2lOW5QgY2ub3yiv1Qnz7qTSb5mgKZSLz48kuoUAFOPZW9e2HYMLj4YmjS5PBbMskE1INMwncWZ/EyL/MFX3Af4T0FecJu/TmpVg2eeMI209Pl45KPKZCJxAPHsb0tLrwQjGHUKFi3zrZyOtIGbFNYbVlKJG7lVm7hFp7gCcYxLuTP53ifZW769LG9yQYMCPGDIrFDgUwkHvz2G6xdC23bkp1tW100aQJt2x79Nq2QiVte5mXO5mxu5EYWsSikz4a0QgZ2y/LNN+3ByCP34EXyEQUykXgwY4Z9bdOGTz+FxYvtn1vmmN6vCmTilqIUZRzjOImTuJIr2UZOdyHlLOQVMoDTT7f/Ub/7LnzxRYgfFok+BTKReDB9OtSpg1OrNo88ArVr26crjxXYslQgEzdUpSpjGctylnM91wfdNDbkFbKABx6ABg2gd2/Yvj2MAUSiR4FMpKDLyrJPWLZpwxdfwLx5MHgwJCYe/9bAPZYlKen7NKVgOp/zeZEX+ZRPeZAHg/pMcrK9XDwrK8RiSUl2hWzNGm1dSr6jQCZS0P3wA2zbhtOmLY88Yh9Ku+GGnN+qeyzFC33py83czCM8wihG5fn+QLf+LVvCKHb22TaMDR+urUvJVxTIRAq6g+fHZie145tv4O67oWjRnN+qLv3iBYPhdV6nNa25mZvz7OQfdLf+E3noIbt1ecst2rqUfEOBTKSgmz4dzjiDR18pQ8WK0KvXid+qLv3ilSIUYRzjqEUtruRK/uCPE7436PssTyQpCd55B1avPr63i0iMUiATKch27YJvvuG7U29i+nS7k1Os2InfrhUy8VIyyXzKpzg4XM7lbCHnPcmIV8gAzjkHBg607TCmT49gIBF/KJCJFGTffAP79vHoH11JTob/+7/c365AJl6rRz0mMIE/+INruIb97D/uPRGvkAU8/DDUr2+XhbcF33ZDJBoUyEQKsunTmZ/YnMlzKvCvf0GJEid+6y52sZOdCmTiuQu4gOEMZyYzuY3bcDj6DkpXVsjALgcHnrrs3z/CwUS8pUAmUpDNmMFjZZ6mdGno1y/3twaawuoMmfihJz25n/t5h3cYytCjvle6tG2+H/EKGdity3vvhfffhzFjXBhQxBsKZCIF1ebNLPxhD+M2tqZfP/uHXG7UpV/89iAPciu38gRPMIxhh75uDJQt6+K1lEOHQosWds9+zRqXBhVxlwKZSEH15Zfcz0OULJ7FgAF5v11d+sVvBsNrvMZVXMUABvAhHx76Xtjd+nNSuDCMHAl798KNN0J2cDcGiPhJgUykgPrho98ZT2cGDjSHDknnRluWEg0JJPAhH9KKVvSkJ1OZCthzZK4FMoB69eDFF+0Tly+95OLAIu5QIBMpoO6dch7Jhbcz8M6EoN6vLUuJliSS+IRPaEQjOtOZOcyhXDkXtywDeveGK66Ae+6BX391eXCRyCiQiRRAX49bz+e7LuDuDgsoVSq4z2xgA0UoonssJSpKU5opTKECFehAB0je7O4KGdjDacOH2wOV3bvDnj0uFxAJnwKZSAHjODB0iENF1nHHfWWD/lwmmVSggu6xlKipTGVmMIOTOInp5UaTuemA+0UqVLBd/BcsUBd/iSkKZCIFzPTpMPv3Sgwt+TLFmzUM+nNqCiuxoDa1mclMEpO3s+vvBH7d97v7RS69FAYNgldegfHj3R9fJAwKZCIFiOPA0KEONQr9xa2XrbFbNEHawAYFMokJ9ajHwHI3AHDx5tRc770M2+OP21YYN98My5e7P75IiBTIRAqQiRNh7lzD/dkPUvTyi0P6rFbIJJY0TK4MwO7NxbiIi1jJSncLFCkCaWn2f3ftCvv2uTu+SIgUyEQKiOxsuO8+qJu8iZ68D+3ahfT5wBkykVgQaNXy3KZ32M52LuRCVrDC3SK1a8Nbb8GcObZ5rEgUKZCJFBAffQS//AIPlXuZxOZnQkrwq126x1JiTdmDz6OkbKnPNKaxla2cz/ksZam7ha65Bvr2hWefhcmT3R1bJAQKZCIFwO7dMGQInHl6Fl2XPQodOoT0+UCXfq2QSawItGvZvh2a0Ywv+ZK97OUCLmARi9wt9txz0KQJ9OwJq1a5O7ZIkBTIRAqAYcPsnyPPdfqKQs4BuOSSkD6/nvUAVKKSF9MTCVng7tVt2+xrE5qQQQYGQyta8RM/uVcsKclePL5vn10xU38yiQIFMpF8LjPTPjDWsSNc+Nf7dq+nRYuQxljHOkCBTGLHsYEMoCENmc1silOcC7mQOcxxr+App8B778HcudC/v3vjigRJgUwkn3voIdi1C55+yoHPP7eH+ROCuy4pIBDIKlLRiymKhKxoUfsg5PbtR3+9LnWZzWySSaYNbZjGNPeKXnWVvVZp+HAYMcK9cUWCoEAmko/99hu88Qbcdhs02PszrFsX8vkxOLxlqTNkEktKlz56hSygJjX5mq+pQx0u4zI+4iP3ij76KLRtC7ffDvPmuTeuSB4UyETysbvvhuLF4YEHsKtjAO3bhzzOOtZRjnIUoYi7ExSJwIkCGdhrlmYzm3/wD67jOl7gBXeKJiTYR5YrVoTOnWHjRnfGFcmDAplIPpWRYRvBDhlir+djyhQ44wyoXDnksdaxTtuVEnNKlTp+y/JIpSnN53xOZzozkIHczd04OJEXLl/eXqm0fj106wYHPLhTU+QYCmQi+VB2tr2Kr0YN+Oc/scsI33wT8tOVAetZrwP9EnNyWyELSCKJ0YymD314mqfpSU/24ULX/aZN4bXX7OWwuoRcfJAY7QmISOg++AB+/BFGjoRixYApM+zf4sM4PwZ2hewcznF3kiIRKl0a/gjiGssEEniVV6lCFe7jPlaykvGMpxzlIpvAzTfDzz/DCy9Aw4bQu3dk44nkwtMVMmNMB2PMEmPMMmPMPTl8/wJjzI/GmCxjzDVezkWkoNi6Fe66C845x+6mAHa7slQpOPfcsMbUlqXEory2LI9kMNzLvYxiFN/zPedwDktYEvkknnvOnsvs08eeExDxiGeBzBiTALwKXAI0BLoZYxoe87ZVwI3/3959x0dZZWEc/92ETgCRovQiICiiQEBQUQRsrGBZXbDg2kAXRF2VtaCuBVfEXlARC7ZVlLIgYEFBUYoEQxcQUMTQpEsL9e4fJ8EQAqTMzDuTeb58sjNJ3nnf87KSnLn33HOB/4YrDpHC5t//tt5jAwdCQgLgM9pddOgARYvm+Xxb2cp2tmvKUqJObqYss7uSK5nABDazmVa0YgITChZEkSK2CXm9elbkn5shO5F8COcIWUtgiff+Z+/9LuBD4KKsB3jvl3nv5wD7whiHSKExeza89JK9WW/WLOOL8+dDWlq+68fUFFaiVeYImc9jnf5pnMZ0plONapzHeQxmcMECOeoo+OQTe96pU96zRJFcCGdCVg34LcvnaRlfE5F88N5aIx19tLVK2i+z3UUB6sdACZlEn3LlbAHLtm15f21tajOZyXSgAz3oQS96FazYv149GD4cFi+GLl1gz578n0skB+FMyFwOX8vXemTnXA/n3Azn3Iy1a9cWMCyR2PTuu7aQ8oknbHek/T79FBo3hurV83XezKawqiGTaJPT9kl5ej3l+IRP6EMfXuZl2tGOVazKf0Bt29rKy88/t+2V8jp0J3IY4UzI0oAaWT6vDqzMz4m8969575O998mVKlUKSXAisWTTJlt5f+qpcO21Wb6xeTNMmpTv6UrQCJlEr7Jl7bEgM4RFKMIABjCUocxkJs1pzhSm5P+E3bvbqppXXoH+/fN/HpFswpmQpQD1nXN1nHPFgK7A6DBeT6TQyizkf/nljEL+TJ99ZlMnF110yNceyWpWk0ACFalY8EBFQihzhCy3Ky0P52/8jWlMoxSlaEtbXuGV/DeRffxxuOoq68r87rsFD06EMCZk3vs9wC3A58AC4CPv/Xzn3CPOuc4AzrkWzrk04HJgkHNufrjiEYlVmYX8N9+cpZA/0yefWFfxVvnvIbaGNVSiEonkbUNykXAr6JRldidxEimk0IEO9KQnV3M1W9iS9xMlJMCbb0K7dtarbHwINziXuBXWPmTe+3He+wbe++O8949lfO1B7/3ojOcp3vvq3vvS3vsK3vsTwxmPSKzZu9dmSCpUyFbID7B7N4wdCxdeaPvv5dNqVmu6UqJSKKYssytPecYwhn7040M+pDnNmcWsvJ+oWDHbXumEE6wdxqx8nEMkC22dJBLFnn8eUlLgxRdtdeUBJk+24rJOnQp0DSVkEq1COWWZVQIJ9KUvE5nINrbRilb5m8IsVw7GjbPHCy6AZctCG6jEFSVkIlFq6VK4/37Lt/72txwOGD3a3qWfe26BrrOGNVphKVEp1FOW2Z3JmcxiFu1oR0960oUubGRj3k5SrZrVcqanW3PmVQVYxSlxTQmZSBTy3qYqixa1Qn6XvYmM95aQtW8PSUn5vw5eI2QStZKS7L/9cPZhrUQlxjCGJ3iCkYykCU3y3t3/xBNtpGz1anuDtH59eIKVQk0JmUgUeuMNmDgRnnzyEO3FFi60IbTOnQt0nU1sYhe7lJBJVEpIgDJlQj9ledB1SOBf/IupTKU0pWlPe+7kTtJJz/1JWreGUaOscewFF8CWfCwWkLimhEwkyqxcCXfdZT0ob7zxEAeNGmWPF15YoGupKaxEu/zsZ5lfySSTSio96ckzPENLWjKXubk/Qfv28PHHkJpqtQY7doQvWCl0lJCJRBHvoWdP2LkTBg/O1nMsq+HDrUtsPrvzZ1JTWIl2ZctGduvIUpRiIAMZy1h+53eSSeZxHmcPudwqqVMn6002aRJcdhnsKsB2TRJXlJCJRJGPP7bBr0cfta3zcrR8OcyYAZdeWuDrKSGTaFeuXPinLHPSkY7MZS6d6cx93MepnMoc5uTuxVdcAa++anVlXbsqKZNcUUImEiVWrLDmry1awO23H+bAkSPt8ZJLCnzNzClLJWQSrSI5ZZldJSrxMR8zjGGkkUZzmvMQD+Vuk/IePaxvzciRthm5kjI5AiVkIlFg3z7bo3LnTnjvPShS5DAHjxgBJ50E9esX+LqrWU1RilKe8kc+WCQAkZ6yzMlf+Ss/8iNd6MLDPEwyyUxn+pFfeOut8MIL8L//We8aJWVyGErIRKLACy/Al1/Cs89CgwaHOXDNGvj225BMV4IlZJWpjCN7Xw2R6BDUlGV2FajAe7zHaEazgQ20ohU96ckmNh3+hb17295no0ZZTdnOnZEJWGKOEjKRgM2dC/fcYx0sunc/wsGjRlnlf4gSshWsoDoFWxggEk7RMEKWVSc6sYAF3MqtDGIQDWnIB3xw+C7/vXrBwIG296ySMjkEJWQiAUpPh6uugqOOgtdfz6EBbHYjRsBxx9mUZQikkaaETKJauXL27ySaZvvKUIbneI4UUqhBDa7kSs7jPH7ip0O/qGdPeOUVGDMGLroItm2LXMASE5SQiQTo/vtthOzNN6FSpSMcvH49fPWVjY4dMXPLnTTSqEa1kJxLJBzCtZ9lKDSjGdOYxku8xPd8T2Ma04c+/MEhgr35Zuv6PH68dfTfmMdtmqRQU0ImEpCvvoKnn7Y3zh075uIFI0bAnj22jD4E/uAPtrBFI2QS1cqWtcdomrbMKpFEetGLn/iJbnTjKZ6iAQ0YwhD2se/gF1x/PQwdCikp1v159eqIxyzRSQmZSABWrbKpyoYNbXukXBk61JqTNW0akhhWsAJACZlEtXBvMB4qx3AMb/AG05lOHepwHdfRmtZMZvLBB192GYwdC0uWQJs2sGxZxOOV6KOETCTC9uyxvpFbtsCwYVCqVC5etGaNbW7ZtWtIpysBTVlKVIvmKcuctKAFk5nMO7zDb/zGGZzBxVzMAhYceOA559jS6nXr4PTT4ccfgwlYooYSMpEIe+AB+OYba+R94om5fNGwYdasrEuXkMWhETKJBdE+ZZmTBBLoRjcWs5h+9GMCE2hMY7rTff+/O8A2JJ80yf5tn366/WCQuKWETCSCPvkE+ve3Jt7duuXhhUOHWvbWuHHIYskcIatK1ZCdUyTUYmXKMielKU1f+rKUpfSmN2/zNvWpz33cx2Yybuikk2DqVKhSxUbN3n8/2KAlMErIRCLkl1/gmmusBOz55/PwwrQ0awYbwtExsISsIhUpQYmQnlcklGJtyjInlajEczzHIhZxCZfwOI9Tl7o8wzPsYAfUrg2TJ9so2dVXw2OPWb9BiStKyEQiYOdOuPxy+xk7bBiUyEsO9NFH9hjihExNYSUWxOKU5aHUoQ7v8z6ppNKc5tzJndSlLs/yLNvLF4fPPrOE7P77bRh99+6gQ5YIUkImEmbe2+4pP/wA77wDdevm8QTvvms7jh92T6W8Uw8yiQXFi9tHYUjIMjWlKV/wBV/zNY1oxB3cQR3q8FTxF9n2zitWaPr669YPZ8OGoMOVCFFCJhJmL74IgwfDvffa9kh5MncuzJqVx4Kz3NEImcSKaNnPMtTO4iwmMIFJTKIJTehDH+q4ugx4JIn0Ia9awX/LljB/ftChSgQoIRMJo88+g3/+Ey6+GPr1y8cJ3n0XihQJWTPYTOmks5a1SsgkJkTbfpah1oY2jGc8k5lMU5pyN3dT/e99GfL1tezbthVatbJ9bKVQU0ImEiY//mhlX02aWF6VkNd/bXv32oqrCy7Ixb5KebOSlYB6kElsKFeucCdkmU7jND7nc6YylZa05LrWr9FgxmZ+bVTK3tU9+qiK/QsxJWQiYbBuHXTqBCVLwujRkJSUj5NMmAArV9rSzBBTDzKJJYV1yvJQWtGKcYxjDnM4o1oXGk/awLvdgAcfZMNl7eMjO41DSshEQmzXLtv/e8UKm2WoUSOfJ3rnHftNdOGFIY0P/uxBpoRMYkFhn7I8lJM4iSEMYWGJZcx7uw/3PlOcsqMmkpZchUmzXsh5r0yJWUrIRELIe1ut/u238NZbcOqp+TzRH3/YZuJ/+1see2TkjrZNkliSlARbtwYdRXCqUY0n3ADu/efvDPvmFhJ37KJlq9t4YPCxPO2fYgNaiVkYKCETCaF774W334aHHrL9KvPtgw9g+3a44YZQhXaAFaygDGUoS9mwnF8klMqUie+ELFNZytL19BepOHM5G89qwmM91lLp732ov60q13M9P/BD0CFKASghEwmRp5+GJ56Af/wDHnywgCcbPNhWA7RsGZLYsksjTdOVEjPifYQsu6KVqlJlXCo88gjd3nPMSy7JwpkfkEwyrWjFu7xLOulBhyl5pIRMJATefhvuusu68b/4IjhXgJOlploX2e7dC3iiQ1NTWIklSUmQng579gQdSRRJTIQHHsB9+SVV/ijF5FP3Mumpzmzat4FruIbqVOd2bmcOc4KOVHJJCZlIAY0ZYzOLHTpYe4vExAKecPBgqxu76qqQxJcTNYWVWJK5SlmjZDlo1w7mzMF16kSbPqNZcE4Nvkl7n3a042Ve5mROpgUteIVX2MSmoKOVw1BCJlIA331no2JNm1oNfvHiBTzhtm3We+zyy6F8+ZDEmN0e9rCKVRohk5hRpow9KiE7hAoVbJPcN97Aff89Zza5hY8+uoyVrOQ5nmMnO+lJT6pQhau5molM1ArNKKSETCSfpk2zreZq1oRx4/78pVEgH30EW7bYdGWYrGENe9mrETKJGRohywXn4PrrYeZMqF8funSh4uX/4Lbfr2A2s5nBDK7nesYwhna0ow51uId7mMvcoCOXDErIRPJh6lQ491w45hj46qsQNdL33grQTjgBzjgjBCfM2a/8CkBNaobtGiKhpIQsD+rXt6H7xx6zrtQnnID7cCjNfTMGMpBVrOJ93qcxjXmKp2hCE07iJPrTf//PBgmGEjKRPJoyxZKxY4+Fr7+G6qEaaJoyxd7d9u4dtmJ+gKUsBeA4jgvbNURCKTMh27Il2DhiRtGicN99tkDouOOsB8+ll8Lq1ZSkJFdyJWMZyypW8RIvUZay3Mu91KY2bWjDK7zCGtYEfRdxRwmZSB5MngznnQdVqsDEiVAtlGVYL7xgnfm7dQvhSQ+2hCU4HLWpHdbriISKasjy6cQT7YfWgAHw6afQqBEMGgT7rH6sEpXoRS8mM5mf+Zl+9GM96+lJT6pSlba05SVe2r/3rYSXEjKRXPruO0vGqla1kbGQJmNpaTB8ONx4I5QuHcITH2wpS6lBDYpT0BUIIpGhKcsCKFIE+vSB2bPhlFPg5pvh9NNh1qwDDqtDHfrSl/nMZzazuZ/7WctaetObalTjDM7gWZ5lOcsDupHCTwmZSC6MHWvJWPXqloxVrRriC7z6qr1r7dUrxCc+2FKWarpSYooSshA4/niYMMH2yF26FJo3hzvuOGge2OFoQhMe5mHmM58f+ZFHeZStbOUO7qAWtWhBCx7hEWYyE48P6IYKHyVkIkfw1ltw0UXQsCF8841NV4bUjh02jdC5M9SpE+KTH0wJmcQaJWQh4pyVRCxaZCu5n3vOpjHfe2//NGZ2jWjE/dzPLGaxmMX0pz9FKcpDPEQzmlGLWvSkJ5/yqXYHKCAlZCKH4D08/ritJD/7bBsZO+aYMFzorbdg3Tq4/fYwnPxAW9jCWtZSj3phv5ZIqKioP8TKl7dR+SlT7B1mt25w2mnWy+cw6lGPu7mbKUxhNat5kzdJJpm3eZuOdKQiFbmUSxnMYK3YzAclZCI52LcPbrvNFipdcYVNWYakz1h2e/bAk09Cq1Zw1llhuMCBtMJSYlGRIrZ5hUbIQqxVK/j+exgyBJYvh9at4eqrrab1CCpTmeu4jhGMYD3rGcc4ruEaUkihBz2oTW0a0pDbuI2xjGUb28J/PzFOCZlINtu3Q9eu1hLsn/+00fxixcJ0saFDYdkyuOeesLa6yKSETGKVNhgPk4QE+Pvf4aef7B3osGHQoAHcey9syt1WSyUowQVcwMu8zHKWM5/5PMuz1KEOgxnMhVzI0RxNe9ozgAHMYpZqz3KghEwki99+gzZt7GfSk0/CM8/Yz6uw8B7697dGsJ06hekiB1JCJrFKCVmYJSVZM9mFC+GSS+xnU9261jJjx45cn8bhOIETuJ3b+ZRP2cAGxjOeW7mVtazlbu6mKU2pQhWu5Epe4zV+4iclaCghE9lv8mRITobFi63B9V13hfmCY8fCvHk2Oha2rO9AS1lKRSpSlrIRuZ5IqCQlqYYsImrXtv10Z860Kc2774Z69WDwYNi9O8+nK0EJOtCBJ3mSOcxhBSsYwhDa056v+ZqbuInjOZ5qVONKrmQQg1jEorhM0Jz3sXXTycnJfsaMGUGHIYXM669Dz55Qq5YlY40ahfmC3kPLlrB2rWWARYuG+YKmPe3ZxjamcfjiXZFoc9pp1qJv/PigI4kz33xj05dTp1qydt99NsUZgjoOj2cJS/g6489EJrKKVQAcy7G0zfjThjY0pCEJMTqG5Jz7wXuffKTjYvPuREJk1y7bqah7d1tJOX16BJIxgFGjYMYMePDBiCVjoJYXErs0ZRmQs86y6YMxY2zT3h49rMZs0CD7AVoADkd96tOd7rzP+6xgBT/xE6/xGu1oxzd8w83czImcSAUq0JGO9KMfE5jAVgrffwwaIZO4tWSJFe//8APceaeVTBQpEoEL791rHbN37YL58yN0UdjFLkpSkr705REeicg1RULl0kttMHnu3KAjiWPew2efwcMP2+rMGjXsh+cNN/zZmySUl8sYQZvMZKZk/JnPfAASSOBkTua0LH9qUQtH+BdH5VVuR8gi85tAJMp88AHcdJPlQv/7nzV+jZihQ6127MMPI5aMASxjGfvYpxEyiUkaIYsCzsEFF8D559vccb9+1j/x4Ydtl5HevaFy5dBdLmMErT71uZZrAdjIRr7n+/0J2hCGMJCBABzDMSSTTIssfypRKWTxhJtGyCSubNsGt94Kb75p27n9979Qs2YEA9i921ZVli4NqakRK+YH+JRP6UhHvuVbzuCMiF1XJBR69bL3MuvWBR2JHGDaNFuSPnKk1ZVde601cYxI7QfsYQ/zmMdkJjOd6aSQwkIW7l8UUJOa+5Oz5Iw/5SgXkdgyaYRMJJvp060WddEiuP9++Pe/IzpAZV55xeZKP/kkoskY/NnyQl36JRZphCxKtWoFw4dbH7Onn7Yms4MGQbt2lkV37hzWH7RFKMIpGX96YXsBb2ELqaSSQgozmEEKKQxn+P7X1Kc+zWhGU5rSkY6cxElhiy8vlJBJoZeebsnXU0/ZpuDjx0P79gEEsm6dBXLuufCXv0T88ktZSmlKcwzh2P9JJLySkmDnThtkjuA6GMmtzEL/fv3gjTfszedf/wrVq1t9SPfuYdp77mBlKMNZGX8yrWc9P/DD/iRtGtMYylCKUjRqEjKtspRCbepUq58fMMD2pJw3L6BkDOCBB6yR0rPPRqQrf3aLWUxd6kZl0avIkWTWjG/TDjzRrVIl6634889WoNuokf3sq1EDrrrKVmwGUCpVgQqcy7n0pS8jGckylrGe9VzHdRGP5VCUkEmhtHWrLf45/XRrMv3FF9bXsFxkSwf+NHs2vPaaDeGfcEIgIcxkJidzciDXFimozL1k1Rw2RiQm2mqpL76w7v//+Ie1zjjjDGjYEB5/PFd7ZobT0RxNecoHGkNWSsikUPHeFi82bGjbHt10k42KnXNOgEHt22eJWPny8NBDgYSwhjWsZCXNaBbI9UUKKnOETHVkMej44+H552HFCltRdeyx1mC2Zk047zxb9p6H7ZkKKyVkUmjMm2fNXa+4wkoVpkyxMobMd9aBGTjQhumfecaSsgDMZCaAEjKJWUrICoGkJLjuOuv+v2SJTWUuWgRXXmlJ2l49CgMAAA/jSURBVE03waRJ1qsxDikhk5i3YYO1wjnlFGsa+eqrtqKydeugIwOWLbNtR84/H7p1CyyMVFIBOIVTAotBpCCUkBUyxx1n/ct+/hkmTICLL4b33rOdAWrUsJ5m335rMwxxQgmZxKxt2+Cxx6BuXXjhBVvE89NP9iYrMTHo6LD50x49rIB/0KBACvkzpZJKPepFvP+OSKiohqyQSkiwqY2334Y1a6zmpHVr22D4zDMtObv1Vvjuu0KfnCkhk5iza5fNAh53nPUTO+ssmDPHpicrVAg6uixeftl6bDzxRIS7zx4slVSa0jTQGEQKQiNkcSApCbp0sb5mv/9utWWnnmoLotq0seTspptscUAhrDlTQiYxY+dOe9PUqBHccovViU6ZYvt0N24cdHTZzJoFd9wBHTvCzTcHGspGNvILv6h+TGKaErI4U6aMbTY8YgSsXWvbqrRubY+dOtm7706dLFlbsSLoaENCjWEl6m3dav/mnn4aVq6E5s3h009tcU6As4CHtnWrvcurWNG6Vke4I392KuiXwkAJWRwrU8ZWa11xhb0znzTJdjv55BMbLQNo1szeAJ9zju0eUKxYsDHng0bIJGqtXWs1n7VqWU+x44+3ljYpKVYjH5XJmPc2pL5kCbz/vjVJDFhmQqYpS4llpUvboxKyOFe8uCVdL7xgCwLmzYP+/aFkSfjPf6yGpXx5S86eecbqWWJkz26NkEnUmT7dasQ+/NDqxTp3toWKrVoFHVku/Oc/NqTerx+0bRt0NIDVj9WgBpUIPjkUya8iRex3ror6ZT/n4MQT7ePuu2HzZvj6a6vd/fJLeycPULmybdHStq0tFDj++Kh8R6+ETKLCjh3w8cfw0ks2ApaUZKsme/YMrLF93n38sa0yuOoqa3oYJVTQL4WFNhiXwypXznYHuOgi+/y33+CrryxB++orWyQAlqC1aWPJWceOUK9ecDFnoYRMAuM9TJtmq50//NDe3DRsaElZt25QtmzQEebBd9/BNdfAaafZyoMoefe1la0sYhFd6Rp0KCIFpoRM8qRGDbj2WvvwHhYvtvqzzI/hw60mrU+foCMFlJBJAH791Wb1hgyxvmElS8Jll9m/mbPPjppcJve+/97eZdWsCSNHQokSQUe0XwopeLwK+qVQUEIm+eYcNGhgHzfeaF/79Vf7BRQllJBJRPz8MwwbZh8pKfa1M8+0af/LLoux0bCsfvjBlntWrmzdpitXDjqiA4xiFMUpztmcHXQoIgVWpoxqyCSEatUKOoIDKCGTsPDeWnGNHWuDRqm2cw/JydYn9bLLrMN+TJswAS65BI4+2p5XqxZ0RAfweEYykvM4jySSgg5HpMCSkqy0QaQwUkImIbNli9VNjh0L48ZZzzCAli1hwABLwurUCTbGkHnvPbj+ehv+Hjcu8E78OUklleUs52EeDjoUkZBISio0PUBFDqKETPJtxw7rlD9hgn2kpMDevTb9eO658Je/wAUXwDHHBB1pCO3eDQ88YMN8Z59tXaSPOiroqHI0ghEkkkgnOgUdikhIqIZMCjMlZJJra9fC1KmWhE2ZYrXsu3bZRt4tW8I990CHDnD66VC0aNDRhsHy5dYpesoUa/76/PPWpDBKjWAEbWlLBaJpg0+R/CtTRgmZFF5KyCRHmzZZDVhqKsycacnX4sX2vaJFbZeK3r2hXTtr51KmTLDxhtW+fbZ30z332PMPPrA91qLYAhawkIX0pnfQoYiETFKSivql8FJCFuf27bPeefPmWeKV+fHLL38eU7WqFePfeKO12WrePKpWCodXSgrcdpsNDbZrB4MGRU0TwcMZwQgALubigCMRCZ2kJBuV37UrJrcqFDksJWRxYscOG+FauNA+Fiywx0WL7HuZ6tWz5Kt7d2ja1D4KVQ1Ybs2dCw8+CP/7n20S/vbb1q02BpqkbWc7b/AGrWlNVaoGHY5IyGRuML5tmxIyKXyUkBUSO3dCWhosW2ajW5mPmc9XrfrzWOes/UqjRlaX3rChPW/SxHaeiFt799oS0RdftH3Qypa13c1vvz2mGqU9xEP8wi+8xVtBhyISUpmlEVu32v7RIoWJErIo5r3VS6xdC+vWwe+/W2K1YoW1lMj6uG7dga9NTLRdI+rUgfPPt8cGDSz5ql8fSpUK5p6izp49ViA3dCh89BGsWQPVq8Njj1nhfoXYKohPJZWneZrudOcszgo6HJGQyhwhUx2ZFEZhTcicc+cDzwOJwOve+/7Zvl8ceAdoDqwHunjvl4UzpkjzHtLT4Y8/rKHh5s0HP9+0yRKqzMQr83HdOquVyM45m0asWtXaX7VqZc+rV7fEq04de15E6fbBdu+GOXNsE82vv7aRsE2bbLXkhRfaxuCdOsXkX95udnMjN1KZygxgQNDhiIRcZkKmlZZSGIXtt45zLhEYCJwDpAEpzrnR3vsfsxx2A7DRe1/POdcVeALoEq6YcuOPP2ybnx07YPt2+8h8nv3xUF/7448Dk67du4983fLlrVSpUiWoXRtatPjz86yPVavCscfGZL4QOenptlJh+fI/P5YssZULCxbY/C5YZ/1LL7WmaeefH9Pztd/wDbdwC/OYx3CGcxTR2RtNpCCUkElhFs5f6y2BJd77nwGccx8CFwFZE7KLgIcyng8DXnLOOe+9D2NchzXs46XccONxRzwuMWEvJUrsonjRPZQovovixXZTovhuihfbTemS6VSqkk7puukklUonqeQOSpdKp3TJdJJK7bCvlUqndIk/v14kcR94yPgfk/WvYRP4jbBiMawAnPf4A449MD6X7a8w+19p9u8fcK3DfQ9wPvOS/rDfz/p6v/97hz/3Qefa50nYuZvE9N0k7NxDYvouEnfuISF9N4npuyiydSfFNm6zj03b7XHzdrLbVv1oNjWuzub27Vjfoi5rW9dje40KGUX6u4FPDrynXMjtseE4J1jx/mIWM5vZTGQitajFCEZwCZfk+hwisSQzIfvqK3vTK1JQJ55oZTzRIJwJWTXgtyyfpwGnHuoY7/0e59xmoAJwQEWUc64H0AOgZpi3qKnshjKCFEqyg1JspxTb9z/P+lh03x44+Pe+hNE+BzuLQ3oJe9xZHLYkwcbysLE6bGpsz3+vDMtr2sdvNSCtOuwqvgHYAMwJ+jZCqhSlaEADHuIh+tCHUqg4UAqvKlXsPdR//hN0JFJYPPEE/OtfQUdhwpmQ5dQfIPvb/9wcg/f+NeA1gOTk5LCOnrX569Wsbp6tz5Rz7AJ2AZszPt//LRw++11kbY2Q7diD2iZk+9y7jOMO8f3s58567GHjAJxLOOz3DxX34T7ff/283NcRjj0gTufYV6IYvngx60ib7dhiOI4BctOZw+X4n1v+jwvXsbk9rhjFOJZjSSDhyAeLFAJVq9rK8Y0bg45ECosqVYKO4E/hTMjSgBpZPq8OrDzEMWnOuSJAOWwYIzDlytWk3MnRt1G0iIhYy55atYKOQiT0wvnWOgWo75yr45wrBnQFRmc7ZjTw94znlwETgqwfExEREQlC2EbIMmrCbgE+x9pevOm9n++cewSY4b0fDbwBvOucW4KNjEX3BoEiIiIiYRDW5gne+3HAuGxfezDL83Tg8nDGICIiIhLtVA0sIiIiEjAlZCIiIiIBU0ImIiIiEjAXa4sanXNrgV/DfJmKZGtOG2fi+f7j+d4hvu9f9x6/4vn+4/neITL3X8t7X+lIB8VcQhYJzrkZ3vvkoOMISjzffzzfO8T3/eve4/PeIb7vP57vHaLr/jVlKSIiIhIwJWQiIiIiAVNClrPXgg4gYPF8//F87xDf9697j1/xfP/xfO8QRfevGjIRERGRgGmETERERCRgSsiOwDl3l3POO+cqBh1LJDnnHnXOzXHOzXLOfeGcqxp0TJHinHvSObcw4/5HOueOCjqmSHHOXe6cm++c2+eci4qVR5HgnDvfObfIObfEOXdP0PFEinPuTefc7865eUHHEmnOuRrOuYnOuQUZ/83fFnRMkeScK+Gcm+6cm51x/w8HHVOkOecSnXMznXNjgo4FlJAdlnOuBnAOsDzoWALwpPe+iff+FGAM8OCRXlCIjAcae++bAD8B9wYcTyTNAy4FJgUdSKQ45xKBgcAFwAnAFc65E4KNKmKGAOcHHURA9gB3eu8bAa2AXnH0/zvATqCd9/5k4BTgfOdcq4BjirTbgAVBB5FJCdnhPQv8C4i7Qjvv/R9ZPi1NHP0deO+/8N7vyfh0GlA9yHgiyXu/wHu/KOg4IqwlsMR7/7P3fhfwIXBRwDFFhPd+ErAh6DiC4L1f5b1PzXi+BfvFXC3YqCLHm60ZnxbN+Iibn/POuerAX4DXg44lkxKyQ3DOdQZWeO9nBx1LUJxzjznnfgOuIr5GyLK6Hvg06CAkrKoBv2X5PI04+sUs4JyrDTQFvg82ksjKmLKbBfwOjPfex9P9P4cNuOwLOpBMRYIOIEjOuS+BY3P4Vl/gPuDcyEYUWYe7f+/9KO99X6Cvc+5e4Bbg3xENMIyOdO8Zx/TFpjXej2Rs4Zabe48zLoevxc1IQbxzziUBw4Hbs80MFHre+73AKRl1siOdc42994W+ntA5dyHwu/f+B+dc26DjyRTXCZn3vkNOX3fOnQTUAWY758CmrFKdcy2996sjGGJYHer+c/BfYCyFKCE70r075/4OXAi094WsN0we/n+PF2lAjSyfVwdWBhSLRJBzriiWjL3vvR8RdDxB8d5vcs59jdUTFvqEDDgd6Oyc6wiUAMo6597z3l8dZFCassyB936u976y976297429gO7WWFKxo7EOVc/y6edgYVBxRJpzrnzgbuBzt777UHHI2GXAtR3ztVxzhUDugKjA45JwszZu+03gAXe+2eCjifSnHOVMleQO+dKAh2Ik5/z3vt7vffVM36/dwUmBJ2MgRIyObT+zrl5zrk52NRtPC0JfwkoA4zPaPvxatABRYpz7hLnXBrQGhjrnPs86JjCLWMBxy3A51hh90fe+/nBRhUZzrkPgKnA8c65NOfcDUHHFEGnA92Adhn/zmdljJjEiyrAxIyf8SlYDVlUtH+IV+rULyIiIhIwjZCJiIiIBEwJmYiIiEjAlJCJiIiIBEwJmYiIiEjAlJCJiIiIBEwJmYiIiEjAlJCJiIiIBEwJmYjELedcC+fcHOdcCedcaefcfOdc46DjEpH4o8awIhLXnHP9sP3sSgJp3vvHAw5JROKQEjIRiWsZ+1emAOnAad77vQGHJCJxSFOWIhLvjgaSsP1LSwQci4jEKY2QiUhcc86NBj4E6gBVvPe3BBySiMShIkEHICISFOfcNcAe7/1/nXOJwBTnXDvv/YSgYxOR+KIRMhEREZGAqYZMREREJGBKyEREREQCpoRMREREJGBKyEREREQCpoRMREREJGBKyEREREQCpoRMREREJGBKyEREREQC9n+qtkGwdudyuQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np \n",
    "import matplotlib.pyplot as plt \n",
    "from scipy.stats import genextreme \n",
    "\n",
    "plt.figure(figsize=(10,8)) \n",
    "x = np.linspace(-4, 4, 201) \n",
    "#Please notice that in genextreme the shape parameter (gamma) equals to -gamma in our equations here\n",
    "plt.plot(x, genextreme.pdf(x, -0.6), color='#00FF00', label='Fréchet (\\u03B3=0.6) : \"Heavy tails\" distributions\\\n",
    " (Cauchy, Student, Pareto..)') \n",
    "plt.plot(x, genextreme.pdf(x, 0.0), 'r', label='Gumbel (\\u03B3=0): \"Light tails\" distributions\\\n",
    " (Gaussian, Log-normal, Gamma, Exponential)') \n",
    "plt.plot(x, genextreme.pdf(x, 0.6), 'b', label='Weibull (\\u03B3=-0.6): \"Finite tails\" distributions\\\n",
    " (Uniform, Beta)') \n",
    "plt.ylim(-0.01, 0.51) \n",
    "plt.legend(loc='upper left') \n",
    "plt.xlabel('x') \n",
    "plt.ylabel('Density') \n",
    "plt.title('Generalized Extreme Value Distributions') \n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use again our dataset of Nidd River, we will fit a GEV on it and estimate extreme quantiles.\n",
    "In pratice we will estimate the three parameters $\\gamma$, $a_{n}$ and $b_{n}$ such that:\n",
    "$$P(X_{n,n} \\leq x) \\simeq H_{\\gamma}(\\frac{x-a_{n}}{b_{n}})$$\n",
    "\n",
    "So first we need to turn our trimestrial data, that is our observations of $X$ into observations of X_{n,n}, lets say yearly maximums, we will create a vector of observations by taking the maximum flow for each year:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i in range(0, int(df.shape[0]/4)):\n",
    "    if i == 0:\n",
    "        df_max = df.iloc[4*i:4*i+4].max()\n",
    "    else:\n",
    "        df_max = np.append(df_max, df.iloc[4*i:4*i+4].max())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([ 9., 10.,  3.,  7.,  4.,  1.,  0.,  2.,  1.,  1.]),\n",
       " array([ 76.74 ,  99.641, 122.542, 145.443, 168.344, 191.245, 214.146,\n",
       "        237.047, 259.948, 282.849, 305.75 ]),\n",
       " <a list of 10 Patch objects>)"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAADD5JREFUeJzt3W+MZXddx/H3x93SUoq2tQOpbddpTUNsjKHNhFQxxLQQoTUuJn2wJGg1JJto0GI0ZomJ4LNilKgJgaxQrdq0aGlCY/FPU9oQE13c/qG0rLUrVCis3SWEAj6wVL4+mLPsOOzs7txzZm7nO+9XMrn3nnt2zm9+c/Lec8+9d26qCknS1vd98x6AJGkaBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhM7N3NjF110US0uLm7mJiVpy3v44Ye/WlULp1tvU4O+uLjIwYMHN3OTkrTlJfnPM1nPUy6S1IRBl6QmDLokNWHQJakJgy5JTZw26EluS3I0yRMrll2Y5P4kTw+XF2zsMCVJp3MmR+h/Drx51bJ9wANVdSXwwHBbkjRHpw16VX0K+NqqxbuB24frtwNvnXhckqR1mvUc+qur6gjAcPmq6YYkSZrFhr9TNMleYC/Arl27Nnpzk1vcd9/ctv3MrTfObduStp5Zj9CfS3IxwHB5dK0Vq2p/VS1V1dLCwmn/FIEkaUazBv1e4Obh+s3Ax6cZjiRpVmfyssU7gX8GXpPk2STvAG4F3pTkaeBNw21J0hyd9hx6Vb1tjbuun3gskqQRfKeoJDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUxIZ/puhU5vnZnpK0FXiELklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJamJUUFP8htJnkzyRJI7k5wz1cAkSeszc9CTXAL8OrBUVT8G7AD2TDUwSdL6jD3lshN4eZKdwLnAV8YPSZI0i5mDXlVfBv4A+CJwBHi+qv5xqoFJktZnzCmXC4DdwOXADwGvSPL2k6y3N8nBJAePHTs2+0glSac05pTLG4EvVNWxqvo2cA/wk6tXqqr9VbVUVUsLCwsjNidJOpUxQf8icG2Sc5MEuB44NM2wJEnrNeYc+gHgbuAR4LPD99o/0bgkSeu0c8w/rqr3AO+ZaCySpBF8p6gkNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJkZ9wIV6Wtx339y2/cytN85t29JW5xG6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCZGBT3J+UnuTvJvSQ4l+YmpBiZJWp+xn1j0x8DfV9VNSV4GnDvBmCRJM5g56Em+H3gD8EsAVfUC8MI0w5IkrdeYUy5XAMeAP0vyaJIPJ3nFROOSJK3TmKDvBK4BPlhVVwP/DexbvVKSvUkOJjl47NixEZuTJJ3KmKA/CzxbVQeG23ezHPj/p6r2V9VSVS0tLCyM2Jwk6VRmDnpV/RfwpSSvGRZdD3xuklFJktZt7Ktcfg24Y3iFy+eBXx4/JEnSLEYFvaoeA5YmGoskaQTfKSpJTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1IRBl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqYnRQU+yI8mjSf52igFJkmYzxRH6LcChCb6PJGmEUUFPcilwI/DhaYYjSZrV2CP0PwJ+G/jOBGORJI2wc9Z/mORngaNV9XCSnz7FenuBvQC7du2adXPaJhb33TeX7T5z641z2a40pTFH6K8Hfi7JM8BdwHVJ/mr1SlW1v6qWqmppYWFhxOYkSacyc9Cr6t1VdWlVLQJ7gE9W1dsnG5kkaV18HbokNTHzOfSVquoh4KEpvpckaTYeoUtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqYpIPuNDGmNcHJkvamjxCl6QmDLokNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktTEzEFPclmSB5McSvJkklumHJgkaX3GfGLRi8BvVtUjSV4JPJzk/qr63ERjkyStw8xH6FV1pKoeGa5/EzgEXDLVwCRJ6zPJZ4omWQSuBg6c5L69wF6AXbt2TbE5SROY52fWPnPrjXPbdmejnxRNch7wMeBdVfWN1fdX1f6qWqqqpYWFhbGbkyStYVTQk5zFcszvqKp7phmSJGkWY17lEuAjwKGqev90Q5IkzWLMEfrrgV8Arkvy2PB1w0TjkiSt08xPilbVPwGZcCySpBF8p6gkNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktTEJB8SLW11fmDy5prnfM/DZv2OPUKXpCYMuiQ1YdAlqQmDLklNGHRJasKgS1ITBl2SmjDoktSEQZekJgy6JDVh0CWpCYMuSU0YdElqwqBLUhMGXZKaMOiS1MSooCd5c5KnkhxOsm+qQUmS1m/moCfZAXwAeAtwFfC2JFdNNTBJ0vqMOUJ/HXC4qj5fVS8AdwG7pxmWJGm9xgT9EuBLK24/OyyTJM3BmA+JzkmW1feslOwF9g43v5XkqRHbnJeLgK/OexAvAc7DCZPNRd43xXeZG/eJE9aciwl+xz98JiuNCfqzwGUrbl8KfGX1SlW1H9g/Yjtzl+RgVS3Nexzz5jyc4Fwscx5OeCnMxZhTLv8KXJnk8iQvA/YA904zLEnSes18hF5VLyZ5J/APwA7gtqp6crKRSZLWZcwpF6rqE8AnJhrLS9mWPmU0IefhBOdimfNwwtznIlXf8zymJGkL8q3/ktSEQQeS3JbkaJInViy7MMn9SZ4eLi8YlifJnwx/7uDxJNfMb+TTWmMe3pvky0keG75uWHHfu4d5eCrJz8xn1NNLclmSB5McSvJkkluG5dtxn1hrLrbVfpHknCSfTvKZYR5+b1h+eZIDwz7x0eEFIiQ5e7h9eLh/cVMGWlXb/gt4A3AN8MSKZb8P7Buu7wPeN1y/Afg7ll+Hfy1wYN7j3+B5eC/wWydZ9yrgM8DZwOXAfwA75v0zTDQPFwPXDNdfCfz78PNux31irbnYVvvF8Ls9b7h+FnBg+F3/NbBnWP4h4FeG678KfGi4vgf46GaM0yN0oKo+BXxt1eLdwO3D9duBt65Y/he17F+A85NcvDkj3VhrzMNadgN3VdX/VNUXgMMs/zmILa+qjlTVI8P1bwKHWH4X9HbcJ9aai7W03C+G3+23hptnDV8FXAfcPSxfvU8c31fuBq5PcrI3Y07KoK/t1VV1BJZ3auBVw/Lt+CcP3jmcSrjt+GkGtsk8DA+Vr2b5iGxb7xOr5gK22X6RZEeSx4CjwP0sP/r4elW9OKyy8mf97jwM9z8P/OBGj9Ggr98Z/cmDRj4I/AjwWuAI8IfD8vbzkOQ84GPAu6rqG6da9STLus/Fttsvqup/q+q1LL8r/nXAj55steFyLvNg0Nf23PGHzcPl0WH5Gf3Jgy6q6rlhR/4O8KecePjceh6SnMVywO6oqnuGxdtynzjZXGzX/QKgqr4OPMTyOfTzkxx/P8/Kn/W78zDc/wOc+enMmRn0td0L3Dxcvxn4+Irlvzi8suFa4PnjD8M7WnUu+OeB46+AuRfYMzybfzlwJfDpzR7fRhjOdX4EOFRV719x17bbJ9aai+22XyRZSHL+cP3lwBtZfj7hQeCmYbXV+8TxfeUm4JM1PEO6oeb97PFL4Qu4k+WHjd9m+X/Wd7B8vusB4Onh8sI68Wz3B1g+f/ZZYGne49/gefjL4ed8nOWd9OIV6//OMA9PAW+Z9/gnnIefYvnh8ePAY8PXDdt0n1hrLrbVfgH8OPDo8PM+AfzusPwKlv/DOgz8DXD2sPyc4fbh4f4rNmOcvlNUkprwlIskNWHQJakJgy5JTRh0SWrCoEtSEwZdkpow6JLUhEGXpCb+D41qBHx4CqJTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(df_max)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "shape, loc, scale = genextreme.fit(df_max) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 212,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.3284055812905212"
      ]
     },
     "execution_count": 212,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This shape parameters tells us that the GEV is a Fréchet type, with a heavy tails, which means our river's flow can sometimes be very big, we better be careful and construct strong flood protections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the fitted distribution to see that it makes sense"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAD8CAYAAAB3u9PLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xt0nPV95/H3d266WlfLxviCDXZJRQg345Bb2w3ZYkgTsik5MT1tScseetLQJJtNt9C0SZc2m7BNQ5OWJEtLWpomAUK7rU9CQlIg3bYJvgEBjHEQtsG2bFm2Zcm6zP27f8zj0Wg8ksa6jaT5vM7R0fP85vc8+j2PZ+br3/Uxd0dERCRU6QKIiMj8oIAgIiKAAoKIiAQUEEREBFBAEBGRgAKCiIgACggiIhJQQBAREUABQUREApFKF+BcLF261NeuXVvpYoiILCi7du067u4dk+VbUAFh7dq17Ny5s9LFEBFZUMzs1XLyqclIREQABQQREQkoIIiICKCAICIiAQUEEREBFBBERCSggCAiIsACm4ewkL00NMTf9fTwzOAgw5kM58Vi/EJLC+9ftoy2aLTSxRMRUUCYbX2pFB/u6uLve3rOeu3h3l5+95VX+Njq1XzigguoCanCJiKVU9Y3kJltNrO9ZtZlZneUeL3GzB4KXt9mZmuD9HYze9LMBs3sLwvy15vZd8zsJTPbbWafnakLmk/2jYxw1a5dJYPBGUPZLH/86qtcvWsX+0dG5rB0IiJjTRoQzCwM3AtcD3QCN5tZZ1G2W4E+d18P3APcHaTHgT8EPl7i1J9z99cBVwBvMbPrp3YJ89Or8Tg/98wz7I/H82nvam/nwc5OHr/sMr6wfj2XNjTkX3t+aIhrnn6aHQMDlSiuiEhZNYRNQJe773P3JPAgcGNRnhuBB4LtR4Brzczcfcjd/51cYMhz92F3fzLYTgJPA6umcR3zykgmw3tfeIHDySQAtaEQ3+rsZOull/L+Zct4e2srH161imc3buQvN2ygxgyAY6kUv/jcc/xkcLCSxReRKlVOQFgJHCzYPxSklczj7mmgH2gvpwBm1gK8C3i8nPwLwR379vF08KUeMePbl17KTcuWnZUvZMaHVq7k8csvpy2S6845lU7ziz/5CQfUfCQic6ycgGAl0nwKec4+sVkE+CbwRXffN06e28xsp5nt7O3tnbSwlfZUfz9/cfhwfv+eiy7i2tbWCY95S3Mz/3LZZTSHw0CupvDe3bsZyWRmtawiIoXKCQiHgNUF+6uA7vHyBF/yzcDJMs59H/Cyu//5eBnc/T533+juGzs6Jl3Ou6Ky7vz2yy/nI+HmtjY+tLK4MlXaFUuW8O1LLyUaNB89MzjIb/30p7hPGldFRGZEOQFhB7DBzNaZWQzYAmwtyrMVuCXYvgl4wif5JjOzPyEXOD56bkWev/6ht5dngqaiulCIL2/YgFmpylNpb21p4Yvr1+f3v9bTwzeOHZvxcoqIlDJpQAj6BG4HHgP2AA+7+24zu8vM3h1kux9oN7Mu4GNAfmiqmR0APg98wMwOmVmnma0CPkFu1NLTZvasmf3XmbywuZbOZvnkgQP5/Q+vXMnaurpzPs9vnX8+v3Heefn9219+mcOJxEwUUURkQraQmiQ2btzo8/WJad/s6eFX9uwBoCkcZv8110x5BvLpdJrLdu7MD1m9rrWV777hDedU2xAROcPMdrn7xsnyaWrsDPliQUfyR1etmtZyFEsiER543evyPfWP9fXxwNGj0yyhiMjEFBBmwI6BAZ4KJpTFzPjtMjuSJ/K2lhb+26rRqRn/Y98++lKpaZ9XRGQ8CggzoHCY6fuXLWN5LDYj571r3TpW19QA0JtK8Yn9+2fkvCIipSggTFN/Os3DBSOBfmcGagdnNITDfKFg1NFXurvZdfr0jJ1fRKSQVjudpkd6e0kEHfOXNzZydVPTmNfX3vGdaZ3fgdqraoh3RHDgLT/YzvLt8ZIzAYsd+Ow7p/W3RaS6qIYwTYUrmf7a8uUzfn4D2vYkIZsLOom2MCMd4Rn/OyIiCgjT8Fo8zg9PnQJyN/LmEusVzYTosLPkYDq/f+riGK4RqCIywxQQpqGw7+Adra2sCDqAZ0NzVxJL52oJqcYQgyvV2iciM0sBYRr+6fjx/Pb7Z6l2cEY4Bc37RoedntoQJauWIxGZQQoIU9STTPKjYO5BiNzDb2bbkgMpwvEsANmaEKfX6FnMIjJzFBCmaOvx4/lVTd/W3EzHDM09mEgoC81do7WEgXWqJYjIzFFAmKLC5qL3LF06Z3+38XCa8HBQS4gZpy9QLUFEZoYCwhQMZzI83teX379xDgOCOTS/UlBLWKtagojMDAWEKfi3/v78ZLTO+nrWTWGZ6+lo7E4TKaglDKiWICIzQAFhCn5wcvRhcP95ksdjzobiWsLpdVGyGoUqItOkgDAF3y9oLvrFtraKlKGhO01kKKglRE0jjkRk2hQQztGRRILnh4YAiJrx8y0tFSmH+dh5CQMXRMnqX1NEpkFfIefoXwpqB29pbqYhXLke3YbuNOGRM/MSjMFVajcSkalTQDhHPygICJXoPyhkDk37x85L0BpHIjJVCgjnwN15MljMDiofEAAaD6UJJXIjnjJ1IYZWqJYgIlOjgHAOXo3HOZRIANAYDnNFY2OFS5Sbvdz0akEt4cJofga1iMi5UEA4B//W35/fflNTE5HQ/Lh9S15LjVkJdWSZZqqJyLmbH99oC0RhQHhbc3MFSzJWKJ0LCmf0X6Ragoicu6ppcJ7uoywBDr+1DhpzMfTeh/dw/8nd0z7nTGk6kOL0BVE8bCSbw8TbFOtF5NzoW6NMmSikg2BA1omdyla2QEXCSWg4NPpUtYF1mqgmIuemrIBgZpvNbK+ZdZnZHSVerzGzh4LXt5nZ2iC93cyeNLNBM/vLomOuMrPng2O+aGbzesBkonW0Xb6mP0tofsUDIFdLIFhjKd4R4YXBwQqXSEQWkkkDgpmFgXuB64FO4GYz6yzKdivQ5+7rgXuAu4P0OPCHwMdLnPrLwG3AhuBn81QuYK7ECwNCX6aCJRlfdMSp7xkt258dOlTB0ojIQlNODWET0OXu+9w9CTwI3FiU50bggWD7EeBaMzN3H3L3fycXGPLMbAXQ5O4/dncH/g54z3QuZLYlW0ZvVc08ay4qVDhR7es9PRwOhsmKiEymnICwEjhYsH8oSCuZx93TQD8w0TMlVwbnmeic84YbJJtGb1Wsf/4GhJr+LDUnc7WElDt/oVqCiJSpnIBQqm2/eFRjOXmmlN/MbjOznWa2s7e3d4JTzp5UYwgP54ocHskSSczvQZ1NB0ZrCV/p7uZ0Oj1BbhGRnHICwiFgdcH+KqB7vDxmFgGagZOM71BwnonOCYC73+fuG919Y0dHRxnFnXmJ5oLmonlcOzij7lgmvzR2fybDXx85UuESichCUE5A2AFsMLN1ZhYDtgBbi/JsBW4Jtm8Cngj6Bkpy9yPAaTO7Jhhd9OvAP59z6edIsnlhNBedYYztS/jzQ4dIZed/uUWksiYNCEGfwO3AY8Ae4GF3321md5nZu4Ns9wPtZtYFfAzID001swPA54EPmNmhghFKHwT+GugCXgG+OzOXNPPG1hDm5wijYg3daTqiubkIryUSfKtCzW0isnCUNVPZ3R8FHi1K+2TBdhx43zjHrh0nfSfw+nILWinZcK4PAQD3BVFDgNyid7+zciWfPHAAgM8dPMjNy5Yxz6d7iEgFaabyJJJLQhDKfYlGB53QwqggAPDB88+nLliA75nBwTFLd4uIFFNAmETh/IPYAmkuOmNpLMZvnHdefv9PDx6cILeIVDsFhEkkmscuWbHQfGz16vwY3++dPKnlLERkXAoIk1goE9LGc1FdHe9dujS//znVEkRkHAoIE8iGIV0f/P8660QHF15AAPj46tFpJN84dkzLWYhISQoIE0guCUEwKic65PNyhdNyXNPczFuDB/qk3PmClrMQkRIUECaQWlLQXHR6gUaDwO8W1BL+T3c3A1rOQkSKKCBMYEz/wcDCGmFU7Jfa27m4rg6AgUyGv9JyFiJSRAFhAsmCGkJ0YGHXEEJm/PeCWoKWsxCRYgoI43BbXE1GAL+2fDnLg+UsDiUSPHjsWIVLJCLziQLCOFL1NmbJ63BqkgMWgNpwmN9ZNbrI7J8ePMgEaxCKSJVRQBhHqmlx1Q7O+OD559MQLGfx/NAQ3+/rq3CJRGS+UEAYR2H/QWyB9x8UaotGuXXFivy+JqqJyBkKCONINo0uWRFdRDUEgI+uWpX/h/+Xvj6eOX26ouURkflBAWEcyUXWoVxoXV0d7yt4+pxqCSICCgglZaKQrcl1KFvaiQwvvo7X312zJr/90LFjvBqPV7A0IjIfKCCUUDjcNDqUZTE+UuaqJUv4Ty0tAGTIzUsQkeqmgFBCsrEgICzQBe3KUbjo3V91d9OXWgRja0VkyhQQSkgVBoRF1n9Q6Pq2Ni6prwdgKJvl/3R3V7hEIlJJCgglFAaE2ODi6z84w8zG1BK+cPgwCS1nIVK1FBCKOEU1hEXcZATwK8uXc34sBsDRZJKv9/RUuEQiUikKCEUyNUY2FowwSjnh+OKtIQDEQiE+UrCcxecOHiSr5SxEqpICQpFU4+iYotjg4hxhVOy2FStoDOcm4u0ZHuY7J05UuEQiUgkKCEWqqbnojJZolNsKlrP4zGuvadE7kSqkgFCkGgMCwMdWryYaPC70xwMD/OupUxUukYjMtbICgpltNrO9ZtZlZneUeL3GzB4KXt9mZmsLXrszSN9rZtcVpP83M9ttZi+Y2TfNrHYmLmi6xjwUp4oCwsqaGj5w3nn5/U+/9loFSyMilTBpQDCzMHAvcD3QCdxsZp1F2W4F+tx9PXAPcHdwbCewBbgE2Ax8yczCZrYS+DCw0d1fD4SDfBVVPMIodrq6mk1+b82aMYvebR8YqGh5RGRulVND2AR0ufs+d08CDwI3FuW5EXgg2H4EuNbMLEh/0N0T7r4f6ArOBxAB6swsAtQDFZ8Vlak1PJJrNgklnVCyugLCRXV1bFm2LL//GdUSRKpKpIw8K4HC5TAPAW8cL4+7p82sH2gP0p8qOnalu//YzD4HvAaMAN939+9P7RJmTnH/QTWMMCp255o1fCN4tOY/HT/OC4ODvL6x8ZzOsfaO78xG0cpy4LPvrNjfFlnoyqkhlPpeLP6v83h5SqabWSu52sM64Hygwcx+teQfN7vNzHaa2c7e3t4yijt11dp/UOj1jY28Z+nS/P5nVUsQqRrlBIRDwOqC/VWc3byTzxM0ATUDJyc49h3AfnfvdfcU8I/Am0v9cXe/z903uvvGjoI1/GdD8RyEavX7BUtjf/PYMV4ZGalgaURkrpQTEHYAG8xsnZnFyHX+bi3KsxW4Jdi+CXjCcwPZtwJbglFI64ANwHZyTUXXmFl90NdwLbBn+pczPamG0dsRqeKAcHVTE/+5tRWALPC/VUsQqQqTBgR3TwO3A4+R+9J+2N13m9ldZvbuINv9QLuZdQEfA+4Ijt0NPAy8CHwP+JC7Z9x9G7nO56eB54Ny3DejV3aOHEgXBIToUHV1KBcrrCX87dGjHE4kKlgaEZkL5XQq4+6PAo8WpX2yYDsOvG+cYz8NfLpE+qeAT51LYWdTNgrZ6OhT0sKJ6g4IP9/SwpubmvjRwABJd/7s4EE+v359pYslIrNIM5UDqTG1g+ocYVTIzPj9Cy7I73+lu5ueZLKCJRKR2aaAEBjTf1DlzUVn3NDWxhXBkNORbFZ9CSKLnAJCIN0wWieIDlVvh3IhM+OP1q7N73+pu5sj6ksQWbQUEALFTUaS8672dq4KagnxbJa7VUsQWbQUEAIpjTAqqbiW8JXubrpVSxBZlBQQADdI1482GUWGVUMo9M72dq5esgSAhLtmL4ssUgoIQLrOIJQLCOF4llCmwgWaZ4prCfd1d2tegsgipICAmovKcX1bG28sqCV85tVXK1wiEZlpCghAqmCEUUQdyiUV1xL+6sgRDsbjlSuQiMw4BQSKl6xQQBjPdW1tvKmpCYCkO3+sWoLIoqKAgJqMymVm3FVQS/jqkSPsHR6uXIFEZEYpIFA8S1k1hIm8o62Na1taAMgAf7B/f2ULJCIzpuoDQjYC2ZqgDyHjREZUQ5jMZy68ML/9SG8vO/TsZZFFoeoDwpjmomGv+kXtynF1UxM3FTys6M59+ypYGhGZKQoI6lCekj9Zt45wsP34qVP84OTJipZHRKZPAUFDTqfk4vp6fnPFivz+nfv2kXU1t4ksZFUfEPSUtKn71Nq11IZy92/X4CCP9PZWuEQiMh1VHxDUZDR1K2tq+MjKlfn9T+zfTzKreyiyUFV1QHAgVa/nIEzH761ZQ0sk9yTWrpERvtzdXeESichUVXVASNcZhHMBIZRwQukKF2gBao1G+YOCR23+0YEDZKIVLJCITFl1BwQ1F82I21eu5KLaWgBOpdP0r49VuEQiMhVVHRBSemzmjKgJhfjTiy7K759eHRlzb0VkYajugFBfuGSFRhhNx3uWLuXnm5tzOyGj72LVEkQWmqoOCGoymjlmxj3r1+dneo8sizDSXtVvL5EFp6o/sWOajPTYzGm7YskSfuO88/L7fa+rwdVyJLJglBUQzGyzme01sy4zu6PE6zVm9lDw+jYzW1vw2p1B+l4zu64gvcXMHjGzl8xsj5m9aSYuqFzZMGTqgsvPOpFhNRnNhD9Ztw5L5+5lakmI06sjFS6RiJRr0oBgZmHgXuB6oBO42cw6i7LdCvS5+3rgHuDu4NhOYAtwCbAZ+FJwPoAvAN9z99cBlwF7pn855UsX9h+MOKZ4MCNW1NTQvC+V3z+1IUZG3QkiC0I5NYRNQJe773P3JPAgcGNRnhuBB4LtR4BrzcyC9AfdPeHu+4EuYJOZNQE/B9wP4O5Jdz81/cspn0YYzZ6m/an8ulAeVQezyEJRTkBYCRws2D8UpJXM4+5poB9on+DYC4Fe4G/M7Bkz+2szayj1x83sNjPbaWY7e2dwrRwtWTF7zKHtxWR+f2hllHhrVXdXiSwI5XxKS3ULFjewjJdnvPQIcCXwZXe/AhgCzuqbAHD3+9x9o7tv7ChYg3+6xj4lTe1FM63uRIb6o6NTv092qoNZZL4rJyAcAlYX7K8CihesyecxswjQDJyc4NhDwCF33xakP0IuQMyZtJqMZl3rnuTYDuY16mAWmc/KCQg7gA1mts7MYuQ6ibcW5dkK3BJs3wQ84e4epG8JRiGtAzYA2939KHDQzC4OjrkWeHGa11I2R01GcyGScJpfGdvBnK5RNUFkvpr0v2zunjaz24HHgDDwVXffbWZ3ATvdfSu5zuGvmVkXuZrBluDY3Wb2MLkv+zTwIXfPBKf+HeDrQZDZB/zGDF/buDI1hkeCRe1STig5yQEyZU0HUgytjJBqDOERo+9nY3Q8m6h0sUSkhLLq8O7+KPBoUdonC7bjwPvGOfbTwKdLpD8LbDyXws6U4qek6f+ssyfXwZygZ1MdAMPnRRhelqb+WGaSI0VkrlXl0A8tWTG3ak9maTg02nR0sjNGVt0JIvNOVQaElB6bOedaX0oSjueCb6Y2pLkJIvNQlQaEsU1GMvvC6dyoozMGV0eJt1Xl209k3qrKT6SajCqjoSdDXcHchBOX1JCtynegyPxUdR9HDwWPzgRwLWo319r2JLFU7p6nG0Kc+hk1HYnMF1UXEFL1BpYLCOERJ6QKwpyKJJzWl0abjk5fEFHTkcg8UXWfRE1Iq7zGw2lqe4OmIzOOX1pDNjzxMSIy+6ouIKQ1wqjiDGh/IUkombv/mboQJ39WTUcilVZ1AUHLXs8PkYTT9uLojOWhVVGGO1RNEKmkKgwIhaucKiBUUsPRDPVHCkYdvb6GTLSCBRKpclUVEBw1Gc03bS8m8hPWsjXGiUtrzlpbXUTmRlUFhGwMstFck5GlnXBCXz2VFk7l+hPOGFkW4fQFWtdCpBKqKiCk6seOMNKidvND3fEMS/aPrnXUd3GMRFNVvTVF5oWq+tSN6T/QhLR5pfWnSWL9wQqoIeP4ZRqKKjLXqiog6Clp85c5LP1JIv+EtXRDiBOXqD9BZC5VVUBINWpS2nwWHXbad48ORR0+P8LgKvUniMyV6goIhSOMBvV/z/mo4UiGxoNjn52QaK6qt6lIxVTNJ82tYFE7IDKsGsJ81bonSXRgtD+h93LNTxCZC1VTH0/XG4TOLGqXJaQnOM5boSx0PJPg6JvryEaNTF2I45fXsmxnHJukYrf2ju/MTSGLHPjsOyvyd0VmUtXUEPSUtIUlOuK0PzfanxBvD3Nqg6oJIrOpSgOCmosWgvreDM0vj05aG7gwxtB5GosqMluqJyA0asjpQtT8Soq6Y6PrHR2/tEadzCKzpGo+WWMmpQ0qICwUBix9LjH6bxY2jl1ZQ7pW88xFZlpVBAR316J2C1goDcuejuefn5CtCXHsSs1kFplpVREQelMpLWq3wEWHnY5n4pDN/dulmsIcv0wzmUVmUlkBwcw2m9leM+sysztKvF5jZg8Fr28zs7UFr90ZpO81s+uKjgub2TNm9u3pXshEXhoezm9rUbuFq7YvS/vusSujnuyMKSiIzJBJA4KZhYF7geuBTuBmM+ssynYr0Ofu64F7gLuDYzuBLcAlwGbgS8H5zvgIsGe6FzGZsQFBXx8LWePhNE37RoPC4Joo/RdpOKrITCinhrAJ6HL3fe6eBB4EbizKcyPwQLD9CHCtmVmQ/qC7J9x9P9AVnA8zWwW8E/jr6V/GxPYWBAQ9JW3ha/lpivru0ZFH/RtinF5dNXMsRWZNOQFhJXCwYP9QkFYyj7ungX6gfZJj/xz4H8CE39BmdpuZ7TSznb29vWUU92xjaggaYbTgGbD0+QS1x0eDwsnOGMPL1cssMh3lBIRSTe7F7S7j5SmZbma/BBxz912T/XF3v8/dN7r7xo6OjslLW0JxH4IsfOa55S1ip4I1SMzofUMN8baqGCchMivK+fQcAlYX7K8CusfLY2YRoBk4OcGxbwHebWYHyDVBvd3M/n4K5Z9UPJPhQDye23HXg3EWkVAGlu2KjzYDho1jV9YSb1FQEJmKcj45O4ANZrbOzGLkOom3FuXZCtwSbN8EPOHuHqRvCUYhrQM2ANvd/U53X+Xua4PzPeHuvzoD13OWrpGRfJtUZMQJqYKwqIRTsGxnnHA89w/rEePYxlrNZhaZgkk/NUGfwO3AY+RGBD3s7rvN7C4ze3eQ7X6g3cy6gI8BdwTH7gYeBl4Evgd8yN3ndJ3RvSMj+W11KC9O0RFn+Y44ocRoUOjZWKvnMouco7KGZrj7o8CjRWmfLNiOA+8b59hPA5+e4Nw/BH5YTjmmQkNOq0N0yFm+PU7PpjqyNYZHjWNX17J8e5zYaf1HQKQci/6/UOpQrh6xIWf5jpHRJS6iRs8mNR+JlGvRf1Iua2jg2pYWwvGshpxWgdigs2xnnFCqIChcXavRRyJlWPSfko+vWcO/XH45q344Qm2fAkI1qBnIsnz76GJ4HjGOXVXLyFLNUxCZyKIPCFKdYqezLN82Mjr6KFg2e0iT10TGpfn+i1ilni88X8SGnOXb4vRcXUumPgQh4/jlNWT2JGl6LT35CUSqjGoIsqhFR5zztsVHH7BjRl9nDX0Xa5VUkWIKCLLoRRLOedtGRpe5AAbWRXPPU9AnQCRPHwepCuEULN8ep65ntKloeEWEno21ZLR6tgiggCBVJJTNLYi35NVUPi3RFubIm+tILtFHQUSfAqkqBrTuSdL6UiKflqkLcfSNtRqBJFVPAUGqjgFNB9J07Ipj6dG5CsevqKVvQ1SdzVK1FBCkatX3Zljx45Exix4OXBTj2FU16leQqqSAIFUtOuSs+PEItb2jnc3xjghH3lyn5ypI1dE7XqpeKA3LdiVo2pfMp2XqQvRsqqV/nZqQpHooIIgQdDb/NEXHrtE1kAgZpy4OmpBiFS2eyJxQQBApUN+bYcWPRqjpG53EFu+I0P3WeoaXaRSSLG4KCCJFIvHcw3YKm5CyMaP3ylqOXxojqxXAZJFSQBApwTzXhLRsx+iKqQBDK6N0v6WOET1fQRYhvatFJlB3IsuKfx+hoXt0FFKmLsSxTXUcf31Mw1NlUVHlV2QS4TQsfS5BXU+ak5fUkI0ZAEOroox0RGh7KYG7Y2YVLqnI9KiGIFKmhp4MK/5jhPqjo7WFbI1x/LJaNj/3HPtGRipYOpHpU0AQOQeRhNPxbIKOXXHCI6N9C9/v66Nz+3Y+sW8fg2k9fEcWJgUEkSmo781w/r+PsORACjw3byHhzv967TUu3r6dr/f04K4pbbKwqA9BZIpCGWh7KUlDd5rV1y1jx+nTAHQnk/zqnj3ce/gw96xfzxubmipc0tlRqUe0HvjsOyvyd6tBWTUEM9tsZnvNrMvM7ijxeo2ZPRS8vs3M1ha8dmeQvtfMrgvSVpvZk2a2x8x2m9lHZuqCROZazUCWp668kr993es4LzY6pfnHAwNc8/TTvPeFF3hxaKiCJRQpz6QBwczCwL3A9UAncLOZdRZluxXoc/f1wD3A3cGxncAW4BJgM/Cl4Hxp4L+7+88C1wAfKnFOkQUjZMYt553HTzdt4vdWryZWMOLo/x4/zqU7dvCbL73Ea/F4BUspMrFyagibgC533+fuSeBB4MaiPDcCDwTbjwDXWm4M3o3Ag+6ecPf9QBewyd2PuPvTAO5+GtgDrJz+5YhU1pJIhM9edBG7r76aLcuW5dOzwN8cPcqGbdv48Msvc0iBQeahcgLCSuBgwf4hzv7yzudx9zTQD7SXc2zQvHQFsK38YovMb+vr6/lmZydPX3UVm9va8ulJd/7i8GEu3LaN39q7V0NVZV4pJyCUmm1TPHxivDwTHmtmjcA/AB9194GSf9zsNjPbaWY7e3t7yyiuyPxxxZIlfPcNb+DJyy7jmoLO5ZQ79x05ws9s28av79nDHvUxyDxQTkA4BKwu2F8FdI+Xx8wiQDNwcqJjzSxKLhh83d3/cbw/7u73uftGd9/Y0dFRRnFF5p9faG3lR1dcwXcvvZS3FASGDPC1nh4u2bGDdz3/PI/39Wm4qlRMOQHo01fRAAALp0lEQVRhB7DBzNaZWYxcJ/HWojxbgVuC7ZuAJzz3rt4KbAlGIa0DNgDbg/6F+4E97v75mbgQkfnOzNjc3s6/XXEFP7z8ct7R2pp/zYFvnzjBO37yEy7buZOvHjlCPJMZ/2Qis2DSgBD0CdwOPEau8/dhd99tZneZ2buDbPcD7WbWBXwMuCM4djfwMPAi8D3gQ+6eAd4C/BrwdjN7Nvi5YYavTWReMjN+vqWFH1x2GT++4gre1d4+pm31+aEhbt27lzVPPcUf7NvHAfUzyByxhVQ93bhxo+/cuXNKx1ZqEo1Uh+lOlnp5eJgvHj7M3xw5wlA2O+Y1A65ra+O2FSv4pfZ2oqH5scCAJqYtHGa2y903TpZvfryzRKrchvp6/mLDBg6+6U386YUXsqamJv+aA987eZL37t7Nmqee4vf37aNreLhyhZVFSwFBZB5pjUb5+Jo1vPLGN/LIJZdwXWvrmOako8kkn3ntNTZs386bnn6aew8fpjeZHPd8IudCaxmJzEORUIhf7ujglzs62D8ywv1HjvDVo0c5UvDl/9TAAE8NDPDRri42t7Xxq8uX8672durDevazTI0Cgsg8t66ujj+58EI+tXYt3zlxgq8ePcp3T54kHfT/pd359okTfPvECepDIW5ob+eXly7lhvZ2miL6iEv59G4RWSCioRDv6ejgPR0dHE8m+VZvL3/f08OPBkbndA5nszzS28sjvb3EzPjFtjZ+eelS3r10KW1RPe9TJqaAILIALY3F+ODKlXxw5UpeGRnhGz09fOPYMV4q6GxOFtQcwnv38tbmZm5ob+eGtjYuaWjQIz/lLAoIIgvcRXV1/OHatfzh2rW8ODTEP/b28g/Hj/Ps4GA+Twb41/5+/rW/n9/bt4/VNTVc39bGDe3tXNvSQqOalgQFBJFFpbOhgc6GBv5g7VpeGRnh/wbB4amBsUuFHUwkuO/IEe47coSYGW9qauLtra38p5YW3tjURGyezHWQuaWAILJIXVRXx8fXrOHja9bQk0zyvZMnefTECR47eZL+gmUxku752sOngPpQiLc2N+cDxJWNjUQUIKqCZiqLLHDnOnM3nc3y44EBHg0CxHOTrLTaGA7zxiVLeHNzM29uauKapiZaolF9pubQdGdnlztTWTUEkSoTCYV4W0sLb2tp4TMXXsiRRIIfnjrFE6dO8WRfH68UPbxnMJPh8VOnePzUKSC3lMYlDQ2cuCRGTV+Wmv4MkSEvuda9LCwKCCJVbkVNDTcvX87Ny5cD8Go8zpN9fTwZBIlDicSY/A68MDQEq6MMBovbW9qJDWRzP/0ZagayChILkAKCiIxxQW0tH1ixgg+sWAHAwXicHw0M8KP+fv6jv59nBwcpXpjbI0aiLUyiLQzk5juMCRKns0RPZ4kOZQlpVe95SwFBRCa0uraW99fW8v7gGdFDmQw7Bgb4L9/aQaIlRKI5RLbm7E7nsUHiTKITGXaig0GQCH5Hhh1bON2Zi5YCgoick4ZwmF9obaV5XwrINSFlaoxkU4hkc4hk0/hBAjPSDUa6IcTI8oL0rBMdciJDWaLDWaJDTnQoS2QoSzg1J5clKCCIyDQZEEk4kd4M9b259qAxQaIpRGpJiGRjiHSDQakZ0iEjtcRILQlR/DigUHI0OESHnMhIrkYRGckSSpV+cLtMjQKCiMy4UkECwEOQagiRXBIi1WjB7xCZuvHnOWRjRiIWJtF69iqulgoCxIjng0Q0+B0ecULZEieUcSkgiMicsSzETuf6DQplI5CqD5EKmpNSDaPbHh6/DuBRIxUNk2oq/Xoo6YTjWSJxJxz3/O/CNAWNUQoIIlJxoTTUDGSpGQAKxjA5kKm10QBRHyJdb6TrQqTrDI9M3GCUjRnZ2PgBA0aDRjjhoz9JH7ufcCyz+JunFBBEZN4yIBJ3IvEMdSfGvuZANkYuOBQEiXR98LvWIDT5V3g+aExWlsxocAidCRpJzwWUVPA7CaFgeyHWPBQQRGRBMiCchHAyS00/UDQ7IhcwcoEhU1v4O5T7XWdkasoLGgAetlzgqS+zfGknlAqCRsoJJSkIHE4o7YRSQQAJ8oZSVHT4rQKCiCxKuYCR+/JloHSeM6OhzvqJGdmitIn6MkqeO2JkIkam7hzLnT47WNyyZw/rg2XOZ5MCgohUrfxoqMTE/y13wMOMCRiZWiMbNTJRC5qdjEw0VyvJxMqveZz1t84EktrRtL/r6eHyxkYFBBGRSjPAMhAadqLDk7fp5ANILBc0zgSJbDRIixnZSPBalDHbJedpAK1z8BAjBQQRkRmWDyAjDiPldwqcCSS54DAaLD5782UsnYNnYpf11Asz22xme82sy8zuKPF6jZk9FLy+zczWFrx2Z5C+18yuK/ecIiLVxoBQJjeyKnY6S+3JLPXHMvzmihW8e+nSWf/7kwYEMwsD9wLXA53AzWbWWZTtVqDP3dcD9wB3B8d2AluAS4DNwJfMLFzmOUVEZA6VU0PYBHS5+z53TwIPAjcW5bkReCDYfgS41swsSH/Q3RPuvh/oCs5XzjlFRGQOlRMQVgIHC/YPBWkl87h7GugH2ic4tpxziojIHCqnU7lUl3dxL8l4ecZLLxWISva8mNltwG3B7qCZ7R2nnPPZUuB4pQsxD+g+jJqxe2F3z8RZKkbviVHj3osZ+De+oJxM5QSEQ8Dqgv1VQPc4eQ6ZWQRoBk5Ocuxk5wTA3e8D7iujnPOWme0s5wHXi53uwyjdixzdh1Hz4V6U02S0A9hgZuvMLEauk3hrUZ6twC3B9k3AE+7uQfqWYBTSOmADsL3Mc4qIyByatIbg7mkzux14DAgDX3X33WZ2F7DT3bcC9wNfM7MucjWDLcGxu83sYeBFIA18yN0zAKXOOfOXJyIi5bLcf+RlNpnZbUHTV1XTfRile5Gj+zBqPtwLBQQREQHKnKksIiKLnwLCDDCzr5rZMTN7oSCtzcx+YGYvB79bg3Qzsy8GS3Y8Z2ZXVq7kM2uc+/BHZnbYzJ4Nfm4oeK3ksiYLnZmtNrMnzWyPme02s48E6VX1npjgPlTje6LWzLab2U+Ce/E/g/R1wXI/LwfL/8SC9HGXA5pV7q6faf4APwdcCbxQkPa/gTuC7TuAu4PtG4DvkpujcQ2wrdLln+X78EfAx0vk7QR+AtQA64BXgHClr2GG7sMK4Mpgewnw0+B6q+o9McF9qMb3hAGNwXYU2Bb8Wz8MbAnSvwJ8MNj+beArwfYW4KG5KKdqCDPA3f8fudFVhQqX83gAeE9B+t95zlNAi5mtmJuSzq5x7sN4xlvWZMFz9yPu/nSwfRrYQ24mflW9Jya4D+NZzO8Jd/fBYDca/DjwdnLL/cDZ74lSywHNKgWE2bPc3Y9A7oMBLAvSq3HZjtuDppCvnmkmoUruQ1DVv4Lc/wir9j1RdB+gCt8TwcKezwLHgB+QqwGd8txyPzD2esdbDmhWKSDMvXKWAllMvgxcBFwOHAH+LEhf9PfBzBqBfwA+6u7jPMQxl7VE2qK5FyXuQ1W+J9w94+6Xk1uZYRPws6WyBb8rci8UEGZPz5lqf/D7WJBezlIgi4a79wQfhCzwV4w2ASzq+2BmUXJfgl93938MkqvuPVHqPlTre+IMdz8F/JBcH0JLsNwPjL3e/L0oWg5oVikgzJ7C5TxuAf65IP3Xg5El1wD9Z5oRFqOitvD/ApwZgTTesiYLXtDWez+wx90/X/BSVb0nxrsPVfqe6DCzlmC7DngHuT6VJ8kt9wNnvydKLQc0uyrd+74YfoBvkqv6pshF9lvJtfc9Drwc/G7z0dEG95JrP3we2Fjp8s/yffhacJ3PkXuTryjI/4ngPuwFrq90+WfwPryVXPX+OeDZ4OeGantPTHAfqvE98QbgmeCaXwA+GaRfSC7odQHfAmqC9Npgvyt4/cK5KKdmKouICKAmIxERCSggiIgIoIAgIiIBBQQREQEUEEREJKCAICIigAKCiIgEFBBERASA/w95PH93tJvmwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(df_max, normed=True); \n",
    "\n",
    "x = np.linspace(df_max.min(), df_max.max(), 1000) \n",
    "y = genextreme.pdf(x, shape, loc, scale) \n",
    "plt.plot(x, y, 'c', linewidth=3);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**We can now calculate the probability that the river's flow exceeds 350 $m^3/s$:**\n",
    "Indeed, because the observations of $X$ that we have are independant, hence:\n",
    "$P(X_{n,n} \\leq x) = F^{n}(x)$\n",
    "and from the extreme value theorem this means $F(x) \\simeq H_{\\gamma}^{1/n}(\\frac{x-a_{n}}{b_{n}})$, taking the log and developping around $\\log(1+u)$ gives :\n",
    "$P(X\\geq x) \\simeq -\\frac{1}{n}\\log H_{\\gamma}(\\frac{x-a_{n}}{b_{n}})$\n",
    "Therefore we now have an approximation of the survival function in tail:\n",
    "\n",
    "$$ \\hat{F}(x) = P(X\\geq x) \\simeq \\begin{cases}\n",
    "                \\frac{1}{n}[1+\\gamma(\\frac{x-a_{n}}{b_{n}})]^{-1/\\gamma}        &\\text{for } \\gamma != 0\\\\\n",
    "                \\frac{1}{n}\\exp(-(\\frac{x-a_{n}}{b_{n}}))   &\\text{for }\n",
    "                                                        \\gamma = 0\n",
    "              \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 185,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def surv(x, n, shape, loc, scale):\n",
    "    if shape != 0:\n",
    "        surv = 1/n * (1 - shape*(x-loc)/scale) ** (1/shape)\n",
    "    else:\n",
    "        surv = 1/n * np.exp(-(x-loc)/scale)\n",
    "    return surv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 207,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.003669861860990634"
      ]
     },
     "execution_count": 207,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "surv(400, 4, shape, loc, scale)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 209,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "62.5"
      ]
     },
     "execution_count": 209,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "1/0.004/4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the probability that the trimestrial flow will be greater than $400$ $m³/s$ is  about 0.004, which means one trimester in 250, which means one time in 62.5 years! *Instead of \"never\" before*\n",
    "\n",
    "We can also answer questions like : what is the river's flow that we expect to surpass every century, because we also have the inverse of the survival function :\n",
    "\n",
    "$$ \\hat{F}^{-1}(x)\\simeq \\begin{cases}\n",
    "                a_{n} + \\frac{b_{n}}{\\gamma}[(np)^{-\\gamma} - 1]         &\\text{for } \\gamma != 0\\\\\n",
    "                a_{n} - b_{n}\\log(np)   &\\text{for }\n",
    "                                                        \\gamma = 0\n",
    "              \\end{cases}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 203,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def invsurv(p, n, shape, loc, scale):\n",
    "    if shape != 0:\n",
    "        invsurv = loc + scale/-shape * ((n*p)**shape - 1)\n",
    "    else:\n",
    "        invsurv = loc - scale * np.log(n*p)\n",
    "    return invsurv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 205,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "451.6079380637543"
      ]
     },
     "execution_count": 205,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "invsurv(1/400, 4, shape, loc, scale)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we computed the flow that we would get every 400 trimester, that is every century : $452$ $m³/s$!\n",
    "It's too bad these guys didn't follow my tutorial before:  <center> \n",
    "<img align=\"center\" src=\"http://www.planat.ch/uploads/pics/p0032_01.jpg\" style=\"width: 600px;\"/>\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Conclusion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**What's new :**\n",
    "- You now understand the problem with extreme quantiles in classical statistics.\n",
    "- You now know that there exist an extreme value theorem that looks like CLT that tells us there is 3 types of distributions for the max of any random variable, depending on a shape parameter $\\gamma$ : Gumbel, Fréchet, Weibull.\n",
    "- To fit such a distribution we need a dataset of maximum, and to estimate $\\gamma$, $a_{n}$ (loc) and $b_{n}$ (scale).\n",
    "- We can do that with scipy.stats.genextreme very quickly.\n",
    "- With it we can compute extreme quantiles and extreme value propabilities where we have no data to learn !!!\n",
    "\n",
    "**What's not new yet:**\n",
    "- It's sometimes hard to form a dataset of maximums, there is also a similar alternative with threshold levels instead of maximums, which is easier in general. And the distribution to fit is called a GPD (Generalized Pareto Distribution), that is in scipy.stats.pareto.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Thanks for your attention !** **God bless you!**"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}