{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Ensemble Kalman Filters"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"#format the book\n",
"%matplotlib inline\n",
"from __future__ import division, print_function\n",
"import matplotlib.pyplot as plt\n",
"import book_format\n",
"book_format.load_style()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"> I am not well versed with Ensemble filters. I have implemented one for this book, and made it work, but I have not used one in real life. Different sources use slightly different forms of these equations. If I implement the equations given in the sources the filter does not work. It is possible that I am doing something wrong. However, in various places on the web I have seen comments by people stating that they do the kinds of things I have done in my filter to make it work. In short, I do not understand this topic well, but choose to present my lack of knowledge rather than to gloss over the subject. I hope to master this topic in the future and to author a more definitive chapter. At the end of the chapter I document my current confusion and questions. In any case if I got confused by the sources perhaps you also will, so documenting my confusion can help you avoid the same.\n",
"\n",
"\n",
"The ensemble Kalman filter (EnKF) is very similar to the unscented Kalman filter (UKF) of the last chapter. If you recall, the UKF uses a set of deterministically chosen weighted sigma points passed through nonlinear state and measurement functions. After the sigma points are passed through the function, we find the mean and covariance of the points and use this as the filter's new mean and covariance. It is only an approximation of the true value, and thus suboptimal, but in practice the filter is highly accurate. It has the advantage of often producing more accurate estimates than the EKF does, and also does not require you to analytically derive the linearization of the state and measurement equations. \n",
"\n",
"The ensemble Kalman filter works in a similar way, except it uses a *Monte Carlo* method to choose a large numbers of sigma points. It came about from the geophysical sciences as an answer for the very large states and systems needed to model things such as the ocean and atmosphere. There is an interesting article on it's development in weather modeling in *SIAM News* [1]. The filter starts by randomly generating a large number of points distributed about the filter's initial state. This distribution is proportional to the filter's covariance $\\mathbf{P}$. In other words 68% of the points will be within one standard deviation of the mean, 95% percent within two standard deviations, and so on. Let's look at this in two dimensions. We will use `numpy.random.multivariate_normal()` function to randomly create points from a multivariate normal distribution drawn from the mean (5, 3) with the covariance\n",
"\n",
"$$\\begin{bmatrix}\n",
"32 & 15 \\\\ 15 & 40\n",
"\\end{bmatrix}$$\n",
"\n",
"I've drawn the covariance ellipse representing two standard deviations to illustrate how the points are distributed."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAEWCAYAAABRx5AbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl0Ved97//3s/c+0tGEEKABBAihATHZggMSFlg2cewm\nTeobX8eN46Rxkvbn9sZN47T3l9use1O7ies2ua5XfllOnI6pXTu1YydxZod4xCAjhLAwk0ASCARC\nMmLQfKSz9/P8/nh0NABisIUQ9ve1FstCZ599traWxec8+j7frzLGGIQQQgghhBCXlXOlL0AIIYQQ\nQoj3AwneQgghhBBCTAIJ3kIIIYQQQkwCCd5CCCGEEEJMAgneQgghhBBCTAIJ3kIIIYQQQkwCCd5C\nCCGEEEJMggkN3t/97ne59tprSU9PJz09nYqKCn7961+POeaBBx4gNzeX5ORk1q9fz549eybyEoQQ\nQgghhJiSJjR4z5s3j29961u8+eab1NbW8oEPfICPfexj7NixA4BvfvObPPLIIzz66KPU1NSQlZXF\nzTffTE9Pz0RehhBCCCGEEFOOutyTK2fOnMk//MM/8Cd/8ifMmTOHv/iLv+CrX/0qANFolKysLB5+\n+GHuueeey3kZQgghhBBCXFGXrcY7CAKefvppotEolZWVHDx4kPb2dm655ZbhY8LhMJWVlVRVVV2u\nyxBCCCGEEGJK8Cb6hDt37uS6665jYGCApKQkfvSjH7Fo0aLhcJ2dnT3m+KysLFpbWyf6MoQQQggh\nhJhSJjx4l5SU8NZbb9HZ2cmzzz7LnXfeySuvvHLe5yilxvy9s7Nzoi9LCCGEEEKISZOenn7W5ya8\n1CQUCrFw4UJWrFjBQw89xJo1a/jud7/L7NmzAWhvbx9zfHt7Ozk5ORN9GUIIIYQQQkwpl72PdxAE\naK3Jz88nJyeHDRs2DD8WjUbZtGkTFRUVl/syhBBCCCGEuKImtNTkr//6r/noRz/K3Llz6e7u5oc/\n/CGvvfYaL7zwAgD33XcfDz30ECUlJRQVFfHggw+SlpbGXXfdNe45z7VML96btm3bBsCqVauu8JUI\nMbXI/xtCnJv8vyGmmguVS09o8G5vb+fTn/40bW1tpKenc+211/LCCy9w8803A/CVr3yF/v5+7r33\nXk6dOsWaNWvYsGEDKSkpE3kZQgghhBBCTDkTGrx/8IMfXPCY+++/n/vvv38iX1YIIYQQQogp77LX\neAshhBBCCCEkeAshhBBCCDEpJHgLIYQQQggxCSR4CyGEEEIIMQkkeAshhBBCCDEJJHgLIYQQQggx\nCSR4CyGEEEIIMQkkeAshhBBCCDEJJHgLIYQQQggxCSR4CyGEEEIIMQkkeAshhBBCCDEJJHgLIYQQ\nQggxCSR4CyGEEEIIMQkkeAshhBBCCDEJJHgLIYQQQggxCSR4CyGEEEIIMQkkeAshhBBCCDEJJHgL\nIYQQQggxCSR4CyHEe4Dva3xfX+nLEEIIcR4SvIUQ4irn+5rNm2HzZq5o+JbwL4QQ5yfBWwghxLs2\nVcL/aPJGQAgx1XhX+gKEEEK8O57nsHatHv74cooH2cv9Ou9W/I0AwNq1espfrxDi/UGCtxBCvAdM\nRrA8X5idzPAvhBBXKwneQgghLpnvG+Ds8D1VyBsBIcRUJMFbCCHERYmHWd83VFcrYGqXcUzV6xJC\nvH9J8BZCCHHRbJi9uA2LV0s9uBBCTBYJ3kIIIS7JxZRxyOZGIYQ4mwRvIYQQl+x8gXsqkVV3IcRU\nIsFbCCHEhBi7ym3/wJULvbLqLoSYaibsp9Df//3fs3r1atLT08nKyuLWW29l9+7dZx33wAMPkJub\nS3JyMuvXr2fPnj0TdQlCCPG+NFUHxXieI2FXCCFGmbCfiK+99hp//ud/zhtvvMHLL7+M53l88IMf\n5NSpU8PHfPOb3+SRRx7h0UcfpaamhqysLG6++WZ6enom6jKEEOJ9ZSImRk5UcLe133aleyID9zu9\nvst1PUII8U5NWKnJCy+8MObv//mf/0l6ejpVVVV85CMfwRjDt7/9bb761a9y2223AfD444+TlZXF\nD3/4Q+65556JuhQhhBAXaaLLMSY64J55fZf6OhK4hRBTyWX7idTV1YXWmoyMDAAOHjxIe3s7t9xy\ny/Ax4XCYyspKqqqqLtdlCCHEe9p7dVX3XKvcvm/e9eq+EEJcSZdtc+WXvvQlVqxYwXXXXQdAW1sb\nANnZ2WOOy8rKorW19XJdhhBCvOe9m8A9FSc8jqxyG8rLYe1aNfSIOs+zhBBi6rsswfsv//Ivqaqq\nYtOmTSh14R+U5ztm27ZtE3lp4iog33Mhzu1K/L9hjA3jSk3eCrMxDocPp3PoUBIHDvSzbFnn8Osn\nJ9vrqauTFW8xQv7dEFNFUVHReR+f8OWNL3/5yzzzzDO8/PLLLFiwYPjzOTk5ALS3t485vr29ffgx\nIYQQU4cxDrt3T2P37mnDAXwyKKVZtKib3Nx+lDJnPTaZbwKEEGIiTeiK95e+9CWeffZZXnnlFYqL\ni8c8lp+fT05ODhs2bCASiQAQjUbZtGkTDz/88LjnXLVq1UReopjC4isW8j0XYqwr9f+G72v6+uzH\nK1ZMfinKypVTqwRGTD3y74aYajo7O8/7+IQF73vvvZcnn3yS559/nvT09OGa7rS0NFJSUlBKcd99\n9/HQQw9RUlJCUVERDz74IGlpadx1110TdRlCCCEmyJWu/74SrymTLoUQl9OEBe/HHnsMpRQ33XTT\nmM8/8MAD/M3f/A0AX/nKV+jv7+fee+/l1KlTrFmzhg0bNpCSkjJRlyGEEGICvZ8CqEy6FEJcbhMW\nvLW+uJq7+++/n/vvv3+iXlYIIcQ7NNGru+c637t9DVmBFkK8l1y2doJCCCEu7EoFy4le3T3X+S72\nNca7B5O9An2lS2uEEO99EryFEOIKmexgeTlC/rkG2fi+AS7ut6BTrbzjSr++EOK9TYK3EEK8D5wr\n4L7b1d2x57R/fN9QXa3GfO6dvsb5rlFKUIQQVyMJ3kIIcYW8k/B7ocB5KYF0okJrEGh8XxEOu5y5\n0n2h14jfA983+D7A2FXvcz1/qq2SCyHExZLgLYQQV9ClhMYLBc7zPX456pc9zyES8amthepq9a5W\n0qurob7eUFRkqKy8/CvZsmIuhLgSJHgLIcT7xKWG/As9x/c1W7dCQwOUlBhAjXlONBoM/X3s58cT\nBIbzNcgafU3v5k2ErJgLIa4UCd5CCDFFXCjsXihwXmwgvZhylXN1KDnXc1zXoahIU14+9rFoNODx\nxw1gV7Fd1x035HqeQ3l5gA3u6qKvSQghrjYSvIUQYgq42GB5MTXTZ5539OffSYAd7zkjQd85z/Ua\n6uvBdc1QOD/3+a3JCdPSNlAIcaVI8BZCiPeod7pKPHrD43gr0KOPPZdw2OXuuwN832Hr1vixZ59r\n9DVGIppwePwQX15u8Dw1IWFZArcQ4kqQ4C2EEFPAZK3Cnut1Rq+Kj+7LHW8LWF4eX62+tNBru5xA\nRUW81nu85xrq6w1BAJWVZz86Es7VUMtCWa0WQlydJHgLIcQUMdFBcrwwHw/Y8QAbX3EuLw/GhG1Q\nBIGmuhrioRds8O3tDTAYcDWOUqQlJ6DUuVe0R/p6n73qHq/vDgJbL34hY/uES623EOLqIsFbCCHe\nwy7UBzsesM9+nuK66wzb97Xzzz88TntXBwm/O8H+Ix0cfbuXnugggQmGjw95DjPSkpg5LYmFc2ZQ\nVjKH1SW5rCjMAZLHvPaZ1xUOu1RWXtym0QuVvgghxFQmwVsIIa5i77bswvPUqOmSLrkFJ3lp+0Ee\ne+ggL20/SEdn37jPVdiV7kBreqMx2k/10n6qlz2HOvjlG/uHj1u7bB6f+uA1LOlezJ63kuznLtAt\n5ezrHHlcNkYKIa5WEryFEOIqMroG+52UXcSfP3qUe/Ox0zzzyi6efmU3dY1tY47PnDaNxblzKF8+\nkwQ/k7xZs/jg9anU73ZI8FxWrQpRV+ewd98gJ7v7yZzTT9LME2xvOEbNvla27Wtl864WNu9q4b7v\nvkBlyWI+t/4GYMa72vwphBBXIwneQghxlRgJqoYgMAC4Llxs+cXooLtmjeY3NY38089r2bCtCW3s\n+aYlJ1Kal09kYT6fvy2fjpYZKKUoLx8J+ZnTNS8f0ezYAfv2aUpKDMqEWLPCo6IinXA4l0/fcg2+\nr+nuG+AXb+znP3/3Fi9tP8CLO3fx2t49VB9ZxVfvup7RZSijrxMuLmDLRkshxNVEgrcQQlyVDK7r\nUF5uP46H73N1KBkdSmN+wO/e2smf/lsV+1pO2Mcdh2tzi/nYdcv50qcLqXszBEBJHpDH0Dnc4RIP\n31fMmwfNzXZDZFmZoqLCdjw5u194InfdtJzP/N61HGo7zf3/8SpPbNjBd36ylSc2vMX3vvQR7rhx\nyTvqMy5DdYQQVxsJ3kIIcZWIdwDxfcPWrWpo1VuN6kQSnPPjtWs1g37A93++jUeefYOjHd0A5GWn\nc+9/K6MwbQnHDqdQUgKJCc6YMpQzXz/epSQx0eV//k9NYqIz3DZwPLYfuCYvZzr/8dcf4y/+ezlf\n+f7veOnNg9z1dz/m9Z2H+f+++HuEvPOfRwghrnYSvIUQ4ioxujVffMXbrnaPzxjDjzfu5Sv//CKH\n2zsBWLogk/+2ooIPLFuC0YqmJli61BCL2RAfH1Rz7muIv54iJcW94Hj73t6AqiqGR8YDdB/L5v/c\nehclmVv5pxdf5rGf19Bw5AQ/f+hOkhJDF9w8OXolXzZaCiGuJhK8hRBiirhQvbLv29pu13WoqLDj\n10eXgJz5cfqcNr78vd/yal0zANcWZPPgH3+AWyIFbNpkwHbiJgjsH8cZ3bcbIhEfz1PDK9rRaEB1\nNQSBGXr98VeofV8TjWqefhqOHoUbbwyAkeOVcviHL5axeO4cvvbMc7y4/QD//Ws/4md/dycJofOf\nV8pLhBBXKwneQggxBVwoUBrjUF2tcF3be3t0eceZw3FOdfdz/w9e5Xs/qyHQhhlpSXz9c+v5s1sj\nuK4tF4lvyoxENKGQAlzKyw2+b3jqKUhJ8dm1y+B5LnffbdsOVlUZGhqgpETheeqcbxTsYB67ETMW\nM3ge5OZCWdnI1MuRNwcef/rJ+VSUfZobv/wEL9Q0ctc3fszT938c7xzDdOLnHqlntyUsl9KWUAgh\nriQJ3kIIcVVReOP85DbG8IPf1PG//vlFOjr7cBzF/7h1FR8uuZFpyUkYM/Y8YKittSG2osKGed/X\nFBUFHDyo2btXMWOGDeOep3Bdh4KCgEjEPn/0xEuwQbi2Vg2tykMo5PKJT+jhVfNzhWPPc1icl8mf\nX/9JHnnlKX78+l6+9m8v8/f3fHDM1zZ6bLwd+sO49ewSvoUQU5UEbyGEmAIuVK+slB530yNAT/8g\n9zz8C/7r5V0A3HBtHt/54odZkpc5HJDPfC3fh+rq+LkMvm9Da1mZpqXFYd48w4oVEA7bbiWRiE9t\nLdTW2pXy+MpzVZWmpcXQ22sAh6VLbbcV+wbhwt1KPE9xU/lsZmR9nP/3yaf41tNVfLi8iMpr88a5\nVzK9UghxdZLgLYQQV8B4q7/nYoxzzsfjpRdNx07wia8/x+7m46SEQzz25Y/w6ZuvQSlbDhLfLHnm\na3keQwHcjFk9BohGHaZPZ7hNoO9rtm41vPoq5OZqggDAEIlAXZ2ipcUwMKBYtw7Kyxle4Y6H7UhE\nE4sZHEfh+4p4iUj8WiorIRLJo6bpOp554w0+982fUfevf0pacuLwMWe+MTl3bbusdgshpi4J3kII\nMckutVf17t3TAFixYuRY39ds3Kh55uV6nqz5BX2Dg5TMn8VPvv6HLM7LPON1RsbCxz9/Nlu+UVVl\ncF1FQYENyeHwyLW5rkNOTsD8+fbve/eC1orVqw3RqCIUMpSVqbPaCwaBprYWDh2C/HxNVZWD66rh\nsGyvyfD004r1CyrZ09bEzoNv840nNvKtP7t5+JhztTcUQoiriQRvIYS4TC51w98lTWwMNP/04kv8\n6I0tANy2djHf+fPfZ9b0pFHnM8RiGsdRgDPcaaS21j4eBKC1wXEAFGVlDNVoa9asUcMlJvFrKi8P\niMXiZzd0dCjCYfv8cFgxMABbthjWrQsIh91RJS2K6mrw/TO/3pGV9khkqADdeDx230dZ96V/59Hn\nt/JXf3gd2TNSL3jfpNOJEOJqIMFbCCEug/PXNJ9dNjHe8Z7nsHRp15hj20728Im/fY6Nbx3Ccx0e\n/NxNzOyP8K2HFCUlAZ/5jBnuQnLoEBQW2pBbVWVobNQcOwY5OYr8fENDgyEUUhQWAihiMU1TE0Mr\n0me/GWhqsgF+8WJYsACiUWhstMe1tSkSEgwQsG4dw+Hb82wJSyTCGavhI7Xa4bDD3XfbjZzhcC5/\nUFHML6r2862nN/OPX/g96VoihHhPkOAthBBXwPgB0uD78Y2JllIjJRmbdh7mD//2WY6d6GH2jFR+\n+LXbWVMyj6eeCobOa9iyxfbodhyD7zuMDriuq8jMNCxYAImJcOqUYs4cOzzHdjgxQ329NdGoYssW\nW3ISidjzFxVBS4s91x/8AdTV2SE+AwPguhAOG/bvdwiFGDViPr6ybcP86K997BuQkQFB/+fTlfyi\naj/f+9k2/uqOChr2pAwff66SE6nxFkJcDSR4CyHEBIuvzp6vC8mZ4qUcVVWG6mrnnAHzid/u4I//\n78/xA831y+dz3y23ozpT8Tz4+McVt96qSUhwqKmBpiZFUZHhzjsNKSnO0AZGTUWFGZ4+WVvrMGdO\nbLhmOxrVrF5tCIXs61ZVaTZtguxs266vv9/hi1/UNDQ4vPqqQmvDmjX22quroahIE4mooTBuz7d1\nK9j68Xjv7ZF2iOdbxV5ZNJsPlxfym+pGnn1tD6WZq4feDEA4LPXeQoirk/ykEkKICRQvGYmXjVxK\nIIz3yrbnMUNdSzTGOPy69gif/ebz+IHmvtvX8KuHPsX05OShMGonRD77rKKuThEKKYqKNA0Nipqa\nkXPZsg9Fba1Dba1iyRKflhbFli2KIAioqTH84hcQi2nAlqGkpGgWLQqGrzEhwaGgALKzDYcOqeE+\n4GCvPTHRobRUEwQBW7faczQ0KCIRu2mzuloNf13x+zR2BLztgOL7hk9+YBkAz7yya6h9oR3us3nz\nyL0590ZRIYSYmiY0eG/cuJFbb72VuXPn4jgOjz/++FnHPPDAA+Tm5pKcnMz69evZs2fPRF6CEEJc\ntUYHz+pqGzA3btT86y9O8LfP7MAYePDz6/nY8pup2ergurZ0JAgMR44YWlvtJEeA0lK7eXL/fti0\nyefVV4NRkx819fWGN9+05R/xkfGDgzZEaw39/Zr0dI3nGfbu9Vi2zPDFL2oSEx3WrXP49KcV+fkA\nCt83LF/uE4tpqqo0Tz1lePVVW7ZSUqIoKbHXObrkZTzxyZmPP66ZoQtJ9Dze2HOEX73UidZjj7PB\n3RCNBuOfUAghppAJLTXp7e3lmmuu4e677+Yzn/kMSo39IfvNb36TRx55hMcff5zi4mK+/vWvc/PN\nN7Nv3z5SU8+/a10IIa4GF6o3vtAmQd+3NdpNTVBSAi/v2s2/bHoDY+BvP3sj/+uT60YNxBkJswsW\n2M/EYpr9+xXRqGLpUhugq6sN/f2alSsNO3d6xGIBRUXguh5f+EKMHTsgFPK47z6N6yqeegpOntQk\nJEB/v8LzDNXVcOQIJCRoCgshFrMlKaWlAc88Y+joMKSlQXGxDdnZ2YZIRJGYaDd51tY6Z/UTP9cG\n03gZDEBSQoiKRcW8snsPr+/dx8P3lbNu3egBOob6eluTXlERnNWrXAghppoJ/Qn14Q9/mAcffJDb\nb78dxxl7amMM3/72t/nqV7/KbbfdxtKlS3n88cfp7u7mhz/84URehhBCXFG2pOPcofvM8orR5RLR\naEBVVcD+/VBQYDgc3cPf/+xnGGP4WOk13LDgesDWjldUKMrLobzcjmlPTDR4Hhw4YFfCGxsNAwN2\ndPvp0xAKaerqNLGYz/79Cq0VpaUBe/a4HDrksGuXZvt2e52FhYZw2J5v1Sr7OikpitZWAENra8CO\nHQFNTQFBYDunHD2qKCkxrFvncscdhmXL7MZL3zfD5TOjg3E0GgyNoh/5++bNhq1b4Y47NEVFAA53\n3rIQgLa+I4TD7qhOKXY6pn0DYcP96PsqhBBT0aRtrjx48CDt7e3ccsstw58Lh8NUVlZSVVXFPffc\nM1mXIoQQU8LYyY4+mzcbNm+2w2uO+fv43P99Hq0NHytdzq3XXDPmuSOTJm3tdFeXwhhYvFhRWmr4\n6U8VBw9CQYFdOW9rg/p6h+LiAGMUQeBSWwsNDYaSEs2+fbBli2LvXsP8+ZreXojFFNddp0hO9giF\n7Er08uWa55+HPXsUs2fbUF1ZydCqs20V+PTTiuPHNQUFPkHgUlamh0K3fbynx+eJJ+x13323LRPZ\ntMm+4UhMtKv4Nqwr1iyda7/e+qNn3T/PU1RU2N8S2E2c57/X9jmyIi6EuHImLXi3tbUBkJ2dPebz\nWVlZtNpllHPatm3bZb0uMfXI91y8lyUn2+BXV2c3TR48OA3XVRw6pGhuTiYlBd48spdHn3idQBs+\nunwZf7D8GpYu7UapurOel5TUTVIShMN2umVqahdNTZqEhHQSEmDatG6WLXNQKhXXddmwIY3+fpes\nrNN0dwe0tyeTkTFIRoahqyuRAwcSOHbM0Nfn0d8Pv/71KYqLTxMO238umpoGSUzMIC0tnYSEQfbt\nayM93V7Tzp0arT1Onszi1KkEjh4dxHH6ePZZg1JmuB95Q8N0Dh5MJz09xo4dHezenczOnelkZcVY\nuLCfw4cDFi/uQylNf0dAUijEkeNd/OrFzWRPtyPkjXHYvdveg3j7xUWLuqmr88+65/FjAZYu7RrT\nnlG8N8i/G2KqKLK/rhvXlGgneGYtuBBCvFeNDn1KaZYu7RoKhmkkJmqO+Q08+rIN3Z9dX8j6Bdei\nlEIpPfxcpTRLlvSwa9c0du2axtKlPcMr0vHjiot7AXAcG0SXLz+N1mH270+luxsGBgxLlw6iVJhX\nXkknJyfKhz50kubmMAcPhpkzZ5DBQUVrayLJyRm8/noaShk+8YmjFBeforAwChiUsqF73740jFEo\nZViyZADfjxGL+RQU9LF7dzKuq9DaQylNEBiys2MEATQ1JQMwa9YgeXn9LFjQA8DevfbzixdHmZ2W\nxYGTR9lzuHs4eJ8pCIwEaiHElDdpwTsnJweA9vZ25s6dO/z59vb24cfOZdWqVZf92sTUEF+xkO+5\neD+IbySM1z1fc03AY/91kEef3EigDX91RwXf+Nx6tm5VHD7cDEBp6UrAlktEowFvvqk5ftwwb54h\nN9fBdV2WLzcMDGieecYG8bvuckhN9fB9zcaNmrw8zU03abSeza9/DT09htRUg+clEo3O4CMfsV1S\nXFcxOKh59lnFqVMaYwxB4FBSshjXVXzrW3ZIz403gus6dHfbYxYv1mzfrmhutvXXs2fbmvP+fjuO\nfulSh1tv1WzZAk1NkJMzk4qKkfKP6mo7sj4/315Hb29AfnYTB04epeV4OqWlK/G8+PTMeAmLGr4v\n41mxQkpN3ovk3w0x1XR2dp738UkL3vn5+eTk5LBhwwYikQgA0WiUTZs28fDDD0/WZQghxBUzeuPf\nxo2ahgZbf712rabtdDd/9/xPCLTm42vKKHRuZONGQyik8X1bLrFxo31+ZaWtby4p0bS2Gl5/3eGm\nm2xLv6oq2LtXc+gQDA7C66/DzTfb57muQ2EhlJU51NaC42iSkqCy0nDwoEN9ve2KcuAALFyoCYUc\nlDLMneuQn2/7aIOivz8gISH+lShisYCBAcPJk5CebjueDAwoTpwArTV2H78aer6mtlYRChmWLgVw\nht98xO+P69oOKAMDhn/8Rwe/Nw2Akz29w/dx7P27cKCWwC2EmAomvJ1gQ0MDYH/YHjp0iLq6OmbO\nnMm8efO47777eOihhygpKaGoqIgHH3yQtLQ07rrrrom8DCGEmHJGb6QsLzdjHuuLxrjta89woquf\nD60u5H/c8kHq90JTkyExUZGQYIP30I9XKirsqvb+/QalYM4cTWOjw+7dmr4+28WksNBw/Lji8GHb\nZzscdolEbNlJaqpHeXlAaandxOh5Hp5niMUMELB3r8OJE4Zp0wI8z2XNGnBdly1bNN/+tu3tnZlp\nmD/frrD/9rcBJ0/aYH3kiENpqWHOnICEBIeBAYecHMXq1Ybt2yEWg1DI9g4f3Q4xLt5y0Pdt+cjc\nuTC3Pw0Og0rqHhPQR+6tDfQSroUQU92EBu+amho+8IEPALZu+/777+f+++/ns5/9LP/+7//OV77y\nFfr7+7n33ns5deoUa9asYcOGDaSkpEzkZQghxJTmeYrKSkVZmcZ1DX/yj7+irrGNwtwZ/Nff3E5q\n2KVstaa21iEW06Sl9eI4PiUlNqQODGiefDLgrbdg5kzD/Plw4ICmuRlmzdLk5ioWLIBFi0aC7ehO\nInfeGaO21l5LWZm9ntLSgO3bNU1NDmlptowEFPPnm6ER8DZYnzgBiYmQmalwHMWvf21X3BcuDFDK\nITUVkpOhrc0lGoUgUOTmQigETU0OBQWGSMQQDo8NydFoMNSpRRGJ2FKZQ4cU111nuH3RNJ7aBm0n\ne4bun+2kUlFh38BUVwMY1q69tPAtnU6EEJNtQoP3jTfeOPRrxfHFw7gQQrwfjB2HPvKx79tg/aOq\nLfzwpZ2kJiXwzNfuIOyF8DyH1FSH8vKA6mqH/ftTWLSom7VrFdGoprbWloNkZBgWLjSEQh7Lltmh\nNz09DklJAb/8pUtxsaG4WLNpkyIUsi0Fs7M1dXXwyit2tbylBebMsePd7Uq57RHe3q7wfUUopGls\nhMJCh9WrHQ4c8Hn7bVi4ULFsWcDx41BSYstSQJGYCK++6uB5hpwc222ksBCCwKGoyAb/2lo7nTNe\nex7/TUAQ6OGe30Fg0NrgOA4z0sOA/c1AnO3lbQN7fb0N4OXl9vUu9vsSf91LCewS1oUQ78aU6Goi\nhBDvRfERb/VMAAAgAElEQVRwFwSaigpb7jHymGHHwcN8/3cvAvAvf3Urb26ayZubNJ/6FEODYhRB\nENDcnIzvw5Ilms2bDVu2QGGhJjkZBgYctIbly+3I9wMHDH19toZ6cNDQ0gIDAy4lJZrKSsPKlQ47\ndzo4ToDjGI4dAwjYv98lHDbccotPZmaY2loFBOzZozh61JCf7xMEdgOn52lqaqC62hCNKlJTFd3d\niqwsSEgwxGKKVasUN9zgjpoyCaCoqgp45RXYvVvx+c8HYx4vK4NwGMAZ+vrM0N/H53mKoiIz/PF4\n3wf7+LsLy+80rAshRJwEbyGEuIyCwG4CdN2RsOb7mg0v9/P1536GNoYvf3wNHy1fxLdeM7iurd+2\ndc52hdhxwBi7UtzcbOukHcchNxdee83h2DEboLu6FIWFdtU5FjOcOuWQkGBXpBsb4dgxh1DIIRKx\nw3F8H1pbDQkJmmnTFLEY/OAHHtddpykvh5oaxbFjMHu2z+bNisZGw8KFkJdnV8/37TPEYvba5s9n\n6JodFi40pKa6w+UkGzcGaG1YuRJc14xpIWsnUAbDY+XXrrVvSmIxxf79LgMDkF9qxrm7I2Un8Y/P\nZMOyff7osHzmbyCEEGIySPAWQoiLdKkrp57nUFFhx7aP3kRojOHhX/yKE72dLJozmwc//wESQg5r\n1wY0NkJtLbiuHm7rV1gYZcGCPqJRO5a9oEDhOC6RiA3pR44Y2tsNnZ0OxcUOZWUeg4MBv/mNXTFe\nvFhz+LBDEIxcg+MoXNeQlGTo6VHMmaM5edIhIcGwfbvm7bcdcnIgJ0eTkRGwZ0+IaNQwd25AU1OI\n3NyA3l7o7FRUVgasWRNi504H0EOv4xCNavr7A157DU6eNOzYYRgcdLjpJk1pqTv8GwDPU7iuGrrH\nhupqWxpz442GMzdfjnefx/+emXFLUS41cEtYF0K8WxK8hRDiIrzTMgPPU0OdOkbC2r/95k1er68n\nLSmBHz1wG8lJIQBuuMEG5VgsoL5eA4qFCzX79iUxMGCYO9d+bt06B9831NZCR4eth45GISsrhut6\neF6IsjLF7t2aUMjwk5/YH/Vf+EJAUpJLba2hpUXT2KhITzccOuQybZrhC1/Q7NkTYts2xcAArFyp\naWsz1NS4pKZq0tLgxRcVvm83hYbDkJioaWx0UQpCoYADB+zGzuXLY/z4x4q2NoPvw6xZoBT09Rl2\n7bJhOhLxSU31xgTaeNB2XZdIRLN1K7z1lh0r7ziXPmxtvFKUd1p+IoFbCPFuSPAWQojLZCSsK9au\ntZ/b33KC+x59AYA7rv19Nv1uOsXzAsJhuwJsyy6grU2Rk6M5elTR3R3CGDW08dB2F9myxbB5swYC\nenrsJsjWVpeODgiFbKjs7QXXVWRmBrS3K377W8PJkwGHDyuyszVz5zocP+4wc2bA6dMue/c6JCY6\nlJUZkpNhzx6HxkZDb6/ixht9WltdXNewenWMwUGXjAxDd7fiyBFDX59m9my7WbO9HVpaNF1dmo4O\nl4ULYckS20Kwv1/R1qZ49VXD7t3w+c8HQ/XsI4F2JIQ7QEBVTT8A01MuUPB9DucqRZFabSHElSI/\nbYQQ4jzshEk9tCp7ccNaxoq35rMlJl969AUGYgGf/MByckNLaGuz7QGj0WDodWyZRl6eGmofqMjO\njlFS0k8kYohEbCh1HEVODoTDCs+D1FRIT7edRACamzUNDTbkam1fv6HBbrZMTTWEQjasz55tWLAA\nMjJg+3aF42hCIU1Li+2NnZWlWLZMsX69R3a27Scei9nzNjW5nDzpEA7bkpV58wzXXWeYP1+ze7d9\nvLjYDLU1dGho8EhKcrj1VsOsWQyXl5xp9P0tK1MkptrgPSMt6RLu+9jzSbgWQkwFsuIthBDjmIiV\n0SCIbwxU/GpLAy9sbWRaciIP/9nNNOyxq9c1NXZYjp3CqKisdIaG5Nga5ePHDS+/PINt2zThsKKw\nECIRQ2mp4plnHDo7DWlpLrm5sGCBS2mpLenIzLSBuKPDJQggJQWmT4dPftLwy1+6tLZCerrm1CmX\nWMwMT6N84w1Fe7uhqwumTbN9uHftcunq0syeramvt/chK0uTnW1IT7elKfX1CmMUqamalBQDOCxY\n4KC17bxiSz4Ue/Z4LFumKS01hMNn/zM09r475C/qgRqYMyt1QjqUSK22EOJKkeAthHhfuFL9l+N9\nqQcGfb786G8B+KPrbyArI4WsoQ4er702Mo49fo3RqM9zz9n66IKCQQ4ejM9oN7S2BuzbBwsW2K4l\nixYZcnMVhw45BIGmrs7l8GGH6dM1M2bY8O37EAo5zJ9vJ1cWFdke4NEoJCTYoWd5eS5LlgS0tto2\nf4mJDHdVicWguNhh2TLDo486hEKGnBzD8eMuixfDwYMOx4/DjBkMtR3U5OcHLFnisG+fOzwGHqC6\nWtPYaDuXrFs30lLw3F1JDC3Hu+w97JrGxo221/e53ghdyvdYArcQ4kqQ4C2EeM975xsjL25ldLzA\nN/r5//jsVhpbT5I3axa3la0aftz3A44cAaXsREfPc4db4O3aZcs3PE8xc+YAn/iEw+Cg5tFH4fhx\nQ0JCjMFBO6a9p0dx4oR9rcJCzfHj4PuKxESHo0cVGRmKggLD0aOKQ4c0AwO2deHp07ZtYVmZQ2mp\n5l//VXHypCIpSXHkiO1q0txsN3MuXGhD+403Gvr7A958UzE4CKWltvY8P98QiSgGBxXf+Y5DXR30\n9EB5uSYSAVAMDGh27dJD/cMhFDJD3VtGb7BkuK1gdTVsqTsFQGZa2oR/j4UQYjJJ8BZCiPMYHeBs\nvbcNwmM36p3dJzr+GEBrRzffeGIjAN/98i3cUOYO9/P2fYPv21ru0V03XNfWbc+cqTl+PIFo1NZ/\nDw6C1oaEBDtS/e23XWbO9AFNV5eH4zh86EMBzc02zM6ZY3tux9votbUZWltBa8jLM/i+pq/P4fBh\nTUKCoqPDlpgoZUhKsl+n60I0Cvv3KxzHsHKlYft2GBiwdeWW7QNeU2OnTa5Zo6mutiPjtYbaWtvW\nLxTStLcrZs82FBbar1dr23Ix3krQ3kvbhcQYTcuptwG44/dzWDB7pA+3EEJcbSR4CyHe8yaiptf3\nNRs3BtTXQ0mJorIyvmJ97j7RIyuwmgd//Dt6ozEWzSwmNZY/5vEggDvvNCQmOkMDaGx4X7fO4DiG\nvXs9YjFYuDBKX1/Ab3+rmTlTAw6+r0hIgK4ulyDQpKUFrFql2bxZ8dpr8T7coLUtF3EcRUaG7Vji\nupCUBG+95ZKYqOnqgqYmh+uvN2zfbuu277zT0Nhov6AlSwJ27NAcPGhobbU9wH0fTp6EuroA34fm\nZjNUbqL4zGdc1q0zw73In3sOjh415OQobrpJUVpqe3dv3WpLYGzLxdHt/uzf80t66B3sJyMtTF7O\ntLO+h6N/2zDe91jGvAshpgoJ3kKI94WLDV3jhTTfNxw8aFePi4tHJimef2S54ZevH+bFnbvxHJeb\nF9405tGRqZYO5eWMKZUIh10qKkBrzeHDfRQU9PHUU1BXB4sXQ1aWorhY0dOj2brVcOKEIi3Nrog3\nNTlkZhqysmwXEtdVhMO2u0piomLhQkMoZGhpcUhOhvR0RVeXXdkuKFCkptpQ3thoA73Wmp/8xNDa\najh50sHzIBKxmy7b2gxvvKGYNk2jtSEUcigu1oBLTY1dJS8pgfx8TV6eYfVqdzhw2+Ptan+888tI\niYm9l90hu9q9PD+bqio1fH/ivzGI37NIxMfz1PBQntHfTylBEUJMFRK8hRBiyPlCmucpCgsdCgsN\n5eUjpSbjjSz3PIfS0kH+n+//DoAvf3wN//uuDBITR547dqrlSJj3fUNvbwywq8HhMICD48CcOTAw\n4JCR4XDttYaf/1yhlA3ic+YYjh/3WLFC09xsR8a7rn0dgL17FadPG5TS1NQ4ZGTA+vV29by/33D4\nMLS0KHp74dQpgzE+mZlQV6fo6FAkJNjx8tOnK3JyHH7v9+C55+DYMejuNjgOlJUZtHbZssVw9Kjm\n5EnbznDvXsPx44pQyBCLaZqaHIIASkoMkUhAVZXCdeO/NRh5A/PWgXYAli/MGnWPznyDo3niCXsv\n7747OCt8CyHEVCHBWwghLkK8lMH37d/jJSHxx87k+5oH/+0tGtrayJ2Vxv/+9DrqtttAOFISoYY6\nfRjCYZdIxCcI7IbCn/40YNo0uPFGTXNzMgcPprBuHdx+uy0PcRzFnj3Q2WlISFCUlkJamu0eUloK\nhw6B79vV9FiModeAgQFbKhIOg+cZtmxxSU+HadNsUF+0SNPZadi3Dw4etG8IcnIC+vocolHFvHma\n1asVWjvs3AkJCYbkZJ+2Nodp0+xGy+3b7Wp2NOqQmQmLFwf8/OeKvj5De3vAzp0ehYUBubkKcKmr\ng4MHDQMDhiBQVFaq4ZXvv//VIQBWFc8Zas1o8H1neBLo2rWaaNS2Mjzf922875MQQkwmCd5CCDHk\nYkJadTXU1xuKisxwnTecXaIS8wP+41W7ofKz199EYig0fI54KUUQ6KEhMopIxOeJJwzGaBYv9uno\ncOnuhqKigBdeSKStLcTcuT67drl4nuKOOww//ani9GnF7Nma7m7F9OmG0tKRc86cqVm0KGDXLttH\n+0MfCvj+9xXt7Yo/+INBtE5g715wXbvhMjPTloKEQj5Hj9qyktZWxfz5huxsW84yfbqL646sKOfl\naTo6FP398VaChqNH7X/nz2foDYLDggWGmTM1M2Y45OdruroccnMVZWVQW6soLNREo2q4/SLAxk0B\nr75pg/fNqxbSuEdRXw9gqKwc6YCSmupx9912rPy5VrslcAshpgoJ3kIIMcp4IS3egWS8x84sUflV\ndQMd3d3MmzmDD16zDM8bWcWNi3fzAEUQGDo7NSdPQlGRw/r1BnBJSXEIhfqZMSPA9x06OhSZmZCQ\nYANsNGpbBIbDAdu3u9TVBSxdamhvVxw+DCdPGsJhn/Z2h2PHbLlILAYvvRRi8WJITjZ4nqGnB5SC\n/v6AN95QxGKG7OwAY+DAAYeUFKishPR0+8/GwIDGcRShkENOjqGwMCAjQ7F9u8uBA5q0NENBgQ3g\nr77qMGuWZsECl/p6B6UCEhJGVqiDQFNfb2vB7W8AFNGo5s0DzQz4PtcWZJObOY2Z5QFBYFfxx3ZA\n0VJeIoS4KkjwFkKICxgJ1rY0xNZ4q/OupH7nx9UAfGz1KlatGhmOE1/pjkQ0dXUOQWAoK9MEwciK\ncXOzy4IFUF5u+3DPmDFIXt4geXkzyM+HSMSWqcRitqtJUpIhCGyLwIYG6OuD8vIAx3Ho6gLHsa0C\nfd9l1aqA11+Hzk4F+HieS3e33ZQZjWpefFHh+w7Ll9vr+d3vHIyBjAzDqVMeN9xgqK011NZCZqZP\nURHk5ir27vU4cQJmzw7IyID+fmhpsV1V7ERNh3gtu++7LF0KsZhhyxZwnPj9sfdo40ZNQ4Ph1Z1N\nANwcWQjY1eyRlW57Hy1Z0RZCXB0keAshxCU4V+A+s0TlzYZjvL7zMIleAjdfs5ytWxnaOGjb69n2\ng/E/sHWr7ZgSBIoVK8Dz7NAaz4Ply32UUuzfn8icOZq0tBB1dRAEAU1N0NamCIUCTpxwyMsLiEbh\n9Gm7WXL6dE1np8Pp0y7XX2+orFSAR3JyQGtrjAMHPHxfMXeurdV+6y27GbOyMqC93WX2bJgxw9Ze\nJyUZenvjAdle+969isZGw8qVtsf4wAA0NBjefhvy8qCvTwEOn/rUSI/y8nL7X983PPWUPVdBgWHp\nUmf4DU38vuw8ehCAW1YVjLnX9vnxkhoz1Hpwwr/VQggx4eRHlRBCXIR4n+nzT7A0gOaxn20D4MMr\nrmH99WFqa+ObMBVlZTZkNzQ43HmnGW6t53mGwsKRVWDQJCYaqqoUp04lEot5bNumWbdukKamEAMD\nhsWLYelSaG5WQ5sNFT09duhOba3L6dN2wM2JEw47dxocR3PokO3BvXt3CGM0M2fakAx2UmVGBmRk\nKDo6HIwxRCKa3bvh+HGH2bM1tbV206VSAVrbY958U5GXp1m82K7ouy6kpQV0do6Uf9iyEDudMxy2\nPctLSgJaWjQHDigKCw2e5wx3iVmwqIs/e/o4yeEQldfmjXvP6+vtbwkqK6VVoBBi6pPgLYQQ5zG6\nzGTt2vMdY1eys3P7eOrFnQB8489Wk5rqDXVDsavF4bBtSQg2xIJhcFDT1weDg3DgABw8CLNmQVKS\n5u23PTwPMjM1jY0h8vM1sVhAS4vdDDl9uktvr8uMGZrmZkUQKJKTDd3dirQ0zcKFAQMDcOKEornZ\n0NNjR8UnJsLAgA3hbW22hWBuru18snevQ2qqz4kTir4+RWamRmtFe7tDZqYmFApQCpYutbXcr7wS\nIhaDxYsNK1YE1NQ4NDd7rFmjAcXmzQ7hsKax0bB7N5SUaCoq7KTLri5FYqItk4nzPIeNO5sBuPHa\nBbjOuaeHlpczXPMthBBXAwneQoj3hcmZXmh48sW36BuIcdPKfErmzyIatd02Ro9Cj0QMmzdDUxNM\nm6bZt8/Q2alIS4P58x1AEw7DrFkeH/2o4tixLk6fDjNzZiJBoCgp8TlyxKOpyWHFCjsaXilDOKwo\nKDDk5LgsXKjp6AjYujWE50FBgWHRIti/H8AOzpk1y6C1Yu9eG/p933ZGyczUtLXZlevTp+0QnrVr\nAzo7FQ0NcPSoS1eXYv58zdq1kJenmDfPsGeP4tAhl+xsQ0KCPd+2bXal+0//VBOLKfbti4dsRSjk\nUlwc4DgMd0qJf5+efnkXAIUzFrJ5s6G8PMDzFJs366GuMorKSmeofEY6lwghrg4SvIUQ73nvZnrh\nSP/ukQ2S4x2zYgV843O1APzpR1cNbRKEkhKGRqfbbhy1tYrmZhuGw2GHri7D3LmaefMUkYgDOGzZ\nYl+vogJ27+6mqSmgtHQ6u3bBxo0e8+YZlFIkJdnWgtu3O0MbMu3XVltrqK938Ty44QYoK4PGRofj\nxw0zZhh27HC59tqA7GyYPt1ucPR9zdGjdjrnjBkuS5fGeO65ECdPwsyZimuusd1P2tshMREyM20r\nvy98AUDx9NO23jojQ5GXZ0fVFxRo9u932L7dIRQypKbaaZXxLi8wthvJ5s1wuKOD39Y0keh55IWX\n0NioCQJbptPYaNsVFhSY4XsvhBBXCwneQghxEUa3rjtX2PM8h1d2NHD05Cmy0qdxy8pCampseI7F\nNCtXQnLyyCbDkpKR7iRLlvg884zhd78DWwdtV4ODwJajvP76DObNG2DRIsPzzyu6uuzKeSikaGiw\n5SIHDkBamqaqyg7IaWuzZSIrV2qOH1d8//suy5cHHD3qUlPj8KtfhWhrU5SVaW6/PSAWs1Mljx9X\n7N3rUlRkOHRIMW1awPTp0N3tAoq5cxVZWZqUlIDNmz1279bEYg6OoyguDlDKYe5ch1DI1l/Pn++w\ndq3t5+37hqYmhe/bmvlzjXcPAs1PqmsA+OA1y7lmSdLQiriD5xny8yE/31BR4UjoFkJcdSR4CyHe\n8yZjeqHva/7hP+2mylsjEWprHUIhxR13BPz0p7Bxo+LGGzUVFQ7l5bYPdm2tLcNoadHs2GHbAu7b\nFxAE3tCYd0V1NSQnw7FjYfbscVi+XJOUZFizxsN1FS0tmpYWjedpduyw/cBnzAjo6PCYNg2ysmDz\nZpfkZENHh2LZMp/58w2dnfCFLwzQ1+fw4x+7zJxpmDULUlIMxcV2mE1dnUdysmb9+gDHCREKKfbt\nsyPmFy+2LQi7uw1tbZpTpxx6ehzmzYM1a6CuzgEMoRCsWWPfbNTWKubO1cydO/IGJH7v4voGovy2\nztbIf+3zESKLRsK55420ExxvcNGZ35PxHhNCiCtBfhoJId4X4h0z3ulz1661q8zjneNUdz9bG5vw\nHIf//fnSoQ1/ioQEB9+Ph0xDVZXhtdc0P/6x5qWX7CZLOwrdDrOZP18Nre7acBqJGNLTY/i+DbJz\n5ig6Olxqamz5im3hBwMDkJKiOXnSQSlFWZnPrFmao0cV118Pf/InthZ7/36H9nYbgNvb7Up1VpZG\nKcO+fS6nT7vccINPfr6iuBgKCxWHDyfgeQE9PZpTpyA11bB8ORQXG9raHHp7FSUlhv5+RW+v7T2+\ndq3iU5+CUMihtlaxZYuhvh4SExWjhngOlwFt3my7wry4czf9sQGWz5/H0rzss753oz+ORgM2btRD\nzx21O/Os8459TAghrhRZ8RZCiIswOnCfayX1tbcOoY3h+mvmsyB3GnOz48fYceYDA7bDx5Ythv37\n7WbKuXMNpaW2t/aRI5qMDId161w8z+D78WE7ho6ORFpawrS0aPr6HNrbDcnJEArZke0DAwHNzQ4Z\nGXbD5JEjLn19mtOnFampUFio2bHDJTVVk5kJxijmzw/Yts1h2TLNjBlw5IjD/Pm2BeArrySwcKFm\n2jSXefMMr75qaGhwmDfPoDWcOAE7digWLHDp7TUsWKAoLdXs328IAsXAgMbzXMLh+P0xOI7d+Ok4\naszI+dFcFzbstr81+OiKCE89BSUl5pzlPb6vqa62bzqKiuwqvBBCTHUSvIUQ4hKMt1HzN9WNAHyo\nrBA4e2Xc9r8GpTSJiYpIBCoqFCkpNoSuWxdjyRIIh52h85vhyYzFxYPEYiGCIJHk5IBrrtGcOuXx\n+uswd67Pjh0evg8LF2oGBqC723Yj6euzGyBbW6G/X+P7UFAQkJ2tOHXK0NsbDG90dBxDTo5BqYA3\n37T13MuWBbS0OEN9uTUzZihycjSNjbBvn8vddysiEcXWrYaf/1wxY4amr8/w3HO2ZWJ5ue3gYstK\n7J94icno+xPvkf5q3SH2Hu4gPSmVuYmLhn4TcO4NrVZ8zPzZ93syyouEEOJSSfAWQohLdOao8lgs\n4Jeb7Xjzm1bmjzk2vjLb0mJQyhCNOpSU2NKSujpbwrJ8uc8jj8BLLzn81V8F2E4ftovH1q2wb1+I\nIPCZOdOugkejDuGw4eRJOy3ScQxz5sDv/77iyScdMjI0OTm2JWF+vsO+fbbdYEaGz6ZNHrm5NoQf\nPeowe7YhJcUG323bPDIzDWvXBnR1ubz9tiIpCYqLNceP25KS5mb7J/41hMMOWvscOmRHwweBnaaZ\nnx9QXe1QXx9fkbb17GVlitRUb/jejO6R/r2f2U2VNxatYOkSl0hEEQ6fe2jRSLAef6iRBG4hxFQj\nP5WEEOISua4aGldu1R/u4O2uLqYnJ7N0fvY5a4qjUcXs2Yo777RdSxoaHBobNb29AXV1drOjMXaS\nZbyePBx2cN2xNdEZGXY8ekJCABg6Ox0WLdJUVMC0aSFWrzYUFdkpkJ2dDkeP2tHzJ08qWlocwmHI\nyNAcOuQQjSpcV3PwoO3ZnZ4OmZkOt93mEok4JCU59PZCU5Pi1CmYOxdmznSYMcOGbLB12aGQQ3a2\nLTOx9d1QVhavyTZDx2leecXw1FOGnh5/zL0JAs2B1tP8rGofIc/hrvUrCAI1NOFy/H+m3k3dvhBC\nXAlX5CfW9773PfLz80lKSmLVqlVs2rTpSlyGEEK8Q/HSCeulNw8C8KHyhdTUOGM29HmeQyRiWLAA\nWlsVzz6rqKmxq8ADA7BlC7S0KCoqDF/+su2LPdKxwxCNGmIxh+xsO2inrw8SEw1tbS4ZGZqCAk1C\ngkdrq8uTT2pqaqCxERoaFMeOKfr7DSUlhoICzeHDDomJkJWlmTUrICND09XlkpamKSyEe+4xzJ8P\nu3a5rFwJhYUBsZgN3p6nue46xR/9kcv69Q7Tpzts3Wqorrar3IsWKZYuVRQW2o2VALGYIS8PysoU\niYmKxETDkSOGrVvN8P2JRDSuC3/7L7Vobbj9+sVkZ6TLNEohxHvSpJeaPPPMM9x333089thjrFu3\nju9+97t8+MMfZs+ePcybN2+yL0cIIca4UAu6c9UO/7bGlpn83uqCs85lB+Y4hEIBhYVQX69wHCgt\nNQwM2Bpm34do1MV1GZ50WVVlO540Ndl67unTBzl+XNHTY/ttG6NYuhQ++EGPpCSHujpobg5ISPA5\ncsSlu1uxYIGht1eze7dDQYEmN9dl1ixNerpLUZFm3z5FKKRIS7MdWJ5/HvbvN8ya5dPWZmvDQyFN\ndjYcO+awc6ftPR4KKZqaIBw2zJunKStTbN0KTU2G+fMDamtdGhttx5WSknhdt8P119sw7rouvm+o\nrlbEYoY393Tyo6qtANz7sTLWLGHM/RVCiPeKSQ/ejzzyCJ/73Of44z/+YwC+853v8MILL/DYY4/x\n0EMPTfblCCHEsIudcDmmw0mgeW1HMwAfKi9k1rSR4zZvjk+stCPRy8sNkYghCAxPPaU4edJw/fW2\nPltrzXPPOQRBwKJFdiNmT09AW5sdvLNiRSeOk0FXlyEU0iQn///t3Xt8VeWd9/3PtdZOsnMgHAIh\nHJNAEgOBCgYJBMFDVfBUnbYeCgPt9KA96dRO2/txbmfKTL3bTuv0eUZtp3q3t7d1arVnrbWiU0sV\nEgKCUc4JIUA4JJxDEgjJXms9f1zZO9kQQtAQTt/36+XLsLP22tfecZvvvvit389j69YEXnzRMH68\nDbS1tXDwYIgRIzwyMiA9PeDgQTh+PKC+3pCdbcN2dbXhsssMl19uGDkSbE15QHV1QCjkc/gwbNxo\nOHbM1lhfe63P7t2GTZsMra328VNToaXF6Rj/boflbNpkh+RcdpkN3QUFdtBNOOx2fFjp2s3E7oo7\njuGVzW8Q8SNcWzSRGRNHK3CLyEWrX4N3W1sba9as4Rvf+Ebc7TfeeCNlZWX9uRQRkT5RVXeAY8cj\nZA8fyKCUZEIhh9ZWr2PEvK3Rtp096Nj9Jta3u63NAD61tXbQjOME7NgREAQBmZke77zjUF9vKCiw\nF1vm5bkdxxs2b3Y5cMCWnXh2kxxjDEEQcOSIQ0YGtLcb9u93GDUqQlqaS1ZWwNtvG44etfXco0c7\nHQsEbzgAACAASURBVB8I6JgqGSEhwaehwQZkG97t1MhRowKqqnxqamwv7uPHAwoKAhISQoRCAQUF\n9rHz8gxTpoAN1k7coByI/9Aya5bPm+/WsXTDRpKTQvzkH68H7Aeg07VvFBG5EPVr8N6/fz+e5zF8\n+PC42zMzM6mvr+/2Pm+//XZ/LE3OI/qZy7mUkmLDXWVl74auvFa5G4BkbxCPPtrADTcc5PXXhwBw\nww0HWb8+paMsA3bsSGbv3kSGDWvnuusOAj5btw6gvDyD5GSPz3xmL4cOpbNpUwKRSBtJScmMG2eH\n4/zud0OZOnUnOTlNpKUNICEhnZEjoaDgKIMGtbBpUxIFBQ5jxiTyzjsp7N/vM2pUO8ePp7Bvn8fc\nuXuAdvLzh1Bfn0RWVjOFhY04Thtr14Lvh9i9ezh794YZNaqNgQPbSUryaWuDl19OZPhwj9RUH7BD\nc3zfZ8CAfRjjs3atz8CBMG1aIuvWpVFXFxAOQ1ubT0pKIwDr1g0EYNKkRoyxr217BO577C0AFl09\njv07t7P8zwPwvICioiMY4xMEDuvX279GiN4mciL93pDzRX5+fo/fVztBEZEuosEuCJy4P59K9e4m\nALLShsTdbjuROOzeHcYYQ05OK8YEDBsWIRSCmpoweXlH8H2fyy47ztChbbz6ajrHjyeQn3+c1NQ2\n7rrrEODw/PPDOHDAjog3xicnp4kgMOzbF2bp0jQikXSGDImQk9NEenobU6dCY6NDU1MCAwfaspDf\n/jaDAQPaGD/+OEePJrBzZwpJSQH5+c1s2JBEQoJLcrIhHPbIy2uhqiqJ/fuTGDfuOKNGedTWJuE4\nMHbscZKTfbKyjgGwebMNyhMnNrN5cwrV1QPwPBgyJMLYsUdjr+WuXckAFBU1xW77z9/vZsueI4wY\nnMyCq/NYv34Au3Ylk519rC9+lCIi551+Dd5Dhw7FdV0aGhribm9oaGDEiBHd3mfatGn9sTQ5D0R3\nLPQzl3Mtvtb71CUOkYjPpqeqAPi7u3K469rhhMMjKSqKsGJFQEtLJrNnAxhKSuD4cZ/ycp+KCsOx\nYxlMnGhobrbTJq+4wuexxwzNzQHp6ckcPerQ2uoweXLA9OkNHWF2FM3No3Fdw9ChHkeO2A8FrguZ\nmYnccksqv/yloa7O58orfTZudBk7NuDyyz1eeimZhISA0tIIrgvvvedQUZHKli3DGDUqYOtWh/Z2\nmDYtYMiQFNraHMaNC5g0KZl16wIiEYfc3ICmpiSCwOFTnxrCmjW2bKWw0DBxIjQ3QxD4RCI+27YZ\nmpuHMHnyWEIh210FYOLETFavNjQePcrv310CwGMP3MrM6YX47XYSZ2mpIRzurAefOlWlJtI9/d6Q\n801jY2OP3+/X4J2YmEhxcTGvvfYaH/vYx2K3v/7669x55539uRQRkT6xbd9eAPymYbF65nDY6djx\nthdThsP29rIyw9atdmhNbpc5O64LiYkO114bUFcHnuewaxc0NHiUlcGOHelkZ7d2XAhpz+s4toVf\nKOThugHgsnq14dgxOyp+3z6HnByf/fsNvu8yb55PXV3AU085NDfb8e2HDkFdHSQmBiQmBqSlwciR\nLp5nyMnxGTMGfN9h3z7D9Ok+c+fCT37SGX5ragyeZ1sOrl6dwPTpPsXFhtWrHY4fh/HjIRQysU4w\nkUjAsmUBNTUBf9ryV5paW7luai4fu3oCnhd0TLA8uTe3AreIXCz6vdTkq1/9KgsXLmT69OmUlpby\n4x//mPr6ej7/+c/391JERE4SvZBv1iz751DIOeXFfY5j2N9kSyea9g0iErEXUYZCDiUlHmVltpVg\n9Fyu61BYaNvvhUKGigqACNu3G2pqDB/9aOeI9MxM2LEDDh4MOHbM5dChEJMn+yQnu6xY4VFdHdDa\naqdEHjnikJUV0NYWkJLiM2yYw969hpEjI6xbF6KuznDrrXaSZWOjbWeYnW0YNizCxo3Q3OwSChky\nMmDyZLAdThyqq+loWRgwebJDYqIN0I5jSEmxEzi3bDFUVRmKinxWrrT39X0oKnLjRrnb19Gjpsaw\n/UA9f1izBtcxPHb/PDwviJtgKSJyser34H3XXXdx4MABHnnkEfbs2cPkyZN55ZVX1MNbRM657kpM\nemoxuO9wCxHPJz05mQ9NSozr4BEKdU63tIHcUFJig3U47NLcHGH9ep8gsBdeHj7s8ctf2vsVFhpK\nS2H2bIe2Np+//OUQW7cm8NJLhrw8j/Jyn4YGGDTIMHSoTyRiaGx0SEz0aWoyZGX5ZGd7vPuuS12d\nITUVDh600zNzc+Hyyz0SEhJJTTUcORLQ2urR1OTi+/DCCwG+b5gwITrm3bBli2HduoBQKKCw0KGk\nxD6HkhKP9nb7HIuLYeXKgKVLDaNGwYIFAeFw56+Y6IeXyy6D//y/r+P7AX//sRKKcjO7nfQpInIx\nOicXV37hC1/gC1/4wrl4aBGRDywaFPcebgFgdGYqs2aZuFDetbyiosK2DbQ9vW1gXbkyYN++gKws\nwy23+Dz9dMA77xiysz08z7B+vUthYYDrOuzZEyIx0bB9e0AkYqdXpqTA4MEe9fUO+fkRDh2CQ4cM\nKSk+dXUOu3YZmpvtTvjo0baXuOPA0KGGqio7rGf8eMjP9/jv/3YYNMgnMzPgrbdcBg+GSZNsUPd9\nyM0N2LjRPi/bCtHEXoOEBAPYDxbTp0N1tf1+1w8h9sOLrfHecOgd3t2+g6EDU1j8qWtix9gyE6Oy\nEhG5qKmriYhIhxOnUp5YdgJ2sqQN0tCeaoP30IGp3QZGe5ttwbdpk+23XVjoEYnYsBoO2+mPiYkO\nQ4YYPM8nJyegsRE8z8f3DeDx7rsppKTAiBEenueQl2dwHJ9t2xx27jS0t4coKYmQnu7S2upw7Jjt\n4X355T633GJwXYcf/hBaWmDyZI89e0IMH+6xdWvAgQMOBw4Y2toMGRkBOTl2GuWWLSFaWwMuu8ww\nfbqD60J7u4/rOh1j4u3z2rLFp67O0N7uc801Lp/8ZNDx3E2sH3ckYp//7sZ9fP+N1wD4jy/PIy2c\n2OX1VJmJiFz8FLxFRLqIBugTy04gfhIlGJqOHgcgPTUp7hxda8Kj9d6eZ/B9e9/Vq20wHjQooKXF\nUFnpsGCBvaiyrMywa5dPQoIt+UhIcBg82KehIYEhQwwFBXDFFQ5r1kBqqkdycsDevS6e5zJsGBw7\nBgMG2HUEQYgtWxw8L8LRowFHjhj27TOMHu2Tng41NQ6u6zF2rGHoUMO8eQEvvuiwf78hLw/q6uxI\nebt7bWvDbfmJHWcPkJRkh/hEa9PDYbeb8hzD2Jw2vvvT33HseIRPzZvCXdcUxY7xPBvoRUQudgre\nIiK9EIkEHaHb6SiLgFcqbFgM/KDjmGhpib1PtCY8HHaZM8entdXn5z+HUCggOxvy820Lv64cx9Da\n6tDWFp1I6TBjxhG2bUsiNXUIjgNr1xreestw4IBLbq7PjBkBhw87bN0Kx48bxo61ZSRVVbbryMaN\nBmNg/HiPMWMcBg922L7dYfZsn3XrHJKTITPT58UX7c52To7puDDSPj50dhyJRALKygI2bTLk5QW0\nthqGDYMZM0yXDy1Bx7OxYTwUcvjt2v9m2759XDYmg8cfuKnLM7b17NGLUkVELmYK3iIi3ehadgJQ\nUWFwXVuLHO0x3dZuv19XB83NEVavdvA8n02b7A5wSYnTEV7pCOCQn28vWiwtNR3nBfD52c8C6uth\n9uyA++4L2LzZAQzt7R7NzR6zZu0nP38Ia9e6QMDw4eD70NZmy1lGj/Zob3fxPBg/3nDllaaj/tol\nLy/AGJ+hQx2GD3fwPJf8fJgyxWHTJtt7e+VK20HlyisDEhIMlZW2DnzKFI+yMhP7wBEN1a5rcJwg\nVvMdremORGyf8vZ2ryOMh/jlX9bzkz+uISnB5YV//jhpyYkAXcp6Ont2i4hczBS8ReSid6p2gL25\nj/26s5SiO0EQvejQinYzOVEo5FBaao+LlmR08mlsNLzzjq2ZjkQgKamdbdscUlNTcByHdesMubnt\nFBe7TJnisHKlT0NDZ825MQEFBQEJCSFWrQo6eom7zJrlc/nlAb/5jWH/fodFi2yv8UgkID8/gu/7\n7NljSE+3z9N2KrHlH6tXQ3W1vQhz2TLbuzs/P+Cee2x7xO3bA07keT41NTb4DxxZz6e/9yIA//6F\nG7k8Lyvu9eiN1lYv9pqJiFzIFLxF5KLWUzvA093H1h7b+ubiYrsrXFFhKCnxCIddkhLtuVLTAior\n3Y52gU5HDfjJgTS6G2xLN7xYS76yMkNBgUNhoUdZmUNtrWHOnAivvRad+OiwY0cY3/eproZt22xd\ndlWVYeBAW3vd2gotLYYdOyAvz2PFCofMTJ+774YVK+DNN+HoUcjPj64lYNkyn7/+1ZCY6DBxosex\nYw5ZWdGdcsOUKT6VlYbCQiguNqxcaT9g1NXZDxclJSY2CCi64x39cOG60HSslc8s/iUtre0suH4y\nX7zjyjP++bW2ejzzjH0tP/lJT+FbRC5oKqgTETmBreeO7y0d7cu9aZOtcY5EfIYMSAag8eix2DHh\nsMvkyR7Hj9uQ3drqxe1st7dHeO89j2XL/I6uHzawOo7D9OkJFBcbZsywpSGHD9td6NGjW8nKOsqV\nVwY0N8O2bbBrFyQk2B3qsWN9cnI8RozwyM6GyZMNbW2wZw+0tfm0t9sa8GHD4KMftR8gysp8PM8O\n6hk7FpqbHRobDZ5na7rBZ/VqWL/e5+hRu87SUsM999i67k2b7POZM8dhzpz4aZPhsMuMGQGP//l3\n1Ow+xNT8LJ76h9swpvu/CbC18erlLSIXP+14i8hFq7splL25T2c9d9edXLejO4ktsYhEAoYPTgOg\nqa0pNnCnuTnC44/DgQNw660RKircjjXYCy937AjYuhUSEw1XXRXEDdqBgN27DcZ4bNwYYuBAn2uu\nMQwadBCA3NxRvP56wJEjMGaM3QmPRAxHjhhaWgKSkmDsWJfkZENWVoQjR3x++1tDbq7Lxz4G4JKY\nCLW1fsf54EMfMhw9CrW1hoSEgCDwWb/erqaw0A7eOXAAkpI8XNcFXMaP93GcznHw3fnmM0t57e0a\nhg5M4Xf/ejeJITfWXvDE1zv6NxIlJd5J5wyHXT75SZWaiMjFQcFbRC5K3U2hfD+6thcMhQxz5pjY\nUJy2iA3eew404Ti2b7XnBXie7Yk9bRpUVgbYHWQb1iMRO5q9sLAz1Fv2QsVIxGCMw8CBAeBSUuJQ\nVeWzbt1AGhoCEhMDHMeQn+9TXe2QlOSyYIEN0q5rOvpnBwweDJs3O1RVge8HzJ/vsGoVlJf7HDhg\nx7pnZzu4rkNSUkBurr29oQH274cgcLjzTjs4Z//+rus0XHWV02Po/uHvVvJvv1iO6xh+tfhORg1N\nP225j+edXIYTpcAtIhcLBW8RkS6ifbfLygIqKjo7m7z5pv23LauwQTQxFCIjPZkDR46xe38TNRvT\nAJevfCWC6xpSU11cN9pqMGD1aodx4wLuvhtSU93YDnl7u0dVlUNzs8/ttwekpYU6WhcawmEH3w+x\nc2cKO3YYrrzS48gRw09+4hIEcM01tpuK67pMnhxhxQqfhASHjIyAjAzbnjAx0cfzDDU10NYWYccO\nB8cxXHEFpKcbWlsDli0LWLnSYcQIn1GjbCvD5GSXCROgvT1g+nRDONy5+38qP/r9Kr782J8A+OHf\n38w1U3J6LCPpnPBpqKggNmhozpze1eOLiFxIFLxF5KJ04hTKU+na8aTza9PlAkm7g1xdbY8vLbXt\nBKPnzv9lBgc27GRd7V6SScPz/I6Wf6bjYsvOMhKwQTgpye4Wt7Z6/PznAfX1hqFDI7z0kuGtt+C6\n63ySkkKUlAQsXw5bt6aSn3+c2lrYuzfEuHF2+E12dsCxY7B9O+TlRXj2WY933jEMGeIzfHhAOGyY\nODHC4cMua9Z4HD5sqKtzSUqCoiJISbEBeuXKgJoaKCryGD8eIITj2NBfWmp391evNnF/c9Bdp5gn\nX3qbL/3HKwA8/sBNfObmK2LlJSUlXkd3mO4nfIZCxJXyiIhcjBS8ReSi1V3gPjFod60vjrblKykJ\nurQHtOUbhYXRc3Z27wCYMXEUKzbsZNXmXTw0f1xs5xZs/XZnfbkbe4yKCtOlR7ghKwtuvRW2bbOl\nFsbEd0QxJiA3t5nGRjDGMHOmw/Tpdmf4V7+yu+of+pDPunXRXtrRke2GY8cctm83JCfDkCFw5Igd\nWHP11S7hsEtrq0dtLTQ1BRw/bvt+L1hgQ3fnyPvoaxcQidgAvXKlDcjR0pH//fJqPv///hGwbQM/\nf9u0uNe2rMx+eCksDE5ZbhIdNHSqn52IyIVOwVtELhknthaM/178tMWuu6497Z6XTBgNVFCxcVcs\nqHZ2RImvg+5a0x19vAULosN1knnoobaOCzud2KCeWbN8kpObCYIQ4MR6hq9aZc9fWOhTVxewYYPL\n3XdHSEkJSEhw+MhHHN5+27B1q6GoCHJy7ACfggKYMqXrBwhDbm5AJAK7djmx26Lr7iwFCTpKQWwd\ne1KSITfXBxx++sc13PvvLwPwhRuuZ8rQK7vpfR50DBM6dT/07l5fEZGLiYK3iFyy4kOlnTpZWkpc\nKUnXANqdGRNHA7Biw07a2jwcx8RCuw2fftw5oiUXZWV07AAbZs2yYTQ9PURKSudud3R3fv36Aezd\nm8yAAT65uYbKSsObb8Lo0QFjxwZs3mxobPRpaDDU1NhR7q7rkJBguOwy0zFB0nTs6Pv8/OcBrhvE\nuoW4rkNhoeHOOwMSEw3xF33SZec76Pgz5OVBSYnh2dff5XP//gcAvvPZ6xncXEJ1NZSWdnaTsWHb\nkJ0ddHSKUbgWkUuTgreIXHDezyTK6PHdB+poV5D4Ee+9eczs4QPJyRrEtvrD/PgXO/jigmxKSgKO\nH/cpK7O7vnPmdJa2RDt3+L7fJeB2Bl3P8/F9j7/+1ZaCFBfbIT6RiCE11cR6d2dn23aAvm9LVcaN\nMxw4AImJPjt3Gp5/3iccNh073U4s9B896rNli52MGS0Zqa62u9GeZ/D9AMfxY2s+8bUrLjasWGGP\nfe6N9/jsoy8RBPBv917PV++cyfLl0XDeuWseidgpmK5L7PXtzWsrInKxUfAWkQvK+5lE2VV3x5/u\nQsyeHtMYw13XTOR7z5fx53XruTcylooK2LbNZ+3agIwMQ3ExDBwYvU/Apk0BubkwejR4XtfHCdiy\nJaC62pCcHDBzZoSVKw2eF3D99QdoaRnKpk22NOTjHzckJblUVMCkSQHFxXZYTksL1NXZloDG2OE4\n0dpsWyrikJMTMGuWHRvvujB+vIfTsbyaGvvv0lIbxruGYlsSAwkJsKTyXb77og3d3/nch/nGJ2YR\nifix3ufdhfbuXt8P+vMUEbmQKHiLiPDBdlsXXP8hvvd8GcurNhAwFzBs22YvcvQ8WLUKrrnGBs8p\nU+yAmro6u2sdHcYTidgEPm4c7NxpyzISEgKqq0MkJoLjtDFlSkBtbUB5uUNSki3niFq5EmxNucuA\nAT6hUMDIkXD55XaXvb3d7xgFD65LbBfcBmLb3SQSCTq6uVitrZ0XnMZCsYE/bvwz3/99GQCPfOZa\n/p/5V3UJ0KZLiUnfvL4iIhcLBW8RuaD0tk1gXz9mSUl0a9qcNIFx8rhMJmYPY8P2fby8oorbZhTS\n3u5w/LgNuwkJbkddt73wsqjIdg8pLrY14BUVdhR9fj5cdZXLlVface3V1S75+QGpqS34fojnnzcc\nPBgwfnxnWPa8AN+3F1S2twfk5QE4eF7A0aOGX/0qYPduj3DYUFrqM368XU/X59b5NbEa9GhJjO0n\nbo/Ze6iFe771a/7yzjZcx/D9z9/Ag3fO/MCvbX//PEVEzhUFbxG54PR3QOtam+15fqyNXtf1fPlv\nruSL/98r/Mszf+X20kKuusqhosKus6QkIBIJeOstnx07YPp0n3vuMaxe7dLe7sXqtn3f/jspyQEC\nCgttmcqvf52C49iykCFDHO64A1JTnVhnFMcxTJ7ssXatg+cZ7r7brmvlSlu2Mnx4QEoKuK6L7we0\nt/tEIk5cvXUk4scCt+f5VFfbtSxYAOEwrNq8izsX/4pd+5sYPjiVF/7541w9JSfuNXi/AVqBW0Qu\nFQreInLJ+iAX9XXuCtsyjE/fNJXv/HwZ67ft4zs/2cj/+PQEug7P8byArKyAI0eI3e55Htu324sk\n77zTp7LSUFZmJ1a6rqG4OIhdgOl5QUfXEXuxZPTxXdd+GPj1r8F1ffLy3NiEydLSIFaOEg3pK1dC\nba2t046Wj0TLROzuti1/ycnxcd2ApCSH/3xpFV/90WtEPJ9JY8aw+M6PMWvSgJNeEwVoEZGeKXiL\nyCXpTC7q69zNNXSG6fiWe67j8PW7Z/HAE3/i//xlKV/7VAGzZoW69L+297nxxoBw2KGy0nD8uM/x\n4+D7sHatobYWsrNt8AY6htQYJkw4yubNafzhD5CT4+N5dJSXRNcTsHOnYfhwOspXTKze2k7PhNWr\n7S667wdEIvFrj0Q6S0qix1dUOBxra2PRd/7I839ZD8ADHy3hI0UfJuS6cXXpJ15MKSIi3VPwFhHp\nhe6CZbQe2k6rNEwYNJXxI1dQs/sA//Kzv/K9z99A1/7XnmdiFzhu2+bT3GxDq+P4bNpkb58xwxAO\nm1hvcQBjbEvBigqoqwuYOtVj+3aHwsKA4uKAVatsS8GCAjtxslPQ0dKQjosmDVdd5XDVVZ1hOVpG\n47rEhvZEIj51+w/wTy/8im379pEaTuCnX/8Id183qUtJCmzZYj84FBaak9oPiojIyRS8ReSSdGJN\n8omj5KNfn0ok4lNWZqdGjhkT4PsOCQkOP/3a7Xz4a8/w6C/LuGVGPldPyYn1v7adR6C1FVpaDBMm\nRHe3bdg1pnMnOhSyu9WhkKGyMkJeXgsHDwIYkpN9kpIC2tvtMJ2aGkNhYcD06YbWVr9jVL3t1d21\nzjw6xr7752U6enn7/O8/ruZ//PS/aTrWRlZ6Bv9r/sf42JzhXV6Tzg8TIiLSewreInLJ6jrgJVp2\nUlLSTQu9bkQith/3xo0+9fUwY4a9aDHhyBgemn8Vj/zXW/ztt39H+Q8/w+hh6YTD0V1nG4ALC6G0\n1CEUMh311y6e57N69YkXcUIQOGzcGKaoCFw3oK7OYEy09Z9DYSEUF9uLKZcuhVGjAhYssLvfXadH\ndn2O0QDe9QNI2fo6Hnj8T7xb0wDA1RMn8DcTbiFnWNJJr9usWT4lJSePmBcRkVNT8BYReR9CIcP4\n8Ybdu6G93dZaO44N7P+4YDZvvLONsvV1XP8PP+PN//g7MgennlQn3jk9k47R9U5srPumTTZkl5TY\n4G13xg2bNhl27QoYNSogIcGNDayx5+ycxrNypQ3ms2bF7+J7nn/SB4v6Q818/cev8/wb6wAYO3wg\nP/jCjdw8vQDP6+z5feIwnRO7onR9TiIicjIFbxG55MWXnbi9aosXCjlcfXVAOGwnS5aUOLEOJKFQ\nAr//1t18+Gs/Y+3Wvdz49Wd54wefZEh6co/n6wzgJjbRMhIJWL/edhApLrYXVXqeIT/fUFIC4XBn\nT+45c2yrQs8LqKx0Tjp/9NzR4N3aFuGx5yv4Xz9/i6Ot7YQTQ3zjnlL+4c5SUsIJXS4+jQ7TATsg\nJ/5vAjR9UkSkd/R/RxEROndxT/z6RPbiQhvMw2GXWbMc5sxxCYfd2P0iEZ8N7yXzLx9dQMHoIbxb\n00DJF3/C+tq9Pa6h63lLSw3Tp8Pq1YadO20f79WrDY7jsGiRiT3myfd1qKx0Oz4MBCeNbg+HXaYW\nt7O6fiVFf/cj/udP3+Boazsfv3oi6/7PF7hu3NW8szoh1n4QAlpbfcrKbGnNB63t7vr6iYhcarTj\nLSLSS93t7Pa0uzskLY3Xvr+Qv/nnF3inup7pX/wJj315Hp++eSqeZwPs6erMW1t9EhLsEBuwPbZD\nofjwG3/fzmE8kUh8CcihpmP88Per+I/fVLC/8SgARTnDeOz+m7juilwiEZ+dW+g43l7cWVYWrTuH\n/Hw6Slu6303v+nx6+/qJiFxKFLxF5JLS21rkD1KzHB9EB7HssU/z+R+8zLOvv8dnH/0DP3nlHf52\nxvVMGjsmdlzXHebjx214bm/3CYfBmICjR32mTwc7Yt6uqbg4Ejdgx4ZyQ3GxzwsvQE1NQF5eO1X1\nDazZ8y7Pv7GO5mNtABSOHMmC2bP4xqcvIzHR7WbdDpGIR3W1LaVZsACSkjrLabp7ziIi0rM+C95P\nPfUUv/jFL3jnnXc4cuQI27ZtY+zYsXHHHDp0iAceeIA//OEPAHzkIx/h8ccfZ+DAgX21DBGRU+rt\njuupjuupBWHX+0Yi8W37UsIJ/Owf/4YPX5HLVx5/nRUbdrJiw/9l1mUF/M/E6YRbcwiCAPDZsgXW\nr7dBOikp6JgieZycnGOxITie5+P7AS+8ABCQnQ0JCU5HaYmL6/rsOLCfd+o2svKP62loOhhb3/VX\njOMb95RyzZQcjDm5G0l8aYph/HifmhqorDQd4f7kGu/e+iBj5UVELgZ9FryPHTvGvHnzuOOOO3jw\nwQe7PWb+/Pns3LmTJUuWEAQBn/3sZ1m4cCEvvfRSXy1DRKTX3s+udnelIV1Hr7/5pkdtLeTlOScF\n1AXXf4hM5zKeX17OrytWsHxzFTf/YxWZA9P58KRJzJ5QgCETsDvQrmvIzfWprXXYvDmF0aNti8Hp\n0wM8D2pqDKFQQBAE1B9u5Lk3tvOXym38efVW9hxsjj3ugKQU5k2dzEOfmsKR+kw4Csac/nmHQg5X\nXRWQkBC95YP37lbgFpFLWZ8F77//+78H4O233+72+xs3bmTJkiUsX76ckpISAJ588klmz55NVVUV\nBQUFfbUUEZFudd1xBU65+93ZAcSWfPSW7e0N9fWQl3fyfUMhhxuuS+KG667hkQPF/OMPV/H6lM7f\nUwAAEddJREFUe+toaGzkF8vL+MXyMhJCDpNzhjN60AiGpKax9kCY6g2GjEEOc4cf5VhLOy/+uJ2m\nY620sJ+qnQeo+u1+mo62xT1WejiVGQXjueFDE0lpy6Eg32FCtkPZrujz710AtheQRu8T3wZRRETO\nTL/VeJeXl5OWlsbMmTNjt5WWlpKamkp5ebmCt4j0i6471qfT0yCdU5VNFBZCQQGUlHQ/VCZ625jh\nA3n6n6/D96+lfEMdz/15LW+t3cGmHftZs2UPa9gTf8ft8Oy7p15rRnoyMyaO5rqpuQwxubQeyCAv\nDxISDFu2BDgdS7H9wM+MgraISN/ot+BdX1/PsGHD4m4zxpCZmUl9ff0p73eqHXS5eOlnLv0lJcUG\nysrKk0N4EDjU1qZ3HHcEY3oO6kHgsH59OkFgKCpqYt26yEnfB7o9TyrwuatH8LmrR9Dc2s6mnY1s\n3tXE4ZY2Go+2cbilnZbWCG3Hk0hJDDF8KKQmuYwZmkx2ZjrHD40gPRymqMiuMwh2E2TXd3zt4DgD\n2LEjIC2tmdraNFzXkJzchOP0fo0i5zP93pDzRX5+fo/f7zF4P/zww3z729/u8QRLly5lzpw5Z74y\nEZGzoDfhsTfHGONTVHTktMedfL/gpOOjoRyIheNTSQsnMC1vKNPyhp50Dt8PsWFDKsYEXUK2w/r1\n4ZPWHn0MY3wmTWqMfT1xYjPr1w9gw4a0uLWcyRpFROT96TF4P/jggyxatKjHE4wZM6ZXD5SVlcW+\nffvibguCgL1795KVlXXK+02bNq1X55cLX3THQj9zeb/iL3jsvkSiN8ec6WOCPc/Uqd1frBmJ+By1\nbbOZOvXkcpfTrWHVqjWsX5/O6NHZjBtnACfuPKd63FOt99ix7tfS3Rq7e54i5wv93pDzTWNjY4/f\n7zF4Z2RkkJGR0ScLmTlzJs3NzZSXl8fqvMvLy2lpaaG0tLRPHkNEpD/1tj1hdxdrnqorSvT4EwWB\noabGkJ8fUFpq2wZ2PT/Yse4QP0b+VGs58XF6avV3poNvFNJFRLrXZzXe9fX11NfXU1VVBcD69es5\nePAg2dnZDB48mAkTJjBv3jzuu+8+nnrqKYIg4L777uO22247bT2MiEhPuga90/WJPt0xZx4aoy32\nTn/RYteLNbt73DfftLfPmRP/+LbspYljx+xud6ib/3O3tno884xdy4IFEcLhU0/VPNPbz4SmU4qI\nnFqfBe8f//jH/Ou//itgL5q85ZZbMMbw9NNPx8pVnnvuOe6//37mzp0LwO23384TTzzRV0sQkUvQ\n+wl6Zzo4pyfR0e9n0nYwuoauHwBaWz2qq+337I52/PGOE2HWrNO38wuFAlasgISEvgu+GnwjItI3\n+ix4L168mMWLF/d4zKBBg3j22Wf76iFFRM4aWxZycnA9cUfcdXsXRLsLrydOiSws7Pz6TIXDLp/8\npEckQseEy97p7Q5/bwO3QrqIyKn1WztBEZGzoS+DXtda7O56eHe3I37iY/cUZHveqe75eQSBc9rd\n+Ghtd29fj7NVFqLALSLSPQVvEbng9WXQs+fqfSu9ro99phdM9nSuD0LBV0Tk/KTgLSJygp52n0tK\nAkKh7qdSnqgvd5SN8Zk1q3N9fUFlISIi/UvBW0SkG6duqWdiAbi7+3QNsna3+/1dfNmbNZ2JU+28\nK3CLiPQfBW8RkT50YpD1vADfD4hEum8DeKKz0QO7tdWjrCzAdR21+BMROYcUvEVEOpzuwsj3W5ZR\nU2N61d7vbAyqiUR8Kiqguhry831AoVtE5FxR8BYRoXeh90wDdyjkUFoa4LrQF6UmXZ1ZSLetCktK\nVFoiInIuKXiLiGD7dnte0Ou+3L0VDru93ik/Gxc7dp6zdxeEiojI2aPgLSKXPFuOYXDdaNcS95yt\nA3rfevBMAr2IiJx7Ct4iIjGmVxdAnoneloS8n9aDCtQiIhcWBW8RueRdCv2sz0a3FBEROTMK3iIi\nnL1A2ttQfzbD/9kaDS8iImdGwVtE5CzrbdBVIBYRubgpeIuIXOQuhVIaEZELgYK3iMglQIFbROTc\n0/+JRURERET6gYK3iMgFKhLxY91KRETk/KfgLSJylpzNYBztVLJ8OQrfIiIXCAVvEZGzQMFYRERO\npIsrRUQuQOpUIiJy4VHwFhE5C7oLxn09PVKBW0TkwqLgLSJylnQNxpoeKSIi+j+/iIiIiEg/0I63\niMgZej8lI6rJFhERBW8RkTPwQUpGFLhFRC5t+i0gIiIiItIPtOMtInIGVDIiIiLvl4K3iMgZUuAW\nEZH3o09+exw6dIj777+fCRMmkJKSwtixY/niF7/IwYMHTzpu4cKFDBo0iEGDBrFo0SIaGxv7Ygki\nIiIiIue1Pgneu3fvZvfu3Xz/+99n3bp1/Nd//Rdvvvkmn/jEJ+KOmz9/PpWVlSxZsoRXX32VNWvW\nsHDhwr5YgoiIiIjIea1PSk2Kior4zW9+E/vzuHHj+P73v8+tt95Kc3MzaWlpbNy4kSVLlrB8+XJK\nSkoAePLJJ5k9ezZVVVUUFBT0xVJERERERM5LZ61QsbGxkaSkJFJSUgAoLy8nLS2NmTNnxo4pLS0l\nNTWV8vLys7UMEZF+F4n4sV7fIiIiUWcleB8+fJh/+qd/4t5778Vx7EPU19czbNiwuOOMMWRmZlJf\nX382liEi0u+ifb6XL0fhW0RE4vRYavLwww/z7W9/u8cTLF26lDlz5sT+3NzczG233caYMWP43ve+\n94EX+Pbbb3/gc8iFRT9zuZAFgUNtbToAKSlHMKbvwrfeGyLd03tDzhf5+fk9fr/H4P3ggw+yaNGi\nHk8wZsyY2NfNzc3cfPPNOI7Dyy+/TGJiYux7WVlZ7Nu3L+6+QRCwd+9esrKyenwMEZELhTE+RUVH\nYl+LiIhE9Ri8MzIyyMjI6NWJmpqauOmmmzDG8Kc//SlW2x01c+ZMmpubKS8vj9V5l5eX09LSQmlp\n6SnPO23atF49vlz4ojsW+pmLxNN7Q6R7em/I+eZ0bbL7pKtJU1MTN954I01NTfz+97+nqamJpqYm\nwIb3hIQEJkyYwLx587jvvvt46qmnCIKA++67j9tuu+202/IiIiIiIhe6Pgneq1evpqKiAmNMXFtA\nYwx/+ctfYjXgzz33HPfffz9z584F4Pbbb+eJJ57oiyWIiIiIiJzX+iR4X3PNNfj+6WsZBw0axLPP\nPtsXDykiIiIickE5a328RURERESkk4K3iIiIiEg/UPAWEREREekHCt4iIiIiIv1AwVtEREREpB8o\neIuIiIiI9AMFbxERERGRfqDgLSIiIiLSDxS8RURERET6gYK3iIiIiEg/UPAWEREREekHCt4iIiIi\nIv1AwVtEREREpB8oeIuIiIiI9AMFbxERERGRfqDgLSIiIiLSDxS8RURERET6gYK3iIiIiEg/UPAW\nEREREekHCt4iIiIiIv1AwVtEREREpB8oeIuIiIiI9AMFbxERERGRfqDgLSIiIiLSDxS8RURERET6\ngYK3iIiIiEg/UPAWEREREekHCt4iIiIiIv1AwVtEREREpB/0WfD+3Oc+R15eHikpKWRmZnLHHXew\ncePGuGMOHTrEwoULGTRoEIMGDWLRokU0Njb21RJERERERM5bfRa8r7zySp555hk2bdrEkiVLCIKA\n66+/nkgkEjtm/vz5VFZWsmTJEl599VXWrFnDwoUL+2oJIiIiIiLnrVBfnejee++NfT127Fi+9a1v\nMWXKFGpra8nPz2fjxo0sWbKE5cuXU1JSAsCTTz7J7NmzqaqqoqCgoK+WIiIiIiJy3jkrNd4tLS08\n/fTT5Ofnk5ubC0B5eTlpaWnMnDkzdlxpaSmpqamUl5efjWWIiIiIiJw3+jR4/+hHP2LAgAEMGDCA\nl19+mT/+8Y+EQnZTvb6+nmHDhsUdb4whMzOT+vr6vlyGiIiIiMh5p8dSk4cffphvf/vbPZ5g6dKl\nzJkzB4C//du/Ze7cuezevZtHH32Um266iTVr1jBgwID3vUBdfHnpyM/PB/QzFzmR3hsi3dN7Qy40\nPQbvBx98kEWLFvV4gjFjxsS+Tk9PJz09nfHjxzNjxgwGDx7M7373OxYtWkRWVhb79u2Lu28QBOzd\nu5esrKwP8BRERERERM5/PQbvjIwMMjIy3teJfd8nCAI8zwNg5syZNDc3U15eHqvzLi8vp6WlhdLS\n0vf1GCIiIiIiFwoTBEHwQU9SU1PDr3/9a2644QaGDh3Kzp07+e53v8vy5cvZtGlTrLb75ptvZufO\nnTz11FMEQcC9997LuHHjePHFFz/wExEREREROZ/1ycWVSUlJ/PWvf+Wmm24iPz+fe+65h4EDB1Je\nXh53QeVzzz3H5Zdfzty5c5k3bx5Tp07l2Wef7YsliIiIiIic1/pkx1tERERERHp2Vvp4i5yJa665\nBsdx4v6ZP39+3DGHDh1i4cKFDBo0iEGDBrFo0SJdxS6XhB/96Efk5uaSnJzMtGnTWLZs2blekki/\nWrx48Um/I0aOHHnSMaNGjSIlJYVrr72WDRs2nKPVivRMwVvOOWMMn/70p6mvr4/98+STT8YdM3/+\nfCorK1myZAmvvvoqa9asYeHChedoxSL944UXXuArX/kKDz/8MJWVlZSWlnLTTTdRV1d3rpcm0q8K\nCwvjfkesXbs29r1/+7d/4wc/+AFPPPEEq1atIjMzkxtuuIHm5uZzuGKR7qnURM65a6+9lkmTJvH4\n4493+/2NGzdSVFTE8uXLYx1xli9fzuzZs9m0aRMFBQX9uVyRflNSUsKUKVPiPogWFBTw8Y9//LQz\nFkQuFosXL+Y3v/lNXNiOCoKAkSNH8sADD/DQQw8B0NraSmZmJo8++ij33ntvfy9XpEfa8ZbzwvPP\nP8+wYcOYNGkSX//61+N2KsrLy0lLS4uFboDS0lJSU1MpLy8/F8sVOeva2tpYs2YNN954Y9ztN954\nI2VlZedoVSLnxtatWxk1ahTjxo3jE5/4BLW1tQDU1tbS0NAQ9z4Jh8PMmTNH7xM5L/XYx1ukP8yf\nP5+cnBxGjhzJunXreOihh3jvvfdYsmQJAPX19XHdccCWp2RmZlJfX38ulixy1u3fvx/P8xg+fHjc\n7frvXi41M2bM4JlnnqGwsJCGhgYeeeQRSktLWb9+fey90N37ZPfu3ediuSI9UvCWs+Lhhx8+7V+F\nL126lDlz5vC5z30udltRURHjx49n+vTpVFZWMmXKlLO9VBEROY/Nmzcv9vWkSZOYOXMmubm5PPPM\nM5SUlJzyfsaY/lieyBlR8Jaz4sEHH2TRokU9HjNmzJhub7/iiitwXZfq6mqmTJlCVlYW+/btizsm\nCAL27t1LVlZWn61Z5HwydOhQXNeloaEh7vaGhgZGjBhxjlYlcu6lpKRQVFTEli1buOOOOwD7vhg9\nenTsmIaGBv1+kPOSarzlrMjIyKCgoKDHf5KTk7u979q1a/E8LxYuZs6cSXNzc1w9d3l5OS0tLZSW\nlvbL8xHpb4mJiRQXF/Paa6/F3f7666/rv3u5pLW2trJx40ZGjBhBbm4uWVlZce+T1tZWli1bpveJ\nnJfcxYsXLz7Xi5BL19atW3n88cdJS0ujra2NsrIy7r33XrKzs/nWt76FMYZhw4ZRUVHBc889x9Sp\nU6mrq+O+++5jxowZfOlLXzrXT0HkrElPT+eb3/wmI0eOJDk5mUceeYRly5bx9NNPM3DgwHO9PJF+\n8bWvfY1wOIzv+1RVVfHlL3+ZrVu38uSTTzJw4EA8z+O73/0ul112GZ7n8dWvfpWGhgaeeuopEhMT\nz/XyReKo1ETOqcTERN544w0ee+wxmpubGTNmDLfeeivf/OY34+rznnvuOe6//37mzp0LwO23384T\nTzxxrpYt0i/uuusuDhw4wCOPPMKePXuYPHkyr7zyyinLtEQuRrt27eITn/gE+/fvZ9iwYcycOZMV\nK1bE3gff+MY3OHbsGF/60pc4dOgQM2bM4LXXXiM1NfUcr1zkZOrjLSIiIiLSD1TjLSIiIiLSDxS8\nRURERET6gYK3iIiIiEg/UPAWEREREekHCt4iIiIiIv1AwVtEREREpB8oeIuIiIiI9AMFbxERERGR\nfqDgLSIiIiLSD/5/+8wmciNNB54AAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"from numpy.random import multivariate_normal\n",
"from filterpy.stats import (covariance_ellipse, \n",
" plot_covariance_ellipse)\n",
"\n",
"mean = (5, 3)\n",
"P = np.array([[32, 15],\n",
" [15., 40.]])\n",
"\n",
"x,y = multivariate_normal(mean=mean, cov=P, size=2500).T\n",
"plt.scatter(x, y, alpha=0.3, marker='.')\n",
"plt.axis('equal')\n",
"\n",
"plot_covariance_ellipse(mean=mean, cov=P,\n",
" variance=2.**2,\n",
" facecolor='none')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Algorithm"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As I already stated, when the filter is initialized a large number of sigma points are drawn from the initial state ($\\mathbf{x}$) and covariance ($\\mathbf{P}$). From there the algorithm proceeds very similarly to the UKF. During the prediction step the sigma points are passed through the state transition function, and then perturbed by adding a bit of noise to account for the process noise. During the update step the sigma points are translated into measurement space by passing them through the measurement function, they are perturbed by a small amount to account for the measurement noise. The Kalman gain is computed from the \n",
"\n",
"We already mentioned the main difference between the UKF and EnKF - the UKF choses the sigma points deterministically. There is another difference, implied by the algorithm above. With the UKF we generate new sigma points during each predict step, and after passing the points through the nonlinear function we reconstitute them into a mean and covariance by using the *unscented transform*. The EnKF just keeps propagating the originally created sigma points; we only need to compute a mean and covariance as outputs for the filter! \n",
"\n",
"Let's look at the equations for the filter. As usual, I will leave out the typical subscripts and superscripts; I am expressing an algorithm, not mathematical functions. Here $N$ is the number of sigma points, $\\chi$ is the set of sigma points."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Initialize Step"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\\boldsymbol\\chi \\sim \\mathcal{N}(\\mathbf{x}_0, \\mathbf{P}_0)$$\n",
"\n",
"This just says to select the sigma points from the filter's initial mean and covariance. In code this might look like\n",
"\n",
" N = 1000\n",
" sigmas = multivariate_normal(mean=x, cov=P, size=N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Predict Step"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{aligned}\n",
"\\boldsymbol\\chi &= f(\\boldsymbol\\chi, \\mathbf{u}) + v_Q \\\\\n",
"\\mathbf{x} &= \\frac{1}{N} \\sum_1^N \\boldsymbol\\chi\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"That is short and sweet, but perhaps not entirely clear. The first line passes all of the sigma points through a use supplied state transition function and then adds some noise distributed according to the $\\mathbf{Q}$ matrix. In Python we might write"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" for i, s in enumerate(sigmas):\n",
" sigmas[i] = fx(x=s, dt=0.1, u=0.)\n",
"\n",
" sigmas += multivariate_normal(x, Q, N)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The second line computes the mean from the sigmas. In Python we will take advantage of `numpy.mean` to do this very concisely and quickly.\n",
"\n",
" x = np.mean(sigmas, axis=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now optionally compute the covariance of the mean. The algorithm does not need to compute this value, but it is often useful for analysis. The equation is\n",
"\n",
"$$\\mathbf{P} = \\frac{1}{N-1}\\sum_1^N[\\boldsymbol\\chi-\\mathbf{x}^-][\\boldsymbol\\chi-\\mathbf{x}^-]^\\mathsf{T}$$\n",
"\n",
"$\\boldsymbol\\chi-\\mathbf{x}^-$ is a one dimensional vector, so we will use `numpy.outer` to compute the $[\\boldsymbol\\chi-\\mathbf{x}^-][\\boldsymbol\\chi-\\mathbf{x}^-]^\\mathsf{T}$ term. In Python we might write\n",
"\n",
" P = 0\n",
" for s in sigmas:\n",
" P += outer(s-x, s-x)\n",
" P = P / (N-1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Update Step"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the update step we pass the sigma points through the measurement function, compute the mean and covariance of the sigma points, compute the Kalman gain from the covariance, and then update the Kalman state by scaling the residual by the Kalman gain. The equations are\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\boldsymbol\\chi_h &= h(\\boldsymbol\\chi, u)\\\\\n",
"\\mathbf{z}_{mean} &= \\frac{1}{N}\\sum_1^N \\boldsymbol\\chi_h \\\\ \\\\\n",
"\\mathbf{P}_{zz} &= \\frac{1}{N-1}\\sum_1^N [\\boldsymbol\\chi_h - \\mathbf{z}_{mean}][\\boldsymbol\\chi_h - \\mathbf{z}_{mean}]^\\mathsf{T} + \\mathbf{R} \\\\\n",
"\\mathbf{P}_{xz} &= \\frac{1}{N-1}\\sum_1^N [\\boldsymbol\\chi - \\mathbf{x}^-][\\boldsymbol\\chi_h - \\mathbf{z}_{mean}]^\\mathsf{T} \\\\\n",
"\\\\\n",
"\\mathbf{K} &= \\mathbf{P}_{xz} \\mathbf{P}_{zz}^{-1}\\\\ \n",
"\\boldsymbol\\chi & = \\boldsymbol\\chi + \\mathbf{K}[\\mathbf{z} -\\boldsymbol\\chi_h + \\mathbf{v}_R] \\\\ \\\\\n",
"\\mathbf{x} &= \\frac{1}{N} \\sum_1^N \\boldsymbol\\chi \\\\\n",
"\\mathbf{P} &= \\mathbf{P} - \\mathbf{KP}_{zz}\\mathbf{K}^\\mathsf{T}\n",
"\\end{aligned}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is very similar to the linear KF and the UKF. Let's just go line by line.\n",
"\n",
"The first line,\n",
"\n",
"$$\\boldsymbol\\chi_h = h(\\boldsymbol\\chi, u),$$\n",
"\n",
"just passes the sigma points through the measurement function $h$. We name the resulting points $\\chi_h$ to distinguish them from the sigma points. In Python we could write this as\n",
"\n",
" sigmas_h = h(sigmas, u)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The next line computes the mean of the measurement sigmas.\n",
"\n",
"$$\\mathbf{z}_{mean} = \\frac{1}{N}\\sum_1^N \\boldsymbol\\chi_h$$\n",
"\n",
"In Python we can compute that with\n",
"\n",
" z_mean = np.mean(sigmas_h, axis=0)\n",
" \n",
"Now that we have the mean of the measurement sigmas we can compute the covariance for every measurement sigma point, and the *cross variance* for the measurement sigma points vs the sigma points. That is expressed by these two equations\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\mathbf{P}_{zz} &= \\frac{1}{N-1}\\sum_1^N [\\boldsymbol\\chi_h - \\mathbf{z}_{mean}][\\boldsymbol\\chi_h - \\mathbf{z}_{mean}]^\\mathsf{T} + \\mathbf{R} \\\\\n",
"\\mathbf{P}_{xz} &= \\frac{1}{N-1}\\sum_1^N [\\boldsymbol\\chi - \\mathbf{x}^-][\\boldsymbol\\chi_h - \\mathbf{z}_{mean}]^\\mathsf{T}\n",
"\\end{aligned}$$\n",
"\n",
"We can express this in Python with\n",
"\n",
" P_zz = 0\n",
" for sigma in sigmas_h:\n",
" s = sigma - z_mean\n",
" P_zz += outer(s, s)\n",
" P_zz = P_zz / (N-1) + R\n",
"\n",
" P_xz = 0\n",
" for i in range(N):\n",
" P_xz += outer(self.sigmas[i] - self.x, sigmas_h[i] - z_mean)\n",
" P_xz /= N-1"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Computation of the Kalman gain is straightforward $\\mathbf{K} = \\mathbf{P}_{xz} \\mathbf{P}_{zz}^{-1}$.\n",
"\n",
"In Python this is the trivial\n",
"\n",
" K = np.dot(P_xz, inv(P_zz))\n",
"\n",
"Next, we update the sigma points with\n",
"\n",
"$$\\boldsymbol\\chi = \\boldsymbol\\chi + \\mathbf{K}[\\mathbf{z} -\\boldsymbol\\chi_h + \\mathbf{v}_R]$$ \n",
"\n",
"Here $\\mathbf{v}_R$ is the perturbation that we add to the sigmas. In Python we can implement this with\n",
"\n",
" v_r = multivariate_normal([0]*dim_z, R, N)\n",
" for i in range(N):\n",
" sigmas[i] += dot(K, z + v_r[i] - sigmas_h[i])\n",
"\n",
"\n",
"Our final step is recompute the filter's mean and covariance.\n",
"\n",
" x = np.mean(sigmas, axis=0)\n",
" P = self.P - dot3(K, P_zz, K.T)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation and Example"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"I have implemented an EnKF in the `FilterPy` library. It is in many ways a toy. Filtering with a large number of sigma points gives us very slow performance. Furthermore, there are many minor variations on the algorithm in the literature. I wrote this mostly because I was interested in learning a bit about the filter. I have not used it for a real world problem, and I can give no advice on using the filter for the large problems for which it is suited. Therefore I will refine my comments to implementing a very simple filter. I will use it to track an object in one dimension, and compare the output to a linear Kalman filter. This is a filter we have designed many times already in this book, so I will design it with little comment. Our state vector will be\n",
"\n",
"$$\\mathbf{x} = \\begin{bmatrix}x\\\\ \\dot{x}\\end{bmatrix}$$\n",
"\n",
"The state transition function is\n",
"\n",
"$$\\mathbf{F} = \\begin{bmatrix}1&1\\\\0&1\\end{bmatrix}$$\n",
"\n",
"and the measurement function is\n",
"\n",
"$$\\mathbf{H} = \\begin{bmatrix}1&0\\end{bmatrix}$$\n",
"\n",
"The EnKF is designed for nonlinear problems, so instead of using matrices to implement the state transition and measurement functions you will need to supply Python functions. For this problem they can be written as:\n",
"\n",
" def hx(x):\n",
" return np.array([x[0]])\n",
"\n",
" def fx(x, dt):\n",
" return np.dot(F, x)\n",
" \n",
"One final thing: the EnKF code, like the UKF code, uses a single dimension for $\\mathbf{x}$, not a two dimensional column matrix as used by the linear kalman filter code.\n",
"\n",
"Without further ado, here is the code."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAt4AAAEPCAYAAAB1HsNIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd8VFXawPHf9Jlk0ishkAQSSgi9iwFRWkCayCosKqiw\nuoiwrour6yq4oqIs9rXtIigC6uoriqw0qQLSewkQSnqflMnUe+/7x0BgTAIBEhL0fD8fTObec889\nM7kmz5x57nNUiqIoCIIgCIIgCIJQr9QNPQBBEARBEARB+C0QgbcgCIIgCIIg3AAi8BYEQRAEQRCE\nG0AE3oIgCIIgCIJwA4jAWxAEQRAEQRBuAG1DD6A6JSUlDT0EQRAEQRAEQbhmAQEBVbaJGW9BEARB\nEARBuAFE4C0IgiAIgiAIN0CjTDW5VHXT9IJwwa5duwDo1q1bA49EuBmI60WoLXGtCLUlrpUGlpcH\nERFgNkNJCagbdk75SunSYsZbEARBEARBuDnt2+f52rFjgwfdtdH4RygIgiAIgiAI1dm71/O1c2fP\nV1mG3NyGG88ViMBbEARBEARBuDnZbODnB506eWa/W7aEsWMbelQ1EoG3IAiCIAiCcHOaNQssFrj/\nfk/QnZ8PmzfD8eMNPbJqicBbEARBEARBuHmp1aDTeWa+773Xs23BgoYdUw0afVWTK5FlGafT2dDD\nEBpITEwMAHa7vYFHcv30ej3qm+DGEEEQBEFotB56CP7zH1i4EF580ROQNyI3deCtKAoOhwOj0YhK\npWro4QgNwGg0NvQQ6oSiKNjtdnEtC4IgCML16NULEhPhyBH4/nsYNaqhR+Tlpp5eczqd6PV6EagI\nNz2VSoVerxef3giCIAjC9VCpPLPePXuCydTQo6nipg68FUVBo9E09DAEoU5oNBoURWnoYQiCIAjC\nzWHdOjh1Cn75t3PGDNi+HQYPbphxXcZNHXgLgiAIgiAIv0FuN9x5J8THe6qaXOoG3i8lV1hxf/oW\nDkcRFRXnrtheBN6CIAiCIAjCzSU1Fex2iImBoKA67VqR5crvZVlG/m4JktuJ01mK3ZqD++FBlJce\no6RkL5aKHdi0uZSVbqGiYv8V+76pb64UBEEQBEEQfoN+uWLlNbiQ3il9/E/k3z2ArAVZdqKffC+2\n199ENqhQFDu6jB3Y8s2gUwMKqr9ORHGkevLJAcfAXrU+pwi8BUEQBEEQhJvLvn2er506XbGpUpCL\n4h8IWh1udxmqv07G8dg0pBAjbncJuqAKnCWbweApPVjxzz+DOhNcnuPdY++40JPnv4F+1zxskWrS\nyCxcuBC1Wl3jv9WrV19Tfzt27PDaXl5eTnJyMnq9nq+//hqAWbNm1Xje+fPn19lzFARBEARBuC4X\nAu/LzHgrsox7xlTku2/DunMZxcVrKS3dTOljo7H7FeNy5aAoNpzJnSuDbgC09Ve4Q8x4N1KzZ8+m\nZcuWVbZ36NDhuvu2Wq0MHTqUHTt2sGzZMu666y6v/e+++y4BAQFe27p27Xrd5xUEQRAEQagTXbtC\nWVmVwFtRFNzffoJb78beqzW++7eh33gMft6B0nqIp42/uSFGDIjAu9EaPHgwPXr0qPN+LwTdP//8\nM0uXLq0SdAOMGTOG8PDwOj+3IAiCIAhCnXjllcpvFUshctoxHG3icLmykVppUExGFKkY57Be6Dfs\nRb/pAI77hzTceM+rdarJpk2bGDFiBNHR0ajVahYtWlSlzaxZs2jatCk+Pj7079+fI0eOeO13OBxM\nmzaNsLAwzGYzI0eOJDMz8/qfxW+QWq3m0Ucf5ZtvviEpKQmj0UhSUhKrVq2q8ZiKigqGDRvG9u3b\nawy6BUEQBEEQbhay7MaWvR/boRVUVOzG5cpCDgtAMfsA4Ez2ZArofjpYtd53A6h14G21WunQoQNv\nvvkmJpOpymqRc+fOZf78+bzzzjvs3LmT8PBwBg4cSHl5eWWbGTNm8PXXX7Ns2TI2b95MaWkpd955\nJ/IlZVsED4vFQkFBQZV/l9q2bRuPPfYY48eP59VXX8VutzNmzBiKioqq9Ge1Whk2bBjbtm27YtBd\nWFjodc7q+hMEQRAEQWgIiqIg/d8irPnHsFg2UxFmxT70lmrbyi2ikCKDUReUoDl25Trb9a3WqSYp\nKSmkpKQAMHHiRK99iqLwxhtv8PTTTzN69GgAFi1aRHh4OEuWLGHKlCmUlJSwYMECFi5cyB13eO4O\n/fTTT4mJiWHt2rUMGjSojp7Sr8OQIdV/HGK329Hr9QAcO3aMI0eO0KJFCwD69+9Px44dWbp0KVOn\nTvU6btKkSWRlZVWb0/1L7dq183ocGhpKXl7etT4VQRAEQRCEOiHLMnZ7FrKciaN0D0rAFSqMqFS4\nkjug+XID2oNpSG1j6m9wtZhRr5Mc79OnT5Obm+sVPBuNRvr27cvWrVuZMmUKu3fvxuVyebWJjo6m\nbdu2bN26td4D75Xbl/LDz5/XW/9Det7D0F7j6qy/t99+m7Zt21bZrtNdvOu2f//+lUE3QPv27fH3\n9+f06dNVjsvLy8NoNNK8efMrnvvLL78k6JJi9BcCfUEQBEEQhBtFuSSQdW9bjZx1koo7OiFJRdCv\nfa37qXh6AtYXH0YJ9q+PYVYyLliJ9Ei/y7apk8A7JycHgIiICK/t4eHhZGVlVbbRaDSEhIR4tYmI\niCA3N7cuhvGr0r179yveXFldEB0UFERxcXGV7R988AFPPvkkKSkpbNy4kcTExBr7TU5OFjdXCoIg\nCIJwwymKgiS5kHauRbVlFbYpk5CkYpTQQlT+PijF6Zg+/A535wRcd9Su4prc7MbENOq8YqQrtKn3\nqia/zAW/Wrt27apxX0xMDEaj8br6v5lpNNXXmVSq+aijdevWrFq1iv79+zNo0CA2b95MXFxcfQ9R\nuEplZWUcOnSooYfxq3e53yuCcClxrQi1Ja6V66MrLyF69WLOjhmAWyoGdwX0jEI++qNXO79Nu2g3\ndwnlbZtxNKZ+Z7Cv2tieVM1V8FYngXdkZCQAubm5REdHV27Pzc2t3BcZGYkkSRQWFnrNeufk5NC3\nb9+6GMZlDe01rk5TQW5GnTp1YsWKFQwaNIiBAweyefNmmjRp0tDDEgRBEATht+j8RKFW50DyyyO7\ncwROZ45nKXZd9SGqb6qnGp61VdMbNszaKnfYr9imTgLvuLg4IiMjWb16deVCK3a7nS1btjBv3jzA\nswCLTqdj9erVjBvnCYAzMjI4duwYt9xS/Z2oAN26datxn91+5ScoeOvTpw9fffUVI0eOZNCgQWzc\nuJHg4OCGHpZwnp+f32WveeH6XJiREq+xcCXiWhFqS1wr10ZRFNwfzsERH4m9YyRghpg+VzzOnPM9\nAL59OhEfn1DPo7w8t+TmTG4OB0/ns2ZXMAdPdWL7e5c/ptaBt9Vq5cSJE4DnjtKzZ8+yb98+QkJC\naNasGTNmzOCll16iTZs2JCQk8OKLL+Ln58f48eMBCAgI4KGHHmLmzJmEh4cTHBzME088QceOHRkw\nYMC1P+tfqR9++IHU1NQq23v27ElCQs0XWnVpJr80ZMgQFi9ezLhx40hJSWHdunWYzQ23ipMgCIIg\nCL8iNhu89x7cfTdccj+anJ8NORm4W7fBbj+Lq39LlAATUPv62tpDngIS7qQWV2hZDbsT3a5juHq1\nu6Zl4RVFIb/EQlp2JtuPOtlyoDmnMvuQnd8WWbkQUpdeto9aB947d+7k9ttvBzx5288//zzPP/88\nEydOZMGCBcycORObzcbUqVMpLi6mV69erF69Gl9f38o+3njjDbRaLffccw82m40BAwawePHi684D\n/zW58FrMmjWr2n1vv/32ZQPv6l7L6raNHTuW0tJSJk+ezMiRI1m5cmWNbQVBEARBEGrto4/gz39G\nXvM/5Mcm4b5jAJJUipx5ANXR/dib5AIKBF/lpJ/Dheb4ORSVCvc1lAUMHPxntEfOYvnhNdxdW9fq\nmBJrOYdOF7DxgJY9qUGk57Yit2gY2nINy4+N5Nnm8WQFqIgKO0yLpruBUZftT6XUZor0BispKan8\nPiAgoMZ2drv9N31zpfDrI67p+iU+EhZqS1wrQm2Ja+UiudSCvPor3GHhGG8bgRzsR/nHT+G8pUPd\nnMBqx7hkLersQiqee+CqD/d98l+YFv1A8VP3kvXQEMrtFVjtdixlCrnFevItevItJopKjRSV+ZBT\n6M+53FiKy6K9O1IUvjoxmtEF33IiKoYf33mIji1DiQmPxMd0e2Wz6mLYeq9qIgiCIAiCINwAiuK5\nMfFGnU6WUTZ+j7t3X9xyMS5XDmrHKZzNfTECqnIbzu41ly++ar5G7JPvrHVzl9vNxgPprNphJr/E\nSJeMdrzADxz6OI3fHe6JzRFAhT0Qt3T5CS+12kVYYBoRwSeIDj/Do9lrGLbtOyQ/EyHfPMOYuNoX\nqhCBtyAIgiAIQmMhSfDQQxAfD88+W/vjsrKgSxcYMgQWLqy34Sk2K4pGixs7Tmcu6sPrsLWwo/ie\nD14H9wZAah6O5lwemrQspNZXXryvLtkcDj5bV8pna1py7MxQJNmzEOAu5+28wD/pWnCQooI4XGrP\ndo3GgY+hBB+jBZPRgo/Rgo/BgtlURGToabrE20mIjiAuMoombn9Ceq8DwDp/GvJVBN0gAm9BEARB\nEITG4+BBWLTI8/2ECRAbW7vj3nkHcnM9x86fD3VcsUxRFNzuCpT5T2K/vSfOhFBAgXvvqLa9u22s\nJ/A+eu6GBd65xU7e+K+Wb7d0Jd8Se36rTHT4foL8M/AxlHA6rQlxRdn8o+NULJ2iCA9yEWhWMOr1\nGHV6DHodRp0evU6Pj8FAk+DeaDUXw2XfJ/+FuqwC5x1dcY669arHKAJvQRAEQRCExqJTJxg2DL7/\nHv79b3jxxdodt3Pnxe+/+w4euPocaPCkj+B0gMEzg+3evBLJkomjX2dcrjyYPATUaq5UicTdtRXq\n4jIw6q5pHFfjeLrErIXhbN7XEafbBwCjoYTE2HX0StpMSo8IosMiMBtNhJlvxXE8nSl3tcXdpdVV\nn8vVtyP6TfuxvvDgNY1V3FwpCI2IuKbrl7gJSqgtca0ItVUv18rmzdC3L0REQHo66GoRvCoK3H8/\nLF4MI0bA8uWeWeqD21A2r8Rxv2cNFe32nWgOH0b6w5Oo1XrUp06izs5C7jsASaqA1V9BdjqOe0fh\ndpeiys4AnRo5JPCyp9fuOobP3CXYpt+N69Y6upnyCpwuidmf6Fi08rbKgLtJ6FGSWv5Az8TD9OvQ\njqTYlmjU6ro9sSyff/NRlU53cVFIcXOlIAiCIAhCY3frrdC2LRw9CsuXe+phX4lKhfLqq8g4kP2c\n2EoP4nYXIgcWor49AdmRBoCjYxAk9gLrDgDU7gLUmlLcJRsAGXpHgaopuHIAUCJrl7Jievcb9Bv2\n4e7Qsl4Cb936vRiXrMUx/BacI/qwdk85T73fmfRcT1nAuKif6Zm0lA4trNya1JHW0SPrr0TydQTy\nIvAWBEEQBEFoTFQqeOQReOEFuCQLoDqKJCEv/wTHgNtwGLOQXvs9qvJyFOcZTwMfPbLPJcGzWg0G\nfeVDOSoUOSoUkC+e+yqpT2ej/34bik57VVVHaq3Cgc+rS9DtOk5pVCgPHbuVldvGIMs6fIzF9O3y\nIbckpTKwSw9aREY16jVJROAtCIIgCILQ2EyeDFOmwGXSD91uJw5HBootDUexL4qvCbRqlED/GzhQ\nMH3wLSpFwT6mH3JkSN12Lkn4PfpPT9AdHEz/o3/liKM7AIlxa+jfdSkpPdrQLWEk6rpOKakHIvAW\nBEEQBEFoDH7+GZKSwNcXTKbKzYrLCS4nKh+zJ2/7m49xa2xU9GmFothhcI8GG7KquAzj0rUA2B4d\nWSd9Ot0ucooKScuyEDHre5J/2o5F50+f5hs46kgiwJxN/67vMay3nQGdU/A1mq7c6TVSFZRgXLoW\n2+ThYNRf+YArEIG3IAiCIAhCQysvh1tvRdHpkJYvRi7JxT10BOBCteZ7VAW5OO6+C0VxIif5IBv9\n0G3ZiXbrIRy/H3g+XeTG0xw9i2LU4+qZiJQY67VPVVCCbvthFF8Trv6da+xDkiTO5eeyPy2XbYf0\nHDoTS0buLTQ7W8L2/TNxqPSMbPUtx82t6Zzwfwzvs5qRt3SjWVjENY/b8PmP6Dbtxzr3Dyhmnxrb\n+cz9DNPCH9CczKT8zcev+XwXiMBbEARBEAShASmyjLL8S9RuN3LHtpQmyGA1IVfs8TTo0xxoDq50\nz+MgXwCMH3yL4X8/g0GP7fExVTt2uMBQv+X83LckUbR3Aeqi0sptuUUG1uwKJ37PKe785yvkd+/O\nyYRbCTS7CPJzoVErnMiEdXskdh41cDw9jNyinpRaI736zvaV+UPn+QQFZBDUayPPJnxOrzaxdIlP\nue60EuPHK9HtTsUxui+uAV2rbaM5cgbjJ6tRNGpsU0df1/kuEIG3IAiCIAhCA1EUBUdBGppv/4MK\ncPRphWw2gNlw2ePU2YXoV+9E0Wqw33u7985yG/4TX0J75AxF+z8GXT2Hez4GnIZwVv4UwH++j+Tn\nI3HIioY4ux938idc+zLoM7V/ZXOtxolbqpq2oVa7CAk4R1ToURJj0+jZrpCEKF+iQtoRERSEVlPz\nrPnVcvVpj253KrqfDlQfeCsKvn//DypZxvbwnUitmtXJeUXgLQiCIAiCUMfkfdtRmsej+PsDCixf\njNI9GZpEAwqqZR/hTO6DI8CJpC4i4HguKsCVXLUUn+G/GzB+8C3l78yoXAXSsGQNKknGMfwWlPAg\n7wPMJtRZhajzS9BtO4yrb8d6e577Tkp8tCKUVTsSKbV6xqFWuYmL2oFBW4Z1v4koVzZxhkNkK81x\nuHxxS3qMhhLCAk8TGniG0MAzxDXJpVeiljbNmtI8LAK9rgXQot7G7bq1A7z1FbotB6vdb/h0FfpN\n+5EDzVT85d7L9rVu3QF8fCIZOHAEDsflzysC70Zm4cKFPPjgg2zfvp0ePS7eLFFeXk5KSgo///wz\nS5cu5eDBg7zwwgvV9jFv3jyeeOKJGzVkQRAEQfhNurAGoUqlwrXkbaR2SUjxMUiSFe2BlThU8UhN\nPbnXOv9i3M7dKCUnANDGqJBcqSiSPypLOdqDaSg6La4eiVXOo91+BN2+kxg/WYV1zmSQJIyL1wBg\nf2BItWNzDuuN9o0v0a/YVqeBt6XcxZpdEhv2+7HzaDPOZF8cb4A5i3Yt1tA6ZgO+JgsA+VsN+J6z\nMb3do5yMD0VRVLjcBvQ6B83CwmnVtBkJTZsTFtDyhpYBdPVoi6LTott3ElVJOUqA2buByfOJQ8Vf\nxqEEe1eJcbslcnIsxMS0wmCIIiGhCSaTL0ZjGA7H5cs/isD7JmC1Whk6dCg7duxg2bJl3HXXXRw8\n6HmH9u6771ZZGalr1+pzlQRBEARBuDaKooDbDVotkuRE/vg1pKgw1E38kKQSylvpkc05KLZyAJxD\nu58/UgLA1bn1+ceeetnuxIuzuaqyCpx39ga3BD5VU0zs9w/GtOgHDJ//iPVv96M9cApNRj5SbGS1\nM+QAjjt74/PGl+hXbsf6ypRrXvSloETmhx2wab8f+042IT0vDsWtRa3IuNU6NBoH8dFbSWyxlqZh\nR2keHkHr6Fa0atocX6ORgO2lcG49Y0NakD/kNlyShCzLRAQG4/OLUon6VTtQ9DpcvdvVSQWRy/I1\n4u7SCnVWAariqoG3K7kDpR/+BeeoW722q1R6UlNL+frrXbz++gOoVCq6dImt9WlF4N3IXQi6L8x0\n33XXXV77x4wZQ3h4eAONThAEQRB+nRSbFcVajiokHJerHOW//0aWbdiH9UOSSlANaoViMuBMOwWA\nFJ9wzeeSm4VT9p+natwvdWiJq0srdHtSMSzfgmPcHRRvfAt1XnGNAbXUoSVSs3A06Xlod6fi7t6m\nShubQ01usZH0PA0nMyEtW8u5XD3ZhT7kWfwoLg3E5qhaE/zx0leZnv4OHw4YSkkfG+1igmgd3ZxW\nTcdXCaaVAd2x6fX4dGpNVEjYZV8Hnzmfoj16lpKv/4Eruf7SYy6wPToS03vLkWMjq+yTI0Nwjk4G\nPDPcL7/8NbNnP4vZHENysp6+fcde0zlF4N2IVVRUMGzYMLZv315t0C0IgiAIQt1Q8nNQ8jKRWrXB\n7S6Fjd+j5GVgG5aMotigfyxoNSAVedr71l/t6OrYHxiCbk8qxk9+wDHuDqTE2Crl+7yoVNhSemP4\nfD0710h8fSCWMzk6cooM5Ft8KC4zY3PUXEbvArXaRWjgGaJCDxMVdpQ4/4P845/f4ldawUM9XPiO\nH41OW3M46Rx5K86Rt9a4v/I8Gfloj55F9jXh6lk13eZ6HTp0hsjIYEJDPW8kjh/PIPyWJIKG9Qbg\n2LF0wsMDCQ72A2D//tM0bRpOkyZtMJubkpJixGSKQau9/E2vVyIC70bKarUybNgwtm3bdtmgu7Cw\n0KukjlqtJjg4uNq2giAIgiBcpJQWo5j9cToLcWf8DMf3Yo/IAWToFgFEgFLhaVzflUGuwDEqGd+/\n/wd1vgVVcRlKkJ/XflmGnw6FsPt4APtOmjhy1g9LZh/KEv2Qfq5+7Gq1C19jMT7GYsw+hZh9CvAz\nFXi++ni++ppKaBIcRGxEFLERTWj3qYyfpQJXp3gCJ4285hSWX9Kv2w2A67aOoL/+EoiKomCzuQgI\naIHRGMPevcdITo4hLi4JUNiy5V369WtL8+btARXbtr1P375JxMZ6ZtoPH04lKKgl/v4dUalUDBlS\nfS791frNBN7qPkq99i//VLc3BEyaNImsrKzKnO6atGvXzutxaGgoeXl5dToWQRAEQfg1kWUZZ9ZR\n1K8+Q9lz01CwQTM1NOvKhRzsRsfHQMmq15BaRHkFu1a7hsWrw/ng21gy8n+xXLtaJsCcS0jAWYL9\nz+FvzsXXVISvsRhfUxFGfTkqlYJGrSbI7EeQnz/Bfv4Em/0J8mtBsF8nAnzNaM6fT5VbjN9733rO\nO+vBOgu6AXRrdwHgvKNbnfS3adNptm1L47XX7kKlUvGnP/3Fa/9TTz3r9fgvf3nG6/Hjj/+5Tsbx\nS7+ZwPtmk5eXh9FopHnz5pdt9+WXXxIUdLGMkF5fzzcjCIIgCEIjpRTloQSEoACK4kZ5ZzbypMdQ\nDFpk2YV+6v3YXpmLS1eKbLTCsw8CtoYedq1J8dGV35/J0fH6FyEs39KWCodnQR2zTz4tm24nJOAs\nIQHnCA44h07rqW/nazQRFhBIoNmPQN8QAs2xBPr6EWg242fyqVVFEZ9Xl6CqsOMY0gN3n6S6e2IO\nF/pNBwBqXMymNqxWO35+wRgMCQwfPpgRI1Q3tFJKbfxmAu+6npGubx988AFPPvkkKSkpbNy4kcTE\n6vOdkpOTxc2VgiAIwm+SknkWJSQMSQNutwXNP6Zhm/YQLn81IKNLDMRVsR1cntSFin88CuoskM/H\nBHU4Y3utfJ/5CHebZjh+d/sVK3mUVjhYu1vDgu9j2XmsPYqiASAi5DidWn1Ly6bbUatl/Ew+NAkJ\npUlQIpHBoUQFh2CuZXB9Oa6+HdFvOUjFcxOvq58qZJnyfzyI9lg6cpOQK7evloYZMxbyyivzaN36\n8pOWDek3E3jfbFq3bs2qVavo378/gwYNYvPmzcTFxTX0sARBEAThhlAUBWnFUujQBSUqGkWRUH36\nHlLvvihxLVEUJ5rPXsM+pA+uKH9Ahr9NxFO+73wJvy5tvTs1Xt+NcXVNnVWA6aPvkM0mHOMHVm53\nud3klxSTW1zMwTQV2w434WBaS85kJ+J0eWa3VSqJhGab6dhqBZEhqTQLi6BDXG9aRzfHbLryTZPX\nwjnyVpzDb7nqNyzqc7kYvlyPEmDG/vCdVRuYDDjuH8IV1p6plizL6PUR+Pi0ZenS2zH+oqpKYyMC\n70asU6dOrFixgkGDBjFw4EA2b95MkyZNGnpYgiAIglDnPIH2EtwRwbjaxuJ2F6E2ZCLZZeSSkwBo\nOgUh+6SjWAs9Bz04+PzRjTQv+wourJro7t0ORaPmbG42u0+ksnJ7KKcyupCeO4yyCu9PtQPM2bRs\nuo328f8jJtJJh7h42seOJcivatm/C1QFJRi++wnF3xfHmH7XN+hr+JRAnVuM7ytLcCfFVR94X6Oz\nZwuYN28Fixd/gUajqbN+65MIvBu5Pn368NVXXzFy5EgGDRrExo0bRdUSQRAE4VdB2r4Ot7UQV89O\nuN1FKJEVKP4yssPt2d/lQu1pT4EEKSaqgUZas+Zvf4fKLcE7T8JVpnLoNnvymlMTwvhi+df8fKQj\nu44+TnFps8o2Rn0p0REHaBaxn9gmh4lv6qRZWARJsb1pGhJWq/QR3a5jmGe+j2LQoWg0VRaFqW9S\nG0/qh+ZEhmeRIO31B8labShJScnMmzfopgm6oQ4Db7fbzXPPPceyZcvIzs6mSZMm/P73v2fWrFle\nL8isWbP46KOPKC4upmfPnrz77rs15i8LHkOGDGHx4sWMGzeOlJQU1q1b19BDEgRBEITLk2VYtgz6\n9YOmTQFQiguQz5zA1TYepzMPyZSNYrQj2497jomrupBJY6bdcZTQJRsAKLl3UK0XfbE5HJzLy6HD\n2h0YgbeKO/DFF9MptXqev59PHt0TN9K9zTm6JFiJDA4iPDCYIPMQrxLCteUc2A37mH4Yv9qI/+RX\nse06hvX5iTesRKLi51O5mI/mdDZSQvSVD6rBunX70WiacNddQ1GrNTRtGliHI61/dfaKv/TSS3zw\nwQd88skntG/fnv379zNx4kQMBgPPPusp2TJ37lzmz5/PokWLaNWqFS+88AIDBw7k+PHjmM3mK5zh\nt6O6d69jx46ltLSUyZMnM3LkSHr06NHo7tQVBEEQBADWrYOZM1H27IGJE1D+s8hTKzt3L6oD67E1\nu8XTrmkAENCgQ70epn/9X+X3hi82VBt4y7JMfomFzMI8MvLzyCjIp6DUQkCukz4FZRTpgvjPmTko\nKjVBflnc1e8nHh1ZQfPwEKB1lf6uiUZD+XtP4O7SCt/nF2D64Fu0+05SuugZlJCaU1Q0JzPQ/7AD\n20PDwFT1mR1FAAAgAElEQVT7/Pjvv99Fjx7dCA83oCgOTxpR2xhP4H3kjHfgrSiX/aTA6XSRkVFA\nixbRaLVBtG07CIMhALX65pnlvlSdBd47d+5kxIgRDBs2DIDmzZtz55138vPPPwOe3K033niDp59+\nmtGjRwOwaNEiwsPDWbJkCVOmTKmrodzUJk6cyMSJE6vd99BDD/HQQw9VPn755Zdv0KgEQRAE4cqU\nvXtRnnkS9Q8/AiBHRyJLmZQVr0dWKiAUGHZLww6yjqjTstCv9MQ4p/80Cs2MeykuyKe4vJTi8jKK\ny0opKi8jp6gQh8uNpSyKvOJ48oqSySuOp6wgkt1ttxPuyiMiJJtJw/bxyHAnJoMJqIdVMVUq7FOG\n4+7YEr+HXwVJQvGreh5VYSmG5Vsw/HcDup3HAE8ueumy52vsOjOzELVaS/PmbdHpwtDpSjAY2hEY\n2Byns4gRI+7h3y1jiWEn2mPncI48f2CFg8Bej1DSKxHlvSdAo6GiwsHatfsYPXowGk0AFksF77//\nDR9++BEajZ4OHW7uScc6C7xTUlKYO3cux48fp3Xr1hw5coT169fzzDOeguSnT58mNzeXQYMGVR5j\nNBrp27cvW7duFYG3IAiCINyEFJcLRa3GnZWGtmdPlB4xSAG+2B8bg23KcPAxXFz98QKnq05WJ2wo\nJdZyis6eIaprC3KNsDhOjfXz5VhtwZTbQrDag7FWxGK1hVBY2oz84ha43FUrjaS2a8/QsSd59ba9\n5+9ZrP/XxN0zEcu611FJcpWfgWHZOsx/eseTsw4oPkYcw3pRMXO8V7ucnCLKyx0kJMSi1Yaxc+dp\nIiJiSErqgkqlYvz4CZVtTaZwli79kuATJ5Ba98DWpQ333fc6n3/+Jtr1W7DnFDF0Xwbr/bqjUmkx\nGBQyM48QENAHtVqNvz8sWNCr3l+XG6XOAu8//vGPZGRk0LZtW7RaLW63m2effZZHHnkEgJycHAAi\nIiK8jgsPDycrK6vGfnft2lXjvpiYmEZfNkYQrkZZWRmHDh1q6GH86l3u94ogXEpcK5enUkHLJfPI\n6dcNS7iWmFE9QaUi89n7cAeaIetclWOCNhwk5p3vOPr2IziaNM5iAbKi4JLcON0X/rmw2CooKC+l\noKyUCqeLvKJ4znR6kfTsjhR8GYtbunw8EmAupEVkJonN8+ncooh2zQsI8LUDkJZ2I55VNU4WeT00\nBevpoChk9mzNmb7t0aV0A18fjh4+x9mVBxg+fCBgYNcuC3l5FnS6riiKRKdOnk8xdu/eXeOpzuh0\nkHwbsizzyCOPk5pqp/kX6whXFP5vwACOHMmubJuSMpQ9e/bUy1OubwkJCZfdX2eB91tvvcXHH3/M\nsmXLaNeuHXv37mX69OnExsby4IMPXvZYkassCIIgCDcHQ0E2emc5tuaBuN1ZpI7sjKz3LFhz9k+j\nLl/ZQ1GI+PonjJmFtJnxIYc/eMwToDcQtyxRWF5GXmkJ+WUllDnsOF0uXLJUpa3D6cO53E6cze7K\n2ewu2BzeN/XptDb8fCwEmi2E+JURGlBBZKCVlpFWOsQWERrQeFfIPHHCwrFjudx992AOrv6Unadz\n2LfvOA8G9kVRtEREFGI0liBJnqoynTv3ATxpxFdLrVbTtGlTUBQCtm4FoKRPn7p7Mo2cSrmWV60a\nERERPPvss0ybNq1y25w5c1i4cCEnTpwgLS2N+Ph4du7cSdeuF5cDHTZsGOHh4Xz88ceV20pKSiq/\nDwio+aYLu90uZryFXxVxTdevC7OX3bp1a+CRCI2duFa8KW4XaLQ4nRZcu39ALj6Ls2e7a+pLVWol\nYMTTaA+fwdWlFSVfvwi+N+b3XoXDTmZBPufycjibl0NWUQGyfLEGeJP0MsZ+cYg3h9zFEWMSJdZI\nSsojKS2PoKAktnKlSIBAv0K6t0mja4uj9G5TSK/O8TfNRGJZWQU//niYsWPvRqeLwGbTkpZ2+sZe\n70ePQmIihIZCbm6jWEW0Llwphq2zGW9FUaqUuFGr1ZXvhuLi4oiMjGT16tWVgbfdbmfLli3Mmzev\nroYhCIIgCEIdURQF19nDqN58nvLnZiBJJdDKDFxb0A2g+PtSuux5AoY+hW5PKgF3/52Sla9WmSlX\nWcoxvf4FFU9PuOJS6peSZRmLtZyCUguFpSUUlFgoKC2hsLSECofdq60kackuaMe53I5k5HZgzq5X\niM/9kQlL93Br+7eQVBfDJLVKol2LdAZ0zWXUrWW0jbGiUsHJk57FfFQqFaryCgyfrkYJMOMYP6Dq\n4FxuzE+9j+2Po5Dir72k3rVIT88nOjoSvT6SwMAgCgvPYDYnoVarMZkgOPhal2q/RgcPgk4HQ4b8\naoLu2qizwHvUqFG88sorxMXFkZiYyN69e3n99dd54IEHAM8FOWPGDF566SXatGlDQkICL774In5+\nfowfP/4KvQuCIAiCcCPIbhfS/L/inPIwTqUIybcInhwPkqXuzhEZQukXswgY9hSaMznVpqeo7A58\n/vUNSpAfthlja+zL6XaRWZBPen4u6fm5ZBTk4XC5qm2rKCoKS5qTntuR9NyOZOW3Q+90U6HxLMP+\n55j5DClZTc/yHXxqfJb9YyYRE1FBbBMrraLL8Tc5URWWooRVXztau/Uw5ucWIEWF4hh7W5U62aY3\n/4vx09Votx/BsuWdegw41SxatI7RowcREtIEjcaHl15awPz5rxMaGoNKpWLmzKfq6dy19LvfQUoK\nlJY27DhusDoLvF9//XX8/f2ZOnUqubm5NGnShClTpvDcc89Vtpk5cyY2m42pU6dSXFxMr169WL16\nNb6+vnU1DEEQBEEQrpIiSSiSG4dkweE4DT2icduOeQJHFWCofsZZu/UQPm/+F/uEQTiHX12ZQCk+\nGsvK1zD8b3u1+zWnPIUXTOf7V0IDcLpdlbPYWUUFpOfnklNUiHxJ1qzbrcdS3pxSqydFpNQaTok1\nklKr53tJuliPur31AJuP9GXDHQ9j/9NYerQtwrx9Cvzuee7ZMp/Bs6KREmMr2+vW7MF/0svY/jia\nimcuVu64wDWgK+5WzdCmpmP4ZjOOsf0vPp/9J/H55+cAWF99pE6DbpXKwKxZnzFhwjg6duyKRuND\ns2ZW/Py6EhAQBsCXX/7fFXq5AU6cgL//HQID4f33wc/P8+83pM5yvOuSyPEWfqvENV2/RN6uUFs3\n/bViscC778LDD8Mvqon9kqIouBa8gitEh+3Wq1uwxfevH2D6z/dUPD6Gir8/cD0j9iLJMjnFhURN\nepWwbcfZN6Q9X43pQIm1/JJxg9UWTGFJLAWW8/9KYrGURXnlYv9SVIiN5I4F9OtYwPgPp+G/Ziu2\nKcOxzpl88Xk9+S9Mi37A3b4FllXzKmeu/Uf/Df2Wg1hnTcI21bMmycmTJwCIj/dUszB8tga/GW/j\nbheHZf0bntl8u5PAgU+gPXYO2+ThWF+azPXR8tFH60hMTGLw4GHodEGcO5dOVFQUBkPtF7q54U6d\ngvh4iIqCzMyGHk29uGE53oIgCIIgNBJz5sC8ebBoEaSmVtmtKApKbiauIB/s9jRcQ1qD/ipDAllG\nv9IzW+288/oWxSm32cgo8KSJZOTnkVVUgFuSaJIcxVPbU0lafZi3ooZwTNseS3kTSsqbUFzWFLuj\namCjVsm0iContkkFMREVxERaz3+tICa8ArOPp2KJZv8p/NdsRTHpqXj8bq8+KmZNRLfruGfFRq2m\nsr1+y0Fkswn7fYOqnPcCx9234fvSYrSHT6PbuA/XbZ3xmfsZ2mPncLdsivXZ+6/hFVKzYsU+JMnA\n+PET0OmCGTeuLcHBwRiNnpUn4+LirqHfGywuDnx8ICsLiooguHGWk6xPIvAWBEEQhF+bNWs8X0+c\ngJ07oXt3FMkN55fZduacQDX/b5T95T5Qq8Bw9Qu3aPeeQJNdiBQVirtT/FUdK8syZ/NyOHL2NGk5\nmRSXl1Xus9oCOZPdnzNZ3ckpbE1E2BM8mPcx/b/LZmHCXK9+zCY7SXFltG9RTmJsKe1iS2ndvAyT\nQf7lKavweW0pALZJQ1Eigrz2KWYfLD++7pUOYnrvGwAc9w1C8b9MiqxBh23ynfjO+RT96l24buuM\n3CwC2Wyi/J3pngWFrkjFzp1nOHQom6lTH0WnC2LAgI7odDpMJs8nGLGxsbXop5FRq6FdO881efAg\n9OvX0CO64UTgLQiCIAi/JrIMp0+jqFQozz8PHTrgdlhQP3YP1uefxW1youhsMPO+y9fcvgL999sA\ncA7rXat+ZFnmXH4uR86e5mj6aax2T4URRYHCkhhOZ3XnTFZ3cotaeR33You/o/bT8v1tf+DRVido\nEWUjromVFlFWmobar+kpaPekYli1A8XHiG3amOobXRJ0qzPyMXyzGUWjxjZlxBX7t09MwdWjLe7e\nnuov9geH4rirL8plapbn55fw449Huf/++9DrI+jatQfx8SX4+Hiqn0RH/0ruh7uQ0719uwi8BUEQ\nBEG4uSnHj6MqLUW+ox0Vd3XCWb4FRbHBrClgsMCFO7t+GbE6XN4z37J82RsAdZsPeA4b1rvmsSgK\nGQV5HD6bxpFzZyi3VZzvWk12YSKnMnpzOrMHZRXhlcfodW76JOUxtFcB/TvnEx1mQ60ezp3UvMr1\n1VLMJpyDuuNu0xwltOZ7yS7lGNsfFAU5OuzK/Qeacd+SVGWb12NF4ejRdJKS2qLTRRIW5oPZLOHn\n1waA8HDP6t6/OrfdBj/+CCtXwlMNXFmlAYjAu5FZuHBh5UqfmzZt4tZbb63SJj4+nrS0NPr168f6\n9etv9BCF87Zu3cqaNWuYMWPGZW8CFgRBuBHkw7txu8pxtAzDffAzVGfTcTdRg+IJdmuqTAKgW7sb\n81/+hfXFh3Gm9MT0/nKMi1ZhWTWvxlnakpWvott6CHevtlX2FZRYOHjmFIfOnKpMI5EkLRn5nc4H\n2z2xXZKfHRZoZ3D3PAZ1zyW5YwG+xqorR14N9blcfF5bhu2REUjtquY+S62aUfrZ3z1vLmpBjg6j\n/O3pnun5OqFFo4ng/fc/4733JuDnF4CfHzzwQGwd9d+IPfUUNG0Kd93V0CNpECLwbqRMJhNLliyp\nEnhv376dtLQ0jEbjTbNC1q/V1q1bmT17NpMmTRKBtyAIDUKRJBSVCocjF6n0CG5XEW5HC4j0g8jE\n2nXicuP7/AI0GfloTmeDWo1u7W40aVmY3vvGs4BNdXRaXP06VT4ss1Vw+EwaB8+cIruowDM+BTLz\n2nP0zO2czuqO03UxXSI20sqdt2QzrFcOnRMsdVrS2vTRCozL1mFctg5nv47Y/jAS1x1dqs7gX+1J\nr/Pv7iefbKJTp54kJ9+BVmtiyZLPr6u/m5JeD+cnGH+LRODdSKWkpPDll1/y1ltvodVe/DEtWbKE\nNm3aoNHUXCrpZmC1Wn819dsbYUVOQRB+hRSXEyX9NKq4VrjdNqRT+9F8NJ+yv01FlssgIRCofmGX\nyzEu+gFtajpSbCS2ycMBqPjr79FvPoDxg++wTRmBEuJf/ZgUhfT8XLYePciJzPTK34dut57jZ/ux\n/+QwikpiKtu3aV7KnbfkMKx3NokxZdcbx9bI9uBQkGSMn61Bv3E/+o37cSdEU/76Y7h71vINSR1w\nuyWKisqJikpAr4/mjjuaExUVhU7nc8PGIDQuv501Om8y48aNo6ioiFWrVlVukySJL774gt///vdV\n2iuKwttvv0379u0xmUxERETw8MMPU1hY6NXu22+/Zfjw4TRr1gyj0UhsbCwzZ87E4XB4tcvNzeXh\nhx+ubBcZGcnQoUM5cuRIZRu1Ws3s2bOrjCU2NpZJkyZVPl64cCFqtZr169fz+OOPExERgd8lBfN3\n7tzJ0KFDCQwMxMfHh+TkZDZs2ODV56xZs1Cr1Rw7dowJEyYQGBhIWFgYf/vb3wBIT09n5MiRBAQE\nEBkZybx586qMy+FwMHv2bBISEjAajURHR/PEE09gs9m82qnVah599FG++eYbkpKSMBqNJCUlef0s\nZs2axcyZMwFPCSe1Wo1arWbTpk0A7Nmzh6FDhxIeHo7JZCI2Npb7778fu917uWJBEISayPYK5O+X\n4nKVY7NlY83bh2vRHIqL11NSsp7y4DxK/jLeE3RfidOFccFKdD/u8dqsKi7D51VPdQ/r7Acrc7zd\nPdrivL0LaqsN0ztfVR2bLHP03BkWrP6OhWu+JzXjHIqiUF4RwtYDE1i44t+s3/1HikpiCA+08+S9\nqWz713o2vb2JmeNSaRd77UG3dsdR/O6fA+U2NIdPo/3pUNXxxTXB+tJkig4swPr8RKSoUDQnM5HD\ng6rp8dq5XG4KC8tQqTSAhhMnClm8eBs6XVMMhhbs3l3O4sUHCAjogY9PUzp27EhY2JVzxIVfLzHj\n3UhFR0eTnJzMkiVLGDZsGABr164lLy+PcePGsXTpUq/2jz76KAsWLGDixIk8/vjjnDt3jrfffpsd\nO3awc+fOyoL6CxcuxGQyMX36dAICAti2bRuvv/466enpXn3efffdHDp0iGnTphEXF0deXh6bNm3i\nxIkTJCZenC2oLt1FpVJVu33atGkEBwfz97//vbLA/MaNGxk8eDBdunTh+eefR6vV8umnnzJo0CDW\nrFlDv1/c8Txu3Djatm3L3Llz+f7773n55ZcJCAjg3//+NwMGDODVV19l8eLFzJw5k65du9K/v2fV\nMEVRGD16NJs2bWLKlCkkJiZy5MgR/vWvf3H48GGvoBpg27ZtfPfdd/zxj3/EbDbz1ltvMWbMGM6d\nO0dwcDBjxozhxIkTLF26lDfeeIPQ0FAA2rZtS35+PgMHDiQ8PJynnnqKoKAgzp07x3fffUdFRYVY\nIEcQ6tr69Z5KCTfhYjdK2nGUmHhUajWyy4H8yp9w/ekpJGxIrhJ0GTuxWcznV48E+7SxIFs9B6vV\ntU6VMHz+I+an3sfdNgbLbZ0qj/N5bRnq4jKct7bHmdLT65iKp3+P/sc9mP7zPbZHRqFEBOFyuzlw\n+iTbjh6kqMyz1LeiQE5hGw6cGMqpjFuQzy9e0znBwuQ7TzOiTxZ6XR19Mqgo+D6/AN2u46gmvYxu\n22EUfx+Kf3wTJbJqTWglwIztsbuw/WEE2p3HkOOaVNOlUvk3y+Vyk5qaSbt2MVXaAZSX2zh6NJ2e\nPTuj0QRy9Og5Pv/8GH/603SCgroSF1eIorQiIKALAMOGteP8n3BBAETg3WipVCrGjx9fOSNrMpn4\n7LPP6NWrFy1atPBqu3XrVj788EM+/fRTr9nwIUOGkJyczCeffMLkyZ5Vsj777DNMJlNlm8mTJ5OQ\nkMCzzz7La6+9RnR0NBaLhZ9++ol58+bxxBNPVLZ96jrvPvbz82PDhg2oz//CVxSFP/zhD/Tt25fV\nq1dXtnvkkUfo3LkzzzzzDD/99JNXH926deOjjz6qHHtsbCx//etfmTNnDk8//TQA9957L1FRUSxY\nsKAy8F66dCmrVq1iw4YNJCcne/U3YcIE1qxZw8CBAyu3Hzt2jCNHjlS+1v3796djx44sXbqUqVOn\n0r59ezp37szSpUsZNWoUzZs3rzx2+fLlFBcXs2bNGrp06VK5fdasWdf1+gmCUI2iIrj9ds/3Tifo\nrr4e9dWQt6xC6X4baM8HvJ9/CCMmgMkHUMGX/0a5cxwYDSiKBO/OQZk0DeX8Y80zU3E9/TyySY+i\nuDC8PxvbE1ORDSokyYq2f2vc1n2Vi7a4x/SvcSzVUeVbPDdD6rz/vDt+dzs+//wc7dGz6JdvwTm6\nLwDOO7qg27wf64sPV8lfdnZoSfnArkgFFvbt3sOZIB0ns9JpuSMNZ1wILpM/x8/25dCpFAosnhsY\nNWqZUX0ymTz8NN1aW+o+lUSlwvrcRAJHPI1+wz7Pcxvc8/J1tQF0Wty3JOFwuNi48RCDB/dErfaj\nrEzi3nufZO3a/6JWm7BYyvj00+W8/fYYFMVFVlYWkyc/yfLlH6BWaykrs7Flyx6GDPH8XJKTu2Ay\nRSNJoNHoiIyMJDIyso6ftPBr8tsJvGv6v7+m/NyrbV8Pxo4dy7Rp0/jmm28YNWoU33zzDS+//HKV\ndl988QVms5lBgwZRUFBQub1169aEh4ezfv36ysD7QtAtyzJlZWW4XC769OmDoijs3buX6OhoTCYT\ner2e9evXM2nSJIKC6uajucmTJ1cG3QD79+8nNTWVp556ymvcAAMGDOCdd96psoT6ww8/XPm9Wq2m\na9euZGZm8tBDD1VuDwgIoHXr1pw+fdrrNWrVqhWJiYle5+rbty8qlYr169d7Bd79+/f3eoPTvn17\n/P39vfqsSWCgJ8fyu+++o0OHDl45+oIg1LHDhy9+b7VC4NXnOF9KcTpAozn/yZ2C642ncY8djxRg\nRJIq0O37H/ZmVmRfz+8lgzYHR8l6sOsBFQYyzj/2vAHQJQXhqtgGLs9j9cPDkOUTYPcE1o6ZE4AS\nOF/Ew9225XWN32/6W+i2HKD0k7/huq3zxR0GHRV/vhe/J97BZ+5SnMP7gFaD646uWPp3BrUap9vF\n0XNnOJWVQX6phcLSEtQDInHqm0LxGSiGAIuNhxfsxKo1EtU7gzLJM8scGuBgwsBzTEw5S1Ro/abU\nuXu3wz5hEPqV27C+NAXHmJprQSuKwmefbWTChLswGEIwGMzs3buF3/3uDjQaDYGB8MMPa/Dx8eSw\nR0aG8u9/L6w8vmXLZqxYsaryniSzGebM6VTdqQShVkRE0IgFBQUxePBgFi9ejFqtxmazcc8991Rp\nl5qaSnl5OREREdX2k5+fX/n9oUOHmDlzJhs3bqyS23wh/cNgMDB37lyefPJJIiIi6NmzJ0OHDuW+\n++4jOjr6mp9Py5bef1BSzy9jfGnQfCmVSkVhYSFNmzat3HbpzDJ4gmydTlel1qm/v7/X805NTeX4\n8ePV5tapVCqvttWdBzw/j+Li4mrHeql+/fpx9913M3v2bObPn0+/fv0YMWIE48ePx8dH3FAjCHXq\n6FHP1/vvv6agW96/HaVZSySz0TPj/PwM7BPGoyjFSFI55T2jkZTjYPek67nuveP8kZ5I2XF7D+/H\nAy6ka3jK1Lk6tfE+X3g9LpGtKGh3H0dlcyK1iKqy23Hv7fi89V+0pzIxfLkBx7g7UBSFzKIC9p1K\n5dCZNJxul/dBBi12h5l8Sxz5xS3ov+04sJpV/imUScF0b1PEpJSzDO+TjUFXu9J8daF8/lSYP7Xa\nSbJPP93AqFEphIXFoNUGolKlodMl4efnqT712mvzvdr7+1d/4yh4/j78WgoBCI3DbyfwvtqZ6kZS\nqWL8+PHcf//9lJaWMnDgwMpc4kvJskxISAiff159WaILM9YlJSX0798fPz8/XnrpJeLj4zGZTGRk\nZDBx4kTkS+qZTp8+nZEjR7J8+XLWrFnDP/7xD1566SVWrFhRJe/6l9xud7XbL01xuTBugLlz59K1\na9dqj/nl862umktNZRUvrTYiyzLt2rXjzTffrLZtVJT3H6maqsbUtoLJF198wc6dO1mxYgVr1qxh\nypQpvPzyy2zfvl3cWCMIdelC4N22ai3p6sjbf0QKC8PVJARJsqA58AMOJQ4pOgRQ4C/jQaXCfdJz\nY7qUkFBPA6976tPZqIvKkMMCkZtVs/CKTkvFzHH4/fF11P/bztaO4ew7lUpBqaWyicutJzOvPXnF\nLSmwxFFgaUmp9eLvrNlnBgCQ168v6/66ifYtSuv9eVWr8ve+moUL19Ov3620a9cJjcafpk0r8PXt\ngNnseQ2mT5/RMGMUhGr8dgLvm9TIkSMxGAxs3bqVRYsWVdumZcuWrF27lp49e172nfn69espLCzk\n66+/9spzXrNmTbXtY2NjmT59OtOnTyczM5NOnToxZ86cysA7KCgIi8XidYzT6SQ7O7tWz+3CDLjZ\nbOb2Czma9SQ+Pp7du3fX6XmuVEe9e/fudO/endmzZ/PDDz8wdOhQPvroI5555pk6G4Mg/OZVE3hf\nuFlOURTknRuR1BLuxATcbguUHEDS+yIFnM/DTbnwpv/8m+qbeH0E3a7jALi6ta7yPNySm4yCfM7E\n+8HMFLZE6ZH37gBAUVRk5idy/OxtnMrog9PlPUli0kskxpbSO/I0/bdvRNFqGPtaBEpAQwTdKr78\ncjuxsQn06zcAnS6AW24JolmzFvj4eH6mv/td1U+GBaGxEIF3I2cymXjvvfdIS0tj1KhR1ba59957\nee+993jhhReYO3eu1z5JkigrKyMwMLByFvfSmW1Zlpk/3/tjtwspKJfOUDdt2pSwsLDKdBTwBM4b\nN270OvbDDz/06v9yunXrRnx8PPPnz+e+++7DbPZeHS0/P79Ws8O1WUjonnvuYeXKlbz33ns8+uij\nXvscDgcul6vK+a/kwpucoqIir9QUi8VCQECA17g6d/bkWl76+glCnfrkE/jwQ/j6a89a078Ryktz\nUEYNR+7WBclhgW8Xo1gtOEYPR5atqKTjKCoFqeL86o1dry+HujHTng+83V1bI8sy6QV5nMnN5kxu\nNhn5eUjy+UTyqP9n777DoyqzB45/79T0QkISEgIhIaF3FKSJIEXsDXvDLro/K65lF9eGXVcWdWXX\ntnZxXXfXhquIIKD0TkhIb5Myvd6Ze9/fHxMCMRUI1ft5nnkCM7e8V0Ny7nvPe064g6XNmU5+6WTy\nSyfj8u77WTsi186YgQ0MzXYyJNtBToYHg15g/ngZEoLQiDxE/IH9vDwUq1btxus1cPbZ52A0dmPi\nxCy6detGVFT4SeW4ceOO2Fg0mkOlBd7HgSuvbL1r2N60h4kTJzJ37lyeffZZtmzZwvTp0zGbzRQW\nFvLpp5/y2GOPcfXVVzNhwgSSkpK45ppruOOOOzAYDCxZsgSPx9PsuPn5+UyZMoXZs2czcOBAzGYz\nX375Jbt27eL5559v2u6GG27glltu4aKLLuL0009n8+bNLF26lOTk5E6lZEiSxN///ndmzpzJwIED\nmR8dR98AACAASURBVDNnDhkZGVRVVTUF9N9//32Hx2nrXPu/f+WVV7JkyRLmzp3L8uXLmxaU5ufn\n88knn7BkyRImTZp0QOc56aSTAHjggQe47LLLMJlMTJ06lffee49FixZxwQUXkJ2djc/n480338Rg\nMHDRRRd1eD0azUG55prw1yeegDZSqk4kQggURSZk3YFqqsJr3gouBcb1AGMmBMvCG+a2zHU+YUmg\nxETyU6zCss8+wPOrvgEhxUhNQx6VtYMpqxmJxZrX9FlmipeLJ1dw8eRKcjI8vz4yAGpCDEQY8c+Z\n1eVDD4UUDIZwLez168vYurWcW265FoMhnr59+xAMKkRHhyunDB6sdQrWHL+0wPsY1JkZ3F/Xyl64\ncCEjR47ktdde4+GHH8ZgMNC7d28uueSSpvSKxMREvvjiC+655x7mz59PbGwsF154IbfccgtDhw5t\nOlavXr248sor+e6773j//feRJIl+/fo11Qnf68Ybb6S4uJi///3vfP3110yaNIlvv/2WqVOntriG\ntq5p4sSJrFmzhscee4xXXnkFp9NJjx49OOmkk5pVMGmrNnhn35ckiX/+85+89NJLvP3223z++edE\nRkaSk5PTVB6wI78+z6hRo1iwYAGvvPIKc+bMQQjBsmXLmDx5MuvWrePjjz+mpqaGuLg4Ro4cyaJF\ni5qCdY3msOlE5Z1jlRACVVWRfvgCJkxHZw5XDhFKCEkf/nWl1Fah/n0B3tuuJxSqQQxPgOFjwecl\n6s9LMGwtwvnuw0fzMrqcrroBNSUB2lh7YnU52Vqyh63jk7ANnga4EH6JUMhETUM/KusGUVk3GEtD\nHoq6r9xiTGSQc8ZXM/u0CsYOtHZYEjw4ZST2fy9AGdp1Tw10umjWrq3kgw++4rXXXsZgiGX4cDtZ\nWXZiYsL59QMH/nae4GhOfJI4Bvtd7/84Pj6+7TvbX5ea02iOd9r39OG1bt06IJzmdMLZuBFGjoTs\nbNiz52iP5oCohTsIRRnwRTsJBuuJ+uAL5MsuQTJGo5OMRN50E/Kr76AaFWR/BdKefJTsns3zmFWV\nboOvRVdnx7byLyj9WlYmOhCFhQUA9O17dBdXmv/xDbF3L8Jz23lU33cxDo8bh9eDw+PG6fFQY2ug\nsqF5VSaLtS9bCmZRWD4eRTU1vS9JgkFZTsYNbmDc4AYmj6gjynzkKpHs1dAg88ADb/Dee29iNieH\nb7Ykqc1F7ce6E/rniuaAdRTDajPeGo1GcyIYMgSio6GoCCwWaKO86LEmFAoQ3PQ1gRQDoYHhVALv\nZdMBGwTD5TsDL98Pys6mWtfkZLY8kE6HPHk4EZ/8gHHZxkMOvI8mIQQNTgf1P21g7P2vAWB6/d+8\nlerDGd/6jbmiGCisGMfWwjOpaQinkEiSYEi2g3GDGhg3pIGxA60kxgZb3f9wUlWVRYu+Yu7cG4mL\nyyQxsRvPPz+SyMjwTLauk903NZoTgRZ4azQazYnAYIC334asLEhKOtqjaZcIBlE3r8E/IBO/vxAx\nuYNZ5U4GZsHJI4j45AdMyzbiv+XcLhjpkSGEoMHloNRSQ4mlmlJLNbLDxb3PL8cQVAiY9JhlhSHb\nqvlpfJ9m+7p9iewomsHO4jNwecP1qOOjZa6YVs61Z5SSleY9GpcEhK/LYEggKiqLIUNkIiJyiYgI\nj7FPnz4d7K3RnJi0wFuj0WhOFBdeeLRH0Kb9sxoDrgrE9+/iyzwjnLcsBEgS+q1FRHy8jFD/XgSu\nmNbO0cLizn8IDHpcC+9EpHVDPjXcUdC4ehv4ZYgwdXCEI0sIgdvvo9Zuo9Zupc5uo9Zho85uJ6g0\n73+gM+jYOrgH8jYrj06+haAvmi1RY5B/ScQvx+P1x+HyxOANmJv2GdDbyQ1nlXDhqRVHJYVkf++9\ntxqzOZmbbroTnU7HRRcdv08gNJqupAXeGo1Gozks9gbboZAfHp1L4IJzkXsnoKoumHNW03Yxd/0F\n5CDK4D5EvvY58mkjOg685SDG1dtBgEgIl7YTqYmEBvXBsL0Y47pdBCcMbf8Yh3BdwVAIORQiGAqG\nvyr7/u71+3H7fXgaX26/H4/Ph9vvxS/L7RwX7O50KixDqKofzt8YgtpHh7+isbRrQ8t9jAaV2YO2\nc2fSP8k6PxPR6+gsRFRVldLSOvLyhmE29+aqq8YTHR2tpZFoNL+iBd4ajUZzPLPZ4NNPYdQoaKwX\nfzSFy/yFUF97hFD/vgRG5aEodqTrZyCiJVCbN13RVdVj/uh7UAXOy8NdEQ1b9jTNgrdFX1SNpKgo\nWWnNZrbdC25CJMai9GslD/wgrsXhcTfNUNfabdQ5bNQ5HPtqYh/0scHjS6LB0QurMxObsy9VdYOx\nuxObbdczzcuEIeVk9fCQFCeTFC+THCeTHB8gKV4mPjpIxOcriLvpWeTVo3B+OP+QxnVwJOrrzSxc\nuIy33pqDTqdDWyOu0bSuSwPv6upqfv/73/PVV1/hcrnIzs7m1VdfbVYf+ZFHHmHx4sXYbDbGjBnD\nokWLGDhwYFcOQ6PRaH471q2DG2+E8eNh5cqjMgQhBMpP36LYK/BPHEkoVAtnDUVERYBiDW8TE9Xq\nvhF/+y9SSCFw7gSC44egdotF1+BEV1mP2rPtBlr63eUAKHnNA+zQKYMO+jo8fh8VdbWU1VkorCjF\n7vUQ/KVlgC062eRSCAlfIA6XNxm3NxmXtzsOdy+c7j7U2jLwBSJb7JMUF2DCkAYmDqtn4tB6stK8\nHZ7HuL6xY+VJ/Tp1nV3lyy/XM3HiaaSnDyMpKZ533pl5RM+v0RyPuizwttvtjB8/nkmTJvHll1/S\nvXt3ioqKSNmvg9rTTz/NCy+8wNtvv01eXh6PPvoo06ZNIz8//4C7Bmo0Go0G2Lw5/HXYsObvBwJg\nNrfcvouodTWo+VsIjhqMLNehdqtHJEiowXBATBuB9v4kl5eIt74GwDf3fJAkQkP7YvphI4ZNBcjt\nBd4F4fOE8noe3PhVlXqng4r6WsrrLJTXWbC6ms/GB+Qo7K4MbHtfzgzsrgzs7h4gJExGPyajnwjT\n3pdMhFnGH4jF7k7E6owjpLT9a7ZbrEz/3i7693LRr5eLMVnV9O8nd3YtaRPD+n0dK48MPSZTTxSl\nCiEyMZsTjtB5NZrjX5cF3s888wwZGRm89dZbTe/17t276c9CCF566SUeeOABzj//fADefvttUlJS\neP/997npppu6aigajUZz5MgymI7iIr5fB97LlsG118KYMfDxxwd9WKGq4HFBTLgKhVJRhPjPOwSv\nnYOieFCtRejLf8bfPxDeISUWiD2gc5jfXYrO5SV4yiBCI8KVTULDcjD9sBHD5j3IZ7XdCly/uyI8\nrtyOU0r8cgCLzYrFbg1/tVmpddgIKftmsxXFQJ09j+r6/tQ09MPSkIfbl9zBcWPwyzE4W2/0CEBi\nrEx6so+MZD/pyT76Znjo38tJ/94uusfLTbPZpi/XEH3p67heu6fdWXtdZR2Rf/4UJbtHuHJLIIhh\nyx6EJBEamdfmfp2xfXsF+fn1XHbZhYDE118vZ8uWXdx3X7iZ2Xff/UR+fil33/0ARmMM1113eHLo\nNZoTWZcF3v/6178444wzuOSSS/jhhx9IT0/nhhtuYO7cuQAUFxdjsViYPn160z4RERFMmjSJVatW\naYG3RqM5/nz4IVx5ZbhN+/33H50x/DrwTkuDsjJQO65qIbxuiIxGVRVUaw28t4jQzXcihAy1FRhf\nfAbvo39AUVyokhXD8GRC3i3hnZMgOO3kQxq6zuVFRJjCs92NAudOQOmbQXBs+ykj7hdux3fb+agZ\nzYPjkKKEm8rU11JRX0dlQx12t6vl/t5u1FjzqKnvR01Df2ptOaj7dXUE6KlWkN3DQ2LfWHIy3OT2\nDL+y0z0YDQK3T4/bZ2h6eRq/xkcHyUj20SPZT3REx7ng+oIKYn73Z3QOD4bNhe0G3vrdFUS++SVq\nUhz+q2Zg2FWKJIcI9ctExEV3eK79lZZa+eqrzdx5560YDAn07u1Gr68kNjac/nnmmb2ZNk0mLi48\noz1lSiYnneTCZDqwGyyNRrNPl3WujIiIQJIk7r77bmbPns3GjRu54447eOqpp5g7dy6rVq1iwoQJ\nlJWV0bPnvkeDc+bMoaqqiq+//rrpvf27/hQUFLR5zt69e9O9e9uPIjWa401dXR2lpaVHexiaThpw\n1VVE79qFY+xYChYuPOLnl2SZEZMmIakqG5cvR42MBFVl+LRpGJxONv/nPwTT0pq21/u9KBFR6HQC\nvdxA3itPsfP/riekOCHoJ6asBmefHkf0Ggw2N6H4qE7X6v41VVWptFupdTmodzuxedyoQoAQpFrc\n1KbEIKtm6mw54Zlsay41DXl4GmezjapMYshGrSmVrNR6hvSuZnBWNYN7VzPtpRdI2FBI6e1nU3f2\nyR0ndR+EmM3F9Jv3BkanF9v4geQ/O6f98wjB4BteJmZHGSW/Owf7Kf1JW/ITcnIcVdee3u65XK4Q\nH374C3fcMQeIxelUKCjYw7BfpylpNJqDlpu7ry/BYe1cqaoqJ598Mk888QQAw4YNo6CggEWLFjXN\nerdFOgw/zDQajeZwMtbWEpUfzq3d89RTR2UMOlnGcsUVGOz2cNANoNPhGTyY+FWriN24AeuMGegM\neiQ89Pv7nyi44QoCRh8yMptvPRuCje3G9dIRD7oBQokHt75HDgUprK1ht6USb2OJPlXV4fSkYXVm\nsuCdNxliKeS08f/hR2YiRPN25On6GhbV3MW0sv9SMXQYJY9dTlxUoOlzyR9ELwcxuHzkLPiY5KUb\nKLr/YgKZ7aefHIjEH7aSO/9ddHII64SBFD52VcfBvSRRce3p9J/3Bunv/4DlgnGU3HtBq5sKoWfh\nwm+5664bMRgS0OsjOPnkeILB7kiSRFQUWtCt0RxhXRZ4p6ent6hO0r9/f8rKygBIa5x1sVgszWa8\nLRZL02etGT16dJuf+f3+QxmyRnPMiY2Nbfd7XnNo1q1bB7T/c6XTXn45XN7i/PMZeeqph368gzV5\nMgDdCedlCyUEM2fCqlX02rqU7ucNQU6PQggZ7zN3knGcT3RYXU5+3rWdTUW7CYZCON0pbNx9DjX1\n/bG6eqIo4QWlF4hChlDIyRXbWZE5i0FZDkb1szEqz84Ew0YGz/sDhqIqAPpgpduQzGZBb2FhAfnP\n38DArdXEPLiY+PWFDL/6ebwPXIHvtvNbHduBMm23IAUVfNfMRH3qZrIN+o53AsjpS+jtZZi2FzPo\n5xL8c2Y1ffT6699y2WUXkZqag8mUxA03DGHYsAmYGtchnHTSmC4Zu2afLv25ojnu7Z+10ZouC7zH\njx/Prl27mr23e/dusrKygHB72LS0NJYuXcqoUaOAcOC8cuVKnnvuua4ahkaj0RwZn3wS/nrxxUd3\nHOytnR1AfelB5BGDUAZHEidJhEI6Aj0MIBqbthzjQXcgKOP0elFUBUVVURS16c9yKMjW4j3kV4RT\nsRzuNNbtvJD8ksmoYt+vsh5JXvr3ciOyhsH7MC/5E277YEhTrrXpyzXEXP8iOo+P0KAsXK/di9K/\nja6KkoR8wanYTh1B9B/+RsQnPyDVt/9L9UDI507AkZkSXljazv8bIUTTk+GSEguxsZGY7p5N7E3P\nsnJ9IanT/OTmDsNg6MbQoRKxsQOJigrPzE+ZMqXLxqvRaA5dlwXed911F+PGjePJJ59syvFeuHAh\nCxYsAMLpJHfeeSdPPvkk/fv3Jzc3l8cff5zY2Fguv/zyrhrGCcPj8fDMM8+wdu1a1q5dS0NDAwsW\nLOD+o7WAS6PRNDdnDsTHw1lntfysuBgUBfr2PWynF0Kg/PgFSkMZ/qljCAZr4YqJYDSAnIi14D1E\n/DFYpjUYCo+xkaqqFNVUsXFPPvkVZagdLAq1udLZtO08dlRMQQg9OknlolPLuHpmOQN7u4iLDrde\nl6zJiA8kum3djFC9QHgm3Lh6OzqPj8B5E3G9dAdEd9zpRSTF4X7lbgKXTCF48oCDv/ZWtFeJRJKM\nLFr0DX379uP8888B9Cxf/gbDhw+lz+X3Exx3OfZN28g0DSQmJnycs1r7ftRoNMeMLltcCfDll1/y\n4IMPkp+fT+/evbn99tu5/fbbm23zpz/9ib/+9a/YbDbGjh3bagOd/afpW0tM38vv9xNxgrbHKikp\nITs7m8zMTPr378+3337LU089xbx584720DSH0Yn8PX0sOCKPhHfvDqd/mM3w00+Qnt5lh1Ytlajb\n1yGPHUUwaCFUsxsIIZKO3TrKuvJadBV16Gpt6OrsmP67Ggw6ah64jF+Mfjbt2Y3TG67HZ3f1wO7u\nQYQc4P5/fUyK087tN9yGJIXzlfNLJ1FQPoEVWybRO1DKX65ZzPl3mMnu4W313PHT78G4sQDHh/MJ\nTg0/aSWkYP58JYELJrU7y1xYGF7Y37dvbpvbHA6rVxdQWxvkssuuwmRKxmZzkJCQgNFo7HhnzVGh\npZpo9tdRDNulnStnzZrFrFmz2t1m/vz5zJ9/NFraHl/S09OpqqoiLS2N0tJS+vTpc7SHpNFoOiM9\nHTIz4ZdfYMYMWLECEg4+MBZOOyImlkCgnqBnF1L9egLexgogScfgjPavxNyzCNOyjc3eUyWJtyZ2\npz45GiEkSqpHs6XgTMotw8MbCMGi8tfpHqpny7fXUhqR1bSvQacwTN5GtOzi3v+rR6Q2b7G+v6Zg\nW91vfsmgJ3Dhkc/JN33zC2pyfIsmN+Xl9axfX8rs2RdhNCYzaNBAcnICREWFb9i0yl0azYmlSwNv\nTdcxmUxNi0678KGERqM53GJi4IsvYOJE2LYNXnkFHnzwgA8jhCDoqUO692pcj96FapAhDphyEhDu\nVmj+10rkKSMJnjaiiy+iawghsGenordk0hClp9Ys4Yw1U5Tdjar4ZHbunsqWwlk43OFqKpGmEGMG\nWtHpoLhyCN3LlnF192Uszz4LIaBfpps7xq8lerILNT4akdL+DY133mVw/zGQyigHib73FfQ1Voo/\nfITPqz1cd90VGI2JJCfLJCUlEBc3GIDs7KM8Vo1Gc1hpgbdGo9F0teRkePxxuOgi+N//2g28hdOO\nMJmQzJEoioxY+EeCM2Ygp0cSCjXA47eCTm6xn3HFFiJf+xwU9ZgKvF1eD8U1VRRWV1JUXYl3iBmG\njATCRWAaHL3ZUXQ6O/8zlWAoXAKxZ3cv159ZwhXTykmICQIQZewOL8L9/b/ijof3zfoaV4RTQJS8\nzE6V3jtcJKcH88fLEBEmAlfuawynqioFBVX065cJGJBlA/dctoAPaqyoA/uRNPNmYpd8QVzc0KaS\nfpmZh28tgEajObZogfcJwOFwMG/ePD7++OMWZWxefvnlFnn2mt+wpUvDM7Lj2m7FrelAKASGTvzo\nbCzzx6pVCK8XKSoKIQRi1f9Q0zNR0lNRFA/6VxcgTx5PMLcHqupDN7UfarwNQo19yNtoLGPYXhIe\nzqCstsfg9mHcVECof29EctvrZQ6FXw5QUlNF/DMfsa5PLBvTmq9R8PrjKbcMo6xmOOWWYXj93Zo+\nGze4gRvPKmbmyTXof1VJLzQsBwDD5sJm7+vzy4HGwPsoUFWVJUt+4vKUbsQ88DpKegqXfbuTf7z7\nAkZjJELoeeaZxXzwwT8wm2ORBNxT04AAdL9/iLj4NK6//vqjMnaNRnP0/WYC75deegmAO++887D8\n/WhxOp1MmjSJqVOnsnz5cjZs2MDNN9/M6tWr6d69e7Oa6ZrfuGAQbr0Viorgm29g+vSO99E0V14e\nbs1+6aXhFJJWCI8LEQoiEuJR/7EYUboZZfl7BMePQVW96JxbCEWXo8Y2zuJePzP8VXWHv6R3LqdX\nv70YgNDgttd/xN7+IuYv1uD6y50ELjn0snKqqtLgclBjtVJja6C0toZqaz0zv9zJhG/yyYk0su2h\nMyjyDKW8ZjhlluHU25vnTqQk+pk+upbrZpUwJNvZ5rlCw8KzwIYdpeGp8sbZa12dHSFJKHlH4meb\nhCSZ+OMf/8EDD9xOt25p6HQxWK1biLj2d4jcN9AXFPD0g+cTFzcIfePdw+ef/3ffIf79b0bv3h3O\n+7/00iMwZo1Gcyzr0qomXUWratLc3gonrVU1ufXWWwkEArzxxhtN740ePZo//OEPnHvuuUd6qJpD\ndFi/p996C667DvLyYMcOWkwx/gYccvWBF1+Eu++GCy+EJUsAEC4HwutGTUohFLLDJ39DiY/AP34g\nQsjgD4DJeNAt0VvlC5CUdQlI0FDyMUSYWt0s8i//JPpPb+G7Ziae5247oFMEQyEsdisWm5UaawPV\ntgZq7VZCitJsuynfF3Lev7ejSDpuH/sEbxh/15RCAhBhUjhlUAOnDq9n8vA6BvR2dS4DRAgMmwsJ\n9e/d8vq8ASRVQcREHdA1dYYkRXLXXa9x4YVn0bt3HkOGnMyqVasZPXo0kZGRzTfe+/0wfXr4ZrY1\n06fDt9+Gtz3KEzWaw0OraqLZ3xGtaqI5shoaGnjjjTfYuXNns/eDwSDKr345an7jQqFwzjHAww//\nJoPuLrG3ac7s2QCoaojAqn+iNlTgmz4aCMEZ4UVyTU1rIsxdPgzDrjIkVSXUL7PNoBsgOLo/AMa1\nO9vcBsAb8FNjbaDGFp7JrrE20OBytLmw2x+IwWLN5aQVNZy36nMArst5k3elqyAEA3o7mTKyjtNG\n1HLyABsRpvZrc7dKkggNb6OUX5SZQ50x2teURsfrr3/HqFEnceqp0zCZEnnyyf5UVlYSCunR6w1M\nnDix9YNcey089FA4haugAHJbGe9HH8Hrr8MNNxziiDUazYlAC7yPYz/99BPp6elk77cM3uVysWvX\nLu3OW9Pc++/Dnj3hhi6XXRZ+z2aDxLZLsWl+pawMVq+GyEjUYQOQfXX4/LtQRiUDyUDoiA1F6ZWK\na9Fd0MHMcWhYDsJoQL+zjLrySqySgsPjxu5x4fC4cbjd2D1uvAF/2+dS9TQ4emNpyKWmoR91tn5Y\nnen0DJTz8YZwoHlX3stYzzqd50dtYcqIWjK6t328o2Hz5mIMBhNDhw5Cp4tm4cIPSU5O4ZprrsRg\niGf27H6kpKQQEREHQM+ePampqen4wImJ4fSRN98MPwF54IHWt9Ean2k0mkZa4H0ck2WZ9F8153j3\n3XeZMWMGvXq10QJZ89u0Nx/5oYfCX3v1gqoqcLvhN5Cu1SUaU0vUM2YSeu95vJefjnqYFix2RCTF\nEZh9WuufCUGDy0FVfR2VDXVM6ZVIxp46vl/8Dvn9Uto+ppBwebpjdWY2vnri8mRRa8skGGo+qx5h\nUujRP5q3hrzMeNMW7nukDybj+i69xkOxY0c5VmuQqVOnYTAkoSh6DIZoEhImAHDvvQMxmUzoGtN/\n+h5Kh9H774drroFJk7pi6BqN5gSnBd7HsL/85S/Y7XbsdjsA33//PbIcfnz9u9/9jgkTJjBv3jwU\nRUGv11NWVsarr77KF198cTSHrTkWLV0Kb7wBV1wRrsgRExNuab5jB4wcebRHd1xQqyqQjEa8Zw7F\nf86oA9xZRap3IFIOzxOGkBJi054CdlWUUtVQh1/eV34wsV8SddEGAiYDcjACl7d7+OVJweXpjtOb\ngtPdA5urJ8FQ62kx2eluRubZGZVnY1SenYFZTkxGAaQ3vo7uUqHqaisFBRamTp2K0ZhEdHQysuwh\nNjbcFXn69OaN3bp0HUW/fuGXRqPRdIIWeB/Dnn/+eUpLSwGQJIlvv/2WpUuXIkkSV199Nb169eK5\n557j5ptvJj09HYvFwn//+18yM49OmS3NMSwurvnCrmHDYOdO2LRJC7w7IIRAKdqBaqzBtfMdhPnA\nWnfr88uIP+8h1NRE7D+83KVjC4ZCrC/YxaqdW3H79rVNFwLcviQs1lxWpV6FxZiLdUsWvrWx7R4v\nNdFPXqaLfpnu8NdebvpluugWF+zScR8oqcGJiDChs7sQMZGI+H0dO3W6OCIiUrFYPMTFjUSSJEaN\nyjp6g9VoNJp2aIH3May4uLjDbS644AIuuOCCIzAazQll+HD48MNw4P0bIskywthx4CyEQP38HwSn\nTEXGRjChGnH7hWBueyFjW5TeaUhuH4Z6R3jWuwvSU+RgkLW7d7Jm11Y8/nA+tdOdwu6ySVisudTa\ncvH4Ws6um40KPbv76JniI7O7j8wULz1TfPRK9ZLX001ibOsBtvHbdUhuL/L5Rz6dImbui0R8vAzn\n3+8n4sPvMH27Dsc7D7Jgu4XbbptLz55D6dZNT27uyUd8bG167TVIT4ezzuraajYajea4pwXeGs1v\n0bBh4a+bNx/dcRxB6j8WM+J3d1J5/a2I0aMRe3ZCXAK6lPA6CeFyIPR6gvogsmxB8u4h0GBCxEaH\nD3AQQTcAESaCYwZgWr4Z48qtyOdNOOhr8Ph9bCjMZ82ubfgCgfB1qTo27T6HX7ZfSkjZlyqSECMz\nItfOiFw7I/PsDMl2kproP7A4UAgi/7yEqCffBZMB29Ac1JyMgx7/wVB7huuaGzYVoi+oAMA08DSm\nZ0gkJvZFpzvGKvTYbHDffeH1E5s27fu3ptFoNGiBt0bz2zR8eHgmTm7Zivx4J4QIp4d89T6qWU/o\nlJMIhTzEPPpHdDqJzOefx/rgmZg3rEL07IViHoJOZ8T45lsEh+YQGNwTEDC96yoDBScMbQy8t3Q6\n8BZCUO90UFFnobzOQnldLQ0uB5OW7+HaHRaWT8rmh/RxLF9/BxZrFgDnjK/ijDE1jMiz0yfN23q9\n7GAIDPoO26lLFhsx972C+aufAfDcNRs1O73dfbqSy+XF4fCQPSxcOWXFih1scAT4vV5P5JBJTDEd\n5I3Q4VRRAdnZ4WZVp5+uBd0ajaYFLfDWaE5EqgpXXQUXXwznntsyyEpLA5cLovY1IBHWWoiORTL/\nqknIcUAIgbJqKUrJFvxnTSYUsiFlKQgjCN9O8AbQF+4rD6eqbnxThob/Egy3IA9cvjeNousXGFBM\niAAAIABJREFUCgYnhs9l+nEznjbG7/J5qbY2UGNroKqhjoq6WnxyoMW2ffc0MCC/jk96XMWSnc+g\nqjoyU7w8e+tWpoysa38gQhBz11+QfAFcL/8fRLe+yND4/QZib34Ond2NGhuF+9W7kWcc3lSOwkIL\nW7dWMHv2eej1Uaxfv4X168u577JrgCcYuSWf0YCUlwfHYtAN4cXLwcZ0nV81O9NoNBrQAm+N5sT0\n2Wfh2t0rVsAZZ4B5XwqCEIJQyAvvvYQydjxqTjaqGsD4t5eQp0zGMGgiERGpjc1Fjl2qrQGxaRXy\n2FHIcg1KDyeiew9EMBxgi9RuTdtKwSDe/7uIqD8vQU6MgUAQOrlI0vTVGvAHCU4cetD52aFhfVG7\nxaKmJBJyumlQZOoctqZA22KztltLey+dTkd6Rfgpxae2qyFW4uZzirj/8nxiIjtumqUrs2D6cg06\nlxf9nkqc7zyE2iu1xXZKVhpSQEaeOgr387ehZnSujX3n6SkpsfOPf3zHk08+jF4fQ48eLhyOAuLj\nhwMwY0YeM2ZcFF4pmpJC99ra8K4DBnTxWLrQ3XfDjz/C4MHhGW+NRqP5Fa1lvEZzDOmy7+kRI8L5\npa+8ArfeCjQuGKwpxxNtR5bL0VVUIxLjENG/nuGWiLBGEZlfjf7sKw59LF1I+DwQEUUgUI/cUIDu\n28/xn9X5nOniTVtQYiLp27eNjoitiJ81D+PaXTg+e5zghKGd3i8YClFRX0u9006D04HNaqXW68bh\ncXdqfyEgpKQj1GF4vLnU2npRVxbBhs8ykSUjJ8+u5Nn/28mIXEfHB9uPvqCC2KuewLCnErVbLK43\nfk9w/JCW2+0uR8nt2WFKSmfZ7Qrz5i3m3XcXYzIlEQiEKC8vp3///h3vfNZZsLdM6h/+AI8+2iVj\nas9BtQEXIjzOESMg48jmwmuOHq1lvGZ/Wst4jea3pqoqHHTHxcGcOU1vK+WFKH+dj3znpSBJqD1b\nznSGCQKiEqK8mIMejMboIzPudgghCMkeuPNiPHOvxfjBf1DOPAX5AIJuACXmANNohEC/O5yKEsrt\nuEynEIKqhno27dnNttIiAsGOc+iDIRMOdxoeXyahUDZeXy8aHGmUWlJweZunVJziXAVAQ0Y2X730\nM0bDgc+bKLk9cXzzLLE3P4/pu/XEXfgHnO/+geDpzWuTK3mHVpZUCHjkkY945JF5xMVlkJiYwLPP\nDiEysgeSJGEwmDsXdAN8/jno9fsOfKySpPBNgkaj0bRBC7w1mhPNL7+Ev550UlOKSTDoxh1bg9IY\ndHdEJCXgT0og4FxNVNQgzDYZXXrvwznqVilLPyWYkUKgZzTBYA386WYku5vY97/FvOQH7D8tOqzl\n2iSLDZ3Dg5oQg0hJaHM7j9/HluJCNu3ZTZ3D3vR+KGTC4++Gx5eIx5+I19cNj78bwWAKvkAKVmcK\nVmfbT/WS4gIMzHIyqI+TgVkupq3/FLZB4rh03AcRdO8l4mNwvvcwUU+8i+n79QTHDjyo46iqihCg\n1+uRJBOPPfYBN954JX369MNgiOeCC6KIicnF3LhuICcn5+AGrN+vcskxngKl0Wg07TmuA29Jkpq6\nNmo0xztFUbomr3pv4H3yyQhrHcEfPsM9uS+q6mkRtEh2N4btxQSH57a60E4IH97KH9Av/jf6x99G\nbzy8i9pEwI9w2gjFRyHLFtSIWkK4UINJ4Q0MekR8NLqKeiRVxfTNL8hnjD1s4zE0znYreZmtBnxl\ntRbW7NrK7ooy1P1mYmtt2azbcRFFlWOA9m8MDHqV3qlestM99OkRfuVkeBjQ20VKQqD5aceNxzY9\nE0wH1sSnVXo93j9eg/e+SyGy9Y6Vv1ZRUU9MTBSJiUkYDAncfvszXH/9VYwbNwmDIYY5c9LIysoh\nKir8lGTq1KmHPk6NRqM5gRzXgbfJZGrKiT3WF4JpNO0RQiDLctfkdz/4IMyciUhLw6/YCYkaVCWt\n1cAxbvZ8jBsLsH/+JKFxg1sfW0IMznsvw+DZgNncG1NpLVJcIrr0Xh0ORSgKFOcj9W1/RlUIgaLI\nhH78GLU8H+85pwACBvZsubHRgOfROcQ8/DciF/6z/cBbVcPXfZA/H/T5jYF3v+ZpFzaXk/9tXMvO\n8pJm79c05LF+52yKq8JpG3qdSnqyl5TEAGnd/KR185Pa+OfUbn769PDSs7sPg76Ts9cRJpRhfQ/q\nWtrUTtC9aVMRMTHxDBgwGL0+li+++IXx4yeRkxMOqP/+93cx7VdhZOjQzufAazQazW/RcR14S5KE\n2WwmEGhZckvz2+ByuQCIjW2/FfbxwGw2d80NZEwM6uhR+EQdPt8OOHVEm5sqg/pg3FiAYVtxm4E3\nAJJEKFQXfu3cjJSag6FbBCZTEpLPixQd/u8vVBWx+WfEkNGEFBchtwXj4kfxPXg/Op0ZnVfG9Nh8\nlOf+hk5nQtdgR7z5LN65NxAK1SBGJMKIsTSV9JODoAqIaD7T7r9iGlHPfYhx7S4MP+8gNKb1wN78\n6XIiX1qC7/8uIjD7tPClyCH0u8pQcjOapzC0IjQwC9+cWU2LKv2yzIptm/glfzuKqjZtV1k7iK2F\nV1BYEa64EWUOcfXMUm47t4i0pP1+PnkDmFZsRrI4CUztoOqFqoY7Xaa07EB5uFRXW3G7ZQYMGIrR\nmITD4SYmpg9xcaOQJImHHnqs2famY7Wsn0aj0RyjjuvAG8LltbTKJr9d27ZtA7TV5PsLbliB+tU7\n+G4+p8NtQ4P7AGDYVtTp4wcmNjYFcf+MJEUR+/TrqDfcA72yCAXtGL9YjLt7HSIifBPhu+8KCIY7\nDiKp+O+8BOFeA0igF+inDUaRy1qmwdTZibvuKZSe3XG/enfzz2Mi8V9/JlHPf0TkX/6Jq43A2/T5\nTxh2lyO5fU3vDb/kKcw1Nqy//BW1T4/2//uMH0xo/GBUVWXD7p38sGVDs7J/FmsOG3fdRmFFduOw\nglx/Zgk3n1NMcnzLhZW66gbirnwctVssgUuntJmfbvh5B9EP/Q0kcHzz3EHnsbvdPurrnWRlpSNJ\neoqLq9m5s4wzzjgZEJSUWNi9u5yZM8djMKRQXu6jvh7Gjj0ZSZK48MIunl3XaDSa37jjPvDWaDT7\nujX6fBX4e/kQ153Rqf2aAu/tJQd5Xi/Oey8HqQbcFgD8t5zf9g46HSJh79MJATpQ+rQsuyZZnSTM\nuBd9eS268mR0FitqWlKzbXzXn4lU78B/S+s3GJLDjen7DQidjsCZpzS97++ZjLnGhr6gosPAWw4F\n2VVeyk/bNzdbNCkHI9heNIfVW6aiCh3x0TI3nVPMDWeWkBgbbPN4anYPlIxk9JX16HeUojT+928a\nc72D6PlvEPHxMgCUHknoympRs9LaHede5eVWVq7M5+qrL0KvjyY/fw9ff72B+fOvQq83YbcXkZCw\ng27dpgNgs+0hMnIbiYmnIUkSs2a189RDo9FoNIfssAXeCxYs4KGHHmLu3LksXLiw6f1HHnmExYsX\nY7PZGDNmDIsWLWLgwINbUa/RaBpL7b36KIFBPfEPSgm/2cnFd8qgLAD0u0rDrcSNB/Ej4TBUFTF/\n+iP68lpCg7JwfPgIIq1bi21E9wQ8z93W5jFMX/2MFAwhTxiCSA2nawQVBW/v7sSvK8BQWEFw+kkt\njysEJZZqthQXsrOsBDnUPJCus57CDxtuwWKNQ69Tue3cPdxzSQGxUaGOL0ySCE4Yiv6j7zGu2Nws\n8Db9czkxv38dnc2FMBvx3X4B3jsubLO7ZHN6zOZsevSIIzs7mbi4EUiSxPjx/Rg/flbTVjk5Oc0q\ni+Tm5pKb2/ma5hqNRqM5NIcl8F6zZg2LFy9m6NChzXJWn376aV544QXefvtt8vLyePTRR5k2bRr5\n+fnExMQcjqFoNCes8Cy3is9XQmBaP9S4CPD4OxmoNR4jJgp50jBEfDSS04tIijuMI+4801drAPDN\nPb/VoLszzP9aGT7G2ePYtGc3GwrzqaivZRJuLgJKfvyZzWPTSYqNJykuHrPRyM6yErYWF+Lwtmzs\nHgwms33PnazcMgiAoTl2Xrx9C0OynQc0ruCkYUR89D2mFVvw33pe0/s6pxedzYV86jDcz9yKmp3e\nqeM9+ugnXHXVDYwe3Y/YWIkzz+zcfhqNRqM58ro88HY4HFx55ZW8+eabPPLII03vCyF46aWXeOCB\nBzj//PCj6LfffpuUlBTef/99brrppq4eiuY3wOB2oMYdG8HikSKEIGgpRnrqPtzz70FRGiAxGvwy\nSf2uQMlKw77sz52evXZ++ljHGx1Bks2FcdU2hF6HPK3ljHSnhBSUmnpUncSLohrrmpqmj2pSwzf5\nEcU1rNm5rcNDJcYkYLVfwIdLZ+DwmIgyh5h3+W5uOru489VI9hOcGF6oaVi1vdlTBv/VM1DSkwlO\nG92pKiySZCQiIpc77niYrKw+WmUnjUajOQ50eeB90003cfHFF3Pqqaeyfzf64uJiLBYL06dPb3ov\nIiKCSZMmsWrVKi3w1hywkNNG3w+eJf+6cwkEcjCbD25m9FgnhCDkd8Fdl+F7/DFCOieq3oV0+4UI\npaFpO8P2YqRAMBy0HUzKyLHCqMfz9C3oymsRCQf2JMzr91NQVc76gl1U3DicRFsetoh9AakEWFJi\ncMSZcce0XUYv0mxm6kaVUGE8z3ErP/hHAnDaiFqeuXUrvVN9be7bEbVHEr5bzyXUrxcoKuzNCtLp\nWk19+TWr1cVbb63kj398ErM5ntxcLeDWaDSa40WX/nZevHgxRUVFvP/++wDNZmBqasIzTqmpzdtU\np6SkUFVV1ZXD0JzgVEVBDlrxhnaw+5rTAR/u+p/Q/e0LDPe/iBRxgG3Bj0FCCJS//JHAmTORE1UU\nxYl0z+UIUQGqBBL7LVIMM2woACA04vjO2RUxUfivmdmpbb1+P6W1NZTWVuNdtx21zMKOQfsWItoS\nowAwG3pQabmU1VsH0i3WQ8W9LnqlltMvsQijqRSby47b5yUpLhWX+3S+3zCKaR9dwSzHEl7pfx69\nBnj4/eW7ufDUyi5pnOh59Pp2P1+6dAMTJw4lKioGnc7MQw+9zv3330RiYjd69Ihh/Ph4TKY4bZZb\no9FojjNdFnjn5+fz0EMPsXLlyqZOknsrLXSkvV8e69at66ohao4mIZB0OnQ6hW5rl+PplYc/JdyU\nJKZwG/7kHshxHdcrji/YRHzhGoqnD0YItWlhX0F5ATFDUwnt+oVQKPqwXsrhYrLVIUmgJJtR1TpM\n6SreipWotfstlGywtLl/zo/riQGqeiZQW1hw+Ad8FAghqHM7KbfWY3HYsfvCudi9ymzc+8KPOOLM\n/OmP0wgZ9EiSRIQ0kC0F57Bs80jkUPjHXWltMhv3AIQreERHBOiXUUtaopOfdmTj8IZv3N7z7QDg\n4htKefCU19HpYM+err8mSdLxr39tZMqU8SQnZwDRrF1bT48eccTHJyKE4OSTT6eiQk9dnR/wk5qa\nyoYNG7p+MJoWtN9Bms7Svlc0QIcL1rss8F69ejX19fUMGjSo6T1FUVixYgV//etfm+otWywWevbc\n143OYrGQlta5Ulma449OBz2++5BAchz1Q7KQZTu+KAuyLoSiVAIGzLVbkaOz0CWkI0mR9P7wY+on\nzMTTsx9CGDFVlyEndkcXpeDMUnGkZoeD7l9x94hHCmzFbB6E6jWhmjrXBruzpGCAtNVfU3faTBQ1\nAqGqIMRBV/Uw2eow+D34M7LR6fzEF69CNgdoiAn/e5Azkjo4QnMxO8oAcA/quKPk8cbh9VDSUEtJ\nfS0euWXDrLLMBCrT48iocjJ5Sz2rR1/Kyi0zWbUzFyHCN/bjBhRx0YRNhBQ9uypSyK9IJb8ihXpn\nDBv27OtMmdOjjtlD15C+qhrFbKTfKb5DrtzidvvR63VERUWh10fx5pvLmDDhZPr3HwxEk5FhQJKy\nUZQEAK644hoA1MYmPYMHa2X+NBqN5kQgic5MSXeCw+GgsrKy6e9CCK677jry8vJ48MEHGTBgABkZ\nGdxxxx088MADAPj94Zmb5557jhtvvLHZsfaKj4/viuFpjhAhBKFvPkK1V+I/4zRCoXpwWhFRER12\nCdxL8voQRgMYzeh0kUT97RPkc85GTlZp6mjYqLBxZrdv3313mJLdQ9zCT9A//xE6w6HdW6qFOxAZ\nvZFxIAdq0H/1T/ynjUIyxmGslYlc/AbiqTcwGmPCQbgkNT3BUatKoaYCdehoVNWPuvp/6LZvQL76\nGlQ1gG7bRqipIHDqKIQItri2A+KXSZwwF53FRkPRhweU4y2EILApn+AP6yntm4I6uA9DsnK6LI1B\nqrEi4qIh6sBuhFxeD9tKi9havIcaW0Or2+gkie7xPdBJQ+n7QxkXv7mA4rhc+g7aiZB0mAwKs0+r\n5OZziujXy93q94vFambznniKqqIZM9DK8L4OjOt2kTBrHqGhOdi/e/GAr1kIgSwrREf3QK+P47HH\nXmPKlCmcfvpM9Hoze/bsITU19YTouHqi2jt7qTXn0nRE+17R7K+jGLbLZrzj4+NbnCAqKorExMSm\nOt133nknTz75JP379yc3N5fHH3+c2NhYLr/88q4ahuYoURtqCcVG4veXEOxvAJGBCDbeiMUeWOqH\niNqbo62iqh7cc2YBSuf3T4jGcd/lRAYqidL3OuAAUqgqqhDIch3Sv17FN3MMSlr4ezs4bUx4G+FF\n7g7yPZeCczk6XSzmbeUYV64icMcdqKofLPnoSgoJ9HKEx58jkLKGIfz54RPlJYZfomWHwwMWYcK2\nbjGSw91h0G1zOSmtrcFit2KxWam125jx4VpOW17EtjMH8K0zj417dnPeKZOIjz60Mp+6yjoST7mV\n0NC+OP6zoP1qHaqKz+dnV3U5W0v2UGKpbnUzRelOg+08Ghx5lFlS2FMVi6pKGNQgZ0sv0MdZwGN1\nj2K9/XKum1VKamLLGfL9pXYLML1bbbP39Pnhpwehfpmt7dKhzz7bjNWqcv/95yFJEs8883Kzz/v2\n1TpCajQazW/RYS19IO03+wcwb948fD4fc+fOxWazMXbsWJYuXUp09PGZk3s8E0KAqoCuMR+/shSM\nRnSp4S6CwloLOgNSQvuVQoQQBG1V8KfbcD50A+iBGDPQtWkeByzChM+3FVVxE/HqYnQ3/wF9Qvup\nG0IIQv94nmCcCd/EgQjhg2unt7sPZhMgUFUnvgFx+PqeDoHC8GeZ0ZA5jKabhgjzocxpd4qIbztQ\nFkKwcvtmlm1e3+KzyozwjUXPyvCdeqmlmte++IxZJ41jSJ+cFtt3mhxC8skYf96B6eufkc8Y22KT\nYChEQWU59V/+yIyn/o06PouSM5s31dLr9EQaxrFx9yy+35BHMKTb7zOVfpluBmY5+THtVqb/6yXu\nSv8vnivamX1SFPT55egq6lqtJBI8ZTDuZ25ByWq/s+X+amqsZGTkEhGRzVVXnY7JZNIWP2o0Go2m\nmcMaeC9btqzFe/Pnz2f+/PmH87S/OUIIgss/g93b4drfoddHIdVZwOVAlxvOuRe2eoTPg0jNQFE8\niI8Xo0pB5LNORwgV47oVEBVNaMwYQML47VKIjkadPB2dLhLD/5YiJfdEN+Z0dDodong3SowZX4Qb\nWS2Dh+dA1zcwPESCQGAPyrgsVGUzRlcaRuIxLl2K7oJrkSQJtWQ3SuFWAmOHEQxWoU7KRESaQRxE\nuThJ6nTHyCMtGArxnzUr2FZaBICq6nB5U7C50rG70qlxBrmS24jdHckXKx/gtNGvAA4+W/UDBZVl\nzDp5HBEHkTOv9umB+4kbiXloMVFPv49n6kisHjdWl5MGl4Nam5XdleXIoSDnfb+NaG8QQ2hf/n5m\n9wycrhl88/MprMtPBkCnE8waW83Mky0MzHKSl+kmwtS4jzoZ10w98oSh7Q8spJIw5U4Q0FD6MUSY\nmo87Jx1/Tucb0QSDUdx330KWLPmMyEhtIkGj0Wg0rTuOi/3+tglVRWxbh5zXh0CghGB/HVLPHIRz\nBaDDWFiOodaJkirQ6cxIa5eDoxb/lFEIIcPUHDDoIRgu8xia0C984Mb0kODUxhnHxrQIXRZgqkPY\n/4dOF4N51f+Q+ySHaxHDYWkb3iUkidCAHBAeAoE9yG4vZrkE1bUDgyGaoD8fQoUE/Y2BV0zU0R3v\nYeDyevjox/9R1VCPzx/H/375HRW1w1DUff/8jaqMLP0f2d4iassH8on9BWae8hSpSQVsKy2irM7C\nuadMok9a54NRvyxTYqmidFAC53SLJnZ7CV899DSbhrdyDCEYujWcWlI1rj8n5U5me9FkXv4oj8r6\ncOpRXHSQy08v4/ozS9quo63TEbjw1I4HZzai9k5DX1SFvrgaZUDvTl/XXvn5FcTEdCc3dwwRERl8\n+eUUbYZbo9FoNO3SAu/jkBCCgLcW6d9/xX3NGRBpBoO+sa6zABSC/dMJ9k8HuTi808hUIHVfPvEB\nNlhR01MaTy6jKFa8M0d21eUcUSImCv+00SAXIctAMpA8sJ0dBHgD6BxuJIcHye5GxEWhDOpzpIZ8\nSKoa6vlo+be4fF6szp78d8VDOD3hqinpST5yMtzkZHjom+HBU5tFYlEBl3T/iTesM/jshyeZNOI1\nBmZ/h9Pr4R/ffcWw7FzSEruREBNLYkwsCTGxmAzhWX5VValsqKOoupI91ZVUNtQ1lRM1Tslh9pIt\nzPx6F5uG9WiR693fqZLc4MUbn8jH7tf45qkehJTwzVxOupsbzirmkikVxER2Pte/I6HcnuHAu6Di\ngANvvT6JwsJKMjO7MWzYweWBazQajea3Rwu8jyPqpjUEdTL+3lEEgzVwy3lHe0gnPH1hJYnjbmvx\nvu27F1GGHkLucxfSFVWhc3gIDcwC875Ulx2lxfxr9XJCikJF7WC++ul+AsEYhve1886Da0lLar7o\n0Fg9Ga9tFE9cWIH4XwlvfpXF9+tup8HRj3FDX0evD7G5qIDNvzp/dEQE8dExNDidBIKtLxRdPbY3\nmRV2VkzMJjE2jqS4eLrFxtEtNh5FyUG3YAUA75sv4Iu1Geh1KjNPruHqGWVMGVl7WB6oKLkZ8A3o\nCyo6tX11tY3//Gczd901j4iIFK65ZlzXD0qj0Wg0JzQt8D5OCCHweUoISvWEggf+WFxzcNS4aITZ\niIiPQU2IQXJ50Vc3YPpuPb5jJPCOeOcbohZ9hvfeS/HefzlCCJZv3ciPWzcCsLPkNJatvQ1VGDhj\nTA2v3L2R6IiWM8f+m84Gwutjn+67jeF97cx7bQibC6bhcGczefQTxETaWuzn8fvx+P2tji05rgf+\nwCmohp58dnkcwVAMzmoz2wuNOD0GLLYIthbF8+ei5QQlAz/nzuDhq3ZyyWkVpHZrvxrJoVL6hvsJ\n6AvbDrzdbh/R0dGYTGlkZg5nyJBEIiNT29xeo9FoNJr2aIH3cUAIgcu1A7lfNKAt3DqSREoCDRWf\nNv1dv3kP+opagqcMamevI8u4YTcAwRG5+OUA/1r9I7sryhACftl+GWt3zAbg1vP28Merd3a2nDqX\nnV5B/94urlswmpLqHP6zfCFXz/yW7olFGE0luLwO7B53U5MXgLioaLrF9qe69hROe/1NNrmTeTX1\nbByGhDbPE2UOsfqG+8gdO4snh/uQzIehPWSj8nIP5eUW+vbNJTS4D/LEoSgDs5ptE/Ha5xjX7EC+\naTbX/PVLFi5cSJ8+eUiSxNlna2klGo1Gozl4WuB9jFNWLiVQsRH59HbykI8huoo6oh/4K77bLyA0\n5vgYc7t+lYusDMtBGXZszHQDEFIwbA6XL6zO6c6HX/8bq8tJSDHy/dq57C47FZ0kWHDzNq47o/SA\nDz8i18H/s3ff4VGV2QPHv3f6TCa9Q0IKSUjoJYQqoNLBgoqsrmL7iYVV1F117ViQVXddRXFX3ZW1\nrn1VFBVUEOkd6aGHkEJ6Mplk6v39MRKMhHRIAufzPHkS5r733jPDhZx559zzLnn+J25+bgCrtofy\n9w995U16nZekzjaSYyqIjy4mMqSA3MJQFq2IY8u+IALcZTyzYyyjFR1fDL2NpORyAvxcBPq5CbK6\nfvnZRaDVRY+4cqyW4zPwrd8VRlVVFEWHyZSEXu/GYgnFah2CZ3APvvmTH8XFeUxGAVQWLtxA1yXb\nGLF8HYbr7ubzzxeha+EiTEIIIcRxHe43iqqqeDxOdLo27hN9Bni9XirjLbiDO85H2+aXPsH4zTqM\n36yj6OD7qB2sS4h2636M367FOaw37mHtf5lubeYRFLuDqk4hvLZ+GW6PhzJbJIvX3k1+UTf8TG5e\nv3cTo9OPNXywUwgPcvLR42v45xeJrNoeSuYRK0eOWdh1OIBdhwOAzrXGmw0e7gj/Cg0q7j5JfDnv\n5L7hZ9K9977F/fc/RGpqMjEx5cTExGIyhQFh9Oo1BrvdTmBgHG53JT16hBOa/Q4ASvfuknQLIYRo\nVR3qt4r3wG7sAS7cmnL89zrQDrqgrUM6bTyOamzVmbhMJWCqf+GX9kS37UDNz6b/fEPVHy5rw2ia\nzvDTVizPvY9SYusQibdmo6/d485IMy63h10HL+SnLTfhcpvpFFrFO4+so2dCRYvPo9ep3HHZfu64\nzFcGYrNr2XvUyp4j/mQesbIv20pEsIOxA/M5r3choc995tvxvLS6487Kx/yvL6l8aHqtG0Jbk6KY\nMJm68cgjz5CQkFBnq7+EhBPdafR6Kxm9AuDgQdBqITn5pPFCCCFES3SYxFtVVVxLP8KdFoy7SxSu\nn5ZC70Foz8LFKjw/r8W18N+4brmorUNpmmpnTdmD/U+/o+rWS9o4oKbT7vKVY7jTurRxJA2rsFey\nKXsvaakR7IiP5auVD3AoJwOAi4fl8Oyt2wgJcDXpmIZFa9Av34rj8pG4B6aecpzV4qFfchn9ksvq\n3K5fsxMA1+A6auFVlYAb/oLu5/144qOpvnFik2JsSFZWAe+/v545c55HrzfTtSmVQXsKq5s2AAAg\nAElEQVT2+FpIJieDwdDweCGEEKIJOkzi7XAUYLtkAMeX37ZdPx6D6wD+pp5n1aIVHo8bWxc97t9f\n2NahNJmmtALn6HSU8krs91/d1uE0i253FgCebo1IvKudYNCd1sWDyu2V5JcUU26vpKzSRlmlzfez\nvZLyShveLmYWXDSLpRtup8oRSIDFyV9u2cHlI4/+tjy9UfQ/bsH8xiI88VH1Jt71qnai25yJqii4\nM+qY8VYU7HdNJeDGv2C9/594QwJwjs84afXIU/F6vezenU337vEoihGHQ8stt8zlvfdeRqv1IznZ\nwKWXpqLXm5se+07fGwbS6p6pF0IIIVqi3SfeapUdZ9ZOKsOKOZ50H+d0HqLyaDWWSiOabvUvEe3d\nuhYCg9HEp5zGaJtHVVXc772EK9RM1cBEVNXeIVdQ9EaFUvHmg20dRvN5PGj3HvH9mFp/4m299W8Y\nF66k9Ju/4umV2KphVNgrOfzTWtyL17Kwbziq5kQGHX+omMNxwaAoOJwWVmy9gV0HRwMwvFcBL83a\nSufwulv7NYY31rdQkvZI82vCMegoXfQcut2HUYOsdQ5xThqMO6kzun1HCbj5WYq3/Btv5/BTHtLt\n9qDRKOh0AWg0Icyb9x/efvstTCYrqqry3HPRWK3JNW/Chw499bFqKSuDRYuguhpuuAEmTYIffwST\nqclPWwghhGhIu0+83fu24lnxAerUumeA3YfW4MxRMab0qnPmW1VVXK5yPJsX4RozAavXi6YdLG+u\nlhbjzTmEMyEKpzMHd3ooqtUEqv2MxqHddgDj12twDe2Ja3j9b17OdprD+ShVTjzRoaiBdSeMNbQa\nFKcb/artrZJ4V1ZXsSvrEDsOH+Bwfi53vLyS5P1F2PN68u6g8ZRVRjJw6yHu+f5z/ttlCjO7vkyp\nIxLQYNC5efT6PfzfpIMtnnz3/JJ4a7ILmn8Qjabh7i8aDfZHriPguqdxDUytN+nWaPy5/fa/8uij\nD9K792A0Gg0ff/y/mu2KopCa2szZ+eJiuPpq6NTJl3gHBsKIEc07lhBCCNGAdp9422I8eE6RdAO4\nU+OxpWpRq4/VWtjCm30Qz8512Ael4nJlw+R0oICKiu1YNXFo/QPPQPR183icVB9eC5uXUxXxy+p3\noQFnNAbt1v1Y7/sH2oO5aEoqqLphwmlNvJUKO6rVfFJ7vvZEDbRie+ZW8Hjr3L47y8qLHyexeW8Q\ntxTt416WkvPZQXYNDSKhUxlBftWoeNFqtGg1GnRaba03g6qqUlFlp6SinBJbBcW/fC+pKCe3pAhV\nVVFV6L6siuT9RRQagrmpbAlF3/uua1Px59g1b3BV1v8IKq/i6tT/ktrTzV9v20a3LrZWeQ1aZca7\nkZwTB1P26ZN4unb+zRaFd99dSWJiN8aOvQSDIYh33snAYjkNnwJ16eKb3c7JgfJyCDiz/w6FEEKc\nW9p94u3xnLxSXh2jsNu3oN2pQ5eagVsPDuch1Op9uFy1f1m7XIdxv/w3mHYX2rR+pyfo+iL1uLHZ\nfsbV2Q2d227JacPKn9FvysSd1BlNSQX6VdtP37n+9xOm+17jk2lP8J71CsKDHCRGV9K1cyWJ0ZV0\nibSj16lNOqbLrbBqeyjvfRNLXo6Zqwdnc+nUaoz6upPmxlBDA+q80W/bgQCe/zCZr1ZH1zz2WtVk\n7uUegrdsZcpDw0BRMBoqCAnIJiQgi5DALEIDswgLPIK/XyVajQa3x4Pbc6JcSlUVqp1WqhyBFJQM\nJyuvH1VHOjN7je+6uDvuRUp0EXSPKycxupLYqJ58ZnuVy1+4kwml33DMNoDyWY+iRoU0+zn/lifG\nN/PcohnvJnCd1weAXbuOUFzs5sILx6PXRzBhQgohISEYjcEApyfpBl/3kpQU+Pln2L0bMjJOz3mE\nEEIIOkDi3Viq6sS1+TscRgeOKCDAA0PrnsGtmHEJOosLq9uJTnfmOhd4d27G9fV/cF03+oyd81R0\nq3YAUHX7FKwPvYZuzxGUglLU8FOvMNgUZTYdP24N54dN4XT5dB9Pl5Zx3oJnubXvdVToomuN1Wq8\ndImsomsnG0kxNlJiTnz/dVeOaqeGZZvD+WpNFN+ui6TUZuCKwo/4MfNqPlxyJRnf/Ys7L9/P78dk\nYTI0PwE/bsOeIP7+YTJLNkT+EqeT7onfkRq3jMqqIAr2BBFuL+Q8wyLWMgqH05/cwjRyC2vfmGc2\nlhESmIXFWEq10x97dSBVDt+XqtZeRnLB3usJdRezLWE4o55J5dG+iwn2/3VnkiAqJ/4F3bTZ6Lcf\nwPL8h1Q+eVOrteRTw4OofHg63phwX3eP0/QJhdfrJTe3hC5dEtHro7FYQrHbK7BaffdgpKS0znXY\nKKmpkngLIYQ4I86axBug6tIR/PYGzDoZ9Ljd+VRUrMeapaDrPhBFf3p6CR+nqipVMf44xp/5WfaT\neDzoV/sSb9eovrgGpmFYvhX96h04Lx7WrEOa/vUlqtHAutTxPPLpINbvDsbj9RUca/zv4Iqg9+hf\nuoGlQTP4aupDHMzz50COHwdy/DhaaOZgrh8Hc/34bmPtxYJCAxwkx9gI9HPx07Yw7NUnLtlOYfno\nIw6iy/QwpHIluUVmHnitJy9+nMQdl+3jmrFZmI1NT8DX7Qrm2f+msHyrb/ZXp62mZ9dv6dvtc6zm\nE5/AZHXzw7Knkqu6zaF/zwU4nKEUlcVSUNKZgtIYispiKSrrQpUjkKPHetV5Lqu5mpCAauKjqhnT\nK5cpr25ErTDQ+f1ruCQxr859vImdKP36Ofxv+xvazCOga617FjQ4XRbeCQ3n6innoWlO0l3tbLA7\niaLoycvz8uyz3/DOO/9FURQGDGjDntnH68N37267GIQQQpwTzqrEu6k87iLc338DIWHoY09vt5Pq\n6nyqnHsgsu0Xw9HuOISmvBJPlwi8sRG4hvbwJd6rtjcv8VZVzM9/iLaglDv738YuUyg6rZehPYu4\noP8xLuh3DPZfgPeajfRd8iFLe+aT0CueQT2DCQ8MJtAvjGpHFwpKIth31J+92daar6JyI0U7T6xS\n2rtrKcN77cWr+QSdbg8ajxfHIi0JVVlM6/sI3x+6kbziBB76V0/mfZLEzCn7uXr0EQL83A0+jQM5\nfjz5VmpNSYlBZ6dX8iL6Ji/EbCoHoGt0DCN79yPEGoBu3GXYgvy5/De13D4OVHUvHs8eso4Z2HXY\nn1KbkehQL+FBDsICnYQEOE8qjXFePofS7QfxJnaq/yUPC6T8oyca8ZdTHwVF8eOqq55gwYL5REQk\nUFFRhZ/fHozGWFyuo00+ov8fXkC3cQ+2eXfWlJGA743nE098yIMP3k94eDIhIUbefbd1+3c326hR\nUFICK1b4kvAHH4Tp09s6KiGEEGehczrxRlGonD4BjSYHS3UQBrcJxWxB0bbey6J6vXievpOqaedB\ncPtoEahf75vZcw7txYqfQwlJHUefJ/zwXti3Wccr3lREWEEp+foIdhlTuf3S/fxx2l6sZhe7jhzi\np21byC8t5pKRiVy4dD+Tv9zB/C5BFJSVAgdPxKXTERcRxfA+XUiJ6YK/2Y/cIhN7s63klxhJ71bA\nvtyfWLN7O8fneL1aDblxIcTvLWCk6zvCxmznYE4G63dcSX5JVx59owdz30ll8tBcfj8miyE9ik+q\nnigsM/C3D5J585s43B4NOq2Dvilf0Lfb55gMlQB0Dgvnwr4DiY/8VZlMAy3nFEVBp9OS2MlDYqfS\nxr2YGg2e3k1Z8aXp5s79lKuuuorevQej1/vz9tvdiYiIQFEUQkLM3Hzzrb5PaKoiWL36S7p1C8fU\nmB7bqopuzU60+cV4I0PYuvUA0dGdiYlJw2CIZOrUIIKCktDp2lmrvvPP93397nfw00++EhshhBDi\nNDi3E+9feL0V2GxrsXy1Hl1ADPpLb261RXncnmoqJw7CG9SMxTx+Rbv9IKb3v8c1vBfO8YNadKzq\nGydSMbQ/T/4nlX8+MgQYgtlwNf2zSxjco5jB3YtJTy3Bz9Rw2c6SDRGsuTeHF4H1IUN4f/Y6RvXN\nZ1fWIZZv3/xLcu2zeEw3ovJtLB1Vd2LpcrvZl5PNvpxsFq1fRVRwKCkxXUiJjaVrjMoXq5dTWO5b\nKbH3zzkoRiN+YwdjGFgKe39kcLmG3VoNiZ3XkdBpHYdy09maeRHZx3rz0bIYPloWQ0J0JVePzmLa\n+dkEWl28tjCBeZ8kMThnKf/Nv5c1qSkUTcyrKSkJDwzigj7ppMR0OSsWatLrO3P99XeRkpKKweD7\nJCEyMvKkcYqiYLHE8MMPhwkN7UxMTMPHduzOoriglNDQQHQ9z2fve0sIDIwnIMC3euWI9t6m7/ji\nOd27t20cQgghzlqSeP+KfWI6eL3oytZhMiVi2J+FktoXRatteOc6eL1uKit/xt2lmTeKeb3of9iM\n+R+fYVi+FQDNoTyc4zJadNNbYbmR6QsuZ8OeEPwtLqJDqsnM9mfl9jBWbg8DfDc89u5aTo/4crp2\nttG1UyVdO9uIi7Rj0Ks4XBqefDOV1xYm8u/8fwIw4KYoDkVu4p9fraOwvPYMr16ro++AXlivvY4p\nWi0FZaUcKy2hoLSE/NJiCspKqKyuvfBLXkkReSVFLN+2ufYTUFUuX5RJcF4Zu7p1x9YjDm9EMMlx\n8dxx8UWs2rmNTfv2kNBpAwmdNlBmi2TXwQvJzLqQg7khzHk7jb+8240gq4uicl/yOUn9hMuLPiXA\nlMxCc3esJjMX9E2nd0JSu+j73lJbtx5k4MBxWK3J9OrV+OczZ87cX2a/cyks3MSxY3nEx0cBGrZu\nPUhhYSVjxw5HozHy6Zvvo7NaueG88/EPSOPWWzvQ6o9ut2+5eDhR8y2EEEK0Mkm8f01RQKvF7T6G\nrSQH68efw6wnMAbGorpdKFV2lICGk2hVVfGu+YHq3K24RzTvpjHtz/vxv+15dJm+lRRViwnXwG5U\nvHF/i5LuAzl+XPVEBgdz/egcVsV7j64jLa6CwjIDa3eGsGZnCGt3hvDzgUA27w1i897az/d4BxKv\nCofz/NBpvVzC9wCsDavkh2WLa4036PQMTEljcFpP/EwnZv1jwiKICYuoNbakopzMo0fIPJrF4WN5\neL0n3xhp0OmZqo8gOK8MT2QIZenJoEDxzZeDohAITBg4hPN69mH1ru1syNxFoDWfwb3eI6PH+2Tl\n9+XAkYnsyepLUbmRTmFHGNjjdUZ89SUAudH+9Ouawuh+GZiNxpPO39q0P+/H2ykMNez09ZVXVT2f\nfLKN1NTLCQhoRNJ95Ag8/jj4+cGLL/4y+92JzZvXcvSojT59RqEoBqKi9mI0lhAUNASAGarJ1wv7\nvPNO23M5bQ4cAKfT19fb37+toxFCCHGWksT7VPQ6bDMvB/fP2Ev2ocuzYXr7A9yPPItOZ0VbUY1y\naC9K/2F4vW48G5fC2qU4rr8Gj8eG6p+DGlT3jZSGhavQ7TgIXhW8XqpuufikNn7emHC0R/LxRIdS\nffNkqqePa3g1xQas2xXM9DkDKa4w0LtrKe88tJ6oUAcAYYFOJg3JY9IQXycNm13Lpr1B7M22sv+o\nlf05fuzP8ePIMQsHc/0AiI+q5NU/bqRgwyS2LfqRZWoZ/FJ9bdTryejWg0GpPbAYG1fTG+wfwKBU\n3z4Ol5P9OUfJPJrF3pwjVDkcxEVGc/Hg84iZ/TYAjqkjQftLIvmbNyNWs4Ux/TMY2r0Xq3dtY/2e\nXbhwEx+9ifjoTQzuHUiFPYyI4AMoikp0ru/myf6XTiRqcBNayqkqmsP56DbtwXnZyMbvB+ByE3Dz\ncyjHSij7bA6ePklN278RtNpg/Pz68Mor4xu/k6rCv/8N0dHw4os1D19yyZRaw3r06FF7v+Ji399D\nR0y8MzN936XMRAghxGkkiXeDVLzeSpwRCs4//g7smwAFbW4xxl0HcCRq8HjKIboCZWJvVMcB325R\ndc+Ma44cI+DGv9R6zDF1FJ7fJN5qSABln8/F3TMB9PX/NTldCsdKjfiZPARZXXVOiH+xMpqZf++L\nw6Vl9IB8Xrt3E161nNW7DuBwudBpteg0WnQaDTqNBq1OT0SIlpRYf8IDg9D+Um5T7dRwKM/CsRIj\nvRKPsWLHKhaquTDhRFeYnnGJjE8fgqWBmw/rY9Qb6B6XQPe4BLxeL5XVVVjNFhSHC+NnP/letysv\nAJz1HsfPZGZ0vwyGpPVi9a7trM/cicvtxmIqw2Ly1Yvr3V7Ci+yoGg1Rg5p4g6mqEjT6bjRllRQP\nTKtZ+bExTG99i/ZADp7ETni6xzftvA2w2x388Y9v8sYb72AwNHEGt1Mn38IyubngcEBjZ/7/9z9f\nd5COOGM8eTIUFEBFRVtHIoQQ4izWaon33Llz+fTTT8nMzMRoNDJ48GDmzp170qzY7Nmzef311ykp\nKWHQoEHMnz+f7h1ulknFEx2MPXoAeIp8D5lNqI24f1K3wddRxJ0Wh+OS4aBR8IbWXWbg7neiTEVV\nYfH6CLYdCCS32ERekYncIhN5xSYKy04kRhaTm5iwKjqHVxET7vteYdcz/39dUVQvDwxYwu1/rmLT\n/h38+PNmnO4Ti7OMWZJJxspDfHxFb7b3jKp5XKPREB4QRGRwyC9foUSFunnzu1WU2ytPnNtoYlLG\nUNK6JDTqVVQKyzB8vQbHtePqHafRaPC3+GbZDd+uQ1NWibt3VzxpcbBvb6PO5UvABzIkrWetBDw6\nJIypAV3QeBfiTurUYA/qOoLDNbgHxm/XoV+9A0cjE2+lwo7lr+8DUPnI9AbfXDVVSEgvnnjiWfz9\nm1HCotNBTAwcPgxZWZDchHKp4OCmn6+9CAvzfQkhhBCnSav9tv/xxx/5wx/+wMCBA/F6vTz66KOM\nHj2anTt3EvzLL+NnnnmG559/njfffJOUlBSeeOIJxowZw549e7BaW1ZG0VHoN/kSRcfFw6j647RG\n7eNwabj3lV68/0NsrcfDncco1Aeg1XoJD3Jiq9Jiq9KTme1PZvbJs44vjVnI7U9MIXdhCEvuPrkc\nwOD0EFJaRdK+wlqJt9frJb+0mPzS4l93/6slrUs8EwcOrVXHXS+Xm+CRd6I5VoInLQ53euNuaHNe\nOICKl+9C9W9ea8bjCfjwHr0psVUQGRSCpspJ2WdzwO5o1jHdQ39JvFdtx3Hl+Y3axzzvEzSFZbgy\n0nBOGtKs89bl2LEy4uKGYrEkkJbWgi4scXG+xPvw4aYl3kIIIYQ4pVZLvL/55ptaf3777bcJDAxk\n1apVTJo0CVVVeeGFF3jggQeYMsVXK/rmm28SERHBe++9x4wZM1orlHat8pHpOC4bgbeRN9MdKzVw\nw9x01u8OwWJ0c8PEw8RG2Lnwx3/R993XOPK3+6kcl0xJRQl6nR6jLowyWyg5RRayC8wcLTBTUGpg\ndPphun31EgDZYSeS1vDAIFJj4/F4vejsFliSSa9sG3vjEnG6XBSUl1JqO/XH72aDkYkDh9I9LqFp\n7fb0OqqvHo3lhY+wPPc+5R/Mbtx+VjOOaRfUuUm76zC6DXtwTB3V4My1yWAkOuSXTwr8TLiG1b2y\nZGO4hvQEqFkNtCFKeSXm1xYCUPn4Da22LHt1tYd7732b//734pa3PoyPh+XLfYm3EEIIIVrFaavx\nLi8vx+v11sx2Hzx4kPz8fMaOHVszxmQyMWLECFatWnXOJN4Y9LVKSOqz7UAA0+cM5Gihmc5hVbzx\n5zWEBx8kr7gI96pdaB0OzE/+nRcco/BqT3Sr0Gt1hAYEEhoQSGLnQLRaLat2/kzftb4yl31Joei1\nOkb06sfg1B419duk9EJ99EPCDhRwRe8M1ABfeUe10+mb8S4pJr+kiPzSYirsduIjoxnTPwOruXmz\nz1W3XYLp9S8x/LAJ3cY9uAd0a9ZxjvOf8Ry63Vl40ro0ega9Nbh7JeL1M6M9kIOSV4waFVLveDXA\nj9KvnsHw/cZWi1NR9ISGDuTzzyej07XCP+tbb4VLL4WMJtxoKoQQQoh6nbbEe9asWfTr148hQ3wf\no+fl+bpl/HaxjoiICHJyck5XGB3WwlVR3PFCX+wOHendivnzNZ+wePMSHC5fTbY22UximB+R+RUM\nWpfF6iHxNfu6PO6aHtg1VJWk/YUAeIf34faLxhPo95vyHj8T7n7J6NfvRrdmJ66xAwEwGQzERUQR\nFxFVa7jl2ffQz5mN/Y/TcF04oMnPUQ0JoPr/JmF58WPfrPf7jzX5GL/m7p+CbncWuk17z2jijU6L\n4+rRoNWgeL00Zt1DT88Eqno2rha+Pqqq8vXXW5gyZQYmU2TrLfIzpAnlL1lZvnZ8/fpB4OlriyiE\nEEJ0dKcl8b7nnntYtWoVK1asaFQiUN+YfY28ee5soarwn+8G8fo36QBMTN/BNaM/5MdtP9dK6Dw6\nDV9OTOOGtzYw6Zs9ZJ/fk2qdhvJqOw63+6TjRuVX4G9zUhVipeuA3hTk5lJQx/lj0zrTacMeStZu\nJS+x/p7l3b9dg+XnQ+RmZ1O6L6BZz1c3vhf9XvsCw/cbObJ0NY7Y5t3ctm/fXiJig0gEqn/cwL4L\nzvAiKDeO8n23l8C+klY9tM1WjcPhJjTU90bJ5fKgqioGgw6wsGlTPrGxu/D3z27V8zZW5LvvEvvC\nCxRMmcLhBx9skxiaasOGDW0dgugg5FoRjSXXigBIbuC+qFZPvO+++24+/PBDli5dSnx8fM3jUVG+\n2dL8/HxifrX+dH5+fs22c53TreXJ/47j+y3dUBSVmZN+YvSAH1i6Z3tN0m3WG4gICCTYYiUouQfl\n63IJ2H2UaZuLybl+NADVLhcV1XbKq6oor7ZjdzpIqTZR3isOR0x4vTXFub8bQc7vR+EJqL98RHG4\nsO46gqooVPRu/sytO8jKgT9PpSoxqt6k27I3h6rYcFST/pRjbN27AGDdkdXseNqD0lI72dml9O+f\njkYTzNatO8nM3M811wwDYN26NWzfvofbbrsaCOSaa9LbNF7Lbl8JU6Ws+CiEEELUq1UT71mzZvHR\nRx+xdOlSUlJSam1LSEggKiqKxYsXM2CAryyhurqaFStW8Ne//vWUx0xKOns6KijF5aghdc8M2x0a\nbpibztItEVjNLl7902b6JB3iP0t24fllBccgqz83jJ2M/69qqtW5Fhyvf4n/dReRlBRT57GPc970\nexRVJakVyhF0q7ajcXlw94gnoV/vlh0sKZk6m9BVVoOfCdweQi55CqoclC59AW/ciTdqxz8RSUpK\nhvhEVMt8TEeLSA6ORA1teBbe+OFSzC9+TPX0cVTfcnHLnsevaLILUEpteLrFNtgq0O32kJl5lJ49\n09DpwsnOtrFr1yr69ZuGoij07Vu73WK/fhe2Wpyt4pcbMOMvu4z49LZ9E9CQ4zNS6e08TtH25FoR\njSXXivi1srKyere3WuI9c+ZM3nnnHT777DMCAwNrarr9/f3x8/NDURTuuusunn76aVJTU0lOTuap\np57C39+fq6++urXCaLeU8kpCUq/Fk9iJ0lXzQXPiZkhblZZrnxrIyu1hhAY4+PDxtcRGHGXB4m+p\ndvoWiPEzmbjmgvG1km4A1/DeuIY3IfFtpRpg/ZqdvvMPOn092INH3Yliq8IbHeprO5jYCW+XyFPv\noNNSfc1YVKMePJ5GnUO3/QC6zCMoldUtiFTDE098yC233EhiYiper53pN1/G65t+JtSgw909nqzc\nIjr1S6JqzgyqwkN47bVvueeeG9FoLHi98Oqr7/PWW7ei1Wrp3h26d2+9FoOnVWUl7N7t6/3ds2db\nRyOEEEK0a62WeP/jH/9AURQuvLD2bNzs2bN59NFHAbjvvvuoqqpi5syZlJSUMHjwYBYvXoyfn19r\nhdFu6bbsQ1FV1EC/Wkl3eaWO3z2ewYY9IUQGV/PxE2uIiSjgP4u/paLKDoBBp+fq88cR4t+8OurT\nQbfNt0Kna3CPBkY2k9uDUlyBprwSTaHv3WP1VRc2+Mahcs7NTTqNdrevLMWT2qVJ+y1duo3OnRPo\n02cwBkMEN9/chcTEJMxmXx/zZ0ePJ6i4HOXAATSbMpnZqROfL1qL/9Tr8f/9OBISSrFa+6PVavH3\nh3fffb9J5z8jnngCliyBv/3t1N1Ntm713ZjQowe0YKVSIYQQ4lzQaom395dyiIY89thjPPZYy7pX\ndES6TZlA7dUoi8v1TJs9iK37g4gJt/Pxk2voHFbK298tpqjCl2xqNRqmjRxNdEjTbjpUCst8pRat\n1eXiNyr+fR/2zGy8nUJPy/HRaSne9x6a3CK0e7NRiitwTm79WeDjibe7kYm3opgxm5MIDNTi7x+L\nn5+vvr1Xr9p9wLvOnQtz50JpKdpNm/hm40ZfYjp9Bmg0XH/99a36PE6LnTthxQrYs+fUibfRCFdc\nAb8pLRNCCCHEyU5bO0FRW03iPcCXoBSUGrji0cHsOhxAXFQlnz65hujQCv677Htyigtr9psybBQJ\nUZ2adrIqB4GXPYy3UxgV8+5EjWjGMt52B/oNu3Glp4LFePJ2jabJs8RNpih4O4Xh7XR6lvFWymxo\nc4tQzQa8cfWUsAAlJRW8/voynnjib5hM/owdG9+4kwQFwQUX+L46mrg43/f6FtEZMAA++ujMxCOE\nEEJ0cJJ4nwmqin6jL/F29U8ht8jEFY8OYm+2P8kxFXw4ezWFFTv4dOVGyuyVNbtNHDiU7l2a3jFE\nuzcbTV4xul2HCR51J56YCKr+OA3n6AFwfLGcBgRe9hD6jZmU/+s+0OvQbdiN/eHptcpkOjptpq/9\nnic5toHXRUtU1EDOPz8Yo9Faz7izTGMSbyGEEEI0miTeZ4BSYccbZAWXm+KIOC69dwgHc/3oHl/G\ns7f9j4XrlnOstHbv55G9+5Oektas83l6d6X0x3lYb38ew4ptaArK0N30DEVHPm70MdwZaeg3ZhLw\nf8/WPOa4YhSe7vHNiqk9cg9MpWjX22hKKurcXl5u5+hROxkZl2AwBDNhwtnTYT++k/AAAB00SURB\nVKdRJPEWQgghWtXZM33ZjqkBfpSunE/RljeY9VJfDub60a1LEdPGPMk3Gz+vlXRbjCYmZQxjRM++\nLTqnNzqU8o+foPKR61D1OqqvHtOkem/HJF89tarX4RrSA/u9v/PdGNoBaDOPYHnsDUyvLWxwrBoW\niCe57jaMWVle1q8vxmgMab0VITuS4334JfEWQgghWoXMeDeDUlxOwO+fxHH5SKqnj8PwwyZcw3uh\nWutfdOaNpd34em0UZmMVg3o+SEHZsZptep2OIWm9GJLWE6Pe0DqBarVU3Xk5VTdNqrtOux7uQd0p\n2vYf1AC/Ju/b1jR5xVhe+QxXv2SqZ1zUjCNosVh6MXJkDKNGnYMJ93Fdu8KiRZDQ8qXthRBCCCGJ\nd7MYP1+BfsMe1AA/DN9vwvDdBirm343jyvNPuc+2AwE89oavdGRk/5cIsPqSbkVR6J/UjZG9+mE1\n15+4N5tf89q8qVEhrRwIKIoFVa3G7Xah0zWu3rw+Xq8XrVYPaFBVla+/Xs/E4d1RFQXt9oPs2LyP\nHv2SGn28hQs3UFFh5Y47Jpybs9y/ZjLBhAmn3v7WW2CzwaWXQqcm3gAshBBCnIOk1KQZjB8tA8Ax\ndRTO8Rm1HquLza7l5mf743Rr6dn1G5JiVwOQFhvP7ZMvZ1LGsNOXdLcTO3Zk8dJLywgKGklg4Cge\nf/wL1qzJQ1GaP7uvqipz5iyiuDiewMALCAgYxcaNRQTETERN7QYuNy/99Utcrsa98dBoLEyadCNX\nXXW9JN2NMW8ezJwJe/e2dSRCCCFEh3DOzHhrDuailNvx9OnasuMcykO/fjeqxYhj/CAUlxu/B15D\nv3wrSl5xnbPE97/akwO5VkIDDzG8zwIAxqcPIaPb6Vv1sT3RaKz07DmZoKBStFodWq2OefNexev1\notNpcDqLufHGW3nwwauIjjY36piKYsRiSeOuux4gODgYvd6334svvuwbMHgI7NrNxxMvxh0+Ers9\nk/LyTEwmvW+73QE6DRj0OJ0uqqr0xMZmYDD4n46X4OzjdMK2bb6f+/Vr21iEEEKIDuLcmPH2eAic\n8hBB4/6I9pcVF5vL+MmPADgmDgarGTXYH+eYdBSvF+P/lp80fskCO/mfHyCQQsYN+Rs6nZNh3Xuf\nE0n3tm0HOXxYJSBgKKGhcfTp06dmm16vx2g0otXqMZsjmTv376SmTsHPLx0I4/vvt6Cq6knHLCws\n46WXluDvPwSLJZbg4FP0KB80yPd97Vp0OgP+/j148smv2LjxCACmd5cQGncllqffYenSfXzwwVZJ\nupti505f8p2cDAHtZ0VVIYQQoj07J2a8dRsz0R71LUrjN3sB5R8/0ewVHfXrdgFQMvlCfv9YBmWV\nBu6LncY01mD6aBnVt11aM3Zvth9Vf/+MZblv8I/zLmVXgErvhCQu6Jve8ifVzimKkcJCf7TaIHS6\nhm/OjP+lg4ZeH01pKeze7WTChERcrmxU1fXLMXVERWUwYEAAen0D/bTHj4c33oDhw3/ZV+H55/6O\nYcsmvM/9F90Hn6O4PRDRmSuvvBOttpVuaD1XbNrk+96/f9vGIYQQQnQg50Tibfhm7Ymfl29F/+MW\nXKOa9/F4+fuP4Vq1lys+v4bVeyIAuN57C7qQn/g5dCyWL7tw8fA8rGY305/uxcdlGwCo6ltIYnRv\nLhp83lldP2yzVREU1AWLpQdTpzZvBjk6OprHHnsSVVXxeJJ4993X8ffXMHnyNej1AVx8cSOWJ4+L\ngxtuqPWQ38yZ8M47ACy2WMgLD+Oaydej0Tfv5tNzwtq1vjru7t19N1MeJ4m3EEII0WTnROJdPX08\nalggmqxjOJNjKB/YDZe9EpfbjcvjxuV242cyE+Lf8EfmDreWaxf7ku6IYBs3TFjJ+t1JXGN8H4dT\nB6/DQ//uReewcopztfS0b8ejUfD0SWTaeRegPYtWfqzLww9/wJ/+9DB9+7a8bENRFHQ6ExdffA37\n9+/HYAhs2QFHjIBVq2DiRAaffz7bg4PR9Ozd4jjPano9bNzoKyv5tWnTICQExoxpm7iEEEKIDuic\nSLyLwvxYmGLkcJAeVc2Dz/5b57hOoeH0T0qhR1xinb203R6F/3u2N8u2RGAxlTFm0INUunPongRJ\ncSYO5mSwN2s4WXn9OHIskJH2H9DiJScmmKljJ7Zef+52ScFkSuX119/BbG7cDZKNFRwcTHp6K5Tn\n3Hgj/N//gaIQAAxt+RHPfsdXrzx0CFT1RInWeef5voQQQgjRaGd94l1iq+CtJV9RZq9scGxOUQE5\nRQUs3riOHnEJ9E9KpVNoGACH8vK486W+rN0Rg1Fv4+IRswkOyKnZ16CvplvccrrFLafaYeVQbjrX\n/PwFAJbh/TCcpe0CVVXlyy83MHnytfj5dW3fZTTalvcNP+eEhICfH1RUQGkpnOpmViGEEEI06KxO\nvEsqynnru0WU2SupqAyjpKIzQf6FhAYUYzIo6HU69FodWq2G/JJiPF4vAE63i837M9m8P5PIoBC8\nXpWPl13Btn0D0euqmHzeU4QFHcJiNBEXEUW1y0mVoxq7w0GVwwFGG6nxy1CrnZRq+qG5cGAbvxKn\nptH4odH48dZbXzBgQAppaV1QVQ/ffLOabt1iiI8Pa3D/8vIAPB7/9p10i+ZRFN+s986dvqXjJfEW\nQgghmu2sTbx/nXRXVgXz4Xd/pcrhqxHWaFRiwqqIj64kPspOfFQlqZ3KsDv2kFuymSpXVs0n6vml\nxWzafAWTlm2lMiKJXiPfISOtlIEpo0iNjUdXxyyqp9yG3evB/Ds/3Nr2+RK7XG5WrjzAJZfchk5n\nJiNDpUuXLgQHRwNgsVQQEJCGv38EHk85jz02l6lTh5OYGAyoeL1e9PpwrNY+zJp1Yds+GXF6/Trx\n7tu3raMRQgghOqz2mRXWQ3MgB212Aa4RfU45pqSinDe/W4SrsATVZOD79XdS5Qikc1gViqJytNBM\ndr6RsdvepYd9B3cnvPDLnkOA6zHo3AT5l2I05KLVVDNo236eznqIGfwb2w2PEBWcccpzm//2PuaX\n/oft5Vk4J7ffKuLKSh3bt5dy2WUmFEVh8ODBtbZPmTLlV38KYfr0PxAfH4/JpOBylXHJJdewYMFb\nBAf7ndnAxZn34oswfz7ExLR1JEIIIUSH1uESb93Kbfg/+Drl/30U1/CTO1IcT7rL7ZXc8vZGwo+4\n2JhQha2Tk6+fXUFUqAOHS0P+1gr6XHwHWpcT+6Uj+ck4hOwCM0cLzZRUGDhWEgb4yixeKngcgLDb\nRmMNDq03PtViQlNZhfHjZe028dZoAoiLS2f27IsavU+vXr1qftbrrXz88WcEBrawy4joGJKTT/zs\n9cLYsZCaCn//u6/riRBCCCEapUMl3oXlpRxYt5aJ1U5K//h3PnriMkIDgwgNCCI0IBCzwcinK5dS\nbq/EWO0iJbMIrcfLfmNXXrhjK1GhDgCMei9d0v1w3HYRlnmf8KcND3HzV8/UdGywVWnJKTSTXWCm\n+EAVE/7wLapWg+PShrs4OKechzr7PxiWbEApqUANbj+rIXo8Xt55ZwU33vgn9PqWzVRL0n2O2rcP\nvv8edu+Gl19u62iEEEKIDqXDJN7HSkt4+/uv2R81jsGmTOIOFGFdvo3N3SPrHJ+6qxS9x8MK/2FM\nvNjO+EH5J42pumsqpve+Q79+N4aFq3BePAwAq9lDSqyNlFgbps1fofF4cI5ORw1rONn0RoXiGtEb\nw7ItGBauxDF9fIue968pigmNxsjGjT9jMGjo3j22SftrtX74+ydgMskS36KZZOEcIYQQotk6xGou\neSVFvPXdIjbtSefjdY/zVNRjAIz74gCKVz1pvE6rJWmDr/vCqoTxPH7jzjqPq/pbsN9/NQABNz1T\n5xjFXg1A9dSRjY7XcZlvrP8fX/F9NN8KNm06gMnUh6Cg8/B4klHVblitQ7BY+jF//nK++movOl04\nimJg375cCgvLau2v0VgICBjIjBm3t3qfbXEOkcRbCCGEaLZ2n3jnFhfy9ndf8/O+fny3bhag4dNu\n13LU0ImEvGN0W9+P/kndiY+Mxt9sIdjqT4z1BnrvXgnAiEcSsBhPnfxWXzMWV/8UvOFBdW+/8gK8\noQE4xw1qdMzOi4bg6RyGa0AKtMJKlaoKX321m9JSF4qiMG7cOEaOvBCTKQyLJYb773+SK664hcDA\nQQQFjWbjxiqysvwwm9PQ62P517+WkZVlbfnKj+Lcpqq+VSwBBgxo21iEEEKIDqjdl5rseXw+oeVd\nWVxyD6qq5eaLNvPI9FyW3XYrfZe8z2dbh+IM6sv8u7YQGeIgr8jI1TMSGWTpTnLAMZJGNVBWodNS\n9vWzKOX2OjergX6ULn0R/EyNjlm1WihZ9Q/Qtc77GpMphRdemHzKPtkm04nYtFotM2fOqrV9+HAb\niYlprRKLOEdNngw//gg2m+/PMuMthBBCNFmbzHi/8sorJCQkYDabSU9PZ8WKFacce9GHW7nvyw8Z\nUL6J6yds46mbjmIyeBn3ai+2ffQ6m7uMYPnWcEbNGsG36yKY+UJftru78si1n8LqustHTqLRoAZZ\n695mMuCNrr+TSZ0sRjC0rOODzVZFQYEeP7/kFi1OM3z4cAICpK5btEB1tS/pnjcPPvoIOnVq64iE\nEEKIDueMJ94ffPABd911Fw8//DBbtmxh6NChTJgwgSNHjtQ53uh2syD8elIus/LMLYdrFrZBp2X0\noCKWvbickX0LKCo3cu2cDH76OZzQAAcvzdqKxtyxW51t3pzDF19sQdMK5SpCtEhcnO+7Xg9XXAGy\nSqkQQgjRZGc8o3v++ee54YYbuOmmm+jWrRvz5s0jOjqaf/zjH3WOP2SMY/l1M3l+5oE6f9dHhjj4\n4LG1PHb9TvQ6Xy33i3duJTLEcTqfxmmnKGYmTbqJ++//c1uHIsSJxPvw4baNQwghhOjAzmiNt9Pp\nZNOmTdx33321Hh87diyrVq2qc58Flz3FM/fm1jvBptHAzCkHGDcwn+IKAxlpJa0Z9hmXn19O166D\nMBjaTw9wcY6TxFsIIYRosTM6411YWIjH4yEysnbv7YiICPLy8urc5/a/BzauMYititSCLR0+6Xa7\nPdx331tUVmrbOhQhTjieeB892rZxCCGEEB1Yu+9qcvDg3gbHGHOL6XnTi+hLbBy7KIOCiQOp6Jt4\nBqJrXRqNAYMhhdmzn+LgwYMcPHiwrUPqMDZs2NDWIZzVFL0e7ddf4w4JgbPgtZbrRTSWXCuiseRa\nEQDJycn1bj+jM95hYWFotVry82uvIpmfn090dHSzj+uICsYRHQJAxMJ1RH6yskVxnml2u5M33lgN\n9MbtDmlRBxMhTgfVaMQdFtYqfemFEEKIc9UZnfE2GAwMGDCAxYsXc/nll9c8vmTJEqZOnVrnPklJ\n9b9zOM7z+E0w1beipWHqhY3er+1pMRqTqa7uxYABw9BqpcSkKY7PMKSnp7dxJKIjkOtFNJZcK6Kx\n5FoRv1ZWVlbv9jNeanLPPfdw7bXXkpGRwdChQ/nnP/9JXl4et956a4uO6xrZF8fFw9DuOoxrzMBW\nivb0ycsrITu7ivPPn4bBEMjEiR3ljYIQQgghhGiOM554X3nllRQVFfHUU0+Rm5tLr169WLRoEbGx\nsS07sKJQ8e/7WyfI006L3R7K0aO5GI11L1UvhBBCCCHOLm1yc+Vtt93Gbbfd1hanbgcULJbeDBzY\nmYwMqeUWQgghhDhXyJ1SZ9DPPx/ixReXYTZ3lhsohRBCCCHOMe2+neDZpHfvC4mMVCXpFkIIIYQ4\nB0nifQaoqorZ3A0/vxTCwiTpFkIIIYQ4F0mpSQuoKjT0EpaX25kx45/odHEy0y2EEEIIcQ6TxLuZ\n9u3L5c9//oDAwJGYzT1YtGgHe/ZknzQuLKw7f/vby5hMpjaIUgghhBBCtBdSatIsCj17juaRR8ai\n11vR661ERw8gPDwKs9mM05nP1q2b6dt3GFZrDwIC5P2NEEIIIcS5ThLvJqqudhEc3A+LJR6r9UTp\nyJgxY2p+NpsT+eSTd+nSZQJBQZJ0CyGEEEIIKTVpEputmhtueAmNJqreem2NRsNLL71CYmLXMxid\nEEIIIYRozyTxBg4fPsaKFTsaGKUQHt6fTz/9ArPZfEbiEkIIIYQQZw9JvIGAgBRycixYLL3Q6zuR\nnV3Kvn05Ndtzc4sxGrthsSTi7+/fhpEKIYQQQoiO6pxPvI3GJJKShnP77XdiscQTENCfgoJIMjMV\njMYkNJpA5sz5nMJCnbQDFEIIIYQQzXbO3lyZl1fM/Pk/MH/+gloJtaIojBkztubPqprKu+8OQ6vV\ntkWYQgghhBDiLHHOJt6xsQO4+eZeaDT1T/oriiJJtxBCCCGEaLFzLvH2eLyYzUn4+3enf38pHRFC\nCCGEEGfGOVXj7XC4uP76eahqjNRrCyGEEEKIM+qcSrz9/VNYsOBdAgMD2zoUIYQQQghxjjknSk1U\nVUWvj8Jq7S7LtwshhBBCiDZxTmSh77+/gv/9b2uDN1IKIYQQQghxupwTM95Tp16HVhvR1mEIIYQQ\nQohz2FmfeOv10YSGdpebKYUQQgghRJs6q2svvvxyI4WFfpJ0CyGEEEKINndWJ94QhtUa0tZBCCGE\nEEIIcfaWmuj1nbjuusky2y2EEEIIIdqFs3LG+/DhEiyWVEm6hRBCCCFEu9EqiXdJSQl33HEHaWlp\nWCwWunTpwu23305xcfFJ46699lqCgoIICgpi+vTplJWVtUYINRwOF08//QlVVZ5WPa4QQgghhBAt\n0SqJd05ODjk5OTz33HNs376dd955h+XLl3PVVVfVGnf11VezZcsWvv32W7755hs2bdrEtdde2xoh\n1LBa4/noo/8REBDQqscVQgghhBCiJVqlxrtHjx588sknNX9OTEzkueeeY/LkydhsNqxWK7t27eLb\nb79l5cqVDBo0CIBXX32V8847j8zMTFJSUloUQ1lZJVptEJ07p8lCOUIIIYQQot05bRlqWVkZRqMR\ni8UCwOrVq7FarQwZMqRmzNChQ/Hz82P16tUtPt/ChTtYubIAvd7c4mMJIYQQQgjR2k5LV5PS0lIe\neeQRZsyYUTP7nJeXR3h4eK1xiqIQERFBXl7eKY+1b9/eBs6mwWRKYdiwi1EUhQ0bNrQ0fNEByd+7\naAq5XkRjybUiGkuuFQGQnJxc7/Z6Z7wffvhhNBpNvV/Lly+vtY/NZuOiiy4iNjaWZ599tuXPoB5L\nluxi/34jbneEdDARQgghhBDtWr0z3nfffTfTp0+v9wCxsbE1P9tsNiZOnIhGo+HLL7/EYDDUbIuK\niqKgoKDWvqqqcuzYMaKiok55/KSkut85aDRWvN6+xMV1o3PnzvXGKM5ex2cY0tPT2zgS0RHI9SIa\nS64V0VhyrYhfa6hbX72Jd2hoKKGhoY06UUVFBRMmTEBRFL7++uua2u7jhgwZgs1mY/Xq1TV13qtX\nr6ayspKhQ4c26hwApaU2QkLi8Pfvy9ChloZ3EEIIIYQQoh1olZsrKyoqGDt2LKWlpSxYsICKigry\n8vLIy8vD5XIBkJaWxvjx47nllltYs2YNq1ev5pZbbuGiiy5qsB7m1/7610Xs3q2i10vSLYQQQggh\nOo5Wubly48aNrF27FkVRarUFVBSFpUuXMmLECADee+897rjjDsaNGwfAJZdcwssvv9zIs2gxm9OY\nN29crRIWIYQQQgghOoJWSbxHjRqF1+ttcFxQUBBvv/12k46tqipvvrmM3/3uZkJD4+UmSiGEEEII\n0SG1+5VmdLpgunQZgMUSKUm3EEIIIYTosE5LH+/W5O8/iCuvPK+twxBCCCGEEKJFFFVV1bYO4rca\nasUihBBCCCFEexYYGHjSY+2+1EQIIYQQQoizgSTeQgghhBBCnAHtstRECCGEEEKIs43MeAshhBBC\nCHEGSOIthBBCCCHEGdAuE+9XXnmFhIQEzGYz6enprFixoq1DEm1s7ty5DBw4kMDAQCIiIrj44ovZ\nsWPHSeNmz55N586dsVgsnH/++ezcubMNohXtydy5c9FoNNxxxx21HpdrRQDk5uZy3XXXERERgdls\npkePHixfvrzWGLlWhNvt5sEHHyQxMRGz2UxiYiKPPPIIHo+n1ji5VkRD2l3i/cEHH3DXXXfx8MMP\ns2XLFoYOHcqECRM4cuRIW4cm2tCPP/7IH/7wB1avXs0PP/yATqdj9OjRlJSU1Ix55plneP7553n5\n5ZdZv349ERERjBkzBpvN1oaRi7a0Zs0aXn/9dXr37l1rAS65VgRAaWkpw4YNQ1EUFi1axO7du3n5\n5ZeJiIioGSPXigB4+umnefXVV3nppZfYs2cPL774Iq+88gpz586tGSPXimgUtZ3JyMhQZ8yYUeux\n5ORk9YEHHmijiER7ZLPZVK1Wq3755Zeqqqqq1+tVo6Ki1KeffrpmTFVVlerv76+++uqrbRWmaEOl\npaVq165d1WXLlqmjRo1S77jjDlVV5VoRJzzwwAPq8OHDT7ldrhVx3OTJk9Xrr7++1mPTp09XJ0+e\nrKqqXCui8drVjLfT6WTTpk2MHTu21uNjx45l1apVbRSVaI/Ky8vxer0EBwcDcPDgQfLz82tdOyaT\niREjRsi1c46aMWMGU6dOZeTIkai/at4k14o47rPPPiMjI4Np06YRGRlJv379mD9/fs12uVbEcRMm\nTOCHH35gz549AOzcuZOlS5cyadIkQK4V0Xjtasn4wsJCPB4PkZGRtR6PiIggLy+vjaIS7dGsWbPo\n168fQ4YMAai5Puq6dnJycs54fKJtvf766xw4cID33nsPoFaZiVwr4rgDBw7wyiuvcM899/Dggw+y\nefPmmnsBZs6cKdeKqHH77beTnZ1NWloaOp0Ot9vNww8/zK233grI/yui8dpV4i1EY9xzzz2sWrWK\nFStW1EqoTqUxY8TZY8+ePTz00EOsWLECrVYLgKqqtWa9T0WulXOL1+slIyODOXPmANCnTx/27t3L\n/PnzmTlzZr37yrVybpk3bx4LFizg/fffp0ePHmzevJlZs2YRHx/PjTfeWO++cq2IX2tXpSZhYWFo\ntVry8/NrPZ6fn090dHQbRSXak7vvvpsPPviAH374gfj4+JrHo6KiAOq8do5vE+eG1atXU1hYSI8e\nPdDr9ej1epYvX84rr7yCwWAgLCwMkGtFQKdOnejevXutx1JTU8nKygLk/xVxwpw5c3jwwQe58sor\n6dGjB9dccw333HNPzc2Vcq2IxmpXibfBYGDAgAEsXry41uNLlixh6NChbRSVaC9mzZpVk3SnpKTU\n2paQkEBUVFSta6e6upoVK1bItXOOmTJlCtu3b2fr1q1s3bqVLVu2kJ6ezlVXXcWWLVtITk6Wa0UA\nMGzYMHbv3l3rsczMzJo39fL/ijhOVVU0mtopk0ajqfkkTa4V0Vja2bNnz27rIH4tICCAxx57jE6d\nOmE2m3nqqadYsWIFCxYsIDAwsK3DE21k5syZvPXWW3z00UfExMRgs9mw2WwoioLBYEBRFDweD3/5\ny1/o1q0bHo+He+65h/z8fF577TUMBkNbPwVxhphMpv9v5w5RFovCMACfq4jXBSjcdMWizaLd5Bps\nLsIkGnQNgtniEoxicA0mt2ExfZNGRoZhfhi4/gzPA6d94QtveMM5J7Xb7dfpdDrpeDymsizTfD6X\nFV7Kskzb7TbV6/VUFEU6n89pvV6n5XKZxuOxrPByv9/T4XBIg8EgNRqNdLlc0mq1SrPZLE2nU1nh\n6z76p8of7Pf76Ha70Ww2YzQaxfV6/fRKfFiWZVGr1SLLsrez3W7f5jabTRRFEXmex2Qyidvt9qGN\n+U5+/U7wJ1khIuJ0OsVwOIw8z6Pf78dut/ttRlZ4PB6xWCyi2+1Gq9WKXq8Xq9Uqns/n25ys8DdZ\nxBdeHAEAAP/kW93xBgCA/5XiDQAAFVC8AQCgAoo3AABUQPEGAIAKKN4AAFABxRsAACqgeAMAQAUU\nbwAAqMAPmgLIcozeIMoAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from numpy.random import randn\n",
"from filterpy.kalman import EnsembleKalmanFilter as EnKF\n",
"from filterpy.kalman import KalmanFilter\n",
"from filterpy.common import Q_discrete_white_noise\n",
"import book_plots as bp\n",
"\n",
"np.random.seed(1234)\n",
"\n",
"def hx(x):\n",
" return np.array([x[0]])\n",
"\n",
"def fx(x, dt):\n",
" return np.dot(F, x)\n",
" \n",
"F = np.array([[1., 1.],[0., 1.]])\n",
"\n",
"x = np.array([0., 1.])\n",
"P = np.eye(2) * 100.\n",
"enf = EnKF(x=x, P=P, dim_z=1, dt=1., N=20, hx=hx, fx=fx)\n",
"\n",
"std_noise = 10.\n",
"enf.R *= std_noise**2\n",
"enf.Q = Q_discrete_white_noise(2, 1., .001)\n",
"\n",
"kf = KalmanFilter(dim_x=2, dim_z=1)\n",
"kf.x = np.array([x]).T\n",
"kf.F = F.copy()\n",
"kf.P = P.copy()\n",
"kf.R = enf.R.copy()\n",
"kf.Q = enf.Q.copy()\n",
"kf.H = np.array([[1., 0.]])\n",
"\n",
"measurements = []\n",
"results = []\n",
"ps = []\n",
"kf_results = []\n",
"\n",
"zs = []\n",
"for t in range (0,100):\n",
" # create measurement = t plus white noise\n",
" z = t + randn()*std_noise\n",
" zs.append(z)\n",
"\n",
" enf.predict()\n",
" enf.update(np.asarray([z]))\n",
" \n",
" kf.predict()\n",
" kf.update(np.asarray([[z]]))\n",
"\n",
" # save data\n",
" results.append (enf.x[0])\n",
" kf_results.append (kf.x[0,0])\n",
" measurements.append(z)\n",
" ps.append(3*(enf.P[0,0]**.5))\n",
"\n",
"results = np.asarray(results)\n",
"ps = np.asarray(ps)\n",
"\n",
"plt.plot(results, label='EnKF')\n",
"plt.plot(kf_results, label='KF', c='b', lw=2)\n",
"bp.plot_measurements(measurements)\n",
"plt.plot (results - ps, c='k',linestyle=':', lw=1, label='1$\\sigma$')\n",
"plt.plot(results + ps, c='k', linestyle=':', lw=1)\n",
"plt.fill_between(range(100), results - ps, results + ps, facecolor='y', alpha=.3)\n",
"plt.legend(loc='best');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can be a bit difficult to see, but the KF and EnKF start off slightly different, but soon converge to producing nearly the same values. The EnKF is a suboptimal filter, so it will not produce the optimal solution that the KF produces. However, I deliberately chose $N$ to be quite small (20) to guarantee that the EnKF output is quite suboptimal. If I chose a more reasonable number such as 2000 you would be unable to see the difference between the two filter outputs on this graph."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Outstanding Questions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"All of this should be considered as *my* questions, not lingering questions in the literature. However, I am copying equations directly from well known sources in the field, and they do not address the discrepancies.\n",
"\n",
"First, in Brown [2] we have all sums multiplied by $\\frac{1}{N}$, as in \n",
"\n",
"$$ \\hat{x} = \\frac{1}{N}\\sum_{i=1}^N\\chi_k^{(i)}$$\n",
"\n",
"The same equation in Crassidis [3] reads (I'll use the same notation as in Brown, although Crassidis' is different)\n",
"\n",
"$$ \\hat{x} = \\frac{1}{N-1}\\sum_{i=1}^N\\chi_k^{(i)}$$\n",
"\n",
"The same is true in both sources for the sums in the computation for the covariances. Crassidis, in the context of talking about the filter's covariance, states that $N-1$ is used to ensure an unbiased estimate. Given the following standard equations for the mean and standard deviation (p.2 of Crassidis), this makes sense for the covariance.\n",
"\n",
"$$\n",
"\\begin{aligned}\n",
"\\mu &= \\frac{1}{N}\\sum_{i=1}^N[\\tilde{z}(t_i) - \\hat{z}(t_i)] \\\\\n",
" \\sigma^2 &= \\frac{1}{N-1}\\sum_{i=1}^N\\{[\\tilde{z}(t_i) - \\hat{z}(t_i)] - \\mu\\}^2\n",
"\\end{aligned}\n",
"$$\n",
"\n",
"However, I see no justification or reason to use $N-1$ to compute the mean. If I use $N-1$ in the filter for the mean the filter does not converge and the state essentially follows the measurements without any filtering. However, I do see a reason to use it for the covariance as in Crassidis, in contrast to Brown. Again, I support my decision empirically - $N-1$ works in the implementation of the filter, $N$ does not.\n",
"\n",
"My second question relates to the use of the $\\mathbf{R}$ matrix. In Brown $\\mathbf{R}$ is added to $\\mathbf{P}_{zz}$ whereas it isn't in Crassidis and other sources. I have read on the web notes by other implementers that adding R helps the filter, and it certainly seems reasonable and necessary to me, so this is what I do. \n",
"\n",
"My third question relates to the computation of the covariance $\\mathbf{P}$. Again, we have different equations in Crassidis and Brown. I have chosen the implementation given in Brown as it seems to give me the behavior that I expect (convergence of $\\mathbf{P}$ over time) and it closely compares to the form in the linear KF. In contrast I find the equations in Crassidis do not seem to converge much.\n",
"\n",
"My fourth question relates to the state estimate update. In Brown we have \n",
"\n",
"$$\\boldsymbol\\chi = \\boldsymbol\\chi + \\mathbf{K}[\\mathbf{z} -\\mathbf{z}_{mean} + \\mathbf{v}_R]$$ \n",
"\n",
"whereas in Crassidis we have\n",
"\n",
"$$\\boldsymbol\\chi = \\boldsymbol\\chi + \\mathbf{K}[\\mathbf{z} -\\boldsymbol\\chi_h + \\mathbf{v}_R]$$ \n",
"\n",
"To me the Crassidis equation seems logical, and it produces a filter that performs like the linear KF for linear problems, so that is the formulation that I have chosen. \n",
"\n",
"I am not comfortable saying either book is wrong; it is quite possible that I missed some point that makes each set of equations work. I can say that when I implemented them as written I did not get a filter that worked. I define \"work\" as performs essentially the same as the linear KF for linear problems. Between reading implementation notes on the web and reasoning about various issues I have chosen the implementation in this chapter, which does in fact seem to work correctly. I have yet to explore the significant amount of original literature that will likely definitively explain the discrepancies. I would like to leave this here in some form even if I do find an explanation that reconciles the various differences, as if I got confused by these books than probably others will as well."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## References"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"- [1] Mackenzie, Dana. *Ensemble Kalman Filters Bring Weather Models Up to Date* Siam News, Volume 36, Number 8, October 2003. http://www.siam.org/pdf/news/362.pdf\n",
"\n",
"- [2] Brown, Robert Grover, and Patrick Y.C. Hwang. *Introduction to Random Signals and Applied Kalman Filtering, With MATLABĀ® excercises and solutions.* Wiley, 2012.\n",
"\n",
"- [3] Crassidis, John L., and John L. Junkins. *Optimal estimation of dynamic systems*. CRC press, 2011."
]
}
],
"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.4.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}