lectures/lecture-11/lecture-11.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 11: From dense to sparse linear algebra"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "Randomized matmul, trace, linear least squares."
   ]
  },
  {
   "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"
   ]
  },
  {
   "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": 1,
   "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",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x10415f880>"
      ]
     },
     "execution_count": 1,
     "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": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x10ebb9d00>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAAGiCAYAAAB+sGhNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAvs0lEQVR4nO3dfXRU5bn38V9CSAjiJAJmkmhiKb4Ax1jeNEZRew5ZIqDngLQHNAK1Ao80WFC0yvER26MSoNUqPQpFq3AqQut6xCIV2jQoqMQAURQQEas1qE1CG5MJWsJL7ucPmjEDyWR2sve87Pl+1prVMvvOzp69ID8n131dk2CMMQIAwMUSI30BAAA4jbADALgeYQcAcD3CDgDgeoQdAMD1CDsAgOsRdgAA1yPsAACuR9gBAFyPsAMAuJ5rw+7xxx/XN77xDfXo0UP5+fnatm1bpC/JcSUlJbr44ot1+umnKyMjQ+PGjdO+ffsC1hw+fFjFxcXq06ePevXqpQkTJqimpiZgTVVVlcaOHauePXsqIyNDd911l44dOxbOlxIWCxcuVEJCgubMmeN/Lt7vz2effaabbrpJffr0UWpqqvLy8rRjxw7/cWOM5s+fr6ysLKWmpqqwsFD79+8POEddXZ2Kiork8XiUnp6uW265RYcOHQr3S3HE8ePHdd9996lfv35KTU1V//799cADD6j11MV4ukdbtmzRddddp+zsbCUkJOjFF18MOG7XvXj33Xd1xRVXqEePHsrJydHixYutX6xxoTVr1pjk5GTz9NNPmz179pjp06eb9PR0U1NTE+lLc9SoUaPMM888Y3bv3m127txpxowZY3Jzc82hQ4f8a2699VaTk5NjysrKzI4dO8yll15qLrvsMv/xY8eOmQsvvNAUFhaat99+27z88sumb9++Zt68eZF4SY7Ztm2b+cY3vmEuuugiM3v2bP/z8Xx/6urqzDnnnGO+973vmYqKCvPRRx+ZP/zhD+bDDz/0r1m4cKFJS0szL774onnnnXfMv//7v5t+/fqZf/zjH/4111xzjfnWt75l3nzzTfPaa6+Zc88919xwww2ReEm2e+ihh0yfPn3M+vXrzccff2yef/5506tXL/PYY4/518TTPXr55ZfNvffea1544QUjyaxduzbguB33oqGhwXi9XlNUVGR2795tVq9ebVJTU80vf/lLS9fqyrC75JJLTHFxsf/Px48fN9nZ2aakpCSCVxV+tbW1RpLZvHmzMcaY+vp60717d/P888/71+zdu9dIMuXl5caYE395ExMTTXV1tX/N0qVLjcfjMU1NTeF9AQ5pbGw05513niktLTVXXXWVP+zi/f7cfffdZsSIEe0eb25uNpmZmeanP/2p/7n6+nqTkpJiVq9ebYwx5r333jOSzPbt2/1rNmzYYBISEsxnn33m3MWHydixY833v//9gOeuv/56U1RUZIyJ73t0ctjZdS+eeOIJc8YZZwT8+7r77rvNBRdcYOn6XPdrzCNHjqiyslKFhYX+5xITE1VYWKjy8vIIXln4NTQ0SJJ69+4tSaqsrNTRo0cD7s2AAQOUm5vrvzfl5eXKy8uT1+v1rxk1apR8Pp/27NkTxqt3TnFxscaOHRtwHyTuz7p16zR8+HB997vfVUZGhoYMGaInn3zSf/zjjz9WdXV1wP1JS0tTfn5+wP1JT0/X8OHD/WsKCwuVmJioioqK8L0Yh1x22WUqKyvTBx98IEl655139Prrr2v06NGSuEet2XUvysvLdeWVVyo5Odm/ZtSoUdq3b5+++OKLkK8nqasvKNr87W9/0/HjxwN+GEmS1+vV+++/H6GrCr/m5mbNmTNHl19+uS688EJJUnV1tZKTk5Wenh6w1uv1qrq62r+mrXvXcizWrVmzRm+99Za2b99+yrF4vz8fffSRli5dqjvuuEP/9V//pe3bt+uHP/yhkpOTNXXqVP/ra+v1t74/GRkZAceTkpLUu3fvmL8/knTPPffI5/NpwIAB6tatm44fP66HHnpIRUVFksQ9asWue1FdXa1+/fqdco6WY2eccUZI1+O6sMMJxcXF2r17t15//fVIX0rUOHDggGbPnq3S0lL16NEj0pcTdZqbmzV8+HAtWLBAkjRkyBDt3r1by5Yt09SpUyN8ddHht7/9rVatWqXnnntO//Iv/6KdO3dqzpw5ys7O5h5FOdf9GrNv377q1q3bKTvoampqlJmZGaGrCq9Zs2Zp/fr1euWVV3T22Wf7n8/MzNSRI0dUX18fsL71vcnMzGzz3rUci2WVlZWqra3V0KFDlZSUpKSkJG3evFlLlixRUlKSvF5vXN+frKwsDRo0KOC5gQMHqqqqStLXry/Yv63MzEzV1tYGHD927Jjq6upi/v5I0l133aV77rlHkyZNUl5eniZPnqzbb79dJSUlkrhHrdl1L+z6N+e6sEtOTtawYcNUVlbmf665uVllZWUqKCiI4JU5zxijWbNmae3atdq0adMpb/2HDRum7t27B9ybffv2qaqqyn9vCgoKtGvXroC/gKWlpfJ4PKf8IIw1I0eO1K5du7Rz507/Y/jw4SoqKvL//3i+P5dffvkprSoffPCBzjnnHElSv379lJmZGXB/fD6fKioqAu5PfX29Kisr/Ws2bdqk5uZm5efnh+FVOOurr75SYmLgj81u3bqpublZEveoNbvuRUFBgbZs2aKjR4/615SWluqCCy4I+VeYktzbepCSkmJWrFhh3nvvPTNjxgyTnp4esIPOjWbOnGnS0tLMq6++av7617/6H1999ZV/za233mpyc3PNpk2bzI4dO0xBQYEpKCjwH2/ZWn/11VebnTt3mo0bN5ozzzzTFVvr29J6N6Yx8X1/tm3bZpKSksxDDz1k9u/fb1atWmV69uxpnn32Wf+ahQsXmvT0dPO73/3OvPvuu+Y//uM/2txKPmTIEFNRUWFef/11c95558Xktvq2TJ061Zx11ln+1oMXXnjB9O3b1/zoRz/yr4mne9TY2Gjefvtt8/bbbxtJ5pFHHjFvv/22+eSTT4wx9tyL+vp64/V6zeTJk83u3bvNmjVrTM+ePWk9aPGLX/zC5ObmmuTkZHPJJZeYN998M9KX5DhJbT6eeeYZ/5p//OMf5gc/+IE544wzTM+ePc348ePNX//614Dz/OUvfzGjR482qamppm/fvmbu3Lnm6NGjYX414XFy2MX7/XnppZfMhRdeaFJSUsyAAQPM8uXLA443Nzeb++67z3i9XpOSkmJGjhxp9u3bF7Dm73//u7nhhhtMr169jMfjMTfffLNpbGwM58twjM/nM7Nnzza5ubmmR48e5pvf/Ka59957A7bFx9M9euWVV9r8mTN16lRjjH334p133jEjRowwKSkp5qyzzjILFy60fK0JxrRq/QcAwIVcV7MDAOBkhB0AwPUIOwCA6xF2AADXI+wAAK5H2AEAXI+wAwC4XlSHXVc+bbypqUk//vGP1dTU5OAVxi7uT3Dcn45xj4Lj/gQX7vsTtU3lv/nNbzRlyhQtW7ZM+fn5evTRR/X8889r3759p3wkRFt8Pp/S0tLU0NAgj8cThiuOLdyf4Lg/HeMeBcf9CS7c9ydq39k98sgjmj59um6++WYNGjRIy5YtU8+ePfX0009H+tIAADEmKsOOTxsHANgpKj+8tTOfNt7U1BTwu9+WzyRraGhw7Dpjmc/nC/hfBOL+dIx7FBz3J7iO7o8xRo2NjcrOzj7lY5U6IyrDrjNKSkr0k5/85JTnc3NzI3A1sSMnJyfSlxDVuD8d4x4Fx/0JrqP7c+DAgYAPoe6sqAy7znza+Lx583THHXf4/9zQ0KDc3FydNXOFElN6Bqx9espwXdK/j//Pg3/yRx1rbnufzgszC3R+JsVlAAgnn8+nnJwcnX766bacLyrDrvWnjY8bN07S1582PmvWrDa/JiUlRSkpKac8n5jS85Swm/ab97R6Wr4Kzu0rSVo189u64amKNs/7naffCVgLAAifhIQEW84TlRtUJOmOO+7Qk08+qZUrV2rv3r2aOXOmvvzyS9188822nP+GpypU/uHfJEkF5/bV8HPa/3j31msBALEnasNu4sSJ+tnPfqb58+dr8ODB2rlzpzZu3HjKppWuaB1iS28aJk+P9t/o3vBUhd7/K4VmAIhFUdtU3lUtDYtLfv+2Ht7yWdC1Lb+mPNjYpFt/XanKqi/aXbtx9hUakEUNDwCcZHfTedS+s7PL1BHf1P8dPSDompZ3eGeenqL/94PLtHpafrtrr3nsNX6lCQAxxvVhJ0nTruofcuBJJ2p4wQKPGh4AxJa4CDvpROAFCzDpRIgdbDzRmB7KphVqeAAQG1wfdv/n1zsCAqyjd3j5C/7kD7GlNw3TsNz2A++ax14j8AAgBrg+7N748O8qKCnzh1JHv9JsNl/X5UKt4RF4ABDdXB92knSs2QSEkt01PDatAEB0i4uwa3Fy4J2f0Svoemp4AOAOcRV2UmDgrZp+qQZmBp+7dvFDfwq58ZxfaQJAdIq7sJO+DqUzT0/RhjlXWurDK5v7bTatAECMicuwk079lSaN5wDgXq4PuyE56e0eax1KodbwGB4NALHH9WH384mDQx7wvGr6pTrvzNOCno/h0QAQe1wfdn3/WWcLZWPJmaenqHTut6nhAYDLuD7sJFkOJWp4AOAucRF2kixPQwk18FrWMzwaAKJX3IRdi442lpwceB0Njx6z5DUazwEgysVd2EnWmsM7Gh7dbE5tPKeGBwDRJS7D7kwLm1YkZ2p4BB4AhE9chp3k7KYVieHRABBN4jbspK83rVip4Q3K8gQ9J8OjASD6xHXYtbBSw1v5/Us6bDxneDQARBfCTtZqeDSeA0DsIez+icZzAHAvwq6VUGt4DI8GgNhC2LXByoBnhkcDQPQj7NpADQ8A3IWwa4dTNbyWgKSGBwDhQ9gF4cTw6NYhxvBoAAgPwi4Edg+PpvEcAMKLsAuRncOjJYZHA0A4EXYh6szw6EXj84Kek+HRABAehJ0FVjetTMzPZXg0AEQBws6izgyPDqXxnBoeADiHsOskKzW8UBrPGR4NAM4h7DqJxnMAiB2EXRcwPBoAYgNh10UMjwaA6EfY2YTh0QAQvQg7m1DDA4DoRdjZKBpqeAQeAJyKsLOZE8OjaTwHgK4h7BzC8GgAiB6EnYOsDo9eeuPQoOej8RwAOoewc5DV4dGjL8pi0woAOICwc1g0bFqhhgcg3hF2YUDjOQBEFmEXRjSeA0BkEHZhROM5AEQGYRdm1PAAIPwIuwiwGkqhBl7rNoZg56aGByDeEHYRZKU5PJTG8zFLXqPxHADaQNhFmNXG82Dv8JrNqY3n1PAAgLCLuGio4RF4ANyOsIsCnRkevWh8XtBzMjwaAL5G2EURK8OjJ+bnhtR4Tg0PAAi7qGOlhhdK4znDowGAsIs6NJ4DgP0IuygUDZtWqOEBcBPCLkoxPBoA7EPYRTmGRwNA1xF2UY4aHgB0HWEXA6jhAUDX2B52JSUluvjii3X66acrIyND48aN0759+wLWHD58WMXFxerTp4969eqlCRMmqKamJmBNVVWVxo4dq549eyojI0N33XWXjh07ZvflxozONJ6HGngSw6MBuJvtYbd582YVFxfrzTffVGlpqY4ePaqrr75aX375pX/N7bffrpdeeknPP/+8Nm/erM8//1zXX3+9//jx48c1duxYHTlyRFu3btXKlSu1YsUKzZ8/3+7LjTlWGs9DGR5N4zmAeJBgjDFOfoODBw8qIyNDmzdv1pVXXqmGhgadeeaZeu655/Sd73xHkvT+++9r4MCBKi8v16WXXqoNGzbo2muv1eeffy6v1ytJWrZsme6++24dPHhQycnJHX5fn8+ntLQ0NTQ0yOPxOPkSw+5gY5NGPvyqfIfbf6e7cfYVGpB14nU/tfnPenDD+0HPuXpavgrO7auDjU269deVqqz6IqRzA4AT7P4Z7njNrqGhQZLUu3dvSVJlZaWOHj2qwsJC/5oBAwYoNzdX5eXlkqTy8nLl5eX5g06SRo0aJZ/Ppz179rT5fZqamuTz+QIebmVl04p04h3e0huHBj0nw6MBuJmjYdfc3Kw5c+bo8ssv14UXXihJqq6uVnJystLT0wPWer1eVVdX+9e0DrqW4y3H2lJSUqK0tDT/Iycnx+ZXE12sbloZfVGWrTU8Nq0AiCWOhl1xcbF2796tNWvWOPltJEnz5s1TQ0OD/3HgwAHHv2ekhdp43vodHsOjAcQjx8Ju1qxZWr9+vV555RWdffbZ/uczMzN15MgR1dfXB6yvqalRZmamf83JuzNb/tyy5mQpKSnyeDwBj3jB8GgACM72sDPGaNasWVq7dq02bdqkfv36BRwfNmyYunfvrrKyMv9z+/btU1VVlQoKCiRJBQUF2rVrl2pra/1rSktL5fF4NGjQILsvOebReA4Awdm+G/MHP/iBnnvuOf3ud7/TBRdc4H8+LS1NqampkqSZM2fq5Zdf1ooVK+TxeHTbbbdJkrZu3SrpROvB4MGDlZ2drcWLF6u6ulqTJ0/WtGnTtGDBgpCuw827MdtjdSellV2aklT+4d90w1MVIa0FgK6w+2e47WGXkJDQ5vPPPPOMvve970k60VQ+d+5crV69Wk1NTRo1apSeeOKJgF9RfvLJJ5o5c6ZeffVVnXbaaZo6daoWLlyopKT23720Fo9h1+I7S7dqxyftB17rULr6kc36oPZQ0PO1Xm/l3ADQWVEfdtEinsPOSh/ewcYm3bi8XPsPftnu2pPXW+nxA4DOiLk+O4SfEzW8llYDangAYhFh51IMjwaArxF2LsbwaAA4gbCLAwyPBhDvCLs4YaU5vODcvlo0Pi/o+U5uPKeGByCaEXZxwurw6In5ubbX8Ag8AJFC2MURJzetSAyPBhC9CLs4w/BoAPGIsItTDI8GEE8IuzjF8GgA8YSwi2M0ngOIF4RdnAu1htcSSqHW8FpvWumohkfgAXAaYQdJHdfZWm8sCaWG1zrErJwbAJxA2EESNTwA7kbYwY8aHgC3IuwQoDPDo5feODToORkeDSDSCDu0ycrw6NEXZTE8GkBUI+zQLqvDozv6lSbDowFECmGHdlkdHu1EDY/AA2AHwg5BMTwagBsQdugQw6MBxDrCDiFjeDSAWEXYIWRONZ63rKfxHIBTCDtY4kQNr6UuR+M5AKcQdrCM4dEAYg1hh05jeDSAWEHYodMYHg0gVhB26JLO1PAWjc8Lek4azwHYjbBDl1kNpYn5uTSeAwgrwg62sTI8etpV/RkeDSBsCDvYyunh0TSeA+gMwg62cnp4NJtWAHQGYQfb8YnnAKINYQdH0HgOIJoQdnCUE43nLetpPAcQKsIOjnKi8bz1LE1qeABCQdjBcdTwAEQaYYewsBpKdn/iOTU8IL4RdggrK83h067qr42zrwh6PhrPAYSCsEPYWWkOH5Dlsdx4Tg0PwMkIO4RdNNTwCDwgvhB2iAiroWR3DY9NK0B8IewQUVaHR4fSeE4ND8DJCDtEnJUaXiiN5wyPBnAywg4RxyeeA3AaYYeo4NSmlZaApPEciG+EHaKGE8OjW69neDQQvwg7RB0nhkeHWsNj0wrgToQdog41PAB2I+wQlaKh8ZxfaQLuQdghanWm8XzR+Lyg52R4NBCfCDtEPSuN5xPzc4MGmETjORCPCDvEBCvN4QXn9mV4NIAAhB1igpVNKxLDowEEIuwQM5zctCIxPBpwM8IOMSXUxnOGRwNojbBDTGJ4NAArCDvEJKcaz1vW03gOuAthh5jlRA2vpS5H4zngLoQdYpoTw6NP3rTC8Ggg9hF2cAWGRwMIxvGwW7hwoRISEjRnzhz/c4cPH1ZxcbH69OmjXr16acKECaqpqQn4uqqqKo0dO1Y9e/ZURkaG7rrrLh07dszpy0WM6kwNL9TRYtTwgNjnaNht375dv/zlL3XRRRcFPH/77bfrpZde0vPPP6/Nmzfr888/1/XXX+8/fvz4cY0dO1ZHjhzR1q1btXLlSq1YsULz58938nIR46yG0sT8XIZHA3HCsbA7dOiQioqK9OSTT+qMM77+4dPQ0KBf/epXeuSRR/Rv//ZvGjZsmJ555hlt3bpVb775piTpj3/8o9577z09++yzGjx4sEaPHq0HHnhAjz/+uI4cOeLUJcMFOjM82s7Gc2p4QHRyLOyKi4s1duxYFRYWBjxfWVmpo0ePBjw/YMAA5ebmqry8XJJUXl6uvLw8eb1e/5pRo0bJ5/Npz549bX6/pqYm+Xy+gAfil5Xh0dOu6s/waMDlHAm7NWvW6K233lJJSckpx6qrq5WcnKz09PSA571er6qrq/1rWgddy/GWY20pKSlRWlqa/5GTk2PDK0EsY3g0gBa2h92BAwc0e/ZsrVq1Sj169LD79O2aN2+eGhoa/I8DBw6E7XsjOjE8GkAL28OusrJStbW1Gjp0qJKSkpSUlKTNmzdryZIlSkpKktfr1ZEjR1RfXx/wdTU1NcrMzJQkZWZmnrI7s+XPLWtOlpKSIo/HE/AAGB4NQHIg7EaOHKldu3Zp586d/sfw4cNVVFTk///du3dXWVmZ/2v27dunqqoqFRQUSJIKCgq0a9cu1dbW+teUlpbK4/Fo0KBBdl8yXI7h0QASjDHG6W/y7W9/W4MHD9ajjz4qSZo5c6ZefvllrVixQh6PR7fddpskaevWrZJOtB4MHjxY2dnZWrx4saqrqzV58mRNmzZNCxYsCOl7+nw+paWlqaGhgXd5kCQdbGzSyIdfle9w+/2aG2dfoQFZHh1sbNKNy8u1/+CXQc+5elq+Cs7ta+ncADpm98/wiExQ+fnPf65rr71WEyZM0JVXXqnMzEy98MIL/uPdunXT+vXr1a1bNxUUFOimm27SlClT9N///d+RuFy4hFPDo2k8B6JfWN7ZRQLv7NCeg41NuvXXlaqs+qLdNa3fhT21+c96cMP7Qc/Z8g5Pkso//JtueKoipLUA2uaKd3ZAJDE8Gog/hB3iltXh0RedlRb0fAyPBqIXYYe4ZbWGt+62EdTwgBhF2CGuOdmHx/BoIHoQdoh7DI8G3I+wA/6J4dGAexF2QCsMjwbcibADWmF4NOBOhB1wEqc2rbR+R8imFSC8CDugDU4Mjx6z5DVqeECEEHZAEFZqeKumX6rzzjyt3bXN5tQaXqjnBtA1hB0QBMOjAXcg7IAOdKaGt2h8XtBz0ngOhBdhB4TA6vDoifm5DI8GoghhB1hgdXh0sBpey3qGRwPOI+wAC6jhAbGJsAMsiobh0QQeYA1hB3RCpIdHs2kFsIawA7qA4dFAbCDsgC5yeng0jedA1xF2QBc5PTyaTStA1xF2gA2iYdMKNTygfYQdYBOrjeehDI+m8RywB2EH2IzGcyD6EHaAzTrTeL70xqFBz0kND+gawg5wgNVQGn1RFjU8wEGEHeAQq6Fkd+M5NTzga4Qd4DArzeE0ngPOIOyAMHC68ZwaHhAcYQeEQTT04RF4iGeEHRAmDI8GIoewA8LM6vDoUBrPQz03NTzEK8IOiAArNbxQGs+/Oe/3DI8GgiDsgAiw+xPPmw2N50AwhB0QIdGwaYUaHuIFYQdEUGeGRy8anxf0nAyPBk5F2AFRwMqA54n5uZYCj+HRAGEHRAWrHwA7MT+XD4AFLCDsgChBDQ9wDmEHRJFIN55Tw4NbEXZAFLLaeM7waCA4wg6IUgyPBuxD2AFRyuqmFYZHA+0j7IAo5uSmFYnh0YgfhB0Q5UJtPLc6PJoaHuIJYQfECLuHR59cw2N4NNyMsANihN3DoyUazxE/CDsghnSmhrf0xqFBz0njOeIBYQfEGKvDo0dflBVSDY/h0XAzwg6IUVYGPIdSw2N4NNyMsANiFDU8IHSEHRDDGB4NhIawA2Icw6OBjhF2gEswPBpoH2EHuIjdw6PzF/zJv57h0YhlhB3gInYPj242X9flGB6NWEbYAS7D8GjgVIQd4EIMjwYCEXaAizE8GjiBsANcrDON54vG5wU9J43niEWEHeByVkNpYn4ujedwHUfC7rPPPtNNN92kPn36KDU1VXl5edqxY4f/uDFG8+fPV1ZWllJTU1VYWKj9+/cHnKOurk5FRUXyeDxKT0/XLbfcokOHDjlxuYDrWR0eHWoNj+HRiBW2h90XX3yhyy+/XN27d9eGDRv03nvv6eGHH9YZZ3z9D2Hx4sVasmSJli1bpoqKCp122mkaNWqUDh8+7F9TVFSkPXv2qLS0VOvXr9eWLVs0Y8YMuy8XiCsMj0a8SjDGGDtPeM899+iNN97Qa6+91uZxY4yys7M1d+5c3XnnnZKkhoYGeb1erVixQpMmTdLevXs1aNAgbd++XcOHD5ckbdy4UWPGjNGnn36q7OzsDq/D5/MpLS1NDQ0N8ng89r1AIMYdbGzSyIdfle/wsXbXbJx9hQZknfh389TmP+vBDe8HPefqafkqOLevDjY26dZfV6qy6ouQzg20x+6f4ba/s1u3bp2GDx+u7373u8rIyNCQIUP05JNP+o9//PHHqq6uVmFhof+5tLQ05efnq7y8XJJUXl6u9PR0f9BJUmFhoRITE1VRUdHm921qapLP5wt4ADgVw6MRj2wPu48++khLly7Veeedpz/84Q+aOXOmfvjDH2rlypWSpOrqakmS1+sN+Dqv1+s/Vl1drYyMjIDjSUlJ6t27t3/NyUpKSpSWluZ/5OTk2P3SANdgeDTije1h19zcrKFDh2rBggUaMmSIZsyYoenTp2vZsmV2f6sA8+bNU0NDg/9x4MABR78f4AYMj0a8sD3ssrKyNGjQoIDnBg4cqKqqKklSZmamJKmmpiZgTU1Njf9YZmamamtrA44fO3ZMdXV1/jUnS0lJkcfjCXgA6Jjdw6NPbjynDw/RwPawu/zyy7Vv376A5z744AOdc845kqR+/fopMzNTZWVl/uM+n08VFRUqKCiQJBUUFKi+vl6VlZX+NZs2bVJzc7Py84P/lyUAa+weHi1Zr+EReHCa7WF3++23680339SCBQv04Ycf6rnnntPy5ctVXFwsSUpISNCcOXP04IMPat26ddq1a5emTJmi7OxsjRs3TtKJd4LXXHONpk+frm3btumNN97QrFmzNGnSpJB2YgKwhuHRcDvbw+7iiy/W2rVrtXr1al144YV64IEH9Oijj6qoqMi/5kc/+pFuu+02zZgxQxdffLEOHTqkjRs3qkePHv41q1at0oABAzRy5EiNGTNGI0aM0PLly+2+XAD/1Jnh0d86Oy3oOanhIVrY3mcXLeizAzrHSh/ewcYm3bi8XPsPfhn0nK378Kz0+CF+RX2fHYDY1pnh0VZqeAyPRiQQdgBOQeM53IawA9AmhkfDTQg7AEExPBpuQNgBCIoaHtyAsAPQIadqeC0BSeM5nEbYAQiJE8OjW9f8aDyHkwg7AJYwPBqxiLADYJnTw6NDPTcQKsIOgGVOD49m0wrsRtgB6JTObFpZND4v6DlpPIdTCDsAnWa18Xxifi6N54gIwg5Al9F4jmhH2AHoMhrPEe0IOwC2YHg0ohlhB8A2VkPJ7k88p4aH9hB2AGxnpTmcxnOEA2EHwBF2N57nL/iTf/3Sm4ZRw4MlhB0AR9hdw2s2X/8KlOHRsIqwA+AYJ4ZHW6nhsWkFLQg7AI6zOjw6lMZzaniwgrADEBZWanirpl+qgZmnBz0fw6NhBWEHICysNp5vmHMljeewDWEHIGxoPEekEHYAwsrq8OhQa3gMj0YwhB2AiGB4NMKJsAMQEQyPRjgRdgAihhoewoWwAxBRTjWetx5FxvBoEHYAooLVxvOOhkePWfIajefwI+wARA07h0c3m1Mbz6nhxS/CDkDUsLJpRXKmhkfguRNhByCqOLlpRWJ4dLwi7ABEnVAbz1sH3qAsT9BzMjw6vhF2AKKWlRreyu9f0mHjOcOj4xdhByBq0XgOuxB2AKIajeewA2EHIOoxPBpdRdgBiBkMj0ZnEXYAYgY1PHQWYQcgpjhVw2sJSGp47kTYAYg5TgyPbh1iDI92H8IOQMyye3g0jefuRdgBiGl2Do+WGB7tVoQdgJjWmeHRi8bnBT0nw6Pdh7ADEPOsblqZmJ/L8Og4Q9gBcIXODI8OpfGcGp47EHYAXMVKDS+UxnOGR7sDYQfAVWg8R1sIOwCuw/BonIywA+BKDI9Ga4QdAFdjeDQkwg6Ay1HDg0TYAYgD0VDDI/Aii7ADEBecGB5N43nsIOwAxBWGR8cnwg5A3LE6PHrpjUODno/G8+hH2AGIO1aHR4++KItNKzGOsAMQl6Jh0wo1vPAh7ADELRrP4wdhByDu0XjufraH3fHjx3XfffepX79+Sk1NVf/+/fXAAw/IGONfY4zR/PnzlZWVpdTUVBUWFmr//v0B56mrq1NRUZE8Ho/S09N1yy236NChQ3ZfLgDQeB4HbA+7RYsWaenSpfqf//kf7d27V4sWLdLixYv1i1/8wr9m8eLFWrJkiZYtW6aKigqddtppGjVqlA4fPuxfU1RUpD179qi0tFTr16/Xli1bNGPGDLsvFwAkUcNzuwTT+i2XDa699lp5vV796le/8j83YcIEpaam6tlnn5UxRtnZ2Zo7d67uvPNOSVJDQ4O8Xq9WrFihSZMmae/evRo0aJC2b9+u4cOHS5I2btyoMWPG6NNPP1V2dnaH1+Hz+ZSWlqaGhgZ5PB47XyIAlyv/8G+64amKdo+vnpavgnP7SpKe2vxnPbjh/aDn2zj7Cg3I8lg+dzyz+2e47e/sLrvsMpWVlemDDz6QJL3zzjt6/fXXNXr0aEnSxx9/rOrqahUWFvq/Ji0tTfn5+SovL5cklZeXKz093R90klRYWKjExERVVLT9l6SpqUk+ny/gAQCdYaU5PJTG8zFLXqPxPMJsD7t77rlHkyZN0oABA9S9e3cNGTJEc+bMUVFRkSSpurpakuT1egO+zuv1+o9VV1crIyMj4HhSUpJ69+7tX3OykpISpaWl+R85OTl2vzQAccRq43mwX2k2m1Mbz6nhhZftYffb3/5Wq1at0nPPPae33npLK1eu1M9+9jOtXLnS7m8VYN68eWpoaPA/Dhw44Oj3A+Bu0VDDI/DsY3vY3XXXXf53d3l5eZo8ebJuv/12lZSUSJIyMzMlSTU1NQFfV1NT4z+WmZmp2tragOPHjh1TXV2df83JUlJS5PF4Ah4A0BWdGR69aHxe0HMyPDoybA+7r776SomJgaft1q2bmpubJUn9+vVTZmamysrK/Md9Pp8qKipUUFAgSSooKFB9fb0qKyv9azZt2qTm5mbl5wf/3TgA2M3K8OiJ+bkhNZ5Twwsv28Puuuuu00MPPaTf//73+stf/qK1a9fqkUce0fjx4yVJCQkJmjNnjh588EGtW7dOu3bt0pQpU5Sdna1x48ZJkgYOHKhrrrlG06dP17Zt2/TGG29o1qxZmjRpUkg7MQHAblZqeKE0njM8Orxsbz1obGzUfffdp7Vr16q2tlbZ2dm64YYbNH/+fCUnJ0s60VR+//33a/ny5aqvr9eIESP0xBNP6Pzzz/efp66uTrNmzdJLL72kxMRETZgwQUuWLFGvXsH/i6kFrQcA7HawsUkjH35VvsPH2l3Tus0glLaEllaDg41NuvXXlaqs+iKkc7ud3T/DbQ+7aEHYAXCC1VCyEngSfXgtor7PDgDcjOHRsYmwA4BOYHh0bCHsAKATGB4dWwg7AOikaGg851eaoSHsAKALOtN4HmrgSR03nlPDCw1hBwA2sNJ4HsrwaBrP7UXYAYBN7BweLTE82k6EHQDYxMqmFenEO7ylNw4Nek6GR9uDsAMAG1ndtDL6oixba3hsWmkbYQcANgu18bz1OzyGRzuLsAMAhzA8OnoQdgDgEBrPowdhBwAOovE8OhB2AOAwhkdHHmEHAGHixPDolvUMjw6OsAOAMHGihtfyjpAaXnCEHQCEETW8yCDsACDMGB4dfoQdAEQIw6PDh7ADgAiyOjx60fi8oOdjeHTbCDsAiCCrw6Mn5ufaXsOLh8Aj7AAgwpzctCIxPFoi7AAgKjA82lmEHQBEEYZHO4OwA4AowvBoZxB2ABBlaDy3H2EHAFGI4dH2IuwAIIo5MTw61BqemzatEHYAEMWo4dmDsAOAKEcNr+sIOwCIAZ0ZHr30xqFBzxlPw6MJOwCIIVaGR4++KIvh0f9E2AFAjLE6PLqjX2nGw/Bowg4AYozV4dFO1PBiLfAIOwCIQQyPtoawA4AYxfDo0BF2ABDjGB7dMcIOAGKcU43nLevd0HhO2AGACzhRw2upy7mh8ZywAwCXYHh0+wg7AHAZhkefirADAJdhePSpCDsAcKHO1PAWjc8Les5Ybjwn7ADApayG0sT8XNc2nhN2AOByVoZHT7uqvyuHRxN2ABAHnB4eHe2N54QdAMQBp4dHR/umFcIOAOJEPH/iOWEHAHEkXhvPCTsAiENONJ63rI/GxnPCDgDikBON561naUZbDY+wA4A4FU81PMIOAOKY1VCy+xPPw1XDI+wAAJaaw6dd1V8bZ18R9HzR1nhO2AEAJFlrDh+Q5bHceB7JGh5hBwCQFB01PKcCj7ADAPhZDSW7a3hObVoh7AAAp7A6PDqUxnMrNbwPqu19h0fYAQDaZKWGF0rjuZXh0dcvLe/EFbePsAMAtCnSn3huJ8IOANAupzattARkR/VBu1gOuy1btui6665Tdna2EhIS9OKLLwYcN8Zo/vz5ysrKUmpqqgoLC7V///6ANXV1dSoqKpLH41F6erpuueUWHTp0KGDNu+++qyuuuEI9evRQTk6OFi9ebP3VAQC6zInh0a3Xd1TDs4PlsPvyyy/1rW99S48//nibxxcvXqwlS5Zo2bJlqqio0GmnnaZRo0bp8OHD/jVFRUXas2ePSktLtX79em3ZskUzZszwH/f5fLr66qt1zjnnqLKyUj/96U/14x//WMuXL+/ESwQA2MGJ4dGh1vC6KsEYYzr9xQkJWrt2rcaNGyfpxLu67OxszZ07V3feeackqaGhQV6vVytWrNCkSZO0d+9eDRo0SNu3b9fw4cMlSRs3btSYMWP06aefKjs7W0uXLtW9996r6upqJScnS5Luuecevfjii3r//fdDujafz6e0tDQ1NDTI4/F09iUCAFo52NikkQ+/Kt/hY+2u2Tj7Cg3IOvFz96nNf9aDG4L/3F49LV8F5/bVwcYm3frrSlVWfaHmpq904NH/tO1nuK01u48//ljV1dUqLCz0P5eWlqb8/HyVl5/YWVNeXq709HR/0ElSYWGhEhMTVVFR4V9z5ZVX+oNOkkaNGqV9+/bpiy++sPOSAQAWhKvxPCkxwdbrtjXsqqurJUlerzfgea/X6z9WXV2tjIyMgONJSUnq3bt3wJq2ztH6e5ysqalJPp8v4AEAsF9nGs8Xjc8Les6TG8933n+1fRcsF+3GLCkpUVpamv+Rk5MT6UsCAFez0ng+MT+3w12XrRvP7WZr2GVmZkqSampqAp6vqanxH8vMzFRtbW3A8WPHjqmuri5gTVvnaP09TjZv3jw1NDT4HwcOHOj6CwIABGWl8bzg3L6Whkfbydaw69evnzIzM1VWVuZ/zufzqaKiQgUFBZKkgoIC1dfXq7Ky0r9m06ZNam5uVn5+vn/Nli1bdPToUf+a0tJSXXDBBTrjjLb/KyIlJUUejyfgAQBwlpXGcyn0Gt62P//d1uu0HHaHDh3Szp07tXPnTkknNqXs3LlTVVVVSkhI0Jw5c/Tggw9q3bp12rVrl6ZMmaLs7Gz/js2BAwfqmmuu0fTp07Vt2za98cYbmjVrliZNmqTs7GxJ0o033qjk5GTdcsst2rNnj37zm9/oscce0x133GHbCwcA2MOJTSvf/98dtl6j5daDV199Vf/6r/96yvNTp07VihUrZIzR/fffr+XLl6u+vl4jRozQE088ofPPP9+/tq6uTrNmzdJLL72kxMRETZgwQUuWLFGvXl83Ib777rsqLi7W9u3b1bdvX9122226++67Q75OWg8AIPy+s3SrdnzS/q751m0JVz+yWR/UHmpznd2tB13qs4tmhB0AhJ+VPryDjU26cXm59h/88pQ1Ud1nBwCIb6HU8K79xWs62NgU8vBoOzg3myXCWt6w0m8HAOGVIunF6UM1Z81O7fy0/pTjRyT98H/f0C8nnxgu8p9DztQ/vmzUT0u/nqPc3PSVpK9/lneVa3+N+dFHH6l///6RvgwAQBccOHBAZ599dpfP49p3dr1795YkVVVVKS0tLcJXE318Pp9ycnJ04MABappt4P50jHsUHPcnuI7ujzFGjY2N/l36XeXasEtMPFGOTEtL4y9aEPQkBsf96Rj3KDjuT3DB7o+db1TYoAIAcD3CDgDgeq4Nu5SUFN1///1KSUmJ9KVEJe5PcNyfjnGPguP+BBfu++Pa3ZgAALRw7Ts7AABaEHYAANcj7AAArkfYAQBcj7ADALgeYQcAcD3CDgDgeoQdAMD1/j9uFPRqfvM3aQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "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 = 32\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": [
      "785 µs ± 18.2 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)\n",
      "4.62 ms ± 7.62 µ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 = 512\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": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  2.6e-02, Efficiency:  6.2e+01 Gflops\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "n = 4000\n",
    "k = 50\n",
    "a = np.random.randn(n, n)\n",
    "v = np.random.randn(n, k)\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,  ((k*2 * n ** 2)/t) / 10 ** 9))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Time:  4.0e-04, Efficiency:  1.5e+00 Gflops\n"
     ]
    }
   ],
   "source": [
    "n = 4000\n",
    "k = 50\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, k)\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 * k) / 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 (Least recently used)\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": 18,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'Reshaped solution')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHFCAYAAADSY6wWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB5r0lEQVR4nO3deVxU5f4H8M+AMCAwg8iubC6JC26oOOaaJCpW3NQ0ydRQbwWWSyqm4lJdTa3cM1vE7tUs7WZlRpJrNwkTrYSEzDQ1G7QQxiVZz+8Pf5wY2WFmzjkzn/frNS9lzjMz33POLJ95znOeUQmCIICIiIiI6sVO6gKIiIiIlIghioiIiKgBGKKIiIiIGoAhioiIiKgBGKKIiIiIGoAhioiIiKgBGKKIiIiIGoAhioiIiKgBGKKIiIiIGoAhiiQxceJEuLq6Sl1Ggxw6dAgqlQqHDh2y2GMuXrwYKpXK5PcbHByMiRMnmvx+iarD13798LUvbwxRMpGcnAyVSgUnJyf89ttvlZYPHDgQnTp1atB9b9y4EcnJyY2skJTq6NGjWLx4MfLz86UuhWSk/D2n/NKkSRO0aNECEydOrPI9iJSHr33zayJ1AWSssLAQy5cvx7p160x2nxs3boSnpye/ddioo0ePYsmSJZg4cSLc3d2NluXk5MDOjt+lbNnSpUsREhKC27dv45tvvkFycjL+97//ITMzE05OTlKXR43A1775cQvKTNeuXfHmm2/i8uXLUpciKikpQVFRkdRlkBmo1Wo4ODhIXQZJaNiwYXjssccwefJkvPXWW3juuedw9uxZfPLJJ1KXRmbE175pMETJzPPPP4/S0lIsX7681rYlJSV44YUX0Lp1a6jVagQHB+P5559HYWGh2CY4OBhZWVk4fPiw2G0/cODAau/z/PnzUKlUWLVqFVavXi3e948//ggAyM7OxqhRo+Dh4QEnJyf06NGj0pttcXExlixZgrZt28LJyQnNmzdH3759kZqaWunxfvvtN8TExMDV1RVeXl547rnnUFpaatRm1apV6NOnD5o3bw5nZ2eEh4dj165dle5LpVIhISEB27ZtQ7t27eDk5ITw8HAcOXKkysd94okn4OPjA7VajY4dO+Kdd96p1O7SpUuIiYmBi4sLvL29MWPGDKPtW5Pr169j+vTpCA4Ohlqthre3N+6//36cOHHCqN3OnTsRHh4OZ2dneHp64rHHHqv1cEr5fqrqMK1KpcLixYsB3BlPMXv2bABASEiI+Bw4f/48gKrHRfzyyy8YPXo0PDw80LRpU/Tu3RufffaZUZvysSEffPABXnrpJbRs2RJOTk4YPHgwfv755zptH5Knfv36AQDOnj1rdD1f+3ztA3zt342H82QmJCQEjz/+ON58800kJibC39+/2raTJ0/G1q1bMWrUKMyaNQvp6elYtmwZTp8+jY8++ggAsHr1akybNg2urq6YP38+AMDHx6fWOrZs2YLbt29j6tSpUKvV8PDwQFZWFu699160aNECiYmJcHFxwQcffICYmBh8+OGH+Mc//gHgzot32bJlmDx5Mnr16gWDwYDjx4/jxIkTuP/++8XHKC0tRVRUFCIiIrBq1Sp8+eWXeOWVV9C6dWs89dRTYrs1a9bgwQcfRGxsLIqKirBjxw6MHj0ae/bsQXR0tFHdhw8fxvvvv49nnnkGarUaGzduxNChQ3Hs2DFxTFlubi569+4tvvF6eXnh888/R1xcHAwGA6ZPnw4A+OuvvzB48GBcuHABzzzzDPz9/fHvf/8bBw4cqMOeBJ588kns2rULCQkJ6NChA/7880/873//w+nTp9G9e3cAd8alTJo0CT179sSyZcuQm5uLNWvW4Ouvv8bJkycrdcHX18MPP4yffvoJ7733Hl577TV4enoCALy8vKpsn5ubiz59+uDWrVt45pln0Lx5c2zduhUPPvggdu3aJe7jcsuXL4ednR2ee+45FBQUYMWKFYiNjUV6enqj6ibplH/INmvWTLyOr32+9vnar4ZAsrBlyxYBgPDtt98KZ8+eFZo0aSI888wz4vIBAwYIHTt2FP/+7rvvBADC5MmTje7nueeeEwAIBw4cEK/r2LGjMGDAgDrVce7cOQGAoNFohCtXrhgtGzx4sBAWFibcvn1bvK6srEzo06eP0LZtW/G6Ll26CNHR0TU+zoQJEwQAwtKlS42u79atmxAeHm503a1bt4z+LioqEjp16iTcd999RtcDEAAIx48fF6/79ddfBScnJ+Ef//iHeF1cXJzg5+cn/PHHH0a3Hzt2rKDVasXHW716tQBA+OCDD8Q2N2/eFNq0aSMAEA4ePFjjOmq1WiE+Pr7a5UVFRYK3t7fQqVMn4a+//hKv37NnjwBASEpKEq9btGiRUPHlWr6ftmzZUul+AQiLFi0S/165cqUAQDh37lyltkFBQcKECRPEv6dPny4AEL766ivxuuvXrwshISFCcHCwUFpaKgiCIBw8eFAAILRv314oLCwU265Zs0YAIJw6dara9SZ5KH/P+fLLL4WrV68KFy9eFHbt2iV4eXkJarVauHjxotiWr32+9vnarxoP58lQq1atMH78eGzevBm///57lW327t0LAJg5c6bR9bNmzQKASl2w9TVy5Eijbyx5eXk4cOAAHnnkEVy/fh1//PEH/vjjD/z555+IiorCmTNnxG5od3d3ZGVl4cyZM7U+zpNPPmn0d79+/fDLL78YXefs7Cz+/9q1aygoKEC/fv0qdY0DgE6nQ3h4uPh3YGAgHnroIXzxxRcoLS2FIAj48MMP8cADD0AQBHE9/vjjD0RFRaGgoEC8371798LPzw+jRo0S769p06aYOnVqretVvh3S09OrHd92/PhxXLlyBU8//bTRAN7o6GiEhoY2eh82xN69e9GrVy/07dtXvM7V1RVTp07F+fPnxcO65SZNmgRHR0fx7/JDQXfvQ5KvyMhIeHl5ISAgAKNGjYKLiws++eQTtGzZEgBf++X42udrvyoMUTK1YMEClJSUVDs26tdff4WdnR3atGljdL2vry/c3d3x66+/NurxQ0JCjP7++eefIQgCFi5cCC8vL6PLokWLAABXrlwBcOdsn/z8fNxzzz0ICwvD7Nmz8cMPP1R6DCcnp0pdy82aNcO1a9eMrtuzZw969+4NJycneHh4wMvLC6+//joKCgoq3Wfbtm0rXXfPPffg1q1buHr1Kq5evYr8/Hxs3ry50npMmjTJaD1+/fVXtGnTptIcLe3atatx25VbsWIFMjMzERAQgF69emHx4sVGbzDl+6iq+wsNDW30PmyIX3/9tcp62rdvLy6vKDAw0Ojv8kNAd+9Dkq8NGzYgNTUVu3btwvDhw/HHH39ArVaLy/na/xtf+3/ja/8OjomSqVatWuGxxx7D5s2bkZiYWG07c0zCBhh/AwSAsrIyAMBzzz2HqKioKm9THuj69++Ps2fP4uOPP8a+ffvw1ltv4bXXXsOmTZswefJksb29vX2tdXz11Vd48MEH0b9/f2zcuBF+fn5wcHDAli1bsH379nqvV/l6PPbYY5gwYUKVbTp37lzv+63KI488gn79+uGjjz7Cvn37sHLlSrz88sv473//i2HDhjXqvqvb73cPzDW36vahIAgWrYMarlevXujRowcAICYmBn379sW4ceOQk5MDV1dXvvYbgK9928EQJWMLFizAf/7zH7z88suVlgUFBaGsrAxnzpwRvykAdwYH5ufnIygoSLzOFEGrVatWAAAHBwdERkbW2t7DwwOTJk3CpEmTcOPGDfTv3x+LFy82eiOtiw8//BBOTk744osvjL4db9mypcr2VR1G+Omnn9C0aVPxm6+bmxtKS0trXY+goCBkZmZCEASjbZiTk1Pn+v38/PD000/j6aefxpUrV9C9e3e89NJLGDZsmLiPcnJycN999xndLicnx2gf3q38W9/dk+hV9Q22Pvs/KCioyvXLzs4Wl5P1sre3x7JlyzBo0CCsX78eiYmJfO3ztS8up8p4OE/GWrdujcceewxvvPEG9Hq90bLhw4cDuHP2XUWvvvoqABidueLi4tLoGWu9vb0xcOBAvPHGG1WO07p69ar4/z///NNomaurK9q0aVPn04Mrsre3h0qlMvqWdf78eezevbvK9mlpaUbjJS5evIiPP/4YQ4YMgb29Pezt7TFy5Eh8+OGHyMzMrHE9hg8fjsuXLxudUn3r1i1s3ry51rpLS0srHXLw9vaGv7+/uB169OgBb29vbNq0yWjbfP755zh9+nSls48q0mg08PT0rHQK98aNGyu1dXFxAVD5Tbcqw4cPx7Fjx5CWliZed/PmTWzevBnBwcHo0KFDrfdByjZw4ED06tULq1evxu3bt/na/3987fO1XxX2RMnc/Pnz8e9//xs5OTno2LGjeH2XLl0wYcIEbN68Gfn5+RgwYACOHTuGrVu3IiYmBoMGDRLbhoeH4/XXX8eLL76INm3awNvbu9K3n7rYsGED+vbti7CwMEyZMgWtWrVCbm4u0tLScOnSJXz//fcAgA4dOmDgwIEIDw+Hh4cHjh8/Lp7uW1/R0dF49dVXMXToUIwbNw5XrlzBhg0b0KZNmyrHWnTq1AlRUVFGpzkDwJIlS8Q2y5cvx8GDBxEREYEpU6agQ4cOyMvLw4kTJ/Dll18iLy8PADBlyhSsX78ejz/+ODIyMuDn54d///vfaNq0aa11X79+HS1btsSoUaPQpUsXuLq64ssvv8S3336LV155BcCdb/Yvv/wyJk2ahAEDBuDRRx8VT3MODg7GjBkzanyMyZMnY/ny5Zg8eTJ69OiBI0eO4KeffqrUrnyw7fz58zF27Fg4ODjggQceEN9gK0pMTMR7772HYcOG4ZlnnoGHhwe2bt2Kc+fO4cMPP+QMxzZi9uzZGD16NJKTk/Hkk0/ytc/XPl/71ZHqtEAyVnGKg7uVnxJccYoDQRCE4uJiYcmSJUJISIjg4OAgBAQECPPmzTM6DVkQBEGv1wvR0dGCm5ubAKDG6Q7KT59duXJllcvPnj0rPP7444Kvr6/g4OAgtGjRQhgxYoSwa9cusc2LL74o9OrVS3B3dxecnZ2F0NBQ4aWXXhKKioqM1snFxaXS/d99Oq8gCMLbb78ttG3bVlCr1UJoaKiwZcuWKtsBEOLj44X//Oc/Yvtu3bpVeTpybm6uEB8fLwQEBAgODg6Cr6+vMHjwYGHz5s1G7X799VfhwQcfFJo2bSp4enoKzz77rJCSklLrac6FhYXC7NmzhS5dughubm6Ci4uL0KVLF2Hjxo2V2r7//vtCt27dBLVaLXh4eAixsbHCpUuXat0ut27dEuLi4gStViu4ubkJjzzyiHDlypVKpzkLgiC88MILQosWLQQ7OzujU57vPs1ZEO7s41GjRgnu7u6Ck5OT0KtXL2HPnj1GbcpPc965c6fR9TWdfk3yUtN7TmlpqdC6dWuhdevWQklJiSAIfO3ztX8HX/vGVIJgY6PAyGqpVCrEx8dj/fr1UpdCRBbE1z5Jhf1zRERERA3AEEVERETUAAxRRERERA3AMVFEREREDcCeKCIiIqIGYIgiIiIiagBOtllHZWVluHz5Mtzc3Mz2e3VEVDNBEHD9+nX4+/srZvI/vncQScuc7xsMUXV0+fJlBAQESF0GEeHOT3q0bNlS6jLqhO8dRPJgjvcNhqg6cnNzA3BnJ2g0GomrIbJNBoMBAQEB4utRCcprPZTuBVdXZfSeEVmTGzfKMDDiqlneNxii6qi8G16j0TBEEUlMSYfFymt1dbWDqxtDFJFUzPG+wVc0ERERUQMwRBGR4nz99dd44IEH4O/vD5VKhd27dxstFwQBSUlJ8PPzg7OzMyIjI3HmzBmjNnl5eYiNjYVGo4G7uzvi4uJw48YNozY//PAD+vXrBycnJwQEBGDFihXmXjUiUhCGKCJSnFu3bqFLly7YsGFDlctXrFiBtWvXYtOmTUhPT4eLiwuioqJw+/ZtsU1sbCyysrKQmpqKPXv24MiRI5g6daq43GAwYMiQIQgKCkJGRgZWrlyJxYsXY/PmzWZfPyJSBs5YXkcGgwFarRYFBQUcE0UkkapehyqVCh999BFiYmIA3OmF8vf3x6xZs/Dcc88BAAoKCuDj44Pk5GSMHTsWp0+fRocOHfDtt9+iR48eAICUlBQMHz4cly5dgr+/P15//XXMnz8fer0ejo6OAIDExETs3r0b2dnZ9a75eJYPx0QRSeDG9TL06Jhrls9vvqKJyKqcO3cOer0ekZGR4nVarRYRERFIS0sDAKSlpcHd3V0MUAAQGRkJOzs7pKeni2369+8vBigAiIqKQk5ODq5du1bt4xcWFsJgMBhdiMg6MUQRkVXR6/UAAB8fH6PrfXx8xGV6vR7e3t5Gy5s0aQIPDw+jNlXdR8XHqMqyZcug1WrFC+eIIrJeDFFERCY0b948FBQUiJeLFy9KXRIRmQlDFBFZFV9fXwBAbm6u0fW5ubniMl9fX1y5csVoeUlJCfLy8ozaVHUfFR+jKmq1WpxPjvPKEVk3TrZJRFYlJCQEvr6+2L9/P7p27QrgzuDu9PR0PPXUUwAAnU6H/Px8ZGRkIDw8HABw4MABlJWVISIiQmwzf/58FBcXw8HBAQCQmpqKdu3aoVmzZpZfMRuxy9Bd6hIUY5TmhNQl2DyGKCJSnBs3buCXX34R/z537hy+++47eHh4IDAwENOnT8eLL76Itm3bIiQkBAsXLoS/v794Bl/79u0xdOhQTJkyBZs2bUJxcTESEhIwduxY+Pv7AwDGjRuHJUuWIC4uDnPnzkVmZibWrFmD1157TYpVVhyGIfNr6DZm+DIdmwtRGzZswMqVK6HX69GlSxesW7cOvXr1krosIqqHkydPYsSIEeLfM2fOBABMmDABycnJmDNnDm7evImpU6ciPz8fffv2RUpKCpycnMTbbNu2DQkJCRg8eDDs7OwwcuRIrF27Vlyu1Wqxb98+xMfHIzw8HJ6enkhKSjKaS8pWMSApW132H4NW3djUPFHvv/8+Hn/8cWzatAkRERFYvXo1du7ciZycnEpn6tyN80QRSU+Jr0MlzxPFsER3U2K4Muc8UTYVoiIiItCzZ0+sX78eAFBWVoaAgABMmzYNiYmJRm0LCwtRWFgo/l3+6/G17YQ1X57Ba1/+VKd6nBzs0LmlO9r7uiHE0wUtmjWFt5saHi6OcFU3gbOjPRzs7WCnUtYPrpL1EwQBZQJQVFKGv4pLcf12MfJuFiHXUIjf8v/CL1dv4PTvBnx/qQClZXV7i1k1ugtGhbessQ1DlPkwMFFDyT1YmTNE2czhvKKiImRkZGDevHnidXZ2doiMjBQn4Kto2bJlWLJkSb0fp64BCgBuF5fh2Lk8HDuXV+/HsQUqFeBgbwcHOxUcmtihiZ0dHOxVaGKvgr1KhSb/HzCb2NlBpQKa2NtB9f+3a2JXv9BZUnYnFEAQxP+XlpWh9P//X1JWhpJSAcWlAkrKylBcUobiUgFFpWVmWXdb9NzO72sNUWRaDE5kChWfR3IPVKZmMyHqjz/+QGlpaZWT51X1Ew7z5s0Tx1kAf/dEkeUI/9/TUQQARaVSl0OkeAxNZG53P8esPVTZTIiqL7VaDbVaXe/btfV2xZkrN2pvSERGZke1k7oEq8XwRFIpf+5Za5iS7wF6E/P09IS9vX2NE/CZwpT+rUx2X0S25IHO/lKXYHV2GbozQJEsWOtz0WZClKOjI8LDw7F//37xurKyMuzfvx86nc5kj3NvG0+T3ReRLWnZzFnqEqyGtX5gkfJZ23PTpg7nzZw5ExMmTECPHj3Qq1cvrF69Gjdv3sSkSZNM9hgt3PlBQNQQdvU8GYAqs6YPJ7Ju1nKYz6ZC1JgxY3D16lUkJSVBr9eja9euSElJqTTYnIhIaRigSIl2GborOkjZVIgCgISEBCQkJEhdBhGRydhygErVh0pdgsnc71v5THFboOQgZXMhiojImlhrgLKmcFRXdV1nawxbSg1SDFFEJLlhnUx3hqwtUXqAssWgZAo1bTclBywlBimGKDOIjQjEtvQLUpdBpBjjIgKlLkFxlBSgGJYsp7ptrZRwpbQgxRBlBuN1QQxRRPXQl1OD1IvcAxRDk/zcvU/kHKqUFKQYoswguLmL1CUQKQp/YLvu5BigGJqUR+6hSilBiiHKDJwc7KUugYjIrBicrEvF/Sm3QCVnDFFERAohdS8Ug5NtkEugUkJvlM387AsRkZJJGaBS9aEMUDZK6n0v9ReH2rAniogkFdS8qdQlUDUYnKhc+XOBh/qMsSfKTHoENZO6BCJFmHZfW6lLkD1LfxuXuveB5EuK54ace6MYosxkzlC+ARHVRXSYn9Ql0P9jeKK64nPlDh7OM5O23q5Sl0CkCM6OPJtVavwwpIay9cN8DFFm0szFUeoSiMgKmPtQhi0FqPOXvCz+mMEtr1r8MaWQqg81a5CS65l6DFFERDbKWgKUFOGorupamzWELXMHKTliiCIisjFKDE9yDkqmUNP6KSlg2drhPYYoIiKZMsehPCUEKGsPTPVV1faQe7AyR6+UHA/pMUSZUSsvF/xy9abUZRDJ1pIHO0pdgk2Ra4BiaKq/u7eZHEOVLRze4xQHZvTsYM5/Q1STqI6+UpdgM+QWoM5f8hIv1Hhy3Z5ye96ZGnuizEjXqrnUJRDJmq/WSeoSZEvOEww2lNw+4K1Vxe0sxx6qxpDbIT2GKDPy1vADgoikJ2VvAIOTtOQQqKz5sB4P5xERWTGpApQcDy3ZOin3ibUe1mOIIiKyUlJ8cDE8yZ9U+8gagxQP5xERyYwSx0MxOClP+T6ztnFTlsSeKDPz58BZoio9P9z6vpXKiaW+9bPnSfksuQ9N8byU05cMhigzm8ZpDoiqxOkNzMcSAYrhyfpYap9a02E9higzG9TOW+oSiGQpqLmL1CVQAzE8WTfu37pjiDIzzoNDRJZk7m/5/IC1Debez9bSG8WB5UREVCtrDk/qC46Nvo/CwCITVCIvHHheO4YoIiIZacygWXN9u1dygDJFQDLV4yg1aJ2/5GWWINWYSTjlMnM5QxQREVVLKQHKUmGpMaqrUQnhylxBSukYoiyghbszfsv/S+oyiGQjcZh1jIewdnIOUEoITXV197rINVQxSFXGgeUW8MzgNlKXQCQrwzpxegNTM/WhPLkFKPUFR6OLNZPzupr6eaH0AeZmC1EvvfQS+vTpg6ZNm8Ld3b3KNhcuXEB0dDSaNm0Kb29vzJ49GyUlJUZtDh06hO7du0OtVqNNmzZITk6udD8bNmxAcHAwnJycEBERgWPHjhktv337NuLj49G8eXO4urpi5MiRyM3NNdWq1mogpzkgMsLpDeRNLgFKrkHC0uS2HeTy/JADs4WooqIijB49Gk899VSVy0tLSxEdHY2ioiIcPXoUW7duRXJyMpKSksQ2586dQ3R0NAYNGoTvvvsO06dPx+TJk/HFF1+Ibd5//33MnDkTixYtwokTJ9ClSxdERUXhypUrYpsZM2bg008/xc6dO3H48GFcvnwZDz/8sLlWvRIfDac5ICLzMeW3eak/IOUWGORGLtvHlM8TJfdGqQRBEMz5AMnJyZg+fTry8/ONrv/8888xYsQIXL58GT4+PgCATZs2Ye7cubh69SocHR0xd+5cfPbZZ8jMzBRvN3bsWOTn5yMlJQUAEBERgZ49e2L9+vUAgLKyMgQEBGDatGlITExEQUEBvLy8sH37dowaNQoAkJ2djfbt2yMtLQ29e/eu03oYDAZotVoUFBRAo9HUezsEJ35W79sQWavzy6MbdLvGvg6lUF7z8SwfuLrV/r21IWfnmepDSMoAJXUoUDIpx1CZaoxUQ87Sq+vZeTeul6FHx1yzvG9INiYqLS0NYWFhYoACgKioKBgMBmRlZYltIiMjjW4XFRWFtLQ0AHd6uzIyMoza2NnZITIyUmyTkZGB4uJiozahoaEIDAwU21SlsLAQBoPB6EJERKYlh14VpeM2lI5kIUqv1xsFKADi33q9vsY2BoMBf/31F/744w+UlpZW2abifTg6OlYal1WxTVWWLVsGrVYrXgICAhq0nkREdWVLvVD84Dc9KbapqZ43DXkey+GHiOsVohITE6FSqWq8ZGc3bOIsuZk3bx4KCgrEy8WLFxt1f228XU1UGZGyLXmwo9Ql0F0sGaAYnszP0ttY6nF0UqrXPFGzZs3CxIkTa2zTqlWrOt2Xr69vpbPoys+Y8/X1Ff+9+yy63NxcaDQaODs7w97eHvb29lW2qXgfRUVFyM/PN+qNqtimKmq1Gmq1uk7rUhez7r8HT22TfnZVIqkNC+P0BraIwcnyyre5XOedsgb16ony8vJCaGhojRdHx7q9UHQ6HU6dOmV0Fl1qaio0Gg06dOggttm/f7/R7VJTU6HT6QAAjo6OCA8PN2pTVlaG/fv3i23Cw8Ph4OBg1CYnJwcXLlwQ21hCrxAPiz0WkZx5u/FsVTmxRC8CA5S0LLH9bbU3ymwzll+4cAF5eXm4cOECSktL8d133wEA2rRpA1dXVwwZMgQdOnTA+PHjsWLFCuj1eixYsADx8fFiD9CTTz6J9evXY86cOXjiiSdw4MABfPDBB/jss7/PdJs5cyYmTJiAHj16oFevXli9ejVu3ryJSZMmAQC0Wi3i4uIwc+ZMeHh4QKPRYNq0adDpdHU+M88UmruarleLiAho/Hgoc3/wMTzJhyV6pRo7o3ljfktPKmYLUUlJSdi6dav4d7du3QAABw8exMCBA2Fvb489e/bgqaeegk6ng4uLCyZMmIClS5eKtwkJCcFnn32GGTNmYM2aNWjZsiXeeustREVFiW3GjBmDq1evIikpCXq9Hl27dkVKSorRYPPXXnsNdnZ2GDlyJAoLCxEVFYWNGzeaa9WJiGweA5Q8qS848vCeCZl9nihrYYr5aThXFFHD54gCrH+eqPqebdSYnihz9kIxQMmfOYNUY3qj6tsTVZe5oqxynigisj3ebjysbSpyneWZAUoZuJ9MgyHKgnoFc3A52bY5Q+X5wW9rzNULxQ9mZTHX/mrM80uuXw6qY7YxUVRZ4vBQPLzxqNRlEElmRGc/qUsgM7GGAOX2a/1Ht1wPUpmhEsvhGKnGYYiyoHt83KQugUhSTg72UpcgW5aafdkcvVBKCVANCUmNvU8lhCxzBKnGnqlXV7sM3ev8G3rmwBBlQa5qbm4isi5yDlDmCE2NrUGuoYo9Ug3DMVFERArT0HEjpu6FkluAcvtVMLrIkZxrNPX+bOjzTUnjotg1QkRE9SaXACW3IFJfFeuXQy8Ve6TqhyGKiCxiYDvb/FkIayR1gFJ6cKpO+XpJHaYYpOqOh/MsbGKfYKlLIJLE1P51+3FyMg9THcqTMkDJ8RCYOcjhcJ+p9rO1/6YeQ5SFjdcFSV0CkSR0rZpb9PFKS0uxcOFChISEwNnZGa1bt8YLL7yAij/SIAgCkpKS4OfnB2dnZ0RGRuLMmTNG95OXl4fY2FhoNBq4u7sjLi4ON27csOi62DqpA4WUbHndlYAhysJCmrtIXQKRJFQqyx6iePnll/H6669j/fr1OH36NF5++WWsWLEC69atE9usWLECa9euxaZNm5Ceng4XFxdERUXh9u3bYpvY2FhkZWUhNTUVe/bswZEjRzB16lSLrotcWLoXigHib1JsC6kP2yoBQ5SF2dlJP3CQyBYcPXoUDz30EKKjoxEcHIxRo0ZhyJAhOHbsGIA7vVCrV6/GggUL8NBDD6Fz58549913cfnyZezevRsAcPr0aaSkpOCtt95CREQE+vbti3Xr1mHHjh24fPmyJOvVkDOXTHFIRYoARZUpMUg15PmnlDP0GKKIyCr16dMH+/fvx08//QQA+P777/G///0Pw4YNAwCcO3cOer0ekZGR4m20Wi0iIiKQlpYGAEhLS4O7uzt69OghtomMjISdnR3S09OrfNzCwkIYDAajC9Ude59qx20kHzw7j4isUmJiIgwGA0JDQ2Fvb4/S0lK89NJLiI2NBQDo9XoAgI+Pj9HtfHx8xGV6vR7e3t5Gy5s0aQIPDw+xzd2WLVuGJUuWmHp1JGWpXigGg/px+1WwyJl8PFuveuyJkkDLZs5Sl0BkUfOGWb5r/oMPPsC2bduwfft2nDhxAlu3bsWqVauwdetWsz7uvHnzUFBQIF4uXrxo1serjRLOjmLPSsMpZdsp4XnYEAxREpjLX7InG/NAF3+LP+bs2bORmJiIsWPHIiwsDOPHj8eMGTOwbNkyAICvry8AIDc31+h2ubm54jJfX19cuXLFaHlJSQny8vLENndTq9XQaDRGFyUzdy+UEgKAEph7O3KQedUYoiTQt42n1CUQWZS/u+V7X2/dugU7O+O3OHt7e5SVlQEAQkJC4Ovri/3794vLDQYD0tPTodPpAAA6nQ75+fnIyMgQ2xw4cABlZWWIiIiwwFpIiwFKWbg9LY9joiTQzIWJnsjcHnjgAbz00ksIDAxEx44dcfLkSbz66qt44oknANyZcmH69Ol48cUX0bZtW4SEhGDhwoXw9/dHTEwMAKB9+/YYOnQopkyZgk2bNqG4uBgJCQkYO3Ys/P0t37umlDOW6oIf+OZhznFSlh4blaoPxf2+2RZ7vIZgiCIiq7Ru3TosXLgQTz/9NK5cuQJ/f3/885//RFJSkthmzpw5uHnzJqZOnYr8/Hz07dsXKSkpcHJyEtts27YNCQkJGDx4MOzs7DBy5EisXbvWpLXuMnQ36f2Va8w4FHP2QjFAmZelBpzX1/lLXghueVXqMkxKJVScvpeqZTAYoNVqUVBQYJIxDsGJn5mgKiJlOL882iT3Y+rXoSWU13w8yweublWPoKhriKpvT5QcQxQDlOWYK0g1pjeqviGqrj1RozQnql1243oZenTMNcv7BsdEEZFZdQ90l7oEkgkGKMvi9jY/hiiJjO0ZIHUJRBbxzOC2UpdA9WSOXih+oEvDHNudZ+r9jSFKIpP7hUhdApFFDLjHOueHkTs5zcvDACUtOW1/OT0vTYEhSiKtPF2lLoHIIiz9w8PUOKbuZZDTB7gtM/V+YG/UHQxREuEPERORtWOAkhfuD9PjFAdERGSztGcLTXZfBa3VJrsvUgb2RBERKUB9pjdo6LgTUx6ikWuvh/ZsodFFKfdtKqbcLw19vtTn+Sn3CWbZEyWhdj5uyMm9LnUZRGaz9KGOUpdAEpBbgJIq0FR8XDn1Usl1Mk4lYk+UhBaO6CB1CURm9VCXFlKXQDZKbj1CcquHTIMhSkJdArRSl0BkVtqmDlKXQHVkqkN5UvdCKSGoyKFGU+0nWz9Lj4fzJOTmxA8YIrIeUgYoqUNJQ5TXLNWhPh7WazyGKCIiUiwlhqe7SR2mqOHMdjjv/PnziIuLQ0hICJydndG6dWssWrQIRUXGP1z4ww8/oF+/fnByckJAQABWrFhR6b527tyJ0NBQODk5ISwsDHv37jVaLggCkpKS4OfnB2dnZ0RGRuLMmTNGbfLy8hAbGwuNRgN3d3fExcXhxo0bpl9xIiIJNeTMPFMckpGiF8oaAlRFUqyPKfZbQ54/1jJzudlCVHZ2NsrKyvDGG28gKysLr732GjZt2oTnn39ebGMwGDBkyBAEBQUhIyMDK1euxOLFi7F582axzdGjR/Hoo48iLi4OJ0+eRExMDGJiYpCZmSm2WbFiBdauXYtNmzYhPT0dLi4uiIqKwu3bt8U2sbGxyMrKQmpqKvbs2YMjR45g6tSp5lp9Ips3sU+w1CWQlZLDmCJzseZ1s0YqQRAs9vVh5cqVeP311/HLL78AAF5//XXMnz8fer0ejo53kmxiYiJ2796N7OxsAMCYMWNw8+ZN7NmzR7yf3r17o2vXrti0aRMEQYC/vz9mzZqF5557DgBQUFAAHx8fJCcnY+zYsTh9+jQ6dOiAb7/9Fj169AAApKSkYPjw4bh06RL8/f1rrd1gMECr1aKgoAAajcZk2+Tfaeex8OMsk90fkVzsnzUArb1M+/NG5nodmlN5zcezfODqVvX31l2G7rXeT13ny5GiJ8qSvVC2FDAseXivsWOjCgOLam90l+CWV+vU7n7f7FrbjNKcqHbZjetl6NEx1yzvGxY9O6+goAAeHh7i32lpaejfv78YoAAgKioKOTk5uHbtmtgmMjLS6H6ioqKQlpYGADh37hz0er1RG61Wi4iICLFNWloa3N3dxQAFAJGRkbCzs0N6enqVtRYWFsJgMBhdzGFYmJ9Z7pdIaqYOUGQeSjq7ypYCFKCs9VXS88iULBaifv75Z6xbtw7//Oc/xev0ej18fHyM2pX/rdfra2xTcXnF21XXxtvb22h5kyZN4OHhIba527Jly6DVasVLQEBAvda3rjxdOZCQiJTLUr1QSgoUpmSp9ZZ6agqlqneISkxMhEqlqvFSfiiu3G+//YahQ4di9OjRmDJlismKN6d58+ahoKBAvFy8eFHqkojIRsn9py/MzVYDVDlbX385P//rPcXBrFmzMHHixBrbtGrVSvz/5cuXMWjQIPTp08dowDgA+Pr6Ijc31+i68r99fX1rbFNxefl1fn5+Rm26du0qtrly5YrRfZSUlCAvL0+8/d3UajXUavYSEZFyWPqMJ0v0Xth6gCinPVto9jFSlp436vwlrzqPi5KrevdEeXl5ITQ0tMZL+Rin3377DQMHDkR4eDi2bNkCOzvjh9PpdDhy5AiKi4vF61JTU9GuXTs0a9ZMbLN//36j26WmpkKn0wEAQkJC4Ovra9TGYDAgPT1dbKPT6ZCfn4+MjAyxzYEDB1BWVoaIiIj6bgIiqsWgdtZx+rK1k/s4FgYoY3LfHnJ/PpmD2cZElQeowMBArFq1ClevXoVerzcagzRu3Dg4OjoiLi4OWVlZeP/997FmzRrMnDlTbPPss88iJSUFr7zyCrKzs7F48WIcP34cCQkJAACVSoXp06fjxRdfxCeffIJTp07h8ccfh7+/P2JiYgAA7du3x9ChQzFlyhQcO3YMX3/9NRISEjB27Ng6nZlnbk8PbC11CUQmNW1wW6lLIIWTe2CQCreLvJhtxvLU1FT8/PPP+Pnnn9GyZUujZeWzKmi1Wuzbtw/x8fEIDw+Hp6cnkpKSjOZv6tOnD7Zv344FCxbg+eefR9u2bbF792506tRJbDNnzhzcvHkTU6dORX5+Pvr27YuUlBQ4OTmJbbZt24aEhAQMHjwYdnZ2GDlyJNauXWuu1a+X2N5B2HjorNRlEJlM98BmUpdAZmbOQ3kMCjUz56E9/hRM/Vh0niglM/f8NMGJn5n8Pomkcn55tFnu11bniTLXHFGNOfxirhDFAFV35gpSjQlR9Z0vylRzRdnEPFFERKR8DFDyYK7txekO6o4hiojIBsltEDADVMPIbbvJ7XllbgxRRGRSXQPcpS5BUeryky9EJE8MUTIxQRckdQlEJjEnqp3UJdgkS80RZY5DPXLrTVEac2w/Sx3SM9XzVqovIwxRMjG5X6vaGxEpgK51c6lLIAVhgDINbkdpMETJRMtmzlKXQGQSKhVPjyYi28AQJRP84CEiS2no4F9TH+Jh74lpmXp7NnR/29LgcoYoIiIiogZgiCIikwlroZW6BFII9kKZB7erZTFEycjEPsFSl0DUKPOG1W1mbaq7us5WTkSWxxAlI5P7hUhdAlGj8Mw862XK8VDsLTEvU25fucxeLtcvEwxRMtKyWVOpSyBqFJ4gQUS2hCGKiMiGSH3mFHuhLEPq7Sz188xSmkhdABERNY6lZiu3JY7Zl6pdVhTa0oKVWL/zl7wQ3PKq1GU0CEMUEZlEv7aeUpdAZmKqcTFS947UpKbQVFtbuYYq7dlCFLRWN/p+3H4VcD2Ih+qrwhAlM8/c1wZrD/wsdRlE9TZrCH8zj5SlPsGprvcj10BF5sExUTLzWG/+EDEpU9cAd6lLIBmTUy+UY/YlkwUoS953Q8hpu1sj9kTJjLfGSeoSiIiskiXDTfljsWfKurEnioiIqiWXeYIaS6reITn1SjWGtTwPTI0hiojIRkh12rmUh5TkcHhN6hqk2v62MM0BQxQRNdrYngFSl0BUidTh6W5yq4cajyFKhl4eGSZ1CUT18uSA1lKXQGREroFFrnVRwzBEydDQjn5Sl0BUL8GeLlKXQDIlxaEkuQcVKerjWXrmwRAlQ9qmDlKXQESkyMHEcg9Q5ZRSZ0VKfD6YG0MUERFZBaUFE6XVS5UxRBERkVlY8hCSUgOJJevmIT3TY4iSqdZeHGNCyrBwRAepSyAbp9QAVU7p9dsyhiiZWvZwZ6lLIKqTMZzeQFLnL3lJXQJRoyn1ecwQJVNdArRSl0BUJ65q/nqUElh64kNLHTqyll4cS60HD+mZFkOUTKmb2EtdAhHZMCWciWUtAaqcta0PYP2zljNEERERUZ0oIVxbEkMUETWYn9ZJ6hLIRlljrw1gvetlrRiiZOzRXhywS/K2fCRPgKDKOO5G3rh/TMesIerBBx9EYGAgnJyc4Ofnh/Hjx+Py5ctGbX744Qf069cPTk5OCAgIwIoVKyrdz86dOxEaGgonJyeEhYVh7969RssFQUBSUhL8/Pzg7OyMyMhInDlzxqhNXl4eYmNjodFo4O7ujri4ONy4ccP0K21C0+5rK3UJRDXq39ZT6hLIBll7b421r581MWuIGjRoED744APk5OTgww8/xNmzZzFq1ChxucFgwJAhQxAUFISMjAysXLkSixcvxubNm8U2R48exaOPPoq4uDicPHkSMTExiImJQWZmpthmxYoVWLt2LTZt2oT09HS4uLggKioKt2/fFtvExsYiKysLqamp2LNnD44cOYKpU6eac/Ubzd/dWeoSiGqkUqmkLoGISDIqQRAsNkrsk08+QUxMDAoLC+Hg4IDXX38d8+fPh16vh6PjnRH8iYmJ2L17N7KzswEAY8aMwc2bN7Fnzx7xfnr37o2uXbti06ZNEAQB/v7+mDVrFp577jkAQEFBAXx8fJCcnIyxY8fi9OnT6NChA7799lv06NEDAJCSkoLhw4fj0qVL8Pf3r7V2g8EArVaLgoICaDQaU2+aagUnfmaxxyKqr/PLoy36eFK9DhujvObjWT5wdav8vXWXoXu1t03Vh9Z6/3WdX6e+Z0k1ZgCxOQ8X2VIvTVFoS7Pdd0FrdYNvez2ofl+eCgOL6tQuuOXVGpff75td4/JRmhNVXn/jehl6dMw1y/uGxcZE5eXlYdu2bejTpw8cHO78wG5aWhr69+8vBigAiIqKQk5ODq5duya2iYyMNLqvqKgopKWlAQDOnTsHvV5v1Ear1SIiIkJsk5aWBnd3dzFAAUBkZCTs7OyQnp5eZb2FhYUwGAxGFyIiIpJGXb5UWJrZQ9TcuXPh4uKC5s2b48KFC/j444/FZXq9Hj4+Pkbty//W6/U1tqm4vOLtqmvj7e1ttLxJkybw8PAQ29xt2bJl0Gq14iUgQJpB3k0dOV8UydM/+7eSugQyE/ZCyYM517cx+4nTHPyt3iEqMTERKpWqxkv5oTgAmD17Nk6ePIl9+/bB3t4ejz/+OCx4BLHB5s2bh4KCAvFy8eJFSepY92g3SR6XqDZx/UKkLoGISFL1DlGzZs3C6dOna7y0avX3N1RPT0/cc889uP/++7Fjxw7s3bsX33zzDQDA19cXubm5Rvdf/revr2+NbSour3i76tpcuXLFaHlJSQny8vLENndTq9XQaDRGFyn0ac2zn0ievN3kP0fUb7/9hsceewzNmzeHs7MzwsLCcPz4cXG5tZ7Za41srReqnK2ut1LUO0R5eXkhNDS0xkvFMU4VlZWVAbgz3ggAdDodjhw5guLiYrFNamoq2rVrh2bNmolt9u/fb3Q/qamp0Ol0AICQkBD4+voatTEYDEhPTxfb6HQ65OfnIyMjQ2xz4MABlJWVISIior6bwKKceTiPqEGuXbuGe++9Fw4ODvj888/x448/4pVXXhHfWwDrPbOXiCzDbL8cmp6ejm+//RZ9+/ZFs2bNcPbsWSxcuBCtW7cWw824ceOwZMkSxMXFYe7cucjMzMSaNWvw2muviffz7LPPYsCAAXjllVcQHR2NHTt24Pjx4+I0CCqVCtOnT8eLL76Itm3bIiQkBAsXLoS/vz9iYmIAAO3bt8fQoUMxZcoUbNq0CcXFxUhISMDYsWPrdGYeESnPyy+/jICAAGzZskW8LiTk70OQgiBg9erVWLBgAR566CEAwLvvvgsfHx/s3r1bPLM3JSXF6MzedevWYfjw4Vi1alWV7x+FhYXiF0UAsjgpxdp/v4zkTX3Bsc5n6CmN2QaWN23aFP/9738xePBgtGvXDnFxcejcuTMOHz4MtfrOqZVarRb79u3DuXPnEB4ejlmzZiEpKcnoW16fPn2wfft2bN68GV26dMGuXbuwe/dudOrUSWwzZ84cTJs2DVOnTkXPnj1x48YNpKSkwMnp78MN27ZtQ2hoKAYPHozhw4ejb9++RvNREVHd+Wrkfyjvk08+QY8ePTB69Gh4e3ujW7duePPNN8Xl5jqzVy4npVgTWz+kZevrL2dm64kKCwvDgQMHam3XuXNnfPXVVzW2GT16NEaPHl3tcpVKhaVLl2Lp0qXVtvHw8MD27dtrrUeOJuiCsDXtV6nLIBK9+kgXqUuo1S+//ILXX38dM2fOxPPPP49vv/0WzzzzDBwdHTFhwgSzndk7b948zJw5U/zbYDDYTJDiz4koi/ZsYaPmiyIzhigynacHtWGIIlnRtW4udQm1KisrQ48ePfCvf/0LANCtWzdkZmZi06ZNmDBhgtkeV61Wi73tRGTd+APECuCjgEMnZFuU8HMvfn5+6NChg9F17du3x4ULFwCY78xeMi0eyrqD20GeGKKIyCrde++9yMnJMbrup59+QlBQEADrP7OXiMyPh/OIyCrNmDEDffr0wb/+9S888sgjOHbsGDZv3swze2vB2aiJ6o49UQoR6usmdQlEAIDFD3SovZEM9OzZEx999BHee+89dOrUCS+88AJWr16N2NhYsQ3P7CVqGIbtO9gTpRCrRnfBiHX/k7oMIjwaESh1CXU2YsQIjBgxotrl1n5mr9JxHJAxx+xLKAptKXUZVAF7ohSivZ80PztDdDd1E86iLxfnL3lJXYKI0xsoE/db4zBEKYS9nfzPhiIiImooOX0pqCuGKCIiIqIGYIgiojqLVdB4KFI2joeqGreLvDBEKcimx8KlLoFs3LT72kpdAhGRbDBEKcjAdso7XkzWxVfL2fOJiMoxRCmIkwPPiiIiIpILhigiIiKiBmCIIqI66d3KQ+oSiIhkhSFKYZY+1FHqEshGvRjTSeoSyEbwDLSacfvIB0OUwjzcnVP+kzTaePP3G4mIKmKIUhhXNX/ukIiISA4YooiIqFH4+2vKxv3XcAxRRFSr1l4uUpdARCQ7DFEKNDuqndQlkI1ZPaab1CUQEckOQ5QCjdcFSV0C2ZiwllqpSyAikh2GKAXSODlIXQIREZHNY4giIiIiagCGKCKqUSsOKiciqhJDlEJxcDlZytqxHFRORFQVhiiFepyDy8lCOrXgoHIioqowRCmUGweXE5GJuf0qSF0Cfxeujrid5IEhioiIiOpNDqFbagxRRFStrgHuUpdARCRbDFEK9kJMJ6lLICu3anQXqUsgIpIthigFG9m9hdQlkJVr4+0qdQlERLJlkRBVWFiIrl27QqVS4bvvvjNa9sMPP6Bfv35wcnJCQEAAVqxYUen2O3fuRGhoKJycnBAWFoa9e/caLRcEAUlJSfDz84OzszMiIyNx5swZozZ5eXmIjY2FRqOBu7s74uLicOPGDZOvqyU1dWwidQlERGSjrgeppC5BchYJUXPmzIG/v3+l6w0GA4YMGYKgoCBkZGRg5cqVWLx4MTZv3iy2OXr0KB599FHExcXh5MmTiImJQUxMDDIzM8U2K1aswNq1a7Fp0yakp6fDxcUFUVFRuH37ttgmNjYWWVlZSE1NxZ49e3DkyBFMnTrVvCtORERkBkWhLaUugWCBEPX5559j3759WLVqVaVl27ZtQ1FREd555x107NgRY8eOxTPPPINXX31VbLNmzRoMHToUs2fPRvv27fHCCy+ge/fuWL9+PYA7vVCrV6/GggUL8NBDD6Fz58549913cfnyZezevRsAcPr0aaSkpOCtt95CREQE+vbti3Xr1mHHjh24fPmyuTcBkSI93I2Hi22NHHoWGA5IScwaonJzczFlyhT8+9//RtOmTSstT0tLQ//+/eHo6CheFxUVhZycHFy7dk1sExkZaXS7qKgopKWlAQDOnTsHvV5v1Ear1SIiIkJsk5aWBnd3d/To0UNsExkZCTs7O6Snp1dZe2FhIQwGg9FFjt6Z2KP2RkQNMG94e6lLICKSNbOFKEEQMHHiRDz55JNG4aUivV4PHx8fo+vK/9br9TW2qbi84u2qa+Pt7W20vEmTJvDw8BDb3G3ZsmXQarXiJSAgoNZ1lkLfNl5Sl0BWystNLXUJRESyVu8QlZiYCJVKVeMlOzsb69atw/Xr1zFv3jxz1G128+bNQ0FBgXi5ePGi1CVVybEJT7AkIiKSQr1P75o1axYmTpxYY5tWrVrhwIEDSEtLg1pt/G22R48eiI2NxdatW+Hr64vc3Fyj5eV/+/r6iv9W1abi8vLr/Pz8jNp07dpVbHPlyhWj+ygpKUFeXp54+7up1epKtRMRERGVq3eI8vLygpdX7YeQ1q5dixdffFH8+/Lly4iKisL777+PiIgIAIBOp8P8+fNRXFwMB4c7vwWXmpqKdu3aoVmzZmKb/fv3Y/r06eJ9paamQqfTAQBCQkLg6+uL/fv3i6HJYDAgPT0dTz31lHgf+fn5yMjIQHh4OADgwIEDKCsrE2tRslBfN2Trr0tdBlmRhSM6SF0CEZHsme1YUGBgIDp16iRe7rnnHgBA69at0bLlnbMvxo0bB0dHR8TFxSErKwvvv/8+1qxZg5kzZ4r38+yzzyIlJQWvvPIKsrOzsXjxYhw/fhwJCQkAAJVKhenTp+PFF1/EJ598glOnTuHxxx+Hv78/YmJiAADt27fH0KFDMWXKFBw7dgxff/01EhISMHbs2CqnXlCaDbHdpS6BrMzjuiCpSyAikj1JZ2vUarXYt28f4uPjER4eDk9PTyQlJRnN39SnTx9s374dCxYswPPPP4+2bdti9+7d6NTp7588mTNnDm7evImpU6ciPz8fffv2RUpKCpycnMQ227ZtQ0JCAgYPHgw7OzuMHDkSa9eutej6mktrL84qTablYM+xdkREtbFYiAoODoYgVP7F586dO+Orr76q8bajR4/G6NGjq12uUqmwdOlSLF26tNo2Hh4e2L59e90LJiIiIqoBv24SEVGjFLTmSThKxv3XcAxRVuKJe0OkLoGsxPbJyj/ZgojIEhiirMSzkW2lLoGsRJ82nlKXQESkCAxRVkLr7CB1CUREJsHfz6sZt498MEQRERERNQBDFBGJHuyi/HnTiIgshSHKimzjgGBqpMUPdpS6BCIixWCIsiK9WzWXugRSOA8XR6lLICJSDIYoK2Jvp5K6BCIiIpvBEEVERLLDM9Cqxu0iLwxRVuaRHnyBUcNsmdhT6hKIiBSFIcrKPD+8vdQlkEINCvWWugSqp+CWV6UugRROTj/5osTnM0OUlXFvyoHBRGR5cvowJrIUhigiIpIljv8xxu0hPwxRRITxvYOkLoFk4noQz/Kl2vF5cgdDlBX6NKGv1CWQwiQOC5W6BCIixWGIskKdWmikLoEUxkXdROoSiIgUhyHKCqlU7GYlIuvAcUB3cDvIE0MUkY3zdOUZnUREDcEQZaVeeIg/JEt1896U3lKXQFaC0xwoC/dX4zFEWalHegZIXQIpRFsfN6lLIDMrDCySuoRGsfVDWUpff6U//2rCEGWl1E3spS6BiIjIqjFEERGR7Cm9N6ahbHW9lYIhyoo92itQ6hJI5vijw0REDccQZcWeH84JFKlmtvKjw8uXL4dKpcL06dPF627fvo34+Hg0b94crq6uGDlyJHJzc41ud+HCBURHR6Np06bw9vbG7NmzUVJSYuHqieSFs5X/jSHKirk5OUhdApHkvv32W7zxxhvo3Lmz0fUzZszAp59+ip07d+Lw4cO4fPkyHn74YXF5aWkpoqOjUVRUhKNHj2Lr1q1ITk5GUlKSpVfB4hrzIWnOM75s7dCWOdeXZ+aZBkMUEVmtGzduIDY2Fm+++SaaNWsmXl9QUIC3334br776Ku677z6Eh4djy5YtOHr0KL755hsAwL59+/Djjz/iP//5D7p27Yphw4bhhRdewIYNG1BUZJmzje73zbbI4xApgRxfDwxRVq5vG0+pSyCZWjGqc+2NFC4+Ph7R0dGIjIw0uj4jIwPFxcVG14eGhiIwMBBpaWkAgLS0NISFhcHHx0dsExUVBYPBgKysrGofs7CwEAaDwehCpmMrvVG2sp5KxxBl5V4b01XqEkimRnW37jfpHTt24MSJE1i2bFmlZXq9Ho6OjnB3dze63sfHB3q9XmxTMUCVLy9fVp1ly5ZBq9WKl4AAztlGZK0YoqyclxuPe1PV7Oysd3DoxYsX8eyzz2Lbtm1wcnKy6GPPmzcPBQUF4uXixYsWfXw5MPd4G2vvpTH3+nE8lOkwRBGR1cnIyMCVK1fQvXt3NGnSBE2aNMHhw4exdu1aNGnSBD4+PigqKkJ+fr7R7XJzc+Hr6wsA8PX1rXS2Xvnf5W2qolarodFojC5EZJ0YomxAUPOmUpdAMjNnaDupSzCrwYMH49SpU/juu+/ES48ePRAbGyv+38HBAfv37xdvk5OTgwsXLkCn0wEAdDodTp06hStXrohtUlNTodFo0KFDB4uvExmz1t4oa10va9VE6gLI/N6b0ht9lh+QugySkX/2by11CWbl5uaGTp06GV3n4uKC5s2bi9fHxcVh5syZ8PDwgEajwbRp06DT6dC7950fZB4yZAg6dOiA8ePHY8WKFdDr9ViwYAHi4+OhVlv/4ZDrQSq4/SpIXUaNikJbwjH7ktRlmIwSAhTniDJm1p6o4OBgqFQqo8vy5cuN2vzwww/o168fnJycEBAQgBUrVlS6n507dyI0NBROTk4ICwvD3r17jZYLgoCkpCT4+fnB2dkZkZGROHPmjFGbvLw8xMbGQqPRwN3dHXFxcbhx44bpV1qG/N2dpS6BZMbeisdD1dVrr72GESNGYOTIkejfvz98fX3x3//+V1xub2+PPXv2wN7eHjqdDo899hgef/xxLF26VMKqG87SPwLLcTfyZOn9Ys0/PgxYoCdq6dKlmDJlivi3m9vfvxhvMBgwZMgQREZGYtOmTTh16hSeeOIJuLu7Y+rUqQCAo0eP4tFHH8WyZcswYsQIbN++HTExMThx4oT4jXLFihVYu3Yttm7dipCQECxcuBBRUVH48ccfxUGlsbGx+P3335Gamori4mJMmjQJU6dOxfbt2829CYhIBg4dOmT0t5OTEzZs2IANGzZUe5ugoKBKX9pIXqylN0oJvVBUmdnHRLm5ucHX11e8uLi4iMu2bduGoqIivPPOO+jYsSPGjh2LZ555Bq+++qrYZs2aNRg6dChmz56N9u3b44UXXkD37t2xfv16AHd6oVavXo0FCxbgoYceQufOnfHuu+/i8uXL2L17NwDg9OnTSElJwVtvvYWIiAj07dsX69atw44dO3D58uUq67a2uV68eZYe/b/pkW2lLoFMKLjlValLkJzSA4jS6zcFpT6PzR6ili9fjubNm6Nbt25YuXKl0e9OpaWloX///nB0dBSvi4qKQk5ODq5duya2uXuivKioKHFCvHPnzkGv1xu10Wq1iIiIMJo0z93dHT169BDbREZGws7ODunp6VXWbW1zveyOv1fqEkgmpt3HEEXWR6lBRKl10x1mDVHPPPMMduzYgYMHD+Kf//wn/vWvf2HOnDni8rpMZlddm4rLK96uujbe3sY/tNqkSRN4eHhUO2metc31wnFRVI7jochSOC5KXrg/TK/eY6ISExPx8ssv19jm9OnTCA0NxcyZM8XrOnfuDEdHR/zzn//EsmXLZH92i1qtln2NRETmpIQz9CpS2vgopfVC8cy8yuodombNmoWJEyfW2KZVq1ZVXh8REYGSkhKcP38e7dq1q9NkdtW1qbi8/Do/Pz+jNl27dhXbVJzrBQBKSkqQl5dX46R51qa9nwanf1f22C5qnIUjOL8RWTelBCmlBSiqWr0P53l5eSE0NLTGS8UxThV99913sLOzEw+t6XQ6HDlyBMXFxWKb1NRUtGvXTvzFdZ1OZzQhXnmb8gnxQkJC4Ovra9TGYDAgPT3daNK8/Px8ZGRkiG0OHDiAsrIyRERE1HcTKNa7T/SSugSS2KQ+wVKXQDZGikNIcg8oUtTHQ3nmYbYpDtLS0pCeno5BgwbBzc0NaWlpmDFjBh577DExII0bNw5LlixBXFwc5s6di8zMTKxZswavvfaaeD/PPvssBgwYgFdeeQXR0dHYsWMHjh8/js2bNwMAVCoVpk+fjhdffBFt27YVpzjw9/dHTEwMAKB9+/YYOnQopkyZgk2bNqG4uBgJCQkYO3Ys/P39zbUJZIe/o0fW/Ht5RBWVBxU59UrJPdxR/ZktRKnVauzYsQOLFy9GYWEhQkJCMGPGDKNxUlqtFvv27UN8fDzCw8Ph6emJpKQkcY4oAOjTpw+2b9+OBQsW4Pnnn0fbtm2xe/duo9mI58yZg5s3b2Lq1KnIz89H3759kZKSYvTDo9u2bUNCQgIGDx4MOzs7jBw5EmvXrjXX6hMRyU5hYBHUF6o+UmCt5HJ4zxYDlLVPtAkAKkEQlDNqUEIGgwFarRYFBQWK/kHR+G0n8Nmp36UugySwYVx3RHf2q72hjCnxdVhe8/EsH7i6VR5BscvQvcbbp+pDa32M85e86lxPQ0KUKQaXa88WNvo+GkPKICV1gDLFobyGDCqvT4iqyzxR9/tm17h8lOZEldffuF6GHh1zzfK+wR8gtjHLR4ZJXQJJZHiY7ZxEQXS3otCWFg8zUjwmWRZ/gNjGuDk5SF0CSUSl4ngoIkuMlZJTcOKAcvNiiCIiIrMraK2W/JBeReYIU3IKT2QZDFE2aOb99+DV1J+kLoMsiD/7Q1S1isGnIYGKwcm2MUTZoCcHtGaIsjFdA9ylLoHMKLjl1XoNLq8vU81cLrfeqLtZWyAy1aE8c89UrtQfHwY4sNwmOTbhbiciImosfpoSWblWni5Sl0AyIvXcPRzobBlSb2epn2eWwhBlo/4TZzs/d2PrPnhSJ3UJRERWiSHKRt3bprnUJZCFeLrym7+S1TbBoKWYclyM1L0k1s6U29fc46HqSi6vg7sxRNkozhlERETUOAxRNiyoeVOpSyAzmzO0ndQlEFWJvVHmwe1qWQxRNuyjpzl3kLV7akBrqUsgKyKXQztkGdzftWOIsmEeLrb1a+62iIdtqSpyOXOKvSamJZftKZfnlyUwRBERkWTk8sGvdNyO0mCIsnEvxHSSugQyk30z+ktdAlmQpWZ95iEe22Cp/azk2coBhiibF9srUOoSyEzu8XGTugSqg1GaE1KXIDn2ojQOt590ryOGKBtnZ8dvlUQkPQaBhuF2kxZDFEHN39KzOnF9Q6QugWSuMYN/zXWoh4Ggfsy1vRqzf21pUDnAEEUA0uYNlroEMrH5w9tLXQIRkdVjiCJOdWCFeJiWlIq9UXXD7SQPDFFERFRv5jx7iwGhZubcPjz7sn4YoggAsDG2u9QlkIl8NWeQ1CWQQsh5/AqDVNXkvF3k/HwyF4YoAgAM6+QrdQlkIgEe/E1EW2XpOXfM3Wsh58AgBXNvD0v3Qil9jiiAIYr+H38ehIjkiEHqDm4HeWKIItHD3VpIXQI10ubx4VKXQGZwv2+21CVIytYDhK2vv5yf/wxRJFo2MkzqEqiRhnTkYVmqn8aOY7HUISBbDRKWWu/G7kdbHA8FAE2kLoDkQ93EXuoSiIiqVR4otGcLJa7E/Gw1NCoNe6KIrMTI7i2lLoFslKUHJFt7wLD0+nFag4ZjiCIjGQsipS6BGmjlqM5Sl0A2jEHKNBiglIUhiow0d7XONyZbwFnKCWjYaeNKHc9S0FptNWFKyevSkOePNUxvAHBMFJFVCPBwlroEIlwPUsHtV8Hij6vksVJSBif2QjUee6KokpTp/aQugerp82f7S10CkeSU1JujpFqpemYNUZ999hkiIiLg7OyMZs2aISYmxmj5hQsXEB0djaZNm8Lb2xuzZ89GSUmJUZtDhw6he/fuUKvVaNOmDZKTkys9zoYNGxAcHAwnJydERETg2LFjRstv376N+Ph4NG/eHK6urhg5ciRyc3NNvbpWI9RXI3UJVE+uanYqU+OY6pCeHHo35BxQ5FKbqfaTUg8Fm4rZQtSHH36I8ePHY9KkSfj+++/x9ddfY9y4ceLy0tJSREdHo6ioCEePHsXWrVuRnJyMpKQksc25c+cQHR2NQYMG4bvvvsP06dMxefJkfPHFF2Kb999/HzNnzsSiRYtw4sQJdOnSBVFRUbhy5YrYZsaMGfj000+xc+dOHD58GJcvX8bDDz9srlUnIrJpcghSwN+BRerQIpc6ysll/1gDlSAIJj+AXVJSguDgYCxZsgRxcXFVtvn8888xYsQIXL58GT4+PgCATZs2Ye7cubh69SocHR0xd+5cfPbZZ8jMzBRvN3bsWOTn5yMlJQUAEBERgZ49e2L9+vUAgLKyMgQEBGDatGlITExEQUEBvLy8sH37dowaNQoAkJ2djfbt2yMtLQ29e/eu0zoZDAZotVoUFBRAo7H+nprDP13FhHeO1d6QJJf+/GD4aJykLsMilPg6LK/5eJYPXN2q/t66y1D7D4Cn6kPr9bjnL3nVqz0AqC841vs21ZFibFRdWWLslFwCU1VMGaIsMai8LjOWj9KcqHbZjetl6NEx1yzvG2bpiTpx4gR+++032NnZoVu3bvDz88OwYcOMwlBaWhrCwsLEAAUAUVFRMBgMyMrKEttERhqfch8VFYW0tDQAQFFRETIyMoza2NnZITIyUmyTkZGB4uJiozahoaEIDAwU21SlsLAQBoPB6GJLBtxT/zdgkoatBChbZ4mfvjDloRk593ZU7BlqbC+RKe/LEqQOUPUl5598Acx0dt4vv/wCAFi8eDFeffVVBAcH45VXXsHAgQPx008/wcPDA3q93ihAARD/1uv14r9VtTEYDPjrr79w7do1lJaWVtkmOztbvA9HR0e4u7tXalP+OFVZtmwZlixZUv+VJyIiANKdrddQcg4/piDnYKtU9eqJSkxMhEqlqvGSnZ2NsrIyAMD8+fMxcuRIhIeHY8uWLVCpVNi5c6dZVsTU5s2bh4KCAvFy8eJFqUuyuNdjaz/EQNL6as4gqUsgIrJZ9eqJmjVrFiZOnFhjm1atWuH3338HAHTo0EG8Xq1Wo1WrVrhw4QIAwNfXt9JZdOVnzPn6+or/3n0WXW5uLjQaDZydnWFvbw97e/sq21S8j6KiIuTn5xv1RlVsUxW1Wg212rq/ldRmWJif1CVQLQI8mkpdAlmZwsAik46NUlpvlLUydS+UrZ+VV65ePVFeXl4IDQ2t8eLo6Ijw8HCo1Wrk5OSIty0uLsb58+cRFBQEANDpdDh16pTRWXSpqanQaDRi+NLpdNi/f79RDampqdDpdAAgPlbFNmVlZdi/f7/YJjw8HA4ODkZtcnJycOHCBbENEZG1kdOM0DyMJC05bX85PS9NwSwDyzUaDZ588kksWrQI+/btQ05ODp566ikAwOjRowEAQ4YMQYcOHTB+/Hh8//33+OKLL7BgwQLEx8eLPUBPPvkkfvnlF8yZMwfZ2dnYuHEjPvjgA8yYMUN8rJkzZ+LNN9/E1q1bcfr0aTz11FO4efMmJk2aBADQarWIi4vDzJkzcfDgQWRkZGDSpEnQ6XR1PjPPlm16LFzqEqgaPJRH5mKOXgY5fZDbEnNsd/ZC/c1sM/StXLkSTZo0wfjx4/HXX38hIiICBw4cQLNmzQAA9vb22LNnD5566inodDq4uLhgwoQJWLp0qXgfISEh+OyzzzBjxgysWbMGLVu2xFtvvYWoqCixzZgxY3D16lUkJSVBr9eja9euSElJMRps/tprr8HOzg4jR45EYWEhoqKisHHjRnOtulUZ2qn6Q54kLR7KI6XhoT3LYnA1P7PME2WNlDg/jakEJ34mdQlUhfPLo6UuweKU+DqsyzxRgHzmiqrIlGOjKmKQMj9zBajG9kJZeo4oQIHzRJF12T45QuoS6C7H5g+WugSSQH3nzJHr+BP2kJiXXLevOQKU1BiiqFZ92nhKXQLdxduNE2yS+Zlz7ItcP+iVzpzblWOhKmOIIlKYpo72UpdAZBIMUqZzPUjF7SkBhiiqk4PPDZS6BPp/GQvul7oEsiHm7n3gh3/jWWL7sReqagxRVCchni5Sl0D/z5k9UVQPch0XdTcGqYZRynZTyvOwvsw2xQERmd7wME45QZZn6lnMq1MeCHj2Xu0sGZ7YC1U99kRRnWUuiaq9EZnVhnH8PUOShiU/SJXSuyIVBij5YIiiOnNVs+NSaioVP1xsXUNO+1bioRSOlapMqdukIc8/JUxvADBEUT092MVf6hJs1vpx3aQugWycFL0SSg0OpiTVNmAvVO0YoqheVo/pKnUJNmtEZwbY+li2bBl69uwJNzc3eHt7IyYmxuhH0QHg9u3biI+PR/PmzeHq6oqRI0ciNzfXqM2FCxcQHR2Npk2bwtvbG7Nnz0ZJSYklV0VWpPpgtcUwJeU6M0DVDUMU1YudnW29iZFyHT58GPHx8fjmm2+QmpqK4uJiDBkyBDdv3hTbzJgxA59++il27tyJw4cP4/Lly3j44YfF5aWlpYiOjkZRURGOHj2KrVu3Ijk5GUlJSVKsUqMo8ZBeVWwhTFnTOlrL8646HORC9fafuAg89na61GXYlPTn+TMv9ZWSkmL0d3JyMry9vZGRkYH+/fujoKAAb7/9NrZv34777rsPALBlyxa0b98e33zzDXr37o19+/bhxx9/xJdffgkfHx907doVL7zwAubOnYvFixfD0dH8Z6zJkaXO1qtJxZBhDWfzySk0sReq7tgTRfXWty1/BsbSfDT8mZfGKigoAAB4eHgAADIyMlBcXIzIyEixTWhoKAIDA5GWlgYASEtLQ1hYGHx8fMQ2UVFRMBgMyMrKqvJxCgsLYTAYjC5yYcpeATl90Jb33MgpiNSFHOs25X5t6PNNKYPKAYYoItkb2pFzQzVWWVkZpk+fjnvvvRedOnUCAOj1ejg6OsLd3d2orY+PD/R6vdimYoAqX16+rCrLli2DVqsVLwEBASZeG/l8yMgpSJWrGEzkFE4AedcGyHN/yh1DFDXI6aVDpS7BZrz+GOeGaqz4+HhkZmZix44dZn+sefPmoaCgQLxcvHjR7I8pJbl/8N4dXCwRYKR4zMaS+36UK46JogbhT49YDueGapyEhATs2bMHR44cQcuWLcXrfX19UVRUhPz8fKPeqNzcXPj6+optjh07ZnR/5Wfvlbe5m1qthlqtrnedozQnsMtg/sAc3PIqzl/yMul9ymGMVH3JPdRYkjkClKUGlI/SnLDI41SHPVHUYCtGdpa6BKv3xfT+UpegWIIgICEhAR999BEOHDiAkJAQo+Xh4eFwcHDA/v37xetycnJw4cIF6HQ6AIBOp8OpU6dw5coVsU1qaio0Gg06dOhgmRVRCPZkKBP3W+MwRFGDPdLT9GM9yFg7XzepS1Cs+Ph4/Oc//8H27dvh5uYGvV4PvV6Pv/76CwCg1WoRFxeHmTNn4uDBg8jIyMCkSZOg0+nQu3dvAMCQIUPQoUMHjB8/Ht9//z2++OILLFiwAPHx8Q3qbTKlxoyLMlcvAT+QlcVc+6sxzy+5jPerKx7OI5KpQe1Me8jF1rz++usAgIEDBxpdv2XLFkycOBEA8Nprr8HOzg4jR45EYWEhoqKisHHjRrGtvb099uzZg6eeego6nQ4uLi6YMGECli5daqnVUBwlHtqzRQy8psEQRY2S/cJQhC5Mqb0h1ds7E3tKXYKiCULtcwc5OTlhw4YN2LBhQ7VtgoKCsHfvXlOWJgvmGBtVrvwDmmFKfswdnqx9cs278XAeNYqTAweYmwsHlJO5mfsDj70d8sIAZXoMUdRoH/xTJ3UJVufYfM5QTrVTwvgRBil5UMJ+UMLz+W4MUdRovUI8pC7B6ni7cYZysgxL9B4UBhYp4kPcGllq29tiLxTAEEUm8lBXf6lLsBprH+0mdQlEZsEwZTnc1pbBEEUmsXpMV6lLsBoPdmEgpbozxSEQS/ci8MPdvCy9fU3x/FHioTyAIYpMhIOgTaNrgLvUJZCEpJx9WYogxd4S05Fqe0p5GE/q2coBhigyoewX+Ht6jfXR032kLoFsmFQfiAxTDSfltrPVcVAVMUSRyXC6g8Zjjx41hCkPhUj5wcgwVXdSbytTPk+UeigPYIgiE/tmHk/Nb6icF9mTRwTwUF91uF3khzOWk0n5anlqfkOpm7Anj+TBnLOZ11fFwGCLM6DLMTDxMN7f2BNFJvdJwr1Sl6A4JxbeL3UJpHCmPiQixw/Kij0xcgwXpiD3dTT180LJh/IA9kSRGXRu6S51CYrj4WJ737BJ/uTUI1WVqkKGknqr5BiSaiLHYC01higyi7cn9EDc1uNSl6EIR2YPkroEshL3+2YjVR9q0vss/+CUc5iqqKZgIkXAUlpQqoq5wlNjeqHkML0BYMbDeYcOHYJKpary8u2334rtfvjhB/Tr1w9OTk4ICAjAihUrKt3Xzp07ERoaCicnJ4SFhVX6RXVBEJCUlAQ/Pz84OzsjMjISZ86cMWqTl5eH2NhYaDQauLu7Iy4uDjdu3DDPyhMGt/eRugTFCGzeVOoSSEbk8uFwN2vohbj7UJklLkpnDfvdnMwWovr06YPff//d6DJ58mSEhISgR48eAACDwYAhQ4YgKCgIGRkZWLlyJRYvXozNmzeL93P06FE8+uijiIuLw8mTJxETE4OYmBhkZmaKbVasWIG1a9di06ZNSE9Ph4uLC6KionD79m2xTWxsLLKyspCamoo9e/bgyJEjmDp1qrlWnwBsjO0udQmyt3/WAKlLICtjzjEm/EC1Lebc30ofC1VOJQiCYIkHKi4uRosWLTBt2jQsXLgQAPD6669j/vz50Ov1cHS8082amJiI3bt3Izv7zgYeM2YMbt68iT179oj31bt3b3Tt2hWbNm2CIAjw9/fHrFmz8NxzzwEACgoK4OPjg+TkZIwdOxanT59Ghw4d8O2334oBLiUlBcOHD8elS5fg71/5ZzYKCwtRWFgo/m0wGBAQEICCggJoNBrzbCQrFJz4mdQlyNr55dFSl6AoBoMBWq1WUa/D8pqPZ/nA1a1u31t3GRr3BcTUh/TuppRDe9QwlgjLjQ1R9emxvXG9DD065prlfcNiZ+d98skn+PPPPzFp0iTxurS0NPTv318MUAAQFRWFnJwcXLt2TWwTGRlpdF9RUVFIS0sDAJw7dw56vd6ojVarRUREhNgmLS0N7u7uYoACgMjISNjZ2SE9Pb3KepctWwatViteAgICGrkFbNOWiT2lLkG2vprDsVBkHub+lh/c8ip7payUEgKUnFgsRL399tuIiopCy5Ytxev0ej18fIzHzpT/rdfra2xTcXnF21XXxtvb22h5kyZN4OHhIba527x581BQUCBeLl68WK/1pTsGhXrX3shGBXhwLBQpG8OU9eC+bJh6h6jExMRqB4yXX8oPxZW7dOkSvvjiC8TFxZmscHNTq9XQaDRGF2qYvc/0k7oE2fk+aYjUJZCMmWJwuSW/7fMDWLksve9M8byU08kX9Z7iYNasWZg4cWKNbVq1amX095YtW9C8eXM8+OCDRtf7+voiNzfX6Lryv319fWtsU3F5+XV+fn5Gbbp27Sq2uXLlitF9lJSUIC8vT7w9mU8HfwbQu2mbOkhdAtkAc0x5UBOlTYdgy6QIvdZ0GK9cvXuivLy8EBoaWuOl4hgnQRCwZcsWPP7443BwMP7g0Ol0OHLkCIqLi8XrUlNT0a5dOzRr1kxss3//fqPbpaamQqfTAQBCQkLg6+tr1MZgMCA9PV1so9PpkJ+fj4yMDLHNgQMHUFZWhoiIiPpuAmoAzsj9tzMvDZO6BCKzYs+UfHHfmJbZx0QdOHAA586dw+TJkystGzduHBwdHREXF4esrCy8//77WLNmDWbOnCm2efbZZ5GSkoJXXnkF2dnZWLx4MY4fP46EhAQAd371fvr06XjxxRfxySef4NSpU3j88cfh7++PmJgYAED79u0xdOhQTJkyBceOHcPXX3+NhIQEjB07tsoz88j0OCP3He5NHeBgz19bIsuR8ts/P7DlQ+p9YY29UIAFQtTbb7+NPn36IDS0cpeyVqvFvn37cO7cOYSHh2PWrFlISkoymr+pT58+2L59OzZv3owuXbpg165d2L17Nzp16iS2mTNnDqZNm4apU6eiZ8+euHHjBlJSUuDk9PeP4W7btg2hoaEYPHgwhg8fjr59+xrNR0Xmd/Zfw6UuQXIn2SNHEpD6A6z8A5yByrLkst1N+fyT03gowILzRCmdEuenkaMp7x5H6o+5tTe0Qi/9oxNiI4KkLkPRlPg6bMg8UeUaO1/U3Sw5PqouOHbK9KQOTHczdYBvSIgy5zxR/O08sqg3H+9hsxNwMkCR1Cw90Lw2d3/gM1TVn9xCU0VS94BaAkMUWdyeaX0xYt3/pC7DojiwnuRCbkGqIoaq2sk5NFVkCwEKYIgiCXRqoZW6BIvjwHpqiFGaEyY/pAfIO0hVVF1gsIVwpZSwVBVzBSi5jYcCGKJIIr/8azhaPb9X6jIs4twyDqgn+VFKkKpKTQFDSQFLyUGpOrbSA1WOIYokYWenwkv/6IT5H2VKXYpZ7Y6/FyqVSuoyiKqk5CBVnYYEE1MEL2sMRPVlawEKYIgiCcVGBFl9iOoa4C51CaRw5jqkV678g8/awlR9MAA1jiXCkxwP5QEW/AFioqr8bMWzd/MwHinJ/b7ZNtmTQA3H5wxDFEmsib0d3pnYQ+oyTO7gcwN5GI9MxpLfwm39Q5HqxpLPE7n2QgE8nEcycF+oj9QlmFSXllqEeLpIXQZRg/EQH1WHIdsYQxTJwvnl0VYzCefHCX2lLoGskLnHRlWl4gcmA5XtkjI4ybkXCuDhPJKRH5dGSV1Co/3C3wckK8XxL7aH+7x27Iki2Wjq2ATvTemNR9/8RupSGuR/cwfBzo7joMh8pOiNuht7p6ybnEKT3HuhAIYokhld6+Z4sIs/Pvn+stSl1Mu//hGGls2aSl0GkUUxUFkHOQUnpWGIItlZ+2g3RYWoLgHuGBcRKHUZZCPk0BtVlbs/iBmq5EsJoUkJvVAAQxTJlJIGmn8cf6/UJZCNkWuQqqiqD2oGK8tTQmC6m1ICFMAQRTKmhCB1fnm01CWQjVJCkLpbdR/oDFeNp8SwVBUlBSiAIYpkTs5BigGKpKbEIFWV2gIAQ5b1hKSaKC1AAQxRpAByDFIMUCQX1hKkalLfAKGE0GULoag+lBigAIYoUgi5BKnOLbX4hJNpkszYQpCqDwYUZVFqgAI42SYpyPnl0ZjSL0Syx18ztisDFMnWKM0JRX8Yke2xhucsQxQpyvzoDvg68T6LP27mkig81LWFxR+XqL6U/qFE1s8awlM5hihSnBbuzhYbk9TcxRHnl0fDVc0j36Qc1vQhRdbF2p6X/GQgxTq/PBpXrt9Gr5f2m+X+Ty8dCmdHe7PcN5ElVPzA4pgpkoq1BaeKGKJI0bzdnHB+eTRO/27AsDVfmeQ+058fDB+Nk0nui0guGKjIkqw5OFXEEEVWob2fBueXR6OsTECnxV/gVlFpvW7fPdAd/32aM4+TbWCgInOwleBUEUMUWRU7OxV+XDpU/PtmYQm2fH0Oq/b9JF6nUgFLHuyIR3oEwMmBh+vItt39wcdQRXVli6HpbgxRZNVc1E2QcF9bJNzXVupSiBShqg9GBitiYKoaQxQREdWopg9QBizrwaBUfwxRRETUYHX94GXYkg7DkfkwRBERkdmZ4oPcFoMYA5C8MUQREZEiMFCQ3HDGciIiIqIGMFuI+umnn/DQQw/B09MTGo0Gffv2xcGDB43aXLhwAdHR0WjatCm8vb0xe/ZslJSUGLU5dOgQunfvDrVajTZt2iA5ObnSY23YsAHBwcFwcnJCREQEjh07ZrT89u3biI+PR/PmzeHq6oqRI0ciNzfX5OtMREREtsNsIWrEiBEoKSnBgQMHkJGRgS5dumDEiBHQ6/UAgNLSUkRHR6OoqAhHjx7F1q1bkZycjKSkJPE+zp07h+joaAwaNAjfffcdpk+fjsmTJ+OLL74Q27z//vuYOXMmFi1ahBMnTqBLly6IiorClStXxDYzZszAp59+ip07d+Lw4cO4fPkyHn74YXOtOhEREdkAlSAIgqnv9I8//oCXlxeOHDmCfv36AQCuX78OjUaD1NRUREZG4vPPP8eIESNw+fJl+Pj4AAA2bdqEuXPn4urVq3B0dMTcuXPx2WefITMzU7zvsWPHIj8/HykpKQCAiIgI9OzZE+vXrwcAlJWVISAgANOmTUNiYiIKCgrg5eWF7du3Y9SoUQCA7OxstG/fHmlpaejdu3ed1slgMECr1aKgoAAajcZk24qI6k6Jr8Pymo9n+cDVjSMoiCztxvUy9OiYa5b3DbO8ops3b4527drh3Xffxc2bN1FSUoI33ngD3t7eCA8PBwCkpaUhLCxMDFAAEBUVBYPBgKysLLFNZGSk0X1HRUUhLS0NAFBUVISMjAyjNnZ2doiMjBTbZGRkoLi42KhNaGgoAgMDxTZVKSwshMFgMLoQERERlTPL2XkqlQpffvklYmJi4ObmBjs7O3h7eyMlJQXNmjUDAOj1eqMABUD8u/yQX3VtDAYD/vrrL1y7dg2lpaVVtsnOzhbvw9HREe7u7pXalD9OVZYtW4YlS5bUf+WJiIjIJtSrJyoxMREqlarGS3Z2NgRBQHx8PLy9vfHVV1/h2LFjiImJwQMPPIDff//dXOtiUvPmzUNBQYF4uXjxotQlEZFEajt5hYhsU716ombNmoWJEyfW2KZVq1Y4cOAA9uzZg2vXronHHzdu3IjU1FRs3boViYmJ8PX1rfRGVH7GnK+vr/jv3WfR5ebmQqPRwNnZGfb29rC3t6+yTcX7KCoqQn5+vlFvVMU2VVGr1VCr1TWuKxFZv/KTVzZt2oSIiAisXr0aUVFRyMnJgbe3t9TlEZGE6tUT5eXlhdDQ0Bovjo6OuHXr1p07tzO+ezs7O5SVlQEAdDodTp06ZXQWXWpqKjQaDTp06CC22b9/v9F9pKamQqfTAQAcHR0RHh5u1KasrAz79+8X24SHh8PBwcGoTU5ODi5cuCC2ISKqzquvvoopU6Zg0qRJ6NChAzZt2oSmTZvinXfekbo0IpKYWcZE6XQ6NGvWDBMmTEBSUhKcnZ3x5ptvilMWAMCQIUPQoUMHjB8/HitWrIBer8eCBQsQHx8v9gA9+eSTWL9+PebMmYMnnngCBw4cwAcffIDPPvtMfKyZM2diwoQJ6NGjB3r16oXVq1fj5s2bmDRpEgBAq9UiLi4OM2fOhIeHBzQaDaZNmwadTlfnM/MAoPwkRg4wJ5JO+evPDCcVV6n85JV58+aJ19198srdCgsLUVhYKP5dUFAAALhxo8y8xRJRlcpfe2Z53xDM5NtvvxWGDBkieHh4CG5ubkLv3r2FvXv3GrU5f/68MGzYMMHZ2Vnw9PQUZs2aJRQXFxu1OXjwoNC1a1fB0dFRaNWqlbBly5ZKj7Vu3TohMDBQcHR0FHr16iV88803Rsv/+usv4emnnxaaNWsmNG3aVPjHP/4h/P777/Van4sXLwoAeOGFFxlcLl68WK/Xb0P99ttvAgDh6NGjRtfPnj1b6NWrV5W3WbRokeTbhxdeeKl8Mcf7hlnmibJGZWVluHz5Mtzc3KBSqaptZzAYEBAQgIsXLypmHpu7cR3kgetQmSAIuH79Ovz9/SsNFzCHy5cvo0WLFjh69KjR4f85c+bg8OHDSE9Pr3Sbu3ui8vPzERQUhAsXLkCr1Zq9ZkuyhudodbhuylTVupnzfYM/QFxHdnZ2aNmyZZ3bazQaxT85uQ7ywHUwZskg4unpWevJK3er7qQUrVar+P1YHWt4jlaH66ZMd6+bud43OH0uEVE16nLyChHZLvZEERHVoLaTV4jIdjFEmZharcaiRYsUPccU10EeuA7yMGbMGFy9ehVJSUnQ6/Xo2rUrUlJSKv1SQnWsYRtUh+umTFw30+HAciIiIqIG4JgoIiIiogZgiCIiIiJqAIYoIiIiogZgiCIiIiJqAIYoIiIiogZgiKqDDRs2IDg4GE5OToiIiMCxY8dqbL9z506EhobCyckJYWFh2Lt3r9FyQRCQlJQEPz8/ODs7IzIyEmfOnDHnKtRrHd58803069cPzZo1Q7NmzRAZGVmp/cSJE6FSqYwuQ4cOlc06JCcnV6rPycnJqI3c98PAgQMrrYNKpRJ/xBuw7H44cuQIHnjgAfj7+0OlUmH37t213ubQoUPo3r071Go12rRpg+Tk5Ept6vv6UhIlrtuyZcvQs2dPuLm5wdvbGzExMcjJyTFqU9Vz88knnzRqc+HCBURHR6Np06bw9vbG7NmzUVJSYslVqWTx4sWV6g4NDRWX3759G/Hx8WjevDlcXV0xcuTISrPVy3G9ACA4OLjK94v4+HgAytpntb3X1OW9Oy8vD7GxsdBoNHB3d0dcXBxu3Lhh1OaHH35Av3794OTkhICAAKxYsaL+xZr81/iszI4dOwRHR0fhnXfeEbKysoQpU6YI7u7uQm5ubpXtv/76a8He3l5YsWKF8OOPPwoLFiwQHBwchFOnToltli9fLmi1WmH37t3C999/Lzz44INCSEiI8Ndff8liHcaNGyds2LBBOHnypHD69Glh4sSJglarFS5duiS2mTBhgjB06FDh999/Fy95eXlmqb8h67BlyxZBo9EY1afX643ayH0//Pnnn0b1Z2ZmCvb29kY/wm3J/bB3715h/vz5wn//+18BgPDRRx/V2P6XX34RmjZtKsycOVP48ccfhXXr1gn29vZCSkqK2Ka+20RJlLpuUVFRwpYtW4TMzEzhu+++E4YPHy4EBgYKN27cENsMGDBAmDJlitHzrqCgQFxeUlIidOrUSYiMjBROnjwp7N27V/D09BTmzZsnxSqJFi1aJHTs2NGo7qtXr4rLn3zySSEgIEDYv3+/cPz4caF3795Cnz59xOVyXS9BEIQrV64YrVdqaqoAQDh48KAgCMraZ7W919TlvXvo0KFCly5dhG+++Ub46quvhDZt2giPPvqouLygoEDw8fERYmNjhczMTOG9994TnJ2dhTfeeKNetTJE1aJXr15CfHy8+Hdpaang7+8vLFu2rMr2jzzyiBAdHW10XUREhPDPf/5TEARBKCsrE3x9fYWVK1eKy/Pz8wW1Wi289957ZliD+q/D3UpKSgQ3Nzdh69at4nUTJkwQHnroIVOXWq36rsOWLVsErVZb7f0pcT+89tprgpubm9GHmaX3Q7m6hKg5c+YIHTt2NLpuzJgxQlRUlPh3Y7eJnFnLul25ckUAIBw+fFi8bsCAAcKzzz5b7W327t0r2NnZGX1xef311wWNRiMUFhaas9waLVq0SOjSpUuVy/Lz8wUHBwdh586d4nWnT58WAAhpaWmCIMh3vary7LPPCq1btxbKysoEQVDuPrv7vaYu790//vijAED49ttvxTaff/65oFKphN9++00QBEHYuHGj0KxZM6N1mzt3rtCuXbt61cfDeTUoKipCRkYGIiMjxevs7OwQGRmJtLS0Km+TlpZm1B4AoqKixPbnzp2DXq83aqPVahEREVHtfVp6He5269YtFBcXw8PDw+j6Q4cOwdvbG+3atcNTTz2FP//806S1l2voOty4cQNBQUEICAjAQw89hKysLHGZEvfD22+/jbFjx8LFxcXoekvth/qq7bVgim0iV9a0bgUFBQBQ6fW/bds2eHp6olOnTpg3bx5u3bolLktLS0NYWJjRrO5RUVEwGAxGr0MpnDlzBv7+/mjVqhViY2Nx4cIFAEBGRgaKi4uN9lloaCgCAwPFfSbn9aqoqKgI//nPf/DEE09ApVKJ1yt1n1VUl/futLQ0uLu7o0ePHmKbyMhI2NnZIT09XWzTv39/ODo6im2ioqKQk5ODa9eu1bke/uxLDf744w+UlpZW+nkHHx8fZGdnV3kbvV5fZXu9Xi8uL7+uujam1JB1uNvcuXPh7+9v9KQdOnQoHn74YYSEhODs2bN4/vnnMWzYMKSlpcHe3l7ydWjXrh3eeecddO7cGQUFBVi1ahX69OmDrKwstGzZUnH74dixY8jMzMTbb79tdL0l90N9VfdaMBgM+Ouvv3Dt2rVGPzflyhSvOzkoKyvD9OnTce+996JTp07i9ePGjUNQUBD8/f3xww8/YO7cucjJycF///tfANXv+/JlUomIiEBycjLatWuH33//HUuWLEG/fv2QmZkJvV4PR0dHuLu7G93m7vdvOa7X3Xbv3o38/HxMnDhRvE6p++xudXnv1uv18Pb2NlrepEkTeHh4GLUJCQmpdB/ly5o1a1anehiiqEbLly/Hjh07cOjQIaOB2WPHjhX/HxYWhs6dO6N169Y4dOgQBg8eLEWpRnQ6HXQ6nfh3nz590L59e7zxxht44YUXJKysYd5++22EhYWhV69eRtfLfT+QssXHxyMzMxP/+9//jK6fOnWq+P+wsDD4+flh8ODBOHv2LFq3bm3pMuts2LBh4v87d+6MiIgIBAUF4YMPPoCzs7OElZnW22+/jWHDhsHf31+8Tqn7TO54OK8Gnp6esLe3r3R2Rm5uLnx9fau8ja+vb43ty/+tz302RkPWodyqVauwfPly7Nu3D507d66xbatWreDp6Ymff/650TXfrTHrUM7BwQHdunUT61PSfrh58yZ27NiBuLi4Wh/HnPuhvqp7LWg0Gjg7O5tkv8qVNaxbQkIC9uzZg4MHD6Jly5Y1to2IiAAAo9dXVetevkwu3N3dcc899+Dnn3+Gr68vioqKkJ+fb9Tm7vdvua/Xr7/+ii+//BKTJ0+usZ1S91ld3rt9fX1x5coVo+UlJSXIy8sz+b5kiKqBo6MjwsPDsX//fvG6srIy7N+/36iXoyKdTmfUHgBSU1PF9iEhIfD19TVqYzAYkJ6eXu19WnodAGDFihV44YUXkJKSYnRcuTqXLl3Cn3/+CT8/P5PUXVFD16Gi0tJSnDp1SqxPKfsBuDNlRmFhIR577LFaH8ec+6G+anstmGK/ypWS100QBCQkJOCjjz7CgQMHKh3yqMp3330HAOLzTqfT4dSpU0YfZKmpqdBoNOjQoYNZ6m6IGzdu4OzZs/Dz80N4eDgcHByM9llOTg4uXLgg7jMlrNeWLVvg7e1tNBVKVZS6z+ry3q3T6ZCfn4+MjAyxzYEDB1BWViaGR51OhyNHjqC4uFhsk5qainbt2tX5UB4ATnFQmx07dghqtVpITk4WfvzxR2Hq1KmCu7u7eAbD+PHjhcTERLH9119/LTRp0kRYtWqVcPr0aWHRokVVTnHg7u4ufPzxx8IPP/wgPPTQQ2Y/tb4+67B8+XLB0dFR2LVrl9HpsNevXxcEQRCuX78uPPfcc0JaWppw7tw54csvvxS6d+8utG3bVrh9+7Ys1mHJkiXCF198IZw9e1bIyMgQxo4dKzg5OQlZWVlG6ynn/VCub9++wpgxYypdb+n9cP36deHkyZPCyZMnBQDCq6++Kpw8eVL49ddfBUEQhMTERGH8+PFi+/IpDmbPni2cPn1a2LBhQ5VTHNS0TZRMqev21FNPCVqtVjh06JDR6//WrVuCIAjCzz//LCxdulQ4fvy4cO7cOeHjjz8WWrVqJfTv31+8j/LT5YcMGSJ89913QkpKiuDl5SX5VACzZs0SDh06JJw7d074+uuvhcjISMHT01O4cuWKIAh3pjgIDAwUDhw4IBw/flzQ6XSCTqcTby/X9SpXWloqBAYGCnPnzjW6Xmn7rLb3mrq8dw8dOlTo1q2bkJ6eLvzvf/8T2rZtazTFQX5+vuDj4yOMHz9eyMzMFHbs2CE0bdqUUxyYw7p164TAwEDB0dFR6NWrl/DNN9+IywYMGCBMmDDBqP0HH3wg3HPPPYKjo6PQsWNH4bPPPjNaXlZWJixcuFDw8fER1Gq1MHjwYCEnJ0c26xAUFCQAqHRZtGiRIAiCcOvWLWHIkCGCl5eX4ODgIAQFBQlTpkwx+4dDfdZh+vTpYlsfHx9h+PDhwokTJ4zuT+77QRAEITs7WwAg7Nu3r9J9WXo/HDx4sMrnRXnNEyZMEAYMGFDpNl27dhUcHR2FVq1aGc1xVa6mbaJ0Sly3qvYxAHHfXbhwQejfv7/g4eEhqNVqoU2bNsLs2bON5hwSBEE4f/68MGzYMMHZ2Vnw9PQUZs2aJRQXF0uwRn8bM2aM4OfnJzg6OgotWrQQxowZI/z888/i8r/++kt4+umnhWbNmglNmzYV/vGPfwi///670X3Icb3KffHFFwKASu9jSttntb3X1OW9+88//xQeffRRwdXVVdBoNMKkSZPEjoBy33//vdC3b19BrVYLLVq0EJYvX17vWlWCIAh177ciIiIiIoBjooiIiIgahCGKiIiIqAEYooiIiIgagCGKiIiIqAEYooiIiIgagCGKiIiIqAEYooiIiIgagCGKiEzmyJEjeOCBB+Dv7w+VSoXdu3fX+z4EQcCqVatwzz33QK1Wo0WLFnjppZdMXywRUSM1kboAIrIeN2/eRJcuXfDEE0/g4YcfbtB9PPvss9i3bx9WrVqFsLAw5OXlIS8vz8SVEhE1HmcsJyKzUKlU+OijjxATEyNeV1hYiPnz5+O9995Dfn4+OnXqhJdffhkDBw4EAJw+fRqdO3dGZmYm2rVrJ03hRER1xMN5RGQxCQkJSEtLw44dO/DDDz9g9OjRGDp0KM6cOQMA+PTTT9GqVSvs2bMHISEhCA4OxuTJk9kTRUSyxBBFRBZx4cIFbNmyBTt37kS/fv3QunVrPPfcc+jbty+2bNkCAPjll1/w66+/YufOnXj33XeRnJyMjIwMjBo1SuLqiYgq45goIrKIU6dOobS0FPfcc4/R9YWFhWjevDkAoKysDIWFhXj33XfFdm+//TbCw8ORk5PDQ3xEJCsMUURkETdu3IC9vT0yMjJgb29vtMzV1RUA4OfnhyZNmhgFrfbt2wO405PFEEVEcsIQRUQW0a1bN5SWluLKlSvo169flW3uvfdelJSU4OzZs2jdujUA4KeffgIABAUFWaxWIqK64Nl5RGQyN27cwM8//wzgTmh69dVXMWjQIHh4eCAwMBCPPfYYvv76a7zyyivo1q0brl69iv3796Nz586Ijo5GWVkZevbsCVdXV6xevRplZWWIj4+HRqPBvn37JF47IiJjDFFEZDKHDh3CoEGDKl0/YcIEJCcno7i4GC+++CLeffdd/Pbbb/D09ETv3r2xZMkShIWFAQAuX76MadOmYd++fXBxccGwYcPwyiuvwMPDw9KrQ0RUI4YoIiIiogbgFAdEREREDcAQRURERNQADFFEREREDcAQRURERNQADFFEREREDcAQRURERNQADFFEREREDcAQRURERNQADFFEREREDcAQRURERNQADFFEREREDfB/v8h022z2/4gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 1024\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": [
    "## Take home message\n",
    "\n",
    "- About parallel matrix-by-vector product and different blocking.\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"
   ]
  },
  {
   "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 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "nav_menu": {},
  "rise": {
   "scroll": true,
   "theme": "sky",
   "transition": "zoom"
  },
  "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
}