lectures/lecture5/lecture-5.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 5: Linear systems"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Previous lecture\n",
    "- Matrix rank\n",
    "- Skeleton decomposition\n",
    "- Low-rank approximation\n",
    "- Singular Value Decomposition (SVD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture\n",
    "\n",
    "- Linear systems\n",
    "- Inverse matrix\n",
    "- Condition number\n",
    "- Gaussian elimination"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear systems\n",
    "\n",
    "- Linear systems of equations are the basic tool in NLA.\n",
    "\n",
    "- They appear as:\n",
    "\n",
    "    - Linear regression problems\n",
    "    - Discretization of partial differential/integral equations\n",
    "    - Linearizations of nonlinear regression problems\n",
    "    - Optimization (i.e., Gauss-Newton and Newton-Raphson methods, KKT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear equations and matrices\n",
    "- From school we know about linear equations. \n",
    "\n",
    "- A linear system of equations can be written in the form\n",
    "\n",
    "\\begin{align*}\n",
    "    &2 x + 3 y = 5\\quad &\\longrightarrow \\quad &2x + 3 y + 0 z = 5\\\\\n",
    "    &2 x + 3z = 5\\quad &\\longrightarrow\\quad &2 x + 0 y + 3 z = 5\\\\\n",
    "    &x + y = 2\\quad &\\longrightarrow\\quad  & 1 x + 1 y + 0 z = 2\\\\\n",
    "\\end{align*}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrix form\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "2 & 3 & 0 \\\\\n",
    "2 & 0 & 3 \\\\\n",
    "1 & 1 & 0 \\\\\n",
    "\\end{pmatrix}\\begin{pmatrix}\n",
    "x \\\\\n",
    "y \\\\\n",
    "z \n",
    "\\end{pmatrix} = \n",
    "\\begin{pmatrix}\n",
    "5 \\\\\n",
    "5 \\\\\n",
    "2\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "\n",
    "or simply\n",
    "\n",
    "$$ A u = f,  $$\n",
    "\n",
    "where $A$ is a $3 \\times 3$ matrix and $f$ is right-hand side"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Over/under determined linear systems\n",
    "\n",
    "If the system $Au = f$ has\n",
    "\n",
    "-  more equations than unknowns it is called **overdetermined** system (generically, no solution)\n",
    "\n",
    "- less equations than unknowns it is called **underdetermined** system (solution is non-unique, to make it unique additional assumptions have to be made)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence of solutions \n",
    "\n",
    "A solution to the linear system of equations with a **square** matrix $A$\n",
    "\n",
    "$$A u = f$$\n",
    "\n",
    "exists, iff \n",
    "* $\\det A \\ne 0$\n",
    "\n",
    "or\n",
    "\n",
    "* matrix $A$ has full rank."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Scales of linear systems\n",
    "\n",
    "In different applications, the typical size of the linear systems can be different. \n",
    "\n",
    "- Small: $n \\leq 10^4$ (full matrix can be stored in memory, **dense matrix**)\n",
    "- Medium: $n = 10^4 - 10^6$ (typically, **sparse** or **structured** matrix)\n",
    "- Large: $n = 10^8 - 10^9$ (typically **sparse** matrix + parallel computations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear systems can be big\n",
    "- We take a continious problem, discretize it on a mesh with $N$ elements and get a linear system with $N\\times N$ matrix.  \n",
    "- Example of a mesh around A319 aircraft\n",
    "(taken from [GMSH website](http://geuz.org/gmsh/)).  \n",
    "<img src=\"a319_4.png\" width=50%>\n",
    "\n",
    "The main difficulty is that these systems are big: millions or billions of unknowns!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Linear systems can be structured\n",
    "\n",
    "- Storing $N^2$ elements of a matrix is prohibitive even for $N = 100000$.  \n",
    "\n",
    "**Q:** how to work with such matrices?  \n",
    "\n",
    "**A:** fortunately, those matrices are **structured** and require $\\mathcal{O}(N)$ parameters to be stored.  \n",
    "\n",
    "- The most widespread structure are **sparse matrices**: such matrices have only $\\mathcal{O}(N)$ non-zeros!  \n",
    "\n",
    "- Example (one of the famous matrices around for $n = 5$):\n",
    "\n",
    "$$\n",
    "  \\begin{pmatrix}\n",
    "  2 & -1 & 0 & 0 & 0 \\\\\n",
    "  -1 & 2 & -1 & 0 & 0 \\\\\n",
    "  0 & -1 & 2 & -1 & 0 \\\\\n",
    "  0 & 0 &-1& 2 & -1  \\\\\n",
    "  0 & 0 & 0 & -1 & 2 \\\\\n",
    "  \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "- At least you can store such matrices\n",
    "- Also you can multiply such matrix by vector fast\n",
    "- But how to solve linear systems?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Main questions about linear systems\n",
    "\n",
    "1. What is the accuracy we get from the solution (due to rounding errors)?\n",
    "2. How we compute the solution? (LU-decomposition, Gaussian elimination)\n",
    "3. What is the complexity of the solution of linear systems?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to solve linear systems?\n",
    "**Important**: forget about determinants and the **Cramer rule** (it is  good for $2 \\times 2$ matrices still)!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to solve linear systems?\n",
    "\n",
    "The main tool is variable elimination. \n",
    "\n",
    "\\begin{align*}\n",
    "    &2 y + 3 x = 5 \\quad&\\longrightarrow \\quad &y = 5/2 -  3/2 x \\\\\n",
    "    &2 x + 3z = 5 \\quad&\\longrightarrow\\quad &z = 5/3 - 2/3 x\\\\\n",
    "    &z + y = 2 \\quad&\\longrightarrow\\quad  & 5/2 + 5/3 - (3/2 + 2/3) x = 2,\\\\\n",
    "\\end{align*}\n",
    "\n",
    "and that is how you find $x$ (and all previous ones).  \n",
    "\n",
    "This process is called **Gaussian elimination** and is one of the most widely used algorithms. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian elimination\n",
    "Gaussian elimination consists of two steps:\n",
    "1. Forward step\n",
    "2. Backward step"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Forward step\n",
    "- In the forward step, we eliminate $x_1$:\n",
    "\n",
    "$$\n",
    "   x_1 = f_1 - (a_{12} x_2 + \\ldots + a_{1n} x_n)/a_{11},\n",
    "$$\n",
    "\n",
    "and then substitute this into the equations $2, \\ldots, n$. \n",
    "\n",
    "- Then we eliminate $x_2$ and so on from the second equation. \n",
    "\n",
    "- The important thing is that the **pivots** (that we divide over) are not equal to $0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Backward step\n",
    "In the backward step:\n",
    "- solve equation for $x_n$\n",
    "- put it into the equation for $x_{n-1}$ and so on, until we \n",
    "compute all $x_i, i=1,\\ldots, n$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Gaussian elimination and LU decomposition\n",
    "\n",
    "- Gaussian elimination is the computation of one of the most important matrix decompositions: **LU-decomposition**.\n",
    "\n",
    "**Definition**: LU-decomposition of the square matrix $A$ is the representation\n",
    "\n",
    "$$A =  LU,$$\n",
    "\n",
    "where \n",
    "- $L$ is **lower triangular** (elements strictly above the diagonal are zero)\n",
    "- $U$ is **upper triangular** matrix (elements strictly below the diagonal are zero)\n",
    "\n",
    "This factorization is **non-unique**, so it is typical to require that the matrix $L$ has ones on the diagonal."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Main goal** of the LU decomposition is to solve linear system, because\n",
    "\n",
    "$$ A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, $$\n",
    "\n",
    "and this reduces to the solution of two linear systems **forward step**\n",
    "\n",
    "$$ L y = f, $$\n",
    "\n",
    "and **backward step**\n",
    "\n",
    "$$ U x = y. $$\n",
    "\n",
    "Does $LU$ decomposition always exist?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Complexity of the Gaussian elimination/LU decomposition\n",
    "\n",
    "- Each elimination step requires $\\mathcal{O}(n^2)$ operations. \n",
    "\n",
    "- Thus, the cost of the naive algorithm is $\\mathcal{O}(n^3)$.  \n",
    "\n",
    "**Think a little bit**: can Strassen algorithm help here? "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Block LU-decomposition\n",
    "\n",
    "We can try to compute **block** version of LU-decomposition:\n",
    "\n",
    "\n",
    "$$\\begin{pmatrix} A_{11} & A_{12} \\\\\n",
    "A_{21} & A_{22}\n",
    "\\end{pmatrix} = \\begin{pmatrix} L_{11} & 0 \\\\\n",
    "L_{21} & L_{22}\n",
    "\\end{pmatrix} \\begin{pmatrix} U_{11} & U_{12} \\\\\n",
    "0 & U_{22} \n",
    "\\end{pmatrix} $$\n",
    "\n",
    "- There are two basic operations: compute LU-factorization of half-matrices + matrix-by-matrix product."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Existence of the LU-decomposition\n",
    "- The LU-decomposition algorithm does not fail if **we do not divide by zero** at every step of the Gaussian elimination.\n",
    "\n",
    "**Q:** when it is so, for which class of matrices?\n",
    "\n",
    "**A:** it is true for **strictly regular matrices**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Strictly regular matrices and LU decomposition \n",
    "\n",
    "- **Definition.** A matrix $A$ is called *strictly regular*, if all of its **leading principal minors** (i.e, submatrices in the first $k$ rows and $k$ columns) are non-singular. \n",
    "\n",
    "- In this case, there always exists an LU-decomposition. The reverse is also true (check!)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Corollary:** If $L$ is unit triangular (ones on the diagonal), then $LU$-decomposition is unique. <br>\n",
    "\n",
    "**Proof:** Indeed, $L_1 U_1 = L_2 U_2$ means $L_2^{-1} L_1 = U_2 U_1^{-1}$. $L_2^{-1} L_1 $ is lower triangular with ones on the diagonal. $U_2 U_1^{-1}$ is upper triangular. Thus, $L_2^{-1} L_1 = U_2 U_1^{-1} = I$ and $L_1 = L_2$, $U_1 = U_2$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LU for positive definite Hermitian matrices (Cholesky factorization)\n",
    "\n",
    "- Strictly regular matrices have LU-decomposition. \n",
    "\n",
    "- An important **subclass** of strictly regular matrices is the class of **Hermitian positive definite matrices**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "\n",
    "**Definition.** A matrix $A$ is called <font color='red'> positive definite </font> if for any $x: \\Vert x \\Vert \\ne 0$ we have\n",
    "\n",
    "$$\n",
    "(x, Ax) > 0.\n",
    "$$\n",
    "- if this holds for $x \\in \\mathbb{C}^n$, then the matrix $A$ has to be hermitian\n",
    "- if this holds for $x \\in \\mathbb{R}^n$, then the matrix $A$ can be non symmetric"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- **Claim:** A Hermitian positive definite matrix $A$ is strictly regular and has **Cholesky factorization** of the form \n",
    "\n",
    "$$A = RR^*,$$\n",
    "\n",
    "where $R$ is a lower triangular matrix.\n",
    "\n",
    "- Let us try to prove this fact (on the whiteboard).\n",
    "\n",
    "- It is sometimes referred to as \"square root\" of the matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Computing LU-decomposition\n",
    "\n",
    "- In many cases, computing LU-decomposition once is a good idea!\n",
    "\n",
    "- Once the decomposition is found (it costs $\\mathcal{O}(n^3)$ operations), then solving linear systems with $L$ and $U$ costs only $\\mathcal{O}(n^2)$ operations.\n",
    "\n",
    "**Check:**\n",
    "\n",
    "- Solving linear systems with triangular matrices is easy (why?). \n",
    "- How we compute the $L$ and $U$ factors?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## When LU fails\n",
    "\n",
    "- What happens, if the matrix is not strictly regular (or the **pivots** in the Gaussian elimination are really small?). \n",
    "\n",
    "- There is classical $2 \\times 2$ example of a matrix with a bad LU decomposition.  \n",
    "\n",
    "- The matrix we look at is  \n",
    "\n",
    "$$\n",
    "    A = \\begin{pmatrix}\n",
    "    \\varepsilon & 1 \\\\\n",
    "    1 & 1 \n",
    "    \\end{pmatrix}\n",
    "$$\n",
    "\n",
    "- If $\\varepsilon$ is sufficiently small, we **might** fail. In contrast the Cholesky factorization is **always stable**.\n",
    "\n",
    "Let us do some demo here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "L * U - A:\n",
      " [[0. 0.]\n",
      " [0. 0.]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray([[1.e+00, 0.e+00],\n",
       "             [1.e+04, 1.e+00]], dtype=float64)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "import jax\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "\n",
    "eps = 1e-4#1.12e-16\n",
    "a = [[eps, 1],[1.0,  1]]\n",
    "a = jnp.array(a)\n",
    "a0 = a.copy()\n",
    "n = a.shape[0]\n",
    "L = jnp.zeros((n, n))\n",
    "U = jnp.zeros((n, n))\n",
    "for k in range(n): #Eliminate one row   \n",
    "    L = jax.ops.index_update(L, jax.ops.index[k, k], 1)\n",
    "    for i in range(k+1, n):\n",
    "        L = jax.ops.index_update(L, jax.ops.index[i, k], a[i, k] / a[k, k])\n",
    "        for j in range(k+1, n):\n",
    "            a = jax.ops.index_add(a, jax.ops.index[i, j], -L[i, k] * a[k, j])\n",
    "    for j in range(k, n):\n",
    "        U = jax.ops.index_update(U, jax.ops.index[k, j], a[k, j])\n",
    "\n",
    "print('L * U - A:\\n', jnp.dot(L, U) - a0)\n",
    "L"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## The concept of pivoting\n",
    "\n",
    "- We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over.  \n",
    "\n",
    "- The simplest but effective strategy is the **row pivoting**: at each step, select the index that is maximal in modulus, and put it onto the diagonal. \n",
    "\n",
    "- It gives us the decomposition \n",
    "\n",
    "$$A = P L U,$$\n",
    "\n",
    "where $P$ is a **permutation matrix**.\n",
    "\n",
    "\n",
    "**Q.** What makes **row pivoting** good?  \n",
    "\n",
    "**A.** It is made good by the fact that  \n",
    "\n",
    "$$ | L_{ij}|<1, $$\n",
    "\n",
    "but the elements of $U$ can grow, up to $2^n$! (in practice, this is very rarely encountered).\n",
    "\n",
    "- Can you come up with a matrix where the elements of $U$ grow as much as possible?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Stability of linear systems\n",
    "\n",
    "- There is a fundamental problem of solving linear systems which is independent on the algorithm used.\n",
    "\n",
    "- It occures when elements of a matrix are represented as floating point numbers or there is some measurement noise.\n",
    "\n",
    "Let us illustrate this issue on the following example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "11.1598010212843\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEDCAYAAAA4FgP0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABR20lEQVR4nO29e3hcd3nv+3nnKs1oJNm62bKTyHF8iZ3ECTG5EYhcSAm07AAbWnJoS3fZO5uzyz5tn7YHTvuc7p52d7dP2S277FJoaHlKSwuFFkoIISEXK4mTEGInjhPHduzYji+yJY3uM6O5/84fa9ZoJM1trVkjzUi/z/PosTSzlub1SFrvem/fV5RSaDQajWbt4lppAzQajUazsmhHoNFoNGsc7Qg0Go1mjaMdgUaj0axxtCPQaDSaNY52BBqNRrPGaVpHICJfFZFREXmtimPfJSIviUhaRD6y6Lk/FZGjInJMRL4gIlI/qzUajabxaFpHAPwdcE+Vx54Dfhn4p8IHReQO4B3ADcB1wNuBuxyzUKPRaJqApnUESqmngYnCx0Rkq4g8IiKHROQZEdmZO/asUuoIkF38bYAWwAf4AS8wUn/rNRqNpnFoWkdQggeA/6qUuhn4LeCvyh2slHoe2A9cyn08qpQ6VncrNRqNpoHwrLQBTiEibcAdwLcL0vz+CudcA1wLbM499JiIvCsXbWg0Gs2aYNU4AozoZkopdaOFcz4E/FgpFQEQkR8CtwHaEWg0mjXDqkkNKaVmgDMi8lEAMdhT4bRzwF0i4hERL0ahWKeGNBrNmqJpHYGIfAN4HtghIhdE5JPAx4FPisgrwFHg3tyxbxeRC8BHgb8WkaO5b/MvwJvAq8ArwCtKqe8v839Fo9FoVhTRMtQajUaztmnaiECj0Wg0ztCUxeLu7m41MDBg69xoNEowGHTWoGWkme1vZtuhue1vZttB2+8Uhw4dCiulehY/3pSOYGBggIMHD9o6d2hoiMHBQWcNWkaa2f5mth2a2/5mth20/U4hIm8Ve1ynhjQajWaNox2BRqPRrHG0I9BoNJo1jnYEGo1Gs8bRjkCj0WjWONoRaDQazRpHOwKNRqNZ42hHoNFoqiaeyvDtg+fR0jSrC+0INBpN1fzo9RF++1+O8MZIZKVN0TiIdgQajaZqpmNJAGbjqRW2ROMk2hFoNJqqmYmnAYgk0itsicZJtCPQaDRVM5tzBLFkZoUt0TiJdgQajaZqZnIpIR0RrC60I9BoNFWTjwi0I1hVaEeg0WiqxiwSR3VqaFWhHYFGo6kaMyKI6ohgVaEdgUajqZqZuVxEoB3BqsIRRyAiXxWRURF5rcTzIiJfEJFTInJERN5W8Nw9InIi99xnnbBHo9HUh3xEoFNDqwqnIoK/A+4p8/z7gG25j/uBLwGIiBv4Yu75XcB9IrLLIZs0Go3D5GsEOiJYVTjiCJRSTwMTZQ65F/h7ZfBjoFNENgK3AKeUUqeVUkngm7ljNRpNg5HOZPORgI4IVhfLtbx+E3C+4OsLuceKPX5rsW8gIvdjRBP09fUxNDRky5BIJGL73Eagme1vZtuhue13wvZIcl5obnh0fFnfi2Z+76Hx7V8uRyBFHlNlHl/6oFIPAA8A7N27Vw0ODtoyZGhoCLvnNgLNbH8z2w7Nbb8Ttp+fiMGT+wHwtAQZHHyXA5ZVRzO/99D49i+XI7gAXFHw9WZgGPCVeFyj0TQY5lRxwOcmmtQ1gtXEcrWPPgj8Uq576DZgWil1CXgR2CYiW0TEB3wsd6xGo2kwZuaMi/+GjhaiCV0jWE04EhGIyDeAQaBbRC4A/w3wAiilvgw8DLwfOAXEgP+Qey4tIp8GHgXcwFeVUkedsEmj0TiL2TG0ob2Fi5NzK2yNxkkccQRKqfsqPK+AXy3x3MMYjkKj0TQw5gzBho4WEuks6UwWj1vPpK4G9E9Ro1lmRmfjTdmHb0YEGztaAN1CuprQjkCjWWbue+DH/NmP3lhpMyxjLqXZ0J5zBE3ozDTF0Y5Ao1lGlFKcn5zj4lRspU2xzGw8RYvXRUfAB0BMdw6tGrQj0GiWkVgyQzKdZSrWfDt/Z+Np2lu8BH1uACK6c2jVoB2BRrOMTOaWv0/PNZ8jmImnCLV4CPqNHhO9nGb1oB2BRrOMTEYNB9CsEUGoxUtbzhHodZWrB+0INJplxIwIpuaSK2yJdWbiadpbvQRyqSG9wH71oB2BRrOMmI4gnsoSTzXXhXQ2lxrSEcHqQzsCjWYZmYzORwIzTVYnmJlL097iIWDWCHTX0KpBOwKNZhmZLKgNTDWZIzAiAi8Br+4aWm1oR6DRLCNmagiaq2CcTGdJpLO0t3hwuYSAz627hlYR2hFoNMvIgogg1jwFY1NeItTiBSDo92gp6lWEdgQazTIyFUvS3eY3Pm+i1JApLxFqMeoDQZ9bp4ZWEdoRaDTLyEQ0yZbuAADTTZQaMiOC9oKIQKeGVg/aEWg0y8hULMUV6wO4XdJU08WzSyICj24fXUVoR6DRLCMT0STrAz46Wr1NNVRmtrrO1wjceqBsFeGIIxCRe0TkhIicEpHPFnn+t0XkcO7jNRHJiMj63HNnReTV3HMHnbBHo2lE4qkMc6kM64I+Olu9TdU1tDgiCPg9WoZ6FVHzhjIRcQNfBO7GWFL/oog8qJR63TxGKfU54HO54z8A/IZSaqLg2+xTSoVrtUWjaWTMC/+6gI/2Vm9TpYbMxfXtrUZE0ObTXUOrCScigluAU0qp00qpJPBN4N4yx98HfMOB19VomoqJ3FTxuoCXzkBzRgSmvETA79YL7FcRTuws3gScL/j6AnBrsQNFJADcA3y64GEF/EhEFPDXSqkHSpx7P3A/QF9fH0NDQ7aMjUQits9tBJrZ/ma2HWq3//Vx48L51snXScykuDyVXbb3o2bbTyVoccMzTz8FQPhykmgizf79+xERh6wszVr/3ak3TjiCYr8FqsSxHwCeXZQWeodSalhEeoHHROS4UurpJd/QcBAPAOzdu1cNDg7aMnZoaAi75zYCzWx/M9sOtdsfOTIML77Mvjtu4fJPznH0pQvL9n7UavtDY6+wbiqc/x7H5U2+/+Zxbn3HOwn4nLiMlGet/+7UGydSQxeAKwq+3gwMlzj2YyxKCymlhnP/jgLfxUg1aTSrDnOqeF3QS0erl5l4mky21D1TYzEbT+VnCID8ljKdHlodOOEIXgS2icgWEfFhXOwfXHyQiHQAdwHfK3gsKCIh83Pgp4HXHLBJo2k4pnI1gs5WH50B46LaLAqkM3PpfMcQkN9SpjuHVgc1x3RKqbSIfBp4FHADX1VKHRWRT+We/3Lu0A8BP1JKRQtO7wO+m8sxeoB/Uko9UqtNGk0jMhFL0ub34PO46Mh130zNpVgX9K2wZZWZTaToyUljAPl0kO4cWh04ktxTSj0MPLzosS8v+vrvgL9b9NhpYI8TNmg0jc5ULMW6oOEAzIjAEJ4LrqBV1TEbT7O1py3/dVs+ItCpodWAnizWNB0vnp3gU/9wqGny6yYT0STrAsbdf0er8W+zzBIY+4rn7xsDfrNGoCOC1YB2BJqm46+fOs0jRy8zOhtfaVMsMRWbdwRmRNAMjkApxcxcKi8vAQURgU4NrQq0I9A0FdOxFE+9MQrA2Gxiha2xxmQsxbqcA+g0awRNMFQWT2VJZ9XCiMCnI4LVhHYEmqbi0dcvk8oYKaFwpMkcQTRJZy4iaG8iR7BYghp0jWC1oR2Bpql46MglQrmLUDNFBKlMltlEmvW5DiGv20Wb39MUCqSLl9JAQdeQjghWBdoRaJqG8UiCZ0+F+bm3G/OLzeQIzF3FZmoIoKPV2xTLaWaKRAQ+jwuf20VUS1GvCrQj0DQNP3ztMpms4iM3bybU4iEcafy7aZO88mjBzEBnoDkUSE3BufbWhd3mhvCcjghWA9oRaJqG778yzNaeIDs3hOgJ+ZsrIsgrjy50BM2wt3jx4nqToJaiXjVoR6BpCi5Px/nJ2Qk+sKcfEaGnrckcQS411FmQGups9eUGyhqbmbmlNQIwCsbLERH83vde4zsnG/99ambqLxuo0TjAD169hFLwszf0A9Ad8nNseGaFraoeU3BufUFqqFmW05SKCALLtK5y6MQYybiOPOqJjgg0TcFDR4bZtbGda3oNmYNmiwgmSqWGYimUauwJ6dl4GpfMK46atPmXZ4F9OJJgJKZ0PaKOaEegaXjOT8R4+dwUH9jTn3+sJ+RnNpEmnmqOrpWpWJJWr5sW7/zFtLPVSzqrGn4J/GzcmCpevIAm4HMTq/McQSyZJpbMoIDjl5snAmw2tCPQNDwPHbkEwM/esDH/mKmE2SxRQeFUsUleeK7B00Mzi3SGTILLEBGEZ+drA683USqw2dCOQNPwfP+VYW68opMr1gfyj/WEco6gSaaLC6eKTUzhuUYvGC9eSmOyHF1D4ej8z/eodgR1QzsCTUPz5liE1y/NLEgLQYEjaJqIILmgUAwFwnMNPlRWLiKod2oonPv5tnnh9UvaEdQL7Qg0Dc1Dr1xCBH7m+o0LHu/OpYaaRW9oMpZa0DoKLFhO08gYEtTFIgI3yUyWZDpbt9c2hwZ3d7k5fnmWVKZ+r7WWccQRiMg9InJCRE6JyGeLPD8oItMicjj38XvVnqtZuyilePCVi9wysJ4NHS0LnutqM+6umykiWBcoHhE0uvDczFyK9hIRARgF3XphOvrrut0k01lOj0UrnKGxQ82OQETcwBeB9wG7gPtEZFeRQ59RSt2Y+/gDi+dq1iDHL8/y5liUn12UFgJDtG190NcUjiCTVUwXWUnZuUzLaZRSRJL2W1Rn46m8WmohwdxymnoWjMcjCdpbPFzdabzW0eHpur3WWsaJiOAW4JRS6rRSKgl8E7h3Gc7V1Jm5ZIYz4ZW7A/v+K8O4XcL7rttQ9PnuNl9TpIam51IoxZKuoRavC5/HVXcF0qfeGOPX9scYnbG+yCebVUQSpWsEQF3bX8ORJN0hPxuDQovXpTuH6oQTk8WbgPMFX18Abi1y3O0i8gowDPyWUuqohXMRkfuB+wH6+voYGhqyZWwkErF9biOwnPb/4HSSfzuV4i/2BQh4pfIJFbBiu1KKb78wx7XrXLx28Pmix3jSc7x5MbZs74fd9/5SxMhrX37rTYZSby14LuBWHHvzHENDI06YWJRHzqTIKPju48+yY7278gkFzKUVWQWjF88xNHRpwXNvjhmRwNPP/4ThTmvft1pOXZjDC8SiGfoDbp59/S2G2kbr8lr1pNGvO044gmJXiMVx6EvAVUqpiIi8H/g3YFuV5xoPKvUA8ADA3r171eDgoC1jh4aGsHtuI7Cc9j86cYRU9jy+TbsY3NVX8/ezYvvh81OMPfosv/0zuxnce0XRY/7t8sscOje5LO/H/3r8DVxz5/i/ftb6ax08OwEHnueOvXu4a3vPgud6X36KQEcbg4M3O2TpUp6fOwYnTrPpmmsZvGFpmq0cw1Nz8PiT3Lh7B4O3XLngucCZCT5/6Hl27t7Dndu6nTQ5zx8eGmLHhhBtbbPctrOLHxy5xF133bVkuK3RafTrjhOpoQtA4V/qZoy7/jxKqRmlVCT3+cOAV0S6qzlXs3KMzBhpl+dPjy/7az/0yjBet/DeXcXTQmB0DoVnk8si0fDlp97kuWF7uXBTZ2hxagiMzqF6p4bGc503ozPW02h5CepiXUPLUCMIR5L5DrFdG9uZnktxcWqubq+3VnHCEbwIbBORLSLiAz4GPFh4gIhskJwLF5Fbcq87Xs25mpXDXA7/3JvL6wiyWcVDRy5x1/YeOopcPE16Qn7mUpm6L0eJJdPEU1nCc/YczvxSGt+S5zpafXXvGjLrKHaG7+YF54rUCHz17RpKprNMz6XyjmB3fzugJ4zrQc2OQCmVBj4NPAocA76llDoqIp8SkU/lDvsI8FquRvAF4GPKoOi5tdqkcYaRmQQicOzSTF40bTk4+NYkl2fiS4bIFrNcQ2XmHfVE3KYjMAXngksdQWfAy0ydu4ZMR2AnIpgp5wj89V1XOZ6bKjZbhXduaMclesK4HjgiQ51L9zy86LEvF3z+l8BfVnuuZuVJZ7KEIwnesbWbA6fCvHB6nPctGuqqF99/ZZgWr4v3XFu+LlE4VLalO1g3e0wnOBlXpDJZvG5r90+TsRQ+t2uJeicYwnP1Higz9XrMCM8K89vJSqeG6hWRmQ64u80Pc9Dqc3N1T5ueMK4DerJYU5TxaBKl4D3X9hL0uZctPZTOZHn41Uu8e2df/o6zFMsVEZiOQAEjNlowDZ2hpeqdYEQEsWSGRLo+F1OlVP7O2s77VGxxvUmr141I/SICM5VlOnww6gQ6NeQ82hFoimJe8Po7W7lly3qeezO8LK/7/OlxxqPJBUqjpVi21FBBWmx4yoYjKDJVbGLKTNRrqGx6LkUqoxBg1I4jmFu6uN5ERAj66qdAauoM9RQ4gt397Vycmmt4ob5mQzsCTVHMfHJvewu3b+3izbGorbthqzz0yiWCPjf7dvZWPHZdwIdL6q83NF7w/YdtdKxMxVKsCxYvenfkHES9hOfM96YvKExEk5Z1gWbjabxuwe8pfqkI+uu3k8B0wGaNAGCXLhjXBe0INEUZyeWT+9r93LHV6BH/8TK0kT52bIT37OpbsMClFG6X0LUMm8omokncLiOtY6d1caJMRNBZ54hgLFcfuDJk/KlbdZqmBHWpvv2gz0OkTl1D4dkErV73ghThro05R9AkdYJ4KsNfPH6SRKaxt9BpR6ApymiuY6i7zc+1G9vpaPXy3Kn6OoJEOsNENMm23DrKaliOlZXj0SQ9bX7avHYjgmTRjiGov/CcWR+4st34U7f6Xs2WkKA2MaSo6+QIIgm6Qwvft642PxvaW5qmc+iZk2E+//gbvDza2Fvo9PJ6TVFGZ+N0BX35Dpnbrl7Pc6frWycwL4aLF7iUozvkr3tqaCKapKvNh5+kZUeglCq6nczEFJ6rV+eQmWc3IwKrdYKZ3JrKUgR8bqJ1Sg2FI0m6gv4lj+/qb56C8dmcVteZ6cZ2BDoi0BRldCZBb2he+vn2q7s4PzHH+YlY3V6z3OBVKZYrIlgf9NHVIpaLxTPxNJmsKl0szkcE9Sl+hiNGWmtTm+kIrNk/G0/T3lr6frHNX78tZeFIYkHHkMnu/nZOjUWaYl/16ZwjODvd2HsUtCPQFGVkNk5v+/wf4R3XGHWCespNTEZLSzGUoifkJxypr8zERDRBV9BHV6tYjgimKji3kN+DSP1qBOFIgvVBHx1+QcT6UNlsPEXIXyYi8Hvq1j4ajiTpCS1933ZtbCeTVbwxMluX13USMyJ4ayZLJtu4dQLtCDRFGZlJ0FcQEWzrbaO7zcfzdZwnyEcEJfLpxehu85HMZJmZq5/ezUQkyfqgn/UtLmYT6fy0bVXn5qeKi19MXS4x9Ibq2DXU3ebH4xLWB3yWU0OVagRtfnddBsoyWcVEtFRE0AE0x4Tx2fEorV438QycHoustDkl0Y5As4R0Jst4JLEgIhARbt/azXNvhut2920rNZRfYl+f1tZ4Tsuoq82ICMBawXgqLzhX+v/U2eqtX9dQJEl3rv2yJ+RnzGJqaGauUo2gPhHBZCxJVkFXkZuCzetaCfk9DV8nmEtmuDQd5+6ccu+RC427VEc7As0SxqPGH2Fv+8L1kLdf3cXITCKf93Sa+WKxtdQQ2BuWqgbzjt6sEYA1R5CPCMo4go6Ar67FYvOuure9xdL7lMkqoslM2RpB0O8hlsyQdTjtkZeXCC2NCFwu4dr+9obfVnZ23Pg7efe1vfjdcOTC1MoaVAbtCDRLMPPIfYv+CO/Y2gVQt/TQZDRJwOeuaobApCevN1SfYmsxR3DRQsG4miins9XLdB2KxUqpXGrIeO3ekLXCeiQvL1HaMbfl9IZiDhduw0XkJQrZtbGd45dnGzrvbtYHtva0MdDu4sjFxnVc2hFolmBOEC+OCK7qCtDf0VI3R1Bu8KoU9ZaZMC9IXUEf7X7B67ZWMJ6KpXC7pGyevTNQH+G5aDJDIp2djwhyjqDau/dyyqMmAVOK2uH0UCVHsLu/nVgyk7/rbkTO5Gwb6A4y0GGs2UxlGrN7SDsCzRLM9EFf+8I/QrNO8PzpccdTAWBcNK2khcDQ6vG6pW6zBIURgUuEjR2t1lJDsSSdrV5crtIbtepVLDZnCMyLaU/ITzqr8lFKJUxHUExnyKQtN/XrtN5QOK88WvzGoBmkJs6Go/SE/LT5PWzpcJNIZxu200k7As0SRmbi+anixdy+tYuJaJITdfiFLifOVgoRobuOswSmIzAHm/o7WyxGBKWnik06W73MxFOOO9f8XXXIjAiMCK/aOsH8drJyEUFOitrhobJwJIHXLXlRvsVs6w3hdUtDdw6dCUfZ0mXIo2/JTXa/2qAFY+0INEsYnU0smCou5PY61gkmo5UvmsXosZj7tsJ4NInHJfmCaX9nq6WhsolosuJcREfAh1LzF16nmE+v5GoE7dYK67NV1Qhyy2kcHioLzyboCvpLahz5PC629YYaWnPoTDjGQHcAgN6A0N7i4ZXV7AhE5B4ROSEip0Tks0We/7iIHMl9PCciewqeOysir4rIYRE56IQ9mtoYnYnTE2op+tymzlYGugJ12U9QToqhHD1t9ZOZMGYIfPkL0qbOVi7PxElXmes10l2VIwLA8d3FY7n0Sk9BjQCMn281mBLUZWsEddpSFo4kFqiOFmN3fzuvD08vy85qq8zGU4QjCQZyC5NEhBs2d/LqxamVNawENTsCEXEDXwTeB+wC7hORXYsOOwPcpZS6AfhD4IFFz+9TSt2olNpbqz2a2hmZjS+pDxRy+9ZuXjg9XvXFsBoyWcVMvPJFsxj1TA2Z8hIm/Z2tZLKKkSpfbzKWZH0lR1An4TmzRmDabz01lKsRlEjPwHzXkNNDZePRZMlCscmu/nbCkWTdJUbs8Na4IcVydcHmvOs3d3D80mxDSmM4ERHcApxSSp1WSiWBbwL3Fh6glHpOKTWZ+/LHwGYHXldTJ0YXTRUv5vatXcwm0o7mZ6fnUigF6+1EBCG/MftQhwL2RHThnWl/ZytQ3SyBUorJaIrOElPFJnlH4HDnUDiSYF3AiyeX4mv1uQn5PVVfOGfLbCczMbuGHI8IZotPFRfSyBPG5qzNQIEj2LO5g3RWcfxy4xWMnVAf3QScL/j6AnBrmeM/Cfyw4GsF/EhEFPDXSqnF0QIAInI/cD9AX18fQ0NDtoyNRCK2z20E6m1/VinGZhPEJi8zNDRR9BiVMC64//jYT5i8uvo7+HK2D0eM6OLSW28ylHrLks1Tl1NksoqHHh+i3Ve6O8cOF8MxBtpdDA0NEYlEmHnjVQAef+4lomfL//nE04pkJsvk5fMMDY2UPM78vz9/8BXUsHOCwMfPxml1ZfO2Dw0NEfRkOHr6PENDYxXPP3oygc8Fzz7zdMljoinjd+HI6yfYGDvtiN1KKUZn4sxNjuR/X4r97syljdd+6NnDyGXrkWQ92X/KSMu9dfQQl48LkUiExNwxAP51/4tMXWn9hqeeOPFbV+wvr+itmYjsw3AEdxY8/A6l1LCI9AKPichxpdSS37ycg3gAYO/evWpwcNCWsUNDQ9g9txGot/2jM3HUo09wy/U7GLztqpLHffH1p7hMK4ODt1T9vcvZfvDsBBx4njv27uGu7T2WbI4eucTXj73E9htuZueGdkvnViI29CjXXr2ZwcHdDA0Nse/2O/mdA4/S0T/A4OA1Zc89PxGDx/ez97prGXz7FSWPG5tN8DsHHqd/yzUM3j7gmO3/+9hzXBV0MTh4W/69v+rE82SVYnDwjornPzJ+hI7x0bK/b+lMFp74IRs2DzA4uM0Ru6fnUqQf/RE3XXsNg++6Gij9u3PVy/uZa2lncPBmR17bKR4cOUx/xzjvffc+wLD/Z+66iz8+9Dhzrb0MDu6p8B2WFydSQxeAwt/yzcDw4oNE5Abgb4B7lVL5SqNSajj37yjwXYxUk2aFGDFXVBYZ7S/k9qu7ePHMhOXVh6WYjFlXHjWp11BZMp1lNp5eUCMI+j10BrxVpYaqlczI7y12uEYwHkkskWiwIjMxG0+XbR0F8Lhd+D0uYg52DY3n214r3+U36jL70+HogrQQGAXj6zd1NGQLqROO4EVgm4hsEREf8DHgwcIDRORK4DvALyql3ih4PCgiIfNz4KeB1xywSWOT0fyKytI1AjAKxnOpDK84pJ8yWYUmTylMR+B051DhMFkh/R3VtZCag1uLz1+Mz+Mi6HPXoUaQXDKQ1RvyMzqTqKrTptJSGpOg39kF9vPDZOVvRsDoHDo7HssXthuFs+NLHQHA9Zs7OTk666jjdIKaHYFSKg18GngUOAZ8Syl1VEQ+JSKfyh32e0AX8FeL2kT7gAMi8grwE+AHSqlHarVpLfPEsRH+6zdett1SV21EcNvV6xFxbp7AjgS1iXmxczoiMNc8LlbANGYJKkcE5v+pmk6ozoDP0a6heCpDJJFecjHtDfmZyz1XiZkKEtQmQb+bmINdQ5XkJQoxJ4wbqQA7FUsyFUvlh8kK2bO5g6xqvAK3I3MESqmHlVLblVJblVJ/lHvsy0qpL+c+/49KqXW5FtF8m2iu02hP7mO3ea7GPn/xxEm+/8pw/m7WKqbOUE8FR9AZ8LG7v53n3nRmfeVkLIXXLQR91QvOmbT5PbR4XY47gvxU8aIL0qbOlqqW2M9HOZXvqjtavUw7OEdgvhc9ix1Be/VptNl4qmzrqEnQ52xEYKaGKs0RQEHnUJWCbnPJTN3vxs/kOoa2FI0IDHtfOT9VVxusoieLVxEnLs/mNc/tinGNzhpqlcWmihdz+9VdvPTWlCN90ZNRQ16i1CRpOUQkv6nMSUqlhjata2U2XnlBzUQshQglZRIKcVpvKFwiz25llqCaGgGYUtTOXVzHIklEqDh/AUaE0xX0VTVh/OqFad75p0/ym996xQkzS3KmSOuoSW+ohY0dLbzaYEqk2hGsIv7l0HwX75mwvd3C5aaKF3PH1m6SmSyH3pqsfHAF7OgMFVKPoTJTE79YagjgUoU6wVQsSXvLfB9/OToDzi6nCUcWaiSZWNnfMGupRuBsamh9wFfV+yYi7Opvr5hqGToxys8/8DzhSNKR39dynA1HcQlcuT5Q9PlGLBhrR7BKSGWyfPfli7x7Zy9ul+S10K0yOpsoO1VcyNu3rMftEkfqBHaURwupxxL7iaix+H3xHX21Q2WTsVTFQrGJ01LUiwXnTKqVmUims8RTWUL+KiICn9vRgbLwbGV5iUJ29bdzciRSsoPt2wfP88mvHWSgK8gn79zC6Gwiv0u6HpwZj7F5XQCfp/jl9YbNHZwOR+u2lc4O2hGsEp46MUY4kuS+W65k87rWvBa6VUZm4hULxSZtfg97Nnc4UieoNSIwUkNOF4sNwbjFEtKbco6gUp1gMpqs2rl1tPqYjqUc080x5SUWRzMdrV58nsr1lGrkJUyCfo+j+wiqkZcoZHd/B8lMllOjC3cCK6X430+c5Lf/5Qh3bO3in//zbdx5TTcAb4zUb3/w2SKto4XcsLkTqL6usRxoR7BK+Pah83S3+bhrRw8DXUFbEUEma2y0qtQ6WsjtW7t45cJ0zcXCySrkmsvR3eZnIpZ0dPHHRDRR9I6+p81f1YIaK86tM+Almcky55AOzXg0SajFs2Tbm4jQ0+avmBqqRl7CJOhzO9w+WlleopBdG3O7CQrqBOlMlt/57mv82WNv8OGbNvG3n3g7oRYv2zeEAOq2F0AplZOfLp4WAiM1BDSUEql2BKuA8UiCJ46N8qGbNuF1u9jSbTgCq3eX45FE0V3F5bhjazeZrOLFs8XlKKpBKcWUTeVRk56QH6Ww3S1VjIlFgnMmLpewoaPyXgKzAF4NZvrJqYLxWCSxpGPIpLfdn58XKcX8drIqI4JkxtFoxooj2NIdpNXrzu8wnktm+NTXD/GNn5zjvwxu5c9+bk8+TdPf0UKb31M3RxCOJIkk0mUjgnVBH1euDzSUEql2BKuA7x0eJp1VfORmY8B7oCtANJlhzGKqxLxLrDY1BHDzVevwuISXaijAzSbSpLOq5tQQODtLMB5NLim2mlQzVGZFVrvTYUdQ7mJaze7iapbSmAT9HtJZRcKBKfO5ZIZoMmOpRuB2CTs3hnh9eIbxSIL7vvJjnjw+yh9+8Dr+73t2LuhEExG29bXVzRGY3XrFWkcLuX5zB6+c1xGBxkG+fegCN2zuYEcu7DXvRs5a7BwyZwispIZavG56Qn4uT1e/rGUx+X77GlNDgGXnV45SEQEYdYJyNYJ4KsNcKlP1/6kj5zCcKiCGI4mSEg29ocoyE7NWIoLc7IcTQ2VmnadUNFOK3bnOoY98+XmOXZrhS79wM79YQitre2+objWCcjMEhezZ3MHFqbn8zMRKox1Bk3N0eJpjl2b46M3zyt5b8o7AWp2g2qnixfSEKuecy1GLzpBJr8MRQSqTZapM109/bkFNpoT0dX5SutoaQatxnFNDZYa8RPGfY0/Iz1QsRSJd+sI9Y6FG4ORymlLzD5XYtbGDSCLNZCzJP/2nW3nv7g0lj92+IcRENFmXZUZnwlE8Lsk3FJTi+k2dAA0zT6AdQZPz7YMX8LldfGBPf/6xTZ2teFxiuXPIzBtXmipeTG/NjqB6KYZSmBc9p/64TZtKLU83F9SUyrVPRq05NyeX0yTTWabnUiXTWtU4TXM7WbnF9SZOrqu0ojNUyHt39/HRmzfzL5+6g5uvWl/22O19bQC8UQdZirPhKFeuD1ScgbhuUzsi5AdAVxrtCJqYZDrL9w5f5O7dfQsuoh63iyvXB2xFBKV2FZejJ9TCWIXiYzmmYtVLMZSi1eemzcLSlUrMTxWXqBF0GumzUgVjq9pJTi6nMTWSSqaGqthdbNYI2qqJCPIL7Gt3BPPyEtYcQVebn899dA/X9LZVPHZHX/06h86EoxXTQmCk3K7uDmpHoKmdJ4+PMBlLLUgLmQx0B/P5ymoZm41b6hgy6c1tCLO7unIid/dc7fBVKZxcYj8RKa8cOj9LUCIisJgaavW68brFkYggPFv+rjovMzFT3hG0+T24XZUlP/IRgQPTxWZEt3j+wUl6Qn46Wr2ccLhOkM2qkqqjxbhhcydHHFLvrRXtCJqYbx+8QF+7n3duW7rIZaAryFvjMUstfSMz1U8VF2K2bo7bbN2ciiVxSXVpiLJ2OLjEfjwvOFf8grTRdASTpSKCXGqowppKExExhsocqBGEzYigTNcQUDaKM+QlqttbFXS0RlB8/sFJRIQdfSFOOhwRjMzGiaeyFhxBB6OziXyTxkqiHUGTMjobZ+iNMT78ts1F79q2dAeYS2XyBeBqv6fVQjEUyhbYuwhPxpJ0tC6d4LVKd8hXh9RQcUfQ5vfQ0Vp6QY3ZCWUWgavBKb2hcAnlUZOuNj8uqVAjiKeqdsxBc2+xA11DYxaHyeyyra+NEyOzjs0+wHzH0NUWHAE0hhKpdgRNyr+9fJFMVvGRImkhmG8hrTY9lMkau4qttI6amOmkSkNKpZiMpmpqHTVxUm9oPGooYJZL7ZTbSzAZS9Lm95TUmylGp0MKpPmCa4kagdslrA+WL/DPVrmLAIx9BOBcjaBUgd5JdmwIMRtPW7pRqoTZrl1tRLBrYwdulzRE55B2BE2IUopvH7zA267sZGtP8eLYQG4pRrVy1OPR3FRxDRGB3YtwrTpDJj0hPzPxdNm2yGqZiCbobPWWzZGX20swGU1WnRYy6Qw45QgSBHxuAr7SF/JKnV7WHIFxnBMyE+XaXp1kW6/zBeMz4Qh+j4uNVd5MtfrcbOtta4iCsSOOQETuEZETInJKRD5b5HkRkS/knj8iIm+r9lzNUo5cmObkaISP7i29EL2/sxWf21V155CZ1rFTLDb/cO22kFqZwK3GDif2EpQbJjMpHxGkLDs3o0bgjCOodDGtJDNR7ZpKAL/HhdsljuwksKozZJd8C6mjjiDGVV0BSynOGzZ3cOTClKMpKjtU5/LLICJu4IvA3RiL7F8UkQeVUq8XHPY+YFvu41bgS8CtVZ67ajA1dcKRBGORBGOzCcKRJGOzCUItHn75joH83VU5vn3oPC1eFz9zw8aSx7hdwpVdgapTQ2bByk5E4PO4WBfw2k4NTcWS7M6tHKyFQpmJSgM9lRiPlJaXMOnvbGUmni6q228nyjGW0zhQLI5UlnHuDfnLLn2fjadpb63u8iAiBHzumruGzCE+K/ISdulq89Pd5nPUEZwdj7K1p7q0kMkNmzv51sELXJic44oS+wuWg5odAXALcEopdRpARL4J3AsUXszvBf5eGW7vxyLSKSIbgYEqznWMLw29ybeem+PPXztQ9jgrztklxh+CS4yLr/m5SwSXCCLGkNDYbILxaIJUZuk397qFVEbxTy+c4w/u3c27r+0r+XrJjOLBw8Pcs3tDxWLeQFew6tSQeTdvp0YAOdkCm/nWau6+qyG/xN6BOsF4NMk1JdJuJvkFNdPxoo6g2qKhSWfASzSZIZXJWp7lKCQ8m+TKMuqXYPy8wpEEmaxakv5SSlW9lMakze+puUZgFuiXIyIA2N4XcqyFNJNVnBuP8e5rey2dZxaMj1yYbnpHsAk4X/D1BYy7/krHbKryXABE5H7gfoC+vj6GhoYsGzp8LkXAnYFE+R++5D4qoTCcRgbjjyer5h9TYHytIOgVtoWEm7s8dPiFdr/Q4RPjc58Q9MKpqSx/dzTOJ792kL19bj5+rY91LUsvBs+fizITF7Z7Jyq+B565BGfG0jy5fz+uCisgXzhl/BG+/tKPecNG944nPcebF6NlbYpEIkueT2QMsbKJy+cZGhqx/LqFjM8ZcwwHDh3BM1pbqmlkMsoV/vgCexfbPzpp3AE/8vQLDPcs/FMKz8wRnUxZ+j0dvWCkhX74+FO0++13UA1PRtnonStr+9TlFFkF339sP53+hb9nyYwilVGMDZ9jaOhydS+aTnDmwiWGhuyLD741Y7yfI2+dZCh+ZsFzxX53aiWQSvDScHV/H5UYi2VJZrIkwxeK/h6Xsj+VVXgEfvD8qwQnTtRkQy044QiKvYOLb3tLHVPNucaDSj0APACwd+9eNTg4aMFEg0FgaGgIO+fWm33AJz6Q5SvPnOYLT5zk/30+xW+/dwe/cNtVC+7Y/ufBR9jU6eNTH9pXMRd5sfUtHjn7Gjtuuq1iquRHk6/Sdeky7/mpfbbsf3DkMC+cmSj73hZ774en5uCxJ7n5up0M3nKlrdc2SaQz/OZTj7C+f4DBwW22v08mq4g8+jDXbxtgcHBH/vHF9u+cjvPfX3iCriu3MXjrvMBZMp1l7pEfcv32qy3ZMX34Il8/dphdN729qgnZcrbfsL287fHXLvH1Yy+x7fqb8wvgTUZn4vDYE+zZtYPBEsJti+l97QDBgI/BwVts2Q3GOkmee5G7bn0bewcWykTU4+92uPUcj597lW17bq35bvzpN8bg6Z/w3ne8jduu7lryfDn7d71+gEnxMDh4W8XXGZ2N09Pmt7XbuxxOFIsvAIVVy83AcJXHVHPumsHncfGr+67hR7/xLm66spP/9uBRPvxXz+Z11i9PxzkazvDv37apqoLUlq7qxedGZ+xNFZv0tButm1aLXmY6wImuIb/HTWfAW/NQ2VQsiVKVJ517Qn48rqULaqbm7ElmmDIhtQyVTUQN2xevqFxMud3FMxYkqE0CvtpTQ+M2dYbsYhaMT47WXicwU7BW04FgLKp57eI02RIChtms4tlTYf7LPx7ijj9+kh+ftr/7oxROOIIXgW0iskVEfMDHgAcXHfMg8Eu57qHbgGml1KUqz11zXNUV5O9/5Rb+4mM3cnFqjn/3l8/yRz94na//+C0U8O9LzA4sxsoswehswlah2KQ31EIyk7Xc9TLlgPJoIU4ssc8Pk1W4ILnzC2oWFsnzgnMW6x5O7CTIq3dW6hrKyUyMFanr5NdUWqgRBP2emgfKSu1ZrhfbcppDJy7XXic4E44S9LktCzYC7NncyWwivaSeNxlN8jfPnObdf/4UH/+bF3j+zXF+5c4tFes/dqg5NaSUSovIp4FHATfwVaXUURH5VO75LwMPA+8HTgEx4D+UO7dWm1YDIsK9N25icHsvf/LIcb7yjJEz3bHOxVVd1d11bGhvwe+proV0ZCbOztw+AzsUduxYURG1Ks5W0Q4HHEFeXqIKm/qL7CWwqjNk4sSWsmodwXxEsLTTy4oEtUnQX/sC+3AkQYvXld9vUG86Wr1saG9xRGriTDjKVV1BWymb6wsKxlu6g7x0bop/fOEtHjpyiWQ6y96r1vFr797GPddtqJv0hhM1ApRSD2Nc7Asf+3LB5wr41WrP1czTEfDyxx++nn//tk18/vE3uL2z+rsXl0uq6hwydhUn83eJdugtSDWYd1rVMJWXoHYmIugJ+WsW8qokL1HIps7WJWs67f6fOh1YTjPvCMrb3uJ1097iKZoasrK43sRYV1mrIzBadp3Of5dj+4YQJxxwBGfDUXZv6qh8YBG29bbR4nXx9R+/xV8/fZpjl2Zo83v42Nuv4P+49Up2bqi9tboSjjgCTf3ZO7Cef/yPt1nunBjoDnBqtLzzGI8abYR2BOdMesvcYZZjIq/b70xE4ERqyFpE0MLl6fiCNky7aqqhFq/RblyLI5g1xfIq/yx724u3/FpZXG/S5vfUPFlsbFVbnrSQyfbeNv7h9HjRNtpqSWWynJ+c42dv6K98cBE8bhc3bO7kJ2cm2N3fzh9/+Hr+3Z7+qmaKnEI7glXOQHeQ/cfHyv6i1zJVbGKea/UiPBlLEvJ7auqbL6Qn5CeazBBLpstKLJTDlKCuJl3V39lKOqfTtKHDeA/spobcLqG9xct0DUNl4UgCn9tVVaHXkJlY6ritrKk0CfjcxFPZmi6o4UiSTZ32fwftsH1DiEQ6y7mJWFV7BIpxYXKOTFZVrTFUjM///I1MRo3ByuWMiEy01tAqZ0tXkGQmW1IKAebv4mspFgd9blq9bstDZVOxJJ0WNXnKMT9UVkvnTYL2luqcU39+L8H8+zsZTdLqddvK53YGvDVFBGM50bZqLiY9IX/RHc8zc2lcgqVcvRNbypZLXqKQ7Q4sqTkTNiLuLd32i7ibOlu5blPHijgB0I5g1VNN55CpwGh3qhiM4rahX2PNEUzEUqx3KC0E87nxsYh9jfdwtHrhM3M+o9DR1qKd1FGjAul4JFl1eqU35Gd0ZmnLrzlVbOWiZEZfMZsyE9msYiKaXBZ5iUK29da+tvJMTnV0S7e92Y9GQDuCVU5+kX2ZgrF5F1/r3VipVEM5pmLJmnYVL6anRiVUMFJD1eb3N3YsXVlZy/+po7W2iMDKXXVvqIVEOpvvEjKxojxqYkpR260TTM2lyGTVskcEQb+HzetaeaNCHa0cZ8NR2ls8jrVArwTaEaxyekN+Aj53+YhgNk5X0GdJO7/4a7XYqhE4+QeUdwQ1KJBa0T4KtXhpb/EscAQTMfvaSZ0BX35xvB3CFvT8zd3FizeVzcTTluoDULCcxqYjqLbttR5s7wvVGBEYe4pXKq3jBNoRrHJEhKu6gmVnCUZnErYGYRbTU0HjvhiT0ZSjEcH6gA+psH2rEuMWUxTGLMH8xXQqlrLdDttZgwJpNquM1FCVF9P8LMGiuo6xncxqRFBbjcAUClzu1BAYjuB0OELK5s7tahfWNzLaEawBtnQHODseK/n86Gy8pvqASU/Iz2w8TTxVXZ44mc4SSaQdUR418bhddAV9tmUmslnFpMU7+k2LhsomovYX7ZjrKkvJDZRjei5F2kJ6Jb/EfpHTnLUTEeS3lNmrEYRzLbul1mvWk+19baQyqurdHYXEUxmGp+dq6hhqBLQjWAMMdAU5PxEjXeKOZ2QmXtMMgYnV3cV2NXkqUcsswUzcyFWvr7CLoJDCBTWZrGImbn/1Zkerl6yCWRspFtP5VXtXbaaGFtd1ZmuICOwOlZkRwUqlhgDesCFJfX4ihlLoiEDT+Ax0B0lnFRcml7aQOjFVbJKfJaiyY8fsjnEyNQS5tkibjsDKMJlJf2cr03MpIok003MplLLv3EyZiWkbnUNmK2i1d9Uhvwe/x7XkvZqZS1kvFvtqW1cZjiTwuCT//19OrultwyXYmjA+nYsiBqqUfWlUtCNYA5h3K2eKdA45MVVsYl6Aqo0InFQeXWyH3dSQFXkJk/7cENSlqTlb5xcyr0Bq3RHML62v7mdZrOVXKUUkkbYkLwHzqSG77aPhSIL1QZ+lNY9O0eJ1c1VX0JbmkJlO0qkhTcMzUEaO2rxo9zgSEVjbXTyVF5xz9i7QjAjs7IE1pZCt1gjAGCqb1xmyXyOA+bSZFcZtdN4s3iwXTWbIKmvyEjA/R2A3IrBS5K4H23rbbA2VnR2P0hX0rUgk4yTaEawButt8tPk9xR1BLj/sRESwPuDD7ZKqZwkmY87qDJl0t/lJpLO28uzmHb3VriGA4al4QZRjv2sI7CmQhiMJ3C7Jf49qWDz7YUdeAgx5jFav236NYAV0hgrZsSHE2fEYibS1iOZMONr00QBoR7AmEBEGugOcKdI5NOrAVLGJyyV0t/mqzs/XLTVUw+7iiahxjpWIoDfkx51bUDNVo3PryEcENhzBbJIui+mV3kUtvzNz5lIa644s6HcTsZ0aSlY9/1APtvWFyGQVp8esdQ6dCUebvj4A2hGsGQZKzBKMODRVbNIbarGUGmrxumh1WH++lunicCSZK6JWb5PH7WJDewvDU3M171eYLxZbTw3Z0erpbW9Z0PI7HxFYF+yzK0WtlFoRnaFCdtjQHIol04zMJLi6RzsCTZOwpTvIhckYyfTCFtJRh6aKTUz9mmowNHmcvws0LyjFBNUqMRFNst7Gnak5SzARS+Jz21+u4vcY4n12U0NWB7IWF/jtSFCbBG2uq4wk0iTS2RWNCLZ0B/G4xJIjOJvTGFrzEYGIrBeRx0TkZO7fdUWOuUJE9ovIMRE5KiK/VvDc74vIRRE5nPt4fy32aEoz0BUkq+DcxML00IhDU8UmVoTnnNYZMqktNWRPHqK/s4Xh6TmmosZUcS1yA+ZQmVXCkaTlgaweU2Yi1/I7Y7NGAOaWMuupofAy7youhs/jYkt30NLaSlO/a6AG1dFGodbbwM8CTyiltgFP5L5eTBr4TaXUtcBtwK+KyK6C5z+vlLox96E3ldUJs6C1OD3k1FSxSU+bn4lcS2oljAlc57stOlu9eFxiKyIYjyYtzRCY9He2cnk6zng0UXOUY0d4TillSFBbdOqLhwDzi+tb7aWG7EhMzA/CrZwjAGOwzMoi+zOrZIYAancE9wJfy33+NeCDiw9QSl1SSr2U+3wWOAZsqvF1NRYppUI6OlPb0vrF9LS3kFXzrYzlmIrZn8Ath8sldFkoWhcyEU3YjAhaSWUUb4xEam6H7Qx4LQ+UzSbSJG2kVxbLTNhZXG8S9NnbUjZe5XrNerO9L8S5iRhzyeqimjPhKH3t/mXdJFYvav0f9CmlLoFxwReR3nIHi8gAcBPwQsHDnxaRXwIOYkQOkyXOvR+4H6Cvr8/yykaTSCRi+9xGwK79SikCHnj2yEmuyZwDIKsUo7Nx4pMjjr0noyPGheDh/c8y0LEwT77Y9tHpKAOtibr8PFpIceKtSwwNFf11KopSivBsgsh48fej3Hs/Pmb8v89NxOj11vZ/SkXjXI5mLX2Py1Gj9jN2/jRDQ+eXPF/K9qxSuARefO0Nrkqe5eiJJG6B5w88bTm9NTORYHImY/n//tw5w/mcPHKIsTeK35sux99tMpxGKfjnHw4t+d1dzHRC8dhrc2xuk6rsavTrTkVHICKPAxuKPPW7Vl5IRNqAfwV+XSk1k3v4S8AfAir3758Bv1LsfKXUA8ADAHv37lWDg4NWXj7P0NAQds9tBGqxf9vRAyT9HgYHbwOMrprso49zy/XbGbx9wBH7Os5N8r9ffo4rd1zH4M6+Bc8V2p7JKqKPPszubQMMDu5w5LULufrMTwhHkgwO3ln1OTPxFJlHf8RN117D4LuuXvJ8ufd+4+VZPn/oaQC2D2xicPB6W3YDPDJ+hPPHRy39nF88OwHPPM87334j79res+T5crb3PP84ret6GBzcw2OTr9Ixepl9+/ZZtnto5iiHwxcs/36+/Ngb8PpJfvbuQTwltsItx9/tFWMRvnj4Kdo272Dw5s0lj0tlsnz8b14gkY3zuY+/g139lZfLN/p1p6IjUEq9p9RzIjIiIhtz0cBGYLTEcV4MJ/CPSqnvFHzvkYJjvgI8ZMV4jTUGuoMcPDt/hzwyYxQInZgqNqm2dXMmr8lTn3RAd5ufY5esTYpO2JgqNukv2LVba93DTo2gFtG2wpbf2XjasuCcSVuuRqCUshRNhCMJ1gW8JZ3AcnHV+gA+t6ui1MSf/PA4Pzkzwf/6+RurcgLNQK3v/IPAJ3KffwL43uIDxPiN+FvgmFLqzxc9t7Hgyw8Br9Voj6YMA11Bhqfn8j3j5sXaialik1Ia94uZrJO8RKEd49GEJTnncRtTxSahFm++5bLmYnHASzKdrVrOGwoWu4Ssv3ahSJ+5ptIOAb+brIJ4ypqu/0rLS5h43C629raVFZ/73uGL/O2BM/yHdwzwwZtWT6mzVkfwJ8DdInISuDv3NSLSLyJmB9A7gF8EfqpIm+ifisirInIE2Af8Ro32aMqwpTuIKmghNSOCXge7hvweN50Bb8UW0sk6KY+a9IT8pDLKUhtmXl7CggR1IabmUK2OoLPVON/KLMFYJIkItvY/F04X21lTaWJ3gf1KD5MVsr2vjZMl5KiPXZrhM/96hFsG1vM77792mS2rLzUVi5VS48C7izw+DLw/9/kBoGicqJT6xVpeX2ONwkX22/tC+T9+p5eBVLO7eLJO8hImhUNl1XYmmd0rdgbKwOgcOn551pGuITCE5zZ0VOekjfSKz1Z6pTcXPaUzWWbiKa62uYQ9ULCu0sqFPRxJcN2mDluv6TTb+0J87/DwkshoOpbiU18/REerl7/8+E14VziN5TSr63+jKcuWRSqkIzNx1js4VWxSzT4AMzVk5w62WhvA2lCZnV0EhZh1gtojAuvCc+HZ6ncVL6anvQWljP9/bRGBvS1ljZIagvklNScLltlns4pf/+eXGZ6a468+frMjuzsaDe0I1hAdAS/rAt78LMGIwzMEJtXoDeWX0tSxRgDWZCYmokkCPjctXnvyEP0OpYba7TiCGtIrhUNldtZUmuQjAgupoXgqw2wi7eh0ey3kNYcKltn/xRMn2X9ijN/7wG5uvmqJeMKqoPknITSWGOgO5icixxyeKjYxc87lukcmYkk8LiFUp2GcfGrIQkRgV17C5Kd3beDMWJTN61ptfw+YTw3NWKhvjEeT7Nncaev1TEdwaXqOSMJ+RJBfYG9hqCzcIMNkJpvXtdLqdefXVj5xbIS/eOIkH7l5M79w65UrbF390BHBGmNLVzAvllWviKAn5CeZzublCoph6gzVoslTjvYWDz6PK18Qrwa78hIm1/S28bmP7qm5DdIsoFtZTmOkhuz9LM27cfMGwep2MpN8sdhCasjUGbJboHcal0vY1mcsqTkTjvLr/3yY6za1898/eF3dflcbAe0I1hgD3UEuz8SJJtKMRRJ1iQjmZwlKX4Qno6m66AyZiAjXb+rghTMTVZ9jV17CaYI+Nx6XVJ0amktmiCYztlpHYf7n9eaYcRdsNyII5BRXraSG8vISDZIaAtjWG+LYpRk+9Q+H8LiEL//CzbbThc2CdgRrDLNz6OVzU2SyKr9e0kny+jVlZgkmY8m6dQyZ/NTOXo5cmK56Y9pEJMn6BrgzFRE6A9UPlYVtrKgsxGz5fTO3lKWWgTJo7tQQwI4NbYxHk5wcneUL993E5nXNry5aCe0I1hhm59ALZ8YB6tIBUc3u4slYMp8LrxeDOwyphaETYxWPVUoZqaEGuSB1tFYvPDfmwMW0N+TnVK5Txo7gHBgDZWDVEay8BPVirus3Wll/6707eOe2pXIdqxHtCNYYpnb6C6eNlEk9IoJqZCYmY6m6p2F2bWxnQ3sL+48XVT5ZQCyZIZHO1lQjcBJDZqK6GkEt8hImvaGW/PCd3a4hv8eN1y1Eq1TvBON3pM3vaajUy+1bu/jB/3Un/+ddW1falGVDO4I1RqjFS3ebj8PnpwBndhUveQ2/hxavq2RKRilVt6U0hYgI+3b28MzJMKlMedkDc6q4EWoEYBSMq52KduKuurBpwG6NAIwWUisRwXh0ZXcVF0NE2N3fsaqLw4vRjmANMtAVJJm7MDo9VQzGH1K5WYJIIk0qo+paLDbZt6OXSCJtqHOWYX45SmNclDpbvVUXi52wvafdGUfQ5vdY6xqqodtJ4xzaEaxBzIJxPaaKTcrtLjYvcPVYSrOYd1zTjc/tqpgemo8IGuOi1GFhOc14JEF7iwe/x356pfCGwG5qCIzOIavF4kZxvmsZ7QjWIOa2snrMEJj0hPwlp3rzyqN1Tg2BMeR069XrebKCI6hVXsJpOlt9zCbSFVNaYKSGam2/NIUHW7yumm4OrK6rNFJDjeF81zLaEaxBzB2rTqqOLsaICIrXCCbygnP1Tw2BkR56cyzKufFYyWMar0ZgvDem0yzHmAPqneZNQS3RAJgL7KtzBOlMlsmYdgSNgHYEaxCzc6ivjhFBb3sLM/F0UU39qTpLUC/mp3YaG1SfPD5S8piJaBK/x5UfilppbtjcgQj80Q+OoVT5nQrhSKLmWo/pCOzOEJgEfR5iVXYNTUSTKNVYw2RrFe0I1iADXUFcAhs7a9PEKUe5FtK88ugy3X0PdAe5ujvIk2XmCcYjhrxEo3SK3HTlOn7z7u187/Awf3vgTNljw7O159nN6LD2iKD6Bfb5+YcGicLWMtoRrEGCfg9f/+St/PIdA3V7jfymsqKOIIWI0Su/XOzb2cuPT48TK5G/nogmbO8hqBe/uu8a7tm9gf/x8DGePRUuekwinWEmbk3/vxhtfg8Bn7umjiEwUkPVRgT5tlcdEaw4NTkCEVkvIo+JyMncv0U1WkXkbG4T2WEROWj1fI3z3HFNd13vyHvL6A1NRpO0t3hxu5bv7nvfjl6S6SzPnRov+vxENNkwwmcmIsL//Lk9bO1p49P/9BLnJ5bWOMYdnMy9cn2g5klzKxHBkdwsy9Yee4twNM5Ra0TwWeAJpdQ24Inc16XYp5S6USm11+b5miYirzdUIjW03EXZW7asJ+hz8+SJ4t1DtSqP1os2v4cHfmkv6aziP//DIeYW3W07qdXzlV/ay++8f2dN3yPo85BMZ6vqdjpwKszu/vaGKdCvZWp1BPcCX8t9/jXgg8t8vqZBWR/04XZJ0VmCqViq7jpDi/F5XNy5rZv9x0eLFl9r3UVQT7Z0B/nCfTdx7LKxM7fQ/nEH0ytXrA/QVWNkYe4kiFUYKosm0rx0bpI7t3XX9HoaZ6h1K0ifUuoSgFLqkoj0ljhOAT8SEQX8tVLqAYvnIyL3A/cD9PX1MTQ0ZMvgSCRi+9xGoJnsD3nhyMmzDPkvAfO2nxuZY12LLPv/o19SPDqd5OsP7eeK0Pw9UCKjiCUzTI1eYGio9LzBSr73Anz4Gi//+sowrfEw79tiONIDF4wOrDdfe5mZ06Xv65bL9gvnDXsef+oZulpL23NkzJgub48OMzRUupvLpJl+74vR8PYrpcp+AI8DrxX5uBeYWnTsZInv0Z/7txd4BXhX7uuqzl/8cfPNNyu77N+/3/a5jUAz2f8zX3ha/fJXX8h/bdp++/94XP3GP7+87PZcnp5TV33mIfWXT55c8Pj5iai66jMPqW/+5K2y56/0e5/NZtWn/uGg2vLZh9Qzb4wppZT64v6T6qrPPKRiiXTZc5fL9gcPX1RXfeYh9cblmbLH/cH3j6ptv/uwmkuWt9tkpd/7WmkU+4GDqsg1tWJqSCn1HqXUdUU+vgeMiMhGgNy/RW+nlFLDuX9Hge8Ct+Sequp8TXNSSm9oMpaq29L6cvS1t7C7v32J3ESjyUuUQkT43EdzxeNvGMXj8GySoM9Na4PMPwRzUtSVCsYHToa5ZWB9Q6mOrmVqrRE8CHwi9/kngO8tPkBEgiISMj8HfhojoqjqfE3zYu4uLiSeyjCXyiyLzlAxfmpnLy+dm2SqYGJ3vMGmisthFo8zWcX9/3CI85OxmvP6ThLMLbAv10I6OhvnxMgs77hG1wcahVodwZ8Ad4vISeDu3NeISL+IPJw7pg84ICKvAD8BfqCUeqTc+ZrVQU/Iz3gkQSY7X9w0h8mWu1hssm9nL1kFT70xP1w2EWksnaFKmMXj45dneOz1kYaScTaLxeUiAnMm4p26UNww1FQsVkqNA+8u8vgw8P7c56eBPVbO16wOekN+sgrGo4l8O+lk1CgmrkRqCGDP5k7WB33sPz7KvTduAgpSQw10Qa3Evh29/NZP7+Bzj55oKK2efNdQGeG5AyfHWRfwsmtj+3KZpalArV1DGk1Jegp2F5uOYCofEazMRdftEga397D/xCiZrMLtEsajSbxuIeRvrj+H/zK4lUgizY1XdK60KXmCPrNGUDw1pJTiwKkx7rimG9cyDhRqyqMlJjR1I683VCBHPWFKUAdXJjUEMLizl8lYisPnJw2bogm6gv6G0RmqFhHhM/fs5L27N6y0KXnm5wiKRwSnRiOMzCR4p64PNBTaEWjqRl5momCobNJcSrNCEQHAXdt6cLskv6OgkYfJmo1Wb/kF9gdy9QFdKG4stCPQ1I154bl5vaGp6MoWi8HY/nXzlevYf9woGI9Hk3pLlkO4XELQ5y65wP7AyTADXQGuWB9YZss05dCOQFM3WrxuOlq9C1pIJ2Mpgj53TWsVnWDfzl5evzTD5em4jggcJugvvsA+lcny49PjWlaiAdGOQFNXekL+BTsJJmPJFSsUF7JvZw8A+0+MMh7RjsBJjHWVSyOCw+eniCYz3KnTQg2HdgSaurJ4qGwllEeLsaMvRH9HC4+8dplIIt00MwTNQKl1lc+cDOMSuP1q7QgaDe0INHXFcATzNYLJFVAeLYaIsG9nL8+cNOoEjS4v0UwEfMV3Ejx7Ksz1mzvpaICfv2Yh2hFo6kpvewujM4m8dPJkNLmiHUOF/FRuyhiaQ16iWWjze5YMlM3EUxw+P6XbRhsU7Qg0daWnzU8inWU2d4c4GUuyrkHuCO/Y2o3PY/wJ6K4h5wj43EQXDZS9cHqCTFbpQnGDoh2Bpq70tudaSGcSpLOK2Xh6xQTnFtPqc3P71V1A8+gMNQNtRbqGDpwco9Xr5qYrO1fGKE1ZtCPQ1JXCWYKczFDDpIYA7r2xn5Dfw4aO2nb1auYJ+Io4glNhbr16/Yq3DWuK01ziKpqmw9QYGptNEEkZCflGKBabfPhtm/mZGzbqC5SDtPndxFIZslmFyyVcmp7jzbEo991y5UqbpimBjgg0dSWvNzSbIJI0HEEjRQSAdgIOE/B7UArmUkad4JmThqyErg80LtoRaOpKe4sHv8fFaEFEoDt0Vjem8Fw01zn07Kkw3W1+dvSFVtIsTRm0I9DUFRGht93P6Ew8HxE0UmpI4zymFHU0YaSHnj0V5s5ruppO3XUtUZMjEJH1IvKYiJzM/buuyDE7RORwwceMiPx67rnfF5GLBc+9vxZ7NI1JT5t/QUTQaKkhjbPkI4JEmhMjs4QjSe7c1rPCVmnKUWtE8FngCaXUNuCJ3NcLUEqdUErdqJS6EbgZiGEssDf5vPm8Uurhxedrmp/eUEuuWAw+j4tAgyxa19QHc29xNJHmgFkf0INkDU2tjuBe4Gu5z78GfLDC8e8G3lRKvVXj62qaiN72XESQVKwLeHWKYJUT9BuOPpbM8MypMNf0tun23AZHzNF/WyeLTCmlOgu+nlRKLUkPFTz/VeAlpdRf5r7+feCXgRngIPCbSqnJEufeD9wP0NfXd/M3v/lNWzZHIhHa2tpsndsINKP9D76Z5DsnU1zbqZhNu/jvdzanFn0zvvcmy2n7xUiW3z0wx3+63sfXjia56woPH7+2Ni2nZn7voXHs37dv3yGl1N7Fj1ecIxCRx4Fiu/B+14oBIuID/h3w/xQ8/CXgDwGV+/fPgF8pdr5S6gHgAYC9e/eqwcFBKy+fZ2hoCLvnNgLNaP9I8BzfOfkqI3EX2/vXMTh4+0qbZItmfO9NltP2i1NzcOBJxtzdJLPD/PxdNzK4q6+m79nM7z00vv0VHYFS6j2lnhORERHZqJS6JCIbgdEy3+p9GNHASMH3zn8uIl8BHqrObE0zYQ6VTcSVbh1dA7TlagRPHB/F7RJu29q1whZpKlFrjeBB4BO5zz8BfK/MsfcB3yh8IOc8TD4EvFajPZoGxBwqAxpiKY2mvgRyNYLZeJqbruikza8FDBqdWh3BnwB3i8hJ4O7c14hIv4jkO4BEJJB7/juLzv9TEXlVRI4A+4DfqNEeTQNiCs8BDaM8qqkfXrcrr+qqp4mbg5pctVJqHKMTaPHjw8D7C76OAUviQ6XUL9by+prmoCvoxyWQVXqGYK0Q9LlJprO8UzuCpkBPFmvqjtsldLUZUYF2BGuDoN9Dm9/DDZs7V9oUTRVoR6BZFnpzdYJ1QZ0aWgv0d7YyuKMHr1tfYpoBXcXRLAu9IT9H0cXitcLffGIvHpceHGwWtCPQLAtm59B67QjWBO0tOvJrJnTcplkWzFkCXSPQaBoPHRFoloUPvW0Tly68RXur/pXTaBoNHRFoloWtPW18YKtPC85pNA2IdgQajUazxtGOQKPRaNY42hFoNBrNGkc7Ao1Go1njaEeg0Wg0axztCDQajWaNox2BRqPRrHG0I9BoNJo1Tk3L61cKERkD3rJ5ejcQdtCc5aaZ7W9m26G57W9m20Hb7xRXKaV6Fj/YlI6gFkTkoFJq70rbYZdmtr+ZbYfmtr+ZbQdtf73RqSGNRqNZ42hHoNFoNGuctegIHlhpA2qkme1vZtuhue1vZttB219X1lyNQKPRaDQLWYsRgUaj0WgK0I5Ao9Fo1jhryhGIyD0ickJETonIZ1faHiuIyFkReVVEDovIwZW2pxIi8lURGRWR1woeWy8ij4nIydy/61bSxnKUsP/3ReRi7mdwWETev5I2lkJErhCR/SJyTESOisiv5R5v+Pe/jO3N8t63iMhPROSVnP3/X+7xhn7v10yNQETcwBvA3cAF4EXgPqXU6ytqWJWIyFlgr1KqEYZSKiIi7wIiwN8rpa7LPfanwIRS6k9yjnidUuozK2lnKUrY//tARCn1P1fStkqIyEZgo1LqJREJAYeADwK/TIO//2Vs/zma470XIKiUioiIFzgA/BrwYRr4vV9LEcEtwCml1GmlVBL4JnDvCtu0alFKPQ1MLHr4XuBruc+/hvEH3pCUsL8pUEpdUkq9lPt8FjgGbKIJ3v8ytjcFyiCS+9Kb+1A0+Hu/lhzBJuB8wdcXaKJfMIxfph+JyCERuX+ljbFJn1LqEhh/8EDvCttjh0+LyJFc6qihwvtiiMgAcBPwAk32/i+yHZrkvRcRt4gcBkaBx5RSDf/eryVHUGxrejPlxd6hlHob8D7gV3OpC83y8iVgK3AjcAn4sxW1pgIi0gb8K/DrSqmZlbbHCkVsb5r3XimVUUrdCGwGbhGR61bYpIqsJUdwAbii4OvNwPAK2WIZpdRw7t9R4LsYqa5mYySXAzZzwaMrbI8llFIjuT/yLPAVGvhnkMtP/yvwj0qp7+Qebor3v5jtzfTemyilpoAh4B4a/L1fS47gRWCbiGwRER/wMeDBFbapKkQkmCucISJB4KeB18qf1ZA8CHwi9/kngO+toC2WMf+Qc3yIBv0Z5AqWfwscU0r9ecFTDf/+l7K9id77HhHpzH3eCrwHOE6Dv/drpmsIINdy9r8AN/BVpdQfraxF1SEiV2NEAQAe4J8a3XYR+QYwiCG/OwL8N+DfgG8BVwLngI8qpRqyIFvC/kGM1IQCzgL/2cz7NhIicifwDPAqkM09/DsYufaGfv/L2H4fzfHe34BRDHZj3Gh/Syn1ByLSRQO/92vKEWg0Go1mKWspNaTRaDSaImhHoNFoNGsc7Qg0Go1mjaMdgUaj0axxtCPQaDSaNY52BBqNRrPG0Y5Ao9Fo1jj/P06IffLy2ov1AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "import jax\n",
    "%matplotlib inline\n",
    "n = 33\n",
    "a = [[1.0/(i + j + 1) for i in range(n)] for j in range(n)]\n",
    "a = jnp.array(a)\n",
    "rhs = jax.random.normal(jax.random.PRNGKey(0), (n,)) #Right-hand side\n",
    "x = jnp.linalg.solve(a, rhs) #This function computes LU-factorization and solves linear system\n",
    "\n",
    "#And check if everything is fine\n",
    "er = jnp.linalg.norm(a.dot(x) - rhs) / jnp.linalg.norm(rhs)\n",
    "print(er)\n",
    "plt.plot(x)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "- As you see, the error grows with larger $n$, and we have to find out why.  \n",
    "- **Important point** is that it is not a problem of the algorithm: it is a problem of representing the matrix in the memory. \n",
    "- The error occurs in the moment when the matrix elements are evaluated approximately."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Questions from the demo\n",
    "\n",
    "- What was the problem in the previous example? \n",
    "\n",
    "- Why the error grows so quickly?  \n",
    "\n",
    "- And here is one of the main concepts of numerical linear algebra: the concept of **condition number** of a matrix.  \n",
    "\n",
    "But before that we have to define the **inverse**. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inverse: definition\n",
    "\n",
    "- The inverse of a matrix $A$ is defined as a matrix $X$ denoted by $A^{-1}$ such that  \n",
    "\n",
    "$$ AX = XA = I,  $$\n",
    "\n",
    "where $I$ is the identity matrix (i.e., $I_{ij} = 0$ if $i \\ne j$ and $1$ otherwise).\n",
    "- The computation of the inverse is linked to the solution of linear systems.  Indeed, the $i$-th column of the product gives  \n",
    "\n",
    "$$ A x_i = e_i,$$\n",
    "\n",
    "where $e_i$ is the $i$-th column of the identity matrix. \n",
    "- Thus, we can apply Gaussian elimination to solve this system. Moreover, if there are no divisions by zero in this process (and the pivots do not depend on the right-hand side), then it is possible to solve the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inverse matrix and linear systems\n",
    "If we have computed $A^{-1}$, the solution of linear system  \n",
    "\n",
    "$$Ax = f$$\n",
    "\n",
    "is just $x = A^{-1} f$.  \n",
    "\n",
    "Indeed,  \n",
    "\n",
    "$$ A(A^{-1} f) = (AA^{-1})f = I f = f. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Neumann series \n",
    "\n",
    "- To study, why there can be such big errors in a solution (see the example above on the Hilbert matrix) we need an important auxiliary result.  \n",
    "\n",
    "**Neumann series**:  \n",
    "\n",
    "If a matrix $F$ is such that $\\Vert F \\Vert < 1$ holds, then the matrix $(I - F)$ is invertible and\n",
    "\n",
    "$$(I - F)^{-1} = I + F + F^2 + F^3 + \\ldots = \\sum_{k=0}^{\\infty} F^k.$$\n",
    "\n",
    "Note that it is a matrix version of the geometric progression. \n",
    "\n",
    "**Q**: what norm is considered here? What is the \"best possible\" norm here?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Proof\n",
    "The proof is constructive. First of all, prove that the series $\\sum_{k=0}^{\\infty} F^k$ converges.  \n",
    "\n",
    "Like in the scalar case, we have  \n",
    "\n",
    "$$ (I - F) \\sum_{k=0}^N F^k = (I - F^{N+1}) \\rightarrow I, \\quad N \\to +\\infty $$\n",
    "\n",
    "Indeed, \n",
    "\n",
    "$$ \\| (I - F^{N+1}) - I\\| = \\|F^{N+1}\\| \\leqslant \\|F\\|^{N+1} \\to 0, \\quad N\\to +\\infty. $$\n",
    "\n",
    "We can also estimate the **norm of the inverse**:\n",
    "\n",
    "$$ \\left\\Vert \\sum_{k=0}^N F^k \\right\\Vert \\leq \\sum_{k=0}^N \\Vert F \\Vert^k \\Vert I \\Vert \\leq \\frac{\\Vert I \\Vert}{1 - \\Vert F \\Vert} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Small perturbation of the inverse\n",
    "\n",
    "- Using this result, we can estimate, how the perturbation of the matrix influences the inverse matrix. \n",
    "- We assume that the perturbation $E$ is small in the sense that $\\Vert A^{-1} E \\Vert < 1$. \n",
    "- Then\n",
    "\n",
    "$$(A + E)^{-1} = \\sum_{k=0}^{\\infty} (-A^{-1} E)^k A^{-1}$$\n",
    "\n",
    "and moreover, \n",
    "\n",
    "$$ \\frac{\\Vert (A + E)^{-1} - A^{-1} \\Vert}{\\Vert A^{-1} \\Vert} \\leq \\frac{\\Vert A^{-1} \\Vert \\Vert E \\Vert \\Vert I \\Vert}{1 - \\Vert A^{-1} E \\Vert}. $$\n",
    "\n",
    "As you see, the norm of the inverse enters the estimate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Condition number of a linear system\n",
    "Now consider the **perturbed** linear system:\n",
    "\n",
    "$$ (A + \\Delta A) \\widehat{x} = f + \\Delta f. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Estimates\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\widehat{x} - x &= (A + \\Delta A)^{-1} (f + \\Delta f) - A^{-1} f =\\\\ \n",
    "&= \\left((A + \\Delta A)^{-1} - A^{-1}\\right)f + (A + \\Delta A)^{-1} \\Delta f = \\\\\n",
    "& = \\Big[ \\sum_{k=0}^{\\infty} (-A^{-1} E)^k - I \\Big]A^{-1} f + \\Big[\\sum_{k=0}^{\\infty} (-A^{-1} E)^k \\Big]A^{-1} \\Delta f  \\\\\n",
    "&= \\Big[\\sum_{k=1}^{\\infty} (-A^{-1} \\Delta A)^k\\Big] A^{-1} f + \\Big[\\sum_{k=0}^{\\infty} (-A^{-1} \\Delta A)^k \\Big] A^{-1} \\Delta f,\n",
    "\\end{split}\n",
    "$$   \n",
    "\n",
    "therefore\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \n",
    "&\\frac{1}{\\|A^{-1}f\\|} \\Big[ \\frac{\\|A^{-1}\\|\\|\\Delta A\\|}{1 - \\|A^{-1}\\Delta A\\|}\\|A^{-1}f\\| + \\frac{1}{1 - \\|A^{-1} \\Delta A\\|} \\|A^{-1} \\Delta f\\|  \\Big] \\\\\n",
    "\\leq & \\frac{\\|A\\|\\|A^{-1}\\|}{1 - \\|A^{-1}\\Delta A\\|} \\frac{\\|\\Delta A\\|}{\\|A\\|} + \\frac{\\|A^{-1}\\|}{1 - \\|A^{-1}\\Delta A\\|} \\frac{\\|\\Delta f\\|}{\\|A^{-1}f\\|}\\\\\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "Note that $\\|AA^{-1}f\\| \\leq \\|A\\|\\|A^{-1}f\\|$, therefore $\\| A^{-1} f \\| \\geq \\frac{\\|f\\|}{\\|A\\|}$\n",
    "\n",
    "Now we are ready to get the final estimate\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "\\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \n",
    "&\\frac{\\Vert A \\Vert \\Vert A^{-1} \\Vert}{1 - \\|A^{-1}\\Delta A\\|} \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big) \\leq \\\\\n",
    "\\leq\n",
    "&\\frac{\\Vert A \\Vert \\Vert A^{-1} \\Vert}{1 - \\|A\\|\\|A^{-1}\\|\\frac{\\|\\Delta A\\|}{\\|A\\|}} \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big) \\equiv \\\\\n",
    "\\equiv &\\frac{\\mathrm{cond}(A)}{1 - \\mathrm{cond}(A)\\frac{\\|\\Delta A\\|}{\\|A\\|}} \\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big)\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "The crucial role is played by the **condition number** $\\mathrm{cond}(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Condition number \n",
    "\n",
    "- The larger the condition number, the less number of digits we can recover. Note, that the condition number is different for different norms.\n",
    "\n",
    "- Note, that if $\\Delta A = 0$, then\n",
    "\n",
    "$$ \\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \\mathrm{cond}(A) \\frac{\\|\\Delta f\\|}{\\|f\\|} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- The spectral norm of the matrix is equal to the **largest singular value**, and the singular values of the inverse matrix are equal to the inverses of the singular values.  \n",
    "- Thus, the condition number in spectral norm is equal to the ratio of the largest singular value and the smallest singular value.\n",
    "\n",
    "$$ \\mathrm{cond}_2 (A) = \\|A\\|_2 \\|A^{-1}\\|_2 = \\frac{\\sigma_{\\max}}{\\sigma_{\\min}} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hilbert matrix (again)\n",
    "\n",
    "- We can also try to test how tight is the estimate, both with ones in the right-hand side, and with a random vector in the right-hand side. \n",
    "- The results are strickingly different"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 8.981202871503108 Log Condition number: 21.341344339184218\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "DeviceArray([-2.59701783e-03, -3.62615992e-01, -4.92682651e+00,\n",
       "             -1.39317220e+00, -8.60041880e+01,  3.03460550e+02,\n",
       "              9.20961751e+03, -3.41310628e+04, -8.24928647e+04,\n",
       "              5.27596100e+06,  3.30516690e+07,  5.81349493e+08,\n",
       "              3.30494319e+09,  5.91765636e+10,  1.54937154e+11,\n",
       "             -1.03870175e+13,  6.97324076e+13, -1.97562227e+13,\n",
       "             -5.47702944e+13,  1.94494972e+14,  2.37939836e+13,\n",
       "             -6.35816689e+13, -1.25422856e+14,  5.42975038e+13,\n",
       "              7.27578993e+12, -3.22040681e+13,  1.57513412e+14,\n",
       "             -1.41688598e+14, -1.34967174e+14, -3.54499321e+13,\n",
       "              3.86080313e+13, -3.02898541e+13,  5.14933533e+13,\n",
       "             -1.09592987e+14, -1.92757750e+14,  2.72745301e+13,\n",
       "              7.99814548e+13,  6.06250212e+13, -1.03983362e+14,\n",
       "             -8.07437378e+13,  7.07436338e+13, -1.43980174e+14,\n",
       "              4.39149056e+13, -8.82456882e+13, -1.58302961e+13,\n",
       "             -1.13785201e+14,  1.38974003e+11, -8.10241429e+13,\n",
       "             -6.89874139e+13,  7.33904252e+13], dtype=float64)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "\n",
    "n = 50\n",
    "a = [[1.0/(i + j + 0.01) for i in range(n)] for j in range(n)]\n",
    "a = jnp.array(a)\n",
    "rhs = jax.random.normal(jax.random.PRNGKey(10), [n]) \n",
    "#rhs = jnp.ones(n) #Right-hand side\n",
    "f = jnp.linalg.solve(a, rhs)\n",
    "\n",
    "#And check if everything is fine\n",
    "er = jnp.linalg.norm(a.dot(f) - rhs) / jnp.linalg.norm(rhs)\n",
    "cn = jnp.linalg.cond(a, 2)\n",
    "print('Error:', er, 'Log Condition number:', jnp.log10(cn))\n",
    "\n",
    "u1, s1, v1 = jnp.linalg.svd(a)\n",
    "cf = u1.T@rhs\n",
    "\n",
    "cf/s1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "And with random right-hand side..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error: 17.674315761144477 Condition number: 4.073996146476839e+19\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZoAAAEXCAYAAACUKIJlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZwcV3ku/Jzuqt5mpme0jhZLli3b8m5jGxMbsMd2TJILSchyiROSLyQQsl1yA+T7Em5I4EJ+IQmQ5XKzkRtIIMQQIECSLxgbrLHl3ZJl2fKi1VpHM1pm67W6lnP/OHWqTlWdWrq7emak6ef3009Sd3X16apT5z3P+z7v+xJKKfroo48++uijV8gs9gD66KOPPvq4sNE3NH300UcfffQUfUPTRx999NFHT9E3NH300UcfffQUfUPTRx999NFHT6Es9gAWG6tXr6Zbtmzp6LO1Wg0DAwPpDug8wHL83cvxNwPL83f3f3My7Nq16yyldE2SY5e9odmyZQt27tzZ0WfHx8cxNjaW7oDOAyzH370cfzOwPH93/zcnAyHkaNJj+66zPvroo48+eoq+oemjjz766KOnOK8NDSHkCkLIxwghTxFCzhBCKoSQ5wkhv0sIWV5O1j766KOPJYrz2tAA+EUA7wdwCMDHAPy/APYB+AMATxBCios4tj766KOPPnD+iwG+BuATlNI54bW/IYQcAPC7AN4N4H8vysj66KOPPvoAcJ4zGkrpTp+R4fiK/fe1CzmePvroo48+gjivDU0ELrL/nlrUUfTRRx999AFyobUJIIRkATwG4BYA11JK90mOeS+A9wLA6OjozV/+8pc7+q5qtYrBwcEuRnt+Yjn+bv6bT1YsVHWKbSuziz2kBcFyvtfLCZ385rvuumsXpfSWRAdTSi+oPwA+A4AC+FCS42+++WbaKbZv397xZ89nLMffzX/zr/3TLnrPp8cXdzALiOV8r5cTOvnNAHbShOvyBeU6I4R8HMB/A/BZSuknFns8fVx4mG200NTNxR5GH32cV7hgDA0h5KMAPgzg8wB+ZXFH08eFikrTgG5aiz2MCxrb953GL31hJ/dQXNCYmm/CtC7833lBGBpCyEcAfATAFwC8hy6HGdrHoqDSNGCY/enVSzz72jQeenkK9daFzRznmzru+JPt+P9fPLXYQ+k5zntDQwj5fQAfBfBFAL9AKe1vN/voGSpNHa0+o+kpNINd3/mmvsgj6S1mai1ohoWzFW2xh9JznNcJm4SQXwfwPwEcA/BdAD9DCBEPmaKUPrQYY+vjwkSlacA7xfpIG5rBmMxcQ8f64Qu3uEdNY79zObhiExsaQsgdAOYppc/3cDzt4vX235sB/KPk/UcA9A1NH6mgZVjQDAvZTN/S9BKazhbeufqFzWhqLQPA8jA07bjOtsPOPVkqoJS+i1JKIv6MLfYY+7hwULFdOaZFYS2DAO5igbsm55vGIo+kt6hp7Pe1lkHMrx1DcxZAo1cD6aOPpY6KsPDp1oW/C10sOIymcWEzGi526DMaL8YB3N6jcfTRRyS+98oU3vG3Ty4qk/AYmmWwC10siDGaCxmc0ehG39CI+DCAbYSQjxNC1F4NqI8+ZHh0/xk889q0o0iKQ6Nl4n3378bkXDO1MVQ0d+FbDovDYsFRnS0XQ7MMGE07qrMPAdgL4H8AeDchZA+ASbByLyIopfTdKY2vjz4AAMdnmNe2ZVooIr7O2KEzVfz7ngn84DXr8Nbr16cyhr7rbGHADc0Fz2hs19lyiNG0Y2jeJfx7nf1HBgrWB6aPPlLD8ek6AKb8SgK+SzRSNAh919nCgLvOLnRGU19GqrN2DM0lPRtFH0sWp+Ya+L1v7sWn33EjhouL4zGllOKEzWiSPpTcECQ1TElQafZdZwsBLga40BM2+3k0ElBKj/ZyIH0sTTxx8By++8ppvHhiDm+6fPWijOFstYWG3t5DyY9Lk3l4Gc2FvzgsFhbaddayc6MWOj9qOcVozvsSNH30FqfmGJM4XUkvqN4uTszUnX8nfShbPXCdVbW+62whsNCqs5/7+6fx8f94eUG+SwSXN7eMC38utV2ChhDyfQDeA+B1AEYAzAHYBeDzlNIn0h1eH4uNk7PMwEzNL149Ji4EAJBYdcZdWz1znS2DXehiwVWd9T5hs2VY2HV0BgP5ha/GtZwqA7R1dQkhfwCmPvNzzBsB/CIh5I8ppf8jrcH1sfhYCoyGCwGA5EyCH2ekmHcz33edLQgWMmHzyLkaDIvG9hiqagZqmoHRciG1764voxhNYtcZIeS/gkmbj4ExmksBFO2/32O//tuEkHf0YJx9LBImZrmhWTxG04nrzInRpMpoXEPTr+DcO7RMCxkCNHQzVUYqw77JCgA4McAw/PlD+/HO//N0qt9d7cdopHgfgCkAr6eUfo5SeoRSqtl/fw6swOUZAL/ei4H2sTg4ZbvOziym62y64VRMTmo4Wo4YIMUYTVNHUWU5PP2eNL2BYVowLYpVg3kAvVee7Z+yDU1M75uzVQ1nq+k+A1zevBzyaNoxNDcA+Bql9KzsTfv1r4K50fq4ADDf1FGxd11Ti+k6m6lj4wgrF6+1y2hSdJ1VmgZWDuQ85+8jXfD4zBrb0PTafcYZTZzrrGVajksvLfCEzeUglW/H0CgA6jHH1HGe97jpwwVnM+vKBZye1xalta5pUUzMNnDpmkEAyR9KflzarrO+oektHEMztDCGxmE0MYZG0y00DTPVZ6Avb5bjIIC3EUKkn7Ff/y8ADqUxsD4WHxO2EOCGTcNo6KZH3rtQmJxvQjcpLl09AGBxxQCVpu4YmuXg7lgMcGnzWtvQ9LI6QKNl4qgtNIlznWmGBUrTk7VbFu1Xbw7B/QCuAvAtQsjl4huEkK0AvgbgagD/nN7w+lhMcCHADZtGACyOIIArzrautRlNm3k0aQXtLUpRa5mOoTGWweKwGODuqbXl3jOag6eroBTYsqqEZoxbjIsSmka0QUoKkUEth5ysdgzNnwJ4FMBbAbxCCDlGCHmaEHIUwD4AbwfwuH1cH0sUTd3EuYRBzVOzTWQzBNdtHAYATM2nF6f50L++iLs/NY6/Hj+EMxEGzDE0NqNJqkJqpew64ykdfddZb+GP0fSS0eyz3WY3bBpBy7QiNw+caaUVp+E5NMDyUDAmNjSU0haAewH8LoDXAFwEpjTbZP//dwHcYx/XxxLFJ7+zD+/42ycTHTsx28C6cgHrh1nuQJRBaBfPH5/FydkG/viBV3HbJ76HP31wn/S4EzNMcbZ5VQlA8ofSLaqZzm6xYWdv911nvYXjOrPzVXrJaPZPVZBTMti2bggA0IzYlHADGCcaSApe56xcUJbFpqWtEjSUUp1S+glK6eUAymBGpkwpvdx+/cKugncBYPexGZytJtsLTMw1sH64gDVD7KE/naLEudEy8APXrMN3P3Anrlw/hG/vnZQed3ymjnXlAgbtzO1282jS2i0GGM0yUAotBviCXi6oyCuZnrZz3jdZwWVrBjFkz60oI8IZctLKFHHgQoCRUm5ZzKV2EjbvIIRs5v+nlFYppScppVXhmE2EkDvSHmQf6YBSigNT1cSL9am5JjaMFFEuKCiomVSrA9RaJkq5LC5bO4ir1pU9yZAiTkw3sGlFCWqWTdXkbQIY40jrIa7rXkaTZg21Plxw11RezWC4qGKu3ltGs23dEAp2blSUICBtRsOFACMltR+j8WE7vD1pZPh/7OP6WII4NddERTMSGRrLojg128T6kQIIIVg7VEi13lmjZaKUYzvJoYIaqmg7PlPHRSuLyClsqrZfVLM3rrOoxaHeMrDr6HQq37vc0DLZApzLZlAuqj1znc03dZyaa+KK0SEUc8zQRBkRJ0aTFqOxYzTDRRUt01qU1IGFRDuGJkkNbYJgx80+lgh4zoBu0tiJfa7WQsu0sGGYJUquHcqnxmgopai3DJTsB3ywoKCqGTB9RkEzTEzON7FpRQmKXcI9aWzEyaNJ2XW2omTHaCIWnK8/dxLv+NunLvh+Kr2An9H06hoesJ+FbesGnWoPUbk03MBoqcVo2ITi8ylNGf5SRNptAjYDqKR8zj5SAjc0QPzE5tLmDXZG/tpyPjV5s2ZYsChQyrMHvFxgzEZU4rAxNEEpsGllCYQQ5LKZDvrRpGVo2PUqFxUoGRJ53vmGDtOioe7APsLBF/S8kmWusx4xmn2TzON/xeiQa2gSuM7SYjS8oOZIiTUTvNAFAZFZ/ISQ3/e9NEaIlNhkwYzMfQAeS2dofaSN/VNOOA26aTlxDxl41WauOFs7VMCO/dLqQ22D+6dL9gPOA/2VpoFywe3iyaXNm1YwY6dmSfsxmpT833VuaAoq1Gwm0lBzF0x9ERJcz3dwF1VeyaBcULB/qjeGZv9UBQO5LDaOFB2XcBijoZS6eTRpMRp7UzVid63VDQrkUjn1kkRcuZiPCv+mAMbsP2E4CeB3uhpRHz2DyGjiJjbvQ7NRYDQVzXAKAXYD7jYQYzQAUPUxgON21eZNK5m0OackZzRpF9Vs6MzQ5ZUMlBiDx3e9tZhs816jqZvQDGvRWnB3ApfR2K6znjGaCq5YNwRCiMNowoyIqFxMK2HTFQNwufwyZjQA7rL/JgAeBvAPAP5RcpwJ4ByAfZTSC/uKnaewLKY4yysZaIYFPUY1dWq2gYKacaj92hQlznznyF1ngwXOaLyLyomZBtQscXqAqB25ztITAwzmlUQuvKXCaD75nX344lNH8f7vvwK/9OZLoEQw2KUCN0bDXGcVzYBlUWRSbLNMKcW+qQruvWoUABwxQBijEd1laSVsVjUDapZgwH4GlrXrjFL6CP83IeQfAXxTfK2P8wcnZhpo6CZu3DSC54/Pxk5sLm3mrtJRuyRIGnEal9Gwh2yIGxrfwjxTa2FFKef0clezmcRtb3sRoxkq5JxxJDI0i8xoXpqYA6UUf/zAq/iPFybwJz95Pa7ZMLyoY4qDx3VWVEEpc6kOl9JjZWeqGqZrLVy5niVqujEa+T0VjUtq8mbNQCmnOO7rC93QtFMZ4Bcopf/Wy8H00TvwchvXbCgDsF1nETg523AUZ4DAaFJQnvGgq+M6E2I0IuYausftk1MyySsD2L8vrZpkdcM1iKpCIvvRuK6zxWU0J2Ya+C/Xrcdfv/MmTM1reNfnn13U8SSBZrCmZ0qGoGzf+7SVZ7w1AK8IEKc6E+dcevJmEwO5bN/Q+EEIuZsQ8jlCyIaQ9zfY74+lNro+UgOPz1xtG5q4BfuUXRWAg1fTTcN15ogBHEYjj9EEDE0200Hjs/RcZ46hyUQbvKXAaAzTwqk5Jg3/oevW4xfftAVnKlpsleLFhmZYyCtZEEKce5+28owbmivXsWehkGPLYBhbESXNccU3k6LeMjCQV4RE5L68meM3ANxOKZ2QvWm/fhtYJ84+lhj2T1WwYbiAlaX4zPaWYeF0RXOkzQCTYeaymVQaoPGdvphHAwRjNHMN3dnVAoxJLGYJGm4Q411nNqNZxBjNqbkmTIviIluxt8pONJ2uL+1ShJpuIq+yZYkrENM2NK9OVrBmKO8k3+ayGWRIhKERYzQpiQGqmolSXkFOYW7hPqNxcROAJ2KOeQzALZ0Ppz0QQj5ECPkqIeQwIYQSQo4s1Hefb9g/VcUV64Zcqh6xg5qaZ/krG0ZcRkMIwZqhfCotnf2us4FcFoQgUB3Az2jUbBuuM14ZIC3XmU49rrMopsQXo8VkNCdmmDz9ohVMsccTA6cT1rlbLLRMCzl7jvJ7343y7ImDZwMbmH2TFVxpu80AOMqzMLYnKgxTYzSa0XedhWAtACmbETBlH7dQ+EMAd4M1W5tZwO89r2CYFg6drmLb6BCULM+wD5/Y/mRNjrSSNms+1xkhBIN5JT5Gk80sWh5N06ROLCkpo1lcQ8Ol4TajGWSG5lxt4XsKffTfXsL/+t6BRMdquuUwGi4A6JTRHDlbw8/8n6fxd48edl4zLcpqnI0OeY4t5rLJVGcpMZqaXYLJcZ31DY2DObBqzVHYBKDW+XDaxlZK6SpK6b2IN4LLFken62iZFi4fHXJ2i1E7/Um774wYowHSK0PTaHnzaADmJhENDc+sL/vEAInzaFIsQUMp9brOMklVZ4vnOjtut1dYbws6Vg6wGNvMIrjOnjp8Dk+/di7RsTxGA6DrGM1DL08BAHYcdBONj5yrQTMsRwjAUVCjDE36MZqaZmAwLzKafoyG4xkAbyeErJO9aYsE3m4ftyCglB6OP6qP/VxlMzoEVYmf2LxXxlDBKykdLadTWLPeMqFkiFMoE2DVAaqau6BwYYDfdZb0gUwzYbPeMmFRtOE64zGaxWU068oF5xrz2Ny5RXCdtUwrMbvTDBN5e8wDuSyyGdKx6owbmj3HZ51zcCHAVevLnmMLajY8YbMHjKbeMliMxnFl9xkNx2cADAHYQQj5EUJIHgAIIXlCyI+Cdd8cBPC/0h9mH91g/1QVhACXrR1M5BPmD1xe8U6PtUN5zDX0rpt+1VumkyTHMVTwus74LtZraNoXA1gUgWKd7YKPK6kYQFsCjObEDGuvwFEuKshmyKIwmpZhObW94sAYDZt3hBCUC0pHjOZcVcPOo9N4wyUrYVHg6cOsmvarkxVk7GdBRFSMhm8cshmSIqOx5c3LRAwQVxnAAaX0QULIxwH8HoBvAKCEkBkAK8AqBxAAH6OUPtCTkaYIQsh7AbwXAEZHRzE+Pt7RearVasefFdEyKXLZ9DKf/Xh8bxOrCwRPP7EDR+fZw/TcnhdAJuW3/+XDbDF65snHPeOanmAP/KmZWle/+9BRDQo1Pedo1Zo416LOa0fm2DiPH3oV45WDAIDZc03MzluJvrspsKPvbR/v6vqerLJF4Pjh/RhvHMbcTBMzDRo6jkqduRdPTJ5JZX50goOn6rhyZdbz/QMKsPfAUYzn5U3mZEhjjlfrTWhNJDrP1NkGssQ9NgcDB46cxPh4Mtcbx44TOiwK3Dtax+5jwFfGn4d6Oo/H9zaxtkTw1OM7PMdr9QZONdj3+n/z8xNsw1DMUkyePtv19bAoRUM3cXriBHZTdi/2vPgSiufkXWYXAmmtZWFIbGgAgFL6EULI42AS5jcAGAEwDeApAJ+hlD6U/hDTB6X0swA+CwC33HILHRsb6+g84+PjGBsbw+RcE7ONlqPLbwf7pyr4ob/YgYfefwcuXTMY/4EO8OcvPY5tGxWMjb2B5dM88SiuvOoajF2/Xnr8bn0/sP8A7r17DJ4iqvtO43N7n4WeLaLTawYAX514Diu0ec85vn5qN/aenHNee+zAWeDJp/HG19+EWy9ZCQD4t9PP46Q2nei7rYe+jVyWuW1ue+ObAm7AdvDcsRngsSdw603XY2zbWvzLyV2oTFUxNnan9Hjz4QcAmCgMljE2dnvH39spWoaF2e98G7dcuQVjY9uc19fvfhT5cgljY8mFoXyOdwPy6IOgGZLoPH+29zGsGMhhbOxWAMDo3seQL7n/T4ovfWEnNgzP4d0/ejd2TD+LI7MNjI3diY8+ux03XVrG2NjNnuM/d/gZzDV0jI29MfCbp549BrzwIlaXSygNFTA2dltbY/Gj0tSB7zyIa7ddhjdduw54bDsu23Ylxm6+qKvzdoM07nMU2i5+RCl9kFL6w5TStZTSnP33j5wvRqYX+LOH9uPXvvRcR589MVOHaVGctJVevcDJ2YZTHDOR68wwkVMy8FfqXj3IAsrzre5cUY2W6dQ54/CrzmSus6RtAiilaJmW456LyuJPAj4u3s4g1nW2yDGaybkmLApctLLkeX3lQG5RXGe6SRNfC9F1BqCjwpqNlokdB87g3qtHQQjBmy5bjYOnq3jtbA1Hp+vYNhrcEBbVTGivGae9dFFNJUbjJCzn+/LmPtrAfFPvWOvP6yj1KmNbM0ycEZIv1QTyZk23UFCCU4Mv+jW92xiNgZLqJdPlguLJd+CGplx0j0sqBuAl/Ady6RQs5OPyxGhCgre6aTkxocWK0fCq1zxZk2PlQA7naosTo2noJqwEsTJRdQawxb3dZ2vHgTNo6hbecg3TLd1+2SoAwD88/hooRUBxBtgxmhgxQLmgphKj4Ym8AznFeR77hsYHQsj1hJA/IoR8ixDyXeH1LYSQdxBCVqQ7xKUPzbA6NhROL/II1YlpUfzEXz+B7ftOt33+U7zc/4rkjEYzTOTVbOB1ntdQ6zJRux7CaDTDch5qKaNRkuXR8N/GGU23OQquGEBgNCGLpqhcWqw2AU4OzQoJo1lgQ8PZJZCsxL6mmx41YrnQfvOzB1+eQrmgOC7Xq9aVsXIgh6/uOgEAnmRNjmIuXgxQLirpMppc1lGBJs0PO1/RlqEhhHwMwHMA/j8APwy3jQA/1/0Afja10Z0n0AwTDd3sqO83n7jNiEWpqhnYdXTGUc60A+6S87vOotxJTd1CQQ1OjcGcggwBal3WZaq3TCdZk4Mv4rw6wFxDh5p1e4XwsScxGrzqwYCdYNm968zPaMLVb3xRymUzi1ZX7MRMA9kMCeRBrRjIYdbu/rlQEBloEomz1HXW1BM/W6ZF8fCrp3H3lWuduZ7JENy+dRVTO6pZbPa5FIGYPBr79cG8kgqj4XN8UJQ39/NoGAgh9wH4MICHANwI4BPi+3ZOy04AP5LmAGPG9HOEkA8TQj4MYA2AYf5/QsjPLdQ4NJ21JvYvgueqGh7Zfyb6s0b8bo8bo7PV9nNYgoYmnqprhomCEmQ0mQwrdFjv1nWmGSj6XGeDvsKavCqAGCfK2Qt83KLD70MpJddZtWmAwO0IGuU644xm5UAOtZbR0eajWxyfZjk0/v4zqwZyoHRhkzbFZyKJxLnlc50NF1XoJg01An6cnGlgutbCbVtXeV5/02WrAQBXrBuS9rYpRuTRaCYzfgU1mxKjsROW8/02ATL8BoCDAH6UUvoCANlsfQXA5WkMLCHeDeDj9p+1YCo4/v93L9QgHNeAr5/Fl589jl/4/DORPSz4e1HH8DhOR4bGzhBfN+w2DxPHLB+TWwbEj+Gi2n2MRjedhk8cnNHwxLp5X0FNgI2dJsiL0R1Do9j/7268800DBQXOAhUVK+I73hX2op5W3kU7ODHTCMRnADhFJBfSfSYa5LoeH7PSDO/c4zG6+UayeBffsA3mvXPnjbahuXI06DYDmKHRTSpd8DXdNTTpxGjYGHlCaob0DY2I6wB8h1IaNUsnAIx2N6TkoJSOUUpJyJ+xhRoHNwR+VjLfYFr+MxH1wVwxQDTDADpnNGuH8o7fO5nrTM5oAG5o2h6GB3VNkrCZ97rO5pt6oAUx/w1x7jM9ZUZTaRooKT5mZcmZlcto7HjWIggCTsw0nPbXIrihWUhBgIfRxLjOLIvFc/K+GA2QvCeN06HTJ2bZtLKE37jncvz0GzZLP1eIaOesGRZyStbuTtuZi1yEyGiA9orFhuEvtx/Ee7+ws6tz9BLtGBoCIO5qjALovhjWeQZuCPw+eU73owxEEtcZ30VFGawwTAjSZgCJdlCaYTkPnh/DpVxXjMYwLbRMCwM5r+uMxz8qPteZiCSVp4FeGBodgvgNSgSz4veT1xZb6DiNZpiYqjSXDKMRg9xxrjO+2PpVZ0DyCs5Oh04JI//AvVfgxk0j0s8VIto5twyX0Vi0e4Zcta/DoP0MsD5L3Z1z19EZ7DhwdlFctUnQjqE5ACA0+4wQkgXwJgAvdTuo8w18cfFPUr6DOxtRXyrMSMnOf67aSiQRFXFytoGNPvVR3A6qqZuBHSFHt66zuu6t3Mwx6IgB2IIiNTQJGQ1vIpWW66yqGSgKjCaqECIPHK8sLQ6jmZhlLR78ijNgKTCa6GshYyNlp1dRsuvIn5V8CCMPAxed+N3f7JysRw4fVxL1XBTqNmvnrF5to1hsGGbrLTR0M3KtWUy0Y2j+BcBNhJAPhrz/IQCXAfjnrkd1niHM0PD/n0vAaKKCjHzxMizaltTTsihOzTY9fWUAHsyOLgoZymiKSldiAL6rldU6A6IZTS5BDhDgMhgeB+paDBAwNOHj4IvQCntRX+ikzRMhOTSA25NmsRhNXEBfxkbabefMz5EL2SiFwTE0kudQM1iPHC7517qM09RaJnLZjODOTl7DLwyzdXZ9eA7VUkM7d+PPAewB8CeEkKcB/BAAEEI+Zf//f4KVovls6qNc4uCGwC9RbrQSuM54fCdi8or9MNqJ05ytamiZFi7y9ZVRsySyw2YUoxkp5lAz0DFF57tav+tsMO8aGsuiTAxQkMdo4irdBsUA3RuagjBcPg5ZqwV+P3lHy4VO2jw+bTc8k8RockoGQwVlQRmN3kaMRsZGnBhNUtdZSIwmDkW7nbPMs9AyLOTVrMtoEirgwsAqN7cv248CVxIenz7PDQ2ltAGWN/NFsG6bt4LFbT4A4GYA/wTgBymli1eydpEQymhSc52577UTpznBpc0r/IYmvnGXLGETYK4ziwa7YSYFX2z8jKagZpHLZlBpGqi2DFgU4TGaONcZZzROjKY711lNM1AQGI2SCXed8R0xj9EsdPOzEzN1KBmCdeWC9P2VAzlMLxKjiWtt7TAawUi4asRk882N87RnaDiDl7EuzTCRz2acY7QukyurmuHZaOXaaH8hgyV4Onhn1aWGdotqzgF4FyHkAwBeD2AVWEO0Zyil0QkjFygM03JKngRiNPb/z6QkBog7lx9hnTLVbMaJY8jHZEoTNgFvM6pOClXyhdfPaAAWp6lqOubqwaoAfNxAEteZHaPJp8NoapqJgmAXo3KRXHkzG/tCM5oTMw2sHykgK8kVARa+3pnHdRZjdPm1E91eBTWLnJJpn9GEbJTCUBQMjf/KaYaFwbySHqPRvPL+qLysJKg02cYMWLqMpi1Dw0EpnQbwnZTHcl5CXPT87i/uSouK0fBJm5TRtBPsOznjTdbkiPMJs7yBMNUZW0Bn6zou6qDYEF94/YwGcHvSuHXO/DGaZOU6+EPLBQdRbsI4WBZFrWWgqHhL4QByg6cJCZtA72I077t/N65YO4j33eNNW5tr6E4sRoaVpRxOzS2cMNQjBoiN0cjZSLmgth2jad91xsUAJvzRrZZhIT+QHqOptQxPd1nWSC/ZOXXTQoYQz0ZituGuCUuV0UtCFuYAACAASURBVLRTGeArhJAfIoT0C3EKEAODQUbDFtVo11l8rbNOYzQnZxsoF5QA81CzmdDF17RzGeIYTadFRMU6T34M5hVUm4Zz7rA8mjg3gyMG4DGaLqSjdd0EpfC4zqJykZqOvLm3MZpdR6bx4sm5wOtcihuGxWQ09RjXWStEMVYuKoldZ2HGKg7FSNeZhZziqs7CqjwnRb0VZDRJYzS/9IWd+J2vv+B5bcb2AAzmlQtCDPBfAfwHgJOEkE8SQq7r0ZjOK4hGICgGiM/od1RnCaoHDBUUnG0jRjMx2wi4zYBo1xl/2MNVZ+31cf/W8yfx6uS88/8o15mf0XQbo0mjqCaPK4iuMyUT5TqzVWel3jKaimZIf1fLtCIVVysHWQXnhcq3aKfWWVgOTLmQvIKzU2su7RiNkk2P0fhiNHExU46WYeGJg+dYTykBs/bG4ZoNZUzMNha0ll1StHM3bgPwtwByAD4I4HlCyE5CyPsIIat7MrrzAFHyzYa9m52t6xFFGE3pZz3H2Kxp40ixrRhNWCmSKNdZWBtnjhHuOkv44H/4m3vx+ceOOP+Pcp0N5lVUNMHQlPyGJpm8md8T/jDL1GFJwUUPHkYT5TozLKhZAjWbiSw93w0opahphtSF2LKluGFYWcqhZVgLVlm6ZbLvyWZI566zopqc0djfEXUNZHAMTZjqTGA03cZoai3DKfgKJBcDvDo5j5ZpBTwkXNp8/UXD0E2KyfmllzPfjursaUrprwFYD+AdAP4TwPUA/gKM5fwrIeTthJCO4j7nK8T4ibioULtdK5e5ngtxn7ny5ghGY5jIZgjWDRfadp354zNAtOtMS5HRUEpRb5k4XXEnvsNo8sHz85403B8fKgaIlTdzMUD3qjOH0YjyZqdCgZzR8PI9pVw2VmnVCZp2EVfZzprvvsOw0NUBuNtypKjGigHCki2HCgoqSRmNXcLG37QvDsXYEjTtxWiOnqvhd7/xovR8dc1bvTxpHs2e47MAmCBIZKSc0Vy7cRjA0hQEdNJhs0Up/Rql9IcBbARjN68AeDuAr4PVO1s2ECed+CBphuXpchhmIBxpdJQYwG5EtmYwj7OVZAvEfFNHpWmEus7C4hZxjKaoZpEl7i4qCprBmoCJLIwbGlktNaY6Y4wmmyGOPJkjryRznQVK0HQhBuDVpGWVAQyJi0KUhpfy2Z7ImznLCmU0MTEaYOGqA2j2vRguqvHy5pC515YYQI+OUYVBzbIAu9R1prfHaKqagff840586elj2DdZkb4vMhrmyo6fo88fZzG5lmGhIlxLHqO5oAyNCErpGUrpnwF4HYDfAmCASZ6XDURGI05A/u9Ntusq3NDYyZ4xYoC8msXqoTzO1bREZWgmQnJoAOb6CXM/cZl1GKMhhGBAJYkYDTeeYu4PaxGQlZZqF2M0/hYBQPIYjWNo1O7FAFVZjCbChacJya4DOaUnjKaWgqFZKEbDxzhcUhNUBuDSZL/rTElcvZk/K+2CEGK3CpDHvZLGaCil+K1/2YMDp6sAgkbJMC1ohrfWX9ISNHtOzII/NqKHZK6ho1xQsGlFCYQAx5eg8qwrQ0MI2UYI+UMARwF8EoAK1kpg2UBUnYmTiu9kNzmMJsR1Zk/almGFGhCeqb96MA/dTFaGZmJWLm0GADUT0bhL566z8KkxoCZTnXGf/FmhRpusRQDHYF6FaVFMzmlOjSvPuJ3KAMn60eTVTNcl2HmtsoKnenO464yV72Hvl3K9idFw4ycrW6TFGJpVdiLpQjEafu1HimpHlQEAxmhappUoNsJch50ta7LmZ4bdmpu5zuIZzV+NH8IDL03irdetBxCeWyc+A0liNJWmjkNnqrj5YpZTIG5cZ+otjJRyyCkZrC8XcOJCYDSEkBWEkF8lhDwF4GUAvwOgDODvAbyZUrot5TEuaXhcZ8Kk4v/eHOc6Ew1VSNImrz22ejAXeS4RYTk0gB2jCe2nwt0X4bvCAZV4tPth4GII06KOpLauGVIhAOBmgZ+cbQTiM2zcbLHX4hiNbYjUbMZuu9yNGIC7+sRxRFQG0E1n51vqEaOJc51FLbQ8kXS61n4l8E7AxzhSysXKm8NyYNqpdxZnaKNQzGUCylFRoMCfiTBG8/zxWXzqwX344Rs24Nfu2gogaJR4rT9PHk2CGM2LJ+dAKXDPVawLi6g+nanrWGELZy5aWVqSuTTt5NG8jRDyVbAYzP8GcAuA74K1bl5HKX0vpfTx3gxz6UJzXE0ZNMScGnvCrhrIoahmpbJkSimahukssGH1zvgubc0Q240mUZ6dmG0gl2UsyI8oqu6KAaIYTTLXmbiD5WOut0yptBlwDc2JmXogWRMA8lk75pKg1lk2w3zu3ZZg54ZCVlRTJqhoCjvqUq43MZqaw2gkjCpG3szbB09321QoIViCIfveWNUZrwyQ9cdokjc/i0o2joNMJdgSDE02Q6BmSSijefHELCgFPvzWqxxD4n+mOUMOVAaIMTR77PjMPVeuBeDdbM7VWxi25fSbVpSWZC5NO6b/3wD8BJib7MMALqaU/gCl9J8ppUtPT7dA4A/7cFH17Ib4hC3mslhl5y74oZsUlLqS4TA3S9MORq6xjUaS6gATs02sHylIYyFqloTHaBIwmpKaTHUm5pCcnmcPRkMPNj3jECs4SxmNEt+Gmr/PjYESU0A0DtWmgQwBcsKTEqV+03S38vVAXgltEzDf1PHeL+zE6Q6kqGGMhlLKGE2EtJcQghUD6oIympySSSSM4FWS/XO2neZn/sZp7UBmaNy8HHZP80o2lNFw9lsuqKEJoHwDWlS9hiZODLDn+Cy2rCrhktUDIAQ4I6wBIqPZtLKIyflmKi2n00Q7d+SzAG6nlF5JKf0EpfRkrwZ1PoHvwoaL3mCnmAG/ejAvdXfxyTBSZLuR0J7lhsnEANzQxCRtWhbFc0dncNmaQen7aiZ8B9VMwGgGVZJIddYQWvdyQYA/Wc1zXqH9rtx1ljxhkx+bNBkuDFwhJAoTIl1nhug6y4Y2+9o3WcGDL0/h2SMzbY+JG3C/C9GNTUXv6FcO5BeM0bC8ogxKqoKWYUXmNIW5/dx2zglcZxGVx+NQULMB9affnVdQM6HPab3FNiUF1Y3n+M/nrgtCHo0SH6PZc2IWN2wagZLNYGUp5ylrNVtvOQnCm1aUQCnbaC4ltJNH8yuU0qd6OZjzEaKx8MRouIzXNhCyqssiGwKiDA17AIeLKpQMiXWdPX7oLE7ONvBjN22Uvq8qJDRGoyVhNApBpWnEZiCHuc7iGA0gNzQ8Iz8+j8ZNWowrIBqHmmY4LQw44opqJnGd8fnRSTkY3hyuZXjbSfPrEpesuGogt2CMRrcZBpeaR7nPeIMxP1xGk8B11qHqDGDPqv8ZdFxn9riiGE1NY25hQkhopQHX0+H+zrgYzdR8E6fmmrjhItYddNVgztm4GqaFecEDwMVHS03i3K9b1iUcY1HyJqTx3Xwpp2DNUE7q7hI/C4Qbmqbtd85kCJtkMYzmy88ex4qSinuvHpW+H1Vbqel7sGQYUNlCW4lxZXgMTcU1NLI6ZwA8C7rM0BBCkFMyaMXVOjOowGi6c535s7gBQf0mrQzgFQM0dFNqkPm97kRmXBVYkngfk5ZfWbGArQJanNHkwzPvOdiGSpLIW+RtvpOJAXriOrPnUz6G0fDfyZJGg6WlxA0oB0ugpqGqU56oecMmlifDPCTs/nHj64gB7HSGpRan6RuaLsEn4khR9VYJsOucFdUsVg3kMS3Jf+GT0GU04WIATsXXDMndcBzTtRYefGkSP/a6i0JZSZQ7iY8pLI8GYPJmID5pk6uMhgqKz9DIXWdiozOZoQG4FDRBjMaO53TrOqs0JYwmqh+N7sqbecBXFnvjBn0mgQvSD1HJJrI7MXAdhVULaGh0W5zgMJpYQxPBaJKIAQyzC9WZzNDw+mtJYjSuW5jn5QQZjbsB5XBcsSEboj0nZpHNEFyzQTQ07HnijHjEdp2NlgtQs8RpgLdU0Dc0XYJPurKvxIZY02v1YA4WDbpJRCMFhO/2RCWNuJuR4Ru7T0I3KX7q9ZtCj1GzEa4zHqOJkTcD8YIA7ibZsmpAMDRGKKMRlThhhiaJFFSM0ShdNpVq33XmloDhi4lM1uswmg5cZ6Kh0SSGJpbRlHKYbxpd9+lJgpbtxizaybNRcm9NlxuJgpqBmiXJ5M0dVgZg3xNM2PTXX4uO0Zge9ivLyxE3oBy5iJgfwBRnV64bcjZ/qwfzjleDb/a4oCibIdg4UuwzmgsNXHrMk/O4z5xPxqKd0Q8E1WKBGE1oHo3ruw4TFgBMdfSVZ4/hxk0j2LZuKHTMUVS9qZsgxF1MZeCGJq6wZqPFznXRiiJOV5pO/Td/aRkOJevufGXyZj72dmI0uS77sde0YIJpNkNAQhJBNd1y7hX/nGwX342hqYYwmqSuM36Nuy0OmQQt240Zxe44wuIrhJDEFZx5Fn8nKKrZ0Dwafk3zSsaT+yaipnk3UUU16xgWDllRWWfjEjKvXztbw7ZR93lePZRDrWWi0TKdOmcjQg+iTStLSy5ps29ougTfQRXULCzq+szrLRNKhsUUHLWYz0A0k7rOdMthGNzQyMq87z4+i/1TVdwXwWaAaKquGey7oooSJmY0LRMlNYu1Q0wM0dQtUAoUQ1xngBunCXWdRZTP4dBN6iwMaanORBBC7PN674HFe/nY98rZxUskzt3EaOIYTdxC6zRu67LcfRK02nKdhSvGhgrJetJ0ozor5jJB15nuZzTZUOlwkNEE2Y+4AeWIivmx83rn4OoBdz3hrtcVQqXzTStLOHy2tiCMNSlC7wgh5E8JIW8R/r+ZEFJemGGdP+C7ME5rm/YOpqGbzmQKy+h3XGcxeTTsOzijyYWWofnKM8dRymXxths2RI7Zdf2EZbbH7IjtOR1vaAwUcwrWDOUx3zQwbe++wlxngKs8i47RxDc+c11npDvXWcvAUF5SDkdSxsdf+Tqa0XQeowljNLwkfxyjieoQ2g3OVLTAItwyTI/rrBHRCC4qkF8uJmM04rPSLopqFoZFPcVSHcm4k0eTCd0QBhhNLqhia+imk/jJEdei3C+gWT3kricyRvP9V61FpWnguy9PxfzihUPUHflNAN8n/P81AP+9t8M5/8B3Yf4ErYYg410dkmjpigFynv+L4B0v+QO4ZkjOjgBg17EZ3L51dSCm4IfbITI6xhAGh9HEuH14J8G1QwUAwLFzjM5HGZpBO/Ab5TqLqwzQMtyEzW4YDaUU1WaQ0QDy6gr+ytc8RiOLS3SrOuM5jR7XWUhmvR/5HjCaMxUNd39qHH89fsjzOmeX/J5HNYKLiq+UC2qs6oxSygxNm71oOPgGQdwX+CtKRzGaWssbz5OJATjLFz0GUTEa02K/SXS1ievJbF1HhsCzGbrzirXYOFLEl54+5jkXpRSTC9jGW0TUHakCKAn/J/afPgTwXRjXxfMFRMyAHy6qULMklNEMR4gB/B0veXWAM5J2AdWmgZUD8gVaRNQOSiwKGfr5DFPUJHGdFdWsYxyPnqsBQKjqDGDlRojvwRGRzHXmTdgMEz7EQTMsGBaVGxqJAWuH0fAFqKK1H5SvaYaToCcuejyBM25H3wvX2V98bz8qmhHIF+NGn8t+Y/NoQjY5Sdo5J01YDYPbkdWdL34xQBSjYX1mosUATd1EwbfRikpE5p8vSQzNuaqG2QYrqClWU8hmCH761k147OBZvHa25rz+R99+FXf8yfbELRfSRNSMPAjgxwkhbyCEbLZfG7FdaJF/FmDcSwZcEeZnNHyRBZhPf9VAPpD/wifxYF5BNkOkYgB/ZvLqiHpnVc3wZNeHIWoHlYTRAMw4xhmahk35HUNjByhLIdWbAXYtygVVWjoHSK46cxM2OxcDcCYiY4gyFx7fZDjVm213UZTrDGhfEFDTDKfcv0zeHMdo+PvdtiTmOHSmivufOQ4g6P515c1B19kXnzqKj3xrr/P/SNdZAjFAWIfOpODPqxi/98e9whgNpdTOuXLntqzSQENYFziczrGS++GKB9w5uEpwxc/UdUe1KuIdt2yCkiG4/xnGap4/Pou/23GYdehsox18Woi6I58EsBXAE2BuM4C5zl6L+XO4V4NdiuC6fX8msOg6A5hfNcho3IWpoGQCChXAXZBEeTMQLENjWZQZGkl5fT94PxW56yye0QDM0MTm0bQMO2HVx2gidpzXbhzG6zaPhL6fSHXmSdiMZ0Bh4G4eOaMJGrCmsylwG58B7mIhQnSTJinnI6IqGBqxDE3ShTbtGM0nH9jHGvMN5aWZ9TlhIya6zh58aRJff+6kI2xpRcRXWDvn6OuUNI8oDE6/GdF1JlGdyRgN73o64HOd+Y25LGE5Sgwgq42WV7IYKig4W21hrq47MV4Ra8sF3Hv1KL668zgqTR2//bUXwLdFSeoUpo3QVYlSej8h5DUAbwXrpPkuAC8AeH5hhpYMhJAMmAH8ZQBbAJwB8C8Afp9SWov4aCpwXGeOGMB1nYkTatVAMP9FE4xIMZeNZDR88R+xy9D4jRZXNoW5nEREUXVeVy0Ow6V4RlNvmVg9mMfKgRwIAY6c5TGa8DH++l2XRZ4zp2Q8wXAZWMJm+66zidkGBguKkyBYFRmNr1GiInOd+Xr5DDgxGlnCpvtaO8mTut04yzE0epDRLKTqbNfRaTzw0iQ+cO8VePDlSWkeCu9eySqcu797ar6JqmZgutbCqsF8aGUAgM3rpm5FutfC+tkkhcNohPniz03ijIZS6omzOFWZA/LmoBjAnwwd5WEQayaKWDOYx5mqhpl6C6PlgvT3vPMNF+Pbeyfxs3//DPZNVfBrY1vxV+OHEqn30kbkqmTXNnsKAAgh7wLwDUrpxxZgXO3gzwD8BoBvAPg0gKvs/7+OEPL9lNKeavw0w8JwUQ0wmnrL9EgOVw/mcWCqEvgswHzqeSWoUPEcYz88mQyR7u74YpaE0bjVh+WZ7UMJzjFcVGPrKfHdm2oXAjyWwHUWhySVAVqmKAZI7jr7ub9/GrdtXYU/ePt1ALyGxv9oymqoua4z7mZhZUhkjKbRMu2xUUc5lATcnbeCu87MoKGJU53xHX8arrM/+varWDOUx3vefAl2HDgTWFh10/KII8RrMWVX9D46XbcNTbg02S1DYyA/GGJoeOC+Q9WZVAxgsDSFrO3KzSsZWJSLHFxDI+sz4zesgOtOFhElzuGGxl8fkCdtztb10Jy527euwpZVJew5PosfvXEDfux1G/FX44cWhdG0c0d+AcA3ezWQTkAIuQbA+wD8K6X0xymlf0cp/QCADwC4C8B9vR4D1+3zicAnVlM3PX7V1Xa9MzH/hS9MuWwmNOPYr2QC3JbHInihxTjFGRCd2c6brMVhJEGMhhXQZONZM5R3Fu4o1Vkc1AT9ZcSETRnzCMPUvIaj51zjyRd1WUfQnKSGWtPnuiGEoKTKC2s2dQvrhtlOtJ1KyvwarhoIKhU5+42VN9t9fbplNE3dxLNHZvDTt25GKaewzHq/vFkQZhSFa9HUTWf+HJ+uu4qxUEPjtpAIg78uWbtwxACWVwwgjslt5+z9nbI+MwVJSRsx7YHDidFEuM787ubVQ6z1iFi52Y9MhuBX7tyKS1YP4CM/fI0jOlrShoZS+o+U0hd6OZgO8NNgSrg/973+dwDqYE3ZegouPXbEAC3OaAwUhZ3VmsE8WnalVQ6x/wbT3MsXfsBbe2yooAQCo/wBbIfRyIpNahJqL0MyMYCbV8DjNEC06ywOSfqr66Ybo0mSdwMwGWlVMzAl9IepRogBpKozSZ24Ul6RMpqmYWJ9mRVAbEcMwJkrX1xkjCZxjCbC0EzMNvChf30BNT382vHv5o3JZMFv3XCN/kDebZsgXuej5+pOb6Ywt61b7yx8zokegk5QlDCalq9aAb+2gYZmzqbEG6Pxt2hvRKrOgtfaVZ155+DqwTwm55qotUypGIDjvls3Y/tvjWHlQM7tVLqUDQ0HIeQ+Qsh3CSHnCCEGIWSaEPIQIaTn7EGC1wOwADwjvmg3Ynvefr+n4KozJ2HTnuwNX/HIYclNFl0FBSX4kPLzA96HZyivShhN+zEauessWWb1sN0DPmyxopSiLsSpvIamG0ZDYl0+uiHKm5O5zviOlLtzAEF1JjHeSpYEmJWsl89ALiuP0egWykXWIKsd15nDaAYjVGeJxQByqXGlqeMX/+FZ3P/McRyaDZcj+w1bUVZmX+j4Wcy5XTbFfI6j5+ouGwthI0naOSeNUYVBFqPR7IRTjnwoowkaBCdu6ym2awbYSS5CDCArWQOwmC+fCyMDckbjR0HNIqdklrahIQxfBPAlAHcDKIMF3ocA3APgS4SQf+7JKMOxAcBZSqlMr3cSwGpCSLK70CF4JjKfCKIYQNzZ8l1x1Vc+hE/cMDFAYtdZG4wmF9GpMkkeDeBWMwhjNbzcTElwnQFs4Vc7dG0A7DokitEovMNmMjEAv55zDd255lXJLpVDpmaTdSf1xyXEY4u5LFYO5DpynTliAF+tM0Lcvj1hiErYNEwL77t/N16dZPHEekTcOBgo9yqyKKUedllSs468ecpWTa4cyOHYdC2WjSSp4OxPBWgXBTsXzq86E8fkbCh9jKYucbM67nRf512/0YjMowkRA/DqAAAiGY0fSTwRvUA7PoxfBvBOALsA/DaARyilJiEkC+BOAH8E4KcIIY9SSv8m/aFKUQIQJgpvCsd4toyEkPcCeC8AjI6OYnx8vKMvr1arqDUJzkxO4OnHzwIAXt5/EN81jkI3KaZOHsP4+CQA4NBZNmF2PPUsplawSXP0uAZqmBgfH0d1romzDRoYy3OTbAK/uPs5TB9kE7I+p+HMnOk5ducJNnn27t6JyWL0g3Z4jo3luef3wJrwToFas4UzpyYwPn4u8nefnDgIAHjokcexYTD4ffMttrifPHoY4+PHMT/FxpfLBH9jOzg9qaHeNCLP0TIsTJw4jvHxKUwcb6FlWti+fXtk/bbjFfch//eHHsGaUgZ7D7Jps/PJx1Cv1TzfOT/bxKzm/S0vHmO/cdezT+FQ3l48Gg2c1KqB8c5V65g9q0GxLBw8fgrj48k6bT5rz4fDrzAv9r4DhzBOWQ7LwddaUAnwyCOPRJ5jXmP3Zu8r+zBac7MRKKX4p1daGD9m4O2XqfjmQR0z1WbotZ6qsWt2cP8+jNcO49xpDZWGe29022V08tgRjI9PoF5pYqbJrtnjr7FrdcmgiQOnZvHIjscBAK8dOoBx7Ujgu8412Hc9u2cvBqb3Scez+zS7Nnv37Eb1SPuspmEz1EpDc37DyVNNGJrl/P/AFPuOx596GifK7nfw5+/F53ZissTu/RH7te07Hsdq+5msNXWcmZzA+PhZ57Nn6uy3vbD3ZQzPHvCM6YUj7BzPPfMU9ufc+Ts55RrcYwdewfjM/kS/UbVaOHhsAuPj057Xq9XgHE0T7RiaXwRwBMAdlFKn2QGl1ATwMCHkTgB7AbwbwEIZmjqAtSHvFYRjPKCUfhasNTVuueUWOjY21tGXj4+Pw0QDW7dcjHvuvhK5730b6y7ajFtv3wo8+CCuvuIyjN1xKQBg+NgMPrXzCVx+1XUYu5IN+WsTz2FYn8fY2Bi+NvEcZifYv0VMP3cCeH4P3nT7G3DxqgEAwCOVl7Dn7AnPsYcfew3Y+zLuvfPNTiO1MKydmAee3IErr74WY9eu87ynP/ifuOzSizE2dmXk7/6+y6/G37zwLLZddyNuvnhl4Jjj03Xg4e24/porMXbLJszvmcD9r+5GuVQI/MZ2sKP6Mp6eOh56DtOioA/8Jy6/9BKMjV2OF80DwKH9eNMdd0YyqZ1HpoHHnwQAXHL1jbhly0o8XnsZxaPHcPddd2F8fNzznV8+vgvNszWMjd3hvHZwx2Hg5Vdw951vdnbgnzv8DOYaOsbG3uj5Pvrog9iyaQNwroaqZgTeD8PpnceB51/A3W+6DR97ahwbN12MsbFtAIDx+ZdQmDwRe33nmzqw/UFcfMlWjL35Uuf1pw+fw/e+8xR+6c2X4INv2YZv/t4DMLO50PPtn6oAOx7FDdddg7HrN+DJ+ivYMXHEOb7S1IEHH8S2y9lz8PVTuzF/cg5jY2PYUX0ZxcPHcNcNW7Hrwf248oabgUceww3XXI2x1wU7w1Y1Ax985DtYv/lSjN25VTqe2gungOeew+3fdyuuGA2vXh4G3bSA734bRHF/8xeOPIum0sTY2JsBAGT/GWD3M7j2hptw88UrnM8efeIIsPcl3HPnmxy2Ob9nAti7Gzfe/HpctnYIpkWhP/Cf2LZ1C8bGrnA+OzXfBB79HrZevg1jb/Dmu7+0/SDw6j7ce/cdHqY8dHQGn9n9BADgzttuwbUbhxP9xvWvPIGcksHY2Pd5XvfP77TRDse8GkzeLO2oY7/+TTB58UJhAsw9lpe8txHMrdazDk+UUk+f87yaQaNlOu4zkSJzyXDF7zpzKv3K5c3+hE12LhUVzdtKuRqhkPIjzHVmmBZMi0b2ouHgRfzCaLhf/89L53QTnwHiEzb5b3JiNAqXjka7z0RXJI/TVDVT6jYDeLFOuetMvH4DuWxIPxrmohwp5dpK2BSDzrlsxluCxjCRS3DveMzB7/qbtAP09926GXklg1w2g6ih+SsRFHzBbx7c5q41psDjsbAmRst5bLY3TwdPVwGEu70GcllkSJzqLDrOEwc1m4GSIQF5s/jsudJwuerM3yYAcJ9hWeVm/r1AeIwmmyGB38SfJwDShM0wDCdIfO0F2rkjFPG1zha6FtqzYL/hVs8gCCkAuBHAzl5+OXfTch8uNxaO9t0To7GTAJt+Q+M+pLLqzf6ETcBV+YjxnqpmoKhmoSR4yJSMfGInaePMwcUNYYtk3Ulg88ZoulGcAW6tM1mbBMBdPMWimuLrYRAfPq6IqmlGaE5RThqjsQK9fFiM8Dcc3AAAIABJREFUxntfeV+eoprFypLaVsKmKLnOq16jm7SNsWNofAZbrCxACEG5qKAeISXnv98N9nuD3/z8jrw558qbmaEpYLPd494xNCFzj43Hu0geOlOVj79D1RlgZ/P7Ejal8uZAjIZVZRaPlZWlAoKbrah0A17Kyu/25WIQwFu5OQ7lgrK05c0AXgGrfVaUvWm//nYAL6cxsIT4CpgB/E3f678EFpv5Ui+/nM+1nPAgNXRTWghv0DEOgupMUHjJSooD8mxnhx0JD12lmaz8DBC+y0/SxpljJEaT3/CxOtfQdMdocrx8Tkh/dd0XoFYjyu2I8DCaCjM0rBeNfLyyigOaYQZ6+Qzks4F+NO6CmMVIKYf5ph47PmecmoFcliX4MkbjVZ0lMTQZu0y9X73nn2vlgop6lLzZLwbwSX91nyEayDNDQynF1LyG0XIBF9uGZr+dzBylGBPrnT2w9xTu+fQjeHVyPjCeTlVnAMt90X1iAFHFF8ZoqnaLAPHe881hQ8gdYq/LGU1YHo1fPAAwRltUs1CzJLSRoAzDRRVzHbSm6BbtGJrPAdgM4FFCyD2EEAUACCFZQshdALYDuNg+bkFAKX0RwF+CGcB/JYS8hxDyaQB/CuARAD1VwfFgp6Mc8zEaUS/PSoMzdwxHU1CdFewaSv6dukx1xv3/4uJY1eR9U2QISxBz5LkJHtRyLKPxGttyQfGUi+8UUW4G9jr1HBeVoyBCVHOdcVxnbg/4wDgUmessqNgTd/EcbqkapjqjNHkSXU0wfjklE5A3x0mbOXISF6S70WDnGCqqiVRneR+j4Rstf52wUk5x2l5MzjexbriAkZKKoYKCA1PRrjPAW8H5y88yAYQoR+9WdQbYjEZM2NTljCagOmsFW34H6h/q3s0XhzNHJexRVhuNY/VQDiOlXKTIxY/hInO7y7rr9hLt3JG/BXA/gJsBPAigQQiZAlN3fRfMffXVBVSccfwmgN8CcA2Y0bkPwGcAvK3X5Wf4zsfr/rIc4yDq5TMZgsGc4nWdCYzG1ecHd5k8qZNjSGZomnpiRpMLWawdo5bA9ZDNkEgaXvexOkII1pULztg7hRri9uEIxGgkv/WLTx4JlM+pNHVkMwQXryo5jKamBRcPDiUjlzf7d9MDOQUtw/J8f1Nwh3L/etIGaDUhbpRXMoGimkkXWb+R4p9n53U3B4kYTZaX3PHKed332dzlrqSJ2SZahoW1Q3kQwq75EbvgahJGMzXfxKP7zwDwdRv1dcPsBEU1603Y9LWGDo/RBA2C40qMcZ1l7RI3Ua4zGVYP5tuSNgNsg0ipN1a8EGinMgCllL4TTOL8MIA5ACvtvx8G8E5K6YInbVJKTUrppyml2yileUrpRkrpByil1fhPdwcnRiMmrLXM0PpEgwXF4zprCeVe3MBhcPfrf3C460xMvKpGLIp+KCGGRpMID6KwYiAXmmzIA+BiTObP77sR77/3CunxSRFXedgfN/D7v6uagd/71kv42q4Tns9Vmuz6rSsXnF1yTdLGWRyHrB+Nn9HIWhiLFXm5QilpdQDxPueUbKCoZmJGIzM0uumJMZWLaqIYDc9ZKvjmsMx1BgCvnWWPJi/Bc/HKAfANdtQmh7Vz1vGN3Sed4/0xz2yGJIpThqGQ8xoaTTc91zQqj8Y/V/zPNL/vMtd0WGJxQzdCGc1brl6HH7hmnfS9MCxWdYC2I7OU0vvBmM2yB48TOMqxXBZnKpo0RgOwpM1AwqbAhgBGr8Ui+U3DDDx8roLNG6PhgdU4hLVydsrcJwymjhTV0J24bPd20+YV0mPbQVSlW/a6dxftd51VhcRMEZUmC/yPlgt47CDLcYhqu6BmSSBG05SU7+GLT71lOAIKl9FknVIySTttiiwr72sC1zKTJduyz2YDxprPR+6KYTGa8HP4VWf+hbXlY5e87t1rdhXvdXbV4c2r3HkbpRhjjMbA13adwGVrB3HwdNX3PJkdK87c71Bwcj681plbgsbHaDQJo/ExvDDVGRDezqLeMkPdt786Jpd5R0Gsd7ap7U93ju7uyjKHX3XGq7Xy7Gf/ojPoy+gXS9C4HTqDLMPPMKSus4S9aABADVOdSeS5URgu5TAbJgYI8Ud3C757DmvnzP3cfHHjWfIuo2Hj9TMIZmhUrBnKo9I00GiZkSxRzWZgWNTj625KXFduC2P3XjWdGE3GqcLcDqMZcBhNBi2fvDkpGw1znYlzNlZ1JimhDwhiAJ8h4q5kzmh4eXtxgxS1ySkXVUzON3HwdBU/f9vFAGSVNrpb0jYMFzHd9KvOZK4z77Xzt3EGxM0jO9bdfIU10ksuBugUi1VYs29ouoDfdcaLCjZCJlSA0QhGhC/uAddZFKPpUAyQyRAoEp+wW8AzOaOZC3Odhej/u0WcXJnX7/Ln0biGhr3vZ2KVpu4wGgCYmGugqVvhYgDOlITCpE092MvHLT1keo4D2D1f0WaMpupjNH7VWdLr7VesOeP3iU4MKzgnOTSfa8wv59Uc1xoXA3BDw+Ixa8tMiXixaGhiYjTsmAx+9HUbkVcyHgOeVHUXhQ0jRcxq1InB+FVnip1r478mdV9tQz5OwNveHQhnNDIxgL+vVbeQ1VxcCPQNTRfg7hh/0mU9ZEIN5pVgHo3q3Q0Ge4wHGY1THM+WN1NKUW1D3gzwhEO5vDnprnikpIYyGu5KaEcRkwRhOSAcLR+j8bva+PX3x5YqTWaoR+3F77UzbDEMlzcH3Y+yFgvcKIgLIr/HhRzrPJlXMm25zviY8t2ozpSgq8bfeEzsASMDZyz5LHcd23JeHqPxMxr7Wrx2poYVJdX5LtF1Fqc6A4AfvHYdygWmVgtLgO4UG0bYRmNqToNlK+T8YypIOmeK94UjYzd7c2M08gKZgFzFCESrzjpBuc9ozj9IxQC6WxnAzwxERuOfxP5AKodmmFKGUS4oToFBzbBgWDQ0cC2DtMx9J4ymoUulkrIGT2kgqu2t+DqvfsBdZ0aM66yqGR5Gc9h274QlbMqaVWl6sPL1oIR9agKjIYRgRSmX2HUmqs5yPkbTvupMwp49jMYWnYRkkvuFF3kfK+fv532MZmKu6ekKuX646BjuqPFz4cRP3HQRABb/qmlyV3Sn2DjC0gRPzjbc8fueB9bO2R+jMaQusaKQiB3lTg6L0TRaJopqd0nOIvqus/MQ3NDwhZknbIZl8w4WXEbjl5LycwQNjXzxGCqoTsImX8SSus4AuU84LKEsDMMllgMi2/GyFgHpPSAcScUAftcZf4j5WGdrMteZitEhtgAeOs0ZTbShERcHmRhAVrVbjNEATL2XpIIzpRS1lusi9efCtMNo/G43wHbleipQRLtZ3Mx/W77sk/MGKgMI14YrzgAm7920ohSrGHvL1evwNz97M958+WoAEg+BnkKMxjY0E7ON0EZqfkZjWawlhmyuiD16Gi32GZnrTPY8UkpZX6tcesv0QC6LbIb0Dc35BN2nOiuoWVh28p1sNz+UV1BtGXY3QW9ymfuQBiWnMneA2CrAadDVruss0Iq4vTwEruGfbQR3442WEar/7wZt59FkOPNgv5XvgCua4RxLKXVUZ+UiSyzlcYRwQxPiOvMzGonrzG/QV5TURD1pWFY9hDyabKoJm027sgEHd1WF9ZhvGRYyxJXL+1VWQXmzey25QefYtLIUO++KuSx+8Np1zgZuQKLi7DYmyA0gMzRchemdx35G0zTs+yJ55kVGU9cN5JSM0xZahCpp0KcZFizafdkmEYSQRWkV0Dc0XYDPNX8wdLrWkrKCwYICStmC4e+cGS4GkEtWywKjcXrR5JMnb6nZjCeQzb6rPUbDkw1l1QFkcs80ENUkCnCbVrmMJphHw8EftqbOXI+DBQWEEIyW847rLEp1BnjVb0wM4L1XAxJG0/AbmoEcphMYGn9/HOY6E1Vn3ebRtMloTCsyx8TPaMT5MDrsNTTXbRx25M5JMRQwNMlVd2EoqFmUcwQTc43QjqV5H6Phje1KIYzGKaoZkXwpy6MR863SBCusuUQTNsNACLmSEPJ+QsgvE0KS1aq+QCBTnQHAdL0lXWSdwpqaEchiDhMDhAU4RUbD82mSJmwC8hbHbTMabmgkC1FdN6UPXreI6q8OBAPQfheXGDzmLIJfPy4bHx0q4GyVvRenOjM8qjMrIA3PKxmoWSJ1nRU9jCZ+h+lvLS2KASilgSz2KLQjBgiN0fgYBFcZNpwYjbd6M8vRYcdy0QXH++65DN/49WStEjj8MZpWCvJmAFhVJDg52/QUGRXhZzROoVMZoxFqGEYpyGRVyf3VNdLCYhTWTHxXCCG/Twg5RQhZKbz2/QB2A/gUgL8C8BwhZFX6w1ya8LvOuC91uhZiaITAsOZLjvSXFOcIa63Ms6QBl9GEBa5lYHLK4HcpbWRWDxdZcFbm9mm0jEDL2jQQVj6Hw3Gd2UzG7zoTffpcUswNNg9+i4HqODEAV7lxd6ifDRJCArEEf/26lSVWYSGu/pS/L71oLPyB9zj4FWuATAwQ3dWSMSifYRVUVv6ETkKIY7j97CWvZJ1AdVKwShvyBOhusKpAmOssZONVUL3xLV40VeZm9bjOIhiNrNJElEqtG5SXuOvshwC8SikVW7N9Aqx68kcA/DWASwD89/SGt7ThVG9O6DobEtwofjFAPkIMID1XQQ3GaNpgELJ+KmHfFYaods5pyzI52o7RRLjOuKS44jPUa4XddnyMxk5ONCksKlfs+XfeTYOVNeH160ZKOVg0nDlw+O8zW5xY0mhY4DoMSRI2C2oGWRI+Ll0i/RX7KvkTOgF30Rxt000mg7zSRvdzzjE0ThzVz1KzjnIQcBMxZey3YPeoAuRiEQ5ZjCYqwbMbDBfVJZ1HswWsVQAAgBCyEazA5l9RSv+AUvrfwGqevT3VES5h6CZbcHhwj0+iSlNen8hpFSAyGp9bIVjrLJzR1FsmDNPqSAzAYjTBEirt7AijetLUU85o5mg3RuP03rFcMYB/3DzWxV2b4iIY1SZAHEczZFFi5/XmezRbpkc0wGW7cX1peCxgUBADsN9sufGEhK6jXDYbFAP47j8hBCU1WnXmjwlx5SXArg0vGMlRStnQNHXLkZiHPSvtYmUxg3rLxOkKq3nn/41+RsOfv5JkrhTUrDM3ojZfshhNWBHObrHUDc0KACKbeSMYm/kP4bVdYK0ElgUMi3oWFpEWS3tI5NyeNH5aTghBQQn2pGmG+J15PKGqGc6OvO0YTcyONg5qNoPBvBJiaMILUnYDN+YSV+vMl7BpuPLmi1YwCSvPXfG7HtfavXPULAndIftrqGk+ybKIoOvM8syPqArO9ZZb0p3nAIltAvh3+91Uccj5Kj8DctdTSSGRqjOxyRsAzxxumcH3i2oWSoZg1UDyZl1hGHAUffIs/k6xqsDGzJWHsoRN8Tmta+GMhhfaBViMJmzzJcuj6VUZJ+46C2se2Au0c1fOgLVH5rgLgA7gaeG1XJvnPK+hW95JKE4IWZKVWDqmKZFO8lppHFGtlcvCuaqaYS+KyS+9vJ9K+zvC4aIakDdbFmWLaQ9jNKGuM19uh8x1tq5cgJIhjojB7zrju+0ow+1vh+22WJAzWbH5WdMXy3EqOPsYzVxdxxv/6GH8zaOH7LF7GY1jaExXxdiu6kxcbJjqzDv+kkoSq84AVu2gIajO/IZvIK9g7VDe0/aiU3BXNBdzpOY6K7KxHbENjf83+nOQ3BiNXAzgJGy2wl1nsjyaRg8ZjWHRQJ+kXqKdVeV5AD9CCLmWEHIZgJ8C8BiltCEcswXAqRTHt6QRMDQeRiPf2QJy1Rn/vCgGiGpNyxnNXENn5WfySlvlXpSMzHUWXGjiMFIKduwLq16dBvyGww/dZO2UubuGu84MhxWwUj0jJbfFAY9BDDmuM8ZoohiZvx12VNMtmRhA3DysHy4iQ4CnXzvn+dwXnzqCmbqOb+2eABAUA/DvahmC6yzhQut8VsglkmXWl5TkqjMAKKoZZwcvM0TrhgvYunYw0Rjj4Gc0qanOCuwch8/Ke+QEGY23bbn/2EZC1Zk/r81xnaVYGQBYnOoA7dyVPwEwDGAPgH32vz/N3ySEFACMAdiZ4viWNHSL+hiJ+29ZAM/JqWi6YoCC7/Mio4nK1PczmnbiM4BcdRZW7iYKsnpnvfItA0FXmB8tk0LNuqXuHTm0fTwvs7+ipGKm5mU0g44YIJ7RBGI0evB+cviD1g1fJ841Q3m8/XUb8YUnj2Jqvmmfz8TnHz+CvJLBvqkKjpytoaYZIMS9rmIlYVngPQp+ZuiKGYKMJqzWmSxGI8YkdIkh+uOfuB5/+c6bEo0xDmJ79LC6ZJ1gKMeubZjrLK94PQ81Pt/DYjR259xGVB6NxMPQK9XZcIxsvRdIfFcopTsAvA3ANwF8A8BPUkq/LRxyO4Aj9nvLAn5G4zcafuSUDPJKxladBXfAedWrZgnT8QNiqwDdrujbnjQ0J5nYsiZrcRgpBpuf1VvBpmdpIZshICSa0YiLGyHEE2it2MVHxfpirGVz1mFBQ3Y/9ihG47rO7BhNRLKr39DI1Efv//4rYFGKzzx8AADw1V0ncK7Wwh/+2HUAgAdfnmT3OecyV5HR8O9vx3XGPyuOXxqjCdn5amZQ3lwUSq7IGM1gXnFk091i0F7Yq5opyLu7X5QJIdg4UsSZEDHAaLmApm45rs6aZkAJqVTODYtmWG3HaMIaKHYLh9EkrBieBtpaVSilD1BKf4JS+pOU0m/43nuYUvo6SunX0h3i0oVuhsdownbzQ7b2X2ZEimpG7joLSdgEbEbTTN4igEPJyIpqhvuQwzBcCmrye8lomOEIBrI5dEkAWsmw3jFsQbYwlFcYExNUZ2KLaUII1pbzkYYmjNHIDPVAnikETYsLB4LXedPKEu57/WZ8+ZnjOHymir979DBu3DSCH79pI65eX8Z3XppCtekVWOQkjKadopqA6zoL29SUVBIub5YwFpHRMLFA70K2TgJ003VFpyEGANyaZ0Dwmly6ZgCAW3i13mJ1zmSu66LNXHn7kNA8GkmMpt6jygDlwhJ2nRFCPkcI+ZGYY95GCPlc98M6P6Bb1FuCIyReI4Lvbp0YTSLXWfA2uSXc9c5dZ5LKAEmbnnGMFNmCLQaVe7UT48iF9O4AbEbjWxjULEHLsDwxDpHRcJYj4lfv3IqfuTVcQOlvhx3l5uSbAh40buryskLvu/syKFmCd33+WRybruNX7twKQgh+4Jp1eO7YDI6eq3sCzjm7PH/LsAK9YeKQDzAa+aamqLDxakYwcBxWQp8Xj5TdizTBr0UtxEPQDXi7ACAo8Lh0NYsxHbJbSdRsRiwDfwYqTQOGRSNK0ASfx4YtzpHVRusGSz1G8y4AN8YccwOAn+94NOcZmOvMnThKNuO2tQ2ZeLyCs+zBKPqCjEkYzTyP0bTJaGSuM1nb6DiMlJiCpSYoWMIav6UFWRY1R8uggV0064ZpeRIeOaOhlDotAkTcd+tm/OC14f3Y3VI43HUWLW8GXBl1IyRxb225gJ+/bQuOTddx6ZoBvOXqUQDAD1w7CkqBZ45Me+4zv1ctw3J39G3Im/lnAVE1F3SdAfIK3TJ5sziHtR4zGi7eqIR4CLqByGj81/SiFaytAY/h1FpGaLklsSwVEL4uqNkMTIs6rBdgLuheeAXcGM3C1TtLexbkASycZm6R4Y/RAELLgAhGwx+MDHH7pbDPehmN24gseJvUbAYFNYNKU5fuyOMg7UfTCaOR9LyvOTGa3jAaWXIbhz9Gw45nDIgbmiFbddYymd983m7j3A74dxg+RiPbFLjqKMM5Nmx+/MqdW3HpmgF88N5tjgR42+iQ0+7Y4zrLcteZG6NIKuZwP8uTHeWbmpLKxiCL08jFAG4JmoVlNEEPQTfghoYQBF2x2Qw2ryzh8Jmq/f3yFgGAa2j48xFqaCRqSlnXzjQwVFBAyNJlNABL0JSCEJIHcAeAya5GdB5BN2lgB8gnUtgiO5hXHdVZXvH2rMmrGeeBB1iyJns9LN6j2qozvcMYjc911gmjkdDwXun/OWQFCDlkcQGu6HEZjeppoczbOLc7BkCmOpMwGh5PEwxNWCxsxUAOD39wDG+9fr3zGnOfjdpjD8ZoRHkzd6fFIRijCWM07G/Z7lcW7C+qWRgWhW5XK0i7lbcIxd5shYlrugFvgMYqdgRdV5euGcRh23VWb0W4zoSyVOL//ZDV8Gv0qLpGJsPq7y1kdYDIu0IIOcz/2C+9X3xN+HMUwAyANwP4914PeqnA7zoD3B1M2ITiYgBZSfmA6yyC0fBzzdRbaOpW266zsITNdsUAnNGI1QF6VaOJIycodHYdncF1H/2OIwvWTcvZHXKods4Qd10N5LMeJlbpREyRDVGdhZSgAQRG00G+x1uuWec5F+DNheledSZ3PQ3EMRqfYRObn8kMUdpwYp49cp2FGcpL1wzg6Lk6TIvaLTHk84dfD25oovJoAG9/o17VCwSw4D1p4p6uDFwWQwEQ+48fOoAXAXwPwB+kNrolDpnrjBuYsJ3IQD7riAH8i1LAdRbTWrlcUHFqji2w7brO/CoXlrAXbNwVB7dVgOs6q/dI/88hxmiePHQWlaaBl0/NY7RcsMueyFxnlsMohgqKs8DO1lnCa9uMJuOPcYQn14oxGtNWv7Xrorxp8wpcvKqETbYLDZCXoGmnerM4/rB4II/RyJRnMtcZZ98N3YRu0J4yGsBNhk1bdbbe7pcT5k3YunoQLdPCiZk6ai0jtCYev888RhNVVBPwMZoIF2u3WOh6Z5FPF6V0C/83IcQC8GeU0o/1elDnC3SLSlU3QIQYwHGdhTMaSikIIZF+f4AtmK+cmrfP277rx6KAaVFkMwQt0wKl7fu4RySFNXspbwa8Cp19U8xPfny6DoDLm73XVckSGJarOvu/7Z15mF1Vleh/6041plKVkQyEhIQEDEOEBAEZikG7bUR90k79bAm2Tbct+kD7dbdDN3nYT1tBkBZbBVsmFRUbefKUQYQClDEyx0AYMhEyD6QqlbpVdWv3H3ufe0+dOufce+5Y99b+fV99N9lnn3P2PsNeZ6299lrtTUmcOdddfWkODmUiz9HEYkIiJtl8NGH3yrk3vS4TT1QhHI8J91xy+qiB2zlXOlPMgs3c+g53+70fNeY7YkyqACf/jZ/pDGBgcEQL/QprNE5k7KBIy8XSnIwzrT0VKLizLs47D4TO0TgRQpy1ZkGaj3dhMWjT2bT20mPC+THeNBo3Z6IXZFoMQyNjB2bnRQsKGzGpOcFgZoTegWFfR4IR5biNxgkLQQNao3ESdEX9Ik+4wtzHY/GiTQ8dPnM0/YMZUvFYxTyOHHdlgHXbegG3oBnrQqoXw+VMZ+3NCYxCkt0v6vVzjusIvNf3HmT6pCZfV1S36Sw7l1PEAOz9Gs5pNMXFOoPiNRqn30HOMAPDmYrP0cBo5xq/9pTC7M6WUQtt3Rw+Xbs4v7brQOgcjXPPdvflmaPxiUquvc5afeuXSkdzkleNM0M1iBIZ4EGl1MZKNqaeUEqNWbAJuS/VZp9YZ5AbdHYdGAyc38mmfi1Ao8kdtzivqULWgYTRnIzTkoyPig5wcHC4YmYzyCX8Ghweyb4sm0ZpNKMH+5THdNaajNNpkrY5+0XVCEELa2egfnHbfo48ZJJvPXfooWKvsx/uOZpBHy/GMHLOAMYVOWA+MBXXx/SaWZw5Mj/3ZtBf41rjKe8aEC/tRqMZzGNmLoalsydnnQK8dLUmmdyiB+sw7zDneuzNujf7t89vjiYsCGepjGeNBgARWQ6ciE4b4HcVlFLqy6U2bLwzPKJQBM/RBD14zoC250CamZNG5+RwHqr0UAZaknnnaEYJmqK9pkaHuS/mi7DTk4q4kpOYoNvelx5m/a4DDI8oYgKb9ujYrn5eZ06SNyf4aCwmpGJCWyrO5r2ORhM9LIozzzWcGeHlHX2sPGW+f71EjFQiRt/gcHYOrhwDiDtemWPGKjSwauAcjedZExE6WpJj1tEEpSVocc/R+Lial5v25gR9O8cmEiwHl793KUGR9EWEw6e3sWbLm0Bw3iKvM0BLntTgozSakCCcpTK5NVnVWGcFj04i0gHcjjahhT3NCmh4QRP0YOfzOnMEwu6+wezaCO++B10L3iDY88UdMyrqF3nCkyEyLFZXPnSqAJegCYnpVA6cAf6l7dpstnz+FNa+sd9omWPnBZLxGAcGM9mAmg6drSk2GwHVUaTpbDij2LD7AIPDI4EaDej4aaM1mtIH4FhMx3FLD4+QHspEGtSDFmwGBXD1Dkq5OaHR9ZuyWnmm4iFowG+Opnzny9f2w6e1c+ezb2Tb4YfjDODkGgqODOAsAPauo6mcRuNEfCincA4iytt1BXAW8DBwA7AZqN7S0nGGY2rwW7CWCgkb4bjR9g+OvcEtHtOZkzEw6CvVrdGUax1IsRqNO0Bff3rYN2R6uXDW0azb1ks8JnQvmc4T6/ewr3+IocxYTyctEEbGhOrpakvyxzeMM0Uxgsa4iK/dqgXekYd0BNZtz7q1B0d5LoaUuRaDmWgpHtxx0ty/fve/w8dDKcj5oMUjaCrt3jypKUFvBbzOCuHw6W1ZwRD0vMdiQioRy5nO8q2jMdfV8U6s1Aebe251xqTxJWjeCzwFnKmU8l8tV2VE5G/Qi0RPAI4AYkqpyhqFDUEv5okLpoSm5HUPaEETqW6NJmxAmlSCRpPymM5KmTvobEllAwxC5dI4O2j3ZsVL23tZMK2NRWZidtOeft85GieSQG96dFDKrtZU1vusGNNZMqbnil7ctp94TFg4oy2wblvKfHmX0XQGWoNwgoVG0miymUpHCxq/Y3Q0J8cs2HTmdvxSOYN+BoY9sQArQVtTgvTwSNalvhpf5w4Lp+fud5jm0ZKM551DSiZGv4/9FY6u4Wjw+w8OMWNS6Wm18xFldJoM3DJehIzh88BU4GlnumR+AAAgAElEQVSgDZhbrRMH2bTfu2wO7102x28XYLRA8A42za6vQX2O8IyXjhbjzlFSKF6Nxi8/TqF452gODmWyWSMrQVaj2d7L0bMnM2+qNkFu3tvv697seIf1DYyOoOAs2oRSvM5GeHFrLwunt4UOcu3N+su7nHM0oAVD2nh4RdFGvflonGfNT3vuaElkF8Q6BAkmZyB1NKBKm86c92m3+birpkazYFougVtYpO+WZJw3Dw7RkowHWie87+PBbGDaylgGTjtiOnd86u3M7aqMV5uXKHflZWBmpRpSJN3AZKXU6eiEbFWjWL/9cI1mtKDRGS/DBI3+Co+aXRNyczRjgioW8aJONsnPnAjOlXYGSJkcKZv29LN45iQONS/Lpj39oc4AB9KZUYLeCUMDRQqahDCcUby4rTfUbAbGO2pwODRUTTE0JWM5jSbCvXPmd7KCJiQXkdZo/N2bvV5ljmnI0YDKOWfiR865ZpBkXMoe6TiMw6a24rx2YYLGuddh74R3jiaXXbMy79GUthTLDu2smFeblyhPwbeB80Qk+HO9yiilNtRKwyrWS2uSyw3Zu+/klpyjAJj8MCGCzBkco4ZPAVdQyBFv9OHiTGeDwyPZQbQ/PVyx8DOgv/5608MoBUsOaaetKcHUNj2xP5QZa65xnAe8czTOYtNUPFaUySUZj7H7wCBb9h3kyFnBjgCQW8HuCPRyrfjWGs1IZI3GvS+E5yKa1JwYs2AzKLaac4yqaTQu55pqms1A93Vul3Z/DlpH49Rz//rhXW5Q6UXP1SbKU3AXcC/wexG5UESOFZF5fn8Vauu4otiMfs3JGM5Hl3fydsG0dtpScZ7ZvA8wX5khX75OuO+iJrID19EU5wwAuTA0lXTLhNHmmsUz9QB/6JRWNgfO0WjvsN6BoTFeZ1CcNgN6jsaJzHBUPo3GcQYowbvPj1Qilg2qGdVs5OwL4c9aR3OSg2Zy3yHIGcARds4ajWrM0YA2nVVae/LDyU0TlCYAcvNWYfOWY0xnQ8VFkBivRHnDNpCLd/b9kHoq4nGrjohcBFwEMHPmTHp6eiIfY80u/SD88YVnGd4SUdjEoX8Ytm7ZTE/P9lHb5rUrHvrjZno6d7F150EGMwS2r29QayMj6f7IfXhpj27/k394mgMb4jy3WQ8MTz35OOubw1/Yvr6+Ued7fZv+2r3voUeZ3Sb0HRxi17Y36OnZFalNhbJ1ixZoiRhseOFJNonQNDTAS7tGGB5RbNm8iZ6eXBDxbVvT9KeHGRiG3du20NOzU5e/odudUEN5r5+3zwB9vQezGsGe9S/Qsy34uu3ZPsj+/iFeWLsOgNWPP5oNWFkK6f6DbEv3kc5AUzz4WfFDZYbZ+PoWenp2sXnrAJn0yNg+9vWxffcGAO6+/0E6UrrNz+zQ1+6FZ5+mf6MnYkEcXnt9KwCvvfwSPf2vUSle2auf4617eolJtP4H4Xevg0ildbrnp594jJdT/vdzoE+70GcGgt/T3Qf1c/T8mrV0vflKdnx5ac3zsLXywiZKn4shikC4mZA0AcUiIp3AJRF2+Xel1J5SzqmUug64DmD58uWqu7s78jEya7fD6tWctOIEjp3bGWnfrsfup3/fQZYsXEB39xGjtq1Ov8R3HnyVE085lda1T9CVjNHdfZLvcYYzI3D/XcyeMZXu7hMjtWHSxr3wxCMsPeZYzlg8nfW/Xw9r/siZp586apLcj56eHtzXLPXqLr79zOMcuvhobn1iE8Oqn3eddDTdx82O1KZCeWpoHax/mSWHdHDWmacB8GT6RZ7oeRWAIxYeTnf3omz9R/vXct+m11DA0iUL6T5jIQDqpR1c99yTzOjqoLv71NBzevsMcOP6J/jj7p10NCd4/5+eGTpP9lzmZe7esI5D5s6HF9dxdvfpZdFqvrvuUUZGQIaGmTGpme7uFQXv2/HEA0yd3kV39zJ+uPFJ+mMDdHefNqpOT08PJ8xZxA/XPssxx5/Igmna02rgha3w1FOcdOIK3jJ7tDbX9tC9pNo6gN0cc/TSij0HALO29cLjD3FgWIf2L+Zd9uJ3r4Nonrcb+d16/uycE7L5g7zcsuFJ1u7ZwfSpnXR3n+xbZ0fvADz4Ww5ftJjukw5jcM02WP0HTjlxOcfMnVxsVwomSp+LoWBBo5RaWaE2dAKXRaj/Q6AkQVMOSlmJ7Jhv/PY9/rBOMiOK515/k/TwSNbf3Y9EPEZrKl7SHI3jt7+3fwiR4kKxOOFc/v62Z9nVN8hl572F8yo4uKTiuYRgDvOmtGZXcXs9oRJxyW5r97g3e8ui4Jg7jpzVkdcZI2fiSSNSvknyVEJ7NBUTV8xtOhvI4wwAo1MFhMVWa0nGs84D1YgMACY3VJXnaABOOnwqJx0+NbROvkXcMHaOZiKbziqCUmoD4ZEGxiWlrER2Xg4/m/hbD+0C4KlNe3XOmjzHn9PZkg1pHgVvRr+dvQNMbWsiUcTAkJ2j6R/imx9axvveWll/EWdwW+JaiX9o19jw+Q7uCWk/r7Oi52iMwDsqJCKAg/MxsKsvTXMi2M01Kk2JnNdZ1Bw3hToDdGRT/+YETVhaguaUS9BUOtaZy+kkav+rhXNdw73OGtsZoKg3TESOBI4C2pVSt5S3SfVBOiT/SD6cwc7Po6yrLcWCaW08tXFf3gWbALdedFJRD2PCyadiHuzt+9PMmNQU+TgAMzuaOf/4uZx33Cy6l8wo6hhRcF7KxW5B4wrn47eOxsHfGSD6Yk33cZfkcQSA3MfFrt7BsgZ+TCVy62iK0mhc66iCNDtHELvjneWCavoImkQ86zlZaMbPYnHHGKuFM0AhOIE0wzQa5zo640qjCZpId0ZElonIamAN8HPgRte2M0SkX0TOK28Txyf54pCFEabRALx1XifPbC5Mo5nW3lSUK3Eu5722Ke3oHWBmR3GCJh4TvvHB46oiZECvAUjEhLfMyg3wsyY3Z9dQ+EUGcHB76E1qShCPSQkajWM6y6/RtLk1mjKuXWgq2evMid48Emh6csKruEPmD4WZzlJxegcc9+bKajROOmeoblSAKORLhgj6Os7pbGGNCYl0sMLJA6tNlKCai4EedMTma4DFwLtcVR5Cz538OVVK52yE2nHmv4tM2ZfM//cppa6t1LlzOdajPwiTsnM0/gPD8fO6uP2pLWW15Xvxms527E+zdFblJx3LwZ8dM4vj5nYysyNnMkzE9Yu6aU9/waazWEy47Ly3cMJhXUW1I6vRzMwvaNpdgmZae3EC3Y+mhDZ/pYeDBUXYvrn00mMT8Tk4be9Pj9Vo/ASNk1cpaHu5aW9KMjCUrmpUgCgUMkcDcOqiadz1wlYyI4r+wQzxmFR8jqtaRPmUuwxIAScopdaKyGW4BI1SSonIo0Dhbi+lcz5wgafMiRy9EaicoCkhCGWYMwBoQQOgVPnWW3hx24QzI4pdfemiNZpqk4zHmD9tbFyxeVNa2bSn3ycygL+gAfjYyfOLbsepi6ahlApdFe497/6BYeaUMexHNqhmMRpNPMYe1zqaoMXBrcY8dWAwl2Y8KE0AjB5QqyNo4uzqG7+ms3xZdx3efsQ0frp6M89veVNH1wgJWVNvRBE0ZwO3K6XWhtTZBLyjtCYVjvGEW1mt87lJD48gFJ5oyk3WdBbwYiw5ZBKtqbiJ8FwhjSY7R6PY3ZdmRMH0jsoH16skh07Rq7S9giYVYDorlXOPncW5x84qqK77vOWco2lKxvUcjU9a5XyMWrAZotE0JeIk45LVfiAnaPxMY+6Po2p8ked7n2pNIaYzgFMWau+137+yi4MVDkxbbaLcmU7g9QKOV7loiuOI9HCGZJyivjiyGk3Aix2PCceZtTnFmOYKwTGdDWdG2L5fLzor1hlgvOA4BPhFBnAo1pW5VEYFUy3jXEIqHis6xYPXGSBs/9ZUIjtBDZAOSbTWXGWNxplDGrdzNKnCTGfT2pt4y6wOHn55Z8Wja1SbKE/BDsw8SAhL0XlqGp7B4RGK/TDNZzoDvZ5G16mQRuMyne3o1ZF5Z9a5RuMkkhu7jkb/PxmXmn31umNhlfNL1d2fqH1zHAmUUsbxJLhdban4KGeAweERmgK0FfeAWulYZ5Dzihu/7s35g2o6nHrENJ7auI89B9IVi9xcC6LcmfvRQTWX+G0UkRVo89o95WjYeCc9PEKyyEixc7taiceEqe3Byp8zT1MpjcYx+Q1mFDt6G0OjeduCqZy5ZDpHzhrtbuyYztqKiHJdLhLxWHYALrd7s9+/C913cFiH7RlR4YKqrSmRzZEC5kMroL5bkFZFo8njXFNrWgoIqulw6qJpDGZGeHLDXlrGqeAshig9+So6o+ZDIvJJYDaAiCw1/78T6AWuLHsrxyHpEjSaty+ayiP/dBazJrcE1lk+fwqLZrSHpgcuBREdJl6bzrRGU05vqFowfVITN1x44phcOM5Xda3MZg5tIeunimWUoIm6jiaeS5oG4QNha1OCvnTOdDaUCV630+xqUzU0Gue+1rvXGcCK+VOyDh6VjIBebaKEoHlJRM4HbiXnzSXAc+Z3H/B+pdSmsrdyHOLM0RSDiOQ1U01uSXLfZ88o7gQF4iTu2ndwiKltqXH7opZKYpwImknNCXb1pcuqpbrNXZEjAyRipDMj2ayfYfu3N8VHuzeHeLk1uzSaamgZhZiia4mzMLirgGSALak4y+d38ciruxvKGSDSm6eUultEFqBdik9CZ7d8E3gMuKHUYJf1RHqoeNPZeMHJPLlj/wDT69xsFobjHFDswsxy4QyI5cpFA16NJtpxHdPZQEg4GYfWVII9Bw5m/x/m5ebW2Kqp0YxX09lxcyfz40+8jeUFrtd6+6JpPPLq7oZyBoj85iml9qEXbF5T/ubUD91LpjNlZG+tm1ESybgwmBlhR2+67h0BwnAGu0LWu1QSJ1zKeJmjcQbmvmw2zHBnAK97c5DpzPkSj8eqk/FyvM/RiAinLJpWcP3TjpjGFfe81FCCZnzemTrgL0+ez7mH17cnt04INsKOEuKc1QPjZY6m3WRXLXcIGr9/F4IjKJwAmGEC0OsMEJY62tHYqrWqPRfSqTEG5qWzJzO3q2VU/L56J/KbJyLT0AE15wK+0QiVUjeX2C5LFUiaNRg7+9LMqJOoAMUwXkxnzvnHi0bj1HfikoVqNE2JMe7NgaazZM6dvBpknQEaJFxLPCbc99kzGqY/EC3WWRNwFfBxghdlCjo5mhU0dUAyLmzfP0BmRE0I01mtNZqc6awyGk2xprPerOksRKNJJRgY0uGK4jFtcm0L8Ipy+peq0uR8vgXQ9UilQk/Viihv3pXAJ4G1wE+BLWh3Z0udkozH2LJPT/BOBNNZredoxp3pLOGYzoygCTWdOfHOhuloTmr35jyms2rNmTj5chppTqPRiPLmfRDtyrxCKTWUr7Jl/JOMx9j2pl5DM6OBNRpncWqtNZr2Cmg0bk+zogXNwfymM2dNR386Q0dzMtQZwOlftUxnx86ZzNfOP4a3R5hwt1SXKG9eG/AbK2Qah2RcGDbx3BtZo+lqSzGpKcHCGe01bUcu4V05g2qW4N4cH206C3cGyGk0ED5H43idVWtdViwmfGjFvKqcy1IcUQTNGqCwULWWusAdPr+R19G0NyV49rJ3UuuI6+3N5TedubWKYhZsQoHOAEajcVycQ50BEo5G0zhzJpbSiPIkXAn8D5MAzdIAOINUV2ty3K6qLhexmNQ8t4djOivniu+SQtBEcQZocgSNjiIQumDTpC5u1EgTluhECUFzm4jMAh4Wkf8AnkJHBfCr+1CZ2mepII4Nfcakxp2fGU8cP6+Ld7xlJkvKGL+uHF5n+wtybzamM6PRpMMWbCatRmMZTdTZ0S70XM2/5KnX2J/HDYIzEDTyGprxxIyOZq7/2PKyHjNV0oJN/ZpmnQHyLNiEwuZomqvsdWYZ/0RZR/N5dDrn3Wj35jew7s11TVbQWI2mbnG0kJiMnnMrhEimM8frbDCDUkqbzgLOl4zHSDRQvntL6UTRaC4CXgNOUEr5msws9UXWdGY1mrrFuYfFzIe4F2wGZct0aHWZzjIjCqXCz9mSjFvTmSVLlCfhEOCXVsg0Ds5AMLOBPc4aHREhlYgV5czh9jrLZ+bKeZ1lsumfwwRNUzJunQEsWaI8Ca8BnZVqiKX6OBkSG3mx5kSgKREralB39jkwGJ7GGXT8reZkjAODwwyatAJhprGjZk1iUY3XLVnGD1FMZ98BLhORQ5RS2yrVIEv1cPLpNPJizYlAUyJW1HyIWzgVEuizvSnBgbRL0IQIt1v+6m2R22NpXKIImjuBbuAREbkc+APB7s0TIstmvZM1nVmNpq5JxWNFBZQctdizAI2oNZWgfzCTTf1sJ/sthRJF0KxHR2YW4D9D6qmIx7XUCMd01shRASYCTcl4cRrNKEGTf46nNRWnLz1c0ByNxeImikC4GS1ELA3CucfMor0p0XAhyScaqXisqDUrsZiQjAtDGVWQRtRukp8VYjqzWNxEiQywsoLtsNSAo+dM5ug5k2vdDEuJNCWL8zoDLaSGMplsfLIwWpsSvHlwiKGMNZ1ZomGfFIulzjl0Sitzp7QUta+jlRSm0cTpL9AZwGJxU9RciogciU7n3K6UuqW8TbJYLFG45kPLig4YqjWh/OtoIOcMYAWNJSqRnhQRWSYiq9EpA34O3OjadoaI9IvIeeVtosViCSMRjxGPFSdoshpNAaazNuMMkLbOAJaIFPykmPQAPcAS4BrgLk+Vh4A9wJ+Xq3EWi6Wy5ARN/qGgzesMYOdoLAUS5Um5DEgBJyqlPgs86d6olFLAo8CK8jUvGBGZIyKfF5EHRWSriBwQkTUicoWITK1GGyyWescRFoV4HrY1JRjKqGyqAKvRWAolypNyNnC7UmptSJ1NwOzSmlQw5wGr0NGkrwAuAR4xv0+LyCFVaofFUrdE0mhMwrY9Bwb1vlajsRRIFGeATuD1PHViaK2nGjwMHOYJh3O9iDwOXA/8vfmzWCwBRPE6azU5afb1D43a12LJR5QnZQewKE+dpcDm4ptTOEqpNQEx135qfo+uRjsslnqmKZIzgBY0e/uNRmMFjaVAojwp9wPnicgSv40isgJtXrunHA0rgbnmd3tNW2Gx1AG5OZpCnAG0MLIajSUqoufwC6ioBcxTQB96bmQZ8AngWOB0tLNAM3BsLYNqisjPgA8AZyul7g+ocxE6kRszZ8484Sc/+UlR5+rr66O9feKFQp+I/W7UPl/79ACrt2f46FEpzjksOWa7u9/r9mb4yuMDHDUlxto9I1z3jlZS8eLcqsczjXqvwyimz2eeeeYflFIF5SaPEoLmJRE5H7gVuNYUC/Cc+d0HvD+qkBGRTvQEfqH8u1JqT8CxPocWMtcFCRkApdR1wHUAy5cvV93d3RFOn6Onp4di961nJmK/G7XPd2x7mtXb3+Doo5bQfeK8Mdvd/Z7xxn54/GFUqg3o5Zwzu4kVuX5nPNOo9zqMSvc5UmQApdTdIrIAuAA4CZiKThXwGHBDkADIQydaGyqUH6LX64xCRD6B9j77FXBxEe2wWCYcUZwB3KazREwaUshYKkPkEDRKqX3oBZvXlKMBSqkNaI2oaETk42gN5V7gfKXUUBmaZrE0PJEiAzTlnAHs/IwlCnX/tIjIhWh35vuA9yml0jVuksVSN6TiWsAU5AxgvM7SwyNW0FgiEajRiMjpxR5UKfVQsftGQURWAt8HHgDeq5QaqMZ5LZZGIYpG05yMERMYUXaxpiUaYaazHopPdFbxTFoi8h50ps/96LUz53si2PYppe6odDsslnomSmQAEaEtlaA3PWw1GkskwgTN5YzvjJrHo01/nRgPMg8bAStoLJYQoizYBGhtimtBYzUaSwQCBY1SalUV2xEZ075VNW6GxVLXOAKjEK8zcBwC0lajsUTCPi0WywTGERiFpHKGnEOAFTSWKIQ+LSJyuoiMXcUVXP84EflY6c2yWCzVYE5nC62pOJ1tY6MC+NFqIjhb05klCvmelgeAle4CEflHEdkdUP99wA1laJfFYqkCZx81gye/eA4dzYUJmvYmq9FYopPvafFbSNmMnoC3WCx1johkF2IWQqsVNJYisE+LxWIpmPYmazqzRMc+LRaLpWBajTNA0mo0lgjYp8VisRSMk865yWo0lgjYp8VisRRMm52jsRRBIU/LeI4OYLFYqoh1BrAUQyHuJqtEZJW3UEQy5W+OxWIZz1hnAEsxFCJoouaKsRqQxdKgtNrIAJYiCBU0Sin7NFksliw2BI2lGOzTYrFYCsZJ55y0pjNLBOzTYrFYCsbxOiskf43F4mCfFovFUjDWvdlSDPZpsVgsBTN7cjOfOfsIzjlqZq2bYqkjCo+mZ7FYJjwiwmffsbjWzbDUGVajsVgsFktFsYLGYrFYLBXFChqLxWKxVBQraCwWi8VSUaygsVgsFktFsYLGYrFYLBXFChqLxWKxVBQraCwWi8VSUUSpiR3VX0R2AhuL3H0asKuMzakXJmK/J2KfYWL22/a5MA5TSk0vpOKEFzSlICKrlVLLa92OajMR+z0R+wwTs9+2z+XHms4sFovFUlGsoLFYLBZLRbGCpjSuq3UDasRE7PdE7DNMzH7bPpcZO0djsVgslopiNRqLxWKxVBQraCwWi8VSUaygsVgsFktFsYImIiISE5FLReRFERkQkc0i8g0Raat120pFRBaLyOUi8piI7BSRXhF5RkS+6Nc/EVkiIneIyF4ROSAiD4vIWbVoezkRkVYRWS8iSkSu9dneEP0WkSkicqWIvGKe5Z0i8oCInOap9zYRuc88D/tF5G4RWVardpeCiLSLyBdE5HnTn10i8oiIrBQR8dStu36LyOdF5DYRec08vxvy1C+4jyIyW0RuNs/JQRFZLSIfKKhhSin7F+EPuAZQwO3AXwNXAUPA/UCs1u0rsW//BvQCPwI+Dfwt8FPT32eBFlfdhcBuYDvweeDvgKfNtTin1n0p8Tpcaa6DAq71bGuIfgOHAeuBnea+fxy4FLgB+LCr3knAAPCq2X6p+XcvcEyt+xGxzzHgYSAD/AC4CLgEeNzc66/Ve79NP3YDvwH2ABtC6hbcR2AK8BrQB1xurl2POd+FedtV6wtTT3/AUmAE+C9P+afNBf+LWrexxP4tByb7lP+r6d/FrrKfmRd2mausHR3O5yWMR2O9/QHHA8PAZwMETUP02wy4m4FZeeo9AewH5rjK5piye2vdj4h9Ptnc06s95SkziO6r934Dh7v+/UIeQVNwH4Gvm2t3nqssbo6xG2gPa5c1nUXjI4AA3/SUXw/0Ax+teovKiFJqtVLqTZ9NPzW/RwMYM9p7gB6l1DOu/fuA7wOLgRUVbm7ZEZE4+l7ejdZYvdsbot8icjpwKvB1pdRWEUmKSKtPvUXo/tymlNrilJt/3wacIyKHVKvdZaDD/L7hLlRKDaLjfB2A+u63Uuq1QuoV0ce/AF5VSt3pqpsBvoXWdv4s7HxW0ERjBVqjecJdqJQaAJ6hDgaZIplrfreb32OBJuBRn7qPmd96vBaXAkcCFwdsb5R+O4PCJhG5EzgIHBCRdSLi/lhy+hLUXwFOqFwzy84TwD7gH0TkAyIyz8y3fRXdj1WmXqP124+C+ygis9CazmMBdd3H88UKmmjMBnYppdI+27YA00QkVeU2VRTzlf8vaHPSj03xbPO7xWcXp2xOhZtWVkRkAfB/gMuVUhsCqjVKv5eY3+vRX6MXAH8FDAK3iMiFZnuj9BcApdRetEa6B20C3Qi8CHwKOF8pdb2p2lD9DiBKH0u+HolITbO0An5CBvSkmlNnsDrNqQrfRE8afkEp9ZIpc8wsftdiwFOnXvgOenL8qpA6jdLvSea3FzjTmI4QkV+g5yq+IiI30Tj9ddOHnrv4JfAIWtB+CvixiLxXKfUbGrPfXqL0seTrYQVNNPqBGQHbml11GgIR+TLajHSdUuqrrk1OH5t8dqu762DMRe8ETldKDYVUbZR+HzS/tzpCBvQXv4j8EvgYWutplP4CICLHoIXLpUqp77rKb0ULn+tFZCEN1u8AovSx5OthTWfReANtHvO74HPQZrWG0GZEZBXwJbS76996NjuTqX7qslPmp2aPO8y9vAr4NbBNRBaZidLDTJXJpqyTxun36+Z3m8+2rea3i8bpr8Ol6IHxNnehUqof+BX6ns+n8frtR5Q+lnw9rKCJxpPoa3aiu1BEmoFlwOpaNKrciMhlwGXAzcAnlPFldPE8Wo0+2Wf3k8xvvVyLFmA6cC7wsuuvx2z/qPn/J2icfjvOLHN9tjllO9DPOwT3VwF/KG/TKoozKMZ9tiVcv43Wbz8K7qNSaitakJwUUBfyPfe19vuupz/gGMLX0Xy01m0sQx//xfTlZkIWoKK/CjPAca4yZz3JOupnPUkS+HOfv0+a63CX+f/iRuk3WlvZj9Zs2l3ls9BzGOtcZU+aurNdZbNN2X217kvEfl9t7uk/eModbXUPkGiUfpN/HU3BfQSuIHgdzV5gUlhbbJqAiIjIt9DzFr9Am1uOAj4D/B44Syk1UsPmlYSIfAq4FtgE/DNaqLrZrvRkqeOH/wR6RfzV6Ifzr9HC+Fyl1D3VanclEJH5aOeAbyulLnaVN0S/ReQi4HvAGvQq+RRauM4C3q2UutfUOwV4AC2UvmV2/zQwE3i7UurZKje9aETkMOAptKD9EfqdnYK+f/OBTyml/sPUrct+i8hfkjP7fhp9X79h/r9RKXWLq27BfRSRqWgNZyra1LwFva6wG231+M/QhtVa6tbbH1qKfw69CjxtLvhV5FkZWw9/wI3or5agvx5P/aOA/4dem9AP/I46CsOS51rMxycyQCP1G3g/eh3EAbQH2r1mgPHWOxn4LVrb6QXuAY6vdfuL7PNC4CYzuA6hPxQeAt7fCP0mFxYm7/sbtY9o0+Mt6MWtA2ih/aFC2mU1GtFl4cAAAArzSURBVIvFYrFUFOsMYLFYLJaKYgWNxWKxWCqKFTQWi8ViqShW0FgsFoulolhBY7FYLJaKYgWNxWKxWCqKFTQWi8ViqShW0FjGICLzRUSJyI21bkspiMhyEfmNiOwy/Xkm/16WciAiPSJiF+lZAJsmoO5wXl6llNS6LeMZEelAR+RtJrea2S9a8bjF3OsHlVLdtW6LpTBcoYtuUkqtrGljxhFW0FgalRPRuYO+qJT6Sq0bY7FMZKzpzNKoOOln3witZbFYKo4VNA2Ae07F/PsnZl5iQERWi8i7A/abJCJXicjrpu6LIvJZQp4LEWkVkc+LyDMickBE+kTkURH5iKdeSkSeNO16j89xbjHbvhShn2eLyN0isse0d52I/JuITPZeC3TgRIAbzHmUiKwMOfZHTB3fVM4i0iQie0Vkm4gkfPZ9wGwfEJG1IvKlgAR5iMiRIvIDEdkgImkR2SEiD4vIJ832la75jTNc7VcmIZ37WB8UkYdE5E0ROSgiz5v7M+bc5nwbRKTD3PcNIjLkPWbINXqbiPzcXINBEdksIt8Tkdn59x51nD8RkV+bZzQtIq+KyBUmuVxQm9tF5GpzzoPm+XufqZMQkS+IyMvm+r8qIhePPXNJ5281dTaZfV4RkX8UEXHVXYU2mwFc4LlvK00dEZELROQREdlp2rtZRO4RkQ9FuY51Ra2jjdq/yNFZlb5to8rmm/IH0AmrHkOHsL8JHWU1g84N796nCR3uXgHPAF9Dh43fi45MrIAbPft0oiO2OkmRvgV8G3jFlP2rp/7hwJvAbuBQV/mFpv5vCcl54znW36DTFvSiw9r/m+mnQoe673S1cRVwh9l2h/n/KmBZyPGb0dGYt2Nykni2f9Ac70pP+X+a8s3m399Ah5937kfCU/9cdMTnDHoO6avAd9AphtebOstMexWwwdX+VUC361hfMXV2mmNcgc5BotBRfJOec29Aa3irgdeA64ArgQsKuP4XAsPoSM+3Al9Hp8rImGPO89TvwfOcmnIn39Fu9PN5BTpisHMfO3zavAV4FB0x/VrT7l5z7rOB/0JHY77ebN9ujjcmsnAJ5/+duWbfQz/zW8w+l7nqdgPfJPdOue/bMs89e80c5yvoLLYvAD+v9fhSqb+aN8D+Rbxh4YJm1INvtv2JKf+1p/wLpvy/cA32wAJ0Aig/QXMj/omjmoG70YJgmWebM0A/jE6xcJQZrLYDhxTY58PQKRn2A0d6tv2HOf51nvKVpnxlhGv7PbPPu322/cpsO8bnHLcDLZ76q8y2/+Uqm4YWvIPAGT7nmOtzr3sC2nqy2b7JfR3R8653mm1f8OyzwZTfB7RFuC6LTZtfAeZ4tp2FHvB/4Snv8XlOzzTnfwTzYeBzLa8OaPOdQJOr/DRTvgedwKvTte1w096ny3j+X7vvMXr+b5/5S7rK5+Pz7ri270YLxVafbdMKvSf19lfzBti/iDcsXNBsAOI++2wEdnnKXjYDxEKf+qu8Lws64dEw8GRAu44z+3zdZ9t3zbYr0emQR4B3RujzF83+X/HZ5mSLPOgZiJyBY2WE85xi9rnNU36I6ftTnvKn0TlNOn2OFUd7uj3hKvucOf41Ee51T8C26832i3y2LTb39jVPuTNoHlfI+V37OZkpzw3Y/gtzfSa5ynp8ntNfmOMsDTjO08COgDb7PaevmW1n+Wx7wLQpXqbzL/Kpf5PZdrSrbL733fHssxttXmvy296of9brrLF4RimV8SnfjCs3uIhMAhYBm5VSr/rU7wEu85StQA+eY+YJDEnze5TPtkvQg/jnzP+/qkwGxwI53vze792glNorIk8DpwNHAkVnPlRKPSIi64DzRKRLKbXXbPqf6L7f6NQVkVa0cN0FXOIy1btJM/p6OPnV7yq2jS7Crsk6EXkdWCAinUqpfa7NA8BzEc/lPDtniMgKn+0z0NdnMSbPfMhxhoAPiMgHfLangOkiMlUptdtVvi/gOX0DrYH7nXOLadMh5t+lnP9NpdQrPvU3m98un21B/AidxXKNiNwGPAg8qpR6M8Ix6g4raBqLfQHlw4ye4Hcmz7cH1PdbbzLV/K4wf0G0ewuUUgMi8it0uuNhtG06Ck57twZsd8rHTOYWwU3A/wU+jJ73ALgAPUDd6qrXBQgwnbFCOQinfVtCaxVGIddknqnnfi52KPNpHQHn3v/vPPXG3Huf4yTIf73a0V/+DkGD8DBAwCA9bH6TrrJizx/2XoEWaIVyKfAq8HHgn8zfsIj8GvhcgECre6zX2cTEeTFnBmw/JGSfq5VSEvJ3pndHETkVPUjtQr/oP5AAFSBPe/3aBTrPvbteKdyCNu1dACAib0ULyF8rpXb6tOnpPNfD3U9nwJpThnYWe02iChn3MSbn6euDBRxnb77rpZTaWEQbC+1HLc+PUiqjlLpGKXUc+v07H23Sew9wd5CnYr1jBc0ERCnVi5nYFZGFPlW6fcqeQA/Ap0U5l4hMQWsCQ+iJ4x8B7wT+McJhng5ql3FJXYY2Ca2N0jY/lFKb0eaot4nIEozAIecu7dTrQ3spLTV9LITHzO+7Cqw/QvDXctg1WQTMRXuxBX2NR8Fpd6R7H3CcLhFZWuJxxvP5HdN1Xi1HKbVDKXW7UuqD6GduIXB0BdtWM6ygmbjcgL7/XxOR7HMgIguAz3grK6V2oIXEchH5Z/GsJTH7LjT7u7kRPehdqpR6HvhbtCPCl0XklALb+kO0oPq0GUTdfBnoAH6olEoXeLx83Gh+/wr4CNqM8v996l2Ftuv/IGANRpeIHO8qugntuPBJETndp/5cT9Fu4NCANv7A/H5JRKa7jhFHO13E0O7W5eBa9PW/WkQWezeKXjNViBC62vxe77f2RkTaROQkb3kZqcb596K1xnk+x28SvRZMPOVJwPlY6S/x/OMSO0czcfkG8D606v6UiNyDtud/CHgIrcp7uRg4Argc+EsR+R16nmc2etJ7BXpgXg8gIpcA5wG3K6W+C1oTEJEPo9dF3Coiy1yT7r4opTaYY33btPVn6LUjZ6AneF8kmoaUj9vRAuEStI3/W0qpIZ92/UBETgD+DnjVXMNN6EFjAdpB4Qa0cEUptUtE/gL4OfCAiNyFnpjvAI5FCxW3oP4t8GERuRM94T0MPKSUesg4Lnwd+AfgBRH5Odpt/F3or+LfodeIlIxS6kUR+ThauK0RkbuBdebazENrOjvRzhhhx/mtiPwTeu3Qy2ZeYj16TuQw9P38HfCn5Wh3Lc5vnu/HgdNE5Efo65QBfol+Nu4DNpg6G9FLA96Bfn9+qZQqWSsfl9Ta7c3+Rfsj3L35xoB9erz7mPIO9Ff5FrTp6UW0Z9jhQcdDf8FfjF6L8Cbas2oTelC8BJhq6p1gtm3A3/33M+Ycd0To+zuBe9FfjWm0+e/rAcdfSUT3Zs/+33euNXBCnrrvRms8O9DrN7ahTY3/imfdj6m/FLjZXPdBtLB+EI+rMtqb68dme8a0ZZWnzofRg2OvuYdr0O7gzT7n3QBsKOHZOwat7W00138PeqHh9/C4GAc9c2bbqcDP0F5jg2gh9Yx5FpcX2uY857jRXK/5FTz/KnOObk/5IvS6n91o86cyz2MS/WFwF/qdGTDnfgz9MZIq9t6M9z8xF8ZisVgslopg52gsFovFUlGsoLFYLBZLRbGCxmKxWCwVxQoai8VisVQUK2gsFovFUlGsoLFYLBZLRbGCxmKxWCwVxQoai8VisVQUK2gsFovFUlH+G8MK8CToTXqmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "n = 100\n",
    "a = [[1.0/(i + j + 1) for i in range(n)] for j in range(n)]\n",
    "a = jnp.array(a)\n",
    "rhs = jax.random.normal(jax.random.PRNGKey(-1), (n, )) #Right-hand side\n",
    "f = jnp.linalg.solve(a, rhs)\n",
    "\n",
    "#And check if everything is fine\n",
    "er = jnp.linalg.norm(a.dot(f) - rhs) / jnp.linalg.norm(rhs)\n",
    "cn = jnp.linalg.cond(a)\n",
    "print('Error:', er, 'Condition number:', cn)\n",
    "\n",
    "u, s, v = jnp.linalg.svd(a)\n",
    "rhs = jax.random.normal(jax.random.PRNGKey(1), (n, ))\n",
    "# rhs = jnp.ones((n,))\n",
    "plt.plot(u.T.dot(rhs))\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Index of vector elements\", fontsize=20)\n",
    "plt.ylabel(\"Elements of vector\", fontsize=20)\n",
    "plt.xticks(fontsize=18)\n",
    "_ = plt.yticks(fontsize=18)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "**Can you think about an explanation?**  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overdetermined linear systems\n",
    "\n",
    "- Important class of problems are **overdetermined linear systems**, when the number of equations is greater, than the number of unknowns. \n",
    "\n",
    "- The simplest example that you all know, is **linear fitting**, fitting a set of 2D points by a line.\n",
    "\n",
    "- Then, a typical way is to minimize the residual (**least squares**)\n",
    "\n",
    "$$\\Vert A x - b \\Vert_2 \\rightarrow \\min$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Overdetermined system and Gram matrix\n",
    "\n",
    "The optimality condition is $0\\equiv \\nabla \\left(\\|Ax-b\\|_2^2\\right)$, where $\\nabla$ denotes gradient. Therefore,\n",
    "\n",
    "$$ 0 \\equiv \\nabla \\left(\\|Ax-b\\|_2^2\\right) = 2(A^*A x - A^*b) = 0. $$\n",
    "\n",
    "Thus,\n",
    "\n",
    "$$ \\quad A^* A x = A^* b $$\n",
    "\n",
    "- The matrix $A^* A$ is called **Gram matrix** and the system is called **normal equation**. \n",
    "\n",
    "- This is not a good way to do it, since the condition number of $A^* A$ is a square of condition number of $A$ (check why)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Pseudoinverse\n",
    "\n",
    "- Matrix $A^* A$ can be singular in general case. \n",
    "- Therefore, we need to introduce the concept of pseudoinverse matrix  $A^{\\dagger}$ such that solution of the linear least squares problem can formally be written as\n",
    "\n",
    "$$x = A^{\\dagger} b.$$\n",
    "\n",
    "- The matrix $$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^*$$ is called **Moore-Penrose** pseudoinverse of the matrix $A$.\n",
    "\n",
    "- If matrix $A$ has full column rank, then $A^* A$ is non-singular and we get \n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^*. $$\n",
    "\n",
    "- If matrix $A$ is squared and non-singular we get standard inverse of $A$:\n",
    "\n",
    "$$A^{\\dagger} = \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = (A^* A)^{-1} A^* = A^{-1} A^{-*} A^* = A^{-1}$$\n",
    "\n",
    "* If $A$ has linearly dependent columns, then $A^\\dagger b$ gives solution that has minimal Euclidean norm "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Compute pseudoinverse via SVD\n",
    "Let $A = U \\Sigma V^*$ be the SVD of $A$. Then,\n",
    "\n",
    "$$A^{\\dagger} = V \\Sigma^{\\dagger} U^*,$$\n",
    "\n",
    "where $\\Sigma^{\\dagger}$ consists of inverses of non-zero singular values of $A$. Indeed,\n",
    "\n",
    "\\begin{align*}\n",
    "A^{\\dagger} &= \\lim_{\\alpha \\rightarrow 0}(\\alpha I + A^* A)^{-1} A^* = \\lim_{\\alpha \\rightarrow 0}( \\alpha VV^* + V \\Sigma^2 V^*)^{-1} V \\Sigma U^* \\\\ & = \\lim_{\\alpha \\rightarrow 0}( V(\\alpha I + \\Sigma^2) V^*)^{-1} V \\Sigma U^* = V \\lim_{\\alpha \\rightarrow 0}(\\alpha I + \\Sigma^2)^{-1} \\Sigma U^* = V \\Sigma^{\\dagger} U^*.\n",
    "\\end{align*}\n",
    "\n",
    "* One can check that $\\Sigma^{\\dagger}$ contains just the inversion of nonzero singular values.\n",
    "* If singular values are small, one can skip inverting them. This will result in a solution which is less sensitive to the noise in the right-hand side.\n",
    "* The condition number for the Euclidean norm is still just the ratio of the largest and the  smallest non-zero singular values."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A canonical way to solve linear least squares\n",
    "\n",
    "Is to use the $QR$ decomposition.\n",
    "\n",
    "- Any matrix can be factored into a product \n",
    "\n",
    "$$ A = Q R, $$\n",
    "\n",
    "where $Q$ is unitary, and $R$ is upper triangular (details in the next lectures).\n",
    "\n",
    "- Then, if $A$ has full column rank, then\n",
    "\n",
    "$$ x = A^{\\dagger}b = (A^*A)^{-1}A^*b = ((QR)^*(QR))^{-1}(QR)^*b = (R^*Q^*QR)^{-1}R^*Q^*b = R^{-1}Q^*b.  $$\n",
    "\n",
    "- Thus, finding optimal $x$ is equivalent to solving \n",
    "\n",
    "$$ Rx = Q^* b. $$\n",
    "\n",
    "- Since $R$ is upper triangular, the solving of this linear system costs $\\mathcal{O}(n^2)$. \n",
    "- Also it is more stable, than using the pseudo-inverse matrix directly. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Padding into a bigger system\n",
    "\n",
    "- Instead of solving $A^* A x = A^* b$, we introduce a new variable $r = Ax - b$ and then have\n",
    "\n",
    "$$A^* r = 0, \\quad r = Ax - b,$$\n",
    "\n",
    "or in the block form\n",
    "\n",
    "$$ \\begin{pmatrix} 0 & A^* \\\\ A & -I \\end{pmatrix} \\begin{pmatrix} x \\\\ r \\end{pmatrix} = \\begin{pmatrix} 0 \\\\ b \\end{pmatrix}, $$  \n",
    "\n",
    "the total size of the system is $(n + m)$ square, and the condition number is the same as for $A$ \n",
    "- How we define the condition number of a rectangular matrix?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example of LS\n",
    "Consider a two-dimensional example. Suppose we have a linear model \n",
    "\n",
    "$$y = ax + b$$\n",
    "\n",
    "and noisy data $(x_1, y_1), \\dots (x_n, y_n)$. Then the linear system on coefficients will look as follows\n",
    "\n",
    "$$\n",
    "\\begin{split}\n",
    "a x_1 &+ b &= y_1 \\\\\n",
    "&\\vdots \\\\\n",
    "a x_n &+ b &= y_n \\\\\n",
    "\\end{split}\n",
    "$$\n",
    "\n",
    "or in a matrix form\n",
    "\n",
    "$$\n",
    "\\begin{pmatrix}\n",
    "x_1 & 1 \\\\\n",
    "\\vdots & \\vdots \\\\\n",
    "x_n & 1 \\\\\n",
    "\\end{pmatrix}\n",
    "\\begin{pmatrix}\n",
    "a \\\\\n",
    "b\n",
    "\\end{pmatrix} =\n",
    "\\begin{pmatrix}\n",
    "y_1 \\\\\n",
    "\\vdots  \\\\\n",
    "y_n \\\\\n",
    "\\end{pmatrix},\n",
    "$$\n",
    "\n",
    "which represents overdetermined system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deVhV1frA8e8SGQ6KkIqYs5ZSmhXJNRRTHHHKyMxsULNuZmWZlhVNNnn1Xq20X6OleS3Luko4C06ImppjOaWZGgEqzoIys35/bECPgCIc2Gd4P8/jo+yzz94vC3hdrL3etZTWGiGEEI6nitkBCCGEKBtJ4EII4aAkgQshhIOSBC6EEA5KErgQQjioqpV5s9q1a+smTZqU6b3nz5+nWrVqtg3IgUl7XCRtYU3aw5oztMfWrVtPaK39Lz9eqQm8SZMmbNmypUzvjYuLIywszLYBOTBpj4ukLaxJe1hzhvZQSv1V3HEZQhFCCAclCVwIIRyUJHAhhHBQksCFEMJBSQIXQggHJQlcCCEclCRwIYRwUJU6D7y0zp49y4kTJ8jKyio85uvry969e02Myr5Ie1zkyG3h4eFB7dq18fX1NTsUUQGitycxa9k6+p6PYma1xxnTsxURQfVtdn27S+AZGRkcO3aMBg0aYLFYUEoBkJqaio+Pj8nR2Q9pj4sctS201qSnp5OYmIinpydeXl5mhyRsKHpbArui3+cb9T24QdS5DkRG5QHYLInb3RDK8ePH8ff3x9vbuzB5C+GMlFJ4e3tTu3Ztjh8/bnY4wpZS9nLjwvt4vcpMNufdRI/M/7BLNyM9O5dJMftsdhu77IHXrVvX7DCEqDQ+Pj6cPHnS7DCELeRkwtoPYO371MvzZFT208zPCwUudkaTz6Tb7HZ2l8BzcnKoWtXuwhKiwlStWpWcnByzwxDl9fcvsOBZOP47tB7I4P292Z3pUeS0en4Wm93yqkMoSqkZSqkUpdSuS47VVEotV0r9kf/3dTaLyLi+LS8nhF2T73cHl5kKS16C6T0gMw0engv3fckTPdticXezOtXi7sbY8ECb3bo0Y+AzgZ6XHXsFWKm1bg6szP9YCCFcy/5Y+CQEfpkGbYfDMxuheXfAeFA5oX9r6vtZUEB9PwsT+reu3FkoWut4pVSTyw7fA4Tl//u/QBzwss2iEkIIe3b+BCx7BXb+D/xvgsdjoWHbIqdFBNW3acK+XFkHmwO01kcAtNZHlFJ1SjpRKTUcGA4QEBBAXFzcFS/s6+tLampqkeO5ubnFHndV0h4XOUNbZGRkXPVno7TS0tJsdi1nYNP20JqAY2u48cBXuOWm81eTB0lodB/6zwvwp43ucQ0q/Gmh1noaMA0gODhYX21h9b179xY7p9dR5/pWFGmPi5yhLby8vAgKCrLJtZxhAwNbsll7nEmARaPhwApo8A/o9380rXMzTct/5TIr6zzwY0qp6wHy/06xXUiuZ9euXVStWpXly5ebHcpVRUdH4+HhwR9//GF2KEJUjrxc2PiZMdb91wbo9R94LAbq3Gx2ZGVO4AuAofn/HgrMt004rmnMmDGEhobSvXt3s0O5qoiICFq3bs3LL8sjD+ECju0xZpcsewUat4dnNsGdT0IVt6u/txJcdQhFKfU9xgPL2kqpRGAcMBH4USn1OJAA3F+RQTqzDRs2sHz5cqKjo80OpdRGjRrF0KFD2b17N61atTI7HCFsLycT1r5vFOV41YD+X0HrAWBnUz6v2gPXWj+otb5ea+2utW6gtZ6utT6pte6qtW6e//epygjWGX366afUqlWL3r17mx1KqfXv3x9vb28+//xzs0MRokTR25MInbiKR5edJ3TiKqK3J5XujQmb4PO7YM2/4Zb+8MxmuPV+u0veYIdroVSkgi9o01cWX9sXtILk5OQQHR1N9+7dcXd3t3otPT2dBg0a0KhRIzIzM61e++c//4mfnx9z5swp031HjBiBUork5OQir+3btw8PDw9GjRpV4vurV6/OXXfdxf/+978y3V+Iiha9PYnIqJ0k5ZetJ51JJzJq55V/5jNTYfGLMCMcsi/Aw/Og/zSoVquSor52LpPAL/2Cakr5Ba1gW7duJS0tjbZti84ftVgsvP322/z99998+umnhccjIyOZPn06kyZNYtCgQWW6b7t27QD45Zdfirw2evRoatSowVtvvXXVaxw7dozff/+9TDEIUZEmxewjPTvX6tgVF5LaHwOf3Ambv4I7R8DTG6F5t0qItHxcJoFf8xe0EuzZsweAG264odjXH330UVq1asWECRNIS0tjypQpTJw4kbfffpsnnniizPcNCQkBiibwxYsXs3TpUt555x2uu+7KqyMUxLx79+4yxyFERSlpwagix9OOw9zH4buB4FkDHl8OvSaCZ/VKiLL8XCaBl/oLWokKlhCtWbNmsa+7ubkxceJEjh8/TkREBGPGjOHZZ5/lzTfftDqvbdu2xMTElPq+gYGB1KxZ0yqBZ2dnM2bMGG655RaefPLJq16/Vi3j18qUFJlBKuxPSQtGFR7XGn6dA5/8A/YugM6vwZPx0PAflRhl+bnMsn/1/CyF42GXHzdLwSJGWusSz+nbty933HEHK1euZNCgQUydOrXIOcUNhVxNSEgI69evR2uNUoqpU6eyf/9+VqxYgZub9RSp4q5fELMsxCTs0djwQCKjdlr91l24kNTpv2DR8/DnKmh4J9z9EdS5ycRoy85leuBjwwMrfGWwa+Xv7w/AqVMlT+L58ccf2bFjB2CsG22rhBkSEsLZs2fZt28fKSkpvPvuu0RERNC1a9dSvb8g5oLPQQh7culCUpC/kNS9LYnImA+fhhhLv/aeDMOWOWzyBhdK4JWxMti1uuWWWwBKrGqMjY1l8ODB3HvvvQwaNIgZM2YU2ftx4cKF3Hrrrdd870sfZL766qtkZmby/vvvFzmvpOsfOHDA6nMQwt5EBNVn/StdmNmzGuuHBRCx9VGIiYQmHYyHlG2fgCqOnQJdZggFKn5lsGsVFBREjRo12LhxY5HXNm3aRP/+/QkNDWX27NkkJiYyb948IiMjrYp+tm3bVqY1NO68806qVKnC9OnTWbduHWPHjqVZs2ZFzivp+hs3biQgIIDAQPN+gxHiqnIyaXJoNsRHgZcv3DcdbrnPLud0l4Vj//fj4Nzc3Ojfvz+rVq2ymuu9d+9e+vTpQ4sWLYiOjsbT05MbbriBxx9/nPnz57N+/frCc4tLsE2aNLnqUIuPjw8tW7YkPj6eOnXq8NprrxV7XnHXT0tLY+3atdx/vxTgCjuWsBE+70CTv36E1vcbBTl2WE1ZHpLATfbUU09x+vRpFi1aBEBCQgI9evTA19eXpUuXUqNGjcJz33zzTSwWCy+99FLhse3btxebYOvVq3fVexfMP58wYUKJq/kVd/158+Zx4cKFIrNVhLALGedg8Qv5BTkZ/HrrOLj3c7suyCkrlxpCsUdt27YlPDycKVOmcN9999GoUSP+/vvvYs+9/vrruXDhAmAsoXrixAkSExO5/fbbC8/57bffOHnyJDNmzLjifbOzs4mLiyM4OJihQ4cWe05x1wf46KOPuPfee2X8W9iffctg8Rg4lwwhT0Pn1zi9YYvZUVUYSeB24P333+e2224jNjaWHj16lPp927Zto2nTpvj6+hYei4mJ4bbbbisxKReYPHkyhw4dYvbs2SUOtxR3/ejoaHbu3FnmMn4hKkTacVj2MuyaB3VawsBZ0CDY7KgqnCRwO9CqVasy7Upe3PDG2LFjGTt2bLHnnzp1ipiYGH777TcmTZrEmDFjCqsyS3v9iIgIsrKyrjlWISqE1vDr9xDzKmSdNwpyQp+HqkV3g3dGksAd2LWuyR0TE8NDDz1EnTp1GD16NBMnTrTp9YUAY92hSTH7SD6TTj0/C2PDAytm9tfpw7DweTi4GhqGQL+PwN+1ZkVJAnchDz74IA8++KDZYQgnVrBoXEEFZMGicYDtknheLmz6HFa9B6qKUZAT/LjDz+kuC0ngQgibudKicTZJ4Ed3wYJnIXkbtOgJfd4H3wblv66DkgQuhLCZCls0LjsD4ifB+ing5QcDZkCr/k41p7ssJIELIWymQhaN++tnWPAcnPwDbnsQwv8F3sWv4OlqXG/QSAhRYWy6aFzGOVg0Br7uBbmZ8EiUUZAjybuQ9MCFEDZTMM5d7lko+5YayTvtKIQ8A11eA49qFRCxY5MELoSwqXItGpeWAktfgt0/GQU5D3wLDdrYNkAnIglcCGE+rWHHd0ZBTvYF6PI6tB/lMgU5ZSUJXAhhrlOHjB1yDsZBo3bGDjn+LcyOyiFIAhdCmCM3xyjIWT0elJsxp7vNYy5ZkFNWksCFEJXv6C5YMBKSt0OLXvkFOfaz2YqjkAQuhKg82RkQ/x9YP1UKcmxAflcxSVxcHEopJk+ebHYoVzRz5kymTJlidhjCGRxeD5+Hwtr34dYHYORmp9rezAySwMUVSQIX5ZZxFhaNhpm9ITcLBv8EEZ9KQY4NyBCKEJdJTU0tcYs5cY1+X2xsb5Z2DNqNhM6vSkGODZWrB66UGq2U2q2U2qWU+l4p5WWrwMRFmZmZ/Otf/6JVq1Z4eXnh5+fHwIED2b59u9V5eXl5jB8/no4dO1K3bl08PDxo1KgRTz31FCdPnixy3VmzZtG2bVv8/PyoVq0azZo14+GHH+b48eOAsTnymjVr+Ouvv1BKFf6Ji4u7YryLFy+mU6dO1K5dG4vFQqNGjejfvz/79++3Om/dunWEhoZisVgICAhg5MiR7Nq1C6UUb731VuF5M2fOLPG+YWFhRbZ2i42N5YEHHqBZs2ZYLBb8/Pzo0aMHa9asKfb9TZo04eDBgwwYMICaNWta7UOqteazzz6jTZs2eHt74+PjQ+fOnVm9evU1t6dLSUuBH4fCnIfAUhMeXwHh4yV521iZe+BKqfrAc0BLrXW6UupHYBAw00axCYy9K3v27MnPP//M4MGDGTlyJGfPnuWLL74gNDSU+Ph4goONraOysrKYNGkS9913H/fccw/VqlVj8+bNTJ8+nXXr1rF161Y8PIzCiG+//ZahQ4dy11138c4772CxWEhISGDp0qWkpKTg7+/PlClTiIyM5MSJE3z44YeFMd18880lxrtmzRr69etH69atiYyMxM/Pj+TkZFasWMGBAwdo0cKY37tp0ya6deuGj48PL7/8Mn5+fsyZM4chQ4aUu81mzpzJqVOnGDJkCA0aNCApKYmvvvqKrl27snr1au666y6r89PS0ujUqROhoaGMHz+elJSUwtcGDx7M999/z4ABAxg2bBiZmZnMnj2b7t27ExUVRb9+/Urdni5Ba9gxG2Jeu1iQE/o8uLmbHZlTKu8QSlXAopTKBryB5PKHVDzP1ePg5L6KunzZ1G0Nva68q015ffzxx8TFxbFs2TLCw8MLjw8ePJh27drx4osvFvZMPT09OXLkCBbLxZXfRowYQfv27fnnP/9JdHQ0AwcOBCAqKgofHx9WrVpF1aoXvw3efffdwn9HREQwZcoU0tPTeeSRR0oV7/z588nLyyM2NpY6deoUHn/jjTeszhs9ejR5eXmsX7++MKk//fTTdOjQoZQtU7Ivv/ySatWse3ojRoygVatWTJgwoUgCP3nyJK+99hrvvfee1fGffvqJ2bNn88UXXzB8+PDC46NGjSIkJIRRo0Zx9913o5QqVXs6vVOHYOEoOLQGGrWHu6dKQU4FK3MC11onKaUmAwlAOhCrtY69/Dyl1HBgOEBAQMBVf/329fUlNTW1yHGPPE1O7rXvG1mR8rKzyCwm1tIo2F0+MzOz2M+3wKxZs2jRogWBgYEcPny48Hhubi5hYWF89913pKSkWCXt1NRUcnNzC/9u27YtAGvXrqVXr14AeHt7c+HCBebOnUvv3r1L3Ng4NzcXrfUVY7yUl5cxijZ79mweffRRq2RW4Pjx42zYsIGIiAiuv/56q2uPGDGCzZs3W7VLRkYGYLTZ5XGUFF/Bx2lpaWRmZgLQpk0bNm7caHVubq6x+cCTTz5Z5BozZ87Ex8eHbt26WbU9QI8ePZgwYQLbt2+nefPmpW7PkmRkZFz1Z6O00tLSbHat0lB5udRPWkjTQ7PRyo2DzUeQXC8cdidTgX26Uqvs9qhUWusy/QGuA1YB/oA7EA08cqX3tGnTRl/Nnj17ij1+7ty5q77XkaxevVoDetKkSVc8z2KxaOCKfxISEgrP/+GHH3Tbtm21u7t7kfOGDRtWeN7+/ft1s2bNNKBr1aql+/fvr7/88ssi7dypUyfduHHjUn9eJ06c0EFBQRrQPj4+ulevXnrq1Kk6JSWl8JwNGzZoQL/++utF3r9jxw4N6HHjxhUe+/rrrzWgV69eXeT8Tp066UaNGlkdO3DggH7ggQe0n59fkTZQShV5v7+/f7Gfy80333zVto+Pj9dal749S1LS931ZFNdOFSb5V60/76j1uBpafzdI6zOJlXfvUqrU9qggwBZdTE4tzxBKN+CQ1vo4gFIqCmgPfFuOa4rLaK1p3bo1H3zwgdXxCxcu4O3tDVA4vhoVFcUDDzxA27ZtmTp1Kg0bNsTLy4vc3Fx69uxJXl5e4fubN2/Onj17WLlyJStXrmTNmjU88cQTjBs3jvj4eG644YYyxVurVi02b97M2rVrWb58OfHx8YwePZpx48axZMkS2rVrV9ABKLaXWvDapa7Um83Jsf6tLC0tjY4dO3L+/Hmef/55WrdujY+PD1WqVGHChAmsWrWqyDUK2rG4WPz9/fnuu+9KvH/BA9SKak+7lZ0Oa/ILcrxrwv0zoWWEzOmuZOVJ4AlAiFLKG2MIpSuwxSZRiULNmzfn+PHjdOnShSqXrBFR3FS3b775Bi8vL1avXm2VlH7//fdir+3p6Unv3r3p3bs3AEuWLKFPnz588MEHfPLJJ8CVk2dJ3NzcCAsLIywsDIDffvuNNm3a8N5777F48eLCZLZ3794i7y3uWM2axnzhU6dOFXnt0KFDVsM0K1euJDk5mRkzZjBs2DCrc19//fVr+jyaN2/O/v37CQkJoXr16lc9vzTt6RQOrzN2yDn1J9z+MPR4T+Z0m6TM0wi11puAucA2YGf+tabZKC6Rb8iQIRw9erRID7zAsWPHCv/t5uaGUsqqp621LvJwDuDEiRNFjt1xxx2AdaKsXr06p0+fLrZnXJzirnvTTTdhsVgKr1unTh1CQkKYP3++1dTCrKwsq9kuBQoecq5YscLq+Pfff09ysvUYq5ubsRvM5fHGxsayadOmUn0OBYYMGUJeXh6RkZHFvn5p25e2PR1axlnjIeXMPpCXA4OjpSDHZOWahaK1HgeMs1EsLmnlypWFD+kuVbt2bUaMGMGoUaNYvnw5Y8eOZdWqVXTp0oUaNWrwxx9/sG7dusIeN8CAAQOYN28eXbp0YciQIWRnZxMdHV34wPRSPXr0wNfXl44dO9KwYUPOnDlTON968ODBheeFhISwaNEiRo4cSfv27XFzc6NLly5WM0wu9cQTT5CYmEiPHj1o3Lgx6enp/PDDD6SmplpNEfzggw8ICwsjNDSUZ555pnAa4eVDIgCBgYF069aNL774Aq01t99+Ozt27OCnn37ixhtvJCsrq/DcDh06ULduXV544QUOHz5MgwYN2LFjB9988w2tW7dm586dpf7aFEwd/Pjjj9m2bRt9+/aldu3aJCYmsmHDBg4cOMDBgwevqT0d1t5FRkHO+RQpyLEnxQ2MV9QfeYh5UcFDzJL+BAYGFp6bnZ2tp06dqoODg7W3t7f29vbWzZo10w899JCOiYmxuu60adP0zTffrD09PXXdunX1E088oU+ePKkBPXToUKvzunXrpgMCArS7u7uuW7eu7tWrl161apXV9dLS0vRjjz2m69Spo6tUqVLiw8QC8+bN03fffbeuX7++9vDw0LVr19YdO3bUc+fOLXLumjVrdLt27bSnp6f29/fXTz/9tN65c2eRh5haa33kyBE9YMAA7ePjo6tVq6Z79uyp9+zZU+xDzF9//VWHh4drPz8/Xb16dd2pUycdHx+vhw4dqo1v+YtK85B21qxZukOHDtrHx0d7enrqxo0b63vvvVfPmTPnmtuzJHb7EPPcUa1/GGw8pPy0vdaJW2137UrizA8xlS7lr8a2EBwcrLdsufIw+d69e4stFJHyZmvO2h6HDx+madOmjBs3zqoaszinL2Rx7GwGWbl5eLhVIcDXi+u8HXMHl5K+78siLi6u8PlDmWkN27+B2NeNFQTDXob2zzlkQY5N2sNkSqmtWuvgy4/LWijCIZ2+kEXS6XTy8jsgWbl5JJ1OB3DYJG43Th3ML8iJh8ahRkFO7eZmRyWKIQlcOKRjZzMKk3eBPK05djZDEnhZ5ebAxk9g9QSjp933Q7jjUdkhx45JAhcOKSs375qOi6s48puxQ86RXyGwD/SZDDXqmR2VuApJ4MKuNGnSpFRTFj3cqhSbrD3cpLd4TbLTYc2/Yf1H4F0L7v8vtLxHCnIchCRw4ZACfL2sxsABqihFgK+saFxqh9bCwueMMe+gR6D7uzKn28FIAhcOqWCc21lmoVSq9DOw/E3Y9l+4rgkMmQ/NwkwOSpSFXSZwrXWZSriFa7nO24PrvD0cfkplZU7lZe9CWPyiUZDT/lkIexU8il8LRtg/u0vgVatWJScnB3d3x5tvKkRZ5OTkFLvsrk2lHoUlLxoJPKA1PDQH6gVV7D1FhbO7BO7l5UVaWhrXXXed2aEIUSlSU1ML11G3Oa1h2yyIfQNyMqDrOKPn7YAFOaIou0vg/v7+JCQk4OnpicVikaEU4bS01qSnp3PixAkaNWpk+xuc/NMoyDm8Fhp3yC/IudH29xGmsbsE7uXlRUBAAEePHi3cSQWMHUsqrJfigKQ9LnLktvD09CQgIMCm8au8XFj3IcRNBDcPI3EHDZGCHCdkdwkcjG3VfH19rY7FxcURFCRjdgWkPS6StrjEkV+5Y9uLkHYQbuoLvSdDjevNjkpUELtM4EKIa5SdDnET4OeP8XCvAQNnGQU5wqlJAhfC0R2KN3bIOX0Iggaz2bsnHVr2NTsqUQlkUEwIR5V+GhY8C/+92/h4yAK452Ny3K++/ZtwDpLAhXBEe+bDJ3fC9tkQOgqe+pnoszcSOnEVjy47T+jEVURvTzI7SlHBZAhFCEdy7ohRkPP7IqjbGh76EerdTvT2JCKjdpKenQtA0pl0IqOM7eMiguqbGbGoQNIDF8IR5OXB1plGr/vACuj2NjyxGurdDsCkmH2FybtAenYuk2L2mRCsqCzSAxfC3l1akNPkLmNed60brE5JPpNe7FtLOi6cgyRwIexVbjb8/H9GQU5VL7j7I7hjSLFrddfzs5BUTLKu52epjEiFSWQIRQh7lLwdvuwMK9+GFj1g5C/QZmiJGy2MDQ/E4u5mdczi7sbY8MDKiFaYRHrgQtiTrAtGQc6Gj6FaHRj4DbTsd9W3FTyonBSzj6Qz6dT3szA2PFAeYDo5SeBC2IuDa4yx7tOH4I6h0P0dsPiV+u0RQfWJCKpPXFwcYWFhFRensBuSwIUwW/ppiH0dtn8LNZvB0IXQtKPZUQkHIAlcCLNobRTkLBkLF05C6PMQ9gq4y4NHUTqSwIUww7lkI3H/vgiuvw0emWv8LcQ1kAQuRGXKyzM2E17+JuRmGePcIc+Am/woimsn3zVCVJYTB4yHlH+tM8a4+04pUpAjxLUoVwJXSvkBXwG3ABp4TGu9wRaBCeE0crPh548g7t/g7gX9PoagR0qc0y1EaZW3Bz4VWKa1HqCU8gC8bRCTEM4jaZuxVvexncYGC73+Az51zY5KOIkyJ3ClVA2gI/AogNY6C8iyTVhCOLisC7B6PGz81CjIeWA23CybLAjbUlrrsr1RqduBacAe4DZgKzBKa33+svOGA8MBAgIC2syZM6dM90tLS6N6dVmovoC0x0X21hZ+p38lcN+nWDKOknx9OAebDanUTRbsrT3M5gzt0blz561a6+DLj5cngQcDG4FQrfUmpdRU4JzW+o2S3hMcHKy3bNlSpvtJdZk1aY+L7KYtLpyC2Ddgx7dQ8wbo9xE06VDpYdhNe9gJZ2gPpVSxCbw8Y+CJQKLWelP+x3OBV8pxPSEck9awJxqWvGQU5HQYA51ekoIcUeHKnMC11keVUn8rpQK11vuArhjDKUK4jnPJsPhF2LcYrr8dBkcZO+UIUQnKOwvlWWB2/gyUg8Cw8ockhAPIy4OtX8OKt4xpgt3fhZCnpSBHVKpyfbdprXcARcZlhHBqJ/4wpgYm/AxNO8HdU4xFqISoZNJdEKK0crNh/VRY8x+jIOeeT+D2h6UgR5hGErgQpZG0Nb8gZxe0jMgvyAkwOyrh4iSBC3ElWedh9b+MgpzqATDoO7ipj9lRCQFIAheiZH+uNhafOvMXtBkG3d8GL1+zoxKikCRwIS534ZSxQ86O2VDrRnh0CTQJNTsqIYqQBC5EAa1h90+w9CVjm7O7XoCOLxkPLIWwQ5LAhQA4mwSLX4D9S6FeEAz+SQpyhN2TBC5cW14ebJ0By9+CvBzoMR7uHCEFOcIhyHepcF3H98PC5yBhAzQLM3bIqdnU7KiEKDVJ4ML15GTBzwUFOd5wz6dw+0NSkCMcjiRw4VqStsL8ZyFlN7S61yjIqV7H7KiEKBNJ4MI1ZJ2HVeNh02f5BTnfw029zY5KiHKRBC6c34GVsOh5OJMAwY9Dt3FSkCOcgiRw4bwunIKY1+DX74yCnGFLoXF7s6MSwmYkgQvnozXsmgdLX4aMM3DXi9BxrBTkCKcjCVw4l7NJsHgM7F8G9e6AfvOh7i1mRyVEhZAELpxDXh5smQ4r3r5YkBPyFFRxMzsyISqMJHDh8LzPJ8LXveDvjdCsM/T9UApyhEuQBC4cV04WrJ9K8JaJ4FkdIj6D2x6UghzhMiSBC8eUuAUWPAspezjhH0qdoTOlIEe4HEngwrFkpsHq8bDxM/C5Hh6cw54jFupI8hYuSBK4cBwHVsDC0XA2Af7xT+g6DrxqwJE4syMTwhSSwIX9u3AKlkXCb3OgVnMYtgwatzM7KiFMJwlc2K/LC3I6jjWKcqQgRwhAEriwV2cTYdEY+NFwcXcAAA62SURBVCMG6reBfgsgoJXZUQlhVySBC/tSWJDzFug8CJ8Adz4pBTlCFEMSuLAfx/cZUwP/3gQ3dDF2yLmusdlRCWG3JIEL8+VkwboPYe1k8KgG934Btz4gBTlCXIUkcGGuvzcbve7je+GWAdBzIlT3NzsqIRxCuRO4UsoN2AIkaa37lj8k4RIy02DVe7Dpc6hRDx78AQJ7mh2VEA7FFj3wUcBeoIYNriVcwR8rjB1yziYaBTndxoGnj9lRCeFwqpTnzUqpBkAf4CvbhCOc2vmTEDUcZt9n7Ab/2DLoM1mStxBlpLTWZX+zUnOBCYAP8GJxQyhKqeHAcICAgIA2c+bMKdO90tLSqF69epljdTYO1R5aUyclnhsPfEXVnAskNLqPvxrfj67ibpPLO1RbVAJpD2vO0B6dO3feqrUOvvx4mYdQlFJ9gRSt9ValVFhJ52mtpwHTAIKDg3VYWImnXlFcXBxlfa8zcpj2OPO3sUPOH7FQPxj6/R9NAlrSxIa3cJi2qCTSHtacuT3KMwYeCvRTSvUGvIAaSqlvtdaP2CY04dDycmHzV8YOOWDMLmk7XApyhLChMidwrXUkEAmQ3wN/UZK3ACDld2NqYOIvcENXY4ccKcgRwuZkHriwnZxMoyAnPv/B5L3T4NaBUpAjRAWxSQLXWscBcba4lriy6O1JTIrZR9KZdOpvXMXY8EAiguqbHRb8/Ut+Qc7v0Pp+Y8ikWm2zoxLCqUkP3IFEb08iMmon6dm5ACSdSScyaieAeUk8MxVWvgu/TIMa9eGh/0GLHubEIoSLKdc8cFG5JsXsK0zeBdKzc5kUs8+cgP5YDp+2M5J32+HwzEZJ3kJUIumBO5DkM+nXdLzCnD9h7JCz80eoHQiPxUCjOys3BiGEJHBHUs/PQlIxybqen6VyAtAafvsRlr1iDJ10egXuGgNVPSvn/kIIKzKE4kDGhgdicbeeR21xd2NseGDF3/xMAsweAD8Nh1o3wIi10DlSkrcQJpIeuAMpeFBZOAvFz1Lxs1DycuGXL2HlO8bHvf5jLEAlBTlCmE4SuIOJCKpPRFD9yikPTtmbX5CzGW7sDn0/AL9GFXtPIUSpSQIXReVkwtoPYO37RkFO/y+Nud1SkCOEXZEELqwlbDJ63Sf2QeuB0HOCFOQIYackgQtDZqoxzv3Ll+DbAB6eC827mx2VEOIKJIEL2B8Li0bDuSS480no8rpssiCEA5AE7srOn4ClL8OuueB/EzweCw3bmh2VEKKUJIG7Iq3htx+MasrMVAiLhA6jZU63EA5GErirOZMAC5+HP1dCg7bQ7/+gzk1mRyWEKANJ4K4iL9dYdGrlu8Z0wF6T4B+PS0GOEA5MErgrOLbHmBqYtCW/IOdD8GtodlROo2CN9uQz6dSrjOpYIfJJAndmOZnG7jjrPgSvGtD/K2g9QApybMgu12gXLkMSuLNK2JhfkLMfbn0AwidAtVpmR+V0rrRGuyRwUdEkgTubjHNGQc7mr/ILcuZB825mR+W07GaNduGSJIE7k/0x+QU5yfkFOW+AZ3Wzo3Jqpq/RLlyarAfuDNKOw9zH4LuB4FkDHl8Ovf4tybsSmLpGu3B50gN3ZFrDr3MgJhIy0yDs1fyCHA+zI3MZl67RLrNQRGWTBO6gvNKPwbf94c9V0PBOuPsjKcgxScEa7UJUNkngjiYvFzZ9wT82vwVV3aH3ZAh+HKpU3miYzHsWwj5IAnckx3bnF+Rs5UzNNtQa8t9KL8iRec9C2A95iOkIsjNg1XvwRUc4fRjum87O1m+YUk15pXnPQojKJT1we/fXBlj4XH5BziAI/5dRkBMXZ0o4Mu9ZCPshCdxeZZyDlW/nF+Q0gkfmwY3mF+TIvGch7IcModijfcvg0xDYPB1CnoanN9hF8gaZ9yyEPSlzD1wp1RCYBdQF8oBpWuuptgrMJaUdh6Uvwe4oqNMSBs6CBsFmR2VF5j0LYT/KM4SSA7ygtd6mlPIBtiqllmut99goNtehNfz6PcS8ClnnofNrEPq83RbkyLxnIexDmRO41voIcCT/36lKqb1AfUAS+LU4fdjYIefgamgYAv0+An8ZjhBCXJ3SWpf/Iko1AeKBW7TW5y57bTgwHCAgIKDNnDlzynSPtLQ0qld3orU9dC4NEhfR9NBstFIcbDaU5Ho9QZXusYTTtUc5SFtYk/aw5gzt0blz561a6yLjqeVO4Eqp6sAaYLzWOupK5wYHB+stW7aU6T5xcXGEhYWV6b125+guoyAneRu06Al93jeWfr0GTtUe5SRtYU3aw5oztIdSqtgEXq5phEopd2AeMPtqyVtgFOTET4L1U8DLDwbMgFb9ZYccIUSZlGcWigKmA3u11h/YLiQn9dfPsOA5OPkH3PYQhI8H75pmRyWEcGDl6YGHAoOBnUqpHfnHXtVaLyl/WE4k4yyseAu2zAC/RvBIFNzY1eyohBBOoDyzUNYB8rv/lfy+BBa/AGlHod1I6PwqeFQzOyohhJOQUvqKkJaSX5DzE9RpBQ98Cw3amB2VEMLJSAK3Ja1hx3dGQU72BejyOrQfZbcFOUIIxyYJ3FZOHYJFz8PBOGjUztghx7+F2VEJIZyYJPDyys2BTZ/BqvFQpSr0+QDaDKvUHXKEEK5JEnh5HN2ZX5CzHVr0yi/IkTVChBCVQxJ4WWRnQPx/YP1UsFwHA76GVvdKQY4QolJJAr9Wh9cbO+ScPAC3Pww93pOCHCGEKSSBl1bGWVg+DrZ+DX6NYXA03NDZ7KiEEC5MEnhp/L44vyDnmBTkCCHshiTwK0k9BkvHwp75EHALDPoO6t9hdlRCCAFIAi+e1rD9W4h9zXhg2fVNaP8cuLmbHZkQQhSSBH65Uwdh4Sg4FA+NQ+HuqVC7OdHbk2QfSCGEXZEEXiA3BzZ+Cqv/ZfS0+34IdzwKVaoQvT2JyKidpGfnApB0Jp3IqJ0AksSFEKaRBA5w5DejIOfIDgjsA30mQ416hS9PitlXmLwLpGfnMilmnyRwIYRpXDuBZ6fDmn/D+o/Auxbc/19oeU+RgpzkM+nFvr2k40IIURlcN4EfXmfskHPqTwh6BLq/W2JBTj0/C0nFJOt6fpaKjlIIIUrkeisupZ8xHlLO7AM6F4bMh3s+uWI15djwQCzublbHLO5ujA0PrOhohRCiRK7VA9+7EBa/COdToP2zEPYqeHhf9W0F49wyC0UIYU9cI4GnHoUlY2HvAghoDQ/NgXpB13SJiKD6krCFEHbFuRO41rD9G4h9Pb8gZ5zR85aCHCGEE3DeBH7yT2Os+/BaaNwhvyDnRrOjEkIIm3G+BJ6bAxs+hrgJ4OZhJO6gIbJDjhDC6ThXAj/yK8wfCUd/g5v6Qu/JUON6s6MSQogK4RwJPDvd6HH//DFUqw0DZxkFOUII4cQcP4EfWmvskHPqIAQNhh7vGtucCSGEk3PcBJ5+Bpa/AdtmwXVNYcgCaNbJ7KiEEKLSOGYC37MAlrwI509A6Cjo9EqpCnKEEMKZOFYCP3fE2CFn70Ko2xoe+hHq3W52VEIIYQrHSOBac31yLGwYArmZ0O1taPeMFOQIIVxauSZHK6V6KqX2KaUOKKVesVVQl4rensSad8MJ3P8J27IbsrxTFHR4XpK3EMLllbkHrpRyAz4BugOJwGal1AKt9R5bBVewE06v3CCWqFv5ITcMS2wqE6onybokQgiXV54eeFvggNb6oNY6C5gD2HTydcFOOFF5HfkhtzOgCnfCEUIIV1eeMfD6wN+XfJwI3Hn5SUqp4cBwgICAAOLi4kp9g+I2USg4fi3XcUZpaWku3wYFpC2sSXtYc+b2KE8CV8Uc00UOaD0NmAYQHBysw8LCSn2D+htXFZvE6/tZuJbrOKO4uDiXb4MC0hbWpD2sOXN7lGcIJRFoeMnHDYDk8oVjTXbCEUKIkpWnB74ZaK6UagokAYOAh2wSVb5Ld8JJOpNOfdkJRwghCpU5gWutc5RSI4EYwA2YobXebbPI8hXshOPMvwYJIURZlKuQR2u9BFhio1iEEEJcA9nlQAghHJQkcCGEcFCSwIUQwkFJAhdCCAeltC5Se1NxN1PqOPBXGd9eGzhhw3AcnbTHRdIW1qQ9rDlDezTWWvtffrBSE3h5KKW2aK2DzY7DXkh7XCRtYU3aw5ozt4cMoQghhIOSBC6EEA7KkRL4NLMDsDPSHhdJW1iT9rDmtO3hMGPgQgghrDlSD1wIIcQlJIELIYSDcogEXhmbJzsCpVRDpdRqpdRepdRupdQos2OyB0opN6XUdqXUIrNjMZtSyk8pNVcp9Xv+90k7s2Myi1JqdP7PyS6l1PdKKS+zY7I1u0/gl2ye3AtoCTyolGppblSmyQFe0FrfDIQAz7hwW1xqFLDX7CDsxFRgmdb6JuA2XLRdlFL1geeAYK31LRhLXg8yNyrbs/sETiVsnuwotNZHtNbb8v+divHD6dK7WyilGgB9gK/MjsVsSqkaQEdgOoDWOktrfcbcqExVFbAopaoC3th4xzB74AgJvLjNk106aQEopZoAQcAmcyMx3RTgJSDP7EDsQDPgOPB1/pDSV0qpamYHZQatdRIwGUgAjgBntdax5kZle46QwEu1ebIrUUpVB+YBz2utz5kdj1mUUn2BFK31VrNjsRNVgTuAz7TWQcB5wCWfGSmlrsP4Tb0pUA+oppR6xNyobM8REniFb57sSJRS7hjJe7bWOsrseEwWCvRTSh3GGFrropT61tyQTJUIJGqtC34rm4uR0F1RN+CQ1vq41jobiALamxyTzTlCAi/cPFkp5YHxIGKByTGZQimlMMY392qtPzA7HrNprSO11g201k0wvi9Waa2drpdVWlrro8DfSqnA/ENdgT0mhmSmBCBEKeWd/3PTFSd8oFuuPTErQ2VtnuwgQoHBwE6l1I78Y6/m700qBMCzwOz8zs5BYJjJ8ZhCa71JKTUX2IYxe2s7TlhSL6X0QgjhoBxhCEUIIUQxJIELIYSDkgQuhBAOShK4EEI4KEngQgjhoCSBCyGEg5IELoQQDur/AQfauSyvg09JAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "a_exact = 1.\n",
    "b_exact = 2.\n",
    "\n",
    "n = 10\n",
    "xi = jnp.arange(n)\n",
    "yi = a_exact * xi + b_exact + 2 * jax.random.normal(jax.random.PRNGKey(1), (n, ))\n",
    "\n",
    "A = jnp.array([xi, jnp.ones(n)])\n",
    "coef = jnp.linalg.pinv(A).T.dot(yi) # coef is [a, b]\n",
    "\n",
    "plt.plot(xi, yi, 'o', label='$(x_i, y_i)$')\n",
    "plt.plot(xi, coef[0]*xi + coef[1], label='Least squares')\n",
    "plt.legend(loc='best', fontsize=18)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Lacking for structure\n",
    "\n",
    "- A typical 3D-problem requires a $100 \\times 100 \\times 100$ discretization  \n",
    "\n",
    "- This gives a linear system with $10^6$ unknowns, right-hand side takes $8$ megabytes of memory\n",
    "\n",
    "- This matrix has $10^6 \\times 10^6 = 10^{12}$ elements, takes $8$ terabytes of memory.\n",
    "\n",
    "Fortunately, the matrices in real-life are not **dense**, but have certain **structure**:\n",
    "\n",
    "- Sparse matrices\n",
    "- Low-rank matrices\n",
    "- Toeplitz matrices (shift-invariant property)\n",
    "- Sparse in certain bases"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "- Linear systems can be solved by Gaussian elimination, complexity is $\\mathcal{O}(n^3)$.\n",
    "- Linear systems can be solved by LU-decomposition, complexity is $\\mathcal{O}(n^3)$ for the decomposition, $\\mathcal{O}(n^2)$ for each solve\n",
    "- Linear least squares can be solved by normal equation (bad)\n",
    "- Linear least squares can be solved by QR-decomposition (good) or by augmentation (not bad)\n",
    "- Without structure, we can solve up to $10^4$ linear systems on a laptop (memory restrictions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- Eigenvectors & eigenvalues \n",
    "- Schur theorem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 8,
     "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.7.7"
  },
  "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
}