lecture-8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#Lecture 8. Multigrid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Motivation\n",
    "Elliptic problems are typically **ill-conditioned**.\n",
    "For instance, complexity to solve $N\\times N$ linear system from the Poisson equation:\n",
    "* Sparse LU <font color='red'>$\\mathcal{O}(N^{3/2})$ </font>\n",
    "* Direct FFT solver  <font color='red'>$\\mathcal{O}(N \\log N)$</font> (works only in rectangular domains)\n",
    "* Iterative methods <font color='red'>$\\mathcal{O}(Nk)$</font>, k is number of iterations\n",
    " * $k=\\mathcal{O}(N^2)$ - Jacobi or Gauss-Seidel\n",
    " * $k = \\mathcal{O}(N)$ - Kyrlov subspace methods (cg, gmres)\n",
    "* **Multigrid method** <font color='red'>$\\mathcal{O}(N)$</font>. Optimal complexity\n",
    "\n",
    "The trick is that multigrid uses additional information about the grid structure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Idea\n",
    "Iterative methods <font color='red'> smooth error </font> (works for **any method**, inculinding Jacobi, Gauss-Seidel or ILU).\n",
    "\n",
    "As a result, error can be approximated on a coarser grid and so on **recursively**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A bit of history\n",
    "* Was initially proposed by Fedorenko (1961) (многосеточный метод)\n",
    "* Became popular after works of Brandt (1973, 1977)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence issues and  smoothing property of iterative methods\n",
    "\n",
    "In the currect section we will show that iterative methods does not have optimal complexity, but have some nice property (smoothing) that will help us to build an $\\mathcal{O}(N)$ method (multigrid)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#####  Model problem \n",
    "\n",
    "Consider a discretized 1D diffusion equation with zero boundary conditions\n",
    "\n",
    "$$A_h u_h = f_h,$$ \n",
    "\n",
    "where\n",
    "$$\n",
    "A_h = \\frac{1}{h^2} \n",
    "\\begin{pmatrix}\n",
    "2 & -1 \\\\\n",
    "-1 & 2 & -1 \\\\\n",
    "& -1 & \\ddots & \\ddots \\\\\n",
    "& & \\ddots & \\ddots & -1 \\\\\n",
    "& & & -1 & 2\n",
    "\\end{pmatrix}\\quad \\text{of size $N\\times N$, $\\ $ $ h = N^{-1}$}\n",
    "$$\n",
    "As we know so far $$\\lambda_k(A_h) = \\frac{4}{h^2} \\sin^2 \\frac{\\pi kh}{2}>0, \\quad \\text{cond}(A_h) = \\frac{\\lambda_\\max}{\\lambda_\\min} = \\mathcal{O}\\left(\\frac{1}{h^2}\\right)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Iterative process\n",
    "Consider Jacobi iterative process with relaxation parameter $\\tau$ (fixed-point iteration with diagonal as a perconditioner):\n",
    "\n",
    "$$\n",
    "u^{(k+1)} = u^{(k)} - \\tau \\ D^{-1}(A_h u^{(k)} - f_h), \\quad D = \\text{diag}(A_h)\\equiv \\frac{2}{h^2} I. \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Convergence. Bad news\n",
    "Let us show that the number of iterations required to get accuracy $\\epsilon$ is $\\mathcal{O}(N^2 \\ln \\epsilon)$.\n",
    "\n",
    "For the error $e^{(k)} = u^{(k)} - u_h$ we have\n",
    "$$\n",
    "e^{(k+1)} = \\left(I - \\frac{\\tau h^2}{2}A_h \\right) e^{(k)}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "It is known that (**show why**) the optimal tau is $\\frac{\\tau_{\\text{opt}}h^2}{2} = \\frac{2}{\\lambda_{\\min} + \\lambda_{\\max}}$.\n",
    "\n",
    "Therefore,\n",
    "$ \\|e^{(k)}\\|_2 \\leqslant \\left(\\frac{1 - \\lambda_\\min/\\lambda_\\max}{1 + \\lambda_\\min/\\lambda_\\max}\\right)^{k}\\|e^{(0)}\\|_2 = \\left(\\frac{1 - \\mathcal{O}(N^{-2})}{1 + \\mathcal{O}(N^{-2})}\\right)^{k}\\|e^{(0)}\\|_2\\leqslant \\epsilon \\quad \\Longrightarrow$ <font color='red'>$\\quad k = \\mathcal{O}(N^2 \\ln \\epsilon)$</font>."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For conjugate gradients $ \\|e^{(k)}\\|_2 \\leqslant \\left(\\frac{1 - \\sqrt{\\lambda_\\min/\\lambda_\\max}}{1 + \\sqrt{\\lambda_\\min/\\lambda_\\max}}\\right)^{k}\\|e^{(0)}\\|_2 \\quad \\Longrightarrow$ <font color='red'>$\\quad k = \\mathcal{O}(N \\ln \\epsilon)$</font>.\n",
    "\n",
    "The overall algorithm consists of $k$ matvecs $\\quad \\Longrightarrow \\quad$ the total complexity is $\\mathcal{O}(N k)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Convergence example**\n",
    "\n",
    "The following example shows that even for $N = 20$ Jacobi method requires $1000$ iterations to reach $\\epsilon=10^{-4}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false,
    "scrolled": true,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x10fa550d0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAZEAAAESCAYAAAA8BeghAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAHStJREFUeJzt3XmYXFWZx/FvpZOQsAQIsiUEimEzLkBEWWQrRSCoDCK4\n",
       "oIkQx3FlABklCTjQLsOg4DCSkS0YZVADuAKOk8A4BEGUPRAIgSQIJAQC6DjqPAMi/uaPc5qu2+mt\n",
       "uqvr3Fv393meelL3dtWtt1/o99S559xzwczMzMzMzMzMzMzMzMzMzMzMzMzMzMzMzMzMzCyndgau\n",
       "AL6XOhAzMysuNyJmZjk1qsWftwBYDyzrsX86sAJYCcxucUxmZlYQBwPTyDYiHcAqoAqMAZYCU+t+\n",
       "7p6ImVlOtboncivw3z327UtoRB4HXgKuBo4BJgKXAnvj3omZWS6NTh0AMBlYU7e9FtgP+C3w8SQR\n",
       "mZnZoOShEdEw3rsK2KVZgZiZlcRqYNfUQQxVleyYyP7AorrtuQz+9JWATqDWhLiKrjN1ADnSmTqA\n",
       "HOlMHUCOdKYOIAdqhDwM58t7RqvHRHpzN7AboXEZC7wPuL6B93cCS5odlJlZG1pCkxvTVjciC4Hb\n",
       "gd0J4yCzgD8DJwOLgeXANcDDDRyzE/dEzMwGo4Z7ZBlN65K1gVrqAHKkljqAHKmlDiBHaqkDyBHX\n",
       "zshjImZmg1ejyWMiRedEmJk1rq0G1s3MrKDaoRHpxKezzMwGo4YH1jN8OsvMrHE+nWVmZum1QyPS\n",
       "iU9nmZkNRg2fzsoQaPPUQZiZFYyHAiKBFoMmpA7EzKxAPCbS7cSJsOgh0L6pIzEzy7kaPp2VEVtT\n",
       "vR/0NOhi0BZpQzIzyz2fzorqEqEtQZeCngF9CjQmXVhmZrnmRiTqJRHaC3Qj6BHQu0CV1odlZpZr\n",
       "bkSifhKhI0HLQLeC9mtdSGZmuedGJBpgFV91gD4Megp0Daiwt4M0M2uCGl7FN2OQidAmoLNAz4Pm\n",
       "gbYZ2bDMzHLNjUjUYCK0NehroN+AzgZtOjJhmZnlmhuRaIiJ0F+BvgtaB/qEZ3KZWcm4EYmGmQjt\n",
       "A/rPOJPrOM/kMrOScCMSNSERqoCOAC0F/Qp0yPCPaWaWa25EoiYmQqNAM0CPg24AvbZ5xzYzyxU3\n",
       "ItEAU3yHdMhxoE+DngV9A7RD845tZpZUDU/xzRjBRGgL0D/FmVzneU0uM2sjbkSiFiRCO8QeybOg\n",
       "00EbjfxnmpmNKDciUQsTodeCro9jJjPDGIqZWSG5EYkSJEIHg34ZZ3Md6WnBZlZApWhENgGuBC4H\n",
       "PtDHaxIlQhXQu+P1JT8L15uYmRVGKRqRmcA74vOr+3hN4kRoDOjj8cr3heFKeDOz3Cvs7XEXAOuB\n",
       "ZT32TwdWACuB2XHfZGBNfP5yS6JrWOUlqFwK7A4sB+6Ma3NtnTgwM7O2dDAwjWwj0gGsAqrAGGAp\n",
       "MBWYQXdPZGEfx8tZl0zbxFWCnwd9LqwebGaWOzmrnY2pkm1EDgAW1W3PiY+NCT2Xi4ET+jhWThOh\n",
       "XUFXx/uYfBQ0OnVEZmZ1clo7B6dKthE5Hphftz0DmDfIY+U8EXoT6L9AK0DHeiaXmeVE02pnHr4h\n",
       "D/eX6ax7viQ+cqJyF+gwwpjPl4HPgM6Ayi8SB2Zm5VKjqctDpVUl2xPZn+zprLl0D64PJOc9kXrq\n",
       "AH0I9AToOtDU1BGZWWkVqHZuqEq2ERkNrI77x9I9sD4YI7AA40jTONBnQM+B5oMmp47IzEqjRpMX\n",
       "YGz1FN+FwO2EKbFrgFnAn4GTgcWEabLXAA+3OK4WqrwAlQsIOfgt8ADoXNDmiQMzMyudQnfJAk0B\n",
       "fRO0HnSaF3g0sxZog9rZHG2UCL0e9BPQY6APeIFHMxtBbVQ7h6eAYyID0aGgO0D3gg5PHY2ZtZUa\n",
       "vilVRpsmQhXQe0ArQTeCpqWOyMzaSpvWzsa1YU+knsaAPgl6GvQd0M6pIzKzQqvhnkhGSRKhzUDn\n",
       "xFv1Xgh6VeqIzKzQSlI7B1ayRGhb0NfjAo9ngjZOHZGZFVLJamff2vx0Vl+0G+ha0FrQR7zAo5kN\n",
       "Ug2fzsooeSK0L2gJaDnoGC/waGaDVPLa2c2JCDO53gFaBroN9ObUEZlZ7rl2Rk7EK9QBOgn0JOiH\n",
       "oFenjsjMcsu1MyrpmEh/ND4sN6/nQJeBtk8dkZnlRg2PiWQ4EX3SRND5cVrwF0ETUkdkZrnh2hk5\n",
       "EQPSTqAr4wKPp4DGpo7IzJJz7YyciEHTnqCfglaD3u8FHs1KzbUz8phIw/QW0F2gu+Ote82sPGp4\n",
       "TCTDiRgSjQK9F7QKtAi0V+qIzKylXDsjJ2JYNBZ0MugZ0FVh/MTMSsC1M3IimkITQF+IM7m+Ctoq\n",
       "dURmNqJcOyMnoqm0HeiSuMDjnHDNiZm1IdfOyIkYEdoD9H3QGtCHw9XwZtZGXDsjJ2JEaX/Qz0EP\n",
       "gt7pBR7N2oZrZ+QpviNOFdDRoIdAt4SGxcwKqoan+GY4ES2jjnhqa0081bV76ojMbMhcOyMnouW0\n",
       "cRx0fz4Owm+XOiIza5hrZ+REJKOt4nTg34A+D9osdURmNmiunZETkZyq8ULFZ+KFi17g0Sz/SlE7\n",
       "dwauAL7Xz2tKkYhi0N5xCZVVcUkVz+Qyy69S1U43IoWiw0D3xEUe35I6GjPrVaFq5wJgPbCsx/7p\n",
       "wApgJTC7n/e7ESkcjQKdAHosLj+/Z+qIzCyjULXzYGAa2UakA1gFVIExwFJgKjATuBCYVPdaNyKF\n",
       "pbHxRljrQd8C7Zg6IjMDClg7q2QbkQOARXXbc+Kj3kTgUvrvqRQuEeWkCaAvxZlc54Mmpo7IrOQK\n",
       "VzurZBuR44H5ddszgHlDOG7hElFumgS6DPQc6Awv8GiWTNNq5+hmHahBzSz+nXXPl8SH5VJlHfAx\n",
       "0IXAPwEng84GroLKy2ljM2trNQq+PFSVbE9kf7Kns+bS/+B6X9wTKTS9GXQb6AHQ2z0t2KxlClc7\n",
       "q2QbkdHA6rh/LN0D643yAoyFpwroGNDDoJtB+6aOyKyN1SjgAowLgXXAi8AaYFbcfxTwCGGW1twh\n",
       "HrtQibD+aDToI6C1oGtBu6aOyKyNuXZG7om0HW0MOjMu8PivoG1TR2TWRmoUsCcykpyItqVXhQF4\n",
       "PQ86B7Rp6ojM2ohrZ+SeSNvTzqDvgJ4GfRI0JnVEZgVWwz2RDCeiNPQG0E2gR0HHeyaX2bC4dkbu\n",
       "iZSODgfdC/oV6JDU0ZgVTA33RDKciFLSKNAHQb8G/QT0utQRmRWMa2fkRJSaNgKdFhd4XADaIXVE\n",
       "ZgXh2hn5dJYB2hz0j3GBxy+DtkwdkVlO1fDprAwnwupoMmg+6FnQ34PGpY7ILKdcOyMnwnqh14Cu\n",
       "Az0B+hCoI3VEZjnj2hk5EdYPHQS6HXQ/aLqnBZu9wrUz8piIDUAV0LGgFaD/Ar0xdURmCdXwmEiG\n",
       "E2GDpNGgj4KeAl0N2iV1RGYJuXZGToQ1SJuAPhfX5LoItE3qiMwScO2MnAgbIm0N+lpsTP4hNC5m\n",
       "peHaGXlMxIZJu4AWgtaBPuYFHq3N1fCYSIYTYU2iN4J+BnoE9G7P5LI259oZORHWRKqAjgQtjVOD\n",
       "D0odkdkIce2MnAgbARoFmgl6HHR9uHjRrK24dkZOhI0gjQOdHpdRuSIsq2LWFlw7IyfCWkBbgM6L\n",
       "CzyeG7bNCs21M3IirIW0A+gbsWfy6bAUvVkhtax2VoAprfqwIfAUX0tArwPdEG+K9cEwhmJWCDVa\n",
       "PMW3AjzYqg8bAvdELCEdQrhN732gI1JHY9aAltbOK4F9W/mBDXAjYompAjoO9CjoJtAbUkdkNggt\n",
       "rZ2PAC8DjwHL4uOBVgbQDzcilhMaA/oE6GnQd0E7p47IrB8trZ3V+NgpPrq288CNiOWMNgWdHdfk\n",
       "+pewRpdZ7rS8du4N/B1wMrBXqz+8H25ELKe0DWhebEzO8gKPljMtrZ2nEgbXvwB8kXA665QWffYx\n",
       "wOXA1cDhvfzcjYjlnHYFXRPvY/K34b4mZsm1tHYuA+q/RW0S97XSFsAVvex3I2IFoTeBbgY9DHqX\n",
       "F3i0xFreiIyv2x5P443IAmB9L++bDqwAVgKz+3n/BYRTaj25EbECUQV0FOgB0C9AB6aOyEqrpbXz\n",
       "dMJsrE7g88D9wKcbPMbBwDSyjUgHsIowSD8GWApMBWYCFwKTCNepfBk4rI/juhGxAlIH6ETQk6Af\n",
       "g6amjshKp2W1cxRwILAPYWzkFEJjMBRVso3IAcCiuu058VHvFOBu4BLgY70c042IFZjGgT4Leg50\n",
       "OWhS6oisNJpWOwca5PsL8HXCqaR7mvWh0WRgTd32WmC/Hq+5KD7601n3fEl8mBVA5QXg/LBCMHOB\n",
       "ZaBLwr7K/6SNzdpMjYTLQ10AHE84tTQcVbI9keOA+XXbM4B5DR7TPRFrI9oR9E3QetCpXuDRRlDT\n",
       "audgFo77OHAt8CfgD/Hx+yZ89lNkF3ecQuiNNKoTL8BobaHyJFRmAW8DjgAeBp3gBR6tiWpkz96M\n",
       "uK4xkWaoku2JjAZWx/1j6R5Yb4R7ItbGVAPdCboH9LbU0VhbaWntXNqEYywE1gEvEsZBZsX9RxHW\n",
       "5lpFOCfcKC8Fb21OFdB7QCtBi0G9TXU3G6waLV4KHpo3JjIS3BOxktAY0KdAz4C+DaqmjsgKraW1\n",
       "84+EWVov0dwxkWZwT8RKRpuBOuOtev8ZtFXqiKxQaiToiXQAHwLOids7seFU3FTcE7GS0nagi+MC\n",
       "j3NBG6eOyAqlpbXzUsK1Ig/H7YnAXa0MoB/uiVjJaXfQ90BrQX/jBR5tADUS9ETu6/EvhKVP8sA9\n",
       "ETMAtB/oFtBDoL/2Ao82gJbWzjsIp7S6GpGtyTYoKbkRMXuFKqB3gB4E3Qo6IHVEllstrZ0zgOsJ\n",
       "FweeCzwKvLeVAfTDp7PMNqAO0CzQGtAPQHukjshyo0aC01kQLgI8OT7ytOKoeyJmfdJ40BlxgcdL\n",
       "Qdunjshyw7UzciLMBqSJoAvitOAvgiakjsiSc+2MnAizQdNOoH+LFyz+HWhs6ogsGdfOyGMiZg3T\n",
       "XqD/AK0Cvc8LPJZKjURjInnlRJgNmd4Kujs+3po6Gmsp187IiTAbFo2KvZHVsXeyZ+qIrCVcOyMn\n",
       "wqwpNDaOkzwDujKMn1gbc+2MnAizptIE0BfiTK4Lwswua0OunZEH1s1GhLaP15Y8B5odrjmxNlDD\n",
       "A+sZToTZiNIe8ar3J+NV8B2pI7KmcO2MnAizltABcT2uZXF9Li/wWGyunZETYdYyqoCOjisF3xJW\n",
       "DraCcu2MnAizltPoeO+StfFeJruljsga5toZORFmyWhj0Jw4+H4xaNvUEdmguXZGnp1llpy2Itzv\n",
       "/TeE+79vljoi61MNz87KcCLMckNV0LdBT4M+BRqTOiLrk2tn5ESY5Y6mgRaDVoLe45lcueTaGTkR\n",
       "Zrmlt4HuBd0JqqWOxjJcOyMnwizXNAp0Augx0L+DXp86IgNcO1/hRJgVgjYCnQpaD/omaMfUEZVc\n",
       "KWrnq4FLgGuBv+njNaVIhFn70OagL8WZXF8BbZk6opIqVe0cRWhIelOqRJi1D00CXQ56FvQZ0LjU\n",
       "EZVMoWrnAmA9sKzH/unACmAlMLuP9x4N/Afw7j5+XqhEmFlPmgr6EegJ0Ile4LFlClU7DwamkW1E\n",
       "OoBVQBUYAywFpgIzgQuBST2OcV0fxy5UIsysLzoQ9AvQA6C3e1rwiCtc7aySbUQOABbVbc+Jj3qH\n",
       "Al8DLgNO6+O4hUuEmfVFFdC7QA+Dbga9KXVEbaxptXN0sw7UoMnAmrrttUDPFUFviY+BdNY9XxIf\n",
       "ZlY4FQE/Bv0EmBWf3wacBZVVaWMrvBoFXx6qSrYnchwwv257BjBvCMd1T8SsbWkT0Fmg50HzQNuk\n",
       "jqiNNK12jmrWgRr0FDClbnsKoTcyFJ0UvIU1s95U/hcq/0iY7v8ysBx0NmjTxIEVWY3s2ZvCqJLt\n",
       "iYwGVsf9Y+keWG+UeyJmpaG/An0XtA70CS/wOCyFqp0LgXXAi4RxkFlx/1HAI4RZWnOHeGwvBW9W\n",
       "OnoD6CbQI6DjPJOrITW8FHyGE2FWWjoCdB/oV6BDUkdTMK6dkXsiZqWmUaAZoMdBN4BemzqinKvh\n",
       "nkiGE2FmxAUePx2XUfkGaIfUEeWca2fkRJhZHW0BOjcu8Hhe2LZeuHZGPp1lZr3QZNAVsWdyeuip\n",
       "GD6dtQEnwsz6odeCrotjJjPDGIrh2vkKJ8LMBkEHg34JWgo60tOCXTu7+HSWmQ2SKqB3x+tLfgba\n",
       "J3VECdTw6awMJ8LMGqQxoI/FK98XhivhS8e1M3IizGyItAnoH+JMrotAW6eOqIVcOyOfzjKzYdI2\n",
       "sRF5HvS50Li0rRo+nZXhRJhZk2gX0NWgp0AfBaW631IruHZGToSZNZneGAfeV4CObdOZXK6dkRNh\n",
       "ZiNAFdB00P2Ee78fmDqiJnPtjJwIMxtB6ogXKT4RL1ocyn2P8si1M/LAupm1gMaB/j4uozI/LKtS\n",
       "SDU8sJ7hRJhZC2lL0JfjtOBzQZunjmiIXDsjJ8LMEtAU0ALQetBpBVzg0bUzciLMLCG9Lt4M69eg\n",
       "DxZogUfXzsiJMLMc0KGgO0D3gg5PHc0guHZGToSZ5YQqoONBj4JuBE1LHVE/XDsjz84ys5zRGNAn\n",
       "QE+DvgPaOXVEdWp4dlaGE2FmOaVNQefEmVz/AnpV6ojquHZGToSZ5Zy2BX09LvB4Jmjj1BHh2vkK\n",
       "J8LMCkK7ga4FrQV9JPECj66dkRNhZgWjfUFLQMtBxyRa4NG1M3IizKyAVAG9HbQMdBvoza0OoMWf\n",
       "l8wmwF3AO/r4eWkSYWbtSB2gk0BPgn4EenWrPrhZB8r71ZVnANekDsLMbGRUXobKt4A9gNuBW0GX\n",
       "gbZPG9fgtaIRWQCsB5b12D8dWAGsBGb38r7DgeXAcyManZlZcpX/g8r5hMbk98CDoC+BJiQOLBcO\n",
       "BqaRbUQ6gFVAFRgDLAWmAjOBC4FJwJfi88XAj4HeBp98OsvM2pB2BH0L9AzoFNDYZn9Ak4834qpk\n",
       "G5EDgEV123PiozcnAm/v42eFS4SZ2eBpT9BPQatB72/iAo+Fq51Vso3I8cD8uu0ZwLwhHLdwiTAz\n",
       "a5zeAroTdDfosGYcsAnHACDVxS7NLP6ddc+XxIeZWRup3AzaD3gPcBloFTAbKvcP8gA1RmiNwVSz\n",
       "s54CptRtTwHWDuN4SwiNyZJhHMPMLMcqgsq1wGuAG4DFoKtA1UG8eQkFr5FVsqezRgOr4/6xdA+s\n",
       "N8qns8yspLQZ6PNxgcevgrZq5M0jFtYIWAisA14E1gCz4v6jgEcIs7TmDvHYXgrezEpO24EuiQs8\n",
       "zgGN7+fFNbwUfIYTYWYGgHYHfR+0BvThcDV83y9uWVg5556ImVmG9gf9HPQQ6OgeCzzWcE8kw4kw\n",
       "M9uAKqB3gh6MDcr+PV+QJKwcciLMzPqkjnhqa0081bVH1w+ShpUjPp1lZjYgjQfNhpt+Bx+6Czci\n",
       "r3AizMwGTVuF6cCunV2cCDOzxpXmfiKD0YlPZ5mZDUaN7FJRpeeeiJlZ49wTMTOz9NqhEenEp7PM\n",
       "zAajhk9nZfh0lplZ43w6y8zM0nMjYmZmQ+ZGxMzMhqwdGpFOPLBuZjYYNTywnuGBdTOzxnlg3czM\n",
       "0nMjYmZmQ+ZGxMzMhsyNiJmZDZkbETMzG7J2aEQ68RRfM7PBqOEpvhme4mtm1jhP8TUzs/TciJiZ\n",
       "2ZC5ETEzsyHLcyNSA24FLgEOTRuKmZn1Js+NyF+APwAbAWsTx2JmZoksANYDy3rsnw6sAFYCs3t5\n",
       "XyX+uw3w7T6O7dlZ3WqpA8iRWuoAcqSWOoAcqaUOIEcKNTvrm4QGo14H8K9x/2uAE4CpwEzgQmAS\n",
       "3b/k7wi9EetfLXUAOVJLHUCO1FIHkCO11AG0o9Et+IxbgWqPffsCq4DH4/bVwDHAecBVcd+xwJHA\n",
       "FsC8kQ7SzMwa14pGpDeTgTV122uB/Xq85kfxYWZmOZWqEWnW+bjVTTxWOzgndQA54lx0cy66ORfB\n",
       "6mYdKFUj8hQwpW57CkObgbVrc8IxM7M8q5KdnTWa0BJWgbHAUsLAupmZWcZCYB3wImEcZFbcfxTw\n",
       "CGGAfW6a0MzMrKwGus6k3UwBbgYeAh4ETon7JwI3AY8CNxJms3WZS8jPCuCIlkXaGh3AfcANcbus\n",
       "eYDwu34feBhYTpikUsZ8zCX8fSwDvku4NKAseejteryh/O77xGOsBL42gvEm10HowVSBMZTjdNh2\n",
       "wN7x+aaEXtxU4CvAGXH/bMI0aQjX3ywl5KdKyFeeVyho1OnAd4Dr43ZZ8wBwJfDh+Hw0sDnly0cV\n",
       "eIzua8quAU6kPHk4GJhGthFp5Hfvurj7TsIlGAA/ZcNr/NrGAcCiuu058VEmPwbeRvgmsW3ct13c\n",
       "hvBNo76HtgjYv2XRjawdgP8E3kJ3T6SMeYDQYDzWy/6y5WMi4YvVloSG9AbgcMqVhyrZRqTR3317\n",
       "Qm+2y/uBSwf60KK2vL1dZzI5USwpVAnfOu4g/E+yPu5fT/f/NJPIznhrpxxdCHyWsL5alzLmAWBn\n",
       "4DnCyhD3AvOBTShfPn4LfBV4kjAG+zvCqZyy5aFeo797z/1PMYicFLURKfO1IZsCPwBOJSxQWU/0\n",
       "n5t2yNs7gWcJ4yGVPl5Thjx0GQ28Abg4/vu/bNgrL0M+dgFOI3zBmkT4O5nR4zVlyENfBvrdh6yo\n",
       "jUizrjMpmjGEBuQqwuksCN8wtovPtycUWNgwRzvEfUX3ZuCvgV8TZv69lZCPsuWhy9r4uCtuf5/Q\n",
       "mDxDufLxRuB24DfAn4EfEk57ly0P9Rr5m1gb9+/QY3+75eQVZbzOpAL8G+FUTr2v0H1+cw4bDp6N\n",
       "JZzyWE3f39yL6lC6x0TKnIefA7vH552EXJQtH3sRZi2OJ/w+VwKfolx5qLLhwHqjv/sdhNl9Fdp8\n",
       "YB3Kd53JQYQxgKWEUzn3Ef4DTyQMMvc2je9MQn5WEBazbDeH0j07q8x52IvQE7mf8A18c8qZjzPo\n",
       "nuJ7JaHnXpY8dF2P9ye6r8cbyu/eNcV3FXDRiEdtZmZmZmZmZmZmZmZmZmZmZmZmZmZmZhb8Mf67\n",
       "E3BCk499Zo/tXzT5+M12EjAvdRBWTkVd9sSsax2gnYEPNPjegW4L3fPi1QMbPH6rDXdNJNcBGzL/\n",
       "z2NFdx7hXgr3ERalHAWcT7gvwv3AR+PrasCtwHWE5TEgrD92d9z+27rjjY/Huyru6+r1VOKxlwEP\n",
       "AO+tO/YS4HuEpbS/3UesS+Lx7yCstnBQ3H8S2Z7ET4BD6j77KzHGmwhLdt9CWKri6Lr3dN207FHg\n",
       "7Lr9M+Ln3UdY1rvrb/6PwAWEFRCKvgS6mVnDulYwrl8/C0KjcVZ8vhFhOZAqodD/kXD6q8uW8d/x\n",
       "hIaha7vn6shd28cRlo+oANsATxAWuKsRlh6fFH92O733Xm4mNEIQlu25KT4/iWwjcgPdjchf6F6W\n",
       "4ofx8zuAPQkNQ9f718X4x8XfZR/CenLXx9dDWOl3Zt1xj+8lRrOGDNStN8u7novmHQG8nu4COQHY\n",
       "lbCy652Ewt/lVOBd8fkUYLf4mr4cRLjtqggrot4CvAn4fXzfuvi6pYSGq7exlB/Gf++NrxnIn4DF\n",
       "8fky4AXgZULPpP79NwL/XfcZB8XX7UPobUFoLJ+Jz18mrAhtNixuRKwdnUz3t/wuNcK9Nuq3DyOc\n",
       "ynmB0EsYN8BxxYaNVtd4xIt1+16m77+tF3t5zZ/Jnlquj+Oluud/ITQqXc/7+oxKXVxXsuFEAQi/\n",
       "czvfP8NaxGMiVnR/ADar214MfJLuArs7sHEv75tA+Ob+AvBqsuMCL9F7gb4VeB/h72ZrwimnOxn+\n",
       "EuKPA3vH40yh+x7XjTiccDprPHAMcBvwM0KPbOv4monAjsOM1SzDPRErqq5v0fcTvtUvJdwi9iLC\n",
       "aZ57CUX5WeBYNryz2yLg48BywiD3L+t+djlh4PwewhhC1/t+RLjR0f1x32fj8aey4bf6wXzL73rN\n",
       "bYSbbC0nDMzf089x1MtzERqzHxBuJHQV4fcH+BzhVNcoQuP4ScItZN0LMTMzMzMzMzMzMzMzMzMz\n",
       "MzMzMzMzMzMzMzMzM7Ph+3/miZisbbJ6hwAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f94ba90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg as spla\n",
    "import scipy\n",
    "from scipy.sparse import csc_matrix\n",
    "\n",
    "n = 20\n",
    "ex = np.ones(n);\n",
    "lp1 = n**2*sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "ev1, vec = spla.eigs(lp1, k=2, which='LR')\n",
    "ev2, vec = spla.eigs(lp1, k=2, which='SR')\n",
    "lam_max = ev1[0]\n",
    "lam_min = ev2[0]\n",
    "\n",
    "tau_opt = 2.0/(lam_max + lam_min)\n",
    "\n",
    "niters = 1000\n",
    "x = np.zeros(n)\n",
    "res_all = []\n",
    "for i in xrange(niters):\n",
    "    rr = lp1.dot(x) - rhs\n",
    "    x = x - tau_opt * rr\n",
    "    res_all.append(np.linalg.norm(rr))\n",
    "    \n",
    "#Convergence of an ordinary Richardson (with optimal parameter)\n",
    "plt.semilogy(res_all)\n",
    "plt.xlabel('Iteration number')\n",
    "plt.ylabel('error')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Smoothing property. Good news \n",
    "Iterative processes remove high-frequency components in the error! Let us show that.\n",
    "\n",
    "From the previous slide\n",
    "$$\n",
    "e^{(k+1)} = T e^{(k)} = T^{k} e^{(0)}, \\quad \\text{where}\\quad T = \\left(I - \\frac{\\tau h^2}{2}A_h \\right).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Let $$T \\psi_i = \\mu_i \\psi_i,$$ \n",
    "where $\\psi_i$ and $\\mu_i$, $i=1,\\dots,N$ are eigenvectors and eigenvalues of $T$.\n",
    "Since $\\{\\psi_i\\}$ is orthonormal basis (**explain why**), we have\n",
    "$$e^{(0)} = \\sum_{i=1}^N c_i \\psi_i,$$\n",
    "with some coefficients $c_i$. Therefore,\n",
    "$$\n",
    "e^{(k+1)} = T^{k} e^{(0)} = \\sum_{i} c_i\\lambda_i^k \\psi_i\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "Now the goal is to show that $\\lambda_i^k$ decays fast enough with respect to $k$ when, for instance, $i>N/2$ (high-frequecy part of the spectrum). Indeed, \n",
    "$$\n",
    "| \\lambda_i(T)| = |1 - \\frac{h^2}{4}\\lambda_k(A_h)|<({\\bf\\text{show why}})<\\frac{1}{2},\\quad \\text{when}\\quad k>N/2 \n",
    "$$\n",
    "So, the high-frequency part of the error ($k>N/2$) becomes at least <font color='red'> 2 times smaller on each iteration </font>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": [
       "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\n",
       "AAALEgAACxIB0t1+/AAAIABJREFUeJztnXe4HGX1xz+TCumNVCCXIr2EIKEpBBEMoIAIooKAImKh\n",
       "KKAIqLx0aQqCP0qkKkiRjiggEKkSUWoInUAgEBJaep/fH+edu7OzU3dnd2c35/M898m9O2XfbPnO\n",
       "me8573lBURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFKQATgJeAV4HjQ7bvBTwLPA38F/iC\n",
       "b9s04Dm7bXJdR6koiqLUTFfgNaAD6A48A2wY2Ke37/dN7f4ebwKD6jg+RVEUJQNdEraPQ0R8GrAU\n",
       "uAGJ7P3M9/3eB5gd2O7UMD5FURQlR5JEfxQw3ff3O/axIHsDU4G/A0f5HneBfwJPAYdVP0xFURQl\n",
       "D7olbHdTnud2+/N54E/A+vbx7YH3gNWA+5HcwCPZh6koiqLkQZLovwus4ft7DSTaj+IRe87BwIeI\n",
       "4APMAm5D7KKg6L8GrJNyvIqiKIrwOrBu3iftZk/cAfQgPJG7DiXffqzdH6AX0Nf+3ht4DNg15DnS\n",
       "3k2sDJhmD6BAmGYPoECYZg+gQJhmD6BAVKWdSZH+MuAI4F6kkucKxLs/3G6/DPgacBCS6J0HfMNu\n",
       "Gw7c6nue64D7qhmkoiiK0j5opF/CNHsABcI0ewAFwjR7AAXCNHsABaIq7Uyq3lEay6RmD6BATGr2\n",
       "AArEpGYPoEBMavYAlNrRSF9RFCU7dfH0FUUpBh8BA5s9CKUpfEybdTbQSF9RktHvycpL1Huvnr6i\n",
       "KIoSj4q+oijKSoSKvqIoykqEir6iKEo6DNJbLIoXgB0aM5TqUdFXFCUPpgEzkfYrHt8DHmrKaGpn\n",
       "POUdhiE5cboJ8HBdRpMjKvqKouRFF+DoZg+ijrTF2iAq+oqi5IELnAccB/SP2Gc74D/AJ8jyqdvG\n",
       "nO94pKPvHKQlu7cMqwFuRmyWOchyrJ8BTkDuNN4CdvGdZyRwJ9L191Xk7sOjJ3AB0k34XeB3SGPJ\n",
       "3sjaICOBufZ5Rtj/Yw/gGvvYC8CWvvNNC4zzpph9xyLLyM6x+90InBbzerQVWn+sKMkU/XvyJrAz\n",
       "cAsl8fLbO4OQSUYHIMHmN5AJZ2GTjtYH3kaaNgKsCaxtfzfAQkTYuyKiOg0R/a72Od/wneth4GJE\n",
       "rDcHPgB2sttOBR4Hhtifx+xjADtSae94zz0BifrPBJ4IvAZfSLFvD+TidKQd81eBxb7nDpJrnX4R\n",
       "aNmBK0oDSfE9cd18fqrCE7yNkUh+COWi/23g34FjHgcODjnXukjUvjOyNrcfg3T99fgKEo171ktf\n",
       "YAXQD1n/Yxnl63ifCVxlf38dEWWPXe3/A8I9fUN5p+CNgAW+v4OiH7XvDlSuS/IIDRJ9tXcUpW1w\n",
       "nHx+amIKcDfwC8pFaSQSvft5i/DlV18DfoII50zgL4i94vGB7/eFyLrcru9vkPW6RyJ3E/51vN+2\n",
       "j2PP+VbEtihm+n5fAKxCtI5G7TsSsZP8TKdBOQMVfUVR8uZkZE1sv6C/C4wO7Dea6JX4/oIsvzoa\n",
       "EfSzqxjHDMQ+6uN7bE1KgjsDWSDKv22G/T0sis7LlXiPyovdmjmePxYVfUVR8uZ1JDHpr+T5O7Ae\n",
       "8E2k0eP+wAbIXUGQ9RCbpCfidS8CllcxjumIhXSWPddmwHeBP9vtfwF+ScnT/zWlOvyZyLKv/Xzn\n",
       "yysSfwL5/xyBvBZ7AVvldO5EVPQVRakHpyI1+170+iHwZeBYxI45zv79UcixPRGhnoVExUOQRC32\n",
       "fMGIOO7vbyLR/AxkJb9fAw/abacDTyEVQM/Z30+3215CLgpv2DF61TtJz+1/PGrfJcA+wKGUktt3\n",
       "28dXCjSRqyjJ6PekvXmS8KQ2aPWOoqyU6PekvdgBKUnthoj9fGBYxL4q+oqyEqLfk/biMOB9pNz0\n",
       "GWC3mH1V9BVlJUS/JysvWqevKIqiVIeKvqIoykpEGtGfgJQvvYo0QQqyF/As0jzov5SmIac5VlEa\n",
       "jLs5uGOaPYp8cNcFN6q5maJURVdkSnQH0gPjGWDDwD7+vhab2v3THgvqVSoNxb0M3N80exRVEPI9\n",
       "cW8A96eNH4rSYBrq6Y9DhHsasBS4AYns/fj7WvRBJl6kPVZRGs2mlE/Lb2X6AVs0exBKa5Ek+qMo\n",
       "7zT3DuENkvYGpiJTrY/KeKyiNAjXQVY36tvskeTEyij6a1LeVbMZHEB5p8+WIkn0094+3I5YN19B\n",
       "eldkfUOM72d8xmMVJS2jEcFvF9HvC2wEbq/EPetLD+AK5K5+DpLfmxDYZ2ckvzcfaYOwZmD72YhL\n",
       "MBuIs9/eRv7fnjZNQtoZ1IsOpFWzXyuvA75Ux+eMYjzlWlkXtgH+4fv7BJITsq8jjYrSHquevtIg\n",
       "3C+Duxjc+5L3LRxhnv6b4H4E7taNH04ZvZDOmp6Q74GIv9dVcwjSY/9ryAXiHMoXHzkcuSCMtD9T\n",
       "7GNpeIjaRD8p8O1ARL9rDc9RKw2dnNUNEfEO5M0KS8auQymyH2v3T3ssqOgrDcM9AdzHwH0ied/C\n",
       "ESb6s8G9A9wfNn44iTyLrAgF8H3gUd+2Xkh/+fXs349Tvozhdyi/KPjpoCTCZyCLpCxELJ/f2302\n",
       "AO5Hmry9BOznO/5q4BLgHmAeUm24B3J38ilyJ3Gyb/+37fN5yyZuAxyCLHriEbcM5CSk+dyj9vh7\n",
       "kaAYpL/+n5G7m4/tsUND/s8Nn5G7G/AykpT1Ot0dTulK/HNk/cenkRdiq4Rjg6joKw3CvR7c88B9\n",
       "odkjqYIw0V8C7jHgXt744cQyDBFiT9QvBP4Q2Oc5SheFTyjXjS0RgQyjg3K75SGkXbJHbySXeLDd\n",
       "ZwzSrdMLOK+2z+cJc09kacSN7d+bIu0RvKKT0VTaO4dQEv2oZSAH2u2TkJL1dRGRfwjpIAqioXfa\n",
       "xx0kPxNmPeYq+t1S7PN3++PnMt/v59iftMcqSh1x9wbuACfsC7Ep8Fdg38aOqR64PZEv/RPAtwAw\n",
       "OQVQpqYkaXfE874aeMU+1hsRXj9zKAlcHyTK9m/LUmHlH++XkWULr7F/P4O0VN6P0nKEt1O6k1gM\n",
       "/Mt3/PNIpeGOwB0k5yf3QALb6+zfNyDFLHvaMbjI8oxeKftNdhtIK+XByMLuzyOBc91JI/qK0iK4\n",
       "XZAv+KaIL+zf1gOJtibTHoncfog4Pockc3tgal7qsFa6IIUci5AFQjzmUb4YCUB/xDIJ297fPpYW\n",
       "/8VuNLA1En17dAOu9e0bXK1rayR5vDFiRfdExDkNUctA+pddfN/3+0JKF7Q/Iev43gAMQKyekxDL\n",
       "qm5oGwalneiHRGY7h2xbD/lyzqI9RL8vMBec+UjVzEbNHQ4OUsGzGpKw9a90NQXY3Pd3byQXOMW3\n",
       "3T9LenPEMk5D8O7mbSRyH+j76Qv8OOYc1yPR/+qI+F5KSRuT7p6iloEMroEbxjLk7mNjJC/wZeCg\n",
       "FMfVhIq+0k54LQnCRH9T5BZ6CeDYyL+VsaIPwP9ofr3+JUgCdU/EMvFzGzI/Yh/Evz4ZsV08++da\n",
       "4BgkOh5lf7865fPORC4gHncjF/gDEaupO5Iv2MBuD7sb6oPcGSxBJpV+i5LYz0I8/XVCjoN0y0BG\n",
       "3YHthHwuuyLv5VKqWxYyEyr6SjsxALmV3hHcoHW5CfC89frn0vrRvmfvgIj+2CaOZTRSobM5pR7x\n",
       "cxEhBKlO+RpSbfMR8Fkk4elxGXAXclF+zv4el5z2R98XIjmaj4ALEFtoV3v+d5HlFs9CbBvv2GD0\n",
       "/iMk4p4D/ApZ39djgR33Y/Y5tg6cI80ykG7gd+/vYcDNSD7jRSTp+ydWArR6R8kJdwdwHwb3ucra\n",
       "dfdOcPexv78FbvCWvOgEvifu7uDeY3/fSUpRlTZF++krSgQDkKjpASotnk0p+cTzaK9I/2lgc3Cb\n",
       "OYFIaRFU9JV2oj9Sgx0QfbcvMunFmzjYDvaOz9N3PkG87c80cTxKi6Cir7QTXqT/MDAO3FXt4xsD\n",
       "U8HxkmRzaf1Om/5IH5rv6ystgoq+0iDcz4H7gzo/iY30nTlIUnA7+7hXuePRDvaOv3oHVPSVlKjo\n",
       "K41iW8CEVNXE4E4E94sZnsOL9KHc4vH7+dAe9k4w0n+ahpVtup+zbaqVFkRFX2kUA5EStSwtadcC\n",
       "Nsuwv+fpQ6XoByP9Vrd3gpH+08DY+oux6yCTn7as7/Mo9ULbMCiNYiBSg30I8LeUx/SlcrZjHP5I\n",
       "/wmkPcEAKkW/FSP9jykr0XNAOlNe7vsbZCJRHel8nv/U93kUHx8n75IejfSVRjEQmd6+C7iDk3a2\n",
       "9EG6KqbFF+k7ixHh/7rd5u9/0oqJ3EGI4tof9+/g7hF47C5w9y1/LO8fdyi4S21b5571fS79sT+D\n",
       "yBEVfaVRDER6xNxD+WzMOPqSTfT9kT6IxXM08ALlXTfbIZHbj3J7BzotnrrSG5nl+iLSYVJpMVT0\n",
       "C0fbJsgGIrep1yAWTxr6kM3e8Xv6IKK/EeXWDrSmvROkL5U959+k/utQ90aWPLwG6VmvtBgq+oXC\n",
       "PZbk5ShblYFIP5J/AiPB3STFMX2AHtaXT0Mw0n8auQiEiX4Ge8c9Fdzgmq7NJpjIBbkIBFsY500f\n",
       "5E7pr8BO4K5W5+dTckZFv1iMRCpW2hEb6TvLka6KCVGi29P+8hrpLZ5ApO8sB85EIn4/Gewdtwti\n",
       "EX0h5RgaRbBkE+QiUG/Rt5G+MwfpJPnNhP2VgqGiXyy8dgFthusgUbgnyNcABybU7HsR5Vuksnhc\n",
       "u+Scs6j8cedccF4P7JzF3lkLEdKiTXxqVqTv2TuQ2uJxx4J7df2GpGRBRb9YtKno0xdYCM5S+dN5\n",
       "CUnqxtXse6I2jXSRftDPjyNLnf4WSPvcAom+dxfkBPvW+5cgrBd+0X8AGJbCqtuQbPMtlDqiol8s\n",
       "+iITmFoQ91zb2CwML4nr52rgOzEn9CL9aaQT/aCfH0eWSH8s0uO8SF0sw6J8aKynj7XP/kxytD+M\n",
       "nMsOlepR0S8WLRrpuwOQxSOiVhcaRKXoP0h89JzR3skU6WcV/QeBDyhOF8swPx8a6ul3ksaqG46K\n",
       "fmFQ0S8WfeWnsztkq7CN/TeqkiMs0v+EeOHNau9kifRT2juug4j+0xSroVlUpD8P6GWTz/UiIPrO\n",
       "VKRuf9uYY4Yjn+vudRyXkhIV/WLhiWCrlcF53SyHRGwPE/05QP+YeQlZ7Z0skf4CoGeK5m8j7b/v\n",
       "UizRj4j0nRWIINdztnEw0geYSvx75FmWaUtvlTqSRvQnAC8BrxJeQ34A8CzSV+UxyhM20+zjTwOT\n",
       "axnoSkJfYAatZ/Fsh7zXGSJ9ZzHSS6ZnyP5QimY/lH3cJNsiQ6TvuMgFpXfCjmOB/9n9iyT6UZE+\n",
       "1N/X70Ol6M+gdIEMY7j9Vy2eApAk+l2BixHh3wipyd0wsM8bwA6I2J9G+YLGLjAeqYAYV/tw256+\n",
       "yOpOLST6blfkvb2LbPYOxAuUjfQdF7mgJPn6/Ulv70C6Wn3P2oGGdbFMRVgLBo96+/q96UzkdpJG\n",
       "9N9GPgdKk0kS/XHI5JhpwFLgBmCvwD5PUPqyPQmsHthehC9JC+B2BVZBptK3UgXPJsiX/mWy2TuQ\n",
       "SvSBdBaPfx5AGtIkc7dAInzA+cAes3aG56gXYS0YPOod6YfZOzGi73ZFIvxX0Ei/ECSJ/ihguu/v\n",
       "d4jv7XEo0lDLw0Wm3T8FHFbNAFcivNvmmbRUpM92wOPALOIj/Y9CHo8TKL+FkaaCJ2ukn6YVg7V3\n",
       "OimKxZNk79SzVj9K9EdE7L8acsH/ABX9QpCUyHITtvvZCfgusL3vse2RzP5qwP1IbuCRkGON7/dJ\n",
       "9mdlw/sif0D0F6iIbIusSZsk+tVE+l7AMY38I/0Ee8cdglxI3vQ96In+zRmepx5ElWxC8Tz94Uhb\n",
       "649Qe6dWxtufmkgS/XeBNXx/r4FE+0E2AyYi3r//y/2e/XcWcBtiFyWJ/sqKX/Q3b/JYsrAd8Bsk\n",
       "/5OnveOPZqeRnBOqJtKPi4i3AJ62FTEe/wOOyPAc9cJL+IfRDE//PaSJnhNoYQ1iVb6PvP8ZIn13\n",
       "FWAFOEuqH2rbMYnygPjkak6SZO88hUxI6QB6APsDdwb2WRO4FTgQ8f89elH6UvUGdqWy26FSwi/6\n",
       "LWLvuMOAwcgdXD0ifU9c0tg71UT6cfZO0NqBzki/6cncZkb6IfaOswBYSHgkPxyxLD8im71zKsW4\n",
       "wLYdSZH+MuSFvxeJ5K5AanIPt9svA36NvNmX2MeWIlHZcORi4D3PdcB9eQ28DfFEv5U8/W2Bf0s0\n",
       "7H4IDJKJQU5wyb5GJHLzjvTHIl0k/cxAliNcnfJcV6MpmqcPJYsnmLvx2ztZ8iGrAzqZqw6kWSP3\n",
       "7/bHz2W+379nf4K8AYypclwrI/5Iv1Wqd7ZFkrhIMzV3LiLwHwb2C2vDACLSaeydD4De4PYGJ0xw\n",
       "IP/qnbFItOnDccH1fP1min5cpD8XufuqF2GePpRE/4XA48MQmzijvdMygU/LoTNyi4MnctYmqetU\n",
       "+rzwKnc8QiwetwvRs2XtrNxQ/I29XJItnmrq9CPsHbcfImAvh2wsQgVPMydnhXn6EJ3MrdbeWQ0V\n",
       "/rrQCsJSUNwTwN0ixxPaL7KzBPlSFXzKutsDSXb6Z1qH+fqBtsplpE3kgoh+R8RYuhBfux5GXKQ/\n",
       "BngenGUh24og+k3y9N0uyFyShSEb40S/muqdobTOHW9LoaJfPV+gcnZyLQTtjKJHOWOA18DxC/Ns\n",
       "Kit4ovx8SO/pQ7yv3xdZzWl5zHiDxNXphyVxPYog+s3y9HshF/BgzgaiRb+K6h23CxI8qOjXARX9\n",
       "UNz+4P4sYae+yJcgL1pN9IPWDoRH+nmKfpS9k9XPh/g6fd9M3AqmIZ0sbT8Zd1NwLwXXZHz+WmhW\n",
       "G4aoJC50lm1W4EX6VvRTVT4NRJriDSzQGgZtg4p+OFsAP07Ypy/JDbuy4Bf9mRQ/yqmj6LsOld5x\n",
       "jL2T2c+HeHvH33MnQGfztZ+D+whS5NAH2CXj89dCs9owBC/EfkIifbeHHctHtsHeYtJ9Z4YiF5FP\n",
       "iJ77oVSJin44o0m+Ra6n6LdbpB/WggGiBaoXsDhg10wjWvSrjfRD7B23BzI3JViF4uceZALdBcga\n",
       "uqfSsPfL7YmsBRxcKtGjnqIfF+mH2TtDgQ98dlBai2c1Oa4Vgh93VXBbanF4Ff1wRgP9Em5F+7DS\n",
       "2juuV0P9ZmBDXp5+mGc9jWh7J89IfyQwK0ZUAed8cHYG5xaboG6kOCUlrOvp6SfZO8MDVWeeteOR\n",
       "toLHXixaYs7KpsDVKVp/FwYV/XBGI3MYInq9uw4rd6Tv7zPvJy9PP8xGmAkMiFhVrJpIP0r0RyF1\n",
       "5VmYA/Ro0IpncUlcaJqn7yxGXgf/Rd8r1/RIW8HjiX4rzFkZjHQr2K3ZA0mLin44HfbfqIipJ3JR\n",
       "yFP0/cm5oou+7UtTQR1F31mB9GQPi/arifSj6vSrEH3HpXHvWVy5JsAioKu1gfImztOHSovHq9zx\n",
       "SGvv+CP9oov+IGAJ8NVmDyQtKvrhjEam20eJvvf4SmrvRIp+VntnIdDd+uh+oqLZaYSLfrMjfWic\n",
       "QCVE+o5L/SyeOHsHKkW/FntnFq1h7wwG7gAm1OlCmzsq+hW4XZC+H68RP3EIVt7qnYRIvywXEtWC\n",
       "AStQYa0YoiLKqAqeaiL9+Uhrh+B3oFrRL0qkD8US/WrsHS+R2yr2zovAFGTuTuFR0a9kOCIgs9BI\n",
       "PwR3MBJZv1G5zZmPrMHgvxjGRfoQbvH0ITrSXyvk8SoifWc5cqcRfA9bPNIH6ufrZxX9lcHeGYz0\n",
       "mroN2LvJY0mFin4lo5GIMq6Ouy+V4lYr/i/zJ8gEoCLeLo4Bno2YlQmVFk81ot+X8Ej/FWD9kMer\n",
       "ifQh/D0uSKTv9gH3jJANaSP9eoh+VLM1j+AKWrVW7xQ4+OlkMPL/uh3YqxUmk6noV9KBRJRxt8h9\n",
       "EXHLSfS92mtsmaDjEt+fvplEWTsewXFXG+mHif5UwltfVOPpQ3grhqJE+p9BJoAF2wunifTrJfpR\n",
       "zdY8wiL9Wqp3WijSd15DPvvbNHk8iajoV+KP9OM8/ZnkZ+94zdb8JZBFjXIaIfpRwvYqMDok8Vtt\n",
       "pB9oxeA6iGhVG+nnKVAjkAqxoJ3Vap5+RnvH7Y68nx/R+R1o+qI1cXj2Dki0X3iLR0W/kjT2Th/k\n",
       "w5yXvRMmcq0q+j57J7atskeGSN9ZjJRtrhvYUEuk73+PByEzgeOELYq8K0088QzaWS3i6burIqXN\n",
       "/otxGnvHi5xXgLMQufuNar9dBPyifxvw1YJfpFT0Q0jr6ddb9At4a+v2QuyvF2N28kf6fYEFEW2V\n",
       "PbIkckEsno0Cj9US6fvtnVFErz2bRD0ifagU/bhmax7N8vRnItVbXem0dsruXtPYO5614z9nwb4H\n",
       "ZfhF/2lkpvrGzRtOMir6lXiin+Tp18HeKaOIkf5mwNQEEfeLfpK1A9kSuRDu6+cV6Vdr7UD+kf4I\n",
       "5OIaFuk3K5Gb4Ok7SxFhH0qltQPpIv0w0S/a98Di9kDWF7Dvh+MiFk+hJ2qp6JfhOkgkm8bT/xDo\n",
       "Bm6aJSeTaBXR3wJ4JmGf2dQu+nEzPwOi766CfI4XJTxPGEHRrzaJC/L/HpjT5wHkAjQJWC/weNpE\n",
       "bjM8fSi1WA4T/blIVVrc2rdB0S9yrf4gpIOo/27GWjzFRUW/nMHAEnA+JdnemYt8AfKI9ltJ9OP8\n",
       "fJBI3yvZrCXSj7N3/JG+zRlU9AFKQ5i9U6XoO8uR/2te69OOAB4i3N5JivSb5elDydcPVu5g36OP\n",
       "iV8VzpuN61Fke8dv7Xg8CnwG3MKufKeiX45n7UA60V9APr5+u4l+lkj/UyoTdXGR/kvA+r566Gr9\n",
       "fMg30od8BWoE8F+gjyzq00kzSzaTPH0oiX5YpA/JFo83G9ejwPZOmOg7yxCrsbBdN1X0y/GLfpKn\n",
       "70X6K4nou16C6tmEHf2iH9OCoZOMiVxnLmKleD14svn5hg4Ma9u/gnX6o1j/jiUYpmDKJhkFz3E2\n",
       "hmNDtogVYXAwtSyl6XZBRHMGlRPSmjk5K6lOH5JFP6lss5XsnbBIH/LThbqgol/OaGRiFiR7+vW2\n",
       "d4p2W7sBMB2cpC+9f0ZuPRK5UG7xlEf6JmalJUMX4CbgdgzdqFwycRRf/uG+9rFTI86xCfBd4ERM\n",
       "xeQ5LyrdB5iCianiMHTBEFXaNxiZt7EYeJly0S+6px9j7wDJFTytVL3jzcYN0vKiPwG5pX4VOD5k\n",
       "+wFI9Pcc8BhS4ZH22KLRQXHsnVkUa2LKGJKtHRCR72PvDOqRyIVy0S9F+obNgPcxkUsXeisczQSO\n",
       "JPger/74aHrP3BzYFtgTw+ZlR4tInw+cBlwPnBg4/wesOnsUcB7SedGEjsLgsLD//cwbGtZmAcTa\n",
       "ec/+/grlydxW8fSrtXdaqHonNtIPa9tdCJJEvytwMSLeGyFfmuBt6xvADojYnwZcnuHYopHW0/fs\n",
       "hzraO84ipCFYUSampPHzsT15PkKi/XwTuYZRGNYjLNIXQf4tcDPwR0zgnIZewFnAMcARwEms9c/u\n",
       "eF/OAW/2ZNefDcDt8ksM7yKR/vmBaHw3JDC4BDgdOAhT1up5JhOO2QPx4g8EPodhTMj/eW+WrzKO\n",
       "nnOPxoQmNUcCMzAMpizST1wq0aMO9o7bDalBT6qSysPe8SdyW9HemUcLR/rjkBbD04ClwA3AXoF9\n",
       "nqB0e/0k0pY47bFFoxpPv172DhTL108p+kDJ4qlC9N1uyEpEC8v2EvG9DniQte+fQWWk/2VEbL4N\n",
       "3IdE235+CjyJ4VEMLwOXsdf3DsR7j/fb97us8ukyui672u5/OZLY3d0+f3ckyj8Ow1IMM4E/AKd0\n",
       "PsMajy5hg9u2Bn6GYT5wdtl2OU9v4ALuvHwab2//DvCTkNdkBNud1xN4m41unkXJ3klTow/18fRt\n",
       "lJ9YJTUDed2qtXfCErmtJvotbe+MAqb7/n7HPhbFocii0dUcWwT8or8QWQIvrO66EfYOFEb0XYf0\n",
       "9g6UkrnVRPo2WVghLl+z5/sL39zzaJzlG9lx9af7vHmIIB+DYRlwLPAlDLsCYBiORPi/8J3vDPq8\n",
       "vwHr3z4KQ0+GTvklD5z+Gobl9pilwHHAeVbwv49U9tztO8f5wG7W54e9v7MHz3/rHUzn2sGXAlti\n",
       "2Mp3zEnAo7yy5yDuPX82cAQmIIKDXl2bHU/dEniCvQ/ZAykB7EI6Px8671JztQbTWDsgn9lBgBuR\n",
       "/4mxdzpbN/gvbHOQhXbybGOeF4NoQ9HPUvu8E5Lg8rz7LMca38/4DMfliNsP+cDNlr8dl8pEH/aL\n",
       "1AB7ByiOn9mBRHmzkna0ZBH9+cAqHD5mMIb/8e1ddiLo5xtWRSL3nwAn0H3RKuxwRjfEQujPVw7f\n",
       "BHgdwz/s/nOAw4CJGPojtuOVGF73nXMB03Y8gy8duxZwDPOHvc/Le08NjO0eJFg5Hvg1cCzG97k2\n",
       "fIpE86dj2Jb+b2/M/efM8W1fBJyBlxQ2bAAcxowtTwRG8MFmI4E7kbuQEl/5/n7MGPsUcDA9Fnyb\n",
       "Pu99CqxBsAWD4RgM21W+pM4ypGdNnkKZUvSd5cjnNszagXh7x0b5/gt+Q5eizEqjI/3xlGtlVSSJ\n",
       "/rvIh81jDeRLEGQzYCKwJ6Uvedpjofw/MilhTPXCRvllEWaYr78qsNR+sRZQf3unCLe2Y0kf5UO5\n",
       "vRNW3eDDcYG5DHrtDGA+HZMmMuTFoGd9DPAUhodsJP8Ntr6wB7se93X6TR/KhrfuAoESSsN9wL1I\n",
       "tc6eiPiW8+d77+KjdVcAp3LPRfcRrNEXgT8GsWjuxISWq/4f8vpcz4frn8ni/sGKniuADTFsD1wE\n",
       "nMHlT3W3zzWST1f/DfAjjBVCwwRWm7oGd1x1hc0tXMfOJy1Fkrkle8fwJeRidkFEFVDeFk+aGn2P\n",
       "GYRbOxBv7wT9fI+iWjyNFv1JNED0n0L6encgPuv+SGTiZ03gViRx9VrGY4uE39rxCBN9fzlhvSP9\n",
       "okQ444HmxsESAAAgAElEQVRHMuyfJdKH4U8vpPuC/YB9eGfriRy4+0ibxJTkrQjvzzr3N7zDA2f+\n",
       "ky0vM+x52PZ8uN4jmNAmcMchfvgpmNBa/rncOXEJcAivfMUT4nIMzyMVaieEjl2i+ROAN7j1TxMJ\n",
       "VlwZliDifBMiXBcj34lXgRn8bvpypF/LTzH0AS7lbxdP45O1PIvobDa+cSjDn/4sXqQvpaJXIhez\n",
       "VYE9QkaWt+intXdARD8q0o+r3glW7ngU5XsQpC09/WVIpcO9SPOnG5HKicPtD8ht70CkouFpYHLC\n",
       "sUUlTPSTKktWlkTuF4F/Ztjflpv62iobDsZwDYbyXvgGhz1+1Jf3x1yGYRZXPXo/r+/yAVJLvwrw\n",
       "G+BSn08u/Pfw+3h1j9dZ47ER/OOCywlDbJ6NkM9mGPOYs2ZvjHs9cbNxDTdgPNsvdPt1wBeZOWYB\n",
       "4a2Ar0Wq3H5o71Q6kAKHN4B1kLuQHyEXhH8xdd9eeCWbhneZsdVkdjh9P6AvrJiD3FVfj+EB4GTg\n",
       "1JBofy6bX7s2JmAdVeBeB+4X4/cB0k3M8ogT/Th7J0r0Cxjpuw6dvXcqCLb3KBRhScogf7c/fi7z\n",
       "/f49+5P22KLin5jlERXpewK9gHQrASURJ/o75nD+GnBHIV/GpEZrfmYDayNtlZfZyPRc4AXgrxj2\n",
       "w3irhLEvvWY7TJx8r9WCPtx96bNs+ceFSBXOOoQvkTiVW/80h/5vv8zH6wQv1iUMC6KH6SwFdxnS\n",
       "KbG2Fgwlr98TqE9825YCn/ft3YEEGMuAtTHch+FWYC+W9N7QjuO9zr0f+/lE9v3GlYx4ahQb3j4S\n",
       "EZT97dbbgV8hlXG3dx7Tff48Jhz9W2AtDJMxPBYx8jHQOUM5jiyR/mMQOfEszt4JVu54FFD06Yv0\n",
       "6Qorn23pSH9lIou944/0a3xzyxLDQd6j+RVPOwMPxqyJG8YsxIP2rJ3fAH8CvoREwrdjWNWWL57P\n",
       "/ec+z6KB3uvYF7frPOBgpILqp5iIpRNXdN+Qj9fpTfW9d6DUiqHWvjseae7OvADDi/QBfg7sxJnz\n",
       "XGAROKWL1au7P84L31jIXocexNYXbgkc0HnRNKxAon2Dsd9ng8N++6/FwkEfAj8Gzgj3/Tu7yqZZ\n",
       "ljOD6Dt/BudPERttpB9aWdRK9k6UtQMq+i1DB+lE3y/Qedg7vZCIYVnItpeADZo8K3dn4IHEvQyr\n",
       "+YRlFiJsH2PYFpmgd4qNeL+JfPHvQhKkj/PS3m9RstFkNq5hIYYvYbgp4hmn231HUF0vfQ+vQksm\n",
       "RNVOmqi0g5LoS5Rt+BjDFMpn43q8xaRf92C1FzfhuQMnheQv7kLuGryWvkcw9IW+XPWvq4Fr7DnD\n",
       "LJyhyOcvunVFiSyJ3BicJciFP8z+aHAi1x0GbpY7WD8q+m1AlKefZO/U+ubG1F47s4DlNO3W1nVI\n",
       "4+ebzsj1GlteORvoQpdlHyMTmH5u/XWsp/1tJKr+IZKg9edOUtaiOy5yUexBuglLUcxF3vslKfoK\n",
       "pSFNVNqBvF6vU2mthFx8nGXMXf0NLnnuce65uPICLNbSyciF9QvAL7np5r8xd/We9vX+NWHR/qDX\n",
       "NmC//WCLKzevOGclWTz9JKIsnkZ7+sOBjaoMqlT0Wxt3FeRDGIywwnqY5GzvJIrci1QuD9go1kci\n",
       "yNcT9rsA+D2SI3qUH20sdz87njoA+b9dX7a3TID6DrARhumUi35S3x0/U+X8zvKU+4cxD2kml4e1\n",
       "A4kC5faw29+l094pE52wSB/gZWZvuAlu16jPyj3I5/Eu4ABmbPUOpdf0ZuTiWJoRb+jJt794IQPe\n",
       "Ws6OZhwmMb+XxdNPIqqCp9H2zgCktUQ1eTkV/RZnTaSDZNC3TpPIrdXeKbLo2yg/Zuq9YQ+k5fIp\n",
       "SGnjdQx98WE2+csCtrlgA2TGaeXxhhWYspYXftFPM+sURPRr8fOxz7U++Yl+kkCtDswQO8/5GLmT\n",
       "8y+8Ei36Ik7hdzXyGv8I+B6Gf+J/TcX3/xVwGoauyAzjG1g8oDtXPjKJeSNXIBMr48hT9KMqeBqd\n",
       "yPWqrKo5d9RsXKhL9Y7bKy+bV0VfCLN2oCGJ3ETRD1sIvFF8kTg/X6ycixBhX4TBxfBb4ED2+PEq\n",
       "TNtxiq1zTyJo72SJ9Gvx80Fe+wZG+p3WjkfQ4onKLbxs/43+rBj+i+Ev9q9gufHdyOv6LaSEtAd/\n",
       "fHwyy3s+yb3nLwBOton1KHLy9IFQe8d1iPb0PwT6JyyzWA1es7vhVRzb6Ej/NuBzeZxIRV+wkX4F\n",
       "SZ5+Honcgkb6bjekXPQBDHtjuA3TWWnicQIyU/a+skcND3DRS//j5pvuSPlk1do7k5DOmbUwj3wj\n",
       "/ayi76/ggehI/xX7b9q7oPKARe4ETkImdQ0BvsayXh3AZKZv3w9Z5q+y+ZthqJ39m7enH4z0+wLL\n",
       "yqqWOnGWU+rcGoM7zreiWhpqifQbLfprkS7hnoiKvjCA8EkWSZ5+nRO5QPPsnS2BdzDOYqTVwHSk\n",
       "S+VJGHogLY5/hMyWrWTB0PdZtmr0hKZyqkjkAjifgHN98n6xeIncRtk7HcRH+nH2DqRPWldOLDQ8\n",
       "iNg4e9uZxGshcycc3t3qdKQ8dohv/y2B/wA3M/I/w6ivvRPl53uksXhuRLrBpqWVIv0R5GQZqegL\n",
       "UULTiDYMSSL3HtLtM5erfAZ2Rqp2TgfuxnAU8FlgG2Si1rXAWZjIfko3QOSEoCCfUl2knwdzkYlE\n",
       "jYr0g1ZiqWxTiLB3nNlIVVRaOyu8DYPhTxjmW6tkBPA2MIuJk+ci79kv7X4HAf9ALuqn86Vjt6Cu\n",
       "9k4q0Y+5mLoO8tplSfh6M8arjfSj+kpZXcir1Nrtg3wvVPRzJKvoNzCR67hItN/oBWi+yPbnvAXs\n",
       "i9eS2DAN6ffyS6TP0u+jD3euA+e5lM81h9KtdpZEbh54F5i8RP9ToKdtExxGB5H2jusQHekDbAXO\n",
       "tIhtQZJ676wpz+MsRXz0IUg30G9juAIp89wJwy3AhQx8ox97/DBsZnQ1hNk7UUlcj6Tmg4OQCqUs\n",
       "oj8AuYPKOdJ3liFVbz2rOG8Y3nrNKvo5EiW8YZ6+X5QWIb2+s/iIaZ/bT4MtHrcXzrJx7HziQcDx\n",
       "GF9EI8naWzEcaCdb5UG1idw88F77nEQ/sRVwB9H2Tj9gefR8gdSCD/Erv2Gf8w37+2xgNQwfIDmS\n",
       "1YCtMLwAgGExD5z+Hptf+0MMtXzWPeph73jCmDXSfznhvFHE2TuQbwVPrqKfpvfOykDUuqMJnr7j\n",
       "guv5+tVOEOpLdHMqj/qJvtRn34JUwlyEtPPdns/9ZiZdls9FbJx6U20iNw/mIWWTcYKTFS8qDVSE\n",
       "ud2RqNJvib2DdObsidgTUVF+VpIi/bWgs4md1xUVDOdRueoYPHvQUnY/8lMkJzCxbJssEvMCJrDa\n",
       "WTRR9k7ceg1JuZJqRH8A0iRy0wzHeCSJvmf9ps1rxeH933JZ7F4jfaFaewdqt3jSRPr1LNs8Col4\n",
       "VgWex3Atn/2/7/D53wxHukJmWQynWqqt08+DuYjNUcsEryBR/vPqss3x3SE5y5AkeQfy5c6jFQQk\n",
       "i74/0vfsnRi69OaNXbyOnmLFGTbBcA9S+XNqhrGF2Tu1RvojEUulAZG+2x0R9Lg5Inkmc0cgroLa\n",
       "OzkSJbx2geOyhExw31rf3ObZO4Y1gBOBwzAcjXjLU/jCr77Ke2PvDenxUi/mAb2sTdZoe+dTohf3\n",
       "qZYo/zmskyuUkrlxfn5WZEWy8OU+oTzSF3snnj7ceOuTSNfcczD8EZnDcS/yuTkEw8Ypx1aNvfM+\n",
       "pYg3jBFIcJRFwAcgpbBD7XKUabHrRMQ2Icxb9F9H7Z1ciVhw2lkO7kLKa5SbIfrTkckpA6RMMTcu\n",
       "RCydVwFp+gVn43ywJ86K3xFofV8/nBXgzkdumVfYplyN4l/AITmfMyrS7yBa9NdB7rZyEn3HBdfz\n",
       "lcM+M8FIvyPhhN534CTgSaS1xvqYzvUSTgEuxvCFFHeHHwHDwD0T+f7MQ1bfuzLmmLcpX4kvyEjg\n",
       "WbDrFaejP/J/n4NchNJaMUnWDuQv+q+gop8r5WuPluP5+vNsxB+cpNIAe8dZAe5UpILniRqeq4Th\n",
       "K8gX5FsV29yuo3C75lXNkpY5SHvjRkb52AvMy4m7ZWMmUh0TpINw0fcnc/Oyd6Bk8YSJfoZI37VX\n",
       "f2cJhvcI/79dChwKfAM6ZwVHMRdptufVng9DLKK4jpfT5XldJ6ItyAjg38AXEp7bzwDktZmJ5FqK\n",
       "LPr/A7bK42Rq7whxwuv39Xsjvc79/m8jIn2otmzT0A0T6Mkv0+0vQtZmXVR+QGLZYL2Yg0RrDRb9\n",
       "uhCVdOwgvN1HPewdiPT13X7IXYVnp5QSueEk992Rbp4/Bs7FJCUcHRecq8A5E5wTwTkKnMNtV9mo\n",
       "Y+YAS4hedcuL9FdLVx/vrmLPuwixjrLYQmlEP+/qnVfJKZGrkb6QVvTD9mtEIheq9/V/DZyA4TVk\n",
       "ss0/gF2Bx2xjriCDgfngpK3EyAtP9BuZxK0XUUnH0chiMkFeR+ydD8lX9MOqz6Azyu+MmJMSuen6\n",
       "7hgety05TkbWJ/b6M+0ILLEzgmthOmLxhAnuCOTOZSESwSetzexF+VCK9NMSNzHLI+9I/1XU3skL\n",
       "tzvSXjVK5Py1+mFJxkZG+j/MdGbDBkirhLWRyHM3pNviOkRPV89rBamsfEpT7J26EJXI7SDc3nkT\n",
       "eY9WIX97Jyw69Pv5kJzIzdJh8xdI+eYcYDtge6QsciMMn8fUtE7224i1FLCByu5OvbusJNHvT6n6\n",
       "ph6Rfk6i7/ZE3sO3UNHPDSu6ke2D/dFSWDlho0Q/W9mmLJjxf0g73elIlPRfpK1CHHmtIJUVz9Nv\n",
       "h0j/fWB1sRAca5+53ZDXNqSxnzPHzvdYm4bYO2V+PkjU2k8CICdswl36ZmuGDzAcibTxuBzYH8On\n",
       "GI4AJmLYwbZ6rgZP9IMMQBbBWQCuJ/pJeRqvBQNUF+k3ytMfjoxvDir6uZEkunW0d9yuiLeaJop6\n",
       "Eykt65NyhacDkC/DHzIOahTNFf02iPSdmeDej1hrJ9oHRwEfxFQmvU5tk/zCiBL9tcFWbAG2UMCr\n",
       "nZ8Zsn+2XvqGG5HmZ34uQT6ThwGXpT5XOTaZW4E/UEno0dPJAMoj/Sz5ssGUXzTDyFP036MzRxCZ\n",
       "yE6NJnJrF/1a3lzrlaZZdNxZjpRtbZC4q2EgcC7wA5tgy8JImmPvtFMiF+BI4FBwPRutg/Akrscb\n",
       "yCSxPCfDJXj6ZcRZPLX30pfV0g4DTscwssqzRJVt+hPgaVfZapVIfwTwvg0WVpBDPx8V/cga/U7q\n",
       "KfoZ2ggD6ZO5ZwG3YphcxZiaae+0SyIXcN4Hfg5cafNGUROzPF4n/9c9racP8cncfFbNkl4+lyKV\n",
       "Y9UQZe/4P7NpRT8Y6WcR/bhVszzsxM6a8V/QcqkISmPvTEDWQO0K/BE4O7B9A+AqJDF4EnC+b9s0\n",
       "5IO3HFgKjKttuHUhrkYfknu9LyB+pmActYm+oQvy2q+JRJFvIV/OPal+Bu8oZNZlo5mDRJrtEumD\n",
       "9C36JlLJ0p140X8JicDzZA4VCUq3C+EJ5bhIP88FVM4AnrUL89ye8dgo0Q9G+mk++8FIvx6J3Dw8\n",
       "+DDRr6mfT5LodwUuRpbNexdZUOFOKMvAf4jcyu4dcrwLjCe5vKmZpLF3hsTs2+hI/2Df3wcjXuSJ\n",
       "SCQ5GvnAH9w5UzI7zYz0oW0ifbCzYg9HEuhTgWtidr4BWRIvT+YAnwk8Nlwer8gLxdXq57c+riyr\n",
       "+X3gDmSN5K6I4yD2j+HfMUfPQGbydrM9izxGUrqIfYBoThL+SH8WMFhybMEeTO5ngNl2PWOPRts7\n",
       "T9nfGxLpj0P6pk+zf98A7EW56M+yP3tEnCOnhQTqRp0SuW5fYBdwbq3huYOUKngMgxEbZ3cM/8tw\n",
       "jiSaLfrtFOkDzlvgGsTSiGlK5ixHPkt5Eubph/n5EG/v5Lk+Lhj+hWGsPe8K+7MtcA2GLTBRr4Oz\n",
       "1FbnjESifo8RlGaqZ/H0X/Wd9xPk/x9MZF8MLAd3D3sRd2hsnb4/0p9LA0R/FOUlZu8AW2c4v4us\n",
       "vrQcydhPjN+9KWT19IMldVFv7lZ43nr8c2cR/deAIeCOBOcU4KZkwXc3QS7Il4IT1xUQ6z2HffAb\n",
       "QZuKPiCls+si9eqNJMzTD/PzQSyDdSPOk1+k72EqxvAihp0R++enMUd6Fk9Q9LMmcv2Ts6Dk6/s+\n",
       "+66D2NZeS+krkNdieYrJi/UQ/XnkMCs3SfRrrSTYHhnwasD9iG/5SMh+xvf7JPvTKLJ6+q8Etke9\n",
       "ucOBUQklVhlF31kK7i2MnXgCsDuJ3qW7DtIF8X/AseD+FrgInKgv8DBgVuDWuVG0ob3j4awgbNHx\n",
       "+hNWshkX6W8bcZ7eNKYtx5FIe+/bMDwcsY83K9dPNYlc/+QsCPf1hyNOxT7Ag7YM1yHZ2oG6iP5f\n",
       "V4HfH4roatUkVe+8S/kLvAbZ2tB6H5RZiF8Zlcg1vp9JGc6fEvdC228kjCz2TtgCH1F1+sORNz2u\n",
       "p3nWSB96fnId257/PVznGExcP293BHAfcBo4X0Gmwm8BvAbuoREHNatGH0pfwHaM9JvFHGC0vfh7\n",
       "rEV4pF//6p0kDB8is86vsv2hwggkcyt6RXkTzZJaxEZF+n7GAM+A8zzwOyTSH0I60c9QveMOB/eY\n",
       "kMe7IgGzvfvYdxo8/DdKWlkVSaL/FJII6kD67O6PJHLDCHr3vShvVLYr8HxVo6wJ1wG+T3Tr2HrV\n",
       "6XtRQ1xNcnbRP37wWOYPczltUcyMQ3cgEuFfCc6l8pgzFZz9kVYMfyg1nCqjWTX60NaRftN4Hkke\n",
       "PwnuDXbOwNqER/pJ1Tv1F30Awx3AY8BvIvYIVvD0B5aVEtPOCuT/krAoTKpIfwtKLR/OscecSPpI\n",
       "P63//gXgVyE9/VdD+vZ7s6QbkshdBhyBCEhX5Eo3FTjcbr8MuTr+B4loVwBHI7bDUEp+djfgOiTy\n",
       "zIh7HnAOONUuZzcE6WkyOGJ7ver0vahhFET2G4kXfcMWSAKwL3IR7QUrevPQKVeyosdBhLaidXsB\n",
       "dyELXJxZud15Btw3gPWRroR+mpXEhfb29JuEsxwRk3OQwOcuJCrOGunnm8hN5mjE5nkbCTJf8fXo\n",
       "fxsJID3COpN6s3LjPstpI/075FdnGbgHI3mZqMDXTxZ7Z4wdz7qU28d2YlYnDUnkgtRsB+u2/dOo\n",
       "3yd8ltw85D9TA64DfA+4m+rXMPXGFtWSNY86/Sh7ZzYE2hqX05eo2m3TOS/ieiQZvsD+fMpb40cA\n",
       "j4D78xD//XIkkjs2JpfwArAxlaLfrGZrUHpdVfRzx5kLnA/uRYgfHDYz2EbHoTmoPOv0kzF8jGFv\n",
       "pFXzkUAXDA8gtk/Q0w8LVNL4+mGR/uaBfcYAp5T+dKaC+xPSifl8ZDW4NG0TxiAN4sZRKfr+C1pD\n",
       "ErnNZiDy5lSzWr1HkuintHdcJ2LfOHvnf1Rv7/zAbvtt5UpEzqs2Wt8VuKf0uPtV5IMzJqG1wwuE\n",
       "rzA0EiITaHXGWW5Xz1J7p244S4CHIrYtAncJEuAEc0WNs3c8DE8B37GNAz8D7ALczDf2+jI33OG3\n",
       "d8Ii/QTRd7tQeYcfiPTdvkgQFLBRPbs0CWc5uEsRlyGm0sd1ENG/Fvnu/tm3MUz044LIVBRd9L3V\n",
       "hGoRfe8DUqXoO0vA9XpeZI307yE50q98bsNwJFGzY8zSc9cCB9Ep+u4QpLnaftJtMJYplE/y8mim\n",
       "vQPwHcpvZ5XG4lk8zRd9D/n8v4Jn8ax/58V0WbIKKzobD0ZF+nGa0RfpeeWfiBX09DcFptRYyeYF\n",
       "hHHlncOR3OptVHY7CBP9mu2dovfeyUP010C8uzjRT+ps6Pn6YSK9EOhpM+0WtyuSQ3iWeNGPspbO\n",
       "Ba5IWJj8JmCCrJsLiPd/AziPxRzj4dk7QZpp7wDOzemazyl1IiqZ22hPP4pLcJjFLr9YSOkOvopI\n",
       "v8LPh0pP31bu1ESaCh7vef4LbBqoOlppRf9T0tXdRrEGIr5RidwkTx/kojCA0DbIjktltD8E8eje\n",
       "It7eqRRZw3hgB+C0+CE5HyFe/77g7gNsCfwy/phOXgdGgBv8QDY70leaS1QrhuZF+n4k6v8um1/T\n",
       "h51P3M0+OpLsoh/080EueAPsBEUor9ypljQVPF5Z6Dzke7mZb5vXVtkjl0Ru0UV/HWSx41rtnWeo\n",
       "3tPHbh8OLIiIRIOiPxyJHGYQGem7DsHyOUMPZPbmTzCpvmTXIitj/QH4Tgpbx+IsQ26ZfZO73F7I\n",
       "Ra3IfZKU+hJVwdPYRG4chvd56NSH2eoPJ2Loh0TDWRO5IZG+sxz57HsXvTwi/TQVPJtTKqiYTPlc\n",
       "prokcosu+msjPTVqtXciRN/tgoh10gd6LvFtf4Nvrjed+z1k4ZOuIccMgxXzMc5qGPbDcBaSZHsT\n",
       "Uncf/Afy/7s+pa3jJ2jx2Cg/137uSmsRYu+4DkWJ9D3+8+NHeW/sW8BldJ8/inwifehcNtHthgRE\n",
       "z9U40jSi77+4pBH9to/0U4i+uwG4+0Vs64YI8POER/q9iY7e/XiiH3VxCBP99+2kio8I+wD2mLsu\n",
       "h2+5CjID+QAkN3AW8PWY5G0AZwliBZ2YtGcIwQoetXaUMHtnFWTyUzNac0QxnRtvmYpLT45cfy2O\n",
       "H/QlDN1926vx9KG0mMr6wDspV6iLI0H03d6IE/GSfeBJOkXfdai0dxrWT79JuN0RIZqMXH2j6l2/\n",
       "CHwduDlk2wjkgzyTcNFP4+dD8gIfQXtnGKUqFM/iKY9G9jz0ULosWwismV7kw3CqXWh6ClIH7dHk\n",
       "JK5SAGZTuTJbsaJ84W0WDRrJKe5XWf3fE/jetl8GjsZwMnAdxhP9SM1IiPQZRO3WDiRH+psAL/lm\n",
       "3E4B1gS3PxKQLwlYtm0f6a+J2A2fIF06o7ys1YleQnANpEGTt/5nkLRtEKq1d0CEtDyZaxjKevfs\n",
       "x/3n3F6b4NdEhL2jrMSERfrF8fNLeMsmjuCdbd7BsCvSBfMXwNG2oeAKogUyKdLPw8+H5OqdwPM4\n",
       "y5AZv1tSORsXVoJE7jqUpovHrWwzClgN3LDqHE/0FwJdwF01sD1NuSYki/4CQu0dQEQ/mMw9m1f2\n",
       "mMZruz2Z4rnrxVvAQF/Jp4q+EpbILWKk/w7y3S7dQRv+hawYd5Lt1R+3QHpSpJ9H5Q4kV++EXVw8\n",
       "Xz+sFHU+pYmiVVNk0ff3/Y4T/dWRq3pYtG/7bjsu4dF+lkh/VMy+dsp1J2H2jmDYHtiFuy/9lPAe\n",
       "KA3CWYHcTnrRvto7Slid/gjSBUYNxFmAfBc3xy+M0qP/SOAGVvl4NtGa0ahIP8neySj6zhLE9ahp\n",
       "cfRWEf24GXarIy/UhiHbvEgfpDNeUPTz8vTT2TuGbkhJ5nEsGjiapoo+IKLvJXM10lfC7J2fIqXB\n",
       "ReNtZEGn8s+s4QbgEfY5cATVRfpbIsKax8zwGNF3uyKzfoMVQnGRPuTg6xdd9F+3v0dE+q6DRKj/\n",
       "JDzS94t+rZF+75h9o+r0oTzS/zEwm9++fTvy5cqyNkE98Pv6GukrnwKrgGsjSXccEk3/sYljisIT\n",
       "/TBhPIrhz/Rh9x99JeJY/6LofmYifX6ezql0OS7SXwdZsCg4jmlIJL8lK6no++2dsKv2QGAJ0vc/\n",
       "xt4BokU/rafv/zeI7811e9jzSs/tPjNmsuEt62P4M9JP5wjmrNEh4wouwtxwbNmm6xA+s1FZqXBc\n",
       "5HPr+fqnAGeCs7h5Y4pkOjCasLtTw3zuvvQWxlzzDUzo6nL+RdH9eIFaHtYOxIt+hIXkuEi0vxvh\n",
       "38eak7kFFX3XIV0id3UkWp5Ksr3zEZWtGLJE+v5/g/jf3KHALIzTF8NFHLv6vXzu7FFIDe4GGKYS\n",
       "vU5po/HsnQFIeVjRqjSUxmMtHndbZILSlU0eTxReMBceqLzylRd47OePApMw7BvYGhXpf4hYO3mJ\n",
       "flz1zuYxzzMZsZ6jIv2aZuUWVPQZhKzP+7H9O0r0PUviTfndvxqUuyry5np9+MMi/bSefpLo++0d\n",
       "z88/ExjB8u5bM3Hycoz7R0ynz18U0Z8BdEeiDrV2FCitOnUKcEZBo3woiX5UHuoD/nXyx8ha0udg\n",
       "+D2mMwEaEek7K+x5/5vTGJMi/eB6Fh6T7b8rlb1jRbHTV0uI9J2lsj/rhWzzZtvW4uknLeXnf3OH\n",
       "0fHQAmBf4PucvuQNKnvw+PMVTcRxEYtnVzSJqwizgL0Rb/vq5g4lFu8OPsqSlAla0pd/LHLX/wiG\n",
       "DqIjfYDNZL2KXIgr2YyrEPqP/XelE32/KCbZOyAWj9/X91s7EF69k9XTj7I/SnX6XZaM4Ms/WA+p\n",
       "F/aalwUnaBUl0gcVfaWcWcji5KfbEsGiMg35XkZ9f0utGAyfAPsAN+M6/6TnJ12BReGH5WpxRkT6\n",
       "7mr28bAVzABnNrAzOB+HbGxb0ff7+ZBs74D0r/D7+kHRr7V6x/9vkFKd/q4/342uS5ZR7oWGRfpF\n",
       "Ef0pSCSk9o4CYu9Mo5hlmj6c94BNY6psyvvvGFwM57K476Ps/d0VmJrmN6Ulyt6xnTXjKoScByM2\n",
       "tGsit0IU5wDdbftfP/5I/yXKI31/5Q6EJ3Lz8vTlzTUMYsxVE3jo1Osw+Ju4+WblhrRUbi4v2H81\n",
       "0ldAelgd7OsHU2CciEgZkIvXwIoOt1c+ch5DXnKAw+s5MkuU6G9G9R082zaRGxB9xyV8glYWe6eW\n",
       "SH8BMus3KZF7Gm/u/B7PffupwPYZlOyd1YBF4ISVjDWDKfZfjfQVwJkCzqPNHkXtOMsQ3768rcQH\n",
       "m63KX69/BTgNU7EQet5EVe+sA1SbN2hbeycs0Rlm8YyiJPovA+vZHvmQXvRTePqOi3yAooR6Puvd\n",
       "ORzYl7v/7wNKs3E9/P13imTtAM4s5IKqNfpKuxHWYnkAM8fMBH4C3ISpfVGSGKIi/Q7EQquGdhR9\n",
       "t7crI/YAABCsSURBVAcyBXl6YENA9N0+SK9vm+xw5iLC7i2EHmbvVFuyCbANOB9UPGoYw/GDzmHP\n",
       "768NfJf5w4dQOYW7wKIPSGvq4N2JorQ6YaIvLRgM1wGPAJdj6qaDtsCjokHaWlRv7zZE9Ccgfvmr\n",
       "wPEh2zdAFjpZBByb8dgwRgPvhniKwVm5NsovS4b4k7nBSH8+0K28lj+1vUNFGZdhfQy3AH9n4aCH\n",
       "ueCNGRj+RnmzNQ+/vVNA0Xf+1RoerqJkIjzSL5VrHoVYxFdgCFvdrkac5cBiZBlSi+sgkX61ol/3\n",
       "RG5X4GJEvDcCvknlzNcPkc5251VxbBhRohi0d/zWjodN5rr9AYcyO6az0+ZA+dt1yCT6Pgy9gfuQ\n",
       "6Hgdfv/axSzr1dsmmntSaQPZks3OJG7BRF9R2pI3oKINQ6nZmmEBok9rAH8OrL6VF0GLZ5g8VnVp\n",
       "aN0TueOA1xD/aSlwA7BXYJ9ZiPgFI8U0x4aRVvRXpzL56LVj8LdU9uOv4FlFxlVVhHsS8CiGs+wH\n",
       "x6vTHybjDD6vs9DuMxgVfUVpFPcjc1D8lLdVNswHvoJYvTf6Zu3mRVD0O6itcq/u9s4oyi2Sd6hc\n",
       "ECTvY6Nmq4aJfkSkX2HtePh9/Sx+fgnD+sD3gZ/5Hl2IRPhhq914eBaPir6iNIbHgA0DCyxVtlU2\n",
       "LEQmbwHchqFHjmMIVvCsRfVJXO98dRX9WtqLZjnWlH5uHkfN9k4q0c9u7Rgc4PfAmRh/XbuzAhH+\n",
       "tYkW/Xft9mEh41YUJXecxcDDyDraHuELqBgWA/sjDdfOznEQwUi/liTueFj3YPhZB6KXVZEk+u8i\n",
       "AuqxBukFK8uxpvSzX3+qt3feQyLuLSiv3PEIin7WFYH2QS42F4VsW4CIerBc02MGsB0w3dYQK4pS\n",
       "f+4FvuT7O2oBFTAsBQ4Gvoph75yeP9h/pxbRnwSv/QbOXUQdRf8ppPFSB9ADuRLeGbFvsCwpy7GW\n",
       "2ERnCnvHcZFofxfCI31//51skb4kb38L/Nh+OILMJznS/xxq7ShKI7Gi31k2GbVUoiD9sr4BXGab\n",
       "s9VKnpE+NCCRuww4AnnhXgRuRJKlh1OaxjwcEdifAr9EIuw+McfGMRhYHtFo6GM5r7eqT6i9g32O\n",
       "tYm2dzx/L6un7yVv/xWxfQEy0y5O9D+Lir6iNJJXkYWWvBXioiN9D8O/EYvnxhz8/bBE7rQaz9en\n",
       "lsXRu6XY5+/2x89lvt/fp9zGSTo2jjWJ7jy3AtxZwFBwZyIRe+VkKYn0IdreGW1/TxfpG/ogH4A9\n",
       "keXZovAi/Th7pzsq+orSQBwX3HuRKp4XSIr0S/wOGA+cReX8oyz4V9XrimhlXM+gBJwl4Na0OHrR\n",
       "ZuSOJv4F8SyekcD7EcsNeqIfdheQzdM37IQ0RuoFbFaevK1gfmlcoXj5BxV9RWksfl8/OdIH6coJ\n",
       "hwBfC1l5Kwv+6p1RwGxwIto6Zzpn1RU8aSL9RjKa8Ajdw5uV25PopPALyAVhYci2dNU7ckv3OyS6\n",
       "PxzDPYkjF3sH4ks2QUVfURrNg8C1dvJkP9IWcBg+wvBV4D4MMzA8XsVz++2dWv18j5pm5bZqpB9W\n",
       "uWNxXkNWpQkjbZ3+Mcgkr01TCj7Im+uNMYyZyIVBRV9RGorzKbJK1R7AgkzVc4angW8Dt2LKVuZL\n",
       "i796p4N8RL+mSL9ooh/j6QPloh9TOupECe+HlBK54ZG+YU3gOOBQu+JOWuYTO73aWQ6sDU6WcyqK\n",
       "kg/3IY0Fs3//DP9ACjnuwVT08kkiGOlPy/z8ldRUwVM00U9j7wwjunIniTSe/u+A32MyX5EXEG3t\n",
       "WCIvRoqi1Jd7kUi/unUsDFcA1wN32fLttNTD3mmrSD8HeyeWeUBPW/ZZGekbJiDW0DlVnHs+0daO\n",
       "oijN5b/IrPla7rRPRkrCr7Kz89Ogoh+N2wsR4rAyTI+U9k4UZZ02yz19abR0EXAkJmrR5FhSRPqK\n",
       "ojQHZznSgK36FeukoucHyKTTH6Q8yl+904EmcstYE2lRsCJmn1rtHShZPMFI/2fACxkSt0HmoOvM\n",
       "KkqRuZXwZo7pkYBwf+BUTGTBiB8b6bs9yK/vVrtE+onWDshdwAhkFnC1Auu1Yih5+oZ1kOXTflrl\n",
       "OQEmAqfUcLyiKHXFuQmco2o+jeEV4Ghkxm6S+Hr2zprAjJz6brVNIjepcgc6V7jnE5mZVhVeKwaJ\n",
       "9A3jgH8Bv8LUkll35oAzu/rjFUVpGQzXI8stXpLg73slm3n5+bByRfrOMkT4a7lF8uydfhy19peA\n",
       "u5EmapfUcE5FUVY+jgLGIjN3o/AifRX9EJLKNT1mUqvod108mF2PG8jAaccBX8BwRw3nUxRlZURW\n",
       "zfs6cA6GzSP28kS/g/xEv60SuWkaEc2kunJNj484cMJBjPhvF5b22grDCzWcS1GUlRnDFCTivwXj\n",
       "rb9dhle9k9fELO+cbSH6aRK5UGukv925Qxg6ZR2uv+sTzpz/YdXnURRFATD8BbGJ/4yp0NQFwKpI\n",
       "B9487Z1WT+S63ZAOlWnEfCJwe1VPY+jJDmfsx73nv8PSPtnXx1UURQnnZ4gQ/6r8YWcFsBjp5aWe\n",
       "vo8RwKx0FTnOJHBerPJ5jmVJ77d57oABZF8qUVEUJRxZTe/rwGEYdg9snY+sHpjX5M22EP201k71\n",
       "GNYCjuHRX5wJXYaSdVF0RVGUOAzvI8J/VWCpxfnAWwkTT7PQFonctJU71SF1tBcBv2XykV7iVkVf\n",
       "UZR8kZ775wLXYzrXK5lPftYOtEmkn7Zyp1r2RNavPQ+p0wcVfUVR6sNvEfvY2L/nkb/ot3oit472\n",
       "juEzwP8hE7CWIG/GctTTVxSlHhhWAAcD37VLruYd6fsXZslMkUQ/f3vHsBHwEPBrDA/Kg52dNjXS\n",
       "VxSlPhhmIjN1r6Xf9HnAK/md3FmCBK5VUSTRzzfSlxly/wR+YRdA8KOiryhKfTHcB/yFo9fpSsek\n",
       "u3I+e9X6VRTRz9fTN3wWWR7taAx/DtlDRV9RlEbwS7ouHcohO/0o5/NGLMuaTBrRnwC8BLwKHB+x\n",
       "z+/t9meBLXyPTwOeA54GJsc8xzLpUpkD0ib5b8D3MdwcsddHqKevKEq9kTziNxGLOao/TzVULfpJ\n",
       "dAVeQ5oFdUdWlN8wsM/u0LnwyNbAv33b3qS0Jm0ULrjP1jxSkNJMw30YfpbwlCeDu3Muz6koipKE\n",
       "4QAML2VcXzcG99+AW82RSZH+OET0pwFLgRuAvQL77AlcY39/EhiArBDjkWYtybysnW8CQ4EL4ndz\n",
       "TgHngZyeU1EUJR7DdYg+XpjTGetm74wCpvv+fsc+lnYfF0mmPgUcFvM8tVfuGAYB5yO2ztKaz6co\n",
       "ipIvRwA7YNg/h3NVnZPslrA97e1DVDT/OWRZw9WQRYlfQlabCbD3epQmMkyyP1k5G/grJjZ3oCiK\n",
       "0hwMczF8E/g7hsmYzLX74+0P8L21cx2bj22Af/j+PoHKZO6lwDd8f79Eub3jcTJwbMjjLrhfr2WQ\n",
       "GD6P4R0M/Wo6j6IoSr0x/BTDE742DVXgXkKdPP2ngM8gidweyCrwdwb2uRM4yP6+DfAJ0vO+F6Wp\n",
       "wr2BXYHnI56nek/f0BO4HDgKoxU5iqIUnguR6sFf1nCOunn6yxAf6l7gReBGYCpwuP0Bqdx5A0n4\n",
       "XgZ49ajDESvnGSSBcTdSOx9GdZ6+NFK7EHgZuK2qcyiKojQSadNwCPADDNtVeZaq2zqkqaypNy64\n",
       "XatqO2r4NbA3MF6jfEVRWgrDXsDvgDFV6pdLFRpekBm5VQn+YUhTo91V8BVFaTkMdyDux8WNfNqC\n",
       "iH5G5Ap5KjDBLlygKIrSihwLjLNVPQ2hIPZOhnGIB3Y7EuE/Va9BKYqiNATDWKRK8rOYTPnNquyd\n",
       "1hJ9Q1+kAugIDHfXc1CKoigNw3Ai8EXgizbRm4ZW9vRTcxbwkAq+oihtxtlIWfxP6v1ErSP6hs8D\n",
       "XwWOafZQFEVRcsWwHPg2cAKGTev5VK0h+oZVgSsQW+fjZg9HURQld6Qtw8+B6+yk07rQGqIvfXme\n",
       "wegELEVR2pqrkYmup9XrCYov+rIK1iHAkU0eiaIoSn0xuMD3gQMw7FiPpyi26Bv6IFe+Y+1Cw4qi\n",
       "KO2NYTbSiv4aDP3zPn1xRV/66lyFLLN4XZNHoyiK0jgM9wB/R5aizZXiir60cV4T+JG95VEURVmZ\n",
       "OA7YDsO+eZ60mJOzDHsg7ZLHYXi3GYNSFEVpOoatkfb1YzC8F9jaJjNyDesjLZn3xvB4swalKIpS\n",
       "CAynIOuV7x5wPdpgRq4kbm8HTlLBVxRFAeB0YAjwgzxOVizRh/OAyRgmNnsgiqIohcCwFJmtexqG\n",
       "dWs9XXFE37AbsBtwVLOHoiiKUigMLyETtq7B0LWWUxVD9A2DgT8C38HwabOHoyiKUkAuAhYjVT1V\n",
       "UwzRhz8AN2N4sNkDURRFKSTScvk7wHEYNqv2NEUR/c2RunxFURQlCsNbSFO2P1V7iqKI/kEYFjZ7\n",
       "EIqiKC3A1cC0Jo+hJnS2raIoShYMg6ijdk4AXgJeBY6P2Of3dvuzwBYZj1XRVxRFyU5dtLMr0tu5\n",
       "A+gOPANsGNhnd+Ae+/vWwL8zHAsq+n7GN3sABWJ8swdQIMY3ewAFYnyzB1AgqtLOJE9/HCLc04Cl\n",
       "wA3AXoF99gSusb8/CQwAhqc8VilnfLMHUCDGN3sABWJ8swdQIMY3ewCtTpLojwKm+/5+xz6WZp+R\n",
       "KY5VFEVRGkiS6Ke9fShC4zZFURQlgW4J298F1vD9vQYSscfts7rdp3uKYwFeR319Pyc3ewAFQl+L\n",
       "EvpalNDXQni9HiftZk/cAfQgOZG7DaVEbppjFUVRlIKxG/AykpT1Zs0ebn88LrbbnwXGJhyrKIqi\n",
       "KIqiKEo7Ucskr3Yj6bU4AHkNngMeg+qbK7UAaT4XAFsBy4B9GjGoJpHmtRgPPA28AExqyKiaQ9Jr\n",
       "MQT4B2IbvwAc0rCRNZYrgZnA8zH7FFI3a5nk1W6keS22Bfrb3yewcr8W3n4PAncDX2vU4BpMmtdi\n",
       "ADAFKZYAEb52JM1rYYCz7O9DgA9JLkxpRT6PCHmU6GfWzUY1XKt2ktewBo2vkaR5LZ6AznUFnqT0\n",
       "JW830k7gOxL4KzCrYSNrPGlei28Bt1CqgpvdqME1mDSvxXtAP/t7P0T0lzVofI3kEeDjmO2ZdbNR\n",
       "ol/tJK92FLs0r4WfQyldyduNtJ+LvYBL7N/tWt6b5rX4DDAIeAh4CllCrx1J81pMBDYGZiC2xtGN\n",
       "GVrhyKybjbodqnaSVzt+wbP8n3YCvgtsX6exNJs0r8UFwC/svg7tOxEwzWvRHamO2xnohdwR/hvx\n",
       "c9uJNK/FiYjtMx5YB7gfWZdjbv2GVVgy6WajRL/aSV7v1nlczSDNawGSvJ2IePpxt3etTJrXYkvk\n",
       "9h7Eu90NueW/s+6jayxpXovpiKWz0P48jAhdu4l+mtdiO+AM+/vrwJvA+sgd0MpEYXWzlkle7Uaa\n",
       "12JNxNPcpqEjazxZJ/BdRftW76R5LTYA/okkOnshyb2NGjfEhpHmtfgtpZm5w5CLwqAGja/RdJAu\n",
       "kVs43axlkle7kfRa/BFJTD1tfyY3eoANJM3nwqOdRR/SvRbHIRU8zwNHNXR0jSXptRgC3IVoxfNI\n",
       "krsd+QuSt1iC3Ol9l5VXNxVFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFURRFUZrH\n",
       "/wPWmDNEEKgBiQAAAABJRU5ErkJggg==\n"
      ],
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10f500e90>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import scipy.sparse\n",
    "import scipy.sparse.linalg as spla\n",
    "import scipy\n",
    "from scipy.sparse import csc_matrix\n",
    "\n",
    "n = 80\n",
    "ex = np.ones(n);\n",
    "lp1 = n**2 * sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); \n",
    "rhs = np.ones(n)\n",
    "ev1, vec = spla.eigs(lp1, k=4, which='LR')\n",
    "ev2, vec = spla.eigs(lp1, k=4, which='SR')\n",
    "lam_max = ev1[0]\n",
    "lam_min = ev2[0]\n",
    "\n",
    "tau_opt = np.real(2.0/(lam_max + lam_min))\n",
    "\n",
    "sol = spla.spsolve(lp1, rhs)\n",
    "niters = 1000\n",
    "x = 0.2*np.random.random(n) #random x0 to add high frequencies to the error\n",
    "err_all = []\n",
    "grid = np.linspace(0, 1, n)\n",
    "for i in xrange(niters):\n",
    "    rr = lp1.dot(x) - rhs\n",
    "    x = x - tau_opt * rr\n",
    "    err = spla.spsolve(lp1, rr)\n",
    "    err_all.append(err)\n",
    "init, = plt.plot(grid, err_all[0])\n",
    "smooth, = plt.plot(grid, err_all[200])\n",
    "plt.legend([init, smooth], ['No smoothing', '200 iterations'])\n",
    "\n",
    "#x = np.linspace(0, 1, n)\n",
    "#def set_cursor(i):\n",
    "#    plt.ylim(0, 1.1*np.max(err_all[0]))\n",
    "#    plt.plot(x, err_all[i])\n",
    "#    plt.xlabel('space')\n",
    "#    plt.ylabel('error[i]')\n",
    "#    plt.title('Error smoothing')\n",
    "#from IPython.html.widgets import interact\n",
    "#interact(set_cursor, i=(0, niters-1, 5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Geometric Multigrid\n",
    "\n",
    "Now we know that in an iterative processes \n",
    "* High-frequecies (k>N/2) in error decay twice at each iteration \n",
    "* Low-frequencies in error decay really slow\n",
    "\n",
    "But, we can approximate **smooth error vector on coarser grids**!\n",
    "\n",
    "Equation on error: $A_h e_{h} = r_h $, where $r_h = f_h - A_h u$ is residual. Given $e_h$ we can simply get $u_h = e_h + u$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### A Two-Grid Cycle (almost multigrid :)) \n",
    "\n",
    "1. **Smoothing.** Do several iterations (say 2-3 Jacobi iterations) to reach $u$ with smooth $e_h = u - u_h$\n",
    "2. **Restriction.** Restrict the residual $r_h$ to the coarser grid $r_{2h} = Rr_h$\n",
    "3. Solve $A_{2h} e_{2h} = r_{2h}$\n",
    "4. **Prolongarion** (Interpolation). $e_h \\approx P e_{2h}$\n",
    "5. $u:= e_{h} + u$, then repeat from 1 if necesarry\n",
    "\n",
    "Here \n",
    "* $R$ - restriction matrix\n",
    "* $P$ - prolongation matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "For the 1D case, N = 7:\n",
    "$$\n",
    "   P = \n",
    " \\frac{1}{2}\n",
    "\\begin{pmatrix}\n",
    "1& \\\\\n",
    "2& \\\\\n",
    "1& 1& \\\\\n",
    " & 2& \\\\ \n",
    " & 1& 1 \\\\\n",
    " &  & 2 \\\\ \n",
    " &  & 1\n",
    "\\end{pmatrix},\n",
    "\\quad\n",
    "R = \\frac{1}{2} P^{T}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "In the 2D case:\n",
    "$$\n",
    " P_{2D} = P\\otimes P, \\quad R_{2D} = R\\otimes R\\equiv \\frac{1}{4}P_{2D}^{T}, \n",
    "$$\n",
    "which corresponds to the following stencils\n",
    "\n",
    "$$\n",
    "\\text{Restriction:}\\quad\n",
    "\\begin{bmatrix}\n",
    "  \\frac{1}{16} & \\frac{1}{8} & \\frac{1}{16} \\\\\n",
    "  \\frac{1}{8} & \\frac{1}{4} & \\frac{1}{8} \\\\\n",
    "  \\frac{1}{16} & \\frac{1}{8} & \\frac{1}{16} \n",
    "\\end{bmatrix}.\n",
    "\\quad\n",
    "\\text{Prolongation:} \\quad\n",
    "\\begin{bmatrix}\n",
    "  \\frac{1}{4} & \\frac{1}{2} & \\frac{1}{4} \\\\\n",
    "  \\frac{1}{2} & 1 & \\frac{1}{2} \\\\\n",
    "  \\frac{1}{4} & \\frac{1}{2} & \\frac{1}{4} \n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "$A_{2h}$ can be chosen either as an operator on a coarser grid (from differetial formulation) or as Galerkin projection $R A_h P$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#####  Finally... Multigrid\n",
    "\n",
    "Multigrid need not stop at two grids! \n",
    "\n",
    "To find error on $2h$ we can go to $4h$ and so on. However, there are several different ways to go through grid levels.\n",
    "\n",
    "<img width=\"700\" src=\"fig/mg-cycle.jpg\" />\n",
    "[<a href=\"http://cvsp.cs.ntua.gr/research/gac/\">source</a>]\n",
    "\n",
    "1. **V-cycle**. Straightforward generalization of the two-grid cycle\n",
    "2. **W-cycle**. On each coarse grid level V-cycle is called twice\n",
    "3. **F-cycle** (full multigrid, FMG). Starts from the coarsest grid. Then prolongates to the next level to get a good initial approximation to run V-cycle\n",
    "\n",
    "<!---\n",
    "4. **S-cycle** (sawtooth cycle). No smoothing during restriction. Then starts from the coarsest grid and does prolongation procedures as it was in the V-cycle\n",
    "-->\n",
    "\n",
    "**W-cycle** is typically superior to a **V-cycle**. At the same time the **F-cycle** is asymptotically better than $V$ or $W$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Pseudocode (V- and W-cycles)\n",
    "\n",
    "```\n",
    "def MGM(l, u_old, rhs):\n",
    "\n",
    "    if l==0:    \n",
    "0.      u_new = solve(A0, rhs)     \n",
    "        \n",
    "    else (l>0):    \n",
    "1.      u = u_old\n",
    "        u=several iterations of smoother\n",
    "        \n",
    "2.      res = R (A_k * u - rhs)  (residual on a coarser grid)\n",
    "\n",
    "3.      e = 0 (initial guess for error)\n",
    "        for i = 1..gamma   #gamma=1 is V-cycle, gamma=2 is W-cycle\n",
    "            e = MGM(l-1, e, res)\n",
    "            \n",
    "4.      u_new = u + P e\n",
    "\n",
    "    return u_new\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Complexity of one cycle\n",
    "\n",
    "Let $d$ be dimension of a problem. Then comutational cost of one cycle is\n",
    "* V-cycle: $\\mathcal{O} \\left(N (1+2^{-d}+2^{-2d} + 2^{-3d} + \\dots  )\\right) = \\mathcal{O} \\left(\\frac{N}{1 - 2^{-d}}\\right)$\n",
    "* F-cycle: $\\mathcal{O} \\left(\\frac{N}{1 - 2^{-d}} (1+2^{-d}+2^{-2d} + 2^{-3d} + \\dots  )\\right) = \\mathcal{O} \\left(\\frac{N}{(1 - 2^{-d})^2}\\right)$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algebraic Multigrid\n",
    "\n",
    "What we have considered before is called **geometric multigrid** (GMG). \n",
    "\n",
    "There is an alternative **algebraic multigrid** (AMG), where only sparse matrix is given (with no information concerning a grid structure). Then the grids and operators have to be constructed **automatically**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### Why AMG?\n",
    "\n",
    "1. Difficult to build hierarchy of grid levels on **unstructured grids** (e.g. on irregular triangulation)\n",
    "2. Uniform coarsering with linear interpolation do not take into account operator structure. As a result, smoothers sometimes can not produce good approximation to be restricted on coarse grids. That is why geometric version may fail on\n",
    " * Poisson equation with **discontinuous coefficients**\n",
    " * **Skewed operators**\n",
    "3. Some problems are purely discrete with **no geometric background**\n",
    "4. AMG is **black-box**, so one does not need to know how to programm multigrid (AMG is more expensive than GMG though)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### AMG Idea\n",
    "#### Smooth vectors\n",
    "Vector $u$ is smooth if $\\|A u\\| \\approx \\|u\\|$. This is due to the fact that $A$ amplifies high-frequencies.\n",
    "\n",
    "####How to find neighboring points? \n",
    "Nodes $i$ and $j$ are said to be **strongly connected** if $A_{ij}$ is big enough (e.g. greater than $\\alpha A_{ii}$ say with $\\alpha = 0.1$ ).\n",
    "\n",
    "#### Coarse grid\n",
    "1. Remove all weak connections\n",
    "2. The fact that $i$ and $j$ are strongly connected means that $u_j$ strongly influences $u_i$. So, we do not need to have both $i$ and $j$ nodes in the coarse set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "##### PyAMG\n",
    "\n",
    "To install just run\n",
    "```\n",
    "conda install pyamg\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Detected a Hermitian matrix\n",
      "    maxiter = 400\n",
      "    iteration 1\n",
      "    iteration 2\n",
      "    iteration 3\n",
      "    iteration 4\n",
      "    iteration 5\n",
      "    iteration 6\n",
      "    iteration 7\n",
      "    iteration 8\n",
      "    iteration 9\n",
      "  Residual reduction ||r_k||/||r_0|| = 5.26e-08\n"
     ]
    }
   ],
   "source": [
    "import numpy as np  \n",
    "import pyamg\n",
    "\n",
    "n = 200\n",
    "A = pyamg.gallery.poisson((n, n), format='csr')   \n",
    "b = np.array(np.arange(A.shape[0]), dtype=float)\n",
    "x = pyamg.solve(A, b, verb=True, tol=1e-8)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary\n",
    "\n",
    "1. Multigrid is a strong tool to solve elliptic problems (but not limited to them!) with **optimal $\\mathcal{O}(N)$ complexity**\n",
    "2. Key idea: \n",
    "  * **smooth error** with iterative method\n",
    "  * smoothed error can be well approximated on **coarser grids**\n",
    "3. Sometimes it is not possible to use geometric interpretaion of grid hierarchy. In these cases one can use the **algebraic version**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 100%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 100%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 100%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/alex.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}