{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Kalman filtering 1.\n",
"\n",
"## Brief history\n",
"\n",
"* Rudolf Emil Kalman (1930-[2016](https://en.wikipedia.org/wiki/Rudolf_E._K%C3%A1lm%C3%A1n)) American electrical engineer and mathematician, born in Hungary, publishes his theory (1960) \n",
"* First application: unmanned Lunar mission (1963)\n",
"* Applications in navigation and engineering\n",
"* Sensor fusion (radar, laser scanner, camera, INS,…)\n",
"* Integral part of GNSS receivers, smartphones, computers, games, etc.\n",
"* Forecasting stock prices, time series analysis, approximation of functions\n",
"\n",
"The original [publication](https://www.cs.unc.edu/~welch/kalman/media/pdf/Kalman1960.pdf) received more than 30 thousand citations. Let's have a look at it:\n",
"\n",
"\n",
"## Do you understand?\n",
"\n",
"Kalman filter filter yields the exact conditional probability estimate in the special case that all errors are Gaussian. The underlying model is similar to a hidden Markov model except that the state space of the latent variables is continuous and all latent and observed variables have Gaussian distributions.\n",
"\n",
"The above definition is correct but does not make it understand for most of us **what** Kalman filter really is. To understand it, we will consider its properties one after the other:\n",
"\n",
"* recursive\n",
"* optimal\n",
"* makes estimation of indirect (hidden) state\n",
"\n",
"Finally, Kalman filter equations will be derived which may be easier to grasp from geodetic perspective and discuss a simple navigational example to show how Kalman filter works."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recursive estimation\n",
"\n",
"A very simple example is measurement of the state $x$ (which is a scalar variable) of a system at the epochs $t_1, t_2, ... , t_n$. \n",
"\n",
"Our measrurements are $x_1, x_2, ... , x_n$.\n",
"\n",
"The best state estimate (assuming Gaussian distribution) in this case is the well-known arithmetic mean or *average*:\n",
"\n",
"$$ \\mu_n = \\frac{1}{n} \\sum_{i=1}^n x_i.$$\n",
"\n",
"Suppose we have a new measurement $x_{n+1}$, and our goal is to calculate a new state estimate by using this new measurement. This is the average based on $n+1$ measurements:\n",
"\n",
"$$ \\mu_{n+1} = \\frac{1}{n+1} \\sum_{i=1}^{n+1} x_i.$$\n",
"\n",
"A much more efficient calculation is to \"update\" our previous average $\\mu_n$ by using our new measurement:\n",
"\n",
"$$ \\mu_{n+1} = \\frac{n}{n+1} \\left( \\frac{1}{n} \\sum_{i=1}^n x_i + \\frac{1}{n} x_{n+1} \\right) = \\frac{n}{n+1}\\mu_n + \\frac{1}{n+1} x_{n+1}.$$\n",
"\n",
"Our calculator's `STAT` function works similarly: it \"updates\" the average at each step by using the input data at each step. This is a *recursive estimation*: estimation is updated after getting more data. In our example recursive state estimation can be written in a slightly different form:\n",
"\n",
"$$\\mu_{n+1} = \\frac{n}{n+1}\\mu_n + \\frac{1}{n} x_{n+1} = \\mu_n + K(x_{n+1} - \\mu_n),$$\n",
"\n",
"where \n",
"\n",
"$$K = \\frac{1}{n+1}$$\n",
"\n",
"denotes *gain*, which tells us the correction to the previous average $\\mu_n$, and\n",
"\n",
"$$x_{n+1} - \\mu_n,$$\n",
"\n",
"is *innovation* (i.e. new information), since it decides whether our previous estimate changes at all. For when the new measurement $x_{n+1}$ and the previous average $\\mu_n$ are the same, our estimate does not change at all.\n",
"\n",
"Let us have a numerical example:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEICAYAAAB25L6yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABCAklEQVR4nO2deXxU1dn4v0+SCQlbwhIQAkhcEdkNClKVRUVFFLX21R+18qr1VastXaz42iq1VrHSurS1vNa9Vaq4IFVbFHHFNWxu4IKsAdnDEhKynd8f906YDLPcWe69Zybn+/nkk5k7d+597plznvOc5zznOaKUwmAwGAz6kuO3AAaDwWCIjVHUBoPBoDlGURsMBoPmGEVtMBgMmmMUtcFgMGiOUdQGg8GgOUZRGxwjIieJyBd+yxELEdkrIod5dK9ZIvJrL+4Vdt+rRWSz/axdvL6/wXvExFHrg4isAboDjcBe4D/AtUqpvX7KpSsi8gbwD6XUgx7cawpwhVLqO27fK44cAWA3MEIptdxPWQzeYSxq/ZiolGoPDAGGAjd6eXMRyfXyfoaE6Q4UAJ/5LYjBO4yi1hSl1LfAfCyFDYCIjBCRd0WkSkSWi8jokM86i8gjIrJRRHaKyFz7+BQReSf02iKiROQI+/WjIvJXEXlZRKqBMSJyloh8LiJ7RKRSRH5hnztaRDbYr6eJyDNh171XRO6zXxeJyEMissm+xm3ROgERybGvt0pEtovI0yLS2f6sQET+YR+vEpGPRKS7iPwOOAn4s+0C+HOUZ7tfRP5tn7NIRA4RkXvsMlopIkND5AjKsMd+/vPs48cAs4CR9nWqQq5/W8j3fygiX4vIDhGZJyI9w8r8KhH5yr73X0REopRHG1vGjfbfPfaxo4Cg66lKRBZG+G5f+17/LSLr7XtdJSLDReRjuwz/HPady0RkhX3ufBE5NOw3XS8iu0VksYicFPLZdPu3etwus89EpDzSMxlSRCll/jT5A9YAp9qvewGfAPfa70uB7cBZWB3safb7Evvzl4CngE5AADjFPj4FeCfsPgo4wn79KLALGGVftwDYBJxkf94JGGa/Hg1ssF8fCuwDOtrvc+3vjbDfzwX+D2gHdAM+BP4nynNPBd63n7mN/b3Z9mf/A/wLaGvf47iQe76B5Y6I9Wzb7O8UAAuB1cAP7GvdBrwe8t0LgZ52OfwXUA30iFGOjwK32a/H2vcaZj/Dn4C3wuR6ESgG+gBbgTOilMetdnl0A0qAd4Hf2p/1ta+VF+W7wc9n2c98OlBr/x7dsOrRFg7Uj0nA18AxQB7wK+DdkOt9H+hif/Zz4FugwP5sun3ts+zyvAN43+92lI1/vgtg/kJ+DEtR7wX22I3tNaDY/uwG4O9h588HLgV6AE1ApwjXjKRgwpXZ42Gfr8NSkB3Djo/GVtT2+3eAH9ivTwNW2a+7A/uBwpBzLw5VimHXXQGMC3nfA6i3lcNltqIaFOF7bxBfUf8t5LPrgBUh7wcCVTF+j2XAuTHK8VEOKOqHgN+HfNbefoa+IXJ9J+Tzp4FpUe67Cjgr5P14YI39ui/OFHVpyLHtwH+FvH8WmGq//jdwechnOVgd8KFRrr8TGGy/ng4sCPmsP1DjdzvKxj/j+tCPSUqpDlhKsR/Q1T5+KHChPXStsoff38FSar2BHUqpnUnec33Y+wuwrKS1IvKmiIyM8r0nsRQwwP+z3wdlDQCbQmT9PyyLLhKHAs+HnLsCa0K1O/B3rA7pn7Yb4PdiTag5ZXPI65oI79sH34jID0RkWYgcAzhQ/vHoCawNvlHWBPB2LAs2yLchr/eF3jvWtezXPaOcGw2nz30ocG/IM+8AJCi3iPzcdovssj8vomWZhD9TgYjkJSirIQ5GUWuKUupNLIttpn1oPZZFXRzy104pNcP+rLOIFEe4VDWW2wAAETkk0u3C7v2RUupcLMU6F8v6i8QcYLSI9ALO44CiXo9lUXcNkbWjUurYKNdZD5wZ9mwFSqlKpVS9Uuo3Sqn+wInA2Viui4PkTgXbL/s34Fqgi1KqGPgUS2k5uddGLKUXvF47LJdBZRLitLgWlqtkYxLXccJ6LJdUaNkXKqXetf3RNwDfwxqtFWO5ySL61g3uYRS13twDnCYiQ4B/ABNFZLyI5NqTbKNFpJdSahPWEPZ+EekkIgEROdm+xnLgWBEZIiIFWMPVqIhIvohMFpEipVQ9VihYY6RzlVJbsdwPjwCrlVIr7OObgFeAP4hIR7EmCw8XkVOi3HYW8LvgJJaIlIjIufbrMSIyUKyJyN1Y7oSgPJuBdMVMt8NSxlvt+/43lkUdZDPQS0Tyo3z/SeC/7XJuA9wOfKCUWpOELLOBX9nl0BW4Gev3d4NZwI0iciw0TwJfaH/WAWjAKpM8EbkZ6OiSHIYYGEWtMbYifBz4tVJqPXAu8L9YDWc9cD0HfsNLsJTYSqzJoqn2Nb7EmpxaAHyF5VeOxyXAGhHZDVyFNaEUjSeBUzlgTQf5AZAPfI7l13wGy00TiXuBecArIrIHayLtBPuzQ+zv7sZyibzJAaV1L/BdO1rhPgfPFRWl1OfAH4D3sJTyQGBRyCkLsULivhWRbRG+/xrwayz/7ybgcOCiJMW5DagAPsaaUF5iH0s7SqnngTuxXEu7sUYRZ9ofz8cyAL7Ecr/UcrCbzOABZsGLwWAwaI6xqA0Gg0FzjKI2GAwGzTGK2mAwGDTHKGqDwWDQHFcC07t27ar69u3rxqUNBoMhK1m8ePE2pVRJpM9cUdR9+/aloqLCjUsbDAZDViIia6N9ZlwfBoPBoDlGURsMBoPmGEVtMBgMmmMUtcFgMGiOUdQGg8GgORmVN3bu0krumv8FG6tq6FlcyPXjj2bS0NL4XzQYDIYMJmMU9dylldz43CfU1FsZLiurarjxuU8AjLI2GAxZTca4Pu6a/0Wzkg5SU9/IXfO/iPINg8FgyA4yRlFvrKpJ6LjBYDBkCxmjqHsWFyZ03GAwGLKFjFHU148/msJAbotjhYFcrh9/tE8SGQwGgzc4UtQi8lMR+UxEPhWR2fbee54yaWgpd5w/kNLiQmuL5OJC7jh/oJlINBgMWU/cqA8RKQV+DPRXStWIyNNYe8E96rJsBzFpaKlRzAaDodXh1PWRBxSKSB7QFve2rjcYDAZDGHEVtVKqEpgJrMPaXXmXUuqV8PNE5EoRqRCRiq1bt6ZfUoPBYGilxFXUItIJOBcoA3oC7UTk++HnKaUeUEqVK6XKS0oi5r42GAwGQxI4cX2cCqxWSm1VStUDzwEnuiuWwWAwGII4UdTrgBEi0lZEBBgHrHBXLIPBYDAEceKj/gB4BlgCfGJ/5wGX5TIYDAaDjaOkTEqpW4BbXJbFYDAYDBHImOx5BkNrwqT0NYRiFLXBoBkmpa8hnIzJ9WEwtBZMSl9DOEZRGwyaYVL6GsIxrg+DIU2ky6/cs7iQyghK2aT0bb0Yi9pgSANBv3JlVQ2KA37luUsrE76WSelrCMcoaoMhDaTTr2xS+hrCMa4PgyENpNuvbFL6GkIxFrXBkAbMVnEGNzGK2mBIA8avbHAT4/owGNJA0E1hVhMa3MAoaoMhTRi/ssEtjKI2pIzJS2EwuItR1IaUMHkpDAb3MZOJhpQweSkMBvcxFnUrI91uCpOXwmBwH6OoWxG/mvsJT7y/DmW/T4ebwuSlMBjcx7g+Wglzl1a2UNJBUnVTmPhhg8F9jEXdSrhr/hcHKekgqbgpTPxw9mOievzHKOpWQixlnKqbwsQPZy8mqkcPjOujlRBNGQsYN4UhKiaqRw/iKmoROVpEloX87RaRqR7IZkgjkXzJAkwe0cdYRoaomKgePYjr+lBKfQEMARCRXKASeN5dsQzpxviSMw8dfMMmqkcPEvVRjwNWKaXWuiGMwV1aiy9ZBwWXKrr4hq8ff3QLOcBE9fhBoj7qi4DZkT4QkStFpEJEKrZu3Zq6ZAZDEqRzSyw/0cU3bHab0QNRKlrQVtiJIvnARuBYpdTmWOeWl5erioqKNIhnMCTGqBkLIw7VS4sLWTRtrA8SJUfZtJcihlMKsHrGBK/FMXiAiCxWSpVH+iwR18eZwJJ4Stpg8JNsmfxK1TecDe4fwwEScX1cTBS3h8GgC9myJVYqKz6zxf1jOIAjRS0ibYHTgOfcFcdgSI1sWdKeim9YF/+2IX04cn0opfYBXVyWxWBImWwKQ0w2Sidb3D+GA5gl5Iaso7WEIUYjFf+28W3riVHUhlZPtimnZGOf0x27nW3l6icm14ehVZONE2/J+rfT6dvOxnL1E2NRG1o1sZRTNMWWCZZiMu6fdPq2kylXQ3SMovaYTGjkrYlElZMuS7vdIJ15PVrbhKbb7dq4PjzEDAf1I9G462wOfUtnaGMi5Tp3aSWjZiykbNpLjJqxMOPagxft2ihqD8nmRp6pJKqcstlSTGdeD6flmg3Gixft2rg+bLxwSSQzzM4GN4nOz5Fo3LWuaT/TVcbpCm10Wq7Z4Mv2ovPWXlF70ci98jsm0sizxReaCc+RiHLSMe2nrmXspFz9HqGkQ7940Xlr7frwaljklUsikWG2zm6SRHyKOj9HMuiY9lOXMk7G1+xnbpZ06Rcv0hZobVF7NSzyqldPZJidiExeuhYStd78tpjcQLeVjzqUcbJWfbpGKMm0gXTpFy/SFmitqL2qgF76HZ02cqcyeT3sTbRy6+rTzSZilbFXnXiySi8dSi7ZNpBO/eJ2562168OrYZGOGdecyuT1sDfRyq1j2WYb0cp4TL8SzyIqUlF6k4aWsmjaWFbPmMCiaWMTVnjJtoFMSomrtaL2qpHr6Hd0KpPXw95olVhBRL+kjmWbbUQr49dXbvWsE/dT6SXbBjLJiHC8FVcipHMrLp1Du3TA662nwoeZ4RQGcn1RxKaeHIyX23lFqhde1YVU2kBovSkqDCACVfvqfalD6dqKyxfc9v1kegP3Olws1KcYqXH4EQOra3ia33g99wL+5AFPpQ0E9YvudUh7i9pN/LQC0olfnY0uG7Bmy4a26SZb6rcTUm0DOtShjLao3SQbVkWBf+FiukR0eOGnz8SRlxMrNxOfKxLR2oDT59MhxDEWrVpR6/7j6I4uq/Tc7jB0HxbHIlYnnsnP5YREnk8XoyMa2kZ9eJFRK5PCc3REl4gOt2fvdVn5l26y9bmCJPJ8ukeAaGlRe9XT62IRZjI6rNJzeyIrW0de2fpcQRJ5vlh1SAf3kJaK2ivfcTbtWN3acbPD0H1YnCzZ+lxBEn2+SHVIF/eQI9eHiBSLyDMislJEVojISDeF8rKnT3VVlCE2mZ4UHvQfFidLtj5XkHQ8ny7uIacW9b3Af5RS3xWRfKCtizJlfU/fWtDFGkmVbB15ZeNzhbspLjiulNdXbk36+XRxD8VV1CLSETgZmAKglKoD6twUaky/Ep54f12LGN1s6ukziVT8c9kS/gjuulb89IHqMMeQLiIZBs8urkxpglsXo9GJ6+MwYCvwiIgsFZEHRaRd+EkicqWIVIhIxdatW5MWaO7SSp5dXNlCSQtwwXHZU6EyhVTz9epijehMNmxFpQtuuCl0cQ85UdR5wDDgr0qpoUA1MC38JKXUA0qpcqVUeUlJSdICRSpsBby+Mnnlb0iOVCt+ogmcWiO6+ECzATcMA11CUJ34qDcAG5RSH9jvnyGCok4XxgpriZ/D4lR+i7lLK6ne3xD180z1V6eb1rqPphtEc1MUFQZSuq4O7qG4FrVS6ltgvYgEbf1xwOduCZTMIpREIwsyJRLB72FxsguCgnJX1dTHPM9YjomVsd/1QXeuH380gRw56Hh1XUPGl5HTlYnXAU+IyMfAEOB2twRK1CeUaOWNdb5uCtzvYXGy/rlIckejtY6UgkQqY8Gql+F10O/6oDuThpbSvuBgJ0F9o8r4MnIUnqeUWgZEzOqUbhINGUo0siDa+dPnfcb+hiatQsn8dgMlG76ViHytPeQyPG2sQPNEengd9Ls+ZAJV+yKP4jK9jLRcmZiITyjRyhvteKRhut+hZDqEBiXjn4smd6gSAhNyGSRYxpFSbYbWQR3qg05E8tdnaxlpm5TJKYn6URP9wfzsiXUJDXJK0HUUtAxDKQzkMnlEH99nz3UmnnGRafXBTaK5MMf0K/GsjLx0lWppUSdCoomVop1fEMhhZ4Rhk589cSatHAtfbKA4YEGXaiy3TsSzBjOpPrhNNBfm6yu3csf5A10vI69X3Wa8ok608kY7H9Ayk54OoUFOiBb/rusuKzqGuTkxOjKlPrhNrNGHF2Xk9arbjFfUkHjljXV+JAU+asZCrRq0jmTSRJcb1lA6FL+xmJ3jty/a6/qeFYo6XYQr8GxJKuQFfjecRIgV+ZPM75rOemIsZmf4nUve6/qecZOJXjrwTdyqczJpoitW5E8y9cmveuJ2W9BtXUEofi/t9rq+Z5RFPXdpJc8990++2/QZu3Lb8UTVqa5auJk0nPebTBq2R7OGgKR8jH7UE7dHe06v31oz/3ld3zNKUd81/wsekMc4NrAWgG9VZ/5Tf7xrDvxMGs7rQKYM268ffzRTn1oW8bNklKsf9cTtySwn12/trkEv63tGuT42VtVQKtuY3TCGfaoNJ+SsaD7uBl4Ob3QeZmYbk4aW0qlt5EQ9yShXP9w+blvxTq5vXIPekVGK+rAioViqWae6s6TpCEbYitoty8UrP5hJtuM9t0w8Nm3K1Q9/aboWdKVyfeMa9I6Mcn3cOKoDLISNqjMfNB3DT/Oe5ZBADdePH+LaPbMxJrM1EM93mm4fo9dun3RGPUQqKyfXN65B78goRX1qqbVysL59Tz7c04UcUdwzch8jfIx9TQfGMkkvTn2nmeJTj0S6OppoZXXH+QPjrvDzO0Qu9Bl0aMduklGKml2WK+D+qydC+0Ngxu8ZkbPS+kwpaKyDvDaOLqXTRIixTNJLaxmhpKOjiVVWi6aNjXl9HSJ9dGrHbpJZinp3JSDQoSfk5UOfEfDxUzDgfPj3DVBfA1e+YX0WB7caczK9u9uWSWuwOEIxIxTnpFpWfo9KWkunrI2idqRMdm2A9t0OKOKzZsIjZ8KD40ByQTVCxUMw4uq493OjMSfbu7tpmbQWiyOUWFsyeZUOILQ+FxUGELFyJSdzXzc7Wi9Hc248R2vplLVQ1I6Vye5K6BjyvuQo+MFcePmXcPIv4N0/wRszYOD3oF2XmPdMtII6qWSp9O5OLJNkKnprsThCiTRCCeQI1XUNzXnH3eywwutzaK7zRO/rdkfrlZ/ZredoLW5DLcLzgsrkvJy3GSSrgCjxmLsqoSjsRz1kIFz2bzhiHIz/HdTthcfOhqp1LU4Lj1NOJG+t0+27oq12S0fvHkmGqU8tY+itr8QM40unxaFLrHc8OSKFy7UvyKO+UbU4z62Y33hbkSVyX7djlb0KLXTrOTIpdUEqaGFRB5XGbYGH+UfjqXzccHiL43OXVnLXf1Yyv3YtL+88kvyllZErUvdjYfIz8PSl8MhZcM170KZDxN782cWVXHBcKa+v3BrXQnW6fVc00tG7R2v8O/fVN1smwfPc2PFCFxdKshEdZdNeing9N4bITq65sarG0QjJi6G9F35mt55DhwlNL9BCUQeVSR0B2lDf4niwYQbqd9O+oJYva4t4IpaCOHwMfP8ZeOh0WPAbmDAzZpLxg3IlNzXCtq+g5GgQa5+SRLbvCiddvXusCh1rz8cLjivl2cWVKQ9tdXGhJCuHl0PkWLlEghQVBhx1OMnIrePksZvl7/eEphdo4foIDl/qyCPfVtRBZRJsmD1kOwCbVJf4Q6bex8MJ/wMfPQir33bem1cuhr+NhftPgCf/y5q8JLnKlO5hZDwZqmrqIyqwF5dvSsvQVpdJm2Tl8HKIHOle4fcVwZErIFG5vVjlmowLTDcXhS5uPKdoYVEHlUbDvDYUNNW32Lrpp08t477Anzgj50MANiprkjCughj7a/h6ATx9CSd0/C3v7+580CktlN+2r+Cxc6BNBxh5LVQ8DH85AUZczfPt3ya/5lPaUMe/GkfydOMp7Msr5tu8UnbUHOyOcGNXk0iTPk4IWv2pyqPLpE2ycng5RA6/V6Soj586TAqVqNxuj3x0jGxKFF3ceIngSFGLyBpgD9AINCilytMtyKShpbCoiPO7l3D+hQeUSs+iAsbVLuETdRgvNoxgubL813EVRJv2MHkOPHgafwvMZGTgdvaGeCpa9OZ1+yy/dm4+XPGaNWF5/A/hX1PhrbvoVtSH1WUT+HT9Ns7hLb6X9yYAG5q68lGgHwXsZ2nTEbzUOIIdgUNcsRKCFWj6vM8OcrnE2vMRkkvdGY7fq9CCw/ngxrnJ7Gbu5RA53r2CzxJOpHodruSCVnek67s98nE7sskLnGYG1KFTCZKIRT1GKbXNNUnAWlXYsL/FoV+N7kK7+ft5oeFEHm8cDySgIDofBmffTYenL+GhEVv52Sd9Di54peCFH8GWz62JyGBUSae+cMnzsHMNFPehLCeXMuDf7y1n7ssv07FpB2fnvM/xOSuoUwHODHzE/wZmU922lHYLm+A/++DQE+G4S+GoM5r93akQrOiRKhGQ1tSdke4NSVpETU2wdhGseo3NX37E+m27WFnXjafaTebyM0Y4CkvMto1zE+n4ErEA3R756OICS4V4z6Cjxa2F66OZ3IMV9Zk9qgHY07YvsofEe7d+E6C4DydsmcOi8Jn/vVvgrZnw2XNw6m/gyFNbfi4CnctaHLrtze1U1g8CYE7j6Objx3es4ukxu2i3/gPIbwc5ubDqdZh9EfQ9CU7/LfQcap28YzU01EJJv6QUeDTL5Df/+szZTuorX4LKJVan1H8StD3YLZTIfVuwbwfU7rImZTd8ZCnob96AXetpkjy2N/UClc+FuW9yTu27PPf8WF7bfSXjThkd9ZKZtnGuU9rk5TQ/V6e2AW6ZeGxCUUeRrNi0jHyaGq10DAG73uzdykevz2Xh8lWcnZPLcnU4Vaod+TSSRwMBaUDI5YXFazl3QBcrNDa/rdW+Nn9mtYW8QggUQF4BdDgEuvW3jvtAvM5Ml4nzUJwqagW8IiIK+D+l1APhJ4jIlcCVAH369ElSmjZWBQll+9cA3H3N+dzdqW/i18zJheE/hFd/Dc9dCZJjKeXls+G1W63VjMdNgVE/cXS5aL3xR7uLYeRkGHnNgYON9bD4UXjjDnhgNAy6CEqHwfyboKkeuh5tLX8/9jwryqTO6pTIbxdfEKUOUvK3TDyWG5/7mJz6anJQ7KFty0basB/+c6O1ejPIK7+GoZfA0WdA6XGWjz5RlIIlj8Fbf4BdLePXKSi2RhbjbubUlzrwza4mAPrKJq7Pe4r/l/Mfcl7/NwR+CyN/FLHjilbm8SIrksGLIW+4xQZQW98U9Xzr+RUF1NGRfXSQfexW7aiqqrWUYrsSy8jZ/jWTSvZQPEb4+/uVFO9dxcmFqxjT9hs6vpPD5ve688W2OrbX57M3t4gG8mjfsJN+gc30bVdPh8I2UFAEW1bA/l3Qvjs0NcC+7QwHhgPEyM5QNy8X9a8mBBX9pCD5HaBXOfQYZMm/c43VwfccCp0Pt5R5hx7QrmtyCr2hDlST1Tk0NR3IAyRyUGcmNFEQyOOXpx0G21cR2LWaYqw2WE0h9eRRSC3bq+pi14+6ati9Cboekbi8cRCl4heqiPRUSm0UkW7Aq8B1Sqm3op1fXl6uKioqEpfm7+fB/j1wxYIDx169Gd7/K9z0bfI9cM1OuHeInbhpv+WL3r8b+p8Lo2+Ebsc4vtSoGQsjKoiY1l3tLnjnbnjvfuv+ZSfDMefAZ3MtixMFHXvBnk1WR9L3O3DmndD1KPh8LlQ8AttXWdZ9/T7Y8y1Ub4UjToPRN1iN672/wCdzUDW7EJpoUsLynGOQgd9lyNjvQfUWePFnsGmZ1SmNuQm2fQlv/xFWvnigg+xUBocMgO4DLTkPHXnw8zQ1wuo3Leu5bi8sedyKmOk9whrBtCuxGknPIVByDORYwUVl0146qAl3Yjd3BB7ijNyPoG1XKxa+TQc48nQYMhly86KWuQCTR/Q5OBZ+SE/rOavWWbJ26mtZcIGCsOdogr3fwsalsOJF1u6sZfaa9lQ35rCyqQ8fqaMpDOQlFiWjlNVBx8g346gO1e6Cj5+GFfPYuXoZ7VQ1+RJlIllyINAO6vYc/FlBEfQ+gY17FTs2riJXNdCeGjrJHvJoZDft+LqpJ7tzOjCoZwd65tdY9a5jT6haC7n5/HVZHS9X92OLKqaz7GFwzirasp868qgnj3qVS540UibfEihoxxWTTrfmfQqL4ZBBlnwNtVYunoZaSymv/wDWfQBbV1idQX4H63ffszHs2XKt+lTYyXqWtp2h+FCrLXToAfu2W+1mzyZLSe751rrGvu1YeYEOsdp/Qy3kBKx7FBSxnSLWV9VR0LiHspxv7bDg8NkPi72qgPZSC8BWVURdiH0rCMVtA7TNw7p3++7wi+QW8YjI4mjzf44UddjFpgN7lVIzo52TtKJ+8iJrmfhVbx849s/JlkK59qPErxdK3T6rR93yuWVZHzYaTv9dsxJxSiRrqDCQ66wxV62Db96EwRdBrr3DyO5NsGIerHnHsqob9sOyJy1Fd8hASyEW94HeJ0DVesva7tDDUjrL/2kpSrAq9IALoNOh0KajdfyzubAtpNK0KYLzZkG/s1rKtX+v1WFs+hg2fwLffmK5Z1DWaKNtF2sIW1BkVfiNy6xGHKTLkTDqxzDk+5CTE9XqiKqgigpYNGEHfPM6bP3CamQ7V1sKtrScL2uLmbtiD7k0ki8N5FNPPg1sVx35XPWlhnyOlA0cn7OSJslndNEmOu75uuVNcvOhqLdVLnkFkJNnhV822q62wk7srGmkE7ubv7JBdWW76siWvJ6cNuFCKzIov52VomDLZ1Zn2eVIKOho/W5bPocPH4StK+2OodAasXXoYV1QNUGXw7n+zf2sbepOjigCWIrzMNlEkVRz5UllVj1ZtdCStetRrG43mFdW17OjsZDdtKVaFVKSt4/zBnRiwOF9rBW7NTug57ADnaRqsupNt/6QkxO17Fv8DlGMjUgdbDQEWD1jgsOzsTrL/busupmTA3s2W79LUPnu+Rb2brY6rtoqqN5uKfr66pZ3bd/NKucOPSzl3LGnVQY711rpJAqKLYt3/x5LcVdvtT7Pbwddj4RAW+tSxYeyeN0u5i9eSX1jEx2ooVj2UpVTTJtc6Fy/mVya7LtapdK2TR4TBva02l5JPzhmYlIuzZQUtYi0A3KUUnvs168Ctyql/hPtO0kr6qcvtYZd13544Nj9I60Ge/HsiF/xY3bW9Xvu+AYen2T5+E7/LZRfHrlD2buFxa89w8JP1jC/+ghqio5oKYtSsPlTWLMI2pdAnxOhYw9nMuzfA2/+Ht69z+oESvpZFltegaXwhk6GbscCqoWvPVZHBjjr5JSCFf+CpX+HLSutBttk+d4blVBHgHry6Cj7Woi8QXVFKWF3bjHHTvyxPScgsGMVbKiwFGBBR6ivta5X1NtqXF2Pgj4jKbtpPkXsJUADp+R+zLicJbRlP8fmrKGr7KYxJ59c1WA18GiUHGPNdWz+zD5PLGUjOVZZbV91oHMIo5Z8CvJyoKgX9DoeTriyeV4j1TrnRNlGU7JOlHwQT+YNlLIU7Z5N1iisfXfITe90W6Ty/ulTy6KW4T3/NSRlHZCqoj4MeN5+mwc8qZT6XazvJK2on7vSGhL9ZLn1vqkJbu8Bw6+w8niEkZJ1qzu1uy03R4dDop7iyfPvWH1gyOmAeMP6eAon4udDejJmxn9Yt6ueRg64vzqwj6NkPXk08S2dWKusskrYqosjew5NlMkmtub1YOb4bpyeu9ga7XQus+ZQ6vZZFntxH8s6i2VNNTXy6qIP+Ocr77CvEfarALXkszmvlDOPO9xRSoNkSMWijlTPAjkCQov8KVnT9qIQqwzT8eyxFHXcbkgp9Q0wOOm7J0JuvjUJEGR3pTXU7nJ4xNN1nJ1NGwUdrb8YePL8YVEv8YgX+hQrciRWWNRPzhh0kLLYS1sWq4OjGZINRYu2qKiJHFapUqiH37xTzenTDkwYz/1GQjqWDVw/vl3sss/J5bSTTqS6/aEtOqQz+5W0WOrvRZa8UGJFhkQLzYx0LOPbXQxilaHbekev8Ly8AksxB9lhZdKjc2RFHSsaoGzaS1lfeXSMaU0ljjfebiPBc4KKYUyYcoPUFuGEKiQnmRBTibcN77BGzViYcqcba7TiZLVkvJWFkT7P1rYVieCzurleIRqaKeqw8Lwdq63/Uay6WMlvQvMcQHZWqFSUYqo+z2jfTyWONxlrvPzQzmm16oL3iDbMDS3bdI5oUu10nXQauqwM9IN0zStNGlqa0KrSdKGfog61qHdvtCZhOkSeAHOS/8JLV4jXE5vJKsV4jTraysdQa6y6rqHZPxlJKSRTDsl0PG4pHydlm84RTaorCrPaDZgiydT3WGXmRzoFvRR1bhsrprKpyYpy2L3RntENRDw9XClEmxb1whXgx7LTZJVivCTu4c9x/ZzlLSaOIqV3DVUKySpPv/OJhOKkbNO5XDvVZ9fRDaYLidb3eO3WjwRTeinq4CKBxv2QU3jw1lsRCFUKToarbhGtMkx9ahl3zf8i5R8yWq+fjFKM1agjPUd9k7Mo2mxLAh+vbNPZsaT67LpkNwxHh+RGidZ3JyMRr91Imilqe+VYQ621WGD3RmsRiEP8tMhiKalUret0W+uxGnUqytaNJPDBvMFOGroTpZBOxZHujiWVxq/TaCSILsmNkqnvuo1EtNg4oJlc26IOhujt3hjXog5l0lBv9n+DlonHh/zmlbgLkVLZHy7d+83FSuKerLJ1c0NUJ0nwnZzrRlL9SUNLWTRtLKtnTGDRtLG+Wv9e1X2nuL3fo1OSqe9+j0TC0deirt1trYTr2DOhS3gxJIm1y3Qsku2l093rx7MEnSxuCOQI7QvyHId3JUMiw1In52b7hJtuUR26WKuJ1ne/RyKR0ExRt7H+N9ZZ/mlIWFF7QbxdpqORbC/thv8xXlysDosbEmnoTs7VRXG0FnTymyda33Xq8EA3Rd3s+qi1ErGAlfdAE0J3GUmUVHppr/2PuixuSKShOzlXJ8WRzaRjN57QuYREF+ckim4jkUjo5aNudn3UWf5p0MaiDvVvOqG4MJA2f6GO/kcvSGRDVCfn6rbBajYS3k6Cu/GA83obPpdQVVPPzn31rm3WmwnoZVGHhuft3oiVvjB6UiK3CLUIckVoVKr5vxMKA7lMPyfybh2h109kqJUJvX66SWRY6uTcTBnmZjLp2I0nnmsxm+YVnKKZog6ZTNxdaeWYjZGA3Q3CJwqDyjmWki5OYGiW7lVS6UKHeNdIJNJBOTm3NXZ44bj5W6djHsDJuW7s7pMsXrQdvRR1aHhegqF5cHCBjelXknDayEQnChPNvxst8uDnTy+nYu0OVzOoRUOXeFeD+8T6rSH10Ua0eYCiwoDjePhYOXyCiP0sftdPr9qOZj5qO+qjodbatSIB/3SkGNl/vL/uoJjZX839pDn+edSMhQf5uhLp+ZPxb0a7fqNSPPH+Ol/iTnWJd3WD0Hj3SL93JpLKM0X7rafP+ywtMeaR5gECOUJ1XYPja0e6RjjKfha/8artaKaobddHY521x1/77o6/6sQSrqlv5IkIyju0wsSLAMgVSWlCL9b1/cpVkq1ha24scElGhnR2FKk+U7TftKqmPi0KJ9LEd/uCvBYx+PGuHX6NaOhQP71qO5q6PmqtPdIKix1/1WnBhCvD8ImJWBn50rGLg5OMf+G4HT6ma9haqr4/vxe4uDEsTvWZnLgVQklG4YTPA5RNeynha7uRw8eNFANetR09LeqaKiuLXpvYO5yEkkrBhFaY0N4cLAsa0hcSF7x+rsPNL70IH9MxbC0d1rDfIwU3hsWpPlO037pT28gZKtOhcGIt03Yy4khH/XQrxYBXbUcvizoY4VG9xfpfUOT4q04s1cibwR9ckdyODIi2dDWcUo+iL3QMW0uHNez3SMGNjiLVZ4q18tStRVXRFmyN6VfiaMSRjvrpVooBr9qOXoo6155M3BtU1M4t6tACi7QiSoATD+/MknW7tFjXH5T3508vjxj658luzmHyhIYH/jRN6VmTJR1Kzu+Mcm50FOl4pliGiBsKJ5oyS0Qxpmo8xdq2Lxg9kmyd8yLkUy9FHYz6qN5q/U/AoobY2ygpYM32Gu44f6A2lqMuSWGiLfn1M0wvHUrO75GCGx2Fm8/kpsKJdO2ferj3YCzffLCO+z0Ci4VjRS0iuUAFUKmUOtsVaUSsCcW9QUVdnNRlYvWMXi54iLTCMdyd4WbDczp5EqpM4k22ukWkGPh0bFzr5wIXt37bbFm046VidLKDuN8jsFgkYlH/BFgBOPdHJENewYGETAlMJoaiQ88YbYWjV5uOOo04cBLW6PbkWyRZn11cyQXHlSa8YEk3skWpuoGXijH4G8TaQdzvEVgsHClqEekFTAB+B/zMVYly82HfNut1gq6PIDr0jLEUoBdWqlP/nxMl7HYHF03W11du9dRPnyjJhA/6sVRf5/QA4J1inDQ0/g7ioR1r6HyN3+Xm1KK+B/gl0ME9UWzy2oBqsl4nMJkYig49YzwFqMsilnhxtV50cMlO4vipgJKJkfZjqb7u6QFijTjc+H2dGnG6lVtcRS0iZwNblFKLRWR0jPOuBK4E6NOnTwoS2ROKufkH4qqTwO8hZzwF6LaVWtw2wM599XHvG6niBicUvQoPTMZVle6GlKhSSCaUy48FOMne0+/kYG5NbDs14vxeLBWOE4t6FHCOiJwFFAAdReQfSqnvh56klHoAeACgvLzcWT7QSARD9AqKiLsRocbEmrwQYEy/EtfuPXdpJXtrGw46HsiVgywHHUYfybiq0tmQklH6yYwC/FiAk8w9/bImvZrYdmLE+b1YKpy4ilopdSNwI4BtUf8iXEmnVyJbUSc5kagL4XHdoSjg2cWVlB/a2ZWKf9f8L6hvOrivbJefF/F+fo8+kuks0tmQklH6yYwCnI5y0kkycvplTeowsR1Eh4CEUPRaQg4HFHWSE4k6MWmotUN1aYQf183sdNEq8y6Hm/D6QbCsnO7mnc7do5NR+okuHU5klJNOklni7Jc1qcPEdhDd0ioktOBFKfUG8IYrkgQJJmZKciJRR7yu+LpZA26QiLsknr81mfJKdBSQ6CgnXSQzWkl2zsCtXNZBgr+vF/5zHVyCoei1MhEOTCDGsah1DTmKhNeKU4fwRLdx2pCc+FuTLa9EXEZujXKctINEXVuJlke6fNpOJrYBz/znfrsEQ9FQUcd3fegWOhMPP3YRB32sAbdw0pCc+Fu9KC83Omu32kEyo4V0+LSd3Hfora9oFY3hFfoq6hiTibqFzsTDD8WpkzXgJ07dTm6XlxudtZvtIB2jhXTksg5l7tLKiJOxyd4rk9BPUTeH5xVHPUW30BknGMXpD7r4693orHVpB16VcazJdy9+Tz/drfop6mBO6hiuD10an0F//PDXR2vQ6e6sdWkHXpVxvCgcN/Hb3apheF5wMjG660O30BmDvkwaevAefunYqScaXu7TqEs78KqMo3VAxYWBiPdK536Vfm8ArZ9FnRvfom4tk2WG9OCl28nL+ROd2oEXZRzNcp9+zrEHnZtuC9hvN5N+ijpoUcdZmdjafb6ZFJ4YjWx4hnC8btCtqR0k0jGlu8P0282koaKOb1G3dvz2l6WDbHiGSPjdoLMdpx1TujtMv9cm6OejDk3KZIhIuvxl6fThJYrfPj+30MVv3NpJZ4oB8H6uIxz9LOru/aHLkdDOvexymU46rAW/LVq/fX5uoZPfuDXj1n6VZgl5kCNOhesq/JZCa9IxvPZ70VA2uwhak99YV7Ktw9RPURvikg5rwW+L1m+fnyH7yaYO0yjqDCQd1oLfFm22WTwGg5uIUslvxhKN8vJyVVFh3Bc6E+6jBsui9XKCxJB5ZGNIpS6IyGKlVHmkz4xF3UoxFq0hUfyegG7NtApFbayAyGSTD8/gPn5PQLdmsl5RGyvAkCqmo7fwewK6NaPfgpc0k60LKwze4GWSJd1J9yKScPxcgKU7WW9RGyvAkApmuH+AdIVURhqhgHdbbGUiWa+o/Q5DM2Q2pqM/QDomoKO5IgsCOaZDjEHWK2qzsMKQCqajb0mqE9DRRijhx4K0xg4xEnF91CJSICIfishyEflMRH7jhWDpwu9kKobMxiRZSi+JKt7W2iGG48Si3g+MVUrtFZEA8I6I/Fsp9b7LsqUNE4ZmSBYTb55eoo1QAAQIXX5nOsQDxFXUylq6uNd+G7D/0r+c0WDQFNPRp4/rxx/N1KeWRfxMYY14TYd4MI581CKSCywGjgD+opT6IMI5VwJXAvTp0yedMhoMhixh0tBSfvOvz9i5r/6gz0qLC1k0bawPUumPozhqpVSjUmoI0As4XkQGRDjnAaVUuVKqvKTE5JI2GAyRuWXiscbvnyAJRX0opapE5A3gDOBTVyQymJVwhqzG+P0TJ66iFpESoN5W0oXAqcCdrkvWSjFL3g2tAeP3TwwnFnUP4DHbT50DPK2UetFdsVovZiWc/5gRjUE3nER9fAwM9UAWA2YlnN+YEY1BR7I+KVOm4XbiG0NsTBIvg44YRa0ZZiWcv5gRjUFHjKLWDLPk3V/MiMagI1mflCkTMTPi/mGSeBl0xChqgyEEE+Nr0BGjqA2+oHMIXLpHNPX19WzYsIHa2tq0XdOQuRQUFNCrVy8CgYDj7xhFbYiKW8o0Xgiczko8GTZs2ECHDh3o27cvIuK3OAYfUUqxfft2NmzYQFlZmePvGUVtiIib8cTxQuB0jmNOphOpra01StoAgIjQpUsXtm7dmtD3TNSHISJuxhPHCoHTOY45lY1ujZI2BEmmLhhFbYiIm/HEsULgdI5j1rkTMWQ3RlEbIuJmPHGsRT06xzF71YnMXVrJqBkLKZv2EqNmLHRksRvi88Ybb/Duu+/6LUZSGEVtiIibKyRjLerReWWmF51IKu6VTKKhocHzexpFbcg63F4hOWloKYumjWX1jAksmja2+bo6r8z0ohNxy72yZs0a+vXrxxVXXMGAAQOYPHkyCxYsYNSoURx55JF8+OGHVFdXc9lllzF8+HCGDh3KCy+80Pzdk046iWHDhjFs2LBmZbdp0yZOPvlkhgwZwoABA3j77bcBaN++ffN9n3nmGaZMmQLAlClT+NnPfsaYMWO44YYbWLVqFWeccQbHHXccJ510EitXrmw+7+qrr2bMmDEcdthhvPnmm1x22WUcc8wxzdcCeOWVVxg5ciTDhg3jwgsvZO9ea8fAvn37cssttzBs2DAGDhzIypUrWbNmDbNmzeLuu+9myJAhvP3228yZM4cBAwYwePBgTj755JTK121M1IchKn6tkNR1ZaYXi2HcdK98/fXXzJkzhwceeIDhw4fz5JNP8s477zBv3jxuv/12+vfvz9ixY3n44Yepqqri+OOP59RTT6Vbt268+uqrFBQU8NVXX3HxxRdTUVHBk08+yfjx47nppptobGxk3759cWX48ssvWbBgAbm5uYwbN45Zs2Zx5JFH8sEHH3DNNdewcOFCAHbu3MnChQuZN28eEydOZNGiRTz44IMMHz6cZcuW0atXL2677TYWLFhAu3btuPPOO/njH//IzTffDEDXrl1ZsmQJ999/PzNnzuTBBx/kqquuon379vziF78AYODAgcyfP5/S0lKqqqpSLl83MYraYEgAtzuRaLt0p8O9UlZWxsCBAwE49thjGTduHCLCwIEDWbNmDRs2bGDevHnMnDkTsMIK161bR8+ePbn22mtZtmwZubm5fPnllwAMHz6cyy67jPr6eiZNmsSQIUPiynDhhReSm5vL3r17effdd7nwwgubP9u/f3/z64kTJzbL1r179xZyB2X9/PPPGTVqFAB1dXWMHDmy+fvnn38+AMcddxzPPfdcRFlGjRrFlClT+N73vtd8vq4YRW0waISbuUbatGnT/DonJ6f5fU5ODg0NDeTm5vLss89y9NEt7zV9+nS6d+/O8uXLaWpqoqCgAICTTz6Zt956i5deeolLLrmE66+/nh/84Actws/CV2O2a9cOgKamJoqLi1m2bFlMWUPlDJf1tNNOY/bs2TG/n5ubG9UfPmvWLD744ANeeuklhgwZwrJly+jSpUvEc/3G+KgNBo3w00c/fvx4/vSnP6GUAmDp0qUA7Nq1ix49epCTk8Pf//53GhutTmTt2rV069aNH/7wh1x++eUsWbIEgO7du7NixQqampp4/vnnI96rY8eOlJWVMWfOHMBasbd8+XLHso4YMYJFixbx9ddfA7Bv375mSz8aHTp0YM+ePc3vV61axQknnMCtt95K165dWb9+veP7e42xqA0GzfDLR//rX/+aqVOnMmjQIJRS9O3blxdffJFrrrmGCy64gDlz5jBmzJhmq/iNN97grrvuIhAI0L59ex5//HEAZsyYwdlnn03v3r0ZMGBA8yRfOE888QRXX301t912G/X19Vx00UUMHjzYkawlJSU8+uijXHzxxc0uk9tuu42jjjoq6ncmTpzId7/7XV544QX+9Kc/cffdd/PVV1+hlGLcuHGO7+0HEuw900l5ebmqqKhI+3UNhkxkxYoVHHPMMX6LYdCISHVCRBYrpcojnW9cHwaDwaA5RlEbDAaD5hhFbTAYDJoTV1GLSG8ReV1EVojIZyLyEy8EMxgMBoOFk6iPBuDnSqklItIBWCwiryqlPndZNoPBYDDgwKJWSm1SSi2xX+8BVgD6re81GAyGLCUhH7WI9AWGAh9E+OxKEakQkYpEdy8wGAz68Oijj7Jx48bm91dccQWff576AHrNmjU8+eSTze8rKir48Y9/nPJ1YzFnzhyOOeYYxowZ4+p93MaxohaR9sCzwFSl1O7wz5VSDyilypVS5SUlJemU0WAweEi4on7wwQfp379/ytcNV9Tl5eXcd999KV83Fg899BD3338/r7/+uqv3cRtHKxNFJIClpJ9QSkXOcGIwGOLz72nw7SfpveYhA+HMGTFP+cc//sF9991HXV0dJ5xwAvfffz8Al19+ORUVFYgIl112Gb1796aiooLJkydTWFjIe++9x5lnnsnMmTMpLy+nffv2/OhHP2LBggV06tSJ22+/nV/+8pesW7eOe+65h3POOYc1a9ZwySWXUF1dDcCf//xnTjzxRKZNm8aKFSsYMmQIl156KUOHDmXmzJm8+OKL7Nixg8suu4xvvvmGtm3b8sADDzBo0CCmT5/OunXr+Oabb1i3bh1Tp06NaIXPnj2b22+/HaUUEyZM4M477+TWW2/lnXfeYfXq1Zxzzjncddddzee/8cYb3HLLLXTv3p1ly5Zx/vnnM3DgQO69915qamqYO3cuhx9+OFu3buWqq65i3bp1ANxzzz2MGjWKDz/8kKlTp1JTU0NhYSGPPPIIRx99NI8++ijz5s1j3759rFq1ivPOO4/f//73qf/GSqmYf4AAjwP3xDs3+Hfccccpg8Fg8fnnnx948/INSj18Vnr/Xr4h7v3PPvtsVVdXp5RS6uqrr1aPPfaYqqioUKeeemrzeTt37lRKKXXKKaeojz76qPl46HtAvfzyy0oppSZNmqROO+00VVdXp5YtW6YGDx6slFKqurpa1dTUKKWU+vLLL1VQH7z++utqwoQJzdcNfX/ttdeq6dOnK6WUeu2115qvdcstt6iRI0eq2tpatXXrVtW5c+fm5whSWVmpevfurbZs2aLq6+vVmDFj1PPPPx/xWULvXVRUpDZu3Khqa2tVz5491c0336yUUuqee+5RP/nJT5RSSl188cXq7bffVkoptXbtWtWvXz+llFK7du1S9fX1SimlXn31VXX++ecrpZR65JFHVFlZmaqqqlI1NTWqT58+at26dRF/k3CAChVFpzqxqEcBlwCfiMgy+9j/KqVeTr2baEkyOzwbDBlFHMvXDV577TUWL17M8OHDAaipqaFbt25MnDiRb775huuuu44JEyZw+umnx71Wfn4+Z5xxBmDlc27Tpg2BQKA5VSpAfX19xLSosXjnnXd49tlnARg7dizbt29n165dAEyYMIE2bdrQpk0bunXrxubNm+nVq1fzdz/66CNGjx5N0OU6efJk3nrrLSZNmhTznsOHD6dHjx4AHH744c3PP3DgwGZXyYIFC1r453fv3s2ePXvYtWsXl156KV999RUiQn19ffM548aNo6ioCID+/fuzdu1aevfuHbcMYhFXUSul3sGyql0luAVRML1jcAsiwCjrLMd00O6ilOLSSy/ljjvuOOiz5cuXM3/+fP7yl7/w9NNP8/DDD8e8ViAQaE5jGilVKsDdd98dMS1qPBnDCd4nNM1ppLSlkb7rhHhpX8FKx/ree+9RWNgyH/h1113HmDFjeP7551mzZg2jR4+OeN1YaVYTQZuViWaH59ZJa9kj0E/GjRvHM888w5YtWwDYsWMHa9euZdu2bTQ1NXHBBRfw29/+tjlNaXg60ESJlhY11nVPPvlknnjiCcDyH3ft2pWOHTs6ut8JJ5zAm2++ybZt22hsbGT27NmccsopScsfyumnn86f//zn5vfB/Nm7du2itNQyJh599NG03CsW2ihqr3Z4NuiF6aDdp3///tx2222cfvrpDBo0iNNOO41NmzZRWVnJ6NGjGTJkCFOmTGm2uKdMmcJVV13FkCFDqKlJvP1dc801PPbYY4wYMYIvv/yyOS3qoEGDyMvLY/Dgwdx9990tvjN9+nQqKioYNGgQ06ZN47HHHnN8vx49enDHHXcwZswYBg8ezLBhwzj33HMTljsS9913X7Nc/fv3Z9asWQD88pe/5MYbb2TUqFHNHZGbaJPmdNSMhRG3ICotLmTRtLHpEs2gGWXTXiJSDRRg9YwJXovjCibNqSGcjE1z6sUOzwb9iLYXYDr2CDQYsgVtFLWfWxAZ/MN00AZDfLTaisuvLYgM/hH8vbM96kMp1WLTV0PrJRl3s1aK2tA6yfYOuqCggO3bt9OlSxejrFs5Sim2b9/uKGQxFKOoDQaX6dWrFxs2bMAkKzOA1XGHLthxglHUBoPLBAIBysrK/BbDkMFoM5loMBgMhsgYRW0wGAyaYxS1wWAwaI4rKxNFZCuwNsGvdQW2pV2Y9KCrbEauxDByJY6usmWjXIcqpSLuuuKKok4GEamItnzSb3SVzciVGEauxNFVttYml3F9GAwGg+YYRW0wGAyao5OifsBvAWKgq2xGrsQwciWOrrK1Krm08VEbDAaDITI6WdQGg8FgiIBR1AaDwaA5WihqETlDRL4Qka9FZJqPcvQWkddFZIWIfCYiP7GPTxeRShFZZv+d5YNsa0TkE/v+FfaxziLyqoh8Zf/v5LFMR4eUyTIR2S0iU/0qLxF5WES2iMinIceilpGI3GjXuS9EZLzHct0lIitF5GMReV5Eiu3jfUWkJqTsZnksV9TfzufyeipEpjUissw+7mV5RdMP7tcxpZSvf0AusAo4DMgHlgP9fZKlBzDMft0B+BLoD0wHfuFzOa0BuoYd+z0wzX49DbjT59/xW+BQv8oLOBkYBnwar4zs33U50AYos+tgrodynQ7k2a/vDJGrb+h5PpRXxN/O7/IK+/wPwM0+lFc0/eB6HdPBoj4e+Fop9Y1Sqg74J5CenSkTRCm1SSm1xH69B1gB6Jwo+VwguAvoY8Ak/0RhHLBKKZXoitS0oZR6C9gRdjhaGZ0L/FMptV8ptRr4GqsueiKXUuoVpVSD/fZ9ILG8ly7JFQNfyyuIWAm9vwfMduPesYihH1yvYzoo6lJgfcj7DWigHEWkLzAU+MA+dK09TH3YaxeDjQJeEZHFInKlfay7UmoTWJUI6OaDXEEuomXj8bu8gkQrI53q3WXAv0Pel4nIUhF5U0RO8kGeSL+dLuV1ErBZKfVVyDHPyytMP7hex3RQ1JG2vPA1ZlBE2gPPAlOVUruBvwKHA0OATVhDL68ZpZQaBpwJ/EhETvZBhoiISD5wDjDHPqRDecVDi3onIjcBDcAT9qFNQB+l1FDgZ8CTItLRQ5Gi/XZalBdwMS0NAs/LK4J+iHpqhGNJlZkOinoD0DvkfS9go0+yICIBrB/hCaXUcwBKqc1KqUalVBPwN1wa8sVCKbXR/r8FeN6WYbOI9LDl7gFs8VoumzOBJUqpzbaMvpdXCNHKyPd6JyKXAmcDk5Xt1LSHydvt14ux/JpHeSVTjN9Oh/LKA84Hngoe87q8IukHPKhjOijqj4AjRaTMtswuAub5IYjt/3oIWKGU+mPI8R4hp50HfBr+XZflaiciHYKvsSaiPsUqp0vt0y4FXvBSrhBaWDl+l1cY0cpoHnCRiLQRkTLgSOBDr4QSkTOAG4BzlFL7Qo6XiEiu/fowW65vPJQr2m/na3nZnAqsVEptCB7wsryi6Qe8qGNezJY6mE09C2sGdRVwk49yfAdraPIxsMz+Owv4O/CJfXwe0MNjuQ7Dmj1eDnwWLCOgC/Aa8JX9v7MPZdYW2A4UhRzzpbywOotNQD2WNXN5rDICbrLr3BfAmR7L9TWW/zJYz2bZ515g/8bLgSXARI/livrb+Vle9vFHgavCzvWyvKLpB9frmFlCbjAYDJqjg+vDYDAYDDEwitpgMBg0xyhqg8Fg0ByjqA0Gg0FzjKI2GAwGzTGK2mAwGDTHKGqDwWDQnP8PPQ9J63Lr5N0AAAAASUVORK5CYII=",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"\"\"\" Kalman filter: recursive estimation of mean\n",
"\"\"\"\n",
"\n",
"npt = 200\n",
"x = 5 + np.random.randn(npt)\n",
"n = np.arange(1,npt+1)\n",
"K = 1.0/(n+1)\n",
"mu = np.zeros((npt))\n",
"mu[0] = x[0]\n",
"\n",
"for i in range(1,npt):\n",
" mu[i] = mu[i-1] + K[i-1]*(x[i]-mu[i-1])\n",
"\n",
"plt.plot(n,x,'o',label='measurements')\n",
"plt.plot(n,mu,label='estimation of mean')\n",
"plt.legend()\n",
"plt.title('Recursive estimation of mean')\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This example showed an important feature of the Kalman filter: it provides **recursive estimates** of the system state in terms of measurements.\n",
"\n",
"Next let us consider another important feature of the Kalman filter: it provides *best* estimates."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Best estimates\n",
"\n",
"Suppose $x$ is measured with *two different tools* (e.g. when $x$ is a distance, it is measured both with a measuring tape and with a laser distance meter). Our first tool provides the measurement $x_1$, our second $x_2$. These tools do not have the same precision. Their precisions can be expressed by the standard deviations (mean errors) $\\sigma_1$ and $\\sigma_2$.\n",
"\n",
"Our goal is to provide *best* estimate of $x$ in terms of measurements $x_1$, $x_2$ and their standard deviations $\\sigma_1$, $\\sigma_2$. Let us express this as a *weighted average* of the two measurements:\n",
"\n",
"$$\\hat{x} = wx_1 + (1-w)x_2.$$\n",
"\n",
"Our problem is now what weight $w$ makes standard deviation $\\hat\\sigma$ of the estimate $\\hat{x}$ minimal? \n",
"\n",
"Since we suppose *Gaussian distribution* of both $x_1$ and $x_2$, according to a [well-known result](https://www.statlect.com/probability-distributions/normal-distribution-linear-combinations) in probability their linear combination is Gaussian as well with variance \n",
"\n",
"$$\\hat\\sigma^2 = w^2\\sigma_1^2 + (1-w)^2 \\sigma_2^2.$$\n",
"\n",
"The best (optimal) weight $w_{opt}$ is that minimizes the above expression, i.e. where the derivative with respect to $w$ is zero:\n",
"\n",
"$$\\frac{\\partial \\hat\\sigma^2}{\\partial w} = 2w_{opt}\\sigma_1^2-2(1-w_{opt})\\sigma_2^2 = 0.$$\n",
"\n",
"The optimal weight is\n",
"\n",
"$$ w_{opt} = \\frac{\\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}.$$\n",
"\n",
"Optimal estimate $\\hat{x}$ of $x$ is:\n",
"\n",
"$$ \\hat{x} = \\frac{\\sigma_2^2}{\\sigma_1^2+\\sigma_2^2} x_1 + \\frac{\\sigma_1^2}{\\sigma_1^2+\\sigma_2^2} x_2.$$\n",
"\n",
"Variance of the optimal estimate is:\n",
"\n",
"$$ \\hat\\sigma^2 = \\frac{\\sigma_1^2\\sigma_2^2}{\\sigma_1^2+\\sigma_2^2}.$$\n",
"\n",
"For us in geodesy this is the well-known weighted least squares estimation, since $\\hat{x}$ minimizes the sum of weighted residuals:\n",
"\n",
"$$\\hat{x} = \\underset{x}{\\arg \\min} = \\left[ p_1 (x_1 - x)^2 + p_2 (x_2 - x)^2 \\right],$$\n",
"\n",
"where $p_1$, $p_2$ denotes weights proportional to reciprocal variances. In fact if we make derivative of the above expression with respect to $x$ and make it equal to zero the equation \n",
"\n",
"$$ p_1(x_1-\\hat{x}) + p_2(x_2-\\hat{x}) = 0$$\n",
"\n",
"is yielded. Its solution for $\\hat{x}$ is\n",
"\n",
"$$ \\hat{x} = \\frac{1}{p_1+p_2}(p_1 x_1 + p_2 x_2)$$\n",
"\n",
"which is identical with the previous solution when $p_1=1/\\sigma_1^2$ and $p_2=1/\\sigma_2^2$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recursive best estimate\n",
"\n",
"Let us recast our best estimate in the following way: first we use measurement $x_1$ for estimation, then update (correct) our estimate by using the next measurement $x_2$. Obviously in the first step $\\hat{x_1} = x_1$, $\\hat\\sigma_1^2=\\sigma_1^2$, where subscript refers to the step. In the second step we update our estimate as\n",
"\n",
"$$\\hat{x}_2 = \\hat{x}_1 + \\frac{\\hat{\\sigma}_1^2}{\\hat{\\sigma}_1^2 + \\sigma_2^2} (x_2 - \\hat{x}_1),$$\n",
"\n",
"$$\\hat{\\sigma}_2^2 = \\left( 1 - \\frac{\\hat{\\sigma}_1^2}{\\hat{\\sigma}_1^2 + \\sigma_2^2} \\right) \\hat{\\sigma}_1^2.$$\n",
"\n",
"The factor in parentheses $(x_2 - \\hat{x}_1)$ is *innovation* (new information), and its multiplier, \n",
"\n",
"$$\\frac{\\hat{\\sigma}_1^2}{\\hat{\\sigma}_1^2 + \\sigma_2^2}$$\n",
"\n",
"is the *gain* $K$. Let us summarize our recursive estimation scheme:\n",
"\n",
"$$\\hat{x}_2 = \\hat{x}_1 + K(x_2 - \\hat{x}_1),$$\n",
"\n",
"$$\\hat{\\sigma}_2^2 = ( 1 - K ) \\hat{\\sigma}_1^2,$$\n",
"\n",
"$$K = \\frac{\\hat{\\sigma}_1^2}{\\hat{\\sigma}_1^2 + \\sigma_2^2}.$$\n",
"\n",
"As we can see it has the same form as our previous recursive estimate\n",
"\n",
"$$\\mu_{n+1} = \\mu_n + K(x_{n+1} - \\mu_n),$$\n",
"\n",
"but now we have found the *best* (optimal) value of the gain $K$ that minimizes estimation variance. This best estimate can always be updated when a new measurement is available.\n",
"\n",
"At this point we remark that if state $x$ is not a scalar but a vector $\\mathbf{x} = (x_1, x_2, ...)$, the procedure is the same, but we must use covariance matrices instead of variances. Among others, we will form inverses of covariance matrices instead of reciprocals of variances."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Indirect state estimation\n",
"\n",
"Until now it was supposed that we can *directly* measure the system state $x$. In most cases this cannot be made - only a quantity $\\mathbf{z}$ can be measured, which is related to the system's state $x$. This problem is well-known in geodesy: least squares adjustment with parameters. A similar problem is solved by the Kalman filter: estimation of unobservable (hidden) system states with measurements. \n",
"\n",
"First we deal with *static* systems, the state of those do not change with time. But we do have a key assumption for both static and dynamic systems: system state $\\mathbf{x}$ is connected with measured $\\mathbf{z}$ with a **linear** model:\n",
"\n",
"$$ \\mathbf{z} = \\mathbf{H} \\mathbf{x} + \\mathbf{r}.$$\n",
"\n",
"In this equation matrix $\\mathbf{H}$ is the measurement model (the well-known design matrix), moreover measurement residuals $\\mathbf{r}$ are zero mean ($\\bar{\\mathbf{r}}=\\mathbf{0}$) and Gaussian distributed. These residuals (noise) have noise *covariance matrix* $\\mathbf{R}$ \n",
"\n",
"$$ \\mathbf{R} = \\mathbb{E}\\left[ (\\mathbf{r} - \\bar{\\mathbf{r}}) (\\mathbf{r} - \\bar{\\mathbf{r}})^T \\right]$$\n",
"\n",
"where $\\mathbb{E}[.]$ denotes expectation value. We want best (optimal) estimation $\\hat{\\mathbf{x}}$ of the system state that minimizes measurement residuals (noise) weighted by the inverse covariance matrix. This is the very well known adjustment with parameters (now system state $\\mathbf{x}$ are the parameters), where weight matrix $\\mathbf{P}$ is proportional with inverse $\\mathbf{R}^{-1}$ of the noise covariance matrix."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Recursive estimation**\n",
"\n",
"In contrast with traditional adjustment with parameters let us suppose a continuous flow of measurements, i.e. they are not available at once. At each measurement epoch $k$ a new measurement $\\mathbf{z}_k$ is available and we need best estimates of the state $\\hat{\\mathbf{x}}_k$ and state covariance matrix $\\hat{\\mathbf{P}}_k$. Our minimization problem is therefore:\n",
"\n",
"$$\\hat{\\mathbf{x}}_k = \\underset{x}{\\arg \\min} = \\left[ (\\mathbf{z}_k-\\mathbf{H}_k \\mathbf{x})^T \\mathbf{R}_k^{-1} (\\mathbf{z}_k-\\mathbf{H}_k \\mathbf{x}) \\right].$$\n",
"\n",
"*Static Kalman filter* provides solution to this minimization problem. It provides new estimate $\\hat{\\mathbf{x}}_k$ and its covariance matrix $\\hat{\\mathbf{P}}_k$ from previous estimates $\\hat{\\mathbf{x}}_{k-1}$, $\\hat{\\mathbf{P}}_{k-1}$ and measurement $\\mathbf{z}_k$ and its covariance matrix $\\mathbf{R}_k$. \n",
"\n",
"A typical example of static Kalman filtering would be position estimation of a static GNSS receiver. In this case the hidden system state (parameter vector $\\mathbf{x}$) contains unkonwn position and clock bias of the receiver. Measurements $\\mathbf{z}$ are the satellite-receiver pseudoranges. When all measurements are available this is a problem in traditional least squares estimation. When there is a continuous flow of measurements epoch by epoch, the positioning requires recursive least squares adjustment, i.e. static Kalman filtering.\n",
"\n",
"Next we derive equations of the static Kalman filter by *adjustment in groups* used in geodesy. It will be recognized that equations of the static Kalman filter are identical with that of parametric adjustment in groups.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Equations of the static Kalman filter\n",
"\n",
"This derivation follows that of Section 5.8.2 found in textbook of Detrekői: Adjustment calculations (Kiegyenlítő számítások, in Hungarian). Our notation, however, will be more close to the notation used in Kalman filtering literature, to ease comparision of the two procedures. We restrict our attention to the linear case, hence no linearization of equations for the parameters is required.\n",
"\n",
"Within the procedure for adjustment in groups there are measurements at different epochs and measurements at a later epoch are processed by using results obtained from previous epochs. Let us asume that at epoch $i=1$ by using $n$ measurements $\\mathbf{z}_1$, adjusted parameters $\\mathbf{x}_1$ and their covariance matrix $\\mathbf{P}_{11}$ have already been determined by using the covariance matrix $\\mathbf{R}_{11}$ of measurements from linear measurement (residual) equations \n",
"\n",
"$$ \\mathbf{r}_1 = \\mathbf{z}_1 - \\mathbf{H}_1\\mathbf{x}_1 $$\n",
"\n",
"through the normal equations\n",
"\n",
"$$ \\mathbf{N}_1 \\mathbf{x}_1 - \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 = \\mathbf{0}.$$\n",
"\n",
"If the coefficient matrix $ \\mathbf{N}_1= \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{H}_1$ of the normal equations is nonsingular, then the unknown parameters are\n",
"\n",
"$$\\mathbf{x}_1 = \\mathbf{N}_1^{-1}\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 .$$\n",
"\n",
"Covariance matrix of estimated parameters is\n",
"\n",
"$$ \\mathbf{P}_{11} = \\mathbf{N}_1^{-1} = (\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{H}_1)^{-1}.$$\n",
"\n",
"In the next epoch $i+1=2$ additional $s$ measurements $\\mathbf{z}_2$ are made, their covariance matrix is $\\mathbf{R}_{22}$. Let us assume that these measurements are independent of the measurements at the previous epoch. Let us calculate the updated parameter vector $\\mathbf{x}_2$, updated measurement residual vector $\\mathbf{r}_r$ and the covariance matrix $\\mathbf{P}_{22}$ of the updated parameter vector."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Updated normal equation is\n",
"\n",
"$$ (\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{H}_1 + \\mathbf{H}_2^T\\mathbf{R}_{22}^{-1}\\mathbf{H}_2)\\mathbf{x}_2 - ( \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 + \\mathbf{H}_2^T\\mathbf{R}_{22}^{-1}\\mathbf{z}_2) = \\mathbf{0}.$$\n",
"\n",
"Let us introduce a size $s$ new parameter $\\mathbf{y}$ vector:\n",
"\n",
"$$ \\mathbf{y}= \\mathbf{R}_{22}^{-1}\\mathbf{r}_2 = \\mathbf{R}_{22}^{-1}\\mathbf{z}_2 -\\mathbf{R}_{22}^{-1}\\mathbf{H}_2 \\mathbf{x}_2 .$$\n",
"\n",
"Then, with the help of $\\mathbf{y}$, the equation of residuals\n",
"\n",
"$$ \\mathbf{r}_2 = \\mathbf{z}_2 - \\mathbf{H}_2\\mathbf{x}_2 $$\n",
"\n",
"can be transformed into a parametric constraint equation:\n",
"\n",
"$$ \\mathbf{R}_{22}\\mathbf{y}= \\mathbf{z}_2 - \\mathbf{H}_2\\mathbf{x}_2,$$\n",
"\n",
"from which we can derive the usual form of constraint equations:\n",
"\n",
"$$ \\mathbf{z}_2 - \\mathbf{H}_2\\mathbf{x}_2 - \\mathbf{R}_{22}\\mathbf{y} = \\mathbf{0}$$\n",
"\n",
"Using vector $\\mathbf{y}$ we transform normal equations:\n",
"\n",
"$$(\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{H}_1)\\mathbf{x}_2 - \\mathbf{H}_2^T\\mathbf{y} - \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 = \\mathbf{0}$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Constraint equations together with transformed normal equations lead to the following system of equations with blocks of size $n+s$:\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{H}_1 & - \\mathbf{H}_2^T \\\\ - \\mathbf{H}_2 & - \\mathbf{R}_{22}\\end{bmatrix} \\begin{bmatrix} \\mathbf{x}_2\\\\ \\mathbf{y}\\end{bmatrix} = \\begin{bmatrix} \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1\\\\ -\\mathbf{z}_2\\end{bmatrix}. $$\n",
"\n",
"To solve for $\\mathbf{x}_2$ we invert its coefficient matrix $\\mathbf{A}$ with blocks by using Gauss-Jordan elimination. Our starting point is the augmented matrix\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{A} & | & \\mathbf{I} \\end{bmatrix} $$\n",
"\n",
"and with Gauss-Jordan elimination we transform it into the form\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{I} & | & \\mathbf{A}^{-1} \\end{bmatrix} $$\n",
"\n",
"where $\\mathbf{I}$ denotes identity matrix.\n",
"\n",
"Starting with the augmented matrix (for simplicity of notation $\\mathbf{N}=\\mathbf{N}_1$)\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{N} & - \\mathbf{H}_2^T & \\mathbf{I} & \\mathbf{0} \\\\ - \\mathbf{H}_2 & - \\mathbf{R}_{22} & \\mathbf{0} & \\mathbf{I}\\end{bmatrix} $$\n",
"\n",
"we multiply its first row from left with $\\mathbf{N}^{-1}$-el, and add to second row first row multiplied from left with $\\mathbf{H}_2$. Since $\\mathbf{H}_2\\mathbf{I}=\\mathbf{H}_2$ if $\\mathbf{I}$ is an $r\\times r$ identity matrix it yields the augmented matrix\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{I} & - \\mathbf{N}^{-1}\\mathbf{H}_2^T & \\mathbf{N}^{-1} & \\mathbf{0} \\\\ \\mathbf{0} & - \\mathbf{R}_{22}- \\mathbf{H}_2\\mathbf{N}^{-1}\\mathbf{H}_2^T & \\mathbf{H}_2\\mathbf{N}^{-1} & \\mathbf{I}\\end{bmatrix}.$$\n",
"\n",
"Let us introduce the notation $\\mathbf{S} = (\\mathbf{R}_{22}+ \\mathbf{H}_2\\mathbf{N}^{-1}\\mathbf{H}_2^T)^{-1}$. If we multiply from left by $-\\mathbf{S}$ its second row we get matrix\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{I} & - \\mathbf{N}^{-1}\\mathbf{H}_2^T & \\mathbf{N}^{-1} & \\mathbf{0} \\\\ \\mathbf{0} & \\mathbf{I} & -\\mathbf{S}\\mathbf{H}_2\\mathbf{N}^{-1} & -\\mathbf{S}\\end{bmatrix}.$$\n",
"\n",
"Finally, add to the first row the second row multiplied from left by $\\mathbf{N}^{-1}\\mathbf{H}_2^T$. It yields the augmented matrix\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{I} & \\mathbf{0} & \\mathbf{N}^{-1} - \\mathbf{N}^{-1}\\mathbf{H}_2^T\\mathbf{S}\\mathbf{H}_2\\mathbf{N}^{-1} & -\\mathbf{N}^{-1}\\mathbf{H}_2^T\\mathbf{S} \\\\ \\mathbf{0} & \\mathbf{I} & -\\mathbf{S}\\mathbf{H}_2\\mathbf{N}^{-1} & -\\mathbf{S}\\end{bmatrix}. $$\n",
"\n",
"On the right side of this matrix there is the inverse that was looked for.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we multiply with this inverse from left the equation to be solved we find the solution\n",
"\n",
"$$\\begin{bmatrix} \\mathbf{x}_2\\\\ \\mathbf{y}\\end{bmatrix} = \\begin{bmatrix} \\mathbf{N}^{-1} - \\mathbf{N}^{-1}\\mathbf{H}_2^T\\mathbf{S}\\mathbf{H}_2\\mathbf{N}^{-1} & -\\mathbf{N}^{-1}\\mathbf{H}_2^T\\mathbf{S} \\\\ -\\mathbf{S}\\mathbf{H}_2\\mathbf{N}^{-1} & -\\mathbf{S}\\end{bmatrix} \\begin{bmatrix} \\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1\\\\ -\\mathbf{z}_2\\end{bmatrix}. $$\n",
"\n",
"In the following, if we use the notation $\\mathbf{K} =\\mathbf{N}^{-1}\\mathbf{H}_2^T\\mathbf{S} $ then the solution $\\mathbf{x}_2$ is\n",
"\n",
"$$ \\mathbf{x}_2 = \\mathbf{N}^{-1}\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 - \\mathbf{K}\\mathbf{H}_2\\mathbf{N}^{-1}\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 + \\mathbf{K}\\mathbf{z}_2.$$\n",
"\n",
"This solution can be simplified by taking into account the relation $\\mathbf{x}_1 = \\mathbf{N}_1^{-1}\\mathbf{H}_1^T\\mathbf{R}_{11}^{-1}\\mathbf{z}_1 $. Finally we get\n",
"\n",
"$$ \\mathbf{x}_2 = \\mathbf{x}_1 - \\mathbf{K}\\mathbf{H}_2\\mathbf{x}_1 + \\mathbf{K}\\mathbf{z}_2 = \\mathbf{x}_1 + \\mathbf{K}(\\mathbf{z}_2 - \\mathbf{H}_2\\mathbf{x}_1 ).$$\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Substituting matrix $\\mathbf{S}$ into the equation (since additionally $\\mathbf{N}^{-1}=\\mathbf{P}_{11}$)\n",
"\n",
"$$ \\mathbf{K} = \\mathbf{P}_{11}\\mathbf{H}_2^T(\\mathbf{R}_{22} + \\mathbf{H}_2 \\mathbf{P}_{11} \\mathbf{H}_2^T)^{-1}.$$\n",
"\n",
"Covariance matrix of adjusted parameters is the upper left block of the inverted hypermatrix (the block belonging to $\\mathbf{x}_2$)\n",
"\n",
"$$ \\mathbf{P}_{22} = (\\mathbf{I} - \\mathbf{K}\\mathbf{H}_2) \\mathbf{P}_{11}.$$\n",
"\n",
"These are identical with the equations of the static Kalman filter. Kalman gain matrix is $\\mathbf{K}$, measurement matrix is $\\mathbf{H}$, system state is $\\mathbf{x}$. The filter updates state estimate and its covariance matrix $\\mathbf{P}$ after each new measurement.\n",
"\n",
"Equations of the *dynamic* Kalman filter are slightly different since also system state $ \\mathbf{x}$ is variable from epoch $k$ to epoch $k+1$:\n",
"\n",
"$$\\mathbf{x}_{k+1} = \\mathbf{F}\\mathbf{x}_k + \\mathbf{B}\\mathbf{u}_k + \\mathbf{q}$$\n",
"\n",
"where $\\mathbf{F}$ denotes state transition matrix, $\\mathbf{u}$ is system forcing, $\\mathbf{B}$ denotes system control matrix and $\\mathbf{q}$ is a zero-mean system noise with covariance matrix $\\mathbf{Q}$.\n",
"\n",
"Covariance matrix $\\mathbf{P}^-$ after state transition of the dynamical system is simply obtained with the well-known covariance propagation rule:\n",
"\n",
"$$\\mathbf{P}^-_k = \\mathbf{H} \\mathbf{P}_{k-1} \\mathbf{H}^T + \\mathbf{Q}.$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is interesting to note that equations of the Kalman filter can be derived directly from a probabilistic perspective (see explanation by [Pradu](https://math.stackexchange.com/questions/840662/an-explanation-of-the-kalman-filter)). A slightly modified explanation is [here](https://nbviewer.jupyter.org/github/gyulat/Korszeru_matek/blob/master/KF_derivation.ipynb)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Kalman filter example: vehicle navigation\n",
"\n",
"Consider a vehicle moving along a straight path. Our model of the system (vehicle) state $x$ contains position $p$ and velocity $v$ of the vehicle. System input $u$ are known control accelerations and output $y$ are position measurements. The system is controlled by accelerations and its position is measured at $T$ seconds. Basic kinematics provide the following equation of the velocity $v$:\n",
"\n",
"$$v_{k+1} = v_k + Tu_k.$$\n",
"\n",
"But our model is not yet complete. It is because velocity may have change randomly due to unavoidable external factors. Velocity noise $\\tilde{v}_k$ is a time-dependent random variable, hence \n",
"\n",
"$$v_{k+1} = v_k + Tu_k + \\tilde{v}_k$$\n",
"\n",
"is our complete velocity model. A similar equation holds for the position $p$:\n",
"\n",
"$$p_{k+1} = p_k + Tv_k + \\frac{1}{2}T^2u_k + \\tilde{p}_k,$$\n",
"\n",
"where $\\tilde{p}_k$ denotes position uncertainty (noise). State vector contains position and velocity:\n",
"\n",
"$$x_k = \\left[ p_k \\atop v_k \\right].$$\n",
"\n",
"Since we know that $z$ is vehicle's measured position the following system equation can be assembled:\n",
"\n",
"$$x_{k+1} = \\begin{bmatrix}1 & T\\\\0 & 1 \\end{bmatrix} x_k + \\begin{bmatrix}T^2/2\\\\T \\end{bmatrix} u_k + w_k$$\n",
"\n",
"$$y_k =\\begin{bmatrix}1 & 0\\end{bmatrix} x_k + \\nu_k.$$\n",
"\n",
"Measurement noise $\\nu_k$ comes from sensor imperfection. We want best (most accurate) estimates of vehicle's state $x$.\n",
"\n",
"It is relatively straightforward in Kalman filtering to specify measurement noise, in this case to specify covariance \"matrix\" $R$ (here a scalar only) of position noise, since generally sensor noise is known. Let us specify the standard deviation of position noise in this example as ± 3 m. A more difficult task is to define the covariance matrix $O$ of system noise. For this purpose we apply the so called Discrete White Noise Acceleration model (DWNA). If acceleration noise $w_{acc}$ is known (we specify ± 0.2 m/$\\mathrm{s}^2$ for this problem) than covariance matrix of system noise is\n",
"\n",
"$$Q = \\begin{bmatrix}¼T^4 & ½T^3\\\\½T^3 & T^2 \\end{bmatrix} w_{acc}^2.$$\n",
"\n",
"The function `kalman(duration,dt)` makes estimation of vehicle's position and speed with Kalman filtering. It makes plot of estimated and measured vehicle positions by simulation. The plot demonstrates well why this process is called filtering."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAGDCAYAAADd8eLzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABza0lEQVR4nO3dd3hUVf7H8feZlknvPQRCDSX0jiBNQEAs2BtFxV5W197ddcW1r22XnwWsqFgAwQKIgoD03ksSSO9t0iYz5/fHDDFIC5BCwvf1PHkyc9ucewP55Jx77jlKa40QQgghmjZDYxdACCGEEGdOAl0IIYRoBiTQhRBCiGZAAl0IIYRoBiTQhRBCiGZAAl0IIYRoBiTQhThNSqkkpdTIs6Ac4UqpZUqpYqXUK0qpx5RS77nXtVJKaaWUqbHLWVt/PZ/GLo8QTUWT+U8uRH1QSiUBN2utF7vfXw28C1yitf6tMct2CqYBOYCfPsnAEkqpX4FPtNbvNUTBTlOtz0cI8SepoQvhppSaBLwNjGtCYQ7QEtjREOGnlDLW47GVUsrAGZxPU2qJEKKuSaALASilpgGvAKO11ivdy9oopX5RSuUqpXKUUp8qpQKOs/8zSqmvlFKfuJuKtyql2iulHlVKZSmlDimlRtXYfopSaqd72wNKqVtrrBuqlEpRSj3g3jddKTXlOJ87E5gEPKSUKlFKjXSX5ZNjbPs8MBh4y73tW+7l8UqpRUqpPKXUbqXUlTWPr5R6Vym1UCllA4Yd47i/KqVeUEqtUUoVKqXmKqWCaqzvr5RaqZQqUEptVkoN/cu+zyulVgClwEfHOB8PpdTrSqk099frSimPv1yrh5VSGcCHDfmzUEp5um9zJLvP/XellOfJzluIeqG1li/5Ome/gCTgayAT6PaXdW2BCwAPIBRYBrz+l31Hul8/A5QDo3HdyvoISAQeB8zALUBijX3HAW0ABZyPK8x6utcNBaqA59z7jnWvDzzOOcwE/lnj/TO4mtUBWgEaMLnf/4rrFsPhbb2BQ8AUd7l74mru7lzj2IXAIFwVAOsxPv9XIBXo4j7e1zU+PxrIdZ+DwX09c4HQGvseBDq7P998jPN5DvgDCHP/HFYC//jLtXrR/XPybMifBa4WnV/d52kEBrrLccLzli/5qo8vqaEL4fpl+wewteZCrfU+rfUirXWF1jobeBXXL/zjWa61/klrXQV8hSt8pmut7cBsoNXhGr7WeoHWer92+Q34GVft+TA78JzW2q61XgiUAB3q5GyPNB5I0lp/qLWu0lpvwBXIl9fYZq7WeoXW2qm1Lj/OcT7WWm/TWtuAJ4Er3c3z1wMLtdYL3fsvAtbhCrrDZmqtt7s/336MY1+H61pkuX8OzwI31FjvBJ52/5zK3Mvq/Wfhvj0wFbhXa52qtXZorVdqrStqed5C1CkJdCHgNqA98J5SSh1eqJQKU0rNVkqlKqWKgE+AkBMcJ7PG6zIgR2vtqPEewMd97AuVUn+4m7kLcP2ir3nsXHcYHVZ6eN861hLo524WLnCX5TogosY2h2pxnJrbJOOqzYa4j3/FX45/HhB5CsePch+z5vGjarzPPsYfGg3xswgBrMD+Y5S5NuctRJ2SQBcCsoARuGpl79RY/gKu5uquWms/XLUudfTup8Z9//dr4GUgXGsdACysi2PXwl87mh0CftNaB9T48tFa336CfY6lRY3XsbhqtTnu43/8l+N7a62nn8Lx03AFZM3jp51i+Y7pDH8WObia9tscY11tzluIOiWBLgSgtU4DhgNjlFKvuRf74mpeLVBKRQMP1tHHWXDdZ80GqpRSFwKjTrxLnckEWtd4/z3QXil1g1LK7P7qo5TqeIrHvV4p1Ukp5YXrfvMcd434E+AipdRopZRRKWV1dzSLOYVjfw48oZQKVUqFAE+5j1sXTvtnobV2Ah8AryqlotznN8D9R0JdnLcQp0QCXQg3rfUhXKF+uVLqBVz3anvi6hS2APimjj6nGLgH+BLIB64F5tXFsWvhDVznl6+U+o+7LKOAq3HVejP4s4PZqfgYV2e2DFzN0PdA9TW9GHgMV2gewvWH0an87vknrvvPW3D1c9jgXnbG6uBn8Xd3mdYCebiunaGOzluIU6K0lnEbhBCnTzWNwWqEaPbkr0UhhBCiGZBAF0IIIZoBaXIXQgghmgGpoQshhBDNgAS6EEII0Qw06ZmJQkJCdKtWrRq7GEIIIUSDWb9+fY7WOvSvy5t0oLdq1Yp169Y1djGEEEKIBqOUSj7WcmlyF0IIIZoBCXQhhBCiGZBAF0IIIZqBJn0P/VjsdjspKSmUlx9v2mZxOqxWKzExMZjN5sYuihBCiGNodoGekpKCr68vrVq1osbU1uIMaK3Jzc0lJSWFuLi4xi6OEEKIY2h2Te7l5eUEBwdLmNchpRTBwcHS6iGEEGexZhfogIR5PZBrKoQQZ7dmGeji9Pn4+DR2EYQQQpwGCfRzQFVVVWMXQQghRD2TQK8HSUlJxMfHc/PNN9OlSxeuu+46Fi9ezKBBg2jXrh1r1qzBZrMxdepU+vTpQ48ePZg7d271voMHD6Znz5707NmTlStXApCens6QIUPo3r07Xbp0Yfny5cCRNeo5c+YwefJkACZPnsz999/PsGHDePjhh9m/fz9jxoyhV69eDB48mF27dgGQmJjIgAED6NOnD08++WQDXiUhhBB1qdn1cq/p2fnb2ZFWVKfH7BTlx9MXdT7pdvv27eOrr75ixowZ9OnTh88++4zff/+defPm8a9//YtOnToxfPhwPvjgAwoKCujbty8jR44kLCyMRYsWYbVa2bt3L9dccw3r1q3js88+Y/To0Tz++OM4HA5KS0tPWoY9e/awePFijEYjI0aM4L///S/t2rVj9erV3HHHHfzyyy/ce++93H777dx44428/fbbdXGJhBBCNIJmHeiNKS4ujoSEBAA6d+7MiBEjUEqRkJBAUlISKSkpzJs3j5dffhlw9c4/ePAgUVFR3HXXXWzatAmj0ciePXsA6NOnD1OnTsVut3PJJZfQvXv3k5bhiiuuwGg0UlJSwsqVK7niiiuq11VUVACwYsUKvv76awBuuOEGHn744bq8DEIIce5yVMGO76DLRGiAjsXNOtBrU5OuLx4eHtWvDQZD9XuDwUBVVRVGo5Gvv/6aDh06HLHfM888Q3h4OJs3b8bpdGK1WgEYMmQIy5YtY8GCBdxwww08+OCD3HjjjUf0Pv/rY2Xe3t4AOJ1OAgIC2LRp0zHLKj3YhRCijlUUw1dTYN8i8I2AVufV+0fKPfRGMnr0aN5880201gBs3LgRgMLCQiIjIzEYDHz88cc4HA4AkpOTCQsL45ZbbuGmm25iw4YNAISHh7Nz506cTifffvvtMT/Lz8+PuLg4vvrqK8A1UMzmzZsBGDRoELNnzwbg008/rb8TFkKIc0VROnxwIez/Bca/3iBhDhLojebJJ5/EbrfTtWtXunTpUt0h7Y477mDWrFn079+fPXv2VNeyf/31V7p3706PHj34+uuvuffeewGYPn0648ePZ/jw4URGRh738z799FPef/99unXrRufOnas74b3xxhu8/fbb9OnTh8LCwno+ayGEaObyDsAHoyE/Ea77EnpPabCPVodriHV+YKU+AMYDWVrrLu5lQcAXQCsgCbhSa53vXvcocBPgAO7RWv90ss/o3bu3/ut86Dt37qRjx451dyKimlxbIYQ4gcwd8PEl4LDD9V/jiOzGbym/MTx2eJ1+jFJqvda691+X12cNfSYw5i/LHgGWaK3bAUvc71FKdQKuBjq793lHKWWsx7IJIYQQp2V3RjHvLT9w5MKUdfDhhaAMMOUHcoJiuXXRrdy79F7WZqxtkHLVW6BrrZcBeX9ZfDEwy/16FnBJjeWztdYVWutEYB/Qt77KJoQQQpyuz1Yn888FO9mXVeJacOBXmDUBPAOomDSPrwq2cvm8y9mQtQl7xuWo8jYNUq6G7uUerrVOB9BapyulwtzLo4E/amyX4l4mhBBCnFWS81zjgCzcms4t4VvZ+MPdbA6LZF9cf9YtuYW88jy6BHfBN/8y8nQw3VoENEi5zpbH1o713NQxb+4rpaYB0wBiY2Prs0xCCCHEUQ7mloKqZOXuvzPz0H7sYUEo7MQUJdI7vDdXx1+Ns7Q1V/3+B89OaI3Z2DD9zxs60DOVUpHu2nkkkOVengK0qLFdDJB2rANorWcAM8DVKa4+CyuEEELU5HBq0mx7adnmv2w3V3KR05exQ5+jR/RAvM3e1dtN+mANIT4WrurT4gRHq1sN/djaPGCS+/UkYG6N5VcrpTyUUnFAO2BNA5dNCCGEOKGk9INERb2FQZVxaUorolp8wHmtLjgizLelFvLbnmymnheH1dxw/bvrLdCVUp8Dq4AOSqkUpdRNwHTgAqXUXuAC93u01tuBL4EdwI/AnVprR32Vrb4VFBTwzjvvNHYxhBBCnIKv16fw4YrE6veZReVHzgdSlM6sBZeQY3Zwb8gl7Ax+lO+3ZfPXx7/f/XU/vlYT1/dv2VBFB+q3l/s1WutIrbVZax2jtX5fa52rtR6htW7n/p5XY/vntdZttNYdtNY/1Fe5GsLxAv3wqG9CCCHOPh/9kcwbS/bidLoC+um527nk7RX8vj+Vnzd/wL8/u4BvPaoIze9I7yFPcGmPGHZlFLNwa0b1MfZnl7BwWzo3DmiJn9XcoOWXkeLqwSOPPML+/fvp3r07ffr0YdiwYVx77bXVE7N06dKletuXX36ZZ555BuC4U5wKIYSofyl5pRSU2tmXXYLTqfkjMRe7yuP+pVfxwKbX+MwKrYztOZh7A5H+nlzZO4aEaH+enreNgtJKAP776348TAamDIpr8PKfLb3c68cPj0DG1ro9ZkQCXDj9hJtMnz6dbdu2sWnTJn799VfGjRvHtm3biIuLIykp6bj7TZs27ZhTnAohhKhftooqcm2uUF6TmIfDqSmqOkR02//hcNp4MsvIRVO/5e8/FmELKMJoUIBi+sQEJry1gofmbOG8diF8uzGV6/u3JMTH48QfWA+ad6CfJfr27Utc3In/WjvRFKdCCCHq16H80urX65LyKC3LJjr2LUyOCp6vbMXtuTfRvSqc5Lx0WgR5VW/bOcqfO4e24T+/7OPnHZmE+Fi4ZUjrxjiFZh7oJ6lJN5TDE6wAmEwmnE5n9fvDU56ebIpTIYQQpyclv5RALwveHsePvEN5ZQDEBHpyIHEfJQXPU2KtYlbkWCIGPEfZC0tZsCWd5NxSesYGHrHv/aM6cEXvFpiNBgK8zA3as70muYdeD3x9fSkuLj7muvDwcLKyssjNzaWiooLvv/8eOPEUp0IIIU7fVf/7g1s+WndUb/SaDrlHf7u5fSn9PZ9gtaeDwc6BdB39MmF+XvRpGcSc9SkUl1cRW6OGfliLIC8i/K2NFuYggV4vgoODGTRoEF26dOHBBx88Yp3ZbOapp56iX79+jB8/nvj4+Op1x5viVAghxOmxVVSRWlDGyv25zFmfctztDuaV0t9zEyuy/s7ngRa8C9tzXsIT1esvTIggtcBVi28Z7H28wzSq5t3k3og+++yz46675557uOeee45aHhcXx48//lifxRJCiHPK4RD2shh5fuFOhsWHHbPDWkjKQpyRH7HJYsE/ezQpOUMZ0Cakev2YLhE8O38HAC2Dj66hnw2khi6EEKLZOtyU/syEztgqqnjrl31HbqA1/PEuEaWvs9PDwjP9n6B9yFW0CPI6ovNbpL8nPWIDAI7Z5H42kBq6EEKIZisl31VDH9YhjIFtQvjjQO6fKyttMP8+SrZ9xSsxsYQY2zK2w5X0jqrAVlF11LHuHNqWZXuzG/U++YlIoAshhGi2DuWV4mEyEOJjoXuLAN78ZS+2iiq8S5LhixtIz93Ni/GDKCk/yEURt6GUItzPesxjjewUzshO4Q18BrUnTe5CCCGarZT8MmICPVFK0T02AKeGQ6u+oWLGMJ4llwtjY/ilPIWKnGH0jerW2MU9I1JDF0II0Wwdyi+tvhfePcqX+01fErR8PlNiYtlqcNA3aDxDIi7n6a/Tjrhn3hRJoAshhGi2UvLLXJ3ZSvNwfDuJqtAtjPNrgTJbsB+6giU7O1He1jUqZ0ygZ+MW9gxJk3sjmzlzJmlpadXvb775Znbs2HHGx01KSjrho3PHM3nyZObMmXPGny+EEPWhoLSSQdN/YW1S3km3LSq3U1hmp4cpiaIZ53OVfT8z/f2xl3fjquhXKC/sRHSAJyv25RLic+KR5JoCCfRG9tdAf++99+jUqdMZH/d0A10IIc5muzOKSS0o49fdWSfdNiWvjCuMv3LJhqm87KnJNZu5qsWLFCRfxdd/VNI1xp9Pbu6Hr4fprH0U7VRIoNeTTz75hL59+9K9e3duvfVWHA4HkydPpkuXLiQkJPDaa68xZ84c1q1bx3XXXUf37t0pKytj6NChrFu3DgAfHx8efvhhevXqxciRI1mzZg1Dhw6ldevWzJs3D3AF9+DBg+nZsyc9e/Zk5cqVgGsK1+XLl9O9e3dee+01HA4HDz74IH369KFr167873//A1xDzN5111106tSJcePGkZV18v8kQgjRULTW/Lo7q3qO8sOPoW1PKzrxjo4qfBf/nZfMM/g1MoFvrYrJnacwtl1/ANIKy7mqTwviQryZfWt/nr80oV7PoyE07faFk3hxzYvsyqvbOcXjg+J5uO/DJ9xm586dfPHFF6xYsQKz2cwdd9zBP//5T1JTU9m2bRsABQUFBAQE8NZbb/Hyyy/Tu3fvo45js9kYOnQoL774IpdeeilPPPEEixYtYseOHUyaNIkJEyYQFhbGokWLsFqt7N27l2uuuYZ169Yxffp0Xn755eqx4mfMmIG/vz9r166loqKCQYMGMWrUKDZu3Mju3bvZunUrmZmZdOrUialTp9bpNRNCiNO1dHcWU2eu44PJvRkeH1498tsJA91eBnOm0uLAQt6qGs/8wFxaGVtxe/fbQZuwGA0YDYoJ3aIA14xpzUGzDvTGsmTJEtavX0+fPn0AKCsrY8yYMRw4cIC7776bcePGMWrUqJMex2KxMGbMGAASEhLw8PDAbDaTkJBQPa+63W7nrrvuYtOmTRiNRvbs2XPMY/38889s2bKl+v54YWEhe/fuZdmyZVxzzTUYjUaioqIYPnx4HVwBIYSoG4t2ZAKwL6vEFejuGnp2cQVZxeWE+f7lmfHyQvj8WkhewcIWD/BunsJg+5K3R7yNh9E15Ou4rpFE+lvxtZob9FzqW7MO9JPVpOuL1ppJkybxwgsvHLH8+eef56effuLtt9/myy+/5IMPPjjhccxmM0opAAwGAx4eHtWvq6pcoxi99tprhIeHs3nzZpxOJ1brsQdE0Frz5ptvMnr06COWL1y4sPozhBDibOJ0apbsdN0GTMyxAZBSUIrFaKDS4WR7WhFhHf78nVdWcBDr7GtQWbtg4nt8vSEWc/CjtAuKZ3D04OrtXruqe4OeR0ORe+j1YMSIEcyZM6f6fnReXh7Jyck4nU4mTpzIP/7xDzZs2ACceKrV2igsLCQyMhKDwcDHH3+Mw+E45nFHjx7Nu+++i91uB2DPnj3YbDaGDBnC7NmzcTgcpKens3Tp0tMuixBC1KWtqYVkFVeg1J+BnppfRv82wQDscDe7b8jcwAOLbmfgd+N4zZ4B13wBCZezt+R3HMYspnWddk5UXJp1Db2xdOrUiX/+85+MGjUKp9OJ2Wzm1Vdf5dJLL8XpdAJU194nT57MbbfdhqenJ6tWrTrlz7rjjjuYOHEiX331FcOGDcPb2zWtX9euXTGZTHTr1o3Jkydz7733kpSURM+ePdFaExoaynfffcell17KL7/8QkJCAu3bt+f888+vuwshhBBnYMnOTAwKhnYIY0daEU6nJq2gnNGdI0jMKWFzahqPLf+A+Qfm4+fUdK5y8KGfF/28PGljyyDfMh9fQzQjYkc09qk0CHWiCd/Pdr1799aHe4QftnPnTjp27NhIJWre5NoKIRrShW8sx9fDxOB2IbyyaA+/PTiU81/6lX9c3Jnv9/7KzqoZGE2l3FRUyk3lCnXtbK5Z+xx55Xk4tJOC0lIuinycFy68rLFPpU4ppdZrrY/qSS1N7kIIIc46aQVl7EwvYkTHMFqFuFoef9+XAzhZV/wRO/QreFQZ+Cwjh9urvPCY8iPWqB78e8i/sdltBFpCsCXdTfeQo58gaq4k0IUQQpx1tqQUADCgTTBx7kBfsS8Hk+9Wfkn/ihFeCSxJ24GlPIR+GQ9x6ew0tqQU0C6wHT9O/JFnes1AV4YS7O3RiGfRsCTQhRBCnHUODyDTMsi7Rg09C4/QRbSxBPPK9oXs1u240/IPLjmvO2kFZVz89gqW7MwkxDOEglJXf6UQ33Mn0Jtlpzit9TnRo7EhNeW+FkKIpudQXim+Hib8PE0opQj19aDYsAyzRw53H8rG2H4M7cfNYKGvLwaD4p6R7ej7/GJW7s9lRMdwcksqAQj2tjTymTScZldDt1qt5ObmSgDVIa01ubm5x33GXQgh6tqh/DJigryqK2cdA/OJCP2WLhUVDE+YAld9jL+/HwaDa72f1UykvycZReUAZJe4ZlAL8ZEaepMVExNDSkoK2dnZjV2UZsVqtRITE9PYxRBCnCNS8ktpFeyNUzv5YN3r7PL8kEqlOb98POrCF465T4SflYxCV6DnllTibTHiaTE2ZLEbVbMLdLPZTFxcXGMXQwghxGnSWnMor4zz2oYyf/N7vLHjQ4aWV5CScRWmPtccd78IfytrEl3TquaUVJxT98+hGTa5CyGEaNpybZWU2R30caznvfWvE293cE2n/7Cxoj/RAZ7H3S/C30pmUTlOpybXVnFO3T8HCXQhhBBnmZS8Um4yLoS9D5NkMjJtwBO0SBiC2ajoEn38mdEi/KxUOTW5tkpyiivPqfvn0Ayb3IUQQjRhVRUELb6fx83fcGlYO1r7RTKi09UYlIFNT43C2+P4sRXh7+q4m1FYTq6tgp4tAxuq1GcFqaELIYRoNJVVTn7fm4PDqaEsHz6ZSOzBb7jPej4HqODmrtMwKFdUnSjMASLdgZ5aUEqerZJQn3OryV1q6EIIIRrNJ38k89z3O7gwpoI3HM9jKTrIR7GPslTPp1NwJ8bGja31sSL8XIG+M70Yp4bgc6zJXWroQgghGs2ve7I53zuZf+bcR2leOukTZjNTZYOxmCf6PYHRUPvHzoJ9PDAZFNvTCoFz6xl0kEAXQgjRgFbuz+GG91dTbndQbnfgmfgz7zmfwd/Pn8tMd3HHjqVkq8VEGIaQEJpwSsc2GhRhvh5sS3XNkx4iTe5CCCFE/VixL4fle3NYsCWdLjk/8LbhFQoDO/Fmt8FkJ31OVoUBZ2lbBsbecFrHj/C3suFgAXDuNblLoAshhGgw2cWuIVnzf3mDDrYZLFad+LhNazYkLWBi6+v56Mc4nA5v2vaMOK3jH+7pDhAqgS6EEELUj+yicv5m+oqbbd/yg+rDs7FeVOZt54XBLzC+9XhyDm3g+y3pxAQefwCZE4nwc+1nNir8PM+tiDu3zlYIIUTjqSzl+swXGWFawhfOoTwbbsJk2MerQ15hZMuRANw3sh15tkq6xwSc1kdE+Ltq5cHeHufcrJsS6EIIIc7YT9szyCwq58YBrQCoclZRUllCgDXAtUHufvjyRoZWbOeliLF87ZWLUadzU/yD1WEO0DbMl89u6X/a5Yjwd9XQg8+xDnEggS6EEOIM5dsq+ftXG3AY8rmwqw9rM9fw5sY3SSlJIdI7ksGe0fx9yyI8DCYuCBpKtuc2Wvq05kL/J7mn7xV1WpbDz6Kfa4+sgQS6EEKIM/Ti4jU4Il/DaM1g+JzpALT0acuo8OshfxFflaxld0QYHePGkZ38Hb0DLuO9i546pWfMa+vwaHFSQxdCCHFOeumnXZzfPoy+cUGntN+K5F3Mz3ocs4eN0ozxXNQtmrEd4/llUSHXbn6OLoYkfu4xkUeKNrE5+TsqC3pzabdp9RLmAGF+HhgUhJ5jU6eCBLoQQpzzbBVVvL10Pyn5ZccN9AOFB7A77LQLbFc9tvr3B77nyeXPoYyK/wz7L/d/lI9HoR9D9/zEkENvUKo8WD/wXUaNupaQzA18tX0pn+9sR5if9ZifURc8TEbeua4nCafZqa4pk0AXQohzjN3h5LJ3VvK3C9oxPD6cQ/mlAGxLLTxq2y3ZW3h387v8nvo7AEHWIKJ9oimuLCapKAlV0YpBvndzfsveXBL+KZN33YNJp/GN8zxesF/DBbYEegE9w3tyKD2cz9lU7/e3x3SJrNfjn60aJdCVUn8DbgY0sBWYAngBXwCtgCTgSq11fmOUTwghmrOckgq2phaybE8Ow+PDOZjrCvQDOTZsFVV4e5jQWvPxjo95Zf0r+Fv8uafHPYR5hbE6fTW55bmEeYVxfuR43voukhGXtoFlL/FY+gtkOf35X6tXeWV/FP3aBbFsTzZaa5RS1YPKnIvN4Q2hwQNdKRUN3AN00lqXKaW+BK4GOgFLtNbTlVKPAI8ADzd0+YQQornLLakEYH92CQAH81yBjrGI19e+R5Uxk8TCRNZnrsfL3h2fwhu5st1Q/D3NXNz24urjfPJHMj3VN1yy5l+Qu5P8luMYs/tiSvb6MrBNEKM6hfPk3O0k5thoHepDdkkFFpMBP6s0DteHxpqcxQR4KqVMuGrmacDFwCz3+lnAJY1TNCGEaN7ybK5AP5BtAyA5twSv6C/wbjud2fvf5peDv1BUWcS0LneTue9KdqRWMuXDNdgqqv48SOZ2ui6/jW88nsFcWQBXf4bntR9RrHxwODVDO4QxpH0oAMv2ZAOuYV9Dfc69AV8aSoMHutY6FXgZOAikA4Va65+BcK11unubdCCsocsmhBDngvxSV6CnFZZRbnewLn8uRr+NmEoGMcjjJX676je+mfAN7TwmAAZuPb81m1MKeWbedshLhK9vQb87iDjbJhaE3oy6ax3Ej8Pbw0SHCD8AhnYIpWWwNy2DvVi2NwdwBXqINLfXm8Zocg/EVRuPAwqAr5RS15/C/tOAaQCxsbH1UUQhhGjWDje5aw0rk3eRor4liB7EeU8iMb2ierv1yXlYTAbuv6A9JcWFtNz+OnrXApQyUNDjDoau6sZjfQeAh0/1PsPjQzEaoHWINwBD2oUyZ30KFVUOckoqiQ6ovx7u57rGaHIfCSRqrbO11nbgG2AgkKmUigRwf8861s5a6xla695a696hoaENVmghhGgu8myVKGMJJr9NvLjhCbTTyJCgW0mIDmBvVjFllQ4A1iXnc0G4DY+lz/L0vquZxjfkxo6BezawIPxWCvGhf+vgI4794Oh45t91XnWz+pD2oZTZHaxPync1uUsNvd40Rs+Eg0B/pZQXUAaMANYBNmASMN39fW4jlE0IIZqskooqXvxhFw+O6YCf1Xzc7Q4Vp+Ld9mWUoZzsMi/K0ycS3zmGYG8PnBp2puXRpXgl92e8zGDDFlhpRLUbw1U7+tA9bAyP+kXx+971RPpbiQ3yOur4Ne+RD2gTjMmgWLo7izxbxTk3pWlDavBA11qvVkrNATYAVcBGYAbgA3yplLoJV+jX7QC/QgjRTNgqqvhxWwaX9Yw+Ijw/37SMb/MeZsP3rRgZN5DhscPpGtL1qE5om2wfYVBOrDn3YLLHkVdcQWyQF238nNxt/Ib4z+/FUpFNWxXE/s730Gb07Zj9ojD+3x8s3ZXFxJ4x/Lwjg1sGtz5pBzcfDxO9WgYyb3MaTo3cQ69HjfLsgNb6aeDpvyyuwFVbF0IIcQLzN6fxyDdb6RDhS5dofwDKq8r5eP+/0U4LVQ7FR9s/4oNtHxDmFUaEVwReZi/Gxo0lxDOEPNYTWnUJkf4JrNiXiyfldEn8gICN7/KAOZ8/7D3Y0/JvPLu7BevGjgFv17jow+PD+OeCnTw0ZwveHiZuH9qmVuUd0j6U1T/lAUgNvR7Jw4BCCNHEpOSXAbAvq6Q60N/d/C759lTK027ivJ6juPfSGH499CvLU5ZTXFlMui2dp1Y+hUJhdITQ3jqOcG8rbRJ/4m7TdwSuLIS2F7Cvy73cNq+Mgt122oR6E+j95yQnQzu4An3ToQIeHN2BAK/aTYByfvtQXvppNyCDytQnCXQhhGhi0gpdgf7LwaX8kLOc3LJc9hbsJdJwPntK25GaX4afxY8JbSYwoc0EALTWLDm4hM92fcaGzb3pHp3IFftfJdi8h42GzoRO/hJi+9MW+CysiBs/WM3IjuFHfG6bUNdjaKWVDqYMalXr8naK9CPY20KurfKcnNa0oUigCyFEE3OwIBvPmJksLdhFTFUMbQPb0ieiD6s39AbKqwO/JqUUI1uOZFhIT774/Sau2buUSs8wbqu8j+K4MXwa2796205Rfqx6dAR/vTuulOLNa3pgUAovS+3jw2BQDG4Xwneb0uQeej2SQBdCiCYmyfEtRp89eJdcyrwbnsRscPVoH/TLLwCkFRwd6DjssH4maukLXGnIZ0fsdfiPfZofX1/LNcE+R21uNh77qeaupzmL2W1D29A+whcfD4md+iJXVgghmpD0kgzKrStxFPYmO7M/TqcBDOBwajKKyrEYDeSUVFJud2A1u+ccT1oB398HOXsoj+rP5QUTuK3XpYwLC6VbjD+D2gaf8DPrQnyEH/HuUeRE/WissdyFEELUktaa0krXOOrvbnwPlKa9x8U4nJqkHNfEKlnF5Ticmq4xrk5yaQVlUF4E3/8NZo6Fqgq4+nO2jviUHboVwd4WjAbF3LvOY3zXqEY7N1F3JNCFEOIssyYxj6fnbkNrDbhmNev3/BL25R3i+6RvsRf2ZGS7jgDszSoGIK2gHIDerYIAKN3xI7wzANZ9CP3vhDtWQfxY8krtAATWsoe6aDok0IUQ4iyzcGs6s1Yls2xvDk6n5v9+P0CZdTXX/3A1WkNlzjAGtQtBKdib6ZoC9fB98wFRBl4y/ZcuS28CizfctAjG/Mv1Gsh1z7QW7COB3tzIPXQhhDjLHA7d95YfACBd/4Rn1EJ8jR24JOJO/rujlDYhPsQGebEv689Av8CwjsE/34fTmMPq6Mn0m/wimI+cDOXw1KlSQ29+JNCFEOIsk1PsmvFs+d4cMm15WEN/wVjeifamB6ko9cTbcgg/TxPtwnxcTe62HAZueohbLYvBpwtTKh8kxL8f/cxHz2yWZ6vE18OExSQNtM2N/ESFEKKRzd+cxtNzt1W/z7VV0L91EJ5mI8lVc1GGCrp5Xcf2tGLSCsqIDPBEKUXbUB865i5Gv92Xjvm/MsvjOrhlKSVBXUgtKK0+ntOp+WrdIQrL7OTZKgmS5vZmSWroQgjRyOZuSmPZ3myemdAZpRQ5JZX0aRVEdEgpPxX9wejYCbSiE79tdw2fGhvsDcUZ3HDwCaJNiyn17sbfTc9SGtiBSSYL0QGebDpUUH38hdvSeXDOFvZmlbgC3VsCvTmSGroQQjSyxJwSKquc5JfaqXI4yS91DZFa4T8Xq8nCQ/3uJcE9ZntSbikjDOvh7X5EZS/nFX09D/q9zB+2CKICPAGICvAkvbAMp1PjcGpeX7wXgM/XHCQlv5RgCfRmSQJdCCEaUZXDycE8V/N4RmE5eaWVaA2FbOG3lKXc0f02Qr1C3YGuudU4nxuTHoXAlqjbV6AG3cOC7dnk2SqJdgd6dKAndocmu6SC77eksS+rhJvPi6O4vIqk3FLpENdMSaALIUQjSi0ow+5wPW+eWVRObkklKDtLc2bQ2r81N3S8AYBAD8073u/xqPlzUiJHwZQfIaQdNw1ujZ/Vdfc00t/VCS46wPX9tUV7ePnn3cRH+PLY2I50bxEAIPfQmykJdCGEaEQHcmzVrzOKyskpqcAStIK8ynQe6/cYZqMZSrJh1gTGOpbymn0ih0a8DRYvAPw9zdwyuDUAMYGuZV2i/Inws/LlukOk5pfx4OgOGAyKqefFAUiTezMlneKEEKIRJdUM9MJyLCYwB66iW3Af+kX2g8zt8NnVYMvi507TeWNDLMsCvY84xrTzW9MyxJveLQMBCPOz8sdjI9DadQ/d5J5o5cIuEUwb0poLOkU03AmKBiOBLoQQDeRQXimVDidtQv+c3Swxx4av1YSHyUhmUTmZjrUYzIVcFX817P4Bvr4ZLD4wZSH9grrySptMYoO9jjiuh8nIhG5Hj8eulMJk/HMSVLPRwGNjO9bfCYpGJYEuhBAN5InvtlFQZmfunYOqlyXm2Ggd4o1Du2ZLS61YgLb7MyZ9Oyx5FiK7wTWfg18U/sDEXjGNdwLirCaBLoQQDeRgXilFZfYjlh3IttGnVSAlFVUkFiaSZdjK6JIgzEuegU6XwCXvVt8vF+JEpFOcEEI0AK01aQVl5Npcc5UDlNsdpBWWERvsQZl1BTleb2LS8GjRFhj6KFwxU8Jc1JoEuhBCNIA8WyUVVU7A9XgaQHJuKVrD9oqP2FL+Pm0debyRkccXQY/B0EdAqRMdUogjSKALIUQDSC8sP+p1Yo4NDOVszl/EeJud/0sr4D/FD3IgbFRjFVM0YRLoQgjRAFLd85WD6/E0cAV6vP8CKpyVXF6uuKLyadY52hEiA7+I0yCBLoQQDSC9RqCnFbpe++6bi1/gCto7DARM+I79OhqAEB+PRimjaNok0IUQogGkF5ZjMRnws5pcNfRV79Ar61/s9LAwsfc9hES3rt42WGro4jTIY2tCCNEAUgvKiPK3YjUb6Zv4Lmz8hJdC2mPEwbgOl+NnMeFpNlJmd0gNXZwWqaELIUQDSC8sJ9Lfk6l8x/iCT1gWP5a1vuX0CrwIfw9/lFJEuCdXkRq6OB0S6EII0QDSCsq4VP/MlQXv871hAM+Qh7MyiMmdbq3eJtzPVTMPlRq6OA3S5C6EEPWsyuEkoeR3Lq94ne2BA3hQBWOq2EJ5+s0kRIVUbxfh56qhB8psaOI0SKALIUQ9y9+zkjdMb7E4uAPPBJVhtG8lpGo8JcbOBHj9Gd79WweTWVSB2SiNp+LUSaALIUR9ytlHwHfXc0gH8kKIN1aDmYy913CoMpb+rX2P2PTqvrFc3Te2kQoqmjr5M1AIIepLSRZ8chlODdebryWnMovr2t+Cs7wFVU5Nu3Cfkx9DiFqSQBdCiPpgL4fPr4aSLBZ0eZ08vwP4mH24uP3o6k3ah/ue4ABCnBoJdCGEqGtaw8IHIHU9XDaDdc4ozH7buDDuQkK8ffC1uu52tpcauqhDEuhCCFHX1n0AGz+BwX+HThPYVrgMDHYubXspAFH+ngC0DZMauqg70ilOCCHq0sHV8MPD0PYCGPYYeWUFJNoX4GWJpktIFwAi/K0Ultnx9zQ3cmFFcyKBLoQQdaU4A768AfxjYOL/4QBu//kBtCmPKe1fRrnnN79vZDvySysbt6yi2ZFAF0KIuuB0wFeToaIYbvgWbQ3glXUvsaNgDSr3Mqb0Hl69aY/YwMYrp2i2JNCFEKIurHgdDq6CS2dgD23PP1Y+zbf7vkUXDmJ8q4l4mIyNXULRzEmgCyHEmUrfAktfoKTjRSzyNPPVD5PYmrOV4RHXMXdnFyaMi2rsEopzgPRyF0KIE3A4Nf/9bT8p+aXH3qCyFL69lf2+wUwghadWPU1RZRHPn/c8ZZkjCfO10i8uuGELLc5JEuhCCHECX6w9xPQfdvHxquRjb/DDg+zO38fU8CC0UswcM5P5l8xndMtxLN+bw+jOERgNqmELLc5JEuhCCHEc+bZK/v3TLgD+OJB71PqqjZ8wa/9cbmgRg9nszfP93qG1TwJKKdYl5VNmdzC0Q2hDF1uco+QeuhBCHMdLP++muLyKMZ0j+HlHBsXldnytZt77fQ84F/D9nhnsDw5kSNQAHu7zBJe9uZ3OUcXMmtqX3/ZkYzEa6N9amttFw5BAF0KIY8gpqeDzNQeZNKAVozqF8+P2DNYl5eNlreD1nTejLPm0UwZe7/c0wztMZPHOLHJKKvltTzb7skr4bXc2feIC8faQX7OiYci/NCGEOIakHBtaw/kdQukRG4jFaGDVgRx+y38FkzmfVzOyGXb5bFSboQB8uzGFAC8zpRUO/v3jLnZnFjOxV3zjnoQ4pzTKPXSlVIBSao5SapdSaqdSaoBSKkgptUgptdf9XUZeEEI0mtSCMgBiAjzxtBjp3iKAH5MWkFa1hrvyC9hUNI7MkP4AFJbZWbwzi0u6RzO+WyQ/78gEYEh7uX8uGk5jdYp7A/hRax0PdAN2Ao8AS7TW7YAl7vdCCNEoUvJdgR4d6JpIpX2LEgqsn9OjrIILjB15y3EJa5LyAFi4NZ3KKieX9YxmysA4AML9POgg06OKBtTgga6U8gOGAO8DaK0rtdYFwMXALPdms4BLGrpsQghxWEp+GUHeFrwsJnLKcvi9aDoBupKHc6qInvwxXh4W1iS6er5/vT6FNqHeJET7kxDjz7iESK7pG1s9drsQDaEx7qG3BrKBD5VS3YD1wL1AuNY6HUBrna6UCjvWzkqpacA0gNjY2IYpsRDinJOSX0pMoCdVziru/eVebJU5zMzMIKXHu3T2C6d3q0DWJObxy65M1iXn89T4TtUB/vZ1PRu59OJc1BhN7iagJ/Cu1roHYOMUmte11jO01r211r1DQ+X+lBCifqQWlBEd4Mni5MVsydnCU9lZdOx+E6PGXg5A37gg9mSW8OR322kd6s31/Vs2conFua4xAj0FSNFar3a/n4Mr4DOVUpEA7u9ZjVA2IYRAa01qfhnRAVY+3PoeLR2asR6RGEY+XV0L7xcXBLiC/8nxnbCYZJwu0bhq9S9QKRWolOqslGqtlDqjf7Va6wzgkFKqg3vRCGAHMA+Y5F42CZh7Jp8jhBCnK6ekkooqJ9p6gB35u7kxvwDjJf8Di1f1NgnRAfh4mBgeH8awDse8QyhEgzruPXSllD9wJ3ANYMF139sKhCul/gDe0VovPc3PvRv4VCllAQ4AU3D9cfGlUuom4CBwxWkeWwghzsjhiVh2535AkMPBhITJ0KLPEdtYTAa+vWMgEf7WRiihEEc7Uae4OcBHwGB3L/RqSqlewA1KqdZa6/dP9UO11puA3sdYNeJUjyWEEHUttaAMX0sim8v3cUeVFevwJ4+5XTt5LE2cRY4b6FrrC06wbj2u3ulCCNHspOSX0T30E3Y5nVw96nUweTR2kYQ4qVo9tqaU6gq0qrm91vqbeiqTEEI0rqRv2eZbwkTfDgS2GtLYpRGiVk4a6EqpD4CuwHbA6V6sAQl0IUTzY8ulsHAGDj8PJg1/qbFLI0St1aaG3l9r3aneSyKEEGeBQ9/fxXc+ZqJ0PC0C2zR2cYSotdo8grZKKSWBLoRo9laseJGrbZupxEzHgCmNXRwhTkltauizcIV6BlABKEBrrbvWa8mEEKKBaK15f91r/Gfvx7RVFrYcuJcu7Ts2drGEOCW1CfQPgBuArfx5D10IIZoFm93G48sfZcmhpVxYVkmbiH+xwW6kZ2xAYxdNiFNSm0A/qLWeV+8lEUKIBqS1ZnXKbl5Y/zDJhYn8PTefC/r9m/MXmrisZxQ9YgMbu4hCnJLa3EPfpZT6TCl1jVLqssNf9V4yIYSohcIyOzvSio5aPm9zGv3/tYSSiqpj7vfl1hXcvOh6MgtTmZGewfWdp3Df1ji8PUw8Plaa20XTU5tA98R173wUcJH7a3x9FkoIIWrr3z/u4or/rsTp1NXLnE7N64v3kFFUzs50V9iXVZWx5OASbHYbh4oO8eqWh7E6THydnIiPuTdjtg5lTVIej4/tSLCPDCQjmp6TNrlrraWrpxDirOR0ahbtyMRW6SDHVkGYr2tc9aW7sziQbQNgV3oRHaMs3LnkTjZkbcDb7I23yRvtsPNpZjJVlnZcln0zEcEG3p/UmxEdwxvzlIQ4bSeanOUJXBOw5B1n/XDAS2v9fX0VTgghTmRbWiFZxRUApBeUVwf6/y0/QKS/FVtFFRtSD/Jj3uNsz93O33v/nV15u1ibtorXMrPxqgog+t55zMw20iM2AA+TsTFPR4gzcqIa+lZgvlKqHNjAn7OttQO6A4uBf9V3AYUQ4ngW78yqfp1WUEa3FgFsSy3kjwN5PDA6ljkH3mNJyTJMZYqn+k5n14GWPDvgAiwbfqCo3MHfA/7BDN8w+sscK6IZOO49dK31XK31IOA2XMO+GoEi4BOgr9b6b1rr7IYpphBCHG3JzkzahfkAkFZYDsCiHZkoBSmG2eSbfsVR3JPvJszlYEprflq+Ev3+KCjJ4gnr4xDUujGLL0SdOmmnOK31Xq31TK31C1rr17XWP2mtyxqicEIIcTzphWVsTytiYq8YrGYDaQWuX0uJOTYiQnP4IWk+/YMvxZZ6GU57EPu3rWGO5RmoKIbJ8/nF1oqoAM/GPQkh6lBterkLIcRZZ4m7uX1kxzCiAjxJL3QF+oGcEnTQXAKtgUzpfAsAmzat5an8R3FgZFbH/1EU3JWSiioi/a2NVn4h6poEuhCiSfppewatgr1oE+pDlL8nqQXlaK1JKluJTe3lrh530T06gpYqk4G/T0UBN+knWVcSQnqBq3k+UmroohmRQBdCNDl5tkpW7s9lbEIkSimiAqykF5SRlJ+HDppLqCWOy9pehmfhfr62/gOzruBvlmeIatuNxBwbae7afJTU0EUzUpv50EOBW4BWNbfXWk+tv2IJIcTx/bw9A4dTMzYhEoBIf0+ySyp4ee2rKFMJkzv8G2PGVvhkIhaD5oryJ+jdvS8+VhO/7c4mJd8V6FJDF81JbcZynwssx/WYmqN+iyOEECe3YGs6LYO96BzlB7hq2kafzSzLmI89fyBjdD58eCV4BTGv0+vs/r2ChzuGkVNcSaXDybqkPAwKwn1lRDjRfNQm0L201g/Xe0mEEKIWDje33zqkNUopduXt4sPkh/GMOYDFGcWYPE9Cv78RIhLg2i8ZXOnLtfYDDGobwpaUQgBW7MslzNeKySh3HUXzUZt/zd8rpcbWe0mEEKIW/trc/uKaFympyqcs7QrG7u/AS8YPUK2HweSF4BtBy2Bv/nVpAh4mI61DvAHIKakgQu6fi2amNoF+L65QL1dKFbu/jp7aSAghGsDqxDzC/TzoHOXH3vy9rMtcxw0db+TmklT+YZrFRq+BcM3n4OFz1L5B3hb8Pc0ARAVIoIvmpTaTs8igiEKIs0ZmUTkxgV4opZi9azYeRg+uSt1BgHk2cx0D2dnp3/QwHfveuFKKuBBvNh0qINJfOsSJ5qVWN5CUUhOUUi+7v2TqVCFEo8kqriDUx4PiymLmH5jPhfgQsPZ9vvW4mPvsd9AyzP+E+7cOdTW7y6Ayork5aaArpabjanbf4f66171MCCEaXHZxBWF+Hszd9QVlVWVcc3AbXPAcCyLvQmOgVbD3CfdvE+pqipdhX0VzU5te7mOB7lprJ4BSahawEXikPgsmhDg3/bIrE6PBwPntQ49aV253UFhmJ8qjiJkb/kPPigo6jXsLul1NVM424M8a+PF0CHfdRWwZ7FX3hReiEdX2mY2AGq9P3J4lhBBn4KWf9vDsvO3HXJddXEEY+Xjuu51Mpbm1x93Q7WoAhnUIY2THcMJO8mz58Pgwvr59IJ2j5FeZaF5qU0N/AdiolFoKKGAI8Gi9lkoIcc7KKion11ZJakEZ0X9pFi9M38fnlme53dNMV982DOhzV/W6YfFhDIsPO+nxDQZFr5aBdV5uIRpbbaZP/RzoD3zj/hqgtZ5d3wUTQpx77A4nubZKAFbszTlyZfZu2s6fyCq/KtJNRm7t+yBKqUYopRBnp+MGulIq3v29JxAJpACHgCj3MiGEqFPZxRXVr3/fVyPQk1fB+6OwaQcvBYQTH9iZwdGDG6GEQpy9TtTkfj8wDXjlGOs0MLxeSiSEaNYe/Gozw+LDqkd6y7NV4mEy4O1hIssd6EHeFlbsy8Hp1Bh2fAPf3g4Bsfw9bDDO8sU81u8RqZ0L8RfHDXSt9TT3ywu11uU11yml5AFOIcQpKy6389X6FDRUB/qkD9bQOcqP6RO7klXk+lUzoVsUs1YeIHf+k5RtfZek6ATMwx5jw7KHMdh60yO8e+OdhBBnqdp0ilsJ/LWJ/VjLhBDihPZmlQCuWvlhSbk2DAZXbTvTXUOf2NmPnutf5oPkJGa3iKaKXFj2AAbtQUTVZQ1fcCGagOMGulIqAogGPJVSPXD1cAfwA+QBTiHEKdubWQxQ3fGtospBcXkVh/JKAcguKidOZRD44zM83KKEQyY/Lm1zGZe0m0B2WTYvL8gi0je80covxNnsRDX00cBkIAZ4tcbyYuCxeiyTEKKZ2pN5uIbuqonn2+zu95UUl9vxS1/B6z7PMMnTjxKTL6WJk8gy9qfbwB4YDYqnChfTI1zmMBfiWE50D30WMEspNVFr/XUDlkkI0UztcdfQ80pcNfScksO92jW2397ihsTnmBATg8MrgI9Hv8/v28089/0O/rVwJ4+N7UhOSSWhJxk4Rohz1Yma3K/XWn8CtFJK3f/X9VrrV4+xmxBCHNc+9z10W6WDcruDPFslFuz80/QBEat+44WABNLMhbx93vO0D2xP+/NgW1ohn65OZvLAVjicmjBf6ZMrxLGcaGCZwwMi+wC+x/gSQohaKyq3k15YXj2Gep6tkpLcFD63/JMrTb+xotUUPvNzEmSIP+IZ88t6xFBud/LNhlSAkw7tKsS56kRN7v9zf3+24YojhGiu9rrvn/ePCyY5t5TSpLWc/+tN5JiKudpyBYVmO1QVMyDwxiOeMe/XOgh/TzOfrzkIIE3uQhxHbaZP/bdSyk8pZVZKLVFK5Silrm+Iwgkhmo/DPdz7tQ5igmEFcfMvJw8jE6Jj2R65mpSqZdgLetI5pOsR+5mNBkbEh5HhfkZdmtyFOLbazLY2SmtdBIzHNfxre+DBei2VEKLZ2ZNZgrdZMezQW/zH8jb5AQncHTMcp7GcLsYHsO19nPL0Kwg/Rg18VOeI6tdSQxfi2GozsIzZ/X0s8LnWOk+GXBRCnExKfikRflZMRle9ISUjgw+sLzM9JYnUoC7Etp3E/qyX8ak8nx4R/VlVtR+AML+ja+BD2ofgYTJgMRrwtBgb9DyEaCpqU0Ofr5TaBfQGliilQoHyk+wjhDiHlVRUMfLV3/i/5YkA6PTNPJ56J4XmXSz08WaTXzHzs17CqH1oZbiM2KA/x6o6Vqc3L4uJCzqF0zrU+6h1QgiX2kyf+ggwAOittbYDNuDi+i6YEOLsl1tSwaDpv7A1pfCI5an5ZZTbnXy/ORVW/w/9fyOx6FJejupEK79WeKQ/STvLJXgW3EiYT8ARgX68JvV/X96VWVP71uv5CNGU1aZTnBm4AfhCKTUHuAnIre+CCSHOfrszi0ktKGPDwfwjlqcWlBJCIQ/kPAk/PMQBv76M9ZhCuiOTW7vdSog1nBD7xRTnxxHkbaGFO9CDvS2Yjcf+teRlMRHgZan3cxKiqapNk/u7QC/gHfdXT/eyM6KUMiqlNiqlvne/D1JKLVJK7XV/DzzTzxBC1K9Md8/ztMKyI1fs+ZkfPR5moGE7K9o/wuTKOzFFraKVXysubHUhQd4W0gvLKamoIsTHg6gAT0wGJR3ehDgDtQn0PlrrSVrrX9xfU4A+dfDZ9wI7a7x/BFiitW4HLHG/F0LUoyqHk9zq4VdPXUZhhfu7u1tNaR4seIDhG+4iR/tzt+9r3JvYg1zv96hQGTzY50GMBiNBPhb2Z7ueSw/ytmA0KFoEeRHhL4+kCXG6ahPoDqVUm8NvlFKtAceZfKhSKgYYB7xXY/HFwCz361nAJWfyGUKIk/vkj2QGvfgLSTm209r/cA09O78Ifn0R3ugGa9/n18CJ3On9Mm0SulDi/yEmn7082udJhsQMAVxN66WVjurXANMvS+DB0R3q4KyEODfV5rG1B4GlSqkDuKZQbQlMOcPPfR14iCOHkA3XWqcDaK3TlVJhx9pRKTUNmAYQGxt7hsUQ4tx2MM/Vee2fC3by3qTep7x/RmE5bVUKz2a9A5lJED8ehj3O298WYPXexg/5/8bkk0MrruOaTpdX7xfk/ee98GAf1+t+rYPP+HyEOJedNNC11kuUUu2ADrgCfZfW+rTb6JRS44EsrfV6pdTQU91faz0DmAHQu3dvfbrlEEJAfqlr1rPFOzNZtiebIe1DT2n/+KwFvG75DzanJ85rZmOIvxCAZNv/KAt5mzaerenvfT839Bx8xH7BNQPdW+6bC1EXThroSikrcAdwHqCB5Uqp/2qtT/dZ9EHABKXUWMAK+CmlPgEylVKR7tp5JJB1mscXQtRSnq2S+AhfyuwOnvt+Bz/fNwSDoRYDRzns8POT3Ff8Lr8Sz98r7mZh9DDCgLzSAmz+HxNgjOCzsZ/hZfY6avegGiEe5CM914WoC7Vpcv8IKAbedL+/BvgYuOJ0PlBr/SjwKIC7hv53rfX1SqmXgEnAdPf3uadzfCFE7eWXVhLpb2V4fBhPzt1OWmEZMYFHB/ARSrLhq8mQ/Dv/NJ/PF9EHcVa9wxMr1tM1vA1r0zehTMVc1uKZY4Y5/NnkbjEa8PWoza8hIcTJ1OZ/Ugetdbca75cqpTbXQ1mmA18qpW4CDnKafzAIIWovz1ZJ2zAfWoW4RmBLyT9JoKdugC9ugNIccka/wWdbP8bb4E9xaQt25G3mj8wlAFRkXUjvgd2Oe5jD982DvC3IUNJC1I3a9HLfqJTqf/iNUqofsKIuPlxr/avWerz7da7WeoTWup37e15dfIYQ4vjybZUEeVmIDvAEXCO8Aazan0uP534mz1b558abPocPxlCiFFWTF/K2vRCjRxaXxNxDedo13NTyfTbcsIHHOs/Bnjek+pjHcriGHizN7ULUmdoEej9gpVIqSSmVBKwCzldKbVVKbanX0gkh6k253YGt0kGgt4Uod/imuAN9TWIe+aV2dme4pjxl2Utkzb+Tx1q0ZkCQoueSKXyd+D72wq5c1H4EHiYD6YXlmAwmcotdv1aiAo7/THmglwWljuztLoQ4M7Vpch9T76UQQjS4wz3cg7wtWM1Gwnw9SMkvBWB91ho8Y7/ku32HaLE/ma92fsrnsS2oMlRyVdtrCfT0Y21yGr/t6UqEv5VIfyvp7sFl0grKCPQy42U5/q8Xo0ER6GUhxEd6uAtRV2rz2FpyQxRECNGwDjenB7rHR48J9CS1oIw9+XvYVPk6Rk8HCzPeZiGgAvwZ3fICormMN+cXMHtaf2wZWfyuDxDi7UGEv5X0AlftPq2grLrGfyL/uLgLLYNP0gFPCFFr0r1UiHNUvs0OuGroGbYMKny/52BxOXcs3orRYeT/UoowWQpY1/ECRg99jt2pXtz2yXq0hrmbUimtcBDuZ8VgUET5e7I60dXtJbWgjFbBJ5/mdFzXyHo9PyHONbW5hy6EaOJmrkhk2kfrjliWV3q4hm7ikeWPcMi5gHLPpTjKCvgsYz8tnJV85/Uct17yCVSFc+/sjXSNCWBYh1B+3p5JWmEZ4X6uJvPIACuZReU4nJrU/NrV0IUQdUsCXYhzwEerkvl5RybpNWZFy3c3ua/LXcL6zPVMCLyGLxKdLD2wl41lg7jW/Do/lbimcVi+N5tyu5P/XN2di7tHk1Vcwfrk/OrJVCL8PalyapbtzcZW6ThhD3chRP2QQBeimduXVcIB9+Qrv+/NqV6eZ6tEGcv539Y3SPCO5qktbxKtcvgwdjqPVd1Cz/YtySquoNzuYFtqEYFeZmKDvBgWH4bJoLA7NOF+rkCPcgf7bR+vJ9zPg/HdpDldiIYmgS5EM7doRyYAPh4mVuz7M9DzSyvxD/uNvPJcHt+zHmdwB8ZWvMD/ZbbHZFAMcE+WkpJfxvb0QrpE+6OUwt/TzIA2rnUR7kCP9HfVyL09THx6c7/q90KIhiOBLkQzVG53UFjm6vS2aEcGXaL9GB4fxu/7ctHaNaeRI38Pnr5LGW4rpXOvaTBlIWmEkFZYTmywF3Hu0eMOZJewO6OYzlH+1ccf3TkCoLrJvV24D9f3j+Xjm/rSNqzmJIpCiIYigS5EM/TYt1sZ9vKvLN2VxcZDBVzQMYLz2oaQU1LBnswS2P0j8XkPUGKEG/vcD6Ofx2r1rH4uvHWINy2CXI+U/bIrC7tD0znKr/r4F3WL4rKe0QxsEwKA2Wjgn5ckHBH6QoiGJY+tCdEMJeXYyLNVMmXmWgAu6BSOv5cZ0BQtfgnnvv/wSUwL/HUMPXrfUb1fTKAnOSUVtAr2JtTHAw+TgZ/dTfZdov8Ma39PM69e2b0hT0kIcRIS6EI0Q5lFFQxuF0JWUQVOrekY6YuqquDekNf5rmA/n7ZNIM1RQHePS46YHCUm0JNNhwqIC/XGYFDEBHqyP9uGj4eJlkEyCIwQZzMJdCGaGa012cUVjO8WyYeTO1Be5aSi8CAvf3MFX4SWYXUE4O/hgTMnjo6h5x2xb3SgqzPb4fvnLYK82J9to1OkX+3mSRdCNBq5hy5EM1NYZqfS4STM14rJaMAn8Wce+2IMXxjLGOfdh+w9z/Bgp0+wJd9KiM+RE6i0D/PFZFC0c3dsa+GeSrVztN9RnyOEOLtIoAvRhGUWlZNbUnHEsqxi1/twbwMsfIhlc6ewyGrirnZX8dSEGVhMFr7ZkAL8OY77YZf0iGbJA+cT6uvqHNciyFVj7yKd3YQ460mTuxBN2G2frMfPambW1L7Vy7KKKgigmMGrb6Us8w/+1bo9rX0imNrvYcxGE/1bB/PLrizg6OlLjQZFyxrjsHeJ8sdkUPRqGdgwJySEOG0S6EI0UVpr9mQU49RgdzgxG10NbqWp25hneQK/rEJe6j2R1Jw1fDDgScxGMwDDOoSybE82AIEnmY98YNsQ1j9xgbuHvBDibCZN7kI0Ubm2SmyVDsrsDralFroW7lrA0OXXYFV25lzwCB/nrOHqDlfTJ6JP9X7DOoRVvw7yOnGgAxLmQjQRUkMXoolKzrVVv15/IIMe21+E1e+S7RXPVabLsO37hP6R/Xmo70NH7NcqxJu4EG8Sc2wnraELIZoOqaEL0UQl55YC0NWcyqiV18Pqd6HfbfwrcjrF4d8T4xPDK0NfwWw4uoZ9QadwfD1M+Fnlb3ohmgv53yxEE3Uwp5hbTfN50DSHokpPnFd9hqHjOLa8fxdOYyH/Ou8d/CzHftzs/gvac23f2CMGlRFCNG0S6EI0ovSSdIoqi3BqJ06cGDDQLrAdJsNJ/mvmHWDCxlvw8djBfyK7MrPKQpvdH9GveCsFpt+IMYwhITThuLtbzUZahXgfd70QoumRQBeiEWTYMnh9w+ssOLDgqHXdQ7vz6tBXCfUKPXrHqkpY+R9Y9jIKExfFtKHUkIujKobSCgef7fwMXRlM34hrGuAshBBnEwl0IRrYogMruP+3ezAZNTd1uYkuIV1QSmHAQLotndc3vM7V31/Nc4OeY1D0INdOlTbY9Bmlv72Bl+0QlfHjubwoB4chl9njZjPlf+m09w9mxsVRnP/vZUR3kIFghDjXSKAL0YA2ZW3i4d/vw2EPwDv7Vm656jK8PY78b9grvBf3/3o/ty2+jQGhPbhD+9Nt6zwOVpXwj8AY9gX1JsBaRmVFGqNDH6JzSGfOa2dnyc4s7hrWFpxehPlaj1MCIURzJb3chWggPyf9zB2L78BZ5Ytf/p1k5Hnz+uI9R23XATPfeiXwUIWZHRnruCHnVy6NCuPS2Fj+8DSTWe6Hw2GmPGMCw1uMBGBM5wgKy+zM3ZQKQJh76FYhxLlDauhC1LPiymKeW/UcPyb9SBu/eDbvuIT7x/VkZ3oxH6xI4rKeMXSM9IPM7fDrdNg5H4vBxA1xg5nYciDf+/owL30FPXQEv6zqha7yw2oMxJ6fT6x7StMh7UPxshiZvfYQAGF+EuhCnGsk0IWoY3aHk2V7shnWIYyDxcncs/QeDhUd4u4ed5OdMpBtzoOMS4hkQrco5m1K5bsly+no+S1snQMevjD4fug7DXwj8AKuBK7scTtvLtnLkqo9tA71Zm1SPgAtg12BbjUbGRYfxoIt6QCESpO7EOccCXQh6ti8TWk88NVmbh5TxMK0NzEZTDzf/y2CTfHcu2kT57cPJdjHA7J2MTPkE3rsXYA2e6DOuw8G3gNeQcc87q6MYloGezGxZwwv/bSbIG8LvtY/B425sEsEC7ak42EyyIAxQpyD5H+9EHXst30pWKM+54vkzXQJ7sLYiIe444M04A9A85+++TDrOUhcRi+DhVmOC7AOeZBrRvQ54XF3ZhTRMcKPcQmRvPTT7urm9sOGdQjDw2Qg3M8qA8YIcQ6SQBeiDjmdTpblv4XZbwsVWRdwUfzdvLRwL11j/PlHjxLabZ6O1+8bwTcKRj6DoceNzJ25i+JNNrp3KuKTP5K5rl9LOkUdOcJbWaWDpBwbE7pF0SrEm/Pbh9L5L9t4e5i4rGc0FXZnQ56yEOIsIYEuRB16de0MHJ6bGRoylU1Z3Xjs2x20MWXxSdBy/BZ97Qry8a9D92vB5Oq4dl2/WB6cs4UL31gOgFPDC5cdOcrbnkzXNKnxEa4Qrzn/eU0vXNa1/k5OCHFWk0AX4jTllFSwbE82l/WMQWvNl7u/5KNd72Iv7Mrfxt7CLvNynMv/w1jjGgz7TXDe/TD4AfDwOeI447tGsWBrOp0i/dhwMJ81iblHfdaujCIAOkb6Nsi5CSGaHgl0IU6D1pq/fbGJ35N2U2kJZXnm9yxLWUag6kLnkm7E/XA9rQ8sxWH1RfW5BwbcDr4RxzyWp8XIzCmuGvd/f9vP9B92kVNSQYjPn4+e7UwvxttipEWg1zGPIYQQEuhCnIa5m9JYnT8bn7aLeX49eBgtPBw5gm4rfyOBhZAVDiOfwdh7KlhrPwxr3zhXD/e1iXlcmBBZvXxnehEdInwxGKSzmxDi2GSkOCFOUZ6tkmeWfIFH6GI8iuN5vKwTyw6mc/3KD/FylrG+69Nw7xY472+nFOYAXaL88TQbWZ2YB0BFlYMXf9zF2qQ8urUIqIezEUI0F1JDF+IUvbV8BY6gT+movJmZuxiLNmDqdgVzTaO4d4WFlcNHgPn0BnaxmAz0bBnAmsQ8bBVVXDVjFdtSi7iydwz3X9C+js9ECNGcSA1diFOQZ8tiUdpDBFLOf1IS2RRyBRc438R58bvMTo8iPsKPqADPM/qMvq2C2ZlRxD2fb2RHWhH/vb4X/7682xGDyAghxF9JoAtRSxVpm7j3i1GUqVLut7Ug4vY1pPZ7ksRKfzanFLA2KY+hHcLO+HP6xgWhNSzZlcWDo+MZ0+XYnemEEKImCXRxznvn1308v2DHcddrp5Nlvz7FVQuvYZPRQcv0/liHzYTAltX3td/9dT9VTs2wDqFnXJ4esQH4eJgY3Tmc285vfcbHE0KcG+QeujjnfbMhlYJSO4+P63TUOntpHk9/exnzq3JpafZkQsCtfLozhv9z90ZvE+qDp9nIzzsy8bWa6NUy8IzLYzUbWXT/EEJ8PGQIVyFErUkNXZzTyu0ODmSXkFNSQVG5/Yh1toMrue+z85lflctNvl359rpVHMrtQftwH4K8LQAYDYou0a7R24a0D8VkrJv/UpH+npjr6FhCiHOD1NDFOW13hmtIVYCkHBvxkV58sWs2C7fNYldZJg4TjLaM478bhzKwh431SXlc2jP6iGMkRAewNimfYXVw/1wIIU6XBLo459gdThRgMhrYme4aUlWZipi960s2rv6GlJJUupZXMNkrmvPPf5an5xsptxdww/urqahy0i8u+IjjDY8PY+HW9Dq5fy6EEKdLAl2cc6bOXIuPh4l3r+/FjvQivCN/wBDwG/NToZ028W5mNucNehgG3Ud5lZPtqT8xvmskfxzIpaKkkn5xR85Xfl67EP54bEQjnY0QQrhIoItzSnphGcv35mAyKApL7WzM2IbBfxm+ZZ15tXAb/ezpqMtmQKcJAGxJKaTKqbmkezT3jGjHuqR8wvxOb9AYIYSoTw0e6EqpFsBHQATgBGZord9QSgUBXwCtgCTgSq11fkOXTzQ/lVVOjAaF0aD4cVsGAFVOzY/b00nUs7EqLz7J2UiEoxA15TuI7V+97/pk1z/Bni0DCfK20D5cZjsTQpydGqMbbRXwgNa6I9AfuFMp1Ql4BFiitW4HLHG/F+KMXf7fldz12Qa01vywNYMO4b5EB3jyzup5KM/93FOQT6S2McXxBLpFvyP2XZ+cR+sQ7+pe7UIIcbZq8EDXWqdrrTe4XxcDO4Fo4GJglnuzWcAlDV020fxkFZezJaWQH7ZlMGtlEmuT8xibEMmA+EqKPGfSsrKKiZVGfu4/k9WVcWQVV5BWUMb65Hy01qxPzq+TZ8uFEKK+Neo9dKVUK6AHsBoI11qngyv0lVLHfAZIKTUNmAYQGxvbQCUVTdUa96xlob4ePDPfNRpcx5bFfLnqCQKxcU+GGcOdPxCUawHWsD+7hBd/2MWW1ELuH9me/FK7BLoQoklotJErlFI+wNfAfVrrotrup7WeobXurbXuHRoqjwmJE1t9IA8vi5H/u7E3RoOibaiFV1ffjqfdxt0Zobzu+SKeQdHEhXgD8NHKZDanFBLs7cEri/YASKALIZqERgl0pZQZV5h/qrX+xr04UykV6V4fCWQ1RtlE87ImMY9eLQPp3iKA167sxkXh75DqKOUfAb0IvHQ2D4zrCUCUvyceJgM/bs8gyt/KT/cNJiHanwg/K21CfRr5LIQQ4uQao5e7At4HdmqtX62xah4wCZju/j63ocsmmpc8WyW7M4uZ0D0KnE6GHHyB6ZV7OM8znAGXfQw1xkk3GBRxId7syijmtqFtCPbx4OvbB1JcbsdgkPHUhRBnv8a4hz4IuAHYqpTa5F72GK4g/1IpdRNwELiiEcommpHD98/7tfSD727n3ZSfsPn78sCY/x0R5od1jPSjoNTOlb1bAGAxGQj28WjQMgshxOlq8EDXWv8OHK/KI8NtiTqzOjEXH5ODnmvuZ2PiT3weFcGV7a+kbWC7Y27/7MWdKa90YDUbG7ikQghx5mSkONEsFZfbWbUrhU+9X6d89wYeb9eZKA9f/tb7/uPu42c142c1N2AphRCi7sj8jKLZOZRXyvXvLOa54qfoWrGB13pOIMVezD8G/QNvs3djF08IIeqFBLpoVrTW3PbeEp4repI+xr0cHP9vvszfzNXxV9Mnok9jF08IIeqNBLpoVg6mpvJCyZN0UYmoK2cxozwZi8HCtK7TGrtoQghRryTQRfNhy8H/y0vpoFLIuPB9DkZ3Y8GBBVzR4QpCPEMau3RCCFGvJNBF81CcCTPH4V2cxD3qISJ7X8yMLTMwGUxM7TK1sUsnhBD1Tnq5i6avKA1mXQRFaTzm9RQVIf3YlL2RefvncUOnG6R2LoQ4J0gNXTRtBQfhwwuhOJOSK77kq9zWdGvhxVMrnyLKJ4o7u9/Z2CUUQogGIYEumpSKKgc3z1rL1pRCyEuED8dCaT7c+B1rdQcADupvSC5K5tmBz+Jl9mrkEgshRMOQQBdNyp6MEhbvzGLDxrWuMK8sgUlzIaY365PyMXslsyRtDld1uIp+kf0au7hCCNFgJNDFWam43M70H3ZRWGo/YvnerGLaqhQu3XwLOCph0vcQ1QOAtcmZ+MR8TaR3JH/r9bfGKLYQQjQaCXRxVnrrl33897f9LNyWfsTyvORtzLb8E4dTw+QFENEFgMoqJ9tKv6DKmMWzg56VEeGEEOccCXRx1knKsfHBikQANiTn/7mi4BCXbrsbDdxmfA7C4gHX6HBPLnsZQ+AyBoaOp39k/0YotRBCNC55bE2cdV74YSdmo4HOUb5sOOgO9OJM+PhSPBwlXKrvYn8FFJSVsKdgB7N3z2bRoUVUFfThxSufatzCCyFEI5FAF2eVLSkF/LQ9kwdHu3qsv/TTbgozk/D/ciK6KJ2rjdeR0eoLvJVm8JcvAWA2mPEvv4hI63gCPD0bs/hCCNFopMld1ClbRRWH8kprvX1llZPf9mSjtQZg/uY0zEbFDQNa0iM2gGiysXw8Hooz2T36PRIjVuNtDKcsbSIXxkzi9aGv8934JaQkDmJo+7D6Oi0hhDjrSaCLOvXKz3uY+O7KWm+/ZGcmkz5Yw+KdWWit+WFbBue1DcHPaqa7dz5fevwDVVaA84ZveSF9McpcyLSOj1NV2IcEr8sZ0XIE65NsAAxpF1pfpyWEEGc9CXRRp9Ym5ZFVXIHd4azV9hlF5QB88Hsi29OKSMkv48IukZCzF69PL8JqqGRq+JVcvO45NuQtwp47lGu6DsbDZOCguyVg+Z4cgrwtdI7yq7fzEkKIs50Euqgz5XYHuzKKAMgvrazVPnk213arDuTy+uI9GA2K0WH58OFYdisHE+M6stVjEb4WX+KNtxHmuBhPi5EWQV4czCvF6dQs25vDeW1DMBhUvZ2bEEKc7STQRZ3ZmV6E3eG6F17wlwFhjifXVomPhwkvi5HFO7O4IiYf/y8uYbvZyDWhPpQbyyhLuZ4HE96mKDuB9mGuWnhskBcH88pYuT+XnJIKLugUXm/nJYQQTYEEuqgzW1IKq18frnmfTG5JBZH+Vi7vFUOCOsCzeY/gNFl5Pq4T/h4BzBr1JR4V3bjuvdXszy6hbbgP4Ar0lLxSvlp/CD+rSQJdCHHOk0AXdWZzSkH16/xaBnqerZJgHwt3t8vnS88XMHkHMG/E/Wwt2Mvfev2NDqGR/HjfEAa1DaHKqekWEwBAiyAviiuqWLg1nYu7R2M1G+vhjIQQoumQ59BFndl8qIDOUX5sTysi/xSa3C/yTyT0uycgIJzia2fz2pJpdAvtxvjW4wFXeP/fjb1JyS8lOsD1nHmLQNd3u0NzRe+Y+jkhIYRoQqSGLupEcbmdAzk2hnZwPTpW205xbYvXcWfqw+AXBZMX8sa+ryioKODRfo9iUEf+84wJ9EIpV8e32GDXtKgdwn1JiPavwzMRQoimSQJd1ImtqYVoDX1aBeFlMdbqHnrV7p94U0+n0LMFTF7IpvJMvtz9JdfGX0vn4M4n3LdlkDf+nmYmDWxVHfJCCHEukyZ3UScOd4jrGhNAoJfl5DX0XQswfjWZXTqa7f3e5zLPAJ5dMo0wrzDu6nHXST/P02Jk3RMjMRvlb1IhhACpoYs6oLXmh63ptA3zIcjbQpC35cSd4rZ/C1/eSHlwJ66tfAzvgDBeWPMC+wr28Xi/x2s99amEuRBC/El+I4oztmp/LptTCpkyqBUAAV5m8o7VKU5rWPN/MGcqRPdm89CZFOHDxqJ5fLXnK6Z2mcqw2GENW3ghhGgmJNDFGXv3t/2E+Hgwsaert3mQt4WCvza5VxTD1zfBwr9D25Fw/ddk2T0weh3gqwPvMDJ2JPf2vLcRSi+EEM2D3EMXZ2RbaiHL9+bw0JgO1c+CB3pZqjvF7c4oZt2a37k2+UlU3n4Y8RQM+hsYDOQUZ+IRPp8wrwj+NfhfR/VqF0IIUXsS6I1gV0YR6QXlDItv+tN9zlqZhI+Hiev7t6xeFuhlobi8CrvDyc4f/8tlif/G6e2P8cZ5EDcYp3ZiANbmLMZoTee+ntPxNMk85kIIcSYk0BvB20v3szYxjz8eG9HYRTkjWmuW783h/A6h+FnN1cuDvM14UEnVt3dhzv6WW8OiCW01gS62ffww/032FexjWIthrCtai6qIZVzrsY14FkII0TxIoDeCrKJyCstqN5La2exAjo2MonIGtQmpXvbFri/4OOkDrgtP4+X0Yr4MC0FXeWLI+Y6fczQdAjuS4D+C1emrKNdFBJVPkufIhRCiDkigN4LskgrK7A4qq5xYTE33vvHK/bkADGwTDECFo4J31r9GVWUxXwUZAF8qc4dQkTWafq39eeO6Dnz8ewFvL93Pt3c+xDMLV2D1kElVhBCiLjTdNGnCsosrANdwqU3Zqv05RPlbaRnsBVWVLJx3E3lVNv5VEYD33lu4OOQVKrLGEurrxfbUcsI8w1m0IxOANQeKKC4JIMTHo5HPQgghmgcJ9AZWbndQXF4FQJH7e1PkdGpW7c9lYOsg1PZv0e8O4KOctbQ3+dLuigVkVLVhf6pr7vKLu0VRUlHF0t1Z7MksAeD3fTnkllQS5G1pzNMQQohmQwK9gR2unQMU1dN99HK7g/tmb+RAdkmdHVNrfcT7nRlFeJam80jmAzBnCqvMBvZZLNzY72GC/HwB2HSoAKXg4u7RALy6aA8AIzuGsyYxj8IyuwS6EELUEQn0evDjtozqpuW/yi6pEegnaXJ3OjW5NbavrY0HC/huU1p1gJ4prTU3z1rHRW/+zo60IgBSV3zODx6PEFSym6qL/sN/YloT6hnKhXEXYjUb8bIYKbM7aBHoRacoP7wsRranFdEh3Jdr+7WgosoJQIiPBLoQQtQFCfR68MrPu7l39kbSCsqOWndkDf3ETe5fb0hhwPRfSM61ndLn78ksBmDh1nQO5ZWe0r7HsmBrOkt2ZbE3q5ir31rMj89PZNT2h0kztcBw2+98aoXtuTt4qM9DWIyugA70cn1vG+aD0aDo4p7i9IJO4fSLC8ZsdPVsD/KWe+hCCFEXJNDrmNaatIIySisdPL9g51Hrjwj0k9TQt6cVUVnlZMayA6dUhj2ZxXhZjBgNivd/Tzylff+qtLKK5xfspGuEJ+tGJ7Pc6yFG2ZewoeVNBN21hEMmE29tfIuhMUMZ3Wp09X6B3q7n0tuEuiZa6RbzZ6B7e5joERsIIE3uQghRRyTQ61hRWRW2SgfRAZ4s2JrO73tzjlh/KvfQk9w186/Wp5BVXF7rMuzJLKZzlB8Xd4/mi7WHTjzz2Un8b+luziv5ka/sd+Gz+GH8I9tgmPojPae8SpC/lcdXPI7RYOTx/o8f8Tz54Rp6m1AfAK7r15K/j2pPgrumPrit69l1aXIXQoi6IYFex1LdzewPjGpPy2Avps5ay6PfbK1enl1SQZC3BaNBnbSGnpRjIyHaH7vDyYcrkmr1+Vpr9mSW0C7cl2lDWlNe5eDln3ef+ok4HbD5Cy5fNZGXzDPw8A+D67+GqT9CbH8A/rPhP2zM2shT/Z8iwjviiN2rAz3MFeitQry5a3g7DAZX6F/XvyWPjY2vDnwhhBBnRgK9jqUXuoI7LsSbz27pz+W9Yvh6Qwp3froBcNXQw3w98LOaTngP3e5wkpJfxpD2IYztEsknq5IprTz5Y25ZxRUUltnpEO5L+3Bfbj4vjk9XH2Tp7qzanYDT6Zqv/J0B8O00SpxmlnR7HW5Z6polzV0LX3hgITO3z+TqDlcz9hhDtx5uSj9eYAd5W5g2pE11wAshhDgzEujHobXm1Z93szuj+JT2O9wRLjrAk+gAT/51aQK3DI5jW2oh5XYH2cUVhPp64OdpPqqGXlpZxa4Mdy/y/DKqnJpWwd5c1jOa4ooqtrt7mJ/I4Q5x7cJdQfrAqA50CPfloTlbqmdAO6byItdc5e8OhK8mA7Cu72uMrfwXYX0uqw7yUnsp//zjnzy8/GG6hXbjwT4PHvNww+PDuLxXjNwjF0KIBiKBfhyH8sr4zy/7uOuzDVRUOWq9X2pBOWajOmIEtC5R/lQ5NXsyi12B7uOBn9V8xD10h1Nzy0frmPDWCmwVVdX3z1uFeBMbCgbPZLanFpz08w//AdIh3PUsuNVs5NWrulFYaueq/606ste7ww57foY5U+Hl9q65yk0WuHQG3LGKBY5+eJhNxEe6jrUlewtXfn8lX+7+kkmdJvH+6Pere7X/1ZD2obx8RbdaXzchhBBnRsZyP469WcXu7yW8/cs+7h/VoVb7pRWUEeFvPaIpuXOUqyPYttQiskvKSeJjMn32YCy/t3qbN3/Zy4p9rrHRNxzMJynHhjLl83nii/yWshjvVpW8f2ARvdo/R0JowvHLnVlCsLeF4Bp/UHSK9OPf10Tx5IKlTJj5G//sF8Oogh2Yt38LpTngGQQ9roPu10J0r+r9NiTn0y0mAJNB8d7W93hr41uEeoXy3qj36BvZt1bXQwghRMM46wJdKTUGeAMwAu9prac3Rjn2ZblGWRvZMYx3ft3P2K6RxEf4nXS/tIIyovyPnNu7RZAnvlYTqw7kon1Wsbf8ZzDAIdP/qHIOZV1SIW/8so3endPZmpHEkr2+5Jbn4R33NivSnFza7lJ+32Ek3fk91y68lqExQ7m56810Cz26Brw7o4iEUBMUplJYnMrMvXNYlLmG5IpciIAq4JF98E+nk6iIYHx9OnNpl0lc3H7iEccptzvYnlbEzYNb8cKaF/h81+eMaTWGJwc8iZ/l5NdBCCFEwzqrAl0pZQTeBi4AUoC1Sql5WusdDV2WvVklhPp68NLl3ejz/GLmb05zBXpVJRu3bKClTiWoNBm0E7yCwDMQPIOw5u0nrlVrsJfjBP7I2sAXu7/Aq+VBfs0JxyPiD9r79sGjohtbeY+pP00lKdeGT7tEdjvtWMLg6+wfUdqIWQUxe/z7xPnH8UL+Tj5b2Zrbxuzh86QFXJ/yK/eFDeImcwQUpkJhCrooha/yUzErBz/9z5MXgoPINxoYUFbO9aVlxCsrZea2vFsSSkivaOzGUrZnJ/LEqmfYmLOVh/o8hJfZC4AtKYVUqRK2VL7J1l2/M6nTJO7vfT8GJXdphBDibKT+OkZ3Y1JKDQCe0VqPdr9/FEBr/cKxtu/du7det25dnXy2vaKElOTf0NoJzire/XkjAaqEKzv7sHDdTrqZCznPnIsxPwnlPHFv81KlmOvjzWd+viRZzAQ5IcxhZJ/JQWCVYlZlKGWVVl5xVpAcUY693AejPZIXhoxj2eZCvk/7DZO1gOt0W25sqSF3P6UZu/Eqy6g+/jMhQfzg483kwhL+5vTB4N8Cm2cEs7bZKY6HTxyb6OwVxdPtrqVjWDfwDgW/KArLqxg0/ReGx4fRLsyHVxbtokuXVSQ75qNQxPjG0DagLfmF/mzIW4zFXM59ve7jxk43yrzlQghxFlBKrdda9/7r8rOqhg5EA4dqvE8B+jXEB2dkbWbCqkf+XBDg+jY7BXA/Yq0As18MVU4wOD3pEBqLv4cfiRkltPAMZGRYS2Zt/IO8oIOUYaezJZh/ecYyusJJQU4O2w8WotAEx5goLcngDVsy1qIaPd2/+IZewN8AigB+g+IACG5LVYuBvLLdwOC+fejbqS3Trf4EJM1l5oG57Irsz/PnPc8fe6t4df/HeDk+ZViLYbw69FVMhiN/xP6eZq7rF8uM5QfQGiwmEzkHRzBz6kTWZK5hV+5eduTsIaP0EGYdxezxr9MhqHb9B4QQQjSesy3Qj1UFPKIJQSk1DZgGEBsbW2cfHBTUnultr0EpA7ZKzX9XpnN+Qhv6xbfio7UZ7M7OZMr5wSTnFTF3UwoYS9EhkF1WQkp5LunOvaxJ/h0daKB38FD+1m8qXUO6Vtdqi7KKmfLqMgA2XncBCzan8cy8ray/pxN/f3cOE9qauaR3a0ocJm6bvY0S7cnki4ZzyUBXBzhvp+b/nv4Rm6Elfdt1wgA8GtOb9uHd+ffaf3PRtxfhcJjxjC4iIaQrLw558agwP+ym8+L4aFUyvVoGMjw+jOe+30GIOZ5bE3oy5o1l7M0sARzcfF4bCXMhhGgizrZATwFa1HgfA6TV3EBrPQOYAa4m97r6YG/vUMYNegyA3/fmsL9oNf/o1I+BbULYnrSHVVv3Mi3hQuZvTuPLrM0AdIxrBQo2JCYBDi7saeDHzSU8duE42of6HnH8uBAfPM1GqpxO/D3N+Hma0Bg4UBnIL/YuDG7bCTrH4QNkL/Fjd2YxUZHR1fsbDYr4CD92pBdWL1NKcXn7y+kV3ouZ22eycEsmUdZg3hn5EJ6mIzvm1RTmZ+W3B4cS5G3hQI7r8bjVB/LICq1gT2YJD47uwBW9Ywj1kYlThBCiqTjbejitBdoppeKUUhbgamBeQxfi8CNrbd3DlsYEeOLUkFFYzqG8MpSC89qGsGhHJt9uTOWCTuH4eniweLMB7fAl0t961DGNBkXHSF9CfDwwGBR+VtfkJTvdz41HB/wZwH3jggBoFex1xDE6RfmxPa2IV3/ezZQP11SP0R7nH8ftnR8lO+kirmhzE/4e/ic9xzA/KyajgXZhPgR5W/gjMZd5m9LwMBmYNLAVYb5WuWcuhBBNyFlVQ9daVyml7gJ+wvXY2gda6+0NXY59WSX4e5qra6jRga6wPZRfysG8UsJ9rVyYEMHj324DYPLAVvh7mpmzPgU/qwlfd1j/1bQhbchwDw3r5+kO9HTX6G8xgX+G99Tz4ogJ9CTU98gackK0P5+tPsibS/ehtWt61ZsHtwbgjwOuZ9j7tw4+pXNVStG3VRB/7M+losrJyI7h+HicVf8shBBC1MJZ95tba70QWNiYZdiXVULbMJ/qGurh2nNqfhmH8ktpEeTJiPhwHmcbsUFeDGgdjFNr5qxPISrg+E3dY7r8OYHJ4Rr6Dvdwrof/aADXOPC3nt/mqP0v7RGNxWhgUNsQbv1kPXPWHxnoflYTHSNP/Rnxfq2D+HG7qwf9Rd2iTnl/IYQQje9sa3I/K+zLKqFtjUlFIgNcTeipBWWk5JXSIsiLCH8r04a05qExHTAYFAPbhBDm60FskNfxDnsEP0/X31K7M4rx9TDh73nsWn1NVrORib1iiPC3cnnPaHZlFLM9zXVPfdWBXPq1DsZ4GpOd9Itz1ep9PUwM7RB6yvsLIYRofGddDb2xLdmZSa6tsnpyEwAPk5FwPw+ScmykF5XTwt08/tjYjtXbGA2KT27uh5fFWKvPOVxDL7M7iI/wPcnWR7uoWxT/+H4nX69PpcqhSc4t5cYBrU75OADxEb6E+XowomMYVnPtyi+EEOLsIoFewxdrD/LYt9voEu3H5b1ijlgXHeDJmsQ8tIYWx6mFtw+vfTB7WYwYDQqHUx/RIa62ArwsjOgYxhdrD/LxH0mE+Fi4sEvEyXc8BoNBseCewfha5Z+DEEI0VdLk7rY1pZCHv97KoLYhzJ42gACvI2cRiw70Iq2wHKDWzeonopTCzx2gMYGnHugAV/eNxVbpYFTnCH7+2/knvH9/MqG+HlI7F0KIJkyqZG4JMf68e11PRnYKx2w8+u+cmrXoFkGnH5w1+XmayS+1H9Eh7lSc3z6UNY+PIMz36MfkhBBCnFsk0Gu4MCHyuOsOh67FaCC8jgL08H306IDTr/FLmAshhABpcq+1GHcNPSbQ84i5zs/E4Z7up1tDF0IIIQ6TQK+lw6EbUwf3zw87XEM/3XvoQgghxGES6LV0+B56izoMXz+rGQ+TgWBvy8k3FkIIIU5A7qHXkreHiXuGt2VofFidHfOGAS3pGxckY6YLIYQ4YxLop+D+UXU7lWiXaH+6RJ98IhUhhBDiZKTJXQghhGgGJNCFEEKIZkACXQghhGgGJNCFEEKIZkACXQghhGgGJNCFEEKIZkACXQghhGgGJNCFEEKIZkACXQghhGgGJNCFEEKIZkACXQghhGgGJNCFEEKIZkACXQghhGgGlNa6sctw2pRS2UByHR82BMip42Oea+Qa1g25jmdOrmHdkOt45uryGrbUWof+dWGTDvT6oJRap7Xu3djlaMrkGtYNuY5nTq5h3ZDreOYa4hpKk7sQQgjRDEigCyGEEM2ABPrRZjR2AZoBuYZ1Q67jmZNrWDfkOp65er+Gcg9dCCGEaAakhi6EEEI0AxLobkqpMUqp3UqpfUqpRxq7PE2FUqqFUmqpUmqnUmq7Uupe9/IgpdQipdRe9/fAxi7r2U4pZVRKbVRKfe9+L9fwFCmlApRSc5RSu9z/JgfIdTw1Sqm/uf8vb1NKfa6Usso1PDml1AdKqSyl1LYay4573ZRSj7rzZrdSanRdlEECHdcvUuBt4EKgE3CNUqpT45aqyagCHtBadwT6A3e6r90jwBKtdTtgifu9OLF7gZ013ss1PHVvAD9qreOBbriup1zHWlJKRQP3AL211l0AI3A1cg1rYyYw5i/Ljnnd3L8jrwY6u/d5x51DZ0QC3aUvsE9rfUBrXQnMBi5u5DI1CVrrdK31BvfrYly/QKNxXb9Z7s1mAZc0SgGbCKVUDDAOeK/GYrmGp0Ap5QcMAd4H0FpXaq0LkOt4qkyAp1LKBHgBacg1PCmt9TIg7y+Lj3fdLgZma60rtNaJwD5cOXRGJNBdooFDNd6nuJeJU6CUagX0AFYD4VrrdHCFPhDWiEVrCl4HHgKcNZbJNTw1rYFs4EP3rYv3lFLeyHWsNa11KvAycBBIBwq11j8j1/B0He+61UvmSKC7qGMsk+7/p0Ap5QN8DdyntS5q7PI0JUqp8UCW1np9Y5eliTMBPYF3tdY9ABvSNHxK3Pd4LwbigCjAWyl1feOWqlmql8yRQHdJAVrUeB+Dq5lJ1IJSyowrzD/VWn/jXpyplIp0r48EshqrfE3AIGCCUioJ1+2e4UqpT5BreKpSgBSt9Wr3+zm4Al6uY+2NBBK11tlaazvwDTAQuYan63jXrV4yRwLdZS3QTikVp5Sy4OqsMK+Ry9QkKKUUrnuWO7XWr9ZYNQ+Y5H49CZjb0GVrKrTWj2qtY7TWrXD92/tFa309cg1PidY6AziklOrgXjQC2IFcx1NxEOivlPJy/98egatfjFzD03O86zYPuFop5aGUigPaAWvO9MNkYBk3pdRYXPcxjcAHWuvnG7dETYNS6jxgObCVP+//PobrPvqXQCyuXxJXaK3/2mFE/IVSaijwd631eKVUMHINT4lSqjuujoUW4AAwBVfFRa5jLSmlngWuwvUEy0bgZsAHuYYnpJT6HBiKa1a1TOBp4DuOc92UUo8DU3Fd5/u01j+ccRkk0IUQQoimT5rchRBCiGZAAl0IIYRoBiTQhRBCiGZAAl0IIYRoBiTQhRBCiGZAAl0IIYRoBiTQhThHuKcWvaPG+yil1Jx6+qxLlFJPnWB9glJqZn18thDnKnkOXYhzhHvynO/d02LW92etBCZorXNOsM1iYKrW+mB9l0eIc4HU0IU4d0wH2iilNimlXlJKtVJKbQNQSk1WSn2nlJqvlEpUSt2llLrfPWvZH0qpIPd2bZRSPyql1iulliul4v/6IUqp9kDF4TBXSl2hlNqmlNqslFpWY9P5uIa6FULUAQl0Ic4djwD7tdbdtdYPHmN9F+BaXPMyPw+UumctWwXc6N5mBnC31roX8HfgnWMcZxCwocb7p4DRWutuwIQay9cBg8/gfIQQNZgauwBCiLPGUq11MVCslCrEVYMG1zj9Xd1T5A4EvnLN2wGAxzGOE4lrXvLDVgAzlVJf4pq967AsXFN0CiHqgAS6EOKwihqvnTXeO3H9rjAABVrr7ic5Thngf/iN1vo2pVQ/YBywSSnVXWudC1jd2woh6oA0uQtx7igGfE93Z611EZColLoCXFPnKqW6HWPTnUDbw2+UUm201qu11k8BOfw5D3R7YNvplkcIcSQJdCHOEe5a8Qp3B7WXTvMw1wE3KaU2A9uBi4+xzTKgh/qzXf4lpdRWdwe8ZcBm9/JhwILTLIcQ4i/ksTUhRJ1TSr0BzNdaLz7Oeg/gN+A8rXVVgxZOiGZKauhCiPrwL8DrBOtjgUckzIWoO1JDF0IIIZoBqaELIYQQzYAEuhBCCNEMSKALIYQQzYAEuhBCCNEMSKALIYQQzcD/A6Vqh0i8XRzdAAAAAElFTkSuQmCC",
"text/plain": [
""
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"\n",
"def kalman(duration, dt, accelnoise = 0.2):\n",
" ## function kalman(duration, dt) - Kalman filter simulation\n",
" ## duration = length of simulation (seconds)\n",
" ## accelnoise = acceleration noise (m/sec**2)\n",
" ## dt = step size (seconds)\n",
"\n",
" measnoise = 3 ## position measurement noise (m)\n",
"\n",
" F = np.array([[1, dt],[0, 1]]) ## state transition matrix\n",
" H = np.array([[1, 0]]) ## measurement matrix\n",
" x = np.array([[0, 0]]).T ## initial state vector\n",
" xhat = x ## initial state estimate\n",
"\n",
" Q = accelnoise**2 * np.array([[dt**4/4, dt**3/2],[dt**3/2, dt**2]]) ## process noise covariance\n",
" P = Q ## initial estimation covariance\n",
" R = measnoise**2 ## measurement error covariance\n",
"\n",
" dimt = int(duration/dt) ## number of epochs\n",
" pos = np.zeros(dimt) ## true position array\n",
" poshat = np.zeros(dimt) ## estimated position array\n",
" posmeas = np.zeros(dimt) ## measured position array\n",
"\n",
" k = 0\n",
" for t in np.arange(0,duration,dt): \n",
" ## Simulate process\n",
" ProcessNoise = accelnoise * np.array([[(dt**2/2)*np.random.normal(), dt*np.random.normal()]]).T\n",
" x = np.dot(F,x) + ProcessNoise\n",
" ## Simulate measurement z at epoch t\n",
" MeasNoise = measnoise * np.random.normal() \n",
" z = np.dot(H,x) + MeasNoise \n",
" ## Innovation: z - H.xh\n",
" Inn = z - np.dot(H,xhat) \n",
" ## Covariance of Innovation: S = H.P.H' + R\n",
" S = np.dot(np.dot(H,P),H.T) + R \n",
" ## Gain matrix K = F.P.H'.inv(S)\n",
" K = np.dot(np.dot(np.dot(F,P),H.T),np.linalg.inv(S))\n",
" ## State estimate: xh = F.xh + K.Inn \n",
" xhat = np.dot(F,xhat) + np.dot(K,Inn) \n",
" ## Covariance of prediction error: P = F.P.F' + Q - K.H.P.F'\n",
" P = np.dot(np.dot(F,P),F.T) + Q - np.dot(np.dot(np.dot(K,H),P),F.T) \n",
" ## Save some parameters in vectors for plotting later\n",
" pos[k] = x[0] \n",
" posmeas[k] = z \n",
" poshat[k] = xhat[0]\n",
" k = k + 1\n",
"\n",
" return pos,posmeas,poshat \n",
"\n",
"## Simulate\n",
"\n",
"duration = 100\n",
"dt = 0.5\n",
"p,pm,phat = kalman(duration,dt)\n",
"\n",
"## Plot the results\n",
"t = np.arange(0,duration,dt)\n",
"plt.figure(figsize=(8,6))\n",
"plt.plot(t,pm,label='measured')\n",
"plt.plot(t,p,label='true')\n",
"plt.plot(t,phat,label='estimated')\n",
"plt.legend()\n",
"plt.xlabel('time (s)')\n",
"plt.ylabel('position (m)')\n",
"plt.title('Kalman filter performance')\n",
"plt.show() "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This simple example gave a good demonstration of basic Kalman filtering. In practice, however, often we have a nonlinear dependence of measurements on the system state, that cannot be characterized with a measurement matrix $\\mathbf{H}$. This will be discussed further and will lead to the Extended Kalman Filter (EKF) and Unscented Kalman Filter (UKF)."
]
}
],
"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.9.16"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
![](\"REK.jpg\")
\n",
"* First application: unmanned Lunar mission (1963)\n",
"* Applications in navigation and engineering\n",
"* Sensor fusion (radar, laser scanner, camera, INS,…)\n",
"* Integral part of GNSS receivers, smartphones, computers, games, etc.\n",
"* Forecasting stock prices, time series analysis, approximation of functions\n",
"\n",
"The original [publication](https://www.cs.unc.edu/~welch/kalman/media/pdf/Kalman1960.pdf) received more than 30 thousand citations. Let's have a look at it:\n",
"![](\"KalmanP.png\")
\n",
"\n",
"## Do you understand?\n",
"\n",
"Kalman filter filter yields the exact conditional probability estimate in the special case that all errors are Gaussian. The underlying model is similar to a hidden Markov model except that the state space of the latent variables is continuous and all latent and observed variables have Gaussian distributions.\n",
"\n",
"The above definition is correct but does not make it understand for most of us **what** Kalman filter really is. To understand it, we will consider its properties one after the other:\n",
"\n",
"* recursive\n",
"* optimal\n",
"* makes estimation of indirect (hidden) state\n",
"\n",
"Finally, Kalman filter equations will be derived which may be easier to grasp from geodetic perspective and discuss a simple navigational example to show how Kalman filter works."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Recursive estimation\n",
"\n",
"A very simple example is measurement of the state $x$ (which is a scalar variable) of a system at the epochs $t_1, t_2, ... , t_n$. \n",
"\n",
"Our measrurements are $x_1, x_2, ... , x_n$.\n",
"\n",
"The best state estimate (assuming Gaussian distribution) in this case is the well-known arithmetic mean or *average*:\n",
"\n",
"$$ \\mu_n = \\frac{1}{n} \\sum_{i=1}^n x_i.$$\n",
"\n",
"Suppose we have a new measurement $x_{n+1}$, and our goal is to calculate a new state estimate by using this new measurement. This is the average based on $n+1$ measurements:\n",
"\n",
"$$ \\mu_{n+1} = \\frac{1}{n+1} \\sum_{i=1}^{n+1} x_i.$$\n",
"\n",
"A much more efficient calculation is to \"update\" our previous average $\\mu_n$ by using our new measurement:\n",
"\n",
"$$ \\mu_{n+1} = \\frac{n}{n+1} \\left( \\frac{1}{n} \\sum_{i=1}^n x_i + \\frac{1}{n} x_{n+1} \\right) = \\frac{n}{n+1}\\mu_n + \\frac{1}{n+1} x_{n+1}.$$\n",
"\n",
"Our calculator's `STAT` function works similarly: it \"updates\" the average at each step by using the input data at each step. This is a *recursive estimation*: estimation is updated after getting more data. In our example recursive state estimation can be written in a slightly different form:\n",
"\n",
"$$\\mu_{n+1} = \\frac{n}{n+1}\\mu_n + \\frac{1}{n} x_{n+1} = \\mu_n + K(x_{n+1} - \\mu_n),$$\n",
"\n",
"where \n",
"\n",
"$$K = \\frac{1}{n+1}$$\n",
"\n",
"denotes *gain*, which tells us the correction to the previous average $\\mu_n$, and\n",
"\n",
"$$x_{n+1} - \\mu_n,$$\n",
"\n",
"is *innovation* (i.e. new information), since it decides whether our previous estimate changes at all. For when the new measurement $x_{n+1}$ and the previous average $\\mu_n$ are the same, our estimate does not change at all.\n",
"\n",
"Let us have a numerical example:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAWoAAAEICAYAAAB25L6yAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAABCAklEQVR4nO2deXxU1dn4v0+SCQlbwhIQAkhcEdkNClKVRUVFFLX21R+18qr1VastXaz42iq1VrHSurS1vNa9Vaq4IFVbFHHFNWxu4IKsAdnDEhKynd8f906YDLPcWe69Zybn+/nkk5k7d+597plznvOc5zznOaKUwmAwGAz6kuO3AAaDwWCIjVHUBoPBoDlGURsMBoPmGEVtMBgMmmMUtcFgMGiOUdQGg8GgOUZRGxwjIieJyBd+yxELEdkrIod5dK9ZIvJrL+4Vdt+rRWSz/axdvL6/wXvExFHrg4isAboDjcBe4D/AtUqpvX7KpSsi8gbwD6XUgx7cawpwhVLqO27fK44cAWA3MEIptdxPWQzeYSxq/ZiolGoPDAGGAjd6eXMRyfXyfoaE6Q4UAJ/5LYjBO4yi1hSl1LfAfCyFDYCIjBCRd0WkSkSWi8jokM86i8gjIrJRRHaKyFz7+BQReSf02iKiROQI+/WjIvJXEXlZRKqBMSJyloh8LiJ7RKRSRH5hnztaRDbYr6eJyDNh171XRO6zXxeJyEMissm+xm3ROgERybGvt0pEtovI0yLS2f6sQET+YR+vEpGPRKS7iPwOOAn4s+0C+HOUZ7tfRP5tn7NIRA4RkXvsMlopIkND5AjKsMd+/vPs48cAs4CR9nWqQq5/W8j3fygiX4vIDhGZJyI9w8r8KhH5yr73X0REopRHG1vGjfbfPfaxo4Cg66lKRBZG+G5f+17/LSLr7XtdJSLDReRjuwz/HPady0RkhX3ufBE5NOw3XS8iu0VksYicFPLZdPu3etwus89EpDzSMxlSRCll/jT5A9YAp9qvewGfAPfa70uB7cBZWB3safb7Evvzl4CngE5AADjFPj4FeCfsPgo4wn79KLALGGVftwDYBJxkf94JGGa/Hg1ssF8fCuwDOtrvc+3vjbDfzwX+D2gHdAM+BP4nynNPBd63n7mN/b3Z9mf/A/wLaGvf47iQe76B5Y6I9Wzb7O8UAAuB1cAP7GvdBrwe8t0LgZ52OfwXUA30iFGOjwK32a/H2vcaZj/Dn4C3wuR6ESgG+gBbgTOilMetdnl0A0qAd4Hf2p/1ta+VF+W7wc9n2c98OlBr/x7dsOrRFg7Uj0nA18AxQB7wK+DdkOt9H+hif/Zz4FugwP5sun3ts+zyvAN43+92lI1/vgtg/kJ+DEtR7wX22I3tNaDY/uwG4O9h588HLgV6AE1ApwjXjKRgwpXZ42Gfr8NSkB3Djo/GVtT2+3eAH9ivTwNW2a+7A/uBwpBzLw5VimHXXQGMC3nfA6i3lcNltqIaFOF7bxBfUf8t5LPrgBUh7wcCVTF+j2XAuTHK8VEOKOqHgN+HfNbefoa+IXJ9J+Tzp4FpUe67Cjgr5P14YI39ui/OFHVpyLHtwH+FvH8WmGq//jdwechnOVgd8KFRrr8TGGy/ng4sCPmsP1DjdzvKxj/j+tCPSUqpDlhKsR/Q1T5+KHChPXStsoff38FSar2BHUqpnUnec33Y+wuwrKS1IvKmiIyM8r0nsRQwwP+z3wdlDQCbQmT9PyyLLhKHAs+HnLsCa0K1O/B3rA7pn7Yb4PdiTag5ZXPI65oI79sH34jID0RkWYgcAzhQ/vHoCawNvlHWBPB2LAs2yLchr/eF3jvWtezXPaOcGw2nz30ocG/IM+8AJCi3iPzcdovssj8vomWZhD9TgYjkJSirIQ5GUWuKUupNLIttpn1oPZZFXRzy104pNcP+rLOIFEe4VDWW2wAAETkk0u3C7v2RUupcLMU6F8v6i8QcYLSI9ALO44CiXo9lUXcNkbWjUurYKNdZD5wZ9mwFSqlKpVS9Uuo3Sqn+wInA2Viui4PkTgXbL/s34Fqgi1KqGPgUS2k5uddGLKUXvF47LJdBZRLitLgWlqtkYxLXccJ6LJdUaNkXKqXetf3RNwDfwxqtFWO5ySL61g3uYRS13twDnCYiQ4B/ABNFZLyI5NqTbKNFpJdSahPWEPZ+EekkIgEROdm+xnLgWBEZIiIFWMPVqIhIvohMFpEipVQ9VihYY6RzlVJbsdwPjwCrlVIr7OObgFeAP4hIR7EmCw8XkVOi3HYW8LvgJJaIlIjIufbrMSIyUKyJyN1Y7oSgPJuBdMVMt8NSxlvt+/43lkUdZDPQS0Tyo3z/SeC/7XJuA9wOfKCUWpOELLOBX9nl0BW4Gev3d4NZwI0iciw0TwJfaH/WAWjAKpM8EbkZ6OiSHIYYGEWtMbYifBz4tVJqPXAu8L9YDWc9cD0HfsNLsJTYSqzJoqn2Nb7EmpxaAHyF5VeOxyXAGhHZDVyFNaEUjSeBUzlgTQf5AZAPfI7l13wGy00TiXuBecArIrIHayLtBPuzQ+zv7sZyibzJAaV1L/BdO1rhPgfPFRWl1OfAH4D3sJTyQGBRyCkLsULivhWRbRG+/xrwayz/7ybgcOCiJMW5DagAPsaaUF5iH0s7SqnngTuxXEu7sUYRZ9ofz8cyAL7Ecr/UcrCbzOABZsGLwWAwaI6xqA0Gg0FzjKI2GAwGzTGK2mAwGDTHKGqDwWDQHFcC07t27ar69u3rxqUNBoMhK1m8ePE2pVRJpM9cUdR9+/aloqLCjUsbDAZDViIia6N9ZlwfBoPBoDlGURsMBoPmGEVtMBgMmmMUtcFgMGiOUdQGg8GgORmVN3bu0krumv8FG6tq6FlcyPXjj2bS0NL4XzQYDIYMJmMU9dylldz43CfU1FsZLiurarjxuU8AjLI2GAxZTca4Pu6a/0Wzkg5SU9/IXfO/iPINg8FgyA4yRlFvrKpJ6LjBYDBkCxmjqHsWFyZ03GAwGLKFjFHU148/msJAbotjhYFcrh9/tE8SGQwGgzc4UtQi8lMR+UxEPhWR2fbee54yaWgpd5w/kNLiQmuL5OJC7jh/oJlINBgMWU/cqA8RKQV+DPRXStWIyNNYe8E96rJsBzFpaKlRzAaDodXh1PWRBxSKSB7QFve2rjcYDAZDGHEVtVKqEpgJrMPaXXmXUuqV8PNE5EoRqRCRiq1bt6ZfUoPBYGilxFXUItIJOBcoA3oC7UTk++HnKaUeUEqVK6XKS0oi5r42GAwGQxI4cX2cCqxWSm1VStUDzwEnuiuWwWAwGII4UdTrgBEi0lZEBBgHrHBXLIPBYDAEceKj/gB4BlgCfGJ/5wGX5TIYDAaDjaOkTEqpW4BbXJbFYDAYDBHImOx5BkNrwqT0NYRiFLXBoBkmpa8hnIzJ9WEwtBZMSl9DOEZRGwyaYVL6GsIxrg+DIU2ky6/cs7iQyghK2aT0bb0Yi9pgSANBv3JlVQ2KA37luUsrE76WSelrCMcoaoMhDaTTr2xS+hrCMa4PgyENpNuvbFL6GkIxFrXBkAbMVnEGNzGK2mBIA8avbHAT4/owGNJA0E1hVhMa3MAoaoMhTRi/ssEtjKI2pIzJS2EwuItR1IaUMHkpDAb3MZOJhpQweSkMBvcxFnUrI91uCpOXwmBwH6OoWxG/mvsJT7y/DmW/T4ebwuSlMBjcx7g+Wglzl1a2UNJBUnVTmPhhg8F9jEXdSrhr/hcHKekgqbgpTPxw9mOievzHKOpWQixlnKqbwsQPZy8mqkcPjOujlRBNGQsYN4UhKiaqRw/iKmoROVpEloX87RaRqR7IZkgjkXzJAkwe0cdYRoaomKgePYjr+lBKfQEMARCRXKASeN5dsQzpxviSMw8dfMMmqkcPEvVRjwNWKaXWuiGMwV1aiy9ZBwWXKrr4hq8ff3QLOcBE9fhBoj7qi4DZkT4QkStFpEJEKrZu3Zq6ZAZDEqRzSyw/0cU3bHab0QNRKlrQVtiJIvnARuBYpdTmWOeWl5erioqKNIhnMCTGqBkLIw7VS4sLWTRtrA8SJUfZtJcihlMKsHrGBK/FMXiAiCxWSpVH+iwR18eZwJJ4Stpg8JNsmfxK1TecDe4fwwEScX1cTBS3h8GgC9myJVYqKz6zxf1jOIAjRS0ibYHTgOfcFcdgSI1sWdKeim9YF/+2IX04cn0opfYBXVyWxWBImWwKQ0w2Sidb3D+GA5gl5Iaso7WEIUYjFf+28W3riVHUhlZPtimnZGOf0x27nW3l6icm14ehVZONE2/J+rfT6dvOxnL1E2NRG1o1sZRTNMWWCZZiMu6fdPq2kylXQ3SMovaYTGjkrYlElZMuS7vdIJ15PVrbhKbb7dq4PjzEDAf1I9G462wOfUtnaGMi5Tp3aSWjZiykbNpLjJqxMOPagxft2ihqD8nmRp6pJKqcstlSTGdeD6flmg3Gixft2rg+bLxwSSQzzM4GN4nOz5Fo3LWuaT/TVcbpCm10Wq7Z4Mv2ovPWXlF70ci98jsm0sizxReaCc+RiHLSMe2nrmXspFz9HqGkQ7940Xlr7frwaljklUsikWG2zm6SRHyKOj9HMuiY9lOXMk7G1+xnbpZ06Rcv0hZobVF7NSzyqldPZJidiExeuhYStd78tpjcQLeVjzqUcbJWfbpGKMm0gXTpFy/SFmitqL2qgF76HZ02cqcyeT3sTbRy6+rTzSZilbFXnXiySi8dSi7ZNpBO/eJ2562168OrYZGOGdecyuT1sDfRyq1j2WYb0cp4TL8SzyIqUlF6k4aWsmjaWFbPmMCiaWMTVnjJtoFMSomrtaL2qpHr6Hd0KpPXw95olVhBRL+kjmWbbUQr49dXbvWsE/dT6SXbBjLJiHC8FVcipHMrLp1Du3TA662nwoeZ4RQGcn1RxKaeHIyX23lFqhde1YVU2kBovSkqDCACVfvqfalD6dqKyxfc9v1kegP3Olws1KcYqXH4EQOra3ia33g99wL+5AFPpQ0E9YvudUh7i9pN/LQC0olfnY0uG7Bmy4a26SZb6rcTUm0DOtShjLao3SQbVkWBf+FiukR0eOGnz8SRlxMrNxOfKxLR2oDT59MhxDEWrVpR6/7j6I4uq/Tc7jB0HxbHIlYnnsnP5YREnk8XoyMa2kZ9eJFRK5PCc3REl4gOt2fvdVn5l26y9bmCJPJ8ukeAaGlRe9XT62IRZjI6rNJzeyIrW0de2fpcQRJ5vlh1SAf3kJaK2ivfcTbtWN3acbPD0H1YnCzZ+lxBEn2+SHVIF/eQI9eHiBSLyDMislJEVojISDeF8rKnT3VVlCE2mZ4UHvQfFidLtj5XkHQ8ny7uIacW9b3Af5RS3xWRfKCtizJlfU/fWtDFGkmVbB15ZeNzhbspLjiulNdXbk36+XRxD8VV1CLSETgZmAKglKoD6twUaky/Ep54f12LGN1s6ukziVT8c9kS/gjuulb89IHqMMeQLiIZBs8urkxpglsXo9GJ6+MwYCvwiIgsFZEHRaRd+EkicqWIVIhIxdatW5MWaO7SSp5dXNlCSQtwwXHZU6EyhVTz9epijehMNmxFpQtuuCl0cQ85UdR5wDDgr0qpoUA1MC38JKXUA0qpcqVUeUlJSdICRSpsBby+Mnnlb0iOVCt+ogmcWiO6+ECzATcMA11CUJ34qDcAG5RSH9jvnyGCok4XxgpriZ/D4lR+i7lLK6ne3xD180z1V6eb1rqPphtEc1MUFQZSuq4O7qG4FrVS6ltgvYgEbf1xwOduCZTMIpREIwsyJRLB72FxsguCgnJX1dTHPM9YjomVsd/1QXeuH380gRw56Hh1XUPGl5HTlYnXAU+IyMfAEOB2twRK1CeUaOWNdb5uCtzvYXGy/rlIckejtY6UgkQqY8Gql+F10O/6oDuThpbSvuBgJ0F9o8r4MnIUnqeUWgZEzOqUbhINGUo0siDa+dPnfcb+hiatQsn8dgMlG76ViHytPeQyPG2sQPNEengd9Ls+ZAJV+yKP4jK9jLRcmZiITyjRyhvteKRhut+hZDqEBiXjn4smd6gSAhNyGSRYxpFSbYbWQR3qg05E8tdnaxlpm5TJKYn6URP9wfzsiXUJDXJK0HUUtAxDKQzkMnlEH99nz3UmnnGRafXBTaK5MMf0K/GsjLx0lWppUSdCoomVop1fEMhhZ4Rhk589cSatHAtfbKA4YEGXaiy3TsSzBjOpPrhNNBfm6yu3csf5A10vI69X3Wa8ok608kY7H9Ayk54OoUFOiBb/rusuKzqGuTkxOjKlPrhNrNGHF2Xk9arbjFfUkHjljXV+JAU+asZCrRq0jmTSRJcb1lA6FL+xmJ3jty/a6/qeFYo6XYQr8GxJKuQFfjecRIgV+ZPM75rOemIsZmf4nUve6/qecZOJXjrwTdyqczJpoitW5E8y9cmveuJ2W9BtXUEofi/t9rq+Z5RFPXdpJc8990++2/QZu3Lb8UTVqa5auJk0nPebTBq2R7OGgKR8jH7UE7dHe06v31oz/3ld3zNKUd81/wsekMc4NrAWgG9VZ/5Tf7xrDvxMGs7rQKYM268ffzRTn1oW8bNklKsf9cTtySwn12/trkEv63tGuT42VtVQKtuY3TCGfaoNJ+SsaD7uBl4Ob3QeZmYbk4aW0qlt5EQ9yShXP9w+blvxTq5vXIPekVGK+rAioViqWae6s6TpCEbYitoty8UrP5hJtuM9t0w8Nm3K1Q9/aboWdKVyfeMa9I6Mcn3cOKoDLISNqjMfNB3DT/Oe5ZBADdePH+LaPbMxJrM1EM93mm4fo9dun3RGPUQqKyfXN65B78goRX1qqbVysL59Tz7c04UcUdwzch8jfIx9TQfGMkkvTn2nmeJTj0S6OppoZXXH+QPjrvDzO0Qu9Bl0aMduklGKml2WK+D+qydC+0Ngxu8ZkbPS+kwpaKyDvDaOLqXTRIixTNJLaxmhpKOjiVVWi6aNjXl9HSJ9dGrHbpJZinp3JSDQoSfk5UOfEfDxUzDgfPj3DVBfA1e+YX0WB7caczK9u9uWSWuwOEIxIxTnpFpWfo9KWkunrI2idqRMdm2A9t0OKOKzZsIjZ8KD40ByQTVCxUMw4uq493OjMSfbu7tpmbQWiyOUWFsyeZUOILQ+FxUGELFyJSdzXzc7Wi9Hc248R2vplLVQ1I6Vye5K6BjyvuQo+MFcePmXcPIv4N0/wRszYOD3oF2XmPdMtII6qWSp9O5OLJNkKnprsThCiTRCCeQI1XUNzXnH3eywwutzaK7zRO/rdkfrlZ/ZredoLW5DLcLzgsrkvJy3GSSrgCjxmLsqoSjsRz1kIFz2bzhiHIz/HdTthcfOhqp1LU4Lj1NOJG+t0+27oq12S0fvHkmGqU8tY+itr8QM40unxaFLrHc8OSKFy7UvyKO+UbU4z62Y33hbkSVyX7djlb0KLXTrOTIpdUEqaGFRB5XGbYGH+UfjqXzccHiL43OXVnLXf1Yyv3YtL+88kvyllZErUvdjYfIz8PSl8MhZcM170KZDxN782cWVXHBcKa+v3BrXQnW6fVc00tG7R2v8O/fVN1smwfPc2PFCFxdKshEdZdNeing9N4bITq65sarG0QjJi6G9F35mt55DhwlNL9BCUQeVSR0B2lDf4niwYQbqd9O+oJYva4t4IpaCOHwMfP8ZeOh0WPAbmDAzZpLxg3IlNzXCtq+g5GgQa5+SRLbvCiddvXusCh1rz8cLjivl2cWVKQ9tdXGhJCuHl0PkWLlEghQVBhx1OMnIrePksZvl7/eEphdo4foIDl/qyCPfVtRBZRJsmD1kOwCbVJf4Q6bex8MJ/wMfPQir33bem1cuhr+NhftPgCf/y5q8JLnKlO5hZDwZqmrqIyqwF5dvSsvQVpdJm2Tl8HKIHOle4fcVwZErIFG5vVjlmowLTDcXhS5uPKdoYVEHlUbDvDYUNNW32Lrpp08t477Anzgj50MANiprkjCughj7a/h6ATx9CSd0/C3v7+580CktlN+2r+Cxc6BNBxh5LVQ8DH85AUZczfPt3ya/5lPaUMe/GkfydOMp7Msr5tu8UnbUHOyOcGNXk0iTPk4IWv2pyqPLpE2ycng5RA6/V6Soj586TAqVqNxuj3x0jGxKFF3ceIngSFGLyBpgD9AINCilytMtyKShpbCoiPO7l3D+hQeUSs+iAsbVLuETdRgvNoxgubL813EVRJv2MHkOPHgafwvMZGTgdvaGeCpa9OZ1+yy/dm4+XPGaNWF5/A/hX1PhrbvoVtSH1WUT+HT9Ns7hLb6X9yYAG5q68lGgHwXsZ2nTEbzUOIIdgUNcsRKCFWj6vM8OcrnE2vMRkkvdGY7fq9CCw/ngxrnJ7Gbu5RA53r2CzxJOpHodruSCVnek67s98nE7sskLnGYG1KFTCZKIRT1GKbXNNUnAWlXYsL/FoV+N7kK7+ft5oeFEHm8cDySgIDofBmffTYenL+GhEVv52Sd9Di54peCFH8GWz62JyGBUSae+cMnzsHMNFPehLCeXMuDf7y1n7ssv07FpB2fnvM/xOSuoUwHODHzE/wZmU922lHYLm+A/++DQE+G4S+GoM5r93akQrOiRKhGQ1tSdke4NSVpETU2wdhGseo3NX37E+m27WFnXjafaTebyM0Y4CkvMto1zE+n4ErEA3R756OICS4V4z6Cjxa2F66OZ3IMV9Zk9qgHY07YvsofEe7d+E6C4DydsmcOi8Jn/vVvgrZnw2XNw6m/gyFNbfi4CnctaHLrtze1U1g8CYE7j6Objx3es4ukxu2i3/gPIbwc5ubDqdZh9EfQ9CU7/LfQcap28YzU01EJJv6QUeDTL5Df/+szZTuorX4LKJVan1H8StD3YLZTIfVuwbwfU7rImZTd8ZCnob96AXetpkjy2N/UClc+FuW9yTu27PPf8WF7bfSXjThkd9ZKZtnGuU9rk5TQ/V6e2AW6ZeGxCUUeRrNi0jHyaGq10DAG73uzdykevz2Xh8lWcnZPLcnU4Vaod+TSSRwMBaUDI5YXFazl3QBcrNDa/rdW+Nn9mtYW8QggUQF4BdDgEuvW3jvtAvM5Ml4nzUJwqagW8IiIK+D+l1APhJ4jIlcCVAH369ElSmjZWBQll+9cA3H3N+dzdqW/i18zJheE/hFd/Dc9dCZJjKeXls+G1W63VjMdNgVE/cXS5aL3xR7uLYeRkGHnNgYON9bD4UXjjDnhgNAy6CEqHwfyboKkeuh5tLX8/9jwryqTO6pTIbxdfEKUOUvK3TDyWG5/7mJz6anJQ7KFty0basB/+c6O1ejPIK7+GoZfA0WdA6XGWjz5RlIIlj8Fbf4BdLePXKSi2RhbjbubUlzrwza4mAPrKJq7Pe4r/l/Mfcl7/NwR+CyN/FLHjilbm8SIrksGLIW+4xQZQW98U9Xzr+RUF1NGRfXSQfexW7aiqqrWUYrsSy8jZ/jWTSvZQPEb4+/uVFO9dxcmFqxjT9hs6vpPD5ve688W2OrbX57M3t4gG8mjfsJN+gc30bVdPh8I2UFAEW1bA/l3Qvjs0NcC+7QwHhgPEyM5QNy8X9a8mBBX9pCD5HaBXOfQYZMm/c43VwfccCp0Pt5R5hx7QrmtyCr2hDlST1Tk0NR3IAyRyUGcmNFEQyOOXpx0G21cR2LWaYqw2WE0h9eRRSC3bq+pi14+6ati9Cboekbi8cRCl4heqiPRUSm0UkW7Aq8B1Sqm3op1fXl6uKioqEpfm7+fB/j1wxYIDx169Gd7/K9z0bfI9cM1OuHeInbhpv+WL3r8b+p8Lo2+Ebsc4vtSoGQsjKoiY1l3tLnjnbnjvfuv+ZSfDMefAZ3MtixMFHXvBnk1WR9L3O3DmndD1KPh8LlQ8AttXWdZ9/T7Y8y1Ub4UjToPRN1iN672/wCdzUDW7EJpoUsLynGOQgd9lyNjvQfUWePFnsGmZ1SmNuQm2fQlv/xFWvnigg+xUBocMgO4DLTkPHXnw8zQ1wuo3Leu5bi8sedyKmOk9whrBtCuxGknPIVByDORYwUVl0146qAl3Yjd3BB7ijNyPoG1XKxa+TQc48nQYMhly86KWuQCTR/Q5OBZ+SE/rOavWWbJ26mtZcIGCsOdogr3fwsalsOJF1u6sZfaa9lQ35rCyqQ8fqaMpDOQlFiWjlNVBx8g346gO1e6Cj5+GFfPYuXoZ7VQ1+RJlIllyINAO6vYc/FlBEfQ+gY17FTs2riJXNdCeGjrJHvJoZDft+LqpJ7tzOjCoZwd65tdY9a5jT6haC7n5/HVZHS9X92OLKqaz7GFwzirasp868qgnj3qVS540UibfEihoxxWTTrfmfQqL4ZBBlnwNtVYunoZaSymv/wDWfQBbV1idQX4H63ffszHs2XKt+lTYyXqWtp2h+FCrLXToAfu2W+1mzyZLSe751rrGvu1YeYEOsdp/Qy3kBKx7FBSxnSLWV9VR0LiHspxv7bDg8NkPi72qgPZSC8BWVURdiH0rCMVtA7TNw7p3++7wi+QW8YjI4mjzf44UddjFpgN7lVIzo52TtKJ+8iJrmfhVbx849s/JlkK59qPErxdK3T6rR93yuWVZHzYaTv9dsxJxSiRrqDCQ66wxV62Db96EwRdBrr3DyO5NsGIerHnHsqob9sOyJy1Fd8hASyEW94HeJ0DVesva7tDDUjrL/2kpSrAq9IALoNOh0KajdfyzubAtpNK0KYLzZkG/s1rKtX+v1WFs+hg2fwLffmK5Z1DWaKNtF2sIW1BkVfiNy6xGHKTLkTDqxzDk+5CTE9XqiKqgigpYNGEHfPM6bP3CamQ7V1sKtrScL2uLmbtiD7k0ki8N5FNPPg1sVx35XPWlhnyOlA0cn7OSJslndNEmOu75uuVNcvOhqLdVLnkFkJNnhV822q62wk7srGmkE7ubv7JBdWW76siWvJ6cNuFCKzIov52VomDLZ1Zn2eVIKOho/W5bPocPH4StK+2OodAasXXoYV1QNUGXw7n+zf2sbepOjigCWIrzMNlEkVRz5UllVj1ZtdCStetRrG43mFdW17OjsZDdtKVaFVKSt4/zBnRiwOF9rBW7NTug57ADnaRqsupNt/6QkxO17Fv8DlGMjUgdbDQEWD1jgsOzsTrL/busupmTA3s2W79LUPnu+Rb2brY6rtoqqN5uKfr66pZ3bd/NKucOPSzl3LGnVQY711rpJAqKLYt3/x5LcVdvtT7Pbwddj4RAW+tSxYeyeN0u5i9eSX1jEx2ooVj2UpVTTJtc6Fy/mVya7LtapdK2TR4TBva02l5JPzhmYlIuzZQUtYi0A3KUUnvs168Ctyql/hPtO0kr6qcvtYZd13544Nj9I60Ge/HsiF/xY3bW9Xvu+AYen2T5+E7/LZRfHrlD2buFxa89w8JP1jC/+ghqio5oKYtSsPlTWLMI2pdAnxOhYw9nMuzfA2/+Ht69z+oESvpZFltegaXwhk6GbscCqoWvPVZHBjjr5JSCFf+CpX+HLSutBttk+d4blVBHgHry6Cj7Woi8QXVFKWF3bjHHTvyxPScgsGMVbKiwFGBBR6ivta5X1NtqXF2Pgj4jKbtpPkXsJUADp+R+zLicJbRlP8fmrKGr7KYxJ59c1WA18GiUHGPNdWz+zD5PLGUjOVZZbV91oHMIo5Z8CvJyoKgX9DoeTriyeV4j1TrnRNlGU7JOlHwQT+YNlLIU7Z5N1iisfXfITe90W6Ty/ulTy6KW4T3/NSRlHZCqoj4MeN5+mwc8qZT6XazvJK2on7vSGhL9ZLn1vqkJbu8Bw6+w8niEkZJ1qzu1uy03R4dDop7iyfPvWH1gyOmAeMP6eAon4udDejJmxn9Yt6ueRg64vzqwj6NkPXk08S2dWKusskrYqosjew5NlMkmtub1YOb4bpyeu9ga7XQus+ZQ6vZZFntxH8s6i2VNNTXy6qIP+Ocr77CvEfarALXkszmvlDOPO9xRSoNkSMWijlTPAjkCQov8KVnT9qIQqwzT8eyxFHXcbkgp9Q0wOOm7J0JuvjUJEGR3pTXU7nJ4xNN1nJ1NGwUdrb8YePL8YVEv8YgX+hQrciRWWNRPzhh0kLLYS1sWq4OjGZINRYu2qKiJHFapUqiH37xTzenTDkwYz/1GQjqWDVw/vl3sss/J5bSTTqS6/aEtOqQz+5W0WOrvRZa8UGJFhkQLzYx0LOPbXQxilaHbekev8Ly8AksxB9lhZdKjc2RFHSsaoGzaS1lfeXSMaU0ljjfebiPBc4KKYUyYcoPUFuGEKiQnmRBTibcN77BGzViYcqcba7TiZLVkvJWFkT7P1rYVieCzurleIRqaKeqw8Lwdq63/Uay6WMlvQvMcQHZWqFSUYqo+z2jfTyWONxlrvPzQzmm16oL3iDbMDS3bdI5oUu10nXQauqwM9IN0zStNGlqa0KrSdKGfog61qHdvtCZhOkSeAHOS/8JLV4jXE5vJKsV4jTraysdQa6y6rqHZPxlJKSRTDsl0PG4pHydlm84RTaorCrPaDZgiydT3WGXmRzoFvRR1bhsrprKpyYpy2L3RntENRDw9XClEmxb1whXgx7LTZJVivCTu4c9x/ZzlLSaOIqV3DVUKySpPv/OJhOKkbNO5XDvVZ9fRDaYLidb3eO3WjwRTeinq4CKBxv2QU3jw1lsRCFUKToarbhGtMkx9ahl3zf8i5R8yWq+fjFKM1agjPUd9k7Mo2mxLAh+vbNPZsaT67LpkNwxHh+RGidZ3JyMRr91Imilqe+VYQ621WGD3RmsRiEP8tMhiKalUret0W+uxGnUqytaNJPDBvMFOGroTpZBOxZHujiWVxq/TaCSILsmNkqnvuo1EtNg4oJlc26IOhujt3hjXog5l0lBv9n+DlonHh/zmlbgLkVLZHy7d+83FSuKerLJ1c0NUJ0nwnZzrRlL9SUNLWTRtLKtnTGDRtLG+Wv9e1X2nuL3fo1OSqe9+j0TC0deirt1trYTr2DOhS3gxJIm1y3Qsku2l093rx7MEnSxuCOQI7QvyHId3JUMiw1In52b7hJtuUR26WKuJ1ne/RyKR0ExRt7H+N9ZZ/mlIWFF7QbxdpqORbC/thv8xXlysDosbEmnoTs7VRXG0FnTymyda33Xq8EA3Rd3s+qi1ErGAlfdAE0J3GUmUVHppr/2PuixuSKShOzlXJ8WRzaRjN57QuYREF+ckim4jkUjo5aNudn3UWf5p0MaiDvVvOqG4MJA2f6GO/kcvSGRDVCfn6rbBajYS3k6Cu/GA83obPpdQVVPPzn31rm3WmwnoZVGHhuft3oiVvjB6UiK3CLUIckVoVKr5vxMKA7lMPyfybh2h109kqJUJvX66SWRY6uTcTBnmZjLp2I0nnmsxm+YVnKKZog6ZTNxdaeWYjZGA3Q3CJwqDyjmWki5OYGiW7lVS6UKHeNdIJNJBOTm3NXZ44bj5W6djHsDJuW7s7pMsXrQdvRR1aHhegqF5cHCBjelXknDayEQnChPNvxst8uDnTy+nYu0OVzOoRUOXeFeD+8T6rSH10Ua0eYCiwoDjePhYOXyCiP0sftdPr9qOZj5qO+qjodbatSIB/3SkGNl/vL/uoJjZX839pDn+edSMhQf5uhLp+ZPxb0a7fqNSPPH+Ol/iTnWJd3WD0Hj3SL93JpLKM0X7rafP+ywtMeaR5gECOUJ1XYPja0e6RjjKfha/8artaKaobddHY521x1/77o6/6sQSrqlv5IkIyju0wsSLAMgVSWlCL9b1/cpVkq1ha24scElGhnR2FKk+U7TftKqmPi0KJ9LEd/uCvBYx+PGuHX6NaOhQP71qO5q6PmqtPdIKix1/1WnBhCvD8ImJWBn50rGLg5OMf+G4HT6ma9haqr4/vxe4uDEsTvWZnLgVQklG4YTPA5RNeynha7uRw8eNFANetR09LeqaKiuLXpvYO5yEkkrBhFaY0N4cLAsa0hcSF7x+rsPNL70IH9MxbC0d1rDfIwU3hsWpPlO037pT28gZKtOhcGIt03Yy4khH/XQrxYBXbUcvizoY4VG9xfpfUOT4q04s1cibwR9ckdyODIi2dDWcUo+iL3QMW0uHNez3SMGNjiLVZ4q18tStRVXRFmyN6VfiaMSRjvrpVooBr9qOXoo6155M3BtU1M4t6tACi7QiSoATD+/MknW7tFjXH5T3508vjxj658luzmHyhIYH/jRN6VmTJR1Kzu+Mcm50FOl4pliGiBsKJ5oyS0Qxpmo8xdq2Lxg9kmyd8yLkUy9FHYz6qN5q/U/AoobY2ygpYM32Gu44f6A2lqMuSWGiLfn1M0wvHUrO75GCGx2Fm8/kpsKJdO2ferj3YCzffLCO+z0Ci4VjRS0iuUAFUKmUOtsVaUSsCcW9QUVdnNRlYvWMXi54iLTCMdyd4WbDczp5EqpM4k22ukWkGPh0bFzr5wIXt37bbFm046VidLKDuN8jsFgkYlH/BFgBOPdHJENewYGETAlMJoaiQ88YbYWjV5uOOo04cBLW6PbkWyRZn11cyQXHlSa8YEk3skWpuoGXijH4G8TaQdzvEVgsHClqEekFTAB+B/zMVYly82HfNut1gq6PIDr0jLEUoBdWqlP/nxMl7HYHF03W11du9dRPnyjJhA/6sVRf5/QA4J1inDQ0/g7ioR1r6HyN3+Xm1KK+B/gl0ME9UWzy2oBqsl4nMJkYig49YzwFqMsilnhxtV50cMlO4vipgJKJkfZjqb7u6QFijTjc+H2dGnG6lVtcRS0iZwNblFKLRWR0jPOuBK4E6NOnTwoS2ROKufkH4qqTwO8hZzwF6LaVWtw2wM599XHvG6niBicUvQoPTMZVle6GlKhSSCaUy48FOMne0+/kYG5NbDs14vxeLBWOE4t6FHCOiJwFFAAdReQfSqnvh56klHoAeACgvLzcWT7QSARD9AqKiLsRocbEmrwQYEy/EtfuPXdpJXtrGw46HsiVgywHHUYfybiq0tmQklH6yYwC/FiAk8w9/bImvZrYdmLE+b1YKpy4ilopdSNwI4BtUf8iXEmnVyJbUSc5kagL4XHdoSjg2cWVlB/a2ZWKf9f8L6hvOrivbJefF/F+fo8+kuks0tmQklH6yYwCnI5y0kkycvplTeowsR1Eh4CEUPRaQg4HFHWSE4k6MWmotUN1aYQf183sdNEq8y6Hm/D6QbCsnO7mnc7do5NR+okuHU5klJNOklni7Jc1qcPEdhDd0ioktOBFKfUG8IYrkgQJJmZKciJRR7yu+LpZA26QiLsknr81mfJKdBSQ6CgnXSQzWkl2zsCtXNZBgr+vF/5zHVyCoei1MhEOTCDGsah1DTmKhNeKU4fwRLdx2pCc+FuTLa9EXEZujXKctINEXVuJlke6fNpOJrYBz/znfrsEQ9FQUcd3fegWOhMPP3YRB32sAbdw0pCc+Fu9KC83Omu32kEyo4V0+LSd3Hfora9oFY3hFfoq6hiTibqFzsTDD8WpkzXgJ07dTm6XlxudtZvtIB2jhXTksg5l7tLKiJOxyd4rk9BPUTeH5xVHPUW30BknGMXpD7r4693orHVpB16VcazJdy9+Tz/drfop6mBO6hiuD10an0F//PDXR2vQ6e6sdWkHXpVxvCgcN/Hb3apheF5wMjG660O30BmDvkwaevAefunYqScaXu7TqEs78KqMo3VAxYWBiPdK536Vfm8ArZ9FnRvfom4tk2WG9OCl28nL+ROd2oEXZRzNcp9+zrEHnZtuC9hvN5N+ijpoUcdZmdjafb6ZFJ4YjWx4hnC8btCtqR0k0jGlu8P0282koaKOb1G3dvz2l6WDbHiGSPjdoLMdpx1TujtMv9cm6OejDk3KZIhIuvxl6fThJYrfPj+30MVv3NpJZ4oB8H6uIxz9LOru/aHLkdDOvexymU46rAW/LVq/fX5uoZPfuDXj1n6VZgl5kCNOhesq/JZCa9IxvPZ70VA2uwhak99YV7Ktw9RPURvikg5rwW+L1m+fnyH7yaYO0yjqDCQd1oLfFm22WTwGg5uIUslvxhKN8vJyVVFh3Bc6E+6jBsui9XKCxJB5ZGNIpS6IyGKlVHmkz4xF3UoxFq0hUfyegG7NtApFbayAyGSTD8/gPn5PQLdmsl5RGyvAkCqmo7fwewK6NaPfgpc0k60LKwze4GWSJd1J9yKScPxcgKU7WW9RGyvAkApmuH+AdIVURhqhgHdbbGUiWa+o/Q5DM2Q2pqM/QDomoKO5IgsCOaZDjEHWK2qzsMKQCqajb0mqE9DRRijhx4K0xg4xEnF91CJSICIfishyEflMRH7jhWDpwu9kKobMxiRZSi+JKt7W2iGG48Si3g+MVUrtFZEA8I6I/Fsp9b7LsqUNE4ZmSBYTb55eoo1QAAQIXX5nOsQDxFXUylq6uNd+G7D/0r+c0WDQFNPRp4/rxx/N1KeWRfxMYY14TYd4MI581CKSCywGjgD+opT6IMI5VwJXAvTp0yedMhoMhixh0tBSfvOvz9i5r/6gz0qLC1k0bawPUumPozhqpVSjUmoI0As4XkQGRDjnAaVUuVKqvKTE5JI2GAyRuWXiscbvnyAJRX0opapE5A3gDOBTVyQymJVwhqzG+P0TJ66iFpESoN5W0oXAqcCdrkvWSjFL3g2tAeP3TwwnFnUP4DHbT50DPK2UetFdsVovZiWc/5gRjUE3nER9fAwM9UAWA2YlnN+YEY1BR7I+KVOm4XbiG0NsTBIvg44YRa0ZZiWcv5gRjUFHjKLWDLPk3V/MiMagI1mflCkTMTPi/mGSeBl0xChqgyEEE+Nr0BGjqA2+oHMIXLpHNPX19WzYsIHa2tq0XdOQuRQUFNCrVy8CgYDj7xhFbYiKW8o0Xgiczko8GTZs2ECHDh3o27cvIuK3OAYfUUqxfft2NmzYQFlZmePvGUVtiIib8cTxQuB0jmNOphOpra01StoAgIjQpUsXtm7dmtD3TNSHISJuxhPHCoHTOY45lY1ujZI2BEmmLhhFbYiIm/HEsULgdI5j1rkTMWQ3RlEbIuJmPHGsRT06xzF71YnMXVrJqBkLKZv2EqNmLHRksRvi88Ybb/Duu+/6LUZSGEVtiIibKyRjLerReWWmF51IKu6VTKKhocHzexpFbcg63F4hOWloKYumjWX1jAksmja2+bo6r8z0ohNxy72yZs0a+vXrxxVXXMGAAQOYPHkyCxYsYNSoURx55JF8+OGHVFdXc9lllzF8+HCGDh3KCy+80Pzdk046iWHDhjFs2LBmZbdp0yZOPvlkhgwZwoABA3j77bcBaN++ffN9n3nmGaZMmQLAlClT+NnPfsaYMWO44YYbWLVqFWeccQbHHXccJ510EitXrmw+7+qrr2bMmDEcdthhvPnmm1x22WUcc8wxzdcCeOWVVxg5ciTDhg3jwgsvZO9ea8fAvn37cssttzBs2DAGDhzIypUrWbNmDbNmzeLuu+9myJAhvP3228yZM4cBAwYwePBgTj755JTK121M1IchKn6tkNR1ZaYXi2HcdK98/fXXzJkzhwceeIDhw4fz5JNP8s477zBv3jxuv/12+vfvz9ixY3n44Yepqqri+OOP59RTT6Vbt268+uqrFBQU8NVXX3HxxRdTUVHBk08+yfjx47nppptobGxk3759cWX48ssvWbBgAbm5uYwbN45Zs2Zx5JFH8sEHH3DNNdewcOFCAHbu3MnChQuZN28eEydOZNGiRTz44IMMHz6cZcuW0atXL2677TYWLFhAu3btuPPOO/njH//IzTffDEDXrl1ZsmQJ999/PzNnzuTBBx/kqquuon379vziF78AYODAgcyfP5/S0lKqqqpSLl83MYraYEgAtzuRaLt0p8O9UlZWxsCBAwE49thjGTduHCLCwIEDWbNmDRs2bGDevHnMnDkTsMIK161bR8+ePbn22mtZtmwZubm5fPnllwAMHz6cyy67jPr6eiZNmsSQIUPiynDhhReSm5vL3r17effdd7nwwgubP9u/f3/z64kTJzbL1r179xZyB2X9/PPPGTVqFAB1dXWMHDmy+fvnn38+AMcddxzPPfdcRFlGjRrFlClT+N73vtd8vq4YRW0waISbuUbatGnT/DonJ6f5fU5ODg0NDeTm5vLss89y9NEt7zV9+nS6d+/O8uXLaWpqoqCgAICTTz6Zt956i5deeolLLrmE66+/nh/84Actws/CV2O2a9cOgKamJoqLi1m2bFlMWUPlDJf1tNNOY/bs2TG/n5ubG9UfPmvWLD744ANeeuklhgwZwrJly+jSpUvEc/3G+KgNBo3w00c/fvx4/vSnP6GUAmDp0qUA7Nq1ix49epCTk8Pf//53GhutTmTt2rV069aNH/7wh1x++eUsWbIEgO7du7NixQqampp4/vnnI96rY8eOlJWVMWfOHMBasbd8+XLHso4YMYJFixbx9ddfA7Bv375mSz8aHTp0YM+ePc3vV61axQknnMCtt95K165dWb9+veP7e42xqA0GzfDLR//rX/+aqVOnMmjQIJRS9O3blxdffJFrrrmGCy64gDlz5jBmzJhmq/iNN97grrvuIhAI0L59ex5//HEAZsyYwdlnn03v3r0ZMGBA8yRfOE888QRXX301t912G/X19Vx00UUMHjzYkawlJSU8+uijXHzxxc0uk9tuu42jjjoq6ncmTpzId7/7XV544QX+9Kc/cffdd/PVV1+hlGLcuHGO7+0HEuw900l5ebmqqKhI+3UNhkxkxYoVHHPMMX6LYdCISHVCRBYrpcojnW9cHwaDwaA5RlEbDAaD5hhFbTAYDJoTV1GLSG8ReV1EVojIZyLyEy8EMxgMBoOFk6iPBuDnSqklItIBWCwiryqlPndZNoPBYDDgwKJWSm1SSi2xX+8BVgD6re81GAyGLCUhH7WI9AWGAh9E+OxKEakQkYpEdy8wGAz68Oijj7Jx48bm91dccQWff576AHrNmjU8+eSTze8rKir48Y9/nPJ1YzFnzhyOOeYYxowZ4+p93MaxohaR9sCzwFSl1O7wz5VSDyilypVS5SUlJemU0WAweEi4on7wwQfp379/ytcNV9Tl5eXcd999KV83Fg899BD3338/r7/+uqv3cRtHKxNFJIClpJ9QSkXOcGIwGOLz72nw7SfpveYhA+HMGTFP+cc//sF9991HXV0dJ5xwAvfffz8Al19+ORUVFYgIl112Gb1796aiooLJkydTWFjIe++9x5lnnsnMmTMpLy+nffv2/OhHP2LBggV06tSJ22+/nV/+8pesW7eOe+65h3POOYc1a9ZwySWXUF1dDcCf//xnTjzxRKZNm8aKFSsYMmQIl156KUOHDmXmzJm8+OKL7Nixg8suu4xvvvmGtm3b8sADDzBo0CCmT5/OunXr+Oabb1i3bh1Tp06NaIXPnj2b22+/HaUUEyZM4M477+TWW2/lnXfeYfXq1Zxzzjncddddzee/8cYb3HLLLXTv3p1ly5Zx/vnnM3DgQO69915qamqYO3cuhx9+OFu3buWqq65i3bp1ANxzzz2MGjWKDz/8kKlTp1JTU0NhYSGPPPIIRx99NI8++ijz5s1j3759rFq1ivPOO4/f//73qf/GSqmYf4AAjwP3xDs3+Hfccccpg8Fg8fnnnx948/INSj18Vnr/Xr4h7v3PPvtsVVdXp5RS6uqrr1aPPfaYqqioUKeeemrzeTt37lRKKXXKKaeojz76qPl46HtAvfzyy0oppSZNmqROO+00VVdXp5YtW6YGDx6slFKqurpa1dTUKKWU+vLLL1VQH7z++utqwoQJzdcNfX/ttdeq6dOnK6WUeu2115qvdcstt6iRI0eq2tpatXXrVtW5c+fm5whSWVmpevfurbZs2aLq6+vVmDFj1PPPPx/xWULvXVRUpDZu3Khqa2tVz5491c0336yUUuqee+5RP/nJT5RSSl188cXq7bffVkoptXbtWtWvXz+llFK7du1S9fX1SimlXn31VXX++ecrpZR65JFHVFlZmaqqqlI1NTWqT58+at26dRF/k3CAChVFpzqxqEcBlwCfiMgy+9j/KqVeTr2baEkyOzwbDBlFHMvXDV577TUWL17M8OHDAaipqaFbt25MnDiRb775huuuu44JEyZw+umnx71Wfn4+Z5xxBmDlc27Tpg2BQKA5VSpAfX19xLSosXjnnXd49tlnARg7dizbt29n165dAEyYMIE2bdrQpk0bunXrxubNm+nVq1fzdz/66CNGjx5N0OU6efJk3nrrLSZNmhTznsOHD6dHjx4AHH744c3PP3DgwGZXyYIFC1r453fv3s2ePXvYtWsXl156KV999RUiQn19ffM548aNo6ioCID+/fuzdu1aevfuHbcMYhFXUSul3sGyql0luAVRML1jcAsiwCjrLMd00O6ilOLSSy/ljjvuOOiz5cuXM3/+fP7yl7/w9NNP8/DDD8e8ViAQaE5jGilVKsDdd98dMS1qPBnDCd4nNM1ppLSlkb7rhHhpX8FKx/ree+9RWNgyH/h1113HmDFjeP7551mzZg2jR4+OeN1YaVYTQZuViWaH59ZJa9kj0E/GjRvHM888w5YtWwDYsWMHa9euZdu2bTQ1NXHBBRfw29/+tjlNaXg60ESJlhY11nVPPvlknnjiCcDyH3ft2pWOHTs6ut8JJ5zAm2++ybZt22hsbGT27NmccsopScsfyumnn86f//zn5vfB/Nm7du2itNQyJh599NG03CsW2ihqr3Z4NuiF6aDdp3///tx2222cfvrpDBo0iNNOO41NmzZRWVnJ6NGjGTJkCFOmTGm2uKdMmcJVV13FkCFDqKlJvP1dc801PPbYY4wYMYIvv/yyOS3qoEGDyMvLY/Dgwdx9990tvjN9+nQqKioYNGgQ06ZN47HHHnN8vx49enDHHXcwZswYBg8ezLBhwzj33HMTljsS9913X7Nc/fv3Z9asWQD88pe/5MYbb2TUqFHNHZGbaJPmdNSMhRG3ICotLmTRtLHpEs2gGWXTXiJSDRRg9YwJXovjCibNqSGcjE1z6sUOzwb9iLYXYDr2CDQYsgVtFLWfWxAZ/MN00AZDfLTaisuvLYgM/hH8vbM96kMp1WLTV0PrJRl3s1aK2tA6yfYOuqCggO3bt9OlSxejrFs5Sim2b9/uKGQxFKOoDQaX6dWrFxs2bMAkKzOA1XGHLthxglHUBoPLBAIBysrK/BbDkMFoM5loMBgMhsgYRW0wGAyaYxS1wWAwaI4rKxNFZCuwNsGvdQW2pV2Y9KCrbEauxDByJY6usmWjXIcqpSLuuuKKok4GEamItnzSb3SVzciVGEauxNFVttYml3F9GAwGg+YYRW0wGAyao5OifsBvAWKgq2xGrsQwciWOrrK1Krm08VEbDAaDITI6WdQGg8FgiIBR1AaDwaA5WihqETlDRL4Qka9FZJqPcvQWkddFZIWIfCYiP7GPTxeRShFZZv+d5YNsa0TkE/v+FfaxziLyqoh8Zf/v5LFMR4eUyTIR2S0iU/0qLxF5WES2iMinIceilpGI3GjXuS9EZLzHct0lIitF5GMReV5Eiu3jfUWkJqTsZnksV9TfzufyeipEpjUissw+7mV5RdMP7tcxpZSvf0AusAo4DMgHlgP9fZKlBzDMft0B+BLoD0wHfuFzOa0BuoYd+z0wzX49DbjT59/xW+BQv8oLOBkYBnwar4zs33U50AYos+tgrodynQ7k2a/vDJGrb+h5PpRXxN/O7/IK+/wPwM0+lFc0/eB6HdPBoj4e+Fop9Y1Sqg74J5CenSkTRCm1SSm1xH69B1gB6Jwo+VwguAvoY8Ak/0RhHLBKKZXoitS0oZR6C9gRdjhaGZ0L/FMptV8ptRr4GqsueiKXUuoVpVSD/fZ9ILG8ly7JFQNfyyuIWAm9vwfMduPesYihH1yvYzoo6lJgfcj7DWigHEWkLzAU+MA+dK09TH3YaxeDjQJeEZHFInKlfay7UmoTWJUI6OaDXEEuomXj8bu8gkQrI53q3WXAv0Pel4nIUhF5U0RO8kGeSL+dLuV1ErBZKfVVyDHPyytMP7hex3RQ1JG2vPA1ZlBE2gPPAlOVUruBvwKHA0OATVhDL68ZpZQaBpwJ/EhETvZBhoiISD5wDjDHPqRDecVDi3onIjcBDcAT9qFNQB+l1FDgZ8CTItLRQ5Gi/XZalBdwMS0NAs/LK4J+iHpqhGNJlZkOinoD0DvkfS9go0+yICIBrB/hCaXUcwBKqc1KqUalVBPwN1wa8sVCKbXR/r8FeN6WYbOI9LDl7gFs8VoumzOBJUqpzbaMvpdXCNHKyPd6JyKXAmcDk5Xt1LSHydvt14ux/JpHeSVTjN9Oh/LKA84Hngoe87q8IukHPKhjOijqj4AjRaTMtswuAub5IYjt/3oIWKGU+mPI8R4hp50HfBr+XZflaiciHYKvsSaiPsUqp0vt0y4FXvBSrhBaWDl+l1cY0cpoHnCRiLQRkTLgSOBDr4QSkTOAG4BzlFL7Qo6XiEiu/fowW65vPJQr2m/na3nZnAqsVEptCB7wsryi6Qe8qGNezJY6mE09C2sGdRVwk49yfAdraPIxsMz+Owv4O/CJfXwe0MNjuQ7Dmj1eDnwWLCOgC/Aa8JX9v7MPZdYW2A4UhRzzpbywOotNQD2WNXN5rDICbrLr3BfAmR7L9TWW/zJYz2bZ515g/8bLgSXARI/livrb+Vle9vFHgavCzvWyvKLpB9frmFlCbjAYDJqjg+vDYDAYDDEwitpgMBg0xyhqg8Fg0ByjqA0Gg0FzjKI2GAwGzTGK2mAwGDTHKGqDwWDQnP8PPQ9J63Lr5N0AAAAASUVORK5CYII=",
"text/plain": [
"