{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Table of Contents](table_of_contents.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Topic 16.  Recursive Least Squares\n",
    "Author: Seth Nielsen - sethmnielsen@gmail.com\n",
    "\n",
    "Forgetting Factor Section: Brian Nelson - brnnelson4@gmail.com\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Introduction\n",
    "\n",
    "Least squares approximation is a method used to approximate the solution to an overdetermined system of equations. Linear least squares can be applied to a set of data to find a model that minimizes the sum of the squared errors between the data and their corresponding modeled values, thus creating a model of best fit. This method can further be applied to create a least-squares filter; a type of filter that, given a sequence of input data, uses least squares to find the coefficients that minimize the error between the filter's output and a desired sequence. \n",
    "\n",
    "The traditional form of a least squares filter is performed in *batch* form, where the minimizing solution is computed on an entire block of data after the data has been collected. This method of calculating a solution can be effective for applications suited to offline computation, but is not particularly useful for systems that require real-time or online parameter estimation for parameters that are unknown in advance or are changing. In this case, an adaptive filter is required. The recursive least squares filter modifies the least squares filter to be adaptive by continuously updating the estimated parameters (coefficients) as new data arrive. Thanks to the matrix inversion lemma (see the following derivation), the recursive least squares filter does not require the calculation of any matrix inverse, which can greatly lower computation requirements and precision errors when filtering large amounts of data. The computation can also have a \"forgetting factor\" added that allows the filter to follow a system with changing parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Explanation of the theory\n",
    "\n",
    "We begin with the batch least squares filter with output $y_t$, minimizing coefficients $h_t$, and input signal $f_t$:\n",
    "\n",
    "\\begin{equation}\n",
    "y_t = \\sum_{i=0}^{m-1} h_{i}\\ f_{t-i}\n",
    "\\end{equation}\n",
    "\n",
    "The coefficients minimizing the least-squares error between filter output $y_t$ and a desired signal $d_t$ satisfy the normal equation (with data matrix $A$)\n",
    "\n",
    "\\begin{equation}\n",
    "Rh = A^H d, \\text{ where } R = A^H A = \\sum_{i=1}^N q_{i}\\ q_{i}^H \\\\\n",
    "\\text{with } q_i = \\begin{bmatrix}f_{i}\\\\f_{i-1}\\\\...\\\\f_{i-m+1}\\end{bmatrix}. \\\\\n",
    "\\end{equation}\n",
    "\n",
    "Let\n",
    "\n",
    "\\begin{equation}\n",
    "p = A^H d = \\sum_{i=1}^N q_{i}d_i\n",
    "\\end{equation}\n",
    "\n",
    "and thus the coefficients can be computed as\n",
    "\n",
    "\\begin{equation}\n",
    "h = R^{-1}A^{H} = R^{-1}p\n",
    "\\end{equation}\n",
    "\n",
    "The algorithm can now be made to be adaptive. From this point forward, all terms with subscript $_t$ are computing using the data only up to time $t$. \n",
    "\n",
    "Now define the Grammian matrix $R_t$ as\n",
    "\n",
    "\\begin{equation}\n",
    "R_t = \\sum_{i=1}^t q_{i}q_{i}^H\n",
    "\\end{equation}\n",
    "\n",
    "and \n",
    "\n",
    "\\begin{equation}\n",
    "p_t = \\sum_{i=1}^t q_{i}d_i = p_{t-1} + q_{t}d_{t} \\\\\n",
    "\\hat{h}_t = R_{t}^{-1}p_t \\\\\n",
    "\\end{equation}\n",
    "\n",
    "where $\\hat{h}_t$ is the estimated filter coefficients at time $t$.\n",
    "\n",
    "The algorithm is now adaptive, but it still needs to be recursive ($R_{t}^{-1}$ is currently calculated at each time step). We can break $R_t$ up to be\n",
    "\n",
    "\\begin{align}\n",
    "R_t &= \\sum_{i=1}^{t-1}{q_i}{q_{i}^H} + {q_t}{q_{t}^H} \\\\\n",
    "    &= R_{t-1} + q_{t}q_{t}^H\n",
    "\\end{align}\n",
    "\n",
    "By the _matrix inversion lemma_ (eq. 4.32),\n",
    "\n",
    "\\begin{equation}\n",
    "R_{t}^{-1} = R_{t-1}^{-1} - \\frac{R_{t-1}^{-1}{q_t}\\:{q_{t}^H}\\:R_{t-1}^{-1}}{1 + {q_{t}^H}\\:{R_{t-1}^{-1}}{q_t}}\n",
    "\\end{equation}\n",
    "\n",
    "Let\n",
    "\n",
    "\\begin{equation}\n",
    "P_t = R_{t}^{-1}\n",
    "\\end{equation}\n",
    "\n",
    "and the _Kalman gain_\n",
    "\n",
    "\\begin{equation}\n",
    "k_t = \\frac{R_{t-1}^{-1}{q_t}}{1 + {q_{t}^H}{R_{t-1}^{-1}}{q_t}} = \\frac{P_{t-1}^{-1}{q_t}}{1 + {q_{t}^H}{P_{t-1}^{-1}}{q_t}}.\n",
    "\\end{equation}\n",
    "\n",
    "Then we arrive at\n",
    "\\begin{equation}\n",
    "P_t = P_{t-1} - {k_t}\\:{q_{t}^H}{P_{t-1}}. \\quad \\text{(eq. 4.38)}\n",
    "\\end{equation}\n",
    "\n",
    "The coefficients for the filter, or the estimated parameters, are computed by\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{h}_t &= P_t p_t = P_t (\\ p_{t-1} + q_t\\ d_t\\ ) \\\\\n",
    "&= P_t\\ p_{t-1} + P_t\\ q_t\\ d_t \\\\\n",
    "&= P_{t-1}\\ p_{t-1} - k_t\\ q_t^H\\ P_{t-1}\\ p_{t-1} + P_t\\ q_t\\ d_t \\ \\Leftarrow\\ \\text{using (eq. 4.38)} \\\\\n",
    "&= \\hat{h}_{t-1} - k_t\\ q_t^H\\ \\hat{h}_{t-1} + k_t\\ d_t \\qquad \\qquad \\quad \\Leftarrow\\ k_t = P_t\\ q_t \\\\\n",
    "&= \\hat{h}_{t-1} + k_t (\\ d_t - q_t^H\\ \\hat{h}_{t-1}\\ ).\n",
    "\\end{align}\n",
    "\n",
    "The quantity multiplied by $k_t$ respresents the filter error\n",
    "\n",
    "\\begin{equation}\n",
    "\\varepsilon_t = d_t - q_t^H\\ \\hat{h}_{t-1},\n",
    "\\end{equation}\n",
    "\n",
    "which allows us to write the update of the filter coefficients in the simple form of\n",
    "\n",
    "\\begin{equation}\n",
    "\\hat{h}_t = \\hat{h}_{t-1} + k_t \\varepsilon_t.\n",
    "\\end{equation}\n",
    "\n",
    "### Initialization\n",
    "\n",
    "To initialize the algorithm, the initial condition $P_0 = R_{0}^{-1}$ is needed. A slight perturbation of some small $\\delta$ to the matrix $R_t$ gives the following:\n",
    "\\begin{align}\n",
    "R_t &= \\sum_{i=1}^t q_{i}q_{i}^H + \\delta I \\\\\n",
    "R_0 &= \\delta I \\\\\n",
    "P_0 &= \\delta^{-1} I\n",
    "\\end{align}\n",
    "\n",
    "Computing $P_0$ in this way allows the recursive algorithm to be free of calculating any matrix inverses.\n",
    "\n",
    "The initial filter coefficients can be assumed to be zero, and thus $\\hat{h}_0 = 0$."
   ]
  },
  {
   "source": [
    "## Adding a Forgetting Factor\n",
    "\n",
    "The derivation above allows us to approximate a system with static coefficients--as the algorithm runs and more data is processed, it converges to a set of coefficients that minimize the error between the output of the system that it is generating, and the actual outputs. Well it is doing this, each input and output pair is given the same amount of weight. The first input and output will be weighted the same as the most recently obtained inputs and outputs. This works well if the coefficents of the system are known to be static, but if the system is changing, then we need a way for the system to adapt to place greater weight on new inputs and \"forget\" the old ones.\n",
    "\n",
    "### Modifying the Equations shown above\n",
    "\n",
    "To add a forgetting factor, we can look at the initial equation for how $R_t$ is calculated. As shown above, $R_t$ is the sum shown below.\n",
    "\n",
    "\\begin{align}\n",
    "R_t &= \\sum_{i=1}^{t-1}{q_i}{q_{i}^H} + {q_t}{q_{t}^H} \\\\\n",
    "    &= R_{t-1} + q_{t}q_{t}^H\n",
    "\\end{align}\n",
    "\n",
    "If we pick a number $\\lambda$, where $0 < \\lambda < 1$, and weight the previously computed values with this, the formula will place less weight on older inputs. This is shown below:\n",
    "\n",
    "\\begin{align}\n",
    "R_t &= \\lambda\\sum_{i=1}^{t-1}{q_i}{q_{i}^H} + {q_t}{q_{t}^H} \\\\\n",
    "    &= \\lambda R_{t-1} + q_{t}q_{t}^H\n",
    "\\end{align}\n",
    "\n",
    "Because this is redone at each step, the older an input gets, the less weight it will have in the current $R_t$ matrix. This ends up causing an exponentially decreasing focus on old inputs, that fits exactly what we are looking; as the system slowly changes with time, we still want to use the old inputs in the current computation, but we want more weight on more recent input/output pairs and less focus on inputs and outputs that happened far in the past.\n",
    "\n",
    "If the derivation shown in the previous section to find the formulas for $P_t$, $k_t$, and $\\hat{h}_t$ are followed using the formula for $R_t$ with a forgetting factor, then these formulas end up being as shown below.\n",
    "\n",
    "\\begin{align}\n",
    "P_t = \\lambda ^{-1}(P_{t - 1} - k_t q_t^H P_{t - 1})\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "k_t = \\frac{P_{t-1}q_t}{\\lambda + q_{t}^H P_{t-1}}\n",
    "\\end{align}\n",
    "\n",
    "\\begin{align}\n",
    "\\hat{h}_t = \\lambda^{-1}\\hat{h}_{t-1}+k_t(y_t - \\lambda^{-1}q_t\\hat{h}_{t-1})\n",
    "\\end{align}\n",
    "\n",
    "$\\lambda$ is generally chosen to be between $0.98$ and $1$. The higher the value of $\\lambda$ is chosen, the more the old input/output pairs will affect the current estimates of the system coefficients.\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simple Numerical Examples\n",
    "\n",
    "We can start by estimating the parameters of a simple system using a recursive least squares filter. Let's choose a system with impulse response $h = [\\ 5,6,7,8,9\\ ]$ . The following filter written in Python will give normally-distributed random values $f_t$ as inputs and attempt to match the system's true output $d_t$ by computing the error $\\varepsilon_t = d_t - y_t$ , where $y_t$ is the filtered output evaluated as $y_t = q_{t}^H \\ \\hat{h}_t$, and $\\hat{h}_t$ is the vector of $m$ estimated parameters that is updated by computing $\\hat{h}_t = \\hat{h}_{t-1} + k_t \\varepsilon_t$ at each timestep. The vector $q_t$ is simply the $m$ most recent input values taken from the signal $f_t$. \n",
    "\n",
    "The variable $m$ can be set at the beginning of the script to choose how many filter coefficients will be used, thus defining the size of $\\hat{h}$ and $q_t$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hhat: [5. 6. 7. 8. 9.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "np.set_printoptions(precision=4)\n",
    "\n",
    "\n",
    "m = 5      # number of parameters\n",
    "n = 10000   # iterations\n",
    "\n",
    "h = np.array([5,6,7,8,9])  # impulse response\n",
    "hhat = np.zeros(m)         # initial estimated parameters\n",
    "f = np.random.randn(n)     # normally-distributed random input\n",
    "fn = np.hstack(([0,0,0,0],f))  # convenience array used to shift through input (f)\n",
    "q = np.zeros(m)  # input data for one time step\n",
    "delta = .0001\n",
    "P = 1/delta * np.eye(m)  # initialize P[0] without calculating inverse of R[0]\n",
    "\n",
    "\n",
    "for i in range(n):\n",
    "    d = fn[i:i+5] @ h  # true output\n",
    "    q = fn[i:i+m]      # q update\n",
    "\n",
    "    k = P @ q / (1 + q.T @ P @ q)  # kalman gain vector\n",
    "    y = q @ hhat  # filter output\n",
    "    e = d - y     # error\n",
    "    \n",
    "    hhat = hhat + k * e  # update of estimated parameters\n",
    "    P = P - k * q @ P    # P update\n",
    "\n",
    "print('hhat:',hhat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Engineering Application - Mass-Spring-Damper System Indentification"
   ]
  },
  {
   "attachments": {
    "figure_1.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAF3CAYAAABT8rn8AAAABmJLR0QA/wD/AP+gvaeTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAAB3RJTUUH4gwOBCMDkt4EHAAAIABJREFUeNrt3Xd8FNX+//H37G42FRJKRHoAEZFmbi7KTymilIAUAQlN4KLXixcBIfgF5VoAUa9YQESxAQIWFAEVLs2CIEWULh1CNQSQFEjf7O78/vCSa0yAIAlkJ6/n4+FDsju7Zz8zZ+Y9Z3Z2xjBN0xQAAPBpNmYBAAAEOgAAINABAACBDgAACHQAAAh0AABAoAMAAAIdAAAQ6AAAEOgAAIBABwAABDoAACDQAQAg0AEAAIEOAAAIdAAAQKADAECgAwAAAh0AABDoAACAQAcAgEAHAAAEOgAAINABAACBDgAAgQ4AAAh0AABAoAMAAAIdAAACHQAAEOgAAIBABwAABHrheDweS9fndrupzUeZpklt9Etqoz4CvbCWLVtm6frGjRtHbT4oMzNT+/fvpzYf9Oabb1Ib9RHo12rjYmUpKSnU5qMj2JycHGrzQefOnaM26iPQAfzGMAwZhkFtPloftVHfteQorRvOHTt2aP369Zb8Ts/hcGjTpk2KjY2V1+ulNh/i8XiUmpqqsLAwavOxQNi+fbtSUlIsdxTCyrVJkt1uV0JCgjwej+x2O4Hui86dO6cePXqoRo0algt1m82mzp07q2bNmtTmY7xer7KzsxUYGEhtPiYhIUFVq1a15I6mlWtLTk7W448/bolReqkNdMMwVKNGDUVERFiyPo/HQ20+yDRNZWVlWTL0rFybJPn5+alatWrU5mNCQ0OtM+BRKcbPg6itJNZm1fqsXBvrHLUR6AAAgEAHAAAEOgAABDoAACDQAQAAgQ4AAAh0AAAIdAAAQKADAAACHQAAEOgAABDoAACAQAcAAAQ6AAAg0AEAINABAACBDgAACHQAAECgAwBAoAMAAAIdAAAQ6AAAgEAHAIBABwAABDoAACDQAQAAgQ4AAIEOAAAIdAAAQKADAAACHQCAUshhlUJM09TEiRO1efNmBQcHq2rVqpo0aRJLGABAoPuSb775RqdPn9bnn3/OUgUAlDqWOeQ+c+ZMPfzww+rXr58eeughpaSksHQBAAS6rzlw4IAGDBigESNGqHnz5urSpQtLFwBAoPuanJwcvfDCC2ratKn69+8vf39/7dix44LT2+12GYZBD0CJYhiGZfullWuD7/dLm83349Ay36HXrl1bFStWzP3b6XQqIyMj92+Px6Nly5YpMzNThmEoLi5OSUlJKlOmjEzTtFwHzc7O1pkzZyy58lm1Nknyer3KycmRv78/tfmYzMxMJSYmWm57YvXakpKSlJycrE8++UQ2m002m00dOnRQUFAQgX6t9OnTR9OmTdPMmTOVkpKi48ePq1mzZnlG5J06dcr9e+3atSpXrpwqVKhgyY1LSkpKnh0cavOd0MvKyvLJjUlpru186Fl1e2Ll2gzDUFhYmHr27Onzo3TLBHq3bt20efNmNWvWTMHBwXrvvfcEAEBpYZlAdzgc+ve//80SBQCUSlwpDgAAAh0AABDoAACAQAcAAAQ6AAAEOgAAINABAACBDgAACHQAAAh0AABAoAMAAAIdAAAQ6AAAEOgAAIBABwAABDoAACDQAQAg0AEAAIEOAAAIdAAAQKADAECgAwAAAh0AABDoAACAQAcAgEAHAAAEOgAAINABAACBDgAAgQ4AAAh0AABAoAMAAAIdAAACHQAAEOgAAIBABwAABDoAAAQ6AAAg0AEAwDXkYBYAwOUxTVP79++XYRiSpJo1a+rkyZPatWuXqlatqkaNGkmSfvjhB507d04tW7ZUQEAAMw4EOgCUJGfPnlXXrl21b98+SdJzzz2nf/3rX7nPz507V88//7z27dsnr9ergIAApaamyuFgk4viwyF3ALhMYWFh2rFjR+7fN998s44dO6b69etLkiZMmKD169dr8+bNkiSXy6WjR48y40CgA0BJ89lnn0mSKleurA4dOkhSbmj37t1bYWFhOnHihCTJMAyVKVOGmQYCHQBKmkWLFkmSKlasKKfTqXPnzikjI0OS1LRpU0nS119/nRvo1113HTMNBDoAlDQrV66UJEVGRsowDO3atSv3ucjISEnS2rVrJUnNmzdnhoFAB4CSJjk5Wenp6ZKk6OhoSdLOnTtzn69WrZpycnL066+/SpI6derETAOBfrkWLVqkatWqadu2bSxdAMXi+PHj8ng8eQJ9y5YtkqRmzZpJkrKzs3XkyBFJ0pgxY7R161ZmHIqV5X5D8cUXX+jOO+/M/X0oABS1xo0byzTNPI9Nnz5d1atXl/Tb79RDQkLyTQMwQi+kJ554QjExMcrOzmbJArhmGFCAQL8CR48e1datW9WxY0eWKkDQURu1lTqWOeT+73//W3Pnzi10B7Tb7XRSlMgNp1X7pWEYiouLk9PptGR9CQkJysrKktfrtVxtZ8+eVbVq1Sy9ztlsvj++tUSg//rrr1q2bJnKly8vh8Oh+Ph4TZs2TYMHD1ZUVJQkyePxaNmyZcrMzMzdsCQlJalMmTKW+57LMAxlZ2frzJkzllz5rFqbJHm9XuXk5Mjf399SddlsNm3YsEGdOj0m6RZJVvxu2ZDk/e//rcSU9KkOHjyo0NBQyy21pKQkJScn65NPPpHNZpPNZlOHDh0UFBREoF8LZcqU0QcffCC73a6AgAD9+OOPatmypWrUqJFnRP77n46sXbtW5cqVU4UKFSwZDCkpKapYsSK1+WCgZ2Vl+eTG5FKCggIlDZY0QvA1W1WmTBlLrneGYSgsLEw9e/b0+VG6JQI9ICAgz4UbTNNUw4YNFR4eznoIlLjRHnxsN5PlRqBfo82FaWr58uUsWQBAqWK5C8twohsAgEAHAAAEOgAAINABAACBDgBA6eVgFqA4bNiwQfHx8bLb7QoNDdVdd93FTAEARujwNRs3btSjjz6q7t27a86cOcwQACDQ4YtGjBihkJAQSVKXLl2YIQBAoMNXHTx4UJJyr+KXlJSkxMREZgwAEOjwFXv37pXX61XZsmUVEBAgt9utyMhIValShZkDAAQ6fMWXX34pSQoNDVVISIj27t2rqKgo7dy5k5kDAMWAs9xRLL744gtJUnh4uL766isNHjxYR44cYcYAACN0+JJNmzZJkg4fPqzo6GgdPXpUK1euZMYAAIEOX3Hy5El5vV5J0qRJk9SnTx9JUkxMjFwuFzMIAAh0+IKEhITcQH/wwQf10UcfqXz58jp79qw++eQTZhAAEOjwBXFxcfJ6vQoODs69ne24ceMkSQMGDNDcuXOZSQBQxDgpDkWuW7duysjIkM32v/3FRx55RH//+98lSXa7nZkEAAQ6Sjq73a7AwMA8j9lstnyPAQCKDofcAQAg0AEAAIEOAAAIdAAAQKADAECg4+oZtXKUPKbnqrebmJGo8avHX5Oad57eqYV7Fl6TtmdsnaGkzCQ6HgACHUUjNTtVjac31pL9S/TZrs+uevtTf5yq8d+N10/xP13Vdr2mV61nt9YT3zxx1WvO8eRo1MpRaju3LR0QAIGOohkdV51cVRPvmqi1g9aq36J+SnWlXrX2T6Se0Htb3tPmwZs1YsWIq1r7kv1L1L5Oe/Vt2FePffXYVW27w0cd9Fr71+T2urXpxCY6IgACHX/ertO7FPVOlOb1mKcu9booPDhcE+6coCe+vjojVlOmun3STa+2e1WR10cqyBGk749+f9XqH7ZsmN7q9JZG3zFa0zZO08m0k1el3T2/7tH+M/s18JaBWtxnsR5e8jCdEQCBjj9nS8IWNXm7iRb3WayOdTvmPj62xVgtPbhUbo+72D/DsZRjsht29WrYS5L0cruXFf1htNze4m87dmWs7q51t0KcIQr0C9T8nvP1/PfPX5V5P3LFSH3V/ytJUo3QGgpxhmjGlhl0SgAEOi7P7G2z1eXjLjo8/LAaVWqU7/lxrcap5fsti/1z3DX7Lk3rOC337ybXN9HAJgP1ya7ivVOay+PSxz9/rJldZ+Y+1rFuRy3cs1B7z+wt1raX7F+ihLQE3VjxxtzHPu/9uR7/5vGrsiMDAAS6RXyw4wNNWDNBB4YdUPXQ6gVOM6DJAJ1KO6UtCVuKdaeiUkgl/aXyX/I8PiV6ih5bWbzfZ49YPkKvtHslz2N2m11vdXpLvT7rJVNmsbTrNb16eMnD+vS+T2XIyH08LCBMQ5sO1cytM+mgAAh0XNror0brxXUvaveQ3Qr0u/gNTL4d+K2GLh1abJ9l/JrxWvO3Nfked9qd6t2wtx744oFiaXf7qe2avX22YhrE5Huu042dVDm4svaf2V8sbS/dv1Td63dXvYr18j03/LbhGrJ0iM5knKGjAiDQkd/50eb9C+/XvsR9+vmfP8vf4X/J19UIqyGv6dUnO4v+8Pfz3z+ve264Rw57wTfhm9x+sr4+9LXSXGlF3vZjKx/T8n7L5bAV3Par0a+qzdw2Rd5uZk6m7pt/n8a2GFvg8+UCy2lq9FSN/mo0nRYAgY78DBnq/HFnBfsF64veX1zW6xb2WqgRK0bI5XEV2ec5nX5aT3/3tF5o88JFp5vWcZoGLxlcpPNi96+7lZmTqRY1W1xwmpvDb1a9ivX0zpZ3irTtCWsmaNhtw3R9yPUXnGZI0yFad3ydMlwZdFwABDryqv1abbWq2Upvd377sl9bpUwVDbplkObtnFdkn2fSukn6sPuHCnGGXHS6Tjd20tIDS/Vj/I9F0q7b61aLWS00uf3kS067tO9STVo3SUF+QUV2lOTTXZ/qpbYvXXLaia0nqtmMZnRcAAQ6fnMq7ZSqvlpVz7Z+Vo/d/udPMhvZbKQe+OIBpbvSr/gznUw7qTnb5+i++vddusMYNn3d/+siOwS9eN9itanVRk2rNr3ktE67U/fVv0/Pbnq2SNru9GEn/d/t/1eoaXs26ClJ2nB8A50YAIFe2iVmJKrWa7W0tO9S9Wvc74reKzw4XM/d9dwVXx7VNE11/birpnWcJrvNXqjXRFWJUoBfgNYeXXvF82TI0iF6t8u7hZ5+4l0T9e7OdxWfGn9F7e45s0ebT27Ww38t/MVjVvZfqYf/w8VmABDoJU6wX/BVa+v7o9+rztQ6WjNojZpc36RI3nNM8zFaEbfiik5SO5B0QE6Hs8Czyy/m5ba/XWzmSr7Hj10Rq843dlZZ/7KFfo3D5tBrzV/TC9+/cEXzbuTykVrZb+VlvaZSSCVVCKygd7e8W3wrpM2mIGeQJde3IHuQABQfR6kNc2ew7ph5h1Su+NuyG3YlZiQqPjZewc6i3YkY22Ks2s1tp/UPrv9Tr4/+MFqf9/r8sl/X8LqG6tOwjxbtWZR7RbnLkeZK0wc7PtDp/zt92a+9q/pdem7xcxr1/0apVrlal/36ZQeWKSEtQQ0rNbys1xkytKjXItV9va4GNB5QqF8lXO48aTWrlVKyUmSzWWtf25ChjPR0qeJYiV8AAgR6UcrMydSMLjNUPaK6TNMs1rZM01SFoApFHuaSNLDJQL3w/QvacWqHGldqfFmvnb1ttqqXrX7ZrztvcvRk1X297p8K9JErRmpK9JQ/1a7T5tTMrjN17yf3auvgrbIZhQ8/j9ejB798UOseWHdZrzsvNCBUI5qN0Ltb3tXQW4vuegBnMs6owRsNNKvLLEWUjVBAQICl1je73a61a7/X/WdIc4BAL2Je06vqodVVK6yWz9eytN9SDfp8kFYPWn1Zr3vu++f08z9//tPthjhDFHNzjP6x5B96p1Phf0628/ROzds5T9M6TPtzO0gy1b5Oe7207iXt+XWPGlzXoNCvXRG3Ql3rdf1TI/vzHv7rwwp/KVx9G/VV+cDyV7z8tp3cpuYzm+s//f6jltVbKjMrU0FB1js8fSgwTtKvbHWBYsJJcRZQq1wtZXmy9PHOjwv9molrJqpj3Y5XfNj4tQ6vacm+JcrMySz0a0YsH6Fv+n9zxW1PiZ6i1rNbF3r6HE+Ous7rqidbPnlF7ZYPLK83O75ZJGf6rzi4Qi1mtVDc8Di1qtmq2C5vC4BAhw8wZGhBzAKNXD5SXtN7yenPZJzRuNXjNKH1hCJp//WOrxf6krDbTm6Ty+PSrdVuveJ2G17XUFFVojT9p+mFmv6pVU8ptlmsqpatesVtD/7rYK3/Zb0SMxL/9HvM2jpLsStjtfeRvaoUUomODIBAh1StbDUNaDJA83fPv+S0k3+YrNldZ1/W2eUX0+2mblp2cJk2xm+86HSmaerO9+/UpLaTiqzuZf2W6cX1LxZq2o93fqwX275YZG2PbzVed86+80+99qV1L+mpVU9p8z82F8kOBgAQ6BYy+o7R6r+wv1JdqRec5vjZ45q1dZZ6N+xddJ3IsGntA2sVuzz2otN9vvdzdb6xs5pVK9orrvVu0PuSd4JrN7ednmj+RJG227NBTzntTv0U/9Nlve7+hfdrzbE1OjLiiAIcAXRcAAQ68qoYVFFPt3paT616qsAzuE2Zum/+fXq1/auFvohMYTW8rqEC/AK04ZcLX0nt0RWP6p0u7xR53c+0ekZv/PSGElITCnx+1+ld2nl652VdRKawFvdZrEeWPlKoaT2mRw988YCy3Fla3GfxBW9EAwAEOvRkyye19MBSuU13gaPz87dBLQ6T209WmzltCrzYzPBlw9W+dnsFOgKLvN1Av0AtjFmoCWsKPidgxPIR+m7gd8VSc5UyVRTsDNY7my+9o9J4emNVDKqo+THz6agACPQLiY+P14cffqjXX39dq1evLtUL9ckWT6rfV/kvLdt2Tlu93uH1Ymu3caXG6t+4vz7d9Wmex7PcWVqwZ8FlXeL1crW/ob2+3PelDiQeyPP4soPLdCr9lOpWqFtsbX/W8zM9+e2T8ng9BT7/a/qvuvmNmzU4arAmtZ0kQwZbHgAE+oWMGjVKv/zyi2rXrq2nnnpKH3zwQaldqAOaDFB8ary2ntya+9ic7XMUFhimW66/pVjbnhI9RSNXjMzzWOyKWL3a7tXi7ciGTW/f87bum39f7oWC3F63Bi8erAUxC2QYxReiFYIqaPhtwzVz28x8z6W50lT/jfp6pd0rGn7bcLY4AAj0S5kzZ47GjBmje+65R4sXL9bMmTNL9YKdfffsPN/tjl89vtgOO/9egCNAAxoP0INfPChJ2nFqh2Zvn63uN3cv9rY71eukKmWqaNevuyT9die3bvW7Fevo/Lxhtw7TP5f8U6fT/3cp252ndyr8pXAt6btEHep2YGsDgEAvDKfTmfvvc+fOqWLFiqV6wUaUjZApU4v2LNKL617UPXXvUaBf4FVp+5X2r2hF3Aqlu9I1+qvRWtp3qfxsflel7Sntpyj6w2h5Ta96Leilx+94/Kq0GxoQqinRUzT2m7GSpNVHVuuOmXfo4LCDRX5WPwAUxHKn2Xq9XnXv3l3vvPNOqV6whmHoy95f6vqXr5fX41X6k+lXtf3p90xX27ltVda/rFpFtLpq7darWE8Nwxuq6qtVNarZKFUuU/mqtT301qFq9GYjzdk+R2O+HqOTj50slpMAAaBUBHrfvn310EMPKTIyMs/jHo9Hy5YtU2ZmpgzD0JEjR5ScnKyyZcsW+81ZrkWYZ2dnK9AVqFGRoxQRGiF3mltnPFfvxhhRYVHKyMjQuFvH6UwR3pDjfG0Xe88ZbWao8/zOeuIvTxRp24UxOmq0XvzuRW3os0HpKelK1+XtSHm9XuXk5CgjI8NSfdJut+vs2bNscX13q6LExETL3QVQkpKSkpSSkqL58+fLZrPJZrOpQ4cOPnk/BUsFev/+/dW4cWP94x//KHCD0qlTp9y/N2zYoPLly6tChQqWXP1SUlIUEhaiSZ0nXbPPsG3YtmKr7VJfqWwZtuXa9MGK/dX/tv5/+vWmaSorK0uBgdYb2YeFhZGLPstUxYoVLflVps1mU7ly5dSrVy/fr8UqC2XChAnKzs5Wp06dtG3bNm3evPmi03s8HsuNzGGBzaZpWrZfsr6x/EryOuf1en2+FsuM0FNSUuR0OjVx4kRJUlpampYuXcpaCAAoFSwT6K+++ipLEwBQanHpVwAACHQAAECgAwAAAh0AABDoAAAQ6AAAgEAHAAAEOgAAINABACDQAQAAgQ4AAAh0AABAoAMAQKBbjWEY1EZt1AbAEhyltXCbzaZDhw4pJydHpmlaLhSOHj0qt9ttycCzam2S5PV6VaVKFbZMAAj0wvL399fUqVMVEBBgyfrS0tIUEhJCbT62s/LLL78oJiZGw4YNY+sEgEAvjMzMTL322muqWbOmJes7fPiwatWqRW0+ZtWqVdq2bRtbJgCXrVR/h+71eq9KOzk5Oerfv/9VPbR/tWqz8nK7Fqz29Q8AAt1SkpKS9MMPP7CxBgAQ6L6MM5cBAAQ6AAAg0AEAINABAACBDgAACHQAAECgAwBAoAMAAAIdAAAQ6AAAgEAHAIBAL03F265O+Xa7nZ6GEtcvqQuwllJ7+1SHw6H169fr0KFDxX7TlOTkZDkcDnobLsk0Te3Zs0erV6+Wx+OxVJhv3rypNG9yAAK9ODcwe/fu1alTp4o90FNTU2Wz2bjbGi7J4/EoPj5e27Zts9RtYm02m+Li4iTVZSEDBHrRcrlcGjp0qGrUqFHsbSUmJmru3LncdQ2XXiEdDrVp00aPPvqo5Wpbs2a1pk/fwkIGimvHubSPhqzUDqzBSiNz1gOAQAcAAAQ6AAAEOgAAINABAACBDgAACHQAAAh0AABAoAMAAAIdAAAQ6AAAEOgoRk6nUxkZGcwIAACBXhiTJk1So0aNFBMTo+Tk5BLzucLCwnTo0CHuBw0AKDaWudval19+qfj4eK1bt04///yz/vnPf2revHkX3pOx2ZSYmKiQkBDL3db0fG1ly5alNh9it9tL1I4oAAL9mvj444/11ltvqWzZsrrttts0atQopaWlKSQkpMDpAwMD1adPH0suVMMw5HK55OfnR20+GOqvvPIKWyYApTfQT548mbuRt9lsqlixovbs2aOmTZsWOH1GRoaWL1+u2rVrW3IUe+DAAd1www3U5mM8Ho9cLhdbJgClN9DdbneeUZzNZrvohtFut8vhcMgwDBmGYcmRHrX55g7L7/uytRiSsiRlSrLivdGzJKX8t87f11xJkvmHx0794TFTUpikgD8OPSSlFrDZrvDHXUFJZwr4TNf94fNIUpKknD+0HSyp7B8+k/u/72lKciktLU1BQUEFHvXMzMzMd7/7gIAAORyOfNvprKysfOtzYGBgvvdMTU3Ns46bpqkyZcrkm87lcuXb1l/oPdPS0i643lnhHCfLBHpYWJi8Xq8kyev1Kj09XZUrV87TkcaNG6eUlBRJ0qZNm9S5c2d5PB7LjfQMw1BaWpri4uKsFwkWru38Rsvj8eTbEFphublcLgUGPiOH4+kCQsb3eb1emaaRL7ztdqOAIzHmH6aTbLb8O6mmacrrNfPNS5vtj+9p/vc9/7jzayvgc5r5tnkFvef/2jYk+al+/fpyuVzatWuX/P39/xciDod69OihHTt25JkX48ePV+/evfO857x58/TMM8/kCc/GjRtrwYIFeXZks7Oz1aBBgzxfreXk5Gj16tWqVq1a7uf39/fX+PHjNWPGjNx1xuPxKCoqSgsWLMgT9GlpaWrWrFmeHQ+n06l58+bp2LFjeuSRR2QYhhwOh5555hmVK1fO99Yz0yJp9vLLL+uGG27QvffeK4/Ho9tvv10bN2684PSxsbEaPny4IiIiLBkMcXFxqlOnDrX5YKBnZmYqKCjIkrWdPXvWkrVJ0pEjR3TdddflDizOK2hE+8eRommaCggIyHduSE5OjrKzs/ONJgMDA/OEv9frLfCnscHBwfl2EjIyMvJ9Rj8/vzwh/fuBkWEYSkhI0E033SS3251vuvOf849R4nA48o16vV5vviNQhmEUeE5MdnZ2vhF6QW17PB653e7caU3TlM1mK/A9Czpqm5qaqr59+2rFihWM0EuKgQMHqlu3btq4caN++uknTZw4sRB706aAkhZ6Vq7N6XTK6XRasr7AwECFhYUVatry5csXajqn06ng4OBCTRsQEFDo9yys8++Znp4uu90uu91e4HSFPUnVZrMVuv2CwrsgF/tcha39t6MRXp8/7G6ZQA8PD9e3336rrKwsORwOy44CAACwdKCf3/uy6t6/L5oxY0bu4bUmTZqoWbNmzBQAINDha06cOKFJkyYpLS1N8+fPZ4YAQDHiWqQoNk899VTuCL1z587MEAAg0OGL4uPjlZWVpfLly8vpdOrbb7/Vd999x4wBgGLAIXcUm4ULF0r67aczoaGhys7OlsvlUlBQkM6cOVPghR8AAIzQUcIsWrRI0m8nK54+fTr397QZGRlauXIlMwgACHT4gvOH17t3757vN7InTpxgBgEAgY6S7tSpU7kXSenQoYMkaf/+/bnP16xZk5kEAAQ6Srrfj8AjIyMl/e87dUlq2LAhMwkACHSUdMeOHcv99+7du5WcnKxZs2ZJ+u0yvTVq1GAmAQCBjpJu1apVkqRu3bpp0KBBCg8P1/79+/X000/r/fffZwYBQBHjZ2soFlOmTNGUKVPyjNgZlQMAI3T4OMIcAAh0AABAoAMAQKADAAACHQAAEOgAAIBABwCAQAcAAAQ6AAAg0AEAAIEOAACBDgAACHQAAECgAwAAAh0AAAIdAAAQ6AAAgEAHAAAEOgAABDoAACDQAQAAgQ4AAAh0AAAIdAAAQKADAAACHQAAEOgAAIBABwCAQAcAAAQ6AAAg0AEAAIEOAACBDgAACHQAAECg57dz50797W9/0y233KIePXroxIkTLFkAQKnisEIRJ06c0JgxY1SpUiUdOHBAbdq00dq1a1W+fHmWMACAQPcV7dq1y/33bbfdJqfTqaSkJAIdAFBqWO479H379iknJ0eVK1dm6QIAGKGXJCdPntQ777wjp9OZ53GXy6WRI0eqTJkykiSv16uhQ4fqrbfeUnBw8MULdzhks1n3nEDDMKjNF/ewbTbL1mfl2ljnfLtfnv+PQL8KgoKC1LT3y92CAAAI7UlEQVRpUzkceT+u2+3ODXm326127dpp0KBBatGiRb73cLvdGjdunFJSUiRJmzZtUufOneXxeGSapuVWvrS0NMXFxVlyw2LV2iTJNE15PJ58fZ3aSr6MjAy5XC7LbU+sXltKSoqOHTumRx55RIZhyOFw6JlnnlG5cuV8b/toWmQJ9e/fXx06dFDfvn0LNX1sbKyGDx+uiIgIS25c4uLiVKdOHWrzwdDLzMxUUFAQtfmY48ePq3r16tTmY5KSktSnTx+tWLGCEXpJMH78eG3dulX16tXThAkTZBiGBg0apGrVql3wNV6v15J7m/D90KM24Or2S9M05fV6ff6wuyUCffTo0Ro5cmSexy71HToAAFZiiUAPDAxUYGAgSxMAUGpx6VcAAAh0AABAoAMAAAIdAAAQ6AAAEOgAAIBABwAABDoAACDQAQAg0AEAAIEOAAAIdAAAQKADAECgAwAAAh0AABDoAACAQAcAgEAHAAAEOgAAINABAACBDgAAgQ4AAAh0AABAoAMAAAIdAAACHQAAEOgAAIBABwAABDoAAAQ6AAAg0AEAAIEOAAAIdAAACPTSwzAMGYZh6fqojX5JbfRLais9tTlKa6BnZmbq2WefVdmyZWWapuU6aHZ2tvz9/anNB5mmadkNqJVry8nJkdPptGS/tGpthmEoKytLpmnK6/XKZvPtMW6pDfTIyEi1bdvWkhuXgIAAjR07Vk888YTcbje1+diO5vHjx3XjjTdSm4+ZPn26YmNjlZmZSW2+FIIOhyZPniyHw/fjsNQGevny5VW7dm3L1hccHKzq1atTmw+GnmmaqlWrFrX5mNDQUF1//fXU5oMqVKhgiTpK7XfoXq/X0vV5PB5q89F+adW+aeXaWOeoj0C/hm666SZL19e2bVtq80F+fn4KDw+nNh90++23Uxv1XVOGadUziwAAKEX4HToAAAQ6AAC+6/cHqX39gHWpC3S3261hw4apWbNmqlOnTr6fYUyfPl3NmzdX27ZttXHjRp+vd9KkSWrbtq1atGihZ599Vjk5OZZaEWfNmqW7775brVq10ssvv2y5/rpgwQJVrVpVe/futUxNp0+f1qBBg9S6dWu1b99es2bN8ul60tLSNHDgQDVv3lx9+vTR2bNnLbOsnn/+eUVHR+vuu+/WkCFDlJ6ebrl1zDAMdevWTU2bNvX9nzGbpYzb7TaPHz9uHjx40GzUqJGZkZGR+9zKlSvN1q1bm26324yPjzcbNGhgJicn+2yt8+bNM3v27GlmZ2ebpmmaXbp0MX/66SfLLMvly5ebI0eONE3TNF0ul5mVlWW5/jpo0CDzzjvvNPft22eZmhITE83k5GTT5XKZ6enpZps2bcyVK1f6bD2PP/64OXPmTNPr9ZqrVq0ye/fubZlltXv3bjMrK8t0u93m3LlzzdjYWMutY++//74ZExNjtmvXzudrKXUjdLvdrmrVqikgICDf4ZVVq1bp0Ucfld1uV5UqVVS3bl2tXbvWZ2t1OBxKTU2Vx+NRZmam0tPTFRoaapll+dZbb2nIkCGaNGmSZs+ebbmLBI0dO1a9evWy1FEV6bdrQISFhcnPz09BQUGKiIjQoUOHfLaejRs3atCgQTIMQ5GRkdq6datlRrL169eXv7+/bDaboqKitG/fPkv1RY/Ho08//VQTJ060xFXw+A79dxISElSvXr3cvxs2bKgTJ074bD09evRQRESE/vKXv6hBgwYaNmyY6tata5nl9fPPP+vee+9VVFSUDhw4oMGDB1umtsOHD2vLli1q3769pde5jIwMbd++Xffee68l6gkJCZHT6VRCQoLlDkt/8MEHGjhwoKXqmjp1qsaMGWOZq05a7kpxqampuv/+++Xn5/fHrxY0dOhQtW7d+mJfP+S5lq/NZivRe21xcXF69NFHFRAQkOfxzMxMzZ07VwkJCbLb7Vq/fr0Mw1CPHj1Ur149n/kN/tSpU7VmzZp8jwcFBWnOnDlKTU3VRx99pFtvvVWtW7dWs2bNlJycrHLlyvl8bS+++KLef//9PH3Riutgv379NGrUKFWqVMkywWcYhuWuef7GG29o165dmjhxomVq2rdvnz7//HONHDnSMueoWC7QAwMD9fTTT+fbAJqmqZo1a170teHh4Tp06FDutaYPHDigW265pcTWWrlyZY0bN052uz3P4x6PR2XLltXo0aPVsmXL3IDr1auXtm/f7jOB3qVLF7Vo0SLf4+eXbfXq1XMvI2qapkJCQpSYmOgTgX6h2ux2u/bu3asff/xR69atk7+/v7Kzs7Vs2TLFxMT4RPAVdh0cPny46tevr169ellm+5Oeni6Xy2WpC+isXr1aixYt0tdff22prHj77bfVtWtXLVq0SL/88ovS09O1du1aNW/enEAvMQU5HIqKiirUdJLyjCLuuOMOvfnmm4qOjlZ6erq2bdumyZMnl9hag4KC9Ne//vWCz1eoUEE//vijBgwYIElasmSJxowZ4zPLMiIiQhERERd8vkePHpozZ45GjRolr9erU6dO6YYbbvD52pKSkvSvf/1Ldrs9d2fNbrf7zDkChVkHhw0bpjJlyuj555/3+W1OgwYNtGTJEnXq1EmHDx9WzZo1FRYWZpkwHzlypFavXm21qFC/fv10/Pjx3PXr9//3VaXySnExMTFKSUnR7t271bhxY0VGRuq5556TJMXGxurIkSNyuVzq16+f+vTp47N1pqSkaOjQoXK5XDp37pzq1KmjV155Jd8hel917ty53O/00tLS9Pjjj+vuu++2XH+NjIzUp59+apnzH7Zt26bo6Gi1atVKDodDpmmqefPmGjJkiE/WEx8fr7///e8KCwtTSkqK3nvvPVWtWtUSy6pJkyYKDw9X5cqV5fV65XK5NH/+fMutY3v27NFDDz3k0ydBl9pABwD8OVa+p72v4yx3AEDhR4GEOYEOAAAIdAAAQKADAECgAwAAAh0AABDoAACAQAcAgEAHAAAEOgAAINABAACBDgAAgQ4AAAh0AABAoAMAAAIdAAACHQAAEOgAAIBABwAABDoAAAQ6AADwUf8fRTkeWTQohQsAAAAASUVORK5CYII="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![figure_1.png](attachment:figure_1.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "In this example, the following Python code will attempt to find the optimal filter coefficients to match the output of a mass-spring-damper system given a horizontal force $F$ as input. The force varies over time sinusoidally. The output of the system is the position of the mass. The output is calculated using the Runge-Kutta RK4 algorithm to integrate the first order ODE \n",
    "\\begin{equation}\n",
    "m\\ddot{x} + b\\dot{x} + kx = F.\n",
    "\\end{equation}\n",
    "\n",
    "and propagating the dynamics at each time step $T_s$.\n",
    "\n",
    "Gaussian noise has been added to the output as well as some uncertainty in the system's physical parameters (mass $m$, spring-constant $k$, and damping coefficient $b$). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class msdParam:\n",
    "    \n",
    "    def __init__(self):\n",
    "        # Physical parameters\n",
    "        self.m = 5.;  # kg\n",
    "        self.k = 3.;  # N/m\n",
    "        self.b = 0.5; # N m s\n",
    "\n",
    "        # Simulation Parameters\n",
    "        self.t_start = 0.0;  # Start time of simulation\n",
    "        self.Ts = 0.1;       # sample time for simulation\n",
    "\n",
    "        # Initial Conditions\n",
    "        self.z0 = 0.0;\n",
    "        self.zd0 = 0.0;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "P_ = msdParam()\n",
    "\n",
    "\n",
    "class msdDynamics:\n",
    "    '''\n",
    "        Model the physical system\n",
    "    '''\n",
    "\n",
    "    def __init__(self):\n",
    "        # Initial state conditions\n",
    "        self.state = np.matrix([[P_.z0],    # initial position\n",
    "                                [P_.zd0]])  # initial velocity\n",
    "\n",
    "        alpha = 0.02  # Uncertainty parameter\n",
    "        self.m = P_.m * (1+2*alpha*np.random.rand()-alpha)  # Mass of the msd, kg\n",
    "        self.k = P_.k * (1+2*alpha*np.random.rand()-alpha)  # Spring constant, N/m\n",
    "        self.b = P_.b * (1+2*alpha*np.random.rand()-alpha)  # Damping coefficient, Ns\n",
    "        self.Ts = P_.Ts  # sample rate at which the dynamics are propagated\n",
    "\n",
    "    def propagateDynamics(self, u):\n",
    "        '''\n",
    "            Integrate the differential equations defining dynamics\n",
    "            P.Ts is the time step between function calls.\n",
    "            u contains the system input(s).\n",
    "        '''\n",
    "        # Integrate ODE using Runge-Kutta RK4 algorithm\n",
    "        k1 = self.derivatives(self.state, u)\n",
    "        k2 = self.derivatives(self.state + self.Ts/2*k1, u)\n",
    "        k3 = self.derivatives(self.state + self.Ts/2*k2, u)\n",
    "        k4 = self.derivatives(self.state + self.Ts*k3, u)\n",
    "        self.state += self.Ts/6 * (k1 + 2*k2 + 2*k3 + k4)\n",
    "\n",
    "    def derivatives(self, state, u):\n",
    "        '''\n",
    "            Return xdot = f(x,u), the derivatives of the continuous states, as a matrix\n",
    "        '''\n",
    "        # re-label states and inputs for readability\n",
    "        z = state.item(0)\n",
    "        zd = state.item(1)\n",
    "        F = u[0]\n",
    "\n",
    "        # Equations of motion\n",
    "        zdd = (F - self.b*zd - self.k*z)/self.m\n",
    "        # build xdot and return\n",
    "        xdot = np.matrix([[zd], [zdd]])\n",
    "        return xdot\n",
    "\n",
    "    def outputs(self):\n",
    "        '''\n",
    "            Returns the measured outputs as a list\n",
    "            [z] with added Gaussian noise\n",
    "        '''\n",
    "        # re-label states for readability\n",
    "        z = self.state.item(0)\n",
    "        # add Gaussian noise to outputs\n",
    "        z_m = z + random.gauss(0, 0.001)\n",
    "        # return measured outputs\n",
    "        return [z_m]\n",
    "\n",
    "    def states(self):\n",
    "        '''\n",
    "            Returns all current states as a list\n",
    "        '''\n",
    "        return self.state.T.tolist()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzIAAAI8CAYAAAAqSJ7pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAIABJREFUeJzs3Xl8VOX1+PHPncm+73tCAgQS9k02QVHctS7ghlqtijv219pK1bbWrWqt1q/FpWqtdauiIipuIIvsO4EkJIHshIQlLIGQPZP7++PJTIIECMnM3JnkvF8vXr1hJvc5yDTMmeec82i6rusIIYQQQgghhBsxGR2AEEIIIYQQQpwpSWSEEEIIIYQQbkcSGSGEEEIIIYTbkURGCCGEEEII4XYkkRFCCCGEEEK4HUlkhBBCCCGEEG5HEhkhhBBCCCGE25FERgghhBBCCOF2JJERQgghhBBCuB1JZIQQQgghhBBuRxIZIYQQQgghhNuRREYIIYQQQgjhdiSREUIIIYQQQrgdSWSEEEIIIYQQbsfD6ABcyaFDh3j77bdZtmwZe/bswcfHh/j4eMaPH8/s2bONDk8IIYQQQgjRStN1XTc6CFeQnZ3NnXfeydGjR+nfvz8DBgzg2LFjFBQUsG/fPrKzs40OUQghhBBCCNFKSstQOzEzZ86ksbGRN954gwULFvDSSy/x5ptvsmTJEj7++GO7rWWxWJgzZw4Wi8Vu9xQ9k7xWxJmQ14voLHmtiDMhrxfRWUa8VmRHBnjiiSeYO3cuf/nLX7jxxhsdulZTUxNDhgwhOzsbT09Ph64l3Ju8VsSZkNeL6Cx5rYgzIa8X0VlGvFZ6/Y5MQ0MDCxYswNfXl2nTphkdjhBCCCGEEKITen2zf1ZWFjU1NYwZMwYvLy+WL1/O2rVraWhoIDk5mUsvvZSoqCijwxRCCCGEEEK00+sTmcLCQgDCwsJ44IEHWLJkCZqmAaDrOi+//DJ//etfufzyy40MUwghhBBCCNFOry8tO3LkCABLly5l1apVPPHEE6xZs4alS5dy5513Ul9fzyOPPEJeXp7BkQohhBBCCCGsev2OjHXWgcVi4aGHHjqu2f/hhx+mvLychQsX8s477/D3v/+9W2u1tLRw8OBBAGpqavDy8urW/UTP1tTUBEBtba00WIrTkteL6Cx5rYgzIa8X0Rm6rlNdXQ2o97vO0usTGT8/P9v1Nddcc8Lj06ZN44cffmDjxo2dvqfFYunwL/HgwYOce+65AIwbN64L0YreaOzYsUaHINyIvF5EZ8lrRZwJeb2IzqqsrCQ6OrrDx0wmE2az2W5r9fpEJi4uDgAfHx9CQ0NPeDw+Ph7AtpPSGa+//jqvvvqqfQIUQgghhBDCTUydOvWkj82aNYsHH3zQbmv1+kRm0KBBgBrD3NTUdMK2qbWHpv3Ozencf//93HvvvSf8fk1NjW0nZs2aNfj6+nY1bCGEEEIIIVxCXV0dEydOBGD9+vX4+/t3+DyTyb7t+b0+kYmNjSUtLY0dO3awceNG21+C1YYNGwAYPHhwp+9pNps73DZr3xPj6+t7RsmREEIIIYQQrs7Ly0sOxHSmmTNnous6f/vb36isrLT9fm5uLu+++y6aph03BEAIIUTn6DrMmwezZ8Njj0F+vtERCSGE6Cl6/Y4MwBVXXMHq1av58ssvufzyyxk5ciT19fVkZGTQ1NTE9ddfz0UXXWR0mEII4VaOHIFbb4Wvv277vRdfhKeegkceMS4uIYQQPYMkMq2ee+45Ro0axdy5c9mwYQOapjFkyBBuuOEGrrrqKqPDE0IIt9LSArfcAt98A15ecMcdUFwMCxfCo49CWBjcfbfRUQohhHBnmm49SEU4XG1tLSNHjgQgIyNDemSEED3WCy/AH/4AXgHH+MU//0B2zRKONR4jpfqXrPrbw5gbw1i3DsaMMTpSIYQQ3WXUe1zpkRFCCGFXe/bAE08AQWWEPzKWebteZ8fBHZRXl7OK5wn4zTgs3pXMmqV2boQQQoiukERGCCGEXT3zDNQ1NBFw+w3sac4lLjCOL67/gnnXzyMpOIlj3gWYbrmS9Zvree89o6MVQgjhriSREUIIYTe7dsFbbwHnPs2x0LUEeQex6vZVXJN+DdPSp7HolkWE+YbRErcOJvyD55+XXRkhhBBdI83+LkzXderq6owOQ5wBX19fNE0zOgwhDPPvf0OzbwXapBfQgbeueIuU0BTb4wMjBjLn0jnc/MXNcM5f2TnnVhYvTkAGQwohhDhTksi4sLq6OlvjlHAPMsRB9GbNzfCf/wCTnkM3NzApaRLXD77+hOfNGDKD1za+xpqyNTDlCebM+bckMkIIIc6YlJYJIYSwix9+gPKjFTDmLQCemvJUhzuUmqbxwgUvqC+Gfcg3P+2jvNyZkQohhOgJZEfGTaxZswZfX1+jwxAdqKurY+LEiUaH4ZIaG6GiAhITwWw2OhrhaO+9B4z8D5gbmZg4kfNSzjvpcycmTmRc/DjWl6+HMW8wd+4TPPSQ82IVQgjh/iSRcRO+vr5SsiTchq7Dm2/C00+rRMbPD+69F559Fry9jY5OOEJNDXz7XQvc9Q4A946+95TP1zSN347/LTfOuxHOep2PPnmMhx7yckaoQggheggpLRNC2N3TT8N996kkRtOgthb+8Q845xz1hlf0PD/8AHUxSyC0hBCfEK4ddO1pv2f6oOlE+kaDfyVbqhaTn++EQIUQQvQYksgIIezqgw/gL39R108/rRKXL7+EsDDYsAEeeMDY+IRjfP45MEIdCnPL0Fvw9Tx9KayHyYMbh7YOAxjyCZ995sAAhRBC9DiSyAgh7KayEv7f/1PXv33sMPrkp7l63sV8armZP7+zDJNJ9VHMnWtsnMK+GhthwfcNMGABADcNvanT33vjkBvVRdqXfLNQxs0LIYToPOmREULYzezZcPgwDJicyQchUznw0wHbY59on3DBn15i0VO/4dFHYdo08PQ0MFhhN6tXQ03kUvA5SmxALOMSxnX6e8cnjCfOP4kKdrH+4PccPjyN0FAHBiuEEKLHkB0ZIYRdFBS0Tq0KLeLApRdxoPYAaRFpzLl0DrcMu4UWvYVFpt8SNOFTiovh3XeNjljYy8KFQPoXAFyTdg0mrfP/tJg0EzcMnQ5AS/8F/PijIyIUQgjRE0kiI4Swi5dfBl3XCbv9Dg417mNY9DDW3LGGWWNn8f7V7/PwxIcBaLh4JoQW8uyzYLEYHLSwi4WLWiDtK0A18J+py1MvVxep3/Pd9y32DE0IIUQPJomMEKLbDhxo3WEZ9iGHgpbj6+HLVzd+RaivqhHSNI1npz7L5KTJNFCN5+V/oLQUFi0yNm7Rffv2wda9W8G/kgDPQCYnTT7je0zuMxlfcwAE7OPbLVvQdQcEKoQQoseRREaIn3nkkUdIS0vj0UcfNToUt/HBB1DX2IDHZbMBePzcx0kOST7uOR4mD964/A00NJr6z4Pobbz1lgHBCrtasgTouxiAKSnn4mk+88YnL7MXF/a9EIADod9RWGjPCIUQQvRUksgI8TOapqFpmtFhuJWPPgKGzKXZZy8JQQk8NKHjI9oHRw3m+sGt43bPfYoFC2DPHufFKexvxQogZQkAU1Omdvk+v0i7TF2kfsdPP3U/LiGEED2fJDJCdECX2pZOy8uDzZt1GP9/ADxw1gN4mU9+Qvvj5z6uLtK/xOJfps4fEW5rxZoG6LMS6F4ic1G/i9RF3CZ+XF5tj9CEEEL0cJLICCG65aOPgMQ1EJuBr4cvd42665TPHxQ5iCnJU0BrgZHv8uWXTglTOMChQ5B7dB141hHhG8WQqCFdvldScBKxPilgsrAkf7X0yQghhDgtSWSEEN0ybx4w/ANAHW4Y7hd+2u+ZOXKmuhj5Dj+tsHDokAMDFA6zZg2QvByAqX3P63ZJ5gWpUwA4GLCcgoJuBieEEKLHk0RG9Epff/01M2bMYNSoUYwZM4brr7+eTz/91Oiw3E5REeTubIRBnwFw89CbO/V909KnEeITAiG7aElaxnffOTJK4SirVgGJqwGYlDSp2/eb2vdcdZH8k0qShBBCiFOQREb0Oo8++iizZ89m69at1NfXYzab2b59O3/5y1/43e9+Z3R4buXbb4F+i8DvEDEBMapkrBN8PX25Nv1a9UXafL76ymEhCgdaucoCiWsB+yQy5ya3JjJxm1ix7li37yeEEKJnk0RG9Crvv/8+8+fPR9M0brnlFtasWcP69etZv349s2bN4rvvvmPJkiUytayTvv0WGPIJADcMvgGzydzp77067Wp1kfYVi5focjimm2luhs1l2eBdjb9HIEOjhnb7nskhyUR69gFzM8vy19khSiGEED2ZJDKi12hsbOS1115D0zSuvvpq/vjHPxISEgJAQEAADzzwAHfffTdHjx41OFL3UFMDS39qhtRvAbhu0HVn9P1T+07F39Mfgsqp8t1MRoYjohSOsn07NESrsrIJSePPKIk9lbOTJgJQ0ryOahleJoQQ4hQkkekldF298XSVX0ZMJFq1ahVHjhwB4P777+/wOXfddRfe3t7ODMttrVwJTdFrwbeKMN8wxieMP6Pv9/Hw4dLUS9UXafP58UcHBCkcZsMG2vpjEs+2233PS1WvIz1+nVpDCCGEOAlJZHoBXYdJkyAgwHV+TZ7s/GQmOzsbgNjYWBITEzt8TkBAAIMHD3ZmWG5r6VJggNqNuaT/JV36RP7KAVeqi/4/sHixHYMTDrdxI7b+mImJE+12X1tCnLCO1atlBrMQQoiTk0Sml5CWDzh48CAAUVFRp3xeTEyMM8Jxe0uXAqlq3NjlqZd36R4X9L1AXcRmsHLTIWpr7RSccLg1Ww9CaDEAZ8WfZbf7jogZgQfe4HeQnzIL7XZfIYQQPY+H0QEIx9M0VQbkSm8S/fwkuXJnhw/DloLd8IssTJqJi/td3KX7xAbGMihyEDmVOTTFL2Pjxumce66dgxV2V1sLuVWbAEgJSlWjtO3Ey+xFWvAoso+sZUvlOqC/3e4thBCiZ5FEppfQNPD3NzoKY4WHq4Ma9+/ff8rn7du3zxnhuLUVK0BPXgrAmLgxnToE82SmpkwlpzIH+i5h1SpJZNzBtm3QErsRgAlJ9tuNsZrSfzzZm9dyJGAde/bcQmys3ZcQQgjRA0hpmeg1hgwZAsCePXsoKyvr8DnHjh1j+/btzgzLLf30E5CiEpnzk8/v1r2mpkxVFylL1AGLwuVlZABxKpGxZ1mZ1cQ+rfeM28TmzXa/vRBCiB5CEhnRa5x99tkEBwcD8Prrr3f4nLfffpv6+npnhuWWVq3W2xKZlO4lMucmn4sJE0TsZFVmuZwn4wYyMoD41kQmzv6JzOi40eoiehsbNjXb/f5CCCF6BklkRK/h7e3Nfffdh67rfPnllzz77LNUVVUBaifmtdde46233rIlO6JjtbWQUVIEwWV4mjw5O6l7o3dDfEIY0nqY4rGQtciGmOtbn1MBgXswYWZk7Ei7379/WH+8CQTPepbn5Nj9/kIIIXoGSWREr3Lbbbdx9dXqRPn333+fiRMnMm7cOMaNG8err77KZZddxvnnn49uxEE3bmLjRrAkqt2Y8Qnj8fP06/Y9JyZOUBcJa6W8zMU1NUFu1RYA+oek2+Xv/+dMmolBoaMA2FYptWW9QXU1shsrhDhjksiIXkXTNJ5//nn+9re/MWLECHx8fLBYLAwePJinnnqKl156yfY8TcaqdWjNGqDPSgCmJE+xyz1t55AkrmX9ervcUjhIbi40h28D4KzEEQ5bZ1JfVV52xG8zp5nPIdxUSwu8/DKkpEBQkDpj7MILISvL6MiEEO5CppaJXunKK6/kyiuv7PCx5557jueee87JEbmP1auBfmsA+x2EOMG6IxO7mfXzGwBvu9xX2F9GBhCzFYARMcMdts74PqOZsxmI3UxWFkyd6rClhAHq6+GGG+Drr4//vcWLYdQoePNNuOMO4+ITQrgH2ZERQnSarsPqrfshTB1UaDuFvZv6hfYj3CcCPBrZcSSD6mq73FY4wNatQIzakRke7bhEZnRsa8N/zFa2ZUnDf09iscAvf6mSGK+AGm578SNu/exu7p37ZybftIbmZpg5Ez75xOhIhRCuThIZIUSnFRdDVcBaAAZFDLbbQYiapnF2n9bdnYS1MnLXhW3LPQZhBQAMd+COTGp4Kp66H3jWsyavwGHrCOd79ln4/HMw911J2B+H8N6xW3h/+9v8K/cZVg44m76PXIfuW8ntt8POnUZHK4RwZZLICCE6beNGIFGVlZ2dZJ+yMqtx8ePURdwmNmyw662FHWXtywJNJ8I7lij/KIetY9JMJPupaXZb92Q6bB3hXNnZ8PTTQL+FcOv57G0oITEokd9P+D0zhszArJkp8vmcwPsvol4/wl13qV4aIYToiCQyQohOa5/ITEiYYNd7j4kboy7iNql1hMupqoIDHqqsbGSs4xr9rUbEDgNgV0OmvJntAXQd7rkHmsK3Yr5pOhaamZY+jez7s/n7RX/nf9P/x6a7NxHlH0V1wFZMM65lxcoWPvjA6MiFEK5KEhkhRKdt2NQMcaruy179MVa2noiInWzYdsSu9xb2kZMDRKtEZnS848rKrM7urxKZprBMSkocvpxwsMWLYc1aC9pVM7GYa7ig7wV8PP1jgryDbM8ZETOCH27+AX9Pf1qSF8PYV3n8cWhoMDBwIYTLkkRGCNEpFgtsKtkBnnX4eQQwMGKgXe8f7hdOUlAyALuattB6VqlwITk5QJQ6sXRI1BCHrzcyTiUyRG+Tkbw9wNNPA6PfRo/dTLB3MB9e8yFeZq8TnjcydiR/v/DvAGgXPsKu6iLefNPJwQoh3IIkMkKITtmxA+qCMwBVVmTS7P/jY2xCW3nZtm12v73opuztOkRlAzA4arDD1xsW3ZrIhOxiQ6Zktu5s9WpYuf4YnP8nAJ45/xmiA6JP+vx7xtzDecnnoXvUwblP8uKL0CzD64QQPyOJjBCiU7ZsAWLVie6jYkc6ZI0xsa2JTPxGSWRc0Nb8veB7GA0TaRFpDl8vxCeEEC0JgHVFsiXjzl5/HRjzBvgdJDUslXvH3HvK55s0Ey9c+IL6YtiHlNXks2CB4+MUQrgXSWSEEJ2SmUm7RGaUQ9YYHdfaJxO7WZ1XIlxK9n5VVpbg1w8fDx+nrJkapCaX5RySyWXuqrISPvuyFia+CMBjkx/Dw3T687jHxI3higFXgKkFzvkrr77q6EiFEO5GEhkhRKdsy2yBGFVa5qhExnbAYlgRW7JrHLKG6JojR+CgSSUyw2IdX1ZmNTppEAD7mvOk4dtNvfsuNKV9BAH7SQ5J5uahN3f6e/98zp/VxdD/sXTDXnbscFCQQgi3JImMEKJTMkqKwecoXiZv0iPSHbJGpH8k4T7qbJLt+3NoanLIMqIL2jf6j4hzXiJzVop6rekRufIm1g3pOvz3v8AY1a3/wFkP4Gn27PT3j40fqyYkmptg9FvqXkII0UoSGSHEaVVWQqVZlZUNiRx6Rm9EztTwGDUNqzk0m7w8hy0jzlBODhCpEpnBkc5LZAZFtibNEbkyucwNZWZCbtVmiNuMl9mLX4341Rnf49djf60uxrzBex82YrHYN0YhhPuSREYIcVpZWUCsKisbE++YsjIr21jfqO3yxtWFqIllrYmMEyaWWdl2/4Iq2JQl5wu5m//9Dxj9NgDT06cT4RdxxveYPmg6Mf4xELiXPX4LWbzYzkEKIdyWJDJCiNNq3+g/0kETy6xsb5KjssnOduhS4gxszd8PPkfQ0BgQPsBp6wb7BBNsigVgQ5Fs0bmTlhb4eG4TDPocgDtG3tGl+3iZvZgxdIb6Yuj/VHIkhBBIItOhqqoqJkyYQFpaGhdffLHR4QhhuG2ZusMnllm17chIIuNKtu/ZCUCMTx+nTSyz6hekdmV2Hs516rqie7ZsgTLzcvA7SIRvBFOSp3T5XjcNvUldpH3F1z9US/+cEAKQRKZDzz//PEeOHEHTNKNDEcIlbN5ZAf6VmDAzNGqoQ9ey9V8ElZO5Uw5BdAVHj0KlrjrtB0UPdPr6w+NVInNAy6W+3unLiy765htg0GcATEuf1qmRyyczOnY0qWGp4FlHVcyXrFhhpyCFEG5NEpmfWbt2LV9++SXXXXcduq4bHY4Qhmtuhrwjajemf3A6vp6+Dl0v2CeY+IBEAEprt1MjU5gNl5MDhKsdmUHRzisrsxqV1NbwX1Dg9OVFF33zXTOkfwHAdYOv69a9NE1r25UZ/Bnz53c3OtET6DoUFsL330NZmdHRCCNIItNOQ0MDjz/+OKmpqdx5551GhyOESygogKZwlciM6+PYsjKroTFtfTI5OU5ZUpxCbi62RMaZ/TFW7SeXyQhm97BnD2zetw78DxDqHdatsjKr6enT1UW/RXzxbTXyWWPvpetqkMSQIdC/P1x2GSQlwdlnw8aNRkcnnEkSmXbmzJlDeXk5Tz75JGaz2ehwhHAJqtHfehCmYxv9rYZESp+MK9m5E1siMzDc+aVltslloUVk50ltmTv4/nsg9TsALh1wSbfKyqyGRA2hf2gqeDSwx/879bNJ9DpNTXDbbXDzzWq32NMTBg4EsxnWrIHx4+G994yOUjiLJDKt8vLy+O9//8v06dMZNco5nzoL4Q4yM4EolU0Mjx7ulDXbj2CWRMZ4OwuaIUzVdBmxIxMTEIMPwWBqYXOx1Ja5g0WLgNRvAbis/2V2uaemaVw7qHVXZtA8lSyJXqWlRSUwH3ygEpcHntjO//30Xx54bw7/Xvk1V91QRUsL/OpX8P77RkcrnEESGUDXdf70pz8RFBTE73//e6PDEcKlbMmugdAioF2C4WDtJ5dJaZnxcipKwdyEl+ZDYnCi09fXNI1EX7Urk3tAJpe5upYWWLRuN8RkoqFxcX/7Tf+clj5NXaR+x7c/NNrtvsI9/O1v8Nln4Bl8gEkv38xrDOGBH2/n1z/8mtsXXcXq0alc/LuPAbj7bjU5T/RsksgA77//Ptu3b2f27NkEBwcbHY4QLiVjdw5oOqGeUUT6RzplzfTIdDQ08K8kp3S/U9YUHdN1KD6qysqSg1Ixacb8szE4SiUyu+tzpTfCxW3bBofD1XbJ2PhxXToE82RGx40mzDsSvGpYU7aGI3JGaq+xbh388Y9AaBGhv5/I8kPqQKEpyVOYlj6NvqF9OVB7gIWBNzHgrqdpaIDrroPaWmPjFo7V6xOZiooKXnnlFcaOHcvVV19tdDjCQY4ePcrw4cNJS0vjhx9+OOVz/+///o+0tDQuvPBCJ0Xnuo4ehb0WVds1NNo5uzEAfp5+9AnqC8Cu+mwZuWugPXugIUB12A+OcX5ZmdWY5DQA6gNy2S+5rUtbvBjo9yMAl/a/xK73NmkmLh1wEQAtKYtYssSutxcuqrkZ7r0XdN8D+N8/lf2WfPoE92HTXZtYdtsy5l0/j9wHcnls0mMA7Ix/nOBLX6KoCJ5/3uDghUP1+kTmqaeeorm5mSeeeMJu97RYLDQ1NXX4SxgjKCiIyy5TddqffvrpSZ/X0tLC/Pnz0TSN667r3rjQniA3F1t/zIg45yUyAMNjW9eL3C4jdw2Un4+t0T8t0rhEZlhsa8N/pEwuc3U/Lm6BlKUAXND3Arvf/+J+raVq/RaydKndby9c0BtvwLasZjxmXE+NZwn9Qvux9s61jI4bbXuOl9mLv079K3+74G8AHBv3CMRt5IUX1Ihm4Twnew/c1NSExWKx61rdHyPi5n766SeCg4N5/PHHj/v9xkZVe7tv3z5++ctfAuqT+vDw8NPe8/XXX+fVV1+1f7CiW2bMmMH8+fNZt24du3fvJiEh4YTnLF++nH379uHh4cH06dMNiNK1tE9knNUfYzUkaghf7fgKorLZsUON2RTO1z6RMWJimVW6dQRz+A5y8yycc45MlnRFjY2wYsc2mHgQP48AxsaPtfsaF/VTOzLEbWHRl/uBKLuvIVxHTQ088www9lWaE5cR4BXAVzd+RWxgbIfPf3jiw2yq2MRnOZ/h98ubqX0pk2ee8eHdd50bd282duzJ/38/a9YsHnzwQbut1esTGU3TOHr0KJs2berw8YaGBjZt2oSmaTQ0NHTqnvfffz/33nvvCb9fW1t7yr9cR9J1ndom1ykU9fP0Q9M0p645bNgw0tPTycvL47PPPuO3v/3tCc+ZO3cuAFOnTu1U0trT5eVhWCJjG7kbsUM+gTdQQQGGniFjlRKSgln3xuJZz8b8Uu6hr2GxiJPbsgUa4lW915SUc/E0e9p9jeiAaIZGjCDrwFbyLYupqLiJuDi7LyNcxOuvw/763WhT/4wOvHTRSwyOGnzS52uaxptXvMnqstVUkA8T/sGHHz7Gk0+qs2aE423YsAE/P78OHzOZ7FsM1usTmdzcjifglJeXM3XqVJKSkli4cOEZ3dNsNnd4Do2np/1/oHeGrutMencSa8rWGLJ+R85OPJuVt690ejIzY8YMHn/8cb744gt+/etfH/f3tG/fPlauVDHdcMMNTo3LVWXmH4ZhFQCn/IfDEfqH9VcXYQWSyBgot6AGhqojs41MZMwmM9EeqVRYssmqyAdJZFzSihVAikpkLkiZ6rB1LkqdStaBrZC8jGXLbuLmmx22lDBQXR288AIw9VF0z2NMSJjAzFEzT/t9ob6hvHDBC9wy/xZMU/5K87ZbefHFBP75T4eHLFDvd531nrfX98j0FhrOTRhc1RVXXIG/vz8HDhxg2bJlxz32+eefY7FYSEhIYOLEiQZF6Fqy9+YBEOmVSJB3kFPXTg1PVRdB5eTku85uYm+Ts081KAV5hBPuZ+wuZUqwSm6Lj0jTlKv6aUUzJK0C4PyU8x22znnJ56mLlGXS8N+D/e9/cEDLgWEfATDn0jmdnpx409CbmJg4kRZzLUx6nnfegaoqR0YrjNDrd2R6A03TWHn7yl5fWgbg5+fHlVdeyccff8zcuXO54ALViKrrOvPmzUPTNK6//nqnx+WKGhthd71KZNIinN8bEeYbRpBnKEebDrNjfyG6PhQDXjK9WksL7Dqmysr6hRi3G2M1JL4/qw/BAUsBjY3g5WV0RKI9iwVW7twGZx0jwCOYodFDHbbW5D6TMWGmJaxxqmrLAAAgAElEQVSQpT+VAc4/30g4lq6jdlCmPAGaztVpVx/X3H86mqbxzHnPcP7756ONfofa5X/mnXei+d3vHBayMIDsyJyCEW+0HUXTNPy9/F3ml5H/bWfMmAHAmjVrqKhQZVOrVq2ioqICs9nMtGnTDIvNlRQWgh6mEplh8WmGxDAgQu3KVHvlU1lpSAi9WkUFNAapur4hccYnMiOS1I6MHlpAUZHBwYgTZGXBsTC1GzM5+WyHnjkU5B3EyGj1prZUW8a+fQ5bShhkxQrILCuEQZ8D8OSUJ8/4HlOSpzA2fiy6uR7Gv8Krr6qEW/QcksicRHx8PLm5uWfcHyNc34ABAxg5ciQtLS18/rn6AWkdyXzBBRdIk3+r3FwgQiUy6RFGJTLSJ2Ok9hPL0g3Ylfu5VOmbcmmrVwN9VgJwTp/JDl/vgv6t5WXJy1jjOi2gwk7+8x/grDdA07mk/yUMix52xvfQNI1HJz2qrse8RUlZA/K2rmeRREb0SjNmzLCVk+3bt49ly5ZJWdnP5OVhS2TSDEpkUsNa+2TC8uWNqwHy84GwfKBdz5KBbAMgQovIzZOPVV3NmrW6rT9mUtIkh693bp9z1UXSapVEiR7j2DH4/OsaGPkOALPOmtXle/1iwC9ICEpA9z0I6fN5/317RSlcgSQyole65JJLCAkJYf/+/fz+97+nublZmvx/ZnteI4SpU8SMSmRkcpmxCgqAUFXD1S+0n7HBAAlBCZh1L/BoJKOozOhwxM+szC6EgH14al6MiRvj8PXGJ4xXF+H5/LRBak97knnzoDZlLvhW0Te0L5f0v6TL9zKbzNw58k71xei3+PJLafrvSSSREb2Sl5cX11xzDbqus3HjRtmN6cC2skIwWfAxBRAXaMwhDbZEJlx2ZIyQW3gMAvYD0DfU+HHHZpOZaC8VR16lTC5zJfv3QxmqvmtUzBh8PHwcvmaobyj9gtSHLFsPrqOuzuFLCif54ANghDrBcubImZhN3TsA946Rd6jprSnLaPAvoLWaXPQAksiIXmvGjBm2oQPS5H88XYfCKlVW1j84zbDhDLbSsqBycgtcZ+peb5G7pxiAQI8wgn2CDY5GsY5g3lUtiYwrWbsWiF8PwNl9xjtt3XP6ql10S+xaTnKutXAzBw7Asm0F0GcVJs3ErcNv7fY9k4KTuKjfReqLoR/xySfdvqVwEZLIiF4rKSmJ9PR0NE2TJv+fKS+H+oDWiWVxxpSVgRrBHOwVAkBxVSFNTYaF0uuo0cuqrCw52PjdGKvBcSqRqTIVUF9vcDDCZu1aIGEd0K7kywkmJk5QF4lrWLXKacsKB1qwAFqGvgfARf0uIj4o3i73vWXYLepi2If8tFyXSZg9hCQyotc6cOAAO3eqiUxSVna89o3+g6KMS2Q0TbONYLYEy8hdZ9q9G5r81X/wtCjXSWSGxrf1TRUXGxuLaLN2Ux1EZwIwLmGc09adkNCayMRtZOXqZqetKxxn3he67QDM24bfZrf7XpN2Df6e/hBegB67ga++stuthYEkkRG91scff0xzczNJSUnS5P8z7ROZgQaP3ZWGf2MUFWFr9O8blmJsMO20H8FcWGhsLELRddhcsQXMzUR4x5IY5LzDKdMj0wnwCAavWlblZ9LS4rSlhQNUV8OirM0QWoyv2Y8rB15pt3v7e/lzTfo16othHzFvnt1uLQwkiYzolbKysnj33XfRNI077rjD6HBcTk6uDhEqazBqYpmVrU9GGv6dqrgYCFVbHq7Q6G/VNoK5kJ358q7VFRQXQ02IKiub2GecU3vqTJqJiUmqlK06eK36EEa4re+/h6YBqhP/FwOvwM/Tz673v2HwDeoi/QsWL2nhyBG73l4YQBIZ0aucf/75TJo0ieuuu466ujoGDRrEtddea3RYLiezcB/4HEHD1PbG0SDtd2SktMx5VCLTuiPjQolMn5A+mHQP8Kwnq6TC6HAEsHkzEL8BgAmJzisrs2rfJyPnybi3eV/oMOgzAK4fbP+S7wv6XkCAVwAEldMctYlFi+y+hHAySWREr7Jnzx4OHjxIZGQk06ZN4+2338Zs7t5Yx54o74D6WDPON8UpY1RPxXYQY1i+lBI5UWFRC4S43o6Mh8mDCA9V6pazVyaXuYItW4DYDABGx452+voTE1tLgxPXsmGD05cXdlJfDws2bYHQEnzMflyaeqnd1/Dx8OGy1MvUF+lfsGCB3ZcQTuZhdABCOFNubq7RIbi8mho4qKlEZnC0sWVl0G5HJng3BaV1gK+h8fQWOyv2Qmo9JsxO7XnojOSg/uw/nE9JdT4wxehwer31W6thfD4AI2JGOH39cQnj0NDQQ4tZu24vEOP0GET3LVkCdYkqs7g09RK7l5VZTUubxqfbP4W0+Xz38fNYLCCfZ7ov2ZERQhynsBBbo//QWOMTmXDfcNsI5l3VhTTLYCKnKK5SZWUxvkl4mj0NjuZ4g2JUcnvAUoDFYnAwvZyuw+bdalpZlE88kf6RTo8hyDuIASGDAcg9tlYOxnRTX30FDFCJzC8GXOGwdS5NvRRPkydE7ORgS6EaHS7cliQyQojj5OdjS2SMbvQHNYI5NVy9cbUEF7Brl8EB9QJ1dXCwRSUy/cJdZ2KZ1fBE9XpoCSmgrMzgYHq5sjI46qfKysbEjzQsjsl9VZ9MS+x6MjMNC0N0ka7DtyvLIW4LGhqXD7jcYWsFeQcxwdpX1W8RCxc6bCnhBJLICCGOU1BA2+jlcGNHL1v1D5eGf2cqKcHW6D8g0nX6Y6wGRMgIZlexeTMQ09ofY2AiMzbuLHURu5lNmwwLQ3TRjh1Q4f8NAGfFjiPKP8qh613U9yJ10W8RS5c6dCnhYJLICCGOk1dQD8Fq22NA+ACDo1GSg5PVRXCpvHF1gvYTy/q5UKO/VftJdgUFurHB9HLtG/1HxhiXyIyKHaUuYrewcZO8JtzNwoVAf7U1cmWa48rKrC7q15rIpCxl/cZmqqsdvqRwEElkhBDHyS4vBk3HRwt0+KdindUnpI+6CN4liYwTFBfjkhPLrPoE9wFdA69aMgsrjQ6nV9u4pRGisgFjGv2thkQNwUPzBL9DrM0pNSwO0TULF1kgeRkAF/a70OHrjYodRahPKPgcxRKzgRUrHL6kcBBJZIQQxymqUiNtkwL7O/Vgu1NJCk5SFyGlUlrmBK56hoyVt4c3waY4AHIqig2OpvfSddhQnAvmJgI9QkgOSTYsFm8Pb9LChgCQf2wLNTWGhSLOUGMjLM3LAN8qAj2DnTLC22wyc0HfC9QXUl7m1iSREULY1NbCIdQY1fQoYw/CbK9PsOzIOFNBSR0EqcMmXTGRAYj3V0MIig6XGBtIL1ZRAYe9W8vKYkcY/sHHuCRVXqbHbGHbNkNDEWdg40ZoiFsMwPl9p2A2OWcWsq28rO+PLFnilCWFA8g5Mm6iTuZJuqye9HdTWAiEqR2ZwTGpxgbTjm1HxvcwhWXV6HogLrJZ1CPt2FcCQ8HXFEiYb5jR4XSoX3gyOcdWsbe+GF1HXg8G2LIFl2j0txodO5p3Mt6xNfxPnGh0RKIzli0DUlQmMTVlqtPWvbBvawlbwnq2fXSEyspgIp0/PVx0kyQybmKi/EQWTpCfjy2RsTVUu4BA70BCfUI5XH+YatMuDh4cTESE0VH1XLuOqbKypMC+hn/KfjJD4lNYUAqNfsXs3w/R0UZH1PtkZGBr9DeyP8bK1vAft7m14d81X7vieEuW18O4VQBM7eu8RKZPSB8GhA9g58GdkLKMn366muuuc9rywk6ktEwIYVNQAISp0rLUcNfZkYF2uzJSXuZQhw9DrVfr6OUo1ywrA0iNaD3fJqREXg8GycpugZitgLETy6yGRQ/DhBn8K1m3vdzocEQnNDTAmrK14FlPpE8s6RHpTl3ftivTb5GUl7kp2ZFxYb6+vmRkZBgdhjgDvr6+RofQLXn5DRCnRi+70o4MqE/Ptu3bZmv4HzfO6Ih6JtXorxroU13wMEwrW2N5aDGFhVJGZIQtxcUwuBpPzdslDs/19fRlYNggcg9lUVCzherqBAIDjY5KnMqGDdAYrzKIi1LPd/oO8IV9L+S1ja9ByhKW/ODUpYWdSCLjwjRNw8/Pz+gwRC+SXVEMCS34aAFE+7tWrU5SkOzIOENRERCsxtcaOYXqdFJCW5Os4FIKi1qQAgPnqq+Hknr1QVt6+BA8zZ4GR6SMSxpN7qEsiN1CRsaVnHOO0RGJU1m5EkP6Y6wm95msLiJ2UlBRya5dkSQlOT0M0Q3yk18IYVN4WPXHJAa4zuhlq7bSMjkU05GKi7EdiGr7b+6CEoISVBmRRyO5u/cYHU6vk5cHLZGZAIxJNL4/xmpUjPVgzM1s3WpsLOL0lq09CvEbAef2x1iF+Ya1lbMlrmH1aqeHILpJEhkhBNA6ellXiYwrjV62an8oppwl4zjqMEy1I2P7b+6CPEwehHkkArCzUs6ScbasLCByOwBDo4YYG0w7tob/WBnB7OosFlizewWYLCT5pxr2wcnZiWeri6TVrFljSAiiGySREUIArSVFrY3+6TGul8i0PxRTdmQcp6CkDvwrAdfekQFIClTlZburS4wNpBfKzgaiVCIzKHKQscG0MyJmBBoaBFWwecdeo8MRp7B9O9RGrgDgwtQphsVxdlJrIpO4mrVrDQtDdJEkMkII4PjRywNcbGIZtDsUM7Cc8j3N1NcbG09PVVCpysp8zQGE+oQaHM2pDYhUiczBlmKamw0OppfZmt3QduZU5GCDo2nj7+VP3yA1eCDn8BZ5Xbiw1auBJDV2+dzkyYbFYduRidtERlY9NTWGhSK6QBIZIQRgHb3semfIWEUHRONl9gJTCwRWqBIoYVctLbC7WiUyCQFJLtcn9XPpsckA6MHFlJUZG0tvs233DjBZCPAIJi4wzuhwjjMmcTgATaFZ6gMa4ZJ+Wl0HcZsAmJQ0ybA4+of1J9IvEjwaaYnezMaNhoUiukASGSEEAHn5jRBSAkBqmOvtyJg0E4lBqieC4FLpk3GAykpo8lf9MSnhrl1WBtA3rO0smZISQ0PpVaqqYF9LDgCDIge7XMI7LHqouojOIjPT2FjEyf20cyOYmwj3ijN0QqKmaW3lZdIn43YkkRFCAJBTUQKmFrw0P2ICYowOp0PtD8UsLTU2lp6otBTbxLKUUNdt9LdKCbEmMsWyQ+dE2dnYGv2Hx7pOWZnVsOhh6iI6Uxr+XVRZGez3sZaVTTI8GZ6U2LojlCiJjLuRREYIAUBRVevoZX/XG71sZZuiFVIqiYwD7NqFW4xetrJ9ihtcRlGJNEM4i6s2+lsNjWrdkYnIIyOzydhgRIfa98dM6WtcWZlVW8P/Gtas1WlpMTYe0XmSyAghaGqCSosqJk+LdL2yMqv2h2JKKZH9qR2Z1tHLwa6/IxMbGIsH3mCysL1st9Hh9BrtRy+7UqO/VVJwEv7mIDA3saV0h9HhiA6sXKVD/HoAJiROMDgaNbbbx8MH/A9wWNvJzp1GRyQ6SxIZIQRlZaCHtk4ginW9Rn+r9odiyo6M/bnbjoxJMxHppRKu/ANSW+YsmTn1bRPLolwvkdE0jSGtZ9vs1zI5eNDggMQJlmYUgN8hPDXvtlJAA3mZvTgr7iz1hfTJuBVJZIQQrWfIqMNZUsNdN5Fpfyim7MjYX+muFghW479c+TDM9pICkwEorykxNI7eQtdhW/kOMLUQ6BlCbECs0SF1aERca3lZlDT8u5ojRyDvmNqNGRY1Sk2jdAG2McwJayWRcSOSyAghVCITqsaA9Q3ta2wwp9C+2X/fPp26OmPj6WkK9+0FcxMmTC43UvdkBkarhv8qimloMDiYXqCiAqq91cSyIVGuN7HMqq3hP0sa/l3Mxo3YysomJ48zNph2zopv3ZGJ2ySJjBuRREYIQWFRC4So0hxXTmRs45e9j4HvYVUKJexm1xH1HzTKNx4Pk4fB0XROWkyyuggplteDE2zfjq3Rf2i065WVWdka/mVHxuVs3AgkqERmXILrJDJj4saoi6hscvPrOXLE2HhE50giI4QgZ9ce8GjEhJmEoASjwzkpX09fovyj1BfSJ2NXNTVwVGs9Q8YNRi9b9Q1tHcEcKiOYnSE3F1ujvytOLLOy9sgQsovN26uMDUYcZ92mBojZCsC4eNdJZBKDEonwiwBzM0RlsXmz0RGJzpBERgjBzgOqrCzKq4/LfxJvq8n3r5Q+GTtq3+jvDodhWrWNYC6VRMYJ8vKw7ci4YqO/VahvKDF+6kOZ3IPZNMt0bpextmgrmJsI8Yw09CDMn9M0jdGxo9UXcZvYtMnYeETnSCIjhKC8RiUyfYJTDI7k9IK8g9SF91HZkbGj0lIgxH1GL1vZhhIEVlBYImeGONr2HfUQqgaDuOLo5fZGxqk+mabQLPLzDQ5GALBnD1R6t5aVJY51uR4rW3lZ3GZVAidcniQyQvRyVVVQ46U+yk6Pcd3+GKv2iYzsyNiPu41etoryj2o9S6aFnN1yloyjZe/NA1MLQZ6hxATEGB3OKQ2Lbu2Tic5SvT3CcO0b/c9Ocp2yMivZkXE/ksgIh1m/Hm65BdLTISoKhg+Hhx9G3ny6mOJibBPLBka5USLjVS07MnbUPpFxpx0ZdZaMSryKDsoLwpEOHYLDJnXAZFpEmst9mv5ztob/6Eyys42NRSgbNuCSjf5Wth2ZyO2U7K6nstLYeMTpSSIj7O7IEfjlL2H8ePjoI1VTXVkJmZnw4oswcCC88oo6j0AYr6gI28SylBDXLy0L9ApUF7IjY1elpUCwSgTcaUcGICkoGYDdx0oMjaOny8sDwtWR54OiBxobTCcMte7IRGWTlS3/4LiC1RkHbWeWjY0fa3A0J0oISiDSL1I1/EdnSsO/G5BERthVRQWMGwcffggmE9x2GyxaBFu3wmefwZQp0NgIv/kN3HOPJDOuwF3OkLFqX1pWUaFeT6L7isqPgq+a7uRuicyAaLWDVG0qpbbW4GB6sPaJzICwAcYG0wlpEWmYNQ/wOcK2Yik7NJquw6Y9GwDo4z+QEJ8QgyM6kaZpjI5rKy+TPhnXJ4mMsJu9e1WismMHJCTAypXw3//ChReqsrJrr4WlS2HOHJXkvP02/P73RkctdhbVQVAFACmhrr8jY01kzP5H0XUoKzM4oB6i+JAqKwv0CCXQO9DgaM5MamRrKVxIqezSOVBuLrZEZmCE6+/IeJm9SAnqD0Bxda4cmGqwoiI4Ftx6EGaK65WVWY2JbS0vi5WGf3cgiYywi6NH4bLLID8f+vRRSczEiSc+T9Ng1ix45x319T/+AZ984txYxfFy96hyIh8tkHDfcIOjOT1rIuMXehRA+mTswGKBffUqkUkMcq/dGIBk6+Sy4FI5FNOBcnL1th2ZcNffkQEYGpsOQEtYLjt2GBxML7dhAxC7BYCxCWOMDeYU2nZkNkvDvxuQREZ0m66rnpiMDNXUv2QJJCef+nt+9Sv44x/V9V13QWGho6MUJ1NcpcrK4v36unzzLrQlMt5BKpGRT+C7r6ICWgLU1lZKeKLB0Zw521kUISWS2DrQ9pJK8K1CQ6NfaD+jw+mU9AiVyBCZKw3/Btu4EVsiMyp2lLHBnIJtcllUNnsq6ygvNzYecWqSyIhue/11+Ppr8PKCb76Bfp389+2JJ+Ccc+DYMXjwQemXMYLFAnsbVKN//3DXLysDbGVPHn7VgOzI2MOuXUCQ6iFICna/RMZ2lkxwGaW7WowNpoeqr4ddx9RuTHxAEr6evgZH1DnpkdZEJkcSGYOt3rofgsrR0BgeM9zocE4qISiBKP8oMFkgOlN2ZVycJDKiW4qK2vpcXngBzjqr89/r4aH6ZDw94fvv4auvHBOjOLnycmgJUjsy6bGu3+gPbTsymo/syNiLGr2sdmQSg9wvkYkLjMOEGcxN5JVXGB1Oj5SfD3qYSmTSo9yjrAxgUOQgdRGRK2fJGMhigW37MgBIDhxAgFeAwRGdnKZp7c6TkT4ZVyeJjOiW//f/1Cd1558Pv/512++XHSnjlXWv8K9N/2Jx0WKq6qs6/P4BA9oSodmz1Q874TztJ5b1C3OPHRlrItPsIT0y9lJaim1HJiEowdhgusDD5EGYRzwARQdk+oMjHNfoH+76jf5Wtlj9D7At/4CxwfRiO3dCQ5hKZMYmjjQ4mtMbHt26YxSVJYmMi5NERnTZggWqlMzTE157TTXyl1aVcu8399Lvn/34zcLfcN+393HhBxcS82IMH2V+1OF9Hn0UwsLUJ35z5zr5D9HLqURGlZa5w+hlaEtkGpEdGXtpX1rmjokMQFyA2knafUwSGUdon8i4S6M/gL+XP4mBqvSwtDaXmhqDA+qlMjKw9ceMjnPd/hirYdHD1EV0Jlu3GhuLODVJZID6+noWL17MY489xiWXXMKwYcMYOXIkV111Fa+99hq1cjDBCWpr23Zgfvc7SEuD+bnzGfjqQN7c/CZNLU2cnXg2vxjwC5JDkmmwNPDL+b/kv1v/e8K9AgPhoYfU9dNPy66MMxUV6251hgy0JTK1lqOAzu7d0NxsbEzurqRUd/tEJiVMJTKHmsvk9eAAx50h40aJDMDg6NY+mYhclZAJp8vIAGLUjszIWNffkWlLZLLYv19n715j4xEnJ4kMsGDBAmbNmsX8+fPx8PBg6tSpjBkzhvLycubMmcO1117LoUOHjA7TpTz/vPokPDER/vQn+Dznc6777DoaLA2c0+cclv9qOavuWMXXM76m8NeF3DfmPnR07vjqDt7e/PYJ93vwQQgJUf9Yfvut8/88vdXOXYfAWzXN2yY/ubhAL9Xs39TShIdPAxYL7NljcFBurnRvFXipj6rjg+INjqZrUqNVIqMHllEhbTJ2l5NrgbACwP0SGdvksgiZXGaUjZlHIFy9fkbGuH4iMyB8AF5mL/XvY0ip7Mq4MElkAE9PT2644Qa+++47vvnmG15++WXefvttfvjhBwYNGkRxcTHPPvus0WG6jP374cUX1fU//gHr9y/lpnk3YdEt3Dr8VpbeupRz+pxje75JM/HaZa/x4NgH0dG555t7yNqXddw9g4Lg7rvV9SuvOOtPIvIPqN2YUI84fDx8DI6mc9o3icYlqyRMzg7pnt1H1W5MiFc4fp5+BkfTNbZpa0G7pW/KzlpaIG/vLvBowNPkRVKwe5011H4EszT8O5+uw5aKbQDE+CYR7uf655V5mj3bXjdSXubSJJEBrr76ap588klSUo5vdo6IiODxxx9H13V+/PFHmqVeAYC//x3q6mDsWBh5fhHXzL2GppYmrht0Hf+58j+YTeYTvkfTNF655BWuGngVOjovr3v5hOc88ACYzbB0KWRmOuNPIsqOqUQmMcA9Gv0BzCazLZmJ6aP6ZCSR6bqaGjiKSmQSg92zrAzalcQFlcnrwc7KyqAxQJWVpYb17/BnvCuzTS6TEcyGKCuD6gDVH3NWguvvxlhJn4x7kETmNNLS0gBobGykqqrjyVu9yf79qrEf4E+PN3HzFzdxtOEoExIm8P4175/yHzhN03hk0iMAfJT1EXuPHV90mpQE06ap6zfecEj4op2mJjikq0b/1Ej3SWSgrU8mMkESme4qK6PtDJkQ901kEq07MsFlsiNjZ/n5tE0si3CfiWVWtrNkgsvI2nHM2GB6ofb9MWPiXb/R30oSGfcgicxplJWpCTgeHh4EBwcbHI3x2u/GbPF/jvXl6wnxCeHj6R93qjRpfMJ4JiRMoNHSyD/X//OEx++5R/3vJ5+osc7CcXbvBj24BID0GPdMZEJiVCJTJoOqumz3bty+0R/anX8TsJfiXY3GBtPDtE9k3K0/BiDMN4xI3ygAyhvyOCa5jFO1n1jmDv0xVu0TmZ07kYl3LkoSmdN47733ADjnnHPw9PQ0OBpjtd+NueMPOfx15TMAvH7Z620na3fC7LNnA/DP9f+ksqbyuMfOO08NEKiqgq+/tk/comMlJUBICQApockGRnLmrA3/wZGyI9NdZWW49WGYVpH+kXjgBZpOvkx/sCt3T2QABkW1Nfzn5RkbS2+zaWsdRKpxcaNi3XBHJiwf3aNWyhJdlCQyp7B8+XLmzZuHp6cnv25/2mMv9eKLajfmrLE6H1bdQ1NLE1cMuIIbh9x4Rve5auBVjI4dTU1TDS+sfuG4x0wmuO02df3uu/aKXHSkfSLjLhPLrKw7MnpwKXhVSyLTDT1lR8akmYj0VvGXVskWnT3l5wMROwD3TWTaN/xLIuNcG3dlgclCiGckcYFxRofTadH+0UT6RYKpBSJzpLzMRUkicxKFhYU8/PDDAPzhD39g4ED3qwu2pyNH4F//UtcXzJrPqrJV+Hr48tplr6Fp2hndS9M0njrvKQD+tflfWFqOPzjm1lvV//74I1RW/vy7hb2UlOoQrJoJ3DWReaP0Qbh3BLt2yyCOrmrfI+POiQxAfGv8e2rK0HWDg+lB8grqIVh9WuC2iUyknCVjhIMHYb+mJpaNiht5xu8XjKRpmvTJuAFJZDqwb98+7rrrLqqrq7n99tu55ZZbzuj7LRYLTU1NHf5yV2+9BdXVkD64ic+rVMP+7yb8rstjOC/udzEBXgEcazxG3oHjPx5LTYXRo9XBmPPmdTt0cRJ5ZfvAsx4Nk9u9gbUmMgCEFXG4rkrq3ruoJyUy/SNVaVyDTxmHDxscTA/R3AzFh0tA0wnwDFSfULuhth2ZHElknCgjA4hWxy2MjBtqbDBdIIlM15zsPXBTUxMWO5967mHXu/UAR44c4Y477mDPnj1Mnz6d2bNnn/E9Xn/9dV599VUHRGeMxsa2s11G3f0mHx3KJ8o/ytbr0hVmk5mRMSNZuWslmyo2MThq8HGP33ADbN4Mc+fCvfd2J3pxMvmVJZAEYR7x6uAvN3JcIgPgVU1ZWQTp6cbE485K9x6xHYrq7olMcujxZ8mEhRkbT09QWgqWoEIA+oX1datP1NtLi1ATSA3E0/YAACAASURBVAktIndpE9C7e16dJSMDiFKJzNAo90tkbDFHZZG5Qn3Aanav6eOGGDt27EkfmzVrFg8++KDd1pJEpp3a2lpmzpxJUVERF110EU8//XSX7nP//fdzbwfvvmtra0/5l+uqPvkEysshus8RFtY/CcAT5z5BoHdgt+47OnY0K3etZPOezdw24rbjHrv+epg9G5YvV6e2x8Z2aynRgd3VqqwsISDZ2EC6QP953ZB3NWVlSCLTBWWth2EGeYbi7+VvcDTdYxvB3HqWzEj3GZDksvLzgVB13lTf0L7GBtMN8UHx+Jh9qaeO/MoSmppS6eXze5xiS4YOia2JTLT7JTK2HZmYbdTW6hQUaPTyToNO2bBhA35+HR+ubDLZtxhMSstaNTY2ct9995Gdnc3kyZN56aWXuvzJk9lsxtPTs8Nf7kbXVZM/wKCZf+dA3QHSItKYOWpmt+89Om40AJv3bD7hsT59YPx4tf78+d1eSvxMczMcsJQA0D8i2dBYuqKq4WdnOknDf5dUV8Mxzf0Pw7Sy7SjJWTJ20z6R6Rfaz9hgusGkmRgYofp7LCE7KSw0OKBeYvOOveB3EBOmtvI+NzIochAmzQR+ByFwD9u2GR2RezjZe2BPT0/Mdt7SkkQGaGlp4aGHHmL9+vWMGTOGOXPm4OEhm1UAixZBVhb4hR9io0nVlz17/rN4mruflI2JGwPA1r1bT2j4B7j6avW/MobZ/srLQQ8qASA9NtnQWLriT5P/xJCoIW2/4X1UEpkuaD+xzJ0Pw7SyjY9u3ZER3ZefD4Spd/3uvCMD7QYVhO+UyWVO0NAAhdVqNyYlOBVfT1+DIzpzvp6+9A/rr76IyFWlcsKlyLt14MMPP2Tx4sVomkZISAhPPPFEh8/7wx/+QEhIiHODM5i11WfIzH+yoekYw6KHcXXa1Xa594DwAcc1/P+8T+aqq+CRR2DpUjh6FIKCTnIjccZKS3HbM2RAnS6edV8W5/73XFaUrgBv2ZHpit27sZ0h4+79MdCutCxgP0XbGwBvQ+PpCfLzgb7uX1oGxycyubltH5YJx8jNhZYIdfjKiLghp3m260qLSGPnwZ0QkUdW1lSjwxE/I4kMcPToUVsZ2eLFizt8jqZpPPjgg70qkSkthW+/BbyPkhf0CjSpT8Lt1exp0kyMih3FitIVrChdcUIiM3CgmmCWn692hq691i7LCtz7DJn2rAdjqmZ/Y2NxRz1pYhlAuG84XpoPjXo9hfvLAfd+4+0KduzUYXRraVmY+5aWQbtEJmKHTC5zgqws3LrR3yotPI2v+Roic9WfSbgUKS1DTVDIzc095a+cnBzi4tznICd7eOst1aPS98ZXOdpURXpEOtMHTbfrGlekXgHApzmfnvCYpqldGZDyMnsrLtZ7RiJjHTghOzJd0j6RsZVluTFN04jxU3+OsiOS2XZXYyOUHNgLnnWYNFOXx+27ip/vyAjHyszENnrZHRv9rdrOIMpj1y51rp5wHZLIiA41NsK//w14HaOy/z8A+OPkP6qmNzu6fvD1ACwvWU5FdcUJj1+h8hwWLYKWFrsu3avt2L3fbc+Qac+2I+N9lLIy5BDEM9S+R8adXwftWXt9DreUUV9vcDBurrgY9GC1G5MYlOh2Y9p/zpbIBJWTW3hMfl44WGaWBSK3A26+I9M6utscrRqrZFfGtUgiIzo0fz7s3w+BU96m2nKQ/mH9uWHIDXZfp09IHyYkTEBH57Ptn53w+IQJ4O8P+/bJDw972llZAkCo2f3OkGnPdp6MdzUNDVBZaWw87qanlZYB9A1vO0tGyg2757iJZW5eVgYQ5htGhG8EADXeBZSXGxxQD5dRWgie9XibfN26v2pguJq3bPEvB69qeS/iYiSR+f/svXd4XHeZ9v850zSSZjTqXbLlXmJbLqkmm04ghCWhE94UwhLyprCwu2STvCxZYCm/8BJe2BTqkgABNgsbIGSXDU5wAoTguKlbtnq1rG5JozKamd8f3zlnRsWyyswczTnfz3Xl0tFI1nmkjEbnPvdzP49kXp56ClD8WC/9JgCfvuzT2CyxiVSprsx/NfzXnI85HHDlleL4pZdicnpT0j7SAiTmDplIVEcmJV0sdJTtZUuj9fQIOEWfhFGEzOxdMpLlI4RMMwBl6WX6FhMlNmXL9rJ40N8PvYq44t+Wsw2rJXG3SGYkZ5CXmifeya6XQmaVIYWMZA61tWIRpbLl1wwpLWQlZ3Hrzltjdj71Ts3QxNC8H7/uOvH2d7+LWQmmwu+HvukWIDF3yESiZmSSpZBZFu1D4pa0y5624gW3qwUt6yN3yayYU6cwRJYuknBOpl6OYI4hIugvJpbtKkjctjIVtb1MTC7TtxbJTKSQkczhW98SbzPf8XUAPr734zGd/55sE1973Dc+78ff+lbx9rXXYHz+T5Esge5uCLhbANhSsEbfYlaI6sg4XELIyFaixXP2LIxZRC6tyF2kczXRQ3OWpCOzYgwpZDKlIxMPZgT9EzgfozJbyMh81epBChnJDMbG4JlngIIj9Lv+gM1i454L74npOVPsKQB4fV46z3byyO8fmRH837IFiovFcq0//jGmpZiCyNHL6xJwh0wkqotgTTkLSEdmKXR0AO6QkPEU6FtMFNFay6Qjs2JOngQ84odoGCEjJ5fFBaOMXlbZmi0mlyk5dQwPy5tmqwkpZCQz+OlPxZ1a17XfAOAD2z9AUVps79aqbs/49DhPvvkkn3/t83zr8Le0jytKuL1M5mRWTuQyzES/OIkM+4MUMkuhvR1NyBS6jTNaXnNkUvppaZ/Ut5gEZmIC2jr84BG/VGs8ie3eqmzOFsFtsk5Sd0LeVo8Vx6vHIbMBgAtyE3cZporqyDgK5eSy1YYUMhKNYDAU8nd1M77+ZwB88pJPxvy8qiMz7htnYHwAQHuroraXyZzMyhE7ZIxxl1VtLfNbpZBZKpGOTKHLOEImw5mBw5IEQFNvt87VJC6NjYCrC6zT2Cw2w4jd9RnrUVAgeYies30MzR/NlKyAQACqe2rBEiDdkUW+K1/vklaMKmR8aafAMi2FzCpCChmJxptvwtGjYL30Sfz42F+yn32F+2J+XjUj4/V5GfONATAxPXMBxDXXiLcVFXD6dMxLMjQn2nvBPg4o4TacBEVtLZtShJDp6NCzmsTCqI6MoijkpYjvp2u0U+6fWiaR+ZhST2lCT52KJNmeHF7smV0v28tiQFMTTHjElf6ugh0oiqJzRSunxFNCsi2ZgOKD9GaRAZKsCqSQkWg89RRgG8d2iWjriocbA2FHZtI/yejUKCDazCLJyYE9e8TxgQNxKcuwnArtkMlM8B0yEHZkvH6RkenuBp9Pz4oSB6MKGYDSdNEOO+3soqdH52ISFCFkjOHcziYyJyMnl0WfyIllO/MSPx8DYFEs4bZEOblsVSGFjASAgQH42c+AHT9h0trHGs8abtpyU1zOHTkRrX+8H5h/gtm114q3Bw/Goyrj0na2BYDC1MTveVczMuPTXqx2P8GgdOwWy4zWMoMJmSJP6Ptxd8l2w2US6cgYJR+jIgP/saWyEkMF/VXUwD/ZJzhxAqam9K1HIpBCRgKISWUTE0GcV/4/AO6/6P6YLcCcjdpaBtDn7QPmOjIAl18u3srJZcsnEIDeKXGXdX1W4l+cRO4+KVwj3DzZXrY42tqDhhUyWubH3SUnly0TI45eVpFCJrZUVaGNXjZC0F9FzcnYC+uYnob6ep0LkgBSyEgQIf9vfQsoe4UJTzWp9lQ+uuejcTu/1WLVWpw0ITOPI3PZZeJtfT309satPEPR0wN+t7iyS/QdMgBJ1iRNcOeVypzMUujoGwSbmOpV4DbO+GWIEGbuLjkmdZkYcfSyyobMDeIgo1G2lsWAYyf6wS0GbRhRyDiL5eSy1YQUMhJeeUX80bK+RSzA/Ej5R0h3pse1BtWV6feGWsvmcWQyM2HbNnH8+utxK81QtLSgXZyUZSS+kFEURcvJZBeJnIwUMudneBhGFeHGZDozcdqcOlcUXbSR8WmdsrVsGYyNQVcXhnVk1mesFweZjTQ2BZmYWPjzJYvH64WmMXGFX+JeO8M1T3RUITOVJoRMba2e1UhUpJCRiJB/1kn8618E4BMXfyLuNaiBf3/QD4SnlgVnrc99y1vEW9letjzEMkwhZNakJ76QgXBOJjNfOjKLZUY+Js1YbWUw05GRQmbpNDQASiD8WmGwjMza9LViBLNjjGDyGTFqWhIVamrQ8jHlBcbJxwCsy1gHwKRlEJKGxfcq0R0pZExOVxf88pfARf8KwDs2voONWRvjXkdk4B9Ea9nnX/08ef83j6PdR7XHpZBZGa2taI6MUS5O1Dt+r6XdCxv+m85OnQtKAIw8sQykkFkpTU2AqxusPqyKNeZLkeNNki0pPHo+s1FmHaJIVRWQI6yK7Tnb9S0myrgcLrKSs8Q76a3SkVklSCFjcr73PfDbhrHsfRqI38jl2aiOjMr49DiPHHyEXm8ve7+zV3Nm9u8XHz9yBMbndp9JzkN96xA4RQuWtkshwVFbyzr8R+BDfy0dmUVgGiGTNELr6RF9i0lAmprQ3JgST0ncBr/EE629LEMKmWhSWQlknQTCrVhGQmuz9LTS0ACTk7qWI0EKGVPj88F3vgPs/jcCtlG252znmrJrdKklcnIZCEcmw5mhvf+2Z9/GQwceYs3aAAUFovY334x3lYlPfY+4OHFZskl1pOpcTXQYmohYzW2dpr1DbkA8HzNay1zGEzIuhwu3Q7Qc9k92yZseS0QImRbAePkYlcicjBQy0aOqCsgWP1BtOpyBUFuykwtaCATk5LLVgBQyJuZXv4LOLj+WS0Vb2Scu/oRuG3jntJZNj+NxerT3X2p8ia/86Sv8oe012V62AtqGhZApTDFGWxkwZ6ln59AZuc39PBjdkQEoSpOTy5ZLczOG3SGjsj5TOjKxoPLEKKSJ/l5tgaSBUH8fMsrE31LZXqY/UsiYmMcfBzb9hoCnmczkTP7Xzv+lWy1zWst844xMipaQO8vv5Io1VwDw48ofSyGzTIJB6Jk0zsQylceuf4wbNt6gve93tXLmjI4FJQBGXoapInMyy8dcjkyDGDUtWTG9vdAXOAVAVnI2mcmZOlcUfdTfB0duCyCFzGpAChmTUlUFr74KXPINAD6252NzxEQ8md1aFiSotQz90xX/xOeu/BwAP6/9OfsuERPNXn8deed9CZw5A9Mp4opuU74x8jEAV5ddzYu3vMhlJaFFQ+mtMidzHszgyMhdMssjEAg5MgbdIaOiOTKZjQwMQF+fvvUYgdpatLayzdnGayuDsCOj7mOTk8v0RwoZk/LEE0BeJZT9Hqti5d4L79W1nvlElDqK2eVwcfmayylJK2F4cpj25N+Qmip2YVRXx7vSxKW1FS3Auy7TOI6MitYC42mTQmYBgkFEjsglFtYZVsio2R+5S2ZJdHXB1BSQ0QIYWMiojkxqLzhGZHtZFKitBbJCQibLeG1lEM7InFVka9lqQQoZEzI0BD/6EXDxNwF499Z3h0dR6sRsRyaSVHsqFsWitb49Xfl9Lr5YfOwvf4lHdcYgchmmEfvew0JGOjILcfYsjAX7wDoNQL4rX+eKYoM2Mli2li2JpiZACaAYdIeMisfpCY/SlYH/qDDDkTGokFGF/fB0L9jHOHUqJPwluiGFjAl55hnw0oey61kA/vbiv9W5ovkdGQCLYtG2jn9090cB+J+G/2HzxS0AHDoUl/IMgRGXYUaijZOWrWUL0tmJ1laWm5qL3WrXt6AYITMyy6OpCUjtIWidxKJYKE4r1rukmCED/9GlpgZt9LIRJ5YBpDvTtSXMqUVt+P3IjJXOSCFjMgKBUFvZRY8TtE6wr3BfOFugI7Onlqm4HC5tktr6zPVcu+5aggQ5U/JdQAqZpdDYOgGuHsCYd1k1cSZbyxYkUsgYta0MpJBZLmJimbjhUZxWbFihC3IEc7SpqQ2GW8sMOLFMRXVlii9oAWR7md5IIWMyfvc7ONUyhnLx4wB8+rJP6zZyOZJztZal2mfuOvn43o8D8KexZ4Ag1dUwNhbr6ozBidPiai5JSTXkNBnNkZGtZQtiSiHTHiS0U1dyHswwsUwlcimmvKu+Mvr74Yy3G5JGsSiW8M/WgJSllwGQsa4JkEJGb6SQMRmPPw7s/gHB5H7WZazjPVvfo3dJwLlby1wO14z3b9x0I8m2ZE6PdZK7vY5AAI4ciUeFiU/LoLjLmu9csyrEa7TRXKbkIdp6zupbzCpmhpAx4DJMlQJXgTiwTTLBAP39+taTKJhKyERMLmtogOlpfetJZOrq0NrKytLLSLIl6VtQDNmQuQEAe54YNS0nl+mLFDIm4uRJ+M1/TcOlXwPgHy79B6wWq85VCc7VWjZ7+7zT5uTyNZcDkHfZ7wDZXrYYgkE4PW7cfAyAO8mNx5EBQMdoq7wDfw7MsEMGIMmWRHZKtnhHtpctmkghY8QW1EhU10DJbMTnC+UIJctC5GNEW5lR8zEqGzM3AjCRKoSMdGT0RQoZE/HYY8C2n0NGCzkpOdxRfofeJWks1pEBuLbsWgB8JQcAKWQWQ38/TCULIbM5z7gXJ2tDIs2X0sLAgM7FrFLM0loGMiezVLxeOH0aw++QUVEdmaCnDaxTMiezAswwsUxlY5YQMuryz5MnwefTsyJzI4WMSThzBn7wdBD2/38A3H/R/ed0QfRgsRkZgGvXCSHTqhwEi08KmUUQObFsXZZxhcz6rHXiIKNZ5mTOgSmFTFqnXIq5CJqbxVtLVgtgfCFT4CoQf3uUAHhapZBZATN2yBg46A9hR6Z9tJlU9zTT03DqlM5FmRgpZEzCE0/AVNEBKDhOij2Fey68R++SZrAUR2ZX/i6yU7IZ949C8V9obYWenlhXmNi0tmLoHTIqagiTdClkzoWphIxLOjJLoakJIEgwzRyOjKIorMsI3fzIlIH/lWCG0csqRWlFOG1OpgPTrN/bAsj2Mj2RQsYEeL2hkcv7HwXgb3b/DVkpWfoWNYuFxi/PxqJYuKbsGgCyLxLtZW++GbvajIDRd8ioaEImo0kKmXmYmoKe3mlIFcrf6EJGLsVcGs3NQGovQds4Coqhd8ioyF0yK2dwELrPTEGGsPSM3lpmUSxa4D93q8zJ6I0UMibg6aeh33EU1h/Aqlj51KWf0rukOSyltQzC7WWWDTInsxiaWqYhTVzZG9mR0e6uytayeenuBlLPgCWARbGQm5qrd0kxRWZklkZk0L8orQiH1aFrPfFgrWetOPC0SSGzTOrqgIxGsPhxOVyGv0EC4fay1BIpZPRGChmD4/PBV78KXPE5AD54wQdXZbvAUlrLAK5bdx0Afc43IOmsFDLnob6rCyx+LNjId+XrXU7MKMsIOzLtHXJs2Wwi28ryXfmrZmphrNCEWkqvzMgsAjONXlYp8ZSIA0873d1wVk5uXzIiHxNuKzPieP/ZqELGny5HMOuNFDIG54c/hJapw7Dl11gUC5/5q8/oXdK8RLaWOW1O7Xj2+GWVNelr2JC5gQB+WPMqhw4hx+0uQPOAaCvLTSox9MWrdvHlGKP5dJ+utaxGzJSPAXA73OLAMUpXl5wsdD5MKWTShJCxZwkLV+Zklk5NDdrEMqPnY1TUyWUj9gYA6uvlHiK9kELGwPh88MUvAlc+AsCHd3yYLdlb9C3qHEQ6MlnJ4fzOuRwZCLsylo2/Y3BQDapKZhMMQpdXCJlSA7eVgRDBWQ5xgd56tlnnalYfZhMy2utH0gjBYOj7l8xLMBh6DTXBUJBI1ByQJUNYdlLILJ0ZE8sMno9RyUnJAWCCQVJTxfVWY6PORZkUKWQMzI9+BM2+N2DTf2FVrHz2is/qXdI5iczIRA4iWEjIXF12NQBJm/4AwLFjMSouwRkagokkcXGyKdf4FydrPSInc3qySbp0s+jsBFzdQMRELwPjThKOjCV5BEDmZBagpwfGxwH3aUCMJjYDamuZz9kJSkDmZJaBmXbIqKjdImO+MbZuFY/J9jJ9kELGoPh88C//gubG3L7rdm3Kxmok0pHRtnFz7rA/wGUllwEwkVYJjhEpZM5BSwvgEVdw6zKNL2Q25YiczGRys+x3n0VHB6ZyZMKtZVLInA/V0U7K6AUw/CAIlQJXAQoKAcUHqWekkFkiZ8+GXldMMnpZRb02GZ0aZds28ZgM/OuDFDIG5Uc/gubgQdjwEjaLbdVmY1QiMzKLbS0rdBeyxrOGoBKAokNSyJwDs4xeVtmYrQb+5eSy2ZittUx1ZAKWSbD4ZOB/AdRlmFa3EDI5qTk6VhM/7FY7Be6Q+5TWLoXMEqmtBZyDkCqeN6YRMqojMzXG9u3iMSlk9EEKGQPi9cI/PeKH68WY5bv23BWe5rRKsVls2Cw2YKaQOVfYX0V1ZSh5naNHY1ZeQmOWZZgqmljztEkhMwuzCZkZN0Ico9KRWQDVkfE7zwDhDIAZ0PblpHVw8iQEAvrWk0hEtpUVugu1mwdGR3VkxnxjmiMjW8v0QQoZA/K1r0FX9o+g4DieJA//fOU/613SolD/cGrjMFnYkYEIIVP6Oj09oT0Zkhk0twRN5chEXpRIIRNGC7ubSMg4rI7wLpSkESlkFqCpCbBMM2kdAMzjyEB4cpklox2vVw6FWAqzRy+bBfUmq9fnZctWoXzl5DJ9kELGYHR3w1ceG4VrHgbgM3/1mYT5g/TT9/yUZ9/97Iwsz2KFjKX0z6AEZHvZPJzs6AP7OBD+g21ktO8xrV1ekETQ3w+TPh+4xB13MwgZmJmTkULm3DQ1Acn9ACgoM5xxo6O+ZnhK5OSypVJTg+kmlsHM/G5e0TjJyTA5GW7RlMQPKWQMxmc/C97yr4K7m7L0Mu6/6H69S1o0V6y9glt23DJjgtlCYX+AnXk7SbGnEHAMQ3adbC+bh6Z+4cZk2gtIsiXpXE3s0Rw951maOmXaX0VMLBMTqewW+4zpgEZGa3VxjMqMzAI0NaHlHLJSsgy9b2o2qoubnC8sXJmTWTxmnFgGMwcUjfvl5DI9kULGQFRUwPd/VQ/7HwXg0eseTcgL18jg//kcGZvFxsVFF4t3Sl6Xjsw8dIwKIVPiNn5bGYjnTIqSDkBjn7xyVYlsKytwF2BRzPHyH7lLZngYhof1rWc1MjkZen6kmi8fA+GbH4pHvF5IIbM4RkZCkwBVRybbPELGarFqy7vHpsbk5DIdMcdfMhMwNQW3f2Sa4LtuB/sE1627jvdsfY/eZS2LGY7MecL+MDPwL4XMTIaHwWsXQmZDjjmEDECuU9xh7TgrhYzKDCFjkh0hEG4tc2WKEczSlZlLa6vIUCVlmmv0soraWjbuEL2HUsgsjhMnACUAWacAc2VkYGbgX04u0w8pZAzCv/wLVKR+FYr/gtuRxvf/+vsoiqJ3WctCdWQUlBmi5lxECpnmZhgcjGV1iUXkxLL1WeYRMuod1jMTMu2vIlrLegDId+XrW0wcUVvLMgvkLplzoU4syyox1+hlFXWq56C/A6xTUsgskpoaxI4y2yR2i5216Wv1LimuRI5glpPL9EMKGQNw+DB88XtVcJVYfvmvb//mjMlfiYbae+pOci9KjF1SfIk4yD4JKX0cPx7L6hKLyB0ypZ5SXWuJJ+tzxPPfa2tnbEznYlYJHR1AqhAyeal5+hYTR9TWMk/OKCCFzGyCwSB1jeKXxJVnztayvNQ8km3JBAmCp5XWVpiY0Luq1Y+YWCZU34bMDdoKBbOgvrZEjmA+cQL8fh2LMiFSyCQ4Xi/c+tFRAjd/CKw+3rnpndy26za9y1oRm7I2cfuu2/nM5Ytb4pmZnMnW7FDSrvgN2V4WQUsL4R0yJhi9rLI+OyTkPXJymUpkBsJMrUOytWxh7nrhLv6uzwVZ9ThDrWVmEzKKomiuTEpRM8EgNDToXFQCYNbRyypaa9nUGGVlkJQkBHBLi751mQ0pZBKYYBA++jdBTmy+E3JryE3J5zvv/E7CtpSpWBQLT9/0NJ/e/+lF/5tLiy8VB8V/lkImgtZWwjtkTLAMUyVyBLPcJSOIbC3Lc5nHkVGFTHK6bC2bj+8d+5442P8oFrc5MzIA6zLWAZC7WfTZyfay82PWiWUqWmuZbwyrFbZsEY/X1elYlAmRQibE5OQk3/jGN7j++uvZuXMnl19+OQ8//DA9PT16l3ZOHnsMftb2Vdj+H9gUO//5gZ+bqvc9kn2F+8RBwVE5gjmCU20jkCxCQ2ZyZORSzLmY1pEJZWQcLilkFiRphOkkc2ZkAMrShSOTWiQWgUghszBjY6GdKSacWKYS6cgAcnKZTkghA0xNTXHbbbfx1FNPMT4+zrXXXkthYSH/+Z//ybvf/W46VuGV0O9+B5/+zn/DNQ8B8M23f4P9pft1rko/9hTsEQeFRzhRH5T9zSEaeoUbk2pNJy0pTedq4oeWEfO0094e1LeYVcD4OAwMYOqMjDVFZmQWJOksXsyZkYGwI0OmdGQWw4kT4q0lVzoylT2V3PXCXeRtaQSkkIk35kpmnYMnnniCiooK9uzZw/e//32Sk8WkrKeffpqvfOUrPPzww/zwhz/Uucowhw/DzZ/6PcH3vRssAT5Sfid377tb77J0ZWfeTqyKFX9qL4HUTurqitm9W++q9EfdIVPsMo8bAxGOjGOMpq4hIEPXevRGywm5TOjIONSFmMKR6ewUYVyrefY9Lo6kswxPy9ay8SThyJw8qWc1q5/aWsDuJeAWoTMzZ2S+eeibAJQmvQFUSiETZ0zvyPh8Pn7yk5+gKAqf/exnNREDcMcdd7B582befPNNalfJM/PPf4Yr7/kFYzfdAPYJblh/I9+68amEz8WslGR7MttyQr5uwVEqK/WtZzUwMgIjFvONXgYx+S5VEZvr5VLMkJCxToJzCDBZRibUWjZtGcFqBZ8PVnHHsG7Y3AMMjA8A5m4t6/OHHZmgNHPPiRi9LP6+C1atOAAAIABJREFUuB1uslOy9S1IB1Qho9I2WQWIjIx87sQP0wuZo0ePMjIyQmlpKVvUpFYE119/PQCvvPJKvEubw09+EuTyB77O2A3vA/sEb1/3Tn7xwf/AYXXoXdqqQGsvk0IGmLVDJttcQgYgzynay9qHV19raLyJzMfYLDYynOZxqFRHZtQ3QlGReEy2lwn8gfCc2On0eoIESbYlm7K1TJ1adtY3CM4hBgehr0/nolYxtbVoC3aL04pNeTN19sJum8WGzQajo8hsZhwxvZA5EWr03KamtGaxfft2gsEgJ3X0mf1++D+f8/Lhn30c/7V/B0qQu8rv5YUPP4/T5tStrtXG3oK94qDgiBQyqKOXxRWbmSaWqajtZb2T0pGZHfQ300WHmpEZnRqlNLRKSQoZwfj0+JzH1qavNdXzQ8XlcGkCrmCrDPyfj9paIE30rBa6C/UtRidmOzKp9lQ2hTrsVkkTjykwvZDp7u4GID9//mlfeXmiBaNTp2UUFVV+tn/gp3ypvxz2fheAr177f/nWX/8rVots8o5EOjIzmTF62UQTy1TUpZgjlnYmJ3UuRmdmjF42UdAfwq1lI1MjFJUIB0LukhGM++YKGS30bkLUqZ+FG4QVI4XM/IyPQ1MTmiNTlFakb0E6MduRcTlcbA2ttJNCJn6YXsh4vV4URcHpnN/ZSEkRW+bH4rgePBCAF/5nhH13fZ/yb19A/Y5bIOsUHms+v7v1d/zD/r835R2z87ErfxcKCqR1ccZ72vR98DOWYZrQkdmQG94l09Wlby1609GBKUcvQ7i1rGmwif/Y5IQtz0tHJoTX553zmJoVMSPqhWleifh7L4XM/Jw4ITIgzhxxg7fIbVIhM9uRcaTKEcw6IKeW6czh1loO1h+hq3+U0wMj1HW1c2KwgomsQ1AkbiPbpzO4f9/f8dnr78fj9Ohc8erF5XCxJXsLdX11IVfmBq67Tu+q9KOpdQq2CcfRjI5MacQI5o4OKDPv9Vmotcx8yzAh7MgABJRpuO4B2tpu1rGi1cO8QibDvL8o6oVpdqEQMnJy2fyoF+mpBZ1MYOLWslmOjNPm1ISMXIoZP0wvZFJSUggGg0ycY/GI1yte6FNTU+f9+Hz4/X4CgcCcx30+n3ac+cDF+Kw+Apmzbvk4gNB1Rvr0Ju4o/yifu/FuU+0AWQl7CvZIIRPi1Jl22B7EoThNGd6NXIpp9laizk5ga8iRSTGXI6NmZDT6N0lHJsTIxNzWMunIQEaudGQWQhUy1vRQa5l0ZACYmJ6Y4cgEg2DW5hmfzzfjmjcSi8WCNYrz700vZAoKCgA4ffr0vB/vCfUnFRUt/hf1ySef5PHHH1/wc6YyagjaguC3YencTzJZpNpdFLjyuXj9Vv73jZexq3ijbCFbInsK9vBs1bMyJwO0DYu2soKUUlM+j0rSwq1lHR1BwHw/AxDDQrq7gQtN6sg43DMf8KVKIROipVM6MpGoF6Zp2ULINDbC9DTYTH+lNBNVyEw5OyEgMzIqY1NjbNoEFgsMDoox7+eIXxueiy666Jwfu++++7j//vujdi7T/3qqI5fPtSempqYGgE2bFr/s6Z577uHuu+cuqPR6vdr/3C9d8ALZuTau3FLO+vw806r2aDNjctlBXUvRlbExGEYImXUm2yGjojky9gkauvoB8+05ADhzRlyMqWF/s2VkkmxJMx9wjNDXJwLLEWvDTElTu8zIRKIKGXvKGMnJ4jnS3AwbN+pc2CqjpgZQApwNiNZl07aWzXJkxnxjOJ2wbh00NAjBZ1Yhc+jQIS1jPhuLJbrxfNOH/ffs2YPb7aatrU0bxRzJb3/7WxRF4eqrr17017Rardjt9nn/U/nEO6/ib666ng0FUsREk/L8cnGQ3kZNcx/ncDYNT1sb2sQysy3DVEmyJeFWxEV7Y695e8vUgYu2dOE6F7gKdKxGHy4pvkQ7tiSPAnJyGUBb98zWsgxnhqlzmOoddq9vjM2bxWPzXBaYlkAwwHcO/YCGoROQeoYAfiyKRZv2Zjbmc2QAGfiHc14D2+32qLaVQZSFzNatW3n44YfP+3mf+cxnzrm3Jd7Y7XY+/OEPEwwG+fznP8/4ePiF/Qc/+AEnT57koosuWjX1ShbG4/SwIXMDAL6sY6YNa86YWGbCoL9KbmgpZseIebeTaZPjXeLuaYHbfELm5dte5tl3PwuAPXUEkLtkADp6ZjoyZh69DOE77GMRQkbmZML8R81/8PH/vpPgvVtxFYoXlrzUPGwWczb3zHZkfAEfPr9PBv7jTFSffcFgkGAwuOjPXS3cc889vPHGGxw7doy3vvWt7Nu3j66uLioqKsjOzuaLX/yi3iVKlsDegr00DDRA4REqK69j+3a9K4o/Zl+GqVKUVkzj+BF6xs17+72zE7BOMe3oBzDl3dMUewpr09cCYA05MlLIQHevFyLiDZqjbVLUO+xjU2OEus6lkIng2Olj2nHJ1i7qMG9bGcwzSAQhgrdtSwfM7cjEE11ay0ZGRnA4HHqcel4cDgc//OEPueeee0hJSeHll1+mq6uL97znPfziF7+guLhY7xIlS0AuxpTLMFXWZwtHZph2kRMxIZGjl+0WO5nJmfoWpBPqRUfQLhwZ2VoGZwZFB8Kl6e/hldte4Wtv/ZrOFenLfI6MbC0LE/na4SlrBMwb9Ie5rWUgRLBcihlfVuzIdM3aNOf1euc8puL3+2lqauJPf/oTpaWlKz11VHE4HNx///1RnaQg0Ydw4N+8Qqa5JQAbxZWamR2ZTXklUAOktdPTA0sYPmgYOjoAt2gry3PlYVHMGY1Up5dNW6Ujo9I7JFrLcjwuriq7Sudq9EdzZHzSkZkPuyWc8x3PeR185h29DHNby2Dmc+fMGejvh6ysOBdmMlYsZK6++uoZo11feuklXnrppQX/TTAY5H3ve99KTy2RzMvugt3iILOR4yeGgHRd69GDU92nYesUFqymvmO2Jj1yBLM5hUxnJ+Ayb9BfRXVkfHhB8dPWFt3AaaIxNATj0yEhk27y8W0hNEcmNEYXoLcXBgYg05xG5gxGp0a146bgywBsylr8RFejYbXMfQ0ZmxrDlQVr1ojOiLo6eMtbdCjORKxYyFx44YXa8ZtvvklWVhZl51ih7XA4yM3N5eqrr+Y6M28qlMSUzORMStPW0na2ha7AMQYGrjLdH6HW0A6Z3OQi0wYxAUo8qpDpoKMDLr5Y33r0QAgZ4ciYMR+j4k6K2CfjGKOtzdxLhpubAZtoLfMkzz8m1WxEOjKpqVBcLBzN+nq49FKdi1sFDI+HhczI9AAAW7K36FXOqmBf4T4qTleQ7kyn19tLz1gPR7qOsG3bXlpbRXuZFDKxZcVXOD/60Y+04y1btnD55Zfz5S9/eaVfViJZERcW7aXtbAsUHKWq6iquuELviuLH+DgM+IWQWZth3rYyiNglk9ZBW3sAM06c7+wE9khHJsmahFWx4g/6wTFCW1uaqTdvNzUBduHIpNilkIGZjgzAli1CyJw4IYUMQEfvyJzHtmZv1aGS1cPrd77OpH+St/zbW+j19vL2Z98OwE07XoD/vlHmZOJAVP+qv/zyyzzwwAPR/JISybLQAv+FR0yXk2lrQxu9vD5rdWXR4k2RuwiCCtimaOjq1bucuHP2LIyMoGVkzOzIKIoSdmUco0xMiP51sxIpZJLtsrUMZjoygBzBPIuu/tEZ76fYU8Kut0mxW+24HK45wf+a9EcBGfiPB1EVMkVFRWRkZETzS0oky0IbI5pXYTohM2NimYmD/iD+yKRZxMV7gwmXYs5ZhmnCHTKRqDmZrEK5S6a5GbCL1jLpyAjmc2RATi5T6R2aKWQ2Z2027fCQ2cwO/nutYuiVFDKxJ6rN848//viiP1dRFO69995onl4i0diVt0scZNdz/KUJwKlrPfFkxg4ZE49eVslJKuHsRDftwx3APr3LiSuakMnoZhpzOzIQnlyWUzRKf5UQMnv26FyUTjQ1ARmytSySSEfGH/CzYVMQsElHJsTA2AhErE4xez4mktmOzIBPCJnOTuGMp5k7khdToi5kFEU557JLdbpZMBiUQkYSUwrdhXgcmQxPDVDdU4vfvwerSYYUCSEjHRmVYncxjROH6JkwryNDqszIQNiRySyQjkxTE3BJqLXMJlvLIHxX3evzUv7tcvw+CyjHaGy04POB3X6eL2BwRiZGZwgZs+djIpntyIxPj5NX1kdPczZ1deYcNBMvoipkzhXyDwQCdHd38/rrr3P06FE+/OEPc8EFF0Tz1BLJDBRFYXfBLg62/p4JTwVNTXvYuFHvquJDc0sQ1splmCrrskt4tReGAu2mC3cLIRNkKkkIGdM7MqGMjCdHtMiYdSmm3x+64fEW2VoWSeRd9eoz1QAkZwwxPpBJczPaSGYz4vPBeGBma9nWHClkVObbKVOw5xg9zddRWyuFTCyJqpC5+eabF/z4fffdx3e/+12efPJJ3v/+90fz1BLJHMrzhZAhv4LqakwjZBrah2CzuONc6jF32B9gc0EJ1EHA1U5fH+Tk6F1R/OjsBJIHCShTgBQyqiPjyjK3I9PZKS5MccjWskjm+zms2zxOzZ9F4N/MQubUKcAeFjIKSniojmROaxlAyrqjwHUyJxNj4p7S+tjHPkZeXh5f//rX431qicnYlR/KyeRVUFWlby3xpGlAuDEZjhx5gQKUZc7cJWMmOjqAlD5A5EOSbEn6FqQzakYmxSMuyMwqZJqaxFt7ipxaFolFscxps1uzQQT/zR74r60FksQNgBc+9AJ/+MgfWJexTt+iVhHzOTITmUcAsRRTEjt0GTexadMmjhw5osepJSZCC/znV1BZNX9uy2h4vRE7ZGRbGQAlaSEh42k3nZDp7AScQwBkJMuJkqojk+Q2tyOjChlrkmwtm83sO+sl64WQMXvgv7YWcIgbADvzdrK/dL++Ba0yIp83qsBrCf4BCEpHJsboImTa29uZnp7W49QSE7E1ZysWrJA8yPHGzvP/AwPQ0oI2enldlhQyELEU091Je4df32LiTKSQSXem61vMKkB1ZCwpQsh0d8PUlJ4V6UNzc+hAtpbNYfad9fwS6cgAVNf6wDYJhG8ISMJEPm9u2nwTTpuTganTkH2ClhYYG9OvNqMTVyEzPDzMV77yFerq6ti5c2c8Ty0xIU6bk40ZYjxks7eCiQmdC4oDzc3IiWWzKHAXoAQtYJ2mvrNH73Lihs8HPT1IIROBegEWsI6SlATBIHR16VyUDqiOTMAip5bNZrYjk1ssHRmA6pPhfIwUMnOJfN6UekrZXyIcq9QdrxAMyudPLIlq2P+aa64558e8Xi9DQ0MEg0GcTid///d/H81TSyTzsrdoF/WDNQRyK6irewe7d+tdUWxpakJzZGTQX2Cz2EhTChmmI7QUs1DvkuLC6dPiQt2SOkgAyHDK1jJ1atmob4SSEmhoEO1la9fqW1e8EUImyDSytWw2sx2Z9BwhZPr6oL8fsrL0qEpfpqfhVKsQMg6LA4fVoXNFq4/I36F8Vz5Xl13Ny80vk7z1FcZ+fy+1tebdWRVrourIdHZ2nvO/kZERCgoKeNe73sXPf/5z6chI4kJk4L+6Wt9a4oFwZOQyzNnkJImcTPuQeUIyah4oLVc6MirqneSRyRFKSgOAOXMyTU2AdYoA4mcgw/5hZjsyAauXklDMzqx31RsbYdoihIwrSbox8xEpgPNceVxddjUAI9m/ByUgA/8xJKqOzAmzN5FKVh0780KCOd8ck8uam4ENsrVsNsVpJTT0/tlUSzHVZZipmUMMIYUMhDMyL5x8Adtb3HD8OG1tJpnLHmJ0FM6cAZxe7THpyISZ7ciMTY2xZYvYOVRfD5ddplNhOhI5sUy2lc1PpADOd+WzIXMDqfZUxnyDkF1Hbe12HaszNrqE/SWSeKFNLss8xfEa78KfbAAaW8fBdQaQjkwkZZki8D/gF0sxzYAqZJLSpSOjEnkRNm3xwo6fmm4pphr0T88RbWVWxYrdYvKV9RHMdmTGfGNs3iyOzXqvtqYGbWKZFDLzY1Ws2nG+Kx+bxRa+kZpXKSeXxZCYC5nh4WGGh4cJmuXqQbKqyHflk27PAUuA4501epcTU4JBaBoQfTIpNpfMRESwpVD0hkyntHP2rM7FxAlVyNhcUsioqBkZDdu46VrLVCFTXBbeIaMoio4VrS78gZmTDVVHBszbWhY5ell1NSUzCQQD2rEnyQPM7AhpaIDJST0qMz5RbS1Tefnll3n22Wc5duwYE6FRUU6nk927d3PLLbdw7bXXxuK0EskcFEWhvGAXB9sO0GutYHDwQjIMen0/OAhjtvAOGXlxEmZ99sylmB6PvvXEAzUjoyQPgl+G/WGeu8npLbSZ7OJUnVhWuMZLNXJi2WxGpkZmvD/mG+NSkzsykUJGOjLzc/may7mm7Bp25O7Q/vaqHSG2okqmA3DqFFxwgZ5VGpOoOjLBYJCHHnqI++67j9dff53x8XHcbjdut5vx8XFef/117r//fh588EHp0Ejixp7CsL1r5MB/UxPa6OW1GbKtLJISjypk2jWnwuio36ffLh0ZldzU3JkPpLeYzpFRhUx6iZg7PednYnLOTs60bCMdmcZGMdbcTPj9IQHnkBmZhbBZbBy47QBff9vXtcdUR0YpqACQ7WUxIqpC5plnnuH5558nJyeHf/7nf+bw4cMcOnSIQ4cOcfjwYT73uc+Rk5PDr371K5555plonloiOSfhyWWVhg78NzejjV6WQf+ZhJdidtHabo5lvKqQmbJIIaNS6inl+Q88zxM3PCEeSG/h7FkYHta3rniiChlrtjhYn7lex2pWHyOTcx2ZoiJITRVjiNWfn1loahItUbbUUGvZ7PZMyTlRhYzP2QUpfVLIxIioCpnnnnuO5ORknn32WT74wQ/icoWVu8vl4gMf+ADPPvssTqeT5557LpqnlkjOSThwV0FVtXGdQLkM89zkpeahBG1gCXCio1vvcmLKlH+Kh1/+P7RbXgPAG5BCJpKbttzEBy/4oHjHfRps46YK/KsX4j6XOFiXvk7HalYf+wr3zXjf6/OiKGiBf7PlZNQuhpzCUGuZXToyi8Wd5GZdRuj3Swb+Y0ZUhUxHRweXXHIJJerQ9XkoKSnhkksuoaPDPPscJPqyNXsrVmyQPMSRk8Z93skdMufGarGSphQBhJZiGpfvHPkOX/7jl5j68BUAjPikkJlNhjMjHFr2tJmmvSwQCIf9z1qlIzMfX3vr1/jH/f/I313yd4BwZADTTi6rCc3IyciXrWXLIfJGqhQysSGqQiYzMxO7/fxjHO12OxlGTVxLVh1JtiTK0kSTc01/pWHH78rWsoXJdYSWYg4bV8wCNA40asc5BeNM+sWonIxk+ZqroigKa9PXindMlJM5fRomJsBqhe5J8TzR7hhLAMhKyeIr135Fc2bGpoSQMevkMtWRcWfJ1rLloK2AyKvk5EnRniiJLlEVMtdeey1/+ctfGF6g4XhoaIg33nhDTi6TxJW9xeKuiNdVadiwd1PLNKSJi/RST6nO1aw+Ct0iJ3Paa2xHJnJLe/5a4cZYFIu8kzoLMwoZ1Y0pKQ3SPBRqLZNCZl7UfTJjvjEaBxpZs1FcyJtNyKiOTLJHTi1bDqqQUQor8PnEwAhJdImqkPnkJz9JcXExt99+O3/+85/nfPyNN97gzjvvpKSkhE996lPRPLVEsiC7C8Lz3I0Y+A8EoKW/Cyx+7BY7Be4CvUtadazLEo7MgN/YjkzkON2s4gFA7DWwKHL/cSSRQsYsGRk1H1OyuY/RKXFhqv0cJDNItQshc6z7GBv/dSPfHnwPYK7WMp8vQrgli9eStKQ0/QpKQLTWsuwasPhke1kMiOoemXvuuQe73U5NTQ133nknHo+HwsJCALq7uxkaEncHd+3axT333DPj3yqKIieZSWJG5Ibdqip4+9v1rSfadHWBL7UFEG6MvGidy5aCEjgJk0ntjI9DskHXZ0Q6MqmFQrTJfMxctAv4jGbTODKqkMlYJw6K3EU4bU4dK1q9qI6MLyDmLR8bfA0s0/T32+jrg+xsPauLD6dOCTHjckHLWB0Am7I26VxVYlGWUYbL4RI3DrJOUlu7nZtv1rsqYxFVIXPo0CHtOBgMMjQ0pImXSI4fPz7nMbm8TxJLtBHMWfVU1EwAxvrjPWNimQz6z8uGvNAI5tAumQ0b9K0nViiEX0ttWS0QkEJmPsrSy8RBRhNth/WtJV6oQia5oAn8sq1sIVRHRmVieoKC7Q10V22hvt4cQkbNx2zZOcLhoRYAtuds16+gBMSiWNiRu4M/d/wZ8iuoq5M/v2gTVSHz8ssvR/PLSSRRo8BVgNuaxQj9HG6tBfboXVJUEUH/FgDWetbqWcqqpVRbitlBR4dxhYwa7geYdjfDsAz6z8fGrI3iIPMU7R1B/H4Fq1XfmmKNtgMlown6pJBZCNWRiSTnggpNyOzfr0NRcUbNx+TvFAcFrgKyUrJ0rCgx2ZW3SwiZvApqa2/RuxzDEVUhU1RUFM0vJ5FEDUVR2J69kzd6fk/TWAXT03uwRfXZry8zJpZJR2ZeStJCQsbdTUu7Dzj/hMVEZHI6LGTG7CLdLR2ZuazPCI0dTh7C7+inpyebUCe0YVGFTLdFdE/INqFzk2JPmfOYvaQC+IBpcjKqkElZWw0TsCNvh74FJSiRS7nr/gh+P4a/aRJPotpI//jjjy/KlXnllVd4/PHHo3lqieS8XLRG5GSmsyo5dUrnYqJMZGuZDO/OT05qDpagHZQgJzq79C4nZkxMT2jHA8EWANKTpJCZTbI9OSxuM08ZPiczMSGydCQN80bfbwH4681/rW9Rq5jZrWUAk+mVgHkC/2pr2VS6OLgg5wIdq0lctIxufgUTE9Daqm89RiPqQubAgQPn/bxXXnmFJ554IpqnlkjOS3nEXRH1BdooRLaWyR0y82NRLKQhcjKneow7pmpkPOzItI6KkUPZKSZo6F8GWntZlvGFTEuLeOss/zVTgSm2Zm+VeYcFmM+R6VEqAEwxeWpiAhoaxHGfJSRkcqWQWQ47ckNOlrsbUnpN8fyJJ7qMNvL7/VgscqqSJL5E3hWprDLWVszGpgB4xJWYbC07N9mhpZhtQ8Ydwdw/HHZkhifFTq+iNNn2Ox8bM8M5GaMLGbWtzLH7OQDev/39csjOAlgt1jkT3XonOyB5gOZmGB/XqbA4UV8vWqDS0+HUsBQyK8Gd5A63suZVUlenbz1GQxc10dDQQFqanEUuiS/bcrahYIGUft6s69a7nKgxOQmdw6fBNoVVsVKcVqx3SauWQpfxl2IOjkzOeUw+J+ZHEzImcGSamgDLNKO5LwHw3m3v1begBMCqhIMMmcmZALg3HyYQgJMn9aoqPqj5mHX7GukZ60FBYVvONn2LSmAib6RKRya6rDju/NBDD814/+jRo3MeU/H7/TQ3N1NdXc2111670lNLJEsi2Z5MSfJm2sbrqOypBIyR7G1rQ8vHFKUVYbMYaIpBlCnLKuG1QRjwG1fIDI9NzHmsyC0dmfmYMbnMuE8JICRkPK0ELFM4bU55UboIxnxj2vENG2/gx5U/xrPzNUaOv5W6Oti1S8fiYozafu0v/zYA162/bt5JbpLFsStvF8+feB7yKqWQiTIrvuJ5/vnntWNFUWhtbaX1PEmmzZs388ADD6z01BLJktldtJO2hjq6A5WMjb2NVAO8LkdOLJNB/4XZUlACDeC1deDzgd2Ag8tGvHMdGdlaNj+Rjkzr4SBg3FarpiYgS9gIGzI3yKW5S+TKNVfy48of4yt6FTB+TqamBrCNc8r9fQjCvRfeq3dJCU14KXcFdS9DMAiyszM6rFjI/PCHPwTEAszbb7+dyy+/nI997GPzfq7dbic3N1eOaZboxkWlO/lVw79rd0UuvFDvilaODPovns0F4aWY3d1QWqpvPbFgdGJixr5Xi2Ih35WvX0GrmHUZ67BgIZA0QmvvGSBP75JihhAyYlyjJuAki+aKtVcA0Jf0F7B7qaubOwzASNTUANt+gTc4QKmnlHdsfIfeJSU02gjmnFpGxnx0dtoplh2/UWHFQuaiiy7Sjm+++Wb27t074zGJZDUReVekutpAQibUWiaFzMKsSQ+N2/W009FhTCHjnZycIWTyUvNku+E5SLIlUegupmOkjYFgE15vHikGvD4NBkOvE5dLIbNc1mesp9BdSNdIFxS/QW3t1XqXFDO83pDwvVG4Tx+64ENYLXLxyUpYm74Wt8PNCCOQXU9t7QVSyESJqHrLX/7yl3nve2WAULJ62ZUXuiuSfYJjVXNbcBKRpiY0R0a2li2MtjfE1UNrx5S+xcSIcd/MjIxsK1uY9Zll4iCjmQ6DDrPr64PRUTRHRi7CXBouhwtFUbhy7ZXigTWvcuoU+Hy6lhUz6uqE+LWWHAXgwkID3PHTGYtiCS8UzZOB/2gim2QlpqI4rZgUJR2s0xxqMsZWs8iMjBy9vDDZKdlYAkkA1LR36lxNbJj0zxToMui/MGUZISGT3mzYyWXq6GVrrsjIaEMOJItCFX77S/YDYC19E58PGhv1rCp2VFcD1kkC2VUA7C3cq29BBkG7kSoD/1Elqv0GW7duXfTnKopCrfw/KYkziqKwOX0nxwZfo26gEkj8sTMNjUG4vgWQrWXnQ1EU0ihmiMbQUswyvUuKKl4v+JVZjowUMgtSlh52ZAwtZKxT+N3ihodsLVscL3zoBb74hy/y9LueBmB3/m4ALIXH8SOciy1b9KsvVtTUALk1BC0+MpwZ8u9KlNCEjBzBHFWiKmQKCgrmfTwYDNLb28v09DQAhYXGGHsrSUwuWiOEzNnkCnp7byUnR++Kls/AAAxN9oFdbGcr9Rgw9BFlsu0lDPkbaR00Xh9RZydgnenIyB0yC6O1YxrdkcloAiWAy+GSwx8WyY2bbuTGTTdq7+/I24GCgs/ZDalnqK0GCRB8AAAgAElEQVTN5eabdSwwRlRXAwVHAOHGyMWp0SEyo1v7opxcFi2iKmReeeWVc34sEAjwxhtv8MUvfpFNmzbx2GOPRfPUEsmi2Veyi28fB/Iqqa6Gq67Su6Ll09iI1lZW4CogyZakb0EJQIGrmIZhOD1uvMUhHR2ATWZklkKkI2PUXTKRo5c3Zm6UF6bLxOVwsSFzA6cGTokxunXX6V1STKipAXYKIbMnf4++xRgILSPjPs3g1BnOnMklz7iDEuNG3DIyFouFyy67jG9961scPHiQ7373u/E6tUQyg/BdkUqqqvStZaU0NCCD/kukLFME/gemjXfVKoSMzMgsBS0j42mjpW1a32JihBAy9YAM+q+U8vxycZB/3JDtQWfPhpYsF4igv8zHRA+Xw8X6jPXiHZmTiRpxD/uXlJRQXl7OL37xi3ifWiIBYHvOdggq4OrhUG2P3uWsiIYGwqOXZdB/UWzOF0Jm1NJBIKBzMVEm0pEpcBXgtDnDdwEl81LoLsSuOMDip7nfeO2GEHJuc8RV07acbfoWk+BECpkTJzDca0htLWDxQV4lAHsLpJCJJto+GSlkooYuU8ucTifd3d16nFoiIdWRSr5jAwBHOip1rmZliNayFkAG/RfL1tDw/qC7nTNndC4mynR0oGVkXr3jVVo/2Upuaq6+Ra1yLIqFIpf43ekYayYY1LmgKDM5GXpeSCETFdTAtlJQwfg4tLbqXFCUqapCPFdsk3iSPKzLWKd3SYYiPLmsgro6fWsxCnEXMp2dnRw+fJjs7Ox4n1oi0diZK15MmsYqE/qOmmgtE39JZWvZ4lBby/C0i3C8gWjvCIBN7MdJd6ZLEbNINmSL9jJfajN9fToXE2Wam8XAHXLEVZMUMitDdWSCWSfANm64i9HKSqAwlI8p2CPzVFFGa22Xk8uiRlTD/r/85S/P+bGxsTFaWlr49a9/zejoKLfccks0Ty2RLIlL1+/kpc6fM5VeSVsbrF2rd0XLo6EB2BFqLZOOzKLQlmKm9tLUNsHevU59C4oi7V1TEMrmysEPi2d9ZhkHmoF0EfhP5EmGs2lsBNI6IWkEm8XGhswNepeU0BS6C8lOyabP2we5NdTU7OOGG/SuKnpUVqLlY/YUyKB/tNEcmZxaan7lA+y61mMEoipkHnzwwQXVezDk2d90003cf//90Ty1RLIkdhfMDPwnopAZGYGeHsAjZsbK0cuLIzM5E0vAScAyQU17J+9jvd4lRY2O0+GJZU6bcQRarJm9S2aPga7fIvMxGzI34LA69C0owVEUhV15u3i5+WXIq6CmZp/eJUWNYDDUWvbu0OhlmY+JOmvT1+J2uBmZGuGM/wT9/TvIytK7qsQmqkLm3nvvPaeQsdvt5OTkcOGFF1JSUhLN00okS0YL3OXUcrzSxzvfmXh3RZqaAMcIJA8BUOKRv1eLQVEU0oIlDHEqtBTTGEJmagr6BsMTy+yWxHtO64U2ucyAu2Rk0D/6lOeXCyGTf1zsXDEInZ0wODwN+RWAdGRigaIo7MzbyZ/a/wT5FVRX7+CKK/SuKrGJqpBJRJelqamJAwcO8Mc//pGTJ08yOjpKeno6u3fv5vbbb2ffPuPcbZGEWeNZQxJuJm0jvH6yHrhA75KWjJhYJkYIe5I8pCWl6VtQApFlL2EocMpQSzG7utAmliVZk2Rv+xKY7cgYiRlCJlsKmWgwYwTzQfD7wWrVtaSoUFkJZJ8A+zguh4uNWRv1LsmQ7MrbJYRMqCNECpmVEVUhE8mxY8c4fPgwPT1ivG1eXh779u1j9+7dsTrlsrjjjjvo7e0lJSWF8vJyPB4PDQ0NHDhwgAMHDvDQQw9x22236V2mJMooisJ6105qR/9EdW8liStkZFvZcihILaZxBLrHjLNLJnKHjGwrWxqaI+PuprluHEjWtZ5o0tgIXCoS6VtztupbjEHQcg75FYxPBGhutrDBANEjkY8RbWW783djUXQZbGt4wrvsKsTPXLIioi5kmpubeeCBB6gO+a1qLka9O3jBBRfw1a9+lbWrJJSwYcMGPv3pT3P99dfjcIR7h5977jk++9nP8uijj7J//37WrzdG+4kkzJ7indSe+BOd05VMTd2CI8Fax0WIV1yISyGzNNZmlPDHEeg30FLMyB0yMui/NLKSs3BaXEwERmnsawW26F1SVPD7Qy2oN54E5DLMaLElewsOq4OppBFIb6G6ep2BhExoEabMx8SMyF0yVX/StxYjEFW5febMGW699VaqqqrIycnh1ltv5aGHHuLhhx/mtttuIzc3l6qqKm699VbOrJIFDv/2b//GO9/5zhkiBuD9738/+/fvx+/389vf/lan6iSxZP968WISzK2gvl7nYpZBpCOjTeKSLIrIpZhG2RsSuUNGOjJLQ1EUilOFK9M20qxzNdGjsxOmlLPgEn9vN2bKVqFoYLfauSA35OLnH6emRt96okVVFdro5b2FUsjEih25O1BQwH2ayqaehF4BsRqIqpB56qmn6Ovr44477uDAgQM8/PDD3H777dx222089NBDHDhwgI985CP09vby7W9/O5qnjglbtmwhGAyuGtEliS678sKTyxIxsClby5bP9lKxFNOf2s7QkM7FRIkZjoxVOjJLZUOWEDKDwWYmJs7zyQlCYyOQ2QBATkoOHqdH34IMRHleKCeTV5GQfz9mMzkJdfV+yD8OyKB/LEl1pLIhdFPB666gpUXfehKdqAqZV199lbKyMh588ME5DgeIyWUPPPAAZWVl/P73v4/mqWNCe3s7iqKQY6SlAhIN7Y5aWhd/qUqsLXgTE6ELV49sLVsO67PDSzE7DJL3lxmZlbEpz3iTy4SQOQUgg9tRRmsPMsjkshMnwO85CY4xUuwpbM7arHdJhmZ3QXhgRFWVvrWslDNn4MYbQa/mpagKmd7eXrZtW3gqiqIobNu2jd7e3mieOuq0tbVx8OBBAK6++mp9i5HEBHeSm2zLOgD+0pJYryRiWzco6aHWMjl6eUlorXgp/TS2jetbTJSQGZmVsS5iclmzQbrLGhuBrJCQkW1lUSVyctmJE2L8eSITmY8pzy/HajHAGLZVjOboGUDIPPssvPgi/Pu/63P+qAoZl8vF6dOnz/t5p0+fxuVyRfPUUcXv9/Pggw/i8/m44YYbzivOJInL1kxxV61+uELnSpZGQwOgBGTYf5mkO9Ox+lMBqG4zhiUjMzIrI3KXjFFaPSIdmQ2ZBkijryK0yWXpbUzbBjl1St96VkplJeF8jAz6x5ywED6W8EJGrV+vgRdRnVpWXl7OwYMHOXjwIFdeeeW8n/Pqq69y9OhRrrrqqqic895776WpqWlJ/+bRRx9lx44d5/z4F77wBY4ePcqaNWt45JFHVlqiZBVzSdlO/tD3PIP2SkZGwO3Wu6LF0dAApPQStE6hoFDkLtK7pIRCURTcwWKGqKe+ux1I7LvV09PQ3Q1kyozMcgnvkmkyliNTLjIy0pGJLh6nh7L0MpqHmkOLDa9k+3a9q1o+VVVojozMx8QeTchk13P84BiQqms9K0Ftrdy+HX7+8/ifP6pC5q677uK1117jvvvu4+1vfzs33ngjxcUiVNvZ2cmLL77Iiy++iMVi4a677orKOTs7O2lZwu0zRVGYWCDJ+dRTT/Gzn/2MnJwcvve975GWtvQlg36/n8A8Yyh8Pt+Sv5YktlxSthPeRAv8X3qp3hUtjsigf74rH7tVbnFfKln2EoaC9bQMJL4jc/o0BAJgcUwSQDoyy0FzZJKHqD8xBKTrWs9KCQZDQuYamZGJFeX55SEhc5zq6iv5wAf0rmj5VFQGYI8cvRwvCtwF5CTn0TveQ8PZaiYmLsaZgC/bgQDa1L5IIe/z+c55zWuxWLBGcYNsVIXM7t27+dKXvsQjjzzCCy+8wG9+85sZHw8GgzidTj73uc9RXl4elXP+8pe/jMrXAfjpT3/KN77xDTweD9/73vcoKVle7uDJJ5/k8ccfj1pdktihtQfk1FBRNc2ll8ZsR2xUOXkSyBC3jbULMMmSKHKV0DgCXaOJv0umPfQtpGdPMIDMyCwHl8OFx5bL8PQZTvU1A6trefNSGRiA4YlhSBV5VOnIRJ9debt4/sTzCR/47+2F01MNkDSC0+aUi1PjxJ7Ccv6n8X8I5B6nru5iVtm++EXR1AReLyQlQVnEpchFF110zn9z3333cf/990ethqhftb3rXe/i4osv5rnnnuPIkSPa6OLc3Fz27dvHe9/7XgoKCqJ92hXz4osv8oUvfIHk5GS+/e1vs3nz8id23HPPPdx9991zHvd6vQv+z5XEn7KMMuzBVHz2Mf5Yd4q7SYwX8Pp6oDQkZNKlkFkOazOLeW0E+nyJL2TUKVuerEkGkI7Mcil1r6Nq8AztI00kupCJDPrnpebhTkqQvtkEQmsPyqug+mV9a1kJx46htZXtytuFzZIYN/QSnfJ8IWTIP05lJQkpZNR8zLZtYIt42hw6dIiUlJR5/43FEtV4fvSFDEB+fj6f+MQnYvGlY8Krr77KP/7jP2K323nyySdX7BZZrdZ5bTO7Xbb/rDYsioXSpB00Tr3Bsc5KSAAhMzYWugO/SwqZlbCloARaw0sxFUXvipaP6sikZciMzErYnLuOqsE3OGttYmwMUhO3bT0kZMSm301Zm/QtxqBoQia3hobmKcbHHSQn61vTchBCRgT9ZT4mfmjPn4LEDfyrdc+Ondvt9rhd80ZXFiUgR44c4W//9m8BeOyxx7g0UUISkqixM18sxmz2VibElnd1Oo4tR7aWrYQLSkXraMDVziqfBn9eVEfGlS6nlq2EzbnhwH+iTy4TQuYkgNwJEiNKPaWkO9PB6oPsOmpr9a5oeRw/jubIyHxM/Ag7epVUVvn1LWaZVFUByf0c3fgeDjQd0KUG0/uHd999N5OTkxQXF3PgwAEOHJj7P2LPnj28733v06E6STz4q827eL4NxtMq6emB/Hy9K1qYk+LaBFt2M9NIR2a5lGWJQSR42mlvh9xcfetZCaojk5I2AePSkVku6zLEXil1clkiT6FqaEATMtKRiQ2KolCeX87BloOQf5yKil3sTUAdcPRYEG4KCZnCBPwGEpSNmRtxWlOYwMvR1lPAFr1LWjLV1UD5M1T7/5OfVjt0qcH0QmZ0dBSAjo4OOhZY8S2FjHHZVywcGfIrqKpa/UKmvh5QAkyltAIRF1+SJaEtxUwe5GTLGHv3Jm4fkSpkkt2TMC4dmeUSFjKJv0umsRHYLlrLNmdLRyZW7MrbFRIyFcLZSDBGR+HkmWZIHsJhcbAtR+7NixdWi5WduTs51P0G/fbjnD69ZdVff0QyMRHqELnpTQC2ZG/hKEfjXofphUxdXZ3eJUh0ZkduqLnT085fKge57roMfQs6D/X1gLuLgDKFzWKjOK1Y75ISEo/Tg83vZto6ElqKmbgXe2prmSMllJGRU8uWhSZk0ltobPYDibvdvKExCH8lHZlYE15seDwhhUxlJVo+Zmf+ThxWfe6qm5U9heUc6n4j9Pz5IG97m94VLZ66OvD7wVLyJgFEW+JP+Enc6zB9RkYi8Tg9eIJrAHi9YfUn7urrgXSRjyn1lGK1JO7Flt6kIUTgydOJO7lsfBwt4+NIlhmZlVDkLsKKHaw+TnR26l3OsvF64fRoFzjGsCpW6drGkEghc+x4kHlWyK1qjh8HCkNB/3wZ9I834efPsYQTwiIfM0AgvRHQb1CEFDISCbAxTeyTqemr1LmShQkGQ0ImowmQ+ZiVkuMQ7WWtg4m7FFPtiE1JgaBVTi1bCVaLldwkcVOjob9J52qWT1MTWj6mLKNM3mWPIdtytmG32CF5kFFLW8K1JEaOXpb5mPgTObns2PEEmDYUQVUVmpu3IXODGHyhA1LISCTAhSUiJ9Ppr8C/ioeH9PTAyAiQKUcvR4PiUE7mtDdxHRk1H1NSApN+6cislDKPcC86vYkrZBobgexQPkZOLIspDmtEriT/OBUV+tazVI4eC8rRyzqyI28HFiyQ2suRE6f1LmdJVFcDRSIfs69wn251SCEjkQBXbBFCxp9VSXOzzsUsQL24NsFVJEIRa9PX6leMAVAnl/VPJ76QKS2FPm8fINolJctja4EQMuNJzQwN6VzMMokcvSzzMbFHEwAFRxOqPcjn4/9n787Doyrv/o+/ZyYz2feQfQVCQgj7IouIgqBoa9G61eWxrV0son36tM9l21/Lo7/+Wq32ad1ta61aLdpWAaWl7sgmS8KaQBIISVjCFtYkhOzz++MkQxCQhMzMmUk+r+vqxSGZOed7lWMyn3Pf3/umaPceCDlGgCXgTL+oeE2IPYRB0cbDhl0Nm+lYf8ovFBUByUaQGZ883rQ6FGREgFFJHSuXxRezeavvDsl0BpmgWOPJTVJ4konV+L+hKcaITEPAPlpaTC7mEnU2+qelwb5aY56Za0U26bGcAWeWYPa3aUKdum6GqREZz/PXIFNSAi0DjNGY/Ph8LRJiknEpnfvJbPabjTGPH4fqalwjMgoyIiYbHDMYW3swOBpYUbTL7HIuqDPIWMMPAZAQmmBiNf4vryPIELmH/fvNreVSdY7IpKa1nwkykQoyl8q1wWzHXjL+SCMy3uWvQcbYCNMIMuqPMU/XPhl/uX+2bAHCDkBENVaLldFJo02rRUFGBKPJNzkgH4CCPb7b8N8ZZJodHUEmTEGmN7KijcZuInezZ49/NVp26hyRiU45QnNbMxYspISnmFuUH+u6Kaa/jsiUVzZDtJHCtIeM541IGIEFC0TsZ8/RQxw/bnZF3bNpEz7xRL2/O2vlu03m1tJdmzYByYUADI0bSpgjzLRaFGREOuTFGtPLdtb6epBxUu88DGhEprdcIxeB9ZRW+WdDROeIjD3OOEgMS8Rus5tYkX9zBZmww5RV+tGE9Q4tLVB1sgKsbYQGhJEUpumnnhbmCDsTGJM2+U3D/8ZNTteHUTObtfs7V5CJKWdDUZ25xXTTWSE4xdwQrCAj0mHKYGMJ5qP2rTQ2mlzMeTQ3Y0x1CT5Oq9No6IgPjTe3KD8XYg8hsG0AAMX7dptcTc85nWdGZIgwgow2SO2dqKAoQizGprglB6pMreVSVFVBW1RHf0zcECwWi7kF9RP+Nr2svR02Vu6C4OPYrQ6GJ6jR3yzxofHEByeDxUlxzVZaW82u6OKMERnfGM1TkBHpMDW7s+F/CyUl5tZyPhUVxi66wXHGtLKooCg1Z7pBtCUdgJ2H91zklb7n5Elcq9w0Bao/xl1SQoxRmYrj/rcE886dnOmPiVN/jLe4NpP0kyCzYwfURxijMaMSR2mvIZONSzV6TJpjN7qmkPuq06dhe4nTZ6YlKsiIdBiR2PFEKrqK9VtPmlvM51Qer+TldW8BTlJy1OjvTknBRp/M3lr/G5HpnFYWEwOHTht/0YplvTco1mj4P9BU4Xc7te/YgfaQMYG/jcgUFuJ6oj7B5KlBAuM7p/YlF/r8/VNcDO3huyHkKHarnREJI0ytR0FGpENMcAyhbca0nOUlxSZXc7ZvL/k2j1fdAumriR+oRn93yogygsyhJv8bkem6h8zeWk0tc5dhycaITFt4hd+tZnfWiIxWLPMaV59DdCXbKo7T1GRuPRdTUIDPPFGXLj1KfhBkuk4rG5EwwvSZIQoyIl0MDDH6ZDbt962G//11HZ+moncRlWrsIaMRGffITjCmlp3E/0ZktIeMZ2THdq5ctstYytiP7NiB9pAxQXRwNFlRxkhea9wmn98PZH1BGyRtBNTo7wvGJnUsfz2ghHWbfHuREV9b7U5BRqSLsanGEGllwxacPrQa76mWU8ZB6GFXj4yCjHuMSDdGZJqDd3PqlMnF9FDniExa2pkRGfXI9N7gmMHGQYz/BZmy3ScgzFjVMDs22+Rq+peu08sKC82t5Yu0tsLGvaXgOEVIQCi5cblml9TvJYUnER+UAhYnG6o3+/SU1rMa/X1gWqKCjEgXV+UZQaYpaqvrQ6IvONV8JshYwjW1zJ2GJhkjMkTuObMCmJ/orDc1rZ3q2mrjWFPLes0VZKIr2LmrzdxieqCxEfY27AQgISSJiMAIkyvqX/wlyGzbBs1xxgfRccljsVltJlckABPTjZGxhqhCY2TVB7W1wZat7ZBsbKTqC6N5CjIiXYxL61y5rIjCDb7zSKS+I8jYIg7TYNGIjDt19sgQfpDyKh9cd/sLdIZta8J2WtpbCLAGaN8QN0iNSMWGA2wtFO/xoScaF7FrF65pZbkD1B/jbf4SZAoK8Kkn6mIYl9wxvSy50Gfvn7IyaAzdAYF1BAcEkzcgz+ySFGREuhoSOwSbMxAC61m2qcrscgBoa2+jqc34gB2acIjDDRqRcafY4FhsbSEAbKncZ3I1PdM5IrPJ+QoAXxryJW2G6QY2q43EQKNPZseRcpOr6T6jP8Z4lKv+GO9zBZnYHRTtqOX0aXPruZDCQnyqx0EMXRv+CwrMreVCNm/GFYLHJI0hwBpgbkEoyIicJcAaQHLAMAA+q/CN7ZkbWhpcxwGRhzlUrxEZd7JYLIQ7jell26v9p+G/ra1jRMbWzEc1fwHgm6O+aW5RfcigaGN62d5T/hNkjBXLjBEZrVjmffGh8WREZoDFSXtiIVt841fIOdYVNkOCUZwvTA0Sw9jOEZm4MtZsrDW3mAvwtUZ/UJAROcfwjjXRd5zwjZXLXI3+QGvgIQ6fMhp5NSLjPvEOY3pZxTH/aZI5cMBo2rXmLOVoYw2JYYnMzp5tdll9xvAUI8icDi7n2DGTi+mms0Zk4jQiY4YJKROMg5R1Pjk9qLERig8XQUAzUYExDIweaHZJ0iE+NJ7kEOOh2pZDm2hpMbmg8zAa/Y0b21dCsIKMyOdMyzGWYK4P3cqBAyYXQ5dGf6CW/TS1NWGz2EgOTzaxqr4lLcL45bH/lP+MyOzuKDVs6GcA3JR7k08M8/cVQ+M7Vy4r95uVy3bsbNceMia7LOUy4yB1HRs2mFvL+WzdCq1JawCYmDoBi8VickXS1WUdDf/NcYVs22ZyMZ/jdMLGLS2QuAnwnf4qBRmRz3E1/CdsZeNGc2sBONnQcM7XBscMxmFzmFBN3zR4gDEic7TNf4JMVZXxpyOhAtATeHc7swSz/wSZ0ur94GjAZglw7Wki3nVZakeQSVlHQaEPreHfoaAASDWCzKS0SeYWI+fo2vDva30yFRVwPGAb2BuJDIw88zPSZAoyIp8zomNqGbHlfLbB/HmqpRXnbm6idf/dKz/VCDIN9j00N5tcTDd1jsi0RxpBRlNE3OvMEsy72FnuOysYXkhdHdS0G/0xWZEDteiDScYkjcFmsUH4Qbbv2+dze1MVFABpxiju5LTJ5hYj5/Dlhv9163A1+o9LHofV4hsRwjeqEPEhcSFxRFuMqUbLd5o/JFOy69zfhEPjhppQSd+Vn9a5l8xu9vnJwmWdQeaUQ0HGEzKiMrASAPZGinfvN7uci9q5E4jrWHo5XtPKzBJiD3E9DHMmrzNWefIhyzcchOgqLFjO9POIzxib1DEiE1vOms3HzS3mc9avx+ca/UFBRuS8hsca/5EWHzP/kciOyvMEmQEKMu6UFd2xl0zEXiqrfP/pO3RMLQs6TpPlJACZUZlmltPnBFgDGGDPBKD0sO+vXGYEmVJADzrMdqbhf71PNfwfOgRVrca0sry44dow1QfFhsSSGWGMBm8/uc6nlvA2RmR8q9EfFGREzmt6rhFkToYWUFNjbi2V1Zpa5mnJ4cngtEJAM0UVh8wup1t27waijdGYpLAkQuwh5hbUB2VFGh8o9tT7fpDZsQNXkNHPB3O5Gv59bOWyNWtwTSu7PEP9Mb5qaqbxb9Oe8pnPjOg1N8PGogaILwJ8p9EfFGREzmvqwI7/SJMLjOUGTbTngIKMp9ltdkLbUgAo9oO9ZJzOs4OMppV5xrAkI8icsJb71JPR8+k6IqOfD+ZyNfwnF7JmXau5xXSxejWQ1tHon6og46tcvUupa3ymT2brVmiOKwBbKynhKaRFpJldkouCjMh5uOapRlfx6XrzhmQaG+Hw8bODTEp4iqYEeEBcgDG9bFeN7+8lc/iwcW8QYwSZrGitUOUJ+clnVi6rqDC3lovZXl4PkXsBBRmz5cblEu6IAEcDu2q3c8hHBnlXftbsmhqkRn/f5QqZqev4bG2bucV0WL+esxaJ8KVluxVkRM4jMiiSeKuxnO3HpeY9EikpAafdCDIWjB8c6o/xjOQwo+F/b63vj8h0NvqHpHSMyERpRMYTsv1kCWanE0prjP1jYgLjiQmOMbmi/s1qsTKhc+pNyjpjJMRkjY2wYf8mCGgiOjDOZ5bOlXPlx+cTbA2DwDpWlfnGZjLr1uGzq90pyIhcwOiEjob/4wU4TdoOYONGwG7sI3Nz3s3ckncLP7n8J+YU08cNijVGZGpafD/IdO4hY4/X1DJP6rqXTHm57+0J0ungQTgVXAJAXrxGY3yBq+E/dR2rVplbC0BhIbQmGtPKpqRP9Kkn6nI2m9XGhI4+q2rLGp8Y0Vu33qkgI+JvZuYZQaYhqsD1wdHbNmwAOkZkBkUP4u+3/J3pWdPNKaaPG5aSCUCdrYp2H1+4rHNEpjXcGCbQ1DLPyIzKxOK0guMUW3b5wKeJCygpwdUfkzdAQcYXTEydaBykrvWJIGP0x/jmB1E5V2fDP2lrWLvW3FpOnICymh0QcowgWxCjEkeZW9DnKMiIXMDkjDMN/2vWmPM0duNGwGEEmVBHqCk19BejMoww4Iys5MABk4u5iN27AVsTDXYj0QyJ1b4hnhAYEEis3ZhyWLzfd1cuKy1Fjf4+xhVk4rexcfsJ0zfGXL0aSO1o9E9To7+vc/0bpa4xVpszUWEhkG7MjxyfMh6HzWFuQZ+jICNyAaMSR2F1BkDYYT4q2Ov167e2wpYtuEZkQu0KMp40KDbTOIiqoqrKd6cRQcfUsugKnJZ2whxhJF9X+IcAACAASURBVIQmmF1Sn9W5BHPFCd8NMl1HZBRkfEN8aDzZMdkAtCWtNZqlTeJ0wsqi3RC5D5vF5lObGcr5uYJw3A5WFB4xtRZf7o8BBRmRCwq2B5MWmA/AqgrvN/yXlBgNmgEhRpDRPiGelR6ZDk4LOBooqjB586CL2L0biN0JGKMxmu/uOZ0rl52w7qS+3uRiLqCktA1ijWZ/BRnf4frQl/aZqdPLduyAE5GfAjAuabxG9/1ATHAMA8ON/5YLD66lqcm8WhRkRPzYxHTjydWuxgKv/yDZsMH4MyxGU8u8ITAgkJC2ZAC27qk0uZoLc+0h0/HBVdPKPGtEsrF6IXFlxqaTPmhbdRUENOOwBpERlWF2OdLhTJBZbWqQWbUKyPwUgKuyrjSvEOmRKwYZ08taEtaYNqLX3g4rCo7BAGMxEQUZET8zPdcIMu2JBV7fYbczyARHaGqZt8TajD6ZHTVV5hbyBY4fh7o6zgSZGAUZT8qJ7QwypZSVmVvL+dTVwcFWY1rZkJgcrBb9WvcVZzY2XMdna1tpM2lLkE8+ATKXAzAtc5o5RUiPTekyovfpp+bUUFwMJ8ON1QayY4YQFxJnTiFfQD/xRL7AJNfKM+tYvca7OzRv3Gj8aQ/RiIy3pIYaQWZvne+OyHSuWGZP1IiMN+TEdQSZmHJKd/jG5nRdlZXhelo6LEHTynxJ3oA8IgMjwXGK+rDNFBV5vwanEz5cvxuiK7FiY0raFO8XIZfEFYRT1vPJ8hZTali+HNe0sinpvjcaAwoyIl9oWPwwgjB+Eb2/eavXrtvWhmsEyOnQiIy3DIzJBOBwc5WpdXyRziCjqWXekRGZQQCBENDEpgrf22NIjf6+y2qxckXGFcZfsj4xZXpZSQnUhBqjMWOTxhEeGO79IuSS5MblEumIAUcDayo3mtIns3w5kG7cuJNSfXO1OwUZkS9gtVgZHmU8hSg87L3fQqWl0NAAoaHQirEhpkZkPC8/1RiRqbVVmrYJ6sVUVQGOOlqCjDWis2OzTa2nr7NZbSQHGf8fbz9canI159LSy75tRtYM4yDrE+NDoZd9/DFd+mM0rcyfWC1WpmVeDkBT4goKvLzmkNMJy1efdi3bPS3DN+8fBRmRi7gmzxiKPxa62ms77HZOKxs9Gk61aETGW8YMzASgPaKSw4fNreVCdu8GYoylgAeEDCAqKMrcgvqB3I7pZXtPl/lcwFWQ8W2uDYzTV7JsRbPXN9v9+GMgw0hQV2Ze6d2LS6+5RvQyVnq9T6akBI4EfwYBzSSHJfvs6L+CjMhFXJ1tPBEhfRWrVnnnU0xno//oMU5ONatHxluGDDBGZIjcTfkuL3/i6KaqKjStzMvGpBsBoTm8jP37TS7mc4p2HYGQo1iw6H7wQfnx+QwIGQCOBo4GraO42HvXbm2FTzbsgZgKbBYbU9LVH+NvpmZMNQ7SV7HsU+/+Tlq+HMhaBsD0gdN9dpl/BRmRixifMt7YGDNiP++urPLKNTtHZEaOaabNaTQYax8Zz0uNSMXitEFAMxt2HDC7nPOqqEBBxsvy4jsa/mPLfGrlspYWqKjfBkBKWIZ+Rvggi8VyZlRm4MfGCmJesnEj1MUYozFjksYSERjhvYuLW4xOHE2wLQSCj7N6xzav9smsWAFkGTfsVZlXee/CPaQgI3IRIfYQBoWMBWBZ+WqPX6+tDTZtMo5zh59yfV1TyzwvwBpAaFsaAFt3V5lbzHk4nQoyZnBN2Yor9am9ZCoqoC3GCDKjkoabXI1ciCvIZH3i1SBj9Md0Tivzzf4G+WJ2m52pGR19Mikfe61PxumEZZ/VQbJxQdc97IMUZES6YcYQY0h+r2UVR4969lrbtkF9PYSHQ0qmEWTsVjt2m92zFxYA4u3G9LKyw763BHNNjXFvKMh4l2sJ5vCDbC2rNbeYLkpLgXhjrlJ+/DBzi5ELcjX8p67l09WnaPHSSroffeyErI8B323UloubNWiWcTDofa8tGFFeDoccq8DWSmZkFplRmd658CVQkBHphlk5Z/pkVqzw7LXWGAuEMGECNLarP8bb0sIzAd/cS6aiAsCJZYARZLJjtGKZN0QERhBhTQRg8z7fmVtWUkKXIJNvbjFyQQOjB5IekQ62FuqiV7Funeev2dgIq0pLILoKhzVQjf5+zBVkMpfz/seNXrnmWf0xWb47rQwUZES6xdUkGb+N95Yf9+i1PjP2nmLSJKg5VQOglam8KCfBGJGpaa0yt5DzqKgAQo7iDDwBwOCYweYW1I8MjDCml5Wf8J0lmIu3ORVk/IDFYmHGwI5RmYEf8957nr/mypXQnPEvAK7KulIPw/xYfnw+8cFJYD/N6r2rqavz/DU/+ghXf4wvTysDBRmRbokPjSfJYTz9/rDkM49eq3NEZtIk2HzQ2BVzeLzmv3vLyIxMABoclZw+bW4tn7drF65pZemR6QTbg80tqB8ZnmRML6txlpmyMd35bNpxEIKPY8V2Zvqb+KSufTLeCDKLFgHZSwG4Pvt6z19QPMZisTB7iDEq0571PsuWefZ6bW0YD2wTjWbdqzQiI9I3TMsyppdVtq/k2DHPXOPIEdi50zieOBE2HjSWLxuTNMYzF5RzDO/YFJOoKmOpYx+iRn/zjErrCAoxZZSXm1sLGCuW7ThhjMZkRgwmKCDI5Irki7j6ZJI2sqGkxqN7krW3w8KlJ107sl8/REHG353pk/nA40G4oABORq4Aazs5sTkkhyd79oK9pCAj0k3X5HZsTJX1icf6ZNauNf7MzYWYGNh4wAgyoxNHe+aCco6s6EzjIHIPO3e1mlrL550VZGIUZLxpaJeVy0pKzK0FjAcerTFGkBmVrGllvi4pPImRCSPB4oTB73v0w+j69XAo9AOwtZITm8vA6IGeu5h4xcyBM42DxC0s+fSgRzfmfe89YOBHgO9PKwMFGZFuc/0gSS7k38s9s3RZ54okU6ZAY2sj22u2AxqR8abk8GSs7Q6wtrGxfJ/Z5Zyl69Qyjch4l2vqVuxOira1mVsMUFSEqz9meIKCjD9wTfHKXsq773ruOosWAUP+1XHN6zx3IfGaAaEDGBVvfA7YZ//Qow9T/v2e03X/XDv4Ws9dyE0UZM7jueeeIzc3l9zcXJYsWWJ2OeIjUiJSSHXkgcXJBzs8sxnAp58af151FRQfLqa1vZXY4FhSI1I9cj05l9ViJZIMAIqrq8wtpovGRqiuxhVksmO1Ypk3ZURmEEAgBDRRULbH7HI6goyxh4wa/f3DdZ2hYvB7vPdBG40eWIDK6YSFi9ph8L8BTSvrS67NPjO9zFNB+NAhWF9RCtGVOKyOM1MifZiCzOdUVFTwhz/8AavVisViMbsc8THXDjFGZapsH1JT495znzxp7MQMMG0arNy9EoDRSaN1L3pZQmAmALuOVphbSBdVVYClHWKMBg2NyHiXzWojNdgIj0UHzZ9btrWoHQYYQWbYAO0h4w8uS72M6KBoCD5OQ/RaY8NKNyspgfJTGyDsMOGOcC5Pv9z9FxFTXDP4GuNg0AcsfqfdI9dYsgQYbCwScaWfrHanIPM58+fPJyIigunTfX9eoHjfjSO6Nty5d5LqypVGk+bgwfDg6pv4rw/+C4CxSWPdeh25uM5ljfc17DK5kjN27QLCq8F+mgBrgE9vUNZXjUgyAkN1S7FHnqb3xKaKPRBYj93i0DLcfiLAGnBmqk7uOyxe7P5rdJ1WNmvQLBw2h/svIqaYnDaZMHs4hB1m3b4CDh50/zXeeYcu0xL9YzRPQaaLv//972zYsIEf//jHhIeHm12O+KBpGdOwOu0QtZu/f+zepYs6p5WNmbGLRaWLsGDhxtwbmTdhnluvIxc3PMX4YHicco82VfZE10b/QdGDCLAGmFtQPzQ+w1gG3TmgmDIT98Wsq4N9TcZozJCYXOw2u3nFSI/cmHujcZC7iLcXOmludu/5Fy0Chr4NwJeGfMm9JxdTOWwOrh/SMT0xZ7Hbp5fV18MHK49BujEbREHGzxw5coTf/OY3TJ48mS99Sf/xy/mFOkIZGW1sjvnJ7g9oc2PPb+c0g+AR7wMwNWMqC29bqP4YE4wbaEwhaoss9+gyqT2hpZfNN7yzFyW+mOJi8+ooLsbV6D9SK5b5ldnZswm0BUJsOccDit26ell5OWzYXQIJxditdr6S8xX3nVx8wpzcOcbB0EX87W/uPfd770Fz5j/B1kp+fD6DYga59wIeoiDT4Re/+AVNTU08/PDDZpciPu6mUUafTEPihxQUuOec+/fDZmPvSw6GGkHmmkHXuOfk0mND4zum6sTuZNcu3xiS0Ypl5nM11Q/YztZi85bm7rpimfpj/EuYI+zMniBDF/H66+4792uvAcP+AcDMQTOJDo5238nFJ1yXfR0OqwPiyvikqIQDB9x37jffBHIXAXBT7k3uO7GHKcgAy5Yt4/333+e+++4jLS3N7HLEx7lWDslcxuIlLW45Z+dTubETmlm931gRzR+WPeyrsqKzwGmBwDo273Tzqg6XqOuITHaMViwzQ1Z0Fg5CIKCJ9TvN65/qGmS0Ypn/uWlox4fEvH/w7rtw4kTvz9neDn/5i3FOgFvybun9ScXnRARGMGNgx0piQ9/mrbfcc96TJ2HJ+6dgsPEg9cahN7rnxF7Q74NMQ0MDjzzyCAMHDuRb3/qW2eWIHxidOJowawwE1bJg+Xq39FD821gpk/zZn1HfXM+AkAGMShzV+xPLJQkKCCK0zXiosXG3+du4O52dQWYnoBEZs1gtVgaGGyMg244UmVbHluIWGGCsnKYg43/m5M4xmvATimmK2sqCBb0/56pVUNW4GRKKcdgcmlbWh9067FbjIP8NXv+re2YMLF4Mzen/BvtpMqMyjc1b/YTfd4vef//9VFT0bInUxx9/nOHDjabN3/72txw6dIhXX30Vu10Nk3JxNquNWYNnsnDH39gb9E9KSqaQl3fp52tpgQ8+MI6Pp70O+4zhY6ul3z9nMFWCfTAVzj2UHS4HJptay8GDcLqpBaKNn3UKMuYZmzqc0pICaizF1NXdjLfXhXE6YUt1CQQ0ERYQSVZUlncLkF6LCori+uzrWVS6CIb/lT/+cQTf+x70ZpX9554DRr0MwA05N2haWR82J3cO313yXZrjt7O+qoht20YwrJczTF9/HRj+VwBuzbvVr7Z88PtPStXV1VRVVXX7f7t376axY93MrVu3smDBAr7yla8wYcIEt9XU1tZGS0vLef8nfcMtwzsb7hby9tu9eyKyYgXU1kJMUi0fHXoDgG+N0eig2bIijT6Z3fXmj8hUVADRlWBtI8QeQnJ4stkl9Vtj0zob/ovYvt371z9wAGpDNgEwOmmUX33gkDPuHH6ncTD8DbZsbaew8NLPtXcvvLWoGUYYH0S/MeobbqhQfFVUUBTXda5elv8GL77Yu/NVVMBHq49BtrF/zF0j7uplhVzwM3BLSwtt7lwliT4wIrO4Fwuxr1ixgvb2dsrKyrj77rvP+l5lZSUAL7zwAn//+9+ZOnUq3/nOd7p13ueff55nn332kusS33dd9nUE4KA1bgevvVfCz39+6UMynSuP5N22gFUtDQyNG8qUtCluqlQuVX7SYD4+ATWtPhJkujT668OrefK7rFxWVASXXebd6xcVAUlGkBmbPNq7Fxe3uX7I9UQGRnIyci8M/JBnn72GV1+9tHM9/zy0D14CIUdJCks6s5iA9Flfy/8ai0sXw4i/8uqr/4/HHrMRFHRp53rpJSDvLQhoZkTCCIYnDO91fV80ODBv3jweeOCBXl+jk98Hmd6yWCyUlpZe8PuVlZVUVlaSmtr9JXDnzp3Lfffdd87XGxoa3DryI+aJCIzgqoyr+XD3UnbaFrFlSx4jL2FKaUsLLFxoHNek/hnq4dtjvq0Pqj5gQvZgKIGm0HKOH4doE2dq7NyJVizzEa5f8jHlbC4+DQR79fobNwJJGwEYnaQg46+CAoL4+qiv89S6p2DCcyxYcA2/+AWkp/fsPCdPwgsvAHOMh6ffGPUN7THVD9yQcwPRQdEcZy8nYt9nwYLr+OY3e36elhZ4+WVgtpGiXSOFvbR+/XpCQkLO+z2r1b2Twfx+allvzJs3j5KSkvP+b84cY+rQE088QUlJCY8++mi3z2uz2bDb7ef9n/Qdt43oWHlm2D+MZS8vwbJlcPQoxAwup6y+AJvFxp0j3PODRHpneHLHEswx5UaQMFFZGVqxzEckhCYQZo0FazvrK0u8fv2CwnZINNZqH52oIOPP5o6faxwM+SetYZX87nc9P8czz8DJwGLI+hSbxcZ94859iCp9T1BAEPeMvMf4y9g/8L//yyUtPPTmm3CgfQukf0aANYD/GPkfbqnvQp+B7XY7NpvNLdfo1K+DTHc4fWVbb/E5Nw69kQCLHRK38MrSYlovYVuJzj0EBn3FmF82Y+AM4kPj3VilXCrXZmDBx9mw/ZiptezYAcRoxTJfYLFYGBJljMqUnfD+ymVryyogsA6HNYihA4Z6/friPkNihxjTwCxOuOwZ/vAH2Lev+++vrcUIPxOM0Zg5uXNIi9QWEv3Fd8Z2tDsM+Sfb9+3t8eaq7e3w6KPA+BcAuDH3RhLDEt1bpBcoyIhcopjgGGYPNhrujib/lSVLevb+48fhH8aS/xxLehOA24fd7s4SpRdC7CGEtqUAULDLvD6Z9vaOIKOpZT5jfLrRJ1MbWMzhw967bk0NHMCYVpY/YLimEPUBP5j4AwCsE/7AaY7y0592/70PPwzHWvZjGf0KAA9e9qD7CxSfNXTAUK7MvBKs7TDxKR5+2Ph90V3vvAMllSdci0R8b9z3PFKnpynIiPTCf4zqWN1jxF958qke/ATBGI1pbITBV2xgV30xdqvdrzah6g8SHcb0su0HzQsy1dXQ0HIKIo1HtQoy5hudcqbhf/Nm7113wwYg0Wj0H5eiaWV9wTWDrmF04mjabQ1w2dO89hqsXn3x9xUVwdNPA5OfwGlr4vL0y5maPtXj9Ypv+e/J/20cjPsD67eecC0edDHNzfDQQ8CEZ8BRz7ABw4xQ5IcUZC7g0UcfpaSkhC9/+ctmlyI+7EtDvkSEIxIi97Ji30fd/lDjdMIf/2gcR137WwBuGXYLUUFRHqpULsXAKCPIVNWaF2TKyoAY4/qxwbHEBMeYVosYXA3/CUVs2uS96xYW4lqxTI3+fYPFYuGnU41hGMcVT0FIDXfeCSdOXPg9dXVw++3QFlKN9bLfAzD/ivlaJKYfmj14trGSoqMexj/HQw8ZUw4v5qmnYOeeWiyTjcasn13xM7+9fxRkRHohKCCIe0Z1NMeNf45f/7p773v3XSguhpCkPWxqMR6h/GjSjzxUpVyqEalGkDnSXt6jIXt36tror9EY3zBsQMfucxHVrN3ivf6pgkLnmRXL1OjfZ9w09CZGJ46m2XqS8K/MZ/duuO02OH363Ne2tMCdd8L27RB40wO0WxuZkjaFqwde7f3CxXQWi4WfXP4T43jq4+w9eoQf/OCL37NtG/zP/wATf4cz6Di5cbnckneL54v1EAUZkV5yrTyTs4Q336u86BNap9OY2wyQf+9TtDnbmJ41XU9YfdC4gUaQaYssZ/9+c2ooKwPiygDIjtWKZb4gMiiSxMAsAAr2em9u2frt+yG0BitWt+z1IL7BarHy5LVPAnAq548EZm7ggw9g9mzYs+fM6+rrjZGYJUsgIH8xTVmLCLAG8ML1L/jt03Tpvdvzb2dU4iicjlqY9gv+/Gf461/P/9raWrjlFjjt2I31iscAeOTKR7BZ3buSmDcpyIj0Um5cLjMHzuxYeeZpHnroi5dBfPNN2LwZQmNPsC3YmF+m0RjflBt/ZgnmHTvMqWHHDiDOWOZ3aJxWqfIV41LGAFDdvpG6Os9f7+BBOGgxnpLkxA4lxH7+PRrEP12RcQW3DruVdtqJ+97thEbXsXw55OXBTTcZozCZmca+Y/boA4Tc/m3A+N2hUNu/WS1Wnpj5BACWCc9Byjq+/nV4++2zX3fiBFxzDZSUOAmc8wPabY1My5jm16MxoCAj4hY/nPRD42Dc7/lwzQFjp9zzOHQIOje0nfzAi5xqqSdvQB7XDr7WO4VKjwyK7liCObSGzSUnTamhrAwYsB2AvAF5ptQg55qUaQQZkjayZYvnr9e10X+sGv37pBeuf4G0iDSqT5cz7rE7mHR5E6dOwaJFsGCBsefYoNxTDP3516htPcKoxFE8fOXDZpctPuDqgVdze/7tOC1thN1zB622Wm6+GWbOhGefNRr7Bw+GtWshZOofacpahM1i4+nZT/v9aJ6CjIgbzBo0i0mpk8DeCJc/xve/37EDdxfNzfAf/2H8Msofe5KiUKPJ/0eTfuT3P0j6qvDAcELbjXX1Cyq9PyRz+jRU7W53TS3TiIzvGJN0Jsh8/r91TzAa/dUf05fFBMfw5s1vEmgLZPmBfxL6nS/z9sdVPPEEPPEEPPe3EqJ/cBVba5cTYg/hja++QWBAoNlli4944foXyIjMoN5RQep/34AtqIGPPjIenj7+uPHZI3nqB7TO/D4Aj139GCMSRphcde8pyIi4gcVi4f9e9X+N4wkv0BBWxLRp8Pe/G9PMjh415jZ/8AGEhEDuAw9x8NRBBscM5o7hd5hcvXyRlKBcALYf9v4u7uXlQORusJ8m0BZIVnSW12uQ83OFidgdrC6o9/j1CguBlPUAjEse5/HriTkmp03mn3f8k+CAYD6q+JDbVw/hndipLIyezAOl+RQeKCA2OJaP7v6I3Lhcs8sVHxIVFMXC2xYSERjBvoDlDP/NDB54eAfXXgv33NvI7c/+msMzr6e5vYk5uXPOzCTxcwoyIm4yI2sGN+TcgNPaQvjd36C+oZXbboPkZMjIMKYH2O3w85c/4q2qPwDw4pdf1BM1Hzc0zpjOtee094NM12llQ2KHaANEH5IQlkCsPRksTlbv8uzcMqcT1pXsg4j92Cw2xiaN9ej1xFxXD7yatd9ay4ysGbS0t7BqzyrW7FtDu7OdG3NvZP231zMpbZLZZYoPGpM0hqV3LCUiMILNR9byDDlsmpHI37KiefPIj2ltb+XO4Xfyt5v/1mdmgijIiLiJxWLhhetfICooirrwDYz9+TxCw5wcPAinTsGIEfDGv6r57W5jBOa7Y7/rtxtQ9SfjM43pXCcCSmhu9u61jSBjBCj1x/ieMw3/Gzh61HPXqa6GGsc6APLi8gl1hHruYuITRiSM4MO7P6T0/lJenfMqb3z1DcrmlbHwtoUMjB5odnniw6akT2HrfVu5ZtA1ABw6dYjG1kZSI1J5+Ssv89qNr+GwOUyu0n30eE/EjZLDk3nphpe4+e83s8HyB+78czN3JjxGYng8bfGbuHPR7dQ01DAyYSS/u+Z3Zpcr3XDZwKHwGRBXQkUF5HpxNkfXERn1x/ieSRnjeL/qn5BcwPr1xnK5nrB6NZBqBJlJaZd55iLicywWCzlxOeTE5ZhdiviZjKgM3rvrPU42nqT8WDnRwdGkR6b3yVF9jciIuNlNQ2/i+eufB+Cv21/my58mc9OKLMb/aQw7ju4gLSKNt259i2B7sMmVSnfkxXcEiJhdFJc0efXaXZde1oiM77kstSNUpK5j3TrPXWfVKiDFuMDE1Imeu5CI9CmRQZGMTR7LwOiBfTLEgIKMiEfcN+4+lt2zjPHJ42lztlF1ogqrxcpXh36VTd/dxOCYwWaXKN2UFJaEvS0CrG18VrbTa9d1OqG0zHlmRGaARmR8zfjk8cZB7E5WFh7z2HVWrm6F5EKgS3gSERFNLRPxlCszr2T9t9ezr3YfJTUljE4aTVxInNllSQ9ZLBYSbEPZxzo27S0B8r1y3ZoaONl2AIJqsVlsZMdke+W60n2xIbGkhgxiX8Mu1u0roL39GqxufjxYWwtbDxWBo4Fwe4RWqhIR6UIjMiIelhqRysxBMxVi/Fh2tDEaUnbMeyuXde2PGRQzSKvb+ajLs4wRklOR69m2zf3nX7sWnKmrAJiSMRmrRb+2RUQ66SeiiMhFjMswgszh9hLa2rxzzW3bUH+MH5iYOsE4SF1n9LK42cqVQLpx4qnpU91/ARERP6YgIyJyEZOHGEGmLcZYucwbiovRimV+4EzD/1pWrGx3+/k/WeaEjJUAXJ5+udvPLyLizxRkREQuIj+hI0jElrF5q3eGZLZtQ3vI+IExSWMItAZDyFE+LS5167lPnYL1Oyoh/AB2q4MJKRPcen4REX+nICMichFZUVlYnYFgb2RV8W6vXLPr1DKNyPguh83BZSnGksgHHSvZ7cbbY/VqaE02ppWNTx5HUECQ+04uItIHKMiIiFyEzWojwTYEgMIqzzf8Hz4MNfVHIewwgFaq8nFXZnX0rmSs4NNP3XfeTz4BMpYDmlYmInI+CjIiIt3QuXLZzuOeDzLGtDKjPyY9Mp1QR6jHrymXbmpGZ5BZyUcfue+8y5YBWZ8AMD1ruvtOLCLSRyjIiIh0Q+fKZTWUcPq0Z6+1bRuQuBmAkQkjPXsx6bVJqZOwWQIgci/vr92N09n7cx49Cut3VkJ0FQGWAI3IiIich4KMiEg3jM/q6FOJK6HEw4MyxcW4gszoxNGevZj0WqgjlDGJYwGoCfvELffHBx8AmcZozGWpl2lUTkTkPBRkRES6IW9AR5AZsJ2tW93wyP0LbN0KJG0CYFTiKI9eS9xj1uCrjYNBH7pletnSpWhamYjIRSjIiIh0Q05sDlZnAASd5LNtez12nbY22FLcDAOMbeIVZPzDzIEzjYOBH7H0373bT6a9Hf79XjsM/BiAqzKv6m15IiJ9koKMiEg3BAYEkmQ3RmUK9m7x2HV27YKGkFIIaCYyMJLMqEyPXUvcZ1LaJIIDQiG0ho+Lt1Jbe+nnKiiAo/ZNEHaIMEcYk9Mmu69QEZE+REFGRKSb8gcYjfe76j0XZDZvBhLPTCuzWCweu5a4j8PmlGsOFwAAHl9JREFU4KrMaQC0pn/Ie+9d+rkWLQKylwJw9cCrCQwIdEOFIiJ9j4KMiEg3XT7YCDJ1IVs4dswz1zCCjNHor2ll/mXWoFnGQfa/Wbz40s7hdMLbb+MKMtcNvs49xYmI9EEKMiIi3TQhY4RxkLiFoiLPXKNrkNGKZf7lS0O+ZBxkrOCfHx2nqann59i6Fcr3H4HUdQDMzp7txgpFRPoWBRkRkW5y7ekSU85nhac8co1Nm50akfFTg2IGkReXB9Y26hLeY8mSnp/jrbeAIUvA4mRkwkhSI1LdXqeISF+hICMi0k0JYQmEkQAWJx8XFbv9/AcOwMHTuyH4BHarnaGdSz6L37gh5wbjIOddXn21Z+9ta4O//AXIfxOAm/Nudm9xIiJ9jIKMiEgPDI02Rkk2Hdjk9nMXFuIajcmPz8dhc7j9GuJZX875snGQ/W+WftDEoUPdf+9778Geo4ddyy7fNuw2D1QoItJ3KMiIiPTA1MFjADgWuIGDB9177oICNK3Mz12WchnJ4ckQdJL2Qf/q0ajM738PDH0brG2MTRpLdmy2x+oUEekLFGRERHpgStY44yB5A+vWuffcRpA5s/Sy+B+b1cZdw+8y/jLqFZ57DlpbL/6+nTth6VJg5F8AuD3/ds8VKSLSRyjIiIj0wNikscZBfBGr1zW67bxO59kjMlqxzH/dM+oe4yB7KXuOHuKddy7+nl/9CtoHbIW0tQRYA7hrxF2eLVJEpA9QkBER6YH0yHRCLbFga2XZNvetwVxVBUdP10DUHgBGJIxw27nFu/IG5DEhZQJY22DUKzz2GLS3X/j1lZXw2mvA2D8AMCd3Dolhid4pVkTEjynIiIj0gMViYcQAY1Sm6OgG2trcc96CAlx7hwyNG0pkUKR7TiymuG/sfQBYJj5N4aZm/vSn87/O6YT//E9oC6glYMzrAHx37He9VaaIiF9TkBER6aFpg40+maaYDZSWuueca9YAqWsAmJQ6yT0nFdPcMfwOksKScIbvh/w3eegh2LPn3Ne98Qa8+y5YJz1Da0AtuXG5TM+a7v2CRUT8kIKMiEgPTUjtaPhPXeu2hv9Vq4zzAUxMneiek4ppAgMCefCyBwEIuvpRTtS2Mns2HD9+5jXLlsG3vw0E1hJ45f8CMP+K+Vgt+tUsItId+mkpItJDk9MmGwfx21ix/kSvz1dfDxs3t0HKekBBpq/43rjvERscS2N4KZHTX2T7dhg2zGjs/8EP4LrroKEBsu75f5zmOLlxudw67FazyxYR8RsKMiIiPZQQlkCiYzBYnCyvWNPr861bB+2x2yCwnnBHOHkD8txQpZgtMiiSR658BADr9PlkDDvMgQPwf/4PPPkkNDbClK+tpir5NwA8fvXj2Kw2M0sWEfErCjIiIpdgasYUAKraV1Nf37tzGdPKjEA0PmW8Psz2Id8Z+x3yBuRxvPkIcQ9+icd+W88dd8Ctt8JLCyvYPf52nDi5Z+Q9fDnny2aXKyLiVxRkREQuwcwcI8iQtsoIIr2wciWQsQKAy9Mu793JxKfYbXYW3rqQmOAYNhwo4GX7OGb996tM+9HzzK+8gn21+8iNy+XJa580u1QREb+jICMicgmmpHcEmZT1vPdhyyWfp7ERVq12QuZyAKZlTnNHeeJDcuJy+Ncd/yIhNIGyo2V8/Z2vc//S+6muqyYnNodP/uMTooKizC5TRMTvBJhdgIiIP8qNyyXMGkO9/RhLVhXyJJe2ZPKqVdAUXAER1ditdjX691ETUydSOq+UX638FQX7C2hua+bmoTdz75h7iQiMMLs8ERG/pCAjInIJrBYrV2RMY2nlIiraPuXgwUkkXsJm7B9+iGs0ZnzKeELsIe4tVHxGVFAUj8983OwyRET6DE0tExG5RNfmXGUcZC3jo48u7RwffICrP2ZahqaViYiIdJeCjIjIJboy80rjIG01S5Y29/j9hw/D5s1OyPoYgCsyrnBjdSIiIn2bgoyIyCUaFj+MKHscOBr41+YCmnuYZf75TyC+GCL3ERQQpCAjIiLSAwoyIiKXyGqxMmPwlQCciv+Y5ct79v633gKG/AuA6VnT1R8jIiLSAwoyIiK9MHvwtcbBkCW8807333fiBEZfTfZSAK7Pvt79xYmIiPRhCjIiIr3w5ZwvY8ECKYW89cFe2tq6974lS6DFdhzSPgPguuzrPFiliIhI36Mg08VHH33Evffey6RJkxgxYgRXXnkl8+bNY+PGjWaXJiI+Kj40nkmpxuaYh6LeNVYh64bXXwdyloC1jWEDhpEZlemxGkVERPoiBRnA6XTy05/+lHnz5rFp0yby8/OZOXMmSUlJrFixgvXr15tdooj4sJuGzjEOhi7ipZcu/vqKio5ll4f9DYBb8m7xXHEiIiJ9lDbEBJ599lkWLlzIjBkzePTRR4mIOLPLcl1dHcePHzexOhHxdXNy5/CjD38EmctY/PQBDh9OIj7+wq9/8UUg+BiWwR/gBG7Lv81bpYqIiPQZ/X5E5tChQ/zxj38kJSWF3/3ud2eFGIDw8HDS09NNqk5E/MGgmEFMTpsM1nba8l7nyScv/NqGBvjzn4HcxTitrYxIGEFuXK7XahUREekr+n2QWbhwIa2trdxyyy04HA6zyxERP3XPyHuMg1Gv8tTTTmpqzv+6J580NsIMnPQnAG4bptEYERGRS9Hvg8y6desAGDVqFDU1Nbz00kv8z//8D7/5zW9YuXKlydWJiL+4dditBNoCIX4bDVEFPPLIua+pqYHHHgNS1tEUvwaHzcE3R3/T67WKiIj0Bf2+R6a8vNz15wMPPEB9fb3re3/605+YMGECzz//PGFhYWaVKCJ+ICooitvyb+MvW/4CU3/Fc88t5uqrYU7HOgBtbXDXXVBXB9G3P8lx4Gv5XyMxLNHUukVERPxVvx+Rqa2tBeCxxx4jNzeXRYsWsWHDBl5++WXS0tIoKCjgZz/7mclViog/+MnlPzH2lMl9BxI3c9dd8PLLsHMn3HOPsVJZUNo2alP/AcB/TvxPkysWERHxX34/InP//fdTUVHRo/c8/vjjDB8+HID29nYAIiMjefHFFwkKCgJg4sSJPP/889xwww28//777N69m4yMjG6dv62tzXXerlpaWnpUp4j4l9y4XG7Pv503it8g5vYfceypD/jmN7s+L3Iy5MH/YuupNm7MvZFRiaNMq1VERMQTWlpaLviZ12q1YrPZ3HYtvw8y1dXVVFVVdfv1FouFxsZG199DQ0Opra3l2muvdYWYTtnZ2QwfPpyioiIKCgq6HWSef/55nn322W7XJCJ9xyNXPsLi0sUci/qY6x95lsLnHuTIERg7Fmb812s8WvoBdqudJ2Y+YXapIiIibjdhwoQLfm/evHk88MADbruW3weZxYsX9+r9ycnJ1NbWkpKSct7vp6SkUFRUxLFjx7p9zrlz53Lfffed8/WGhoYv/McVEf+XHZvNb2b9hvuX3s/7/JAnFrfzjVHf5G/b3mTuv+YC8NCUhxgUM8jkSkVERNxv/fr1hISEnPd7Vqt7u1r8Psj01tChQyktLeXkyZPn/X7n1y/0D3I+NpvtvMNmdrv90ooUEb/yvXHfY+2+tby29TV+8P4P+MH7P3B9756R9/DIVedZ0kxERKQPsNvtXvvM2++b/adPn47T6aSgoOCc7zU0NLBt2zbACDwiIt1hsVh4dc6rPDP7GdeqZPGh8Tw24zH+dMOfsFr6/Y9eERGRXuv3IzLTp09n0KBBbNq0iQULFnDHHXcAxiIAjz76KCdPniQnJ4exY8eaXKmI+BOLxcK8CfOYN2EeRxqOEBEYgcOmTXdFRETcxeJ0Op1mF2G20tJS7r77burq6sjNzSU9PZ2SkhL27t1LTEwMf/nLXxg8eHCvr9PQ0MDo0aMB2LRpU4+mq4mIiIiI+CKzPuNqfgOQm5vL4sWLufHGGzl69CjLli2jtbWVW2+9lbffftstIUZERERERNyn308t65SSksKjjz5qdhkiIiIiItINGpERERERERG/oyAjIiIiIiJ+R0FGRERERET8joKMiIiIiIj4HQUZERERERHxOwoyIiIiIiLidxRkRERERETE7yjIiIiIiIiI31GQERERERERv6MgIyIiIiIifkdBRkRERERE/I6CjIiIiIiI+B0FGRERERER8TsKMiIiIiIi4ncUZERERERExO8oyIiIiIiIiN9RkBEREREREb+jICMiIiIiIn5HQUZERERERPyOgoyIiIiIiPgdBRkREREREfE7CjIiIiIiIuJ3FGRERERERMTvKMiIiIiIiIjfUZARERERERG/oyAjIiIiIiJ+R0FGRERERET8joKMiIiIiIj4HQUZERERERHxOwoyIiIiIiLidxRkRERERETE7yjIiIiIiIiI31GQERERERERv6MgIyIiIiIifkdBRkRERERE/I6CjIiIiIiI+B0FGRERERER8TsKMiIiIiIi4ncUZERERERExO8oyIiIiIiIiN9RkBEREREREb+jICMiIiIiIn5HQUZERERERPyOgoyIiIiIiPgdBRkREREREfE7CjIiIiIiIuJ3FGRERERERMTvKMiIiIiIiIjfUZARERERERG/E2B2Ab6gubmZV155hffee4/KykpaW1uJj49n8uTJfOc73yEtLc3sEkVEREREpIt+H2Sam5u5++672bJlC5GRkVx22WU4HA62b9/OP/7xD5YuXcprr71GXl6e2aWKiIiIiEiHfh9k3nzzTbZs2cLIkSP585//TGhoKABOp5Nf/vKXvP766zz66KO89tprJlcqIiIiIiKd+n2PTGFhIRaLhXvuuccVYgAsFgsPPvggAMXFxWaVJyIiIiIi59Hvg4zD4bjoa6KiorxQiYiIiIiIdFe/DzKXX345TqeTV155hfr6etfX29vbeeqpp7BYLNx8880mVigiIiIiIp/X73tkbrjhBlauXMnSpUuZPn06Y8aMweFwsG3bNo4dO8a3vvUt5s6da3aZIiIiIiLSRb8PMlarlSeeeIKkpCReeuklli9f7vpeXl4eEydOxGKxmFihiIiIiIh8nt8Hmfvvv5+Kiooevefxxx9n+PDhANTW1nL//fdTXFzMz372M2bOnElwcDAFBQX84he/4Nvf/jb/+7//y+zZs3tdq9PpdB2fPn261+cTERERETFb18+1XT/veprfB5nq6mqqqqq6/XqLxUJjY6Pr77/61a8oLCzkpz/9KXfeeafr69OnTyc+Pp5bbrmFX//618yaNQubzdata7S1tdHe3n7O1+vq6lzHkydP7nbNIiIiIiL+oK6u7oKLaVmt1m5/nu4Oi9ObscnHtLe3M3LkSFpbW/n0009JSEg45zUzZ85k3759/Pvf/yYzM7Nb533mmWd49tln3VytiIiIiIj/mjdvHg888IDbzuf3IzK9cfToUVpaWrBYLISHh5/3NZ1fr62t7fZ5586dy3333XfO19vb26mpqWHGjBmsW7euW0s/S//V0tLChAkTWL9+PXa73exyxMfpfpHu0r0iPaH7RbrD6XRSV1fHtGnT2LBhA4GBged9ndXq3gWT+3WQiYyMxG6309raSnFxMRMmTDjr+/X19VRWVgKQnJzc7fPabLYLDpt1jvqEhobqB4J8oZaWFgBCQkJ0r8hF6X6R7tK9Ij2h+0W6q/MBfWBgoNfulX69j4zD4WDq1Kk4nU4ee+wxampqXN9rbm7m4Ycf5vTp04wdO5a4uDgTKxURERERka769YgMwE9+8hO2bt1KSUkJ1157LaNGjSIoKIiioiIOHz5MdHQ0jzzyiNllioiIiIhIF/16RAYgLS2Nd999l69//eskJiayYcMGVq1aRXBwMHfddReLFy9m0KBBZpcpIiIiIiJd9PsRGYCYmBgeeughHnroIbNLERERERGRbuj3IzIiIiIiIuJ/bA8//PDDZhfR31gsFiZMmOD2Jeik79G9Ij2h+0W6S/eK9ITuF+kub98r/XpDTBERERER8U+K1iIiIiIi4ncUZERERERExO8oyIiIiIiIiN9RkBEREREREb+jICMiIiIiIn5HQUZERERERPyOgoyIiIiIiPidALML6C+ampr4/e9/z9KlSzlw4ACRkZFMnTqV73//+yQkJJhdnnjZtm3bWL16NUVFRWzdupVDhw5hsVgoKSn5wvctXLiQBQsWsGvXLhwOByNHjuR73/seo0eP9lLl4k2NjY2sWrWKTz75hI0bN7J//35sNhvp6enMmjWLb3zjG4SEhJz3vbpX+qeXX36ZjRs3smPHDo4ePUpTUxMDBgxg/Pjx3HvvvQwZMuS879P9IidOnGD27NkcP36cjIwM3n///Qu+VvdL/3P33XdTUFBwwe//6U9/4vLLLz/n656+V7Qhphc0Nzdz9913s2XLFuLj4xk3bhzV1dVs2bKF2NhY/va3v5Gammp2meJF999/Px9//DEWiwUAp9N50SDzy1/+ktdee43g4GCmTJlCU1MTa9aswel08vTTTzNjxgxvlS9e8o9//IOf//znWCwWBg0aRHZ2NvX19WzatIn6+noGDhzI66+/TkxMzFnv073Sf02cOJHGxkZycnJcD8l27txJZWUlAQEBPPfcc0ybNu2s9+h+EYAf//jHvPvuuzidTtLT0y8YZHS/9E933303hYWFzJo165wHaBaLhW984xtkZ2ef9XWv3CtO8bjf/va3zpycHOfXvvY1Z0NDg+vrL7/8sjMnJ8d59913m1idmOHFF190Pv30085PP/3UeeTIEefw4cOdubm5F3z96tWrnTk5Oc6JEyc69+zZ4/r65s2bnfn5+c4JEyY46+rqvFG6eNGiRYuc8+fPd1ZUVJz19ZqaGueNN97ozM3Ndf7whz8863u6V/q3jRs3Opuams75+oIFC5w5OTnOKVOmONva2lxf1/0iTqfT+dlnnzlzcnKc8+fPd+bk5DhnzZp13tfpfum/7rrrLmdubq6zurq6W6/31r2iHhkPa2lpYcGCBVgsFubPn09wcLDre1//+tfJycmhoKCA7du3m1ileNu3vvUtHnjgAaZNm0ZsbOxFX//yyy9jsViYO3cuaWlprq+PHDmS22+/ndraWt566y1PliwmmDNnDo888ghZWVlnfT0uLo758+fjdDr58MMPaW1tdX1P90r/Nnr0aBwOxzlf/9rXvkZ6ejpHjx6lvLzc9XXdL9LU1MT8+fPJzs7m3nvv/cLX6n6R7vLWvaIg42EbN26krq6O9PR0cnNzz/n+NddcA8Ann3zi7dLETzQ1NbFu3ToAZs2adc73r7nmGpxOp+6hfqbz50lzczMnTpwAdK/IFwsIMNpi7XY7oPtFDM888wzV1dU88sgj2Gy2C75O94t0lzfvFTX7e1hpaSkAeXl55/3+sGHDcDqd7Nixw5tliR+prKykubmZ2NjY8y4MMWzYMADKysq8XZqYaO/evYDx4TQyMhLQvSIXtnjxYiorK8nMzCQzMxPQ/SLGZ5RXXnmFr371q4wZM4bq6uoLvlb3i4DRu3nixAmsViuZmZlcffXVJCUlnfUab94rCjIeduDAAQASExPP+/3Of+Av+uEh/dv+/fsBLri6XXBwMBEREdTW1tLQ0HDBVaykb3n11VcBuOKKK1xP2HWvSKeXXnqJ8vJyGhoaqKioYOfOnSQmJvLb3/7WtciI7pf+zel08rOf/YyIiAh+9KMfXfT1ul8E4Pe//73r2Ol08utf/5q5c+cyd+5c19e9ea8oyHhYQ0MDFouFoKCg836/8x/v1KlT3ixL/EhDQwPAWf1VnxccHExdXR2nTp3SL49+YPny5bz99tvY7XYefPBB19d1r0in/9/enYVE+ehhHH8GHBPT1BIv0rSizRQqlaEFC1KShMAWgkgFowKNughK6iLCoMEgKB0Rr4zJICIwoosEoaKFSkkMFxAMMteEESfInNT3fxHO/1hjnAM6c96Z7we8eDf4XTyMPO/66tUrvX371rucmJioyspKpaameteRl9DmdDrV2dkpu93uvar7N+QltNlsNh09elTbtm1TQkKChoaG1NTUpNraWlVXVys6OlpFRUWS/JsVnpEBABPp7e3VhQsXJEnl5eXauHFjgCfC/6P6+np1d3erpaVFDQ0NWr16tQoLC+ecTUXoGhwc1O3bt2Wz2VRQUBDocWACZ8+e1YEDB5SUlKTw8HClpKTo9OnTcjgcMgxDDodDHo/H73NRZBZZZGSkDMPQjx8/fG6fba1Lly7151gwkdkzFRMTE/PuM7uNHAW3kZERnTp1St++fVNJSYkKCwvnbCcr+F1UVJQyMzNVV1entLQ0VVVVqaOjQxJ5CWUVFRWamprS1atX/+tjyAt82bVrl9LT0+V2u9Xe3i7Jv1nh1rJFNvsA1PDwsM/tIyMjkn5d9gd8WblypaR/s/K7iYkJud1uxcTEcCk/iI2Pj+vEiRMaGhrS4cOHdfHixT/2ISuYT1hYmPLz89XV1aVnz54pPT2dvISw58+fKyYmRleuXJmzfvaM+sjIiPc2oVu3bmnFihXkBfNKSUlRZ2enRkdHJfn3fxFFZpHNviJ1vu/EdHZ2SpI2bNjgt5lgLmvWrFF4eLhcLpe+fv2qhISEOdtnM8QtRsHr+/fvOnnypD59+qR9+/bp2rVrPvcjK/ibuLg4GYYhl8slibyEMovFIrfbrdbWVp/bJycn1draKovFosnJSUnkBfNzu92S/n0mxp9Z4dayRZaRkaHo6Gj19fV5X8X8n54+fSqLxaK9e/cGYDqYwZIlS7R9+3ZJv/LyOzIU3Dwej0pLS9XR0aHs7GzdvHnT+9ap35EV/M379+9lsViUnJwsibyEsu7ubp9/zc3NkqTk5GR1d3erq6vLe3advMAXl8vlLcSznxrxZ1YoMovMarXq+PHjMgxDFRUVc+4XrK+vV09Pj2w227zfmQEkqaSkRIZhqLa2Vp8/f/aub2tr04MHD7Rs2TIdOXIkgBNiMczMzOj8+fN69+6dsrKyVF1d7f2o4XzISuj68OGDXr58KcMw5qyfmprS3bt39fjxY0VERCg/P9+7jbzgf0FeQlNbW5uam5s1MzMzZ31/f7/OnDmjiYkJ5eTkzHndsr+yYjF+/8XDgvN4PCouLlZ7e7vi4+OVlZWlwcFB7/L9+/eVlJQU6DHhRy9evFBNTY33zPrHjx9lGIa2bNni3aesrEx79uzxLtvtdjmdTkVERGjnzp36+fOn3rx5I0mqqqriLFgQcjqdun79uiwWi3JzcxUVFeVzv/LycsXGxnqXyUpoamxs1KVLlxQXF6e0tDTFxsZqbGxMPT09Gh0dVUREhCorK5WXlzfnOPKCWQMDA8rJyVFKSoqampp87kNeQs/sb0t8fLzS0tIUHR2twcFBdXZ2yuPxaP369bpz546WL18+5zh/ZIUi4ycej0d1dXV68uSJhoaGFBMTo927d+vcuXPzfjAIwauxsVGXL1/+6z52u/2P12I+evRIDQ0N6u3tVXh4uLZu3aqysrI5BQjBw+FwqKam5q/7WCwWNTc3e2//mEVWQk9/f78ePnyolpYWffnyRWNjY7JarUpMTNSOHTtUVFSkVatW+TyWvED6VWRyc3OVnJw8b5GRyEuo6e3t1b1799Te3q7h4WGNj48rMjJSa9eu1f79+3Xs2DGFh4f7PHaxs0KRAQAAAGA6PCMDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAAAwHQoMgAAAABMhyIDAAg6AwMD2rRpk4qLiwM9CgBgkVBkAAAAAJgORQYAEHQMwwj0CACARWYx+LUHAAQRh8Mhh8Mhi8XyR6E5ePCg7HZ7gCYDACyksEAPAADAQkpNTVVeXp6ampoUHx+v7Oxs77bMzMwATgYAWEhckQEABJ2BgQHl5OTIZrPJ6XQGehwAwCLgGRkAAAAApkORAQAAAGA6FBkAAAAApkORAQAAAGA6FBkAAAAApkORAQAEHavVKkmanp4O8CQAgMVCkQEABJ24uDiFhYWpr6/vj49iAgCCA9+RAQAEpdLSUj1//lzr1q3T5s2bZbValZGRoUOHDgV6NADAAqDIAACCksvl0o0bN/T69WuNjY1penpaBQUFstvtgR4NALAAKDIAAAAATIdnZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOlQZAAAAACYDkUGAAAAgOn8AwBR53UPxCILAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7eff9f3cf828>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "hhat: [-2.8695 -0.4514  3.3125]\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_style(\"white\")\n",
    "\n",
    "msd = msdDynamics()\n",
    "param = msdParam()\n",
    "\n",
    "# Input signal (sine wave)\n",
    "amplitude=5\n",
    "frequency=0.1\n",
    "\n",
    "m = 3     # number of parameters (k, b, m)\n",
    "time = 50 # secs\n",
    "n = int(time//param.Ts)  # iterations\n",
    "\n",
    "hhat = np.zeros(m)         # initial estimated parameters\n",
    "\n",
    "q = np.zeros(m)\n",
    "delta = .00001\n",
    "P = 1/delta * np.eye(m)  # initial P\n",
    "\n",
    "# Arrays for plotting\n",
    "d_arr = np.zeros(n)\n",
    "y_arr = np.zeros(n)\n",
    "t_arr = np.zeros(n)\n",
    "\n",
    "t = param.t_start  # time starts at t_start\n",
    "for i in range(n):\n",
    "    u = [amplitude*np.sin(2*np.pi*frequency*t)]  # input (force on mass)\n",
    "\n",
    "    msd.propagateDynamics(u)\n",
    "    d = msd.outputs()[0]          # system output (position of mass)\n",
    "    q = np.hstack((u[0],q[:-1]))  # q update\n",
    "\n",
    "    k = P @ q / (1 + q.T @ P @ q)  # kalman gain vector\n",
    "    y = q @ hhat  # filter output\n",
    "    e = d - y     # error\n",
    "\n",
    "    hhat = hhat + k * e  # update of estimated parameters\n",
    "    P = P - k * q @ P    # P update\n",
    "\n",
    "    # Saving values for plot\n",
    "    d_arr[i] = d\n",
    "    y_arr[i] = y\n",
    "    t_arr[i] = t\n",
    "\n",
    "    t = t + param.Ts  # advance time by Ts\n",
    "\n",
    "\n",
    "fig = plt.figure(dpi=150)\n",
    "plt.plot(t_arr, d_arr, label='d')\n",
    "plt.plot(t_arr, y_arr, label='y')\n",
    "plt.xlabel('t')\n",
    "plt.ylabel('output')\n",
    "plt.legend(loc='upper left')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "print('\\nhhat:',hhat)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The filter's output $y$ very quickly converges to match the system's output $d$ with only minor deviations thereafter. Note that the filter coefficents $\\hat{h}$ differ from the system's parameters, yet the filter output is still very close to that of the system."
   ]
  },
  {
   "source": [
    "## Homework Problem\n",
    "\n",
    "### Background\n",
    "Recursive Least Squares is a powerful tool that can be used to model an unknown system based solely on the inputs and outputs of that system. One application of this is in filtering audio. In this problem you will be taking an input audio clip that is all of the notes of the song “Mary Had a Little Lamb” being played on the piano at the same time for the length of the song, and an output that is an audio clip of the whole song “Mary Had a Little Lamb”. Given this input and output data, build a system that filters the input clip to sound like the output clip. \n",
    "\n",
    "\n",
    "### Other information\n",
    "There is some starter code shown below. The input audio clip is called mary_had_a_little_lamb_in.wav, and the output is mary_had_a_little_lamb_out.wav.\n",
    "The starter code given will read in the .wav files, set up a few helper arrays and variables to make things a little easier, and will write the output to a .wav file called generated_out.wav.\n",
    "\n",
    "HINT: The filtering will need to be adaptive. This is because the audio output changes over time while the input to the system stays consistent. This means that the recursive least squares algorithm will need to include a forgetting factor to place greater weight on newer input/output pairs.\n",
    "\n",
    "### Input and output files:\n",
    "[mary_had_a_little_lamb_in.wav](files/mary_had_a_little_lamb_in.wav)\n",
    "\n",
    "[mary_had_a_little_lamb_out.wav](files/mary_had_a_little_lamb_out.wav)\n"
   ],
   "cell_type": "markdown",
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io.wavfile as wav\n",
    "\n",
    "# get the input and output files\n",
    "[rateIn, inWav] = wav.read(\"mary_had_a_little_lamb_in.wav\", 'r')\n",
    "[rateOut, expectedOutWav] = wav.read(\"mary_had_a_little_lamb_out.wav\", 'r')\n",
    "\n",
    "# number of taps in the filter.\n",
    "numberOfParameters = 5\n",
    "\n",
    "n = expectedOutWav.size   # length of the expected output\n",
    "\n",
    "# convenience array used to shift through input. The beginning is padded with zeros,\n",
    "# and the end is also padded with zeros.\n",
    "fn = np.hstack((np.zeros(numberOfParameters - 1),\n",
    "                inWav, np.zeros(numberOfParameters - 1)))\n",
    "\n",
    "# impulse response of system we are generating\n",
    "h = np.zeros(numberOfParameters)\n",
    "\n",
    "# preallocate the output that we will be generating\n",
    "generated_output = np.zeros(n)\n",
    "\n",
    "'''\n",
    "WRITE YOUR CODE HERE.\n",
    "\n",
    "- fn is an array of inputs that has already been zero padded\n",
    "- h is the array for the system coefficients that your code will be using to \n",
    "  approximate the system\n",
    "- generated_output is an array that you should write each output of your system to.\n",
    "  This will be saved to a .wav file that you can listen to once the system is done\n",
    "  processing.\n",
    "\n",
    "You will need to:\n",
    "1. set up P\n",
    "2. loop through the inputs and outputs to compute P, k, h, and the generated_output\n",
    "   at each time step. Save your guessed outputs to the generated_output array for\n",
    "   it to be written to a wave file at the end of the processing\n",
    "'''\n",
    "\n",
    "# write the output to a wav file.\n",
    "generated_output = generated_output.astype(np.int16)\n",
    "wav.write(\"generated_out.wav\", rateOut, generated_output)\n"
   ]
  },
  {
   "source": [],
   "cell_type": "markdown",
   "metadata": {}
  }
 ],
 "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}