{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "# 无迹卡尔曼滤波器Unscented Kalman Filter\n", " 本节主要介绍UKF这一滤波器,在书中介绍到,这一方法相比于EKF这一传统的非线性卡尔曼滤波器有诸多有点,比如说计算量的降低。因此,我本着做工程的需要,便没有过多研读EKF这一滤波器,有兴趣深入研究的同学,可以回到原著中一窥究竟,我就不讲了,哈哈哈哈哈~\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Compared UKF to KF \n", "我们先给出这一方法的算法过程:\n", " \n", " 先看UKF-KF对照表!!!\n", "\n", "$$\\begin{array}{l|l}\n", "\\textrm{Kalman Filter} & \\textrm{Unscented Kalman Filter} & \\\\\n", "\\hline \n", "& \\boldsymbol{\\mathcal Y} = f(\\boldsymbol\\chi) & State Transition Function \\\\\n", "\\mathbf{\\bar x} = \\mathbf{Fx} & \\mathbf{\\bar x} = \\sum w^m\\boldsymbol{\\mathcal Y} & StateEstimation \\\\\n", "\\mathbf{\\bar P} = \\mathbf{FPF}^\\mathsf T+\\mathbf Q & \\mathbf{\\bar P} = \\sum w^c({\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})(\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})^\\mathsf T}+\\mathbf Q \\\\\n", "\\hline \n", "& \\boldsymbol{\\mathcal Z} = h(\\boldsymbol{\\mathcal{Y}}) & Measurement Function \\\\\n", "& \\boldsymbol\\mu_z = \\sum w^m\\boldsymbol{\\mathcal{Z}} & MeasurementEstimation \\\\\n", "\\mathbf y = \\mathbf z - \\mathbf{H\\bar x} &\n", "\\mathbf y = \\mathbf z - \\boldsymbol\\mu_z \\\\\n", "\\mathbf S = \\mathbf{H\\bar PH}^\\mathsf{T} + \\mathbf R & \n", "\\mathbf P_z = \\sum w^c{(\\boldsymbol{\\mathcal Z}-\\boldsymbol\\mu_z)(\\boldsymbol{\\mathcal{Z}}-\\boldsymbol\\mu_z)^\\mathsf{T}} + \\mathbf R \\\\ \n", "\\mathbf K = \\mathbf{\\bar PH}^\\mathsf T \\mathbf S^{-1} &\n", "\\mathbf K = \\left[\\sum w^c(\\boldsymbol{\\mathcal Y}-\\bar{\\mathbf x})(\\boldsymbol{\\mathcal{Z}}-\\boldsymbol\\mu_z)^\\mathsf{T}\\right] \\mathbf P_z^{-1} \\\\\n", "\\mathbf x = \\mathbf{\\bar x} + \\mathbf{Ky} & \\mathbf x = \\mathbf{\\bar x} + \\mathbf{Ky}\\\\\n", "\\mathbf P = (\\mathbf{I}-\\mathbf{KH})\\mathbf{\\bar P} & \\mathbf P = \\bar{\\mathbf P} - \\mathbf{KP_z}\\mathbf{K}^\\mathsf{T}\n", "\\end{array}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Predict Step\n", "\n", "The UKF's predict step computes the prior using the process model $f()$. $f()$ is assumed to be nonlinear, so we generate sigma points $\\mathcal{X}$ and their corresponding weights $W^m, W^c$\n", "according to some function:\n", "\n", "$$\\begin{aligned}\n", "\\boldsymbol\\chi &= \\text{sigma-function}(\\mathbf x, \\mathbf P) \\\\\n", "W^m, W^c &= \\text{weight-function}(\\mathtt{n, parameters})\\end{aligned}$$\n", "\n", "We pass each sigma point through $f(\\mathbf x, \\Delta t)$. This projects the sigma points forward in time according to the process model, forming the new prior, which is a set of sigma points we name $\\boldsymbol{\\mathcal Y}$:\n", "\n", "$$\\boldsymbol{\\mathcal{Y}} = f(\\boldsymbol{\\chi}, \\Delta t)$$\n", "\n", "We compute the mean and covariance of the prior using the *unscented transform* on the transformed sigma points. \n", "\n", "$$\\mathbf{\\bar x}, \\mathbf{\\bar P} = \n", "UT(\\mathcal{Y}, w_m, w_c, \\mathbf Q)$$\n", "\n", "These are the equations for the unscented transform:\n", "\n", "$$\\begin{aligned}\n", "\\mathbf{\\bar x} &= \\sum_{i=0}^{2n} w^m_i\\boldsymbol{\\mathcal Y}_i \\\\\n", "\\mathbf{\\bar P} &= \\sum_{i=0}^{2n} w^c_i({\\boldsymbol{\\mathcal Y}_i - \\mathbf{\\bar x})(\\boldsymbol{\\mathcal Y}_i-\\mathbf{\\bar x})^\\mathsf{T}} + \\mathbf Q\n", "\\end{aligned}\n", "$$\n", "\n", "This table compares the linear Kalman filter with the Unscented Kalman Filter equations. I've dropped the subscript $i$ for readability.\n", "\n", "$$\\begin{array}{l|l}\n", "\\text{Kalman} & \\text{Unscented} \\\\\n", "\\hline \n", "& \\boldsymbol{\\mathcal Y} = f(\\boldsymbol\\chi) \\\\\n", "\\mathbf{\\bar x} = \\mathbf{Fx} & \n", "\\mathbf{\\bar x} = \\sum w^m\\boldsymbol{\\mathcal Y} \\\\\n", "\\mathbf{\\bar P} = \\mathbf{FPF}^\\mathsf T + \\mathbf Q & \n", "\\mathbf{\\bar P} = \\sum w^c({\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})(\\boldsymbol{\\mathcal Y} - \\mathbf{\\bar x})^\\mathsf T}+\\mathbf Q\n", "\\end{array}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Update Step\n", "\n", "Kalman filters perform the update in measurement space. Thus we must convert the sigma points of the prior into measurements using a measurement function $h(x)$ that you define.\n", "\n", "$$\\boldsymbol{\\mathcal{Z}} = h(\\boldsymbol{\\mathcal{Y}})$$\n", "\n", "We compute the mean and covariance of these points using the unscented transform. The $z$ subscript denotes that these are the mean and covariance of the measurement sigma points.\n", "\n", "$$\\begin{aligned}\n", "\\boldsymbol\\mu_z, \\mathbf P_z &= \n", "UT(\\boldsymbol{\\mathcal Z}, w_m, w_c, \\mathbf R) \\\\\n", "\\boldsymbol\\mu_z &= \\sum_{i=0}^{2n} w^m_i\\boldsymbol{\\mathcal Z}_i \\\\\n", "\\mathbf P_z &= \\sum_{i=0}^{2n} w^c_i{(\\boldsymbol{\\mathcal Z}_i-\\boldsymbol{\\mu}_z)(\\boldsymbol{\\mathcal Z}_i-\\boldsymbol{\\mu}_z)^\\mathsf T} + \\mathbf R\n", "\\end{aligned}\n", "$$\n", "\n", "Next we compute the residual and Kalman gain. The residual of the measurement $\\mathbf z$ is trivial to compute:\n", "\n", "$$\\mathbf y = \\mathbf z - \\boldsymbol\\mu_z$$\n", "\n", "To compute the Kalman gain we first compute the [cross covariance](https://en.wikipedia.org/wiki/Cross-covariance) of the state and the measurements, which is defined as: \n", "\n", "$$\\mathbf P_{xz} =\\sum_{i=0}^{2n} w^c_i(\\boldsymbol{\\mathcal Y}_i-\\mathbf{\\bar x})(\\boldsymbol{\\mathcal Z}_i-\\boldsymbol\\mu_z)^\\mathsf T$$\n", "\n", "And then the Kalman gain is defined as\n", "\n", "$$\\mathbf{K} = \\mathbf P_{xz} \\mathbf P_z^{-1}$$\n", "\n", "If you think of the inverse as a *kind of* matrix reciprocal, you can see that the Kalman gain is a simple ratio which computes:\n", "\n", "$$\\mathbf{K} \\approx \\frac{\\mathbf P_{xz}}{\\mathbf P_z} \n", "\\approx \\frac{\\text{belief in state}}{\\text{belief in measurement}}$$\n", "\n", "Finally, we compute the new state estimate using the residual and Kalman gain:\n", "\n", "$$\\mathbf x = \\bar{\\mathbf x} + \\mathbf{Ky}$$\n", "\n", "and the new covariance is computed as:\n", "\n", "$$ \\mathbf P = \\mathbf{\\bar P} - \\mathbf{KP_z}\\mathbf{K}^\\mathsf{T}$$\n", "\n", "This step contains a few equations you have to take on faith, but you should be able to see how they relate to the linear Kalman filter equations. The linear algebra is slightly different from the linear Kalman filter, but the algorithm is the same Bayesian algorithm we have been implementing throughout the book. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "怎么样,我们发现,其实UKF方法和KalmanFilter有着很多相似之处,可能唯一有出入的地方大概是每个状态的建立过程,其中包含了:$\\bar x $ 、 $\\bar P$ 、 $\\mu_z$ 和 K。这些不同也几乎都最后指向了$w^m、w^c$是个什么鬼了。 \n", "这里我就要对这两个权系数引入点介绍:这些权系数其实是作用在描述当下状态的几个sigma点上,通过对sigma点的状态加权而得到当下状态。 \n", "那么,问题来了\n", "\n", "**What is sigma point?** \n", "\n", "我们看个图:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWQAAADuCAYAAAAOR30qAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3WmQJFd97/1v1r4v3V29r7OvkkYCSdhIw2aZTSErANkY\nMBhCCptwgM315TEYYTCYCzfsiy6+F/uFA+MwRkjGCEEI0ANIDxgQEpKFNKPZemZ6767q7tq3rMrK\nzOdFTY9mRj0zPaNasrv/n4gJIkRP1ek5mb88+c+T5yimaSKEEKL9bO1ugBBCiDoJZCGEsAgJZCGE\nsAgJZCGEsAgJZCGEsAgJZCGEsAgJZCGEsAgJZCGEsAgJZCGEsAjHlfxwV1eXOTo62qSmCCHExvTM\nM88sm6YZu9zPXVEgj46O8vTTT199q4QQYhNSFGVqLT8nJQshhLAICWQhhLAICWQhhLAICWQhhLAI\nCWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQh\nhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAI\nCWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQh\nhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAI\nCWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAICWQh\nhLAICWQhhLAICWQhhLAICWQhhLAICWQhhLAIR7sbIITYGEzT5Oc//zn3338/iUSCzs5O3vGOd/C6\n170Om03GfmshgSyEeNmeffZZ3v72t5NIJCiVSpimCcDXv/51wuEw999/P7fcckubW2l9EsgNZBgm\nhXKVbFGlpGpUazoVTaeq6Zimid1uw25TsNts+DxOAl4XQa+LgNeF3S4jiGZQqzWyBZV8uUqlWqNa\nq/dHTTewnekLm03B5bDX+8PnIuhz43HJqbFWzz77LLfeeiuFQuEl/1+hUKBQKPDGN76RRx55hFtv\nPUiuVCFbUClXa1Q1/WyfAGf6pN4v/pVzxOcm4HVhsymt/tVaTo66l8EwTJK5Eol0kVSuTK5UoVQy\nKRZBrUBNg1qt/scwwWYDu63+v24PeD3g9YLXqxCLeOmJBuiJ+vF7Xe3+1datkqoRTxVYyhTJFisU\nyjVKRSiVQKvV+0TTQDfAptT7wmYDhwO8vhf7JBxw0tNR74/OkG9ThMHVME2TO+64Y9UwrlPA6aOk\nBLj93X/Cv/zbA2iag2IJKuqZPjnTL/Bif9js4HGvnB/g8yl0R/z0RP30dAQ27AVzY/5WTWQYJvFU\ngYVknsVMkUzWIJOBXK5+0jvtLvxuNx6nC4/djtPpwOmxA/WD1zBNdMNALVVJZ6vMV6tU9SqhUIlI\npEQkAj0dXrYOdNAT9aMoEgSXky9VmFvOE08VSGUrZDKQyUKxAHrNjs/txu9y43Q48NrtOLx2HDY7\nJvW+MEwTrVajlKqSrFYoV6u43BrRaJpIJE00YmesN8JobwT3Bg2Cq/X444+TTqfP/4+KDVwBcAXB\n6QMlBESoGjG+/vUTXH/tAXxnzhGfzY7DaT97jhimiWEY6KZBpaSRzFQoVavUzCrhUIFotEA4nKC/\ny8/WgQ66wr7W/9JNJEfXGpVUjcl4hpnFLMspnWQSMhlwKG6i/gCDIR/+mAf7VTy8qOk62VKR9FKB\n6eki4XCZ6fgcvV0utg92MhgLNeE3Wt8Mw2Q+mWcyniGRLLOchEwayiUbYZ+fqC/AcK8Xt9N5xZ9t\nmibFikqmVGRiucApW4W53iQnulOM9ITZOdQpwXzGV7/61RdHx3YXeKLgDoLSAXQAEagBVY1qtszJ\np57lvW+984q/p1qrkS0VScULTE4UmeksMhUv0tflYcdgJz0dgUb+Wm0jR9Vl5EsVjs8kmUnkWVqC\nxSWwG266gmH6+wJXdcJfyGG30xkM0RkMoRsGS7ksJ4+nmZmuEl9eYLg3w7XbeglIKQNdN5iMZzg1\nn2ZxucbiIuRyNjr8IQZDQQLdXmwv865CURQCHi8Bj5fBji7y5TILiynm5grM92aYWcyxdyzGSE94\n09/BLCwsgMMLvk5whoFuIAaaDSpV0Ir1et0Z2Uzmqr7H5XAQC4WJhcJous5iLsPRF9LMhlQWlucY\n6w+wf0vPui9lrO/WN1GxXOX4TJLJhRwLC7C0pBDxBdnaGSXg8TTte+02G72RKN3hCMlCjvHjyywv\nlUnmJtm/JcZYX7Rp321lhmEyvZhlfDbJfLzG3BwoupuecJQtw8GrujNZq6DXS9A7QKlaYTa5zK+X\nCuTyCRYG8ly/vW/TjpazBRXd3wfhrcAAmB1QqYGq1ov0q/B6vS/7e512OwPRTnrDUZZyWY6+sEwy\nWWApU+a6bT30dwVf9ne0y+Y8ki5Bq+kcnVrm9HyWhQWTREKhKxBh/2AHLkfr/rlsikIsGCbqCzCd\nXOL557OUy4vkS1X2b+neVCOzhWSeI5NLzC9qzM4CNQ9DHV1EfP6WtsPncrOjb4BUIc/4iQTZbImi\nOs1NuwcI+twtbUs7qdUaL0wsMrmQZ3THm3E9aaeazYNaANO86N9zuz0cPHhrw9qxMnjpCASZXErw\nXKZASZ3n2m1dbB/sbNj3tJIE8jniqQKHTieYmq0xN6sQ9YXZN9DZkLLE1XLY7Wzp7iWZ93HsWBxN\ny1CuaNywsx/HBp8qV6nWODSxyOm5PJMTUKu4GezoIupvb72wIxAk6PEynpjncKVMVZvhxt39G+4B\n02qmE1kOTywxPaMTjytct+MgDxQeplquXPbvmqbJa1/72oa3yeVwsKNvgHgmzZEXFqlpy5RUjf1b\netbd7BgJZOon/uEzJ/7pCUDzsquvB5/LOqOezmAIl8PJyVNzVKtFdGOOm/cMrrsDbq1mFs858Rds\nDERjdMesU7N1Ohzs7Bvk9GKcwy/kqemzvHr/EB2hl39LbkUlVeO5U3Em50tMToCLAPsGunE7nfzp\nhz7M3/7t31KtXjyUXS4399xzD15v8y5avZEobqeT48cX0GpZdMPkwPZeyxwza7HpA3l2Kceh04tM\nz+oszNsYiHbRE4tYshODXi+7+4c5NjOD3V7C44qvuwPucsoVjedOJZicKzIxAU787B3oaetdysXY\nbTa29fQxuWzn+IkMDvsct147vKEevpqmycRChhcml5meMUgu2Rnu6KYz+OLMn5tvvpk/+7M/44tf\n/CKKolCpqGf/P5fLjWma3H333dx2221Nb2/UH8BpH+TEyVnsthxet4PdI7Gmf2+jbNpA1nWDZ0/G\nOTVTHxU7DD/7LHrin8vjdLG9d4DjEzM4nTl8Hie7hrva3ayGWEjm+a8TcaamDZYW7Qx3dtMVtPaU\nP0VRGO3q5mSixomTBVzOWW7ZP7whHvRVNZ1fHZtjcr7MxAQEHCH2DcRwrvIs5Td/8ze5/voDPP74\n4/zoRz8mn8/j9/s5ePAgb3jDGwgGW/egLeDxsrW7n5Mn53A6U/jcTkZ6Iy37/pdj/R81V6FSrfHU\nsTlOTqrMTNsZ7IgRC4bb3aw187s9Zw84tztJV9i37uuXp+fTPHtikRPj4CbI/sHuVU98K1IUhS3d\nfRybn+HUhErIl+DG3QPtbtbLUixXefLoHOOnqyTiDka7ei5bu/d6fbz5zW/hzW9+S4taeXFhn5/B\nSA/HT8RxuxfpCHnXxYPXjf1UaBWFcpWfHZrmheMqczNOdvUNr6swXhH2+ekJdjIxAb8+Gad2kWlG\nVmeaJocnFnn66CJHjkCHu4vtvf3rJoxX2G02tvcOkIjbmZgrMLuUa3eTrlo6X+Y/D01z+FiV1KKb\nvQMjbX+QejVioTAhV5jJKZNfn4yfXfDIyjZVIKdyZX763DSHj2pklj3s6R/G61q/9b7eSAea6mF6\nVuPo1FK7m3PFdN3g6ePzPHc8zbFjCkORPvqj63O6EtSf9g92xDg9AYdOL6JWa+1u0hVbSOb5z+dn\nOPyCTrXgZ/fAcEunezbaUEeMTMrB9LzKqfn05f9Cm22aQJ5fzvOzQzO8cFSnVgqwq39o3Y3CLmRT\nFLbEepmdVRifzZAvXX7qkVVUqjV+8cIMh08UOHXSztbY4HkPitarWDCM0/QzPatzYibZ7uZckdPz\naX5xaJ4Xjpi4zDA7egea+sJNKzjsdrbEepmYgGPTSSoWv0iu73/tNTo5l+KJw/O88IKJx4ywvbd/\n3R9oK3xuN1FvmIUF1k0ArFY2CjVxOlSrDXd2E19QmFjIUq5o7W7OZa1WNhqLbZzZO2GfH689QDxh\nWH6UvDFS6SJM0+TQ6QTPHF3ihSPQ6Y0xGuvZMAfaiv5oJ0uLClOJvOVHyel8mf98fuOUjVbjdbkI\ne4LE4ybjs6l2N+eSDMPcUGWji+mPdjI/D6cXMpYeJW/oQD4+k+TwyQzHjimMRPvpi3S0u0lN4XI4\n6PSHiS/U7wasqlCu8sQLc7xwZOOUjS6mP9pJIlEfJVu5lvzcqThHT22sstFqAh4Pfkd9lDwRv7oF\njlphwwbyzGKWQyeTjI8rbIkN0BFYvwuOrEVfpIOlJZhbLlhyxkVV03nq6BwnT+nYaoENVTZajdfl\nIuDyk0yZzC/n292cVZ2YSXJsIsfEhI0dvYMbqmy0mt5IlKXF+stgVp1xsSHPiOVsiWeOJzh+AgYi\n3S1fhKYd3E4nHoeXdNognrrY7g3tYRgmvzpWn9NazHnY0t234cpGq+kMhEguY8kpcLNLOZ4/ucz4\nuMJYVx9+d/NWMLSKoMdLreokmdZI5crtbs6qNlwg50sVnjo6z4kTJlFPB92h9fGGTiN0BUMsWywA\nTNPk2fEFxifLLMadG+LJ/VpF/AFKRTtLadVStf1ktsQzx+McPw794e51Ocf4aiiKUj9HktY6R861\noc6MSrXGU0fnOHFSx2EEGerYGK8Ur1WHP0gmo5BIl9BqerubA9Tr+Mcn80xN2tjRO7Cu57ReKZui\nEPUHSKWwzF1LoVzlyTMDlog7Sk948wxYoL5SXypZ7w8rli02TCDrusFTx+Y4cVpDzXvY0r1xpu2s\nlcNux+fyUMibpPPq5f9Ck00n6nX8kycVtnb343Nb/9XVRgt5/eRykLTALXJV03nyyCzjJ3XseoCh\nzvWz6E6j+FxuTMNBrqBTKFfb3ZyX2BCBbJomz56Mc3JSJbm4uW6LLxT0eMnnIZkrtbUdS5ki/3Xi\nxTp+eBPU8VcT8njJF+pvibZzRGYYZn39lgmNUn7z1PFXE/L4zvaJ1WyI1Dpx5rZ4ZtrO9t6BDTuV\nai2CXh/5fHsPtpKq8atjC5uyjn8hp8OBQ3GRLxhki+2rIz9/OsHJTVjHX03A4yVvkbuWC637XskU\nVI5Mpjh1ErbE+iy1qHw7BNweCoX6v0s7RmSmWV/IZWJKx2EENl0dfzVBj/dsn7RDPFVgfDrL1NTm\nq+Ovpt39cSnrOpB13eDZ8QUmJ006/R2b9rb4XA67HbviQK2YbXkhYTKeYXKuRHLJwegGev325fA4\nXaiV+pKWrVbVdJ47meD0BPRHujZlHf9CbqeTShWKqma5B3vr+lJ5fCbJ1FyVUt7F2Drd1LAZ3E4n\nlUqNkqrhdbduwf1iucrhiWVOT8BoVw9Ou71l321lbqeTYgVKbVjX4vnTCSZnaig1Hz2xzVs6WqHr\nOocPH+aFyVMUC1n29HvZOjrU7madtW4DOZUrc2wqxfSUwvbuvk1dE7uQ2+FErZQpVTRadZlaKVVM\nTRkEnKFNM7d1LTxOJ5V8vbbeSvPLeU7P5onP29g7sLnvVnRd55vf/CYPP/wwuq6j+z3Y3Kf5P5/7\nCG84+Bv83d/9Hbt27Wp3M9dnIBuGWb/yT0Es0EHAs/HfMroSbqeTaqW1ATCdyDI1XyaVdLB/sLsl\n33nq1CkeeughDh06hGEYjI6Ocuedd3Lddddhs9AF2u2o3yK3coSs1XQOnV5kYgIGO2KW35qsmWq1\nGp/61Kc4duzYixux2k0wDaiZfP/73+enP/0pP/7xj7nxxhvb2tZ1GcgTC2lmFiqU8i62Dkmp4kIO\nm52qDlqL1rSoVGscmVpmaqq+9KSjyaUK0zT5yle+wve//300TcM067/nc8/9muPHj7Nz507uvfde\nXBZZRc5us6HXoKYbmKbZkpHq8Zkks/M1qHnX5Y44jfTVr371/DAGMEzADooN0zQpFAr89m//NtPT\n0y3d/+9C1hlGrFG5onF0OsnUFIx0dWPbxLdhF2O32dB1WrbI0JGpJWZmdRymn84WLOL0ne98h+9/\n//tUq5WzYbxCVcscOXKE++77YtPbsVaKorS0T7IFlfGZDLOzCqNdG2+52Suhqio/+MEPzg9jAPPF\nQF6haRpf+9rXWtvAC6y7QD48scjMrIHPHtwUiwZdDZuiYBj10k6zJbMlTs/lWFhQGI01v1RRq9X4\nxje+8dIT7ByaVuXJJ59icXGx6e1ZK6VFfWKa9XLe9LRJVyC66WdVPPHEE5coXylwzsWqWCxy3333\ntaZhF7GuShbpfJmphQKLcRv7WlSnXKHrOr/61a/47ne/y+LiIl6vl9e+9rX81m/9FoGAxR5gKQqG\nCUYLpvQcmVpiehq6gx14nM0vETz33HMYxuXX6TBNk8cee4zf+73fa3qb1kJRFMwW9Mn8cp7ZuEo2\n4+AaKeeRSCRQ1VXmG5tQH4+ef/cwNzfXimZd1LoK5JNzKRYWoDsUbenk9kwmzcc+9nGSySSq+uLb\nPf/2b/Pcf//9fOxjH+PAgQMta89lmSY2haaXc5LZEvFllXzOzthQaxb/TyaT6Gu47a/VNBKJRAta\ntDb12nHz++TUfJr5eRiIdsrMI8Bms2Gz2V56EVegnsrnXyDtbZ6quW56rFCuMrNYIJlU6Gnhq7i6\nrvOxj32chYWF88IYoFqtoKpl/uZv/oaJiYmWtelydNPAZgeHvbnde+4FslUnv9/vx76G30tRbIRC\n1tmUQDea3ydLmSLxZZVi3kHXJn+Qt2Lnzp2rP9xVFECvz7Q4xw033NCahl3EugnkU3MpEgno8Idb\nulbFk08+STKZvORtcrVa5Rvf+EbL2nQ5umFgb/LJnytWmFsqkkrZWrpWxXXXXUdtDUuLOp1OXv3q\nW1rQosurlylMHHZlTReTq3VyLsX8AvSEI/Kw+4z9+/fj862yE8oqgRwIBPjoRz/ausatYl0Eslqt\nMZXIkUhAXyTa0u/+7ne/+5KR8UuZ/OpXv6Jcbu8Kayt0w8Bua24gn5pPEY9DVyDc0jfy/H4/r3vd\n63BdYs0Sm81Of38/27dvb1m7LsVoweg4W1CZWyqRTbf2Aml1NpuND3zgAy89XhQ4N5BdLhc7d+7k\ntttua3kbz7Uuasin59MkEiYhd7AlD47OtdY6pMPhIJ3O4LXAvmQ1XcfraF4AlCsa04k8S0sKe/tb\ne4EEuPvuu5mcnOT06dNo2vnrQ9jtDsLhMJ/85CdX/bvFYpGf/OT/49Sp0zidDq6//gZuuOGGptYO\na4aOo8mBfHIuRSIOXcFI0+eBrze33HILmUyGf/mXf6FWq9XvdhUbUAPTwO/3s337dn74wx+2/YUi\nyweyVtOZjGeIx2FbrLW7Rl/JwiOapp2d75lMJllaWsLj8TA8PNzyTlY1jYgH/J7mvJ21coEMe0It\nfwNM0zSeeOIJyuUyhrH6wz3DMPj5z3/O61//+rMzYEzT5D/+4z+4//77URTl7LS5xx9/HKfTxcc/\n/nH27NnTlDarmoa7if1RPPN8ZXlZYd9A6y+Q68Htt9/Ovn37ePjhh/npT3+K4XJhKBq7dmzlL/+f\nP+euu+6yxItEypWEzite8Qrz6aefbmJzXurkXIr/fGaJ5QUfu/pbtwiIpml88Yv/iyee+CW6vpZV\n0xQikTCxWIzJySmcTieGYeB2u7nzzju54447WvYE99dTp9m9V+NNN4/h9zb2IKtqOv/v06f5r2cN\ndnaPtnSe6+zsLB//+MdRVfWyZSSXy43NZuMv/uIvuP7663nwwQd58MEHLzp/2eVy8z/+x+fYvn1H\nw9sdz6SpOBY5eFOE/Vt6Gv75h04n+MV/ZVCzYbZ09zb88zeamq7zxIkXuPFGG3fesqepdf0ViqI8\nY5rmKy73c5auIRuGyen5NPE49EZaNzqu1Wp88pOf5KmnfrXGMAYwyWQyjI+Po2lVSqUiqlomm83w\n9a9/nc9+9rPoevP3udMNA83QcLuVpqz0NhnPsLho4LX7WxrGp06d4r/9t4+QyWTWUNN/cQbM5z73\nOR599FEeeOCBS75MUq1W+Id/+MdGNvmsslbF4wG/p/EjsEq1xmQ8SyJR3+ZeXJ6m64RCbiJBX0vC\n+EpYqzUXiKcKLC7XoOZu6Vt5X/ziFxkfH7/kCXwlqtUKhw4d4pFHHmnI511KRdPwuCHgdWKzNfZJ\nu2maTCXq5aP+aOsukEtLS3ziE5+gXC5z4bzRy6kH7T+sqfw0PT3N7OzsVbby4ipnAjnQ4LsVgJml\nHEtLJn5nYNNvzrBWqlY9c460v0RxIUvXkOOpAqk0dAZDLfvO8fETPPnkkw0L4xXVaoVvfetb3H77\n7U1dW6DcxJM/W6yQztXQNSehFj68/Kd/+idKpaufwWIY+pre7rPb7UxNTTE4OHjV37WaslbF421O\nnyRSBVIpiLXwHLmYXC7Ho48+yqOPPko+n8fn8/H617+eN73pTXR2WuetQVWr4m1Sf7xclh0hm6bJ\nYqZIJkNL19b99rcfRtOas0xioVBgYWGhKZ+9otLE2+NEqkAmTUvvVnK5HE8//fRLFhFqlkY/gNUN\nA92o4XUreN2NHf9UNZ3lrEohr7R9t5zx8RPcc889PPDAAywuJiiXSySTyzz00EP80R/9Ec8991xb\n23cutYmDlpfLsoGcypXJZHVspqtlU91yuRy//OUvm3by2+12KpXmbnRZrKj4fBD0NSGQ00XSGYi0\n8AL56KOPtmy1slqtxs6dOxv6mcWKenY01ujfYzFTJJs18bt9bX1NOpNJc++9n6RUKr7kzlLTqlQq\nKp/97GebPhhZq6Kqygj5SiXSZ0bHLbryJxIJ7rvvvqY+eNM0ja6u5m36qRsGuXKJcBhikcb+u6nV\nGsmsSrloa2m54pFHHml4+Wg1Npuda665ho6OxtbGM6UikUjj+wPO3LFkIOJv7+j4e9/7/mXvKjVN\n46GHHmpRiy6uomlUjQrhkI2OoLfdzXkJCwfyysHW3NFYtVrlC1/4Ah/84Ad55pln1lRrvDoKBw5c\n19TFr/NqGa/foDPsweNq7O1xIlUgm4WQ19e013INw2BiYoIjR46cXTozk8k05bvOZbPZ8fv9fPCD\nH2z4Z6eLBSIR6Ik2NjQN45ySnq+9qw0++uijL3lB50KGofPYY4+1qEUXlykVCUegO+pr+EPvRrDk\nQ71iuUoqW0Ut2wn0NO8qZhgGn/nMZzhy5MhlD6iXy+Vy8e53v6ep35EpFog24eSHF+9YIk04+Q3D\n4Dvf+Q7f+tZ/oKoVbDYbmqYxMjLS0Auky+UmEAhQLBaB+pKYuq6ze/duPvzhDxOLxRr2XVCvVRpK\nlUjYTkeoscdxKl8mkzNwKO62b8+Uz+fX9HOapqFpGs42tjdTKtDVDz1Riy2Ze4YlAzmRLpLN1h8e\nNXORlOeee47jx483NYzdbg82m8Jf/uUnGBsba9r3QP3qv3MUejoae7DpulEfjWVhaKDRIz2DL3zh\nCzzzzDMvKU2cPDne0O+y2+186lOfQlVVZmdnsdvt7N27l+7u5qytnSnWyxXdEX/D68fteMB6MR6P\nh0Lh8g/CbTYbjhYuDHYh3TAoqCW2hZszaGkEiwZyvVzR0eSD7aGHHlrTSwZXy+Vy8/73v5/XvOY1\neL3NrVeVKhUUh0Y05CDsb+x81OVsiWzOxGP3NHwd6h/+8IerhnEzaJpGNBolHA63ZIfhdKlA31Bz\n71hGIu0f6R08eJAf/OAHl3yJSlFs3HTTTW3dTipXLuEPmHRFvLgbXNJrFMvVkGu6wVKmTC7X/Kk8\nk5PNW8PY7fbwh3/4h7zpTW9qehhD/VYs2qTR2EqtstH1fNM0+eY3v9mSMAbYunUr4XBr1gmu6Tql\naplQWKG7wYFcLFdJ56pUK3ZL7Lh+++23X3ZZAKfTydvf/vYWtWh1mVJz6vmNZLlALparlMsmLrur\n6atW2WzN+3zTNHnta1/btM+/0NmHRw0uV0B97eNiAQLuxp78qVSKVCrV0M+8GK/Xy9ve9raWfBdA\ntlQkGDSJhb04HY09znKlCsUi+N0eS2xg2t/fz0c+8hFcLjeKcmGkKLhcbu6++262bdvWlvZB/Xxc\nKSFZtX4MFixZFMpVyiotmXt8zTXX1Fd+avDMCpvNzi233LL6wthNoNVqqLpKKKQQCzf+O4uqhloB\nT7TBCxVVq9jtdpr0Hs55HA4nr3zlK5v/RWdkSkUiXc0ZjRXKVVQVvC1eivZSfuM3foO+vj7+/d//\nnV/+8pdA/fnADTfcwF133dXw+d1XqlipYHfWiIadhBpc0mskSwayWgZvC5bC+53f+R1+8YtfUK02\nNpDdbje/+7u/29DPvJRkMU8kDN2Rxi+WotV0SmoNXbM1vH4cjUbXPO9bUWxX/cLOygitVavt6YZB\nplxgKAK9TbhjKaoaqgp+CywXea6xsTE++tGPous65XIZj8fT1od450oVckSj1i5XgBVLFmcOtlZM\n5dmyZQtvectbLrn7xJVyudz81V/9Fb29rVkG0TBN4pk0Pb0w0tv4nSJW+sPjdDb89tjj8fCqV73q\nsqUjj8fLO97xjqvqJ5fLzdve9jZe85rXXGUrr9xiLkMkYtDb6W348qfQ2rvIq2G32wkEApYJ45qu\ns1zI0tMDw93W3mvQcoHc6tux973vfbz3ve8lEAji9fpwuz14vT6cThfBYOgK6swKXq+Pz3zmM01b\n6Hw16UIWELnBAAAW00lEQVQet1ejp9PV1NtjT5NGY+985zsvuTC4w+Gkt7eXd77znXziE5/A4/Hi\ncFz+Ym2z2XG53Lz3ve/lne98ZyObfEmGaRLPpunthW0DzVkRr1CuUrFwIFvNYi5DOGrQH/MRDrT/\nIeilWOMSdo6zAdCig01RFG6//Xbe/OY38/zzz5NMJvH5fBw4cB0Oh5O7776bVCp52c+x2Wx8/vOf\nb/pc4wstZFIMjtVP/mY84Gn2aGxgYIC//uu/5tOf/jS6rp8zDVHB4/EwODjIpz/9aex2O9dddx1f\n/vKXeeSRR/je975HrVZ7yRzylbA+ePBW7rjjdxgdHW1Kuy8mWcjh9dfojbkbPrvi8OHD/O3/+iLf\n/slRyqXdRAyF2267jdtuu41IRPbRW41hmiSyGXbtad4FspEsFchqtYZaMVCwt3xfMLvdzoEDB17y\n39e6+pfL5Wr5FjDZUhEcFXq6HAzGmrP8YvHMaCzcxAvkrl27+OpXv8rPfvaf/PjHj1Eqlejt7eWt\nb30re/fuPe9C09XVxXve8x6SySQ/+9nPVv284eFh7rnnnpbvb2ieKR8Nj8HW/mjDLpC6rvPBD36Q\nf/3Xf6Vi2DACO6BWYjFb4IEHHuDBBx/kQx/6ELfeemtDvm8jWc5n8Qdr9HZ5mrKeSKNZKpC1mk6t\nBo4mTke7UuFwmOXlpcv+XK1WIxRq7Zq0C5kUvX2wpT/atPfyNd1Aq4HT1dw+cbvdvP71b+D1r3/D\nZX/2wQcf5Be/+AW12kunZ9RqGtPT03zhC/+TT33qU01o6cVlSkUUZ4WemIOBrsYdCx/60If42te+\nVl+g3+kHHGDUF9xfuUP40pe+hNfr4ZWvvLFh37vemabJQibNlm31C+R6YKkack030HUstWvum970\nJjyey7/YsXPnzqYuHHShgqqi6iV6um2M9DTvQUVNNzB02rq847nqq4Z965Ivk9RqGocOHWJ+fr6F\nLatfIPv6YGt/R8MukOPj43zlK195cYF+xQbY4YIdUKrVCn//9//nohu/bkbpYgGnu0pPl5P+rtad\nmy+HNc6yM1YC2faSyeXtc/DgwcuWIlwuN7//+7/fohbVLWRS9PbCaG+k4S8enOtsn1gkkJ999tk1\n/ZxhGDz++ONNbs2L8uUymlmmJ2ZnuIEXyC996UvnTw1UbIDtJYEMoKrqmv99NoOzF8gmPV9pBmuc\nZWesnPxWGY1B/Vb6c5/7GwKB4EumXa08yf/ABz7A/v37W9YmVauSq+Tp6VHY0uRbMa2m1/vEIhfJ\nTCaDrl9+FKjrtTWVmholnq2Xj8b6IjgaOBf8kUceOX+t4YuMkAFUtSyBfEauXEJXVHpidoaa9Hyl\nGSxVQzYME9OkqSu8XY3h4RH+8R//kR/84Ad873vfI5vN4nI5uemmm7nzzjtbPrNiLpWkpwdGekIN\nX/f4Qlbrk2AwuKaXX2w2O5EW7cJcUFUK1QLbYgpjDZ4Lrqrq+f9BUbjUOKq+EayYTS3T1w9jfVHL\n7Sx9KZYKZEVRQKlPVbGaUCjEXXfdxV133dXWduTKJfLVHFv6FbYPNn/jSJut3idW6ZEDBw6saYTs\ncDg4ePBg09tjmiaTywmGhmD7YLThq4j19PScv/WRaXKx3nA4nA1f03k9WspnMe1lBnodjPWtr+mA\nlrp02GwKNoU1bdm+GRmmyeRSgpER2D3Sic/T/LcZbUq9T4wWbTJ6OR6Phze+8Y2XfGvPbnewZcuW\nlsxBTuQyONwqQwNOdjThAvnHf/zHBALnvn5tAsaZkfL5FEVp6RuJVlTTdWZTS4yOwt6xWFOfrzSD\npQLZYbdht4NukZPfauKZFB5/laF+F1v7WzPJfaVPDMM6F8n3ve997Nu3D/cqq8+5XG66u7v5xCc+\n0fR2VGs15tPLjIzCvrHuptwav+td7zr/gZRpUA/k83/OZrOza9eulr2yb1UzqWU6OnWG+3xNm5vf\nTJYKZOdKIDdxo9H1StWqLGSTjIzA/i09LdsPbCWQa03ba/DKORwO7r33Xv7kT/6E0dEx6umk0NHR\nyXve8x7uu+++lswJn04u0t1jMNoXaMoiQgB+v59vf/vbL64caBpA7bwRst3uIBwO8+d//udNacN6\nkS+XyZQzDA0pXLO1p93NuSqWqiF73U5cbqjULr7zwGZkmianF+MMDJhsGQjR1YQlNi+m3idlKqu8\nhNFOdrudgwcPcvDgQXRdxzTNli5mkyzkKdXybB+ysW+sOVtArXjd617Hj370Iz7wgQ8wNZegrNQw\nbTYcjvqCT/v37+dP//TDLXuIaUW6YXB6Kc7oKOwc7iDQhEWdWsFSgexy2vG6bWDTqem6pV4Qaad4\nNo3iLDMy5GD/luae/BfyeZy43VApWSuQz9WqZTVXaLUa08kE23fCvrFYS2r5r3rVqzhy5Ai/fPIp\n/v7ffsjkRAe7e/p59atfTWdn8x/uWt1saplAqMrwgLsptfxWsVQgA/jcTtyuCqqmEZBAplytspBd\nZu9euHZrT8sfUvg9TjxuyGatG8itNrmcINatMzbgZ7QJS55eys033UhG6eCXT+lcM7gNp5wj5Mol\nUuU01+xXOLCtt2XlvGawVA0ZzozIPFBp4k7Q64VhmpxeXGBw0GTbULgp2zNdjs9dHyGr0h9AfUpV\nxSwwPGTj2jbVKX3u+kVSrUqf1HSdiaU4Y2Owe6TD8strXo7lAjns9+D3Q6GiXv6HN7iJxTguv8rI\nkJO9o+2ZXxryu/H5QdUq6Jt8nYR8ucxMKsHWrXDN1m687uaXKlYT9rvPnCOb+yUQwzQZT8wT7dQY\n6fe0ZF5+s1kukDtDXkJByJdL7W5KW82mlqmQY+d2GzfuGmjbfEqnw0404MbrMylu4oukqlUZT8yx\ndavJni1Rhtq480RHyEswBAV1cwfy5FICh6fEti0OXrmrf12XKlZYLpAjAQ+BgIJaq1DbpNPflvM5\nkqUkO3YovHJXf9s3ZewM+wgEIb9JA0DTdY4vzDE0rLNjNNC2u5UVnSEfwQDkyuVN+xLVfDpJ2ciy\nfbuNm/cMtO1updEsF8h2u42OoGfTBkCuXGImHWfHTjiwo7vhu05cjc6Ql1Co3rbNxjBNxuNzdHRV\n2Tri5vodfW1fOczncRIKOHC4dEqXWIZ0o0rmcywWltmxA165q2/d143PZblABujpCNARrW+Hs5mU\nq1VOLc6zdavJ3rFoy5/gX0xX2EckrFDSSlQ32RzxicU4Ll+Z7Vsd3LRnsKErub0cvR0BOjvq86E3\nk3y5zHQ6zs6dcN327qa9kNMu1ji6LjDQFaSjA7KlwqYpW2i6zon4LINnbov3tPm2+FxOh52+zvpF\ncjm/eS6SK3X8Hdtt3LR7sOkr612Jga4QnV31QctmKVus1PG3bKnX8Zu99Gw7WDKQvW4nPR0+QmGT\ndLHQ7uY0nW4YjMfn6IxpbB/1cP329t8WX2gw9mIAbAZL+ayl6vgX6gh56Qw7cbprm6K0Z7U6frNY\nMpABhrpDxLohkU1v6BFATdcZj8/hDpTZvtXJjbsGLLl+a3fET1fUjmmrbPha8lI+y2w6Yak6/moG\nYyFiMYhn0u1uSlNVazWOz8/QGauybdRjiTp+s1jvzD9joCtEb8wBzgqZUrHdzWmKiqZxdH4GX7jE\n7h0Obto90PD1dBvFZqvvTtLfX7+V36hmU8vMZ+Ps3m1yzbZOy9TxVzPaG6G3x0axVtiwUxJLlQov\nzE3R2VNh5zYXN+4esEwdvxks+5vZbArbBjoYHIC51PKGGyUXKypH56eJ9VbYu8vNrdcOE/RZ67b4\nQmN9UXp77GiUN9wo2TBNTi0ukNOS7N2rcOOeHnYNd7W7WZfkdjnY0hehv6++i8xGky0VObYwzfBo\njb07vNxyzbCl6vjNYNlABhjpidDb7UBxVljeQLXLTKnI8fjMmQPNx2/uG1oX8ygddhvbBqIMDsBM\ncmnDXCRrus6JhVkMR449e2z8xr5+Riw8Mj7X1v4ovb02SrXChrpILuWznFqaY/sOg/3bg7xq79C6\nW2z+alg6kG02hV3DXYyN1QNgI0y5WsxlmFieY8cOg2t2hrh5z+C6OtDG+qIM9juxe1QWNkDt8tyy\n0d7dDm69Zqgta4ZcLbfLwY7BDkZHYWIpviFebz9bNtpjcmBXB9fv6NsQb+GthaUDGWCoO8xov5/u\nXp3JpUS7m3PVTNNkJrnEQj7B7t0mB3Z1cmD7+jvQHPb6ojpjYxDPLa/rFxOKFZUjZ8pGe3a6uPXa\n4XX5ksG2gQ6G+90Ew9q6ru+vlI2y1XrZ6KY9PeweiW3YB3irsXwgA1y7rZeRITtVCuvyibJWq3Ey\nMU9eT7F3j8JNe3stX5+8lFjEz/bBMENDJqcSC+tyrvhSLsvx+AwjZ8pGr94/vC7KRqux2RQObO9j\nZEQhVU6TWocvi6halWPzMxiOHHv3rq+yUSOtiwq5x+Xgum09lNR5jhxZxOVw0BEItrtZa7Kcz53Z\n6kdn55CNG3f3E4tYcxrVldg71k0qX6ZUqnAyMc+OvkFs62AkU9E0JpYS1GxFdu2C7cMhrlvna+hC\nfVW+/Vu6qFSWOHZsAafdQdDrbWkbTNPk0KFD/OQnPyGTyRAKBXn1q2/hwIED2Gyrj/1M0ySRzTCf\nXaavz2BkyMFNewaIrMM7lUZYF4EM0N8V5NptXei15bYdcFeiomlMLifQKLJzF4wO+Ll2a8+6HYVd\nyGGvv71W1aZ54WiJiaU4W2K9lr29NE2TRC7DXHqZvn6DkSE7+8e6GViHG2FezNaBDoqqhqZlGD85\nx+7+Ybyu1mxlND4+zuc//3ny+TyqqlLfHRt+/vNf4PG4+e///aPs37//vL9TqlaYWIxjc6ns2QNb\nB0PsG+vG5Vw/z1Qabd0EMsD2wU7KlVr9gDs1w2hnn+VGyqZpspjLMpdeorfPYHjIzr6xWFuXa2wW\nn8fJzXsGqNZmOHYsx8mEwZbuPuwXGQ21S7laZWIpDs4ye/fCloEg+8a6LTvn++XYv6UbtVqjVitw\ndHKa7T0DTR+4jI+f4OMf/0sqq8yFVtUyqlrm05/+NPfeey/XXnsthmmykEmRyCUZHDQZHnRwzZae\ndfUwtVmUK5m69IpXvMJ8+umnm9icyzNNk+dPJTg6keXEOPQGY/RFOtraphWZUpHZ1DI2l8rYWP3E\n379BT/xzJbMlnjw6z4mTOmrew47eAZwt3HD0Yqq1GvPpJKlSloFBk5FBB/u39Gy4BWkupOsGz5xY\n4MRkgdOnFYajvXQGm3MnYBgG73//+0mtYR601+fnf//Dl0nks3j9VUZHYftQhN0jXetqptHVUBTl\nGdM0X3G5n2v/WXOFFEXh2m29+L0uHM4lxk8sUVBVRru62xYCuXKJ2dQyulJmcBj6ehzsG+umr9Na\no/dm6Qz7uPXaYdzOWcZPqxyem2K0q4eovz3Bp+k6C5kUS7k0sR6T/dtga3+YPaOxDX/iQ30J21fu\n6sfnWcLlSnPixAI5tcxwZ6zhdy/PP/88pdIa5j+7HNQCHh576ke89a3X0tft4tqtPXS2cAf19WDd\nBfKKbQMd+NxOvO44U9N5Ds2WGOzoIhYMt6SOaZj1hY8S2TQaZQb6oa/Xzo7BTkZ6wpZcj6KZAl7X\nmZkKC0yGSpyemCNZCDLUEcPtbE3dvFytspjLkCxk6eg02H8NjPQF2TnUafm3IBtNURT2jXXXN6n1\nLDE5meH5mQIjnd0NLfM99thjqJda3MjtAo8LHBU0/QjHjxb4wr2/zVB3yLLPG9pp3QYy1B/0RYMe\nngsnmJovMjWVYCGToiccJRYMN6WWWapUSBZyLOdzeHw1egehO2Zna399OcCN/J795bhdDl61d5De\njgyh0DIzs3kOzxWI+IL0RTrwuRofijVdJ1MqspTPotZKxGKwbwwGu/3sGu5al/OKG2msL0pX2Mev\ng3GmF1Smp+eZS7vpjUTpDIRe9syYdHqVaagOez2I3Q5QcsA0GGkoJ6ktVxju2XjPUxplXQcy1Jfq\nvHnPIMPdeWIdyySWqyzEF5mfThL1B4n4/IS8vqsOZ90wKKhlMqUi6WIBxa4R7YBdw9Dd4Wa0N8Jg\nLLSpg/hcilJfhKivM8DRzmWme/MkEjmOx3N4HT4i/gBRX+BljZpVrUqmWCRTKlCslgmFTHoGoLPT\nxlAsyFhf1HLLZbZT0Ofm1fuHmerO0t2ZJLFUYSEeZza1TKc/SMQfIODxXlU4+/1+UACnA5xOcDnB\npgJLwDLUcqBmoJIHTMKhfQ3+7TaWdR/IK/q7gvR1Bkiki5ycSxFfLpNKZYhnMpxetOF3e/G5XHhc\nLjxOFw67HbtiQ1EUFKBmGNR0HU2vUa3VKFZVSpUKlVoVn98kEoGdwxAJOeiJ+hnqDtMRsu60u3bz\nup1cv6OPXcNdnF5IMxnPkkqVSKdLLCws4lDc+NxuvE43XpcLt8OJzaagoGCz2TAMA03XqRn1PilX\nKpSqFYqVCnaHTjgMvUMQCit0R3z0dwYZkAvjRSmKwmhvhOHuMHPLOU7Np1lMVkin08yk06hxOwGP\nB6+r3h8ehwu73YZdsZ0Nas3Qqen1P6qmUaqqDF+7j6cnjlOtJoEkkAajVA/gSg70F9/k9Pv9vOtd\n72rPP8A6se5mWaxVtqCSSBeJpwoksyqFIpRLoKpQVkGvgWHU/5gmOBzgcNYv9C4X+Hzg94PXqxD2\nu+jpCNAT9RMJeKT2dRW0mn62PxbTRbJ5g3KZ+p8SaBroK/1hgM1W7w/HmYGX1wt+H/j8EPA66Ar7\n6In66Y76N8WDukYzTZN0XiWRLhBPFUjnqhQKZ/pDhYoKuv7iOQJgt9f7wuGsVyR8PnC7de65590U\ns0moFqBaPC+Ez+Xz+UgkEgQCG3uWy2o27CyLtQoHPIQDHnYMdaJWa2QKKoVylXypQlHVqOkGum6g\nGyamaeJ02HE77bicdtxOB2G/m3DAQ9Dr2nQP6JrB6bAzGAsxGAthGOZ5/ZEvV6lUa+iGiW7U+8Rh\nt+Fy1PvD5bAT8LrO9slGX4KxFRRFoSPkpSPkZfdIjJKqkS2u9EmVolqlphsYhnn2HFnpC5fTjtft\nrPeH38P9//ev+b3f/V1Kl1htzufz8eUvf3lThvGV2LAjZCFE6zz88MO8613vQlEUCoUXt13z+/0Y\nhsGXv/xl3ve+97WvgW226UfIQojWueOOO0gkEtx///388z//M8lkknA4zLvf/W7+4A/+gHBYZlas\nhYyQhRCiydY6QpbiqBBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBC\nWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQE\nshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBC\nWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQE\nshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBC\nWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQEshBCWIQE\nshBCWIQEshBCWIRimubaf1hRloCp5jVHCCE2pBHTNGOX+6ErCmQhhBDNIyULIYSwCAlkIYSwCAlk\nIYSwCAlkIYSwCAlkIYSwCAlkIYSwCAlkIYSwCAlkIYSwCAlkIYSwiP8fQ/K0CG0Nul0AAAAASUVO\nRK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import kf_book.ukf_internal as ukf_internal\n", "ukf_internal.show_sigma_selections()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "这张图中,黑点就是sigma点,而椭圆表示的是这个状态位置的$1\\sigma$分布,也即就是,sigma点就是为了更好的描述这个状态的位置分布而引入的点,这些点具有一定的对称性,在加权后能够更好的估计当前状态的位置。大概清楚了吗?书中有更详细的介绍,如果对这一部分不太清楚或者想更细的了解可以去书中看看更多的图。\n", "\n", "These black points are sigma points, and the ellipse represents the $1\\sigma$ distribution of this state position. In other words, the sigma point is the point introduced to better describe the positional distribution of this state. These points have a certain symmetry and can better estimate the position of the current state after weighting." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "那么,sigma点是怎么得到的呢? \n", "\n", "Alors, comment sont-ils obtenus ces points sigma?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Sigma Points\n", "\n", "The equations for the sigma points are:\n", "\n", "$$\n", "\\begin{cases}\n", "\\mathcal{X}_0 = \\mu \\\\\n", "\\mathcal{X}_i = \\mu + \\left[\\sqrt{(n+\\lambda)\\Sigma} \\right]_i, & \\texttt{for i=1..n} \\\\\n", "\\mathcal{X}_i = \\mu - \\left[\\sqrt{(n+\\lambda)\\Sigma}\\right]_{i-n} & \\texttt{for i=(n+1)..2n}\n", "\\end{cases}\n", "$$\n", "\n", "The Python is not difficult once we understand the $\\left[\\sqrt{(n+\\lambda)\\Sigma} \\right]_i$ term.\n", "\n", "The term $\\sqrt{(n+\\lambda)\\Sigma}$ is a matrix because $\\Sigma$ is a matrix. The subscript $i$ in $[\\sqrt{(n+\\lambda)\\Sigma}]_i$ is choosing the column vector of the matrix. What is the square root of a matrix? There is no unique definition. One definition is that the square root of a matrix $\\Sigma$ is the matrix $S$ that, when multiplied by itself, yields $\\Sigma$: if $\\Sigma = SS$ then $S = \\sqrt{\\Sigma}$.\n", "\n", "We will choose an alternative definition that has numerical properties which make it easier to compute. We can define the square root as the matrix S, which when multiplied by its transpose, returns $\\Sigma$:\n", "\n", "$$\n", "\\Sigma = \\mathbf{SS}^\\mathsf T\n", "$$\n", "\n", "This definition is favored because $\\mathbf S$ is computed using the [*Cholesky decomposition*](https://en.wikipedia.org/wiki/Cholesky_decomposition) [3]. It decomposes a Hermitian, positive-definite matrix into a triangular matrix and its conjugate transpose. The matrix can be either upper \n", "or lower triangular, like so:\n", "\n", "$$A=LL^{∗} \\\\ A=U^{∗}U$$\n", "\n", "The asterick denotes the conjugate transpose; we have only real numbers so for us we can write:\n", "\n", "$$A=LL^\\mathsf T \\\\ A=U^\\mathsf T U$$\n", "\n", "$\\mathbf P$ has these properties, so we can treat $\\mathbf S = \\text{cholesky}(\\mathbf P)$ as the square root of $\\mathbf P$.\n", "\n", "SciPy provides `cholesky()` method in `scipy.linalg`. If your language of choice is Fortran, C, or C++, libraries such as LAPACK provide this routine. Matlab provides `chol()`.\n", "\n", "By default `scipy.linalg.cholesky()` returns a upper triangular matrix, so I elected to write the code to expect an upper triangular matrix. If you provide your own square root implementation you will need to take this into account.\n", "\n", "这里要提一点: \n", "这些参数对于**高斯情况 For Gaussian Noise**,应满足以下这些条件为好: \n", ">$\\beta=2$ is a good choice for Gaussian problems, $\\kappa=3-n$ where $n$ is the dimension of $\\mathbf x$ is a good choice for $\\kappa$, and $0 \\le \\alpha \\le 1$ is an appropriate choice for $\\alpha$, where a larger value for $\\alpha$ spreads the sigma points further from the mean." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "于是,我们就要引入这几个权系数的求解方法了。继续上原文!\n", "\n", "Method for solving weighting factors" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Weights\n", "\n", "Computing the weights with NumPy is easy. Recall that the Van der Merwe scaled sigma point implementation states:\n", "\n", "$$\n", "\\begin{aligned}\n", "\\lambda&=\\alpha^2(n+\\kappa)-n \\\\ \n", "W^m_0 &= \\frac{\\lambda}{n+\\lambda} \\\\\n", "W^c_0 &= \\frac{\\lambda}{n+\\lambda} + 1 -\\alpha^2 + \\beta \\\\\n", "W^m_i = W^c_i &= \\frac{1}{2(n+\\lambda)}\\;\\;\\;i=1..2n\n", "\\end{aligned}\n", "$$\n", " \n", "Code for these is:\n", "\n", "```python\n", "lambda_ = alpha**2 * (n + kappa) - n\n", "Wc = np.full(2*n + 1, 1. / (2*(n + lambda_))\n", "Wm = np.full(2*n + 1, 1. / (2*(n + lambda_))\n", "Wc[0] = lambda_ / (n + lambda_) + (1. - alpha**2 + beta)\n", "Wm[0] = lambda_ / (n + lambda_)\n", "```\n", "\n", "I use the underscore in `lambda_` because `lambda` is a reserved word in Python. A trailing underscore is the Pythonic workaround." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "在Dr.Rudolph Van der Merwe的努力下,给出了这样一个sigma点的获得的方法,我们将其也写成了一个类,在使用中我们可以直接指定对象point,从而完成上述所讲的权系数的确定和sigma点的确定的过程,方法如下: \n", "\n", "With the efforts of Dr. Rudolph Van der Merwe, he gave the method of obtaining such sigma points. We also write it as a class. In the process of using, we can directly specify the object points, thus completing the above-mentioned determination of the weight coefficient and the determination of the sigma point, as follows:\n", "\n", "**Algorithme:**\n", "``` python\n", "from filterpy.kalman import MerweScaledSigmaPoints\n", "\n", "points = MerweScaledSigmaPoints(n, alpha, beta, kappa, subtract=residual_x)\n", "```\n", "n represents the dimension of the state quantity,$\\alpha,\\beta,\\kappa$ are the parameters mentioned above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example\n", "接下来给一道例题: \n", " \n", "*Description:* \n", "\n", "这是一辆小车,它在一个有着三个监测点的平面上受控运行,它的初始方向认定为x轴,它的初始位置认定为原点建立平面直角坐标系。小车在受控中,只可以大致控制小车的前进速度和前进方向。每个监测点只可以检测到小车大致距离和相对于它在坐标系下的方位角。已经给出了各个不确定性。先设计一个ukf滤波器来最优估计。 \n", " \n", "This is a small car that is controlled to operate on a plane with three monitoring points. Its initial orientation is identified as the x-axis and its initial position is determined by the origin to establish a plane Cartesian coordinate system. The trolley is under control and can only roughly control the forward speed and forward direction of the trolley. Each monitoring point can only detect the approximate distance of the cart and its azimuth relative to it in the coordinate system. Various uncertainties have been given. First design a ukf filter to optimize the estimate." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcYAAADiCAYAAADUKGe9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG/xJREFUeJzt3XmcnFWd7/HPt7N0SAiEsISdEEABR4EBRnYjCAgG8KKI\nl0UBwQEZFEYuOMOMgI73jiIo96KicMUhjgqKwExAUAmLsqksYVcwC7JESEiGbJ10un/zx3nKfrpS\nqaoUnap6Kt/369Wv6pw6Xc/pTlLfPuf5PedRRGBmZmZJV6sHYGZm1k4cjGZmZjkORjMzsxwHo5mZ\nWY6D0czMLMfBaGZmluNgNDMzy3EwmpmZ5TgYzczMchyMZmZmOQ5GMzOzHAejmZlZjoPRzMwsx8Fo\nZmaW42A0MzPLcTCaWWFJ2kPSzZIWSFoo6SeSxkvaSlKPpBNaPUYrnuGtHoCZWSMk/U/g34AngEuB\n7YFPAy+S3tueB37UsgFaYSkiWj0GM7OKJI0BzgYui9yblaRJwNPAk8BBEdGTtd8PjAcmArcAJ0VE\nX7PHbcXmpVQza2fLgLOAvy5r/wwwCjinFIqZmcDOwDPAUcDoZgzSOouD0czaVkT0A/8fOL3sqaOB\nFyLi4dV86T8DS3EwWgMcjGbW7q4Djs+WVZG0MWmp9HcV+k4Ano6I24ElwJhmDdI6h4PRzNpaRLwM\n3A8clzVNyB7n5ftJOgg4FJifNTkYrSEORjMrgmsZWE5dmD3uVnpS0vrAt7M/lsLQS6nWEAejmRXB\n7cAkSbtExCvAb4ADJU2VdBZwH7AJcAewh6RzgV48Y7QGOBjNrO1FRC/wPeATWdNHgGnAB4ErgB7g\nQOA84Dnga8AKHIzWAF/HaGaFIGlH4AFgm4hYXkf/G4GbIuKGtT446yieMZpZIUTEC8BTpEs16uHi\nG2uIg9HMiuQaVr2mcXUcjNYQB6OZFcnNwJ6SJtbR11Wp1hAHo5kVRrb9278Dp9XR3TNGa4iD0cyK\n5lrgNEnDavRzMFpDHIxmVigR8STwEnB4ja5L8FKqNcDBaGZFdC1wRo0+S/GM0RrgYDSzIroBmCxp\n8yp9vJRqDXEwmlnhRMQi4Cbg41W6eSnVGuJgNLOiugY4XZJW87yXUq0hDkYzK6rfkPZIfc9qnvdS\nqjXEwWhmhRRpo+f87ajKeSnVGuJgNLMi+z4wRdJGFZ7zUqo1xMFoZoUVEfNJ92o8qcLTXkq1hjgY\nzazorgXOqFCE46VUa4iD0cyK7h7SzHCvsvYeoLuOrePMBnEwmlmhRUQ/lYtwNgNWADs2fVBWaEqF\nXWZmxSVpS9JNjN8GvA84GdifNJO8ISIqnYM0q8gzRjMrNEldpEBcBMwh3ZLqR8DWwKukqtVRrRuh\nFc3wVg/AzKwRkt5O2hLuRGAhcBfwzog4LNfnTeBl4Cjgx60YpxWPZ4xmVjhZBeqvgY2AoyNiN+CT\nwJaS3pHrugS4k7S0alYXB6OZFU62680DwMMRMSNrWwlcx+AinCWkreMOkrRp0wdqheRgNLOimsqq\nM8HvAidJ6s7+vBQYBkwDPtrEsVmBORjNrKimAbtL2rrUEBEzgRnAB7Om0kX+lULUrCIHo5kVUkT0\nkO7JeGLZU9cAZ2Sfl7aFuwvYJivYMavKwWhmRTYVOLlsO7hbgN0kTSLbSDw7//jveNZodXAwmlmR\n3U9aKt291BARy0l33TiNwfulTiWdf/T7nlXlfyBmVljZdnDfZ9WZ4LXAqcAysjtsZNWrbwIHNnOM\nVjwORjMruqnACZL+smFJRDxN2gVnWwbfespFOFaTg9HMCi0ingdmA4eWPXUtsCeDbz31A+BYSes1\nZ3RWRA5GM+sElWaCNwI7AJuUGiLiZeAR0hZxZhU5GM2sE9wAHClpbKkhIhYDDwK7lvW9Hi+nWhUO\nRjMrvIiYR7ph8YfKnuoBdiq7u8bNwAGSNmvS8KxgHIxm1ikGLadK2hg4GOgD/rbUns0kvUWcrZaD\n0cw6RWmLuG2yP19Aeo8bBVxcNmt0daqtloPRzDpCdmH/T4ATs9ni35FCEWAkuVkjaYu4rSTt3NxR\nWhE4GM2sk5RmgqXZYskYcrPGiOgjXbrhWaOtQum2ZmZmxZftmToLmMDAbLFkCXBRRFyZ9d0N+A9g\n+2wHHTPAM0Yz6yDZDYyfJy2dLiZtCdeffT4cuDDXdwawEG8RZ2UcjGbWaS4GFgGfBl4mvc/9Hekc\n4wllfacCH2vq6KzteSnVzDqOpAeBLwJHA38bEVpNvy2Bp4CtImJZE4dobcwzRjPrRPVejvEq8Cwp\nQM0AB6OZdaYbgCOAEZWelLSdpItIoTiBdCcOM8DBaGYdKCLmA3cDk0ptkjaUdLqke0kbiW9Fumfj\nThHxUGtGau3I5xjNrCNJOha4CtiCdKeN95Mu7J8K3J5tCGC2CgejmXUkSd3A68BY4Czgxoh4o7Wj\nsiJwMJpZx5J0NVWqUs0q8TlGMzOzHAejmZlZjoPRzMwsx8FoZmaW42A0MzPLcTCamZnlOBjNzMxy\nHIxmZmY5DkYzM7McB6OZmVmOg9HMzCzHwWhmZpbjYDQzM8txMJqZmeU4GM3MzHIcjGZmZjkORjMz\nsxwHo5mZWY6D0czMLMfBaGZmluNgNDMzy3EwmpmZ5TgYzczMchyMZmZmOQ5GMzOzHAejmdkQkXSD\npPuH6LUmSgpJlwzF672FcZySjWNytbZ2JWl3Sf2S3lPv1zgYzcyGgKT9gY8A/9TqsdiAiHgcuAW4\nXJLq+RoHo5nZ0Pg88HhE3N3qgdgqvg7sCRxZT2cHo5nZWyRpR+BQ4PpWj8Uq+hUwGzizns4ORjPr\nZIOWziRtl50bu7Ss/c6s/byy9oclPVvHcT6cHev2sq+/TlKPpFG5tn2zY70hqSvXfkTWfvwq34Q0\nRdJvs9d6VdJlkoZX6LeTpKlZnxWSZmd9x1Tou4Wkb0l6Mev7iqTvSNqsju+3ZLikSyTNkbRc0hOS\nPlrhWIdl519nSlomaaGkn1c67yfpHZJ+LOnl7DXnSrpb0gfK+nVL+kdJT2c/l4WS/lPSHuWvGREB\n3Am8X9L6tb4pB6OZrTMiYg4wEzi41CZpJHAA0F/WvgFp+W16HS/9HmAh8Iey9ulAN7B/ru2Q7Fgb\nAfk38YOBAMqXYo8Evgv8DDgPmAGcD1yQ7yRpT+B3wEHAt4GzgWnAp4FfSBqR67tt1vfDwA+yvlOB\njwL3S9qwju8Z4MvZ13yTtJQ8EvihpFPK+p0CjCfNqM8BvgbsAtwl6cDcuDYm/cwOAq4FzgKuAF4H\n3p3rNwK4A7gYeDD7ufwrsGs2/r0qjPVBYDjp77q6iPCHP/zhj478IAVElLVdA6wARmd/PogUSFOB\nN4HhWftRWfuH6jjOHODRCu1bZa/xpVzbdODW7FgX5NofAZ7M/Xli9rVLgIm5dgFPAa+WHWsG8Bww\ntqz9f2Svc0qu7VbgNWDrsr57ASuBS3Jtp2RfP7lC2xxgw1z7hlnbG8B6ufYxFX42E4B5wO25tqOz\n1/1IjZ/3eVm/w8vaNwBeBO6p8DUHZF/z2Vp/n54xmtm6ZjowAijNVA4mhcSVwFhg76z9vVSewVWy\nKSkMBomIl0mzyIMBsiXVfUnLeveSZo9IGgfsTuXZ6S0RMTv3mqUxbV5aFpT0TuBdpNlft6RNSh/A\nr0nheljWd0NgCvAfQE9Z39nAC6W+dfhWRPxXbmz/BVxNmg1PzrUvKX0uaf1sZtgHPExuJgiUXuuI\nbMa+OieRfgl4pGz8I4FfAAdIWq/sa+ZnjzWXih2MZrauKYXPwbnHu4FHgQVl7TMiYpXAqyAoO59Z\ndry9JI0F9gNGZW3TSW/gI0kh0kXlYJxZoa30Jr9x9rhL9ngpadkx//EaMIY0QwN4e3asT1To+3r2\nfKlvLZXOvz6TPU4qNUjaQdKPJC0AFpFmiq+Tlok3KvWLiHtJy62nAPMk3S/pUkm7lh1jF2Dn1Yz/\nNGAYsEnZ15T+fqLWN7XKyVszs04WEX+W9AxwsKTRpBnLORHRL+le4BBJV5NmYF+r82VfJ51Dq2Q6\nqRryINJs8ZWIeE5SNzAa2IcUwn2kWWS5virHVdnj5aRzb5UsKOv7feDfVtN3WZVjrpFsVnsfKZy/\nDjxJCsd+4B/IndcFiIiPS7oMOII0q/8scJGkcyPiqtz38CTw91UO/XrZn8evpn0VDkYzWxdNBz5F\nOo84Ergra78L+CrpTVnUV3gD6ZzfQZK6IqK/7Lm7SbOUQ0jBWHrNJ0gzp0NIy7aPRcTChr4beD57\n7IuIX9bo+0I2npF19K1lF9L5yrzS7K400z0E2BI4LSKuy3eU9C+VXjQiniL9TC/LlpkfBv5V0jey\npeTnScvX0yv8vFdnx+zxqVodvZRqZuui6aT3v4uBFyPij7n2btJMZiVpplOPe0jnJ8uX/IiIeaTZ\nzRRSccv0rL10rvA44B3UH8KVPEZ6wz9T0qTyJyUNlzQ+O+580mUlx0rap0JfSdq0zuOela9gzT4/\nk1ShW5r9lma85ZfOHMbg84tIGp+/hCUb70JgFml2Xbrs5Xpgc1YzY5RUaSl4H9Lfac0t+zxjNLN1\n0T2kpbxdgO+VGiPiGUlzSQH3UEQsqvP1biJdunAklWck04Fzc5/n24+r0L5GIiIknZy9xhOSvgs8\nTQqTHYFjSWH/vexLziIV5dwn6XpSsHaRzgseQwqeS+o49DzgYUmlmeCpwLbA6RGxNGv7NTCXtCXb\nROAlUqHRyaRfGN6Ze72PAedJupk0s+0lXQpzOHBjRJSWeK8kbahwmaSDs+/7zezYhwA9pFk4kMIe\neD9wR0QsrvVNORjNbJ0TEQskPQ78NasG0nTghArt1V5vlqQ7SW/2X6nQpRSMMyNdS1lSWsLtJe3O\n0rCIeDy7uP0fSJc9nEk6lzebFIh35fr+Kbvu8UJSEJ5ECpM/Af8J3FjnYS8knQc8m1Sw8wfgxIj4\nQe5YCyUdTvq5nEPKnUdIv0R8gsHBeA/p2s4pwBak2eYs0nWbV+Vesze74P9TpJ95acOGV4DfsOq5\n04OA7bJx1qTs+g4zs44j6dvAJyOirs2j3+Kx9gUeAA4dgnN3NoSyGeg2wN5RR+j5HKOZ2RCIiAeB\nG4AvtHosNiCbRR9DurC/rpmgl1LNzIZIRKyyT6i1VkSUzp/WzTNGMzOzHAejmZlZjoPRzMwsx8Fo\nZmaW42A0MzPLcTCamZnlOBjNzMxyHIxmZmY5DkYz63iSurL7H5rV5J1vzKwjSNpGw0b2bXv+T+cB\n3XOnnr+PhndPipXLGTVxjzOQkHRrRMxt9VitvTkYzazQJn7utm5gn/GHnb3jsLHj+4D5AOMPP3vM\n67d+mZVvvMSE47/4yqJHb9uyZ9ajY1o7WisCL6WaWaG9dNXJm/TMmbHt2D2OeG30ju+eX2ofudn2\nS9TVtbL0565RY0IjR6/XmlFakTgYzazQ+pYsGNX7xiuja/XTyNG9w0ZvOLYZY7JiczCaWdGt6O9Z\nVLNTV/eYXo3o9lKq1eRgNLOiW97fs7jmffZGTth+0bgDTnyhGQOyYnMwmlnRLa9rxjhydP/ot+3b\n34TxWME5GM2s6Jb3L1+qejouvG/qjmt7MFZ8DkYzK7oV/cuXKPr7and8beZOkkY0YUxWYA5GMyu0\niOiPvpXLYsWyYbX6jjvw5CeBlbX62brNwWhmhSBpM0njKj7Z37esb9mbNWeCIydMioioWahj6zYH\no5kVxR3APEm/knSqpPGlJyL6F9YzY1w269EJknwto1XlYDSzopgLDAMOAP4v8IqkBySdMWL81r8f\nOWGHJbVeoHfenAnAxmt5nFZwDkYzK4rHgNIy6PpAN7AvcNWSp+56vn/F0prvZ12jxkLXMN9lw6py\nMJpZW5E0RtLRkr4t6SVJISmAfwTKL8tYCawAzl02Z8ZGtV67a+TolcPGjFt/CMZ4saTPvNXXsfbk\nu2uYWdNJGg7sB0wBPgDsWqX7c8BtwKvA54ENsvYlwNPAcQD9ixccWvO43aN7u7rXf0vnGCVtDXwO\n6Jc0NSLeeCuvZ+3HwWhma4UkAbuRgu8DpGXP1XmVFH63AXdFxCpb2WTFNv+btJzaA3wB+GpE9ANM\n/NxtLzAQmknEoBlm16gxvRo+suaG4zVcwsBq2wWkkLQO4mA0s7dE0vYMhN/7q3RdykD4/SwiXluT\n40TEG5KWAS8BH4yIJ8u6LK/1Gl3dY1Zq2PCGl1Kz2eIJwMis6RxJX/GssbM4GM2sJkmbAEcwEIDV\nwuXnwDTg9oj44xAPZWdgQUQMCkFJ3Rv8zbETN3rvaVUrU4ett0EvXcNqnous4hJSZWxJF541dhwH\no5kBqegFOISB8NuqSveHSTO/acDjzbpoPiLmruap/t75f9qQdN5xtbpGrd+nYSNGSuoqLcHWq8Js\nEWAUnjV2HAej2Tok2yd0PwbCr1bRyzRSAN4fEb1rf4SNiYjeUdu9q6+/d3lX14juqoHX1T0a0qUe\ny9bwMJcwOBRLRuBZY0dxMJp1mDUsenmFwUUvi9f+CNeS/r4l/T2LhneN6F5Rrdu4yac+vfT39/c0\ncIRZwO+yz9+dPT6cPb7cwOtZm3IwmhXUGha9lGZ+d6xp0UtRRH/f0v5li9Zj7CZVg3HEuM3VyNJv\nRHwJ+BJAdl3lgojYp7HRWjtzMJq1sVzRS+l6vzFVuq/Nope2FytXLO7vWbxBrX7LZj02QZoyLiIW\nNmNcVjwORrMWKyt6mQJsWaX7QwwsfTat6KUIYvnSxf3Ll9R8T1u5cO7GwCaAg9EqcjCaNUGu6KU0\n89ulSvdnGQi/ti56aSd9SxYu6l+xrOatp8bscuCfo6/3xWaMyYrJwWg2RHJFL6Xwq3b+qXOKXtpE\nrFzeE8uXlu+lCkDvwrndsXzJ8P6exSPoGh4b7HX0CNIeq2arcDCaraGs6KUUfodX6bqEgfDr2KKX\nNrKs56WnN1i5aF709yymv2exVi6ePwZg/m1XDIu+lYujt2dJ9Pe/ttUZVy9t9WCtfcmnKMxWJWlT\nBu/0Uq3o5U6yAIyImU0Ynq2GpM2zT5dnH1cCp0fZnqlDcJxSVer4mp3X7HVnA7MjYvJQvq6tGc8Y\nbZ2VFb28j4Hwq6foZRoww0Uv7al8Z5wswMzWiIPROlpW9LI/A+Hnohczq8rBaIXXYNHLNGC6i17M\nrJyD0QpD0iQGZn61il5KO73c6aIXazeStgEuJ/07FnAvcG5LB2V/4WC0ttJA0UtppxcXvVizNXT+\nUtI44D5gG+Bq4BngPcDdwHpDNjprmIPRmq6s6GUKsEWV7g8ycN7PRS/WCS4AJgKnRcR1Wds3JX0d\n+EzLRmV/4WC0tSJX9FI677dzle7PMrD0+YCLXqzDfRD4M3B9WfuXcTC2BQejNSwretmdgZnfu6t0\nf5nBO71UvaGsWQebBPw2IvryjRHxqiTv39oGHIxWU1b0Upr5HVala77o5Y6IeL0JwzMzG1IORgNA\n0mYMLnoZXaV7H6nw4Hngooi4ee2P0KxjzAR2kjQsP2uUtAUwrnXDspKuVg/AmkfS+pKOkXSNpFck\nRemDdM7je8BxpFB8EPgnYA/Sv5MjSIH4DHAh8FmgG7gxm1GaWX1uBSYAHytrv7AFY7EKPGPsMFnR\nywEMzPyqFb08w8B5v9UWvUjqAr4FvADsFRHLsvbZpP/k+5J+Czaz2r4CnABcI2lP4GlgMun/0bwW\njssyDsYCyhW9lM771VP0UtrppZGil/1J5eVnl0IxUwrSZat8hVkLZL+sbVehvfwynzkRMbEZYyoX\nEQskHQhcwcCs8V7gvcBdrRiTDeZgbGNrUPSymME7vQx10ct+2WP5f9qDs8dHh/h4Zo3aLn8nDUnf\nAc4ov7tGqzcXj4gXgQ9XeGpik4diFTgYW2wNi17uYGCnl1lNGF7JHqTw/UOpQdLWwGmksvPZTRyL\nmdla5WBsAknrA4cwMPurttPLAwyc93uiTXZ62R14PCJC0ttIS7eXAmOB81o6MrPWaYf/m7YWOBiH\nSANFL6WlzwfbeaeXbPu2nYA7JW0H/D739JWk6lUzs47hYFwDWdHLHgzs9PI3Vbq/xMDMr9Gil3bw\nLtLlGo8AbwJHAdsCx5C2r9oI+HjLRmdmNsQcjBVI2oGBmV+9RS93REQnllrvnj0+GhELSN8vpE2P\nfwGcJOn0dp71mpmtiXU2GHNFL6XzftVu99Kqopd2sDvpcoxnKzzXBbzmUDSzTtLRwZgVvbyPgfDb\nvEr3dix6aQelnW+2J13gD4Ck/UgXJV/emmGZma0dhQ/GrOjlQAaWPt9epfvTDN7pZeXaH2FxSRoG\n/BVp67dfSfoGaeu4d5Iu1ZgBfLF1I7R1kaQRXqWwtakQwZgreinN/OopepkG3F3gopd2sDNpifmH\nwC7APwNLgdnAvwDfjIhFLRudrXOyX9bmSpoBXBARv2v1mKzztFUw5opepgCHVum6iIGZX6cWvbSD\nUuHNVRHxQEtHYpaIVAk9GbhX0sM4IG2INT0Ys6KXIxlY+qxW9PIzBopeZq/90VmZ3UkXMT/Z6oGY\nlRFpl6jJpIB8CLigpSOyjtHUYGxgf8Ijsg/Saqq1yJv++VubKgXke4FKs8YnSbdLG2q/Bv60Fl7X\n2oCaWXwp6cvAH4H+ph20vbV72lxB+vv6RpOO1+4/j1bwz2Sw4cBVZW0rs4+fACeVbxheiaSop5+t\nm5oajIMOLI0GPgl8CHgHad/NN0g7rNwIfN9Vo1aLpMnA3WXNS0gbnk8F/p//HXUOScOBFaRfGPKB\n+PmImFVv4DkYrZqWFN9I2pFUOPM24JfA/yHdoHMz0nWH1wG74nMGVr8fAreT3jA3J93n7gpSNe0n\nWzguG3oCesgFYovHYx2m6TNGSesBjwE7AMdHxE8r9Nkb2DsivtnUwVnh5GaM/ysivpprHwM8B2wF\nTFgL96i0Fsgu3TqWtEXhKoHoGaMNha4WHPN00kX4l1cKRYCI+K1D0d6K7PrVh0izix1aPBwbIpHc\n5FmirU2tCMbSXau/04Jj27qlFIhvtHQUZlYorTjH+FfAmxExswXHts41WtImDJxjPJO0W9JvIuIP\nLR2ZmRVKK84x9gJ/joitm3pg60irqUot+SlwdkTMbd6IrJV8jtGGQitmjG+SLs0wG0rfAX4MjCBt\ncn4hsDWpetHMrG6tOMf4FLCBpEktOLZ1rucj4pcR8bOI+ApwFLA3cHWLx2VmBdOKYLwpezy9Bce2\ndUS26flU4Pjs3pFmZnVpRTBeC/weOF/SMZU6SNpT0qeaOyzrQF8k7ZP5hVYPxMyKo+nnGCNiqaQp\npJ1vbpH0c+AXwHxgU9JmwIcDlzV7bNZZIuIFST8CTpR0YET8qtVjMrP214oZIxHxAqmU/u+BMcBF\npOKJ87MxnZq1mb1VXyJtWu9Z47phjqSo9QHMafVArX21bBNxMzOzdtSSGaOZmVm7cjCamZnlOBjN\nzMxyHIxmZmY5DkYzM7McB6OZmVmOg9HMzCzHwWhmZpbjYDQzM8txMJqZmeU4GM3MzHIcjGZmZjkO\nRjMzsxwHo5mZWY6D0czMLMfBaGZmluNgNDMzy3EwmpmZ5TgYzczMchyMZmZmOQ5GMzOzHAejmZlZ\njoPRzMwsx8FoZmaW42A0MzPLcTCamZnlOBjNzMxyHIxmZmY5DkYzM7McB6OZmVmOg9HMzCznvwG5\nG4bAQgGyYQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import kf_book.ekf_internal as ekf_internal\n", "ekf_internal.plot_bicycle()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "那么,一下是设计过程,在设计过程中的注释里面会一步一步的解释这个过程,可以不断浏览不断回顾之前学过的思想与方法,从而掌握ukf的方法。\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from math import tan, sin, cos, sqrt, atan2\n", "\n", "def move(x, u, dt, wheelbase):\n", " \"\"\" 生成实际位置和方向的状态量,用于拟合出一组用于分析的测量量,现实中不存在。\"\"\"\n", " hdg = x[2]\n", " vel = u[0]\n", " steering_angle = u[1]\n", " dist = vel * dt\n", "\n", " if abs(steering_angle) > 0.001: # is robot turning?\n", " beta = (dist / wheelbase) * tan(steering_angle)\n", " r = wheelbase / tan(steering_angle) # radius\n", "\n", " sinh, sinhb = sin(hdg), sin(hdg + beta)\n", " cosh, coshb = cos(hdg), cos(hdg + beta)\n", " return x + np.array([-r*sinh + r*sinhb, \n", " r*cosh - r*coshb, beta])\n", " else: # moving in straight line\n", " return x + np.array([dist*cos(hdg), dist*sin(hdg), 0])" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def normalize_angle(x):\n", " \"\"\" 用于处理转动的周期性,使超过2pi的角度自动正规回2pi内\"\"\"\n", " x = x % (2 * np.pi) # force in range [0, 2 pi)\n", " if x > np.pi: # move to [-pi, pi)\n", " x -= 2 * np.pi\n", " return x" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def residual_h(a, b):\n", " y = a - b\n", " # data in format [dist_1, bearing_1, dist_2, bearing_2,...]\n", " for i in range(0, len(y), 2):\n", " y[i + 1] = normalize_angle(y[i + 1])\n", " return y\n", "\n", "def residual_x(a, b):\n", " y = a - b\n", " y[2] = normalize_angle(y[2])\n", " return y" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def fx(x, dt, u):\n", " \"\"\" 分别初步确定每个sigma点下一个状态的预估值,相当于KalmanFilter中的 X_bar = FX + Bu。不过还需要state_mean过程,对每个sigma点进行加权\"\"\"\n", " return move(x, u, dt, wheelbase)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def Hx(x, landmarks):\n", " \"\"\" 分别输入每个sigma点得到的状态预估值。\n", " 输出:分别初步确定每个sigma点的这个新状态,相对于这组landmarks,测量时所需要描述的维度上的值,相当于KalmanFilter中的 HX_bar。\n", " 不过还需要z_mean过程,对每个sigma点进行加权,\n", " \"\"\"\n", " hx = []\n", " for lmark in landmarks:\n", " px, py = lmark\n", " dist = sqrt((px - x[0])**2 + (py - x[1])**2)\n", " angle = atan2(py - x[1], px - x[0])\n", " hx.extend([dist, normalize_angle(angle - x[2])])\n", " return np.array(hx)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def state_mean(sigmas, Wm):\n", " \"\"\" 对每个sigma点的状态预估值进行加权,最终确定下一状态的预估计值,即:KalmanFilter中的X_bar\"\"\"\n", " x = np.zeros(3)\n", "\n", " sum_sin = np.sum(np.dot(np.sin(sigmas[:, 2]), Wm))\n", " #print(sum_sin-np.dot(np.sin(sigmas[:, 2]), Wm))\n", " sum_cos = np.sum(np.dot(np.cos(sigmas[:, 2]), Wm))\n", " x[0] = np.sum(np.dot(sigmas[:, 0], Wm))\n", " x[1] = np.sum(np.dot(sigmas[:, 1], Wm))\n", " x[2] = atan2(sum_sin, sum_cos)\n", " return x\n", "\n", "def z_mean(sigmas, Wm):\n", " \"\"\" 顺Hx()函数的结果,对每个sigma点加权,最终确定这一新状态,相对于这组landmarks,测量时所需要描述的维度上的值,即:KalmanFilter中的HX_bar\"\"\"\n", " z_count = sigmas.shape[1]\n", " x = np.zeros(z_count)\n", "\n", " for z in range(0, z_count, 2):\n", " sum_sin = np.sum(np.dot(np.sin(sigmas[:, z+1]), Wm))\n", " sum_cos = np.sum(np.dot(np.cos(sigmas[:, z+1]), Wm))\n", "\n", " x[z] = np.sum(np.dot(sigmas[:,z], Wm))\n", " x[z+1] = atan2(sum_sin, sum_cos)\n", " return x" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from filterpy.stats import plot_covariance_ellipse\n", "from numpy.random import random, randn\n", "\n", "def run_localization(\n", " cmds, landmarks, sigma_vel, sigma_steer, sigma_range, sigma_bearing, ellipse_step=1, step=10):\n", "\n", "## found the MSSP & UKF modle\n", " plt.figure()\n", " points = MerweScaledSigmaPoints(n=3, alpha=.00001, beta=2, kappa=0, \n", " subtract=residual_x) \n", " # 生成2n+1个sigma点,并生成对应每个sigma点的mean权系数W_m和covariance权系数W_c\n", " # n = dim_x, 0<=alpha<=1, beta = 2, kappa = 3-n\n", " ukf = UKF(dim_x=3, dim_z=2*len(landmarks), fx=fx, hx=Hx,\n", " dt=dt, points=points, x_mean_fn=state_mean, \n", " z_mean_fn=z_mean, residual_x=residual_x, \n", " residual_z=residual_h) # 生成dim_x个状态变量,dim_z个测试变量,作用在sigma点集points上的UKF模型。\n", " # 其中由函数fx表达状态传递方程y=f(x),由函数hx表达状态测量方程z=h(y)。\n", " # x_mean_fn表示状态变量期望μ(x)计算函数,即wm对y加权求和,也即:KalmanFilter中的X_bar;\n", " # z_mean_fn表示测试变量期望μ(z)计算函数,即wm对z加权求和,也即:KalmanFilter中的HX_bar。\n", " # residual_x和residual_z都是用来解决变量做差运算时出现的不合题设要求(如周期)的问题。\n", "\n", " ukf.x = np.array([2, 6, .3]) # 初始化状态量,分别为x、y位置量range和当前方向量bearing\n", " ukf.P = np.diag([.1, .1, .05]) # 初始化状态协方差,分别描述x、y、bearing各自独立方向方差\n", " ukf.R = np.diag([sigma_range**2, sigma_bearing**2]*len(landmarks)) # 描述len(landmarks)个标记点(测量点),各自测量测量的不准确性\n", " ukf.Q = np.eye(3)*0.0001 # 描述状态变量在变化过程中出现的过程不确定性\n", " \n", "##\n", " sim_pos = ukf.x.copy() # 对状态初值的浅层复制\n", " \n", " # plot landmarks\n", " if len(landmarks) > 0:\n", " plt.scatter(landmarks[:, 0], landmarks[:, 1], \n", " marker='s', s=60)\n", " \n", " # plot the points\n", " track, xs, covs, Qs = [], [], [], []\n", " for i, u in enumerate(cmds): \n", " sim_pos = move(sim_pos, u, dt/step, wheelbase) # 真实状态的记录,实际上不用。其目在于生成之后的测量值,而一般这是实际测量的\n", " track.append(sim_pos)\n", "\n", " if i % step == 0: \n", " if i % ellipse_step == 0:\n", " plot_covariance_ellipse(\n", " (ukf.x[0], ukf.x[1]), ukf.P[0:2, 0:2], std=6,\n", " facecolor='k', alpha=0.3)\n", "\n", " [x, y, theta] = sim_pos \n", " z = []\n", " for lmark in landmarks:\n", " dx, dy = lmark[0] - x, lmark[1] - y\n", " d = sqrt(dx**2 + dy**2) + randn()*sigma_range # 生成测量值中的目标距离测量点的距离\n", " bearing = atan2(lmark[1] - y, lmark[0] - x)\n", " a = normalize_angle(\n", " bearing - theta + randn()*sigma_bearing) # 生成测量值中的测量点-目标-原点之间的夹角\n", " z.extend([d, a]) # creat the measurement variate\n", " ## ukf滤波器求解,也可以用ukf.batch_filter()函数\n", " ukf.predict(fx_args=u) # ukf.predict(u) 给一个输入的控制量,预估计一下接下来的状态\n", " ukf.update(z, hx_args=(landmarks,)) # ukf.update(z) 再用一组各测量点贡献的测量值,后验估计这一新状态\n", " \n", " xs.append(ukf.x)\n", " covs.append(ukf.P)\n", " Qs.append(ukf.Q)\n", " \n", " if i % ellipse_step == 0:\n", " plot_covariance_ellipse(\n", " (ukf.x[0], ukf.x[1]), ukf.P[0:2, 0:2], std=6,\n", " facecolor='g', alpha=0.8)\n", " \n", " track = np.array(track)\n", " plt.plot(track[:, 0], track[:,1], color='k', lw=2)\n", " plt.axis('equal')\n", " plt.title(\"UKF Robot localization\") \n", " #Means, P, K = ukf.rts_smoother(np.asarray(xs), covs) # 在这里我们无法使用,是由于fx()函数需要控制量u的输入,而rts_smoother需要fx()函数却不需要控制量u的输入,出现了矛盾。\n", " #ukf_internal.plot_rts_output(np.asarray(xs)[:,0], Means[:,0], range(200))\n", " plt.show()\n", " return ukf" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd0XMX58PHvbF/trnq3imVbcsUGdzCmmW5KaKEnIRAI\nSSCQX0IJJAQIqYQWQwJ5IUBoIfTeDLYxYIN7L7Lqqq5Wu9re5/1j5SKvZAsjGVuezzk6lu7cO3fm\n2n40O3eKkFKiKIqiHPw033YBFEVRlIGhArqiKMoQoQK6oijKEKECuqIoyhChArqiKMoQoQK6oijK\nEKECunJQEEL8QAixeBDzrxNCnDiI+R8nhLDv8vN6IcRxg3AfnxBixEDnqxwcVEA/xAkhpBBi1G7H\nfieEeKb7+90DkUEI8YoQ4jMhRHr3udHuQLL966Y+7lUnhAh2n9MqhHhSCGEd3Br2rM+BQko5Xkq5\n4JvkIYRYIIS4ard8rVLKmm9UOOWgpQK60m9CCCPwCpAJnCyl9HQn/bc7kGz/+ssesjlTSmkFDgeO\nAG4d3FIryqFDBXSlX4QQacCbgA6YK6X0f5P8pJStwPskA/v2e2QIIZ4WQjiEEPVCiNuFELv+GxVC\niHlCiC4hxCYhxJxdEoqFEG8IITqFENVCiB91Hz8V+DVwYfcng9X9qKtRCPGAEKK5++uB7l9m29PP\nFkKsEkJ4hBDbuu+BEOIKIcRGIYRXCFEjhLhmD/fY0cUjhHDv8unG3/2pabgQIksI8Vb383B1f1/S\nfc09wGxgXvd187qP7/jEtafnub0LSwhxb3fetUKI0/b2bJQDmwroSn8YgXeBEHC2lDL4TTPsDkyn\nAdW7HP47kAGMAI4FvgdcsUv6DGAbkAvcAbwihMjuTnsBsAPFwPnAH4QQJ0gp3wP+wM5PEZP6Ubzb\ngJkkf9lMAqYDt3eXezrwNPArkp9UjgHquq9rB84A0rvLfb8QYvLebialzNz+6QZ4EPgUaCL5//Pf\nQDlQBgSBed3X3NZ93s+6r/1ZL1n353luJvk8/wI8LoQQeyuvcgCTUqqvQ/gLkMCo3Y79Dnim+/vj\nSAbyCHBeL9f/rjvNvctXcR/3qgN8gLf7vvOBzO40bXc+43Y5/xpgQff3PwCaAbFL+pfA5UApEAds\nu6T9EXhy9/rs4TnUASd2f78NOH2XtFOAuu7vHwXu7+ezfQ34+S7P0d7b/XY5dmH38bw+8jsccO3y\n8wLgqt7+Pvv5PKt3SUvrvrbw2/43qb72/Uu10JU4oN/tmB6I7vJzB3AR8JQQ4pRe8nhRJluZ27+a\n93C/70gpbSQD3BiSrUO6/9QD9bucWw8M2+XnJtkdfXZJL+7+6pRSevdw7ddR3Es5iru/LyUZ8FMI\nIU4TQizp7vZxA6ezs357JIQ4gmTr+xwppaP7WJoQ4tHu7hIPsAjIFEJo+5Flf55n6/ZvpJSB7m8H\n/SW1MnhUQFcagOG7HaugZyBASvkK8CPgJSHE8d/0plLKhcCTwL3dhzpI/hIp3+W0MpJdD9sN261L\noIxkq70ZyBZC2Pq49usuKdrcSzm2/5JqBEbufkF3H/vLJOtTIKXMBN4B9tqFIYTIJ9ma/6mUcuUu\nSf8HjAZmSCnTSXbvsEuee6pXf56nMsSogK78F7hdCFEihNB0v6g7E3hp9xOllM8DPwNeF0LMGoB7\nPwCcJISYJKWMAy8C9wghbEKIcuAXwK7DDfOB64UQeiHEBcBY4B0pZSPwOfBHIYRJCDERuHKXa9uA\n4bu9YN2T50k+kzwhRC7w213yehy4Qggxp/t5DRNCjAEMJN81OIBY9wvGk/d2IyGEjuSzfkZK+eJu\nyTaS/ebu7ncFd+yW3kayfzxFP5+nMsSogK7cRTIYLgZcJF+OXSqlXNfbyVLKp0i2HN/ufkG4z7q7\nFp4mGTABrgP8QE13eZ4DntjlkqVAJcnW5z3A+VJKZ3faxSQ/aTQDrwJ3SCk/6k77X/efTiHEin4U\n7ffAMmANsBZY0X0MKeWXdL/wBLqAhUB5d3fP9SSDqAu4BHijH/cqITla5QbRcyx/GclfeObu+i4B\n3tvt2geB87tHqTzUS957e57KECN6dkkqiqIoByvVQlcURRkiVEBXFEUZIlRAVxRFGSJUQFcURRki\ndPvzZrm5uXL48OH785aKoigHveXLl3dIKfP2dt5+DejDhw9n2bJl+/OWiqIoBz0hRP3ez1JdLoqi\nKEOGCuiKoihDhAroiqIoQ4QK6IqiKEOECuiKoihDhAroiqIoQ4QK6IqiKEPEfh2HriiKcrBrd/nZ\nWO/Almbk8FGFaDQHzjasKqAriqL004otLfzsoXfo9Heh1+gZVZzPRceP54Ljxh8QgV0FdEVRlH56\n7K3l1HY2EDY0Ew/padnSwuamFhasruOeK+eQnW7+VsunArqiKEo/VTd14o8EqJi5Hp0hgrc9l6aN\no3h3RYDGdg8PXXcaI4qz+p2f0+kkIyMDnW5gQrF6KaooitJP2/coF0IiNJBe2EHFUSvw6bexon4L\nP/zL66yrbe9XXtXV1UyfPp2f/vSnDNTOcSqgK4qi9JPFpEen0RILG3cc0xkjlE1bSyKjho1t27jh\n4fewOzx7zGfZsmUcddRR1NTUsHz5cnw+34CUTwV0RVGUfppcWUSaPg2/M7PHcY02wbBJG5DpDWxs\nqeP6v7+Lxx/uNY8PPviA4447DofDwcknn8yCBQuw2WwDUj4V0BVFUfrpyPElWA1WvG05KWlCAyWT\nNhHUN7Kyrob/+8f7RKLxHuc888wzzJ07F7/fz6WXXsqbb76J1WodsPKpgK4oitJPRx9WRo41g4gn\ni0jAmJKu0cUpnbyezpidReu2ctfTC5FSIqXk3nvv5fLLLycWi/GrX/2Kp59+GoPBMKDlU6NcFEVR\n+sls1HPcpOE0dDbjaiimYExtyjl6U4TSyeto+ErHa58bmFpVxKJX/8X9998PwH333ceNN944KOVT\nAV1RFOVruPzkSby9dDM19mHkVNjRGaMp55jS/RSM3UrTWj1XX3UFji1L0ev1PPXUU1x88cWDVjYV\n0BVFUb6GMWW5zDliJM7PnDhrSykYU9MjfXOjE4CROXEalq4k0tSK3mjirTfe4OSTTxrUsqmArijK\nIe29pVt5+oM1RGJxZo0v5aI5EyjK2fOok2vOnMr8ldvY1ljc3UqP9EhPeGNseWozkaYgwmSm5JSf\n4jOVDGY1ABXQFUU5hC3ZYOfa+9/A6awnIRMsWpnJI69/yVVnTOGX3z0Krbb3cSNVpTmcNGUUzsVO\nOmpKKRy7bUfL3Nfgg1c6wZdAZOkov3wq7hYtD7y8hGMmlZObkTZo9VGjXBRFOWQ9/f5q3F3NSG0N\nluyVBEJraWhczUMvfcLFd79Mly/U57VXnzGFfEseHnsJIU9y6GGsJggvOMGXQFNiJO2KIrLHBdBn\ntWB3tfHvd1cOan32GtCFEE8IIdqFEOt2O36dEGKTEGK9EOIvg1dERVGUweH0BIhGQ1iyOrDltVFY\nuZ7sknV0ODbyybKVnPOb/9LQ1tXrtZUlOVw6ZxIFaYU0r6skpx5CL7RBRKIbm8akW8YxZnQ+QkBe\nVR0dAScvf7qBdpd/0OrTnxb6k8Cpux4QQhwPnA1MklKOB+4d+KIpiqIMrlA4hpQJhGbnBCBLVicF\no9bi81ezatM6fvCn1/psqV979lTGDBtG+KtN1P+nDhKgPzId47l5aPQ7w6vJ5ictv5U2bwdPvrdq\n0Oqz14AupVwEdO52+FrgT1LKcPc5/VuNRlEU5QBSnGtDpzMSC/Xs1zaYgxRWriUYqmdt9WauuvcN\nYrF4yvVaIYmve5nAhsWAoPCcSib+oJIxZbkp5+aNrKczOLit9H3tQ68CZgshlgohFgohpvV1ohDi\naiHEMiHEMofDsY+3UxRFGXiTRhViMloJ+TJS0rT6KHkVG/F01bF49UZ+8cgHPVZFdLlcnHLKKbz7\n5isYjGaK51xByHIC8WjvY02MtgBp+a20+zp46v3BaaXva0DXAdnATOBXwIti+7qSu5FSPialnCql\nnJqXl7ePt1MURRl4J00ZgTktk3DARiKRGsIM5iC5ZZvodG7j1UWrePClpQDU1dUxa9YsFixYQGFh\nIYsWLeTYY0/HKkuwrxxHIt777kU5FY24gi7e+7KaWDwx4PXZ14BuB16RSV8CCSD1M4aiKMoBbOSw\nbCoK89DrbfhdqQtuAZgzPGQWVeNo38rDry3h2VfeY+bMmWzcuJHx48ezdOlSZkyfxv0/PYXxxcPR\n+ktoXV9Fb0ucG21+NOYuWrvcLN1gH/D67GtAfw04HkAIUQUYgI6BKpSiKMr+ct4xY7GlF+B1lNDX\nNhO23HaM1haa1n3I9y48m7a2Nk444QQWL15MWVkZAAXZVh742amMzC0n4iinY1t5Sj5CgK3QgTfs\n5bN1jQNel/4MW3we+AIYLYSwCyGuBJ4ARnQPZXwB+L4cqC03FEVR9qMrTjucorwiEvFMAu6e28fV\ntriobXEhpUTb8RXBNe+QiEWYcsxpvPvuu2Rm9lwXfUxZLn+++iTKMkvpqh1JZ31xyv0sOW78UT9L\nNzYNeF32OlNUStnXSjKXDXBZFEVRvpF4PMEDLy/h83WN6LQaTpoygnNmjyVnD7MzLWYDFxw7jgfa\n7HjanFgyXT3SZVzSPr8d35bkrkLa4VPxlp5Co8PHyGHZKfkdfVgZt11yLHc9k6BhqyAWMpJXVcv2\nt4wmm59wLEJr58DsUrQrNfVfUZQh444nF/Dvdz6jy92KRqNh4Yp1PPjKEn59yTFcfOJhfV53zZlT\nefajNXS5mwl50mnx1wOQCMVhqQ9fZxy0UHhSIT4KcDjt/P6ZRfz75u/0mt95x47DoNdy59MLqG/S\n0ei1UDxxMzpDFIQEJLLPDp59pwK6oihDQjgS46UFG+hw1JKeXw1IPO48PJ5cbnnUx1ebm/nLNSei\n02lTrs3JSOPc2eN4rLOVziYPiYwGhCcKX/ggmACzBjHTimWEBUO4htbNOSxaXcum+g7GlPc+HuTM\no0ZTkGXl1/9vPlva0ti2KIP0ojZ0higGrYEMi2nAn4Fay0VRlCHB7vAQjARABMgoaCajoIWi0Wuw\nZG+mtW0DL3y0lMv+8GrKtnDb3XzxLCpKhiNEHra2QvjUmwzmWVrKLyxlxNjksGu9MYLJ5sDtbuGR\n17/aY5mmjx3Gc7efx3emT2Zkxmj0nROJNU2kyFbEubPHDPgzUC10RVGGhBanj3gsilYb23FMABmF\nTRjMPjoaoixcKbnxYTPzfn46u0+dsaYZuemio7jyutdwr9kASCgxICZb0KX1DJXphXbaqwt478st\n2B2zKMlL77Nc+VkWHrr+NGpbXLz+2WbcvhATKvI5d/bYgaw+oAK6oihDRHlhBjq9kXjMgCQZzLcz\nZ3SRV74JR72O1xfrKS/I4KaLj+5xfSQS4b1nHsC95m0AdBVjKDklglaTOknIaA5gMDvpdLfyyGtf\n8Ycfzdlr+SqKsrjh/JnfpIp7pQK6oiiMv+M9/OHeuyIALEYt6+88tc/0A0FJXjqZVistmIgGzRjM\nwR7ppnQPmcVbcDTreOQ1A2X5GVw0J/mi1Ol0csEFF/DJJ59gMpkomnUZ7VJPZ2MNeeXVvd4vPc+O\n057P/BW1SClTWvzfBtWHrijKHoN5f9IPBEIIJo4owGi0EvJm9nqOLacDa04d7e1buOvpBWxr6mTN\nmjVMmzaNTz75hMLCQhYuXMgT9/6aosIqIr5SulqH9ZqXMd2DTPhp7XCyfHPzYFat31RAVxRlyDj6\nsFLS0rLwufpeNyqzuAGtsZmW9jouvv5ujjzySGpra5kyZQpffvkl06dP57jDK7jjByeQX1CJxzES\nX2fqSBaNAKO1i0Cwi8VrB37W575QAV1RlCHjohMOIy+ngEQkk5C3576g22d9CiC7dBue9e/z1SsP\nEAgEuOyyy/j0008pLS3dcf5lJ03kR2fMJC9vJK6mSoJdqSsymixewmEfq2vaBrtq/aICuqIoQ0a6\nxcjcmVVY0vPwOFKn3QMkIgk6Pmog0bgeEBTPOI8H5v0Ts9mccu5tl83mjFmTyMmtpKNhAt6Ogh7p\nWn2URCKOPxhJufbboF6KKopywJJS8sZnm1lZ3UpGmpHTZlT2OZFnux/Nncz/FqzF3thKNFyLvTPZ\nek5ICb44tR/VgzeBxqhBN/ZIInmjeeClpdx95QkpeQkhmHf9aZiNOl5aoMfRoiPks5FdWotWGyce\n1aLRaLGYDYNS/69LBXRFUQ5YNz/6Ec/PX4bf34lWa+DBl7OYXFXKX398Uq/rqEBySdxZE4bzpsuO\nu3k4mLq7Q9qi8JUfYhJsWkrOGkZMBHA2tPH6Z5u55ZKjew3MOp2W+35yCiMKs3jwFTPtDitN6/Mw\nWrqIhszYbDYqh2WlXPdtUF0uiqIckNZsa+PFBatpb9tEQrOaUGwVLS0r+WTZl5x123O8tnhTn9f+\n+tLZ5GaXEvIVUphWRmYLyWn8MYmlwkLFhWXoM/SY0z1odG4cnW3858M1feYnhOC682bw3O3nc/zU\nmQwbdjhmw0QyM8cwqmwE15w5dRCewNenWuiKomAxavc6Dn1/W7rRTiDgxpDmJLe8BoBY1ICzfhS1\n9T5unBeitdPHj89KDaZjynO58ISJ/OvlNlo/dJJoT26LLMaaKTiuoMeYcVtuC56OPF74eB3XnDll\nj+PJp48t4fU/XMyWRicfr6jBbNRzzuyxpFuMA1z7faMCuqIoB+SkoW3NLqLRIHpTYMcxnT5C/qgN\nuJp8tDTH+MvzUJxj46xZo1Ouv3BmEfff/AwJpx20WgpPzsNSYUk5z5LdjrvVzbamFuavqOXEKSP2\nWraq0hyqSnvf4ejbpLpcFEU5IFlNejRCi0z0/HQggOxhDZiz6mlt3crNj37AV5t6bhbx1ltvMXvW\nkQScdnTWXHSTzkRmpO4gBKDRSNLSnfj9nXyysnawqrNfqICuKMoBqawgE73eRCycOpwQILukDp3J\njr1lMz++7y0a27uIx+P89re/5cwzz6Srq4uzzz6bm+99ioLhU3E1VRLyWXvNy2RzEw77WLG1ZTCr\nNOhUl4uiKAekKaOLMJqsuFyWlMW26P45d/gW2rYaqLOnccP9r+Jf+QLvv/8+Go2Ge+65h5tuugkh\nBB3+BK8tiuCojVNYuQa9KdQjL5PNg7MxSEObh3g8gVZ7cLZ1VUBXFOWANKY0l7zMTNrbrIS6MjBn\ndKWco9FIcodvpnl5iNcX/ot4wE1ubi4vvPACc+bsXAHxwZ+dSmunj8WrIrTXxMir2Nhj8S6tLgZE\nCYRDdHQFKMjuvSV/oDs4fw0pijLkabUaTp9RicWai8dZmJK+fSp/sKaD2Jr3iQfcWPKGs2jx5z2C\nOYDRoOPxX53F4VVjsVoraas+DL9r5wSlRFwgpRa9Vk/uHvYfPdCpgK4oyn7h8YeZ9+qX3PTPD/nn\nG8uobXHt9ZofnDqJdFseYV8Osai+R5qMSxIr/Tg+cUBCoi2uwDLtYl5a2vvKhzkZabxy94XMmXYE\nuXljcdnH0VY9Dm9HPu7WUvR6M/lZloO2uwX60eUihHgCOANol1JO2C3t/4B7gTwpZcfgFFFRlINd\npyfAWbe9QHVDHeGIH4M+jT8/l8npM8dw9w+PJzu991bxiOJspo0ppcNZj8dRRHZxA7UtLmQgjlzq\nA3ccNCAOt5A/NodOezuvfLqRX104C5MxNbylW4w8e9u53P/SEh5940vcHgf+DjcIyMkt4bKTJg72\noxhU/elDfxKYBzy960EhRClwMtAw8MVSFGUo+dNzn7G5rhqfrxpzhpNA0ILblcWL852s2NrKwz8/\nnclVRb1e+/1TJvHFuhra2jpIz2tFtkaQy/0QkZCmQcywIjJ1mDNciGY3bc52Xl60gUv7CM5arYZf\nXngU5x8zlpc/3cii1fVoNILTZ1Ry1dzJg/kYBt1eA7qUcpEQYngvSfcDNwGvD3CZFEUZYr7a1Izf\n10HmsFosGcmulnAgDWeDn/VbfFz6+wD/vuUcZo4rSbn2tBmVzBhfwYddrbTMDyNrfMmEAj1iqoUR\n5Tsn+FizHXi72nnxk/V9BvTthhdl8X/fPYr/++5RA1fRb9k+dRYJIc4GmqSUq/tx7tVCiGVCiGUO\nh2NfbqcoykEsFI7R0J6c9Wm2unccN6YFKKxci9TWYG/ezE8eeKvXfnUhBD8+pYrwqjeI1NSAgOwZ\n2YgjrQhDzxBmy20lEnKzelsT6+vaB71uB5qvHdCFEGnAr4Hf9ud8KeVjUsqpUsqpeXl97yKiKMrQ\nZDRo0Wq1gEgZTK7Rxsmr2EyCJmobtvLDP79Bl6/nGPH333+fi846iZCjBo3Jin7i8dgmFjCiOJuK\noqyU/IxWF/6Ai09W1A1uxQ5A+9JCHwlUAKuFEHVACbBCCJE6rkhRlEOeEIJMqwmtRkcsbEpJ12gk\n+SM3Ego3sq56Mz998B2klMRiMW6//XZOO+00Ojo6mDPnRGZf8SesRVNo3zaWRLz3BcOMaT4iYf8B\ns4vQ/vS1A7qUcq2UMl9KOVxKORywA5OllK0DXjpFUYaEsvwMdAYz4UDvE3a0uuRkH4+ngYWrtvLY\nyws48cQTueeeexBCcPfdd/PBB+/zzJ2XUzm8Cp1mGI66ShIyNS+9OUA0FqK5wzPItTrw7DWgCyGe\nB74ARgsh7EKIKwe/WIqiDCVTqoowmzMIuFNXKNw+QchgDpKe14hj66f87PKzWbhwIYWFhXz00Ufc\nfvvtaDQaSvLSeez/zqS0pIpEpARnXSWJRM9+HK0+QiIRw71b182hYK8BXUp5sZSySEqpl1KWSCkf\n3y19uBqDriiHFikl1XYnK7a0EIv1vY76dhedMAGbNZewP4t4rPfBdTIhidWvJbz6bWIhLyPHT2HV\nqlUcf/zxPc6bNKqQ+649lWHFY4iHR9C2dQLR0M71yCNBMzqdkWG56d+skgchtZaLoihfy6b6Dq57\n6F222ttIJOJkWG1MH1vCHd8/lrKCjF6vqSjKYtKoYjo6avE588koaN4xoiUhJQQT1LzUAB0xALTl\nR2Cc8n1s6b1v7XbK9FE8cdO5XPf3d2hostC61YLJ6sRk8eDtLMKSlsGYsj3vPToUHbxzXBVF2e/C\nkRg/vv8tlq1fRUvLCtodK9iybRmvLviU027+D89/tLbPa78zawwWay4+Z0HPvu+WCHziSQZzo6Do\nrCIMI0fg8rh4dQ/bzB09sYy3/3gJpx41g8LCiejERILuyZhNIykpLOO6c6YPYM0PDqqFrihKv32y\nspbqxiaCITvF41ah1caJBM24miqore/ilscC1La6+PVlx6Rce+HxE/jHG8vwdLXgcxRSnpfA+YUT\nz1o/AOZSM/lz8tGl6Yi2OvB5Onjts017nCBUnJvOM7edy1ebmnhnaTU1zS6Kcqz89DvTyctK3Z1o\nqFMBXVGUfvt4ZR2BgBtzeidabbLv3GAOkj9qA57WLtpaozz6BgzLTef7px7e41qTUcdPzp7Gr9va\naK9pxP3xV0Q7I8nh6ePNFM0u2rGfpzWnna52Nyu3NNHY3kVpfu9dOZAcFjl9bAnTx6bOMj3UqC4X\nRVH6zRMIE4tH0BuDPY4LIKOwCVteDe3tW7n7Pwv4eEXqdm6XzJlAdqCW6KoPiHZG0GfoKTmvhBHH\nFPfYnFmrj2Iwu/H5XHzw1bbBrtaQoQK6oij9ZtBp0QhNyj6f22UUNmGwNtDauoUb5r1LQ9vOTSlc\nLhcXXXQha997HBIxNAWV2I6dijHf2GtehjQf4aifzY3OQanLUKQCuqIo/VaUY0WrMxCN9B6EAXLK\ntiE1zTS11nHHk58AsHjxYiZNmsTLL7+MzWbjR7/8PSWzr8LvGo3fld1rPnpjiHgsQmunb1DqMhSp\ngK4oSr+NKcvFaLQQCaSO8d4+QUgjILesGp+vhfnLNvGDH9/IscceS2NjI9OnT2fVqlU89tfb+P5p\n08jNG0mnvYqQ15aSn0YbI5GI4wtG9kfVhgT1UlRRlH47ZepIstJzcDpsRIImDObeZ2PqTSFMxm00\nfzSfp9wtCCG45ZZbuOuuu9DrkzsP3fH9Y2l0eHhvSRxHnYaMgm3Y8tt2rN8VDljR682UF2bup9od\n/FRAVxSl39LMBmaOK8XeshW/Kx+DuaHnBCGgprkTGiKwthYZTaAx2fjTfQ/zq2sv75GXVqvhHzfM\n5ReP6HjzcwMdHUb87jysmQ40+ih+Vz45OZlMH1283+t5sFIBXVGUr2XuzEreXbIalyuXjKLdNiwL\nJ5CrA9AcBUBflI957Hk0xXtfOttk1PHwDaczpaqI+1/+gvaOFgJuN4lEnPT0bA4bWcF5x44b7CoN\nGSqgK8ohzh+M8MLH67A7PEyoyOeIykJGFPf+ohLg9BmVlBYU4XbZ8TsLqOjeOa5mVVtya7iwROgF\nebPz0JcU0lEb5tM19cTjiV43YBZCcOXcyZw2YxQvLdzI4rUNBMNRxg3P4+aLZ2HQ9z6iRkmlArqi\nHMJWV7fy0wffYZu9kUjYj8FowWS0ctwRI7nriuN7ndBjMur40dwj+G1bC51tbszWdlxftiPXJUej\nmIpM5M/JR5+uR+JHaL20dTr5ZFUdJ04Z0WdZinPTuf68GVx/3oxBq+9Qp0a5KMoh7A/Pfsr6retx\nu9cT06zB41tJc/Mq3lj4OWfc+hxvf7Gl1+t+cOrhjBtejjYIjS8241nnAQ1kz8ym+Oxi9OnJF58C\nMFrdBANuPl/XuB9rdmhSLXRFOUTZHR6Wb7bj9zsoHr0OnTE5PDAaMtHR4KW6zst1DwWRci5nHDV6\nt6slFZG1zP/qJZAJtOkmik7JxZiXOj7dYAoS7opgdxx6G07sbyqgK8oh6q0vtuD1OTGY3TuCOSSH\nHBZUrqOzIUBra4Jf/VNPYbaVqWOGAbBt2za+973v8fnnnwNgGXkkkYJySNsIpE4C0mijJBIxNZ58\nP1BdLopyiAqGosQTcXSGcEqaRkBOeQ1akx17y2auue8tapqcPPLII0ycOJHPP/+coqIi3nvvPS75\n8a1k51aBNPLyAAAgAElEQVTRXjuWSNCckpdMaFM2h1YGh2qhK8ohymDQIoRAxntv1wkgr3wLbdUG\nqrdEmHn033HUrQfgwgsv5OGHHyYnJ4djj4vR1unj87Ux2qq1ZBbXYM1p3xHDA95MTMZ0jqgs2j8V\nO4SpgK4oh6iibBs6nZFAqO91WYQmgSmyhM7POyAew5aRxeP/epQLLrhgxzkmo47/3HYuP3vQxMfL\n03C2mvF2FGKyuknEdUT8BeQW5zB35qj9Ua1DmgroinKImlxZiMlkw9lhJRHXotHu3Bu0tsWFDCYw\nrQ8TbEwulavNq6Bq7nWcdfY5KXmlW4w8det3ePSNEv7+2lI63U5CIQ8CDYVFeVxz5nTGDc/fb3U7\nVKmAriiHqOFFWYwuzcfRno7fnYUtJ7nXu5QSWR9Grg0QjEo0Rg05s3PxRKfQ5Pbw+Dsr+cl3pqXk\nJ4Tgx2dP44LjxvPh8hq+3NiEUa9lzpQRexx/rgwcFdAVZQjZ1tTJO0u20ukNMqEin6MmlFKUk7qS\n4XZzJlewfOMGAl052HI6qKl1Ilf6oTU5dZ9CPfIIC+kVNkRHM96OPN78YnOvAX27nIw0LjphAhed\nMGGgq6fsxV4DuhDiCeAMoF1KOaH72F+BM4EIsA24QkrpHsyCKoqyZ898uIa7nvoYV5eDRDyKwWjB\nmpbB+cdO4NZLZ5NuSe0rP3vWGB5+7Qu6mrJwbwgiP+uCqASdgIlmRJlxx05ClqwOXE0eNjW009DW\nRVlB39vCKd+O/gxbfBI4dbdjHwITpJQTgS3ArQNcLkVRvoZN9R3c9dTHNNrXE4qsJaFbice7iobG\n1Tzx1gLm3vosK7a0pFxXVZrDlOFZsGUpzgXNEJWYy8yIEzPQlJsYUZxNRVEWABptHKPFg9/fyWuL\nN+3vKir9sNeALqVcBHTuduwDKWWs+8clgNqdVVG+Rf/9ZB1Odws6cwsFozaQU1pL0ei15A1fhce7\nmdUbV/G9P77C6urWHddIKXniiSf44JFfEG2vBq2erKNKKJpbhDD3HhrSMpwEA24Wrq7bTzVTvo6B\nmFj0Q+DdvhKFEFcLIZYJIZY5HI4BuJ2iKLv7dE0DQb8ba1Zbj+Mmm5ei0WuQmjoa7Ju5+m9v0tje\nRU1NDSeddBJXXnklPq+HignTyTruGkKmGciEjoqirB0t810ZLF6isRBtLv/+qpryNXyjgC6EuA2I\nAc/2dY6U8jEp5VQp5dS8vN7XRFYU5ZvxhSIk4lH0pmBKmkYbJ69iM7GEnc21G5jz3Ws57LDDmD9/\nPjk5OTz77LMs++wTKisnotMU46irJCF7v4/OECERj9HlS51dqnz79nmUixDiByRfls6RUvbx168o\nyv6g1QjQaJB9tNE02jhZWato+WAhXd5kD+oll1zCAw88wPaG1j9/cSaX/SFMXX0EV2OEnLKa1Ixk\n8gWpRP2XPxDtUwtdCHEqcBNwlpQyMLBFUhTl68qymdBqdMTCqSNZZFzS+VUnza/WIr2dYEgjf/YV\nXHfbn9n1U/MRlYXcd+2pFBWOJugto6N+JIlEz0VYAu5sDAYLJXmpm0Qr3769BnQhxPPAF8BoIYRd\nCHElMA+wAR8KIVYJIf45yOVUFGUPRhRlYzBaCPl6DiUMtYWoeb4e11cuSED6+HQyjp9OyJTOH5/5\nlN0/XJ8yfRS//f4JlAwbTzxUScumw/F15hILG4hH9Xg6irBacjgrZTld5UCw1y4XKeXFvRx+fBDK\noihKtyUb7Ly0cAPtLj/DCzOZOrqYuTMre93CDeC0GaN4ZeEKnJ3ZQB2JaILOLzvpWt2VPMGiofjE\nQszDzMRjzTRvLGLl1gY+XFbDydNG9sjre6dMYmRxFr/854fU2uvxtufiigRJyARWax45OflcePz4\nQX4Cyr5QM0UV5QDz8Gtf8tcXFtHV1UY8FkKnN5NmzmDs8DJ+d8VxzJpQmnLNCUdUkJ2ZjcNhpWtz\nFPfSZmK+7pHFlUYYY6ZVE4KWEBVFWaRltePxtPHKpxtTAjrArMPKmH/v5Tzw8hK+WG+nttVNMBxh\nclUx91w5h5yMtMF+DMo+UAFdUQ4ga7a18bf/Lqa1eSPG9CbMNi+RUBqOjhw8nja+94cOrj17Br+8\n8Kge1xkNOsYXGdn01iI62huSBzO1iMMtyMzUTZat2e20b3PzxfpGItF4rxsxp5kN/PqyYwCIxxO4\nfSEVyA9wKqArygHk+flr6epqxWBtIa+8esfxrOIGOu3Dsdt9PPhSCJ1Gww0XzAQgFosxb9483rj/\ndqIBP2h1ZE3LIetwC0IjqG1xAfQYV25MC6DR+nC6O1mwqq7XVvqutFqNCuYHARXQFeUA8unaBgJB\nF9nDWnsc12jj5JZvw+vw0d4KD76spzjXysj0MFdffTUrVqwAoHTcDLryJxIzNSE0m/d4L5PNRSDg\nZv6K1H505eCkArqiHED8oQiJWAy9KdRrui2vjVjUSGtTgmuufRPP1sUkEgnKysqYN28eVZOO5Kzb\nnqWuPobf1Y4ly9XrjE8Ag9lPIBikxZm6D6hycFIBXVEOIAIBou8NOKWU6IMbCC/7nFAoDELDDTfc\nyN1334XVagXg+6ccwbxXPDjtEXSGtRgtvU/T1+iiJBJxuny9//JQDj4qoCvKAcSaZkCj0RKPGdAZ\nIz3Sop4oHYs7CNQl5/IJWzbZUy6geMZ5O4I5wM0XH82mRicffBGlvTZOYeVq9LvlBbB9sqfQqB2c\nh4qBWJxLUZQ96PKFCAR7Cai9GFmchcFgIeTdOUEoEUvQuayThucaCNQF0Bg05M7OJf/0EfhkjP98\nsAaXd+caLlqthsd+cSbTxo/BaiujfdsEwgFLyr1C/kwMhjTK8tW65kOFaqEryiBpcXq5498LWLi6\nllhCUpBpZURxFteePa3XseQAxx9ewbtfrMTblU1GYRP+ej8dn3YQ83SPKS8xUDqnCJ1Fh8SFt8NJ\na0cT/3xjGbdeOntHPiajjidv+Q7n3xFlfY2e9m0mLNl2Mgqa0epiRIImAq488vJzOWuWmvU5VKiA\nriiDwOMPc8atz1HXVIPP00ZCxmlqMrK+2sriNTWcNWs8f/rRHNLMhh7XzZ1ZyV1PZeG0Q/PbrQTr\nu/u/bVqYaIY8PY0eL3iSwxBteS142wr4eEVdj4AOya3gXr/nIm569EPeXZqOqzODpg1lCE0IGTeS\nkTmMKaOHc/zhw/fTU1EGmwroijII/vX2cpramggEaiio3ITOECIaTsPnzKel2cWLH3nYUOfgxd+d\nT3b6zvHdFqMGnf0zYivfJJaII/QCRptgpGn7Qoc9pGU46bT7qG529LotXLrFyD9/cQYfLR/H319e\nyuZGJ8FIAJ3GwOxJw/nLNSf1uZyAcvBRAV1RBlgoHOOFj9fj8bSSWdSIwZzs3zam+TGm1WLNdtBR\nH2TVlihX/vUNXvjN+RgNOt5++22uv/56amqSy9Zq8soZdooWQ3oykvc2QUijlRjTPAT8nbz1xZY+\nN28+ccoITpwygkg0Tk2zi+JcW697jCoHNxXQFWWAbW7soM3ZSSLRRVpWR0q60eIjf8QGWrfqWLLW\nyJV3PYVnzeu8+eabAIwfPx7TuLPZ4vLg82wgO71+j/czp7sIuj18uamJn9B7QN/OoNcypjx33yun\nHNBUQFeUAeb0BIknImj0UfoaEKg3hcgpWU/bZ+t59s0NkIhhs9m48847+dnPfsana+1c8eeXaG3p\nwpLVgdHi73OCkNYYIp6I4vKk7lakHFpUQFeUAabZS5e0lBJftQ/n53Uk/HEAMkdO49M3n2fC2OQU\n/BMmV3DajLG8tsCLozZCQeWa3seSA1pNDJlIEIjEek1XDh3qbYiiDLAMiwmtxkAipk9JC7WHqPlv\nA+0fthP3xzHkGjBOPhrdmBN55tPaHufed+3JTBk3Gou1nPaacSTiqSsiAsTjOoQQGHTqv/OhTrXQ\nFWUvwpEYT72/io+W12JLM1JRlMnhIwuYe2QVopdp+odV5JNpS6etzUIkaMJgDhHzx+hc2ol3kzd5\nklGQd2QutjE2IkEv7TWtvL54I7+44Ejys5KTgNLMBv5983c4744wazeHad0qyRu+OWWdl0BXNkZT\nOuOGq03YD3UqoCvKHtQ0d/LDv7zB5vp6/L4OQKDXmzCnZTJpVDl/uGoOk0YV9rhGp9MycWQBdY1W\ngu4M/BtbcC13IaMSBDDSCKPNOA1xnG1uKooEOpOTTncb//14HdedN2NHXvlZFh7/1dlcdk+cbXYj\nrVvN2HLtWLLb0RsjBLsyCHnyKSzK5YJjx+3fh6MccNRnNEXZg3ue+ZS1Wzbidq/HlLkKU+YK4tpV\nONrXsnjlV3z3zv/y9PurU66bMXYYwt2M871VdC7pTAbzIj3ixAyYkAb6ni17S4aTYLCLLzbYU/Kq\nKs3h7T9ewqkzp1NQcBhh7wRaN8+gfvVMOhomkZ0zipOnjWXamGGD9hyUg4NqoStKH5ZssDN/+VZ8\nvhYKKzf06OpIxOvpaBhJQ6Of3z0ZJSfdzNwjqwBYs2YNz/3tl3iWLQZAn2kkd3Y2aaXJCUS9jSc3\nZ3TS1epjTU1brzsI5WVZeOa2c3nz8y08+d4qtto78YeCWExmzppVxW8uP7bX7h/l0KIC+kFi/B3v\n4Q/H+0y3GLWsv/PU/Viioe+ZD9fQ5WnBnN6a0m+t0cbJr9iC0x6htU3LTY99gIkQ//33PB5//HES\niQTGNBua4dPQlepJK63u4y5JemMEjc6P2+tmyYZGjpk0POUcIQRnzRrNWbNGE48ncHqC5GWmqUCu\n7LDXLhchxBNCiHYhxLpdjmULIT4UQmzt/rP3AbLKgNlTMO9PuvL1NbZ3EY0EMKe7+zwna1gdGk0d\n1Qv/w8nHTOdf//oXQgiuu+46Pv7sK4omnUEkUETYv3N524qirF7HlOtMfqKRIHWtXXstm1arIT/L\nooK50kN/+tCfBHZv+t0CzJdSVgLzu39WlCHFF0xuAKEz9D7+WyYkvo0eQl98RKT2SxKxMKOPmMX6\n9et56KGHOOrw0Zwzezy2jGI66quIx3ofdridRkgSMkEgHB2M6iiHgL12uUgpFwkhhu92+GzguO7v\nnwIWADcPYLkU5VsXTySQ23eB2IWUktpV7ch1QfAmPxnpc9KQw45FTDiN9NziHefedtlsvtrUzOpN\nQdq3RSmoXIdGk5onQEJq0Aot1t1WYFSU/trXUS4FUsqW7u9bgYK+ThRCXC2EWCaEWOZwOPbxdoqy\n/5XkpaPXmQn7d24OEXaEaXmjBfmFD7xxdDYdBScVUPrdQowFejpcrTz+9ood52fZzDxx01lUlFUh\nKKajropEIvVeibggErChN5opybPtj+opQ9A3HrYopZTQSzNmZ/pjUsqpUsqpeXlq4oPy7YnHE3y4\nbBt3P72QJ95ZweI1DXj84T7PnzSyEIPRQtifTtQbpe2jNuz/sxNsCiaHHU4wEz/ehsMaRQiBLa8Z\nn7edN77YQvK/RdLIYdnM+/lchhVVkYiU07p1IpGguce9vB2F6HVZVJYUc2wvL0QVpT/2dZRLmxCi\nSErZIoQoAtoHslCKMtCWbLBz2//7mM0NdoLBLnRaPXq9mayMLK48bTLXnTs9ZV3w6WOLMQgt7nW1\neJsbkXGZbAKNMEGVEQyaHi8lzRkuOhv9tLu6qGt193jxOWtCKf+48SxuevRD6ppqadtqQ292ozMG\nkQkNIU8BefnD+OHpR6gXnco+29eA/gbwfeBP3X++PmAlUpQB1uUL8fO/v8fG6jWEI+2YbC4iUR0x\nbxpOZyZ/ed7JJ6tqefiGuZTkpQMQCARY9OZztH1wH4lIchVD6ygr2TOy0Wfoe1+bXCRHqoRCXpZu\nsKeMZDlhcgVv//ES7nxqIe8u3UQw5CEaDYGUFBXn853ZE7jo+PH76akoQ9FeA7oQ4nmSL0BzhRB2\n4A6SgfxFIcSVQD3w3cEspJIcZ763cehK73731AJq7LXEEs0MG7sWjTb5HCXg78zF2eRn8UovP/hj\nlP/dcS7/ff4Z7rrrLlpakq+JdNllGKpKKJjcttd7Gc0+wmEfK7a2ctGcw1LSC7KtPHLjXKrtM1lZ\n3cpWuxNfMMq5s8cwVc30VL6h/oxyubiPpDkDXBZlD9SkoX3T0NbFa4s34Omykz9i645gDsllVazZ\nHZgsHlqro3y1cCvlI2/G62wFYMqUKfz0xlu5+616mppWEw74MKYl9/jsa21yvTlIIBCm2endY7lG\nleQwqiRnYCqpKN3UTFFlSFuyoZFQ0IPO2IXR4k9Jl1ISbnWRWPMWHldyNmhuUSn/eOg+zjvvPIQQ\nrHa9xzPvOXA2+CisXING2+cYAJACgUCn9ulUvgUqoCtDWk2zi0g0iN4USEkLtgTpXNJJqCUZyDVm\nA5qSaWRPP49Zx5+y4+XkrZcczZINdjZs9eGoD5E/YnOfOxHFIga0GgPZNnMfZyjK4FHNCOWg1Nje\nxWdrG+j0pAbqXXV6QyQSMbS6nbv5hDvCtLzdQvOrzYRaQmhMGnKOyqH88mEYh2fT6Xbw91e+3HF+\nQbaVf9yYHHYYD5Xgspf3Ok5XAkFPFkajhRHFajUMZf9TLXTloPLhsm08+PISNtS1EYmGMOhN5GTY\nOPqwUm69ZPaOzSG2y8tMQ6s1EIvqiXRG6PyqE/+27q4XLVBpovzoIjSGZNsms7CRjro83l6yhd9+\n71hMxuR/kSMqi/jDj07kFw+HaW1LEK0xk1Nag06/c1kAnzMXErkUFhRz+UkT98vzUJRdqYCuHDTe\nW7qVa+9/HaezgVCwE60+TCxqpLnJSK29hvnLa/j9lXM4a9boHddUluQgQl5869bT1daUPKgBKoxQ\nZQKjhnpncjGsiqIsjDYvGl0XDpeDt5ds4bxdNo04Z/ZYAqEo9zxjpKUtnZZNGRjSutAZQyRiOkKe\nAnLySrnq9MlkWE3789EoCqACunKQqGtxccu/5tPu2IbOVEvJyFo02jgJCdGAlU57BdV1Lq5/KIA3\nEObSkyZSXV3Nk/ffQcfH/wOZnBSUPi4db4kGYdaQkKkdJwIwWrsIh3ysr3Nw3rE90y89aSJTqor5\nzb8/5suNdYSDPiLRIDoERcX5nD5zLFefMWX/PBRF2Y0K6MpB4fYnPqGhaRtoWsgpr0bT/VZSI8Bo\n8VFQtRZXk5e29ji/ecTL/x77I2+9+j/i8ThCo0UUjCJ/lgZLQYTtC1D0NjkIwGAKEHSFqOlO392Y\n8lxevOMClm9uZkN9B9VNncTiCc4/dhyTq4oG6xEoyl6pgK4c8PzBCF9tsuPzOSgavWVHMN+VRoAt\nYyu+FVuobWqiVibQarVcddVVaIYfw0ufrSUQWIeFLXu9n94UwBsLUdfa9zroQgimjhmmJgMpBxQV\n0JUD3rItzfgDXrR6H3pj6trkMW8M1woXno0e6F7J0FQ6kSt/8ivm3XIZm+o7+GCtg+amNqKhhh27\nD/U1OUirjxKPx/AFe18HXVEOVCqgKwe8zQ0dybHkxmCP49GuKK6VLrybvMlALsBaZSVtdCldXdNY\ntMWLLxBmTHkuJ0wexaueVjoaRlEwah2aPQzYjUf1aLU6Mi3qxaZycFEBXdnvAsEIz3y0hk9W1WF3\neCjNy2DUsGyOHF/CqdNHpaw2qNdpEexcozniiuBa4cK3xbfjoHWUlaypWRiyDUAEr9+Nw9XO659t\n5tKTJvKby2ezfGsz1TU+OhsD5JbX9Fm+cMCKTm+mrCBjUOqvKINFBXRlv9rS6OSye16hsa0Zv6+D\naDTEKp0RvcHMv9/NZeb4Cv527Sk9gmlhthWtzkCsLUjr1tad48gFUGZAVJkpqMrtcR9LpgO/x8m7\nX1Zz6UkTGV6UxX3XnsLV9waxN0XoaguTUdCUUr5EXIvPWURmZjZTRqsXnMrBRQV0Zb9pcXq54s+v\nsaV2PeFIE7bcVjIsXiLBNCJBC+1t+cz3dXCWvZN515/O0RPLAHDat+Fa+gKh5o3JjDRAmRFRZUKm\naZCkjlgx2brwOYM0tnt23P+EyRXc+N1Z/OXZCO0OQdhvJXtYLbrufvmEBKe9Ar2ugNHl5Vx52uT9\n9mwUZSCogK7sF1JKrvrrG2ys2Uwkaqdo9M5lbI0WHwAZ+U046v1U1wW44WENd3yngocf/Bvvv/9+\nMhONFtvYdLKnpNPo9e7ItzdaQ4REPIrL17Pf/SdnT8Og03Lvf020d6TTvDkbrS6EzhAmFjGhEdnk\n55fz+x+esGOWqKIcLNS/WGW/WLymgbXb7Ph9doqqNvRYxnY7nTFC/si1NH/ZweqVb3LWY40AWK1W\nSibNoVGTizFvKzprExXWZEu8r7HkGl2EhIzjD4YJhWM7grMQgh+dMYXZh5Vz59MLWLqhkUg0SCwS\nQmc1Miy/gF9eeOSOTweKcjBRAV3ZL177bBNeXwdpGR07ujh2JRMSf40f90o3EUfyhaXQmzj1nMv4\nzyN/Yv6aFn4x7zU6nF5sBU29jkXfVfK9qkTKBP5QJKW1PaY8l+d/cz4ub5BNDR1saXRSlGPjhCOG\no9OpzUKUg5MK6Mp+sXJrK+GQh8xhPWdfJmIJvJu8uFe5iXmSKyJqzVpMI4cRzzwBT/6RZGZmcd4x\nWdz/0lLc7ibczeVkD6sH+h5LHg2a0WgMZNks5GSk9VmuLJuZI8eXcuT40gGqqaJ8e1RAV762eDzB\ne19W8+riTayvc6DXahhRnEVVSTZnzRrNhIqClGu8wQjxeBS9MbncbTwUx7POg3utm0QwORtIl64j\n8/BMbGNsoNXQtD5IfWsbn66t57jDK/jN5bP56QMumpvDGNM8WLJ6n5oPEOzKxmi0MrY8r89zFGWo\nUQFd+Vp8gTCX//FVvtpYg8/bTijsQyBYtdGM3mjmiXeWc/5xh/G77x/Xo5tDr9UghIaoN0bX8g48\n6z3IWPcLzUwtospM2eR8xC59KWabi4DfxYfLajju8ApOn1nFD0+fxj9eC+K0h9AZNux4obqraMiI\nxzGMvLx85kyuGPRnoigHChXQlX7zBcJc9odX+XzNOrrcNVhyWskr6gQg4rcQ9GXSaG/nybfdrNza\nypO3nE1Rjg2AqKeV6OZFtHy6decMoXwdotKMzNUihaCuLbl2yvZuFHN6F942HxvqHTvK8OtLZ7O2\npp1FK2O0b9NjyWkks7Bpx0vWaNiAo24MVlspMydU8b2TJ+2np6Mo3z4V0JV++8kD7ySDedc28ket\nx2jeuVuQ0eLDlt9GyNtMR32A5Rsj/PzvBq6ZncP999/PsrffTp4okrM6M4/IpDmWvL6voYcaTZSE\nTBCJ7hwRo9VqePrWc7jpUSuvf2ahs9NG47oStLowGm2MWMSK1ZrPqLJRPHLDXLRqb0/lEPKNAroQ\n4kbgKpJtrrXAFVLK0EAUTDmwrNzawqLV2+hy16cE812ZbF7yR6ym+fN6Xlv8JP+7sxUAnd6Atmg8\n+gobBeOTMzQrMAJ9Dz0U2jjIBJFYouc9jDoeuv40zjyqir+9+AU1zZ2Ew0Gi8QhppnRmTyrnlouP\nTtm9SFGGun0O6EKIYcD1wDgpZVAI8SJwEfDkAJVNOYA8N38tHm8bJpujz2AeD8bpWt+FZ62HeDBO\nHNCarPz42p9wzncv54f3v0NryyoS8dZex6HvTgiQSKKx3s89aepITpo6knAkxuZGJ3ZHFzPGluxx\nVIuiDGXftMtFB5iFEFEgDWj+5kVSDkSL1zYSDLjJLm1PSYt0RnCvcePb7EPGk90nhhwDmqIxmIpP\nwDjqBObMnMCYslU4nQ3/v707j5Lrqg88/r2vXu17dVX1qlZ3a99syRK28IbBSyyPgwGPSSYkEMg5\nTkjCgWEyGWclk5zJTjJkYYizkoQZEgIBBxIgECAE40UW2tXae9+79r3qvTt/VLXUavW+luT7OaeO\nqurd997tW0+/unXfXUgMtxFq6722/1xdDwsZH1Z94Umy7DadO7Y0cseWm3vXKMrrybIDupRyUAjx\nu0AfkAe+KqX86qrlTFlzlYrBN0/08p1TfQgh6GwO0NUSYl9nFJ/bfkPaRCZPpVLE7q5OjCWlJN+f\nJ3EiQb7/+vB612YX/jv9OFudlHKC8Sspvn7sKoZh8uFnDnOhf4ThoSxObwKnPzlv/nKJBrzuAPeq\nPuKKsigraXIJAk8BnUAC+IwQ4oellH87I92zwLMA7e0bN5y6XDFIZoskMwWS2SLFcgXTlBimRNME\nmhBYNIHHaSPgceD3OHA7rDdN5Xq7+MQXXuWvvnKCgdFR8oUUSIlVd6DbHIR8Qd79fXfygbffjd1W\nvUQMs1rzNo0KmfNpEicTlGPl6sEsQLudTYej2IK2a+ewuTOgZZhMxnntwhDfd/dWfuiRA/zFl7JM\n9udpdJycdcEKgPRElEopSKAlwjsf2j1rGkVRbrSSJpdHgKtSynEAIcTngHuBGwK6lPJ54HmAQ4cO\nzd6dYQ2YpmQklmF4Mk0yWySVK5DL5sjlsmSyWcrlMqZpIk2JECCEhsWi4XS6cLldeNxu3C4nfreD\nkM/JpogPt9O28IlvAb/xqW/z8c+/yOTEZSQJ7O4kwmKSLTmppB1MTvj56KfH+OrRy/zBTx9hR3sD\nxdQERs9r9B/twSxOLQskEF0OZIcNbBqDhSwMZ681oQjA5sxQKGQ4en6Yu3e18ZH3PMSxi8O8eqbA\nyAWNQMsVPA0TTP/azMaDJIa3EW3cxo89cRchn2oTV5TFWElA7wMOCyFcVJtcHgaOrkquVqBQqtA7\nkqB3NMn4ZIyJiUmy2QzFYhGnVcNl1/DZdWyuaq18qgJuShPDNMgVY8RS4/QVDUw03G43fp+fcCRM\na9hPR1OAaNB9y9bcP/r3L/KJz7/I2Oh5/I0X8USGb5gXRQK5eIj4UIqXT8R58v2v0mVepO8rX6F6\nixLsUTv+fX7GPWWEJubsdggghFnbr9a2brXwyefezvt/38pLZ1zERtykRpPotgK6rUAp76ZSCtIQ\n7uDKYa4AACAASURBVOQH3nKQn/mBe9eyOBTltrKSNvSXhRD/ABwDKsD3qNXEN0I8nefSYIyhiRTj\n4+OMjo4hzBIRr43GBh2nzYa2iCAcdFuvPS9XTDLFCvHECCcHB+gLhbgSjRJpCNLRFKCzKbBh/Zyl\nlHzjez288OJ5Lg/GSGQLhDxOwgE3B7c38+7H7sDjurEd/NXuQf7wsy8yMnoef+MFfNGRm44rAIdr\nHHe5h/ireVK5DJcBzaKjRbbg3WsltL3aZu6t7TNXt8PrRxRYpq35Fg26+fuPPMMn/ukof/LCa8TT\nCcrFPOVKAZ/XTcAb4u0P7ObX3vfmW/aLU1E2wop6uUgpPwJ8ZJXysiyGYdLdN8H5vjEGBgeJTU7i\nd1roCNnwOb0LH2AeVl0jqGsE3VbKhslEOsvli9309DoZaG6mr6WJ/VubCPmcq/TXLE7PcJwP/uGX\n+d7FPtKZMUqlHIZRxmKxYtGsfOnFIJ944SjvfXw/H3jH3ddmD/yTf3qNZGoYp29w1mBeipVInkqS\nPp++Pizf7sK75Y3c+cBTdA+OUCidAs4uOq9SCoQQ6DOmR7RYNH7qbXfzY0fu4mzvOKevjnFlOM7O\nTQ08emiL6nqoKMtwS48UjaXyHL80Qt/gML29vUQ8Fu5oc2PVV7/WbLVoNAfsNPltJPMV+vuuEpuc\nJJHOsqM9ys728LrU1k9eHuVHfv1z9A1dIpcdxhUcxxdKoNtKVEo2KkUHqXiURGKQ3/n0OF966SJ/\n/fNvI+B28B+n+shlJmncdn3pNWlKslezJE8nKQxeHxPmbHXi2+dDBNuIDbTRO1nA6QySSgQpZt3X\nervA3N0OJVDKe/C5XXS1hGZN47Dr3LW9mbu2q+XeFGWlbsmAPr1W3tPTQyGbYnujE49j7f8cIQQB\nlxWfU2conuf06dOkUu2MxrNrXltPZgr81Me+xNW+cxiyn+ad59Gt13uJ2F3VIOtvHCYbDxMbSnOs\nkOE9vyn48ScPksrEsdgy2JwFyuky6XNpUudSGNnqwB2hC7w7vfj3+muLLQMkyMbHyeVjdLR2Uihs\nYqI3R/OO42iW+e9xF5I+wEskGOJBtWCEoqy5Wy6glysGr5wb5GLvIL09PUS9OlvaPItqH19NmhC0\nhRwEXBWu9vcQi8VIZfO8YWcbrRHfvPsahsnnvn2Orx69wpmrYySyBTQhCHod7OmI8sC+dp5+cPdN\nizJ8+ONfofvqJcqVYZp3zL7qzxR3cAK7K8XIJTh+3spvZ6ur8shEL8NfHCbXd320pzVgpdJuhXYb\nkfaGm47laRgjPjiBbuliR0cXJ84lmezfQrjjEnOVumlYiA934vVG+f437lCLRijKOrilAnqxVOGl\nswNcuNLLyNAA2xtdK6qVSynpj+U5N5hjOFms3b4Dm66xpdHJjmYPPuf8x/c4dPa0eRiYzHH27FkM\nw+TQzjY2NwVmTX/swjA/+4l/pbu3n3RmnFIxS6VSQgiBxWLl5Hk3L/xHiN/7zHd575EDvP+th7BY\nNPpGk3zz+GWSiQGatnUvaui8bi8R7TrLyJkCJ7u/Q7H/NSjWppvVwLPFg2+3D0eLg56RxJzHcfji\nmH1ZRiaT/NaPP84vTCQZGi4ydlkn3H4Zi7V8Q3rTsDB6aTdWyyY6Wjt4/1OHFsyroigrd8sE9HLF\n4LtnB+i+eJWJsWF2tXiwW5fXZh3LlPjS8XFe60mRyuUpl7JUKtfbj4VmwWp1YrO6ifodHOz08fDu\nBvyuag+YYsXg/HCOC0NZ8mUDIQROq0bQrVM8dRqJRAhx05D1r7xyiQ9//Mv0D16gWBrDExrD1xzD\n6sgjqa6yU0gHiMcjJOID/ObkOF8/doWP/fQRPvuts6TSk9hdCWzOGxc+no00Jbm+HKkzw5T7Tl+b\nsla43IT2O/Du8GJxWqo9VEYKmLWuh7P1WNEE6LY8pXKOoNfJbzz7GL/yV1aGRrwMdftweBLYnFks\ntiKlnJtcMozd1sz2rt38zc+949oUuoqirK1bIqBXDJOXzw1y4XIvk2PD7G5Z3o1PKSUvfG+Mr54c\nJ5GeoJCLIcmhO9LojhTVqCcwDSvZvI9U0kM87qJvNMC/ngpyV7sfAzjRlyJfzFEu5TDNMgjQhBWr\nzYWuO2g5NsLTb97P+544SEu4Gsy6eyf44B/9C/2DZ9Fs/bRuuXhTLdviyeDwZPA3D5CZDDM+lOFb\nR5M8/ZEUTUEPhUISV3By/rLKVEidS5E6e71tHA0s0XaMxnsQgRyBvd1LLrvpZfjMQ3vY19nIz//Z\n13mtu4dcIUk5kydfKWK1OQk3hNjc3MSf//en2Lk5vOxzKYqyNLdEQD9xaYQLV/oYHR5gd6tnWcHc\nME3+9BsDvHRxhHRqEM02hrflCjZXgrma301Do5gJk0+0kR5r44vxCMJII4w4Fj2Dbk+iOQqAwCjb\nKWb8GBU36bSPj346wZdfucLz/+1JdrSH+a8f/zIjo1fQbP1EOrvnXeRYAN6GCZyeFGM9BS5cLXO2\nN0qlmCPoSd2UXhrVniqjJ2Iwer35w+q34tvtw7vTS6kSZKh3K2aph0pFoOvVGvlUTXz+vuTTzlWr\nye/cHOazv/pOvnOqn+OXRujun2BgPMWWliBv3t/BkXu2qbnIFWWd1X1AH55Mc3lwnL6+Xva0urEt\nI5iXKiYf/1of37s6TDrVhydyFmfg5n7YM2kWE6d/DKHbSUxsRRYGkNoIItiN3ZPBa9fRZkTmSslB\nLtZOKtnCsXNxnvz5FE+/aRenLveRL4zQsvPigivWT9HtJZq2nmHkoqAovciKgSmvzw1eHC+S6k6R\nuZC5PhxfA3enG98eH85W57WBOU5SaFoRU3MRGw0SbY0tLhM18todhuuEENx/Rzv3qx4silIX6jqg\nl8oGJy6PcPXKFdqCdhzW5fWU+IdXRjh2dZhMqgdf00ns3sUHs0rJSzJ+D7IYA/tlCLyGlFYKJZ2K\nIQm6rDcEdd1WwNd0AVfoKomhvQwM5vijz6Uwy2O4Q6NY9IVvZk6nWQwat56h90InEkF8rIgzlSB9\nLk1pctrEVn4LtNtgk4283UKeIp3ixsE5dscIeaORXKwNszmOpl3vdjhfzdw0LFSKHlwOLwe3tywp\n/4qirJ+6DuinrozS1z+EZhSI+pa3+szl0SzfODtOJjWEd4nBXEpIJQ5jFrOg96E1vFxddEEamBU7\nZSCe46agDqDbygTbj5Eci1HIRqCcxOnoZWYtdzE0rYyeOU35Sje5eC+5WhzWHBrebV5SYYkI6Ndu\nbM7F6RyiWOhEGBHiA500tF9Z1PnTE03Y7T72dDQSUasAKUrdqtuAPtXUMjQ0wN5Wz7Lm9JBS8tf/\nMUQqPYzNNYBjCcEcoJhvo1zwQXkY0XT0Wlu7EBJNL2JWmDeoWzSB05egmM+ByJIpJ7BmqzM4LkYp\nXiLdnSZ9Po2Rux58LS02wvuCuDvcCItg6rbjQu3gui2PTWSxuFrJJ8tkJlN4GibmzUO5aCM13kw4\n3MgPvmXvovKtKMrGqMuALqXkXO8EPVd7aAval9098cxghv6JJOXSBKG280vev5DrxCznEO6raPqN\nS6VOD+oVIJUX+F36zV880guiBJYcRqVMIiOw23Qc1tmLvpKrkLmYIXMhQ3G8eO19zeNCRvYiww5E\n+xnsDU6EZWlfchbdQCPPpsYtpDMWxofKmIaONzoy6++GctHG2OW9eLyb2dvVyTNqXnJFqWt1GdAn\nkjnGJhMU8xmikeX3Yf7OhQT5fBybdwzLIgbiTGcaOsVCI1QmITh708T1oC4oVgS5ksA9Y3SnQCCQ\nSGEgtAqGIZhI5GgJe6+NbjXLJtmrWcZOxWGsfK3PuGbTcG9x493pRW8IMNhzL0Z6jErxIqlciYYZ\n0wws1EPFKNvQLDp3dDUS8nbwxe9qTE44ycYjeMLD2Jx5rPYcZsVKNhEmPdGM27WZA9v38H9/8R1q\ntKei1Lm6DOg9IwnGxkaJ+mzLnj61WDE41Z+iWEzhjwwuvMMM5XIIDBNEGs2WnjOdEBLNUsI0BNmi\nwGrRZumJI0AKhKWMWdEoVyrEU3mcSUhfSJO9kr0+u6EAV4cL73Yvrg4X2rVjFXF6xskU3ZDeRsm9\njF8cGS82u5v9W5t47ofu5+5drfz+Z15ieHyA7EQLyUoRo1JCCA27w0MoFGVH+2Y+9YvvULMfKsot\noO4Cer5YZnA8SSwW44625d+AuziSI5PPIrQ0Nmdmyfsb5QBSVsB6c7/vmYRmgFnBMAXpgiDkvr50\nnW5NIjQdaXiQ0kSki5g9JRKDCRKFaTcxQzq0WaHVRsFhoUCZzhlfDIFQN7n0A5iZTRRTk9BwYzPQ\nfEwTCukgTY0h3nygEyEE7z1ygKfu28Eff/5VTl0ZY2gyzWgsg1XX2LU5wtsf2MUzb9p9bRk6RVHq\nW939T+0dTTI+Pk7AaVnRNLi9E3kqlQK6fe7a9XwMw1mtoVsXHmYPVGvpFY2KYaFQNnHaqs0Tmp5C\n5OPQfwl5GkhXru/jseDf5SMRMBEey4K9VOyOBL5gL4lSJ0b8DjKhM3j82Xn3mZIc2YRuDdDeHOGe\nXa3X3g/5XPzSu9907XWpbCClVEFcUW5Bdfe/tn+sGtA7gitbv3MsVcKoFLE4cgsnno2QS+phKAQI\nrYJpamSLGnrRpNhTJH81jxGbtsyqDWjToNWBpdFFKOpnaqbwxYzWDEVOkp5oQ4oWYr1g3Xoau2v+\nvzGXDJKZ3ExTcycfevrwvM1YtmX29VcUZePVVUAvlCqkcwVKxQIex/xT0C4kmatgmmVs1sU3S0wn\nRAUQYC4+wMl8Efol5aEcE/FpN2F1HUKdEHUgtp9Bs2gYJYFpmuSKZVx269wHvSlfoDOEsO/DoXcw\ndlnHHRzE3zxw041fU0JmvInkaCcN4S7e9chBnn6T6qmiKLerugroyUyBXDaH225Z8VqSK50eXROV\n6jSD5vy/FGRBIgeqDyYAaqM3LeBod+DscGJpChAbP4KZGUeWB8CSRtNMpDTJFa4H9IV6qUwRaAQ9\nOnfv3sH3LgaIxwIMnWnG5kmgW4tolgqVsp1S1ocgQLSxk8fv2cuvvveh5ReIoih1r64CeiJTIJvN\n4LKv/Gd/dQ3Lau+SZe1vm0RoVmTJf9M2WZDIoVoQH5u2QQOagGY7NDnwNThqiyMXsDuGyBcDkLwD\nHN8BzUAaBqWKedPx52OagkrZhtPu5pPPvY2j54f4g8+9zPFL/RTyKSqVEmbJwKrb8IY8tESj/MT3\nH+S9Rw6oBZcV5TZXVwE9mS2SzeUI2VYe0C2aACGQcnnHmkz3V/uDm25Mw4oolpCDEjk4VROvEUAj\n0Ao0A1YBZvVmbqFs4rZXn3v8JygWHsfMtmJmm7G4h6lUqvO8SxbfXJ9PBrBa3bRFA/g9Dh4+2MXD\nB7s4dmGY45eGGYllmEzlaQp52dcZ4ZGDXar/uKK8TtRZQC+QzWTZ1LSyG6IALpsFTdMxKss7ltBM\nZKEfhvuRZyUyNa0mrQFREG0C0SIQtmo4NoxqGs1iYJomxYrEba/uolszuDyXyFQ6kPF7MPV/B5EF\nTCqGgdWyuKCbmWzG7Wng0UNdN7yvFlpWFGVFAV0IEQD+DNhLtT77Pinld5dzLNOU5AplyuUSdn1x\nc53MZ1ODA6vuIF+cf6TpRLLaQyTsdyGlZKQnBcOV6iP9hesJLSCaBbSCaBII63wTmpsgq3OwT+f2\nnaBUDFKSYczxByD4IlIvIs3q8ReSTQYp5cM0NkT50cf3L7yDoiivKyutoX8M+LKU8j8LIWzAsocT\nGqaJaZpUW0pW3tbbGXWi6y6M7MJZMmMVUhdSFPoKkJ7eO0VAeAd4w7DzOJp9/j7pUws6SCRIE9OE\nimmia9X3NU3iDX2DNI9SERHM2H2YgdOYoYX7upfyTmL9Wwg1bOZHHttP2wILUSuK8vqz7IAuhPAD\nDwI/CiClLHGti8fSGaZESvPa/CYrtSnkwGm3EzedGGUbFuuNWZuqmZfLBpwtkB2vBnLNoeFod5AL\nSQhbcLCfQtJVDfT2VxZ1bgGgSarNKZLp46M0i4mn4Rvk4g9TMFqRCScJxmhoGcDqmL2LZTYeJDaw\nDZ+/gwM7tvJz73pgqcWhKMrrwEpq6J3AOPCXQog7gdeAD0opbxi6KIR4FngWoL19/Va2sWgam8NO\nJmJu8slmPOHeuRNv0tF8FkSLlUinF6EJ8rHqdAFu90lK+Ucws5sxS+fRbMlF5kAiqX5RzaRpFYKN\n32Q09STC1YluRhi50IhuT2Fzp7Ha8gigVHRRzPoxjQANDe0c3ruDv/wfT6nBP4qizGolAV0H7gI+\nIKV8WQjxMeA54JemJ5JSPg88D3Do0KE5x7ZbNIEQ2oLD35fiDV1+TvcFyGUabwro4dpkUxPJHGy2\nXHs9pSnkqT2LY3cNUig3IGN3Y0a/dsNKP3OZ+qFxU8palxaz5MZKji2bdnJ490G+dfwq2XyCQiFD\nIVedNle32gn4vYQDEX740X188OnDaki+oihzWkl0GAAGpJQv117/A9WAviyaEGiaxiwV2mW7Z0uA\nz7zsJ5P2Uym40R2Lm/dkJm/gGOXiESrZRmTiAISOLW5HyU0RvfpSkI834nQGuG9vO8//zPczGsvw\n7yd7OXp+iFi6gJSSaMDNG3a28PBdXfimussoiqLMYdkBXUo5IoToF0LskFKeBx4Gzi73eBaLhsOm\no1utFMvmshe1mM5lt7C71UM86Scb24S/pfumNDNr5rPnrYAv9CKJygOY+W2Y6RSa99ICe80Szaku\naweCQrqBSKSBpx/cBUBjyMMzD+3hmYf2LPyHKYqizGKlUfMDwKeEECeB/cCvr+Rgfo8Dt8tFtri0\nxSjm8/idYdzuCKVMC0Zp+d0h7Y5R3P5uNGcDMnUQM7V93vRSCoQQaDNKWAKFWBe6Hqa9sZFHDnbN\nur+iKMpSrSigSymPSykPSSnvkFK+TUoZX8nxAh4HLrebXGn1AvrWRje72/zYHSHSY1tXdCyP/zQu\n/8VqUE8fwIgdwDTn6JVTm3LAMi2iSykpZf0UE5sJN3TwS+9+8FpXR0VRlJWqq2jid9vxuN1kCqsX\n0AHeeU8zPm8j5XwzxXRo4R3m4fWfqNbUXQ1Q3IMcfQyz2HBDGgkgNUDDOm3h6EI6TG70IL5AB08c\n3s2Re7atKC+KoijT1VWXiYDHgcu1ujV0gLaQgwd3hvnKiQyp0Twhxys39UtfCo//NFb7OOnYYSrF\nNuREAEMfRXiugH2kuoKRENX5ZBAUsz7y8XZK+WY83jYe3L+Hj33gyOr9gYqiKNRZQHfarXhcDiy6\nlWzRwL0Ksy5OeefhJi6P5TjXmyU5tJfApmM3tW8vhd0xit74z2SSByjkNiPLHmSqDYwSUmRBmJia\nYMJ0o2lObA4f/kAD9+3dzO++/zHVl1xRlFVXVwEdoC3ioycSYSw1Rmdk9RYmtlo0fvwtm/i1LxQZ\nGc2THt2Br+n8iuZNt1jK+EOv4PGdJJfZSanYRKUYxSzrCM2C16ljszoIepx0hp1saXLz2JvuVsP2\nFUVZE3UX0Dc3BTgfiXByeIhNIYluWb05vBv9dt5zfyvP/1uZeNwkNSLwNXWveDEMi17AGzgOQKxv\nH1h2saW1mXfd20zUZ6PRb2cgVsC0B2iL+NXgIEVR1kTdRRaP00Zz2E9fIMhEJkeTf3UH1LyhK0Cp\nsplPflsQiwlSQwJv87kVNb9MKaQiGMVmopEIP/lIG+0N1V8YppSMpUrs3tNIR1Ng5SdSFrTnI1+e\nt/ur227hzP98fB1zpChrr656uUzpaAoQjTYyllz+jcv53Lc9yPvetJmGUCeV0hbifW+gUnCu6JjF\nrI/06C68vhYe3Ru9FswB4pkyLo+XaMhHyLey8yiLs9BYhtUc66Ao9aIuA3pTyEOkIYBmtTOZXpug\nfnhrgJ98tJP2pq3Y9a3E++8hG2vDXNqKcAAU0yFSQ3fh8W5mV1sjbz0YvbbNlJKhRJGmRlU7VxRl\nbdVlQBdCsL2tgY7OTvpiBcpLXHdzse5s9/KRt2/l0LYO/IGtFOJ3ELv6RnKxNkxj4aIxyjaSQztJ\nDR/A6+vgrq42PnykA+u0wUJD8SIOt4/Wpqi6Gaooypqquzb0KZubAgzHosRjcXomxtnW5F6T8/hd\nVj70+Ga+ec7Ll457GU/GySXDZCa2oDtSWB1JrM4UQitXdzAtlPIBygU/RiGIzRGkIRzlvu1hfuT+\nlhtGhmaLBmPpCvv2dbB/a5MaFaooypqq24AOcOeWRiaTOU6cjDOZKdHgWflao7MRQvDm3Q3ctz3I\nt7pjfPt8mKFYllIpQ6mQI5suMH2iLd3qxGl1YXN72Nni522HomxtvPELx5SSK2M52ts72NEeVW3n\niqKsuboO6E67lb2dUVKZLi6eP4fXoWPT166Wa9M1Ht0b5tG9YUaTRc4OZrg0mmMoXsSUcmoqc1qC\nDrY1utjT5qFxjl44U00tm9ua2dkeXrM8K4qiTKnrgA7Xm15SySQXhofZ2eJZ1b7pc2n022n023nz\n7oaFE88wnioxkZXs2aOaWhRFWT+3RKQ5sLWJbV0d+EIRLoxkZ13WrV5MZkoMJCrs3LmTA9tbVVPL\nBllo2ojVnFZCUepF3dfQAew2nTfuacOUkrPdBueH42xvcq9LTX0pJtIl+uNlduzYwZ3bWulsDm50\nll631KAh5fXolqihQ7U9/Y2729i9czueYITuoQzF8tp0Z1wqKSXDiSL9iQo7d+3izu2b2L5p6U01\niqIoK3FL1NCnuJ027t/XjkUTXLpq5czgIG0hG1Hfxq23WSybXB3PYVoc7N69g7u2t9HVomrmiqKs\nv1sqoAM4bDr37W3H67LjDwS4cvkKsUyGzohrVdYhXSxZm59lIF6kpaWN9k0t3LmlieYG77rlQVEU\nZbpbLqAD2KwWDu1ooaXBi8/jord/kDODQ7QEbUS8ttrCEmsnVzTom8xXa+V79rKlNcK+rkY1x7mi\nKBvqlgzoU1rCXhp8Tk753VwJBunv62OoL0XYayXqs+FYxQBrSkk8U2Y0VaJoajQ3tbCprYU7uppo\nCatauaIoG++WDuhQ7QFzaEcLrWEvl8IBRieTjI2Nc3ZoHJcVoj4bXqflhvlVFktKSb5kEsuWGU+V\ncLo9NLW1EgmH2BT1s72tQc1trihK3VhxNBJCWICjwKCU8smVZ2l5mhu8NDd4SWQK9IxE6R9LMjEx\nydjkBFcns1gwcdstuO0WXHYLNl3DIgRCgJQgkRgm5EoG2YJBrmSQK5nYbHYCgQA790RpDPnpaArQ\nGvGhq8FCiqLUmdWoXn4QOAfUxVSCAY+D/Vub2NMRoX8syvBkmmS2SDaXJ5vLks1mGclmKZXKSGli\nmiYIgaZpWDQNp9OJO+Cmwe3B5XLidTkI+ZxsbgyoQUKKotS1FQV0IUQb8J+A/wV8eFVytEqsuoWu\nliBdLUGklGQLZRKZAslMgUSmQKliYJgS05QIARZNw6IJPE4bfo+DgMeB323HqqsbnYqi3BpWWkP/\n38DPAnPeFRRCPAs8C9De3r7C0y2PENVA7XHa1JzkiqLctpbdECyEeBIYk1K+Nl86KeXzUspDUspD\nkUhkuadTFEVRFrCSO3v3AW8VQvQAnwbeIoT421XJlaIoirJkyw7oUsqfk1K2SSk7gB8E/k1K+cOr\nljNFURRlSVTfO0VRlNvEqoyKkVJ+E/jmahxLURRFWR5VQ1cURblNqICuKIpym1ABXVEU5TYhpFy/\n9TmFEONA7xybw8DEumVm8VS+lkbla2lUvpauXvO2lvnaLKVccCDPugb0+QghjkopD210PmZS+Voa\nla+lUflaunrNWz3kSzW5KIqi3CZUQFcURblN1FNAf36jMzAHla+lUflaGpWvpavXvG14vuqmDV1R\nFEVZmXqqoSuKoigroAK6oijKbWLdA7oQ4nEhxHkhxCUhxHOzbBdCiD+obT8phLhrHfK0SQjxDSHE\nWSHEGSHEB2dJ85AQIimEOF57/PJa56t23h4hxKnaOY/Osn0jymvHtHI4LoRICSE+NCPNupSXEOIv\nhBBjQojT094LCSH+VQhxsfZvcI59570W1yBfvyOE6K59Tv8ohAjMse+8n/ka5OtXhBCD0z6rJ+bY\nd73L6++m5alHCHF8jn3XsrxmjQ31cI3NSkq5bg/AAlwGugAbcALYPSPNE8C/AAI4DLy8DvlqBu6q\nPfcCF2bJ10PAF9ezvGrn7QHC82xf9/Ka5TMdoTrwYd3LC3gQuAs4Pe293waeqz1/Dvit5VyLa5Cv\nxwC99vy3ZsvXYj7zNcjXrwA/s4jPeV3La8b2jwK/vAHlNWtsqIdrbLbHetfQ7wYuSSmvSClLVBfG\neGpGmqeAv5ZVLwEBIUTzWmZKSjkspTxWe56muuh161qecxWte3nN8DBwWUo51wjgNSWl/HcgNuPt\np4BP1p5/EnjbLLsu5lpc1XxJKb8qpazUXr4EtK3W+VaSr0Va9/KaIoQQwDuB/7da51useWLDhl9j\ns1nvgN4K9E97PcDNgXMxadaMEKIDOAC8PMvme2s/l/9FCLFnnbIkga8JIV4T1fVZZ9rQ8qK6uMlc\n/9E2orwAGqWUw7XnI0DjLGk2utzeR/WX1WwW+szXwgdqn9VfzNF8sJHl9QAwKqW8OMf2dSmvGbGh\nLq8xdVN0GiGEB/gs8CEpZWrG5mNAu5TyDuAPgc+vU7bul1LuB44APyWEeHCdzrsgIYQNeCvwmVk2\nb1R53UBWf/vWVd9cIcQvABXgU3MkWe/P/P9QbRbYDwxTbd6oJ/+F+Wvna15e88WGerrG1jugDwKb\npr1uq7231DSrTghhpfqBfUpK+bmZ26WUKSllpvb8nwGrECK81vmSUg7W/h0D/pHqz7jpNqS8ao4A\nx6SUozM3bFR51YxONTvV/h2bJc1GXWc/CjwJvKsWCG6yiM98VUkpR6WUhpTSBP50jvNtVHnpy4N5\ndgAAAYBJREFUwDuAv5srzVqX1xyxoS6vsfUO6K8C24QQnbXa3Q8CL8xI8wLw7lrvjcNActpPmzVR\na6P7c+CclPL35kjTVEuHEOJuqmU3ucb5cgshvFPPqd5UOz0j2bqX1zRz1pw2orymeQF4T+35e4Av\nzJJmMdfiqhJCPA78LPBWKWVujjSL+cxXO1/T77m8fY7zrXt51TwCdEspB2bbuNblNU9sqMtrbM3u\nts5z1/gJqneKLwO/UHvvJ4CfqD0XwB/Xtp8CDq1Dnu6n+pPpJHC89nhiRr5+GjhD9U71S8C965Cv\nrtr5TtTOXRflVTuvm2qA9k97b93Li+oXyjBQptpG+WNAA/B14CLwNSBUS9sC/PN81+Ia5+sS1TbV\nqWvsEzPzNddnvsb5+pvatXOSasBprofyqr3/V1PX1LS061lec8WGDb/GZnuoof+Koii3CXVTVFEU\n5TahArqiKMptQgV0RVGU24QK6IqiKLcJFdAVRVFuEyqgK4qi3CZUQFcURblN/H/xyOAv0FH2ggAA\nAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Final P: [ 0.00972677 0.0187833 0.00070503]\n" ] } ], "source": [ "#format the book\n", "%matplotlib inline\n", "from __future__ import division, print_function\n", "import kf_book.ukf_internal as ukf_internal\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from filterpy.kalman import MerweScaledSigmaPoints\n", "from filterpy.kalman import UnscentedKalmanFilter as UKF\n", "\n", "dt = 1.0 # 采样时间\n", "wheelbase = 0.5 # 汽车前后轮轴距\n", "landmarks = np.array([[5, 10], [10, 5], [15, 15]]) # 给一组标定点位置,其实就是实际中观测对象的观测点,传感器的位置\n", "cmds = [np.array([1.1, .01])] * 200 # 给一组控制,其实就是实际中不断地输入控制量\n", "\n", "ukf = run_localization(\n", " cmds, landmarks, sigma_vel=0.1, sigma_steer=np.radians(1), sigma_range=0.3, sigma_bearing=0.1)\n", "\n", "print('Final P:', ukf.P.diagonal())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "好了!大概就这样了,其实思想很简单,就是通过描述sigma点来优化对状态的描述,基本过程还是KalmanFilter那一套。" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.3" } }, "nbformat": 4, "nbformat_minor": 2 }