12-Fast-Fourier-Transform.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.01 Discrete Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Represent a continous function as linear combination of sines and cosines.\n",
    "* The sine and cosine components are referred to as **frequencies**.\n",
    "* For some problems, frequency domain is more efficient than time or space domain.\n",
    "\n",
    "**Euler's formula**\n",
    "\n",
    "aka **complex exponential**\n",
    "$$\n",
    "e^{i \\theta} = \\cos \\theta + i \\sin \\theta\n",
    "$$\n",
    "where\n",
    "* $i = \\sqrt{-1}$ is complex\n",
    "\n",
    "**Roots of Unity**\n",
    "\n",
    "aka **twiddle factors**\n",
    "\n",
    "For a given integer $n$ the primitive nth root of unity is given by:\n",
    "$$\n",
    "\\omega_n = \\cos(2\\pi/n) - i \\sin(2\\pi/n) = e^{-2\\pi i/n}\n",
    "$$\n",
    "\n",
    "##### Example: $w_4^k$ \n",
    "$k = [0, 1, 2, 3]$ represents $\\pi/2$ steps counterclockwise around the real-complex plane.\n",
    "* $w_4^0 = 1$\n",
    "* $w_4^1 = -i$\n",
    "* $w_4^2 = -1$\n",
    "* $w_4^3 = i$\n",
    "\n",
    "#### Discrete Fourier Transform (DFT)\n",
    "DFT $y$ of the sequence $x = [x_0, \\cdots, x_{n-1}]^T$ is given by:\n",
    "$$\n",
    "y_m = \\sum_{k=0}^{n-1} {x_k \\omega_n^{mk}}, \\qquad m = 0, \\cdots, n-1\n",
    "$$\n",
    "where\n",
    "* $\\omega_n^{mk}$ is the kth element of the nth-root of unity\n",
    "\n",
    "Expressed in matrix notation:\n",
    "$$\n",
    "y = F_n x\n",
    "$$\n",
    "where\n",
    "* $F_n$ is the Fourier matrix with entries $F_n = \\omega_n^{mk}$.\n",
    "\n",
    "#### Inverse DFT\n",
    "Inverse DFT $x$ of the sequence $y = [y_0, \\cdots, y_{n-1}]^T$ is given by:\n",
    "$$\n",
    "x_k = \\frac{1}{n} \\sum_{m=0}^{n-1} {y_m \\omega_n^{-mk}}, \\qquad k = 0, \\cdots, n-1\n",
    "$$\n",
    "\n",
    "Expressed in matrix notation:\n",
    "$$\n",
    "x = F_n^{-1} y\n",
    "$$\n",
    "where\n",
    "* $F_n^{-1} = (1/n) F_n^H$ is the inverse of the Fourier matrix\n",
    "  * $F_n^H$ is the conjugate transpose of $F_n$\n",
    "    * Note: $F_n$ is almost but **not** Hermitian\n",
    "\n",
    "\n",
    "#### Significance\n",
    "DFT gives a trigonometric interpolant using only matrix-vector multiplication with $O(n^2)$ work.\n",
    "\n",
    "DFT of purely real sequence is in general complex.\n",
    "\n",
    "$y_0 = \\sum_{k=0}^{n-1} x_k$ is called zero frequency or **DC** component\n",
    "\n",
    "$y_{n/2}$ is highest frequency representable aka **Nyquist frequency**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate computing the roots of unity for $n = 2, 4, 8$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "omega2:\n",
      " [ 1.+0.0000000e+00j -1.-1.2246468e-16j]\n",
      "omega4:\n",
      " [ 1.0000000e+00+0.0000000e+00j  6.1232340e-17-1.0000000e+00j\n",
      " -1.0000000e+00-1.2246468e-16j -1.8369702e-16+1.0000000e+00j]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "def roots_unity(n):\n",
    "    \"\"\"\n",
    "    Return the nth roots of unity $e^{(-2k \\pi i)/n}$.\n",
    "    \"\"\"\n",
    "    k = np.arange(0, n, dtype='d')\n",
    "    return np.exp((-2.*k*np.pi*1j)/n)\n",
    "\n",
    "\n",
    "omega2 = roots_unity(2)\n",
    "print(\"omega2:\\n\", omega2)\n",
    "np.testing.assert_almost_equal(omega2, np.array([1, -1]))\n",
    "\n",
    "omega4 = roots_unity(4)\n",
    "print(\"omega4:\\n\", omega4)\n",
    "np.testing.assert_almost_equal(omega4, np.array([1, -1j, -1, 1j]))\n",
    "\n",
    "omega8 = roots_unity(8)\n",
    "# Compare to $\\omega_n = \\cos(2\\pi/n) - i \\sin(2\\pi/n)$.\n",
    "np.testing.assert_almost_equal(omega8, np.array(\n",
    "    [1, \n",
    "     np.cos(1.*np.pi/4.) - 1j*np.sin(1.*np.pi/4.),  # 45d\n",
    "     -1j,\n",
    "     np.cos(3.*np.pi/4.) - 1j*np.sin(3.*np.pi/4.),  # 135d\n",
    "     -1,\n",
    "     np.cos(5.*np.pi/4.) - 1j*np.sin(5.*np.pi/4.),  # 225d\n",
    "     1j,\n",
    "     np.cos(7.*np.pi/4.) - 1j*np.sin(7.*np.pi/4.)]))  # 315d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate computing the Fourier matrix for nth roots of unity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def fourier_matrix(n):\n",
    "    \"\"\"\n",
    "    Return the Fourier matrix for nth roots of unity $F_n$.\n",
    "    \"\"\"\n",
    "    # Compute the nth roots of unity.\n",
    "    omega = np.exp((-2.*np.pi*1j)/n)\n",
    "    # Entries of Fourier matrix given by a_{ij} = \\omega^{i*j}.\n",
    "    i, j = np.meshgrid(np.arange(n), np.arange(n), indexing='ij')\n",
    "    return np.power(omega, i*j)\n",
    "\n",
    "\n",
    "F4 = fourier_matrix(4)\n",
    "np.testing.assert_almost_equal(F4, np.array([[1, 1, 1, 1],\n",
    "                                             [1, -1j, -1, 1j],\n",
    "                                             [1, -1, 1, -1],\n",
    "                                             [1, 1j, -1, -1j]]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the DFT for the sequence $x = [4,0,3,6,2,9,6,5]$ from the Fourier matrix with $n=8$ and verify that the inverse DFT returns the same sequence. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "y:\n",
      " [35.        +0.00000000e+00j -5.07106781+8.65685425e+00j\n",
      " -3.        +2.00000000e+00j  9.07106781+2.65685425e+00j\n",
      " -5.        -2.13162821e-14j  9.07106781-2.65685425e+00j\n",
      " -3.        -2.00000000e+00j -5.07106781-8.65685425e+00j]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "n = 8  # 8th roots of unity.\n",
    "\n",
    "# Compute DFT as y = F_n x.\n",
    "F8 = fourier_matrix(n)\n",
    "x = np.array([4,0,3,6,2,9,6,5], dtype=np.float64)\n",
    "y = np.matmul(F8, x)\n",
    "\n",
    "print(\"y:\\n\", y)\n",
    "np.testing.assert_almost_equal(y, np.array(\n",
    "    [35, \n",
    "     -5.07+8.66j, \n",
    "     -3+2j, \n",
    "     9.07+2.66j, \n",
    "     -5, \n",
    "     9.07-2.66j, \n",
    "     -3-2j, \n",
    "     -5.07-8.66j]), decimal=2)\n",
    "\n",
    "# Confirm that DC component equals the sum of the sequence.\n",
    "np.testing.assert_equal(y[0], np.sum(x))\n",
    "\n",
    "# Compute the inverse F_n^{-1} from complex conjugate transpose.\n",
    "F8inv = (1./n)*F8.conj().T\n",
    "\n",
    "# Verify the inverse DFT is $F_n^{-1} y = x$.\n",
    "xx = np.matmul(F8inv, y)\n",
    "np.testing.assert_almost_equal(xx, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.02 Fast Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "FFT uses symmetry and redundancy to compute the DFT more efficiently.\n",
    "* DFT of n-point sequence can be computed via 2 DFT of $n/2$ (half) length.\n",
    "\n",
    "A pattern emerges from the Fourier matrices.\n",
    "* The matrix $F_4$ is composed of diagonally scaled blocks of $F_2$.\n",
    "* The permutation matrix $P_n$ is used to rearrange the matrix so that even numbered columns appear before the odd numbered columns.\n",
    "$$\n",
    "F_1 = 1, \\qquad\n",
    "F_2 = \n",
    "\\begin{bmatrix}\n",
    "1 & 1 \\\\\n",
    "1 & -1 \\\\\n",
    "\\end{bmatrix}, \\qquad\n",
    "F_4 = \n",
    "\\begin{bmatrix}\n",
    "1 & 1 & 1 & 1 \\\\\n",
    "1 & -i & -1 & i \\\\\n",
    "1 & -1 & 1 & -1 \\\\\n",
    "1 & i & -1 & -i\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "#### FFT\n",
    "Recursive divide-and-conquer algorithm for computing DFT.\n",
    "* Work required $O(n \\log_2 n)$.\n",
    "  * Compare to $O(n^2)$ for matrix-vector product form of DFT.\n",
    "* Can be implemented in-place and using no additional storage.\n",
    "  * Ignore implementation below.\n",
    "* FFT can be reworked to compute inverse DFT.\n",
    "\n",
    "Assumptions of Input Sequence\n",
    "1. Equally spaced.\n",
    "2. Periodic.\n",
    "  * Transforming non-periodic sequence may introduce spurious noise.\n",
    "3. Power of 2 in length.\n",
    "  * Padding a sequence may introduce spurious noise.\n",
    "  * **Mixed-radix FFT** is a work-around for this."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the FFT of the sequence $x = [4,0,3,6,2,9,6,5]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def fft(x, n, omega):\n",
    "    \"\"\"\n",
    "    Compute the DFT of the sequence x using the FFT.\n",
    "    \n",
    "    x is the input sequence.\n",
    "    n is the number of points in the sequence, must be 2^n.\n",
    "    omega are the nth roots of unity.\n",
    "    \n",
    "    Returns y as the result of applying the DFT to x.\n",
    "    \"\"\"\n",
    "    if n == 1:  # Base case.\n",
    "        return np.array([x[0]])\n",
    "    halfn = n//2\n",
    "    # Split x into odd and even sequences.\n",
    "    odd, even = np.zeros(halfn), np.zeros(halfn)\n",
    "    for k in range(halfn):\n",
    "        odd[k], even[k] = x[2*k+1], x[2*k]\n",
    "    # Recursively compute DFT of each sequence.\n",
    "    yodd = fft(odd, halfn, omega*omega)\n",
    "    yeven = fft(even, halfn, omega*omega)\n",
    "    # Combine results.\n",
    "    y = np.zeros(n, dtype=np.complex)\n",
    "    for k in range(n):\n",
    "        y[k] = omega[k]*yodd[k%halfn] + yeven[k%halfn]\n",
    "    return y\n",
    "\n",
    "\n",
    "# Initialize the input sequence and compute FFT.\n",
    "x = np.array([4,0,3,6,2,9,6,5], dtype=np.float64)\n",
    "n = x.size\n",
    "omega = roots_unity(n)\n",
    "y = fft(x, n, omega)\n",
    "\n",
    "# Compare the result of FFT with numpy.fft.fft.\n",
    "expected = np.fft.fft(x, n)\n",
    "np.testing.assert_almost_equal(y, expected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.03 Applications of DFT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Digital Filtering\n",
    "Remove High Frequency Noise\n",
    "1. Compute DFT of the sequence.\n",
    "2. Set high-frequency components of the transformed sequence to zero.\n",
    "3. Compute the inverse DFT of the sequence to transform data back to time domain.\n",
    "\n",
    "Similar applications include separating overlapping cyclic phenomena such as:\n",
    "* diurnal vs. seasonal signal in weather data\n",
    "* daily vs. secular trend in economic data\n",
    "\n",
    "#### Discrete Convolution\n",
    "1. Transform the input sequence and kernel to frequency domain using DFT.\n",
    "  * Zero-pad the kernel to be of the same size as the input sequence.\n",
    "2. Compute the pointwise product of the two sequences in frequency domain.\n",
    "3. Transform the product back to the time domain using inverse FFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate use of FFT to remove Gaussian noise from a DC signal with $c = 5$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAEWCAYAAABsY4yMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAB3t0lEQVR4nO2dd3hUxfrHP7M9yaYnhCQkhN5CL4ICgorlotjF3uvv6rV3r+3qteu9tqvYe0NRREEBCYKg9J5QAoGEJKSXTbLZNr8/ziYkkJ5sGvN5njzZc86cMzO7Z9+d884731dIKVEoFApF90PX0Q1QKBQKhW9QBl6hUCi6KcrAKxQKRTdFGXiFQqHopigDr1AoFN0UZeAVCoWim6IMvKLNEEK8JYT4ZwfVHS+EsAkh9E0omyCEkEIIQ3u0TaEhhFgohLiqo9txLCFUHHzXQAiRBkQBbsAGLAJulVLaOrJd7YEQIgHYBxillK4a+z8EMqSUj7TF9Roo/yFwKeCosfs6KeVXzam3q+O9B/2BPlLKMu++64HLpZTTOrBpinpQI/iuxVlSSiswChgNPNixzTmmeF5Kaa3xV8u4C41j4fukB27v6EYomsaxcEN2O6SU2cAvaIYeACHERCHEKiFEkRBisxBiWo1jSUKIp7zHbUKIH4UQ4UKIz4QQJUKItd5RbVX54737ir3/j/funy2EWFezLUKIO4UQ872vPxRCPOV9PU0IkSGEuFsIkSOEyBJCXFPjvHBvO6rqf0oIsbKl78mRbhchRB8hxO9CiFIhxBIhxBtCiE+POO0yIcQBIUSeEOLhFtSZJIR4WgjxB1AO9BVCDBZCLBZCFAghdgohLjqiz/O9fV4jhPhXVZ/rcht5r399je1rhRDJQohCIcQvQojeNY5JIcTNQojd3nvgDSGEqHH8Bu+5pUKIHUKIMUKIe4UQ3x7Rp1eFEP9toNsvAPcIIULqeU/qvHeO7I8Qor8QYrm3XJ4Q4qsa5ep9DxXNQxn4LogQohdwBrDHux0L/AQ8BYQB9wDfCiEia5x2MXAFEAv0A1YDH3jLJwOPea8V5r3Wq0A48DLwkxAiHPgRGCSEGFDjupcCn9fT1J5AsLfO64A3hBCh3mNvAGXeMld5/9qSz4E13j48jtb3I5kMDAJOBh4VQgxpQT1XADcCgUAusNhbdw+09/xNIcRQb9k3ADsQDVzr/WsSQoizgYeA84BIYAXwxRHFzgTGAyOAi4DTvOdeiPYeXAkEAbOAfOBT4PQqY+39cbkY+LiBpqwDktDusSPb2NC9cyT/An4FQoFewGveawTQ8HuoaA5SSvXXBf6ANDTfeykggaVAiPfY/cAnR5T/BbjK+zoJeLjGsZeAhTW2zwI2eV9fAaw54lqrgau9rz8FHvW+HuBtj793+0PgKe/raUAFYKhxnRxgItpjvhMYVOPYU8DKevqe4O1z0RF/jhr1VZUxAPGAq6pdNdr96RFle9U4vga4uJ76P0QzzFX15tV4X5+sUW42sOKIc99G+/Gs6vPgGsf+XdXnmu2vcTwJuN77eiGa37/qmA7tqaG3d1sCk2sc/xp4oMa9cHs9fVsI3OB9fSawo5F78BQgEShG+6G5Hkhq4r1Tsz8fA3NqfgaNvYcd/R3sin9qBN+1OEdKGYhmPAcDEd79vYELvY/mRUKIIrTRaXSNcw/VeF1Rx7bV+zoG2H9EvfvRRuGgjawu8b6+FPheSlleT3vzZe1JzHJvPZFohji9xrGar+sjQkoZUvVH/U8OMUDBEe2q6/rZdbStPl6sUXdEjf01r9sbOO6Iz+EytKeUuvp85PvcEL2B/9a4bgEgOPy5NNSfOCC1nut+BFzufX058EljDZFSbgMWAA8ccaixe6cm96G1f40QYrsQoupppqH3UNFMlIHvgkgpl6ONKl/07kpHG8GH1PgLkFI+24LLZ6J9yWoSDxz0vl4MRAohRqEZ+vqMbEPkoo2we9XYF9eC69RHFhAmhPD30fVrUjMMLR1YfsTnYJVS3sLhPtdsR3yN12Xe/zXbXNOopQM3HXFtPynlqia0MR3NLVcX3wMjhBCJaCP4z5pwPdCeSm6gtvFu7N6pRkqZLaW8QUoZA9yE5obpT8PvoaKZKAPfdfkPMEMIMRLN/XCWEOI0IYReCGHxTnL2avgSdfIzMFAIcakQwiCEmA0MRRuxIaV0At+gTbaFoRn8ZiGldAPfAY8LIfyFEIPR/MNtgpRyP5qv+HEhhEkIMQnNDeVrFqC9d1cIIYzev/FCiCF19HkoNeYdpJS5aIbwcu9neC21jfJbwINCiGEAQohgr2+9KbyLNjE6Vmj0r5qglVLagbl45yyklAeackEp5R7gK+AfNXY3eO/URAhxYY37sxDth9JDA+9hE/uqqIEy8F0Ur0H4GM0fng5UTcLloo2C7qUFn6+UMh9tJHc32kTcfcCZUsq8GsU+R/PFfiObEEdeD7eiTcBmo7kFvgAqW3iturgMmITWh6fQjFFbXv8opJSlwKloE4OZaH17DjB7i9yK5jbJRnsC++CIS9yA9rnlA8OA6tG5lHKe91pfCiFKgG1oE+1Nadc3wNNon1sp2qg9rEaRj4DhNME9cwRPAgE16mnKvVPFeOAvIYQNmI82R7C3Ce+hohmohU6KToEQ4jmgp5TSJysdvWF4KVLKx3xx/ZYghLgabdJxcge3Ix5IQXv/SzqyLYq2RY3gFR2CN9Z5hNdlMAEtjHJeG15/vBCinxBCJ4Q4He0J5/u2un53QWiLs+4CvlTGvfuhtDgUHUUgmlsmBi2i5yXghza8fk80n3c4kAHcIqXc2IbX7/J4Y84PoUW6nN7BzVH4AOWiUSgUim6KctEoFApFN6VTuWgiIiJkQkJCi84tKysjICCg8YLdCNXn7s+x1l9QfW4u69evz5NSRtZ1zGcGXggxCC00rYq+aCF9/6nvnISEBNatW1ff4QZJSkpi2rRpLTq3q6L63P051voLqs/NRQhR74ponxl4KeVOvGqHQkvCcJA2jJJQKBQKRcO0lw/+ZCDVu8JQoVAoFO1Au0TRCCHeBzZIKV+v49iNaHKrREVFjf3yyy9bVIfNZsNqbUgrqvuh+tz9Odb6C6rPzWX69OnrpZTj6jrmcwMvhDChLTkeJqU81FDZcePGSeWDbzqqz92fY62/oPrcXIQQ9Rr49nDRnIE2em/QuCsUCoWibWkPA38JR2eeUSgUCoWP8amB9y6FnoG2ZFyhUCgU7YhPDbyUskxKGS6lLPZlPQqFQtFueNyw4WNwOzu6JY2ipAoUCoWiORxYDfNvg33LO7oljaIMvEKhUDQHmzdepKKoQ5vRFJSBVygUiuZQ5k1QVVnase1oAsrAKxQKRXMoy9X+V3b+/CjKwCsUCkVzqDbwagSvUCgU3YsqF41djeAVCoWie6FcNAqFQtFNUS4ahUKh6KZUu2g6//pNZeAVCoWiqTjth10zagSvUCgU3YjyvMOvlQ9eoVAouhFV7hlrlBrBKxQKRbfCa+APGWORKkxSoVC0N6//tpvz3vyjo5vRPfFG0CTlBiHcleCq7OAGNYwy8ApFN2PR9mw2phdR6XJ3dFO6H14Dv19Gadud3E2jDLxC0Y2wVbrYkVmClHCwsKKjm9P9KMvFpTORLUO17U4eKqkMvELRjdh0oAiP1F6nKwPf9pTlUWYIpRR/bVuN4BUKRXuxbn9B9ev0gvIObEk3pSyXEl1IDQPfuSdalYFXKLoR6/cXMigqEKNekF6oDHybU5ZLPkGUSj9tW43gFQpFe+D2SDYeKGJCnzBiQ/zIKFAumjanLI8cd+DhEXwnD5VUBl6h6CakZJdgq3QxLiGUuDB/NYJva6SE8jyyXIE1RvDKwCsUinZgXVohAGN7h9Ir1F/54Nsahw1cdjIcAdiUD757s+tQKe+u2NvRzVAoqlm3v5CeQRZiK/dyons1heVObJWujm5W98EbA58vg3BiwIFJuWi6K68s3sVTPyVzqMTe0U1RKABYn1bA2IRQxNInOWXnY4BUo/i2xCtTkE8wBp2gTPirSdbuiK3SxW8pOcDhx2KFoiPJLKogs9jO+Dgr7P8Dg6ucYMqUgW9LqkfwgcSH+2PD79h20QghQoQQc4UQKUKIZCHEJF/W114s2XGISpcHqB13rFB0FOv2awONKX4HNF8x0EvkqcVObUm1gQ+mT3iANtHayUfwBh9f/7/AIinlBUIIE1TNTHRtFmzJJDrYQnyYP+v3qxG8ouNZn1aAv0lPQum66n19jQVqBN+WeA18AYEkRARQtNfv2PXBCyGCganAewBSSoeUsshX9bUXxRVOlu/KZebwaMYnhLE9s4Ryh5rIUnQsa9MKGR0fgn7fcgiOB2CofzEZKlSy7SjLw64LwOIXQFiAiVLph6eTa9H4cgTfB8gFPhBCjATWA7dLKctqFhJC3AjcCBAVFUVSUlKLKrPZbC0+tzmsyHDidEtiXFnYKiRuj+SjH5czJFzv87qPpL363Jk41vrclP5WuCTJWeVc0MeN5+AaMnqdRWxpNjHuTD5Lz+1y71dn/YyH7NuOh0ACdC6yDuwjCn/sxftZ0wZt9VWffWngDcAY4DYp5V9CiP8CDwD/rFlISjkHmAMwbtw4OW3atBZVlpSUREvPbQ7vv7+GuDAb15w9nRK7i1c2/IorJJ5p0wb4vO4jaa8+dyaOtT43pb8rduciWcOVgxzoslzEn3gFLNrOAFcZBQU6TjzxRIQQ7dPgNqDTfsb7XyalMIzekWGMTuxFyW4/LFS2SVt91WdfTrJmABlSyr+823PRDH6XpaDMwR978pg5PAYhBMF+Rgb2CGSt8sMr6qGgzIG7St7RR6xLK0QnYFD5etCbIH4SBPci0pNDucNNQZnDp/UfM5TlkesJJDLQTIDZQCl+CEeptsK1k+IzAy+lzAbShRCDvLtOBnb4qr72YNG2bNweyVkjo6v3jU0IZeP+Qp9/iRVdj53ZpUx8Zimfrzng03rW7y9kUM8gTAdWQNxxYPKH4DiCK7MBJRvcZpTlkuWyEhloJtBioFT6I5DVUUudEV/Hwd8GfCaE2AKMAv7t4/p8yoItmfSNCGBodFD1vnG9QymtdLHrUOcOl1K0Lx6P5MHvtuBweVif5rtQWpfbw8YDhZzYS0D2VuhzonYgJA5TZQEWKlUkTVvg8SDL88h2ayN4q9mgxcFDpw6V9KmBl1JuklKOk1KOkFKeI6Xssr6MnFI7f+7N58wR0bX8meN6hwGH45AVCoDP1hxgw4Eigv2M7MjyXShdSnYpZQ43J1tStB19p2n/vZE0sSJPiY61BRWFCOkhXwYTafW6aGTnV5RUK1mbyMKt2XgknDkyptb+uDA/IgPNPh2lKboWh0rsPL8whRP6h3PlpN6k5pZhd/omP+o67303tGITmIMgZrR2ICQOgEGWItKVbHDrqYqBl4H0CNJcNIdH8MrAd3kWbMlkYJSVgVGB2o78VNjwMUIIxieEqhG8oprH52/H4fbw9DnDGRYTjNsj2Zntm8f4dfsLiQ62EHBwJfQ+AfTewLhgzcAPC2h9LPyenFIlWuY18HkEV0+ylsjOryipDHwTyCquYG1aIWeNqDF6T3oW5t8GJZmM7R1GRmGFEh7rAOxON/M2ZuDpJJPci3ccYuG2bP5x8gASIgIYFqPN12zP9I0RWL+/kBkxlVC477B7BiAwGoSevsaiVvngXW4P57yxinu/2dz6xnZlaihJRlrN+Bv12IRy0XQLftqSBdRwz7idsOsX7fXe5YzrrWVYV8Jj7c/8TZnc+dVmlu/K7eimYKt08egP2xgUFcgNU/oC0CvUj0CLge2Zbb/i8WBRBVnFdk7zq/K/n3j4oN4AQbHE6fI4WFTR4iivvXll2CpdLNyWzab0otY3uqviVZIsEsGE+pvQ6QQeo1U7dqxOsnYXFmzJYlhMEH0iArQdaSuh0vuF3becoTFB+Bn1SnisA9iUUQRoIawdzYu/7CS7xM6/zxuOyaB9tYQQDI0O8slEa5X/fVjlJrBGQeTg2gW8sfBOt2zx02Wyt91mg47nFqYgO3HMt08py8WDQB8Qhk6nBVl4TN5oOuWi6bqkF5SzKb2IM2u6Z1J+AqM/DDwD9i7HqBOMjAtWwmMdwBavgV+cfAiX29Nh7diUXsRHq9O4/LjejPU+0VUxLCaYlKzSNl8rsS6tkACTjuDsVVp45JGrVUNqxMK30E2zI6sEo15w72mDWL03n99357W22V2T8jxsumDCgw7rJeotVjwINYLvyiyocs+M8C5uklIz8P1OgoGnQWkm5O1mXG9NeKzsWJ+MakfsTjc7s0vpGxlAQZmDtR3kInO6PTz43VYirWbuPX3QUceHxgRR4XSzL6+sjrNbzrr9hZwVXYwoy63tnqkiOA5zxSH0uFu82Ck5q5QBPQK5clICcWF+PLcwpdPMd7QrZbkUCs3/XoW/xYRddG5FSWXgG2HBlkxGxYUQF+b95c7cqBn1wWce/lLtW87YhFDcHsnmY9lP2c6kZJfidEv+Pq0/ZoOOX7Z3jJvm/ZX7SM4q4cmzhxFkMR51/PBEa9v54UvtTnZml3C6v9f/3qcOAx8Sh5BuokXLZYN3ZJYwJDoIk0HH3TMGsSOrhAVbs1rR8i5KWR65niAiAw8b+ECzgXLhr1w0XZW9uTa2Z5YcHr2DNnoXem30HtpHW1CyN4kx8aEIoRY8tSdbMooQeJgScICpAyNZtC273UeX6QXlvLJkFzOGRnHasJ61D+5cCK8kMqBgOSa9rk398BsPFOGRMLxyE4T1q457r4U3VDIxoKRFi51ySyvJs1UyJFoLDZ41MoYh0UG89OtOHK6Oc4d1BLIsl2xXYC0Dr61mVQa+y1LlnplZy8AvgN7Hg3+Y5vPsOxXSVhBs1jGwR6Ay8O3I5vRiZvtvoMeXZ3BV5B6yS+xsOdh++txSSh7+fht6IXhi1rDDK5zdTvjlYfjiYihOx7B7IQN7WtnRhqGS6/YXYhIuwvLW1O2eAQjRVrMOCygmowWLnaomWMeb0iBvDzqd4L7TB7E/v5yv1vpWX6ezIW255MnAWi6a6lh45aLpmizYksn4hFCig70r1vL2QG6K5p6pou90sBdD1mYlPNbObMko4pSAvQAcl/05Bp1o12ia+Zsz+X1XLveeNoiYEO89UpQOH5wBq1+H8ddDn6mQuUmLpMksabMolHVpBZwVcQjhKKsd/16T4F4A9DMVtWgEn5xVgh43I5ZcBq+PhU/PZ5rYxHEJIfx36Z5jZ77J5UBXWazFwAdaqncHWgyUSIuaZO2K7DpUyq5DNs6qKU2w8yft/+C/Hd7XZ6r2f28S4xOU8Fh7Yat0sSfXRqLcDQiM+5czO66IRduy2iWUr6jcwZM/7mBkXAhXTErQdu76Bd6eAjkpcMH7MPMliJsIuSmMiDKRX+bgUEllq+t2uT1sSi9iZkAKICBhSt0FjX4QEEkvXR7ZJXYqXc2TS0jOKmFsYBHCWQ4DToXsbYjPL+Sj8r8zs2I+nyRtbXVfugTlWuRQvncVaxUBZj1Fbj+kctF0PRZszkQn4IzEI/zv0SOrH30BsPaAHkNh33IlPNaObD9YjFE66VG+C8ZcAcYArjUsJC2/nF2HfC/f+tyiFIoqnDxz7nD00gWLH4XPL9JGzTcth8TztYIxo0C6GWM+CMCOrNa7kJKzSil3uBnh3AzRIzR3YX0Ex9HDk4OUkFnUvFj45KxSpoR4wyKnPQh3bIXz38MSFMETxo+4fNUZ2Offo8l2dGdqrmKt5YM3aiN4e+cd0CkDXwdSShZsyWJi3/DDH2jpIUhfU9s9U0WfE+HAn/QKFPRQwmPtwpaMYoaK/eg8Tuh/Coy+nL5ZC4kShT5300gp+WlLFueOjmVoQAl8OBP++C+MuxauWwLh/Q4Xjh4FQD/nbgC2H2z9aG/d/gL8sBNeuLl+90wVwb1aFAtf6XKTmmtjpCVH2xExEAwmGH4BXL+EA+f9yGLPWIwbP4DXxsBnF8KepS3sUSenSofmSANvMVCqJlm7HtszS9ibV1Z7cdOuhYCEwTOPPqHvNHDZERlrGaeEx9qFzRlFnBiwX9uIHQcTb0FIN/eH/c7Cbb4N48u1VVJid3GaaQu8NQUObYfz34MzXwGjpXbhoBgI6IE5ZwsJ4f5tEkmzbn8hpwfuRXicdYdH1iQkHnN5FiCb5YfffciGyyPpS4YWjWO21joeP2Iqq0b8mymO1ymZeA9kboJPz4Pdi5vfoc5OWT4A5cZQrObDWU6tZj026YdwVWgT650QZeDrYMGWLPQ6wemJNcLeUn6C0ATNHXMkvY/XQif3JlULj2UXK+ExX7Ilo5gTLGmaqFZwLIT1gcFn8rfKhezPziWtjRcV1SQ1p4yb9D8yY+OtmgG/MUkb2daFEJqbJmsTw2KCWy06JqVkXVoBM627D6fna4jgOHSuCqL0pc2SDa6KoImo2AeRRy/eArhzxkDyRAhPls6Cf2zUdmZubHIdXQbvCF5YI2vttpqN2ggeOu1E6zFj4KWUfL02nU/+3F8d/bA5vYi0vDKKyg/nzdTcM5lM7h9BWIBJO7myFPYmae6ZupIXW4Igdmxt4TGlS+MzCsscHCgoZ5Brp/a+VzHpViyuEi7UL/fpoqfUnBJuM8zD3ns6XL8EIhpJuB4z2jvRauRAQTkl9paP9jTV0kpGuzYdTs/XEN74+FGBpc0awSdnleJnBHNR6tEaN15iQvy4alJvvtuQwa4iCdaeUJjW5Dq6DGW5ODHib60tQWGtlfSj/cJzm4Oh8SLdg43pRdz37ZZ6jwuhrUwLtBg5WFTB7SfX+NLuXgxuR93umSr6nggrXmJomNSEx9IKa7t4OilSyloZqroCWw8WE0oJwfYM6HX94QPxx0Gv8dyS+Qu3bruAm07sV/9FWkHBwV1YhR05/BwtUqUxokeB9DDOok20JmeWcFzf8BbVvX5/IaGUEF66E8Zd2PgJ3sVOQ/2L+a0ZPvjkrBImR9oRBRWa/70e/m9af75ck84Lv+zkndAEKNzf5Dq6DGV5FIpgIoNqu9+6Qtq+Y8bAr07V/GiL79TCGosrnBSVOymu0P6KKpyUVDg5YEvBFPs9vxf/SsrqCML8wghLXkhYaBRhRh1hRamEWcIINgejEzUegPqcCL+/gDF9NSPjQju18JiUkuW7cpnz+16Ss0r4+fYph2P9uwBbMooYqfNGbsSOq31w0q1Ef3MVEQd/I7t4Aj2DLUdfoJXIrG0AiJ7Dm3ZCzCgABrp3A33ZkdVyA79ufwHTzTu1jcYmWKF6BN/fXMhHOU1z0UgpSc4u4Y74XCig3hE8QGiAiZtO7MuLv+4if1g04blrm1RHl6LMu8ipxgQraJOsJXTupB/HjIFflZrH4J6BDKjKyFQH2/O3c8OvL2E0GUkvLWVL3iYK7YVIJISYYfEN1WV1QseAkAG8PO1l4oPiIW4CGPy84ZLX8r/lqZRVuggwd5632OHy8OPmTN5ZsZeU7FKigsyU2l28sWwPT53TRGPVCdicUcz0gAPg0h1OUVfFkLNwBsVzfdFP/Lrjaq6silFvQwKLUvCgQ9djSBNPiAZrFIEF24mwDmmVH35dWiH3WXeCK+jovteFJQRMgfTS5VNQ5mjSPZldYqeo3Mlwk9fNFVn/CB7g2sl9+HDVfn7L9ueCsoMIVyUYzA2e05XwlOWS6669ihW0OHib7Nwj+GPCB1/pcrMurZDj+0XUWyalIIUbf72RQGMgX8z8gu/P+Z7ls5ez8fiXSNqfwXcj7ubdU9/l+anP88CEB7gu8ToOlR/iioVXkJyfrN3QvSdpE62dTHisxO7k7eWpTH1+GXd/sxkp4cULR7LivpO4aHwcX61Nb3Vat/ZkS0YRE0x7IXLIUdEd6PQYj/8743W72L0hqc3rLne46OVIpdgvrnH/dxVCaG6azE0MiwlqsYEvsTvZeaiUMe4ttdPzNVZ3SBw9PNpEYVP88FUTrPGeDM2v7hfaYHl/k4HbT+7P6sJABFJbzduN8JTmkEfQUSP4QLOR0ioXTWvkCla/wbBtz7aihfVzTBj4jQeKqHR5OL5f3Y/Fuwp3ccOvN+Bv9Ofd094lxnrYd67fuZBwvR8DEmdzXPRxnNHnDC4bchn/GPMPPjrjI0x6E9f8cg1rs9dqbprcFMaGVXYK4bGs4gr+/XMyJzzzG88sTKFPRAAfXDOeRXdM4YKxvTAZdNw6vT8CwRvL9nRoW5tKod3DoRI7fSt3Qq+xdRcafRl2vZWJ2Z9TUOZo0/r35pYxRBygIryOaKqGiBkNeTsZGWVkT05pi8S6NuwvJJZcQuwZTXPPVBHci+BKLXS0KZE0VZo5YeV7Gx29VzFrZCzpsoe20Z0mWqVEV55Hvgw+ysBbjDpNTRJa56LJ3ITVtrcVjayfY8LAr0rNRydgQt+jV/ztKdzDDb/egEln4r1T3yMusIYqn8cDO3+G/ifXOZnWN7gvn5zxCVH+Udy8+GaWWrWMT0FZqxkU1XHCY2WVLu79ZjNTn1/Guyv2cuKgSH68dTJf3DiR6YN61JpUjQnx45IJcXyzLoMD+Z1/FL+v2EMfkY3ZVXK0/70KcyClwy7ndN0aVq9f36b178/MIk6XiylmRPNOjBkF0sNx/gdxumWL5CyWJudwomG7tlGfwFhdBMdhKc8Eai92cnqcOOuI307OKiUu1II+b1eD/veaBPkZyDN6V30XpTW9bZ0dRxk6t/2oVaygZevymLwu39YY+PI8nMbgVjSyfo4JA786NY/hscFHaXXvLd7L9b9ej07oeO+09zRfek0yN0JpVt2rV730DOjJR6d/xODwwdy15XW+DesBe5cztnfHCY+99OsuvlmfwaUT4ll+73Rev3QMw3vVfwP93/T+6HSCV3/b3Y6tbBn7SjyM0XsnWHvVY+CBiJNvQwodxnVz2rR+2/5NAAT3aYL/uybeFa2D3NqTUnMXPDlcHhZsyeTckD11p+driJA4dPYiwo0ODhTYWJu9lkf/eJSpX05l8peTuSvpLn5M/ZFibxrK5KwSJkU6wFFabwz8kQghMIVE4xTG7jWCr9ahOdrAA5jN/rgwtM5FU5aHw+QbA+/TGUAhRBpQCrgBl5Sy/m+kjyh3uNiUXsR1k/vW2p9WnMb1v2ghdu+d9h4JwQlHn5yywKv9fmqDdYRYQnhnxjvctfwuHpd/UJiZxNjj7uezvw6w61ApQ6KD2qo7jbLtYDEfrtrHZcfF88TZiU06JyrIwuXH9ebDVfv4+/T+h3PPdkL2FXu40X8/YG3QyIngXuwIm8Hx+T9jK87HGtyyqJWjrntIG0EbmzuCD4oGa0/CS3bgbxrcbOng33flUljuYIRxMww8ue71GPURHMdeo4HIHj/zY+Fuvv0lH3+DPzN6z8CkN5GUnsTi/YvRCz0jI0eR4enJJQFeN2UzfkiiQwPILosirj0NfFE6OCua7EpqNt5k23kyiPCAow18gMVAhTOAwNZMspbl4Qxo4oR9M2mPEfx0KeWojjDuoEUdON2SSTX87+kl6Vz363W4pZt3T32XvsF96z45ZQEkTG50kgnA3+jPaye9xszgQfw3QMemgtcBT7u6adweyUPzthIWYOb+yaFgywF30yRdb5nWD5NBx6tLO+8oXkrJvmI3o/Wpmk9bp2+wvO74W7EKOwd+fbPN2hBUnIxNF6itYG0uMaMQmZsYEh3Etsxc0kvT2ZSziaX7lzJv9zy25W3D6al7EdS8jQeZ4J+NqTK/ye6Z/Ip8Pkv+jIt3f8TZvWLI9P8L4YziuSnPkTQ7iacmP8Wjkx5lyYVL+GLmF1w3/DpyywoxR/3Eq/Z3OCe2J//NXsHm3M1NUuiMDfEjzR3ZfiN4Ww68ewq8NwPKfbSw0LuK1WkOr06kXhOr2UC5CGi5i0ZKn7poOk8Mn49YlZqPQSeYEFIKL00hY9AMrq1IxuF28O6p79I/tH/dJ+bugrxdMOHGJtdl1Bn59wn/JvSLGXzKQkLic1izrydXTOzdRr1pmE9Wp7Elo5jvJ+0h6A2vmiEC/ELAPwICIsA//PB//wiIGgp9pxEZaOaqSQnMWbGXv0/vR/8e9YeT1suWb0C6YeTFbdmtatILKnA5HcTY90Ds6Y2WHzJmMmt/SqT/zo/A/QDoj06n1xzcHklsZSp5wYOwNmEEXemu5Ntd37KveB/59nzydbnkWUrINP4DF5X87bujz7HoLQwNH8rIHiMZFTmKMncZJXYni5MP8Vr8PsikTv2Zcmc5GbYM0kvTySjNYE32Gv44+Adu6WZIcH/uzS/EGHwlTx2Yyhl9Tqs1D6MTOhIjEkmMSCTccTYPL1jOI0O/ZXXlTj7Y9RXv7vycc/ufyxPHP9HgoriYED/2uiKYXPgnQsrmPWU0F7cL5l4L9iJtEeLy5+EMH0SieA08R8gUVGG1GLEJP6JaOoKvLAW3w2cuGl+P4CXwqxBivRCiUUu5N7eMb9ZpIVZOt4fZb69m3sYMACocbma/vZofN2uTRSV2J7PfXs0ir7BUqUMy++3VLNlxCICcUjuz317NL9uzGB0fQuHK9zg3/3ou3ruf8rJDPNn/bh79poA/92oLoFJzbcx+ezXrvRIDO9cuZnblI2wO1L5M2zOLmf326uq8mpvTi5j99mp2Zmsf7Pr9Bcx+ezX7ZDT3uaycYx9Had5Yluf/j3JnOSt35zH77dVkFmlRDEk7c5j99mpySjXNmiU7DjH77dXVUR+LtmUx++3V1cvaf9ycyey3V1Ph0DS9523M4Jm/KnC6tWiMd1fs5ckFO7glPp2Rm5/ki5CbuCzwfTjxfkg8n090s7jy4Cz2Fuxm2d6F3PxnKlN+KeGhn69l14b3mPN7KrtzSvEz6vnPkt28mbSHWz/fUP3ZvLp0N3d8eVhn5OVfd3LPN5urt5+bu4IHv/oL5t0EvzzM0wt28M/vt1Uff+LH7Tzx4/bq7X9+v42nf9pRvf3gd1t4blFK9fY932zm5V93Vm/f8eVGnluUwjCRhl66uDV1PG8mHY78ufmT9cz5/bBs7fUfreWjVWls730loa5crnjtZz5ZnVZ9/LJ3/+SLNYezEs1+e3Wj997Hf6QygHQKQkbUuvcKyhxH3Xtnvr6E0z+5i2fWPMP8nSv45Y8hFNsHM6zSwSj9ZGy7H+Dc6Ef53yn/4+XjP6Vf2avcMugFLhx0IcWlAbz5k4lbf36BhzIeYsZnV+HUZbPRsYO/oofwavIGprz0LTf//CSX/3w5E9+7lOHPvMe539zEHcvu4Nnl37J09XDO7n0d3836jn8MfZuf8u4k0u2izOFmwZaseu+95KwSzDKCX3edwctiIMtnL2dS4K18tjSST7d/XX3vzX57dfW99826dGa/vZpeoX4ckD34smwcl835o/q9/WR1Gle9v6Z6+/2V+7j+o8MLoub8nsrNnxyeDF+w19H4vff6p5C2As78D8+FPcmDK11aQh7g6Z/a8N4ry+UOx/9x0HXYAN/6+Ybqe89q1vNwxWXMyTy8avr6j9by/sp91dtXvb+m/nuvLJfZlY/wc4nmRWiu3WssSszXI/jJUsqDQogewGIhRIqU8veaBbyG/0YA/6g+pKSkkGRLxeWRFBXZSU4uI7R4D5VubXvHjjICC3dR7tS2t23bjiVvJ2VlZRQV2dm6bSuGnGSKKj0UFFayr8jD8dY8irZ9S4r5VswGMx/mVhDy4/3YdQ+zeYMN+wEjWTYPRUWVbNiwkdJ9eqzbfsZtOIP1u9IpPJTJ/hI3RUUO1q1bR26Qnr3F2vbatWvJCtSxu1DbXrNmLf7+AzktO40fOA2XKZXzvr6QaYbrKSqysHr1asL9dGzJdVFU5GTVqlWEmHVszdG2//jjDwJNgm3Z2vbKFSvxNwp2ZGnbv6/4HbNekJzpwu12s3z5cgw6wSeryjHLSu7Mf4JSvxhWWoeyN9fBrYUHyXHmkOr2UISVswMqADMOjx+Um/gtIICftrxCjC0PsyORk3oZ+WlLFpXFuRRXSpKSkgDYt8/BoTJP9XbafgcFdklSUiFIN47ti3HrDByMOZ3Y1a9TbrCQHTKWpCTNh5mRoSW6SErSRkQHMysx6SApSZOjzcyspMQoSErSFtdkZ1fiKBIkJWk38qFDdnIrJNP02hcrr9zD3r17SUL7IuTm2Ul1FZDk0Yx0Xp6dPbKIaP8+7PHE4Mrby65dfiRVpgFQWFjBzp2lJJVr4WlFRRWkpNgavvc2p2ERTrJdQRQVFVXfe6UOWX3vebI3823uYpIL4gnVl3JLwi1EMZg52ZWcYY7gprTP+FH4scwVQskBIy6TizxbPhWllRizTEwInUCM3zjKzXamhF1CqdjG4iwbOn05HxvL+MxShnvzW1SW/Q1X5mairZJ4/QA8hnBODT2PgUHBFJSEs7BUx7AyMwc3H2R73gE8woC5OA04kdUbt1FU5K7z3vszxU6oyYO0l5HtspK5eiN97X1YrUvnxbUv4jpYTlZ+HEVFzup7LyXDSVGRi0P7ykmXPegvMikpyKm+V3YdcFJQ4K7e3pPmJC//8HbqPie5RYe3HZUOckpy6r33DiXvg5xkMuNPY1dRNBmGAPxFOnlf/B/bhj9EenolDg9tcu/127MRJwm4Xc7q+nNy7FjseSSRQWlBJU6po7K8pPp41b2X5NJkGwoK7OzaXVznvRdUrP2wVEgjSUlJzbZ7pY5GXGdSynb5Ax4H7mmozNixY2VLWbZs2VH7Fm/Plr3vXyDT5z4kH30tQY78aITceGijlOUFUn5zrZSPBUn5zslS5u2pfWJJlnZs+fMta8ymL6V8LEju3PC77Puv5+XYjyfIyV9MliszVlYXKat0ys3phfLrtQfk0z/tkNd8sEZ+sHKv9Hg8ze7z4u3ZcvT9n8uiZ4bI5Jf7y8vnXygTP0yUiR8mygmfTpAXzr9Q3pN0j3xtw2ty/p75ckvOFllkL5JSSlm4N0k+9WpvOeLDRHnCFyfI9zd/Ioc9+pO86eN1Te/vn29r79fmr7XtVa9L+ViwlO/OkNKW1/TrNMJFb62SS584XcqXhjT5nEqnWz7+2D1a+/b+3qr6l3z9hpSPBcni1Lrfm/XZ6+XM72bKxA8T5ROrnpCllaVHF3pxkHR9c73s++BP8oVFKY3WOffnpbL3/Qvk4o+flRn/CpWrt30h04rTpMPlaF7jP5gpy96cLnvfv0D+uPlgnUXcbo8c9ugi+ew3v2vv1+o3q48V2YvkGd+eIU/66iSZW55b5/kZheXytAfe1M7d+m3z2leDur7L1eSnSvnvOCnfmiKlo+Lw/t9f8n7Gy1tcb1145l4n9z/aXz790446jz8+f5v84dEzpXwlsWUVJP8k5WNBcu38d1rcRmCdrMem+mwEL4QIAHRSylLv61OBJ31VX12sSs0nzFBJ5v7P+S7cyjXDrmJUj1HawQveg0FnwE93wVuT4bR/w9irNb/hzp+1Mg2ERzaIN41fP9t6TJXDGa1PZJ/4HzcvuYVYZlGeM530AjtV81Ymg44egWZ+S8nhQEEFj8wcgk7XNP9lucPF0z9s5G3rK7xpreTLwABCyg/xxPFPMCV2ChF+EQ36TUP6nMjDAy7hwvXv8FziIF7e+BxhA+NZnHoa2w72JzG2Ed9g8UFY+iT0O4nKIeciXB5Mk/4OQbHw3Y3aBNjlcyGsnonsJuL2SLYdLGa4YQ/ETmzyeSaDjvJBF1Cw80tCVr2Grk896e2agC5nGy70BMXXjk4qd5bz3w3/5YuUL4ixxvDOqe8wMbqeNsaMRp+9iQE9rqh29zXE6ixtkvwE15/4BcYRO3R2y3zbIfFY8rSEHPUtdsoorMBW6WJcgDfJR40QyWBzMK9Me4XLf76cu5Pu5t3T3sWoqz2nERVoJlP4cLGTswK+vlLr/0Uf19ben/h/sO4D+OUhuHF5oxPwTcVdmqMl+rDWLb0QaDZQ5DEj7SW0aMahahK3C06yRgHzvMbFAHwupVzkw/qOYlVqHreHr+BfVjOxlkhuGXlL7QLDL9D0tH/4P1hwB+xcCLNe07Tfw/o2L9a4JkHREDEIfdrvjI6fyOIt+SCuwS/6Bw4G/0BoVCo3j7iHkbExDIgKpHegQP/78/yy38Otf7goKKvkhQtHYtQ3PkXyyq87OVH/b+7tYaPQEMBFg2Zz6+hbCTY344Y56REGpizg3QP7WDrzOZ7f8F/8e7/D/y3ZyBfnP0WsNbb+cxfeBx4Xvw98iHueTyLQYuDDayYQN+wcCOwJX1yiRTpc+nWDcetH4vK4+DPrT0oqS3BLNweLbBj9N/KbsRyX1Q/39o9xSzdu6cakMzGz70zC/eoOhTxpRG8+234St+7+QYu2aCjFXQMEF+/koCGO3jV0Vv7K+ovHVj3GQdtBLh18KbePuR1/YwMSBtGjYOdCRg808Ftaw5EXUkpWZbqYEmfCL2OlNuHf0onL4Dh0tkNE+tUvV1AVmz9Qp6leHnn/DwobxOPHP84DKx7g5XUvc/+E+2sdN+h1BAaFYnMGY/WFgf/5Xsjeqt1LoQmA5rN2eyQWowVOeQy+vQ42fwGjL2+TKt2lOXWuYq3CajFQIv21ydKWTCx74+ydRt+EUvvMwEsp9wIjfXX9xigoc7A3u4D82IWkmYy8Pfkp/Ax1KCYGx8Ll82DtO1pezTcnah/WxJtbFwXQdxps/ISnrx/E1uwKBkZZSQg/iwX7vueZv55hSekDnBH1Ev30FfDhlZC9ldOBP8P7c+Pmy7m+3Mn/Lh+Dv6n+j2h9YTrLDj3IoSgbI8xRvDnjDYY2dwk9gCkAznwF8en5nJK2gcnn/MAtP77C2qJvOGveLK5JvJrrEq872nAl/wgpC5gXfiN3zsthcM9AMosqOPfNVXxw9XiGx0+E6xbDZ+fDh2fC+e/CkIafitweN4vSFvHW5rdIK0mrfTAaniIM8lZrfzV4fdPrzB40m6uHXX2UoZ86IJIPxRhu43tIWwlDZzX/PQJ6OVI5GDKB3oDNYePl9S/zza5viA+M58PTP2RsVD3SCTWJGQVITrAe5IsSf/JslUTUMzrcnllCpk3y3LBUyHW0/IkSvKqSktEh5fWm7kvOKkEnINpxAEyBmkjaEczsO5Ntedv4NPlThkUM48y+tdsUG+JHVn5PBhS1sWzwhk9g4ycw5R4YeFr17hs+Xkdqro3v/+8EwhPPh7/egqX/gqHnHK1T1AJEeR55MpH4egx8gNlAvvRHSDc4y7XvUnMoyweTFY/eN+Js3XYl65978znJbxGfBBo4M3Icx8ceX39hnQ6Ouwlu+l37Inic2g3SGvqeCM5y+lTsYNbIGAb3DMJiNHDBwAv4+G8foxd6rlp4JZ99eiqyKF0blVz8ORF6O9+ZH+e0fc9ww9tLKKxjltzmsPHcmuf5qPgFKk3FPGLqxycXLWmZca+i/ykw4mJY+QqW/FRePeMe9JkPEOgezZwtczhr3lnM2z0Pt0eL4pH2Ysp/uJud9Obh7Knce9ogfrxtMt/93/GYDTpmz1nNspQciOiv5SmNGgZfXQ5/vlVn9R7p4de0Xzl//vk8sOIBjHojL534EvPPmc/P5/7MaUGvctX+Efy2P5MV5/3CqktW8delf7Hu8nX8cPYPnBx/Mh/v+JgzvjuDl9a9RH5FfvW1/Ux6wgZOpBwLct/vddbfGIW5WfSgkIrwwczdNZdZ38/i293fctXQq5g7a27TjDtUr2hNFFqURUMLnr7feBC9gOMqV0NApKZY2lK8uvDD/IvJKKzbRZOcVUJCRACGgp2ae6aeAc5d4+5ibNRYnlj1BDsLdtY6FhNiYb+njWPhs7bAz/do4aHTH6re/ceePJJ25pJeUMEtn23A6ZGaq9WWDatebX29Hg8Ge2G9q1ihDTThy3K1kGUf0W0N/Krd2eT3XEkAOu6d/mLTToocpBmjW1Y1y51QJ71PAKHTMkEdwbCQQXwVMJITymw8G+zH/ePOpLzPFC2hyN//gkm3crFhOa/m3cBbr/2bvbkFpJemszFnI1/v/JpZ38/is+RPObfUxgflPZl94dfo2sLnePozmrzs/NsINOm4+YSx7E85l4dHv0FPa08eXfUoF/x4Ad+lLGbx67diqcjh/bA7mX/7dP4+vT9GvY7+PQKZ93/H0zcygOs/XqeFg1kj4aofYdDfYNH9sOghTecHzQ2x7MAyLvrxIu5efjcePLxw4gvMPWsupyacSp/gPsQFxbE708QMQyZ+fnGEBMYQaArE3+iPWW+mb0hfnpnyDN+f/X21oT/929N5ce2L5FVoj8BTB0fzl3sQzj1Hfx5NIXvXWpb6+/G4bgVPrH6CGGsMn5zxCfeMv6fuJ8P6CIyCwBhiyrXoifqUJV1uDz9szmRshAfT3qXae9eaz9irC9/PXMjBwgo8dUho7Mgq0VZd5+5s0D1p1Bl58cQXCTIFcceyO6olDkCLhd/liNAGLU1cZNcgFUXw9RXgFwYXvF/9Hkgpef6XncQEW3j+/BGs2VeghULGTYDE8+GPV7X5IWBr7lbOn38+Fy+4mEX7FlUPUhrFXoROujQdmnqesmpndWrBYqfyPO3H20d024VO+zNeIyVY8K/eswirxzdbJwaTNtpsLX4h2mrLfcuBhw/vL82Gb64h+MAqXh17Le/HDeK1Lf8j5aeLObvf2eTb88mzOMkdfQo5eTspEN/w+c/f1rr0oMB+/OtgKglOCzH3ztXa3Bb4h8EZz2l+zL/e5spJN/Luir38tNbMp9d+yuK0xfxr9Ys89tddjLfacYTN5Jmrrz5qQrhHkIWvbpzE3z/fwIPfbeVgYQV3nzoQMfsTWPQA/PkG0ujHH4On88bGN9iWv434wHj+Pfnf/K3P39AfYcgcLg8pWUUMMu+iIPwE6nvw7hPch2emPMNNI25izpY5fJL8CV/t/IrZg2YzPfoiFnmGMb3wc+0zCOxZz1WOZv2h9fwr9WVSoyKJM1n475T/Mj1uesszYcWMxpSzhV6hl9erSbMqNZ/c0kruj9oBGaUw5KyW1VVFUC9AEKfLx+H2cKjUXivJS4ndSUZhBVePDoFdhxrVoInwi+Dl6S9z9aKreWDFA7xx8hvohI7YUD82eyI1l0VJRrWvvEV4PPD9LVCcAdcs1Bboefll+yE2pxfx/PkjuGh8HKm5Nt7+fS9DooO47OTHIHkB7qVP8sGA8byx8Q0i/SMx683c+/u9xAXGcfWwq5nVbxYWQwMJYbwyBUUimGC/uhfJWc2Gw5LBLRrB57VsVXQT6ZYj+O3ZB9hpXcuYSsnZkx/ruIb0nQYH1x/+4Pf9Dm9NgaxNcN476M56hetH3cycGXMorizmPxv+w9xdc9mevx2P0Y/BfU/hlMAR3FxQxmO5hfzbfxJzT57Df1NSGVrpYmPiPxEtnDCsl8TzYcBp8Nu/CCjP4OYT+7Fidx7frMtgziIrWVtu5uZCwW6zmQf1W3hg5f1klGYcdZkAs4F3rxzHxePjeH3ZHu76ejOVbsieehdLB03nytTPuGXJLRRWFvLk8U/ywzk/cFa/s44y7gC7DpXSy30Qi9tGSVDjmiMJwQn8e8q/+eHsH5jRewafJH/CzUnn82t0IfOtAezcMbdeSYCa7C7czW1Lb+PqRVdT4C7m7lw7886dz0nxJ7UuzWHMKMjfw5gofb2RNPM2HiTIYmCia63mD/dGZrUYgwkCex7WhT8ikiYlS7tHR/sfHUFTHyMjR/LghAdZeXAl/9v8P0AbwR+WDW6lH37Vf7WItlOfruWecnskL/26k36RAZw3RgsAuO/0wZw4MJLHftjOmqJADk24hhtzfuO/G/7Lyb1PZu6suXx/9ve8Mu0VQswh/OvPf3Hat6fxzpZ3KHHUM/L2Rri4/CLqjWqzWgyUVif9aEFe1rI8bUW5j+iWI/jnl9+HW0j+HncJwtC65emtoo+Wp5W0lZCzA357CsL7a+6KHocfgY+LPo4lFyzB6XHWGYGRnp5Gyoe3MWP7V7h3fIvbI5k/8n9EhEW1fZuFgJkvaZPNC+7k8ou+Yc6Kvdz37RaC/YzMG72F4cn7ueqU93nfcZBPdnzCkv1LuHjwxdw4/EZCLCGApoOyp2gPw4fsYa9czy+5u/jtsxw8QjMsUXrBP4NGcu6sDzA2IiGwOaOIUUJbpdoUA19FlaG/ccSNvLP1HRakLuThwHBIeRvjrvfpH9KfwWGDGRQ2iMFhgxkYOpBAUyDZZdm8sekN5qfOJ8AQwO1jbmfqj69T7InD3Bb3U8xoQDI1MIsfdwZS7nDVmkwvq3SxaFs2546KokfKGhg4o20yJAX3IqRaF76cCX0ODw6qknz09y4ca6qK5IUDL2Rr3lbe2vwWw8KH0Stk7BG68M2QNa7JgT+18Nth52rzYzX4fuNBdufYePOyMRi8kWZ6neDVS0Zzzht/cNO3H+IXswKX2cyTnlDOmfI8QqeVO6X3KZwcfzJrs9fy/rb3eXXjq7y37T0uHHghlw+5nKiAGt+pKpmCBlwomg++hS4arw5NzSeTtqbbGfhlB5axwb6d64oqGXvZXR3bmLjjwGDRlu/bi2HYeTDrVTAfrfNi1BvrNXRxcQmY7/iGB99+h/NLPmaR/1ncO+sCVq9c4Zt2h8TBKY/Dz/fglzKXp8+ZxvJdudw1zkj4R1fA4DOxJp7PP4DZg2bz5uY3+Sz5M77f/T2DwgaRWpRKYeVhkbVgczB9I+NIzYgh3BTPI6dO56RN72Pa8TPMKARrjwabsyW9mAmmfUhTIOX+DYRs1kNCcAJPT36a3p6rsSRdT2h4AbuPu5yUghSS0pOYt2deddle1l7klOcgkVwx5AquH349IQZ/nN/ez7LQpsffN4h3onW4bh9SjiA5q5SxvQ8L2i3ecYgKp5vLYg5h2lrcaORRkwmOw5K5ESGODpVMzioh1N9IkC1VSz0ZHF/PRWojhODh4x5mV+EuHlrxEO+d+glZMgy3MKBvzUTr5i+1J5dZr9Wa7HW4PLyyZBeJsUGcPqy2m81kdDFmzFJ+Tf+eSns8X/Y/jwG//UsTDawROSWEYEL0BCZETyA5P5kPtn3Axzs+5tPkTzmr71n8Y8w/iPCLqDbwhsD670/NB99CF41Xh4aACGjbvDTVdCsDX+Ys4+lVj9Pf4SBGfyZ6UwcnkjZaNDXKvcvhjBdgwg0tDr3sEWjhwVtv5uVfT+KCsb0wG9pmIUe9jLsOtn4Dix7g1L+v5dShifDJuaAzwt9eqC4WFRDFE8c/weVDLufNTW+SU5HD9Pjp9A/pT/+Q/gwIHUC4JRwhBCt353Hzp+t5/Gs7A8/7B323fwd//BdOe7rBpmzOKOJm015E7GhNvrmFTOgTybxfRvJE7kfQ73wY1xspJbkVuaQUpLCzYCfJBclMipnE9cOvr87sVZmxGTMuXJFtMDcD2qRzUCxx9p3ACHZkFtcy8N9tPEhsiB9Di5fgEQZ0/We0Tb0hcYiUBfS0mo5y0SR7J1hF7k5NelfXdO+txWDhlWmvcPGCi3lg5V0E+F9HkSGK8NYY+KzNED3iqMHQl2sPkFFYwdPnDq/lNtlZsJP7f7+f1OJUToq+iPnLRvC6JYb/RM5FLH5UC62s4yloSPgQnj/xeW4rvY2Ptn/EvN3z2FW4iw9P/xCL1wfvF1z/CNtqMVDa0sTbVU8I/srAN4nXNr5GTmUBT+aWkXb81R3dHI1z3gKHDcL6tPpSQRYjj89qIyPTGDodnPUqvD1FmxgdcCrsXQZ/e7HOSaEBoQN4ZforDV5y8oAIvrl5Eld/sIaL5uaSNPA8rGvfg+P/oUWX1EGFw82BnALizXshtnUTjcNigvinzptcPG0FhPZGCEEP/x708O/B1F51+7nzU9cTA1ji2nBZR8xo/HK3EuJ/aa2J1pxSOyt353LLiX0RKQsoCB1JuKWNFsEEx4HbwfAwe60RvNsj2XmolMuO6w27d2mL/5pJjDWGl6a9xA2/3oA55ksO5rXcwAuPW3Npjr++1v5yh4tXl+5hQp8wpg7QjK6Uks9TPufldS8TbA5mzow5TIqZxCDTHp5ftJNTJtzKWVv+DmvegeNvrbfOuMA4Hpn4CMfHHM/ty27nidVP8JTNTrG0Eh5Uf2y72aCnUl9l4Js5gi/3hvIGREJR805tKt1mkjWtMo3Pkz/nohIbW8qnMW5I65bGtxnWyDYx7h1Cj8HawpJtc2HBndBrPIy7tlWXHBIdxGfXH4fb4+H6tOlItwP++E+95XdkFTNY7kMv3a0OXTXqdQTGJVIoQrQJ7yZSkb6ZSmkkqk8b/rhGj0Lk72FcT32tUMn5mzLxSJgdXwJF+8mLaCO3EECI5nYZ5l9MRo3FTvvyyrA7PQyP1ENxepP970cyvud47p9wP3bjVt4LkNDCxU7+5RngskN07R/UD/5II89Wyf2nD0IIgdPj5N7f7+XZNc9yfMzxzJ01l0kx2o/TLSf246yRMfxjbSj50SdqcsJl+XVVV4uT4k/i76P+zoK9C/iwaHuDq1ir22s2Uanza74P3vuEQICKg28Qp8fJl/lfEilM/L3QxjzzLAb0aP0qNgUw+U6IHALuSjjrv22i8dG/RyAfXTuBrRUR/GqYhlz3PpRk1Vl2c3oxo3VeWeD6crA2g3F9wlnpGoxn73KqxYAawZi7g52yF317hLS6/mpitJR/JwZmkZJdissrvfv9poMMjw0m/tBvgCAvohWLm44kuCoWvoisEnt14u+qCdbhZk3uuMUSHcDFgy4mwXQSvwXl8wvlLYoNt9q8ss81DHxxuZO3l6dy8uAejO0dhtPt5N7l9/JL2i/cOfZOXj3pVcIshyeNhRA8f/4IhkYHcW3m2UiHrcGBRE1uHHEjp8SfwqvuLFb4BdCjEQMfYDZg17Ug6UdNF42P6BYG/pMdn3DQeZAHcw6RJKcysP+A1oWxKQ5jMMEV38HVP7fN+gAvI3qF8M6V43i+/Cw8LgfO5S/VWW5LRhETzfs041SPG6c5jE8IZZVnGDpbNuTvafwEKQmz7eSAsS9+pjac94gZBcBowz4cLg+puWXsPlTKtoMlnDs6VtNDip+I0xTSdnV6FzvF6/KRkurcBMlZJRh0gni3V6O8FQZeCMHfYm4hqCKMf0aGs/NA81cOW217weivRZx5efv3VErsLu4+dRBOt5N7lt/D0gNLuX/8/VybeG2d33c/k545V44jwxDPGt1I3DubJoWlEzqenvw0fdzwVlQlHn1uw+01GygT/s038OVVI3hl4OuluLKYtza/xUQZysm2Uv5rP4NJfX33yHNMEhQD8ce1+WWP7x/BfZeewbfuqcj1H1JZkH5UmS0ZxYzWpUJsE6UAGmF0fCh/Sa8aZFPcNKXZWN3FFAW1zG1RLwERENSL+EotReKOrGLmbTyIXic4O8EJh7a2TnumLsyBYAmhh9Ri3av88MlZJfTvYcVYsBv0ptYtTgLiQ4MJyvgbgR4P/1j/HIX25qWtDCzdC1GJ1U+LOaV2PvgjjVkjYxgQ5cfdy+/mt/TfeGDCA1w+tGFRsdgQP/53+ViSHIPR5++C0kNNaoO/0Z+XckuQCN7Y8QhlzrL622sxUIZ/y1w0Jis/JRcxb7dvZlm7vIEPNgfzzon/4fHMvRyIOoU0Gc3x/ZSB7yqcNqwnATMeQEgPf370CO4aS+iLK5wU52US4cpuvXSEF6vZQEDPAeTpI5tk4D1ZWwFwt1UETU1iRhFYsA2zQce2gyX8sCmTyf0jCD/wq3a8rcIjaxIcR4hDS2xRFUmTnFV6WKIgvD/oWxd7ERNi4aArgVcP5ZHvLOXu5Xc3aWEZAB4PVts+LYLGyxu/7cHh9nDbyX24e/ndLEtfxoMTHuSyIZc16ZIT+oTh6u2ViU5rYmix20k/exEnZg7gYNkBHlzxIB7pqbNogNkbSdPcSdayPPAPZ9nOHFYcbANZhzro8gYeYOT+tcQ6SvnCdB4xwRZ6hzcg16rodMycOonU2LOZWLSAF75eUp3gefvBYka1of+9inF9wljhHIJMW1GtiVMfpfu1VHH+8T4QRo0ZhShIZUwPHd9uyOBgUYW2MjNlAUQNb/VIuk5C4jCXZWLUC9ILyyksc5BdYmdIdCDkprR4grUmsSF+lGBlgMfMYwFDWZu9lhfWvtD4iQCF+zC4y6GnZuDTC8r5fM0Bzh/bk1e3/pNl6ct46LiHuHTIpc1qU6/BEyiW/thSljbtBG+ES5CzN/eOv5dl6cuqV+seidVsoET6tcxFExBJvq2SQJNvXMpd38C7KmH1mxSEjOTrzEgm9Ws4wYWiczL4oicw6KDXtv/x/C+aQuHmjGJG6VKRQn9UREVrGJ8QxgrXUER5vhaO1wCVB7eQISOIj/GBXoh3ovWkkGyKyp34m/TMiBfaKk5fjN4BguMQRenEBFtILyivnmBNjDRrK08jWm/gI6xmTHodBaZozqpwcNXQq/gi5Qu+3fVt4ydnb9H+ez/vV5bsQujcFAa8Q1JGEg8f9zCXDL6k2W06rn8P/vQMbXr0VA2ZgksHX8o5/c/hrc1vsXj/4qOKWs0Git2WFozgcyEggoIyhzLw9SIlTL6TNVEXUVDmYJJyz3RNQuLRjbmCi43LmZ/0J28vT9UmWE37EFHDwNR2T2XjEkJZ7fG6XBr5whvzdpDs6U0/X0RlRWsGfqxRkw4+PbEn/vt+BWTb+9+rCIkDRymDQjykF1ZUx+APNR/S6m2DEbxOJ4gOsZCl6wmFadwx9g6Ojzmep/56ik05mxo+OWsLHqGHHkPYdaiUeZv202fIXP48tJJHjnuEiwdf3KI2DewRyGbDcKzlGVB0oPETakS4CCH458R/MiJyBA+vfPgoiWSr2UChx9ICH3w++EeQX+YgSBn4ejBaYOLNLK/Ubkxl4LsuYsrd6HU6notazDMLU1iWkk0ie9rM/15Fj0AL5vA4DhliGzbwzgqCy9LYq08gPKCNFDtrEhAOwfH0d6di1AsuHh+vRc+E9G7TiKVaeEMlEwO0WPgdWSX0CDQTUqb9yLQmgqYmsSF+pHkioWg/BnQ8P/V5YgJiuGPZHWSXZdd/YtZmygLiwWDmhV+2EdDrMzKdG/jnxH8ye/DsFrdHpxNU9poM0LScAN4YdX2gpkNj0pv4z7T/EGgM5PZlt9eaOLZaDBS6LOAsg6ZKEdfQodFG8M3rT1Pp+gbeS3K+m4Rwf2JDOlieQNFyQuIQY67khNJFnNfXTaz7IH6esjb1v1cxLiGMFa4hyP1/1K9bnpOMDg/FQYN85/aLGUlQ/la2Pn4aE6INWv6AIWe1LptYQ1TpwpuKyC9zsGF/oXeCNUWTgQjv1ybVxIT4sasyQtNaKc0i2BzMqye9SoWrgv9b+n+8tvE13t36Lp/u+JS5u+ayYO8ClqQtZmX+VlYFRfPd9tWsLH0JEZDMo5Me5aJBF7W6Tb2HjCVPBlG+c1njhb0jeHPw4dDcSP9I/jP9P+SW53LP8nuqJ46t5hbIFVSWgNuB0xJGucPtMxdNt5AqcLk97Cx0c84Y38WTKtqJyXchNnzM85G/8FtAX9hNm4/gQYuHT9o4hAtYAtmb6w7DPLQNABmVePSxtiJ6FCT/iMVVCqlLNYPoK/cMVIuIxevzgZ6k5ZdzemK0ZuDD+rSNaiWagd9cEQImtBWtwbH0C+nH81Of55E/HuGdLe8gqWOhWagFyIB1N2Kwwv3jHuHCgRe2SZsm9ovgT89QTkpb0Wj+VFdpDlLqCQypbVOGRw7n0UmP8sgfj3DLklt46cSXamvC20vAL7SOKx6B9wmhVK/lTlYGvgG2Z5ZQ4VLumW5BcCyMvRrD2vc4te80MAdD+IA2r2Z8QhgveLwpDvf9XqeBr8zYjFNaCIltukRxs/FOtJK1GZIXtD41X2MERIDBQg9PDqC5gYZEB0LqrjZzzwD0CvFjfk3Z4N5ayswT405kxcUrkFLi8Diwu+zan9uOfc9S7L/cz9YBN/J5YTzSHcjlw1ruljmSAT2sfGMYwZn2P6Fgb4NPK/aiHGwEERl4dEKQs/ufjUTyxOonuGLhFZwX8+jhrE5NnWj1RukUiRAA5YNviFWp2pulFjh1EybfBTqDNqKNHd0sZcOm0iciABkQSZa5T71+eMfBLeyUcfTr4ZuM98BhA5/+F+xe3PrUfI0hhKYL7zjsBx/awwIFqW0ywVpFTIgfB2UkElFnflYhBGa9mWBzMFEBUfQO6s2g0jxGVjqJDzmB9IMJjIka3mbtqarTGa/Fw8u9yxss6yrN0VL11SNTcE7/c5gzYw4F9gLe2n072X5ew95UF43XBZSPdm8FGpWBr5dVqXnEWkWjokCKLkJQ9GFRMx/430H7so9LCOUP9zDYvxpcR6wklBJzfjLJnnjfRNBU4R+miYCteQccpb51z1QRHIfJdhA/ox6TQUcf3SHwuNp0BB8b6ocTA+V+PZuegDtrM4T355DDQqndxdCYtv9h7Td4BJkyjPKdvzVcsCy3QQMPmrjaZ3/7DKsxiNReS/kxoBmLnbwumly3JoccZFYGvk4cLg/r0goZHOZjfXRF+zL5Ts1t4qt4cDQ3za/lA8BVAQfX1T5YdACTq5RdIoG4UB9P3EePgrIcLcFF3xZmQGoOIXGI4nTiw/wZFBWIoUCTSyCi7VxR0cGaa6PQFNMMA78FokdwoFRbfDYsJrjN2lPFpH7hrPYMRX/gjwbF5vQVeeTRuJJk76DePDPxHYz2aB7qEcHr+36oXqjXIF4dmmynNnhQcfD1oNcJPr1+AifHd2BqPkXbExgFN/x22IXhA8YnhPGnZ4jmRjjSTeOdYC0OGlSdFs5nVPVxQBul5muM4Hgoy+XR0/vw8MwhmkQBok0NvMWoJ8JqIlsX1bTcrOUFUHwAokdyoMSDTsCgqKMzn7WWfpFWthpHYnEUQE5yveXMlQUUyKAmhcf2tIah238J55TaeDt7Bff/fj+V7sqGTyrLo9JkZXPJVvyiFvFt0WfN7UqT6BYGfmzvMGKsXb4rinZmaEwQTmMwWf6Djjbw2dvwILRFVr6mKkqoRlo5n+INlTwh0s7EvuFaBE1IfJsuJgNvLLw7EmzZ4ChvuHDVCtae2gi+b6S1bdU7vQghcPeuioevxw/vKMPkqaDCFNqkH3ctL2sgT+YVcEfYOBamLeTaX64lryKvVjkpJbsKd/HR9o+4OW8FJ8SE8VvRvzCErqDQXYjL0/Z6ND6PohFC6IF1wEEpZTs4GBWKpmHU6xgdH8Kq/KFckLFAM0JeI+fJ3sp+GUWvqPoTLrcZCVPgqgVaesf2wLvYieJ0iBigjeDb0P9eRUyIHykZ3sCHogO1Es0fRdZhiYIDJX9ywiDfTWwPHjyMtNQoIncuI2DiLUcX8PrHneamBW0EmA3YMSHRc51fAvHTLuWhFQ9x2U+X8e8p/ybTlsnqzNWszlpdbfT7YuQCjx9rdTdSVtyb2xKMGHRtb47bI0zydiAZ8GEogkLRMsYlhLFg3wAuMDm0SJZ+0wFwZW5hhyeefj3qT9fWZggBfab4vp4qvCN4itK1lZd5u6v73ZbEhvixdWco6NH88A0a+M0QHEcRVvLt0icTrFVM7Kv54c9PX6X1/8ioJa+Bl03UaTcZdN7UfQH4VZYyo/cMYgJiuPW3W7l60dUAhJhDmBQ9iUkx2l/PTy6AoFjOKRxCpNUAVDRYR0vxqYEXQvQCZgJPA3f5si6FoiVMSAjjXc8gLbn1vt81Q1dZiqlkP8meCZwa2Q0zgwXGgNBpI/jCNC1bl49G8POcEYcNfENkb4GeIw5r40T7zsD3jQjgffNILnEt0+o9cp7HG8KoD+zR5GtazQYq9AH4ecMkh0UM48uZX7I8YzmJEYkMDhuMTtRw95TlQ/RICjIcJIT70yUNPPAf4D6g3tkSIcSNwI0AUVFRJCUltagim83W4nO7KqrPrafCJbELC2nGfkRsXsAGw4kEFSczBkiW8SSmbKRgT8epk/rqM55oCqNo51pyC4wMB9anl1Fa3Lb1FB5ykU8QTp2Z7G0rSbXX/SOid1UwOW83adax/LhCk2cu2reNpEzfve951iFQDHuWfEBGfHGtY1FZKxkC5NncTX7v9dJJiduE5+A+ttU4pwc9yMnKIYecw4WlZKoth4y8cnKKyykvdGCzOn3yOfvMwAshzgRypJTrhRDT6isnpZwDzAEYN26cnDat3qINkpSUREvP7aqoPrcNw5JXsqlyDOfZvmLaxNGwdQ9shHzrQE4/pe1dF83BZ5/x3gH0FE56RhlhG4w99RKwtO2oOeJgMa9tXInd2pu4ADdx9fXjwF+wUtJn0tlUbookxJzJrNN8+75n+h1g90+xxLoO0P+IdlUsWw87of+w0Uyb1rRJ9sgtK6gsDybBamr887IXw3IX0QNHYd8NIwf3wyoyfPI5+zL05ARglhAiDfgSOEkI8akP61MoWsS4hFB+KO4P0qMtesrehk0EYI1M6Oim+Y7gOM0Hn7dLc9m0sXEHzUUDUNBYLHyNCJodWSXEB/k+Im5Sv3BWeYZiPPgnuGtnm7IXZVMmzYSGhDT5egFmA7amJv3w+vhtBk2zJswXSqVeGn0nhRBD69g3rbHzpJQPSil7SSkTgIuB36SUDSdQVCg6gAkJYfzp7IdHb4J9vyMPbWOHrzTgOwshcVByUIv3b0OJgpqE+hvxM+rJ0kVpBr6+BUBZm8A/ArtfFHtybMQH+t7AJ4T7s8MyCqO7Ag5uqHXMWdKwTEFdBDYnq5NXh6ZYaAu5OtTAA18LIe4XGn5CiNeAZ3zWIoWinRmbEEolJrKCRsHeJGT2Nra747q3gQ+OA+mG7G0+mWAFLeY8JsSi6cI7yw8n0TgS7wrWPblluDyyXUbwQgj0fSbjQRwdD1+WSz7B9GiGgQ8wGyhuatKPI3RofJJrwEtT3snjgDhgFbAWyERzvzQZKWWSioFXdFZ6BFpICPdnDYmQsx2dq4Idsjf9u2METRVVoZJIiPSdWmZsqD+7HN548rpWtLoc2orS6JHsyNSMY3uM4AFGDOhLsieeil219eF15XnkySAirUcrSdaH1WKg0O2nadE0JlVQpUPj0e6vcKvvVi835Z10osXw+AEWYJ+U9aQXVyi6KOMSwvihqG/1ts9Fxjoary484LMRPEBsiIWtZSHaRl1++Nxk8Dir/e8BJj09/Nsnaknzww/DnLUOnPbq/ebKAopEMEF+TY9BCTQbKHSbtb647A0X9urQZDm0+6ujXTRr0Qz8eGAKcIkQ4huftUih6AAmJISxsiIej8EfN3oyTQnNekTvcgT3OvzahwY+JtivYQNfYwXrjswShkQHofNVNqsjiA/zJ8UyCr3HARlrtJ1S4u8swG4MbVYWrwCzgQK3V5SuMUXJsjwwWcmxC4x6QZDFd9HqTTHw10kpH5VSOqWUWVLKs4H5PmuRQtEBjEsIxYWBzPCJpBn7ERfZvC94l8PkD/4RWoIR/zCfVRMb6ocdMy7/qHoM/GYwBeIJSWBHVolPV7AeiRACU9/JuNAd1oe3F6HHjcPSvNwSVrOBUlkjq1NDlHlzsdochPqbfHqfNWrgpZTr6tj3iW+ao1B0DH0iAggPMPFmyN38nQfp153971WE9wdfpiPkcKhkWUAvLXXfkWRvgZ7DSS+yY6t0+XQFa12MHhDPVk9f7LuTtB1lWoSL9G9e+k9NcKxqBN+IgS/PA/8I8sscPnXPQDdQk1Qo2oKqBCCL99pJKTV3b/97Fee/A2e/4dMqYr0GvtAUffQI3uOG7K21JljbcwQPh+PhzYc2QaWtOsJF1wyZAtAmWQ+n7WtsBJ+rjeDLKgm3KgOvULQL4xPCyC3VdLyPiRF8SLyWA9eH9Ay2oBOQpesJxRm1M2flp2rhk9HaBKteJxjoAw34hugV6scuvzHopAsO/ImrVJMUMAY108DXHME36qLJ9xp4B+EBvp3nUQZeofAyPuGwL7p/e6hIHgMY9TqigiyaLjxSEzirImuz9t87gu8XGYDF2L6Z2YQQ+PU7HicG5N7llBdquWr9Q3o26zoBZgMlNCHxtpTaCF65aBSK9mVoTBB+Rj16nSA+TBn4tiImxI+dlVWx8PsOH8jeDHozRAzUJljb2f9exdj+MWzw9KdyTxIVRZqBDwyPatY1Ai1eqQJo2EVTWQIeJy6/MErtLp8ucgJl4BWKaox6HeMSQukXGYDJoL4abUVMiB9byjTdlVqLnbK2QNRQCuySrGK7T3KwNoWJfcNZ5R6GOXcrMm83hdJKRHDzXHQBtSZZGxjBV+nQ6LWnxTAf++DbI+GHQtFlePb8EVQ43B3djG5FbIgfv2yzIC1mRNVEq5Sai2bYOR02wVpFXJg/uwPGIBzfEn7wN/bLYCKbubrUajbgwoBTZ8FoL66/oNfAF+t8L1MAagSvUNQiNsSP/sdCBE07EhtiweEWuIPiDkfSFKeDvci7glUziEM6yEUDENTvOCowYXSVkU/zhMZAM/AADn1AwyN47yrWAqn1NUxNsioUiq5MbGiNWPgqA189wTqKHZklRAdbfD7h2BDj+0ez1q2pahbrgps92avXCfyMeir0AQ374L0j+ByPFi2kwiQVCkWXJqY6Fj7msA8+awsIPUQN7dAJ1iom9gtntUdL7lFhbNnKXqvFQIVoZATvjbPPcnqFxpSLRqFQdGWqDHym6AmVxVBRqK1gjRyEHROpuWUd5n+vIjbEj73WMQA4zC0z8IFmA2XCv+E4+PJ8MFnJtWuj/iCLsUV1NRVl4BUKhU8JshgJtBjY74nUdhSmaS6aniPYdagUt0d2+AgeIKz/eL5xTWVf2OQWnV8dSdOYi8a7yCnU34RO51u9I2XgFQqFz4kN8WNnpXdknL4WSrMgegTbvRE0HRUiWZMJ/Xtwr+tmKiJHtuh8TXDMv3EXjX8E+TaHz90zoAy8QqFoB2JD/Nha7o2FT/aK0XpXsAaaDfTyTsR2JJP6RiDEYf2c5hJgNlDk8WvERXN4BN8ek8rKwCsUCp8TE+LHnmIB/uGw/w9tZ8/h7MjyasD72FXRFHoGW5h78yQunhDfeOE6CLQYKHKbwVEKnnpyItXUofFxBA0oA69QKNqBmBA/iiucuIN7g/RAaB88piCS21kDvjHG9g6rjmlvLtaaST8cdbhpaujQ5NkqlYtGoVB0Dw7HwntzwUaPYH9BOeUOd6eYYG0LAswGCpzehUt1+eG9OjRuv3BK7C6fL3ICZeAVCkU7EBuiJbAuNEZrO3qO6HCJgrYm0GKg0ONN1F2XH967yKnMoM1F+FqHBpSBVygU7UBsiCalm6XzqjRGj2JHVjEGnWBAVPeQhggw6RvO6uQ18EWifXRoQBl4hULRDkQGmjHoBOvME2DMVdD7eLZnltC/hxWzoX014H2F1WKskdWpDhdNlQ4NysArFIpuhF4niA6xsLs8AGa9CiZ/dmR2rgnW1mI1GyitzupUh6KkdwR/yN0+OjSgDLxCoWgnYoL9yCyqACC3tJKc0spuM8EKNRY6Qd0jeK8OzSGX5pLq0pOsQgiLEGKNEGKzEGK7EOIJX9WlUCg6P7Ehfhws1Ax8clb3mmAFb+JtGki8Xa1DI9AJCPHzrQ4N+DbhRyVwkpTSJoQwAiuFEAullH/6sE6FQtFJiQ31I7vEjsvtYUeVge9WI3g95ZiR6BB1juC1Vaz57aRDAz4cwUsNm3fT6P2TvqpPoVB0bmJC/PBIyC6xsyOzhNgQP0L8O04Dvq2xmo2AwGkIqCdMUlvkVGBrH5kC8HHKPiGEHlgP9AfekFL+VUeZG4EbAaKiokhKSmpRXTabrcXndlVUn7s/3am/+XkuAH5atpq1eyrpGaCrs29dtc8VLm38WibNFOzfRcoRfRibk0alOYK9BTnoBbX66Ks++9TASyndwCghRAgwTwiRKKXcdkSZOcAcgHHjxslp06a1qK6kpCRaem5XRfW5+9Od+tsrx8aL65YTGDuA7LVbuWhiX6ZNG3hUua7aZ49HwpKfcZlD6BkaQM8j+7C+ksD4Qbj2+DGkZxDTpo2pPuSrPrdLFI2UsghYBpzeHvUpFIrOR4x3NetvKTlI2b0mWAF0OoHVbKBC+B8dJllDh6a9lCTBt1E0kd6RO0IIP2AGkOKr+hQKRefG32QgLMDEyj1auGB3mmCtIsCs17I6HTnJWqVD4x9BUbmz6xt4IBpYJoTYAqwFFkspF/iwPoVC0cmJCbFgd3oItHQODfi2xmo2YMP/6DDJah2aEKB9FjmBD33wUsotwGhfXV+hUHQ9YkP82HZQS7ItRMdrwLc1VouR0gq/o0fwXgNfLLTMVd1hBK9QKBS1qErA3d3871VYzXqK68rqdJQOje9XsYIy8AqFoh2pSofXGXKw+gKr2UCx2wLuSnBVHj7glSnIaUcdGvBxmGRb4HQ6ycjIwG63N1guODiY5OTkdmpV56Az9tlisdCrVy+MRt8vw1Z0PQb3DEInYHR8SEc3xSdYzcbDWZ0qS8HgHal7XTTZ1To0ysADkJGRQWBgIAkJCQ367EpLSwkMDGzHlnU8na3PUkry8/PJyMigT58+Hd0cRSfkhP7h/PngyfQIsnR0U3yC1awn3+U16vZiCIjQXpfngymQXLtACAhtpxW8nd5FY7fbCQ8P75YTMt0NIQTh4eGNPm0pjl2EEN3WuIMmOJZfV9q+slwICKegrJIQPyP6dkoy3ukNPKCMexdCfVaKY5kAs4EiWUdWp7K8dl/kBF3EwCsUCkVXINBsoFTW8MFXUZ4HAZHk2RyEW9snggaUgVcoFIo2o5YmvP2IEXxAOAVljnZJ1VeFMvAKhULRRgSYDNiOdNFIqVw0XZGKigpOPPFE3G43RUVFvPnmm9XHHA4HU6dOxeVydWALFQpFe2K1GLBxhIH36tB4/CMoLFcj+C7D+++/z3nnnYderz/KwJtMJk4++WS++uqrBq+RlJTE1Vdf7eOWKhSK9sBqNuDAiFtnOuyiqaFDI2X7xcBDF4iDr8kTP25nR2YdmVIAt9uNXq9v9jWHxgTx2FnDGiyzdu1arrvuOtasWYPb7WbChAl89dVXfPbZZ3z++ecAPPDAA6SmpjJq1ChmzJjBCy+8wDnnnMODDz7IZZdd1ux2HcnmzZu57bbbyMvLIyUlBSkl9913H88991yrr61QKNoGq1kzqU5jIPqqSVavgS/ReXVo2nGStUsZ+I5i/PjxzJo1i0ceeYSKigouv/xyBg4cyN69e0lISADg2WefZdu2bWzatKn6vMTERNauXdvq+u12O7Nnz+bjjz9mwoQJ/POf/8Rut/Pwww+3+toKhaLtsFo0k+rQW7FUuWi8OjT5MhiwEaFG8HXT0Ejb16s6H330UcaPH4/FYuHVV1/l0KFDhISENHiOXq/HZDLV2bbjjjuOyspKbDYbBQUFjBo1CoDnnnuO0047rVbZJUuWMGbMGCZMmADAiBEjWLRokYo5Vyg6GVUjeLven6DqEbymQ5PrCQRshLWTDg10MQPfkeTn52Oz2XA6ndjtdvz8/Jq0YrOyshKL5eiVe3/9paWnTUpK4sMPP+TDDz+s9xrbtm1j+PDh1dsbNmxgzJgx9ZZXKBQdg59Rj05AhQg4ygef7QoA2tcHryZZm8hNN93Ev/71Ly677DLuv/9+QkNDcbvd1UY+MDCQ0tLaGtD5+flERES0WngrPDycLVu2ALBr1y6+++47Lr744lZdU6FQtD1CCALMBspFjaQfNXRooP10aEAZ+Cbx8ccfYzQaufTSS3nggQdYu3Ytv/32G6eeeiorV64ENCN8wgknkJiYyL333gvAsmXLmDlzZqvrv+SSS7DZbCQmJnLjjTfyxRdfEB4e3urrKhSKtifQbMBWM21ftQ6Ng2A/I0Z9+5ld5aJpAldeeSVXXnkloPnVq9wrISEhvPLKK5xyyikA1RE1VXz++ec8++yzDV572rRpjWZTt1qt/Pjjjy1svUKhaE+sFu9ip5ouGv8I8tt5FSuoEXyrGDNmDNOnT8ftdh91zOFwcM455zBw4MAOaJlCoegoAswGiqWf5qKRslqHJt9W2W6JPqpQBr6VXHvttXXG35tMpupRv0KhOHawmg0Uuv0ACQ5bLR2a9pxgBWXgFQqFok2xmg0UuauSfpQcoUPTfoucQBl4hUKhaFOsZgP5Lm9odMnBGjo0TuWDVygUiq6M1WI4nLavYC8AFcZQ3B6pXDQKhULRlbGaDeQ4ahv4Yl0QgJpkVSgUiq6MtWZWJ6+B13RoILy7+OCFEHFCiGVCiB1CiO1CiNt9VVdHkZWVxZlnngnApk2b+Pnnn6uPLViwgEcffbTZ19y0aRMTJ05k1KhRjBs3jjVr1gCapMGqVauqy1199dV8//33reuAQqFoc6wWA6XSm9XJa+DzPJoWVXdy0biAu6WUQ4GJwN+FEEN9WF+78/LLL3PDDTcARxv4mTNn8uOPP1JeXt6sa95333089thjbNq0iSeffJL77rsPONrAKxSKzonVXCPpR8E+ALJcVqD9XTQ+W8kqpcwCsryvS4UQyUAssKPFF134AGRvrfOQn9sF+hZ0p+dwOKP+1ab1acEnJiby7bff8tRTT+FwOHj00UepqKhg5cqVPPjgg8yePZtp06axYMECLrrooiY3RwhBSYm2Aq64uJiYmBjS0tJ466230Ov1fPrpp7z22msA/PHHH/zvf/8jOzub559/ngsuuKD5/VcoFG2KZuAtSASiogBMgeR1gA4NtJNUgRAiARgN/FXHsRuBGwGioqJISkqqdTw4OLhaxMvsdKBz15MCT4KrvmMN4HE6qDxCJKwmgwcP5rTTTuO+++6joqKCCy+8kN69e7N161aCg4NxOBwAPPTQQ2zYsIGXXnoJ0OSLhw0bxtKlSznjjDOa3J6nn36ac889l7vvvhuPx8PixYsJDw/nmmuuwWq18o9//AMAp9NJUVERCxcuZNeuXcyePfsomeGOwm63H/U5thU2m81n1+6MHGv9ha7f590FbiQ6HDo/zJ5yKnQBbNm5Fz8DrFr5e53n+KrPPjfwQggr8C1wh5TyqHRMUso5wByAcePGySN1WZKTkw9rqc96ud56WqMH39hv6tNPP12tBV81ki4tLSUqKqq6TovFgslkqtWG3r17s3Dhwma16+OPP+Y///kP559/Pl9//TW33347S5YswWw2Yzabq69lNBo588wzCQ4OZvz48eTm5vpUD785WCwWRo8e7ZNrJyUlNard05041voLXb/PEQeLYc1KPOYgqCjHL7wXloBIosqK6+2Xr/rs0ygaIYQRzbh/JqX8zpd1+ZIqLfjS0tJqeeCm6MFX6cY3xDXXXMOoUaP429/+BsBHH33EeeedB8CFF15YPclaF2bz4Rl5KWWT+qJQKHxLYFVWJ4PmdycgkoIyB+HtmKqvCl9G0QjgPSBZSln/0LsLcKQWPMDAgQNJS0urLlOXHvyuXbtITEwEYN68eTz44INHXfuDDz6oNUEbExPD8uXLAfjtt98YMGBAvddXKBSdjwBvVqdKfYB3R8fo0IBvR/AnAFcAJwkhNnn//ubD+nxCfVrwAQEB9OvXjz179gAwffp0duzYwahRo/jqq6+A2nrwqampBAUFNVrfO++8w913383IkSN56KGHmDNnDgBnnXUW8+bNY9SoUaxYscJHvVUoFK2lKm1fhc5r4DtIKhh8G0WzEujySUPr04IHuPXWW/nwww956qmnCAsLq5Vg+9ChQ1RUVFSn2tu0aROvvPJKo/VNnjyZ9evXH7V/4MCB1VmdAKZMmVJrRG+z2ZrfOYVC0eaYDToMOkEFWiy8xz+Cwg4awauEH63g3HPPJT8/v85jBw4cqI6oAfj000/bq1kKhaIDEUJoST+EZuArTKG4OkCHBpSBbzXXX399nfvHjx/fzi1RKBSdBavZQKl3BF8svDIF7bzICZQWjUKhULQ5VrOBEo8mGVwoOkaHBpSBVygUijbHajZQKLVJ1hyPZuC7WxSNQqFQHJMEmA0s1k2B894hyxMKKBeNQqFQdAusFgNZTj8YcREFZZWAGsF3OXwhF3zvvfcyePBgRowYwbnnnktRUVGt4wcOHMBqtfLiiy8CUFRUxJtvvll9PCkpqbpNCoWiYwg0Gyir1LSx8mwOrGYDZoO+3duhDHwr8IVc8IwZM9i2bRtbtmxh4MCBPPPMM7WO33XXXbXEy4408AqFouMJMBuw2TUD31GrWKGLhUk+t+Y5UgpS6jzmdrvR65v/Czk4bDD3T7i/3uPtLRd86qmnVr+eOHEic+fOrd7+/vvv6dOnDwEBAdX7HnjgAVJTUxk1ahQzZsxg5syZ2Gw2LrjgArZt28bYsWP59NNP0ZQjFApFe2A1GyhzuPF4pFeHRhn4Tsn48eOZNWsWjzzyCBUVFVx++eUkJiayb98+QkNDqwW/nnzySdatW8frr79efe64ceNYsWJFswx8Td5//31mz54NaCtVn3vuORYvXlztngF49tln2bZtG5s2bQI0F83GjRvZvn07MTExnHDCCfzxxx9Mnjy5he+AQqFoLlWCY2UOF/llDmJDLB3Sji5l4BsaabdGLrgxHn300Wq54FdffRXQ/O+RkZENntejRw8yMzNbVOfTTz+NwWDgsssuA+Dxxx/nzjvvxGq1NnruhAkT6NWrFwCjRo0iLS1NGXiFoh2pEhyzVbooKKtkeGzjOlS+oEsZ+I6iSi7Y6XRit9sJCAhoU7ngjRs3EhMTU+3D//DDD1mwYAFLly6tdq389ddfzJ07l/vuu4+ioiJ0Om36pK4sTjVlhPV6PS5X8xOhKBSKllMlOGazu7w++PZf5ATKwDeJKrngffv2cf/99/P666+3SC54zZo1R02afvDBB7W2Fy1axPPPP8/y5cvx9/ev3l9TQfLxxx/HarVy00034XA4lIywQtHJqDLwmcV2nG7ZIUqSoKJoGqW95YJvvfVWSktLmTFjBqNGjeLmm29usHx4eDgnnHACiYmJ3Hvvva3srUKhaAusXh/8gQItik5F0XRS2lsuuOoHoyEef/xxgOqR++eff17reM3UXzUnfRUKRftQNYI/kF8GdMwqVlAGvlUouWCFQlEX1QbeO4LvCKExUAa+1Si5YIVCcSRVBn5/vtdF00EjeOWDVygUijYm4KgRvDLwCoVC0S0wGXSYDDrKHW78TXosxvbXoQFl4BUKhcInBHpH8R0VQQPKwCsUCoVPqHLThFs7ZoIVlIFvEnq9nlGjRjFs2DBGjhzJSy+9hMfjqT6+Zs0apk6dyqBBgxg9ejTXX399s1UkuxpKxVKhaJiqidaO8r+DMvBNws/Pj02bNrF9+3YWL17MwoULeeKJJwAt3v3CCy/kueeeY+fOnWzcuJHTTz+9268uVQZeoWgYq3LRNJ/Zb6/mm3XpADjdHma/vZp5GzMAqHC4mf32an7crAl8ldidzH57NYu2ZQGaLvPst1ezZMchAHJKG9aSqYsePXowZ84cXn/9daSUvPHGG1x11VVMmjSpuswFF1xAVFRUrfPcbjf33HMPiYmJjBgxgtdeew2ApUuXMnr0aIYPH861115LZaWW/SUhIYEHH3yQUaNGMW7cODZs2MBpp51Gv379eOuttwBNvmDq1KnMnDmTQYMGcfPNN1c/WXzxxRcMHz6cxMRE7r//sEib1Wrl4YcfZuTIkUycOJFDh7T3Ijc3l/PPP5/x48czfvx4/vjjD0BbVHXttdcybdo0+vbtWy22VlOmWK2gVSiOpmo1a7ccwQsh3hdC5Aghtvmqjo6ib9++uN1ucnJyqjXXG2POnDmkpaWxadMmtmzZwmWXXYbdbufqq6/mq6++YuvWrbhcLv73v/9VnxMfH8+mTZuYMmUKV199NXPnzuXPP//kscceqy6zZs0aXnvtNXbs2EFqairfffcdmZmZ3H///fz2229s2rSJtWvX8v333wNQVlbGxIkT2bx5M1OnTuWdd94B4Pbbb+fOO+9k7dq1fPvtt7Xi+1NSUvjll19Ys2YNTzzxBE6nk2effZZ+/fqxadMmXnjhhTZ6ZxWK7kNnGMH7cqHTh8DrwMdtedGvbjo8UjbqddXbpaWl+Jn0tY4HWYy1tsMCTLW2ewS2n0bzkiVLuPnmmzEYvB96WBibN2+mT58+DBw4EICrrrqKN954gzvuuAOAWbNmATB8+HBsNhuBgYEEBgZiNpurU/lNmDCBvn37AnDJJZewcuVKjEYj06ZNq5Yzvuyyy/j9998555xzMJlM1Sn9xo4dy+LFi6vbt2PHjur2lpSUYLPZAC07ldlsxmw206NHj+pRv0KhqJ+A7mzgpZS/CyESfHX9jmTv3r3o9Xp69OjBsGHDWL9+PWeffXab11Ml+6vT6WpJAOt0umoJ4CMzNTWWucloNFaXqSkl7PF4+PPPP7FYjv7RU/LDCkXzqUr6EdGBUTQdLlUghLgRuBEgKiqKpKSkWseDg4ObNGHpdrt9OrFZde28vDyuv/56brjhBmw2G1dffTXTp09n2rRp1fIE8+fPZ+LEifTo0aP6/ClTpvDGG28wbtw4DAYDBQUFxMTEsG/fPjZt2kS/fv14//33Oe644ygtLUVKic1mw2w2Y7fba8kCVx3zeDysWbOGrVu3Eh8fz2effcY111zD0KFDue2220hLSyMkJIRPP/2Um266qfr8qv8VFRU4nU5KS0uZPn06L774IrfffjsAW7ZsYcSIEVRWVmI0GqvP8Xg82Gw2rFYrJSUldb7ndrv9qM+xrbDZbD67dmfkWOsvdJ8+52Y6ANi3cysiu+GFTr7qc4cbeCnlHGAOwLhx42RNJUSA5OTkJmVq8mVGp4qKCqZMmYLT6cRgMHDFFVdw1113odPpCAwM5KuvvuK+++4jJycHnU7H1KlTOffcc2vpud96660cOHCAE044AaPRyA033FCtRnnNNdfgcrkYP348d9xxB2azGSEEVquVwMBALBYLJpOpun9Vx3Q6HePHj+eBBx5gz549TJ8+nUsvvRSdTsdzzz3HWWedhZSSmTNncvHFF1e3peo6fn5+GI1GAgMD+d///sff//53TjjhBFwuF1OnTuWtt96qds1UnaPT6bBarSQkJDB58mQmTZrEGWecUcsPb7FYGD16tE8+i6SkJI68R7ozx1p/ofv0OdWwj3l7djBjyiTiwvwbLOuzPkspffYHJADbmlp+7Nix8kh27Nhx1L66KCkpaVK57sRPP/0kZ86c2dHNOIqmfmYtYdmyZT67dmfkWOuvlN2nz+kFZfKlX3dKj8fTaNnW9BlYJ+uxqR0+glcoFIruSK9Qf+6aMbBD2+DLMMkvgNXAICFEhhDiOl/VdawyZcoUFixY0NHNUCgUnRRfRtFc0obXajQ6RNE50J4YFQpFZ6DTr2S1WCzk5+crw9EFkFKSn59fZ6ilQqFofzq9D75Xr15kZGSQm5vbYDm73X7MGZbO2GeLxUKvXr06uhkKhYIuYOCNRiN9+vRptFxSUpLPQvM6K8dinxUKRdPp9C4ahUKhULQMZeAVCoWim6IMvEKhUHRTRGeKThFC5AL7W3h6BJDXhs3pCqg+d3+Otf6C6nNz6S2ljKzrQKcy8K1BCLFOSjmuo9vRnqg+d3+Otf6C6nNbolw0CoVC0U1RBl6hUCi6Kd3JwM/p6AZ0AKrP3Z9jrb+g+txmdBsfvEKhUChq051G8AqFQqGogTLwCoVC0U3pcgZeCHG6EGKnEGKPEOKBOo6bhRBfeY//1dUTfzehv3cJIXYIIbYIIZYKIXp3RDvbksb6XKPc+UIIKYTo8iF1TemzEOIi72e9XQjxeXu3sa1pwr0dL4RYJoTY6L2//9YR7WwrhBDvCyFyhBDb6jkuhBCvet+PLUKIMa2utL5UT53xD9ADqUBfwARsBoYeUeb/gLe8ry8Gvurodvu4v9MBf+/rW7pyf5vaZ2+5QOB34E9gXEe3ux0+5wHARiDUu92jo9vdDn2eA9zifT0USOvodreyz1OBMdSTxhT4G7AQEMBE4K/W1tnVRvATgD1Syr1SSgfwJXD2EWXOBj7yvp4LnCy6braQRvsrpVwmpSz3bv4JdHWt3qZ8xgD/Ap4D7O3ZOB/RlD7fALwhpSwEkFLmtHMb25qm9FkCQd7XwUBmO7avzZFS/g4UNFDkbOBjqfEnECKEiG5NnV3NwMcC6TW2M7z76iwjpXQBxUB4u7Su7WlKf2tyHdoIoCvTaJ+9j65xUsqf2rNhPqQpn/NAYKAQ4g8hxJ9CiNPbrXW+oSl9fhy4XAiRAfwM3NY+Teswmvt9b5ROrwevaBpCiMuBccCJHd0WXyKE0AEvA1d3cFPaGwOam2Ya2lPa70KI4VLKoo5slI+5BPhQSvmSEGIS8IkQIlFK6enohnUVutoI/iAQV2O7l3dfnWWEEAa0R7v8dmld29OU/iKEOAV4GJglpaxsp7b5isb6HAgkAklCiDQ0X+X8Lj7R2pTPOQOYL6V0Sin3AbvQDH5XpSl9vg74GkBKuRqwoIlydVea9H1vDl3NwK8FBggh+gghTGiTqPOPKDMfuMr7+gLgN+mdweiCNNpfIcRo4G00497V/bLQSJ+llMVSyggpZYKUMgFt3mGWlHJdxzS3TWjKff092ugdIUQEmstmbzu2sa1pSp8PACcDCCGGoBn4hnN3dm3mA1d6o2kmAsVSyqzWXLBLuWiklC4hxK3AL2iz8O9LKbcLIZ4E1kkp5wPvoT3K7UGb0Li441rcOprY3xcAK/CNdy75gJRyVoc1upU0sc/diib2+RfgVCHEDsAN3Cul7KpPpk3t893AO0KIO9EmXK/uwoM1hBBfoP1IR3jnFR4DjABSyrfQ5hn+BuwByoFrWl1nF36/FAqFQtEAXc1Fo1AoFIomogy8QqFQdFOUgVcoFIpuijLwCoVC0U1RBl6hUCi6KcrAKxQNIIQIEUL8X0e3Q6FoCcrAKxQNE4KmUKpQdDmUgVcoGuZZoJ8QYpMQ4oWOboxC0RzUQieFogG8CWMWSCkTO7otCkVzUSN4hUKh6KYoA69QKBTdFGXgFYqGKUWTKFYouhzKwCsUDeBVbPxDCLFNTbIquhpqklWhUCi6KWoEr1AoFN0UZeAVCoWim6IMvEKhUHRTlIFXKBSKbooy8AqFQtFNUQZeoVAouinKwCsUCkU35f8B1HvWLg3DVdsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Create a dc sequence with added Gaussian noise.\n",
    "n, c = 32, 5.\n",
    "x = c + np.random.standard_normal(n)\n",
    "\n",
    "# Setup the k highest frequency components to remove.\n",
    "kths = [8, 24]\n",
    "\n",
    "# Store the inverse DFT associated with each k.\n",
    "xks = []\n",
    "\n",
    "for k in kths:\n",
    "\n",
    "    # Transform the sequence to frequency domain.\n",
    "    xfreq = np.fft.fft(x, n)\n",
    "\n",
    "    # Remove the kth high frequency components.\n",
    "    xfreq[-k:] = 0.\n",
    "    \n",
    "    # Transform the data back to the time domain.\n",
    "    xinv = np.fft.ifft(xfreq, n)\n",
    "    xks.append(xinv)\n",
    "\n",
    "# Plot approximate solution.\n",
    "t = np.linspace(0, 1, n)\n",
    "plt.plot(t, x, label='x(t) + $\\sigma$')\n",
    "for (k, xk) in zip(kths, xks):\n",
    "    plt.plot(t, xk.real, ls='-', label='x(t), -{0}th'.format(k))\n",
    "plt.axhline(c, linestyle=':', label='DC component')\n",
    "plt.title(\"Removing High Frequency Noise\")\n",
    "plt.ylabel('x')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Demonstrate use of FFT to perform 1D convolution of a random signal with a 4-point kernel $[1, 2, 2, 1]$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "nx1, nx2 = 32, 4\n",
    "x1, x2 = np.random.random(nx1), np.array([1,2,2,1], dtype=np.float64)\n",
    "\n",
    "# Zero pad the window to equal size of x1.\n",
    "x2pad = np.hstack((x2, np.zeros(nx1-nx2)))\n",
    "\n",
    "# Compute the DFT of the input sequences.\n",
    "x1fft, x2fft = np.fft.fft(x1, nx1), np.fft.fft(x2pad, nx1)\n",
    "\n",
    "# Compute the pointwise product in frequency domain.\n",
    "x1x2freq = x1fft * x2fft\n",
    "\n",
    "# Transform back to time domain.\n",
    "x1x2 = np.fft.ifft(x1x2freq, nx1)\n",
    "\n",
    "# Compute convolution only where signals overlap eg mode='valid'.\n",
    "expected = np.convolve(x1, x2, mode='valid')\n",
    "\n",
    "# Compare only the overlapping part of the sigals.\n",
    "np.testing.assert_almost_equal(x1x2[nx2-1:], expected)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12.04 Wavelets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Wavelet basis functions are improvement over sine and cosine for representing functions that change abruptly (eg sharp corners).\n",
    "\n",
    "**Mother Wavelet** (def.) Single function $\\phi(x)$ generated by dilation and translation, i.e., $\\phi((x-b)/a)$, where a and b are real numbers with $a \\neq 0$.\n",
    "\n",
    "Sparse Wavelet Basis \n",
    "* Wavelets are localized in frequency (by dilation) and time (by translation).\n",
    "* Fourier basis are localized in frequency but not time.\n",
    "\n",
    "#### Discrete Wavelet Transform\n",
    "* Analogous to the FFT.\n",
    "* Computed by pyramidal, hierarchical algorithm.\n",
    "* DWT requires $O(n)$ work compared to $O(n \\log n)$ for FFT.\n",
    "  * Sparsity of wavelet basis functions make it faster than FFT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary: Fast Fourier Transform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Sines and cosines are appropriate for representing cyclic phenomena.\n",
    "* DFT transforms input sequence to linear combination of sines and cosines.\n",
    "* FFT performs DFT in $O(n \\log n)$ work, provided input $n$ is power of 2.\n",
    "  * Compare to $O(n^2)$ for matrix-vector product form of DFT.\n",
    "* Wavelets are more effective for rapidly changing functions than DFT.\n",
    "* Discrete Wavelet Transform aka DWT requires $O(n)$ work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}