{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4ec40081b048ce2f34f3f4fedbb0be10",
     "grade": false,
     "grade_id": "cell-98f724ece1aacb67",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "# CDS: Numerical Methods Assignments\n",
    "\n",
    "- See lecture notes and documentation on Brightspace for Python and Jupyter basics. If you are stuck, try to google or get in touch via Discord.\n",
    "\n",
    "- Solutions must be submitted via the Jupyter Hub.\n",
    "\n",
    "- Make sure you fill in any place that says `YOUR CODE HERE` or \"YOUR ANSWER HERE\".\n",
    "\n",
    "## Submission\n",
    "\n",
    "1. Name all team members in the the cell below\n",
    "2. make sure everything runs as expected\n",
    "3. **restart the kernel** (in the menubar, select Kernel$\\rightarrow$Restart)\n",
    "4. **run all cells** (in the menubar, select Cell$\\rightarrow$Run All)\n",
    "5. Check all outputs (Out[\\*]) for errors and **resolve them if necessary**\n",
    "6. submit your solutions  **in time (before the deadline)**"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "team_members = \"Koen Vendrig, Kees van Kempen\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "2900c663e4345c7f0707be39cc0cc7f4",
     "grade": false,
     "grade_id": "cell-b6b5c93e38567117",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "## Linear Equation Systems\n",
    "\n",
    "In the following you will implement the Gauss-Seidel (GS), Steepest Descent (SD) and the Conjugate Gradient (CG) algorithms to solve linear equation systems of the form \n",
    "\n",
    "$$A \\mathbf{x} = \\mathbf{b},$$ \n",
    "\n",
    "with $A$ being an $n \\times n$ matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "55d777356d474cb7327bb087dfe5e644",
     "grade": true,
     "grade_id": "cell-bcf2dd8f9a194be8",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as linalg\n",
    "from matplotlib import pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "9717dea99ca91d38963fc25b2ef95e03",
     "grade": false,
     "grade_id": "cell-595bb99f65f79ca7",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 1\n",
    "First, you need to implement a Python function $\\text{diff(a,b)}$, which returns the difference $\\text{d}$ between two $n$-dimensional vectors $\\text{a}$ and $\\text{b}$ according to \n",
    "\n",
    "$$ d = || \\mathbf{a} - \\mathbf{b}||_\\infty = \\underset{i=1,2,\\dots,n}{\\operatorname{max}} |a_i - b_i|. $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "437477eff36ed8c083b4903a6921591c",
     "grade": true,
     "grade_id": "cell-bcd26d8a9f0e6447",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def diff(a, b):\n",
    "    return np.max(np.abs(a - b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "a10d859fa1b26c11ae00c7b846fdd1e8",
     "grade": false,
     "grade_id": "cell-19ae2e6fe5c9264c",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 2 \n",
    "\n",
    "The Gauss-Seidel iteration scheme to solve the linear equation system \n",
    "\n",
    "$$A \\mathbf{x} = \\mathbf{b}$$\n",
    "\n",
    "is defined by \n",
    "\n",
    "$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right].$$\n",
    "\n",
    "Note especially the difference in the sums: the first one involves $x_j^{(k)}$ and the second one $x_j^{(k-1)}$.\n",
    "\n",
    "\n",
    "Give the outline of the derivation in LaTeX math notation in the markdown cell below. (Double click on \"YOUR ANSWER HERE\" to open the cell, and ctrl+enter to compile.) \n",
    "\n",
    "Hint: Similar to the Jacobi scheme, start by seperating the matrix $A$ into its diagonal ($D$), lower triangular ($L$) and upper triangular ($U$) forms, such that $A = D - L - U$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "3c0c1bc83919a390739fa5802fc75743",
     "grade": true,
     "grade_id": "cell-3a8403fee85d9582",
     "locked": false,
     "points": 2,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "source": [
    "---\n",
    "\n",
    "We start from our linear equations:\n",
    "\n",
    "$$Ax = b$$\n",
    "\n",
    "We separate A into different components (diagonal, strictly lower triangular and strictly upper triangular):\n",
    "\n",
    "$$A = D - L - U$$\n",
    "\n",
    "\n",
    "We write $D - L$ as $L'$ to get:\n",
    "\n",
    "\n",
    "$$(L' - U)x = b$$\n",
    "\n",
    "\n",
    "We take the iterative process of the Gauss-Seidel method to write:\n",
    "\n",
    "\n",
    "$$\n",
    "L'x^k = b + Ux^{k-1}\\\\\n",
    "x^k = L'^{-1}(b + Ux^{k-1})\\\\\n",
    "$$\n",
    "If we write every component of the matrix $A$ as $a_{ij}$, we can use forward substitution to rewrite our previous equation to:\n",
    "\n",
    "\n",
    "$$x^k _i = \\frac{1}{a_{ii}}\\left[-\\sum_{j=0}^{i-1}a_{ij}x_{j}^{k} -\\sum_{j=i+1}^{n-1}a_{ij}x_{j}^{k-1} + b_i\\right].$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "dfbc2086b4294a7d942d35d9e1a22d5b",
     "grade": false,
     "grade_id": "cell-4ebd4c774509912d",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 3\n",
    "\n",
    "Implement the Gauss-Seidel iteration scheme derived above\n",
    "\n",
    "$$x_i^{(k)} = \\frac{1}{a_{ii}} \\left[ -\\sum_{j=0}^{i-1} a_{ij} x_j^{(k)} -\\sum_{j=i+1}^{n-1} a_{ij} x_j^{(k-1)} + b_i \\right],$$\n",
    "\n",
    "where $a_{ij}$ are the elements of the matrix $A$, and $x_i$ and $b_i$ the elements of vectors $\\mathbf{x}$ and $\\mathbf{b}$, respectively.\n",
    "\n",
    "Write a Python function $\\text{GS(A, b, eps)}$, where $\\text{A}$ represents the $n \\times n$ $A$ matrix, $\\text{b}$ represents the $n$-dimensional solution vector $\\mathbf{b}$, and $\\text{eps}$ is a scalar $\\varepsilon$ defining the accuracy up to which the iteration is performed. Your function should return both the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. \n",
    "\n",
    "Use an assertion to make sure the diagonal elements of $A$ are all non-zero. Initialize your iteration with $\\mathbf{x}^{(0)} = \\mathbf{0}$ (or with $\\mathbf{x}^{(1)} = D^{-1}\\mathbf{b}$, with $D$ the diagonal of $A$) and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "65b01c4e4c650dc6971b7d6dde62044b",
     "grade": true,
     "grade_id": "cell-ec2390c3cafa2f8e",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def GS(A, b, eps, k_max = 1000000):\n",
    "    \"\"\"\n",
    "    Return the Gauss-Seidel algorithm estimate solution x to the problem\n",
    "    Ax = b and the number of iterations k it took to decrease maximum\n",
    "    norm error below eps or to exceed iteration maximum k_max.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Assert n by n matrix.\n",
    "    assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
    "    n = len(A)\n",
    "    \n",
    "    # First we decompose A = D - L - U.\n",
    "    D = np.diag(np.diag(A))\n",
    "    U = -np.triu(A) + D\n",
    "    L = -np.tril(A) + D\n",
    "    \n",
    "    # We need non-zero diagonals elements.\n",
    "    assert np.all(np.diag(D) != 0)\n",
    "    \n",
    "    x_prev = np.zeros(n)\n",
    "    x_cur  = np.dot(linalg.inv(D), b)\n",
    "    \n",
    "    k = 1\n",
    "    while diff(x_cur, x_prev) > eps and k < k_max:\n",
    "        k += 1\n",
    "        # We will have to copy, as the array elements will point to the same\n",
    "        # memory otherwise, and changes to one array will change the other aswell.\n",
    "        x_prev = x_cur.copy()\n",
    "        for i in range(n):\n",
    "            x_cur[i] = 1/A[i, i]*(-np.dot(A[i, :i], x_cur[:i]) - np.dot(A[i, i + 1:], x_prev[i + 1:]) + b[i])\n",
    "    return x_cur, k"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "badd9e6110882657c17c42e3d9ad8370",
     "grade": false,
     "grade_id": "cell-78d8cad45f5350b9",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 4\n",
    "\n",
    "Verify your implementation by comparing your approximate result to an exact solution. Use $\\text{numpy.linalg.solve()}$ to obtain the exact solution of the system\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\begin{pmatrix}\n",
    "        10 & -1 & 2 & 0 \\\\ \n",
    "        -1 & 11 &-1 & 3 \\\\\n",
    "         2 & -1 & 10&-1 \\\\\n",
    "         0 &  3 & -1& 8\n",
    "    \\end{pmatrix} \\mathbf{x}^*\n",
    "    =\n",
    "    \\begin{pmatrix}\n",
    "    6 \\\\\n",
    "    25 \\\\\n",
    "    -11\\\\\n",
    "    15\n",
    "    \\end{pmatrix}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Then compare you approximate result $\\mathbf{\\tilde{x}}$ to the exact result $\\mathbf{x^*}$ by plotting $|| \\mathbf{x}^* - \\mathbf{\\tilde{x}}||_\\infty$ for different accuracies $\\varepsilon = 10^{-1}, 10^{-2}, 10^{-3}, 10^{-4}$. \n",
    "\n",
    "Implement a unit test for your function using this system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "b54e0939053f9afee0f0ba46e490ba17",
     "grade": true,
     "grade_id": "cell-0585da22f1d304ab",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEKCAYAAADXdbjqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAQxUlEQVR4nO3df2id133H8c9nipOoGZVL445YiWcXBdGwEAzCgXgEN1mmFKLEhKxLukFLTEwysrF/RG3S9b/WYWJjLQlkLg1uyvLDpEbYiVcxSF2HEDYr0aiTGXUmECI5zMkWqU25UNv97g9dx7IiyffqPI+ee+95v0D4Pue599xvOFw+Oec897mOCAEAkOL3qi4AAND+CBMAQDLCBACQjDABACQjTAAAyQgTAECyy6ouoApXX311bNy4seoyAKCtvPHGGx9GxLrFzmUZJhs3btT4+HjVZQBAW7H97lLnslrmsj1ke+/s7GzVpQBAR8kqTCLiUETs7OnpqboUAOgoWYUJAKAchAkAIFlWG/C2hyQN9fX1VV0KAKyq0YlpjYxN6tRMTevXdmt4sF/bN/cW1n9WMxP2TADkaHRiWrsPHNf0TE0haXqmpt0Hjmt0Yrqw98gqTAAgRyNjk6qdOXdRW+3MOY2MTRb2HoQJAHS4UzO1ptpXgjABgA63fm13U+0rkVWY8KVFADkaHuxX95qui9q613RpeLC/sPfIKkzYgAeQo+2be7Xn3hvVu7ZbltS7tlt77r2x0Ku5sro0GABytX1zb6HhsVBWMxMAQDkIEwBAMsIEAJCMMAEAJMsqTLg0GADKkVWYcGkwAJQjqzABAJSDMAEAJCNMAADJCBMAQDLCBACQjDABACTLKkz4ngkAlCOrMOF7JgBQjqzCBABQDsIEAJCMMAEAJCNMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAECyrMKE26kAQDmyChNupwIA5cgqTAAA5SBMAADJCBMAQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMkIEwBAMsIEAJCMMAEAJCNMAADJCBMAQLKOCBPbX7L9lO0XbT9SdT0AkJvKw8T207ZP235rQfudtidtn7S9a7k+IuJERDws6auStpZZLwDg0yoPE0n7JN05v8F2l6QnJX1F0g2SHrB9g+0bbb+04O8L9dfcLellSYdXt3wAwGVVFxARR21vXNC8RdLJiHhHkmw/L+meiNgj6a4l+jko6aDtlyU9u/C87Z2SdkrShg0bivsPAABUHyZL6JX03rzjKUk3L/Vk29sk3SvpCi0xM4mIvZL2StLAwEAUVCcAQK0bJk2JiCOSjlRcBgBkqxX2TBYzLem6ecfX1tuS2B6yvXd2dja1KwDAPK0aJsckXW97k+3LJd0v6WBqpxFxKCJ29vT0JBcIALig8jCx/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE21XWCQBYWuV7JhHxwBLth8VlvgDQFiqfmawm9kwAoBxZhQl7JgBQjqzCBABQjqzChGUuAChHVmHCMhdQvtGJaW19/BVt2vWytj7+ikYnkr8ihjZQ+dVcADrH6MS0dh84rtqZc5Kk6Zmadh84Lknavrm3ytJQsqxmJgDKNTI2+UmQnFc7c04jY5MVVYTVQpgAKMypmVpT7egcWYUJG/BAudav7W6qHZ0jqzBhAx4o1/Bgv7rXdF3U1r2mS8OD/RVVhNXCBjyAwpzfZB8Zm9SpmZrWr+3W8GA/m+8ZIEwAFGr75l7CI0NZLXMBAMqRVZiwAQ8A5cgqTNiAB4ByZBUmAIByECYAgGSECQAgGWECAEiWVZhwNRcAlCOrMOFqLgAoR6FhYvuzRfYHAGgPhd1OxfYjks7ZvjUi/rKofgEAra/Ie3Odqf97tsA+AQBtoMhlrvcl9Uo6XWCfAIA2UGSY3CzpNUmbCuwTANAGCguTiPi2pP+T9FBRfQIA2kOhV3NFxJsRMVNkn0XieyYAUI6mN+Btb2jwqTMR8atm+y9TRBySdGhgYIDZEwAUaCVXc/2ogeeEpH2SnllB/wCANtN0mETEl8soBADQvrJa5gIAlINlLgBAsqbCxHYXy1wAgIWWvTTYdq/t79r+TL3p720/UT/3vdKrAwC0hUt9z+RGSTdJGqgffyTp3frjX9s+ZPsqSbI9aPu1csoEALSyZZe5IuKntr8eEUfrTVsk/Vv93Ldsf03SEdu/lfSxpF2lVgsAaElN7ZlExN22PydJtm/X3K1TfiPpGkkPRsRk8SUCAFpdI7dTeXP+QUR8VH/4mKS/i4htku6T9ILt24otDwDQDi45M4mIkSXab5v3+Ljtr0j6iaRbiiuvWLaHJA319fVVXQoAdJQi7xr8vqTbi+qvDPwGPACUo+i7BteK7A8A0B6Sw6S+dAQAyFgRM5PvFNAHAKCNFREmLqAPAEAbKyJMooA+AABtrNANeABAnggTAECyIsLkfwroAwDQxpLDJCLuKKIQAED7YpkLAJCMMAEAJCNMAADJmgqTeb+q+PvllAMAaEfNzkw+Z/tRSX9cRjEAgPbUbJjcLukbkr5o+wvFl7Nytq+yPW77rqprAYDcNBsm/yHpQUnvRsTpIgqw/bTt07bfWtB+p+1J2ydtN/Lb8t+UtL+ImgAAzWn2N+BP1B/+osAa9kl6QtIz5xtsd0l6UtIdkqYkHbN9UFKXpD0LXv+gpJsk/ZekKwusCwDQoKbCpAwRcdT2xgXNWySdjIh3JMn285LuiYg9kj61jGV7m6SrJN0gqWb7cET8rsy6AQAXNLzMZbvX9ndtf6Z+/D3bZd1+vlfSe/OOp+pti4qIxyLibyU9K+kHiwWJ7Z31PZXxDz74oOh6ASBrzeyZ3Ki55aSB+vGvJR2cd7nwoO3XCq6vKRGxLyJeWuLc3ogYiIiBdevWrXZpANDRGl7mioif2v56RBytH3/L9tckHbH9W0kfS2pko7wR05Kum3d8bb0NANCCVvwNeNu3S3pI0m8kXS3pbyLi1YLqOibpetubbF8u6X5JB1M7tT1ke+/s7GxygQCAC5oNkzfnPX5M0rcjYpuk+yS9YPu2Zguw/Zyk1yX1256yvSMizkp6VNKYpBOS9kfE2832vVBEHIqInT09PaldAQDmcUQxv7pr+xpJP4mIWwrpsEQDAwMxPj5edRkA0FZsvxERA4uda/rSYNsbljm9Y975mYj4VbP9l8n2kKShvr6+qksBgI7S9MzE9s+WOR2SXP93X0Q8s8xzK8PMBACaV+jMJCK+nF4SAKCTFL3MNV/LLXMBAMqxktup/KiB54Tm7rnVUstc7JkAQDkKu5qrnbBn0jlGJ6Y1MjapUzM1rV/breHBfm3fvOSddwAkKHTPBGgVoxPT2n3guGpnzkmSpmdq2n3guCQRKMAq4zfg0bZGxiY/CZLzamfOaWRssqKKgHxlFSbcTqWznJqpNdUOoDxZhQm3U+ks69d2N9UOoDxZhQk6y/Bgv7rXdF3U1r2mS8OD/RVVBOSLDXi0rfOb7FzNBVSPMEFb2765l/AAWkBWy1xswANAObIKEzbgAaAcWYUJAKAchAkAIBlhAgBIRpgAAJIRJgCAZFmFCZcGA0A5sgoTLg0GgHJkFSYAgHIQJgCAZIQJACAZYQIASEaYAACSESYAgGRZhQnfMwGAcmQVJnzPBADKkVWYAADKQZgAAJIRJgCAZIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEiWVZhwOxUAKEdWYcLtVACgHFmFCQCgHIQJACAZYQIASEaYAACSESYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJgCAZIQJACAZYQIASNYRYWJ7m+1XbT9le1vV9QBAbioPE9tP2z5t+60F7XfanrR90vauS3QTkj6WdKWkqbJqBQAs7rKqC5C0T9ITkp4532C7S9KTku7QXDgcs31QUpekPQte/6CkVyPi57b/QNI/SvqLVagbAFBXeZhExFHbGxc0b5F0MiLekSTbz0u6JyL2SLprme4+knTFYids75S0U5I2bNiQWjYAYJ7Kl7mW0CvpvXnHU/W2Rdm+1/Y/S/qx5mY5nxIReyNiICIG1q1bV2ixAJC7ymcmRYiIA5IOVF0HAOSqVWcm05Kum3d8bb0tie0h23tnZ2dTuwIAzNOqYXJM0vW2N9m+XNL9kg6mdhoRhyJiZ09PT3KBAIALKg8T289Jel1Sv+0p2zsi4qykRyWNSTohaX9EvF1lnQCApVW+ZxIRDyzRfljS4SLfy/aQpKG+vr4iuwWA7FU+M1lNLHMBQDmyChMAQDkIEwBAsqzChEuDAaAcWYUJeyYAUI6swgQAUA7CBACQjDABACTLKkzYgAeAcmQVJqkb8KMT09r6+CvatOtlbX38FY1OJN97EgA6QuW3U2kXoxPT2n3guGpnzkmSpmdq2n3guCRp++Ylf2oFALKQ1cwkxcjY5CdBcl7tzDmNjE1WVBEAtA7CpEGnZmpNtQNATrIKk5QN+PVru5tqB4CcZBUmKRvww4P96l7TdVFb95ouDQ/2F1UeALQtNuAbdH6TfWRsUqdmalq/tlvDg/1svgOACJOmbN/cS3gAwCKyWuYCAJSDMAEAJMsqTLidCgCUI6sw4fdMAKAcWYUJAKAcjoiqa1h1tj+Q9O4ip3okLVwDW6ztakkfllDapSxWy2r108hrLvWc5c4vda7Vx0QqZlzKGpNGnlfWuLT7mKy0n07+rPxhRKxb9ExE8Ff/k7S3wbbxVqlvtfpp5DWXes5y55c61+pjUtS4lDUmVY5Lu49JmePSiZ8VlrkudqjBtqoUVctK+mnkNZd6znLnlzrX6mMiFVNPWWPSyPM6cVz4rDReSyGyXOZKZXs8IgaqrgMXMCathzFpTWWNCzOTldlbdQH4FMak9TAmramUcWFmAgBIxswEAJCMMAEAJCNMAADJCJOC2b7K9rjtu6quBXNsf8n2U7ZftP1I1fVAsr3d9g9sv2D7T6uuB3Nsf9H2D22/2OxrCZM620/bPm37rQXtd9qetH3S9q4GuvqmpP3lVJmfIsYlIk5ExMOSvippa5n15qCgMRmNiIckPSzpz8usNxcFjcs7EbFjRe/P1VxzbN8q6WNJz0TEH9XbuiT9UtIdkqYkHZP0gKQuSXsWdPGgpJskfV7SlZI+jIiXVqf6zlXEuETEadt3S3pE0o8j4tnVqr8TFTUm9df9g6R/iYg3V6n8jlXwuLwYEfc18/780mJdRBy1vXFB8xZJJyPiHUmy/bykeyJij6RPLWPZ3ibpKkk3SKrZPhwRvyuz7k5XxLjU+zko6aDtlyURJgkK+qxY0uOS/pUgKUZRn5WVIkyW1yvpvXnHU5JuXurJEfGYJNn+huZmJgRJOZoal3rI3yvpCkmHyywsY02NiaS/lvQnknps90XEU2UWl7FmPyufl/QdSZtt766HTkMIkxJExL6qa8AFEXFE0pGKy8A8EfF9Sd+vug5cLCL+V3P7WE1jA35505Kum3d8bb0N1WJcWg9j0ppWbVwIk+Udk3S97U22L5d0v6SDFdcExqUVMSatadXGhTCps/2cpNcl9duesr0jIs5KelTSmKQTkvZHxNtV1pkbxqX1MCatqepx4dJgAEAyZiYAgGSECQAgGWECAEhGmAAAkhEmAIBkhAkAIBlhAgBIRpgAAJIRJkCLqP+I0X/W//7dNp9PtA2+AQ+0CNv/LenWiHi/6lqAZvF/PkDrOCzpF7b/qepCgGbxeyZAC7B9iyRLuqZ+cz6grTAzAVrDn0n6ZUSc9ZzPVl0Q0Az2TIAWYHuLpB9KCkk1SX8VEW9UWxXQOMIEAJCMZS4AQDLCBACQjDABACQjTAAAyQgTAEAywgQAkIwwAQAkI0wAAMn+H8is0a3Xr+ppAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[ 10, - 1,   2,   0],\n",
    "              [- 1,  11, - 1,   3],\n",
    "              [  2, - 1,  10, - 1],\n",
    "              [  0,   3, - 1,   8]])\n",
    "b = np.array( [  6,  25, -11,  15] )\n",
    "x_exact = linalg.solve(A, b)\n",
    "\n",
    "eps_list  = [1e-1, 1e-2, 1e-3, 1e-4]\n",
    "diff_list = []\n",
    "for eps in eps_list:\n",
    "    x, k = GS(A, b, eps)\n",
    "    diff_list.append(diff(x_exact, x))\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.scatter(eps_list, diff_list)\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xlabel(\"$\\epsilon$\")\n",
    "ax.set_ylabel(\"$||\\\\vec{x}^* - \\\\vec{\\\\tilde{x}}||_\\infty$\")\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# As the three algorithm functions will have the same signature,\n",
    "# it makes sense to only write the test function once.\n",
    "\n",
    "def test_alg(alg, alg_name):\n",
    "    \"\"\"\n",
    "    Check that function alg returns solutions for the example system Ax = b\n",
    "    within the error defined by the same eps as used for the iteration.\n",
    "    \"\"\"\n",
    "    \n",
    "    A = np.array([[ 10, - 1,   2,   0],\n",
    "                  [- 1,  11, - 1,   3],\n",
    "                  [  2, - 1,  10, - 1],\n",
    "                  [  0,   3, - 1,   8]])\n",
    "    b = np.array( [  6,  25, -11,  15] )\n",
    "    x_exact = linalg.solve(A, b)\n",
    "\n",
    "    print(\"Starting with A =\")\n",
    "    print(A)\n",
    "    print(\"and b =\", b)\n",
    "    print(\"We apply the {} algorithm to solve Ax = b.\".format(alg_name))\n",
    "    print()\n",
    "\n",
    "    eps_list = [1e-1, 1e-2, 1e-3, 1e-4]\n",
    "    for eps in eps_list:\n",
    "        x, k = alg(A, b, eps)\n",
    "        print(\"For eps = {:.0e}\\tafter k = {:d}\\t iterations:\".format(eps, k))\n",
    "        print(\"x =\\t\\t\\t\", x)\n",
    "        print(\"Ax =\\t\\t\\t\", np.dot(A, x))\n",
    "        print(\"diff(Ax, b) =\\t\\t\", diff(A @ x, b))\n",
    "        print(\"diff(x, x_exact) =\\t\", diff(x, x_exact))\n",
    "        print()\n",
    "        \n",
    "        assert diff(x, x_exact) < eps\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "d8e54a45261b013b9dd07a401d86401c",
     "grade": true,
     "grade_id": "cell-6524fb6d322a6ea0",
     "locked": false,
     "points": 0,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting with A =\n",
      "[[10 -1  2  0]\n",
      " [-1 11 -1  3]\n",
      " [ 2 -1 10 -1]\n",
      " [ 0  3 -1  8]]\n",
      "and b = [  6  25 -11  15]\n",
      "We apply the Gauss-Seidel algorithm to solve Ax = b.\n",
      "\n",
      "For eps = 1e-01\tafter k = 4\t iterations:\n",
      "x =\t\t\t [ 0.99463393  1.99776509 -0.99803257  1.00108402]\n",
      "Ax =\t\t\t [  5.95250909  24.98206671 -10.98990699  15.        ]\n",
      "diff(Ax, b) =\t\t 0.04749090616931895\n",
      "diff(x, x_exact) =\t 0.005366066491359844\n",
      "\n",
      "For eps = 1e-02\tafter k = 5\t iterations:\n",
      "x =\t\t\t [ 0.99938302  1.99982713 -0.99978549  1.00009164]\n",
      "Ax =\t\t\t [  5.99443213  24.99877578 -10.99900762  15.        ]\n",
      "diff(Ax, b) =\t\t 0.005567865937722516\n",
      "diff(x, x_exact) =\t 0.000616975874427883\n",
      "\n",
      "For eps = 1e-03\tafter k = 6\t iterations:\n",
      "x =\t\t\t [ 0.99993981  1.99998904 -0.99997989  1.00000662]\n",
      "Ax =\t\t\t [  5.99944928  24.99993935 -10.99991498  15.        ]\n",
      "diff(Ax, b) =\t\t 0.000550717702960668\n",
      "diff(x, x_exact) =\t 6.018928065554263e-05\n",
      "\n",
      "For eps = 1e-04\tafter k = 7\t iterations:\n",
      "x =\t\t\t [ 0.99999488  1.99999956 -0.99999836  1.00000037]\n",
      "Ax =\t\t\t [  5.99995255  24.99999971 -10.99999375  15.        ]\n",
      "diff(Ax, b) =\t\t 4.744782363452771e-05\n",
      "diff(x, x_exact) =\t 5.11751035947583e-06\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def test_GS():\n",
    "    \"\"\"\n",
    "    Check that GS returns solutions for the example system Ax = b\n",
    "    within the error defined by the same eps as used for the iteration.\n",
    "    \"\"\"\n",
    "    \n",
    "    return test_alg(GS, \"Gauss-Seidel\")\n",
    "    \n",
    "test_GS()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "f88ad91685d2bd73471acc5b600b5177",
     "grade": false,
     "grade_id": "cell-59514f17e337138f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 5\n",
    "\n",
    "Next, implement the Steepest Descent algorithm in a similar Python function $\\text{SD(A, b, eps)}$, which calculates\n",
    "\n",
    "\\begin{align*}\n",
    "    \\mathbf{v}^{(k)} &= \\mathbf{b} - A \\mathbf{x}^{(k-1)} \\\\\n",
    "    t_k &= \\frac{ \\langle \\mathbf{v}^{(k)}, \\mathbf{v}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)}\\rangle } \\\\\n",
    "    \\mathbf{x}^{(k)} &= \\mathbf{x}^{(k-1)} + t_k \\mathbf{v}^{(k)} .\n",
    "\\end{align*}\n",
    "\n",
    "Initialize your iteration again with $\\mathbf{x}^{(0)} = \\mathbf{0}$ and increase $k$ until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$. Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n",
    "Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "215218647aa6a6a314c093e82a4a24fb",
     "grade": true,
     "grade_id": "cell-a105ce7e1a34ce3c",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def SD(A, b, eps, k_max = 1000000):\n",
    "    \"\"\"\n",
    "    Return the Steepest Descent algorithm estimate solution x to the problem\n",
    "    Ax = b and the number of iterations k it took to decrease maximum\n",
    "    norm error below eps or to exceed iteration maximum k_max.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Assert n by n matrix.\n",
    "    assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
    "    \n",
    "    n = len(A)\n",
    "    \n",
    "    x_cur  = np.zeros(n)\n",
    "    x_prev = np.zeros(n)\n",
    "    \n",
    "    k = 0\n",
    "    while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
    "        k += 1\n",
    "        x_prev = x_cur.copy()\n",
    "        \n",
    "        v     = b - A @ x_prev\n",
    "        t     = np.dot(v, v)/np.dot(v, A @ v)\n",
    "        x_cur = x_prev.copy() + t*v\n",
    "        \n",
    "    return x_cur, k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "f554d4f402c9e1a4cbff1f4e9018e9d0",
     "grade": true,
     "grade_id": "cell-62100bf55800145b",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting with A =\n",
      "[[10 -1  2  0]\n",
      " [-1 11 -1  3]\n",
      " [ 2 -1 10 -1]\n",
      " [ 0  3 -1  8]]\n",
      "and b = [  6  25 -11  15]\n",
      "We apply the Steepest Descent algorithm to solve Ax = b.\n",
      "\n",
      "For eps = 1e-01\tafter k = 4\t iterations:\n",
      "x =\t\t\t [ 0.99748613  1.98300329 -0.98904751  1.01283183]\n",
      "Ax =\t\t\t [  6.01376302  24.84309304 -10.89133793  15.04071202]\n",
      "diff(Ax, b) =\t\t 0.15690696195356324\n",
      "diff(x, x_exact) =\t 0.016996711694861055\n",
      "\n",
      "For eps = 1e-02\tafter k = 6\t iterations:\n",
      "x =\t\t\t [ 0.99983175  1.99716093 -0.99850509  1.00217552]\n",
      "Ax =\t\t\t [  6.00414638  24.97397012 -10.98472385  15.00739201]\n",
      "diff(Ax, b) =\t\t 0.02602988052923294\n",
      "diff(x, x_exact) =\t 0.002839069878101119\n",
      "\n",
      "For eps = 1e-03\tafter k = 9\t iterations:\n",
      "x =\t\t\t [ 0.99991029  1.99994247 -0.99999645  1.00023877]\n",
      "Ax =\t\t\t [  5.99916754  25.00016961 -11.00032515  15.00173404]\n",
      "diff(Ax, b) =\t\t 0.0017340408511579142\n",
      "diff(x, x_exact) =\t 0.00023877424067331177\n",
      "\n",
      "For eps = 1e-04\tafter k = 11\t iterations:\n",
      "x =\t\t\t [ 0.99998551  1.99999065 -0.99999949  1.00003874]\n",
      "Ax =\t\t\t [  5.9998655   25.00002733 -11.00005329  15.00028137]\n",
      "diff(Ax, b) =\t\t 0.0002813662634846281\n",
      "diff(x, x_exact) =\t 3.874139053650083e-05\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def test_SD():\n",
    "    \"\"\"\n",
    "    Check that SD returns solutions for the example system Ax = b\n",
    "    within the error defined by the same eps as used for the iteration.\n",
    "    \"\"\"\n",
    "    \n",
    "    return test_alg(SD, \"Steepest Descent\")\n",
    "    \n",
    "test_SD()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "695745368fdf84d75691832839fb2834",
     "grade": false,
     "grade_id": "cell-3277aa3696a5c87f",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 6\n",
    "\n",
    "Finally, based on your steepest decent implementation from task 5, implement the Conjugate Gradient algorithm in a Python function $\\text{CG(A, b, eps)}$ in the following way: \n",
    "\n",
    "Initialize your procedure with:\n",
    "\n",
    "\\begin{align*}\n",
    "    \\mathbf{x}^{(0)} &= \\mathbf{0} \\\\\n",
    "    \\mathbf{r}^{(0)} &= \\mathbf{b} - A \\mathbf{x}^{(0)} \\\\\n",
    "    \\mathbf{v}^{(0)} &= \\mathbf{r}^{(0)}\n",
    "\\end{align*}\n",
    "\n",
    "Then increase $k$ and repeat the following until $|| \\mathbf{x}^{(k)} - \\mathbf{x}^{(k-1)}||_\\infty < \\varepsilon$.\n",
    "\n",
    "\\begin{align*}\n",
    "    t_k &= \\frac{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle }{ \\langle \\mathbf{v}^{(k)}, A \\mathbf{v}^{(k)} \\rangle } \\\\\n",
    "    \\mathbf{x}^{(k+1)} &= \\mathbf{x}^{(k)} + t_k \\mathbf{v}^{(k)} \\\\\n",
    "    \\mathbf{r}^{(k+1)} &= \\mathbf{r}^{(k)} - t_k A \\mathbf{v}^{(k)} \\\\\n",
    "    s_k &= \\frac{ \\langle \\mathbf{r}^{(k+1)}, \\mathbf{r}^{(k+1)} \\rangle }{ \\langle \\mathbf{r}^{(k)}, \\mathbf{r}^{(k)} \\rangle } \\\\\n",
    "    \\mathbf{v}^{(k+1)} &= \\mathbf{r}^{(k+1)} + s_k \\mathbf{v}^{(k)}\n",
    "\\end{align*}\n",
    "\n",
    "Return the solution vector $\\mathbf{x}^{(k)}$ from the last iteration step and the corresponding iteration index $k$. Implement a unit test for your implementation by comparing your result to the exact solution of the system in task 4.\n",
    "Use $\\text{numpy.dot()}$ for all needed vector/vector and matrix/vector products.\n",
    "\n",
    "How do you expect the number of needed iteration steps to behave when changing the accuracy $\\epsilon$? What do you see?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "f709b7bfc8af2ba2d607f60ed03e1074",
     "grade": true,
     "grade_id": "cell-8ac3d07c9d237e47",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def CG(A, b, eps, k_max = 1000000):\n",
    "    \"\"\"\n",
    "    Return the Conjugate Gradient algorithm estimate solution x to the problem\n",
    "    Ax = b and the number of iterations k it took to decrease maximum\n",
    "    norm error below eps or to exceed iteration maximum k_max.\n",
    "    \"\"\"\n",
    "    \n",
    "    # Assert n by n matrix.\n",
    "    assert len(A.shape) == 2 and A.shape[0] == A.shape[1]\n",
    "    \n",
    "    n = len(A)\n",
    "    \n",
    "    x_cur  = np.zeros(n)\n",
    "    x_prev = x_cur.copy()\n",
    "    r_cur = b - A @ x_cur\n",
    "    v     = r_cur\n",
    "    \n",
    "    k = 0\n",
    "    while diff(x_cur, x_prev) > eps and k < k_max or k == 0:\n",
    "        k += 1\n",
    "        x_prev = x_cur.copy()\n",
    "        r_prev = r_cur\n",
    "        \n",
    "        t     = np.dot(r_prev, r_prev)/np.dot(v, A @ v)\n",
    "        x_cur = x_prev + t*v\n",
    "        r_cur = r_prev - t*A @ v\n",
    "        s     = np.dot(r_cur, r_cur)/np.dot(r_prev, r_prev)\n",
    "        v     = r_cur + s*v\n",
    "        \n",
    "    return x_cur, k"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "0d5b479bee6c0297974c003b4ebd99f3",
     "grade": true,
     "grade_id": "cell-4f5bc81f40ddb3fa",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Starting with A =\n",
      "[[10 -1  2  0]\n",
      " [-1 11 -1  3]\n",
      " [ 2 -1 10 -1]\n",
      " [ 0  3 -1  8]]\n",
      "and b = [  6  25 -11  15]\n",
      "We apply the Conjugate Gradient algorithm to solve Ax = b.\n",
      "\n",
      "For eps = 1e-01\tafter k = 4\t iterations:\n",
      "x =\t\t\t [ 1.  2. -1.  1.]\n",
      "Ax =\t\t\t [  6.  25. -11.  15.]\n",
      "diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
      "diff(x, x_exact) =\t 2.220446049250313e-16\n",
      "\n",
      "For eps = 1e-02\tafter k = 5\t iterations:\n",
      "x =\t\t\t [ 1.  2. -1.  1.]\n",
      "Ax =\t\t\t [  6.  25. -11.  15.]\n",
      "diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
      "diff(x, x_exact) =\t 2.220446049250313e-16\n",
      "\n",
      "For eps = 1e-03\tafter k = 5\t iterations:\n",
      "x =\t\t\t [ 1.  2. -1.  1.]\n",
      "Ax =\t\t\t [  6.  25. -11.  15.]\n",
      "diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
      "diff(x, x_exact) =\t 2.220446049250313e-16\n",
      "\n",
      "For eps = 1e-04\tafter k = 5\t iterations:\n",
      "x =\t\t\t [ 1.  2. -1.  1.]\n",
      "Ax =\t\t\t [  6.  25. -11.  15.]\n",
      "diff(Ax, b) =\t\t 1.7763568394002505e-15\n",
      "diff(x, x_exact) =\t 2.220446049250313e-16\n",
      "\n"
     ]
    }
   ],
   "source": [
    "def test_CG():\n",
    "    \"\"\"\n",
    "    Check that CG returns solutions for the example system Ax = b\n",
    "    within the error defined by the same eps as used for the iteration.\n",
    "    \"\"\"\n",
    "    \n",
    "    # TODO: Do we need to plot the error over epsilon like in task 4?\n",
    "    \n",
    "    return test_alg(CG, \"Conjugate Gradient\")\n",
    "    \n",
    "test_CG()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "65dd782269d01314b02643d0610e3348",
     "grade": false,
     "grade_id": "cell-308595a088486388",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 7\n",
    "\n",
    "Apply all three methods to the following system\n",
    "\n",
    "\\begin{align*}\n",
    "\\begin{pmatrix}\n",
    "0.2&  0.1&  1.0& 1.0&   0.0 \\\\ \n",
    "0.1&  4.0& -1.0& 1.0&  -1.0 \\\\\n",
    "1.0& -1.0& 60.0& 0.0&  -2.0 \\\\\n",
    "1.0&  1.0&  0.0& 8.0&   4.0 \\\\\n",
    "0.0& -1.0& -2.0& 4.0& 700.0\n",
    "\\end{pmatrix} \\mathbf{x}^*\n",
    "=\n",
    "\\begin{pmatrix}\n",
    "1 \\\\\n",
    "2 \\\\\n",
    "3 \\\\\n",
    "4 \\\\\n",
    "5\n",
    "\\end{pmatrix}.\n",
    "\\end{align*}\n",
    "    \n",
    "Plot the number of needed iterations for each method as a function of $\\varepsilon$, using $\\varepsilon = 10^{-1}, 10^{-2}, ..., 10^{-8}$.\n",
    "\n",
    "Explain the observed behavior with the help of the condition number (which you can calculate using $\\text{numpy.linalg.cond()}$). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "3d94dc67b052b82dae655b7ff562c7d2",
     "grade": true,
     "grade_id": "cell-490cd96dc58cbff8",
     "locked": false,
     "points": 4,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEMCAYAAAAvaXplAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA93UlEQVR4nO3dd3hUVf7H8fdJJyEJhCQk9NBJA0LoAhFEUBEUUURUsGFZbPuzgLq2tcuuLpZVdkVEEVRAhQUVFTGggNIJNYAgoaQBCenJzPn9ccOQCglk5s4k39fzzJOZOyWfGcL5zj3n3nOU1hohhBCiLDezAwghhHA+UhyEEEJUIsVBCCFEJVIchBBCVCLFQQghRCVSHIQQQlQixUEIIUQlUhyEEEJU4mF2gIqUUgnA34EdwAKt9arzPSc4OFi3a9fOrrmEEKK+2bhxY4bWOqSq+xxSHJRSs4FRQJrWOrrM9pHAvwB34L9a61cADeQAPkBKTV6/Xbt2bNiwoc5zCyFEfaaUOlTdfY7qVpoDjCy7QSnlDrwDXAFEAhOUUpHAaq31FcDjwHMOyieEEKIMhxQHrXUicKLC5j7APq31Aa11EbAAGKO1tpbefxLwdkQ+IYQQ5Zk55tASOFzmdgrQVyk1FhgBNAHeru7JSqkpwBSANm3a2C+lEEI0QE43IK21XgwsrsHjZimljgFXe3l59bJ/MiGEaDjMPJT1CNC6zO1WpdtqTGu9VGs9JTAwsE6DCSFEQ2dmcfgd6KSUilBKeQE3AktMzCOEEKKUQ4qDUmo+sBboopRKUUrdobUuAaYC3wG7gM+11jtq+bpXK6VmZWVl1X1oIYRwRkV5cPIQpGyA3cvh1OHzP+cCqPqwElx8fLyW8xyEEC5JayjIgtx045KTVvl62W1FOeWfP/otiLv1gn61Umqj1jq+qvucbkC6NpRSVwNXd+zY0ewoQghxltUCeZnVN/YVr1uKqngRBb7NwC8EGodAy15nr/uFnr0e1N4ub0H2HIQQoiZKCmve2Odlgu2UrTLcPKtu4Ku67tsM3O37/b3e7jkIIcRFKzwNWUcgO8X4efpYFQ1/OhRWM7bp6Xe2UW8aAa37VN/w+zQBpRz69i6USxcH6VYSQpxTcT5kH4WsFMg+Ur4InLldVaPfqKnRqDcOhbCYc3zLDwEvP8e/LweQbiUhhGuyFBsNf5WNfmkxyMus/DzfYAhsCQGtSn+2hMBWpT9bQuMw8PBy/PsxgXQrCSFci9UCOanVN/pZR4z7qfDl1jvQaOgDWxoDuBWLQEBL8PQx5S25GpcuDtKtJIQL0tr4Rl9VV8+ZbaePgbWk/PM8fc9+u+94WZlv/WUKgLe/Oe+pHpJuJSFE3SvKhYxkyNhr/Mw6fLbhzz4KJQXlH+/uBQEtKnT1VPjW36ipywzmugrpVhJC1D2tjSN6MvZWuJQWgzOUG/iHGw18eHfociUEti7f3+8bDG6yarEzkeIghDg3SwmcOgTpe842/hml1wvKHOnj6QfBnaBNfwieBCGdIbizcZKWhyzN4mpcujjImIMQdagwBzKTjca/bCE4sb/8GbyNmxuNfvQ44+eZIuDfQr791yMy5iBEQ1KuK2hPmUKQbAwKn6HcISjCaPSDO0Fwl7PXGzUxLb6oWzLmIERDYymBkwerGA+opiuo3cDSItDZKARBEdIV1MBJcRDClZ3pCkqvUAAy94O1+OzjGocZjX/0OAjpcrYQBLSUI4BElaQ4COEqigsg5Tc48DMc2WgUgewyiyfauoK6QOeRpeMBXaBZR+kKErXm0sVBBqRFvWa1wLEtRjH442f4c51xfoByh7BoaHdJaTdQ2aOCGsa0D8L+ZEBaCGehtTE4/MfPRkE4uObspHChURAxGNoPgbYDwSfA3KyiXpABaSGc1anDZ4vBH4mQc9zY3qQtRI2BiCFGUWgcam5O0eBIcRDCkXIzjCJwpiCc/MPY7hdiFIGIIcbeQdN2psYUQoqDEPZUmAOHfj1bDFK3G9u9/I0xg753GwUhtJscNSScihQHIepSSSGk/H52EPnIRmN2UXdvY4WwoU9BRAK06Gn3JSCFuBjy1ynExbBa4Pi2s8Xg0FooyTcmm2vREwY8YHQTte4Lno3MTitEjbl0cZBDWYXDaW1MNfHHz3BglXFEUcEp476QrhB369kjiuTcAuHC5FBWIc4n60j5I4pOHzW2B7aB9oONbqKIweDf3MyUQtSaHMoqRG3knYCDq892FWXuM7b7NqtwRFGEDCKLekuKgxBgHFG0Z7lREI5vBzR4NTa6h+JvLz2iKFKmpBYNhhQHIX75F3z/tLFUZas+cOkTRjFoGQfunmanE8IUUhxEw6U1rHwBVs+AqLEw5m3w8jM7lRBOQYqDaJisVvh2Gvz2vnGE0ag3wc3d7FRCOA0pDqLhsZTAkvth66fQfypc/oIMLAtRgRQH0bCUFMKiO2DXUrj0SRj8qBQGIarglIdeKKX8lFIblFKjzM4i6pGiXJh/o1EYRr4CQx6TwiBENRxSHJRSs5VSaUqppArbRyql9iil9imlppW563Hgc0dkEw1E/in4eKxxVvOYd6DfvWYnEsKpOWrPYQ4wsuwGpZQ78A5wBRAJTFBKRSqlhgM7gTQHZRP1XU46fDTKmARv3IfQ82azEwnh9Bwy5qC1TlRKtauwuQ+wT2t9AEAptQAYAzQG/DAKRr5SarnW2lrxNZVSU4ApAG3atLFjeuHSso7A3DGQlQITFkCny8xOJIRLMHNAuiVwuMztFKCv1noqgFJqMpBRVWEA0FrPAmaBMbeSfaMKl5S5H+ZeY0yMd8tiaDvA7ERCuAynPVpJaz3nfI+RWVlFtVJ3wsfXGGspTFoKLXqYnUgIl2Lm0UpHgNZlbrcq3VZjWuulWuspgYGBdRpMuLiUjTDnSmNNhdu+kcIgxAUwszj8DnRSSkUopbyAG4EltXkBpdTVSqlZWVlZdgkoXNAfq2HuaPAJhNu/hZAuZicSwiU56lDW+cBaoItSKkUpdYfWugSYCnwH7AI+11rvqM3ryp6DKGfPt/DJdRDYCm77Fpq2MzuREC7LUUcrTahm+3JguSMyiHpu+0L48m4Ii4GJi8CvmdmJhHBpTnmGdE1Jt5IAYMOHsOhOY53mW5dIYRCiDrh0cZBuJcEvM+F/D0Gn4XDzIvAJMDuREPWCSxcH2XNowM6sxfD93yDqWhg/DzwbmZ1KiHrDpYuD7Dk0UFYrfPM4JL5urMVw3Qfg4WV2KiHqFac9CU6IKllKYOkDsGWerMUghB1JcRCuo6TQGHjetQQSnpApt4WwI5cuDjJ9RgNSlAuf3Qz7VxprMciU20LYlYw5COdXkHV2LYbRb0thEMIBXHrPQTQAuRnwyVhjIr1xs40jk4QQdufSew5yKGs9l30UPrwC0vcaazFIYRDCYVy6OEi3Uj124gDMHgHZx4y1GGSRHiEcSrqVhPNJ3QkfXwuWIpi0BFrGmZ1IiAbHpfccRD10pHQtBjDWYpDCIIQppDgI5/HHavhoNHgHGGsxhHY1O5EQDZYUB+Ec9n4H88YZazHc/i0ERZidSIgGzaWLgxytVE8kLYIFN0FIV5i8HAJamJ1IiAbPpYuDHK1UD2ycAwvvMNZimLRU1mIQwkm4dHEQLu7Xt2Dpg9DxMpi4UNZiEMKJyKGswvG0hp9egsTXIPIaGPsfmXJbCCcjxUE4ltUK302H9e9Bz1vg6n+Bm7vZqYQQFUhxEI5Tdi2Gfn+BES/KlNtCOCkpDsIxSopg8Z2w82tImA5DHpfCIIQTc+niIOs5uIiiPPj8Ftj3A4x4Cfr/xexEQojzcOmjleRQVhdQkAWfXAf7foTRb0lhEMJFuPSeg3ByuZnwybWQusNYiyF6rNmJhBA1JMVB2Ef2MZg7Bk4dghvnQ+fLzU4khKgFKQ7CPpY/AlkpcPMiaHeJ2WmEELXk0mMOwkkd2QS7/weXPCSFQQgXJcVB1L2fXoRGQdD3HrOTCCEukBQHUbcOrTUOWb3kIZkrSQgX5nTFQSnVTSn1nlJqoVLqXrPziFrQGla+AH6h0Psus9MIIS6CQ4qDUmq2UipNKZVUYftIpdQepdQ+pdQ0AK31Lq31PcANwEBH5BN15I+f4dAaGPwIePmanUYIcREctecwBxhZdoNSyh14B7gCiAQmKKUiS+8bDSwDljson7hYZ/YaAlpBr8lmpxFCXCSHFAetdSJwosLmPsA+rfUBrXURsAAYU/r4JVrrK4CJjsgn6kDyCkj53dhr8PA2O40Q4iKZeZ5DS+BwmdspQF+lVAIwFvDmHHsOSqkpwBSANm3a2C2kqIEzew1N20HPm81OI4SoA053EpzWehWwqgaPmwXMAoiPj9f2TSXOaddSOL4NrnkP3D3NTiOEqANmHq10BGhd5nar0m01ppS6Wik1Kysrq06DiVqwWoxV3YI7Q+wNZqcRQtQRM4vD70AnpVSEUsoLuBFYUpsXkFlZnUDSYkjfZazRICu6CVFvOOpQ1vnAWqCLUipFKXWH1roEmAp8B+wCPtda76jl68qeg5ksJbDqJWgebawFLYSoNxwy5qC1nlDN9uVcxOGqWuulwNL4+Hg548oMW+fDiQPGrKtuTnc+pRDiIrj0/2jZczBRSSH8/Cq0iIMuV5idRghRx1y6OMiYg4k2zYWswzD0KVkLWoh6yKWLgzBJcT4kzoA2A6DDULPTCCHswKWLg3QrmeT3DyDnOAx9UvYahKinXLo4SLeSCQpzYM0b0D5BFvIRoh5z6eIgTPDb+5CXAZc+ZXYSIYQduXRxkG4lB8s/Bb/8CzqPhNa9zU4jhLAjly4O0q3kYOvehYIsuPQJs5MIIezMpYuDcKDcTFj7LkSOgfDuZqcRQtiZFAdRM7/+C4pyIEH2GoRoCFy6OMiYg4OcToX1s4xZV0O7mp1GCOEALl0cZMzBQdb8EyxFMORxs5MIIRzEpYuDcICsFNgwG3rcBM06mJ1GCOEgUhzEuSXOMJYBHfKY2UmEEA4kxUFU78QfsPlj6DUZmsg63UI0JC5dHGRA2s5+fg3cPGDQ/5mdRAjhYDUuDkqpF6rYZuq6kDIgbUfpe2HbAuh9JwSEm51GCOFgtdlzaKmUsq3oppQKBX6o+0jCKax6GTwawSUPm51ECGGC2iwTejfwnVJqP6CBDwE5trE+Op4EOxbDoEfAL9jsNEIIE5y3OCil5gKbgM3AX4BPgRLgGq31PvvGE6b46SXwDoQBU81OIoQwSU26leYACrgN+ARoB5wEblZKjbNbMmGOIxthzzIYcD80amp2GiGESc6756C1XgmsPHNbKeUBdAO6A32BhXZLJxxv5YvQKAj63WN2EiGEiWoz5gCA1roE2F56+aTOE9WCUupq4OqOHTuaGaP+OLQW9v8Iw58Hb3+z0wghTOTS5znIoax1SGtY+QL4hULvu8xOI4QwmUsXB1GH/vgZDq2BwY+Al6/ZaYQQJpPiIM7uNQS0MqbKEEI0eFIcBCSvgJTfYcij4OFtdhohhBOQ4tDQWa3GXkPTdtBjotlphBBOotZHK4l6ZvdSOL4Nrn0f3D3NTiOEcBKy59CQWS3G2dDBnSHmerPTCCGciOw5NGRJiyB9N1w/B9xMnWBXCOFknLI4KKWuAa4CAoAPtNYrzE1UD1lKjJlXm8dAtzFmpxFCOBmHdSsppWYrpdKUUkkVto9USu1RSu1TSk0D0Fp/pbW+C7gHGO+ojA3K1vlw4gBc+gS4Se+iEKI8R7YKc4CRZTeULhb0DnAFEAlMUEpFlnnIU6X3i7pUUgg/vwot4qDLFWanEUI4IYcVB611InCiwuY+wD6t9QGtdRGwABijDK8C32itNzkqY4OxaS5kHYahT4FSZqcRQjghs/sTWgKHy9xOKd12P3AZME4pVeX0oEqpKUqpDUqpDenp6fZPWl8U50PiDGgzADoMNTuNEMJJOeWAtNZ6JjDzPI+ZBcwCiI+P147IVS/8/gHkHIdxs2WvQQhRLbP3HI4ArcvcblW6rUaUUlcrpWZlZWXVebB6qTAH1vwT2l8K7QaanUYI4cTMLg6/A52UUhFKKS/gRmBJTZ8sU3bX0vr3IC/TGGsQQohzcOShrPOBtUAXpVSKUuqO0oWDpgLfAbuAz7XWO2rxmrLnUFP5p+DXmdD5CmgVb3YaIYSTU1q7fnd9fHy83rBhg9kxnNtPLxmHr969GsJjzU4jhHACSqmNWusqvy2a3a10UWTPoYZyM2HtuxB5jRQGIUSNuHRxkDGHGvr1X1CUAwnTzU4ihHARTnkoq6hDp1Nh/SyIvQFCu5qdRghRBatVk1tUQm6hhZzCEnIKS8gt89O4brFtK3v/lMHtGdQppM4zuXRxUEpdDVzdsWNHs6M4rzX/BEsRDHnc7CRC1CuFJRZyCs426LlF5Rvz06X3VdxuXC9fBPKKLDX6nW4KGnt70NjbA7/SS1GJ1S7vz6WLg9Z6KbA0Pj7+LrOzOKWsFNgwG3pOhGYdzE4jhKlKLFajUS4621DbGvbCEvKKyn87zy1t8I1GvnwRyC0sodhSs4N5vD3c8Pcpbcy9jIY9uLEXbZv5Gtu9jPvONvju5QpA2Z8+nm4oB5286tLFQZxH4uvGz8GPmZtDiAtQbLHaGuq8Iku5b9+Vv5FbbI15TqGFvDMNfJmumpp+w3ZTGI24z9lv54293Qlu7G002j5lGm0v9/INuE+Z615GQ+/h7ppDuy5dHKRb6RxOHIDNn0D87dCk9fkfL0Qd0FqTW2QhK7+YrLxiTuUXkZ1fvo88r6j8N/bcIku5LpczhaBWjbn3mW/g7rbGuZWvL429yzfevl7u5b6J+1Zo3P283Wnk6e6wb+cXIy0vje3p24kKjiLML6zOX9+li4N0K53Dz6+BmwcM+j+zkwgXo7Umv9ho4E/lFRsNfeklu/RnVdtPlf4ssZ67u+VMY16239zPy50gP9/Sbe62rpYz39p9vSp3u5zZ5siuFrMUlBSw68QutqVvMy4Z2zieexyAp/o+xfiudb/sjUsXB1GN9L2w7TPo/xfwr/tvFMI1FJQ28FU15ll5ReVvl2ncs/KLz9mfrhQE+HjSxNeTwEbGpWXTRgQ28qRJo7PbzlwCGnna+twbe3vg7VH/G/OLobXmUPYhtmVssxWD5JPJlOgSAFo2bknPkJ7ERMYQGxJL1yD7HIUoxaE+WvUyePrCwIfMTiLqQEGxhRO5ReUuJ0sb91N5Zxv0so18Vn7xebtlAnw8CCzTwIcF+hDYyKtS4162CAQ08sTf2wM3N2nc60pWYRbbM7bb9gi2p28nuygbAD9PP6KDo7kt+jZigmOICYkhuFGwQ3K5dHGQMYcqHE+CHYth0CPg55g/IlFzVqsmu6C4UmOfmVvEyTO388rfd67DHP29PQgo05B3DG1criGv2Lifufj7eOIuDbzDFVuL2XtyL9vSjSKwLWMbh7IPAeCm3OjQpAPD2w4nNiSW2OBYIgIjcHdzNyWrzK1U38y/CQ6ugYe2QqOmZqep9wpLLJzMLdPY5xVxIqewykbe+MZfjKWaPvlGnu4E+XnV6NLU14sAHw+XPRKmIdBaczz3uK17aHvGdnZm7qTQUghAcKNgYoNjiQmJITY4lqjgKPw8/Rya8VxzK7n0noOo4MhG2LMMLn1KCsMF0FpzurCEEzlnGvmqG/iyl5zCkipfSylo0sjT1phHBPvRq23T0tveBPl5Gj99vQhq7EWQrxeNvMz5hijqRl5xHjsyd9jGCbZnbCc931il0tvdm25B3RjfZTwxITF0D+5OmF+YU4+9SHGoT1a+CI2CoF+VK6s2WBarJjO3kLTsQtJOF5CabVxPPV1AWnYBaacLSc0u4ERuUbUDsV4ebjQr8829bTNfmvp60czPi6Z+XuXuC/Lzoomvl3Tb1GNWbeXAqQPGWEHpnsG+U/uwamOcp21AW/qG97V1D3Vu2hlPd0+TU9eOFIf64tCvsP9HGP538PY3O41DWK2aE3lFpGYXGI19mYY+NbuQ9NJCkJ5TWGVXTpCfF6H+3jQP8KFLc3+aNfautrH39XKNY9+FfWTmZ5YbNN6RsYOc4hwA/L38iQ2OZVibYcagcXAMTXyamBu4Drh0cZAB6VJaw8oXoHFz6H2n2WkumtWqOZlXRGrZb/el11OzC0k7XUhadgHppwurPKa+qa8nzQN8CA3woVNzf5oHGAUg1N+H0NLrIY298fKQ/npRWZGliN0ndtsKwbb0bRzJMVYvdlfudG7amavaX0VsSCwxwTG0DWiLm6p/f0suXRzkJLhSB1bBoV/gitfBy9fsNNXSWnMyr7j0m32BrZFPLfOt/0w3T1WNfhNfT5qXNvAdQ4JtjX7zAG9C/M/89MbbQ/ruzZaRn8G6Y+tYd3QdG1M3UmApMDtSjWUVZlFsLQYgzC+MmOAYbuxyI7EhsXRr1o1GHo1MTugYLl0cBGf3GgJaQa9JpsUoLLFwPKuAY1kFHMvK51hWAcezCsr07Rv9/VX16Qc28rQ19O1DmhkNfml3T2iAN6H+PoT4e+PjKY2+s8orzmND6gajIBxbR/LJZAACvALoE9bHpbpZArwCbEcRhfqGmh3HNFIcXF3yCjiyAa6eCR7edvkVhSUWUrMKbY1+2QJwLCufY6cKyMwtqvS8AB+P0m/2PvSN8CO09Fu+0cVT2r0jjb5LKrGWsCNzB+uOrmPtsbVsTd9KibUELzcvejbvyYNxD9K/RX+6Nu1q2nH64uJIcXBlVqux19A0AnrcdEEvUVRiJTW7QoN/qmwRKCAjp7DS8wJ8PAgPbER4Ex9iWjYhPNCn9GJsCwvwwc9b/rzqC601B7MP2rqKfj/+O6eLTwPQLagbt0TeQr/wfsSFxuHj4WNyWlEX5H+vK9u9FI5vg2vfhyoOkyu2nG34j57Kr9Ttc6bhr3gepL+Ph62hj2oRYGvwbY1/oDT8DUFGfgbrj623dRWdmeithV8LLm93Of1a9KNvWF+a+sg5NfWR/A93UcXFxfDDC5QEduAHPZBjifs5eqqgtAAYjX96VQ2/twdhgT6EN2lEt7AAwpv40CKwEWGBPrRo4kNYYCMaS8PfIOUV57EpbZOtq2jvyb2A0QffN7wvd8XcRf/w/rTybyWH9TYALt0K1PdDWbXWpJ0uZOfRbHYey2bXsWwOn8zneFY+A3J/5A3PvTxY9ADLF2wDwM/LnfAmxjf7LmH+tm/54U0a0SLQh7BAH/x9XOtEHGE/FqvFGDco3TPYkraFYmsxnm6exIXG8WDcg/QL70e3oG4ybtAAydxKTqLEYuVARi47jxpFYOexbHYezS430Ns6qBHtmvnRwt+DaQcmgacvW69aQngTP8Kb+BAgDb84B601f57+k7VH17Lu2Dp+O/4bp4uMcYOuQV3pF96P/uH96dm8Z4M5XLOhk7mVnMzpgmJ2Hz9tFIHSvYLdx0/bplj2cnejc1hjhnULJTI8gMgWgXQN9z/b+G+aCzsPw7ULSOgi6zWI6p0oOGEbN1h7dC3Hco8BEO4XzvC2w+kX3o++4X0J8gkyOalwNlIc7EhrzbGsgnJFYOexbA5l5tke09TXk8gWAUzq35bIFgFEhgfSPsQPz+pm2ywpNFZ5a9kLOo900DsRriK/JJ9NqZtsXUW7T+wGjCke+ob15Y7oO+jXoh9t/NvIuIE4JykOdaTYYmVfWk6lQnAqr9j2mHbNfIlqEcD1vVoR2SKAbuEBhAX41O4/6aa5kHUYRs80pv4UDZrFamHXiV22rqLNaZspthbj4eZBz9Ce3N/zfvqH9yeyWaSMG4hakeJwAbLyi9ldZlxg57FsklNzKLIY3ULeHm50DfPniugwIsONItA1PODijwIqzofEGdB2ILS/tA7eiSi2FnM85zgpOSm2KRNcwfHc46w9upb1x9fbxg26NO3CTV1vol8L43wDX0/nnUpFOD8pDuegtSblZH65AeKdx7JJOZlve0wzPy8iWwRw2yXtjPGB8AAigv3sswjL7x9AznEYN1v2GmohrziPw6cPk3I6hcOnD5e7HMs9hkVXv9KaM2vu25xhbYbRP7w/fcL7OGz5SNEwSHEoVVRiJTnt9NkuodKjhrILjMVclIKIYD+6t27ChD5tiGwRQFR4ACH+3o7pu83cD6v/AR2GQruB9v99LkRrzYmCE7YGv2IRyCzILPf4QO9AWjduTUxwDFdEXEFr/9a08m/lUkfoBHgF0Nq/tYwbCLtp0MXh573pfL3lCDuPZrMvLcc2E2gjT3e6hPkzqnuL0qOFAuga5o+vl0kf16k/Ye4Yo0KNfNWcDCYrsZZwPPd4tQUgr+TsIL9C0dyvOa39WzOk9RBb49/avzWt/VsT4BVg4jsRwjU4XXFQSrUHngQCtdbj7Pm7klNPszo5g8jwAC7tGmorBO2a+TnPKl7Zx+Cj0VCYDZOWQkhnsxPZTX5JfqVG/8ztozlHKdFnl+T0dPO0NfjxYfG2hr+VfytaNm6Jt7t9JiEUoqFwyElwSqnZwCggTWsdXWb7SOBfgDvwX631K2XuW1jT4nChJ8FZrNp5ikBVcjPgwysh+wjc8hW07m12oouitSarMKtSv/+ZIpCWn1bu8f5e/rZGv+Il1De0Xi6wUheKi4tJSUmhoMB11lAQ9uXj40OrVq3w9Cx/oqwznAQ3B3gbmFsmlDvwDjAcSAF+V0ot0VrvdFAm5y4M+Sdh7jVGl9LNC12qMBRZithzYg97T+6tVADOzOR5RmijUFr5t2JAywGVCkCgd6BJ78C1paSk4O/vT7t27WRMQqC1JjMzk5SUFCIiImr8PIcUB611olKqXYXNfYB9WusDAEqpBcAYoEbFQSk1BZgC0KZNm7oL6wwKsuGT6yBjD0yYD+0uMTtRtSxWC/uz9rMjYwdJGUkkZSax9+ReSqxGF5CHmwctG7eklX8ruod0L9f4t/JvJdM720FBQYEUBmGjlKJZs2akp6fX6nlmjjm0BA6XuZ0C9FVKNQNeBHoqpaZrrV+u6sla61nALDC6lewd1mGK8uDT8XB0C4z/GDpeZnYiG601KadTSMpMMgpBRhK7Tuwiv8Q4tLexZ2OigqOYFDmJ6OBougZ1JdwvXE6+MoEUBlHWhfw9ON2AtNY6E7inJo+td7OyFhfAgpvg8Dq47r/Q9SpT46Tnpdv2BnZk7CApM4mswiwAvN296RrUlbGdxhLVLIro4Oh6u9C6EA2RmcXhCNC6zO1WpdtqTGu9FFgaHx9/V10GM4WlGL6YDAd+gjHvQvR1Dv312UXZ7MjYwY5Mo3toe8Z20vKMAWJ35U7HJh25rM1lRAVHEd0smo5NO+LpJrPAiuqlpqby8MMPs27dOpo2bYqXlxePPfYY1157rcMyrFu3jgcffJDCwkIKCwsZP348zz77bLWP37BhA3PnzmXmzJmV7mvXrh0bNmwgOLj6kw1r8hhXYWZx+B3opJSKwCgKNwK1Wuuy3uw5WC2w+C7Y+w1cOQN6TrTrr8svyWfPiT1sz9hOUkYSOzJ3cCj7kO3+tgFtiW8eT3RwtK17yJVOEBPm01pzzTXXMGnSJD799FMADh06xJIlSxyaY9KkSXz++ed0794di8XCnj17zvn4+Ph44uOrPHinwXFIcVBKzQcSgGClVArwjNb6A6XUVOA7jENZZ2utd9TmdevFnoPVCl9PhR1fwvC/Q5+6fSvF1mL2n9rP9ozttkHjfaf22aaMCPUNJbpZNGM6jCEqOIqoZlFylJC4aCtXrsTLy4t77jnbQ9y2bVvuv/9+Dh48yC233EJubi4Ab7/9NgMGDGDVqlXMmDGD//3vfwBMnTqV+Ph4Jk+ezLRp01iyZAkeHh5cfvnlzJgxgy+++ILnnnsOd3d3AgMDSUxMrJQjLS2N8PBwANzd3YmMjAQgNzeX+++/n6SkJIqLi3n22WcZM2ZMuQyZmZlMmDCBI0eO0L9/f8oe9v/JJ58wc+ZMioqK6Nu3L++++y7u7vVrbM1RRytNqGb7cmC5IzI4Ja1h+SOw9VNIeAIGPnBRL2fVVg5lH7LtDSRlJLH7xG4KLYWAMeVCdHA0g1sNtu0VhPqG1sU7EU7quaU72Hk0u05fM7JFAM9cHXXOx+zYsYO4uLgq7wsNDeX777/Hx8eH5ORkJkyYwLnOU8rMzOTLL79k9+7dKKU4deoUAM8//zzfffcdLVu2tG2r6OGHH6ZLly4kJCQwcuRIJk2ahI+PDy+++CJDhw5l9uzZnDp1ij59+nDZZeUP/njuuee45JJLePrpp1m2bBkffPABALt27eKzzz7jl19+wdPTk/vuu4958+Zx6623nvMzcTVONyBdGy7draQ1rHgKNnwAAx+EIY/V8uma1LxU21FDSZlJ7MzYaTuPoJFHI7oFdeOGLjcQ3SyamOAYWftXmOYvf/kLa9aswcvLix9++IGpU6eyZcsW3N3d2bt37zmfGxgYiI+PD3fccQejRo1i1KhRAAwcOJDJkydzww03MHbs2Cqf+/TTTzNx4kRWrFjBp59+yvz581m1ahUrVqxgyZIlzJgxAzAO//3zzz/LPTcxMZHFixcDcNVVV9G0aVMAfvzxRzZu3Ejv3sa5R/n5+YSG1r8vWS5dHFy6W2nVy7D2begzBS57rkazrG5J28K6Y+tsBeHMhHIeyoNOTTsxMmKkbY+gfWB7PNxc+p9X1IHzfcO3l6ioKBYtWmS7/c4775CRkUF8fDxvvPEGzZs3Z+vWrVitVnx8jHNdPDw8sFqttuecOcPbw8OD3377jR9//JGFCxfy9ttvs3LlSt577z3Wr1/PsmXL6NWrFxs3buSRRx5h8+bNtGjRguXLjU6JDh06cO+993LXXXcREhJCZmYmWmsWLVpEly5dyuVOTU0973vTWjNp0iRefrnKo+zrDZc+7lApdbVSalZWVpbZUWpnzRvw86vQ82ZjIr3zFIaCkgL+vvbv3PLNLby75V3+PP0nA1oMYFqfacy7ch7rJq7j86s/5+n+TzO201g6N+0shUGYaujQoRQUFPDvf//bti0vz5gcMSsri/DwcNzc3Pj444+xWIzxr7Zt27Jz504KCws5deoUP/74IwA5OTlkZWVx5ZVX8sYbb7B161YA9u/fT9++fXn++ecJCQnh8OHDfPjhh2zZssVWGJYtW2YbK0hOTsbd3Z0mTZowYsQI3nrrLdt9mzdvrvQeBg8ebBtM/+abbzh58iQAw4YNY+HChaSlGUfznThxgkOHDlV6vqtz6RbEJfcc1r8PPzwL0ePg6pngdu76vO/kPh5NfJR9p/YxOWoyU2Kn4O/l75isQlwgpRRfffUVDz/8MK+99hohISH4+fnx6quvEhcXx3XXXcfcuXMZOXIkfn5+ALRu3ZobbriB6OhoIiIi6NmzJwCnT59mzJgxFBQUoLXmn//8JwCPPvooycnJaK0ZNmwY3bt3r5Tj448/5uGHH8bX1xcPDw/mzZuHu7s7f/vb33jooYeIjY3FarUSERFhGwg/45lnnmHChAlERUUxYMAA20wMkZGRvPDCC1x++eVYrVY8PT155513aNu2rT0/UodzyMR79nahE+853Ka5sOR+6DoKrp8D7tWfJ6C15ou9X/Da76/h5+nHS5e8xMCWso6DOL9du3bRrVs3s2MIJ1PV34UzTLwntn0BSx4wpsMYN/uchSGrMIvn1j7H94e+p394f14a9JKs8iWEcCiXLg4uc7TSrqXw5d3GBHrjPwGP6tca2JS6icdXP05GXgZ/7fVXJkVNkikphBAO59KtjtZ6qdZ6SmCgE5+0lfw9fHEbtOxlzLDqWfWZxharhfe2vsdt392Gh/Lg4ys/5rbo26QwCCFM4dJ7Dk7vj0T47GZoHgkTvwDvqgeSj+ceZ/rq6WxI3cCVEVfyt35/o7FXYweHFUKIs6Q42Muf6+HTG6FpBNz8JTRqUuXDfvrzJ/72698oshTxwsAXGN1htJyoJoQwnUsXB6cdczi6GeaNA/8wuPUr8GtW6SGFlkL+seEfzN89n25B3Xht8Gu0C2zn8KhCCFEVl+7Qdsoxh9Sd8PG14NMEJi0xCkQFB04d4KZlNzF/93xuibyFT678RAqDqFdefPFFoqKiiI2NpUePHqxfvx6AN99803YynFlOnTrFu+++W+397u7u9OjRg6ioKLp3784//vGPcmdum+Gll15y/C/VWrv8pVevXtoppCdr/VpHrWd00Tpzf6W7rVarXrR3ke79SW89aP4g/fPhn00IKeq7nTt3mvr7f/31V92vXz9dUFCgtdY6PT1dHzlyRGutddu2bXV6erqZ8fQff/yho6Kiqr3fz8/Pdj01NVUPGzZMP/30046IVq2ymS5UVX8XwAZdTbvq0nsOTuXkIZg7GrQVbl0CQe3L3Z1dlM1jiY/xzK/PEBscy8LRCxncarBJYYWwn2PHjhEcHIy3t3HIdnBwMC1atGDmzJkcPXqUSy+9lEsvvRSAFStW0L9/f+Li4rj++uvJyckBYOPGjQwZMoRevXoxYsQIjh07BkBCQgIPPvggPXr0IDo6mt9++w0wpuC+/fbb6dOnDz179uTrr78GjNlh+/TpQ48ePYiNjSU5OZlp06axf/9+evTowaOPPnrO9xIaGsqsWbN4++230VpjsVh49NFH6d27N7Gxsbz//vu29zx48GBbrtWrVwPw7bffEhcXR/fu3Rk2bNg5s86ZM4exY8cycuRIOnXqxGOPGZNxTps2jfz8fHr06MHEifZd66UsOUO6LmQfhdkjoSALJv8PwmLK3b0lbQvTVk/jeO5xpvacym1Rt8m6ysJuyp0J+800OL69bn9BWAxc8Uq1d+fk5HDJJZeQl5fHZZddxvjx4xkyZAhQfqW0jIwMxo4dyzfffGObWqOwsJDp06czZMgQvv76a0JCQvjss8/47rvvmD17NgkJCXTq1In//Oc/JCYmct9995GUlMQTTzxBZGQkN998s20K7s2bNzNt2jT69evHxIkTKSoqwmKxkJqayqhRo0hKSqoyf+PGjW1F6owmTZqwZ88evv76a9LS0njqqacoLCxk4MCBfPHFFyxevJiCggKefPJJLBYLeXl5FBQUEBcXR2JiIhEREZw4cYKgoKBqs37xxRc8//zzbN68GW9vb7p06cKaNWto3bp1lZlqq0GdIe0UA9I5afDRaMg7AZO+LlcYLFYLH+74kLc3v02YXxgfXfER3UMqz/8iRH3SuHFjNm7cyOrVq/npp58YP348r7zyCpMnTy73uHXr1rFz504GDjSmhSkqKqJ///7s2bOHpKQkhg8fDoDFYrEt2AMwYYKxPMzgwYPJzs7m1KlT1U7B3b9/f1588UVSUlIYO3YsnTp1uqj3tmLFCrZt28bChQsBYxLB5ORkevfuze23305xcTHXXHMNPXr0YNWqVQwePJiIiAgAgoKCbK9R3XThw4YN48wYamRkJIcOHaJ169YVYziESxcHbfbEe3knYO41kH0Ebl5snOhWKi0vjSdWP8H64+sZ2W4kT/d/WibME453jm/49uTu7k5CQgIJCQnExMTw0UcfVSoOWmuGDx/O/Pnzy23fvn07UVFRrF27tsrXrniot1Kq2im4u3XrRt++fVm2bBlXXnkl77//Pu3bl+/yPZ8DBw7g7u5OaGgoWmveeustRowYUelxiYmJLFu2jMmTJ/PXv/7Vtv5DRdVlXb9+va0rDozPsKSkpFZZ65KMOVyogiz4ZCxk7oMbP4W2/W13JaYkMm7JOLZlbOP5Ac/z2uDXpDCIBmPPnj0kJyfbbm/ZssU2Y6m/vz+nTxsLUvXr149ffvmFffv2AUZf/N69e+nSpQvp6em24lBcXMyOHWdXEP7ss88AWLNmDYGBgQQGBlY7BfeBAwdo3749DzzwAGPGjGHbtm3lMpxPeno699xzD1OnTkUpxYgRI/j3v/9NcXExAHv37iU3N5dDhw7RvHlz7rrrLu688042bdpEv379SExM5I8//gCMqb2BGk0XXpGnp6ftdzqKS+85mKYoF+bdYPTljp8HHYzBtSJLEW9sfINPdn1C56adeX3w67RvUrtvKUK4upycHO6//35OnTqFh4cHHTt2ZNasWQBMmTKFkSNH0qJFC3766SfmzJnDhAkTKCw0lrJ94YUX6Ny5MwsXLuSBBx4gKyuLkpISHnroIaKijIWLfHx86NmzJ8XFxcyePRug2im4P//8cz7++GM8PT0JCwvjiSeeICgoiIEDBxIdHc0VV1zB66+/Xi7/mcHf4uJiPDw8uOWWW/jrX/8KwJ133snBgweJi4tDa01ISAhfffUVq1at4vXXX8fT05PGjRszd+5cQkJCmDVrFmPHjsVqtdqWR63JdOEVTZkyhdjYWOLi4pg3b16d/ntVRwaka6u4AD69AQ6uNmZXjboWgINZB3ks8TF2ndjFTV1v4q/xf8XbvfoJ9oSwl/o8ZXdCQgIzZswgPr7KMVRxDg1qQNrhSorg81uNOZOufQ+irkVrzdf7v+al9S/h7e7NW0PfIqF1gtlJhRDiokhxqClLCSy6A5K/g1FvQPcbySnK4e/r/s7yP5bTO6w3L1/yMs39mpudVIh6a9WqVWZHaDBcujg47FBWqxW+vg92LYERL0P87WxP385jiY9xLPcYU3tM5c6YO+XcBSFEveHSRytpR8ytpDUsexi2fQZDn8La7x5mJ83m1m9uxaItfDjyQ+7ufrcUBiFEveLSew52pzV8Ox02zoFB/0dG79t44vt7WHtsLcPbDueZ/s8Q6O1Ek/4JIUQdkeJwLiv/Duv/DX3v5ZcuQ3liyXXkFufydP+nGddpnKy7IISot1y6W8muEl+H1f+gOO5W/hEczD0/3kuQTxALrlrA9Z2vl8IgxDkcP36cG2+8kQ4dOtCrVy+uvPJK9u7de0Gvdeedd7Jz5846zXe+aburk5OTw7333kuHDh2Ii4ujV69e/Oc//7moLHPmzGHq1KkAvPfee8ydO/eCXufgwYN8+umnF5WlLCkOVVn7Lqx8gT+jR3OLWzpzdn7E+C7jmX/VfDo2dbKFhYRwMlprrr32WhISEti/fz8bN27k5ZdfJjU19YJe77///S+RkZF1mvFCi8Odd95J06ZNSU5OZtOmTXz77be2M5/LutBpL+655x5uvfXWC3quFAd72zAbvpvO0s6XcH3hXg6fPsybCW/yVL+n8PHwMTudEE7vp59+wtPTk3vuuce2rXv37gwaNAitNY8++ijR0dHExMTYpsJYtWoVCQkJjBs3jq5duzJx4kTb9BIJCQmcOcm1ceOza6svXLjQNl/T/v376devHzExMTz11FO2x+Xk5DBs2DDi4uKIiYmxTY9d1bTdr7/+um0q7meeeabS+9q/fz+//fYbL7zwAm5uRtMZEhLC448/bnsPgwYNYvTo0bZids0119CrVy+ioqJsZ4kDfPjhh3Tu3Jk+ffrwyy+/2LY/++yztgn59u/fz8iRI+nVqxeDBg1i9+7dAEyePJkHHniAAQMG0L59e9skgNOmTWP16tX06NGDN954o5b/apXJmENZWxeQu+z/eKl9DEuK/yQuNI5XBr1CeOPw8z9XCCf06m+vsvvE7jp9za5BXXm8z+PV3p+UlESvXr2qvG/x4sVs2bKFrVu3kpGRQe/evRk82FjXZPPmzezYsYMWLVowcOBAfvnlFy655JIaZXrwwQd58MEHmTBhAu+9955tu4+PD19++SUBAQFkZGTQr18/Ro8ezSuvvEJSUhJbtmwBjJlSk5OT+e2339BaM3r0aBITE23ZwFgbonv37rbCUJVNmzaRlJRkm4l19uzZBAUFkZ+fT+/evbnuuusoKirimWeeYePGjQQGBnLppZfSs2fPSq81ZcoU3nvvPTp16sT69eu57777WLlyJWCsH7FmzRp2797N6NGjGTduHK+88gozZsw471QcNSXF4YwdX7Jj+QM81q4dKZzmvu73cVfsXXi4yUckRF1Zs2YNEyZMwN3dnebNmzNkyBB+//13AgIC6NOnD61atQKgR48eHDx4sMbFYe3atXz11VcA3HTTTTzyyCOA0cX1xBNPkJiYiJubG0eOHKmye2vFihWsWLHC1kjn5OSQnJxcrjhU9OKLL/LFF1+QlpbG0aNHAejTp4+tMADMnDmTL7/8EoDDhw+TnJzM8ePHSUhIICQkBIDx48dXGo/Jycnh119/5frrr7dtOzP/FBh7JG5ubkRGRl5wd935OF3Lp5TyA94FioBVWmu7zzJl3b2cj1c8yJvhzWnWKIgPBr9KfJjM3SJc37m+4dtLVFSUraujNmoyXXXZA0EKCgrO+5rz5s0jPT2djRs34unpSbt27ap8ntaa6dOnc/fdd1f7WpGRkWzduhWr1YqbmxtPPvkkTz75ZLmuLj8/P9v1VatW8cMPP7B27Vp8fX1JSEioUWYAq9VKkyZNbHs2FZX9rOw1P55DxhyUUrOVUmlKqaQK20cqpfYopfYppaaVbh4LLNRa3wWMtne2zF1L+Muqh5gRFMjgVoNYNGaxFAYhLsLQoUMpLCws18e+bds2Vq9ezaBBg/jss8+wWCykp6eTmJhInz59avzazZs3Z9euXVitVts3cjCm/160aBEACxYssG3PysoiNDQUT09PfvrpJw4dOgRQadruESNGMHv2bNtqa0eOHCEtLa3c7+7YsSPx8fE89dRTWCwWwChQ1TXOWVlZNG3aFF9fX3bv3s26desA6Nu3Lz///DOZmZkUFxfzxRdfVHpuQEAAERERtvu01mzduvWcn01tpiKvCUcNSM8BRpbdoJRyB94BrgAigQlKqUigFXC49GEWe4Zat2U249ZO5zcfb57s+RBvDntHTmoT4iIppfjyyy/54Ycf6NChA1FRUUyfPp2wsDCuvfZaYmNj6d69O0OHDuW1114jLCysRq8J8MorrzBq1CgGDBhQbnW4N998k3/+85/Exsayb98+22pqEydOZMOGDcTExDB37ly6du0KQLNmzWzTdj/66KNcfvnl3HTTTfTv35+YmBjGjRtXZUP73//+l8zMTFuhGD58OK+99lqVmUeOHElJSQndunWzLVcKEB4ezrPPPkv//v0ZOHBgtTPozps3jw8++IDu3bsTFRVlG0yvTmxsLO7u7nTv3r1OBqQdNmW3Uqod8D+tdXTp7f7As1rrEaW3p5c+NAU4qbX+n1Jqgdb6xmpebwowBaBNmza9znwjqI1FKx7m45QfeW3ELDq37Ffr5wvhjOrblN0xMTEsWbKkXF9+RXl5eTRq1AilFAsWLGD+/PnnbUwbGleasrslZ/cQwCgKfYGZwNtKqauApdU9WWs9C5gFxnoOFxJg7PB/MqowG28f2VsQwhkNHz6cmJiYcxYGgI0bNzJ16lS01jRp0sS2CJC4cE43IK21zgVuq8ljL3ZWVqWUFAYhnNj3339fo8cNGjTovH3yonbMPAnuCNC6zO1WpdtqzCGzsgrhgurDCo+i7lzI34OZxeF3oJNSKkIp5QXcCCypzQsopa5WSs3KysqyS0AhXJGPjw+ZmZlSIARgFIbMzEx8fGo3w4NDBqSVUvOBBCAYSAWe0Vp/oJS6EngTcAdma61fvJDXd+ga0kI4ueLiYlJSUmp8TL2o/3x8fGjVqhWenp7ltp9rQNphRyvZQ5kxh7uSk5PNjiOEEC7lXMXBpSfekzEHIYSwD5cuDkIIIezDpYuDDEgLIYR9uPSYwxlKqXTgEBAI1LZSBAMZZZ5b9jVqcp1abqvN/bXJe65857vt6KwVX08+W3M+2+p+r3y2Deezbau1DqnyVbXW9eYCzLqA52wo+9yyr1GT67XdZq+858p3vtuOziqfrXN8ttX9XvlsG95nW9XFpbuVqlDtdBu1eO7SKrad63ptt9Xm/po893y5a3K7Jlnks636dV35s63u98pnW/65DeGzraRedCtdDKXUBl3NoVzOyJXyulJWcK28rpQVXCuvK2UF++Wtb3sOF2LW+R/iVFwprytlBdfK60pZwbXyulJWsFPeBr/nIIQQojLZcxBCCFGJFAchhBCVSHEQQghRiRSHc1BKtVdKfaCUWmh2lppQSl2jlPqPUuozpdTlZuc5F6VUN6XUe0qphUqpe83Ocz5KKT+l1Aal1Cizs5yPUipBKbW69PNNMDvPuSil3JRSLyql3lJKTTI7z/kopQaVfq7/VUr9anaec1FKtVFKfaWUmq2Umlbb59fb4lD6gaQppZIqbB+plNqjlNp3vg9Ma31Aa32HfZPactVF3q+01ncB9wDjnTzrLq31PcANwEBnzlrqceBz+6Qsl6su8mogB/DBWH7XmbOOwVjoq9ieWUtz1cXf7erSv9v/AR85c1YgBliotb4d6FnrELU9G89VLsBgIA5IKrPNHdgPtAe8gK1AZOmH+L8Kl9Ayz1voYnn/AcQ5e1ZgNPANcJMzZwWGYyxGNRkY5ex/B4Bb6fOaA/OcPOs04G5H/D+r4/9jnwP+zpwVaAb8BKwEbqt1Bnv+Y5h9AdpV+HD7A9+VuT0dmF6D17F7caiLvIACXgUuc/asFV5rmTNnBV7EWJRqBfD1mcbXWfOWeZyXAxrci/1sbwZuKL3+mT2z1tVnC7QB/uPsWYFHgMGl12v9d+BBw9ISOFzmdgrQt7oHK6WaYTQMPZVS07XWL9s5X0W1ygvcD1wGBCqlOmqt37NnuApq+9kmAGMBb2C5PYNVoVZZtdZPAiilJgMZWmurXdNVVtvPdiwwAmgCvG3XZJXV9m92MfCWUmoQkGjPYNWobV6AO4AP7ZaoerXN+i3wrFLqJuBgbX9ZQysOtaK1zsTov3cJWuuZwEyzc9SE1noVsMrkGLWitZ5jdoaa0Fovxmh0nZ7WOg+jsXUZWutnzM5QE1rrJGDchT6/3g5IV+MI0LrM7Val25yVK+WVrPbjSnldKSu4Vl6HZm1oxeF3oJNSKkIp5YUxyLjE5Ezn4kp5Jav9uFJeV8oKrpXXsVntPahi1gWYDxzj7CFyd5RuvxLYizHq/6TZOV0xr2SVvK6W1dXyOkNWmXhPCCFEJQ2tW0kIIUQNSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIeykdGGWLaWX9Uop+f8mXIacIS2EnSilkjHm0z9mdhYhaku+yQhhP8uBbUqpN80OIkRtyXoOQtiBUmoAxsp84VrrErPzCFFbsucghH1cD+zVWpcoQ4DZgYSoDRlzEMIOlFJ9gA8ADeQD92mtN5qbSoiak+IghBCiEulWEkIIUYkUByGEEJVIcRBCCFGJFAchhBCVSHEQQghRiRQHIYQQlUhxEEIIUYkUByGEEJX8PwHirA85xCEuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "A = np.array([[  .2,   .1,   1.0,  1.0,    0.0],\n",
    "              [  .1,  4.0, - 1.0,  1.0, -  1.0],\n",
    "              [ 1.0, -1.0,  60.0,   .0,     .0],\n",
    "              [ 1.0,  1.0,    .0,  8.0,    4.0],\n",
    "              [  .0, -1.0, - 2.0,  4.0,  700.0]])\n",
    "b = np.array( [ 1  ,  2  ,   3  ,  4  ,    5  ] )\n",
    "x_exact = linalg.solve(A, b)\n",
    "\n",
    "eps_list  = np.logspace(-8, -1, 8)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\")]:\n",
    "    k_list = []\n",
    "    for eps in eps_list:\n",
    "        x, k = alg(A, b, eps)\n",
    "        k_list.append(k)\n",
    "    ax.plot(eps_list, k_list, label=alg_name)\n",
    "    \n",
    "ax.set_xscale(\"log\")\n",
    "ax.invert_xaxis()\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_xlabel(\"$\\epsilon$\")\n",
    "ax.set_ylabel(\"$k$\")\n",
    "ax.legend()\n",
    "\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The condition number for A is K(A) = 12269.877667702964 >> 1,\n",
      "so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\n",
      "This explains the great difference in order of required iterations k as observed above.\n"
     ]
    }
   ],
   "source": [
    "print(\"The condition number for A is K(A) =\", linalg.cond(A), \">> 1,\")\n",
    "print(\"so A is ill-conditioned, so the Conjugate Gradient method is highly susceptible to rounding errors.\")\n",
    "print(\"This explains the great difference in order of required iterations k as observed above.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "deletable": false,
    "editable": false,
    "nbgrader": {
     "cell_type": "markdown",
     "checksum": "4d7421993a0fef53e4f0f83ffae975b1",
     "grade": false,
     "grade_id": "cell-c1b8b533da043a26",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "source": [
    "### Task 8\n",
    "\n",
    "Try to get a better convergence behavior by pre-conditioning your matrix $A$. Instead of $A$ use\n",
    "\n",
    "$$ \\tilde{A} = C A C,$$\n",
    "\n",
    "where $C = \\sqrt{D^{-1}}$. If you do so, you will need to replace $\\mathbf{b}$ by \n",
    "\n",
    "$$\\mathbf{\\tilde{b}} = C \\mathbf{b}$$\n",
    "\n",
    "and the vector $\\mathbf{\\tilde{x}}$ returned by your function will have to be transformed back via\n",
    "\n",
    "$$\\mathbf{x} = C \\mathbf{\\tilde{x}}.$$ \n",
    "\n",
    "What is the effect of $C$ on the condition number and why?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "deletable": false,
    "nbgrader": {
     "cell_type": "code",
     "checksum": "213ac0f9b57d67bb90bbd75dd3808955",
     "grade": true,
     "grade_id": "cell-de3d1bd5ccfa8dfb",
     "locked": false,
     "points": 3,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The number of iterations is brought down by a lot due to the conditioning, bringing down the runtime.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAF3CAYAAAB65myjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACeRElEQVR4nOzdd3zV9b3H8dc3J3svQhICIYS9SQKIEzcO1GodaK17tLXz1ta22t3beXtva22tA6nWWbVVKyiodYOQAKKAKBvCCGTv5JzzvX/8Tk4S9kjyy0neTx95nHN+53dOPgch37x/32WstYiIiIiIiIiEojC3CxARERERERE5Vgq1IiIiIiIiErIUakVERERERCRkKdSKiIiIiIhIyFKoFRERERERkZClUCsiIiIiIiIhK9ztArpCenq6HTp0qNtliIhIH1FSUrLXWjvA7TpCkTFmNjA7ISHhlpEjR7pdjoiI9BGHaptNX9intqioyBYXF7tdhoiI9BHGmBJrbZHbdYQytc0iItKVDtU297rhx8aYmcaYd4wx9xtjZrpdj4iIiIiIiPRePRJqjTFzjTFlxpiP9zk+yxizzhiz3hhzV+CwBeqAaGB7T9QnIiIiIiIioamnemrnAbM6HjDGeID7gPOAscAcY8xY4B1r7XnAd4Gf9FB9IiIiIiIiEoJ6JNRaa98GKvY5PA1Yb63daK1tAZ4CLrbW+gPPVwJRPVGfiIiIiIiIhCY3Vz8eBGzr8Hg7MN0YcylwLpAM/OlgLzbG3ArcCjBkyJDuq1JERERERER6rV63pY+19nng+SM47wHgAXBWWOzuukRERERERKT3cXP141JgcIfHOYFjIiIiIiIiIkfEzVC7DBhhjMkzxkQCVwEvHs0bGGNmG2MeqK6u7pYCRUREREREpHfrqS19ngQWA6OMMduNMTdZa73AHcCrwFrgGWvt6qN5X2vtS9baW5OSkrq+aBEREREREen1emROrbV2zkGOzwfm90QNIiIiIiIi0ve4Ofz4uGn4sYiIiIiISP8W0qFWw49FRAQAvw9qd4HP63YlIiIiAuD3Q80Op43uZr1uSx8REZFOfF6o2+U0jDWlgdt97tfuBL8X7iiG9BFuVywiItL3WQt1u6Fqa+BrC1RuaX9cvQ18LfC1lZCa162lhHSoNcbMBmYPHz7c7VJERORYeFucQHrAwBp4XLcbrL/z6yJiITHb+Rp6cuD+IIhJdedziIiI9DXWQv1eJ6xWdQirbcG1eht4mzq/JjYdUnIhayKMuRCScyG6+0fVhnSotda+BLxUVFR0i9u1iIjIPlobO4TUg/Sy1pft/7qoxPbAmjHGCattobXteHQyGNPjH6m/MsbEAX8GWoA3rbWPu1ySiIgcL2uhoWL/0NoxuHobO78mJtUJrQPHwqhZTmhNzoXkIZA8GCLjXPkoIR1qRUTEJS31+wTVDoG1OtDL2lix/+uik9vDadakzkG17X50Yo9/nP7IGDMXuBAos9aO73B8FvAHwAM8ZK39FXAp8Ky19iVjzNOAQq2ISG9nLTRWdh4e3DG4Vm2FlrrOr4lOdgLqgJEw4uxAWB0SCK6DISrBlY9yOAq1IiLSWVP1oeev1pQ65+wrNs0JpUk5MHjaPr2rgyAxy7UruHJA84A/AY+2HTDGeID7gLOB7cAyY8yLQA7wUeC07l/xQ0REjkxTded5rPsG1+aazudHJToBNSUPhs3sEFoDXz0wVLg7hHSo1ZxaEZGj5G1x5sBUbITq7QcOrS21+7zIQHyGE05Th3Wew9rWy5qQDRHRrnwkOTbW2reNMUP3OTwNWG+t3QhgjHkKuBgn4OYAKznEzgnGmFuBWwGGDBnS9UWLiPQ3zbWdhwMHg2vg8b4XmSPj24cD557kDBXuGFpjUtz5HN0spEOt5tSKiBxAaxNUbnaC675f1ds6L7pkwiA+0wmmA0ZB/hmB3tZB7aE1PhPCI137ONKjBgHbOjzeDkwH/gj8yRhzAfDSwV5srX0AeACgqKjIdmOdIiJ9g98PNduhfD2Ub3Da7449ro2Vnc+PiG0PrYNPcG6DwTXXCa39cM2JkA61IiL9VnMdVG7qEFg3td/WlAId8kR0EqTmQ85UmHil09uamgdJgyF+IHjUFMihWWvrgRvcrkNEJCRZCw3lgeDa8WuD03Z3XEE4PLo9oA4qau9hTQksyBSb1i9D6+HoNxkRkd6qscpp7Co37RNcNzrb3HQUN8AJq3mnOPNkUoe1h9dYbXMjR6wUGNzhcU7g2BHT1CAR6bea66BiQ3tg7RhgOw4TDotw2ue04TD8TOc2bbhzATohU6H1GCjUioi4pW0p/QMNE67YuP/qwQnZTiM44uwOoXWYE2K1YrB0jWXACGNMHk6YvQq4+mjeQFODRKRP87Y4Q4T37XEtXw91uzqcaJwRUWn5MOGK9uCaNgyShmiUVBcL6T9NXQ0WkV7PWqdX9YDBdTM0d1zgIdAApubB2Iv3Ca5DITLWpQ8hfZEx5klgJpBujNkO/Mha+7Ax5g7gVZwtfeZaa1e7WKaISM/z+52pPPuG1vL1zlzXjmtTxKYHelzPcgJssNc1DyJi3PsM/UxIh1pdDRaRXqGt8esYWCs3tQ8Xbm1oP9d4nHkxqcNg8PT2ntbUYc7x8Cj3Pof0K9baOQc5Ph+Y38PliIj0rLbRUgec57qh8zzXiDgnsGZPgQmXd+517aOrCYeakA61IiI9xudt3wpn3/mtlZvB19x+ricyEFTzIO/U9rmtqcOcnlhPhGsfQ6S7aRSViPQqzXVOW33Aea5V7eeFhTttd9pwyD+9Q3AdrnmuIUChVkRkX03VsOsj2Plh+1f5evB7288Jj3FCavoIGHlu56HCidkQ5nGvfhEXaRSViPQ4b4szLHi/ea4boHZH53OD81w/3744U1q+s7Kw5rmGLP2fE5H+rb4cdn3YOcBWbGx/PiEbsibBqPOdRq9tqLCu2oqIiPQsa502eudK2LES9nziBNjKLWB97efFpnXoce0wzzUlT+tT9FEKtSLSf9Tu6hxed37oDCluk5zrBNjJ10DWZMiaCPEZrpUrEoo0/FhEuoTf76xPsWNFe4jduap9gUVPJKSPctrt8Zd1WKBpmLay64dCOtSq4RSRA7LWCav7Btjg3q7GafgGT4dptzoNYtZELfYg0gU0/FhEjprf36EHdkV7u91c4zzviYKB42DCZc5F5+zJMGAMhEe6WLT0JiEdatVwikjwSu7OlZ0DbGOl87zxwIDRkH9mILxOgszxEJXgatkiIiL9kt/vrC68Y2V7D+yuVZ0DbOZ4Z5Xh7MlOiM0Yo0UW5ZBCOtSKSD/j80L5Z/v0wK6CllrneU8kZIyFMRcFAuxkGDhW+8SJiIi4we935rwGhw+v7Nxuh0fDwPEw8YoOPbCjFWDlqCnUikjv5G2BPWs7B9hdH4O30Xk+ItZpCCdd1d4DO2C0hiKJuExTg0T6qU4BdkV7D2xLnfN8eDRkTnDa7bYe2AGjFGClSyjUioj7Whth9+rOQ4h3rwF/q/N8VKITWqfe1B5g04Zr2xyRXkhTg0T6Ab/PCbD7DiEOBtgYZwjxpDmQPcUJsemjtGWOdBv9zRKRntVcu/8esHvWtS/FH5PqNH4n3tEeYJOHQliYm1WLiIj0T34f7P1s/yHErfXO8+ExTg/s5KvbhxArwEoP0982Eek+DRXOlduOAbZ8ffvzCVlOaB0zuz3AJg7S/q8iIiJu8Ptg76f79MB+1B5gI2KdADvlC+1DiNNHKsCK6/Q3UES6Tu1u+OTfsPE/ToCt2tr+XPIQJ7ROusppBDMnQsJA10oVERHp13xeJ8B27IHd9RG0NjjPR8Q6bXXBtR16YEdq6o/0SiEdarUYhUgvULUN1r4Ea1+ErUsAC8m5kDMVpt4c2EJnojZCF+kn1DaL9EJ+P+z5pD3A7ljhBNiOiy9mToSCLzpzYLMmQ/oIBVgJGcZa63YNx62oqMgWFxe7XYZI/1G+Ada84ATZHSucYwPHO8OIx1zk7CenIcQSwowxJdbaIrfrCGVqm0VcZi1sL4bVz8Pqf0HtDud4RBxkTWzvfVWAlRBxqLY5pHtqRaSHWAtla2DNi06vbNlq53h2AZz1YyfIpuW7WqKIiEi/Zy2ULneC7JoXoHqbs4f78LNg9N2QU6TdA6RPUqgVkQOzFnYsd0LsmhehYgNgYMgMmPUrGH0hJA92u0oREZH+zVpnWPHqfzpfVVshLAKGnwln3A2jzoPoJLerFOlWCrUi0s7vg21LnWHFa19yrvAaD+SdCjO+4gRZLe4kIiLiLmud3QXagmzlZggLh2Gnw2l3wejzISbF7SpFeoxCrUh/52uFze86QfaTl6FutzNUKf8MmPk95wqvFnkSERFxl7Wwe3V7kK3Y4Fx4HnYanPJtGH2B2mvptxRqRfojbzNs+I8TZNfNh8ZKZ+XDEWc782NHnAPRiW5XKSIi0r9ZC2Vr24Ns+WdgwpwRVCd9DUbPhrg0t6sUcZ1CrUh/0VIP619z5sd++iq01EJUEoya5axanH8mRMa6XaWIhDht6SPSBfasg4+fd4Ls3nVOkM09CWZ82Qmy8QPcrlCkV1GoFenLmqqdALvmBVj/urMfXUwqjLsExl4MeadBeKTbVYpIH2KtfQl4qaio6Ba3axEJKXs/a++RLVsDGCfITrvFGUWlNS1EDkqhVqSvqS+HdS87PbIb3wR/K8RnwpQvwNiLYMiJ4NE/fREREdeVb2jfR3b3x86xITPgvN86bXZCpqvliYQK/WYr0hfU7IRP/u3Mkd38HlgfJA+B6bc5PbKDiiAszO0qRUREpGKjE2JX/9NZwRhg8HRnu7yxF0NitqvliYSikA61mrcj/VrlFmfbnbUvOtvwYCFtBJz8DWeYUtYkMMbtKkVERKRyc3uQ3bnSOTaoCM79byfIJuW4WJxI6AvpUKt5O9Lv7P3MmR+79kXY+aFzbOAEOP37TpDNGO1ufSIiIuKo2gZr/uUE2dIS51h2AZz9MyfIpuS6Wp5IXxLSoVakz7PWmWOz5kWnV3bPWuf4oCI4+6fOqsWpw9ytUURERBzVpe1Bdvsy51jWZDjrJ84ijSlD3atNpA9TqBXpbax1ruiufdEJs5WbnKX8h5wI5/3G2Vxdw5REpJfS1CDpd2p2OqOoVv8Tti1xjmVOgDN/CGMvgbR8V8sT6Q8UakV6A78Pti5xguzal6CmFMLCA5urfx1GX6g96UQkJGhqkPQLtbvbg+zWxYCFgePhjLth7OcgXRd1RHqSQq2Im1obYdGPnOX86/eAJwqGnwln3AOjZkFMitsVioiICEBdmXPxefW/YPO7gIUBY2Dm95yhxQNGuVygSP+lUCviFr8fnr/V6Zkde7GzH92IcyAqwe3KREREBKB+byDI/tMJstYP6SPhtO86QTZjjNsViggKtSLuef0nTkN5zi/gxDvcrkZERETAuei86ilY9QxsetvZ+z01H075Lxj3OcgYqy3zRHoZhVoRN5TMg/f+D4pughlfcbsaERERafP2b+DNX0JKnrP3+7jPOfNlFWRFei2FWpGetuEN+Pe3YPhZzmrGaiRFRER6h9X/dALtpDlwyV/URouEiDC3CxDpV3avgWeuc+bgfP4R8Oi6koiISK+wYwX880uQMw1m/0GBViSEKNSK9JTa3fDElRARC1c/DdGJblckIiIiALW74MmrIS4drnocwqPcrkhEjoK6iUR6QksDPHkVNOyFG+ZDUo7bFYmIiAg42+s9dTU0VcFNCyE+w+2KROQoKdSKdDe/H56/xRnWdNUTkD3F7YpEREQEwFp48atQWgJX/h0yJ7hdkYgcg145/NgYE2eMKTbGXOh2LSLH7bUfwif/hlm/hNHnu12NiIiItHn39/DRP+CMu2HMbLerEZFj1COh1hgz1xhTZoz5eJ/js4wx64wx640xd3V46rvAMz1Rm0i3Kp4L798LU2+B6be7XY2ISLczxsw2xjxQXV3tdikih7b23/D6T2H85+GUb7tdjYgch57qqZ0HzOp4wBjjAe4DzgPGAnOMMWONMWcDa4CyHqpNpHusfw1e/jaMOAdm/UqrKIpIv2Ctfclae2tSUpLbpYgc3K6P4PlbIbsALv6T2miRENcjc2qttW8bY4buc3gasN5auxHAGPMUcDEQD8ThBN1GY8x8a62/J+oU6TK7V8Mz10PGWPj8XG3dIyIi0lvUlcGTcyA6CeY8CRExblckIsfJzd+0BwHbOjzeDky31t4BYIy5Hth7sEBrjLkVuBVgyJAh3VupyNGo3QWPXwFR8c7WPVEJblckIiIiAN5mePoLUL8XblwACZluVyQiXaDXdh9Za+cd5vkHgAcAioqKbE/UJHJYLfXO1j2NlYGtewa5XZGIiIiAs9Lxv78J2z6Azz+i3QhE+hA3Q20pMLjD45zAMZHQ5PfBc7fAzg/hqiche7LbFYmIiEib9++FlY/DaXfB+EvdrkZEupCbW/osA0YYY/KMMZHAVcCLR/MGWmFRepVFP4R1LzuLQo2adfjzRUREpGd8+qrTTo+9GE77rtvViEgX66ktfZ4EFgOjjDHbjTE3WWu9wB3Aq8Ba4Blr7eqjeV+tsCi9xrKHYPGfYNptMP02t6sRERGRNmVr4dmbIHMCXPIXCHOzT0dEukNPrX485yDH5wPze6IGkW7z2SKYfyeMnAWzful2NSIiItKmvhyeuBIiY2HOUxAZ53ZFItINQvpSlYYfi+t2fQT/uB4GjofLHoYwj9sViYiICIC3BZ75orMrwVVPaPFGkT4spEOthh+Lq2p2Old/oxIDW/fEu12RiIiIgLPS8fxvw5Z34eI/QU6R2xWJSDfqtVv6iPRqLfXw5JXQWAU3vgKJ2W5XJCIiIm0++Css/xuc/C2YeIXb1YhINwvpnloNPxZX+H3w3M3O0OPLH4GsiW5XJCIiIm3Wvw6vfg9GXQBn3ON2NSLSA0I61Gr4sbhi4d2wbj6c9xsYea7b1YiIiEibPZ/CP26AjLFw6QNa6Vikn9C/dJGj8cEDsOTPMP1LMO0Wt6sRERGRNg0VztQgTwTMeVJrXYj0Iwq1Ikfq01fhle/CyPPg3F+4XY2ISI8yxgwzxjxsjHnW7VpE9uNrdXYjqNoGVz0OyUPcrkhEelBIh1rNqZUes3MVPHujs3H7ZQ9p6x4RCSnGmLnGmDJjzMf7HJ9ljFlnjFlvjLnrUO9hrd1orb2peysVOUavfA82vQWz/wBDTnC7GhHpYSEdajWnVnpEzQ5n657oJJijrXtEJCTNA2Z1PGCM8QD3AecBY4E5xpixxpgJxph/7/OV0fMlixyhZQ/Bsgdhxh0w5Rq3qxERF2hLH5FDaa6DJ66A5hq48VVIzHK7IhGRo2atfdsYM3Sfw9OA9dbajQDGmKeAi621vwQu7OESRY7Nxrdg/ndgxDlw9k/drkZEXBLSPbUi3crvg+dugt2r4fJ5kDne7YpERLrSIGBbh8fbA8cOyBiTZoy5H5hijPneQc651RhTbIwp3rNnT9dWK7Kv8g3wzBchfQRc9rCmBon0Y+qpFTmYV78Pn74CF/wPjDjb7WpERFxlrS0Hbj/MOQ8ADwAUFRXZnqhL+qmmanjyKjBhMOcpiE50uyIRcVFI99RqoSjpNkvuhw/ud+bnTL3Z7WpERLpDKTC4w+OcwDGR3s3ndRZvrNgIVzwKqXluVyQiLgvpUKuFoqRbrFsAr34PRl2g+Tki0pctA0YYY/KMMZHAVcCLx/umuuAs3W7RD2H9a3D+7yDvFLerEZFeIKRDrUiX2/khPHsTZE6Eyx7U/BwR6ROMMU8Ci4FRxpjtxpibrLVe4A7gVWAt8Iy1dvXxfi9dcJZutfxRWHIfTL8dim5wuxoR6SU0p1akTXWps3VPTApc/TRExrldkYhIl7DWzjnI8fnA/B4uR+TYbH4P/v0tyD8DzvmF29WISC+inloRgOZaJ9A218E1z0BCptsViYiEJA0/lm5RuRmeuRZScuHzj4BH/TIi0k6hVqRtwYmyNXDFPBg4zu2KRERCloYfS5drroUn54DfC3OehphktysSkV4mpEOtrgbLcbMWXrkLPlsIF/wOhp/ldkUiIiLSxu+D526BPevg8r9B+nC3KxKRXiikQ62uBstx++B+WPYgnPhVKLrR7WpEREKeLjhLl3r9p/DpAjjv15B/utvViEgvFdKhVuS4fPIyvPI9GH0hnKWte0REuoIuOEuXWfkkvPd/zkVn7RkvIoegUCv9044V8NzNkD0FLn0QwvRPQUREpNfYthRe+hoMPQXO+w0Y43ZFItKL6Td56X+qt8MTV0FsGsx5CiJj3a5IRERE2lRtg6euhsRBcMWj4IlwuyIR6eW0Hrr0L0018PgV0NoA174KCQPdrkhERETatNTDU3PA2wzXvwyxqW5XJCIhQD210n/4vPDsDbDnE7jibzBwrNsViYj0OVooSo6Z3w//vA12r4bPz4UBo9yuSERChEKt9A/WwoI7Yf1rcOHvIf8MtysSEemTtFCUHLM3fwlrX4Jzfg4jzna7GhEJISEdanU1WI7Y4vugeC6c9HUovN7takRERKSjj56Ft38DU66FE77sdjUiEmJCOtTqarAckbX/hoV3w5iL4Mwfu12NiIiIdFRaAi98BYacCBf8Xisdi8hRC+lQK3JYpcudrXsGFcClD2jrHhERkd6kZgc8eTXEZ8CVj0F4pNsViUgI0m/40ndVbYUnr4K4Ac7WPRExblckItLnaWqQHLGWBmfrnpY6p52OS3e7IhEJUQq10jc1VcMTV0JrE1zzD+cKsIiIdDtNDZIjYq0z5HjHSrjsIRg4zu2KRCSEaZ9a6Xt8rfCP62Hvp3DNs5Ax2u2KREREpKO3fwurn4ezfgKjznO7GhEJcQq10rdYC/PvhA1vwOw/Qv7pblckIiIiHa15Af7zC5h4lbMrgYjIcdLwY+lb3r8XSh6Bk78Jhde5XY2IiIh0tPND+OftkDMVZv9BKx2LSJdQqJW+Y82LsOiHMPYSOOOHblcjIiIiHdXuhifnQEwqXPk4RES7XZGI9BEKtdI3bC+B52+FnCL43P3aukdExCVa/VgOqLUJnr4GGithzpOQMNDtikSkD9Fv/hL6KrfAk1dC/AC46klt3SMi4iKtfiz7sRZe+jpsXwaf+ytkTXS7IhHpY0I61OpqsNBYBU9cAd4WZ6Xj+AFuVyQiIiIdvfd/sOopOP1uGHuR29WISB8U0qFWV4P7OV8r/OM6KF8PVz4GA0a5XZGIiIh09Ml8eO0nMP4yOPXbblcjIn2UtvSR0GQtvPwt2PgmXHwfDDvN7YpERESko92r4flbIHuy01ZrpWMR6SYh3VMr/dh7f4Dlj8Ip/wVTvuB2NSIiItJR/V544iqIStB6FyLS7dRTK6Fn9b/gtR/BuEud+TkiIiLSe3hb4OkvQH0Z3LAAErPcrkhE+jiFWgkt25bBP2+DnGlwyV+0dY+IiEhvYi28/E3Yuhg+PxcGFbhdkYj0A0oEEjr8PnjuRogf6Oxxp03bRUR6He1M0M8t+TOs+Duc+h1ncSgRkR6gUCuh49NXoWornPMziEt3uxoRETkA7UzQj322CBbeDWMugpnfc7saEelHFGoldCx7EBKyYNQFblciIiIiHe1ZB8/eCAPHwefu1/QgEelRmlMroaF8A2x4A2Z+Hzz6aysi/Y8xJgyYBGQDjcDH1toyd6sSARoq4IkrITwa5jwFkXFuVyQi/YzSgYSGZQ9DWDgUXud2JSIiPcoYkw98FzgL+AzYA0QDI40xDcBfgb9Za/3uVSn92vO3Qs0OuP5lSMpxuxoR6YcUaqX3a2mAlX935ugkZLpdjYhIT/s58BfgNmut7fiEMSYDuBq4FvibC7VJf7d7NaxfBGf9GAZPdbsaEemnFGql9/v4WWiqhqk3u12JiEiPs9bOOcRzZcD/9Vw1IvsofgQ8UVCgkVQi4h7N4pfezVpY+iBkjIXcE92uRkTENcaYy40xCYH79xhjnjfGaBNQcU9zHXz4FIy7BGJT3a5GRPqxXhdqjTFjjDH3G2OeNcZ8ye16xGXbi2HXKph6ExjjdjUiIm66x1pba4w5GTgTeBhnWLKIOz5+DlpqoehGtysRkX6uR0KtMWauMabMGPPxPsdnGWPWGWPWG2PuArDWrrXW3g5cAZzUE/VJL7bsQYhMgIlXul2JiIjbfIHbC4AHrLUvA5Eu1iP9XfFcZyTV4OluVyIi/VxP9dTOA2Z1PGCM8QD3AecBY4E5xpixgecuAl4G5vdQfdIb1e+F1f+EyXMgKsHtakRE3FZqjPkrcCUw3xgTRS8ccSX9ROly2LkSCm/QSCoRcV2PNIbW2reBin0OTwPWW2s3WmtbgKeAiwPnv2itPQ+4pifqk15q+aPga4Gim9yuRESkN7gCeBU411pbBaQCd7pa0QEYY2YbYx6orq52uxTpTiWPQEQsTNJIKhFxn5tXeAcB2zo83g4MMsbMNMb8MXA1+qA9tcaYW40xxcaY4j179nR3rdLT/D5nRcWhp0DGaLerERFxjTEm1RiTirM37ZtAeeBxM1DsZm0HYq19yVp7a1JSktulSHdpqoaPnoXxl0G0/j+LiPt63ZY+1to3cRrtw533APAAQFFRkT3M6RJqPlsI1VvhnJ+5XYmIiNtKAAsYYAhQGbifDGwF8lyrTPqnVc9AawMU3eB2JSIigLs9taXA4A6PcwLHRJxtfBKyYPQFblciIuIqa22etXYY8Bow21qbbq1NAy4EFrpbnfQ71jojqbImQbZ2lBKR3sHNULsMGGGMyTPGRAJXAS8ezRto3k4fVb4BNrwOhdeDJ8LtakREeosTrLXBaTnW2gWANvCWnrVtKZStdrbx0QJRItJL9NSWPk8Ci4FRxpjtxpibrLVe4A6cRS/WAs9Ya1cfzftq3k4fVTwXwsKdUCsiIm12GGPuNsYMDXz9ANjhdlHSz5Q84my1N/7zblciIhLUI3NqrbVzDnJ8Ptq2RzpqaYAVf4cxsyEh0+1qRER6kznAj4B/Bh6/HTgm0jMaKuDj52HKFyAq3u1qRESCet1CUUfDGDMbmD18+HC3S5Gu8vFz0FQFU292uxIRkV7FWlsBfN3tOqQf+/BJ8DVrgSgR6XVCetN2DT/uY6yFZQ/CgDGQe5Lb1YiI9CrGmJHGmAeMMQuNMW+0fbldl/QTbQtE5UyDzAluVyMi0klI99RKH1NaAjs/hPN/p8UnRET29w/gfuAhwOdyLdLfbH4Xyj+DS/7idiUiIvsJ6VCr4cd9zNIHncUnJl3ldiUiIr2R11qrRCHuKJ4L0Ukw7nNuVyIish8NP5beob4cVj/vBNqoBLerERHpjV4yxnzZGJNljElt+3K7KOkH6vbA2pdg8jUQEeN2NSIi+wnpnlrpQ1Y8Cr4WmHqT25WIiPRW1wVu7+xwzALDXKhF+pOVfwd/q7baE5FeS6FW3Of3OcOahp4CGWPcrkZEpFey1ua5XYP0Q36/s0BU7skwYJTb1YiIHFBIDz+WPuKzRVC1Vb20IiKHYIyJMMZ8zRjzbODrDmNMhNt1SR+38Q2o2qJtfESkVwvpUGuMmW2MeaC6utrtUuR4LHsQErJg9IVuVyIi0pv9BSgE/hz4Kgwc6xHGmEuMMQ8aY542xpzTU99XXFb8CMSmw5jZblciInJQIR1qtVBUH1CxEda/5szT8ajDQUTkEKZaa6+z1r4R+LoBmHokLzTGzDXGlBljPt7n+CxjzDpjzHpjzF2Heg9r7b+stbcAtwNXHvOnkNBRswPWLYAp10B4lNvViIgcVEiHWukDlj0MYeFQcN3hzxUR6d98xpj8tgfGmGEc+X6184BZHQ8YYzzAfcB5wFhgjjFmrDFmgjHm3/t8ZXR46d2B10lft/wxsD4tECUivZ4WihL3tDbCir87w44Ts9yuRkSkt7sT+I8xZiNggFzgiCY6WmvfNsYM3efwNGC9tXYjgDHmKeBia+0vgf3mgxhjDPArYIG1dvkxfwoJDT4vLP8b5J8BqVpgW0R6t5AOtcaY2cDs4cOHu12KHIuPn4OmKph6s9uViIj0etba140xI4C2JWjXWWubj+MtBwHbOjzeDkw/xPlfBc4Ckowxw6219+97gjHmVuBWgCFDhhxHaeK69YugphRm/crtSkREDiukhx9rTm0IsxaWPggDxsDQk92uRkSk1zPGfAWIsdaustauAmKNMV/uqe9vrf2jtbbQWnv7gQJt4JwHrLVF1tqiAQMG9FRp0h2K50J8Jow6z+1KREQOK6RDrYSw0uWwc6WzjY8xblcjIhIKbrHWVrU9sNZWArccx/uVAoM7PM4JHJP+rmqrs91ewRe1iKOIhASFWnHHsgchMh4magFNEZEj5AnMawWCCz1FHsf7LQNGGGPyjDGRwFXAi8dZo7bb6wtK/uZccC74otuViIgcEYVa6Xn15fDx8zDpKohOdLsaEZFQ8QrwtDHmTGPMmcCTgWOHZYx5ElgMjDLGbDfG3GSt9QJ3AK8Ca4FnrLWrj7dITQ0Kcb5WWP4ojDgHkgcf/nwRkV4gpBeKkhC14jHwNWuBKBGRo/Nd4DbgS4HHi4CHjuSF1to5Bzk+H5jfJdVJ3/DJy1BfBkU3ul2JiMgRC+lQq9WPQ5Df5yw+kXsyZIxxuxoRkZBhrfUbY+YBb1hr17ldz8GobQ5xJY9A0mAYfpbblYiIHLGQHn6sIU4haP1rULXFWSBKRESOmDHmImAlgSHHxpjJxpjjngPb1dQ2h7DyDbDxTSi4DsI8blcjInLEQjrUSgha+qCzRcCY2W5XIiISan4ETAOqAKy1K4E8F+uRvqbkETAeKLjW7UpERI6KQq30nIqNTk9t4fXaIkBE5Oi1Wmv3XVLYulKJ9D3eZljxOIy+ABIy3a5GROSoKNRKzymeCyYMCq9zuxIRkVC02hhzNc7WPiOMMfcC77td1L60pU+IWvMiNFZA0Q1uVyIictQUaqVntDbCir/DmAshMdvtakREQtFXgXFAM852PjXAN9ws6EA0pzZEFc+FlDzIm+l2JSIiRy2kVz+WEPLx89BYqW18RESOkbW2AfgB8ANjjAeIs9Y2uVyW9AVln8DW9+Hsn0KY+jtEJPToJ5f0jGUPwYDRMPQUtysREQlJxpgnjDGJxpg44CNgjTHmTrfrkj6g5BHwRMLka9yuRETkmIR0qNW8nRBRWgI7lju9tMa4XY2ISKgaa62tAS4BFuCsfNzrlqlV2xxiWhpg5ZMw5iKIS3e7GhGRY3LEodYY8/MDHHN1EzPN2wkRSx+CyHiYeKXblYiIhLIIY0wETqh90VrbSi9c/Vhtc4hZ/U9oroaiG92uRETkmB1NT+0gY8yctgfGmAzgta4vSfqUhgr4+Dkn0EYnul2NiEgo+yuwGYgD3jbG5OIsFiVy7IrnQvooyD3R7UpERI7Z0SwUdRvwqjFmA86V4UeA73ZLVdJ3rHgMfM1aIEpE5DhZa/8I/LHtsTFmK3C6exVJyNu5CkqLYdavND1IRELaYUOtMeZRYDmwAvgK8ATgBS6x1q7v3vIkpPn9sOxhyD0JBo51uxoRkZBkjPkC8IS11t/xuLXWAl5jTD6QZa1915UCJXSVPALh0TDpKrcrERE5LkfSUzsPmATcAEwEhgLLgC8YYz621j7bbdVJaFv/GlRtgbN+7HYlIiKhLA1YYYwpAUqAPUA0MBw4DdgL3OVeeZ0ZY2YDs4cPH+52KXIozbWw6hkYfxnEpLhdjYjIcTlsqLXWvgG80fbYGBMOjMEJutMBhVo5sGUPQvxAGH2h25WIiIQsa+0fjDF/As4ATsK5wNwIrAWutdZudbO+fVlrXwJeKioqusXtWuQQPvoHtNRB4Q1uVyIictyOZk4tANZaL87+eB8Bf+/yiqRvqNgEny2C074D4ZFuVyMiEtKstT5gUeBL5PhY6ywQNXAC5BS5XY2IyHEL6X1qpRcrngsmDAqvd7sSERER6ah0Oez6CIpu0AJRItInKNRK12ttdFY9Hn0BJGa7XY2IiIh0VDwXIuJgwuVuVyIi0iUUaqXrrf4nNFZqGx8REZHeprEqsH/85do/XkT6jJAOtcaY2caYB6qrq90uRTpa9pCzkXveqW5XIiLSZxhjBhpjHjbGLAg8HmuMucntuvaltrmXW/U0eBuh6Ea3KxER6TIhHWqttS9Za29NSkpyuxRpU7ocSkucXlrN0xER6UrzgFeBtnkdnwLfcKuYg1Hb3Iu1LRA1qBCyJrldjYhIlwnpUCu90LKHnHk6k650uxIRkb4m3Vr7DOCH4G4EPndLkpCydTHs+UTb+IhIn6NQK12nocKZpzPpSojWFXoRkS5Wb4xJAyyAMeYEQGN85cgVPwJRSTD+UrcrERHpUke9T63IQa34O3ibtECUiEj3+BbwIpBvjHkPGAB83t2SJGTUl8Oafzm9tJFxblcjIn2QtZaqhlb21DWzp7aZstom9tQ2c/X0XOKjujd2KtRK1/D7ofhhGHIiDBzndjUiIn2OtXa5MeY0YBRggHXW2laXy5JQsfJx8LU4e9OKiByFxhYfe2qb2VPnhNS2r7K2+4EQu7eumVaf3e/1p43MYFRmQrfWqFArXWPD61C5Gc78oduViIj0ScYYD3A+MBSn/T7HGIO19veuFia9n98PJfNgyAzIGON2NSLSC/j8lvL6fcLpviE18Fxds3e/1xsDaXFRDEhwvkYOTHDux7cfa/tK6OZeWlCola6y9EGIy4DRs92uRESkr3oJaAI+IrBYlMgR2fw2VGyAmXe5XYmIdCNrLbXN3s4BNRBSy2raw+qe2mYq6pvx79+pSkJUOAMSokhPiGJMdiKnxkeRkbh/WE2NjSTc03uWZ1KoleNXuRk+Wwin3gnhkW5XIyLSV+VYaye6XcThGGNmA7OHDx/udinSpnguxKTCmIvcrkREjkGL18/eugP1qDbtNxS42bv/Nc8IjwmG0kHJ0UwenNQhpEYzICGKjIQo0uOjiIn0uPAJj59CrRy/4rlgwqDwercrERHpyxYYY86x1i50u5BDsda+BLxUVFR0i9u1CFC7Gz55GabfDhHRblcjIkCz10dVQ2vgq4XKhlaqG53bvbWde1T31DVT1XDg5RNS4yKD4XTq0LgDDv/NSIgiKSYCY0wPf8qepVArx6e1CZY/BqPPh6RBblcjItKXLQH+aYwJA1pxFouy1tpEd8uSXm3FY+D3am9akW7Q7PVR3dBKVaMTUCsbWgKPnYBa1RZW651zqgMBtrH14FuMR0eEkZEQTUZCFMMz4pmRn3bAeappcVFEhvee4b9uU6iV47P6n9BYAVN1QV5EpJv9HpgBfGStPcBMKJF9+H1Q8jfIOw3SNRxc5GBavH6qGluCAbWyviUQQp2g6oTWlk69q1WNrTS0HDycRngMSTGRpMRGkBwbwaDkGMZnJ5IcG0FybKRzG3g+qe1YTASxkZ4+36vaHRRq5fgsewjSR0LeqW5XIiLS120DPlaglSO2/nWo3grn/NTtSkR6RKvPH+wddXpOO4TRxg6htLFzQK0/RDgNDzPtQTQmguzkaMZmJ5IcE7FfQE0OBNiU2EiF0x6mUCvHbscKKC2G837jrOstIiLdaSPwpjFmAdDcdlBb+shBlTzi7Eww6gK3KxHpEj6/5dPdtazcVsWH26rYXtnoDPWtb6W6sfWAW8+08YQZkmOcXtGU2EgyE6MZnRnoOY2JIDkuMhhUU2IjSQrcj48KVzgNAb0y1BpjLgEuABKBh3v7ohj91tKHICIOJl3ldiUiIv3BpsBXZOBL5OCqt8Onr8DJ39TOBBKydlY3snJrFSu3VbFiWxUfl1YHh/wmxUQwbEAcGQnRjMxICPaYOsN5nYCaEjiWFBtBgsJpn9ZjodYYMxe4ECiz1o7vcHwW8AfAAzxkrf2VtfZfwL+MMSnA7wCF2t6moQI+fhYmzYHoJLerERHp86y1P3G7Bgkhyx8Fa6HgOrcrETkidc1eVm13AuzKrVV8uL2K3TXOoJQIj2FsViKXF+YweUgykwenMDQtViFVgnqyp3Ye8Cfg0bYDxhgPcB9wNrAdWGaMedFauyZwyt2B56W3Wfk4eJtg6s1uVyIi0qcZY/7PWvsNY8xLwH7zaa212nxUOvN5nVA7/CxIyXW7GpH9eH1+Pt1d5wTYbZWs3FbFZ2V1tK0YMDQtlhOGpTF5cDKTByczNjuRqPDQ3D9VekaPhVpr7dvGmKH7HJ4GrLfWbgQwxjwFXGyMWQv8ClhgrV3eUzXKEfL7YdnDMGQGZI4//PkiInI8Hgvc/s7VKiR0fPoK1O6ECzTdWtxnrWVndVMgwDq9sB+VVge3tUmOjWDy4GTOG5/l9MLmJJMSpyHzcnTcnlM7CGc1xzbbgenAV4GzgCRjzHBr7f37vtAYcytwK8CQIUN6oFQJ2vAGVG6CM+52uxIRkT7PWlsSuDvZWvuHjs8ZY74OvNXzVUmvVjwXErJhxDluVyL9UG1TKx9tr2ZFW4jdVsWeWmcYcaQnjLHZiVw5dXCwFzZXw4ilC7gdag/IWvtH4I+HOecB4AGAoqIibW/Qk5Y95KymOEYj3kREetB1OGtQdHT9AY65yhgzG5g9fLj2RXVFxSbn4vPMu8DTK3/Nkz7E6/OzLrAacduCTuv3tA8jzkuP4+Th6UwenMykwcmMyUrQMGLpFm7/tCsFBnd4nBM4Jr1V5RZnWNOp39ZqiiIiPcAYMwe4GsgzxrzY4akEoMKdqg7OWvsS8FJRUdEtbtfSLy3/G5gwKPii25VIH2OtpbSqkQ+3VQfnwX5UWk1Tqx+AlMAw4gsnZjN5SDKTcpJIjtXvitIz3A61y4ARxpg8nDB7FU7DfUR0NdgFxXOdPWkLr3e7EhGR/uJ9YCeQDvxPh+O1wCpXKpLeydsCyx+DUedBYrbb1UiIqwkMI165rYoVgV7YvXWBYcThYYzLTmTOtCFMHpzMlMEpDE6N0TBicU1PbunzJDATSDfGbAd+ZK192BhzB/AqzpY+c621q4/0PXU1uIe1NsGKx2DU+ZCU43Y1IiL9grV2C7AFmOF2LdLLffISNOyFwhvcrkRCTKvPz7pdte2LOW2rYkOHYcTD0uM4dUR6YDudZEZnJhIZHuZu0SId9OTqx3MOcnw+ML+n6pDjsOZf0FAO03QNQUSkpxljLgV+DWQAJvBlrbWJrhYmvUfxI5A8BPLPcLsS6cWstWyvbGTltio+DATYj3e0DyNOjYtk8uBkLpqU7cyFzUkmKTbC5apFDs3t4cfHRcOPe9iyhyBtBOSd5nYlIiL90W+A2dbatW4XIr3Qnk9h8ztw5o8gTD1o4oTX8voWSisbKa1qZOOeumAv7N66FgCiwsMYPyiJq6flMnlIMlMGJ5OTomHEEnpCOtRq+HEP2rESti+DWb925tSKiEhP261AKwdVMg/CwmHKF9yuRHqIz2/ZXdNEaVVjMLhur2xge+D+jqrGYO9rm/wBcZw2MoPJg5OYPDiF0VkJRHh0EURCX0iHWulByx6EiFiYdJXblYiI9FfFxpingX8BzW0HrbXPu1aR9A6tjfDhEzBmNsRnuF2NdJFmr4+dVe2hdXvbbWUDpVWN7KpuwuvvvKtlWlwkg1JiGDUwgTNGZTAoJYaclFgGJccwODWGhGgNI5a+SaFWDq+xEj561gm0McluVyMi0l8lAg3AOR2OWUChtr9b84LTVhfd6HYlchTqm737BVbnsdPbuqeuObhQEzgD5TIToxmUHENhbgqDkmM6hdZByTHERGoPWOmfQjrUak5tD1nxOHibYOrNblciItJvWWu1pK0cWPFcSBsOQ09xuxIJsNZS3djK9srG4HBgJ7Q2BO9XNrR2ek2Ex5AdCKenjRzAoJSYYHAdnBJLZlK0hgqLHERIh1rNqe0Bfj8UPwyDT4DMCW5XIyLSbxljHsHpme3EWqvuuf5s92rY9gGc8wutedGD/H7L3rpmtu0bWDs8rm/xdXpNbKQnGFIn5SQHQ2tOoLd1QHwUYWH6fyhyLEI61EoP2PgGVGyE03/gdiUiIv3dvzvcjwY+B+xwqRbpLYofAU8UTL7a7Ur6lFafn13VHeaz7tPLuqOqiRZf50WYkmMjGJQcw9C0OE4anh4MrIOSYxmUEkNKbIRWFRbpJgq1cmjLHoa4Ac7iEyIiPazF66emqZWaxlZqmrxUN7bdbw3c9wafr25s5Tefn0hWUozbZXcLa+1zHR8bY54E3nWpHOkNWuph1dMw7nMQm+p2NSGjbaubHVVOON1Z7awUvKO6iZ2BY2W1TeyzBhMDEqLISYlh/KAkzh2fSU6g17UttMZH6ddqEbeE9L8+zantZlVb4dNX4ORvQXiU29WISAjy+S11bWG0Q/h07ns7hFMntO77fGOr75DvH+ExJMVEkBgdQWJMxH7bV/RxIwAtddufffwcNNdAkaZbd1Tb1MqOqiZ2VDeys6opEFgD96sb2VndRIu388+KqPAwspNjyEqK5uQR6WQnRXcKrFlJ0URHaBEmkd4qpEOt5tR2s+K5zm3h9a6WISLusdZS3+JrD5sHCZ8HCqc1ja3UNnsP+f5hBhKDoTScpJgIhifEd3rc9rxzPzwYYJNiIogKD+s3w/mMMbV0nlO7C/iuS+VIb1A8FzLGwuDpblfSY5pafeyq3jewtve27qxq2u/njifMMDAhiqzkGCbmJDNrXDRZSdFkJ8cEvzQ0WCS0hXSolW7kbYblj8Ko8yF5sNvViMgxstbS2OoLBs/aps5Ddg83pLemyYtv3zF4+4iPcsJnQnQ4iTHOnLKxWYnBANoeTMODYbTtcXxUuH6RPALG+UMaZ63d6nYt0kuULocdK+D83/WZBaJ8fktZbVPnIcHB+87t3rqW/V6XFhdJVnI0Q9PiODE/vUNgjSYrKYaMhCjCtWqwSJ+mUCsHtvpf0FCubXxEXOb3W+pavNR26P2safIGwmmHXtGmVuecfUJrbZMX72FCaXREWKfwmR4fybABcR2G9e4bTtt7TeOjwvXLYg+w1lpjzMuAK8vQG2PGAF8H0oHXrbV/caMO6aDkEYiIhYlXuF3JEbHWUtnQ6vSmVu8zJDhwbFdN034X0eKjwoMhdfygRLKSAr2rSdFkJWtYsEhvt7dxL+kx6d3+fUI61GpObTda9pCz513eaW5XIhLSvD4/dc3e9qB5gNB58GPO8F176ExKbKQnGD4TottDacdjHcNpQrRzrC2YRoXrF8IQsdwYM9Vau+xoXmSMmQtcCJRZa8d3OD4L+APgAR6y1v7qYO9hrV0L3G6MCQMeBRRq3dRUDR89B+Mvg+gkt6sBoL7Zy87qRkqrAostVbeF1fa5rPvOeY/0hJGZFE12cjTT81LJSg70sCbFBO8nRke49IlE5Fg0+5op2V3C+6Xv896O91hftZ5XL3uV7Pjsbv2+IR1qNae2m+z8ELYvhVm/gjD1wEj/5vdbqhtbqeowRLemMdBT2iGM1nboMe14rO4wc0oBEqKcYbltw3ezk2MYHZ3QachuWxjtOP80IXAsQj2lWGvZVLOJvMS8vjyceTpwjTFmC1APGJxO3ImHed084E84YRQAY4wHuA84G9gOLDPGvIgTcH+5z+tvtNaWGWMuAr4EPNYFn0WOx6pnoLUeinp2i+Lqxla2lNezaW89m/c2sDlwf0t5PZUNrZ3ONQYyEqLITo5hTFYiZ4zO6DQkODs5hrS4SO3LKhLirLVsqt7Eezve470d71Gyq4QmXxMRYREUDCzgovyLiA6P7vY6QjrUSjdZ9pAzpGnSHLcrEelyfr+lqrGVivpm9ta1UFHfQnl9CxV1LZTXNwfvO8ebqWxoPeSc0o4LHSVEOz2hQ9Njg4sZtR3rGFA7HouPCsejX+qOyZ6GPSzZucT52rGEssYyXrjkBYYlDXO7tO5y7rG8yFr7tjFm6D6HpwHrrbUbAYwxTwEXW2t/idOre6D3eRF4MTAM+okDnWOMuRW4FWDIkCHHUq4cjrXO3rRZk2BQQZe/fW1TazCwbt5bz6bA7ebyBirqO89nzU6KZmh6HLPGZzEkNbZDYI1mYGK0LriJ9FE1LTUs2bGE93e8z/s73mdn/U4AhiYO5bKRl3Fi9okUDSwiNiK2x2pSqJXOGith1T+cOToxyW5XI3JYPr+lqqE9nJbXtVDRFk4Dj8vrm4P3Kxta9tt7sE1STARpcZGkxkWSmxZLQW5K8HFKXESHcNreaxoX6enLPYO9SkNrA8W7i1mycwmLdyxmfdV6AJKjkpmeNZ0ZWTNIi05zucruY63d0oVvNwjY1uHxdpye4AMyxswELgWigPkHO89a+wDwAEBRUdFhBs7LMdm+DMpWw+w/HPNb1DV7A0G1PbC2Pd53IabMxGiGpsdy7riBDE2LY2h6HEPT4shNi9VcVpF+wuf3sbp8tdMbW/oeH+39CL/1Ex8RzwlZJ3DLxFs4MftEBsUPcq1GhVrpbOUT4G3UAlHiGp/fUtlw4EDa1nvadr+i/ghCanwkaXGR5KXHUZibSnq8E1JT4yJJi4sKPp8SF6lehV6mrRFdvGMxS3YuYeWelXj9XiLDIikYWMCFwy5kRvYMRqeOJszo/113sta+CbzpchkCzjY+kQkw/vOHPK2hxdtpiPDmvfVsKW9gU3k9e2qbO507MDGK3LQ4zhw9kKHpceSlxzI0PY7c1DhiIhVcRfqj3fW7eX+HMy92yc4lVDdXYzCMSxvHzRNu5qTsk5gwYAIRYb1j3rtCrbTz+2HZw85+d1mHm6IlcmS8Pj+VDa3BQNoeVp0e1Yr6luAw4LaQerCFkZJjIwJhNJL8AfFMzXPup8VFkhofFexVTYuPJCVWITXUWGvZVruNxTsWs3jnYpbuWkptSy0AY1LHcO3Ya5mRNYMpGVN6ZH5OH1cKdNyvLSdw7LhpEcdu1FABHz8PBddCVDyNLT62VASGCe9tCA4X3lJez+6azsF1QEIUeWlxzBw5IBBc4wI9r7HERurXQZH+rsnbxPLdy3lvx3u8v+P94GioATEDmJkzk5MGncQJWSeQEp3icqUHFtI/xdRwdrGN/4GKDTDze25XIr2ctZaaJi97apvYXdNMWdttTTO7a5vYU9tMeZ0TWKsaWw8aUlOCITWK4QPiSQuEVCeYBkJqoGdVIbVvqmyq5IOdHwSHFO+o3wFAVlwWZ+eezYysGUzLmkZqdKrLlfY5y4ARxpg8nDB7FXB1V7yxFnHsWk2tPraUOz2uCSse4ERfM9/ZXMDb//06u2qaOp2bHh/J0LQ4ThkxgLx0Z4hw25Dh+KiQ/pVPRLqYtZaN1Rt5r9QJscW7i2n2NRMRFkHhwEIuzr+YEwedyIjkESExzSqkf8Kp4exiyx6G2HQYe5HblYhLrLXUNHrZXdvkBNSaJspqnds9te2Py2qb9tuaAZytZQYmRjMgPoqRAxM6h9NAj2pafFQgpEZof9N+qNnXzPLdy1m8czFLdizhk4pPsFgSIhKYljWNG8ffyAnZJzAkYUhINKKhwBjzJDATSDfGbAd+ZK192BhzB/AqzorHc621q10ss19ravWxraLBGSZc3mGO6956dtY0BS4MWl6PfIIPw0ay3gzlxOFx5HWY4zo0PZYEbX8jIodQ3VzNkp3tCzztqt8FQF5SHpePvNxZ4CmziJjwGJcrPXohHWqlC1Vtg08XwMnfhPAot6uRLmatsy1Np17VQHDd93Gzd/+wGh8VTkZCFBmJUUwenMzAxCgyEqLJCNwOTIwiIzFaPQGyH7/1s65iHYt3LmbxjsWsKFtBs6+Z8LBwJg2YxFcmf4UZ2TMYmzaW8DD9/ekO1toDLmVvrZ3PIRZ9OlYaRXVw9c1ePtxWxZqdNYGtcJwgu6O6sdOIlpTYCHLT4pg+LC0YWMe3rCJ//k645Ec8P/kk9z6EiIQMn9/Hx+UfB/eMbVvgKSEigROyT+C2ibdxYvaJ3b6HbE/QbxDiKJ7r3Bbe4G4dclSstVQ1tO7Xs1rWoYfV6VltpuUAYTUhKjwYTAuHpJCRGB0Ir87twMBtnMKqHIUddTuCw4k/2PkBlc2VAAxPHs7lIy9nRvaMHl/qX3qORlE5rLVsKW9g+dZKlm+tpGRLFet21QQXtkuKiWBoehxTh6aQm5bjzHFNd3pfk2IP0OP67D0QnQzjLunJjyEiIWZX/S5ngadSZ4GnmpYaDIbx6eO5ZcItnDToJCakT+hzF5L71qeRY+NthuWPwsjzIHnw4c+XbmetpbKh9cDDfwPzVstqmtlT20yL7wBhNTo8GEinDk09YFDNSIzS4iDSJWpaali2c5kzpHjnErbUODvPDIgZwCk5p3BC1gmckHUCA2IHuFypSPdpaPGyans1JVsqWbG1khVbqygP7OsaHxXO5MHJ3HH6cApyU5iYk0xqXOSRv3ndHljzIky7BSJCb1igiHSfJm8TJbtLnAWeSt9nQ/UGADJiMjhjyBmclO0s8JQcnexuod1Mv9EKrHkBGvbCNG3j0xN8fktZbROllY1sr2yktKqRXdVNwWHAewJzVlt9+6+ulBQTEQym0/PiOvSsdgirCdHagkG6VauvlQ/3fBicF/tx+cf4rZ+Y8BimZk7lqlFXMSN7BsOShmlerPRJ1lq2VTQGe2GXb61k7c5afIFu2GHpcZw+OoOCISkU5CYzIiMBT9hx/FtY+Xfwt2o0lYhgrWVD1YbgKsUlu0to9jUTGRZJ4cBCPjfic5yYfSLDk4f3qzZYoVZg2UOQmg95M92upE/w+S27a5rYXtnI9soGJ7hWNrK9yrm/o6pxv8CaHNseVocNiAsG1I63AxKitNG9uMJay/qq9cH9Yot3F9PobcRjPMHhTDOyZzAxfSIRHi1U09/1xTm1Ta0+Vm2vdgLslkqWb61ib52zZU5spIfJg5P50mn5FOQmM2VwCilH0wt7OH4/lMyD3JNhwMiue18RCRnVzdUs3rmY90udBZ52N+wGYFjSMC4feTknDTqJwoGFIbnAU1dRqO3vdq6CbR/Aub+EMK1EeyS8Pj+7gqE1EFgD4XV7VQM7q5rw+juH1oyEKHJSYpiYk8z5E7LISYkhJyWWnJQYBiXHKKxKr1PWUBacF7tk5xL2Nu4FYGjiUC7Ov5gZ2TOYmjmVhMgElyuV3ibU59RaaymtamT51iqWB4YSr95RE/y5PjQtllNHpDMlN4WCIcmMGpjQvSu5b/wPVG6GM+7pvu8hEoJafC1UN1cTZsIIDwsnzIThMR48YZ7g/TATmr/bev1ePt77sTM3dsd7fLzXGRGVEJnACVkncFL2SZyYfSJZ8Vlul9prhHSo7YtXg3vcsocgPAYmH3BxzH6p1ednV3Xnnta2+6VVjeysbgoOMQMwpi20xlIwJIWcSU5gHZQcQ05KDNkKrRIC6lvrKdldwuIdzirFbXNyUqNTmZ45nRnZMzgh6wQ1oNLnNLX6WL2jmuVbqgILOlVSVuv0wsZEeJiYk8Qtpw6jYEgKU4Ykkx7fwzsEFM91ttsbM7tnv69IL9Tqb2XJjiUs2LSAN7a9QX1r/SHPN5hguPWEeYL320JwmAkj3ITv93xbOPYYz36vP+D7mXDCwsI6P9/hPQ54/gFq8vl9rNyzkiU7l1DbUkuYCWN82nhunXgrJ2WfxPj08X1ugaeuEtJ/KqF+Ndh1jVXw0T9g4uUQk+J2NT2mxdsWWhs6B9cqp9d1Z3UjHTtajYHMxGhyUmKYOjQ1GFbbelqzkqOJCldoldDSdhW4rTd21Z5VeK2XKE8UhQMLuWT4JZyQfQIjU0aG7JVukQPZWd1IyZbKYIhds6MmuODe4NQYZuSnUTAkhcLcFEZlJhDh5n7aNTth3QI48avabk/6Lb/1s3z3chZsWsCiLYuobK4kISKBc3LPYVzaOPz48fl9+KwPv/Xjsz58fue+13r3O9Z2v9P51offf4jzA49bbWv7+daP1++c3/F+x+f3+16B+5b9101pkxGbwVlDzuLEQScyI2sGSVFJPfinHbpCOtTKcVr5BLQ2wNS+dU2g2etjZ1VTp97VjuF1V3Aje0dYMLTGMj0vtfPQ4JQYspJiiAzXL/X9jbWWtRVrWbRlEa9teY3ttdvdLqlL+XEaYYNhTNoYrht3HTOyZzA5YzJRHv3yLMeuN42iavb6WL2jJjCM2AmxO6ubAIgKD2NiThI3nDSUglynFzYjIdrlivex4jGwPii8zu1KRHqUtZY15WuYv2k+r2x+hbKGMmLCY5iZM5NZebM4edDJRHq6cO56D2sLvB2Dbls4To5K7lcLPHUVhdr+yu93hh7nTIOsiW5Xc1SaWn3sqGpfObjjEOHSykZ21+4fWrOSnN7VE/PTGZTS1tMaw+CUWDKTot29Ei+9RlsjunDLQhZuXsj2uu14jIdpmdM4K/csDH2rkRmVOorpmdP7/DL/0rPcHEW1u6YpsJCTs5jTR6XVwT26ByXHUDQ0lYIhyRQMSWFMVmLvvmDp8zoLROWfAanD3K5GpEdsrNrI/E3zWbBpAVtrtxIeFs7J2SfzX4X/xczBM/vM/uZtQ58B0GC/LqFQ219tehMqNsDMu9yu5ICqG1rZUlHP5vIGtpa33TawpaKe3TXNnc71hBmyk6MZlBzDySPS91uESaFVDsVay8d7P2bhloUs2rKI0rpSwk0407Omc8vEWzhj8BkKfSK9UKvPz5odNcEAu3xLJaVVjQBEhocxYVAS183IDWyrk8LAxF7WC3s46xdBTSmc92u3KxHpVqV1pSzYtIAFmxbwaeWnhJkwpmZO5aYJN3HmkDM1/FaOiEJtf7XsYWfhibEXu/LtrbWU1TazpbyBLeX1zm1Fe4CtbmztdH5GQhS5abGcPHwAQ1Jjgz2tOamxDEyI6t6VJ6XP8Vs/H+39iIWbnSC7s34n4WHhzMiawW0Tb+OMIWeoERXpZfbUNrfvC7ulklXbq2kO9MJmJUVTMCQlOJR4XHZi6K91UPwIxGfCyFluVyLS5fY27uXVza+yYNMCPtzzIQATB0zkrml3ce7Qc0mPSXe5Qgk1CrX9UdU2WDcfTvpGty484fX52VHVxOby+k6BdWt5A1srGmhs9QXP9YQZBiXHkJsWy4UTsxiaFseQtFhy02IZkhpLbKT+qsrx8Vs/q/as4tXNr7JoyyJ2N+wmIiyCE7NP5CuTv8LMwTMVZEV6mTU7avjr2xtYvrWSbRVOL2yExzAuO4lrpudSmJtCQW4yWUl9bG/Gqq3w2UI49U7Q3s/SR1Q3V/P61teZv2k+y3Ytw2/9jEwZydcLvs6sobPISchxu0Q5Av6WFvzV1fg6fdXgq67CV13tPFdVja+mJvh87t/mEZHVvbsnKCn0RyWPOLdFNxz3WzW1+tha0bBfj+uW8npKKxs77dcaFR7GkFQnqJ48Ij0YWIemxTEoJUZDhKXL+a2fFWUrWLRlEYu2LKKsoYyIsAhOGnQSXy/4OjMHz9Q+qyJdrCsXimr2+li8oZyCISl88YShFOQmMy47qe9vk1byN2fp/YIvul2JyHFpaG3gre1vMX/TfN4tfRev38vghMHcPOFmzs87n/zkfLdL7JestfjrG/AHgqivpsYJosGQ2hZQa/YJr9XYxsaDv3FYGJ7ERDxJSYQlJeFJTiYyN9f5edbNFGr7G28zLH/UGc6UPOSIXlLd2BoMrFsrGti8t63ntYFdNU2dzk2IDmdoWhzjByVxwYTOPa4DE6IJC+tbC+1I7+Pz+1hetjy4avGexj1EhkVy8qCTOafwHE7LOY34yHi3yxTps7pyoajJg5P54Ptn9q+VQH2tzqrHI86B5MFuVyNy1Fp8LbxX+h4LNi3gze1v0uhtJCM2g6tHX835eeczNm1s//o33Y2s14uvthZfVRX+Dj2jwYBa06HntO25wHl4vQd9XxMZiScpCU+yE04jcnKIHjfOCazJSc5zbcE1MSl4LCw+HhPmTieVQm1/s+ZFqN8DU28OHrLWsqe2OdDDun+Pa1VD5/mtAxKiyE2N5cThaQxNi+vU45ocG6EfVNLjvH4vy3cvZ+GWhby25TXKm8qJ8kRxyqBTOGfoOZyacypxEXFulykiR6lftifr5kPdbii60e1KRI6Yz+9j6a6lvLL5FRZtWURtSy3JUcnMHjabWXmzKBxYqD3PD8F6vfiqqvBWVOwztLdmn57Tjseq8dfWHvJ9w+LjgwHUk5xEeGZm++O20BroWfUkJbeH0+gQW1gPhdp+w+vzs7O6icR3/oKJHcKf1mWy+f3i4NDhjvNbwwxkJ8cwNC2O8ydkkZsaS26H8BoXpb824j6v30vx7mIWbl7I61tfp6KpgmhPNKfkBILsoFP7zNL/ItKPFM+FpMEw/Cy3KxE5JGstH+75kAWbFvDq5lcpbyonNjyWM4ecyXl553FC9glEhPXPOeHW58NXVYWvogJveQW+ygq8FRX4yivwVjq3vooKvJWV+MrLnZ7TjvtRduTxtIfQxEQ86WlE5g9zQmjbsX17T5OT8SQkYML7z+/sIf1Je9MG771JdUMrJVsrKN5cyeodNWwpr2d7ZSMj7GYWRJXws9Yv8Njirc781tRYZuSntQ8TTo0lJyW2d+/dJ/1Wq7+VZbuWsXDzQt7Y+gaVzZXEhMdwWs5pnJ17NicPOllBVkRCV/kG2PgmnHE3hPXxecMSkqy1fFr5KQs2LeCVza9QWldKZFgkpw0+jVlDZ3FqzqlEh4deL9/hWJ/P6R3tGFLLy/FVVAbuB0JqhXPrq6o6cEg1xgmfqamEp6YSlZ+PZ9pUwlNSA8dS8KSkOOE0MKw3LC6uf45aOUohHWrd3OC9t7DWsrWigeLNlRRvqaRkSwWf7q4DIDzMMGJgAmOzEzlvQhaXlT6Hb0c0N99+Nz/IyNT8VgkJrf5WPtj5AYu2LOL1ra9T3VxNbHgspw0+jXNyz+GkQScRE97HVj4Vkf6pZB6EhcOUa92uRKSTrTVbmb9pPgs2LWBj9UY8xsMJ2Sfw5clf5ozBZ4TcWhXW73dCanl5IIhW4q0IhNSKcrwVlYGQGjhWVQV+/wHfy5OUhCctDU9qClHDhuGZWkR4aiqe1DQnpKamBkOsJzm5X/We9iT9qYaYVp+f1TtqKN5cQckWJ8juqW0GnEWaCoakcNGkbApzU5k8OJmYyMCV3sYq+P18mHg5WZndu6S2yPFq9bWyeOdiFm5eyH+2/YealhriIuKYOXgm5+Sew4nZJ/bJK8EifYFGUR0jbzOs+DuMOh8SMt2uRoRd9buCe8muLl8NQOHAQu4Zcw9n5Z5FanSqyxW2C4bUQE/pQYf8VlY4gbWy8qAhNSwpKRBKU4nKy8NTWIQnNYXw1LTAbYfAmpyMieifQ6x7G4XaXq66sZXlWysp2VxJ8ZYKVm6roqnV+UeYkxLDycPTKcxNoWhoCiMzEg7e+/rhk9Da0GmBKJHepMXXwuIdi1m4ZSH/2fofaltriY+I5/TBp3PO0HOYkT2DKE/37assIl1Do6iO0ZoXobFCC0SJqyqbKlm0ZRHzN81n+e7lWCxj08by7aJvc+7Qc8mM67kLLv6WFnx79+ItL3eG+paX491bjrd8rxNU2461Dff1+Q74PmGJicGQGjl0KDFTCjqE1FTC05znPCkphKekKKSGKIXaXsRay/bKRoq3OPNhizdX8mlZLdaCJ8wwLjuROdOGUJSbStHQFAYmHmFPlbWw7CHImQrZk7v1M4gcjWZfM++Xvs/CLQt5c9ub1LXWkRCZwOlDTufcoedyQtYJRHoi3S5TRKT7lTwCKXmQd5rblUg/U9dSxxvb3mD+pvks2bEEn/UxLGkYX578Zc7LO4/cxNwu+T5te6P6KjqG030Da3nw2MFW9g2LjcWTnk54WhoRQ4YQM3kynrTA8N6UDiE1NVUhtR9RqHVRq8/P2p01gfmwTpAtaxtKHBXOlNwULpiYRVFuCpMGJx/7qsMb34Ty9fC5B7queJFj1ORt4r3S91i4ZSFvbX+L+tZ6EiMTOTv3bM7OPZsTsk4gwqMGSET6kbJPYMt7cPZPwaU9HqV/afI28U7pOyzYtIC3t79Ns6+Z7Lhsrht3Hefnnc/IlJFHtDiRtRZ/dbXTm7q3/ACBtcK5HwistqnpgO/TNi81PC2NqDGjiUtLd8JpWhrhgQDb9nxYjNbRkP0p1PagmqZWVmytonizE2BXbqsKbqUzKDmGGflpFOWmUJibyqjMBDxdtZDTsocgNg3GXdI17ydylBq9jbxb+i6LNi/ire1v0eBtIDkqmVlDZ3F27tlMy5rWb5f9FxGh5BHwRMLka9yuRPqwVn8rS3YsYcGmBbyx7Q3qW+tJi07jshGXcV7eeUwaMAljjLPSb6C31Lt3b7AXtT2wdhgCXFEBra37f7OwsMDQXieIRubmEp6aRnh6Gp60dOc2NdUJrCkpmEiNypLjo1DbTdqGEjuLOTkhdt1uZyhxmIGx2YlcOXUwRUNTKMpNJTOpmxa9qd7ubOR+0tchXPMRpec0tDbwTuk7LNy8kHdK36HR20hKVArnDzufc3LPoSizSEFWRKSlwVn3YsxFEJfudjX9WlVTFeur1rOhagNbarfQ4msBAr2R+LGBLVosFmstFovf+oPntP3nt36wBB8Hn9v31jkp+N6djll/p8cd33vfOvx0+H4d3id43+8nosVPXcUuwipryWyJ4faoMUz0DCGrJQb/0nK8FX9gU9vw38rKA25HYyIigsN+wwcMIHr0GOd+epqzcFJ6WrBn1ZOUhPFoWyrpOQq1XcTr87N2Z60TYLc4CzvtqnGGWMRHhTNlSDLnjc+iaGgKk49nKPHRKpnn/GDSwhO9Sl1LHRuqN7C+cj1barbQ6j/AVc4Qtqt+F++WvkuTr4nU6FRmD5vNOUPPoXBgIeFh+rEjB2d9PmxTE/6mpuBtx/udbhubsM2db9Nvv43wdAUDCSGr/wlN1Wqne1B1c3UwvK6vWs/Gqo2sr1pPeVN58JyY8BiiPFEYDMaYzrcYMBBmwoKP24bqGoxzPHAuELwfZsI6PT7YMWOcrzAfRDf7iWr2E93kJ7LJF3wc1eQnqslHVLOPqCYfkU3ObUSHx5FN3uCt6ZRR64BlwDJqYmODw3ojcocQU1DQPuy3rUe1bdhvQoL2S5VeS79dHqPatqHEgb1hV2ytoqHFGUqcnRTNtDxnMafC3BRGZyZ23VDio+FtgZK/wchZkDyk57+/0NDawKbqTXxW9Vmw8VxftZ5d9buC50SGRfa5VX3jIuO4ePjFnDv0XAoyCvCE6WptKLN+vxMim5uxjY34m5qxTQe+9Tc1Yhub8Dc3HeY28F7B93SCqj3QMLYjYCIiMDExpFx1pUKthJbiuZA+CnJPdLuSPqe6uTrY9m6o2sCG6g1sqNrA3sa9wXNiw2PJT87nlJxTyE/KJz85n+HJw8mMyzzmAOdvacFfVxf88tXV4a+rx1+/z+O25+sP/Ng2Nh7R9wuLiyMsPp6w+CTC4uPwpMQHHsc7j+PjCYuLJywh3hnum5oa7HHV/FTpKxRqj1BpVWNwb9hlmytZt6sGf2Ao8ZisRC4vzKFwaCpFuSlkJ/eSHxBr/gX1ZTBN2/h0t2ZfsxNeK53wuqFqA59VfcaOuh3O0CGc8JqXlEfhwEKGJw8nPymf4SnDGRQ/KHilVqQrWZ+vw+bybRvKB/bqq6rEX1/fHkQPeNuEbWzEtrQcWwHh4YRFR2NiogmLiiYsJhoTFU1YdDSexETCBmY4j9uOd7yNds472O2+xzTMrffQPrVHYecqKC2GWb8G9YAds5qWms7hNfC1p3FP8JyY8Bjyk/I5Kfsk8pM7h9d922BfXR1Nq9fgr6nuHEDrDx9I/fX1R3ZxzuNxwmZ8fDCUelJTiBwy2AmgHQNpWyiNj8cTH9chsMYTFhuL0eJiIgq1B+L1+flkV20gwDpBdme1M5Q4LtLDlCEpfPWMERQNTWHKkBTie2oo8dGoLoVXvgcZ42DYGW5X02e0+lrZXLM52OPa1ohuq90WnFcTbsIZmjSUCekTuGT4JQxPHs7w5OHkJORo6K0cl+Dm8pWVBwyqnTaXLw/s23eQzeU9SUmExcVhYmIIi4rCxMTgiU/ApA8IBMUowqJjDnp72CAaFaVtFPop7VN7FEoegfBomHSl25WEhNqW2gOG17LGsuA5MeExDEsaxozsGcHgmp+cT1Zc1n7h1VtZSVPJcprXb6B54wZa1m+geeNGvLt27futg0xUVHvgDATNiOzszj2iRxBITVSUhvKKdCH9hh2wraKBZ0u2U7KlkhVbK6kPDCXOSoqmMDeFotwUioamMjozgXBPL78i5m2Bf1wH3ia4fJ62BzgGXr+XrbVbWV+5vtOw4a01W/FaLwAe42FwwmBGpozkvLzzyE/OZ0TyCIYkDtECSHJErLX4a2uDC3N4y8vxtQXV8gp8Fc7KksHbysqDby6flER4SgqetDRnc/mCwg6by6c42yG07d+XnIwJ149/EVc118KqZ2D8ZRCT4nY1vUpbeA3Oea125ryWNXQOr3lJeZyQfYLT8xoYOpwdn90pvFpr8ZaV0bhqyX7h1VdRETzPxMYSlZdH3PRpRA7LJzJvKOGpqe09onFxeOLitEqvSC+l32oCdtc08cc3PmN0ZiKXFuQ4qxIPTWVQbxlKfDQW/gC2L4PL/wYDRrpdTa/m8/sorSvdb87r5urNwcWbDIachByGJw/nzCFnBq/85iXlEelR4ybtnI3l6wO9qB3CaMfe045BtbLywFshAGEJCcFQGjF4MDGTJgW2R9Dm8iJ9wkfPQktdv14gqm3RxI69ruur1rO7YXfwnGhPNHlJeUzPnB4cNpyfnL/f1B3r99NaWkp98Vu0bNzYHmA3bMRfVxc8Lywxkaj8fBLOPIPIYflEDc8natgwwrOyNIxXJIT1ulBrjBkG/ABIstZ+vqe+7+TByaz60TkkRIf4L4ar/gFLH4AZd2hf2g781s/O+p2sr+w8bHhj9Uaafc3B87LjshmeMpyTB50cHDacl5RHTHgIXtyQLuFvacG3Zw/evXsPGlS9FYEe1vLyg86lCouLcwJoaioRWVlEjx9HeEoqnjTnmCc1jfBUp6fVk5JCmHoDRPoua50FogZOgEGFblfT7epb6zsH12qnHe64aGKUJ4phScOYmjm1fdhwktPz2nGxQdvaSsvWrdS9/xotGzcEwutGWjZuxDa3t+eeAelEDcsn6aKLiMwfRlT+cKLyh+FJT9ewX5E+qEdCrTFmLnAhUGatHd/h+CzgD4AHeMha+ytr7UbgJmPMsz1RW5twTxgJvX1Y8eGUrYWXvgZDZsBZP3a7GldYa9ndsLtTcF1fuZ4N1Rto9LavIpgRm8GI5BFMzZwaDK/DkocRFxHnYvXSk2xrqxNUy8rw7tlDa1mZc79sT+DWOe6rrDzg601MTCCMphIxIIPoUaODvaidgmqgRzUsqm+tcC0ix6F0OexaBRf8vk8tENXQ2rDfkOENVRvYWb8zeE6UJ2q/RRPbel47hld/UxMtmzZRt35BpyHDLVu2gNcbPC8iO5vI4fnETZ9O1PB8p/c1fxiepKQe/ewi4q6e6qmdB/wJeLTtgDHGA9wHnA1sB5YZY1601q7poZr6lqYaePoLEBnvzKP1hHiP82FYaylvKg+uNtxx4aa61vZhRmnRaQxPGc6lIy4NznkdljyMxMhEF6uX7mS9Xqf3tKwM756yYEBtDYTUttDqq6jYf3N5j8fZ7iAjwxnyW1hAREYG4QMGBLc/8KSkEp6aQlhsrDsfUERCX/Fcp72eeIXblRwzv/XzaeWnLNu1jOJdxXxS8Qk76ncEn29b8X9KxhQuT748OGw4Jz6nU3j11dbSsmEDtRuKO4XX1u3b239GezxEDh5M5PB8Es48sz28DsvTz2IRAXoo1Fpr3zbGDN3n8DRgfaBnFmPMU8DFgELt0bIWXvgKVGyC616ChEy3K+pyNS01FO8qZtmuZawpX8OG6g1UN1cHn0+OSmZ48nAuGHZBsOd1ePJwkqOT3StaupT1+fBVVLT3qHYIqMHguqcMX3nF/iv+hoURnpbmhNWBA4mZMIHwjAzCMwY4x9qCa2qqtoYRke7VWAUfP+eseByV4HY1R6xjiF22axklu0uoaakBYHDCYCZlTOKy5MuCQ4f3Da/eigqaP1lPzcbFNK/fEBw67C1rX/zJREYSmZdHzITxJF18cSC8DiNy6FBNyRCRQ3JzTu0gYFuHx9uB6caYNOAXwBRjzPestb880IuNMbcCtwIMGTKku2vt3RbfB2tfhLN/BkNPcruaLtHQ2sDysuUs3bmUpbuWsrZiLX7rJ9oTzejU0Zyde3YwuOYn55MWnaY5MiHK+v34qqo6h9NgcN3bfnzv3gOu/OsJhNXwjAFEjxtL+ID2sOrczyA8LVWr/YpI71AyD7yNvX6BKL/181nlZ8EQW7y7uFOIPSv3LIoGFjE1cyqZcc7FdGst3t27aV61geqNb3VabdhXVRV877DYWCLz84mbMYPI4flE5TtfETk5urAoIsek1/2WZ60tB24/gvMeAB4AKCoqsoc5ve/a/B4s+iGMmQ0nftXtao5Zs6+ZD8s+ZOkuJ8R+tOcjvNZLeFg4kwZM4raJtzEtcxoTB0zUisMhwlobCKt7OvSslgWHBbe2zV/ds6fT/Kg2npQUwgc44TRqxIj9elXDMzIIT0/Xqr8iEjo2vAFv/AxGnAtZk3rkW1qvF9vSgm1tdW4D9/0tLdiWVmyrc+tvaaK0cgvryz5h897P2FqxEW9zI+E+GBSexIkxg8iJmkJW1ABibQR2VQu2dQm+lnfY3tJC686dtGzYgL++Pvi9PUlJRA4fTsLZZ7cPGR6eT3hmpi5Ei0iXcjPUlgKDOzzOCRyTI1W7C569AVKGwsV/DqnFJlr9razeu9oJsTuXsqJsBS3+FsJMGOPSxnHduOuYljWNKRlTtPKwi/zNzfhravDV1uGvq8VXW4u/7X5NbeBYHf7aWnx1znO+2hr8VdV49+w54ErAYUlJRGQMIHxABlHT8gK9rBnB3taIjAw8AwZoqJmIOPMtN29un1vZcR584L4NPhd8ovO5B3jtfq8JPnfg19hO3/coX9N2W7kFFv0I4odByrXY+fM7hMu20Nl6wAC63/22INraOZge6Pz9pmMcRn7gq7PKwNfHtERG0hoZiYmMxEREBG/DB2aQdMklncKrJzVV4VVEeoSboXYZMMIYk4cTZq8Crj6aNzDGzAZmDx8+vBvK6+V8rfCPG5zN26/9F0T37oWPfH4f6yrXsXTnUj7Y9QHLdy+nwdsAwKiUUVw5+kqmZ06nYGABCZGhM8eoN/O3tOCvrXUCZ6dQuk84bQulbc8Fwqm/tvag29MEGUNYXJyzp2pCAmEJCUQMyCBs+HCnR7Vjr2rgflh0dM/8AYhIyGtc+SHbbrnF7TK6UCxQBy9959CnhYe3B8bICExEBGERbffbA6UnPqFzuOx0PyJ4PywyEhsezl5vFZsbS9lQv4VP6zdRbRvwhkFKwgBGZIxl9MDxjM2cxICk7MD3iSCsw3sSHq6QKiK9Uk9t6fMkMBNIN8ZsB35krX3YGHMH8CrOlj5zrbWrj+Z9rbUvAS8VFRX1pRbvyLz2Y9j6Plz6IAwc63Y1+7HWsqFqAx/s+oClO5d2mouTl5TH7PzZTM+aTtHAIlKiU1yutvexra346vYJm52CaIde0X3DaV0d/poabEvLYb9PWGysE0gTEwiLT8CTlkpkbm6nY2EJ8U5gjQ/cJiTiSYgnLCGBsLg4bVYvIp105QXn6HFjybn/L23v2/GbdL5ln8fBh/se7xjIjvI1B3ztEbymuQ5e+jrU7YaL/ghpwzFh5gAhNHA/IqJL5pW2zYkt3l0cnBNb3VwNEZCTncPUzPM4K3MqRQOLyIrPOu7vJyLipp5a/XjOQY7PB+b3RA19yup/weI/wdRbes12ANZattVu44NdH7Bs5zI+2PUBFU0VAAyKH8RZuWcxNXMq0zKnkRGb4XK17vM3NtK46iMaSoppWvUR3sqKTiHVNjUd9j1MTEywd9QTH48nKYmInEF44hMISwz0nMYntAfQ+Hg8iYntx+LjtSCHiHS5rrzgHJ6aSsLMmcdflFu8LfD458GzAb7yHAyb2W3fym/9rK9aH9xip3h3MVXNVYDTDp8++HSmZU5TiBWRPqnXLRR1NPrl8OO9n8ELd8CgIjj3v10tZVf9LpbuWsoHOz9g6a6l7KrfBcCAmAHMyJ7B9MzpTM2cSk5Cjqt19gbeykoaV6ygoaSExuISGtesgdZWMIbI/GFEZAwkIjMr0Cua6NzGBwJrQrwTRBMT2sNpfLwWSBIR6c2shZe+Bpvegkvu7/JAe7gQO3PwTKYGemKz47O79HuLiPQ2IR1q+93w4+Y6ePoLEB4JV/zNue1Bexv3UryrODikeGvtVgBSolIoyizi5vE3My1rGkMTh/b7OTetpaU0lJTQULKchpJiWtZvcJ6IiCBm/HjSrr+OmMJCYqdMwZOU5G6xIiLS9d78JXz4JJz+A5h8wAFrR8Vv/Wyo2sDSXUsVYkVE9hHSobZfsdaZk7P3U/jC85DU/b2f1c3VFO8uDu4Vu75qPQDxEfEUDSziqtFXMS1zGiNSRhBm+u+8Suv30/zZehqXl9BQXELD8uV4d+4EICw+npgpU0i6cDaxhQVET5jQIwsltba2sn37dpqOYBizSH8WHR1NTk4OERr5IF1p+WPw1q9hyrVw6p3H9BZtIbZtPmzxrmIqmysBhdhQpbZZ5MgcS9usUBsqlj4IHz8LZ9wD+ad3y7doaG2gZHdJcK/YteVrsViiPdEUDCzgwmEXMj1rOqNTRxMe1n//6vhbWmj6+OPgUOKGFSvw1ziLYIUPGEBMUSGxN91EbGEBUSNHujJvdfv27SQkJDB0qHrNRQ7GWkt5eTnbt28nLy/P7XKkr1j/mnMROv9MuPB/j3i7vcOF2FNzTnVCbGYRg+IHdecnkG6itlnk8I61bQ7pZNJv5tRuWwqvfh9GzoKTv9Vlb9vkbeLDPR8G94r9eO/HeK2XiLAIJg2YxJcmf4lpmdOYmD6RCE//7cXw1dbSuHKl0wsbWNipbWXhyLw8Es89h5iCQmKLConIyekVDVVTU5MaTZHDMMaQlpbGnj173C5F+oqdq+CZ65xdCa74Gxyi7WzbJWDZ7mXBebEKsX2b2maRwzvWtjmkQ22/mFNbv9dpIJMGwefuh+PYPqXV38rqvauDCzutLFtJi78Fj/EwLn0c14+/nmmZ05icMZmY8Jgu/BChpXV3WaehxM3r1jmb13s8RI8dS8rVVxNTWEBsYSHhqalul3tQajRFDk//TqTLVG2Dxy+H6GS4+h8QdeA915ftWsaTnzxJye6S4C4B2XHZCrH9hH7miBzesfw7CelQ2+f5ffDsjdBYATctgpij28/Vb/2srVjL0p1L+WDXByzfvZxGbyMGw6jUUVw1+iqmZ02nIKOA+Mj4bvoQvZu1lpZNm9qHEi9fTuu2bYCzZU7M5Emkf/nLxBYWEDNpEmGxsS5XHFp2797NN7/5TZYsWUJKSgqRkZF85zvf4XOf+1yP1bBkyRK+/vWv09zcTHNzM1deeSU//vGPD3p+cXExjz76KH/84x/3e27o0KEUFxeTnp5+0NcfyTki0sc0VjmBtrURbnwFEvffMsdv/Tz80cP8aeWfSItO4+RBJzM1cypTM6cqxEqPUtssfZFCbW/2n184WwFcfB9kTTyqlza0NvDtt77NO6XvADAsaRgX51/M9KzpFA0sIjk6uRsK7v1saytNn3xCQ3GJ0xtbshxfhXOl3JOaSmxhASlXX01sUSHRo0dr25zjYK3lkksu4brrruOJJ54AYMuWLbz44os9Wsd1113HM888w6RJk/D5fKxbt+6Q5xcVFVFUVNRD1YlIyPO2ODsTlK+HLzznDD3eR3VzNd9753u8U/oO5+edz49m/IjYCF0klZ6ntln6qpBestYYM9sY80B1dbXbpXS9dQvgnf+Bgi/ClC8c1Usrmiq46dWbeG/He/xX4X/xxuVv8MIlL/CDE37AWbln9atA66+vp37xYvbc+ye23HAD66afwObLr6Ds17+m6ZN1xJ96Kpk/+ynD5s9nxHvvknPvvaTdcD0xEyYo0B6nN954g8jISG6//fbgsdzcXL761a+yefNmTjnlFAoKCigoKOD9998H4M033+TCCy8Mnn/HHXcwb948AO666y7Gjh3LxIkT+fa3vw3AP/7xD8aPH8+kSZM49dRTD1hHWVkZWVlOr4nH42HsWOcXzvr6em688UamTZvGlClTeOGFF/aroby8nHPOOYdx48Zx8803Y60Nvu/f//53pk2bxuTJk7ntttvw+Xxd8ccmIqHEWnjxDtj8jnMBethp+53y8d6PueKlK1iycwl3T7+bX53yKwVacY3aZumrQrqnts/Oqa3YCM/fBlmT4LzfHtVLS+tKuW3Rbeyq38X/zfw/Th/SPSsl91beiopOQ4mb1qwBnw+MIWr0aJIvvdQZSlxQSMTADLfL7dNWr15NQUHBAZ/LyMhg0aJFREdH89lnnzFnzhyKi4sP+l7l5eX885//5JNPPsEYQ1VVFQA//elPefXVVxk0aFDw2L6++c1vMmrUKGbOnMmsWbO47rrriI6O5he/+AVnnHEGc+fOpaqqimnTpnHWWWd1eu1PfvITTj75ZH74wx/y8ssv8/DDDwOwdu1ann76ad577z0iIiL48pe/zOOPP84Xv/jFo/+DEpHQ9Z9fwKqn4Yy7YdKVnZ6y1vL0uqf59bJfkxGTwaPnPcr49PEuFSriUNssfVVIh9o+qbURnv6iswXAFY9CxJHvabquYh23v3Y7Lb4WHjznQaZkTOnGQt1nraV1+/b2ocTFJbRs2gSAiYwkZuJE0m6+mdiiQmImT8aTcOBFO/q6n7y0mjU7arr0PcdmJ/Kj2eOO6jVf+cpXePfdd4mMjOS1117jjjvuYOXKlXg8Hj799NNDvjYpKYno6GhuuukmLrzwwuDV2pNOOonrr7+eK664gksvvfSAr/3hD3/INddcw8KFC3niiSd48sknefPNN1m4cCEvvvgiv/vd7wBnVcqtW7d2eu3bb7/N888/D8AFF1xASoozr/3111+npKSEqVOnAtDY2EhGhi6SiPQrJfPg7d86I6pO+XanpxpaG/jx4h+zYNMCTs05lf8++b9Jikpyp07pldQ2q22WrqVQG+CrqqLpk0/cLcJaeO+PsOFTOOOHsG4XsOuIXrquYh33rriXMeExfLPwv8je2Ez9xiXdW68b/H6aN26ioaSYxpLleMvKAAhLTCS2oICkSz9HbGER0ePHERYZ6XKx/du4ceN47rnngo/vu+8+9u7dS1FREf/7v//LwIED+fDDD/H7/URHOxdvwsPD8fv9wde0bVAfHh7O0qVLef3113n22Wf505/+xBtvvMH999/PBx98wMsvv0xhYSElJSV8+9vfZsWKFWRnZzN//nwA8vPz+dKXvsQtt9zCgAEDKC8vx1rLc889x6hRozrVvXv37sN+Nmst1113Hb/85S+P+89JRELQZ4vg39+C4WfBBZ33ot1QtYFvvfktNtds5usFX+fG8TcSZkJ6tpf0IWqbpa8K6VDblfvUNn3yCVuvv+H4i+oSafCf+47qFTHAdwCox/vID9h66NNDXnhmJrFTpzq9sAWFRI0YjjmO7Y76sqO9attVzjjjDL7//e/zl7/8hS996UsANDQ0AFBdXU1OTg5hYWH87W9/C855yc3NZc2aNTQ3N9PY2Mjrr7/OySefTF1dHQ0NDZx//vmcdNJJDBs2DIANGzYwffp0pk+fzoIFC9i2bRuPPPJIpzpefvllzj//fIwxfPbZZ3g8HpKTkzn33HO59957uffeezHGsGLFCqZM6Ty64dRTT+WJJ57g7rvvZsGCBVRWOntInnnmmVx88cV885vfJCMjg4qKCmpra8nNze3WP1MR6QV2rAzsRTsOLp8HnvZfpV7e+DI/WfwTYsJjePDsB5mWNc21MqV3U9ustlm6VkiH2q6cUxs9diy5jz3aBVUdo72fwsvfdlY5PvuncIRXdV/b+jqPr/k7I1JG8LWCrxMfEdfNhbovIjubiEHa/qC3M8bwr3/9i29+85v85je/YcCAAcTFxfHrX/+agoICLrvsMh599FFmzZpFXJzz93bw4MFcccUVjB8/nry8vGBDVltby8UXX0xTUxPWWn7/+98DcOedd/LZZ59hreXMM89k0qRJ+9Xx2GOP8c1vfpPY2FjCw8N5/PHH8Xg83HPPPXzjG99g4sSJ+P1+8vLy+Pe//93ptT/60Y+YM2cO48aN48QTT2TIkCEAjB07lp///Oecc845+P1+IiIiuO+++9RwSp9mjIkD3gJ+bK399+HO75OqtsITV0BsKlz9THAv2hZfC79Z9hueXvc0BRkF/Pa035IRq2GP0vuobZa+ynRcMSxUFRUV2UNNZO/1Girgr6cBFm59C+LSDvsSay33rriXBz96kNMHn85vTv0N0eFHPv9W+ra1a9cyZswYt8sQCQkH+vdijCmx1vaJ/SOMMXOBC4Eya+34DsdnAX8APMBD1tpfHeZ9fgrUAWuOJNSGfNu8r8ZKmDsLanbCTa9ChvN3prSulP96879YXb6aG8bfwNemfI3wsJDuM5BuorZZ5Mgdbdusn7pu8/vh+VugbpezYfsRBFqv38tPF/+Uf67/J58f+Xl+MP0HakBFRORg5gF/AoLDkYwxHuA+4GxgO7DMGPMiTsDdd0LajcAkYA3QP6+eepvh6WuhfANc+89goH1r21t8/93vY63lD6f/gTOGnOFyoSIi/ZOSkNve/i2sfw0u+D0MKjzs6Y3eRu58607e2v4Wt0+6nS9P+jKmwwIVIiIiHVlr3zbGDN3n8DRgvbV2I4Ax5ingYmvtL3F6dTsxxswE4oCxQKMxZr611r/veX2StfDCV5y9aC99CPJOwev3ct/K+3joo4cYkzqG/5n5PwxOGOx2pSIi/ZZCrZs+ew3e/CVMmgNFNx729KqmKu544w5W7VnFPSfcwxWjruiBIkVEpA8aBGzr8Hg7MP1gJ1trfwBgjLke2HuwQGuMuRW4FQjOcwt5b/wMPvoHnPlDmHg5exv38p23v8OyXcv4/MjPc9e0u4jyRLldpYhIvxbSobYrVz/ucVVb4fmbIWOs00t7mN7WnXU7ue212yitLeX3M3/PWblnHfJ8ERGRrmatnXeY5x8AHgBnTm1P1NStih+Bd/4HCq+Hk79F8a5i7nz7Tupa6vjFyb/govyL3K5QRESAkN4DxVr7krX21qSkENvQ3NsMz3wR/D648jGIjD3k6Z9WfsoX5n+BvQ17+evZf1WgFRGR41UKdBwvmxM4dtyMMbONMQ9UV1d3xdu559OF8PK3YMQ52PN+x9zVj3DzwpuJj4jn8QseV6AVEelFQjrUhqwF34UdK+Bz90Na/iFPLd5VzPULrgdg3nnzKMrsE4txioiIu5YBI4wxecaYSOAq4MWueOOQveDc0Y4V8I/rIXMC1Rf9ga+9/S3+t+R/OSv3LJ684ElGpox0u0IREelAobanrXwCSh6Bk74Boy845Kmvb3md2xbdRnpsOo+d/5gaUQkpv/jFLxg3bhwTJ05k8uTJfPDBBwD83//9X3Cjd7dUVVXx5z//+aDPezweJk+ezLhx45g0aRL/8z//g9/v7po4//3f/+3q95fQZYx5ElgMjDLGbDfG3GSt9QJ3AK8Ca4FnrLWr3ayz16jcAk9cCbFprDnvv7ly4Y28W/oud027i9+e+lviI+PdrlDkmKlt7lpqm3sPhdqetOsj+Pc3YegpcMY9hzz1mXXP8K23vsXotNE8OutRsuOze6hIkeO3ePFi/v3vf7N8+XJWrVrFa6+9xuDBzkjHUGg4Y2JiWLlyJatXr2bRokUsWLCAn/zkJz1Y4f7UcMqxstbOsdZmWWsjrLU51tqHA8fnW2tHWmvzrbW/cLvOXqGxEh6/HOtt4h+n3sa1b30dn/Uxb9Y8rhlzjXYbkJCmtrnrqW3uPRRqe0pjlbPHXUwKfH4ueA68Rpe1lj+v/DM/W/IzThl0Cg+d8xDJ0ck9WqrI8dq5cyfp6elERTkrgqanp5Odnc0f//hHduzYwemnn87pp58OwMKFC5kxYwYFBQVcfvnl1NXVAVBSUsJpp51GYWEh5557Ljt37gRg5syZfP3rX2fy5MmMHz+epUuXAlBfX8+NN97ItGnTmDJlCi+88AIAq1evZtq0aUyePJmJEyfy2Wefcdddd7FhwwYmT57MnXfeecjPkpGRwQMPPMCf/vQnrLX4fD7uvPNOpk6dysSJE/nrX/8a/MynnnpqsK533nkHgFdeeYWCggImTZrEmWeeecha582bx6WXXsqsWbMYMWIE3/nOdwC46667aGxsZPLkyVxzzTVd8z9JpJuE7JxabzM8dQ0NVZv5weRz+OnHf2Vq5lSeufAZJg2Y5HZ1IsdNbbPa5j7NWhvyX4WFhbZX8/msfeIqa3+Sau2WxQc9zevz2h+//2M7ft54e8+799hWX2sPFil9yZo1a1z9/rW1tXbSpEl2xIgR9ktf+pJ98803g8/l5ubaPXv2WGut3bNnjz3llFNsXV2dtdbaX/3qV/YnP/mJbWlpsTNmzLBlZWXWWmufeuope8MNN1hrrT3ttNPszTffbK219q233rLjxo2z1lr7ve99zz722GPWWmsrKyvtiBEjbF1dnb3jjjvs3//+d2uttc3NzbahocFu2rQp+LoDiYuL2+9YUlKS3bVrl/3rX/9qf/azn1lrrW1qarKFhYV248aN9ne/+539+c9/bq211uv12pqaGltWVmZzcnLsxo0brbXWlpeXH7LWRx55xObl5dmqqirb2NhohwwZYrdu3XrQmqRrHOjfC1Bse0H7Fspfvb5t7sjns/YfN9iNP0+1lzx1up0wb4L9y8q/WJ/f53Zl0oeobVbbLEfuaNvmkN7SJ2S893+wbj7M+hUMOeGApzR5m/ju29/ljW1vcMuEW/jqlK9qmJN0jQV3OUPfu1LmBDjvVwd9Oj4+npKSEt555x3+85//cOWVV/KrX/2K66+/vtN5S5YsYc2aNZx00kkAtLS0MGPGDNatW8fHH3/M2WefDYDP5yMrKyv4ujlz5gBw6qmnUlNTQ1VVFQsXLuTFF1/kd7/7HQBNTU1s3bqVGTNm8Itf/ILt27dz6aWXMmLEiOP66AsXLmTVqlU8++yzAFRXV/PZZ58xdepUbrzxRlpbW7nkkkuYPHkyb775Jqeeeip5eXkApKamBt/jQLUCnHnmmbQtsDN27Fi2bNkSHB4mIt3k9Z/wyqYF/GjwEKKsl/vPvp8Ts090uyrpy9Q2q22WLhXSoTYk9qnd+Jazcfu4S2H67Qc8pbq5mq+98TVWlK3ge9O+x9Vjru7hIkW6nsfjYebMmcycOZMJEybwt7/9bb+G01rL2WefzZNPPtnp+EcffcS4ceNYvHjxAd973ws+xhistTz33HOMGjWq03Njxoxh+vTpvPzyy5x//vn89a9/ZdiwYUf1WTZu3IjH4yEjIwNrLffeey/nnnvufue9/fbbvPzyy1x//fV861vfIiUl5YDvd7BaP/jgg+CwMHD+DL1e71HVKuK2kGibO2hd+ld+t/ZvPJGRzuT08fz2tN+SGZfpdlki3UJts9rmviqkQ6219iXgpaKiolvcruWAanbAszdC2gi46F44QM/rrvpdfOm1L7GlZgu/Oe03zBo6y4VCpU87xFXb7rJu3TrCwsKCV15XrlxJbm4uAAkJCdTW1pKens4JJ5zAV77yFdavX8/w4cOpr6+ntLSUUaNGsWfPHhYvXsyMGTNobW3l008/Zdy4cQA8/fTTnH766bz77rskJSWRlJTEueeey7333su9996LMYYVK1YwZcoUNm7cyLBhw/ja177G1q1bWbVqFZMmTaK2tvaIPsuePXu4/fbbueOOOzDGcO655/KXv/yFM844g4iICD799FMGDRrE3r17ycnJ4ZZbbqG5uZnly5fzgx/8gC9/+cts2rSJvLw8KioqSE1NPWithxIREUFraysRERHH8X9GpPv1+ra5g52rnuTbK/+HVUkJfHHMF/hG0beICNO/MekBapvVNkuXCulQ26t5W+CZ68DbBFf+HaL23wJgQ9UGbn/tdmpbarn/rPuZljXNhUJFul5dXR1f/epXqaqqIjw8nOHDh/PAAw8AcOuttzJr1iyys7P5z3/+w7x585gzZw7Nzc0A/PznP2fkyJE8++yzfO1rX6O6uhqv18s3vvGNYMMZHR3NlClTaG1tZe7cuQDcc889fOMb32DixIn4/X7y8vL497//zTPPPMNjjz1GREQEmZmZfP/73yc1NZWTTjqJ8ePHc9555/Hb3/62U/1tCz+0trYSHh7Otddey7e+9S0Abr75ZjZv3kxBQQHWWgYMGMC//vUv3nzzTX77298SERFBfHw8jz76KAMGDOCBBx7g0ksvxe/3k5GRwaJFiw5a66HceuutTJw4kYKCAh5//PEu/f8l0h+9++E87ir5Ld7IKH5/8n9zdv5st0sS6VZqm9U292XGmXMb2oqKimxxcbHbZXS24Lvwwf3w+Udg/KX7Pb2ybCVfef0rRHoi+ctZf2F06mgXipS+au3atYwZM8btMrrFzJkz+d3vfkdRUZHbpUgfcaB/L8aYEmut/pIdh17ZNgM+v4+/fPBLHlj3FCN8ht+fP4/crEK3y5J+QG2zyJE72rZZPbXd4aNnnUB7wpcPGGj/s/U/3Pn2nWTGZXL/WfeTk5DjQpEiIiL9S3ljOd9987/4oKyESxpb+cHnniU6c4LbZYmIyHFSqO1qZWvhxa/C4BPg7J/u9/Rznz7HT5f8lLGpY7nvrPtIjU51oUiR0PXmm2+6XYKIHEJvXShqRdkKvv3mt6luKOOnFdV87rKnnNViReS4qW0Wt4W5XUCf0lwLT18LkXFw+TzwtE8at9by1w//yo8X/5gZ2TN4+NyHFWhFRKTPsda+ZK29tW37C7dZa/nb6r9x4ys3EtVUzd937ORz5/wv5GrLHhGRvkI9tV3FWnjhK1CxEb74AiS279vl8/v45dJf8vS6p7ko/yJ+fOKPtbqiiIhIN6ttqeWe9+7h9a2vc1Z0Fj/9ZCkJZ/0Exl/mdmkiItKFFGq7ypI/w5oXnCHHeacEDzf7mvneO99j0ZZF3DD+Br5Z8M399vESERGRrrWuYh3fevNblNaV8u2Mk/niB09gpt4MJ37N7dJERKSLKdR2hS3vw8J7YPSFnRrLmpYavv7G1yneXcx3pn6Ha8de62KRIiIi/cM/P/snv/jgFyRFJjF3zC0UzL8bRp4H5/3mgHvGi4hIaNOc2uNVuxv+cQOkDIVL/hxsLMsayrj+letZuWclvzn1Nwq00u/s2rWLq666ivz8fAoLCzn//PP59NNPj+m9br75ZtasWdOl9VVVVfHnP//5qF9XV1fHl770JfLz8ykoKKCwsJAHH3zwuGqZN28ed9xxBwD3338/jz766DG9z+bNm3niiSeOqxaR42WMmW2MeaC6urrHv3ejt5F73ruHH77/QyZnTOaZqXdT8OrPIGsyfP5hCPP0eE0ivYna5iOntjm0hHSodbPhBMDnhWdvgKZquPIxiHYWxdhUvYlr519LaW0pfz7zz5yXd5479Ym4xFrL5z73OWbOnMmGDRsoKSnhl7/8Jbt37z6m93vooYcYO3Zsl9Z4rA3nzTffTEpKCp999hnLly/nlVdeoaKiYr/zvF7vMdV1++2388UvfvGYXquGU3oDtxaK2lKzhS/M/wL/Wv8vbpt4G38tvIu0Z2+B+Ay4+mlnEUeRfkxts9rmviykQ63rKyy+/hPY8h7M/gMMHAfAqj2r+OKCL9Lka+KRWY8wI3uGO7WJuOg///kPERER3H777cFjkyZN4pRTTsFay5133sn48eOZMGECTz/9NOBsBzBz5kw+//nPM3r0aK655hqstYCzqXtxcTEA8fHxwfd89tlnuf766wHYsGEDJ5xwAhMmTODuu+8OnldXV8eZZ55JQUEBEyZM4IUXXgDgrrvuYsOGDUyePJk777wTgN/+9rdMnTqViRMn8qMf/Wi/z7VhwwaWLl3Kz3/+c8LCnB+fAwYM4Lvf/W7wM5xyyilcdNFFwYb+kksuobCwkHHjxvHAAw8E3+uRRx5h5MiRTJs2jffeey94/Mc//jG/+93vgt9v1qxZFBYWcsopp/DJJ58AcP311/O1r32NE088kWHDhvHss88GP9M777zD5MmT+d///d+j/L8mErpe2/IaV/37KnY37ObPZ/6ZO0ZdjeeJK8HvhS885wRbkX5ObbPa5r5Mc2qP1ZoX4f0/QtFNMOlKAN7e/jbf/v/27j06qvr89/j7mzRykwiKXCoiEERISAImhAAGIhpAQQoKReBYPIiIlqrnFH5ysQu6BKWAqHgDRFAQxHIJRdsqxYIBAYEoUAQUsHCQiyTBpiL35Dl/BMcACZmESSYz+bzWmrUye/Z89/Ps7+z95Ltn7z2fDKdWlVrMuHMGN4bf6OcgRfxj+/btxMXFFfja0qVL2bJlC1u3biUzM5PWrVvToUMHAL744gu+/PJLfvnLX9K+fXs+/fRTbrvtNq+W+cQTT/DEE0/Qr18/pk+f7pleuXJlUlNTCQ8PJzMzk8TERHr06MHEiRPZvn07W7ZsAWDFihXs3r2bjRs3Ymb06NGDtLQ0T2wAX375JbGxsZ6iWZDPP/+c7du306hRIwBmz57Ntddey8mTJ2ndujX33XcfZ86cYezYsaSnp3PNNddw++2306pVq0vaGjJkCNOnT+fmm2/ms88+47HHHuOf//wnAIcPH2bt2rXs2rWLHj160Lt3byZOnMiUKVP44IMPvFpnIoHubO5ZXkh/gXk75hFdK5rnOz5PvUo1Ye6v4D8H8n6NoNbN/g5TpFxQbVZtDmYa1JZE5h5Y9hjcEAddnwNg2Z5ljFs3jqY1m/Lana9Rq0otPwcpkudPG//ErmO7fNpms2ub8VTCUyV679q1a+nXrx+hoaHUqVOHjh07smnTJsLDw0lISKB+/foAtGzZkn379nldONevX8+yZcsA6N+/P8OHDwfyTrcaPXo0aWlphISEcPDgwQJPtVqxYgUrVqzwFLDjx4+ze/fuCwrnxSZMmMCiRYs4evQohw4dAiAhIcFTNAGmTZtGamoqAAcOHGD37t0cOXKE5ORkrr/+egD69u17yTVNx48fZ926dfTp08cz7fTp056/e/bsSUhICJGRkSU+dUwkkH3343cM/2Q4WzK20L9Zf4bHDyfMheZdFnRgQ97vxd+ks6WkfFJtVm0W39KgtrjO/Ah/fgBCw6DP21joVbz5r1m89PlLJNZL5MXbX6RamK7bkYotKirKc9pNcVSqVMnzd2hoaIHXvuT/SaxTp04V2eb8+fPJyMggPT2dsLAwGjZsWOD7zIxRo0bxyCOPFNpWZGQkW7duJTc3l5CQEMaMGcOYMWMuOO2qWrWft//Vq1ezcuVK1q9fT9WqVUlOTvYqZoDc3Fxq1KjhOVp9sfzr6qdTwUQqivWH1jNyzUhOnTvF5A6T6dqoa94LH42BHcug83iI6uXXGEXKG9XmPKrNwUmD2uIwg/efgKM74YGl5F5zA5M2/Yn5O+dzV6O7mNB+AmGhYf6OUuQCJT1qeyU6derE6NGjmTlzJkOGDAFg27ZtZGdnk5SUxIwZMxg4cCDHjh0jLS2NyZMne65JKUqdOnXYuXMnt9xyC6mpqVSvXh2AxMRElixZQt++fVm4cKFn/uzsbGrXrk1YWBirVq1i//79AFSvXp0ffvjBM1+XLl34wx/+wIABA7j66qs5ePAgYWFh1K7987V4TZo0IT4+nqeffppnnnmG0NBQTp06VWjhys7OpmbNmlStWpVdu3axYcMGANq0acMTTzxBVlYW4eHhLFq0iNjY2AveGx4eTqNGjVi0aBF9+vTBzNi2bdsl8+V3cU4i/uCcuwe4p0mTJj5vO9dymbFtBq9veZ2IGhE8n/w8ja9pnPfiZzNg/SuQMATaDvP5skV8SbVZtVl8K6BvFFXmNs2Cfy2C28dwpuFt/E/a/zB/53weiHyAiUkTNaAVOc85R2pqKitXriQiIoKoqChGjRpF3bp16dWrFzExMcTGxtKpUycmTZpE3bp1vWoTYOLEiXTv3p127dpRr149z+svvvgiU6dOJSYmhj179vDTDeQGDBjA5s2biY6OZu7cuTRr1gyA6667jvbt29OiRQtGjBhB586d6d+/P23btiU6OprevXsXWIRmzZpFVlaWp4impKQwadKkAmPu2rUr586do3nz5owcOZLExEQA6tWrx7hx42jbti3t27enefPmBb5//vz5vPnmm8TGxhIVFeW5kUZhYmJiCA0NJTY2VjejEL8prZs4fn/qex5b+RivbXmNbo27Mf/u+T8PaHf9Ff7+FNzSDbpO1G/RihRAtTmPanNwcsHw1Xh8fLz9dPe1UnNgE8y5CyI6cfy+WTz5yf/lsyOf8fu43/NgiwdLd9kixbRz585Cd8aBKDo6muXLl19wPczFTpw4QZUqVXDOsXDhQt59990iC40IFLy9OOfSzSzeTyEFBV/W5q0ZWxn+yXCyTmYxqs0oet/c++fTHb/dDG91hzqRMPADuKqqT5Yp4muqzarN4r3i1madfuyNHzNh0UAIr0fmXc/y6D8eYs/3e3j2tme5J+Ief0cnEtRSUlKIjo6+bNEESE9PZ9iwYZgZNWrUYPbs2WUUoYiUFjNjwa4FTNk8hTpV6zDv7nlEXRf18wzHvoEFv4bqdaDfexrQipQR1WYpbzSoLUpuDix5CH7MZP/9c3lk1TCOnTrGy3e8zG03eHfnNxEpuX/84x9ezZeUlMTWrVtLORoRKUsf7fuIiRsnklw/mfG3jeeaSvlOaf4xC97pnXe/iwFL4Orr/ReoSAWj2izljQa1RVn1LHyzmu13juax9GcBmN1lNi1qtfBzYCIiIsEt5aYUJiZN5K5GdxHi8t0G5OxJWNgPsr+Fge9DLd/flEpERAKHbhR1OV99CGum8GmLbgzav4SqYVWZd/c8DWhFRETKQGhIKN0ad7twQJubA0uHwIGNcN8b0KCN/wIUEZFyQYPawhz7N6QO4f0bmjHsxA4aVG/AO3e/w03hN/k7MhERkYprxR9g53LoMgEif+XvaEREpBzQoLYgZ0/Cnx/g7WpXMfqqE9xaJ445XedQq0otf0cmIiJSrjnn7nHOzczOzvZ94xtehw2vQpuhkPiY79sXEZGAVO4Gtc65as65t51zbzjnBvgjhty//p7JZw4wJbwKXRp24fU7X6f6VdX9EYpIwDpy5Aj3338/ERERxMXFcffdd/P111+XqK3BgwezY8cOn8b3n//8h9dee63Y7zt+/DiPPvooERER3HrrrcTFxfHGG29cUSxvvfUWw4YNA2D69OnMnTu3RO3s27ePBQsWFPr64cOH6d69e4na9lbDhg3JzMwEoF27dgXGtXnzZh5//PFSjSP/On3llVd0x80yVFq/U8uO5fDhKGjWHbo8q9+iFSkB1WbvqTb7XmnW5jIZ1DrnZjvnjjrntl80vatz7ivn3B7n3Mjzk+8FFpvZw0CPsogvv7ObZjP64IfMvSac/s36M6nDJK4KvaqswxAJaGZGr169SE5OZu/evaSnp/Pcc8/x3Xfflai9WbNmERkZ6dMYS1o4Bw8eTM2aNdm9ezeff/45H374IceOHbtkvnPnzpUorqFDh/Kb3/ymRO8tqnBOnTqVhx9+uERtl8S6deuAS+OKj49n2rRpZRbHoEGDePnll8tseVIKDmyEpQ9D/Xi49w0ICfV3RCIBR7VZtRmCtzaX1Te1bwFd809wzoUCrwJ3AZFAP+dcJFAfOHB+tpwyig+AH//fOoalT+SvV1fjiZa/Y2TCyAtvTiEiXlm1ahVhYWEMHTrUMy02NpakpCTMjBEjRtCiRQuio6N57733AFi9ejXJycn07t2bZs2aMWDAAMwMgOTkZDZv3gzA1Vdf7Wlz8eLFPPjggwDs3buXxMREoqOjefrppz3zHT9+nDvuuINbb72V6Ohoz4++jxw5kr1799KyZUtGjBgBwOTJk2ndujUxMTGMHTv2krz27t3Lxo0bGT9+PCEhefuG66+/nqeeesqTQ1JSEj169PAU+p49exIXF0dUVBQzZ870tDVnzhyaNm1KQkICn376qWf6uHHjmDJlimd5Xbt2JS4ujqSkJHbt2gXAgw8+yOOPP067du1o3Lgxixcv9uS0Zs0aWrZsyQsvvHBJ/EuWLKFr17xdcU5ODsOHD6dFixbExMR4CsvHH39Mq1atiI6OZtCgQZw+fRrIO8o7duxYz3r8KZasrCw6d+5MVFQUgwcP9vRZ/r66OK7Vq1d7jkofO3aMnj17EhMTQ2JiItu2bfOsh0GDBpGcnEzjxo0vKLTvvPMOCQkJtGzZkkceeYScnJzLrtOqVavSsGFDNm7ceMk6kQCQtRcW9IXq9aDfQv0WrUgJqTarNufvq6CrzWZWJg+gIbA93/O2wEf5no86/3gA6H5+2kJv2o6Li7MrlZm1x349q4XFzomy1O3zrrg9EX/asWOHX5f/0ksv2ZNPPlnga4sXL7Y777zTzp07Z0eOHLEbb7zRDh06ZKtWrbLw8HA7cOCA5eTkWGJioq1Zs8bMzDp27GibNm0yM7Nq1ap52lq0aJENHDjQzMy6detmCxYsMDOz119/3TPf2bNnLTs728zMMjIyLCIiwnJzc+3f//63RUVFedr66KOP7OGHH7bc3FzLycmxbt262SeffHJB7H/5y1+sZ8+ehea9atUqq1q1qn3zzTeeaVlZWWZmduLECYuKirLMzEw7dOiQ3XjjjXb06FE7ffq0tWvXzn7729+amdnYsWNt8uTJZmbWqVMn+/rrr83MbMOGDXb77bebmdnAgQOtd+/elpOTY19++aVFRER4lt+tW7cCY/vmm2/s1ltv9Tx/7bXX7L777rOzZ8964jx58qTVr1/fvvrqKzMze+CBB+yFF14wM7ObbrrJpk2bZmZmr776qj300ENmZva73/3O/vjHP5qZ2QcffGCAZWRkmNnPfXVxXPmfDxs2zMaNG2dmZh9//LHFxsZ61kPbtm3t1KlTlpGRYddee62dOXPGduzYYd27d7czZ86Ymdmjjz5qb7/99mXXqZnZ+PHjbcqUKQWum4K2F2CzlVF9DNaHL2qzHc8we6ml2Z8amWXuufL2RPxItVm1+WKqzb6rzf78ndob+PkbWYBvgTbANOAV51w34P3C3uycGwIMAWjQoMEVB/Pv/as4EGJMa/l/6BD1v664PZHy4sizz3J65y6ftlmpeTPqjh5doveuXbuWfv36ERoaSp06dejYsSObNm0iPDychIQE6tevD0DLli3Zt28ft912m1ftrl+/nmXLlgHQv39/hg8fDuQduBs9ejRpaWmEhIRw8ODBAk+1WrFiBStWrKBVq1ZA3lHk3bt306FDh0KXOWHCBBYtWsTRo0c5dOgQAAkJCTRq1Mgzz7Rp00hNTQXgwIED7N69myNHjpCcnMz1118PQN++fS+5pun48eOsW7eOPn36eKb9dGQW8o4yh4SEEBkZ6dWpY4cPH/YsD2DlypUMHTqUX/wirwxce+21bN26lUaNGtG0aVMABg4cyKuvvsqTTz4JwL333gtAXFwcS5cuBSAtLc3zd7du3ahZs2aRseS3du1alixZAkCnTp3Iysriv//9r6e9SpUqUalSJWrXrs13333Hxx9/THp6Oq1btwbg5MmT1K5dm88+++yy67R27dqeI9gSQI5sgxNZMGAxXBfh72hEfEa1WbUZVJt9WZv9OagtkJn9CPxvL+abCcwEiI+PtyJmL1J8q8F8FHE31cN/eaVNiVR4UVFRntNuiqNSpUqev0NDQwu89sXluznMqVOnimxz/vz5ZGRkkJ6eTlhYGA0bNizwfWbGqFGjeOSRRwptKzIykq1bt5Kbm0tISAhjxoxhzJgxF5x2Va1aNc/fq1evZuXKlaxfv56qVauSnJzsVcwAubm51KhRgy1bthT4ev51lXfw8vKqVKni9bIL89MyC+sbXyvo82BmDBw4kOeee+6CeX/6p6kwp06dokqVKqURppSmiE7w5L+gso9vOiVSAak251FtvjLltTb784LRg8CN+Z7XPz/NbzSglWBUd/Robpo316ePoo4Ed+rUidOnT19wncq2bdtYs2YNSUlJvPfee+Tk5JCRkUFaWhoJCQle51OnTh127txJbm6u5ygrQGJioueo4sKFCz3Ts7OzqV27NmFhYaxatYr9+/cDUL16dX744QfPfF26dGH27NkcP34cgIMHD3L06NELlt2kSRPi4+N5+umnPdeKnDp1qtDClZ2dTc2aNalatSq7du1iw4YNALRp04ZPPvmErKwszp49y6JFiy55b3h4OI0aNfK8ZmZs3br1suvm4pzya9q0Kfv27fM8T0lJYcaMGZ4CeOzYMW655Rb27dvHnj17AJg3bx4dO3a87DI7dOjgudHE3//+d77//vtixZWUlMT8+fOBvH80atWqRXh4eKHLu+OOO1i8eLGnb44dO8b+/fuLXKdff/01LVq0uGwuUk5pQCtBSLVZtRlUm31Zm/05qN0E3Oyca+Scuwq4H1henAZK9bfwRKTEnHOkpqaycuVKIiIiiIqKYtSoUdStW5devXoRExNDbGwsnTp1YtKkSdStW9erNgEmTpxI9+7dadeuHfXq1fO8/uKLLzJ16lRiYmLYs2cPP/2cyIABA9i8eTPR0dHMnTuXZs2aAXDdddfRvn17WrRowYgRI+jcuTP9+/enbdu2REdH07t37wJ39rNmzSIrK8tTRFNSUpg0aVKBMXft2pVz587RvHlzRo4cSWJiIgD16tVj3LhxtG3blvbt29O8efMC3z9//nzefPNNYmNjiYqK8txIozAxMTGEhoYSGxt7yc0oqlWrRkREhKcoDh48mAYNGnj6YsGCBVSuXJk5c+bQp08foqOjCQkJueCGIgUZO3YsaWlpREVFsXTp0gIvB7lcXOPGjSM9PZ2YmBhGjhzJ22+/fdnlRUZGMn78eDp37kxMTAwpKSkcPny4yHX66aefkpKSctm2RUSCmWpzHtXmouMKxNrsvPlq/IoX4ty7QDJQC/gOGGtmbzrn7gZeBEKB2WY2oSTtx8fH2093XxMR2LlzZ6E740AUHR3N8uXLL7ge5mInTpygSpUqOOdYuHAh7777bpGFpqJJTU0lPT2d8ePH+zuUMvXFF18wdepU5s2bV+DrBW0vzrl0M4svi/iClWqzyIVUm1WbC6La7JvaXCbX1JpZv0Km/w34W1nEICKBKSUlhejo6MsWTYD09HSGDRuGmVGjRg2f/qB3sOjVqxdZWVn+DqPMZWZm8swzz/g7DBGRoKHa7Duqzb5RJt/Ulhbn3D3APU2aNHl49+7d/g5HpNwItqPBIqVJ39T6lmqzSMFUm0W8V9za7M9raq+Ymb1vZkN+Oj9fRERE/Eu1WUREylpAD2pFpHCBfBaGSFnRdiIiZUn7HJGilWQ70aBWJAhVrlyZrKwsFU+RyzAzsrKyqFy5sr9DEZEKQLVZpGglrc1lcqOo0pLvuh1/hyJSrtSvX59vv/2WjIwMf4ciUq5VrlyZ+vXr+zsMEakAVJtFvFOS2hzQg1ozex94Pz4+/mF/xyJSnoSFhRV5R0IREREpO6rNIqVHpx+LiIiIiIhIwNKgVkRERERERAJWQA9qnXP3OOdmZmdn+zsUERERERER8QMXDHdgc85lAPvPP70GKOkotxaQWUA73vxd2LIvF483sV5JPlBwTsXJoTznVNJ+KU5O3sZaljl5k09Jphd3nsvR5y74cvI2fn9uS77ch99kZteXMA5BtbkI2keqNpckVn3uLqXarNqcx8yC6gHMvIL3bi6oHW/+LmzZl4vHm1ivJJ/CcipODuU5p5L2S3Fy8jbWsszJm3yuNCd97vz7uSuPOXkbvz+3pdLYh+vhm0d569fysu3lb6si7U+uNKey2J8UNydv9ydXkpM+d/793JXHnLyN35/bUmnswwt6BPTpx4V4vxTa8ebvwpZ9uXi8idVX+eRvqzg5lOecStovxcnJ21jL2+euJNOLO4+39LkrfPmBlJO38QfjtiRXrrz1a3nZ9vK3VZH2J948L2za5aaXdL7itKPaXPQ0b14rzjzeUG0u/rTLTS/pfMVpx6e1OShOP/YV59xmM4v3dxy+pJwCQ7DlFGz5gHIKBMGWj+QJxn5VToEh2HIKtnxAOQWCssonGL+pvRIz/R1AKVBOgSHYcgq2fEA5BYJgy0fyBGO/KqfAEGw5BVs+oJwCQZnko29qRUREREREJGDpm1oREREREREJWBrUioiIiIiISMDSoFZEREREREQClga1XnLONXbOvemcW+zvWHzFOdfTOfeGc+4951xnf8fjC8655s656c65xc65R/0djy8456o55zY757r7OxZfcM4lO+fWnO+nZH/H4wvOuRDn3ATn3MvOuYH+judKOeeSzvfPLOfcOn/H4wvOuQbOuWXOudnOuZH+jkd8Q7U5MKg2l3+qzeWfavPlVYhB7fkVddQ5t/2i6V2dc1855/YUtSLN7Bsze6h0I/Wej3JaZmYPA0OBvqUZrzd8lNNOMxsK/BpoX5rxFsUX+Zz3FPDn0omyeHyUkwHHgcrAt6UVq7d8lNOvgPrAWfyck4+2ozXnt6MPgLdLM15v+KiPooHFZjYIaFVqwYrXVJsLptpculSbC6XaXIpUmwvls9pcIe5+7JzrQN6GOtfMWpyfFgp8DaSQ90HfBPQDQoHnLmpikJkdPf++xWbWu6xiL4yPc3oemG9mn5dR+AXyVU7OuR7Ao8A8M1tQVvFfzBf5ALHAdeQVmUwz+6Bsoi+Yj3LKNLNc51wdYKqZDSir+Avio5wGAd+b2Qx/7yN8vG/4M/CQmf1QRuEXyEd9lAMsJu8ft3lmNqdsopfCqDarNvuDarNqsz+oNpd+bf5FSd8YSMwszTnX8KLJCcAeM/sGwDm3EPiVmT0HlPtTSXyRk3POAROBv/u7aILv+snMlgPLnXN/BfxWOH3UR8lANSASOOmc+5uZ5ZZm3Jfj423pe6BSqQRaDD7qp2+BM+ef5pRiuEXyVR855xoA2f4umuCzPhoOjD3f1mJAg1o/U21WbfYH1WbVZn9QbS792lwhBrWFuAE4kO/5t0CbwmZ2zl0HTABaOedGne+c8qZYOQG/A+4ErnHONTGz6aUZXAkVt5+SgXvJ2yH/rTQDK6Fi5WNmYwCccw9y/ihqqUZXMsXto3uBLkAN4JVSjazkirstLQVeds4lAWmlGVgJFTcfgIco3wO/4ub0ITDOOdcf2FeKccmVUW1WbfYH1WbVZn9QbfZhba7Ig9piMbMs8q5vCRpmNg2Y5u84fMnMVgOr/RyGz5nZW/6OwVfMbCl5hSZomNkJ8gpN0DCzsf6OwZfMbDvg99NTxbdUmwODanP5p9ocGFSbC1chbhRViIPAjfme1z8/LZApp/Iv2PIB5RQIgi0fCM6cJDj7VTmVf8GWDyinQBBs+YAfc6rIg9pNwM3OuUbOuauA+4Hlfo7pSimn8i/Y8gHlFAiCLR8IzpwkOPtVOZV/wZYPKKdAEGz5gD9zMrOgfwDvAof5+ZbeD52ffjd5d+jaC4zxd5zKKbhyCrZ8lJP/Y62I+QRrTnoEZ78qp/L/CLZ8lJP/Y62I+ZTHnCrET/qIiIiIiIhIcKrIpx+LiIiIiIhIgNOgVkRERERERAKWBrUiIiIiIiISsDSoFRERERERkYClQa2IiIiIiIgELA1qRUREREREJGBpUCsiIiIiIiIBS4NaERERERERCVga1IpUAM65rs65LecfnznntO2LiIiISFBwZubvGESklDnndgMdzOywv2MREREREfElfVsjUjH8DdjmnHvR34GIiIiIiPjSL/wdgIiULudcO8AB9czsnL/jERERERHxJX1TKxL8+gBfm9k5lyfc3wGJiIiIiPiKrqkVCXLOuQTgTcCAk8BjZpbu36hERERERHxDg1oREREREREJWDr9WERERERERAKWBrUiIiIiIiISsDSoFRERERERkYClQa2IiIiIiIgELA1qRUREREREJGBpUCsiIiIiIiIBS4NaERERERERCVga1IqIiIiIiEjA+v+/pOZl9fXihAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1152x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def CG_cond(A, b, eps, k_max = 1000000):\n",
    "    \"\"\"\n",
    "    Return the Conjugate Gradient algorithm estimate solution x to the problem\n",
    "    Ax = b, after diagonal conditioning A, and the number of iterations k it\n",
    "    took to decrease maximum norm error below eps or to exceed iteration maximum\n",
    "    k_max.\n",
    "    \"\"\"\n",
    "    \n",
    "    D = np.diag(np.diag(A))\n",
    "    C = np.sqrt(linalg.inv(D))\n",
    "    A_tilde = C @ A @ C\n",
    "    b_tilde = C @ b\n",
    "    \n",
    "    x_tilde, k = CG(A_tilde, b_tilde, eps, k_max)\n",
    "    x = C @ x_tilde\n",
    "    return x, k\n",
    "\n",
    "# Sorry for copying.\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, sharex = True, figsize = (16,6))\n",
    "\n",
    "for alg, alg_name in [(GS, \"Gauss-Seidel\"), (SD, \"Steepest Descent\"), (CG, \"Conjugate Gradient\"), (CG_cond, \"Conjugate Gradient (conditioned)\")]:\n",
    "    k_list = []\n",
    "    t_list = []\n",
    "    for eps in eps_list:\n",
    "        t_start = time.time()\n",
    "        x, k = alg(A, b, eps)\n",
    "        t_list.append(time.time() - t_start)\n",
    "        k_list.append(k)\n",
    "    ax[0].plot(eps_list, k_list, label=alg_name)\n",
    "    ax[1].plot(eps_list, t_list, label=alg_name)\n",
    "    \n",
    "ax[0].set_xscale(\"log\")\n",
    "ax[0].invert_xaxis()\n",
    "ax[0].set_yscale(\"log\")\n",
    "ax[0].set_xlabel(\"$\\epsilon$\")\n",
    "ax[0].set_ylabel(\"$k$\")\n",
    "ax[0].legend()\n",
    "\n",
    "ax[1].set_yscale(\"log\")\n",
    "ax[0].set_xlabel(\"$\\epsilon$\")\n",
    "ax[1].set_ylabel(\"runtime (seconds)\")\n",
    "ax[1].legend()\n",
    "\n",
    "fig.show()\n",
    "\n",
    "print(\"The number of iterations is brought down by a lot due to the conditioning, bringing down the runtime.\")"
   ]
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}