{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical Methods\n",
    "\n",
    "## Lecture 7: Numerical Linear Algebra III\n",
    "\n",
    "### Exercise solutions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import scipy.linalg as sl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:blue\">Exercise 7.1: matrix norms</span>\n",
    "\n",
    "Write some code to explicitly compute the two matrix norms defined mathematically above and compare against the values found above using in-built scipy functions.\n",
    "\n",
    "Based on the above code and comments, what is the mathematical definition of the 1-norm and the 2-norm?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "\n",
    "def frob(A):\n",
    "    m, n = A.shape\n",
    "    squsum = 0.\n",
    "    for i in range(m):\n",
    "        for j in range(n):\n",
    "            squsum += A[i,j]**2\n",
    "    return np.sqrt(squsum)\n",
    "\n",
    "\n",
    "def mars(A):\n",
    "    m, n = A.shape\n",
    "    maxarsum = 0.\n",
    "    for i in range(m):\n",
    "        arsum = np.sum(np.abs(A[i]))\n",
    "        maxarsum = arsum if arsum > maxarsum else maxarsum\n",
    "    return maxarsum\n",
    "\n",
    "\n",
    "A = np.array([[10., 2., 1.],\n",
    "                 [6., 5., 4.],\n",
    "                 [1., 4., 7.]])\n",
    "\n",
    "\n",
    "print(frob(A) == sl.norm(A,'fro') and mars(A) == sl.norm(A,np.inf))\n",
    "\n",
    "\n",
    "print(np.allclose(frob(A), sl.norm(A,'fro')))\n",
    "print(np.allclose(mars(A), sl.norm(A,np.inf)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:blue\">Exercise 7.2: Ill-conditioned matrix</span>\n",
    "\n",
    "Consider a range of small values for $\\epsilon$ and calculate the matrix determinant and condition number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0 singular\n",
      "0.0019999999999997797 [ 1501.5 -3000. ]\n",
      "0.0039999999999995595 [  751.5 -1500. ]\n",
      "0.005999999999999339 [  501.5 -1000. ]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[2.,1.],\n",
    "                 [2.,1.]])\n",
    "b = np.array([3.,0.])\n",
    "print(sl.det(A), 'singular')\n",
    "\n",
    "for i in range(3):\n",
    "    A[1,1] += 0.001\n",
    "    print(sl.det(A), sl.inv(A) @ b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <span style=\"color:blue\">Exercise 7.3: Implement Gauss-Seidel's method.</span>\n",
    "\n",
    "Generalise the Jacobi code to solve the matrix problem using Gauss-Seidel's method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "499\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def gauss_seidel(A, b, maxit=500):\n",
    "    m, n = A.shape\n",
    "    x = np.zeros_like(b)\n",
    "    for k in range(maxit):\n",
    "        for i in range(m):\n",
    "            x[i] = 1/A[i,i] * (b[i] - np.dot(A[i,:i], x[:i]) - np.dot(A[i,i+1:], x[i+1:]))\n",
    "    print(k)\n",
    "    return x\n",
    "\n",
    "\n",
    "A = np.array([[10., 2., 3., 5.],\n",
    "                 [1., 14., 6., 2.],\n",
    "                 [-1., 4., 16.,-4],\n",
    "                 [5. ,4. ,3. ,11.]])\n",
    "b = np.array([1., 2., 3., 4.])\n",
    "\n",
    "# check gauss_seidel solution agrees with multiplying through by the inverse matrix\n",
    "np.allclose(sl.inv(A)@b, gauss_seidel(A, b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "def jacobi(A, b, maxit=500, tol=1.e-6):\n",
    "    m, n = A.shape\n",
    "    x = np.zeros(A.shape[0])\n",
    "    residuals = []\n",
    "    for k in range(maxit):\n",
    "        x_new = np.zeros(A.shape[0])\n",
    "        for i in range(m):\n",
    "            x_new[i] = (1./A[i, i]) * (b[i] - np.dot(A[i, :i], x[:i]) \n",
    "                                   - np.dot(A[i, i + 1:], x[i + 1:]))\n",
    "        x = x_new # update old solution \n",
    "        residual = sl.norm(A@x - b)\n",
    "        residuals.append(residual)\n",
    "        if (residual < tol): break       \n",
    "    return x, residuals\n",
    "\n",
    "def gauss_seidel(A, b, maxit=500, tol=1.e-6):\n",
    "    m, n = A.shape\n",
    "    x = np.zeros(A.shape[0])\n",
    "    residuals = []\n",
    "    for k in range(maxit):\n",
    "        for i in range(m):\n",
    "            x[i] = (1./A[i, i]) * (b[i] - np.dot(A[i,:i], x[:i]) \n",
    "                                   - np.dot(A[i,i+1:], x[i+1:])) \n",
    "        residual = sl.norm(A@x - b)\n",
    "        residuals.append(residual)\n",
    "        if (residual < tol): break\n",
    "        \n",
    "    return x, residuals\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[-0.16340807 -0.01532701  0.27335259  0.36893548]\n",
      "[-0.16340812 -0.01532701  0.27335261  0.36893553]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEWCAYAAAB42tAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeZxN9f/A8dd7dmvWSvYtjMyYMZaQfd+G7CVLtqwhRQlZslSyC2XJEFHWLGOJH0kxZCeGqNG3bMlSQvP5/XEuTdMYM+bOnDt33s/H4zzce+45n/ue23Tf89nFGINSSil1Px52B6CUUsq1aaJQSikVL00USiml4qWJQimlVLw0USillIqXJgqllFLx0kShlEJECoiIEREvu2NRrkcThXJ5IvKciESIyHUR+Z+IrBeRynbHlVaJyFsistDuOFTK0UShXJqIDAAmAWOAx4B8wAwg1M64YtK/wpW700ShXJaIPAKMBHoZY5YbY24YY24bY9YYY151XOMrIpNE5GfHMUlEfB2vVRORKBF5RUTOO2ojnRyvVRCRX0TEM8b7NRORg47HHiIyWEROicglEVkqItkcr91tpuksIj8CXzrOtxeRs47rh4rIGRGplYjyOojIjyJyUUSGxIjLU0TecNx7TUT2ikhex2vFRWSTiFwWke9FpFU8n+c2ERkrIrtF5HcRWXU3hjiufUJEVjvKjRSRro7z9YA3gNaOGt6Bh/qPq1IVTRTKlT0N+AEr4rlmCFABKA0EAuWAN2O8/jjwCJAb6AxMF5GsxphvgBtAjRjXPgd84njcF2gKVAWeAH4Dpsd676pACaCuiPhj1XSeB3LFeM+7ElJeZaAYUBMYJiIlHOcHAG2BBkBm4EXgDxHJAGxyxPyo45oZIlLyvp8WtHfc/wRwB5hyn+sWA1GO61oAY0SkpjFmA1bt7lNjTEZjTGA876XchTFGDz1c8sD60v3lAdecAhrEeF4XOON4XA34E/CK8fp5oILj8WhgruNxJqzEkd/x/BhQM8Z9uYDbgBdQADBAoRivDwMWx3ieHrgF1EpEeXlivL4baON4/D0QGsfP3hrYEevcLGD4fT6rbcC4GM/9HTF6xojBC8gL/A1kinHtWGC+4/FbwEK7fz/0SLlD21aVK7sE5BARL2PMnftc8wRwNsbzs45z98qIde8fQEbH40+Ar0WkB/AssM8Yc7es/MAKEYmOce/fWP0kd/0UK457z40xf4jIpRivJ6S8X+4TZ16shBhbfqC8iFyJcc4LCIvj2rhiPgt4AzliXfMEcNkYcy3WtSHxlKvcmDY9KVe2C7iJ1WRzPz9jfWHelc9x7oGMMUexvgDr8+9mJ7C+UOsbY7LEOPyMMediFhHj8f+APHefiEg6IHsiy7ufn4DC9zn/f7HKzGiM6RFPWXljPM6HVau5GOuan4FsIpIp1rV3Y9Ulp9MYTRTKZRljfsdq0pkuIk1FJL2IeItIfRF5x3HZYuBNEckpIjkc1ydm6OYnWP0HVYBlMc7PBN4WkfwAjvLjG2n1GdBYRCqKiA8wApAklBfTR8AoESkqlgARyQ58ATwpIi84PhdvESkbo28jLu1ExF9E0mMNFPjMGPN3zAuMMT8BXwNjRcRPRAKw+ncWOS75FSggIvr9kUbof2jl0owx72N15r4JXMD6K7o3sNJxyWggAjgIHAL2Oc4l1GKsvowvjTEx/7KeDKwGNorINeAboHw8cR4B+gBLsGoX17D6Q/56mPJieR9YCmwErgJzgHSOpqE6QBusWsAvwHjAN56ywoD5jmv9sJJkXNpi9Vv8jDWYYLgxZpPjtbsJ9ZKI7Evgz6BSMTFGa5FKOZuIZASuAEWNMT/YHQ9Yw2OxOqE/sjsWlbpojUIpJxGRxo7msQzAe1g1nDP2RqVU0mmiUMp5QrGaan4GimINb9Uqu0r1tOlJKaVUvLRGoZRSKl5uOeEuR44cpkCBAnaHoZRSqcrevXsvGmNyxj7vVolCRBoDjYsUKUJERITd4SilVKoiImfjOu9WTU/GWlW02yOPPGJ3KEop5TbcKlE4hifO/v333+0ORSml3IZbJQqtUSillPO5bR+FUsp13L59m6ioKG7evGl3KArw8/MjT548eHt7J+h6t5xHERISYrQzWynX8cMPP5ApUyayZ8+OiDz4BpVsjDFcunSJa9euUbBgwX+9JiJ7jTH/WU7erZqetI9CKdd08+ZNTRIuQkTInj17omp3bpUotI9CKdelScJ1JPa/hVsliqQKDw/ngw8+IDo6+sEXK6VUGuFWiSKpTU+ffvopPXv2pGrVqhw/ftzJ0Sml7JQxY8YHX5RAHTt25LPPPvvP+YiICPr2vd8WH6mXWyWKpDY9zZkzh3nz5nHkyBECAwN5++23uX37tpOjVEq5q5CQEKZMmWJ3GE7nVokiqUSEjh07cvToUUJDQ3nzzTcJCQnR5UCUchPXr1+nZs2aBAcHU6pUKVatWnXvtQULFhAQEEBgYCAvvPACAGfPnqVmzZoEBARQs2ZNfvzxx3vXb968mWeeeYYnn3ySL774AoBt27bRqFGjlP2hUoDOo4jD448/ztKlS1m1ahU9e/akfPny9O/fn5EjR5I+fXrnBKtUGtWvXz/279/v1DJLly7NpEmTHnidn58fK1asIHPmzFy8eJEKFSrQpEkTjh49yttvv83OnTvJkSMHly9fBqB37960b9+eDh06MHfuXPr27cvKldYuvGfOnOH//u//OHXqFNWrVycyMtKpP5MrcasaRVKbns6dg6+++ud5aGgoR44coUuXLkyYMIFSpUrx5ZdfOilapVRKM8bwxhtvEBAQQK1atTh37hy//vorX375JS1atCBHjhwAZMuWDYBdu3bx3HPPAfDCCy/wVYwviFatWuHh4UHRokUpVKiQW/drulWNIql694bwcPjiC6hRwzqXJUsWZs2aRdu2benatSs1a9akc+fOvPvuu2TNmtXegJVKhRLyl39yWbRoERcuXGDv3r14e3tToEABbt68iTEmQUNGY14T+3p3Hv7rVjWKpJo1CwoVgkaNYMuWf79WrVo1Dh48yKBBg5g/fz7+/v58/vnn9gSqlHoov//+O48++ije3t5s3bqVs2etVbVr1qzJ0qVLuXTpEsC9pqeKFSuyZMkSwEoylStXvlfWsmXLiI6O5tSpU5w+fZpixYql8E+TcjRRxPDoo7B1KxQpYiWLzZv//Xq6dOkYN24cu3fvJleuXLRo0YJnn32W//3vf/YErJRKkDt37uDr68vzzz9PREQEISEhLFq0iOLFiwNQsmRJhgwZQtWqVQkMDGTAgAEATJkyhXnz5hEQEEBYWBiTJ0++V2axYsWoWrUq9evXZ+bMmfj5+dnys6UIY4zbHEBjYHaRIkVMUpw/b0xAgDF+fsaEh8d9za1bt8y4ceOMn5+fyZw5s3n77bfN9evXk/S+Srmro0eP2vr++/fvN2XLlrU1BlcT138TIMLE8d3qVjUK46QlPHLmtJqeihWDJk2sfovYvL29GTRoEAcOHKBq1aoMGTKEwoULM23aNG7dupWk91dKOc/MmTNp27Yto0ePtjuUVMutEoUz5chhJYsSJSA0FNavj/u6J598ktWrV7Nz506KFy9Onz59KFasGAsWLODvv/9O2aCVUv/x0ksvcfToUerUqWN3KKmWJop4ZM9uJQt/f2jaFNauvf+1FStWZOvWrWzYsIFs2bLRoUMHAgMDWbly5d1mMaWUSpU0UTxAtmxWp3apUtCsGaxZc/9rRYS6deuyZ88eli5dyu3bt2nWrBlPP/00W7duTbmglVLKiTRRJEC2bLBpE5QuDc2bw+rV8V/v4eFBy5YtOXLkCB9++CHnzp2jRo0a1KlTR5cDUUqlOi6fKEQkg4h8LCIfisjzdsWRNSts3AhBQdCiBThm8cfLy8uLLl26cPLkSSZMmMC+ffsoW7YsLVq0cOtZnEop92JLohCRuSJyXkQOxzpfT0S+F5FIERnsOP0s8JkxpivQJMWDjSFLFitZBAdDy5awYkXC7vPz82PAgAGcPn2a4cOHEx4eTsmSJfnggw+SN2Cl1D2//vorzz33HIUKFaJMmTI8/fTTrEjo/8RO8s0331C+fHlKly5NiRIleOutt+K9Pr5lywsUKMDFixfjvT8h1yRIXGNmk/sAqgDBwOEY5zyBU0AhwAc4APgDrwOlHdd8kpDyy5Qpk+QxxvH5/Xdjnn7aGC8vYz77LPH3nz9/3jRs2NAAZtq0ac4PUCkXY/c8iujoaFOhQgXzwQcf3Dt35swZM2XKlBSN48knnzT79+83xhhz584dc+TIkYcuK3/+/ObChQsPfY3Lz6MwxmwHLsc6XQ6INMacNsbcApYAoUAUkMdxzX3jFZFuIhIhIhEXLlxIjrDvyZwZNmyAcuWgdWtYtixx9+fMmZPPP/+cJk2a0Lt3b6ZNm5Y8gSqlAPjyyy/x8fHhpZdeuncuf/789OnThzNnzvDMM88QHBxMcHAwX3/9NfDfJcN79+7N/PnzARg8eDD+/v4EBAQwcOBAwFrS46mnniIwMJAqVarEGcf58+fJlSsXAJ6envj7+wNw48YNXnzxRcqWLUtQUNC95c9jxnDp0iXq1KlDUFAQ3bt3/9doyoULF1KuXDlKly5N9+7dnT4035UWBcwN/BTjeRRQHpgCTBORhsB9xxwZY2aLyP+Axj4+PmWSNVL+SRb160PbtnDnjvVvQvn6+rJs2TJatWpFnz59iI6OdsudsZSKrV8/cPIq45QuDfGtNXjkyBGCg4PjfO3RRx9l06ZN+Pn5cfLkSdq2bRvvoJPLly+zYsUKjh8/johw5coVAEaOHEl4eDi5c+e+dy62/v37U6xYMapVq0a9evXo0KEDfn5+vP3229SoUYO5c+dy5coVypUrR61atf5174gRI6hcuTLDhg1j7dq1zJ49G4Bjx47x6aefsnPnTry9venZsyeLFi2iffv28X1kieJKndlxLb1ojDE3jDGdjDE9jDGL4ivAOGlmdkJlymQli0qVoF07WLAgcff7+PiwdOlSmjZtyssvv2zrqppKpSW9evUiMDCQsmXLcvv2bbp27UqpUqVo2bIlR48ejffezJkz4+fnR5cuXVi+fPm9PWoqVapEx44d+fDDD+/7F/2wYcOIiIigTp06fPLJJ9SrVw+AjRs3Mm7cOEqXLk21atW4efPmvzZJAti+fTvt2rUDoGHDhvdWr96yZQt79+6lbNmylC5dmi1btnD69OkkfT6xuVKNIgrIG+N5HuDnxBTgrI2LEiNjRli3zpq93bEj3L4NnTsn/P67yaJNmzb0798fYwz9+/dPtniVspsdfw+VLFnyX6s9T58+nYsXLxISEsLEiRN57LHHOHDgANHR0fcW9/Py8iI6OvrePTdv3rx3fvfu3WzZsoUlS5Ywbdo0vvzyS2bOnMm3337L2rVrKV26NPv372fgwIF89913PPHEE6xbtw6AwoUL06NHD7p27UrOnDm5dOkSxhg+//zz/6xA++uvv/7reVxLmRtj6NChA2PHjnXOhxUHV6pR7AGKikhBEfEB2gAPmLHwbyldo7grQwZrIl7dutClC8yYkbj7vb29WbJkCc2bN2fAgAFMmDAheQJVKo2qUaMGN2/e/NdIwz/++AOwlh7PlSsXHh4ehIWF3asN5M+fn6NHj/LXX3/x+++/s8Wx98D169f5/fffadCgAZMmTbq3W9+pU6coX748I0eOJEeOHPz000/MmzeP/fv330sSa9euvde3cPLkSTw9PcmSJQt169Zl6tSp91777rvv/vMzVKlShUWLrEaV9evX89tvvwHWEumfffYZ58+fB6ymsbvLpzuLLTUKEVkMVANyiEgUMNwYM0dEegPhWCOg5hpjjiSy3BSvUdyVLp01t6JlS+jVC27dstpiE8rb25vFixfz/PPPM3DgQIwx9zrJlFJJIyKsXLmS/v37884775AzZ04yZMjA+PHjCQ4Opnnz5ixbtozq1auTIUMGAPLmzUurVq0ICAigaNGiBAUFAXDt2jVCQ0PvbXg0ceJEAF599VVOnjyJMYaaNWsSGBj4nzjCwsLo378/6dOnx8vLi0WLFuHp6cnQoUPp168fAQEBGGMoUKDAvX247xo+fDht27YlODiYqlWrki9fPgD8/f0ZPXo0derUITo6Gm9vb6ZPn07+/Pmd9wHGNRQqtR44aZnxpPjrL2OaNzcGjBk/PvH3375927Ru3doAZvzDFKCUC7J7eKz6r8QMj3WlPookM8asAdaEhIR0tSsGHx9YsgReeAEGDYK//oKhQxN+v5eXFwsXLkREGDRoENHR0QwePPjBNyqlVDJxq0ThKry8YOFC8PaGYcOsZqiRIyGhW+p6eXkRFhaGiPD6668THR3NG2+8kbxBK6XUfbhVorCzjyI2T0+YN8+qYYwebSWLceMSlywWLFiAh4cHQ4YMwRjDkCFDkjdopZKRMSbOUTsq5ZlEbn3gVonCFZqeYvL0hNmzrWTxzjtWM9TEiYlLFh9//DEeHh68+eabREdHMzQx7VhKuQg/Pz8uXbpE9uzZNVnYzBjDpUuXErXHt1slClfk4QHTp1vJYvJkq2YxbZp1PiE8PT2ZN28eIsKwYcM4f/4848ePvzfJR6nUIE+ePERFRZHcy+uohPHz8yNPnjwPvtDBrRKFKzU9xSRi1SR8fa2axe3bMGtW4pLF3LlzyZo1K5MnT2bDhg3MmzePypUrJ2/gSjmJt7c3BQsWtDsM9ZBcacJdkhmbJtwlhIjVRzF0KHz0kTWLOzHrdnl6ejJp0iS+/PJL/v77b6pUqUK/fv3uTRpSSqnk4laJwtWJWKOfRo6EsDB48UWIsUJAglSvXp2DBw/Ss2dPJk+eTEBAADt27EiegJVSCjdLFCLSWERm//7773aHEq+hQ61ksWAB9OwJiRyAQMaMGZk2bRpbt24lOjqaqlWr0q9fP27cuJE8ASul0jS3ShSu3PQU25tvwuuvW30V/fsnPlkAVKtWjYMHD9KrVy8mT55MYGCg1i6UUk7nVokiNRGBt9+21oOaPBneeOPhkkXGjBmZOnUq27ZtwxhD1apVefnll7V2oZRyGk0UNhKB99+Hl16yOrpHjXr4sqpWrcrBgwfp3bs3U6ZMITAwkO3btzsvWKVUmqWJwmYi1jyLjh1h+HB4992HLytDhgxMmTLlX7WLvn37au1CKZUkbpUoUktndmweHtaQ2TZt4LXXYOrUpJV3t3bRt29fpk6dytNPP82lS5ecE6xSKs1xq0SRmjqzY/P0tEZBNWsGfftaS38kRYYMGe5Nzjtx4gT169fn6tWrzglWKZWmuFWiSO28vWHxYqhf3+q3CAtLepl169Zl2bJlfPfddzRu3Fgn6CmlEk0ThYvx9YXPP4fq1a1+i6VLk15m48aNCQsLY8eOHbRo0YJbt24lvVClVJqhicIFpUsHq1dDxYrw/POwalXSy2zTpg2zZs1i/fr1tGvX7t6+wEop9SCaKFxUhgywdi0EB0OrVhAenvQyu3btyoQJE1i2bBldu3YlOrHrhyil0iSXTxQiUkhE5ojIZ3bHktIyZ4YNG8DfH5o2ha1bk17mgAEDGD58OPPmzaN///6J3sBEKZX2JGuiEJG5InJeRA7HOl9PRL4XkUgRiXdDaGPMaWNM5+SM05VlzQobN0KhQtC4MezcmfQyhw8fTv/+/ZkyZQrDhw9PeoFKKbeW3PtRzAemAQvunhART2A6UBuIAvaIyGrAExgb6/4XjTHnkzlGl5czJ2zZAlWrQoMGsHkzlC378OWJCBMmTODq1auMGjWKzJkzM3DgQOcFrJRyK8maKIwx20WkQKzT5YBIY8xpABFZAoQaY8YCjR72vUSkG9ANIF++fA9bjMt6/HErWVSpAnXqWM1QpUs/fHkiwqxZs7h27RqvvvoqmTNnplu3bs4LWCnlNuzoo8gN/BTjeZTjXJxEJLuIzASCROT1+11njJkNjAD2+fj4OCtWl5InD3z5JWTKBLVqweHDD74nPp6enoSFhdGwYUNeeuklPvnkE+cEqpRyK3Ykirh2Vr9vj6ox5pIx5iVjTGFHreO+UvPM7IQqUMBKFr6+VrL4/vuklefj48OyZcuoWrUq7du3Z5UzxuIqpdyKHYkiCsgb43ke4GdnFJxa13pKrCJFrGYoY6BGDTh1KmnlpUuXjtWrV1OmTBlatWrFli1bnBOoUsot2JEo9gBFRaSgiPgAbYDVNsSRqhUvbiWLv/6yksXZs0krL1OmTKxfv54nn3yS0NBQdu3a5ZxAlVKpXnIPj10M7AKKiUiUiHQ2xtwBegPhwDFgqTHmiDPeLy00PcX01FOwaRNcvWoli6iopJWXLVs2Nm3aRK5cuahduzYjRozQhQSVUog7TbgSkcZA4yJFinQ9efKk3eGkmN27rf6KXLng//7PGiGVFFFRUbz88sssX76c7NmzM3jwYHr16kW6dOmcE7BSyiWJyF5jTEjs8y4/Mzsx0lqN4q5y5WD9ejh3DmrWhAsXklZenjx5+Pzzz9m9ezdlypTh1VdfpUiRIsycOZPbt287J2ilVKrhVokiLatUCdasgdOnoXZtuHw56WWWLVuW8PBwtm7dSoECBejRowfFixdn4cKFuqigUmmIWyWKtDLq6X6qV7dWmj12DOrWBWd9DNWqVeOrr75i7dq1ZMqUiRdeeIHAwEBWrlypa0UplQa4VaJIq01PMdWpY+1nceCAtQHStWvOKVdEaNCgAfv27ePTTz/l9u3bNGvWjPLly7N582ZNGEq5MbdKFGm9RnFXo0awZInVyd2oEThzUzsPDw9atWrFkSNHmDNnDr/88gu1a9emZs2aOqRWKTflVolCaxT/ePZZWLgQvvoKQkPh5k3nlu/l5cWLL77IyZMnmTx5MkeOHKFixYpUqlSJRYsW8ddffzn3DZVStnGrRKH+rU0bmDvXWm22dWtIjgFLvr6+9O3bl1OnTvH+++9z4cIF2rVrR548eXj99dc5c+aM899UKZWiNFG4uQ4dYPp0a2vV9u0huQYrZcyYkf79+3P8+HE2btxI5cqVeeeddyhUqBCNGjVi3bp1OlJKqVTKrRKF9lHErWdPeOcdq9+iWzdIzh1QPTw8qF27NitWrODMmTO8+eabRERE0LBhQ4oWLco777zDxYsXky8ApZTTudXM7LtCQkJMRESE3WG4nOHDYeRI6NsXJk0CiWsd32Rw69YtVq5cyYwZM/i///s/fH19adWqFT179qR8+fJISgWilIpXmpiZreL31lswYABMmQJvvply7+vj40OrVq3Ytm0bhw8fpmvXrqxcuZKnn36aMmXKsHHjxpQLRimVaJoo0hAReO896N4dxoyxjpRWsmRJpk6dys8//8zMmTO5ceMG9erVY+jQodqHoZSL0kSRxojAjBnQrh0MGWLVLuyQMWNGunfvznfffUenTp0YPXo0tWvX5pdffrEnIKXUfblVotDO7ITx8IB586BZM3j5ZZgzx75Y0qdPz5w5c5g3bx7ffPMNpUuXZuvWrfYFpJT6D7dKFDrhLuG8vGDxYqhXD7p2tR7bqWPHjuzevZssWbJQq1YtRo8eTXRyDs9SSiWYWyUKlTi+vta6UFWqwAsvWAsK2umpp54iIiKCNm3aMHToUBo0aMCFpK6ZrpRKMk0UaVz69Nby5CEh0KqVtWOenTJmzMjChQuZOXMm27ZtIygoiJ07d9oblFJpnCYKRaZM1sZHJUpA06bW+lB2EhG6d+/Orl278PPzo2rVqrz33nu6Qq1SNtFEoQDImhU2boS8eaFhQ3CF+YpBQUHs3buXpk2b8uqrrxIaGsplZ+zIpJRKlFSRKESkqYh8KCKrRKSO3fG4q0cfhS1bIHt2a+OjQ4fsjggeeeQRli1bxuTJk9mwYQPBwcHs3r3b7rCUSlOSPVGIyFwROS8ih2Odryci34tIpIgMjq8MY8xKY0xXoCPQOhnDTfNy57aSRbp01paqJ07YHZHVFNW3b1927NiBMYbKlSszYsQIzp8/b3doSqUJKVGjmA/Ui3lCRDyB6UB9wB9oKyL+IlJKRL6IdTwa49Y3HfepZFSwoLU0eXQ01KwJrrJSePny5fnuu+9o1KgRb731Fnny5OG55567l0CUUskj2ROFMWY7ELthuRwQaYw5bYy5BSwBQo0xh4wxjWId58UyHlhvjNkX1/uISDcRiRCRCB1SmXTFi1vJ4sYNK1mcO2d3RJZs2bKxfPlyjh49So8ePVi3bh1VqlQhICCAGTNmcPXqVbtDVMrt2NVHkRv4KcbzKMe5++kD1AJaiMhLcV1gjJkNjAD2+fj4OCvONC0gADZsgAsXoFYtcKWWnhIlSjB58mTOnTvHRx99hI+PD7169SJ37tz06NGDgwcP2h2iUm7DrkQR17rS9207MMZMMcaUMca8ZIyZGc91OjPbycqVg7Vr4exZqFMHfvvN7oj+LUOGDHTu3JmIiAi+/fZbmjdvzvz58wkMDKRy5cq6LatSTmBXoogC8sZ4ngf4OamF6lpPyeOZZ2DlSjh2DOrXh2vX7I7ov0SEcuXKMX/+fM6dO8eECRP49ddf723LOnjwYN2WVamHZFei2AMUFZGCIuIDtAFW2xSLSoA6dWDpUmt+RaNG8Mcfdkd0f9myZWPAgAF8//33bNy4kWeeeYZ3330Xf39/1q9fb3d4SqU6KTE8djGwCygmIlEi0tkYcwfoDYQDx4ClxpgjSX0vbXpKXqGhsHAh7NhhrTzr6i06d7dlXb58OT/88APFixenSZMmLFmyxO7QlEpVvOJ7UUSuEXffgQDGGJP5QW9gjGl7n/PrgHUJCTKhRKQx0LhIkSLOLFbF0KaNVZvo3Blat4Zly8Db2+6oHixfvnxs3bqVJk2a8Nxzz3HlyhVeeinOcRFKqVjirVEYYzIZYzLHcWRKSJJIaVqjSBkvvghTp1qrzXboAKllY7pHHnmEDRs20LBhQ3r06MGYMWN0/oVSCRBvjSI2x+Q3v7vPjTE/Oj0ilSr07m3NsRg82FqBdvZsa0MkV5cuXTqWL19Op06dGDJkCJcvX+bdd99FJK6BeEopSGCiEJEmwATgCeA8kB+rbyiu9TcAACAASURBVKFk8oWWeNr0lLIGDbKSxahRVrKYPNnaatXVeXt7s2DBArJly8aECRP47bffmDVrFl5eifq7Sak0I6H/Z4wCKgCbjTFBIlIdiLPvwU7GmDXAmpCQkK52x5JWjBgB16/DxImQMSOMGWN3RAnj4eHB5MmTyZYtGyNGjOC3337jk08+wc/P78E3K5XGJDRR3DbGXBIRDxHxMMZsdSyp4VK0RpHyRGDCBKuDe+xYyJABhgyxO6qEERHeeustsmbNSr9+/WjUqBErVqwgU6ZMdoemlEtJaKvyFRHJCGwHFonIZOBO8oX1cLQz2x4iMGMGtGsHb75pNUGlJi+//DIff/wx27Zto1atWly6dMnukJRyKQlNFKHAn0B/YANwCmicXEGp1MfDA+bNs+ZX9OsHc+faHVHitG/fnuXLl3PgwAGqVKnCOVdZBVEpF5CgRGGMuWGM+dsYc8cY87Fj7SX9s0v9i5cXLF4M9epBly6Q2ua1NWnShA0bNvDTTz9RuXJlIiMj7Q5JKZeQoEQhItdE5KrjuCkif4uIy63nrGs92c/XFz7/3Fof6oUXYM0auyNKnGrVqrF161auX79O5cqVOXDggN0hKWW7hNYoYk688wOaA9OSN7TE0z4K15A+vZUggoKgZUtrX4vUpEyZMuzYsQNvb2+qVKnCmDFj+M3Vls1VKgU91BQpY8xKoIaTY1FuJHNmay+LJ5+01ojaudPuiBKnePHi7Ny5k0qVKjFkyBDy5cvHwIEDiYqKsjs0pVJcQpueno1xtBCRccSzf4RSANmywaZNkCcPNGgA++Lcm9B15cuXj3Xr1rF//35CQ0OZNGkShQoVolOnThw7dszu8JRKMQmtUTSOcdQFrmGNhFIqXo89ZjU9Zc1qLVV+JMlrBKe8wMBAFi5cSGRkJC+99BKffvop/v7+hIaG8vXXX9sdnlLJTtxpUbQYE+66njx50u5wVAyRkVClivV4xw4oXNjeeJLi4sWLTJs2jalTp3L58mUqV67MoEGDaNCgAR6pYcErpe5DRPYaY0L+cz6+RCEiU4l/i9K+zgnPuUJCQkxERITdYahYjh61kkXGjFayyJv3wfe4shs3bjBnzhwmTJjAjz/+SMmSJXnttddo27Yt3qlh7XWlYrlfonjQnz8RwF6sFWODgZOOozSQShaXVq7C3x82brT23a5VC3791e6IkiZDhgz07duXyMhIFi5ciIeHBx06dKBw4cKsWLHC7vCUcpoH7UfxsTHmY6AoUN0YM9UYMxWoiZUslEqU4GBYvx6ioqB2bbh82e6Iks7b25vnn3+eAwcOsG7dOnLmzMmzzz7LwIEDuX37tt3hKZVkCW1QfQKIuVJaRsc5pRKtYkVr06Pvv4f69eGqy03dfDgiQv369fn666/p1asXEyZMoHr16rociEr1EpooxgHfich8EZkP7ANSyYLSyhXVqgWffWYNmW3c2Fp91l34+voybdo0Fi9ezP79+wkKCmLLli12h6XUQ0vozOx5QHlgheN42tEklexEpISIzBSRz0SkR0q8p0oZjRvDwoVWx3bz5vDXX3ZH5Fxt2rRhz5495MyZk9q1azN69Giio6PtDkupRIs3UYhIcce/wVhNTT85jicc5+IlInNF5LyIHI51vp6IfC8ikSIyOL4yjDHHjDEvAa2A//TGq9StdWv46CNrFnfbtnDH5RavT5oSJUqwe/dunnvuOYYOHUqjRo10GXOV6jxo46IBQDesbVBjMzx4GY/5WGtCLbh7QkQ8gelAbSAK2CMiqwFPYGys+180xpx3bMU6GBdcX0ol3YsvWrvkvfwydOoEH3+cOvbfTqgMGTIQFhZG5cqVefnllwkKCmLZsmWUL1/e7tCUShhjTLIeQAHgcIznTwPhMZ6/DryewLLWxvNaN6zhvBH58uUzKvUZM8YYMKZ7d2Oio+2OJnns2bPH5M+f33h7e5upU6eaaHf9QVWqBESYOL5fE7rWU0sRyeR4/KaILBeRoIfMTbmxmq/uinKcu997VxORKSIyC1h3v+uMMbOBEcA+Hx+fhwxN2en1161j1iwYOBDcaNGAe0JCQti3bx9169alT58+tG3blmvXrtkdllLxSmgFf6gx5pqIVMZa6+ljYOZDvqfEcS6+2d/bjDF9jTHdjTHT4yvY6DLjqd7bb0OfPvD++zBihN3RJI9s2bKxatUqxo4dy7JlyyhbtixHUuMiWCrNSGiiuDsLuyHwgTFmFfCwf7ZHATEXb8gD/PyQZf2LblyU+onApElWv8WIEfDuu3ZHlDw8PDwYPHgwW7Zs4cqVK5QrV46ZM2fy99+64IFyPQlNFOccTT+tgHUi4puIe2PbAxQVkYIi4gO0AVY/ZFnKDXl4wOzZ1oio116DDz6wO6LkU61aNb777jvKly9Pjx49CAoKYtOmTXaHpdS/JPTLvhUQDtQzxlwBsgGvPugmEVkM7AKKiUiUiHQ2xtwBejvKOwYsNcY4pd6tTU/uw9MTwsKsuRY9e8KCBQ++J7XKlSsXW7ZsYenSpVy/fp06derQsGFDjh49andoSgGJWGbc0T9R1BgzT0RyAhmNMT8ka3SJpMuMu5+bN6FRI9i6FZYutSbmubO//vqLqVOnMmrUKG7cuEG3bt146623ePTRR+0OTaUBD7t67N2bhwODsIayAngDC50XnnNojcL9+PlZ60JVqGBNyFt333Fv7sHX15eBAwfe2yRp9uzZFC1alPHjx3Pz5k27w1NpVEKbnpoBTYAbAMaYn/n3IoFKJZsMGWDtWihVyqpRbNtmd0TJL2fOnEybNo3Dhw9TpUoVBg8eTPHixfn0009JaCuAUs6S0ERxyzEZwwCISIbkC+nh6agn95UlC4SHQ6FCVlPUN9/YHVHKKF68OGvWrGHz5s1kyZKFNm3aULFiRXbt2mV3aCoNSWiiWOoY9ZRFRLoCm4GPki+sh6NNT+4tRw5r/+3HH7eWJ9+/3+6IUk7NmjXZu3cvc+fO5ezZs1SsWJHWrVvzww8u1U2o3FRCV499D/gM+BwoBgwzxkxJzsAehtYo3F+uXLBlC2TKBHXqwPHjdkeUcjw9PenUqRMnTpxg2LBhrFmzhmLFitGnTx9++eUXu8NTbizBo57+dZO1sF8bY8wi54eUdLpntvs7ccLaf9vLy1qmvGBBuyNKeefOnWPkyJHMnTsXb29v+vbty2uvvUa2bNnsDk2lUg816klEMovI6yIyTUTqiKU3cBprboVStnjySdi0ydrwqGZNSIubyOXOnZtZs2Zx/PhxmjdvzjvvvEPBggUZNWqUrh+lnOpBTU9hWE1Nh4AuwEagJRBqjAlN5tiUilepUlYH98WL1o55Fy7YHZE9ChcuTFhYGAcPHqRmzZoMGzaMQoUKMWHCBP7880+7w1Nu4EGJopAxpqMxZhbQFmvjoEbGGJfsRtQ+irSnbFn44gs4e9bqs7hyxe6I7PPUU0+xfPlydu/eTXBwMAMHDqRIkSLMnDmTW7du2R2eSsUelChu331gjPkb+MEY47J1Wh31lDZVqQIrVsCRI9ZoqOvX7Y7IXmXLliU8PJxt27ZRsGBBevToQfHixQkLC9NFB9VDeVCiCBSRq47jGhBw97GIXE2JAJVKiLp14dNPYc8eCA21lv5I66pWrcqOHTtYt24dWbJkoX379gQEBLB8+XKdtKcSJd5EYYzxNMZkdhyZjDFeMR5nTqkglUqIZs1g/nxrXagWLUBbW0BEqF+/PhERESxbtgxjDM2bN6d58+ZcvHjR7vBUKuFGOxMrBe3aWcuSr10LL7wA2tJi8fDwoEWLFhw6dIh33nmHL774glKlSrF+/Xq7Q1OpgFslCu3MVgDdu8N771mrzXbtCtHRdkfkOjw9PXn11VfZs2cPOXLkoEGDBvTq1Ys//vjD7tCUC3OrRKGd2equV16B4cNh3jzo1889999OisDAQPbs2cOAAQOYMWMGQUFB7Nmzx+6wlItyq0ShVEzDh1sJY+pUGDLE7mhcj5+fHxMmTGDLli388ccfVKxYkVGjRnHnzh27Q1MuRhOFclsi1p7b3bvD2LHWof6rRo0aHDx4kFatWjFs2DCeeeYZIiMj7Q5LuRBNFMqticCMGVYn9xtvWLUL9V9Zs2Zl0aJFLF68mOPHj1O6dGk+/PBDHUargFSSKEQkg4jsFZFGdseiUh8PD6uvolkz6NsX5s61OyLX1aZNGw4dOkSFChXo1q0bTZs25fz583aHpWyWrIlCROaKyHkRORzrfD0R+V5EIkVkcAKKGgQsTZ4oVVrg5QWLF1sT87p2tSbnqbjlyZOHjRs3MnHiRMLDw3nqqadYs2aN3WEpGyV3jWI+UC/mCccS5dOB+oA/0FZE/EWklIh8Eet4VERqAUeBX5M5VuXmfH1h+XKoXNlqitLvvvvz8PCgX79+RERE8MQTT9CkSRMaNGjAjh077A5N2SBZE4UxZjtwOdbpckCkMea0MeYWsARrNdpDxphGsY7zQHWgAvAc0FVEUkVzmXJN6dNbCSIoCFq2tDZBUvf31FNP8e233zJ27FgiIiKoUqUKzzzzDOvXr9f+izTEji/d3MBPMZ5HOc7FyRgzxBjTD/gE+NAYE+f0KRHpJiIRIhJxIa2uN60SJHNm2LDB2tOiSRPYudPuiFybr68vgwcP5syZM0yZMoWzZ8/SoEEDgoODWbp0qS40mAbYkSgkjnMP/NPEGDPfGPNFPK/PBkYA+3x8fJIQnkoLsmWzNj7KkwcaNIB9++yOyPWlT5+ePn36EBkZybx58/jzzz9p3bo1JUqUYM6cObqUuRuzI1FEAXljPM8D/OyMgnVmtkqMxx6DzZsha1ZrL4sjR+yOKHXw8fGhY8eOHDlyhM8++4xMmTLRpUsXChcuzOTJk7lx44bdISonsyNR7AGKikhBEfEB2gCrnVGwrvWkEitvXitZ+PhA7dpw6pTdEaUenp6eNG/enIiICMLDwylcuDD9+vUjf/78jB49mt9++83uEJWTJPfw2MXALqCYiESJSGdjzB2gNxAOHAOWGmP0bzllmyJFrGRx65a1//ZPPz34HvUPEaFOnTps27aNr776igoVKjB06FDy58/P+PHjtQ/DDYg7jlwICQkxERERdoehUpl9+6B6dXj8cdi+3WqaUg/nwIEDDB8+nFWrVlG9enXCwsLInfu+Y1aUixCRvcaYkNjn3WqoqTY9qaQIDoZ16yAqymqGuhx7YLdKsMDAQFasWMGcOXP49ttvCQwMZNWqVXaHpR6SWyUK7cxWSVWpEqxaBd9/D/XqwVXd8PehiQgvvvgi+/btI1++fDRt2pRevXrx559/2h2aSiS3ShRao1DOUKsWfPYZfPcdNGoEuqdP0hQrVoxdu3bxyiuvMGPGDMqVK8fhw4cffKNyGW6VKLRGoZylcWNYuBC++gqefRb++svuiFI3X19f3nvvPTZs2MCFCxcICQlh+vTpOrs7lXCrRKGUM7VuDR99BOHh0LYt6H4+SVe3bl0OHDhA9erV6d27N02bNuXixYt2h6UewK0ShTY9KWd78UWYPBlWrIBOnXT/bWd47LHHWLt2LRMnTmTDhg0EBgby5Zdf2h2WiodbJQptelLJoW9fGD3aaorq1Uv333aGu6vTfvPNN2TKlIlatWrx+uuvc/v2bbtDU3Fwq0ShVHJ54w0YPBhmzoRXX9Vk4SxBQUHs3buXzp07M27cOCpVqqTbsLogTRRKJYAIjBkDvXvDhAkwcqTdEbmPDBky8OGHH7J06VJOnjxJyZIl6du3L7/88ovdoSkHt0oU2kehkpOI1V/RsSO89ZaVMJTztGzZkkOHDtG+fXtmzJhBoUKFeO2117Sz2wW4VaLQPgqV3Dw84MMPoUULGDgQZs2yOyL3kidPHj788EOOHz9O8+bNee+99yhYsCDDhg3jypUrdoeXZrlVolAqJXh5waJF1j4WPXpYndzKuYoUKUJYWBiHDx+mXr16jBo1ikKFCjF27FiuX79ud3hpjiYKpR6Cj481e7taNaspasUKuyNyT/7+/ixbtox9+/ZRqVIl3njjDQoVKsTEiRN1KZAUpIlCqYeULp21LlRIiDU5Lzzc7ojcV1BQEGvWrGHXrl0EBgYyYMAAihQpwgcffKA766UAt0oU2pmtUlqmTLB+Pfj7Q7Nm1vLkKvlUqFCBTZs2sXXrVgoWLEjPnj158sknmT9/PtE6GzLZuFWi0M5sZYesWWHjRsif31pEcM8euyNyf9WqVWPHjh1s2LCBnDlz0qlTJ6pXr873339vd2huya0ShVJ2efRR2LQJsmeHunXh0CG7I3J/IkLdunXZvXs3c+bM4eDBgwQGBjJmzBid4e1kmiiUcpI8eWDLFqvvonZtOHHC7ojShrv7Xhw7dozGjRszZMgQQkJC0F0unUcThVJOVKiQtf/2339b+1qcPWt3RGnH448/zrJly1ixYgUXL16kfPnyDBw4kBs3btgdWqrn8olCRKqJyA4RmSki1eyOR6kHKVHC6rO4ehVq1oT//c/uiNKWpk2bcvToUbp27cqECRMoVaoUmzdvtjusVC1ZE4WIzBWR8yJyONb5eiLyvYhEisjgBxRjgOuAHxCVXLEq5UxBQdZoqF9+sZqhdBWKlPXII48wc+ZMtm3bhpeXF7Vr16ZTp05c1o3QH0py1yjmA/VinhART2A6UB/wB9qKiL+IlBKRL2IdjwI7jDH1gUHAiGSOVymnefppWL0aIiOt/bd11HbKq1q1KgcOHOD1118nLCyMEiVKsHTpUt1ZL5GSNVEYY7YDsVN4OSDSGHPaGHMLWAKEGmMOGWMaxTrOG2PuDo7+DfBNzniVcrYaNeDzz+HAAWjYELS5POWlS5eOMWPGEBERQd68eWndujWhoaFERWkDRULZ0UeRG/gpxvMox7k4icizIjILCAOmxXNdNxGJEJGICxcuOC1YpZKqYUP45BPYtcualHfzpt0RpU2lS5fmm2++4d1332Xz5s34+/szceJEHUqbAHYkConj3H3rgcaY5caY7saY1saYbfFcNxuraWqfj49P0qNUyolatoQ5c6y5Fm3agH432cPLy4uBAwdy6NAhKlasyIABAwgMDGTjxo12h+bS7EgUUUDeGM/zAD/bEIdSKapjR5g2zVofqkMHawitskfhwoVZv349q1ev5tatW9StW5emTZty6tQpu0NzSXYkij1AUREpKCI+QBtgtTMK1iU8lKvr1QvGjYPFi+Gll3RLVTuJCI0bN+bIkSOMHTv2XnPUkCFDdCnzWJJ7eOxiYBdQTESiRKSzMeYO0BsIB44BS40xR5z0froooHJ5gwbBkCHw0UcwYIAmC7v5+voyePBgTpw4QatWrRgzZgzFihVj0aJFOjrKQdzxgwgJCTE6fV+5MmOgf39ra9WhQ3UPblfy9ddf07dvX/bu3UulSpWYMmUKwcHBdoeVIkRkrzEmJPZ5l5+ZnRhao1CphQhMnAidO8OoUfDOO3ZHpO6qWLEiu3fv5qOPPuLEiROEhITQrVs30vJoSrdKFNpHoVITEWvP7TZtrOaoGTPsjkjd5eHhQefOnTlx4gT9+vVj3rx5FC1alEmTJqXJ4bRulSi0RqFSG09PWLAAmjSxOroXLLA7IhVTlixZeP/99zl48CDly5enf//+BAQEsG7dujTVf+FWiUJrFCo18vaGTz+1Vpvt1Mmaya1cS4kSJdiwYQOrVq3i77//pmHDhtSrV48jR5wyDsfluVWiUCq18vODlSut9aHatoV16+yOSMUmIjRp0oTDhw/z/vvvs3v3bgICAujZs6fb91+4VaLQpieVmmXIAGvXQqlS0Lw5bNtmd0QqLj4+PvTv35+TJ0/Ss2dPZs+eTdGiRZkwYQK3bt2yO7xk4VaJQpueVGr3yCMQHm5tgNS4MXz7rd0RqfvJkSMHU6dO5eDBg1SsWJGBAwfi7+/PypUr3a7/wq0ShVLuIEcOa5e8xx6zlic/cMDuiFR8/P39WbduHRs2bMDX15dmzZpRo0YN9u/fb3doTqOJQikXlCuXtf92pkzWxkfHj9sdkXqQunXrcuDAAWbMmMGhQ4cIDg6mS5cu/PLLL3aHlmRulSi0j0K5k/z5rZqFh4c1IuqHH+yOSD2Il5cXPXr0IDIykgEDBrBgwQKKFi3K+PHjU/X8C7dKFNpHodzNk09aS5P/8YeVLM6dszsilRBZsmThvffe4+jRo9SsWZPBgwdTpkwZvk2lnU5ulSiUckelSlkd3BcuWMnCzUdiupUiRYqwcuVKVq5cyeXLl3n66afp06cPV69etTu0RNFEoVQqULYsfPEFnD0LderAlSt2R6QSIzQ0lGPHjtGnTx+mT59+b3RUaqGJQqlUokoVWLECjhyB+vVBt0xIXTJlysTkyZP55ptvyJ49O82aNePZZ5/lXCpoT3SrRKGd2crd1a1rLfexZ4+1PtSff9odkUqscuXKERERwfjx41m/fj0lSpRgxowZREdH2x3afblVotDObJUWNGsG8+dbM7dbtgQ3nQzs1ry9vXnttdc4fPgwFSpUoFevXlSqVIlDhw7ZHVqc3CpRKJVWtGsHH3xgLfnxwgu6/3ZqVbhwYcLDwwkLCyMyMpLg4GCGDBnCny5WVdREoVQq1b07vPceLF0KXbqAC7dcqHiICO3atePYsWO0a9eOMWPGEBAQwObNm+0O7R5NFEqlYq+8AsOHW01RL7+s+2+nZjly5GDevHls2bIFgNq1a9O4cWOOHj1qc2SpIFGIiIeIvC0iU0Wkg93xKOVqhg+3Esa0aTBkiN3RqKSqUaMGBw8eZNy4cWzfvp1SpUrRrVs3fv75Z9tiStZEISJzReS8iByOdb6eiHwvIpEiMvgBxYQCuYHbQFRyxapUaiUC775rNUWNHQtjxtgdkUqqdOnSMWjQIE6dOkXfvn2ZP38+RYsWZdiwYVy7di3F40nuGsV8oF7MEyLiCUwH6gP+QFsR8ReRUiLyRazjUaAYsMsYMwDokczxKpUqiVh7bj//vFWrmDLF7oiUM+TIkYOJEydy/PhxmjRpwqhRoyhSpAgzZsxI0bWjkjVRGGO2A5djnS4HRBpjThtjbgFLgFBjzCFjTKNYx3msWsRvjnt1bIdS9+HhYfVVNGtm9VfMnWt3RMpZChUqxOLFi/n2228pXrw4vXr14qmnnmLFihUpsveFHX0UuYGfYjyPcpy7n+VAXRGZCmy/30Ui0k1EIkQkwt23JVTqfry8YPFia2Jely7W5DzlPsqVK8e2bdtYs2YNnp6ePPvss1SuXJmvv/46Wd/XjkQhcZy7b0o0xvxhjOlsjOljjJkez3WzgRHAPh8fHyeEqVTq5OsLy5dD5crWfIs1a+yOSDmTiNCoUSMOHjzI7NmzOX36NJUqVaJFixacOHEiWd7TjkQRBeSN8TwPYF93vlJuKH16axHBoCBr9rZjxKVyI15eXnTt2pXIyEhGjhxJeHg4JUuWJDw83OnvZUei2AMUFZGCIuIDtAFWO6NgXcJDqX9kzgwbNlh7WjRpAjt32h2RSg4ZMmRg6NChREZG8sorr1ClShWnv0dyD49dDOwCiolIlIh0NsbcAXoD4cAxYKkx5oiT3k8XBVQqhmzZrI2PcueGBg1g3z67I1LJ5bHHHmPcuHGkS5fO6WVLSvSYp7SQkBATERFhdxhKuYwff4RnnoEbN2D7dvD3tzsi5YpEZK8xJiT2eZefmZ0YWqNQKm758ln9FN7e1i55p0/bHZFKTdwqUWgfhVL3V6QIbN4MgYGg/4uoxPCyOwBnEpHGQOMiRYrYHYpSLqlkSVi/3u4oVGqjNQqllFLxcqtEoZRSyvncKlFoZ7ZSSjmfWyUKbXpSSinnc6tEoZRSyvk0USillIqXWyUK7aNQSinnc6tEoX0USinlfG651pOIXADOPuTtOYCLTgwntdLPwaKfwz/0s7C48+eQ3xiTM/ZJt0wUSSEiEXEtipXW6Odg0c/hH/pZWNLi5+BWTU9KKaWcTxOFUkqpeGmi+K/ZdgfgIvRzsOjn8A/9LCxp7nPQPgqllFLx0hqFUkqpeGmiUEopFS9NFDGISD0R+V5EIkVksN3x2EVEzojIIRHZLyJpZvNxEZkrIudF5HCMc9lEZJOInHT8m9XOGFPCfT6Ht0TknON3Yr+INLAzxpQgInlFZKuIHBORIyLysuN8mvud0EThICKewHSgPuAPtBWRtLwFfXVjTOk0Nl58PlAv1rnBwBZjTFFgi+O5u5vPfz8HgImO34nSxph1KRyTHe4ArxhjSgAVgF6O74Q09zuhieIf5YBIY8xpY8wtYAkQanNMKgUZY7YDl2OdDgU+djz+GGiaokHZ4D6fQ5pjjPmfMWaf4/E14BiQmzT4O6GJ4h+5gZ9iPI9ynEuLDLBRRPaKSDe7g7HZY8aY/4H1xQE8anM8duotIgcdTVNu39wSk4gUAIKAb0mDvxOaKP4hcZxLq2OHKxljgrGa4XqJSBW7A1K2+wAoDJQG/gdMsDeclCMiGYHPgX7GmKt2x2MHTRT/iALyxnieB/jZplhsZYz52fHveWAFVrNcWvWriOQCcPx73uZ4bGGM+dUY87cxJhr4kDTyOyEi3lhJYpExZrnjdJr7ndBE8Y89QFERKSgiPkAbYLXNMaU4EckgIpnuPgbqAIfjv8utrQY6OB53AFbZGItt7n4xOjQjDfxOiIgAc4Bjxpj3Y7yU5n4ndGZ2DI4hf5MAT2CuMeZtm0NKcSJSCKsWAeAFfJJWPgcRWQxUw1pG+ldgOLASWArkA34EWhpj3Lqj9z6fQzWsZicDnAG6322nd1ciUhnYARwCoh2n38Dqp0hbvxOaKJRSSsVHm56UUkrFSxOFUkqpeGmiUEopFS9NFEoppeKliUIppVS8NFEoFQ8Rue74t4CIPOfkst+I9fxrZ5avlLNoolAqYQoAiUoUjhWJ4/OvRGGMqZjI8KpbswAAAfJJREFUmJRKEZoolEqYccAzjr0Y+ouIp4i8KyJ7HAvldQcQkWqOPQw+wZqohYisdCyweOTuIosiMg5I5yhvkePc3dqLOMo+7NgXpHWMsreJyGciclxEFjlmDyuVrLzsDkCpVGIwMNAY0wjA8YX/uzGmrIj4AjtFZKPj2nLAU8aYHxzPXzTGXBaRdMAeEfncGDNYRHobY0rH8V7PYs2CDsSaHb1HRLY7XgsCSmKtQ7YTqAR85fwfV6l/aI1CqYdTB2gvIvuxlnTIDhR1vLY7RpIA6CsiB4BvsBaeLEr8KgOL/7+9O8SJGIjiMP49hyEorkC4AGINQXAHBBcAgeUeWBQ3wLIOkIgN1OEREASCQEKAPESnSWlg0mxg1fdTnUk6SUXz77wmb0oTvgfgAtjorX1XmvNd05bEpH/ljkKaTwAHmTn9NhmxBbwMxtvAJDNfI+IcWBqx9m/eetef+A5rAdxRSOM8A8u98RTYL22oiYi10m13aAV4KiGxTnukZue9u3/gEtgp/0FWgU3g6k+eQpqDXyPSOA3wUUpIJ8ARbdlnVn4oP/LzkZhnwF5ENMAtbfmpcww0ETHLzN3e/CkwAW5ou7UeZuZ9CRpp4eweK0mqsvQkSaoyKCRJVQaFJKnKoJAkVRkUkqQqg0KSVGVQSJKqvgClCMKF4ydP7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "A = np.array([[10., 2., 3., 5.],\n",
    "                 [1., 14., 6., 2.],\n",
    "                 [-1., 4., 16.,-4],\n",
    "                 [5. ,4. ,3. ,11.]])\n",
    "b = np.array([1., 2., 3., 4.])\n",
    "# an initial guess at the solution - here just a vector of zeros of length the number of rows in A\n",
    "x = np.zeros(A.shape[0]) \n",
    "\n",
    "tol = 1.e-6 # iteration tolerance\n",
    "it_max = 1000 # upper limit on iterations if we don't hit tolerance\n",
    "\n",
    "x_j, res_j = jacobi(A,b)\n",
    "x_gs, res_gs = gauss_seidel(A,b)\n",
    "print(x_j)\n",
    "print(x_gs)\n",
    "\n",
    "# plot the log of the residual against iteration number \n",
    "fig = plt.figure(figsize=(6, 4))\n",
    "ax1 = plt.subplot(111)\n",
    "ax1.semilogy(res_j,'k',label='Jacobi') # plot the log of the residual against iteration number \n",
    "ax1.semilogy(res_gs,'b',label='Gauss-Seidel')\n",
    "ax1.set_xlabel('Iteration')\n",
    "ax1.set_ylabel('Residual')\n",
    "ax1.set_title('Convergence plot')\n",
    "ax1.legend(loc='best')\n",
    "plt.show()\n",
    "\n",
    "# check our solutions agrees with multiplying through by the inverse matrix\n",
    "# [0] as our implemntations also return the residuals\n",
    "print(np.allclose(sl.inv(A)@b, jacobi(A, b)[0]))\n",
    "print(np.allclose(sl.inv(A)@b, gauss_seidel(A, b)[0]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "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.3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": false,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "213.809px"
   },
   "toc_section_display": false,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}