From 6d9c500e50e305124fb6739c9608fffe58fea24d Mon Sep 17 00:00:00 2001 From: Roger Labbe Date: Sun, 10 May 2015 19:12:44 -0700 Subject: [PATCH] Restructuring for print. I have moved less pertinent material to appendixes. The idea is that a print edition should contain only the most important material, and online appendixes allow the reader to access supplemental material. --- 10_Unscented_Kalman_Filter.ipynb | 543 ------- Appendix_B_Symbols_and_Notations.ipynb | 854 +++++----- Appendix_C_Walking_Through_KF_Code.ipynb | 2 +- ...pynb => Appendix_D_HInfinity_Filters.ipynb | 720 ++++----- ...> Appendix_E_Ensemble_Kalman_Filters.ipynb | 0 Appendix_F_FilterPy_Code.ipynb | 1379 +++++++++++++++++ pdf/merge_book.py | 111 +- table_of_contents.ipynb | 38 +- 8 files changed, 2240 insertions(+), 1407 deletions(-) rename 16_HInfinity_Filters.ipynb => Appendix_D_HInfinity_Filters.ipynb (99%) rename 17_Ensemble_Kalman_Filters.ipynb => Appendix_E_Ensemble_Kalman_Filters.ipynb (100%) create mode 100644 Appendix_F_FilterPy_Code.ipynb diff --git a/10_Unscented_Kalman_Filter.ipynb b/10_Unscented_Kalman_Filter.ipynb index 2d39c69..0819351 100644 --- a/10_Unscented_Kalman_Filter.ipynb +++ b/10_Unscented_Kalman_Filter.ipynb @@ -2989,549 +2989,6 @@ " self.P = self.Pp - dot3(K, Pz, K.T)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Full Source from FilterPy" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Without further explanation, here is the full source from FilterPy.\n", - "\n", - "** author's note ** this is somewhat dated, but while authorship of the book is in progress I am not going to update this section every time I make a minor change to the filterpy code." - ] - }, - { - "cell_type": "code", - "execution_count": 44, - "metadata": { - "collapsed": false, - "scrolled": false - }, - "outputs": [], - "source": [ - "# -*- coding: utf-8 -*-\n", - "\"\"\"Copyright 2014 Roger R Labbe Jr.\n", - "\n", - "filterpy library.\n", - "http://github.com/rlabbe/filterpy\n", - "\n", - "Documentation at:\n", - "https://filterpy.readthedocs.org\n", - "\n", - "Supporting book at:\n", - "https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", - "\n", - "This is licensed under an MIT license. See the readme.MD file\n", - "for more information.\n", - "\"\"\"\n", - "\n", - "# pylint bug - warns about numpy functions which do in fact exist.\n", - "# pylint: disable=E1101\n", - "\n", - "\n", - "#I like aligning equal signs for readability of math\n", - "# pylint: disable=C0326\n", - "\n", - "from __future__ import (absolute_import, division, print_function,\n", - " unicode_literals)\n", - "\n", - "from numpy.linalg import inv, cholesky\n", - "import numpy as np\n", - "from numpy import asarray, eye, zeros, dot, isscalar, outer\n", - "from filterpy.common import dot3\n", - "\n", - "\n", - "class UnscentedKalmanFilter(object):\n", - " # pylint: disable=too-many-instance-attributes\n", - " # pylint: disable=C0103\n", - " \"\"\" Implements the Unscented Kalman filter (UKF) as defined by Simon J.\n", - " Julier and Jeffery K. Uhlmann [1]. Succintly, the UKF selects a set of\n", - " sigma points and weights inside the covariance matrix of the filter's\n", - " state. These points are transformed through the nonlinear process being\n", - " filtered, and are rebuilt into a mean and covariance by computed the\n", - " weighted mean and expected value of the transformed points. Read the paper;\n", - " it is excellent. My book \"Kalman and Bayesian Filters in Python\" [2]\n", - " explains the algorithm, develops this code, and provides examples of the\n", - " filter in use.\n", - "\n", - "\n", - " You will have to set the following attributes after constructing this\n", - " object for the filter to perform properly.\n", - "\n", - " **Attributes**\n", - "\n", - " x : numpy.array(dim_x)\n", - " state estimate vector\n", - "\n", - " P : numpy.array(dim_x, dim_x)\n", - " covariance estimate matrix\n", - "\n", - " R : numpy.array(dim_z, dim_z)\n", - " measurement noise matrix\n", - "\n", - " Q : numpy.array(dim_x, dim_x)\n", - " process noise matrix\n", - "\n", - "\n", - " You may read the following attributes.\n", - "\n", - " **Readable Attributes**\n", - "\n", - " Pxz : numpy.aray(dim_x, dim_z)\n", - " Cross variance of x and z computed during update() call.\n", - "\n", - "\n", - " **References**\n", - "\n", - " .. [1] Julier, Simon J.; Uhlmann, Jeffrey \"A New Extension of the Kalman\n", - " Filter to Nonlinear Systems\". Proc. SPIE 3068, Signal Processing,\n", - " Sensor Fusion, and Target Recognition VI, 182 (July 28, 1997)\n", - "\n", - " .. [2] Labbe, Roger R. \"Kalman and Bayesian Filters in Python\"\n", - "\n", - " https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", - " \"\"\"\n", - "\n", - " def __init__(self, dim_x, dim_z, dt, hx, fx, kappa=0.):\n", - " \"\"\" Create a Kalman filter. You are responsible for setting the\n", - " various state variables to reasonable values; the defaults below will\n", - " not give you a functional filter.\n", - "\n", - " **Parameters**\n", - "\n", - " dim_x : int\n", - " Number of state variables for the filter. For example, if\n", - " you are tracking the position and velocity of an object in two\n", - " dimensions, dim_x would be 4.\n", - "\n", - "\n", - " dim_z : int\n", - " Number of of measurement inputs. For example, if the sensor\n", - " provides you with position in (x,y), dim_z would be 2.\n", - "\n", - " dt : float\n", - " Time between steps in seconds.\n", - "\n", - " hx : function(x)\n", - " Measurement function. Converts state vector x into a measurement\n", - " vector of shape (dim_z).\n", - "\n", - " fx : function(x,dt)\n", - " function that returns the state x transformed by the\n", - " state transistion function. dt is the time step in seconds.\n", - "\n", - " kappa : float, default=0.\n", - " Scaling factor that can reduce high order errors. kappa=0 gives\n", - " the standard unscented filter. According to [1], if you set\n", - " kappa to 3-dim_x for a Gaussian x you will minimize the fourth\n", - " order errors in x and P.\n", - "\n", - " **References**\n", - "\n", - " [1] S. Julier, J. Uhlmann, and H. Durrant-Whyte. \"A new method for\n", - " the nonlinear transformation of means and covariances in filters\n", - " and estimators,\" IEEE Transactions on Automatic Control, 45(3),\n", - " pp. 477-482 (March 2000).\n", - " \"\"\"\n", - "\n", - " self.Q = eye(dim_x)\n", - " self.R = eye(dim_z)\n", - " self.x = zeros(dim_x)\n", - " self.P = eye(dim_x)\n", - " self._dim_x = dim_x\n", - " self._dim_z = dim_z\n", - " self._dt = dt\n", - " self._num_sigmas = 2*dim_x + 1\n", - " self.kappa = kappa\n", - " self.hx = hx\n", - " self.fx = fx\n", - "\n", - " # weights for the sigma points\n", - " self.W = self.weights(dim_x, kappa)\n", - "\n", - " # sigma points transformed through f(x) and h(x)\n", - " # variables for efficiency so we don't recreate every update\n", - " self.sigmas_f = zeros((self._num_sigmas, self._dim_x))\n", - "\n", - "\n", - " def update(self, z, R=None, residual=np.subtract, UT=None):\n", - " \"\"\" Update the UKF with the given measurements. On return,\n", - " self.x and self.P contain the new mean and covariance of the filter.\n", - "\n", - " **Parameters**\n", - "\n", - " z : numpy.array of shape (dim_z)\n", - " measurement vector\n", - "\n", - " R : numpy.array((dim_z, dim_z)), optional\n", - " Measurement noise. If provided, overrides self.R for\n", - " this function call.\n", - "\n", - " residual : function (z, z2), optional\n", - " Optional function that computes the residual (difference) between\n", - " the two measurement vectors. If you do not provide this, then the\n", - " built in minus operator will be used. You will normally want to use\n", - " the built in unless your residual computation is nonlinear (for\n", - " example, if they are angles)\n", - "\n", - " UT : function(sigmas, Wm, Wc, noise_cov), optional\n", - " Optional function to compute the unscented transform for the sigma\n", - " points passed through hx. Typically the default function will\n", - " work, but if for example you are using angles the default method\n", - " of computing means and residuals will not work, and you will have\n", - " to define how to compute it.\n", - " \"\"\"\n", - "\n", - " if isscalar(z):\n", - " dim_z = 1\n", - " else:\n", - " dim_z = len(z)\n", - "\n", - " if R is None:\n", - " R = self.R\n", - " elif np.isscalar(R):\n", - " R = eye(self._dim_z) * R\n", - "\n", - " # rename for readability\n", - " sigmas_f = self.sigmas_f\n", - " sigmas_h = zeros((self._num_sigmas, dim_z))\n", - "\n", - " if UT is None:\n", - " UT = unscented_transform\n", - "\n", - " # transform sigma points into measurement space\n", - " for i in range(self._num_sigmas):\n", - " sigmas_h[i] = self.hx(sigmas_f[i])\n", - "\n", - " # mean and covariance of prediction passed through inscented transform\n", - " zp, Pz = UT(sigmas_h, self.W, self.W, R)\n", - "\n", - " # compute cross variance of the state and the measurements\n", - " '''self.Pxz = zeros((self._dim_x, dim_z))\n", - " for i in range(self._num_sigmas):\n", - " self.Pxz += self.W[i] * np.outer(sigmas_f[i] - self.x,\n", - " residual(sigmas_h[i], zp))'''\n", - "\n", - " # this is the unreadable but fast implementation of the\n", - " # commented out loop above\n", - " yh = sigmas_f - self.x[np.newaxis, :]\n", - " yz = residual(sigmas_h, zp[np.newaxis, :])\n", - " self.Pxz = yh.T.dot(np.diag(self.W)).dot(yz)\n", - "\n", - " K = dot(self.Pxz, inv(Pz)) # Kalman gain\n", - " y = residual(z, zp)\n", - "\n", - " self.x = self.x + dot(K, y)\n", - " self.P = self.P - dot3(K, Pz, K.T)\n", - "\n", - "\n", - "\n", - " def predict(self, dt=None):\n", - " \"\"\" Performs the predict step of the UKF. On return, self.xp and\n", - " self.Pp contain the predicted state (xp) and covariance (Pp). 'p'\n", - " stands for prediction.\n", - "\n", - " **Parameters**\n", - " dt : double, optional\n", - " If specified, the time step to be used for this prediction.\n", - " self._dt is used if this is not provided.\n", - "\n", - " Important: this MUST be called before update() is called for the\n", - " first time.\n", - " \"\"\"\n", - "\n", - " if dt is None:\n", - " dt = self._dt\n", - "\n", - " # calculate sigma points for given mean and covariance\n", - " sigmas = self.sigma_points(self.x, self.P, self.kappa)\n", - "\n", - " for i in range(self._num_sigmas):\n", - " self.sigmas_f[i] = self.fx(sigmas[i], dt)\n", - "\n", - " self.x, self.P = unscented_transform(\n", - " self.sigmas_f, self.W, self.W, self.Q)\n", - "\n", - "\n", - " def batch_filter(self, zs, Rs=None, residual=np.subtract, UT=None):\n", - " \"\"\" Performs the UKF filter over the list of measurement in `zs`.\n", - "\n", - "\n", - " **Parameters**\n", - "\n", - " zs : list-like\n", - " list of measurements at each time step `self._dt` Missing\n", - " measurements must be represented by 'None'.\n", - "\n", - " Rs : list-like, optional\n", - " optional list of values to use for the measurement error\n", - " covariance; a value of None in any position will cause the filter\n", - " to use `self.R` for that time step.\n", - "\n", - " residual : function (z, z2), optional\n", - " Optional function that computes the residual (difference) between\n", - " the two measurement vectors. If you do not provide this, then the\n", - " built in minus operator will be used. You will normally want to use\n", - " the built in unless your residual computation is nonlinear (for\n", - " example, if they are angles)\n", - "\n", - " UT : function(sigmas, Wm, Wc, noise_cov), optional\n", - " Optional function to compute the unscented transform for the sigma\n", - " points passed through hx. Typically the default function will\n", - " work, but if for example you are using angles the default method\n", - " of computing means and residuals will not work, and you will have\n", - " to define how to compute it.\n", - "\n", - " **Returns**\n", - "\n", - " means: np.array((n,dim_x,1))\n", - " array of the state for each time step after the update. Each entry\n", - " is an np.array. In other words `means[k,:]` is the state at step\n", - " `k`.\n", - "\n", - " covariance: np.array((n,dim_x,dim_x))\n", - " array of the covariances for each time step after the update.\n", - " In other words `covariance[k,:,:]` is the covariance at step `k`.\n", - " \n", - " \"\"\"\n", - "\n", - " try:\n", - " z = zs[0]\n", - " except:\n", - " assert not isscalar(zs), 'zs must be list-like'\n", - "\n", - " if self._dim_z == 1:\n", - " assert isscalar(z) or (z.ndim==1 and len(z) == 1), \\\n", - " 'zs must be a list of scalars or 1D, 1 element arrays'\n", - "\n", - " else:\n", - " assert len(z) == self._dim_z, 'each element in zs must be a' \\\n", - " '1D array of length {}'.format(self._dim_z)\n", - "\n", - " n = np.size(zs,0)\n", - " if Rs is None:\n", - " Rs = [None]*n\n", - "\n", - " # mean estimates from Kalman Filter\n", - " if self.x.ndim == 1:\n", - " means = zeros((n, self._dim_x))\n", - " else:\n", - " means = zeros((n, self._dim_x, 1))\n", - "\n", - " # state covariances from Kalman Filter\n", - " covariances = zeros((n, self._dim_x, self._dim_x))\n", - " \n", - " for i, (z, r) in enumerate(zip(zs, Rs)):\n", - " self.predict()\n", - " self.update(z, r)\n", - " means[i,:] = self.x\n", - " covariances[i,:,:] = self.P\n", - " \n", - " return (means, covariances)\n", - "\n", - " \n", - "\n", - " def rts_smoother(self, Xs, Ps, Qs=None, dt=None):\n", - " \"\"\" Runs the Rauch-Tung-Striebal Kalman smoother on a set of\n", - " means and covariances computed by the UKF. The usual input\n", - " would come from the output of `batch_filter()`.\n", - "\n", - " **Parameters**\n", - "\n", - " Xs : numpy.array\n", - " array of the means (state variable x) of the output of a Kalman\n", - " filter.\n", - "\n", - " Ps : numpy.array\n", - " array of the covariances of the output of a kalman filter.\n", - "\n", - " Q : list-like collection of numpy.array, optional\n", - " Process noise of the Kalman filter at each time step. Optional,\n", - " if not provided the filter's self.Q will be used\n", - "\n", - " dt : optional, float or array-like of float\n", - " If provided, specifies the time step of each step of the filter.\n", - " If float, then the same time step is used for all steps. If\n", - " an array, then each element k contains the time at step k.\n", - " Units are seconds.\n", - "\n", - " **Returns**\n", - "\n", - " 'x' : numpy.ndarray\n", - " smoothed means\n", - "\n", - " 'P' : numpy.ndarray\n", - " smoothed state covariances\n", - "\n", - " 'K' : numpy.ndarray\n", - " smoother gain at each step\n", - "\n", - "\n", - " **Example**::\n", - "\n", - " zs = [t + random.randn()*4 for t in range (40)]\n", - "\n", - " (mu, cov, _, _) = kalman.batch_filter(zs)\n", - " (x, P, K) = rts_smoother(mu, cov, fk.F, fk.Q)\n", - "\n", - " \"\"\"\n", - " assert len(Xs) == len(Ps)\n", - " n, dim_x = Xs.shape\n", - "\n", - " if dt is None:\n", - " dt = [self._dt] * n\n", - " elif isscalar(dt):\n", - " dt = [dt] * n\n", - "\n", - " if Qs is None:\n", - " Qs = [self.Q] * n\n", - "\n", - " # smoother gain\n", - " Ks = zeros((n,dim_x,dim_x))\n", - "\n", - " num_sigmas = 2*dim_x + 1\n", - "\n", - " xs, ps = Xs.copy(), Ps.copy()\n", - " sigmas_f = zeros((num_sigmas, dim_x))\n", - "\n", - "\n", - " for k in range(n-2,-1,-1):\n", - " # create sigma points from state estimate, pass through state func\n", - " sigmas = self.sigma_points(xs[k], ps[k], self.kappa)\n", - " for i in range(num_sigmas):\n", - " sigmas_f[i] = self.fx(sigmas[i], dt[k])\n", - "\n", - " # compute backwards prior state and covariance\n", - " xb = dot(self.W, sigmas_f)\n", - " Pb = 0\n", - " x = Xs[k]\n", - " for i in range(num_sigmas):\n", - " y = sigmas_f[i] - x\n", - " Pb += self.W[i] * outer(y, y)\n", - " Pb += Qs[k]\n", - "\n", - " # compute cross variance\n", - " Pxb = 0\n", - " for i in range(num_sigmas):\n", - " z = sigmas[i] - Xs[k]\n", - " y = sigmas_f[i] - xb\n", - " Pxb += self.W[i] * outer(z, y)\n", - "\n", - " # compute gain\n", - " K = dot(Pxb, inv(Pb))\n", - "\n", - " # update the smoothed estimates\n", - " xs[k] += dot (K, xs[k+1] - xb)\n", - " ps[k] += dot3(K, ps[k+1] - Pb, K.T)\n", - " Ks[k] = K\n", - "\n", - " return (xs, ps, Ks)\n", - "\n", - "\n", - " @staticmethod\n", - " def weights(n, kappa):\n", - " \"\"\" Computes the weights for an unscented Kalman filter. See\n", - " __init__() for meaning of parameters.\n", - " \"\"\"\n", - "\n", - " assert kappa >= 0.0, \\\n", - " \"kappa cannot be negative, it's value is {}\".format(kappa)\n", - " assert n > 0, \"n must be greater than 0, it's value is {}\".format(n)\n", - "\n", - " k = .5 / (n+kappa)\n", - " W = np.full(2*n+1, k)\n", - " W[0] = kappa / (n+kappa)\n", - " return W\n", - "\n", - "\n", - " @staticmethod\n", - " def sigma_points(x, P, kappa):\n", - " \"\"\" Computes the sigma points for an unscented Kalman filter\n", - " given the mean (x) and covariance(P) of the filter.\n", - " kappa is an arbitrary constant. Returns sigma points.\n", - "\n", - " Works with both scalar and array inputs:\n", - " sigma_points (5, 9, 2) # mean 5, covariance 9\n", - " sigma_points ([5, 2], 9*eye(2), 2) # means 5 and 2, covariance 9I\n", - "\n", - " **Parameters**\n", - "\n", - " X An array-like object of the means of length n\n", - " Can be a scalar if 1D.\n", - " examples: 1, [1,2], np.array([1,2])\n", - "\n", - " P : scalar, or np.array\n", - " Covariance of the filter. If scalar, is treated as eye(n)*P.\n", - "\n", - " kappa : float\n", - " Scaling factor.\n", - "\n", - " **Returns**\n", - "\n", - " sigmas : np.array, of size (n, 2n+1)\n", - " 2D array of sigma points. Each column contains all of\n", - " the sigmas for one dimension in the problem space. They\n", - " are ordered as:\n", - "\n", - " .. math::\n", - " sigmas[0] = x \\n\n", - " sigmas[1..n] = x + [\\sqrt{(n+\\kappa)P}]_k \\n\n", - " sigmas[n+1..2n] = x - [\\sqrt{(n+\\kappa)P}]_k\n", - " \"\"\"\n", - "\n", - " if np.isscalar(x):\n", - " x = asarray([x])\n", - " n = np.size(x) # dimension of problem\n", - "\n", - " if np.isscalar(P):\n", - " P = eye(n)*P\n", - "\n", - " sigmas = zeros((2*n+1, n))\n", - "\n", - " # implements U'*U = (n+kappa)*P. Returns lower triangular matrix.\n", - " # Take transpose so we can access with U[i]\n", - " U = cholesky((n+kappa)*P).T\n", - " #U = sqrtm((n+kappa)*P).T\n", - "\n", - " sigmas[0] = x\n", - " sigmas[1:n+1] = x + U\n", - " sigmas[n+1:2*n+2] = x - U\n", - "\n", - " return sigmas\n", - "\n", - "\n", - "def unscented_transform(Sigmas, Wm, Wc, noise_cov):\n", - " \"\"\" Computes unscented transform of a set of sigma points and weights.\n", - " returns the mean and covariance in a tuple.\n", - " \"\"\"\n", - "\n", - " kmax, n = Sigmas.shape\n", - "\n", - " # new mean is just the sum of the sigmas * weight\n", - " x = dot(Wm, Sigmas) # dot = \\Sigma^n_1 (W[k]*Xi[k])\n", - "\n", - " # new covariance is the sum of the outer product of the residuals\n", - " # times the weights\n", - " '''P = zeros((n, n))\n", - " for k in range(kmax):\n", - " y = Sigmas[k] - x\n", - " P += Wc[k] * np.outer(y, y)'''\n", - "\n", - " # this is the fast way to do the commented out code above\n", - " y = Sigmas - x[np.newaxis,:]\n", - " P = y.T.dot(np.diag(Wc)).dot(y)\n", - "\n", - " if noise_cov is not None:\n", - " P += noise_cov\n", - "\n", - " return (x, P)" - ] - }, { "cell_type": "markdown", "metadata": {}, diff --git a/Appendix_B_Symbols_and_Notations.ipynb b/Appendix_B_Symbols_and_Notations.ipynb index c049a21..61a27df 100644 --- a/Appendix_B_Symbols_and_Notations.ipynb +++ b/Appendix_B_Symbols_and_Notations.ipynb @@ -1,427 +1,427 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/html": [ - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "#format the book\n", - "%matplotlib inline\n", - "from __future__ import division, print_function\n", - "import matplotlib.pyplot as plt\n", - "import book_format\n", - "book_format.load_style()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Symbology" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "This is just notes at this point. \n", - "\n", - "\n", - "## State\n", - "\n", - "$x$ (Brookner, Zarchan, Brown)\n", - "\n", - "$\\underline{x}$ Gelb)\n", - "\n", - "## State at step n\n", - "\n", - "$x_n$ (Brookner)\n", - "\n", - "$x_k$ (Brown, Zarchan)\n", - "\n", - "$\\underline{x}_k$ (Gelb)\n", - "\n", - "\n", - "\n", - "## Prediction\n", - "\n", - "$x^-$\n", - "\n", - "$x_{n,n-1}$ (Brookner) \n", - "\n", - "$x_{k+1,k}$\n", - "\n", - "\n", - "## measurement\n", - "\n", - "\n", - "$x^*$\n", - "\n", - "\n", - "\n", - "Y_n (Brookner)\n", - "\n", - "##control transition Matrix\n", - "\n", - "$G$ (Zarchan)\n", - "\n", - "\n", - "Not used (Brookner)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#Nomenclature\n", - "\n", - "\n", - "## Equations\n", - "### Brookner\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "X^*_{n+1,n} &= \\Phi X^*_{n,n} \\\\\n", - "X^*_{n,n} &= X^*_{n,n-1} +H_n(Y_n - MX^*_{n,n-1}) \\\\\n", - "H_n &= S^*_{n,n-1}M^T[R_n + MS^*_{n,n-1}M^T]^{-1} \\\\\n", - "S^*_{n,n-1} &= \\Phi S^*_{n-1,n-1}\\Phi^T + Q_n \\\\\n", - "S^*_{n-1,n-1} &= (I-H_{n-1}M)S^*_{n-1,n-2}\n", - "\\end{aligned}$$\n", - "\n", - "### Gelb\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "\\underline{\\hat{x}}_k(-) &= \\Phi_{k-1} \\underline{\\hat{x}}_{k-1}(+) \\\\\n", - "\\underline{\\hat{x}}_k(+) &= \\underline{\\hat{x}}_k(-) +K_k[Z_k - H_k\\underline{\\hat{x}}_k(-)] \\\\\n", - "K_k &= P_k(-)H_k^T[H_kP_k(-)H_k^T + R_k]^{-1}\\\\\n", - "P_k(+) &= \\Phi_{k-1} P_{k-1}(+)\\Phi_{k-1}^T + Q_{k-1} \\\\\n", - "P_k(-) &= (I-K_kH_k)P_k(-)\n", - "\\end{aligned}$$\n", - "\n", - "\n", - "### Brown\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "\\hat{\\textbf{x}}^-_{k+1} &= \\mathbf{\\phi}_{k}\\hat{\\textbf{x}}_{k} \\\\\n", - "\\hat{\\textbf{x}}_k &= \\hat{\\textbf{x}}^-_k +\\textbf{K}_k[\\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{}x}^-_k] \\\\\n", - "\\textbf{K}_k &= \\textbf{P}^-_k\\textbf{H}_k^T[\\textbf{H}_k\\textbf{P}^-_k\\textbf{H}_k^T + \\textbf{R}_k]^{-1}\\\\\n", - "\\textbf{P}^-_{k+1} &= \\mathbf{\\phi}_k \\textbf{P}_k\\mathbf{\\phi}_k^T + \\textbf{Q}_{k} \\\\\n", - "\\mathbf{P}_k &= (\\mathbf{I}-\\mathbf{K}_k\\mathbf{H}_k)\\mathbf{P}^-_k\n", - "\\end{aligned}$$\n", - "\n", - "\n", - "### Zarchan\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "\\hat{x}_{k} &= \\Phi_{k}\\hat{x}_{k-1} + G_ku_{k-1} + K_k[z_k - H\\Phi_{k}\\hat{x}_{k-1} - HG_ku_{k-1} ] \\\\\n", - "M_{k} &= \\Phi_k P_{k-1}\\phi_k^T + Q_{k} \\\\\n", - "K_k &= M_kH^T[HM_kH^T + R_k]^{-1}\\\\\n", - "P_k &= (I-K_kH)M_k\n", - "\\end{aligned}$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Wikipedia\n", - "$$\n", - "\\begin{aligned}\n", - "\\hat{\\textbf{x}}_{k\\mid k-1} &= \\textbf{F}_{k}\\hat{\\textbf{x}}_{k-1\\mid k-1} + \\textbf{B}_{k} \\textbf{u}_{k} \\\\\n", - "\\textbf{P}_{k\\mid k-1} &= \\textbf{F}_{k} \\textbf{P}_{k-1\\mid k-1} \\textbf{F}_{k}^{\\text{T}} + \\textbf{Q}_{k}\\\\\n", - "\\tilde{\\textbf{y}}_k &= \\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{x}}_{k\\mid k-1} \\\\\n", - "\\textbf{S}_k &= \\textbf{H}_k \\textbf{P}_{k\\mid k-1} \\textbf{H}_k^\\text{T} + \\textbf{R}_k \\\\\n", - "\\textbf{K}_k &= \\textbf{P}_{k\\mid k-1}\\textbf{H}_k^\\text{T}\\textbf{S}_k^{-1} \\\\\n", - "\\hat{\\textbf{x}}_{k\\mid k} &= \\hat{\\textbf{x}}_{k\\mid k-1} + \\textbf{K}_k\\tilde{\\textbf{y}}_k \\\\\n", - "\\textbf{P}_{k|k} &= (I - \\textbf{K}_k \\textbf{H}_k) \\textbf{P}_{k|k-1}\n", - "\\end{aligned}$$" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Labbe\n", - "\n", - "$$\n", - "\\begin{aligned}\n", - "\\hat{\\textbf{x}}^-_{k+1} &= \\mathbf{F}_{k}\\hat{\\textbf{x}}_{k} + \\mathbf{B}_k\\mathbf{u}_k \\\\\n", - "\\textbf{P}^-_{k+1} &= \\mathbf{F}_k \\textbf{P}_k\\mathbf{F}_k^T + \\textbf{Q}_{k} \\\\\n", - "\\textbf{y}_k &= \\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{}x}^-_k \\\\\n", - "\\mathbf{S}_k &= \\textbf{H}_k\\textbf{P}^-_k\\textbf{H}_k^T + \\textbf{R}_k \\\\\n", - "\\textbf{K}_k &= \\textbf{P}^-_k\\textbf{H}_k^T\\mathbf{S}_k^{-1} \\\\\n", - "\\hat{\\textbf{x}}_k &= \\hat{\\textbf{x}}^-_k +\\textbf{K}_k\\textbf{y} \\\\\n", - "\\mathbf{P}_k &= (\\mathbf{I}-\\mathbf{K}_k\\mathbf{H}_k)\\mathbf{P}^-_k\n", - "\\end{aligned}$$" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.4.3" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#format the book\n", + "%matplotlib inline\n", + "from __future__ import division, print_function\n", + "import matplotlib.pyplot as plt\n", + "import book_format\n", + "book_format.load_style()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Symbology" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is just notes at this point. \n", + "\n", + "\n", + "## State\n", + "\n", + "$x$ (Brookner, Zarchan, Brown)\n", + "\n", + "$\\underline{x}$ Gelb)\n", + "\n", + "## State at step n\n", + "\n", + "$x_n$ (Brookner)\n", + "\n", + "$x_k$ (Brown, Zarchan)\n", + "\n", + "$\\underline{x}_k$ (Gelb)\n", + "\n", + "\n", + "\n", + "## Prediction\n", + "\n", + "$x^-$\n", + "\n", + "$x_{n,n-1}$ (Brookner) \n", + "\n", + "$x_{k+1,k}$\n", + "\n", + "\n", + "## measurement\n", + "\n", + "\n", + "$x^*$\n", + "\n", + "\n", + "\n", + "Y_n (Brookner)\n", + "\n", + "##control transition Matrix\n", + "\n", + "$G$ (Zarchan)\n", + "\n", + "\n", + "Not used (Brookner)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "##Nomenclature\n", + "\n", + "\n", + "### Equations\n", + "#### Brookner\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "X^*_{n+1,n} &= \\Phi X^*_{n,n} \\\\\n", + "X^*_{n,n} &= X^*_{n,n-1} +H_n(Y_n - MX^*_{n,n-1}) \\\\\n", + "H_n &= S^*_{n,n-1}M^T[R_n + MS^*_{n,n-1}M^T]^{-1} \\\\\n", + "S^*_{n,n-1} &= \\Phi S^*_{n-1,n-1}\\Phi^T + Q_n \\\\\n", + "S^*_{n-1,n-1} &= (I-H_{n-1}M)S^*_{n-1,n-2}\n", + "\\end{aligned}$$\n", + "\n", + "#### Gelb\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\underline{\\hat{x}}_k(-) &= \\Phi_{k-1} \\underline{\\hat{x}}_{k-1}(+) \\\\\n", + "\\underline{\\hat{x}}_k(+) &= \\underline{\\hat{x}}_k(-) +K_k[Z_k - H_k\\underline{\\hat{x}}_k(-)] \\\\\n", + "K_k &= P_k(-)H_k^T[H_kP_k(-)H_k^T + R_k]^{-1}\\\\\n", + "P_k(+) &= \\Phi_{k-1} P_{k-1}(+)\\Phi_{k-1}^T + Q_{k-1} \\\\\n", + "P_k(-) &= (I-K_kH_k)P_k(-)\n", + "\\end{aligned}$$\n", + "\n", + "\n", + "#### Brown\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\hat{\\textbf{x}}^-_{k+1} &= \\mathbf{\\phi}_{k}\\hat{\\textbf{x}}_{k} \\\\\n", + "\\hat{\\textbf{x}}_k &= \\hat{\\textbf{x}}^-_k +\\textbf{K}_k[\\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{}x}^-_k] \\\\\n", + "\\textbf{K}_k &= \\textbf{P}^-_k\\textbf{H}_k^T[\\textbf{H}_k\\textbf{P}^-_k\\textbf{H}_k^T + \\textbf{R}_k]^{-1}\\\\\n", + "\\textbf{P}^-_{k+1} &= \\mathbf{\\phi}_k \\textbf{P}_k\\mathbf{\\phi}_k^T + \\textbf{Q}_{k} \\\\\n", + "\\mathbf{P}_k &= (\\mathbf{I}-\\mathbf{K}_k\\mathbf{H}_k)\\mathbf{P}^-_k\n", + "\\end{aligned}$$\n", + "\n", + "\n", + "#### Zarchan\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\hat{x}_{k} &= \\Phi_{k}\\hat{x}_{k-1} + G_ku_{k-1} + K_k[z_k - H\\Phi_{k}\\hat{x}_{k-1} - HG_ku_{k-1} ] \\\\\n", + "M_{k} &= \\Phi_k P_{k-1}\\phi_k^T + Q_{k} \\\\\n", + "K_k &= M_kH^T[HM_kH^T + R_k]^{-1}\\\\\n", + "P_k &= (I-K_kH)M_k\n", + "\\end{aligned}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Wikipedia\n", + "$$\n", + "\\begin{aligned}\n", + "\\hat{\\textbf{x}}_{k\\mid k-1} &= \\textbf{F}_{k}\\hat{\\textbf{x}}_{k-1\\mid k-1} + \\textbf{B}_{k} \\textbf{u}_{k} \\\\\n", + "\\textbf{P}_{k\\mid k-1} &= \\textbf{F}_{k} \\textbf{P}_{k-1\\mid k-1} \\textbf{F}_{k}^{\\text{T}} + \\textbf{Q}_{k}\\\\\n", + "\\tilde{\\textbf{y}}_k &= \\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{x}}_{k\\mid k-1} \\\\\n", + "\\textbf{S}_k &= \\textbf{H}_k \\textbf{P}_{k\\mid k-1} \\textbf{H}_k^\\text{T} + \\textbf{R}_k \\\\\n", + "\\textbf{K}_k &= \\textbf{P}_{k\\mid k-1}\\textbf{H}_k^\\text{T}\\textbf{S}_k^{-1} \\\\\n", + "\\hat{\\textbf{x}}_{k\\mid k} &= \\hat{\\textbf{x}}_{k\\mid k-1} + \\textbf{K}_k\\tilde{\\textbf{y}}_k \\\\\n", + "\\textbf{P}_{k|k} &= (I - \\textbf{K}_k \\textbf{H}_k) \\textbf{P}_{k|k-1}\n", + "\\end{aligned}$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Labbe\n", + "\n", + "$$\n", + "\\begin{aligned}\n", + "\\hat{\\textbf{x}}^-_{k+1} &= \\mathbf{F}_{k}\\hat{\\textbf{x}}_{k} + \\mathbf{B}_k\\mathbf{u}_k \\\\\n", + "\\textbf{P}^-_{k+1} &= \\mathbf{F}_k \\textbf{P}_k\\mathbf{F}_k^T + \\textbf{Q}_{k} \\\\\n", + "\\textbf{y}_k &= \\textbf{z}_k - \\textbf{H}_k\\hat{\\textbf{}x}^-_k \\\\\n", + "\\mathbf{S}_k &= \\textbf{H}_k\\textbf{P}^-_k\\textbf{H}_k^T + \\textbf{R}_k \\\\\n", + "\\textbf{K}_k &= \\textbf{P}^-_k\\textbf{H}_k^T\\mathbf{S}_k^{-1} \\\\\n", + "\\hat{\\textbf{x}}_k &= \\hat{\\textbf{x}}^-_k +\\textbf{K}_k\\textbf{y} \\\\\n", + "\\mathbf{P}_k &= (\\mathbf{I}-\\mathbf{K}_k\\mathbf{H}_k)\\mathbf{P}^-_k\n", + "\\end{aligned}$$" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.3" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/Appendix_C_Walking_Through_KF_Code.ipynb b/Appendix_C_Walking_Through_KF_Code.ipynb index 71b8e0c..8e43ee1 100644 --- a/Appendix_C_Walking_Through_KF_Code.ipynb +++ b/Appendix_C_Walking_Through_KF_Code.ipynb @@ -11,7 +11,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Walking through the Kalman Filter code" + "# Walking through the Kalman Filter code" ] }, { diff --git a/16_HInfinity_Filters.ipynb b/Appendix_D_HInfinity_Filters.ipynb similarity index 99% rename from 16_HInfinity_Filters.ipynb rename to Appendix_D_HInfinity_Filters.ipynb index 1ff0f0a..c0d399b 100644 --- a/16_HInfinity_Filters.ipynb +++ b/Appendix_D_HInfinity_Filters.ipynb @@ -1,360 +1,360 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# H Infinity filter" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "text/html": [ - "\n", - "\n" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 1, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "#format the book\n", - "%matplotlib inline\n", - "from __future__ import division, print_function\n", - "import book_format\n", - "book_format.load_style()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "I am still mulling over how to write this chapter. In the meantime, Professor Dan Simon at Clevant State University has an accessible introduction here:\n", - "\n", - "http://academic.csuohio.edu/simond/courses/eec641/hinfinity.pdf\n", - "\n", - "In one sentence the $H_\\infty$ (H infinity) filter is like a Kalman filter, but it is robust in the face of non-Gaussian, non-predictable inputs.\n", - "\n", - "\n", - "My FilterPy library contains an H-Infinity filter. I've pasted some test code below which implements the filter designed by Simon in the article above. Hope it helps." - ] - }, - { - "cell_type": "code", - "execution_count": 2, - "metadata": { - "collapsed": false - }, - "outputs": [ - { - "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAEjCAYAAABzSrSOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNXdP/DPnX0mmZlkshMCBBuIIJuEVVlEoWpbHq21\nCiqLWGqLVuHxUcGFgAiVUmoXsBZbjFLc2ufRujwUrQqkwPMjAiqraDAoYYZsM1lnMsv5/TGTSyb7\nMjfr5/16zevOnHvunDPndYWvh+89RxJCCBARERERUZdQdXcHiIiIiIj6EwbgRERERERdiAE4ERER\nEVEXYgBORERERNSFGIATEREREXUhBuBERERERF2IATgRUR/y4osvQqVSIScnp13XqVQqXHPNNQr1\nioiI6mMATkTUh0iSJL/qU6lUSE9Pb/VaIiJSnqa7O0BERJFz8803Y8qUKUhOTm50rqUA+9SpUzCZ\nTEp2jYiIQhiAExH1IRaLBRaLpd3XDRs2TIHeEBFRU5iCQkTUSV9//bWcQ33+/HnccccdSEhIgMlk\nwoQJE/D66683ukYIgRdeeAGTJ0+G2WxGVFQUrrzySmzevBk+n69R/c8++wzz589Heno6jEYjEhIS\nMGbMGCxbtgzl5eVyvYY54B9//DFUKlVYP+teixcvlq9rLge8oqICjz/+ODIzM2E0GhEbG4trr70W\n//jHP1och5KSEixduhQpKSkwGAy44oor8OKLL7Z7bImI+iLOgBMRRUhZWRmuvvpq2Gw23HPPPSgt\nLcXrr7+O22+/HYWFhXjwwQflugsXLsSOHTswcOBALFmyBFqtFv/4xz/w0EMPYffu3Xj33XehVqsB\nBIPvSZMmQa1W4/vf/z4uu+wyVFZWIj8/Hzk5OfjP//zPRrPedekm6enpWL16NdasWQOr1Yrly5fL\ndcaOHdvkNXVcLheuvvpqHD9+HFdeeSUefPBBlJWV4Y033sBNN92ENWvW4Iknnmg0Dk6nE1dddRX0\nej1+/OMfw+Px4PXXX8fdd98NlUqFBQsWdG6giYh6O0FERJ1y9uxZIUmSkCRJ3H777WHnvvrqKxET\nEyP0er0oKCgQQgjx6quvCkmSxNixY0VFRYVct7a2Vlx77bVCkiSxadMmuXzFihVCkiTx1ltvNWq7\nsrJSeDwe+fP27duFJEkiJycnrJ4kSSI9Pb3Z3yBJkrjmmmvCyu69914hSZJYsmRJWPm3334rUlJS\nhEqlEocOHWpyHH7yk5+IQCAgnztx4oTQaDRixIgRzfaBiKi/YAoKEVGEaDQa/PKXvwwrGzp0KH7+\n85+jtrYWO3bsAAC88MILAIANGzYgOjparqvVavGb3/wGALBt27ZG328wGBqVRUVFQafTRew31PF6\nvXjppZcQFRWFjRs3hp1LTU3FqlWr5DSapvq0efPmsBn1yy+/HFOnTsWpU6dQXV0d8f4SEfUmDMCJ\niCJk0KBBGDx4cKPy6dOnAwCOHj0KADh8+DAkSWoy53rUqFFISEjAmTNn5ED19ttvh1qtxk033YS7\n7roL27dvx6lTpxT8JcFVUWpqajBq1CjYbLZG56+77joAwJEjRxqdy8jICPsfizppaWkQQqCsrCzy\nHSYi6kUYgBMRRUhSUlKL5S6XSz5arVbo9fom66ekpITVnzBhAnJzczF79mz8z//8D5YsWYIRI0Zg\n6NCheP755yP9M8Labmo5w/rlTqez0bmYmJgmr9Fogo8d+f3+SHSRiKjXYgBORBQhDoejxXKr1Sof\nXS4XPB5Pk/UvXLgQVh8AJk2ahLfeegtOpxMHDx7EU089hZqaGvzsZz/Dyy+/HMmfEda23W5vcx+J\niKhtGIATEUXIuXPnUFBQ0Kh8z549AIBx48YBAMaPHw8hBD766KNGdY8dO4aioiIMGzasyY1xNBoN\nJk6ciMceewwvvfQSAODNN99stW+SJLVr5vnyyy+H0WjE559/jpKSkkbn//Wvf8m/hYiI2ocBOBFR\nhPh8PjzyyCMQQshlX331FbZu3QqdToc77rgDALBkyRIAwKpVq1BVVSXX9Xq9WLFiBQDgnnvukcv3\n798Pt9vdqL26Wei27GAZFxeHoqKiJr+nKRqNBgsWLEBVVRVWrlwZdq6wsBAbNmyASqXC3Xff3abv\nIyKiS7gOOBFRhIwePRr/93//h6ysLMyePVteB7yiogKbN2/GoEGDAAC33XYb3n77bezcuRMjRozA\nTTfdBK1Wi7fffhtnzpzBddddF7Ze98aNG/Hhhx9i2rRpGDJkCCwWC7744gu88847MJlMYeuLN2fO\nnDnYuXMnrr/+ekybNg16vR5jx47F97///Wav+eUvf4l9+/bhhRdewJEjR3DttdfC6XTijTfegNPp\nxJNPPokJEyZ0fuCIiPqZds+Ab9iwARMmTIDVakViYiLmzp2L48ePN6qXnZ2N1NRUmEwmXHPNNThx\n4kREOkxE1FPZbDb8+9//xuWXX46//OUv2LFjBzIyMvDKK6/ggQceCKv78ssv449//COSk5Px5z//\nGVu3boXJZMKvfvUrvPfee/LulQCwbNky3HrrrSgoKMDOnTvxu9/9DseOHcOiRYtw+PDhsDQQSZIa\nbagDAM8++yzuuusufPHFF9iwYQNWr16N//7v/27x91itVuzfvx8rV65ERUUFnn32Wbz66qsYM2YM\n/v73v2P16tXtGp/m+kZE1N9Iov6/lbbB9ddfj3nz5mHChAkIBAJ48sknceDAAZw4cQKxsbEAgGee\neQZPP/00cnJyMGzYMKxduxa5ubk4ffp0k0tTERH1Zl9//TWGDh2KmTNn4sMPP+zu7hARUQ/X7gC8\noaqqKlitVrz11lv43ve+ByEEBgwYgF/84hdy3qDb7UZiYiI2bdqEpUuXRqTjREQ9BQNwIiJqj04/\nhFleXo5AICDPfp89exYOhwNz5syR6xgMBkyfPh379+/vbHNERERERL1apwPwBx54AOPGjcOUKVMA\nXFoztuGGFImJic2uJ0tERERE1F90ahWUFStWYP/+/cjNzW3TgzUN69TttEZE1JvFxsbK26vzzzUi\nov6lIxuSdXgGfPny5Xjttdfw4YcfYsiQIXJ53fbEDXeEczgczW5pTERERETUX3QoAH/ggQfk4HvY\nsGFh59LT05GcnIzdu3fLZW63G7m5uZg6dWrnektERERE1Mu1OwVl2bJl2LFjB958801YrVY5r9ts\nNiMqKgqSJOHBBx/E+vXrkZmZiYyMDKxbtw5msxnz589v9ns7Mn1PTcvLywMAZGVldXNP+haOqzI4\nrsrguCqD46ocjq0yOK7K6Gy6YbsD8Oeeew6SJOHaa68NK8/OzsaTTz4JAHj44YdRU1ODZcuWoays\nDJMnT8bu3bsRFRXVqc4SERERUd8XEAE4K0pQ7LLLrxKXHd+bcgcSYwd0d/c6rd0BeCAQaFO91atX\nt3uXNCIiIiLqH2p9HpS4LqLYdQElLkd4sF3ugN/va3TN+OHT+mcATkRERETUGiEEqt0VoaA6PMAu\ndtnhqixp93cWu/rGktYMwImIiIioQwIBP5yV9VNFHCh2XQjOYjvtqKmtjmh7xU4G4ERERETUxwVT\nRUIz2M5L+djFLjtKKi42mSoSCVEGM+KtycFXTPCYlniZIm11NQbgRERERP2YEAJVdakizgsNHnx0\nwFVVqki7kqRCbHScHGDHWVNCAXcS4q3JMOr77uIdDMCJiIiI+rhAwI+yymJ5BrsuuK57745wqkgd\nrUaHeGsy4qzJiLckyTPZ8dZk2CyJ0Ki1irTb0zEAJyIiIuoDar2eRsH1V+dOo8Jdhr8eKIc/oFCq\niNFyKVUkNHsdfKXAEhULSZIUabc3YwBORERE1AsIIVBZU46SckejVJFilx3lVWWKtNt8qkgw4O7L\nqSJKYQBORERE1EM0lSpSf1a7S1JFGrxsloR+myqiFAbgRERERF2ofqpIw1zs0vKLXZQqkhy2wojF\nxFSRrsQAnIiIiCiCWlpVRPFUEXN8WHDtKq5CtCEW06fMglFvUqRdaj8G4ERERETt1N2ritS94lpJ\nFcnLywMABt89DANwIiIioiY0tapIV6SKRButiAtbTYSpIn0NA3AiIiLql3pSqoi8woglmbPV/UCH\nAvC9e/di06ZNOHz4MAoLC7F9+3YsXLhQPr9o0SK89NJLYddMnjwZ+/fv71xviYiIiNqhu1cVCQ+w\ng8v3xZrjuapIP9ehALyqqgqjR4/GwoULsWDBgkb/FCJJEmbPno2XX35ZLtPpdJ3rKREREVETGq4q\nUj9lROlUkWAedhJTRahdOhSA33DDDbjhhhsABGe7GxJCQKfTITExsVOdIyIiImKqCPU1iuSAS5KE\n3NxcJCUlISYmBjNmzMDTTz+NhIQEJZojIiKiXi4gAqj2lOP0uU+bnMnmBjTUlygSgF9//fW45ZZb\nkJ6ejrNnz+Lxxx/HrFmz8MknnzAVhYiIqJ+q9XkurSbibJgyYkdABIBPIt8uN6ChnkYSQojOfIHZ\nbMaWLVuwYMGCZutcuHABgwcPxmuvvYabb75ZLne5XPL7M2fOdKYbRERE1M2EEPD4alDpLkNFo5cT\nNbUVirQrQUKU3oJoQyzMDV7RhljoNHpF2qX+KyMjQ35vtVrbfX2XLEOYkpKCgQMH4ssvv+yK5oiI\niEghdaki9QPr+gG31+9RpF21SiMH1A2D7Ci9FWqVWpF2iZTQJQF4UVERzp8/j5SUlGbrZGVldUVX\n+oW6Xa84ppHFcVUGx1UZHFdl9JdxDaaKXESx64KcHlKXMlJScRF+vzKriug1JiTHD2SqSAT1l3u2\nq9XP4uiIDi9DWJcyEggEUFBQgKNHjyIuLg42mw2rV6/Gj370IyQnJ+Prr7/GypUrkZSUFJZ+QkRE\nRN1DCIHqulVF6l5OO4rLg/nZrsoSRdqVJBVio+MurSRiTZGD7G/yC6HTGBgoUr/QoQD80KFDmDVr\nFoDgiierV6/G6tWrsWjRImzduhXHjh3Dyy+/DKfTiZSUFMyaNQt/+9vfEBUVFdHOExERUdMCIgBn\nRUmDJftCM9pOO2q6fFWRJNgsic2uKuI4V6pIf4h6og4F4DNnzkQgEGj2/K5duzrcISIiImobr68W\nJeWOpnd5LHcolioSZTCHpYdcCrJTYIliqghRa7okB5yIiIg6pspd0WSA3RWpInEN8rDr3hv1/Bdt\nos5gAE5ERNSNAiIAV2VpvVzs8J0eazxVirSrVevCt1CPqb8BTfOpIkTUeQzAiYiIFOb1eVEaesCx\nUaqIywGf36tIu6b6qSINAm1LVCxUkkqRdomoZQzAiYiIIqDaU9l0LrbTDmdlCQQ6te9dkyRIiDHH\ny8F1wxltkz464m0SUecxACciImqD+hvQHDzuahRoV7uV2eVRo9aGB9ZygJ0CmzkRWg1TRYh6Gwbg\nREREIT6/F6XlF8PXxq6XKuL11wYrHotsuyZ9dFh6SP2HH63RNqaKEPUxDMCJiKhfqfFUN8jBviAH\n2mWVJRCi+WV2O0qChJjoOMTFJCPekiTPYNcF2SYDU0WI+hMG4ERE1KcIIVBeXdbs0n1VNeWKtKtW\naxBnuZQqklAvwLZZEqHV6BRpl4h6HwbgRETU6/j9PpRWFDW5bF+Jy4Fan0eRdrVqPcyGWAwaMBTx\n9bZRj7cmIybaBpVKrUi7RNS3MAAnIqIeyVNbg2JXw6X7gsF2WXkRAgqkigCANcrWzC6PyThx7DQk\nSUJWVpYibRNR/8AAnIiIuoUQApU1riYfeCx22VFR7VSkXbVKA5slscm1seOsSdBp9M1eyy3WiSgS\nGIATEZFi/AE/nBXF4TPY9QJtj9etSLt6nbHxsn2hQDs2Op6pIkTUrRiAExFRp9R6PU3v8Oi0o6Ti\nIgIBvyLtWkyxclAd1yDQjjZaOFtNRD1WhwLwvXv3YtOmTTh8+DAKCwuxfft2LFy4MKxOdnY2tm3b\nhrKyMkyaNAlbtmzBiBEjItJpIiLqOkIIVLsrwgNspx1FoXzs8qoyRdpVqdSwmROazMeOsyZDrzUo\n0i4RkdI6FIBXVVVh9OjRWLhwIRYsWNBoluGZZ57B5s2bkZOTg2HDhmHt2rWYPXs2Tp8+jehornVK\nRNTTBEQAzoqSJh94LHHaUVNbrUi7Oq2h+VQRcwLUTBUhoj6oQwH4DTfcgBtuuAEAsGjRorBzQgg8\n++yzWLlyJW6++WYAQE5ODhITE7Fz504sXbq0cz0mIqIO8fq8KC0PX1WkKLSEX0m5A36/T5F2o43W\nJh54DC7hZzZZmSpCRP1OxHPAz549C4fDgTlz5shlBoMB06dPx/79+xmAExEpqNpTiWKnHV8Xn0CF\nuwxfOA/Ks9jOyhIIiIi3KUkqxJrjm5zFjrMkw6g3RbxNIqLeLOIBuN1uBwAkJSWFlScmJqKwsDDS\nzRER9StCCJRXlcnpIQ2X76tyVyjSrlatQ5w1qV5wXX+XxwRo1FpF2iUi6ou6dBWUlv6ZMS8vrwt7\n0j9wTJXBcVUGx/WSQMCPSo8LFe4yVLjLUBk61r38AWVSRXQaI8yGWJgNMaHjpZdRZw7/M9wLVBcH\ncK64EOfQ/yZXeL8qh2OrDI5rZGVkZHTq+ogH4MnJyQAAh8OBgQMHyuUOh0M+R0TU33n9tY0C64qa\n4LHK41IkVQQATDrzpcDaaAsLsnUaripCRNQVIh6Ap6enIzk5Gbt378b48eMBAG63G7m5udi0aVOz\n13Fb38ip+79cjmlkcVyV0VfHNbjLY3mDpfsupY0otsujWoM4SxK0CM5mj8gYXS8fOwlajU6RdvuL\nvnq/9gQcW2VwXJXhcrk6dX2HlyE8c+YMACAQCKCgoABHjx5FXFwc0tLS8OCDD2L9+vXIzMxERkYG\n1q1bB7PZjPnz53eqs0REPUkg4IezsqTRiiLyLo+1NYq0a9CZmtxGPd6ajJjoOKhU6kt/6Y7jX7pE\nRD1NhwLwQ4cOYdasWQCCed2rV6/G6tWrsWjRIvzlL3/Bww8/jJqaGixbtgxlZWWYPHkydu/ejaio\nqIh2nohIaV5fLUrKHZeW66v30GNJxUXFlu5rapfHhNCDj1EGM5fuIyLqxToUgM+cOROBQKDFOnVB\nORFRT1e3dF/DNJFilx2uylJF8rFVkgo2SyJ3eSQi6oe6dBUUIqLu0HDpviKnPSzIrlZo6T6dRh8W\nXMfVSxuxmROgVvOPYCKi/oh/+hNRn+D3+1BaURQ2i11UlzLissPrq1Wk3SijRQ6qE6wpYTPZZlMM\nU0WIiKgRBuBE1Gt4vO6wTWeCr2CwXVZehIBoOTWuIyRIiGm4y2O9TWi4yyMREbUXA3Ai6jGEEKhy\nV4StKFL/ocfy6jJF2lWrNYi3NM7FDu7ymASthrs8EhFR5DAAJ6IuFRABOCsuLd137OsjqHCX4cMv\ndqLYZYe7tlqRdo06E+Lk4DolLMiOibZBpVIr0i4REVFDDMCJKOK8Pi9Kyx0NNqGxo8h1ASXlDuWW\n7ouKDU8VqZcuwqX7iIiop2AATkQdUuOpvpSD3SAv21lRzKX7iIiImsEAnIiaJIRARbWz0Qx23eeq\nmnJF2tVqdJdWFYlJ4dJ9RETU5/BvMqJ+zB/wo6yiqPHKIs4LKC53oNbrVqTdKINZDqprqwXMhlhM\nGDsF8THJsJhimSpCRER9GgNwoj6u1udBicvRaBa7xBncSj0Q8Ee8TQkSYqLj6j30GL6VulEfJdfN\ny8sDAFyWOiLi/SAiIuqJGIAT9QHV7sqwpfvqz2S7qkoVaVOt0iCuLh+73rrY8THJiLMkQavRKdIu\nERFRb8cAnKgXqNtKvajhA4+hgLvaU6lIu3qtITy4rjeLHRMdx6X7iIiIOoABOFEP0XAr9aJQkK30\nVupmo1VOFWm4lXq00cp8bCIioghTJADPzs7G2rVrw8qSk5NRWFioRHNEvYbH65YD6iJnF22lLqkQ\nGx3XaCY7GGinwKAzRrxNIiIiap5iM+CZmZn4+OOP5c9qNf+pmvo+IQSq3RXB2etG+djKbaWuUWsR\nZ01qlCYS3Eo9ERo1t1InIiLqKRQLwNVqNRITE5X6eqJuExABuCpL5KC62BW+skiNglupN57FDn62\nRtugklSKtEtERESRpVgAnp+fj9TUVOj1ekyaNAnr169Henq6Us0RRZTP70Vp+cWwlUW+KjiNCncZ\ndh50wef3KtJu/a3UExoE2yZupU5ERNQnSEKIiO8XvWvXLlRWViIzMxMOhwPr1q3DqVOncPz4cdhs\nNrmey+WS3585cybS3SBqkdfnQYW7rMlXtadcka3UJUiIMlhhNsTCbLCFjpdeTBUhIiLq+TIyMuT3\nVqu13dcrEoA3VF1djfT0dDz66KNYvny5XM4AnJQkhIDbW40Kd6kcWFe6y1DhdqLCXQq3V5lUEbVK\nExZURxtiYQm9j9JbuXQfERFRL9fZALxLliE0mUwYOXIkvvzyy2brZGVldUVX+oW6nQX7w5gGAn6U\nVRbXWxs7fJ1sj0JbqZvqbaUeTBepe58CSxS3Um+P/nS/diWOqzI4rsrh2CqD46qM+pPIHdElAbjb\n7cbJkycxa9asrmiO+hivrxYl5Q55G/WSeg8/lpRfhD/gU6Rda3TcpQDbmgxXSTXMhlhMnzILJkO0\nIm0SERFR36dIAP7QQw9h7ty5SEtLw8WLF/HUU0+hpqYGCxcuVKI56gOqPZVN7vBY7LLDVVmqSD62\nSqVGnDmx0Tbq8dZkxFmToNPow+rXzSIw+CYiIqLOUCQAP3/+PObNm4fi4mIkJCRgypQpOHjwINLS\n0pRojnqBuq3U6zadKaofbLvsqHZXKNKuTqNvfit1czzUzMcmIiKiLqZIAP7KK68o8bXUw9XfSr2o\n3gy20lupRxktl7ZRrzeLHW9NhtkUw3xsIiIi6lG6JAec+g5PbQ2KXQ55Jrt+2khZhUJbqUNCjDm+\n0eYzdS+j3hTxNomIiIiUwgCcwgghUFnjCts+vX6qSEW1U5F21WoN4i3JjWaw42NSYDMnQqvh+thE\nRETUNzAA74f8AT+cFcX1Aut6S/eVO+CprVGkXYPOVC+4Dp/Fjom2cX1sIiIi6hcYgPdRXn8tCou/\nDqWLXJrBLnHaUVJxEYGAX5F2LabY8BnseikjUdxKnYiIiIgBeG8VniriuPSwo9OOC8XfoMZbCRyM\nfLsqlRo2c0KjZfuCS/clQ681RL5RIiIioj6EAXgP5g/4UVZR1CgPu25VEaV2edRpDeEz2PWC7Vhz\nApfuIyIiIuoEBuDdrMZTjZJyO0rkVJFLK4yUlSuzqggAmI1WxDVMFQnlZZtNVqaKEBERESmEAbjC\nAiIAV2VpaObaETaDXVzuQFVNuSLtSpBgsyQ2mypi0BkVaZeIiIiIWsYAPAI8XjdKXA6UlNcPsEPv\nyx3w+32KtKvT6MO3T7ckIT4mBRfOFSFKb8WkiZMUaZeIiIiIOo4BeBsEZ7FLUOxyoLTcgWKXIzib\nHUodUWptbAAwm2LqzVwntSlVpKooT7H+EBEREVHnMAAPqZ+LXVIvyC5xBZftU2oWW63SwGZJrBdc\nJ4Vms4Pv9UwVISIiIupT+k0AXuv1oLTiYijAvojS8osoKXeEjhdR7a5QrG2TwYx4SxLiY4JpInHy\nQ49JiImO4wY0RERERP1InwnAfX4vSsuLGgXWde+VTBNRSSrEWhKCOdihhxzlmWxrEkz6aMXaJiIi\nIqLeRdEAfOvWrfjVr34Fu92OkSNH4tlnn8XVV1/d7u8RQqDaUwlnRTFKK4pQVlEMZ0UxykLvSysu\nwlVZCgGhwK8IijKY5cDaZklCvDVJDrhjzPFcG5uIiIiI2kSxAPy1117Dgw8+iOeeew5XX301tmzZ\nghtuuAEnTpxAWlpao/oXywrhrAwG1aUNAuyyymLUKrTpTB21WoM4cyLi5Icdk0LpIsGjUR+laPtE\nRERE1D8oFoBv3rwZixcvxpIlSwAAv/vd77Br1y4899xzWL9+faP66176uVJdARBME4kxxwcfeLQk\nIc6SKL+3WRJhjbZBJakU7QMRERERkSIBeG1tLQ4fPoyHH344rHzOnDnYv3+/Ek0CAKxRtmBAbU0M\nBdhJcrAdEx0HtbrPpLwTERERUS8lCSEinjhdWFiIgQMHYu/evWE532vXrsXOnTtx6tQpAIDL5Yp0\n00REREREXcZqtbb7GuZcEBERERF1IUUC8Pj4eKjVajgcjrByh8OBlJQUJZokIiIiIuoVFEmK1ul0\nGD9+PHbv3o1bbrlFLn///fdx6623yp87MmVPRERERNSbKfZU4ooVK3DXXXdh4sSJmDp1Kv74xz/C\nbrfj3nvvVapJIiIiIqIeT7EA/Mc//jFKSkqwbt06XLhwAaNGjcJ7773X5BrgRERERET9hSKroBAR\nERERUdO6bRWUrVu3Ij09HUajEVlZWcjNze2urvQZ2dnZUKlUYa8BAwZ0d7d6lb1792Lu3LkYOHAg\nVCoVcnJyGtXJzs5GamoqTCYTrrnmGpw4caIbetr7tDa2ixYtanT/Tp06tZt62zts2LABEyZMgNVq\nRWJiIubOnYvjx483qsd7tv3aMra8Z9tvy5YtGDNmDKxWK6xWK6ZOnYr33nsvrA7v1/ZrbVx5r0bG\nhg0boFKpcP/994eVd+Se7ZYAvG6b+scffxxHjx7F1KlTccMNN+Cbb77pju70KZmZmbDb7fLr888/\n7+4u9SpVVVUYPXo0fvvb38JoNEKSpLDzzzzzDDZv3ow//OEPOHToEBITEzF79mxUVlZ2U497j9bG\nVpIkzJ49O+z+bfgXM4Xbs2cP7rvvPhw4cAAffvghNBoNrrvuOpSVlcl1eM92TFvGlvds+6WlpWHj\nxo04cuQIPvnkE8yaNQs33XQTPv30UwC8XzuqtXHlvdp5Bw8exLZt2zB69Oiwv786fM+KbjBx4kSx\ndOnSsLKMjAyxcuXK7uhOn7F69WpxxRVXdHc3+ozo6GiRk5Mjfw4EAiI5OVmsX79eLqupqRFms1k8\n//zz3dHFXqvh2AohxMKFC8X3v//9bupR31BZWSnUarV45513hBC8ZyOp4dgKwXs2Umw2m/jTn/7E\n+zXC6sZsTX8NAAAgAElEQVRVCN6rneV0OsVll10mPv74YzFz5kxx//33CyE692dsl8+A121TP2fO\nnLBypbep7y/y8/ORmpqKoUOHYt68eTh79mx3d6nPOHv2LBwOR9i9azAYMH36dN67ESBJEnJzc5GU\nlIThw4dj6dKlKCoq6u5u9Srl5eUIBAKIjY0FwHs2khqOLcB7trP8fj9effVVuN1uTJ8+nfdrhDQc\nV4D3amctXboUt956K2bMmAFR79HJztyziq2C0pzi4mL4/X4kJSWFlScmJsJut3d1d/qUyZMnIycn\nB5mZmXA4HFi3bh2mTp2K48ePw2azdXf3er26+7Ope7ewsLA7utSnXH/99bjllluQnp6Os2fP4vHH\nH8esWbPwySefQKfTdXf3eoUHHngA48aNw5QpUwDwno2khmML8J7tqM8//xxTpkyBx+OB0WjE66+/\njuHDh8sBC+/XjmluXAHeq52xbds25OfnY+fOnQAQln7SmT9juzwAJ+Vcf/318vsrrrgCU6ZMQXp6\nOnJycrB8+fJu7Fnf1zCfmdrvtttuk9+PHDkS48ePx+DBg/Huu+/i5ptv7sae9Q4rVqzA/v37kZub\n26b7kfds2zU3trxnOyYzMxOfffYZXC4X3njjDdx+++346KOPWryG92vrmhvXrKws3qsddPr0aTz2\n2GPIzc2FWq0GAAghwmbBm9PaPdvlKSjcpr7rmEwmjBw5El9++WV3d6VPSE5OBoAm7926cxQ5KSkp\nGDhwIO/fNli+fDlee+01fPjhhxgyZIhcznu285ob26bwnm0brVaLoUOHYty4cVi/fj0mT56MLVu2\nyDEA79eOaW5cm8J7tW0OHDiA4uJijBw5ElqtFlqtFnv37sXWrVuh0+kQHx8PoGP3bJcH4PW3qa/v\n/fff55I4EeZ2u3Hy5En+j02EpKenIzk5OezedbvdyM3N5b2rgKKiIpw/f573byseeOABOUAcNmxY\n2Dnes53T0tg2hfdsx/j9fgQCAd6vEVY3rk3hvdo2N998M44dO4ZPP/0Un376KY4ePYqsrCzMmzcP\nR48eRUZGRofvWXV2dna2wv1vxGKxYPXq1RgwYACMRiPWrVuH3NxcbN++HVartau702c89NBDMBgM\nCAQC+OKLL3DfffchPz8fzz//PMe1jaqqqnDixAnY7Xb8+c9/xqhRo2C1WuH1emG1WuH3+/HLX/4S\nw4cPh9/vx4oVK+BwOPCnP/2JeXStaGlsNRoNVq1aBYvFAp/Ph6NHj+Kee+5BIBDAH/7wB45tM5Yt\nW4aXXnoJb7zxBgYOHIjKykpUVlZCkiTodDpIksR7toNaG9uqqiresx3w6KOPyn9PffPNN3j22Wex\nc+dObNy4EZdddhnv1w5qaVyTk5N5r3aQwWBAQkKC/EpMTMRf//pXDB48GAsXLuzcn7HKLdrSsq1b\nt4ohQ4YIvV4vsrKyxL59+7qrK33G7bffLgYMGCB0Op1ITU0VP/rRj8TJkye7u1u9ykcffSQkSRKS\nJAmVSiW/X7x4sVwnOztbpKSkCIPBIGbOnCmOHz/ejT3uPVoa25qaGvHd735XJCYmCp1OJwYPHiwW\nL14svv322+7udo/WcCzrXmvWrAmrx3u2/VobW96zHbNo0SIxePBgodfrRWJiopg9e7bYvXt3WB3e\nr+3X0rjyXo2s+ssQ1unIPcut6ImIiIiIulC3bUVPRERERNQfMQAnIiIiIupCDMCJiIiIiLoQA3Ai\nIiIioi7EAJyIiIiIqAsxACciIiIi6kIMwImIiIiIuhADcCIiIiKiLsQAnIiIiIioCzEAJyIiIiLq\nQgzAiYiIiIi6EANwIiIiIqIuxACciKgXmDlzJlSqrv0je8iQIUhPT+/SNomI+gMG4EREvYQkSV3e\nXsM2s7OzoVKpkJOT06V9ISLqSzTd3QEiIuqZPvzww2bPdfX/DBAR9SUMwImIqEktpZ8IIbqwJ0RE\nfQtTUIiIOungwYNQqVSYO3dus3WysrKgVqtRUFAgl3300UeYO3cuEhISoNfrMWTIECxbtgwOh6PN\nbQsh8MILL2Dy5Mkwm82IiorClVdeic2bN8Pn8zV5zfnz5/Hggw9i2LBhMJlMsNlsyMrKwurVq8Ou\naZgDPnPmTKxduxYAsHjxYqhUKvlVUFCAlStXQqVS4aWXXmqy3VOnTkGlUmHatGlt/n1ERH2RJDiN\nQUTUaSNGjMCXX36J8+fPIyEhIezc8ePHMWrUKMycOVNO63jmmWewcuVKxMXF4Xvf+x6Sk5Px6aef\n4p///CdSU1Nx8OBBpKamyt8xc+ZM7Nu3D36/P+y7FyxYgB07dmDgwIH44Q9/CK1Wi3/84x84c+YM\n5syZg3fffRdqtVqun5eXh+uvvx6lpaWYPn06Jk+eDLfbjZMnT+Ljjz9GUVERLBYLgGAArlKpkJ+f\nDwDIycnBiy++iD179uCmm27C2LFj5e994IEH4HK5MHToUEyaNAn//ve/G43R8uXL8dvf/hY7duzA\n/PnzOzniRES9mCAiok575plnhCRJ4je/+U2jc//1X/8lJEkSOTk5Qggh9uzZIyRJElOnThUulyus\n7ssvvywkSRK33HJLWPmMGTOESqUKK3v11VeFJEli7NixoqKiQi6vra0V1157rZAkSWzatEku93g8\nYsiQIUKlUomXX365UT8dDofw+Xzy58GDB4v09PSwOqtXrw77LQ394Ac/EJIkic8++yysvKamRsTG\nxoqEhARRW1vb5LVERP0FU1CIiCLgrrvuglqtbrQ6iN/vx44dOxAdHY0f/ehHAIDf/va3AIDnn39e\nnm2uc+edd2Ls2LF46623UFlZ2WKbL7zwAgBgw4YNiI6Olsu1Wi1+85vfAAC2bdsml7/99tsoKCjA\njTfeiDvvvLPR9yUmJobNlnfEz3/+cwDB31bf66+/DqfTiUWLFkGr1XaqDSKi3o4PYRIRRUBKSgpm\nz56NXbt24dNPP8WYMWMAALt374bdbseiRYtgMpkAAP/+97+h0Wjw97//HX/7298afZfH44Hf78cX\nX3yBK6+8stk2Dx8+DEmScM011zQ6N2rUKCQkJODMmTOorq6GyWTCwYMHAQA33HBDJH5yk66//noM\nHToUO3bswMaNG+Xf/Pzzz0OlUuGnP/2pYm0TEfUWDMCJiCJk8eLF2LVrF3JycrB582YAkGfEFy1a\nJNcrKSmB3+/HmjVrmv0uSZJQVVXVYnsulwtWqxV6vb7J8ykpKSguLobL5YLJZILT6QSAsNxyJdx7\n7714+OGH8corr2DJkiX4/PPPceDAAVx33XW47LLLFG2biKg3YAoKEVGE/Md//AdiY2Oxc+dOBAIB\nOJ1OvPXWWxg6dCimT58u17NarbBYLAgEAs2+/H5/q6uFWK1WuFwueDyeJs9fuHBBrgcAMTExAIBv\nv/02Ej+3WXfffTcMBoOchlJ3vPfeexVtl4iot2AATkQUITqdDrfffjsuXryId999F6+99ho8Hg8W\nLFgQVm/q1KkoLy/HZ5991qn2xo8fDyEEPvroo0bnjh07hqKiInmpQQCYMmUKAOB///d/O9xmXY54\nw9VY6rPZbLjtttuQl5eH3Nxc7NixAykpKbjppps63C4RUV/CAJyIKIIWL14MIJh6kpOTA5VKhYUL\nF4bVWbFiBQBg6dKlOH/+fKPvcLvdyM3NbbWtJUuWAABWrVoVlq7i9XrlNu655x65/Ac/+AGGDBmC\n9957Dzt27Gj0fQ6Ho8XAGgDi4uIAIGw986bUPYw5b948lJeXY8mSJVCp+FcOERHAdcCJiCJu1KhR\nOH36NHw+X9ja3/X9+te/xiOPPAKdTocbb7wR6enpqKmpwblz57B3714MHToUhw8fluvPnDkTe/fu\nRSAQCPueO++8Ezt37kRaWhpuuukmaLVavP322zhz5gyuu+467Nq1Kyzw/eSTT/Dd735XXgd80qRJ\nqK2txenTp/Gvf/2rxXXAAeD06dMYOXIkoqOjcddddyEpKQkA8Itf/KLRii4TJkzAJ598ArVajbNn\nz2LgwIGdH1wior6gLWsVFhYWigULFoiEhARhMBjEiBEjxJ49e1q85rPPPhPTp08XRqNRpKamirVr\n13Z+0UQiol5g06ZNQpIkoVKpml0vWwghDh48KObNmycGDhwodDqdiI+PF6NHjxb33Xef2LdvX1jd\nmTNnNloHXAghAoGAeP7558XEiRNFVFSUMBqNYuzYsWLTpk3C6/U22e4333wj7rvvPjF06FCh1+tF\nXFycmDBhglizZk3YNUOGDGm0DrgQQrzyyiti/PjxwmQyyb+zoKCgUb0tW7YISZLED37wg2bHgIio\nP2p1BtzpdOLKK6/E9OnTcd999yEhIQH5+flISUlBZmZmk9eUl5dj2LBhmDlzJp588kmcPHkSixcv\nRnZ2tvzPokRE1Lf95Cc/wZ///Ge88847uPHGG7u7O0REPUarAfiqVauwb98+7Nu3r81f+txzz2Hl\nypVwOBzy8lhPP/00nnvuOcWfviciou53/vx5fOc738HAgQNx5syZ7u4OEVGP0uoTMW+++SYmTpyI\n2267DUlJSRg3bhy2bNnS4jUHDhzAtGnTwtamnTNnDgoLC1t9cIeIiHqvnTt3Ys2aNbjuuutQW1uL\ntWvXdneXiIh6nFY34snPz8fWrVuxYsUKrFq1CkeOHMH9998PAFi2bFmT19jtdgwaNCisrO5BHbvd\njsGDBwMIbiJBRER9x3PPPYf9+/djwIABWLNmDW688Ub+WU9EfVrdXgvt0WoAHggEMHHiRDz99NMA\ngDFjxuDMmTPYsmVLswG4JEnt7ggREfV+77zzTnd3gYiox2s1BWXAgAEYMWJEWFlmZibOnTvX7DXJ\nycmw2+1hZQ6HQz5HRERERNRftToDftVVV+HUqVNhZV988QWGDBnS7DVTpkzBI488Ao/HI+eBv//+\n+0hNTZXTTxrqyPQ9NS0vLw8AkJWV1c096Vs4rsrguCqD46oMjmv7CSEQEAH4Az74/X4EAj74A/7g\n54A/+PL7cOzY5wiIAIYNH4aACJYFRAB+f7Be82WXzvlFAIHQdwbqXsIPfyAgf/YLf4M6Abmsyc8i\ngEDo+kvlgVB5vbJQPSECrQ8KddiPr7kXV4++vru70enUulYD8OXLl2Pq1KlYv349fvzjH+PIkSP4\n/e9/jw0bNsh1Vq5ciUOHDuGDDz4AAMyfPx9r1qzBokWL8Pjjj+P06dN45plnkJ2d3anOEhER9Xf+\ngB8+vxd+vw8+vw/+gDd09NUr88Hn99Z774Pf7730Xi7zwRfwISCXhwLZ0Hf5A374Al75fd05X8CH\ngN8PXyC8bl1QHfDXC64Dvnb9vvc+U2jgqE/oK/+D02oAnpWVhTfffBOrVq3CU089hcGDB2PdunX4\n2c9+Jtex2+1hO6VZLBa8//77WLZsGbKysmCz2fDQQw9h+fLlyvwKIiIihfgDfvh8tfD6vajylMMf\n8KGwuEAOcH3+2tD74GevL/xz8+eDQXGwjhc+OWiuC64bvg8Gvn0lACHqCIG+sYF7qwE4ANx4440t\nbqKwffv2RmVXXHEF9uzZ0/GeERERhQgh5GDW6wu+an218Po88Pm9oc8eObgN++zzwuv3hK7zwuuv\nK2t89IYCbV+9Y6CpgPdw148B9R4qSQVJpYJaUkOlUkMlqeRjsPzS+frH+vXkY6hu3XuVSgVJUrW5\nzsWLRZAgYcCA1LB6kiQ1/j4p9D6sjUvlkhR+rSRJUKnUkCDVu0aCJKkb1FFBwqX3l77vUp363938\n+2B7fUGbAnAiIqLm+P0+eHxueL218Hjd8Po88Hg9oWPwc63XIwfEwWPovdcjB9K1Pg+83lrU+oPH\nYNmlQJozv72XJKmgVqlDLw3UKjVUak3YZ4+nFipJBbPZArUUrKsKnVepVKFj/fJLx0vnVFCpNFDX\nBacN60vq0HfVBcYNv09Vr44m1G9VWN2636JSNQyWg+3XBZY9BZ9b6JkYgBMR9QP+gB+1Xjc8Xnej\nY/h7D2q9NfDUulHrc8NT64bH54HX6wkF2ZeC6VqvG7W+2nbn+FLnSJCgUWuhUWugUWuhVmugVmug\nUWlDR029skt16so1ag3UYXVDQbBaC7VKHTofetWdCwXJ6tA5jTr8sxxINwiqVaHvUEmtLrrGQJH6\nlVYD8Ozs7EY7mSUnJ6OwsLDJ+l9//TWGDh3aqHzXrl2YM2dOB7tJRNR/CCHg9dXCXVsDjzf0qq2B\nx+sOlblDn4OBsrveeU9tDdzeGrjKy+D11+LvnwTg8brh83u7+2f1WhIkaDU6aDQ6iACgltSIioqG\nVq0LC4Q1Gp38vv45rUYHtVobKteGysMD6EZlqsZlcrDdR/4Jnqg/a9MMeGZmJj7++GP5s1rd+n/8\n//znPzFmzBj5c2xsbPt7R0TUiwREAO7aarg9weDYXVsNd23o6Kn3vn55gzJP6HOTecf9nEqlhlaj\ng06tg1ajg1ajDwXG2lBZvc9159XBz1qNHjpNMCgOXquTA2WtRgtNk0cdtGptMMc1lFLAWVoiioQ2\nBeBqtRqJiYnt+mKbzdbua4iIuosQIhQMV6HGU4UaTzAorvFUoSYUQNeEPrvDyqrkc+7a6u7+Gd1C\nklTQafXQawzQanXQawzQaQ3QaXTBo1YfCoD10Gl10Kr10Gr1cpAsH+uVaevVDV6v48wvEfUZbQrA\n8/PzkZqaCr1ej0mTJmH9+vVIT09v8Zof/vCHcLvdyMjIwPLly3HLLbdEpMNERM3x+32oqa1GtbsS\nNZ5KVHuCwXTwcxWqPfXK3XWfq+R6ffkhPwlSMEjWGqEPBcV6rRE6nQF6jR56nRE6jR56nQE6rRH6\nuvNaPXQaQ+iolwNqnUYPvdYArUYPjVrTox46IyLq6VoNwCdPnoycnBxkZmbC4XBg3bp1mDp1Ko4f\nPw6bzdaovtlsxq9//WtcddVV0Gg0eOutt3DbbbchJycHd9xxhyI/goj6Fq/Pi2pPBardlah2V6DK\nXYlqdyWq3BWodgfLqzwVqK6pQJWnUq7n8bq7u+sRo1ZrYNAaodcZYQgFynWf9VoDDDqjHEzrdUYY\ndEbotEa5/KszX0Gj1iHryomhQFnHIJmIqIeQhBDtWtG8uroa6enpePTRR9u8sc59992Hffv24dNP\nPw0rr7+N55kzZ9rTDSLqJXx+L9y+ani8wZfbV9P4va8aHp8bHm8Nan018AV67wODGlUonUKtD6Za\naEIpF+pgmUYdSquoVx6sF8xXrvvMdAsiop4rIyNDfm+1Wtt9fbuXITSZTBg5ciS+/PLLNl8zYcIE\n/OUvf2lvU0TUwwgh4PHVwO2tCr2q4Q4F08HyUDDtvfS+Ny1Rp1ZpoFMboNPooW10DKVgaAyhYNoQ\nyl82BIPtUJ22LLdGRET9W7sDcLfbjZMnT2LWrFltvubo0aMYMGBAi3X4RHnk8Cl9ZfTFcRVCwON1\no6LaicqaclTWuFBR7UJltRMVNS5UVrvCjlU15T16dQ4JEgx6E0z6aBj1UTDpo2A01L2/VGYKlRn1\n0cH3uigY9VHQarTd/RMipi/erz0Bx1U5HFtlcFyVUT+LoyNaDcAfeughzJ07F2lpabh48SKeeuop\n1NTUYOHChQCAlStX4tChQ/jggw8AADk5OdDpdBg7dixUKhXefvttbN26FRs3buxUR4mo7by+WpRX\nl6G8yomK0LG8ugwVoWN5tRPlVWWorHbB66/t7u42IkkqmAzRiNJHw2Q0I0pvhskQDJajDGaYDGZE\nGaJhMphh0kcjyhg8GvQmzkATEVGP12oAfv78ecybNw/FxcVISEjAlClTcPDgQaSlpQEA7HY78vPz\n5fqSJGHdunUoKCiAWq3G8OHDsX37dsyfP1+5X0HUDwghUO2phKuyFK6qUpRXBQPpitCxfoBd46nq\n7u7K1CoNoo0WRBktiDaYEWW0ICp0jDZa4LhQDL3GiHGjx8sBtl5nZCBNRER9VqsB+CuvvNLi+e3b\nt4d9XrBgARYsWNC5XhH1M7U+D1yVpSivKoWzshSuqhI50HZVlsJZVYLyyrIeMVtt1JkQbYqB2WhF\ntMmKaKMZ0UYrogwWRBnNwWDbYJGDbr3W0OLqG3m+4D+PDkr6Tlf9BCIiom7V7hxwImqfGk8VyiqK\nUVZRBGdlCcoqikOBdYkcYFd7Krutfxq1FuZ6AfWlwNoKc4NjtNHap/KkiYiIugMDcKJO8Pq8cFYW\nw1lZHAqyQ4F2RTHKQmXdsTuiSlIFg+qoGFhMsbCYYmCJioXZFANLlA0WU0zofWyrM9REREQUWQzA\niVpQ46lGabkDJeUOnCj8BFVuFz61f4CyyhKUVRShotrZpf3RafSwRsfBGhULa5QN5qjw4NoaFQuz\nKRZRRjNzqImIiHqoVgPw7OxsrF27NqwsOTkZhYWFzV7z+eef47777sOhQ4dgs9nw05/+FE888UTn\ne0sUYV5fLUrLL6Kk3IESlwMloffBsouodld0ST9UkgqWUFBtjY5DTLQNligbYqLjQmU2WKNsMOhM\nnK0mIiLq5do0A56ZmYmPP/5Y/qxWN79DW3l5OWbPno2ZM2ciLy8PJ0+exOLFixEVFYUVK1Z0usNE\n7REQAbgqS1DktKPEZW8QYDtQXlWmeB80ai1io+MRa45HjDl4tEbFyUF1THQcoo0WqLjzIRERUb/Q\npgBcrVYjMTGxTV/417/+FW63Gzk5OdDr9RgxYgROnTqFzZs3MwAnRQghUF5VhovOQhQ5L6Co3rHY\nZYfXp9zKIRIkWKJtcoAdfCUgpt7naKOVs9ZEREQka1MAnp+fj9TUVOj1ekyaNAnr169Henp6k3UP\nHDiAadOmQa/Xy2Vz5szBE088gYKCAgwePDgyPad+RQiBimpXWHAtH1121HrdirSrVmtgMycizpKI\ngEeFKIMVozLHyoG2NcoGtZqPUhAREVHbtRo5TJ48GTk5OcjMzITD4cC6deswdepUHD9+HDabrVF9\nu92OQYMGhZUlJSXJ5xiAU0sCAT+KXQ7YS7+BveQc7KXfwl72DYqcF+CprYl4e5KkQkx0HOIsiYiz\nJMFmTQq9T4TNkgRrtE1+mFHezjeT2/kSERFRx0lCCNGeC6qrq5Geno5HH30Uy5cvb3T+u9/9LtLS\n0vDCCy/IZefOncOQIUNw4MABTJo0SS53uVzy+zNnznSk/9RLBQJ+VLjL4KwugqumOHisLoarpgQB\n4Y9oWzq1AWajDWZDDKINMYjWB49mfQxMeivUzL0mIiKidsjIyJDfW63Wdl/f7n87N5lMGDlyJL78\n8ssmzycnJ8Nut4eVORwO+Rz1L/6AHxXuUjiri+oF2cUorylBQAQi1o5WrYPZYIPFaJOPFoMNZqMN\neo2ROdhERETUY7Q7AHe73Th58iRmzZrV5PkpU6bgkUcegcfjkfPA33//faSmpraYfpKVxX/WjxQ5\nVaKLx7TaXYlvi87i26J8fFuUj/NFZ+Eo/TZigbZOa0BCTAoSYlKQGDMACTEpiLemICFmAMwm5R90\n7K5x7es4rsrguCqD46ocjq0yOK7KqJ/F0RGtBuAPPfQQ5s6di7S0NFy8eBFPPfUUampqsHDhQgDA\nypUrcejQIXzwwQcAgPnz52PNmjVYtGgRHn/8cZw+fRrPPPMMsrOzO9VR6jmEEHBWFoeC7bM4X5SP\nby/mo7SiKCLfbzZakRw3CMm2NCTbBiLJloYkWyospljOZBMREVGv12oAfv78ecybNw/FxcVISEjA\nlClTcPDgQaSlpQEIPliZn58v17dYLHj//fexbNkyZGVlwWaz4aGHHmoyX5x6vkDAj4vOQnx7MT8U\nbAdnuKsisEGNNcoWDLLj0sKC7WijJQI9JyIiIuqZWg3AX3nllRbPb9++vVHZFVdcgT179nS8V9Qt\nhBAorbiIs4WncPbCaXxz8SsUFn+NWp+nU98bGx2PJDnIvhRsmwzREeo5ERERUe/BBYz7Ma/Pi2+L\nvsLZC6eCQbf9dKd2hlSp1Ei2pWFgQjpSE9IxMGEoUhOGwKRnoE1ERERUhwF4P1JeVYazF04HA+4L\np3Du4pfw+30d+i6d1oDU+CGhIDsdAxPSkRI3CFqNLsK9JiIiIupbGID3UQERwLdF+XI6ydkLp1BS\n7ujQd0UbrRgYmtEemBgMuBOsyVBx/WwiIiKidmtXAL5hwwY89thjWLZsGX7/+983Wefrr7/G0KFD\nG5Xv2rULc+bM6VgvqVVCCBQWf42TBUdw6Ng+FFV8C99+b7u/R6vRYXBSBtJTMjEkZTgGJX4Hliiu\nPkJEREQUKW0OwA8ePIht27Zh9OjRbQrG/vnPf2LMmDHy59jY2I71kJpVUe3C6XNHcercUZwqOIry\n6vbnb8eaE5Cekon0lOFIT8lEavwQqNX8hxEiIiIipbQp0nK5XLjzzjuxffv2Nq/nbbPZkJiY2Jm+\nUQN+vw9n7adxquAIThYcwbcX8yEg2ny9WqXBwIT0YMA9IBNDkocj1hyvYI+JiIiIqKE2BeBLly7F\nrbfeihkzZkCItgV8P/zhD+F2u5GRkYHly5fjlltu6VRH+6sSlwMnC47g1LkjOP3NZ/DU1rT52mij\nVZ7ZTk/JRFrSZdBp9Ar2loiIiIhaI4lWIupt27bhT3/6Ew4ePAi1Wo1rrrkGo0aNwu9+97sm65eU\nlOCll17CVVddBY1Gg7feegtPP/00cnJycMcdd4TVrb+N55kzZyLwc3o/r78WDlcBzju/QmFZPirc\npW2+VqPSItk6BANiL8OAmHSYDTbmbhMRERFFWEZGhvzearW2+/oWZ8BPnz6Nxx57DLm5uVCrgyte\nCCFanAWPi4sL2/XyyiuvRElJCTZu3NgoAKcgn9+Lb8vOIL/oGArLvkJA+Nt8rS0qGQNihmJA7FAk\nmNOg5sokRERERD1aizPgL774Iu6++245+AYAv98PSZKgVqtRVVUFrVbbaiM5OTn42c9+hurq6rDy\n+vSBSUAAABb+SURBVDPgHfm/h94sEPDjzLfHkHd6Lz798gDctdWtXwTAbLRi+OCxuHzwOAxPGwtL\nVEyjOnl5eQCArKysiPa5v+O4KoPjqgyOqzI4rsrh2CqD46qMzsawLc6A33zzzZg4caL8WQiBxYsX\nY9iwYVi1alWbgm8AOHr0KAYMGNDuzvU1QgicLz6LvFN78cnpvXBVtZ5eolZpkD4gE5cPGofMweOQ\nmjAEKknVBb0lIiIiIiW0GIBbrdZGUb3JZEJsbCxGjBgBAFi5ciUOHTqEDz74AEBwtlun02Hs2LFQ\nqVR4++23sXXrVmzcuFGhn9DzlZZfRN7pvcg7tQf20m9arR9nTcKIweOROXgsMgaOgkFn7IJeEhER\nEVFXaPeCz5IkhT3YZ7fbkZ+fH3Z+3bp1KCgogFqtxvDhw7F9+3bMnz8/Mj3uJarcFTh6Zj/yTu3B\nV4UnWq1vNsVg/LBpyMqcgbTEy/jwJBEREVEf1e4A/KOPPgr7vH379rDPCxYswIIFCzrXq17K66vF\nsbN5yDv1MU58fRj+gK/F+nqtAWO+MwXjh0/HsLTRfICSiIiIqB/glocRUFZRhA/y/geHTn3c6sOU\nKpUalw8ah6zMGRg1dCJ0Wq7LTURERNSfMADvhNLyi3j/0N9x8MS/Wp3tHpI8HFmZMzAu4yqYTf1r\nxRciIiIiuqRdy2ls2LABKpUK999/f4v1Pv/8c8yYMQMmkwkDBw7EU0891alO9jQl5Q68+q8teCrn\n5/j3sX82G3wnxgzADZPn4YmFz2HFbc9g+pgbGXwTERER9XNtngE/ePAgtm3bhtGjR7f4gGB5eTlm\nz56NmTNnIi8vDydPnsTixYsRFRWFFStWRKTT3aXYZcfuQ3/D/zv5EQKBpjfLMRutuHL4NGQNn4FB\nSd/hw5REREREFKZNAbjL5cKdd96J7du3Izs7u8W6f/3rX+F2u5GTkwO9Xo8RI0bg1KlT2Lx5c68N\nwIucF7D70N9w6ORHCIhAk3XiLEmYM+FHmHD5TGjUbVsfnYiIiIj6nzYF4EuXLsWtt96KGTNmtLgN\nPQAcOHAA06ZNg15/6eHCOXPm4IknnkBBQQEGDx7cuR53oYtlhdh96A3kndrTbOAdb03GnAm3YkLm\nDKjVTKknIiIiopa1GjFu27YN+fn52LlzJwC0mlJht9sxaNCgsLKkpCT5XG8IwB1l5/HP//c6Pjm9\nD6KZwDvBmoI5E29FVuYMLh9IRERERG3WYgB++vRpPPbYY8jNzYVaHQwyhRAtzoJ3NOc5Ly+vQ9dF\nkrO6GJ9/sw9fF5+AQNO/0WKwYXTaNAxJGAlVtQpHDh/p4l62XU8Y076I46oMjqsyOK7K4Lgqh2Or\nDI5rZGVkZHTq+hYD8AMHDqC4uBgjR46Uy/x+P/bt24fnn38eVVVV0GrD852Tk5Nht9vDyhwOh3yu\nJ3JWF+GzUODdHKsxDqPSpmFI/AiopHYtHkNEREREJGsxAL/55psxceJE+bMQAosXL8awYcOwatWq\nRsE3AEyZMgWPPPIIPB6PnAf+/vvvIzU1tcX0k6ysrI7+hg4LiAD+lfc/ePfoX5vN8U62peG7E3+M\ncRlToeolqSZ1/5fbHWPal3FclcFxVQbHVRkcV+VwbJXBcVWGy+Xq1PX/v717D2ryXvMA/k2AcKuk\nUA0XUUBPhCo1UrEaqngFbXWol/aI1FE4p+tOCw6FdXekdSpOOVo5W9dOBUdtF6P1gmyd9thSBA+I\nsNCKVpDipbSopWpSbS0WKqjJu390zPrKJRcCMfj9zGQm+b2/N3l45hl5eP3l9/bYgMvlcsjl4n2r\nPTw84O3tjdGjRwMAMjIyUFNTgyNHjgAAEhISsG7dOiQmJmLNmjU4f/48Nm7caHL3lP72e3srdhdv\nRsOFrv9Lxv+J4Zj9zJ8xThnFK95EREREZDMWb9shkUhE67y1Wi2ampqMr728vFBSUoLk5GRERkbC\nx8cHq1atQlpamm0itoHmn77Hf3+ejZ9v6jodC3giCLMnLobqT5PYeBMRERGRzVncgJeVlYle5+Xl\ndZoTHh6O8vJy66PqI4IgoOqbYnxc/gHu6u+Ijnm4DcKfp/8rr3gTERERUZ96ZDauvn2nA/mlW1Fz\n7minY8N9lfjL8/8OHy9F/wdGRERERI+UR6IB/+nGZXz4+UZc/fmHTsemjH0e86ckwcWZd68kIiIi\nor434BvwU41V2HvkfXTcviUal7m4YcnM1zA+NNpOkRERERHRo8jkYuecnByoVCrjjihRUVEoLCzs\ndv7FixchlUo7PYqLi20auCl6/V0cLP8QeYXZnZpvX59A/Nviv7P5JiIiIqJ+Z/IK+LBhw5CdnQ2l\nUgmDwYCdO3di/vz5qKmpgUql6va8w4cPi457e3vbJmIz3PjtOnZ+8Z+4cPVcp2NPj5qCJTNfg6vM\nvd/iISIiIiK6x2QDHhcXJ3qdlZWFrVu34vjx4z024D4+PlAo+v9Ljed/qIOmaBNab4k3SHeSOmNB\n9F8wZexzom0UiYiIiIj6k0VrwPV6PQoKCtDe3o7o6J6XbyxcuBDt7e1QKpVIS0vDokWLehWoKQbB\ngOLjBfjiy/0QIIiOeT82GElz/wPBfqP6NAYiIiIiIlMkgiAIpibV19dDrVajo6MD7u7u2LdvH+bO\nndvl3J9//hm7du3Cs88+C2dnZ3z66af429/+Bo1Gg5dfflk09/7beDY2Nlr9Q7Tf+R3/2/gpLt/4\nvtOxgMdHYvKoF+Dm4mH1+xMRERER3aNUKo3PH7xrvDnMasDv3LmD5uZmtLS0oKCgAO+//z7KysoQ\nGRlp1oekpKSgoqICdXV1onFbNODXf7uM8vMfo63jZqdj44ZPxVOBk7nkhIiIiIhspl8a8AfFxMQg\nMDCwy7tgdkWj0eDVV1/F77//Lhq/vwG3JvhTjVXYVbQJesNd0binuxeWz05HWNA4i99zIDhx4gQA\nmP0HEpmHee0bzGvfYF77BvPad5jbvsG89o3e9rBW7QOu1+thMBjMnl9bW4uAgABrPqpb1369ij3F\n73VqvoP9QpH0/Cp4Dxpi088jIiIiIrIFkw346tWrMW/ePAQGBuK3337D3r17UV5ejqKiIgBARkYG\nampqcOTIEQB/XO2WyWQYN24cpFIpDh06hNzcXGRnZ9ssaL1Bj93Fm3H7bodofOq4eXhh8nI4O/Gu\nlkRERET0cDLZgOt0OixduhRarRZyuRwqlQpFRUWIiYkBAGi1WjQ1NRnnSyQSZGVl4dKlS3ByckJo\naCjy8vKQkJBgs6D/eeIgLl49Lxp7cdq/IFrV9RdDiYiIiIgeFiYbcFPrvB88vmzZMixbtqx3UfWg\n+acmFH61XzQWHjIBU8Y+32efSURERERkKyZvRf8wuXP3NnYf/i8YDHrjmKe7F+JnJnOnEyIiIiJy\nCA7VgB+q+gjaX5pFY/EzXoOX5+N2ioiIiIiIyDIO04B/21yPo6f+IRqb+OQMqP40yU4RERERERFZ\nzmQDnpOTA5VKBblcDrlcjqioKBQWFvZ4Tn19PaZOnQoPDw8EBgbi7bff7lWQtzrasKf4PdGYz6Ah\nWDj1lV69LxERERFRfzP5Jcxhw4YhOzsbSqUSBoMBO3fuxPz581FTUwOVStVp/s2bNxETE4Np06bh\nxIkTOHv2LJKSkuDp6Yn09HSrgvyfoztwo/W68bUEErwcmwp3V95enoiIiIgci8kGPC4uTvQ6KysL\nW7duxfHjx7tswPfs2YP29nZoNBq4urpi9OjROHfuHDZt2mRVA17bWIWac0dFY9OffgHKwHCL34uI\niIiIyN4sWgOu1+uxf/9+tLe3Izo6uss51dXVmDJlClxdXY1jsbGxuHLlCi5dumRRcC1tvyC/dKto\nzP+J4Zirtt2e4kRERERE/UkiCIJgalJ9fT3UajU6Ojrg7u6Offv2Ye7crm96Exsbi+HDh+ODDz4w\njv3www8IDg5GdXU1Jk6caBxvaWkxPm9sbBS9jyAIKD2bj8s3vjOOSSVSPK/6K3w8fc3/CYmIiIiI\nbEipVBqfy+Vyi8836wp4WFgYTp8+jePHjyMlJQXx8fE4ceJEl3NttR93o+5rUfMNAKrhU9l8ExER\nEZFDM7kGHABcXFwwYsQIAEBERARqamqQk5PT5V0y/fz8oNVqRWM6nc54rDuRkZHG59d+vYr9X/1d\ndHyE/5NYHpcCqdTJnJAfaff+OLo/p9R7zGvfYF77BvPaN5jXvsPc9g3mtW/cv4rDGlbtA67X62Ew\nGLo8plarUVFRgY6ODuNYSUkJhg4diqCgINPvbdBjd/Fm3L77/+e7urhh6exUNt9ERERE5PBMNuCr\nV69GZWUlLl68iPr6emRkZKC8vBxLly4FAGRkZGDWrFnG+QkJCfDw8EBiYiIaGhpw8OBBbNy40ewd\nUP554iAuXj0vGlsQ/VcMlnd/9ZyIiIiIyFGYXIKi0+mwdOlSaLVayOVyqFQqFBUVISYmBgCg1WrR\n1NRknO/l5YWSkhIkJycjMjISPj4+WLVqFdLS0kwG0/zT9yj8ar9oLDxkAtRjZnVzBhERERGRYzHZ\ngHe1ztvU8fDwcJSXl1sUyJ27t7H78GYYDHrjmKe7F+JnJtvsi51ERERERPZm1RrwvnCo6iNof2kW\njS2Z+Rq8PB+3U0RERERERLb30DTgR0/9Q/R64pMzMHbkJDtFQ0RERETUN0w24Bs2bMCECRMgl8uh\nUCgQFxeHhoaGHs+5ePEipFJpp0dxcbFZQfkMGoKFU18x7ycgIiIiInIgJhvw8vJypKSkoLq6GqWl\npXB2dsasWbNw48YNk29++PBhaLVa42P69Okmz5FAgpdjU+Hu6mHeT0BERERE5EBMfgmzqKhI9Hr3\n7t2Qy+Woqqrq9nb09/j4+EChUFgU0PSnX4AyMNyic4iIiIiIHIXFa8Bv3rwJg8EAb29vk3MXLlwI\nX19fTJ48GR9//LHJ+f5PDMdcdYKlIREREREROQyLG/DU1FRERERArVZ3O2fQoEF49913UVBQgC++\n+AIzZ87E4sWLsWfPnm7PcZI6Y9nsNLg4yywNiYiIiIjIYUgEQRDMnZyeno4DBw6gsrISwcHBFn1Q\nSkoKKioqUFdXZxxraWmx6D2IiIiIiB4mcrnc4nPMvgKelpaG/Px8lJaWWtx8A8CECRPQ2Nho8XlE\nRERERAOJyS9hAn8sOykoKEBZWRlGjRpl1QfV1tYiICDAqnOJiIiIiAYKkw14cnIyPvroI3zyySeQ\ny+XQarUA/ljn7enpCQDIyMhATU0Njhw5AgDQaDSQyWQYN24cpFIpDh06hNzcXGRnZ4ve25pL9kRE\nREREjsxkA75161ZIJBLMnDlTNJ6ZmYm33noLAKDVatHU1GQ8JpFIkJWVhUuXLsHJyQmhoaHIy8tD\nQgJ3OCEiIiKiR5tFX8IkIiIiIqLesXgbQlvJzc1FSEgI3N3dERkZicrKSnuFMmBkZmZCKpWKHlx3\nb5ljx44hLi4OgYGBkEql0Gg0neZkZmZi6NCh8PDwwPTp03HmzBk7ROp4TOU2MTGxU/1GRUXZKVrH\nsGHDBkyYMAFyuRwKhQJxcXFoaGjoNI81azlzcsuatVxOTg5UKhXkcjnkcjmioqJQWFgomsN6tZyp\nvLJWbWPDhg2QSqVYuXKlaNyamrVLA56fn4/XX38da9asQW1tLaKiovDcc8+hubnZHuEMKGFhYdBq\ntcZHfX29vUNyKG1tbRg7dizee+89uLu7QyKRiI5v3LgRmzZtwpYtW1BTUwOFQoGYmBi0trbaKWLH\nYSq3EokEMTExovp98BcziZWXlyMlJQXV1dUoLS2Fs7MzZs2ahRs3bhjnsGatY05uWbOWGzZsGLKz\ns3Hq1CmcPHkSM2bMwPz5841bFLNerWMqr6zV3vvyyy+xY8cOjB07VvT7y+qaFezgmWeeEVasWCEa\nUyqVQkZGhj3CGTDWrl0rhIeH2zuMAeOxxx4TNBqN8bXBYBD8/PyE9evXG8du3bolDBo0SNi2bZs9\nQnRYD+ZWEARh+fLlwrx58+wU0cDQ2toqODk5CZ999pkgCKxZW3owt4LAmrUVHx8fYfv27axXG7uX\nV0FgrfbWr7/+KowcOVI4evSoMG3aNGHlypWCIPTu39h+vwJ++/ZtfP3114iNjRWNx8bGoqqqqr/D\nGXCampowdOhQjBgxAkuWLMGFCxfsHdKAceHCBeh0OlHturm5ITo6mrVrAxKJBJWVlfD19UVoaChW\nrFiBa9eu2Tssh3Lz5k0YDAZ4e3sDYM3a0oO5BVizvaXX67F//360t7cjOjqa9WojD+YVYK321ooV\nK/DSSy9h6tSpEO776mRvatasfcBt6fr169Dr9fD19RWNKxQK4xaHZJ1JkyZBo9EgLCwMOp0OWVlZ\niIqKQkNDA3x8fOwdnsO7V59d1e6VK1fsEdKAMmfOHCxatAghISG4cOEC1qxZgxkzZuDkyZOQyWT2\nDs8hpKamIiIiAmq1GgBr1pYezC3AmrVWfX091Go1Ojo64O7ujgMHDiA0NNTYsLBerdNdXgHWam/s\n2LEDTU1N2Lt3LwCIlp/05t/Yfm/Aqe/MmTPH+Dw8PBxqtRohISHQaDRIS0uzY2QD34Prmclyixcv\nNj4fM2YMxo8fj6CgIHz++edYsGCBHSNzDOnp6aiqqkJlZaVZ9ciaNV93uWXNWicsLAynT59GS0sL\nCgoKEB8fj7Kysh7PYb2a1l1eIyMjWatWOn/+PN58801UVlbCyckJACAIgugqeHdM1Wy/L0EZPHgw\nnJycoNPpROM6nQ7+/v79Hc6A5uHhgTFjxuC7776zdygDgp+fHwB0Wbv3jpHt+Pv7IzAwkPVrhrS0\nNOTn56O0tBTBwcHGcdZs73WX266wZs3j4uKCESNGICIiAuvXr8ekSZOQk5Nj7AFYr9bpLq9dYa2a\np7q6GtevX8eYMWPg4uICFxcXHDt2DLm5uZDJZBg8eDAA62q23xtwmUyG8ePHo7i4WDReUlLCLXFs\nrL29HWfPnuUfNjYSEhICPz8/Ue22t7ejsrKStdsHrl27hsuXL7N+TUhNTTU2iKNGjRIdY832Tk+5\n7Qpr1jp6vR4Gg4H1amP38toV1qp5FixYgG+++QZ1dXWoq6tDbW0tIiMjsWTJEtTW1kKpVFpds06Z\nmZmZfRx/J15eXli7di0CAgLg7u6OrKwsVFZWIi8vj7en74VVq1bBzc0NBoMB3377LVJSUtDU1IRt\n27Yxr2Zqa2vDmTNnoNVq8eGHH+Kpp56CXC7HnTt3IJfLodfr8c477yA0NBR6vR7p6enQ6XTYvn07\n19GZ0FNunZ2d8cYbb8DLywt3795FbW0tXnnlFRgMBmzZsoW57UZycjJ27dqFgoICBAYGorW1Fa2t\nrZBIJJDJZJBIJKxZK5nKbVtbG2vWCqtXrzb+nmpubsbmzZuxd+9eZGdnY+TIkaxXK/WUVz8/P9aq\nldzc3DBkyBDjQ6FQYM+ePQgKCsLy5ct7929s323a0rPc3FwhODhYcHV1FSIjI4WKigp7hTJgxMfH\nCwEBAYJMJhOGDh0qvPjii8LZs2ftHZZDKSsrEyQSiSCRSASpVGp8npSUZJyTmZkp+Pv7C25ubsK0\nadOEhoYGO0bsOHrK7a1bt4TZs2cLCoVCkMlkQlBQkJCUlCT8+OOP9g77ofZgLu891q1bJ5rHmrWc\nqdyyZq2TmJgoBAUFCa6uroJCoRBiYmKE4uJi0RzWq+V6yitr1bbu34bwHmtqlreiJyIiIiLqR3a7\nFT0RERER0aOIDTgRERERUT9iA05ERERE1I/YgBMRERER9SM24ERERERE/YgNOBERERFRP2IDTkRE\nRETUj9iAExERERH1IzbgRERERET96P8AHFK+m1F6fwcAAAAASUVORK5CYII=\n", - "text/plain": [ - "" - ] - }, - "metadata": {}, - "output_type": "display_data" - } - ], - "source": [ - "from __future__ import (absolute_import, division, print_function,\n", - " unicode_literals)\n", - "\n", - "from numpy import array\n", - "import matplotlib.pyplot as plt\n", - "\n", - "from filterpy.hinfinity import HInfinityFilter\n", - "\n", - "dt = 0.1\n", - "f = HInfinityFilter(2, 1, dim_u=1, gamma=.01)\n", - "\n", - "f.F = array([[1., dt],\n", - " [0., 1.]])\n", - "\n", - "f.H = array([[0., 1.]])\n", - "f.G = array([[dt**2 / 2, dt]]).T\n", - "\n", - "f.P = 0.01\n", - "f.W = array([[0.0003, 0.005],\n", - " [0.0050, 0.100]])/ 1000 #process noise\n", - "\n", - "f.V = 0.01\n", - "f.Q = 0.01\n", - "u = 1. #acceleration of 1 f/sec**2\n", - "\n", - "xs = []\n", - "vs = []\n", - "\n", - "for i in range(1,40):\n", - " f.update (5)\n", - " #print(f.x.T)\n", - " xs.append(f.x[0,0])\n", - " vs.append(f.x[1,0])\n", - " f.predict(u=u)\n", - "\n", - "plt.subplot(211)\n", - "plt.plot(xs)\n", - "plt.title('position')\n", - "plt.subplot(212)\n", - "plt.plot(vs) \n", - "plt.title('velocity')\n", - "plt.show()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.4.3" - } - }, - "nbformat": 4, - "nbformat_minor": 0 -} +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# H Infinity filter" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "#format the book\n", + "%matplotlib inline\n", + "from __future__ import division, print_function\n", + "import book_format\n", + "book_format.load_style()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "I am still mulling over how to write this chapter. In the meantime, Professor Dan Simon at Clevant State University has an accessible introduction here:\n", + "\n", + "http://academic.csuohio.edu/simond/courses/eec641/hinfinity.pdf\n", + "\n", + "In one sentence the $H_\\infty$ (H infinity) filter is like a Kalman filter, but it is robust in the face of non-Gaussian, non-predictable inputs.\n", + "\n", + "\n", + "My FilterPy library contains an H-Infinity filter. I've pasted some test code below which implements the filter designed by Simon in the article above. Hope it helps." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuAAAAEjCAYAAABzSrSOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl8VNXdP/DPnX0mmZlkshMCBBuIIJuEVVlEoWpbHq21\nCiqLWGqLVuHxUcGFgAiVUmoXsBZbjFLc2ufRujwUrQqkwPMjAiqraDAoYYZsM1lnMsv5/TGTSyb7\nMjfr5/16zevOnHvunDPndYWvh+89RxJCCBARERERUZdQdXcHiIiIiIj6EwbgRERERERdiAE4ERER\nEVEXYgBORERERNSFGIATEREREXUhBuBERERERF2IATgRUR/y4osvQqVSIScnp13XqVQqXHPNNQr1\nioiI6mMATkTUh0iSJL/qU6lUSE9Pb/VaIiJSnqa7O0BERJFz8803Y8qUKUhOTm50rqUA+9SpUzCZ\nTEp2jYiIQhiAExH1IRaLBRaLpd3XDRs2TIHeEBFRU5iCQkTUSV9//bWcQ33+/HnccccdSEhIgMlk\nwoQJE/D66683ukYIgRdeeAGTJ0+G2WxGVFQUrrzySmzevBk+n69R/c8++wzz589Heno6jEYjEhIS\nMGbMGCxbtgzl5eVyvYY54B9//DFUKlVYP+teixcvlq9rLge8oqICjz/+ODIzM2E0GhEbG4trr70W\n//jHP1och5KSEixduhQpKSkwGAy44oor8OKLL7Z7bImI+iLOgBMRRUhZWRmuvvpq2Gw23HPPPSgt\nLcXrr7+O22+/HYWFhXjwwQflugsXLsSOHTswcOBALFmyBFqtFv/4xz/w0EMPYffu3Xj33XehVqsB\nBIPvSZMmQa1W4/vf/z4uu+wyVFZWIj8/Hzk5OfjP//zPRrPedekm6enpWL16NdasWQOr1Yrly5fL\ndcaOHdvkNXVcLheuvvpqHD9+HFdeeSUefPBBlJWV4Y033sBNN92ENWvW4Iknnmg0Dk6nE1dddRX0\nej1+/OMfw+Px4PXXX8fdd98NlUqFBQsWdG6giYh6O0FERJ1y9uxZIUmSkCRJ3H777WHnvvrqKxET\nEyP0er0oKCgQQgjx6quvCkmSxNixY0VFRYVct7a2Vlx77bVCkiSxadMmuXzFihVCkiTx1ltvNWq7\nsrJSeDwe+fP27duFJEkiJycnrJ4kSSI9Pb3Z3yBJkrjmmmvCyu69914hSZJYsmRJWPm3334rUlJS\nhEqlEocOHWpyHH7yk5+IQCAgnztx4oTQaDRixIgRzfaBiKi/YAoKEVGEaDQa/PKXvwwrGzp0KH7+\n85+jtrYWO3bsAAC88MILAIANGzYgOjparqvVavGb3/wGALBt27ZG328wGBqVRUVFQafTRew31PF6\nvXjppZcQFRWFjRs3hp1LTU3FqlWr5DSapvq0efPmsBn1yy+/HFOnTsWpU6dQXV0d8f4SEfUmDMCJ\niCJk0KBBGDx4cKPy6dOnAwCOHj0KADh8+DAkSWoy53rUqFFISEjAmTNn5ED19ttvh1qtxk033YS7\n7roL27dvx6lTpxT8JcFVUWpqajBq1CjYbLZG56+77joAwJEjRxqdy8jICPsfizppaWkQQqCsrCzy\nHSYi6kUYgBMRRUhSUlKL5S6XSz5arVbo9fom66ekpITVnzBhAnJzczF79mz8z//8D5YsWYIRI0Zg\n6NCheP755yP9M8Labmo5w/rlTqez0bmYmJgmr9Fogo8d+f3+SHSRiKjXYgBORBQhDoejxXKr1Sof\nXS4XPB5Pk/UvXLgQVh8AJk2ahLfeegtOpxMHDx7EU089hZqaGvzsZz/Dyy+/HMmfEda23W5vcx+J\niKhtGIATEUXIuXPnUFBQ0Kh8z549AIBx48YBAMaPHw8hBD766KNGdY8dO4aioiIMGzasyY1xNBoN\nJk6ciMceewwvvfQSAODNN99stW+SJLVr5vnyyy+H0WjE559/jpKSkkbn//Wvf8m/hYiI2ocBOBFR\nhPh8PjzyyCMQQshlX331FbZu3QqdToc77rgDALBkyRIAwKpVq1BVVSXX9Xq9WLFiBQDgnnvukcv3\n798Pt9vdqL26Wei27GAZFxeHoqKiJr+nKRqNBgsWLEBVVRVWrlwZdq6wsBAbNmyASqXC3Xff3abv\nIyKiS7gOOBFRhIwePRr/93//h6ysLMyePVteB7yiogKbN2/GoEGDAAC33XYb3n77bezcuRMjRozA\nTTfdBK1Wi7fffhtnzpzBddddF7Ze98aNG/Hhhx9i2rRpGDJkCCwWC7744gu88847MJlMYeuLN2fO\nnDnYuXMnrr/+ekybNg16vR5jx47F97///Wav+eUvf4l9+/bhhRdewJEjR3DttdfC6XTijTfegNPp\nxJNPPokJEyZ0fuCIiPqZds+Ab9iwARMmTIDVakViYiLmzp2L48ePN6qXnZ2N1NRUmEwmXHPNNThx\n4kREOkxE1FPZbDb8+9//xuWXX46//OUv2LFjBzIyMvDKK6/ggQceCKv78ssv449//COSk5Px5z//\nGVu3boXJZMKvfvUrvPfee/LulQCwbNky3HrrrSgoKMDOnTvxu9/9DseOHcOiRYtw+PDhsDQQSZIa\nbagDAM8++yzuuusufPHFF9iwYQNWr16N//7v/27x91itVuzfvx8rV65ERUUFnn32Wbz66qsYM2YM\n/v73v2P16tXtGp/m+kZE1N9Iov6/lbbB9ddfj3nz5mHChAkIBAJ48sknceDAAZw4cQKxsbEAgGee\neQZPP/00cnJyMGzYMKxduxa5ubk4ffp0k0tTERH1Zl9//TWGDh2KmTNn4sMPP+zu7hARUQ/X7gC8\noaqqKlitVrz11lv43ve+ByEEBgwYgF/84hdy3qDb7UZiYiI2bdqEpUuXRqTjREQ9BQNwIiJqj04/\nhFleXo5AICDPfp89exYOhwNz5syR6xgMBkyfPh379+/vbHNERERERL1apwPwBx54AOPGjcOUKVMA\nXFoztuGGFImJic2uJ0tERERE1F90ahWUFStWYP/+/cjNzW3TgzUN69TttEZE1JvFxsbK26vzzzUi\nov6lIxuSdXgGfPny5Xjttdfw4YcfYsiQIXJ53fbEDXeEczgczW5pTERERETUX3QoAH/ggQfk4HvY\nsGFh59LT05GcnIzdu3fLZW63G7m5uZg6dWrnektERERE1Mu1OwVl2bJl2LFjB958801YrVY5r9ts\nNiMqKgqSJOHBBx/E+vXrkZmZiYyMDKxbtw5msxnz589v9ns7Mn1PTcvLywMAZGVldXNP+haOqzI4\nrsrguCqD46ocjq0yOK7K6Gy6YbsD8Oeeew6SJOHaa68NK8/OzsaTTz4JAHj44YdRU1ODZcuWoays\nDJMnT8bu3bsRFRXVqc4SERERUd8XEAE4K0pQ7LLLrxKXHd+bcgcSYwd0d/c6rd0BeCAQaFO91atX\nt3uXNCIiIiLqH2p9HpS4LqLYdQElLkd4sF3ugN/va3TN+OHT+mcATkRERETUGiEEqt0VoaA6PMAu\ndtnhqixp93cWu/rGktYMwImIiIioQwIBP5yV9VNFHCh2XQjOYjvtqKmtjmh7xU4G4ERERETUxwVT\nRUIz2M5L+djFLjtKKi42mSoSCVEGM+KtycFXTPCYlniZIm11NQbgRERERP2YEAJVdakizgsNHnx0\nwFVVqki7kqRCbHScHGDHWVNCAXcS4q3JMOr77uIdDMCJiIiI+rhAwI+yymJ5BrsuuK57745wqkgd\nrUaHeGsy4qzJiLckyTPZ8dZk2CyJ0Ki1irTb0zEAJyIiIuoDar2eRsH1V+dOo8Jdhr8eKIc/oFCq\niNFyKVUkNHsdfKXAEhULSZIUabc3YwBORERE1AsIIVBZU46SckejVJFilx3lVWWKtNt8qkgw4O7L\nqSJKYQBORERE1EM0lSpSf1a7S1JFGrxsloR+myqiFAbgRERERF2ofqpIw1zs0vKLXZQqkhy2wojF\nxFSRrsQAnIiIiCiCWlpVRPFUEXN8WHDtKq5CtCEW06fMglFvUqRdaj8G4ERERETt1N2ritS94lpJ\nFcnLywMABt89DANwIiIioiY0tapIV6SKRButiAtbTYSpIn0NA3AiIiLql3pSqoi8woglmbPV/UCH\nAvC9e/di06ZNOHz4MAoLC7F9+3YsXLhQPr9o0SK89NJLYddMnjwZ+/fv71xviYiIiNqhu1cVCQ+w\ng8v3xZrjuapIP9ehALyqqgqjR4/GwoULsWDBgkb/FCJJEmbPno2XX35ZLtPpdJ3rKREREVETGq4q\nUj9lROlUkWAedhJTRahdOhSA33DDDbjhhhsABGe7GxJCQKfTITExsVOdIyIiImKqCPU1iuSAS5KE\n3NxcJCUlISYmBjNmzMDTTz+NhIQEJZojIiKiXi4gAqj2lOP0uU+bnMnmBjTUlygSgF9//fW45ZZb\nkJ6ejrNnz+Lxxx/HrFmz8MknnzAVhYiIqJ+q9XkurSbibJgyYkdABIBPIt8uN6ChnkYSQojOfIHZ\nbMaWLVuwYMGCZutcuHABgwcPxmuvvYabb75ZLne5XPL7M2fOdKYbRERE1M2EEPD4alDpLkNFo5cT\nNbUVirQrQUKU3oJoQyzMDV7RhljoNHpF2qX+KyMjQ35vtVrbfX2XLEOYkpKCgQMH4ssvv+yK5oiI\niEghdaki9QPr+gG31+9RpF21SiMH1A2D7Ci9FWqVWpF2iZTQJQF4UVERzp8/j5SUlGbrZGVldUVX\n+oW6Xa84ppHFcVUGx1UZHFdl9JdxDaaKXESx64KcHlKXMlJScRF+vzKriug1JiTHD2SqSAT1l3u2\nq9XP4uiIDi9DWJcyEggEUFBQgKNHjyIuLg42mw2rV6/Gj370IyQnJ+Prr7/GypUrkZSUFJZ+QkRE\nRN1DCIHqulVF6l5OO4rLg/nZrsoSRdqVJBVio+MurSRiTZGD7G/yC6HTGBgoUr/QoQD80KFDmDVr\nFoDgiierV6/G6tWrsWjRImzduhXHjh3Dyy+/DKfTiZSUFMyaNQt/+9vfEBUVFdHOExERUdMCIgBn\nRUmDJftCM9pOO2q6fFWRJNgsic2uKuI4V6pIf4h6og4F4DNnzkQgEGj2/K5duzrcISIiImobr68W\nJeWOpnd5LHcolioSZTCHpYdcCrJTYIliqghRa7okB5yIiIg6pspd0WSA3RWpInEN8rDr3hv1/Bdt\nos5gAE5ERNSNAiIAV2VpvVzs8J0eazxVirSrVevCt1CPqb8BTfOpIkTUeQzAiYiIFOb1eVEaesCx\nUaqIywGf36tIu6b6qSINAm1LVCxUkkqRdomoZQzAiYiIIqDaU9l0LrbTDmdlCQQ6te9dkyRIiDHH\ny8F1wxltkz464m0SUecxACciImqD+hvQHDzuahRoV7uV2eVRo9aGB9ZygJ0CmzkRWg1TRYh6Gwbg\nREREIT6/F6XlF8PXxq6XKuL11wYrHotsuyZ9dFh6SP2HH63RNqaKEPUxDMCJiKhfqfFUN8jBviAH\n2mWVJRCi+WV2O0qChJjoOMTFJCPekiTPYNcF2SYDU0WI+hMG4ERE1KcIIVBeXdbs0n1VNeWKtKtW\naxBnuZQqklAvwLZZEqHV6BRpl4h6HwbgRETU6/j9PpRWFDW5bF+Jy4Fan0eRdrVqPcyGWAwaMBTx\n9bZRj7cmIybaBpVKrUi7RNS3MAAnIqIeyVNbg2JXw6X7gsF2WXkRAgqkigCANcrWzC6PyThx7DQk\nSUJWVpYibRNR/8AAnIiIuoUQApU1riYfeCx22VFR7VSkXbVKA5slscm1seOsSdBp9M1eyy3WiSgS\nGIATEZFi/AE/nBXF4TPY9QJtj9etSLt6nbHxsn2hQDs2Op6pIkTUrRiAExFRp9R6PU3v8Oi0o6Ti\nIgIBvyLtWkyxclAd1yDQjjZaOFtNRD1WhwLwvXv3YtOmTTh8+DAKCwuxfft2LFy4MKxOdnY2tm3b\nhrKyMkyaNAlbtmzBiBEjItJpIiLqOkIIVLsrwgNspx1FoXzs8qoyRdpVqdSwmROazMeOsyZDrzUo\n0i4RkdI6FIBXVVVh9OjRWLhwIRYsWNBoluGZZ57B5s2bkZOTg2HDhmHt2rWYPXs2Tp8+jehornVK\nRNTTBEQAzoqSJh94LHHaUVNbrUi7Oq2h+VQRcwLUTBUhoj6oQwH4DTfcgBtuuAEAsGjRorBzQgg8\n++yzWLlyJW6++WYAQE5ODhITE7Fz504sXbq0cz0mIqIO8fq8KC0PX1WkKLSEX0m5A36/T5F2o43W\nJh54DC7hZzZZmSpCRP1OxHPAz549C4fDgTlz5shlBoMB06dPx/79+xmAExEpqNpTiWKnHV8Xn0CF\nuwxfOA/Ks9jOyhIIiIi3KUkqxJrjm5zFjrMkw6g3RbxNIqLeLOIBuN1uBwAkJSWFlScmJqKwsDDS\nzRER9StCCJRXlcnpIQ2X76tyVyjSrlatQ5w1qV5wXX+XxwRo1FpF2iUi6ou6dBWUlv6ZMS8vrwt7\n0j9wTJXBcVUGx/WSQMCPSo8LFe4yVLjLUBk61r38AWVSRXQaI8yGWJgNMaHjpZdRZw7/M9wLVBcH\ncK64EOfQ/yZXeL8qh2OrDI5rZGVkZHTq+ogH4MnJyQAAh8OBgQMHyuUOh0M+R0TU33n9tY0C64qa\n4LHK41IkVQQATDrzpcDaaAsLsnUaripCRNQVIh6Ap6enIzk5Gbt378b48eMBAG63G7m5udi0aVOz\n13Fb38ip+79cjmlkcVyV0VfHNbjLY3mDpfsupY0otsujWoM4SxK0CM5mj8gYXS8fOwlajU6RdvuL\nvnq/9gQcW2VwXJXhcrk6dX2HlyE8c+YMACAQCKCgoABHjx5FXFwc0tLS8OCDD2L9+vXIzMxERkYG\n1q1bB7PZjPnz53eqs0REPUkg4IezsqTRiiLyLo+1NYq0a9CZmtxGPd6ajJjoOKhU6kt/6Y7jX7pE\nRD1NhwLwQ4cOYdasWQCCed2rV6/G6tWrsWjRIvzlL3/Bww8/jJqaGixbtgxlZWWYPHkydu/ejaio\nqIh2nohIaV5fLUrKHZeW66v30GNJxUXFlu5rapfHhNCDj1EGM5fuIyLqxToUgM+cOROBQKDFOnVB\nORFRT1e3dF/DNJFilx2uylJF8rFVkgo2SyJ3eSQi6oe6dBUUIqLu0HDpviKnPSzIrlZo6T6dRh8W\nXMfVSxuxmROgVvOPYCKi/oh/+hNRn+D3+1BaURQ2i11UlzLissPrq1Wk3SijRQ6qE6wpYTPZZlMM\nU0WIiKgRBuBE1Gt4vO6wTWeCr2CwXVZehIBoOTWuIyRIiGm4y2O9TWi4yyMREbUXA3Ai6jGEEKhy\nV4StKFL/ocfy6jJF2lWrNYi3NM7FDu7ymASthrs8EhFR5DAAJ6IuFRABOCsuLd137OsjqHCX4cMv\ndqLYZYe7tlqRdo06E+Lk4DolLMiOibZBpVIr0i4REVFDDMCJKOK8Pi9Kyx0NNqGxo8h1ASXlDuWW\n7ouKDU8VqZcuwqX7iIiop2AATkQdUuOpvpSD3SAv21lRzKX7iIiImsEAnIiaJIRARbWz0Qx23eeq\nmnJF2tVqdJdWFYlJ4dJ9RETU5/BvMqJ+zB/wo6yiqPHKIs4LKC53oNbrVqTdKINZDqprqwXMhlhM\nGDsF8THJsJhimSpCRER9GgNwoj6u1udBicvRaBa7xBncSj0Q8Ee8TQkSYqLj6j30GL6VulEfJdfN\ny8sDAFyWOiLi/SAiIuqJGIAT9QHV7sqwpfvqz2S7qkoVaVOt0iCuLh+73rrY8THJiLMkQavRKdIu\nERFRb8cAnKgXqNtKvajhA4+hgLvaU6lIu3qtITy4rjeLHRMdx6X7iIiIOoABOFEP0XAr9aJQkK30\nVupmo1VOFWm4lXq00cp8bCIioghTJADPzs7G2rVrw8qSk5NRWFioRHNEvYbH65YD6iJnF22lLqkQ\nGx3XaCY7GGinwKAzRrxNIiIiap5iM+CZmZn4+OOP5c9qNf+pmvo+IQSq3RXB2etG+djKbaWuUWsR\nZ01qlCYS3Eo9ERo1t1InIiLqKRQLwNVqNRITE5X6eqJuExABuCpL5KC62BW+skiNglupN57FDn62\nRtugklSKtEtERESRpVgAnp+fj9TUVOj1ekyaNAnr169Henq6Us0RRZTP70Vp+cWwlUW+KjiNCncZ\ndh50wef3KtJu/a3UExoE2yZupU5ERNQnSEKIiO8XvWvXLlRWViIzMxMOhwPr1q3DqVOncPz4cdhs\nNrmey+WS3585cybS3SBqkdfnQYW7rMlXtadcka3UJUiIMlhhNsTCbLCFjpdeTBUhIiLq+TIyMuT3\nVqu13dcrEoA3VF1djfT0dDz66KNYvny5XM4AnJQkhIDbW40Kd6kcWFe6y1DhdqLCXQq3V5lUEbVK\nExZURxtiYQm9j9JbuXQfERFRL9fZALxLliE0mUwYOXIkvvzyy2brZGVldUVX+oW6nQX7w5gGAn6U\nVRbXWxs7fJ1sj0JbqZvqbaUeTBepe58CSxS3Um+P/nS/diWOqzI4rsrh2CqD46qM+pPIHdElAbjb\n7cbJkycxa9asrmiO+hivrxYl5Q55G/WSeg8/lpRfhD/gU6Rda3TcpQDbmgxXSTXMhlhMnzILJkO0\nIm0SERFR36dIAP7QQw9h7ty5SEtLw8WLF/HUU0+hpqYGCxcuVKI56gOqPZVN7vBY7LLDVVmqSD62\nSqVGnDmx0Tbq8dZkxFmToNPow+rXzSIw+CYiIqLOUCQAP3/+PObNm4fi4mIkJCRgypQpOHjwINLS\n0pRojnqBuq3U6zadKaofbLvsqHZXKNKuTqNvfit1czzUzMcmIiKiLqZIAP7KK68o8bXUw9XfSr2o\n3gy20lupRxktl7ZRrzeLHW9NhtkUw3xsIiIi6lG6JAec+g5PbQ2KXQ55Jrt+2khZhUJbqUNCjDm+\n0eYzdS+j3hTxNomIiIiUwgCcwgghUFnjCts+vX6qSEW1U5F21WoN4i3JjWaw42NSYDMnQqvh+thE\nRETUNzAA74f8AT+cFcX1Aut6S/eVO+CprVGkXYPOVC+4Dp/Fjom2cX1sIiIi6hcYgPdRXn8tCou/\nDqWLXJrBLnHaUVJxEYGAX5F2LabY8BnseikjUdxKnYiIiIgBeG8VniriuPSwo9OOC8XfoMZbCRyM\nfLsqlRo2c0KjZfuCS/clQ681RL5RIiIioj6EAXgP5g/4UVZR1CgPu25VEaV2edRpDeEz2PWC7Vhz\nApfuIyIiIuoEBuDdrMZTjZJyO0rkVJFLK4yUlSuzqggAmI1WxDVMFQnlZZtNVqaKEBERESmEAbjC\nAiIAV2VpaObaETaDXVzuQFVNuSLtSpBgsyQ2mypi0BkVaZeIiIiIWsYAPAI8XjdKXA6UlNcPsEPv\nyx3w+32KtKvT6MO3T7ckIT4mBRfOFSFKb8WkiZMUaZeIiIiIOo4BeBsEZ7FLUOxyoLTcgWKXIzib\nHUodUWptbAAwm2LqzVwntSlVpKooT7H+EBEREVHnMAAPqZ+LXVIvyC5xBZftU2oWW63SwGZJrBdc\nJ4Vms4Pv9UwVISIiIupT+k0AXuv1oLTiYijAvojS8osoKXeEjhdR7a5QrG2TwYx4SxLiY4JpInHy\nQ49JiImO4wY0RERERP1InwnAfX4vSsuLGgXWde+VTBNRSSrEWhKCOdihhxzlmWxrEkz6aMXaJiIi\nIqLeRdEAfOvWrfjVr34Fu92OkSNH4tlnn8XVV1/d7u8RQqDaUwlnRTFKK4pQVlEMZ0UxykLvSysu\nwlVZCgGhwK8IijKY5cDaZklCvDVJDrhjzPFcG5uIiIiI2kSxAPy1117Dgw8+iOeeew5XX301tmzZ\nghtuuAEnTpxAWlpao/oXywrhrAwG1aUNAuyyymLUKrTpTB21WoM4cyLi5Icdk0LpIsGjUR+laPtE\nRERE1D8oFoBv3rwZixcvxpIlSwAAv/vd77Br1y4899xzWL9+faP66176uVJdARBME4kxxwcfeLQk\nIc6SKL+3WRJhjbZBJakU7QMRERERkSIBeG1tLQ4fPoyHH344rHzOnDnYv3+/Ek0CAKxRtmBAbU0M\nBdhJcrAdEx0HtbrPpLwTERERUS8lCSEinjhdWFiIgQMHYu/evWE532vXrsXOnTtx6tQpAIDL5Yp0\n00REREREXcZqtbb7GuZcEBERERF1IUUC8Pj4eKjVajgcjrByh8OBlJQUJZokIiIiIuoVFEmK1ul0\nGD9+PHbv3o1bbrlFLn///fdx6623yp87MmVPRERERNSbKfZU4ooVK3DXXXdh4sSJmDp1Kv74xz/C\nbrfj3nvvVapJIiIiIqIeT7EA/Mc//jFKSkqwbt06XLhwAaNGjcJ7773X5BrgRERERET9hSKroBAR\nERERUdO6bRWUrVu3Ij09HUajEVlZWcjNze2urvQZ2dnZUKlUYa8BAwZ0d7d6lb1792Lu3LkYOHAg\nVCoVcnJyGtXJzs5GamoqTCYTrrnmGpw4caIbetr7tDa2ixYtanT/Tp06tZt62zts2LABEyZMgNVq\nRWJiIubOnYvjx483qsd7tv3aMra8Z9tvy5YtGDNmDKxWK6xWK6ZOnYr33nsvrA7v1/ZrbVx5r0bG\nhg0boFKpcP/994eVd+Se7ZYAvG6b+scffxxHjx7F1KlTccMNN+Cbb77pju70KZmZmbDb7fLr888/\n7+4u9SpVVVUYPXo0fvvb38JoNEKSpLDzzzzzDDZv3ow//OEPOHToEBITEzF79mxUVlZ2U497j9bG\nVpIkzJ49O+z+bfgXM4Xbs2cP7rvvPhw4cAAffvghNBoNrrvuOpSVlcl1eM92TFvGlvds+6WlpWHj\nxo04cuQIPvnkE8yaNQs33XQTPv30UwC8XzuqtXHlvdp5Bw8exLZt2zB69Oiwv786fM+KbjBx4kSx\ndOnSsLKMjAyxcuXK7uhOn7F69WpxxRVXdHc3+ozo6GiRk5Mjfw4EAiI5OVmsX79eLqupqRFms1k8\n//zz3dHFXqvh2AohxMKFC8X3v//9bupR31BZWSnUarV45513hBC8ZyOp4dgKwXs2Umw2m/jTn/7E\n+zXC6sZsTX8NAAAgAElEQVRVCN6rneV0OsVll10mPv74YzFz5kxx//33CyE692dsl8+A121TP2fO\nnLBypbep7y/y8/ORmpqKoUOHYt68eTh79mx3d6nPOHv2LBwOR9i9azAYMH36dN67ESBJEnJzc5GU\nlIThw4dj6dKlKCoq6u5u9Srl5eUIBAKIjY0FwHs2khqOLcB7trP8fj9effVVuN1uTJ8+nfdrhDQc\nV4D3amctXboUt956K2bMmAFR79HJztyziq2C0pzi4mL4/X4kJSWFlScmJsJut3d1d/qUyZMnIycn\nB5mZmXA4HFi3bh2mTp2K48ePw2azdXf3er26+7Ope7ewsLA7utSnXH/99bjllluQnp6Os2fP4vHH\nH8esWbPwySefQKfTdXf3eoUHHngA48aNw5QpUwDwno2khmML8J7tqM8//xxTpkyBx+OB0WjE66+/\njuHDh8sBC+/XjmluXAHeq52xbds25OfnY+fOnQAQln7SmT9juzwAJ+Vcf/318vsrrrgCU6ZMQXp6\nOnJycrB8+fJu7Fnf1zCfmdrvtttuk9+PHDkS48ePx+DBg/Huu+/i5ptv7sae9Q4rVqzA/v37kZub\n26b7kfds2zU3trxnOyYzMxOfffYZXC4X3njjDdx+++346KOPWryG92vrmhvXrKws3qsddPr0aTz2\n2GPIzc2FWq0GAAghwmbBm9PaPdvlKSjcpr7rmEwmjBw5El9++WV3d6VPSE5OBoAm7926cxQ5KSkp\nGDhwIO/fNli+fDlee+01fPjhhxgyZIhcznu285ob26bwnm0brVaLoUOHYty4cVi/fj0mT56MLVu2\nyDEA79eOaW5cm8J7tW0OHDiA4uJijBw5ElqtFlqtFnv37sXWrVuh0+kQHx8PoGP3bJcH4PW3qa/v\n/fff55I4EeZ2u3Hy5En+j02EpKenIzk5OezedbvdyM3N5b2rgKKiIpw/f573byseeOABOUAcNmxY\n2Dnes53T0tg2hfdsx/j9fgQCAd6vEVY3rk3hvdo2N998M44dO4ZPP/0Un376KY4ePYqsrCzMmzcP\nR48eRUZGRofvWXV2dna2wv1vxGKxYPXq1RgwYACMRiPWrVuH3NxcbN++HVartau702c89NBDMBgM\nCAQC+OKLL3DfffchPz8fzz//PMe1jaqqqnDixAnY7Xb8+c9/xqhRo2C1WuH1emG1WuH3+/HLX/4S\nw4cPh9/vx4oVK+BwOPCnP/2JeXStaGlsNRoNVq1aBYvFAp/Ph6NHj+Kee+5BIBDAH/7wB45tM5Yt\nW4aXXnoJb7zxBgYOHIjKykpUVlZCkiTodDpIksR7toNaG9uqqiresx3w6KOPyn9PffPNN3j22Wex\nc+dObNy4EZdddhnv1w5qaVyTk5N5r3aQwWBAQkKC/EpMTMRf//pXDB48GAsXLuzcn7HKLdrSsq1b\nt4ohQ4YIvV4vsrKyxL59+7qrK33G7bffLgYMGCB0Op1ITU0VP/rRj8TJkye7u1u9ykcffSQkSRKS\nJAmVSiW/X7x4sVwnOztbpKSkCIPBIGbOnCmOHz/ejT3uPVoa25qaGvHd735XJCYmCp1OJwYPHiwW\nL14svv322+7udo/WcCzrXmvWrAmrx3u2/VobW96zHbNo0SIxePBgodfrRWJiopg9e7bYvXt3WB3e\nr+3X0rjyXo2s+ssQ1unIPcut6ImIiIiIulC3bUVPRERERNQfMQAnIiIiIupCDMCJiIiIiLoQA3Ai\nIiIioi7EAJyIiIiIqAsxACciIiIi6kIMwImIiIiIuhADcCIiIiKiLsQAnIiIiIioCzEAJyIiIiLq\nQgzAiYiIiIi6EANwIiIiIqIuxACciKgXmDlzJlSqrv0je8iQIUhPT+/SNomI+gMG4EREvYQkSV3e\nXsM2s7OzoVKpkJOT06V9ISLqSzTd3QEiIuqZPvzww2bPdfX/DBAR9SUMwImIqEktpZ8IIbqwJ0RE\nfQtTUIiIOungwYNQqVSYO3dus3WysrKgVqtRUFAgl3300UeYO3cuEhISoNfrMWTIECxbtgwOh6PN\nbQsh8MILL2Dy5Mkwm82IiorClVdeic2bN8Pn8zV5zfnz5/Hggw9i2LBhMJlMsNlsyMrKwurVq8Ou\naZgDPnPmTKxduxYAsHjxYqhUKvlVUFCAlStXQqVS4aWXXmqy3VOnTkGlUmHatGlt/n1ERH2RJDiN\nQUTUaSNGjMCXX36J8+fPIyEhIezc8ePHMWrUKMycOVNO63jmmWewcuVKxMXF4Xvf+x6Sk5Px6aef\n4p///CdSU1Nx8OBBpKamyt8xc+ZM7Nu3D36/P+y7FyxYgB07dmDgwIH44Q9/CK1Wi3/84x84c+YM\n5syZg3fffRdqtVqun5eXh+uvvx6lpaWYPn06Jk+eDLfbjZMnT+Ljjz9GUVERLBYLgGAArlKpkJ+f\nDwDIycnBiy++iD179uCmm27C2LFj5e994IEH4HK5MHToUEyaNAn//ve/G43R8uXL8dvf/hY7duzA\n/PnzOzniRES9mCAiok575plnhCRJ4je/+U2jc//1X/8lJEkSOTk5Qggh9uzZIyRJElOnThUulyus\n7ssvvywkSRK33HJLWPmMGTOESqUKK3v11VeFJEli7NixoqKiQi6vra0V1157rZAkSWzatEku93g8\nYsiQIUKlUomXX365UT8dDofw+Xzy58GDB4v09PSwOqtXrw77LQ394Ac/EJIkic8++yysvKamRsTG\nxoqEhARRW1vb5LVERP0FU1CIiCLgrrvuglqtbrQ6iN/vx44dOxAdHY0f/ehHAIDf/va3AIDnn39e\nnm2uc+edd2Ls2LF46623UFlZ2WKbL7zwAgBgw4YNiI6Olsu1Wi1+85vfAAC2bdsml7/99tsoKCjA\njTfeiDvvvLPR9yUmJobNlnfEz3/+cwDB31bf66+/DqfTiUWLFkGr1XaqDSKi3o4PYRIRRUBKSgpm\nz56NXbt24dNPP8WYMWMAALt374bdbseiRYtgMpkAAP/+97+h0Wjw97//HX/7298afZfH44Hf78cX\nX3yBK6+8stk2Dx8+DEmScM011zQ6N2rUKCQkJODMmTOorq6GyWTCwYMHAQA33HBDJH5yk66//noM\nHToUO3bswMaNG+Xf/Pzzz0OlUuGnP/2pYm0TEfUWDMCJiCJk8eLF2LVrF3JycrB582YAkGfEFy1a\nJNcrKSmB3+/HmjVrmv0uSZJQVVXVYnsulwtWqxV6vb7J8ykpKSguLobL5YLJZILT6QSAsNxyJdx7\n7714+OGH8corr2DJkiX4/PPPceDAAVx33XW47LLLFG2biKg3YAoKEVGE/Md//AdiY2Oxc+dOBAIB\nOJ1OvPXWWxg6dCimT58u17NarbBYLAgEAs2+/H5/q6uFWK1WuFwueDyeJs9fuHBBrgcAMTExAIBv\nv/02Ej+3WXfffTcMBoOchlJ3vPfeexVtl4iot2AATkQUITqdDrfffjsuXryId999F6+99ho8Hg8W\nLFgQVm/q1KkoLy/HZ5991qn2xo8fDyEEPvroo0bnjh07hqKiInmpQQCYMmUKAOB///d/O9xmXY54\nw9VY6rPZbLjtttuQl5eH3Nxc7NixAykpKbjppps63C4RUV/CAJyIKIIWL14MIJh6kpOTA5VKhYUL\nF4bVWbFiBQBg6dKlOH/+fKPvcLvdyM3NbbWtJUuWAABWrVoVlq7i9XrlNu655x65/Ac/+AGGDBmC\n9957Dzt27Gj0fQ6Ho8XAGgDi4uIAIGw986bUPYw5b948lJeXY8mSJVCp+FcOERHAdcCJiCJu1KhR\nOH36NHw+X9ja3/X9+te/xiOPPAKdTocbb7wR6enpqKmpwblz57B3714MHToUhw8fluvPnDkTe/fu\nRSAQCPueO++8Ezt37kRaWhpuuukmaLVavP322zhz5gyuu+467Nq1Kyzw/eSTT/Dd735XXgd80qRJ\nqK2txenTp/Gvf/2rxXXAAeD06dMYOXIkoqOjcddddyEpKQkA8Itf/KLRii4TJkzAJ598ArVajbNn\nz2LgwIGdH1wior6gLWsVFhYWigULFoiEhARhMBjEiBEjxJ49e1q85rPPPhPTp08XRqNRpKamirVr\n13Z+0UQiol5g06ZNQpIkoVKpml0vWwghDh48KObNmycGDhwodDqdiI+PF6NHjxb33Xef2LdvX1jd\nmTNnNloHXAghAoGAeP7558XEiRNFVFSUMBqNYuzYsWLTpk3C6/U22e4333wj7rvvPjF06FCh1+tF\nXFycmDBhglizZk3YNUOGDGm0DrgQQrzyyiti/PjxwmQyyb+zoKCgUb0tW7YISZLED37wg2bHgIio\nP2p1BtzpdOLKK6/E9OnTcd999yEhIQH5+flISUlBZmZmk9eUl5dj2LBhmDlzJp588kmcPHkSixcv\nRnZ2tvzPokRE1Lf95Cc/wZ///Ge88847uPHGG7u7O0REPUarAfiqVauwb98+7Nu3r81f+txzz2Hl\nypVwOBzy8lhPP/00nnvuOcWfviciou53/vx5fOc738HAgQNx5syZ7u4OEVGP0uoTMW+++SYmTpyI\n2267DUlJSRg3bhy2bNnS4jUHDhzAtGnTwtamnTNnDgoLC1t9cIeIiHqvnTt3Ys2aNbjuuutQW1uL\ntWvXdneXiIh6nFY34snPz8fWrVuxYsUKrFq1CkeOHMH9998PAFi2bFmT19jtdgwaNCisrO5BHbvd\njsGDBwMIbiJBRER9x3PPPYf9+/djwIABWLNmDW688Ub+WU9EfVrdXgvt0WoAHggEMHHiRDz99NMA\ngDFjxuDMmTPYsmVLswG4JEnt7ggREfV+77zzTnd3gYiox2s1BWXAgAEYMWJEWFlmZibOnTvX7DXJ\nycmw2+1hZQ6HQz5HRERERNRftToDftVVV+HUqVNhZV988QWGDBnS7DVTpkzBI488Ao/HI+eBv//+\n+0hNTZXTTxrqyPQ9NS0vLw8AkJWV1c096Vs4rsrguCqD46oMjmv7CSEQEAH4Az74/X4EAj74A/7g\n54A/+PL7cOzY5wiIAIYNH4aACJYFRAB+f7Be82WXzvlFAIHQdwbqXsIPfyAgf/YLf4M6Abmsyc8i\ngEDo+kvlgVB5vbJQPSECrQ8KddiPr7kXV4++vru70enUulYD8OXLl2Pq1KlYv349fvzjH+PIkSP4\n/e9/jw0bNsh1Vq5ciUOHDuGDDz4AAMyfPx9r1qzBokWL8Pjjj+P06dN45plnkJ2d3anOEhER9Xf+\ngB8+vxd+vw8+vw/+gDd09NUr88Hn99Z774Pf7730Xi7zwRfwISCXhwLZ0Hf5A374Al75fd05X8CH\ngN8PXyC8bl1QHfDXC64Dvnb9vvc+U2jgqE/oK/+D02oAnpWVhTfffBOrVq3CU089hcGDB2PdunX4\n2c9+Jtex2+1hO6VZLBa8//77WLZsGbKysmCz2fDQQw9h+fLlyvwKIiIihfgDfvh8tfD6vajylMMf\n8KGwuEAOcH3+2tD74GevL/xz8+eDQXGwjhc+OWiuC64bvg8Gvn0lACHqCIG+sYF7qwE4ANx4440t\nbqKwffv2RmVXXHEF9uzZ0/GeERERhQgh5GDW6wu+an218Po88Pm9oc8eObgN++zzwuv3hK7zwuuv\nK2t89IYCbV+9Y6CpgPdw148B9R4qSQVJpYJaUkOlUkMlqeRjsPzS+frH+vXkY6hu3XuVSgVJUrW5\nzsWLRZAgYcCA1LB6kiQ1/j4p9D6sjUvlkhR+rSRJUKnUkCDVu0aCJKkb1FFBwqX3l77vUp363938\n+2B7fUGbAnAiIqLm+P0+eHxueL218Hjd8Po88Hg9oWPwc63XIwfEwWPovdcjB9K1Pg+83lrU+oPH\nYNmlQJozv72XJKmgVqlDLw3UKjVUak3YZ4+nFipJBbPZArUUrKsKnVepVKFj/fJLx0vnVFCpNFDX\nBacN60vq0HfVBcYNv09Vr44m1G9VWN2636JSNQyWg+3XBZY9BZ9b6JkYgBMR9QP+gB+1Xjc8Xnej\nY/h7D2q9NfDUulHrc8NT64bH54HX6wkF2ZeC6VqvG7W+2nbn+FLnSJCgUWuhUWugUWuhVmugVmug\nUWlDR029skt16so1ag3UYXVDQbBaC7VKHTofetWdCwXJ6tA5jTr8sxxINwiqVaHvUEmtLrrGQJH6\nlVYD8Ozs7EY7mSUnJ6OwsLDJ+l9//TWGDh3aqHzXrl2YM2dOB7tJRNR/CCHg9dXCXVsDjzf0qq2B\nx+sOlblDn4OBsrveeU9tDdzeGrjKy+D11+LvnwTg8brh83u7+2f1WhIkaDU6aDQ6iACgltSIioqG\nVq0LC4Q1Gp38vv45rUYHtVobKteGysMD6EZlqsZlcrDdR/4Jnqg/a9MMeGZmJj7++GP5s1rd+n/8\n//znPzFmzBj5c2xsbPt7R0TUiwREAO7aarg9weDYXVsNd23o6Kn3vn55gzJP6HOTecf9nEqlhlaj\ng06tg1ajg1ajDwXG2lBZvc9159XBz1qNHjpNMCgOXquTA2WtRgtNk0cdtGptMMc1lFLAWVoiioQ2\nBeBqtRqJiYnt+mKbzdbua4iIuosQIhQMV6HGU4UaTzAorvFUoSYUQNeEPrvDyqrkc+7a6u7+Gd1C\nklTQafXQawzQanXQawzQaQ3QaXTBo1YfCoD10Gl10Kr10Gr1cpAsH+uVaevVDV6v48wvEfUZbQrA\n8/PzkZqaCr1ej0mTJmH9+vVIT09v8Zof/vCHcLvdyMjIwPLly3HLLbdEpMNERM3x+32oqa1GtbsS\nNZ5KVHuCwXTwcxWqPfXK3XWfq+R6ffkhPwlSMEjWGqEPBcV6rRE6nQF6jR56nRE6jR56nQE6rRH6\nuvNaPXQaQ+iolwNqnUYPvdYArUYPjVrTox46IyLq6VoNwCdPnoycnBxkZmbC4XBg3bp1mDp1Ko4f\nPw6bzdaovtlsxq9//WtcddVV0Gg0eOutt3DbbbchJycHd9xxhyI/goj6Fq/Pi2pPBardlah2V6DK\nXYlqdyWq3BWodgfLqzwVqK6pQJWnUq7n8bq7u+sRo1ZrYNAaodcZYQgFynWf9VoDDDqjHEzrdUYY\ndEbotEa5/KszX0Gj1iHryomhQFnHIJmIqIeQhBDtWtG8uroa6enpePTRR9u8sc59992Hffv24dNP\nPw0rr7+N55kzZ9rTDSLqJXx+L9y+ani8wZfbV9P4va8aHp8bHm8Nan018AV67wODGlUonUKtD6Za\naEIpF+pgmUYdSquoVx6sF8xXrvvMdAsiop4rIyNDfm+1Wtt9fbuXITSZTBg5ciS+/PLLNl8zYcIE\n/OUvf2lvU0TUwwgh4PHVwO2tCr2q4Q4F08HyUDDtvfS+Ny1Rp1ZpoFMboNPooW10DKVgaAyhYNoQ\nyl82BIPtUJ22LLdGRET9W7sDcLfbjZMnT2LWrFltvubo0aMYMGBAi3X4RHnk8Cl9ZfTFcRVCwON1\no6LaicqaclTWuFBR7UJltRMVNS5UVrvCjlU15T16dQ4JEgx6E0z6aBj1UTDpo2A01L2/VGYKlRn1\n0cH3uigY9VHQarTd/RMipi/erz0Bx1U5HFtlcFyVUT+LoyNaDcAfeughzJ07F2lpabh48SKeeuop\n1NTUYOHChQCAlStX4tChQ/jggw8AADk5OdDpdBg7dixUKhXefvttbN26FRs3buxUR4mo7by+WpRX\nl6G8yomK0LG8ugwVoWN5tRPlVWWorHbB66/t7u42IkkqmAzRiNJHw2Q0I0pvhskQDJajDGaYDGZE\nGaJhMphh0kcjyhg8GvQmzkATEVGP12oAfv78ecybNw/FxcVISEjAlClTcPDgQaSlpQEA7HY78vPz\n5fqSJGHdunUoKCiAWq3G8OHDsX37dsyfP1+5X0HUDwghUO2phKuyFK6qUpRXBQPpitCxfoBd46nq\n7u7K1CoNoo0WRBktiDaYEWW0ICp0jDZa4LhQDL3GiHGjx8sBtl5nZCBNRER9VqsB+CuvvNLi+e3b\nt4d9XrBgARYsWNC5XhH1M7U+D1yVpSivKoWzshSuqhI50HZVlsJZVYLyyrIeMVtt1JkQbYqB2WhF\ntMmKaKMZ0UYrogwWRBnNwWDbYJGDbr3W0OLqG3m+4D+PDkr6Tlf9BCIiom7V7hxwImqfGk8VyiqK\nUVZRBGdlCcoqikOBdYkcYFd7Krutfxq1FuZ6AfWlwNoKc4NjtNHap/KkiYiIugMDcKJO8Pq8cFYW\nw1lZHAqyQ4F2RTHKQmXdsTuiSlIFg+qoGFhMsbCYYmCJioXZFANLlA0WU0zofWyrM9REREQUWQzA\niVpQ46lGabkDJeUOnCj8BFVuFz61f4CyyhKUVRShotrZpf3RafSwRsfBGhULa5QN5qjw4NoaFQuz\nKRZRRjNzqImIiHqoVgPw7OxsrF27NqwsOTkZhYWFzV7z+eef47777sOhQ4dgs9nw05/+FE888UTn\ne0sUYV5fLUrLL6Kk3IESlwMloffBsouodld0ST9UkgqWUFBtjY5DTLQNligbYqLjQmU2WKNsMOhM\nnK0mIiLq5do0A56ZmYmPP/5Y/qxWN79DW3l5OWbPno2ZM2ciLy8PJ0+exOLFixEVFYUVK1Z0usNE\n7REQAbgqS1DktKPEZW8QYDtQXlWmeB80ai1io+MRa45HjDl4tEbFyUF1THQcoo0WqLjzIRERUb/Q\npgBcrVYjMTGxTV/417/+FW63Gzk5OdDr9RgxYgROnTqFzZs3MwAnRQghUF5VhovOQhQ5L6Co3rHY\nZYfXp9zKIRIkWKJtcoAdfCUgpt7naKOVs9ZEREQka1MAnp+fj9TUVOj1ekyaNAnr169Henp6k3UP\nHDiAadOmQa/Xy2Vz5szBE088gYKCAgwePDgyPad+RQiBimpXWHAtH1121HrdirSrVmtgMycizpKI\ngEeFKIMVozLHyoG2NcoGtZqPUhAREVHbtRo5TJ48GTk5OcjMzITD4cC6deswdepUHD9+HDabrVF9\nu92OQYMGhZUlJSXJ5xiAU0sCAT+KXQ7YS7+BveQc7KXfwl72DYqcF+CprYl4e5KkQkx0HOIsiYiz\nJMFmTQq9T4TNkgRrtE1+mFHezjeT2/kSERFRx0lCCNGeC6qrq5Geno5HH30Uy5cvb3T+u9/9LtLS\n0vDCCy/IZefOncOQIUNw4MABTJo0SS53uVzy+zNnznSk/9RLBQJ+VLjL4KwugqumOHisLoarpgQB\n4Y9oWzq1AWajDWZDDKINMYjWB49mfQxMeivUzL0mIiKidsjIyJDfW63Wdl/f7n87N5lMGDlyJL78\n8ssmzycnJ8Nut4eVORwO+Rz1L/6AHxXuUjiri+oF2cUorylBQAQi1o5WrYPZYIPFaJOPFoMNZqMN\neo2ROdhERETUY7Q7AHe73Th58iRmzZrV5PkpU6bgkUcegcfjkfPA33//faSmpraYfpKVxX/WjxQ5\nVaKLx7TaXYlvi87i26J8fFuUj/NFZ+Eo/TZigbZOa0BCTAoSYlKQGDMACTEpiLemICFmAMwm5R90\n7K5x7es4rsrguCqD46ocjq0yOK7KqJ/F0RGtBuAPPfQQ5s6di7S0NFy8eBFPPfUUampqsHDhQgDA\nypUrcejQIXzwwQcAgPnz52PNmjVYtGgRHn/8cZw+fRrPPPMMsrOzO9VR6jmEEHBWFoeC7bM4X5SP\nby/mo7SiKCLfbzZakRw3CMm2NCTbBiLJloYkWyospljOZBMREVGv12oAfv78ecybNw/FxcVISEjA\nlClTcPDgQaSlpQEIPliZn58v17dYLHj//fexbNkyZGVlwWaz4aGHHmoyX5x6vkDAj4vOQnx7MT8U\nbAdnuKsisEGNNcoWDLLj0sKC7WijJQI9JyIiIuqZWg3AX3nllRbPb9++vVHZFVdcgT179nS8V9Qt\nhBAorbiIs4WncPbCaXxz8SsUFn+NWp+nU98bGx2PJDnIvhRsmwzREeo5ERERUe/BBYz7Ma/Pi2+L\nvsLZC6eCQbf9dKd2hlSp1Ei2pWFgQjpSE9IxMGEoUhOGwKRnoE1ERERUhwF4P1JeVYazF04HA+4L\np3Du4pfw+30d+i6d1oDU+CGhIDsdAxPSkRI3CFqNLsK9JiIiIupbGID3UQERwLdF+XI6ydkLp1BS\n7ujQd0UbrRgYmtEemBgMuBOsyVBx/WwiIiKidmtXAL5hwwY89thjWLZsGX7/+983Wefrr7/G0KFD\nG5Xv2rULc+bM6VgvqVVCCBQWf42TBUdw6Ng+FFV8C99+b7u/R6vRYXBSBtJTMjEkZTgGJX4Hliiu\nPkJEREQUKW0OwA8ePIht27Zh9OjRbQrG/vnPf2LMmDHy59jY2I71kJpVUe3C6XNHcercUZwqOIry\n6vbnb8eaE5Cekon0lOFIT8lEavwQqNX8hxEiIiIipbQp0nK5XLjzzjuxffv2Nq/nbbPZkJiY2Jm+\nUQN+vw9n7adxquAIThYcwbcX8yEg2ny9WqXBwIT0YMA9IBNDkocj1hyvYI+JiIiIqKE2BeBLly7F\nrbfeihkzZkCItgV8P/zhD+F2u5GRkYHly5fjlltu6VRH+6sSlwMnC47g1LkjOP3NZ/DU1rT52mij\nVZ7ZTk/JRFrSZdBp9Ar2loiIiIhaI4lWIupt27bhT3/6Ew4ePAi1Wo1rrrkGo0aNwu9+97sm65eU\nlOCll17CVVddBY1Gg7feegtPP/00cnJycMcdd4TVrb+N55kzZyLwc3o/r78WDlcBzju/QmFZPirc\npW2+VqPSItk6BANiL8OAmHSYDTbmbhMRERFFWEZGhvzearW2+/oWZ8BPnz6Nxx57DLm5uVCrgyte\nCCFanAWPi4sL2/XyyiuvRElJCTZu3NgoAKcgn9+Lb8vOIL/oGArLvkJA+Nt8rS0qGQNihmJA7FAk\nmNOg5sokRERERD1aizPgL774Iu6++245+AYAv98PSZKgVqtRVVUFrVbbaiM5OTn42c9+hurq6rDy\n+vSBSUAAABb+SURBVDPgHfm/h94sEPDjzLfHkHd6Lz798gDctdWtXwTAbLRi+OCxuHzwOAxPGwtL\nVEyjOnl5eQCArKysiPa5v+O4KoPjqgyOqzI4rsrh2CqD46qMzsawLc6A33zzzZg4caL8WQiBxYsX\nY9iwYVi1alWbgm8AOHr0KAYMGNDuzvU1QgicLz6LvFN78cnpvXBVtZ5eolZpkD4gE5cPGofMweOQ\nmjAEKknVBb0lIiIiIiW0GIBbrdZGUb3JZEJsbCxGjBgBAFi5ciUOHTqEDz74AEBwtlun02Hs2LFQ\nqVR4++23sXXrVmzcuFGhn9DzlZZfRN7pvcg7tQf20m9arR9nTcKIweOROXgsMgaOgkFn7IJeEhER\nEVFXaPeCz5IkhT3YZ7fbkZ+fH3Z+3bp1KCgogFqtxvDhw7F9+3bMnz8/Mj3uJarcFTh6Zj/yTu3B\nV4UnWq1vNsVg/LBpyMqcgbTEy/jwJBEREVEf1e4A/KOPPgr7vH379rDPCxYswIIFCzrXq17K66vF\nsbN5yDv1MU58fRj+gK/F+nqtAWO+MwXjh0/HsLTRfICSiIiIqB/glocRUFZRhA/y/geHTn3c6sOU\nKpUalw8ah6zMGRg1dCJ0Wq7LTURERNSfMADvhNLyi3j/0N9x8MS/Wp3tHpI8HFmZMzAu4yqYTf1r\nxRciIiIiuqRdy2ls2LABKpUK999/f4v1Pv/8c8yYMQMmkwkDBw7EU0891alO9jQl5Q68+q8teCrn\n5/j3sX82G3wnxgzADZPn4YmFz2HFbc9g+pgbGXwTERER9XNtngE/ePAgtm3bhtGjR7f4gGB5eTlm\nz56NmTNnIi8vDydPnsTixYsRFRWFFStWRKTT3aXYZcfuQ3/D/zv5EQKBpjfLMRutuHL4NGQNn4FB\nSd/hw5REREREFKZNAbjL5cKdd96J7du3Izs7u8W6f/3rX+F2u5GTkwO9Xo8RI0bg1KlT2Lx5c68N\nwIucF7D70N9w6ORHCIhAk3XiLEmYM+FHmHD5TGjUbVsfnYiIiIj6nzYF4EuXLsWtt96KGTNmtLgN\nPQAcOHAA06ZNg15/6eHCOXPm4IknnkBBQQEGDx7cuR53oYtlhdh96A3kndrTbOAdb03GnAm3YkLm\nDKjVTKknIiIiopa1GjFu27YN+fn52LlzJwC0mlJht9sxaNCgsLKkpCT5XG8IwB1l5/HP//c6Pjm9\nD6KZwDvBmoI5E29FVuYMLh9IRERERG3WYgB++vRpPPbYY8jNzYVaHQwyhRAtzoJ3NOc5Ly+vQ9dF\nkrO6GJ9/sw9fF5+AQNO/0WKwYXTaNAxJGAlVtQpHDh/p4l62XU8Y076I46oMjqsyOK7K4Lgqh2Or\nDI5rZGVkZHTq+hYD8AMHDqC4uBgjR46Uy/x+P/bt24fnn38eVVVV0GrD852Tk5Nht9vDyhwOh3yu\nJ3JWF+GzUODdHKsxDqPSpmFI/AiopHYtHkNEREREJGsxAL/55psxceJE+bMQAosXL8awYcOwatWq\nRsE3AEyZMgWPPPIIPB6PnAf+/vvvIzU1tcX0k6ysrI7+hg4LiAD+lfc/ePfoX5vN8U62peG7E3+M\ncRlToeolqSZ1/5fbHWPal3FclcFxVQbHVRkcV+VwbJXBcVWGy+Xq1PX/v717D2ryXvMA/k2AcKuk\nUA0XUUBPhCo1UrEaqngFbXWol/aI1FE4p+tOCw6FdXekdSpOOVo5W9dOBUdtF6P1gmyd9thSBA+I\nsNCKVpDipbSopWpSbS0WKqjJu390zPrKJRcCMfj9zGQm+b2/N3l45hl5eP3l9/bYgMvlcsjl4n2r\nPTw84O3tjdGjRwMAMjIyUFNTgyNHjgAAEhISsG7dOiQmJmLNmjU4f/48Nm7caHL3lP72e3srdhdv\nRsOFrv9Lxv+J4Zj9zJ8xThnFK95EREREZDMWb9shkUhE67y1Wi2ampqMr728vFBSUoLk5GRERkbC\nx8cHq1atQlpamm0itoHmn77Hf3+ejZ9v6jodC3giCLMnLobqT5PYeBMRERGRzVncgJeVlYle5+Xl\ndZoTHh6O8vJy66PqI4IgoOqbYnxc/gHu6u+Ijnm4DcKfp/8rr3gTERERUZ96ZDauvn2nA/mlW1Fz\n7minY8N9lfjL8/8OHy9F/wdGRERERI+UR6IB/+nGZXz4+UZc/fmHTsemjH0e86ckwcWZd68kIiIi\nor434BvwU41V2HvkfXTcviUal7m4YcnM1zA+NNpOkRERERHRo8jkYuecnByoVCrjjihRUVEoLCzs\ndv7FixchlUo7PYqLi20auCl6/V0cLP8QeYXZnZpvX59A/Nviv7P5JiIiIqJ+Z/IK+LBhw5CdnQ2l\nUgmDwYCdO3di/vz5qKmpgUql6va8w4cPi457e3vbJmIz3PjtOnZ+8Z+4cPVcp2NPj5qCJTNfg6vM\nvd/iISIiIiK6x2QDHhcXJ3qdlZWFrVu34vjx4z024D4+PlAo+v9Ljed/qIOmaBNab4k3SHeSOmNB\n9F8wZexzom0UiYiIiIj6k0VrwPV6PQoKCtDe3o7o6J6XbyxcuBDt7e1QKpVIS0vDokWLehWoKQbB\ngOLjBfjiy/0QIIiOeT82GElz/wPBfqP6NAYiIiIiIlMkgiAIpibV19dDrVajo6MD7u7u2LdvH+bO\nndvl3J9//hm7du3Cs88+C2dnZ3z66af429/+Bo1Gg5dfflk09/7beDY2Nlr9Q7Tf+R3/2/gpLt/4\nvtOxgMdHYvKoF+Dm4mH1+xMRERER3aNUKo3PH7xrvDnMasDv3LmD5uZmtLS0oKCgAO+//z7KysoQ\nGRlp1oekpKSgoqICdXV1onFbNODXf7uM8vMfo63jZqdj44ZPxVOBk7nkhIiIiIhspl8a8AfFxMQg\nMDCwy7tgdkWj0eDVV1/F77//Lhq/vwG3JvhTjVXYVbQJesNd0binuxeWz05HWNA4i99zIDhx4gQA\nmP0HEpmHee0bzGvfYF77BvPad5jbvsG89o3e9rBW7QOu1+thMBjMnl9bW4uAgABrPqpb1369ij3F\n73VqvoP9QpH0/Cp4Dxpi088jIiIiIrIFkw346tWrMW/ePAQGBuK3337D3r17UV5ejqKiIgBARkYG\nampqcOTIEQB/XO2WyWQYN24cpFIpDh06hNzcXGRnZ9ssaL1Bj93Fm3H7bodofOq4eXhh8nI4O/Gu\nlkRERET0cDLZgOt0OixduhRarRZyuRwqlQpFRUWIiYkBAGi1WjQ1NRnnSyQSZGVl4dKlS3ByckJo\naCjy8vKQkJBgs6D/eeIgLl49Lxp7cdq/IFrV9RdDiYiIiIgeFiYbcFPrvB88vmzZMixbtqx3UfWg\n+acmFH61XzQWHjIBU8Y+32efSURERERkKyZvRf8wuXP3NnYf/i8YDHrjmKe7F+JnJnOnEyIiIiJy\nCA7VgB+q+gjaX5pFY/EzXoOX5+N2ioiIiIiIyDIO04B/21yPo6f+IRqb+OQMqP40yU4RERERERFZ\nzmQDnpOTA5VKBblcDrlcjqioKBQWFvZ4Tn19PaZOnQoPDw8EBgbi7bff7lWQtzrasKf4PdGYz6Ah\nWDj1lV69LxERERFRfzP5Jcxhw4YhOzsbSqUSBoMBO3fuxPz581FTUwOVStVp/s2bNxETE4Np06bh\nxIkTOHv2LJKSkuDp6Yn09HSrgvyfoztwo/W68bUEErwcmwp3V95enoiIiIgci8kGPC4uTvQ6KysL\nW7duxfHjx7tswPfs2YP29nZoNBq4urpi9OjROHfuHDZt2mRVA17bWIWac0dFY9OffgHKwHCL34uI\niIiIyN4sWgOu1+uxf/9+tLe3Izo6uss51dXVmDJlClxdXY1jsbGxuHLlCi5dumRRcC1tvyC/dKto\nzP+J4Zirtt2e4kRERERE/UkiCIJgalJ9fT3UajU6Ojrg7u6Offv2Ye7crm96Exsbi+HDh+ODDz4w\njv3www8IDg5GdXU1Jk6caBxvaWkxPm9sbBS9jyAIKD2bj8s3vjOOSSVSPK/6K3w8fc3/CYmIiIiI\nbEipVBqfy+Vyi8836wp4WFgYTp8+jePHjyMlJQXx8fE4ceJEl3NttR93o+5rUfMNAKrhU9l8ExER\nEZFDM7kGHABcXFwwYsQIAEBERARqamqQk5PT5V0y/fz8oNVqRWM6nc54rDuRkZHG59d+vYr9X/1d\ndHyE/5NYHpcCqdTJnJAfaff+OLo/p9R7zGvfYF77BvPaN5jXvsPc9g3mtW/cv4rDGlbtA67X62Ew\nGLo8plarUVFRgY6ODuNYSUkJhg4diqCgINPvbdBjd/Fm3L77/+e7urhh6exUNt9ERERE5PBMNuCr\nV69GZWUlLl68iPr6emRkZKC8vBxLly4FAGRkZGDWrFnG+QkJCfDw8EBiYiIaGhpw8OBBbNy40ewd\nUP554iAuXj0vGlsQ/VcMlnd/9ZyIiIiIyFGYXIKi0+mwdOlSaLVayOVyqFQqFBUVISYmBgCg1WrR\n1NRknO/l5YWSkhIkJycjMjISPj4+WLVqFdLS0kwG0/zT9yj8ar9oLDxkAtRjZnVzBhERERGRYzHZ\ngHe1ztvU8fDwcJSXl1sUyJ27t7H78GYYDHrjmKe7F+JnJtvsi51ERERERPZm1RrwvnCo6iNof2kW\njS2Z+Rq8PB+3U0RERERERLb30DTgR0/9Q/R64pMzMHbkJDtFQ0RERETUN0w24Bs2bMCECRMgl8uh\nUCgQFxeHhoaGHs+5ePEipFJpp0dxcbFZQfkMGoKFU18x7ycgIiIiInIgJhvw8vJypKSkoLq6GqWl\npXB2dsasWbNw48YNk29++PBhaLVa42P69Okmz5FAgpdjU+Hu6mHeT0BERERE5EBMfgmzqKhI9Hr3\n7t2Qy+Woqqrq9nb09/j4+EChUFgU0PSnX4AyMNyic4iIiIiIHIXFa8Bv3rwJg8EAb29vk3MXLlwI\nX19fTJ48GR9//LHJ+f5PDMdcdYKlIREREREROQyLG/DU1FRERERArVZ3O2fQoEF49913UVBQgC++\n+AIzZ87E4sWLsWfPnm7PcZI6Y9nsNLg4yywNiYiIiIjIYUgEQRDMnZyeno4DBw6gsrISwcHBFn1Q\nSkoKKioqUFdXZxxraWmx6D2IiIiIiB4mcrnc4nPMvgKelpaG/Px8lJaWWtx8A8CECRPQ2Nho8XlE\nRERERAOJyS9hAn8sOykoKEBZWRlGjRpl1QfV1tYiICDAqnOJiIiIiAYKkw14cnIyPvroI3zyySeQ\ny+XQarUA/ljn7enpCQDIyMhATU0Njhw5AgDQaDSQyWQYN24cpFIpDh06hNzcXGRnZ4ve25pL9kRE\nREREjsxkA75161ZIJBLMnDlTNJ6ZmYm33noLAKDVatHU1GQ8JpFIkJWVhUuXLsHJyQmhoaHIy8tD\nQgJ3OCEiIiKiR5tFX8IkIiIiIqLesXgbQlvJzc1FSEgI3N3dERkZicrKSnuFMmBkZmZCKpWKHlx3\nb5ljx44hLi4OgYGBkEql0Gg0neZkZmZi6NCh8PDwwPTp03HmzBk7ROp4TOU2MTGxU/1GRUXZKVrH\nsGHDBkyYMAFyuRwKhQJxcXFoaGjoNI81azlzcsuatVxOTg5UKhXkcjnkcjmioqJQWFgomsN6tZyp\nvLJWbWPDhg2QSqVYuXKlaNyamrVLA56fn4/XX38da9asQW1tLaKiovDcc8+hubnZHuEMKGFhYdBq\ntcZHfX29vUNyKG1tbRg7dizee+89uLu7QyKRiI5v3LgRmzZtwpYtW1BTUwOFQoGYmBi0trbaKWLH\nYSq3EokEMTExovp98BcziZWXlyMlJQXV1dUoLS2Fs7MzZs2ahRs3bhjnsGatY05uWbOWGzZsGLKz\ns3Hq1CmcPHkSM2bMwPz5841bFLNerWMqr6zV3vvyyy+xY8cOjB07VvT7y+qaFezgmWeeEVasWCEa\nUyqVQkZGhj3CGTDWrl0rhIeH2zuMAeOxxx4TNBqN8bXBYBD8/PyE9evXG8du3bolDBo0SNi2bZs9\nQnRYD+ZWEARh+fLlwrx58+wU0cDQ2toqODk5CZ999pkgCKxZW3owt4LAmrUVHx8fYfv27axXG7uX\nV0FgrfbWr7/+KowcOVI4evSoMG3aNGHlypWCIPTu39h+vwJ++/ZtfP3114iNjRWNx8bGoqqqqr/D\nGXCampowdOhQjBgxAkuWLMGFCxfsHdKAceHCBeh0OlHturm5ITo6mrVrAxKJBJWVlfD19UVoaChW\nrFiBa9eu2Tssh3Lz5k0YDAZ4e3sDYM3a0oO5BVizvaXX67F//360t7cjOjqa9WojD+YVYK321ooV\nK/DSSy9h6tSpEO776mRvatasfcBt6fr169Dr9fD19RWNKxQK4xaHZJ1JkyZBo9EgLCwMOp0OWVlZ\niIqKQkNDA3x8fOwdnsO7V59d1e6VK1fsEdKAMmfOHCxatAghISG4cOEC1qxZgxkzZuDkyZOQyWT2\nDs8hpKamIiIiAmq1GgBr1pYezC3AmrVWfX091Go1Ojo64O7ujgMHDiA0NNTYsLBerdNdXgHWam/s\n2LEDTU1N2Lt3LwCIlp/05t/Yfm/Aqe/MmTPH+Dw8PBxqtRohISHQaDRIS0uzY2QD34Prmclyixcv\nNj4fM2YMxo8fj6CgIHz++edYsGCBHSNzDOnp6aiqqkJlZaVZ9ciaNV93uWXNWicsLAynT59GS0sL\nCgoKEB8fj7Kysh7PYb2a1l1eIyMjWatWOn/+PN58801UVlbCyckJACAIgugqeHdM1Wy/L0EZPHgw\nnJycoNPpROM6nQ7+/v79Hc6A5uHhgTFjxuC7776zdygDgp+fHwB0Wbv3jpHt+Pv7IzAwkPVrhrS0\nNOTn56O0tBTBwcHGcdZs73WX266wZs3j4uKCESNGICIiAuvXr8ekSZOQk5Nj7AFYr9bpLq9dYa2a\np7q6GtevX8eYMWPg4uICFxcXHDt2DLm5uZDJZBg8eDAA62q23xtwmUyG8ePHo7i4WDReUlLCLXFs\nrL29HWfPnuUfNjYSEhICPz8/Ue22t7ejsrKStdsHrl27hsuXL7N+TUhNTTU2iKNGjRIdY832Tk+5\n7Qpr1jp6vR4Gg4H1amP38toV1qp5FixYgG+++QZ1dXWoq6tDbW0tIiMjsWTJEtTW1kKpVFpds06Z\nmZmZfRx/J15eXli7di0CAgLg7u6OrKwsVFZWIi8vj7en74VVq1bBzc0NBoMB3377LVJSUtDU1IRt\n27Yxr2Zqa2vDmTNnoNVq8eGHH+Kpp56CXC7HnTt3IJfLodfr8c477yA0NBR6vR7p6enQ6XTYvn07\n19GZ0FNunZ2d8cYbb8DLywt3795FbW0tXnnlFRgMBmzZsoW57UZycjJ27dqFgoICBAYGorW1Fa2t\nrZBIJJDJZJBIJKxZK5nKbVtbG2vWCqtXrzb+nmpubsbmzZuxd+9eZGdnY+TIkaxXK/WUVz8/P9aq\nldzc3DBkyBDjQ6FQYM+ePQgKCsLy5ct7929s323a0rPc3FwhODhYcHV1FSIjI4WKigp7hTJgxMfH\nCwEBAYJMJhOGDh0qvPjii8LZs2ftHZZDKSsrEyQSiSCRSASpVGp8npSUZJyTmZkp+Pv7C25ubsK0\nadOEhoYGO0bsOHrK7a1bt4TZs2cLCoVCkMlkQlBQkJCUlCT8+OOP9g77ofZgLu891q1bJ5rHmrWc\nqdyyZq2TmJgoBAUFCa6uroJCoRBiYmKE4uJi0RzWq+V6yitr1bbu34bwHmtqlreiJyIiIiLqR3a7\nFT0RERER0aOIDTgRERERUT9iA05ERERE1I/YgBMRERER9SM24ERERERE/YgNOBERERFRP2IDTkRE\nRETUj9iAExERERH1IzbgRERERET96P8AHFK+m1F6fwcAAAAASUVORK5CYII=\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from __future__ import (absolute_import, division, print_function,\n", + " unicode_literals)\n", + "\n", + "from numpy import array\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from filterpy.hinfinity import HInfinityFilter\n", + "\n", + "dt = 0.1\n", + "f = HInfinityFilter(2, 1, dim_u=1, gamma=.01)\n", + "\n", + "f.F = array([[1., dt],\n", + " [0., 1.]])\n", + "\n", + "f.H = array([[0., 1.]])\n", + "f.G = array([[dt**2 / 2, dt]]).T\n", + "\n", + "f.P = 0.01\n", + "f.W = array([[0.0003, 0.005],\n", + " [0.0050, 0.100]])/ 1000 #process noise\n", + "\n", + "f.V = 0.01\n", + "f.Q = 0.01\n", + "u = 1. #acceleration of 1 f/sec**2\n", + "\n", + "xs = []\n", + "vs = []\n", + "\n", + "for i in range(1,40):\n", + " f.update (5)\n", + " #print(f.x.T)\n", + " xs.append(f.x[0,0])\n", + " vs.append(f.x[1,0])\n", + " f.predict(u=u)\n", + "\n", + "plt.subplot(211)\n", + "plt.plot(xs)\n", + "plt.title('position')\n", + "plt.subplot(212)\n", + "plt.plot(vs) \n", + "plt.title('velocity')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.3" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/17_Ensemble_Kalman_Filters.ipynb b/Appendix_E_Ensemble_Kalman_Filters.ipynb similarity index 100% rename from 17_Ensemble_Kalman_Filters.ipynb rename to Appendix_E_Ensemble_Kalman_Filters.ipynb diff --git a/Appendix_F_FilterPy_Code.ipynb b/Appendix_F_FilterPy_Code.ipynb new file mode 100644 index 0000000..b530a2e --- /dev/null +++ b/Appendix_F_FilterPy_Code.ipynb @@ -0,0 +1,1379 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "[Table of Contents](http://nbviewer.ipython.org/github/rlabbe/Kalman-and-Bayesian-Filters-in-Python/blob/master/table_of_contents.ipynb)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# FilterPy Source" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For your convienence I have loaded several of FilterPy's core algorithms into this appendix. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## KalmanFilter" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "# %load https://raw.githubusercontent.com/rlabbe/filterpy/master/filterpy/kalman/kalman_filter.py\n", + "\n", + "\"\"\"Copyright 2014 Roger R Labbe Jr.\n", + "\n", + "filterpy library.\n", + "http://github.com/rlabbe/filterpy\n", + "\n", + "Documentation at:\n", + "https://filterpy.readthedocs.org\n", + "\n", + "Supporting book at:\n", + "https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", + "\n", + "This is licensed under an MIT license. See the readme.MD file\n", + "for more information.\n", + "\"\"\"\n", + "\n", + "from __future__ import (absolute_import, division, print_function,\n", + " unicode_literals)\n", + "import numpy as np\n", + "import scipy.linalg as linalg\n", + "from numpy import dot, zeros, eye, isscalar\n", + "from filterpy.common import setter, setter_1d, setter_scalar, dot3\n", + "\n", + "\n", + "\n", + "class KalmanFilter(object):\n", + " \"\"\" Implements a Kalman filter. You are responsible for setting the\n", + " various state variables to reasonable values; the defaults will\n", + " not give you a functional filter.\n", + "\n", + " You will have to set the following attributes after constructing this\n", + " object for the filter to perform properly. Please note that there are\n", + " various checks in place to ensure that you have made everything the\n", + " 'correct' size. However, it is possible to provide incorrectly sized\n", + " arrays such that the linear algebra can not perform an operation.\n", + " It can also fail silently - you can end up with matrices of a size that\n", + " allows the linear algebra to work, but are the wrong shape for the problem\n", + " you are trying to solve.\n", + "\n", + " **Attributes**\n", + "\n", + " x : numpy.array(dim_x, 1)\n", + " state estimate vector\n", + "\n", + " P : numpy.array(dim_x, dim_x)\n", + " covariance estimate matrix\n", + "\n", + " R : numpy.array(dim_z, dim_z)\n", + " measurement noise matrix\n", + "\n", + " Q : numpy.array(dim_x, dim_x)\n", + " process noise matrix\n", + "\n", + " F : numpy.array()\n", + " State Transition matrix\n", + "\n", + " H : numpy.array(dim_x, dim_x)\n", + "\n", + "\n", + " You may read the following attributes.\n", + "\n", + " **Readable Attributes**\n", + "\n", + " y : numpy.array\n", + " Residual of the update step.\n", + "\n", + " K : numpy.array(dim_x, dim_z)\n", + " Kalman gain of the update step\n", + "\n", + " S : numpy.array\n", + " Systen uncertaintly projected to measurement space\n", + "\n", + " \"\"\"\n", + "\n", + "\n", + "\n", + " def __init__(self, dim_x, dim_z, dim_u=0):\n", + " \"\"\" Create a Kalman filter. You are responsible for setting the\n", + " various state variables to reasonable values; the defaults below will\n", + " not give you a functional filter.\n", + "\n", + " **Parameters**\n", + "\n", + " dim_x : int\n", + " Number of state variables for the Kalman filter. For example, if\n", + " you are tracking the position and velocity of an object in two\n", + " dimensions, dim_x would be 4.\n", + "\n", + " This is used to set the default size of P, Q, and u\n", + "\n", + " dim_z : int\n", + " Number of of measurement inputs. For example, if the sensor\n", + " provides you with position in (x,y), dim_z would be 2.\n", + "\n", + " dim_u : int (optional)\n", + " size of the control input, if it is being used.\n", + " Default value of 0 indicates it is not used.\n", + " \"\"\"\n", + "\n", + " assert dim_x > 0\n", + " assert dim_z > 0\n", + " assert dim_u >= 0\n", + "\n", + " self.dim_x = dim_x\n", + " self.dim_z = dim_z\n", + " self.dim_u = dim_u\n", + "\n", + " self._x = zeros((dim_x,1)) # state\n", + " self._P = eye(dim_x) # uncertainty covariance\n", + " self._Q = eye(dim_x) # process uncertainty\n", + " self._B = 0 # control transition matrix\n", + " self._F = 0 # state transition matrix\n", + " self.H = 0 # Measurement function\n", + " self.R = eye(dim_z) # state uncertainty\n", + " self._alpha_sq = 1. # fading memory control\n", + "\n", + " # gain and residual are computed during the innovation step. We\n", + " # save them so that in case you want to inspect them for various\n", + " # purposes\n", + " self._K = 0 # kalman gain\n", + " self._y = zeros((dim_z, 1))\n", + " self._S = 0 # system uncertainty in measurement space\n", + "\n", + " # identity matrix. Do not alter this.\n", + " self._I = np.eye(dim_x)\n", + "\n", + "\n", + " def update(self, z, R=None, H=None):\n", + " \"\"\"\n", + " Add a new measurement (z) to the kalman filter. If z is None, nothing\n", + " is changed.\n", + "\n", + " **Parameters**\n", + "\n", + " z : np.array\n", + " measurement for this update.\n", + "\n", + " R : np.array, scalar, or None\n", + " Optionally provide R to override the measurement noise for this\n", + " one call, otherwise self.R will be used.\n", + " \"\"\"\n", + "\n", + " if z is None:\n", + " return\n", + "\n", + " if R is None:\n", + " R = self.R\n", + " elif isscalar(R):\n", + " R = eye(self.dim_z) * R\n", + "\n", + " # rename for readability and a tiny extra bit of speed\n", + " if H is None:\n", + " H = self.H\n", + " P = self._P\n", + " x = self._x\n", + "\n", + " # y = z - Hx\n", + " # error (residual) between measurement and prediction\n", + " self._y = z - dot(H, x)\n", + "\n", + " # S = HPH' + R\n", + " # project system uncertainty into measurement space\n", + " S = dot3(H, P, H.T) + R\n", + "\n", + " # K = PH'inv(S)\n", + " # map system uncertainty into kalman gain\n", + " K = dot3(P, H.T, linalg.inv(S))\n", + "\n", + " # x = x + Ky\n", + " # predict new x with residual scaled by the kalman gain\n", + " self._x = x + dot(K, self._y)\n", + "\n", + " # P = (I-KH)P(I-KH)' + KRK'\n", + " I_KH = self._I - dot(K, H)\n", + " self._P = dot3(I_KH, P, I_KH.T) + dot3(K, R, K.T)\n", + "\n", + " self._S = S\n", + " self._K = K\n", + "\n", + "\n", + " def test_matrix_dimensions(self):\n", + " \"\"\" Performs a series of asserts to check that the size of everything\n", + " is what it should be. This can help you debug problems in your design.\n", + "\n", + " This is only a test; you do not need to use it while filtering.\n", + " However, to use you will want to perform at least one predict() and\n", + " one update() before calling; some bad designs will cause the shapes\n", + " of x and P to change in a silent and bad way. For example, if you\n", + " pass in a badly dimensioned z into update that can cause x to be\n", + " misshapen.\"\"\"\n", + "\n", + " assert self._x.shape == (self.dim_x, 1), \\\n", + " \"Shape of x must be ({},{}), but is {}\".format(\n", + " self.dim_x, 1, self._x.shape)\n", + "\n", + " assert self._P.shape == (self.dim_x, self.dim_x), \\\n", + " \"Shape of P must be ({},{}), but is {}\".format(\n", + " self.dim_x, self.dim_x, self._P.shape)\n", + "\n", + " assert self._Q.shape == (self.dim_x, self.dim_x), \\\n", + " \"Shape of P must be ({},{}), but is {}\".format(\n", + " self.dim_x, self.dim_x, self._P.shape)\n", + "\n", + "\n", + " def predict(self, u=0):\n", + " \"\"\" Predict next position.\n", + "\n", + " **Parameters**\n", + "\n", + " u : np.array\n", + " Optional control vector. If non-zero, it is multiplied by B\n", + " to create the control input into the system.\n", + " \"\"\"\n", + "\n", + " # x = Fx + Bu\n", + " self._x = dot(self._F, self.x) + dot(self._B, u)\n", + "\n", + " # P = FPF' + Q\n", + " self._P = self._alpha_sq * dot3(self._F, self._P, self._F.T) + self._Q\n", + "\n", + "\n", + " def batch_filter(self, zs, Rs=None, update_first=False):\n", + " \"\"\" Batch processes a sequences of measurements.\n", + "\n", + " **Parameters**\n", + "\n", + " zs : list-like\n", + " list of measurements at each time step `self.dt` Missing\n", + " measurements must be represented by 'None'.\n", + "\n", + " Rs : list-like, optional\n", + " optional list of values to use for the measurement error\n", + " covariance; a value of None in any position will cause the filter\n", + " to use `self.R` for that time step.\n", + "\n", + " update_first : bool, optional,\n", + " controls whether the order of operations is update followed by\n", + " predict, or predict followed by update. Default is predict->update.\n", + "\n", + " **Returns**\n", + "\n", + "\n", + " means: np.array((n,dim_x,1))\n", + " array of the state for each time step after the update. Each entry\n", + " is an np.array. In other words `means[k,:]` is the state at step\n", + " `k`.\n", + "\n", + " covariance: np.array((n,dim_x,dim_x))\n", + " array of the covariances for each time step after the update.\n", + " In other words `covariance[k,:,:]` is the covariance at step `k`.\n", + "\n", + " means_predictions: np.array((n,dim_x,1))\n", + " array of the state for each time step after the predictions. Each\n", + " entry is an np.array. In other words `means[k,:]` is the state at\n", + " step `k`.\n", + "\n", + " covariance_predictions: np.array((n,dim_x,dim_x))\n", + " array of the covariances for each time step after the prediction.\n", + " In other words `covariance[k,:,:]` is the covariance at step `k`.\n", + " \"\"\"\n", + "\n", + " try:\n", + " z = zs[0]\n", + " except:\n", + " assert not isscalar(zs), 'zs must be list-like'\n", + "\n", + " if self.dim_z == 1:\n", + " assert isscalar(z) or (z.ndim==1 and len(z) == 1), \\\n", + " 'zs must be a list of scalars or 1D, 1 element arrays'\n", + "\n", + " else:\n", + " assert len(z) == self.dim_z, 'each element in zs must be a'\n", + " '1D array of length {}'.format(self.dim_z)\n", + "\n", + " n = np.size(zs,0)\n", + " if Rs is None:\n", + " Rs = [None]*n\n", + "\n", + " # mean estimates from Kalman Filter\n", + " if self.x.ndim == 1:\n", + " means = zeros((n, self.dim_x))\n", + " means_p = zeros((n, self.dim_x))\n", + " else:\n", + " means = zeros((n, self.dim_x, 1))\n", + " means_p = zeros((n, self.dim_x, 1))\n", + "\n", + " # state covariances from Kalman Filter\n", + " covariances = zeros((n, self.dim_x, self.dim_x))\n", + " covariances_p = zeros((n, self.dim_x, self.dim_x))\n", + "\n", + " if update_first:\n", + " for i, (z, r) in enumerate(zip(zs, Rs)):\n", + " self.update(z, r)\n", + " means[i,:] = self._x\n", + " covariances[i,:,:] = self._P\n", + "\n", + " self.predict()\n", + " means_p[i,:] = self._x\n", + " covariances_p[i,:,:] = self._P\n", + " else:\n", + " for i, (z, r) in enumerate(zip(zs, Rs)):\n", + " self.predict()\n", + " means_p[i,:] = self._x\n", + " covariances_p[i,:,:] = self._P\n", + "\n", + " self.update(z, r)\n", + " means[i,:] = self._x\n", + " covariances[i,:,:] = self._P\n", + "\n", + " return (means, covariances, means_p, covariances_p)\n", + "\n", + "\n", + "\n", + " def rts_smoother(self, Xs, Ps, Qs=None):\n", + " \"\"\" Runs the Rauch-Tung-Striebal Kalman smoother on a set of\n", + " means and covariances computed by a Kalman filter. The usual input\n", + " would come from the output of `KalmanFilter.batch_filter()`.\n", + "\n", + " **Parameters**\n", + "\n", + " Xs : numpy.array\n", + " array of the means (state variable x) of the output of a Kalman\n", + " filter.\n", + "\n", + " Ps : numpy.array\n", + " array of the covariances of the output of a kalman filter.\n", + "\n", + " Q : list-like collection of numpy.array, optional\n", + " Process noise of the Kalman filter at each time step. Optional,\n", + " if not provided the filter's self.Q will be used\n", + "\n", + "\n", + " **Returns**\n", + "\n", + " 'x' : numpy.ndarray\n", + " smoothed means\n", + "\n", + " 'P' : numpy.ndarray\n", + " smoothed state covariances\n", + "\n", + " 'K' : numpy.ndarray\n", + " smoother gain at each step\n", + "\n", + "\n", + " **Example**::\n", + "\n", + " zs = [t + random.randn()*4 for t in range (40)]\n", + "\n", + " (mu, cov, _, _) = kalman.batch_filter(zs)\n", + " (x, P, K) = rts_smoother(mu, cov, fk.F, fk.Q)\n", + "\n", + " \"\"\"\n", + "\n", + " assert len(Xs) == len(Ps)\n", + " shape = Xs.shape\n", + " n = shape[0]\n", + " dim_x = shape[1]\n", + "\n", + " F = self._F\n", + " if not Qs:\n", + " Qs = [self.Q] * n\n", + "\n", + " # smoother gain\n", + " K = zeros((n,dim_x,dim_x))\n", + "\n", + " x, P = Xs.copy(), Ps.copy()\n", + "\n", + " for k in range(n-2,-1,-1):\n", + " P_pred = dot3(F, P[k], F.T) + Qs[k]\n", + "\n", + " K[k] = dot3(P[k], F.T, linalg.inv(P_pred))\n", + " x[k] += dot (K[k], x[k+1] - dot(F, x[k]))\n", + " P[k] += dot3 (K[k], P[k+1] - P_pred, K[k].T)\n", + "\n", + " return (x, P, K)\n", + "\n", + "\n", + " def get_prediction(self, u=0):\n", + " \"\"\" Predicts the next state of the filter and returns it. Does not\n", + " alter the state of the filter.\n", + "\n", + " **Parameters**\n", + "\n", + " u : np.array\n", + " optional control input\n", + "\n", + " **Returns**\n", + "\n", + " (x, P)\n", + " State vector and covariance array of the prediction.\n", + " \"\"\"\n", + "\n", + " x = dot(self._F, self._x) + dot(self._B, u)\n", + " P = self._alpha_sq * dot3(self._F, self._P, self._F.T) + self._Q\n", + " return (x, P)\n", + "\n", + "\n", + " def residual_of(self, z):\n", + " \"\"\" returns the residual for the given measurement (z). Does not alter\n", + " the state of the filter.\n", + " \"\"\"\n", + " return z - dot(self.H, self._x)\n", + "\n", + "\n", + " def measurement_of_state(self, x):\n", + " \"\"\" Helper function that converts a state into a measurement.\n", + "\n", + " **Parameters**\n", + "\n", + " x : np.array\n", + " kalman state vector\n", + "\n", + " **Returns**\n", + "\n", + " z : np.array\n", + " measurement corresponding to the given state\n", + " \"\"\"\n", + "\n", + " return dot(self.H, x)\n", + "\n", + "\n", + " @property\n", + " def alpha(self):\n", + " \"\"\" Fading memory setting. 1.0 gives the normal Kalman filter, and\n", + " values slightly larger than 1.0 (such as 1.02) give a fading\n", + " memory effect - previous measurements have less influence on the\n", + " filter's estimates. This formulation of the Fading memory filter\n", + " (there are many) is due to Dan Simon [1].\n", + "\n", + " ** References **\n", + "\n", + " [1] Dan Simon. \"Optimal State Estimation.\" John Wiley & Sons. \n", + " p. 208-212. (2006)\n", + " \"\"\"\n", + "\n", + " return self._alpha_sq**.5\n", + "\n", + "\n", + " @alpha.setter\n", + " def alpha(self, value):\n", + " assert np.isscalar(value)\n", + " assert value > 0\n", + "\n", + " self._alpha_sq = value**2\n", + "\n", + " @property\n", + " def Q(self):\n", + " \"\"\" Process uncertainty\"\"\"\n", + " return self._Q\n", + "\n", + "\n", + " @Q.setter\n", + " def Q(self, value):\n", + " self._Q = setter_scalar(value, self.dim_x)\n", + "\n", + " @property\n", + " def P(self):\n", + " \"\"\" covariance matrix\"\"\"\n", + " return self._P\n", + "\n", + "\n", + " @P.setter\n", + " def P(self, value):\n", + " self._P = setter_scalar(value, self.dim_x)\n", + "\n", + "\n", + " @property\n", + " def F(self):\n", + " \"\"\" state transition matrix\"\"\"\n", + " return self._F\n", + "\n", + "\n", + " @F.setter\n", + " def F(self, value):\n", + " self._F = setter(value, self.dim_x, self.dim_x)\n", + "\n", + " @property\n", + " def B(self):\n", + " \"\"\" control transition matrix\"\"\"\n", + " return self._B\n", + "\n", + "\n", + " @B.setter\n", + " def B(self, value):\n", + " \"\"\" control transition matrix\"\"\"\n", + " self._B = setter (value, self.dim_x, self.dim_u)\n", + "\n", + "\n", + " @property\n", + " def x(self):\n", + " \"\"\" filter state vector.\"\"\"\n", + " return self._x\n", + "\n", + "\n", + " @x.setter\n", + " def x(self, value):\n", + " self._x = setter_1d(value, self.dim_x)\n", + "\n", + "\n", + " @property\n", + " def K(self):\n", + " \"\"\" Kalman gain \"\"\"\n", + " return self._K\n", + "\n", + "\n", + " @property\n", + " def y(self):\n", + " \"\"\" measurement residual (innovation) \"\"\"\n", + " return self._y\n", + "\n", + " @property\n", + " def S(self):\n", + " \"\"\" system uncertainty in measurement space \"\"\"\n", + " return self._S\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## ExtendedKalmanFilter" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# %load https://raw.githubusercontent.com/rlabbe/filterpy/master/filterpy/kalman/EKF.py\n", + "\n", + "\"\"\"Copyright 2014 Roger R Labbe Jr.\n", + "\n", + "filterpy library.\n", + "http://github.com/rlabbe/filterpy\n", + "\n", + "Documentation at:\n", + "https://filterpy.readthedocs.org\n", + "\n", + "Supporting book at:\n", + "https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", + "\n", + "This is licensed under an MIT license. See the readme.MD file\n", + "for more information.\n", + "\"\"\"\n", + "\n", + "from __future__ import (absolute_import, division, print_function,\n", + " unicode_literals)\n", + "import numpy as np\n", + "import scipy.linalg as linalg\n", + "from numpy import dot, zeros, eye\n", + "from filterpy.common import setter, setter_1d, setter_scalar, dot3\n", + "\n", + "\n", + "class ExtendedKalmanFilter(object):\n", + "\n", + " def __init__(self, dim_x, dim_z, dim_u=0):\n", + " \"\"\" Extended Kalman filter. You are responsible for setting the\n", + " various state variables to reasonable values; the defaults below will\n", + " not give you a functional filter.\n", + "\n", + " **Parameters**\n", + "\n", + " dim_x : int\n", + " Number of state variables for the Kalman filter. For example, if\n", + " you are tracking the position and velocity of an object in two\n", + " dimensions, dim_x would be 4.\n", + "\n", + " This is used to set the default size of P, Q, and u\n", + "\n", + " dim_z : int\n", + " Number of of measurement inputs. For example, if the sensor\n", + " provides you with position in (x,y), dim_z would be 2.\n", + " \"\"\"\n", + "\n", + " self.dim_x = dim_x\n", + " self.dim_z = dim_z\n", + "\n", + " self._x = zeros((dim_x,1)) # state\n", + " self._P = eye(dim_x) # uncertainty covariance\n", + " self._B = 0 # control transition matrix\n", + " self._F = 0 # state transition matrix\n", + " self._R = eye(dim_z) # state uncertainty\n", + " self._Q = eye(dim_x) # process uncertainty\n", + " self._y = zeros((dim_z, 1))\n", + "\n", + " # identity matrix. Do not alter this.\n", + " self._I = np.eye(dim_x)\n", + "\n", + "\n", + " def predict_update(self, z, HJacobian, Hx, u=0):\n", + " \"\"\" Performs the predict/update innovation of the extended Kalman\n", + " filter.\n", + "\n", + " **Parameters**\n", + "\n", + " z : np.array\n", + " measurement for this step.\n", + " If `None`, only predict step is perfomed.\n", + "\n", + " HJacobian : function\n", + " function which computes the Jacobian of the H matrix (measurement\n", + " function). Takes state variable (self.x) as input, returns H.\n", + "\n", + "\n", + " Hx : function\n", + " function which takes a state variable and returns the measurement\n", + " that would correspond to that state.\n", + "\n", + " u : np.array or scalar\n", + " optional control vector input to the filter.\n", + " \"\"\"\n", + " if np.isscalar(z) and self.dim_z == 1:\n", + " z = np.asarray([z], float)\n", + "\n", + " F = self._F\n", + " B = self._B\n", + " P = self._P\n", + " Q = self._Q\n", + " R = self._R\n", + " x = self._x\n", + "\n", + " H = HJacobian(x)\n", + "\n", + " # predict step\n", + " x = dot(F, x) + dot(B, u)\n", + " P = dot3(F, P, F.T) + Q\n", + "\n", + " # update step\n", + " S = dot3(H, P, H.T) + R\n", + " K = dot3(P, H.T, linalg.inv (S))\n", + "\n", + " self._x = x + dot(K, (z - Hx(x)))\n", + "\n", + " I_KH = self._I - dot(K, H)\n", + " self._P = dot3(I_KH, P, I_KH.T) + dot3(K, R, K.T)\n", + "\n", + "\n", + " def update(self, z, HJacobian, Hx, R=None):\n", + " \"\"\" Performs the update innovation of the extended Kalman filter.\n", + "\n", + " **Parameters**\n", + "\n", + " z : np.array\n", + " measurement for this step.\n", + " If `None`, only predict step is perfomed.\n", + "\n", + " HJacobian : function\n", + " function which computes the Jacobian of the H matrix (measurement\n", + " function). Takes state variable (self.x) as input, returns H.\n", + "\n", + "\n", + " Hx : function\n", + " function which takes a state variable and returns the measurement\n", + " that would correspond to that state.\n", + " \"\"\"\n", + "\n", + " P = self._P\n", + " if R is None:\n", + " R = self._R\n", + " elif np.isscalar(R):\n", + " R = eye(self.dim_z) * R\n", + "\n", + " if np.isscalar(z) and self.dim_z == 1:\n", + " z = np.asarray([z], float)\n", + "\n", + " x = self._x\n", + "\n", + " H = HJacobian(x)\n", + "\n", + " S = dot3(H, P, H.T) + R\n", + " K = dot3(P, H.T, linalg.inv (S))\n", + "\n", + " y = z - Hx(x)\n", + " self._x = x + dot(K, y)\n", + "\n", + " I_KH = self._I - dot(K, H)\n", + " self._P = dot3(I_KH, P, I_KH.T) + dot3(K, R, K.T)\n", + "\n", + "\n", + " def predict_x(self, u=0):\n", + " \"\"\" predicts the next state of X. If you need to\n", + " compute the next state yourself, override this function. You would\n", + " need to do this, for example, if the usual Taylor expansion to\n", + " generate F is not providing accurate results for you. \"\"\"\n", + "\n", + " self._x = dot(self._F, self._x) + dot(self._B, u)\n", + "\n", + "\n", + " def predict(self, u=0):\n", + " \"\"\" Predict next position.\n", + "\n", + " **Parameters**\n", + "\n", + " u : np.array\n", + " Optional control vector. If non-zero, it is multiplied by B\n", + " to create the control input into the system.\n", + " \"\"\"\n", + "\n", + " self.predict_x(u)\n", + " self._P = dot3(self._F, self._P, self._F.T) + self._Q\n", + "\n", + "\n", + " @property\n", + " def Q(self):\n", + " \"\"\" Process uncertainty\"\"\"\n", + " return self._Q\n", + "\n", + "\n", + " @Q.setter\n", + " def Q(self, value):\n", + " self._Q = setter_scalar(value, self.dim_x)\n", + "\n", + "\n", + " @property\n", + " def P(self):\n", + " \"\"\" covariance matrix\"\"\"\n", + " return self._P\n", + "\n", + "\n", + " @P.setter\n", + " def P(self, value):\n", + " self._P = setter_scalar(value, self.dim_x)\n", + "\n", + "\n", + " @property\n", + " def R(self):\n", + " \"\"\" measurement uncertainty\"\"\"\n", + " return self._R\n", + "\n", + "\n", + " @R.setter\n", + " def R(self, value):\n", + " self._R = setter_scalar(value, self.dim_z)\n", + "\n", + "\n", + " @property\n", + " def F(self):\n", + " return self._F\n", + "\n", + "\n", + " @F.setter\n", + " def F(self, value):\n", + " self._F = setter(value, self.dim_x, self.dim_x)\n", + "\n", + "\n", + " @property\n", + " def B(self):\n", + " return self._B\n", + "\n", + "\n", + " @B.setter\n", + " def B(self, value):\n", + " \"\"\" control transition matrix\"\"\"\n", + " self._B = setter(value, self.dim_x, self.dim_u)\n", + "\n", + "\n", + " @property\n", + " def x(self):\n", + " return self._x\n", + "\n", + " @x.setter\n", + " def x(self, value):\n", + " self._x = setter_1d(value, self.dim_x)\n", + "\n", + " @property\n", + " def K(self):\n", + " \"\"\" Kalman gain \"\"\"\n", + " return self._K\n", + "\n", + " @property\n", + " def y(self):\n", + " \"\"\" measurement residual (innovation) \"\"\"\n", + " return self._y\n", + "\n", + " @property\n", + " def S(self):\n", + " \"\"\" system uncertainty in measurement space \"\"\"\n", + " return self._S\n", + "\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## UnscentedKalmanFilter" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "# %load https://raw.githubusercontent.com/rlabbe/filterpy/master/filterpy/kalman/UKF.py\n", + "\"\"\"Copyright 2014 Roger R Labbe Jr.\n", + "\n", + "filterpy library.\n", + "http://github.com/rlabbe/filterpy\n", + "\n", + "Documentation at:\n", + "https://filterpy.readthedocs.org\n", + "\n", + "Supporting book at:\n", + "https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", + "\n", + "This is licensed under an MIT license. See the readme.MD file\n", + "for more information.\n", + "\"\"\"\n", + "\n", + "# pylint bug - warns about numpy functions which do in fact exist.\n", + "# pylint: disable=E1101\n", + "\n", + "\n", + "#I like aligning equal signs for readability of math\n", + "# pylint: disable=C0326\n", + "\n", + "from __future__ import (absolute_import, division, print_function,\n", + " unicode_literals)\n", + "\n", + "from numpy.linalg import inv, cholesky\n", + "import numpy as np\n", + "from numpy import asarray, eye, zeros, dot, isscalar, outer\n", + "from filterpy.common import dot3\n", + "\n", + "\n", + "class UnscentedKalmanFilter(object):\n", + " # pylint: disable=too-many-instance-attributes\n", + " # pylint: disable=C0103\n", + " \"\"\" Implements the Unscented Kalman filter (UKF) as defined by Simon J.\n", + " Julier and Jeffery K. Uhlmann [1]. Succintly, the UKF selects a set of\n", + " sigma points and weights inside the covariance matrix of the filter's\n", + " state. These points are transformed through the nonlinear process being\n", + " filtered, and are rebuilt into a mean and covariance by computed the\n", + " weighted mean and expected value of the transformed points. Read the paper;\n", + " it is excellent. My book \"Kalman and Bayesian Filters in Python\" [2]\n", + " explains the algorithm, develops this code, and provides examples of the\n", + " filter in use.\n", + "\n", + "\n", + " You will have to set the following attributes after constructing this\n", + " object for the filter to perform properly.\n", + "\n", + " **Attributes**\n", + "\n", + " x : numpy.array(dim_x)\n", + " state estimate vector\n", + "\n", + " P : numpy.array(dim_x, dim_x)\n", + " covariance estimate matrix\n", + "\n", + " R : numpy.array(dim_z, dim_z)\n", + " measurement noise matrix\n", + "\n", + " Q : numpy.array(dim_x, dim_x)\n", + " process noise matrix\n", + "\n", + "\n", + " You may read the following attributes.\n", + "\n", + " **Readable Attributes**\n", + "\n", + " Pxz : numpy.aray(dim_x, dim_z)\n", + " Cross variance of x and z computed during update() call.\n", + "\n", + "\n", + " **References**\n", + "\n", + " .. [1] Julier, Simon J.; Uhlmann, Jeffrey \"A New Extension of the Kalman\n", + " Filter to Nonlinear Systems\". Proc. SPIE 3068, Signal Processing,\n", + " Sensor Fusion, and Target Recognition VI, 182 (July 28, 1997)\n", + "\n", + " .. [2] Labbe, Roger R. \"Kalman and Bayesian Filters in Python\"\n", + "\n", + " https://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n", + " \"\"\"\n", + "\n", + " def __init__(self, dim_x, dim_z, dt, hx, fx, kappa=0.):\n", + " \"\"\" Create a Kalman filter. You are responsible for setting the\n", + " various state variables to reasonable values; the defaults below will\n", + " not give you a functional filter.\n", + "\n", + " **Parameters**\n", + "\n", + " dim_x : int\n", + " Number of state variables for the filter. For example, if\n", + " you are tracking the position and velocity of an object in two\n", + " dimensions, dim_x would be 4.\n", + "\n", + "\n", + " dim_z : int\n", + " Number of of measurement inputs. For example, if the sensor\n", + " provides you with position in (x,y), dim_z would be 2.\n", + "\n", + " dt : float\n", + " Time between steps in seconds.\n", + "\n", + " hx : function(x)\n", + " Measurement function. Converts state vector x into a measurement\n", + " vector of shape (dim_z).\n", + "\n", + " fx : function(x,dt)\n", + " function that returns the state x transformed by the\n", + " state transistion function. dt is the time step in seconds.\n", + "\n", + " kappa : float, default=0.\n", + " Scaling factor that can reduce high order errors. kappa=0 gives\n", + " the standard unscented filter. According to [1], if you set\n", + " kappa to 3-dim_x for a Gaussian x you will minimize the fourth\n", + " order errors in x and P.\n", + "\n", + " **References**\n", + "\n", + " [1] S. Julier, J. Uhlmann, and H. Durrant-Whyte. \"A new method for\n", + " the nonlinear transformation of means and covariances in filters\n", + " and estimators,\" IEEE Transactions on Automatic Control, 45(3),\n", + " pp. 477-482 (March 2000).\n", + " \"\"\"\n", + "\n", + " self.Q = eye(dim_x)\n", + " self.R = eye(dim_z)\n", + " self.x = zeros(dim_x)\n", + " self.P = eye(dim_x)\n", + " self._dim_x = dim_x\n", + " self._dim_z = dim_z\n", + " self._dt = dt\n", + " self._num_sigmas = 2*dim_x + 1\n", + " self.kappa = kappa\n", + " self.hx = hx\n", + " self.fx = fx\n", + "\n", + " # weights for the sigma points\n", + " self.W = self.weights(dim_x, kappa)\n", + "\n", + " # sigma points transformed through f(x) and h(x)\n", + " # variables for efficiency so we don't recreate every update\n", + " self.sigmas_f = zeros((self._num_sigmas, self._dim_x))\n", + "\n", + "\n", + " def update(self, z, R=None, residual=np.subtract, UT=None):\n", + " \"\"\" Update the UKF with the given measurements. On return,\n", + " self.x and self.P contain the new mean and covariance of the filter.\n", + "\n", + " **Parameters**\n", + "\n", + " z : numpy.array of shape (dim_z)\n", + " measurement vector\n", + "\n", + " R : numpy.array((dim_z, dim_z)), optional\n", + " Measurement noise. If provided, overrides self.R for\n", + " this function call.\n", + "\n", + " residual : function (z, z2), optional\n", + " Optional function that computes the residual (difference) between\n", + " the two measurement vectors. If you do not provide this, then the\n", + " built in minus operator will be used. You will normally want to use\n", + " the built in unless your residual computation is nonlinear (for\n", + " example, if they are angles)\n", + "\n", + " UT : function(sigmas, Wm, Wc, noise_cov), optional\n", + " Optional function to compute the unscented transform for the sigma\n", + " points passed through hx. Typically the default function will\n", + " work, but if for example you are using angles the default method\n", + " of computing means and residuals will not work, and you will have\n", + " to define how to compute it.\n", + " \"\"\"\n", + "\n", + " if isscalar(z):\n", + " dim_z = 1\n", + " else:\n", + " dim_z = len(z)\n", + "\n", + " if R is None:\n", + " R = self.R\n", + " elif np.isscalar(R):\n", + " R = eye(self._dim_z) * R\n", + "\n", + " # rename for readability\n", + " sigmas_f = self.sigmas_f\n", + " sigmas_h = zeros((self._num_sigmas, dim_z))\n", + "\n", + " if UT is None:\n", + " UT = unscented_transform\n", + "\n", + " # transform sigma points into measurement space\n", + " for i in range(self._num_sigmas):\n", + " sigmas_h[i] = self.hx(sigmas_f[i])\n", + "\n", + " # mean and covariance of prediction passed through inscented transform\n", + " zp, Pz = UT(sigmas_h, self.W, self.W, R)\n", + "\n", + " # compute cross variance of the state and the measurements\n", + " '''self.Pxz = zeros((self._dim_x, dim_z))\n", + " for i in range(self._num_sigmas):\n", + " self.Pxz += self.W[i] * np.outer(sigmas_f[i] - self.x,\n", + " residual(sigmas_h[i], zp))'''\n", + "\n", + " # this is the unreadable but fast implementation of the\n", + " # commented out loop above\n", + " yh = sigmas_f - self.x[np.newaxis, :]\n", + " yz = residual(sigmas_h, zp[np.newaxis, :])\n", + " self.Pxz = yh.T.dot(np.diag(self.W)).dot(yz)\n", + "\n", + " K = dot(self.Pxz, inv(Pz)) # Kalman gain\n", + " y = residual(z, zp)\n", + "\n", + " self.x = self.x + dot(K, y)\n", + " self.P = self.P - dot3(K, Pz, K.T)\n", + "\n", + "\n", + "\n", + " def predict(self, dt=None):\n", + " \"\"\" Performs the predict step of the UKF. On return, self.xp and\n", + " self.Pp contain the predicted state (xp) and covariance (Pp). 'p'\n", + " stands for prediction.\n", + "\n", + " **Parameters**\n", + " dt : double, optional\n", + " If specified, the time step to be used for this prediction.\n", + " self._dt is used if this is not provided.\n", + "\n", + " Important: this MUST be called before update() is called for the\n", + " first time.\n", + " \"\"\"\n", + "\n", + " if dt is None:\n", + " dt = self._dt\n", + "\n", + " # calculate sigma points for given mean and covariance\n", + " sigmas = self.sigma_points(self.x, self.P, self.kappa)\n", + "\n", + " for i in range(self._num_sigmas):\n", + " self.sigmas_f[i] = self.fx(sigmas[i], dt)\n", + "\n", + " self.x, self.P = unscented_transform(\n", + " self.sigmas_f, self.W, self.W, self.Q)\n", + "\n", + "\n", + " def batch_filter(self, zs, Rs=None, residual=np.subtract, UT=None):\n", + " \"\"\" Performs the UKF filter over the list of measurement in `zs`.\n", + "\n", + "\n", + " **Parameters**\n", + "\n", + " zs : list-like\n", + " list of measurements at each time step `self._dt` Missing\n", + " measurements must be represented by 'None'.\n", + "\n", + " Rs : list-like, optional\n", + " optional list of values to use for the measurement error\n", + " covariance; a value of None in any position will cause the filter\n", + " to use `self.R` for that time step.\n", + "\n", + " residual : function (z, z2), optional\n", + " Optional function that computes the residual (difference) between\n", + " the two measurement vectors. If you do not provide this, then the\n", + " built in minus operator will be used. You will normally want to use\n", + " the built in unless your residual computation is nonlinear (for\n", + " example, if they are angles)\n", + "\n", + " UT : function(sigmas, Wm, Wc, noise_cov), optional\n", + " Optional function to compute the unscented transform for the sigma\n", + " points passed through hx. Typically the default function will\n", + " work, but if for example you are using angles the default method\n", + " of computing means and residuals will not work, and you will have\n", + " to define how to compute it.\n", + "\n", + " **Returns**\n", + "\n", + " means: np.array((n,dim_x,1))\n", + " array of the state for each time step after the update. Each entry\n", + " is an np.array. In other words `means[k,:]` is the state at step\n", + " `k`.\n", + "\n", + " covariance: np.array((n,dim_x,dim_x))\n", + " array of the covariances for each time step after the update.\n", + " In other words `covariance[k,:,:]` is the covariance at step `k`.\n", + " \n", + " \"\"\"\n", + "\n", + " try:\n", + " z = zs[0]\n", + " except:\n", + " assert not isscalar(zs), 'zs must be list-like'\n", + "\n", + " if self._dim_z == 1:\n", + " assert isscalar(z) or (z.ndim==1 and len(z) == 1), \\\n", + " 'zs must be a list of scalars or 1D, 1 element arrays'\n", + "\n", + " else:\n", + " assert len(z) == self._dim_z, 'each element in zs must be a' \\\n", + " '1D array of length {}'.format(self._dim_z)\n", + "\n", + " n = np.size(zs,0)\n", + " if Rs is None:\n", + " Rs = [None]*n\n", + "\n", + " # mean estimates from Kalman Filter\n", + " if self.x.ndim == 1:\n", + " means = zeros((n, self._dim_x))\n", + " else:\n", + " means = zeros((n, self._dim_x, 1))\n", + "\n", + " # state covariances from Kalman Filter\n", + " covariances = zeros((n, self._dim_x, self._dim_x))\n", + " \n", + " for i, (z, r) in enumerate(zip(zs, Rs)):\n", + " self.predict()\n", + " self.update(z, r)\n", + " means[i,:] = self.x\n", + " covariances[i,:,:] = self.P\n", + " \n", + " return (means, covariances)\n", + "\n", + " \n", + "\n", + " def rts_smoother(self, Xs, Ps, Qs=None, dt=None):\n", + " \"\"\" Runs the Rauch-Tung-Striebal Kalman smoother on a set of\n", + " means and covariances computed by the UKF. The usual input\n", + " would come from the output of `batch_filter()`.\n", + "\n", + " **Parameters**\n", + "\n", + " Xs : numpy.array\n", + " array of the means (state variable x) of the output of a Kalman\n", + " filter.\n", + "\n", + " Ps : numpy.array\n", + " array of the covariances of the output of a kalman filter.\n", + "\n", + " Q : list-like collection of numpy.array, optional\n", + " Process noise of the Kalman filter at each time step. Optional,\n", + " if not provided the filter's self.Q will be used\n", + "\n", + " dt : optional, float or array-like of float\n", + " If provided, specifies the time step of each step of the filter.\n", + " If float, then the same time step is used for all steps. If\n", + " an array, then each element k contains the time at step k.\n", + " Units are seconds.\n", + "\n", + " **Returns**\n", + "\n", + " 'x' : numpy.ndarray\n", + " smoothed means\n", + "\n", + " 'P' : numpy.ndarray\n", + " smoothed state covariances\n", + "\n", + " 'K' : numpy.ndarray\n", + " smoother gain at each step\n", + "\n", + "\n", + " **Example**::\n", + "\n", + " zs = [t + random.randn()*4 for t in range (40)]\n", + "\n", + " (mu, cov, _, _) = kalman.batch_filter(zs)\n", + " (x, P, K) = rts_smoother(mu, cov, fk.F, fk.Q)\n", + "\n", + " \"\"\"\n", + " assert len(Xs) == len(Ps)\n", + " n, dim_x = Xs.shape\n", + "\n", + " if dt is None:\n", + " dt = [self._dt] * n\n", + " elif isscalar(dt):\n", + " dt = [dt] * n\n", + "\n", + " if Qs is None:\n", + " Qs = [self.Q] * n\n", + "\n", + " # smoother gain\n", + " Ks = zeros((n,dim_x,dim_x))\n", + "\n", + " num_sigmas = 2*dim_x + 1\n", + "\n", + " xs, ps = Xs.copy(), Ps.copy()\n", + " sigmas_f = zeros((num_sigmas, dim_x))\n", + "\n", + "\n", + " for k in range(n-2,-1,-1):\n", + " # create sigma points from state estimate, pass through state func\n", + " sigmas = self.sigma_points(xs[k], ps[k], self.kappa)\n", + " for i in range(num_sigmas):\n", + " sigmas_f[i] = self.fx(sigmas[i], dt[k])\n", + "\n", + " # compute backwards prior state and covariance\n", + " xb = dot(self.W, sigmas_f)\n", + " Pb = 0\n", + " x = Xs[k]\n", + " for i in range(num_sigmas):\n", + " y = sigmas_f[i] - x\n", + " Pb += self.W[i] * outer(y, y)\n", + " Pb += Qs[k]\n", + "\n", + " # compute cross variance\n", + " Pxb = 0\n", + " for i in range(num_sigmas):\n", + " z = sigmas[i] - Xs[k]\n", + " y = sigmas_f[i] - xb\n", + " Pxb += self.W[i] * outer(z, y)\n", + "\n", + " # compute gain\n", + " K = dot(Pxb, inv(Pb))\n", + "\n", + " # update the smoothed estimates\n", + " xs[k] += dot (K, xs[k+1] - xb)\n", + " ps[k] += dot3(K, ps[k+1] - Pb, K.T)\n", + " Ks[k] = K\n", + "\n", + " return (xs, ps, Ks)\n", + "\n", + "\n", + " @staticmethod\n", + " def weights(n, kappa):\n", + " \"\"\" Computes the weights for an unscented Kalman filter. See\n", + " __init__() for meaning of parameters.\n", + " \"\"\"\n", + "\n", + " assert n > 0, \"n must be greater than 0, it's value is {}\".format(n)\n", + "\n", + " k = .5 / (n+kappa)\n", + " W = np.full(2*n+1, k)\n", + " W[0] = kappa / (n+kappa)\n", + " return W\n", + "\n", + "\n", + " @staticmethod\n", + " def sigma_points(x, P, kappa):\n", + " \"\"\" Computes the sigma points for an unscented Kalman filter\n", + " given the mean (x) and covariance(P) of the filter.\n", + " kappa is an arbitrary constant. Returns sigma points.\n", + "\n", + " Works with both scalar and array inputs:\n", + " sigma_points (5, 9, 2) # mean 5, covariance 9\n", + " sigma_points ([5, 2], 9*eye(2), 2) # means 5 and 2, covariance 9I\n", + "\n", + " **Parameters**\n", + "\n", + " X An array-like object of the means of length n\n", + " Can be a scalar if 1D.\n", + " examples: 1, [1,2], np.array([1,2])\n", + "\n", + " P : scalar, or np.array\n", + " Covariance of the filter. If scalar, is treated as eye(n)*P.\n", + "\n", + " kappa : float\n", + " Scaling factor.\n", + "\n", + " **Returns**\n", + "\n", + " sigmas : np.array, of size (n, 2n+1)\n", + " 2D array of sigma points. Each column contains all of\n", + " the sigmas for one dimension in the problem space. They\n", + " are ordered as:\n", + "\n", + " .. math::\n", + " sigmas[0] = x \\n\n", + " sigmas[1..n] = x + [\\sqrt{(n+\\kappa)P}]_k \\n\n", + " sigmas[n+1..2n] = x - [\\sqrt{(n+\\kappa)P}]_k\n", + " \"\"\"\n", + "\n", + " if np.isscalar(x):\n", + " x = asarray([x])\n", + " n = np.size(x) # dimension of problem\n", + "\n", + " if np.isscalar(P):\n", + " P = eye(n)*P\n", + "\n", + " sigmas = zeros((2*n+1, n))\n", + "\n", + " # implements U'*U = (n+kappa)*P. Returns lower triangular matrix.\n", + " # Take transpose so we can access with U[i]\n", + " U = cholesky((n+kappa)*P).T\n", + " #U = sqrtm((n+kappa)*P).T\n", + "\n", + " sigmas[0] = x\n", + " sigmas[1:n+1] = x + U\n", + " sigmas[n+1:2*n+2] = x - U\n", + "\n", + " return sigmas\n", + "\n", + "\n", + "def unscented_transform(Sigmas, Wm, Wc, noise_cov):\n", + " \"\"\" Computes unscented transform of a set of sigma points and weights.\n", + " returns the mean and covariance in a tuple.\n", + " \"\"\"\n", + "\n", + " kmax, n = Sigmas.shape\n", + "\n", + " # new mean is just the sum of the sigmas * weight\n", + " x = dot(Wm, Sigmas) # dot = \\Sigma^n_1 (W[k]*Xi[k])\n", + "\n", + " # new covariance is the sum of the outer product of the residuals\n", + " # times the weights\n", + " '''P = zeros((n, n))\n", + " for k in range(kmax):\n", + " y = Sigmas[k] - x\n", + " P += Wc[k] * np.outer(y, y)'''\n", + "\n", + " # this is the fast way to do the commented out code above\n", + " y = Sigmas - x[np.newaxis,:]\n", + " P = y.T.dot(np.diag(Wc)).dot(y)\n", + "\n", + " if noise_cov is not None:\n", + " P += noise_cov\n", + "\n", + " return (x, P)\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.4.3" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/pdf/merge_book.py b/pdf/merge_book.py index 88e9e4d..9a630f4 100644 --- a/pdf/merge_book.py +++ b/pdf/merge_book.py @@ -1,56 +1,55 @@ -from __future__ import print_function -import io -import IPython.nbformat as nbformat -import sys -from formatting import * - - -def merge_notebooks(outfile, filenames): - merged = None - added_appendix = False - for fname in filenames: - with io.open(fname, 'r', encoding='utf-8') as f: - nb = nbformat.read(f, nbformat.NO_CONVERT) - remove_formatting(nb) - if not added_appendix and fname[0:8] == 'Appendix': - remove_links_add_appendix(nb) - added_appendix = True - else: - remove_links(nb) - if merged is None: - merged = nb - else: - merged.cells.extend(nb.cells) - #merged.metadata.name += "_merged" - - outfile.write(nbformat.writes(merged, nbformat.NO_CONVERT)) - - -if __name__ == '__main__': - f = open('book.ipynb', 'w', encoding='utf-8') - '''merge_notebooks( - ['../00_Preface.ipynb', - '../01_g-h_filter.ipynb', - '../Appendix_A_Installation.ipynb'])''' - - merge_notebooks(f, - ['../00_Preface.ipynb', - '../01_g-h_filter.ipynb', - '../02_Discrete_Bayes.ipynb', - '../03_Least_Squares_Filters.ipynb', - '../04_Gaussians.ipynb', - '../05_Kalman_Filters.ipynb', - '../06_Multivariate_Kalman_Filters.ipynb', - '../07_Kalman_Filter_Math.ipynb', - '../08_Designing_Kalman_Filters.ipynb', - '../09_Nonlinear_Filtering.ipynb', - '../10_Unscented_Kalman_Filter.ipynb', - '../11_Extended_Kalman_Filters.ipynb', - '../12_Designing_Nonlinear_Kalman_Filters.ipynb', - '../13_Particle_Filters.ipynb', - '../14_Smoothing.ipynb', - '../15_Adaptive_Filtering.ipynb', - '../16_HInfinity_Filters.ipynb', - '../17_Ensemble_Kalman_Filters.ipynb', - '../Appendix_A_Installation.ipynb', - '../Appendix_B_Symbols_and_Notations.ipynb']) +from __future__ import print_function +import io +import IPython.nbformat as nbformat +import sys +from formatting import * + + +def merge_notebooks(outfile, filenames): + merged = None + added_appendix = False + for fname in filenames: + with io.open(fname, 'r', encoding='utf-8') as f: + nb = nbformat.read(f, nbformat.NO_CONVERT) + remove_formatting(nb) + if not added_appendix and fname[0:8] == 'Appendix': + remove_links_add_appendix(nb) + added_appendix = True + else: + remove_links(nb) + if merged is None: + merged = nb + else: + merged.cells.extend(nb.cells) + #merged.metadata.name += "_merged" + + outfile.write(nbformat.writes(merged, nbformat.NO_CONVERT)) + + +if __name__ == '__main__': + f = open('book.ipynb', 'w', encoding='utf-8') + '''merge_notebooks( + ['../00_Preface.ipynb', + '../01_g-h_filter.ipynb', + '../Appendix_A_Installation.ipynb'])''' + + merge_notebooks(f, + ['../00_Preface.ipynb', + '../01_g-h_filter.ipynb', + '../02_Discrete_Bayes.ipynb', + '../03_Least_Squares_Filters.ipynb', + '../04_Gaussians.ipynb', + '../05_Kalman_Filters.ipynb', + '../06_Multivariate_Kalman_Filters.ipynb', + '../07_Kalman_Filter_Math.ipynb', + '../08_Designing_Kalman_Filters.ipynb', + '../09_Nonlinear_Filtering.ipynb', + '../10_Unscented_Kalman_Filter.ipynb', + '../11_Extended_Kalman_Filters.ipynb', + '../12_Designing_Nonlinear_Kalman_Filters.ipynb', + #'../13_Particle_Filters.ipynb', + '../14_Smoothing.ipynb', + '../15_Adaptive_Filtering.ipynb', + '../Appendix_A_Installation.ipynb', + '../Appendix_B_Symbols_and_Notations.ipynb', + '../Appendix_C_Walking_Through_KF_Code.ipynb']) diff --git a/table_of_contents.ipynb b/table_of_contents.ipynb index 5af597d..6184da8 100644 --- a/table_of_contents.ipynb +++ b/table_of_contents.ipynb @@ -15,7 +15,7 @@ "Motivation behind writing the book. How to download and read the book. Requirements for IPython Notebook and Python. github links.\n", "\n", "\n", - "[**Chapter 1: The g-h Filter ($\\alpha$-$\\beta$ Filter)**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/01_g-h_filter.ipynb)\n", + "[**Chapter 1: The g-h Filter**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/01_g-h_filter.ipynb)\n", "\n", "Intuitive introduction to the g-h filter, also known as the $\\alpha$-$\\beta$ Filter, which is a family of filters that includes the Kalman filter. Once you understand this chapter you will understand the concepts behind the Kalman filter. \n", "\n", @@ -92,25 +92,6 @@ "Kalman filters assume a single process model, but manuevering targets typically need to be described by several different process models. Adaptive filtering uses several techniques to allow the Kalman filter to adapt to the changing behavior of the target.\n", "\n", "\n", - "[**Chapter 16: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/16_HInfinity_Filters.ipynb)\n", - " \n", - "Describes the $H_\\infty$ filter. \n", - "\n", - "*I have code that implements the filter, but no supporting text yet.*\n", - "\n", - "\n", - "[**Chapter 17: Ensemble Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/17_Ensemble_Kalman_Filters.ipynb)\n", - "\n", - "Discusses the ensemble Kalman Filter, which uses a Monte Carlo approach to deal with very large Kalman filter states in nonlinear systems.\n", - "\n", - "\n", - "[**Chapter XX: Numerical Stability**](not implemented)\n", - "\n", - "EKF and UKF are linear approximations of nonlinear problems. Unless programmed carefully, they are not numerically stable. We discuss some common approaches to this problem.\n", - "\n", - "*This chapter is not started. I'm likely to rearrange where this material goes - this is just a placeholder.*\n", - "\n", - "\n", "[**Appendix A: Installation, Python, NumPy, and FilterPy**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/Appendix_A_Installation.ipynb)\n", "\n", "Brief introduction of Python and how it is used in this book. Description of the companion\n", @@ -128,6 +109,23 @@ "\n", "A brief walkthrough of the KalmanFilter class from FilterPy.\n", "\n", + "\n", + "[**Appendix D: H-Infinity Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/16_HInfinity_Filters.ipynb)\n", + " \n", + "Describes the $H_\\infty$ filter. \n", + "\n", + "*I have code that implements the filter, but no supporting text yet.*\n", + "\n", + "\n", + "[**Appendix E: Ensemble Kalman Filters**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/17_Ensemble_Kalman_Filters.ipynb)\n", + "\n", + "Discusses the ensemble Kalman Filter, which uses a Monte Carlo approach to deal with very large Kalman filter states in nonlinear systems.\n", + "\n", + "\n", + "[**Appendix F: FilterPy Source Code**](http://nbviewer.ipython.org/urls/raw.github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python/master/Appendix_D_Filterpy_Code.ipynb)\n", + "\n", + "Listings of important classes from FilterPy that are used in this book.\n", + "\n", "### Github repository\n", "http://github.com/rlabbe/Kalman-and-Bayesian-Filters-in-Python\n" ]