jmbarr/01_smoother_notebook.ipynb

## 01_smoother_notebook.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "# 1 - Kalman smoother\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook introduces the most basic Kalman smoother.\n",
    "\n",
    "## Notation\n",
    "\n",
    "Let $p(\\mathbf{x}_{k})$, the probability distribution over $\\mathbf{x}_{k} \\in \\mathbb{R}^n$, represent a hidden state at some discrete point $k$. The measurement is given by $\\mathbf{z}_{k} \\in \\mathbb{R}^m$.\n",
    "\n",
    "Our goal is to infer $p(\\mathbf{x}_{k})$, given a sequence of measurements $\\mathbf{z}_{1}, ..., \\mathbf{z}_{K}$, (written as $\\mathbf{z}_{1:K}$), where $k < K$. \n",
    "\n",
    "### Smoothing\n",
    "\n",
    "Smoothing is a forward-backward process. In the forward step we undertake filtering from $k = 1...K$. We take advantage of the Stone Soup filtering algorithms (take your pick) for this. The backward step is undertaken using the Smoother classes.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A nearly-constant velocity example\n",
    "\n",
    "This is the same simple scenario from the first tutorial notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from datetime import datetime, timedelta\n",
    "\n",
    "from matplotlib.pyplot import plot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate a target\n",
    "\n",
    "We consider a 2d Cartesian scheme where the state vector is $[x \\ \\dot{x} \\ y \\ \\dot{y}]^T$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.types.groundtruth import GroundTruthPath, GroundTruthState\n",
    "from stonesoup.models.transition.linear import CombinedLinearGaussianTransitionModel, \\\n",
    "                                               ConstantVelocity\n",
    "\n",
    "# And the clock starts\n",
    "start_time = datetime.now()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fix random number generator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(1991)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "transition_model = CombinedLinearGaussianTransitionModel([ConstantVelocity(0.05),\n",
    "                                                          ConstantVelocity(0.05)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A 'truth path' is created starting at 0,0 moving to the NE\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "truth = GroundTruthPath([GroundTruthState([0, 1, 0, 1], timestamp=start_time)])\n",
    "\n",
    "for k in range(1, 21):\n",
    "    truth.append(GroundTruthState(\n",
    "        transition_model.function(truth[k-1], noise=True, time_interval=timedelta(seconds=1)),\n",
    "        timestamp=start_time+timedelta(seconds=k)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate measurements\n",
    "\n",
    "As in tutorial 1\n",
    "\n",
    "\\begin{align}H_k &= \\begin{bmatrix}\n",
    "                    1 & 0 & 0 & 0\\\\\n",
    "                    0  & 0 & 1 & 0\\\\\n",
    "                      \\end{bmatrix}\\\\\n",
    "          R &= \\begin{bmatrix}\n",
    "                  1 & 0\\\\\n",
    "                    0 & 1\\\\\n",
    "               \\end{bmatrix} \\omega\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.types.detection import Detection\n",
    "from stonesoup.models.measurement.linear import LinearGaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The linear Gaussian measurement model is set up by indicating the number of dimensions in the\n",
    "state vector and the dimensions that are measured (so specifying $H_k$) and the noise\n",
    "covariance matrix $R$.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "measurement_model = LinearGaussian(\n",
    "    ndim_state=4,  # Number of state dimensions (position and velocity in 2D)\n",
    "    mapping=(0, 2),  # Mapping measurement vector index to state index\n",
    "    noise_covar=np.array([[5, 0],  # Covariance matrix for Gaussian PDF\n",
    "                          [0, 5]])\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate the measurements\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "measurements = []\n",
    "for state in truth:\n",
    "    measurement = measurement_model.function(state, noise=True)\n",
    "    measurements.append(Detection(measurement,\n",
    "                                  timestamp=state.timestamp,\n",
    "                                  measurement_model=measurement_model))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Construct a Kalman filter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.predictor.kalman import KalmanPredictor\n",
    "predictor = KalmanPredictor(transition_model)\n",
    "\n",
    "from stonesoup.updater.kalman import KalmanUpdater\n",
    "updater = KalmanUpdater(measurement_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Run the Kalman filter\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.types.state import GaussianState\n",
    "prior = GaussianState([[0], [1], [0], [1]], np.diag([1.5, 0.5, 1.5, 0.5]), timestamp=start_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.types.hypothesis import SingleHypothesis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.types.track import Track\n",
    "track = Track()\n",
    "for measurement in measurements:\n",
    "    prediction = predictor.predict(prior, timestamp=measurement.timestamp)\n",
    "    hypothesis = SingleHypothesis(prediction, measurement)  # Group a prediction and measurement\n",
    "    post = updater.update(hypothesis)\n",
    "    track.append(post)\n",
    "    prior = track[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the resulting track, including uncertainty ellipses\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create and run a Kalman smoother\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "from stonesoup.smoother.kalman import KalmanSmoother\n",
    "smoother = KalmanSmoother(transition_model)\n",
    "smoothed_track = smoother.smooth(track)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the figure that follows the truth is represented by the dashed blue line, the measurements by large blue dots, the filtered result is orange and the smoothed is green. It's evident that the smoothed result is less affected by measurement scatter and more closely follows the ground truth. One could use metrics to quantify this if one chose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAFzCAYAAACZwbV4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeXxV1bn4/886Y0YSMhIIIQnzlISQgIioYIGqiAgqKLRFaina2ur96tVbvIoDt1q8ra21+qNapwLlipWiCIogIgUJCYQYBoFAwhTIPOckZ1i/P3YIUyIIyTkBnvfrldc5e5219n7OSQhP1lp7LaW1RgghhBBCeI/J1wEIIYQQQlxtJAETQgghhPAyScCEEEIIIbxMEjAhhBBCCC+TBEwIIYQQwsskARNCCCGE8DKLNy6ilMoHqgE34NJapymlwoClQDyQD9yttS7/rvNERETo+Pj4do1VCCGEEKItZGVllWitI1t6zSsJWJPRWuuS046fANZqrV9QSj3RdPz4d50gPj6ezMzM9oxRCCGEEKJNKKUKWnvNl0OQtwPvND1/B5jkw1iEEEIIIbzGWwmYBj5TSmUppWY3lUVrrQsBmh6jvBSLEEIIIYRPeWsIcqTW+phSKgpYo5Tac6ENmxK22QBxcXHtFZ8QQgghhNd4JQHTWh9reixSSn0IDANOKKVitNaFSqkYoKiVtguBhQBpaWmycaUQQogrmtPp5MiRIzgcDl+HIi6Qn58fsbGxWK3WC27T7gmYUioQMGmtq5uejwOeBVYAPwFeaHr8V3vHIoQQQnR0R44cITg4mPj4eJRSvg5HnIfWmtLSUo4cOUJCQsIFt/NGD1g08GHTD5EFWKy1Xq2U2gr8n1Lqp8Ah4C4vxCKEEEJ0aA6HQ5Kvy4hSivDwcIqLi79Xu3ZPwLTWB4DkFspLgZva+/pCCCHE5UaSr8vLxXy/ZCV8IYQQQpzhxIkT3HvvvSQmJjJ06FBGjBjBhx9+6PU44uPjKSkpOaNs+PDhpKSkEBcXR2RkJCkpKaSkpJCfn3/B5123bh1ff/118/GMGTNYvnx5W4V9Qby5EKsQQgghOjitNZMmTeInP/kJixcvBqCgoIAVK1acU9flcmGxeDeV2LJlCwBvv/02mZmZ/PnPf26xntvtxmw2t/jaunXriIiI4Jprrmm3OM9HesCEEEII0WzdunXYbDbmzJnTXNajRw8eeughwEh87rrrLm677TbGjRuH1prHHnuMQYMGMXjwYJYuXQrA+vXrmTBhQvM5fvnLX/L2228DRs/W008/TWpqKoMHD2bPHmN1qtLSUsaNG8eQIUP4+c9/jtYXvviBy+UiNDSUJ598kmHDhpGRkUFsbCwVFRUAfP311/zgBz8gLy+PN954gwULFpCSksKmTZsA+OKLL7j22mtJTEz0Sm+f9IAJIYQQHdiNN954Ttndd9/Ngw8+SF1dHbfccss5r8+cOZOZM2dSUlLCnXfeecZr69ev/87r7dy5k9TU1O+ss3nzZnJycggLC+ODDz4gOzubHTt2UFJSQnp6Otdff/1531dERATbtm3jL3/5Cy+99BJvvPEGzzzzDNdddx1PPfUUK1euZOHChec9z+kqKytJTU3l+eefb7VOz549uf/++4mIiODhhx8G4C9/+QtFRUX8+9//5ptvvuHuu+/mjjvu+F7X/r6kB0wIIYQQrfrFL35BcnIy6enpzWVjx44lLCwMgI0bN3LPPfdgNpuJjo7mhhtuYOvWrec97+TJkwEYOnRo8/ytDRs2MGPGDABuvfVWOnfu/L1itdlsF504TZo0CaUUSUlJHD169KLO8X1ID5gQQggA3G74yU9gzhy47jpfRyNO+q4eq4CAgO98PSIi4rw9XmcbOHAgH3zwQfPxq6++SklJCWlpac1lgYGBzc9bGya0WCx4PJ7m47MXlrXb7QCYzWZcLldz+aXcAerv739G+9NjON/CtifjgdbfU1uSHjAhhBAALF8OixbBiRO+jkT40pgxY3A4HLz22mvNZXV1da3Wv/7661m6dClut5vi4mI2bNjAsGHD6NGjB7t27aKhoYHKykrWrl173mtff/31LFq0CIBVq1ZRXl5+Se8lPj6erKwsgDOSyuDgYKqrqy/p3JdKesCEEEKgNSxYAD17wqRJvo5G+JJSiuXLl/PII4/wu9/9jsjISAIDA3nxxRdbrH/HHXewefNmkpOTUUrxu9/9ji5dugDGXLWkpCR69+7NkCFDznvtp59+mnvuuYfU1FRuuOGGS94Det68efzsZz+jS5cuDBs2rLn89ttv56677uKf//wnr7766iVd42Ipb3SztZW0tDSdmZnp6zCEEOKK43LB//wPJCTAj37k62iubrt376Z///6+DkN8Ty1935RSWVrrtJbqSw+YEEIILBZ46ilfRyHE1UPmgAkhxFVu/3744ANjEr4QwjskARNCiKvcggUwYwaUlvo6EiGuHpKACSHEVayoCN55x1h+IirK19EIcfWQBEwIIa5if/4zNDbCI4/4OhIhri6SgAkhxFWqrg5efRUmToS+fX0djRBXF7kLUgghrlDVNcc4UrSdwvLdOF0NmE1WwoJj6dElnc6dEjhwwEJYGDz6qK8jFR2NUooZM2bw3nvvAcZG1zExMQwfPpyPP/7Yx9G1v/z8fDZt2sS9997bbteQHjAhhLjCOBoq2LrrbTbmvsHRslwC7WGEB3enU0AkVXUnyNizhA3b/0xsXAF79sDIkb6OWFyKRYsgPh5MJuOxaSH5SxIYGEhubi719fUArFmzhm7dul36iS/C6dsUeUt+fj6LFy9u12tIAiaEEFeQ+voyNuf+jaq6YqJCE+kc1A2LxYYymTCbrQQHRBIVmsiJwgi+2LaY8spvuYSt94SPLVoEs2dDQYGxm0FBgXHcFknYzTffzMqVKwFYsmQJ99xzT/NrtbW1zJo1i/T0dIYMGcK//vUvwEhcRo0aRWpqKqmpqWzatAmAwsJCrr/+elJSUhg0aBBfffUVAEFBQc3nXLZsGTNnzgRg5syZ/Md//AejR4/m8ccfb/V6b7/9NpMmTeK2224jISGBP//5z/z+979nyJAhXHPNNZSVlQGQl5fHD3/4Q4YOHcqoUaPYs2dP83V+9atfce2115KYmMiyZcsAeOKJJ/jqq69ISUnhD3/4Azt37mTYsGGkpKSQlJTEvn37LvnzlQRMCCGuEG5XI9v2/gONJjQoptV6WsPzc8fz5K9mkbVvGbV1RV6MUrSluXONuXynq6szyi/VtGnT+Mc//oHD4SAnJ4fhw4c3vzZ//nzGjBnD1q1b+eKLL3jssceora0lKiqKNWvWsG3bNpYuXcqvfvUrABYvXsz48ePJzs5mx44dpKSknPf6e/fu5fPPP+d///d/W70eQG5uLosXLyYjI4O5c+cSEBDA9u3bGTFiBO+++y4As2fP5pVXXiErK4uXXnqJBx98sPk6hYWFbNy4kY8//pgnnngCgBdeeIFRo0aRnZ3NI488wuuvv86vf/1rsrOzyczMJDY29pI/X5kDJoQQV4iyqjyq6kqI6pz4nfV2bIskNyeSx/87A5PJzOETmfRLuMVLUYq2dOjQ9yv/PpKSksjPz2fJkiXccsuZPx+fffYZK1as4KWXXgLA4XBw6NAhunbtyi9/+Uuys7Mxm83s3bsXgPT0dGbNmoXT6WTSpEkXlIDdddddmM3m77wewOjRowkODiY4OJiQkBBuu+02AAYPHkxOTg41NTVs2rSJu+66q/ncDQ0Nzc8nTZqEyWRiwIABnGhlJ/oRI0Ywf/58jhw5wuTJk+ndu/cFfYbfRRIwIYS4Qhw8vplA/87nrffeWwMICXVw2x15WG3RHCrOpmfsjVitAV6IUrSluDhj2LGl8rYwceJEHn30UdavX0/paSv1aq354IMP6HvW7bPz5s0jOjqaHTt24PF48PPzA+D6669nw4YNrFy5kh/96Ec89thj/PjHP0adNv7tcDjOOFdgYOB5r7dlyxbsdnvzsclkaj42mUy4XC48Hg+hoaFkZ2e3+B5Pb9/a/tj33nsvw4cPZ+XKlYwfP5433niDMWPGtFj3QskQpBBCXAEcDRWUVR8j0K/lBOzbgk9ZlvEim7dvZsO6WO6+dy9+/m7MZitu7aKiuoX/xUWHN38+BJyVNwcEGOVtYdasWTz11FMMHjz4jPLx48fzyiuvNCcs27dvB6CyspKYmBhMJhPvvfce7qb9rQoKCoiKiuJnP/sZP/3pT9m2bRsA0dHR7N69G4/Hw4cffthqHK1d70J06tSJhIQE3n//fcBIsnbs2PGdbYKDg6murm4+PnDgAImJifzqV79i4sSJ5OTkXPD1WyM9YEIIcQVwuRyoVv6mztn/AbOylmMC7CqXpJGDueve8ObXFSZc7oYW24qObfp043HuXGPYMS7OSL5Oll+q2NhYfv3rX59T/t///d88/PDDJCUlobUmPj6ejz/+mAcffJApU6bw/vvvM3r06OZerPXr17NgwQKsVitBQUHNc7NeeOEFJkyYQPfu3Rk0aBA1NTUtxtHa9S7UokWLeOCBB3j++edxOp1MmzaN5OTkVusnJSVhsVhITk5m5syZOBwO/v73v2O1WunSpQtPtcHO9aq17raOKC0tTWdmZvo6DCGE6HBq64r46puFRIbEn/mCdvHcqtn8q9oJgBl4KCaZGdefWvyrpLKA5MTb6BKZ5L2ARat2795N//79fR2G+J5a+r4ppbK01mkt1ZchSCGEuALYrEEoFB7PmWsmdTv6fxxtcKIwki+bgmEJqWfU8eDBZg1ECOE9MgQphBBXAKs1gK5h/SiuyickMBqAkIosOLaKbY1we7dkPHtD6Wa5nj7d+zS3c7oc+FkCCe3Uw1ehC3FVkh4wIYS4QnSPTqPBaayNZGsoIr5gIe81hOEBrguezYoF6/HzTDyjTUVtIQnR6SxZYmnz1dSFEK2THjAhhLgcedzQWAbueuO5yUaIXzhdO/ejqGwX1x5/B4C/1wfTOySQrR+PwWp1c+vtB5tPUV1XTIA9hK/WD2HOz08t6HlyNXVou8ncQogzSQ+YEEJcTtwNUL0fTqyDskyo2gM1+6HyG1TJVwwM7UZy1RcE1ufzZfid5JQXMKbrSD5ZkcCYcYcI7dyAx+OivPoYoEjrey///WRQu62mLoRomfSACSHE5cJVC6WZ4HaALRRMZ/0K1xrLkX8RVbqBii63sfiEscWQ2jOe6io7P5y4jZLKAjSa2PCB9Oo+Gj97aLuupi6EaJkkYEIIcTlwNxjJFxr8IlquU1sAe/4XQpMJ7f1TMnf/F2ldkpnQMxbPL75hzI3QOfgHRIX1x27v1NysvVdTF5eX0tJSbrrpJgCOHz+O2WwmMjISgIyMDGw220Wd94033iA3N5eXX365zWK9nEkCJoQQl4OaA+BpAHt4y6+7HXi2PwYmG2V9H2dvxQlySnbzp/F/4NZrfsCtPwAY3GLT+fONOV+nD0O25Wrq4vISHh7evG3PvHnzCAoK4tFHHz2jjtYarTUmk8xkuljyyQkhREfncULtIWPYsQV1jfWUb/sNqjafrC7Tyaws4c+5/8SEoixjOF+sd7XY7qTp02HhQujRA5QyHhculAn4l5PNmzfz29/+ls2bN7fbNfbv38+gQYOYM2cOqampFBYWMnv2bNLS0hg4cCDPPvtsc90tW7YwYsQIkpOTGT58OHVnTTJcsWIFI0eOpKysrN3i7eikB0wIITo6RzGgQZnPeam2eAs1excSXb2Do10mYYocQYTWfHVsO0Mi+vG/c1PYMKqS60Z1wmq2tnqJ6dMl4eqIHn744VY3kT6psrKSnJwcPB4PJpOJpKQkQkJCWq2fkpJy0cOAu3bt4q233uL1118HjK2EwsLCcLlcjB49mjvvvJPExESmTZvGBx98QGpqKpWVlWdseL1s2TL+9Kc/8cknn3xnnFc6ScCEEKKjc1aBycqqVfDqq3D8OHTpAo8/uI1rzb8mQLvQKKo7GUOMu8oOcqSmiHTnLLIq/Rl3Zy7fFJlJjUk9z4XE5aiyshKPxwOAx+OhsrKy3RKbnj17kp6e3ny8ZMkS3nzzTVwuF8eOHWPXrl00NDQQFxdHaqrx83Z6LGvWrCEjI4PPPvuMoKCgdonxciEJmBBCdHTaw5q1Jp5/HhwOo6iwEHZs3Mh1o1woBRpFcM0eaoP68OmhzVhNFg58NJNu3SoZM0ZzouYENY01BNmu7v/0LjcX0lO1efNmbrrpJhobG7HZbCxatIgRI0a0SzwnN9cG2LdvH3/84x/JyMggNDSUGTNm4HA40FqjlGqxfa9evdi/fz/79u1jyJAh7RLj5ULmgAkhREdntvHWm+7m5AtAo9lwPAaNQgNaWagO6k928V7+dWA9vQJ6s2NTXyZN2YXJZMJisnC06qjP3oJoPyNGjGDt2rU899xzrF27tt2Sr7NVVVURHBxMp06dKCws5NNPPwVg4MCBFBQUsG3btuZ6brcbgISEBN5//32mT5/O7t27vRJnRyU9YEII0dHZIygv23dmmXKzrbgblY4AbOGRHI79MZsd8MD63+L0uNhbs4/QwWu57c4aULEEWAOocFT4Jn7R7kaMGOG1xOuk1NRUBgwYwKBBg0hMTGTkyJEA2O12lixZwgMPPIDD4cDf359169Y1txswYADvvfceU6ZMYeXKlSQkJHg17o5CEjAhhOjorCEEhQRSWeugweXXXBxsq6Ozfy1HOt9GbVBvsg6twOkx7njUeLj3yZeJiDu1fIBGez10cXmbN29e8/NevXqdcUOAUor33nuvxXbXXHMNW7ZsOaPs/vvvb34+dOhQdu3a1bbBXmZkCFIIITo6pZg2uzddwivgZBKlzfQLLwfAYY8BICWyb3MTq8lMWkwyWIw5O/XOekLsV+8dZ0J0NNIDJoQQl4Ep93bBRgIfvJPPzrxIoqNNPDjRmBTm8DMSsAZ3IwDhx+4gPO8+ku4LNxb2ApweJ906dfNN8EKIc0gCJoQQlwOluO3e/tx2m9XYfFuZcR0owVNgosrcCbvWfFawkQCzP6V/W8QdvygBk7EXZIWjgsjASDqdtv2QEMK3JAETQojLhVIQ0M0Yhaw7hKVqD257FFU1xzGZLHxxNJvuDeP51uXPrXdV4va4KW8oJ9ASyOColrchEkL4hswBE0KIy0FjBZRmQdGXUJsHSkNDMeagHoyM7seu2kqqnXWUfjmdAanl2KMOU+4oZ8eGOKaPHoa/zU58PCxa5Os3IoQAScCEEKLjqzsKxZvAVQ32SPCLBEsw1B+DTv0IDkkgq3AbweZgSjZNZOZMxdCuQzn+9Rie+WV/Dh20oTUUFBibbksSJoTvSQImhBAdWX0RlO8AezhYg5sn1VN/DLQLAuNxuF38q+BrpvQcxSfv5DL7XitRgVE8/aSVs/ZApq4O5s71/tsQbcTjhMZyqC+Eqr1Qts1Izov+bTyWfG38vFTng6PI2MZKe773ZfLz8xk0aNAZZfPmzeOll15qozdyruXLl1/Q0hSvv/4677777nfWyc7O5pNPPmmr0NqF1+aAKaXMQCZwVGs9QSmVAPwDCAO2AT/SWjd6Kx4hhOjwPG6o+AbsncF01q/r2nzjMSieVQWbqHbWcm//WxjbpRzUXmAIhw61fNrWykUH1lAKtQVGUgWAMn4mTDZQplOJuXYbw9WOIuM5gLJCYHcIiAVLgE/CPx+Xy8Xy5cuZMGECAwYM+M66c+bMOe/5srOzyczM5JZbbmmrENucN3vAfg2cvu/Ai8AftNa9gXLgp16MRQghOr7GMtCNxn+yZ/HUHASgyhrB4m9XE2IO4/O3J1LT0Bnqj4Ornri4lk/bWrnogLQHKvdAyRajN8seYQxB+0WALdRIqMx+YLY3ffmBNQjsYU31Io3j2gIo3mgkcpfoxhtv5PHHH2fYsGH06dOHr776CgC3282jjz7K4MGDSUpK4pVXXgEgKyuLG264gaFDhzJ+/HgKCwubz/Ob3/yGG264gRdffJEVK1bw2GOPkZKSQl5eHn/9619JT08nOTmZKVOmUNfUnXt6T1xLsTQ2NvLUU0+xdOlSUlJSWLp0Kb1796a4uBgwNizv1asXJSUll/xZXAqvJGBKqVjgVuCNpmMFjAGWNVV5B5jkjViEEOKyUZvfvJDqSW6Pm/zKY5QeW49L2cg6+Ckf5X9F4OFbePcfUfj5a1BmqD/B/PkQcFaHR0AAzJ/vtXcgLlXtYajOa0qkThuC/j5MFiMhswRB6VZw1Z2/zXm4XC4yMjJ4+eWXeeaZZwBYuHAhBw8eZPv27eTk5DB9+nScTicPPfQQy5YtIysri1mzZjH3tDHwiooKvvzyS+bOncvEiRNZsGAB2dnZ9OzZk8mTJ7N161Z27NhB//79efPNNy8oFpvNxrPPPsvUqVPJzs5m6tSpzJgxg0VNkx8///xzkpOTiYiIuOTP4VJ4awjyZeA/geCm43CgQmvtajo+AsgKgUIIcbrGSrCdWr3e7XGzvfhbGsuyGFGTC2iO73mNBreb42tnMnn8Luo8NXQy28FVyfTpRru5c41hx7g4I/k6WS4uA9oJJrMxzHipTFajR+3k0OR3UK0keifLJ0+eDBhbCuXn5wNGYjNnzhwsFiO1CAsLIzc3l9zcXMaOHQsYvWQxMTHN55s6dWqrMeTm5vLkk09SUVFBTU0N48ePb7FeS7GcbdasWdx+++08/PDD/O1vf+O+++5r9bre0u4JmFJqAlCktc5SSt14sriFqi1uUqaUmg3MBoiTfnMhxNVEuzn91+XBqmMU15czpH4/oFHA+1VuOqsAyg/cyITfLWdbUSnXR/fB1DTxevp0SbguawGx4Cg2vuxhRu/mxXA7jLlhwb2NnrDzCA8Pp7y8/IyysrKy5o2z7XY7AGazGZeraf9Rrc9J3LTWDBw4kM2bN7d4ncDAwBbLAWbOnMny5ctJTk7m7bffZv369S3WaymWs3Xv3p3o6GjWrVvHli1bmnvDfMkbQ5AjgYlKqXyMSfdjMHrEQpVSJxPAWOBYS4211gu11mla67TIyEgvhCuEEB2EydbcW+H2uDlYdYwwvxD8HEcBWFMLK2ohsCiVfv0rGDSgnnqXg/L68hbnjYnLkNkPwtMhMAEaKoy7YhvLwVVr3BHZEq3B3QDOanCUGMmbdkPnVOjU54KGMYOCgoiJiWHt2rWAkXytXr2a6667rtU248aN4/XXX29OgsrKyujbty/FxcXNCZjT6WTnzp0ttg8ODqa6urr5uLq6mpiYGJxO5/dOmM4+Fxibgc+YMYO7774bs/kiE9k21O4JmNb6v7TWsVrreGAasE5rPR34ArizqdpPgH+1dyxCCHFZCYg1Jl4D1c463NpNYOMJQqp28JllILcWmnADx8IzSL/7QwDsZhtFNUXgF+3DwEWbMlkgpC90GQMRwyAw3ujFcjec6h1r/iqBxlIjCbOFQ0h/iBwJkaMgIOa8lzrdu+++y/PPP09KSgpjxozh6aefpmfPnq3Wv//++4mLiyMpKYnk5GQWL16MzWZj2bJlPP744yQnJ5OSksKmTZtabD9t2jQWLFjAkCFDyMvL47nnnmP48OGMHTuWfv36fa/YR48eza5du5on4QNMnDiRmpqaDjH8CKC0bnHkr30uZgxBPtq0DEUip5ah2A7M0Fo3fFf7tLQ0nZmZ2f6BCiFER+CqhaIN4BdFmaOSrSd2kV74HiFVOTxiGccrO1cAYFYm5gyawn0DJlLTUEWYxc7gAfe1zbwh4XW7d++mf//+F1bZ427qJdWAMnq3lOXiJutf4TIzM3nkkUea79psay1935RSWVrrtJbqe3UvSK31emB90/MDwDBvXl8IIS4rlkDwi4GGYuxmf4LqDxFWkcGxLpPop/oARgJmMVkYGtUftMbhKCUodrQkX1cLkxnw/XBaR/fCCy/w2muvdYi5XyfJv1AhhOjIQgeCJYBATz0DSz/FaQ7kRPTNlDYYQ5Ns/Tm/CFtAUngvdEMpblsEXSJTfRuzEB3ME088QUFBwXfOYfM2ScCEEKIjM1khLB1q8wmtziEvdBT12sTaIxnYG7vQaePvuXOYFd1QRpH2JzH2RvxtHXO1cyHEKZKACSFER2eywsF3wa8L4QP+gxOVh9h0LAfnjjsZM24vNfYAiv3iSIi9nj4RfX0drRDiAnh1DpgQQoiLUPipsY1M+l+I7D6WmrJCnNoF30xl0qs24rteS9fgrgTaWl9TSQjRsUgCJoQQHZn2wI7fGOtAJRpb5n64Zzn2xhi626/lx7ea5IY3IS5DkoAJIURH4aoFVz3gMVY8twTB0Y+gfDuMeBfMNmoaa1i1fxU/Tr+fB++S5Eu0vdLSUm666SYAjh8/jtls5uRC6BkZGdhsF7fI7xtvvEFubi4vv/xym8V6OZMETAghfElraCiF2gPgKD1z/Sbtge3/CZ36QY97AVi5dyUOl4PpKXeSEu+bkMWVLTw8nOzsbADmzZtHUFAQjz766Bl1tNZorTGZZCr5xZJPTgghfMXjhspcKNkC7nrwjwK/yFNfpVug9iB0vxPqjO2H3t/1PgGeaA5+2XFupxe+t3nzZn7729+2uudiW9i/fz+DBg1izpw5pKamUlhYyOzZs0lLS2PgwIE8++yzzXW3bNnCiBEjSE5OZvjw4dTV1Z1xrhUrVjBy5EjKysraLd6OTnrAhBDCF7SGyp1QfxT8o89dudzTCPv/Cp0G0BA9lqrCjdQEJLBy3yc0bLuPXfWy+ObV4OGHH27ujWpNZWUlOTk5eDweTCYTSUlJhISEtFo/JSXloocBd+3axVtvvcXrr78OGAuchoWF4XK5GD16NHfeeSeJiYlMmzaNDz74gNTUVCorK5s3zAZYtmwZf/rTn/jkk0++M84rnSRgQgjhC44iqDts7NnY0kSuw8vBUUhBj5+yu3AH2uNi464VOFz1kHsnE/7LAfh5PWzR8VRWVuLxeADweDxUVla2W2LTs2dP0tPTm4+XLFnCm2++icvl4tixY+zatYuGhgbi4uJITTUWBD49ljVr1pCRkcFnn31GUFBQu8R4uZAETAghfKHmAFg7tZx8uerReW9QE9SfPSqScL8QTIHcsHYAACAASURBVMrE1uJ9WBvDCWkchqvLZhyuEfhZJAm7kl1IT9XmzZu56aabaGxsxGazsWjRIkaMGNEu8QQGnlrqZN++ffzxj38kIyOD0NBQZsyYgcPhQGuNauXukF69erF//3727dvHkCFD2iXGy4XMARNCCG9zVkNjBVhaWLG+PAeyH0M1lvFN2DgiAsIwKRMOVwNfHd+JO3cK426txqkbOVx52Puxiw5nxIgRrF27lueee461a9e2W/J1tqqqKoKDg+nUqROFhYV8+umnAAwcOJCCggK2bdvWXM/tdgOQkJDA+++/z/Tp09m9e7dX4uyopAdMCCG8zVXbcs9XeQ5sfQA8DWgU/hY/PE0vLdn7KQ53A91tXfnBLUWE2kPJr8inZ1hPTLLx9lVvxIgRXku8TkpNTWXAgAEMGjSIxMRERo4cCYDdbmfJkiU88MADOBwO/P39WbduXXO7AQMG8N577zFlyhRWrlxJQkKCV+PuKJTW2tcxXLC0tDSdmZnp6zCEEOLS1B2DihzjTsfT5b0F+/4CaDSKw10mU9x1Ejkl+7h/3fN4tAe7ycprN/+RpK7DKK4rZkzCGGzmi1uXSXRMu3fvpn///r4OQ3xPLX3flFJZWuu0lurLn01CCOFtqpU7GMOGNr/mURbKAnoDkHFiJx5t9IW5PG6yju/Aoz2YMGFu7VxCiA5NEjAhhPA2axDrvtBMmABpaTBhAqxaBXROgi7jABPlg16g0NYVAH9L0y38HhNmZWFo12FUNVTRrVM3zCZJwIS4HMkcMCGE8LJFSwN57flIzLoGrYMoLITnnzdeuznOD2yhhMaMpNPxXCoaaiioPo7JbceW8Rh/fiqE3uF9qHPVER8a79P3IYS4eNIDJoQQbWTRIoiPB5PJeFy0qOV6c+fC7iMJBNhrUcoYWnQ44NVXMe6QtAZjMVkYGtWfQIuddYczYd+t/CDwDrqG98alXQzvNpxAW2DLFxBCdHjSAyaEEG1g0SKYPRtO7rhSUGAcA0yffmbdQ4dA63D2n+hFry77Ka6KRGsTx48DzkqwGgtX+lnseNBUNFbBzinc/WQgw7pfR2f/znLnoxCXOUnAhBCiDcydeyr5Oqmuzig/OwGLizMStP0neqNRJMftICSggqgIF9QfNxZorckHW2eW712NyWMlpORmpt3RGavVa29JCNGOJAETQog2cOjQhZfPn2/0jil3LRFBJWitsFhMTLnDYawRZo+Eqm/RwIf5GxnZdRS/+ZskXwKcbic1jTU4XA6qG6upaajB4T61+rxJmfC3+BNiDyHAFoCfxY8gW9D37jHNz89nwoQJ5ObmNpfNmzePoKAgHn300bZ+WwAsX76cPn36MGDAgO+s9/rrrxMQEMCPf/zjVutkZ2dz7NgxbrnllrYOs81IAiaEEG3gZK9WS+Vnmz4dLFTx0VtbOHzURqM1njvug/TxLljXAIE9IGoUuyqOsr+qkP+XPJUfXl8HtLByvrgqlNaVUlBZQFFNEShQKCwmCzazDZMyoVCgwa3dVDgqKKotwu1xgwKryUr3Tt2JDYklwNoxf4ZcLhfLly9nwoQJ503A5syZc97zZWdnk5mZ2aETMJlEIIQQbWD+fAg46/+2gACj/BzuBqbelMnfF/vx5b+D+PhjuPlmjNXxXbVgDwdLAB8e/BIA24GJUJYFHle7vw/RsXi0hz0le9hydAtVjioiAiKIDIgkIiCCUL9QAqxGL5fdYsdusTf3eIX5hxEZGElkQCRBtiAKKgvYeGgjpXWllxzTjTfeyOOPP86wYcPo06cPX331FQBut5tHH32UwYMHk5SUxCuvvAJAVlYWN9xwA0OHDmX8+PEUFhY2n+c3v/kNN9xwAy+++CIrVqzgscceIyUlhby8PP7617+Snp5OcnIyU6ZMoa5pjH/evHm89NJLrcbS2NjIU089xdKlS0lJSWHp0qX07t2b4uJi4zP1eOjVqxclJSWX/FlcCknAhBCiDUyfDgsXQo8eRh7Vo4dxfPb8LwDqCsHjPHcvSGcNoI05YMCHeeuxFA7n89XJxmuO4nZ/H6JjOVx5mLyyPCIDIgm2B7e6yfV3sZgshPmHEWQNYuvRrdQ5687f6DxcLhcZGRm8/PLLPPPMMwAsXLiQgwcPsn37dnJycpg+fTpOp5OHHnqIZcuWkZWVxaxZs5g7d27zeSoqKvjyyy+ZO3cuEydOZMGCBWRnZ9OzZ08mT57M1q1b2bFjB/379+fNN9+8oFhsNhvPPvssU6dOJTs7m6lTpzJjxgwWNd2W/Pnnn5OcnExERMQlfw6XQoYghRCijUyf3krCdTrtgdqDYAs55yWPsxITcMhRxzcHv2Jb8R7IfZGJMyvAGgy1B8C/S8v7SIorktPtxKzMbXLXq9VsxYPHGJo8j9YSvZPlkydPBmDo0KHk5+cDRmIzZ84cLBYjtQgLCyM3N5fc3FzGjh0LGL1kMTExzeebOnVqqzHk5uby5JNPUlFRQU1NDePHj2+xXkuxnG3WrFncfvvtPPzww/ztb3/jvvvua/W63iIJmBBCeJOzCjwOMHU6o7jR7eRw/kf0BJz1J/jsRDkApr0TSRqeB5YwcBSBuw4ssv7X1SI2JJbiumKKa4sJ8w+76J0PHC4HFY4Keof1JsgWdN764eHhlJeXn1FWVlbWvHG23W7szmA2m3G5jKHxkzcCnE5rzcCBA9m8eXOL1wkMbP1neebMmSxfvpzk5GTefvtt1q9f32K9lmI5W/fu3YmOjmbdunVs2bKluTfMl2QIUgghvMnjBM7tXThQsIqEI+8BkFC0gq+PbMZa0Y8hCSEUOHdT0VBttPM4vRuv8Ck/ix/p3dJJ6JxARYMxub68vpzaxlqc7pZ/FrTWNLgaqG6opqS+hOK6YtweN6kxqfSJ6HNBw5hBQUHExMSwdu1awEi+Vq9ezXXXXddqm3HjxvH66683J0FlZWX07duX4uLi5gTM6XSyc+fOFtsHBwdTXV3dfFxdXU1MTAxOp/N7J0xnnwvg/vvvZ8aMGdx9992Yzb7fwksSMCGE8Kpz//OrbqxFlW9DaeM/rtXVLjJKC3AdGkJZ5VF2ZNk5WHWs1fbiymYxWegb0ZcxCWMY1m0Y8aHxBNmCaHA3NPeOFdcZXyV1JZTWl6LRhAeE0z+iPyO7j2RUj1HEBMec/2Kneffdd3n++edJSUlhzJgxPP300/Ts2bPV+vfffz9xcXEkJSWRnJzM4sWLsdlsLFu2jMcff5zk5GRSUlLYtGlTi+2nTZvGggULGDJkCHl5eTz33HMMHz6csWPH0q9fv+8V++jRo9m1a1fzJHyAiRMnUlNT0yGGHwGU1trXMVywtLQ0nZmZ6eswhBDi4jVWQskm8ItqLjpWU8yhw58xvOAPbK6DG4+CE1BuC903zSWkrid3z6xk7q2JqOgbz528L64ou3fvpn///hdU1+1x49bu5uG/k8tTXMxk/StdZmYmjzzySPNdm22tpe+bUipLa53WUn3pARNCCG+ydjLmcLkdZxRrpVDAW4XdcOqTZR7qI3bR2AifflLfvDyFECeZTWZsZht2ix2b2YbVbJXkqwUvvPACU6ZM4be//a2vQ2kmCZgQQniTUhCYaGy63STYFkCXym14lIWvvm66K0wr8FgJKO2Px1KHu8wfFZToo6CFuLw98cQTFBQUfOccNm+TuyCFEMLb/KOhrsAYjrSFEGz1x16zg9LA/liCLMY0r+wfEV2Ugl9jGPaQI4Q4rzF6wIQQVwTpARNCCG8zWaFzKigzNJRBeTY2ZzkVYdcQlLYZnH6YVv+BYG0lIOQwgZUD+PHsUdAGa0GJy8PlND9bXNz3S/41CyGEL1j8IXw4+HWBox+ByUZi/O0cU99iLriJmMgTROtAOlXfzC8fnsA9022+jlh4iZ+fH6WlpZKEXSa01pSWluLn5/e92skQpBBCeIv2GD1etfnQWGYco6B4M0SOYk9dHYdri2D3RF7+nyjuvGc6mCXxutrExsZy5MiR5r0LRcfn5+dHbGzs92ojCZgQQnhDYzmUZRt3P1oDwNbZGFIszYDGEghP5+P89QA88INbGTcxHHy/VqTwAavV2rzivLhySQImhBDtzVEMZZnGEhS2M7cg4vhaMPtB11v5aMsDpEb05i8/7yy/nYW4wskcMCGEaE+uWijbDtYQI9E6nccFJ9ZB5ChKnA1sPr6LgWo0zhPb4AI2TBZCXL4kARNCiPZUe9hY+8tsP/e1skxjaLLLWFYXbMKDh/fm/4yjR1zQWOr9WIUQXiMJmBBCtBePE2oPgS2k5dePr0GbA2gIS2fFwQ3YGqIZFNqP+AQzVB/wbqxCCK+SBEwIIdpLYyXgNtb7OourZCueo6so9evBmqPZrDy4EefOCdw8ttzYbshZAa5678cshPAKScCEEKK9aBct/Zp1lW7HlPUQJt1IWN0+Sou3Ued2oL+9jehB26h11p/WXghxJZIETAgh2o0Czl1Ms+rEBlRTcqW0m38f+RqTtmE+PJqBScXsLjtwWnshxJVIEjAhhGgvZhtnJ2Baaw6Yu6CVsc6EVmY+Kj1ObHAYz7z5/xHZyY/iunLqXA1gkkVYhbhSSQImhBDtxRoCJj/wNDYXubWbEr/uHOl2DwB/sVzLcUc1h2uO89z+J8kp2Ydy19Foi5BV8IW4gkkCJoQQ7UWZILhX02R8g0mZMCsT1f49APhn8RHA6Cdzul1kFe1Ge1xYgnr4ImIhhJe0ewKmlPJTSmUopXYopXYqpZ5pKk9QSm1RSu1TSi1VSsmfekKIK49/F7AGNydhJmWiR3BXKlzGHLDyhjrQCtxmrCYLA4KjCAmJJyiwmy+jFkK0M2/0gDUAY7TWyUAK8EOl1DXAi8AftNa9gXLgp16IRQghvMtkhc6pxqOjBLSbHp26YLEEUe+B3TXFdC6cQOTOJ/j9NT8lMbwfA3rcYizeKoS4YrV7AqYNNU2H1qYvDYwBljWVvwNMau9YhBDCJyz+ED4cAuOgoRy7s5Lk8AQ21EODx031lz9npG0iQ3uMZUT/6YQGhPs6YiFEO/PKdq9KKTOQBfQCXgXygAqtmxe5OQJIf7sQ4spltkFIf2NOmKMYv9pDrKoDu7LQkDeGqc/B8D7DfB2lEMJLvDIJX2vt1lqnALHAMKB/S9VaaquUmq2UylRKZRYXF7dnmEII0f5MVhZ92JXe16SyuhbiqhMICfTnphv9fR2ZEMKLvHoXpNa6AlgPXAOEKqVO9sDFAsdaabNQa52mtU6LjIz0TqBCCNFOFi2C2bNhf9kRvnVC7+p4Ghvhs898HZkQwpu8cRdkpFIqtOm5P/ADYDfwBXBnU7WfAP9q71iEEMLX5s6Fujqg16cAJDV0o77eKBdCXD28MQcsBninaR6YCfg/rfXHSqldwD+UUs8D24E3vRCLEEL41KFDTU96rSbOZCauabX75nIhxFWh3RMwrXUOMKSF8gMY88GEEOKqERcHBUcaIHEtoywBBNnrmsuFEFcPWQlfCCG8xOl28sSzJdgGfQy2Wq7RnQmw1REQAPPn+zo6IYQ3eWUZCiGEuJpprTlYcZD9pfvpMcJN8ox/sLXawkBXBM6AWhYuhOnTfR2lEMKbJAETQoh2lleWx96yvUT4R2A2mSmzbWdQ2FAo9qd37zLG3ONCfh0LcXWRIUghhGhHDa4G9pfvJzIgErPJzPr89eSV59E9MJFOERaCAqsoqinydZhCCC+TBEwIIdpRhaMCrTUmZSLnRA5PrH0CgLXHP6AkroIAawOHqw77OEohhLdJAiaEEO3Ird0ojI21swqzcHmMHdhcbhdflRVhdlXjKt0G1XngqvNlqEIIL5IETAgh2pG/xR/dtNNaSuSpXdg8TjsJdT0weRrp7BcMNXlQ9CWUZYO7wVfhCiG8RBIwIYRoR6F+oQTbgqmpO4G55gAA/fkBvLOWUd1CMelGYkO6gz0c7JHQUAKlW8Dt8HHkQoj2JAmYEEK0I6UUyZH9cFXsYl3hNygUUVkvEt6QBIEOzB4nnayBJyuDvTN4XFC2HTxu3wYvhGg3koAJIUQ7C3ZXMTKqHzsrDtG3cw++2dyXYcOL6R7aA4UHtPPMBrYQaKyAxjLfBCyEaHeSgAkhRHvyuKHmIA6znR0l+xgdNo6ykkDuuEnhbwsx6rQ03GgNhNqD3o1VCOE1koAJIUR7clWBbmTt0e14tIc7+g7no8X7mTCuwpjvBcZw49ksgdBQBq5678YrhPAKWXpZCCHayaJF8Pv/cRJpUexK/5qAHoHcGD8Qa89KKM+B/EVGxR3/BemvQ+ekc0+inYC/V+MWQrQ/6QETQoh2sGgRzJ4Nh4+AB82Jzl9jL0jngV93Y/e3flCWBbppkr3HZRy3SHktZiGE90gCJoQQ7WDuXKirA6fLSn3QERqDCrEXjOLNJbFs3BIEYUNBNf0KNlmN49Npfeo1IcQVRxIwIYRoB4cOGY+V9SEcj8gBQOWNBeDGkdXGcGNoClhDIf21c4cfXbXGumBmP2+GLYTwEknAhBCiHXSPb4BOBejobA5Fb8G/tgv1RQOwWxvpldi00r27BkIGtDz3y1ULQfFejVkI4T2SgAkhRBsrrSvlnv/agD12D/iV0Rj1DX7lCdQ4/RieVolCG8tT1B2DgG7nnqChDOwRYOvs/eCFEF4hd0EKIUQbqnfWk1WYxcTxQQSZ7Pz+/7ZSbnVgOjwZtI1rb/gUXXwE5a4HVzVowFFsJFvKBA3lYA2Czimn5ogJIa44koAJIUQbKqwpRKOxW+zcfDNs8lvB6v2KZx8vIuXpVzlRVUBVYxQh7gqjgS0UKr4Bd6Mx5Bg2BIJ7y+R7Ia5w8ueVEEK0oaLaIoKsQQDknMhh9f7VaDT/8e+X+dZdTkBUT6qC+p5q4N8N/LtA6GAwBxpfknwJccWTBEwIIdqQCRO6aQmJDQUb0BjPG5welmw4jkajzH7gaZqIH32DMRE/sBsEdoXKXGMOmBDiiiYJmBBCtKGY4BhqnbUAmJsWWjWhwG0jqn4YWmtC7cFQfwwswWANPtVYmY35XzV5vghdCOFFkoAJIUQbig6KxmqyUttYy5GKPDpZAxhtuR/eWcvIQYF0CYggyBYAdUchoOu5J7AEgaPEWIZCCHHFkgRMCCHakM1sI71bOk53I5sLsxkePQi/7F/hXz6UoYPcDArvaVSsOQAep7En5NmUSYYhhbjCSQImhBBtLNgeTKDFTpWzjpvjr2P/N90ZkVbDsJh+WM0WKM8Gx3EjCdv6wLlJmMlyao6YEOKKJAmYEEK0g88OrkWhuLf3LQzo5WHi+NOGFI9/fup5Sxtxaw/y61mIK5usAyaEEO1gdd5npEf1JTogkH++e+DMF5v3dzQZvV1nb8TtcYEl0CtxCiF8Q/7EEkKINlZWX8aWo1v4YeI4Gqqrzq3QWAnmAOg959yNuD0u425Ie7j3AhZCeJ30gAkhRBv7/MDneLSHH/a7gx/eM5jwCMWytw+eqlCZC6FJ0HPWuY2dlRDYw+gZE0JcsaQHTAgh2tjq/avp7NeZ5Ojr2JLdmdioyqZ5XYCrHqrzIHTQuQ2dNaCsENTDuwELIbxO/sQSQog2Uu+sp9HdyKr9qxjbcyzf7LBQXw/XXW8HRxFYQ6BqN+CBkNMSMO2Gxgoj+QpPP22OmBDiSiUJmBBCXKKaxhr2lOyhuLaY/Ip8jtccp294X9Z92QDYGTm2K4SaoWY/lHxtNArobiRdHiegIDAOghIk+RLiKiEJmBBCXILaxlq+PvI1ZmUmMiCSVftXATA4ajB/WldBfHwkMTEmIBr8omDX74zkK6S/MdneEgR+kbIBtxBXGUnAhBDiEuwv249C0cneCYA1eWsI9w+nprGGG28twjS2EehuVFYKyrdD5EgjARNCXLVkEr4QQlykRncjhTWFhNhDANhyZAt7SvdQWl/KAysfIGHkVkbcvpNGd6PRoL4Q6g5D+HAfRi2E6AgkARNCiIvk8rgAUEoB8PG+j5tfc7pdrN2Zg9an6lGyxXgMH+bVOIUQHY8kYEIIcZFsZhsmTLg9bsCYDwZgUiaUx8oHz90JHhM2s81oULoFlAU6D/FVyEKIDkISMCGEuEgWk4XuId2paKjAoz3sKtlFetd0Hkh7gO5fraRv2ADiw7pjObmoamkGdE4Gi79vAxdC+JwkYEIIcQkSOyfib/Hn6yNfU1JXwoQ+E5ja98cc3XgjA1NqSeycaFT0uKF0q8z/EkIAkoAJIcQlsVvsDI8dzt7SvSgUfcP6kvMNOJ0mJo7pgt1iNypW7QFXtSRgQghAEjAhhLhkNrONzUc2MyJ2BJMHTEYfSQPg2mtOW9urVCbgCyFOOW8CppT6XCmV7I1ghBDicnSs+hhZhVlM6DMBq9nKtGkmPvoIunc/rVJphrEVUac+PotTCNFxXEgP2H8Cf1BKvaWUimnvgIQQ4nLzyb5PAJjQZwIAEREwYYKx7mqz0i1G75eSgQchxAUkYFrrbVrrMcDHwGql1NNKKbmFRwghmny892N6hPRgUNQgqqrgpZcg/6AHGsqg/gRUH4CKbyAs3dehCiE6iAv6U0wZqwx+C7wGPATsU0r9qD0DE0KIy4HD5WDNgTVM6DMBpRSZW5w89hjs3ZJj9HqVZ8Phf4J2GxttV+eB2+HrsIUQPnYhc8A2AkeBPwDdgJnAjcAwpdTCC2jfXSn1hVJqt1Jqp1Lq103lYUqpNUqpfU2PnS/ljQghhDd5tIfaxlpW7VtFnbPOGH501pDx5WEA0oY6jc23/SKg/qjRKCwNavKgeCM0VvoweiGEr13IZtxzgJ1aa31W+UNKqd0X0N4F/D+t9TalVDCQpZRag5HIrdVav6CUegJ4Anj8e8QuhBBep7XmWPUx9pbuxeF28Lftf8PP4kdicDd0yRYytifRK9FBWPhpf99W5oJ/V/CPMo5d9UbvWMQIsAb75o0IIXzqQuaA5baQfJ106wW0L9Rab2t6Xg3sxuhJux14p6naO8CkC4pYCCF86ED5AbKPZ2M324n0jySzMJP0runsP7qebysOkbE9hGGptWc2qsiFkIGnji3+YLYb88KEEFelS7odR2t94PvUV0rFA0OALUC01rqw6TyFQNSlxCKEEO2tzlnH3tK9RAVGYbfYWbV/FcdrjtMzpAeRJjfZxxo5UWxhWGrdqUYNJeA4DqGDzjyZJcgYhpShSCGuShcyBNkmlFJBwAfAw1rrKnXG/dnf2W42MBsgLi6u/QIUQojzKKopwmQy8elqE39YkkPZtc+CGf7+zVKu7xxFbHgi23d8RHzQab+rjn5sPJr8zj2h2QZ1R8AW4p03IIToMLyyII1SyoqRfC3SWv+zqfjEyXXFmh6LWmqrtV6otU7TWqdFRkZ6I1whhGhRvaueL9dZmT8fymxZYHIB4PK4+TD7IDaTFY+lnqAgj9GgPAf2vWY83/N74/h0lkBoLPPiOxBCdBTtnoA1LWHxJrBba/37015aAfyk6flPgH+1dyxCCHEpAqwBvPmWC4cDKEsEBWiF0haOft2PZ5+4gc8+Om2l+6IvjeUnADwuKMs664zKKBdCXHW80QM2EvgRMEYpld30dQvwAjBWKbUPGNt0LIQQHVZ0UDRFxR5QLuhUCEBI3g8ZtvtHRJQHsW5VX+pKgozKrno48UVTSxOYLBA29MwTajeYbN57A0KIDqPd54BprTdi/J3Ykpva+/pCCNFW/Cx+dLX1p7HTWhwDluOujSYx73YsASXUaWMpw/FDtkN53f/f3p3H11UW+B//PHfLzc2+J13Tlu5t6EYBQQSEEUYFHdRBMsggQ0HRn7uCnfmJCqIoo6MiM50BEQ2iAv5EFgXLvpSlLV1DaZs2adOkSbM2uTe5y3l+f5xb0rQpQ0t7b5bv+/XK695z7r3nPs/rvJp+86xQ90t3fNfML4ITg8LFPLaqittvh+ZmKC+HL17Xw3kfnpLmWolIOmhTMhGRd8o6/Nc3u6nKi9GTv52MvSfT7emnp7+EXLuU7KwYixbEoeEP0PocTL8WpvwTTLuSx1ZVcdNN0NQE1kJTk+Xn/5Hgvj+PT3etRCQNFMBERN6p7q186OxdVH28EYyldPNl5Lecgnf3Yl58ZTrzT2rE3/kcND0GJWdBziyIu0tS3H477tixpPysTupbyrh+eShNlRGRdFIAExF5J+IR6KmDYClbQ09R5htHUc9cbCyAYw3WJphf9DjO1jshdzZMvRI8frcbErfb8YDczC6i8QCbG+fS0JCm+ohIWimAiYi8E5EmMB66Y7080fAyoW3n0N/nDm/1GMsZM17gln+8nraeIphyOTj94A25C7Em+ikvh6yMHkpzW4hEM3m1binReAZa3lBkbFIAExH531gLvTshkM+jO18g6sQIbD7nrZcXV77CPddeQYa/j2vv/DkUnQq+HIh3QbQb9m/ny9e1gDeTV+tO4ZW6U+mPBQmF4Oab01ctEUmflK2ELyIyYtkEOFHw5PHg9qcoDxUxzVtFMzB/4np+8enP4fU4xBJ+pkzYD4Ep7ur2ThT62iCzgnMvnUsTQZYvh/ZemDzZDV/V1emunIikg1rARETehrWW9kg7G9p38FzjWh7e+RwfrDyDz13nIRiE9858Dp/XwRjwGIerLzlosVVPAHwhd8V7b5Dqati5ExzHfVT4Ehm71AImInIECSfBxpaNNHbtJrOvh1VtO4jE+5mSO54Z03fxr/86kabnMpLv9WA8Pmaeduhiq/Gh94EUkTFNAUxE5AjqOurYs38PZTnlYPp4Yf395AayOHfCEt7oqGfp2dk07G2ktbuYooWfwFO8BAqqBl/ESUCm9rEVkcHUBSkiMoS4E2dH5w6KMosAWNu5hyca1zC/cBoBb4DcQBbbO+o5Ke85Xm+5AM/0Tx8evuI9ECx2uyBFRA6iACYiMoRwLIxjHbweL+v3ruezf/0KUSfOKy2bWL9vKyF/EKdjI3mZAhGZegAAIABJREFUnZTMXHL4BZw4xHoh56TUF15Ehj0FMBGRIRgMFgvA6qbVxJwYAAnHsrppLdZxKI28gcXDgvfNH/xhJ+qu/5VfBYGCVBddREYAjQETERlCViCLoDdINBGlqnSga9Hv9bN43Cns793DxL11RLPmkOHPcVu8nD6Ihd3Zj4WnQGZpGmsgIsOZWsBERIbgMR5mFM2gPdJOZ18nAB+eei53nPUlZoVyifV1UWy28cia86G/E+K94M2CoiVQfo7Cl4i8LbWAiYgcwfjc8URjPXz/6X8lP5DFZ2eci/Vl0oOH8R2ZeD0OGeULwZeR7G7MT3eRRWSEUAuYiMiRJPooitTxWstmPjrtXGaVVrG4bD7nTFyKbdpMR28+VWfNBTzQ+hL07Ut3iUVkhFAAExEZinWgfS3/r+4Z+hJRrp77USrzxlEaKsRnvFR4n2XVjnOYOMGCL9Pdeqh9tdsVKSLyv1AAExEZSn87RDu5d/uzVOaO47TygZmO8a7tFIWaaPOcOfB+bwZ4fNBbn4bCishIowAmIjKU3jpaYv38bdcrfHL6BzDGvPWSr2MVAOd+dN7gz/hzoXcXJKKpLKmIjEAKYCIih7DxCLFIC7/f8SIJm+CymR8Y/IamxyFQwLj8xsHnjQewEO1IWVlFZGTSLEgRkaS4E2d3127q9m0k2v46d2x8gOn5kzgpb+LAm/a9iu3aDBjMq5+BU+4YvAWRMeD0p7zsIjKyqAVMRAQ3fK3Zs4bafbWEfCFiiTib2+s4s+JkXmreQCTeB0Cs7g8YA8ZYd/HV9tVDXM0McU5EZIACmIgI0NjdSFukjdKsUvz+TB5vXAPAR6aejWMdtnQ0ALC/pQHHMVi87qD7wsWDL+Q47oB8EZG3oS5IERnzrLXUddRREEzu2+gJ8Mf6lynPLKCtr4txWSU09+4j6uun0LeNe1Z9jssvt274Orj70SbAeLX/o4j8r9QCJiJjXsIm6I/34/f6AXh066M09rawN9LBZ57+PhvatoGBjg0PEe7PpHfcMsy0KweHL4BoF2RNBo8/DbUQkZFEAUxExjyv8eL1eIk7cQDur70fAAvEE3FWt9TiiUcoCv+VJ7dewscvSRx+kYQ7RoysiYe/JiJyCAUwERnzjDFU5lfS2d9JNBFle/t2PMaD13jweb3MySllemQjPsJ86F/Oo7jokAAW63FXwC86BXyh9FRCREYUjQETEQEm5U2iaX8Tj2x9hN5YL18+9cv0O/3MLprBtECAsq0/JJE5DW/WVHCi7lZFiQg4MQgUQuEi8OekuxoiMkIogImIABm+DJZOWMoNT95ASaiEc6eei8d4KAmVMNt0Eaqt48Y//ZxvnZmBcaLg8ULmeAhNUPASkaOmLkgRGVNqaqCyEjwe97GmZuC1lt4WXtr1EssWL+OsyWdxduXZLBm/BGfTb9gfySY051OYkqVQeiYUnw55sxW+ROSYqAVMRMaMmhpYtgzCYfe4vt49BqiuhrtfvxuL5aqFV5GTkQxW0Q6Ce+/jrhev4LJvK2yJyPGhFjARGTOWLx8IXweEw+55x4nzy7V38v4JS5nStx32/BWanySx8Vb8nj42Ra9hwoT0lFtERh+1gInImNHQMMRJb5T9+/fw6Cv3srOrge8tuQICeYABJ07flt+xYdupnHt+FsQj4MtMdbFFZBRSABORMWPSJLfbEQBvPxS+SW75RhbN2saPNjxIjj/EORVVENvvLqba/QZZdgeTTvs0i2fvhn2NUHyqlpoQkXdNXZAiMmbcfDOEQoA3CuNeJbNoO4srdnPh+ZaXmjdwXtlJrG1YSWfrq9CxFrvtTvBmMe7ks/Bn5YEx7ubbzhALsYqIHAUFMBEZM6qrYcUKGDdzD2T0cPrUfVz5id2s9j9E1ElQVTiN7MwiNod7wfhI7FvN47UfxeneCol+d8ZjvAf6WtNdFREZ4RTARGRMOLD8xOWXQzy7jpu+keDWb6wmPr6b39avBuD7G//Mtv176YpG2L3pJXyeOC/uuQyPjULXJkhEwZcDvXXprYyIjHgKYCIy6h1YfqK+Hqy1tLRFeP7BDdRuzeDR5k3Y5PviToLVbTsp6Gskt+Nh1u+q4jP/0uu2fDlR6NnhDsKPdblbD4mIHCMFMBEZ9QYvP2HI8ycwvi5WPhNkd287AB4MPo+XM7IyWdLyB3KDncwdX0tZcIP7MV8uRPdBvI8DMyRFRI6VZkGKyKh36PITlT5DVyKBJxzmtbadnFc+l5l5FSwuquSU2Jt4ccOVx5OA7lrIme4OwMcD/fvAG0x9JURkVFELmIiMepMmDTz3eWOUGi9O9yR2TlxF3Cb47KzzuHzqGYzPLCQz3oUxYDEYjw9yZx/04ZAbwHDcZSpERI6RApiIjHpvLT8BeD0JDB58nfNoLNzCacXTCHkDdMciTLNhxvWsB8AYD1RWu61fBxgvxHogUKi1wETkXVEAE5FR78DyE5MnAxjKSuH0656n0+nhxnkXcl7FPM6rmMeO5xNgk0PyrXWXnBjEgo1C9tQU10BERhsFMBEZE6qrYedO6Ov3c89vAjxjf8vM/EmcXzYdL3Ea9mTzh5UL3KFeGDi0+xHc7sdgmdsCJiLyLiiAiciY0RfvozvawzPtrbzSspnPn3wpnqzJsH8HN9w8LRm+gPLzYM71A92P1kKsGxwHKs4HjzdtdRCR0UGzIEVk1OuN9vJG2xu09ror2P/glV+Q5cvgkqIS6K1nVe0sfvuXKl645Vtg/JA71w1d/e1gPIDjLsCaOwOC5emtjIiMCie8BcwYc5cxpsUYs/Ggc4XGmCeMMVuTjwUnuhwiMjaFY2FW7V5FV6SL4sxisPBcwwtcUDGXja1vEvZkctuvz6S8JMzSWW9AZjn4siAehp7tYAKQX+WeL1zMQDOZiMixS0UX5N3ABYecux5Yaa2dDqxMHouIHHc7OnZgrSUvmIcxhjtW3UbCJlgy/jRM5ni2t23lnluf4rG7nsQX3Q3Z0yBvFhQugIKFEG2H3l1QsAgC+ltRRI6PEx7ArLXPAu2HnL4Y+FXy+a+Aj5zocojI2BN34jR2N5IXzANgdeMqHqr7G1j45gv/wxd+aPnj6iL8gX4WTF7rznoMlrnbDEW7IN4NmWXgy4YMDbwXkeMnXYPwy6y1TQDJx9I0lUNERrGEk8DBwWPcX3U/e/I37gsGrIlTF+nipq9dw/3PVA2sbp8zAwL5EJrgtnoVVLn7P/a1pKkWIjIaDftZkMaYZcaY14wxr7W2tqa7OCIygvi9fvweP3EnTjwR5c2uWrAGHA/G8RHZ+DEycxv59S+90P6a+6FgCeScBFkT3OAF7mbcPdsH1ggTEXmX0hXA9hpjKgCSj0f809Jau8Jau8Rau6SkpCRlBRSRkc9jPEwpmEJHXwePb/8L0Ywuims/TvEbH8N/35/w7l3AuKKdlPo3QcP97ofWfg061g++kDfoDsq3idRXQkRGpXQFsIeAK5LPrwD+lKZyiMgoNylvEvkZ+fz3mjvJjlRQtP2DmFVfIrrtAkonvkZmXzFnV60GHPcDThzaVw99MQUwETlOUrEMxW+Bl4CZxpjdxpirgO8D5xtjtgLnJ49FRI47n8dHfVc9u/Y38pFx52CyuujtGkd2VgtlTpDMDB+zT5uVfHdyBfzCxUNfzGgBVhE5Pk74QqzW2k8e4aX3n+jvFhGx1nLL87cwo2gGd19yFU8UZvCz21tp3ONlfIWP666DpYvjsAYYfzFMuMgdeH+weAT8uW44ExE5DvTbRERGJWstFssjbz7Cur3ruPviu3ltyzwWzNvII48EgYO6E7esBeODOV8dmA15sPh+yF+YsrKLyOinACYio0p3fzc7O3fStL8Jx3G4YeUNTM6bzIcrL+Pk+T6mjF/Isw9vAY9/4EMdr0PenKHDV6Lf3Z4oWJy6SojIqDfsl6EQEXmn2sJtvNDwAvt691GUWURdZx2bWjfx4Zkf5ms3RGhsNNz6fcfd49GJux9K9EHXZnfV+0M5MYh2QsHJ6n4UkeNKAUxERoWEk+D15tfJy8h7a9uhu9beRWlWKScn/pm7V+Rw9TVxTju7yA1U/e0Q64bOjWDjgwOYddzg1d/uDsgPDiyBU1MDlZXg8biPNTUpr6qIjAL6k05ERoWOvg6iTpR8Xz4Af9j8B9Y0r+GTcy/jZzfOJK8wynXXNwOTITTe3V6odwe0Pg8YyJrsBq4DS01kjofsye7g+6SaGli2DMJh97i+3j0GqK5OXV1FZORTABORUaEv1ocn2ai/rnkdt75wKwAPvvEAC2dfztWfqMSbedDA+0AeBBZApAlyZ0F+cuajNwgZReDNOOw7li8fCF8HhMPueQUwETkaCmAiMioEfAEc6y6m+kDtA1jcbYPiTpzFn3yY82ZeTNA3dfCHnDi0vQxTroCcqYde8jANDUd3XkTkSDQGTERGhcLMQrweL33xPtY0rcFgMHjxGj8LyxcSd+KUZ5cP/lDH6xDvgdL3vqPvmDTp6M6LiByJApiIjAo+j495pfOoWV9Dc28z/zzrC3if+TbzNj5IeXY5JxWeRFYgy531GO2CaAc0P+F+uOSdBbCbb4ZQaPC5UMg9LyJyNNQFKSKjRpY/i99v/j2LyhdR/7uvwPPj+PrK11hYPp6KQBDa10KkGYxxP9D4MAQrIBFxN9v2hd72+gfGeS1f7nY7Tprkhi+N/xKRo6UAJiKjxi3P30JbpI0fL3qcT/1xAl/9epxPvGcBpmsD7G8BX6a7pIQxYK27/lfJe6BnO3S/CfnzIGvi235HdbUCl4i8ewpgIjJi9UR7aOltIZqI0hXp4ierfsLlVZfzn99aRFkZ/NsNYDrWQLQbMksHf7h3J8Q6oXCRO+vRiUPnenASkFOZjuqIyBiiACYiI461li1tW9jRsQOvx4vXeLnp2ZuwWL649Gv88Vy45hrIZYs73muobYQ61rqPBxZg9fjc1rHuzZCRB4GC1FVIRMYcBTARGXF2de2irqOO0lApxhg2tmzk6fqn+VTVp9gXaeKG/zuNkMcDe3dBxhH2cOxYC4FCCB00hdF43XFgPTugUAFMRE4czYIUkRHFsQ7b2rdRGCzEGMO65nXcsPIGcgI5VOz4Gi+tLKGxuxHCzYBnYMD9odpfh4IFh7/uy4a+ve6gfBGRE0QBTERGlEgsQtSJ4vf6Wb93Pdc+ci1NPU1EYhH+4z87+dsDk2npbYFIAwRyhr5I81PQ1wTBssNfMwYw7rZEIiIniAKYiIwo5qAWq5d2v0TMiQEQdyz9Fc/ymet3uu9J9IHxH36BjvWw7gb3+a4H3ONDeXzg9J2I4ouIAApgIjLCZPoyyQ5k0xfvY2vbVgAMHogHeN+MBRRNamVC7oShPxyPQO2tYOPusZOA9tUpKrmIyAAFMBEZUYwxzCicwYu7XuTp+qc5f+r5TNm1nMDvHufzV+cT8AbcLYe8meBEBz64fzusugK633AH2+N1W7oKFx/+JU4CPJkpq5OIjD2aBSkiI05eMI8Va1ZQllXGNYuv5fnd0+kf56GiopuqsqUEvAF3dmN3LXiKofEh2HyrO8B+ye1uOGtf7YavgqrBF7cWrAMZhempnIiMCQpgIjKs1dS4W//U70owuTLGrd/uZXPZ99jWvo2HL/k1p05YwpmfzyDoC5KbkTswRiyz3A1ZW/4Nmv4CRUuh6jsDy1IcGrwOiO+HUIW7ar6IyAmiACYiw1ZNDVz92Qix0BuUVG0gr2APtzzYwoaTf8onp50Na2Zz7x97+MzVrfgLpgz+cPcbsPYr0NsA06+FqVcmux7fhhOHeB8UVJ6oKomIAApgIjKM3fB/w9iKlcwpf4Ogk0FPTyl18++EeAafnv1JPvuJufj9ls9++hXYtwqyp0DODNj+37D6i+4WQ6f8wj3/TsJX/z7IPxkC+ampoIiMWQpgIjJstcXWM7tyM05vEXtzdtE251f0FW6jfO3VPBCdzta6IH++dxu+zGywWdC5EdZ82e1yrLgATr/H3VKocxNE9oA3AP5cMAfNP3Ki7l6R4G5LFBqXnsqKyJiiACYiw1Is3s+ixc+xvymfzpxdNJx+C9YTA2vIdyq55/YzmDWnmTt+2sWNy+HMqlq+9/EbCJkmmLMcTv7OQNAqqILsSujdBZFGwA58kScDcme5Y8a8GemoqoiMQQpgIjIsJSJ7+dCFvdx7dzm9RZvd8GUABzqC3cSiHmzGizQ3T+LS03/HFy/4CW3dRWzI/Tmnjq8a3MplDATy3J/cGcnlKRy3W9KbOfi9IiIpoAAmIsNHPOJ2FfbuJNCxgcXTGyn+9H5uWd/IPgNYg9/r55/PL+cPbz5NoTfMzdVf45w5T/NM7Xv59gPfIpSfz8Nnt0C0yw1ch/IG3B8RkTRSABOR9HPi0FUL4V1uq5QvC4/Hy9TCWbwRfp36PS+zIG887ymZzpIJp7E0cx/z//5G/m5SLfmhbm575Evc++JlgKG7D3eB1f7WoQOYiMgwoHZ3ETnhamqgshI8HvexpuagF50YtL/mtnwFSyFY/NaG2CXBXG6r/Qt5/hA3LvgY23/9I+ruK2P21lu4bO4LFOW0c8tD13Pvi9W4/ZNQXo4bwBL9Ka+niMg7pRYwETmhampg2TIIh93j+nr3GKD6Mgud691ZiMFi+hNROvq6cRJ9ZEUjfH3Dg9T3tnH/+z5Pw9rTePGZ6dz10Y/hwQEDiYSHwqyOt74rGITrrkseaFyXiAxjCmAickItXz4Qvg4Ih93z1R/vhMhebLCUHV27ebOzwZ2gaBO8ULeSFVuf4itzLuS8wvdwyT2VrP/BIiqy691uSmsxHj879y/GGLfl67rr4MILgf4YeILpqK6IyDuiACYiJ1RDw9uc760HXyaNPS3Utu+gNFSIx3hoCbfzs63PMiOnlIsnLOb+3zTwp89fiS8jiJn7TcBAxzq8Ey7iRxcOsaWQk4BgyYmslojIu6IAJiIn1KRJbrfjoaZP7YNIM06gkK2dWygM5uExHl5vfZNvvfxfRBJR/qdyPJN2PsB7F9SyvXMJ095zxcAq9RklEJp4+IXjPRAsAn/2ia2YiMi7oEESInJiJKLQ18ptNzcxuayFDH/fWy+FQvDdb7vHfYkY/U6cgNfP+n1bufap79HY24J1Eng71jLZ1tLkPZ0J5/yfwVsEeXyQiAz+TpuAWA9kT0tFDUVEjplawETk+IpHoKcOwrsByyXnQP5tcM/dlrVbxhPLmMa/3pjFJz7mQJs74dFaB4A/bn+KuE0A4ADPROD0TMgszCYjaA75IuMGrgMO7OWYO9vdA1JEZBhTABOR4yfWA22vAA5kFL41E/H9F8L7L7AQawW7F4qWAu7m2Jm+ILmBbLZ2NLBy1ysYDB4sAQNnBg1x6ydQMB+i7e6q9d7M5Jc5YHxu8Ip1g41D3jzInpyWqouIHA0FMBE5PpwYtL0GHi/4hlgA1Ri3CzEecUNazkw3VNkEE3x+ql/8GV6Pl9sWXEp4x+/IalrEjp2nMe4fSpk2+QNuyOrdDdHkshPRTgiUuOezKt1NtH1ZKa2yiMixUgATkeOjrxWcCARKD3vpscfg9tuhrTXKjGktfPzSVyme9SyOE6XMxPnG5pXs6tnLDxd9gk/tf56Ev5gP//Ue7vr1m0wpG+e2pAUK3J94BBK97mPpWeDPdceDiYiMIPqtJSLHR0+dG4YO8dJj66l/ajUTMudQOjlOe9Y27njUw0UxP/OWTuKWdffywO51fOfkf+ALWWG8+3dRXfN77rmjj5nlEw7/Hl+mO9OxYKHbzSkiMgIpgInIu+fEk8s/HLL21r6XOcX5AkvPdkic5WX5k5dSv288gXiIV1/oJj51B7/c8TLnlk7nwvwCPE13U993Lldd5WFm+U5g0uHfFe0Ef47b5SgiMkIpgInIcWCHPr3lP/B5428dzijbycbm6QBsYgu/efEvVGQVc2XBZxm/4zZsVojJZ36Myb5eCHdARv5Aq5p13PFf3hAULgaP/0RXSkTkhNE6YCLy7hnfwIzEA/Y8CvvfxLHu8hGO9fB6y0SC/gg9JetYVfln+p0YLeFOnvvreioydrOv5KPgRN2NtI0Xehog0Qd9+6C/HTLHuzMovdpmSERGNrWAicixs9adhRjbD54A9OyE7Eroa4ZN34eChWxs/yBTIj8hYaHf08e4op08NfEvkFzWK5ZwmLjwV2z3LiV/+j+Bz0CkCeK418mqhLw5ECwFb0b66ioichypBUxkDKmpgcpK8Hjcx5qad3Gx/nbY9yK0vghdmyDWCfvfhJYXYM1XwXix877NjAV5NAT+gfxQD5dMqeXx8pfo8vXixQuOh6AnwXuzfDSWzKfAJNyZjnlzoGgRFJwMebMga6LCl4iMKmoBExkjampg2TIIh93j+nr3GKC6+igvFm6GjrXuYPjMg5ad8Pih9kfQs422SVewqa2ecPNL2BJDYd8sVubUsi0B/zLlVJy172PKlB/x/iw4NTNOT24env1bwBuAwIF1xMzgbk0RkVFCAUxkjFi+fCB8HRAOu+ePKoDFe6HzdcgoOHwgfE8d7H2SzrxFrArHyXfepMTEIVjED/dm8ONOuDovxJfmXoSv8C+c1O1gcIfw5/btgoI50L0Viha7C7eC21wnIjLKpPU3mzHmAmPMFmPMNmPM9eksi8ho19BwdOePqHe3O0D+0PAV7YAN3yIRmszLpf9IaekpBGPtrO8N82+bn+XHO9cxvaeKO0rCTGhfx/ScguQwMIPx+Nw9HD0ZYPvdcWXgznz0BI6ygCIiw1/aWsCMMV7gduB8YDfwqjHmIWvt5nSVSWQ0mzTJ7XYc6vw7Zh0INxy+4GrHOtj4HYh2sXvecjyJAB5/BuvDvVyz6SliTgKsYduD/87GKd+kyjzgzmT0F0D5eZA3G3Lc5SnwZLiD730hN3wFCo65ziIiw1U6W8CWAtustXXW2ihwH3BxGssjMqrdfDOEQoPPhULu+XfMiYFNDN76p2M9vHIN9LrpLhrtJuANgBPj8ZY6N3wBWA8fuva/cE6aACSS2wl1Dw5f4I4Bi4ch2g3Z097a0FtEZDRJ52+28cCug453J8+JyAlQXQ0rVsDkye7wqsmT3eOjG4Bv3KUnDta+GmvdkGWtQ37vVvoTURp79/FY8xZ3gJfjJeD18NEFeWTFWnlrDQprobv28O+Ih92QF6o4tsqKiAxz6RyEb4Y4d9hy2saYZcAygElH1VciIoeqrj6GGY8H8/jBn+UujuoN0hfvp8lXziTjw9gE1nh501PEprbt3Lb21wAsn38VeyK9vLc8lzmZWbSbCe51nLgbsnJnD/6OaDdgoGiJlp4QkVErnQFsNzDxoOMJwJ5D32StXQGsAFiyZMkR9jsRkZQwBrKmQtdGIhZWNW8kZgqITL+evJ432J89mx5vCbet+i6tkS5+uuAilkw9HbwZhKNh9vXUc3LRLEzBddDbAPlz3e5Ha93ZlU7U7eac8JEhN/YWERkt0hnAXgWmG2OmAI3ApcBlaSyPiLwTmWXQW8fmlk041qEomEeEPCLZM4jE+/ja0z+gJdLOl+deTHHxQlr3N2MDOeRl5LB00jmUBLOgfx+Em8Dpd2dPYiCj2P1xEpA9Od21FBE5odIWwKy1cWPM54C/Al7gLmvtpnSVR0TeIY+fcM5sWupfosQfeqsrcXVLLd999U4ae/byg6VX8p4pFxD3F3JqRhzj9ZKRUThwjcxx7o913B/jdc/37YXcOYMH+YuIjEJp/S1nrX0UeDSdZRCRo9fjWEzeXIzHQmQPrze/zmee+ykOFp/xUly8kGDuNFrDrVC0hIyu9RDtBH/ewAKr4M5wNB6327GvDXKmuntJioiMcprfLSLHxuOD0DgoXMSrPT04yTk0FljdumXgff5sKD7dXc+rvxX69rnjveJhd8HVSAvEeiF/PuTOGhzQRERGKbXzi8hRy/JnYQ8sR2E8nDrxDO5e/xviThyfx8fiisX0x/sJ+oJkeDPcUFW4yA1dkWaIdSXXE8uAvArIKNR6XyIypiiAichRywpkURwqpru/m9yMXKrKqrjjg3ewumk1iysWU1VWxd7evcwvnY85uEXLF3K7GUVExjgFMBE5JnNK5vDy7pfp7OskLyOPqrIqqsqqiDtxWsOtlGaVMi5nXLqLKSIyLCmAicgxyQpkcdrE09jatpWm/U3JRfItXo+XaQXTmFowFa/Hm+5iiogMSwpgInLMQv4QJ5efzMzimURiEYwxZAey8WkZCRGRt6XfkiLyrgV9QYK+YLqLISIyYmjakYiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIikVE0NVFaCx+M+1tSku0QiIqnnS3cBRGTsqKmBZcsgHHaP6+vdY4Dq6vSVS0Qk1dQCJiIps3z5QPg6IBx2z4uIjCUKYCKSMg0NR3deRGS0UgATkZSZNOnozouIjFYKYCKSMjffDKHQ4HOhkHteRGQsUQATkZSproYVK2DyZDDGfVyxQgPwRWTs0SxIEUmp6moFLhERtYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJix1qa7DO+YMaYVqD9OlysG9h2na0nq6L6NTLpvI5Pu28ijeza8TLbWlgz1wogKYMeTMeY1a+2SdJdDjo7u28ik+zYy6b6NPLpnI4e6IEVERERSTAFMREREJMXGcgBbke4CyDHRfRuZdN9GJt23kUf3bIQYs2PARERERNJlLLeAiYiIiKTFmAxgxpivGmOsMaY4eWyMMT81xmwzxqxa1nHvAAAEcUlEQVQ3xixKdxllgDHmh8aYN5L35o/GmPyDXrshed+2GGM+kM5yyuGMMRck7802Y8z16S6PDM0YM9EY85QxptYYs8kY84Xk+UJjzBPGmK3Jx4J0l1UOZ4zxGmPWGmMeTh5PMca8nLxvvzPGBNJdRjncmAtgxpiJwPlAw0GnLwSmJ3+WAXekoWhyZE8A86y1VcCbwA0Axpg5wKXAXOAC4BfGGG/aSimDJO/F7bj/vuYAn0zeMxl+4sBXrLWzgdOA65L36npgpbV2OrAyeSzDzxeA2oOOfwD8OHnfOoCr0lIqeVtjLoABPwa+Dhw8+O1i4B7rWgXkG2Mq0lI6OYy19nFrbTx5uAqYkHx+MXCftbbfWrsD2AYsTUcZZUhLgW3W2jprbRS4D/eeyTBjrW2y1q5JPt+P+5/5eNz79avk234FfCQ9JZQjMcZMAD4I/E/y2ADnAvcn36L7NkyNqQBmjLkIaLTWrjvkpfHAroOOdyfPyfDzaeCx5HPdt+FN92cEMsZUAguBl4Eya20TuCENKE1fyeQIfoLbqOAkj4uAzoP+aNW/u2HKl+4CHG/GmL8B5UO8tBz4JvB3Q31siHOaHppCb3ffrLV/Sr5nOW5XSc2Bjw3xft234UP3Z4QxxmQDDwBftNZ2u40pMlwZYz4EtFhrVxtjzj5weoi36t/dMDTqApi19ryhzhtj5gNTgHXJXyoTgDXGmKW4fyFMPOjtE4A9J7iocpAj3bcDjDFXAB8C3m8H1k7RfRvedH9GEGOMHzd81VhrH0ye3muMqbDWNiWHZbSkr4QyhDOAi4wxfw8EgVzcFrF8Y4wv2Qqmf3fD1JjpgrTWbrDWllprK621lbj/OSyy1jYDDwGfSs6GPA3oOtDsLulnjLkA+AZwkbU2fNBLDwGXGmMyjDFTcCdRvJKOMsqQXgWmJ2dkBXAnTDyU5jLJEJLjhu4Eaq21/37QSw8BVySfXwH8KdVlkyOz1t5grZ2Q/D/tUuBJa2018BTwseTbdN+GqVHXAnaMHgX+HncQdxi4Mr3FkUP8HMgAnki2Xq6y1l5rrd1kjPk9sBm3a/I6a20ijeWUg1hr48aYzwF/BbzAXdbaTWkulgztDOByYIMx5vXkuW8C3wd+b4y5Cnfm+MfTVD45Ot8A7jPG3ASsxQ3XMsxoJXwRERGRFBszXZAiIiIiw4UCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiIiEiKKYCJiIiIpJgCmIiMWcaYp4wx5yef32SM+Wm6yyQiY4NWwheRsexbwHeMMaXAQuCiNJdHRMYIrYQvImOaMeYZIBs421q7P93lEZGxQV2QIjJmGWPmAxVAv8KXiKSSApiIjEnGmAqgBrgY6DXGfCDNRRKRMUQBTETGHGNMCHgQ+Iq1thb4LnBjWgslImOKxoCJiIiIpJhawERERERSTAFMREREJMUUwERERERSTAFMREREJMUUwERERERSTAFMREREJMUUwERERERSTAFMREREJMX+P5sPXuddmE5gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from stonesoup.plotter import Plotter\n",
    "plotter = Plotter()\n",
    "plotter.plot_ground_truths(truth, [0, 2], color='b')\n",
    "plotter.plot_measurements(measurements, [0, 2], color='b')\n",
    "plotter.plot_tracks(track, [0, 2], uncertainty=True, color='orange')\n",
    "plotter.plot_tracks(smoothed_track, [0, 2], uncertainty=True, color='g', label='smooth')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Key points\n",
    "\n",
    "1. Smoother be, like, smoothing.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    ".. [#] Bayesian filtering and smoothing, Sarkka, S.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}