Lesson_6_NumPy/Lesson_6_NumPy.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/How-to-Learn-to-Code/python-class/blob/master/Lesson_6_NumPy/Lesson_6_NumPy.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# NumPy\n",
    "\n",
    "The python package **NumPy** is an indespensible tool for efficient data manipulation and numerical processing. \n",
    "It is a highly curated and tested collection of functions centered around the manipulation of thier proprietary data type, the `numpy.array`. \n",
    "The `numpy.array` is similar to the standard `list` type, but it has some important differences that make it particularly useful for generic vector and matrix algebra and element-wise data manipulation.\n",
    "\n",
    "## Learning objectives \n",
    "* What is scientific computing?\n",
    "* What is NumPy?\n",
    "* What are the benefits of NumPy?\n",
    "* What is a `numpy.array`, and why use it?\n",
    "* Where is the NumPy documentation located? Right [here](https://numpy.org/doc/stable/reference/index.html)!\n",
    "\n",
    "## What's the deal with NumPy?\n",
    "\n",
    "[**NumPy**](https://numpy.org/) is the foundational Python package for scientific computing.\n",
    "Scientific computing encompasses any tools, functions, or algorithms that are used to generate or simulate mathematical models and problems. \n",
    "Python in its base form is fine for some data manipulation with support for generic algorithm construction (e.g. loops and function declaration), but as soon as you need to do any mathematics, especially on a large scale (e.g. large data sets), generic Python falls short with respect to efficient built-in functionality.\n",
    "NumPy is a heavily curated and tested package that contains thousands, if not hundreds of thousands, of useful functions for generic mathematical manipulation of your data.\n",
    "This includes general mathematics, like find the natural logarithm, as well as general statistics, like the variance and quanitles.\n",
    "These included functions are also very effienct with respect to compuational time and memory, so they should be used whenever possible.\n",
    "\n",
    "So what makes NumPy so efficient? \n",
    "Part of the efficiency comes from the the fact that the algorithms used behind the scenes are tuned to be as efficient as possible by dedicated computer scientists.\n",
    "Another reason that the code is efficient is that large chunks of the code are written in *C* and *C++* rather than in Python.\n",
    "A pre-compiled *C* or *C++* function is called by NumPy to execute your code, and the results are passed back to you in Python. \n",
    "It is not necessarily important to understand why this is more efficient, but it is important to point out that it is done in this manner because it is more efficient than strictly using Python.\n",
    "They are playing 4-dimensional chess while you are playing checkers, metaphorically, so there is no use in reinvinteing the wheel with respect to creating homemade functions.\n",
    "\n",
    "A final reason that NumPy is so efficient is the introduction of the `numpy.array`.\n",
    "The `numpy.array` is a generic n-dimensional array object.\n",
    "The easiest way to think about it is that a `numpy.array` is an n-dimensional matrix which can be used to store and manipulate your data or just do simple matrix and vector algebra.\n",
    "The reason this object is more efficient than the standard Python `list` is how the memory is stored in each case.\n",
    "In a standard `list` all of the members are stored in random spots in memory, so when they need to be accessed, the computer has to search around for them.\n",
    "In a `numpy.array` all of the data points are stored in sequential memory.\n",
    "This means that if you are at a point in your array, `my_array[a][b][c]`, and want to access the adjacent piece of data, `my_array[a][b][c+1]`, the computer can just accesses the adjacent memory address instead of searching through all of the memory for the correct piece of data.\n",
    "In the previous example `my_array` designates a 3-dimensional array as seen by the three different indexes following the `my_array` call.\n",
    "Most functions in Numpy are built around manipulating these objects, and are coded in such a way as to take advantage of the memory allocation in an optimal way, and are therefore very efficient.\n",
    "\n",
    "In order to demonstrate the power of NumPy, we can start with the standard dot product. In the code block below we generate two vectors of length 1000 and store them in `numpy.array` objects and `list` objects. We create a generic function, `my_dot_0`, that takes the dot product of two input vectors using a standard for loop. We then time the output.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "list time for my_dot_0:\n",
      "3.3 ms ± 305 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "np.array time for my_dot_0:\n",
      "4.89 ms ± 493 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%pip install numpy\n",
    "# np is the standard abbreviation for numpy\n",
    "import numpy as np \n",
    " \n",
    "# This example uses the random package so setting the seed ensures\n",
    "# we get the same result every time we run this block.\n",
    "np.random.seed(5)\n",
    "        \n",
    "# generic code block to run a dot product with any array like object\n",
    "def my_dot_0(a0, a1):\n",
    "    dot = 0\n",
    "    for i in range(len(a0)):\n",
    "        dot += a0[i] * a1[i]\n",
    "    return dot\n",
    "\n",
    "# create the vectors as numpy arrays\n",
    "# np.random.randn is the normal distribution with the number\n",
    "# in the parentheses indicating the size of the array you want\n",
    "A_arr = np.random.randn(10000)\n",
    "B_arr = np.random.randn(10000)\n",
    "\n",
    "# copy the vectors as lists\n",
    "A_list = list(A_arr)\n",
    "B_list = list(B_arr)\n",
    "\n",
    "# %timeit is a built it jupyter notebook function for timing functions easily\n",
    "# -n 10 is a statement about how many runs it will use for the time output\n",
    "# tmp is a placeholder variable to prevent output clutter\n",
    "print(\"list time for my_dot_0:\")\n",
    "%timeit -n 10 tmp = my_dot_0(A_list, B_list)\n",
    "print(\"\\nnp.array time for my_dot_0:\")\n",
    "%timeit -n 10 tmp = my_dot_0(A_arr, B_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a surprise result, seemingly contrary to everything previously stated, the `list` manipulation came out faster than the manipulation using the `numpy.array`. \n",
    "The small caveat is that the `numpy.array` object is built for speed when used with NumPy's built in functions. \n",
    "We can try the test again, but this time we use some built in functions from NumPy to clean up the function definition and improve our results, hopefully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "list time for my_dot_1:\n",
      "1.49 ms ± 171 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "np.array time for my_dot_1:\n",
      "21.8 µs ± 14.5 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "# using built in NumPy functions\n",
    "# np.multiply does element-wise multiplication\n",
    "# np.sum adds the elements of the array together \n",
    "def my_dot_1(a0, a1):\n",
    "    return np.sum(np.multiply(a0, a1))\n",
    "\n",
    "print(\"list time for my_dot_1:\")\n",
    "%timeit -n 10 tmp = my_dot_1(A_list, B_list)\n",
    "print(\"\\nnp.array time for my_dot_1:\")\n",
    "%timeit -n 10 tmp = my_dot_1(A_arr, B_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's more like it! \n",
    "Notice that the `list` speed sees some improvement using the built in NumPy functions as well, but certainly not as much as the `numpy.array`. \n",
    "In general the NumPy functions can handle generic `list` types, but the full potential is unlocked with the `numpy.array` type.\n",
    "For this final speed test, we can use the the totally pre-built `numpy.dot` function and check out the speed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "list time for np.dot:\n",
      "2.97 ms ± 399 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "np.array time for np.dot:\n",
      "The slowest run took 14.49 times longer than the fastest. This could mean that an intermediate result is being cached.\n",
      "15.5 µs ± 24.2 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "print(\"list time for np.dot:\")\n",
    "%timeit -n 10 tmp = np.dot(A_list, B_list)\n",
    "print(\"\\nnp.array time for np.dot:\")\n",
    "%timeit -n 10 tmp = np.dot(A_arr, B_arr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Will you look at that! \n",
    "The totally pre-built function is even faster than our last test with the two different functions. \n",
    "Notice that the `list` test really struggles, comparitively, with this function. \n",
    "Together, all of these speed tests should highlight:\n",
    "1. NumPy functions are very efficient.\n",
    "2. NumPy functions are relatively easy to use and clean up your code.\n",
    "3. NumPy functions work best with `numpy.array` objects.\n",
    "4. If you can find a NumPy function to handle your `numpy.array` objects, you should use it rather than trying to build your own function.\n",
    "\n",
    "## What's the deal with NumPy arrays?\n",
    "\n",
    "Now that we know why we should use the humble `numpy.array`, we should probably learn how to use this object.\n",
    "As mentioned previously is, it is very similar to you standard `list` with a few built in features that make it easier to use in a a numerical sense.\n",
    "The first feature that is crucial is that unlike a list all of the rows/columns should have the same length, much like any matrix you would normally deal with.\n",
    "This means that `numpy.array` are best suited for comparing tests or doing data manipulation on data that would fit into a traditional n-dimensional matrix. \n",
    "Below is an example of a data set that is fine in a `list` but should be avoided in a `numpy.array`.\n",
    "Python will throw a warning that this functionality is depricated (outdated, no longer supported)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0, 1, 1], [2, 3]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/marshall/.local/lib/python3.6/site-packages/ipykernel_launcher.py:3: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "my_list = [[0, 1, 1], [2, 3]]\n",
    "print(my_list)\n",
    "my_matrix = np.array(my_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we know what data works best in a `numpy.array`, we can look into how to build these arrays.\n",
    "In the first example we will look at how to make a (2,3,4) matrix that is filled with 0s for simplicity.\n",
    "In this case we use the built in `shape` attribute of the array to check its dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "a0 shape = (2, 3, 4) , a1 shape = (2, 3, 4) , a2 shape = (2, 3, 4) \n",
      "\n",
      "the above arrays have dimension 3 \n",
      "\n",
      "[[ 0  1  2  3]\n",
      " [ 4  5  6  7]\n",
      " [ 8  9 10 11]]\n"
     ]
    }
   ],
   "source": [
    "# can make an array of any size filled with \n",
    "a0 = np.zeros((2, 3, 4)) # just outputs 0s for the shape that you tell it\n",
    "a1 = np.array([[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]], [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]])\n",
    "a2 = np.array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0])\n",
    "a2.resize((2, 3, 4)) # resize changes the shape of the original array, this can be called on any array\n",
    "\n",
    "#check that the shapes are equivalent\n",
    "print('a0 shape =',a0.shape,', a1 shape =', a1.shape,', a2 shape =', a2.shape,'\\n')\n",
    "\n",
    "# sanity check on the dimensions of the above arrays\n",
    "print('the above arrays have dimension', a0.ndim, '\\n')\n",
    "\n",
    "# this how the resize works with respect to rows/columns in that it starts by filling rows\n",
    "# np.arange creates an array with 12 numbers from 0-11\n",
    "a3 = np.arange(12)\n",
    "a3.resize((3, 4))\n",
    "print(a3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At some point it might also be helpful to sort data within an array, find specific instances of things within an array, delete rows or columns, or add data. \n",
    "This can all be done using NumPy's built in functions as well. \n",
    "In the following example we will try to use NumPy to manipulate a data set in a way that would be potentially be useful for your own data. \n",
    "We can also run some simple statistics for fun!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "the location of the NaN is (3, 5) \n",
      "\n",
      "row 2 mean = 0.02053209316073448 and variance = 1.0315188879373587\n",
      "row 4 mean = -0.10027662190124971 and variance = 0.8912458143138923\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAASmElEQVR4nO3de4yddZ3H8fcXWhnlLh1rpcA0u1yK5dJ2uCklC5UFhXCrq2IrxQUbolUuuxEWkoXdjREBuW02sECFikQkIEKgC7KAsZDQZUoBgaJM3GKnqdBWoVTSUPC7f8zTbqm9zJznTM/Mj/crmcx5nvM8v/P9TSef/uZ3nvP8IjORJJVlm1YXIElqPsNdkgpkuEtSgQx3SSqQ4S5JBRrW6gIARowYkR0dHa0uQ5KGlPnz5y/PzPaNPTcowr2jo4Ourq5WlyFJQ0pEvLqp55yWkaQCGe6SVCDDXZIKNCjm3CV9MK1Zs4aenh5Wr17d6lIGtba2NkaPHs3w4cP7fI7hLqllenp62HHHHeno6CAiWl3OoJSZrFixgp6eHsaMGdPn85yWkdQyq1evZrfddjPYNyMi2G233fr9143hLqmlDPYta+RnZLhLUoEMd0kqkG+oami6bOcBaPPN5rcptYgjd0nqp4cffpjbb7+91WVsluEu6QNt8eLFHH300ey///588pOf5LrrrgNgyZIlHH300VxzzTV88Ytf5J133mHmzJlcfPHFXHHFFRxyyCEtrnzznJaRNGh0XPRgU9tbdPkJWzxm2LBhfP/732fChAm89dZbTJw4kWOPPZZFixbx5S9/ma997WtMnTqVG264genTp3PIIYdw5JFHsu+++26yzcwkM9lmm9aNnx25S/pAGzVqFBMmTABgxx13ZOzYsSxZsoTnnnuOSZMmAb2XIi5YsIADDjiAt956ixEjRvzF5YmLFi1i33335YwzzmDcuHEsXryYq6++mnHjxjFu3DiuvfbadcdeeeWVXH/99QCcf/75HHPMMQA89thjTJ06tSn9MtwlqbJo0SIWLFjAYYcdRnd3N/vssw/Lly/n4x//OMcddxznnHMOF154Ifvss89Gz3/llVf4+te/zosvvsjy5cu59dZbmTdvHk899RQ333wzCxYsAGDSpEnMnTsXgK6uLlatWsWaNWuYO3cuRx11VFP64rSMJAGrVq1iypQpXHvttey0007MmjULgBEjRnDVVVcBcPrpp2+2jb322ovDDz8cgCeeeIJTTz2V7bffHoDTTjuNuXPnMn78eCZOnMj8+fNZuXIl2223HRMmTKCrq4u5c+euG9HXZbhLDWj23DD0bX5YA2PNmjVMmTKFqVOnctpppzXcztog35Lhw4czZswYbrvtNj71qU9x4IEH8vjjj9Pd3c3YsWMbfv31OS0j6QMtMznrrLMYO3YsF1xwQdPanTRpEj/72c94++23+dOf/sS99967bg5/7fNXXXUVRx11FJMmTeLGG29k/PjxTbsdg+Eu6QPtySef5Pbbb+exxx7j4IMP5uCDD2bOnDm1250wYQJnnnkmhx56KIcddhhnn30248ePX/f8pEmTWLp0KUcccQQjR46kra3tfeFfl9MykgaNVkxNHXnkkWRm7XY6Ojp44YUX3rfvggsu2ORfA5MnT2bNmjXrtn/zm9/UrmF9jtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBdpiuEfEDyLi9Yh4Yb19H42IRyLiler7rtX+iIjrI6I7Ip6PiAkDWbwkDaQ//vGPrS6hYX0Zud8GHL/BvouARzNzb+DRahvgs8De1dcM4IbmlClJW9/555+/7vHZZ5/dwkr6b4vhnpm/BP6wwe6TgdnV49nAKevt/2H2egrYJSJGNalWSRow7733HuPHj+fEE08E4KGHHuLll1/myiuv5O2336a7u5tLLrmEU089tcWV9k2j95YZmZlLq8e/B0ZWj3cHFq93XE+1bykbiIgZ9I7u2XPPPRssQ1JRLtu5ye292edDr7vuOsaOHcvKlSuB3vu4T5s2jZkzZ/LEE08wZcoUvvnNbzJt2rTNtjMYltiDJryhmr133On3XXcy86bM7MzMzvb29rplSFLDenp6ePDBB9839fL8889z0EEHAfD0008zefJkALbddtu/OL+vS+xtjeX11mo03F9bO91SfX+92r8E2GO940ZX+yRp0DrvvPO44oor3jfaHjFiBLfccgsLFy7kpZdeYr/99mP58uVsajDalyX2tsbyems1Gu73A9Orx9OB+9bbf0Z11czhwJvrTd9I0qDzwAMP8LGPfYyJEye+b/9JJ53E7NmzGTt2LDfffDPbbLPN+5bc29CmltjbYYcd1i2xt+HyekccccS65fWaeS936MOce0T8GPgbYERE9ACXApcDd0XEWcCrwBeqw+cAnwO6gbeBrza1WklqsieffJL777+fOXPmsHr1alauXMm0adP40Y9+1K92+rLE3tZYXm+tvlwtc3pmjsrM4Zk5OjNnZeaKzJycmXtn5mcy8w/VsZmZ38jMv8rMAzKzq6nVSlKTffe736Wnp4dFixZx5513cswxx/Q72De0uSX2Bnp5vbX8hKokNdnmltgb6OX11opmLC9VV2dnZ3Z1OchXPzT7kjno12VzHRc92PSXb8USc622cOHCpk9HlGpjP6uImJ+ZnRs73pG7JBXIBbI14AZklNvW9CZbr8V/jagsjtwlqUCGu6SWGgzv+w12jfyMDHdJLdPW1saKFSsM+M3ITFasWEFbW//mIp1zl9Qyo0ePpqenh2XLlrW6lEGtra2N0aNH9+scw11Sy6z9xKaaz2kZSSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUC1wj0izo+IFyPihYj4cUS0RcSYiJgXEd0R8ZOI+FCzipUk9U3D4R4RuwPfAjozcxywLfAl4HvANZn518AfgbOaUagkqe/qTssMAz4cEcOAjwBLgWOAu6vnZwOn1HwNSVI/NRzumbkEuAr4Hb2h/iYwH3gjM9+tDusBdt/Y+RExIyK6IqJr2bJljZYhSdqIOtMyuwInA2OATwDbA8f39fzMvCkzOzOzs729vdEyJEkbUWda5jPA/2bmssxcA/wU+DSwSzVNAzAaWFKzRklSP9UJ998Bh0fERyIigMnAS8DjwOerY6YD99UrUZLUX3Xm3OfR+8bpM8CvqrZuAi4ELoiIbmA3YFYT6pQk9cOwLR+yaZl5KXDpBrt/Cxxap11JUj1+QlWSCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFqhXuEbFLRNwdES9HxMKIOCIiPhoRj0TEK9X3XZtVrCSpb+qO3K8DHsrM/YCDgIXARcCjmbk38Gi1LUnaihoO94jYGTgKmAWQme9k5hvAycDs6rDZwCn1SpQk9VedkfsYYBlwa0QsiIhbImJ7YGRmLq2O+T0wcmMnR8SMiOiKiK5ly5bVKEOStKE64T4MmADckJnjgT+xwRRMZiaQGzs5M2/KzM7M7Gxvb69RhiRpQ3XCvQfoycx51fbd9Ib9axExCqD6/nq9EiVJ/dVwuGfm74HFEbFvtWsy8BJwPzC92jcduK9WhZKkfhtW8/xvAndExIeA3wJfpfc/jLsi4izgVeALNV9DktRPtcI9M58FOjfy1OQ67UqS6vETqpJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKlCtBbI1hFy28wC0+Wbz25TUFI7cJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqUO1wj4htI2JBRDxQbY+JiHkR0R0RP4mID9UvU5LUH834ENO5wEJgp2r7e8A1mXlnRNwInAXc0ITXkTRAOi56sOltLrr8hKa3qb6rNXKPiNHACcAt1XYAxwB3V4fMBk6p8xqSpP6rO3K/Fvg2sGO1vRvwRma+W233ALtv7MSImAHMANhzzz1rllGWARlFtTW9SUmDWMMj94g4EXg9M+c3cn5m3pSZnZnZ2d7e3mgZkqSNqDNy/zRwUkR8Dmijd879OmCXiBhWjd5HA0vqlylJ6o+GR+6Z+U+ZOTozO4AvAY9l5lTgceDz1WHTgftqVylJ6peBuM79QuCCiOimdw5+1gC8hiRpM5pyP/fM/AXwi+rxb4FDm9GuJKkxfkJVkgrkSkySBoarf7WUI3dJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFchwl6QCGe6SVCDDXZIKZLhLUoEMd0kqkOEuSQUy3CWpQIa7JBXIcJekAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFajjcI2KPiHg8Il6KiBcj4txq/0cj4pGIeKX6vmvzypUk9UWdkfu7wD9k5v7A4cA3ImJ/4CLg0czcG3i02pYkbUUNh3tmLs3MZ6rHbwELgd2Bk4HZ1WGzgVNq1ihJ6qemzLlHRAcwHpgHjMzMpdVTvwdGbuKcGRHRFRFdy5Yta0YZkqRK7XCPiB2Ae4DzMnPl+s9lZgK5sfMy86bM7MzMzvb29rplSJLWUyvcI2I4vcF+R2b+tNr9WkSMqp4fBbxer0RJUn/VuVomgFnAwsy8er2n7gemV4+nA/c1Xp4kqRHDapz7aeArwK8i4tlq38XA5cBdEXEW8CrwhVoVSpL6reFwz8wngNjE05MbbVeSVJ+fUJWkAhnuklQgw12SCmS4S1KBDHdJKpDhLkkFMtwlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgQx3SSqQ4S5JBTLcJalAhrskFajOSkyDQsdFDza9zUWXn9D0NiVtXR/0bHDkLkkFGvIj9wFx2c4D0OabzW9TkjbBkbskFchwl6QCGe6SVCDDXZIK5BuqktRXQ+hiC0fuklQgw12SCmS4S1KBDHdJKtCAhHtEHB8Rv46I7oi4aCBeQ5K0aU0P94jYFvgP4LPA/sDpEbF/s19HkrRpAzFyPxTozszfZuY7wJ3AyQPwOpKkTYjMbG6DEZ8Hjs/Ms6vtrwCHZebMDY6bAcyoNvcFft2H5kcAy5tY7mBQWp/sz+BXWp9K6w/0vU97ZWb7xp5o2YeYMvMm4Kb+nBMRXZnZOUAltURpfbI/g19pfSqtP9CcPg3EtMwSYI/1tkdX+yRJW8lAhPvTwN4RMSYiPgR8Cbh/AF5HkrQJTZ+Wycx3I2Im8DCwLfCDzHyxSc33axpniCitT/Zn8CutT6X1B5rQp6a/oSpJaj0/oSpJBTLcJalAQy7cI+LfIuL5iHg2In4eEZ9odU11RMSVEfFy1ad7I2KXVtdUV0T8XUS8GBF/joghe4laabfRiIgfRMTrEfFCq2tphojYIyIej4iXqt+3c1tdU10R0RYR/xMRz1V9+peG2xpqc+4RsVNmrqwefwvYPzPPaXFZDYuIvwUeq96I/h5AZl7Y4rJqiYixwJ+B/wT+MTO7WlxSv1W30fgNcCzQQ+9VYKdn5kstLayGiDgKWAX8MDPHtbqeuiJiFDAqM5+JiB2B+cApQ/zfKIDtM3NVRAwHngDOzcyn+tvWkBu5rw32yvbA0PrfaQOZ+fPMfLfafIrezwUMaZm5MDP78onjway422hk5i+BP7S6jmbJzKWZ+Uz1+C1gIbB7a6uqJ3utqjaHV18NZdyQC3eAiPhORCwGpgL/3Op6mujvgf9qdRECekNi8XrbPQzx4ChZRHQA44F5LS6ltojYNiKeBV4HHsnMhvo0KMM9Iv47Il7YyNfJAJl5SWbuAdwBzNx8a623pf5Ux1wCvEtvnwa9vvRJ2hoiYgfgHuC8Df6yH5Iy873MPJjev+IPjYiGptAG5QLZmfmZPh56BzAHuHQAy6ltS/2JiDOBE4HJOUTeBOnHv9FQ5W00hoBqXvoe4I7M/Gmr62mmzHwjIh4Hjgf6/Sb4oBy5b05E7L3e5snAy62qpRki4njg28BJmfl2q+vROt5GY5Cr3nycBSzMzKtbXU8zRET72ivmIuLD9L6h31DGDcWrZe6h9xbBfwZeBc7JzCE7ooqIbmA7YEW166mhfPUPQEScCvw70A68ATybmce1tKgGRMTngGv5/9tofKe1FdUTET8G/obe28m+BlyambNaWlQNEXEkMBf4Fb15AHBxZs5pXVX1RMSBwGx6f+e2Ae7KzH9tqK2hFu6SpC0bctMykqQtM9wlqUCGuyQVyHCXpAIZ7pJUIMNdkgpkuEtSgf4PkqFwBTUfjmgAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# generate a 5x30 array full or of random numbers\n",
    "a0 = np.random.randn(5, 300)\n",
    "\n",
    "# pretend that one of the data points was unfilled, or a NaN\n",
    "a0[3][5] = np.nan\n",
    "\n",
    "# we want to be able to find that nan\n",
    "# np.isnan sets all instances of nan to true within the array, everything else is false\n",
    "# np.where generates a tuple of arrays that show the indices of where the trues are located\n",
    "# the tuple that np.where generates will always have the same dimension as the array queried\n",
    "where_nan = np.where(np.isnan(a0))\n",
    "# we are using NaN in this example but you could use any generic conditional to check\n",
    "where_negative = np.where(a0 < 0.)\n",
    "\n",
    "# here is a way you could print out the ordered pairs of where the nans were located\n",
    "for i in range(len(where_nan[0])):\n",
    "    print('the location of the NaN is',(where_nan[0][i],where_nan[1][i]),'\\n')\n",
    "    \n",
    "# since the NaN would mess up our data we want to delete the row it is located in\n",
    "row_delete = np.delete(a0, where_nan[0], 0)\n",
    "# we could do the same thing if we wanted the column gone instead\n",
    "col_delete = np.delete(a0, where_nan[0], 1)\n",
    "# we could also replace the nan with a 0 if we knew what that value should be \n",
    "# the tilde (~) flips the boolean values since np.where replaces false not true\n",
    "# unlike the previous np.where usage, this outputs our original array with replacements\n",
    "# instead of a tuple. This has to do with the arguments fed into the function. \n",
    "# The NumPy documentation is wonderful and should be consulted regularly.\n",
    "replace_nan = np.where(~np.isnan(a0), a0, 0.)\n",
    "\n",
    "# let's plot some of the data we generated\n",
    "# first we import the plotting package\n",
    "%pip install matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# we then want to plot the distribution of the 2nd and 4th rows\n",
    "# we can take extract the rows using standard indexing\n",
    "# we then have to take the transpose of the data so that the plt.hist function will read the data correctly\n",
    "hist = plt.hist(np.transpose(replace_nan[[1, 3]]), [-3, -2, -1, 0, 1, 2, 3])\n",
    "plt.legend([\"$2^{nd}$ row\",\"$4^{th}$ row\"])\n",
    "\n",
    "# we can output the mean and variance as well to see if it matches the expected 0,1\n",
    "mean = np.average(replace_nan[[1 ,3]], axis = 1)\n",
    "variance = np.var(replace_nan[[1, 3]], axis = 1)\n",
    "rows = [2, 4]\n",
    "for i in range(2):\n",
    "    print('row', rows[i], 'mean =',mean[i], 'and variance =', variance[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In class exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this exercise we will walk through some of the stuff we have covered previously and introduce some new helper functions that should be useful for most people.\n",
    "In general we will have an \"expected\" response curve, we will generate some data with noise, and then we will compare the noisy data to the expected response.\n",
    "We are going to pretend that we are doing some sampling every 0.01s between 0.0s and 1.0s, and this will be our independent variable space, or our time domain.\n",
    "We will start with generating the independent variable space using two different methods, using `np.arange` and `np.linspace`.\n",
    "They are both useful, but sometimes one is better than the other depending on the context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The step size is 0.01\n",
      "The final step is at 1.0\n"
     ]
    }
   ],
   "source": [
    "# we will start with generating an array that goes from 0-1 with a step size of 0.01\n",
    "# np.linspace, as well as np.arange, but the work in slightly different ways\n",
    "# np.linspace always includes the start and stop values, but the step size is variable\n",
    "t_linspace = np.linspace(0., 1., 101) \n",
    "\n",
    "# np.arange always includes the start value with the specified step value, but may not include the \n",
    "# wanted last step\n",
    "t_arange = np.arange(0., 1.0001, 0.01)  \n",
    "\n",
    "# make sure the outputs are correct\n",
    "print('The step size is', t_linspace[1] - t_linspace[0]) # uses np.linspace output\n",
    "print('The final step is at', t_arange[-1])   # uses np.arange output\n",
    "t = t_linspace # specify t for later use now that we have checked both"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have generated our independent variable space, we can specify what our expected response function should look like. \n",
    "In our case we would like to generate the following curve that we can pretend is expected position, *y*, given a time value, *t*, from our time domain. \n",
    "\n",
    "$y = 0.5 e^{t^2} + 1.5$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqtElEQVR4nO3dd3yV9d3/8dcnGwg7YRPClCEiQ8FR966tt1rv1gHivFtt1dZardreHXet/qxaa2sVZ521rjqqVVHcAgKC7L0JJKwkJGR/fn+cC3uaJuQAObmSnPfz8Thwje91XZ/vSXI+5/v9XsPcHRERSVxJYQcgIiLhUiIQEUlwSgQiIglOiUBEJMEpEYiIJDglAhGRBKdEIAnHzN43s8vDjkOkuVAikBbFzHLNzM0sJexYRFoLJQKRvVDC2TuL0OdIC6cfoDTIzHqZ2YtmVmBmq83smmB5FzPbYGbfCOYzzWyFmU0K5h83swfM7B0zKzazD8ysX9R+hwbrtpvZUjP776h1bczsLjNba2aFZvaxmbUBPgyK7DSzXWZ2RFD+UjNbbGY7zOytWsc52cyWBPv5I2B7qesvzOwFM3vKzIqAyWbW0cweMbM8M9toZv9nZslB+UFBvQrNbKuZPRe1Lzeza8xsVbDuzj0fmmaWZGa3BvXLN7MnzKxjsG5Pq+diM1sXbHtL1H4PN7NZZlZkZlvM7O6odRPM7FMz22lm88zsuL3Uta+ZvRT8XLcF782e9+CpqHL/1goLutZ+Y2afAKXADWY2q9a+f2hmrwbT6Wb2u6AuW4LfiTb1xSUhcHe99Kr3ReTLwmzg50AaMABYBZwarD8F2Ax0Ax4CXoja9nGgGDgGSAfuBT4O1rUD1gOXACnAaGArMDxY/yfgfaA3kAwcGewjF3AgJeo4ZwErgGHBvm4FPg3WZQUxfAtIBX4IVAGX11PfXwCVwH8FdW8DvAw8GMTcDZgJ/E9Q/lnglqBsBnB01L4cmAZ0AXKAZXuOC1waxDwAyAReAp4M1u2p40PB8UcB5cCwYP1nwMRgOhOYEEz3BrYBZwTxnBzMZ9dRz2RgHnBPUK+vYg/eg6eiyv7bex78XNYBI4L3u2PwHg+O2uZz4DvB9D3Aq8H70B54Dfht2L/bekX9PoQdgF7N+wWMB9bVWvZT4LGo+fuA+cBGoGvU8seBv0bNZwLVQF/g28BHtfb7IPC/wYfYbmBUHfH824dSsOxN4LKo+SQi31T7AZOA6VHrDNjA3hPBh1Hz3YMP4TZRy84HpgXTTwBTgD517MuB06LmrwLeDabfBa6KWncQkQSUElXHPlHrZ0Z9sH4I/BLIqnW8GwmSSdSyt4CL64jtCKAg+n2s9R40lAh+VWubp4CfB9ODg8TQNni/S4CBtY69Ouzfbb3+9VLXkDSkH9Ar6GrYaWY7gZuJfEDuMQU4GHjc3bfV2n79ngl33wVsB3oF+x1fa78XAj2IfIvPAFbuQ4z3Ru1nO5EPoN7BsaJj8Oj5ekSv70ekJZEXtf8HibQMAH4SHGummS00s0v3sq+1QTwE/6+ttS6Ff39fN0dNlxJJpACXAUOAJWb2uZmdGRXrebXe06OBnnXUsS+w1t2r6lgXi9rv4TNEEiTABcDf3b0UyCaSEGZHxfTPYLk0ExoIk4asJ/LtbXBdK4O+8ilEvhlfZWaPufuKqCJ9o8pmEuke2BTs9wN3P7mOfSYBZcBAIt0X0eq6Xe564Dfu/nQd+xpcKwaLnq9H9DHWE2kRZNX1oenum4Ergn0fDUw1sw+j3oO+wMJgOodI3Qn+7xe1qxwiXVZbgD57Dc59OXB+8D6dA7xgZl2DWJ909ysaqN+eeuWYWUod9Soh8uG9R4+6wqg1/w6QbWaHEkkIPwyWbyXSuhvh7htjiEtCoBaBNGQmUGxmNwYDuMlmdrCZHRasv5nIh8KlwJ3AE3sGUgNnmNnRZpYG/JpIN8164HVgiJlNNLPU4HWYmQ1z9xrgUeBuiwxUJ5vZEWaWTqQ7o4ZI3/oeDwA/NbMRAMHg7nnBun8AI8zsnGCw8xrq/mCrk7vnAW8Dd5lZh2CQd6CZHRsc6zwz2/PBvSN4L2qidnGDmXU2s77AtcCeweRngR+aWf8gQd4GPBfLN3Qzu8jMsoP3aWewuIZI98w3zOzU4D3LMLPjouKLNhPIA243s3ZB2aOCdXOBY8wsJxjA/mkM71Ml8DyR34EuRBIDQYwPAfeYWbcg/t5mdmpD+5Smo0Qge+Xu1cCZwKHAaiLf8B4GOprZWOBHwKSg3B1EPghvitrFM0T6/bcDY4GLgv0WExlo/g6Rb8ebg+3Tg+1+TGTc4fNg2zuApKC74TfAJ0FXwwR3fzlY/1eLnOmzADg9OM5W4DzgdiIDp4OBT/bxbZhEZKB8EZEP+xf4V3fLYcAMM9tFZED0WndfFbXtK0QG2+cSSUqPBMsfBZ4k0t+/mkgL6AcxxnMasDA45r1Exg52Bwn2LCLJuYDIt/4bqOPvPPh5fQMYRGTgdwORcRvc/R0iCevLIPbXY4zrGeAk4PlaCe1GIgPj04Ofz1QiYyLSTFiky1Sk8ZnZ48AGd7817FjCYGZO5EyaFQ0WFgmRWgQiIglOiUBEJMGpa0hEJMGpRSAikuBa3HUEWVlZnpubG3YYIiItyuzZs7e6e50X8rW4RJCbm8usWbMaLigiIl8xs7X1rVPXkIhIglMiEBFJcEoEIiIJTolARCTBKRGIiCQ4JQIRkQSnRCAikuCUCEREWoDfT13GrDXb47LvFndBmYhIolmwsZDfT11Okhnjcrs0+v7VIhARaebue2857TNSuPjI3LjsX4lARKQZW5xXxFsLt3DJUf3p2CY1LsdQIhARacbue285mekpXHZU/7gdQ4lARKSZWrq5mDfmb2bykbl0bBuf1gAoEYiINFv3vbecdmnJXHZ0/FoDoEQgItIsLdtSzD/m5zHxiFw6t0uL67GUCEREmqF7py6nbWoyVx4zIO7HUiIQEWlmFucV8Y/5eVx6dH+6xLk1AEoEIiLNzj3vLKN9egqXHx3/1gAoEYiINCvzNxTy9qItXPa1/nE9Uyha3BKBmfU1s2lmtsjMFprZtXWU6Whmr5nZvKDMJfGKR0SkJfj91GV0bJPKpXE+UyhaPFsEVcD17j4cmABcbWbDa5W5Gljk7qOA44C7zCz+HWIiIs3QnHU7eHdJPlceM4AOGU3TGoA4JgJ3z3P3OcF0MbAY6F27GNDezAzIBLYTSSAiIgnF3bnzn0vJykxjcpzuKVSfJhkjMLNcYDQwo9aqPwLDgE3AfOBad6+pY/srzWyWmc0qKCiId7giIk3ukxXb+GzVNq4+fhDt0pv2xtBxTwRmlgm8CFzn7kW1Vp8KzAV6AYcCfzSzDrX34e5T3H2cu4/Lzs6Oc8QiIk3L3bnzrSX06pjBBeNzmvz4cU0EZpZKJAk87e4v1VHkEuAlj1gBrAaGxjMmEZHm5u1FW5i3oZDrThpCekpykx8/nmcNGfAIsNjd766n2DrgxKB8d+AgYFW8YhIRaW6qa5y73l7KgOx2nDOm9jBq04hnR9RRwERgvpnNDZbdDOQAuPsDwK+Bx81sPmDAje6+NY4xiYg0Ky9/sZFlW3bxpwvGkJIczqVdcUsE7v4xkQ/3vZXZBJwSrxhERJqzsspq7n57KaP6dOSMkT1Ci0NXFouIhOSJz9awqbCMG08fSqQ3PRxKBCIiISgsreRP01Zy7JBsjhyYFWosSgQiIiH48wcrKSqr5MbTwj9RUolARKSJbdq5m8c+Wc3Zh/ZmeK//uHSqySkRiIg0sd+9vRQHfnTKkLBDAZQIRESa1IKNhbw0ZyOXHtWfPp3bhh0OoEQgItJk3J3/+8ciurRL46rjB4YdzleUCEREmsi7i/OZvmo71500uElvM90QJQIRkSZQWV3Db99czIDsdpx/eNPfWG5vlAhERJrAMzPWsbKghJ+ePozUkG4lUZ/mFY2ISCu0s7SCe6Yu46hBXTlpWLeww/kPSgQiInH2+6nLKdpdyc/OHB7qrSTqo0QgIhJHK/KLeXL6Ws4/PIehPcK/eKwuSgQiInH069cX0zYtmR+d3DwuHquLEoGISJy8t2QLHywr4NoTB9M1Mz3scOqlRCAiEgflVdX86rVFDMhux6QjcsMOZ6/i+YQyEZGE9fBHq1mzrZQnLj2ctJTm/Z27eUcnItICbdq5mz++t4JThnfnmCHZYYfTICUCEZFGdtsbi6lx52dnDg87lJgoEYiINKLPVm7j9S/z+J9jB9K3S/O4u2hDlAhERBpJRVUNP3tlAX06t+F7xzafu4s2RIPFIiKN5NFPVrMifxePXDyONmnJYYcTM7UIREQawcadu7l36nJOHt6dE4d1DzucfaJEICLSCH712kIc53+/0TIGiKMpEYiIHKBpS/J5a+EWfnDC4Gbz+Ml9oUQgInIASiuquPXvCxjULZMrvjYg7HD2iwaLRUQOwL1Tl7Nx526eu3JCs7+CuD4tM2oRkWZg0aYiHv54Nd8e15fxA7qGHc5+UyIQEdkP1TXOzS/Pp1ObVH56xtCwwzkgSgQiIvvh6Rlrmbt+Jz87czid2qaFHc4BUSIQEdlHG3fu5o43l/C1wVmcdWivsMM5YEoEIiL7wN255eX5OHDb2SOb5TOI95USgYjIPnhl7ibeX1rAj085qMXcVK4hSgQiIjHatqucX762kEP7duLiI3PDDqfRKBGIiMToF68tYld5Ff/vW4eQnNTyu4T2UCIQEYnBPxfk8dq8TVxzwmCGdG8fdjiNSolARKQBO0oquPXvCxjRqwPfPa7lPGcgVnFLBGbW18ymmdkiM1toZtfWU+44M5sblPkgXvGIiOyvX7y2kJ2lldz5rVGkJre+78/xvNdQFXC9u88xs/bAbDN7x90X7SlgZp2A+4HT3H2dmXWLYzwiIvvsrYWbeWXuJq47aTDDe3UIO5y4iFtqc/c8d58TTBcDi4HetYpdALzk7uuCcvnxikdEZF9t21XOLS/PZ1jPDlx13KCww4mbJmnjmFkuMBqYUWvVEKCzmb1vZrPNbFI9219pZrPMbFZBQUGcoxUR2XPh2AKKdldxz7dHtdg7i8Yi7jUzs0zgReA6dy+qtToFGAt8HTgV+JmZDam9D3ef4u7j3H1cdnZ2vEMWEeGVuZv458LN/PDkIQzt0Tq7hPaI6/MIzCyVSBJ42t1fqqPIBmCbu5cAJWb2ITAKWBbPuERE9mZzYRk/f2UBY/t15spjWubDZvZFPM8aMuARYLG7311PsVeAo80sxczaAuOJjCWIiISipsa54YV5VFY7d503qlVdOFafeLYIjgImAvPNbG6w7GYgB8DdH3D3xWb2T+BLoAZ42N0XxDEmEZG9evzTNXy0fCu/OftgcrPahR1Ok4hbInD3j4EGU6m73wncGa84RERitXRzMbf/cwknDevGBYfnhB1Ok2m9w+AiIvugvKqaa//6BR0yUrj93ENaxe2lY6WH14uIAL97aylLNhfz6ORxZGWmhx1Ok1KLQEQS3gfLCnjoo9VMnNCPE4Z2DzucJqdEICIJraC4nOv/NpeDurfnlq8PCzucUKhrSEQSVk2N8+Pn51FcVsXTl08gIzU57JBCoRaBiCSsRz9ZzQfLCrj168M4qEfresbAvlAiEJGENHf9Tu745xJOGd6diyb0CzucUCkRiEjCKdxdyfefmUO39hnc+a1RCXWqaF00RiAiCcXd+ckL89hcWMbz3z2Cjm1Tww4pdGoRiEhC+cuna3hr4RZuPG0oo3M6hx1Os6BEICIJ44t1O/jNG4s5cWg3Lju6f9jhNBtKBCKSEHaUVHD103Po3iGDu/57FEkJcFfRWGmMQERavZoa57rn5rJ1VwUvfO8IOrVNCzukZkUtAhFp9f44bQUfLCvg598YziF9OoUdTrOjRCAirdq0pfncM3UZZ4/uzYXjE+fW0vtCiUBEWq2120q49tkvGNqjA7edPTLhrxeojxKBiLRKuyuq+Z8nZ2NmPHjRWNqkJeZ9hGIR02CxmWUDVwC50du4+6XxCUtEZP+5Oze99CVLtxTz2OTDyOnaNuyQmrVYzxp6BfgImApUxy8cEZEDN+XDVbwydxM3nHoQxx3ULexwmr1YE0Fbd78xrpGIiDSC95fmc/s/l/D1kT256riBYYfTIsQ6RvC6mZ0R10hERA7QqoJd/CAYHL7zvMR67vCBiDURXEskGZSZWXHwKopnYCIi+6JwdyVXPDGLlCRjysSxtE3T9bKxiumdcvfEfWKDiDR7VdU1fP+ZOazbXsqTl42nbxcNDu+LmFOmmX0TOCaYfd/dX49PSCIi++b//rGYj5Zv5Y5zRzJhQNeww2lxYuoaMrPbiXQPLQpe15rZb+MZmIhILJ6cvpbHP13DFV/rz7cP05XD+yPWFsEZwKHuXgNgZn8BvgB+Gq/AREQa8v7SfH7x6kJOGNqNm04fFnY4Lda+XFncKWq6YyPHISKyT5ZsLuL7z3zBQd3bc9/5o0nWbaX3W6wtgt8CX5jZNMCIjBXcFLeoRET2Ir+ojEsf+5zM9BQenXwY7dJ1htCBiPWsoWfN7H3gsGDRje6+OW5RiYjUo6S8isv+Mouduyt5/rtH0KNjRtghtXh77Roys6HB/2OAnsCG4NUrWCYi0mQqq2u4+pk5LMor4o8XjGZEL/VSN4aGWgQ/Aq4E7qpjnQMnNHpEIiJ1cHdufXkB7y8t4PZzRnLC0O5hh9Rq7DURuPuVweTp7l4Wvc7M1B4TkSZz77vLeW7Weq45cTDfOVyniTamWM8a+jTGZSIije7pGWv5/dTlfGtsH3540uCww2l19toiMLMeQG+gjZmNJnLGEEAHQNdwi0jcvTk/j5/9fQEnDO3Gb8/RU8bioaExglOByUAf4O6o5cXAzXGKSUQEgE9XbuXav85ldE5n/nTBGFKT9VDFeGhojOAvwF/M7Fx3f7GJYhIRYf6GQq58Yja5WW155OJxetRkHDXUNXSRuz8F5JrZj2qvd/e769hMROSALN9SzKRHZ9CpbSpPXDqeTm3Twg6pVWuondUu+D8TaF/Hq15m1tfMppnZIjNbaGbX7qXsYWZWZWbf2ofYRaQVWr+9lImPzCQlOYmnLhuvC8aaQENdQw8G//9yP/ZdBVzv7nPMrD0w28zecfdF0YXMLBm4A3h7P44hIq3IlqIyLnpkBrsrq3nufyaQm9Wu4Y3kgMV6G+r/Z2YdzCzVzN41swIzu2hv27h7nrvPCaaLgcVEzkCq7QfAi0D+PsYuIq3I1l3lXPDQdLYWl/PYJYcxtEeHsENKGLEOwZ/i7kXAmcAaYBBwQ6wHMbNcYDQwo9by3sDZwJ8b2P5KM5tlZrMKCgpiPayItBA7Syu46OEZbNy5m0cnH8aYnM5hh5RQYk0Ee7qQvg487+6FsR7AzDKJfOO/Lkgm0X5P5AZ2NXvbh7tPcfdx7j4uOzs71kOLSAtQVFbJpEdnsqqghIcmjWO8njDW5GK9d+vrZrYE2A18z8yygbIGtsHMUokkgafd/aU6iowD/hpcIJIFnGFmVe7+9xjjEpEWrKiskkmPzGRxXhF/vnAsXxusL3phiPU21DeZ2f8DCt292sxKgLP2to1FPt0fARbXd5qpu/ePKv848LqSgEhiKC6r5OJHZ7JgYyH3XziGk4brJnJhiSkRBN/sLwKOCb69fwA80MBmRwETgflmNjdYdjOQA+DuDW0vIq1UcVklkx/7nPkbCvnjBWM4ZUSPsENKaLF2Df0ZSAXuD+YnBssur28Dd/+Yf92bqEHuPjnWsiLSchXu/ldL4L7zR3PawUoCYYs1ERzm7qOi5t8zs3nxCEhEWq+dpRVMejQyJnD/hWoJNBexnjVUbWYD98yY2QCgOj4hiUhrtL2kggsfnsGSvGIeuGiskkAzEmuL4AZgmpmtCuZzgUviEpGItDr5RWVc+PAM1m0vZcqksRx3ULewQ5IosbYIPgEeBGqA7cH0Z/EKSkRajw07Sjnvwc/YuHM3j19yuJJAMxRri+AJoAj4dTB/AfAkcF48ghKR1mFVwS4mPjKTorJKnrp8vK4YbqZiTQQHu/vwqPlpZrao3tIikvAWbCzk4kdnAvDsFRM4uHfHkCOS+sTaNTTHzCbsmTGz8cCs+IQkIi3djFXbOH/KdDJSk3n+u0coCTRzsbYIxgKfmtm6YD4HWGpm8wF390PiEp2ItDhvL9zMD579gj6d2/DU5ePp2bFN2CFJA2JNBKfFNQoRaRWenbmOW16ez8g+nXhs8mF0aacni7UEsd5raG28AxGRlsvd+cO7K7hn6jKOOyib+y8cQ9u0WL9nStj0kxKRA1JVXcPPXlnAszPXc+6YPtx+7khSk2MdfpTmQIlARPbbrvIqrn56Dh8sK+D7xw/i+lOGENyYUloQJQIR2S9bisq49PHPWbK5mNvOHskF43PCDkn2kxKBiOyzRZuKuOwvn1O4u5KHJ43j+KG6WrglUyIQkX3y7uIt/ODZL+jYJpXnv3sEI3rpGoGWTolARGLi7jzy8Wpue2MxI3p15OGLx9G9Q0bYYUkjUCIQkQaVV1Vz68sLeH72Bk4b0YO7vz1Kp4e2IvpJishebd1VznefnM2stTu45sTBXHfiYJKSdGZQa6JEICL1mr+hkCufnMX2kgruO3803xjVK+yQJA6UCESkTi9/sYGbXpxPVmY6L37vSN04rhVTIhCRf1NZXcNtbyzmsU/WML5/F+6/cAxdM9PDDkviSIlARL6SX1TG1c/M4fM1O5h8ZC63fH2YbheRAJQIRASAmau3c/Uzc9hVVsW93zmUsw7tHXZI0kSUCEQSXE2NM+WjVdz51lJyurTlycsOZ2iPDmGHJU1IiUAkge0sreDHz89j6uJ8vj6yJ7efO5L2GalhhyVNTIlAJEHNXruda56dS35xGb/85ggmHdFPdw5NUEoEIgmmpsZ54MOV3PX2Mnp3asML3z2SUX07hR2WhEiJQCSBbCkq4/q/zePjFVs585Ce3HbOSDqoKyjhKRGIJIh3Fm3hJy/Mo6yyht+eM5LvHNZXXUECKBGItHqlFVXc9sZinpq+juE9O/CH80czqFtm2GFJM6JEINKKzVu/kx8+N5dVW0u4/Oj+3HDaQaSnJIcdljQzSgQirVBldQ33T1vJH95bTvf26Txz+XiOHJQVdljSTCkRiLQyy7YUc/3f5jF/YyH/dWgvfnnWwXRsowFhqZ8SgUgrUVVdw8Mfr+but5eRmZHCny8cw+kje4YdlrQASgQircDSzcX85IV5zNtQyKkjuvObs0eSpTuGSoyUCERasIqqGh78YCX3vbeCzIwU7jt/NGce0lOnhco+USIQaaG+WLeDm16cz9ItxZx5SE9++c0Rem6A7Je4JQIz6ws8AXQHHJji7vfWKnMhcCNgQDHwPXefF6+YRFqDXeVV3PX2Uh7/dA3d22fw0KRxnDy8e9hhSQsWzxZBFXC9u88xs/bAbDN7x90XRZVZDRzr7jvM7HRgCjA+jjGJtFjuzlsLt/CLVxeypbiMi8b34yenHaS7hcoBi1sicPc8IC+YLjazxUBvYFFUmU+jNpkO9IlXPCIt2frtpfzytYVMXZzPsJ4d+PNFYxid0znssKSVaJIxAjPLBUYDM/ZS7DLgzXq2vxK4EiAnJ6exwxNptsoqq5ny4Sr+NG0FyUnGLWcM45KjcknR4yOlEcU9EZhZJvAicJ27F9VT5ngiieDouta7+xQi3UaMGzfO4xSqSLPy3pIt/Oq1RazZVsrXR/bk1jOH0bNjm7DDklYoronAzFKJJIGn3f2lesocAjwMnO7u2+IZj0hLsHprCb96bSHTlhYwILsdT1x6OMcMyQ47LGnF4nnWkAGPAIvd/e56yuQALwET3X1ZvGIRaQkKd1fyx/eW8/ina0hPSeaWM4Zx8ZG5pKWoG0jiK54tgqOAicB8M5sbLLsZyAFw9weAnwNdgfuDC2Cq3H1cHGMSaXaqqmt49vP13PPOMnaUVvDfY/ty/alD6NY+I+zQJEHE86yhj4lcH7C3MpcDl8crBpHmzN15b0k+t72xmJUFJRzevws/P3M4B/fuGHZokmB0ZbFICOat38ntby7hs1XbGJDVjgcnjuWU4d11awgJhRKBSBNavbWE3721lH/Mz6NruzR++c0RXDA+h1SdDiohUiIQaQJ5hbv5w7sr+Nus9aSnJHHNiYO58pgBZKbrT1DCp99CkTgqKC7ngQ9W8uT0tbg7Eyf046rjB2ogWJoVJQKRONheUsGDH67kiU/XUl5VzTlj+nDtiYPp26Vt2KGJ/AclApFGtHVXOQ99tIonP1vL7spqzhrVi2tOHMyA7MywQxOplxKBSCPILyrjoY9W8dT0dZRVVXPmIb249sRBDOrWPuzQRBqkRCByANZvL+WBD1by/OwNVFXX8F+H9uaq4wcxqJtaANJyKBGI7IdFm4qY8uFKXvsyj2Qzzh3bh+8eO4B+XduFHZrIPlMiEImRu/Ppym08+OEqPlxWQLu0ZC45MpfLvzaAHh11FpC0XEoEIg2oqKrh9S838dBHq1mcV0RWZho3nHoQF43vR8e2ejqYtHxKBCL12LarnGdmrOPJ6WvJLy5ncLdM7jh3JGcd2puM1OSwwxNpNEoEIrUs2FjI45+u4dV5m6ioquGYIdnc8a1cjh2cTVKS7gUkrY8SgQiRR0K+uSCPJz9by5x1O2mblsx5Y/twyVG5OgVUWj0lAkloa7aW8Ozn63h+1ga2l1TQP6sdPz9zOOeO7UPHNur/l8SgRCAJp7yqmrcWbuGvM9fx6cptJCcZJw/rzkUT+nHkwK7q/pGEo0QgCWNxXhHPfb6ev8/dyM7SSnp3asOPTxnCeeP60r2DTv+UxKVEIK3a9pIKXpm7kRfnbGDBxiLSkpM4eUR3vj2uL0cPytK3fxGUCKQVKqus5r0l+bz8xUbeX5pPZbUzolcHfn7mcM4e3ZvO7dLCDlGkWVEikFahusaZvmobr87dxBsL8iguq6Jb+3QmH5nLuWP7MLRHh7BDFGm2lAikxaqpceas28HrX+bxj/l5FBSX0y4tmVNH9ODsMb05cmAWyer6EWmQEoG0KDU1zux1O3hz/mbeXJBHXmEZaSlJHH9QNt8c1ZsTh3XTVb8i+0iJQJq9yuoapq/axtsLt/DWws3kF5eTlpzEMUOyuPG0oZw0vLue/StyAPTXI81ScVklHywrYOqiLby3JJ+isioyUpM4bkg3Th/ZgxOGdqN9hi74EmkMSgTSbKzZWsK7S/KZtiSfGau3UVntdG6bysnDe3DqiO58bXA2bdLU7SPS2JQIJDRlldV8tmobHywt4INlBazeWgLAoG6ZXHJUf04e3p0xOZ014CsSZ0oE0mRqapwlm4v5eEUBHy7bysw126moqiE9JYkjBnbl4iP6ccLQ7uR0bRt2qCIJRYlA4sbdWbOtlM9WbuOTlVuZvnIb20oqABjSPZOJE/pxzJBsxvfvojN9REKkRCCNxt1ZvbWEGau3M3P1dj5buY3NRWUAdO+QzrFDsjlyUBZHD8rSox1FmhElAtlvldU1LNpUxKy1O5i1Zjufr9nB1l3lAGRlpjNhQBcmDOjKhAFdGZjdDjP19Ys0R0oEErPNhWXMXb+DL9bv5Iu1O/ly407KKmsA6NulDV8bnMXh/bswvn8X+mfpg1+kpVAikDrtKKlgwaZCvtxQyLz1O/lyQ+FX3TypycaIXh254PB+jOnXiXH9uqirR6QFUyJIcO5OXmEZizYVsSiviIWbClmwsYiNO3d/VaZ/VjsO79+FQ/t24tCcTgzv2UGDuyKtiBJBAiksrWRZfjHLthSzbHMxizcXsySviKKyKgDMILdrO0bndGLSEf04uHdHDu7dUY9sFGnllAhamZoaZ3NRGasKSli1dRcr83exPH8XK/J3kV9c/lW5zPQUhnTP5MxRvRjaoz0jenXgoB4ddM8ekQSkv/oWqLK6hs2FZazbXsrabaWs3V7C2q2lrNlWwpptJV8N4ELkA39gt0yOGZLNoG6ZHNS9PYO7Z9K7UxsN5ooIoETQ7Lg7RWVVbC4sI69wN3mFZWzauZuNO3azIfg/r3A3Nf6vbdKSk+jTpQ39u7bjqEFZ5Ga1Y2B2OwZmZ9Ktfbo+8EVkr+KWCMysL/AE0B1wYIq731urjAH3AmcApcBkd58Tr5jC4u6UVlSzvaSC7SUVbCspZ2txBVtLyikojrzyi8vJLypjS1E5uyur/237JIMeHTLo3bkNh+V2pm+X3vTt3JY+XdrQr2s7enTI0P14RGS/xbNFUAVc7+5zzKw9MNvM3nH3RVFlTgcGB6/xwJ+D/0Pn7lTVOOVVNZRVVn/1Kq2oZndFNSUVVewqr6akvIpdZVUUl1dRXFZJ0e4qisoqKdxdSWFpJTtKK9i5u5KKqpo6j9MuLZluHTLIzkxnRO+OnDgsgx4dMujeMYNeHTPo2akN3dqnk5qc1MTvgIgkirglAnfPA/KC6WIzWwz0BqITwVnAE+7uwHQz62RmPYNtG9X7S/P59euRQ3vwT7U71TVOTY1TWeNUVddQVe1UVNdQUV2D+153+R8y01Po2CaV9hkpdGiTSr+ubTm0byc6tU2lc7s0urRLo0vbNLpmppGVmU7XzDTapql3TkTC1SSfQmaWC4wGZtRa1RtYHzW/IVj2b4nAzK4ErgTIycnZrxjaZ6RGHmAe9KAYkJxkJJuRlGSkJhspSUmkJBtpyUmkpSSRlpxERmoy6alJZKQkk5GWTNvUZNqmJdM2PYXM9GTapqWQmZFCZloKSeqeEZEWKO6JwMwygReB69y9aH/24e5TgCkA48aN28fv6RFj+3VmbL/O+7OpiEirFteOZzNLJZIEnnb3l+ooshHoGzXfJ1gmIiJNJG6JIDgj6BFgsbvfXU+xV4FJFjEBKIzH+ICIiNQvnl1DRwETgflmNjdYdjOQA+DuDwBvEDl1dAWR00cviWM8IiJSh3ieNfQxXw3N1lvGgavjFYOIiDRMJ6eLiCQ4JQIRkQSnRCAikuCUCEREEpz5vt5HIWRmVgCs3c/Ns4CtjRhOS6A6JwbVOTEcSJ37uXt2XStaXCI4EGY2y93HhR1HU1KdE4PqnBjiVWd1DYmIJDglAhGRBJdoiWBK2AGEQHVODKpzYohLnRNqjEBERP5TorUIRESkFiUCEZEE1yoTgZmdZmZLzWyFmd1Ux/p0M3suWD8jeIJaixZDnX9kZovM7Esze9fM+oURZ2NqqM5R5c41MzezFn+qYSx1NrP/Dn7WC83smaaOsbHF8LudY2bTzOyL4Pf7jDDibCxm9qiZ5ZvZgnrWm5n9IXg/vjSzMQd8UHdvVS8gGVgJDADSgHnA8FplrgIeCKa/AzwXdtxNUOfjgbbB9PcSoc5BufbAh8B0YFzYcTfBz3kw8AXQOZjvFnbcTVDnKcD3gunhwJqw4z7AOh8DjAEW1LP+DOBNInd3ngDMONBjtsYWweHACndf5e4VwF+Bs2qVOQv4SzD9AnBi8CCdlqrBOrv7NHcvDWanE3kaXEsWy88Z4NfAHUBZUwYXJ7HU+QrgT+6+A8Dd85s4xsYWS50d6BBMdwQ2NWF8jc7dPwS276XIWcATHjEd6GRmPQ/kmK0xEfQG1kfNbwiW1VnG3auAQqBrk0QXH7HUOdplRL5RtGQN1jloMvd19380ZWBxFMvPeQgwxMw+MbPpZnZak0UXH7HU+RfARWa2gcjDrn7QNKGFZl//3hsU94fXS/NiZhcB44Bjw44lnswsCbgbmBxyKE0thUj30HFEWn0fmtlId98ZZlBxdj7wuLvfZWZHAE+a2cHuXhN2YC1Fa2wRbAT6Rs33CZbVWcbMUog0J7c1SXTxEUudMbOTgFuAb7p7eRPFFi8N1bk9cDDwvpmtIdKX+moLHzCO5ee8AXjV3SvdfTWwjEhiaKliqfNlwN8A3P0zIIPIzdlaq5j+3vdFa0wEnwODzay/maURGQx+tVaZV4GLg+lvAe95MArTQjVYZzMbDTxIJAm09H5jaKDO7l7o7lnunuvuuUTGRb7p7rPCCbdRxPK7/XcirQHMLItIV9GqJoyxscVS53XAiQBmNoxIIiho0iib1qvApODsoQlAobvnHcgOW13XkLtXmdn3gbeInHHwqLsvNLNfAbPc/VXgESLNxxVEBmW+E17EBy7GOt8JZALPB+Pi69z9m6EFfYBirHOrEmOd3wJOMbNFQDVwg7u32NZujHW+HnjIzH5IZOB4ckv+YmdmzxJJ5lnBuMf/AqkA7v4AkXGQM4AVQClwyQEfswW/XyIi0ghaY9eQiIjsAyUCEZEEp0QgIpLglAhERBKcEoGISIJTIhBpgJl1MrOrguleZvZC2DGJNCadPirSgOA25a+7+8FhxyISD63ugjKROLgdGGhmc4HlwDB3P9jMJgP/BbQjchuH3xG5VfJEoBw4w923m9lA4E9ANpELgK5w9yVNXQmR+qhrSKRhNwEr3f1Q4IZa6w4GzgEOA34DlLr7aOAzYFJQZgrwA3cfC/wYuL8pghaJlVoEIgdmmrsXA8VmVgi8FiyfDxxiZpnAkfzr1h4A6U0fpkj9lAhEDkz0XVxrouZriPx9JQE7g9aESLOkriGRhhUTua31PnP3ImC1mZ0HXz1vdlRjBidyoJQIRBoQ3L3zk+Bh4nfuxy4uBC4zs3nAQup+pKZIaHT6qIhIglOLQEQkwSkRiIgkOCUCEZEEp0QgIpLglAhERBKcEoGISIJTIhARSXD/H3ddFB3bFC54AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# this is fairly straight forward as you can use generic python syntax for everything but the exponential\n",
    "# the only numpy specific function you might need is np.exp for the exponential\n",
    "expected_response = 0.5 * np.exp(t**2) + 1.5\n",
    "\n",
    "# plot expected curve\n",
    "fig = plt.plot(t, expected_response)\n",
    "plt.title('expected response curve')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('position');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our expected value we can simulate noise in our data. \n",
    "We will do this through introducing random values such that the average should replicate our data.\n",
    "Below we have *N* represent the normal distribution.\n",
    "The `np.random.randn` function is N(0,1), and the transformation to the correct Normal distribution space is included below since this is not a statistics course. \n",
    "We will generate 50 samples of the same length as our *t* vector.\n",
    "\n",
    "$t_{rand} = t * N_{mult}(\\mu = 1, \\sigma^2 = 0.001) + N_{add}(\\mu = 0, \\sigma^2 = 0.005)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABRmUlEQVR4nO3dd3hUVfrA8e+ZlknvhBYIvUNooqhYcAV1Ze2oWNe6Yl0V2659d3Vtq2t3de1i7z9dCyhWpEqXXgIBUiczmT73/P6YyZBMZpIJZEqS83mePE45c++5Cd53TnuPkFKiKIqidF26RFdAURRFSSwVCBRFUbo4FQgURVG6OBUIFEVRujgVCBRFUbo4Q6Ir0FYFBQWypKQk0dVQFEXpUJYsWVIppSwM916HCwQlJSUsXrw40dVQFEXpUIQQ2yK9p7qGFEVRujgVCBRFUbo4FQgURVG6uA43RhCOx+OhrKwMp9OZ6Koo7cxsNtO7d2+MRmOiq6IonVanCARlZWVkZmZSUlKCECLR1VHaiZSSqqoqysrK6NevX6KroyidVqfoGnI6neTn56sg0MkIIcjPz1ctPUWJsU4RCAAVBDop9XdVlNjrNIFAURRF2T8qELQTvV5PaWlp8Oe+++6L+Tlra2t58skn2/y5O++8kwcffDAGNVIUpT05Pb64nKdTDBYng9TUVJYvXx7XczYEgiuuuCKu51UUJT721Dnpk5cW8y5S1SKIIYvFwpAhQ/jtt98AOOuss3juuecAyMjI4LrrrmPEiBFMnTqViooKADZt2sT06dMZP348hx9+OOvWrQNgz549nHzyyYwZM4YxY8bw448/cvPNN7Np0yZKS0u58cYbAXjggQeYOHEio0eP5o477gjW5W9/+xuDBw/msMMOC9ZHUZTkVufwUu+OfasgZi0CIYQZWACkBM7zjpTyjghlTwXeASZKKQ84kdCRRx7Z7LUzzjiDK664ArvdzvHHH9/s/QsuuIALLriAyspKTjvttCbvffPNN62e0+FwUFpaGnx+yy23MHPmTB5//HEuuOACrrnmGmpqarjkkksAqK+vZ8KECTzyyCPcfffd3HXXXTz++ONceumlPP300wwaNIiFCxdyxRVXMG/ePK6++mqOOOII3n//fXw+Hzabjfvuu49Vq1YFWyJffPEFGzZs4JdffkFKyYwZM1iwYAHp6enMnTuX5cuX4/V6GTduHOPHj4/696koSmJYXR6sTg8ZKbHtvInl0V3A0VJKmxDCCHwvhPhMSvlz40JCiEzgGmBhDOsSc5G6hn73u9/x9ttvM3v2bH799dfg6zqdjpkzZwJwzjnncMopp2Cz2fjxxx85/fTTg+VcLhcA8+bN4+WXXwb84xHZ2dnU1NQ0OdcXX3zBF198wdixYwGw2Wxs2LABq9XKySefTFpaGgAzZsxovwtXFCUmNE3icPuwOr30yI7tuWIWCKSUErAFnhoDPzJM0XuA+4Eb2+vcLX2DT0tLa/H9goKCqFoA0dI0jbVr15KWlkZNTQ29e/cOW04IgaZp5OTk7PdYg5SSW265hcsuu6zJ6//617/263iKoiSOze1FSrA6vTE/V0zHCIQQeiHEcmAv8KWUcmHI++OAYinlp60c51IhxGIhxOKGvvSO4pFHHmHYsGG8/vrrXHjhhXg8HsAfIN555x0AXn/9dQ477DCysrLo168fb7/9NuC/sTe0IqZOncpTTz0FgM/nw2KxkJmZidVqDZ5r2rRpvPDCC9hs/vi7c+dO9u7dy5QpU/jggw9wOBxYrVY+/vjjuF2/oij7p97lDwBWpyfm54ppIJBS+qSUpUBv4CAhxMiG94QQOuBh4PoojvOslHKClHJCYWHYfRUSrmGMoOHn5ptv5rfffuM///kPDz30EIcffjhTpkzh3nvvBSA9PZ1ffvmFkSNHMm/ePG6//XYAXnvtNZ5//nnGjBnDiBEj+PDDDwF49NFHmT9/PqNGjWL8+PGsWbOG/Px8Dj30UEaOHMmNN97Isccey9lnn80hhxzCqFGjOO2007BarYwbN46ZM2cyZswYjjvuOCZOnJiw35OiKNGxBVoCXp+M+TRS4e/BiT0hxO2AXUr5YOB5NrCJfd1H3YFqYEZLA8YTJkyQoRvTrF27lmHDhsWk3rGSkZER/OautKwj/n0V5UAt3V5Dtc0NwOjibLplmg/oeEKIJVLKCeHei1mLQAhRKITICTxOBX4HrGt4X0ppkVIWSClLpJQlwM+0EgQURVG6ioauIYj9OEEsu4Z6APOFECuARfjHCD4RQtwthOjy01ZUa0BRlEjcXg2XRws+j3UgiOWsoRXA2DCv3x6h/JGxqouiKEpHYnM1vfHHesBYrSxWFEVJMvUhgcDl0XB7tQilD5wKBIqiKEkmXFdQLFsFKhAoiqIkgMcX+Rt+aNcQxHacQAWCdrJnzx7OPvts+vfvz/jx4znkkEN4//33E1afb775hh9//LHNnyspKaGysjLs64cffniT10pLSxk5cmSzsoqitG5HtT3ie6FdQxA+OLQXFQjagZSSk046iSlTprB582aWLFnC3LlzKSsri+l5vd7I/zD2NxC0xGq1smPHDsA/t19RlNaFWwxmd3upsLrClre7vfi05uu7HDFcVKYCQTuYN28eJpOJyy+/PPha3759ueqqqwB/Sogbb7wxmB76mWeeAfw36yOPPJLTTjuNoUOHMmvWLBoW+C1ZsoQjjjiC8ePHM23aNMrLywF/ZtVrr72WCRMm8Oijj/Lxxx8zadIkxo4dyzHHHMOePXvYunUrTz/9NI888gilpaV89913VFRUcOqppzJx4kQmTpzIDz/8AEBVVRXHHnssI0aM4OKLL6alBYZnnHEGb775JgBvvPEGZ511VvC9SNdos9mYOnUq48aNY9SoUcGV0lu3bmXYsGFccskljBgxgmOPPRaHw9Eufw9FiaW2LMKtsLpYtLW62Y293OLEFWHwt94Vn81oGut8G9Ncey209wYxpaXQQuK21atXM27cuIjvP//882RnZ7No0SJcLheHHnooxx57LADLli1j9erV9OzZk0MPPZQffviBSZMmcdVVV/Hhhx9SWFjIm2++yW233cYLL7wAgNvtpmF1dU1NDT///DNCCP7zn//wz3/+k4ceeojLL7+cjIwMbrjhBgDOPvtsrrvuOg477DC2b9/OtGnTWLt2LXfddReHHXYYt99+O59++inPP/98xOs49dRTufDCC7nhhhv4+OOPee2113jllVdavMbi4mLef/99srKyqKys5OCDDw5mP92wYQNvvPEGzz33HGeccQbvvvsu55xzTtR/FkVJBIfHR5qp9VunxeFh1U4LPk1SVmOnb3568L3dFicen4aUstmmM/HalayxzhcIksDs2bP5/vvvMZlMLFq0iC+++IIVK1YEk8xZLBY2bNiAyWTioIMOCmYkLS0tZevWreTk5LBq1Sp+97vfAf5v2z169AgevyF9NUBZWRkzZ86kvLwct9tNv379wtbpq6++Ys2aNcHndXV12Gw2FixYwHvvvQfACSecQG5ubsTrys/PJzc3l7lz5zJs2LBgWmsg4jX27t2bW2+9lQULFqDT6di5cyd79uwBoF+/fsE9HMaPH8/WrVuj+v0qSiLZ3a0HAofbx687aoMtgW1Vdopz09DpBLV2N47AZjMur4bZqG/yWZdXBYIDl4CUyyNGjODdd98NPn/iiSeorKxkwgR/Wg8pJf/+97+ZNm1ak8998803pKSkBJ/r9Xq8Xi9SSkaMGMFPP/0U9nzp6fu+WVx11VX8+c9/ZsaMGXzzzTfceeedYT+jaRo///wzZvOB5SuZOXMms2fP5sUXX2zyeqRrfPHFF6moqGDJkiUYjUZKSkpwOp0Aza5ddQ0pHYEjih3DdtTYm8z7d3s1dtY6KM5Lo9ziDL4eLhA4PbFbLxCJGiNoB0cffTROpzOYJhrAbt83I2DatGk89dRTwRTU69evp76+PuLxhgwZQkVFRTAQeDweVq9eHbasxWKhV69eALz00kvB10NTVB977LH8+9//Dj5v2PNgypQpvP766wB89tlnzTa7CXXyySczZ86cZjf8SNdosVjo1q0bRqOR+fPns23bthaPryjJzh5FIHCFuZlvrarH69PYU9c4EDQ/ViJaBCoQtAMhBB988AHffvst/fr146CDDuL888/n/vvvB+Diiy9m+PDhjBs3jpEjR3LZZZe1OOPHZDLxzjvvcNNNNzFmzBhKS0sjzgC68847Of300xk/fjwFBQXB10888UTef//94GDxY489xuLFixk9ejTDhw/n6aefBuCOO+5gwYIFjBgxgvfee48+ffq0eK2ZmZncdNNNmEymJq9HusZZs2axePFiRo0axcsvv8zQoUOj+p0qSrKyu1ufxun2hbnBezRW7LTg9e0bOA63WjgRLYK4paFuL50lDbUSPfX3VZLJjxsrmTywoOUymyqxRzH7p19hOgMKM5q8Nm/dHrQwsSA7zcjEkrw21bWxhKShVhRF6Ww0TeKMousm2rxAoV1Ibq8WNgjEmgoEiqIoUfJqEk1rOT2Epskm3T8tcYccJ5ogEwsqECiKokSpYTpoS9/4Q2/uLXGFrBkIN8gcDyoQKIrS5exvSmdvoN+mvQJBsxZBAhaTgQoEiqJ0Qfs7RTPYImjhZt+WIOP2ak1SVkRKOxFrKhAoitLl7O8N1+OLomuoDceWsmldEtUiiNnKYiGEGVgApATO846U8o6QMn8GLga8QAXwRynlAa84+mrNngM9RBPHDC9qtYwQglmzZvHqq68C/sygPXr0YNKkSXzyySd89NFHrFmzhptvvrnVYx133HE899xznHPOOWzevJlt27YF85GcdNJJfPXVV+2y57Fer2fUqFF4vV769evHK6+8Qk5OzgEfV1GSma8Ng7nhPgstB5K2dju5fftWF3fGFoELOFpKOQYoBaYLIQ4OKbMMmCClHA28A/wzhvWJqfT0dFatWhVMk/Dll18GV/wCzJgxI6og4HA4qKqqCuYfysnJCWYKra2tDWYhbQ+pqaksX76cVatWkZeXxxNPPNFux1aUZOXTZLCvv63ae4wAmg4Qhw4ex0vMAoH0a/jaagz8yJAy86WUDbkYfgZ6x6o+8XD88cfz6aefAs3TNL/44otceeWVAFxwwQVcffXVTJ48mf79+wcTtcG+1NQNzjzzTObOnQvAe++9xymnnBJ8L1KK50WLFjF69GicTif19fWMGDGCVatWtVj3Qw45hJ07dwLwyy+/cMghhzB27FgmT57Mb7/9FryGU045henTpzNo0CDmzJkT/Pzzzz/P4MGDOeigg7jkkkuC1xop/bWiJEp7tAjaa4wg9FidsUWAEEIvhFgO7AW+lFIubKH4RcBnsaxPrDXctJ1OJytWrGDSpEkRy5aXl/P999/zySefNGkpfPbZZ0yfPj34fOrUqSxYsACfz8fcuXObZB41m828//77LF26lPnz53P99dcjpWTixInMmDGDv/zlL8yZM4dzzjmnxZ3EfD4fX3/9dTA99NChQ/nuu+9YtmwZd999N7feemuw7PLly3nzzTdZuXIlb775Jjt27GDXrl3cc889/Pzzz/zwww+sW7cuWP6aa67huuuuY9GiRbz77rtcfPHFbfulKko780mJN8zGL9HwRjF9tK0384ZWgMenhd2QJh5imn1USukDSoUQOcD7QoiRUspmX02FEOcAE4Ajwh1HCHEpcCnQai6cRBo9ejRbt27ljTfe4Pjjj2+x7EknnYROp2P48OHBtMwAP/zwAw8++GDwuV6v57DDDmPu3Lk4HA5KSkqC70kpw6Z47t69O7fffjsTJ07EbDbz2GOPha2Dw+GgtLSUnTt3MmzYsGDaa4vFwvnnn8+GDRsQQgQTyYE/MGVnZwMwfPhwtm3bRmVlJUcccQR5ef7l76effjrr168HIqe/zshouqxeUeLF55P7fcONah1BWwNBoHyiBoohTrOGpJS1wHxgeuh7QohjgNuAGVLKsHu3SSmflVJOkFJOKCwsjGldD9SMGTO44YYbmnQLhdM4BXPD9LHNmzdTXFzcLKHbmWeeydVXX80ZZ5zR5PXXXnstmOJ5+fLlFBUVBVM8V1VVYbPZsFqtwddCNYwRbNu2DSllcIzgr3/9K0cddRSrVq3i448/bvL5cGmzW9KQ/nr58uUsX76cnTt3qiCgJJRPyhZXBrekoUspXFK5Bm0dI2gIHInqFoIYBgIhRGGgJYAQIhX4HbAupMxY4Bn8QWBvrOoST3/84x+54447GDVqVJs/G9ot1ODwww/nlltuaRZcWkrxfNlll3HPPfcwa9YsbrrpphbPm5aWxmOPPcZDDz2E1+ttkto6dN+BcCZOnMi3335LTU0NXq+3yd4MkdJfK0qi+LQDbxFESjMhpcTb1sHiJGgRxLJrqAfwkhBCjz/gvCWl/EQIcTewWEr5EfAAkAG8HZgeuV1KOeNATxzNdM9Y6d27N1dfffV+ffbzzz9vctNsIIQIbjnZ2KxZszjxxBMZNWoUEyZMCKZ4fvnllzEajZx99tn4fD4mT57MvHnzOProoyOee+zYsYwePZo33niDOXPmcP7553PvvfdywgkntFrvXr16ceutt3LQQQeRl5fH0KFDg91Hjz32GLNnz2b06NF4vV6mTJkSTIGtKIlwILOGPI0+5/ZqGPVNv0u7fRptTejcsLgtkS0ClYY6STTs8xt6bR1FQ7+/1+vl5JNP5o9//CMnn3xyuxy7M/x9leSxs9ZBWbWdSf3z2/zZRVursdj9Y2bj++aSm960G9fq9LBwc3WbjikEHD20G2vK6yivDd+NCyoNdZeQkpLSYYMA+DfIKS0tZeTIkfTr14+TTjop0VVSlLC0A+gaarKpTJguoP3JYSSl/1iJbBF0vj2LlYRoPNNJUZKZV5N4DnCMAMLf9Ns6UNz4WJ1+1lA8dLQuLiU66u+qtDf/YPGBrSyG8H36+5vV1OVNbIugUwQCs9lMVVWVuml0MlJKqqqqMJvNia6K0olo0r+5jLYfrYJWWwT7eTO3u3z49nO1c3voFF1DvXv3pqysjIqKikRXRWlnZrM5mHdJUdpDQz+/R9NI0emj/pxPk01mBIXrBtrfb/V1Tk/rhWKoUwQCo9FIv379El0NRVE6AC1wN2/rgHHolNP2HCOwOBIbCDpF15CiKEq0GgKAp41dMaGBoz27hhzuxA0UgwoEiqJ0MQ2J48K1CFpqJYQmqguXZmJ/A0GiqUCgKEqX0tA1FG51sbWFvvrQwdzQNBPyAHIYJZoKBIqidCkNg8Xh9iRoqa/eEyZwNL7x7096iWShAoGiKF1KsEUQEgg8Pg17C3314ZYeNN5drKN2C4EKBIqidDEN4wDhZgE5WljdG67bp/EsIRUIFEVROghfhMHi1tI8hBtI3lFtDwaIts5CSiYqECiK0qVEmj7q8mpNunpChdvestbuYcm2Gpwen2oRKIqidASNv9WHaxH4NBncH6ClzzZmc3pZvLWGGru7/SoaZ51iZbGiKEo0Go8LhM4Caujvd3o0UgzNU0+0tJmN0+NLaPbQA6VaBIqidBmN7+WRVgq7ItzQ93cPg45ABQJFUboMX6OJ/qHTRxu6hCLNHAo3RtBZqECgKEqX0Xh1cKQkcs4IA8aqRbAfhBBmIcQvQohfhRCrhRB3hSmTIoR4UwixUQixUAhREqv6KIrS+Xl9rUwBbdwiaJY7yB8AIrUIOmr6iGjEskXgAo6WUo4BSoHpQoiDQ8pcBNRIKQcCjwD3x7A+iqJ0ck6vRl0LaSKazBryySabWXmCg8VqjKDdSD9b4Kkx8BP6m/wD8FLg8TvAVCGEiFWdFEXp3FweH7VRBgLY1yrw+LTgQLIaI2hnQgi9EGI5sBf4Ukq5MKRIL2AHgJTSC1iA/DDHuVQIsVgIsVjtQqYoSiQur9Zi4jhfSFa4hsDQeDGYzxc+i+j+bG3ZUcQ0EEgpfVLKUqA3cJAQYuR+HudZKeUEKeWEwsLCdq2joiidh9Pjw+r0RLxph6aSbrjhh64KDu0e8nbgzKLRiMusISllLTAfmB7y1k6gGEAIYQCygap41ElRlM7H5fV38UTaAzhiiyCkBRDaPdSZu4UgtrOGCoUQOYHHqcDvgHUhxT4Czg88Pg2YJ2VnjruKosRSw+bxtfYIgSDCGEFojqHQ5515oBhim2KiB/CSEEKPP+C8JaX8RAhxN7BYSvkR8DzwihBiI1ANnBnD+iiK0sk1dOlEGjDWQr5nNiwq6+otgpgFAinlCmBsmNdvb/TYCZweqzooitJ5aJpEp2t5UmFDiyDSgHHoauKGRWWtjRF09haBWlmsKEqHUOf00FLPsaZJPIEbuserUe/yNi8TbYvAHdoi6LyLyUAFAkVROgiHxxf8xh9O6HvhWgWhXTz7xgia3vid3q41RqACgaIoHYLTozX7pt70/abvhRswDr2hR5o15AnsTdAg3Eb3nYkKBIqidAhOj6/FPYVDWwS1juYbxYR2De3bZrJ5S6PxuTr7YLEKBIqidAitB4Km79ldPrwhN/jQb/Y+TTZJL9HYbouzUTk1RqAoipJwrXcNNb9Zh/b1Nxss1rSIew1vrawPBgPVIlAURUkCTq8v4n7C0LxFAGFSRYQOFvtki5vOrym3UGt3J8UYQf47r8P27TE5tgoEiqIkPY9Pw+eTONzRzxqC5tNAQ3MQ+TTZbKC4aXn4tcwSdipqPKVvWEu/m66Bhx+OyfHV5vWKoiS9hpu8y+uLuLAs3D4Coa2E0FlDHq3lFgH4ZxB5WikTU1Iy+P478GVmYfjrX2NyCtUiUBQl6TXc5KX0dxGFkjL8Db3xuEG4tQA+TWtxbUIyKPj2S/J/WsCua+ZAfrMs/e1CBQJFUZJe42/74QaFXd7waaKbTgENt8dA5B3JkoHweBj0wF3U9xvI3nP+GLPzqK4hRVGSXuObf7gppKHZQvd9bl/ZSDNAE93/35Lec18kfesmlj35CtJojNl5VCBQFCXpNb6hh5tCGmk2kcujBccUQvciaGBPphaBlPR58SlylywkbfMG0nZspWrykVRNOQZdnQXIi8lpVdeQoihJr/GNPvygcOR+/oYxhdDdyRpEej0RUndsY/CDd5O+YS22wcPYcuk1rP7Hv/lo7ouceexktm7dGpPzqhaBoihJr7WuoZb6+Z0ejTRT893JklHWquUArPzX81iHjQJg6U8LeOIff2HyEcdQXFwck/OqFoGiKEnvgFoEgfIdIZV01spl+FLM2AYOBWDnti3c8+dL6NN/EPc8+gx6vT4m542qRSCEKAQuAUoaf0ZKGbthbEVRFBrWDjR63qjfv3GZSBpaEB0gDpC1ajnWYSODA8N5Bd04dOpxzLr8OtIzMmN23mi7hj4EvgO+ApJoZEVRlM4u3HRRh8dHeoqhxTL73guMESR515Dweslau4Kdp87C43bj9XhITU/nhnv/FfNzRxsI0qSUN7XlwEKIYuBloAiQwLNSykdDymQDrwJ9AnV5UEr537acR1GUzi100xhoHghaahE0BIlkGhQOJ33TevQOB5aRpTx69xw2rl3FY69/ginFHPNzRztG8IkQ4vg2HtsLXC+lHA4cDMwWQgwPKTMbWCOlHAMcCTwkhDC18TyKonRiYbOKNgoOVTZXi90+rg7SImgYKJ67dhX/e38uhxx1bFyCAEQfCK7BHwycQghr4KeupQ9IKcullEsDj63AWqBXaDEgUwghgAygGn8AURRFAcKnlGgIBG6vxpryFm9F+6aPJnkq6ayVy3CkpvKPl57mqONP5rzZN8bt3FF1DUkpD2iUQghRAowFFoa89TjwEbALyARmSimbxXYhxKXApQB9+vQ5kKooitLBhJsl1JCFdG15XcRVxQ0a0kgkeyAw/vI9PzidjJxwMDfc+wj+78fxEfX0USHEDCHEg4Gf37fhcxnAu8C1UsrQ0D0NWA70BEqBx4UQWaHHkFI+K6WcIKWcUFhYGO2pFUXpBCINFu+otlNhdUV1DJdHS+pAoHM6KCzbTnnf/tz12Itx6xIKnj+aQkKI+/B3D60J/FwjhPhHFJ8z4g8Cr0kp3wtT5ELgPem3EdgCDI228oqidH7hWgT1Li8b9lqjPobD42u2O1myqKutJn3NSnQ+H6P//Bcys3PiXodoZw0dD5Q2dNsIIV4ClgG3RPpAoN//eWCtlDLSbgrbganAd0KIImAIsDnKOimK0slpEfYLaOu3e6fHl9jtJjWN3F9+oOh/H5O6Yyvm8p0Y62pZfck1XPLua9yUnsEkoG5kaUKq15YUEzn4B3MBsqMofyhwLrBSCLE88Nqt+KeKIqV8GrgHeFEIsRIQwE1Syso21ElRlE6svfYKcHoTM0Yg3G76/vdJer37Gqk7d+DNyKS+/yBsQ4ZjLN/F2Ptv53idnsMPOhRXYRGuoh5xryNEHwj+ASwTQszHf8OeAtzc0geklN8HyrZUZhdwbJR1UBSli6l1uNvlOA53YgLB4H/eQfEb/6XqkClsvOZWKo45Di3FjM/r5W9XXcCNwBOaD++vi6k+ZErc69cg2llDbwghvgEmBl66SUq5O2a1UhRFgagHg1vT0srjWCn69H2K3/gv2y64nA033tnkvX/fewvfLfiK8TfdzcG//EDh/P8lrFsIWgkEQoihUsp1QohxgZfKAv/tKYTo2bBOQFEUpb1pmqSqvn1aBE6vD5M+fjk20zf+xrA7rqdm3CQ2Xntbs/ennngqPfuUcMJ5l7LizAvo+f4b7DnupLjVL1RrLYI/45+//1CY9yRwdLvXSFEUBaiqd7dbWgifT+KS8UmTprPXM+q6i/GlpbHyoWea7Cy2Yc0KBg0fzajxBzNq/MEASJOJnTPPj0vdImkxREopLw08PE5KeVTjH/wziRRFUWKivbqFGsQr+2iv994gY/MGVt//JO5u3YOvv/3iU1xx+rEs+m5efCrSBtG2lX6M8jVFUZQDJqWk0ta+gSAupKTXWy9jGTW2yeDvh6+/wLMP3MUR02cwLoGDwpG0NkbQHX9+oFQhxFj2zQLKAtJiXDdFUZKc0+PDbGz/zVLqHN6w6weSXc7ShWRsWs+ae/Ytnfr8vdd5/G+3Mvno6dx83xPoDcm3MWRrNZoGXAD0BhovCrPiXxOgKEoXZnF4YhIIKmzOdj9mPPR66xW8GZnsnv4HAHZs2cgjd9zAhEOP4raHnsHQaLwgmbQYCKSULwEvCSFOlVK+G6c6KYrSQVgcHoqyDjwvToXVhd3tpXduGnqdYG87jw/Eg7Gmim5ffMKu02ahpaUDUNxvIH956FkmHn40JlNKgmsYWWtdQ+dIKV8FSoQQfw59v4XUEYqidAG1dk+bytc5PWSZm38r3lPnZLfFyfZqOz2yU7G7Ot5GiD0+fAu920XZGefyzWcfUNCtByPHT+LwY6PO0ZkwrQ0Wpwf+m4E/TXToj6IoXZjN5WnTit29da6wu4nV2P3rBVweja2V9e1Wv7iRkl5vv0pt6UQ+Xr+Of8y5grnP/zvRtYpaa11DzwT+e1d8qqMoSkfh8Wlomj+zZ0ZKdAOgTo+PKpubnjmpwddsLm+rewokFSmh0V4Bwu2i/1MPk751E3NPn8z9N89m1PiDue2BZxJYybaJNg31P4UQWUIIoxDiayFEhRDinFhXTlGU5NWQEM7hjr4bpyEQNFZta5/Vw/FQ+NX/cdS4EsZccQ7dP36H3J+/Y9Kpx9Dv2UdZOnoc57/9CmMOOpR7n3qV1PT01g+YJKKdx3SslHKOEOJkYCtwCrAA/8bziqJ0QQ3TO8PtFxCJw+PD5vIipQzuwFVV33EGhos++wAtJYXM39ZS+O1XADh69GLZk69ywyfvMjozmzsfe4EUc2orR0ou0QaChnInAG9LKS3x3EZNUZTk09DXH20gkNK/t4CU/kHm3HQTmibbPOCcKMLrJf/HBVQcczxr7n6Y7OWLSd+0nq1HHoupsBtzJh+BpvnivrtYe4h2ZfEnQoh1wHjgayFEIdAxJ/oqitIuGloEjigDgdPjDwJAcNVwraNtg82JlLVqOca6WioPOwp0OmrHTuTvu3dy+YWnUFdbg8Fo7JBBAKIMBFLKm4HJwAQppQeoB/4Qy4opipLcXMGuoegGehu3HCoCgaC6nbKLxkP+D/OROh3VB09B0zSevO+vvPb0I4wYdxDpmc22Wu9QouoaCuw9fA4wJdAl9C3wdAzrpShKkmtzi6DRtFG7y4fD7etYgeD7+VhGj8OZnsGDt1zF15+8y6nnXcalN96BThe/FNexEG3tn8LfLfRk4Gdc4DVFUbqohjECj1eLqnsntOVQbnFgdXaM8QFjbTVZK5dRdehRvPDo3/n6k3e58JpbuGzOnR0+CED0g8UTpZRjGj2fJ4T4NRYVUhSlY2g89z+atQSh00y3VduDYwbJLu/HBQgpqTrsKE7vVUz/ISM45sTTEl2tdhNtKPMJIQY0PBFC9AdabA8KIYqFEPOFEGuEEKuFENdEKHekEGJ5oMy30VddUZREcvkaBYIo1hI4Q1YUt9emM/GQ/sUn1KWkUD14OLn5hZ0qCED0LYIbgflCiM2B5yXAha18xgtcL6VcKoTIBJYIIb6UUq5pKCCEyMHf1TRdSrldCNGtTbVXFCUuQtNN+zTZ5EYezRTStqw3SAhNQ+9w4AtZCLZl3WpGff1//E+nw1K2lZKBQxNUwdiJtkXwA/AMoAHVgcc/tfQBKWV5w57GUkorsBb/3gaNnQ28J6XcHii3N/qqK4oSL6G7hYXmC4rmJp/saSSG/P02jpg8lJE3/onsZYtASpZ+N4/3Z51AkaZR+KfrO2UQgOhbBC8DdcA9gednA68Ap0fzYSFECTAWWBjy1mDAKIT4Bn8Su0ellC+H+fyl+PdOpk+fPlFWWVGU9rLX6qQ4b99eVKGbxrQ2c8gd5YByomT9uoTiN/6LZfQ48r/7mu7/9z629AyOrLcxB9B0enSnnp3oasZMtIFgpJRyeKPn84UQayKWbkQIkQG8C1wrpawLc/7xwFQgFfhJCPGzlHJ940JSymeBZwEmTJiQvP+aFKUTcns1au0eNE2i0/kzCrhCA0ErYwSh4wPJRHi9DLtrDs6iHiz9z1tIIejxybuw4GsWbtnImLMuxF06AXdhUaKrGjPRBoKlQoiDpZQ/AwghJgGLW/tQYP3Bu8BrUsr3whQpA6qklPVAvRBiATAGWB+mrKIoCeDw+JASrC4v2an+vQRCu3mcIYGhwuqiMHPfRizJND6QtWIpOcsWUf6H0/Hk5FH82vNk/raaX//1PHa9nu++/JRjzjgPzjiP3kBVoiscB9EGgvHAj0KI7YHnfYDfhBArASmlHB36AeFfefY8sLaFDWw+BB4XQhgAEzAJeKQtF6AoSmw1fNu3NQoEbl/TG7vHq+H1aRj0OmrtblbtsnDEoMJgC8LpTp7xgUEP3U3u4p8Z8O/72PWHmfT46G0qp0zltzETuOPCU/lt5TL6Dx5O/yHDWz9YJxFtIJi+H8c+FDgXWCmEWB547Vb8QQQp5dNSyrVCiM+BFfgHov8jpVy1H+dSFCVG7G4vQGDxlz+rZri0Eg6PjwydYG25FZ9PUlnvolumP/dOsnQNGSy1ZC9bRPmM05BCR693XkXqDXx+xnn8+azjsFpquOPRF7pUEIAoA4GUcltbDyyl/B5oNUWplPIB4IG2Hl9RlPiwB1oEVqc3+Jrb1zwQOD0a1fV26l3+cnvrGgWCJOkayv/xG3Q+H2VnnI9l7EQ2XXMLy/73MXNuuJysnDweeeUjBg4blehqxl3HXxutKEpMNdzEbU7/PgIQfiqoxeFmc6NtJittLrTATKFoE9PFWsGCr3Dn5GEZPQ4AV1EPKroVMXDYKJ548/MuGQRABQJFUVrR0CLwaTL4OFyLYFuVvckiM69PUh3YizgpWgQ+H/nfzaPqsKOw1dv45buvAThi+h94+OUPyC0oTHAFEyfaMQJFUbogr09rsmbA6vSSatTj8TYPBOHyBu2tc5GXZmq27iARslYtx1RTzdohw7nyzOOo2rubV79cRHZufqdIHHcguvbVK4rSInvIN3mr0xO2NRBJhc3V7BiJUvDtV2hCcM4TD2K3WfnHM2+QnZuf6GolBdUiUBQlImfIQjGry9umVBEer8Zui6O9q9VmUkp4+2V+lJK8YaP4y8PPUtCte6KrlTRUIFAUJSJ7aCBwenH52vYNv6wmtoFA53SAEGgtbBOZUrGH/tVV/Fg6kQf/+y4GozGmdYqFWO4SrwKBoigRhQYCj1ejztG2zWS8MUw3nbZ1E2MvmYnRUsveY45n9+9PpXb8JDRTCgjBwm+/JEunZ9rmDQAMvP1+bB0wCAB0z47dfsgqEChKFxeaDqIxh8fb7LVKW3JsL5m5ZgVjLzsLgL1Tj6Pb15/R88O3AJA6HU69nsM8Hhpun44evbANHpag2h6YNJOeXjmpMTu+CgSK0oXZ3V5+222lIMNEYD/ykPebdwPZnM2DQ7zl/vIDY648H09WNsv+8xb2kgGsczooWPA1vhVLWPjp+9j37qbv4OEMP24GWmF3LGPGQ5hr7Aj6F2aE/fu0FxUIFKULq7S6cXp87K5z0iO76TdOTZNJuYeAsaaK0ddchLOoB8ueexNX954AaOZUlg4expW3/xkpJX9++FlGT5vBzgTX90BlmA0UZYVvsbUXFQgUpQtwuH3odQKToemM8cp6/4YzWyvtzQJBskz7DDXgX/9Ab7ex8pEPg0FASokQgl59+3PizPM5/vRz6NG7b4Jr2j76F6bHtDUAah2BonRq5RYHS7ZV88PGSrZV1Td5z6dJagMrf+td3ma7kDUkm0ukjHWrMe3dHXyetXIZvd59jR3nXEz9wCEArFuxlCtnTmfPrh0IIbjouts6RRBIMerolZsazNcUSyoQKEoH5PZq2Fwt36hdXh+rd9ZRU++f5bOz1tFkl7Dqejdao56f7dVNA0WiU0en7CrjoDOP45AZUyj67APQNIb87Vbc+YVs/tP1eD0eXn7iAa4550Rqqyupre4cOwfkphs5dGABhw8qZFiPrLicUwUCRemA6l1edtW2PD+/ur7p7B6vT1LeaHFXVX3TFkBNvYeyGjt1Tg8urw97mBlD8dT/af/WJPa+/Rl1w+VMPPsEslcuY8MNt7Nl726umfV7XnnyIY4+4WSeeW8eQ0aWJrS+7SU71USqSR/Xc6oxAkXpgGyBQDCwMCO4+UuoqjDTPHdUO+idmxbx/XXl1vat6H5K3baZHh/MpezMC9gw5y76PfMI/Z5+hNqxB7H796fy3l1zKC/bzl8ffo4p005MdHXbVUZK/G/LKhAoSgdU7/bi9Un2Wl1hFxpJKamqb36jr3d5qa53k2LQtbrPcCINePwBpMnE1kuuQRoMbJ59I8vHTsKWlkpPIbjk+r9y7uwbyO8A+wjnpBmptUe/CC8tJb6tAVBdQ4rSIdW7/DfxnRG6h+qc3rAZQgF2VNvDtgaSRfr6tRR99gHbZ12Mu7AbPp+Pd196hrOvOp8HH/6bv0xmVocIAkLAyF7ZGPTRz/pJN6kWgaIoUWjYBaym3o3d7SUt5OZRZXOF+xjg3zAm3EKxpKBpDPzX3/GlZ7DtwivYtmk9D/31Otb+uoRJRxzDNbf/M9E1bJPcdBNmo54e2ansqLa3Wj7FqEMfoasvlmLWIhBCFAsh5gsh1gghVgshrmmh7EQhhFcIcVqs6qMonYUnZI+AcIPG4bqFGki5L5AkE+F2MfKm2RR++yVbLruWZZt+40+nHsPObVu4+f4nuOeJVygMrBvoKIqy/N12vXKjSw+RnoDxAYhti8ALXC+lXCqEyASWCCG+lFKuaVxICKEH7ge+iGFdFKXTCL2J76p10q8gI/hN0uNre2K4uNE0CLMJjLG2mtFX/5HcJT+zevYNlF94BUM9Hk4+9xJOu+BycvM73u5hOh10C+RwykgxkJ1mxNLKWEEiuoUghi0CKWW5lHJp4LEVWAv0ClP0KuBdYG+s6qIonUno+gG3V2PZ9ho8gQ1jquvdYXcLSzTzrh0cfuRoBv/9NoR33zWkb1zHhLN/T9aKJdw36TCOevd16m1WjCYTl1z/1w4ZBABy00wY9ftusY2Txhn0gpKC9GafSYvztNEGcRksFkKUAGOBhSGv9wJOBp5q5fOXCiEWCyEWV1RUxKyeitIRhOvfr7V7WLKtBpfXR2UL4wOJ1O/pf2GqqabPa89TevnZGCy1dP/4HQ468zi0qgp+b07ltsU/ceRxf0CnT8wNsa0yzZG/wYfO5irKMmPQC1KMOiaU5FGSn9YsB16iAkHM2yFCiAz83/ivlVLWhbz9L+AmKaXWUi4NKeWzwLMAEyZMSMLvOooSP5FWFNucXpZsrcGrJd//Iqnbt9Ljg7nsOOtCrMNGMuzOOUw+YTKmmmqWZGTye5uV3DHjefL2fzJg6IhEVzcq3bPNDC7K5MdNlc32XNDpoDCjaaI4vU7QvyCDblkpmI3+G35WatPuos44RoAQwog/CLwmpXwvTJEJwNxAECgAjhdCeKWUH8SyXorSkbU00Juss4H6PfMI0mBk68VX4e7WHVuvvoy643p2nXQmN+0q46zJRzD9lLM6zCbyZqOeId0zMep19CtIZ8MeW5P389NTMOibX0uf/LQmz3PTTMFAoNeLYICIt5gFAuG/uz8PrJVSPhyujJSyX6PyLwKfqCCgKJF5fVpSpoZuSdrWTfT46G22n3MJzoJu/O/d13nt6Yd5+KUP6NazNzcnuoIRDOmeyZ46Z9jFYMN7ZgX7/4tz09hZ4wgGYbNRz8BuGVGdIy/dxNZKf46ntAQFAYhti+BQ4FxgpRBieeC1W4E+AFLKp2N4bkXplBoWkiUzU2UFfV58Ck92Lrahw+n53htophT+N3EyD545nQ2rVzBi7EG4XM5EVzUivV7QOzeV4rw0qmwudtU68WgaPk2Sl24iL90ULKvTCQZ2y2BFmYUMs4HS4pyov9nnpBrR6fyTqRLVLQQxDARSyu9pw37LUsoLYlUXRUl2dU4PJr2u1RuILQlSQ7ckY+0qSq88H1PFbnSNNrmf228gl111PgVFPbj5/ic4+oRTYp5j/0BkmY3B+uVnpJCf0fLGMN2yzJQUeOibn95kplBrdDpBdqqJmnp35wwEiqJEb3uVHZvLy/i+uS3eSOxJuBCsQeFX/8fIm2fjyc5l0ZufU1vQjfxtm0ndsZWPVv7KecefzOkX/glzalrrB0uw7NS2b3A/sFvmfp0rLz0QCBI0YwhUriFFSTivT6PC6sLm9LKirBathVk/re1BkCj5333NmGv+iG3QMH54+QNeWvQTZ5x0FD+ZUig/+Sz+ePt9nHvF9R0iCMD+BYL9lZfm72ZKUy0CRem69lhdwQ1jauo9rNplYXiPrLCzTpJxjMBQZ2HYHTdgGzCYf542ixcuPJXdO3cw7pAppKU1XzTVEcQzEGSlGjDoRacdLFYUJQrlIbmC9ta5qK6vpE9eGsV5aRj1OpweH3UOD84k3Ed40AN3klKxhwtz85l7+/UMGDqSvz/zTyYcemRSjwNEkmrSN9vbOZaEEPTMSY24r0Q8qECgKAlU7/KGnZ7o9Uk2V9SzvdqOECJiSulEq3/pGXq99wZbLr6Kol7F3JKe6V8Z3EHWA4QTz9ZAg95RJqWLFRUIFCWBGm8dGY5/xWoSrRSWEnN5GXu++ozVb73MhVs2srd7TzbPvoHfm1qeWZNIQkBOmn9QtjGzUY9BL7A59429JCIQhKYRjzcVCBQljiqsLuqcHnpmp2I26ii3JO9c+iCfj25f/R/5P8wn65svyKyqBOBUwJ5iZtkDTyGTOAiA/+Y+rk8O6/fYgvsC5KQZGdU7G7vLx5JtNcGyWQkIBImmAoGixNGuWgcVVhdbKurJMBs6xCrhgf/6OyUvPIEnM5v/ed3MT0un54mnMebci6FPf0iCBHE5aUZy003kpZlYU17XbBvOgowUhBAM6Z5JhtlAncPDkKJMdDpBikFPQWYKlVYXOh1kJnD2TqJ0vStWlATxaZLqRl0TjbsjkpGUkvKXnuHoF55g+x/OYMPdD7Nr03qOLO5LahLNBspKNTKhJC/4vCgrha2VTXcDy8/YtxK4V05qk5TQAAO7ZVBlc5FlNiZ00DZRVCBQlDipsu2bJpps0jespe8LT5K+ZSPrr7mFj+ssfPjco7yzdhU79HrmnXEexQYD/YcMT3RVm2mc7gH86Z4bB4IUo45Mc8vdPRkpBrpnm9u0KrgzUYFAUQ6A26uFnWpY5/SQFXLz2WtNrn0CDJZa8n5eQI8P36Lw26/wpabizMph/MVnMB+4PiOTQcDCp1+nuHRCoqsbUUFG00CQaTaSnmIIZmnNT49u/GJAYQbWJG+lxYoKBIpyAFbutDCuT06z+fLrd1vpm59OYWCrQk2TSbNhTN6P39L/8X+SvXIZQtNw5uSx4MTTkDffjdeYgnPmdOZs2QA2K9vPuQTr5CMSXeWI9HoRdpZPUVYKmyv8N/WCTFOz98MxG/UJSwOdaF2zHaQo7aDK5qKm3k1NyDoAp8dHrd3Dhr1WZGDPyGq7u9nmJYmQ9/18Sq84F2NtDctnns9NU4+jwGln2lefYjGakOnppHzyHcufeJmyM85j47W3JLrKLcpLM4VdtNawO5hOty+FgxKZahEoyn7aUeNfA1BucTTpp94dmBJqd/koq3FQnJfG3rrEtwZyf/6OMVdfSHXvPpzdoxdfvvFfjEYTR//+FE49/zLS0vfl0K888lgqjzw2gbWNTuj4QIM0k4GsVCMGvQibqkNpSgUCRdkPDrePqkBXz16ri2GaDM422V23b23A5sp6irLMCe0W0tfbyP3wbUY+eBf2vv34/K/3s/jaizhv9g38/ozzyS1I/ObwOWlGLA4Pso2NpvyMyN/2u2eZI76nNKUCgaK0wuvT+LWsliHds8gIzDEvq7EHb1o+n6TC5qIoy4zN5W0yLdTj1VhRVos7Tiki9DYrOcsWYbTUYLTUYPj+G3r/MJ8Un4/fsrLZ/Z+36J5fyOtfL8VgTJ6FU92zzZiN+mBrqkGqSY/bp+EL062WatK3uCK3W1YKWlsjSxelAoHSYTg9PjQp474cf91uKzX1HhZvrWZscS4ZZgM7QxLFlVucFGWZm93IgLC5hGLBUFvDxHNOJH3LxuBrlcALOh1rDj2KQVdcz/B8/7f/ZAoC4O/iyUkzNfv9jeiZhU4nWL69eTCN1C3UoKsO/O4PFQiUDqPh5jugMLr9YNtDucURvDl5fZKl22sozExpNvBbXe/C7dXYUxf7lBHC7WLovbcivB7Wz7kLb06u/7XLziK1bDsrHn6W//70HV98P4/Jp5/HtNNn8Ye8gpjXa3+lGHXB4F6YmUJFYJptjxwzOYGB3okleSzbUYO9URru/FYCgRK9WG5eXwy8DBThz5r1rJTy0ZAys4Cb8G9paQX+JKX8NVZ1UpKLx6dh0ImoUhVrmmRnjQOdEDELBD5NUlZjJz3FQE6qEbdPY91ua7My4b71axqs32NtltqgvenrbYy++kLyf/4OzWAg78dvefXE0+j+4Zt0r6rkgwv/RMa0GUydcgzH3X5/h8gCmttoVk+/wnQqrC4MesGgRjt+pZr0TCzJY0+dk5p6D7UON7kqELSbWLYIvMD1UsqlQohMYIkQ4ksp5ZpGZbYAR0gpa4QQxwHPApNiWCcliey2OBECeue2vmvVXqsr2DVQXe9utVugrZweH7/uqA0uKBIC9DoRtm86knABoj2ZKisovWIWGetWs/SOf/LSD99w+defc/ELTwDw/iFH4L3gCoAOsxMY0OSGnmU2kpdhojAjpdlCPaNeR+/cNHrnxruGnV8sN68vB8oDj61CiLVAL2BNozI/NvrIz0DvWNVHST5V9W4sDg9FWa0v7S+r2ZcyYFdt0+maHp+GxeGhIGSD8Qqri992W/FqGpqUpJsMjOiVHRzwbVBrd7OizNKkD1pKEj7vX2evZ9BD95C75GdS9uzGWFeLx2Ti13+/SOXhU/nglefYdMJJ3K5pFPQsJvPaW/0RrIMJnec/tHsmqap/P67iMkYghCgBxgILWyh2EfBZPOqjJJ6mSWrsbnw+yZbKegYXRd742+r0NBlwrbC68Pq04PzwVTstVNncFOelMahbBjqdYGtlPZsqbE2mI1qdXn7ZUsXAwkx65aZSYXWxs9bRLEd9MkjbuonR115E+sbfWDNoGCvNKayogwWp6fz1sKPR63Q8+8E36PV66oC6ONbNoBf7FSSNBh0DCtNZv8eKFoi5ZqOe1JBN2xOdm78rivlvXAiRAbwLXCulDPvvVQhxFP5AcFiE9y8FLgXo06dPjGqqxFOtwxPsdimrsdM7NzXiDWBHddMZOj5NstfqomdOKlsq66myuQPl7NTa3aSnGCJ20zT05W+s2HczSiY6ez1FX37KkL/fhtPn4wSdns/Wr6F3yQCmXnk+l/3+VPSBtM/6BKR/1usFg4oyWburbaEn02xgTHEOZqMep8cXTAqXm55cs5e6qpgGAiGEEX8QeE1K+V6EMqOB/wDHSSmrwpWRUj6Lf/yACRMmqInBnUBVowVW/puzjdLinGblnB5f2Jk45RYHKQYdmytsTV63Or1RJQ5LmiAgJWlbN5H1wzeYPnqbkrUrMWkalhFjeOeiKzEv/YV/n3AyQ0aNTYr9f3vnpFKYkcI6QbPFXzpd+N9r92wzw3tkBRfc9SvIoNzixOXRmgwUK4kTy1lDAngeWCulfDhCmT7Ae8C5Usr1saqLknwqbU27YyqtLnbVOugZkid+3W5r2NTNNfUerE5Lm1eiJgOdw07+9/Mp+PJTsr6bR2ZdLQCbgedMKaRdNJs+l/+ZvgYDV0w7MaF1bUyng+K8NEwGHRkphmYBd2TPbHZZnFQ2yrLaPdvMiJ5ZTYKYXicYXJTJyjJLuw/6K/snli2CQ4FzgZVCiOWB124F+gBIKZ8GbgfygScD/1C8UsrkzXertAunxxdMEdzYb3us5KQZg11E5RZHk5tKqEQP5raV3lrHwLvm0H3e5xhdTtzZubxfb+V7cyr1hx/NoFPOZuwhUzAYjcm0S3FQUZY5uEgrL93UJBAY9IKCjBQKM1NYU15Hea0zbBBofKzqXLda9JUkYjlr6Hv86wNaKnMxcHGs6qAkp6oIg7M+n2RlmYWJJXl4NI3fQubwdxTmXTsYcdOVVB12FNsuvILqOgurP3qbWU8/TPd6G29l5zD0yVeonXAIlk3rOaH/oKRb6RtOSf6+Xcny0k1sq9o3k6sgIyXY9TOiZzb56SkUZaW02J01tHvkCQJKfKnheSXuWvqWb3V62VRhw+72dbhv/ODP9TPminNJ37KJ3KULMb/8DLfW1nAPUCAE/zh6OhnnX0a38QcjhEjKHb/CKchMIb3RtNucNFOTMYFuWU2n7jakgW5JMox5KH4qEChxpWmSanvL0zUbf9NMelKCEDgddlb8tIBpf7sN8+6dzLv/CfKzsun3l2t5GbDn5LLimTeYPLI00TVuQq8XDCjIIMWow6ATVNrc7Khu+vsXAvoXNt2jWK/zbwhTU+9BrxcURLkLmJKcVCBQ4qrS5mrTat1Ey1i7it5vv0z2r0vZNPtGKo+eFnyvx4dvMfBvt7JDSr52OuivaYwH7hsygsKhI2HgEGr+9ws9PnyLqsOPxtkruaY+G/SCscW5ZKft65bKSTNRaXM1SZXRKze12babAHnpKdTUeyhs1C2kdEwqEChxo2mSDXttrRdMAjmLfmTQw/eSvWIpvhQzrm5FjLn6Ar6ePoN/paZxXU01I+b/j8oRYyjbvIFz9XrSNI1N51/GhDl3BY+jpaax88wLEnchjYhGUz4NesG4vrnNbvANM3p+3VELgMmgi5jbKS/NxCagW6ZqDXR0KhAoB8Tq9JAZ5ttiOFur6mOelO1AGeosDHzoHnq/8yqOXsX8dtNdPLBrJwsX/cjNUjLrsw8ZIAT9pGTnqbNY99f7kEYjP3m9mHfuwNGnJNGXEKTT+fM45aQayU4zYtLrcHh82N0+Uo36Jn3+jRVmpgSzgA4uyoyY/iMr1UCKUUd+hgoEHZ0KBMp+21JZz6a9NkYXZ9Mts+XBQafHlxR9/8Ljoed7r5Oz9Bcqj/wdFUdPR0sxk7KnnIL33qD4ledIq6vlg2GjyHrlQ7TUNH45/yRERgZfXTmHPtWVTH7vdTbMvpFtF14RzO0jDQYcffvF7Tp0OuielUqqSU+qUU91vZtdjfZIMOgFpcU5wTTODdJMhqhSOAzpnomk5UFfIQQDu2WgV91CHZ6QHWxFzoQJE+TixYsTXY0uTdMka3f754qDP4fMpH55Lc4JX1FWm9h9ezWNwnmfM/CRv5G+dRPe9AwM9TY8WdnsyM2n77bN6IEFwPU6He7xk3jwv+8hhEDTtKbpnH0+SEB6h8b65Kc1y8+01+pkXbkVIWBsn9xmyfXaSkqpZvZ0IkKIJZHWaakWgdJma8rrmuTy8Xg11pTXMa7PvvzAmiapd3uxu33UOTxxDwI6l5MeH75FwTdfkLZ9K6ll29F53Ozp1p1HDj+al/fu4YXZNzDgq08x/rSA54p6sO7wqRROnc5fxk0iPWPfTbZZTv8EBwG9XjSZ09+gW6aZnFQTmpTtslBLBYGuQwUCpU12W5xhE7pV29xsqazHZNBRaXVRHcgsGk/C7cK8u5zCeZ/R579PYa7cS33/QezIyeOL3bv40ePmnb27MdmsDBsznk0Dh+KcehwAAwM/yaZ3XioOty+YWA+gbyDNQziRXleUlqhAoERkd3ub9Cc7PT7W7Y6cdXJTLGcESUlaYC9ee/9BwZdTdpXR74kHyP3mS9Jrq4OvfwVsufgq+l17K9s2reej5x5j2JjxPD72IPoNHpaQzJ1tlZdhYkhRJpqEZdtrqLV7MBp09MnrOJvOKB2DCgRdlNXpweXV8GkSo17XLPlXhdXFyp21ZKeaGNo9kzSTntW7LHFb7au31pFWto3U7VvJWbqQgm+/JG3HNgCq+vZn4eDhFAoY982XICWve9xsAPaaU7EMGY7+4ClMmXYiCEHfgUO4+f4n4lLv1nTLSqHG7sHjbTn9aVqKnlG9shFCoBcwpjiHxVtr6JWTGtyHQVHaiwoEXUzDXP7Q1aONN3WpqXezcmctmgY19W4WbqkiO9VETb0nwlEPqEJkr1hKty8+Jn3jesx7dpGydzfGOkuwiC/FzJL8Aj7u0Qvr3j2csW0zx2/bjA/Yc9JMNs2+kfUrlzJw6CgO71OSFH3b4VIyG/SCYT2ysDq9LNteEzFzqtmop7Q4p8m0TaNex9g+OZhUEFBiQAWCDspi96DTEfUcfvB39awss4TN17+j2o7F4aEkP43V5XVNbmINAaFdSIl553ayf11KzrJfKJz/BebdO/EZjJT36MkWvYGtmVmsBWoLu3PG/U9g7zeAG2afhxCCAdNn8OnwMezMy6d7z2JcgXn7R/RMjl1OhfCvxO1XkM6vOyzUOfYFzwGFGcHW1+CizGZJ9dJTDJQUpFGUaQ67Uldl6lRiRQWCDkJKidXlpdLqYrfFid3tizhXvIGmSaxOL9V2N9X1biwOd4sbstQ5PKwos0QucAB0Lifd3nyZ4hefIntPOQAOgwHroUex8dpbuOrDt/jmpwXoDQaKSwZQMnocI8ZOxDZsJAD/fP7tZsdM1GRUvU4E5++bjXq8mobHJ5FS0r8gI5iyYXTvbBZtrcbl0UhL0dM7d99eC8V5adhcXmwuL3npJgrSU5qkelCUeFKBIMnV1LvZXm2nxu5u1j/v9UmWba9ldO9s8jNSggndqgObwludnrjvxKVpGpW7d1G2fQs1a1dyyuDh5K1YSt7zj5PrsPMT8DKwUAiqe/Xhmcf+i95gYObAIZyl09O73wBMpsSsVM1KNdI9y4zN5aXO6aHe5W3SfZOdZqR3bmrEb+yhzEY9o3vnsGRbNYOLMpt1WQ3rkdXel6Ao+0UtKIuRLZX1ZJkN+7383unxsXGvLeLeu43pdJCbZqLW7gm7m9d+0zT/wUN46+txLPsFuXwRQ4WOzNpq6lYtx75+LSkuJ6lSkgH0aPSZLUOG8+bQkVgPPYq+A4fQu6Q/ppTWUxXHkhBg0OvISNHTNz+dgpC/lZQSj0/iDUTT/d1U3eLwkJ2qvu0riaUWlMXZ+j1WtlfZEQIGdcukT75/up/D7WN7tR2fJslJM5KbZkIIqHf5F145PD7cXg23T8PSaHP31mgaTeaZt4WpYi/ZK5bg7Nmb+pIBSIORnK//j8K3X6XXoh9BCHzpmbgMejyWWlJ8PtJDvjx4srJJTc9gmzkVZ7ciPFk5OHLz2TN0BL6Jk7GNGI0nJ494bz3XMyeVkoI0dEKg1wkEIAlmjsagEy0OLAshMBkEJg5sgFYFASXZddlA4PVp7T4NT0rJ2nJrMOeLlP6gUOf0DxjuqXMGuxoa54WJmcDAbNaq5WStWk7qzh3U5xVSkZ2Nrd5G358W0H/LRnSBb7waUI8gE0k58G+g9MhjGdC9J7byMlasXoEuOxt9bgHukgFoI8eQM+UYTEU9WqpFTOl1/nw3uekmymsdlFucCOHvdgn9hq8oSnhdKhBU2VzsqnVS63Dj8mhkpxnplpnSZC/W/eH1aeyuc7KzxhF2Rk403TutkhJ9vY2Uij2Yd+8iZfdOjLU16LxehMeDlJK9lhrqvB581VX03LSefrvKyHH6A47PaGSz10tPKekfOORvwKcHH0732Tfi3rCOVf99kt4GI+uHj2LPpMMo6N2H6iEj8Ob4U0fkH/hVtJlBLyjKMpNuMpCWokfTJFX1/nEQs1HP8B5ZpJr8f7tBRZkM7JaBT5Nqrr2itEHMxgiEEMX4xwWL8LfIn5VSPhpSRgCPAscDduACKeXSlo57IGME63bXUVbd/Ju4EP7NtEsK0slIMeDxaZTXOtllceD2amiB35F/cY9ApwN9oLtBpxNYwvTN6611oNPhS0sPZqjE58NUVYGxzoLeXo/eYUfndKB3OtA7ncjaamTFHkRlBYaqSow1VWQ67GQ5nRirK9G7o58n8xuwyGDAe+hRDLvyRqr79uPZx+6noLCIPukZFKVnYBgxhm49e2NObZ+VqhlmA05P+20xmWk2MLp3TvBGryjK/kvUGIEXuF5KuVQIkQksEUJ8KaVc06jMccCgwM8k4KnAf9vfli1kv/8pPqeGNBjQjEa8GVl4M7PwZmRi1XysdrvJ0mvY7S7weDH5fKTW2zBY6zDYbWguNz6PG83rwa3X49Yb8aaYKczJQe+wU7NpPRnr19Jzy0YKKvcC4Daa8BZ1R+d2Y6rYgy6KwGsBqoEKYE9eASWHH407r4BXP3yT3UBddi51+YVo3boxatIUjj5pJgDL5n9OjjmNjMJuZJYMoLvZP13RChiB2bfcG4vfLClGHYO6ZdI924ymSSrrXeyxuHB4fPg0iSYlqSY92alGssxGTAadP4gKgv33eiGwurzsqXOyp85JXrqJYd2z1M5XihIHcZs1JIT4EHhcSvllo9eeAb6RUr4ReP4bcKSUsjzScfa7RfDWWzBzZts/10blwELgF8ADDM3N48RDj8RnTuXzRT+x3mrBkZaOLy0NLS2DbgMHc/z5l6GZzXwx/wtsBiPpeXlkZGaRmZVDXmER3Xr2Ch6/YaCz8XPwv6bXCYx6HTrhn1rqk7LJ9FGd8CclMxv1mAw6jHr/PrVGvQ4hoOGwXk3i0yQen4bD48Pp8eH0aKQEPptq0qMXInjOHtmpKie9oiS5hM8aEkKUAGPx3yMb6wXsaPS8LPBak0AghLgUuBSgT5/93Pf197/njfvv5/tvv8UoJQafjwxNI1vTuO6ii8Bo5Mtvv2XNxo1gNKIzGsFgwJCby59uugmysvj488/Ztn07RpOJNJ2OdCA3JYWjpk2D9HTWbt+OXa9nSHo6pamppKenk5GRQWqq/5v5ha1UcdIxB+3ftSmKohyAmLcIhBAZwLfA36SU74W89wlwn5Ty+8Dzr4GbpJQRv/J3lHUEiqIoyaSlFkFMp1YIIYzAu8BroUEgYCdQ3Oh578BriqIoSpzELBAEZgQ9D6yVUj4codhHwHnC72DA0tL4gKIoitL+YjlGcChwLrBSCLE88NqtQB8AKeXTwP/hnzq6Ef/00da60RVFUZR2FrNAEOj3b3EqifQPUMyOVR0URVGU1qnll4qiKF2cCgSKoihdnAoEiqIoXZwKBIqiKF1ch9uYRghRAWzbz48XAJXtWJ2OQF1z16CuuWs4kGvuK6UsDPdGhwsEB0IIsTjSyrrOSl1z16CuuWuI1TWrriFFUZQuTgUCRVGULq6rBYJnE12BBFDX3DWoa+4aYnLNXWqMQFEURWmuq7UIFEVRlBAqECiKonRxnTIQCCGmCyF+E0JsFELcHOb9FCHEm4H3FwZ2UOvQorjmPwsh1gghVgghvhZC9E1EPdtTa9fcqNypQggphOjwUw2juWYhxBmBv/VqIcTr8a5je4vi33YfIcR8IcSywL/v4xNRz/YihHhBCLFXCLEqwvtCCPFY4PexQggx7oBPKqXsVD+AHtgE9AdMwK/A8JAyVwBPBx6fCbyZ6HrH4ZqPAtICj//UFa45UC4TWAD8DExIdL3j8HceBCwDcgPPuyW63nG45meBPwUeDwe2JrreB3jNU4BxwKoI7x8PfIY/u/PBwMIDPWdnbBEcBGyUUm6WUrqBucAfQsr8AXgp8PgdYGpgI52OqtVrllLOl1LaA09/xr8bXEcWzd8Z4B7gfsAZz8rFSDTXfAnwhJSyBkBKuTfOdWxv0VyzBLICj7OBXXGsX7uTUi4Aqlso8gfgZen3M5AjhOhxIOfsjIGgF7Cj0fOywGthy0gpvYAFyI9L7WIjmmtu7CL83yg6slavOdBkLpZSfhrPisVQNH/nwcBgIcQPQoifhRDT41a72Ijmmu8EzhFClOHf7Oqq+FQtYdr6/3urYrlDmZKEhBDnABOAIxJdl1gSQuiAh4ELElyVeDPg7x46En+rb4EQYpSUsjaRlYqxs4AXpZQPCSEOAV4RQoyUUmqJrlhH0RlbBDuB4kbPewdeC1tGCGHA35ysikvtYiOaa0YIcQxwGzBDSumKU91ipbVrzgRGAt8IIbbi70v9qIMPGEfzdy4DPpJSeqSUW4D1+ANDRxXNNV8EvAUgpfwJMONPztZZRfX/e1t0xkCwCBgkhOgnhDDhHwz+KKTMR8D5gcenAfNkYBSmg2r1moUQY4Fn8AeBjt5vDK1cs5TSIqUskFKWSClL8I+LzJBSLk5MddtFNP+2P8DfGkAIUYC/q2hzHOvY3qK55u3AVAAhxDD8gaAirrWMr4+A8wKzhw4GLFLK8gM5YKfrGpJSeoUQVwL/wz/j4AUp5WohxN3AYinlR8Dz+JuPG/EPypyZuBofuCiv+QEgA3g7MC6+XUo5I2GVPkBRXnOnEuU1/w84VgixBvABN0opO2xrN8prvh54TghxHf6B4ws68hc7IcQb+IN5QWDc4w7ACCClfBr/OMjxwEbADlx4wOfswL8vRVEUpR10xq4hRVEUpQ1UIFAUReniVCBQFEXp4lQgUBRF6eJUIFAUReniVCBQlFYIIXKEEFcEHvcUQryT6DopSntS00cVpRWBNOWfSClHJrouihILnW5BmaLEwH3AACHEcmADMExKOVIIcQFwEpCOP43Dg/hTJZ8LuIDjpZTVQogBwBNAIf4FQJdIKdfF+yIUJRLVNaQorbsZ2CSlLAVuDHlvJHAKMBH4G2CXUo4FfgLOC5R5FrhKSjkeuAF4Mh6VVpRoqRaBohyY+VJKK2AVQliAjwOvrwRGCyEygMnsS+0BkBL/aipKZCoQKMqBaZzFVWv0XMP//5cOqA20JhQlKamuIUVpnRV/Wus2k1LWAVuEEKdDcL/ZMe1ZOUU5UCoQKEorAtk7fwhsJv7AfhxiFnCREOJXYDXht9RUlIRR00cVRVG6ONUiUBRF6eJUIFAUReniVCBQFEXp4lQgUBRF6eJUIFAUReniVCBQFEXp4lQgUBRF6eL+HxHIOPw5F3gVAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# create the two random arrays of shape 50 x length of the time domain\n",
    "mult_rand = np.sqrt(0.001) * np.random.randn(50, 101) + 1.\n",
    "add_rand = np.sqrt(0.005) * np.random.randn(50, 101)\n",
    "\n",
    "# add noise to the t variable, this will come directly from the equation above\n",
    "# it's a linear transformation of a*t + b, but using the appropriate NumPy syntax\n",
    "t_rand = np.multiply(t, mult_rand) + add_rand\n",
    "\n",
    "# use the same syntax as before to process the randomized time data\n",
    "# this should be a direct copy from what you have in the previous code snippet\n",
    "random_response = 0.5 * np.exp(t_rand**2) + 1.5\n",
    "\n",
    "# find the min, max, and average response at each time step\n",
    "r_max = np.max(random_response, axis = 0)\n",
    "r_min = np.min(random_response, axis = 0)\n",
    "r_avg = np.mean(random_response, axis = 0)\n",
    "\n",
    "# this plots the output generated above as \n",
    "fig = plt.plot(t, expected_response, 'k--', t, r_avg, 'r')\n",
    "plt.fill_between(t, r_min, r_max, alpha = 0.3)\n",
    "plt.legend(['Expected', 'Generated Mean', 'Min/Max Range'], loc = 'upper left')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('position');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's finally do some procesing so we can finalize our data for output!\n",
    "We can pretend that the 8th and 23rd sample ranges were done by a trainee that had set the sample frequency to be a little too slow.\n",
    "We are also limited by our instruments in that they are unreliable when the position is greater than 3.25 so we will remove any samples that reached more than 3.25 at any time in the experiment.\n",
    "**Side note:** this is not a real example, and data should only be excluded when appropriate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of samples left after removing the bad trainee data is 48\n",
      "The number of samples left after removing the data > 3.25 is 36\n",
      "our final_data for output has shape (52, 101)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABTZ0lEQVR4nO3dd3hUVfrA8e+Zlslk0hNqCAEE6R3EhgUruqxd7NhX3VV37b27rrr2tWDvuPb2s67YRQUEpDcpoSQhdXq5c35/zGSYTCbJBDIlyfk8Tx4y996Ze24S7ntPe4+QUqIoiqJ0X7pUF0BRFEVJLRUIFEVRujkVCBRFUbo5FQgURVG6ORUIFEVRujlDqgvQXkVFRbKsrCzVxVAURelUFixYsENKWRxrX6cLBGVlZcyfPz/VxVAURelUhBAbW9qnmoYURVG6ORUIFEVRujkVCBRFUbq5TtdHEIvP56O8vBy3253qoii7wGw2U1JSgtFoTHVRFKVb6hKBoLy8nOzsbMrKyhBCpLo4SjtIKamurqa8vJwBAwakujiK0i11iaYht9tNYWGhCgKdkBCCwsJCVZtTlBTqEoEAUEGgE1O/O0VJrS4TCBRFUZRdowJBB9Hr9YwdOzb8dc899yT8nHV1dTz++OPtft+tt97K/fffn4ASKYoSze3TUl2ENnWJzuJ0kJmZyaJFi5J6zsZAcPHFFyf1vIqixG9rnYuBxdZUF6NVqkaQQPX19ey5556sWrUKgFNOOYWnn34aAKvVyt///ndGjBjBtGnTqKqqAmDdunUcccQRTJgwgf3335+VK1cCUFFRwbHHHsuYMWMYM2YMP/74I9deey3r1q1j7NixXHXVVQDcd999TJo0idGjR3PLLbeEy3LXXXcxZMgQ9ttvv3B5FEVJLC0g2d6Q/gMhumSN4MADD2y27aSTTuLiiy/G6XQyffr0ZvtnzZrFrFmz2LFjByeccEKTfV9//XWb53S5XIwdOzb8+rrrruPkk0/mscceY9asWVx22WXU1tZy/vnnA+BwOJg4cSIPPvggt99+O7fddhuPPfYYF1xwAU8++SSDBw/m559/5uKLL+arr77i0ksv5YADDuDdd99F0zTsdjv33HMPS5cuDddEPv/8c9asWcMvv/yClJIZM2bw7bffkpWVxZw5c1i0aBF+v5/x48czYcKEuH+eiqLsmjqnF5dXIxCQ6HTpOygiYYFACGEGvgUyQud5S0p5S9Qx/wDOA/xAFXCOlLLFxEjprKWmoUMPPZQ333yTSy65hMWLF4e363Q6Tj75ZABOP/10jjvuOOx2Oz/++CMnnnhi+DiPxwPAV199xUsvvQQE+yNyc3Opra1tcq7PP/+czz//nHHjxgFgt9tZs2YNNpuNY489FovFAsCMGTM67sIVpZvxawEM+vgaU+pcPqQEp0/DmpG+z92JLJkHOFhKaRdCGIHvhRCfSCnnRRzzGzBRSukUQlwE3AucvLsnbu0J3mKxtLq/qKgorhpAvAKBACtWrMBisVBbW0tJSUnM44QQBAIB8vLydrmvQUrJddddx4UXXthk+0MPPbRLn6coSnMOj0auJc5A4PQC4PT60zoQJKyPQAbZQy+NoS8ZdcxcKaUz9HIeEPsu2Yk9+OCDDBs2jNdee42zzz4bn88HBAPEW2+9BcBrr73GfvvtR05ODgMGDODNN98Egjf2xlrEtGnTeOKJJwDQNI36+nqys7Ox2Wzhcx1++OE899xz2O3BH/uWLVuorKxk6tSpvPfee7hcLmw2Gx9++GHSrl9Ruhq71x/XcYGApMEVPNbpSe+RQwntLBZC6IUQi4BK4Asp5c+tHH4u8EkLn3OBEGK+EGJ+Y6dqumnsI2j8uvbaa1m1ahXPPPMM//73v9l///2ZOnUqd955JwBZWVn88ssvjBw5kq+++oqbb74ZgFdffZVnn32WMWPGMGLECN5//30AHn74YebOncuoUaOYMGECy5cvp7CwkH333ZeRI0dy1VVXcdhhh3Hqqaey9957M2rUKE444QRsNhvjx4/n5JNPZsyYMRx55JFMmjQpZT8nRensHJ74AkGD24cWCD77OuIMHgDVds8ulWt3CCll20ft7kmEyAPeBf4mpVwaY//pwF+BA6SUrf4UJk6cKKMXplmxYgXDhg3ruAIngdVqDT+5K53zd6h0Tws31TK+NL/N4zbscLC2Mvh/PNdiZFJZQVyfv3RLPYOKrWSa9LtVzmhCiAVSyomx9iVl+KiUsg6YCxwRvU8IcQhwAzCjrSCgKIqSavHWCGpD/QPteQ9AvctHtaP5rdCnBeL+jPZKWCAQQhSHagIIITKBQ4GVUceMA54iGAQqE1WWdKRqA4rS+QQCEo8vQCDQekuKlJJ6ly/82q9JvP62b+RefwCXV6Pa7m22ryKB8xESWSPoDcwVQiwBfiXYR/CREOJ2IUTj+MX7ACvwphBikRDigwSWR1EUZbf4AsGbudvfeuevzePHrzUNFs44+gkag0eN00tks30gIKm0Ja7BJGHjmaSUS4BxMbbfHPH9IYk6v6IoSkdrvLl7fAEsppaPq3f6mm1zeDXyLK1/fmMg0LRgjSIvdJIquyfc8ZwIKsWEoihKnPyhm3FbNYIGd/NA4IyjnyCyOWlHRPPQ9vrEpqlQgUBRFCVO/lCHrcfXenu/I8a8Aae39eAhpWwSQGocwUDg0wIxO487kgoEHSQ6DfWGDRvYZ599ANiwYQMjR44EYNGiRfzf//1fwspx4IEHEj28tnF7aWlpk3bHY445Bqs1vbMiKko6ibdGEGuUUPRcgugOZ4dXQ4voV7C5fXj9ASoa3AQSN2AIUIGgwzTmGmr8Kisr48cff2x23K4EAr8//qFnrcnLy+OHH34Agimst23b1iGfqyjdRWMgaK1G4PJqMdvz3T4t/CCmBSTzN9biiQgokc1CAFIGawWJbhYCFQgSKvpp2+v1cvPNN/PGG28wduxY3njjDRwOB+eccw6TJ09m3Lhx4ZnEL7zwAjNmzODggw9m2rRpLR7ncrmYOXMmw4YN49hjj8XlcrVYnpkzZzJnzhwA3nnnHY477rgm+1tKYX3MMccwYcIERowYwezZs5tc3w033MCYMWOYMmUKFRUVu/cDU5Q019g01NpiM/YW+gICAXCF3rem0kaDyxeecAaxO5i31Lmoi7G9o6VvFqRddfnl0NELxIwdC20kbotMQz1gwADefffdZseYTCZuv/125s+fz2OPPQbA9ddfz8EHH8xzzz1HXV0dkydP5pBDgoOpFi5cyJIlSygoKGjxuKeeegqLxcKKFStYsmQJ48ePb7GM06ZN4/zzz0fTNObMmcPs2bO54447gJZTWE+dOpXnnnuOgoICXC4XkyZN4vjjj6ewsBCHw8GUKVO46667uPrqq3n66ae58cYb2//zVZROIlwjaGVOQGuTxxweDZdXo7wm+MC2rc5NSZ6FXIuxWY0AoNbRfD5BInS9QJAiu7pC2eeff84HH3wQXjrS7XazadMmIJjCuqCgoNXjvv32Wy699FIARo8ezejRo1s8l16vZ7/99mPOnDm4XC7KysqalCNWCuupU6fyyCOPhAPb5s2bWbNmDYWFhZhMJo4++mgAJkyYwBdffNHu61eUzqRx+KjXH2hxjYGWagQQbP7ZVt+01r5yewPj++fHNc8gUbpeIOhkKZellLz99tvsueeeTbb//PPPZGVltXlce82cOZNjjz2WW2+9tVk5YqWw/vrrr/nyyy/56aefsFgsHHjggbjdwTZLo9GIEMH/CHq9vsP6MhQlXUWmefD4AzHzAbVWI9hY7SA6vZvN7WfFtoZm25NJ9REkWazU0Y8++mi4E+m3336L+b6Wjps6dSqvvfYaAEuXLmXJkiWtnn///ffnuuuu45RTTmn2+bFSWNfX15Ofn4/FYmHlypXMmzcv1scqSrcQ2QnsiTFySErZ6jDRlm72lQ2pTbOmAkGSHXTQQSxfvjzcWXzTTTfh8/kYPXo0I0aM4Kabbor5vpaOu+iii7Db7QwbNoybb765zSUohRBceeWVFBUVNdneUgrrI444Ar/fz7Bhw7j22muZMmVKx/wgFCUNtTV71x8xjtMdY+SQyxd7xFC6S0oa6o7UVdJQK02p36GSDuqdPnItxhb3z1tfjd0dbPoZ3NNK/8KsJvsrbW6WbK5PSNnak8o6lpSnoVYURekMbJ7Wh2o2bRpqXiOINaO4M1CBQFEUheBN3tVGGojIzuJYcwnas+5AOlGBQFEUheBNvrX5AdB2jaC1oaPpTAUCRVEUwKsFWp0x7NcCTUb9RB8bHDGkAoGiKEqn5fO3XiPwR40G8voDTZI4unxawpPDJYoKBIqiKIBPkzHnBjSKDgRSNm0e6qzNQtAVZxYDXy7v2ORnhwzv2eYxQghOO+00XnnlFSCYMbR3797stddefPTRR3zwwQcsX76ca6+9ts3POvLII3n66ac5/fTTWb9+PRs3bgzP4D3mmGP48ssvO2TNY71ez6hRo/D7/QwYMICXX36ZvLy83f5cRemMfFqAQCD4pG8yNH9G1rTmQ+09vgBmY3B2cWcdMQSqRtBhsrKyWLp0aTj75xdffEHfvn3D+2fMmBFXEHC5XFRXV1NSUgIkNnV0Y36kpUuXUlBQwH/+858O+2xF6Wy8jYvOtFAr8MVo93H5NCob3CzeXMcfO3b/4SxVVCDoQNOnT+fjjz8G4PXXX2+SxuGFF17gr3/9KwCzZs3i0ksvZZ999mHgwIG89dZb4eO+/vprDjzwwPDr1lJH2+12pk2bxvjx4xk1alQ4NfWvv/7K6NGjcbvdOBwORowYwdKlS1st+957782WLVsA+OWXX9h7770ZN24c++yzD6tWrQpfw3HHHccRRxzB4MGDufrqq8Pvf/bZZxkyZAiTJ0/m/PPPD19rVVUVxx9/PJMmTWLSpEnhoKYo6cYXDgSxG/qjF6MHWLa1niXl9VTZPJ22fwBUIOhQjTdtt9vNkiVL2GuvvVo8dtu2bXz//fd89NFHTWoKn3zyCUcccUT49bRp0/j222/DqaNPPvnk8D6z2cy7777LwoULmTt3LldccQVSSiZNmsSMGTO48cYbufrqqzn99NPDK6TFomka//vf/5gxYwYAQ4cO5bvvvuO3337j9ttv5/rrrw8fu2jRIt544w1+//133njjDTZv3szWrVu54447mDdvHj/88AMrV64MH3/ZZZfx97//nV9//ZW3336b8847r30/VEVJEp+/9RTT/hh3+k6WmKFFXbKPIFVGjx7Nhg0beP3115k+fXqrxx5zzDHodDqGDx/eZEGXH374IZxqGlpPHS2l5Prrr+fbb79Fp9OxZcsWKioq6NWrFzfffDOTJk3CbDbzyCOPxCxD4xoKW7ZsYdiwYRx66KEA1NfXc9ZZZ7FmzRqEEPh8O2dbTps2jdzcXACGDx/Oxo0b2bFjBwcccEA4ZfaJJ57I6tWrAfjyyy9Zvnx5+P0NDQ3Y7Xa1RKaSdsJNQy0MIY1VI+gqVI2gg82YMYMrr7yyWXbPaBkZGeHvG4egrV+/nn79+mEymZocO3PmTC699FJOOumkJttfffVVqqqqWLBgAYsWLaJnz57hFNHV1dXY7XZsNlt4W7TGPoKNGzcipQz3Edx0000cdNBBLF26lA8//LDJ+yPLHU/q6UAgwLx588JLeG7ZskUFASUt+cKrj7VUI1CBoN2EEGYhxC9CiMVCiGVCiNtiHJMhhHhDCLFWCPGzEKIsUeVJlnPOOYdbbrmFUaNGtfu90c1CjVpKHV1fX0+PHj0wGo3MnTuXjRs3hvddeOGF3HHHHZx22mlcc801rZ7XYrHwyCOP8O9//xu/3099fX24o/uFF15os9yTJk3im2++oba2Fr/fz9tvvx3ed9hhh/Hoo4+GX+/K4j2Kkgy+NjqLYzUNdRWJbBryAAdLKe1CCCPwvRDiEyllZEL7c4FaKeUeQoiZwL+Ak2N9WHvEM9wzUUpKSsIrhrXXp59+2uSm2agxdXS00047jT/96U+MGjWKiRMnMnToUABeeukljEYjp556Kpqmsc8++/DVV19x8MEHt3jucePGMXr0aF5//XWuvvpqzjrrLO68806OOuqoNsvdt29frr/+eiZPnkxBQQFDhw4NNx898sgjXHLJJYwePRq/38/UqVN58skn4/2RKErS7EpncVeRlDTUQggL8D1wkZTy54jtnwG3Sil/EkIYgO1AsWylUF01DbXH42Hfffcl+to6i8Z2f7/fz7HHHss555zDscceG/f7u8LvUOm8/FqAr1dVAWA06DhgSHGzYxZtrmOHLXULyHTaNNRCCL0QYhFQCXwRGQRC+gKbAaSUfqAeKIzxORcIIeYLIeZXVVUlssgpk5GR0WmDAMCtt97K2LFjGTlyJAMGDOCYY45JdZEUJW6+iKd9X2g94mh+TTUN7RIppQaMFULkAe8KIUZKKVsf0B77c2YDsyFYI+jYUiodIXKkk6J0Nt6om3ys9YhVZ/FuklLWAXOB6J7QLUA/gFDTUC5QvYvn2I0SKqmkfndKskX/zfmaBYLmHcZduY8gkaOGikM1AYQQmcChwMqowz4Azgp9fwLwVWv9Ay0xm81UV1erG0onJKWkuroas9mc6qIo3YgtKkFcdCCINYQ0VoqJriKRTUO9gReFEHqCAee/UsqPhBC3A/OllB8AzwIvCyHWAjXAzF05UUlJCeXl5XTV/oOuzmw2h3MrKUqiSSlpcPnIMe9cm7hxVnGj6BqBlDJm0rmuImGBQEq5BBgXY/vNEd+7gRN391xGo5EBAwbs7scoitINeLVAs0yhsfoIInXl/gFQM4sVRekkOqrp169JHFEriXmjbvyeqKYhTQUCRVGU1OuohV98WgBnVI2grc7i6P1djQoEiqJ0Ci6v1iE35Ma1iSOf8psHAlUjUBRFSTsef+uLy8ercfJYZPNQ8z6C6BqBCgSKoigp52ljcfl4+UKfEdk8FH2jb1yyspGqESiKoqQBj1/roBpB8AbfWCOQUsZMHxFZK1B9BIqiKGmgvTUCp9cf8wbf+PTfWCPwaTLmSmPb63euw6GGjyqKoqQBbzv6CAIBye/l9TEDR3SNoKWn/Y3VTqrtwWyjWheeVQwqECiK0km0p0awtsqOze1vNRA4vX6klK02+yzb2oDbp6nOYkVRlFQLBCS+OGsEVTYPm6qdQPOJYrBzhFAgEMwpFD1iqMmx/gDLtjZ06YRzoBavVxSlE/C2sXpYo0BAsnxbQ/h1rCyikU/3Dq+/zaf9WocXXRd/ZO7il6coSlfQmPJB02KP8GlU7/KFh4dC7MAR+X6nR2tyfEuS1UWg87jZ58i96fHp+8k5YeN5k3o2RVGUXeDRdj7Zu1u5cde5fE1eRzcN+bRAkxFCdo8/rYaGZm78A8umP+j/4lNJPa8KBIqipL3IJHCeVvoJ6pzepu9rI2eQ0+tvtY8g2bI2rAMgd8lCLOtWJ+28KhAoipL2Ipt4WqoRSCmpj6oRRGcRjV53wOFNrxFBlo3rAQjo9fR5740m+3o/ej8sXJiQ86pAoChK2ots4mmpRmD3+JuN7vFEPe1HP/37/AGcHZTVtCNYNqzD3aMX1ftPo/eHbyH8wbIVfvMFJQ/8E955JyHnVYFAUZS0F9nEE2sZSYA6p6/ZtujO5Vj9AU7v7qet6CiWDetx9h/I1mNOJqOqgoIfv8bQUM+eN/0D+x57wk03JeS8KhAoipL2PE1GAsW+cUc3CzWKrAWk+3wAy8Z1OMsGseOAQ/HmFdDnvTfY4583YqquYpYMIE2mhJxXzSNQFCXtRTYNtadGAMF+Akvo/plOHcPRDPV1mGprcJYNRJpMbD/6OEpefx6dpnEXMPmqGxFCJOTcqkagKEpak1I27SOIUSNw+1rOTBpZm0inoaLRGjuKnaUDAdh6zEx0msZS4NfjTuWgw49O2LlVjUBRlLTWbCF5TaIFJHrdzqfjlmoD0LQ20SkCQVkwEGzp1Ztrs3OYl5PHFdfdkdBzqxqBoihpLdbs4Oin/zqXt9kxO9/fOdYVsGxYh9TpcJX0B0Dza3w1ZgIzH36WTEtWQs+tagSKoqS1WE1BHn+ArIydr1urEUQGEq8/fTuLLRvX4+pbGu4QLijuwd1PvZ6UcyesRiCE6CeEmCuEWC6EWCaEuCzGMblCiA+FEItDx5ydqPIoitI5RU8Kg6Y1ArdPw9HKXIDIQOBP43UFLBvW4ywbyKqli7jhotOora5K2rkTWSPwA1dIKRcKIbKBBUKIL6SUyyOOuQRYLqX8kxCiGFglhHhVStlyPU9RlG4l1kifyJv76gpbzBXGdh7bCZqGpMSyYR1Vo8Zx15UXovn96PXJa7BJ2JmklNuAbaHvbUKIFUBfIDIQSCBbBMdEWYEaggFEURQFaL1GUGXzUNngafX9jZ3Ffi2QtCyi7WXaUYnB5eT9xfOp3LaFf7/wLjl5+Uk7f1I6i4UQZcA44OeoXY8Bw4CtwO/AZVLKZr8qIcQFQoj5Qoj5VVXJqy4pipJ6LfURaAHJqu22Nt/fOMoonXIKRbOEks19smo55/39RkaMm5TU8yc8EAghrMDbwOVSyoao3YcDi4A+wFjgMSFETvRnSClnSyknSiknFhcXJ7jEiqKkk1irjLl9Guur7HGvYezxa/jStToAZIQyjWbvfQDHn3Vh0s+f0EAghDASDAKvSiljZUs6G3hHBq0F/gCGJrJMiqJ0LrGGj7q8GptqnHF/htcfiGsBmkQSXi9957yAqaqy2b7s8o1oJhNn3v9kwmYPtyaRo4YE8CywQkr5QAuHbQKmhY7vCewJrE9UmRRF6VxaWlxeC8hWO4ijefyBlDcN9X/xSYbdcS2TTz6cnCXBdNKapvHeq8+SuX4trtIBWJPYLxApkd3S+wJnAL8LIRaFtl0PlAJIKZ8E7gBeEEL8DgjgGinljgSWSVGUTsTjD7Trht/i5/gCpOBBOyxj+1YGPPUgtZP2JmPbFiaeeQwrb76HOzasZ86zj3FLz944R41LWfkSOWroe4I399aO2QoclqgyKIrSucUaMbQrvJqWkiaXRoPvvx0CkmV3PoRmzWbklX9h+E3/4CHg8t4lFFZVsKn/wJSVT6WYUBQlbTl9HTOa3O0LpGwOQd6vP9Lrk/fYeO4luEv648sr4INrbuMvpgx+yCtghDkDofmpHzsxJeUDlWJCUZQ01lGLxni1AFIm/7lX+P3sefcNuPqUsOHcvwLBfo9/Xn8ZVVYrU9/4jHl9ShA+H9JoTHr5GqlAoChK2nJ6OiYQeHwBSMF9tucn75G9egVLHniagDkTACEEV975EC6ngx59SgBSGgRABQJFUdKY09sxTUMev4Yu2V0EUtL/+SewDxpC5aFHAbDq99/Yc9Q4Bu45PMmFaZ3qI1AUJW0545ww1ha/JnEneR5BwbzvyF61jE2z/gI6Hd98+j5/nXkkX374VlLLEQ8VCBRFSUtun4bWgWP/kz2hrPT5J/AUFrPt6ONZtXQR915/GSPGTWbq4X9KajniEVfTUCgz6PlAWeR7pJTnJKZYiqJ0Bn4tgEGfmOdJVwd1FKeCddVyin6Yy9pLr6WqrpZb/jaL/MIibnn4WUymjLY/IMni7SN4H/gO+BLovL8dRVE6lMOjkWtJTCDoqGahVCh98Um0zEw2n3g6t11yJg5bAw+/+hH5hemZKy3eQGCRUl6T0JIoitLp2Dw+ci27P+JFStlswpezlcVm0pl11XJ6ffwuW046A62giJPOuQS93pB2HcSR4g0EHwkhpksp/y+hpVEUpVOxd9DNusrmoUeOucm2jppDkCzC76f0+ScY9J/78Ofk8MthR5ML7B8aMZTO4q3TXUYwGLiFELbQV3RKaUVRuhm7u2MCwQ67lwZ303WHO1MgMNbsYMKZf2bwQ3dRdeCh3Hn+ZZx83kn8vmBeqosWl7hqBFLK7EQXRFGUzsfp1dACEv1uDtJvcPuosevJMQebmaSUuDoovUQylD39KDlLF/H7vx7nk4JC7r7oNMZO3o9hoyekumhxiXtCmRBiBjA19PJrKeVHiSmSoiidhU8L4PFrWEy7PjdVC0gcHj/VDi9lRVlAMDdQGq8j04Shrpa+b75MxfRj+XGPPbntzD9TOnAwNz0wG0OKZwzHK66mISHEPQSbh5aHvi4TQvwzkQVTFCW9+bVgimj3bmYItbv9SAn1Li9aIDhvoKNmFHc0c/lGBj10Nzq3K7yt3+vPY3A5WXrCadxw0WlYrNnc9cSrZGU3W2wxbcUbxqcDYxvXExZCvAj8BlyXqIIpipLeGhd6ibWmcHs09g0EAlDr9FJkzUjb/oG+b77CgGceJbN8I0vvfQKdx02/V5+h6oBDYMIUZpxyNnsfeBjFvfqkuqjt0p76XB5QE/o+t+OLoihKZ9K4BvDu1gjqXTs7iWsc6R0I8n/9ES0zk16fvI+zdCDewiJMtTUsOPoEDEJw6gWXpbqIuyTeQPBP4DchxFyCi81MBa5NWKkURUl7jSkbdrdGYIsYeVRt90LP9Gwa0jsc5CxbzMZZF2GqrWbgUw/is2azJDePU+6+nuf3OwhrTud8Ro531NDrQoivgUmhTddIKbcnrFSKoqS9cNPQbtQItIBsctN3ePy4fVpappfIXfQrOr+f2kn7UDt5X8zlmyj8+XtuBE697s5OGwSgjc5iIcTQ0L/jgd5AeeirT2iboijdVOOKX+7dSAVhc/uarUlc7fDiSsP0EvnzfyKg11M/bjLSZOKqYaOYDljPv5RjTz8v1cXbLW3VCP4BXAD8O8Y+CRzc4SVSFKVTaAwEnt3I6tngat4EVF7j7JAF6zta/i8/YhsxBi0ri68+fpfnX3iC6SecztmXdf4xM60GAinlBaFvj5RSuiP3CSHMMd6iKEo30dg05PUHCAQkul2YVBY9mxia9hmkC53TQc7S39h01l8A2HfaEVxw1S0cd8YFzXIkdUbxppj4Mc5tiqJ0E5GLwXt3cWH4BlfzQJCO8hYvQOf387MlC3tDPRnmTE6cdRF6vT7VResQrdYIhBC9gL5AphBiHMERQwA5gCXBZVMUJY1FBgK3T8NsbN9N0a8F0naYaLT8X34koNNx4eP3M7VyO5fe/K9UF6lDtdVHcDgwCygBHojYbgOub+2NQoh+wEtAT4L9CbOllA/HOO5A4CGCS0vvkFIeEFfJFUVJKV/E6mG70k9Q10lqAwCGrz5hvpT0HjGGc/9xY0LPlWHU7dZIrF3RVh/Bi8CLQojjpZRvt/Oz/cAVUsqFQohsYIEQ4gsp5fLGA4QQecDjwBFSyk1CiB7tPIeiKEng9Pqb5RPyR9QI2nvj2lzjZG2lvUPKlmjLvvuKqWtX8X8FRdz91OtkWRObg7Mgy0SVzYO/A5fpbEtbTUOnSylfAcqEEP+I3i+lfCDG2xr3bQO2hb63CSFWEGxmWh5x2KnAO1LKTaHjKtt/CYqiJFqt09csEET2C7jjnFTm9mks29pArcPboeVLFM3nY8Mt/8AEjLj2dty5eQk/p0mvozArg4oGd9sHd5C2moayQv9ad+ckQogyYBzwc9SuIYAxNFktG3hYSvlSjPdfQHAYK6WlpbtTFEVRdkGtw0vfvMzwayllkyfWeGsEW+pcaR0EzFs3k7VuNVnr1pCzdBEFP33L4XU1uPML8B1waFLKYNDrKMgwpE8gkFI+Ffr3tl09gRDCCrwNXC6ljF7MxgBMAKYBmcBPQoh5UsrVUeWYDcwGmDhxYhqOMFaUrksLyGZDOn1RzRbx1gjq07hfoP8zjzL4wbvCr2syLew49Chq9j2Q6n0OQEtwk1Ajg05QmGVKyrkaxZuG+l4hRI4QwiiE+J8QokoIcXoc7zMSDAKvSinfiXFIOfCZlNIhpdwBfAuMac8FKIqSWC6f1mz2sD9qsYDoGoE3RuexlDJtA4HO5aT/c49Ts9d+zLnuTkrMmYzp2Zt519zG9qOPx1dQlLSyGPU6zEY9VvOur/HQXvHOIzgs9DR/NLAB2AO4qrU3iOAsi2eBFa30JbwP7CeEMAghLMBewIo4y6QoShK4fcFVyCKTy/n8TWsEXk1DhqYDawHJ0q31zT7H7vGjJbEDtCWWdasp/O5/Tbb1+ugdTPW1fHbAYcx68C70fUq4/4V3yMkrSHr5DPrgKP1k1griDQSNoeko4E0pZfPfcnP7AmcABwshFoW+pgsh/iKE+AuAlHIF8CmwBPgFeEZKubR9l6AoSiI11gbc3pYnkAUCO4eQbq5xUmP3NqtFpEttYMh9tzL24jPIWbIwuEFKSl95moqS/sx66C569+vP/S+8Q2Fxz5SUz6gL3pYLrRlJO2e8dY+PhBArARdwkRCiGGi1J0NK+T07J6C1dtx9wH1xlkNRlCRrXG/A7dfIJbj0YnTTEAQDgV4n2FDtAILJ4yI7mNMhEAivl7z58xCBACOu+xs/v/UFuYsXYl27im/O+gsDf/uFOx9/mdz8wpSVsbFGkJdpRK8XSalFxZuG+lohxL1AvZRSE0I4gD8ntmiKoqSDxif7yNTQ0U1DEFyXIHL8e409KhA4Ux8Icn7/DYPLyeZTz6Hfa8+xx0N3I9auwltQSOZl1/Kw0YROF29DSWI0BgKdTpBvMbHD5kn8OeM5KNTpezowNZRg6RvgyQSWS1GUNBEOBBFNPbFyCzW4/GyucYZfVzs8SCkRQuD1p0c6iYKfv0cKwbpLrkIKQekrzxAAfpx+DIEMc9xt5YlkjAhEhVnJCQTxXvcTBId5Ph76Gh/apihKF+eKEQhiNQ1trnGGF58H8Gs7RwmlQ7MQQMHP32EbPgp/Xj735BeyCggIgetv16S6aADodaJJFtdeueZwDSGR4u0jmCSljBzW+ZUQYnEiCqQoSvoIBGR4KGhk52+spqHIINCo2uElz2JKi0CgczrIXbSAjWdewNP3385/n3+c6oMO5++nn4coHZDq4gE0u+kb9TrKCrMSno4j3hqBJoQY1PhCCDEQSH09T1GUhPL4A+FFYpoEghg1gliq7cFZxPWu1M8mzlv4Czq/j19ycvnv84/zp5mzOPeR57FP2T/VRQszxOif6FdgIcOY2EareGsEVwFzhRDrQ6/LgLMTUiJFUdJG5M0/ENiZbtoXZ7ZRm9uHx6/FXIks2Qp+/p6AwUjxaedyz4gxjN97atotKmOM0Qyk1wkGFlvZWudK2HnjDTM/AE8BAaAm9P1PiSqUoijpITp1ROMM4ugUEy2REjZVO2M2GyVTQ10NzjdeZPvgoQQsWUzY54C0CwIQzDMUS59cM1mmxM00jjcQvAQMAO4AHgUGAi8nqlCKoqQHV9RIn8YO43ibhgDKaxP3JBuPbZs3csvM6Qx12Fk9YI+UlqUthhaW+xRCMLA4K+a+DjlvnMeNlFIOj3g9VwixvMWjFUXpEtxROYRcPo1AQLZrklOyagPmLZvo+9+X0bKsbDrzAgLmTFYuWchNl5zBdJcLHZA1cxbxpEVIFWMLNQKg3SvAtUe8gWChEGKKlHIegBBiL2B+wkqlKErS1Dt95FqMMfdFNw25vFq7agMJFQhg3roZ6+qV9P7gv/T43ydIIdBpGn3encOXs/7Cjf+6heOtVu7p1x//5o00jBqX6lK3KhlDRWOeN87jJgA/CiE2hV6XAquEEL8DUko5OiGlUxQlobz+AKsrbUwqi51czR3VNOT2a3H3DySMpjH8xsvp8cXHGFzBCWy+nDw2nn0xm085G8vG9Qy77Wpm3H4NFXoDGdU78ErJ2n/ciDQlN71zexlTNKs53kBwREJLoShKStQ5vdQ7fdjcPrLNzWsF0TUCt1eLe8RQouzx4F30+eBNth47k7oxE3HssSe2oSMIZFrw+3w8+d+XOe6Z/zLmo7fJqNhG5WFHUzdxb6QheWmdd5XRkMY1AinlxkQXRFGU5KsN5f/ZUudiaK+mgcDj14huBQrWCFIXCHp+/C5lzz/O5pmzWHXTPU322erruPOKC1j407f06N2X3AsvT00hd4O+hc7iREuH1BqKoiSQTwuwqdoZc3ZvrTM40Wt7vbtZp250RzEE5xLYPamZE5C94neG3/wPaidMYfU1tzfZt2n9Gv4680h+nz+PK+98iKNOOiMlZdxdplY6ixMp/etKiqLsstUVNrbUutACkt55ZnIzc8P7fFoAe2gJSr8mqWhw0yciW2j0egKNGtxJDgSaRp93XmePh+7Gl5fPkgefbtLWv3zRfK678BRMGRnc9/zbjBg3Kbnl60AtzSNINFUjUJROyOPXsLlbz9/j9QeaTOaqtHmaPPU31gYabYmaudpSIGjrvB0pZ8lCJp9yJMNvvRLHoCEsfPoNfIXFTY4pKRvE+ClT+c+cTzt1EICW5xEkmgoEitIJubwaFQ2trg3VrClI0ySVtp3vqYtaH6De6QutJ9CYZC52X0D0+sSJYqivY/y5J2LaUcXv9z7BghffxTlwMAAup4OX/nMfXq+HnLx8bnn4WXr0KUlKuRKptXkEiaSahhSlE3J6NSoaPOzRI7vFY2L1CWyrd9M7N9j8U+tonghu8eY6AIwGXXgN4lTp+9YrGJwO5r/8AfahI8Lbt2z8g9suP4eNa1cxbMxEJu13UApL2XF0OtVZrChKOzi9Gi6v1mp654YYTTi1juBawj4t0Gqnr88fCK80lgrC56Pfa89Rs9d+TYLAj199xiUnH86Oiu3c9eRrnSIIxDsjWJ/CldFUIFCUTqgxB1BlC81DUsqYQULKYK2gzukjxQ/8rerx+UeYt29l05kXhLe99+qz3PK3s+hTWsbjb37GxH0PTF0B22FIL2tcxxlTVBsAFQgUpVNyeoNP8xUNsZcxdHi1FvMBbat3UedM/foALZKS0peewlE2iB1TDwlvHrfXfvxp5iweevkDevUtTWEB45dh1NEj29xiCo9IRoOqESiK0g6NWUDdPi3mTb21JiOnR2s2Qiid5P72K7lLF7HpjPP57ZcfefLeW5BS0n+PPbn0pnswZZhTXcS4FWQFh7n2zG67zKkaMQQJDARCiH5CiLlCiOVCiGVCiMtaOXaSEMIvhDghUeVRlK7CG9V+H6tWUO9sfYhnKtv/29L/xSfx5eRy/9ZyrjnvRH797isctoZUFyss0xR/FtBiawYAPXIy2jw2VSOGILE1Aj9wRSh99RTgEiHE8OiDhBB64F/A5wksi6J0GY3NQo0qbe5mI3zSYY3gdpOSQQ//kx5f/h/PZll57tnHOGTGiTw251OsObltvz8J9HrB5AEF5Ge1nbxOp9tZIzAb9eS10TyUqsyjkMBAIKXcJqVcGPreBqwA+sY49G/A20BlosqiKF2Jy9d81bDNNTubevxaoFmwSAemqgr2mzaeHp992Gyf8PsZdsuVDJj9MK9kZXNFbQ1X//NRrr77ETKzErcgS3v1zDZj1OsY2y+vzXb/3Exjk5nCPXOaNg9F1xJirVecLEk5sxCiDBgH/By1vS9wLPBEG++/QAgxXwgxv6qqKmHlVJTOwOltPuN3XZU9fPNvcPvTckRQ6UuzMW/fyuD7bkXniRjtFAgw/B/n0fftV/njgsvZ8uDTPPbWFxw648TUFbYFvXODN3O9TjC2Xx5ZGS1PxSqyNr3RF2dn0Lg65p69shnSs+kckFjrFSdLwgOBEMJK8In/cilldEPfQ8A1UspWpypKKWdLKSdKKScWFxe3dqiidHnRy0dCcBWwFdtsQHo2Cxka6il540Xse+xJ5rYtlMx5MbzPcs9N9P7fp7xz4GGsu+xaJux7IP3ScEnJTJO+SZOQUa9jSM+Wh4YWRgWCxuahwT2t9CuwYDbqyTDuvAWnKs8QJDgQCCGMBIPAq1LKd2IcMhGYI4TYAJwAPC6EOCaRZVKUzi5WjQCCk8W21LnSMhCUzHkBg8PO0n89TvU+B1D21EOI2hq+v/0aJr76LB9mmFl7+nmpLmarGmsDkQqtGTGbiMxGPdYYtYURfXLpX7izqSs3c+d7UzmPIGEpJoQQAngWWCGlfCDWMVLKARHHvwB8JKV8L1FlUpSuoLX2/zUVtiSWJD46l5PSl2ezY/+DsQ8dwdq/38BeJx7GlmMPZGZVJfUZZsT7XzO+X1mqi9qqxtQc0coKs1jsrGuyrSg7dmdy9CzjvEwTlaFRX121RrAvcAZwsBBiUehruhDiL0KIvyTwvIrSZUUPHY3m12TaDQ3t8+4cTDXVbDjvbwDYho9mxZT9ObOqkiFCsOmJVzCnYRAY3z+fXqFaQH6WscVho8XZGWSbdz5TW80GBhXHN5s4skaQylFDCasRSCm/B+K+MinlrESVRVG6ilj9A+nM8sda+j//OHVjJ7GiR29+e/MVpp94OjW33o/v+GmUn3YutXvtl+piNqPXCwqyTBRkmSjJz2y2aE+0AUVZLCmvx5KhZ3xpftxzArLNBnS64II/qVqvGFT2UUVJCy6vhsmgazP7pNOXfsNCoxnqain+5gv6vP0a+QvmETAYePqAQ7nmuIMwGE3sd+h06Nef775ZQiDTkurixpQd0b6fZ2l7zkCPHDPF2W727JWNqR2pInQ6QbbZSL3T1zVrBIqixG9LnRO3L8DIvq1PnGqpozjVhNdL6cuzKf7qM3KXLEAEAjhLB/DrWX/hqoU/883rzzNp/4O5/Jb7yMkrAEjbIADB5p32GtMvb5fOlZsZCgRdsbNYUZT4VTZ4cHo18ixGSvJbvkGma9PQwMfvZ8DTj1A/Ygx/XHg5O/afxraBgznt0InohI6r736EQ2aciBCpu9m1R7a57SRxHSUv08hWvUjpz0YFAkVJMZvbF37SX11hI9tsbNKJGCl6VnE6yFmykLJnH2PL8aey4vYH2Lh2FaWDhpAlBFff/SjDRo8nv6hzzf+JNfQzUXIyjSnNMwQq+6iipFxk0rhAAJaU11ETY/UwAEcri8mkgs7jZsQNl+Hp0ZuFF/6dh269ivP+fAA/zf0MgH0OPrzTBQEhmvYRJFpLcw6SSdUIFCXFItcRhmDuoIUbaynOzmBwTysmvY4Gt59apzfthoYOfPRestav4ZlzLuGqU46koa6W48+6kHF77Z/qojWh1wmKszPYXt/6Os8AFpMBXZLb64uy285OmkgqEChKCtk9fpye2M09VTYPO+zB2kI65g7q9dHb9H/xSd4v6c/5z/2HoaPH88+n5rDHsJGpLloz2WYDe/SwUmXztDkUNHsXOop3V2Ec2UwTSQUCRUmhlpaabJSOAUD4/ZTeczODX3+O2olT+P240/i7z8sRx52CLoVj4VuTk2nEbNRTWmjhjypHeLteLzDqdLgj+l5ykthR3CjedY0TRQUCRUmhSlvspSbTlaFyO73PO5nB61bx9bjJaM+8yRRj8m+c7dV4cy8rzGJrnQuPL4BBLxhXmo/D42f51p35MHdl6Ghnl57hW1G6AafXj92dPp2/pqoKev7fezGrIcLrxXj7NUycNp6ydau4pV8ZK2+4C5lGQcBo0NEr1xxzUl5OZvDmrtcJBhVbMRl0TCwrIDfTSK8cc5Mn8lQ0DaVa97tiRUkT2+LouEwaTWP05eeSt2g+fquV6ohF4y0b1lF26tH0qa/lswwzSy+9lv3OujBt5gRkmvSM6JNDbqYRIQSLNtexI6KmZdALLKadt7reuWYKskzhm79OJygtsLC6wobZqE/5UM5U6H5XrChpIp4RLMlS+vLTwSCQZWXwA3eCptFQV4u9rpbhN1xOvt/Pv486Ht8Pyxkz6y9pEwQABvewkmcxhctUZG3a8ZoTNSdDCNGsTb5vfiZGg65b1gZA1QgUJSXqnN60mSVs+WMtgx65h6qDDmfb0ccx+ooLKb/mEi794WvuHTqcvEW/suzuRxjz55NSXdRm8rOM9IhaAjK4MtjOdNzxdP7qdYKS/NhpprsDVSNQlN3g9be6uF6LttalSW1A0xh+4+UEMswsvfGfvOp285vRyMGfvMdBe+zJOUsXUz1lf7al4bKRQsDgqOUeITgCJ/LJvrF/oC398i3ktTCju6tTgUBRdsPGakfM7a2NVQ8EZLNJZKmQ+9uvjPnbWeQtms+q6+/k8Zdnc+8Nl/FA336UAv/duA69prHy5nshjZqCGvXOzWzxaT9ygla8w0FNBl2z5SW7C9U0pCi7SErJljpXeP3ZSGsr7fTNz4yZOqDK7knpDOHc335l8P23kbdoPi5LFgvOvIDao4/nqNHjGTpqHPsf9id2XHQaRd/PZc3fb8DVf0DbH5pker1gUI+sFvcXWTP4o8pBhlGX8jH6nYGqESjKLrJ5/Pg1GTMvUKXNzYptDcgYQzFTOVqoz1uvMmHWceg2beDufmUUOR3c5HCAEJT0H8gBR/wZnU7HylvuZc3lN7DprPRcTHBwDysZhpZv8LmZRjKMupRMDuuMVI1AUXZRbSgA1Di89Mnb2dHY4Pbh8QXw+AKU1wZrDI28/gDV9sRNIjNVVZC9chm+3Dx8+YX4LZZws86AJx+k9NVn+TEnj6NqdqDT6znnhrs58oTTmn2Ou08/Np7/t4SVM1rjKl3xKLCaWk3V3agwK6PF5SWVplQgUJQ4ePxasyfQWqcPoFmNoCpiDPvaKjvF2RmY9Dq2N7jZWO1MXNoIKRl7yRnkLFvS4iEfjhzLudu2cMo1t3PUSWeQYU6PkTJ98yxUNLjb7Hw36AXDe+fE9ZlF2Sb0adi3kY5UIFCUNji9fpaU1zNlYGF4m5SSOmcwAHj9AWxuX3gxk8jJTJomWby5Dq8WrCEkUo8vPiZn2RLWXn499sFDMdTVsnnxAuZ//xWT959G8bQjCYyZwPMGQ1IDgEEv2uwTybMY0esEG3Y07XwvtJrIt5ioaHBjc/sZ0jM77jb/wqwMAumYrCkNqUCgdBpOrx8tIJO6ehTAyu027G4/O+ye0Bj1nf0DjWodwUDg9mnYotJGRL9OCE1j0KP/wj5wMOvOvJBvvviYN158ivWrl9OzTwlyv4PYZ98Dabl7NXEGFVv5Y4ej1af93MzgYjwbqx3hGlPj8FBrhoGyoixcXq1dTT16nUCPqhHEQwUCpdNo7GRNZiDYXu+mxh588t9c4wwHgjqHr8lx1Q4PpYWWcNroROs75wUMDjubzvoL0mCg94dvkbV+DYsfeparLjiZ3+fPo3TgYK6662EOPuo4DCnMCVRoNRGQkjUV9pj7I0f2FFkzwk1rPXPMTUZdqfb+xFGBQEkZKWW7UhVU1LsRIpg0LBl8WoDVFTtnqFbbvTi9fiwmA7XOpv0CdU4fgYBs0j+QKD0+fZ9hd1wLQPFXn/HTlTfT856bqR82iqpDpvPngMYJZ/2FKQcelvK00BlGHRaTgZJ8PRurnTFrBZHLcpbkZ1Jl86DTkbTfs5LA4aNCiH5CiLlCiOVCiGVCiMtiHHOaEGKJEOJ3IcSPQogxiSqPkn4ab57xHevF6dVwePw4vcnJ2Lm+qnlzRnmtC6BZINACkmqHt9n2jmZduYwRN15O3bjJzL3sOoxLF3Ho6X+iyFbPp4f/CYTggMNnsM/BR6Q8CADkW4J5f/Q6Qf/C2CN9IgNBoTUDi0lPn7xMVQNIokTWCPzAFVLKhUKIbGCBEOILKeXyiGP+AA6QUtYKIY4EZgN7JbBMShqxuf3ohCDX0nazReTY+4oGDwOKmv7ptrd20RavP8DWOlez7VvrXPTIzojZ+bm20h73EMi4SEn+zz9grKvGMWgIvtx8xlw6C681hxOF4MuH/8mepgzeLO5JzqAhFJ6XvOGe8cqzRD7tW9hQ7cQXFVxzo9I6lBZawk1wSnIkLBBIKbcB20Lf24QQK4C+wPKIY36MeMs8oCRR5VHSj93jRwjaDASBgKQiYiWvKpuHAUU7uz132D1srHYwuiSv3SmEnV4/Xn+APEvTjJXltc6YaSL8mmTldluz7dCxC8tnbtrAkHtupPibL5ts10wZzH/xPRqeuJ9zLr+eo048g4q8fCo67MwdKz/i56rXCfoXWFhbubOvQIjmfT7xzBFQOlZS+giEEGXAOODnVg47F/ikhfdfAFwAUFpa2tHFU1LE4fUTkJJ+bRy3w9E0JUODy4fbp2E26vFrAVZsa8DjC/DrHzWMLc1rkns+lkBAUu/ysanGyQ67B6Nex5SBhZgMuvD+zbXNawONErmYjLGuhtIXn6L0hSeRBgM/nHkBr5Vvwvbd/xgm4YBHnsM5ehx3PfFqwsrQFqvZENfPwGTQkRWVYqNfgYVNNTv7CqwZhpgLySjJlfBAIISwAm8Dl0spG1o45iCCgWC/WPullLMJNhsxceJENTC4i7B7/Pi0tttSYuXtr2wIjtJZXWEPj893ejV++aOG0SV5FEQtBu72aayttGNz+3H5/E2acLz+YKfwyL65AGytdzVrvki0jIptlL7wJH3ffBmDy8nKvQ/ggpodfPfSbMyZmUw7ZiYFM2fhHDoiqeWKZjLoGFicxZLN9W0emxejpqfXCcoKs8Kd8PE0CyqJl9BAIIQwEgwCr0op32nhmNHAM8CRUsrqRJZHSR9un4amSVwBDb8WwNBCk47bp8Uckllpc5OVoW/Wju/XJL9tqmVwj2xKQ52T9S4fS8rrWp3Qtb3eTe9cM4XWDDbVOHfjytqvx6fvM/ymf6Bzu1h/4GFUXXYdi30+Nl97CZdcfxeH/OkErDm5SS1TS/rkmZs090QSoukqly0dV5KfycYaBx5foFn/gJIaCQsEIthz9yywQkr5QAvHlALvAGdIKVcnqixK+rGH2tOlDHYa52fFvmm01AFb7/KxfFvMCiZSwuoKGw1uH4VWEyu32VpNC91oxTYbe/SQOD3JWTBGeL2U3XU9g956hQWZFk4MBBick8tVe+zJHsDs975Oq5XAIJgKwqgPNvlE94kM651Dea2LBldwjkWsGgEEl4YcUJTFym02FQjSRCJrBPsCZwC/CyEWhbZdD5QCSCmfBG4GCoHHQ3/wfinlxASWSUkTkTeRBrcvZiCod/laXM5RStpM2bC93t2u5SDdPo3l29pu8thVwudjwJMPYF27GoOtHu+yxfSy2/g3MLtfGdNPOJVpR5+w8/g0CwIFVlN4SGduprHJ71AIKM7OoMiawYKNtXj8WqsT//rmZbK1zt1mf46SHIkcNfQ9tD6/W0p5HnBeosqgpK/ItAsNrtgdj2sqYo/OSaQOHf4Z9cH9L53FwG//hz00FHRbUQ+eHzeJ/Muu4z/DRiXoxB2nb0SG1TyLsUmzXG6mMTxia1xpXpvNa0IIRvVNj+YuRc0sVlIkukYQbXu9mzpn8+2dgamqkpHXXEz9mIksPPxoPvvhaw549j8cUl/L9Tod/e99gkGhTt9JCS5LdLv9rjIZdBRHjO2PbvaJXNnLbNQzJMYSktHUhLH0oQKBknRSSpwRC7e7vBpefyA8fFMLyCZjzTsT4fMx6orzyV28kPxffqDH7IfIBQ4Bvth7KuPvfZy8gqKklKVxhbRVLcx7iCZE8CtWrahPnhldxDBPi8mA0aALj64qtMbu41E6BxUIlKRz+bRmnbcNbl94NumyrfW4fcnpsO0oDXW1/PjVJ4x7/AHyt5Xz+71PULfHnuhvuIyDVy2n/MTTEDf9i7wktfvnZ5kY2isbIQTVDm+T1NiRMk16RpfkYjbqMep1lNc6WbmtaeDQ60TMSV65mUZ22DyYDGolsM5OBQIl6ewxZuA2uIKB4I8dDiobkpPBMx6WP9ZSMudFchf9yvpLrqR66iHhfXm//kjhrVexxOng1R2V5AQCnAbMnTAF7ahjgwe99QXfOOxolqykLQBvCd3cGzubh/fO4Wd3dbPOdYNeMLZfXpNJX33zMtlU7WxSY+tXkBlzDYC8UCCInrOhdD4qEChJ54gxPLPB7afK5mFdmjQJWdatZs+7b6Bw3ncEDEY8xT0Ye/EZ/Hb+35hTOpBz3C6G3HMT1XoDB3o9/DnUEF87bhKBZ/7b5LO0rORl0bSaDYzqm9sk1YbJoGN47xx+21QX3qbT0SwIQLATd1APK7+XB0dPGfSC/oWxVzFo7CdQeYE6PxUIlN0SCMgmbcfxiJWTp86Z+Myd8Sr48RtG/eN8pMHAmsuu44eRY/n+lx84fs4LHDn7EQzAUKDqgEP49Zb7MOUXkr/oV/IW/sKWE05HmpL/hJxp0jOwOIteOeaYw04LrRnsu0cRNo8Pu9tPttnYLL9So545ZjaYHdjcfvoXZrWYvynHbESnQ9UIugAVCJRdZvf4WbK5jvH98+NePrDxfdHaWsowKaSk7xsvsufdN2AfsAdLnniFTxcv4K7zT0YIwS9jJ+HLzOToed/zx3l/Zd1fr8akD1537eR9qZ28b1KLazToKLKaguP3szLaDMiZJj2ZJj092h7Qwx49rCzb2kBpQcsJ4HQ6Qd88S7iTX+m8VCBQdonN7WPhpjp8/gCrttsY0y8vrvcFRwwlZz2BWEyV2xly320UffsltiHDqR8zHnfPPhjnfUfBgnnk2xr4WOj48ugTOKpPP8ZbrFx550NMnjqN/MJiAL72uAlkmFN2DQDZZgOTBxQkbNJZoTWD0SW5bSaEa2mNAaVzUYFAaTe7xx8OAhBMCx25nm9L3D4ttBhNMkrZlM7poOTNVxj4n/sQPh8Vhx2NZfMGSl95Fp3PyzbgM2B+QREbjzyGg/aZCkBOXj6HHzuzyWelOggA9C/MSvjM45aajiK1pyaopC8VCJR2kVKydEt9s+ycq7bbKBhoCjdP2Nw+6pw+bG4/do8fh8cfV76fjqJ3OCia+xmFP8wlZ9lisv5YiwgEWFQ6gGuzrFTV7OBfr32M8Hp46583YRw8lEn7T+OIfv2TVsZdZTbq6ZmjOmiVjqMCgdIum2tcMXPRu7waqyttGPU6KhrcSUvcBqC328hdvABjbQ2muhpyliyk+H+fYHC78BQWsSzLyic6Hf8LBPhm8waGjBzL5HHBOb3SlMHxt9ybtLK2pKQgE7vb32w29aAeVpxeP9vqduZMKi2wpF0eIqVzU4FAiZvbp7FuR8vDO8trWl7MZXcJr4dht1+Dob6OVdffhad3XwCyly9h9KVnk7ltS/jYeoOBF4VgyBOv4tnvID59bw5rVyxl6l778bdJ+5Cdm5ewcu4KIaB/QRZCwLz11eGO8/wsE2WFFqQMLoZjc/sx6AV98zPb+ERFaR8VCLoxnxbAr0n0OtFs5IeUkl831NIjO4P+hcEn0LWVdrQUjO7RO+yMvvRsCud9h5ZhZsqxB/H75dcj9HpG//Mm7JkWTjEYWeH3sQPI7lfG6Mn7curQERTpdBxx3KlJL3NLcjKN2D1N+0l6ZJvDeXeG985hSXk9RoOOEX1yEEIgBIzpl8fPf9TQNy9TreildDgVCLqhGoeXJeV14SfPTJOevQYUNFkcZlONkwaXjwaXj+0NbvrkZrYrpfMu0TRMNTvIqKpA+P1o5uCT7/AbLiN75TJeP/5UvnQ4uPDrz5l8x7XBa5m0D59efh36Tz/klIlTGDl+ctJy+eyKYb2z2WH3Npk4Vxox8qZHjpm++V4KraYmHbFmo55RfXPJylCds0rHU4Ggm9la52Ll9oYmT6Qur8aKbTZGlQTTArt9Gut3OML77W4/q90dmxJa73BQ9syjZC9bTMaOSkzVVZhqdiBiDCnyZ5j5U0Dj/95+DWtOLssmTuGvOj3DSvpTd81t9DAYuGhsx+TxzDDq2lznoC25FiM9sjNYW2lvkvmzODuDbLMRa4aBaruHOqePPIux2eIsjTmCoqmJW0qiqEDQCdU6vJTXusgw6uJK9wvBm/vmGicbq2Pnia9ocFNQZ6JvXiZrKhLbBFTw07cMu+UKzFvLaRg+GnevPvxmMLAuy8ry6irW2m14gEnjJnHSiWdQP3IsY1Ys5fjhoygpG4ROF6y51HZwuaxmA5PLCvAFAtQ6fGyrd1Ftb3u2sxDBWlVuppGSPEt4Hd6ApMmT/8DirNDxghF9cpn3R3WT2sDOz1NNP0pyqUDQSTS4fVTbvWyvdzdJ0WDU6xhQFDsXjNcfYHOtkyqbJ+ZIn2irt9sIBCQVDR3bBCR8PjJXLiPwzefkfv81o35fyEZzJvdOmMKxL74LwMUnHoZHr2fPg49g8IgxDBkxhj2GjWRbqHno4EFDOrRM0fS64EIpOp0gQ6enV66eXrlm1lXZ+aPK0eTYAquJHLMBi8lAlsmA1WyI2W5fVmihxuGh1uEL1wYaZZr0jO6bq57ylbSgAkGa21rnYl2VvcXminWVdoz6pmmCAwHJphonG6od7UrdoAVk3LnrWyKlpGprORVLFlD44zcc5nRQ9N1XGBzBJ+Mq4F6h46X+AxkydGT4fY++/n/oDan7cxzaO7tZAjaAQcVWskwGVmxvoNiaQVlRFtYYx8USfvJfX82A4ubBulAla1PShAoECVJe68RiMuzyE5/XH2DFtgaqWsgjH2nVdhu1Dh8BKfEHgikcdredu0VSkrlpA8atm3GtWUnPgEbWls045/+EaeN6cjwe8oDGLk13YRHbjzyGX/PyWF1QRP5e+zNq4GAeikrMlqogoNNB79xMeue2PCSzV66Z4uyMXRqtYzbqmVhWEHfwUJRUUH+dCeD0+llTYUciGdk3lx7ZwZQEgYCkvNaFPxCg0JpBjtmAEAK3T8Ph8eP0anj8ATx+jWq7F68/vpu5lOxyc47w+chauwpX6QC0rOBTa8a2LfR95zWKvv4Cb14+7pL+VEhJ4OfvGLxtC3m+ppOe/NZstufksjzLiq90AIaefTCXlCIPPAz/vgeCTkc2MGGXSrjrCqwmBhZlodcJ9DqBTgikhICUCBFsVmsps2a03RmyqYKAku7UX2gHk1KybGtDOJ3C7+X1DOsd/H59lSO88tb6KgfG0Nj96HQNiaB3ODDtqMRbUIhmzUZvt9H3rVfo98ozZG7fihSCioIidhiMDKvYigB+0ukYVDqAHst/p6SuhjXA59k5rB80BN+AwZgH78nQI48hs18ZCIEOaGzsSOZqw0XZGRRnZ7C1zkW904dOF2zSaSmPvqIoTalA0ME21Tipj0gTICUs39oQ89hEBoCMreUUf/YBBT99S/Yfa8ncWh7e59HpkFJilpJ1AwfjvvNB7L/Np+LtVxkgBE/kFzJ3z+EweBjTTziN/nvsieawgzmTAr2egoSVupXrMeqwZhjQApJ6lw8pg6N1BhZbw53lffMysXv8BKRUSycqSjsIKdMgD3w7TJw4Uc6fP3+X3rup2onJoKNX7s7skTa3jx12L33yzGQYmk7Wcfs03D4Nb2gGrhCgFwKdTmAINTcYdDr8geB+jz/A8m31u5ddU0oyKrdjrKuFQAAhAxga6sks30Tmlk0Y62vRELh9XvwN9WTbGsizNyCcTsob6qkRYHU4GOMKDhMtLyhCv/dUqnr15oFn/0Mvg4EBWVayLVn8PGIMA868kJET9sLn9WJvqCevsCiBqY1NDO6ZTb3Lx5oKW1wd2X3zM9mjh7VJE47XH2CH3UOmUU++GnWjKHERQiyQUk6MuS9RgUAI0Q94CegJSGC2lPLhqGME8DAwHXACs6SUC1v73N0JBCu3N1Be48JqNtA3L5NKm5taR/DpXa8X9C+w0L8wixqHl/JaJzUOL7vy49HbGiiY9x1CSty9++Lu1QdTTTXWVcuwrl6O3uVCy8wkYM4ErxdTQz3GuhpYt5q8LZvJ8MRu7/cBdUIgpMQQem3PzSd71Fj8mRZWzP2cIr0OYTTxS2kZS0dNoO+RMxgzaR8CgQD2hjqyc/MTcqPX6wT9Cy04vRqVNnc4GAoRXOS8rCirSZpqrz/AmkpbMH+OLhhcbW5/k1rSgOIsBhUnb5lHRenKWgsEiWwa8gNXSCkXCiGygQVCiC+klMsjjjkSGBz62gt4IvRvx1u+nMLnX0EYs/Dl5FJtzUFkZJCXkUHAaEJoGjWan/qAhvRriIBGoV/DYKvHWF+Hsb6WgK0B6XCA140r04LNmoM7r4CePXthcNipXLaYXosXMGjDOvQtVAs8QmAXOswBjSzALwRaQRG+3DwWbNnEYo+HlcB2AJ2OkeP34s8XXYGrXyn3zn4EYTSSX1hMfmExeYVFlA4cTL8Be4Q/vyL0b2noq5FOpyMnLzGNOnkWI8P75GAxBf+cPH4r2+vdZBj0FGSZYq5gZTLoGNEnt8m2QEBSZfdQXuuiR3YG/VpZHUtRlI6TtKYhIcT7wGNSyi8itj0FfC2lfD30ehVwoJRyW0ufs8s1gjfegJkz2z6uFW6dDmcggA/IY2fHaKTlwAfAR4BDCPYZMIjLTz8Xf04ej33+MQvsDVhy87BYc8i2Wum/x1COPCGYFG3Zb79gMprIzc8nN6+ALGs2hEa6SIK/p8a28UiN2ww6HUZ98OlaC0j8mkQiMYVGxxj1utAIGtDrdGQYdJiNekyGYJ9BQAY7u32aRAtIfFoAp1fD4fXj9mpkGHVkZRiwGA3odKATwSay4uwMNRtWUdJcSpqGogpQBnwLjJRSNkRs/wi4R0r5fej1/4BrpJTzo95/AXABQGlp6YSNGzfuUjmemz2bHz7+mCyfjyyfDzNg1um45vLLwWDgky++YNnKlehMJnQGA8JoxFhUxMU33QS5uXz0ySds2rSJjIwMMs1mcjSNIk1jyv77Q3Y2m+vq0DIysFqtZGVlYTbHXkhcURQl2VIaCIQQVuAb4C4p5TtR++IKBJF2p49AURSlu2otEMQ3m2bXT2wE3gZejQ4CIVuAfhGvS0LbFEVRlCRJWCAIjQh6FlghpXyghcM+AM4UQVOA+tb6BxRFUZSOl8hRQ/sCZwC/CyEWhbZdT2gwi5TySeD/CA4dXUtw+OjZCSyPoiiKEkPCAkGo3b/VnlIZ7KC4JFFlUBRFUdqW0D4CRVEUJf2pQKAoitLNqUCgKIrSzalAoCiK0s11uuyjQogqYNemFkMRsKMDi9MZqGvuHtQ1dw+7c839pZTFsXZ0ukCwO4QQ81uaWddVqWvuHtQ1dw+JumbVNKQoitLNqUCgKIrSzXW3QDA71QVIAXXN3YO65u4hIdfcrfoIFEVRlOa6W41AURRFiaICgaIoSjfXJQOBEOIIIcQqIcRaIcS1MfZnCCHeCO3/ObSCWqcWxzX/QwixXAixRAjxPyFE/1SUsyO1dc0Rxx0vhJBCiE4/1DCeaxZCnBT6XS8TQryW7DJ2tDj+tkuFEHOFEL+F/r6np6KcHUUI8ZwQolIIsbSF/UII8Ujo57FECDF+t08qpexSX4AeWAcMBEzAYmB41DEXA0+Gvp8JvJHqcifhmg8CLKHvL+oO1xw6LpvgMqnzgImpLncSfs+Dgd+A/NDrHqkudxKueTZwUej74cCGVJd7N695KjAeWNrC/unAJwSzO08Bft7dc3bFGsFkYK2Ucr2U0gvMAf4cdcyfgRdD378FTBOde3HhNq9ZSjlXSukMvZxHcDW4ziye3zPAHcC/AHcyC5cg8Vzz+cB/pJS1AFLKyiSXsaPFc80SyAl9nwtsTWL5OpyU8lugppVD/gy8JIPmAXlCiN67c86uGAj6ApsjXpeHtsU8RkrpB+qBwqSULjHiueZI5xJ8oujM2rzmUJW5n5Ty42QWLIHi+T0PAYYIIX4QQswTQhyRtNIlRjzXfCtwuhCinOBiV39LTtFSpr3/39uUyBXKlDQkhDgdmAgckOqyJJIQQgc8AMxKcVGSzUCweehAgrW+b4UQo6SUdaksVIKdArwgpfy3EGJv4GUhxEgpZSDVBessumKNYAvQL+J1SWhbzGOEEAaC1cnqpJQuMeK5ZoQQhwA3ADOklJ4klS1R2rrmbGAk8LUQYgPBttQPOnmHcTy/53LgAymlT0r5B7CaYGDorOK55nOB/wJIKX8CzASTs3VVcf1/b4+uGAh+BQYLIQYIIUwEO4M/iDrmA+Cs0PcnAF/JUC9MJ9XmNQshxgFPEQwCnb3dGNq4ZillvZSySEpZJqUsI9gvMkNKOT81xe0Q8fxtv0ewNoAQoohgU9H6JJaxo8VzzZuAaQBCiGEEA0FVUkuZXB8AZ4ZGD00B6qWU23bnA7tc05CU0i+E+CvwGcERB89JKZcJIW4H5kspPwCeJVh9XEuwU2Zm6kq8++K85vsAK/BmqF98k5RyRsoKvZvivOYuJc5r/gw4TAixHNCAq6SUnba2G+c1XwE8LYT4O8GO41md+cFOCPE6wWBeFOr3uAUwAkgpnyTYDzIdWAs4gbN3+5yd+OelKIqidICu2DSkKIqitIMKBIqiKN2cCgSKoijdnAoEiqIo3ZwKBIqiKN2cCgSK0gYhRJ4Q4uLQ932EEG+lukyK0pHU8FFFaUMoTflHUsqRqS6LoiRCl5tQpigJcA8wSAixCFgDDJNSjhRCzAKOAbIIpnG4n2Cq5DMADzBdSlkjhBgE/AcoJjgB6Hwp5cpkX4SitEQ1DSlK264F1kkpxwJXRe0bCRwHTALuApxSynHAT8CZoWNmA3+TUk4ArgQeT0ahFSVeqkagKLtnrpTSBtiEEPXAh6HtvwOjhRBWYB92pvYAyEh+MRWlZSoQKMruicziGoh4HSD4/0sH1IVqE4qSllTTkKK0zUYwrXW7SSkbgD+EECdCeL3ZMR1ZOEXZXSoQKEobQtk7fwgtJn7fLnzEacC5QojFwDJiL6mpKCmjho8qiqJ0c6pGoCiK0s2pQKAoitLNqUCgKIrSzalAoCiK0s2pQKAoitLNqUCgKIrSzalAoCiK0s39PylK0+yCY95wAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# remove the 8th and 23rd sample rows from the data\n",
    "exclude_data = np.delete(random_response, [7, 22], axis = 0)\n",
    "\n",
    "# store the number of samples left over after removing the bad trainee data\n",
    "# this should be equal to the number of rows left over in the data\n",
    "number_of_samples = exclude_data.shape[0]\n",
    "print('The number of samples left after removing the bad trainee data is', number_of_samples)\n",
    "\n",
    "# we need to find where the remaining data exceeds 3.25, and then delete those rows\n",
    "# we will then re calculate the number of samples to get an idea for how many samples were removed by this filter\n",
    "too_big_indices = np.where(exclude_data > 3.25)\n",
    "exclude_data = np.delete(exclude_data, too_big_indices[0], axis = 0)\n",
    "number_of_samples = exclude_data.shape[0]\n",
    "print('The number of samples left after removing the data > 3.25 is', number_of_samples)\n",
    "\n",
    "# this replots the data to show the results\n",
    "# any large peaks should be gone now\n",
    "r_max = np.max(exclude_data, axis = 0)\n",
    "r_min = np.min(exclude_data, axis = 0)\n",
    "r_avg = np.mean(exclude_data, axis = 0)\n",
    "fig = plt.plot(t, expected_response, 'k--', t, r_avg, 'r')\n",
    "plt.fill_between(t, r_min, r_max, alpha = 0.3)\n",
    "plt.legend(['Expected', 'Filtered Mean', 'Min/Max Range'],loc = 'upper left')\n",
    "plt.xlabel('time')\n",
    "plt.ylabel('position');\n",
    "\n",
    "# in this final step we will at the time domain on top of rest of the data set \n",
    "# this will allow us to save the data in one large array rather than an array for the time\n",
    "# and an array for the pre-filtered data\n",
    "\n",
    "# a subtlety of NumPy is that ranges have dimension 1, while our response has two dimensions\n",
    "# in order for the append function to work we have to make our time domain two dimensional\n",
    "# which is just 1 x number of time points\n",
    "# we use reshape here instead of resize because resize will change the original t data\n",
    "# while reshape will just output a new array of the correct shape\n",
    "t_vec = t.reshape([1, len(t)]) \n",
    "\n",
    "# now we can use append to add the sampled data to the end of the time domain\n",
    "final_data = np.append(t_vec, random_response, axis = 0)\n",
    "\n",
    "# this step we realize it might be nice to include the expected response in the output as well\n",
    "# so it does not have to be generated again later\n",
    "final_data = np.insert(final_data, 1, expected_response, axis = 0)\n",
    "          \n",
    "# so the rows should look like 1: time, 2: expected response, 3-52: sampled data\n",
    "# step just checks that our final data has the correct shape, 52 x number of time steps\n",
    "print('our final_data for output has shape', final_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Homework"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Here is a website that goes through a good sample of NumPy [problems](https://www.w3schools.com/python/numpy/numpy_exercises.asp). \n",
    "It covers some functions that this lesson did not explicitly cover as this lesson tried to introduce some more advanced attributes of NumPy that might be immediately relevant for your data processing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "104px",
    "width": "251px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}