mscroggs/lecture5.ipynb

## lecture5.ipynb
{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyMHb2uYBz0glmyJedvMKy7T",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/mscroggs/0eeca24839cc19da79bb8acbf861d956/lecture5.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Scoping variables in Python"
      ],
      "metadata": {
        "id": "KuNjMCrw6ZFv"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "\n",
        "a = [5]\n",
        "matrix1 = np.random.rand(2, 2)\n",
        "\n",
        "def do_timing():\n",
        "  print(matrix1)\n",
        "\n",
        "def f(n):\n",
        "  global matrix1\n",
        "  matrix1 = np.random.rand(n, n)\n",
        "  do_timing()\n",
        "\n",
        "f(3)\n",
        "print(a)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "LKZ40TO7FFII",
        "outputId": "c400318c-9252-415b-9f15-964ea97889e9"
      },
      "execution_count": null,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[[0.12018935 0.97651145 0.63981664]\n",
            " [0.99614862 0.66402939 0.81235863]\n",
            " [0.05879349 0.85772499 0.84502935]]\n",
            "[5]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Looking for a solution in a Krylov subspace"
      ],
      "metadata": {
        "id": "8RzTNsf-6eY4"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import scipy\n",
        "from scipy.sparse import linalg\n",
        "import matplotlib.pylab as plt"
      ],
      "metadata": {
        "id": "AW-b3xq2GJfb"
      },
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "N = 100\n",
        "\n",
        "w = np.linspace(0.1, 1, N)\n",
        "A = scipy.sparse.diags([w], [0])"
      ],
      "metadata": {
        "id": "IpCa-SIDLHAY"
      },
      "execution_count": 2,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "b = np.random.rand(N)"
      ],
      "metadata": {
        "id": "WeDT6qM1Lc-q"
      },
      "execution_count": 3,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "Niter = 100\n",
        "\n",
        "residuals = []\n",
        "\n",
        "for m in range(1, Niter):\n",
        "    V_m = np.zeros((N, m))\n",
        "    V_m[:, 0] = b\n",
        "    for i in range(1, m):\n",
        "        V_m[:, i] = A**i @ b\n",
        "    # print(np.linalg.norm(V_m[:, -1]))\n",
        "    b_m = V_m.transpose() @ b\n",
        "    A_m = V_m.transpose() @ A @ V_m\n",
        "    y_m = np.linalg.solve(A_m, b_m)\n",
        "    x_m = V_m @ y_m\n",
        "    residuals.append(np.linalg.norm(b - A @ x_m))\n",
        "    print(residuals[-1])"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "CU57EPTeLseK",
        "outputId": "b69673a3-31d2-41ab-c385-2260bc48d16b"
      },
      "execution_count": 4,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2.5436061592773647\n",
            "1.5268525024832769\n",
            "0.7909911169131705\n",
            "0.4696176583131388\n",
            "0.24136040736392161\n",
            "0.12512897553534058\n",
            "0.06089175014831749\n",
            "0.032984590120863536\n",
            "0.017464647397230564\n",
            "0.009580523399923897\n",
            "0.005047419612004152\n",
            "0.008617175470525375\n",
            "0.014783316177311842\n",
            "0.025198894401631873\n",
            "0.0036891340482487853\n",
            "0.005562739997005804\n",
            "0.002479379682888194\n",
            "0.003919095355641206\n",
            "0.0034913902620380024\n",
            "0.004718870348811943\n",
            "0.003927381162396119\n",
            "0.004347369409142133\n",
            "0.0071584231966449795\n",
            "0.003389033975667293\n",
            "0.00403510991747916\n",
            "0.0032925436306968244\n",
            "0.0019120985755934133\n",
            "0.008922360863554501\n",
            "0.002121822396769439\n",
            "0.006224355068486315\n",
            "0.0019758500698531397\n",
            "0.001963828857433013\n",
            "0.0070520293245848725\n",
            "0.00203227604965537\n",
            "0.0026868682602188216\n",
            "0.0035780496878030856\n",
            "0.0028959913703471164\n",
            "0.0038140883124964118\n",
            "0.006495752839612724\n",
            "0.0048662259052545155\n",
            "0.004726850533732949\n",
            "0.0038550919564787974\n",
            "0.0020594552139296\n",
            "0.009681176697094142\n",
            "0.004059547004115637\n",
            "0.0033004224357727615\n",
            "0.07316329312750489\n",
            "0.0033458304053946893\n",
            "0.00811556851128525\n",
            "0.006131606302736748\n",
            "0.0029345632384562563\n",
            "0.002756789832940716\n",
            "0.004777825867425813\n",
            "0.003400581182203216\n",
            "0.005187433991987029\n",
            "0.0031608240602008317\n",
            "0.0040055504811540195\n",
            "0.003237239874247785\n",
            "0.00550471820380855\n",
            "0.005414868084347376\n",
            "0.004455365760456202\n",
            "0.0040602901207817295\n",
            "0.002597693094856829\n",
            "0.003065305480807656\n",
            "0.021700095066164726\n",
            "0.0026889043746681144\n",
            "0.005619051527126164\n",
            "0.005314573214248376\n",
            "0.04298796417498151\n",
            "0.0026243556719718875\n",
            "0.00836044845982334\n",
            "0.056348249898451774\n",
            "0.003775520077539402\n",
            "0.0020734416316075858\n",
            "0.014356350796210008\n",
            "0.002367605961017925\n",
            "0.0057151517855772\n",
            "0.007445044484260813\n",
            "0.0028913160194635433\n",
            "0.07644090077471447\n",
            "0.005017306724747539\n",
            "0.014888283385417408\n",
            "0.009997656013732007\n",
            "0.0020723388204639996\n",
            "0.0037543932257588634\n",
            "4.313020978658278\n",
            "0.024943760455354562\n",
            "0.00801747622377593\n",
            "0.020698039406411163\n",
            "0.05559842930279566\n",
            "0.016947654180791555\n",
            "0.0021510712095071044\n",
            "0.0038086454077604407\n",
            "0.004799665852173488\n",
            "0.002574022384034436\n",
            "0.005931748496574927\n",
            "0.0026713219004615204\n",
            "0.006409179520646468\n",
            "0.008519140241483453\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.plot(range(1, len(residuals) + 1), residuals)\n",
        "plt.xscale(\"log\")\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        },
        "id": "NaWFemS1M6PI",
        "outputId": "28ce1163-9127-4199-aabb-00d7dd0ac757"
      },
      "execution_count": 5,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3yb9bU/8M/R9pA8ZccrcRzHziQDBwiQAYQww4YAvS20FEqBW3ppe0t77/396Kbtr6WU0TaMAm3Zl1BWgDKSEAgkDglkOnGWR5x476H1/f3x6Hn0SJYcy9Z4ZJ3368ULIsnSN4k4Ojrf8z0PCSHAGGNs4tPFewGMMcZigwM+Y4wlCQ74jDGWJDjgM8ZYkuCAzxhjSYIDPmOMJQlDvBcwktzcXFFaWhrvZTDGWMLYtm1bqxDCHuw+TQZ8IloFYFV5eTmqq6vjvRzGGEsYRHQ01H2aLOkIIV4XQtyWkZER76UwxtiEocmAzxhjLPI44DPGWJLggM8YY0lCkwGfiFYR0Zqurq54L4UxxiYMTQZ83rRljLHI02TAZ4yxiehE9yDaeofi9vqa7MNnjLGJ6PRfvg8AOHL/JXF5fc7wGWMsSWgy4POmLWOMRZ4mAz5v2jLGWORpMuAzxhiLPA74jDGWJDjgM8ZYkuCAzxhjSYIDPmOMJQkO+IwxliQ0GfC5D58xxiJPkwGf+/AZYyzyNBnwGWOMRd6EHJ52xz+2Qa/TYXmFHUsr7LBbzfFeEmOMxd2EC/hCCNgsRry39wRe/+IYAGBOkQ3LK/KwvNKO+SWZMOj5iw1jLPlMuIBPRLj/6lPg8QjsaerG+ppmrK9pwaPra/Hwh7WwWQxYUmHH8go7llXakWe1xHvJjDEWExMu4Mt0OsKcogzMKcrAXedOR1e/Ex/VtmB9TQs27G/Bm182AQBmFdiwvNKO5ZV5WDiZs3/G2MQ1YQN+oIxUIy49pRCXnlIIIeTsvwUbalrwl42H8Oj6g7BaDFgyPRfLK/KwrNKOfBtn/4yxiSNpAr4aEWF2YQZmF2bgznPK0TXgxMe1rVhf04wN+1vw1s7jAIAZk6xYXinV/k+dkgUjZ/+MsQSWlAE/UEaKERfPLcDFcwsghMC+4z1YX9OC9TXNePyjQ/jzhoOwmg04qzwXyyul2n9BRkq8l80YY2HRZMAnolUAVpWXl8fjtTGzwIaZBTZ8e/k0dA868Ultq/cDoAVv75ay/8p8qxL8q6Zkw2Tg7J8xpm0khIj3GkKqqqoS1dXV8V6GQgiB/Sd6lc6f6qPtcLoF0s0GnDktRyn/FGZy9s8YG6703jcBRPci5kS0TQhRFew+TWb4WkVEqJxkReUkK761bBp6h1ze2n8LNtQ04909JwAAFfnpWF6Zh2UVdlSVZsFs0Md55YwxxgF/XNLNBlwwexIumD0JQgjUNvdKpZ/9zfjrx4exZuMhpJr0OHNarrf1047irNR4L5sxlqQ44EcIEWF6vhXT8624dWkZ+oZc+ORgm1L+eW+vlP2X56Urh75Om5rN2T9jLGY44EdJmtmA82fl4/xZ+RBC4GBLn9L2+czmo3h802GkGPXe2r908Kskm7N/xlj0cMCPASJCeV46yvPS8c0lZeh3uLD5YJtS/nl/XzOA3Sizp+FrZ0zBVxeXQq+jeC+bMTbBcMCPg1STAefNzMd5M6Xs/1BrHzbUtGDdribc9/oevLK9Eb+8ci7mFPH1ABhjkcPN43FGRJhmT8c3zp6KF7+1GH+8YQGOdQ7isoc34edv7EHfkCveS2SMTRAc8DWEiHDZvEK8f88yrF40GY9vOoyVD2zE+94NX8YYGw8O+BqUkWrEr66ai5dvX4xUkx63PF2NO/6xDSe6B+O9NMZYAuOAr2FVpdl48ztL8IMLKvHe3mas+N0G/G3zEbg92j0dzRjTLg74Gmcy6HDnOeV497tLMa8kE//zz924+k+fYG9Td7yXxhhLMBzwE0Rpbhr+dstp+MPq+ahv78elD23Cr9btxYDDHe+lMcYSRMwCPhGlEdHTRPQYEX0lVq87kRARrlhQhPfuWYarFxbhLxsO4fwHNmB9TXO8l8YYSwDjCvhE9CQRNRPRroDbLySiGiKqJaJ7vTdfBeBlIcStAC4bz+smu6w0E35zzTw8f9sZMBt0uPmvW3HXs5+juYc3dRljoY03w38KwIXqG4hID+ARABcBmAXgBiKaBaAYQL33YVyHiIAzynLw1t1L8B8rKvDu7hNY8bsNePazOnh4U5cxFsS4Ar4QYiOA9oCbTwNQK4Q4JIRwAHgewOUAGiAF/XG/LvMxG/S4e8V0rPvuEswqtOHHa3fi2r9sRs3xnngvjTGmMdEIvEXwZfKAFOiLALwC4Goi+hOA10P9MBHdRkTVRFTd0tISheVNTNPs6Xju1jPw/66dh0Mtvbjkjx/ht+/sw6CTv0wxxiQxm6UjhOgD8PVRPG4NgDWAdMWraK9rIiEiXHNqMc6dkYdfvLkXj3x4EG982YSfXzEHS6bb4708xlicRSPDbwRQovp1sfc2FiPZaSb87rp5ePabp0NHhK8+sQXffX47WnuH4r00xlgcRSPgbwUwnYimEpEJwPUAXgvnCYhoFRGt6erqisLykseZ5blYd/cSfOfccry5swnn/W4DXthaBy1fx5gxFj3jbct8DsBmAJVE1EBEtwghXADuAvAOgL0AXhRC7A7neYUQrwshbsvI4PHA42Ux6nHPykqsu3sJKvOt+OH/7sTqNZ+itpk3dRlLNuOq4Qshbghx+1sA3hrPc7PIKs+z4vnbzsBL2+rxy7f24aIHP8J3V1Tg9mXT+GIrjCUJTbZHckknOnQ6wupFk/H+95Zh5axJ+O07Nbh+zWbUt/fHe2mMsRjQZMDnkk505aab8fCNC/DA6nnY19SDix78CC9va+DaPmMTnCYDPos+IsKVC4qVA1vff+kL3PGPz9HR54j30hibkLSQUGky4HNJJ3aKs1Lx3K1n4N6LZuC9vSdwwR82YsN+PvDGWKRpIN5rM+BzSSe29DrC7cum4dU7z0JGihE3PbkF9722m0/pMhZBGoj32gz4LD5mF2bg9X8/G18/qxRPfXIElz60Cbsa+VsWY5HAJR2mORajHv931Wz87ZbT0DPoxJWPfoxH19fyZRUZGyct/B+kyYDPNfz4WzLdjne+uxQrZ03Cb96uwQ1rPuX2TcbGQQMJvjYDPtfwtSEz1YSHb1yA3183D3ubunHRgx/hf7l9k7ExEaocP17/D2ky4DPtICJctdDXvvm9l77Anc9y+yZj4VLH+HjlTBzw2ajI7Zs/vHAG/rVHat/cyO2bjI1JvL4jc8Bno6bXEb69fBrW3iG1b36N2zcZGzV1Vu/hko4Pb9pq25wiqX3z5jOl9s1VD21CbXNvvJfFmKb51/DjswZNBnzetNU+i1GP+y6T2jc7+h246cktaO4ejPeyGNMsD2f4LNEtmW7HU1+Xgv43nt6KviFXvJfEmCZpobuNAz4btzlFGXj4xgXYc6wb33luOx/SYiwI9f8VnOGzhHbujHz85LLZeH9fM37y+m5NZDOMaYkW2jLHdcWraCGiVQBWlZeXx3spLAxfXVyK+o4BrNl4CJOzU/HNJWXxXhJj2iGC/mdMaTLD503bxHXvhTNw0ZxJ+MVbe/H2rqZ4L4cxzVB36XBJh00IOh3hgdXzMb8kE3c/vwPb6zrivSTGNEELJR0O+CziLEY9Hv9aFfJtFnzz6WrUtfHQNcbUMZ5n6bAJJSfdjKe+vghuIXDzU1vQ2c+zd1hyUwd5zvDZhFNmT8ear1ahoX0At/1tG4ZcPIKBJS9uy2QT3mlTs/Hba0/BlsPt+M+Xv+R2TZa01EE+Xv8XaLItk00sl88vQkPHAH77Tg0mZ6fieysr470kxmJPA6MVNBnwuQ9/4rlj+TQ0dPTjoQ9qUZKViusWlcR7SYzFlAj5i9jRZEmH+/AnHiLCTy+fg6UVdvx47U58dIBn6bPk4j8eOT5r0GTAZxOTUa/DIzcuQHleOu74++fYd7w73ktiLGb8xiPHKcXngM9iymox4q9fX4RUsx7f+OtWnOCRyixJcIbPklJBRgqevHkRugac+MZTPFKZJQc+eMWS1uzCDDz8lYXY09SNn76+J97LYSzq+OAVS2rnVObhzuXleKG6Hm9+yYPW2MTGs3RY0rt7xXTML8nEj175Eo2dA/FeDmMxwSdtWVIy6nV48Pr5cHsE/uP5HXy1LDZhaeGkLQd8FndTctLwsyvmYMuRdjz6YW28l8NYVAgNnLTVZMAnolVEtKarqyveS2ExcuWCIlw+vxB/eP8APucZ+mwC8u/Sic8aNBnw+aRt8iEi/OyKOSjIsODu57ejZ9AZ7yUxFlH+XTqc4bMkZ7MY8eD183GscxD/55+7470cxiJKhPjvWOKAzzTl1CnZ+M6507F2eyPWbm+I93IYixiu4TMWxJ3nTMOi0iz8z6u7+fKIbALhg1eMDWPQ6/DA6vkgAu5+YTucbk+8l8TYuHGGz1gIxVmp+OWVc7G9rhMPvX8g3sthbNy4S4exEayaV4hrTi3Gwx/W4rNDbfFeDmPjwqMVGDuJ+y6bjcnZqfiPF3agq59bNVni8j9pyyUdxoZJNxvw4PUL0NwzhB+t5Yugs8TF8/AZG4V5JZm4Z2UF3tp5HC9Vc6smS0x+V7ziTVvGQvvW0mlYXJaD+17fjSOtffFeDmNh4wyfsVHS6wi/Xz0PBOBX6/bGezmMjRNn+IyNqCAjBbcvm4Z3dp/AlsPt8V4OY2FJqgyfiMqI6AkiejlWr8kmnm8uKUO+zYxfvLkHHp6dzxKISJSTtkT0JBE1E9GugNsvJKIaIqolontHeg4hxCEhxC3jWSxjKSY9vr+yEl80dOGNnXxZRJY4Eumk7VMALlTfQER6AI8AuAjALAA3ENEsIppLRG8E/JMX0VWzpHbVwmLMLLDhN2/vw5DLHe/lMDYqCXPSVgixEUBg0fQ0ALXezN0B4HkAlwshdgohLg34pznC62ZJTK8j/NfFM9HQMYBnPjka7+UwNioiwQ9eFQGoV/26wXtbUESUQ0R/BrCAiH40wuNuI6JqIqpuaWkZx/LYRHb29Fwsq7DjoQ8OoKPPEe/lMHZSCZPhR4IQok0IcbsQYpoQ4lcjPG6NEKJKCFFlt9tjtTyWgH588Uz0Drnw0Ad8HVymff5XvIrPGsYT8BsBlKh+Xey9jbGYqJxkxXVVJfjbp0dwtI0PYzFtS6RN22C2AphORFOJyATgegCvRWJRfBFzNlr3nF8Bg06H37xdE++lMDaihLnEIRE9B2AzgEoiaiCiW4QQLgB3AXgHwF4ALwohInIhUr6IORutPJsFty0tw5s7m7DtaEe8l8NYSAmT4QshbhBCFAghjEKIYiHEE97b3xJCVHjr8r+I7lIZC+62pWWwW8345Vt7eZom0yy/92YC1vCjhks6LBxpZgO+d34Fth3twNu7jsd7OYwBAB7beAjX/WWz8mt1jNd0hh9rXNJh4bq2qgSV+Vbc//Y+OFx8Ddxkdd9ru7Fm48F4LwMAcKi1F7XNvcqv+YpXjEWIXkf40cUzcLStH3//lA9jJatNta3YclgbezkeD+D2BD9sxRm+Cpd02Fgsq7Dj7PJc/PGDA+ga4MshJiOX2wO3Rxvf8NxC+A/4i38JX5sBn0s6bCyIpCy/a8CJRz/kw1jJyOURcGlkiqrHI+D2G6fgw1e8YiwCZhdm4OqFxfjrx0dQ394f7+WwGHO5BVxujQR84f/h40nwk7aMadL3VlZApwN+8w4fxko2UoavlZIO/Eo6SXUBlHBwDZ+NR0FGCm5bUobXvziGz+u0sYHHYsPl8cAZwwx/wOFGc89g0Ps8YoSSTgJOy4waruGz8frWsmnIs5rxszf28GGsJOJyC7/OmGj704aDuOZPm4Pe5/EICOGr16vfhx4B3PPCDvz09T0xWadMkwGfsfFKMxvw/Qsqsb2uE298yVfGShZShh+7kk5r71DI8dxyzV7+AArctH1leyOe/PhwtJfohwM+m7CuXliMWQU23L9uHwadfGWsZBDrDN/l9oTsCpI/d5SyDh+8Yix69DrCf18yE42dAzHPpFjsCW9XTCzbMl1u/zp94HoA6QAWEHARc67h+/CmLYuUM8tzsWJmPh798CBaeobivRwWRXJmH8suHZcn4HCVej1ySUep4fvui1cjkSYDPm/askj68cUzMOh044H39sd7KSyK5Mw+ln34Lk/oko58s1LD55O2jEVfmT0dX108Bc9vqUPN8Z54L4dFiRLwo1zScbo92NvU7f1vuWwz/DXl2zxBNm15lg5jUXT3edNhtRjx8ze5TXOicnl3SV1R7tJZt+s4Ln1oE9r7HL7XDBLw5cxeLul4NJDic8BnSSEz1YS7z5uOjw60Yv3+lngvh0VBrDL8rgEn3B6BfodLea1gGbtHBGT4qoc4Ytg6qqbJgM+btiwa/u2MKZiam4ZfvLk36lkgiz25dh/tGr7b+95xe4TS8x+sFdQTsGmrTuvj1SasyYDPm7YsGkwGHX500QzUNvfiua318V4Oi7CRgm8kyVm92+Mb1BbsW8VIm7ZDcbpIjyYDPmPRcv6sfCwuy8ED/9rPM/MnGDmwOqPc8yhv1LpVPf/BNm3dyn3Sr9WP4AyfsRggIvzXJTPR0e/gmfkTjNx/L0TwABwp8gVW3MI3mTPY4SsxQh8+Z/iMxcicogxc452ZX9fGM/MnCnVZJZpZvlO1V+BSZfuB3MNm6QSv4cdyP4kDPktK37+gEkTAwx8eiPdSWISoN2ujWceXs3qPOMmmrTeOe4Jk+OqAH8txzhzwWVLKt1lww2mT8crnjWjsHIj3clgEqKdkRjOI+m3aekJn+CNNyxx0+tYayxZNDvgsad22tAwA8NjGQ3FeCYsEddCNaoav3rQdqaQT8GGgPvDnn+EnecDnPnwWC4WZKbhqYRGe21LHg9UmAHVWH826uEtVxhlp09YTcMLWr6TjUn8bSfKAz334LFa+vbwcTrcHT2zi8cmJTp1lR/O0bbA+/OAlHf91hdq0dbq4hs9YTEzNTcMlpxTi758eRVc/9+UnMnVnTjRP2ypB/mSbtiNk+EOqgM81fMZi6I7l09A75MLTm4/EdR0ej8BD7x/AoZbeuK4jUamDfDRn4o9209ZXw5d+7d+lo9q0jWFPPgd8lvRmFtiwYmYenvz4MPqGXHFbx7a6DvzuX/vx6o5jcVtDInOrM/wYtGWerKQjhpV0fPqdvvdZ0tfwGYu1O84pR2e/E89+Vhe3Nazd3ghAujA2C5960zaaQVQd5J0jbNoqoxWUko7vMX1D3KXDWNwsnJyFM6flYM1Hh+Iy52TI5cabXzYBAFq5Y2hMYtaW6Q3yTrcYlsWrjdSH3zvoy/C5hs9YHNx1Tjlaeobw8raGmL/2+poWdA04kWrSo4Uz/DGJ1sGrurZ+3PZMtZIIyBn+kMuXGIxqPHKIefh80paxOFg8LQcLJmfizxsOxvRrNgC8ur0ROWkmrJiZz2cCxsgVpQx/W1073t1zAnXt0twlp0cO+L73yMjTMoe3Zao5edOWsdgjIty5vBwNHQN4LYYbp10DTry/txmr5hViUoYFLT1DfBnGMVAH/EgevJIzcF8LpvRvdWtluPPw/Z8/yQM+n7Rl8XLezDzMmGTFo+trozpiV23dziY43B5cuaAI9nQzhlwe9MaxWyhRqYN8JLt0Ajtx5A8A9WnZoCdtAzdtvbfrdeT3uKSv4fNJWxYvRIQ7zynHwZY+vLP7eExec+32RpTlpuGU4gzYrWYA4LLOGPiftI1khu/bpAV8HyxDql56d5A6vG/T1v/XZoN/2OU+fMbi6OK5BZiam4aHP6yNemmlsXMAnx1uxxULikBEHPDHwX+WTuT+3uSA7wo4Veu3aRusLTPEBVDkgC//mzdtGYsjvY7w7WXTsPtYNzYeaI3qa/1zh9R7f8X8IgBAbro34HOnTtiiVtLxPpcrsKTjHHnTVr7JE9CWaTboAQBWi8H7fJzhMxZXly8oRJ7VjMc/it7oZCEE1n7eiKopWZickwoAnOGPg9+mbURr+HJJJ3SGH3TTNnDsgjfFN3kz+1QTB3zGNMFs0OPms0rx0YFW7D4WneaBPU3dONDciysWFCm3ZaYYYdARn7YdA5ff8LTIBVGH6pKGgG9Im19b5gjjkd0iMMOXA76U6UfzkFggDviMhfCV06cgzaTH4x9FZ3Tyq9sbYdQTLplboNym0xFy082c4Y9BtDN8+QPFd/BKtWkb8HpCCKWk0z3ghNPt8dXwjf4BP5pzfwJxwGcshIwUI1YvmozXvziGYxG+DKLL7cE/dxzD8so8ZKWZ/O6zWzngj4UrSpu2gTV8ObgPBvTht/UOKcP31DH852/uxV3Pfq40AJj0UthNM0slnVi1/wIc8Bkb0TfOLoUA8NQnRyL6vBsPtKC5ZwhXLywadl9uuok3bcdAXcZxR7AtU26bdAUcwAo8aXvTX7fgt+/USL8OKPG8s/vEsE3bFCNn+IxpSnFWKi6ZW4BnP6tD92DkLpDy4tYG5KSZcO6M/GH32a1mtPY4IvZaycLlETB4DzVFstXRNyxNLu14SzrDMnwH2vqkv7eRxiXLJR2LUQ8dBa//RwsHfMZO4tYlZegdcuH5LZEZndzWO4T39p7AlQuKlI4NNbvVjNbeoZh+1Z8IXG4BizHyG6FyZq+UdryBfzBg09bpFsp9wWJ44Kat2aCDXkec4TOmJXOLM7C4LAdPbjoSkVORa7c3wuURuG5RSdD77elmuDwCnQN8ycVwuDwCFm/27IzoSVv/QB8sw5cvaO5UXf4wkFBO2kofSmajFPC5hs+Yxty2tAzHuwfx5s7xDVUTQuDF6nrML8lERb416GNyuRd/TFwejxJMg406GKthoxW8AdoR0KXjdgtl7yBYmSbwpK1Jr4eeOMNnTHOWV9oxPS8dazYeHte4hS8burD/RC+uqwqe3QNShg9wwA+XVNKRM/zI1/DdASWdwLZMp8ejBO9gWbv8rUOu4csZ/oTswyeiK4joMSJ6gYhWxup1GYsEIsKtS8uwt6kbm2rHPm7hhep6WIw6XDqvIORj5NO2fPgqPC6PB0a9HESjMB7Z44HH4+uvDzxp6/YI5dtAsBje4d3QNej8a/iaC/hE9CQRNRPRroDbLySiGiKqJaJ7R3oOIcSrQohbAdwOYPXYl8xYfFw+vxB2qxlrNo5t3MKAw43XdxzDxXMLYLMYQz6OxyuMjcstoNcRDDqKbB++29eWqS6/BM7SkTZtA8YpqMgdPOSdjmzQEfQ6XdB6f7SMNsN/CsCF6huISA/gEQAXAZgF4AYimkVEc4nojYB/8lQ/+t/en2MsoZgNetx8pjRuYW9Td9g///buJvQMuUYs5wBAutkAi1HHvfhhcnkEDHodDDqKaFumetNWPb5BneErdX7lYidBAn6vf6stEUGvi+x+w8mMKuALITYCaA+4+TQAtUKIQ0IIB4DnAVwuhNgphLg04J9mkvwawDohxOeR/W0wFhv/dvoUpJr0eGwMQ9Ve3NqAKTmpOH1q9oiPI+LxCqOlvlCMy+OBUUcw6HURLun4grn6g0T930PK4SxvvT9IwG/3ZvjyXURSeUeLGX4wRQDqVb9u8N4Wyr8DWAHgGiK6PdSDiOg2IqomouqWlpZxLI+xyMtINWL1ohK8tuMYTnQPjvrn6tr6sflQG649tRhEdNLHy734LLTGzgHM/8m72Ha0A4AUgPU6glFPEd609WX4oertQwGncYM9Tj2KAQAIBJ1ugg5PE0L8UQhxqhDidiHEn0d43BohRJUQosput8dqeYyN2k2LS+HyCLy8rWHUP/PStnroCLj61OJRPd7OGf5JHe8agMsjcLi1D4AUOJVN22jU8D0i5BROubwjd+IES9oDxyArGX6CBPxGAOpiZLH3tnHja9oyLSvNTcPpU7PxUnX9qFo03d4Ph6UVdhRkpIzqNRJpgFrXgBPzf/ouNh9si+nrDjg8yusDUmCWNm11ET14pR6PHOybg458lzscKcMPPLRH3p9NlIC/FcB0IppKRCYA1wN4LRKL4mvaMq1bvagER9r6seVw4NbWcBv3t6CpaxDXnjryZq2a3WpGe78jphfHGKtjnQPo7Hei5nj4G9njIZdIuvql2rjLI2DUEwz6yLY6qscjB/vmYDbolZKO/LrBDl45An42VIb/9CdHcPavP0C/I/IXsh9tW+ZzADYDqCSiBiK6RQjhAnAXgHcA7AXwohBid8RXyJgGXTSnAFazAS9U15/0sc9sPgK71YzzZw0flBZKbroZQvg2+rSs25thdw9GPkCNZEAO+EqGL2DQ6ZS2zEc+rEVdW/+4X0d9acNg3xxMBp2vpOMOfdJWvk/ewpFq+MNP2jZ2DqC5Z0iZphlJo+3SuUEIUSCEMAohioUQT3hvf0sIUSGEmCaE+EWkFsUlHaZ1KSY9Vs0vxFs7m9AzwhTNo219WL+/BTeeNjnooLRQEqkXv8cb6LtjPPsnMOA7PR7o9VJJ50T3IH77Tg2e3zr+gXfqi5gHZuNmg/QBo2zaKhl+6Oc5c1ouAGDB5EwYghwSa+0Zgj3dPKrN/XBpcrQCl3RYIriuqgSDTg9e/6Ip5GP+/ulR6Ilw4+mTw3puJeAnQKeOPDY6kuOjR2MwIOC7PcLblklo95Z59ozhvEQg/7ZM/+Asn5aVa/hO74fCup3HQz7Pipl5+PK+lagqzYZORwisErX0DinzlCJNkwGfsUQwrzgDlflWvBiirDPgcOOFrfW4YM4k5NssYT13vObpCCHwrz0n/A4VnYxS0hmIcUnHMbyko/eWdOQxBnuOjT/gu9QHrwKis9molwK+98/L5RZ4YtMhPPDe/mHPI/ftE5Fy0jpYht/izfCjgQM+Y2NERLi2qhg76jux/0TPsPv/uaMR3YMu3LS4NOznjldJZ0d9J259phpvjPCtJZBc0ukaRUnnaFsf5v/0XaWVcjyG1fA9Hu+mrU4ZLd3cMzTu8wxyZi6NQB5e0tH7lXQ8+LJheClarwtengk2S6e1d0j5+480TQZ8ruGzRHHlgiIY9YQXtvpn+UIIPLP5KHSe2H0AABcgSURBVGZMsmJRaVbYz2sx6mE1G2J++GrrEanr6GBL76h/JpySTm1zLzr7nahtHv3zh+IL+NIHjjxLR68jvz74sYzBUFNv2gb24QcGfKdbBP07k6/EFViW15N/wHe5PWjrc8Ce7n+d40jRZMDnGj5LFDnpZqyYmY+12xv9+qy3He3AnqZufG1x6Zg33+LRi199RDq1Gk4GLpdyRhPwlQ+HUXwbaO9z4MI/bAwZsOW6edeAA0JI9XWjXgej3v/Pe7xlHaeqLVMO/vJfqdkgzbRXXwzlRPfwvzOj98Llge+EwBbS9n4HhEByZfiMJZLrFpWgvc+B9/eeUG57evNRWC0GXLGgcMzPmxvjgC+EUMYUhBPwe4aG1/Df3X0cP3rly2GPDefDYVdjF/Yd78E7u4dvgAK+Gr7TLTDgdMPtvaatXucLazaLYVwbt0L45uc4VdMylcsUGv0zfABo6hoY9jwGvZzh+4d8XUCG39Uv/blkpCZRhs9YIlk63Y5JNouyedvcPYh1O5twXVUJUk2GMT+v3WqOaZfOkbZ+tPU5kJlqxJG2vlFfek8O4j2DTuVn3tt7Ai9WNww7iRzOBm9DhxQ45Q+hQAOqrLprwAmnR0CvJxhV5ZN5JZnj2i8ILLfIJR35yloWg35YwFePTZbJM/CHZfg6UoanHWzpxdObj3ifNzqhWZMBn2v4LJHodYRrTi3Ghv0tON41iOe21MPlEfi3M6aM63ljPU+n2lu/v2J+EQadHhwf5XA4OVv3CKDXezq0vc8Bt0cMO4wVTr2/vkM6NLWjrjPoydnAgC+1ZeqUDVKr2YDsNBM6+8feLqrepHW6BY56D3KpM3wd0Um7muQyky4ww1fN7l/10Cb8/VPp3EA4ZzbCocmAzzV8lmiurSqGRwDPb63DPz47imUVdkzNTRvXc9qtZvQMuvymLNa19WPpbz6MSLthoG1HO2CzGLBipnQieLSZcY8qqMsZvHxCuLPf/6SwUtIZRQ1fzvB7hlw40Dy8C0r959LZLwV8aVqmFNasFiMyU4zD1hDK2u0N+Pkbe3C0zff7dqg2aXc2duGnb+wBoLpMoUEHg/7k8/flkk5gim/QkXIqt9/h+/3I3yAiTZMBn7FEMyUnDWeUZeORD2vR3DOEm84cX3YP+Hrx1V0fz2w+grr2fmw8EPnR4dVHO3DqlCxMy5M+qEYb8LsHnMj1dpXIAV0O+B0B2bWc2feMYgxDQ0c/SrKlYXPByjoDDjdsFoPf6xn1pGT4thQjMlJN6B50DfuG4PEIpVdf9tAHtXh802Hc/Netym2hrpxl0ssBXz8saw/GGKKkE2y0AuD7QIk0DviMRch1VSVwugVKslOwrCLv5D9wEoG9+INON17+XBrJvDvCGX5nvwO1zb2oKs1GvtUCi1E3qoAvhEDPoAtFWakAfAG9TQn4ARl+OCWd9gGcWZaL3HRT0IA/6HJjUoZ0oK3N+6Fo0OuUbNpqMSAzRTrgFDj+4o2dTTjz/g/8zg7Im8CHW/uUax0EG4d81YIilNnTAfhGK5yMb9M24HYdBd0rMSdTSYexRHTRnAKU5qTi28vKQx60CUduwGnbt3Y2obPfiXybGbsbI7u/JQfUU6dkQacjlOakjSrgD7k8cLg9KM6SMvHuASccLo+SwYcs6Zwk4A863WjtHUJJdgoWTs7CZ4fah+1nDDjcygnm1l75AuGkZNM2iwGZqUbvOvxfb19TNwacbhzrlMpGQgi09TmUMxPb6zoBIOg45N9eO0+psZuNOuhO8ndtsxhUm7b+j9WTlOEHHlpLqoDPm7YsEaWY9Fj/g3PCnpsTipzhy8HsH5/VoSw3DTeeNgWHWvv8Lu83GnubuvHDl79EX5Cfqz7aAYOOMK84EwBQZh9dwJdr8XLA7xpw+mX1HX3BSzon69KR6/fFWam4bH4hjnUN4Kxff+B30ZlBpwd2qxlEQFufN8PXEfTebNpmMSLDm+HLJ2+vX7MZv357H5q6pAxe/hDpHXLB4fJg6XQ7jHrC9nrvVbS83TdyFp9m0isXSgd8ffihvHT7Ymz97xXKpu2wg1feDH/eT971uz2pavi8acsYkOOti7f0DGFvUze2He3AjadPxtxiG4DwT5D++u19eKG6Hr94a++w+7Yd6cDsogykmKRAMzU3DfXt/Sedxy934RRnpii/Vl+sOzDDVyZrniTDlzt0SrJTcOkphXjvnmWYZk/HnzccVB4z4HQjzWSAzWJUavh6vU5py7T6ZfjS/bsau1F9pF3plZcDvvzzBZkpmF2YoWT48kXL5VHFaWZpz0AO8vKmbSgpRj3MBj0MIQ5e6XUU9IM7qTJ8xph0OjM7zYSW3kE8+1kdTAYdrl5YjNmFUiK0K4yyTm1zD9bXtKAoMwXPflaHD2ualfs6+hzY0dCJqim+ERClOWlweYSSaYciB+7CTF9Jxy/DV5VShBCqPnzniFcLU2f4ADDNno7rqopR29yrfPMYcLiRYtIjI8WofAsyqg5eWS1GZKRIH5pdA04MOt3oHXLhcGu/L8P31v7lPYecNBOqpmRhR30n2nqHlO4bi/eDMN27SSxH7pNt2sr3+UYrBJR0dBT0OgJJleEzxiT2dDOOtvVj7fZGXDq3AFlpJuRZzchNN2NX4+gz/Cc/PgKTQYeXbl+Mynwrfvjyl+joc+BgSy+uePRjAMAlpxQojy+zy506I8+8kQN4ZqoJVrMB3YNOJXjqyH/TdsDphssjkJtugkcAfY7QvesNHf0wGXR+UyPldtH3956AENLpWotRCvjqTVu5fGJL8WX4XQNOJZtv7R1CY4d/hi9/K8lJN2H1ohI4XB78/dM65RuOnOGnezN8uevHYtSNuF8j3xdqtEKon02qPnzGmCTXasJHB1rRO+TCV86Q9gaICHOKbNh9bHQZfkefA6983oAr5xehMDMFv189Dx39Dnzr79tw5SMfo3fQheduPR0LJ/sy/Km5UhfK4VaptPLWzib88f0Dw55bLtHYLAbYUozoHnCh3Rt8S7JT/TZL5bq90tEzQi9+Q8cAijNT/DZES7JTMWOSFe/uOaGcbLUYdchMNSqB26DztWVa1TX8fqdfe6vcCtnSM4RnP6tT/iyz00yYnm/FuTPy8MzmI8p8pMCAL7drysPTQpHvCtWHzwEfvGnLmEzOcGdMsvoF5DmFGTjQ3Ot3+AgANuxvweJfve831+fZLXUYdHrwjbOnAgBmF2bguysqsOVwOyZlWPDqnWfh1CnZfs+TlSoFy8OtvTjWOYDvvfgFfv+v/dhR3+n3OLmkY0sxwpZiRNeAE+19DhBJZSF1hi8/VunoGaGO39Dej6Ks4Rd8P39WvrcGL5VkUox6VOZblbKIwTseGZA2bY16HdJMem/AH34Aa0d9J368diceXS/tDeSkSX/eyyvtaOtzKIPQ5L0NuYYvZ/5m48ibtvIHltylE1j+CfWzkejyCrqeqDzrOPGmLWMSuVPnK2dM8av/zi60we0RqDnufwL12c+OoqlrEN98php/3nAQDpcHz2w+giXTc1E5yao87vZl0/CnryzE/377TJRkpw57XSJCaa7UqfOzN/ZAQCAjxYgHAy7sIWftNosRNouvpJOVakJOwFgDpaNHqfcH79TxeAQOtfZhcpB1LS7LgUf4xkCkGPVYUmFX7jeoOmis3np7ZqoJnQOOYWOLi7NSUNcufYNxuDxINemVwC63xMqbu8MyfNUQNf0Im7ZygDcZQnTpjPCz0aDJgM8Yk8wtzkRRZgqumO8/dXNOkXfjVlXW6Xe4sGF/C66rKsbFcwtw/7p9uOpPH+NE9xC+cdZUv5/X6wgXzS2A1XvlpWDKctOw9UgH1u06jrvOKcdtS8vwYU0LvlBl+T2DThh0BItR5y3pSBl+dpoJmammkTN8VUlnb1M3fvduDYQQ2NPUjZ5BFxaV+n/rAKActJIDdYpJj9NKs5USiEHnf/AKADK862r1ln3kMo/cgirLTvNNqJQ/aI91St8kLN6TrxZv4PcFfD3STKE3WPXkn+FbAjZjR/p2EA0c8BnTsMvmFWLTD88ZFpiLs1Jgsxj8TtxuqGnBoNODKxYU4eEbFuCe8yuwq7EbZfY0LFNlwaM1NTcNDpcHZblpuHVpGW46sxSZqUY8qKrldw86YUsxKpft6xl0oc0b8LNSjeh3uJXBYnJGL3feyGOVAeDJTYfx0Ae12NnYhU8PtQEATi8bHvDzvAet5CFmFqOUlZ8+VXqsQe/L8G3ewJ6ZakRnv/TNw2oxoCI/HSa9zu8bDyB16MjsgRm+N6jL7ZK+qZm6oB+agX338odQmtk/4I/mlG4kccBnTOOCXUCFiDC7MMPvxO3bu48jO82E00qzQUT4znnT8dLti7Hmq6ee9DRoMLMKbCACfnL5bJgNeqSbDbh1SRk+2NesZPk9gy5lno0txSC1ZfY5kJNmQqY3gMoz3odn+L6SzmZvkH9zZxM+PdSO0pxUFGQMr+GnmfRIMepxtN0X8AFgyfRcAFImLbc0ypl8RooRnQNOtPRK14pdODkLMwttyLdJQX1eiZTp56g6guQMX94rkFs9h2X4Rp3yTUJN7spRunS8Py+XhGRj+XsZj7EP62aMxdWcIhue3nwUTrcHHiHwwd5mXHJKgbJpCSBoWWS0zpuZh833nqeUUQDga4un4A/v7ce6XccxryQT3QNOJcO1WYzoGXJB3zuERVOzkeVtiezodyLPZlFKOEUBJZ369n40dAzAoCO88UUTegaduHhuAYIhIuTZzKiXSzreAHzJKYVYt+s4puenY2aBFelmg1KHlzP81p4h5Kab8Z8XzoBHCGw60AoAWHVKAera+vwuNJ9mNiDFqFcyfI+qDROA31z84Bm+DoDb14evZPj+ITfWGT4HfMYS1JyiDDhcHhxskTppeoZcuGDOpIg9PxH5BXtAanWcZk9HzXGplNQ96IItxVcrB6QWyJw0E7K8V22S6/jdgy5YjDqkmqRgKmf8cnb/9bNK8dhHhwEAZ5TlhFxXntWslHTkgF+UmYK1d5ylPObqU4uV/85IMUk1/N4hVE6ySte9BWFGgRVFmSlYWmHHWeW5yslmmd1qVvYK5BO3cobvVm3aBsvw5T0FOZ7LGX9gwI91hs8lHcYS1OxCacTCrsZuvL3rOKxmA86alhv1151ZYMM+b3dQz6ATNjnDT/FlutKmrf9Yg+4B9WMNSkln88E25KSZcOc55UrtO1j9XpZn9X0IpZhOHsIyU41wuD2o7xhQsn4AKMhIwcf3nouKfCtmFtj8nhfwlXX0OlKyeLmGL5/ANehJKWmpyeOT5Xn6ciYfWNLhGj64D5+x0Ziam44Uox5f1HfiX3tO4LyZeVE7sKM2Y5IVTV2D6Ox3oHvApWS46sCX7Zfh+2r48oeCzWJE96A0XmHzwTacMS0HmakmnFOZh4r89KD1e5n6At9yxj2Sud6OJofL4xfwT0beuJ1TaFMCc2CGb9QH37RdXiltksvfQOQy27AaPnfpcB8+Y6Oh1xFmFdqwdnsjOvqduHBO8Lp3pM0okL5Z7DveIwXxIBl+Tpp5eElnQL3BKwX8I239ON49iMXeEs7vrpuH5249Y8TXz7OFF/DPKs/Ft5aWAfD/sDgZ+YPs9LIc5YCbkuF7SzxS9j88w7/vstnY8IPlykawMUQNP9a4hs9YAptTaMO2ox1IMerH1Ho5FjO97Yy7GrvQ73D7Ze2y7DQTUkx6mA065fBV96BT6XW3WQxo7XXgk4PSxuniaVLAH+lcgMyvpDOKgA8A/3nhDEyzp+OC2aPf42j0zspfODkLDrcHr+44hlOKpSRUyfB1odoydZiS47vEpTwnLjWgZ189IkK++LkxioexOOAzlsDkyZnLK+1Kr3i02a1mZKeZsNV72lU54JSqyvC9G6BZqSblUoLdA06UeoOg1WLEodY+bDrQioIMC8rCuP5vXpglHUDKxK9bVDLq1wCkjqRPDrbhjLJsZKaasHJWvvJ6j32tCn/bfBTFWSloH8U1c+WzCIFjj4e83xwWlGRh86E2/OCCStx4WmSupxAMB3zGEtjCKVkgAi4POIkbTUSEGZOs2HpEukiIUtJRlTbkDdvMVKOqhu/r6LGlGNDR58AnB9uwclZ+0LMGocglHdNJBpeN14VzCnDk/kuUX6s/XCryrfjZFXMAIGhJJ5CcyQcG/EHv7csr7bj/6rmYnJ0a1p9FuDjgM5bAyvPS8cm95464yRkNMybZ8MlBqZ1SLumkmQzQkfRv+fBTVqoJnf0OZRa+zaLetJW6dM6eHl5nkVzSGW05J9oCZ9fL3VNqyh5AwJrl2y1GvV8JKFo44DOW4GId7AFgRoFvLIGc4eq87YuZqtJOZqoR+0/0KLPwlXq/aoP3rPLwAn5mihFGPSmHoLTk1TvPwtQggVtdq1fzBfzY/F444DPGwjZzki+LVW/W2lIMfkPIMlOliZnqqZrqf88ssIXVKglIHyy56eaoXQZwPOaXZAa9fcjpO5mrJm/UjnYvYrw44DPGwjY9Px06AjwCSl0eACZnp/qNKMhKlebYdA3Ic/MNfv9eEmY5R5ZnNStZcyIItWn7nfOm40hbH5ZX5sVkHRzwGWNhsxj1mJqbhoMtfX5tiX/5apXfyN+sVBPcHoEXq+sB+DJ7+UNBPqAUrmUV9qDXgo2X+1bNgnGEbxz/Z9Vs6HV7h5WvyvPS8dpdZ0d7eQpNBnwiWgVgVXl5ebyXwhgLYUaBDYda+2BVHSYKPEk6o8AKHQFPbJJm5BRmSoG+akoW3vzO2UpbabjuWVk5xlVHx80B1xsINDU3DY/fVBWj1YRGI105Pt6qqqpEdXV1vJfBGAvi3d3H8fau4/j96vkjPm7Q6UZ7nwMC0pAzFl1EtE0IEfTTRZMZPmNM+1bOnoSVozi5ajHqUciBXhO0t83NGGMsKjjgM8ZYkuCAzxhjSYIDPmOMJQkO+IwxliQ44DPGWJLggM8YY0mCAz5jjCUJTZ+0JaIuAAdGeEgGgFBXOs8F0BrxRUXfSL8nLb/WeJ4r3J8d7eNH87iRHsPvL+281kR8f53s/rG+x6YIIYIPKRJCaPYfAGvGej+A6nivPxq/Z62+1nieK9yfHe3jR/O4k7yH+P2lkdeaiO+vk90fjfeY1ks6r4/z/kQUy99TJF9rPM8V7s+O9vGjedxIj+H3l3ZeayK+v8J5rYjQdElnPIioWoQYIMTYePH7i0VbNN5jWs/wx2NNvBfAJjR+f7Foi/h7bMJm+IwxxvxN5AyfMcaYCgd8xhhLEhzwGWMsSSRFwCeiNCJ6mogeI6KvxHs9bOIhojIieoKIXo73WtjERERXeGPYC0S0cizPkbABn4ieJKJmItoVcPuFRFRDRLVEdK/35qsAvCyEuBXAZTFfLEtI4bzHhBCHhBC3xGelLFGF+R571RvDbgeweiyvl7ABH8BTAC5U30BEegCPALgIwCwANxDRLADFAOq9D3PHcI0ssT2F0b/HGBuLpxD+e+y/vfeHLWEDvhBiI4D2gJtPA1DrzbYcAJ4HcDmABkhBH0jg3zOLrTDfY4yFLZz3GEl+DWCdEOLzsbzeRAt+RfBl8oAU6IsAvALgaiL6EybmcXkWO0HfY0SUQ0R/BrCAiH4Un6WxCSJUHPt3ACsAXENEt4/liQ3jX5v2CSH6AHw93utgE5cQog1SbZWxqBBC/BHAH8fzHBMtw28EUKL6dbH3NsYihd9jLNqi9h6baAF/K4DpRDSViEwArgfwWpzXxCYWfo+xaIvaeyxhAz4RPQdgM4BKImogoluEEC4AdwF4B8BeAC8KIXbHc50scfF7jEVbrN9jPDyNMcaSRMJm+IwxxsLDAZ8xxpIEB3zGGEsSHPAZYyxJcMBnjLEkwQGfMcaSBAd8xhhLEhzwGWMsSXDAZ4yxJPH/AUuXJkqY+5SXAAAAAElFTkSuQmCC\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# Orthogonalising the basis"
      ],
      "metadata": {
        "id": "z2AlBeXY60FX"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "residuals2 = []\n",
        "\n",
        "for m in range(1, Niter):\n",
        "    V_m = np.zeros((N, m))\n",
        "    V_m[:, 0] = b / np.linalg.norm(b)\n",
        "    for i in range(1, m):\n",
        "        v = A**i @ b\n",
        "        for j in range(i):\n",
        "            v -= np.dot(v, V_m[:, j]) * V_m[:, j]\n",
        "        V_m[:, i] = v / np.linalg.norm(v)\n",
        "    # print(np.linalg.norm(V_m[:, -1]))\n",
        "    b_m = V_m.transpose() @ b\n",
        "    A_m = V_m.transpose() @ A @ V_m\n",
        "    y_m = np.linalg.solve(A_m, b_m)\n",
        "    x_m = V_m @ y_m\n",
        "    residuals2.append(np.linalg.norm(b - A @ x_m))\n",
        "    print(residuals2[-1])"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "UXv92Xo465sK",
        "outputId": "1d35a012-56e4-4f39-d0e2-dc21f87ea8ec"
      },
      "execution_count": 7,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "2.543606159277365\n",
            "1.5268525024832802\n",
            "0.7909911169132886\n",
            "0.4696176583152166\n",
            "0.2413604073795114\n",
            "0.12512897559019778\n",
            "0.060891748931758964\n",
            "0.03298433313613242\n",
            "0.01746159978324473\n",
            "0.00951152946224309\n",
            "0.004664271554544925\n",
            "0.002467109299810796\n",
            "0.0013873620572100263\n",
            "0.0007026934259215055\n",
            "0.0003475143611268963\n",
            "0.00018605783456961366\n",
            "0.00010247231876058544\n",
            "5.4015093113319406e-05\n",
            "2.5029171385935433e-05\n",
            "1.378869922106686e-05\n",
            "8.373214226616235e-06\n",
            "8.756998135159632e-06\n",
            "6.3045964326279245e-06\n",
            "5.3022346681989085e-06\n",
            "4.2139709328365024e-06\n",
            "3.1175673459982686e-06\n",
            "3.2390429860032925e-06\n",
            "2.7501197239099114e-06\n",
            "2.5174201779388734e-06\n",
            "2.518058816950447e-06\n",
            "2.261635050840596e-06\n",
            "1.9883625161188013e-06\n",
            "2.0147727199586173e-06\n",
            "2.0621679746011067e-06\n",
            "1.912919778281111e-06\n",
            "1.7199806362682747e-06\n",
            "1.58727039377306e-06\n",
            "1.5406818101670758e-06\n",
            "1.5266071926936856e-06\n",
            "1.4853911214713953e-06\n",
            "1.4721716475478571e-06\n",
            "1.5065439497433622e-06\n",
            "1.4931803039115955e-06\n",
            "1.4525001353958056e-06\n",
            "1.3188409935265844e-06\n",
            "1.1947308048483045e-06\n",
            "1.1905974510018879e-06\n",
            "1.1629295673937046e-06\n",
            "1.1633516062367977e-06\n",
            "1.0738226311482194e-06\n",
            "1.0411769227053313e-06\n",
            "1.0533930552334082e-06\n",
            "1.0572800330267503e-06\n",
            "1.0583799864005707e-06\n",
            "1.0583719972756822e-06\n",
            "1.0073814027553532e-06\n",
            "9.838459582829616e-07\n",
            "8.729585970136393e-07\n",
            "8.751426091855014e-07\n",
            "8.792979274790214e-07\n",
            "8.716524329701659e-07\n",
            "8.606928790612869e-07\n",
            "8.587013499844685e-07\n",
            "8.493326992344008e-07\n",
            "8.423371180967301e-07\n",
            "8.136266106483476e-07\n",
            "8.112913568304705e-07\n",
            "8.063708300636296e-07\n",
            "7.522466564749314e-07\n",
            "7.061425100013603e-07\n",
            "6.71593481485072e-07\n",
            "5.390396306699317e-07\n",
            "5.226082226509487e-07\n",
            "5.219812370415829e-07\n",
            "5.224100670405016e-07\n",
            "5.03622965816725e-07\n",
            "4.454997369486936e-07\n",
            "3.7787658596618824e-07\n",
            "3.76907130865984e-07\n",
            "3.77036895106751e-07\n",
            "3.653589834267529e-07\n",
            "3.636006864847608e-07\n",
            "3.40158140608746e-07\n",
            "2.965562608224648e-07\n",
            "2.935051810321643e-07\n",
            "2.9089053582986196e-07\n",
            "2.777397364609529e-07\n",
            "2.7106851966444936e-07\n",
            "2.489357091295347e-07\n",
            "2.45981196884943e-07\n",
            "1.246210214326445e-07\n",
            "1.057093945497795e-07\n",
            "7.429365777218514e-08\n",
            "7.14533128449237e-08\n",
            "7.044169126338894e-08\n",
            "5.913766579960124e-08\n",
            "2.517280171832061e-08\n",
            "8.262043688257606e-09\n",
            "4.791084107529659e-09\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.plot(range(1, len(residuals) + 1), residuals)\n",
        "plt.plot(range(1, len(residuals2) + 1), residuals2)\n",
        "plt.xscale(\"log\")\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 273
        },
        "id": "G3jtTwxC7Mp3",
        "outputId": "feb9e7b8-4260-474b-ff88-13ebddac292d"
      },
      "execution_count": 8,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5xU1d3H8c+ZmZ3tvQFbYOkuHRaWpiKKIrbECmo0imKJicljnjyaHjWJT0yMjy12MdijqIgo2EGpC0hvS1nYZRts7zuz5/njzg6zDbZN2d3f+/Xa1zJ379z7GzP5zplzzj1Xaa0RQgjR+5m8XYAQQgjPkMAXQog+QgJfCCH6CAl8IYToIyTwhRCij5DAF0KIPkICXwgh+ggJfCGE6CM8FvhKqcFKqZeUUu966pxCCCFO6VLgK6VeVkoVKKV2Nts+Vym1TymVqZS6H0BrfUhrvbAr5xNCCNF5li4+fzHwFPDvxg1KKTPwNDAHyAY2KaWWaa13d/TgMTExetCgQV0sUQgh+pbNmzef0FrHNt/epcDXWq9WSg1qtnkKkKm1PgSglHoLuAJoV+ArpRYBiwCSk5PJyMjoSolCCNHnKKWyWtvujj78BOCYy+NsIEEpFa2UehaYoJR6oK0na62f11qnaa3TYmNbfEAJIYTopK526bSb1vokcKenzieEEKIpd7Twc4Akl8eJjm1CCCG8yB2BvwkYppRKUUpZgfnAMjecRwghRAd0dVrmm8A6YIRSKlsptVBrbQPuAVYCe4B3tNa7ul6qEEKIrujqLJ0FbWxfAazoyrGFEEJ0L1laQQghPKSoso680hqvnd9js3SEEKKvm/TwZ2gNRx65xCvnlxa+EEJ4iNbePb8EvhBC9BES+EII0Uf0yj78g9vX0mC3ERk/kMi4BMyWXvkyhRCiQ3plElau+C1jazYDUK/NFKpISi0xVAbEURfUD0L7Y4lIICgmmfD4ZKL7DyIgMNjLVQshhHv1ysAPu+JvbD22l7riHBpKc7BU5hFYk09M1UGiyzcSXNByWlQJIRSZYii3xlIbGI89dADm8AH4RyUSGmt8KIRFxqJM0gsmhOiZemXgDzorjUFnpbX59/LSIopyj1BWkEX1yWzspTmYynPxr84npLaAhJoDxBSXtHhejfbjhCmaUr9YqgPisQXFQ9gArFGJBMUkERE/kJh+yVj8rG58dUII0Tm9MvDPJDQ8itDwKBg5sc196mprOJmXRWn+USpPHKW+OAfKjuNXlU9QTT4DyncQU/o11jxbk+fZtaJQRVBiiaHCGktdUD8aQvtjCU8gMCaJsDjj20JwaIS7X6YQQjTRJwO/Paz+AfQfOIL+A0e0uY9uaKDoRC7FeVlUnDhKzclsdOlxTBW5BNQUEFmTTVTVNsJOVLZ4brkO5KQ5hnK/WGoC47EF98MUnkBo0mgGjzubgKAQd748IUQfJIHfBcpkIiougai4BGB6m/tVVZRyMvcIpQXGh0J9yXFMZTlYq/MJri0gvuQI0cXFmHM07Ia6T83s9RtOScxEAobMZNCE2UTE9PPY6xJC9E4S+B4QFBJO0LBxJA0b1+Y+tvo68vOPkbtvE9WZa4g8sYWJuW9hzXsdvoMsUxL5EeNRydNIGHc+/QcOlwFkIUSHSOD7CIuflfjEIcQnDoHz5wNQU1XB7u3fUrpvNUF5mxhZ9AVhRR/B91BAFMdCxlKfkE5M6rmkjEqX6w2E6CEaGjQmk/L4eSUhfFhAUAipU+fC1LkANNjtHN6bQcGubzAfW09S+ffE7/sa9v0vFUsDORSYSmX8ZEKHz2TI+FkEBod69wUIIVpla9BYJfDF6ZjMZlJGpZMyKt25Le/oAbK3fYn9yFriir9n9JHnMGU9S/0qM/v8hlIcPRH/ITMZOP48x1iDEMLb7A3eWUVNAr+H65c8jH7Jw4A7ACgtKuTI919RlfktEYUZTMz7D9b8N2EtHDUlkBc+HpKnMWDMeSQMTpVxACG8wNbQAJg9fl4J/F4mPCqWcbOvhdnXAlBTXcne7d9RvG8NgbkbGVH8NeHFH8M2YxwgK3wypJzLwMnziEtI8XL1QvQN0sIXbhEQGMzI9Ash/ULAGAc4sn8r+Tu+wu/oGoaWriXy+5Xw/a/JMiWSFzUFv2GzGTJ5LuFRsV6uXojeySaBLzzBZDa7LD3x3zTY7RzctYHC7SsJyv6OMYUfE3RiKfa1igN+QzkRO42Qs2YzLG2OXAwmRDeRFr7wCpPZzJCx0xky1rhwrK62ht1bv6Z01+dE5K0l7fjr+OX+m9ov/Njln0rZgOlEjp7D0HFny5pBQnSStPCFT7D6BzSZClpZXsLujM+o3vsFsSfWM+rIv+DIvyj/KJCDweOpSZxJ/PiLGDRykgwAC9FOdrsEvvBBwaERjDvvGjjvGgCKCnI4nPEptsyvSSjeSOL+dbD/UbJMSeSNvJkx8xYRFBLu5aqF8G3GLB3Pk8AXHRIVl0DUvIXAQgBys/ZxdONyova+Qfruhynb/U/W9/sByXPvZcCgtheeE6Ivkz580SM1riiqG37B3owvqFz9FGm5b6JeeYOtITOwzvgJqVPnSnePEC6kD1/0aMpkYuSUOTBlDnnHMjn8yROMPL6UyFULOPhFCkWjbmHM3IUy00cIvNfCl2aX6Hb9koYybdETBP5qLxvH/AmT1kze/nuq/zaSdS/cS372QW+XKIRXeauFL4Ev3CYgKIQpV/2cQb/dyq45b3AkaCxTsl8l+oU0Nv/9CvZu/AztpcErIbzJLoO2ordSJhOjZlwCMy7h+OG9HP30cVLzPyRsxdUcWDWM8sn3MmHODdLPL/oMm5emZcr/w4RHDUgZydS7nsV83x42pP6GAHslE9fdw4G/pLPjm6XS4he9ltanQl768EWfEhwaQfq1v6L/r7excdxDhNmKGfPVLex+5Bz2bljl7fKE6HauGS99+KJPsvhZmfLDnxF5/3Y2jLyf+LpjjPzkGrb97xwyt33n7fKE6DaurXpp4Ys+zT8giPT5DxD0y+2sG/wzBlXvYuj789jy98vJ2ve9t8sTossaXLp0pIUvBMYN36fd9BDcu531iQsZUb6BxDdmsenx+Rw/ss/b5QnRaU1b+N4Zq/JY4CulBiulXlJKveupc4qeKzwyhqm3PUbtT7awqd98xhZ/Tswr09jw1C0UFeR4uzwhOszeU1r4SqmXlVIFSqmdzbbPVUrtU0plKqXuP90xtNaHtNYLu1Ks6Hui4hKYetezlNy+ga3RlzCx8EN4ZipbVi7xdmlCdEhDD+rDXwzMdd2glDIDTwMXA6nAAqVUqlJqjFJqebOfuG6tWvQ58YlDSP/ZErKvW0mROZaJ6+4h47GrKC0q9HZpQrSLa8j79Dx8rfVqoKjZ5ilApqPlXge8BVyhtd6htb602U9BN9ct+qiU1MkM/J91rEtexLjSr6h7YjLbvnzH22UJcUb2Hj4PPwE45vI427GtVUqpaKXUs8AEpdQDp9lvkVIqQymVUVgorTfRkp/Vn2m3PkrWlcuoMIUxbvXtbHx8AWUlJ71dmhBtch2n9ek+/O6gtT6ptb5Taz1Ea/3X0+z3vNY6TWudFhsrN9EWbRs6biYDfrWedQNuZlLxJ1Q9PoUdqz/0dlnCiyprbZRW13u7DACyTlbywupDzsdNW/g9b5ZODpDk8jjRsU0Ij/EPCGLaoifIvGwpdSZ/xnx5ExuevJnK8hJvlya84Hcf7uSu1zZ7uwwAlm/P5c8r9lBZawOaDtr2xBb+JmCYUipFKWUF5gPLuqcsITpmRNps4n65kfXx85l84kNKHpvC7nWfeLss4WHHS6rJLa3xdhkA1Nbbjd82ozXf0FP68JVSbwLrgBFKqWyl1EKttQ24B1gJ7AHe0Vrvcl+pQpxeQFAIU+96jr0XvwVA6sr5bHh6IVUVpV6uTHhKVZ3d2aL2hMyCcv66Yk+ThdEa1Tlm4tTajOC3N2vhbzlazGe78z1TqEN7Z+ks0Fr311r7aa0TtdYvObav0FoPd/TL/9m9pQrRPqlT5xJ53yY2xF5NeuG7FP9jMrvWrvB2WcIDKmptVNXZPXa+lbvyeW71IUqqWo4b1Dla9rX1rbfwr3xmLbf/O8MzhTrI0gqiVwoKCSf9Jy+x68I30ShGrVrAhqdulb79Xq6q1k5lna3VFnd3agzzCse3iRpbyw+ZOnvTLh27yzitT3fpCNFTjZo+j6j7NrI+7lomFy6l9LEp7PruY2+XJdzECHuoqXffLJj9+eWk/v5TDp+opKox8Fs5n7OF30aXTiN3fzi5ksAXvV5QSDhT736BvRe/RQMmRn12PRueukVa+72M1trZnVNZ575+/CMnKrE1aHKKq6moNc5XU99KC98Z+K116Zz6gKj0YBeUBL7oM1KnziX6l5uMmTyF71P62BS52UovUmtrcLaiq2rdF6KNHyo19acGiFsL/PrGQdv6xi6dU4Ff6/KNwJPXDUjgiz4lMDiUqXc9x/5579CAiSEr5rP+jYfl1oq9gOvsHHe28F377RvP01qXTm3zLh2XFv7Jyjrnv8sk8IVwr5HpFxJ273fsDJ7K1P2PsuWfV0kXTw/nOjunyo2BX+US8qcftG3WpePSwi8sr3X+WwJfCA8Ii4hm3H0fsS7lHsaXfUXhP8/m2IFtXqsnr7SGy578lle+O+y1GnqyCtcWvhu7dBqPXe3SpVPbah9+4yydloO2BeWnLg6TLh0hPMRkNjPt5j+z+/xXCW8oIeK1i9i66jWP11FYXsv1L65nR04p72+VFUo6w7VV350tfJu9gRdWH3L207uGfKVz0PY0s3Qa+/BdunSOFVU7/11W47kLxSTwhQDGnHMFtbd+Sa5fEhPW/oR1z/8UW33dmZ/YDUqq6vjRSxvILalh9sg4duaUUlbjGwuA9SSurfrubOFvOVrCn1fsYe3BE8axXQZtK04zaNuyS+fU36pd9pcWvhBe0C95GMn3fc2GqMuZdvzf7P37HIoLc916zrKaem56eSOHTlTywk1p3HZ2Cg0aNh8pdut5eyN3tfDLHR++5Y6WeKXL3PvTztKxNVtawdHCD/QzN9lP+vCF8JKAwGDSf7aETWMfZFjNLqqeOZfDu92z+mJlrY1bXtnEntwynr1xIjOHxTAxORKr2cT6Q7K2f0dVuLbwu3Fue2MrvvFbQ+OHSVlNvfMCqhpbK1069mZLKzj2jQq2AhDibyHE3+LRb3MS+EK0YvKV95J1+TsE6Fpi3r6ELV907121qups3PZqBluPFvN/8ycwe2Q8AAF+ZsYnRbD+cPMbzIkzadLC78YF1Jq37BuD/2TFqS6/mno7aw4U8s3+UzdtauzD/8dn+/njsl3OQdvGwI8KthJoNbf67cBdJPCFaMPwSbOxL/yCQks/xq1exJolD3bLfP2KWhs/fnkTGw6f5LFrxzNvTP8mf08fHMXOnNIms07EmTUGscWk3NLCd7b0HR8sJypOTa2sqW/gyS8yefzz/c5ttS6t/sVrjzivtHUN/AA/k1uXgWhOAl+I04hLGsqA//qGnaEzOPvgP1j3fz+ipqb6zE9sQ1lNPTe9tIHNjpb9Dya0vCvo1MHR2Bs0GUekld8RlbU2TAoig63d2odf0aKFb/x2vXiqcQDX9eKvumZz8xtn4zQGvr/FRKCfmWpZWkEI3xEQHM7Y/1pGRvKtTC9dzoG/z6Eg/3iHj1NaVc+NL25gR04pT18/kcvGDWh1v4nJkfiZFesPSeB3RGWdjWCr0S/enbN0mrfsGy/wOunSwq+12amqszU5b+PSCo2+P2YMxDcGvsWsCPAzt3rRlrtI4AvRDspkJu3Wf7J9yt8YUb+H2n+dR+a+He1+flFlHQteWM/e3HL+dcMk5o7u1+a+gVYz4xIj2HBYBm5Px2Zv4O7XN7P1qBGkVbV2gvzNBFnN3TxLp7FLx+74bTwudlkDv6a+gco6e5NuuDp7066aHdnGjXgaA9+kFAEW6cMXwmeNnXcHx3/wLmFUYH3zKvYeOHDG5+SX1XD9C+s5WFjB8zdN4oLU+DM+J31wFNuzSz1696bOsjdoZj36Fe9kHPPoeXNLa1ixI4+PtxtTZyvrbAT7Wwi2Gi38qjpbt6w7X1FrBHtlra3JipyNlDq1kFrjPvYG3eLcJxyDvP4WI3YtJkWA1Ux1sz783cfLuP+97eSUdL7rsC0S+EJ00KAJ51F9zVvEUIp6/Sp2ZGa1uW9mQQVXPrOWo0VVvPzjycwaEdeuczT242/O8v35+EeLqjhysorvMk949Lz5ZcbyBHvzygEjkIOtFoL8zVTW2bjwn6t5aPnuLp/HddDWdUXORtHBVqrr7VTV2bE1aGptDc7uG1eN0y/DAv0AGBgdTIDF1GJZhm/2F/LWpmPOD4buJIEvRCf0G3U2lT94hcHkUL/kGrZktlwOYXNWEVc/u5Zam523F01jxtCYdh9/0sBILCbl8fn4pVX1XPbkt6zY0f4Lzg7kG4G7P7/ijPvW1Nu5540tZBaUd7rGRvllRh/63rwywJh7H2Q1E2y1kHWyiuzial7fkNXllrLroG1r37gGRARS7DKA+9nufK7617oW+zV2DU0bHM2TCyZw/8UjW52WueHwSYbGhRAT4t+lulsjgS9EJ8WOn0flvGcYr/ZTvuQG1h/Ic/5t5a48rn9hA5FBVpbeNYMxieEdOnaQ1cLYxHA2eHg+/t9W7mVHTinPfnOw3c/JLDSC/mBhhbP1e7ykmhfXHGpxN6ddx8tYvj2XT3fmtThOcxW1Nmb+75e8veloq3/Pc7TwT1TUUVheS5WjSyfIanYuV1Bv1zz9VWa7X0tryl1m5zR25yhl/C3YaiY80K/JjJ23N7Xs2ooLPRXeZpPisnEDCPAzE2AxO5dZMLqLbGQcKSY9JapLNbdFAl+ILoiYch0VF/yNc9VWCpfcypr9+SxZn8Vdr23mrP5hvHfXdJKjgzp17PTB0Ww7VtJkANJmb+CtjUebLK/bXbYeLeaNjUdJigpke3Ypu46Xtut5mY6WfZ2tgayTlQAsWZ/Fwx/v4dCJyib7Nn4bOFBw5m8DX+zJJ7u4moc/3tPq6y0oO7Xi5L68cipr7UYfvr/FuX3m0BiWfX+8ydLEbTlWVMV7m7M53KzmxhZ+cVU9D39sdBFFOLplokP88beYm1yE9W0rXVsDIgKd/zablPPfrvPwl6zPIvX3K6motTEsLuSM9XaGBL4QXRQ2cxFVM3/NZabvOLzkHn73wQ5mj4znzdunOmdkdMbUwdHYGjRbsox1+mvq7dz52hbuX7qDPy7b1V3lA8YHyW/e30lcqD9v3j4Vq8XEWxvbNwh7oKCCaMfrbOzWaZw5szOntMW+rvudzkfbcokKtlJb38CDy3e3WIIgv6yGsAAj3PfmlTn68I1ZOmAsXXDJ2P5U1No4VlzlfJ7WmqMnq1jw/PomHyQPLt/Nff/ZxtzHVzs/uOBUH35pdT0rd+UDEB8WAEB0iHHxVPUZZtqEBpz6EDIpl8B3dOnY7A28uznbuT2yC++b05HAF6IbBJ3/K2rS7uIm8yoWD/qCZ2+cSKDVfOYnnsakgZGYHf34jYusfbE3n7SBkXy8I7dFmJ7J0i3ZnPO3r3jqywMtpi2+ui6L3bll/OGyUSRGBjFvdD8++D7njBcFNTRoMgsquHCUMfNof3459gbNdscUxOY17ne08F27f1pTVlPP6v2F/GB8AnecO5iPth1n7B9XcceSDGyO6Y55ZTWM6BdKbKg/e/PKqaqzE2Q91cIfEhfCqAFhgDHzBeDnb23l9n9v5r0t2aw7dJIv9xoBrrUxQH72sBjMJsVDy/c4Z9tU1dmbtMof+sFozh0RCxgDtgF+bf/vnBARyNjE8CaB36SFbzFTa2tg2iNfOv+bGcft/v57kMAXonsoRcC8v8D4G5iV9zKWDc90+ZAh/hbGJITzxd4C5j+3ni1ZxTx+3Xhe+vFkwgP9eOyz/a0+r6CspknfeUOD5h+r9vFf72zDZm/g76v2M+vRr1my7ghPf5XJ5U99y0PLd3Pu8FgudlwfMH9KMuU1Nufg7ZETlXy6M7dFn/zx0mqq6+2MSYggKSqQ/fnl7M8vd/Z178wpa7J/ZkEF/hYTdbYGjhVV0ZbPd+dTZ2/gkrH9+cUFw1l8y2Rum5nCyl35PPLJXsfrrCUuLICR/ULZfbzMMS3TTLDjg3ZobAjD40MxKdidW4a9QfP5ngI+35PvbE2vO3iSfXnl7Mgppaiyjrmj+3Hv+cP4fE8+972zzblSZqzLAOqN6ckE+RkBHhFktPDb8sfLR7HsnpkEW10C36WF39goaN5l1ZVvhqdjOfMuQoh2MZngsiegthxW/Qb8Q2HSzV06ZPrgKJ775hABfiZevDnNOa1z0TmDeXTlPjZnFTNpYKRz/5e/PcyDy3czIDyAOanxXJAaz9ubjrF8ey7XpSXx0A9Gsz27hL+s2MPvPjS6hcYlRfCruSO4IX0gyhFG6SlRpMQE8/qGLI4VV/HMVwepszdwYWo8j149jvAgow+7sYtmWHwIw+NCOZBfwdajRhfUzKExbMsuQWuNUoqymnpyS2u4aFQ8K3flsz+/nEExwYDjm0JhBcPjQ9Fa8+7mbBIiApmYHIFSilkj4pg1Ig5bg+bFbw8za0Qc+WU1zBoRx8j4UP7h+PAL9rcQ5GjhD40LIcDPzJDYEPbklrEnt8zZPZNTUo3FpPhmfyErduZhcbS6JyZHMrJfKKXV9Tzz9UFnSz4+PIC8shrCA/1QSuHvCPkQf4vzua6Mi7/s+JmVs65GJpfPh4A2pl5Gh0iXjhC+z2yBq16EoRfAR/fCzve6dLgrxiUwLjGc12+b2mQO/y0zBhETYuUfq/Y5t63clcdDH+9mxtBoRiWE83bGMX700kaWb8/l/otH8shVY7BaTKQNiuK9u6az9O7prH/gfD78yQzunjWUcMdAJIBSiusmJ7HlaAmPf36AuaP78au5I/hybwGXPLmGPblGy71xwHZobAjD4kM5dKKCTUeKiAq2csnY/pTX2DjqaMlnOj4cGheLcx24feTTvVz4z9W8/O1hlm7JYe3Bk9w6M8X5AdTogXkjCfQz8/7WHCrr7MSH+XP7OYM5q7/RdRPsmJYJRuADnNU/jN3Hy5wznkb2CwXgxqkDKa6qp87WQFWdnRB/C8PjQ1FKcfvZgwH4fE8BAP3CjBZ+jCOI6x0LowVZzcSGtux+GRYf6vzvCKf68C0mRZBLa7+t7qCIIL9Wt3eVtPCF6G4Wf7h2Cbx2FSxdBNYQGH5Rpw6VOiCMD++Z2WJ7kNXC3bOG8uDy3azNPEFIgIV739rK2MQIXrxpMoFWY1GubzNPEBVsbfItAIwgmpgc2eK4rhZMTuZAfgVXjB/AOcONlu60wdHc/u/N3P/edj74yQwOFJQTE+JPZLCV4fEh1Ns1q3blkT44mjEJxlTUHTmlDIwOds7QGZ8UwYDwAOfjTUeKeGHNIaKCrTz08W4C/cykDYzkx9MHtajJ32JmckoUn+40upr6hQcQ4GfmyQUTuHXxJobHh5IQGcjMoTFMHhTp/G+4bNtxVu7KIzEykEeuGsuqXXlcPSmRxWuPcH16MtuOldA/PNDZvx4ZbGVkv1C+3mcEfmOoN/6ucIyBBPtbGBLbckZN/7AAtnHq5imNLfzoEGuTPvy2xnn8LV0b/2mLBL4Q7mANguvfhlcvg7d/BDe+CynndOsprk9P5oU1h3j44z0UlNcSHezPizelOUMk0GpmTjuWcWhLeJAf/7h2XJNtE5Ij+cWcYfzm/Z18m3mCzIIKhsYZ3TLDHa3ayjo7E5IiGBYfgp9ZsTOnjEvHDuBAfgUBfiaSIoMYGh/KgYIKqups/PI/20iMDOTDn8zkziWb2ZFTyqPXjGsSjK6mD4lmtWPd+bhQY7bM0LgQVv/qPOc+r92W7vz31MHRKAUbDxdx5YQExidFMD4pAoB/3zqFNMcHg6nZt4mpg6NZvPYIYQEWBkUbrzHWcb4qx7o6wVYzQ1qZQnnuiFg+3ZXnnI4Z0hj4zQZj3RXsbZEuHSHcJSAMblwKUSnw5gLIzujew/uZ+ensYezOLaPWZmfxLZNb7V7obldPSiQ+zJ+nvszkQEEFw+KMoB8aF+K8IGlCciT+FjMj+oU6Z+rsL6hgaFwIJpNieFwImQUV/HHZLo4WVfHo1eOICrby2m3pfPPfs0hx9O23ZsaQU1csx4ed+fWOT4rgg7tnsGBKMjc1+9ZwzvBYgqwWgqyWFt0rMx1XRv/u0lQuGdufIKuZO84xunrSBxsXRk0cGElSZCDNzZ+cxOr/Ps/5LSrEpYXvqnEsYEB4gHPbhOSIM76mzpIWvhDuFBwNP/oAXrkYXrsSbv4I+o878/Pa6Zq0RPbnl3PJ2P7OfmN387eYuf3swTz88R7AGLAF4wNoYFQQWUVVjE0yunNGDwjn0115aK05kF/O1MHRzufU2hp4JyObu2cNcW63WkzEhQW0ctZTUgeEERZgoazG5pwPfybjkiIYl9SxID3/rDg+ufds5/jA7gfnOv926dgBzBwaQ0RQy8HVsYnhKKWaXHDXuHJmdLPZN433s52QHMlTZ6eQ2j/stNM8u0pa+EK4W1h/uHkZWENhyQ+hYE+3HdrPbOKPl49i8iD3XIrfluvTk51TB4e6dGmMT4pg9IBwwgKMQcfRCeGUVNXz6/d3kFta4/xwaPxwGpcYzi/mDO/Quc0mxbQh0YQGNL2qtrsppZxh35rWwn7Pg3N5987pLbY3hvioAU2X2Eh1HP+atEQmJke6NewBVPN5tb4kLS1NZ2R079dgIbzm5EF4ZR6g4ZZPIHqItyvqkn99fZB/frafdQ/MJtoxT72y1ka9vcEZhlknK7nxpQ2UVtWjlGLxLZOZkByJzd7AY5/t5/r0ZBIjO770xOETlRw5Wcl57Vx91N2OFVVxsrLOOTbQXEOD5tNdeVw0ql+LsYnGaavdSSm1WWud1mK7BL4QHlSwFxbPA0sg3LICIgd6u6JO01qTW1rTZJ0Y4RvaCnzp0hHCk+JGGn36deXw78uhrOO3SvQVSikJ+x7GY4GvlDpLKfWsUupdpdRdnjqvED6n/1i48X2oPAlvXAf1NWd+jove98AAABCbSURBVBDdoF2Br5R6WSlVoJTa2Wz7XKXUPqVUplLq/tMdQ2u9R2t9J3AtMKPzJQvRCyROgqtegLztsPLX3q5G9BHtbeEvBua6blBKmYGngYuBVGCBUipVKTVGKbW82U+c4zmXAx8DK7rtFQjRU424GKb/FDJe6vISDEK0R7vmNGmtVyulBjXbPAXI1FofAlBKvQVcobX+K3BpG8dZBixTSn0MvNHZooXoNc7/AxzdAMvuhf7je/zMHeHbutKHnwC43iEh27GtVUqpWUqpJ5RSz3GaFr5SapFSKkMplVFYWNiF8oToAcx+cM0rxqJr79wM9V27/6oQp+OxQVut9dda659pre/QWj99mv2e11qnaa3TYmNjPVWeEN4Tngg/fB7yd8CnD3i7GtGLdSXwc4Akl8eJjm1CiI4afiHM+DlsfgV2vOvtakQv1ZXA3wQMU0qlKKWswHxgWfeUJUQfNPu3kDgFPv6vHj0/X/iu9k7LfBNYB4xQSmUrpRZqrW3APcBKYA/wjta6e++sLERfYvaDHz4LtjpY9jPw4avgRc/U3lk6C9rYvgKZYilE94keAnP+BJ/8CrYugYk3ebsi0YvI0gpC+JrJt8Ogs+HTX0PJUW9XI3oRCXwhfI3JBFc8DWj48CfQ0ODtikQvIYEvhC+KHAgX/RkOrzauxBWiG0jgC+GrJt4MQy+Az34PJzK9XY3oBSTwhfBVSsHlT4LZCu8tNGbvCNEFEvhC+LKwAUbo534PX//F29WIHk4CXwhfl3q5MT3z28eNPn0hOkkCX4ieYO4jxhz9pXdAVZG3qxE9lAS+ED2BNRiuehEqC2H5z+UqXNEpEvhC9BQDJhjr7ez+ELa+5u1qRA8kgS9ETzL9Z46rcO+HokPerkb0MBL4QvQkJpOxwJoyG/35dpu3KxI9iAS+ED1NeCJc+hhkb4RvH/N2NaIHkcAXoicaczWMuQa+fgSyN3u7GtFDSOAL0VPN+zuE9oelt0NdpberET2ABL4QPVVghNGfX3QIvnjI29WIHkACX4ieLOVsSLsVNj4Hx7/3djXCx0ngC9HTnf97CIqG5b+ABru3qxE+TAJfiJ4uMAIu+isc3wKbX/F2NcKHSeAL0RuMuRpSzoXPH4TyfG9XI3yUBL4QvYFScMljYKuGVb/xdjXCR0ngC9FbxAyFmb+AHf+BrHXerkb4IAl8IXqTGT+HoBhY/ai3KxE+SAJfiN7EGgTT74GDX8gVuKIFCXwhepvJt0FABKz5u7crET5GAl+I3sY/FKbeDftWQN4Ob1cjfIgEvhC9UfoisIbCmn94uxLhQyTwheiNAiON0N/1ARTu83Y1wkdI4AvRW029G/wCYbX05QuDBL4QvVVwDExZBDvegewMb1cjfIAEvhC92Tm/NNbM//g+WVhNSOAL0av5h8KFD0Pu97B1iberEV4mgS9Ebzf6Khg4Az7/E1QVebsa4UUS+EL0dkrBvEehphS+fNjb1Qgv8ljgK6VmKaXWKKWeVUrN8tR5hRBA/ChjADfjZVlyoQ9rV+ArpV5WShUopXY22z5XKbVPKZWplLr/DIfRQAUQAGR3rlwhRKed9wCEDYAP7oT6Gm9XI7ygvS38xcBc1w1KKTPwNHAxkAosUEqlKqXGKKWWN/uJA9ZorS8G/gf4U/e9BCFEuwSEw+VPwon98JV07fRFlvbspLVerZQa1GzzFCBTa30IQCn1FnCF1vqvwKWnOVwx4N/xUoUQXTb0fJh0C6x9CkZeCslTvV2R8KB2BX4bEoBjLo+zgfS2dlZKXQlcBEQAT51mv0XAIoDk5OQulCeEaNWFDxnLJ39wF8z+LRQdhrIcGDsfktv8v7DoBboS+B2itV4KLG3Hfs8DzwOkpaVpd9clRJ/jHwpXPAOvXgbv3mpsswTAtrfhpg8habJ36xNu05XAzwGSXB4nOrYJIXxdytlwx2pQJogcBHWV8PJF8PrVcMsKY1aP6HW6Mi1zEzBMKZWilLIC84Fl3VOWEMLt+o+FfqPBPwRC443WvV8gLPkh5O088/NFj9PeaZlvAuuAEUqpbKXUQq21DbgHWAnsAd7RWu9yX6lCCLeKHAg/+gB0Azx3DnzyP1Bd4u2qRDdSWvtuN3laWprOyJBV/oTwqKoi44rcza9AYBT88DkYdoG3qxIdoJTarLVOa75dllYQQjQVFAWXPgaLvobQfvDGNfDdE+DDjUPRPhL4QojW9R8HC1fBWZfDZ78zpnHabd6uSnSBBL4Qom3WYLhmMcx6ALa9CSvuk5Z+D+axefhCiB5KKZh1P9hq4Nt/QugAmPU/RmvfXmt8KIgeQQJfCNE+5/8ByvPg678Yrf3SY4CCs+8zfixWb1cozkACXwjRPkoZi68FRED5cRj1Qyg+At88AnuWwXm/hmEXgkWWyvJVEvhCiPYz+8HFjzTdNvY64565b98IgZGQlG5cwRsQAWm3ylINPkQCXwjRNSPmwtAL4NBXsO0tOLHP2J71HWx7A1LOgTkPwoAJ3q1TSOALIbqB2QLD5hg/jWrLYfNi+O7/4PnzYPz1xro91cVQdhyqTsKgmTBuvrFduJ1caSuEcK+aUvjmb7DhWWiwgV+wcUGXfyjkbjP2GXsdzP4NRMiS6N2hrSttJfCFEJ5RXw0mizEO0KjkGGx6AdY/a6zhM/pKGH8D9BtjXPErOkUCXwjhu0qOwdon4fs3oK7c2GYJNL4FRCRD1GDjFo11lcZ0ULMfBMdCXKrx7SCsv3fr9zES+EII31dbDlnrjIHfinyjO6j4iPFTUwZ+QcYHQIMNynONO3VZQ+H838Pk28AkiwdA24Evg7ZCCN/hHwrDLzR+2uPkQVjxS/jkv2HHf2DQDGOFz+ghYLZCwiTpGnIhgS+E6Lmih8CNS2H72/D5n4ybszfUn/q7fzhMvcu4FqCuCuqrjG8HWhu3dRx2gXHtQB8hgS+E6NmUMqZ2jptvPK4qMrqA6iqMKaHfPNL2c83+kDARpt0DZ13qkXK9SQJfCNG7BEWd6sZJOQcqT0LhHvAPMxZ6M1mMK4ErC2DHe5D5GbzzI5j0Y1BmY7ZQULRx20e/IAiOg5A44z6/JrNXX1pXSeALIXq34GgIntlye0SS0cdf9xv4z49h29uOBeAU1JQYwe8qMNII/+AYOOeXMGS2J6rvVhL4Qoi+zRoMN/yn6Ta7zbgSuL4SKgqhJAsOr4baMuNisSU/hAk/gtm/NS4i6yEk8IUQojmzxejSAeMagOR0GHut8bi+Gr76C6x7GnZ/CHevg/BE79XaATJpVQghOsIvEC58CO5cY1wItvF5b1fUbtLCF0KIzogfBWddBhtfNG4MY7Ya3wYm/dhn5/5L4AshRGfNedCY13/kO9B24+rfXe/DbZ/75I1gJPCFEKKzIgfC/NdPPd6zHN6+AVY/agzo+hgJfCGE6C5nXQqjrzIC314HYQkwcAbEDPeJe/5K4AshRHea/TtjAbjv/u/UtpjhxtRPL9/oRQJfCCG6U1QK3LcHGuzGEg9Z38HK38K7C43QD4w0loPwAgl8IYRwB5PZWNwteogxsLv8F/C3FIgeBtN/CuMWeLybR+bhCyGEu026xVjV88I/Gx8EH/0M3rgGcrZ4tAwJfCGEcDelYOj5MP0euHs9zH0EDn0NSxd5tAzp0hFCCE9Sylijv+gwbHnVWJvfQ3360sIXQghviEoBW42xSJuHSOALIYQ3hCUYv0uzPXZKCXwhhPCGcEfglx332Ckl8IUQwhsaW/hlOR47pccGbZVSZwM3OM6ZqrWe7qlzCyGEzwmOA5OfcXGWh7Srha+UelkpVaCU2tls+1yl1D6lVKZS6v7THUNrvUZrfSewHHi18yULIUQvYDIZ99zd8R+w1XrmlO3cbzEw13WDUsoMPA1cDKQCC5RSqUqpMUqp5c1+4lyeej3wRjfULoQQPdu0u6EiHw585pHTtatLR2u9Wik1qNnmKUCm1voQgFLqLeAKrfVfgUtbO45SKhko1VqXt3UupdQiYBFAcnJye8oTQoieaeBMMFkgJ8NYadPNujJomwAcc3mc7dh2OguBV063g9b6ea11mtY6LTY2tgvlCSGEj/MLgPjRkJ3hkdN5dJaO1voPWuu1njynEEL4tMQ0OLIG1j5p3CDdjboS+DlAksvjRMc2IYQQ7TX1bug/Hlb9Fr76s1tP1ZXA3wQMU0qlKKWswHxgWfeUJYQQfUT0ELjjGxhxCWx9za0zdto7LfNNYB0wQimVrZRaqLW2AfcAK4E9wDta611uq1QIIXqz1MuhuhhKjrrtFO2dpbOgje0rgBXdWpEQQvRFwY5JKpWFEDPMLaeQpRWEEMIXOAP/hNtOIYEvhBC+IDjG+F0lgS+EEL1bULTxW1r4QgjRy1n8wT9cAl8IIfqE4Gjp0hFCiD4hONaYpeMmEvhCCOErgqKhqshth5fAF0IIXxEQATWlbju8BL4QQviKgHAJfCGE6BMCwqG2DBrsbjm8BL4QQviKwAjjt5ta+RL4QgjhKwLCjd8S+EII0ctJ4AshRB8RIF06QgjRNzhb+CVuObwEvhBC+Arp0hFCiD6icZZOtbTwhRCid7OGgF8wlOe65fAS+EII4SuUgsiBUJzllsNL4AshhC+JSIYSCXwhhOj9IhwtfK27/dAS+EII4UsiB0JdOVQXd/uhJfCFEMKXRA2ByBS3rItv6fYjCiGE6LwRc40fN5AWvhBC9BES+EII0UdI4AshRB8hgS+EEH2EBL4QQvQREvhCCNFHSOALIUQfIYEvhBB9hNJuWK+huyilSoEDp9klHGjrTgExwIluL8r9TveafPlcXTlWR5/b3v3bs9+Z9ult7zF5f3Xf/r78/hqotY5tsVVr7bM/wPOd/TuQ4e363fGaffVcXTlWR5/b3v3bs19fe4/J+6v79u+J7y9f79L5qIt/74k8+Zq681xdOVZHn9ve/duzX197j8n7q/v273HvL5/u0ukKpVSG1jrN23WI3kveY8Kd3PH+8vUWflc87+0CRK8n7zHhTt3+/uq1LXwhhBBN9eYWvhBCCBcS+EII0UdI4AshRB8hgS+EEH1Enwl8pVSwUupVpdQLSqkbvF2P6F2UUoOVUi8ppd71di2id1JK/cCRX28rpS7szDF6dOArpV5WShUopXY22z5XKbVPKZWplLrfsflK4F2t9e3A5R4vVvQ4HXl/aa0Paa0XeqdS0VN18D32gSO/7gSu68z5enTgA4uBJnf7VUqZgaeBi4FUYIFSKhVIBI45drN7sEbRcy2m/e8vITpjMR1/j/3W8fcO69GBr7VeDRQ12zwFyHS0uOqAt4ArgGyM0Ice/rqFZ3Tw/SVEh3XkPaYM/wt8orXe0pnz9cbgS+BUSx6MoE8AlgJXKaX+Re9bH0V4TqvvL6VUtFLqWWCCUuoB75Qmeom2MuynwAXA1UqpOztzYEvXa+sZtNaVwC3erkP0Tlrrkxh9q0K4hdb6CeCJrhyjN7bwc4Akl8eJjm1CdAd5fwl3c9t7rDcG/iZgmFIqRSllBeYDy7xck+g95P0l3M1t77EeHfhKqTeBdcAIpVS2Umqh1toG3AOsBPYA72itd3mzTtEzyftLuJun32OyWqYQQvQRPbqFL4QQov0k8IUQoo+QwBdCiD5CAl8IIfoICXwhhOgjJPCFEKKPkMAXQog+QgJfCCH6iP8Hvx0Y7WJceAcAAAAASUVORK5CYII=\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "# GMRES"
      ],
      "metadata": {
        "id": "rUi73jjV63kf"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "from scipy.sparse import linalg\n",
        "\n",
        "residuals = []\n",
        "\n",
        "def f(res):\n",
        "  global residuals\n",
        "  residuals.append(res)\n",
        "\n",
        "solution, info = linalg.gmres(A, b, callback=f)\n",
        "\n",
        "print(solution)"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "44PZG1IBNBew",
        "outputId": "0c9744f7-fc91-4183-80e4-5ce037b95ab5"
      },
      "execution_count": 9,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "[6.2061547  7.37529814 3.80638523 4.41712863 4.76537425 0.94556993\n",
            " 2.84892075 4.77705512 0.07258629 0.03185439 4.34497236 2.0025436\n",
            " 1.78032811 3.82502829 4.1088909  2.13914714 1.99739365 1.18787736\n",
            " 0.15356758 1.2260391  3.19792866 1.522114   1.75730334 1.38212578\n",
            " 0.25951973 2.04132572 0.23701198 1.36317349 0.70856287 1.05702538\n",
            " 1.05337964 0.33249999 0.62775716 2.18985603 1.4359648  1.41745245\n",
            " 0.59344453 0.94941067 1.36112542 0.2613947  0.41750947 0.81029216\n",
            " 0.368578   1.83640797 0.88954292 0.8717927  0.13297308 0.19800423\n",
            " 0.3937471  1.0557808  0.89984746 0.12003587 1.26470125 1.64665169\n",
            " 1.36254533 1.58606573 1.50160286 0.57418062 0.44352974 1.50439107\n",
            " 0.91804037 0.31613971 0.84617016 0.18497274 1.29127099 0.82262941\n",
            " 0.86661596 0.01958306 0.10478585 0.2526019  0.1925097  0.59928046\n",
            " 0.20294674 1.3063265  0.20549465 0.65569566 0.65432446 0.45434763\n",
            " 1.21169446 0.38511625 0.57156817 0.94589106 0.35760295 0.1468299\n",
            " 0.68206223 0.29666252 0.97546308 0.86374882 0.41636435 0.61622003\n",
            " 0.25710019 0.46219372 0.16834214 0.57996767 0.00917777 0.05477064\n",
            " 0.35649646 0.79328252 0.65727736 0.32342666]\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "plt.plot(range(1, len(residuals) + 1), residuals)\n",
        "plt.xscale(\"log\")\n",
        "plt.yscale(\"log\")"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 269
        },
        "id": "QlwLLzk4PrBo",
        "outputId": "26f5fbcd-e264-422f-b13f-1eb85b84f963"
      },
      "execution_count": 10,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXhU5dnH8e89k5UQwpKEfQkkhEVQNIKyKIIoKqittYJa37ZWigta7WurtXaztnaxtS4VUXitS1Vci6Di0iKgaAla2ddASlgT9i1ke94/EjBGwEAmc87M/D7XNdflHGbOuSeX85sz9/OcZ8w5h4iIRL+A1wWIiEh4KPBFRGKEAl9EJEYo8EVEYoQCX0QkRijwRURiRJzXBRxLenq669Kli9dliIhElAULFpQ45zLqbvdl4JvZaGB0dnY2+fn5XpcjIhJRzKzwSNt92dJxzr3unBuXlpbmdSkiIlHDl4EvIiKh58vAN7PRZjZp165dXpciIhI1fBn4aumIiISeLwNfRERCz5eBr5aOiEjo+TLwG9rSWVS0iw9Wl7Bx5wGqqrT8s4gI+HQefkNNnlvAa//ZCEBSfIAurVLompFCl1YpZKVX/3dWelNaNInHzDyuVkQkPKIy8H9yUU++eXpH1pbsY13JPtaW7GP55j28vWQLFbXO+NOS4+mSnkLX9OoPgkO3LukpNE2Myj+NiMQw8+MvXtW60va6VatWhWy/FZVVFO04wNpt+1hbXP1BcOi2cdcBav8pMlMTD38bqP3NoGPLJiTGBUNWk4hIqJnZAudc3pe2+zHwD8nLy3PhWlqhtLySwm37WVuyl4KS6g+EdduqPwxK9pYdflzAoEOLJkf8ZtCueTLBgFpEIuKtowW++hY1kuKD5LZJJbdN6pf+bdeB8sOtodq3lwp3sPdgxeHHJQQDdG7VpPoDICOFrJpvBlkZKWQ0TdR4gYh4SoFfD2nJ8ZzcsTknd2z+he3OOYr3HmRdyZe/GcxaWUxZRdXhxzZNjDs8PpBV69tBl/QU0pLjw/2SRCQGKfAbwMzITE0iMzWJ/lktv/BvlVWOjTsPfOlbwWfrdzJj4UZqzxZtlZLweWvo0DeDmrGDpHiNF4hIaPiyh99Yg7Z+cbCikvXb97O25pvB2pJ9FNQMIm/dc/ALj23fPJku6TVtovSmh78ZdGiRTFzQl5dRiIjHNGgbIfYerDjieEFB8V52l34+XhAXMDq1bHLEbwatU5MIaPBYJGZp0DZCNE2M46T2aZzU/otXGTvn2LG/vHqsoNYMooLifXywpoTS8s/HC5LjgzVjBZ9/Mzg0btAiJSHcL0lEfEKBHyHMjJYpCbRMaclpnb84XlBV5di8u5R1JfuqB45rbss2Hflis0Ph36NtKoOzM+jZNlUziERigAI/CgQCRrvmybRrnszA7PQv/Fv5oYvN6nwzmFewjVc+3QAsJ71pIkNy0hmSk87gnHQyU5O8eSEi0qgU+FEuPhg43Ocf1uOL/7Z5VylzV5cwZ1Uxs1cW8+qnGwDo0Sa15gMgg/5ZLTVTSCRKaNBWgOq20LLNu5mzqvoDYP7aHZRVVpEQF2BAVsvDHwA92qj9I+J3mqUjx+VAWSUfr912+ANg5Za9AGr/iESAiAr8aJ+HH4lqt3/mriph277q9YV6tEnlrO4ZDMlJ5/Quav+I+EFEBf4hOsP3p6oqx9JNn7d/8tdVt38S4wL0V/tHxHMKfGk0R2v/ZKQmMiS7uvWj9o9I+OjCK2k0yQlBhuZmMjQ3E6hu/8xZVcycVSW8v7K4Zvqn2j8iXtMZvjSqr2r/nJWTweCcdLV/REJILR3xhf1lFXy8djtzVpYwd/WX2z9DuqczKFvtH5GGUEtHfKFJQhzn5GZyzhHaP7NqtX96tm12ePqn2j8ioaEzfPGNY7V/huSkc1HftpzbszWpSfrBGJFjUUtHIs6h9s/7K4qZuWQzm3aVkhAX4OzuGYzq25bhPVvTNFFfUkXq8jzwzawrcBeQ5pz7Rn2eo8CXQ6qqHJ+u38mMhZt4Y9EmNu+uDv9zcjO4qG87hvfIJEXhLwI0MPDNbAowCtjqnDup1vaRwF+AIPCEc+6+euzrJQW+NERVleOT/+5gek34b91zkMS4AMN6ZHJR37YM65FJkwSFv8Suhgb+WcBe4KlDgW9mQWAlMAIoAuYDY6kO/9/W2cV3nXNba56nwJeQqapy5BfuYMbCjbyxeDPFew6SFB9geI/WXNS3LefkZpKcoAFfiS0NbumYWRdgeq3APxP4hXPu/Jr7dwI45+qGfd39HDPwzWwcMA6gU6dOpxUWFtarPpHKKsf8dduZsXATby7eRMneMpLjgwzvmcmovm0Zmpup2T4SExpjWmZ7YH2t+0XAgGMU0Aq4F+hnZnce7YPBOTcJmATVZ/gNqE9iTDBgnNG1FWd0bcUvLu7Nx2u3MWPhJt5avJnpCzfRJCHIuT2rz/zP7p6h8JeYE7ZGp3NuGzC+Po+ttVpm4xYlUSsYMAZ2S2dgt3R+eXFv/r12O9MXVYf/tM82kpIQ5NxerbmoT1vOUvhLjAh7S+d4qIcvoVZRWcVHBduZsWgjby3ezI795TRNjGNETfgP6Z5OYpzCXyJbY/Tw46getB0ObKB60PZK59ySEBSr9fCl0ZVXVjFvTU3bZ8lmdh0oJzUxjhG9WzOqb1sGZ2eQEBfwukyR49bQWTrPAUOBdGAL8HPn3GQzuxB4gOqZOVOcc/eGsmid4Uu4lFdW8cHqEmYs3MTMJZvZXVpBs6Q4LurbjvFnd6VzqxSvSxSpN88vvDoRCnzxQllFdfi/vnAj0xduorLKcekp7bnxnG50zWjqdXkiXymiAl8tHfGLLbtLeez9Av7+70LKKqq4+OR23DQsm+zMVK9LEzmqiAr8Q3SGL35RvOcgT8wp4Kl5hZRWVHJhn7ZMGJZNjzbNvC5N5EsU+CIhsG3vQSbPXcvfPlzHvrJKRvZuw4Th2fRul+Z1aSKHRVTgq6UjfrdzfxlT5q7l/z5Yx56DFZzbszU3D8+mb4fmXpcmElmBf4jO8MXvdh0o528frmPy3LXsOlDOObkZTBiew6mdWnhdmsQwBb5II9pTWs5T8wp5Yk4BO/aXMyQnnZuH53B6l5ZelyYxKKICXy0diVT7DlbwzEeFTJpdwLZ9ZZzZtRU3D8/hzG6tvC5NYkhEBf4hOsOXSHWgrJJnPy7ksdkFFO85SP8uLbl5eA6DslthZl6XJ1FOgS/igdLySl6Yv55HZ61h8+5STu3UnJuH53B29wwFvzQaBb6Ihw5WVPJifhGPzlrDhp0HOLlDGhOG5TC8Z6aCX0IuogJfPXyJVmUVVbz6aREP/2s167cfYGC3Vtxz6Ul005INEkIRFfiH6AxfolV5ZRXPz1/P799azsHyKsYP7cYNQ7tpXX4JiaMFvtZ+FfFAfDDAt87ozHs/PJsL+rThwfdWMfKB2cxZVex1aRLFFPgiHspMTeIvY/rxzLUDMDO+Nfnf3Pzcp2zdU+p1aRKFFPgiPjA4J503bxnCLcNzeGvxZobf/z5Pf1RIVZV/W64SeXwZ+GY22swm7dq1y+tSRMImKT7IrSO689YPhtC3Qxp3v7aYrz/6IUs26n0goaFBWxEfcs7xj/9s5NczlrJ9XxnfGZTFrSO60zQxzuvSJAJo0FYkgpgZl/Zrz3u3DWVM/05MnruWEX96n7cWb8bPJ2nibwp8ER9LaxLPb77Wh5evH0hacjzjn1nA9/6WT9GO/V6XJhFIgS8SAU7r3ILpEwZz14U9mVewjRF/ms3E99dQXlnldWkSQRT4IhEiLhjgurO68s5tZzMkJ5373lzOqAfnkr9uu9elSYRQ4ItEmPbNk5l0TR6PX5PHntJyvjFxHne8vJBd+8u9Lk18ToEvEqFG9GrNO7edzffP6sqLC4q46KE5LN6gKZxydL4MfM3DF6mflMQ47rywJy9fP5DKKsdlj37Iq58WeV2W+JQvA98597pzblxaWprXpYhEhFM6Nuf1CYM5pWNzbn3hM34xbYkGdOVLfBn4InL80psm8uz3BvC9wVk8+eE6rnz8I63JI1+gwBeJInHBAD8d1YsHx/Zj8YbdjHpwLgsKd3hdlviEAl8kCl18cjteuWEgyQlBxkyax9MfFeoKXVHgi0Srnm2bMe3GwQzOTufu1xZz+0sLKS2v9Los8ZACXySKpTWJZ/L/nM4tw3N4aUERl0+cp2UZYpgCXyTKBQLGrSO688Q1eawr2cfoh+Yyd1WJ12WJBxT4IjHi3F6tmTZhMBmpiVwz5WMmvr9Gff0YE9bAN7NLzexxM3vBzM4L57FFBLLSU3j1hkFc0Kct9725nBv//gl7D1Z4XZaESb0D38ymmNlWM1tcZ/tIM1thZqvN7I5j7cM595pz7jpgPHDFiZUsIg2RkhjHw2P78ZMLe/DW4s187ZEPKCje63VZEgbHc4b/JDCy9gYzCwKPABcAvYCxZtbLzPqY2fQ6t8xaT/1pzfNExANmxrizuvHMtQPYtq+MSx7+gHeWbvG6LGlk9Q5859xsoO46rP2B1c65AudcGfA8cIlzbpFzblSd21ar9jvgTefcJ0c6jpmNM7N8M8svLi4+0dclIvUwMDud1ycMJisjhXFP5/PIv1arrx/FGtrDbw+sr3W/qGbb0UwAzgW+YWbjj/QA59wk51yecy4vIyOjgeWJyFdp3zyZqd8/k9F92/GHmSv4wQv/0Xz9KBXWX0R2zj0IPPhVjzOz0cDo7Ozsxi9KREiKD/KXMaeQ2yaVP8xcwbqSfUy6Jo/WzZK8Lk1CqKFn+BuAjrXud6jZ1iBaLVMk/MyMG8/J5rFvncaqrXu5+OG5LCza6XVZEkINDfz5QI6ZZZlZAjAGmNbQorQevoh3zu/dhpevH0hcIMDlE+cx7bONXpckIXI80zKfA+YBuWZWZGbXOucqgJuAmcAyYKpzbklDi9IZvoi3erZtxrSbBtG3Qxo3P/cp97+9gqoqDeZGOvPziHxeXp7Lz8/3ugyRmFVWUcXdry3mhfz1nN+7NX/65imkJIZ16E9OgJktcM7l1d3uy6UV1NIR8YeEuAD3XdaHn43qxTtLt3DZox9q8bUI5svAV0tHxD/MjO8OzuL/vtOfDTsPcMnDHzB/Xd1LciQS+DLwRcR/zu6ewWs3DqJZcjxXPv4RU+ev/+onia/4MvDV0hHxp24ZTXnthkEMyGrFj15eyD3Tl1KhH0uPGL4MfLV0RPwrrUk8T37ndL49sAuT567luqfyOVCmK3MjgS8DX0T8LS4Y4BcX9+bXl57ErJXFXDPlY3aXlntdlnwFXwa+WjoikeHqMzrz0Nh+/Gf9TsZO+oiSvQe9LkmOwZeBr5aOSOQY1bcdk67JY03xXr752Dw27jzgdUlyFL4MfBGJLOfkZvLUdwdQvPsgl0+cx9qSfV6XJEegwBeRkOif1ZLnxp3BgfJKLp/4IUs37va6JKnDl4GvHr5IZDqpfRpTv38m8cEAYybNY0GhLtDyE18Gvnr4IpErO7MpL44/k5YpCVz9xL+Zs0q/XOcXvgx8EYlsHVo04cXxA+ncqgnXPpnPW4s3eV2SoMAXkUaSkZrIC+PO5KT2zbjh2U94MV9LMXhNgS8ijSatSTzPfG8Ag7LTuf2lhUyZu9brkmKaLwNfg7Yi0aNJQhxP/E8eI3u34VfTl/LAuyvx8+9wRDNfBr4GbUWiS2JckIev7Mc3TuvAA++u4g8zVyj0PaCfrhGRsIgLBvj9ZX2JDwb466w1xAUD3Daiu9dlxRQFvoiETSBg3HvpSVRVOR58bxVBM245N8frsmKGAl9EwioQMH779T5UOsef311JMAA3DVPoh4MCX0TCLhAwfndZX6qqHH98eyXBQIDrh3bzuqyop8AXEU8EA8YfLj+ZSuf43VvLCQZg3FkK/cakwBcRzwQDxv2Xn0xlleM3bywnYMb3hnT1uqyo5cvAN7PRwOjs7GyvSxGRRhYXDPDAFadQ5Ry/nrGMuIDx7UFZXpcVlTQPX0Q8FxcM8Jcx/Ti/d2t+8fpSnp63zuuSopIvA19EYk98MMBDY0/l3J6tufsfS3j240KvS4o6CnwR8Y2EuACPXNWPYT0yuevVxbww/79elxRVFPgi4iuJcUEevfpUhuZmcMcri5ixUEsrh4oCX0R8JzEuyMSrT+PUTi343xc/Y9km/VxiKCjwRcSXkuKDPHrVqTRLjmPc0/ns3F/mdUkRT4EvIr6V2SyJR68+jS27DjLhuU+prNIKmw2hwBcRXzu1Uwt+dUlv5qwq4Q8zV3hdTkQLW+CbWU8zm2hmL5nZ9eE6rohEvjH9O3HVgE5MfH8N0xdu9LqciFWvwDezKWa21cwW19k+0sxWmNlqM7vjWPtwzi1zzo0HvgkMOvGSRSQW/Xx0b/I6t+D2FxdqEPcE1fcM/0lgZO0NZhYEHgEuAHoBY82sl5n1MbPpdW6ZNc+5GJgBvBGyVyAiMSEhLsBfr9YgbkPUK/Cdc7OB7XU29wdWO+cKnHNlwPPAJc65Rc65UXVuW2v2M805dwFw1dGOZWbjzCzfzPKLi4tP7FWJSFTKTE1iogZxT1hDevjtgfW17hfVbDsiMxtqZg+a2WMc4wzfOTfJOZfnnMvLyMhoQHkiEo36dWrBPZdWD+L+fuZyr8uJKGFbLdM5NwuYVZ/HarVMETmWK07vxMKiXTz2fgEntUtj9MntvC4pIjTkDH8D0LHW/Q412xpMq2WKyFc5NIj7o5cWsnSjBnHroyGBPx/IMbMsM0sAxgDTQlGUmY02s0m7du0Kxe5EJArVHsT9/jP57NinQdyvUt9pmc8B84BcMysys2udcxXATcBMYBkw1Tm3JBRF6QxfROqj9iDuD174D1UaxD2m+s7SGeuca+uci3fOdXDOTa7Z/oZzrrtzrptz7t7GLVVE5Mv6dWrBz0b34v2VxTz0z9Vel+NrvlxaQS0dETkeVw3oxNf7teeB91Yye6Wmcx+NLwNfLR0ROR5mxr1f60Nu61Ruef5TNuw84HVJvuTLwBcROV7JCUH+etWplFc6bnj2Ew5WVHpdku/4MvDV0hGRE9E1oyl/vLwvn63fyb0zlnldju/4MvDV0hGREzXypLZcNySLp+YV8tqnIbk0KGr4MvBFRBriRyN70L9LS+58ZRErt+zxuhzf8GXgq6UjIg0RHwzw8JX9SEmMY/zTC9hTWu51Sb7gy8BXS0dEGiqzWRIPX9mPwu37+fHLC3FOF2X5MvBFRELhjK6t+NH5ubyxaDNTPljndTmeU+CLSFQbd1ZXzuvVmt++sYz56+r+rEds8WXgq4cvIqFiZvzxmyfToUUyNz77CcV7Dnpdkmd8Gfjq4YtIKDVLiufRq09jd2k5E577hIrKKq9L8oQvA19EJNR6tm3GvZf24aOC7dz/zkqvy/GEAl9EYsZlp3XgygGdeHTWGt5estnrcsJOgS8iMeVno3rRt0MaP3zxM9aV7PO6nLDyZeBr0FZEGktSfJBHrjyVYMAY/8wCDpTFziJrvgx8DdqKSGPq2LIJD1xxCiu27OHufyyOmYuyfBn4IiKNbWhuJhOG5fDSgiJemL/e63LCQoEvIjHrluE5DMlJ52fTlrCoKPpbyAp8EYlZwYDxlzH9SE9J4PpnF7Bzf5nXJTUqBb6IxLSWKQn89erT2LK7lNumfhbV/XwFvojEvFM6NueuC3vyz+VbeTG/yOtyGo0CX0QEuObMLvTPask9M5ayZXep1+U0Cl8Gvubhi0i4BQLG7y7rS1lFFXe/Fp1TNX0Z+JqHLyJeyEpP4bYR3Xl76RZmLNrkdTkh58vAFxHxyrWDs+jbIY2f/2MJO/ZF16wdBb6ISC1xwQC/u6wvuw6U86vpS70uJ6QU+CIidfRs24wbzsnm1U838K/lW70uJ2QU+CIiR3DjOd3IyWzKT15dxJ7Scq/LCQkFvojIESTGBfn9N/qyZXcp97253OtyQkKBLyJyFP06teC7g7J49uP/Mm/NNq/LaTAFvojIMfzwvFw6tWzCna8sjPi188Ma+GaWYmb5ZjYqnMcVETlRyQlB7rusD+u27eeBdyP7t3DrFfhmNsXMtprZ4jrbR5rZCjNbbWZ31GNXPwamnkihIiJeGdgtnSvyOjJ57lpWb93rdTknrL5n+E8CI2tvMLMg8AhwAdALGGtmvcysj5lNr3PLNLMRwFIgeuY4iUjMuH1kLsnxQe6dEblz8+Pq8yDn3Gwz61Jnc39gtXOuAMDMngcucc79FvhSy8bMhgIpVH84HDCzN5xzVUd43DhgHECnTp3q/UJERBpTetNEJgzP5jdvLGfWiq0Mzc30uqTj1pAefnug9u+CFdVsOyLn3F3OuR8AfwceP1LY1zxuknMuzzmXl5GR0YDyRERC69sDs+jSqgn3TF9KeeURI8zXwj5Lxzn3pHNu+rEeo9UyRcSPEuIC/PSiXqwp3sfT8wq9Lue4NSTwNwAda93vULOtwbRapoj41fCemQzJSeeBd1eyPcIWV2tI4M8Hcswsy8wSgDHAtFAUpTN8EfErM+PuUb3YV1bJn9+JrGma9Z2W+RwwD8g1syIzu9Y5VwHcBMwElgFTnXNLQlGUzvBFxM+6t07lqgGdePbjQlZs3uN1OfVWr8B3zo11zrV1zsU75zo45ybXbH/DOdfdOdfNOXdv45YqIuIft57bndSkeO6ZvjRifh3Ll0srqKUjIn7XIiWBW8/NYe7qEt5dFhmXF/ky8NXSEZFIcNUZncnObMqvZyzlYIX/19nxZeCLiESC+GCAu0f1onDbfp78YJ3X5XwlXwa+WjoiEinO7p7B8B6ZPPTP1RTvOeh1Ocfky8BXS0dEIsldF/XkYEUlf5y5wutSjsmXgS8iEkm6ZjTlO4OymLpgPYuK/NuZ8GXgq6UjIpHmpmHZtEpJ4JevL/HtNE1fBr5aOiISaZolxXP7+bnkF+5g2mcbvS7niHwZ+CIikejy0zrSp30a9725nP1lFV6X8yUKfBGREAkEjJ+P7sWmXaVMnLXG63K+xJeBrx6+iESqvC4tufjkdjw2u4D12/d7Xc4X+DLw1cMXkUh254U9CJhx35vLvS7lC3wZ+CIikaxtWjLXD+3GjEWb+Khgm9flHKbAFxFpBOPO6kr75sn88vWlVFb5Y5qmAl9EpBEkxQf5yYU9WbZpN8/P/6/X5QA+DXwN2opINLiwTxsGZLXk/rdXsutAudfl+DPwNWgrItHAzPjZ6F7s2F/Gg++t8rocfwa+iEi06N0ujSvyOvK3D9dRULzX01oU+CIijeyH5+WSFB/k3hnLPK1DgS8i0sgyUhOZMCyb95ZvZfbKYs/qUOCLiITBtwd1oXOrJtwzfSkVlVWe1KDAFxEJg8S46mmaq7bu5dmPvZmmqcAXEQmT83q1ZmC3Vvz53ZXs3F8W9uP7MvA1D19EotGhaZq7D5TzwLvhn6bpy8DXPHwRiVY92jRjbP9OPP1RIau37gnrsX0Z+CIi0ey2Ed1pkhDknunhnaapwBcRCbNWTRO5ZXgO768sDus0TQW+iIgHvnVmZ9o3T+aPb68I24+eK/BFRDyQGBfkB+fmsLBoFzOXbAnLMRX4IiIe+Vq/9nTNSOH+t1eEZc18Bb6IiEfiggF+OCKXVVv3Mu2zDY1+PAW+iIiHLjipDb3bNePP76yirKJxl1wIW+Cb2VAzm2NmE81saLiOKyLiZ4GA8b/n5fLf7fuZmr++cY9VnweZ2RQz22pmi+tsH2lmK8xstZnd8RW7ccBeIAkoOrFyRUSiz9DcDPI6t+Chf66itLyy0Y5T3zP8J4GRtTeYWRB4BLgA6AWMNbNeZtbHzKbXuWUCc5xzFwA/Bn4ZupcgIhLZzIzbz89ly+6DPD2vsNGOU6/Ad87NBrbX2dwfWO2cK3DOlQHPA5c45xY550bVuW11zh1qTu0AEo92LDMbZ2b5ZpZfXOzdutEiIuE0oGsrhuSk89dZq9lT2ji/f9uQHn57oHbDqahm2xGZ2dfN7DHgaeDhoz3OOTfJOZfnnMvLyMhoQHkiIpHl9vNz2bG/nClz1zXK/sM2aOuce8U5933n3BXOuVnHeqxWyxSRWNS3Q3NG9m7D43MK2LEv9MsnNyTwNwAda93vULOtwbRapojEqtvO605yQpBVW0P/g+dxDXjufCDHzLKoDvoxwJUhqUpEJEZ1b53KBz8eRkJc6Bsw9Z2W+RwwD8g1syIzu9Y5VwHcBMwElgFTnXNLQlGUWjoiEssaI+wBLFyrtJ2IvLw8l5+f73UZIiIRxcwWOOfy6m735dIKOsMXEQk9Xwa+Bm1FRELPl4EvIiKh58vAV0tHRCT0fBn4aumIiISeLwNfRERCz9fTMs2sGKi9dFwaUN8+TzpQEvKiosvx/D295lWtjX3cUO6/ofs60eefyPP0Xg6tun/Pzs65Ly9G5pyLmBsw6Tgem+91vX6/Hc/f0+ubV7U29nFDuf+G7utEn38iz9N7ObS3+v49I62l87rXBUSZSPp7elVrYx83lPtv6L5O9Pkn8rxI+n8vEtTr7+nrlk5DmFm+O8KVZiISWfReDp1IO8M/HpO8LkBEQkLv5RCJ2jN8ERH5omg+wxcRkVoU+CIiMUKBLyISI2Im8M0sxcz+ZmaPm9lVXtcjIsfPzLqa2WQze8nrWiJRRAe+mU0xs61mtrjO9pFmtsLMVpvZHTWbvw685Jy7Drg47MWKyBEdz/vYOVfgnLvWm0ojX0QHPvAkMLL2BjMLAo8AFwC9gLFm1ovqH1lfX/OwyjDWKCLH9iT1fx9LA0R04DvnZgPb62zuD6yuORMoA54HLgGKqA59iPDXLRJNjvN9LA0QjcHXns/P5KE66NsDrwCXmdmj6LJuEb874vvYzFqZ2USgn5nd6U1pkSvO6wLCxTm3D/iO13WIyIlzzm0DxntdR6SKxjP8DUDHWvc71GwTkcih93EjiMbAnw/kmFmWmSUAY4BpHtckIsdH7+NGENGBb2bPAfOAXDMrMrNrnXMVwE3ATGAZMNU5t8TLOkXk6PQ+Dh8tniYiEiMi+gxfRETqT4EvIhIjFPgiIjFCgS8iEiMU+CIiMUKBLyISIxT4IiIxQoEvIhIjFPgiIqI8zFwAAAAISURBVDHi/wG2bFFARhw1oQAAAABJRU5ErkJggg==\n"
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "code",
      "source": [],
      "metadata": {
        "id": "lNltyFAcQIgv"
      },
      "execution_count": null,
      "outputs": []
    }
  ]
}