EquationSolvers_temp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Linear equation system solving techniques comparison\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%% \n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import datetime as dt\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "<h4> \n",
    "\n",
    "    <i>\n",
    "        index = 172110 </br>\n",
    "        c = 1\n",
    "        </br>\n",
    "        d = 0\n",
    "        </br>\n",
    "        e = 1 \n",
    "        </br>\n",
    "        f = 2\n",
    "        </br>\n",
    "        9cd = 910\n",
    "    </i>\n",
    "</h4>\n",
    "\n",
    "**Exercise A**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "N = 300\n",
    "iter_limit = 100\n",
    "epsilon = 10.0 ** (-6)\n",
    "f = 2\n",
    "e = 1\n",
    "a1 = a2 = a3 = 0\n",
    "\n",
    "def init_a_exercise():\n",
    "    global a1, a2, a3, iter_limit\n",
    "    a1 = 5 + e\n",
    "    a2 = a3 = -1\n",
    "    iter_limit = 50\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Striped matrix generating function**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def generate_striped_matrix(size: int):\n",
    "    array = np.ndarray(shape=(size,size))\n",
    "    # wybór elementu na diagonali (i,i)\n",
    "    for i in range(size):\n",
    "        # wybór 2 elementy przed, diagonal i 2 elementy za diagonalą\n",
    "        for j in range(size):\n",
    "            if 0 <= j < size:\n",
    "                if abs(i - j) == 2:\n",
    "                    array[i][j] = a3\n",
    "                elif abs(i - j) == 1:\n",
    "                    array[i][j] = a2\n",
    "                elif i == j:\n",
    "                    array[i][j] = a1\n",
    "                else:\n",
    "                    array[i][j] = 0\n",
    "    return array\n",
    "                \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   },
   "source": [
    "**Residue vector second norm calculation**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def calculate_residue_vector_norm(A: np.ndarray, x: np.ndarray, B: np.ndarray):\n",
    "    n_size = A.shape[0]\n",
    "    residue_norm = 0.0\n",
    "    for i in range(n_size):\n",
    "        row = 0.0\n",
    "        for j in range(n_size):\n",
    "            row += A[i][j] * x[j]\n",
    "        row -= B[i]\n",
    "        residue_norm += row ** 2\n",
    "    return math.sqrt(residue_norm)\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Timing functions**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "starting_time: dt.datetime\n",
    "def tic():\n",
    "    global starting_time\n",
    "    starting_time = dt.datetime.now()\n",
    "    \n",
    "def toc() -> dt.timedelta:\n",
    "    time_diff = dt.datetime.now() - starting_time \n",
    "    return time_diff\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Function printing details about tested solution.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def print_details(A: np.ndarray, B: np.ndarray, x: np.ndarray, iteration_count: int, time: dt.timedelta):\n",
    "    print()\n",
    "    print(\"Residue vector second norm value: \", calculate_residue_vector_norm(A, x, B))\n",
    "    print(\"Minutes: \", int(time.seconds / 60), \"seconds: \", time.seconds % 60)\n",
    "    print(\"Iteration count: \", iteration_count)\n",
    "    print(\"Time per iteration in seconds: \", time.seconds / iteration_count)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   },
   "source": [
    "**Matrix initialization**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "B: np.ndarray\n",
    "A: np.ndarray\n",
    "x: np.ndarray\n",
    "\n",
    "def generate_equation_system():\n",
    "    global A\n",
    "    global B   \n",
    "    global x\n",
    "    \n",
    "    # Wektor pobudzenia\n",
    "    B = np.ndarray(shape=(N, 1), dtype=float)\n",
    "    for i in range(N):\n",
    "        B[i] = math.sin(i * (f + 1))\n",
    "    \n",
    "    # Macierz systemowa\n",
    "    A = generate_striped_matrix(N)\n",
    "    \n",
    "    # Macierz rozwiązań przybliżonych. Zawiera początkowe przybliżenie a później przyjmuje wartości nadane przez\n",
    "    # metody iteracyjne.\n",
    "    x = np.ndarray(shape=(N, 1), dtype=float)\n",
    "    for i in range(N):\n",
    "        x[i] = 0\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Exercise B\n",
    "\n",
    "**Jacobi's method**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def jacobis_next_iteration(matrix: np.ndarray, b_vec: np.ndarray, estimates: np.ndarray) -> np.ndarray:\n",
    "    this_iteration: np.ndarray = np.ndarray(shape=(estimates.shape[0],1))\n",
    "    for i in range(matrix.shape[0]):\n",
    "        # Copy of the i-th equation coefficients\n",
    "        curr_equation = matrix[i]\n",
    "        # The coefficient of the x which we're currently we're estimating new value for.\n",
    "        curr_calc_coeff = curr_equation[i]\n",
    "        # Inserting the b vector value into the equation.\n",
    "        new_x = - b_vec[i]\n",
    "        for j in range(curr_equation.shape[0]):\n",
    "            # for every variable except the one we're currently calculating -> (a1 * x1 + ... + an * xn - b)/ aj = xj \n",
    "            if j != i:\n",
    "                # eq[j] is the coefficient of x(j) and estimates[j] is the previous value of the x(j)\n",
    "                new_x += curr_equation[j] * estimates[j]\n",
    "        this_iteration[i] = new_x / - curr_calc_coeff\n",
    "    return this_iteration\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Gauss-Seidl's method.**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def gauss_seidl_next_iteration(matrix: np.ndarray, b_vec: np.ndarray, estimates: np.ndarray) -> np.ndarray:\n",
    "    this_iteration: np.ndarray = np.ndarray(shape=(estimates.shape[0],1))\n",
    "    for i in range(matrix.shape[0]):\n",
    "        # Copy of the i-th equation coefficients\n",
    "        current_equation = matrix[i]\n",
    "        # The coefficient of the x which we're calculating new value for\n",
    "        curr_calc_coeff = current_equation[i]\n",
    "        # Inserting the b vector value into the equation\n",
    "        new_x = - b_vec[i]\n",
    "        for j in range(i):\n",
    "            new_x += current_equation[j] * this_iteration[j]\n",
    "        for j in range(i + 1, matrix.shape[0]):\n",
    "            new_x += current_equation[j] * estimates[j]\n",
    "        this_iteration[i] = new_x / - curr_calc_coeff\n",
    "    return this_iteration\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   },
   "source": [
    "**Testing the Jacobi's method on data from exercise A**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "........................\n",
      "Residue vector second norm value:  8.140195485706073e-07\n",
      "Minutes:  0 seconds:  11\n",
      "Iteration count:  24\n",
      "Time per iteration in seconds:  0.4583333333333333\n"
     ]
    }
   ],
   "source": [
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "# Printing details about this method time, iteration count etc.\n",
    "print_details(A, B, x, i, toc())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**Testing Gauss-Seidl's on data from A exercise**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "................\n",
      "Residue vector second norm value:  7.696797973615508e-07\n",
      "Minutes:  0 seconds:  7\n",
      "Iteration count:  16\n",
      "Time per iteration in seconds:  0.4375\n"
     ]
    }
   ],
   "source": [
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "\n",
    "# residue_vector_norm_values = [float]\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while  residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    # residue_vector_norm_values.append(residue_vector_norm)\n",
    "    print(\".\", end=\"\")\n",
    "    \n",
    "# Printing details about this method time, iteration count etc.\n",
    "print_details(A, B, x, i, toc())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The Gauss-Seidl's method consequently performs better than Jacobi's method. This is true both for iteration count\n",
    "and the time of a single iteration. This is partially due to the method of variables caching by the processor and \n",
    "the faster access the ones used more recently. However less efficient as sequential algorithm, the Jacobi's method\n",
    "is possible to run concurrently. \n",
    " \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Exercise C\n",
    "\n",
    "**Changing the coefficients used to generate equation system to match the specification of exercise C**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "def init_c_exercise():\n",
    "    global a1, a2, a3, iter_limit\n",
    "    a1 = 3\n",
    "    a2 = a3 = -1\n",
    "    iter_limit = 30\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   },
   "source": [
    "**Testing the Jacobi's method on data from exercise C**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..............................\n",
      "Residue vector second norm value:  64.92804780786507\n",
      "Minutes:  0 seconds:  14\n",
      "Iteration count:  30\n",
      "Time per iteration in seconds:  0.4666666666666667\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XecXGXZ//HPd0uySTbJEhJCElJJaKFmQxPFBFAQlaIUG82C+hPlEfUBy/MIVhAVsSKIgg9IQLqotJBQVAhZekhCIL1nk91kN2XbXL8/zj1hdtlyNtnZadf79ZrXzClzznXPmTnXnPs+5z4yM5xzzhWuokwH4JxzLrM8ETjnXIHzROCccwXOE4FzzhU4TwTOOVfgPBE451yB80TQAUnzJE3rYNo0SSt7OSQXk6TZkj6b6TjSSdKFkp7p5XX+U9IFvbnO3SVpjKR6ScUZjOGTkh7N1PrjyPlEIGmppO1hY6+VdIuk8t1drplNNrPZPRBixkkySRMzHYfLbWb2ATO7NdNxpOoq6ZvZcjMrN7OWOPP3QDzjwu+tJCWG283s/elaZ0/I+UQQfNjMyoHDgSOAb2Y4nryR+oV2rWXyX6bLjHzd5vmSCAAws7XAI0QJAQBJfSX9VNJySesk3SCpX5g2VNJDkmolbZL0tKSiMG2ppJPC637hSKNG0uvAkanrbfuPO8z7g5ThD0l6Kazn35IObS/+ENtP24x7QNJl4fVISfdI2iBpiaSvpMxXLOlbkt6SVCepStJoSU+FWV4OR03nhvk/J+nNUO4HJY1sU54vSVoELGonzjJJt0naGMr0vKThYdpgSTdLWiNplaQfpP54wnrnhxhflzQljD8w/FurDdVyp7X5PH8j6e/hfc9J2jdl+vskLZC0WdKvAbX3+YZ5r5R0l6Q/h2XNkzQ1ZXpXcfxO0j8kbQWmh3G/DdUm9ZL+JWlvSb8I35cFko7Yxe19Rcr2fF3SmR0s5x3/QtXmn6+kT4fPvUbSI5LGhvGSdJ2k9eHze0XSwR2sZ+cyFaqmFP22asL38QOdfO5LJX09LH+zpDsllaVMb/c3Imnf8B1Nfk9GSqpWVD37Q+A9wK/DZ//rzj6bjuaXdICkx8J6Fko6J+X97W3zD0p6UdIWSSskXZmyyuTvrTas41i1qcaT9C5Fv5nN4fldbT7j74fvUZ2kRyUNDdM6/N3tNjPL6QewFDgpvN4HeBW4PmX6L4AHgSHAQOBvwI/DtB8DNwCl4fEeQO0s92rg6bCM0cBrwMqUdRgwMWX4FuAH4fUUYD1wNFAMXBCW3bedshwPrEiJYQ9gOzCSKGlXAf8L9AEmAIuBk8O83whl359oR3gYsGcH8Z0AVIfY+gK/Ap5qU57HQnn7tRPn58Pn2D+UqRIYFKbdD/weGADsBcwBPh+mnQ2sIkqkAiYCY8Nn/ybwrVC2E4A6YP+Uz3MTcBRQAtwOzAjThgJbgLPCcr4KNAOf7eD7ciWwAzg1xP5j4NkwLU4cm4HjwvYoC+Oqw2dQBjwBLAHOD8v/ATCrg1g63N4pn1dy258LbAVGhGkXAs+E1+PCNitJWfbs5GcAnBHKdWD4/L4D/DtMO5noe1URtsmByXW0E2/qMi8EmoDPhXJ+EVidLEsHv9M5oTxDgPnAF+L8RsI65hN93x4BftpeTB2st9Vn03Z+ou/pCuCi8NlMCdtzcifbfBpwSBg+FFgHnNHJtkjdVkOAGuC8sL6Ph+E9U+J7C9gP6BeGr+7qd7fb+9FM7Lx78hG+MPVEP1gDZgIVYZqIfjz7psx/LLAkvP4e8AApO8k2y00mgsXAKSnTLiZ+Ivgd8P02y14IvLeddQpYDhyf8gN4Irw+GljeZv5vAn9KWebpHXxGbeO7GfhJynA50Y96XMr8J3TymX8a+DdwaJvxw4EGUpJH+KLPCq8fAS5tZ3nvAdYCRSnj7gCuTPk8/5Ay7VRgQXh9PmFHnvIZrqTzRPB4yvBBwPZuxPHnNsu7BbgpZfjLwPyU4UOA2g5i6XB7dzD/S8ltTPcSwT+Bz6RMKwK2ESXhE4A3gGNSy93B+lOXeSHwZsq0/iGGvTv5nX4qZfgnwA1xfyNEf+ZeBV4h5U8Uu58IzgWebvOe3wPf7Wibt7OOXwDXdbItUrfVecCcNu//D3BhSnzfSZn2/4CHO/vd9cQjX6qGzjCzgUSZ+gCif4kAw4i+oFXhUKoWeDiMB7iW6J/So5IWS7qig+WPJPrXkLSsG7GNBb6WXH+IYXRYZisWbe0ZRDtPgE8Q/ftNLmdkm+V8i2jnS1jmWzFjGplaBjOrBzYCo1LmWdH2TSn+j2inPkPSakk/kVTK2//u16TE+HuiI4POYhwJrDCzRMq4ZW3iWZvyehtR8tr53pSyWBext7esslCtEieO9pa9LuX19naG2z15oYvtjaTzU6pLaoGDefu73R1jgetTlrOJKAmNMrMngF8DvwHWSbpR0qCYy935OZrZtvCysxM1OtqGcX4jNxGV/1dm1hAzvjjGAke3Wfcngb1T5mm1zSUdLWmWoirazcAXiL9dWv32grjf9Y5+d7stXxIBAGb2JFEGT9a7VhP9ECebWUV4DLaoYRkzqzOzr5nZBODDwGWSTmxn0WuIvphJY9pM30aUcJLafol+mLL+CjPrb2Z3dFCMO4CzQh3u0cA9KctZ0mY5A83s1JTp+7azvPasJvoBACBpALAnUbVNknX0ZjNrMrOrzOwg4F3Ah4j+ma8gOiIYmhLjIDOb3EWMq4HRCu0zwZg28XSk1baRJFpvq+6IE0eHn8suand7h+GbgEuIqg0qiKok22v/2BqeO/sOfr7Nd6efmf0bwMx+aWaVwGSiKolv9FzxYun0N6LoLMBfEB3JXilpSMp7u7s92s6/AniyzbrLzeyLnbznL0RHKKPNbDBR9bI6mLetVr+9INZ3vZPf3W7Lq0QQ/AJ4n6TDwz+7m4DrJO0FIGmUpJPD6w9Jmhh2HluAlvBo6y7gm5L2kLQP0eF/qpeATyhqsD0FeG/KtJuAL4R/EZI0IDQ2DWwveDN7EdgA/AF4xMxqw6Q5wBZJlytqvC6WdLCkZMP1H4DvS5oU1nOopD3DtHVEbQpJfwEuknS4pL7Aj4DnzGxpRx9qKknTJR2iqBF4C1G1UouZrQEeBX4maZCkIkWNfcnP4w/A1yVVhhgnhh3ec0Q7s/+WVKro+o0PE/1b7srfgcmSPhL+1X+F1jvB7tidOHZJJ9t7ANFOZQOApIuI/hG3t4wNRDuST4XvxadpnXBvIPr+Tg7LGizp7PD6yPDdLCUq+w7a/w2kU1e/keuBKjP7LNH2viHlvW2/211pO/9DwH6SzgvbvDR8Jgd2soyBwCYz2yHpKKIjuaQNQKKTmP4R1vcJRQ3Y5xJVTz7UVeAd/e66el8ceZcIwo/iz8D/hFGXE1X/PCtpC/A4UYMqwKQwXE9UT/dba//agauIDt+WEO3o/q/N9EuJdhjJw8r7U+KZS1T3+2uiRqE3ieoMO3MHcBLRDju5nJawjsNDHNVEO4/BYZafEyWsR4m+JDcTNTZBVC9+azj0PcfMZhJ9PvcQ/aPeF/hYFzGl2hu4O6xnPvAkcFuYdj5RQ+vrobx3AyNCGf4K/DCUq47ocxpiZo3AacAHQrl+C5xvZgu6CsTMqokaVa8mqt6aBPyrG2VJXdYux7Gb2tverwM/I/periNqa+isXJ8j+ie/keif/b9TlnUfcA1RlcIWoiOL5Bk+g4h2xDVE3/GNvH1E3Ss6+41IOh04haj6BeAyYIqkT4bh64mOqGok/TLG6lrNb2Z1wPuJvv+riaplriE6iaIj/w/4nqQ6opM37kopyzai7/i/wu/tmDZl3Uj0T/5rRJ/1fwMfCt/jrnT2u9stybMVnHPOFai8OyJwzjnXPZ4InHOuwHkicM65AueJwDnnClxOdCg2dOhQGzduXKtxW7duZcCAAZkJKA3yrTyQf2Xy8mS/fCvT7panqqqq2syGdTVfTiSCcePGMXfu3FbjZs+ezbRp0zITUBrkW3kg/8rk5cl++Vam3S2PpFi9IHjVkHPOFThPBM45V+A8ETjnXIHzROCccwXOE4FzzhU4TwTOOVfgPBE451yB80TgnHNZ6K0N9dy7qJH1W3akfV2eCJxzLgs9/cYGHnyriaZE+m8V4InAOeey0NxlNQwpE6Mq+nU9827yROCcc1nohWU1TKzonV20JwLnnMsyq2u3s3rzDibtUdwr6/NE4JxzWWbushoAJvkRgXPOFaaqpZvo36eY0QM9ETjnXEGqWl7D4aMrKC5Sr6zPE4FzzmWRrQ3NzF9Tx9Sxe/TaOj0ROOdcFnlpRS0tCWOKJwLnnCtMc5fWIOGJwDnnClXV8hr2Hz6QQWWlvbZOTwTOOZclWhLGi8tqqOzFowFIcyKQVCHpbkkLJM2XdKykIZIek7QoPPduiZ1zLku9sa6Ouobm/EoEwPXAw2Z2AHAYMB+4AphpZpOAmWHYOecKXvJCsqljh/TqetOWCCQNAo4HbgYws0YzqwVOB24Ns90KnJGuGJxzLpe8sKyGYQP7MnpI+juaSyWz9HRxKulw4EbgdaKjgSrgUmCVmVWkzFdjZu84DpJ0MXAxwPDhwytnzJjRanp9fT3l5eVpiT0T8q08kH9l8vJkv1wv0zee3MbYQUVcckQZsPvlmT59epWZTe1yRjNLywOYCjQDR4fh64HvA7Vt5qvpalmVlZXW1qxZs94xLpflW3nM8q9MXp7sl8tlWrd5u429/CG76am3do7b3fIAcy3G/jqdbQQrgZVm9lwYvhuYAqyTNAIgPK9PYwzOOZcTku0Dvd1QDGlsIzCztcAKSfuHUScSVRM9CFwQxl0APJCuGJxzLldULauhb0kRk0cO7vV1l6R5+V8GbpfUB1gMXESUfO6S9BlgOXB2mmNwzrmsN3dZDYeNrqBPSe9f3pXWRGBmLxG1FbR1YjrX65xzuWR7YwvzVm3mc8dPyMj6/cpi55zLsJdX1tKcsF7tcTSVJwLnnMuwqgw2FIMnAuecy7iqZTVM3Kuciv59MrJ+TwTOOZdBiYRRtayGyjGZ63bNE4FzzmXQWxvq2by9icpxngicc64gVe3saM4TgXPOFaS5y2oYMqAP44cOyFgMngiccy6DqpbVMGXMHkjKWAyeCJxzLkOq6xtYUr2VqRlsH4AYVxZLmgp8Gxgb5hdgZnZommNzzrm89kIWtA9AvC4mbge+AbwKJNIbjnPOFY6qZTX0KS7i4FG939FcqjiJYIOZPZj2SJxzrsDMXVbDIfsMpqy0OKNxxEkE35X0B6L7CzckR5rZvWmLyjnn8tyOphZeXbmZC48bl+lQYiWCi4ADgFLerhoywBOBc87totdWbaaxJZGx/oVSxUkEh5nZIWmPxDnnCkimO5pLFef00WclHZT2SJxzroDMXVbD+KEDGFreN9OhxDoieDdwgaQlRG0Efvqoc87tBjPjhWU1TNt/r0yHAsRLBKekPQrnnCsgS6q3snFrY8YvJEvqNBFIKgL+bmYH91I8zjmX97Kho7lUnbYRmFkCeFnSmF6Kxznn8l7VshoG9ytl32HlmQ4FiFc1NAKYJ2kOsDU50sxOS1tUzjmXx+Yuq2HKmAqKijLX0VyqOIngqrRH4ZxzBaJ2WyNvrq/nzCNGZTqUnbpMBGb2pKThwJFh1BwzWx9n4ZKWAnVAC9BsZlMlDQHuBMYBS4FzzKym+6E751zueWF59lw/kNTldQSSzgHmAGcD5wDPSTqrG+uYbmaHm9nUMHwFMNPMJhF1W3FFN2N2zrmcNXdpDSVF4rB9KjIdyk5xqoa+DRyZPAqQNAx4HLh7F9d5OjAtvL4VmA1cvovLcs65nDJ3WQ2TRw6iX5/MdjSXSmbW+QzSq6ldTIRTSl+O0+1EuAithqhvot+b2Y2Sas2sImWeGjN7xzGSpIuBiwGGDx9eOWPGjFbT6+vrKS/Pjhb3npBv5YH8K5OXJ/tle5m2NxuXzNzGyeNKOWf/Pl3Ov7vlmT59elVKbUzHzKzTB3At8AhwYXj8E7imq/eF944Mz3sBLwPHA7Vt5qnpajmVlZXW1qxZs94xLpflW3nM8q9MXp7sl+1leujl1Tb28odszpKNsebf3fIAcy3GvjpOY/E3JH0UOI6oe4kbzey+ONnIzFaH5/WS7gOOAtZJGmFmaySNAGI1PDvnXK6bOX8dFf1LOWJ09rQPQLw2AszsHuCe7ixY0gCgyMzqwuv3A98DHgQuAK4Ozw90K2LnnMtBLQlj1sL1TN9/L0qKs+t28XHuWfwR4Bqi6h3xdqdzg7p463DgPknJ9fzFzB6W9Dxwl6TPAMuJzkZyzrm89uLyGmq2NXHigdnR0VyqOEcEPwE+bGbzu7NgM1sMHNbO+I3Aid1ZlnPO5brH56+npEgcv9+wTIfyDnGOT9Z1Nwk455xr7YkF6zhq/BAGlZVmOpR3iHNEMFfSncD9+D2LnXOu21Zs2sYb6+o598js7L8zTiIYBGwjauxN8nsWO+dcTDPnrwPgxAOyr30A4vU1dFFvBOKcc/lq5oL17DtsAOOGDsh0KO3KrnOYnHMuz9TtaOLZxRs58cDhmQ6lQ54InHMujZ5ZVE1Ti2VttRB4InDOubR6fP56Bvcrzapup9vqsI1A0mWdvdHMft7z4TjnXP5oSRizF65n2v7Dsu5q4lSdNRYPDM/7E92U5sEw/GHgqXQG5Zxz+eClFbVs3NqY1e0D0EkiMLOrACQ9Ckwxs7owfCXw116JzjnnctjM+esoLhLvzcKriVPFOVYZAzSmDDcS3WbSOedcJ55YsJ4jx+3B4H7ZdzVxqjgXlP0fMCd0I23AmcCf0xqVc87luJU121iwto7vfPDATIfSpTgXlP1Q0sPAu8Ooi8zsxfSG5ZxzuW3m/OhWKydk8WmjSbHuRwC8BKxJzi9pjJktT1tUzjmX42YuWM+EoQOYMCx7b52ZFOd+BF8GvgusA1oI9yMADk1vaM45l5vqG5p59q2NXPCusZkOJZY4RwSXAvuH+wg455zrwjOLqmlsSXDCAdl92mhSnLOGVgCb0x2Ic87li5nz1zGorISp47L3auJUcY4IFgOzJf2d1vcj8CuLnXOujUS4N/G0/feiNIuvJk4VJxEsD48+4eGcc64DL6+spbq+MSvvTdyROKePJq8wHhgNWn3ao3LOuRw1c/76nLiaOFWXxy2SDpb0IvAaME9SlaTJ6Q/NOedyz8wF66kcuwcV/XOnAiVOBdaNwGVmNtbMxgJfA26KuwJJxZJelPRQGB4v6TlJiyTdKSl3Pi3nnOvEqtrtzF+zhZNyqFoI4iWCAWY2KzlgZrOB7txv7VJgfsrwNcB1ZjYJqAE+041lOedc1noieW/iLO9ttK04iWCxpP+RNC48vgMsibNwSfsAHwT+EIYFnADcHWa5FTij+2E751z2mblgPeP27M+ELL03cUdkZp3PIO0BXMXbfQ09BVxlZjVdLly6G/gx0b0Nvg5cCDxrZhPD9NHAP83s4HbeezFwMcDw4cMrZ8yY0Wp6fX095eXZf+l2XPlWHsi/Mnl5sl8my9TQbHzpiW2cOLqEjx/Yt0eWubvlmT59epWZTe1yRjNLywP4EPDb8Hoa8BAwDHgzZZ7RwKtdLauystLamjVr1jvG5bJ8K49Z/pXJy5P9Mlmmh19bY2Mvf8j+tWhDjy1zd8sDzLUY++s4Zw09JqkiZXgPSY/ESEbHAadJWgrMIKoS+gVQISl52uo+wOoYy3LOuaz2xPz1DCwr4cjxQzIdSrfFaSMYama1yQGLqoS6bBI3s2+a2T5mNg74GPCEmX0SmAWcFWa7AHig21E751wWSSSMmQvW8979huXM1cSp4kSckDQmOSBpLFHvo7vqcuAySW8CewI378aynHMu415ZtZnq+oacupo4VZwuJr4NPCPpyTB8PKERNy6LTjmdHV4vBo7qzvudcy6b3f/iKvqUFDF9/zxNBGb2sKQpwDFE9yL4qplVpz0y55zLATuaWrj3hZWcMnnvnLqaOFWcxmIBpwBTzOxvQH9J/o/eOeeAR+atZcuOZj525OhMh7LL4rQR/BY4Fvh4GK4DfpO2iJxzLofMmLOCMUP6c8yEPTMdyi6LkwiONrMvATtg51lDuXn845xzPWjZxq38Z/FGzpm6D0VFynQ4uyxOImiSVEw4U0jSMCCR1qiccy4H3DV3BUWCsypzt1oI4iWCXwL3AXtJ+iHwDPCjtEblnHNZrrklwV/nrmT6/nux9+CyTIezW+KcNXS7pCrgRKKzhs4ws/ldvM055/La7IUbWF/XwLk53EicFOesoX2BJWb2G6Kb07wvtcsJ55wrRDOeX8GwgX2ZfkBuXjuQKk7V0D1Ai6SJRN1Jjwf+ktaonHMui63bsoNZC9fz0Sn75GSXEm3F6mLCzJqBjwDXm9lXgRHpDcs557LX3VUraUlYXlQLQfyzhj4OnE/UlTRAafpCcs657GVm3DV3BUePH8L4HLsBTUfiJIKLiC4o+6GZLZE0HrgtvWE551x2enbxJpZt3MbHjsqPowGId9bQ68BXUoaXAFenMyjnnMtWdz6/nIFlJXzg4PypIc/9Vg7nnOslm7c18Y/X1nLG4aMoKy3OdDg9xhOBc87FdP9Lq2hsTuRNI3FSp4lAUrGka3srGOecy1Zmxh1zlnPwqEEcPGpwpsPpUZ0mAjNrASpDV9TOOVewXl21mQVr6zj3yDFdz5xj4tyh7EXgAUl/BbYmR5rZvWmLyjnnssydz6+grLSI0w4bmelQelycRDAE2AickDLOAE8EzrmCsK2xmQdfWs2pB49gcL/8u4wqzumjF/VGIM45l63+8epa6hqa866ROClOp3P7SLpP0npJ6yTdI2mf3gjOOeeywZ3PL2f80AEcNX5IpkNJizinj/4JeBAYCYwC/hbGOedc3ntzfT3PL63h3CNHk6/nzcRJBMPM7E9m1hwetwDDunqTpDJJcyS9LGmepKvC+PGSnpO0SNKdkvy2l865rHXX3BWUFImPTBmV6VDSJk4iqJb0qXBNQbGkTxE1HnelATjBzA4DDgdOkXQMcA1wnZlNAmqAz+xq8M45l06NzQnufWElJxywF3sNzO27kHUmTiL4NHAOsBZYA5xF1BFdpyxSHwZLw8OIzj66O4y/FTijmzE751yveGLBOqrrG/Oqg7n2yMw6n0E6zsz+1dW4Dt5bDFQBE4HfANcCz5rZxDB9NPBPMzu4nfdeDFwMMHz48MoZM2a0ml5fX095eXlXIeSMfCsP5F+ZvDzZr6fL9NO5O1hVl+Cn7+1HcVHvtw/sbnmmT59eZWZTu5zRzDp9AC/EGdfFMiqAWcB7gDdTxo8GXu3q/ZWVldbWrFmz3jEul+Vbeczyr0xenuzXk2V6aXmNjb38Ifvl42/02DK7a3fLA8y1GPvoDq8jkHQs8C5gmKTLUiYNArrV7Z6Z1UqaDRwDVEgqseiuZ/sAq7uzLOec6w3XPf4GFf1LufC4cZkOJe06ayPoA5QTXXQ2MOWxhaidoFOShiVvci+pH3ASMJ/oyCD5/guAB3Y1eOecS4eqZTXMXriBi4+fwMCy/LuSuK0OjwjM7EngSUm3mNmyXVj2CODW0E5QBNxlZg9Jeh2YIekHRP0Y3bwrgTvnXLr84vE3GDKgDxccOy7TofSKOH0N/UHS2WZWCyBpD2CGmZ3c2ZvM7BXgiHbGLwaO2pVgnXMu3Z5fuomnF1XzrVMPYEDfOLvI3Bfn9NGhySQAYGY1wF7pC8k55zLnusfeYGh5X847ZlymQ+k1cRJBQtLODrgljSW6HsA55/LKf97ayL/f2sgXp+1Lvz75cyvKrsQ57vk28IykJ8Pw8YTz+51zLl+YGdc9/gZ7DezLJ4/Ov5vPdCZON9QPS5pCdOqngK+aWXXaI3POuV7077c2MmfJJq46bXJe3Zg+jjjdUAs4BZhiZn8D+kvyxl7nXN4wM37+2BuMGFyWt/cc6EycNoLfAscCHw/DdUTdRTjnXF54alE1Vctq+NL0iQV3NADx2giONrMpkl6E6Kwh7zraOZcvkkcDoyr6cc7UwjsagHhHBE3hojCD6IphIJHWqJxzrpfMWriel1fU8uUTJtKnJM4uMf/EKfUvgfuAvST9EHgG+FFao3LOuV6QPBoYPaQfH60s3Dvwxjlr6HZJVcCJRGcNnWFm89MemXPOpdljr6/jtVVbuPasQyktLsyjAYiRCCRdD9xpZt5A7JzLG4mEcd3jixi3Z3/OPCJ/b0MZR5wU+ALwHUlvSrpWUtc3OXDOuSz36Otrmb9mC5eeNImSAj4agBiJwMxuNbNTiTqKewO4RtKitEfmnHNpkkgY1z22iAnDBnDaYYV9NADxjgiSJgIHAOOABWmJxjnnesE/XlvDwnV1/NdJ+2XkFpTZJs6VxckjgO8BrwGVZvbhtEfmnHNp0JIwfvH4IvYbXs4HDxmR6XCyQpwLypYAx3r/Qs65fHDX3BW8ub6e33xiih8NBHFOH72hNwJxzrl0W7N5Oz/6+3yOnbAnHzh470yHkzUKu6ncOVcwzIxv3vsqzQnjmo8eSpEfDezkicA5VxDurlrJ7IUbuPyU/RmzZ/9Mh5NVYiUCSe+WdFF4PUzS+PSG5ZxzPWfdlh18/6HXOWrcEM4vkBvSd0ecs4a+C1wOfDOMKgVuS2dQzjnXU8yMb9/3Kg3NCa45y6uE2hPniOBM4DRgK4CZrQYGdvUmSaMlzZI0X9I8SZeG8UMkPSZpUXjeY3cK4JxznXngpdU8Pn893zh5f8YPHZDpcLJSnETQaGbG291Qx/0km4GvmdmBRLe5/JKkg4ArgJlmNgmYGYadc67Hra/bwZV/m8eUMRVcdJzXaHckTiK4S9LvgQpJnwMeB27q6k1mtsbMXgiv64D5wCjgdODWMNutwBm7ErhzznXGzPjf++exrbGFn5x1mF8z0AlFf/a7mEl6H/B+om6oHzGzx7q1Emkc8BRwMLDczCpSptWY2TuqhyRdDFwMMHz48MoZM2a0ml5fX095eXl3wshq+VYeyL8yeXmyX2qZ5qxp5rcvN3DOfqWcOiE3b6q4u9to+vTpVWbWdUehZpbWB1AOVAEfCcOsccyhAAATjElEQVS1babXdLWMyspKa2vWrFnvGJfL8q08ZvlXJi9P9kuWqbpuhx3xvUfttF89bU3NLZkNajfs7jYC5lqM/XScs4bqJG0Jjx2SWiRtiZONJJUC9wC3m9m9YfQ6SSPC9BHA+jjLcs65uL774DzqdzRz7dmHFXwX03HE6YZ6oJkNCo8y4KPAr7t6nyQBNwPzzeznKZMeBC4Iry8AHuh+2M45176HX1vDQ6+s4SsnTmS/4V2e4OjYhSuLzex+4IQYsx4HnAecIOml8DgVuBp4X+jR9H1h2Dnndlt9o/Gd+19j8shBfP69+2Y6nJwR51aVH0kZLAKmEk4l7YyZPUPUuNyeE2NF55xz3XD7/AZqtyX486ePLuh7EHdXnG6oU+890AwsJToF1Dnnssbjr6/jP2tauPTESRw0clCmw8kpcbqhvqg3AnHOuV21fOM2/vueVxg9sIgvTZ+Y6XByToeJQNKv6KQKyMy+kpaInHOuGzZva+LCW+bQkjC+NKUvfUq8Sqi7OvvE5hKd/18GTAEWhcfhQEv6Q3POuc41Nif4wm1VrNi0jRvPq2TvAZ4EdkWHRwRmdiuApAuB6WbWFIZvAB7tleicc64DZsa37nuV/yzeyHXnHsbRE/Zk9vJMR5Wb4qTPkbTubbQ8jHPOuYz5zaw3ubtqJf910iTOPGKfTIeT0+KcNXQ18KKkWWH4vcCVaYvIOee68MBLq/jpo29w5hGjuPTESZkOJ+fFOWvoT5L+CRwdRl1hZmvTG5ZzzrXv+aWb+MZfX+Go8UO4+qOHEHVi4HZHh1VDkg4Iz1OIqoJWhMfIMM4553rV0uqtXPznueyzRz9uPK+SviXFmQ4pL3R2RHAZUTfQP2tnmhGvmwnnnOsRNVsbueiW5wH444VHUtE/N7uWzkadnTV0cXie3nvhOOfcOzU0t/D526pYVbOd2z93NOP8lpM9Kk431GdLGhhef0fSvZKOSH9ozjkXnSb6zXteZc6STVx79qEcOW5IpkPKO3FOH/0fM6uT9G7gZKLbS96Q3rCccy7yy5lvcu+Lq/ja+/bj9MNHZTqcvBQnESSvIv4g8DszewDwyjnnXNrd/MwSrnv8DT46ZR8uOcH7EEqXONcRrAo3rz8JuEZSX3bhPgbOORdXImFc8/ACfv/UYk6ePJwff8RPE02nODv0c4BHgFPMrBYYAnwjrVE55wpWY3OCy+56id8/tZjzjhnLbz9Z6R3JpVmcC8q2SVoPvJuo07nm8Oyccz2qvqGZL95WxdOLqvn6+/fjS9Mn+pFAL4hzh7LvEt2VbH/gT0ApcBvRrSidc65HbKhr4KJb5jB/TR0/+eihnHPk6EyHVDDitBGcCRwBvABgZquTp5M651xPWFK9lQv+OIf1dTu46fxKTjhgeKZDKihxEkGjmZkkA5DkV3I453rMyytq+fQtz5Mw447PHcMRY/bIdEgFJ04LzF3hrKEKSZ8DHgf+kN6wnHOFYPbC9Xzsxmfp16eYe774Lk8CGRKnsfinkt4HbCFqJ/hfM3ss7ZE55/La3VUrueKeV9hv+EBu+fSR7DWwLNMhFaw4VUOEHf9jAJKKJX3SzG7v7D2S/gh8CFhvZgeHcUOAO4FxwFLgHDOr2eXonXM5pyVh/PqJN7nu8Tc4buKe3PCpSgaWlWY6rILWWTfUgyR9U9KvJb1fkUuAxUTXFnTlFuCUNuOuAGaa2SRgZhh2zhWI5Ru3ce7v/8N1j0c3lfnThUd5EsgCnR0R/B9QA/wH+CzRRWR9gNPN7KWuFmxmT0ka12b06cC08PpWYDZweXcCds7lHjNjxvMr+P5Dr1Ms8fNzDuPMI0b5NQJZQmbW/gTpVTM7JLwuBqqBMWZWF3vhUSJ4KKVqqNbMKlKm15hZu61Dki4muh8Cw4cPr5wxY0ar6fX19ZSXl8cNJevlW3kg/8rk5dk1tQ0J/vRaIy9vaOHAIUV89pC+7NkvPVcK+zZqbfr06VVmNrXLGc2s3QfwQmfDcR5EbQGvpQzXtpleE2c5lZWV1tasWbPeMS6X5Vt5zPKvTF6e7vvnq6vt8Ksesf2+/Q+7+enF1tKSSOv6fBu1Bsy1GPvYzqqGDpO0JbwW0C8MK8ofNmgXEtQ6SSPMbI2kEcD6XViGcy7LbdnRxJUPzuPeF1ZxyKjB/Pycw5g03K9DzVad3aEsHTcDfRC4ALg6PD+QhnU45zLo329W8/W/vsy6uga+csJEvnziJEqLvdO4bBbr9NFdIekOoobhoZJWAt8lSgB3SfoMsBw4O13rd871rq0Nzfzs0Tf447+WMH7oAO7+wrF+gViOSFsiMLOPdzDpxHSt0znX+xqbE9wxZzm/emIR1fWNnH/sWK74wAH075O23YvrYb6lnHO7JJEwHnx5NT97bCErNm3n6PFDuPH8A5jiRwE5xxOBc65bzIzZCzdwzcMLWLC2joNGDOKWiw7mvfsN8+sCcpQnAudcbFXLNnHNPxcyZ+kmxgzpz/UfO5wPHzqSoiJPALnME4FzrksL19Zx7SMLeXz+OoaW9+X7p0/m3CPH+C0k84QnAudcu8yMqmU1/Pk/y/jbK6sp71PCN07en4uOG+cNwXnGt6ZzrpX6hmbue3EVtz+7jAVr6xjYt4TPvWcCX3zvvuwxoE+mw3Np4InAOQfAgrVbuO3ZZdz3wiq2NrZw0IhB/Pgjh3DaYSMZ0Nd3FfnMt65zBayhuYWHX1vLbc8u4/mlNfQpKeJDh47gU8eM5YjRFX4WUIHwROBcgTEz5q+p468LG/na00+wcWsjY/fsz7dOPYCzK0d79U8B8kTgXAFIJIwXltfwyLy1PDJvHcs3bUPASQcN57xjxvLuiUP9FNAC5onAuTzV2JzgP4s38si8tTz2+jo21DVQWizete9QvjhtX/rXvMXpJ3fdVb3Lf54InMsjWxuaeeqNDTwyby0zF6ynbkcz/fsUM33/vXj/5OFMP2AvBoVbQ86evTjD0bps4YnAuRxWt6OJuctqeG7xJp5bspFXV26mOWFU9C/llMl7c/LkvXn3pKGUlaajV3mXLzwROJdDNm9v4vkl0U7/uSWbeG3VZhIGpcXi0H0quPj4Cbx70lCOGjeEEr8HgIvJE4FzWaq5JcFbG7Yyb/VmXlm5mTlLNjF/7RbMoE9JEYePruCS6RM5esKeTBmzB/36+L9+t2s8ETiXBbY1NrNgbR3zVm/h9dWbeX31FhasraOhOQFAWWkRU8bswX+duB9HTxjC4aMrvLrH9RhPBM71oh1NLSzduJUlG7ayuHorC9fWMW/1ZpZUbyVh0TyD+5UyeeQgzj92LAeNHMTkkYOZMHSAV/W4tPFE4FwPa25JsLJmO0uqo539kup6llRvZWn1NlbVbm8178jBZRw0cjAfOnRk2OkPYlRFP7+i1/UqTwTOdYOZUbutidWbt7O6dgdrwvPq2u07X6/bsoPm5N97YGBZCROGlXPU+CGMHzqAcUMHMCE8l3sfPi4L+LfQOaIqm41bG6mua6C6PvlofPu5roEl67ZRO/NhdjQlWr23tFjsPbiMkYP7cdT4IYysKGPcngMYPzR6DBnQx//hu6zmicDlDTOjoTlBfUMzm7c3RY9tTdRub6R2WxO125p2jq/d1kjt9mhcdV0DdQ3N7S6zvG8JQ8v7sGd5X/YpL+KDk8YwsqIfIweXMaKiHyMryhg6oK93z+ByWkYSgaRTgOuBYuAPZnZ1JuJwmWFmNLUYDc0t7GhKsKOphe1NLWxrbGFbYzPbG98e3t6YfG5mW2MLWxubqdvRTH1DM/XJ54a3h1OrZNozqKyEwf1LqejXh4r+pYyq6MfQ8r4MLe/D0PK+7JnyetjAvq3OzJk9ezbTph2U7o/HuV7X64lAUjHwG+B9wErgeUkPmtnrPb2uG558i7fW1zO4XykV/UsZ3K+Uwf37RMP9SneOH1hWSnEO/aNLJIwWM1oSRiL5nKDVuOT4loTRnHxuMZoTiVbD0fTEzmlNLUZTSzTcFMY3tUTvaW5J0NhiNDYnaGpJtHpuaEnQ1JygseXtcdWbtvOjF5+koTlBQ1OCHc0tNDQlaGhuoYv99TtI0K+0mPK+JZSXlUTPfUsYM6A/5WUlDOxbwoAwbWDfEgaWvb3NK/r3oaJfKYP65dZ2dq63ZOKI4CjgTTNbDCBpBnA60OOJYGn1Vp55s5rN25vY1tjS4XxSVAVQUiSKi4QkigTFCq+LotdFEhKt6nvNWu/RrM2AhXkSBoZhRngYBiQsGrejoZHSpx/buQM3I+zQox18chkt3d2DpkGf4iL6lBRRWqzwHA2/PT56PbCPGDm0nL6lRZSVFNO3tIi+JUWUlRa3eu5bWkz/PtGjX2nJztdlO8eXUFZa5PXszqWJ2u7I0r5C6SzgFDP7bBg+DzjazC5pM9/FwMUAw4cPr5wxY0ar5dTX11NeXh57vU0JY2uTsa0JtjZZyiMa3tZstOzcSUOC8GyQSN2BE41L3Sd1tnsS0bwKA0Vh7uS45HJampvoW1q6c3yRkvNESakoZf7kcNHO+d6eJzk+OV9xSGhRMnt7XJHUari4KPksStoMF4tW4+LukLu7jbKdlyf75VuZdrc806dPrzKzLruYzcQRQXt7kXdkIzO7EbgRYOrUqTZt2rRW06P62mlt35az8q08kH9l8vJkv3wrU2+VJxOXKq4ERqcM7wOszkAczjnnyEwieB6YJGm8pD7Ax4AHMxCHc845MlA1ZGbNki4BHiE6ffSPZjavt+NwzjkXych1BGb2D+AfmVi3c8651rw7Q+ecK3CeCJxzrsB5InDOuQLnicA55wpcr19ZvCskbQCWtRk9FKjOQDjpkm/lgfwrk5cn++VbmXa3PGPNbFhXM+VEImiPpLlxLp3OFflWHsi/Mnl5sl++lam3yuNVQ845V+A8ETjnXIHL5URwY6YD6GH5Vh7IvzJ5ebJfvpWpV8qTs20EzjnnekYuHxE455zrAZ4InHOuwOVcIpB0iqSFkt6UdEWm4+kJkpZKelXSS5LmZjqe7pL0R0nrJb2WMm6IpMckLQrPe2Qyxu7qoExXSloVttNLkk7NZIzdIWm0pFmS5kuaJ+nSMD4nt1Mn5cnlbVQmaY6kl0OZrgrjx0t6LmyjO0P3/T277lxqIwg3vn+DlBvfAx9Px43ve5OkpcBUM8vJC2EkHQ/UA382s4PDuJ8Am8zs6pCw9zCzyzMZZ3d0UKYrgXoz+2kmY9sVkkYAI8zsBUkDgSrgDOBCcnA7dVKec8jdbSRggJnVSyoFngEuBS4D7jWzGZJuAF42s9/15Lpz7Yhg543vzawRSN743mWQmT0FbGoz+nTg1vD6VqIfac7ooEw5y8zWmNkL4XUdMB8YRY5up07Kk7MsUh8GS8PDgBOAu8P4tGyjXEsEo4AVKcMryfGNHxjwqKQqSRdnOpgeMtzM1kD0owX2ynA8PeUSSa+EqqOcqEZpS9I44AjgOfJgO7UpD+TwNpJULOklYD3wGPAWUGtmzWGWtOzzci0RxLrxfQ46zsymAB8AvhSqJVz2+R2wL3A4sAb4WWbD6T5J5cA9wH+Z2ZZMx7O72ilPTm8jM2sxs8OJ7uV+FHBge7P19HpzLRHk5Y3vzWx1eF4P3Ef0Bch160I9brI+d32G49ltZrYu/FATwE3k2HYK9c73ALeb2b1hdM5up/bKk+vbKMnMaoHZwDFAhaTk3STTss/LtUSQdze+lzQgNHYhaQDwfuC1zt+VEx4ELgivLwAeyGAsPSK5wwzOJIe2U2iIvBmYb2Y/T5mUk9upo/Lk+DYaJqkivO4HnETU9jELOCvMlpZtlFNnDQGE08F+wds3vv9hhkPaLZImEB0FQHQP6b/kWpkk3QFMI+oydx3wXeB+4C5gDLAcONvMcqbxtYMyTSOqcjBgKfD5ZP16tpP0buBp4FUgEUZ/i6hePee2Uyfl+Ti5u40OJWoMLib6k36XmX0v7CNmAEOAF4FPmVlDj6471xKBc865npVrVUPOOed6mCcC55wrcJ4InHOuwHkicM65AueJwDnnCpwnAldQJNWH53GSPtHDy/5Wm+F/9+TynUsXTwSuUI0DupUIQu+3nWmVCMzsXd2MybmM8ETgCtXVwHtCn/VfDZ19XSvp+dBh2ecBJE0L/d7/hejiJSTdHzoInJfsJFDS1UC/sLzbw7jk0YfCsl9TdN+Jc1OWPVvS3ZIWSLo9XDHrXK8q6XoW5/LSFcDXzexDAGGHvtnMjpTUF/iXpEfDvEcBB5vZkjD8aTPbFLoBeF7SPWZ2haRLQodhbX2E6GrXw4iuVH5e0lNh2hHAZKL+Y/4FHEfUD71zvcaPCJyLvB84P3QB/BywJzApTJuTkgQAviLpZeBZok4QJ9G5dwN3hM7Q1gFPAkemLHtl6CTtJaIqK+d6lR8ROBcR8GUze6TVSGkasLXN8EnAsWa2TdJsoCzGsjuS2mdMC/6bdBngRwSuUNUBA1OGHwG+GLo2RtJ+oTfYtgYDNSEJHEDUTXBSU/L9bTwFnBvaIYYBxwNzeqQUzvUA//fhCtUrQHOo4rkFuJ6oWuaF0GC7gfZvCfgw8AVJrwALiaqHkm4EXpH0gpl9MmX8fcCxwMtEvWL+t5mtDYnEuYzz3kedc67AedWQc84VOE8EzjlX4DwROOdcgfNE4JxzBc4TgXPOFThPBM45V+A8ETjnXIH7/6smSlyKgPjJAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "init_c_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "i = 0\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "next_norm_values = []\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    next_norm_values.append(residue_vector_norm)\n",
    "    i += 1\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "# Printing details about this method time, iteration count etc.\n",
    "print_details(A, B, x, i, toc())\n",
    "\n",
    "# Plotting the value of residue vector norm\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(range(1, i + 1), next_norm_values)\n",
    "\n",
    "ax.set(xlabel='Iteration', ylabel='Residue vector second norm',\n",
    "       title='Residue vector second norm values in next iterations')\n",
    "ax.grid()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The plot clearly shows that changed input matrix makes Jacobis method not suitable for solving \n",
    "this system. Residue vector increases each iteration, which means that the approximated\n",
    "solution vector won't converge.  \n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false
    }
   },
   "source": [
    "**Testing Gauss-Seidl's on data from C exercise**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "..............................\n",
      "Residue vector second norm value:  1536159.86028818\n",
      "Minutes:  0 seconds:  14\n",
      "Iteration count:  30\n",
      "Time per iteration in seconds:  0.4666666666666667\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEWCAYAAADPZygPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XmcHFW5//HPNzPZ95AwhCQQliD7EiKrywgCAZXggoIiEdGoF5cr1wVcfqDCFZd7EQTBCJGgyCKI5CoSAtIgsoVAwpKAxAQyQ3aSmckkmSwzz++POk1qOt0znZnprume5/169WuqTp2qc05Xdz9TVadOycxwzjnnktAr6Qo455zruTwIOeecS4wHIeecc4nxIOSccy4xHoScc84lxoOQc865xHgQ6qYkvSypOseyakm1Ra6Sy5OklKTPJV2PQpL0GUmPF7nMv0maWswyO0vSXpIaJVUkWIdPSXowqfLb40GokyS9Lmlz+KCtlHSLpEGd3a6ZHWJmqS6oYuIkmaT9k66HK21mdrqZzUy6HnHt/cNhZsvMbJCZNeeTvwvqMz583ypjdbjNzE4tVJmd5UGoa3zIzAYBRwJHAZcmXJ+yEf8yudaS/O/aJaMc97kHoS5kZiuB2UTBCABJfSX9XNIySask3Sipf1g2UtJfJNVJWifpH5J6hWWvS3p/mO4fjrDWS1oIvDNebuaRRsh7RWz+g5Lmh3KekHR4tvqHuv08I+0+SReH6T0l3SNpjaSlkr4ay1ch6TuS/i1pg6R5ksZJeixkWRCOFj8R8n9e0uLQ7lmS9sxoz0WSXgNey1LPfpJ+L+mt0Ka5kqrCsqGSbpa0QtKbkq6If3FDuYtCHRdKmhjSDwr/pdaFU6FnZryf10v6a1jvaUn7xZafIukVSfWSrgOU7f0NeS+XdJekW8O2XpY0Kba8vXrcIOl+SRuB94W0X4VTVY2S/ilpD0m/CJ+XVyQd1cH9fUlsfy6U9OEc29npv29l/Mcv6bPhfV8vabakvUO6JF0taXV4/16QdGiOct7epsLpQEXfrfXh83h6G+/765K+EbZfL+lOSf1iy7N+RyTtFz6j6c/JnpLWKjolfiXwbuC68N5f19Z7kyu/pAMlzQnlvCrp47H1s+3zD0h6XlKDpBpJl8eKTH/f6kIZxyvj1KmkExR9Z+rD3xMy3uMfhc/RBkkPShoZluX83nWKmfmrEy/gdeD9YXos8CJwTWz5L4BZwAhgMPB/wI/Dsh8DNwK9w+vdgLJs9yrgH2Eb44CXgNpYGQbsH5u/BbgiTE8EVgPHAhXA1LDtvlna8h6gJlaH4cBmYE+if1jmAf8P6APsCywBTgt5vxna/g6iH+EjgN1y1O8kYG2oW1/gl8BjGe2ZE9rbP0s9vxDexwGhTUcDQ8KyPwO/BgYCuwPPAF8Iy84G3iQK4gL2B/YO7/1i4DuhbScBG4B3xN7PdcAxQCVwG3BHWDYSaAA+FrbzdWA78Lkcn5fLgSbgjFD3HwNPhWX51KMeODHsj34hbW14D/oBfweWAueH7V8BPJKjLjn3d+z9Su/7TwAbgdFh2WeAx8P0+LDPKmPbTqXfA+Cs0K6Dwvv3PeCJsOw0os/VsLBPDkqXkaW+8W1+BtgGfD6080vA8nRbcnxPnwntGQEsAr6Yz3cklLGI6PM2G/h5tjrlKLfVe5OZn+hzWgNcEN6biWF/HtLGPq8GDgvzhwOrgLPa2BfxfTUCWA98OpR3bpjfLVa/fwMHAP3D/FXtfe869Rua1I93ubzCh7WR6MfCgIeBYWGZiL64+8XyHw8sDdM/BO4j9gOdsd10EFoCTI4tm0b+QegG4EcZ234VeG+WMgUsA94T5j8P/D1MHwssy8h/KfDb2Dan5HiPMut3M/DT2Pwgoh+U8bH8J7Xxnn8WeAI4PCO9CthCLHCFL9kjYXo28LUs23s3sBLoFUu7Hbg89n7eFFt2BvBKmD6fEERi72EtbQehh2LzBwObd6Eet2Zs7xbgN7H5rwCLYvOHAXU56pJzf+fIPz+9j9m1IPQ34MLYsl7AJqJ/AE4C/gUcF293jvLj2/wMsDi2bECowx5tfE/Pi83/FLgx3+8I0T+SLwIvEPsHjs4HoU8A/8hY59fAZbn2eZYyfgFc3ca+iO+rTwPPZKz/JPCZWP2+F1v2H8ADbX3vOvvy03Fd4ywzG0z0H8qBRP8dA4wi+nLMC4evdcADIR3gZ0T/IT4oaYmkS3Jsf0+i/5bS3tiFuu0N/Fe6/FCHcWGbrVj0SbuD6Icb4JNE//Wnt7Nnxna+Q/TDT9jmv/Os057xNphZI/AWMCaWpyZzpZjfEQWUOyQtl/RTSb3ZcVSzIlbHXxMdEbVVxz2BGjNriaW9kVGflbHpTUSB8+11Y22xduqebVv9wqmsfOqRbdurYtObs8xn7SjTzv5G0vmxU1R1wKHs+Gzvir2Ba2LbWUcUAMeY2d+B64DrgVWSpksakud2334fzWxTmGyrU1CufZjPd+Q3RO3/pZltybN++dgbODaj7E8Be8TytNrnko6V9Iii0+L1wBfJf7+0+u4F+X7Wc33vOsWDUBcys0eJ/nNJn2dfS/QjcIiZDQuvoRZ1YsDMNpjZf5nZvsCHgIslnZxl0yuIvhRpe2Us30QU7NIyP8BXxsofZmYDzOz2HM24HfhYOGd/LHBPbDtLM7Yz2MzOiC3fL8v2sllO9OUDQNJAYDeiU2VplmtlM9tmZj8ws4OBE4APEh2R1BAdCY2M1XGImR3STh2XA+MUrscFe2XUJ5dW+0aSaL2vdkU+9cj5vnRQ1v0d5n8DfJnoVM0wotPA2a53bQx/2/oMfiHjs9PfzJ4AMLNrzexo4BCi00Df7Lrm5aXN74ii3q6/IDqCv1zSiNi6u7o/MvPXAI9mlD3IzL7Uxjp/IDoyG2dmQ4lO6StH3kytvntBXp/1Nr53neJBqOv9AjhF0pHhP9rfAFdL2h1A0hhJp4XpD0raP/xwNQDN4ZXpLuBSScMljSU65RI3H/ikos4Bk4H3xpb9Bvhi+O9JkgaGC5uDs1XezJ4H1gA3AbPNrC4segZokPRtRR0lKiQdKindSeIm4EeSJoRyDpe0W1i2iugaUtofgAskHSmpL/DfwNNm9nquNzVO0vskHaaow0ED0am8ZjNbATwI/I+kIZJ6KbqwnH4/bgK+IenoUMf9w4/t00Q/pN+S1FvR/VkfIjpKaM9fgUMkfSQczXyV1j/Au6Iz9eiQNvb3QKIftDUAki4gOhLIto01RD9i54XPxWdpHexvJPr8HhK2NVTS2WH6neGz2Zuo7U1k/w4UUnvfkWuAeWb2OaL9fWNs3czPdnsy8/8FOEDSp8M+7x3ek4Pa2MZgYJ2ZNUk6hugINm0N0NJGne4P5X1SUWeJTxCdEv5LexXP9b1rb732eBDqYuELeSvw/ZD0baJTbk9JagAeIrp4DzAhzDcSnZf9lWW/N+gHRIfMS4l+ZH+XsfxrRD9W6UP5P8fq8yzRuf7riC5ALiY6R9yW24H3EwWL9HaaQxlHhnqsJfrhGhqy/C9RsHyQ6AN6M9GFTYiug8wMpxs+bmYPE70/9xAdSewHnNNOneL2AO4O5SwCHgV+H5adT3RRf2Fo793A6NCGPwJXhnZtIHqfRpjZVuBM4PTQrl8B55vZK+1VxMzWEl3Av4rolOIE4J+70Jb4tjpcj07Ktr8XAv9D9LlcRXRtqa12fZ7oCOYtoiOaJ2Lbuhf4CdFpnAaiI6p0T7YhREFgPdFn/C12nEkoira+I5KmAJOJTnkBXAxMlPSpMH8N0ZHkeknX5lFcq/xmtgE4lejzv5zoVNhPiDrs5PIfwA8lbSDqKHRXrC2biD7j/wzft+My2voW0RHMfxG9198CPhg+x+1p63vXYeleMc4551zR+ZGQc865xHgQcs45lxgPQs455xJTsCAkaYaioTheykj/iqKhKV6W9NNY+qWKhnF5Nd17LKRPDmmLFbuPRtI+ioZPeU3REBx9QnrfML84LB/fXhnOOeeSUbCOCZLeQ9Tr61YzOzSkvQ/4LvABM9siaXczWy3pYKIeOscQ3Uz1ENH9AhDdTX0K0V3oc4FzzWyhpLuAP5nZHZJuBBaY2Q2S/oPojt4vSjoH+LCZfSJXGaHXV04jR4608ePHt0rbuHEjAwcO7OQ71H2UW3ug/NpUbu2B8mtTubUHOtemefPmrTWzUe1m7MrhFzJfRENIvBSbv4swFE1GvkuBS2Pzs4mGtzme6N6FVvmIbsxay46hMN7Ol143TFeGfMpVRnttOProoy3TI488slNaKSu39piVX5vKrT1m5demcmuPWefaBDxrecSJYg+TfwDwbkWjyTYB3zCzuURDRjwVy1fLjmEkajLSjyW6u77OzLZnyT8mvY6ZbVc0rMVu7ZTRiqRpROOzUVVVRSqVarW8sbFxp7RSVm7tgfJrU7m1B8qvTeXWHihOm4odhCqJRuo9jmgk47sk7Uv2oUCM7NesrI38tLGsrXVaJ5pNB6YDTJo0yaqrq1stT6VSZKaVsnJrD5Rfm8qtPVB+bSq39kBx2lTs3nG1RNdxzMyeIRpeYmRIj4+3NZbo7uFc6WuBYdrx/JJ0OvF1wvKhRAMm5tqWc865hBQ7CP2ZaOh2JB1ANLzKWqLB+M4JPdv2IRr65BmijggTQk+4PkRDW8wK5xsfIXqGC0TP/7gvTM8K84Tlfw/5c5XhnHMuIQU7HSfpdqJHG4yUVAtcBswAZoRu21uBqSFAvBx6uy0keiDYRbbjmexfJupEUAHMMLOXQxHfJhqL6grgeaKxygh/fydpMdER0DkAZpazDOecc8koWBAys3NzLDovR/4riQbey0y/n2jk18z0JUTdrTPTm4gGlMy7DOecc8nwEROcc84lxoOQc865nVzz0Gu8tLbwVyw8CDnnnGtle3ML1zz8L/613oOQc865IlvTuIUWg+F9s91e2bU8CDnnnGtlZX0TAMP7eRByzjlXZB6EnHPOJWZlQzoIFT5EeBByzjnXysr6JvpU9GJw78KX5UHIOedcKysbmqga2hfJT8c555wrshX1TYwe0r8oZXkQcs4518qqhiaqhvYrSlkehJxzzr3NzKIjIQ9Czjnniq1u0za2bm+haogHIeecc0WW7p7tR0LOOeeKLn2jqh8JOeecK7r0kdAefiTknHOu2FbUNyHB7oP7FqW8ggUhSTMkrQ6P8s5c9g1JJmlkmJekayUtlvSCpImxvFMlvRZeU2PpR0t6MaxzrcJdVZJGSJoT8s+RNLy9MpxzzkVW1TcxclBfelcU5xilkKXcAkzOTJQ0DjgFWBZLPh2YEF7TgBtC3hHAZcCxRI/yviwdVEKeabH10mVdAjxsZhOAh8N8zjKcc87tsKKheN2zoYBByMweA9ZlWXQ18C3AYmlTgFst8hQwTNJo4DRgjpmtM7P1wBxgclg2xMyeNDMDbgXOim1rZpiemZGerQznnHPBqvqmonVKAKgsWkmApDOBN81sQcaYRGOAmth8bUhrK702SzpAlZmtADCzFZJ2b6eMFVnqOY3oaImqqipSqVSr5Y2NjTullbJyaw+UX5vKrT1Qfm0ql/bUvLWRsX03k0qlitKmogUhSQOA7wKnZlucJc06kN5mFfJdx8ymA9MBJk2aZNXV1a2Wp1IpMtNKWbm1B8qvTeXWHii/NpVDezZt3c6mB2Yz8aD9qK7evyhtKmbvuP2AfYAFkl4HxgLPSdqD6KhkXCzvWGB5O+ljs6QDrEqfZgt/V4f0XNtyzjnHjnuEyuKaUCYze9HMdjez8WY2nigoTDSzlcAs4PzQg+04oD6cUpsNnCppeOiQcCowOyzbIOm40CvufOC+UNQsIN2LbmpGerYynHPOEbtHqByuCUm6HagGRkqqBS4zs5tzZL8fOANYDGwCLgAws3WSfgTMDfl+aGbpzg5fIuqB1x/4W3gBXAXcJelCoh54Z7dVhnPOuciqIt+oCgUMQmZ2bjvLx8emDbgoR74ZwIws6c8Ch2ZJfws4OUt6zjKcc85FN6pCcYOQj5jgnHMOiLpnD+5XyYA+xes47UHIOeccQFGfI5TmQcg55xwQnqhaxE4J4EHIOedc4EdCzjnnErGtuYU1jVuK2j0bPAg555wD1mzYghnsMbR/UctttwuEpElEw+3sHfKLqMfz4QWum3POuSLZ8TC74jxHKC2ffni3Ad8EXgRaClsd55xzSUgP2bPHkG52JASsMbNZBa+Jc865xKxM4EZVyC8IXSbpJqIHxG1JJ5rZnwpWK+ecc0W1qqGJPpW9GD6gd1HLzScIXQAcCPRmx+k4AzwIOedcmVhR38QeQ/qR8ay3gssnCB1hZocVvCbOOecSs7KhqejdsyG/LtpPSTq44DVxzjmXmJX1TUW/HgT5HQm9C5gqaSnRNSHvou2cc2XEzKIjoW4ahCYXvBbOOecSs37TNrZub0nkdFybQUhSL+CvZrbTc3ucc86Vh6S6Z0M714TMrAVYIGmvXd2wpBmSVkt6KZb2M0mvSHpB0r2ShsWWXSppsaRXJZ0WS58c0hZLuiSWvo+kpyW9JulOSX1Cet8wvzgsH99eGc4515OtbNgMdMMgFIwGXpb0sKRZ6Vce693Czqfy5gCHhutJ/wIuBQgdH84BDgnr/EpShaQK4HrgdOBg4NxYJ4mfAFeb2QRgPXBhSL8QWG9m+wNXh3w5y8ijHc45V9ZW1ke3gHa703HBDzqyYTN7LH4UEtIejM0+BXwsTE8B7jCzLcBSSYuBY8KyxWa2BEDSHcAUSYuAk4BPhjwzgcuBG8K2Lg/pdwPXKer4nquMJzvSPuecKxcrG5roJRg1uLjjxkEeR0Jm9ijwCjA4vBaFtM76LPC3MD0GqIktqw1pudJ3A+rMbHtGeqttheX1IX+ubTnnXI+2sn4zIwf1pXdF8R+skM8o2h8HfgakiLpn/1LSN83s7o4WKum7wHaiwVEJ281kZA+S1kb+trbV1jqZ9ZsGTAOoqqoilUq1Wt7Y2LhTWikrt/ZA+bWp3NoD5demUm7PwtebGChL5Lcun9Nx3wXeaWarASSNAh4iOtW1yyRNBT4InGxm6SBQC4yLZRsLLA/T2dLXAsMkVYajnXj+9LZqJVUCQ4F17ZTRiplNB6YDTJo0yaqrq1stT6VSZKaVsnJrD5Rfm8qtPVB+bSrl9vz3848yYexAqqsntUovRpvyOfbqlQ5AwVt5rrcTSZOBbwNnmtmm2KJZwDmhZ9s+wATgGWAuMCH0hOtD1LFgVghej7DjmtJU4L7YtqaG6Y8Bfw/5c5XhnHM9WlKjJUB+R0IPSJoN3B7mPwHc395Kkm4HqoGRkmqBy4h6w/UF5oRB8p4ysy+a2cuS7gIWEp2mu8jMmsN2vgzMBiqAGWb2ciji28Adkq4AngduDuk3A78LHQ/WEQUu2irDOed6qk1bt9PQtL37BiEz+6akjwInEl1XmW5m9+ax3rlZkm/OkpbOfyVwZZb0+8kS9EKPuWOypDcBZ+9KGc4511PteJhdNw1CAGZ2D3BPgevinHOuyJIcLQHyuLYj6SNhVIJ6SQ2SNkhqKEblnHPOFdbKhu5/JPRT4ENmtqjQlXHOOVdcK7r7kRCwygOQc86Vp1UNTQzpV8mAPnldnely+ZT6rKQ7gT8TPU8IADPzx3s751yJW1nfxOih/RMrP58gNATYBJwaSzPAg5BzzpW4lQ1NVCV0Kg7y66J9QTEq4pxzrvhW1jdx4B6DEyu/+KPVOeec6xa2NbewpnELeyR4Os6DkHPO9VBrNmzBLLnu2eBByDnneqx09+zR3fGakKSL21rRzP6366vjnHOuWFaFG1WrEjwSaqtjQvpK1TuAdxKNQg3wIeCxQlbKOedc4XXrIyEz+wGApAeBiWa2IcxfDvyxKLVzzjlXMKsamuhT2YthA3onVod8rgntBWyNzW8FxhekNs4554pmRX0To4f2IzxaJxH53Kz6O+AZSfcS3aT6YeDWgtbKOedcwa2qb0r0ehDkd7PqlZIeAN4Vki4ws+cLWy3nnHOFtrKhiaP2GpZoHfIdsW4+sCKdX9JeZrasYLVyzjlXUGbGyoamRO8RgvyeJ/QVYBUwB/gL8Nfwt731ZkhaLemlWNoISXPC84nmSBoe0iXpWkmLJb0gaWJsnakh/2uSpsbSj5b0YljnWoWTmh0pwznnepr1m7axdXtL4qfj8umY8DXgHWZ2iJkdbmaHmdnheax3CzA5I+0S4GEzmwA8HOYBTgcmhNc04AaIAgpwGXAs0aO8L0sHlZBnWmy9yR0pwznneqIV9ZuBZLtnQ35BqAao39UNm9ljwLqM5CnAzDA9Ezgrln6rRZ4ChkkaDZwGzDGzdWa2nuhobHJYNsTMnjQzI+oocVYHy3DOuR7n7RtVEw5C+VwTWgKkJP2V1s8T6siICVVmtiKsv0LS7iF9DFGwS6sNaW2l12ZJ70gZKzIrKWka0dESVVVVpFKpVssbGxt3Sitl5dYeKL82lVt7oPzaVGrteXTZNgBeX/g8DUuyH48Uo035BKFl4dUnvAohWyd160B6R8rYOdFsOjAdYNKkSVZdXd1qeSqVIjOtlJVbe6D82lRu7YHya1Optee5B1+l16LFnHlKNZUV2YNQMdqUTxft9MgJg6NZa+xEeaskjQ5HKKOB1SG9FhgXyzcWWB7SqzPSUyF9bJb8HSnDOed6nBX1TYwa3DdnACqWfHrHHSrpeeAl4GVJ8yQd0sHyZgHpHm5Tgfti6eeHHmzHAfXhlNps4FRJw0OHhFOB2WHZBknHhV5x52dsa1fKcM65Hqc7dM+G/E7HTQcuNrNHACRVA78BTmhrJUm3Ex3FjJRUS9TL7SrgLkkXEp3iOztkvx84A1hM9CjxCwDMbJ2kHwFzQ74fmlm6s8OXiHrg9Qf+Fl7sahnOOdcTraxvYt9RA5OuRl5BaGA6AAGYWUpSuzU3s3NzLDo5S14DLsqxnRnAjCzpzwKHZkl/a1fLcM65nmZlQxMn7Ldb0tXIr3ecpO8TjSEHcB6wtHBVcs45V0gbt2xnQ9P2RB/rnZbPFanPAqOAP4XXSPxUlnPOlayV4R6hPYb2Tbgm+fWOWw98tQh1cc45VwSr6pN/ompaPr3j5kgaFpsfLml2YavlnHOuUHY8UbU0TseNNLO69Ew4Mtq9jfzOOee6sbdPx5XCkRDQImmv9IykvWl/dALnnHPd1Mr6Job2703/PhVJVyWv3nHfBR6X9GiYfw9hXDXnnHOlp7vcqAr5dUx4IDx75zii8de+bmZrC14z55xzBbGyvok9Eh49Oy2fjgkielbPRDP7P2CApGMKXjPnnHMF0Z2OhPK5JvQr4HggPQLCBuD6gtXIOedcwWxrbmFt45ZucySUzzWhY81sYhjEFDNbL6lQj3RwzjlXQKs3bMGMbhOE8jkS2iapgtAjTtIooKWgtXLOOVcQK+vToyWUThC6FrgX2F3SlcDjwH8XtFbOOecK4u0g1E2uCeXTO+42SfOIRqYWcJaZLSp4zZxzznW57nSjKuTXO24/YKmZXU/0YLtT4sP4OOecKx0r6zfTt7IXwwb0TroqQH6n4+4BmiXtD9wE7AP8oaC1cs45VxArG6KecdHdN8nLa9geM9sOfAS4xsy+DowubLWcc84Vwsr6zd3mVBzk3zvuXOB84C8hrVPHcZK+LullSS9Jul1SP0n7SHpa0muS7kx3A5fUN8wvDsvHx7ZzaUh/VdJpsfTJIW2xpEti6VnLcM65nmJlQ/cZLQHyC0IXEN2seqWZLZW0D/D7jhYoaQzR84kmmdmhQAVwDvAT4GozmwCsBy4Mq1wIrDez/YGrQz4kHRzWO4RoRIdfSaoI3cmvB04HDgbODXlpowznnCt7Zsaq+u5zoyrkEYTMbKGZfdXMbg/zS83sqk6WWwn0l1QJDABWACcBd4flM4GzwvSUME9YfnIYSmgKcIeZbTGzpcBi4JjwWmxmS8xsK3AHMCWsk6sM55wre+s2bmVrc0u3Oh2Xz4gJXcrM3pT0c2AZsBl4EJgH1IVrTwC1wJgwPQaoCetul1QP7BbSn4ptOr5OTUb6sWGdXGW0ImkaYaTwqqoqUqlUq+WNjY07pZWycmsPlF+byq09UH5tKoX2vNHQDMBbNf8mlXqj3fzFaFPRg5Ck4URHMfsAdcAfiU6dZUo/syhbFw5rIz3b0V1b+XdONJsOTAeYNGmSVVdXt1qeSqXITCtl5dYeKL82lVt7oPzaVArteXjRKnjiWU4+4WiO2mt4u/mL0aY2T8eFayw/6+Iy309039EaM9sG/Ak4ARgWTs8BjAWWh+laYFyoTyUwFFgXT89YJ1f62jbKcM65svf2jaqlck3IzJqBo9W1HcqXAcdJGhC2ezKwEHgE+FjIMxW4L0zPCvOE5X83Mwvp54Tec/sAE4BngLnAhNATrg9R54VZYZ1cZTjnXNlbvLqRfr17MWpQ36Sr8rZ8Tsc9D9wn6Y/AxnSimf2pIwWa2dOS7gaeA7aH7U8H/grcIemKkHZzWOVm4HeSFhMdAZ0TtvOypLuIAth24KIQNJH0ZWA2Uc+7GWb2ctjWt3OU4ZxzZW9BTR2HjRlKZUU+HaOLI58gNAJ4i6hnWZoRnUbrEDO7DLgsI3kJUc+2zLxNwNk5tnMlcGWW9PuB+7OkZy3DOefK3bbmFl5a3sD5x+2ddFVayWcA0wuKURHnnHOF88qKDWzd3sKRe3WvoT/zGcB0rKR7Ja2WtErSPZLGFqNyzjnnusb82joAjhhbYkEI+C1RJ4A9ie6r+b+Q5pxzrkQsqKljt4F9GDu8f9JVaSWfIDTKzH5rZtvD6xZgVIHr5ZxzrgvNr6njyHHDus3o2Wn5BKG1ks5Lj8sm6TyijgrOOedKQEPTNv69ppEjxnWvU3GQXxD6LPBxYCXRGG8fIxrU1DnnXAl4sbYeMziyGwahfLpojzOzM+MJkk4kuunUOedcNze/JuqUcPjYoQnXZGf5HAn9Ms8055xz3dCCmjr2GTmQYQO63yPUch4JSTqeaEy3UZIuji0aQjQSgXPOuW7OzJhfU8eJ+49MuipZtXU6rg8wKOQZHEtvYMf4a84557qxlQ1NrN6whSO64ak4aCP1ByMxAAAYx0lEQVQImdmjwKOSbjGz9h884ZxzrtuZvyy6HnRkHo9uSEI+14RukvR2lwpJwyXNLmCdnHPOdZH5tXX0rhAHjR7cfuYE5BOERppZXXrGzNYDuxeuSs4557rKgpo6Dh49hL6V3fNSfj5BqEXSXukZSXuT44mkzjnnuo/mFuPF2vpueX9QWj73CX0XeFzSo2H+PcC0wlXJOedcV1i8upGNW5u75UgJafk8yuEBSROB4wABXzeztQWvmXPOuU5ZEG5S7c5HQvk8ykHAZGCimf0fMECSPxjOOee6uedr6hjSr5Lxuw1Muio55XNN6FfA8cC5YX4DcH1nCpU0TNLdkl6RtEjS8ZJGSJoj6bXwd3jIK0nXSlos6YVwVJbeztSQ/zVJU2PpR0t6MaxzbQik5CrDOefK0YKaOo4YN4xevbrXyNlx+QShY83sIqAJ3u4d19mxH64BHjCzA4EjgEXAJcDDZjYBeDjMA5wOTAivacANEAUUokeEH0v0yO7LYkHlhpA3vd7kkJ6rDOecKyubtzbz6qoN3fpUHOQXhLZJqiD0iJM0CmjpaIGShhB1brgZwMy2hi7gU4CZIdtM4KwwPQW41SJPAcMkjQZOA+aY2boQGOcAk8OyIWb2pJkZcGvGtrKV4ZxzZeWl5fU0t1i3e5Jqpnx6x10L3AvsLulKoiF7vteJMvcF1gC/lXQEMA/4GlBlZisAzGyFpPS9SGOAmtj6tSGtrfTaLOm0UUYrkqYRegBWVVWRSqVaLW9sbNwprZSVW3ug/NpUbu2B8mtTd2vPA0u3AbCxZiGp1Ys6tI1itCmf3nG3SZoHnEzUO+4sM+tYi3aUORH4ipk9Leka2j4tlu1kpnUgPW9mNh2YDjBp0iSrrq5utTyVSpGZVsrKrT1Qfm0qt/ZA+bWpu7Xnj8ufY8ywOqac9r4Ob6MYbcqnd9w1wAgzu97MrutkAILoyKTWzJ4O83cTBaVV4VQa4e/qWP5xsfXHAsvbSR+bJZ02ynDOubKyIDzOu7vL55rQc8D3Qk+zn0ma1JkCzWwlUCPpHSHpZGAhMAtI93CbCtwXpmcB54decscB9eGU2mzg1DCW3XDgVGB2WLZB0nGhV9z5GdvKVoZzzpWNtY1bqF2/uSSCUD6n42YCM0NvtI8CP5G0V+hh1lFfAW6T1AdYQvS48F7AXZIuJHpq69kh7/3AGcBiYFPIi5mtk/QjYG7I90MzWxemvwTcAvQH/hZeAFflKMM558pG+ibV7jxSQlo+HRPS9gcOBMYTHbl0mJnNB7IdUZ2cJa8BF+XYzgxgRpb0Z4FDs6S/la0M55wrJwtq6qjoJQ4dMyTpqrQrn2tCP5H0GvBD4CXgaDP7UMFr5pxzrkOer6njgKrBDOizK8cZycinhkuB4328OOec6/7MjAU1dXzg8NFJVyUv+VwTurEYFXHOOdd5r7+1iYam7SXRKQHy6x3nnHOuRMyvWQ+URqcE8CDknHNlZUFNPQP6VDBh9+75OO9MeQUhSe+SdEGYHiVpn8JWyznnXEfMr6njsDFDqejGI2fH5dM77jLg28ClIak38PtCVso559yu27K9mYXLG0rmehDkdyT0YeBMYCOAmS0HSuM4zznnepBXVmxga3NL2QWhreGG0fSjHLrvI/qcc64Hm19CIyWk5ROE7pL0a6Ln+HweeAj4TWGr5ZxzblctqKlj98F9GT20X9JVyVs+9wn9XNIpQAPwDuD/mdmcgtfMOefcLplfGz3OOxq7uTTkNaZDCDoeeJxzrpuq37SNJWs28tGJY9vP3I20G4QkbWDHQ+H6EPWO22hm3X9kPOec6yFeeDO6HlRKnRIgv9NxrXrCSToLOKZgNXLOObfLFtTUIcFhY4cmXZVdsssjJpjZn4GTClAX55xzHTS/po79Rg1iSL/eSVdll+RzOu4jsdleRM8BshzZnXPOFZmZMb+mnvceMCrpquyyfI6EPhR7nQZsAKZ0tmBJFZKel/SXML+PpKclvSbpzvDUVST1DfOLw/LxsW1cGtJflXRaLH1ySFss6ZJYetYynHOulL1Zt5m1jVs4clxpnYqDPIKQmV0Qe33ezK40s9VdUPbXgEWx+Z8AV4fHhq8HLgzpFwLrzWx/4OqQD0kHA+cAhwCTgV+FwFYBXA+cDhwMnBvytlWGc86VrAU19QAcOW54wjXZdTlPx0n6JW2cdjOzr3a0UEljgQ8AVwIXK+rUfhLwyZBlJnA5cAPRUdflIf1u4LqQfwpwh5ltAZZKWsyODhOLzWxJKOsOYIqkRW2U4ZxzJWtBbR19Knvxjj1Kb0S1tq4JPRv+nkh0RHFnmD8bmNfJcn8BfIsdY9DtBtSZ2fYwXwuMCdNjgBoAM9suqT7kHwM8FdtmfJ2ajPRj2ymjFUnTgGkAVVVVpFKpVssbGxt3Sitl5dYeKL82lVt7oPzalGR7Ui9uZq9B8MTjj3XpdovRppxByMxmAkj6DPA+M9sW5m8EHuxogZI+CKw2s3mSqtPJ2arQzrJc6dlOMbaVf+dEs+nAdIBJkyZZdXV1q+WpVIrMtFJWbu2B8mtTubUHyq9NSbWnoWkbbzz0EJ86dm+qqw9uf4VdUIw25TNiwp5ERyzrwvygkNZRJwJnSjoD6AcMIToyGiapMhypjAWWh/y1wDigVlIlMDTUJZ2eFl8nW/raNspwzrmSdM+8WrZsb+HDR2U9sdPt5dM77irgeUm3SLoFeA74744WaGaXmtlYMxtP1LHg72b2KeAR4GMh21TgvjA9K8wTlv89jOo9Czgn9J7bB5gAPAPMBSaEnnB9Qhmzwjq5ynDOuZLT0mL87sk3OGqvYSV3k2paPr3jfkt0TeXe8Do+faqui32bqJPCYqLrNzeH9JuB3UL6xcAloV4vA3cBC4EHgIvMrDkc5XwZmE3U++6ukLetMpxzruQ8vngtS9ZuZOrx45OuSoe11TvuQDN7RdLEkJS+2L+npD3N7LnOFm5mKSAVppeQZTggM2si6gyRbf0riXrYZabfD9yfJT1rGc45V4puffINRg7qw+mH7ZF0VTqsrWtCFxP1EPufLMsMH7rHOecSU7NuEw+/soqLqvenb2VF0tXpsLZ6x00Lf99XvOo455zLx++ffoNeEp88dq+kq9Ip7V4TknS2pMFh+nuS/iTpqMJXzTnnXDZN25q5c24Npx5cxZ7D+iddnU7Jp3fc981sg6R3EY0dNxO4sbDVcs45l8usBcup27SN80u4Q0JaPkGoOfz9AHCDmd1H9HA755xzRWZm3Prk6xxQNYjj9h2RdHU6LZ8g9KakXwMfB+6X1DfP9ZxzznWx52vqeOnNBj59/HiiYTRLWz7B5ONE99xMNrM6YATwzYLWyjnnXFa3PvE6g/tW8pESHSEhUz43q24CVgPvCknbgdcKWSnnnHM7W7NhC399cQUfPXosA/vmM+pa95dP77jLiEYauDQk9QZ+X8hKOeec29kdzyxjW7Px6eP3TroqXSaf03EfBs4ENgKY2XJ2PILBOedcEWxvbuG2p5fx7gkj2W/UoKSr02XyCUJbw+CfBiBpYGGr5JxzLtOchatY2dBUFt2y4/IJQneF3nHDJH0eeAi4qbDVcs45FzfzydcZM6w/Jx24e9JV6VLtXtkys59LOgVoAN4B/D8zm1PwmjnnnAPg1ZUbeGrJOi45/UAqepV+t+y4vLpXhKAzB0BShaRPmdltBa2Zc845AG598nX6VvbiE5PGtZu31OQ8HSdpiKRLJV0n6VRFvgwsIbp3yDnnXIE1NG3j3uff5ENH7MnwgeU3WE1bR0K/A9YDTwKfI7pBtQ8wxczmF6FuzjnX490zr5ZNW5tL+sF1bWkrCO1rZocBSLoJWAvsZWYbilIz55zr4crh8d3taat33Lb0hJk1A0u7IgBJGifpEUmLJL0s6WshfYSkOZJeC3+Hh3RJulbSYkkvxJ70iqSpIf9rkqbG0o+W9GJY51qFAZZyleGcc91ROTy+uz1tBaEjJDWE1wbg8PS0pIZOlLkd+C8zOwg4DrhI0sHAJcDDZjYBeDjMA5wOTAivacANEAUU4DLgWKJHdl8WCyo3hLzp9SaH9FxlOOdct3Prk6+z28DSfnx3e3IGITOrMLMh4TXYzCpj00M6WqCZrTCz58L0BmARMAaYQvSsIsLfs8L0FOBWizxFdL/SaKJnG80xs3Vmtp6o997ksGyImT0ZbrK9NWNb2cpwzrluJXp892rOPWavkn58d3sSHQFP0njgKOBpoMrMVkAUqCSl78gaA9TEVqsNaW2l12ZJp40yMus1jehIiqqqKlKpVKvljY2NO6WVsnJrD5Rfm8qtPVB+berK9rSY8YvnttAL2KflTVKpFV2y3V1VjH2UWBCSNAi4B/hPM2to47kY2RZYB9LzZmbTgekAkyZNsurq6lbLU6kUmWmlrNzaA+XXpnJrD5Rfm7qyPVf8ZSEvrFnKFWcdykePS26w0mLso0QeTiepN1EAus3M/hSSV4VTaYS/q0N6LRC/Q2sssLyd9LFZ0tsqwznnuoU75y7jpseX8pkTxnNeggGoWIoehEJPtZuBRWb2v7FFs4B0D7epwH2x9PNDL7njgPpwSm02cKqk4aFDwqnA7LBsg6TjQlnnZ2wrWxnOOZe4p5e8xff+/BLvnjCS733goKSrUxRJnI47Efg08KKk9E2v3wGuIhos9UJgGXB2WHY/cAawGNgEXABgZusk/QiYG/L90MzWhekvAbcA/YG/hRdtlOGcc4la9tYmvvj7eYwbMYDrPjmRyopETlQVXdGDkJk9TvbrNgAnZ8lvwEU5tjUDmJEl/Vng0Czpb2UrwznnkrShaRsXzpxLi8HNU9/J0P69k65S0fSMUOucc91Uc4vx1dufZ+najdxw3kT2GdmzHtlWHg8pd865EvXj+xfxyKtruPLDh3LCfiOTrk7R+ZGQc84lJN4T7lPHln9PuGw8CDnnXAKe6oE94bLxIOScc0W27K1NfKkH9oTLpue23DnnEtAQ6wk3o4f1hMvGOyY451yRrG5o4j/vnM/StRu59cJjGN/DesJl40HIOecKbHtzCzOffIOr5/yLrc0t/Pgjh/XInnDZeBByzrkCevb1dXzvzy/xysoNvPeAUfzgzEP8CCjGg5BzzhXA2sYtXPW3V7h7Xi17Du3HjedN5LRD9qCNJwb0SB6EnHOuCzW3GH94Zhk/e+AVNm1t5ovv3Y+vnrw/A/r4z202/q4451wXWVBTx/fve4kXaus5ft/d+NFZh7D/7oOTrla35kHIOec6adlbm5j58hZSs//JyEF9ueacIznziD391FsePAg559wuamkxFtTW8dCiVcxZuIp/rWpEwGdOHM/XTzmAIf169r0/u8KDkHPO5aFpWzNP/Hstcxau4qFFq1mzYQsVvcQx40fw/Q/uxZCGpZx9xiFJV7PkeBByzrkszIzVG7bw2L/WMGfhKv7x2lo2b2tmUN9K3nvAKE45uIrqd4xi2IA+AKRSbyRc49LkQcg516Nt3d7CsnUbWbx6I/9e08iSNem/jTQ0bQdg9NB+fOzosZxycBXH7juCvpUVCde6fPTIICRpMnANUAHcZGZXJVwl51wXMzM2bm1m/catrNu4lXWbtr49vXrDFpasaeTfazaybN0mmlvs7fX2GNKPfUcNZMqRY9h31EDeOX4Eh+w5xDsZFEiPC0KSKoDrgVOAWmCupFlmtjDZmjlXnsyMFoMWM5pbDDNoNmPjNmPNhi1sb2lhe7OxrbmF7S3hb5jfFv5u3tZM07ZmNm9tZtPW5lbzm2N/6zdvY93GrazftJX1G7extbkla536VPZi35EDOWj0YD54+Gj2GzWIfUcNZN9RgxjUt8f9LCaqJ77bxwCLzWwJgKQ7gClAlwahR/+1hiv+UhpxbeOmTQx87tGkq9FlDNi0cRMDunmbzKz9TMHGTZsYMC+VsYEs28yxfXs7LZ7X3p5vlW7WKn86n6Xnw/J0cNkxH003p4NOSzTdZjMffqiNhW2ToH/viujVp4J+vSsY2r83Y4cP4PCxQxkxsC8jBvZm+IA+jBjYh+ED+zBiQPR3SL9KP7LpJnpiEBoD1MTma4Fj4xkkTQOmAVRVVZFKpVptoLGxcae0TIvXNzNU2zpf2yIY1LeFCm1OuhpdanD/Fip7df825fszOLR/C5UVTXmtn+u3VTtNgFCrbaTXzZUmhb8504XUi16CXorSe5FtWmzfuoUB/ftSIaJXL1EpqOgVzVf2ggqJyl7QpwL69BJ9K6O/fSqgdy8yAokBW8MraAEao1cD0ev17G9Pp+Xzu1BqitGmnhiEsn1FW/2vZmbTgekAkyZNsurq6laZU6kUmWmZqoHPdbyORZVPe0pNubWp3NoD5demcmsPFKdNPfGhdrXAuNj8WGB5QnVxzrkerScGobnABEn7SOoDnAPMSrhOzjnXI/W403Fmtl3Sl4HZRF20Z5jZywlXyznneqQeF4QAzOx+4P6k6+Gccz1dTzwd55xzrpvwIOSccy4xHoScc84lxoOQc865xGhXhg7piSStATLHaB8JrE2gOoVSbu2B8mtTubUHyq9N5dYe6Fyb9jazUe1l8iDUAZKeNbNJSdejq5Rbe6D82lRu7YHya1O5tQeK0yY/Heeccy4xHoScc84lxoNQx0xPugJdrNzaA+XXpnJrD5Rfm8qtPVCENvk1Ieecc4nxIyHnnHOJ8SDknHMuMR6EdoGkyZJelbRY0iVJ16crSHpd0ouS5kt6Nun6dISkGZJWS3opljZC0hxJr4W/w5Os467I0Z7LJb0Z9tN8SWckWcddIWmcpEckLZL0sqSvhfRS3ke52lSS+0lSP0nPSFoQ2vODkL6PpKfDProzPP6ma8v2a0L5kVQB/As4hejBeHOBc81sYaIV6yRJrwOTzKxkb7KT9B6ihzjfamaHhrSfAuvM7KrwD8NwM/t2kvXMV472XA40mtnPk6xbR0gaDYw2s+ckDQbmAWcBn6F091GuNn2cEtxPip6TPtDMGiX1Bh4HvgZcDPzJzO6QdCOwwMxu6Mqy/Ugof8cAi81siZltBe4ApiRcJweY2WPAuozkKcDMMD2T6AeiJORoT8kysxVm9lyY3gAsAsZQ2vsoV5tKkkUaw2zv8DLgJODukF6QfeRBKH9jgJrYfC0l/KGLMeBBSfMkTUu6Ml2oysxWQPSDAeyecH26wpclvRBO15XMqas4SeOBo4CnKZN9lNEmKNH9JKlC0nxgNTAH+DdQZ2bbQ5aC/OZ5EMqfsqSVw7nME81sInA6cFE4FeS6nxuA/YAjgRXA/yRbnV0naRBwD/CfZtaQdH26QpY2lex+MrNmMzsSGEt05uegbNm6ulwPQvmrBcbF5scCyxOqS5cxs+Xh72rgXqIPXzlYFc7bp8/fr064Pp1iZqvCj0QL8BtKbD+F6wz3ALeZ2Z9Ccknvo2xtKvX9BGBmdUAKOA4YJin9BO6C/OZ5EMrfXGBC6C3SBzgHmJVwnTpF0sBwURVJA4FTgZfaXqtkzAKmhumpwH0J1qXT0j/WwYcpof0ULnrfDCwys/+NLSrZfZSrTaW6nySNkjQsTPcH3k90nesR4GMhW0H2kfeO2wWhu+UvgApghpldmXCVOkXSvkRHPwCVwB9KsU2SbgeqiYadXwVcBvwZuAvYC1gGnG1mJXGxP0d7qolO8RjwOvCF9PWU7k7Su4B/AC8CLSH5O0TXUEp1H+Vq07mU4H6SdDhRx4MKooOTu8zsh+E34g5gBPA8cJ6ZbenSsj0IOeecS4qfjnPOOZcYD0LOOecS40HIOedcYjwIOeecS4wHIeecc4nxIORckUhqDH/HS/pkF2/7OxnzT3Tl9p0rFA9CzhXfeGCXglAYxb0trYKQmZ2wi3VyLhEehJwrvquAd4fnzXw9DBz5M0lzw8CXXwCQVB2eWfMHopsikfTnMNjsy+kBZyVdBfQP27stpKWPuhS2/ZKi50Z9IrbtlKS7Jb0i6bYwCoBzRVXZfhbnXBe7BPiGmX0QIASTejN7p6S+wD8lPRjyHgMcamZLw/xnzWxdGFplrqR7zOwSSV8Og09m+gjRHfxHEI3AMFfSY2HZUcAhROOB/RM4keg5Ms4VjR8JOZe8U4HzwzD6TwO7ARPCsmdiAQjgq5IWAE8RDag7gba9C7g9DKq5CngUeGds27VhsM35RKcJnSsqPxJyLnkCvmJms1slStXAxoz59wPHm9kmSSmgXx7bziU+Blgz/nvgEuBHQs4V3wZgcGx+NvCl8GgAJB0QRjXPNBRYHwLQgURD7adtS6+f4THgE+G60yjgPcAzXdIK57qA/+fjXPG9AGwPp9VuAa4hOhX2XOgcsIbsj1F+APiipBeAV4lOyaVNB16Q9JyZfSqWfi9wPLCAaGTnb5nZyhDEnEucj6LtnHMuMX46zjnnXGI8CDnnnEuMByHnnHOJ8SDknHMuMR6EnHPOJcaDkHPOucR4EHLOOZeY/w8tUx+0TixnrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "init_c_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "i = 0\n",
    "next_norm_values = []\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)  \n",
    "while  residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    next_norm_values.append(residue_vector_norm)\n",
    "    i += 1\n",
    "    print(\".\", end=\"\")\n",
    "    \n",
    "# Printing details about this method time, iteration count etc.\n",
    "print_details(A, B, x, i, toc())\n",
    "\n",
    "\n",
    "# Plotting the value of residue vector norm\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(range(1, i + 1), next_norm_values)\n",
    "\n",
    "ax.set(xlabel='Iteration', ylabel='Residue vector second norm',\n",
    "       title='Residue vector second norm values in next iterations')\n",
    "ax.grid()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gauss-Seidl's method also does have requirements regarding input matrix for it to converge.\n",
    "In this case those reqiurements are not met and the solution vector won't converge to accurate\n",
    "solution. The chart shows exponential growth of residue vector value in respect to iteration \n",
    "count.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Exercise C\n",
    "\n",
    "**LU Factorization method**\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "def matrix_decomposition(A: np.ndarray) -> (np.ndarray, np.ndarray):\n",
    "    m = A.shape[0]\n",
    "    U = A.copy()\n",
    "    L = np.identity(m)\n",
    "    for k in range(0, m - 1):\n",
    "        for j in range(k + 1, m):\n",
    "            L[j][k] = U[j][k]/U[k][k]\n",
    "            for i in range(k, m):\n",
    "                U[j][i] = U[j][i] - L[j][k]*U[k][i]\n",
    "    return L, U\n",
    "\n",
    "def forward_substitution(L: np.ndarray, B: np.ndarray) -> np.ndarray:\n",
    "    x = np.ndarray(shape=(L.shape[0],1))\n",
    "    for e in range(L.shape[0]):\n",
    "        curr_calc_coeff_val = 0.0\n",
    "        curr_eq = L[e].copy()\n",
    "        # adding every calculated coeff so far\n",
    "        for i in range(e):\n",
    "            curr_calc_coeff_val -= curr_eq[i] * x[i]\n",
    "        x[e] = B[e] + curr_calc_coeff_val\n",
    "    return x\n",
    "\n",
    "def backward_substitution(U: np.ndarray, y: np.ndarray) -> np.ndarray:\n",
    "    n = U.shape[0]\n",
    "    x = np.ndarray(shape=(n, 1))\n",
    "    for e in range(n):\n",
    "        curr_calc_coeff_val = 0.0\n",
    "        curr_eq = U[-(e+1)].copy()\n",
    "        # adding every calculated coeff so far\n",
    "        for i in range(e):\n",
    "            curr_calc_coeff_val -= curr_eq[-(i+1)] * x[-(i+1)]\n",
    "        x[-(e+1)] = (y[-(e+1)] + curr_calc_coeff_val) / curr_eq[-(e+1)]\n",
    "    return x\n",
    "        \n",
    "    \n",
    "def solve_lu_equation(L: np.ndarray, U: np.ndarray, B: np.ndarray) -> np.ndarray:\n",
    "    n = L.shape[0]\n",
    "    # LUx = b -> L(Ux) = b\n",
    "    # Ly = b && y = Ux\n",
    "    # Solving for y\n",
    "    y = forward_substitution(L, B)\n",
    "    x = backward_substitution(U, y)\n",
    "    return x\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "metadata": false,
     "name": "#%% md\n"
    }
   },
   "source": [
    "**LU Factorization method testing**\n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "metadata": false,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Residue vector norm value:  7.245968407401433e-15\n"
     ]
    }
   ],
   "source": [
    "init_c_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "x = solve_lu_equation(L,U,B)\n",
    "\n",
    "print(\"Residue vector norm value: \", calculate_residue_vector_norm(A, x, B))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "LU Factorization is a direct method of solving equations and therefore the requirements \n",
    "posed by the iterative methods do not hold. Equation system from ex.C is solved by LU Factorization \n",
    "method almost perfectly accurate (using floats account for the imperfection).\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Initializing last exercise\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "is_executing": true,
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "100 x 100: \n",
      "........................\n",
      "................\n"
     ]
    }
   ],
   "source": [
    "jacobis = []\n",
    "gauss = []\n",
    "LU_fac = []\n",
    "\n",
    "print(\"100 x 100: \")\n",
    "N =  100\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "print()\n",
    "jacobis.append(toc())\n",
    "\n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "print()\n",
    "gauss.append(toc())\n",
    "\n",
    "tic()\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "solve_lu_equation(L,U,B)\n",
    "LU_fac.append(toc())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "500 x 500\n",
      "........................\n",
      "................\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print(\"500 x 500\")\n",
    "N = 500  \n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "    \n",
    "jacobis.append(toc())\n",
    "print()\n",
    "\n",
    "\n",
    "generate_equation_system()\n",
    "tic()\n",
    "\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "gauss.append(toc())\n",
    "print()\n",
    "\n",
    "\n",
    "tic()\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "solve_lu_equation(L,U,B)\n",
    "LU_fac.append(toc())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1000 x 1000: \n",
      "....................\n",
      ".............\n"
     ]
    }
   ],
   "source": [
    "\n",
    "print(\"1000 x 1000: \")\n",
    "N = 1000  \n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "    \n",
    "jacobis.append(toc())\n",
    "print()\n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "gauss.append(toc())\n",
    "print()\n",
    "\n",
    "\n",
    "\n",
    "tic()\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "solve_lu_equation(L,U,B)\n",
    "LU_fac.append(toc())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2000 x 2000: \n",
      "......................\n",
      "...............\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-18-16957cea6000>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     37\u001b[0m \u001b[0minit_a_exercise\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     38\u001b[0m \u001b[0mgenerate_equation_system\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 39\u001b[1;33m \u001b[0mL\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mU\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mmatrix_decomposition\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mA\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     40\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     41\u001b[0m \u001b[0msolve_lu_equation\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mL\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mU\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mB\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<ipython-input-14-8397e733509c>\u001b[0m in \u001b[0;36mmatrix_decomposition\u001b[1;34m(A)\u001b[0m\n\u001b[0;32m      7\u001b[0m             \u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mU\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mU\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      8\u001b[0m             \u001b[1;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 9\u001b[1;33m                 \u001b[0mU\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mU\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mL\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mj\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mU\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     10\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mL\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mU\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     11\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "\n",
    "print(\"2000 x 2000: \")\n",
    "N = 2000  \n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "# residue_vector_norm_values = [float]\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "    \n",
    "jacobis.append(toc())\n",
    "print()\n",
    "\n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "tic()\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "    \n",
    "gauss.append(toc())\n",
    "print()\n",
    "\n",
    "\n",
    "tic()\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "solve_lu_equation(L,U,B)\n",
    "LU_fac.append(toc())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "is_executing": false,
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "print(\"3000 x 3000: \")\n",
    "\n",
    "N = 3000\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "\n",
    "tic()\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = jacobis_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "jacobis.append(toc())\n",
    "\n",
    "\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "\n",
    "tic()\n",
    "residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "i = 0\n",
    "while residue_vector_norm > epsilon and i < iter_limit:\n",
    "    x = gauss_seidl_next_iteration(A, B, x)\n",
    "    i += 1\n",
    "    residue_vector_norm = calculate_residue_vector_norm(A, x, B)\n",
    "    print(\".\", end=\"\")\n",
    "\n",
    "gauss.append(toc())\n",
    "\n",
    "tic()\n",
    "init_a_exercise()\n",
    "generate_equation_system()\n",
    "L, U = matrix_decomposition(A)\n",
    "\n",
    "solve_lu_equation(L,U,B)\n",
    "LU_fac.append(toc())\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "<h1>Summary</h1>\n",
    "Both iterative methods have uses. While Gauss performs better as a single threaded algorithm, it is not possible to make it concurrent. A significant drawback of these methods is their requirement regarding the input matrix - it's spectral radius has to be smaller than 1. This limits the use of previously mentioned iterative methods. This limitation is not true for direct solving method - LU Factorization. However Factorization is significantly slower for bigger matrices rendring this method almost useless on slower computational units. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the value of residue vector norm\n",
    "fig, ax = plt.subplots()\n",
    "its = [100,500,1000,2000,3000]\n",
    "\n",
    "ax.plot(its, jacobis)\n",
    "\n",
    "ax.set(xlabel='Matrix size', ylabel='Time to meet specified precision',\n",
    "       title='Calculation time in respect to input matrix size (Jacobis)')\n",
    "ax.grid()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(its, gauss)\n",
    "\n",
    "ax.set(xlabel='Matrix size', ylabel='Time to meet specified precision',\n",
    "       title='Calculation time in respect to input matrix size (Gauss)')\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  },
  "pycharm": {
   "stem_cell": {
    "cell_type": "raw",
    "metadata": {
     "collapsed": false
    },
    "source": []
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}