after_lectures/12/cg_convergence.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.sparse.linalg as spsplin\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convergence of CG for different spectrum\n",
    "n = 1000\n",
    "A = np.random.randn(n, n)\n",
    "# A = A.T @ A\n",
    "Q, _ = np.linalg.qr(A)\n",
    "A = Q @ np.diag([1, 10] + [100 + i for i in range(n-2)]) @ Q.T\n",
    "eigvals = np.linalg.eigvalsh(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000000000000953 1097.0000000000016\n"
     ]
    }
   ],
   "source": [
    "print(min(eigvals), max(eigvals))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fc4f94b6190>"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAPgklEQVR4nO3cf6zddX3H8edr/aEFZ1rWO4O3ja2mYTboLLshRRZDxITCjCX8M8gYSnSNCTp/LBrQP8z+WHTTGDAaSAfVMRlsQWSNYcPFHyEzitzSCsXSWUDtbau9hgBGyPjhe3+cr9vxetp7bjn13vvZ85Gc9Hw/n+855/Npy7Pnfu+5pKqQJLXrd+Z7AZKkk8vQS1LjDL0kNc7QS1LjDL0kNW7pfC9gkNWrV9e6devmexmStGjs2rXrZ1U1NmhuQYZ+3bp1TE5OzvcyJGnRSPKjY8156UaSGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGjdr6JPsSHI0yd5jzCfJZ5IcSPJAkrNmzC9JsjvJV0a1aEnS8IZ5R/8FYMtx5i8ENnS3bcD1M+bfB+w7kcVJkl68WUNfVfcAjx/nlK3AzdXzHWBlktMBkqwB/gS4cRSLlSTN3Siu0Y8DB/uOp7oxgGuBDwO/nO1JkmxLMplkcnp6egTLkiTBaEKfAWOV5K3A0araNcyTVNX2qpqoqomxsbERLEuSBKMJ/RSwtu94DXAYOBd4W5IfArcBb07yxRG8niRpDkYR+p3AFd2nbzYDT1bVkaq6pqrWVNU64FLg61V1+QheT5I0B0tnOyHJrcB5wOokU8DHgGUAVXUDcBdwEXAAeBq48mQtVpI0d7OGvqoum2W+gKtmOeebwDfnsjBJ0mj4k7GS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNmzX0SXYkOZpk7zHmk+QzSQ4keSDJWd342iTfSLIvyUNJ3jfqxUuSZjfMO/ovAFuOM38hsKG7bQOu78afB/6qql4LbAauSrLxxJcqSToRs4a+qu4BHj/OKVuBm6vnO8DKJKdX1ZGqur97jp8D+4DxUSxakjS8UVyjHwcO9h1PMSPoSdYBm4B7R/B6kqQ5GEXoM2Cs/ncyeRnwJeD9VfXUMZ8k2ZZkMsnk9PT0CJYlSYLRhH4KWNt3vAY4DJBkGb3I31JVdxzvSapqe1VNVNXE2NjYCJYlSYLRhH4ncEX36ZvNwJNVdSRJgJuAfVX16RG8jiTpBCyd7YQktwLnAauTTAEfA5YBVNUNwF3ARcAB4Gngyu6h5wJ/DjyYZE839pGqumuE65ckzWLW0FfVZbPMF3DVgPH/ZPD1e0nSb5E/GStJjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjVs62wlJdgBvBY5W1ZkD5gNcB1wEPA28o6ru7+a2dHNLgBur6hMjXPuvuXP3IT55934OPfEMSxJeqGJ85Qo+dMEZXLxpfOSv9cF/3sMvR/qsktRz7mtO45a/OGdkzzfMO/ovAFuOM38hsKG7bQOuB0iyBPhcN78RuCzJxhez2GO5c/chrrnjQQ498QwAL1QBcOiJZ7jmjge5c/ehkb7W+428pJPoW488zp/9/bdH9nyzhr6q7gEeP84pW4Gbq+c7wMokpwNnAweq6tGqeha4rTt35D55936eee6FgXPPPPcCn7x7/0hfS5JOtm89crzszs0ortGPAwf7jqe6sWOND5RkW5LJJJPT09NzWsDh7p38ic6P8rUkaaEZRegzYKyOMz5QVW2vqomqmhgbG5vTAl65csWLmh/la0nSQjOK0E8Ba/uO1wCHjzM+ch+64AxWLFsycG7FsiV86IIzRvpaknSynfua00b2XKMI/U7givRsBp6sqiPAfcCGJOuTLAcu7c4duYs3jfPxS17HePdue0l6X0yMr1zBxy953Ug/dXPxpnGu/dM3+LlUSSfNqD91k6pjXk3pnZDcCpwHrAZ+CnwMWAZQVTd0H6/8LL1P5jwNXFlVk91jLwKupffxyh1V9TfDLGpiYqImJydPYDuS9P9Tkl1VNTFobtbP0VfVZbPMF3DVMebuAu4aZpGSpJPDKxCS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNM/SS1DhDL0mNGyr0SbYk2Z/kQJKrB8yvSvLlJA8k+W6SM/vmPpDkoSR7k9ya5KWj3IAk6fhmDX2SJcDngAuBjcBlSTbOOO0jwJ6qej1wBXBd99hx4C+Biao6E1gCXDq65UuSZjPMO/qzgQNV9WhVPQvcBmydcc5G4GsAVfUwsC7JK7q5pcCKJEuBU4DDI1m5JGkow4R+HDjYdzzVjfX7HnAJQJKzgVcBa6rqEPAp4MfAEeDJqvrqi120JGl4w4Q+A8ZqxvEngFVJ9gDvBXYDzydZRe/d/3rglcCpSS4f+CLJtiSTSSanp6eHXb8kaRbDhH4KWNt3vIYZl1+q6qmqurKq3kDvGv0Y8BjwFuCxqpququeAO4A3DnqRqtpeVRNVNTE2Njb3nUiSBhom9PcBG5KsT7Kc3jdTd/afkGRlNwfwLuCeqnqK3iWbzUlOSRLgfGDf6JYvSZrN0tlOqKrnk7wHuJvep2Z2VNVDSd7dzd8AvBa4OckLwPeBd3Zz9ya5HbgfeJ7eJZ3tJ2UnkqSBUjXzcvv8m5iYqMnJyflehiQtGkl2VdXEoDl/MlaSGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGmfoJalxhl6SGjdU6JNsSbI/yYEkVw+YX5Xky0keSPLdJGf2za1McnuSh5PsS3LOKDcgSTq+WUOfZAnwOeBCYCNwWZKNM077CLCnql4PXAFc1zd3HfDvVfUHwB8C+0axcEnScIZ5R382cKCqHq2qZ4HbgK0zztkIfA2gqh4G1iV5RZKXA28Cburmnq2qJ0a1eEnS7IYJ/ThwsO94qhvr9z3gEoAkZwOvAtYArwamgc8n2Z3kxiSnDnqRJNuSTCaZnJ6enuM2JEnHMkzoM2CsZhx/AliVZA/wXmA38DywFDgLuL6qNgG/AH7jGj9AVW2vqomqmhgbGxty+ZKk2Swd4pwpYG3f8RrgcP8JVfUUcCVAkgCPdbdTgKmqurc79XaOEXpJ0skxzDv6+4ANSdYnWQ5cCuzsP6H7ZM3y7vBdwD1V9VRV/QQ4mOSMbu584PsjWrskaQizvqOvqueTvAe4G1gC7Kiqh5K8u5u/AXgtcHOSF+iF/J19T/Fe4JbuH4JH6d75S5J+O1I183L7/JuYmKjJycn5XoYkLRpJdlXVxKA5fzJWkhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcYZekhpn6CWpcamq+V7Db0gyDfzoBB++GvjZCJezkLS8N2h7fy3vDdre32LZ26uqamzQxIIM/YuRZLKqJuZ7HSdDy3uDtvfX8t6g7f21sDcv3UhS4wy9JDWuxdBvn+8FnEQt7w3a3l/Le4O297fo99bcNXpJ0q9r8R29JKmPoZekxjUT+iRbkuxPciDJ1fO9nrlKsjbJN5LsS/JQkvd146cl+Y8kP+h+XdX3mGu6/e5PcsH8rX44SZYk2Z3kK91xS3tbmeT2JA93f4bnNLa/D3R/L/cmuTXJSxfz/pLsSHI0yd6+sTnvJ8kfJXmwm/tMkvy29zKUqlr0N2AJ8AjwamA58D1g43yva457OB04q7v/u8B/ARuBvwOu7savBv62u7+x2+dLgPXd/pfM9z5m2eMHgX8CvtIdt7S3fwDe1d1fDqxsZX/AOPAYsKI7/hfgHYt5f8CbgLOAvX1jc94P8F3gHCDAvwEXzvfeBt1aeUd/NnCgqh6tqmeB24Ct87ymOamqI1V1f3f/58A+ev+BbaUXEbpfL+7ubwVuq6r/rqrHgAP0fh8WpCRrgD8BbuwbbmVvL6cXjpsAqurZqnqCRvbXWQqsSLIUOAU4zCLeX1XdAzw+Y3hO+0lyOvDyqvp29ap/c99jFpRWQj8OHOw7nurGFqUk64BNwL3AK6rqCPT+MQB+vzttse35WuDDwC/7xlrZ26uBaeDz3aWpG5OcSiP7q6pDwKeAHwNHgCer6qs0sr8+c93PeHd/5viC00roB10XW5SfG03yMuBLwPur6qnjnTpgbEHuOclbgaNVtWvYhwwYW5B76yyldxng+qraBPyC3pf+x7Ko9tddq95K77LFK4FTk1x+vIcMGFuw+xvCsfazaPbZSuingLV9x2vofWm5qCRZRi/yt1TVHd3wT7svEel+PdqNL6Y9nwu8LckP6V1We3OSL9LG3qC33qmqurc7vp1e+FvZ31uAx6pquqqeA+4A3kg7+/uVue5nqrs/c3zBaSX09wEbkqxPshy4FNg5z2uak+679TcB+6rq031TO4G3d/ffDvxr3/ilSV6SZD2wgd43hhacqrqmqtZU1Tp6fzZfr6rLaWBvAFX1E+BgkjO6ofOB79PI/uhdstmc5JTu7+n59L6H1Mr+fmVO++ku7/w8yebu9+WKvscsLPP93eBR3YCL6H1S5RHgo/O9nhNY/x/T+7LvAWBPd7sI+D3ga8APul9P63vMR7v97meBfrd/wD7P4/8+ddPM3oA3AJPdn9+dwKrG9vfXwMPAXuAf6X0CZdHuD7iV3vcbnqP3zvydJ7IfYKL7PXkE+Czd/21god38XyBIUuNauXQjSToGQy9JjTP0ktQ4Qy9JjTP0ktQ4Qy9JjTP0ktS4/wFTI4sC6VJ5PgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(eigvals, np.ones(n))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1097.0000000000516\n",
      "0.9981785063752278\n",
      "0.9413850541544568\n"
     ]
    }
   ],
   "source": [
    "q = np.linalg.cond(A)\n",
    "print(q)\n",
    "print((q - 1) / (q + 1))\n",
    "print((np.sqrt(q) - 1) / (np.sqrt(q) + 1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Richardson iteration convergence \n",
    "\n",
    "$$ \\|x_{k+1} - x^* \\|_2 \\leq C \\left(\\frac{\\kappa - 1}{\\kappa+1}\\right)^k $$\n",
    "\n",
    "- Chebyshev iteration convergence\n",
    "\n",
    "$$ \\|x_{k+1} - x^* \\|_2 \\leq C \\left(\\frac{\\sqrt{\\kappa} - 1}{\\sqrt{\\kappa}+1}\\right)^k $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "rhs = np.random.randn(n)\n",
    "x0 = np.random.randn(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [],
   "source": [
    "conv_cg = []\n",
    "\n",
    "def cl_cg(x):\n",
    "    conv_cg.append(x.copy())\n",
    "    \n",
    "x_cg, info = spsplin.cg(A, rhs, callback=cl_cg, x0=x0, tol=1e-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.504019900942985e-05"
      ]
     },
     "execution_count": 77,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(A @ x_cg - rhs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [],
   "source": [
    "def richardson(A, b, tau, x0, tol, max_iter):\n",
    "    x = x0.copy()\n",
    "    conv = []\n",
    "    for i in range(max_iter):\n",
    "        x = x - tau * (A @ x - b)\n",
    "        conv.append(x)\n",
    "        if np.linalg.norm(A @ x - b) < tol:\n",
    "            break\n",
    "    return x, conv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "35.52946857894798\n"
     ]
    }
   ],
   "source": [
    "x_rich, conv_rich = richardson(A, rhs, 2 / (min(eigvals) + max(eigvals)), x0,\n",
    "                               1e-6, 1000)\n",
    "print(np.linalg.norm(A @ x_rich - rhs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$||Ax_k - b||_2$')"
      ]
     },
     "execution_count": 80,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEUCAYAAADwYOuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+QklEQVR4nO3deXxU1d348c83k2SSTEJIZN8REbe6AIKA1aj0UaoW11pbLPhUrVpbfdQK1LagtWqtVUHlsf5UqPVR61ZURK0Lcd/AXUAEBAQFZCfBkJB8f3+cO8kkmck6S2bm+3697mvmnrudc2cy35x7zz1HVBVjjDEmHjISnQFjjDHpw4KOMcaYuLGgY4wxJm4s6BhjjIkbCzrGGGPiJjPRGejounTpogMGDGjTtuXl5QQCgehmqIOzMqcHK3N6aE+ZFy1atElVuzZMt6DTjAEDBrBw4cI2bVtaWkpJSUl0M9TBWZnTg5U5PbSnzCKyOly6XV4zxhgTNxZ0jDHGxI0FHWOMMXFjQccYY0zcWNAxxhgTNxZ0jDHGxI0FHWOMMXGTdkFHRAIiskhETorpgZbMo89Xc2N6CGOMSTZJH3RE5D4R2SginzZIP0FEPheR5SIyJWTRZOCRmGds+Yv0X/0o2HhFxhhTK+mDDjAHOCE0QUR8wJ3AOOAA4GwROUBExgKLgQ0xz1XXIWTtKYPyb2N+KGOMSRZJ3w2Oqr4qIgMaJI8AlqvqSgAReRgYD+QDAVwg+k5E5qtqTUwy1nWIe/12KeR3i8khjDEm2SR90ImgN/BVyPxaYKSqXgIgIpOATZECjohcAFwA0L17d0pLS1udAX/FFkYBy958mq9XxyaudURlZWVtOl/JzMqcHqzM0ZGqQUfCpNXeXFHVOU1trKp3A3cDDB8+XNvU4Z0qu9+/kn2zN7FvGnUSaJ0ipgcrc3qIRZlT4Z5OOGuBviHzfYCv45oDEbYUD4MVL0F1VVwPbYwxHVWqBp33gMEiMlBEsoGfAE/FOxObuoyEiu3w+bPxPrQxxnRISR90ROQh4C1giIisFZFfqOoe4BLgeWAJ8IiqfhbPfO2pruHrwqFQ2A/euA1q0ue+jjHGRJL093RU9ewI6fOB+XHOTq2rHv+YV5dUcvz4qTD3IvjoITjsZ4nKjjHGdAhJX9PpqPyZPiqrgYN/An2PgGevgo1LE50tY4xJKAs6MZKTlUFVjUJGBpw5G7Ly4MEzYeuqRGfNGGMSxoJOjORk+aiq9mY69YKfPgwVO+C+cbD+k4TmzRhjEsWCToz4MzPYo1Bd4z0e1HsYTHoGtAbuGQsLZ1u/bMaYtGNBJ0ZysnwA7N5TXZfY4yC48HXoNwrmXQb3j4dNXyQmg8YYkwAWdGIkJ9Od2oqqBk2l87vChCfghzfD1x/CrFHwzJWwfV38M2mMMXFmQSdG/OFqOkEZGTDifPj1QteMetFsmHmoCz6bV8Q3o8YYE0cWdGIkJytCTSdUfjc4eQb8ehEcfJYLPrcPhQfOgC9esAdKjTEpx4JOjORkuppORVWYmk5DRQNg/B3wP5/B0VNg/cfwf2fAjIPhpT/BpuWxzawxxsSJBZ0Yyc12QWdXZQuCTlBBDzhmKlz2KZxxnxuT5/Vb4I5hrsXbu/8PyjbGKMfGGBN7Sd8NTkcV8LtTu6tyT+s3zsyGg053045v4JNHXTc686+E+b+F/mPggPGw/8nQqWeUc26MMbFjQSdG8ryaTvnuVtR0wunUE8b8Bkb/GjYugcVPuunZ37qudfqOhAN+BEPGQfHeUci5McbEjgWdGAlkt6OmE44IdD/ATcdMhW8/rwtAz//OTV32hX1PcAGozwjw2cdrjOlY7FcpRvL8bbin0xpdh8DRV7lpy0pY9rwbt+ftWfDmTMgtgn1+AENOgH3GQk5hbPJhjDGtYEEnRvKiXdNpSvHecMRFbqrYDitehs+fgy/+A588AhmZrqfrfY6FQcdBj4Pds0LGGBNnFnRiJM97OLSsvfd0WiunEA481U011bD2PVj2HCx/EV661k2BrjDIC0CDjnW9JBhjTBxY0ImRjAzB74Ndu+NQ04mYCR/0O8JNY6fDzg2uFrTiJReEPv6XW6/nIS4A7XOcuxeUmZ24PBtjUpoFnRjy+4TyWN3TaYuC7nDo2W6qqYH1H8Hyl1wgenOmeyYoOx8GfB/2PhoGHg3d9neNGIwxJgos6MRQTiaUJ7Km05SMDOh1mJuOutKN9fPlq64WtLIUlj3r1gt0hYFHedPRUDwwodk2xiQ3CzoxlOOT+DQkiIacTrD/SW4C2LbGBaEvX4WVr8Cnj7v0zv28AFQCA7/velEwxpgWsqATQzmZUNZRazrN6dwPDpvgJlU37s+Xr7hpyTz44AG3Xtf96mpC/UYnNs/GmA7Pgk4M5fik/T0SdAQi0HVfN40437WKW/9xXS3ogwfg3bsBODyvH5T9APqPdt31WDc9xpgQFnRiyJ8Jm5Pl8lprZPjq7geNuRT2VMLX78PqN6j4YB6Bjx+Bhfe6dYsGwoAxLgD1Hw2d+1vDBGPSmAWdGMrxCbt2pUBNpzmZ2bVNsz+pHkbJ94+EDZ/A6jfdtPSZustxnXp7taDR7nJcl33tQVVj0ogFnRjq0K3XYsmXWVcTGvUr1zz726Ww+g0XhL581fWcDZDTGfoc7jou7Xs49B4G/oKEZt8YEzsWdGIoxyeUV1ahqkg6X1LKyKjrrHTE+a5hwpaVLgCtfRe+es89rIqCZEC3A10A6jvSBaTive2SnDEpwoJODOVkQo26IauDg7oZXADZa5Cbhp7j0r7bBusWugD01TvwyWOw8D63LK8L9B3hAlCf4dDzUNfE2xiTdCzoxJDf5/47L9u9x4JOc3I7u96w9xnr5muq3SW5r951/cd99Q58Pr9u/b0GQ++h3mW8odDje5Cdl5CsG2NazoJODOV4Z/e7jtQVTrLI8EH3A900/FyXtmuLayW37gP3+uWrdf3Hic912dPrUBeEeh0G3Q+yfuSM6WCiGnRE5ATgOuAAYBPwMHCtqpY1WO9nwP2qmtL//gdrOuWp2Gw6EfKK69eGwA3n/bUXhL7+AJbOr2spl5HlHl7t8T3ocZALQj2+5/ZjjEmIqAUdERkGzAM2AE8AvYArgNNE5GRVXRKtYyWLHC+kJk1XOMmoU0837fdDN6/quvD5+n34+kPY8KnrT+6jB0O26V0XgHoc5MYXKhpoTbeNiYNo1nSmAyuAkaq6DUBExgAPAq+JyImq+k4Uj9fh5WR6NZ1U6JUgWYhAUX83HXhqXXrZRlj/iQtC6z91r8tfBPU+m6yAG4216371Xzv3t2BkTBRFM+gMBW4JBhwAVX1DRIYDzwIvishpqvpCFI/ZofmtptNx5Hdz4wXtc1xdWlWFa6yw4VMXkL5dCisX1K8VZeZCl8HuflFtMNrPBSOf3RI1prWi+VfTCdjWMFFVvxWREtylt6dF5JwoHrNDs5pOB5eV4zU8OLR++nfbYNMyF4S+/dx7sPXNukYL4O4XFfWH4kEM+i4b8r5wTcCLB0FhH9cQwhjTSDSDzpfAweEWqGqZ18jgcdzltpeieNwOyxoSJKnczu65oL4j6qfv3umC0calsPkL2LwCtqyk17dfwNqn69bz+aFogBeE9navRQNc7aiwD2T641gYYzqWaAadUuBMEblCVRv9yqpqhYiMxwWdMwCN4rE7pLqGBFbTSQn+AtdNT+9h9ZJfW7CAkmFDvCC0ojYYsXm5G5m1enfI2gIFPd3QEUX93Wvnfi4gde7ngpIvK77lMiaOohl07gd6AsOBt8OtoKp7ROQs4G/AYVE8douIyCnAiUA34E5V/U8sj5ftc/e107L/tXQiAp16uWng9+svq6mBHetci7ptq73XNbB1Nax+y/VBpzUh+8pwresK+7jXTr0av+Z3s8t3JmlFLeio6kLgzBasp8Dl0TquiNwHnARsVNWDQtJPAGYAPuAeVb1RVecCc0WkCLgZiGnQEREC2Zl2TyedZWRA575uYkzj5dVVIUEpJCBtXwvrFsGSpxvUlHAPwhb0rAt0oVOBF5QKekB2IC5FNKY1UqH5zRzgDlxNCwAR8QF3Aj8A1gLvichTqrrYW+X33vKYC/h9VtMxkfmy3P2eogHhl6u6nhh2rIMdX4e8eu83fAZf/AeqdjXeNjsf8rt7kxeI8rtBfo/6aXl7Wc3JxI24ikdyE5EBwLxgTUdERgHTVfV4b36qt+qN3vSCqr7YxP4uAC4A6N69+7CHH364TfkqKyvjuvcz6Ncpg4sPzWnTPpJNWVkZ+fn5ic5GXCW8zKpk7inHv3sT2ZVb603+3VvJrtxGduUWsiu3kVndODgpGVRmF1KZXURldhFVWZ2ozC6kKqvQe9+ZqqxOVGUVUpldSI0vJ/FlTgArc+scc8wxi1R1eMP0VKjphNMb+Cpkfi0wEvg1MBYoFJF9VPWucBur6t3A3QDDhw/XkpKSNmWitLSUrkWZBPKzKSkZ0fwGKaC0tJS2nq9klVRlrtwFZRtCpo3IzvX4yzbgD6aVfwGbN8GeivD7yMqjIiOfnL36QKCL6wU8EJy61p/P65IyHbEm1eccJbEoc7NBR0RifUNCVTXawS/c4CuqqjOBmVE+VpPysn2UW+s101Fk50HxQDc1RRUqy6H8W9i12b2Wb6qd37byU3rk+lyQ2vCZW9bw3lNQZo67hJdbDHlF3mtxSFpxSJr3PqfQxlBKUS35sY/1Jx+L/a8F+obM9wG+jsFxmhXwZ7JhR4T/GI3pqETAn++mMAFqaWkpPUL/A1Z1zzHt2gTlXpDatckFo++2wK6tLnh9t8UFqe+2wHdb67fcq3d8H+QWNQhOzQSs3CJrbp4Emg06qpqMHU+9BwwWkYHAOuAnwE8TkZGAP9Oe0zGpT8QNrJfTyT0Q2xI1NVCxzQWfXVu84LSlLjiFpm1d5Tpx3bUlco0KXOOJ3CI3DHpucArOF0We93eymlWcJP09HRF5CCgBuojIWmCaqt4rIpcAz+OaTN+nqp8lIn+BbB9l1nrNmMYyMuouqe01qGXbqLqWeo2C09aQaZt7rdgGm76om28qWElGSKAqChuYenyzHpbuahy4stKjkVC0JH3QUdWzI6TPB+aHWxZPAX8muyzoGBMdIu75o+yA9+xTK1R9Vz8oBQNTw0AVrHltXuHNbwOU/QA+D3NLODOnidpUZ/eaU+i9L6w/n5mTdjWspA86HV3Aa0hQU6NkZKTXl8uYDiUr102derVuu5oa2L2Dt0uf44iD940cqL7b5qZta+Cbj1xauOenQvmymw5KTc3nFCbl81UWdGIs3xuzeldVNfl+O93GJJ2MDMjtTEVuD+g9tHXb7qmEiu3etK2u5lQ7v73+/K7NXg3L20abuR/s71Q/CLUmcGXlJqSWZb+CMZbvd61pyir2WNAxJt1kZkN+Vze1lipUljUOTE3Nb1lZF7Aqy5ref0ZWs0HKX9G79fluhv0KxliwplO2uwqwG47GmBYScT2b+wtcB7CtVV0FFTtaVsOq2O4uB25dVZdWswf/YTdGrzyeVH04tMMo8Go3OyusMYExJo58WRDYy02t5bUS3Pn6W1HPVqo+HNph1NV0LOgYY5KE10pQM6JfH0jVh0M7jHyr6RhjTC0LKDFWkBMMOlUJzokxxiSeBZ0YK8hxrdespmOMMVEKOiKSJyLpNdBECwUvr+0ICTo1NUpFlfXHZoxJP20OOiJSJCIzvP7OdgLbRaRMRF4Xkaki0j162Uxevgwh359Z7/LalY9+xH5/eI5UGEDPGGNao01BR0R6AAuBS4BeuBZoAuQBo4HrgJUicqOIpH1f451yMtlcVsmry74F4IkP1gHwxcZmHt4yxpgU09aazp+BgcAe4HbgaGAfYDhwMbAA9yTkb4E3vCCVtgrzsnnqo6/5+X3v8s7KzbXpG3c00eutMcakoLY2wj4RUODX3tDOod4H7hKRw4B7cYFovoiMVNW0bMJVlFdX2ftk3fba95vLLegYY9JLW2s6nXC1nDmRVlDVD4AjcUHoEFwNKC0V5WXXvl+/vW4U0a3llYnIjjHGJExbg84KoEpVm/zVVNVduPs+ApzTxmMlvaJAXU1n5aby2vdbdqVlxc8Yk8baGnQeBnJF5IjmVlTVt4EycGMgpaPuBXUdfX6xcWft+zJ7dscYk2aaDTpe8+dxItIzJPkWYBkwU0Sa7DpZRHy4e0flTa2XygZ3r3uE6ast39W+dz1PG2NM+mhJTefPwDxgrYhsEJHngem4+zn7Aq+IyIFNbH8mriXb8+3LavI6/sAePH/ZUZx2WP2xKawTUGNMumlJ67UncA0B9ga6Aj8AxoYsHw58KCIvAP8GFgFbgG64Vm5XAl8BU6KX7eQiIgzpUUDnkAYFBTmZ1jWOMSbttKSX6TMAvG5uDvGmQ73pIFwtxgecABwfZhdbcIEn7R+/D2063acor17XOMYYkw5a/JyOqpYBb3gTACKSgWsgEBqIDsHVcoL2wjU80NYcLxUV59fVdHp3zmXV5rS9zWWMSVPtCgKqWgMs9qaHguleo4PQQHQYMKg9x0oFoa3YigNZfLLOGhIYY9JLTGoeqvoN8A3wXDBNRHJjcaxk0r1TXdApyMmyJtPGmLQTt/F0VPW75tdKbXt3DQAwuFs++f5Myiurqa5J+1tdxpg0ktb3WOIt4M/kxcuPIuDP5JmPvwFcs+nC3LTviNsYkyYs6MTZPt0KgPrDWFvQMcakCxuuOkEC3oii5bttBFFjTPqwoJMgwaBjvRIYY9KJBZ0Eya+t6VjQMcakDws6CRLItqBjjEk/UQs6IvInEflztPaX6oINCezymjEmnUSzpjOFNO7Us7Xsno4xJh3Z5bUECfh9gF1eM8akFws6CeLP9JHty6DMmkwbY9JIWj0cKiIBYBZQCZSq6v8lMj8Bv89GDzXGpJWkr+mIyH0islFEPm2QfoKIfC4iy0UkeK/pNOAxVT0f+FHcM9tAwJ9pD4caY9JK0gcd3LDZJ4QmiIgPuBMYBxwAnC0iBwB9cKOYAiT81z7fb6OHGmPSi6hGp5djEakCMlTVF5Udtu7YA4B5qnqQNz8KmK6qx3vzU71V1wJbVXWeiDysqj+JsL8LgAsAunfvPuzhhx9uU77KysrIz8+PuPzPb39HZgZMHpE6oz40V+ZUZGVOD1bm1jnmmGMWqerwhumpek+nN3U1GnDBZiQwE7hDRE4Eno60sareDdwNMHz4cC0pKWlTJkpLS2lq2/tWvsv2XZWUlBzZpv13RM2VORVZmdODlTk6UjXoSJg0VdVy4Nx4ZyaSfL+Pr7cl/CqfMcbETSrc0wlnLdA3ZL4P8HWC8hJRXnamPadjjEkrqRp03gMGi8hAEckGfgI8leA8NZLvt6BjjEkvSR90ROQh4C1giIisFZFfqOoe4BLgeWAJ8IiqfpbIfIaTl+2jvLKaaDXmMMaYji7p7+mo6tkR0ucD8+OcnVYJ+DOprlF276khJyvujf6MMSbuolnTCXfz3jQh3zr9NMakmWjWdK7FAk+rBEIGcuuS709wbowxJvaiFnRU9dpo7Std5Hs9TVtNxxiTLpK+IUEyq6vp2LM6xpj0YEEngUIvrxljTDqwoJNAwYYE5ZUWdIwx6cGCTgLles2kd1Xa5TVjTHqwoJNAwctru+zymjEmTcQk6Hhdz5hm5GV7NZ0qq+kYY9JDrGo698RovynFn5lBhsAua71mjEkTsQo6p4nIVU2tICLjY3TspCEi5GVn2j0dY0zaiFXQ+QXwZ2+wtEZEZDrweIyOnVTysn3sstZrxpg0EZMOP1X1XyJyCPCgiIxS1cUAIhIA/gmMB/4ci2Mnm4A/03okMMakjVj2Mn01cDDwtIgMB4qBuUB/4HRVnRvDYyeNgN9nD4caY9JGVC6vicjVIvIjERkYTFM3SMzZwG7gBeBdIBsYaQGnTiA707rBMcakjWjVdKYBPgARKQc+Az4BPgZmAjOA/wA/U9UdUTpmSsj3Z/LN9opEZ8MYY+IiWkEnH9gfdznte97ricB5gOKGPDgMeFhEPsYFo49V9dMoHT9pBfyZ1g2OMSZtRCXoqGol8JE31RKRYuoHou8BvwICuGCU9sNlBvyZdk/HGJM22hR0RKSrqn7b3HqqugUo9abgtgLsjQtAaS/f77PWa8aYtNHWhgRvi8jgtmyozgprTOAE/JlUVNWwp7om0VkxxpiYa2vQGQi8ISKjopmZdFQ3vIG1YDPGpL62Bp0HgC7ASyJyWnMri8iJIrKojcdKaTaQmzEmnbQp6Kjqz4HrgBzgERG5NNx6IlIiIm8ATwGHtjWTqcyCjjEmnbT54VBV/SOuj7Vq4BYRuSW4TERGisiLwEvAKFxLtUfamdeUlO93DfisMYExJh20q8m0qs4WkbXAY8ClIrI3LpCdiHs2pwb4F/AnVV3S3symonx/FoD1SmCMSQvtfk5HVV8QkQnAk8DJwWRcsLlWVZe29xipLFBb06lKcE6MMSb22tX3mogMEJG7cTUdcLUbAT4ELraA07wCr6ZTZjUdY0waaFPQEZGBInIP8Dnuvk4WsAA4HfgC1+XNGyLSP1oZTVW1NZ0Kq+kYY1JfW2s6S4FzccHmHeA4VT1OVf+NazjwFq4vtrdEZFhUcpqi8nPsOR1jTPpoa9DJwvWzdrKqjlbVBcEFXtc3x+FGBu0BvCIiJ4ffjfFn+sjyCTsrrPWaMSb1tTXonKWqQ1X1mXALVXW3qp4J3AbkAU+IyCVtPFbKC/gzbchqY0xaaOvDoY+2cL3LgctwjQtua8ux0oEN5GaMSRdRGTm0Kao6EzgTN4KoCSMv22c1HWNMWoh50AHwGhgcG49jJSMXdKymY4xJfdEaOTQsERmAG7wtOP04lsdLVnnZdk8n1ezevZstW7awc+dOqqtT7x+KwsJClixJr05G0rnMPp+PgoICiouL8fv97dpnVIKOiARwg7IdQl2A+R5QEI39p7qA38fX2+w5nVSxe/du1qxZQ1FREQMGDCArKws3dmHq2LlzJwUF6fXnna5lzs/Pp6qqih07drBmzRr69evXrsDT6qAjIoNwQSU0wAzANRaot2rI+y3AsrZlMbpE5BRc33DdgDtV9T+JzZGr6ZRbTSdlbNmyhaKiIrp06ZLorBjTbiJCdnZ27fd5y5Yt9OzZs837a9E9HRG5S0TeEpGduODxGPAH4BTc0NMZuCDzFXXB5iJgDNBVVbuo6ug257IuH/eJyEYR+bRB+gki8rmILBeRKU3tQ1Xnqur5wCTgrPbmKRoCfp+1XkshO3fupFOnTonOhjFR16lTJ3bu3NmufbS0pnMBrhPPYED5ClgMfBYyLVbVMhEJjrv8kKruaFfuGpsD3AHcH0wQER9wJ/ADYC3wnog8BfiAGxps/9+qutF7/3tvu4QL2D2dlFJdXU1WVlais2FM1GVlZbX7HmVrL6/dDUxW1e3tOmobqeqrXuOEUCOA5aq6EkBEHgbGq+oNwEkN9yHu4vqNwLOq+n6444jIBbhAS/fu3SktLW1TfsvKylq07bffVLKrspqXFywgI8mv/be0zKmkYZkLCwspKytLXIbioLq6ut3/8SYbK7NTUVHRrr/xlgadK3A1g/OBUSIyRVWfbfNRo6s3ruYVtBYY2cT6vwbGAoUiso+q3tVwBVW9GxdgGT58uJaUlLQpY6WlpbRk22UZK3hyxVJGjP4++f6YNiiMuZaWOZU0LPOSJUtS/oZzut5UtzJDTk4Ohx12WJv32aJ7Oqp6K7AP7tLWfsA8EXlZRIa3+cjRE65qoJFWVtWZqjpMVS8MF3ASIS/bhqw2xqSHFj8cqqpbVfVS4CBgHlACvCMi/xKRfWKUv5ZYC/QNme8DfJ2gvLRJsHZjQ1YbY1Jdq3skUNUvVHU8roeBj3Bd3HwmIneISLdoZ7AF3gMGe2P8ZAM/AZ5KQD7aLOC3mo4xJj20uRscVS0FhuEGcdsEXAwsj062whORh3Bj9QwRkbUi8gtV3QNcAjwPLAEeUdXPYpmPaLOajkllS5cu5de//jUHHXQQhYWFZGdn06tXL0488UTuvfdeKioqGm2zbNkyLr/8coYOHUpxcTFZWVkUFxczcuRIrrzyShYtWpSAkphoaFffa+rMBgYDf2qwv1+JSF579h/meGerak9VzVLVPqp6r5c+X1X3VdVBqvrnaB4zHmqDjo2pY1LMtddey4EHHsgdd9xBQUEBEydO5Morr2TcuHEsXbqU8847jyOPPLJ2fVXlmmuuYf/99+fWW29FRDjrrLO46qqrmDBhArm5udx+++0MHz6cO+/sEE88mFaKSlMpVd0FTBORv+OejZkAXAdcJiJ/BWZ565gw6kYPtaBjUsf111/PtGnT6Nu3L48++igjRzZuVDpv3jz+9re/1c5fe+21TJ8+nb59+/LQQw8xZsyYRtts3LiR2267je3bE/LkhmmnqLbPVdWvgYkiMgO4BTgK+AvwWxH5q6reHM3jpQqr6ZhUs2rVKqZPn05WVhbz58/noIMOCrveSSedxA9+8AMAVq5cyXXXXUd2djbPPvssBx54YNhtunXrxvXXX8+ePfb3koxiMrSBqr6vqiXA6cAKoCsu+JgwgkFnp93TMSli9uzZVFVVcfrpp0cMOEHBziNnz57Nnj17OOOMMyIGnFCZmcn9TFu6iumnpqr/FpGncQ9k/j6Wx0pmOVkZZAjssv7XUt41T3/G4q+j3TtUdB3QqxPTTm7+R78pr7/+OgDHHXdci7d54403ADj2WBt6K5XF/F8Fr3XZrSIyO9bHSlYi4oastns6JkV88803APTp06fF26xfvx6A3r17N1q2atUq5syZUy+tc+fOXHbZZW3Oo0mMuNVPVXVbvI6VjPL8PntOJw20twaRLFRdpyCtGUeoqW1WrVrFNddcUy+tf//+FnSSULP3dESkOsaT/dLiHhAttyGrTYro1asXAGvXrm3xNsExWtatW9doWUlJCaqKqlJVZQMeJrOWNCSQGE8xacyQbALZmVbTMSkj+OzNSy+91OJtgs2jW7ONST7N/uCrakasp3gUtKML2OU1k0LOPfdcsrKyePzxx1m8eHGT6+7evRuASZMmkZmZyWOPPcaSJUvikU2TAPaD30G4mo5dXjOpYcCAAUyfPp3KykpOPPFEFi5cGHa95557jnHjxgEwaNAgfv/731NZWcm4ceN48803w26zbdu2WGXbxIE1dO8gAn4bPdSklt/97nfs2bOHa665hsMPP5zRo0czfPhw8vPz2bBhA6+++ipffPEFw4fXjZDyxz/+EVXlT3/6E2PGjGHYsGGMGDGC4uJitm3bxqpVq3jxxRcBOOqooxJVNNMOFnQ6iIDfR5nVdEyK+eMf/8iZZ57JrFmzWLBgAbNnz6aiooK99tqLQw89lMmTJzNhwoTa9UWE6dOnc/bZZ3PXXXexYMECHnzwQcrLyykoKGDQoEFcdNFFnHPOOQwdOjSBJTNtZUGng7CGBCZV7b///tx+++2t2mbIkCHceuutMcqRSSS7p9NBBPyZfFdVTXVNxEFPjTEm6TVb0xGRWF/zUVVN+xpXwO8DYFflHgpyshKcG2OMiQ17TqeDCI4eusseEDXGpLBmaxj2HE18BLJtyGpjTOqzgNJB5GUHL69ZTccYk7os6HQQtQO5WU3HGJPCLOh0EHm193Qs6BhjUpcFnQ4i4F1es65wjDGpzIJOBxGs6VhDAmNMKrOg00HkB1uvWUMCY0wKs6DTQeQGW69ZTccYk8Is6HQQ2ZkZZPsyrKZjjElpFnQ6kIDfZ63XjDEpzYJOB5KXnWnP6Zi0VlpaWju8QVuVlJQgItHLVJTNmTMHEWHOnDmJzkpCWNDpQPL9NryBSR0iUm/y+XwUFxdTUlLCnDlzULUe1dNR2vfu3JEE/D57TseknGnTpgFQVVXF8uXL+fe//80rr7zCwoULueOOO+qtO2LECJYsWUKXLl0SkVUTBxZ0OpD8nCx2fFeV6GwYE1UNL5W98cYbHHXUUcyaNYsrrriCgQMH1i7Ly8tjv/32i3MOTTzZ5bUOJN/vs3s6JuWNGTOG/fbbD1Vl0aJF9ZY1dU9ny5YtXH311Rx00EHk5eVRWFjIIYccwpQpUygvL2+0/p49e7j++usZPHgwfr+fvn37MnnyZCorKxutO3fuXCZMmMC+++5LIBAgPz+fYcOGMXPmTGpqahqtP2nSJESElStXcvvtt3PwwQeTm5tLSUlJ7TrLly/nzDPPpKioiEAgwOjRo3nmmWcinpePP/6Ys88+mwEDBuD3++natStDhw7lsssuo6qq/j+j27dvZ+rUqQwZMoScnByKioo4/vjjefHFFxvtN/Scfvjhh5x44ol07tyZvLw8jj76aN58882IeYoFq+l0IHZPx6SL4P2crKyWDVj45Zdfcswxx7B69WqGDRvGRRddRE1NDcuWLePWW2/lwgsvJBAI1Nvmpz/9Ka+99hrjxo2jU6dOzJ8/n5tuuomNGzcye/bseutOmTKFjIwMRo4cSe/evdm+fTsvv/wyl156Ke+99x7//Oc/w+br0ksv5bXXXuPEE0/khz/8IT6fe97uiy++YNSoUWzevJlx48Zx6KGHsnz5ck455RTGjRvXaD8ff/wxI0eORET40Y9+xMCBA9mxYwfLly9n1qxZXHfddbXnatu2bYwZM4bFixdz+OGHc9lll7Fp0yYeeeQR/uu//ov//d//5Ze//GWjYyxcuJCbbrqJUaNGcd5557FmzRoef/xxjjvuOD788EOGDBnSos+ivSzodCABfyZlFRZ0UtqzU2D9J4nORdN6fA/G3Riz3b/66qt8/vnnZGdnM2LEiBZtM2HCBFavXs3111/P1KlT6y3btGkT+fn5jbZZsWIFn332GcXFxQD8+c9/5pBDDuH+++/nhhtuoEePHrXrPvPMMwwaNKje9jU1NZx77rncf//9XHLJJRxwwAGNjvH+++/zwQcf1LtECPCrX/2KzZs3c9ttt3HppZfWpj/55JOccsopjfbzj3/8g4qKCubOncv48ePrLdu6dSt5eXm185MnT2bx4sVccMEF3HXXXbUt9SZPnszw4cP5zW9+w/HHH8+AAQPq7eeZZ55h9uzZTJo0qTbt73//OxdeeCEzZsxg1qxZjfIVC3Z5rQPJ92dSVrnHWvWYlDJ9+nSmT5/O1VdfzVlnncXYsWNRVW6++WZ69uzZ7PaLFi3izTff5NBDD2Xy5MmNlnfp0oWcnJxG6X/5y19qAw5AIBDgZz/7GTU1NSxcuLDeug0DDkBGRkZtwHj++efD5u2qq65qFHDWrl3LCy+8wMCBA7nkkkvqLRs/fjxHH310hJJCbm5uo7SioiIyMtxPdVVVFQ888AD5+fnccMMN9ZqGDx48mN/85jdUVlZy//33N9rPmDFj6gUcgP/+7/8mMzOTd999N2Keos1qOh1IwJ+JKlRU1dR2i2NSTAxrEB3VNddcU29eRLj33ns599xzW7T922+/DcDxxx9f++PbEsOHD2+U1rdvX8DVHkJt3ryZv/71r8yfP5+VK1c2uke0bt26sMcIV1P74IMPADjyyCNrL7eFKikp4ZVXXqmXdtZZZzFjxgxOOeUUzjjjDMaOHcuYMWMaBcOlS5eya9cuxowZUy+gBh177LFcd911tXkIFe58ZGVl0b1790bnI5bSrqYjIgERWSQiJyU6Lw0FhzewxgQmlagqqkpZWRkvvPACffv25cILL+Tll19u0fbbtm0DoHfv3q06bufOnRulZWa6/7Orq+seTdi2bRuHH344f/nLX8jNzeXnP/85V199NdOmTaut6ezevTvsMUIv0QVt374dgO7du7d4mxEjRvDaa69x7LHH8thjjzFx4kT22Wcf9ttvPx566KFG+45UQwymB89ZqHDnA9w5CT0fsZY0NR0RuQ84CdioqgeFpJ8AzAB8wD2q2ty/kpOBR2KW0XYI1BvIzZ/YzBgTZYFAgLFjx/L0008zdOhQJk6cyOeff17vfkU4wR/LSLWN9rrnnnv48ssvmTZtWqNWc2+99RYzZsyIuG24ng8KCwsB2LBhQ9ht1q9fHzZ91KhRzJs3j927d7No0SKee+45br/9dn7605/StWtXxo4dW7vvSPv45ptv6uWhI0qmms4c4ITQBBHxAXcC44ADgLNF5AAR+Z6IzGswdRORscBiIPy3IcHygsMb2AOiJoUdfPDBnH/++axdu5Zbb7212fWPOOIIwN1XCdd8ub2WL18OwOmnn95oWcPLYC1x2GGHAfD666+HrUGUlpY2ub3f72f06NFce+21zJw5E3ANEACGDBlCXl4eH374YdhLYgsWLABg6NChrc53vCRNTUdVXxWRAQ2SRwDLVXUlgIg8DIxX1RtwtaJ6ROQYIIALUN+JyHxVbfQtFpELgAvAVZGb+5JEUlZW1qptV2xyX9DX33mPjUXJeU+ntWVOBQ3LXFhYyM6dOxOXoTiorq5ucRnDrXfZZZcxZ84cbr75Zs455xyKiooA2LVrF+AuZwW323fffRk5ciTvvPMO1157LVdccUW9fW3evJlAIFDbmCD4Qx/uuBUVFbWvweXBS1LPPfdcvRZfH330ETfccAPgbuCHljn43ExZWVmj4xQWFnLMMcewYMECbr75Zi6++OLaZc8880xtIAvNw5tvvsmBBx7YqIayZs0awF0CC6774x//mDlz5jBlyhRuvvnm2nVXrlzJzJkzycrK4tRTT61dP9w5DRVsuBRuWbjPuaKion1/48HrrckwAQOAT0Pmz8BdUgvOnwPc0YL9TAJOaskxhw0bpm21YMGCVq2/aPUW7T95nr68dEObj5lorS1zKmhY5sWLFycmI3G0Y8eOZtcB1P3EhHfppZcqoFOmTKlNW7BggQI6bdq0euuuXLlS+/Xrp4AOGzZMr7jiCr388sv15JNPVr/fr19++WXtukcffXTE486ePVsBnT17dm3aunXrtLi4WDMyMvTUU0/Vq666Sk899VTNysrSs846SwGdOHFivTJPnDhRgXrHDbVs2TLda6+9FNAf/vCHOnXqVD3zzDM1MzNTTz755EZ5GD9+vObm5uq4ceP0V7/6lU6ePFlPOukk9fl8WlRUpMuXL69dd/PmzbrffvspoCNHjtQpU6boeeedp506dVIR0TvvvLNeXiKd06D+/ftr//79wy4L9zm39PsNLNQwv6nJdHktnHBdyTbb3lhV56jqvBjkp10C3uW1XXZ5zaSBqVOnkpeXx8yZMyPe/wgaOHAg77//PldddRU7d+7kjjvu4N5772XNmjVcccUVdOvWrc356NWrV+0Dnq+//jp33HEHq1evZtasWdx4Y9taGw4ePJi3336b008/nTfeeIMZM2bw1VdfMXfuXE477bRG61988cWcddZZrFq1in/+85/cfvvtLFu2jIsvvpgPPvigXiu24uJi3nrrLa666io2b97MLbfcwqOPPsqIESN47rnn6tWsOqRwkaijTjSu6YwCng+ZnwpMjeYx41nT+WpLufafPE//9e6aNh8z0aymYzWdVGVldtK9pvMeMFhEBopINvAT4KkE56nN8r3Wa+U2kJsxJkUlTdARkYeAt4AhIrJWRH6hqnuAS4DngSXAI6r6WSLz2R51rdcs6BhjUlMytV47O0L6fGB+nLMTE9mZGWT7MiizezrGmBSVNDWddOEGcrOajjEmNVnQ6WDysjPtno4xJmVZ0Olg8v2Z1mTaGJOyLOh0MHl+n9V0UoBrMWpMaonG99qCTgdjo4cmP5/P12h4YWNSQVVVVdjhGlrDgk4Hk5ftsw4/k1xBQQE7duxIdDaMibodO3ZQUFDQrn1Y0OlgAtmZNp5OkisuLmbr1q1s2rSJyspKu9RmkpqqUllZyaZNm9i6dWvYweNaI2me00kXQ/sX4c9Kzh6mjeP3++nXrx9btmxh1apVcR0gK14qKirCDhGdytK5zD6fj4KCAvr164ff376xvizodDATjuif6CyYKPD7/fTs2TPiCI/JrrS0tHbcmHRhZY4Ou7xmjDEmbizoGGOMiRsLOsYYY+LGgo4xxpi4saBjjDEmbizoGGOMiRsLOsYYY+LGgo4xxpi4Eeuio2ki8i2wuo2bdwE2RTE7ycDKnB6szOmhPWXur6pdGyZa0IkhEVmoqsMTnY94sjKnBytzeohFme3ymjHGmLixoGOMMSZuLOjE1t2JzkACWJnTg5U5PUS9zHZPxxhjTNxYTccYY0zcWNAxxhgTNxZ0YkREThCRz0VkuYhMSXR+okFE+orIAhFZIiKficilXnqxiLwgIl94r0Uh20z1zsHnInJ84nLfPiLiE5EPRGSeN5/SZRaRziLymIgs9T7vUWlQ5v/xvtefishDIpKTamUWkftEZKOIfBqS1uoyisgwEfnEWzZTRKTFmVBVm6I8AT5gBbA3kA18BByQ6HxFoVw9gaHe+wJgGXAAcBMwxUufAvzFe3+AV3Y/MNA7J75El6ONZb8ceBCY582ndJmBfwDnee+zgc6pXGagN/AlkOvNPwJMSrUyA0cBQ4FPQ9JaXUbgXWAUIMCzwLiW5sFqOrExAliuqitVtRJ4GBif4Dy1m6p+o6rve+93Aktwf6zjcT9SeK+neO/HAw+r6m5V/RJYjjs3SUVE+gAnAveEJKdsmUWkE+7H6V4AVa1U1W2kcJk9mUCuiGQCecDXpFiZVfVVYEuD5FaVUUR6Ap1U9S11Eej+kG2aZUEnNnoDX4XMr/XSUoaIDAAOA94BuqvqN+ACE9DNWy1VzsNtwFVATUhaKpd5b+BbYLZ3SfEeEQmQwmVW1XXAzcAa4Btgu6r+hxQuc4jWlrG3975heotY0ImNcNc3U6ZtuojkA48Dl6nqjqZWDZOWVOdBRE4CNqrqopZuEiYtqcqM+49/KPC/qnoYUI677BJJ0pfZu48xHncZqRcQEJEJTW0SJi2pytwCkcrYrrJb0ImNtUDfkPk+uKp60hORLFzA+T9VfcJL3uBVufFeN3rpqXAexgA/EpFVuMukx4rIA6R2mdcCa1X1HW/+MVwQSuUyjwW+VNVvVbUKeAIYTWqXOai1ZVzrvW+Y3iIWdGLjPWCwiAwUkWzgJ8BTCc5Tu3ktVO4FlqjqLSGLngImeu8nAk+GpP9ERPwiMhAYjLsBmTRUdaqq9lHVAbjP8WVVnUBql3k98JWIDPGSjgMWk8Jlxl1WO0JE8rzv+XG4e5apXOagVpXRuwS3U0SO8M7Vz0O2aV6iW1Ok6gT8ENe6awVwdaLzE6UyHYmrRn8MfOhNPwT2Al4CvvBei0O2udo7B5/TihYuHXECSqhrvZbSZQYOBRZ6n/VcoCgNynwNsBT4FPgnrtVWSpUZeAh3z6oKV2P5RVvKCAz3ztMK4A683m1aMlk3OMYYY+LGLq8ZY4yJGws6xhhj4saCjjHGmLixoGOMMSZuLOgYY4yJGws6JumIyHQRURGZk+i8JIKIHC4iT4vIJhGp8c7F9BZsN8BbN22brNo5SDwLOilIROYE/7BEZGFT3Y6LyAPp/AOebERkMFAKnIR7dmYTsAEoa+d+S7xgfkp785goInKKV4aSROfFRJaZ6AyYmBsGnIrr1sMkvwtwPSC/BvxIXe/PLVWFe8gvnBJgGq6X4bltz15CnULdk/WlEdZp6hyYOLCgkx6uFZG5qlrT/KqmgzvQe32klQEHdT0p7xf1HCUROweJZ5fXUtsrwC7cD9VPE5wXEx253mu7LqcZkygWdFLbely/SADTvcGpWizkvtCACMsj3pQVkVJv2SQR6SQiN4nIChH5TkRWisi1IpITsv5xIvK8d3O8XEReFZHvtyCPGeKGGf7I226ziDwlIk0OqOVtd443PO+3IlIpIl+LyL9EZGSEbWobMHjbXyIi74rINi/90Oby2+D4vxCRV0Rki4hUiMiXInK3iOwTZv1V3nku8ZJmh3w+q1p4zEafVzANd2kNYGLIfiN+/iJysog8KSLrvXO30WvcEHbYZu97oCJS6s3/zCv7Zi/9FC/dJyLHiMgMEVkkIhtCPpt/i8ixYfZd4pUheGltWsMyNHUOwuzvMO9e51cistv7Tj4vIqc3sc0qb78l4oZ/vsX7PHeLyDoR+X/i9eQcZtsM7/ws8M5Hlfed/Ezc8NInRDpuUkp0B3Q2xaRTvzm4jjkfBoqB7d78+WHWfcBbNifMMvWmARGOMyC4Tphlpd6y/8H11qu4/84rQ/b7lLfuxbgB0qpD8qrAbmBMmH1P95b/AzfMguKu1W8L2XYPcFaEfBcAL4SsW9PguNXAJc0cd27IcbZ67w9t4eeTBzwfcrzKBnn/DhjfYJv3cP9EBM/fdm9+PfBeC4/b6PPCdV2/3vtsgsde32DqG7J+Vsh3RkPyEjp/U5hjT/KWlQIzQ87zFu/1FG+9gxrsqyIkb8Hpdw32PdrL53ch37N6ZWjJd9ZbfoGXn+CxtnqfcXD+n4QZlhpY5S2fEPK+3Mt/cNsvgaIw2/5fg/Jtw333g/NvJ/o3JZpTwjNgUww+1JCg481P9+bXAP4G68Y66GzD9dx7pJeeDZyHCxIK/AH3Q3o90Nlbpz/wprf83TD7nh6y7z24wBYc234Q8B9v+S5gUJjt/+0t/wjXS3Zw287AVO8PvpoGAS/kuDu9H5OLgDxvWTfcEL4t+Xzuou4H9ZfBzwTYF1gQ8oO1bxPndVIbvhdNfV7BsjX6HjRY79aQH9CzgXwvPR/3gx0MQGc32G5SyLmrAf4Y8nl3ArqFnINHcK3zukNtp8TdgN97n3cNMLKJ7/30Np6D0dQFnEeBPiFl+513XAV+H2bbVdQFqQ+AUV56JvAj6v4xuanBdkdRF4AvAwq8dAF64mpvN8f6NyOeU8IzYFMMPtTGQacTsNlLu7TBurEOOlXAPmGW3xuy//vCLO8f8kfer8Gy6SHbNho2AsjBBToF7mmwbCx1P5rFEcp1lbfOvCaOe0EbP5v+IT9svwyzPA83Fr0C9zdxXie14dhNfV7BsjX6HoSsM9jL+1Zg7wjr/Njbz6cN0ieFnLvr2/Hd/oO3j9lNfO+nt/EcvOQte53wtZnrqQucnRosW+UtWw/sFWbbK7zlKyN8155t6zlJtsnu6aQBdUNK3+TN/k7cePfx8qiqLg+T/mLI+xsaLlTV1bgfX3CXXMLZBdwWZtsK4G/e7Oki9Z5Tmui9zlHVLRH2+6D3eoyI+MIs3wzcF2Hb5pyGu5e6Hrin4UJV3UXdZ3VahOMnys9xeZ+rqisjrPMErqZ4YIR7GNXALWHSW+pp73VMO/bRiIgUA8d4szeoanWY1f6Cq53m42rI4dytqpvDpM/1Xgc2+PsLDvfeTUTS4vc4LQppALgd9xBhN+A3cTzuJxHSg0PiVlAXXBra4L0WRVi+UFXLIyx7xXvtjBv3Pmi09/o/3k3wRhNu8DJwtY69Ihx3T4TjNmeo9/pahB82gJe91wAwJMI6iRA8d2c0ce7W4u77QP2hjoOWq+qmpg4iIrniGoeUeg0UqkJu/n/grdYrKiWqcxjukpZS992pR1W3A4u82aHh1sHdewtnXcj7ziHvX8RdXh4KlIrIBBGJdtk6FHtOJ02o6i4RuR6YAfxWRGZ5f0Sx9k2E9OAP7gb1rjM0sU5WhOXrIqQ3XNYVCP5nHvzvu9CbmpMXJu3bFmwXSVfvtam8rw2zfkcQPHf53tScVp87r3ZUiru3E1SOu6RXA/iALriAHE3B87xdVZtqjh78bCJ9LjvDJapqRUiFOyskfbmIXIRrZfp9b8Jrkfgcrub0ASnEajrp5e/AV7iawxUJzkusRer6J/idH6+q0oJpVZh9RKqhtIa/iWWRgnCiBc/dpS08d6Vh9tHcubsNF3BWAqfj7rvlq2o3Ve0BHBGtwkTQ1OcSE6p6H642fhnwJO7y7QDgQmCRiPwu3nmKJQs6aURVdwN/8mYvE5EuzWwS/IHIibC8JTWFWGrqMkTo/YTQ/66Dl+wOiH52WiSYl/5NrBN6Wao9tapoi+m5E5FsYLw3+zNVfUJVtzZYrXssjk3dec4VkaZql30arB8VqrpBVWeo6im4WtQIXCtLAf4kIgdH83iJZEEn/cwGVuCeVZnSzLrbvNc+EZYfHqU8tdXhIhLuEg7A0d7rNlxLtaC3vNeID/rF2Pve68gm8h58ALKc+PUTFuwiKWLnsNSdu5NFJNIlz/boQl1NI9IlpbFNbN+SMkTyAXU1zGPCrSAihbi+DKHuc4w6dd4DzsRdzssAjozV8eLNgk6a8W6AT/dmL6bp2kKwEcD4hgtExI+7HJBIecClDRO9vF3uzT7W4J7RHO91uIj8vKmdi0ikBgzt8QTux3Ev3HMtDY+ZB/w2uG4TjQ2iLdiKqnMT6/wDl/deuOeZImrjudtB3Q//98Lssyfw62a2h6bLEJbXknGBNzs5Qkuyybhafxkwv7XHCMer3UXKUzXukQNIwGW/WLGgk54eBBbj+vEK+1+d5xHv9XwROdf7MUdEDsT90SW6lc123KWHS0UkF0BE9sZdF98f1zLuxtANVPU56nrcvk9Erglt2isiRSIyXkSepH1Ne8PymoLf7c3eKCIXhJzXfYFngH1wzcGvi/bxm/CZ93qkuOETGlHVJdQ1Ub9GRO70zjcAIpIvIj8QkX/iHq5sFe8G/tve7H3idSvkdRNzHK5VWVO1mGAZTojU5Uwz/oALqkOBh0Wkj3f8fO++SvDKwI3eYwjRcL2IPCZuWIbiYKKIdBeRmbh7PYrrQSM1JPpBIZuiP9Hg4dAI65xO/a435oRZJwv3IxBcp4q6J84342pAzT0cOinC8Uu85auayGPYfVC/O5onvPeV1D31rbgn138SYb8B6nolCE7baNydy+wIx210rlr5+eRR12tCuLxX0KAbnJae12aOO6CJzyuLuodSa3BN2ld5U5+Q9XzArAbnaQd1rcuCaQsa7H+Sl17aTB5H4gJucD9lIfPNfee6UPcQdDWu5eSq0O9YU+fAW/5L6h7ercF10xPaDc4DNN0NTkkTZQvuY0BI2m0NzuV26mp8wel3TZ2zZJusppO+nqCZ69KqWgX8APgr7o+qBnefYQ7u2vZHMc1h8xR33ftyXP9u2bgfv3nAaFV9OOxGquWqeiquq5UncM2Xc73tl+NqgmfgLj9GP9PuAdBxuO6AXsP9qOYBq3EPjH5PVZ+MxbGbyFMVcByub7F1uBaO/b0pM2S9alW9GHeP4QEvz9m487cGF8wn4sa2aUs+3gFG4R6m3IoLhhtxLS8PpYnvnLrnf47Bfabf4m7IB8vQ0uP/HXev8kFc0MrHBYIXgDNVdYJG95Lnrbjn5p4EluFqcn5cK9N/AUep6vVRPF7CBfs1MsYYY2LOajrGGGPixoKOMcaYuLGgY4wxJm4s6BhjjIkbCzrGGGPixoKOMcaYuLGgY4wxJm4s6BhjjIkbCzrGGGPixoKOMcaYuPn/GvMGqUD5QCIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - rhs) for x in conv_cg], label=\"CG\")\n",
    "plt.semilogy([np.linalg.norm(A @ x - rhs) for x in conv_rich], label=\"Richardson\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Number of iterations\", fontsize=24)\n",
    "plt.ylabel(\"$||Ax_k - b||_2$\", fontsize=24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "niters = 256\n",
    "roots = [np.cos((np.pi * (2 * i + 1)) / (2 * niters)) for i in range(niters)]\n",
    "lam_max = max(eigvals)\n",
    "lam_min = min(eigvals)\n",
    "taus = [(lam_max + lam_min - (lam_min - lam_max) * r) / 2 for r in roots]\n",
    "x = np.zeros(n)\n",
    "r = A.dot(x) - rhs\n",
    "res_cheb = [np.linalg.norm(r)]\n",
    "\n",
    "# Implementation may be non-optimal if number of iterations is not power of two\n",
    "def good_shuffle(idx):\n",
    "    if len(idx) == 1:\n",
    "        return idx\n",
    "    else:\n",
    "        new_len = int(np.ceil((len(idx) / 2)))\n",
    "        new_idx = good_shuffle(idx[:new_len])\n",
    "        res_perm = []\n",
    "        perm_count = 0\n",
    "        for i in new_idx:\n",
    "            res_perm.append(i)\n",
    "            perm_count += 1\n",
    "            if perm_count == len(idx):\n",
    "                break\n",
    "            res_perm.append(len(idx) + 1 - i)\n",
    "            perm_count += 1\n",
    "            if perm_count == len(idx):\n",
    "                break\n",
    "        return res_perm\n",
    "\n",
    "good_perm = good_shuffle([i for i in range(1, niters+1)])\n",
    "# good_perm = [i for i in range(niters, 0, -1)]\n",
    "# good_perm = [i for i in range(niters)]\n",
    "# good_perm = np.random.permutation([i for i in range(1, niters+1)])\n",
    "\n",
    "for i in range(niters):\n",
    "    x = x - 1.0/taus[good_perm[i] - 1] * r\n",
    "    r = A.dot(x) - rhs\n",
    "    res_cheb.append(np.linalg.norm(r))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$||Ax_k - b||_2$')"
      ]
     },
     "execution_count": 82,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEUCAYAAADwYOuyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABZWElEQVR4nO2dd5hVxfnHP+/2ysLSpYtIEREBQUp0bVGCirEEO6CxxBg1iRG7YIw9FlRiMAoo9pifBcHOWhHFiggC0kTp7MIWts/vj3POvWdv2XrL3rvv53nOc++ZM3Nm5ty797vvzDvviDEGRVEURYkECdFugKIoitJ6UNFRFEVRIoaKjqIoihIxVHQURVGUiKGioyiKokSMpGg3oKXToUMH07t37yaVLSkpITMzM7QNauFon1sH2ufWQXP6/MUXX+w0xnT0TVfRqYfevXuzbNmyJpXNz88nLy8vtA1q4WifWwfa59ZBc/osIhsDpevwmqIoihIxVHQURVGUiKGioyiKokQMFR1FURQlYqjoKIqiKBFDRUdRFEWJGCo6iqIoSsRodaIjIpki8oWInBjWila9TvefXg5rFYqiKLFGzIuOiDwhIttF5Duf9BNE5AcRWSsi17ouTQNeCHvD1rxNr40vgO5XpCiK4iHmRQeYC5zgThCRROARYDwwCDhLRAaJyLHA98C2sLeq0yCSq0qgaGvYq1IURYkVYj4MjjHmAxHp7ZM8ElhrjFkHICLPAROBLCATS4j2ichCY0xNWBrWaaD1umMltOkalioURVFijZgXnSB0A35ynW8GRhljLgcQkSnAzmCCIyIXAxcDdO7cmfz8/EY3ILliD2OBtZ+8yuaf4sGgbBjFxcVNel6xjPa5daB9Dg3xKjoSIM0zuWKMmVtXYWPMbGA2wIgRI0xTA96VftWNA8wGDmhFQQI1KGLrQPvcOghHn+P1X/DNQA/XeXfgl0g3YmeHkbDhIyjdHemqFUVRWiTxKjqfA/1EpI+IpABnAq9GsgE1poYJSV8wq00GfDY7klUriqK0WGJedETkWWAJ0F9ENovIhcaYKuBy4E1gJfCCMWZFJNv17sotAPyrXQ58+i8o3hHJ6hVFUVokMT+nY4w5K0j6QmBhhJvj4Y3vXK7SlaWw4Cr43VOQEPM6ryiK0mT0FzBMpCa7fBmOnQ6rFsA7N+tiUUVRWjUxb+m0VFKTXXp++GWwez188hBU7oMT7oJEffSKorQ+9JcvTKQmuSwdERh/N6RkwscPwI4f4JRZ0LZn1NqnKIoSDXR4LUykJfs82oQEOG4GnPIv+OUrmDUGls6G6sroNFBRFCUKqOiEiVRf0XEYejb84WPoNgwW/Q1mjYYVL0NNdUTbpyiKEg1UdMJESlKgoAg27XrD+a/AWc9Z5y9OhkdGwpdPQmVZRNqnKIoSDVR0wkRqUj2PVgT6j4c/LoXT50ByBrz6J7hvALxxHWxfFZmGKoqiRBB1JAgTacl1WDpuEhJh8Klw0G9h/fvwxVz47DH4dBZ0HwkHnwEHnQJZncLZXEVRlIigohMm/BwJ6kME9s+zjpKd8M2z8PUz1rzPG9Og969g8Gkw8CTIyA1HkxVFUcKOik6YSE9JbHrhzA4w5k/WsX0lfPeSdbx2BSz4M/QaA/1/Yw3P5fYJXaMVRVHCjIpOmEhP8Vo6VTVVJCU08VF3GghH3whH3WC5Wq9aAKsWwpvXWUfHgTDgN5YI7TdMw+woitKiUdEJE+mu4bXPt37OlpIt/PaA3/LSmpcY32c8mcmZjbuhiOVm3W0YHHMz7F4HPyyyjo8egA//CRkdoO/RcMCx1mtWx9B2SlEUpZmo6ISJ9GTv8NrFb18MwP45+zNjyQw+3/o5dx1xV/MqyN0fRv/ROkp3w5q3Ye078OO7sPwFK0/XQ6DvMZYI9RgJicnNq1NRFKWZqOiEifQUf++1yhor+sC20m2hrSwjFw6ZZB01NbD1G0uA1r4HHz8IH90HKdnQ5wg44Gjokwft+1rWk6IoSgRR0QkTgSISJIiVVlVTFb6KExJgv0Ot44i/QdkeWP8BrH3XOn543crXppslQn2OtF5zuoWvTYqiKDYqOhGkxtQAXosHoKyqjDnfzeHCgy8kJTEl9JWm5Vhu1gNPsrZV2PWjtR5o/Qew+k3LNRsgty/sbwtQ7yMgs33o26IoSqtHRSdMmAD75lQbK76aW3Tmr5zPrG9mkZ6UzpTBU8LbKBHocIB1HHahNRS3fYUlQOveh29fgGVPWHk7H+wVoZ6HW+KlKIrSTFR0wkQNNX5pZVVWXDX38FqlHWW6uLI4Mg1zk5AAXQ62jtF/tCJe//KV1xL67DFY8jBIAnQeDL3GWmuEeo2x1hIpiqI0EhWdMBHI0nFEp9K1nUGy7VHmtn6iRmKy5eXWY6Q1H1S5D376DDYtgY0fWyF6lv7Lytuhvy1AthDpnJCiKA1ARSdMGPxFZ1/VPqC2wCQnWKJTUV0RmYY1huR0a4ht/yOt86oK2PK1JUAbP7GiJHwxx7rWthf0GkuXsvawq4fl0q3ecYqi+KCiEyYCWTp7yvcAtYfXHOeBFmHp1EdSitcSGvdnaw+gbd9ZArTxY1jzJgNKd8EPD0FWZytf95HQY5S1Zig5Ldo9UBQlyqjohIlAls4/v/gnABU1Fewp30NlTSUpCTEkOr4kJFpi0vUQOPwPYAyfLZzPyM6VsHEJbP4MVr5m5U1Mga5DvaLVYxRkd4lq8xVFiTwqOmHCcY8OREV1BeOeGwfAP8b9w5MW84hQmtkDRuTBiAustOLt1rzQT0utV8c5ASCnp1eAeoy0nBUS9SupKPGM/oWHiUCWjkN5dbnnfYue0wkFWZ1g4InWAda80NZvvUK08RP47r/WtaR0y2rqNswKXtptmM4NKUqcoaITLoJrTi2c6NMxObzWFJJSoPsI6xh9mZW2Z7MlQJuXwc9fwrI5UDXLupaWY0dYGOYVozb7qRApSoyiohMm6hpec+M4HMStpdMQcrpbx+DTrPPqKtix0hKgX760Xj+ZCY4DRlZnrwh1HWqtM8ruokKkKDGAik6YqGt4zY2ziLS0qpQzF5zJ/jn7c/uvbg9n01o+iUneRavDJ1tplftg63deEfrlS1i9yFsmo4O3TJch1muHfpazg6IoLQYVnTDRWEunpLKE1QWrWbFrhYpOIJLTocdh1uFQtge2rYCty615oq3LYemj4FiNSenQeVBtMeo4ANLaRKcPiqKo6ISbysJhJLf9Muh1R5xKKks8ab8U/8JT3z/F1SOuJtH+T33nvp2sLljNmP3GhLfBsURajjcsj0N1JexcbQnQlm8tMVrxshVNwaFNd+g0wBKgTgOt3Vc79ofUrEj3QFFaHSEVHRE5AbgNGATsBJ4DbjXGFPvkOwd40hgTt2MfzvBaVUn/OkXHcSBwi86NH9/I51s/57hexzGs8zAApr4xlQ17N/Dt+d8iOncRnMRk6HyQdRxyppVmjOWssHW5NVe0fZX1uv5DcHkSktPTR4z6Q/sDNNipooSQkImOiAwHFgDbgP8B+wF/BU4VkZOMMStDVVcs4Bleq6n7ETuhcdyiI1iiUlHjdS7YsHcDYIlUWLZAiGdEoG0P6xjwG296TTUUbIDtK11itArW5XuH6AAyO1ni076v/WofuX0gKTXSvVGUmCaUls504EdglDGmEEBExgLPAB+KyARjzNIQ1teicSwdQ91WiScIaIB4bO7AoCkJKVTUVFBeXa6iEyoSEm0h6etdRwSW91zBetjxA+xaax8/wpq34KunXDewxGyI5ELpSGjXB9r1suLQtesFqdkR75KitHRCKTrDgPscwQEwxnwsIiOARcA7InKqMebtENbZYvHGXqtbdBxLx02gyNPJicke0clGf8zCSmKS5fnWoZ//tbK9sPtHS4RsQUra8BV8/SxUFNXOm9HeK0C1XntDTg9rzZKitDJCKTptgELfRGPMDhHJwxp6e01EzgthnY1CRE4BJgCdgEeMMW+Fqy636BzfYRpv7rwrYL4i3x8qXFEKXMNrKQkplFBSK5qBEgXS2ni3A7f5Mj+fvCOPhNLdULgBCjZC4UZr6K5gI2z5BlYugFoLgAWyu1oLXXO6Wc4NOd2sc+d9Vmd1+VbijlCKznpgSKALxphi28ngJazhtndDVamIPAGcCGw3xgx2pZ8APAgkAv8xxtxpjHkZeFlE2gH3AuETHXt4LVGEjokjguYLtHlboOE1x/pR0WmhiFhbfGe2h27D/a/XVEPRFq8QFW60nBv2bLbcvle/Bb5WryRawpTTDdp0875mdbYWwzqvKZkR6aKihIJQik4+cIaI/NUYU+V70RhTJiITsUTndBocKKZe5gIPA086CSKSCDwCHAdsBj4XkVeNMd/bWW60r4cNx5EgLTmJkvLaj+OwLofx+dbPASiu8BedQNsdONGoW3XkglgmIdEbeaH3OP/rxsC+Atj7M+z52Xp1v9/yNax6vba3nUNKNmR3tkTIV5Dcr+ntNGqDEnVCKTpPAl2BEcCngTIYY6pEZBLwT+DQQHkaizHmAxHp7ZM8ElhrjFkHICLPARNFZCVwJ7DIGBPcjzmEpCclUewjOm5HgKLK4MNrtRwJ7DKO44ESZ4hARq51dDk4cB5jrCG84q1QtBWKt/m/bvkaVm8Dlzekh8QUyOxobTWe2dGK4pDZwXteK60jpGSEtctK6yRkomOMWQac0YB8BvhLqOoNQjfgJ9f5ZmAU8CfgWCBHRA4wxjwaqLCIXAxcDNC5c2fy8/Mb3YAfy3603lRXsGHzVmjvvVa02ys0S7f4O/T98ssvAKxcvZL8rVbdFaWWhfPZl59RmFbY6PZEiuLi4iY9r1gmOn1OxFqVsB+kYx2dXFerSkmpKCClooDU8t2e9ykVe0gu30Ny8QZSKr4lubKQxJrA1nN1QioVKTlUJluH874iJYec6lS+fWkZVUnZVCZbR1VSJkhCBPoeHfS7HRriNSJBoDEEY4yZCcysr7AxZjYwG2DEiBEmLy+v0Q3I3pYNb0C77CzSk9rWutatcze+3fBt0LJdunaBNdCjTw/yDrbqfuz1x9i8czMDBg/giO5HNLo9kSI/P5+mPK9YJqb7bAxUlEDpTijZCSU7PK+JpbtIL9lBuidtpXU9WER0SbCG8NJti83z2s513t7nWm7MePHF9OfcRMLR53gVnc1AD9d5d+CXSDbA8V5LS0mieF81JHuvOU4BbjKSMiitKgXgg80fAFYQ0BHzR3D+oPPJTLYmiwPNASlKkxGxwv+kZlmu3PVhDJTvZel7Cxk1pJ813LdvN5Tucr23X/dstsIQle72d5Jwk5Jli1A7rzClt4P0tpDW1ue9fZ7W1orHp3NUMUe8is7nQD8R6QP8DJwJnB3JBjjea+nJSeworAJXWK/URP9V7JnJmR7R2bFvBwCllaWUV5fz2PLHOK7XcUBtF+tPfvmEzUWb+V3/34WrG4pSGxFIy2FfRldrT6SGUrnPEp/SXbWFqbTA53y35eFXVgj7CqnT3ygxxStA9QmU73WNJBE16hUdEakOcxuMMabJ4icizwJ5QAcR2QzcYox5XEQuB97EGvx+whizIiStbSCOpZNue69dO/Ja7vzsTiC46Dhi4+BeOOpYOkWVRTyz8hl+3fvXXPL2JQAqOkrLJzndcvnO6dbwMjU1UL7XFqACS4QcMdpX4H3vXC/aYoUz2rcHyvfUfe+k9AYIVFsr7p59pJTvsoYikzPUwmoGDfmxD/fTbdb9jTFnBUlfCCxszr2bg7NPTnqy5b12zsBzuHfZvVTVVAUMY5OR7O8p5I7H5ni0rdq9ijc3vMnr614PU8sVpYWQkGCLQduGDf25qa6yBMsjTnWJ1h4o3GRFJS8rhCBD2GMAlgAJSbXEKPDR1npNbeN/LSWzVYtWvaJjjIlfd5Qw4rF0UhIpKa/CGENKQgpVNVUBLZ3ctFy/NGe4DbzDdc7wmq9VpCiKi8Qkrwt6Y6mqsITIcxRC2R5++PZz+vfsbKWV762dp2ir931lad33l8SGCVawI8ZFK17ndKKOIxIZycnUGNhXWU1KYgqlVaUkin9ok0Ci43YacESs1P5CV9X4rb9VFCUUJKVAVkfrcLFlZy79x+XVX76qwiVKhT4CFuTYud0lWgHWWLmRRCscU32C5bGy2rjO21ivUQyvpKITJhyRyEhNAqopKqtiWKdhvPfTe6QlpQHQq00vNu7dCED7tPZ+99hbsdd7P1vECssLgdrRCowxuseOorQUklIgyV5k2xSqK63Asg0VrLI9sHNtw0ULrCgWgQTJR6xSykO/l5SKTphY8ssSALJSHNGp5M4j7uTHwh89e+O4LZ5Aczpu0XHC6gQSnaqaqoBu2IqixCCJyd44fk3BEa1yR5T2uoYE3e9dolW81dpx17lmj6SkD70jhB2zUNEJE6+vtyb6M1OTgXL2llVxQFI2gzsMZnPxZgDPVtQA6UnpfvfYXrrd896xnDyi4wqRU15drqKjKIpFc0XLGGteqmwve5d9F9q2AeokECYykizLJSPF0vWiMu8cjOOJliRezQ8kOm6Mz3oFt6WjkacVRQkZIpazQpuumITQ/zOrohMmHBHJSnVExysSjiC5hcSZ5wmGd38e/Mq6BWjFzhWs3NWqdgZXFCWGiPnFoS0VZ44mOcl6fG5Lx7nm9kBLS6xbdJx1P4FwWzpnvn4mAMsnL29kixVFUcJPzC8Obak4lo4kWBF89+7zWiOZSVZ0Afc2BY21dNyU1rcuQFEUpYWgi0PDhDOEViPlJCYIe13Da4Ulls6WVXtF58B2B9Z5v9UFq/3SkiSJKlNFcWUxBWUFtEtrF4qmK4qihA0VlDAx8YCJAAzIHUB2WlKt4bUrnrY2MN3nsnS6ZHapc0hsbeFav7Qqe4PWJb8s4Yjnj2DR+kUhabuiKEq4UNEJE3k98nio10P0atOLtunJfLmpgEuf+oKS8ip+KbDmZ9ok1x2iI1DkgkB8ud3aBNVZG6QoitJSCYnoiEiGiGTVn7N10jYjhe9+3ssbK7by5oqtQBL7Np/Dxf3u9st79YirPe993aiDiZAzN1QZbHMtRVGUFkKTRUdE2onIg/Z2AkXAHhEpFpGPROQ6EekcumbGNrmZ3qjSm3bbsdOKDqa6Mtsv76+6/8rz3te5oE1Km4D3d7zXVHQURWnpNEl0RKQLsAy4HGujdrGPDKwI4LcB60TkThFp9Uvl22W4RGeX19NsV4n/3vTJrsfla+m0SQ0sOtXGcst2RylQFEVpiTTV0vkH0AeoAh4CjgQOAEYAlwGLgTTgb8DHtki1WrrmeC2W1du9O3/u2ecvEu7QOL6WTqB9eADW71kP1LZ0qmvCvbxKURSl8TR1UeYErH1k/2SMme1z7UvgURE5FHgcS4gWisgoY0yr/FfcPby2drt3uwK3R5uDe97G19KprK4kJzWHPUF2RXSLTll1GZkJmU1us6IoSjhoqqXTBsvKmRssgzHmK2AclggdgmUBtUomHdaDK44+gGMHdqas0htZoKisisO7Hl4rb1KCKx5bYm3RKasuIzUh+N7ubtEpDrL7oaIoSjRpquj8CFQaY/wnJVwYY0qx5n0EOK+JdcU8malJ/OXX/enW1jtclpKYQFFZJbOOmcWSs7yuzrVEx8fSqampITWpgaJTWUxFtfXxfLntSzYXbW52PxRFUZpLU0XnOSBdRA6vL6Mx5lOgGBjQxLrihrYuh4IeuekUlVn74GSleL3N3cNrvnM6VSbwVtcObkeCl9a8xPD5w1lXuI7Jb0xm/P/Gh6ILiqIozaJe0bHdn8eLSFdX8n3AamCmiNQZNExEErHmjhqwnV180z7LKzo9czNqRZ52cFs6vqJTbarrFh2XpbPgxwUArNi1ImDefVX7PBvDKYqiRIqGWDr/ABYAm0Vkm4i8CUzHms85EHhfRA6qo/wZWJ5sbzavqbFPlzZeEemQlRrYkaCOjd1qamo4pucxtdK6Znr/F3BHrXbiugVau1NdU83Ip0dy+9LbG9kDRVGU5tEQ0fkfsM5+3xE4DrgauB3IxvJO+1pEForIRSIyTER6i8hIEZmB5cH2E3Bt6JsfW3RxuU5npyUHFB33xm5ZybWDPFSZKn5/8O85d+C53vwuy6iyppKUBMuachaMuoXIwVnX88IPLzSlG4qiKE2mIVGmTweww9wcYh9D7WMwlhWTCJwAHB/gFruxRCp4bP5WQt+OlogM7NqG7LQkisurqK4xJCZ4d3cQ8b6/8OAL2bFvBy+vfRmwLBQRqWXdiGtniMrqSjKTM6kor/AMnTnOBG6cDeB8dyNVFEUJNw1ep2OMKQY+tg8ARCQBy0HALUSHAJ1cRdtjOR6YxtQXj2SmJvHeX48kKzWJ177dAkBxWRU5GYGDNmQmZ/L3sX+nb05f/vnFPz0WSnJi4Pzb9233Sws0vFbX3jyKoijhpFkiYIypAb63j2eddNvpwC1EhwJ9m1NXvLC/be1kp1mPfm9ZZVDRcRjQ3nL8c0THGUJrCIGG19TCURQlWoTF8jDGbAG2AG84aSKSHrxE6yM71Xr0JRX+ouCLe54HaofDqU9A1NJRFKUlEbH9dIwx+yJVVyyQaYtOcQBnAl98h9PcbtP1xVira05HURQl0ugmblEiyx5eKypvgKWTUNvSOarHURzc4WAAKmrqDArhGZJzo+tzFEWJFq16Yj+aZNdh6Tz9m6cpqfSupfUdXktOTGbGmBmc+uqpHtfoYAQSGLV0FEWJFmrpRAmPpRNAdIZ0HMLo/UZ7zpMT/B0NnGgFgYbP3LhF5/V1r3P0C0frvjuKokQNtXSiRJbjSNCA4bWAopMYXHQSJMEjNvNXzvekX//R9dSYGnbu29mkNiuKojQXtXSiRGZKw+d0Am3e5lg6gYbKMpMD76PjCFF9Q3KKoijholWJjohkisg8EXlMRM6JZlsSEoSs1KQGea8FEhHfYKD15XejoqMoSrQImeiIyN9F5B+hul8j6n1CRLaLyHc+6SeIyA8islZEnLhvpwL/NcZcBJwc6bb6kpWaFDDStC8ZyRl+aYGG3Dz39YnZ5ot7SG51wWre/+l9AD75+ROd71EUJayE0tK5lugE9ZyLFffNg72dwiPAeGAQcJaIDAK6YwUfBah7gUsEyLLjr9VHggT/mA7peIhf2gFtD6jzfk4EaoDTXj2Ny9+7nBW7VnDJO5dw77J7622PoihKU4n54TVjzAdYQUXdjATWGmPW2bubPgdMBDZjCQ+0gL5npTZMdILx9ulvM/u42Z7zzhmdAeiY0bHOcmVVZX5pzvbWawrXNLk9iqIo9RGv3mvd8Fo0YInNKGAm8LCITABeC1ZYRC4GLgbo3Lkz+fn5TWpEcXFxnWUrS/fxy14adf+68vaRPmxjGz9v/rnOe3y94mu/tOXfLAdgZ8HOJvcX6u9zPKJ9bh1on0NDvIqOBEgzxpgSYGp9hY0xs4HZACNGjDB5eXlNakR+fj51lX32p2Ws31lCXt6R9d7rz9/9mc4ZncnbP8D95lkvvbr34tMfPqVnj54QeMNQAHr27elnGx4y9BB4EzKzMj1tLqsqY86KOfx+8O+DRrb2pb4+xyPa59aB9jk0xKvobAZ6uM67A79EqS1ByUxNoqS8YVNLFwy+oN48jnOBe0+eQOyr8g+D59l/xxVW58nvn2TW17PISMpg8kGTG9RORVGUuoj6vEaY+BzoJyJ9RCQFOBN4Ncpt8iO7gd5r9TG4/WDAG6NNAhp6XgKJjuNG7Y5K7WyLUFxZ3Ow2KoqiQByIjog8CywB+ovIZhG50BhTBVwOvAmsBF4wxtQx4BQdMlOTKKmobvZWA7N/PZvnT3yeREkEAnu7dcrw7qsXSHRKq0qB2vvvOItS6wu1oyiK0lBifnjNGHNWkPSFwMIIN6dRZKYmUV1jKK+qIS05scn3yU7JZlD7Qby76V3AsnQ6pXeqtZNodnI227HOA3mv7au0hMgtMM5wnYqOoiihIuYtnVjG2T20OW7TbtzRqF855ZVa1zJTvFEKHKvGjWP9uIfXHNEJtBGcoihKU1DRiSJZjdjIrSEkJljWksGQlZJFblqu51p6knfj1kXrF/mVDSQ6zvCaio6iKKFCRSeKeHYPDZGl48zpOBu3uQOFZiT5h9Jx88CXDwDeYTZ3eY3VpihKqIj5OZ1YxtnILdCeOk3BER3HMSElwSs6bkunLqpMFdd9eB1bS7Yyqf8kAI3HpihKyFBLJ4pkhXhOx1mf41g64/uM91wLFDQ0GAvWLWDZtmUeL7j6tsRWFEVpKKEUnboXhyh+OMNrpRXhsXQuG3oZQzsOBQJbOnVFqgbvXj3lVTq8pihKaAjl8NqtqPA0Cmcjt4ZGJagPX0snQRJon94eCDynU5+DgCM6JVUlvPDDC3RM78hRPY8KSVsVRWmdhEx0jDG3huperYWMVMsyCbWl44S0AW+UgobO6bhxLKbSylL+/unfAVg+eXlzm6koSitG53SiSIa9IDRUlo4jMIHW2gTa8vqMA8+o836OeJVW+q/rURRFaQoqOlEkKTGB1KSEkFk6zhCa2+25rm2th3QcUuf9nJA4JVUlnrTSylLe3fhuc5qpKEorRl2mo4wVfy00ouMMobkjDjhp7iE3T93JmX5pbpwFoyWVXtG59dNbeX3d67x08ksc2O7AZrc5XikvL2f37t0UFRVRXR31TWpDTk5ODitXrox2MyJKa+5zYmIi2dnZ5Obmkpqa2qx7quhEmczUxJBFJHDcot0BPdMS0/zSHOrzXnPKuIOAbineAsCe8j3Na2wcU15ezqZNm2jXrh29e/cmOTm53u0mYo2ioiKys7Oj3YyI0lr7nJWVRWVlJXv37mXTpk307NmzWcKjw2tRJjMlieIQzel4LB3XHIwjRIGCfDZlCwRPaBxdMBqU3bt3065dOzp06EBKSkrcCY7SuhARUlJS6NChA+3atWP37t31F6qDsIiOvYeN0gCy05IoCdHiUM+cTgBLp6w6gOgE+DHsltXN8z6QA4FnuwNdMBqUoqIi2rRpE+1mKErIadOmDUVFRc26R7gsnf+E6b5xR1ZqUsgiEvTO6c2J+5/IXUfc5Uk7qMNB1mv7g/zyC8KnZ39aKy07xTt8UNeQnNtDbtPeTbyy9hW/vK2V6upqkpMbtr23osQSycnJzZ6jDJfonCoi19SVQUQmhqnumCIrLTlklk5SQhJ3/OoO+uf296Qd2ulQ3j79bU7c/0S//AmS4OdM4F5EGnB4LcF/Y7ezF57NjR/f2OzN6OIJHVJT4pFQfK/DJToXAv8QkQmBLorIdOClMNUdU2SlJrE3RI4EweiS2SXglyXQnI57PU8g0XHWArlFx3Eq0C0QFEWpj7B4rxljnheRQ4BnRGS0MeZ7ABHJBJ4CJgL/CEfdsUZWamLILJ36mHP8HBITEjl/0flA4P9anDkggI1FG/2uO2V8F6BW1lQGnDdSFEVxE06X6RuAIcBrIjICyAVeBnoBpxljXg5j3TFDVmoy+yqrqa4xJCaEd0hmRJcRtc4DiY7b0llTsCbovXw3e6usqdTAoIqi1EtIhtdE5AYROVlE+jhpxhrgPwsoB94GPgNSgFEqOF5Cvb1BYwg0vJaa6O9/747b9uqPrwK1XbCdeR7d7E1RlPoI1ZzOLcD/AWtFZK+ILBGR2cBkYCZwMPAJcJgxpnUt6a2HLDvoZyRFx9knx3l1EyhGW6DIBaVVpcz8ciafb/2c5ETLU0tFRwnEqlWr+NOf/sTgwYPJyckhJSWF/fbbjwkTJvD4449TVuY/LLt69Wr+8pe/MGzYMHJzc0lOTiY3N5dRo0Zx9dVX88UXX0ShJ0ooCJXoZAHDgCnAo8AeYAKW4DwCJAOHAs+JyJ0icraIDA5R3TFNVqr1gx2qqAQNIcHnY3/9t6973tdn6TiUVpby2PLHuODNCzxlVHQUX2699VYOOuggHn74YbKzs5k8eTJXX30148ePZ9WqVfz+979n3LhxnvzGGGbMmMHAgQO5//77EREmTZrENddcw7nnnkt6ejoPPfQQI0aM4JFHHoliz5SmEpI5HWNMBfCNfXgQkVyseZ2DXa9/BDIBAySGov5Yxju8FjnPrwRJAOO1dHq26em5Fkh0spKz/NLc8dh0eE0JxO23384tt9xCjx49ePHFFxk1apRfngULFvDPf/7Tc37rrbcyffp0evTowbPPPsvYsWP9ymzfvp0HHniAPXs0FFMs0iTREZGOxpgd9eUzxuwG8u3DKSvA/lgC1OrJsncPLYqkpWOLjXtOx/FACzS81jmzMyt31x4VrSU6iSo6Sm02bNjA9OnTSU5OZuHChQweHHhg48QTT+S4444DYN26ddx2222kpKSwaNEiDjrIf0EzQKdOnbj99tupqor8PKjSfJo6vPapiPRrSkFj8aM6E1hkR8GRINCcjmPhBBKdLhld/NLc2x04ZXXfHcVhzpw5VFZWctpppwUVHAcneOScOXOoqqri9NNPDyo4bpKSNF5xLNLUT60P8LGITDTGLAllg1objuhE0tJxdhh1k5qYSnFlccD8Oak5fmlFFd74S1kpWZ60HPzzKl5mvLaC73/ZG+1m1Mmg/dpwy0n1/+jXxUcffQTAMccc0+AyH3/8MQBHH310s+pWWjZNFZ35wLnAuyJyrjHmf3VltiMT3GqMGd7E+uIWZ3gtoo4ECfbwmmudzn9+/R9eWvMSnTM6++UPtBHc3nLvD2d2shWvbW/FXhUdBYAtW6wtMLp3797gMlu3bgWgW7duftc2bNjA3Llza6W1bduWq666qsltVKJDk0THGHO+iGwAbgReEJG/GmMe9M0nInlYkQcOb0Yb45rMFMfSiZwjgWPpuL3YDmh3ANNGTmPR+kUA9G7Tmw17NwCBvdcKyws9753tE/ZW7KXaVGOMYf2e9ewq28VhXQ4LUy9ik+ZaELGCE4evMbG66iqzYcMGZsyYUSutV69eKjoxSJNdpo0xN2PFWKsG7hOR+5xrIjJKRN4B3gVGY3mqvdDMtsYlCQlCVmoSRRGc03FEJNAftxNF2j0E5w4C6uAWHWdu6Oein7lq01U88d0TTHxlIhe8eUEom63EEPvttx8AmzdvbnCZrl27AvDzzz/7XcvLy8MYgzGGykqN8RfLNGudjjFmDnAiUAxcKSIvi8irWAtBj8YSm+eAg40xZzW3sfFKRkoipSHayK0heEQnQESCQENpgdLcOFthby62fmD+u/q/zW2iEuM4a2/efffdBpdx3KMbU0aJPZq9ONQY8zbW/I4AJ2EtCnXEZrAx5myNQlA3WalJlFREztJxhsMCRYV2og9UG68IuoOA1oWzm6hGm1amTp1KcnIyL730Et9//32decvLLVf7KVOmkJSUxH//+19WrtSfjHilWaIjIr3tcDfOv7ZiH18DlxljVjWvea2DjAhGmgbvcFkgF2fnmlt0hnQcUuf9HEvH2QqhqkbXT7R2evfuzfTp06moqGDChAksW7YsYL433niD8ePHA9C3b19uvPFGKioqGD9+PJ988knAMoWFheFqthIBmro4tA9WFOnz7HsIsBh4GLgTK+TNxyLyG2OMf3x8pRaZKUmUVERueM0jOlXBRcdtrbRPb8/X533N0KeGBryfwZoAdtyo3WWNMZ65o+KKYkQkYCw3Jf64/vrrqaqqYsaMGRx22GGMGTOGESNGkJWVxbZt2/jggw9Ys2YNI0Z4o5/ffPPNGGP4+9//ztixYxk+fDgjR44kNzeXwsJCNmzYwDvvvAPAEUccEa2uKc2gqS7Tq/CKzafADcaYxQAi8j7wKjAGWCIiJxljNDpfHWSmJrG9KHJ70Rzb61je++k9erXp5XfNGXpzR5EGSEwIHrFoS7HlHltU6S867igHo58dDcDyycub0Xollrj55ps544wzmDVrFosXL2bOnDmUlZXRvn17hg4dyrRp0zj33HM9+UWE6dOnc9ZZZ/Hoo4+yePFinnnmGUpKSsjOzqZv37784Q9/4LzzzmPYsGFR7JnSVJoqOslYQ2g3GWNed18wxuwWkWOw1vKcBrwvImcZY15rVktDhIicgjXv1Al4xBjzVnRbZDkSlETQkeCkvieR1yOP7JRsv2uOk4Gv6NTF0q1LgcDDa+XV5QGjHCith4EDB/LQQw81qkz//v25//77w9QiJZo0dU5nkjFmmK/gOBhjyo0xZwAPABnA/0Tk8ibW5UFEnhCR7SLynU/6CSLyg4isFZFr67qHMeZlY8xFWBGxJzW3TaEgKzUp4vvpBBIcsETnmJ7HMPPomX7Xnhr/lOe941odCLelo/HYFEVx09TFoS82MN9f7EWk92EJ0MNNqc/FXPseTzoJIpKItX3CccBm4HPbbTsRuMOn/AXGmO32+xvtclEnKzUpoo4EdSEiPHDUAwGvtUtr53mflpjWIC+1iuqKUDVNUZQ4IOwR84wxM0XkJ6zhtube6wMR6e2TPBJYa4xZByAizwETjTF3YK0hqoUd5fpOYJEx5stA9YjIxcDFAJ07dyY/P79J7S0uLm5Q2R1bKiitqOa9xYtJaMQK7kjh9GFn5U5PmtQ0rJ0fLPmALsm1A4Y29Xm2VHw/55ycHIqKioIXiAOqq6vjvo++aJ8tysrKmvU3HJEwrcaY/xORcEXx6wb85DrfDPhv3OHlT8CxQI6IHGCMedQ3gzFmNjAbYMSIESYvL69JDcvPz6chZdckrOOVH1dy2OhxZKcFH7aKOPOsF6cPW0u2epzjs9Oz2Vtcf+DKQUMHeV2ufe4XL/h+zitXriQ7O/DwZbxQVFQU9330RftskZaWxqGHHtrke4ZVdGyrZIjr+F04qgmQZoJlNsbMxNrRtMWQaQf9LCmvblmi40NSgvfr0tAFo3vK97Bz3046pHcIV7MURYkhQiI6IpKJtSnbIXgF5mAgEv8WbAZ6uM67A79EoN6QUXv30Ib9mEcDdzy2hnqkffjzh1z27mXc8Svf6TVFUVojjRYdEemLJSpugemNv8XhPt8NrG5aE+vlc6CfvWD1Z+BM4Oww1RUWsp3tDSLoNt0U3Gt16ovH5vDlNmvabOmWpZ60GlNTawM5RVFaDw0SHRF5FEtkBmO5QHsu+WTdBPTEGt76A/AtsNoYs6v5TQUReRbIAzqIyGbgFmPM47Y79ptYHmtPGGNWhKK+SJHl2citZcUsu27kdby7yRt8MUm8Xxdnt1DPtYSkgOFvHJdptxdbWVWZZxGqoiiti4ZaOhdjCYkjMj8B3wMrXMf3xphiEamx8zxrjAnpFonBIlUbYxYCC0NZVyTxbFkdwY3cGsLZA8/m7IFeo7HWnI6PpZOWmEZxjf/Oo44Qud2riyqKVHQUpZXS2OG12cA0Y8yecDSmteLsHhrJPXWagntOx9eRINhwmbPdgROBGmBPxR4SJIGOGR3D0EpFUVoyDR1Y/ytQCFwEfCAi48PWolZIdqrlsVbUwiwdX9zC4ju8ZowhNy03aNmKGu/w2tzv5nL0i0ezctdKtpZspbjC30JSFCU+aZDoGGPuBw7AigYwAFggIu+JyIi6SyoNITPVsiBaSlSCYLh3GvUdXquhpk6PNveczuKfFgOwumA1x/33OE5/7fQQt1RRlJZKg12IjDEFxpgrsZwJFmBN6C8VkedF5IAwta9VkJSYQFpyQsTjrzUHX0unxtT4pblxz+mUVZfVSvu52H97YkVR4pNG+60aY9YYYyZibUf9DXAGsEJEHhaRTqFuYGshGkE/m8PJfU8G4JIhlwDW8Fp9lo4Twbq6ptqTpihu8vPzPdsbNJW8vLxaVnlLY+7cuYgIc+fOjXZTokKTF0sYY/KB4cCFwE7gMmBtaJrV+shsQUE/G0L/3P4sn7ycgbkDAWun0dSEui0dR3ScTd90W+v4RkRqHYmJieTm5pKXl8fcuXMxJmjgECWOaVZEAmN9a+aIyPPANCyHA4c/isiDxhj/7SkVPzJTYkt0HBznAmMM1426jnMWnuO51im9E9v3WUG93Zu5OajotA5uueUWACorK1m7di3/93//x/vvv8+yZct4+OHagedHjhzJypUr6dBBwybFKyEJg2MLyy0i8m+s7QTOBW4DrhKRe4BZKj51EyvDa1ePuJo95V6PeSdKQQ01DOk4hAk5E3h9z+u1rkHgoTQdXmsd+A6VffzxxxxxxBHMmjWLv/71r/Tp08dzLSMjgwEDBkS4hUokCWksEmPML8aYycBhwIdAR+AuYL2IXB3KuuKNzNTI7h7aVCYfNJkrhl3hOXeiFNSYmlrnAOIKWFFRXeHJ405TWh9jx45lwIABGGP44ovaO9nXNaeze/dubrjhBgYPHkxGRgY5OTkccsghXHvttZSUlPjlr6qq4vbbb6dfv36kpqbSo0cPpk2bRkWF//fu5Zdf5txzz+XAAw8kMzOTrKwshg8fzsyZM6mpqfHLP2XKFESEdevW8dBDDzFkyBDS09NrRRtfu3YtZ5xxBu3atSMzM5MxY8bw+usB970E4Ntvv+Wss86id+/epKam0rFjR4YNG8ZVV11FZWXtUYE9e/Zw3XXX0b9/f9LS0mjXrh3HH38877zzjt993c/066+/ZsKECbRt25aMjAyOPPJIPvnkk6BtCgdhiTJt71OTJyK/xRKdA+zXe8NRXzyQlZbMxl2xZwy6rRmoLTo1eP9YK2oqSJf0Wnl1eK314sznJCc3LKr6+vXrOeqoo9i4cSPDhw/nD3/4AzU1NaxevZr777+fSy+9lMzMzFplzj77bD788EPGjx9PmzZtWLhwIXfffTfbt29nzpw5tfJee+21JCQkMGrUKLp168aePXt47733uPLKK/n888956qmnCMSVV17Jhx9+yIQJE/jNb35DYqL197BmzRpGjx7Nrl27GD9+PEOHDmXt2rWccsopjB/vv8zx22+/ZdSoUYgIJ598Mn369GHv3r2sXbuWWbNmcdttt3meVWFhIWPHjuX777/nsMMO46qrrmLnzp288MIL/PrXv+Zf//oXl1xyiV8dy5Yt4+6772b06NH8/ve/Z9OmTbz00kscc8wxfP311/Tv379Bn0VzCevWBvY+Oq9h7WFzYzjrinWyUhNbfESCQPhGIqglOi7LpqiiiCJqbwbla/m0ChZdC1uXR7sVddPlYBh/Z9hu/8EHH/DDDz+QkpLCyJEjG1Tm3HPPZePGjdx+++1cd911ta7t3LmTrKwsvzI//vgjK1asIDfXWrT8j3/8g0MOOYQnn3ySO+64gy5dvJsLvv766/Tt27dW+ZqaGqZOncqTTz7J5ZdfzqBBg/zq+PLLL/nqq69qDREC/PGPf2TXrl088MADXHnllZ70V155hVNOOcXvPvPmzaOsrIyXX36ZiRMn1rpWUFBARoY3bNS0adP4/vvvufjii3n00Uc9nnrTpk1jxIgRXHHFFRx//PH07t271n1ef/115syZw5QpUzxp//73v7n00kt58MEHmTVrll+7wkHYQ/0aY6rsxaV9683ciolVRwJ3PDaoLTrBdzWyaJWi0wqZPn0606dP54YbbmDSpEkce+yxGGO499576dq1a73lv/jiCz755BOGDh3KtGnT/K536NCBtDT/qOd33XWXR3AAMjMzOeecc6ipqWHZsmW18voKDkBCQoJHMN58882Abbvmmmv8BGfz5s28/fbb9OnTh8svv7zWtYkTJ3LkkUcG6Smkp6f7pbVr146EBOunurKykvnz55OVlcUdd9xRyzW8X79+XHHFFVRUVPDkk0/63Wfs2LG1BAfgggsuICkpic8++yxom0JNRHYOBTDGFEaqrlgkMzWJ0opqamoMCQktd42BL+54bBB8eC0Qm4o2ed5X11STIAkten1FSAijBdFSmTFjRq1zEeHxxx9n6tSpDSr/6aefAnD88cd7fnwbwogR/gFTevSwtt4qKCiolb5r1y7uueceFi5cyLp16/zmiH7+OfAC5kCW2ldffQXAuHHjPMNtbvLy8nj//fdrpU2aNIkHH3yQU045hdNPP51jjz2WsWPH+onhqlWrKC0tZezYsbUE1eHoo4/mtttu87TBTaDnkZycTOfOnf2eRzip9xMUkeowH7H3730Y8ESaroitx+ErOgenH0zHdCuQZ6CtDpy1OgAf/fyR5/19X9zHkCeHsK9qX5haqkQLYwzGGIqLi3n77bfp0aMHl156Ke+9916DyhcWFgLQrVu3RtXbtm1bv7SkJOvvrLra67RTWFjIYYcdxl133UV6ejrnn38+N9xwA7fccovH0ikvLw9Yh3uIzmHPHsu7s3Pnzg0uM3LkSD788EOOPvpo/vvf/zJ58mQOOOAABgwYwLPPPut372AWopPuPDM3gZ4HWM/E/TzCTUP+bZAwH7qbF94tq0tjwIPNja8jQWpCKk8c/wQQ2DstK9l/7B3gye+t4QC3O7YSX2RmZnLsscfy2muvUV1dzeTJkyktrd95xvmxDGZtNJf//Oc/rF+/nltuuYWlS5d6Ju6nT5/OpEmT6iwbyDLPyckBYNu2bQHLbN26NWD66NGjWbBgAQUFBXz88cfcdNNNbNu2jbPPPtvjlebcO9g9tmzZUitfS6TeH3xjTEK4j0h0tKWT6dk9NLYtHfAGAw0oOimBRcdB3ajjnyFDhnDRRRexefNm7r///nrzH3744YA1rxLIfbm5rF1rBVI57bTT/K75DoM1hEMPPRSAjz76KKAFkZ+fX2f51NRUxowZw6233srMmTMBywEBoH///mRkZPD1118HHBJbvNgKpjts2LBGtztS6A9+CyEzJTYiTfvia+mAd6+dKuPfF2foLRjqRt06uPHGG0lLS+Pee++tdz5h+PDhjBkzhq+//pq77rrL7/quXbsoKytrclscLy9fMfjqq6+44447Gn2/7t27c9xxx7F+/Xq/iAuvvPJKQCH78MMPPUNnbhxryfFeS0lJ4ZxzzqG4uJibb765Vt4ff/yRmTNnkpyczHnnndfodkeKiDkSKHXjWDoxJzp1WDqB6JPTh8+2BveUUUunddCtWzcuueQSHnzwQe6+++56f9znz59PXl4e119/PS+99BJ5eXkYY1izZg1vvfUWq1at8nMRbijnn38+99xzD1dddRWLFy+mX79+rFmzhgULFnDqqafy/PPPN/qejzzyCKNHj+aqq67irbfe4pBDDvGEADrppJN47bXXauX/5z//yVtvvUVeXh77778/WVlZrFixgkWLFtGuXTsuvvhiT94777yTDz/8kIcffpjPP/+co446yrNOp6ioiIcfftjPo64loZZOCyErRofXfF2mwX/bg4ZeAyivDjxhq8Qf1113HRkZGcycOTPo/IdDnz59+PLLL7nmmms8P6yPP/44mzZt4q9//SudOjU9wP1+++3nWeD50Ucf8fDDD7Nx40ZmzZrFnXc2zduwX79+fPrpp5x22ml8/PHHPPjgg/z000+8/PLLnHrqqX75L7vsMiZNmsSGDRt46qmneOihh1i9ejWXXXYZX331VS0vttzcXJYsWcI111zDrl27uO+++3jxxRcZOXIkb7zxBpdddlmTn0UkEI30WjcjRowwvj79DSU/P79WWIy6WL+zhKPuzef+SYfw20O7N6m+aLC7bDdHPm+tO1g+ebmnzwfPO5jctFx2l+2ulX/KQVOYu2Ju0Ps99uvHOLzr4eFscsjx/ZxXrlzJwIEDo9egCFBUVER2dna0mxFRtM8WDf1+i8gXxhg/P221dFoIzu6hxTHmvZaRlBEw/T+//g/Pn+g/LFHfOpzyKrV0FCWeqXdOR0TC/StojDGtfm4pK0bndIINl43qOqrW+YT9J/D6utdrBQENhLOrqKIo8UlDfuzDvUQ8zpegN4z05EQSJPZEp6ERBByPNt9Ybb6UVanoKEo8U6/o6DqayCAiZKbExp46TcHZwK0+S0cdCRQlvlFBaUHE2pbVjcEZhgtkGXXK8HoeqaWjKPGNik4LIistNi2d+/PuZ87xc+rM47Z03j3j3VrX2qS08bwvLC8MefsURWk5qOi0IKwtq2PLew3g2F7HMqKLfwRbNykJlujUmBo6ZXSia6Y3YGFmsnfzrceWP8Ynv0R2J0NFUSKHik4LIjstieKy+AoD8+/j/s3EvhM94XKc0DjJCd4dI33drpf8siRyDVQUJaKo6LQgslKTKCqLveG1uhiz3xhuG3ebZ58dJ2DjeYO8saEykmuLTqAtERRFiQ9UdFoQ2TE6p9MQfC2dMwecyZj9xgD+lo4G/VSU+EVFpwWRlZpMcZxZOg5OYFC3FeOs2VFLR1FaDyo6LYis1ESKK6qoqYm/eHjJidYcjltQHCHyFR21dBQlflHRaUFkpiZhDOyrjD0PtvpwIhK4t6P2WDpJaukoSmuh1YmOiGSKyBcicmK02+JLrO6p0xAcYXGLjsfS8RGdahN/oqvUz4YNGxARpkyZEtZ6RKTB0d/DTUtqS6SIGdERkSdEZLuIfOeTfoKI/CAia0Xk2gbcahrwQnha2Tyy02JzT52GkJ6UDgSxdHyG1wJtDKfELqtWreJPf/oTgwcPJicnh5SUFPbbbz8mTJjA448/3qxdP5XYI5aiO88FHgaedBJEJBF4BDgO2Ax8LiKvAomA71aEFwBDgO+B4FtbRpHMlDgWnWR/0XE82nwjVavoxA+33norM2bMoKamhsMPP5zJkyeTlZXFtm3byM/P5/e//z3/+te/aOqeVUrsETOiY4z5QER6+ySPBNYaY9YBiMhzwERjzB2A3/CZiBwFZAKDgH0istAYUxMg38XAxQCdO3f22zu9oRQXFzeq7Jpd1rDSx0u/YPfa2PzhDdbnzRWbAdhRuMNzvXBXIQDffv9trbzbt21v8jOPBr59zsnJoaioKHoNigDV1dX19vHee+/l1ltvpXv37sybN4/DDjvML8+iRYt4+OGHKSoqori4GIDKysqwP7+GtD8UZaJ531AQqG1lZWXN+/s0xsTMAfQGvnOdnw78x3V+HvBwA+4zBTixIXUOHz7cNJXFixc3Kv+3PxWaXtMWmLdWbG1yndEmWJ937dtlBs8dbO757B5P2u2f3m4Gzx1snlrxlBk8d7DnuPGjGyPU2tDg2+fvv/8+Og2JIHv37q3z+vr1601ycrJJTk42y5cvrzNvWVmZpwxgJk+ebNavX28mTZpk2rdvb1JTU83w4cPNa6+9FvQezzzzjMnLyzNt27Y1qampZsCAAebvf/+7595uAHPkkUean3/+2Zx77rmmY8eOJi0tzQwbNsw8/fTTtfIuWrTIAGbq1KkB+1xWVmbat29v2rdv76mrvLzcPPjgg+bQQw81bdu2Nenp6aZXr17m5JNPNm+//XbAtuzYscNcdNFFpkuXLiYlJcUMGjTIPPHEE0H7+8Ybb5jx48eb9u3bm5SUFLP//vubq6++2hQUFHjy7Nu3z+Tk5JiOHTuaysrKgPe55JJLDGAWLFgQ8HqgPjf0+w0sMwF+U2NmTicIgeLk1+tvbIyZa4xZEIb2NAtn99B4dCTITcvlndPf4arhV3nS0pKsUU7fjdtq/I1PJcaYM2cOlZWVnHbaaQwePLjOvKmptYdXN27cyMiRI9mwYQPnnXcekyZN4rvvvmPixIksXrzYr/yFF17I2Wefzdq1azn11FP54x//SG5uLjfddBMnnHACVVX+f08FBQWMGTOG5cuXM3XqVM4//3zWrVvHOeecwz333OPJd/zxx9O3b1+ef/559uzZ43efl156iV27djFlyhRPP6ZMmcKVV15JZWUl559/PldccQVHHHEEy5cv54033vC7R2FhIWPHjmXJkiWcfvrpnH/++fzyyy9ccMEFzJs3zy//rbfeygknnMDSpUuZMGECV1xxBQcccAD33nsvY8eOZe/evQCkpaUxadIkduzYwaJFi/zuU15ezgsvvEDnzp05/vjj/a6Hi5gZXgvCZqCH67w78EuU2tJssuLYkQCgc2bnWueO6LjneSC+vdfu+uwuVu1eFe1m1MmA3AFMGzmtWff46KOPADjmmGMaXTY/P5/p06dzyy23eNLOPvtsTjjhBO655x6OOuooT/rcuXN54okn+O1vf8vTTz9Nenq659r06dOZMWMGjzzyCFdeeWWtOr799lvOOOMMnnvuORISrP+9r732WoYPH84NN9zAaaedxv7774+IcOmll/K3v/2N5557jquvvrrWfWbPng3AxRdfDMCePXt47rnnGD58OEuXLiUxsfYw+a5du/z6+80333DhhRfy73//25P/z3/+M0OGDOGuu+5i8uTJnryLFy/mlltuYfTo0SxcuJC2bdvWehZTp07llltu4f777wcsAZw9ezbz5s3jpJNOqlXvq6++SkFBAX/5y19ISoqcFMS6pfM50E9E+ohICnAm8GqU29RknC2r41V0fAnkRt2rTS+qa+JXdFoLW7ZsAaB79+6NLturVy9uvPHGWmnHH388PXv25LPPPquV/uCDD5KUlMQTTzxRS3AAbrrpJtq3b8/TTz/tV0diYiJ33XWXR3AA+vTpwxVXXEFlZSVPPfWUJ33q1KmkpaUxZ07t7Tt++OEH3n//fY466igOPPBAwHKBNsaQmppa694O7du390vLyMjgvvvuqyVQgwYNYuzYsaxcubLWnMrMmTMBeOyxx2oJDlgCM3To0Fr9HT16NAceeCCvvfYau3fvrpXfsaLcohYJYsbSEZFngTygg4hsBm4xxjwuIpcDb2J5rD1hjFkRxWY2i1jdsrqptE1tC1jbHnTL6ubxYotnS6e5FkSsYA3pN3w7czdDhw71sxAAevTowZIl3gjkpaWlfPPNN3To0IEHHngg4L1SU1NZuXKlX3rPnj3p06ePX3peXh4zZszgq6++8qS1b9+e3/3udzz55JN88sknjBljxQx0rJxLL73Uk7dNmzacdNJJvPbaawwdOpTTTjuNX/3qV4waNYqMjNpLAxz69etHmzZt/NJ79LAGcQoLC8nOzgZgyZIlJCcn8+KLL/Liiy/6lamoqGDHjh3s2rXLI3CTJ0/mhhtu4LnnnuOyyy4DYNu2bbz55psceuihDBkyJGC7wkXMiI4x5qwg6QuBhRFuTlgQETLjMNJ0ME7c/0QKygqYNGCSZ67nt6/8Vud04oD99tuPVatWsXnz5kaX9f0P3iEpKckTpRyseRljDDt27GDGjBmNqqNz584B07t06QLgN39z2WWX8eSTT/Lvf/+bMWPGUF5ezrx58+jUqROnnHJKrbzPP/88d911F88884xniDAtLY3TTz+de++916/uuvoLlgeZw65du6iqqqq3v8XFxR7ROf/887npppuYN2+eR3SefvppqqqqIm7lQOwPr8UdWXG8ZbUviQmJTBk8xbNwFKw1OvFs6bQWxo0bB8C7775bT86mk5OTA8Chhx7aEM/XWmzbti3gPbdu3Vrr3g6jRo1i6NChvPDCCxQUFHgcCKZOnUpKSkqtvOnp6UyfPp3Vq1ezadMm5s+fz7hx45g/fz6nn356s/vcrl27evvbq1cvT5nu3btz9NFH89lnn7FqlTWfOG/ePJKTkzn77LOb1Z6moKLTwshMTaKkonWITiASJEEtnThg6tSpJCcn89JLL/H999/Xmbe8vLxJdWRlZXHQQQexYsUKv/mK+ti0aRMbNmzwS3fWnxx66KF+1y688ELKysp48sknmT17NiLCRRddVGc9PXr04JxzzuHNN9+kX79+fPTRRwGdCRrK4YcfTkFBAStWNG4WwQktNG/ePL7++mu+/fZbxo8fT8eOHZvclqaiotPCyIzRLatDhVo68UHv3r2ZPn06FRUVTJgwIWjEgTfeeIPx48c3uZ6//OUvVFRUcMEFF1BYWOh3vaCggC+//NIvvbq6mmnTptUarlu/fj0zZ84kKSmJc88916/MGWecQU5ODnfffTfvv/8+xx13HH379q2VZ8eOHSxdutSvbElJCUVFRSQlJflZRo3hz3/+MwAXXXQRv/zi76hbUlLCp59+6pd+6qmn0qZNG+bPn8/cuXMBwh7jLhgxM6fTWshKTWw1w2uBSEhIqPVDoMQu119/vWf+4bDDDmPMmDGMGDHCEwbngw8+YM2aNYwYMaLJdVxwwQV88cUXzJo1i759+3q83Hbv3s369ev54IMPmDp1Ko8++mitckOGDGHp0qUMHz6cX//61+zZs4fnn3+ewsJC7r77bj8xAcvLbPLkyR4PsksuucQvz88//8zhhx/OwIEDGTZsGD169GDv3r0sWLCArVu3csUVV3icAprCMcccw5133sl1111Hv379+M1vfkOfPn0oLi5m48aNvP/++4wbN85vPVB6ejpnnHEGjz/+OLNmzaJ9+/ZMmDChye1oDio6LYzMlCR2FZdGuxlRI1ESdXgtjrj55ps544wzmDVrFosXL2bOnDmUlZXRvn17hg4dyrRp0wJaFY3hkUceYfz48Tz66KO88847FBYWkpubS8+ePfnb3/4W8P7t2rVj0aJFXHPNNcyZM4e9e/cyaNAgrr766jrnOS644AJmzpxJ165dOfnkk/2u9+7dmxkzZpCfn8/ixYvZuXMnubm59O/fnzvvvJMzzzyzWX0FmDZtGmPHjmXmzJl89NFHvPLKK+Tk5NCtWzcuvvjioO2fMmUKjz/+OJWVlZx11lnNsriag4pOCyMrjresbggJkqDDa3HGwIEDeeihh+rN17t374CT/g51xfs68cQTOfHEhu1W4q5j/vz5DSrj4LhSX3jhhQEXVLZt25abb76Zm2++udFt8WXu3LmeoTBfxo0b53HWaCjjxo2rs75IoaLTwhjZO9ezSLQ1kiiJuomb0iKpqqrivvvuIykpKeDQmtIwWu+vWwvlzJE9o92EqKKWjtLS+Oijj3j//fd59913Wb58OZdffnmTIi0oFio6SotC53SUlsY777zDjBkzaNeuHRdddBF33313tJsU06jLtNKiUEtHaWlMnz4dYwwbN25k9uzZfjHelMahoqO0KNTSUZT4RkVHaVGopaMo8Y2KjtKiSExIjIvFoS3BNVVRQk0ovtcqOkqLIh4sncTERCorK6PdDEUJOZWVlQG3nWgMKjpKiyIeAn5mZ2d7tgxWlHhi7969zQrjAyo6SgsjHgJ+5ubmUlBQwM6dO6moqNChNiWmMcZQUVHBzp07KSgoIDc3t1n303U6SotiSMchZCZnRrsZzSI1NdUTdHLDhg21NuGKF8rKykhLS4t2MyJKa+5zYmIi2dnZ9OzZk9TU1GbdU0VHaVGcNSDgBrExR2pqKl27dqVr167RbkpYyM/PD7jnTDyjfQ4NOrymKIqiRAwVHUVRFCViqOgoiqIoEUNFR1EURYkYKjqKoihKxFDRURRFUSKGio6iKIoSMVR0FEVRlIghGqKjbkRkB7CxicU7ADtD2JxYQPvcOtA+tw6a0+dexpiOvokqOmFERJYZY0ZEux2RRPvcOtA+tw7C0WcdXlMURVEihoqOoiiKEjFUdMLL7Gg3IApon1sH2ufWQcj7rHM6iqIoSsRQS0dRFEWJGCo6iqIoSsRQ0QkTInKCiPwgImtF5NpotycUiEgPEVksIitFZIWIXGmn54rI2yKyxn5t5ypznf0MfhCR46PX+uYhIoki8pWILLDP47rPItJWRP4rIqvsz3t0K+jzn+3v9Xci8qyIpMVbn0XkCRHZLiLfudIa3UcRGS4iy+1rM0VEGtwIY4weIT6AROBHYH8gBfgGGBTtdoWgX12BYfb7bGA1MAi4G7jWTr8WuMt+P8jueyrQx34midHuRxP7/hfgGWCBfR7XfQbmAb+336cAbeO5z0A3YD2Qbp+/AEyJtz4DRwDDgO9caY3uI/AZMBoQYBEwvqFtUEsnPIwE1hpj1hljKoDngIlRblOzMcZsMcZ8ab8vAlZi/bFOxPqRwn49xX4/EXjOGFNujFkPrMV6NjGFiHQHJgD/cSXHbZ9FpA3Wj9PjAMaYCmNMIXHcZ5skIF1EkoAM4BfirM/GmA+A3T7JjeqjiHQF2hhjlhhLgZ50lakXFZ3w0A34yXW+2U6LG0SkN3AosBTobIzZApYwAZ3sbPHyHB4ArgFqXGnx3Of9gR3AHHtI8T8ikkkc99kY8zNwL7AJ2ALsMca8RRz32UVj+9jNfu+b3iBUdMJDoPHNuPFNF5Es4CXgKmPM3rqyBkiLqecgIicC240xXzS0SIC0mOoz1n/8w4B/GWMOBUqwhl2CEfN9tucxJmINI+0HZIrIuXUVCZAWU31uAMH62Ky+q+iEh81AD9d5dyxTPeYRkWQswXnaGPM/O3mbbXJjv2630+PhOYwFThaRDVjDpEeLyHziu8+bgc3GmKX2+X+xRCie+3wssN4Ys8MYUwn8DxhDfPfZobF93Gy/901vECo64eFzoJ+I9BGRFOBM4NUot6nZ2B4qjwMrjTH3uS69Cky2308GXnGlnykiqSLSB+iHNQEZMxhjrjPGdDfG9Mb6HN8zxpxLfPd5K/CTiPS3k44BvieO+4w1rHa4iGTY3/NjsOYs47nPDo3qoz0EVyQih9vP6nxXmfqJtjdFvB7Ab7C8u34Eboh2e0LUp3FYZvS3wNf28RugPfAusMZ+zXWVucF+Bj/QCA+XlngAeXi91+K6z8BQYJn9Wb8MtGsFfZ4BrAK+A57C8tqKqz4Dz2LNWVViWSwXNqWPwAj7Of0IPIwd3aYhh4bBURRFUSKGDq8piqIoEUNFR1EURYkYKjqKoihKxFDRURRFUSKGio6iKIoSMVR0lJhDRKaLiBGRudFuSzQQkcNE5DUR2SkiNfazmN6Acr3tvK3WZVWfQfRR0YlDRGSu84clIsvqCjsuIvNb8w94rCEi/YB84ESstTM7gW1AcTPvm2eL+SnNbWO0EJFT7D7kRbstSnCSot0AJewMB36LFdZDiX0uxoqA/CFwsrGiPzeUSqxFfoHIA27BijL8ctObF1VOwbuyPj9InrqegRIBVHRaB7eKyMvGmJr6syotnIPs1xcaKTgYK5LygJC3KIbQZxB9dHgtvnkfKMX6oTo7ym1RQkO6/dqs4TRFiRYqOvHNVqy4SADT7c2pGoxrXqh3kOtBJ2VFJN++NkVE2ojI3SLyo4jsE5F1InKriKS58h8jIm/ak+MlIvKBiPyqAW1MEGub4W/scrtE5FURqXNDLbvcefb2vDtEpEJEfhGR50VkVJAyHgcGu/zlIvKZiBTa6UPra69P/ReKyPsisltEykRkvYjMFpEDAuTfYD/nPDtpjuvz2dDAOv0+LycNa2gNYLLrvkE/fxE5SUReEZGt9rPbbjs3BNy22f4eGBHJt8/Psfu+y04/xU5PFJGjRORBEflCRLa5Ppv/E5GjA9w7z+6DM7R2i28f6noGAe53qD3X+ZOIlNvfyTdF5LQ6ymyw75sn1vbP99mfZ7mI/Cwij4kdyTlA2QT7+Sy2n0el/Z1cIdb20icEqzcmiXYAOj3CEtRvLlZgzueAXGCPfX5RgLzz7WtzA1wz9tE7SD29nTwBruXb1/6MFa3XYP13XuG676t23suwNkirdrXVAOXA2AD3nm5fn4e1zYLBGqsvdJWtAiYFaXc28LYrb41PvdXA5fXU+7KrngL7/dAGfj4ZwJuu+ip82r4PmOhT5nOsfyKc57fHPt8KfN7Aev0+L6zQ9Vvtz8ape6vP0cOVP9n1nTGutrjP7w5Q9xT7Wj4w0/Wcd9uvp9j5Bvvcq8zVNue43ufeY+x27nN9z2r1oSHfWfv6xXZ7nLoK7M/YOX+KANtSAxvs6+e63pfY7XfKrgfaBSj7tE//CrG++875p9H+TQnlEfUG6BGGD9UlOvb5dPt8E5DqkzfcolOIFbl3nJ2eAvweSyQMcBPWD+ntQFs7Ty/gE/v6ZwHuPd117yosYXP2tu8LvGVfLwX6Bij/f/b1b7CiZDtl2wLX2X/w1fgInqveIvvH5A9Ahn2tE9YWvg35fB7F+4N6ifOZAAcCi10/WAfW8VynNOF7Udfn5fTN73vgk+9+1w/oWUCWnZ6F9YPtCNBZPuWmuJ5dDXCz6/NuA3RyPYMXsLzzOoMnKHEn4Eb7864BRtXxvZ/exGcwBq/gvAh0d/XterteA9wYoOwGvCL1FTDaTk8CTsb7j8ndPuWOwCvAVwHZdroAXbGst3vD/ZsRySPqDdAjDB+qv+i0AXbZaVf65A236FQCBwS4/rjr/k8EuN7L9Ufe0+fadFdZv20jgDQsoTPAf3yuHYv3RzM3SL+usfMsqKPei5v42fRy/bBdEuB6BtZe9AZ4so7nOqUJddf1eTl98/seuPL0s9teAOwfJM/v7Pt855M+xfXsbm/Gd/sm+x5z6vjeT2/iM3jXvvYRga2Z2/EKZxufaxvsa1uB9gHK/tW+vi7Id21RU59JrB06p9MKMNaW0nfbp9eLtd99pHjRGLM2QPo7rvd3+F40xmzE+vEFa8glEKXAAwHKlgH/tE9PE6m1Tmmy/TrXGLM7yH2fsV+PEpHEANd3AU8EKVsfp2LNpW4F/uN70RhTivezOjVI/dHifKy2v2yMWRckz/+wLMWDgsxhVAP3BUhvKK/Zr2ObcQ8/RCQXOMo+vcMYUx0g211Y1mkWloUciNnGmF0B0l+2X/v4/P052713EpFW8XvcKjqpAPAQ1iLCTsAVEax3eZB0Z0vcMrzi4ss2+7VdkOvLjDElQa69b7+2xdr33mGM/fpnexLc78DavAwsq6N9kHqrgtRbH8Ps1w+D/LABvGe/ZgL9g+SJBs6zO72OZ7cZa94Ham917LDWGLOzrkpEJF0s55B820Gh0jX5/5Wdbb+Q9MjLoVhDWgbvd6cWxpg9wBf26bBAebDm3gLxs+t9W9f7d7CGl4cB+SJyroiEum8tCl2n00owxpSKyO3Ag8DfRGSW/UcUbrYESXd+cLcZe5yhjjzJQa7/HCTd91pHwPnP3PnvO8c+6iMjQNqOBpQLRkf7ta62bw6QvyXgPLss+6iPRj872zrKx5rbcSjBGtKrARKBDliCHEqc57zHGFOXO7rz2QT7XIoCJRpjylwGd7Irfa2I/AHLy/RX9oHtkfgGluX0FXGEWjqti38DP2FZDn+NclvCTbDQP853fqIxRhpwbAhwj2AWSmNIreNaMBGONs6zu7KBzy4/wD3qe3YPYAnOOuA0rHm3LGNMJ2NMF+DwUHUmCHV9LmHBGPMEljV+FfAK1vBtb+BS4AsRuT7SbQonKjqtCGNMOfB3+/QqEelQTxHnByItyPWGWArhpK5hCPd8gvu/a2fIblDom9MgnLb0qiOPe1iqOVZVqAnrsxORFGCifXqOMeZ/xpgCn2ydw1E33uecLiJ1WZfdffKHBGPMNmPMg8aYU7CsqJFYXpYC/F1EhoSyvmiiotP6mAP8iLVW5dp68hbar92DXD8sRG1qKoeJSKAhHIAj7ddCLE81hyX2a9CFfmHmS/t1VB1tdxZAlhC5OGFOiKSgwWHxPruTRCTYkGdz6IDX0gg2pHRsHeUb0odgfIXXwjwqUAYRycGKZQjezzHkGIvPgTOwhvMSgHHhqi/SqOi0MuwJ8On26WXUbS04TgATfS+ISCrWcEA0yQCu9E202/YX+/S/PnNGc+3XESJyfl03F5FgDgzN4X9YP47tsda1+NaZAfzNyVuHs0Gocbyo2taRZx5W2/fDWs8UlCY+u714f/gPDnDPrsCf6ikPdfchILYn42L7dFoQT7JpWFZ/MbCwsXUEwrbugrWpGmvJAURh2C9cqOi0Tp4BvseK4xXwvzqbF+zXi0Rkqv1jjogchPVHF20vmz1YQw9Xikg6gIjsjzUuPhDLM+5OdwFjzBt4I24/ISIz3K69ItJORCaKyCs0z7U3ILYr+Gz79E4Rudj1XA8EXgcOwHIHvy3U9dfBCvt1nFjbJ/hhjFmJ10V9hog8Yj9vAEQkS0SOE5GnsBZXNgp7Av9T+/QJscMK2WFijsHyKqvLinH6cEKwkDP1cBOWqA4DnhOR7nb9Wfa8ijMycKe9DCEU3C4i/xVrW4ZcJ1FEOovITKy5HoMVQSM+iPZCIT1Cf+CzODRIntOoHXpjboA8yVg/Ak6eSrwrzndhWUD1LQ6dEqT+PPv6hjraGPAe1A5H8z/7fQXeVd8Ga+X6mUHum4k3KoFzFOIfzmVOkHr9nlUjP58MvFETArW9DJ8wOA19rvXU27uOzysZ76LUGiyX9g320d2VLxGY5fOc9uL1LnPSFvvcf4qdnl9PG0dhCa5zn2LXeX3fuQ54F0FXY3lObnB/x+p6Bvb1S/Au3q3BCtPjDoMzn7rD4OTV0TfnHr1daQ/4PMs9eC0+57i+rmcWa4daOq2X/1HPuLQxphI4DrgH64+qBmueYS7W2PY3YW1h/Risce+/YMV3S8H68VsAjDHGPBewkDElxpjfYoVa+R+W+3K6XX4tliV4OtbwY+gbbS0AHY8VDuhDrB/VDGAj1oLRg40xr4Sj7jraVAkcgxVb7GcsD8de9pHkyldtjLkMa45hvt3mFKzntwlLzCdj7W3TlHYsBUZjLaYswBLD7Viel0Op4ztnrPU/R2F9pjuwJuSdPjS0/n9jzVU+gyVaWVhC8DZwhjHmXBPaIc/7sdbNvQKsxrLkUrG8TJ8HjjDG3B7C+qKOE9dIURRFUcKOWjqKoihKxFDRURRFUSKGio6iKIoSMVR0FEVRlIihoqMoiqJEDBUdRVEUJWKo6CiKoigRQ0VHURRFiRgqOoqiKErEUNFRFEVRIsb/A3wS6VGonozFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy([np.linalg.norm(A @ x - rhs) for x in conv_cg], label=\"CG\")\n",
    "plt.semilogy([np.linalg.norm(A @ x - rhs) for x in conv_rich], label=\"Richardson\")\n",
    "plt.semilogy(res_cheb, label=\"Chebyshev\")\n",
    "plt.legend(fontsize=20)\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Number of iterations\", fontsize=24)\n",
    "plt.ylabel(\"$||Ax_k - b||_2$\", fontsize=24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "8.59 ms ± 402 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "10.6 ms ± 110 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit np.linalg.solve(A, rhs)\n",
    "%timeit spsplin.cg(A, rhs)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}