lectures/lecture9/lecture-9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 9: From dense to sparse linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lectures\n",
    "- Algorithms for the symmetric eigenvalue (SEV) problem\n",
    "    - QR algorithm with shifts\n",
    "    - Divide-and-Conquer\n",
    "    - Bisection\n",
    "    - Jacobi method\n",
    "- Reduction problem of computing SVD to SEV problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Large scale dense matrices\n",
    "\n",
    "- If the size of the dense matrix is huge, then it can not be stored in memory\n",
    "- Possible options\n",
    "    - This matrix is **structured**, e.g. block Toeplitz with Toeptitz blocks (next lectures). Then the compressed storage is possible\n",
    "    - For unstructured dense matrices **distributed memory** helps\n",
    "    - MPI for processing distributed storing matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Distributed memory and MPI\n",
    "\n",
    "- Split matrix into blocks and store them on different machines\n",
    "- Every machine has its own address space and can not damage data on other machines\n",
    "- In this case machines communicate with each other to aggregate results of computations\n",
    "- [MPI (Message Passing Interface)](https://en.wikipedia.org/wiki/Message_Passing_Interface) is a standard for parallel computing in distributed memory\n",
    "\n",
    "- Skoltech owns supercomputer, which is on the 7th place in the list of TOP50 supercomputers from Commonwealth of Independent States (see the rating [here](http://top50.supercomputers.ru/?page=rating), sorry for only in Russian)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Example: matrix-by-vector product\n",
    "\n",
    "- Assume you want to compute $Ax$ and matrix $A$ can not be stored in available memory\n",
    "- Then you can split it on blocks and distribute blocks on separate machines\n",
    "- Possible strategies\n",
    "    - 1D blocking splits only rows on blocks\n",
    "    - 2D blocking splits both rows and columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1D blocking scheme\n",
    "<img src=\"1d_block.jpg\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total time of computing matvec with 1D blocking\n",
    "\n",
    "- Each machine has $n / p $ complete rows and $n / p$ elements of vector\n",
    "- Total operations are $n^2 / p$\n",
    "- Total time for sending and writing data are $t_s \\log p + t_w n$, where $t_s$ time unit for sending and $t_w$ time unit for writing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 2D blocking scheme\n",
    "\n",
    "<img src=\"2d_block.png\" width=400>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Total time of computing matvec with 2D blocking\n",
    "\n",
    "- Each machine has $n / \\sqrt{p} $ size block and $n / \\sqrt{p}$ elements of vector\n",
    "- Total operations are $n^2 / p$\n",
    "- Total time for sending and writing data are approximately $t_s \\log p + t_w (n/\\sqrt{p}) \\log p$, where $t_s$ time unit for sending and $t_w$ time unit for writing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Packages supported distributed storage\n",
    "\n",
    "- [ScaLAPACK](http://www.netlib.org/scalapack/)\n",
    "- [Trilinos](https://trilinos.org/)\n",
    "\n",
    "In Python you can use [mpi4py](https://mpi4py.readthedocs.io/en/stable/) for parallel programming of your algorithm.\n",
    "\n",
    "- PyTorch supports distributed training and data storage, see details [here](https://pytorch.org/tutorials/intermediate/dist_tuto.html) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Summary on large unstructered matrix processing\n",
    "\n",
    "- Distributed manner of storage\n",
    "- MPI\n",
    "- Packages that use parallel computrations\n",
    "- Different blocking strategies"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices intro\n",
    "\n",
    "- For dense linear algebra problems, we are limited by the memory to store the full matrix, it is $N^2$ parameters.\n",
    "\n",
    "- The class of **sparse** matrices where most of the elements are zero, allows us **at least** to store such matrices.\n",
    "\n",
    "The question if we can:\n",
    "\n",
    "- solve linear systems\n",
    "- solve eigenvalue problems\n",
    "\n",
    "with sparse matrices"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan for the next part of the lecture\n",
    "\n",
    "Now we will talk about **sparse matrices**, where they arise, how we store them, how we operate with them.\n",
    "\n",
    "\n",
    "- Formats: list of lists and compressed sparse row format, relation to graphs\n",
    "- Matrix-by-vector product\n",
    "- Parallell processing of sparse matrices \n",
    "- Fast direct solvers for Gaussian elimination (start)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Applications of sparse matrices\n",
    "\n",
    "Sparse matrices arise in:\n",
    "\n",
    "- partial differential equations (PDE), mathematical modelling\n",
    "- graphs mining, e.g. social networks analysis\n",
    "- recommender systems\n",
    "- wherever relations between objects are \"sparse\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices are ubiquitous in PDE\n",
    "\n",
    "The simplest partial differential equation (PDE), called \n",
    "\n",
    "**Laplace equation**:  \n",
    "\n",
    "$$\n",
    "   \\Delta T = \\frac{\\partial^2 T}{\\partial x^2} + \\frac{\\partial^2 T}{\\partial y^2} = f(x,y), \\quad x,y\\in \\Omega\\equiv[0,1]^2,\n",
    "$$\n",
    "\n",
    "$$\n",
    "    T_{\\partial\\Omega} = 0.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Discretization\n",
    "\n",
    "$$\\frac{\\partial^2 T}{\\partial x^2} \\approx \\frac{T(x+h) + T(x-h) - 2T(x)}{h^2} + \\mathcal{O}(h^2),$$\n",
    "\n",
    "same for $\\frac{\\partial^2 T}{\\partial y^2},$\n",
    "and we get a linear system.  \n",
    "    First, let us consider **one-dimensional case**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "After the discretization of the one-dimensional Laplace equation with Dirichlet boundary conditions we have\n",
    "\n",
    "$$\\frac{u_{i+1} + u_{i-1} - 2u_i}{h^2} = f_i,\\quad i=1,\\dots,N-1$$\n",
    "\n",
    "$$ u_{0} = u_N = 0$$\n",
    "or in the matrix form\n",
    "\n",
    "$$ A u = f,$$\n",
    "and (for $n = 5$)\n",
    "$$A=-\\frac{1}{h^2}\\begin{bmatrix} 2 & -1 & 0 & 0 & 0\\\\ -1 & 2 & -1 & 0 &0 \\\\ 0 & -1 & 2& -1 & 0 \\\\ 0 & 0 & -1 & 2  &-1\\\\ 0 & 0 & 0 & -1 & 2 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "The matrix is **triadiagonal** and **sparse**  \n",
    "(and also **Toeplitz**: all elements on the diagonal are the same)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Block structure in 2D\n",
    "In two dimensions, we get equation of the form\n",
    "\n",
    "$$-\\frac{4u_{ij} -u_{(i-1)j} - u_{(i+1)j} - u_{i(j-1)}-u_{i(j+1)}}{h^2} = f_{ij},$$\n",
    "\n",
    "or in the **Kronecker product form**  \n",
    "\n",
    "$$\\Delta_{2D} = \\Delta_{1D} \\otimes I + I \\otimes \\Delta_{1D},$$\n",
    "\n",
    "where $\\Delta_{1D}$ is a 1D Laplacian, and $\\otimes$ is a **Kronecker product** of matrices. \n",
    "\n",
    "For matrices $A\\in\\mathbb{R}^{n\\times m}$ and $B\\in\\mathbb{R}^{l\\times k}$ its Kronecker product is defined as a block matrix of the form \n",
    "\n",
    "$$\n",
    "   A\\otimes B = \\begin{bmatrix}a_{11}B & \\dots & a_{1m}B \\\\ \\vdots & \\ddots & \\vdots \\\\ a_{n1}B & \\dots & a_{nm}B\\end{bmatrix}\\in\\mathbb{R}^{nl\\times mk}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In the block matrix form the 2D-Laplace matrix can be written in the following form:  \n",
    "\n",
    "$$A = -\\frac{1}{h^2}\\begin{bmatrix} \\Delta_1 + 2I & -I & 0 & 0 & 0\\\\ -I & \\Delta_1 + 2I  & -I & 0 &0 \\\\ 0 & -I & \\Delta_1 + 2I & -I & 0 \\\\ 0 & 0 & -I & \\Delta_1 + 2I   &-I\\\\ 0 & 0 & 0 & -I & \\Delta_1 + 2I \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Short list of Kronecker product properties\n",
    "\n",
    "- It is bilinear\n",
    "- $(A\\otimes B) (C\\otimes D) = AC \\otimes BD$\n",
    "- Let $\\mathrm{vec}(X)$ be vectorization of matrix $X$ columnwise. Then \n",
    "$\\mathrm{vec}(AXB) = (B^T \\otimes A) \\mathrm{vec}(X).$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices help in computational graph theory\n",
    "\n",
    "- Graphs are represented with adjacency matrix, which is usually sparse\n",
    "- Numerical solution of graph theory problems are based on processing of this sparse matrix\n",
    "    - Community detection and graph clustering\n",
    "    - Learning to rank\n",
    "    - Random walks\n",
    "    - Others\n",
    "- Example: probably the largest publicly available hyperlink graph consists of 3.5 billion web pages and 128 billion hyperlinks, more details see [here](http://webdatacommons.org/hyperlinkgraph/) \n",
    "- More medium scale graphs to test your algorithms are available in [Stanford Large Network Dataset Collection](https://snap.stanford.edu/data/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Florida sparse matrix collection\n",
    "\n",
    "More sparse matrices you can find in  [Florida sparse matrix collection](http://www.cise.ufl.edu/research/sparse/matrices/) which contains all sorts of matrices for different applications."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"700\"\n",
       "            height=\"450\"\n",
       "            src=\"http://yifanhu.net/GALLERY/GRAPHS/search.html\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x10815f198>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import IFrame\n",
    "IFrame(\"http://yifanhu.net/GALLERY/GRAPHS/search.html\", width=700, height=450)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices and deep learning\n",
    "\n",
    "- DNN has a lot of parameters\n",
    "- Some of them may be redundant\n",
    "- How to prune the parameters without significantly accuracy reduction?\n",
    "- [Sparse variational dropout method](https://github.com/ars-ashuha/variational-dropout-sparsifies-dnn) leads to significantly sparse filters in DNN almost without accuracy decreasing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrix: construction\n",
    "\n",
    "- We can create sparse matrix using **scipy.sparse** package (actually this is **not the best** sparse matrix package)\n",
    "\n",
    "- We can go to really large sizes (at least, to store this matrix in the memory)\n",
    "\n",
    "Please note the following functions\n",
    "- Create sparse matrices with given diagonals ```spdiags```\n",
    "- Kronecker product of sparse matrices ```kron```\n",
    "- There is also overloaded arithmectics for sparse matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0xa18218fd0>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP8AAAD8CAYAAAC4nHJkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAD8RJREFUeJzt3V+IHeUZBvDnaTYJZo0lklWDtdlWgn8QGs1iBYtGrCV6Y3shaKHNhRAvImj1wuCNUih4o7YXRRqrJBfVIlhrEKl/gtYWinRTpMZGUSRadU2OKE0IokbfXuysrsnuzrfnfGfmfed7fhB292R295lD3jO78z2ZoZlBRMrzjbYDiEg7NPwihdLwixRKwy9SKA2/SKE0/CKFam34SW4i+RrJN0huaytHCpL7Sb5M8iWSk23nORbJB0keJLl31mMnk3yG5OvV21VtZpxtnrx3kny3eo5fInlVmxlnkDyD5HMk95F8heRN1eNun99UrQw/ySUAfgvgSgDnAriO5LltZFmEy8xsvZlNtB1kDjsAbDrmsW0AdpvZOgC7q4+92IHj8wLAvdVzvN7Mnmw403yOArjVzM4BcBGArdW/Vc/Pb5K2jvwXAnjDzN40s08B/BHA1S1lCc/MXgDw4TEPXw1gZ/X+TgA/bjTUAubJ65KZTZnZv6r3DwPYB+B0OH5+U7U1/KcD+O+sj9+pHvPKADxNcg/JLW2HSXSqmU0B0/+AAZzScp4UN5L8d/Vrgbsfo0mOAzgfwIuI+fx+TVvDzzke89wzvtjMLsD0rylbSV7SdqAOug/AmQDWA5gCcHe7cb6O5IkAHgVws5kdajtPDm0N/zsAzpj18bcAvNdSllpm9l719iCAxzD9a4t3B0iuAYDq7cGW8yzIzA6Y2edm9gWA++HoOSa5FNOD/wcz+1P1cKjndy5tDf8/Aawj+R2SywBcC2BXS1kWRHKU5MqZ9wH8CMDehT/LhV0ANlfvbwbweItZas0MUuUncPIckySABwDsM7N7Zv1VqOd3Lmzrf/VVSzm/BrAEwINm9qtWgtQg+V1MH+0BYATAQ96yknwYwEYAqwEcAHAHgD8DeATAtwG8DeAaM3Nxkm2evBsx/SO/AdgP4IaZ36nbRPIHAP4G4GUAX1QP347p3/tdPr+pWht+EWmXGn4ihdLwixRKwy9SKA2/SKE0/CKFanX4A1VlAcTKGykrECtvpKwLafvIH+1JjJQ3UlYgVt5IWefV9vCLSEsaLfmsXr3axsfHv/y41+thbGysse8/qEh5I2UFYuX1nHXPnj0fmFlSuJFBvhHJTQB+g+mK7u/N7K6Fth8fH8fkpLsL4Yh0Bsm3Urft+8f+Qa/Gc+STo9jz1kc48snRgbZpiqcsIjkMcuT/8mo8AEBy5mo8/6n7xCOfHMUP7/krDn38GU46YSmeveVSjC4fWfQ2TfGURSSXQU74JV2Nh+QWkpMkJ3u9HgDg1fcP49DHn+HIp5/jfx9/hlffP3zcF0/ZpimesojkMsjwJ12Nx8y2m9mEmU3MnCQ5+7SVOOmEpVixbAm+ecJSnH3ayuO+UMo2TfGURSSXQX527ftqPKPLR/DsLZfi1fcP4+zTVs75I3TKNk3xlEUkl0H+FX95NR4A72L6ajw/Tf3k0eUj2LB24Ws0pmzTFE9ZRHLoe/jN7CjJGwE8ha+uxvNKtmQiMlQD/fxa3VjBy80VRGQRVO8VKVRrw5+r5OOlfOMlh0iqVk5b5yr5eCnfeMkhshitHPlzlXy8lG+85BBZjFYOTzOlGQNqSz6DbtMELzlEFqPR/9I7MTFhM/+r78gnR2tLM7m2aYKXHFI2kntSbyPf2vCLSH6LGX4t9YkUSsMvUigNv0ihXA9/jpKPp/KNpywibk9L5yj5eCrfeMoiAjg+8uco+Xgq33jKIgI4PvLnKPl4Kt94yiICOF/nz1Hy8VS+8ZRFukklH5FCqeQjIrU0/CKFCj/8Odb5Pa2/e8oi3Rb6rFOOdX5P6++eskj3hT7y51jn97T+7imLdF/ow0qOdX5P6++eskj3hV/qy7HO72n93VMWiUfr/CKF0jq/iNTS8IsUSsMvUqjOD3+0kk+KaHnFp06fTo5W8kkRLa/41ekjf7SST4poecWvTh8yopV8UkTLK351fp0/WsknRbS80pzFrPN3/l/O6PIRbFi7auBtPImWV3zq9O/8IjK/gY78JPcDOAzgcwBHU3/cEJH25fix/zIz+yDD1xGRBrX2Y3+0ooruDCRdM+iR3wA8TdIA/M7Mtqd8UrSiiu4MJF006JH/YjO7AMCVALaSvOTYDUhuITlJcrLX6wGIV1TRnYGkiwYafjN7r3p7EMBjAC6cY5vtZjZhZhNjY2MAviqqrFi2JERRpS6vp/3xlEV867vkQ3IUwDfM7HD1/jMAfmlmf5nvc2aXfKIVVXRnIImgqZLPqQAeIznzdR5aaPCPFa2oUpfX0/54yiJ+9T38ZvYmgO9lzCIiDVLDT6RQroff03q17gwkXeP2bJCn9WrdGUi6yO2R39N6te4MJF3k9uXe00UrdGcg6SLXF/PwtF6tOwNJBLpjj0ihdMceEaml4RcplIZfpFDhh99LWSVXycfL/njLIvmFPs3rpaySq+TjZX+8ZZHhCH3k91JWyVXy8bI/3rLIcIR+KfdSVslV8vGyP96yyHCEX+f3UlbJVfLxsj/eskgalXxECqWSj4jU0vCLFErDL1Kozg+/p6JK14pAXnJIfzp9CtdTUaVrRSAvOaR/nT7yeyqqdK0I5CWH9K/TL9WeiipdKwJ5ySH96/w6v6eiSteKQF5yyFdU8hEplEo+IlJLwy9SKA0/fK1XR1rnTxEpa2mKP0vjab060jp/ikhZS1T8kd/TenWkdf4UkbKWqPiXYU/r1ZHW+VNEyloiLfXB13p1pHX+FJGydoHW+UUKpXV+Eaml4RcpVO3wk3yQ5EGSe2c9djLJZ0i+Xr1dNdyYIpJbypF/B4BNxzy2DcBuM1sHYHf1cadFKqtEKwp5ylKS2tOvZvYCyfFjHr4awMbq/Z0AngdwW8ZcrkQqq0QrCnnKUpp+f+c/1cymAKB6e8p8G5LcQnKS5GSv1+vz27UrUlklWlHIU5bSDP2En5ltN7MJM5sYGxsb9rcbipmyyoplS9yXVVKyetofT1lKk7TOX/3Y/4SZnVd9/BqAjWY2RXINgOfN7Ky6rxN5nT9SWSVaUchTluiaWOffBWBz9f5mAI/3+XXCGF0+gg1rV4X4x5mS1dP+eMpSkpSlvocB/APAWSTfIXk9gLsAXEHydQBXVB+LSCApZ/uvm+evLs+cRUQapIafSKE0/Jl4KqrozkCSQmdYMvBUVNGdgSSVjvwZeCqq6M5AkkovoRl4umKN7gwkqXQxj0w8FVV0Z6By6Uo+IoXSlXxEpJaGX6RQGn6RQmn4G+SprJKj5BNtf+TrdOq0IZ7KKjlKPtH2R46nI39DPJVVcpR8ou2PHE8vjw3xVFbJUfKJtj9yPK3zN8hTWSVHySfa/pRAJR+RQqnkIyK1NPwihdLwO+NpvTrHOn+k/SlNuWdGHPK0Xp1jnT/S/pRIR35HPK1X51jnj7Q/JSr7pc8ZT+vVOdb5I+1PibTU54yn9eoc6/yR9qcLtM4vUiit84tILQ2/SKE0/CKF0vAH5KWsEq3kkyJa3kF085Rnh3kpq0Qr+aSIlndQOvIH46WsEq3kkyJa3kF192Wto7yUVaKVfFJEyzsorfMH5KWsEq3kkyJa3mOp5CNSKJV8RKRW7fCTfJDkQZJ7Zz12J8l3Sb5U/blquDFFJLeUI/8OAJvmePxeM1tf/XkybywRGbba4TezFwB82EAWySRaUUV3BmrHIKczbyT5cwCTAG41s48yZZIBRCuq6M5A7en3hN99AM4EsB7AFIC759uQ5BaSkyQne71en99OUkUrqujOQO3pa/jN7ICZfW5mXwC4H8CFC2y73cwmzGxibGys35ySaKaosmLZkhBFlbq8nvbHU5Ycktb5SY4DeMLMzqs+XmNmU9X7vwDwfTO7tu7raJ2/GdGKKrozUD5ZSz4kHwawEcBqAAcA3FF9vB6AAdgP4IaZF4OFaPhFhmsxw1/70mVm183x8AOLTiUirqjhJ1IoDX+hPK1X685A7fB3xkKGztN6te4M1B4d+Qvkab1adwZqj8+XJBkqTxet0J2B2qP/z18oT+vVujNQPrqYh0ihdDEPEaml4RcplIZfpFAafpmXl7JKrpKPl/3xkkVLfTInL2WVXCUfL/vjKYuO/DInL2WVXCUfL/vjKYuO/DInL2WVXCUfL/vjKYvW+WVeXoozuUo+XvZnmFlU8hEplEo+IlJLwy9SKA2/SKE0/NI3D0WVxWSJVARqIoeW+qQvXooqqVkiFYGayqEjv/TFS1ElNUukIlBTOXTkl754KaqkZolUBGoqh9b5pW/RSjORikD95lDJR6RQKvmISC0Nv0ihNPwyVF7WzVOzeMpbZ9CsOtsvQ+Nl3Tw1i6e8dXJk1ZFfhsbLunlqFk956+TI6vNlTTrBy7p5ahZPeevkyKqlPhkqL+vmqVk85a0zV1at84sUSuv8IlJLwy9SqNrhJ3kGyedI7iP5CsmbqsdPJvkMydert6uGH1dEckk58h8FcKuZnQPgIgBbSZ4LYBuA3Wa2DsDu6mORRetasSbK/tSezjSzKQBT1fuHSe4DcDqAqwFsrDbbCeB5ALcNJaV0VteKNZH2Z1G/85McB3A+gBcBnFq9MMy8QJwyz+dsITlJcrLX6w2WVjqna8WaSPuTPPwkTwTwKICbzexQ6ueZ2XYzmzCzibGxsX4ySofNlFVWLFsSplizUNZI+5O0zk9yKYAnADxlZvdUj70GYKOZTZFcA+B5Mztroa+jdX6ZS/RiTT/bDEvWdX6SBPAAgH0zg1/ZBWBz9f5mAI8vNqgIAIwuH8GGtavcDz6QljXK/qSkuxjAzwC8TPKl6rHbAdwF4BGS1wN4G8A1w4koIsOQcrb/7wA4z19fnjeOiDRFDT+RQmn4xT1PpZku3RnI9xkJKZ6n0kzX7gykI7+45qk007U7A+nIL655urpO1+4MpIt5iHueSkDe7wykK/mIFEpX8hGRWhp+kUJp+EUKpeGXTvBQmllMlrptmtgfLfVJeF5KM6lZ6rZpan905JfwvJRmUrPUbdPU/ujIL+F5Kc2kZqnbpqn90Tq/dELXikD97o9KPiKFUslHRGpp+EUKpeGXYkTqAjRxQRCd7ZciROoCNHVBEB35pQiRugBNXRBER34pQqQuQFMXBNFSnxQjUheg366A1vlFCqV1fhGppeEXKZSGX6RQGn6RWbwUgVTyEWmQlyKQSj4iDfNSBFLJR6RhXopAKvmItMBLEaiJko+O/CKzjC4fwYa1q9qOkZRj0KyNHvlJ9gC8Neuh1QA+aCzA4CLljZQViJXXc9a1ZjaWsmGjw3/cNycnU39E8SBS3khZgVh5I2VdiM72ixRKwy9SqLaHf3vL33+xIuWNlBWIlTdS1nm1+ju/iLSn7SO/iLREwy9SKA2/SKE0/CKF0vCLFOr/vZKAwUorgPwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "from scipy.sparse import csc_matrix, csr_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.linalg\n",
    "import scipy.sparse.linalg\n",
    "%matplotlib inline\n",
    "n = 5\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "plt.spy(A, aspect='equal', marker='.', markersize=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparsity pattern \n",
    "\n",
    "- The ```spy``` command plots the sparsity pattern of the matrix: the $(i, j)$ pixel is drawn, if the corresponding matrix element is non-zero.\n",
    "\n",
    "- Sparsity pattern is really important for the understanding the complexity of the sparse linear algebra algorithms.  \n",
    "\n",
    "- Often, only the sparsity pattern is needed for the analysis of \"how complex\" the matrix is."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrix: definition\n",
    "\n",
    "- The definition of \"sparse matrix\" is that the number of **non-zero elements** is much less than the total number of elements. \n",
    "\n",
    "- You can do the basic linear algebra operations (like solving linear systems at the first place) faster, than if working for with the full matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## What we need to find out to see how it actually works\n",
    "\n",
    "- **Question 1:** How to store the sparse matrix in memory?\n",
    "\n",
    "- **Question 2:** How to multiply sparse matrix by vector fast?\n",
    "\n",
    "- **Question 3:** How to solve linear systems with sparse matrices fast?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrix storage\n",
    "\n",
    "There are many storage formats, important ones:\n",
    "\n",
    "- COO (Coordinate format)\n",
    "- LIL (Lists of lists)\n",
    "- CSR (compressed sparse row)\n",
    "- CSC (compressed sparse column)\n",
    "- Block variants\n",
    "\n",
    "In ```scipy``` there are constructors for each of these formats, e.g. \n",
    "\n",
    "```scipy.sparse.lil_matrix(A)```."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Coordinate format (COO)\n",
    "\n",
    "The simplest format is to use **coordinate format** to represent the sparse matrix as positions and values of non-zero elements.\n",
    "\n",
    "```python\n",
    "i, j, val\n",
    "```\n",
    "\n",
    "where ```i, j``` are integer array of indices, ```val``` is the real array of matrix elements. <br>\n",
    "So we need to store $3\\cdot$ **nnz** elements, where **nnz** denotes number of nonzeroes in the matrix.\n",
    "\n",
    "**Q:** What is good and what is bad in this format?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Main disadvantages\n",
    "\n",
    "- It is not optimal in storage (why?)\n",
    "- It is not optimal for matrix-by-vector product (why?)\n",
    "- It is not optimal for removing elements as you must make **nnz** operations to find one element (this is good in LIL format)\n",
    "\n",
    "First two disadvantages are solved by **compressed sparse row** (CSR) format."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Compressed sparse row (CSR)\n",
    "\n",
    "In the CSR format a matrix is stored as 3 different arrays:\n",
    "\n",
    "```python\n",
    "ia, ja, sa\n",
    "```\n",
    "\n",
    "where:\n",
    "\n",
    "- **ia** (row start) is an integer array of length $n+1$ \n",
    "- **ja** (column indices) is an integer array of length **nnz** \n",
    "-  **sa** (values) is an real-value array of length **nnz**\n",
    "\n",
    "<img src=\"https://www.karlrupp.net/wp-content/uploads/2016/02/csr_storage_sparse_marix.png\" width=60% />\n",
    "\n",
    "So, we got $2\\cdot{\\bf nnz} + n+1$ elements."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices in PyTorch and Tensorflow\n",
    "\n",
    "- PyTorch supports sparse matrices stored in COO format\n",
    "- Incompletre backward operation for these matrices, see summary [here](https://github.com/pytorch/pytorch/issues/9674)\n",
    "- Tensorflow also supports sparse matrices in COO format\n",
    "- The list of supported operations is [here](https://www.tensorflow.org/api_docs/python/tf/sparse) and gradient support is also limited"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### CSR helps in sparse matrix by vector product (SpMV)\n",
    "\n",
    "```python\n",
    "\n",
    "   for i in range(n):\n",
    "        \n",
    "        for k in range(ia[i]:ia[i+1]):\n",
    "            \n",
    "            y[i] += sa[k] * x[ja[k]]\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let us do a short timing test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "7.91 ms ± 453 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "10.8 ms ± 536 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg\n",
    "from scipy.sparse import csc_matrix, csr_matrix, coo_matrix\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "n = 1000\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "B = coo_matrix(A)\n",
    "%timeit A.dot(rhs)\n",
    "%timeit B.dot(rhs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "As you see, **CSR** is faster, and for more **unstructured patterns** the gain will be larger."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse matrices and efficiency\n",
    "\n",
    "- Sparse matrices give complexity reduction.  \n",
    "- But they are **not very good** for parallel/GPU implementation.  \n",
    "- They do not give maximal efficiency due to **random data access**.  \n",
    "- Typically, peak efficiency of $10\\%-15\\%$ is considered good.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Recall how we measure efficiency of linear algebra operations\n",
    "\n",
    "The standard way to measure the efficiency of linear algebra operations on a particular computing architecture is to use **flops** (number of floating point operations per second)\n",
    "\n",
    "We can measure peak efficiency of an ordinary matrix-by-vector product."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  7.3e-03, Efficiency:  4.4e+00 Gflops\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "n = 4000\n",
    "a = np.random.randn(n, n)\n",
    "v = np.random.randn(n)\n",
    "t = time.time()\n",
    "np.dot(a, v)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  ((2 * n ** 2)/t) / 10 ** 9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  1.4e-04, Efficiency:  8.4e-03 Gflops\n"
     ]
    }
   ],
   "source": [
    "n = 400\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.random.randn(n)\n",
    "t = time.time()\n",
    "a.dot(rhs)\n",
    "t = time.time() - t\n",
    "print('Time: {0: 3.1e}, Efficiency: {1: 3.1e} Gflops'.\\\n",
    "      format(t,  (3 * n) / t / 10 ** 9))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Random data access and cache misses\n",
    "\n",
    "- Initially all matrix and vector entries are stored in RAM (Random Access Memory)\n",
    "- If you want to compute matvec, the part of matrix and vector elements are moved to **cache** (fast and small capacity memory, see [lecture about Strassen algorithm](https://nbviewer.jupyter.org/github/oseledets/nla2018/blob/master/lectures/lecture-3.ipynb))\n",
    "- After that, CPU takes data from cache to proccess it and then returns result in cache, too\n",
    "- If CPU needs some data that is not already in cache, this situation is called **cache miss**\n",
    "- If cache miss happens, the required data is moved from RAM to cache\n",
    "\n",
    "**Q**: what if cache does not have free space?\n",
    "\n",
    "\n",
    "- The larger number of cache misses, the slower computations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Cache scheme and LRU\n",
    "\n",
    "<img src=\"cache_scheme.png\" width=\"500\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### CSR sparse matrix by vector product\n",
    "\n",
    "```python\n",
    "\n",
    "   for i in range(n):\n",
    "        \n",
    "        for k in range(ia[i]:ia[i+1]):\n",
    "            \n",
    "            y[i] += sa[k] * x[ja[k]]\n",
    "            \n",
    "```\n",
    "\n",
    "- What part of operands is strongly led cache misses?  \n",
    "- How this issue can be solved?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Reordering reduces cache misses\n",
    "\n",
    "- If ```ja``` stores sequential elements, then they will be moved to cache altogether and number of cache misses decreases\n",
    "- This happens when sparse matrix is banded or at least block diagonal\n",
    "- We can convert given sparse matrix to banded or block diagonal with *permutations* \n",
    "\n",
    "- Let $P$ be row permutation matrix and $Q$ be column permutation matrix\n",
    "- $A_1 = PAQ$ is a matrix, which has less bandwith than $A$\n",
    "- $y = Ax \\to \\tilde{y} = A_1 \\tilde{x}$, where $\\tilde{x} = Q^{\\top}x$ and $\\tilde{y} = Py$\n",
    "- [Separated block diagonal form](http://albert-jan.yzelman.net/PDFs/yzelman09-rev.pdf) is a cache-oblivious format for sparse matrix by vector product\n",
    "- It can be extended for 2D, where separated not only rows, but also columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Example\n",
    "\n",
    "- SBD in 1D \n",
    "<img src=\"sbd.png\" width=\"400\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Sparse *transpose* matrix-by-vector product\n",
    "\n",
    "- In some cases it is important to compute not only $Ax$ for sparse $A$, but also $A^{\\top}x$\n",
    "- Mort details will be discussed in the lecture about Krylov methods for non-symmetric linear systems\n",
    "- Transposing is computationally expensive\n",
    "- [Here](https://people.eecs.berkeley.edu/~aydin/csb2009.pdf) is proposed **compressed sparse block** format of storage appropriate for this case"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Compressed sparse block (CSB)\n",
    "\n",
    "- Split matrix in blocks\n",
    "- Store block indices and indices of data **inside** each block\n",
    "- Thus, feasible number of bits to store indices\n",
    "- Ordering of the blocks and inside blocks is impoprtant for parallel implementation\n",
    "- Switching between blockrow to blockcolumn makes this format appropriate to transpose matrix by vector product"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Solve linear systems with sparse matrices\n",
    "\n",
    "- Direct methods (start today and continue in the next lecture)\n",
    "    - LU decomposition\n",
    "    - Number of reordering techniques to minimize fill-in\n",
    "- Krylov methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Let us start with small demo of solving sparse linear system..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5,1,'Reshaped solution')"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEICAYAAABLdt/UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnXm8JGdZ779P72efJbNmMpkskAARAoxhUZRLAgLGoF5UEBH0cnNdUPCqKBdF3L3IVfDjwo2IXmURZZWoQCIGXNiSEEJghpBAkpnMdmY9a+/P/aOquut019pd3V3d5/1+Pudzzumu6nq7uvr3PvV7n/d5RVUxGAwGw+SQGXUDDAaDwZAsRtgNBoNhwjDCbjAYDBOGEXaDwWCYMIywGwwGw4RhhN1gMBgmDCPsCSIirxSRfx91OwBEREXkygRf74D9mrke93+ZiHwiqfYY+kdEHhKRG1LQjjeJyLsSfs2/EpHf6mP/r4jIsxNs0lBJvbDbF99JEZlxPfYqEbkj4v53iMirBtZAQxdenYCqvltVnzfKdo0r9ndgXURWROSELVqzo27XpODVCajqE1T1jhE1qW9SL+w2OeA1wzhQrxGpwTBgvkdVZ4FrgScDrx9xewwpZlyE/feBXxCRLV5PisgzReQLInLB/v1M+/HfBp4F/LEd7fyxx75OdPnfROQR4JP2408Xkf8UkfMi8iX3bZltuXxDRJZF5Jsi8rKO13yLiJyzn3uB6/EfE5FD9n7fEJH/4Xru2SJyVET+l4ictqO0l7meL9qv+4h9B/N2EZlyPf+LInJcRI6JyI8HnUy/9otIRkR+RUQeFpFTIvLXIrLg8xobbuM7bqc/bf8+b5/3Z3TaVH6fmf3cHSLymyLyH3YbPyEiFwW9p82Cqp4APo4l8EDwtSEiF4nIrfZ1fFZE/k1E3N/7a0XkXvtzeJ+IlOz9ttr7LdrX8q0iss91zDtE5HdF5PP2vh8RkW2u54O+P5eJyKfsz/Y2wPezDWq/iDzObsd5sayTm3xeo8sitb/zV4rIzcDLgNfZ1+pH7edb17d9ft9qf7eO2X8X7eec7+3P29+Z4yLyY4Ef4jBQ1VT/AA8BNwAfBH7LfuxVwB3239uAc8DLsSL7l9r/b7efvwN4VcDrHwAU+GtgBpgCLgbOAC/E6vyea/+/w95mCbjK3n8P8AT771cCNeC/A1ngJ4FjgNjPfzdwBSDAdwJrwFPs554N1IE/AIr286uu47wV+Af7/c4BHwV+137u+cBJ4Bq7fe+x39OVHu83qP0/DjwAXA7M2uf8bzrOU879ubhe903Au7y2dZ2bf4/xmT0IPNb+PO4Afm/U1+KovwP23/uALwNvcz0fdG38LvB2IG//PMt1PT4EfB7Ya+97CPgJ+7ntwH8Fpu3X/Hvgw65j3gE86rrmPuD6/H2/P/bzn3Fd598BLDv7erx3z/bbfz8A/C+gADzHfh3nuv4r2nrRuvZcr9v6fri39TnnvwF8FtiJpQH/Cfxmx/f2N+w2vRDre711pNfMqC/aqBe1fQFdsE+sW9hfDny+Y5/PAK90XYBRhP1y12O/hC1orsc+DrzCvojP2xf9VMc2rwQecP0/bb/2bp9jfxh4TccFMuN6/u+AX7Uv5FXgCtdzzwC+af/9TlzChyWIQcLu1/5/AX7K9f9VWB1VjmSFPcpn9iuu534K+Nior8URfwdWsIRL7c9pi/1c2LXxG8BHfK6Fh4Afcf3/ZuDtPm24Fjjn+v+Ojmvu8UAVK6AJ+v7s97jO34O/sHu2H0vgTwAZ12PvBd5k//1XJCfsDwIvdD33XcBD9t/PBtY7rvVTwNNHec2MixWDqt4H3Ar8csdTe4GHOx57GCtqiMMR19+XAj9g3+KdF5HzwLcDe1R1Ffgh4CeA4yLyjyJytWvfE642r9l/zgKIyAtE5LP2LeV5rN7dfRt6zn599/vYi9WZTQN3udrzMftx5xwc6djPk5D2d57Lh7FEfZff6/VIlM/shOvvNexzuIn5XlWdwxKSq2lfN2HXxu9jRbafsO23zu+P53kWkWkR+b+2LbeEZa9tEZGsa/vOay5vt8v3+4P12Xtd5374tX8vcERVmx2vE/d7HwWv78Ve1/9nVLXu+n/k1+vYCLvNr2HZHO4P7xjWheRmP9ZtIlg9cxTc2x3Biji2uH5mVPX3AFT146r6XKwL9TDw52EvbntyHwDeAuxS1S3AP2FFXA5bxZX9Y7+PY8BprKjgCa72LKg1mAZwHLikYz//N+rf/s5z6URXJz1eZhVLUBx2uw8RdHyP4zjHetRjW4MLVf0UVoT5FvuhwGtDVZdV9edV9XLge4D/KSLXRzjUz2PdsT1NVeexLBPYeL12XnM1uz1B35/jeF/nfu/Xr/3HgEs6xgv8rqEN16qI7O54Pu716nwvU8tYCbuqPgC8D/hZ18P/BDxWRH5YRHIi8kNYt4W32s+fxPKM4/Au4HtE5LtEJCsiJXuQZJ+I7BKRm+wLs4J1i9yI8JoFLE9xEaiLNajqlf736yJSEJFnATcCf29HJX8O/KGI7AQQkYtF5Lvsff4OeKWIPF5EprE6QE9C2v9e4Ofswa1Z4HeA93VEIw73AC8RkbyIHARe7HpuEWjif97DPjNDMG8Fnisi14ZdGyJyoz1IKFhjKw2iXa9zWB3GeXtQ1Oua+hHXNfcbwPtVtUHA90dVHwbupH2dfzuWYHsS0P7PYQn26+xr8Nn26/ytx8t8CXiCiFwr1uDwmzqeD9OI9wK/IiI7xBrEf6P9HlPLWAm7zW9g+cQAqOoZLAH8eawBmtcBN6rqaXuTtwEvFmtk/4+iHEBVjwAvwhqYWcSKQH4R63xl7GMdA85iDXL+VITXXMbqkP4Oa6Dwh7EGvNycsJ87BrwbayDrsP3cL2Hdkn7WvjW+HSuiQlX/GevL/kl7m08GNCWo/e8E/gbrtvubQBn4GZ/X+VWsgeBzwK9j+aTOe10Dfhv4D/tW/Okd5yLsMzMEoKqLWIP9v2o/5HttAI+x/1/BGsf4U42Wn/1WrIHr01gDhx/z2OZvsO4eTgAl7IAr5PsD1rX/NKzr79fs9+KHZ/tVtQrcBLzAbuOfAj/q+r60UNX7sXTjduDrQOckwr8AHm9fqx/2aMNvYXVG92INXN9tP5ZanNFxw4ixI453qeq+sG0NhlEj1gTBd6nqO0bdFkM34xixGwwGgyEAI+wGQwAi8hoRuU+sCTCvHXV7DIYo9GXFiMgPYA1EPA64TlXvTKhdBsPIEZFrsAbjrsPK0f4Y8JOq+vWRNsxgCKHfiP0+4PtpTyE3GCaJxwGfVdU1OzPoU8D3jbhNBkMofRW8UtVDAFYmUnQuuugiPXDgQD+HNhh8ueuuu06r6o7wLUO5D/htEdmOlfr3QqzsiA2IVW/kZoDpaXnqZVfE+1qda0yHbzQkluqlUTdhYMznyqNuwga2ZtfCN+rgK1+uRbq2h1bJ0H3x79+/nzvvNK6NYTCISNBMxsio6iER+d/AbVjpdl/CmrDVud0twC0A1zyxoB/4x3j1yt6/9JT+G9sHt524OnyjCeO5u7uyIofOi+fvjr3P1fuPR7q2Q4VdRG5n46xChzeo6keiNsh98R88eNDkWBrGAlX9C6w8Z0Tkd4CjSb7+KER9Mwp5J53nYBRC//6lp/Qk7lEIFXZVHfkKKwbDqBCRnap6SkT2Y40nPSOp1x6WqBshD2dUQj8ocTeLShgMwXzA9thrwE+r6rlRNygMI+T9k4aIvh/6yooRke8TkaNYUcw/isjHk2mWwZAOVPVZqvp4VX2Sqv5LUq87qGjdiPpguO3E1QM7t4O4FvoSdlX9kKruU9Wiqu5S1e8K38tg2NwM4os8SOExtBkXcTdWjMEwxhgxHz7OOU+zPWNKChgMQyTJyMyI+mhJ+i4pyWsjVRH7O/7tG+zbOs3zr/HKruymUm/wkXuO8QNP3Rd7kpQhPXz4i49y/eN2MlfKR9r+0/cvcs+R8/zMc64cq889qS9umgX9oaNJzAvz5sC+xYG9dj8kGcEnlSWTKmH/f595iG+9dFtkYf/Xw4u87v338rjd83zLvoXBNs4wEI6eW+O177uH3/3+b+Gl1wUu/NTiU/cv8r4vHOFnr3/MgFuXLtIg6IMU7n6PPWrhv+3E1amxZ1Il7HPFPEvlWuTtL6xXAThv/zaMH+fXaht+R2FpvcZcKVWXbij9ROujEPRRCnivBLV5WKKfRPSeRNSeqm/H/FSOpXWvVdi8cbaNs48hXTgdeZwOfalcYz6ibZMGxkXUx1HMo+J+b8MQ+X4Fvl9xT5ewl/I8cjZ6YZxeRMGQLtqdc5yIvc78VKou3cQZlqBPspj7MUyRH5U9k6pvx/xUPuYXvLbht2H8aHfOMe7UyjV2z49HFcJeovVBi/pmFHM/hiHyvUbv/UTt6RL2Uj7mF9yO9kzEPrb00jkvlWs8dtfcoJqUGGkSdSPm4QzbrolCr+Keqjz2+akcK5U6jWa04o9tUTAe+7jSS+e8XK4zn/LB0zSI+kNHd7R+DPEYxLnr9fPt5VpKlbA7ecwrEaP2ZROxjz1O57wc8TNXVTsrZnwGT4eNEfNkSfJ8DmvsJFXC7kRhUYXa2S6qKBjSR6tzjmjFrFYbNJVUD56OKlo3gj5YRnlu415Tqfp2zE9ZUdiF9RqXRNjeDJ6OP3Ezm5zPOq3pjr0sc5eUqKeJ4iOFxF6rsj8981QeOrqjb/+910wZS9z/MdK26RJ2+8saPWI3Vsy44wh1udakUm9QzGWDt7c/aycIGHfGXdSTFPC4xxiV4Dvnux+BH3QaZLqE3b69jjIYWm80WamYCUrjjjsLarlcpzgbIuz2Zz2siF1Efg54FaDAl4EfU9VUrIo8TEEfhoDHZdSC32/0PkhxT5nHHj1id0Q9nxUTsY8xS+s18llp/R1lexiOxy4iFwM/CxxU1WuALPCSpF6/n2h90KJefKSw4Wec6Gz7INufNgvMIV3Cbt9eRxkMdba5eMsUa9UGtUZzoG0zDIalco2Lt0wBET/3ytA99hwwJSI5YBo4lsSLplHUx1XIozDI99XPoPWgsmRSJeyzRceKCY/cLtjbXLw1uigY0kWzqaxU6q3PMMqdl2PFDKMImKo+CrwFeAQ4DlxQ1U8M/MA+DCrrZVLF3ItBdl5pEvd+1zz9fRE5LCL3isiHRGRLP6+XzQhzxVy0L7i9zb4tVhaCyYwZP5YrdVTdn2F45+x8zsPIYxeRrcCLgMuAvcCMiPyIx3Y3i8idInLn6rlwf7eXL3LSgj7J0XlUBvH+02LN9Bux3wZco6pPBO4HXt9vg6x6MVG+4LYVYyL2sWW5vPGuK2qHPpXPUsgN5WbzBuCbqrqoqjXgg8AzOzdS1VtU9aCqHpzZmp5I0I/NLOZeTKK497uY9SdU1VHUzwL7+m3QXClmxB5DFAzpwumcW59hpMHToVZ2fAR4uohMi7VU0/XAoWEdPGk2e4QeRNLnJq64J23HJBn2/Djwz35Pum9XFxf9U4SiVnh0ttm31Vgx44rTGe+eL5HNRMtuGmYtdlX9HPB+4G6sVMcMcEs/rxn3C5xU9GcEPRqTcp5CQx8RuR3wWqvuDar6EXubNwB14N1+r6Oqt2B/KQ4ePOhb5Wu+lOfR8+thzWKpXEcE9m4p2f8bYR832qmLeeZL0RZZWSrXhjo5SVV/Dfi1oR3QRRKiPilCNUycc9ZvPnwSs1R7JVTYVfWGoOdF5BXAjcD1qhqtLGMA86UchyOI9HK5xmwxx5Zp60Mwk5TGD2dy0sJUnvmpfMtzD2K5XGfbjBGrKBhR74/iI4WhinuSE5b6zYp5PvBLwE2qGn3powCiWzF15kt5ZgpZMmIi9nHEXfclai3+pfXxWhbPTRwbpt9o3Yh6MozruES/o1B/DBSB26yxJT6rqj/RzwvOl3Is2zXZsxnx3e6CvaCxiDBXirfy0qHjS7z6PXfz/p94JltN9JcI5VqD7/vT/+SNNz6eZ1yxPdI+Tmc8W8oxP5VrzU0IwvncJ5lxFvW5h/u+afdl+VJ/PUgzo7Bk+s2KuVJVL1HVa+2fvkQdYPtsEVU4txZ8C3RmtcKOuSIA04Usa9VG5GMcOr7Eg4ur3H9yua+2GtocObvGoeNLfOXYhcj7rFcbFHMZshlh+0yRMyuVwO1rjSbn12tcNFvst7kTyyhEfe5hbf1MwnG86Pe8DjsFMlUzT4GWWJ9aCv6Sn1qqsMP+gpfyWcr16CUFyjVr21PLwccwRMc5l5VYn0ODUt4q+rVjrhj6eZxZqaLavkYmkTTkQEdhlCLbefxhtWEYnWZSaY+pu6fd6Qj7cpnHM++5jaqyuFxhx7y1bTGXoVyLHrE72y4aYU8M51zG+xyalPJWbLFzrshatcFKpd4qLdHJqeVya9txYxgr5wxSeEYl4FFxt2+Qlk0SA6rDIIXCbqUvBkVvF9ZrVBvN1ralfDaeoNQboccwxMMR3bifgxOx75x37tTKzO6Y9T6GfRe3c77UT1MnkkGIetrF3I/OdqfFmx+m1546K8b5ggdF044gO5FbKZ+hUuvFiklFWe2JwBHdcqzPoUHJXlgjSofe+blPGr3aMEmL+igtlkGQ9HsZtCWTxN1d6oS9lM8yV8pxaslfdB0R2THn9tijR4oVY8Ukzqk+rZjW2EqgsFvXxLgNng5rAeMkmCRBdzNu4t4vqRN2sCKyxYAMicWVjV5rKRfTirG3DRugNUSnZcXEHDwtOlbMXPid2uJyhW0zhWEVABsLkhSYSRV1hzTciQxrcDyV35Cdc6VA0e30Wkv5TCwLwMncMFZMcrSyYmJ57M2Wx74wlaeQywR+JqeWK8aGcWFEvTeSeq9pjtpTKexhqW+nlitMF7Kt7InYg6f2tufWalRjRJgGfxYdjz3G+azUGpTs6FtE2DFbbL2OF6eWKxOd6jgK0hDFjoK0i3u/9l0qhX3nXJFTy2X8Ss90Rm7xhb0tPkGWjyEa69UGy/YatHE7WCdiB2vgPKhDX1wqj52wD8pfT0JQNqOguxnV+x+GHZNOYZ8vUq41W2LRyamOL3gxn4nn7boGWs0Aav+4z2EsK8Y1eAqwY7boa8WoKosrlVb2zCQR94tuRD05krhjSaMlk05hd1LffG7LF5c3fsFLuSzVepNmM9oHVK41mLIjxaDsG0M0HDGeymfjpTvWo0fs59Zq1Bo6sR77sNis1ksYaRT3fu72UirswRkSix1eqyMOUaezl2tN9m+zFugwk5T6xzmH+7dNx0o77bJi5kqcX6tR8XgN51pw5jmMA1G+mMOM1o2gBzNMcR+0HZO6mafgmoXocVvu+LnuL7hzO1+uNZgqZLv26aRca3Dx1inuP7UcS9gvrNf40N1H+dFnHCATUHlynFFV/vYLR7j+6p2RZ3g6dz2XbJvmy4+ej3yccq3ZGjyFdod+eqXKxVumNh6jVU5geFaMiFwFvM/10OXAG1X1rUNrhItxE/WFB/sPmi5cMdyO3DlPaZmt2iupjNh3L1hf6qPnuldSOnrOKvu+Z8FlxdhRX9RosVJvMlPMMV/KcyGkiqSbTx4+yZs++lXui1HBcNx49Pw6r//gl/novccj73PeLre7e6EY2Ypx7q6Kroh9jy3mR892l/Z3rgX35z5oVPVrTuVS4KnAGvChoTXARa+iPmzrZeHBSusnja8XlV7PWVqi9lQK+2wxx+75Eg+eWul67gH7sSt3zLUea0fsUa0YK82ukMtQbcQvRfCAR7smha/b7y1Odku13iSXEWYKucj7OSUg3FbMFTtmAHhg0ftzL+UzXZH8ELkeeFBVH07qBQd9Oz4sQR+W+A5b5Ich7mH06rOn0ooBeMyuWd8vOMAVO2dajzn1RqKKiuPtFrKZWGVmnZz3SRZ2pzONk99frTcp5DIU81kq9Saqir3wii/O3ZU7K2bvwhTThazn+X3g1ApX7JgdpQX2EuC9Xk+IyM3AzQALe6yOJ8k0x16EYhiiPuwo2u/Yg7Rr5h7WsbRlUhmxA1yxY5YHTq105bI/sLjCxVummC60+6SWFRNZ2K00u2IuE1vAoB3VTiJfP2kLe4w7mWrDEnZHpKN0ls5n5XTKAJmMtD73Th44tcKVO72rPg4aESkANwF/7/W8qt6iqgdV9eDM1vSlviXNKKyRIAbdlnEcdO53zdPfFJF7ReQeEfmEiOxNqmFX7pxlrdrg2IWNA6hfP9n9BS/GsGJUtZVmV4gr7LbYeVlEk4JzlxQ7Ys9mYt05lT2sGLA+905hX63UefT8Olf6lPMdAi8A7lbVk6NqQBwGJURpE3Q3aWtb1LusQdlx/Ubsv6+qT7QHl24F3phAmwB4jC3e7i95s6l84/RK6zmHOIOn1UYTVWuffDaex+5Eog+fXZvIUgSq2jrfcYU9n8247pxiROz5jZfglTtnOX6hzIprcto3FlcBy54bES/Fx4YZNGmY/JI20QxiUG0dt6i93zVPl1z/zgCJvfsrPYT90fPrlGvNrojdiRSjzHp0RKfoDJ72YMU0mspDZ1Yj7zcunF6pthaUjnNeKo0mRZcVEy1id4S9O2KHjXdFDywub3humIjINPBc4INJvu6gIrWkBWhcBL2TceqMBkHfHruI/LaIHAFeRoIR+/bZIlun8xuEvZUR0xWxR7diKi5BKWR7E3Zoe9GTxNdPtRf3juWx1x2PPfqdk1MCwitiB7o+91xGuHT7DMNGVddUdbuqTm6OqweTIoyjfA+jvNsKFXYRuV1E7vP4eRGAqr5BVS8B3g28OuB1bhaRO0XkzsXFaMtDXblzlvtPtsXmaye9I7dijMFTt7cbN92x2mgwXcgiAg96ZOz4cd+jF3jJLZ9hqVyLvE+/VOoNfvSdn+c/HzgdeZ8HbctjtpjrKSsmTgfrfFbF3MaI/dJt0+SzsuFzv//kCpdunyafTe1Y/0BIgw0zCSQl7qOyY3rJsAr9pqjqDap6jcfPRzo2fQ/wXwNep5U5sGNHtNvQp122nXuOnOe0XYHxXw6d5Kpdc2yZ3njBO7MXIwm7K80urhVTqyvzpTwzhRxL69FF+t6jF/jsN87ysftORN6nX/7t/tN8+v5Fvngk2kxQoPWedi+U4kfssQdPvT32XDbDUy/dyu2HTqKqrFcb/McDp3na5dsjt2ezkpTwTEKkvtnpNyvmMa5/bwIO99ecjdz4pD00mso/33eC4xfW+cJD57jxiXu6tmtbAPHS7OJH7FZkms9KTOGzjnlrjNmc/XLrvcfsY8e3mmYK2XgdXqOdxw7RhL3SGuvoLgFx4xP38uDiKodPLPPJw6dYqzY8P/e0Mk5L4XUyqaI+qvc1qruufu9tf8+2Ze4Fnge8JoE2tbhq1xxX7pzl1i8d4x9tUbzxSd0ZlXHy2N1WTLEHj90S9gy1GMJea1iR1H88cJozQ6j/Xq41uO2rJ+1jxxPofFYo5uIJe2ceeyQrpu49eArwgmt2k80IH/3SMW699xg75oo87bLNFbEbGyZ5khD3QdgxgxhI72vmqar6Wi9JICLc+MQ9vO1fvs7hE8tcc/E8l13UPYCWzQj5rMROs4trxVRsyyGfzVCtR/+Anei+0VQ+/pWT/PDT9kfa7yvHLvDGj3yFRlN5/QuujmxH3PG1RVar1vuML+zWeVmPWVKg4Ep39KrO2ImfFQPWwPkzr9jOX3/mYSr1Bj983X6yE1p0LSmSEJxJjdY3I6ktKeDw0uv2c/j4sjUY+IwDvttFXdDanWbXixWTt9Mk4womWB3QkXPdBa78uPOhc9z18DkA/vPBM5GF/YhdRCuXkZiWipLPWlbThfU+smL6mKDk8DPPeQz/91MPkskIr3jmgchtGQeGtaBxHDaLqC88WOm7BEHcMgPFRwpU9kcvNpgEqRf2XfMl3v7yp4ZuZ9UpiZdmFz/dsUHRFr64wp7NCMVchlpM79rr79B22ttOFbJUG/HuLJyIPfadTC7jGsSOdudk3Wl5u4HXXbaN6y7bFrkNk0RcG2bcJs+MmiTEPe1MTP5YKZ+JnWbXywSlXj32fNYSsXrEVZ6c/cCK9OPsV7f3my5k47Wz3qSQFQq5bOw7mWIPEbu7FrthdGyWaH0zMTHfrKgLWrsnKDlCG3VJvXZWTCZeJGxPuY9bwsAR5el8/CyVbEYo5WMKu2M19TCoHLukQMeyeIbRsFlFvd/3Hfcuqd/B8LiZVhMk7JmY3q5lOUD0WZbOIGEhG99SKTgWTsz9MmIVOYsr0Dnb5oh/Z2Gdl7jljAvZTHsQO+LgqRH2/unHhtmsou6Qpvef9LjL5Ah7LtpCyu7B02IueplZcFkxufgeuxOxx7ViWvvFuEOoNbTn7J181ilnHCMrxr6TgeiD2JVas1WV09DGpDkOlzSJe5JMzDerlM9GrFHSHrRrRewRhb3W0N499pzEntjkdAi5HgZrLUsl/n6Wx55p+fthNJpKo6mtc1nMR+9gSx6Tk8adpXrw0n1pyYiZVEEbJmketJ4gYY86eNoetCvYGRlRxa9S79Fjd0Xsca0YZ9A1bofQmxVjZ8XEOJ7z+q2IPZ+JVmWz3vDMYTdEp1dhMaK+kUk8H6lPd4xKMZ+NWLa37e3Gjdir9UbbY4+dbdK7FVPowYppdSSx2tn22J1IPGxikGNjOZ1k5DunWtN47BNO4fDRRF6nevW+RF4niF5TINO6dN7EhEzRJyg1u4U96uCpq1ZMbx57b/v1YsVYYwE93FnYdyQQrcNztim6IvbIVowR9p5Jc7ReOHw0MVF3Xm8SGOb4yeQIez4TrQhYvdEatCvEEDBnu0JPloqVx56LmUbYrxVT6CELx/HYIaKwd1oxMWYAj4MVIyJbROT9InJYRA6JyDMGdaxBf/GHJeqDet1BC/yoLZmw8Zc4KY/p/2ZFJE4euzNoV4iRFVNvNGkqvUfCjqUSc6JRr1kxPVkxrpmnAJVGhKUG650ee4w7p/EYPH0b8DFVvRp4EnBoxO1JJcMQXuc4g6QXcU/jIOrEeOxTtqCoKiL+npdlxdgRe4+RaWyPvdFktpgjlxGWyvFmdOZsC6cSwd5wHy/fa/ZONkOxByumkLVEupTPcnY1vC5Gpd5olflNKyIyD3wH8EoAVa0CPRf9SCojphchGVRhkPWzAAAgAElEQVQ0OgqbxDnmMLz3cWViIvbpYpamhs96XKnUmSla/VkxhsdedQ0S9u6xx7diCtneslva+8WfIRurw+uI2GeKWVar9aBdAOtzmC2mW9iBy4FF4C9F5Isi8g4R6Sov6l4drHY+epG3YTFJoj6M44/akkmCiRH2OVus3avbe7FaqTNXsrbtZZCwpzz2em+1YuobLJV4++UyGQq5+HnzhZzEGlSu2nZNPmvdJc0Wc6yGfAb1RpNyrclsMR+5bSMiBzwF+DNVfTKwCvxy50bu1cHyW6aH3caRMGpRd0hLO9LGxAj7TERhX6nUmSlY28aJTCtdwq6oRhNbSzCzsSN9y4qRnvbL92gZOWMBEP+8gCXsy+WwztXqDGbSH7EfBY6q6ufs/9+PJfSGFDEIcY8btUe1x4aVGTMxwj5rC3tYtLhSqTNrR+wtAYsySNhop/U5IhY1iq66slviT1DqbRC00NPxtGcrpugS9kq9GdjeFduqce6c0oqqngCOiMhV9kPXA18dYZNiMwhbwUTJ6WfihD0oWlRVViv11rZxBKw1w9L22N2PRdm3kM2Qy2aoxbRiCn1YMfkYpQFgYz125//QfToGT2cidLAr9mfkbJtyfgZ4t73847XA74yyMaPOwEirqKe1XaMiEWEXkV8QERWRi5J4vV5wovAgQVmvNWgq3cIeQfwcAXMiaIgj7M4M0viDrr1YMe6smGqjGckyUtUNkT5E7fCs125ZMaVwS8x5bnYMhF1V77H98yeq6veq6rlBHMcU/+qfpMV9nAdR+xZ2EbkEeC7wSP/N6Z0oHntnpFi0o8xeBk8h+ozVWr23WjHVHq0Yx/op2HcWUQZsG01FlVZ1R4iaBmrZWG6PHSZH2A1tTFQ8PiQRsf8h8DpgpPeIUbJinOccb7fXtL528bAYHnvOmnkaxxpp13GPv18+E+/Ownn9fK5HK6ZD2IPunJznZlPusY87SUac4yLq49LOQdOXsIvITcCjqvqlCNu2cn0XFxf7OawnkSJ2+7mesmJcE5TyOdtjj1zu11mgQ6g1o1kjYHvlWYldK6beKhNsC3uEmuyOiMfNinHn90P7cwga62jdORWMsI8D4yaW49ZeN0lNYgsVdhG5XUTu8/h5EfAG4I1RDuTO9d2xI/ma1NOFLCIhg3YdkWI2I2QzEi0rZsMEpeiRcKOpNLXtzataj0XBnX4Y34qxBk+d/6McC4hdK6Yz3XGuNdbhf04775wMhqRJStzj3PWMemDbTeg3S1Vv8HpcRL4FuAz4kj2Ffx9wt4hcZ6eJDRURYbYQnEPtRIpubzfq+p7utL44HnvNFQm3H1OilElxr5XatDuEsDK6zjGdwVp3G6K2M94EpY3pju07p5rvPq07J+OxD4ykbJhxjn43Mz1bMar6ZVXdqaoHVPUA1mSOp4xC1B1mS7nYg3aFXDxhj+uxty0OaadJNqNF3/WmUxUyXnplvdGe6Rp1P8eu6TWP3TlWlLTTlUp9QwdpiMawI8JxF/Vxb38/TNQ3K2w6u9egXdSSuO4iYPEE09UhOBObYnjzbs87yvGaTbU7hHjtbHVAro4rStXLar3ZsrTAPXgabMUYG8YwDDaruCcm7Hbkfjqp1+uFmWJwxL7sEbEXc5nIAgYdE5Ri5HnnsxlymeiRvpVX3inQ4fs5dwPu/aIsaL3BY485eFpwRd7ZjDCVzwZbMeX6prNh0rLWaVQ2qyB6MY757BMVsc+FWDGrlTq5jLT8YIhvxcRNB3R713FmrDq5524rph5lv0Z7v0IuvsdeyGXIZIRcJloBMWdVKTeWJeYfsbtn/xqSZxyFaDMxjMloEyXsM4Vca4DUCydSdNdrj5px0rJiYkbQbo+9XWMmbofQa0fS2/HAEvhodyQewh7hzmmzReyG0bEZ7z4mSthnS8Ee+0ql0RUpxh487bVeeUwrZsNgZg8dSS5mh1Cpdwt71P0K2W5hDxvrmDPCnlo2oxAmRVpSHifq2zVbzLV8dC9WKjVvYY8YseezQiYTM9vEFQk37YlJkfZrtiP9XqyYwoasmCgee7sjsfaP3uEVOyL2mWI2+M7JROwGw0CZrIjdjhT9ZnauVhpd09jjCJgTmTpeeSxrJNeeMDQsK6YV6UexVDpmkMa5k+m2YvKhYx2mnECbKJ5rWiJBw3gwUcI+U8zRVKuKoxde3m4vAhYn/bBabw9m5odgxdTcVkwPg6fOPoVcplVGIQjPwdNiNthjL5vBU8Nw6ddeGrcB6YkS9lbJWB8bYKVc6/J2CzHSHR0By8eJhD3quEexVBwrJhfTiql5WDFx8vR7sWK6PPaA7KRao0ml3jTCPiD6FSDjr08GE/XtchZHXqnU2enx/Gql0bUcWxyPvSXsMVZQclsjmYy2XivqfoUkrJgYHruzT7EPKyZoPsHqmJUTEJGHgGWgAdRV9eBoW2QwhDNZEbu9OLKfqKxU6l0LKBeH5bFnMz1bMfEGQbuzYnpNd4xWj71b2OeKOar1puf+TqmBMcuK+S+qeq0R9fFmWHcjaRgPmShhn3FF7J00m8pqtd6K6h3yEYW9YhfkAlwCHUX4nBWGJJbn3U5bdNWYidQBuTz9noqAif074p2M67w4BC2Pt1odr4h9M2FsmMlhooR9zonYPTz2tVoD1e7FHeJYMU5anzMzM06tmLgThuoeVkw9QvEwZ5uCu2xvzJm10G9WjH9t/FaFzfHJilHgEyJyl4jc7LWBe62B2vm1ITfPMGkkUX5ibL5dUQiK2P283agCVusQsKirGrktDsHJY4+RV76hTHA8K6Yfjz3y4GmjSdFjghJ4V3hsV9iMULc4HXybqh4TkZ3AbSJyWFU/7d5AVW8BbgGYu2r36O/DDZueiYrYt0xb+cAX1rsLUJ1fsx5bmNrosZfyVlZM2KpG5XqDoquIej4rERd7dkXsPaQf5jL9WDG9e+ylfJZKPXwBknKtSTG/8TJamLbOsdfn4DzW+TmkFVU9Zv8+BXwIuG60LTIYwpksYZ/KkxE4s1Lteu7MipUGdtFsccPjU/ksjaaGRrXr1QZThbawF3JRa8y0I+GWpdLjBKW4Vkw2I2QkurBnhFb53al81nc+gJtyrcFUfuNdkHOOz6x2p96dtj+bzs8hjYjIjIjMOX8DzwPuG22rDIZwJsqKyWSEbTNFb0FZdQRl4yy/Ut4S6/Vao8srdmMJmDtijybs7qqJ2DoZzVJxBl17s2Li1nxxltNzmCpkWa8GC7uqsl5rMFXYeN62z1jn2K+DzWWE+dJYROy7gA/ZReNywHtU9WOjbZLBEM5ECTtYwn3aQ1BOL1tiv32mI2K3o/ByrRFoD6x7CnuUtMV2tokyeCvGSZPMubJboixmXavrholGpXyWci34eLWG0mjqhvMCliVm3Tl1d7BnVqpsmymQibDE36hR1W8ATxp1OwyGuEyUFQOwfbbgLSirFbIZ6RJvR5TCotNOKyafjVavvNZoIrbFMQwrpuayYpzfUTuSvOuOZSqfpdpoBrbVsWpKHcKezQjbZgqtuyQ3Z1YrbB8DG8ZgGGcmT9hnipzxEhSfSHG60LZigijXmh3CHq1eedVeBUnESpF0HgvDy4qJd4fQLn8QWdiz7XPjnJdywHss2+dsutB947d9pujZwZ5eqXbZYQaDIVn6EnYReZOIPCoi99g/L0yqYb1iReweVsxKteX9unF77H7UG02qjeYGyyHq4Gmt0Z6xKmLZKr1aMdGycDqsmFz0lZDcHnupEH4n4zzX6bGD/+dwZrUyFgOnBsM4k0TE/of2dOtrVfWfEni9vrhotshKpd6KJh38BMUR63KQgNmv1ZPH3hEJ57OZeFZMrh3tx7FiNkbs0e4Q3B57FIvK67w4bJ/1v3Py6mANBkNyTKAVY4nG6Q4b4MxKle0eFsBUBCum5SX36LG7I+E4Qgu46tNEHaztXjAjaj32vJewRzkvXsI+U2gNWDusVeusVRvGYzcYBkwSwv5qEblXRN4pIlv9NnJPu15cXEzgsN60cqg7bIAzK8ER+1pAZFquNjdsC9G962pdO4Q9eocAtHz5XIwJUe589Fgee659Z+HYK4HCXvWP2HfMFVnuuHNyPhOvDtZgMCRHqLCLyO0icp/Hz4uAPwOuAK4FjgP/x+91VPUWVT2oqgd37Oi/FoIfjmi4c9nXqw1Wqw1PQYnisXtZDnGyTTpLEURb4q6dTeMcL6oV00tH0uWxR7FiWh67d8QOcNZlx5xuTRIzwm4wDJLQPHZVvSHKC4nInwO39t2iPnGicncuuyPyF814ROyuPHY/WsJe6LBUouSHe3js0RalbmfTxDpeRz56vKyYbism0nnx8djBitL3bplq/Q3dcwkMBkOy9JsVs8f17/eRgunWrYh9xR0p+lsA03GyP1xT5/MxsmLcgpmL4827UjNzMbJpcq6OxMreiT946qQwRrqT8YrY7XN92nXn1Opg54ywGwyDpN+Zp28WkWuxSps+BPyPvlvUJ9OFHFP57IYcaudvr0G7Ui7ciil7CFh0i2Ojxx55MLNjwlAhYn1078Ha7mJcXvvNu0rpRsmKKQdE7M7dkWcHa7JiDIaB0pewq+rLk2pIkmyfLWxItTsTICiZjFDMZQKFfc1jkDCyx96xJmg+m6HejBZBdwp0PWLk3eWxR6yrvjGP3fp7LcrgaUDEvrGDrTJbzHlm0RgMhuSYuHRHsCJzd7qjYwf4ZWNMFbKDzWPP9Wap9GrFdHv6vZUUgGj5/aVct1BPF7KU8pmNHexqxWTEGAxDYCKFfc98iUfPrbf+f/TcOvOlnOfUdwgvUdvOY3dZI7kMlQhlbTuzTaIuxdcptFGXquu0Ygo5q958lHZ2FgGDcI+9mMt4FvQSEXbPlzh6rr2i0KPn1tk1Xwpti8Fg6I+JFPards/x0JnVllVw6PgSV++e993eEvaAmigeVkwxomCWa40NEW0hohVT9/Dme7FiirlsxHY2KeU3diT5rASPPXQURuvkqt1zHD6+DFhrzh4+sczVu+dC22IwGPpjIoX9cXvmaSrcf3KZZlP52ollHrfHX1BK+Szr1e5l3By8ZliW8lnqTQ3NSe8UzKiWSrXRbE1OirNfpxVTymcCUxbb7dy4QpS1b3BN9rVqw3Pg1OFxe+b5pt3BPnp+nZVKncft8e9gDQZDMkxcPXagJeKHji+xdbrAarXB1QGCMlUIt2LcS80BLbEu15vMZv37x3KtsaFDiGPFdE5sWvVYy9Vrv66IPaSuOkCl1uwa1JzKZ0Pz2IOE/erd86jC104uc3KpbD82fhG7iGSBO4FHVfXGUbfHYAhjIoX9kq3TzBSyHD6x3FoHNShSnC5kA0Vz3SMydaLbSq3RWrzZi0q9SbEjbbEXKybqYG3nfqW85c03m+q7uEWzqdai1B0rSE2HdHjlWrAV83j7nB8+vsTJpQoilj0zhrwGOASY2w3DWDCRVkwmI1y1e46vHl/i8IklROCxu2Z9ty+FeeweAuaO2IPojNh7tWKilvutdkxQco4d5LM7z3VG7GFWTFjEvm/rFDOFLIeOL3Ho+BIHts/4DmCnFRHZB3w38I5Rt8VgiMpECjtYEfqh40t89dgSl4UISpjl4OUllyJMuVdVK2LvTJPs0YqJW/8daEXhQe0st8YQNl4OoRZVyOBpJiNcvWeeQ8eXOXRiKXCcI8W8FXgd4Hvy3QXuaufX/DYzGCJxYF//RRInVtiv3jPPcrnOHV9b5OoQQZmKEJl2RrOOYAb5104kXOxKW4xmqWyM2KNbMb1G7J2Dp+HnpduX7+Tq3XN88cg5Hj6zFpiZlEZE5EbglKreFbSdu8Bdfsv0kFqXPNWr9426CYaEGK/74hh81+N3cddDZ6k2mvzoMw4EbhsWmXpZMU4UXq777+eIfqkjTTJaud/mBqEt5KLlsVteuTt7p4+IPZ/lwrp/OYJyiBUD8NLr9nN+vUZGhO950t7Q9qeMbwNuslcGKwHzIvIuVf2REbfL0APD6riWLx39Qu0TK+w750u89SVPjrRtKWyCkpcVkwu3YhzRdwtmMRcx/bDeoNjLfvakoa52BnRA7XZ2vMcoVkyIsF9z8QJ/8sNPCW13GlHV1wOvBxCRZwO/YETdMA5MrBUTh6l8lmq9ScMnW8VrkNAR3UArptZtcTgTm1SDbZVKbWOWSjEfbUKU5elv3A+sfHo/yrVuywjssYewwdMAj90wfC5cYSpnGoywA+06634R8XqtsWFZPGhHwpVIkbBbaK39wmyVSr3R0SFYnU94h7Bxv5IrLTNoH6udHh57SB77Zinopap3bIYcduOzTwZG2IEpO2PGb3m8crXBdFdWTJRIuLtIVmvQNST6rtQ3zlh1/o67X3ssIKCdrXTH6HnsjaZSrTdb9ewNvVPZ373odydp8G3HmX47rHG7EzLCTvhqQWueeezhHrtXfrgjtGGzQSsdg6etCVEBAl1vNKk3tcv6CWun85xXSYFyzZrc1EnQ6kkGg2G0GGHHtaiEnxXjOfM0PIJuCWa+O688yMKpNyy/f4PHHmE/x97ZMHgaa4JSdx47eA+8OmmQnRaVYfwxdsz4Y4SdtsfulbPdbKptb8SP2B2bJq4V41gjnh1CpMFaj5o2PUTsQasoBa2eZDCkiWF2VGmxzIywE1x73IlWu/LYc+Eee8Uz3THciql4CG3LwgnKm291CB4WTo+Dp+B9XowVk17GzQ/ebEQZU+mXvoVdRH5GRL4mIl8RkTcn0ahhEyhgHrXYAXLZDLmMBOeHe6U7tmrMRBBoDysmymCtd8QeId2xw4pxbBavaL+9LJ6JDaKQxDTxYWLsmPGmrwlKIvJfgBcBT1TViojsTKZZw6XlJXtYDkGRaSkfXBLXa0ZnrFIEnt58vNIApUiRvvcSd20rpvuYXjXqDYa0kUQHNY53QP2GWz8J/J6qVgBU9VT/TRo+03mrf1vxKN27FrBgcymfCYnYncHT/oU21n7uhT0yQkbCI/aMsGGBDqCVyuh9Xur2NhM7eTl1DNu/Hfeofdzb3w/9CvtjgWeJyOdE5FMi8q1+G7or4C0upuu2dH7KEqflcreALdm1Uuan8l3PFXPBVSG9sk0iRd4e1kivEbuI2GmLwR1QKZ9FZKNwzJes97xc7q4Xs7Ret7cxwp5GkooyN7M4jjOh30oRuR3Y7fHUG+z9twJPB74V+DsRuVw9pkeq6i3ALQAHDx4ML1M4ROZsAVvyEjD7sQUvYQ+Z5l+pNRCho4xub9UWow26+njl+eB1TzsXA3Fw3vOSV4cXcF4MhlGTVIc0jjYMRBB2Vb3B7zkR+Ungg7aQf15EmsBFQLpC8hCyGWGumGtFoW6CItNSLhuYbVK2BdMdCbc99ijph901X6JYMZ0iHVZArHMxEAfnTmbJo8Kj85jTKRoml+rV+ygcPjrqZkRmVHcZSVhlSQ2y92vFfBh4DoCIPBYoAKf7bdQomJ/Ke5aovRBgxVgLRQdH7F113COUBvCO2OPUf/eYQRqSN+8l7M6Sf37nZSqf3bAYiGFyMZbMeNHvt/KdwOUich/wt8ArvGyYcWCulPO2Yhxh94hMwzz2cq3b4ohmxXjNWI0+eNpLxO5lxeSyGWaLfuel3oroDelkXG2EfjAdkEVf30xVrQITUZ96YSrvbTmU/SPTUj7jmTHiUK77r7wUOOjqMYM0SqTvl49ejOKx+6QtWufF22M3/npyVPZXKT5SGHUzAkm7JZO0qI9zx2juo23mp/Leg4QBkWlYtkml1uzKDY+T3eLuFEpRIn2PapLW/xE8dh9LxfdOplzzvIsxDJZRT1lPa0Sc1naNCiPsNvMl74j9wrq/gBVzIXnsHasggZV+aC22Ec9SyWcFkZDSAB4Tm8CZSBVcUsAvYg8ae/Aad5gkRKQkIp8XkS/ZM6t/fdRtissgok4jounHCLvN/FTO14rxE7AoM087o2ewV1GKOQja7hDC93OnV1rtDN/PL2L36/CW1uubIYe9AjxHVZ8EXAs8X0SePuI2GToYREcTt0Mc9Z1UJ0bYbRam8ixX6l3L4wV5yaFWTMcydQ5hnrfzmp2+fviEqAa5jJDLxtvPL90R7PPik8c+6R67WqzY/+btn7FMDkiatETtaWlHVIZRAAyMsLfwm2VpWTHekWkxnwlOI6x5pxGGWzFN8lkhm9kYBYRG7B5ZOBCelmm10ydin8p1WTHNprK0CawYABHJisg9wCngNlX9nMc2rVnVtfNrw29kCIMaBBy1qI76+GnGCLuNI1KdGSDW4Kmfx26tQ+q1whA46496ROyRBNqjQ4hgqXh55VYee8hYgMfxwOrwVip16q41WlerdZrqnQI6aahqQ1WvBfYB14nINR7b3KKqB1X1YH7LtO9rJTX5JE23/dWr941EYAd5zHHOhnEwwm7Tnj7fjk6bTWU50IqxTp/fwtReC3RY+4UMZtYbnhF0KZcNXjzbJ7sl1NMPiNid9+5O63SyhybdinGjqueBO4Dnj7gpPTFosRqWuA+6I+nlPKWpo3Uwwm7j2C3ugcKwyNQZGPXzry3vuoeIvR4QsYcMugZF7F5zx1TVM9/ewetOpl0YbbIHT0Vkh4hssf+eAm4ADg/qeIP2X4ch7oMS3WHcGYw6Uk+yZr8RdhtHwNx+8oUQASuGLGJhzej08tiDs2n8inIVc2ETjbytn1I+i6r3nUW10US1e7aqg9PheZ6Xybdi9gD/KiL3Al/A8thvHXGbgN6jxGGIV9ICbLx0i+fujh5TTHbIFYN5DyumXQAsOGL3s0esGixeWTEZzq36R2flWsNzpmtYaQD/DiHjej7btQ/4L5jhfV786+dMEqp6L/DkUbdjHBk3Me61w4vTwQ4rIwZMxN5iwctyCClN217QujsSrjeaNJrqn8fewyBo74Ou/paR12IgbtrnxSXsm9BjnyRGbTmkjUk8H0bYbWYKWTLiZ8UED556CmZAJBxqqfhm0wQPnlY8ZroCrQFVL/vHecx3glKQRTX5Vkxs5nPlwOeT9FH7GbSbRDEbNmkcNHUwwm4jIna9GA/LwbekQJRI2C9LpQdLJUq6o+d+vUfsrUFlj/MyO/kzT4fOMG/XDZPbwRlhd9E5fd6xHPyLgPkX9Gp51x7WSJQVjTzTJKMMunruF6Wd3pfCTCFHRrotqrlirmsClWH4mKi9d9L0/pO8kwMj7BtY6KjwGLZKUClSJNzLIKiPFZMPnrHqV1c9Sjv9Bk8zGa87Gf9JW4bxIk3iNk7E7VD7vROLkxEDRtg3sGU6z5mVSuv/M6sVFqbyvpGpE7GvewjmetVfMMPXSvUZBO0x/91pg2c7Q4QdYOt0gTMr7QvzzGqFrTNG2A3jy6R3aEbYXVyybZoj59Zb/x85u84l26Z8t58tWuK2WukWzFV7puZcsdvGKeay1Ju6YZq+G79BUKeYl98iVX6Drs4Sd6sei4I4j816tNNh39Ypjpxr10A5cnaNS7b6T5039Efc6K7fQbxJF7lONsP77UvYReR9InKP/fOQXSxpbNm/bZqzq9VWIbAjZ9fYv81fwJzBw5VKd1lbZwq+1wCjI76+pQh8inkVcxmaCnW/2jQ+1STn7DZ4VWl0HpsLGAjdv22aR85awt5sKkfOrQeeF8P4sRnEDpJ5n8O2YXqhL2FX1R9S1WvtIkkfAD6YTLNGgyNWR86u02gqR8+tc0mAgE3ns4jAiodgrgREwmELUweVFHCe70RVffdz2uC1jF9QOx32b5vm/FqNC+s1Ti1XqNabgefFEEzSA2VJsVnEfTOQiBUjIgL8IPDeJF5vVDjC/sjZNU4ulak2moGRaSYjzBZyLHsIphMJe0bsef9l7ppNpdrwLykA3qsoOdG/134zjrB7dUD2YzMhwg7WHYwTuZuIfbAM247ZDIwiWh8VSXnszwJOqurXE3q9kXBJDwI2W8oFRuxzxe5BxnaapL9Ae1eFjJC26LFfIZehmPNeeHulUqeYy3iWMHDo5bxsZuJmMKSJSY3aR/W+onTQg7iDC51hIiK3A7s9nnqDqn7E/vulhETrInIzcDPA/v37YzZzOCxM5VmYyvPI2TUWpi1BDhX2Ys5bMMt1shnxqe7oH7E79kxgxB5zP7A8dM87i0o90F8H2L+9fSezWm2QEdi7xX9Q2TAali8V5h7uf4GnC1cUWXiwEr6hIbWECruq3hD0vIjkgO8HnhryOrcAtwAcPHgwtcuLOQOFW6bzkQRsxk/YK3VmClksl2ojjvh6lyIIzn/33c/Jm/cR9tlizjcrJsiGAWvi1pZpq8NbqzbYszAVGOEbDGkiqWh9XGwYSMaKuQE4rKpHE3itkbN/2zRHzq5x5Owae7dMkc8Gn6K5krewL5frvhObokXe8QZPWwtg++Sj+1pG5XrgwKnDfjsV9JGQTCFDcoyyvMCkWDKjfB9JfX69WHtJCPtLGPNBUzeXbJvm6Ll1vnkmmoDNFv089pqvYLYE2qsoV90/8g4aPA3az2mnnxUTRdgvsTs8I+zJMKjMmCSjygtXFMdW4Me57UnQt7Cr6itV9e1JNCYN7N82TbXR5N6j56MLu6fF0fAtkhVkqbQi74C66l4LaAft12qnT8Qe5rGDdV6+eXqVxeVKy3OfdETkEhH5VxE5JCJfEZHXjLpNo2CcRHJQbR0nGwbMQhtdPO8Ju7j/5DK1RpOXPe3S0O39LI7lSp0tPvVUpgtW5L3mIeyODz5d6P5opuz91qvdx1uzZ7967QcBg7wRI/YfPHgJK+U6GYEXXbs3dPsJoQ78vKreLSJzwF0icpuqfnVYDajsr1J8pDCswwXiCGZaB1bHpfNxM6g7NyPsHVw0W+RNNz0h8vZzxRwr1TrNppJx1ZRZKdfYt9V74HUmaIp/1ckrD5po1N0hOKLttR/YHZCfsEeI2C+7aIbf/N5rQrebJFT1OHDc/ntZRA4BFwNDE/ZeSCo7xo+0Zc0MWtB7idZHXX7ZpDb0yQ3HgQgAABFtSURBVGwph2p39L1SqXvWiYFgYXdE2yuKDuwQQmaQzhbzAYOnpqBXGCJyAGuZvM95PHeziNwpIneunrO+0EnmsvciEoO2DhzLY5RR8jCOP24WjIMR9j7xE9uVsn8a4UzBf4r/WsV/JmirmJeHFdOO9L2POVfKUW00N0yKqtQbVBtNZn2ifIOFiMxilcx4raoudT6vqreo6kFVPTizNbptMujSAsMSpWGK/DCP1ev5SzJa7zVAMMLeJ47YugtsNZrKarXhGz1nM8JUPusTsfsLdDGXIZuRwP38I3anA2oL+2rA3YHBQkTyWKL+blUdWS2kXsVi2BHnIER3FHcHwxD1QXbs5hvdJ3Ol7ujbiZ6Dsk2siU1eg6fWYzOF7ihaRJgpZH3LBGczEpgVA9adxLaZQutvgFmzdqkndg2kvwAOqeofjLo9vQ6kOiI1SN+9k3EcyHQYV/vFjYnY+8Txp93+dUswAyLh2WKWNR9LpZTPkPOZGBWUXuk30xXaxciWXSWGnb9NxO7LtwEvB57jKk/9wqg7R7mNHmalx0kQrEGyfKn0dY5GPWDqxnyj+6SdqdIWzKBa7A7TBe8p/lYpgoD9eiwNMOeK2FvHilCLfTOjqv8OpEoN+01/HHTGzLjSb6cXV9SjdOj9DMCbiL1PvBaxiFLj3D/yDhZov9o0q9Xg/WY9LKMo7TSkj34jw34j00ljEs+FEfY+8VrEIkokPFP098oDBbroN+jaCO0QutoZ4c7C0B9ps2PcTKKgxSGpDm4Q0Xq/GGHvE690x3Yk7D8oOeNjqVgzQf3TD2cKOd8OIWi/uQBh98u3N6SXpPzczSruSb3vQfnq/c6DMMLeJwV7kYplj4jdbxYohAyC9mrhBHjzLSvGw2MPK9trGDy9RHFJivtmEvi0i3oSGGFPgLmOAlvLrUg4fsQexWP3mqAUVvNlKp8lI90Ru0i7do1h/EhSXCZd3NPQgQ3LdjPCngCddViiROyWQDdoNjdmKKxU6swGRN69dggiYpXudXdAdi12vxRJQzKM01J5oxa+QZH0+0pztA5G2BOhsyTuSqXGVD7rm4tu7eNd4XGtGmbFZKk1lGpH6d4wCwdgrpTvitiNv54eeo3mkhaZNES2STJuop5EIGC+1QmwdbrAmdX2h31mtcqW6eDZnO5BV8dCUVVWqyGDp679Cjkrn7lab0aq+bJlOs9ZVzvPrlZZmE5HSVhDfwyivK9bEMct9z2NHdMws59MxJ4AuxdKnLhQbv1/4kKZPQulwH280iTXqg1UgwczvdIWHWsmLGLfs1DiuKudxy+U2RvSTkMyjJMd44UTxadRMGFj+wbZxrRbMA5G2BNg70KJU8tl6g3LHjl+ocyesEWwC91pklEE2qvCY1DhMDd7FqY4fmG99f/xC+vs2WKEPU30E9UNS3TSIvLDbscwzm9SAUBfVoyIXAu8HShhrTbzU6r6+SQaNk7s2TJFU+HkcoW9CyWOnV/nhsftDNwncMJQhIh9Q4dQjTaDdM+WEufXaqxXLV///FqNPQvBHZBhvBj2ikudojpIy2bUHUk/DHsSWr8e+5uBX1fVf7aLI70ZeHbfrRozHNvl+Pl1pvNZKvUmu0ME08mY8SqjG5R+6FR9XNmwn7OcXrDH7rTzmCtq32si9tRxYN8iDx3dMepm9ETSvnxaxLyfaH0UM4v7FXYF5u2/F4Bjfb7eWOJEvcculFvrkoZ518EzVuNF7EGrLnm18/j5ts++e95E7MPiubsPc9uJqyNt24+4OyI06rVS0yLK/TIuvrqbfoX9tcDHReQtWH79M/02FJGbgZsB9u/f3+dh04XjU5+4YEXs1mPBguk1eBrHY+9l8HSvI+wmYt8UpGkh7HEkCUGPE60nOcAeKuwicjuw2+OpNwDXAz+nqh8QkR/EWpTgBq/XUdVbgFsADh48OF65UyHMl/LMFnMcO19mKh8vYnfXZA9b3s79nNega1jEvmvBWvzAnRmz22TFpJYkLJm0RO/jxrBFPWlChV1VPYUaQET+GniN/e/fA+9IqF1jh5VKuM50IUsuI2yfDV5Bxons3V55NCvG8ebjR+zFXJaLZgutiP2i2QLFnCkn4IeIvBO4ETilqtck8Zpx7BhIzm830Xs0xtF28aLfdMdjwHfafz8H+Hqfrze27NkyxfELZY5fKLNrvkQ2E+wvZjLOMneuPHZnWbyAiUbFXJZ8VlitugZPq+H7tdq5MMWx82WOnS+bjJhw/gp4/qgbkRSV/dWJEa5BkOS5iRutJz3PoV+P/b8DbxORHFDG9tA3I3sXSnz12BJT+Wxk37qz7ksrHz2gVozffvmsRIq+9yyUeOjMKgAHts9EaudmRVU/LSIHRt2OpLNkTPS+kaQ7u1FaMA59Reyq+u+q+lRVfZKqPk1V70qqYePGnoUpTq9UePjMWuRIuLME72qlznQhSyYk2p8pdO8XtfTu3i1WxH78fJm9IQO8hsHQS3SWtFiY6H0w5yANog6mVkxiHLhoGoATS2UObJ+OtM9sybvaYhhzPe4HcOn26VansH9btHYagnFnfC3sGVxnOYj8drewbZYoflAdWq+iPohyE0bYE+KF37KH7TNFas0m1x3YFmmfHbNFjrkyVE4tl9kxFzzoCrBjrsip5Urs/QBeet1+Lt8xC8DTLovWTkMw7oyvi5+wJVLGV9xBVAdHPAYxgWmSM2gGfXeSJlEHI+yJkc9m+PbHXBRrn10LJb545Hzr/xNLlUhFuXbNl/j6ydPt/S6UeeyuuUjHLOWzfOdjx3NW46TRq7jDYGendorgOAr9sGymtFgvnRhhHyG750ucXa1SqTco5rKcWirz5P1bIu23uFKh0VSyGeHUUoVnPcaIddKIyHuxSmRcJCJHgV9T1b8YbavaDKv0gJ9IpkHwRzlO0K+oD7LipxH2EbJr3rJPTi1V2Dlf5MxqlV1zUSL2Io2mcmalwkwxx3Klzq55M9EoaVT1pYM+Rj9ROwzWmgljmIKftoHeNIs6GGEfKY4Yn1xyzwQN98qd/U4slVvZMFH2Mwyerdm12Pv0K+6QrsJhaRPhJBml9fLi+bv5lYjbmnrsI8SZzn9iqdwS9yiRd2u/C2VOXoi+n2E4vHj+7pEc98C+xdR6vpNAUue2l2g97jVlhH2E7G5F7BVOLllZLlFqt7T2W65wcrm84THDeJLkrbkR92RJssMchqiDsWJGysJUnkIus8GKieKxb58tks0IJy+UWbGtGBOxp4sXz9/N+5eeEmufJCwZB7cQpcWiGSfS0jn2evdnhH2EiAi759vrpRZymdBFsAGyGWHHbJETS2VmiznmirnIM08Nw2PU4u5gRD4agxbzYa57a9RgxOyeL3FiqW2niERbnGDXQomTS2VWCjl2mdK7qSUt4u5gRH4jw4rMh2XBOBhhHzE754vc9+gFhHb6YxR2zRV56MwqM8VcrP0M48Egxd1hs4r8MG2WXqP0fgfgjbCPmN3zJW4/dBKAay5eiL7fQonPfuMMs8UcT79i+6CaZ0iAXqJ2aIvCoAUeJl/kh+2ZD9N28cII+4jZvVCiXGvy8Nk1bnjcrsj77ZovsVSus1ypm4yYMaBXcYfhCjz4i+A4CH4aBj37FfUk0mWNsI+Y5z1+N19+9AKNpvK9T7448n4vuGY3XzuxTFOVm67dO8AWGtLCsAW+kzDRHJbwp0G8vUgiSk9qDoQR9hGzf/s0b3vJk2Pvd/mOWf7opfH3M4yOfqJ2N8Pw33shrYI7aEZtu3hhJigZDEMkqYjsubsPp1JQNhNJfwZJzlg2wm4wjDFG4EdD0uc86TIUfQm7iDxJRD4jIl8WkY+KyHxSDTMYJpVB1JIx4j4cBtGRDuJ66Ddifwfwy6r6LcCHgF/sv0kGw+QzKHE3Aj8Yxu3c9jt4ehXwafvv24CPA7/a52saDIY+cAtQGgdZx4VhCPmgKoH2K+z3ATcBHwF+ALjEb0P3gr/79+/v87AGw3AQkecDbwOywDtU9feSeu2ksmSC8BInI/beDDsiH2R551BhF5Hbgd0eT70B+HHgj0TkjcA/AL4V9t0L/h48eDDSgr8GwygRkSzwJ8BzgaPAF0TkH1T1q0kdYxji3kmngG1WoR8nayUuocKuqjeEbPI8ABF5LPDdSTTKYEgJ1wEPqOo3AETkb4EXAYkJexrYLNZNmoR80IuxiGrvwbOI7FTVUyKSAf4KuENV3xlhv0XgYZ+nLwJO99yoZDFt8SbtbblUVfueBikiLwaer6qvsv9/OfA0VX11x3YtmxG4BsuiHCVp+XxMOzaSRDsiXdv9euwvFZGftv/+IPCXUXYKapiI3KmqB/tsVyKYtnizidriVUO5KxJy24xpODdpaINpx2jb0Zewq+rbsAaWDIZJ5CgbEwL2AcdG1BaDITJm5qnB4M8XgMeIyGUiUgBegpUkYDCkmjQWAbtl1A1wYdrizaZoi6rWReTVWPMzssA7VfUro2pPDNLQBjDt6GRo7ehr8NRgMBgM6cNYMQaDwTBhGGE3GAyGCSM1wi4izxeRr4nIAyLyy0M+9iUi8q8ickhEviIir7Eff5OIPCoi99g/LxxSex6yK2beIyJ32o9tE5HbROTr9u+tQ2jHVa73fo+ILInIa4d1XkTknSJySkTucz3meR7E4o/s6+deERnqdM5RXr+uNnhex6NCRLIi8kURuXWEbdgiIu8XkcP2eXnGCNrwc/bncZ+IvFdEBr+WpaqO/AdrYOpB4HKgAHwJePwQj78HeIr99xxwP/B44E3AL4zgfDwEXNTx2JuxKmkC/DLwv0fwGZ0ALh3WeQG+A3gKcF/YeQBeCPwzVu7504HPDfncjOz6dbXD8zoedjtc7fmfwHuAW0fYhv8HvMr+uwBsGfLxLwa+CUzZ//8d8MpBHzctEXtr6raqVgFn6vZQUNXjqnq3/fcycAjrA0kTL8K6SLF/f++Qj3898KCq+s0YThxV/TRwtuNhv/PwIuCv1eKzwBYR2TOclo72+nVI03UsIvuwSoy8YxTHt9swjxUc/AWAqlZV9fwImpIDpkQkB0wzhLkQaRH2i4Ejrv+PMroL8gDwZOBz9kOvtm/t3zkM+8NGgU+IyF32dHWAXap6HKwvMLBzSG1xeAnwXtf/ozgv4H8eRnkNpeb6dfC4jofNW4HXAc0RHR+sO6hF4C9tS+gdIjIzzAao6qPAW4BHgOPABVX9xKCPmxZhjzR1e+CNEJkFPgC8VlWXgD8DrgCuxfpQ/s+QmvJtqvoU4AXAT4vIdwzpuJ7Yk3NuAv7efmhU5yWIUV5Dqbh+HTyu42Ef/0bglKreNexjd5DDsvL+TFWfDKxi2XdDww56XgRcBuwFZkTkRwZ93LQI+8inbotIHuvL8G5V/SCAqp5U1YaqNoE/x7rlHjiqesz+fQprZarrgJOOtWD/PjWMtti8ALhbVU/a7RrJebHxOw+jvIZGfv06eF3HI+DbgJtE5CEsW+o5IvKuEbTjKHBUVZ27lvdjCf0wuQH4pqouqmoNq6bWMwd90LQI+0inbouIYPlwh1T1D1yPuz3a72MIVftEZEZE5py/scoi34d1Pl5hb/YKrMVNhsVLcdkwozgvLvzOwz8AP2pnxzwd65b3+JDalIrSA37X8bBR1der6j5VPYB1Lj6pqgOPUj3acQI4IiJX2Q9dz/BLLj8CPF1Epu3P53qssY/BMswR4pDR4xdijeI/CLxhyMf+dqxb53uBe+yfFwJ/A3zZfvwfgD1DaMvlWFkVXwK+4pwLYDvwL8DX7d/bhnRupoEzwILrsaGcF6zO5DhQw4q+/pvfecCyQ/7Evn6+DBzcLNevqw2e1/Eo2uJq07MZbVbMtcCd9jn5MLB1BG34deAwVgD0N0Bx0Mc0JQUMBoNhwkiLFWMwGAyGhDDCbjAYDBOGEXaDwWCYMIywGwwGw4RhhN1gMBgmDCPsBoPBMGEYYTcYDIYJ4/8DaqynBLky6L4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 10\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csr_matrix(A)\n",
    "rhs = np.ones(n * n)\n",
    "sol = sp.sparse.linalg.spsolve(A, rhs)\n",
    "_, (ax1, ax2) = plt.subplots(1, 2)\n",
    "ax1.plot(sol)\n",
    "ax1.set_title('Not reshaped solution')\n",
    "ax2.contourf(sol.reshape((n, n), order='f'))\n",
    "ax2.set_title('Reshaped solution')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Now we start discussion of direct methods for solving sparse linear systems based on LU decomposition..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## LU decomposition of sparse matrix\n",
    "\n",
    "- Why sparse linear systems can be solved faster, what is the technique? \n",
    "\n",
    "- In the LU factorization of the matrix $A$ the factors $L$ and $U$ can be also sparse:\n",
    "\n",
    "$$A = L U$$\n",
    "\n",
    "- And solving linear systems with **sparse** triangular matrices is very easy.  \n",
    "\n",
    "<font color='red'> Note that the inverse matrix is not sparse! </font>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-2.  1.  0.  0.  0.  0.  0.]\n",
      " [ 1. -2.  1.  0.  0.  0.  0.]\n",
      " [ 0.  1. -2.  1.  0.  0.  0.]\n",
      " [ 0.  0.  1. -2.  1.  0.  0.]\n",
      " [ 0.  0.  0.  1. -2.  1.  0.]\n",
      " [ 0.  0.  0.  0.  1. -2.  1.]\n",
      " [ 0.  0.  0.  0.  0.  1. -2.]]\n",
      "[[-0.875 -0.75  -0.625 -0.5   -0.375 -0.25  -0.125]\n",
      " [-0.75  -1.5   -1.25  -1.    -0.75  -0.5   -0.25 ]\n",
      " [-0.625 -1.25  -1.875 -1.5   -1.125 -0.75  -0.375]\n",
      " [-0.5   -1.    -1.5   -2.    -1.5   -1.    -0.5  ]\n",
      " [-0.375 -0.75  -1.125 -1.5   -1.875 -1.25  -0.625]\n",
      " [-0.25  -0.5   -0.75  -1.    -1.25  -1.5   -0.75 ]\n",
      " [-0.125 -0.25  -0.375 -0.5   -0.625 -0.75  -0.875]]\n"
     ]
    }
   ],
   "source": [
    "n = 7\n",
    "ex = np.ones(n);\n",
    "a = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "b = np.array(np.linalg.inv(a.toarray()))\n",
    "print(a.toarray())\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## And the factors...\n",
    "\n",
    "$L$ and $U$ are typically sparse. In the tridiagonal case they are even bidiagonal!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.          0.          0.\n",
      "   0.        ]\n",
      " [-0.5         1.          0.          0.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.         -0.66666667  1.          0.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.         -0.75        1.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         -0.8         1.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.          0.         -0.83333333  1.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.          0.          0.         -0.85714286\n",
      "   1.        ]]\n"
     ]
    }
   ],
   "source": [
    "from scipy.sparse.linalg import splu\n",
    "T = splu(a.tocsc(), permc_spec=\"NATURAL\")\n",
    "print(T.L.toarray())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Interesting to note that ```splu``` without ```permc_spec``` will produce permutations which will not yield the bidiagonal factor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.          0.          0.          0.          0.          0.\n",
      "   0.        ]\n",
      " [-0.5         1.          0.          0.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.         -0.66666667  1.          0.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.         -0.75        1.          0.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.          0.          1.          0.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.         -0.8        -0.5         1.\n",
      "   0.        ]\n",
      " [ 0.          0.          0.          0.         -0.5        -0.71428571\n",
      "   1.        ]]\n",
      "[0 1 2 3 5 4 6]\n"
     ]
    }
   ],
   "source": [
    "from scipy.sparse.linalg import splu\n",
    "T = splu(a.tocsc())\n",
    "print(T.L.todense())\n",
    "print(T.perm_c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2D-case\n",
    "In a 2D case everything is much worse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0xa20341160>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQsAAAD/CAYAAADmIfPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJztnX9sZFeV5z/H9aODhiyBDB2l01BOdTluypPdJumNIoFQuZluJ5nxNEggBVtDtMoqqyWsmGV3h2RGWjGrBQ0RTEdoWUZhkiHMOgMsM4gYMVMOsUu9K5ZkuqFJtzvu2FVtQzrZ7owCGdCKtl0++0e9cper68d79d6rn+cjPbnqvfvuva9c/vrec+49R1QVwzCMZgx1ugOGYfQGJhaGYbjCxMIwDFeYWBiG4QoTC8MwXGFiYRiGKzouFiJyl4icE5EVEXkohPpXReS0iJwSkRPOubeJyDMisuz8fKuP+p8QkUsicqbiXM36pcQXnWd9QURuC6CtT4vIBef5TonIPRXXHnbaOiciEx7beoeILIjIiyKyKCKfCOvZGrQV+LOJyDUi8ryI/MRp60+c8zeLyHPOc31DROLO+V3O+xXn+nAAbX1VRM5XPNcB57yv74dTR0REfiwi3w38uVS1YwcQAfJAEogDPwHSAbexCvxm1blHgIec1w8Bn/NR//uA24AzzeoH7gH+DhDgTuC5ANr6NPAfa5RNO5/nLuBm53OOeGjrRuA25/W1wEtOnYE/W4O2An82p39vdl7HgOec/n4TuNc5/+fAv3Vefwz4c+f1vcA3PDxXvba+CnyoRnlf3w+njk8CTwHfdd4H9lydHlncAayoakFV14GvA0fb0O5R4Enn9ZPAB1qtSFWPA6+7rP8o8DUt8UPgOhG50Wdb9TgKfF1VL6vqeWCF0ufttq1XVfVHzutfAi8CNxHCszVoK/Bnc/r3K+dtzDkUOAR8q85zlZ/3W8D7RUR8tlUPX98PEdkL/A7wF857IcDn6rRY3AT8rOL9yzT+krSCAnMiclJEHnDO3aCqr0LpiwrsDrjNevWH9bwfd4atT1RMqQJryxmivpvSf8ZQn62qLQjh2Zyh+ingEvAMpZHJL1R1s0Z92205198Arm+1LVUtP9dnnOc6JiK7gngu4FHgD4Et5/31QT5Xp8WilpIFvf78Pap6G3A38KCIvC/g+r0QxvN+GdgHHABeBb4QZFsi8mbgb4A/UNV/alTUb3s12grl2VS1qKoHgL2URiTvalBfoG2JyG8BDwP7gX8JvA34lN+2ROR3gUuqerLydIP6PLfVabF4GXhHxfu9wCtBNqCqrzg/LwHfpvTluFge3jk/LwXZZoP6A39eVb3ofCG3gK9wZTjuuy0RiVH6451R1b91TofybLXaCvPZnPp/AeQo2QeuE5Fojfq223KuvwX3U8Fabd3lTLtUVS8Df0kwz/Ue4PdEZJXSdP4QpZFGYM/VabH4B2DEsdjGKRlang6qchH5DRG5tvwaOAKccdq4zyl2H/CdoNp0qFf/08BHHav3ncAb5SF9q1TNaT9I6fnKbd3rWL1vBkaA5z3UK8DjwIuq+mcVlwJ/tnpthfFsIvJ2EbnOef0m4Lcp2UgWgA/Vea7y834ImFfHKthiW0sVYiuUbAiVz9XSZ6iqD6vqXlUdpvR3NK+q04E+l1dra9AHJQvwS5TmjX8ccN1JSlbznwCL5fopzc2eBZadn2/z0cZfUxoib1BS6/vr1U9p6Pcl51lPAwcDaOuvnLpecL4AN1aU/2OnrXPA3R7bei+lYekLwCnnuCeMZ2vQVuDPBvxz4MdOnWeA/1zxXXmekrH0fwK7nPPXOO9XnOvJANqad57rDPA/uOIx8fX9qGg3wxVvSGDPJc6NhmEYDen0NMQwjB7BxMIwDFeYWBiG4QoTC8MwXBGaWEjIG8QMw2gvoYiFiEQouYDuprTp5yMikm5Q/oF618Kgne1ZW73VVrvb66W2whpZeN0g1tYvQ5vbs7Z6q612t9czbYUlFu3YIGYYRhsJZVGWiHwYmFDVf+28/33gDlX9dxVlHuCK0t1eeX88HufWW28NvF9lXnvtNd7+9reHVr+11btttbu9drZ18uTJX6nqta3eH21epCWabohR1ceAxwBEZIdira+v8/nPf55MJhNS9wxj8BCRc37uD2sa4nuD2Pj4OE899VQonTMMwzuhjCxUdVNEPg5kKYXOe0JVF73WMz09DcDU1FSwHTQMwzNdsZGsehpSzczMjAmGYfhERE6q6sFW7++KFZzxeLzh9enpaUZGRigUCm3qkWEY1XSFWNx6662kUqmGZVZWVjh06FCbemQYRjVdIRYA2WyW3bsbx81dW1szo6dhdIiuEYtkMsnFixeZmZlpWG56etoEwzA6QNeIRZmpqSkWFhaIxWJ1y5hgGEb76TqxAMhkMqyvrzcsY4JhGO2lK8WiTKPRBZQEY3h42LwkhtEGulos5ubmiEQiDcusra1xyy23mGAYRsh0tVhkMhk2NzebGj2LxaK5VQ0jZLpaLMpMTU3RbKWpuVUNI1x6QizK3HRT45AYZvQ0jPDoKbE4fvw4iUSiYRkTDMMIh54Si2QyyerqatNyJhiGETw9JRZlmhk8wdyqhhE0PSkWU1NTrgRjbW3NvCSGERA9KRbgTTBsSmIY/ulZsYAr+0iaYTYMw/BPT4sFlBZuubVhmGAYRuv0hFgUCgWGh4cREUSEeDxOLpfbvj41NdXUpQomGIbhB19iISKrInJaRE6JyAnn3NtE5BkRWXZ+vtVPG7lcjn379rG2trZ9bmNjgyNHjuwo99nPftZVfeYlMYzW8BWwV0RWgYOq+o8V5x4BXlfVP3USIr9VVT/VqJ6DBw/qiRMnal6LxWJsbm7WvFbZ9127djXd1l5JIpFwtWbDMPqFbgzYexR40nn9JPABrxXkcjni8TgiUlcogB1TEy9CASUvSeVUxjCMxvgVCwXmRORkRYbmG1T1VQDnZ+PAmhWUbRPj4+NsbGw0LV85NWmF8fFxEwzDcIlfsXiPqt4G3A08KCLvc3ujiDwgIidE5MRrr70GwOTkpG8B8IplPjMMd/gSC1V9xfl5Cfg2cAdwUURuBHB+Xqpz72OqelBVD5YTwy4tLblqNxaLufJ+uMW8JIbRnJbFQkR+Q0SuLb8GjgBnKOU0vc8pdh/wnWZ1/fKXvyQej7O1tVW3TDQaRVVRVdbX15mfnw9cMMxLYhj1adkbIiJJSqMJKOVMfUpVPyMi1wPfBN4J/BT4sKq+3qQuV52o1VcR8dLtppiXxOhX/HpDeiLXaZlaffXqMnWD5VY1+pFudJ2GQiKRYGxsjGg0ytjY2PZ0IZvNEo0GmwzebBiGcTVdMbIYGhrSWv1Ip9OcO3eO0dFR1tfXKRQKbG1tMTQ0xP79+1lcXGxYr98pio0wjH6iL6Yho6Oj+tJLL+04l0qlWF5e3n4fjUYpFovb7yORSMMFWxCMPWNhYYFMJuO7HsPoNH0xDbn22mvJ5/Ok02kikQjpdJpsNrujzOjoKENDpe4ODQ0xOjralr6Nj4+bl8Qw6BKxgFJ8zcXFRTY3N1lcXCSZTO64Pjs7y/79+4lEIuzfv5/Z2dmmdbrZuu4GS2RkGF0yDWm0kayd7N27lwsXLtS9bm5Vo5fpi2lIt3D8+PGG1y1EnzHImFhUUD31qcX09LRtPjMGEhOLKtwsIbfNZ8YgYmJRxfz8/LZXppFw2F4SY9Awsaii0iuzurraML+qeUmMQcLEognNjJ7FYpGJiYk29cYwOoeJRROSySSpVKphmZWVFbNhGH2PiYULstksu3c3jg5oXhKj3zGxcEEymeTixYvk8/mGomFeEqOfMbHwQFk0Gk1LbHu70a+YWLRAsxga5lY1+hETixZIJpOcO3euYZm1tTUOHTrUph4ZRvg0FQsReUJELonImYpzNVMUSokvisiKiLwgIreF2flOkkwmm2Zwt70kRj/hZmTxVeCuqnMPAc+q6gjwrPMeSvlDRpzjAeDLwXSzO8lkMtxwww0Ny5gNw+gXmoqFqh4HqqNz10tReBT4mpb4IXBdOYdIv/KDH/ygaRkTDKMfaNVmUS9F4U3AzyrKveyc61uSyaSrzWcmGEavE7SBs1bQy5rRdWqlL+xV3CY8Mi+J0cu0Khb1UhS+DLyjotxe4JVaFdRKX9irJJNJVldXXYXxMy+J0au0Khb1UhQ+DXzU8YrcCbxRnq4MAlNTU009JGBeEqM3ceM6/Wvg/wCjIvKyiNwP/ClwWESWgcPOe4DvAQVgBfgK8LFQet3FZDIZ8vl803JmwzB6jaapvFT1I3Uuvb9GWQUe9NupXqe8BmN8fLxhuenpaQBLZGT0BLaCMyQymYyrKYmNMIxewcQiRMpTEjfb281LYnQ7JhYhU96p2iw/i3lJjG7HxCJkCoUCw8PDrvKumpfE6GZMLEIgl8sRi8UQEfbt28fa2prre82GYXQrTb0hhntyuRyHDx9umt29GeYlMboRG1m0QOXUonwMDw9z5MgR30JRxkYYRrdhYtECk5OTV00t1tbW2NjYCLQd85IY3YSJRQs0i5IVJJbIyOgWTCxaYHR0tO61RrE5W6VYLJpb1eg4JhYuKRQKjIyMICKcPXu2ZplYLMbGxgaqetVRi0gk4rp9c6sancbEwiWTk5OsrKw0LFMsFutei8ViV71vNEKphRk9jU5iYuESN3aK6ilI5Wik2vj5yCOPMDs7y/XXX++pHyYYRqeQZsuQ28HBgwf1xIkTne5GQ3bt2sX6+nrDMkNDQywvL3Po0KGmC7EikQibm5tEo9GGI5J6zMzM2DoMwxMiclJVD7Z6v40smlAeHTQTCoD9+/fXdKvWoiwQrQgFmFvVaD8mFk1wY6uA0kjhS1/6UtvdquYlMdqFiUUTlpaWXJUrFos8+OCDno2WfjEvidEuTCwqyOVyxOPx7SXcN910E1tbW67vP3fuHLOzs64ifR87dmzHTz+Y0dNoC7XWBFStD3iCUvTuMxXnPg1cAE45xz0V1x6mFIPzHDDRrH5V5fbbb9duIBaLKaXUBS0dsViso+3PzMwE9EkY/QhwQl38PdY7Wk1fCHBMVQ84x/cARCQN3AuMOff8dxFxv/Kow7jd21G9ZsLr/X7br4eNMIwwaTV9YT2OAl9X1cuqep7SCOMOH/3rOmKxGHNzc4HUVbkOw01wHDeYl8QICz82i487mdKfKGdRZwDSF66vr5PJZOqOLrzg1tPiFfOSGGHQqlh8GdgHHABeBb7gnO/p9IWpVMp12VqjCy/3g3tPSyuYl8QImpbEQlUvqmpRVbcoJRMqTzV6On1hNpttOmIYGxujUChsR+5Op9NEIhHS6TTZbNZTe148La0wPT1NLpcLtQ1jcGhJLMp5Th0+CJxxXj8N3Csiu0TkZmAEeN5fF9tHMpls+t9+aWmJycnJ7fKLi4tsbm6yuLhIMplsRzc9MT4+biMMIxBaTV/4iIicFpEXgHHg3wOo6iLwTeAs8PfAg6ra2nrmDpFMJhtOJ7a2tgJbpRmPxwOppxk2wjACwY/fNaijW9ZZlMnn85pIJGquZRgaGtJ0Oh1IOwsLC77XVng5EomE5vP5QPpu9B60YZ3FwJFMJlldXd3+kCptE/v372d2djaQdjKZDOvr665/Wfl83tXq0HqYl8Twg21R7zOGh4eb7npdWFggk8m0p0NG12Bb1A3PmNHTaAUTiz7DbfYzWxpueMXEYoAxL4nhBROLPsPrKtLx8XHbS2K4wsSiz8hms549JpbIyHCDiUWfUe32rTwaUSwWmZiYaFMvjV7ExMLYZmVlxYyeRl1MLAYIN8vLzUti1MPEYoDIZrOucrGaYBi1MLEYIDKZzI5crI0wt6pRjYnFANPMzWpuVaMSE4sBJpvNNs3kbpvPjDImFgNMMpnkpZdeYvfu3Q3LWYg+A0wsBp5kMsnFixdZWFhoWM6MnoaJhQGUjJ8mGEYjTCyMbTKZTFOjpwnG4OImBuc7RGRBRF4UkUUR+YRz/m0i8oyILDs/3+qcFxH5ooisOHlFbgv7IYzgcBOh3BIZDSZuRhabwH9Q1XcBdwIPOmkKHwKeVdUR4FnnPcDdlKJ6jwAPUMoxYvQIyWSy6XQEzEsyiLhJX/iqqv7Ief1L4EVKWcaOAk86xZ4EPuC8Pgp8zYkR+kPguqrUAUaXk8lkmJmZaVrOvCSDhSebhYgMA+8GngNuUNVXoSQoQNn/1vcpDAeBqakpV4JhNozBwbVYiMibgb8B/kBV/6lR0Rrnrlpb3I3pC42dTE1NuZqSmGAMBq7EQkRilIRiRlX/1jl9sTy9cH5ecs67SmGoXZi+cNAoFAoMDw9vZ3GvPKLRKLlcznUSaBOM/seNN0SAx4EXVfXPKi49DdznvL4P+E7F+Y86XpE7gTfK0xWju5icnKwb4LdYLHL48GFisRgbGxuu6jMvSX/TNG+IiLwX+F/AaaCcyfePKNktvgm8E/gp8GFVfd0Rl/8G3AX8P+BfqWrDpCCWN6QzRKNRisXgs0smEglWV1cDr9fwh9+8IU2DG6jq/6a2HQLg/TXKK/Bgqx0ywiOXy3HkyBHXI4VWKXtJpqamQm3HaC+2gnOAmJiY8CQUbgLl1MNsGP2HicUAsb6+7qn87t27TTCMbUwsBgg3MTgreeWVV3ZE1iofXnKTWMSt/sHEYoBws+8jDCziVn9gYjFAuNlVWkm9sufPn/fc9traGvv37zfB6GFMLAaMbDZLOp2m5OG+wp49e3ZkMkulUnVHIqOjoy21vbGxYZvPehgTiwEjmUyyuLjI1tbWDjvEhQsXdmQyW15eJplM1qxjdnbWc4rEMrb5rHcxsTA80yhFoqqSz+cbiokZPXsTEwsjcMpi0ojx8XEbYfQYJhZGx7ARRm9hYmGEhhu7hrlVewcTCyM05ufnSafTDA01/ppZiL7ewMTCCI2y56VYLJLP5xsuHTcvSfdjYmG0hWQyyblz5xqWsb0k3Y2JhdE23EQON8HoXkwsjLZiiYx6FxMLo+1YIqPexMTCaDuWyKg38ZO+8NMickFETjnHPRX3POykLzwnIhNhPoDRm1gio97DTRikcvrCH4nItcBJEXnGuXZMVT9fWdhJbXgvMAbsAb4vIreoavCRYY2ephyjc3p6umG58nWL6dlZ/KQvrMdR4OuqellVzwMrwB1BdNboPyyRUe/gJ30hwMedTOlPlLOoY+kLDY9kMhny+XzTZEYmGJ3FT/rCLwP7gAPAq8AXykVr3G7pC42GJJNJlpaWmpYzL0nnaDl9oapeVNWiqm4BX+HKVMPSFxotkUwmXRs9zUvSflpOX1jOc+rwQeCM8/pp4F4R2SUiNwMjwPPBddnoZ6ampsjn8+zevbthOfOStB833pD3AL8PnBaRU865PwI+IiIHKE0xVoF/A6CqiyLyTeAsJU/Kg+YJMbyQTCa5ePEiwFWxQisxL0l78ZO+8HsN7vkM8Bkf/TIGjFZTK5pgtA9bwWl0BV5TK1ZiXpL2YGJhdJRcLkcsFvOcWrEa85KEj4mF0VEmJibY3NwMpC7zkoSLiYXRFgqFAsPDw4jI9jE8POx7RFGNeUnCw8TCaAuTk5Osra3tOFf9PijMhhEOJhZGW2gUUq9RbM5WMcEIHhMLoy00yo+6sbFRM7NZLdLptOs2TTCCxcTCaAuzs7OeRxDVG8tisVjToL/VTE9PMzIyYl6SADCxMNpCObp3tQA0isc5NzdHPB4HIB6PMzc311IG95WVFUZHR00wfGJiYbSN8s7SdDpNJBIhnU43jMeZyWS4fPkyqsrly5fJZDLMzs6STqcbLgOvxebmprlVfSL15obt5ODBg3rixIlOd8PoA/bu3cuFCxfqXp+ZmRnYpeEiclJVD7Z6v40sjL7i+PHjDac2loy5dUwsjL4imUyyvLzcsMz4+Lh5SVrAxMLoSyxEX/CYWBh9ydzcXNPs7bb5zBsmFkZfkslkKBaLTSOH2+Yz95hYGH2Nm9yqtvnMHSYWRl9SjpMhIqysrDQtbzaM5phYGH1JK3EyTDAa4ya69zUi8ryI/MTJdfonzvmbReQ5EVkWkW+ISNw5v8t5v+JcHw73EQzjCk899RQi0nKcDBOM+rgZWVwGDqnqv6CUUOguEbkT+BylXKcjwM+B+53y9wM/V9UUcMwpZxhtoVneVLd1mJfkatzkOlVV/ZXzNuYcChwCvuWcfxL4gPP6qPMe5/r7xetCfsPoMOYluRq3GckiTs6QS8AzQB74haqWJ4WV+Uy3c506198Arq9Rp6UvNGpSaZwsH/F4fHuZdq3r5SNIzEuyE1di4aQpPEApFeEdwLtqFXN+usp1aukLjXrUMk5ubGxw+PBhAA4fPhxYkN9mmA3jCp68Iar6CyAH3AlcJyLlaCaV+Uy3c506198CvB5EZ43epV7A3lp2gXrGybJABCEUiUTCdVkTjBJuvCFvF5HrnNdvAn4beBFYAD7kFLsP+I7z+mnnPc71ee2GffBGR5mYmKgZsHdiYuKqsuWAN2GyurrqqbwJhruRxY3Agoi8APwD8Iyqfhf4FPBJEVmhZJN43Cn/OHC9c/6TwEPBd9voNeotjFpZWblq1FFvZBGJRHb8bDeD7iVxk+v0BeDdNc4XKNkvqs//GvhwIL0zBoJaaQJqUTZg+jVklpd/x2IxzykTy14SryOTfsBWcBqhUCgUGBkZ2R4t1AvWG41GXQfhDcJmkUqltkP5zc3NtZSGYFC9JCYWRihMTEzsmHrU+wPf3Nx0HYS3MnivF1Kp1HZ6geXlZZLJJFDaZLaxsdE09kUtBtGGYWJhhIKbzVtlZmdnm3onYrEY2WyWQqHAjTfe6LruZkGBoTTCMMFojgXsNULBi13By3dweHjYU9pDr9/vVuwhvRIE2AL2Gl1HoVBwbQtoFmuimrDyo/phULwkJhZG4ExOTroyQrqZIvhlbGyMXC531YKwekerrK2tccstt/S1YNg0xHBFoVDg0KFDgf5nb+W7F4/HPbk7h4aGiEajLW9Z90oqlWoaXbxT2DTEaAtu10K4pRWDIuxMaeiGra2ttgkFlAy7/Wr0NLEYIHK5HPF4fHvIXblzs1nyYK8JiRsRi8WYm5tr6d7KlIb1jnQ6vR3Ze2hoqC3LxyvpWy9Jow+9Xcftt9+uRvjEYjGltAO45pFKperem06nG97b7Ggn+Xxe0+m0RiIRTafTurCwoIlEwlf/WzlmZmba+tzNAE6oj79TG1kMEM3m+o3WRrhZC9ENPProo+zbt4+zZ89SLBa3R0SdWJ7dbyMMM3AOEG6s/dXfh6AMm+36ntV6xng8zuXLlxkZGfG0WCwoEokE8/Pz2ytHO4UZOA3XNFv7UOt60IbNTlA2cGazWc/rOoKgX0L0mVj0KdUbuUSk6dqHWteDNGx2irKBs5w0eWFhoaUNZH7oh81nJhZ9SvVGrlZxu8mrWzh27NiO95FI5KqFX+UNZGXDXfU9YdHzNgw/1tGgDvOGBA8tWO9jsdhV9eTz+UA8CYlEQlOp1LaHIp/Pd+BTaZ2FhQWNRqM97SXBpzfEDJx9SitLlyORSMuxIqLRKMVi0XX5bl7pWItYLBZokOBObD4zA6cRGH6mHF7v7YRXwg9BRxPvxc1nftIXflVEzovIKec44JwXEfmik77wBRG5LeyHMPyTSCSYnZ1t+f5eWYfRTfSal8RP+kKA/6SqB5zjlHPubmDEOR4Avhx0p41w+OlPf8rY2BjRaJSxsbGm//Uql4/v27ev512snaCXvCR+0hfW4yjwNee+H1LKL+I+tFGXUb2fopNDx1q5N4Laar22tsaRI0dYWlqiWCyytLTE5ORkw3smJiY8B7wt04n1Dn4I09XaM14SN1ZQIAKcAn4FfM4591XgHPACpQTIu5zz3wXeW3Hvs8DBRvV3szckHo9fZc1Op9Md6Yvf/Rlej0gk0rA/Xuqq3JeSSqUG2htS7wjbS4JPb4gruVTVInDASTb0bRH5LeBh4P8CceAxSnlE/gsu0xeKyAOUpim8853vdNONjlBre3OnFiq1u90gc4e2c5t4GJTXZlQT5GdUzgDfrSH6Wk1feJeqvuoI1mXgL7mSQ2Q7faFDZWrDyrp6ItdprYQ2e/fu7UBP2t9uu/KJGlfoZi9Jq+kLl8p2CClJ6weAM84tTwMfdbwidwJvqOqrofS+DdRaO/DKK1dpX1sIst12L3c23NOtXhI/6QtnROQ0cBr4TeC/OuW/BxSAFeArwMcC73WH2djYaIvBM5fL7QhQ06oxsVoYYrEYzzzzDKpKPp+v6/JsNZpVWPV0I2E9W1d6SfwYPII6usXAubCwsMMQt2fPHleGqUZBY/wQpkGt2khbbcCLxWK6sLDQsH/Hjh1zZSRtVk8vs7CwoENDQz1h9MSWeweH12CwlYTxOQZpPKtmaGjI0/JswztB/f6CWhpuy70DxI1QtDueY1hsbW11uguGS7plHYaJhUdquQDDWuZcyxNjDCbd4CXpCrEohzxzG2k6KKpXRLZKWEauMKchRu+xtrbG6OhoxwSjK2wWb3rTm/TXv/71jnPpdJrFxcVQ2x0bG+Ps2bO+6/GztbsRYYtFN/zu+5mwfn+JRKKlAMR9YbOoFgqApaWl0NsNqo2wokmFaR/pZ3dmt9BvbtWuEItaCly9NqBQKHjaEVmmeiNY5RGEkc/v1u5GZLPZQBZP7dmzhz179my/95Pkx3DP3NxcaIvfOmH07IppiIhc1Ynqof3w8PCOLdBuh2J+3KG1CPLzyuVyHDlyJND+1aIbfsfGTjrhVu2Lacg111xz1bnqoX11rAS3sRPC/kNshXLk7fHx8a7sn9E7tHOE0RVikUqldsQ3SKVSO4b2uVyuA726miDmoLlcjn379vVcWDmje2mXW7UrdhPt2rWL06dPX3U+l8tx+PDhju1+jEaj220HNc+fmJjwXYdhVFPefBZmmsauGFnUY2JioqPbpCtzS6yvr5PJZHzX2etxHTpJLpcjEom4ihQmIjz66KOht9HqEQahe0n8bCwJ6qi3kYwmm2xcbp5p+QiDWpG3wjzC2uTmlloRphKJhKdIWX5yl3ilWab5XjjqbT7D50ayjguF1hCLfD6vqVTApqE9AAADJElEQVSq4QcSjUZd/fJb3bnZLKRcq7QjPFurf5RhUO+Pz4uI+Qkn6JVO/6GHKRh9Jxb5fL6purvZPl2mett5p8WiHVSLbSKR2LHd3svn55cg/pAjkYiJRQCC4VcsunadRS389tXrXLEbPptWcLOMPR6Pc/ny5dD70ugzd/v5trosv5Vl+NXreXqdRCLB/Pw8yWSyP9ZZGMHiJrBvpw2tXlY2tpLAKBaL8f3vf99rt5ifn++rZElBhujrCtepG4L4BSYSib76r1GPm2++uek6jk7vDfHyHz+ZTIbqEnTTVi/vAA7KS+J6ZCEiERH5sYh813l/s4g8JyLLIvINEYk753c571ec68O+e0kwX+5O/4F0E706xTJaw0kz8FY/dXiZhnwCeLHi/eeAY6o6AvwcuN85fz/wc1VNUUo+9Dk/HSwTxOq0bgyvHgbnz59vWqZd61fqjQh7LSNZn/yjSfq52ZVYiMhe4HeAv3DeC3AI+JZT5ElK6QCglL7wSef1t4D3SwBjuCB2iA5KKDk3W+bbFR6wlg0gkUiQzWbb0n5QhLmDtFdwO7J4FPhDoPzXdj3wC1Ut/3t6GbjJeX0T8DMA5/obTnmjTczOzu74z51IJK7aot6uP9ayDaDSBbe6ukoy6eufXNspZyTz4mpcWFjoK4Fp+iQi8rvAJVU9KSKZ8ukaRdXFtcp6t9MXukRF5EceytfiNmr3ryYictJne1DKqfKPAdTTclvVRt2NjQ3Gx8dDaSsk2tlWu9tz09bt7ehIM9zI3nuA3xORe4BrgH9GaaRxnYhEndFDZYrCcvrCl0UkCrwFeL26UlV9jFKOVETkhB//r1fa2Z611Vtttbu9drfl5/6m0xBVfVhV96rqMHAvMK+q08AC8CGn2H3Ad5zXTzvvca7Pq5neDaPn8bMo61PAJ0VkhZJN4nHn/OPA9c75TwIP+euiYRjdgCfri6rmKGVRR1ULXMmcXlnm18CHPfbjMY/l/dLO9qyt3mqr3e31TFtdsTfEMIzux/aGGIbhChMLwzBcYWJhGIYrTCwMw3CFiYVhGK4wsTAMwxUmFoZhuOL/A5CG87aww1K/AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 20\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "T = scipy.sparse.linalg.spilu(A)\n",
    "plt.spy(T.L, marker='.', color='k', markersize=8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "For correct permutation in 2D case the number of nonzeros in  $L$ factor grows as $\\mathcal{O}(N \\log N)$. But  complexity is $\\mathcal{O}(N^{3/2})$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Sparse matrices and graph ordering\n",
    "\n",
    "- The number of nonzeros in LU decomposition has a deep connection to the graph theory.\n",
    "\n",
    "- ``networkx package`` can be used to visualize graphs, given only the adjacency matrix. \n",
    "\n",
    "- It may even recover to some extend the graph structure."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAFCCAYAAADL3BUJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzsnXdYFFfbxg8qti8xebONLgKCqCAqig27EBuxgB1RsaNiw95jL1Fjj4otsWAFjbFHYxcFxS72gihKZ9lld+f+/oCFRSm7O7MFPL/req/3UnfPPCc7M8/MOffz3CYAQCgUCoVCoRgNZQwdAIVCoVAolPzQ5EyhUCgUipFBkzOFQqFQKEYGTc4UCoVCoRgZNDlTKBQKhWJk0ORMoVAoFIqRQZMzhUKhUChGBk3OFAqFQqEYGTQ5UygUCoViZNDkTKFQKBSKkUGTM4VCoVAoRgZNzhQKhUKhGBk0OVMoFAqFYmTQ5EyhUCgUipFBkzOFQqFQKEYGTc4UCoVCoRgZNDlTKBQKhWJk0ORMoVAoFIqRQZMzhUKhUChGBk3OFAqFQqEYGTQ5UygUCoViZNDkTKF8y0REEDJqVPb/UygUo8EEAAwdBIVC0T+ft28n3w0eTCooFASVKhGTvXsJ8fExdFgUCoXQN2cK5ZsDANm9ezc5PGIEqaBQEEIIMcnMJM82bjRwZBQKRQlNzhTKN0RCQgLx8/MjCxYsIC0XLCCkcmVCCCGKChXI8jt3SI8ePciHDx8MHCWFQqHJmUL5Rjh8+DBxdXUl9vb25NatW8Rh/HhC9uwhJCiIlA0LI789fUrs7OyIq6sr+fPPPwnd8aJQDAfdc6ZQSjlJSUlkzJgx5Nq1a2T79u2kadOmRX7+5s2bJDAwkFhbW5MNGzYQa2trPUVKoVCU0DdnCqW0EhFBXvn4kAnVq5Mff/yR3L59u9jETAgh7u7uJDIyknh4eJB69eqRTZs2EebIEarqplD0CH1zplBKIeJ9+0jZfv1IBbmcKCpUIGXDwrRSYt+7d49s69qVLHj+nFRkmOw96j17qKqbQtEx9M2ZQillnD9/nuwfMoRUkMsJIYSUlUoJTp3SaiwTExPSIDk5OzETQohYTJiTJ7kKlUKhFAJNzhRKKUEsFpPg4GDSt29fUnPs2FwldqaJCVlw4wZ59+6d2mMpFAqybNky0rJlS/K/Hj0IcsaSlClDJp85Qx48eKCTOVAolGzKGToACoXCnqtXr5IBAwYQd3d3cvfuXfLTTz8R4u5OyKlTpFyrVkQWE0Pc3NzIokWLSGBgIDExMSl0rNjYWDJgwABSvnx5EhkZSWxtbQnx9ibk1ClSvm1bYh8fT1q0aEGCg4PJ5MmTiampqf4mSqF8I9A9ZwqlBCOVSsns2bPJ9u3bydq1a4mvr2+hn42JiSGDBg0iP/74I9m8eTOpVq1avn9nGIasX7+ezJkzh8yaNYuMGjWKlClT8OLamzdvyLBhw0hcXBwJDQ0l9erV43ReFMq3Ttk5c+bMMXQQFApFc6KiokiHDh0IAPL333+Thg0bFvl5kUhEBg0aRD5+/EgGDBhAKlSoQBrEx5Mya9eSj3Fx5JdJk0hUVBQ5evQo6dixY5Fv1z/88APp06cP+e6770hAQABJTEwkzZo1I+XK0cU4CoUL6JszhVLCkMlkZNGiRWTt2rVkxYoVpF+/fkUm0oJ48uQJCe3Shcx+/JhUYhgiJoQc9/cnXUJDNU6w8fHxZNSoUeTevXtk69atapVrUSiUoqGCMAqlBHH//n3SuHFjcuXKFRIVFUX8/f01TsyEEOLo6EhmNGxIKuWosCsTQrp+951Wb75mZmbkwIEDZMGCBcTPz48EBweT9PR0jcehUCh50ORMoZQAVNXTw4YNI//88w+xsrLSaiyl8cXII0dIVo6YS1KmDJlw8iSJjo7WOsbu3buTe/fukeTkZOLq6krOnDmj9VgUyrcO3SCiUIwcpXra1NSU3Lhx4yshlyYkJCSQESNGkAcPHpCdZ86Q8nFxhJw6RSq0a0fqJCYSb29vMmTIEDJz5kxSsWJFjcf/6aefyI4dO8g///xDAgMDSbt27cjy5cvJjz/+qHXMFMq3CH1zplCMFIZhyNq1a0njxo1Jz549yblz51gl5iNHjhBXV1dSrVo1EhUVRdzd3bM7fa1dS0x++YUMHDiQ3Llzhzx48IDUrVuXXL16VetjtW/fnty9e5eUL1+e1K5dm0TQtp8UikZQQRiFYmxERJDUgwfJ0uhocvb//o9s376dODk5aT1cUlISCQ4OJleuXCE7duwoVrAFgBw4cICMGTOG9OrVi8yfP5/83//9n9bHv3DhAhk8eDBxd3cnGzt0ID9cv06IlxdtAUqhFAF9c6ZQjAiEhxOZnx+psnMnmf3oEbk0aRKrxHzixAni6upKqlSpQu7cuaOWktrExIT4+fmRu3fvkoSEBOLq6krOnTundQwtWrQgd+7cIe0yM4lp//6ErFtH0Ls3NdGgUIqAJmcKxUiIi4sjx4KDiWlWFiGEEFOZjJQ9e1arsdLS0siwYcPIsGHDchuUaPr2y+fzyZ9//kl+//13EhAQQIYOHUpSUlK0iqdy5cqklVxOKuf82UQsJumHD2s1FoXyLUCTM4ViYJTq6bp165K0Ro1y+1hnmpiQRbdukfj4eI3GO3/+PHF1dSUKhYLcvXuXtGnThlV8HTt2JPfu3SMmJiakdu3a5NixYxp9PyUlhQwcOJAsjIwkigoVCCGEZJmakhGHDpGtW7cSurNGoRQAKBSKwfj48SO6d++OmjVrIjIyMvsvw8OBoCBI9+/HtGnTIBAIsH37djAMU+RYGRkZCA4OhoWFBY4dO6aTeM+dOwc7Ozv06dMHCQkJxX7+9OnTsLa2xvDhw5GWlpY7N4SH486dO6hfvz7atm2L58+f6yReCqWkQpMzhWIgDh8+DDMzM4SEhCAzM7PQz0VFRcHNzQ3e3t549epVgZ+5cuUKqlevjt69e+PTp0+6ChkAkJ6ejnHjxsHMzAx79+4t8KEhLS0NI0eOhJWVFU6cOFHoWDKZDEuWLAGPx8Pq1auhUCh0GTqFUmKgam0KRc9oqp4mJLtl57Jly8jKlSvJ3LlzyXALC1LmzBkia9WKzLp5k2zbtq1Y4wuuuXr1KgkMDCSOjo5k/fr1xOLmTUJOnSIPrKyIz5YtpGnTpmT16tVq1Tg/efKEBAYGEoZhyNatW0mNGjX0MAMKxYgx9NMBhfItceLECVhZWWHUqFFIT0/X+PsPHjxAiJMTMsuUAQiB2MQECzw8EB8fr4Noi0cikWDGjBnoV6UKssqXBwhBBiG4Nm2axmMpFAqsWbMGPB4PCxcuRFZWlg4iplBKBvTNmULRA2lpaWTixInkxIkTJDQ0lJVISz5iBCm3cWPunxUjRpCy69dzEabWPPbyIk6nT+f9RVAQIWvXajXWy5cvydChQ8mnT59IaGgocXNz4yhKCqXkQNXaFIqO4VI9ff/+fTLp9GkiyfFZlpQpQyadPk3u3r3LVbgakZWVRWbNmkUW3LhB5OXLE0IIERNCjslkhMkx1dAUW1tbcvLkSTJ69Gji5eVFZsyYQaRSKYdRUyjGD03OFIqOEIvFZOzYsaRv375k7dq1ZMuWLaRKlSpajaVqfFFj0iRS4dAhQoKCSIVDh4jz5MmkdevWZM6cOSQrp0ZaH8TExBAPDw8SHR1Nljx8SMrt309IUBBJXLeOLLp3j7Ro0YI8efJEq7FNTExy24nev3+f1K1bl1y7do3jGVAoRoyh19UplNLI1atX4ejoyIl6OjY2Fk2aNEGLFi0KLTl68+YNOnXqBBcXF9y4cYPV8YpDJpNh4cKF4PP5CA0NLVCtLZfLsXr1avB4PCxevBgymUzr4zEMg3379sHMzAxjx47Vaq+eQilp0ORMoXCIRCLB1KlTIRKJsH//flZjKRQKrF27FjweD6tWrSq2zIhhGPz1118QiUQICQmBWCxmdfyCePjwIRo2bIi2bdsWWtalyvPnz9GmTRvUr18fd+7cYXXshIQE9O3bF3Z2djh37hyrsSgUY4cmZwqFI6KiouDi4oIuXbqwVk+/fPkSrVu3hoeHBx49eqTRdz98+ICePXuievXquHDhAqs4lCgUCvz222/g8/lYt26dRvXIDMNgy5Yt4PP5mDlzJiQSCatYjh49CisrKwwdOhTJycmsxqJQjBWanCkUlmRlZWHevHkQCATYuXNnsZ28ioJhGGzduhV8Ph+LFi1itRx8+PBhWFhYYOTIkUhNTdV6nGfPnqF58+Zo1qwZYmNjtR7n7du36Ny5M2rVqoXr169rPQ4AJCcnY8iQIbCystJZNzQKxZDQ5EyhsOD+/ftwd3eHt7c33rx5w2qsd+/eoUOHDnBzc2O9BKwkMTERAwcOhI2NTZGdugqCYRhs3LgRfD4fy5cvh1wuZx0PwzDYvXs3RCIRJkyYgIyMDFbjnT17FnZ2dujbt69a7UQplJICTc4UihbI5XIsX74cfD4fmzZtYv22vHv3bgiFQsyaNQtSqZTDSLM5efIkqlatioCAAHz+/LnYz7958wZeXl5wd3fHgwcPOI/n48eP6NWrFxwcHHD+/HlWY6Wnp2Ps2LEwMzPDvn37WP0WFIqxQJMzhaIJ4eFI6tcPk2rUQIsWLfDs2TNWw338+BG+vr5wdnbOM77QEampqRg1ahQsLCxw8ODBAj/DMAx27NgBgUCAX3/9ldWyujqEh4fD0tISI0aMQEpKCquxrly5AmdnZ3Tp0gWft23LNdigUEoiNDlTKGqiOHw4t0VllqkpFIcPsxpPaXwxceLEIo0vuObixYtwdHSEr69vPuFafHw8fvnlF7i4uCA6Olpv8SQlJSEwMBA2Njb4559/WI2VmZmJbV27IoMQgBAwlSvTBE0pkdAmJBSKGrx+/ZpEjBpFTHOafJjKZORNaKhWYyUlJZH+/fuTiRMnkv3795Nly5aRihUrchlukTRr1ozcvn2b2NvbE1dXV3J+/HgS+/PPJMTJidSqVYtERkbqtWXmjz/+SLZs2UK2bt1KRowYQQICAkhiYqJWY928eZOUOXuWVM75s4lYTFIPHuQuWApFT9DkTKEUAQCybds2Ur9+fSJv3ZqgcvZtX1GhApl//brGieTkyZPE1dWVVKlShdy5c4c0a9ZMV6EXSaVKlcjixYvJyaAg0nDlSlL95EkSKpGQBR4epEKFCgaJqW3btuTu3bvkhx9+ILVr1yYHNUiqEomETJw4kfTo0YM4jRpFSM7vlGVqSkYePkzWrVundTtRCsUgGPrVnUIxVuLi4tCpU6f86unw8Ny9zLS0NIwePRrm5uY4cOBAkWOlpqZi6NChsLGxwZkzZ/QQffEcPXoUWytXBnKWgEEI7rZoYRSeypcuXYKTkxO6d++O9+/fF/nZGzduwNnZGb6+vnmKbZXf6cGDB2jcuDE8PT3x+PFjPURPobCHJmcKpQD27t0LoVCImTNnFquevnjxYpGJ5Pz586hWrRoGDRpkFE0zkpOTMXDgQFSrVg13FywAchK0omJFjHdwQMuWLfH06VNDh4nMzExMnToVQqEQO3bs+EqFLZVKMXPmTAiFQuzZs6dIlbZcLseqVavA4/GwZMkSnQvdKBS20ORMoaiQkJAAPz8/ODs7a9SjOjMzE1OmTMlLJEeOQDZsGDZ37gwLCwscPXpUh1Grz+nTp2FjY4Nhw4blNSZRectUlojxeDysWLGCk9pmtty6dQt16tRB+/bt8XHLFiAoCC9Wr4abmxs6duyIuLg4tcd69uwZWrduDXd3d8TExOgwagqFHTQ5Uyg5HDlyBObm5qzU0zdv3sToqlUhLlMGIASSsmWRsmsXx5FqTnp6OkaOHAkrKyu1mpHExsaiRYsW8PDwwL179/QQYdFkZWVhb58+uSrsDEJwevRorWqaGYbB5s2bwefzMXv2bJ3UlVMobKGCMMo3T3JyMgkICCATJkwgYWFhrNTTtWvXJj1/+olUyhEfVVAoyPcGtjq8dOkSqVOnDklPTyd3794l3t7exX7HwcGBnDt3jgwYMIC0aNGCzJ8/n8hkMj1EWzCmpqakHSG5KuzKhBD3pCRiYmKi8VgmJiZk8ODBJDo6mty6dYvUr1+fREZGchkuhcIampwp3zSnTp0irq6u5Pvvv2etno6OjiYNGjQg501NCSpVIoQQklmmDJl+/jx5+vQpVyGrjUQiISEhIaRHjx5kxYoVZMeOHeTHH39U+/tlypQhw4cPJ1FRUeTy5cukQYMGJCoqSocRFwzDMGTVqlVkzNGjRFa+PCGEEJmpKRkVHk5WrFhBFAqFVuNaWVmRiIgIMm3aNNK5c2cSEhJCxGIxl6FTKNpj6Fd3CsUQpKamYtiwYbCxscHp06dZjVWg8UXOPq788GGsWLECPB6Ps/7U6hAZGfm1gpkFqp3DpkyZoremKc+fP0eLFi3QtGnTbNMNlf3xp0+folWrVmjYsCHrpXelk5eDgwNnTl4UChtocqZ8cyjV0wMHDmStnlYaX3h5eRVpfBEbG4uWLVtykkiKQhMFsza8f/8e3bt3h5OTEy5dusTp2KowDINNmzaBz+dj2bJlhT7UKBSK3M/NnTuX9f7xkSNHYGFhgREjRrBy8qJQ2EKTM+WbQSwWY+zYsZyop1VVzRs3blQrCSoUilyXJy4SyZfExMTAzc0NnTp10kjBrA0HDhyAubk5xowZg7S0NE7HfvPmDby9vVG/fn3cv39fre+8fv0aHTp0gKurK27evMnq+ImJiRg0aBAn7UQpFG2hyZnyTXDt2jU4OTmhV69e+PTpE6uxYmNj0bRpUzRv3lwr4wtlInFxceHE7EImk2HRokXg8/kIDQ3VmyvT58+f0b9/f1SrVo311gCQf+l83rx5yMrK0vj7u3btglAoxOTJkyEWi1nFc+rUKdja2qrt5EWhcAlNzpRSjUQiwbRp0yASiRAWFsZqLIVCgXXr1oHH42HlypWsOmmpJpJJkyZpnUgeP36MRo0aoU2bNnj16pXW8bDh+PHjsLa2RmBgIJKSkrQaQ2m6Ubt2bURFRbGKJz4+Hn5+fnB0dMTFixdZjaXaBa4wJy8KRRfQ5EwptURHR8PV1RU+Pj753Je04dWrV2jTpg08PDzw6NEjjiLUPpEoFAqsWrUKfD4f69atM3jLzZSUFAwfPhyWlpYI19AFav/+/RCJRJg6dSokEglnMR08eBAWFhYYNWoU66X3wpy8KBRdQZMzpdQhk8nw66+/QiAQFNj2URMYhkFoaCj4fD4WLlyos7aPBw8ehLm5OUaNGlWsEOkrBbMR8e+//8Le3h69evXCx48fi/zs58+f0bt3bzg6OuLq1as6iScxMREDBgyAra0tTp06xWossViMyZMnQygU5qnyKRQdQZMzpVTx4MEDNGjQoFj1tDoojS/q1KmTZ3yhQz5//oyAgABUrVoVJ0+e/OrfVRXM+izL0pSMjAxMmDABIpGoUMX4sWPHYGFhgeDgYGRkZOg8phMnTsDGxgYDBw5EYmIiq7Fu3rwJV1dXtG/fHq9fv+YoQgolPzQ5U0oFSvU0n89XWz1dFErjixkzZui9veM///yTm0jS/voLCArCp9BQeHt7w93dXW0Fs6G5fv06atWqBR8fH7x79w4ID4dkyBCsbtsWtra2+Pfff/UaT2pqKkaOHAkLCwscPnyY1VhSqRTz5s3LPd8Mva1AKX3Q5EwpueQ0pIjbuBHNmjXTWj2tirbGF1yTmpqKDe3b5+slvbdPH40VzIZGIpFg1qxZ6FelCqTlymX3Gy9XDhl79xospgsXLqB69ero0aMHPnz4wGqse/fuoWHDhnlOXipNUigUNtDkTCmZhIeDybE6zCAE4YGBrN9ewsPDWRtfcEnGoEH5vJYzBg0ydEhakZ6ejn8cHPLNBUFBBo1JLBYjJCQEQqEQf/31F6uVFrlcjhUrVqDv998jy9Q0e36VK9METWEF7a1NKZGkHTpETHL6IFcmhHSuWJGUKaPd6aw0vhg/fjxr4wuuOHDgABl+8CDJKleOEEJIVrlyZPjBg2T37t0EgEFj04TLly+TOnXqkIdWVrn9xsWEkKNSKWFyzEEMQaVKlcjSpUvJsWPHyKJFi4iPjw95+/atVmOVLVuWjB8/nszw8CCmSnMQsZiQU6c4jJjyrUGTM6VEAYBs376djDh0iGSZmhJCss0lZly4oJW5hNL44rvvviO3b99mZXzBBYmJiaRPnz5k+vTpZOSJE6T8wYOEBAWR8gcPktGnT7NOJPpCIpGQSZMmEV9fX7Js2TIy7t9/icnevYQEBZGk9evJ0kePiKenJ3n06JFB42zQoAG5desWcXd3J3Xr1iWbN2/W+OFHLBaT4OBgsiQqisgrVMj+O0LIvqQkgzp5UUo4Bn5zp1DU5v379+jcuXOeepqFuURaWhqGDx8Oa2tr1iU2XHHs2DFYWlpi7NixhSqYpVIpZs+eDT6fjz/++MMoy3kiIyNRs2ZNdO/evdByKoVCgTVr1oDH42HRokU6K1HThJiYGLi7u6N169ZqaxeuXLmC6tWro0+fPtldxHLOyY9btqBDhw6oU6cObt26pePIKaURmpwpJQJ11NNKcwkPD48izSUuXLgAOzs7TowvuCAlJQWBgYGoVq2a2gpmbRKJrpFKpZg1axYEAgF2796t1oPDixcv0K5dO9SrVw+3b9/WQ5RFI5PJsHTp0twucIU96EkkEkyePBkikQgHDhwo8DMMw2Dnzp0QCoV6dfKilA5ocqYYNQkJCejRowdq1Kihlnpa1Vziy/7MYrEY48aN48T4givOnDmDqlWrYtiwYRq7IKmbSPSB0nSjQ4cO2WVTGqBs9CIQCDBjxgxOu4Rpy+PHj+Hp6YnGjRvjwYMH+f7t1q1bqF27Nrp06aKW2js+Ph6+vr46d/KilC5ocqYYLUr19IQJEzTuPa3qUhS7YgXe+/piuKUlJ8YXXJCeno5Ro0bBysqKtfNRUYlE16iabmzZsoXVMntcXBy6dOmCmjVr6qxjmCao9lKfP38+JGFhuN6wIfpVqYJdu3ZpPFdlF7jRo0dz7uRFKX3Q5EwxOpKSkhAQEAB7e3v8999/Wo/DMAzOjBmD9JzyHVn58kZR3nL58mU4ODjA39+fdbcqJaqJZMGCBXqph1aabrRu3RovX77kZEyGYbBv3z6YmZlh3LhxeukeVhwvX77EpBo1kGFiAhACRaVKWp9Hyi5wXLQTpZRuaHKmGBUnT56EtbU1Ro4cyfrt4vbt29jN4+Wrr33XrRtHkWpOZmYmQkJCYGZmhkOHDunkGC9fvoS3tzfc3NxYuzsVhtJ0g8fjYc2aNTrpjpWQkIC+ffvCzs4O586d43x8dZHL5Vi6dCk2V6iQ7zySDR/OalxlF7hBgwZp7eRFKd3Q5EwxCrhUT8tkMsyfPx8CgQBng4Nzm5XIypdHwP/+x4lLkabcvHmzWAUzVzAMg+3bt0MgEGDatGmcCpGeP3+Oli1bokmTJnjy5Aln4xbGsWPHYGVlhaFDh+pdvBcbG4smTZqgZcuWiP/jj+zGIoQgs2xZjLC0xOXLl1mNr2wnamlpiSNHjnAUNaW0QJMzxeBwqZ5WGl+0a9cuz5RApaXi58+fOXMpUgelglkoFKqtYOaKuLg4dO3aFTVq1MCVK1dYjaU03eDxeFi6dKlexWfJyckYOnQorKyscOzYMZ0fT1nmxefzsXr16ryVgZzziDlyBGFhYTAzM8OYMWOQnp7O6njnz5+Hg4ODWk5elG8HmpwpBkOpnjY3N0dERASrsZQtFNU1vlC6FOlyWTEmJgZ169ZFx44dNVYwcwXDMLmJJDg4WKtE8vbtW3h7e6N+/fpFlqjpmrNnz8LOzg59+/bVmajv5cuXaN26NRo1alSsb/enT5/Qr18/VKtWDWfOnGF13IyMDEycOBEikUjvD3EU44QmZ4pBuH79OpycnNCzZ0/WN9qnT5/C09MTnp6eGtX86mpZUS6XY/HixeDz+di6datR3Gi1SSQMw2DXrl0QCASYO3euUZhupKenY+zYsTAzM0NYWBhn/20ZhsHWrVvB5/M1bory999/w9raGoMHD2a98nPjxg3Url0bnTt3xtu3b1mNRSnZ0ORM0StSqRTTp0+HSCTCvn37WI3FMAzWr18PPp+P3377TWthktKlqGfPnqyXFXWhYOYSZSIZMmRIkYnkw4cP6Nq1K2rVqmWUHa6uXLkCZ2dndO3aFXFxcazGevfuHTp06IC6desiJiZGqzFSUlIwbNgwWFlZsV4FUu0Ct3nzZqN4uKPoH5qcKXrj9u3bcHV1hY+PD96/f89qrNevX6Ndu3Zo0KABJ7W9SpcikUiklUuRPhTMXKGaSApqxrJ//36IRCJMmTLFKBqCFIZEIsH06dMhEAiwbds2jX8zhmHw119/QSgUYtasWZz4dp87dw729vbo06cPEhISWI2l7ALXpk0bPH/+nHVslJIFTc4UnaOqnt6+fTurNwGGYbBt2zbw+XwsWLCA857MymXFTp06qb2sqG8FM1eoJpKUXbuQGRiIZZ6eqF69OmsBmT6Jjo5G3bp14e3trfZqxcePH9G9e3fUrFkTkZGRnMaTnp6OcePGQSQSYc+ePazOd9UucKtWrTJoFziKfqHJmaJTClRPa4nS+MLV1VWnfZiVy4oCgaDIZUWlgpnP5+tdwcwV6enp2OLjg4yc+l1JuXKQhIUZOiyNycrKwsKFC8Hn87F27doiVy4OHz4MMzMzhISE6LTf9dWrV+Hs7AwfHx/WgkBlF7gmTZrg4cOHHEVIMWZocqboBKV6msfjYcOGDaz3zfbt2wehUIjp06dzsvyoDqrLil8KzYxFwcyWlJQUnHF2ztdgI23AAEOHpTUPHz5EkyZN4OnpicePH+f7t8TERPj7+8Pe3l5vPa4lEglmzJgBgUDAWhyo7ALH5/OxcOFCoxDoUXQHTc4UzlFVTz99+pTVWJ8+fULPnj1Ro0YNXL9+naMI1eerZcVDh/CwbVv0q1LFaBTM2nL27FlUrVoVa728wFSqBBACqakp+lWpYjQqc22Qy+VYvXo1eDwelixZAvmhQ3jZuTMG8ngICgpiXZesDbfyICpDAAAgAElEQVRv30a9evXQrl07vHjxgtVYL1++hJeXF+rWrYvo6GhuAqQYHTQ5UziDYRhs2LCBtXpaidL4Yvz48RobX3DN48ePMdHRMbe/srxiRaPo060NStMNS0tLHD9+PPsvVRq1cJlIDMmzZ88wrXbtvN+sQgWD/mZKkxAej4fff/+d1fWh1F4IBAJMnz7dqIV7FO2gyZnCnvBwpAYEYGadOpyop5OTkzFgwADY2dmxMr7gkgMHDmBzxYr5ln/lI0YYOiyNUZpu9OvXr0jTDdVEYuzq88L4999/se2774zuN1MuvTdr1qzYRifFoewC5+zsnO3kpfKQRSnZ0ORMYQVz5Ei221POkqicpaHDqVOnYG1tjREjRhiFrd7nz5/Rp08fVK9eHQ+XLMnrr1ymDILt7HRmLsE12ppuPHz4EE2bNuUkkeiLjIwMBAcHw8LCAjdmzszXE3uUjY1avuC6RnXpffHixayqDpRd4Px/+AGScuWyH0QqV6YJuoRDkzNFa96/f4+jtrb53kyyhg3Taqy0tDSMGDEC1tbWOHnyJMeRasexY8dgYWGBMWPG5FkXqvRXVi4rcm0uwTVK041u3bpp1WRFoVDg999/5ySR6JqrV6/C0dERvXv3zus8p/Kb7d69GyKRCBMnTjQKO8rnz5+jTZs2qFevHqsKhKtXr2JHlSr5rkUEBXEYKUXf0ORM0Yp9+/ZBJBJhl59fruuTpGxZDDU3x/nz5zUa67///oOdnR0CAgKMwj4vJSUFgYGBsLW1LdaukEtzCa7JysrKLQnTprHKl7x48QJt27ZF/fr1cefOHY6i5AaJRIIpU6ZAJBIhrJhSsI8fP6JXr15wcHDAhQsX9BRh4TAMgy1btoDP52PmzJka7R9LJBJMnToVIpEIlyZNyl0lyDAxwVovL707eVG4gyZnikYUqJ5W2eeKiIiAlZUVhg8fjpSUlCLHEovFGD9+PMzNzY3GMk+pYB4yZAhSU1PV+o6qucTYsWMNogb+EqXpRvv27Tk13WCTSHRFVFQUateujS5duiA+Pl7t74WHh8PS0hIjRoxQ+7fWJe/evYOPjw9q1qyJa9euFfv5qKgouLi45J93zrWYvnt3rpNXQV3gKMYPTc4UtYmIiICFhUWx6umkpCQMHjwY1tbW+Pvvvwv8zPXr11GjRg306NGDdZtDLihQwawhSnMJOzs7nD17luMI1UPVdGPLli06K4dSJpJatWqplUh0QVZWFubOnQuBQICdO3dqNdekpCQMGjQINjY2+Oeff3QQpWYwDIM9e/ZAJBJh/PjxBS69Z2VlYd68eWrN+9y5c7lOXsZwnVHUhyZnSrFoq54+c+YMqlWrBn9//+z9v/BwyEeMwC4/PwiFQuzdu1eHUauPugpmdTl27BisrKyKNZfgGqXpRqtWrfRSAsUwDPbu3VtkItEV9+/fR/369eHt7Y03b96wHu/UqVOwtbVFQEAAPn/+zEGE7Pj48SN69+4Ne3v77G2inDfiV2vWwN3dXaN5K9uJmpmZYd++fSW2fv1bgyZnSpGwVU+np6cjODgYAf/7H7JMTXNVs4nbt+sgWs1QVTAfPHiQ07GTk5OLNJfgEoVCkav8NUTZ01eJRIfI5XIsXboUfD4fmzZt4jTRpKWlYfTo0bCwsOD8fNCW8PBwDPjpp1wVdgYhODFypFbzVrYT7dKlC2snL4ruocmZUiBcqqdlMhnO1aplVEpSVQXzhw8fdHYcLl2KCuLFixdo2bIlGjdubHDTDU30BtoQGxuLJk2aoHnz5hr5dmvKpUuX4OTkBF9fX432sHXFhx49OLt2VNuJhoaG0rdoI4YmZ8pXKNXTAwYMYK2efvjwIRo2bIiZdepAkdPEI4MQnA0ONsiNgWsFszqoLivu3buXk2MyDIM//vgDfD4/u0WlkZhuKPUGNjY2Wu/df4lCocDatWvB4/GwcuVKvawMZGZmYsqUKRAKhVrvZ7NFOe++33+fu+qUYWKC31q2zCsT0xJlFzgvL68S3QWuNEOTMyUXsViMCRMmwNzcHOEsGxgoFAr89ttv4PF4WL9+ffbNLWff7OnKlXBzc8PPP/+MV69ecRR98dy9e1cnCmZ1US4r/vLLL6yO//btW/z888+oV68e7t69y2GE3PGV3kBLXr16hTZt2sDDw8MgTVBu3rwJV1dXtG/fnrWrmiYo592oUaPseedcO+K9ezFmzBiYm5tj//79rI5RGrrAlWZocqYAyPYx5ko9/ezZMzRv3hzNmjUr1PgiKysL8+fPB5/Px/r163V6Y9CXglkd2LgUMQyDXbt2QSAQYM6cOUZvuqHUG2iTSBiGwdatW8Hn87Fo0SKDNj6RSqWYN28e+Hw+Nm7cqNNzVd15K5feu3Xrhvfv37M6pmo70S+dvCiGgybnbxypVIrp06dDJBKxVk8rjS94PB5WrFih1lLr/fv30ahRIzRv3lwne6b6VjCrS3R0NOrWrau2ucSHDx/QtWtX1KpVC7du3dJ9gBxy+fJl1KhRQ+1EEhcXh44dO8LNzc2omp3cu3cPHh4eaNmyJWu3tYJQnXdMTEyxn8/MzMTUqVMhEAiwY8cOVg+dcrk8twvckiVLjLoL3LcCTc7fMLdv30adOnXQuXNn1k/fr1+/hpeXF9zd3TU2vpDL5Vi5ciV4PB6WLVvGyf6pqoKZrQOQrsjKylLLpejAgQMwMzPD5MmTjaLphzZkZmZi2rRpEAqF2L59e4GJhGEY7N69GwKBADNnztSbb7cmqPqUq/sAWhzK2mahUIhZs2ZpPO9bt26hTp06nGwTGXMXuG8Nmpy/QWQyGRYsWACBQFDojVJdGIbBjh07IBAI8Ouvv7J64n727Blat26NBg0asNpLNSYFszoU5lKkarpx+fJlA0bIHVFRUQXqDT5+/AhfX184OzsbhTFFccTGxqJly5bw8PDAvXv3tB4nISEBfn5+qFmzJiIjI7UeR7lNpNR4sLWjNLYucN8iNDl/YyjV023btmUtcHn//j18fHzg6urKmem7qgp5zpw5Gr1FMAyDzZs3G52CWR1UXYoO9O+P5x07YsBPP+U33SglqOoNjg8fjtiff0b/H3/ExIkTjdpA5EsUCgU2btwIPp+PefPmaawBOHz4MMzNzRESEsLZvO/fvw8PDw9OtolUu8Dltuql6A2anL8RlOpppQCLrSgqLCwMIpEI06ZN08mT9Zs3b9CxY0e4uLio9Ubx9u1btG/f3qgVzOrwYvVqZJiYZDt8lS9fqm3/Hi1dioyc2l1ZCZ7r69ev0aFDB7i6uuLmzZvFfj4pKQn+/v6wt7fHpUuXOI9HLpfnVkqw3SZSbSc6YcKEUvegaMzQ5PwNoFRPe3p6shayKI0vnJycdN5TmWEY/PnnnxAKhZg0aVKB/bxLmoK5KM6ePYtt332Xr+HEVXf3UrmseOLECYT+3//lm+utJk1KrBCJYRjs3LkTQqEQU6ZMKfRN+MSJE7CyskJQUJDODVKePn2KVq1asd4mAvK6wDk4OOi8CxwlG5qcSzFK9TSfz+dEvKI0vhg3blyRxhdc8+HDB/To0QOOjo7Zvb1zaj6TduzIVTCr88ZirKSnp2P06NGwtLTEzVmzcm3/FJUqYX7Dhmq7FJUEUlNTMXToUNjY2CB67ty8uVasiOkuLpwkEkMSHx8PX19fODk55XsrTk1NxbBhw2BjY4MzZ87oLR6GYbBp0yattokKIiIiApaWljrrAkfJgybnUopSPd2gQQON1dNfkpycjIEDB6JatWoG9b89fPgwAr7oM7yta9cStU/5JUrTjb59++YZLqhYcKrjUlRSOH/+PKpVq4aBAwfmGYJ8MVdt9QbGxsGDB2Fubo5NHTviZefOGCwUYtCgQQbzV9Z0m6goVF3nuOoCR/kampxLGcyRI3jYpg36VamC+fPns14mPH36NGxsbDB8+HCtjC+4Jtnf36h6dGtLZmYmJk2apLbphj7NJbhGLBZj7NixsLCwUMsEhMtEYkg+hYZCnKMfkJqaGnxPXZ1tIk1Q7QKXsmtX7kMWhRtoci5FJO3YAUnZsgAhkFeowOpCSUtLw8iRI2FlZYUTJ05wGKX2HD9+PL9Dj4kJVrVuzYnNoz65efMmatWqpZXpRnh4eIlaVrx27RqcnJzQq1cvjdp4cp1I9M21a9ew68cf8z1InnF2Zt2rngu+2iZiQXp6OhY3aZL7EILKlWmC5gianEsJYWFh2FKpUr6bwYM2bbRSZV+8eBH29vYICAgwiptJSkoKBg8ejKpVq+Ls2bO5S6EZe/ciKCgIFhYWOHTokKHDLBZV040///xTa8V8UlISAgMDjXpZUSKRYOrUqRCJRAgLC9N6HC4TiT5QnffFSZNy99SZSpWwoX17WFpa4siRI4YOE0D2NpGFhQWCgoKQmpqq8feVvfi35sxR+b/0gQN1EO23B03OJRxV9fSjpUtzbwbyihUxxtZWo65BmZmZucYXxnIDOXfuHGxtbTF48OBC3xT/++8/ODo6ws/Pzygs/gri7t27qFevHtq3b4+3b99yMubp06dha2sLf3//vP1qIyA6OhouLi7w8fHh7Pdgm0j0gXLev/zyS968VfbUgex9dwcHB/Ts2RMfP340YLTZJCYmYsCAAahatapG1rDKXvx+fn7ZS9o59x1puXLw/+EH1u1EKTQ5l2gKVE+r3Aw06Rp048YNODs7w8/PTye+w5qSkZGRq2D++++/i/28WCzG5MmTIRQKsWvXLqO5MaiabmzevJnzuNLS0jhzKWJLVlZWrkGELm7OiYmJGDhwoMaJRNco5y0QCNSyl8zIyMDEiRMhEon0ZltaHCdOnEDVqlUxYMCAIreJpFIpZsyYAaFQmL8Xv8p9h8t2ot8yNDmXQDRVTxfVNajQi82AXL58GdWrV8+vYFaTyMhIuLi4oGPHjnjz5o2OIlSPx48fo3HjxmjZsqXOTTe4dCnShvv378Pd3R1eXl46/+9+8uRJtRKJPmAz7xs3bqB27dro1KkTZ6spbEhNTS1ym+jOnTuoU6cOOnXqhLi4uCLH0qfrXGmFJucShrbq6a+6Bh06hIRevTDG1lati00fSCQSTJ48GSKRCAcOHNB6HKlUirlz54LP52PTpk16fzNRKBS5Dj+rV6/W242JS5cidZHL5Vi+fDl4PB42btyot//WxSUSXaOcN9tzTCqVYvbs2TpbWdGG//77D9WrV0ePHj2QtGMHFCNHIqxfP/D5fISGhmoUI5ftRL81aHIuIXClnn769Cmm1a6dq66UlS8Pxgj2l2/duoVatWqha9euGiuYC+Pu3bto2LAhWrVqpROLv4J48eIFWrVqhUaNGhnMG1dfy4qxsbFo2rQpmjdvjmfPnunsOEWhmki4Om+KQxfzjomJgbu7O9q0aWOw/5aqiMVibOvaNbe9amaZMvi4ZYtWY3HZTvRbgibnEoBSPd2/f3/W6umHDx9in1CYT10pHzGCo0g1JysrC3PmzGGtYC4M1Te73377TWc3BmMz3cjKysKvv/6qk2VFhUKBdevWgcfjYeXKlQZfshSLxZg0aRJEIpFOziElynnz+XydzFsmk2Hp0qXg8XhYtWqVQc8hhUKB6KZN890n0gICWI3JZTvRbwGanI0YLtXTCoUCK1euBJ/Px9/DhoHJUVdmli2L0VWrGqTZw7179zhXMBdGbGwsWrRogUaNGuH+/fucjv3u3Tu0b98edevWNbqbjnJZsUWLFoiNjWU93qtXr9CmTRs0bNgQDx8+5CBC7tCl3kA5bw8Pj3y2nrrg8ePH8PT0ROPGjVl399MGZS/+yc7OUFSsmNtEpe/337PeJlJtJzp37twS3QVO19DkbKQo1dO+vr6s1dPKi61Zs2Z5N+gcdSVz5Ehus4eQkBC9NHuQy+VYsmSJ3vfZFApFbq/xX3/9lbVJhrJRhkAgwOzZs43WdEN1WXH58uVavZExDIPQ0FDw+XwsWLDAaA0quNYbqM574cKFepu36urEggUL9HJuMQyDjRs35j9PVFTYd+/eRYMGDTjZJiotXeB0CU3ORoZUKsXMmTMhFAqxZ88e1k+pSr/Z4owv4uPj4efnh+rVq+u02cOTJ0/0pmAujNevX6N9+/aoU6cObt26pdUYHz58QLdu3UqU6YZyWbFhw4a4d++e2t+Li4tDp06dUKdOHdy5c0eHEXIHF4lEOW83NzeDzfvly5fw8vJC3bp1ERUVpbPjvHnzBl5eXnB3dy9yZUkmk2HZsmWcbBOV9C5wuoYmZyPizp07cHNz40Q9rXqxabI0dujQIZ00ezCUgrkw1LX4K4iDBw/CzMwMkyZNKnGmGwqFQqNlxb1790IoFGLGjBklbgmSTSJRznvmzJkGnzfDMNi2bRsEAgGmTZvG6TnHMAx27NgBgUCAefPmqf2G/uTJEzRv3pyTbSJ9vRiUNGhyNgJkMhkWLlyoVanCl6hebL/++qtWy3Dadg0qjJcvXxpcwVwY79+/R/fu3b+y+CuIxMRE9O3bFw4ODrh8+bKeItQNr1+/RocOHQpdVkxISICfnx9q1KiBGzduGCBC7lBNJMU9qCrn7ezsbHTzjouLQ9euXVGjRg1cuXKF9Xjx8fH45Zdf4OLiotVbuUKhwPr168Hj8TjZJtLVi0FJhSZnA/Pw4UM0bNgQbdq0YV32onqxRUdHs45N3a5BhcEwDLZs2WI0CuaiOHDgAMzNzTFmzJgC68ePHz8OS0tLjB49Gunp6QaIkHu+XFaUhIUBQUG4Nm0azMzMMGHChFKz1KhMJHw+H/Pnzy8wkRw5cgTm5uaYOHGi0a6IMAyDsLAwmJmZYezYsVqfi2FhYRCJRJg6dSokEgmrmF69eoWff/4Zbm5urJfeuX4xKMnQ5GwglOppHo+HdevWsV7m3b9/P0QiEaZNm8b6YlNF22YP7969Q4cOHYxSwVwYnz9/Rv/+/VGtWjWcPn0aCA+HZMgQrGnXLs90oxQSHx+PxU2bIiOn9j3DxAT3Fi40dFg64atEEh6OzMGDsbJVK9jb2+PixYuGDlEtEhIS0K9fP9jZ2Wl0Xn769Am9evWCo6Mjrl27xlk8yhU7oVCIqVOnsn64YftiUBqgyVmf5Cgf4//4Ay1atEDTpk1Zl7coLzYnJydOL7YvUTZ7KM5cgmEY/PXXXxAKhUatYC6K48ePYxCfj8wc+01J2bLI2LPH0GHplFc+PvlqWqVDhxo6JJ2hTCT+P/yQa7EqKVcOmfv2GTo0jTl27BisrKwwZMgQJCcnF/nZo0ePwsLCAmPHjkVGRoZO4lFuE9WoUYP11s9XLwZfmIiUdmhy1hfh4bm1xRmE4NCAAayXeZUXWz7jCx2ibPZQmLnEx48f0b17d9SsWbPEKJgLIiMjA+dq1cqXrBAUZOiwdEJaWhqGDx+OQXx+tgd4zsPIYKEQp06dMnR4OiMtLQ2nnJxKxW+cnJyMoUOHwsrKCkePHi3w3wcNGgRbW1v8+++/eolp//79MDMzK3SbSBMuXLiAYebmuQ9S34pnNE3OeiItICDfjSAzMFDrsZQXW7Vq1XD+/HnuglQT1WYPCVu3AkFBuDJlSolVMKty5coVVK9eHcubNweT44+dYWKCZZ6eRmHxxyUXLlyAnZ0dBg4cmP3WpfJmolxWHDhwYKlbVlTOe3WbNnm/MSH4o1OnEq0nOHfuHOzs7NCnT59sG8egIETPnQsbGxsMGzZM7yKrT58+wd/fP2+biAXxfn757p9MCX2Q0gSanHWMcgmtX5UqkJqa5i6hBfzvf1qZO5w5c8ZgF5sqUqkUu3v1yu29m2FiggeLFxssHrYUaLqRk6wkYWG5Fn9sa8+NAbFYjHHjxsHc3BwRERGFfk51WfHw4cN6jFA3KOdtYWGR94aZ8xun7NqVm0jOnDlj2EBZkJ6eji0+Pvmuy5uzZxs0pr///hvW1tYIDAzUuP2wspKl7/ffQ1a+PEAIxCYmmOvubnDXOV1Dk7MO+apUQeXN5OLFi3B0dISvr69ahvTp6ekICgqCtbU1K+MLLnnZuXO+p9mkfv0MHZJW3Lp1C7Vr10aXLl2KNE+4fv06atWqBR8fH6Ow+NOGa9euwcnJCT179sSnT5/U+s6FCxf0bi7BNcp59+rVq8h5KxPJ4MGDi93DNVbedeuWvyf2gAGGDgkpKSkYPnw4LC0tEa7mkrSykqV169Z4+fJl7v0z68ABzJkzB3w+H3/88UeJf1guDJqcdYQ6pQpisRiTJ0+GUCgs0qT90qVLcHBw4MT4ggtSU1MxdOhQBAoEufuUWTm9d3VpLsE1qqYbBe2hF4REIsGsWbOMyuJPHSQSCaZNmwahUIh9WgifxGIxQkJCIBQK8ddff5XqeSsTiZWVVZErC8ZGZmYmJk6ciP4//ghZznUpNTVFvypVsHXrVqP4zf7991/Y29ujV69ehW4TqVayrF27ttBKlpiYGM7aiRojNDlzjGqpwtWrV9X6zs2bN+Hq6or27dvj9evXuX+vvNi4ML7gin///Re2trYIDAxESkpKvtUALrsG6Rql6cbPP/+s1VvwnTt3UL9+fbRt2xbPnz/XQYTccfv2bbi6usLHxwfv379nNdaNGzfg4uKCTp06Gf3qAdt5KxNJ7969jV5vEBkZCWdnZ3Tv3j07VpXrMjo6GnXr1kW7du0M1jJXlYyMDEyYMKHAbaLnz59rVMnCZTtRY4MmZw5hU6oglUoxb9488Pl8HB8+HO99fTHSyooT4wsuyMjIQHBwMCwtLXHs2LFCP8d11yCukcvlWLp0KSdLYjKZDEuWLDGalqRfIpPJMH/+fAgEAmzfvp2zNyepVGrUy4qq896xYwer+IpKJMaAai/+3bt3FxpfVlYWFi5cCB6PhzVr1hjFuXrt2jXUrFkTPj4+SNi6FXdbtEDf77/XyvP5yZMn8PT0VKsLXEmBJmcO4LJU4dnKlblijixTUzBG8MZ85coVODo6ok+fPvj8+bNa3+GyaxBXPHnyBE2aNOHcdOPx48do1qwZmjRpYjQ2ig8ePECDBg3Qrl27fKsxXKI0l2jdujWePXumk2Noiq7mrZpI3r17x9m4bIiJiYGbmxs6duyodkwPHz5EkyZN0KxZM51bX6qDRCLBH507597zFBUral0mperkVVgXuJIETc4sOX36NGfq6ZiYGOzh8fKJOaKaNjXYUo1EIsGUKVPyK5g1gOuuQdqiUCiwZs0anb7hqh5j0aJFBrNUlMvlWLFiBXg8HjZs2KDzNz3VZcWVK1ca7FyVy+VYvnw5+Hw+Nm7cqJN5SyQSzJw5EwKBwKB7uKq9+LWJQy6XY/Xq1eDxeFi8eLHBzlXl/WFzjme08n/J/v6sxn358iW8vb2N6sVAG2hy1pL09HSMHDkSVlZW+Oeff1iNpXqxnR49OrdZiaJiRUytVQseHh4aWfxxQVRUlFoKZnXgsmuQpqiabujjTeHFixdo164d6tWrh9u3b+v8eKo8ffoUnp6e8PT01LtARrms2LhxY70vKz59+hTNmjWDp6enXt7gb9++jfr16xtkD/fRo0fw8PDIUzCz4Pnz52jTpg3q16+vd0tM1UqWpytXZjcWURGWsn2IZhgG27dv14mTl76gyVkLLl68CHt7e07U08qLLZ/xhYqYQ6FQYMOGDeDz+Zg3b57O7euysrIwd+5cjRTM6sJl16DiUDXd0PebLMMwCA0NhUAgwIwZMzjtdV7Y8ZSmDr/99pvB9hPVMZfgEuW8lWIgfc5bJpNh8eLF4PF4+P3333V+bIVCgVWrVhWrYNYU1etk5syZOj9XgTwfgHyVLCr3vEePHqFp06acbBO9f/8e3bp148zJS5/Q5Kwu4eGQDR+O0C5dYGZmxropg/Ji4/P5ahlfvH79Gu3bt4erq6vOWmPeu3cP9evXh7e3t84K/LnsGlQYqqYbMTExOjmGOsTFxaFLly6oWbOm2sp9TXn9+jXatWuHBg0aGI0QRh96A2OZtzKRNG3aVGcrM0oFc5MmTfDkyROdHOPt27fw8fFBrVq1cP36dZ0c4/Pnz+jdu7dalSxcbxPt378f5ubmCA4Ozu4CVwL6dNPkrA7h4ZDn7Itkli2L5J07WQ2nabmAEoZhsHPnTgiFQkyePJmzftpcKpjVhU3XoMJQNd2YNWuWzlcZ1I1p3759MDMzw7hx4zgzHGAYBtu2bQOfz8eCBQsMtm9YGLpaVjTGeetKb8AwDDZt2gQ+n6+Vglmb4+3ZswcikQgTJkzg1Bzj2LFjsLCwQHBwsEbjvnjxAm3btkW9evVY2+AqXwwGC4W5/RmMuU83Tc7FIJVKcdXdnZMG+VxdbPHx8fD19YWjoyNri7vY2Fg0adIELVq00Hu9rmrXILbNHlRNNyIjIzmKkDsSEhLQt29f2NnZ4dy5c6zGev/+PTp37gxXV1e972trCpfLisp516lTxyjnrZpI2Mb35s0beHt7o379+nrXm3z8+BG9e/eGvb096979XFSyMAyDrVu3gs/nY/r06ayX3rm6n+sampyLQFmqMNfdPVvin9Or9vr06RqPpXqxcdWg4+DBgzA3N8eoUaM03sNVKBRYu3YteDweVq1aZdC6R7bNHg4dOgQzMzOEhIQYvfDj6NGjsLKywrBhw7KbuGjIvn37IBQKMX36dKNYGVCXr5YVNWTv3r0QCoWYMWOGUc9bmUgEAoFWe7jK1TGBQIC5c+catBwoPDwclpaWGD58uFbnqtIHYOjQoZz4ALx79w6//PKL1ttEiYmJ6NevH4aameW+OWt7P9cHNDkXgKp6OjQ0NHuZN2ePImb+fI0SierFNm/ePM4vts+fP6N///6oWrWq2hZ/r169Qps2beDh4WEUtY6Ads0elBebg4MDLl26pIcouSE5ORlDhgyBtbU1/v77b7W+8+nTJ/Ts2RM1atTQ2Z6grlHVG6hrLpGQkIAePXqUuHmrJhJ1fdbj4+PRpUsX1K5d22hKgMwRPL4AACAASURBVJKSkjB48GBYW1vj+PHjan1H6QPARSXLlzAMg71790IkEmm0TfTPP//AysoKo0aNyrfnrOn9XJ/Q5PwF6pQqZGRkYPz48cUmkq+ML3TI8ePHYWNjU6TFn+rykCFrcYtC3WYPX11sJZCzZ8+iWrVq6NevX5FmDBEREbCwsMD48eP14tuta9Q1lwgPD4e5uXmJnbcykZiZmWH8+PFFJhKlgnnKlCl6UUxryunTp2Frawt/f/8iGxFdunQJ9vb28Pf316nVaEJCAvr06VPsNlFqaiqGDBkCGxubQh8I1b2f6xuanHPQplShqERSYLmAjklNTcXIkSNhaWn5VS9upYLZzc1N7zWNmlJUswel6UZRF1tJIj09HWPHjoW5uTn279+f79+Sk5MxYMAA2NnZ4b///jNQhLqhKHOJ0jZvZSKxt7f/at9VEwWzoUlLS8OYMWNgbm7+VVOizMxMhISEwMzMDIcOHdJbTBEREbC0tCxwm+grH4BiMLYucDQ5h4cj2d8fU2vV0qpUQTWRnBo1CuLAQCzz9DToxXb+/Hk4ODigZ8+eSNqxA4+9vOD/ww9Go2BWF9VmD/F//IE3XbpgsFCo9sVWkrhy5QqcnZ3RrVs3JG7fjlc+PhjE52PEiBE6rwk3JEq9QZ8+fZCya1epnrdqIsnYswcvOnZEwE8/aaxgNjSXLl2Ck5MTunfvjsTt2xHv54cga+s80w09o7r0HjlzJrKGDcMfnTrBwsKiSB+Agvjyfs4YsNzq207O4eG53bikpqaQs6hdfr5qFTJMTABCIClbFpKwMA4D1ZyMjAxsaN8+t2etnEXPWkMik8mwt0+f3HnIKlQokfNQh8zMTGzr2vWbmKsqGRkZ2Nix4zcx76SkJKxs1Srvuiyhc/3qXC1f3uA+ALfnzcuNJ7NsWaT++afWYz1ftQrinPu5ocqtypBvmVOniIlYTAghpLxMRp5v3EgAaDxMamoqeb5pE6mc890KCgV5tGYNYRiG03A14dSpU8T0/HlSOefPZSUSkn74sMHi0ZbIyEiS9fffufMoJ5WSpLAwg8akKyIjI0mZs2fzzTXt0CGDxqQPbt68SSpdvJhv3iXxXFWHW7dukR8jI/OuS6mUZEZEGDQmbYiNjSUV//sv7zfLyjLobyaVSknGkSO58VRUKEj8rl1a3c8ZhiFphw6RSsrvisWEnDrFXbBq8m0nZy8vQipn/5yKihXJqvv3Sffu3cn79+/VHuLs2bPE1dWVPLK2JqhUiRBCCFOxItnz6RNp3rw5efz4sU5CL4ykpCTSv39/MnHiRNJ41qzc+WWZmpLhBw+SzZs3a3XC6hupVEqmTp1KunbtSuyGDcudh8zUlIyOiCCLFy8mcrncwFFyQ2ZmJpk4cSLp2bMncRo1Kt9vNuLQIbJ+/XqDPujpiszMTDJhwgTSq1evAucdGhpaIs5VdcjIyCCjR48mAQEBxGXChLy5litHhh88SA4cOGDgCNVDLpeTxYsXk9atWxNe794EX/xmGzZs0Pu5GhUVRdzd3cl5U9Pce7CiQgWy8t49je/nr169Iu3atSO74uMJU7Fi9l9WrpydK/SN3t/VjQ2VNm6ZmZmYPn06BAIBtm3bVqRqLz09HaNGjcpfLvBFT+w1a9boVRl94sQJWFlZISgoKE/BrBJTTEwM3N3djcriryBUTTfi4+Oz/1JlHlw2ezA0169fR40aNdCjR488326VuT548ACNGzdG8+bNdda60RAUN+/bt2+jXr16BjGX4JpLly7BwcEB/fr1y1Mwq8z18uXLqFGjBrp164b3798bNtgiePz4MRo1apS/kkVlHvfv34eHhwdatGihUedDbSnQB+CL+/m0adPUup8XWMli4BafNDkXQHR0NOrWrQtvb+8Cy6kuX74MBwcHtcoF9OFSpImCWSaTYenSpbnNRwxl8VcQqhfbzp071bqY9GUuwTVSqRTTp0+HUCjE3r17i/ysXC7PrSRYunSpUZbAqYsm81aaS/D5fL2YS3CNJgpmZSIRCoXYvn270ZTzAPkrWdasWVPk7yCXy/Hbb7+Bx+Nh+fLlOru/aOIDEBUVBTc3N3h7e+eZC6lgrJUsNDkXQlZWFhYuXJhXWnX4MGTDh2Nb164alwvo0qXo/PnzqFatGgYNGlRkzeiXPH782GAWfwWhremGNs0eDM3t27dRp04ddOrUSaM3pWfPnqF169Zwd3c3qKGHtty+fRuurq4az1sf5hJcExkZiZo1a2qsYFYmkp9//rnARKJvnj9/jpYtW2pcyfL06VO0atUKDRs2xN27dzmLR1sfgKysLCxYsCDPaOjwYTBBQTg/frxR9eJXhSbnYnjw4AFCnJwgLlOGtfHFu3fvcl2K2CaSjIwMBAcHw8LCAkePHtVqDIVCgXXr1oHH42HBggUGaRWovNh4PB42bdqk1RuDJs0eDIlMJsOCBQvUWmYrDIZhsHnzZvD5fMyePdvobigFIZPJMH/+fPD5fK3nrStzCa6RSqWYNWsWBAIBdu/erdVcv0okBlgxUPUBWLp0qVZvwKpjzJ07l/W5yoUPwJf3c7GJCZ4sX84qLl1Bk7MaKEaOzNcoXTFypNZjceFSdPXqVTg6OqJ3795FdpZSl5cvX8LLy0unFn8FoXqxcbEHXlSzB0Pz8OFDNGzYEG3btsXr169Zj/fmzRt06tQJLi4uuHHjBgcR6oaHDx+iQYMGnM3bmPUGMTExqFu3Ljp06MBJEwtD6Q3evn3LqenGmzdv0LFjR7i4uGhlSqP6YMaFD0Dszz9T44tSQ3h4dq0bIcgsUwbjHRxYLytq41IkkUgwdepUiEQihHFcR6204uPa4q8gdG26UVTXIH2jUChy9+DWr1/P6V6i0iJTJBIhJCTEqFpc6nrexqQ3kMlkWLRoEfh8PrZs2cLpXPWpN2AYBrt27dKJ6QbDMPjzzz8hFAoxadIktc/Vly9fonXr1pz4ACQmJsLf3z+f8QW1jCwN5Cj3mCNHOF1WVNelKCoqCi4uLvjll1/yFMw6IC4uDl27duXE4q8g9GW6kZSUpLG5BNc8e/YMzZs3R7NmzfD06VOdHefDhw/o2bMnqlevjgsXLujsOOry9OlTeHp66nzexqA3KFDBrAN0rTf48OEDunbtilq1auHWrVucj68kPj4efn5+xdrdcu0DoOzFn1vJYmAltjrQ5Kwlb9++RefOnVG7dm3Wy4pFuRRpomDmCoZhEBYWBjMzM4wdO5YTYwlDmW6oay7BJQzDYMOGDeDz+VixYoXeFPFHjhyBhYUFRo4cyYlFn6Yo583j8fQ2b0PpDTRRMHOFrvQGhjDdOHToECwsLAq0u1VVT7N9ECnJvfhpcmYBwzDYvXs3RCIRJk6cyHpZ8ezZs7Czs0O/fv2QsmsXPvXujXH29vDy8tJIwcwVnz59Qr9+/WBnZ4ezZ89qPY6hSxXS09MRHBwMMzMzhIWF6fQB5/Xr1/Dy8kKDBg0MooJPTEzEoEGDYGNjgxMnTujtuK9fv0a7du0MNm996g20VTBzBVd6A0ObbiQmJmLAgAGwtbVF1Jw5YIKCcGHCBM7U00rjC00rWYwFmpw54MOHD+jVqxccHBxYLyump6fjj06dcnvEZhlBz9pjx47BysoKQ4YM0egkVz68CIVCzJw50+DKYmWzh65duyIuLo7TsRmGwfbt28Hn8zF//nyDq4lPnjyJqlWrIiAgoEiLP7YY27wjIiLU2ibSBi4UzFzGwkZvcOzYMVhYWBiF6cat2bNzfQm4UE+rVrJoanxhTNDkzCHh4eGwtLRktawYGxuL/WZm+dSE6QMHchyp5iQnJ2Po0KGwsrJSq3Tr48eP8PX1hbOzs1YKTV2hSdcgdXn//j18fHzg6uqK6OhoDqLkhrS0NIwePRoWFhY6sfEz1nkXtU2kLUoFc7169ThRMHOFUm+g7otBSkoKAgMDYWtrq7YQVdc8bd8+3/3uafv2Wo/FdSWLIaHJmWOSkpJylxVz23qqgaqCOTwwMM8tq1w59KtShbNEwpZz587Bzs4Offv2zWu7+AWHDx+GmZkZJk6cqFPVNxuK6xqkLvv27YNQKMS0adMMrhoujIsXL8LR0RF+fn6ciQnDwsIgEomMet5c6A10qWDmEnX0BmfPnkXVqlUxZMgQg2gSvqQg9bS8QgUMMzdHjx498OHDB7XH0mUli6GgyVlHnDp1Cra2tmotKxaoYFZRE6omEl0qQtUlPT0d48aNg5mZGfbt25e97B4UhLS//oK/vz/s7e2LVGIaC2yaPXz69Ak9e/aEk5NTiehMJhaLMXnyZAiFwrw+xFpQ0uadnp6OsWPHwtzcHPv379fou/pSMHOFUm9QtWrVfHqDAn0ADMxXPgAq9zuxWIyQkBAIhUL89ddfxZ6r+qpk0Tc0OesQ5bKiubk5Dh48+NW/a6JgViaS3HaiRtBn+OrVqxhpZYXMsmUBQpBhYoKNHTpwou7WJ8pmD56ennj8+HGxn4+IiICFhQXGjRtnVLXF6nDz5k24urqiQ4cOGjcFiYiIgLm5eYmc95UrV+Ds7Ixu3bqppTc4cOCA3hXMXKF8MVjZqhVedu6Moebm+U03DIgm6ukbN27AxcUFnTp1wtu3b7/6d0NUsugTmpz1gHJZ0dfXN/vJLjwc6QMHYq67u8YKZk0Tia4RBwbm2y9ijLTbTnHI5XKsXLmyyGYPycnJGDBgAKpVq2YU9cTaIpVKMW/ePPD5fGzcuLHYB73SMm91XOc+f/6MPn36oHr16jqp89cXKbt2QZwjspKWK2cU9bza+ABIpVLMmTMnr5d2zirdqzVrUL9+fYNVsugDmpz1hFgsxpQpU+D/ww/ZFwshkJqaIuvAAY3HMhaXovPnz2OwUAhJzpuz2MQEs9zcjGLpXVu+avbw/+3dd1xT5/cH8IMLy+vlt7+aBJBpnTgABQSLW9x11Y17i6Kte9VR6xYrjjqrVG2rFretWne1lrpFELFSHFBxgooyQ/L5/QGJsQokuTfJBc/7v9bkyXPJOPc+9zzn5C23XZ0zB87OzggKCnprX2ZRdf36dfj6+qJp06b5Fgs5duxYsTvut/IN8t7jCzNnwsHBAZ9//rnFM5iFuHTpEn786KM3TpoPVa5ssSVfMfoAREdHY3yVKsjIq4mdRoTDI0cWu6tlXRyczexRjx6i1XW1VJei9PR0jB079vWXLe/HTbl795udvCSw9G4MtVqNDRs2oE+5ctoTqTQrK1yePdvSUxNdTk4Ovvnmm9eFQ/buBYKDkb5jB0aNGgUnJyez7pc2F81toj7lyiG7dGnte3xt7lxLT81o2dnZmD17NmxtbXF6wgRtUqnaxgabu3QRnG9gDE32dK9evQRnTyf37l0kamKLhYOzue3fr/3SaM7+hAQxc3cpOnfuHKpXr17gly02Nhb+/v5o2LChJJbejVVUCuSLIS4uDtNq1dJ260mzskJos2Z49uyZpadmUjHNmxeL9zg6Ovrtphv/KVEpJN/AUGJmT2t2suieSEm5JrZYODhbQt6X5t6qVYUuK+pLUzVIjHKi72Loly0nJwcrVqyATCbD4sWLLV6cwhDp6emYMGECBnz0EZR5WzzSiBDWuXORS4QyRPaIEW8EqpyRIy09JZPRZDAPLF/+jfd494ABReqzmpOTg0WLFunddMPQfANjiJk9/dZOliJQE1ssHJwt7K1lRQFVhzRVg2xtbUUpJ6qh+bJ17NjR4C/b7du3ERAQAG9vb4uU7jTUhQsX4Obmhu7du+fu4877MXi2ZYtBxR6KmgsXLmCko6M2fyCjRAl8XrFikdhCZKg///wTVapUeZ3BnPceP1i/3iK3iYylabrRrFkzg/M89Mk3MFR2dja+/vprUbKndXeyLFiwoEidMImFg7NExMXFoWnTpvDz8xNcgcjQqkH50XzZ5HI5tmzZYvSXTa1WY+PGjVAoFJg1a5Ykt6ZkZWVhxowZsLW1xY4dO/J93N69ey3aXEJsbx23Tve1rVu3wtbWFlOnTpVsMRlDZGRkYNKkSbC3t3/n1kbA/LeJjKFSqbSrUkKaboh5YRATEwMfHx9RsqeTkpLw6aefWqwWv1RwcJYQlUqFdevWQS6X4+uvvxZcjUhIIBHzy6ah6eRVq1YtnD9/XpQxxXDt2jV4enqiffv2eu2BTUlJwaBBgwyuAic1+hz3w4cP0a1bN1SvXh1nz5418wzFc+nSJdSsWRNdunTRq/KUWM0lxHbnzh00bdoUn3zyiWhNN+Li4tCkSROjLgxycnIQEhICmUyGdevWCb5a3r59u2Rq8VsaB2cJSkhIQLt27eDh4YFLly4JGsvQLkU5OTlYunSpKF+2d9Ht5DVhwgSLbllRKpVYsGAB5HI5wsLCDD5WczWXEJtSqdRWRtP3uHfv3o0KFSpgzJgxRWpLlSaDWaFQ6FVtSpfQ5hJi0r2iX7x4sehNN1QqlbbN6dy5c/W6MIiLi0ODBg3QuHFjxMfHC3r9J0+eoHv37qhRo4akToYsiYOzRGlq+tra2mLKlCmCfxj0CSRiftkK8/jxY20nr99//92kr/UusbGx8PX1RYsWLQTV1n758iVGjx6dbxU4qdEcd0BAgMHHnZycjAEDBqBixYo4evSoiWYoHk0Gc9u2bV9nMBtBrNtExvr333/Rtm1bszTdSEhIQNu2beHp6ZlvvoFuH4DQ0FDBSWV79+5FhQoVJF2L3xI4OEvcw4cP0b17d1SrVk1wver8AolKpcLq1atF+7IZQtPJKygoSPQWf++iUqmwbNkyyGQyg+tpF+StKnASI+ZxHz58GC4uLhg8eLAkt1oZmsGsL3PnG+g23fjqq6/M1nRDrVbnm2/wzj4ARnr27FmRqsVvbhyciwjNsuLo0aMFLytqAskif3886tEDMz094evri9jYWJFma5hnz55hyJAhcHZ2xqFDh0z2OvHx8WjcuDEaNmyIuLg40cfXbS4hpVq/pjju1NRUjBo1Co6Ojthn4X7junQzmO/cuSP6+IbeJjKWFJpu6OYbxCxciJjmzdGnXDlRsqc1jS9GjRpVpG6TmBMH5yJEzGXFzPBwbSm8zJIloZTAkuyxY8dQsWJF9OvXT9R7uGq1GuvWrYNMJsPSpUtFv1/3XxcvXoS7uzvatm1r8mIPBTHHcf/++++oUqUKevXqhcePH4s+vr50M5hXrlxp8tUfU+Yb7N69G/b29pgyZYokdjYcHT0aaZq972XLCtpjnJqaihEjRsDFxQXHjh0TcZbFDwfnIkjosmJSUhJ+cXV9o+DEAVdXiwYSDd1OXruMqDv+X4mJiWjVqhV8fHwQExMjwgz1k5WVhTlz5pi02ENBEhISzHbcaWlpmDhxIuzs7LBt2zazrxiYIoNZH2LnG6SkpKBPnz6oWrUq/vzzTxFmKIwme/q7smXf+K1I6NTJqPE0jS8GDRqkd+OL9xkH5yLK2GXFHTt2wNbWFj/26PFG7d3tgYEWCyTvcvbsWVSvXh1du3bFgwcPDH6+Wq3Gli1boFAoMHfuXIsVMYiOjha92ENB/nvc5rpPCeQWMqldu3a+Lf7EZuoMZn2JkW9w8OBBODo6YsyYMZJouaqbPX1r6dLccplEUFpbY7BcjiFDhuh9YfBWLX6mFw7ORdzp06dRtWpV9OzZs8BlRc2Xzc3N7fVWhf+UwjN3IClMRkYGpk2bBoVCYVARlIcPH6JTp05wd3fH1atXTTzLwuluTxNa7KEgUjjurKwszJ49G3K5HN99953JrqJ1M5ijo6NN8hqG0M03MKS5xIsXLzB06FC4urrixIkTJp6lfvbt2/d29rTOb8WLFy8wcuRIODo6Yn8hS9z61OJn78bBuRhIT0/HpEmTYGdn9869nPv27YO9vT0mTJhQ6JYscwUSQ1y+fBmenp5o06ZNodt/du7cCTs7O0yfPl0S9+t0CSn2UJjw8HBJHXdUVBR8fHwQEBCA27dvizaupTKY9aXJN9CnucSJEyfg6uqKYcOGmWWnQmGePXuG/v376509XVC+gaYWv62tLX7++WdTTblY4+BcjOguKz4NC0PG0KEIbdYMlStXxpkzZwway5SBxBjZ2dmYO3cu5HI51qxZ89bSe3JyMgIDA1GtWjWcO3fOQrMsnG6xBzGqwD19+hS9evWS5HErlUosWbIEMpkMy5cvF3yip5vBLLQ4jykVlm+QlpaGMWPGwNHR0aS7Ewxx5MgRODs7G5w9/Va+wb59eNyzJ8a4uhpVi5+9xsG5mMnKysK2Xr202ZWZpUohvYBa0QURO5CIISYmBn5+fmjSpAnur10LBAfjwowZcHBwwNixYy1accwQmmIPQqrAHThwABUqVJD8cf/9999o1KgRPvnkE9y4ccOoMTQZzJMnTy4yhSp0bxMlrVsHBAfjxqJFqFq1Kvr06SOJqnKa7GlnZ2dBO0DOnz+PYGdnbctRZZkyUEtoi11RxMG5OAoOfiO78nm/foKGEyOQiCknJwd7Bw3SnoCkWVkhev58S0/LYAUVeyjI8+fPMWjQIHz88ccWqa5mDN1CN/Pnz9f7RE9qGcyGeuuzSoQ/p0yx9LQAiJs9fePGDexQKIpFb2ypKEGs+GnVisjGhoiIlKVLU/C+fbRixQpSqVRGDefs7EwHDx6kiRMnUrt27WjatGmUmZkp5owNUrJkSfJOSSGbvP+2AajC9esWm4+xrKysqF+/fnTt2jWKi4sjT09POnv2bIHPOX78OHl4eFCZMmXo2rVr1KRJEzPNVpgSJUrQqFGj6PLly3T69Gny8/OjyMjIAp9z+PBhcnd3p/Lly9PVq1fJ39/fTLMVT8mSJck/Le31Z5WIqickWHJKlJGRQePGjaPAwEBasWIFhYWF0YcffmjUWCqVir755htq1KgRfditGyHvd4dsbHJ/h5jxLH12wExEJ7tSjGVFjQcPHqBr164W61L06tUrBAcHY5BMBqW1NUCE7DJl0KdcOcyZM6dId7LZtWtXvlXgXr16hVGjRsHJycmklanMQa1W4/vvv4dCocCXX3751oqBFDOYjaFputH3f/+DskyZNz6r5t7mpqHJnu7Zs6fg7GlNLf5GjRq9rsX/nx0gzHgcnN8TmmVFuVxu0LJifjSBxJxdis6ePYsqVaqgf//+uXssdX4INJ283N3dcfHiRbPMxxSSk5PRv3//11Xg9u/H/S5dMNze/vVxFxNJSUn47LPP4ObmhtjFi4HgYFybOxeurq4YOnSoJDKYjRUdHQ0vLy+0bds2d8+3zmf13r17aNOmTYHNJcSWmZmJ6dOni5I9rXuLYtmyZZKoi1AccXB+z9y9exetW7dG3bp1ceXKFUFj6QYSU5biy8jIwMSJE2Fvb4+9e/fm+zjdTl6TJ0+2aIs/oQ4dOoTBcvnrBBtr62J5NaJWq/HH5Mlv3JO9OHOmpadltJycHCxevLjQfd6aYjGG5hsYIzIyEh4eHujYsaNRBX103bt3Dy1atLBoLf73BQfn95BarcbmzZuhUCgwffp0wT8Mhw4dMlmXoosXL6JGjRro1q2b3rWbNZ28qlatavAWMqm4cOECfvjww/ciwSapa9c3jjPRyPKQlnbr1i188sknaNq0qd5NN0x5m0ipVGLevHlQKBTYvHmzoIIwmlsRmpU3S1Xce59wcH6P6S4rRkRECBpLt2qQGF2KsrKyMHPmTNja2mL79u1G/bDs2bMHDg4OCA4ONkuLPzHoHvfpCRO0JVbTrKwQ0qiRRZtLiC0zMxNTpkxBvw8/1N6T1ZSHHDZsWJGpv6xSqbBy5UrIZDKsWLHCqGVesW8T3bhxA/Xq1UPLli0F18x/8OABOnToAE9PT0RGRgqeG9MPB+f3nFqtRnh4OOzt7TF27FjBdX3F6FIUFRWFOnXq4NNPP0VSUpKg+aSkpGDgwIFwdXXFkSNHBI1lau887rx7lZnh4ZgwYYLFmkuI7dKlS6hVqxY+++wzPHr06K3ykCNGjICTk5PkazHfuXMHzZo1Q/369fH3338LGkuM20S6Ff7Wrl0r+HOiqcU/Y8aMIp1sWRRxcGYAcitN9e3bF5UqVRKcIWtslyKlUokFCxZALpcjLCxM1AD022+/wcXFBQMHDkRKSopo44rBkOM+f/48atWqZbbmEmLTZDArFAr8+OOPBR7ryZMnUblyZfTu3RtPnjwx4ywLZ8qmG5rbRIY0lwCAf/75Bw0bNkSjRo0E18Z/8uQJevToATc3N5w/f17QWMw4HJzZG3799Vc4OTmJsqyoCSQdOnQoNJDcvHkTfn5+aN68Oe7evSvodfOj6eTl4OCAPXv2mOQ1DBUbGwtfX1+DjjszMxOzZs0yeXMJsWkymNu0aaP3icWrV68wbtw42Nvb4+eff5bEsd6/fx9t27ZF3bp1TdZ0w5DmEmq1GmvWrBGtHv7+/ftRoUIFjB8/vkgnVRZ1HJzZW54/fy7asmJhgUSlUmH58uWQyWT49ttvzbItQ9PJq0ePHrlLqhagUqkQGhoq6LivXbsGb29v0ZtLiE3fDOaC/PXXX6hRowY6d+4s+FaHsdRqNX766SfY2tpi9uzZZtmnXNhtooSEBLRs2RL16tUTXMPg2bNnGDBgACpVqlRkEymLEw7OLF+aZcU+ffoIXlZ8K5Ds34/n/fphWq1a8Pf3x61bt0SatX40nbxsbW0LXV4V2+3bt9GkSRM0aNAAcXFxgsZSKpVYvHixaM0lxGZMBnN+MjMzMWPGDCgUCtFvexTm0aNH6NKli0WabqSlpWnzDbZv3w71vn1QBwfj+Oefi5Y9rWl8MXLkSLPVLWAF4+DMCiTmsqKmS1GfcuWQWbJkbsWk0qWRU8DeZVPTdPL69NNPkZiYaNLXUqvVWL9+PeRyOUJCQkQNpH///TcaNmwIf39/Sew/FSODOT9Xr15F3bp10bJlS8EBXx9SabqhaS6Rkbf3Pd3KCreXLxc05suXLxEUFCS48QUTHwdnGhNw6wAAIABJREFUphexlhUTExOxq0IFSe3fzcrKwldffQW5XI7169eb5IosMTERrVu3hre3N2JiYkQfH8gNiN9++y3kcjkWLFhgsS5iYmYw5yc7OxsLFiyATCbDqlWrTHI7RNN0o0qVKpJpuhHbosUb3x21gO/O6dOnUalSJQwcOLBYVZ4rLjg4M70JWVbUdGBSKBTY0bv36/27RAjv108S7SijoqLg4+ODZs2aCc521dBUglIoFGZru3n37l20atUKdevWxdWrV03+ehq6GcyLFi0yyxJ7bGws/P390bBhQ1FPBA4dOgRHR0eMGTNG8PZCMWiyp4McHZFTtmzulXOJEpjp6WlwvkF6ejrGjx+PChUqFJpsxiyHgzMzmGZZsVWrVnplGD98+BCdO3dG7dq1X5cMzdvX+ui77ywSSPKjVCoREhKirRssJMA8fPgQnTp1evO4zUStViMsLEzbXCIzM9Okr2eODOb85OTkYMWKFZDJZFi8eLGg+6+pqakYNmyYpJpuvJU9nffdydmzR5tvsGLFCr0+q+fPn4ebmxt69Oghue1p7E0cnJlRdJcVC8o23rlzJ+zs7DB16tR8A8R/uxSZOpDo49atW2jUqBHq169vVBas5rinTZtm0eNJSkpC586dUaNGDfz111+ij6/JYFYoFGbLYM7P7du3ERAQAG9vb1y7ds3g5586dQoVK1aUTNON58+fY+DAgYVmT+uTb5CVlYUvv/wStra22LFjh6mmzETEwZkJollWbNSo0RvLisnJyQgMDES1atX0DgqacqKmCiSG0u2+M2/ePL0CjzHHbWpqtRo///yzaFXgNB4/foyuXbuiZs2aZs9gzo9arcbGjRshl8sxc+ZMvU6M0tLS8Pnnn8PR0REHDx40wywLd/ToUYOyp1UqFVatWgWZTPZWvoGm8UX79u0FN75g5sPBmQmmu6y4e8AAxLdrhwHly+OLL75AWlqaQWOJXU5UDJpOXnXq1ClwefqXX36Bg4ODUcdtDk+ePEGfPn1QqVIlnDx5UtBYUslgzs+///6LDh06oFatWgVWuIqIiEDVqlXRp08fJCcnm3GG7/by5UuMHDkSzs7ORpWbvXPnDlq2bIm6devi1tKluODnhz7lyuH777+XRAEXpj8Ozkw0Dzds0G7zyC5TRlCLwydPnmjLiQoNJGLQ7eQ1bdq0NwLS8+fPMXjwYFSsWBGnTp2y3CT1pKkCN3z4cIOrwOlmMIvdRUlsarUa27Ztg52dHSZMmPDGCZOm6YadnR127dplwVm+dubMGVSqVAkDBgwQlD2tVqsR3revtg2n6oMPimW70eKOgzMTT3DwG9s8ztWrJ7hYvpBAYgoPHjxAly5d4ObmhhuLFiGhUycMlssxfPjwItP5Csg9oRg+fDicnJzw66+/6vUcTQbz6NGjJbGioa/Hjx+jV69eqFKlCqLnz8ejHj0w2sXlddMNC9PNnhba0U2lUmHZsmX4ztpaUtsVmeE4ODPx7N8P5G2RUn3wAebWq1fosqI+dAOJVLoUnRo3TntlorS2LrJXJidOnEClSpXQp08fPH369J2PkWIGszEipk7VvmfZZcpALUJrU6HEzJ6Oj49H48aN0bBhQyStW6f9LsLGpsh+Pt9nHJyZuHRa/+kuK06cOFHwfdiTJ09qA4klt4GcPXsWW/73vzeuTBI6dbLYfIR69eoVxo4d+84qcJoM5iFDhkgig9lY169fxw6F4o337E779habj5jZ02q1GmvXrn278YXOd5EVPRycmcnpLiuePn1a0FiW7FKUkZGBSZMmwd7eHhFTp2qvTJTW1hgsl2Po0KFFutJSRETE6ypw69bhtIcHBpYvr/eytxTl5ORgyZIlkMvlOBIcrC1+k2NtjaG2tujXr5/ZE8EiIyPh6ekpSva0mI0vmLRwcGZms3//fjg6OmLkyJGCr8LM3aXo4sWLqFmzJrp27fq6O5DOlcmLFy8wYsQIODo64sCBAyafj6lkZmYitFkz7fKvuggnE926dQv+/v5vNt3Qec9evnyJMWPGoEKFCmZJClMqlZg3bx7kcrng7GlNgqJcLse8efMEN75g0sPBmZnVs2fPMGTIELi4uODw4cOCxjJHl6KsrCzMnDkTCoUC27ZtK/Q1Tp06hcqVKyMwMPCdLf6kTJPBvPGDD95Y/n3Rv7+lp2YQzZ5fuVyuV9ONs2fPonr16ujatavJ9gHHxsaiXr16aNGiBe7duydorAcPHqBjx47w8PCQRFU9ZhocnJlFHDt2DBUrVkT//v0FLysaWk5UX1FRUahTpw7atWuH+/fv6/28tLQ0jB8//nWLvyKwv/Ty5cuoXbs2OnfujGdbtmiX7LNKl0afcuXM1mtbqLt376J58+YGN93IyMjA1KlToVAosGXLFtHeM032tEwmw5o1awSPGx4eDjs7O0yfPl0SlfSY6XBwZhYj5rKivuVE9aFUKrFw4ULI5XJs2rTJ6B/Uc+fOoWbNmujYsaNBwd2csrOzMWfOHCgUCvzwww+vj1Vn+Te/KnBSolsZTEjTjcuXL8PT0xNt2rQRfIWrmz0ttJHK06dP0bNnT1SvXh3nzp0TNBYrGjg4M4sTc1lRaCD5+++/Ub9+fTRv3lyUq3DdpXchgd4Url+/Dm9vb7Ru3brQXtZiNpcQ2/3799GuXTvUrVsXUVFRgsfLzs7G3LlzIZfLsWbNGoNP9PLNnjbSgQMH4ODggHHjxuU2vmDvBQ7OTBIyMjIwbdo02NraCl5WNCaQqFQqLF++3GT9gSMjI+Hl5YWWLVu+Tk6yEE0Gs0wmw4YNGwz6W8fHxyMgIAA+Pj6iBEIhNE03bG1tTdJ0IyYmBn5+fmjSpAlu3bql13MSEhLQqlUr+Pj4CM6e1jS++PjjjwXvcmBFDwdnJimaZcW2bdsiISFB0Fj6dim6ffs2mjZtCn9/f71/hI2hVCqxaNEiyGQyrFy50iL3cDUZzE2aNDG4D7CGbt/mWbNmCa4CZwxzNd3IycnR3jMOCQnJ9ypYt7zr3LlzBa8sHDt2DC4uLggKCtKr8QUrfjg4M8nRXVZcu3atoCCmuRepUCjeCiRqtRrr16+HXC7HkiVLBC8/6uvmzZto0KABGjRogJs3b5rlNXW7FumTwawPfZtLiG3Pnj1mb7rxzz//oFmzZqhXr95b/ao12dPu7u6Cs6dfvnyJUaNGGd34ghUfHJyZZOkuK8bFxQkaSzeQ/B0SgpcDB2J23brw8vLC9evXRZqx/lQqFVauXAmZTIaFCxea9B6ubgaz2CcDYleBK0hKSgr69u1rsaYbKpVKezI3Z84cZO/ahVutW6Pfhx+Kkj39xx9/oHLlyoIbX7DigYMzkzTdZcWlS5cKurpVq9U4PWGCtsBGVunSUO7eLeJsDXfnzh20aNECXl5eiIyMFHVs3QxmU58AiFkF7l0OHz4MJycnSTTdSEhIwExPT6RZWeVWGytbVlChlvT0dEyYMEGUxhes+ODgzIoEzbKir6+v0Ve6jx49wsFKld4osJHYubPIMzWcWq3Gpk2boFAoMHPmTFH2r2oymOvUqWPWxC0xq8ABuU03hg8fDhcXFxw/flyEGQp34MABbNI0ldA0PxkxwqixLly4IFrjC1a8cHBmRcZ/lxUNSUTatWsX7OzssKVrV219ZaW1NQaWLy9aIBHq/v376NSpE2rWrGn0XlbdDGZLJWuJVQVOak03dLOno+fP1xZqySxZEiMqVMCZM2f0HisrKwszZswQpfEFK544OLMiJyEhAe3atYOHh0ehmbrJycno3bs3qlatioiIiNz/qVNgQ8xyomJQq9XYsWMH7OzsMH78eIPu4epmMF+8eNGEs9SPsVXg0tLS8MUXX8DBwUEyTTeOHTsGZ2dnjBgx4nXfbp3P0Z49e+Dg4IDg4OBC+3pfu3ZNtMYXrPji4MyKJLVajR9++AG2traYMmXKO4szHDx4EI6Ojvj8888LDXJilhMVw+PHjxEYGIjKlSvj999/L/Txe/fuhb29PSZNmmS2DGZ9GFoFLiIiAtWqVUPv3r0l8T5osqednJzw22+/FfjYlJQUDBw4EK6uru/MtFYqlViwYIEojS9Y8cfBmRVpDx8+RPfu3VGtWjX88ccfwP79yBw2DCtbtICrqytOnjyp91jm7lKkjwMHDsDJyQlBQUHvXNpNSUlBv379ULlyZYtkMOursCpwmZmZmDp1Kuzs7LBz504LzPBtmuzp/v37G5Q9/dtvv8HFxQWDBg3Cy59+AoKDkbh6NXx9fUVpfMHeDxycWbGwZ88eDChfHpklS+beByxVCmnbtxs1ljm6FBni2bNnGDp0KJydnXHw4EHtcuql2bPh5OSE4OBgi2cw6+OtKnD79gHBwfgnNFTbdOPhw4eWniYyMjIEZ0+npqZibdu22p0BaUQ4FBTEV8tMbxycWbGRMXToGxm0CA42fqy8LkVilBMVy/HjxzHU1lZ7ApJmZYXIr7+29LQMdvnyZXxesSIySpTQBq6TY8dK4m984cIF1KhRA927dxecPf28Xz/RPo/s/VOCGCsmynboQGRjQ0RE6VZWtPLmTXr27JlxY5UtSwsXLqRDhw7RsmXL6NNPP6XExEQxp2uwgIAA+tLXl6xVKiIisgHI89Eji87JGF5eXjTDz4/KqtVERGRDRE2VSrKysrLYnLKzs2nmzJnUvn17mj17NoWHh5NcLjdqLAC0fv16Ct63j5SlS+f+TxsbolatRJwxK+44OLPio2NHou3biYKDibZvp1tublS7dm3at2+f0UN6e3vTxYsXyd/fn7y8vGjdunWkzgsq5pSenk5jx46lryIiKMfamoiIMqysaOHly/Tw4UOzz8dYKpWKQkJCaNyhQ6QsU4aIiDJKlKAvf/+d/vnnH4vMKSoqinx9fSkyMpIiIyOpZ8+eRo/177//Ups2bWjjxo305fnzVHrXLu3nkTp2FHHWrNiz9KU7Y6Z0+vRpVK1aFT169MCjR48EjSVmOVFD/PXXX6hWrRoCAwPx9OlT7T3nrJ07MX36dNja2mLz5s2SWBYuSFxcnLbpRnx8vPY4cvbuFa0KnCF0s6fDwsIE/f3UajW2bNkiWuMLxjg4s2IvPT0dkyZNgq2tLX766SfB7Si/+eYbswSSzMxMTJs2rdAM5itXrqBOnTpo06aNJDOBVSoVvv32W8hkMixfvjzfphtiVIHTV2xsrGjZ0w8fPkSnTp3g4eEhuPEFYxocnNl748KFC3B3d0f79u3x77//ChorLi4OTZs2NVkguXLlCtzd3fXOYM7Ozsa8efMgl8uxZs0ai7SjfJd79+4hICAAfn5+ejXdEFIFTh8qlUp7lb569WrBf6fw8HDY2dmJ0viCMV0cnNl7JSsrC1999RXkcjk2bNgg6CpapVJh3bp1ogaS7OxszJkzBwqFAlu3bjV4fjExMahfvz4aN25s0t7UhdHUCze26YYhVeD0FR8fj8aNG6Nhw4aCb0s8ffoUvXr1QvXq1Y0utcpYQTg4s/dSdHQ06tWrh+bNm+fe/xRArEASExMDb29vtG7dGomJiUaPk5OTg9DQUMhkMoSEhJj9/qdu041r164ZPY4+VeD0HUdzEiXGrYhffvkFDg4OGDdunNFzYqwwHJzZe0upVCIkJAQymQyhoaGC21EaG0hycnIQEhICuVyO9evXi5bYFR8fj+bNm6NevXqIjo4WZcyCaHo729raYubMmaItSb9VBc4AiYmJaNWqFXx8fHDjxg1B83j+/DkGDRqESpUqmaQtJmO6ODiz996tW7fQuHFjfPLJJ4J/wA0NJHFxcWjQoAGaNGmC27dvC3rtd1Gr1diwYQPkcjm++uork3Wpevz4Mbp164YaNWqYrOmGprnE6NGj8fLlywIfK3b29LFjx+Di4oKgoKBCX5sxMXBwZgy594/XrFkDuVyOefPmITs7W9B4u3fvRoUKFfINJJoMZrlcXmAGs1gSExPRvn17uLu748KFC6KOrWm6MXHiRJM33dA0l6hYsSKOHj36zsdosqfd3d0FZ0+/evVK78YXjImJgzNjOu7du4c2bdqgTp06uHLliqCxkpOTMWDAALi6uuYGkrx9vY83bjQog1ksml7PdnZ2mDRpkuD7pbpNNwxdbhZKt7lESkqK9m/7x+TJomVPG9v4gjExcHBm7D/UajU2b94MhUKB6dOnC74aPHz4MAbL5a9rYhNhZ79+FitU8ejRI/Ts2RNVq1bFmTNnjBrjt99+s3jTjdTUVIwaNQoDy5eHskwZgAjpVla4uWSJoHHFaHzBmFAcnBnLx4MHD9ClSxe4ubnhzz//FDTWsz59JNcEYe/evXBwcEBwcDBSU1P1ek5qaiqGDx8OFxcXHDt2zMQz1M/F+vVF+9uK2fiCMSG4tjZj+bC3t6fdu3fT3LlzqVu3bvTFF19QWlqaQWMAoB07dtDnv/5K2XlNENKtrGhJZCQ9fvzYFNPWW+fOnen69euUnp5O7u7udOTIkQIff/r0afL09KScnByKioqiFi1amGmm7/bixQsaPHgwrb99m1R59cbTiei0tTUBMGis7OxsmjVrliiNLxgThaXPDhgrCp4+fYp+/frh448/xvHjx/V6zpMnT9C9e/fXGcx590Uzw8MxceJEUcqJiuXIkSNwdXXFgAEDkJyc/Ma/paenY+zYsXBwcMAvv/xioRm+SZM9PWLEiNyr/ry/7a2lS1G7dm2DqsBFRUWhTp06aN++PZKSkkw8c8b0w8GZMQMcPHgQzs7OGDp0KJ4/f57v4/bu3YsKFSoUmMF84cIFgwOJKb18+RKjR4+Gg4MD9uzZA+zfjwdduyLIwQG9evXKbbphYfpkT+tbBU7MxheMiY2DM2MGevHiBYKCguDk5IQDBw688W/Pnj0zKIM5KysLs2fPFqWcqFj++OMPjKhQAeklSgBEuclW+/dbeloGZ09HRUXlWwXu5s2b8PPzQ0BAgCSbhTDGwZkxI506dQqVK1dGYGAgnm/dirsdOmCQTGZUBnNUVBR8fHxEKScq1JUrV/DjRx+9kWSltmACW0ZGBiZOnGhU9vR/q8Apd+/G1YYN0adcOVEaXzBmKhycGRMgLS0Nmzp1QpqVFUCEHGtro68ylUollixZIko5UWNkZ2fj66+/hkKhwMmxY6G2sdFuT5rj7Y2EhASzzgd4nT3drVs3QdnTt27dwvgqVZCe9z6pypaVxGoAY/nhbG3GBLCxsaHBTk5kk5cdXDIri17t3WvUWKVKlaJJkyZRREQE7dmzhxo1akSxsbFiTjdfN27cIH9/fzp79ixduXKFmoWGktX27UTBwVRq506y6tSJvLy8aP369aRWq00+H93s6VmzZgnKngZAp06dohqJifRB3vtUIjOT6OhRMafMmLgsfXbAWJG3fz+Qd5WZVbo0+v7vf9i0aZPgdpSrV6+GTCbD/PnzBZcTzU9OTg6WLl0KmUyGdevWFTjn6Oho+Pr6omnTpvjnn39MMh9A3OzpxMREtG7dGj4+Pri3apX2fYKNDV85M0nj4MyYGPK28mD/fkRGRsLb2xstW7bEnTt3BA179+5dtG7dWpRyov+labrRuHFjve9z6wbzZcuWibr0Lmb2tG7ji6+//vr1yY3O+8SYlHFwZswElEolFi1aBJlMhpUrVwpKPBK7nKjuVXloaKhRc4uLi0OTJk1Qv359xMTECJoPIG72tG7jC7FPaBgzFw7OjJnQzZs30aBBAzRo0EBwk4ukpCR89tlncHNzQ0REhFFj3Lt3DwEBAfD19UVsbKyg+ahUKqxduxZyuRxz5841auldpVIhNDQUcrlclOzpnTt3ws7ODtOmTRPc+IIxS+LgzJiJqVQqrFq1CjKZDAsXLhTU8EKtViM8PBz29vb44osv9N6ypVarERYWBrlcjvnz54vadEPTycvT0xOXL1/W+3m3b99GkyZN0KBBA8TFxQmaQ3JyMgIDA1GtWjWcO3dO0FiMSQEHZ8bM5M6dO2jRogW8vLwQGRkpaKynT5+ib9++epUTTUpKQvv27eHp6Ylr164Jet386N7jnTp1aoFL72q1GuvXr4dcLkdISIjg+9a//vorHBwcMHbsWKSlpQkaizGp4ODMmBmp1Wps2rQJCoUCM2bMELz0qiknOmzYsHeWE92xYwdsbW0xY8YMZGVlCXotfTx48ABdu3ZF9erVcfbs2bf+XTd7Wui96ufPn2Pw4MGoWLEiTp06JWgsxqTGCjCwfQtjTLCkpCQaNWoUxcXFUVhYGPn5+Rk9VmpqKk2ePJkOHjxIa9eupfZqNWUcOEArb96kzcnJtHXrVqpXr56Isy/crl27aMyYMdSjRw9a5O9PZf/4g34vXZp6/vQTjRkzhqZOnUql87p0GePEiRM0ePBgatu2LYWEhFC5cuVEnD1jEmDpswPG3ldqtRo7duyAnZ0dxo8fL3hJ9uTJkxhub4+MvJrYWaVKITM8XKTZGu7p06dY1rSptnpaupUV/gkNFTTmq1evEBwcDCcnJxw+fFikmTImPVwhjDELsbKyop49e1J0dDQ9ePCAPDw86Pfffzd6vLp161KXcuWobF4FrzI5OVTm9GmRZms4mUxGHayttdXTPgDIMSbG6PH+/PNPqlOnDr18+ZKio6OpTZs2Yk2VMcnh4MyYhSkUCtq2bRuFhoZS3759KSgoiFJTUw0a4+jRo+Th4UF3qlQhfPABERFlWFnRgosXKSkpyRTTLlBKSgr17t2bvomOJlXZskRElFWqFI3cs4cOHDhg0FiZmZk0efJk6tatG4WEhNCWLVvo//7v/0wxbcYkg4MzYxLRoUMHun79OqlUKqpduzYdOnSo0Oe8evWKRo4cSUOHDqVNmzZR0KFDZLVjh7YmdmarVlSnTh0KCwsjmCm95ODBg+Tu7k52dnb0TVwclfz5Z6LgYLLevZsG7N5N48ePp8DAQHry5EmhY126dIm8vb3p9u3bFBUVRZ07dzbDETAmAZZeV2eMve348eP4+OOP0bdvXzx9+vSdjzl9+jQqVaqEQYMGvTNTWyMyMhJeXl6ilBMtyIsXL/TKnk5LS8OECRNgZ2eHbdu2vbNMZ1ZWFmbOnAlbW9t8H8NYccbBmTGJevXqFb744gvY29tj586d2rrQmeHhGDduHBwcHPDLL7/oNZZSqcTChQshk8mwatUq0fsYHz9+HK6urhgxYgRSU1P1es758+dRq1YtdOjQAf/++6/2/2saX3z66ae4f/++qPNkrKjgrVSMSVxERAT90KMHhT58SGVVKsqwsqLVDRvSoL17SSaTGTTWzZs3aejQoWRlZUWbNm2iatWqCZpbWloaTZ06lfbt20ffffedwUlaWVlZtGDBAlqzZg390L07ya9epeUxMRSwYgUNHDiQrKysBM2PsaKK7zkzJnH+/v60qkMHKqtSEVFu1vMEDw+DAzMRkZubG505c4Z69OhB/v7+tHjxYsrJyTFqXprs6RcvXlBUVJRR2dPW1tY0Z84c+nX4cGq8di35nDtHW3NyaJBMxoGZvdc4ODNWBJRq25bIxoaIcrOw50REUEJCglFjlShRgsaMGUOXLl2i48ePU/369SkqKkrv5+tmTy9ZsoS2bt1KH330kVFzUavVtGLFCopetoxsNPPLyCA6etSo8RgrLjg4M1YUdOxItH07UXAwld61i0p37UpeXl60Zs0aUuftazZUxYoV6ejRozRy5EgKCAigWbNmUVZWVoHP+W/29GeffWbUaxMR3blzhwICAig8PJzahYZqTz7IxoaoVSujx2WsOOB7zowVUTdu3KAhQ4ZQmTJlaOPGjVS1alWjx0pKSqKRI0dSfHw8hYWFka+v7xv/np2dTfPnz6d169bR8uXLqVevXkYvOwOg7777jr788kuaMmUKjRs3jkqWLEl04EDuFXOrVrknI4y9xzg4M1aEqVQqWrVqFc2bN08b6EqVKmXUWADo559/prFjx1Lfvn1pvp8fWZ8+TXerVaMumzeTg4MDbdiwgRwcHIye7/3792nIkCH09OlT2rJlC9WqVcvosRgrzjg4M1YMxMfH07Bhw+jly5cUFhZG7u7uRo/15MkT2tK1K406e5ZsAEonoj/HjKEWK1YIulr+6aefaPz48aI0vmCsuDPuFJsxJimVK1emEydO0MaNG6l58+YUHBxM06dPpzJlyhg8lkKhoCEuLtqa2DZE1Dgz0+jA/PjxYwoKCqK4uDg6cuQI1a1b16hxGHufcEIYY8WElZUVDRs2jK5evUqXLl0iHx8funjxokFjaLKnxxw4QMq8K9uskiVp1N69dPjwYYPntHv3bvLw8KDq1avTpUuXODAzpide1masGAJA27Zto/Hjx9OAAQNozpw59EFeQ4z83LlzhwYPHkzZ2dm0efNmqhobq03QOm5jQ8OGDaNGjRrR8uXLqXz58gWOlZKSot2utWXLFqpfv76Yh8dYscdXzowVQ1ZWVtSnTx+Kjo6me/fukYeHB505c+adj9VkT/v6+lK7du3ozJkzuZnfHTsSffstUceO1KJFC4qOjqby5ctT7dq1adeuXfm+9qFDh8jDw4MUCgVdvXqVAzNjRuArZ8beA/v27aPg4GDq3LkzLVq0iMqVK0dEudnTQ4cOpcePH9PWrVv1yp6OiIigwYMHU61atWj16tVkb29PRESpqak0fvx4OnHiBIWFhVGzZs1MekyMFWd85czYe6Bz5850/fp1ysjIoNq1a9Plr76imy1b0tSaNal+/fp07tw5vbc1+fv7U2RkJFWvXp08PT3p5NixlNi5M42vUoVKlChBUVFRHJgZE4ivnBl7z1ydM4fc5syhDwBSWVtTyfBwo4t+xC5eTK5Tp5INEeVYW1MpAWMxxl7jK2fG3jN1nzyhD/LOyUtmZdE/a9caNU5ERARdWrhQWxO7VFYW18RmTCQcnBl737Rqpa1jrbK2pqXXrlH37t3p0aNHej09MzOTpkypvzttAAACE0lEQVSZQl27dqUqI0dyTWzGTICDM2PvG50mGiXDwyk0Pp4qVapEHh4e9OOPP1JBd7ouX75M3t7eFB8fT1FRUfTJwoXasWj7dl7SZkwkfM+ZMUZEuR2nBg8eTM7OzrRu3TpydnbW/ptSqaT58+fT2rVrBTe+YIwVjq+cGWNEROTj40OXLl0iPz8/8vLyovXr15N63z5K7t2bJrm50cWLF+nq1asUGBjIgZkxE+MrZ8bYW65fv07ff/YZzbt9mz5Qq0lZpgyVCg8nq06dLD01xt4LfOXMGHtL7dq1KaRVK/pArSYiotLZ2WR17JiFZ8XY+4ODM2PsnUq0bs2Z2IxZCC9rM8byd+CAtvkFZ2IzZj4cnBljjDGJ4WVtxhhjTGI4ODPGGGMSw8GZMcYYkxgOzowxxpjEcHBmjDHGJIaDM2OMMSYxHJwZY4wxieHgzBhjjEkMB2fGGGNMYjg4M8YYYxLDwZkxxhiTGA7OjDHGmMRwcGaMMcYkhoMzY4wxJjEcnBljjDGJ4eDMGGOMSQwHZ8YYY0xiODgzxhhjEsPBmTHGGJMYDs6MMcaYxHBwZowxxiSGgzNjjDEmMRycGWOMMYnh4MwYY4xJDAdnxhhjTGI4ODPGGGMSw8GZMcYYkxgOzowxxpjEcHBmjDHGJIaDM2OMMSYxHJwZY4wxieHgzBhjjEkMB2fGGGNMYjg4M8YYYxLz/3hMvEjMYrUmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import networkx as nx\n",
    "n = 10\n",
    "ex = np.ones(n);\n",
    "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "e = sp.sparse.eye(n)\n",
    "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
    "A = csc_matrix(A)\n",
    "G = nx.Graph(A)\n",
    "nx.draw(G, pos=nx.spectral_layout(G), node_size=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fill-in\n",
    "\n",
    "The fill-in of a matrix are those entries which change from an initial zero to a nonzero value during the execution of an algorithm.\n",
    "\n",
    "The fill-in is different for different permutations. So, before factorization we need to find reordering which produces the smallest fill-in.\n",
    "\n",
    "**Example**\n",
    "\n",
    "$$A = \\begin{bmatrix} * & * & * & * & *\\\\ * & * & 0 & 0 & 0 \\\\ * & 0  & * & 0 & 0 \\\\ * & 0 & 0& * & 0 \\\\ * & 0 & 0& 0 & * \\end{bmatrix} $$\n",
    "\n",
    "If we eliminate elements from the top to the bottom, then we will obtain dense matrix.\n",
    "However, we could maintain sparsity if elimination was done from the bottom to the top."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Fill-in minimization\n",
    "\n",
    "Reordering the rows and the columns of the sparse matrix in order to reduce the number of nonzeros in $L$ and $U$ factors is called **fill-in** minimization. <br> \n",
    "\n",
    "Examples of algorithms for fill-in minimization:\n",
    "\n",
    "- **Minimum degree ordering** - order by the degree of the vertex\n",
    "- **Cuthill–McKee algorithm** (and reverse Cuthill-McKee) - reorder to minimize the bandwidth (does not exploit graph representation).\n",
    "- **Nested dissection**: split the graph into two with minimal number of vertices on the separator (set of vertices removed after we separate the graph into two distinct connected graphs). <br> Complexity of the algorithm depends on the size of the graph separator. For 1D Laplacian separator contains only 1 vertex, in 2D - $\\sqrt{N}$ vertices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Take home message\n",
    "\n",
    "- CSR format for storage\n",
    "- Cache and parallel issues in sparse matrix processing\n",
    "- Reordering and blocking as a way to solve these issues\n",
    "- Sparse matrices & graphs ordering\n",
    "- Ordering is important for LU fill-in: more details in the next lecture"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}