lectures/ludecomposition.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2640a05c",
   "metadata": {},
   "source": [
    "# LU Decomposition\n",
    "$\\newcommand{\\AAA}{\\mathbf{A}}$\n",
    "$\\newcommand{\\BBB}{\\mathbf{b}}$\n",
    "$\\newcommand{\\XXX}{\\mathbf{x}}$\n",
    "$\\newcommand{\\UUU}{\\mathbf{U}}$\n",
    "$\\newcommand{\\YYY}{\\mathbf{y}}$\n",
    "$\\newcommand{\\LLL}{\\mathbf{L}}$\n",
    "\n",
    "LU decomposition is short for lower (L) upper (U) decomposition. It denotes the decomposition of a matrix $\\mathbf{A}$ into a lower $\\mathbf{L}$ and upper $\\mathbf{U}$ matrix, so that the product of the lower and upper matrices gives the matrix $\\mathbf{A}$:\n",
    "\\begin{equation}\n",
    "\\mathbf{A}  = \\mathbf{L}\\mathbf{U}\n",
    "\\end{equation}\n",
    "LU decomposition is a common technique in linear algebra to solve systems of linear equations. It is thus an alternative method to the conjugate gradient method we considered in the previous section.\n",
    "\n",
    "If a matrix $\\mathbf{A}$ can be written as the product of a lower tringular matrix, $\\mathbf{L}$, and an upper triangular matrix $\\mathbf{U} $ the corresponding equation system\n",
    "\\begin{equation}\n",
    "\\mathbf{A}\\mathbf{x} = \\mathbf{b}\n",
    "\\end{equation}\n",
    "can be solved by reducing to two (simpler) systems. First substitute \n",
    "\\begin{eqnarray}\n",
    "  \\mathbf{A}=\\mathbf{LU}\n",
    "\\end{eqnarray}  \n",
    "to get\n",
    "\\begin{eqnarray}\n",
    "  \\mathbf{L}\\mathbf{U}\\mathbf{x} =\\mathbf{b}.\n",
    "\\end{eqnarray}\n",
    "Then introduce \n",
    "\\begin{eqnarray}\n",
    "  \\mathbf{y} = \\mathbf{U}\\mathbf{x} \n",
    "\\end{eqnarray}\n",
    "to get\n",
    "\\begin{eqnarray}\n",
    "  \\mathbf{L}\\mathbf{y}=\\mathbf{b}.\n",
    "\\end{eqnarray}\n",
    "This equation is easy to solve, since $\\mathbf{L}$ is triangular. As we will see below, this can be solved by forward substitution. We therefore first solve for $\\mathbf{y}$ in the equation above. We then solve for $\\mathbf{x}$\n",
    "\\begin{eqnarray}\n",
    "\\UUU\\XXX = \\YYY\n",
    "\\end{eqnarray}\n",
    "This is also an easy system to solve, since $\\UUU$ is triangular. As will be shown below, this can be solved by back substitution.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ee91b60",
   "metadata": {},
   "source": [
    "## Gaussian Elimination\n",
    "$\\newcommand{\\AAA}{\\mathbf{A}}$\n",
    "$\\newcommand{\\III}{\\mathbf{I}}$\n",
    "$\\newcommand{\\BBB}{\\mathbf{b}}$\n",
    "$\\newcommand{\\XXX}{\\mathbf{x}}$\n",
    "\n",
    "Gaussian elimination, sometimes called Gauss-Jordan elimination, is a method to solve systems of linear equations by a sequence of elementary row operations. The method is named after Carl Friedrich Gauss (and Wilhelm Jordan). This method is taught in any linear algebra course.\n",
    "\n",
    "Consider a linear system of equations written as:\n",
    "\\begin{eqnarray}\n",
    "    \\begin{matrix}\n",
    "    a_{1,1}x_1 & + &  a_{1,2}x_2 & + &\\cdots & +a_{1,n}x_n & = & b_1 \\\\\n",
    "    a_{2,1}x_1 & + &  a_{2,2}x_2 & + &\\cdots & +a_{2,n}x_n & = & b_2 \\\\\n",
    "    \\vdots     &   &  \\vdots     & \\ddots & \\vdots  \\\\\n",
    "    a_{n,1}x_1 & + &  a_{n,2}x_2 & + &\\cdots & +a_{n,n}x_n & = & b_n \\\\\n",
    "    \\end{matrix}\n",
    "\\end{eqnarray}\n",
    "We can write this system of linear equations using matrix form as:\n",
    "\\begin{eqnarray}\n",
    " \\AAA\\XXX = \\BBB\n",
    "\\end{eqnarray}\n",
    "where\n",
    "\\begin{eqnarray}\n",
    " \\AAA = \n",
    "    \\begin{pmatrix}\n",
    "    a_{1,1} & a_{1,2} & \\cdots & a_{1,n} \\\\\n",
    "    a_{2,1} & a_{2,2} & \\cdots & a_{2,n} \\\\\n",
    "    \\vdots  & \\vdots  & \\ddots & \\vdots  \\\\\n",
    "    a_{n,1} & a_{n,2} & \\cdots & a_{n,n} \n",
    "    \\end{pmatrix}\n",
    "\\end{eqnarray}\n",
    "and\n",
    "\\begin{eqnarray}\n",
    "    \\XXX =\n",
    "    \\begin{pmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_n\n",
    "    \\end{pmatrix},\n",
    "    \\,\\,\\,\n",
    "    \\BBB =\n",
    "    \\begin{pmatrix}\n",
    "    b_1 \\\\\n",
    "    b_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    b_n\n",
    "    \\end{pmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "#### Example\n",
    "\n",
    "As an example, consider the equation system\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "      &   &  2B & + & 8C & = & 0  \\\\ \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "This can be written as\n",
    "\\begin{eqnarray}\n",
    "  \\AAA \\XXX = \\BBB\n",
    "\\end{eqnarray}\n",
    "where\n",
    "\\begin{eqnarray}\n",
    "\\AAA=\n",
    "\\begin{pmatrix}\n",
    "    1 &  1  &  1  \\\\\n",
    "    0 &  1  &  3  \\\\\n",
    "    0 &  2  &  8  \\\\ \n",
    "\\end{pmatrix},\n",
    "\\,\\,\\,\n",
    "\\XXX = \n",
    "\\begin{pmatrix}\n",
    "A  \\\\\n",
    "B \\\\\n",
    "C  \\\\\n",
    "\\end{pmatrix},\n",
    "\\,\\,\\,\n",
    "\\BBB = \n",
    "\\begin{pmatrix}\n",
    "1  \\\\\n",
    "-2 \\\\\n",
    "0  \\\\\n",
    "\\end{pmatrix}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cd9e6e3a",
   "metadata": {},
   "source": [
    "### Elementary row operations\n",
    "\n",
    "There are three types of elementary row operations. These are:\n",
    "1. **Interchange**: Order of equations changed\n",
    "2. **Scaling**:     Multiplication with a scalar\n",
    "3. **Replacement**:  Equation replaced with the sum of itself and a multiple of any other equation\n",
    "\n",
    "Applying these operations to a linear system of equations does not change its solution. Thus, we can use this operations on the set of linear operations to obtain a system that is easier to solve.\n",
    "\n",
    "#### Example\n",
    "\n",
    "Consider the following matrix:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "    A & + &  2B & + & 4C & = & -1 \\\\\n",
    "    A & + &  3B & + & 9C & = & 1  \\\\ \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "We will now reduce this matrix to an upper triangular matrix using elementary row operations. In fact, we will only use the replacement operation multiple times.\n",
    "\n",
    "Replace second equation with the second equation minus the first equation:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "    A & + &  3B & + & 9C & = & 1  \\\\ \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Replace third equation with: third equation - first equation\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "      &   &  2B & + & 8C & = & 0  \\\\ \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Replace third equation with third equation - 2 $\\times$ second equation\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "      &   &     &   & 2C & = & 4  \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "We have then obtained an upper triangular system. Such a system can easily be solved by back substitution. Let us therefore consider back substitution next."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0dc79df6",
   "metadata": {},
   "source": [
    "###  Forward and back substitution\n",
    "\n",
    "Equations of the form $\\LLL\\XXX = \\BBB$, where $\\LLL$ is an upper triangular matrix, can easily be solved using a method called forward substitution. With this method, if $\\XXX = [x_1, x_2, \\dots, x_n]^T$, you first compute $x_1$ explicitly, since this is the only unknown in the first equation. Then you substitute $x_1$ into the second equation, so that $x_2$ is the only unknown in the second equation, and then solve for $x_2$. Then both $x_1$ and $x_2$ is substituted into the third equation, so that $x_3$ is the only remaining unknown, and thereby can be solved for. This process is continued until the last unknown $x_n$.\n",
    "\n",
    "Similarly, just going the other way, equations of the form $\\UUU\\XXX= \\BBB$, where $\\UUU$ is an upper triangular matrix, can be solved by a method called back substitution. The following example shows the back substitution using elementary row operations:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "      &   &     &   & 2C & = & 4  \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "Divide the third equation by 2 to obtain C:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C  & = & 1  \\\\\n",
    "      &   &  B  & + & 3C & = & -2 \\\\\n",
    "      &   &     &   & C  & = & 2  \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "Insert the value for $C$ in the second equation and solve for B\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A & + &  B  & + & C     & = & 1  \\\\\n",
    "      &   &  B  &   &       & = &  -8 \\\\\n",
    "      &   &     &   & C     & = & 2  \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "Insert value for $C$ and $B$ in the first equation and solve for A\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    A &   &     &   &    & = & 7  \\\\\n",
    "      &   &  B  &   &    & = & -8 \\\\\n",
    "      &   &     &   & C  & = & 2  \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "194d5888",
   "metadata": {},
   "source": [
    "### Augmented matrix\n",
    "It is practical to define an \\emph{augmented} matrix by inserting the right-hand vector $\\BBB$ as an extra (augmented) column in $\\AAA$\n",
    "\\begin{eqnarray}\n",
    "[\\AAA | \\BBB]\n",
    "\\end{eqnarray}\n",
    "Writing out the elements, we have\n",
    "\\begin{eqnarray}\n",
    " [\\AAA | \\BBB] = \n",
    "    \\left(\n",
    "    \\begin{matrix}\n",
    "    a_{1,1} & a_{1,2} & \\cdots & a_{1,n} \\\\\n",
    "    a_{2,1} & a_{2,2} & \\cdots & a_{2,n} \\\\\n",
    "    \\vdots  & \\vdots  & \\ddots & \\vdots  \\\\\n",
    "    a_{n,1} & a_{n,2} & \\cdots & a_{n,n} \n",
    "    \\end{matrix}\n",
    "    \\,\\,\\,\n",
    "    \\left| \n",
    "    \\begin{matrix}\n",
    "     b_1 \\\\\n",
    "     b_2 \\\\\n",
    "     \\vdots \\\\\n",
    "     b_n\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right)\n",
    "\\end{eqnarray}\n",
    "\n",
    "With an augmented matrix, we can solve the linear system $\\AAA \\XXX = \\BBB$ quite efficiently the following way: We use the elementary row operations in Gaussian elimination to reduce the matrix $\\AAA$ to the identity matrix $\\III$. During these elementary row operations we have then transformed the vector $\\BBB$ into the solution vector $\\XXX$.\n",
    "\n",
    "\n",
    "### Example\n",
    "\n",
    "As an example, we will try to find the solution of the equation system\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "    x_1 & + & 2x_2 + & x_3 + & 4x_4 & = & 13 \\\\\n",
    "   2x_1 & + & 0x_2 + & 4x_3 + & 3x_4 & = & 28 \\\\\n",
    "    4x_1 & + & 2x_2 + & 2x_3 + & x_4 & = & 20 \\\\\n",
    "    -3x_1 & + & x_2 + & 3x_3 + & 2x_4 & = & 6 \\\\\n",
    "\\end{matrix}    \n",
    "\\end{eqnarray}\n",
    "\n",
    "This can be written as an augmented matrix in the following form:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{21} = 2               \\\\\n",
    "m_{31} = 4               \\\\\n",
    "m_{41} = -3              \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4  \\\\\n",
    "    2 & 0 & 4 & 3  \\\\\n",
    "    4 & 2 & 2 & 1  \\\\\n",
    "   -3 & 1 & 3 & 2  \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13 \\\\\n",
    "28 \\\\\n",
    "20 \\\\\n",
    "6  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "The first row is called the **pivotal** row, and $a_{11}=1$ the **pivot** element. The multipliers $m_{21}, m_{31}$ and $m_{41}$ are used for multiplying row 1 and subtracting from rows 2-4.\n",
    "These elementary row operations will eliminate all matrix elements below the pivot element.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{21} = 2               \\\\\n",
    "m_{31} = 4               \\\\\n",
    "m_{41} = -3              \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4  \\\\\n",
    "    2 & 0 & 4 & 3  \\\\\n",
    "    4 & 2 & 2 & 1  \\\\\n",
    "   -3 & 1 & 3 & 2  \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13 \\\\\n",
    "28 \\\\\n",
    "20 \\\\\n",
    "6  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "After the elimination we have\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{21} = 2               \\\\\n",
    "m_{31} = 4               \\\\\n",
    "m_{41} = -3              \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4        \\\\\n",
    "    0 & -4 & 2 & -5      \\\\\n",
    "    0 & -6 & -2 & -15    \\\\\n",
    "    0 & 7 & 6 & 14        \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13  \\\\\n",
    "2   \\\\\n",
    "-32 \\\\\n",
    "45  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "We can now continue with the same process, now on the second row:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\,                       \\\\\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{31} = 6/4            \\\\\n",
    "m_{41} = -7/4             \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4        \\\\\n",
    "    0 & -4 & 2 & -5      \\\\\n",
    "    0 & -6 & -2 & -15    \\\\\n",
    "    0 & 7 & 6 & 14        \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13  \\\\\n",
    "2   \\\\\n",
    "-32 \\\\\n",
    "45  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "After elimination we then have\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\,                       \\\\\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{31} = 6/4            \\\\\n",
    "m_{41} = -7/4             \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4        \\\\\n",
    "    0 & -4 & 2 & -5      \\\\\n",
    "    0 & 0  & -5  & -7.5    \\\\\n",
    "    0 & 0 & 9.5 & 5.25        \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13  \\\\\n",
    "2   \\\\\n",
    "-35 \\\\\n",
    "48.5  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "Now, for the third row, we are doing the same again:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\,                       \\\\\n",
    "\\,            \\\\\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{41} = -9.5/5             \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4        \\\\\n",
    "    0 & -4 & 2 & -5      \\\\\n",
    "    0 & 0  & -5  & -7.5    \\\\\n",
    "    0 & 0 & 9.5 & 5.25        \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13  \\\\\n",
    "2   \\\\\n",
    "-35 \\\\\n",
    "48.5  \\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "After elimination we then have\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "\\,                       \\\\\n",
    "\\,            \\\\\n",
    "\\rightarrow \\mbox{pivot} \\\\\n",
    "m_{41} = -9.5/5             \\\\\n",
    "\\end{matrix}\n",
    "\\left[\n",
    "\\begin{matrix}\n",
    "    1 & 2 & 1 & 4        \\\\\n",
    "    0 & -4 & 2 & -5      \\\\\n",
    "    0 & 0  & -5  & -7.5    \\\\\n",
    "    0 & 0 &  & -9.0        \\\\\n",
    "\\end{matrix}    \n",
    "\\,\\,\\,\n",
    "\\left| \n",
    "\\begin{matrix}\n",
    "13  \\\\\n",
    "2   \\\\\n",
    "-35 \\\\\n",
    "-18\\\\\n",
    "\\end{matrix}\n",
    "\\right.\n",
    "\\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "We can now either continue with our elementary row operations to get to $[\\III | \\XXX]$, or we can use back substitution directly to find the the solution: $x_4 = 2, x_3=4, x_2 = -1, x_1 =3$.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa310c1b",
   "metadata": {},
   "source": [
    "## Gaussian Elimination algorithm\n",
    "\n",
    "We will now formulate what was conducted in the example above as an algorithm. For this, assume an augmented matrix of size $N\\times N+1$ with the vector $\\BBB$ stored in column $N+1$:\n",
    "\\begin{eqnarray}\n",
    "\\AAA\\XXX =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & \\cdots & a^{(1)}_{1,N} \\\\\n",
    "    a^{(1)}_{2,1} & a^{(1)}_{2,2} & \\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots  & \\vdots  & \\ddots & \\vdots  \\\\\n",
    "    a^{(1)}_{N,1} & a^{(1)}_{N,2} & \\cdots & a^{(1)}_{N,N} \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_N\n",
    "    \\end{bmatrix}\n",
    "    =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(1)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(1)}_{N,N+1}\n",
    "    \\end{bmatrix}\n",
    "    = \\BBB\n",
    "\\end{eqnarray}\n",
    "\n",
    "We will now present the method to construct an upper triangular system\n",
    "\\begin{eqnarray}\n",
    "\\UUU \\XXX =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    0             & a^{(2)}_{2,2} & a^{(2)}_{2,3}  &\\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    0             & 0             & 0              &\\cdots & a^{(N)}_{N,N} \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_N\n",
    "    \\end{bmatrix}\n",
    "    =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(2)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(N)}_{N,N+1}\n",
    "    \\end{bmatrix}\n",
    "    = \\YYY\n",
    "\\end{eqnarray}\n",
    "\n",
    "The algorithm for this method is as follows:\n",
    "\n",
    "\n",
    "**Step one**: Create agumnted matrix and store coefficients\n",
    "\\begin{eqnarray}\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    a^{(1)}_{2,1} & a^{(1)}_{2,2} & a^{(1)}_{2,3}  &\\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    a^{(1)}_{N,1} & a^{(1)}_{N,2} & a^{(1)}_{N,3}  &\\cdots & a^{(1)}_{N,N} \n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\left.\n",
    "    \\left|\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(1)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(1)}_{N,N+1}\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "The superscript $a^{(1)}_{rc}$ means that the coefficients are stored\n",
    "for the first time in row $r$ and column $c$.\n",
    "\n",
    "\\begin{eqnarray}\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    a^{(1)}_{2,1} & a^{(1)}_{2,2} & a^{(1)}_{2,3}  &\\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    a^{(1)}_{N,1} & a^{(1)}_{N,2} & a^{(1)}_{N,3}  &\\cdots & a^{(1)}_{N,N} \n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\left.\n",
    "    \\left|\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(1)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(1)}_{N,N+1}\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "**Step two**: If $a_{11}=0$, switch rows. Then eliminate\n",
    "the elements below the pivot element.\n",
    "\n",
    "- for r=2:N \n",
    "    - $m_{r1} = a^{(1)}_{r1}/a^{(1)}_{11}$, $a^{(2)}_{r1} = 0$               \n",
    "        - for c=2:N+1        \n",
    "            - $a^{(2)}_{rc} = a^{(1)}_{rc}-m_{r1}*a^{(1)}_{1c}$\n",
    "         \n",
    "The result after step two is:\n",
    "\\begin{eqnarray}\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    0             & a^{(2)}_{2,2} & a^{(2)}_{2,3}  &\\cdots & a^{(2)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    0             & a^{(2)}_{N,2} & a^{(2)}_{N,3}  &\\cdots & a^{(2)}_{N,N} \n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\left.\n",
    "    \\left|\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(2)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(2)}_{N,N+1}\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "**Step three**: If $a_{22}=0$ switch rows with one below, then eliminate elements in the column below the new pivot $a_{22}$\n",
    "        \n",
    "- for r=3:N \n",
    "    - $m_{r2} = a^{(2)}_{r2}/a^{(2)}_{22}$, $a^{(3)}_{r2} = 0$                    \n",
    "        - for c=3:N+1\n",
    "            - $a^{(3)}_{rc} = a^{(2)}_{rc}-m_{r2}*a^{(2)}_{2c}$             \n",
    "\n",
    "The result after step three is:\n",
    "\\begin{eqnarray}\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    0             & a^{(2)}_{2,2} & a^{(2)}_{2,3}  &\\cdots & a^{(2)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    0             & 0             & a^{(3)}_{N,3}  &\\cdots & a^{(3)}_{N,N} \n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\left.\n",
    "    \\left|\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(2)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(3)}_{N,N+1}\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "**General step p+1**: If $a_{pp}=0$ switch rows with one below, then eliminate elements in the column below the new pivot $a_{pp}$.\n",
    "        \n",
    "- for r=p+1:N \n",
    "    - $m_{rp} = a^{(p)}_{rp}/a^{(p)}_{pp}$, $a^{(p+1)}_{rp} = 0$\n",
    "        - for c=p+1:N+1 \n",
    "            - $a^{(p+1)}_{rc} = a^{(p)}_{rc}-m_{rp}*a^{(p)}_{pc}$\n",
    "         \n",
    "The final result will then be\n",
    "\\begin{eqnarray}\n",
    "    \\left[\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    0             & a^{(2)}_{2,2} & a^{(2)}_{2,3}  &\\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    0             & 0             & 0              &\\cdots & a^{(N)}_{N,N} \n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\left.\n",
    "    \\left|\n",
    "    \\begin{matrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(2)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(N)}_{N,N+1}\n",
    "    \\end{matrix}\n",
    "    \\right.\n",
    "    \\right]\n",
    "\\end{eqnarray}\n",
    "\n",
    "This system has the form\n",
    "\\begin{eqnarray}\n",
    "\\UUU\\XXX = \\YYY\n",
    "\\end{eqnarray}\n",
    "where\n",
    "\\begin{eqnarray}\n",
    "\\UUU \\XXX =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,1} & a^{(1)}_{1,2} & a^{(1)}_{1,3}  &\\cdots & a^{(1)}_{1,N} \\\\\n",
    "    0             & a^{(2)}_{2,2} & a^{(2)}_{2,3}  &\\cdots & a^{(1)}_{2,N} \\\\\n",
    "    \\vdots        & \\vdots        & \\ddots         &\\vdots & \\vdots        \\\\\n",
    "    0             & 0             & 0              &\\cdots & a^{(N)}_{N,N} \n",
    "    \\end{bmatrix}\n",
    "    \\begin{bmatrix}\n",
    "    x_1 \\\\\n",
    "    x_2 \\\\\n",
    "    \\vdots \\\\\n",
    "    x_N\n",
    "    \\end{bmatrix}\n",
    "    =\n",
    "    \\begin{bmatrix}\n",
    "    a^{(1)}_{1,N+1} \\\\\n",
    "    a^{(2)}_{2,N+1} \\\\\n",
    "    \\vdots \\\\\n",
    "    a^{(N)}_{N,N+1}\n",
    "    \\end{bmatrix}\n",
    "    = \\YYY\n",
    "\\end{eqnarray}\n",
    "\n",
    "This system can be solved by back substitution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95157037",
   "metadata": {},
   "source": [
    "## LU Decomposition example\n",
    "\n",
    "The matrix $\\AAA=\\LLL\\UUU$ can be decomposed using the Gaussian elimination\n",
    "multipliers:\n",
    "\n",
    " \\begin{eqnarray}\n",
    " \\begin{bmatrix}\n",
    " a_{11} & a_{12} & a_{13} & a_{14} \\\\ \n",
    " a_{21} & a_{22} & a_{23} & a_{24} \\\\ \n",
    " a_{31} & a_{32} & a_{33} & a_{34} \\\\ \n",
    " a_{41} & a_{42} & a_{43} & a_{44} \\\\\n",
    " \\end{bmatrix}\n",
    "=\n",
    "\\begin{bmatrix}\n",
    "1      & 0      & 0      & 0      \\\\ \n",
    "m_{21} & 1      & 0      & 0      \\\\ \n",
    "m_{31} & a_{32} & 1      & 0      \\\\ \n",
    "m_{41} & m_{42} & m_{43} & 1      \\\\ \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "u_{11} & u_{12} & u_{13} & u_{14} \\\\ \n",
    "0      & u_{22} & u_{23} & u_{24} \\\\ \n",
    "0      & 0      & u_{33} & u_{34} \\\\ \n",
    "0      & 0      & 0      & u_{44} \\\\ \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "In equation form\n",
    "first solve the equation system $\\LLL\\YYY=\\BBB$ for $\\YYY$:\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "y_1        &   &            &     &           &          &     & = &b_1 \\\\\n",
    "m_{21}y_1  & + & y_2        &     &           &          &     & = &b_2 \\\\\n",
    "m_{31}y_1  & + & m_{32}y_2  & +   & y_3       &          &     & = &b_3 \\\\\n",
    "m_{41}y_1  & + & m_{42}y_2  & +   & m_{43}y_3 & +        & y_4 & = & b_4\n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "Then solve $\\UUU\\XXX=\\YYY$ for $\\XXX$.\n",
    "\\begin{eqnarray}\n",
    "\\begin{matrix}\n",
    "u_{11}x_1  & + & u_{12}x_2  & + & u_{13}x_3 &+ & u_{14}x_4  & = &y_1 \\\\\n",
    "           &   & u_{22}x_2  & + & u_{23}x_3 &+ & u_{24}x_4  & = &y_2 \\\\\n",
    "           &   &            &   & u_{33}x_3 &+ & u_{34}x_4  & = &y_3 \\\\\n",
    "           &   &            &   &           &  & u_{44}x_4  & = &y_4 \n",
    "\\end{matrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "## Example\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\AAA =\n",
    "\\begin{bmatrix}\n",
    "4  & 3 & -1 \\\\\n",
    "-2 & -4 & 5  \\\\\n",
    "1  & 2  & 6 \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "Construct LU deomposition by Gaussian elimination and\n",
    "saving multipliers\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\AAA = \n",
    "\\begin{bmatrix}\n",
    "1  &  0  & 0 \\\\\n",
    "0  &  1  & 0 \\\\\n",
    "0  &  0  & 1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "4  & 3 & -1 \\\\\n",
    "-2 & -4 & 5  \\\\\n",
    "1  & 2  & 6 \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\AAA = \n",
    "\\begin{bmatrix}\n",
    "1  &  0  & 0 \\\\\n",
    "-0.5  &  1  & 0 \\\\\n",
    "0.25  &  0  & 1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "4  & 3 & -1 \\\\\n",
    "0 & -2.5 & 4.5  \\\\\n",
    "0  & 1.25  & 6.25 \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\AAA = \n",
    "\\begin{bmatrix}\n",
    "1  &  0  & 0 \\\\\n",
    "-0.5  &  1  & 0 \\\\\n",
    "0.25  &  0  & 1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "4  & 3 & -1 \\\\\n",
    "0 & -2.5 & 4.5  \\\\\n",
    "0  & 1.25  & 6.25 \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\AAA = \n",
    "\\begin{bmatrix}\n",
    "1  &  0  & 0 \\\\\n",
    "-0.5  &  1  & 0 \\\\\n",
    "0.25  &  -0.5  & 1 \n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "4  & 3 & -1 \\\\\n",
    "0 & -2.5 & 4.5  \\\\\n",
    "0  & 0  & 8.5 \n",
    "\\end{bmatrix}\n",
    "\\end{eqnarray}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6453507b",
   "metadata": {},
   "source": [
    "# Python libraries \n",
    "\n",
    "Python have several libraries which can be used for solving linear equations. The numpy routine lin.solve gives an exact solution using LU decomposition. The scipy.linalg.lu routine performs LU decomposition and can be used\n",
    "in conjunction with scipy.linalg.solve_triangular to solve a linear equation system.\n",
    "The conjugate gradient solver scipy.linalg.cg can also be used to solve a linear equation system.\n",
    "We will use all of these routines and compare the performance in terms of accuracy and compute time.\n",
    "\n",
    "## Test case\n",
    "\n",
    "For testing the Python solvers we will use the the finite-difference solution of the two-dimensional\n",
    "laplace equation given by:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "\\frac{\\partial f(x,y)}{\\partial x} + \\frac{\\partial f(x,y)}{\\partial y} = 0,\n",
    "\\end{eqnarray}\n",
    "where $f(x,y)$ is the function we want to solve for and $x,y$ are corrdinates.\n",
    "\n",
    "The Laplace equation is discretized on a 2D grid using $x=i\\Delta x$ and $y=j\\Delta x$, \n",
    "where $\\Delta x$ is the size of each grid cell and $i,j$ are integers ranging from $0$ to $N$.\n",
    "Assuming the phyiscal length of the grid in the $x$ and $y$-directions, then $N=L/\\delta x +1$.\n",
    "The finite-difference equations can be derived as for the one-dimensional case.\n",
    "We denote $f_{i,j} = f(x,y)$ and obtain the following equations for the\n",
    "interior points of the grid:\n",
    "\n",
    "\\begin{eqnarray}\n",
    "f_{i+1,j} -2f_{i,j} + f_{i-1,j} +f_{i,j+1} -2f_{i,j} + f_{i,j-1} = 0.\n",
    "\\end{eqnarray}\n",
    "Here $i,j$ range from 1 to $N-2$.\n",
    "The boundary conditions are inserted around the edge of the grid \n",
    "and are zero, except for the top boundary.\n",
    "\n",
    "The finte-difference equation above can be reorganized as a matrix equation of the type\n",
    "\\begin{eqnarray}\n",
    "\\AAA\\mathbf{x} = \\mathbf{b}.\n",
    "\\end{eqnarray}\n",
    "\n",
    "The matrix $A$ is a square matrix :\n",
    "\n",
    "\\begin{eqnarray}\n",
    " \\AAA = \n",
    "    \\begin{pmatrix}\n",
    "    a_{1,1} & a_{1,2} & \\cdots & a_{1,n} \\\\\n",
    "    a_{2,1} & a_{2,2} & \\cdots & a_{2,n} \\\\\n",
    "    \\vdots  & \\vdots  & \\ddots & \\vdots  \\\\\n",
    "    a_{n,1} & a_{n,2} & \\cdots & a_{n,n} \n",
    "    \\end{pmatrix}\n",
    "\\end{eqnarray}\n",
    "\n",
    "\n",
    "where each element is also a matrix of the type:\n",
    "$$\\mathbf{a} = \n",
    "\\begin{bmatrix}\n",
    "-2&1&0&0&\\cdots&0&0&0 \\\\\n",
    "1&-2&1&0&\\cdots&0&0&0 \\\\\n",
    "0&1&-2&1&\\cdots&0&0&0 \\\\\n",
    "\\vdots&\\vdots&\\vdots&\\vdots& & \\vdots&\\vdots&\\vdots \\\\\n",
    "0&0&0&0&\\cdots&1&-2&1 \\\\\n",
    "0&0&0&0&\\cdots&0&1&-2 \n",
    "\\end{bmatrix}$$\n",
    "\n",
    "The size of the $\\mathbf{a}$ matrix is $N \\times N$, so the\n",
    "size of the $A$ matrix is then $N^2 \\times N^2$.\n",
    "The unknown funxtions $f_{ij}$ are contained in the $\\mathbf{x}$ vector, while\n",
    "the boundary conditions are stored in the right-hand side $\\mathbf{b}$ vector.\n",
    "\n",
    "The python code below implements the solution of the above equation using three different solvers:\n",
    "\n",
    "- numpy.linalg.solve\n",
    "- scipy.linalg.lu\n",
    "- scipy.linalg.cg\n",
    "\n",
    "We compare performance in terms of accuracy and compute time by the python code given below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "90f78473",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx :  25\n",
      "cg...\n",
      "0.0036083540035178885 0.00047392613073618664\n",
      "gesv...\n",
      "0.005037247989093885 5.727361557111548e-14\n",
      "lu...\n",
      "0.014883633994031698 4.412204528839466e-14\n",
      "nx :  50\n",
      "cg...\n",
      "0.316160246002255 0.0006355564400891403\n",
      "gesv...\n",
      "0.20946981500310358 1.3573124569578207e-13\n",
      "lu...\n",
      "0.1602477540000109 1.0365857234063024e-13\n",
      "nx :  75\n",
      "cg...\n",
      "2.321929302008357 0.0008198425992008223\n",
      "gesv...\n",
      "1.4760695299919462 2.0792628686756179e-13\n",
      "lu...\n",
      "1.206100494004204 1.495760457267951e-13\n",
      "nx :  100\n",
      "cg...\n",
      "10.448221318001742 0.0009387017108021174\n",
      "gesv...\n",
      "7.570451648993185 3.033929342303371e-13\n",
      "lu...\n",
      "8.93702895999013 2.177598106395483e-13\n",
      "nx :  125\n",
      "cg...\n",
      "34.60822415498842 0.0010157129474658946\n",
      "gesv...\n",
      "34.97466474400426 3.934113612808116e-13\n",
      "lu...\n",
      "32.25540799800365 2.763130514621049e-13\n",
      "nx :  150\n",
      "cg...\n",
      "86.06969190499512 0.0011217912602088352\n",
      "gesv...\n",
      "140.55661413399503 4.886911128027555e-13\n",
      "lu...\n",
      "99.66186075400037 3.4390251297740935e-13\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhIAAAHICAYAAAAbcANrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA97UlEQVR4nO3de3RU9b338c+QQMItqUIJdwwWK0pRGqoCpd7jUg7Wrqr02Apa9MgRRUxrBVG5PJUc7CrLC4LFIq1PwbJqa8tZC4U8tiKKVsV4qXBqq2gQE1NQE66JSfbzB2diLpNk9p59+e2936+1XG2GyWQz7Mt7vnvPJGFZliUAAAAHugW9AAAAILwICQAA4BghAQAAHCMkAACAY4QEAABwjJAAAACOERIAAMAxQgIAADhGSAAAAMcICQAA4JhRIfHcc89p6tSpGjx4sBKJhP74xz92+T1bt25VUVGRcnNzNXLkSD388MPeLygAAIbr6phqWZYWLVqkwYMHq2fPnjrnnHP09ttv2/45RoXEoUOHdNppp2nFihVp3X/37t265JJLNHnyZJWXl+uOO+7QnDlz9Pvf/97jJQUAwGxdHVPvvfdeLV++XCtWrNArr7yigQMH6sILL9SBAwds/ZyEqb+0K5FI6Mknn9Rll13W4X1uv/12bdy4Ubt27Wq+bdasWXrjjTf04osv+rCUAACYr+0x1bIsDR48WHPnztXtt98uSaqrq1NBQYGWLVumG264Ie3HzvZigf3y4osvqri4uNVtF110kdasWaPPP/9c3bt3b/c9dXV1qqura/66qalJn3zyifr166dEIuH5MgMAzGFZlg4cOKDBgwerWzdvhvRHjx5VfX29K49lWVa7Y1VOTo5ycnJsPc7u3btVVVXV6hiak5Ojs88+W9u3b49PSFRVVamgoKDVbQUFBWpoaNC+ffs0aNCgdt9TWlqqxYsX+7WIAIAQ2LNnj4YOHer64x49elTDR/TWv6qbXHm8Pn366ODBg61uW7hwoRYtWmTrcaqqqiQp5TH0gw8+sPVYoQ4JSe3KLHmmpqPpwvz581VSUtL8dU1NjYYPH65lW8erZ59sndPzPe8WFgBghGePjJQkHTnYoNvPflV9+/b15OfU19frX9VN2v7yl9WnT2ZT74MHLU0841/as2eP8vLymm+3O41oKdUx1O50PtQhMXDgwOaqSqqurlZ2drb69euX8ns6GgH934/PVdaBHL365Y8kSZO/9E67+1zQ658uLDUAwA//7/BXUt6+7bOTJElv/Guw9u/uLulVz09t9+mTUN++mZ46OTbVyMvLaxUSTgwcOFDSsclEy+l9dXV1uylFV0IdEhMmTNB///d/t7pty5YtGj9+fMrrIzr1/47Tp2dk6dmqr+q4ga2vWE1GRUcrJYEBAMHpaN+clAwH6Vg8fFp1bPpw3OvdddIL1frQ06UzU2FhoQYOHKiysjKNGzdO0rHpydatW7Vs2TJbj2VUSBw8eFD//OcXB+Xdu3fr9ddf1/HHH6/hw4dr/vz52rt3rx577DFJx96hsWLFCpWUlOj666/Xiy++qDVr1ujxxx+3/bN7VzaqripXdQM/16dVffWGBkuSTvvyR61WwlSTis5WYiIDANzRVTC01HK/LR0LCEnNEZFT1V19Pmp0b+EM1NUxde7cuVq6dKlGjRqlUaNGaenSperVq5euuuoqWz/HqJB49dVXde655zZ/nbyWYcaMGfrVr36lyspKVVRUNP95YWGhNm3apFtvvVUPPfSQBg8erAceeEDf/e53Hf38XlWS1L05Jo4beKB55Tvtf095dBUVbTHFAID02YmFttrGg/RFQEitI6JXVbu7Rk5Xx9Sf/OQnOnLkiG688UZ9+umnOvPMM7Vlyxbb14sY+zkSfqmtrVV+fr7OGzdPdYX9dHBwlg4PlOoGft58n5anOpJB0VY6UZEuIgNA1GUSDC2ligepdUBIqSOi4KVaNb67W898+mvV1NRkfN1BKsljzJs7B2R8jcSBA00ae0q1Z8vqlFETCVO0nExIap5OSGo3oUiyO6noDFMMAFHgViy01VE8SB0HhHQsIpKiflrDT4REG30+atTBwVmSjq10qWJC6jgoJHejoiWuxQBgGq9ioa3O4kFqHxBS6ojoVUVEuI2QaKF3xSEdGt5b0rGV7fDA9jEhKe2gkLyLiraIDABe8SsW2uoqHqSuA0JqHRFwHyGRQnIqkSompPbTCanroJDabxRehkVLRAaArgQVC22lEw9S6oCQuo6I5DSid8Uhh0uItgiJ/5V4v1I6sbDd7Z3FhCRHQZHk17SiM13tPAgNIDpMiYWW0g2HpI4CQuo4IuAtQqKN5OmNltdKJLWNCSn1dEKyFxSSGVGRCtMMIDxMDIVU7MaDZC8gpNYR0XYakZTYvdf2cqA9QiINyamEZC8mJPtBIZkbFW2ls9MiNgD3hCUUUnESD1LnASGlHxEtcVrDXYREJ1pOJdrGhKS0TnUkOQkKKbjrKtzCqRMgPWGOhI44jQep64CQ7EUE79TwDiGRQst3b3TG7nRCch4USWGZVqQr3Z0nwYEwi2IkpJJJOCQ5DQiJayKCQkh0oaOpRFJHMSF1PJ2QMg8KKXpR0RmCA6aJSxx0xY14kNILCCn9iOhoGsFpDfcREjalGxNS19MJyZ2gkMJ/CsQtdnfuhAeSCIP0uRUPUuYBIXUeEfAeIdFCYvdeWYVDJLU+vdH2HRxux4TkXlAkxWlakQmnBw8CxFwEgfvcDIekdANCshcRbXFthPcICRd1FhNS56c6ktwOColphRfcPljFPUw4+JvHi3iQ3AsIKXVEdDaNaHlag7d+uoeQSFM6Uwmp45iQ0p9OSN4ERRLTCvNwIEXQvAqHJDsBIbkTEUwj/EFIdKKrd284jQkpvemE5G1QSEwrgLjyOhyS7AaE5P4kAt4iJGxI9WmXTmJCsjedkLwPiiTCAoguv+JB8iYgpPTf4tl2GsG7NbxDSHjI7ZiQ/AuKJE6DAOHlZzgkOQkIKbOIYBoRLEKijZbv3JDan96wM5WQ0osJKf1THUktN9YgokIiLADTBBEOSUEEhJQ6Irg2wl+EhEsyiQnJ2XQiye8pRRJhAQQryHCQnMdDkhcRkUrb0xq8Y8NdhIQDqaYSkjsxIdmfTiQFFRRJhAXgraDDIcmPgJCcfeQ10wj/ERJpSPXujY5iojPpxISU2XRCCj4oklLt9IgLIH2mhENSpgEhuRcRXBdhDkLCZZ1NJSR7MSE5n05I5gRFS0wtgNRMi4aW/AwIyXlEpJpG8G4N7xESGXByikNKPyakzKcTUjAXZqaLqQXiyuRwkNyJhyQ/IgLBISRSaPvODSn9Xy2e5HZMSJlNJ5JMnFK0xdQCUWN6NLQUVEBImf0acK6NCA4hkSEn10ok2YkJyZ3pRFIYgiKJqQXCJEzR0JKbASF5ExF2pxGpTmvwjg33ERIe6moqITmLCcmd6YRk9mmPznS0syYw4KewRkOS2/Eg2Q8IKfOIYBoRLELCho5Ob3Q2lfAiJiR3pxNJYZpSdITpBbwQ9mBoy4uAkIiIuCIkfOBlTEjuTSeSwjql6AjTC6QrasHQlkkBIXlzOiOJd2v4h5BwSSbXSiQ5iQnJm+lEUhSmFB3p7KBBZERX1GOhLa/iIcnLiOgK0wgzEBIdSPXODcn+uzeS0plKSF9sXKZMJ5KiNqXoCpERbnGLhVRMDQgp/YjgrZ7hQEi4qKupRLoxIWU2nZC8Cwop2lOKdKRzkCI2vEMkdMzreEgyISI6m0Z0dFqDd2x4g5DwmR8xIXl7uiMpblMKO+wc7IgO4iBTYQgIiUlEVBESDnR2eiOdayX8jAnJ2+lEUtynFJlw4yAaZIwQAcHwKx4k/wIiXWG8NuLZIyPVMyuzQ+6RIw2Sqt1ZIBcREiGQSUxI/kwnkphSBIODeTz4GQ9JfkdEptMI3q3hP0LCA25PJSR3YkLyZzqRRFQAmQsiHqTMA0LyJiLCOI2IOkKiEx29c0Pq+t0bJsaEFExQSEQFYFeYA0LyfxKB4BASAQsiJiR/T3e0RVQAqQUVD0mmR0RX04jOTmvwjg3vEBIeSvdDqoKMCcn/6URLRAXiLuh4kNwLCMn9CythPkIiA04/nMoNTj+4KhUTgkIiKhAfJsSD5G5ASM4iwo9pBLxFSHjMq6lEklvTCcmcoJCICkSLKeHQkglTCLciAsEiJAxiQkxIwV4/kQpRgTAyMR4kM6YQEhdXRgkh0YXO3rkhpXd6w84v9DIpJiQzphMtERUwmanxILkfEJI/EZHONILTGsEiJAxkSkxI5gaF1H6nTVjAbyaHQ5IXASGF66JK3rHhLULCJ3Z/zbhJMSGZHRRJhAX8EIZ4kLwLCCmziHB7GoHgERIuSPfdG3Zjwik339HRVhiCIomwgBvCEg4tRSEi0sVpjeAREgZzOpVI8mo6IZl3QWY6CAukI4zhkGRqQEj2I4JpRHgQEi7xaiphekxI4ZhOpEJYIMzR0JKXASEREegcIZGGrt654TWTY0IKf1AkpTqoEBfREZVoaMnrgJD8jwg7OK1hBkIiAE6ulXAjJiRvrptIikpQtERchFMUo6ElPwJCCiYi3J5G8I4N7xESLvL6I7MzjQnJ++mEFM2gaIm4MEvUo6GlsAQE4oWQCIhf7+BIxY+YkKIfFC11djAjMjIXp1hIxa+AkNyLCK+nEZzWMAchEaAgTnEk+RUTUryCIhUiIz1xj4VU/AwIKdiIQHgREmlK94JLP34jqJsxIXl73URLLXeKcY2KttI9eIY5OAgE+8IaEJLziGAaEV6ERMCcnuJwKyYkf6cTSXGfUtjFwTge/A4IKXwRYQcXWvqDkAixsMeERFAAQcRDkgkRgfDrFvQCRJHdsVsmNe7mxptT1T2wq7U/reob6A4V8FvQ67wpEWF3/8dpDfMwkTBEJu/icHMyIQU3nZC4jgLRF3Qwu/1igUkEmEh4xO9qdntjNuF95EG/YgPcklyXg16fTdiuW+KjsKOBkLDB6wt3Mt2ovIgJE3Y8JuyAASdMWXe92pb9PKUhcVrDVJzaQJeCPNXREqc9EAYmhENLXr0YMP2UBu/Y8A8TCQ85qWfTphJJpkwnkkx5pQckmbhOmhoRTCOihYmEgTL9+Gy3L75syZTpRBJTCgTJtHBI8jL6g4gImI2Q8Jgfn3SZSpxiIomogB9MjYckkyMC0URI2JTuR2Vnyo1f6uV1TEj+fby2XUQF3GR6PEjevyPDjYhwOo3gtIbZCAmDmR4TkrnTiZaICjgRhnhICkNE+IkLLf1FSPggqNMbSX7EhGTudKIlogKdCVM8SP58LoRbEcE0IrqMe9fGypUrVVhYqNzcXBUVFWnbtm2d3n/dunU67bTT1KtXLw0aNEjXXnut9u/f79PSes+tC5P8eEVh0rs60mHKhwQhWGFdD+IQEchMQ0OD7rzzThUWFqpnz54aOXKklixZoqamJld/jlEhsWHDBs2dO1cLFixQeXm5Jk+erIsvvlgVFRUp7//8889r+vTpmjlzpt5++2397ne/0yuvvKLrrrvO5yX3VthiImxBIbU+mITtgAJ7wv5v7dc2FrbTGWhv2bJlevjhh7VixQrt2rVL9957r372s5/pwQcfdPXnGBUSy5cv18yZM3Xddddp9OjRuu+++zRs2DCtWrUq5f1feuklnXDCCZozZ44KCwv1zW9+UzfccINeffVVn5e8a6aM5/zaOYQ1KJLCfKBBe1H59/Rrm3JzP5HJCyFT9pth9eKLL+rb3/62pkyZohNOOEGXX365iouLXT9GGhMS9fX12rFjh4qLi1vdXlxcrO3bt6f8nokTJ+rDDz/Upk2bZFmWPv74Yz3xxBOaMmVKhz+nrq5OtbW1rf6zK4gLedwcDfr5SiPMMZEU9lewcRS1fzM/w9yUiEDH2h7D6urqUt7vm9/8pp555hm98847kqQ33nhDzz//vC655BJXl8eYiy337dunxsZGFRQUtLq9oKBAVVWp1+yJEydq3bp1mjZtmo4ePaqGhgZdeumlnY5tSktLtXjxYleXPV1BX3TZktcXYLYUposx05HqwMSFm8GJQih0xO8Qj8LpDFPfsbH9s6+oR0OPjB6j/mC9pJc0bNiwVrcvXLhQixYtanf/22+/XTU1NTr55JOVlZWlxsZG3XPPPfr3f//3jJajLWNCIimRSLT62rKsdrcl7dy5U3PmzNHdd9+tiy66SJWVlbrttts0a9YsrVmzJuX3zJ8/XyUlJc1f19bWtvtHMZUbbwdtyc+YkKIXFC21PZgRFt6Jcji0FPaIyHQawWmNju3Zs0d5eXnNX+fk5KS834YNG/Sb3/xG69ev16mnnqrXX39dc+fO1eDBgzVjxgzXlseYkOjfv7+ysrLaTR+qq6vbTSmSSktLNWnSJN12222SpLFjx6p3796aPHmyfvrTn2rQoEHtvicnJ6fDJ90PmU4lwh4TUjg+eyJTHR3sCIz0xSUY2gridKBpEYHO5eXltQqJjtx2222aN2+evve970mSvva1r+mDDz5QaWlpNEOiR48eKioqUllZmb7zne80315WVqZvf/vbKb/n8OHDys5u/VfIyjp2kLUsy7uFjZigYkKK5nSiMwRGe3ENhraCup7IxNMZTCPccfjwYXXr1vpSyKysLNff/mlMSEhSSUmJrr76ao0fP14TJkzQ6tWrVVFRoVmzZkk6dlpi7969euyxxyRJU6dO1fXXX69Vq1Y1n9qYO3euzjjjDA0ePNjTZfXro7JTcXsqIQUTE1J8g6Ktrg6mYQ4NQqFrUYoIphHmmDp1qu655x4NHz5cp556qsrLy7V8+XL98Ic/dPXnGBUS06ZN0/79+7VkyRJVVlZqzJgx2rRpk0aMGCFJqqysbPWZEtdcc40OHDigFStW6Ec/+pG+9KUv6bzzztOyZcuC+iukxY2LLqMUE1I8TndkwsnB2Kv4IAzcE+S7mkyNiEymEaZeaBmUBx98UHfddZduvPFGVVdXa/Dgwbrhhht09913u/pzElbMzwHU1tYqPz9f5x83Q9nd7F1Rm8lEwo13b7gdEklBxUQSQYGoC/pt0V6dzghjSDQ01euZT3+tmpqatK47sCt5jPneMz9Qjz6Zv2vjt+f/xrNldcqYz5GIGzfOAXo1Qgz6nGnYP8wK6IgJ67bJEYFwIiRCzsuYICgAd5iyLpseEVxkGU6ERAZMOR/n5SuBoGNCMmcnDNhl0rprwraMaCIkAhSW+jZlB2TSThnojEnrqtfTRVOmEaa8sIsjQiIivD4/aUpMSGbtpIGWTFs3vd5uuS4CEiERODenEnGKCcm8nTbiy8R1MUwREZbpLFIz6nMkYL4gP2uiI3yoFYJiWjwkmRb9iDYmEhly47xcmKYSkrk7qeSrQlN37ogG09czP7ZPTmmgJUIigvyKCVODQjJz1IxwC8M6FcaIcOOFFBdaBouQiCi/XjGYHBNSOHb+MFsY1iHTwx7RxjUShnDj928ExcTrJtpqeSDgWgp0xfRwaMnPgDBxGoHgMZFwgaljNT/PY4bp1ZDp57gRnLCtF2GOCEQHIWEQL+qcmOhc2A4ccF9YwzLsEeHW/s7UF3JxwqmNGPDiV453JLlzM/1UR1uc+oiXsEVDS2EMdkQbIWGYMF8r0VIYrpvoCFERTWGOh6QgIsLkaQTMQEjEhJ9TiaQwx0QSURFuUYiHpKhEBKKHayRc4uZ5Oq9qPYidQpTeltbyXHqUDlBRE8V/IyICJmMiETNBTCakaEwn2mJaYYYoBUNbQUW4lxHh5gslLrQ0AyFhKC+vlSAm3Nf2YEZYeCfK4dBSVCZ5iD5CAr4K67s67Ep1sCMunIlLOCQFHRBhmUbAHISEwaI4lUiK8nSiI0wtuha3aGgryhGB6CIkXJTYvVdW4ZCgFyNtJsSEFL+gSOrooBmHwIh7MLQVdEBI3keE29MIro8wByFhOK8/VyLomJDiOZ3oTFcH2bCEBrGQnjhEBKKNkIAxMSERFOmwe4B2IzyIAm/EJSK4NiLaCIkQ8OPTLk2ICYnphBeIAPOYEBCAW/hAKpdx3i5zUfoQK6Al09ZtphFwAyEREn5sjKadJzVphwtkyrT12bTt3Q5esJmFkEArpu1cTHsFB9hl4jrs13bONCIeCIkQ8WujNC0mJDN3xkBnTF1nTdy+EW6EhAeiMHYzdWdj4o4ZaMvU9dTP7ZppRHwQEiHj58ZpckyYuqNGvJm8bpq6PdsVhRdqUUNIoFMm73xM3mkjXkxfF/3ejplGxAufIxFCfnyuREumfMZER/gwKwTF5HhIMvnFQJz8bf8gZR3JyegxGg/XubQ07mIi4ZGojd/CsDMy/VUhoiMs61oQ2y3TiPghJEIqiI01DDEhhWcnj3AKy7oVlu3Vjqi9QIsKQgK2hGnnRFDATaxPXWMaEU+ERIgFtdGGKSYkDgBwLrnuhG394ZQG/ERIeCjKY7iwxYQUzgMCghHmdSWM2ybCjZAIuSBfBYR1hxXmgwS8FfZ1I6ht0o/9UJRfmIUdIYGMhDUmpPAfNOCeKKwLYd4WEW58jkQE+P25Em2Z/jkTXWl5AOGzKOIj7OGQFHRAcG0EmEh4zK9xXNAbc9A7M7dE4ZUpOhelf+OobHcINyYScE3YJxMtMaWIlqiEQ0smRIRfL2C4PsJshESEBH2KQ4pWTCQRFeEUxXhIilNEwHyEBFwXxZhIIirMFuV4SDIhIoCWuEbCB36O5Ux5lRCHnV1YP6woauL072DKdmXKfgZmYCIBz0R5MtEWkwp/xSEa2jIlIvzG9RHmIyQiyIRrJZLiFBNJbQ9yhEXm4hgOLZkUEUwj0BYhEVGmxYSk2AVFEmFhX9zDoSWTIgJIhZDwSWL3XlmFQ4JejEDFcTqRSqqDZJzjgmhIzcSA8HsawWmNcCAkIsykqUQSMZFaXOKCaEgPEYEwISQijpgIr84OuiZHBrGQGRMjAugMIeEjTm98gZjIjJODtZP4IAr8ZWpEMI1AZwiJGDBxKiFxEabfiAKzERGtcX1EePCBVAicqTtQwA99PmpkG0CoERIxYfpokh0p4sj09d70/QbMQEj4LMhxnek7BV6ZIU5MX9eD3F9wWiNcCAkYx/QdLJAJghlRQ0jEjOlTiSR2tIiisKzXYdlPwAyERACCHtuFZSfBKzdERZjW5bDsH2AOQgLGC8sOGEiF9deeoF9owT5CIqbC9qojTK/oACmc62zY9gswAyEREBOqO4w7jbDtmBFPYVxPw7g/gBkIiZgL484jjK/0EA+sm5kx4QUW7CMkEFrssGGSMK+PYXxBAXMQEgj1ToRXgAha2NfBMG//MINxIbFy5UoVFhYqNzdXRUVF2rZtW6f3r6ur04IFCzRixAjl5OToxBNP1KOPPurT0mbGpDFe2HcmYd+ZI5zCvs6ZtN2btD+Mkr179+oHP/iB+vXrp169eun000/Xjh07XP0ZRv32zw0bNmju3LlauXKlJk2apF/84he6+OKLtXPnTg0fPjzl91x55ZX6+OOPtWbNGn3lK19RdXW1GhoafF7yaDD1t4Tawa8nhx/CHhCIh08//VSTJk3Sueeeq6eeekoDBgzQu+++qy996Uuu/hyjQmL58uWaOXOmrrvuOknSfffdp82bN2vVqlUqLS1td/+nn35aW7du1Xvvvafjjz9eknTCCSf4ucgwEL+eHF6JUkCYNI2AN5YtW6Zhw4Zp7dq1zbd5cYw05tRGfX29duzYoeLi4la3FxcXa/v27Sm/Z+PGjRo/frzuvfdeDRkyRCeddJJ+/OMf68iRIx3+nLq6OtXW1rb6L0imjfOitHPhdAfcErV1ybTt3LT9oOnaHsPq6upS3i95jLziiis0YMAAjRs3To888ojry2PMRGLfvn1qbGxUQUFBq9sLCgpUVVWV8nvee+89Pf/888rNzdWTTz6pffv26cYbb9Qnn3zS4XUSpaWlWrx4sevLHyVROMXREhMKOBWleEgyLSLi4rOP+6hbz9yMHqPpSHdJ0rBhw1rdvnDhQi1atKjd/d977z2tWrVKJSUluuOOO/Tyyy9rzpw5ysnJ0fTp0zNalpaMCYmkRCLR6mvLstrdltTU1KREIqF169YpPz9f0rHTI5dffrkeeugh9ezZs933zJ8/XyUlJc1f19bWtvtHQfRiQiIokL4oBgSiY8+ePcrLy2v+OicnJ+X9mpqaNH78eC1dulSSNG7cOL399ttatWqVqyFhzKmN/v37Kysrq930obq6ut2UImnQoEEaMmRIc0RI0ujRo2VZlj788MOU35OTk6O8vLxW/wWNsZ6/ojamhnuivm6YOI1g/2df22NYRyExaNAgnXLKKa1uGz16tCoqKlxdHmNCokePHioqKlJZWVmr28vKyjRx4sSU3zNp0iR99NFHOnjwYPNt77zzjrp166ahQ4d6urxxYOJOx01RP2ggfXFYF6K+PaO9SZMm6e9//3ur29555x2NGDHC1Z9jTEhIUklJiX75y1/q0Ucf1a5du3TrrbeqoqJCs2bNknTstETLccxVV12lfv366dprr9XOnTv13HPP6bbbbtMPf/jDlKc1YF8cdj5xOIggtbj828dhO0Z7t956q1566SUtXbpU//znP7V+/XqtXr1as2fPdvXnGHWNxLRp07R//34tWbJElZWVGjNmjDZt2tRcT5WVla1GMn369FFZWZluvvlmjR8/Xv369dOVV16pn/70p0H9FRxL7N4rq3BI0IuRUhSvl0iFayjiIw7xkGRyRHBaw1vf+MY39OSTT2r+/PlasmSJCgsLdd999+n73/++qz8nYVmW5eojhkxtba3y8/N1/nEzlN2tR6DLYmpIJMUhJtoiKqIlTgGRREh0rqGpXs98+mvV1NR4cs1c8hgz9MHFLrxr46g+vHmhZ8vqlFETCcA0TCnCL47xkGRyRCA6CAmDmHx6Q4rPKY5UWh6MiArzxTkekkyPCBOmEXAHIQFb4hwTSUSFuQiIY0yPCESLUe/aQDgqnZ3UF+Jy1b/Jkv8G/Dscw/YJvzGRgCNMJlprexBjUuEtoiG1sEREGF4wIX2EBBwjJjrG6Q93EQ5dC0tEIHoICQOZftFlS8RE15hW2Ec42ENEIEiEBDJGTNhDWLRHOMQHpzWih5CAK4gJ51IdRKMcF0SDu5hGIGiEhKHCdHojiZhwT0cH2zAFBsHgvbBFBNOIaCIk4CpiwlvpHJz9iA0iIXhhiwhEFyEB1xETweIgH31EBEzCB1IZLMxjQHZ0gDfCum2FeX+GzhES8ExYd3iAqdimYCJCwnBhr3h2fIA72JZgKlvXSGzcuNH2D7jwwgvVs2dP29+H6OCaCSAzYY+IsL8gQudshcRll11m68ETiYT+8Y9/aOTIkba+D9FDTAD2hT0gEA+2T21UVVWpqakprf969erlxTLHTlRqnp0ikL6obC9R2X+hY7ZCYsaMGbZOU/zgBz9QXl6e7YVCdEVl5wh4ie0EYWIrJNauXau+ffumff9Vq1apf//+thcK7UWp6tlJAh2L0vYRpf0WOpbRB1IdPXpUb775pqqrq9XU1NTqzy699NKMFgzRxjUTQHtRigjEh+OQePrppzV9+nTt27ev3Z8lEgk1NvLpem4L4+/f6AwxAXyBiEBYOf4ciZtuuklXXHGFKisr211kSUQgXb0rDrEDRexFcRvgtEZ8OA6J6upqlZSUqKCgwM3lQUxFcUcKdIWQRhQ4DonLL79czz77rIuLgnREufLZoSJOory+R3k/hfYcXyOxYsUKXXHFFdq2bZu+9rWvqXv37q3+fM6cORkvHOKH6yYQB1GOCMSP45BYv369Nm/erJ49e+rZZ59VIpFo/rNEIkFIeChqF122ldzJEhSIoqhHBNOI+HEcEnfeeaeWLFmiefPmqVs3fvcX3Md0AlES9YBAfDkugPr6ek2bNo2ICEhcqp+dL6IgLutxXPZLaM1xBcyYMUMbNmxwc1mAlLiyHWHGuouoc3xqo7GxUffee682b96ssWPHtrvYcvny5RkvHDoX9Wsl2uJUB8KEgEBcOA6Jt956S+PGjZMk/e1vf2v1Zy0vvATcREwgDOIYEZzWiC/HIfGXv/zFzeWAQ3GbSki8qwPmimNAID09Pu6urNzuXd+xE41HzfzUaK6URGix04ZJ4rw+Mo2IN1sh8eabb7b7LZ+defvtt9XQ0GB7oWBPnDdiLsRE0FgHEXe2QmLcuHHav39/2vefMGGCKioqbC8UYBc7cwSBdS7eL2RwjK1rJCzL0l133aVevXqldf/6+npHCwX74nitRCpcjAk/EBDAF2yFxLe+9S39/e9/T/v+EyZMUM+ePW0vFJAJLsaEVwiI1phGQLIZEvy2T7MxlWiNoICbiAggNd61gcjjAIBMcP1NakwjkERIRAwbd2ocDGAX6wyQHscfSAWEEac70BXioWu8YEFLjicSe/bscXM54CI28q7xahNtsU6kh/0L2nIcEieffLLuuusuHTrEhofw4uAB1gEgM45DoqysTFu2bNGoUaO0du1aN5cJLuBVgz0cTOKHf3P72K8gFcchMXHiRP31r3/Vf/3Xf+nuu+/WuHHjeHsoQi95cOEAE038+wLuy/hdG9OnT9c777yjqVOnasqUKfrOd76jf/7zn24sGzLEq4fMcMCJDv4tM8f+BB1x5e2flmWpuLhY//Ef/6GNGzdqzJgx+tGPfqQDBw648fDIABt/5ngVG178uwHec/z2z4cfflivvPKKXnnlFe3atUtZWVkaO3asZs+erdNPP13r1q3TKaecoieffFLjx493c5mBwPD2UfMRDu7jBQk64zgk7rnnHp111lmaMWOGzjrrLI0fP145OTnNf/7DH/5QS5cu1TXXXKO//e1vriwsnOGjs93X8mBFVASPePAOEYGuOA6JdD5HYubMmbrrrruc/gggFIiKYBAPgBk8/WTLAQMG6M9//rOXPwJpYirhD6LCW8SDv5hGIB2ehkQikdDZZ5/t5Y+ADcSEv9oe9AgL+wgHwHz8rg3AJ4RF1wgHczCNQLoIiZhhKmEOwoJwMBURATsICcAQHR1UoxAYBAMQXYREDDGVCJfODsKmRQbBEH5MI2AXIRFTxEQ0OD1wdxYgxEB8ERFwgpAAYohYAOAWV37XBsKJVx8AktgfwClCIubYeQBgP4BMEBJgJwIAcIyQAIAY44UEMkVIQBI7EyCO2O7hBkICzdipAADsIiQAIIZ44QC3GBcSK1euVGFhoXJzc1VUVKRt27al9X0vvPCCsrOzdfrpp3u7gBHHzgWIPrbz+CktLVUikdDcuXNdf2yjQmLDhg2aO3euFixYoPLyck2ePFkXX3yxKioqOv2+mpoaTZ8+Xeeff75PSxpt7GSA6GL7jp9XXnlFq1ev1tixYz15fKNCYvny5Zo5c6auu+46jR49Wvfdd5+GDRumVatWdfp9N9xwg6666ipNmDDBpyWNPnY2QPSwXcfPwYMH9f3vf1+PPPKIjjvuOE9+hjEhUV9frx07dqi4uLjV7cXFxdq+fXuH37d27Vq9++67WrhwYVo/p66uTrW1ta3+Q2rsdADAPG2PYXV1dR3ed/bs2ZoyZYouuOACz5bHmN+1sW/fPjU2NqqgoKDV7QUFBaqqqkr5Pf/4xz80b948bdu2TdnZ6f1VSktLtXjx4oyXFwDChBcGwepVLWX1yOwxGuuP/e+wYcNa3b5w4UItWrSo3f1/+9vf6rXXXtMrr7yS2Q/ugjEhkZRIJFp9bVlWu9skqbGxUVdddZUWL16sk046Ke3Hnz9/vkpKSpq/rq2tbfePgi/wW0KB8CMiomXPnj3Ky8tr/jonJyflfW655RZt2bJFubm5ni6PMSHRv39/ZWVltZs+VFdXt5tSSNKBAwf06quvqry8XDfddJMkqampSZZlKTs7W1u2bNF5553X7vtycnJSPunoGDEBhBcRET15eXmtQiKVHTt2qLq6WkVFRc23NTY26rnnntOKFStUV1enrKwsV5bHmJDo0aOHioqKVFZWpu985zvNt5eVlenb3/52u/vn5eXprbfeanXbypUr9ec//1lPPPGECgsLPV/mOCEmgPAhIuLr/PPPb3eMvPbaa3XyySfr9ttvdy0iJINCQpJKSkp09dVXa/z48ZowYYJWr16tiooKzZo1S9Kx0xJ79+7VY489pm7dumnMmDGtvn/AgAHKzc1tdzvcQUwAQDj07du33bGwd+/e6tevn+vHSKNCYtq0adq/f7+WLFmiyspKjRkzRps2bdKIESMkSZWVlV1+pgS8RUwA4cA0An5JWJZlBb0QQaqtrVV+fr7OP26GsrtleEltTBASgNmIiPQ1NNXrmU9/rZqami6vO3AieYwZ8x9LldUjs4seG+uP6m+r7/BsWZ0y5nMkEB7spABzsX3Cb4QEHGFnBZiH7RJBICTgGDstwBxsjwgKIYGMsPMCgsd2iCAREsgYOzEgOGx/CBohAVewMwP8x3YHExAScA07NcA/bG8wBSEBV7FzA7zHdgaTEBJwHTs5wDtsXzANIQFPsLMD3Md2BRMREvAMOz3APWxPMBUhAU8ldu9lBwhkiG0IJiMk4At2hIAzbDswHSEB37BDBNLHNA9hQUjAV+wYga6xnSBMCAn4jp0k0DG2D4QNIYFAMLYF2mObQBgREggUO06AsEa4ERIIHDtQxBnrP8KOkIAReEWGOGKdRxQQEjAKO1bEAeGMKCEkYBx2sIgy1m9ETXbQCwCkktzZWoVDAl4SwB0EBKKKiQSMxs4XUcB6jChjIgHjMZ1AWBEQiAMmEggNdsoIE9ZXxAUTCYQK0wmYjoBA3BASCCWCAqYhIBBXnNpAqPF+fJiAdRBxRkggEtiRIwiELMCpDUQIpzvgF+IB+AIhgcghKOAVAgJoj5BAZBEUcAsBgUz1rmxUdvfGjB6j4fPMvt8rhAQir+VBgKiAHQQE0DVCArHClAJdIR4AewgJxBJBgbYICMAZQgKxxmkPEBBAZggJ4H8xpYgP4gFwDyEBtMGUIpqIB8AbhATQCaIi3IgHwHuEBJAmoiIciAfAX4QE4EDbgxVhESziAQgOIQG4gGmFvwgHwByEBOAyphXuIxwAcxESgMdSHQSJi44RDUC4EBJAADo6WMYpMAgGIBoICcAgnR1cwxgZxAIQfYQEEBJ2DspeRgdxAKAlQgKIIA72APzSLegFAAAA4UVIAAAAxwgJAADgGCEBAAAcIyQAAIBjhAQAAHCMkAAAAI4REgAAwDFCAgAAOEZIAAAAxwgJAADgGCEBAAAcIyQAAIBjhAQAAHCMkAAAAI4REgAAwDFCAgAAOGZcSKxcuVKFhYXKzc1VUVGRtm3b1uF9//CHP+jCCy/Ul7/8ZeXl5WnChAnavHmzj0sLAICZSktL9Y1vfEN9+/bVgAEDdNlll+nvf/+76z/HqJDYsGGD5s6dqwULFqi8vFyTJ0/WxRdfrIqKipT3f+6553ThhRdq06ZN2rFjh84991xNnTpV5eXlPi85AABm2bp1q2bPnq2XXnpJZWVlamhoUHFxsQ4dOuTqz0lYlmW5+ogZOPPMM/X1r39dq1atar5t9OjRuuyyy1RaWprWY5x66qmaNm2a7r777rTuX1tbq/z8fJ1/3Axld+vhaLkBAOHU0FSvZz79tWpqapSXl+f64yePMWdO/T/K7p6b0WM1fH5Uf/3vuxwv67/+9S8NGDBAW7du1be+9a2MlqWlbNceKUP19fXasWOH5s2b1+r24uJibd++Pa3HaGpq0oEDB3T88cd3eJ+6ujrV1dU1f11bW+tsgQEACEDb41ZOTo5ycnK6/L6amhpJ6vQY6YQxIbFv3z41NjaqoKCg1e0FBQWqqqpK6zF+/vOf69ChQ7ryyis7vE9paakWL16c0bICAGBH7w8PKTurMaPHaGg8KkkaNmxYq9sXLlyoRYsWdfq9lmWppKRE3/zmNzVmzJiMlqMtY0IiKZFItPrasqx2t6Xy+OOPa9GiRfrTn/6kAQMGdHi/+fPnq6SkpPnr2tradv8oAACYas+ePa1ObaQzjbjpppv05ptv6vnnn3d9eYwJif79+ysrK6vd9KG6urrdlKKtDRs2aObMmfrd736nCy64oNP7pjsCAgDARHl5ebaukbj55pu1ceNGPffccxo6dKjry2PMuzZ69OihoqIilZWVtbq9rKxMEydO7PD7Hn/8cV1zzTVav369pkyZ4vViAgAQCpZl6aabbtIf/vAH/fnPf1ZhYaEnP8eYiYQklZSU6Oqrr9b48eM1YcIErV69WhUVFZo1a5akY6cl9u7dq8cee0zSsYiYPn267r//fp111lnN04yePXsqPz8/sL8HAABBmz17ttavX68//elP6tu3b/MxMj8/Xz179nTt5xgVEtOmTdP+/fu1ZMkSVVZWasyYMdq0aZNGjBghSaqsrGz1mRK/+MUv1NDQoNmzZ2v27NnNt8+YMUO/+tWv/F58AACMkfwohXPOOafV7WvXrtU111zj2s8xKiQk6cYbb9SNN96Y8s/axsGzzz7r/QIBABBCfn1MlDHXSAAAgPAhJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADHCAkAAOAYIQEAABwjJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADHCAkAAOAYIQEAABwjJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADHCAkAAOAYIQEAABwjJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADHCAkAAOBYdtALAABA1CXer1SiW4/MHqOp3qWlcRcTCQAA4BghAQAAHCMkAACAY4QEAABwjJAAAACOERIAAMAxQgIAADhGSAAAAMcICQAA4BghAQAAHCMkAACAY4QEAABwjJAAAACOERIAAMAxQgIAADhGSAAAAMcICQAA4BghAQAAHCMkAACAY4QEAABwjJAAAACOERIAAMAxQgIAADhGSAAAAMcICQAA4BghAQAAHCMkAACAY4QEAABwjJAAAACOERIAAMAxQgIAADhGSAAAAMcICQAA4JhxIbFy5UoVFhYqNzdXRUVF2rZtW6f337p1q4qKipSbm6uRI0fq4Ycf9mlJAQAwn93jql1GhcSGDRs0d+5cLViwQOXl5Zo8ebIuvvhiVVRUpLz/7t27dckll2jy5MkqLy/XHXfcoTlz5uj3v/+9z0sOAIB57B5XnUhYlmW59mgZOvPMM/X1r39dq1atar5t9OjRuuyyy1RaWtru/rfffrs2btyoXbt2Nd82a9YsvfHGG3rxxRfT+pm1tbXKz8/X+cfNUHa3Hpn/JQAAodHQVK9nPv21ampqlJeX5/rju3mMcbKsdo+rTmS78iguqK+v144dOzRv3rxWtxcXF2v79u0pv+fFF19UcXFxq9suuugirVmzRp9//rm6d+/e7nvq6upUV1fX/HVNTY0kqcGql5oy/VsAAMKkwaqXJHn9mtqNY0xyWWtra1vdnpOTo5ycnHb3d3JcdcKYkNi3b58aGxtVUFDQ6vaCggJVVVWl/J6qqqqU929oaNC+ffs0aNCgdt9TWlqqxYsXt7t962ePZ7D0AIAw279/v/Lz811/3B49emjgwIHaWuXOMaZPnz4aNmxYq9sWLlyoRYsWtbuvk+OqE8aERFIikWj1tWVZ7W7r6v6pbk+aP3++SkpKmr/+7LPPNGLECFVUVHiyEkVRbW2thg0bpj179ngyCowinjP7eM7s4zmzr6amRsOHD9fxxx/vyePn5uZq9+7dqq+vd+XxUh0TU00jWrJ7XLXLmJDo37+/srKy2lVSdXV1u5pKGjhwYMr7Z2dnq1+/fim/p6MRUH5+PhueTXl5eTxnNvGc2cdzZh/PmX3dunn33oPc3Fzl5uZ69vgdcXJcdcKYd2306NFDRUVFKisra3V7WVmZJk6cmPJ7JkyY0O7+W7Zs0fjx41NeHwEAQFw4Oa46YUxISFJJSYl++ctf6tFHH9WuXbt06623qqKiQrNmzZJ07LTE9OnTm+8/a9YsffDBByopKdGuXbv06KOPas2aNfrxj38c1F8BAABjdHVcdYMxpzYkadq0adq/f7+WLFmiyspKjRkzRps2bdKIESMkSZWVla3e+1pYWKhNmzbp1ltv1UMPPaTBgwfrgQce0He/+920f2ZOTo4WLlzY5TkmfIHnzD6eM/t4zuzjObMv6s9ZV8dVNxj1ORIAACBcjDq1AQAAwoWQAAAAjhESAADAMUICAAA4FouQ4FeT22fnOfvDH/6gCy+8UF/+8peVl5enCRMmaPPmzT4urRmc/qreF154QdnZ2Tr99NO9XUAD2X3O6urqtGDBAo0YMUI5OTk68cQT9eijj/q0tGaw+5ytW7dOp512mnr16qVBgwbp2muv1f79+31a2mA999xzmjp1qgYPHqxEIqE//vGPXX4P+38HrIj77W9/a3Xv3t165JFHrJ07d1q33HKL1bt3b+uDDz5Ief/33nvP6tWrl3XLLbdYO3futB555BGre/fu1hNPPOHzkgfH7nN2yy23WMuWLbNefvll65133rHmz59vde/e3Xrttdd8XvLg2H3Okj777DNr5MiRVnFxsXXaaaf5s7CGcPKcXXrppdaZZ55plZWVWbt377b++te/Wi+88IKPSx0su8/Ztm3brG7duln333+/9d5771nbtm2zTj31VOuyyy7zecmDsWnTJmvBggXW73//e0uS9eSTT3Z6f/b/zkQ+JM444wxr1qxZrW47+eSTrXnz5qW8/09+8hPr5JNPbnXbDTfcYJ111lmeLaNp7D5nqZxyyinW4sWL3V40Yzl9zqZNm2bdeeed1sKFC2MXEnafs6eeesrKz8+39u/f78fiGcnuc/azn/3MGjlyZKvbHnjgAWvo0KGeLaOp0gkJ9v/ORPrURvJXqLb9VeNOfjX5q6++qs8//9yzZTWFk+esraamJh04cMCzX4JjGqfP2dq1a/Xuu+9q4cKFXi+icZw8Zxs3btT48eN17733asiQITrppJP04x//WEeOHPFjkQPn5DmbOHGiPvzwQ23atEmWZenjjz/WE088oSlTpvixyKET9/2/U0Z9sqXb/PrV5FHixq+d/fnPf65Dhw7pyiuv9GIRjePkOfvHP/6hefPmadu2bcrOjvRmmJKT5+y9997T888/r9zcXD355JPat2+fbrzxRn3yySexuE7CyXM2ceJErVu3TtOmTdPRo0fV0NCgSy+9VA8++KAfixw6cd//OxXpiUSS17+aPIqc/trZxx9/XIsWLdKGDRs0YMAArxbPSOk+Z42Njbrqqqu0ePFinXTSSX4tnpHsrGdNTU1KJBJat26dzjjjDF1yySVavny5fvWrX8VmKiHZe8527typOXPm6O6779aOHTv09NNPa/fu3a7+noWoYf9vX6RfCvn1q8mjJJNfO7thwwbNnDlTv/vd73TBBRd4uZhGsfucHThwQK+++qrKy8t10003STp2kLQsS9nZ2dqyZYvOO+88X5Y9KE7Ws0GDBmnIkCHKz89vvm306NGyLEsffvihRo0a5ekyB83Jc1ZaWqpJkybptttukySNHTtWvXv31uTJk/XTn/6UV9htxH3/71SkJxL8anL7nP7a2ccff1zXXHON1q9fH7vzr3afs7y8PL311lt6/fXXm/+bNWuWvvrVr+r111/XmWee6deiB8bJejZp0iR99NFHOnjwYPNt77zzjrp166ahQ4d6urwmcPKcHT58WN26td7NZ2VlSfrilTa+EPf9v2MBXeTpm+TbpdasWWPt3LnTmjt3rtW7d2/r/ffftyzLsubNm2ddffXVzfdPvv3n1ltvtXbu3GmtWbMmdm//sfucrV+/3srOzrYeeughq7Kysvm/zz77LKi/gu/sPmdtxfFdG3afswMHDlhDhw61Lr/8cuvtt9+2tm7dao0aNcq67rrrgvor+M7uc7Z27VorOzvbWrlypfXuu+9azz//vDV+/HjrjDPOCOqv4KsDBw5Y5eXlVnl5uSXJWr58uVVeXt78dln2/+6IfEhYlmU99NBD1ogRI6wePXpYX//6162tW7c2/9mMGTOss88+u9X9n332WWvcuHFWjx49rBNOOMFatWqVz0scPDvP2dlnn21JavffjBkz/F/wANldz1qKY0hYlv3nbNeuXdYFF1xg9ezZ0xo6dKhVUlJiHT582OelDpbd5+yBBx6wTjnlFKtnz57WoEGDrO9///vWhx9+6PNSB+Mvf/lLp/sm9v/u4NeIAwAAxyJ9jQQAAPAWIQEAABwjJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADHCAkAAOAYIQHEyDnnnKNEIqFEIqHXX389o8e65pprmh/rj3/8oyvLByB8CAkgZq6//npVVlZqzJgxGT3O/fffr8rKSpeWCkBYZQe9AAD81atXLw0cODDjx8nPz1d+fr4LSwQgzJhIACH2+OOPKzc3V3v37m2+7brrrtPYsWNVU1OT9uOcc845uvnmmzV37lwdd9xxKigo0OrVq3Xo0CFde+216tu3r0488UQ99dRTXvw1AIQYIQGE2Pe+9z199atfVWlpqSRp8eLF2rx5s5566inb04Jf//rX6t+/v15++WXdfPPN+s///E9dccUVmjhxol577TVddNFFuvrqq3X48GEv/ioAQoqQAEIskUjonnvu0S9/+UstXbpU999/v55++mkNGTLE9mOddtppuvPOOzVq1CjNnz9fPXv2VP/+/XX99ddr1KhRuvvuu7V//369+eabHvxNAIQV10gAIfdv//ZvOuWUU7R48WJt2bJFp556qqPHGTt2bPP/z8rKUr9+/fS1r32t+baCggJJUnV1dWYLDCBSmEgAIbd582b9z//8jxobG5sP9k5079691deJRKLVbYlEQpLU1NTk+GcAiB5CAgix1157TVdccYV+8Ytf6KKLLtJdd90V9CIBiBlObQAh9f7772vKlCmaN2+err76ap1yyin6xje+oR07dqioqCjoxQMQE0wkgBD65JNPdPHFF+vSSy/VHXfcIUkqKirS1KlTtWDBgoCXDkCcMJEAQuj444/Xrl272t3+pz/9ydHjPfvss+1ue//999vdZlmWo8cHEF1MJICYWblypfr06aO33noro8eZNWuW+vTp49JSAQirhMVLDCA29u7dqyNHjkiShg8frh49ejh+rOrqatXW1kqSBg0apN69e7uyjADChZAAAACOcWoDAAA4RkgAAADHCAkAAOAYIQEAABwjJAAAgGOEBAAAcIyQAAAAjhESAADAMUICAAA4RkgAAADH/j+8OlSQk2UxzAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 620x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAHBCAYAAADJgdkTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1P0lEQVR4nO3dd3QUVR/G8e9m00MSCBAgkIReQ+9VioJUK6gogg1REHhRmggiFhAUQQkINiyoWAABUYoivUNoobdQAiEB0uvuvH8sRJGoIW1Tns85OceZvTv7mxGyD3Pv3GsyDMNAREREJI842LsAERERKVoUPkRERCRPKXyIiIhInlL4EBERkTyl8CEiIiJ5SuFDRERE8pTCh4iIiOQphQ8RERHJUwofIiIikqcUPkSKuPnz52MymdJ/HB0dKVeuHA8//DDHjh2zS02hoaFMnDiR06dP3/LagAEDqFixYp7XJCI5R+FDRAD47LPP2LJlC2vWrGHIkCEsXbqUNm3acPXq1TyvJTQ0lNdeey3D8DF+/HgWL16c5zWJSM5xtHcBIpI/BAUF0aRJEwDat2+PxWLh1VdfZcmSJTzxxBN2ru5PVapUsXcJIpJNuvMhIhm6EUQuXbqUvq99+/a0b9/+lrZ/7wo5ffo0JpOJd955h+nTp1OpUiWKFStGy5Yt2bp1679+7vz58+nduzcAHTp0SO8Omj9/foafBWAymRgyZAifffYZNWrUwM3NjSZNmrB161YMw2DatGnpNXTs2JHjx4/f8rlr1qyhU6dOeHl54e7uTuvWrfntt98ycaVE5HYpfIhIhk6dOgVA9erVs3yM4OBgVq9ezYwZM1iwYAHx8fF069aN6Ojof3xP9+7deeutt9Lfv2XLFrZs2UL37t3/9bOWL1/Oxx9/zJQpU/jmm2+IjY2le/fuvPjii2zatIlZs2Yxb948QkNDeeCBB/jrgt5fffUVnTt3xsvLi88//5zvvvsOHx8funTpogAikgvU7SIiAFgsFtLS0khKSmLTpk288cYbtGvXjl69emX5mJ6enixfvhyz2QyAn58fzZo145dffuHhhx/O8D2lS5emWrVqANSuXZsWLVpk6rOSk5NZtWoVHh4egO1uyL333svatWvZvXs3JpMJgMuXLzN8+HAOHDhA3bp1SUhIYNiwYfTo0eOmsSTdunWjUaNGvPzyy2zbti3L10BEbqU7HyICQIsWLXBycsLT05O7776bEiVK8NNPP+HomPV/o3Tv3j09eADUq1cPgDNnzmS73r/r0KFDevAAqFWrFgBdu3ZNDx5/3X+jhs2bN3PlyhX69+9PWlpa+o/VauXuu+9mx44dxMfH53i9IkWZ7nyICABffPEFtWrVIjY2loULFzJ37lweeeQRfvnllywfs2TJkjdtu7i4AJCYmJitWjPi4+Nz07azs/O/7k9KSgL+HNPy4IMP/uOxr1y5clOwEZHsUfgQEcB2R+DGINMOHTpgsVj4+OOP+eGHH9K/mF1dXTMcrxEZGZmnteakUqVKAfDBBx/8YxdPmTJl8rIkkUJP3S4ikqGpU6dSokQJJkyYgNVqBaBixYocPXqU5OTk9HZRUVFs3rw5Rz87N++Q/F3r1q0pXrw4oaGhNGnSJMOfG3dLRCRn6M6HiGSoRIkSjB07llGjRvH111/z2GOP0a9fP+bOnctjjz3GM888Q1RUFFOnTsXLyytHPzsoKAiAefPm4enpiaurK5UqVbqlGycnFCtWjA8++ID+/ftz5coVHnzwQXx9fbl8+TJ79+7l8uXLzJkzJ8c/V6Qo050PEflHL7zwAgEBAUyaNAmLxULr1q35/PPPOXjwIPfccw9vvPEGY8eOzXDuj+yoVKkSM2bMYO/evbRv356mTZuybNmyHP2Mv3rsscdYu3YtcXFxPPvss9x5550MGzaM3bt306lTp1z7XJGiymT89WF3ERERkVymOx8iIiKSpxQ+REREJE8pfIiIiEieUvgQERGRPKXwISIiInlK4UNERETyVKGbZMxqtXLhwgU8PT1vWkxKREREco9hGMTGxuLn54eDw7/f2yh04ePChQv4+/vbuwwREZEi6ezZs1SoUOFf2xS68OHp6QnYTj6np3wWERGRjMXExODv75/+PfxvCl34uNHV4uXlpfAhIiKSxzIz5EEDTkVERCRPKXyIiIhInlL4EBERkTxVaMZ8BAcHExwcjMViyVR7i8VCampqLleV/zk7O//nI1EiIiI5yWQYhmHvInJSTEwM3t7eREdHZzjg1DAMLl68yLVr1/K+uHzIwcGBSpUq4ezsbO9SRESkAPuv79+/KjR3PjLrRvDw9fXF3d29SE9EdmNCtvDwcAICAor0tRARkbxTpMKHxWJJDx4lS5a0dzn5QunSpblw4QJpaWk4OTnZuxwRESkCilRn/40xHu7u7nauJP+40d2S2bEyIiIi2VWkwscN6l74k66FiIjktSIZPkRERMR+FD5EREQkTyl8iIiISJ5S+BAREZE8pfBRQFitVt5++22qVq2Ki4sLAQEBvPnmmwBs3ryZBg0a4OrqSpMmTViyZAkmk4mQkBD7Fi0iIvmOxWph2YllXIy/aLcaitQ8H39nGAaJqfZ5xNTNyXxbT5qMHTuWjz76iPfee482bdoQHh7O4cOHiY2NpWfPnnTr1o2vv/6aM2fOMHz48NwrXERECiTDMFh7di0f7PmA49eO80C1B5jYaqJdainS4SMx1ULtCSvt8tmhk7rg7py5yx8bG8vMmTOZNWsW/fv3B6BKlSq0adOGDz/8EJPJxEcffYSrqyu1a9fm/PnzPPPMM7lZvoiIFCDbw7czc89M9l3eB4CnsycVvSrarZ4iHT4KikOHDpGcnEynTp1uee3IkSPUq1cPV1fX9H3NmjXLy/JERCSfOhh5kJm7Z7IlfAsAbo5uPFbrMQYEDcDL+d/XX8lNRTp8uDmZCZ3UxW6fnem2bm7/+JphGLd03xSytQJFROQ2nbx2klkhs1h9ZjUAjg6O9K7em4H1BlLKrZSdqyvi4cNkMmW668OeqlWrhpubG7/99htPP/30Ta/VrFmTBQsWkJycjIuLCwA7d+60R5kiImJnF+IuMGfvHJaeWIrVsGLCRM8qPXmu/nNU8Kxg7/LS5f9vXsHV1ZXRo0czatQonJ2dad26NZcvX+bgwYP07duXcePGMXDgQMaMGUNYWBjvvPMOoKnTRUSKiqjEKD7a/xHfHfmOVKttHbOO/h15oeELVC1R1c7V3SrfhY+zZ8/Sr18/IiIicHR0ZPz48fTu3dveZdnd+PHjcXR0ZMKECVy4cIFy5coxaNAgvLy8WLZsGc899xwNGjSgbt26TJgwgb59+940DkRERAqf2JRY5h+cz5ehX5KYlghA87LNGdpoKPVK17Nzdf/MZOSzAQLh4eFcunSJBg0aEBERQaNGjThy5AgeHh6Zen9MTAze3t5ER0fj5XXzYJqkpCROnTpFpUqVCvUX84IFC3jiiSeIjo7+1/EiUHSuiYhIYZKUlsQ3h7/hkwOfEJ0cDUCdknUY1mgYLf1a2qWmf/v+/bt8d+ejXLlylCtXDgBfX198fHy4cuVKpsNHUfTFF19QuXJlypcvz969exk9ejR9+vT5z+AhIiIFS6o1lcXHFjN371wiEiMAqOxdmRcavkCngE4Fprs9x2c4Xb9+PT179sTPzw+TycSSJUtuaTN79uz0f2k3btyYDRs2ZHisnTt3YrVa8ff3z+kyC5WLFy/y2GOPUatWLf73v//Ru3dv5s2bZ++yREQkh1gNKytOruDeJffy+tbXiUiMoJxHOV5v/TqLei3izsA7C0zwgFy48xEfH0/9+vV54okneOCBB255feHChQwfPpzZs2fTunVr5s6dS9euXQkNDSUgICC9XVRUFI8//jgff/xxTpdY6IwaNYpRo0bZuwwREclhhmGw4fwG3t/9PkeuHgHAx9WHgfUG0rt6b5zNzlk7cFoyGAY42ae7PcfDR9euXenates/vj59+nSeeuqp9EdGZ8yYwcqVK5kzZw6TJ08GIDk5mfvuu4+xY8fSqlWrf/285ORkkpOT07djYmJy4CxERETsa9elXby/+312R+wGoJhTMQbUGUC/2v1wd3LP2kFjLsCOT2DXfLhjNDQfmHMF34Y8HfORkpLCrl27GDNmzE37O3fuzObNmwFbyhswYAAdO3akX79+/3nMyZMn89prr+VKvSIiInnt8JXDzNw9k43nNwLgYnahb82+PBn0JMVdi9/+AQ0Dzm6DbR/CoWVgTbv+QcuKRviIjIzEYrFQpkyZm/aXKVOGixdtq+tt2rSJhQsXUq9evfTxIl9++SV169bN8Jhjx45lxIgR6dsxMTEaIyIiIgXOmZgzzNozi19P/wqAo8mR+6rdx7P1nqWMR5n/eHcGUpPgwI+20HFx35/7A1tD82ehRvccqvz22eVpl4ymA7+xr02bNlit1kwfy8XFJX1mTxERkYLmYvxFPtz7IUuOL8Fi2FZa71qpK0MaDCHAK+A/3p2B6POw83rXSkKUbZ+jK9TtbQsdZTP+x3xeytPwUapUKcxmc/pdjhsiIiJuuRsiIiJSmF1Nuson+z/hm8PfkGJNAaBdhXYMbTiUGj41bu9ghgFhW2DbXFvXyvUQg7c/NH0KGvUHd5/05lFxyZhMJnw8sjhgNZvyNHw4OzvTuHFjVq9ezX333Ze+f/Xq1dxzzz3ZOnZwcDDBwcFYLJbslikiIpJr4lPj+SL0Cz4/+DnxqfEANPJtxPDGw2no2/D2DpaaCPt/gO1z4eL+P/dXbGu7y1G9K5htX/UxSamsPHCRZfvC2XQ8khc6VmX4ndVz6rRuS46Hj7i4OI4fP56+ferUKUJCQvDx8SEgIIARI0bQr18/mjRpQsuWLZk3bx5hYWEMGjQoW587ePBgBg8enD7DWmHTvn17GjRowIwZM+xdioiIZEGyJZnvjnzHR/s+4mryVQBq+tRkaMOhtCnf5vbm6bh29nrXyueQeMW2z9EN6vWBZgOhbBAAiSkWfjt4gaUhF/jjyGVSLH8OazgWEZdj53a7cjx87Ny5kw4dOqRv3xgM2r9/f+bPn89DDz1EVFQUkyZNIjw8nKCgIFasWEFgYGBOlyIiImJ3adY0lp1Yxuy9s7kYbxt2EOgVyJAGQ+hcsTMOpkzO92kYcGaTrWvl8HIwrgcJ7wBo9jQ07AfuPqSkWVkfeomley+w5tAlElL+7BGo6luMXvX96FGvHJVLF8vpU820HA8f7du357+Wi3n++ed5/vnnc/Rz1e0iIiL5iWEYrD6zmg/2fMDpmNMA+Lr78lz957in6j04OThl7kApCbD/e9g+Dy4d+HN/pXbQ7Fmo0ZU0w8TWk1dYtncfvxwIJyYpLb1ZhRJu9KzvR6/6ftQs65kvZkLNd2u7ZFVh73b5K5PJxOLFi7n33nvT9xUvXpwZM2YwYMAAu9UlIiK20LHlwhZm7plJaFQoAMVdivN03ad5uObDuJgz+YTmtTDY8THs/gISbd00OLpB/Yeh2UCspWuxO+wqy5Yd4uf94UTGpaS/1dfThe71ytGzvh8N/Yvni8DxV4UmfGSJYUBqgn0+28kd8tkfBhERyZ6QiBBm7p7Jzks7AXB3dKd/nf48XvtxijlnopvDMOD0RtvcHEdW/Nm1UjwAmg3EaPAYB686sGzXBZbvW8v5a4npby3u7kTXoHL0rF+O5pVKYnbIv98xRTt8pCbAW372+eyXL4CzVuoVESkMjl49yge7P+CPc38A4OzgzEM1H+Lpuk/j4+rz728GW9fKvoW2rpWI0D/3V24PzZ7lePHWLN1/ieWz93IyMj79ZQ9nM13qlKVnfT9aVy2Fs2OOrxebK4p2+BAREcmGs7FnCQ4JZsXJFRgYOJgcuLfqvTxX/znKepT97wNcPQM7PoLdX0LSNds+J3eo/wjhNfqx+Lwny1aGcyh8Y/pbXBwd6FTLl571/OhQ0xdXJ3PunFwuKjThI0sDTp3cbXcg7CGriwJhG/Px90G9qamp2a1IREQy6XLCZebum8uPR38kzbAN7uwc2JkhDYdQybvSv7/ZMODUettTK0d/+bNrpURFYus9wVJTB34IjWPPxj+/nxwdTLSrXpqe9ctxZ60yeLpmcrBqPlVowkeWBpyaTAWy66N06dKEh4enbx87doyEBDuNXRERKUKik6P59MCnfH3oa5IsSQC09mvNC41eoE7JOv/+5pR4W9fKtnlw+VD67tTA9mwu9SDzwiuzZdU1rMY5wPYV1aJSSXo18OPuOmUpYafZSHNDoQkfRUnHjh2ZNWsWLVq0wGq1Mnr0aJycCnYKFhHJzxJSE1hwaAGfHfiM2NRYAOqXrs+wRsNoWrbpv7/5yinbUyt7voSkaAAMJw9Ole/JJyl3sfCYO2lHDOAaAA0DitOznm0uDl8v11w8K/tR+CiA3n33XZ544gnatWuHn58fM2fOZNeuXfYuS0Sk0Em1pPL90e+Zt28eUUm2RdqqlajG0IZDuaPCHf/8CKthwMk/rnet/ArYusrjPQJY4daTKRcbE3X4RrAwqFXOi571y9Gznh/+Plnvli8oFD4KiD/++CP9v/38/Fi5cuVNr1+7di1vCxIRKcQsVgs/n/qZ2SGzOR93HoAKxSowuOFgulXq9s+zkibHwb5vYftHcPlw+u7DHs14P64jv0QFYWB7b6VSHtcn/ypHVV/PXD+n/KTQhA/NcCoiItllGAa/n/2dWXtmcfyabZ2y0m6lebbes9xf7X6czP/QxX3lJGz/GPZ8Bcm2rpVkBzcWWe/go+Q7OZlkm9bBz9uVnvX96Fnfjzp+Xvlu8q+8UmjCR1Ga4VRERHLetvBtvL/7ffZF7gPAy9mLJ4OepG+tvrg5ut36BsOAE7/D9nkYR1diut61EkY5Pk29ix8s7YjDnVLFnOlf1zbbaKOAEjjk48m/8kqhCR8iIiJZcSDyADN3z2Rr+FYA3BzdeKzWYwwIGoCXs9etb0iOhb3fYmyfhynyKAAmYK2lPvMtd7PeWhdPV2e6N7QFjhaVfXA0F4zJv/KKwoeIiBRJJ6+d5IM9H7AmbA0Ajg6O9K7em4H1BlLKrdStb4g6Ads/wrrnKxxSYjEBsYYbP1ja8YWlMxcdK3BX3TJ8VN+PttVL4eJY8Cb/yisKHyIiUqScjzvPnJA5LDu5DKthxYSJnlV68lz956jgWeHmxlYrnPidpE2zcT39GwAOwAlrOT63dGYZ7WlaI4AXG/jRsaYv7s76Ws0MXSURESkSIhMj+WjfR3x39DvSrLZZSTsFdGJIgyFULVH15sZJMcRt/xLL1rl4J5zhxkOxv1ka8qW1C0blDvSoX56XgsriVcBnG7WHQhM+9LSLiIhkJDYlls8OfMZXh74iMc22Cmzzcs0Z1nAYdUvXvbnt+UNcXP0+Fc4soZhhmzk6xnDje0t7Qsr1plnjJrwbVJaSxVzy/DwKk0ITPvS0i4iI/FViWiLfHP6GT/Z/QkxKDABBJYMY1ngYLcq1SG+XkJzCvj8W4RHyMXUTd3Bjxo3jVj/WeN6DW5NH6daoKk95Z/DEi2RJoQkfIiIiAKnWVBYfW8zcvXOJSIwAoLJ3ZYY2HErHgI6YTCaS0yxsOnCSK5s+o0nEIlqYbOtlWQ0T25yaElHrceq3u49BpYvZ81QKLYWPQu706dNUqlSJPXv20KBBA3uXIyKSa6yGlV9O/UJwSDBnY88C4Ofhx/MNnqdH5R4YhokNxyLZun0L/se/oqexjmKmJDBBLO4cKXcvPu2fp2WNuv/xSZJdCh8iIlKgGYbB+nPreX/P+xy9apt3w8fVh4H1BvJA1QfZdy6e1346SPT+FTyQupxR5v22N5rgsmtFkhs9Tfk7BtDEpWhNcW5PCh8iIlJg7by4k/f3vM+eiD0AFHMqxoA6A2jo3ZPVB6/RY/FK2sWv4knzKio6XAIzWDERXaET3u1foHSVO2xr10ueUvgoIH744Qdee+01jh8/jru7Ow0bNuSnn37Czc2NN954g3nz5nH58mVq1arFlClTuPvuu+1dsohIrjkUdYiZe2ay6fwmAFzMLnQNeBC3hDv5ZlUs311ZzgDzSn4yb8DDKRmAVCcvHBo/jrnZ05TwqWTP8ou8QhM+svKorWEY6Y9d5TU3R7dMLygUHh7OI488wtSpU7nvvvuIjY1lw4YNGIbBzJkzeffdd5k7dy4NGzbk008/pVevXhw8eJBq1arl8lmIiOSt09GnmRUyi5WnbSt7m01manjcyZVz7VhwwEQHhxVMMq+incv+9PdYS9fEofmzONV7CJw97FW6/IXJMAzD3kXkpBuP2kZHR+PldfOc/ElJSZw6dYpKlSrh6upKQmoCzb9ubpc6t/XdhruTe6ba7t69m8aNG3P69GkCAwNveq18+fIMHjyYl19+OX1fs2bNaNq0KcHBwf854PTv10REJD+6GH+RD/d+yJLjS7AYFsBEsdQmXDxzB56pbvQ2/0F/x9UEmGxPtxgmB0w1ukGzgVCpnbpW8sC/ff/+XaG581GY1a9fn06dOlG3bl26dOlC586defDBBzGbzVy4cIHWrVvf1L5169bs3bvXTtWKiOScq0lX+Xj/x3x7+FtSrCkApMXVJDmiC2VSLLzh+BMPuG3E1bB1reBaHBo9jqnp01Ai8J8PLHZVpMOHm6Mb2/pus9tnZ5bZbGb16tVs3ryZVatW8cEHHzBu3DhWr14NcEv3jWEYme7SERHJj+JT4/lo72d8eegLUqy27vG0+EqkXu5Mh+QoXij2JfVMIbbGBuBbB5oPhLp9wDlzd5XFfop0+DCZTJnu+rA3k8lE69atad26NRMmTCAwMJDffvsNPz8/Nm7cSLt27dLbbt68mWbNmtmxWhGRrLEaVubs/JZPQ2eTQjQAliQ/zBF38KJzBH1cgvEyLkAKYHKAmt2h2bNQsY26VgqQIh0+Copt27bx22+/0blzZ3x9fdm2bVv6ky0jR47k1VdfpUqVKjRo0IDPPvuMkJAQFixYYO+yRUQyzTAMPt+9jtn73iXR4TQA1pSSVEtoyWsuF6hreh+HxOsPCLiVgEb9oelTUDzAfkVLlil8FABeXl6sX7+eGTNmEBMTQ2BgIO+++y5du3alS5cuxMTE8OKLLxIREUHt2rVZunSpnnQRkQIhJc3Kgp37mLPvfRJddoADGBYXWhsNmWwNw+fq7D8blwmC5s9C3d7gpHVWCrIi/bSL6JqIiH1EJ6by5dbjfLp/PimeazA5pGAYJhqaqjM17hzlLh+yNTSZbV0rzQdBYCt1reRjetpFRETypfPXEvl0w0kWhv6C4bMMB+8rmICKpjK8GXeNepdtA+lxLgZNnrQ9Klvc3641S85T+BARkVx34Hw0H204yc+H9+DkuxTHsicwASVMHoyKS6R7xA5MAC7etq6VFs+Bu4+dq5bcUmjCR1ZmOBURkdxjGAbrjl5m3vqTbD51FpfSq3GtuBWTycAJM08kpvHUpcO4Gwa4+UDL5213Oly97V265LJCEz4GDx7M4MGD0/ucRETEPpLTLCwNucDHG05x5NI1nEpsp1jV1ZjMCQDclWxlRMR5KqRZwMMXWr1g62JxKWbnyiWvFJrwcTsK2RjbbNG1EJGcEp2YyoJtZ5i/6TQRscmY3Y9TrPJyTC4XAaiWZmX05cs0T0oGTz9oMxwaPa4nV4qgIhU+nJycAEhISMDNTX/YAVJSbNMVm81mO1ciIgXVuasJfLrxNAt3hBGfYsHkdIXiFX/F4rYPAG+rwQtXrvJAbByOxQPgzhHQoC84uti5crGXIhU+zGYzxYsXJyLCtvCQu7t7kZ6G3Gq1cvnyZdzd3XF0LFJ/FEQkB+w/F828DSdZsT8ci9UAUzLlKm4myf130oxUzIZBn5g4Bl+LxrtEZbj3RdscHWYne5cudlbkvnHKli0LkB5AijoHBwcCAgKKdAgTkcyzWv8cRLrlZNT1vQZ1qh8nxnUR11KvggHNE5MYHXWVasWrwn3vQJ37wEF3WMWmyIUPk8lEuXLl8PX1JTU11d7l2J2zszMODg72LkNE8rnkNAs/7bnARxtOciwiDgBHBxPt6iYR47yAI3FHIRXKp6Yx8spVOnpVxXT/dKjRHfQ7Rv6myIWPG8xms8Y5iIj8h+iEVL7adob5m09zOda2bH0xF0fua+pFotM3rLy4FlLAzWpl4LUY+nlWx+XemVDtLs1GKv+oyIYPERH5Z2evJPDJxlN8t/MsCSm2+ZPKebvyeKsKmJyWMP/YN8Rj298zNp5hHtUoc89YqHSHQof8J4UPERFJt+/cNeauP8kv+8OxXn8Sv1Y5L55pWxFvx7VMDxlNmGG7AxKUnMwY1yrU7znetu6KSCYpfIiIFHFWq8HaIxHMW3+SbaeupO9vV700A9tWxo9tTN3Sh02GbaxHqTQLw10C6Nn1dRz8m9qrbCnAFD5ERIqopFQLS/ac56MNJzlxOR6wDSLt1cCPZ9pWpnziHuZsvJ9vrVdIM5lwMgz6OZVhYMe38PBvbufqpSBT+BARKWKuxqfYZiLdfIbIOFsXiqeLI31bBDCgVUV8r+xg0apePGuJ4KrZDCYT7c0lGNn6NQIqdbBz9VIYFJrwoYXlRET+XVhUAp9sPMl3O8+RmGr7Xenn7cqTbSrxUJMKeJ7fwK4fBzM07TyHXZzBbKaygzujmoykda0H7Vy9FCYmo5At7nFjYbno6Gi8vLzsXY6IiN2FnL3GR+tP8suBPweR1i7nxbN3VKZbUFmcTqwifP0Upqec5ddiHgB4mhx5vs6TPNRwEE4OmpFU/tvtfP8WmjsfIiLyJ6vV4LfDEXy0/iTbT/85iPSO6qV5tl1lWlYugenQMhI/msbHyWF86u1FkrMHJuDBit0Z0nwUPq4+9jsBKdQUPkRECpGkVAuLrw8iPXl9EKmT2cQ9DcrzTNvK1CjtBgcXYcyexq+J55juU5xw9+IANC5VjzEtx1PTp6Ydz0CKAoUPEZFC4Gp8Cl9uPcMXW04TGWdbrdrT1ZHHWgQyoFVFyniYYe+38N10DsedZYpPCXb5lgKgnHsZXmw6ks6BnbXOk+QJhQ8RkQLsTFR8+kykSalWAMoXd7MNIm3qTzGHNAj5CjbO4ErseT4o4c2PfmUxTCZczS48WfcpBtQZgJujm53PRIoShQ8RkQJod9hVPlp/kl8PXuTGYwNB5b0Y2K4K3YLK4mhJgl3zYPP7pMaG862XJ3P8yxPrYLuz0bViV0Y0GUFZj7J2PAspqhQ+REQKCKvVYM2hS3y04SQ7Tl9N39+hRmmeaVeZlpVLYkqJgy0zYfMsSIhkk5srb/v7c8rRFjpq+dRidLPRNC7T2F6nIaLwISKS3yWlWvhx9zk+2XCKk5G2QaTOZgfubejH020rU72MJyReg3VTYetsSLrGGUdH3ikfwB/OtmP4uPowtOFQ7q16L2YHregt9qXwISKST12JT+HLLbZBpFHxtkGkXn8ZROrr5QrxUfDbJNj+ESTHEGcyMc+vIl+6QpphxdHkSN9afXm2/rN4OWvuI8kfFD5ERPKZU5HxfLLxJD/sOpc+iLRCCTeealOJPk388XBxhNhLsPJ12PkppCZgBX7yq8ZMD0ei0uLBgNblWzOq6Sgqe1e27wmJ/I3Ch4hIPrHrzFXmrT/BqtBL6YNI61XwZmC7ytxdpyyOZgeIPge/vQ+7P4e0JAD2+tVhio83B+LPQVoygV6BjGo6inYV2tnxbET+mcKHiIgdWawGq0Ntg0h3nflzEGmnmr48064yzSv52ObeuHIKNs2APQvAmgrApQqNmFEukOWRuyA+Fg8nDwbVG8SjtR7Fyawp0SX/UvgQEbGDxBQLP+w+xycbTnI6KgGwDSK9v1F5nm5biaq+nraGkcdgw3TYtxAM22JwyRXb8EXF+nx0fg2JkbswYeLeqvcytNFQSrmVstcpiWSawoeISB6Kikvmiy1n+HLrGa5cH0Tq7eZEvxaBPN4qEF9PV1vDS6Gw4R04sAiw9cEYlTvye+27mHZmKefPLAOgQekGjGk2hjql6tjjdESyROFDRCQPnLwcx8cbT/HjrnMkp9kGkfr7uPF0m8r0blIBd+frv44v7IH178Dh5X++uUY3jjV8mLfDlrHt4BwAfN19GdF4BN0qddOU6FLgKHyIiOQSwzCuDyI9yepDfw4ire9fnGfbVaZLnbKYr884ytntsH4aHFt1/d0mqH0P0S2eZdaFtXy3fQJWw4qzgzMDggbwVNBTuDu52+W8RLIrX4aP++67jz/++INOnTrxww8/2LscEZHbYrEarDp4kXkbTrIn7Fr6/jtrlWFgu8o0rVjCdrfCMODUBlvoOLXO1sjkAHV7k9ZqGD9c28+szaOITo4G4K7AuxjReAQVPCvY4axEck6+DB9Dhw7lySef5PPPP7d3KSIimZaYYuGHXWf5eOMpztwYROrowAONKvBUm0pU9S1ma2gYcPw3W+gI22Lb5+AI9R+GNiPYlhLJlG2vcPzacQCqlajG6KajaV6uuT1OSyTH5cvw0aFDB/744w97lyEikimXY5P5cstpvtx6hqsJtsdgi7s78XiLQPq1rEhpTxdbQ8OAo7/aQsf5XbZ9Zmdo9Di0HsY5swPv7nyXNWFrAPB28eaFBi/wQPUHcHTIl7+uRbIkx/80r1+/nmnTprFr1y7Cw8NZvHgx9957701tZs+ezbRp0wgPD6dOnTrMmDGDtm3b5nQpIiK56sTlOD7ecJIfd58n5fog0sCS7jzdphIPNP7LIFKrFQ79BOvfhUv7bfsc3aDJk9DqBRLcvPl4/8d8fvBzUqwpmE1m+tTow+AGg/F28bbT2YnknhwPH/Hx8dSvX58nnniCBx544JbXFy5cyPDhw5k9ezatW7dm7ty5dO3aldDQUAICAnK6HBGRHBcZl8z01Uf5dnsY1uuDSBtcH0Ta+a+DSC1pcHCR7emVyCO2fc7FoOnT0HIIhkcpfj71M+/tfI+IxAgAmpdrzuimo6lWopodzkwkb+R4+OjatStdu3b9x9enT5/OU089xdNPPw3AjBkzWLlyJXPmzGHy5Mm3/XnJyckkJyenb8fExNx+0SIimZCcZuHzzaf54LfjxCanAXBnLV8G3VGFxoEl/nzkNS3FNinYxulw5aRtn4s3tBgEzQeBuw8HIw8yef3/2Ht5LwDli5VnZNORdPTvqEdnpdDL007ElJQUdu3axZgxY27a37lzZzZv3pylY06ePJnXXnstJ8oTEcmQYRisCr3EWysOpQ8kDSrvxYQedWhWyefPhqlJEPIVbJwB0Wdt+9x8oOVgaPYMuHoTmRjJzE3j+en4TxgYuDm6MbDeQPrV7oeL2SXvT07EDvI0fERGRmKxWChTpsxN+8uUKcPFixfTt7t06cLu3buJj4+nQoUKLF68mKZNm2Z4zLFjxzJixIj07ZiYGPz9/XPnBESkyAm9EMPry0PZcjIKgNKeLozqUoMHGlXA4Ub3SkoC7JoPm9+H2HDbPg9faD3UNq7D2YNUSypfHfiMufvmEp8aD0DPyj0Z1mgYZTzKZPDJIoWXXYZP//2WomEYN+1buXJlpo/l4uKCi4v+tSAiOetybDLTVx/h2x1nMQzbI7MD21bmufZVbEvaAyTHwo6PYfMsSIi07fMqD62HQ6N+4OSGYRisP7uOaTuncSbmDABBJYMY03wM9UvXt8/JidhZnoaPUqVKYTabb7rLARAREXHL3RAREXtITrPw2abTzPr9OHHXx3X0qFeOMV1rUqHEX2YUPbYGFj/7Z+goHghtR0D9vuDoDMDJ6JNM3TGVTec3AVDKrRTDGw2nZ5WeOJgc8vS8RPKTPA0fzs7ONG7cmNWrV3Pfffel71+9ejX33HNPto4dHBxMcHAwFoslu2WKSBFkGAYrD17krRWHCbtiG9dRr4I3E3rUpknFv4zrsFrgj8m2J1gwoGRVaPsS1O0NZtuv1JiUGD7c+yHfHPqGNCMNJwcn+tXux8B6A/Fw8rDD2YnkLzkePuLi4jh+/Hj69qlTpwgJCcHHx4eAgABGjBhBv379aNKkCS1btmTevHmEhYUxaNCgbH3u4MGDGTx4MDExMXh767l4Ecm8A+ejeX15KNtOXQGgjJcLo7rU5L6G5f8c1wEQewl+fApOb7BtN30aOr8JTraVaC1WC4uPL+aDPR9wJcl2rPb+7RnZZCQBXppKQOSGHA8fO3fupEOHDunbNwaD9u/fn/nz5/PQQw8RFRXFpEmTCA8PJygoiBUrVhAYGJjTpYiI/KuI2CTeXXmU73bZxnW4ODrwbLvKPHvHX8Z13HBqgy14xF2yzdXRcybUfTD95V2XdvH29rc5dOUQAJW9KzOq6Shal2+dl6ckUiCYDOPGOouFw407H9HR0Xh5edm7HBHJh5JSLXy66RTBvx8nPsXWVXtPAz9G3V2T8sXdbm5stdrm61j7JhhW8K0Nfb6AUrZJwMLjwpm+azq/nv4VAE8nT55v8DwP1XwIJwenPD0vEXu6ne/fQrNYgMZ8iMh/MQyDFfsvMvmXQ5y7mgjYZiYd36M2jQNL3PqGhCuwaCAcX23bbvAodHsHnN1JtiTz6YFP+XT/pyRZknAwOfBgtQcZ3HAwPq4+tx5LRNLpzoeIFAn7z0UzaflBdpy+CkBZL1fGdK1Jr/p+N4/ruOHsdvj+CYg5B46u0P1daPgYAFsubOGNrW8QFhsGQOMyjRnTbAw1fWrm2fmI5DdF8s6HiEhGLsUkMW3lEX7cfQ7DAFcnBwbdUYWB7Sr/ufDbXxkGbJ0Dq8eDNc32NEvvz6FsEJGJkUzbMY0Vp1YA4Ovmy8imI+lSsYumRBe5DQofIlIoJaVa+HjDSWb/cYKE6+M67mtYnlF316Cct1vGb0q8Bj8NhsPLbdt17ode72N19uCHI98xY9cMYlNjcTA58EjNRxjSYAjFnIvlzQmJFCKFJnxozIeIgG1cx/J94Uz55TDnr9nGdTQMKM6EHrVpGJDBuI4bLoTA9/3h6mkwO0OXt6Dp0xy5epRJv01i3+V9ANQuWZsJLSdQp2Sd3D8ZkUJKYz5EpNDYe/Yak5aHsuuMbVyHn7crY7rVome9cv/cLWIYsOsz+GUMWJKheAD0/pwE35rMDpnNV4e+wmJY8HDy4IWGL/BwjYcxO5jz8KxECgaN+RCRIuVidBJTVx5m0e7zALg5mXm+fRWeblsZN+d/CQrJcbB8OOz/3rZdoxvcO5vfL+9h8k/3cDHethRE58DOjGo6SgvAieQQhQ8RKbASUyx8tOEkc/44QWKqrcv1gUYVGNmlBmW9Xf/9zRGH4LvHIfIomMxw12uE132AyVteZe3ZtQCUL1ael5u/TLsK7XL7VESKFIUPESlwDMNg6d4LvP3LYS5EJwHQJLAE43vUpr5/8f8+QMg3sPx/kJYInn6kPfAxC+KPE7z0XhLTEnE0OTIgaAAD6w3EzfEfBqeKSJYVmvChAaciRcOesKtMWh7KnrBrAJQv7sbYbjXpXvdfxnXckJoIK0bCni9t21U6sveO//F6yHscuXoEgEa+jRjfYjxVS1TNxbMQKdo04FRECoQL1xKZ+uthloRcAMDd2czgDlV5qk0lXJ0yMQA06oStm+XSATA5ENP2RWa6Wfn+6A8YGHi7ePNi4xe5p+o9Wu5eJAs04FRECo2ElDTmrjvJ3PUnSEq1YjLBg9fHdfh6/ce4jhsOLoafXoCUWAyP0vzS9jmmnllGVFIUAL2q9OLFJi9qWnSRPKLwISL5ktVq8NPe87z9yxEuxtjGdTSr6MP4HrWpW8E7cwdJS4ZVr8D2eQCEBTbjjbJ+bDkyH4BK3pUY32I8Tcs2zY1TEJF/oPAhIvnOrjO2cR17z14DoEIJN17uVouuQWUzP4351TPw/QC4sJsU4JN6d/Nx/HFSIi7iYnZhYL2BDKgzAGezc26dhoj8A4UPEck3zl9L5O1fDrN0r21ch4ezmcEdq/Jk60yO67jhyC+weBAkXWO7dyle9wvkdGwoAK38WvFK81fw9/LPjVMQkUwoNOFDT7uIFFzxyWl8uO4E89afJDnNNq6jT2N/XuxSHV/PTI7rALCkwu+vw6aZRDk48G5AdZaZkyDpMqXcSjG66WgtAieSD+hpFxGxG6vVYNGe80z99TARsckANK9kG9cRVD6T4zpuiLkAPzyJNWwLizw9eK90WWKMVEyY6FOjD0MbDcXLWb8TRHKLnnYRkXxvx+krvL48lH3nogEI8HHn5W616FKnzO3fmTjxO/z4DMdSo3ndz489Lo5gpFLTpyYTWkygbum6uXAGIpJVCh8ikqfOXklgyq+H+XlfOADFXBx5oWNVBrSuiIvjbS7YZrXAuqkkrJ/Kh8W9+NK3HGkmcHN0Y0iDIfSt1RdHB/2aE8lv9LdSRPJEXHIac/44zkcbTpGSZsXBBA81DWDEXdUp7emShQNehkVPsz58K2+WL8sFJ9uvs04BnRjTbAxlPcrm8BmISE5R+BCRXGW1Gvyw+xzTVh7h8vVxHS0rl2R8j9rU9sviGIwzm7n44xO87ZrGmrK+AJTzKMfLzV+mvX/7HKpcRHKLwoeI5JptJ6N4/edQDpyPASCwpDvjutXirtpZGNcBYLWStmkG3+x8j1klvEhwcMZscuDxOv0ZVG8Q7k7uOXwGIpIbFD5EJMedvZLA5F8OsWL/RQA8XRwZ2qkaj7cKvP1xHTckXOHAogFMSjjCIZ/iANQvVZcJrSZSvUT1HKpcRPJCoQkfmudDxP5ik1IJXnuCTzeeIsViG9fxSDPbuI6SxbIwruPGcU9v4P2Vz7HQBQwXZzwdXBjRbDT3V39Ai8CJFECa50NEss1iNfh+51neWXWEyLgUANpULcUrPWpRs2zW/x4aVisrfxvF1LAVXL5+x6RHuda81PZNSrqVzJHaRSRnaJ4PEckzW05EMWl5KIfCbeM6KpXy4JXutehY0zdbM4mejTzEm788zSZrDDiaCTS58Modb9MisFNOlS4idqLwISJZciYqnrdWHGLlwUsAeLk6MuzO6vRrEYizY9a7QlItqczfNoW5R78j2QROhsEzpZvzZJdgXBxvY6p1Ecm3FD5E5LbEJKUS/PtxPtt0mhSLFbODiUebBzD8zur4eGRvhdgd4Tt4Y91LnEy+AiZonmrwSrspVKzeI4eqF5H8QOFDRDLFYjVYuOMs7646QlS8bVxH22qlGN+jNtXLeGbr2FeTrvLu9rf56dTPAPhYLIx08qf7I19j8tDYDpHCRuFDRP7T5uORTFoeyuGLsQBULu3B+O61aV+jdLbGdRiGwZLjS3h3x1SiU+MA6B0bz7C6A/FuOwoc9CSLSGGk8CEi/+hUpG1cx+pQ27gObzcnht9ZjcdaBOJkzl4wOHHtBJO2TGJ3xG4AqienMD7BRIN7v4SKrbNdu4jkXwofInKL6MRUZv1+jPmbT5NqMTA7mOjXIpBhnapRIpvjOhLTEvlo30d8dvAz0qxpuFmtPH81mkdLNsSp76dQrHQOnYWI5FeFJnxokjGR7EuzWPl2x1mmrz7KlevjOtrXKM0r3WtR1Td74zoANp7fyBtb3+B83HnbseMTeDnqGuXavAR3jAKHLM5+KiIFiiYZExEANhy7zBvLD3Hkkm1cR1XfYozrXosONXyzfeyIhAje3v42q86sAqCMxcrYyCg64Q73fwRVNXeHSEGnScZEJNNOXI7jrZ8P8dvhCACKuzsx4q7qPNIsINvjOixWC98e+ZYP9nxAfGo8Zkw8Gh3N4KvRuPu3gAc/BS+/nDgNESlAFD5EiqjohFRm/naML7acJs1q4Ohg4vGWFRnWqRre7k7ZPn5oVCiTtkziYNRBAOpaHZkQfpaaKanQehh0HA/m7H+OiBQ8Ch8iRUyaxcrX28OYvvoo1xJSAehU05eXu9eiSuli2T5+XEocs0Jm8c3hb7AaVjzNrgy7co0Hr0RgdvWGhz+Emt2y/TkiUnApfIgUIeuOXuaN5aEci7DNqVG9TDFe6V6bdtWz/4SJYRisCVvDlG1TiEi0deF0dfNn1JEtlLJYwa8h9P4cSgRm+7NEpGBT+BApAo5HxPHmz6GsPXIZgBLuTozoXINHmvrjmM1xHQDn487z5tY32XB+AwD+Hn68EptKq9BNtgbNnoXOr4OjS7Y/S0QKPoUPkUIsJimV6auO8uXWM1iuj+sY0KoiL3Sqhrdb9sdbpFpT+eLgF3y490OSLEk4OjjyVPk7eXrXYlzjL4OzJ/R6H4Luz4GzEZHCQuFDpJA6cD6a5xfsJuxKAgB31irDy91qUjkHxnUA7L60m9e3vs7xa8cBaFqmCa+Y/ai8YTYYVigTZOtmKVU1Rz5PRAoPhQ+RQsYwDL7aeobXlx8ixWKlQgk3ptxfjzbVSuXI8a8lXeO93e+x6NgiAEq4lOCles/Sc+cPmE7a9tGwH3SbBk5uOfKZIlK4KHyIFCJxyWmM+XEfy/eFA7a7He/2rp8jj84ahsGyk8t4Z8c7XE2+CsAD1R5geJk2FF/yAsReAEc36DEdGvTN9ueJSOGl8CFSSBwKj2Hwgt2cjIzH0cHE6Ltr8nTbStladfaGk9EneWPrG+y4uAOAqsWrMr75KzQ6uRm+6gOGBUpVt3WzlKmd7c8TkcJN4UOkgDMMg4U7zvLq0oMkp1kp5+3KrL6NaBxYItvHTkpL4uP9H/PJgU9Is6bhanZlUP1BPF6pJ07LhsORn20Ngx6EnjPBJWfGk4hI4abwIVKAJaSk8criAyzac32hthqlmd6nAT7ZXHkWYPP5zbyx7Q3Oxp4FoG35trzc/GUqxETAx3fCtTNgdoa7p0CTJyEH7rCISNFQaMKHVrWVoubYpVieX7CbYxFxOJjgpS41GNSuCg4O2QsBkYmRTN0+lV9O/wKAr5svY5qP4U7/Tph2fgIrXwZLCpSoaOtm8WuQ/ZMRkSJFq9qKFEA/7jrHK0sOkJhqwdfThQ8eaUjzyiWzdUyL1cL3R7/n/d3vE5sai4PJgb41+zK4wWCKGQYsGwYHfrQ1rtkD7gkGt+LZPxkRKRS0qq1IIZWUauHVnw6ycKetK6RN1VLMeLgBpYplb+bQw1cOM2nLJPZH7gegTsk6jG85njol68Clg/Dd4xB1HBwc4a5J0OJ5dbOISJYpfIgUECcuxzF4wW4OX4zFZILhnaozpGNVzNnoZklITWBWyCwWHFqA1bDi4eTB0IZDeajGQ5gdzLDnK/j5JUhLBK/y0Hs++DfLuZMSkSJJ4UOkAFi69wJjf9xHfIqFUsWcmflwQ1pXzd6kYb+F/cbkbZO5lHAJgM6BnRndbDS+7r6QkgArXoKQBbbGVe+E++aBR/a6dkREQOFDJF9LSrXwxs+hfLU1DIDmlXz44JGG+Hq5ZvmYF+IuMHn7ZP44+wcA5YuVZ1zzcbSt0NbWIPKYrZslIhRMDtBhHLQZAQ7ZX4BORAQUPkTyrTNR8Ty/YDcHL8QAMKRDVYbfWS3Lq9CmWlNZELqA2Xtnk5iWiKPJkQFBAxhYbyBujtenQT/wIywdCilx4OELD34Cldrl1CmJiAAKHyL50q8Hwhn5/T5ik9Mo4e7Eew81oH0N3ywfLyQihElbJ3Hs6jEAGvk2YnyL8VQtcX3Rt7Rk2yO0Oz62bVdsCw98Ap5lsnsqIiK3UPgQyUdS0qxM/uUQn206DUDjwBJ88EhD/IpnbYG22JRYpu+azg9HfwCguEtxRjQewT1V78HBdP0OypVT8P0ACA+xbbcbCXeMAbN+PYhI7tBvF5F84uyVBIZ8s4e9Z68B8Gy7yrzUpQZOWexm2X95PyPXj+R8nG3203uq3MOLTV6khOtfpl0//DMsfg6So8HNB+7/CKrdmd1TERH5VwofIvnA6tBLvPhdCDFJaXi7OfFu7/rcWTtrXR5Ww8qXoV8yY9cM0ow0yhcrz+utX6dp2aZ/NrKkwpqJsGWWbbtCM+j9GXhXyP7JiIj8B4UPETtKtVh5Z+UR5q4/CUB9/+IE921IhRLuWTre1aSrvLLpFdafWw/AXYF3MbHVRLyc/zLbYPQ5+P4JOLfdtt1yCNw5EcxO2TkVEZFMU/gQsZPw6ESGfL2HXWeuAvBE64qM7VoLZ8esdbPsvLiT0RtGE5EQgbODM6ObjaZ39d6Y/joT6fE18OMzkHgFXLzh3mCo1TMnTkdEJNMUPkTs4I8jEfxvYQhXE1LxdHFk6oP16Fq3XJaOZbFa+Gj/R8zZOwerYaWiV0XeueMdavjU+LOR1QJ/TIb17wAGlKtvm63Up3KOnI+IyO1Q+BDJQ2kWK++tOUrw2hMABJX3IrhvIwJLemTpeJcTLjN2w1i2XdwGQK8qvRjXfBzuTn/ptom9BD8+Bac32LabPAVd3gKnrE9UJiKSHQofInkkIiaJF77Zw7ZTVwB4rEUAr3SvjauTOUvH23R+Ey9vfJkrSVdwc3RjXPNx3FP1npsbnd4IPzwJcZfAyQN6zoR6vbN7KiIi2ZIv50tevnw5NWrUoFq1anz88cf2Lkck2zYdj6Tb+xvYduoKHs5m3n+kIW/cWzdLwSPVmsqMXTMYtGYQV5KuUL1Edb7t8e3NwcNqhY3vwec9bcGjdC0Y+IeCh4jkC/nuzkdaWhojRoxg7dq1eHl50ahRI+6//358fHzsXZrIbbNYDT74/RgzfzuGYUDNsp4EP9qIKqWLZel44XHhjFw/kr2X9wLQp3ofRjYdiavjX7pQEq/C4kFw9Ffbdv1HoPu74Jy1rh0RkZyW78LH9u3bqVOnDuXLlwegW7durFy5kkceecTOlYncnsuxyfxvYQgbj0cC8HBTfyb2qpPlbpbfw35n/KbxxKTEUMypGBNbTaRLxS43Nzq/G77vD9fCwOwC3aZBo8fhr0+8iIjYWY53u6xfv56ePXvi5+eHyWRiyZIlt7SZPXs2lSpVwtXVlcaNG7Nhw4b01y5cuJAePAAqVKjA+fPnc7pMkVy17WQU3d/fwMbjkbg5mXm3d32mPFAvS8EjxZLClO1TGLZ2GDEpMQSVDOK7nt/dHDwMA3Z8Ap92sQWPEhXh6dXQuL+Ch4jkOzkePuLj46lfvz6zZs3K8PWFCxcyfPhwxo0bx549e2jbti1du3YlLMy2ZLhhGLe8x6RfnlJAWK0GwWuP88hHW4mITaaqbzGWDmnNA42zNnNoWEwYj614jAWHFgDQv3Z/vuj6Bf6e/n82So6DRQPh5xFgSYEa3WHgOtvjtCIi+VCOd7t07dqVrl27/uPr06dP56mnnuLpp58GYMaMGaxcuZI5c+YwefJkypcvf9OdjnPnztG8efOcLlMkx12JT2HEdyH8ceQyAPc3LM8b9wXh7py1v2YrTq5g0tZJxKfGU9ylOG+2eZN2Ff62vP3lI/Dd43D5MJjMtplKW72gux0ikq/l6ZiPlJQUdu3axZgxY27a37lzZzZv3gxAs2bNOHDgAOfPn8fLy4sVK1YwYcKEfzxmcnIyycnJ6dsxMTG5U7zIv9h15gpDvt5DeHQSLo4OTLqnDn2a+Gfprl1iWiJvb3+bH4/9CEAj30a83e5tynqUvbnh/h9g6VBIjYdiZW1rswS2yonTERHJVXkaPiIjI7FYLJQpc/OCWWXKlOHixYu2ghwdeffdd+nQoQNWq5VRo0ZRsmTJfzzm5MmTee2113K1bpF/YhgGH284xdu/HibNalCplAezH21ErXJe//3mDBy/epyR60dy/NpxTJgYWG8gg+oPwtHhL39V05Jh5cuw4/pj6JXawQOfQDHfHDgjEZHcZ5enXf7+r0HDMG7a16tXL3r16pWpY40dO5YRI0akb8fExODv7/8v7xDJGdEJqbz4/V7WHLoEQM/6fky+vy7FXG7/r5VhGCw+vpjJ2yaTZEmilFspJredTItyLW5uePUMfD8ALuy2bbcbCe3HgkPWnqAREbGHPA0fpUqVwmw2p9/luCEiIuKWuyGZ5eLigouLS06UJ5Jpe89eY/DXuzl3NRFnswPje9bmseYBWepmiUuJY9LWSfxy6hcAWvm14s02b1LKrdTNDY+utA0sTboGbiXg/o+g2l05cDYiInkrT8OHs7MzjRs3ZvXq1dx3333p+1evXs0999zzL+/8b8HBwQQHB2OxWLJbpsg/MgyD+ZtP89aKQ6RaDAJ83Jn9aCOCyntn6XihUaGMXDeSsNgwzCYzQxoO4cmgJ3Ew/eVBNEsa/PEWbHjXtl2+MfT+HIrrDp+IFEw5Hj7i4uI4fvx4+vapU6cICQnBx8eHgIAARowYQb9+/WjSpAktW7Zk3rx5hIWFMWjQoGx97uDBgxk8eDAxMTF4e2fti0Dk38QkpTL6h338csB2565LnTJMfbA+3m5Ot30swzD4+vDXvLvzXVKtqZTzKMfUdlNp4Nvg5oZ/XxSu2bPQ+Q1wdM7m2YiI2E+Oh4+dO3fSoUOH9O0b4zH69+/P/Pnzeeihh4iKimLSpEmEh4cTFBTEihUrCAwMzOlSRHLMgfPRDP56N2eiEnAymxjbtRZPtK6YpW6W6ORoJmyawO9nfwegg38HXm/9Ot4ufwvNpzfBD0/Y1mZxLga93oegB3LidERE7MpkZDSrVwF2485HdHQ0Xl5Ze+JA5AbDMFiwLYxJy0NJSbNSvrgbwY82ooF/8SwdLyQihFHrRxEeH46TgxMvNnmRvjX73hxiDAM2zYTfJoFhsS0K1+cLKF09Z05KRCQX3M73b75b2yWrNOZDclpcchovL9rP0r0XAOhU05d3+9SnuPvtd3lYDSufHfiMD/Z8gMWw4O/pz7Q7plGnZJ2bGyZehSXPw5EVtu16D0OP6VoUTkQKFd35EMnA4YsxPP/Vbk5GxmN2MDGqSw2eaVsZB4fb72aJSoxi3MZxbLqwCYCuFbsyoeUEijn/bWXbCyG22UqvnQGzM3SdCo0HaLZSESkQiuSdD5Gc8t3Os0z46QBJqVbKerkyq29DmlT0ydKxtoVvY8yGMUQmRuJqdmVMszHcX+3+W7tZds2HX0aDJRmKB9q6Wfwa5Mj5iIjkNwofItclpKQxfslBftx9DoB21UvzXp/6lCx2+/PIpFnT+HDvh8zbNw8DgyreVZh2xzSqlah2c8OUeFj+P9i30LZdoxvcO9s2j4eISCFVaMKHxnxIdhyPiOX5Bbs5eikOBxO82LkGz91RJUvdLJfiLzF6w2h2XdoFwP3V7mdMszG4Obrd3PDy0euLwh26vijcq9BqqLpZRKTQ05gPKfIW7znHuMUHSEixUNrThfcfbkjLKv+8ntC/WX9uPeM2juNa8jXcHd2Z0HIC3St3v7Xh/h9g2TBIibMtCvfgp1CxdTbPRETEfjTmQyQTklItvLbsIN9sPwtAqyolmflwQ0p73n43S6ollZm7Z/J56OcA1PKpxbQ7phHo9bf5a9KSYdUrsH2ebbtiW9uicJ5ZW15ARKQgUviQIunk5TgGf72HQ+ExmEwwtGM1hnaqhjkL3SznYs8xav0o9kfuB6Bvzb682ORFnM1/eyT3WphtUbjztu4Y2r4I7V8Gs/4aikjRot96UuQs33eBMT/uJy45jZIezsx4uAFtq5XO0rFWnV7FxM0TiU2NxdPZk9dbv06ngE63Njy2GhY9Y5vHw7U43D8PqnfJ3omIiBRQhSZ8aMCp/JfkNAtv/nyIL7acAaBZRR8+6NuQMl6ut38sSzLTdkxj4RHbUyr1S9dnarup+BXzu7mh1QJr34IN79i2/RpBn8+heEC2zkVEpCDTgFMpEsKiEhj89W72n48G4Pn2VRhxV3UczQ7/8c5bnYw+ych1Izl69SgATwU9xeCGg3Fy+NsCc3ERtkXhTq23bTd9Brq8CY63P6ZERCS/04BTkb/49cBFRv6wl9ikNIq7O/FenwZ0qOmbpWMtPbGUN7a+QWJaIj6uPrzV5i1al8/gKZUzm+H7JyDuIjh52BaFq/tgNs9ERKRwUPiQQislzcrbvx7mk42nAGgUUJwP+jaifHG3/3jnrRJSE3hz25ssPbEUgOZlmzO57WRKu/9trIhhwOYPYM3E64vC1by+KFyN7J6OiEihofAhhdL5a4kMXrCbkLPXAHimbSVG3V0Tpyx0sxy5coSR60dyKvoUDiYHnqv/HM/UfQazg/nmhonX4KfBcHi5bbtuH+g5Q4vCiYj8jcKHFDq/H77EiO/2ci0hFS9XR97pXZ/Odcre9nEMw+D7o9/z9va3SbGm4Ovmy5R2U2hatumtjS+EwPf94erp64vCvQ2Nn9BspSIiGSg04UNPu0iqxco7q44wd91JAOpV8Ca4byP8fdxv+1ixKbFM3DyRVWdWAdC2fFveaPMGPq5/W2DOMGD357Bi1PVF4QKuLwrXMNvnIyJSWOlpFykULkYn8cI3u9lx+ioAA1pVZGy3mrg4mv/jnbc6EHmAl9a9xPm48ziaHBneeDj9avfDwfS3LpuUePj5Rdj7jW27ele4b44WhRORIklPu0iRsv7oZYYvDOFKfArFXBx5+4F6dK9X7raPYxgGX4R+wYzdM0izplG+WHmmtptKvdL1bm0cecy2KFxEKJgcoNP1ReEcbn9MiYhIUaPwIQWWxWowY81RZq09jmFA7XJezH60ERVL3f4Az6tJV3ll0yusP2ebk+OuwLuY2GoiXs4ZpPcDi2DpC9cXhStzfVG4Ntk9HRGRIkPhQwqkiNgkhn0TwpaTUQD0bR7AhB61cXW6/W6WXZd2MWr9KCISInB2cGZU01H0qdEH098Hi6alXF8Ubq5tW4vCiYhkicKHFDibT0Qy9JsQIuOScXc2M/n+utzToPxtH8ditfDx/o+ZvXc2VsNKRa+KvHPHO9TwyWBOjmtnry8Kt9O23WYEdBinReFERLJAvzmlwLBaDWatPc6MNUexGlCjjCfBjzaiqm+x2z7W5YTLjN0wlm0XtwHQq0ovxjUfh7tTBk/GHFsDi57+c1G4++ZCjbuzeTYiIkWXwocUCFFxyQxfGMKGY5EA9G5cgUn3BOHmfPvdLJvPb2bsxrFcSbqCm6Mb45qP456q99za0GqBPybD+ncAw/b4bO/PoURgNs9GRKRoKzThQ/N8FF7bT13hhW92cykmGVcnB16/J4jeTfxv+zip1lSC9wTzyYFPAKheojrT7phGZe/KtzaOu3x9Ubh1tu2mT0OXt7QonIhIDtA8H5JvWa0Gc9ef5J1VR7BYDaqU9mD2o42pUdbzto8VHhfOqPWjCLkcAkCf6n0Y2XQkro6utzY+swV+eAJiw22LwvWcCfV6Z/NsREQKN83zIQXe1fgUXvx+L78fjgDg3gZ+vHlfXTxcbv+P7Nqwtbyy6RViUmIo5lSMia0m0qVil1sbGgZsmQWrX7UtCleqhm22Ut+a2T0dERH5C4UPyXd2h11lyILdXIhOwtnRgdd61eHhpv63Pvr6H1IsKby36z2+OvQVAEElg5h6x1T8PTPosrllUbje0GMGuNz+YFYREfl3Ch+SbxiGwScbTzHll8OkWQ0qlnQn+NFG1PHzvu1jhcWE8dK6lzh05RAA/Wv3Z1ijYTiZnW5tHL7PNlvp1VO2ReHungJNntSicCIiuUThQ/KF6MRURn6/l1WhlwDoXrccUx6oi6drBmHhP/xy6hde2/Ia8anxeLt482brN7nD/45bGxoG7PkSfn7JtiicdwD0+RzKN8ru6YiIyL9Q+BC7O3A+mucW7OLslUSczQ680qMW/VoE3nY3S2JaIm9vf5sfj/0IQCPfRrzd7m3KepS9tXFKwvVF4b62bVe/G+6dA+4+t7YVEZEcpfAhdnU1PoV+n2zjakIqFUq4MfvRRtSrUPy2j3Pi2gleWvcSx68dx4SJZ+o9w3P1n8PRIYM/4pHHry8Kd9C2KFzH8dB6uBaFExHJIwofYlfvrj7C1YRUqpcpxvfPtsLb/fa6WQzDYMnxJby17S2SLEmUdC3JlHZTaFGuRcZvOLgYfnoBUmLBw9e2KFyltjlwJiIiklmFJnxokrGC5+CFaL7eFgbAa72Cbjt4xKfGM2nLJFacWgFAy3IteavtW5RyK3Vr47QUWD0Bts2xbQe2gQc/Ac8MumRERCRXaZIxsQvDMOgzdws7Tl+le71yBPe9vUGeoVGhjFw3krDYMMwmM0MaDuHJoCdxMGXQdRJ9zrYo3Lkdtu02/4MOr2hROBGRHKRJxiTf+ynkAjtOX8XNycy4brUy/T7DMPj68Ne8u/NdUq2plPUoy9R2U2no2zDjNxxfAz8+A4lXwNX7+qJwXXPoLEREJCsUPiTPxSWn8dYK2/wbgztUwa+4W6beF50czYRNE/j97O8AtPdvzxut38DbJYN5QKwWWPc2rJsKGFCuge0x2hIVc+YkREQkyxQ+JM998PsxImKTCfBx5+m2GSzqloGQiBBGrR9FeHw4Tg5OvNjkRfrW7Jvx47jxkbZF4U7+Ydtu8iR0mQxOGazjIiIieU7hQ/LUictxfLrxFAATetTG1cn8r+2thpXPDnzGB3s+wGJY8Pf0Z9od06hTsk7GbwjbCt8/AbEXwMn9+qJwfXL6NEREJBsUPiTPGIbBpGWhpFoM2tcoTadavv/aPioxinEbx7HpwiYAulbsyoSWEyjmnMF6K4YBW4JhzatgTYNS1aHPl1oUTkQkH1L4kDyz5lAE645exslsYkKP2v86g+n28O2M2TCGy4mXcTW7MqbZGO6vdn/G70mKhiXP/7koXNAD0PN9LQonIpJPKXxInkhKtfD68lAAnmpTmcqlMw4GFquFD/d9yNy9czEwqOJdhWl3TKNaiWoZHzh8H3zfH66cBAcnuHsyNH1ai8KJiORjCh+SJz5af5KwKwmU8XLhhY5VM2xzKf4SYzaMYeelnQDcX+1+xjQbg5vjPzwNs/tLWPESpCWBt//1ReEa59YpiIhIDlH4kFx3/loiwX8cB+DlbrXwcLn1j936c+t5ZeMrXE2+irujOxNaTqB75e4ZHzAlAVaMhJCvbNvVOtvm79CicCIiBYLCh+S6t34+RFKqlWYVfehV3++W1+eEzGH23tkA1PKpxbQ7phHoFZjxwaJO2BaFu3Tg+qJwr0Dr/2lROBGRAkThQ3LV5uOR/Lw/HAcTTOxV55YBo7+H/Z4ePPrW7MuLTV7E2eyc8cEOLoGfhlxfFK709UXh2uXyGYiISE5T+JBck2qxMnHZQQAebR5Ibb+b5/qPTIxk4uaJAPSv3Z+Xmr6U8YHSUmyP0G61hRQCWtmCh1e53CpdRERyUaEJH1rVNv/5cssZjl6Ko4S7Ey92rn7Ta4Zh8Mom2xiPGiVqMLTR0IwPEn3ONmnYue227dbDoeN4LQonIlKAFZqO8sGDBxMaGsqOHTvsXYoAkXHJvLfmKAAvdalBcfebu1K+Pvw1m85vwsXswtvt3s64q+X4bzC3nS14uHjDw9/AXa8peIiIFHD6LS65Yuqvh4lNSiOovBcPNw246bXjV48zfed0AEY0HkGV4lVufrPVYlsQbt3b2BaFqw+9PwefSnlUvYiI5CaFD8lxIWev8d3OcwC81qsOZoc/B5mmWFIYvWE0KdYU2pRvwyM1H7n5zfGRsOgZOGFbuZbGT8DdU7QonIhIIaLwITnKajV49acDANzfsDyNA2+ee2Pm7pkcvXqUEi4leL316zc//RK2Db4f8OeicD3eg/oP52H1IiKSFxQ+JEf9sOsce89FU8zFkTFdb17UbcuFLXwR+gUAk1pPopRbKdsLhgFb58Dq8bZF4UpWg4e+BN9aeV2+iIjkAYUPyTHRiam8/ethAIZ2qoqv159dJdeSrvHKxlcA6FO9D+3929teSE2CJc/BwUW27Tr3Q6/3wcUzL0sXEZE8pPAhOWbGmqNExadQubQHA1r9OTjUMAwmbZ1ERGIEFb0q/jmfR+I1+PZROLNRi8KJiBQhCh+SI45cjOWLLWcAmNizDs6Ofz7FveT4ElafWY2jyZEp7abYFoqLCYevHoCIg+DiBQ8v0GylIiJFhMKHZJthGExcehCL1aBz7TK0q146/bWwmDAmb58MwOCGg6lTsg5cPgpf3Q/RZ6FYGXjsRyhb117li4hIHlP4kGxbsf8iW05G4eLowPgetdP3p1pTGbthLIlpiTQu05gn6jwBZ3fA130g8QqUrAqPLYIS/7CInIiIFEoKH5ItCSlpvPlzKADP3lEFfx/39Nfm7ZvHvsh9eDp5MrnNZMzH18B3/SEtEco3hr7fgUcpe5UuIiJ2ovAh2TLnjxNciE6ifHE3nrvjz5lKQyJCmLdvHgDjW46n3LHfYOlQMCxQ9S7o8zk4e9irbBERsaNCs7aL5L2wqATmrj8JwCvda+HmbAYgLiWOMRvGYDWs9Kjcg65nD8JPg23Bo35feOQbBQ8RkSJMdz4kyyYtDyUlzUrrqiW5O6hs+v7J2ydzPu48fh7leDk6EXa+bnuhzf+g06t6lFZEpIhT+JAs+eNIBGsOXcLRwcTEnnXSp0n/9fSvLD2xFAeTA5NTPfE88Blgsq3P0mKQfYsWEZF8QeFDbltKmpVJy2yDTPu3qki1MrbZSC/GX2TSlkkAPG31pNGRNbbJw+6fC0EP2K1eERHJX/LlmI/77ruPEiVK8OCDD9q7FMnAp5tOcTIynlLFXBh2ZzUArIaVcRvHEZsSS12LA4NO7wdnT9scHgoeIiLyF/kyfAwdOpQvvvjC3mVIBi7FJPHBb8cAGH13DbxcnQD4/ODnbL+4HTcDJl84h1OxMvDECqh8hz3LFRGRfChfho8OHTrg6amFxfKjySsOEZ9ioWFAcR5oVAGA0KhQ3t89E4DRkVEEegXCU6ugXD17lioiIvnUbYeP9evX07NnT/z8/DCZTCxZsuSWNrNnz6ZSpUq4urrSuHFjNmzYkBO1ip3tOH2FJSEXMJngtV51cHAwkZiWyJjfXiDNsNAxPoH7varbgkeJivYuV0RE8qnbDh/x8fHUr1+fWbNmZfj6woULGT58OOPGjWPPnj20bduWrl27EhYWlt6mcePGBAUF3fJz4cKFrJ+J5CqL1eDVnw4C8FATf+pVKA7Au78O4lRiBKXT0phYLAhT/+WatVRERP7VbT/t0rVrV7p27fqPr0+fPp2nnnqKp59+GoAZM2awcuVK5syZw+TJtgXGdu3alcVyb5WcnExycnL6dkxMTI4dW/709fYwQsNj8HJ1ZGSXGmAYrF81goVRuwF4w7MeJR74CsxOdq5URETyuxwd85GSksKuXbvo3LnzTfs7d+7M5s2bc/Kj0k2ePBlvb+/0H39//1z5nKLsanwK7646AsCIu6pT0t2JyBX/Y/z5lQA85lGFVr2/VfAQEZFMydHwERkZicVioUyZMjftL1OmDBcvXsz0cbp06ULv3r1ZsWIFFSpUYMeOHf/YduzYsURHR6f/nD17Nsv1S8beWXWEawmp1CzryWNNymL88AQTwpZzxWymmrMPw+9bqFlLRUQk03JlkjHT376IDMO4Zd+/WblyZabburi44OLikun2cnsOnI/m6+228Tqv3x2A47d9+DZyNxtK+eBsMjPl7o9wMev6i4hI5uXonY9SpUphNptvucsRERFxy92QnBYcHEzt2rVp2rRprn5OUWIYBhOXHsQwoG9tZ5qufYyT5zbzTskSAPyvyUtUL1HdzlWKiEhBk6Phw9nZmcaNG7N69eqb9q9evZpWrVrl5EfdYvDgwYSGhv5rF43cniUh59l55iq1nC4xKfJ/pF7az5gyZUk2mWjl14q+tfrau0QRESmAbrvbJS4ujuPHj6dvnzp1ipCQEHx8fAgICGDEiBH069ePJk2a0LJlS+bNm0dYWBiDBmlRsYIkLjmNySsOU890goWu7+IYc43pfoEccjIo7lKcN1q/gYMpX85RJyIi+dxth4+dO3fSoUOH9O0RI0YA0L9/f+bPn89DDz1EVFQUkyZNIjw8nKCgIFasWEFgYGDOVS257oPfjlErfjsfuszALTWZ7eXrMN85DoCJrSZS2r20nSsUEZGCymQYhmHvInJCcHAwwcHBWCwWjh49SnR0NF5eXvYuq0A6cTmOOTPfZLJ5Lk4mC9GV7+AB13guJUbwQLUHmNhqor1LFBGRfCYmJgZvb+9Mff8WmvBxw+2cvNzKsFpZ+P5IHr72sW07qDcjS3qxMmw1gV6BfNfjO9yd3O1cpYiI5De38/2rTnv5k9VK2DfD04PHtYaDWFa/JyvDVuNocmRK2ykKHiIikm25Ms+HFEBpyVgWDSLw2CIAfg8cRvWOA3lr2YMAPNfgOYJKBdmzQhERKSR050MgKQYW9MYcuogUw8wE83AaP/wyYzeMJT41nka+jXgq6Cl7VykiIoVEoQkfmmQsi2IvwfzucGod8YYrT6aOonHPgXx99DNCLodQzKkYb7V9C7OD2d6ViohIIaEBp0VZ1An48j64doYYcwn6JryIe8UmvHyfGwN+HYDFsDC57WR6VO5h70pFRCSfu53vX435KKrO74YFvSEhksRiAfSI+h/nKMP3XSvy8sansBgWulbqquAhIiI5TuGjKDr+GyzsB6nxGGXr83jcCMIMJx5vGchPZ+dwNvYs5TzK8UqLV+xdqYiIFEKFZsyHZNLehfB1H0iNh8rt+bLmbHZEOlHC3YlGtc6y+PhiTJh4s82beDmr20pERHJeoQkfGnCaCZveh8UDwZoGQQ9yuddXTFt7HoBBHUsxbdcbADwZ9CRNy+o6iohI7tCA06LAaoXV42HLLNt2yyFw1+uM/HE/3+86R53yxShX/Uu2XdxGLZ9aLOi2ACezk31rFhGRAkUDTuVPaSnw0/Ow/3vbduc3oNUL7Am7yve7zgHQquEhvj2xDVezK2+3e1vBQ0REcpXCR2GWHGsbWHpyLTg4wj2zof5DWK0GE5ceBOCuBmn8eGoeACObjqSSdyV7ViwiIkWAwkdhFRcBCx6E8L3g5AEPfQFV7wTg+11n2XsummKuBhddPiU1OZX2/u3pXb23nYsWEZGiQOGjMIo6AV/dD1dPg3spePQ7KN8YgOjEVKb+egSABvU3sjfmJCVdS/Jaq9cwmUx2LFpERIoKPe1S2FzYA592sQWP4oHw1Kr04AHw3uqjRMWnUMHvDHtjfgbg9dav4+PqY6eCRUSkqCk04WPw4MGEhoayY8cOe5diPyd+h/k9IP4ylK0HT62GklXSXz5yMZYvt57BZI7DKPUtAH1r9qVthbb2qlhERIogdbsUFvu+gyXP2ebwqHQHPPQVuP75qJNhGLy69AAWq5WAWsu5mnqVqsWr8r/G/7Nj0SIiUhQVmjsfRdrmWbDomeuThz0Aj35/U/AA+Hl/OFtPXsHNZwdXCcHJwYkpbafg6uhqp6JFRKSo0p2Pguzvk4e1eB46vwkON2fKhJQ03vz5ECbny7iW/Zk0A4Y1GkYNnxp2KFpERIo6hY+CKi0FfhoM+7+zbd81CVoNhQyeWJm99gTh0XF4V/mONCOZ5uWa0692vzwuWERExEbhoyBKjoPv+tkGmDo4wj3BUP/hDJueiYpn3vqTOJdeg9X5LF7OXrzZ+k0cTOpxExER+yg04SM4OJjg4GAsFou9S8ldcZfh6962R2qd3KHPl1Dtzn9s/vryUCzOx3EvuQ6Aia0mUsajTF5VKyIicotC88/fIvGo7ZWT8MldtuDhXhL6L//X4LH2SARrjpzBtfx3YDK4t+q93BV4Vx4WLCIicqtCc+ej0LsQYpsuPf4yFA+AxxZDqar/2Dw5zcKkZaG4lv0JB6dr+Hv6M6bZmLyrV0RE5B8ofBQEJ9bCwscgJQ7K1oVHfwTPf+86+XTjac6mbMStdAhmk5nJbSfj4eSRRwWLiIj8M4WP/G7/D7B4EFhToVI7eGjBLXN4/N3F6CQ+WL8N1wpLAHi2/rPUL10/D4oVERH5b4VmzEehtCUYfnzKFjzq3AeP/vCfwQPgrRUHMUp/g8mcTIPSDXim7jN5UKyIiEjm6M5HfmS1wppXYfP7tu3mg6DL5FsmD8vI9lNX+PXcN7j4nsbN7M5bbd/C0UH/m0VEJP/Qt1J+Y0mFn4bAPtvCb9w5EVoPz3DysFveajUY+/PPOJdeDcC4Fi/j7+mfe7WKiIhkgcJHfpIcB989Did+A5MZ7pkFDfpm+u3ztxzmovOnOJistC9/J72q9MrFYkVERLJG4SO/iI+EBb3hwu7rk4d9AdUyPyfHlfgUZuyZjoNnJJ6OpXij7URMmbhbIiIiktcKzYDT4OBgateuTdOmTe1dyu27cur65GG7wc0H+i+7reAB8OLPX4HnVjBMvNt+Mt4u3rlUrIiISPaYDMMw7F1EToqJicHb25vo6Gi8vP77yRC7C98LXz0I8RHXJw9bBKWq3dYhNpw8wXNr+2JyTKBLhYd5p9O4XCpWREQkY7fz/atuF3s6uQ6+fRRSYqFMXXjsB/Ase1uHsFgtvPTHy5gcEyhmCuCt9iNzqVgREZGcofBhLwd+hEXP2ubwqNgWHl4ArrffVTJmzRwSzKEYVkdm3DUNZ7NzLhQrIiKScwrNmI8CZesc+OFJW/CofS889mOWgkfIxVB+vfAJAG1LPknzCrVzuFAREZGcpzsfeckwYM1E2DTDtt1sINw9BRzMt32oZEsyQ9aMBFMaTsl1eK/rczlaqoiISG5R+MgrllRY+gLs/ca23WkCtBmRqcnDMvLaxneItoRhTfNgQsuJuDrpf6WIiBQM+sbKCynx8F1/OL7aNnlYrw+g4aNZPtzm85tZdto2A2oN81PcW7dmTlUqIiKS6xQ+clt8JHzdB87vAkc36PM5VO+S5cNdS7rGyHVjAbBca8n0RzM/A6qIiEh+oPCRm66ehi/vhysnwK0E9P0e/LM+CZphGIzf9CoxqVewJJfmseqDqVjKI+fqFRERyQN62iW3XNwPn3S2BQ9vf3hyVbaCB8Di44v549zvGIYZj+j+DOtYJ4eKFRERyTu685EbTq23TR6WHAO+dWyP0nqVy9Yhz8ScYfK2KQCkRHTmrS534eGi/30iIlLw6Nsrpx1cDIsGgiUFAtvYJg9zK56tQ6ZaUxm7YSxJlkTS4ivTwPseetbLXpgRERGxl0LT7ZIvFpbbNhe+f8IWPGr1st3xyGbwAJi7dy77I/djWFxJCe/Da73qasVaEREpsLSwXE4wDPhtEmycbttu+jR0nZqlycP+bvel3Tyx8gmshpXEc4/waFAvXrsnKNvHFRERyUlaWC4vWVJh2TAIWWDb7vgKtH0py5OH/VVsSiwvb3wZq2El9VojvKxNGXFXjWwfV0RExJ4UPrIjJR6+HwDHVtkmD+s5Exr1y7HDT942mfNx5yHVh6RLvXj13hp4uzvl2PFFRETsQeEjq+Kjrk8ettM2eVjv+VDj7hw7/K+nfmXZyWWAiYTzfahbrgx9mvjn2PFFRETsReEjK66ega8egKhj1ycP+w78m+XY4cPjwpm0dRIAyZc7YEmsyGv31MHsoEGmIiJS8Cl83K6LB2zBI+4ieFWAfougdM6Nw7BYLby88WViU2JxsVQkNrITDzauQKOAEjn2GSIiIvak8HE7Tm2Ab/tenzys9vXJw/xy9CPmH5zPzks7cTK5EnXqQTxdXBh9txaOExGRwkPhI7MOLoFFz1yfPKw1PPx1jszh8VehUaHMCpkFgDXyHozUUgzrXo3Sni45+jkiIiL2VGgmGctVId/YnmqxpECtnvDYohwPHolpiYxeP5o0axp+Ts24FtGAar7F6N+qYo5+joiIiL3pzkdmVGhiG1ha517o9k6OTB72d+/ufJfTMacp4VKK4wfuBkxM7FUHJ7PyoYiIFC4KH5lRqhoM2mgb35EL05qvO7uOhUcWAuAZ+xiWNHe6BpWlddVSOf5ZIiIi9qZ/VmeWd/lcCR6RiZFM2DwBgDal7+fgibK4OjkwrnutHP8sERGR/EDhw44Mw2DCpglcSbpC1eLV2LO3FQDP3VGVCiXc7VydiIhI7lD4sKNvj3zLhvMbcHZwprb5OS5Gp1GhhBvP3lHZ3qWJiIjkGoUPOzlx7QTv7nwXgP41B/P9ljQAxveojatTzg9oFRERyS8UPuwgxZLCmA1jSLYk09qvNSEHg0ixWGlbrRSda5exd3kiIiK5Kt+Fj7Nnz9K+fXtq165NvXr1+P777+1dUo77YM8HHL5ymBIuJbi7zDB+P3wZRwcTr/asgykXBrWKiIjkJ/nuUVtHR0dmzJhBgwYNiIiIoFGjRnTr1g0PDw97l5YjtoVv4/ODnwPwSotXeev7iwA82aYSVX2L2bM0ERGRPJHvwke5cuUoV64cAL6+vvj4+HDlypVCET6ik6N5eePLGBg8WP1Bjp8O5HTUEUp7uvBCx6r2Lk9ERCRP3Ha3y/r16+nZsyd+fn6YTCaWLFlyS5vZs2dTqVIlXF1dady4MRs2bMhScTt37sRqteLv75+l9+cnhmEwacskIhIiqOhVkcerD2XW78cBGNu1Jp6uTnauUEREJG/cdviIj4+nfv36zJo1K8PXFy5cyPDhwxk3bhx79uyhbdu2dO3albCwsPQ2jRs3Jigo6JafCxcupLeJiori8ccfZ968eVk4rfxn6YmlrDqzCkeTI1PaTuG9VadJSLHQOLAE9zUsb+/yRERE8ozJMAwjy282mVi8eDH33ntv+r7mzZvTqFEj5syZk76vVq1a3HvvvUyePDlTx01OTuauu+7imWeeoV+/fv/ZNjk5OX07JiYGf39/oqOj8fLyur0TyiVnY87y4LIHSUhLYFijYdT1uI+H5m3FZIJlQ9oQVN7b3iWKiIhkS0xMDN7e3pn6/s3Rp11SUlLYtWsXnTt3vml/586d2bx5c6aOYRgGAwYMoGPHjv8ZPAAmT56Mt7d3+k9+66JJs6YxduNYEtISaFymMf1q9ufVpQcBeKRZgIKHiIgUOTkaPiIjI7FYLJQpc/NcFWXKlOHixYuZOsamTZtYuHAhS5YsoUGDBjRo0ID9+/f/Y/uxY8cSHR2d/nP27NlsnUNO+2jfR+y9vBdPJ0/eavMWC3ee5/DFWLzdnBjZuYa9yxMREclzufK0y9/nqjAMI9PzV7Rp0war1Zrpz3JxccHFxeW26ssrey/vZe6+uQCMazEOF0ryzso/AHipc3VKeDjbsToRERH7yNE7H6VKlcJsNt9ylyMiIuKWuyGFXXxqPGPWj8FiWOhWqRvdK3fnnVVHiElKo1Y5L/o2D7R3iSIiInaRo+HD2dmZxo0bs3r16pv2r169mlatWuXkR90iODiY2rVr07Rp01z9nMyasn0K5+LO4efhx7gW49h37hrf7rB1Cb3Wqw5mB81kKiIiRdNtd7vExcVx/Pjx9O1Tp04REhKCj48PAQEBjBgxgn79+tGkSRNatmzJvHnzCAsLY9CgQTla+N8NHjyYwYMHp4+2tadVp1ex5PgSHEwOvNX2LYo5evLq0s0YBtzTwI9mlXzsWp+IiIg93Xb42LlzJx06dEjfHjFiBAD9+/dn/vz5PPTQQ0RFRTFp0iTCw8MJCgpixYoVBAYWjW6Gi/EXeW3LawA8FfQUjcs05odd59gTdg13ZzNju9ayc4UiIiL2la15PvKT4OBggoODsVgsHD161C7zfFgNKwNXD2Rb+DbqlKzDl92+JDEFOr6zjsi4ZEbfXZPn2lfJ05pERETygt3m+bCnwYMHExoayo4dO+xWw5ehX7ItfBtujm5MaTsFJwcn3l9zjMi4ZCqV8uDJNhXtVpuIiEh+UWjCh70dvnKYmbtnAjCq6SgqelfkeEQs8zefBmBCz9q4OJrtWKGIiEj+oPCRA5LSkhi9fjSp1lQ6+HfggWoPYBgGE5eGkmY1uLOWLx1q+Nq7TBERkXxB4SMHTN81nZPRJynlVorXWr2GyWRi5cGLbDweibOjA+N71LZ3iSIiIvlGoQkf9prnY8O5DXxz+BsA3mj9BiVcS5CYYuH15YcAGNi2MoElPfK0JhERkfys0IQPeww4jUqMYvym8QA8WutRWpdvDcCH605w/loift6uPN9BT7eIiIj8VaEJH3nNMAwmbp5IVFIUVYtXZXij4QCcvZLAh+tOAPBy91q4O+fK8jkiIiIFlsJHFn1/9Hv+OPcHTg5OTGk7BVdHVwDe+DmU5DQrLSuXpHvdcnauUkREJP9R+MiCk9EnmbZjGgDDGw2nhk8NADYcu8zKg5cwO5iY2KtOplfyFRERKUoKTfjIqwGnqZZUxqwfQ5IliRblWvBY7ccASEmzMnHpQQD6tQikRlnPXK1DRESkoCo04SOvBpwGhwRz6MohvF28ebPNmziYbJfw882nOXE5npIezvzvruq5WoOIiEhBVmjCR17YcXEHnx74FICJLSfi626bOCwiJomZvx0DYNTdNfB2c7JbjSIiIvmdwkcmxaTE8PLGlzEwuL/a/dwZeGf6a1N+PUxcchr1K3jTu7G/HasUERHJ/xQ+MsEwDN7Y8gYX4y8S4BnA6Kaj01/bdeYKi3afB2Birzo4OGiQqYiIyL9R+MiE5SeX88vpXzCbzExuOxl3J3cALFaDV68PMu3duAINA0rYs0wREZECodDMgBUcHExwcDAWiyXHj12nZB1q+dSiY0BH6pWul75/4Y6zHDgfg6eLI6PurpnjnysiIlIYmQzDMOxdRE6KiYnB29ub6OhovLy8cuy4qZZUTCYTjg62vHYtIYUO7/zB1YRUxveozVNtKuXYZ4mIiBQ0t/P9W2jufOQ2J/PNT7BMX32UqwmpVPMtxuMtA+1UlYiISMGjMR9ZEHohhq+2ngHgtV51cDLrMoqIiGSWvjVvk2EYTFx6EKsB3eqWpVXVUvYuSUREpEBR+LhNS/deYPvpK7g6OTCue217lyMiIlLgKHzchvjkNN5acQiA59tXpXxxNztXJCIiUvAUmvCRFwvLzVp7nEsxyfj7uDGwXeVc+xwREZHCTI/aZtKpyHg6v7eOVIvBvH6N6VynbI4dW0REpKC7ne/fQnPnI7dNWnaQVItBu+qluat2GXuXIyIiUmApfGTCb4cusfbIZZzMJl7tWRuTSeu3iIiIZJUmGcuE4u5O1CjjSfsapalSupi9yxERESnQFD4yoXGgDz8PbUOatVANjxEREbELhY9McjQ74Gi2dxUiIiIFn8Z8iIiISJ5S+BAREZE8pfAhIiIiearQhI+8mOFUREREsk8znIqIiEi2aYZTERERybcUPkRERCRPKXyIiIhInlL4EBERkTyl8CEiIiJ5SuFDRERE8pTCh4iIiOQphQ8RERHJU4VuVdsbc6bFxMTYuRIREZGi48b3bmbmLi104SM2NhYAf39/O1ciIiJS9MTGxuLt7f2vbQrd9OpWq5ULFy7g6emJyWQiJiYGf39/zp49q+nW/0LXJWO6LhnTdcmYrkvGdF0yVtivi2EYxMbG4ufnh4PDv4/qKHR3PhwcHKhQocIt+728vArl/+zs0nXJmK5LxnRdMqbrkjFdl4wV5uvyX3c8btCAUxEREclTCh8iIiKSpwp9+HBxceHVV1/FxcXF3qXkK7ouGdN1yZiuS8Z0XTKm65IxXZc/FboBpyIiIpK/Ffo7HyIiIpK/KHyIiIhInlL4EBERkTyl8CEiIiJ5qlCEj8mTJ9O0aVM8PT3x9fXl3nvv5ciRIze1MQyDiRMn4ufnh5ubG+3bt+fgwYN2qtg+Jk+ejMlkYvjw4en7iup1OX/+PI899hglS5bE3d2dBg0asGvXrvTXi+J1SUtL45VXXqFSpUq4ublRuXJlJk2ahNVqTW9TFK7L+vXr6dmzJ35+fphMJpYsWXLT65m5BsnJybzwwguUKlUKDw8PevXqxblz5/LwLHLev12X1NRURo8eTd26dfHw8MDPz4/HH3+cCxcu3HSMonZd/u7ZZ5/FZDIxY8aMm/YXxuvyXwpF+Fi3bh2DBw9m69atrF69mrS0NDp37kx8fHx6m6lTpzJ9+nRmzZrFjh07KFu2LHfddVf6WjCF3Y4dO5g3bx716tW7aX9RvC5Xr16ldevWODk58csvvxAaGsq7775L8eLF09sUxevy9ttv8+GHHzJr1iwOHTrE1KlTmTZtGh988EF6m6JwXeLj46lfvz6zZs3K8PXMXIPhw4ezePFivv32WzZu3EhcXBw9evTAYrHk1WnkuH+7LgkJCezevZvx48eze/duFi1axNGjR+nVq9dN7YradfmrJUuWsG3bNvz8/G55rTBel/9kFEIREREGYKxbt84wDMOwWq1G2bJljSlTpqS3SUpKMry9vY0PP/zQXmXmmdjYWKNatWrG6tWrjTvuuMMYNmyYYRhF97qMHj3aaNOmzT++XlSvS/fu3Y0nn3zypn3333+/8dhjjxmGUTSvC2AsXrw4fTsz1+DatWuGk5OT8e2336a3OX/+vOHg4GD8+uuveVZ7bvr7dcnI9u3bDcA4c+aMYRhF+7qcO3fOKF++vHHgwAEjMDDQeO+999JfKwrXJSOF4s7H30VHRwPg4+MDwKlTp7h48SKdO3dOb+Pi4sIdd9zB5s2b7VJjXho8eDDdu3fnzjvvvGl/Ub0uS5cupUmTJvTu3RtfX18aNmzIRx99lP56Ub0ubdq04bfffuPo0aMA7N27l40bN9KtWzeg6F6Xv8rMNdi1axepqak3tfHz8yMoKKjIXCew/R42mUzpdxSL6nWxWq3069ePkSNHUqdOnVteL6rXpdAtLGcYBiNGjKBNmzYEBQUBcPHiRQDKlClzU9syZcpw5syZPK8xL3377bfs3r2bHTt23PJaUb0uJ0+eZM6cOYwYMYKXX36Z7du3M3ToUFxcXHj88ceL7HUZPXo00dHR1KxZE7PZjMVi4c033+SRRx4Biu6fl7/KzDW4ePEizs7OlChR4pY2N95f2CUlJTFmzBj69u2bvoBaUb0ub7/9No6OjgwdOjTD14vqdSl04WPIkCHs27ePjRs33vKayWS6adswjFv2FSZnz55l2LBhrFq1CldX139sV9Sui9VqpUmTJrz11lsANGzYkIMHDzJnzhwef/zx9HZF7bosXLiQr776iq+//po6deoQEhLC8OHD8fPzo3///untitp1yUhWrkFRuU6pqak8/PDDWK1WZs+e/Z/tC/N12bVrFzNnzmT37t23fY6F+bpAIRlwesMLL7zA0qVLWbt2LRUqVEjfX7ZsWYBbUmRERMQt/4IpTHbt2kVERASNGzfG0dERR0dH1q1bx/vvv4+jo2P6uRe161KuXDlq1659075atWoRFhYGFN0/LyNHjmTMmDE8/PDD1K1bl379+vG///2PyZMnA0X3uvxVZq5B2bJlSUlJ4erVq//YprBKTU2lT58+nDp1itWrV9+0bHxRvC4bNmwgIiKCgICA9N/BZ86c4cUXX6RixYpA0bwuUEjCh2EYDBkyhEWLFvH7779TqVKlm16vVKkSZcuWZfXq1en7UlJSWLduHa1atcrrcvNMp06d2L9/PyEhIek/TZo04dFHHyUkJITKlSsXyevSunXrWx7FPnr0KIGBgUDR/fOSkJCAg8PNvxLMZnP6o7ZF9br8VWauQePGjXFycrqpTXh4OAcOHCjU1+lG8Dh27Bhr1qyhZMmSN71eFK9Lv3792Ldv302/g/38/Bg5ciQrV64EiuZ1AQrH0y7PPfec4e3tbfzxxx9GeHh4+k9CQkJ6mylTphje3t7GokWLjP379xuPPPKIUa5cOSMmJsaOlee9vz7tYhhF87ps377dcHR0NN58803j2LFjxoIFCwx3d3fjq6++Sm9TFK9L//79jfLlyxvLly83Tp06ZSxatMgoVaqUMWrUqPQ2ReG6xMbGGnv27DH27NljAMb06dONPXv2pD+1kZlrMGjQIKNChQrGmjVrjN27dxsdO3Y06tevb6SlpdnrtLLt365Lamqq0atXL6NChQpGSEjITb+Hk5OT049R1K5LRv7+tIthFM7r8l8KRfgAMvz57LPP0ttYrVbj1VdfNcqWLWu4uLgY7dq1M/bv32+/ou3k7+GjqF6XZcuWGUFBQYaLi4tRs2ZNY968eTe9XhSvS0xMjDFs2DAjICDAcHV1NSpXrmyMGzfupi+PonBd1q5dm+Hvk/79+xuGkblrkJiYaAwZMsTw8fEx3NzcjB49ehhhYWF2OJuc82/X5dSpU//4e3jt2rXpxyhq1yUjGYWPwnhd/ovJMAwjL+6wiIiIiEAhGfMhIiIiBYfCh4iIiOQphQ8RERHJUwofIiIikqcUPkRERCRPKXyIiIhInlL4EBERkTyl8CEiIiJ5SuFDRERE8pTCh4iIiOQphQ8RERHJUwofIiIikqf+D9pG2Wyv2RoTAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 620x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAHBCAYAAABOsDAVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABNe0lEQVR4nO3deVzUdf4H8NcwM8wAAgrIfXoi3kLihXQYqKnZoZYb2u7+2ty2DC3FoyztANsOa/Oo1t3OTVPUrNREU0TBCxEP8EoORRDxmEG5hpnP7w9icgJxBo/hy7yej8c8dvl8P9/v5/39RMyr73zn85UJIQSIiIiIJMDO2gUQERERmYvBhYiIiCSDwYWIiIgkg8GFiIiIJIPBhYiIiCSDwYWIiIgkg8GFiIiIJIPBhYiIiCSDwYWIiIgkg8GFiCwik8nMem3fvv2mxwoODsbTTz99x2smotZDYe0CiEhaMjIyTH5+4403sG3bNvzyyy8m7WFhYXezLCKyEQwuRGSRAQMGmPzcvn172NnZNWi/24QQqKqqgoODg1XrIKI7ix8VEdFtd+nSJTz33HPw8/ODvb09OnTogLlz56K6uvqm+2q1Wrz88ssICQmBvb09/Pz8EB8fj2vXrpn0k8lkeP7557Fs2TJ069YNKpUKX3zxBQBg/vz5iIyMhJubG1xcXNCvXz8sX74cf3ymbHBwMEaNGoVNmzahX79+cHBwQGhoKP7zn/+Y9KuoqDDWpFar4ebmhoiICHz77bcm/fbv348xY8bAzc0NarUaffv2xXfffdecKSSiG+AVFyK6raqqqnDffffh119/xfz589GrVy+kpaUhMTERBw8exE8//XTDfSsqKhAdHY2zZ89izpw56NWrF44ePYp58+bh8OHD2LJlC2QymbH/unXrkJaWhnnz5sHb2xuenp4AgPz8fDz77LMIDAwEAOzevRsvvPACioqKMG/ePJMxs7Oz8dJLL2HWrFnw8vLCv//9b/z1r39Fp06dMHToUADA9OnT8dVXX+HNN99E3759ce3aNRw5cgQXL140Hmfbtm0YPnw4IiMjsWzZMri6umLFihWYMGECKioqeC8P0e0iiIhuweTJk4WTk5Px52XLlgkA4rvvvjPpt3DhQgFAbN682dgWFBQkJk+ebPw5MTFR2NnZiX379pnsu3r1agFAbNiwwdgGQLi6uopLly41WZ9erxc6nU4sWLBAuLu7C4PBYDK+Wq0WBQUFxrbKykrh5uYmnn32WWNbjx49xNixY5scJzQ0VPTt21fodDqT9lGjRgkfHx+h1+ub3J+IzMOPiojotvrll1/g5OSExx9/3KS9/orD1q1bb7jvjz/+iB49eqBPnz6ora01vmJjYxv9ptL999+Pdu3aNVrDsGHD4OrqCrlcDqVSiXnz5uHixYsoLS016dunTx/jlRkAUKvV6NKlCwoKCoxt/fv3x8aNGzFr1ixs374dlZWVJsc4deoUjh07hj/96U8AYFL7yJEjUVxcjOPHj9940ojIbAwuRHRbXbx4Ed7e3iYf6QCAp6cnFAqFyccrf3T+/HkcOnQISqXS5OXs7AwhBMrKykz6+/j4NDjG3r17ERMTAwD47LPPsGvXLuzbtw9z584FgAahw93dvcExVCqVSb+PPvoICQkJWLduHe677z64ublh7NixOHnypLFuAHj55Zcb1P7cc88BQIPaiah5eI8LEd1W7u7u2LNnD4QQJuGltLQUtbW18PDwuOG+Hh4ecHBwaHBz7PXbr/fHcAQAK1asgFKpxI8//gi1Wm1sX7dunYVn8jsnJyfMnz8f8+fPx/nz541XX0aPHo1jx44Z65o9ezYeffTRRo/RtWvXZo9PRL9jcCGi2+qBBx7Ad999h3Xr1uGRRx4xtn/55ZfG7TcyatQovP3223B3d0dISEizxpfJZFAoFJDL5ca2yspKfPXVV8063h95eXnh6aefRnZ2NhYtWoSKigp07doVnTt3RnZ2Nt5+++3bMg4RNY7BhYhuq0mTJmHx4sWYPHky8vPz0bNnT+zcuRNvv/02Ro4ciWHDht1w3/j4eCQnJ2Po0KGYNm0aevXqBYPBgMLCQmzevBkvvfQSIiMjmxz/oYcewvvvv4+JEyfib3/7Gy5evIh3330XKpWq2ecUGRmJUaNGoVevXmjXrh1yc3Px1VdfYeDAgXB0dAQAfPLJJxgxYgRiY2Px9NNPw8/PD5cuXUJubi4OHDiAVatWNXt8IvodgwsR3VZqtRrbtm3D3Llz8c9//hMXLlyAn58fXn75Zbz22mtN7uvk5IS0tDQkJSXh008/RV5eHhwcHBAYGIhhw4YhODj4puPff//9+M9//oOFCxdi9OjR8PPzwzPPPANPT0/89a9/bdY53X///Vi/fj0++OADVFRUwM/PD5MmTTLeNwMA9913H/bu3Yu33noL8fHxuHz5Mtzd3REWFobx48c3a1wiakgmxB9WZCIiIiJqofitIiIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgyu43Idg8GAc+fOwdnZudGlxImIiOjOEEKgvLwcvr6+sLO78XUVBpfrnDt3DgEBAdYug4iIyGadOXMG/v7+N9zO4HIdZ2dnAHWT5uLiYuVqiIiIbIdWq0VAQIDxvfhGGFyuU//xkIuLC4MLERGRFdzsVg3enEtERESSweBCREREksHgQkRERJLBe1yaQa/XQ6fTWbsMq1MqlZDL5dYug4iIbAiDiwWEECgpKcGVK1esXUqL0bZtW3h7e3PdGyIiuisYXCxQH1o8PT3h6Oho02/WQghUVFSgtLQUAODj42PlioiIyBYwuJhJr9cbQ4u7u7u1y2kRHBwcAAClpaXw9PTkx0ZERHTH8eZcM9Xf0+Lo6GjlSlqW+vngPT9ERHQ3MLhYyJY/HmoM54OIiO4mBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXG2EwGLBw4UJ06tQJKpUKgYGBeOuttwAA6enp6NOnD9RqNSIiIrBu3TrIZDIcPHjQukUTEVGLc+5KJb7bf8Zq43Mdl1sghEClTm+VsR2Ucou+0TN79mx89tln+OCDDzBkyBAUFxfj2LFjKC8vx+jRozFy5Ej873//Q0FBAeLj4+9c4UREJDlVOj0255zHqv1nsPNUGYQA7gl2Q4iH012vhcHlFlTq9Aib97NVxs5ZEAtHe/P+8ZWXl+PDDz/Exx9/jMmTJwMAOnbsiCFDhmDZsmWQyWT47LPPoFarERYWhqKiIjzzzDN3snwiImrhhBA4dFaDVZlnsP7gOWirao3bIkPccK26tom97xwGFxuQm5uL6upqPPDAAw22HT9+HL169YJarTa29e/f/26WR0RELciF8mqsyyrCqswzOHH+qrHdr60DHuvnh8fC/RHkfvevtNRjcLkFDko5chbEWm1ss/v+tjR/Y4QQDT5yEkI0uy4iIpIend6AX46VYtX+s9h+vBS1hrr3AZXCDsN7eGNceAAGdXSHnZ31Fx1lcLkFMpnM7I9rrKlz585wcHDA1q1b8X//938m20JDQ/HNN9+guroaKpUKALB//35rlElERHfZsRItVu0/i3VZRbh4rcbY3iegLcZF+GNUL1+4OiitWGFDLf9dl26ZWq1GQkICZs6cCXt7ewwePBgXLlzA0aNHMXHiRMydOxd/+9vfMGvWLBQWFuLdd98FwOX8iYhaoysVNViffQ6r9p/F4SKNsb29swqP9vXD4+H+6OzlbMUKm8bgYiNeffVVKBQKzJs3D+fOnYOPjw+mTJkCFxcX/PDDD/j73/+OPn36oGfPnpg3bx4mTpxoct8LERFJl94gkHbyAlZlnkXK0fOo0RsAAEq5DA+EemFchD+iu7SHQt7yV0lhcLERdnZ2mDt3LubOndtg26BBg5CdnW38+ZtvvoFSqURgYODdLJGIiG6z0xeuYnXmWaw5UIQSbZWxvZuPC8aF+2NsXz+4OdlbsULLMbgQvvzyS3To0AF+fn7Izs5GQkICxo8f3+RNvURE1DJdra7FT4fqPgraX3DZ2N7WUYmxfeo+Curh52rFCm8NgwuhpKQE8+bNQ0lJCXx8fDBu3DjjqrpERNTyGQwCu/MuYvX+s9h4pMS4OKqdDIju0h7jIgLwQDdPqBTmfyO1pWJwIcycORMzZ860dhlERGShs5crkJxZhNUHzuDMpUpje4f2ThgXHoBH+/nBy6V13a/I4EJERCQhlTV6/Hy0BKsyzyD914uoX3qrjUqB0b198Hh4APoFtm213wxlcCEiImrhhBA4UHgFqzPP4MfsYpRft9z+oI7uGBfhj+HdfeBgL/2Pgm6GwYWIiKiFOq+twpoDRVideQa/XrhmbA9wc8Dj/eo+Cgpwc7RihXcfgwsREVELUl2rx9bcUqzafwapJy7gt9X3oVbaYWQPH4yLCEBkiFuLWH7fGhhciIiIWoAjRRqszjyL7w8W4XKFztgeHtQO48L98VAvHzirW9by+9bA4EJERGQll67V/PYk5rPILdYa271cVHisnz8eD/dHh/ZtrFhhy8PgYgPuvfde9OnTB4sWLbJ2KURENq9Wb0DqiQtYtf8sth47D52+7rMge7kdHuzuhXHh/ojq3B5yG/0o6GYYXIiIiO6CU6XlWLX/LNZkFeFCebWxvaefK8ZF+GNMb1+0dZTW8vvWwOBCRER0h2irdPjhtycxHzxzxdju7mSPsb89ibmbj4v1CpQgBhcbI5PJsHbtWowdO9bY1rZtWyxatAhPP/201eoiImotDAaB9F8vYlXmGWw6UoLq2ronMcvtZLivqyfGRfjjvq6esFe0/Ccxt0QMLrdCCEBXYZ2xlY5AK10VkYhIigovVmB15hkkHyhC0ZXfl9/v4tUG48IDMLavH9o7q6xYYevA4HIrdBXA277WGXvOOcDeyTpjExERAKCiphYbDpdg1f4z2JN3ydjuolZgTB9fjAsPQC9/11a7/L41MLgQERFZQAiB/QWXsWr/Gfx0qBjXauqexCyTAUM6eWBcRABiwrygVrb+5fetgcHlVigd6658WGvsZpDJZBD1T+T6jU6nu0FvIiKqV6yp/G35/bPIK/t9+f1gd0c8Hu6PR/v5w7etgxUrtA0MLrdCJpPcxzXt27dHcXGx8eeTJ0+iosJK9+kQEbVwVTo9UnLOY1XmWew8+fvy+472cjzUs275/XuC2/GjoLuIwcXG3H///fj4448xYMAAGAwGJCQkQKnkEtJERPWEEDhcpMGq/XXL72urfn8Sc/8QN4wL98fInj5wUvEt1Bo46zbmvffew5///GcMHToUvr6++PDDD5GZmWntsoiIrK7sanXd8vv7z+L4+XJju6+rGo+F1y2/H+QuravsrRGDiw3Yvn278f/7+vri559/Ntl+5cqVu1sQEVELodMbsO1YKVZlnsW2Y6Wo/e2zIJXCDrHdvTEuwh+DOnpw+f0WhMGFiIhszvGScqzafwbrDhah7GqNsb13QFuMC/fH6N6+cHXgx+gtEYMLERHZBE2FDuuz657EfOisxtju0UaFR/v5YVy4Pzp7OVuxQjIHgwsREbVaeoPAzlNlWLX/DDbnnEfNb8vvK+xkeKCbJ8aFByC6a3so5Vx+XyoYXIiIqNX59cJVrDlwFmsOFKFYU2VsD/V2xriIAIzt4wv3Nlx+X4oYXIiISPLKq3TI+PUidpy8gLSTZSi4+Pv6VG0dlXi4ty/GRQSgu68L11yRuGZdG1uyZAlCQkKgVqsRHh6OtLS0JvunpqYiPDwcarUaHTp0wLJlyxr0SU5ORlhYGFQqFcLCwrB27VqT7Tt27MDo0aPh6+sLmUyGdevWNTiGEAKvv/46fH194eDggHvvvRdHjx5tzikSEVELpjcIHDxzBf/aehLjlqWjz4IU/O2rTHy9uxAFFyugsJPh3q7tsXhiP+yZ8wDmP9wDPfz4zKDWwOIrLitXrkR8fDyWLFmCwYMH45NPPsGIESOQk5ODwMDABv3z8vIwcuRIPPPMM/j666+xa9cuPPfcc2jfvj0ee+wxAEBGRgYmTJiAN954A4888gjWrl2L8ePHY+fOnYiMjAQAXLt2Db1798af//xn435/9M477+D999/H559/ji5duuDNN9/Egw8+iOPHj8PZmTdcERFJ2bkrlUg7eQE7TpZh16kyXKkwfVxJiIcTojp7YGjn9hjQ0R1tuEBcqyQTf3xwzU1ERkaiX79+WLp0qbGtW7duGDt2LBITExv0T0hIwPr165Gbm2tsmzJlCrKzs5GRkQEAmDBhArRaLTZu3GjsM3z4cLRr1w7ffvttw6JlMqxduxZjx441tgkh4Ovri/j4eCQkJAAAqqur4eXlhYULF+LZZ5+96blptVq4urpCo9HAxcXFZFtVVRXy8vKMV5qoDueFiO6Uippa7Mm7hLQTZdhx8gJOlV412e6sVmBwRw9EdakLKwFuzXuGG7UMTb0HX8+iOFpTU4PMzEzMmjXLpD0mJgbp6emN7pORkYGYmBiTttjYWCxfvhw6nQ5KpRIZGRmYNm1agz6LFi0yu7a8vDyUlJSYjKVSqRAdHY309HSzggsREVmPEAK5xeW/3adyAfvyLqNGbzBut5PVrbMytHN7DO3igd7+baHgt4FsjkXBpaysDHq9Hl5eXibtXl5eKCkpaXSfkpKSRvvX1tairKwMPj4+N+xzo2PeaJz6/f54nIKCgkb3qa6uRnV1tfFnrVZr9ni2Ij8/HyEhIcjKykKfPn2sXQ4RtTIXyqux89SF366qlKHsarXJdr+2DhjaxQNRndtjcEcPuDpyUThb16wPAP94c5MQoskbnhrr/8d2S495O2pLTEzE/PnzLR6DiIiap7pWj8z8y0g9WRdWcopN/4PRQSnHwI7uGNrZA1Fd2qODhxNvqCUTFgUXDw8PyOXyBldCSktLG1zpqOft7d1of4VCAXd39yb73OiYNxoHqLvy4uPjY9ZxZs+ejenTpxt/1mq1CAgIMHtMIiJqmhACv164ih2/3aey5/QlVOr0Jn16+LkgqnN7DO3cHv2C2kKlkFupWpICiz4ctLe3R3h4OFJSUkzaU1JSMGjQoEb3GThwYIP+mzdvRkREBJRKZZN9bnTMxoSEhMDb29vkODU1NUhNTb3hcVQqFVxcXExerdXq1avRs2dPODg4wN3dHcOGDcO1a9dgMBiwYMEC+Pv7Q6VSoU+fPti0aZO1yyUiCbtSUYMfD51DwupDGJz0C4a9vwMLfszB9uMXUKnTo71z3RL7Hz7RB/tfGYYfX4hCwvBQDOzoztBCN2XxR0XTp09HXFwcIiIiMHDgQHz66acoLCzElClTANRdxSgqKsKXX34JoO4bRB9//DGmT5+OZ555BhkZGVi+fLnJt4VefPFFDB06FAsXLsTDDz+M77//Hlu2bMHOnTuNfa5evYpTp04Zf87Ly8PBgwfh5uaGwMBAyGQyxMfH4+2330bnzp3RuXNnvP3223B0dMTEiRObPUFNEUKgsrbyjhz7ZhwUDmZfPi0uLsaTTz6Jd955B4888gjKy8uRlpYGIQQ+/PBDvPfee/jkk0/Qt29f/Oc//8GYMWNw9OhRdO7c+Q6fBRG1Bjq9AQfPXEHaiQtIPVmGQ2ev4Prvq9or7BAZ4lb3VeUu7dHVy5kf/1CzWRxcJkyYgIsXL2LBggUoLi5Gjx49sGHDBgQFBQGoe5MsLCw09g8JCcGGDRswbdo0LF68GL6+vvjoo49M1mIZNGgQVqxYgVdeeQWvvvoqOnbsiJUrVxrXcAGA/fv347777jP+XP8Rz+TJk/H5558DAGbOnInKyko899xzuHz5MiIjI7F58+Y7toZLZW0lIv8XefOOd8CeiXvgqDTvq3/FxcWora3Fo48+avzn1LNnTwDAu+++i4SEBDzxxBMAgIULF2Lbtm1YtGgRFi9efGeKJyLJK7xY8dt9KheQ8etFlFfXmmzv4tUGQzu3R1SX9ugf7AYHe15JodvD4nVcWjNL13Gp0FVIIrjo9XrExsZi7969iI2NRUxMDB5//HHI5XK4urpi+/btiI6ONvafNm0asrOz8csvv9z0W0Vcx4XINtQvqZ92su5eleuX1AeAdo5KDOncvu6m2s7t4e3KvwdkmTuyjguZclA4YM/EPVYb21xyuRwpKSlIT0/H5s2b8a9//Qtz58413g90u77RRUSth94gcKRIgx0n6p79c6DwMmoNv/93rsJOhvCgdhjapT2iOnugh68r7Oz4d4PuPAaXWyCTycy+6mFtMpkMgwcPxuDBgzFv3jwEBQVh69at8PX1xc6dOzF06FBj3/T0dPTv39+K1RKRNRRrKpF2ogypJy9wSX1qsfhbZwP27NmDrVu3IiYmBp6entizZw8uXLiAbt26YcaMGXjttdfQsWNH9OnTB//9739x8OBBfPPNN9Yum4jusMoaPfbkXcSOE2VIO3kBJ7mkPkkAg4sNcHFxwY4dO7Bo0SJotVoEBQXhvffew4gRIxAbGwutVouXXnoJpaWlCAsLw/r16/mNIqJWqH5J/boHFd54Sf2ozu0RzSX1qYXizbnX4UMWLcd5IWrZuKQ+SQVvziUiskH1S+rvOFmGHScucEl9anUYXIiIJKxuSf1rv3375wJ2N7GkflRnD4QHtePqtCRpDC5ERBJzpaIGu05dNIaVc5oqk+3tnVWI6uyB6C7tMbiTBzzaqKxUKdHtx+BCRNTCXb+k/o7fltQ3cEl9slEMLhbivcymOB9Ed0bhxQrsOHkBO7ikPpEJBhcz1T/JuqKiAg4O5q9a29pVVNQt+10/P0TUPNcvqZ928gLyb7Ckfv0CcFxSn2wVg4uZ5HI52rZti9LSUgCAo6OjTV+KFUKgoqICpaWlaNu2LeRy/tcekSUMBoHDRZq6NVVONL6kfr+gdoj+bUn97r6ukHNJfSIGF0t4e3sDgDG8ENC2bVvjvBDRzZVX6bA68yy+yijA6bJrJtvql9SP6tweA7mkPlGj+G+FBWQyGXx8fODp6QmdTnfzHVo5pVLJKy1EZjpVehVfZuQjOfMsrtXUfV3ZyV6OIb/dUMsl9YnMw+DSDHK5nG/YRHRTeoPAL8dK8UV6PnaeKjO2d/Jsg8kDg/BIP39eVSGyEP+NISK6za5U1GDlvjP4ancBzl6uBFD3HKAHunnh6UHBGNTR3abvkSO6FQwuRES3Sc45Lb7MyMe6g0Wo0tU9vLCtoxIT7gnAU5FB/CiI6DZgcCEiugU6vQGbj57HF+n52Jt/ydjezccFTw8Kwpjeflxjheg2YnAhImqGsqvV+HZPIb7ZU4gSbd2S+3I7GYb38MbTg4IREdSOHwcR3QEMLkREFjh45gq+SM/HT4eKUaOv+zjIo409JvYPxMTIIC4MR3SHMbgQEd1Eda0ePx0qxhcZBcg+c8XY3iegLSYPCsLInj584jLRXcLgQkR0AyWaKnyzpwDf7i1E2dUaAIC93A6jevtg8sBg9A5oa90CiWwQgwsR0XWEENiXfxlfpOdj09ES6H9bht/bRY2nBgTiif6B8GijsnKVRLaLwYWICEBljR7fHyzCFxkFyC3WGtv7h7jh6UHBiAnzgkJuZ8UKiQhgcCEiG3fmUgW+2l2AlfvOQFNZ9ygPtdIOj/T1w6SBwejm42LlConoegwuRGRzhBDYeaoMX6QXYOux8xC/PZQ5wM0BkwYEY1yEP9o62lu3SCJqFIMLEdmMq9W1WHPgLL5Iz8evF35/MnNUZw9MHhiM+0I9Ibfj2itELRmDCxG1er9euIqvMgqwOvMsrlbXAqh7MvPj4f6IGxiMTp5trFwhEZmLwYWIWiW9QWD78VJ8np6PtJO/P5m5Q3snTB4YjEf7+cFZrbRihUTUHAwuRNSqaCp0+G5/3ZOZCy9VAABkMuCBUE9MHhSMwR09YMePg4gki8GFiFqFYyVafJFegHVZRajU6QEALmoFJtwTgLgBwQh055OZiVoDBhcikqxavQEpOefxeXo+9uT9/mTmUG9nTB4UjLF9+GRmotaGwYWIJOfi1Wqs2HcGX+8uQLHm9yczx3b3wuSBwegf4sYnMxO1UgwuRCQZh85ewRfpBfjh0DnU1NY9mdndyR5P9g/EnwYEwsfVwcoVEtGdxuBCRC1aTa0BG48U4/P0fGQVXjG29/Z3xeRBwRjZ0wdqJT8OIrIVDC5E1CKd11bhmz2F+N+eQpRdrQYAKOUyjOrli0kDg9A3sJ2VKyQia2BwIaIWQwiBzILL+Dw9H5uOlKD2tycze7mo8KfIIDzZPxDtnflkZiJbxuBCRFZXpdNj/cFz+CIjH0fP/f5k5nuC22HyoGDEdveGkk9mJiIwuBCRFZ29XIGvdxdi5b5CXK6oezKzSmGHsX38MGlQELr7ulq5QiJqaRhciOiuEkIg49eL+Dw9H1tyz+O3T4Pg19YBkwYGYXxEANo58cnMRNQ4BhciuiuuVddiTVYRvkzPx8nSq8b2IZ08MGlgEB7o5sUnMxPRTTG4ENEdlVd2DV9m5GP1/rMo/+3JzI72cjzWzx+TBwWhk6ezlSskIilhcCGi285gEEg9cQGfp+cj9cQFY3uIhxMmDQzCY+H+cOGTmYmoGRhciOi20VTqsDrzLL7KyEf+xd+fzHxfV09MGhiEoZ3b88nMRHRLGFyI6JadOF+OL9LzsTarCBU1dU9mdlYrMCEiAE8NCEKwh5OVKySi1oLBhYiapVZvwJbcUnyRno+M0xeN7V29nDFpUBAe6esHR3v+iSGi24t/VYjIIpeu1WDFvkJ8s7sQRVcqAQB2MiAmzBuTBwVjQAc+mZmI7hwGFyIyy5EiDb5Iz8f32b8/mdnNyR5P3BOAPw0Igl9bPpmZiO48BhciuqGaWgM2HS3BF+n5yCy4bGzv6Vf3ZOZRvfhkZiK6uxhciKhR6b+WYVbyYRReqvt2kFIuw8iePpg0MBj9Atvy4yAisgoGFyIyUV6lQ9LGY/hmTyEAwKONCk8NCMTE/oHwdFFbuToisnUMLkRklHriAmYnH8I5TRUA4KkBgZg1ohvaqPingohahmY9J37JkiUICQmBWq1GeHg40tLSmuyfmpqK8PBwqNVqdOjQAcuWLWvQJzk5GWFhYVCpVAgLC8PatWstHvfq1at4/vnn4e/vDwcHB3Tr1g1Lly5tzikS2RRNhQ4zVmVj8n/24pymCoFujvj2mQF4c2xPhhYialEsDi4rV65EfHw85s6di6ysLERFRWHEiBEoLCxstH9eXh5GjhyJqKgoZGVlYc6cOZg6dSqSk5ONfTIyMjBhwgTExcUhOzsbcXFxGD9+PPbs2WPRuNOmTcOmTZvw9ddfIzc3F9OmTcMLL7yA77//3tLTJLIZW3LO48EPUrEq8yxkMuAvg0OwKT4KAzu6W7s0IqIGZEIIYckOkZGR6Nevn8mVjG7dumHs2LFITExs0D8hIQHr169Hbm6usW3KlCnIzs5GRkYGAGDChAnQarXYuHGjsc/w4cPRrl07fPvtt2aP26NHD0yYMAGvvvqqsU94eDhGjhyJN95446bnptVq4erqCo1GAxcXF3OnhEiSLl+rwes/HMX3B88BADp4OOGdx3shItjNypURkS0y9z3YoisuNTU1yMzMRExMjEl7TEwM0tPTG90nIyOjQf/Y2Fjs378fOp2uyT71xzR33CFDhmD9+vUoKiqCEALbtm3DiRMnEBsb22ht1dXV0Gq1Ji8iW7DhcDEe/CAV3x88BzsZ8Gx0B2x4MYqhhYhaPIs+vC4rK4Ner4eXl5dJu5eXF0pKShrdp6SkpNH+tbW1KCsrg4+Pzw371B/T3HE/+ugjPPPMM/D394dCoYCdnR3+/e9/Y8iQIY3WlpiYiPnz55t38kStwIXyary2/gg2HK7796aLVxv88/He6B3Q1rqFERGZqVl33f1x/QYhRJNrOjTW/4/t5hzzZn0++ugj7N69G+vXr0dQUBB27NiB5557Dj4+Phg2bFiDumbPno3p06cbf9ZqtQgICLjheRBJlRAC67PP4fX1R3G5QgeFnQzP3dsR/7i/E1QKLiBHRNJhUXDx8PCAXC5vcHWltLS0wdWQet7e3o32VygUcHd3b7JP/THNGbeyshJz5szB2rVr8dBDDwEAevXqhYMHD+Ldd99tNLioVCqoVCpzT59Iks5rqzB37WFsyS0FAIT5uOCf43qhu6+rlSsjIrKcRfe42NvbIzw8HCkpKSbtKSkpGDRoUKP7DBw4sEH/zZs3IyIiAkqlssk+9cc0Z1ydTgedTgc7O9NTksvlMBgMlpwmUasghMB3+89g2Pup2JJbCqVchpdjuuD75wcztBCRdAkLrVixQiiVSrF8+XKRk5Mj4uPjhZOTk8jPzxdCCDFr1iwRFxdn7H/69Gnh6Ogopk2bJnJycsTy5cuFUqkUq1evNvbZtWuXkMvlIikpSeTm5oqkpCShUCjE7t27zR5XCCGio6NF9+7dxbZt28Tp06fFf//7X6FWq8WSJUvMOjeNRiMACI1GY+m0ELUoZy9XiLjle0RQwo8iKOFHMeZfaeJ4idbaZRER3ZC578EWBxchhFi8eLEICgoS9vb2ol+/fiI1NdW4bfLkySI6Otqk//bt20Xfvn2Fvb29CA4OFkuXLm1wzFWrVomuXbsKpVIpQkNDRXJyskXjCiFEcXGxePrpp4Wvr69Qq9Wia9eu4r333hMGg8Gs82JwIanT6w3i6935ovu8TSIo4UfRee4GsWz7KaGr1Vu7NCKiJpn7HmzxOi6tGddxISkrvFiBhORDyDh9EQAQEdQOCx/vhY7t21i5MiKimzP3PZhreRNJnMEg8EVGPt7ZdByVOj0clHLMHN4VkwYGQ27HJzgTUevC4EIkYacvXMXM1Yewv+AyAGBABzcsfKwXgtydrFwZEdGdweBCJEF6g8C/007j/ZQTqK41wMlejtkju2Fi/0DY8SoLEbViDC5EEnPifDlmrD6E7DNXAABRnT2Q+GhP+LdztG5hRER3AYMLkUTo9AZ8kvorPtp6CjV6A5zVCrz6UBjGRfg3uXI1EVFrwuBCJAFHz2kwc/UhHD1X9yDQB0I98dYjPeHtqrZyZUREdxeDC1ELVlNrwMe/nMSS7b+i1iDQ1lGJ10d3x8N9fHmVhYhsEoMLUQuVfeYKZq4+hOPnywEAw7t7Y8HY7vB05lUWIrJdDC5ELUyVTo9FW07i0x2/wiAAdyd7LHi4Bx7q5WPt0oiIrI7BhagFySy4hBmrD+H0hWsAgDG9ffH6mO5wc7K3cmVERC0DgwtRC1BZo8c/fz6O/6bnQQjA01mFN8f2QEx3b2uXRkTUojC4EFlZxq8XMWvNIRRcrAAAPB7uj1cfCoOro9LKlRERtTwMLkRWcrW6Fgs3HsNXuwsAAD6uaiQ+2hP3dvW0cmVERC0XgwuRFaSdvIBZyYdRdKUSADAxMhCzR4TCWc2rLERETWFwIbqLtFU6vP1TLlbsOwMA8G/ngIWP9cLgTh5WroyISBoYXIjukl+OncecNUdQoq0CADw9KBgzYrvCScV/DYmIzMW/mER32JWKGiz4IQdrsooAACEeTlj4WC/0D3GzcmVERNLD4EJ0B206UoJX1h1B2dVq2MmA/4vqgGnDusDBXm7t0oiIJInBhegOuHi1GvPWH8VPh4oBAJ082+Cfj/dC38B2Vq6MiEjaGFyIbiMhBH44VIzX1x/FpWs1kNvJMCW6A6Y+0BkqBa+yEBHdKgYXotukVFuFV9Ydweac8wCAUG9nvDuuN3r4uVq5MiKi1oPBhegWCSGw5kARFvyYA02lDkq5DM/f1xl/v7cj7BV21i6PiKhVYXAhugXFmkrMWXMY245fAAD09HPFP8f1Qqi3i5UrIyJqnRhciJpBCIEV+87g7Z9yUV5dC3uFHeKHdcbfojpAIedVFiKiO4XBhchCZy5VYNaaQ9h16iIAoG9gW/zz8V7o5Ols5cqIiFo/BhciMxkMAl/vKUDSxmOoqNFDrbTDyzFd8efBIZDbyaxdHhGRTWBwITJDftk1zEw+hL15lwAA/UPc8M5jvRDs4WTlyoiIbAuDC1ET9AaB/+7Kw7ubj6NKZ4CjvRyzR4TiT5FBsONVFiKiu47BhegGTpWWY8bqQ8gqvAIAGNLJA4mP9kSAm6N1CyMismEMLkR/UKs34JMdp/HhlpOo0RvgrFJg7kPdMOGeAMhkvMpCRGRNDC5E18kt1mLm6kM4XKQBANzXtT3efrQnfFwdrFwZEREBDC5EAICaWgOWbD+FxdtOQacXcHVQ4rXRYXikrx+vshARtSAMLmTzDp/VYMbqbBwrKQcAxIR54c2xPeDporZyZURE9EcMLmSzqnR6fLT1JD7ZcRp6g4Cbkz3mj+mOUb18eJWFiKiFYnAhm3Sg8DJmrj6EU6VXAQCjevlg/pjucG+jsnJlRETUFAYXsimVNXq8n3Icy3fmwSAAjzYqvDm2B4b38LZ2aUREZAYGF7IZe/MuYebqbORfrAAAPNrXD/NGh6Gto72VKyMiInMxuFCrd626Fu9sOoYvMgoAAN4uarz9aA/cH+pl5cqIiMhSDC7Uqu06VYaE5EM4e7kSAPDEPQGY81A3uKiVVq6MiIiag8GFWiVtlQ6JG47h272FAAC/tg5Ieqwnojq3t3JlRER0KxhcqNXZdrwUc9YcRrGmCgAwaWAQZg4PRRsVf92JiKSOf8mp1dBU6LDgxxwkHzgLAAhyd8TCx3phQAd3K1dGRES3C4MLtQqbj5Zg7rojuFBeDZkM+MvgELwc0xUO9nJrl0ZERLcRgwtJWqm2Cm/8lIsfss8BADq2d8I7j/dGeFA7K1dGRER3AoMLSVJ1rR7/2ZmPj385iWs1etjJgGejO+LFBzpDreRVFiKi1orBhSRFCIHNOefx9oZcFPy2kFzvgLZ44+Hu6OXf1rrFERHRHcfgQpJxvKQcC348il2nLgIAPJ1VmDUiFGP7+MHOjg9FJCKyBQwu1OJdvlaDD7acwNe7C2AQgL3CDn+L6oC/39sRTvyKMxGRTeFffWqxdHoDvtldgA+2nISmUgcAGNHDG3NGdkOAm6OVqyMiImtgcKEWKe3kBSz4IQcnS68CAEK9nTFvdBgGdfSwcmVERGRNDC7UouSXXcObP+ViS+55AEA7RyVeju2KJ+4JhJz3sRAR2Ty75uy0ZMkShISEQK1WIzw8HGlpaU32T01NRXh4ONRqNTp06IBly5Y16JOcnIywsDCoVCqEhYVh7dq1zRo3NzcXY8aMgaurK5ydnTFgwAAUFhY25zTpLiqv0iFxQy4e/CAVW3LPQ2Enw18Gh2D7y/fhT5FBDC1ERASgGcFl5cqViI+Px9y5c5GVlYWoqCiMGDHihuEgLy8PI0eORFRUFLKysjBnzhxMnToVycnJxj4ZGRmYMGEC4uLikJ2djbi4OIwfPx579uyxaNxff/0VQ4YMQWhoKLZv347s7Gy8+uqrUKvVlp4m3SUGg8B3+87gvne345Mdp6HTC0R3aY9N8UMxb3QYXB35FGciIvqdTAghLNkhMjIS/fr1w9KlS41t3bp1w9ixY5GYmNigf0JCAtavX4/c3Fxj25QpU5CdnY2MjAwAwIQJE6DVarFx40Zjn+HDh6Ndu3b49ttvzR73iSeegFKpxFdffWXJKRlptVq4urpCo9HAxcWlWccg8+3Lv4T5PxzFkSItAKCDhxNeHRWG+0I9rVwZERHdbea+B1t0xaWmpgaZmZmIiYkxaY+JiUF6enqj+2RkZDToHxsbi/3790On0zXZp/6Y5oxrMBjw008/oUuXLoiNjYWnpyciIyOxbt06S06R7oKiK5V44dssjFuWgSNFWjirFHjloW7YFD+UoYWIiJpkUXApKyuDXq+Hl5eXSbuXlxdKSkoa3aekpKTR/rW1tSgrK2uyT/0xzRm3tLQUV69eRVJSEoYPH47NmzfjkUcewaOPPorU1NRGa6uuroZWqzV50Z1TWaPHBykn8MB72/FD9jnIZMCT/QOxbca9+L+oDrBXNOuWKyIisiHN+laRTGZ6o6QQokHbzfr/sd2cYzbVx2AwAAAefvhhTJs2DQDQp08fpKenY9myZYiOjm5QV2JiIubPn3/Duun2EELgh0PFSNqQi3OaKgBA/xA3zBsVhh5+rlaujoiIpMSi4OLh4QG5XN7g6kppaWmDqyH1vL29G+2vUCjg7u7eZJ/6Y5ozroeHBxQKBcLCwkz6dOvWDTt37my0ttmzZ2P69OnGn7VaLQICAhrtS81z+KwG8384iv0FlwEAfm0dMPehbhjRw7vJsEtERNQYi67N29vbIzw8HCkpKSbtKSkpGDRoUKP7DBw4sEH/zZs3IyIiAkqlssk+9cc0Z1x7e3vcc889OH78uEmfEydOICgoqNHaVCoVXFxcTF50e5SWV2Hm6myMWbwT+wsuw0Epx0sPdsHWl6IxsqcPQwsRETWPsNCKFSuEUqkUy5cvFzk5OSI+Pl44OTmJ/Px8IYQQs2bNEnFxccb+p0+fFo6OjmLatGkiJydHLF++XCiVSrF69Wpjn127dgm5XC6SkpJEbm6uSEpKEgqFQuzevdvscYUQYs2aNUKpVIpPP/1UnDx5UvzrX/8ScrlcpKWlmXVuGo1GABAajcbSaaHfVOlqxbLtp0T3eZtEUMKPIijhR/HitwfEuSsV1i6NiIhaMHPfgy0OLkIIsXjxYhEUFCTs7e1Fv379RGpqqnHb5MmTRXR0tEn/7du3i759+wp7e3sRHBwsli5d2uCYq1atEl27dhVKpVKEhoaK5ORki8att3z5ctGpUyehVqtF7969xbp168w+LwaX5jMYDCLlaImIfucXY2AZ8680sT//krVLIyIiCTD3PdjidVxaM67j0jwnzpfjjR9zkHay7lti7Z1VSBgeikf7+sGOK94SEZEZzH0P5rOKqNmuVNRg0ZaT+Gp3AfQGAXu5Hf4aFYJ/3NcJbVT81SIiotuP7y5ksVq9Af/bW4j3U07gSkXdIoKx3b0wd2QYAt0drVwdERG1ZgwuZJFdp8ow/4ejOHH+KgCgq5cz5o0Ow+BOHlaujIiIbAGDC5ml4OI1vPVTLjbnnAcAtHVU4qWYrnjyngAo5FzxloiI7g4GF2rS1epafPzLKfxnZx5q9AbI7WSIGxCE+GGd0dbR3trlERGRjWFwoUYZDALJB87inZ+P40J5NQAgqrMH5o0KQ2cvZytXR0REtorBhRrILLiE+T/k4NBZDQAg2N0Rr44Kw/2hnlzxloiIrIrBhYyKNZVI2ngM3x88BwBoo1Jg6gOdMHlQMFQKuZWrIyIiYnAhAFU6PT7dcRpLt/+KSp0eMhkwPjwAL8d2RXtnlbXLIyIiMmJwsWFCCPx0uBiJG46h6EolAOCe4HZ4bXR39PBztXJ1REREDTG42KgjRRos+CEHe/MvAQB8XdWYPbIbRvXik5uJiKjlYnCxMWVXq/Huz8excv8ZCAGolXaYEt0Rzw7tCAd73sdCREQtG4OLjaipNeCL9Hx8tPUkyqtrAQBjevti1ohQ+LZ1sHJ1RERE5mFwaeWEEPjlWCne/CkXeWXXAAA9/Vwxb3QY7gl2s3J1RERElmFwacVOlZZjwY+52HHiAgDAo40KM4d3xeP9/GFnx/tYiIhIehhcWiFNhQ6Ltp7AlxkF0BsElHIZ/jIkBM/f1wnOaqW1yyMiImo2BpdWpFZvwLf7zuD9zcdxuUIHAHgwzAtzR3ZDsIeTlasjIiK6dQwurUT6qTIs+DEHx0rKAQCdPdtg3ugwRHVub+XKiIiIbh8GF4krvFiBtzfkYtPREgCAq4MS0x/sgj9FBkIht7NydURERLcXg4tEXa2uxZJtp/DvnXmoqTVAbifDU5GBiB/WBe2c7K1dHhER0R3B4CIxBoPA2qwiLNx0DKXl1QCAwZ3cMW9Ud3T1drZydURERHcWg4uEHCi8jPk/5CD7zBUAQJC7I+aO7IYHw7y4TD8REdkEBhcJKNFUYeGmY1ibVQQAcLKX44UHOuPPg4OhUnCZfiIish0MLi1YlU6Pz3acxpLtv6JSp4dMBjzezx8zhneFp7Pa2uURERHddQwuLZAQAhuPlOCtn3JRdKUSABAe1A6vjQ5DL/+21i2OiIjIihhcWpij5zRY8EMO9uRdAgD4uKoxa0QoxvT25X0sRERk8xhcWoiLV6vx7uYTWLmvEAYBqBR2mBLdEc9Gd4CjPf8xERERAQwuVldTa8CXGfn4cOtJlFfVAgBG9fLBrBGh8G/naOXqiIiIWhYGFyvadrwUb/yYg9MXrgEAuvu64LXR3dE/xM3KlREREbVMDC5WcKr0Kt78KQfbj18AAHi0sceM2K54PDwAcjvex0JERHQjDC53kaZSh4+2nsQX6fmoNQgo5TL8eXAInr+/E1zUSmuXR0RE1OIxuNwFeoPAin2FeG/zCVy6VgMAGNbNE3MfCkOIh5OVqyMiIpIOBpc7rFZvwKNL03HorAYA0MmzDV4dFYboLu2tXBkREZH0MLjcYQq5HfoHuyG/7BqmPdgFTw0IglJuZ+2yiIiIJEkmhBDWLqKl0Gq1cHV1hUajgYuLy207bnmVDjq9gJuT/W07JhERUWti7nswr7jcBc688ZaIiOi24GcWREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBnNCi5LlixBSEgI1Go1wsPDkZaW1mT/1NRUhIeHQ61Wo0OHDli2bFmDPsnJyQgLC4NKpUJYWBjWrl17S+M+++yzkMlkWLRokcXnR0RERC2TxcFl5cqViI+Px9y5c5GVlYWoqCiMGDEChYWFjfbPy8vDyJEjERUVhaysLMyZMwdTp05FcnKysU9GRgYmTJiAuLg4ZGdnIy4uDuPHj8eePXuaNe66deuwZ88e+Pr6Wnp6RERE1ILJhBDCkh0iIyPRr18/LF261NjWrVs3jB07FomJiQ36JyQkYP369cjNzTW2TZkyBdnZ2cjIyAAATJgwAVqtFhs3bjT2GT58ONq1a4dvv/3WonGLiooQGRmJn3/+GQ899BDi4+MRHx9v1rlptVq4urpCo9HAxcXFvAkhIiKiW2bue7BFV1xqamqQmZmJmJgYk/aYmBikp6c3uk9GRkaD/rGxsdi/fz90Ol2TfeqPae64BoMBcXFxmDFjBrp3737T86muroZWqzV5ERERUctlUXApKyuDXq+Hl5eXSbuXlxdKSkoa3aekpKTR/rW1tSgrK2uyT/0xzR134cKFUCgUmDp1qlnnk5iYCFdXV+MrICDArP2IiIjIOpp1c65MJjP5WQjRoO1m/f/Ybs4xm+qTmZmJDz/8EJ9//nmTtVxv9uzZ0Gg0xteZM2fM2o+IiIisw6Lg4uHhAblc3uDqSmlpaYOrIfW8vb0b7a9QKODu7t5kn/pjmjNuWloaSktLERgYCIVCAYVCgYKCArz00ksIDg5utDaVSgUXFxeTFxEREbVcFgUXe3t7hIeHIyUlxaQ9JSUFgwYNanSfgQMHNui/efNmREREQKlUNtmn/pjmjBsXF4dDhw7h4MGDxpevry9mzJiBn3/+2ZLTJCIiopZKWGjFihVCqVSK5cuXi5ycHBEfHy+cnJxEfn6+EEKIWbNmibi4OGP/06dPC0dHRzFt2jSRk5Mjli9fLpRKpVi9erWxz65du4RcLhdJSUkiNzdXJCUlCYVCIXbv3m32uI0JCgoSH3zwgdnnptFoBACh0WgsmBEiIiK6Vea+ByssDToTJkzAxYsXsWDBAhQXF6NHjx7YsGEDgoKCAADFxcUma6uEhIRgw4YNmDZtGhYvXgxfX1989NFHeOyxx4x9Bg0ahBUrVuCVV17Bq6++io4dO2LlypWIjIw0e1wiIiJq/Sxex6U14zouRERE1nFH1nEhIiIisiYGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhciIiKSjGYFlyVLliAkJARqtRrh4eFIS0trsn9qairCw8OhVqvRoUMHLFu2rEGf5ORkhIWFQaVSISwsDGvXrrVoXJ1Oh4SEBPTs2RNOTk7w9fXFpEmTcO7cueacIhEREbVAFgeXlStXIj4+HnPnzkVWVhaioqIwYsQIFBYWNto/Ly8PI0eORFRUFLKysjBnzhxMnToVycnJxj4ZGRmYMGEC4uLikJ2djbi4OIwfPx579uwxe9yKigocOHAAr776Kg4cOIA1a9bgxIkTGDNmjKWnSERERC2UTAghLNkhMjIS/fr1w9KlS41t3bp1w9ixY5GYmNigf0JCAtavX4/c3Fxj25QpU5CdnY2MjAwAwIQJE6DVarFx40Zjn+HDh6Ndu3b49ttvmzUuAOzbtw/9+/dHQUEBAgMDb3puWq0Wrq6u0Gg0cHFxuWl/IiIiuj3MfQ+26IpLTU0NMjMzERMTY9IeExOD9PT0RvfJyMho0D82Nhb79++HTqdrsk/9MZszLgBoNBrIZDK0bdu20e3V1dXQarUmLyIiImq5LAouZWVl0Ov18PLyMmn38vJCSUlJo/uUlJQ02r+2thZlZWVN9qk/ZnPGraqqwqxZszBx4sQbJrfExES4uroaXwEBATc4cyIiImoJmnVzrkwmM/lZCNGg7Wb9/9huzjHNHVen0+GJJ56AwWDAkiVLbljX7NmzodFojK8zZ87csC8RERFZn8KSzh4eHpDL5Q2ucpSWlja4GlLP29u70f4KhQLu7u5N9qk/piXj6nQ6jB8/Hnl5efjll1+a/JxMpVJBpVI1ccZERETUklh0xcXe3h7h4eFISUkxaU9JScGgQYMa3WfgwIEN+m/evBkRERFQKpVN9qk/prnj1oeWkydPYsuWLcZgRERERK2EsNCKFSuEUqkUy5cvFzk5OSI+Pl44OTmJ/Px8IYQQs2bNEnFxccb+p0+fFo6OjmLatGkiJydHLF++XCiVSrF69Wpjn127dgm5XC6SkpJEbm6uSEpKEgqFQuzevdvscXU6nRgzZozw9/cXBw8eFMXFxcZXdXW1Weem0WgEAKHRaCydFiIiIroF5r4HWxxchBBi8eLFIigoSNjb24t+/fqJ1NRU47bJkyeL6Ohok/7bt28Xffv2Ffb29iI4OFgsXbq0wTFXrVolunbtKpRKpQgNDRXJyckWjZuXlycANPratm2bWefF4EJERGQd5r4HW7yOS2vGdVyIiIis446s40JERERkTQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZDC5EREQkGQwuREREJBkMLkRERCQZzQouS5YsQUhICNRqNcLDw5GWltZk/9TUVISHh0OtVqNDhw5YtmxZgz7JyckICwuDSqVCWFgY1q5da/G4Qgi8/vrr8PX1hYODA+69914cPXq0OadIRERELZDFwWXlypWIj4/H3LlzkZWVhaioKIwYMQKFhYWN9s/Ly8PIkSMRFRWFrKwszJkzB1OnTkVycrKxT0ZGBiZMmIC4uDhkZ2cjLi4O48ePx549eywa95133sH777+Pjz/+GPv27YO3tzcefPBBlJeXW3qaRERE1BIJC/Xv319MmTLFpC00NFTMmjWr0f4zZ84UoaGhJm3PPvusGDBggPHn8ePHi+HDh5v0iY2NFU888YTZ4xoMBuHt7S2SkpKM26uqqoSrq6tYtmyZWeem0WgEAKHRaMzqT0RERLeHue/BCktCTk1NDTIzMzFr1iyT9piYGKSnpze6T0ZGBmJiYkzaYmNjsXz5cuh0OiiVSmRkZGDatGkN+ixatMjscfPy8lBSUmIylkqlQnR0NNLT0/Hss89acqq3jTAYUFl1ySpjExER3SkOajfI7O7+rbIWBZeysjLo9Xp4eXmZtHt5eaGkpKTRfUpKShrtX1tbi7KyMvj4+NywT/0xzRm3/n8b61NQUNBobdXV1aiurjb+rNVqG+13KyqrLiFy1X23/bhERETWtGfcNjg6etz1cZsVlWQymcnPQogGbTfr/8d2c455u/rUS0xMhKurq/EVEBBww3MgIiIi67PoiouHhwfkcnmDqyulpaUNrnTU8/b2brS/QqGAu7t7k33qj2nOuN7e3gDqrrz4+PiYVdvs2bMxffp0489arfa2hxcHtRv2jNt2W49JRERkbQ5qN6uMa1Fwsbe3R3h4OFJSUvDII48Y21NSUvDwww83us/AgQPxww8/mLRt3rwZERERUCqVxj4pKSkm97ls3rwZgwYNMnvckJAQeHt7IyUlBX379gVQd29MamoqFi5c2GhtKpUKKpXKkimwmMzOziqX0oiIiFolS+/6XbFihVAqlWL58uUiJydHxMfHCycnJ5Gfny+EEGLWrFkiLi7O2P/06dPC0dFRTJs2TeTk5Ijly5cLpVIpVq9ebeyza9cuIZfLRVJSksjNzRVJSUlCoVCI3bt3mz2uEEIkJSUJV1dXsWbNGnH48GHx5JNPCh8fH6HVas06N36riIiIyDrMfQ+2OLgIIcTixYtFUFCQsLe3F/369ROpqanGbZMnTxbR0dEm/bdv3y769u0r7O3tRXBwsFi6dGmDY65atUp07dpVKJVKERoaKpKTky0aV4i6r0S/9tprwtvbW6hUKjF06FBx+PBhs8+LwYWIiMg6zH0Plgnx252yBK1WC1dXV2g0Gri4uFi7HCIiIpth7nswn1VEREREksHgQkRERJLB4EJERESSweBCREREksHgQkRERJLB4EJERESSweBCREREksHgQkRERJLB4EJERESSweBCREREkmHR06Fbu/qnH2i1WitXQkREZFvq33tv9iQiBpfrlJeXAwACAgKsXAkREZFtKi8vh6ur6w238yGL1zEYDDh37hycnZ0hk8kA1CXAgIAAnDlzhg9evA7npXGcl8ZxXhrinDSO89I4W5gXIQTKy8vh6+sLO7sb38nCKy7XsbOzg7+/f6PbXFxcWu0vy63gvDSO89I4zktDnJPGcV4a19rnpakrLfV4cy4RERFJBoMLERERSQaDy02oVCq89tprUKlU1i6lReG8NI7z0jjOS0Ock8ZxXhrHefkdb84lIiIiyeAVFyIiIpIMBhciIiKSDAYXIiIikgwGFyIiIpIMBhcAiYmJuOeee+Ds7AxPT0+MHTsWx48fN+kjhMDrr78OX19fODg44N5778XRo0etVLF1JCYmQiaTIT4+3thmq/NSVFSEp556Cu7u7nB0dESfPn2QmZlp3G6L81JbW4tXXnkFISEhcHBwQIcOHbBgwQIYDAZjH1uYlx07dmD06NHw9fWFTCbDunXrTLabMwfV1dV44YUX4OHhAScnJ4wZMwZnz569i2dx+zU1LzqdDgkJCejZsyecnJzg6+uLSZMm4dy5cybHaG3zcrPfles9++yzkMlkWLRokUl7a5sTczC4AEhNTcU//vEP7N69GykpKaitrUVMTAyuXbtm7PPOO+/g/fffx8cff4x9+/bB29sbDz74oPH5Rq3dvn378Omnn6JXr14m7bY4L5cvX8bgwYOhVCqxceNG5OTk4L333kPbtm2NfWxxXhYuXIhly5bh448/Rm5uLt555x3885//xL/+9S9jH1uYl2vXrqF37974+OOPG91uzhzEx8dj7dq1WLFiBXbu3ImrV69i1KhR0Ov1d+s0brum5qWiogIHDhzAq6++igMHDmDNmjU4ceIExowZY9Kvtc3LzX5X6q1btw579uyBr69vg22tbU7MIqiB0tJSAUCkpqYKIYQwGAzC29tbJCUlGftUVVUJV1dXsWzZMmuVedeUl5eLzp07i5SUFBEdHS1efPFFIYTtzktCQoIYMmTIDbfb6rw89NBD4i9/+YtJ26OPPiqeeuopIYRtzgsAsXbtWuPP5szBlStXhFKpFCtWrDD2KSoqEnZ2dmLTpk13rfY76Y/z0pi9e/cKAKKgoEAI0frn5UZzcvbsWeHn5yeOHDkigoKCxAcffGDc1trn5EZ4xaURGo0GAODm5gYAyMvLQ0lJCWJiYox9VCoVoqOjkZ6ebpUa76Z//OMfeOihhzBs2DCTdludl/Xr1yMiIgLjxo2Dp6cn+vbti88++8y43VbnZciQIdi6dStOnDgBAMjOzsbOnTsxcuRIALY7L9czZw4yMzOh0+lM+vj6+qJHjx42M09A3d9hmUxmvJJpi/NiMBgQFxeHGTNmoHv37g222+KcAHzIYgNCCEyfPh1DhgxBjx49AAAlJSUAAC8vL5O+Xl5eKCgouOs13k0rVqzAgQMHsG/fvgbbbHVeTp8+jaVLl2L69OmYM2cO9u7di6lTp0KlUmHSpEk2Oy8JCQnQaDQIDQ2FXC6HXq/HW2+9hSeffBKA7f6+XM+cOSgpKYG9vT3atWvXoE/9/q1dVVUVZs2ahYkTJxofKGiL87Jw4UIoFApMnTq10e22OCcAg0sDzz//PA4dOoSdO3c22CaTyUx+FkI0aGtNzpw5gxdffBGbN2+GWq2+YT9bmxeDwYCIiAi8/fbbAIC+ffvi6NGjWLp0KSZNmmTsZ2vzsnLlSnz99df43//+h+7du+PgwYOIj4+Hr68vJk+ebOxna/PSmObMga3Mk06nwxNPPAGDwYAlS5bctH9rnZfMzEx8+OGHOHDggMXn11rnpB4/KrrOCy+8gPXr12Pbtm3w9/c3tnt7ewNAgwRbWlra4L+cWpPMzEyUlpYiPDwcCoUCCoUCqamp+Oijj6BQKIznbmvz4uPjg7CwMJO2bt26obCwEIDt/r7MmDEDs2bNwhNPPIGePXsiLi4O06ZNQ2JiIgDbnZfrmTMH3t7eqKmpweXLl2/Yp7XS6XQYP3488vLykJKSYrzaAtjevKSlpaG0tBSBgYHGv78FBQV46aWXEBwcDMD25qQegwvq0unzzz+PNWvW4JdffkFISIjJ9pCQEHh7eyMlJcXYVlNTg9TUVAwaNOhul3vXPPDAAzh8+DAOHjxofEVEROBPf/oTDh48iA4dOtjkvAwePLjB1+VPnDiBoKAgALb7+1JRUQE7O9M/KXK53Ph1aFudl+uZMwfh4eFQKpUmfYqLi3HkyJFWPU/1oeXkyZPYsmUL3N3dTbbb2rzExcXh0KFDJn9/fX19MWPGDPz8888AbG9OjKx1V3BL8ve//124urqK7du3i+LiYuOroqLC2CcpKUm4urqKNWvWiMOHD4snn3xS+Pj4CK1Wa8XK777rv1UkhG3Oy969e4VCoRBvvfWWOHnypPjmm2+Eo6Oj+Prrr419bHFeJk+eLPz8/MSPP/4o8vLyxJo1a4SHh4eYOXOmsY8tzEt5ebnIysoSWVlZAoB4//33RVZWlvHbMebMwZQpU4S/v7/YsmWLOHDggLj//vtF7969RW1trbVO65Y1NS86nU6MGTNG+Pv7i4MHD5r8Ha6urjYeo7XNy81+V/7oj98qEqL1zYk5GFxE3dfQGnv997//NfYxGAzitddeE97e3kKlUomhQ4eKw4cPW69oK/ljcLHVefnhhx9Ejx49hEqlEqGhoeLTTz812W6L86LVasWLL74oAgMDhVqtFh06dBBz5841eeOxhXnZtm1bo39PJk+eLIQwbw4qKyvF888/L9zc3ISDg4MYNWqUKCwstMLZ3D5NzUteXt4N/w5v27bNeIzWNi83+135o8aCS2ubE3PIhBDiblzZISIiIrpVvMeFiIiIJIPBhYiIiCSDwYWIiIgkg8GFiIiIJIPBhYiIiCSDwYWIiIgkg8GFiIiIJIPBhYiIiCSDwYWIiIgkg8GFiIiIJIPBhYiIiCSDwYWIiIgk4/8BdVxiM7lfZtIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 620x500 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import scipy.linalg as lin\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "from scipy.sparse.linalg import cg\n",
    "from scipy.linalg import lu\n",
    "from scipy.linalg import solve_triangular\n",
    "from math import *\n",
    "\n",
    "\n",
    "\n",
    "class FDMSteadyFlow2D:\n",
    "    def __init__(self, Nx, Ny, h_top):\n",
    "        self.Nx = Nx\n",
    "        self.Ny = Ny\n",
    "        self.h_top = h_top\n",
    "        \n",
    "    def buildCoeffMatrix(self):\n",
    "        '''\n",
    "        Constructs the coefficient matrix A of the final linear\n",
    "        system of equations of the form Ax = b\n",
    "        '''\n",
    "        # Create the diagonal, upper and lower components of B\n",
    "        Ddiag = -4 * np.eye(self.Nx - 1)\n",
    "        Dupper = np.diag([1] * (self.Nx - 2), 1)\n",
    "        Dlower = np.diag([1] * (self.Nx - 2), -1)\n",
    "\n",
    "        # Add components to create B\n",
    "        D = Ddiag + Dupper + Dlower\n",
    "\n",
    "        # Create a block matrix where the diagonals are each B\n",
    "        Ds = [D] * (self.Nx - 1)\n",
    "        A = lin.block_diag(*Ds)\n",
    "\n",
    "        # Create the identity diagonals\n",
    "        I = np.ones((self.Nx - 1) * (self.Nx - 2))\n",
    "        Iupper = np.diag(I, self.Nx - 1)\n",
    "        Ilower = np.diag(I, -self.Nx + 1)\n",
    "\n",
    "        # Add the identity diagonal to A to complete it\n",
    "        A += Iupper + Ilower\n",
    "\n",
    "        return A\n",
    "    \n",
    "    def buildRHSVector(self):\n",
    "        '''\n",
    "        Constructs the right hand side vector b of the final linear\n",
    "        system of equations of the form Ax = b\n",
    "        '''\n",
    "        b = np.zeros((self.Nx - 1)**2)\n",
    "        b[-self.Nx+1:] = -self.h_top\n",
    "\n",
    "        return b\n",
    "\n",
    "    def solveLinearSystem(self, A, b):\n",
    "        '''\n",
    "        Solves the final linear system Ax = b and concatenates the\n",
    "        boundary conditions (BCs) to the solution\n",
    "        '''\n",
    "        # Solve for h vector and reshape array to 2D\n",
    "        h = lin.solve(A, b)\n",
    "        errv= np.dot(A,h)-b\n",
    "        err = sqrt(np.dot(errv,errv))\n",
    "        h = h.reshape((self.Nx - 1, self.Ny - 1))\n",
    "\n",
    "        # Create empty 2D array with all nodes and insert BCs ans solution\n",
    "        h2D = np.zeros((self.Nx + 1, self.Ny + 1))\n",
    "        h2D[0] = self.h_top         # Insert top BC\n",
    "        h2D[1:-1, 1:-1] = h[::-1]    # Insert solution (::-1 => inverted)\n",
    "\n",
    "        return h2D,err\n",
    "    \n",
    "    def solveLinearSystemcg(self, A, b):\n",
    "        '''        Solves the final linear system Ax = b and concatenates the\n",
    "        boundary conditions (BCs) to the solution\n",
    "        '''\n",
    "        # Solve for h vector and reshape array to 2D       \n",
    "        h,info = cg(A,b)\n",
    "        errv= np.dot(A,h)-b\n",
    "        err = sqrt(np.dot(errv,errv))\n",
    "        #print(\"Accuracy cg: \",err)\n",
    "        #print(\"info: \", info)\n",
    "        #print(\" \")\n",
    "        h = h.reshape((self.Nx - 1, self.Ny - 1))\n",
    "\n",
    "        # Create empty 2D array with all nodes and insert BCs ans solution\n",
    "        h2D = np.zeros((self.Nx + 1, self.Ny + 1))\n",
    "        h2D[0] = self.h_top         # Insert top BC\n",
    "        h2D[1:-1, 1:-1] = h[::-1]    # Insert solution (::-1 => inverted)\n",
    "\n",
    "        return h2D,err\n",
    "    \n",
    "    def solveLinearSystemlu(self, A, b):\n",
    "        '''\n",
    "        Solves the final linear system Ax = b and concatenates the\n",
    "        boundary conditions (BCs) to the solution\n",
    "        '''\n",
    "        # Solve for h vector and reshape array to 2D\n",
    "        L,U=lu(A,permute_l=True)\n",
    "        Y=solve_triangular(L,b,lower=True)\n",
    "        h=solve_triangular(U,Y)\n",
    "        errv= np.dot(A,h)-b\n",
    "        err = sqrt(np.dot(errv,errv))\n",
    "        #print(\"Accuracy lu: \",err)\n",
    "        #print(\" \")\n",
    "        h = h.reshape((self.Nx - 1, self.Ny - 1))\n",
    "\n",
    "        # Create empty 2D array with all nodes and insert BCs ans solution\n",
    "        h2D = np.zeros((self.Nx + 1, self.Ny + 1))\n",
    "        h2D[0] = self.h_top         # Insert top BC\n",
    "        h2D[1:-1, 1:-1] = h[::-1]    # Insert solution (::-1 => inverted)\n",
    "\n",
    "        return h2D,err\n",
    "\n",
    "    \n",
    "    def plotSolution(self, h2D):\n",
    "        '''\n",
    "        Plots the color plot of the solution on a 2D meshgrid\n",
    "        '''\n",
    "        # Create 1D arrays with number of nodes\n",
    "        x = np.linspace(0, 1, self.Nx + 1)\n",
    "        y = np.linspace(1, 0, self.Ny + 1)\n",
    "\n",
    "        # Create 2D mesh grid\n",
    "        X, Y = np.meshgrid(x, y)\n",
    "\n",
    "        # Plot solution on mesh grid\n",
    "        matplotlib.rcParams['figure.figsize'] = 6.2, 5\n",
    "        # plt.clf()\n",
    "        plt.contourf(X, Y, h2D, 10)\n",
    "        plt.colorbar()\n",
    "        plt.xlabel(r'$x$ [$\\mathrm{m}$]')\n",
    "        plt.ylabel(r'$y$ [$\\mathrm{m}$]')\n",
    "        #plt.savefig('2D_Flow_Nx_Ny_' + str(Nx) + '.eps')\n",
    "        plt.show()\n",
    "        \n",
    "if __name__ == \"__main__\":\n",
    "  nruns=6\n",
    "  timesol=np.zeros(nruns)\n",
    "  tolsol=np.zeros(nruns)\n",
    "  timelu=np.zeros(nruns)\n",
    "  tollu=np.zeros(nruns)\n",
    "  timecg=np.zeros(nruns)\n",
    "  tolcg=np.zeros(nruns)\n",
    "  dim=np.zeros(nruns)\n",
    "\n",
    "  nx=0\n",
    "  for i in range(0,nruns) :\n",
    "    nx = nx+25\n",
    "    print(\"nx : \", nx)\n",
    "    fdm = FDMSteadyFlow2D(Nx=nx, Ny=nx, h_top=10)\n",
    "    A = fdm.buildCoeffMatrix()\n",
    "    b = fdm.buildRHSVector()\n",
    "    \n",
    "    dim[i]   = nx\n",
    "    \n",
    "    print(\"cg...\")\n",
    "    t0=time.perf_counter()\n",
    "    hcg,errcg = fdm.solveLinearSystemcg(A, b)\n",
    "    t = time.perf_counter()-t0\n",
    "    timecg[i]  = t\n",
    "    tolcg[i] = errcg\n",
    "    print(t,errcg)\n",
    "    \n",
    "    print(\"gesv...\")\n",
    "    t0=time.perf_counter()\n",
    "    h,errsol = fdm.solveLinearSystem(A, b)\n",
    "    t = time.perf_counter()-t0\n",
    "    timesol[i]  = t\n",
    "    tolsol[i] = errsol\n",
    "    print(t,errsol)\n",
    "    \n",
    "    print(\"lu...\")\n",
    "    t0=time.perf_counter()\n",
    "    hlu,errlu = fdm.solveLinearSystemlu(A, b)\n",
    "    t = time.perf_counter()-t0\n",
    "    timelu[i]  = t\n",
    "    tollu[i]=errlu\n",
    "    print(t,errlu)\n",
    "\n",
    "  fdm.plotSolution(hcg)\n",
    "  fdm.plotSolution(h)\n",
    "  fdm.plotSolution(hlu)\n",
    "  fig=plt.figure()\n",
    "  plt.plot(dim,timecg, label='cg')\n",
    "  plt.plot(dim,timelu, label='lu')\n",
    "  plt.plot(dim,timesol,label='sol')\n",
    "  plt.legend()\n",
    "  plt.yscale('log')\n",
    "  plt.title(\"Run time\")\n",
    "  #plt.savefig(\"runTime.pdf\")\n",
    "  plt.show()\n",
    "  plt.plot(dim,tolcg, label='cg')\n",
    "  plt.plot(dim,tollu, label='lu')\n",
    "  plt.plot(dim,tolsol,label='sol')\n",
    "  plt.legend()\n",
    "  plt.title(\"Toleranse\")\n",
    "  #plt.savefig(\"solvers.pdf\")\n",
    "  plt.show()\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "75d74c51",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}