lectures/lecture12/lecture-12.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 12:  Great Iterative Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Previous part\n",
    "\n",
    "- Concept of **iterative methods** for linear systems:\n",
    "    - Richardson iteration and its convergence \n",
    "    - Chebyshev iteration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Plan for the next part\n",
    "\n",
    "- Lanczos and Arnoldi orthogonalization of Krylov subspaces, optimality result for Krylov subspaces\n",
    "- Main iterative methods: conjugate gradient, GMRES, etc\n",
    "- Convergence estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solution of linear systems and minimization of functionals\n",
    "\n",
    "- Instead of solving a linear system, we can minimize the **residual:**\n",
    "\n",
    "$$R(x) = \\Vert A x - f \\Vert^2_2.$$\n",
    "\n",
    "- The condition $\\nabla R(x) = 0$ gives\n",
    "\n",
    "$$A^* A x = A^* f,$$\n",
    "\n",
    "thus it has squared condition number, so direct minimization of the residual by standard optimization methods is rarely used.\n",
    "\n",
    "- For the symmetric positive definite case there is a much simpler functional."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Energy functional\n",
    "\n",
    "Let $A = A^* > 0$, then the following functional \n",
    "\n",
    "$$\\Phi(x) = (Ax, x)  - 2(f, x)$$\n",
    "\n",
    "is called **energy functional**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Properties of energy functional\n",
    "- It is strictly convex (check!)\n",
    "\n",
    "$$ \\Phi(\\alpha x + (1 - \\alpha)y) < \\alpha \\Phi(x) + (1 - \\alpha) \\Phi(y)$$\n",
    "\n",
    "- Since it is strictly convex, it has unique local minimum, which is also global\n",
    "\n",
    "- Its global minimum $x_*$ satisfies\n",
    "\n",
    "$$A x_* = f.$$\n",
    "\n",
    "Indeed, \n",
    "\n",
    "$$\\nabla \\Phi = 2(Ax - f).$$ \n",
    "\n",
    "and the first order optimality condition $\\nabla \\Phi (x_*) = 0$ yields\n",
    "\n",
    "$$A x_* = f.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Approximation of the solution by a subspace\n",
    "\n",
    "- Given a linear $M$-dimensional subspace $\\{y_1, \\dots, y_M\\}$, we want to find an approximate solution in this basis, i.e. \n",
    "\n",
    "$$A x \\approx f, \\quad x = x_0 +  \\sum_{k=1}^M c_k y_k,$$\n",
    "\n",
    "where $c$ is the vector of coefficients.\n",
    "\n",
    "- In the symmetric positive definite case  we need to minimize \n",
    "\n",
    "$$(Ax, x) - 2(f, x)$$ \n",
    "\n",
    "subject to $$x = x_0 + Y c,$$\n",
    "\n",
    "where $Y=[y_1,\\dots,y_M]$ is $n \\times M$ and vector $c$ has length $M$.\n",
    "\n",
    "- Using the representation of $x$, we have the following minimization for $c$:\n",
    "\n",
    "$$\\widehat{\\Phi}(c) = (A Y c, Y c) + 2 (Y^*Ax_0, c) - 2(f, Y c) = (Y^* A Y c, c) - 2(Y^* (f - Ax_0), c).$$\n",
    "\n",
    "- Note that this is the same functional, but for the **Galerkin projection** of $A$\n",
    "\n",
    "$$Y^* A Y c = Y^* (f - Ax_0) = Y^* r_0,$$\n",
    "\n",
    "which is an $M \\times M$  linear system with symmetric positive definite matrix if $Y$ has full column rank.\n",
    "\n",
    "But how to choose $Y$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Selection of the subspace\n",
    "\n",
    "In the Krylov subspace we generate the whole subspace from a single vector $r_0 = f - Ax_0$:\n",
    "\n",
    "$$y_0\\equiv k_0 = r_0, \\quad y_1\\equiv k_1 = A r_0, \\quad y_2\\equiv k_2 = A^2 r_0, \\ldots, \\quad y_{M-1}\\equiv k_{M-1} = A^{M-1} r_0.$$\n",
    "\n",
    "This gives the **Krylov subpace** of the $M$-th order\n",
    "\n",
    "$$\\mathcal{K}_M(A, r_0) = \\mathrm{Span}(r_0, Ar_0, \\ldots, A^{M-1} r_0).$$\n",
    "\n",
    "- It is known to be quasi-optimal space given only matrix-vector product operation.\n",
    "- Key reference here is \"On the numerical solution of equation by which are determined in technical problems the frequencies of small vibrations of material systems\",  A. N. Krylov, 1931, [text in russian](http://www.mathnet.ru/links/a2431bc65c0764da9a1ae95c73f741f6/im5215.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solution $x_*$ lies in the Krylov subspace: $x_* \\in \\mathcal{K}_n(A, f)$\n",
    "\n",
    "- According to [Cayley–Hamilton theorem](https://en.wikipedia.org/wiki/Cayley%E2%80%93Hamilton_theorem): $p(A) = 0$, where $p(\\lambda) = \\det(A - \\lambda I)$\n",
    "- $p(A)f = A^nf + a_1A^{n-1}f + \\ldots + a_{n-1}Af + a_n f = 0$\n",
    "- $A^{-1}p(A)f = A^{n-1}f + a_1A^{n-2}f + \\ldots + a_{n-1}f + a_nA^{-1}f = 0$\n",
    "- $x_* = A^{-1}f = -\\frac{1}{a_n}(A^{n-1}f + a_1A^{n-2}f + \\ldots + a_{n-1}f)$\n",
    "- Thus, $x_* \\in \\mathcal{K}_n(A, f)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Ill-conditioned of the natural basis\n",
    "\n",
    "The natural basis in the Krylov subspace is very **ill-conditioned**, since \n",
    "\n",
    "$$k_i = A^i r_0 \\rightarrow \\lambda_\\max^i v,$$\n",
    "\n",
    "where $v$ is the eigenvector, corresponding to the maximal eigenvalue of $A$, i.e. $k_i$ become more and more collinear for large $i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Condition number = 595184527.1939951\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.sparse as spsp\n",
    "%matplotlib inline\n",
    "\n",
    "n = 100\n",
    "ex = np.ones(n);\n",
    "A = spsp.spdiags(np.vstack((-ex,  2*ex, -ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "f = np.ones(n)\n",
    "x0 = np.random.randn(n)\n",
    "\n",
    "subspace_order = 10\n",
    "krylov_vectors = np.zeros((n, subspace_order))\n",
    "krylov_vectors[:, 0] = f - A.dot(x0)\n",
    "for i in range(1, subspace_order):\n",
    "    krylov_vectors[:, i] = A.dot(krylov_vectors[:, i-1])\n",
    "    \n",
    "s = np.linalg.svd(krylov_vectors, compute_uv=False)\n",
    "print(\"Condition number = {}\".format(s.max() / s.min()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Solution:** Compute orthogonal basis in the Krylov subspace."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Good basis in a Krylov subspace\n",
    "\n",
    "In order to have stability, we first orthogonalize the vectors from the Krylov subspace using **Gram-Schmidt** orthogonalization process (or, QR-factorization).\n",
    "\n",
    "$$K_j = \\begin{bmatrix} r_0 & Ar_0 & A^2 r_0 & \\ldots & A^{j-1} r_0\\end{bmatrix} = Q_j R_j, $$\n",
    "\n",
    "and the solution will be approximated as \n",
    "\n",
    "$$x \\approx x_0 + Q_j c.$$ \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Short way to Arnoldi relation\n",
    "\n",
    "**Statement.** The Krylov matrix $K_j$ satisfies an important recurrent relation (called **Arnoldi relation**)\n",
    "\n",
    "$$A Q_j = Q_j H_j + h_{j, j-1} q_j e^{\\top}_{j-1},$$\n",
    "\n",
    "where $H_j$ is upper Hessenberg, and  $Q_{j+1} = [q_0,\\dots,q_j]$ has orthogonal columns that spans columns of $K_{j+1}$.\n",
    "\n",
    "Let us prove it (consider $j = 3$ for simplicity):\n",
    "\n",
    "$$A \\begin{bmatrix} k_0 & k_1 & k_2 \\end{bmatrix} = \\begin{bmatrix} k_1 & k_2 & k_3 \\end{bmatrix} = \\begin{bmatrix} k_0 & k_1 & k_2 \\end{bmatrix} \\begin{bmatrix} 0 & 0 & \\alpha_0 \\\\ 1 & 0  & \\alpha_1 \\\\ 0 & 1  & \\alpha_2 \\\\ \\end{bmatrix} + \\begin{bmatrix} 0 & 0 & k_3  - \\alpha_0 k_0 - \\alpha_1 k_1 - \\alpha_2 k_2 \\end{bmatrix}, $$\n",
    " \n",
    "where $\\alpha_s$ will be selected later. \n",
    "Denote $\\widehat{k}_3 = k_3  - \\alpha_0 k_0 - \\alpha_1 k_1 - \\alpha_2 k_2$.\n",
    " \n",
    " In the matrix form,\n",
    " \n",
    " $$A K_3 = K_3 Z + \\widehat k_3 e^{\\top}_2,$$\n",
    " \n",
    " where $Z$ is the **lower shift** matrix with the last column $(\\alpha_0,\\alpha_1,\\alpha_2)^T$, and $e_2$ is the last column of the identity matrix.\n",
    " \n",
    " Let \n",
    " \n",
    " $$K_3 = Q_3 R_3$$ \n",
    " \n",
    " be the QR-factorization. Then,\n",
    " \n",
    " $$A Q_3 R_3 = Q_3 R_3 Z + \\widehat{k}_3 e^{\\top}_2,$$\n",
    " \n",
    " $$ A Q_3 = Q_3 R_3 Z R_3^{-1} + \\widehat{k}_3 e^{\\top}_2 R_3^{-1}.$$\n",
    "\n",
    "\n",
    "Note that \n",
    "\n",
    "$$e^{\\top}_2 R_3^{-1} = \\begin{bmatrix} 0 & 0 & 1 \\end{bmatrix} \\begin{bmatrix} * & * & * \\\\ 0 & * & * \\\\ 0 & 0 & * \\end{bmatrix}  = \\gamma e^{\\top}_2,$$\n",
    " and\n",
    " \n",
    " $$R_3 Z R_3^{-1} = \\begin{bmatrix} * & * & * \\\\* & * & * \\\\  0 & * & * \\\\ \\end{bmatrix},$$      \n",
    "\n",
    "in the general case it will be an **upper Hessenberg matrix** $H$, i.e. a matrix \n",
    "that \n",
    "\n",
    "$$H_{ij} = 0, \\quad \\mbox{if } i > j + 1.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## (Almost) Arnoldi relation\n",
    "\n",
    "Let $Q_j$ be the orthogonal basis in the Krylov subspace, then we have almost the Arnoldi relation\n",
    "\n",
    "$$A Q_j = Q_j H_j +  \\gamma\\widehat{k}_j e^{\\top}_{j-1},$$\n",
    "\n",
    "where $H_j$ is an upper Hessenberg matrix, and \n",
    "\n",
    "$$\\widehat{k}_j = k_j - \\sum_{s=0}^{j-1} \\alpha_s k_s.$$\n",
    "\n",
    "We select $\\alpha_s$ in such a way that \n",
    "\n",
    "$$Q^*_j \\widehat{k}_j = 0.$$\n",
    "\n",
    "Then, $\\widehat{k}_j = h_{j, j-1} q_j,$ where $q_j$ is the last column of $Q_{j+1}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Arnoldi relation: final formula\n",
    "\n",
    "We have \n",
    "\n",
    "$$A Q_j = Q_j H_j + h_{j, j-1} q_j e^{\\top}_{j-1}.$$\n",
    "\n",
    "- This is the crucial formula for the efficient generation of such subspaces.\n",
    "\n",
    "- For non-symmetric case, it is just modified Gram-Schmidt.\n",
    "\n",
    "- For the symmetric case, we have a much simpler form (Lanczos process)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lanczos process\n",
    "\n",
    "If $A = A^*$, then\n",
    "\n",
    "$$Q^*_j A Q_j = H_j, $$\n",
    "\n",
    "thus $H_j$ is hermitian, and thus it is **tridiagonal**, $H_j = T_j$.\n",
    "\n",
    "This gives a short-term recurrence relation to generate the Arnoldi vectors $q_j$ without **full** orthogonalization."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lanczos process (2)\n",
    "\n",
    "$$ A Q_j = Q_j T_j + t_{j, j-1} q_j e^{\\top}_{j-1}.$$\n",
    "\n",
    "\n",
    "In order to get $q_j$, we need to compute just the last column of \n",
    "\n",
    "$$t_{j, j-1} q_j = (A Q_j - Q_j T_j) e_{j-1} = A q_{j-1} - t_{j-1, j-1} q_{j-1} - t_{j-2, j-1} q_{j-2}. $$\n",
    "\n",
    "The coefficients $\\alpha_j = t_{j-1, j-1}$ and $\\beta_j = t_{j-2, j-1}$ can be recovered from orthogonality constraints \n",
    "\n",
    "$(q_j, q_{j-1}) = 0, \\quad (q_j, q_{j-2}) = 0$\n",
    "\n",
    "**All the other constraints will be satisfied automatically!!**\n",
    "\n",
    "And we only need to store two vectors to get the new one."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## From direct Lanczos method to the conjugate gradient\n",
    "\n",
    "We can now get from the Lanczos recurrence to the famous **conjugate gradient** method.\n",
    "\n",
    "We have for $A = A^* > 0$\n",
    "\n",
    "$$A Q_j = Q_j T_j + T_{j, j-1} q_j.$$\n",
    "\n",
    "Recall that when we minimize energy functional in basis $Y$ we get a system $Y^* A Y c = Y^* f,$.\n",
    "Here $Y = Q_j$, so the approximate solution of $Ax \\approx f$ with $x_j = x_0 + Q_j c_j$ can be found by solving a small system\n",
    "\n",
    "$$Q^*_j A Q_j c_j = T_j c_j = Q^*_j r_0 .$$\n",
    "\n",
    "Since $f$ is the first Krylov subspace, then \n",
    "**Note!!!** (recall what the first column in $Q_j$ is)\n",
    "\n",
    "$$Q^*_j r_0  = \\Vert r_0 \\Vert_2^2 e_0 = \\gamma e_0.$$\n",
    "\n",
    "We have a tridiagonal system of equations for $c$:\n",
    "\n",
    "$$T_j c_j = \\gamma e_0$$\n",
    "\n",
    "and $x_j = Q_j c_j$.\n",
    "\n",
    "We could stop at this point, but we want short recurrent formulas instead of solving linear system with matrix $T_j$ at each step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Derivation of the following update formulas is not required on the oral exam!**\n",
    "\n",
    "- Since $A$ is positive definite, $T_j$ is also positive definite, and it allows an LU decomposition\n",
    "\n",
    "- $T_j = L_j U_j$, where $L_j$ is a bidiagonal matrix with ones on the diagonal, $U_j$ is a upper bidiagonal matrix.\n",
    "\n",
    "- We need to define one subdiagonal in $L$ (with elements $c_1, \\ldots, c_{j-1}$), main diagonal of $U_j$ (with elements $d_0, \\ldots, d_{j-1}$ and superdiagonal of $U_j$ (with elements $b_1, \\ldots, b_{j-1}$).\n",
    "\n",
    "- They have convenient recurrences:\n",
    "\n",
    "$$c_i = b_i/d_{i-1}, \\quad d_i = \\begin{cases} a_1, & \\mbox{if } i = 0, \\\\\n",
    "a_i - c_i b_i, & \\mbox{if } i > 0. \\end{cases}$$\n",
    "\n",
    "- For the solution we have\n",
    "\n",
    "$$x_j = Q_j T^{-1}_j \\gamma e_0  = \\gamma Q_j (L_j U_j)^{-1} e_0  = \\gamma Q_j U^{-1}_j L^{-1}_j e_0.$$\n",
    "\n",
    "- We introduce two new quantities:  \n",
    "\n",
    "$$P_j = Q_j U^{-1}_j, \\quad z_j = \\gamma L^{-1}_j e_0.$$\n",
    "\n",
    "- Now we have the following equation for $x_j$:\n",
    "\n",
    "$$ x_j = P_j z_j$$\n",
    "\n",
    "- Due to the recurrence relations, we have\n",
    "\n",
    "$$P_j = \\begin{bmatrix} P_{j-1} & p_j \\end{bmatrix}, $$\n",
    "\n",
    "and \n",
    "\n",
    "$$z_j = \\begin{bmatrix} z_{j-1} \\\\ \\xi_{j} \\end{bmatrix}.$$\n",
    "\n",
    "- For $p_j$ and $\\xi_j$ we have short-term recurrence relations (due to bidiagonal structure)\n",
    "\n",
    "$$p_j = \\frac{1}{d_j}\\left(q_j - b_j p_{j-1} \\right), \\quad \\xi_j = -c_j \\xi_{j-1}.$$\n",
    "\n",
    "- Thus, we arrive at short-term recurrence for $x_j$:\n",
    "\n",
    "$$x_j = P_j z_j = P_{j-1} z_{j-1} + \\xi_j p_j = x_{j-1} + \\xi_j p_j.$$\n",
    "\n",
    "and $q_j$ are found from the Lanczos relation (see slides above).\n",
    "\n",
    "- This method for solving linear systems is called a **direct Lanczos method**. It is closely related to the conjugate gradient method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Direct Lanczos method\n",
    "\n",
    "We have the direct Lanczos method, where we store \n",
    "\n",
    "$$p_{j-1}, q_j, x_{j-1}$$ to get a new estimate of $x_j$.\n",
    "\n",
    "The main problem is with $q_j$: we have the three-term recurrence, but in the floating point arithmetic \n",
    "the orthogonality is can be lost, leading to numerical errors.\n",
    "\n",
    "Let us do some demo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6432256263275\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import scipy as sp\n",
    "import scipy.sparse as spsp\n",
    "from scipy.sparse import csc_matrix\n",
    "\n",
    "n = 128\n",
    "ex = np.ones(n);\n",
    "A = spsp.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "\n",
    "nit = 64\n",
    "q1 = rhs/np.linalg.norm(rhs)\n",
    "q2 = A.dot(q1)\n",
    "q2 = q2 - np.dot(q2, q1)*q1\n",
    "q2 = q2/np.linalg.norm(q2)\n",
    "qall = [q1, q2]\n",
    "for i in range(nit):\n",
    "    qnew = A.dot(qall[-1])\n",
    "    qnew = qnew - np.dot(qnew, qall[-1])*qall[-1]\n",
    "    qnew = qnew/np.linalg.norm(qnew)\n",
    "    qnew = qnew - np.dot(qnew, qall[-2])*qall[-2]\n",
    "    qnew = qnew/np.linalg.norm(qnew)\n",
    "    qall.append(qnew)\n",
    "qall_mat = np.vstack(qall).T\n",
    "print(np.linalg.norm(qall_mat.T.dot(qall_mat) - np.eye(qall_mat.shape[1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conjugate gradient method\n",
    "\n",
    "Instead of $q_j$ (last vector in the modified Gram-Schmidt process), it is more convenient to work with the **residual**\n",
    "\n",
    "$$r_j = f - A x_j.$$\n",
    "\n",
    "The resulting recurrency has the form\n",
    "\n",
    "$x_j = x_{j-1} + \\alpha_{j-1} p_{j-1}$\n",
    "\n",
    "$r_j = r_{j-1} - \\alpha_{j-1}  A p_{j-1}$\n",
    "\n",
    "$p_j = r_j + \\beta_j p_{j-1}$.\n",
    "\n",
    "Hence the name conjugate gradient: to the gradient $r_j$ we add a **conjugate direction** $p_j$.\n",
    "\n",
    "We have **orthogonality** of residuals (check!):\n",
    "\n",
    "$$(r_i, r_j) = 0, \\quad i \\ne j$$\n",
    "\n",
    "and **A-orthogonality** of conjugate directions (check!):\n",
    "\n",
    "$$ (A p_i, p_j) = 0,$$\n",
    "\n",
    "which can be checked from the definition.\n",
    "\n",
    "The equations for $\\alpha_j$ and $\\beta_j$ can be now defined explicitly from these two properties."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CG final formulas\n",
    "\n",
    "We have $(r_{j}, r_{j-1}) = 0 = (r_{j-1} - \\alpha_{j-1} A r_{j-1}, r_{j-1})$, \n",
    "\n",
    "thus\n",
    "\n",
    "$$\\alpha_{j-1} = \\frac{(r_{j-1}, r_{j-1})}{(A r_{j-1}, r_{j-1})}.$$\n",
    "\n",
    "In the similar way, we have\n",
    "\n",
    "$$\\beta_{j-1} = \\frac{(r_j, r_j)}{(r_{j-1}, r_{j-1})}.$$\n",
    "\n",
    "Recall that\n",
    "\n",
    "$x_j = x_{j-1} + \\alpha_{j-1} p_{j-1}$\n",
    "\n",
    "$r_j = r_{j-1} - \\alpha_{j-1}  A p_{j-1}$\n",
    "\n",
    "$p_j = r_j + \\beta_j p_{j-1}$.\n",
    "\n",
    "Only one matrix-by-vector product per iteration.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CG derivation overview\n",
    "\n",
    "- Want to find $x_*$ in Krylov subspace\n",
    "- But natural basis is ill-conditioned, therefore we need orthogonalization\n",
    "- Derive recurrent equation for sequential orthogonalization of the Krylov subspace basis\n",
    "    - Arnoldi process for non-symmetric matrix\n",
    "    - Lanczos process for symmetrix matrix\n",
    "- Clever re-writing of these formulas gives short recurrence "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Some history\n",
    "More details here: https://www.siam.org/meetings/la09/talks/oleary.pdf\n",
    "\n",
    "When Hestenes worked on conjugate bases in 1936, he was advised by a\n",
    "Harvard professor that it was too obvious for publication\n",
    "- CG doesn’t work on slide rules.\n",
    "<img src=\"slide_rule.jpg\">\n",
    "- CG has little advantage over Gauss elimination for computation with\n",
    "calculators.\n",
    "- CG is not well suited for a room of human \"computers\" – too much data exchange.\n",
    "<img src=\"human_computers.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Properties of the CG method\n",
    "\n",
    "- We need to store 3 vectors.\n",
    "\n",
    "- Since it generates $A$-orthogonal sequence $p_1, \\ldots, p_N$, after $n$ steps it should stop (i.e., $p_{N+1} = 0$.) \n",
    "\n",
    "- In practice it does not have this property in finite precision, thus after its invention in 1952 by Hestens and Stiefel it was labeled **unstable**.\n",
    "\n",
    "- In fact, it is a brilliant iterative method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## $A$-optimality\n",
    "\n",
    "\n",
    "Energy functional can be written as \n",
    "\n",
    "$$(Ax, x) - 2(f, x) = (A (x - x_*), (x - x_*)) - (Ax _*, x_*),$$\n",
    "\n",
    "where $A x_* = f$. Up to a constant factor,\n",
    "\n",
    "$$ (A(x - x_*), (x -x_*)) = \\Vert x - x_* \\Vert^2_A$$\n",
    "\n",
    "is the **A-norm** of the error.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence\n",
    "The CG method computes $x_k$ that minimizes the energy functional over the Krylov subspace, i.e. $x_k = p(A)f$, where $p$ is a polynomial of degree $k+1$, so\n",
    "\n",
    "$$\\Vert x_k - x_* \\Vert_A  =  \\inf\\limits_{p} \\Vert \\left(p(A) - A^{-1}\\right) f \\Vert_A. $$\n",
    "\n",
    "Using eigendecomposition of $A$ we have\n",
    "\n",
    "$$A = U \\Lambda U^*, \\quad  g = U^* f,$$ and\n",
    "\n",
    "$\\Vert x - x_* \\Vert^2_A = \\displaystyle{\\inf_p} \\Vert \\left(p(\\Lambda) - \\Lambda^{-1}\\right) g \\Vert_\\Lambda^2 = \\displaystyle{\\inf_p}\n",
    "\\displaystyle{\\sum_{i=1}^n} \\frac{(\\lambda_i p(\\lambda_i) - 1)^2 g^2_i}{\\lambda_i} = \\displaystyle{\\inf_{q, q(0) = 1}} \\displaystyle{\\sum_{i=1}^n} \\frac{q(\\lambda_i)^2 g^2_i}{\\lambda_i}\n",
    "$\n",
    "\n",
    "Selection of the optimal $q$ depends on the eigenvalue distribution.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Absolute and relative error\n",
    "\n",
    "We have\n",
    "\n",
    "$$\\Vert x - x_* \\Vert^2_A \\leq \\sum_{i=1}^n \\frac{g^2_i}{\\lambda_i} \\inf_{q, q(0)=1} \\max_{j} q({\\lambda_j})^2$$\n",
    "\n",
    "The first term is just \n",
    "\n",
    "$$\\sum_{i=1}^n \\frac{g^2_i}{\\lambda_i} = (A^{-1} f, f) = \\Vert x_* \\Vert^2_A.$$\n",
    "\n",
    "And we have relative error bound\n",
    "\n",
    "$$\\frac{\\Vert x - x_* \\Vert_A }{\\Vert x_* \\Vert_A} \\leq \\inf_{q, q(0)=1} \\max_{j} |q({\\lambda_j})|,$$\n",
    "\n",
    "so if matrix has only 2 different eigenvalues, then there exists a polynomial of degree 2 such that $q({\\lambda_1}) =q({\\lambda_2})=0$, so in this case CG converges in 2 iterations.\n",
    "\n",
    "<font color='red'>\n",
    "\n",
    "- If eigenvalues are clustered and there are $l$ outliers, then after first $\\mathcal{O}(l)$ iterations CG will converge as if there are no outliers (and hence the effective condition number is smaller).\n",
    "    \n",
    "- The intuition behind this fact is that after $\\mathcal{O}(l)$ iterations the polynomial has degree more than $l$ and thus is able to zero $l$ outliers.\n",
    "</font>\n",
    "\n",
    "Let us find another useful upper-bound estimate of convergence. \n",
    "Since\n",
    "\n",
    "$$\n",
    "\\inf_{q, q(0)=1} \\max_{j} |q({\\lambda_j})| \\leq \\inf_{q, q(0)=1} \\max_{\\lambda\\in[\\lambda_\\min,\\lambda_\\max]} |q({\\lambda})|\n",
    "$$\n",
    "\n",
    "The last term is just the same as for the Chebyshev acceleration, thus the same \n",
    "upper convergence bound holds:\n",
    "\n",
    "$$\\frac{\\Vert x_k - x_* \\Vert_A }{\\Vert x_* \\Vert_A} \\leq \\gamma \\left( \\frac{\\sqrt{\\mathrm{cond}(A)}-1}{\\sqrt{\\mathrm{cond}(A)}+1}\\right)^k.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Finite termination & clusters\n",
    "\n",
    "1. If $A$ has only $m$ different eigenvalues, CG converges in $m$ iterations (proof in the blackboard).\n",
    "2. If $A$ has $m$ \"clusters\" of eigenvalues, CG converges cluster-by-cluster.\n",
    "\n",
    "As a result, better convergence than Chebyshev acceleration, but slightly higher cost per iteration."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "CG is the method of choice for symmetric positive definite systems:\n",
    "\n",
    "1. $\\mathcal{O}(n)$ memory\n",
    "2. Square root of condition number in the estimates\n",
    "3. Automatic ignoring of the outliers/clusters\n",
    "4. $A$-optimality property"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Non-linear conjugate gradient method\n",
    "\n",
    "- CG minimizes the energy functional, which is quadratic in $x$\n",
    "- CG formulas were used as starting point in developing methods to minimize arbitrary convex function\n",
    "- Most popular CG extensions (so-called non-linear CG method) are \n",
    "    - [Hestenes-Stiefel method](https://www.fing.edu.uy/inco/cursos/numerico/aln/hes_stief1952.pdf)\n",
    "    - [Polak-Ribiere method](http://www.numdam.org/article/M2AN_1969__3_1_35_0.pdf) - original paper in French\n",
    "    - [Fletcher–Reeves method](https://watermark.silverchair.com/070149.pdf?token=AQECAHi208BE49Ooan9kkhW_Ercy7Dm3ZL_9Cf3qfKAc485ysgAAAkMwggI_BgkqhkiG9w0BBwagggIwMIICLAIBADCCAiUGCSqGSIb3DQEHATAeBglghkgBZQMEAS4wEQQMiWM0crDA8bZ0Fj5HAgEQgIIB9tIM8xbxhkhniJQifurRLHIaOcnos9gvgBj8Lkc34eBUtUg4oYzjusGP_FZiHfpBoRnM409FFCs-Y5xxomumnGRQ8UIZJe7CPSmaCcSklAdwHN3Ja-fHfALctP68U3mGGeR11P8vlESoS7HZhZEmVhFb5m7jPk1ENZKLlq9CnVTz2n2bI3HPYCfDIo8ywd8RWo5gw6EwYnZ46l1ExES8AHMqlklpgzdupLhL6HxTzBEHd_qY4B2FjK4m4MpvdaX5imEnOD6YCwZMNC_qLzUdPxumcXl4oTRDyJp1kg_6PcGi98hE10InOp4aHJeIjl5FPJnuEFEBcv28asOq4vngp6N5d-4y1fpeoxcAe4d857TeLtQnGD6LlCAL_3Mh0lQ8iYLLfjw4b_V2jHEuSxKH3BalLPMS5KF2nGqqrzfbDs8x9jW-OlkeL_x68Zidv0rTHqTd0ww3kvSINDjgw6owPMLK9KxEWP5T-4dnu_kn-NXR1W0s52kUHFD1U9AWs5Zzh4RTVOXsdjGQciWzOO-XieLn8Aju-0P8mVEGjUDOjLDZtuPYM-Ep_Z9lp8TNDCi2lB5IQFbKVpBFaAItDDyVUvKNPfBovmKlnPRFLZYjLbgocnU3tcy0YB8RZegC56Gxn7g_WkphptCVPm7IDeCE9l5gcFafuY0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Non-symmetric systems and the generalized minimal residual method (GMRES) [(Y. Saad, M. Schultz, 1986)](https://epubs.siam.org/doi/10.1137/0907058)\n",
    "\n",
    "Before we discussed symmetric positive definite systems. What happens if $A$ is non-symmetric?\n",
    "\n",
    "We can still orthogonalize the Krylov subspace using Arnoldi process, and get\n",
    "\n",
    "$$A Q_j = Q_j H_j + h_{j,j-1}q_j e^{\\top}_{j-1}.$$\n",
    "\n",
    "Let us rewrite the latter expression as\n",
    "\n",
    "$$ A Q_j = Q_j H_j + h_{j,j-1}q_j e^{\\top}_{j-1} = Q_{j+1} \\widetilde H_j, \\quad \\widetilde H_j = \n",
    "\\begin{bmatrix} h_{0,0} & h_{0,1} & \\dots & h_{0,j-2} & h_{0,j-1} \\\\ h_{1,0} & h_{1,1} & \\dots & h_{1,j-2} & h_{1,j-1} \\\\ 0& h_{2,2} &  \\dots & h_{2,j-2} & h_{2,j-1} \\\\\n",
    "0& 0 & \\ddots & \\vdots & \\vdots  \\\\\n",
    "0& 0 &  & h_{j,j-1} & h_{j-1,j-1} \\\\ 0& 0 & \\dots & 0 & h_{j,j-1}\\end{bmatrix}$$\n",
    "\n",
    "Then, if we need to minimize the residual over the Krylov subspace, we have\n",
    "\n",
    "$$x_j = x_0 + Q_j c_j $$\n",
    "\n",
    "and $x_j$ has to be selected as \n",
    "\n",
    "$$ \\Vert A x_j - f \\Vert_2 =  \\Vert A Q_j c_j - r_0 \\Vert_2 \\rightarrow \\min_{c_j}.$$\n",
    "\n",
    "Using the Arnoldi recursion, we have\n",
    "\n",
    "$$ \\Vert Q_{j+1} \\widetilde H_j c_j -  r_0 \\Vert_2 \\rightarrow \\min_{c_j}.$$\n",
    "\n",
    "Using the orthogonal invariance under multiplication by unitary matrix, we get\n",
    "\n",
    "$$ \\Vert \\widetilde H_j c_j - \\gamma e_0 \\Vert_2 \\rightarrow \\min_{c_j},$$\n",
    "\n",
    "where we have used that $Q^*_{j+1} r_0 = \\gamma e_0.$\n",
    "\n",
    "- This is just a linear least squares with $(j+1)$ equations and $j$ unknowns.\n",
    "\n",
    "- The matrix is also upper Hessenberg, thus its QR factorization can be computed in a very cheap way.\n",
    "\n",
    "- This allows the computation of $c_j$. This method is called **GMRES** (generalized minimal residual)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary of the GMRES\n",
    "\n",
    "- Minimizes the residual directly\n",
    "- No normal equations\n",
    "- Memory grows with the number of iterations as $\\mathcal{O}(j^2)$, so **restarts** typically implemented (just start GMRES from the new initial guess)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Figure size 720x360 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAFgCAYAAACmDI9oAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3Xd4FWX2wPHvm14ICUlogUAIRZqUSwJSBIREdBUboO5a1p+FoKuyskiIvWAJurpWDFjWvkAsKwpqAiJ1aaFLTyihQ0gIJf39/TE3EEICSci9c8v5PM99bjIzd+ZEyeHwzpn3VVprhBBCCCGEEAYPswMQQgghhBDCkUiBLIQQQgghRAVSIAshhBBCCFGBFMhCCCGEEEJUIAWyEEIIIYQQFXiZHYAjU0pZgBggBIgFErXWmeZGJYQQ7kPysBDCDFIgV0MpFQLEaK2nWr+PA9KAtqYGJoQQbkLysBDCLMoZ50FWSo0EcrXW6VXsKx9tyASigcyqjqvBNSzATK11W+v3IcAxoJHWOvdS4hdCCGcneVgI4cqcbgTZOoIwDRhVxb5oIFlrHV9h20ylVGZtb8lprTOUUvEVNsVg/GUgSVkI4dYkDwshXJ3TPKSnlIpWSqVgjEbkVHNYApBSaVsKkFyXa1ZK5gnAA3U5jxBCuALJw0IId+GsLRY7gITKt+ys2+MrJtTyW3Jaa2X9fjQX7l9Lq+K8o4EcrXVqff0MQgjhzCQPCyFcmdO1WFTHmoDPG9XQWucqpVBKRWutM8sf9qjFeeOoZf9ceHi4joqKqs1lhBBuYNWqVUe01o3NjsNWJA8LIRxdTfOwyxTIQCgYibia/dEYD4zUmPUBkRytdYb1+5E1Gb2Iiopi5cqVtbmUEMINKKV2mR2DjUkeFkI4tJrmYVcqkEPq82TWB03mAiFKqfLNmUCVidl6+280QKtWreozFCGEcBaSh4UQLsFpHtKzN+ttwEZaa1XhVW3PnNZ6qtY6Rmsd07ixy95BFUIIu5E8LIQwi8sVyNYeOCGEECaRPCyEcHau1GJR3tcWCpzpf6uQqG2+NKlSajgwvF27dra+lBBCOCLJw0KYpLi4mOzsbAoKCswOxSH4+fnRsmVLvL296/R5lymQrU9JZ3J+D1woxsTyNk/MWutZwKyYmBiZp1MI4XYkDwthnuzsbIKCgoiKiqJCz75b0lpz9OhRsrOzadOmTZ3O4WotFukYKy1VZLFuF0IIYXuSh4UwQUFBAWFhYW5fHAMopQgLC7uk0XRnLZBDqfpp6UTOX/o0wbrd5pRSw5VSU/Py8uxxOSGEMJPkYSEcjBTHZ13qfwunKZCVUiFKqWSl1EyMpJyslEpRSo0sP8Y692ai9biRSqkJQLI9butZrz9Laz06ODjYHpcTQgi7kjwshHAXTtODXJ50a3BcBpBh+4iEEMK9SB4WQrgLpxlBFkIIIYQQzik19aILYJ6RmZlJr169iI+PJz4+ntzc6hbntB0pkOtRbXvfSkqKSfn+CbIP7bRtYEII4Sbq0oOcsX0h36S9S2lZqQ0jE8K9TZ8+vcbH5ubmEhcXR1paGmlpaYSE2H9qdadpsXAGtZ1e6Jf/fcm7ebOYMvsHOhR50d67LR2bD6B/95G0adZSmu2FEKKW6jLN24dzk1jol8cbX3zKoDZxxDSLpUfjHkQFR+GhZBxJOJ/nZ23kj33H6/WcnSMa8uzwLtXuz8jIICUlhZycHIYPH86sWbMIDQ0lISGB6dOnk56ezqhRo5g2bVqNCt7U1FQyMzMJDQ0lJSWlPn+UGpEC2UTxvf9McXEhSzJ/YJPazQ8eW/nh4FY8f/mIVkWKpsVhBHl2xD/oKhqF96BJQ1+aBPnSpKEfTYJ8CW/gi5+3p9k/hhBCOLUJN35Ohy/u5ID/fubvmMOszB8BaOjTkO6Nu9OjSQ96NO5B1/CuBHgHmBytEI5rxowZZGVlkZiYSFJSEhaLBQCLxUJGRgYzZ84EjBaK6louJkyYQHR0NMnJyYwcOZJRo0aRnp5OXFyc3X4OkALZVD4+vtx0VQI3XZUAwL5DO1iweiYbDyxlh9rDmoAjFHgshpLFNNlTRrOCQDjVmr0ne7KrsAsaLwJ9PAlt4ENooC9hgT6EBvoQ4u9NkJ83QX5eFV7eld698PWS4loIIaKateH2+35lxWcTeSnrc+Z7R/JJ85tQIafZmbeVhXsXAuClvOgU1glLEws9m/bE0sRCI79GJkcvxPkuNNJrS7feeishISEkJiaSmJhIbm4uKSkpREdHn3NcdHQ0EyZMqPY8ISEhjBxpTI4TGxtLZqZdJsE5hxTI9ehSlziNaNKW24dNPPN9SXEBW3bMZk1WGquObCDD4xhHgzYDm2lRWsZlZYE0UVHg0Yu9pX04eLyMTfuPk3e6mFNFF++l8/HyoGGlojnI15sAX08CfDwJ9PHCv+K7ryf+3l7GPuvXgb6eZ4/x9sTDQ9pChBDmqWsebtYokOFj3+HAqmvoO/tBrtgzjfE7RrOwLIGwoBLaRh4lMHgPx4u38tXmr/j0j08BiA6OxtLUgqWJBUtTCxGBEdIeJ9xedHQ0M2fOJD09nZSUFJKTk8/Zf7ER5IojxitWrCApKcnmMVemtNZ2v6iri4mJ0StXrqz38+qyMvZkL2HVth9YfSiDjNMH2WUdBPYr01yu/IgN6UDvNsPo3OEWCrUv+QUlHC8oJr+gxPoyvj5RWPX2/AKjuDZeJRQUl9UqRj9vj9oX1j7VF+QB3sbXPl7SByicn1Jqlda68ipzwgYuKQ/n7kHPvAe1dyWbo+7iQ/97WL4rn905pwDw9ymjfWQuoWHZnPbczq6Tf3Cy+AQATQOaYmlqoXez3vRp3oeWDeR5EmEfmzZtolOnTqZdv7wHOSUlhcmTJ5OWlgZAcnIyFouF+Ph4gCpHlCvLzMw8MwJtsVjOK7Brqqr/JjXNw1Ig24CtCuSqHDm8idVbviNj31JWndzDZlWCVgr/Mo3FI5DeoV3o0+46OrYfjqeXT63PX1qmOV1cyqnCEk4VlXKyqITTFQpoY1spp4tKOFlYyuniUk4Wnj2m/PgzxxSVWr8voTZ/9Lw9Ff7engT6el2gsPbE38fL+u5JgI+1EK/8OW9PfL098PMy3n08PWTkW9iFFMj2c8l5uKQIfn0KlqdAZB8Y9W8OEsrKncdYuSuHlTuPsXFfHmUaoIzGoccIC9+Lh38WOaWbOVl6DICIwAh6N+9N72bGq2lg03r5+YSozOwC2RFJgexg7FkgV5aXt5uVG79mefZClp/YzXZP4/9vUJkmxjOYPo270bv9jbSLvhrlYd6orNaawpIyo4guLDmnsD5Zofg+VbmwLizhVIWCvbxQP1uEl1JUWrtRbzDaTXy9PPDz9sS30tdnt3ni5228+3pXtd8D3wt8xs/7/P3enjIy7k6kQLafesvD61Phh0fB2x9GfAhtrzqz60RhCWt257J+bx7bDuWz7eAJth86weniEjx8DuMZuIPA4Czw206pMkafm/m34oqIPgyM7Ets01hC/Ow/fZVwTVIgn08KZAdjZoFc2ZHDm1j+x39Yvm8py07vJ9vakhFapunjHUbvpr3o03EELVv0NbVgrk/FpWXnjHCfKrR+XWx8fbq4lMKSUgqKyygsKaWwuIwC63thSRmFxaUUlpRRYH2veOyZz1TYfym/Qp4e6ryCvHIhfqYgr/ReoyK+ms/IqLk5pEC2n3rNw4e3wIy7jfernoQr/wHV5MuyMs3e3NNsPZjPtkMn2Hown60Hj7MjdxslvlvxCtyBp38WyrMIgBDPKDoE96RfRB/io68gMiRUWjJEnUiBfD4pkB1EhYdDHti2bZvZ4VRp797lLN+cyvIDK1hWeJjDnkYijiiF/gEtuLJ1HH0uv5uABk1MjtQ5aK0pKi0vrKsuqgurKKoLis8W5AUXKMiN7ZXOYy3oi0sv7Xe3fNT8bFFdg0K7ihF23+pGy6v5jJeHctsCQApk27NZHi48AT/+HdbPhC43w00fgLdfjT9esXDefDCXVfvXsi1/DTmlf6D8dqE8StDaA1XYikaqG9GBMVwe3ok2jRsQFR5IVFggjQK83fZ3R1ycFMjnkwLZwTjSCPKF6LIysnbNZ/mW71hyKINlpXmc8lB4a41F+XNleA8GdLyV6DZDXWZ02ZWUlukKhfNFCu0Ko+QFVXymsIaj5eWfv5S04aGovgWlmlHzyqPlPp4eZ/rHywv0swW/8fny730qfW9mgS4Fsv3YJA9rDUvehrRnoFVfuP0rCAi9pFOWlWkyjx5jbuYKlh9Yxtb8leSWZhn7SoIoPdGBkhOXUXKyHcG+wbQOC6BVaMCZ91ahgbQOC6BZQz+5K+TmpEA+nxTIDsZZCuTKigtPkrHxSxZlzmHh8R3ssPYvR5TCgIAWDJDRZYExal5cqisU2lUX1ee1qlQzwl5VIW6rUXMApbhoEe1befs5BfnZY6/p0oyo8MBaXFsKZHuxaR7e8A18NwZCWsMdMyG0Tb2e/sjpIyzeu5jf9yxgyb6lnCzJR+FBI8/2+BZ15mReOw4eDqOkwuMWPl4etGzkT6vQACIbGcXzyF4taRRY+4ezhXOSAvl8UiA7GGctkCvbv28VCzd8zqIDK6oYXe7JkMvvIrJl/2p78YSob+Wj5kUlxgh4UUmFgvzM9xW2F5cZLTDFpdb3sjMtMdV99vxjK1+vjFJj6gI+vieGIR1rPiuBFMj2Y/M8vGsJfP1nUB5w+5fQup9NLlNSVsL6I+tZmL2QRXsXsSlnEwDhfuH0CO9Dm8BYGpR15UBuGXtyTrHb+sovKGHZE0Np2rDmbSDCuUmBfD4pkB2MqxTIFVUcXV50fMeZ2TE6FJUwVAUypGF7LmvaE9W4IzTuaIyoeHqbHLUQtlFSahTP3p61m4lECmT7sUseProDvroVju2CPgnQ/+/QoLFNL1k+urxo7yIW71tMflE+vp6+9I3oy9BWQxnUchCN/BqRd6qYhv5e0rPsRtyhQE5MTCQzM5PY2NgLrsRX7lIKZFlJT9SIt28gfSyj6WMZzT+APdnL+G3jF8w9uo4PinOYcnoDLbat4aq1pxh66jQ9i8rwDGsLjS+D8MuM98YdIbwD1GE+ZiEciZenB14yRZ8Iawv3p8PPT8D/3oeVn8AVY6DfI+BvmyWow/3DubHdjdzY7kZKykrIOJjBvD3zmLt7LvP3zMdDedCraS+GRA5hSKshRDSIsEkcQtRWamrqmeWja/OZlJQU0tLSyMjIAGDmzJlnCuWLLThyKWQEuR45wywWtnD09FF+z/6deTt/ZemB5RSVFdPIw4dBBBCff5y+h3fjra3Nch7eRrHctAs07QrNukLTy20+6iKEI5ARZNszLQ8f3grzX4GN34JvMAxOhNgH7DYgoLVmU84m5u6ey7zd89ieux2ATqGdGNJqCNdEXUNUcJRdYhHmcPQR5FGjRjFz5swaH5+ZmUl6ejozZ84kLS2NyZMnY7FYiIuLIzU1lZycHEaPHn3Bc0iLhYNxxRaLmjpVfIpFexcxb888FuxZQH5xPsE+wcQ3ieFa/xb0OnkCz4Mb4eAGyN9/9oMNmp5bMDe7HMLbg4eneT+MEPVMCmT7MS0PH9xozHKxPR3C2sM1r0L7OLuHsev4Ln7b/Rvz9sxjzaE1aDSdQjtxTZtruCbqGhlZdkHnFINzJsKB9fV7gWaXw7WvVru7fKnpnJwchg8fzqxZswgNDSUhIYHp06czdepU4uLimDZtGiEhF18gJyEhgeTkZEaNGkVaWhqJiYncdtttWCwW0tPTSUtLu+gS1NJiIRxGgHcAV0ddzdVRV1NcWsySfUuYs3MOP+2eR2rJacL9wxnWYRjXXP0U3QNaog5tgAMbjIL5wAZY+j6UFRsn82kAzbtDRM+zr9BoYxoCIYRwRE27wB2psO1X+DkJvhwBHa+HP70ODZvbLYzWDVtzT9d7uKfrPRw6dYhfd/7KnJ1zeHPVm7y56k26N+7OtW2u5erWV9M4QO7gifoxY8YMsrKySExMJCkpCYvFAoDFYiEjI+PMCHJmZiapqalVnmPChAmkpqbSq1evcwrp3Nxc2/8AFcgIsg248whydU6XnGZh9kLmZM1hQfYCisqKiAiM4Jo213Bj2xuJDrH2EZUWw5GtsH8t7FsD+zKMfwWXFBj7/YKheQ+jWG5hMd6DI6VoFk5BRpDtxyHycEkR/O89mP8qePrCsEnQ8y5T81V2fjY/7/yZn7N+ZsuxLSgUsc1iuT76eq6OuppA75pPWygci9ktFuUjyCkpKWRmZpKYmEhubi4pKSlER0cTHx9PWlpajc6VkJBATk4OAOnp6SQlJQEQHR3NyJEjSU1NJTMz86IP6kmLhYNxiMTswE4UneC3Pb8xO2s2S/ctpVSXcnn45dzU7iaGRQ0j2Df43A+UFsPhzbBvNezNMN4Pbjw70hzYGCL7QGRv4715j1qtcCWEvUiBbD8OlYePbIdZj8KuxRA9GG54B0JamR0VmbmZ/LzzZ+ZkzWHn8Z34e/kT1yqOG9vdSGyzWDyUPIjqTBypQC5XsRWiYoF8sRHkcrm5uWdaLDIyMpg+fTrJyckkJCSQkJBwZoS6OtJiIZxKA58GDG87nOFth3Pk9BF+yvyJ77d/z4v/e5Hk5ckMaTWEG9vdSN/mffH08DSmi2tm7Uu23G2cpKTQKJL3ZUD2KtizDDb/aOzz8IaIHmeL5pa97XprUwghzhHeDv76I6z6xOhPnjIAbnjLWLLaRNEh0TzU4yEe7P4g646s47/b/8vPWT8zK3MWEYER3NDuBm6IvoHIhpGmximcz+TJk88UwxX7hOPj48+MKNdkmraKLBYL06dPZ9SoUURHR1+0OL5UMoJsAw41cuEkyp/A/n7798zOmk1eYR5NAppwS/tbGNF+BM0Cm138JCcOQ/YKo1jes9wonstbM4JbGcVyqysg6kpjJg1pyxB2JiPI9uOweTgnC765H/auhJ53wrWTwcdx2hoKSgqYt3se/93xX5buW4pGE9sslls73MrQVkPxlvntHZbZI8iOSFosHIzDJmYnUVRaxO/Zv/PNtm9YsncJHsqDQS0Hcdtlt3FFxBU1v+1XUmT0L+9ZdvZVPnNGQLhRLDfpbMyWEdbOKJod6C8q4XqkQLYfh87DpcVGX/LCfxr555pXoe0Qh/tH+4GTB/hhxw98u+1b9p7YS5hfmDFo0WEELRq0MDs8UYkUyOeTAtlBuOs8yLa0J38PqVtT+W7bdxwrPEZkUCSjOozilva3nN+rfDFaw7Es2LnYWCZ2zzLj+/I5mlHGQibNuxstGs17QPNu4BtU7z+XcE9SINueU+XhzN/h+wfh+F6jQL7un8ZMPQ6mtKyUJfuWMGPLDBbsXYDWmitbXsntl91O/xb9pVfZQUiBfD4pkB2MQ49cOKmi0iLSd6Uzfct0Mg5l4O/lz83tbubOzncSGXQJ/XElhcYtzyNb4dAfxswZ+9dC/j7rAcoYWY6wGLNmtLAY8zV7+dbLzyXcixTI9uM0ebik0FiBb94k48HjQYnGSnwO2sqw/8R+Urel8u22bzly+gjRwdHc3flurm97Pb6ekhfNJAXy+aRAdjBOk5id1JacLXz2x2fMzppNmS5jaKuh3N35bno06VF/FzlxyDrN3GrYu8roZz552Njn6WMUyS0s1sK5lyxqImpECmT7cbo8fHw/zJkAm36AZt2MmS4i6jGn1bPi0mJ+3vkzn/3xGZtzNhPqF8rtHW/ntstuI9Qv1Ozw3JIUyOeTAtnBOF1idlIHTx7k681fM2PrDPKL8unWuBujLx/NwJYDUfXdy6c15GWfLZb3ZhgFdFG+sd+ngdGSUT7K3KKXzM8sziMFsv04bR7+4weYPR5OHjFGkgdPBG9/s6Oqltaa5QeW8+nGT1m4dyG+nr7c1O4m7u16r6zWZ2dSIJ9PCmQH47SJ2UmdKj7F99u/57M/PmPvib10Cu1EQrcErmp1lW1748rK4Og2o1guL5wPrIfSImN/QPjZUebylQCDmtouHuHwpEC2H6fOw6ePwa9Pw+rPjQeIb3gHWvczO6qL2pG7g8//+Jz/7vgvaLih3Q3c3/V+mSbOTqRAPp8UyA7GqROzEysuK+anzJ+Ytm4au/N3075RexK6JRDfOt5+D5GUFBnLZpePMu/NgCNbzj4I2LCFtVjuYbRpNO4IIa3BQx5ycQdSINuPS+ThHb8ZC4zk7obYB+CqJyDA8dsXDpw8wMcbPuabrd9Qqku5Lvo6Erol0Kqh+YujuDJXK5BHjRpFbm4uFovlzFzKiYmJZGZmEhsbe2Ye5aq2lZOFQoQAvD28uandTVwffT0/7/yZqeumMv738XQO68xjvR7jiuZX2D4IL5+zbRax1m2FJ4yR5X2rra+Ms4uaAHgHGLNnNOlkFMzl78GRUjgL4c7aXgUPLoV5L8KyFFg3Hfr+Da54CPwamh1dtZoFNuOJPk/wwOUP8O+N/2bGlhnMzpzNiA4jGNN9DOH+4WaHKEyQmprKyJEja3Ts1KlTSUpKwmKxMGrUKDIyMs7smzlz5pmiODc397xt0dH1MxOMjCDbgEuMXLiA0rJSZmfN5p3V77D/5H76RfTj75a/0ynMAf6FXXDcWD770KZz38vnaQbwDjRm0GjcEZp0hMadjHfpbXZaMoJsPy6Xhw9uhN9eNv5x7d8I+o+F3qOdYu72I6ePkLI2hdStqXh7evPXLn/lr53/SgOfBmaH5lIcfQR51KhRzJw5s0bHVix0p06dCnBmNDkuLo7U1FRycnKq3DZ69Ogz55ERZCGq4OnhyfC2w7k66mqmb57O1PVTufXHW7k++noe6/UYTQKamBecX0NjZb/I3uduP30MDm85t3DeMRfWfnX2GJ8GRuEcZl3gJLydMQId2ha8/ez7cwgh7KNpF7j9S+Mu1G8vQ/pzsPQ9uPIf0Ov/HPp3P9w/nCeveJK7Ot/FO6vf4YO1HzB983Qe7vkwI9qPwFNmAKp3ycuT2ZyzuV7P2TG0I4m9E6vdn5GRQUpKCjk5OQwfPpxZs2YRGhpKQkIC06dPJz09nVGjRjFt2jRCQkIueK2Ko8CrVq06c47QUKPFKCQkhBUrVgBUua0+SIEsXJ6vpy93d7mbm9vfzMcbPuazjZ8xb/c8EroncFenuxxr6VT/RsYKf60qtYOcyjl3pPnwZti5ENb9p8JBChq1Norlyq/AMLv+GEIIG4noCXfMhN3L4LdJ8PNEWPw2DBwPPe8y2rwcVKuGrXht0Gvc0+UeXl/5Oi/+70VSt6aS1CeJnk16mh2eqAczZswgKyuLxMTEMy0SABaLhYyMjDMjyJmZmaSmplZ5jop9xBkZGeTk5GCxWEhJSTnv2PIWC1uQAlm4jSCfIMZaxnJLu1uYvGIyb656k++2fcdTVzxFn+Z9zA7vwgJCjafYKz/JXngCcnbAkW3W11bjPWsBlBScPc4/1CiUG1csnNtbHxCU0RshnE6rPvDXWcbv+rxJ8NM4WPwvGDQRut0Gno7713uX8C58POxjftn5C6+vfJ2759zN8OjhPB77OI38Gpkdnku40EivLd16662EhISQmJhIYmIiubm5pKSknNcXHB0dfd4DdZWVf7a8qG7bti2ZmZlYLBZyc3MJCwsjLCzsvG31xXF/g5xQhSVOzQ5FXEBkw0jeGfoOC7MX8uryV7n/1/sZ0X4E42LG0dDHcR98qZJvA2Np7Obdz91eVgp5eyoUzdbCecscyPjs7HFe/tCsq/UcPYzZNRp3dNhVvIS4GLfLw20Gwr2/wPa5xsN8/30IFr0Bg5Og6wiHfV5BKcU1ba5hYMuBTFs/jX9v/DeL9i4iqU8S10RdU/9z2Qu7io6OZubMmaSnp5OSknJmFopyNRlBTkxMPGfUOC4ujunTpzNy5EjS0tJISEgAqHJbfZCH9GzA5R4OcWEFJQW8v/Z9Pt34KeF+4Tzd92kGRw42OyzbOpVztnA+tAn2r4H9684ueuLlZ6zk1aLX2UVPQqMd9i9aZyIP6dmPW+ZhrWHzT0aP8qGNED0Yhr8FjaJMDuzith7byrOLn2XD0Q0MbjmYp654iqaBMm98bZj9kF55D3JKSgqTJ08mLS0NgOTkZCwWC/Hx8QBVjihXNnnyZFJSUs70KiclJTFy5MhzZqqoPPVbxW3lZB5kB+OWidnJbTyykaeXPM22Y9sY1WEUj8c+jr+X465eVe/KyiAn0yiWy5fX3r8Wik8Z+/2CrctqWwvmCAs0bG5uzE5ICmT7ces8XFYGqz6BtGdBl8LQZ4wZLxy8naq0rJQvNn3Bu6vfxdfLlxf6vcCQVkPMDstpmF0gOyIpkB2MWydmJ1ZcWsw7a97hkw2fEB0czeSBk7ks9DKzwzJPaYnxMOA+60qBezOMqaZ0qbE/qLm1WO559t3/wk8muzspkO1H8jCQlw0/PgbbfoWWveHGd40ZcBzczrydTFgwgU05m7jtstsYHzMePy/HnaXDUUiBfD4pkB2MJGbntnTfUp5c9CS5hbkkxiZyW8fbzA7JcRSfNhY9KS+Y964yHhIsF9bOOtJsLZibdjH6pAUgBbI9SR620hrWz4Q5iVB0AgZNgP5/d/jnDIpLi3l79dv8e+O/aRvcljevepM2wW3MDsuhSYF8PimQHYwkZueXU5DDk4ueZNHeRYzqMIqk3kmONR2cIzl9zNqWkXG2aD5xwLpTQWgbaHa59dXNeA9q7pY9zVIg24/k4UpOHIY5j8PG74xl7q9NhqgBZkd1UUv2LmHiwomUlJXw2qDX6N+iv9khOaxNmzbRsWNHecDRSmvN5s2bpUB2JJKYXUNpWSnvrH6HjzZ8RK+mvXhj8BuE+oWaHZZzOL7P6GE+sB4OrIMDG+BY1tn9/qFni+Ymna0rBXZ0ilXBLoUUyPYjebgam3405k7O2wOdhkPc8xDW1uyoLmjfiX08Mu8RtuduZ1yvcdzd+W4pAquQlZVFUFAQYWFhbv/fR2vN0aNHyc/Pp02bc+88SIFsIknMruWnzJ94dsmzhPmFMSV+CtHB9bPOu9spOG70MB9YDwfXW9//gNJC6wHKaNFo3h2adzNGm5t3N+aAdhFSINuP5OELKD4NS9+FhW8LmucrAAAgAElEQVQav39XPAiDEsE3yOzIqnWq+BRPLX6KtF1p3NzuZp7p+wxeHjJTbUXFxcVkZ2dTUFBw8YPdgJ+fHy1btsTb+9y7v1Igm0gSs+vZcGQDf5v7N7TWTImfQpewLmaH5BpKS+DYTji8ySie968zRp6PZ589JjjSmJ85oqfxat7DaYtmKZDtR/JwDZw4ZMydnPGZ0fZ0zSvQ+SaHbX8q02W8v+Z9UtalcFXkVbw26DV8PX3NDks4GSmQTSSJ2TXtOr6L0b+OJq8oj3eHvEtMM6lzbObkUWtrxjrYt8aYfi4n8+z+RlFni+WInsZIsxPMoCEFsv1IHq6F7JXGbBcH1kH7YfCn14xl6x3UV5u+4pXlrxDbLJZ3h7xLgHeA2SEJJyIFsokkMbuugycP8kDaAxw4eYCp8VPp0aSH2SG5j9PHjNHlfavPvnJ3n90f2vbsKHPLWKNo9nasqaGkQLYfycO1VFoCy1Ng3kuAts6dnAAeHmZHVqUfM3/kqUVP0atpL94b+p5MAydqTArkeqCUsgChQAgQDyRrrTMv/ClJzK7u8KnD3PPzPRwrOMaHwz6kc1hns0NyX6dyzi2Y962G43uNfR7exkOALWOtrxhj5NnE28dSINee5GE7y90DP40z5k6OuhJueh9CWpkdVZVm7ZjFk4ueZECLAbx11Vsy05CoESmQ64FS6hjQRmudq5QaDSRorXtd7HOSmF3f/hP7+evPf6WgpIAvr/uSyKBIs0MS5fIPGLeM9660vmdA8UljX0DY2WK5aVcI7wAhrcHTPg/7SIFce5KHTaC10Zf8yxOAgmGToMeddvs9qY0ZW2bw4v9e5E9t/sSrV77q9rM3iIuraR52vD/tNaCUGgnkaq3Tq9hnAWKATCAayKzquBpqo7XOtX6dU8dzCBfUvEFzpsZP5Y7Zd/DovEf5/NrPaeAjC2I4hKBm0Ol64wXWFQE3QfYKyF5lvG/9+ezxnj5Ge0Z4e6NgDu9g/bq9Qz/VbzbJwy5MKej1V4geBN8/BLPGwuK34Mp/QLfbHGqRkVsvu5W8wjzeXv027ULa8UC3B8wOSbgIpxtBVkrFATOBUZUTrlIqGkjRWsdX2DYTSKzJLbmLXDcN49beRZO8jFy4j//t/x9j0sacucXn6eFpdkiiJgry4PBWOFL+2ma852SeXUobICiiisK5AzSMqFOrhquMIEsediNlZbDlJ/h9svEQX0grGPAY9LgDvBxjBgmtNRMXTmR21mz+NfhfDG091OyQhANzuRFka9JNBFZR/ShCApBSaVsKkAyMuoTrjgTSLmEERLioK5pfwcTeE3lp2UtMXT+VB7s/aHZIoib8giEy1nhVVFJkTDt3ZCsc2XK2cF43HQqPnz3Op4ExZ/PVk6DNlXYN3UySh92Qh4exoEjH62HrL7BgsjHjxYLX4Ya3oV2c2RGilOL5fs+zJ38PSYuSSG2USquGjtk3LZyH040gAyildmD0oVUeudgBxFccpVBKhQDHtNbK+v1o4ELLBlWZgK2fG1VxVKQ6MnLhfhIXJPLLzl/44k9f0DW8q9nhiPqmNZw4eP6I89BnjFkzashVRpBB8rDb0hp2zDMK5Oteh6aOMyf8gZMHuOWHW2gT3IZPr/lUFhIRVXLph/SqSszlCRhoVKFfrXyfBtrW5vZe+YiF1npypfNf9DySmN3P8aLjjPhhBH6efky/frrMyymq5OoFsuRhYbY5WXOYsGACf+vxN8Z0H2N2OMIB1TQPO+YEh3UTClA5KVdQ2/WBo4GwSt/nXmoPnXBNDX0a8lL/l9h5fCdvrnrT7HCEMIvkYWGqa9tcy5/a/IkP1n7AxqMbzQ5HODFXKpDrdRkt66jICqXUaOttvSSg2s5/63ErlVIrDx8+XJ+hCCfRu3lv7ux0J9O3TGfDkQ1mhyOEGSQPC9M90ecJQv1CmbR0EmW6zOxwhJNypQK53mmtU7XWU62vUVrrjAscO1VrHaO1jmncuLE9wxQO5G89/kaYfxgv/e8lScxC1APJw6K2gn2DGRczjg1HN/Ddtu/MDkc4KZcrkK09akKYooFPAx7r9Rgbjm4gbVea2eEIYQrJw8Js17W5jp5NevLemvc4XXLa7HCEE3KlArm8Jy204sYKidrmPWtKqeFKqal5eXm2vpRwYNe1uY52Ie14d/W7lJSVmB2OEPYkeVg4BKUUYy1jOXz6MP/Z/B+zwxFOyGUKZOtDIZmc3wMXip0e6tBaz9Jajw4ODrb1pYQD8/Tw5OEeD7Pz+E5+2fmL2eEIYTeSh4Uj6dW0F/1b9OeTDZ/IKLKoNZcpkK3SMZY3rchi3S6E3VzV6ipjLs6Nn+KMUykKcQkkDwuHcX/X+zlWeIxZO2aZHYpwMs5aIIdS9dPSiZy/UlOCdbvNya09Uc5DefDXzn9lU84mlh9YbnY4QtiC5GHh8Ho17UXXsK58uvFTeXBa1IrTFMhKqRClVLJSaiZGUk5WSqUopUaWH2O9vZdoPW6kUmoCkGyvOTPl1p6o6Pq21xPiG8LMrTPNDkWIeiF5WDgbpRR3db6L3fm7ZbBC1IrTrMNYnnRrcFwGUO00QELYi6+nL9dHX8/0LdPJLcglxE8e7BfOTfKwcEZDWw8laFkQ3277liuaX2F2OMJJOM0IshDO6Ob2N1NcVsyPmT+aHYoQQrglX09frmtzHXN3zSWvUFpvRM1IgVyPpPdNVNahUQc6hXZiTtYcs0MRwi1IHhZVuandTRSVFTFv9zyzQxFOQgrkeiS9b6IqV0ddzboj69h/Yr/ZoQjh8iQPi6p0DutMiwYtZAEnUWNSIAthY/Gt4wFI3y2zXAkhhBmUUsS3jmfp/qUcLzpudjjCCUiBLISNtW7Ymg6NOsjIhRBCmCi+dTwlZSX8vud3s0MRTkAK5HokvW+iOoMjB7P28Fp5QEQIG5M8LKrTNbwrYX5hLMxeaHYowglIgVyPpPdNVOfKFldSpstYun+p2aEI4dIkD4vqeCgP+rfoz5L9SygtKzU7HOHgpEAWwg66hnclyCeIRdmLzA5FCCHc1oAWA8grzGPD0Q1mhyIcnBTIQtiBl4cX/SL6sXjfYrTWZocjhBBuqV9EPzyUB4v2ymCFuDApkIWwkwEtBnDk9BG2HNtidihCCOGWgn2DuTz8chbvXWx2KMLBSYFcj+ThEHEh/SP6A8jIhRA2JHlYXEz/Fv3ZcGQDOQU5ZociHJgUyPVIHg4RF9I4oDEdQztKgSyEDUkeFhdzZYsr0WiW7FtidijCgUmBLIQd9Y/oz5pDa8gvyjc7FCGEcEudwzrTyLeRDFaIC5ICWQg7GthyIKW6lAXZC8wORQgh3JKH8mBAiwEszF5IcWmx2eEIByUFshB21KNJD5r4N+HnnT+bHYoQQritYVHDOF50XOamF9WSAlkIO/JQHgxrM4zFexdzvOi42eEIIYRb6hfRjyCfIH7OksEKUTUpkOuRPD0tauLaqGspLiuWxCyEDUgeFjXh7elNfOt45u6eK4MVokpSINcjeXpa1ETX8K50b9ydD9Z+wKniU2aHI4RLkTwsauovHf/C6ZLTTFs3zexQhAOSAlkIO1NKMT5mPIdPH+aNVW/IynpCCGGCy0Iv44a2N/Dlpi9Zf3i92eEIByMFshAm6NGkB3d2upPpW6Yz/vfx7MnfY3ZIQgjhdsZaxtIkoAn3/XofX2/+mqLSIrNDEg5CCmQhTDIhdgKP9HyE+XvmM/y74Tyz+BnyCqVvUggh7KVxQGM+v/ZzOod15uVlL3Ptt9fy/fbvzQ5LOAApkIUwiVKK0d1GM2fEHP7c8c/MypzFHbPvIDs/2+zQhBDCbTQOaMwnwz5havxUmgc25+nFT/Pyspel/c3NSYEshMmaBDQhsXciH139ETkFOUxYMIHiMpm8Xggh7EUpRd+Ivnx6zafc1fkuvt78NTO3zjQ7LGEiKZCFcBCWphae7fss64+sZ+q6qWaHI4QQbsfTw5PxMePp27wvr614jay8LLNDEiaRArkeyfyb4lINixrGDW1vYOq6qaw5tMbscIRwOpKHxaXyUB5MGjAJXy9fkhYmyR09NyUFcj2S+TdFfUjqnUTzwOb8/be/S5EsRC1JHhb1oUlAE57r+xwbj25k/PzxMme9G5ICWQgH08CnAe8PfZ8A7wDu/eVevtv2ndkhCSGE24lrHUdibCLzs+dzx+w7ZDpONyMFshAOKDokmq+v+5peTXvxzJJn+GbrN2aHJIQQbufOzncyJW4Kh04d4p4593D41GGzQxJ2IgWyEA4q2DeYKXFT6B/Rn0nLJkm7hRBCmKBfRD8+HvYx+cX5/P23v1NcKj3J7kAKZCEcmJeHF5MHTaZZQDMmLpzIiaITZockhBBu57LQy5jUfxLrjqxjytopZocj7EAKZCEcXEOfhrxy5SvsP7mfJxY9IaMXQghhgqujrubmdjfz0YaPSNuVZnY4wsakQBbCCfRo0oPE2ER+2/Mbj/72KAUlBWaHJIQQbmdi74l0b9yd8b+PZ9aOWWaHI2xICmQhnMRfOv2F5/o+x+K9i3l47sMykiyEEHYW4B3AB3EfENs0licXPcl/t//X7JCEjUiBLIQTGdFhBC/2f5FlB5bxr4x/mR2OEEK4nQDvAN6Le4/ezXvzwtIX2JKzxeyQhA1IgSyEk7mx3Y3cftntfPbHZ8zfM9/scIQQwu34evoyeeBkQnxDGP/7eE4WnzQ7JFHPpECuR7LEqbCX8bHj6RjakacWP8W+E/vMDkcIhyF5WNhLqF8orw58ld35u3lh6Qtorc0OSdQjKZDrkSxxKuzF19OX1we9TmlZKff9cp8UyUJYSR4W9hTbLJaHuj/E7KzZvLTsJcp0mdkhiXoiBbIQTqp1w9ZMu3oaeUV5JKQlcLzouNkhCSGE2xndbTT/1+X/mL5lOh+s/cDscEQ9kQJZCCfWNbwrb1/1Ntn52UxcMFFu8QkhhJ0ppXis12Pc2PZGpqydwtzdc80OSdQDKZCFcHIxzWJ4PPZxFu5dyMytM80ORwgh3I5Simf6PkPnsM48v+R5jpw+YnZI4hJJgSyEC/hzxz/Tt3lfXl/5OhuPbDQ7HCGEcDs+nj68POBlThafZOKCiTJXvZOTAlkIF6CU4uUrXybUL5SH5j7Envw9ZockhBBup21IW57p+wzLDizjmSXPSNubE5MCWQgXEe4fzvtx71NSVsJD6Q9xrOCY2SEJIYTbubHdjTzc42F+zPyRt1e/bXY4oo6kQBbChUQHR/Pu0HfZe2Ivzy99XkYvhBDCBKO7jWZE+xF8uP5DFmYvNDscUQdSIAvhYno26ckjPR9h7u65/Jj5o9nhCCGE21FK8USfJ2gX0o7nljxHbkGu2SGJWpICWQgXdHfnu7E0sfDC0hfkoT0hhDCBj6cPkwZM4ljhMf7x+z/koT0nIwVyDSmlUsyOQYia8vTw5I3BbxDmH8Yj8x7hwMkDZockxCWTPCycTZewLjzf73mWH1jOpGWTpO3NiUiBXANKqThgtNlxCFEbYf5hvDPkHU6VnOLReY9SWFpodkhC1JnkYeGshrcdzuhuo/l227d8uelLs8MRNeSUBbJSaqQ1WVa1z6KUGq2Uiit/v8RrhQA5gDQQCafTvlF7kq9MZlPOJqasmWJ2OMKFSB4Woub+1uNvDI4czJur3iQzL9PscEQNOF2BbE2006rZFw0ka62naq3TtdZTgQTr9rqK01pnXMLnhTDVoMhB3NzuZj7Z+AmrDq4yOxzhAiQPC1E7HsqDZ/s+i7+3P08sfELu6DkBpymQlVLR1v6zaIyRhKokAJV71FKA5DpeMw5Ir8tnhXAkj8c+TmRQJOPmj2P/if1mhyOclORhIeou3D+c5/s9z8ajG3lh6QvSj+zgvMwOoKa01pkYiRelVGI1h43k/MS8Ekgr/0YpNRpoe4FLpWmt08tv6Wmt5ZaecHpBPkG8PeRt7vjpDsb+NpZPr/0Ufy9/s8MSTkbysBCXZmiroTzU/SHeX/s+lzW6jLu73G12SKIa9V4gK6VKtdae9X3eGlw3hCpGNbTWuUoplFLRWutM6+2+mogDQpVSMdbvQ6xJPd36l4QQTiU6OJrkgck8PPdhJv1vEi8NeMnskIQNmZGLJQ8LcXEJ3RPYemwr/1z1T7qGd8XS1GJ2SKIKtmixUNXuUGqIUup+pVSPCtuGKqWG1MN1Q8FIxNXsr1X/m9Y61dpDN7U8mVu/lqQsnNbAlgO5//L7+WHHD8zfM9/scIRtmZGLJQ8LcREeyoNJAybRPLA5Ty9+mtMlp80OSVTBFgVylU01SqlXgTFAO2CyUuoVAK31XGBmPVw3pB7OcR6lVIhSaoL16wnVPWhifVJ7pVJq5eHDh20RihD14sHuD9K+UXueXfKszI/s2szIxZKHhaiBQO9AXuz/Irvzd/PyspelH9kB2fMhvR1a61u11hO11lcDM5RS4637jtkxjlrRWudqrSdrrZX1vcqRC+uoRozWOqZx48b2DlOIGvP29Oa1ga9RUFLAuPnjKCotMjskYV9Ol4slDwtXFNssloRuCXy//XtmbJlhdjiiErvOYqGUilJK3Q+gtV6ttX5dKfUA0Kger2GTEQwhXEnbkLa8NOAl1h9Zz0vLXpLRCzdj61wseViImnmox0Nc2eJKXl3+KhkHZSZDR1KrAvlS+ta01tMwnlpWVWy/tTZxVKN8RCG04sYKidrmPWtKqeFKqal5eXm2vpQQlyyudRwPXP4A3277lh8zfzQ7HFELDpyLJQ8LUQseyoNXB75KRIMIHv/9cfIK5c+to6hxgVwffWta67nWJHze9prGcYFz52Ik38ojF6FArj0e6tBaz9Jajw4ODrb1pYSoF3/r8TcsTSy8suwV6Ud2Eo6ciyUPC1F7DX0a8tqg18gpyJE7eg6kNtO87dBaTyz/RinVUyk1Xmv9Ouf2rSmlVBS1f1gjV2u9s5afqSwdiAEq3qewIJPMC1ElTw9PJvWfxMhZIxn721g+GfYJAd4BZoclLszRc7HkYSFqqXNYZ8Z0H8O7a96lW3g37ux8p9khub1azYNsTbZxWusPtdargdXWvrXQSodauMAUQ9XQwM4aHhtK1Uk/EWMEpeIcmwnWl80ppYYDw9u1a2ePywlRLyIbRjJ54GQe/e1RnlnyDK8NfA2lavvrK+zJQXKx5GEh6tED3R5gc85mJq+YTHRINP0i+pkdkltTtRnKV0oNBaIr35pTSg0tvzVnq8nprT1sSRjzaI7EuI2XjrHiUmqF4yzAbcAK67EZWmu7jlzExMTolStX2vOSQlyyD9d/yFsZb/Fs32cZ2WGk2eG4JKXUKq11zMWPvOh5TMnFkoeFsK1Txaf4y09/Ibcwl9QbUgn3Dzc7JJdT0zxcqwK5hhc2ZSU9RyKJWTijMl1GQloCaw6t4T/X/4e2IRdaCVjURX0VyDW8llvnYsnDwlltPbaVv/z0F2KaxvB+3Pt4KLtOOObyapqHz2uxqIeeNTN7kIUQdeShPHh5wMuMnDWS8b+P58s/fSn9yCaSXCyEe+rQqAMTYifw4v9e5OMNH3P/5febHZJbqqoHuVcdzlO5Z83WPcgOSXrfhLNrHNCYVwa8woNzH2TCggm8ddVbeHq47SCk2SQX14HkYeEKRnUYxYoDK3g7423aNGzD0NZDzQ7J7UiLhQ3IrT3h7L7e/DUvL3uZh3s8TEJ3uzxb5RakxcJ+JA8LZ1dQUsB9v9zH9tztpA5PJbJhpNkhuYSa5mG7NbZYJ7F/xfr1LRUnuBdCOJbbL7uda9tcy5S1U1hzaI3Z4Yh6JLlYCOfg5+XH64Nex9PDk4kLJ1JcWmx2SG7FFj3IF/KqdSqiBOA/SqlorfW3tbyWEMLGlFI8fcXTrD+8nnHzx/Gf6/9Dk4AmZoflViQXCyGaN2jO8/2eZ9z8cbyy/BWe6fuM2SG5DVv1IFd5jNY6Tyk1EpigtZ6nlBpRh2s5LOl9E64kyCeIt4e8zR2z7+Cx+Y/xybBP8PH0MTssdyK5uA4kDwtXE986nvu63sdHGz6iY2hHbr3sUleEFzVhtx5kawK+DWNy+1DrtiFa63n1GoADkN434UrSd6Xz2PzHuKndTbzQ7wVZROQSOEIPsrvkYsnDwpWUlpXyyLxHWLpvKR8O+5BeTevy72cBDtiDrLX+Rmt9a4WE/Ku9ri2EqLu41nGM6T6G77d/z1ebvzI7HHGJJBcL4Xw8PTxJHphMy6CWjJs/jv0n9psdksurU4GslJpe/mCHUiqqjg95JLjaiIUQrurB7g9yVeRVvLbiNZbtX2Z2OMJKcrEQ7iPIJ4i3hrxFUWkRY38by+mS02aH5NLqOoI8VWu9RinVBhiDcbuutoLreG0hhJ15KA9eufIVohpGMf738WTnZ5sdkjBILhbCjUQHR5M8MJnNOZt5dsmz1HebrDirxgWyUupXpdQUpdQtwA7r+1DgA611Ug0+31ApdX+FEY4s6zlchlJquFJqal5entmhCFHvAr0DeXvI25TqUsb+NpZTxafMDsktSS6+MMnDwtUNbDmQRy2PMidrDp9s/MTscFxWbUaQE4GJGKsyjanwGlnD23pPAI2AJ5RS24BXgdjahevYtNaztNajg4NlQEa4plYNW/H6wNfZnrudpxY/JaMX5pBcfAGSh4U7uK/rfVwTdQ3/WvUvFmYvNDscl1TVNG9V0lqvtn75jfUFgFKqJ9AWuNhqAmla67kVPjcCyKh5qEIIR9CvRT/G9RrH6ytfZ9r6aYzuNtrskNyK5GIhhFKK5/s9T1ZeFokLEvnquq+ICo4yOyyXUqMRZKXUB9XdgtNar9Zaf1PVPqXUNuutwCFU6nOzPkmdVeuIhRCmu7vz3VwffT3vrH6H33b/ZnY4bkNysRCiXIB3AG8NeQsvDy8e/e1R8ovyzQ7JpdS0xSINmHvRo873DZAO3ApMVkodVUqNl6VNhXBuSime7fssncM688SiJ9hzfI/ZIbkLycVCiDNaNGjBPwf/kz3H9/DM4mek7a0e1bRA3oHRs1YrWuuJ1tGJMVrrdhh9bnkYvW/blVJTantOIYRj8PPy443Bb6CUYvyC8RSVFpkdkjuQXCyEOEdss1jGWsaSvjudrzd/bXY4LqOmBfLtQIZS6pdLGXXQWmdqradZJ6lvB0yuy3kclTw9LdxNiwYtmNR/En8c/YM3Vr1hdjjuQHLxRUgeFu7o7i53M7DlQF5f+Tp/HP3D7HBcQo1HkK2rLo3h7KjDL7W5UFW9c67W9yZPTwt3NKTVEO7sdCdfbvqSebtlvQkbk1x8EZKHhTvyUB681P8lQv1CGf/7eE4WnzQ7JKdX0wJ5hlJqiNY6q8Kow7BaXquuvXNCCAc3rtc4Ood15pklz3Dw5EGzw3FlkouFEFUK8QsheWAye0/s5eVlL5sdjtOrUYGstc6rh6VI69Q7J4RwfN6e3rx65asUlRbx5OInKdNlZofkkiQXCyEupFfTXozuNpofdvzAnKw5Zofj1M6bB1kpFQWE1PI8uVrrnRc55nZgtFJqBcYIRrrW+mLzdQohnESb4DZM7D2RZ5c8y6cbP+X/uv6f2SE5NcnFQoi6SOiWwNJ9S3lx6Yt0a9yNFg1amB2SU6pqoZBedTiPBnZe5JgdWutQpVQbIA6jdy64DrcHhRAO6uZ2N7No7yLeznib3s160yW8i9khOTPJxUKIWvPy8OLVK19l1KxRJC1M4uNhH+PlUeN14YSVqu8585RSpVprzyq2BwO96uH2oMOLiYnRK1euNDsMIUyRV5jHiB9G4O/lz4zhM/D38jc7JIehlFqltY6x07XcOhdLHhbu7qfMn5i4cCIP93iYhO4JZofjMGqah2v6kN4lq6feOSGEgwv2DealAS+x8/hO3ln9jtnhiEokFwvhHq6Lvo5ro67lg3UfsCVni9nhOB1b9CArG/bOOTSl1HBgeLt27cwORQhT9Wneh9suu40v/viCoa2G0qtpXboF3Jvk4rqRPCzEWU/0eYLlB5bz1OKn+Oq6r/D28DY7JKdxXouFUmpEHc6jtdbfWj9fBowEVF3P4ezk1p4QcKr4FLf8cAseyoPU4akEeAeYHZLpatNiIbn40kgeFsIwb/c8xv42loe6P8SDPR40OxzT1TQPnzeCrLX+5hKv7RLJVQhxaQK8A3ix/4vc+8u9vJXxFkl9kswOyalILhZC1IchrYZwXfR1TF03lcGRg+kU1snskJyC3XqQhRDuJ7ZZLHd0uoOvNn/FigMrzA5HCCHcUlLvJEL8Qnhy8ZMUlxabHY5TsMW8H27Z9yaEqNqjPR9lYfZCnl78NN/e8K20WtiP5GIhBGA8PP1s32d5ZN4jfLDuAx7p+YjZITk8W02MZ6EOfW9cfP5OIYSTKW+1uOfne/jnyn/ydN+nzQ7JnUguFkIAMDhyMDe0vYGP1n/EVZFX0TW8q9khOTSbFMjS9yaEqMjS1MLdne/m0z8+ZVDkIAa2HGh2SG5BcrEQoqLE3oks27+MpIVJMk/9RdR7D7LWWvqahRDnecTyCO0btee5Jc9xouiE2eG4PMnFQojKGvo0PDNP/Xur3zM7HIcmCVQIYRe+nr680O8Fjpw+wpS1U8wORwgh3FKf5n0Y0X4EX2z6gm3HtpkdjsOSAlkIYTddw7syosMIvtz0JRuObDA7HCGEcEtjLWNp4NOA55c+T2lZqdnhOCQpkIUQdvVYr8cI8w/jiUVPUFBSYHY4Qgjhdhr5NSKpdxJrD6/lk42fmB2OQ5ICuR4ppYYrpabm5eWZHYoQDquhT0Ne7P8iWXlZfLzhY7PDES5G8rAQNfOnNn8ivnU8U9ZMYc/xPWaH43CkQK5HWutZWuvRwcHBZocihEPrF9GPa9tcy0frP2JPviRmUX8kDwtRM0opJvaeiJeHF8krks0Ox+FIgSyEMMU/ev0DLw8vJq+YbHYoQgjhlpoENGFM9zH8nv07vzlrYksAACAASURBVO/53exwHIoUyEIIUzQNbMqY7mOYv2c+C7IXmB2OEEK4pTs73Umb4DYkr0imsLTQ7HAchhTIQgjT3NnpTqIaRpG8PJmi0iKzwxFCCLfj7elNUu8k9uTv4d8b/m12OA5DCmQhhGm8Pb1J6pPE7vzdfLrxU7PDEUIIt9Q3oi/xreP5cP2H7Duxz+xwHIIUyEIIU/WL6Ed863imrpvK/hP7zQ5HCCHc0uMxjwPw+srXTY7EMUiBLIQw3fiY8QC8tvI1kyMRQgj31LxBcx7o9gBpu9JYsm+J2eGYTgpkIYTpIhpEcP/l95O2K42l+5aaHY4QQrile7rcQ6ugVryy7BWKS4vNDsdUUiALIRzCPV3vITIokleWS2IWQggz+Hj6kNg7kZ3Hd/LFpi/MDsdUUiALIRyCr6cvE3tPJCsviy83fWl2OEII4ZYGthzI4MjBTFk7hYMnD5odjmmkQL4ApVSyUkorpY4ppdKUUtFmxySEKxvYciCDWg5iytopHDp1yOxwhAOQPCyE/U2InUBpWSn/XPVPs0MxjRTIF7ZDa6201o201vFa60yzAxLC1SXGJlJSVsI/V7pvYhbnkDwshJ1FBkVy7+X3MidrDisOrDA7HFM4ZYGslBqplIqrZp9FKTVaKRVX/m7v+IQQdRfZMJL/6/p/zM6a7baJ2RlIHhbCtd3X9T5aNGjBy8teprjM/Z4LcboC2Zpop1WzLxpI1lpP1Vqna62nAgmXcEsupPwvAettvpC6xi2EqLn7Lr+PiMAIXl72MiVlJWaHIyqRPCyE6/Pz8uPx2MfZnrud6Zunmx2O3TlNgayUilZKpQDRQE41hyUAKZW2pQDJdbzsVK11qtY6HZgOzK3jeYQQteDv5c+E2Alsz93Ofzb/x+xwhJXkYSHcy5DIIfSP6M97a97jyOkjZodjV0prbXYMtaaU2gEkWBNm5e3n9KhZRxuOaa2V9fvRQNsLnD6t8nkrnEsDjbTWuReKLyYmRq9cubJmP4wQokpaa8akj2H94fXMvmU2IX7OP3ColFqltY4xO476IHlYCPewM28nN/9wMze2vZHn+j1ndjiXrKZ52MsewdiDNQGfN6qhtc5VSqGUitZaZ1pv99XkfBZgmta6V+Xz1VvQQohqKaUYHzOekbNGMm39NB6PfdzskMRFSB4WwvVEBUdx+2W389Xmr7ir8120DbnQv21dh9O0WNRAKFwwcda2/y2TCrcJrT13qXULTQhRF+0btefGtjfy9eavyc7PNjsccXGSh4VwQQndEgj0CuRfq/5ldih240oFcr3ef7Um+EzrE9ijgXjggeqOtx63Uim18vDhw/UZihBu7aEeD+GhPHhn9TtmhyIuTvKwEC4oxC+Eey+/l/nZ891mdiFXKpDrXfkT2NZX4oVu61mPidFaxzRu3NieYQrh0poFNuOuzncxO2s2G49uNDscYWeSh4VwDHd2upOmAU15Y+UbOOPza7XlcgWyTAEkhOu5t+u9hPiG8ObKN90iMTs7ycNCuB4/Lz8e7vkwG45u4Jddv5gdjs25UoFc/sR0aMWNFRK1zVdfUkoNV0pNzcvLs/WlhHArQT5BjOk+hmUHlrF432KzwxHVkzwshAsbHj2c9o3a89aqtygude3FQ1ymQC7vVeP8HrhQINcey5NqrWdprUcHBwfb+lJCuJ1bO9xKywYteWPVG5SWlZodjqiC5GEhXJunhyfjeo0j+0Q2M7bOMDscm3KZAtkqHag8t53Ful0I4cS8Pb0ZaxnLtmPbmJU5y+xwRPUkDwvhwvpH9KdPsz6krE0hvyjf7HBsxlkL5FCqflo6ERhVaVuCdbvNya09IWxrWNQwuoR14YO1H1Bc5tq395yA5GEh3JBSisdiHuNY4TG+3PSl2eHYjNMUyEqpEKVUslJqJkZSTlZKpSilRpYfY729l2g9bqRSagKQbI/betbry609IWxIKcWD3R9k74m9/JT5k9nhuB3Jw0IIgC5hXRjccjCf//E5J4tPmh2OTTjNSnrlSbcGx2UAGbaPSAhhhoEtB9IptBPT1k3j+ujr8fJwmjTm9CQPCyHKJXRP4M8//ZmvN3/N/Zffb3Y49c5pRpCdgdzaE8L2lFIkdEtgd/5u5mTNMTsc4WAkDwthH13Du9K/xf+3d2/BbVz3/cC/hyJlUjfCoCTLDmVJS8kVHd8CUo7jjO3YJmJbEamxRMjTRNOZNhXhhzSTPpT4sy//9qUOOG2fMtNCbNNOG3cKgb6IlGWNAStt49RNRML1LbQsEbQUyZEtE4JkWXfr9GEX1BLEjSSAg939fmYwJBaLPb/dBQ5+e/ac3W/iX97/F1y4ekF1OCXHBLmEeGqPqDIevf1RbLh5AwbeHeAVLWga1sNElfPsPc/izOUziHwYUR1KyTFBJiLLqRE18N/jx8TZCUSPRVWHQ0TkSPetvA8P3PoA/um9f8Kla5dUh1NSTJCJyJK8a7zQGjWE3gnhuryuOhwiIkfy3+PH5KVJvHDkBdWhlBQTZCKypBpRg557enA0dRQHjx9UHQ4RkSO1r2pH+y3t+Om7P8XlLy+rDqdkmCCXEAeHEFXWk2ufxJplaxB6JwQppepwqAqwHiaqPP+9fnx68VO8fORl1aGUDBPkEuLgEKLKWlCzALvu3oUPkh/gP0/8p+pwqAqwHiaqvK+v+jruW3Ef/uG9f8DVL+1xEycmyERkaZu1zWhe0oy/f/vv2YpMRKSAEALP3vssTn1xCkPjQ6rDKQkmyERkaXU1dfjju/8Y70++j19+/EvV4RAROdKDtz2Iu5ruwsC7A7h63fqtyEyQicjyulq6cOviW/F3b/8dW5GJiBRItyKfPH8S+xP7VYczb0yQS4iDQ4jUqFtQh+/f9X28c/od/M/v/kd1OKQQ62EidR5ufhit7lYMvDuAa9evqQ5nXpgglxAHhxCp8/SGp7Fy0UqE3gmpDoUUYj1MpI4QAv57/Dh27hgOfHRAdTjzwgSZiGxh4YKF+KO7/gijn4zi0KlDqsMhInKkR29/FBtu3oCBdwbw5fUvVYczZ0yQicg2tm/Yjqb6JoTeZisyEZEK6Zs4Jc4mED0eVR3OnDFBJiLbqK+txx/e9Yf41alf4f3P3lcdDhGRI3lv92Jd4zr89N2fWnbgNBNkIrKVbRu2YVHtIvxs7GeqQyEicqQFNQuws3UnxpJjeOvTt1SHMydMkEuIo6eJ1Fu6cCme3vA0Dnx0AKcvnFYdDlUY62Gi6tDZ0ollC5dZtrGCCXIJcfQ0UXX47sbv4svrX2LPh3tUh0IVxnqYqDo01Dag+45uvH78dXx8/mPV4cwaE2Qisp3bl92Oh5sfxp7De3D5y8uqwyEicqTf3/j7EBD49w/+XXUos8YEmYhs6Xut30PyUhKvTryqOhQiIkdatXgVHr/9cQweGcSFqxdUhzMrTJCJyJYeuPUBrHetx/Njz1t2FDURkdXtvHMnPr/yOfYl9qkOZVaYIBORLQkh8N3W7+KD5AcY/WRUdThERI5034r7cGfTnZZrrGCCTES2tUXbgsabGvH82POqQyEiciQhBHa27kTibAJvfvym6nCKxgSZiGyrobYB2zdsx8HfHsTJ8ydVh0NE5EhPrH0CTfVN+Nexf1UdStGYIJcQr79JVH2sPIqaZo/1MFH1WbhgIZ7Z+AzeOPkGJs5OqA6nKEyQS4jX3ySqPqsWr0LHmg68cOQFy42iptljPUxUnXx3+FBXU4d/G/s31aEUhQkyEdnezlZ9FPXw+LDqUIiIHGl5w3I8te4p7B3fi3NXzqkOpyAmyERke/euuBdfbfoqnv/geVyX11WHQ0TkSDtbd+LitYt46chLqkMpiAkyEdmeEALfa/0eJs5OYOTUiOpwiIgcqbWpFZ6VHkQ+jFT9Jd+YIBORI3jXeLG0bilePPqi6lCIiBxr+x3bcezcsaq/Pj0TZCJyhPraenxH+w6iH0Vx9jKvcEBEpIJ3jRdL6pbgxSPV3VjBBJmIHGPbhm24cv0KXkm8ojoUIiJHaqhtwOZ1m/HasdeqerAeE2QicozWpla0ulvx4pEXq77/GxGRXW27Yxsuf3kZryZeVR1KTkyQichRtm/YjsNnDuM3yd+oDoWIyJHudN+Jje6NeOHIC6pDyYkJMhE5ylPaU6hfUI8XP6zu/m9ERHYlhMC2DdswlhzD2OSY6nCyYoJcQrzFKVH1W7ZwGbxrvNg/sR8Xr11UHQ6VGOthImvYvG4zFtYsrNrBekyQS4i3OCWyhm0btuH81fOIHouqDoVKjPUwkTU03tQI71ovXkm8gkvXLqkOZwYmyETkOG23tGHNsjV44cPq7f9GRGR32zdsx+dXP6/KxgomyETkOEIIPL3+acQ/jeOjsx+pDoeIyJHab2nH6qWr8dLR6rv1NBNkInKkreu3YoFYwDvrEREpkh6sd+jUIRw7d0x1ONPUqg7AKa5evYoTJ07g0qXq62ejQn19PZqbm1FXV6c6FHKo5Q3L8UjzI9h7dC/+5Gt/groafhaJiCqtq6ULP3nrJ3jpyEv4UduPVIczhQlyhZw4cQJLly7F2rVrIYRQHY5SUkpMTk7ixIkTWLdunepwyMG237EdB397EP/12//C42seVx0OEZHjrFy0Eg81P4S943vxg6/9ALU11ZGasotFhVy6dAlNTU2OT44B/ZRKU1MTW9NJuQdvexArG1aymwURkULb1m/DZxc/wy9O/EJ1KFOYIFcQk+MbuC2oGtTW1GLr+q144+Qb+OSLT1SHQ0TkSA81P4QVDSuqqrGCCTIROdrW9VtxXV7H/on9qkMhInKk2ppabNG24I0TbyB5Kak6HABMkMkwODhY9LyJRAJtbW3wer3wer1IpVJljIyovNYsW4N7V9yLofEhSClVh0NE5EhbWrbgmryGAxMHVIcCgAlyQUKIbvNDdTzlEg6Hi543lUqho6MD0WgU0WgULperjJERlV+n1omjqaM4fOaw6lAoC6fUw0ROdsfNd2CjeyOGx4dVhwKAV7HISwjRCyAhpRwUQrgAvA6g+KbWHP5y+H385uNz847P7M7bluH/d3415+vxeByhUAjJZBKdnZ0YHh6G2+2G3+9HOBxGLBaDz+fDwMBAUQnv4OAgEokE3G43QqFQKVeFqOKeXPckfnzoxxgeH8ZG90bV4ZBJuephIqo+W7Qt+OuRv0bibAJao6Y0FksmyEYLQkpKGcvymgdAO4AEAA16xTpjviLKcAHok1LeDABSyhSAtnkFrtiePXswMTGBQCCAvr4+eDweAIDH40E8HkckEgGgd6HI1eWit7cXmqYhGAyiu7sbPp8PsVgMHR0dFVsPolJrvKkRjzQ/glcSr+BP2/60ai4zVM1YDxNRqX1H+w7+dvRvsW98H37o+aHSWCz3KyCE6AAwAMCX5TUNQFBK6TVNiwghElLKxCyLageQSP8IAPAAGJzDcmbI19JbTjt27IDL5UIgEEAgEEAqlUIoFIKmTT9K0zQNvb29OZfjcrnQ3a2f5dy0aRMSiXlvEiLlOls68frx1/Hmx2/ioeaHVIdT1exQDxNR9VnesBzfuO0b2JfYhx987QeoEep6AlsmQTYq3QCAUQC5hjj6AWSe7w8BCCJLRV6ABr0yjkkpU0KIEaPsllkup+pomoZIJIJYLIZQKIRgMDjt9UItyOYW40OHDqGvr6/sMROV28NfeRiNNzViODHMBDkH1sNEVG5dWhcCvwhg9JNRbFq1SVkclkmQjRYDPwAIIQI5ZuvGzIp5BEA0/UQI0YP8lWvUOBWYgH5aMGWUnxJCaEIIzeqtF/39/YhG9U1iTo69Xu9Ui3K+FmRN0+Dz+ZBKpeDxeKa6ahBZWd2COjy59km8fPRlnL9yHksWLlEdUtVhPUxE5fbo7Y9icd1iDI0PMUEuBaOvmoaMVg2jQkW6QpVS7i5ykdkqX8tez8zj8UwNpuvt7Z2RAKcT5mKkW6CJ7KarpQvhw2FEj0Xx9IanVYdjOayHiWi+Gmob4F3jRfRYFH/+9T9HQ22DkjjsdJk3NzA1iCObWQ2HNFonUkaFn67459KHjogs4u7ld2PtsrUYGh9SHYpVsR4monnr1DrxxdUvcPD4QWUx2ClBLsfFeH0A+ozTgX3I039OCNEjhBgRQoycPn26DKEQUbkJIdDZ0omRT0Zw8vxJ1eFYEethIpq39lXtuHXxrUqviWynBLnkjFOBASnlbuNvzlYLY552KWX7ihUrKhkmEZXQFm0LAGDf+D7FkRDAepjIiWpEDbZoW/Dm797Epxc+VRODklLLKH0qTlHZnUKI3WfPnlUVAhHN021LbsP9q+7HcGKYt56eI9bDRDRfXS1duC6v45XEK0rKt1OCnG5VcJsnmirqsvdZk1IOSyl7Ghsby10UEZVRZ0snjp07hrdPv606FKthPUxEJbG2cS3uWXEPhsaHlDRW2CZBNgaFJDCzD5wb+t2eOKiDiIriXeNFQ20DB+vNEuthIiqlrS1bcTR1FGPJsYqXbZsE2RCDfuclM48xncosEAjA5/Ohv79fdShE87K4bjEev/1xHPjoAC5/eVl1OFbDepiISuKJtU+grqZOyWA9qybIbmQfLR3AzBHOfmN62Vm571uuO+cVeo/Xq99NNh6PAwAikQgmJyd5+2myvK6WLnx+5XP8x2//Q3Uo1Yr1MBGVVeNNjfjW6m9h/8R+XL1+taJlWyZBFkK4hBBBIUQEeqUcFEKEhBDd6XmM03sBY75uIUQvgGClTutZue9bOBye1fyJRALJ5I17AcRisalkedOmTYjF2FhE1nb/qvuxctFKdrMwYT1MRJXW1dKF5KUkfnnylxUt1zJ30ktXukXMFwcQL39E8/Dq/wNOvVvaZa66G3jqxzlfjsfjCIVCSCaT6OzsxPDwMNxuN/x+P8LhMGKxGHw+HwYGBuByFR6AHgwGEQwGp+6oNzk5CbdbH5fjcrlw6NCh0qwXkSILahagU+vEP7//z/js4mdY3rBcdUjK2aoeJiJL+OZXvgl3vRtD40P41upvVaxcyyTIViCE6ATQuX79etWhZLVnzx5MTEwgEAigr68PHo8HgH4b6ng8PpXsJhKJnF0uent7MTg4iLa2tmmJdCrFu7+S/XS1dOEf3/tH7E/sxx989Q9Uh0NFqPZ6mIhmp66mDpvXbUb4cBhnL59F402VOTvEBLmEpJTDAIbb29t35Z0xT0tvOe3YsQMulwuBQACBQACpVAqhUAiaNv3ur5qmobe3N+dyotEokskkotEoRkZG0N/fj5aWFiQSCXg8HqRSKTQ1NZV7dYjKTnNpuKvpLgyNDzFBtoii62Eisoyuli78bOxnODBxAM9sfKYiZTJBdiBN0xCJRBCLxRAKhRAMBqe9XqgFORQKAdBbjX0+H3p7exGPxxEOh9Hd3Y1oNAq/31/29SCqhK71XfirX/0VDicP4/fcv6c6HCIix9no3ogNN2/A0PgQE2Qqn/7+fkSjUQCYlhx7vd6pFuV8LcjZeDwehMNh+Hw+aJo21X2DyOqeWvsU+g/1Y2h8CH/m/jPV4RAROY4QAl1aF/5m9G8wcXYC6xrXlb9M3kq19Nrb2+XIyMi0aWNjY2htbVUUUXXiNiGr+NHPf4T//fR/EfPFUFsz93YFIcSolDLzGsFUBtnqYSKyrtMXTqNjsAPfv+v7+KHnh3NeTrH1sGUu82YFvP4mkT11tnRi8tIk/vvj/1YdChXAepjInlYsWoFv3PYNDCeGcV1eL3t5TJBLiNffJLKnh7/yMFw3uZTczYlmh/UwkX1tbdmKU1+cwqFT5b+ULBNkIqIC6hbU4al1T+Hg8YM4d+Wc6nCIiBzp0dWPYkndkorcwIkJMhFREba2bMWV61fw2kevqQ6FiMiR6mvr8cTaJxA9FsWFqxfKWhYTZCKiItzZdCe0Ro23niYiUqizpRMXr13E68dfL2s5TJBLiINDiOxLCIGuli689elbOH7uuOpwKAfWw0T25lnpQfOSZuwd31vWcpggl5CTBof4fD54vV4EAoGpaYFAAD6fD/39/XmnEVnVFm0LBASGExysV62cVA8TOVG6seLXv/s1Tn1xqmzlMEEmAMh557xsdu/ejb6+PkSjUSQSCcTjccTjcQBAJBLB5OTk1PTMaURWdsviW/DArQ9geLwylxkiIqKZtrRsgYTEvsS+spXBBJkAAOFwuOh5Ozo6pu6U5/V6MTIyglgsBq/XCwDYtGkTYrFY1mlEVte1vgsnz59E/JO46lCIiBxp9dLV8Kz0YGh8COW64R1vNa1A8NdBfJD8oKTL3OjeiMD9gZyvx+NxhEIhJJNJdHZ2Ynh4GG63G36/H+FwGLFYDD6fDwMDA3C5XHnL0jRt6v/R0dGpZbjdbgCAy+XCoUP6NQqzTSOyssdWP4ZFtYswND6E9lW8KR4RkQpdLV34izf/Au999h7uXnF3yZfPBNlB9uzZg4mJCQQCAfT19U21Ans8HsTjcUQiEQBAIpHI2eWit7d36v94PI5kMgmPx4NQKDRj3lQqVYa1IFJrUd0ifHvtt/HasdfQ9/U+NNQ2qA6JiMhxvr3223ju189h7/heJsjVTgjRCaBz/fr1eefL19JbTjt27IDL5UIgEEAgEEAqlUIoFJrWIgzoLcTmRDib9HvTSXVLSwsSiQQ8Hg9SqRSamprQ1NQ0YxqRHXS1dOHloy/j58d/js3aZtXhkEmx9TARWdvShUvx2OrHcOCjAwhsCqBuQV1Jl88EuYSklMMAhtvb23epjiUfTdMQiUQQi8UQCoUQDAanvV5MC3IgEJjWatzR0YFwOIzu7m5Eo1H4/X4AyDqNyOrabmnDTx77CR687UHVoVAGq9TDRDR/z977LHru6Sl5cgwwQXak/v5+RKNRAJiWHHu93qkW5XwtyP39/YjFYmhrawMA9PX1obu7G+FwGD6fD5qmTXXfyDaNyOpqRA0eWf2I6jCIiBxNc2mFZ5ojUa7Rf07W3t4uR0ZGpk0bGxtDa2urooiqE7cJOY0QYlRKyZF9FZCtHiYiKrYe5mXeiIiIiIhMmCATEREREZkwQSYiIiIiMmGCXEJCiE4hxO6zZ89mfZ39vW/gtiCicihUDxMRFYMJcglJKYellD2NjY0zXquvr8fk5CQTQ+jJ8eTkJOrr61WHQkQ2k68eJiIqFi/zViHNzc04ceIETp8+rTqUqlBfX4/m5mbVYRARERHNwAS5Qurq6rBu3TrVYRARERFRAexiQURERERkwgSZiIiIiMiECTIRERERkQlvNV0GQojTAI7N4i3LAXxWpnCqEdfX3ri+ua2RUq4oZzCkYz1cFK6zM3CdpyuqHmaCXAWEECPF3BfcLri+9sb1JSty4n7kOjsD13lu2MWCiIiIiMiECTIRERERkQkT5OqwW3UAFcb1tTeuL1mRE/cj19kZuM5zwD7IREREREQmbEEmKhEhhCfjebcQokMI0ZtvmlVlWd+g8bfHNM0260tkZcV8FzO/00TVrtjfmLn8BjFBVsgpyYMTEichRAeAiOm5BwCklDEAKSGEJ9s0JcGWQOb6GnqEEOMAEsY8dlrfHuMRNE2z9QGQnRXaT3bbj8V8F3N8py2riH084zttdUWsc4fxsMU6F/sbY3y2vbNdPhNkReyUPBTBtolTmrEuCdOkZwCkjP8TADpyTLOkLOsLALuklC3Ga4BN1teoXGNSyt0ANOMHxtYHQHZWaD/ZdD8W/C7m+E5bUhH7eMZ3WkGYJVXkOvuM1z1O+VzPBxNkdWyRPBTJlolTAS4ASdPzphzT7ETLaL2wy/pquPEZTRjPbX0AZHOF9pMd96NdvovFKrQPs32nrS7vOkspY1JKv/FUk1LGKxlcmRT8XAshPKbcY1Zq5xoVzZuTKqz0EbpHStkPZ627Yxj7FkIIrx1aZNKMVqY0D4AwgDY47wDILgrtJ+5H68u7D3N8p62uqM+t0YDhz/aaTbnn+ka2IFPZSSn7jSO4JjslTgWkcOOL6QIwmWOaLRh9+bqNp5PQW2Rstb7GKcm4TVpeyFls9V0sFSd+p42GDL8QwqU6lhLI+7meT+sxwARZJUdUWE5InHII48ZpOw1ALMc0uxjBjfVpMZ7bbX07pJQB439HHQDZTKH9ZMf9mPW7aJMkKZti96H5O211BZNFU7/jBIAeWF+hz7VmDFzsAeCebb9rJsjq2C15yMUJiROMg4D29MFAukXCaDFPSSnj2aYpC3iecqzvDuP5uA3Xt8fUhaQDzjsAspNCP6q22495vouvp+fJ/E5bXMEDgizfaasrtM4dmJ5AW35AZqHPtZRyUEo5aEyb9cEgbxSikHFUk4DeYd62d7ox1jMJfT37TdNsv+5kfabLXyWh/8D4pJSxbJ9hfq6tIce+G5VStuV6nawl3z7O9Z1WF21pFFhnF4AdxqxtpgF7lAMTZCIiIiIiE3axICIiIiIyYYJMRERERGTCBJmIiIiIyIQJMhERERGRCRNkIiIiIiITJshERERERCZMkC1MCNEhhJDGY8ZdcYQQWq7XyhxXUAhR9dcPFEKEhBBnVGyjShFC9AohoqrjyMcKMRJRcYw7l8liH6b3sR6gqlKrOgAqmQAAXtC+SEKICAAPgF3GJMvfVSiHTdDvoFTNrBAjERUnBsCbMc0PoNv4m6uuZT1AVYUJsj3sBtAjhOg23VaR8uuGfvekim8vIUQQQEhKWfakXErpU1l+MeVmi5GIrElKmULG7bmFEB7odW4sV73DeoCqDbtY2EME+lF5n+pArEAIkb5ffUpB2S4AvdBbrytOVfmq15uIiGg2mCDbRxCAxzhSJyIiIqI5YoJsE1LK3dBbRIP55ss2gE4I0WMMmHBlzHdGCOERQowar48aAwM1IUTUmDYuhOjOUZbHmO+MMd+M/mXGgI5py8+Mwfg/ZCzDlbmMPOWml9ljei0IYNx4mp4na/zFLtMUa97tavR7PmO8v4eo6gAABRdJREFUHDFeO5PxnpJuD3Nc5Srf+DxETAMeo6ZW+oLl5th2xWzvM1nm82TM123EKbPFRkTVIbMemO9vUL66bA6x5a1HCpWV8Vt4xvifjVlVjgmyvTwHoKOECYALeveN56APrtCM5+lHus9YJMf7BwCEoA8gBPSEdFqlYrw3ZiwrkTmPMd8ogHbo/VfzdoswljlqLMsLIAwgJIQIGbM8hxsDSAIAWgr1Qy5imcXaZSrbD6AFQFtGOSXdHhUq35+xfDf07VVUuZlmsb3Tn88Qpn8+08tJPx80ykuvU8GDLCKqCnP6DSq2LitGoXqkUFlGspyuD3dB/91JAXhmtrFQhUkp+bDoA/qIXwmgw3juMp6HjOea8bzH9J6gvtunLafHmM+VOV962RnTek3Tuo1pngJlaObYjGlnAAQz5htNTzOVF53FNjljLiPHdkrH0l3CZRa7XXOWXabtMS2uSpQPvZ/xtDIKlJsZY9HbO6OM9DRXxmfTlS9ePvjgo7wP6OMPJAAtzzyZ9cB8foPy1mWzjD1vPVJEvTk+mzqbj+p5sAXZRqTempe+okWpWslGTP+nuyaYRyinRyS7C8SWMN7XDkyNanYB6BXTr4npwcyBXH4UwThSd0FvUTSXHYN+xD7rUdLlWGaOckq+PVSVL6WMG//OpbVmttvb/Fkcz/Ha60K/xipPaRJZz6x+g2ZZlxUjZz1SqCyjFVlDga6PVJ14mTf7CUJvueyBfkpoXmT2U/hzvfpDAkaCjBvJU1uW5SUzYij2cmTpZWabP4E5JGxlWma+ckq5PSpWvnGa8RljOfPZJrPa3jk+n1OvCSG80L8TQSPOmJQy8xqtRFSl5vAbVHRdVmz5eeqRQmW1m/4ni2GCbDNSyoQQYhD6Jd9iheavMA03Ep+pBKiECV96ORqAeMZrGqa3RKhcZr5yypUAl618od/9Kt1Kkt4eo7nfUVQcJdneRstzm3FGpQdAUAjRK6Xsn2N8RFTdSl6X5qpHYPqNzdFwUNQZVqpO7GJhT89BP+1T7CCAnAOmSsWoWDpgVCjGafisV92YR/eQdPI0rQuA0bqZHuxRyWVm265ZWxLKtD2yKWn5xinEDgABKeVuYznZfpSKbUEpxz6ElDJlJMVx6HfsIiIbKmddmlmPFCrLSJoTuDFQvWSxUPmxBdmGpJRxIUQc+sCITOOAft976MlqB/Qj4pIT+uW9ArjRupiCnryn7YJ+2a/0lQhcuHEr0ln3szVOhfmhX/HADf3qB+myB41WgHIts6jtaiwPAPzpv6ZT/iXdHnnWp2TlG2csUtBbVAB9H8/4MShQbuZ8JdmHQr80XMBYl3T3DA8y+jcTke3krcuMA/sggOdMYyayKqIeKVRv+qFf1SICvT5z48ZYCnb3qmJsQbavGUmKYQ/0o98g9Na4JpRh0JchDCBqPJIA2sz9yaR+ebU26BVKFPpl4bIebRdL6teDTvcNi0Bft4Ccx21Mi1zmbLZrP/QEOghTN4JybI8Klb8LeqUfwY31L7rcTCXch3ug98P3m5bTbyyfiGyqiLrMBf3qFO1ZFzBd3nqkUFnp7hm4cQYsfZm3cv3uUokIKWXhuYiIiIhsQggxDsBXqAWZnIstyEREROQ4TI4pHybIRERE5BjGALlSd10jm2EXCyIiIiIiE7YgExERERGZMEEmIiIiIjJhgkxEREREZMIEmYiIiIjIhAkyEREREZEJE2QiIiIiIhMmyEREREREJv8HLLg4stOxeZgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.sparse.linalg as la\n",
    "from scipy.sparse import csc_matrix, csr_matrix\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "\n",
    "%matplotlib inline\n",
    "n = 150\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "\n",
    "plt.figure(figsize=(10, 5))\n",
    "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))\n",
    "for restart in [5, 40, 200]:\n",
    "    hist = []\n",
    "    def callback(rk):\n",
    "        hist.append(np.linalg.norm(rk) / np.linalg.norm(rhs))\n",
    "    st = time.time()\n",
    "    sol = la.gmres(A, rhs, x0=np.zeros(n*n), maxiter=200, restart=restart, callback=callback, tol=1e-16)\n",
    "    current_time = time.time() - st\n",
    "    ax1.semilogy(np.array(hist), label='rst={}'.format(restart))\n",
    "    ax2.semilogy([current_time * i / len(hist) for i in range(len(hist))], np.array(hist), label='rst={}'.format(restart))\n",
    "    \n",
    "\n",
    "ax1.legend(loc='best')\n",
    "ax2.legend(loc='best')\n",
    "ax1.set_xlabel(\"Number of outer iterations\", fontsize=20)\n",
    "ax2.set_xlabel(\"Time, sec\", fontsize=20)\n",
    "ax1.set_ylabel(r\"$\\frac{||r_k||_2}{||rhs||_2}$\", fontsize=20)\n",
    "ax2.set_ylabel(r\"$\\frac{||r_k||_2}{||rhs||_2}$\", fontsize=20)\n",
    "plt.sca(ax1)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.sca(ax2)\n",
    "plt.yticks(fontsize=20)\n",
    "f.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXt8VdWV+L8rJCQkkEQExEIhIgpVRB4pjyIQcMBoaQXrjFZLiVbRVsVpUcHqjODo9GXrgxE7VBRR6mNaRlvFws8HqDysKLUtVRFGUF7K0/DGkPX7Y58bbi735j7PvZeb9f18zufk7uc6+5ycdfbea68tqophGIZhZBt5mRbAMAzDMMJhCsowDMPISkxBGYZhGFmJKSjDMAwjKzEFZRiGYWQlpqAMwzCMrMQUlGEYhpGVmIIyDMMwshJTUIZhGEZWkp9pAYzkadeunVZUVGRaDMMwjJh4++23t6tq+2jpTEHlABUVFaxcuTLTYhiGYcSEiGyIJZ0N8RmGYRhZiSkowzAMIysxBWUYhmFkJaagDMMwjKzEFFSWICLDROQPIrJJRFREajItk2EYRiYxBZU9tAb+DtwIHMiwLIZhGBnHzMyzBFVdACwAEJE5ftd38EAthYWtkTz7RjEMIztptgpKRC4GhgN9gLOBNsA8Vf1OE3k6A3cC1cCJwBbgWWC6qu7yXegUMvWJi/lb3iZOOXwC7QvOpl2Hb9LlS2fQ9cQSurQtpkObQvLyJNNiGobRjGm2Cgq4HaeY9gIbgZ5NJRaRU4FlQAfgOeB9YABuSK5aRIao6g5fJU4hPdtVUrvtc/5RuIs9LV6DXa/RdSu03d+O2r29WHdoKB3KO9ClbbE7PMXVpW0xX27biuKWzfnRMQwjHTTnt8wPcYppLa4n9WqU9DNxymmSqs4IBIrIr7yy7gauDQq/C7gtSpkjVHVx3JKngGvH/SfXAkfqDvOPD55j2do/snznP/hr+Ta+OGExBfoq7Q8XULT3y/xlS18er+1Hvfe4tGyRxx9vOIceHdtkQnTDMJoJoqqZliHjiEgVTkGFHeITkW7AOmA9cKqq1gfFtcEN9QnQQVX3eeHtgHZRqv5YVfeHqW8vcL2qzolF/srKSk2Vq6P9+7ezavVTLN/wEiv2fMQHee5S29Qr/aSUk/J68tLaM7jr0n9h6BldUlKnYRjNCxF5W1Uro6Vrzj2oeBjpnRcFKycAVd0jIkuB0cAg4GUvfDuwPa1SpoDi4nYM+er1DPnq9QDs2L6GN1fPY8XmpSw/sJUl8hac+hYnPPszmJ8HxSdCq7ZQfIJ3bhv5XNwW8luBCCDhz2LzXoZhOExBxUYP77wmQvyHOAV1Op6CihcRaQ10937mAV1EpA+wU1U/DpN+IjARoEsX/3oyJ7Y7nQuGT+cCQOvreXr5PO5e+3NWdr+QXmUnwYFdsH8nHNgJn6525wO7oLEej5MIyosElVdSSs8UpmE0yVUvQcdevhRtCio2yrzz5xHiA+HlSdRRSeN5sOne8RhQE5pYVWcBs8AN8SVRb8xIXh4dOnwN1sLHJw2A4RPDJ6yvh0Ofe4rLU2D7dzjlVXcQVAEFxTvr0bPWHxvWcE6EJJrGhr8NIzrFJ/pWtCmo1BD4zE74jeYZS2T953rrglIA9tXVRk6UlwetTnCHYRhGgtgqzdgI9JDKIsSXhqTLWfLz8tEjhU0rKMMwjBRgCio2PvDOp0eIP807R5qjyin0SIkpKMMwfMcUVGwE5oZGi0ijNvPMzIfg/OetSLdgmUCPFLPvC1NQhmH4iymoGFDVdcAioAK4LiR6OlACzA2sgcplBE9BWQ/KMAyfabZGEiIyFhjr/ezonQcHOWrdrqo3BWX5Ac7V0QMici7wHjAQGIEb2ovmNSJncArq00yLYRhGjtNsFRTOSeyEkLBu3gGwAWhQUKq6TkQqOeos9gKcB4kHcM5id/oucZagR4rZbz0owzB8ptkqKFWdBkyLM88nwBV+yHO8ICLokWIOHNlHXX0d+XnN9hEyDMNnbA7KiBs9UgzA54dy3qreMIwMYgrKiBtTUIZhpANTUEbcBBTU7kO7MyyJYRi5jCkoIy4CZuZgCsowDH8xBWXEjQ3xGYaRDkxBGXFjPSjDMNKBKSgjLkSA+kLypIX1oAzD8BVTUFmCiFwnIn8VkVrvWC4iX8+0XOERSvJLrQdlGIavmILKHjYCU4B+uM0LXwGeFZHeGZUqAsUtSq0HZRiGr5gbgCxBVZ8LCbpNRL4PDAb+mgGRwiLenorWgzIMw2+abQ9KRC4WkRki8ro3pKYi8kSUPJ1F5BER2Swih0RkvYjcJyIp3TpWRFqIyKVAa5yD2qzDFJRhGH7TnHtQtwNnA3txw2s9m0osIqfilEUH4DngfWAAcCNQLSJDVHVHMgKJyFnAcqDIk2ucqv4tmTL9oji/lK37msX+jIZhZIhm24MCfojbIbcU+H4M6WfilNMkVR2rqlNVdSRwL9ADuDs4sYjc5fXKmjqqQur4AOdlfRDwEPCYiPRK6ip9otjrQalqpkUxDCNHabY9KFUN7JKLiDSZVkS6AaOB9cCDIdF3ABOB8SIyOWjTwvuAJocMgY9DZDoMrPV+rhSRr+IU6feilJM2Ak1Vkl/KF/VfcKDuAMUFxZkVyjCMnKTZKqg4GemdF6lqfXCEqu4RkaU4BTYIeNkL3w5sT7LePKAwyTJ8oTi/FHCLdU1BGYbhB815iC8eenjnSJMuH3rn0xOtQER+KiJDRaRCRM4SkZ8AVcC8COknishKEVm5bdu2RKtNmOIWbQDzJmEYhn+YgoqNMu8caeFPILw8iTo64oYEP8D1wr4KnK+qL4ZLrKqzVLVSVSvbt2+fRLWJURLUgzIMw/ADG+JLDYFJrIQtBlS1JjWipIfAEJ8t1jUMwy+sBxUbgbdwWYT40pB0OU+x9aAMw/AZU1Cx8YF3jjTHdJp3bjYLg2wOyjAMvzEFFRsBk/TRItKozUSkDTAEOACsSLdg6SZgZt4iL5/WBa1tiM8wDN8wBRUDqroOWARUANeFRE8HSoC5QWugmgVlhWXWgzIMwzearZGEiIwFxno/O3rnwSIyx/t7u6reFJTlBzhXRw+IyLnAe8BAYARuaO8234XOIlShvLDcFJRhGL7RbBUUzqXQhJCwbt4BsAFoUFCquk5EKoE7gWrgAmAL8AAwXVV3+i5xFiAc9bpRXljO5wdtiM8wDH9otgpKVacB0+LM8wlwhR/yHI+UFZaxoXZDpsUwDCNHsTkoI0GUssIyM5IwDMM3TEEZCVNeWM6eL/ZQV1+XaVEMw8hBTEEZcRHs+L2s0K1btl6UYRh+YArKSIiAFR/A54dNQRmGkXpMQRkJ06CgrAdlGIYPmIIy4iJ4iC+goHYftLVQhmGkHlNQRsIE5qBssa5hGH5gCspICMWG+AzD8BdTUFmCiEwTEQ05tmZarlCCPUmUFJSQL/nWgzIMwxearSeJLOUD3DbvAY5kSI6YEBFzGGsYhm+Ygsou6lQ163pN4VBv7+DywnIb4jMMwxea7RCfiFwsIjNE5HURqfWG1J6IkqeziDwiIptF5JCIrBeR+0TkhBSJ1U1ENonIRyLylIh0i54ls1gPyjAMv2jOPajbgbOBvcBGoGdTiUXkVNx2Gx2A54D3gQHAjUC1iAxR1R1JyPMmUOOV28GTb5mInJlkuSkl2MwcXA/q4z0fZ0YYwzBymmbbgwJ+iNvCvRT4fgzpZ+IUxyRVHauqU1V1JHAv0AO4OzixiNwVxugh9KgKpFfVF1X1GVX9q6q+BIzB3Z/QLUGyivIiG+IzDMMfmm0PSlUD27gjod2CELyhttHAeuDBkOg7gInAeBGZHLSr7n1Ak0OGQMSuh6ruFZHVwGlRysgIipuECgzxqWrUdjQMw4iHZqug4mSkd16kqvXBEaq6R0SW4hTYIOBlL3w7sD3RCkWkCDfs+Gq0tOkkVAWVF5bzRf0XHKg7QHFBcUZkMgwjN2nOQ3zx0MM7r4kQ/6F3Pj3RCkTkHhEZLiKniMhA4HdACfBYhPQTRWSliKzctm1botUmTYO7IzOUMAwjxZiCio0y7xxpsiUQXp5EHZ2BJ3FroeYDh4BBqhp2y1pVnaWqlapa2b59+ySqTYyAmbm5OzIMwy9siC81BEa+NNECVPXSFMmSVqwHZRiGX1gPKjYCPaSyCPGlIelylnBm5mD++AzDSD2moGLjA+8caY4pYGkXaY4qZ7EhPsMw/MIUVGwELOlGi0ijNhORNsAQ4ACwIt2CZYrAWKYpKMMw/MIUVAyo6jpgEVABXBcSPR1nbTc3aA1UDtN4jK8gr4CSghIb4jMMI+U0WyMJERkLjPV+dvTOg0Vkjvf3dlW9KSjLD3Cujh4QkXOB94CBwAjc0N5tvgudpZQXllsPyjCMlNNsFRTQh2PdCHXzDoANQIOCUtV1IlIJ3AlUAxcAW4AHgOmqutN3ibMI1aMGi+Yw1jAMP2i2CkpVpwHT4szzCXCFH/Icz5QXllN7qDbTYhiGkWPYHJQRF+Hc7VkPyjAMPzAFZSSNzUEZhuEHpqCMpCkvLGfP4T3U1ddlWhTDMHIIU1BGXITbUCOwFqr2sM1DGYaROkxBGUlj/vgMw/ADU1BGQgRZmZs/PsMwfMEUlBEX4XbNbehBHbQelGEYqcMUVJYgIutFRMMcL2RatmiYPz7DMPyg2S7UzUK+CrQI+n0y8DbwTGbEaRoN2vrKhvgMw/ADU1BZgqo22rddRL4H1AL/kxmJYqekoIR8ybcelGEYKSUtQ3wiciQd9cSDiFwsIjNE5HURqfWG056IkqeziDwiIptF5JA3LHefiJyQYtkE+B7whKruT2XZyRLOzFxEzJuEYRgpJ109qHDvtUxzO3A2sBfYCPRsKrGInIrzZt4BeA54HxgA3AhUi8gQVd2RItlGAacAD6eoPN8pLyy3IT7DMFJKuowkNHqStPND3A65pcD3Y0g/E6ecJqnqWFWdqqojgXuBHsDdwYlF5K4IRg/BR1WEuq4G3lLVvyR6cX6jIXfUelCGYaSarLDiE5EiEekcJvxMv+pU1VdV9UPV0FftsYhIN2A0sB54MCT6DmAfMF5ESoLC7wO+EuX4c5i6OgAXAr+J85LSQjhnsWD++AzDSD0ZN5IQkXHA/cAuEckHrlTVN73ox4F+GRPuKCO98yJVrQ+OUNU9IrIUp8AGAS974duB7QnUVQMcAp5KWNoMUF5Uzt+2/y3TYhiGkUNkQw/q34H+qno2bgPBR0TkMi8uW+auenjnNRHiP/TOpydTiWcccRXwlKruiZJ2ooisFJGV27ZtayqpL0Qa4ouhQ2oYhhETGe9BAS0DJtaqulJEhgHzRaQ72TN3VeadI1kBBMLLk6ynCjgN+E60hKo6C5gFUFlZmfF2Ki8s54v6LzhQd4DiguJMi2MYRg6QDT2oz0Skd+CHZwk3CjdH0ztiruwi0NNLSlF482KiqsfMTWULEqFTaw5jDcNINQkpqBQbNYwHPgsOUNXDqvptYHgi8vlAoIdUFiG+NCRds6Ospbk7MgwjtcStoDyjhjXACyKyWkQGBkU/Hm95qrpRVbdGiFsab3k+8YF3jjTHdJp3jjRHlXOEdhXNH59hGKkmkTmogFHDNhGpBB4TkbtV9bdEMWoQkY9IbBjsPlV9IIF8qeJV7zxaRPKCLflEpA0wBDgArMiEcOmkKTNzgNpDtmmhYRipIREFlYxRQ00C9YFbf5QxVHWdiCzCmZJfB8wIip4OlAD/rar7MiFfNlBeZHNQhmGklkQU1Gci0ltV/wrOqEFERgGPEcWoQVWXJFCfL4jIWGCs97Ojdx4sInO8v7er6k1BWX6Ac3X0gIicC7wHDARG4Ib2bvNd6Cwi1Jzc5qAMw0g1iSio8UBdcICqHga+LSL/lRKp0kMf3LqrYLp5B8AGoEFBeb2oSuBOoBq4ANgCPABMV9WdvkucxRS0KKCkoMT88RmGkTLiVlCqurGJuGwxaoiKqk4DpsWZ5xPgCj/kyQXM3ZFhGKkkrQt1j2MjCSMGzGGsYRippEkF5YNCqUmgLMiwkYRxLOEeCttywzCMVBKtB1WTYLnrwwVmk5GEkRiRzMzB9aA+2fNJ+oQxDCOnaVJB+alQROR/gP/n+ZRDRHrgrAAXh25/bhwf2ByUYRipJJO++IYBfwEQkROBN3E7yK4WkbMyKJcRC2HG+MoLy9lzeA919XXHRhqGYcRJuueggmmDM9MG+BZuWLASt/D1buCbCdRr+Iw0McYXcHdUe7iWtkVt0yWSYRg5SlrnoEL4GDgV+AS4GJirqnXeQtnjxlzdOEqwR3NTUIZhJEvG5qCAR4AHRWQBzhvDtUEy2YZCWY6G6VgHFJRZ8uUGqsqePXuora1l//79HDlyJNMiGVlGfn4+ZWVltG3blvz81K9aSrhEEblOVR9MNL+q/twbLjoPuElV/8+LGoDz4tCs8JzO/gcwDugArAJuVNW3MipYHDT0oA6aocTxjqry2WefsW/fPtq2bUvHjh1p0aJFk0O8RvNCVTl8+DA7duzgk08+oWvXruTlpdasIZnSRorIeYEfItJaRObHU4Cq/lxVz1XV+4OCTwKeSkKu45WHccp6AnAWsAh4SUQ6ZVSqEJp6PdmWG7nDnj172LdvH127dqW8vJz8/HxTTkYjRITCwkJOPvlk8vPz2bVrV8rrSKZP9l1goYhsBOqB3wK/SlYgVf15smUcb4hIK5yhyLdUdbEXPE1EvgF8H7g9U7LFgw3x5Q61tbW0bduWFi1aZFoUI8sREcrLy9m1axcnnnhiSstOZMPCB0TkSqAH8D3cJoXzgO+qaswbForIcBF5W0Q+FpFXRORXIvJdETlLRHx3wSQiF4vIDBF5XURqRURF5IkoeTqLyCMisllEDonIehG5T0ROSFKcfKAFcDAk/ABwTpJl+4KGse0sKSghX/KtB5UD7N+/n9atW2daDOM4obi4mAMHDqS83EQUwULcgtpbgDOAzriN+r4uIl1V9fkYy5kN/B24F+gO9MVZ83UGDgNFCcgWD7cDZwN7gY1Az6YSi8ipuO02OgDPAe/j5stuBKpFZIiq7khEEFXdIyLLgdtF5O/AVuDbwGBgbSJl+kVTozwiYv74coQjR45Y78mImby8POrr66MnjJOYFJSItMBts/Gcqr4AvBAUVwj0wr3szwViVVAdgdFBxhGB8trilJXf/BCnmNYCwzm6a24kZuKU0yRVbdiwUER+5ZV1N0ctERGRu4i+R9SIoCG98TjLxo3AEeAd4EmgX2yXkx2YP77cweacjFjx61mJSUGp6hERmQksAXaFxB0C3vaOeFgCnAI0UlDevkovx1lW3Khqg0KK1rgi0g23m+56INRy8Q5gIjBeRCYH7ap7H9DkkCFuLVhAnnXAcBEpAUpVdYuIPA18FP1q0k+k1dvWgzIMI1XEM8S3HDfvlPAL09s2/S/Au7gdeP9NRP56HPjeG+mdF6lqo36sNzy3FKfABuEpV1XdDmyPtyJPwe3z5rXOww2lHjeYw1jDMFJFPEYSs4C7ReSUJOpbhZu/ugdnSj4M+EBEHhORa0VkoGfRlm308M5rIsR/6J1PT7QCETlPRM4XkVNEZBRuyPED4NEI6SeKyEoRWbltW/r0uzRpaG5DfIZhpI54elBPeue/e94fluDmSf6iqvtjKUBVpwT+FpGTcHNNfbzjX3HGEioiH6rqGXHI5jdl3jnSmzcQXp5kHT/BGYnsBH4P3KaqX4RL7HmBnwVQWVmZiL9EXwgoKFW1OQzDyGJee+01XnvtNd566y3eeusttmxxrlG3bdtGu3btMiydIx4F9WWcIjnbO9+A86UXt0LxjC7Oxxld/CkovJVXdu845MoGAm/ihBWFqj4DPJMacfwnnJk5uCG+w/WHOVB3gOIC81hlGLGyePFiRowYwfDhw1m8eLHv9U2aNIl3333X93qSIWYFpaqbgE00tuArximsuBRKJKMLVT2Am+taHk95aSDQQyqLEF8aki5nidYpCl6sawrKMLKXUaNG8a1vfYvKykr69+/PSSedlGmRjiGpBbHe0F6iCiVpo4s08oF3jjTHdJp3jjRH1WwI9mh+cuuTMyyNYRiR+MUvfpFpEaKSyQ0LU2F0kS4CJumjRaRRm3lOXofgvD6sSLdgmSKcN3Mwf3xG80FEGuZZZ8+ezcCBAyktLUVE2L3bPf+bN2/m+uuvp3v37hQVFVFcXEyXLl2orq5m1qxZDWVVVVUxYsQIAJYsWdJQtohQVVWV9mvLFjK5YWHSRhfpQlXXeSbyo4HrgBlB0dOBEuC/g9ZANVvMH5/R3LjhhhuYOXMmQ4YMYcyYMaxZswYRYcuWLfTv35+tW7fStWtXqqurKSwsZNOmTaxYsYL169czceJEAKqrqykqKmLhwoWcdNJJVFdXN5Tfs2eTTm5ymkxuWJgyo4tEEJGxwFjvZ0fvPNjbMBFgu6reFJTlBzhXRw+IyLnAe8BA3F5Wa4juNSIniGaXV150dIjPyF2m/3E1/9hcm2kx4uKML5VyxzfOTHm5jz/+OMuXL2fAgAGNwu+//362bt3KNddcw0MPPdTIqvXQoUO8+eabDb+nTp3KoEGDWLhwIT179mTOnDlh65ozZw5XXHFF3DI++uij1NTUxJ0v02Rsw8JUGl0kSB/c1hbBdPMOcHtSNSgorxdVCdwJVAMX4LasfwCY7nnAaPaUtbQhPqN5ccsttxyjnAA+/fRTwPWOQpdcFBYWMmzYsLjr6t69OxMmhL62Yst3POK71/B4SNLoIt66pgHT4szzCRD/50sOEsnMvKBFASUFJTbEl+P40RM5XrnooovChg8YMICZM2cyZYpb/jlq1ChKSkqSquucc87hnHOycoMDX8jkHJRxPBLD2tvywnLrQRnNhq5du4YNHz9+PIsWLeK3v/0t48aNo0WLFvTq1Ythw4Zx6aWX8rWvfS3Nkh5/pHUOyhRe88AcxhrNiVatwntny8vLY968edx66608//zzLF26lKVLlzJjxgxmzJjBlVdeyezZs+Oq64033uDhhx+OW8arrrrquOx5pXsOqibBfOtTKIORApr6yjB/fIZxlF69etGrVy8A6uvrWbBgAZdddhmPPPIIl1xyCaNHj465rLVr1/LYY4/FLUNVVVXuKahU46fRhZEeojmLBfNobhiRyMvLY8yYMVx44YU88cQTvPvuuw0KqmXLlgDU1dVFzF9TU3NcWuMlSiYX6ho5is1BGQbMnTuXd95555jwHTt2sHy5swMLnr/q1KkT4HpJTSmp5kRae1A2B9U8KC8sZ8/hPdTV15Gfl1WGooaRNubPn8+ECRPo1KkTffr0oby8nB07dvD666+zb98+hg4dyrhx4xrSd+3alb59+7Jq1Sp69+5N//79KSwspEePHtx8880pl+/hhx8OO59VXV1Nfr77v+3Xrx8zZ85Med2xku63R02C+danUAYjFUSyM+eou6Paw7W0LWqbLokMI6uYPHkyFRUVLFu2jJUrV7Jr1y7atWtHv379qKmp4fLLL6egoKBRnvnz5zNlyhSWLFnCk08+yZEjRxg+fLgvCmrjxo2NFgsHePvto5ujFxUVpbzeeEiXghJo3nNQIjIMt/C3P/Al4ApVnRNvmkwTyxZPwQ5jTUEZuYo28ZEGMHToUIYOHRpXmRUVFTz99NPJiBUz06ZNY9q0aWmpK1HSMgelqjbXBa2BvwM34hzLJpom6zF/fIZhpAKbIEgTqroAWAAQ5O8v7jTZQjQzc4DdB81QwjCMxMnZno2IXCwiM0TkdRGpFREVkSei5OksIo+IyGYROSQi60XkPhE5IV1yZzuxbOJuW24YhpEKcrkHdTvO8exeYCPQpM96ETkV5628A/Ac8D4wADfcVi0iQ1R1h68S5wgBBWVDfIZhJEPO9qCAH+J2wC0Fvh9D+pk45TRJVceq6lRVHQnci9v59+7gxCJyl9cra+qoSu0lZQ9NzQ+3LmhNvuRbD8owjKTI2R6UqgZ2wT3G1X0oItINtxnheuDBkOg7gInAeBGZHLQp4X1Ak0OGwMdxiJwziAilhaV8fth6UIZhJE7OKqg4GemdF6lqfXCEqu4RkaU4BTYIeNkL3w5sT6uUWUA0ZR/A/PEZhpEsuTzEFw89vPOaCPEfeufTE61ARFqLSB8R6YNr9y7e7y7xpAlKO1FEVorIym3btiUqlm+YuyPDMJLFFJSjzDtH+uQPhJcnUUclsMo7WgHTvb/vjDMNAKo6S1UrVbWyffv2SYiVGNEWKdqWG4ZhJIsN8cVGYFwrET+CLqPqYqJYaceSJtPEKlx5YTmrt6/2VRbDMHIb60E5Aj2ksgjxpSHpjCgEhvii9bQMwzAiYQrK8YF3jjTHdJp3jjRH1eyIpnbKCss4XH+YA3XHrccmwzAyjCkoR8AkfbSINGoTEWkDDMH5xluRbsGOV8wfn2EYyWIKClDVdcAioAK4LiR6OlACzA1aA9VsidHKvJFHc8MwjETIWSMJERkLjPV+dvTOg4OcsG5X1ZuCsvwA5+roARE5F3gPGAiMwA3t3ea70DmE+eMzjOzl888/Z+HChbz44ousWrWKDRs2cPDgQTp37sw//dM/MXnyZLp3755pMXO6B9UHmOAd53lh3YLCLg5O7PWiKoE5OMU0GTgVeAAYbH74GhPN9sGG+AwjPhYvXoyIUFVV5Xtdv/jFL7jkkkuYM2cOBw4cYOTIkZx33nkcPHiQX//61/Tu3ZsXXnjBdzmikbM9KFWdBkyLM88nwBV+yJMrSIyG5uVFNsRnGNlKcXExP/rRj7j22ms57bTTGsK/+OILpkyZwr333st3vvMd1q1bR9u2mdt0NJd7UEYGKWtpQ3yGka38+Mc/5pe//GUj5QRQUFDAPffcw+mnn87u3bsz3osyBWUkRDQz84IWBZQUlNgQn5GziEiDb8rZs2czcOBASktLERF273YfZps3b+b666+ne/fuFBUVUVxcTJcuXaiurmbWrFkNZVVVVTFixAgAlixZ0lB2uob8gsnLy6N3794AbNy4Ma11h5KzQ3yGT8Th58L88RnNgRtuuIGZM2cyZMgQxowZw5o1axARtmzZQv/+/dlvfSNFAAAgAElEQVS6dStdu3alurqawsJCNm3axIoVK1i/fj0TJ04EoLq6mqKiIhYuXMhJJ51EdXV1Q/k9eza5lZ0vfPihcz968sknp73uYExBGb5h/vhymBenwta/ZVqK+Oh4Fpz/05QX+/jjj7N8+XIGDBjQKPz+++9n69atXHPNNTz00EONdgI4dOgQb775ZsPvqVOnMmjQIBYuXEjPnj2ZM2dO2LrmzJnDFVfEP03+6KOPUlNTE1PaP/3pT7z77ru0atWK888/P+66UokpKMM3bMsNozlwyy23HKOcAD799FPA9Y5Ct6kpLCxk2LBhcdfVvXt3JkyYkFC+WNi8eTPf+973ALj11ls56aST4q4rlZiCMhIiFh97ZYVlfLLnkzRIY6QdH3oixysXXXRR2PABAwYwc+ZMpkyZAsCoUaMoKSlJqq5zzjmHc845J6kyIlFbW8uYMWPYvHkz5513Hrfdlvmln2YkYcRFrJ4kwOagjOZB165dw4aPHz+eyy67jDVr1jBu3DjKysro06cPkyZNYtmyZWmWsmn27t3L+eefz6pVqxg6dCjz588nLy/z6sF6UIZvlBeWs+fwHurq68jPs0fNyE1atWoVNjwvL4958+Zx66238vzzz7N06VKWLl3KjBkzmDFjBldeeSWzZ8+Oq6433niDhx9+OG4Zr7rqqog9r3379vH1r3+dZcuWMXDgQF544QWKi4vjrsMP7K2RBkRkGHAT0B/4EnCFqs4JSXMdcA3OHyDAauAuVc38cu4ECbg7qj1cS9uizC32M4xM0qtXL3r16gVAfX09CxYs4LLLLuORRx7hkksuYfTo0TGXtXbtWh577LG4ZaiqqgqroPbv38+YMWN47bXXqKys5E9/+hNt2rSJu3y/yHwfrnnQGvg7cCPOK3o4NgJTgH44l0uvAM+KSO+0SBgj8eymGFBQZihhGI68vDzGjBnDhRdeCMC7777bENeyZUsA6urqIuavqalBVeM+wlnwHThwgDFjxrB48WL69u3LokWLKC9PZtPw1GMKKg2o6gJV/bGq/g6oj5DmOVV9UVXXquoaVb0N2AMMTquwKcT88RnNmblz5/LOO+8cE75jxw6WL18ONJ6/6tSpE+B6SU0pqVRw8OBBvvnNb/Lqq6/Sp08fXnrpJU444QRf60yEnBziE5GLgeE4h7FnA22Aear6nSbydAbuBKqBE4EtwLPAdFXd5bvQjWVpAfwzrueVXbOpHrFslGtbbhjNmfnz5zNhwgQ6depEnz59KC8vZ8eOHbz++uvs27ePoUOHMm7cuIb0Xbt2pW/fvqxatYrevXvTv39/CgsL6dGjBzfffHNKZfvxj3/MSy+9BEDnzp350Y9+FDbd2LFjGTt2bNi4dJCTCgq4HaeY9uKGzppcii0ip+IUQQfgOeB9YABuSK5aRIakw5u5iJwFLAeKcLKPU9XjbDXkUWzLDaM5M3nyZCoqKli2bBkrV65k165dtGvXjn79+lFTU8Pll19OQUFBozzz589nypQpLFmyhCeffJIjR44wfPjwlCuonTt3Nvz9/PPPR0xXUVFhCsoHfohTTGtxPalXm07OTJxymqSqMwKBIvIrr6y7gWuDwu8i+v5QI1R1cZxyf4Dr9ZUD3wIeE5EqVf17nOX4RuiCw6awIT4jl4m2FnDo0KEMHTo0rjIrKip4+umnkxErJubMmRPRW0U2kZMKSlUbFFK0F6qIdANGA+uBB0Oi7wAmAuNFZHLQjrr3AU9EEePjOEQGQFUP45QqwEoR+SpOQX4v3rKygdYFrcmXfOtBGYaREDmpoOJkpHdepKqNDBhUdY+ILMUpsEHAy174dmB7GmTLAwrTUE/caFR/5u7joLSw1BSUYRgJYVZ80MM7r4kQ/6F3Pj3RCkSktYj0EZE+uDbv4v3uEpTmpyIyVEQqROQsEfkJUAXMi1DmRBFZKSIrt23blqhocROPmTmYPz7DMBLHFBSUeedIb9FAeDILBCqBVd7RCpju/X1nUJqOuGHDD3A9ta8C56vqi+EKVNVZqlqpqpXt27dPQjR/MXdHhmEkig3xRSfQaYjBsDo8nrFEk50PVa1JtPxMEIuZOThLvo17M7vpmWEYxyfWgzraQyqLEF8aks6Ig/LCcj4/aE1nGEb8mIJyQ2oQeY7pNO8caY6qWRGPN3M4OsQXy/YchmEYwZiCOrpGarSINGoPEWkDDMH5z1uRbsFygbLCMg7XH+ZAXSQXhIZhGOFp9gpKVdcBi3BexK8LiZ4OlABzg9ZAGcQ+IWeLdQ3DSJScNJIQkbFAwD9HR+88WETmeH9vV9WbgrL8AOfq6AERORd4DxgIjMAN7WV+a8ksQeI0NA/2x3dy65P9EMkwjBwlJxUUzl3QhJCwbt4BsAG3PxPgelEiUslRZ7EX4JzFPoBzFrsTIyHMH59hGImSkwpKVacB0+LM8wlwhR/y5CKx2jzYEJ9hGInS7OegjPiI24qvyLbcMAwjMUxBGb5S1tKG+AzDSAxTUIavFLQooKSgxIb4DCOLqK2t5fbbb+eCCy6gW7dulJaW0rJlS7785S9zySWX8MYbb2RaRMAUlJEgsXgzD2D++AwjOosXL0ZEqKqq8r2uzz77jLvvvps33niDjh07MmrUKMaMGUPr1q155plnGDp0KPfcc4/vckQjJ40kjOyitKVtuWEY2UTHjh1ZsWIFlZWVtGjRolHc008/zeWXX87UqVMZO3Ys3bt3z5CU1oMy0kB5YTm1h2ozLYZhGB6tW7dm4MCBxygngEsuuYThw4dz5MgRXnnllQxIdxRTUEZCxONaz4b4jFxERBp27J49ezYDBw6ktLQUEWH3bve8b968meuvv57u3btTVFREcXExXbp0obq6mlmzZjWUVVVVxYgRIwBYsmRJQ9npGvILJT/fDa4VFRWlve5GcmS09maEiAzDLQ7uD3wJuEJV54SkmYbbZj6YT1W1I1lCvGbm4BbrmoIycpUbbriBmTNnMmTIEMaMGcOaNWsQEbZs2UL//v3ZunUrXbt2pbq6msLCQjZt2sSKFStYv349EydOBKC6upqioiIWLlzISSedRHV1dUP5PXv2TOv1vPjii7z66qu0atWKkSNHRs/gI6ag0kdr4O/AXO+IxAe4nXQDHPFRprRQXlTOnsN7OFJ/hBZ5xw4pGMcfP/vzz3h/5/uZFiMuerbtyZQBU1Je7uOPP87y5csZMGBAo/D777+frVu3cs011/DQQw819LYADh06xJtvvtnwe+rUqQwaNIiFCxfSs2dP5syZE7auOXPmcMUV8fsTePTRR6mpqQkbN2XKFD799FP279/PmjVrePfdd2nTpg2PPfYYnTt3jruuVGIKKk2o6gJgAUCQT8Bw1Knq1rQIlSbKC8tRlNrDtZxQdEKmxTGMlHLLLbcco5wAPv30U8D1jiRk6KGwsJBhw4bFXVf37t2ZMCHUi1ts+SLx+9//nnXr1jX8bteuHb/+9a8ZN25c3PWkmpxVUCJyMTAc55fvbKANME9Vv9NEns4c9cd3Is4f37M4f3y7fBfa0U1ENgGHgTeBH6vq/6Wpbl8I9sdnCio38KMncrxy0UUXhQ0fMGAAM2fOZMoU11ajRo2ipKQkqbrOOecczjnnnKTKCGXt2rUA7N69m/fee4+f/exnXHzxxVx66aU88cQTYQ0p0kUuG0ncDlyPU1CboiUWkVOBt3H++P4M3Av8H3AjsFxETvRP1AbeBGqA84GrcZ7Yl6Wp7piI15s5mD8+I7fp2rVr2PDx48dz2WWXsWbNGsaNG0dZWRl9+vRh0qRJLFu2LM1SRqe8vJzBgwfz7LPP8o1vfIOnnnqKBx98MKMy5WwPCvghsBFYi+tJvdp0cmYCHYBJqjojECgiv/LKuhu4Nij8LqJvwzFCVRfHKrCqvhj8W0RW4JTkBOBXsZaTbQRvuWEYuUarVq3Chufl5TFv3jxuvfVWnn/+eZYuXcrSpUuZMWMGM2bM4Morr2T27Nlx1fXGG2/w8MMPxy3jVVddFVfPa8KECfzxj3/k97//PZMmTYq7vlSRswpKVRsUUuj4bygi0g0YDawHQj8Z7gAmAuNFZHLQxoX3AU9EEePjOEQ+BlXdKyKrObrtfNYQzxbutuWG0Zzp1asXvXr1AqC+vp4FCxZw2WWX8cgjj3DJJZcwevTomMtau3Ytjz32WNwyVFVVxaWg2rdvDziPE5kkl4f44iFgS7lIVeuDI1R1D7AUKAYGBYVvV9X3oxz7kxFKRIqAnri5sKwgETNzG+IzDEdeXh5jxozhwgsvBODdd99tiGvZsiUAdXV1EfPX1NSgqnEfkSz4IhFYoHvaaZn9NjYF5ejhnddEiP/QO5+eaAUi0lpE+ohIH1y7d/F+dwlKc4+IDBeRU0RkIPA73Jbzx3wyichEEVkpIiu3bduWqFhpoXVBa/Il33pQRrNi7ty5vPPOO8eE79ixg+XLlwON5686deoEuF5SU0oqVbK9+uqxsx5Hjhxh7ty5/OxnPwPg6quv9lWOaOTsEF+clHnnSJ/4gfDyJOqopPE82HTveAxnGAHQGXgSaAdsA1YAg1R1Q2hhqjoLmAVQWVkZh1+H1BCPJwkRobTQ/PEZzYv58+czYcIEOnXqRJ8+fSgvL2fHjh28/vrr7Nu3j6FDhzYy5e7atSt9+/Zl1apV9O7dm/79+1NYWEiPHj24+eabUyrbK6+80rDO6eyzz6a8vJxt27axevVqNm3aRF5eHtOmTeMb3/hGSuuNF1NQsREY2EpYEXjGEk0OkKnqpYmWn+2UF5bbEJ/RrJg8eTIVFRUsW7aMlStXsmvXLtq1a0e/fv2oqanh8ssvp6CgoFGe+fPnM2XKFJYsWcKTTz7JkSNHGD58eMoV1NVXX01paSnLly9n5cqV7Ny5k5YtW9KlSxcuuOACrr32Wvr165fSOhPBFJQj8OYsixBfGpKu2ZLAFBRg/viM3COaodDQoUMZOnRoXGVWVFTw9NNPJyNWTAwZMoQhQ4b4Xk+y2ByU4wPvHGmOKTBTGGmOyoiC+eMzDCNeTEE5AnNDo0WkUZuISBtgCHAANydkEP9YZ3lhOZ8fbPYdUMMw4sAUFKCq64BFQAVwXUj0dJwl3dygNVDNlmhryiIRGOKLZ/2UYRjNm5ydgxKRscBY72dgu4rBQY5at6vqTUFZfgAsAx4QkXOB94CBwAjc0F40rxFGE5QVlnG4/jAH6g5QXFCcaXEMwzgOyFkFhfPBF+r2t5t3AGzA7c8EuF6UiFRy1FnsBbgFsg/gnMXu9F3i44h4O0LBi3VNQRmGEQs5q6BUdRowLc48n+CcxRoRSMaKD5y7o5Nbn5w6gQzDyFlsDspIC+aP7/jD5guNWPHrWTEFZaSFgIKyxbrHBy1atODIkeN+M2cjTdTX15OXl3p1YgrKSAiN09DcHMYeXxQXF7N3795Mi2EcJ+zfvz/itiPJYArKiIsErcxtT6jjjNLSUnbu3Gm9KCMqqsru3buT3i04HKagjLRQ0KKA4vxiU1DHCW3atKGkpIQNGzawe/du6urqbE7KaISqcujQIbZs2UJdXR0nnHBCyuvIWSs+w18SeVeZw9jjBxGhQ4cO7Nmzh9raWj777DPrTRnHkJ+fT1lZGR06dPBlDsoUlBEXiXqSABhdMZqOJR2jJzSyAhGhtLSU0tLS6IkNwwdMQRlpY3Ll5EyLYBjGcYTNQaUBERkmIn8QkU0ioiJSEybNei8u9HghAyIbhmFkHFNQ6aE18HfgRpxX9HB8FTg56OiHcxr+TDoEjBebLjcMw29siC8NqOoCYAFAkLPa0DTbgn+LyPeAWuB//JbPMAwjG8nJHpSIXCwiM0TkdRGp9YbKnoiSp7OIPCIim0XkkDfkdp+IpN52MgriLBG+BzyhqvvTXb9hGEY2kKs9qNuBs4G9wEagZ1OJReRU3FYbHYDngPeBAbghuWoRGaKqO3yVuDGjgFOAh9NYZ3zYmhjDMHwmJ3tQwA9x27eXAt+PIf1MnHKapKpjVXWqqo4E7gV6AHcHJxaRuyIYNAQfVUnIfzXwlqr+JYkyfCMJS3PDMIyYyckelKoGtnCPum5HRLoBo4H1wIMh0XcAE4HxIjI5aEfd+4AmhwyBj+MQOVieDsCFHLuzr2EYRrMiJxVUnIz0zotUtT44QlX3iMhSnAIbBLzshW8HtvskTw1wCHjKp/JTgg3wGYbhN6ag3BAeuG3dw/EhTkGdjqeg4kVEWgPdvZ95QBcR6QPsVNWPg9IJcBXwlKruiVLmRFzvDmCviHyQgGjtSFDR3vTToO2IU0fC8viEydM0Jk/TmDyR6RpLIlNQUOadIzmJC4SXJ1FHJfBq0O/p3vEYrscUoAo4DfhOtAJVdRYwKwmZEJGVqlqZTBmpxORpGpOnaUyepsk2eWLBFFR0ApNYCY9qqepiYtgt3Zs7MxMEwzAMcteKLx4CPaSyCPGlIekMwzCMNGAKCgJzN6dHiD/NO0eaozqeSWqI0AdMnqYxeZrG5GmabJMnKpLrm5B565FeBeap6jFzO94i3bU4M/NTgy35RKQNsAWnyNsHmZkbhmEYPtPse1Cqug5YBFRw7Nqj6UAJMNeUk2EYRnrJyR6UiIwFxno/OwLnAf8HvO6FbVfVm4LSh7o6eg8YCIzADe19Lc2ujgzDMAxVzbkDmIazuot0rA+T58vAo7ghvcPABuB+oG2mryfOa/8y8DucUUctMB/oEmPeIuAXXhscAJYDw5KU52Lg9157HsDN+f0EaBND3kj3r08S8lRFKHN3htpncRPX+Sc/2wfoDMzwrmO/l7ci1deNc931Pm4B+gfAtYnKg1uyMcsrbz/OY8s84JQk2/tfk2ifZO9DKttnWhPyHExl+6TjyEkzc1WdhrtR8eT5BLjCD3nShYgUA6/gHvQJuAfrLuBVEemt0YcpZwNfB27G9TivAxaKyGBN3C/gTbiXyI9xjnv74u7NCBH5moZ47wjDHOC/Q8JSYbAyCXgr6HddDHn8aJ8fcNRSNMBg4FfAH2LIP4fE26c78C/A27jRhdER0iV83SJytSffT4CXgHOBmSIiqvpQAvJcCpwJPACsBjoB/wasFJE+3v9xNP4KXBMStj5MuljbBxK8Dz60z8PAn0LCSrywWJ4niL19/CcTWtEOfw6c9/UjQPegsFNwL98fRcl7Nk6hXREUlo/7ovtDEjK1DxP2Xa+ukVHyKnBXituoyiv3n+LM50v7RKhrNu4jo8nee7LtA+QF/X0V4b/IE75uL91nwGMh4Y/gPBoUJCBPuOepK1AP3BnDNS8G3khV+yRzH/xonwj1jPfSfj2V7ZOOo9kbSeQY3wRWqOraQICqfgQsxTmgjZb3C+DpoLx1OJ+A54lIYSICachGjB6BnkunRMrMEL60Tygi0gr4Z+CPqrozFWVGQqP3XiG56x4MtOdYx8qPAycC58QrT7jnSVU3ANtI8fMUY/skQ8rbJwITgE+BhQnmzximoHKLM3Fby4eyGjgjhrwf6bEbJK4GWnLUl2AqGO6d34sh7fe9DST3i8grIjI0RTLME5EjIrJDRH4rIl2ipE9X+1wEtMG5wYoFv9onQDLXfaZ3Dn0mV3vnaM9kTIjIV3AGTrE8TwB9ReRzEflCRP7q7V6dLIncB9/bR0Q644y95nkfFrHgR/skRE7OQTVj2gK7woTvBKLtDNxU3kB80ohIJ+BO4CVVXRkl+RPA88Bm3DDOzcArIjJKnfuoRPgc+CWwBGdE0hc3P7ZcRPqq6mcR8qWlfXDDn58BL8aQ1o/2CSWZ6w7EheZPWZuJSD7wa1wPanYMWV7DGVWswfnX/C7wsIicrKp3JShGovfB9/bBDe/lEfsHjx/tkzCmoHIPDRMWi38/SSJvTHhe3Z/DzYlFNUhR1fFBP18XkedwX5t3ETL8ESuqugpYFRS0REReA/6MM5y4PULWdLTPl4B/Au6P5WvXj/YJJxbJPVNEyJ8q/gv4Gm5+JZwibYSq/ntI0HMi8r/AbSJyn6rujVeAJO5DOtrnu8AqVf1rLIn9aJ9ksCG+3GIX4b+6TiD8V3AwO5vIG4hPGBEpwlkRdQPOU9WN8ZahbguSF4CvJiNLmHLfwX0xNlWur+3j8R3i+9pthE/tk8x1R+oJtA2JTwgR+Qluy5krVXVREkU9iTOlPysZeQLEcR/8bp8BQE8SfJ6CSGn7xIMpqNxiNUfHtYM5A/hHDHlP8UzVQ/MexrmDSggRKcCthRoAXKCqf0u0LCJ/0SdLtHJ9a58gvgu8q6rvJlFGqtsnmesOzKWEPpOBuZVoz2REROQ2YCpwo6o+nmg5geK8cyrbLZb74Fv7eEzAjVb8Nsly0tHTC4spqNziD8Agbxt7AESkAhhC9DUQfwAKcBZkgbz5wCW43YYPJSKQiOThxrTPBS5U1RWJlOOVVYpbj/NmomVEKLcS5yy4qXJ9aZ8QGc4kia9dn9onmetejjOXvjwk/Du43sHSRAQSkUm44bPbVHVGImWEcBluAXIyH04NxHEffGkfT4aWuDVjCyJY0sZDStsnHmwOKrf4DXA9btz4dtwXz38AnxC0iFBEugLrcOtG7gRQ1b+IyNPAfV6P5yPg+7h1VKH/QPHwIO7ldjewT0QGBcVtVNWN4eQRkZtwux2/ytHJ55twrqsSlkdE5uGu7R1gN85I4lZgE26VfrrbJ8B3ifC162f7iMjF3p/9vfP5IrIN2KaqS+K5bhFZC2xQ1XMBVPULEfk33MLTTbiFqCOBK4EbVPVwvPKIyKXAfbiFp6+EPE+1qvqPoLIayeNZ1k3FeVdZj9tiZwLOlH6qhlnIHoM8Md+HdLRPUNIxuKHCiB88qWgf38n0Qiw7UnsAXXDDabXAHuBZjl3sWIFTXtNCwlvhPBhsBQ7ivgCrkpRnPZFdr0yLJA/wDdwX5HbcOpwduK/5AUnKcytupfznXrmf4FznnJyJ9vHKLcBZof0xQrxv7dPEvVkc73V793pxmPBrcHN8h4APgR8kKg/OY0NUmcPJgzOJfxH3MXII2IvzwfntJOSJ+T6ko32C0j3nydIyyv9mUu3j95GTzmINwzCM4x+bgzIMwzCyElNQhmEYRlZiCsowDMPISkxBGYZhGFmJKSjDMAwjKzEFZRiGYWQlpqCMnEFEFotITq2bEJHTROR/RWSriKiI7M60TNmEiFR47TIn07Kkmlx8nuPFFJTRCO+fXUVkg+fgNVya9V4a80TiIyLSArfQ+gLcdg7TgZ9mVCiOD6UgInM8GSsyLUskjgcZM429YIxIdAH+lSx4ITZjTsE5Dv2Nqk7MtDBZyibgKzjPILnGd4FQJ73NCutBGeHYhXNWeauItMu0MM2YL3nnzRmVIotR1S9U9X1V3ZJpWVKNqn6squ9nWo6MkikfS3Zk54Hz67UR13tSYEaYNOu9uPygsCrC+K8LybM+JKzGy1MDjAJex/n/2gY8CpR76frihrh2efF/IMS/oJdusVdeIc7b9Uc4n2LrgDuI4JcMt2fOHJxfvkPApzinrT3CpJ3j1dENuAHn1+8AYXysRairP85X4mdeXRuAmQT5Agy6DxH9F8ZQz7/gdkf9nKOeqG8FCiPc87DyB11vhfd7WhOy1YTkPQ9YgPNTF7gPvwjc13DPB1CK8/u3Hufbrsnr5aifwjkxtF3o89cW+Aluq/gDXlu9DIwOU08NR5/Vau9Z+xzQoDRjcbvrrgH24Z7Vt3EbYebFeH/XB6VZHFx+UHgecC3wllfHPu/v74fWE3x/gXY4v5NbvPuxGrgiTHrBOYldhvtfPIj731gIXJKO91DgsCE+IxIP4jyjXyMiM1R1jY91fRPnffl53PbdX8O9CE4Rkam4l8bruC29z8I56DxVRM5S1fow5T2D2yzud7iX3IW4F2uliHxTvf9CABGpxnlvLgD+iNvfqDNwEfB1ERmhbkPDUO4HhuI2plsAHIl2kSIyBqecxJNtA05hfR+4UESGqOp6L/l03Mt3Am57+sVe+GKiICL/iVNG23GKdi9wPvCfwHneVuRfRCsnAotxW4HfCLyLmyML8JcgGf7du4aduPv6GdAb5+n7AhEZrKq1IWW3BF7BKY5FOIfHHyUg43Scsjgbd58ChiUNBiaeh/jFuDZ+HecdvQT3HP5JRK5R1d+EKftinIJ6EfesVgTF/RSoxznT3YTzBj7Sk+GruO3XY5axCR7HbYHxCfAwTgGNw33onEN4b/blOKe2h3HPXpF3LY+ISL2qBns9vxv3/HyE+1/6HDjZu4Z/Bp6OQcbUkE5taEf2H3g9KO/vi73f80PSrCe1Pag6YHhQeB7w/7y4ncDlIflme3EXhoQv9sLXACcEhRfh9t5RYHxQeGCn4e3AGSFlnYl7sb8TEj7HK2cTcEoc7draq+cIMDQkbopX5qKQ8CbbNEI9g708HwMdg8LzcQpYgR+HueeLI5QXuN6KoLAKQnotIXlGePHLCOktBd3zeyM8Uy8BJXFcb1hZwskd5lmpBy4NCS/HKdoDwElh5K4HqiOUeWqYsMAOyQoMTEBGDQn7tpfnHaB1UHgJsNKLuyzM/VWcMmsRFH4G7n/vHyHpd+BGUYrDyNQu1nuTisPmoIyIqOrvcC/2cSJyjo9VPalBe9mo6xUFdkn9u6rOC0k/1zv3iVDef6hqwxb3qnoQ90UIbq+dAN/FvZDu0KB9hLw8q3H7a/UVkTM4lp+rajxf9xcCJwJPq+rrIXG/xL2gR4lIlzjKDEfg+u5S1a2BQFWtAybjXrBXJVlHNCZ556tVtVGPQFXn4BRApD2rJqvP+w6JyNnAcOD3qvpUiHy7ccPBRcC3wmR/TlX/FK5cVV0XJqwe10MCN+SZLIH7O1VV9wbVsw/3oQPh7+9+4EeqeiQozz9wvaqviEibkPRfEGZUQFW3JyF73NgQnxGNybgv4V+KyCD1PqNSzMowYQHDgLfDxG3yzp0jlLckTNjruJFG1qgAAAXNSURBVK/FvkFhg73z2SIyLUye073zVzh2++0/R6g7Ev288yuhEapaJyKv4XoDfXG9n0Rpqp41IrIRN3RaHqo8Ushg3Avun0Xkn8PEtwTai8iJqrojKPwgbk7PbwL3vSzCfW/vnb8SJi7ifReRE4GbccsCuuF6NcF0ik/MsPTDfWQsDhO3BKdU+oaJ+1CPHVIFN0wI7kNtj/f3PNz86moR+R+v3OWqmnZLSVNQRpOo6nIR+R1uuO9f8Gf8OdyDXxdDXEGE8j4NDVDVIyKyA+gQFHyid746inytw4RtDRPWFGXeOZK1WSC8PM5yE6mni5fOLwV1Iu7dckeUdK1xw0kBPvPpAyiUwH0f5R2RiPm+i0g5zlDhFJwSm4sbnq7j6JxdYYLyBlMG7NQwu+16HzrbafyMB4h0rwP/Sy2Cwn6IM2i5ErfD7lSgTkQW4Hq4axMVPl5MQRmxMBU3RPUTEfnfCGkCxgqRnqky0rdW5SRCeiHeotcTcRPvAQLynK2q8X65x/siDdTVMUL8ySHpEiW4nmOGnCLUo0S+b4kozM9x1mRt48yXDuUER6/9RlV9IM68kWS8CqecpqvqtOAIERmMU1Cp4HOgrYgUaIihi7dwvh2Nn/G48YYB7wfuF5EOOMOLS3EGEmeKyJmqeiiZOmLF5qCMqHhj6zNx/4A3REgWmPP5cmiEiHQn+Z5BPAwPEzYU9xJeFRS2IijObwL1VoVGeC+WwBxfOIvBVNXTHTcs+lHI8N4uwt+3FoSf5wvMTbQIEweuXU8QkTNjlNkPmpLRj/ve3Tv/PkxcuOcRordjOFbh3tvDwsQN88pK9hlqQFU/U9X5qvovuGHjU4FeqSo/GqagjFi5EzdMcBvhhz7ex325Xeh9dQEgIq2AeL9Sk+XfROSEIBmKcOtdwK2vIujv3cAdIjIgtBARyRORqhTJ9CxuyOfbIjIoJO5fcXMWL6lqMvNPAI9459tFJDCXElA29+D+52eH5Pkz0EVERoeE3w50DVPHLlxPIpJBx73e+Tci8qXQSBEpCdMGqSYwdHiMjKq6EjcneZGIXBkaDyAiZwU/xzGw3jtXhZTTl6MGOjHL2ASB+/sTEWnwMuH9HfD6Enp/Y0ZECkXkXBGRkPACnPk/OIOLtGBDfEZMqOpOb33NzyPEfyEi9wP/BqzyhgLzcWP8m0mvN4T3cBO8weugTsWtWQpYB6KqO0TkYuB/gRUi8jJu8WI97qUxGDcsGNYnYTyo6l7vZfg/wBJv8vlj3Dqo0bi5jWtSUM8yEfk5cAvwd68N9uHWQfUC3sAtlg3mHpyF2XMi8jROkX4N12NeTMhL17uWN4GhIjIPZ9Z/BPiDqv5VVV/21q/9BPjQm7v4CPdh0xXXo3gDt57IL17GGSz8xmuDvcBuVf0vL/4yXI9gtohMwq1d2o3rYfbGtdVg3PqtWJjr1XefiIwAPgROw62rmg9ckoCMx6CqvxWRC3HzwatF5Fncx8JY3P16JozVazy0wpn6r/fu8Qbc8z8KZzTyB1V9L4ny4yOdNu12ZP9B0DqoMHGFuBdNYF1Ffki84Oar1uEWBH6MU2jFRPEkEaauKiKsASLy2pfFhPck8X+4CftjvCgElfdfuJfKQVxP8H2cMhsbknYOTaxdiaF9v4pTiNuC2ugh4EvxtEEM9VyKUwJ7vGtajev9FkVI/02cNeVB3Jf9UzhlEvZ6cUNaf/TS1oe7j7hhy2dwHyeHvWv+C85TRGVI2mOejxivM+yz4MX9CPexcojwniTaAD/GWYruxa19+gj3ITORoPVYTT2rQWnOwHk5+Qz3UfA2bm4qIRlp2pPED7z7td873gauowlPEhFkbnR/cYZHt+AWIn/sPQ/bcMOi1xLBG4tfh3hCGYZhGEZWYXNQhmEYRlZiCsowDMPISkxBGYZhGFmJKSjDMAwjKzEFZRiGYWQlpqAMwzCMrMQUlGEYhpGVmIIyDMMwshJTUIZhGEZW8v8BMO16pIJQH54AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.sparse.linalg as la\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "# Example from http://www.caam.rice.edu/~embree/39961.pdf\n",
    "\n",
    "A = np.array([[1, 1, 1],\n",
    "              [0, 1, 3],\n",
    "              [0, 0, 1]]\n",
    "            )\n",
    "rhs = np.array([2, -4, 1])\n",
    "x0 = np.zeros(3)\n",
    "\n",
    "for restart in [1, 2, 3]:\n",
    "    hist = []\n",
    "    def callback(rk):\n",
    "        hist.append(np.linalg.norm(rk)/np.linalg.norm(rhs))\n",
    "    _ = la.gmres(A, rhs, x0=x0, maxiter=20, restart=restart, callback=callback)\n",
    "    plt.semilogy(np.array(hist), label='rst={}'.format(restart))\n",
    "plt.legend(fontsize=22)\n",
    "plt.xlabel(\"Number of outer iterations\", fontsize=20)\n",
    "plt.ylabel(r\"$\\frac{||r_k||_2}{||rhs||_2}$\", fontsize=20)\n",
    "plt.xticks(fontsize=16)\n",
    "plt.yticks(fontsize=20)\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "\n",
    "- Iterative methods continued (BiCG, MINRES, QMR), preconditioners.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.9"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}