lecture-4.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:ab20b32017f664ad9dd21fdef8c975a9c4f8fab35650f3842ba8100ee3a66c55"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Lecture 4: Linear systems"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Syllabus\n",
      "**Week 1:** Python intro  \n",
      "**Week 2:** Matrices, vectors, norms, ranks  \n",
      "**Week 3:** Linear systems, eigenvectors, eigenvalues"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##Recap of the previous lecture\n",
      "- Scalar product\n",
      "- Unitary/orthogonal matrices and norm conservation\n",
      "- Matrix rank definition\n",
      "- Skeleton approximation and dyadic representation of a rank-$r$ matrix\n",
      "- Singular value decomposition"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Today lecture\n",
      "Today we will talk about:\n",
      "- Linear systems, inverse matrices\n",
      "- Gaussian elimination\n",
      "- Sparse matrices\n",
      "- Condition numbers\n",
      "\n",
      "Before that I will spend a few minutes on certain techinal notes on the HW1. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Use numpy arrays in homework!\n",
      "-  Use \n",
      " ```\n",
      "     a = np.zeros((n, n)) #Note for the double brackets\n",
      " ```\n",
      " instead of\n",
      " ```\n",
      "     a = [[0 for i in xrange(n)] for j in xrange(n)]\n",
      " ```\n",
      " In this case you can not add arrays like\n",
      " ```\n",
      "     a = a + a\n",
      " ```\n",
      " and have to write down functions/cycles instead one line\n",
      " \n",
      "- Use numpy arrays indexing like ```a[i, j]``` instead of ```a[i][j]```\n",
      "- Use ```%matplotib inline``` and ```plt.plot()``` (without ```plt.show()```) in the Notebook\n",
      "- For code writing, there is a [PEP8 style guide for Python code](http://legacy.python.org/dev/peps/pep-0008/) - try to follow\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Motivational video on algebra\n",
      "And before we start - a very simple and intuitive video on algebra. It also contains a short ove"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from IPython.display import YouTubeVideo\n",
      "YouTubeVideo('8CX-Q0gtSp8')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "\n",
        "        <iframe\n",
        "            width=\"400\"\n",
        "            height=300\"\n",
        "            src=\"https://www.youtube.com/embed/8CX-Q0gtSp8\"\n",
        "            frameborder=\"0\"\n",
        "            allowfullscreen\n",
        "        ></iframe>\n",
        "        "
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 2,
       "text": [
        "<IPython.lib.display.YouTubeVideo at 0x104457a10>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Linear equations and matrices\n",
      "A linear system of equations can be written in the form\n",
      "\\begin{equation}\n",
      "\\begin{split}\n",
      "    &2 y + 3 x = 5 \\longrightarrow \\quad &3x + 2 y + 0 z = 5\\\\\n",
      "    &2 x + 3z = 5 \\longrightarrow\\quad &2 x + 0 y + 3 z = 5\\\\\n",
      "    &x + y = 2 \\longrightarrow\\quad  & 1 x + 1 y + 0 z = 2\\\\\n",
      "\\end{split}\n",
      "\\end{equation}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "##Matrix form\n",
      "\n",
      "$$\n",
      "\\begin{pmatrix}\n",
      "3 & 2 & 0 \\\\\n",
      "2 & 0 & 3 \\\\\n",
      "1 & 1 & 0 \\\\\n",
      "\\end{pmatrix}\\begin{pmatrix}\n",
      "x \\\\\n",
      "y \\\\\n",
      "z \n",
      "\\end{pmatrix} = \n",
      "\\begin{pmatrix}\n",
      "5 \\\\\n",
      "5 \\\\\n",
      "2\n",
      "\\end{pmatrix}\n",
      "$$\n",
      "The \"matrix form\" is\n",
      "$$\n",
      "A x = f, \n",
      "$$\n",
      "\n",
      "where $A$ is a $3 \\times 3$ matrix."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Linear systems in applications\n",
      "Linear systems are everywhere - in modelling, even if you have nonlinear systems, by **linearization** they are reduced to a sequence of linear systems. They appear in different applications:\n",
      "- Circuit modelling (Kirchoffs law)\n",
      "- Photonics (Maxwell equation, electrodynamics)\n",
      "- Computational fluid dynamics (Navier-Stokes equation)\n",
      "- +$\\infty$ more"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Linear systems are big\n",
      "We take a continious problem, discretize on a mesh with $N$ elements and get a linear system.  \n",
      "Example of a mesh around A319 aircraft\n",
      "(taken from [GMSH website](http://geuz.org/gmsh/)).  \n",
      "<img src=\"a319_4.png\" width=50%>\n",
      "The main difficulty is that these systems are big: millions or billions of unknowns!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Linear systems are structured\n",
      "Storing $N^2$ elements of a matrix is prohibitive even for $N = 100000$. How to work with such matrices?  \n",
      "Fortunately, those matrices are **structured** and require $\\mathcal{O}(N)$ parameters to be stored.  \n",
      "The most widespread structure are **sparse matrices**: in the matrix there are $\\mathcal{O}(N)$ non-zeros!  \n",
      "Example (one of the famous matrices around for $n = 5$):\n",
      "$$\n",
      "  \\begin{pmatrix}\n",
      "  2 & -1 & 0 & 0 & 0 \\\\\n",
      "  -1 & 2 & -1 & 0 & 0 \\\\\n",
      "  0 & -1 & 2 & -1 & 0 \\\\\n",
      "  0 & 0 &-1& 2 & -1  \\\\\n",
      "  0 & 0 & 0 & -1 & 2 \\\\\n",
      "  \\end{pmatrix}\n",
      "$$\n",
      "At least you can store such matrices, and multiply by vector fast; but how to solve linear systems (and that is typically the final goal)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## How to solve linear systems\n",
      "**Important**: forget about determinants and the **Cramer rule** (it is  are good for $2 \\times 2$ matrices still)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "## How to solve linear systems\n",
      "\n",
      "The main tool is variable elimination. \n",
      "\\begin{equation}\n",
      "\\begin{split}\n",
      "    &2 y + 3 x = 5 \\longrightarrow \\quad &y = 5/2 -  3/2 x \\\\\n",
      "    &2 x + 3z = 5 \\longrightarrow\\quad &z = 5/3 - 2/3 x\\\\\n",
      "    &x + y = 2 \\longrightarrow\\quad  & 5/2 + 5/3 - (3/2 + 2/3) x = 2,\\\\\n",
      "\\end{split}\n",
      "\\end{equation}\n",
      "and that is how you find $x$ (and all previous ones).  \n",
      "This process is called **Gaussian elimination** and is one of the most widely used algorithms. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Gaussian elimination\n",
      "Gaussian elimination consists of two steps:\n",
      "1. Forward step\n",
      "2. Backward step"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "## Forward step\n",
      "In the forward step, we eliminate $x_1$:\n",
      "$$\n",
      "   x_1 = f_1 - (a_{12} x_2 + \\ldots + a_{1n} x_n)/a_{11},\n",
      "$$\n",
      "and then substitute this into the equations $2, \\ldots, n$. \n",
      "Then we eliminate $x_2$ and so on from the second equation. The important thing is that the **pivots** (that we divide over) are not equal to $0$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "## Backward step\n",
      "In the backward step, we solve equation for $x_n$, put it into the equation for $x_{n-1}$ and so on, until we \n",
      "compute all $x_i, i=1,\\ldots, n$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "## Complexity of the Gaussian elimination\n",
      "Each elimination step requires $\\mathcal{O}(n^2)$ operations. Thus, the cost of the naive algorithm is $\\mathcal{O}(n^3)$.  \n",
      "**Think a little bit**: Can Strassen help here? "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## A short demo\n",
      "Now we can do a short demo (with a Hilbert matrix)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "n = 500\n",
      "a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]\n",
      "a = np.array(a)\n",
      "rhs = np.ones(n) #Right-hand side\n",
      "x = np.linalg.solve(a, rhs)\n",
      "\n",
      "#And check if everything is fine\n",
      "er = np.linalg.norm(a.dot(x) - rhs) / np.linalg.norm(rhs)\n",
      "print er\n",
      "plt.xkcd()\n",
      "plt.plot(x)\n",
      "plt.title('Yea')"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "3.41412887051e-07\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<matplotlib.text.Text at 0x10fb80ed0>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAELCAYAAADTK53JAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsfXmYXEXV/tvTPTM9Mz37lgUkgQQIewATPllVZBNkk32R\nHVwAl0/lEz8NCsrvU9CwKKBsAio8ogIiBkGEBDHsCAFkCQmJWSaz79M9Pff3R+X0Pbe66t663fdm\nZjL3fZ486enl3rpVp96z1KlTMcuyLESIECFChK0OJePdgAgRIkSIEA4igo8QIUKErRQRwUeIECHC\nVoqI4CNEiBBhK0VE8BEiRIiwlSIi+AgRNBgbG8MjjzyCJUuWIEo2izAZERF8hK0WmUwG3/72t7Fi\nxQoAwOrVq/H1r38d5513HnbaaSeUlZU5/l1zzTW5346NjeGLX/wiPvOZz+CII47AzTffPF6PESFC\nwYhFefARJitGR0dxyCGH4IILLsA555yT9/lVV12F733ve3j77bdRXl6O+fPno7OzE9OmTUNzczOO\nPfZYVFVVAQASiQQuvPBC1NbWAgCuuOIKXH/99TjllFPw1FNPob+/H8uXL8dOO+20JR8xQoTiYEWI\nMEmRyWSsbbfd1komk9batWsdnz388MNWSUmJdeaZZ1qWZVkrV660EomEtXDhQmtgYMD1ul1dXVZ5\nebl10kknWZZlWa+88ooFwDrrrLPCeZAIEUJCqCGaZcuW4YUXXpAVCvr7+/Hqq6/i17/+de7fk08+\nCQBYsmQJ3nzzTeX1hoeHo1hohBwSiQQWL16M4eFhnH322bn3165di8suuwyNjY246aabAACzZ89G\nc3MzFixYgMrKStfrvv766xgZGcGxxx4LANhrr72w//7745577kF/f394DxQhQtAIU3tUV1dbAKw/\n/OEPlmVZ1po1a6yFCxdaACwA1syZM625c+dac+fOtS666CKrs7Mz99mLL76Yu86mTZusiy++2Cot\nLbX23ntv64MPPgiz2REmGU477TQLgPXzn//cymaz1qGHHmoBsB566CHH96ZPn25VVFRYp5xyinX6\n6afn/i1fvtzxvauvvtoCYD3zzDO59/75z39aAKy77rprizxThAhBIFSCP+OMM3JEnk6nrVtuucUC\nYF1zzTXW22+/bY2Ojjq+PzY2Zh1wwAEWAGvBggW59+fPn28BsGpra62zzjrLqq2ttVasWBFm0yNM\nIqxfv96aOXOmVVtbax122GEWAOuAAw6wxsbGct8ZHBy0ampqrNraWmuXXXaxAFiJRMI66KCD8gyG\nBQsWWLFYzOrs7LRGRkasO++80zrvvPMsANYtt9yyhZ8uQoTCESrBP/nkkzmLvLe31/rXv/5lAbB6\nenq0v7njjjssAFZ9fb1lWcJyqqystFauXGmtWrXKGhsbs6644gpru+22y1MQEaYu7r777pyszZs3\nz9qwYYPjc5LFO+64wxobG7Peeecda82aNcpr7bbbbhYAq7m52Zo+fXruugCsY489dks8ToQIgSDU\nGPyBBx6IlpYWAMAHH3yAWCwGANhvv/3wiU98Ap/4xCdw1FFHoa+vL/ebT3/600gmk+jv70dbWxtu\nvvlmHHfccZg9eza22247xGIxnHbaaVi9ejWy2WyYzY8wiXDqqadi2223xezZs/HUU0+htbXV8bm1\nee1m/vz5iMVimDt3LrbZZhvXa3Z3d2OHHXbA4sWL8fDDD6OkpATr168P7RkiRAgaiTAvXlpaip13\n3hltbW144403kEiI25188sl46aWXEI/HccYZZ6C6ujr3m5aWFrS2tmL16tV499138be//Q033nij\n47rLli3D/vvvj7KysjCbH2ESoaysDE1NTejo6MgjdwB4+umnAQjZaWlpwfPPP+/4/OCDD0Z9fb3j\n7yeeeCInswBw+OGH47HHHsP69esxffr0kJ4kQoTgECrBy3j//fdxwAEHYNGiRca/sSwrZ/kDwJtv\nvokrr7wSDz/8cAgtjDDZ8eGHHyrfX7NmDQDg0ksvxWWXXebIxpo2bRqWL1+eI/hsNouFCxc6yB0A\nPvaxj+Gxxx5DJpMJqfURIgQL3wTf3d2N9vb23N9r167F+++/j1NOOQWpVMr1tytXrsTIyAhGRkaw\nfv16jI2NAQAqKysxbdo05W9mz56NX/3qV9htt93w2muv4Wtf+xoWLVqEAw880G/TI0wB1NTUKN8/\n6KCDsGLFChx11FEoKRGRyX322Qfz5s1DY2NjboPTe++9h7feeivPawSAE088Ed/5znccBkeECBMZ\nvgi+s7MTu+yyCzZu3Jj32cqVKx1bvVVYsWIFXnjhBUyfPh1dXV259y+88ELcdtttyt/ccMMNOOmk\nkzB37lzMmjULv/vd77Dvvvv6aXaEKYJ9990XZ555pvKzc889F+eee67nNaqrqzFr1ixH2JAwb948\nPPjgg5g5c2bRbY0QYUvAV6mCsbEx/P73v8fg4CDmzJmD999/H2effTb23XdfPP74444YJuF73/se\nXn/9dVx77bVYsmQJNmzYgO233x4AsOuuu2L77bdHfX19zqoCgMsuuwzt7e248cYb0djYiEwmg6VL\nl+LAAw9EaWlpAI8dIUKECFs/Cq5Fs2TJEpx33nnYsGED1q9fn8uW2ZqQyWTQ39+PoaEh9PX1YWBg\nAIODg+jq6kJPTw96e3vR1dWF3t5eDA0NYWhoCOl0GsPDwxgZGUE6nUYmk0E2m82Fo0pKSlBaWopk\nMomqqiqUl5ejrKwMNTU1qKmpQWVlJaqqqlBbW5v7V1dXh6qqqtxnlZWVW12YwLKsXH8ODAygv78/\n1+cDAwMYGhpCb28vuru7c5/19/cjnU5jdHQURx55JI488kiMjIzg7rvvxrJly3LXpr6KxWKIx+Mo\nLS1FIpFAIpFAWVkZKisrUVFRgdbWVpxwwgnYfvvtMTw8jDfffBPZbBZ1dXWoq6tDS0sLksnkpO97\ny7LQ2dmJjo4O9Pb2ore3Fz09Pbnwa3d3N/r6+jA4OJiTZ5Lt0dHRnCwTdZSUlCAej6OkpCTXp6Wl\npY5/FRUVqK6uRn19PRobG3NyXF9fj4aGBqRSKVRVVaG6unqrNuIymQw2bdqE9evXo6urC+3t7di0\naVOOW4aGhvCjH/0osPsVTPAHHXQQli5dCkBY4vfccw/mz5/v+M5hhx2GxsbG3ABXVVWhsbERra2t\naGhoQE1NDZLJJGpra9HS0oL6+vrAM2Msy8oRNAlye3s7NmzYgM7OTvT19aG7uxsbN27Exo0b0d/f\nn+t0nr7phrKyMlRUVKCiogJlZWVIJpNIJpMoKytDIpFAPB5HPB6HZVkYGxtDJpPB8PAwBgcHc2sS\nvb29GB4eNn6uqqoqNDQ0oKmpKTdpmpqaMGPGDDQ0NKCuri73fzKZRHV1NVKpFJLJJFKpFOLxeKFd\nqkQ6nUZ3dzc6OjrQ09ODwcFBDAwMoL29HT09Pejv78emTZtyBN3f34+Ojg50dHTklKSpKFZWVuae\nR+5j/poTMfV9NptFJpPB6OgoRkdHkU6nMTg4iOHhYfT392NkZMT13hUVFWhqakJTUxOam5tRV1eH\n5uZmpFIp1NfX5+S4vr4eNTU1SKVSqKmpCVy2LcvKyU5fXx86OjrQ1taGdevW5Yib+ritrQ3t7e05\n5dnR0eH6nPF4HKlUCpWVlTnZ5vJcUlLi8Liz2WzOiKE+pT7OZDJIp9MYGhpCf3+/UWpzbW0tmpub\nc7JdV1eHmpoaVFVVoaamBs3NzWhsbMwpjIaGBlRVVaGqqgqpVCpUBZHNZnNKkQh5YGAAPT09OaNj\nw4YN2LBhA3p6etDT04NNmzblxsCLU6qrq/GVr3wFixYtCsSQKJjgBwYGMDIygg8//BDHHXcchoaG\nsHbt2lzn3nbbbbj44ouxZs0a/OQnP8G9996bs8bckEgkcsKVTCZRXl6O0tLSvIkci8Vyk5YEi4SJ\ntOHIyIiRUKVSKbS0tGDatGlIpVJobGzEtGnTcoJVUVGRszCqqqpQV1eH2tra3MQtLy8vpAvzkMlk\ncsRIwsGtWpokZOHSpCWS3Lhxo2Ntww00YcvKylBeXo5kMpmztugzUgKW2BCXm8hk0ZFwU197obq6\nOmetESESSZK3QpO1uroa1dXVOa+FLMDa2trAlRPH6OhozhDo7e3Fpk2b0NPTk+vnTZs25Uizo6Mj\nZwn39fV5ZtcQWaZSKVRUVCAej+feKy0tRUlJSW5SkzIiuabxT6fTSKfT6O/vx+joqOv9yEJuaWlB\nU1NTrm8bGxsxc+ZMNDU1obq6GjU1NTlPsampCalUKhQvhby0zs7OnNx0dnaiq6sL/f39GBwczBlg\nbW1tuc+6u7tzCmtwcNDzPpWVlbk+JvkmT5m4g/5xDiEllc1mkU6nMTIyguHh4dy/oaEhzz4n0MJ9\nTU0Nmpqa0NramhsD4prGxsbc32R83XnnnTj//POxZs0az30aJvCdRbN27Vpcc801OPbYY7Hzzjvn\nCPTjH/+4I62MshJ6enpw3XXX4brrrgMADA4Oor29PTcpRkZG0NXVlWfdkVU1MjKSI+7R0VFHuCMW\ni+UsikQikVME5P6Vl5fnLKiampqcMFOHNzc3o6qqKi8dLmgkEglks1nEYrFc21UoLS3NEd2MGTMK\nutfIyAi6u7tz7nZPT0/OOu3r68u9JgXNBTmTyeS8i3Q67QgrAUB5eTkSiQTKy8tzHgv1dU1NDerq\n6nLWK00y+pus7YmORCKBhoYGNDQ0+P4tWdJETERKfX19jlAekTVZu0NDQ8hkMg7ZoH4muSYlR0qZ\n5JpknIhixowZW0yu/SIWi+Xkm+PWW29FKpVCU1MTjj/+eNdrZDKZnOfX19eHzs7OHG+Q4UN/k0Ik\nTyedTiObzTo8Rc4hFGaKx+M5hUDeeDKZzHkypBDJCE2lUqitrc0ZLs3NzUgmk9pnWLFiBZ588km8\n/PLL6OzsRFtbG/r6+nDJJZfkylf39PQEQvC+LfhVq1bhzDPPxIsvvoiRkRHEYjFccMEF+NnPfuYQ\nqL/+9a847LDD8PTTT+Oggw4quqGTGXKoYCpj9uzZOSU+NDQ03s0ZN6xbty4Xr66rqxvv5owrovlh\nI2je9K3iZ82ahWXLlmHt2rV48sknccQRRyh3DtKia1tbW9GNjLD1YNWqVePdhAkBnmo51Uktgo2g\nebNgH26bbbbB5z73Oe3njY2NAETufIQIESJE8EbQvBlasTHaKBIdkBAhQoQIZgiaN0MjeMosMcmu\niBAhQoQIwfNmqASfSCTQ29sb1i0iRIgQYatC0LwZWh5VLBZDKpWKQjQAvvvd7453EyYMor4QiPrB\nRtQXNoLmzYI3OpmgtbUVxx13HG699dawbhEhQoStEJYFTPKKEAUjSN4M9USnsrIypNPpMG8RIUKE\nrQy33w7U1QHXXAMYbhzdqhAkb4ZK8BUVFVN6M0uECBH8489/Bnp7gW9/GzjvvPFuzZZHkLwZKsFX\nVVV51p6JECFCBA5uvL766vi1Y7wQJG+GSvClpaXR8WYRIkTwBV5U1eOQuK0SQfJmqASfSCSMq69F\niBAhAgBw41WqSzYlECRvhkrw8XjcqP5zhAgRIhB4hqDmiN2tGkHyZqgEX1JSEhVSihAhgi/wGHxz\n8/i1Y7wQJG+GSvARIkSI4Bec4As8FiHCZoRK8JH1HiFCBL/gBN/UNH7tGC8EyZuhEvzY2Jjj7MYI\nESJE8AIPP1dUjF87xgtB8mao7JvJZLbqE9IjRIgQPPiplpur504pBMmboRL80NAQKqaiCo4QIUIg\nKOBo3EmPIHkzVIJPp9OT4qDlCBEiTExMxRh8kLwZKsFns9kJd7J7hAgRJjbicfv15hPsphSC5M0o\nRBMhQoQJBR5+noox+EkTounv70dVVVWYt4gQIcJWBp4lOBXtwyB5M7LgI0SIMKHAi41NxUM/JoUF\nn06nkU6nUTMVi0lEiBChYEzlM4KC5s3QCL6rqwsAUF9fH9YtIkSIsBViKhegDZo3QyN4OhW8eiqu\nkkSIEKFgTGWCD5o3QyN4OhU8NRUr9keIEKFgTOUQTdC8GboFH8XgI0SI4AdTuUZh0LwZxeAjRIgQ\nYYJg0sTgqaENU7GYRIQIESIUgKB5MzSCb2trAwA0TsW9xhEiRCgaU7EQbdC8GaoFX1paGmXRRIgQ\noSBMxTJWQfNmqFk0EblHiBAhgjmC5s3QCH7Tpk1ReCZChAgFYypWGg+aN0Mj+O7u7miBNUKECBF8\nIGjeLIrgR0ZGcMcdd+Af//hH3meRBR8hQoRi0Nc33i3Y8pgwFvz777+PI488Eueffz6OOeYY9Emj\n0d7ejpaWlqIbGCFChKkFqpTLz2adKgiaNwsi+FWrVuGQQw5BW1sbrrrqKnR2dmLTpk25zy3LQnt7\ne2TBR4gQwQif+Qxw4IHAc89NzRrwQDi8WRDBr169GtOnT8fy5ctxzDHHAADKy8tznw8NDWFoaAjN\nzc3BtDJChAhbLbJZ4JFHgGXLgI99DGhvH+8WjQ/C4M2CCP7ggw/G888/j8rKSlxzzTU466yzMHPm\nzNznZM03NjbivvvuQywWU/5btGhRIA8RIUKEyQt+ButUwKJFi/K48JxzznHwZlAoaivB3//+dzzy\nyCNYs2aN4/3u7m4Aop7CunXrirlFhAgRpgDq6oDNtOGAZU2NU51qa2sdvBkUCl5k7e3txaWXXoqL\nLroob1GAGlpbW4uenp7iWhhhwuLXvwb+9a/xbkWErQG6uDs/vm9rRn19vYM3g0LMsgorznnHHXfg\n/PPPVx4Q++CDD+Kzn/0sXn31Vey5556BNDTCxMKLLwIf/ah4PZXLu0YIBjvsAKxcmf/+pk1AU9OW\nb894IAzeLMiC7+npwUUXXQQAmDNnDj7ykY9g4cKFyGazAJBLmYxqwW+96OgY7xZE2Jqg27W6uTz6\nlEAYvFlQDD6VSuGyyy5DPB7HwoUL8d5776G7uxslJUJfUFgmSFcjwsTCvHn266kSJ40QHnRh580H\nHE0JhMGbBRF8PB7H9ddfr/08Os1p6wdLmkJHx9RxoyOEg2nT1O9Ppd2sYfBmKLVo+vv7UV5ejsRU\nrPc5RcBT2957b/zaEWHrgO4I0qlE8GHwZmgEHx22PXWwYcN4tyDCZIeOLgYGtmw7xhNh8GYoBD80\nNIRkMhnGpSNMQEzVnYcRgoOuBPrQ0JZtx3giDN4MheAHBwdRWVkZxqUjTEBI+9wiRPAN3bpiOr1l\n2zGeCIM3QyH4dDqNsqlYrX+KYvMxkhEiFAxdFs3g4JZtx3giDN4MheBHR0ejBdYphPXrx7sFESY7\ndBa8CcF3dwM33wysXh1sm7Y0wuDNUFg4Iviphck+sSKMP3QFFL3y4J94AvjUp8TrV14BfvnLYNu1\nJREGb4ZiwWcyGZSWloZx6QgTEH4s+CefBObMAb7znfDaE2HyQZf67ZVFc/nl9mtVsbLJhDB4MxSC\nz2aziAdcA9SyoponExWmWTTLlgGHHgq8/z5wzz3htinC5EKhBM+LlE32ssNh8GYoBG9ZVq5sQRAY\nGxMHAXz841Nr0WWig1J2N5cg8sRf/2q/3nnn4NsTYfJCd850V5f773gh25GR4NozHgiaN4GQCB4A\nYgEWJ+nuBv75T+Dpp4HojJCJA78Hz3BrLCptEIFDZ8F77WTlimFrSKkMkjeBEC34IME18xtvBHrp\nCEXAL8HzBbPouN4IHMmkumCdV4iGb5AaHQ22TVsaQfMmEKIFHyT4braoQOXEgW73oQ788AZdcakI\nUxOxmDoX3iuLhlv+AUc3tgqE1iVBaiN+KJRfUokQHvyOBXehI0UdQYZK6XvVg+elW7aGJIygrfhQ\nCD4ej+cO/wgC/HCJAI8rjFAk/I5FJmO/nmqKemgIuOuuqDCbG1Rx+I0b3X/T2mq/nuxZNEHzJjAJ\nCV632r414cMPgd//3jw7Zbzgl6THxuzXfGJOBTzwAHDuucD06cDrr493ayYmVIUUvRZZ6+rs1xHB\n5yMUgk8kEhgNcMWDu2lbexXiW24BttsOOPFE4G9/G+/WuMNv4TvufdLEHB11Ev/WCp7ud8UV49eO\niYxCzrngv5nsm+eD5k1gkljwnOC5xt7a8P77wOc/b/890S14v3WROMG/9x5w1FFAaSlw/vnBtmsi\ngmeIeIUdpio4WZum0W6zjf16shP8pLHgg24od9O25vxpeVPHRD/x0O+E4gR/+unAY4+J11TqYGRE\nlDA4/PDJv+1cBu+rqXTOqB/wkJ9pKJbzwWSvjhIGwYei84J2NXia5NYcg5d36U70XHG/MU/d4Q3H\nHCNSKA84AHjpJfHeK6+InctbC7gF77U7M0hks5MnNs3LDkybBrzzjvdvuBE02SuUT5oQTVlZGdIB\nbivj+dNbc4hGVt4TneD9brrTLZgdfDCwcqVN7sDE9178gq9XdHVtmZS+TZvEYvaCBcCLL4Z/v2LB\n89h5CQI3cKXAs7QmI4LmTSAkgk+lUugP0A/llu3WnCYpLzZOdGXmd3FUtysxlcpXbjNmFNamiQpe\ncTOT2TLb6j/4QGSgvfCCIPnHHw//nsWAE7xpMgU3Mjo7g23PlkbQvAmERPDV1dXoC/A4dE4MW5tl\nx8FJLpWa+ItGfi0m3clPZWVOwovHzS24YrF6NXDttd4baorBH/8IfPvbzve2RGEsHpO2LOD558O/\nZzHgIbxCTq6b7MXGguZNIKQYfCqVwuDgYGDlL3k52skeZ3MDt4gnw0YgvxNKZ2GVlzsJfsaMLRc3\n3nlnEQIcHg6vkN0Pf5j/3uBg+MaKbCBM9EqsXAaqquzXzz4LzJoFzJzp/vuJnnXmhaB5EwjJgq/b\nHFvo4TUGisBUTCubDPn+fsMM9P3ycieBJ5POybklw3C0vhPmubKqTV26BecgIRM8X8uaiODt4/J/\nwAEiHfLWW91/P9lLFQTNm0BIBF+z2TQJyt3Y2lLmTDDR4+9A4RZhMilInv/NvZfxePbtt9+y194S\nC4KytzvRQxibNtmvVSGaSy4BXntN//vJXmwsaN4EQiL4is1L20MBmSlbMq1somBLxaCLQaEWoeyd\nyN7oeOx12G678K49e3b+e1uCbHmGCTDx66VzXtPtkl65Uv/7yZ4HHzRvAiERfNXmANqAVzFnQwS8\n7jBhwTMCJkM53UIXJuPx/Hgpd6+31LPze37kI+HdZ7wseLli50RPI+SF2HQE78Z9EzXfv60NOO88\nUV/KDUHzJhAywQ8GvKrDF162RnAXczKkCcqhwu5uYPFip6utQjyen0PPyTbMcAnHltohPWdO/nuX\nXw48+GC4cWM5zDHeFvyHHwLnnAPMnQv8+9/5n/O1Nt0eC7fw9ERdY3joIeDOO0V9qdNP16cXh8Gb\noaVJAsHGkoCtO4MGcFogfk9LGg/Ih22fdhrw5S8DZ5zh/rt4PH8ycuU2fXow7ZOxfDnwmc8A994r\n/uYeSJhZSyqCX7YM+OxngSOOEKT1r38B990X7H1jMadMjeeJR3//O7DjjsDdd4s6RKqUTZ5lpfMO\nJ0qZh7ExUTvKREHzNcTf/AZYt079vTB4M7SNTkDwBD/ZY2xe4FkPkyEGzzfvAMBf/iL+X7bM/Xcq\nV5oTfBg7eK+9FthvP+CRR4Bf/lK8x8mi2IXd9euBL30J2Gsv4N13nZ+VlgLXXaf+3eOPA/ffD+y5\nJ3DmmWbb8/2Ax+HHM43wueec6w5um+Sqq/UL+DLBc4LdkovIV1whFLeJUpbTg3XeSRi8OSli8ISJ\nvvGnWHAPxdSCf+st4JprxifHmSx4OXvBa5MKJxoKu/GxDSMGf8899muy1rklVYx3+OSTwA47ADff\nLLI8/vnP/O/stZf+97/7nf066HRNXrtpPAlezvl38ya22UafOSfH4Pn3tuTz/ehH4v/rr/f+rvws\nPIOMY9LE4Os3JzJ3Brx3eKJb8E8/DXz1q8B//lPY7/nCEi+DqsPGjcAuu4hdkuO5DV2OmXspJz65\nyZILm+B5SGiHHcT/q1YFc+2//tVJPHL2CuDuyvMJH/RC6EQp7SETvFt/bLONPtYu9w/3IscjndqE\nk+Rn0XkvYfBm0QTf39+fd45gTU0NYrFYoAn7gF7zjTc6O0Vt80MOAX7yEzvG6xec4E1I7v777dfj\nGV/daSfn314LltxapjHlQh/GgicnYIrxBxUOkclKlQHiRmhhWtk8JdW0ONxPfiIWgVevDq4dKqWn\nQ0uLnqxlOecEvyU2j8l46y3v75ga5GHwZlEE/9JLL6G6uhpHHnmko45xLBZDTU1N4AQ/UdOgHn3U\nrm0OFE62PEvIZCcrV/TjuTFq222df3sdx6c6m5XvdQhjnPmYUF56UOejyiEqFZm5ETe3boPOBOEy\nZeoBf+tbwA03iPIA3/++/nv/+Y95e+UMOLdNSdOm6TOxPvzQ+feaNfbrLRmmJC/VJFwu06Au1BwG\nbxZF8M3NzUgkEliyZAmuk1aR6uvr0R2wzzRRtyLLE8dvGV0CXyQy2ZXHLZbxrF0jW9xebeEET30X\n9uTkHgKFaIKKd8sTVkXwbqEXvmYRNMH7PdLOspxt+MEP1N/7z3/E5rDp0913lxJkgndry5w5+rpF\nDz8s0ksJXBH09PjjiNdeA1591fz7HCYhVIJswX/rW8DatervBs2bRRH8Rz7yERx77LEAgCVLljg+\na2xsxCavhGifmKh5rrJLXqgi4hkCJtfgpDieVTZnzXL+LW+wkcHzsTMZ4WbLmSdhgjKUdOlqfiEr\neBXBu4UPuPLRTfxCwYnVy2jIZvMXg4eH1bLY1ia+390twpNeXqssn27exKxZ7qV/f/Ur+zU3drNZ\nc0NhcFA86/z5Zt+X4WevhtymW28VXu8dd+R/N2jeLDoGf/LJJwNAXh3jxsZGdAVcY2Ci7mgNynrm\nitsk33c86+Rzi3TbbZ0k5ZWRwn+byYiF4m9+0/03r78OHHoocOCBhR3SzX9D7nVQ80gmK1UWkZvs\n8vBNkHFvwN+JR5SPL0MVXuJKbN069w1IQP4ccVtP+8hH3HdJu9WAN02VfPZZs+/p4MeCl5W7ZYl/\n552X/93cAgkpAAAgAElEQVSgebNogv/kJz8JAHj++ecdO7B48frm5mbEYrG8f4t81medKJscZMjC\nWmiIhhO8SRkAbglvaQuety+VEjv1CF47JrknNjzsnf3w6KPC0nrySZFjX0hohc8Z8riCCgvJxPm7\n3+UrITfZ5RayW62VQsDXZnTb/wk6z0tF8LK8efWlTPByn/E+aG11V4j8OWSCN/XyTcJKbvCT6SUT\n/KJFi5R82NnZGfihH0UTfInG7+MNTRVZ+5YIdKJWw/OzgOQGTpom5OOVmhcm+AT81reACy6w//ZS\nTrzdKtdeDgm88IKTZIrJGOJFxYI65EMmq29/G/jud53vuYVo+PMWmmKrA1+w9lq81n2uCtHIU9qr\nL2WFICsbXqaALFwdePxeNg5UBG9ZwEkniVASbcJTlUrwAy5Hjz4KHHSQqHapgqnS6evrm3gEr0Nd\nXV0un7OhyJOyw9q6HhTkDJZCs0B4yMBkjHked6FeQ6Hg7eOZDIA+ZKVqoyqFTFbksvgUQ/C8bABZ\n2cXW3lcp16uvdv4tT3JOeNwgCPrYOR4+KnQfiYps5TCUVypgMunsZ5nwP/jAfu0VoeAKVSZ41WJ2\ne7vwqh57TIT4vvAFUWagGPB1p6OPBpYuBX7xC/V3TQm+p6fHwZtBIDSCb25uRl9fHzKZDF566SVY\nlpX3zzREM9ErK8oWfKEEz2u7mFiXQS0SFgI+oeXCaLqdrCqiUL0nW7uywiiG4FWle4sNb5kUwZND\nDrrTieT6PsWCk2Ehx+DpkEg4QzomHidPn5X7jMuyvDYiGwZcUcn9qgonyXz5u9/lGyU6tLWJA0cu\nvdQZGlQZnbp1MDlkedVViwBYOPRQJx/uscceDt4MAkUTfDweR2VlJbbbbjuUsp4Psq7CRCd4eeIU\nWlKBd5WJEi8kjzubFQuaPBOhEHBhlz0Y1QKan8wimSxkgi9G9uWcfcA7Nu0Fk9/zZ5LJkRsEQa8z\nce+i0DCebuz4uJhsMuJyIreFy7Ks5OQQGJ9fJIfk5f35z/n3pevttpv4f2DA6TG44bnnxILsTTeJ\n9Nqnnxbvq7w+VQ2lbFafFKAyLIKuR1M0wdfU1GBgYACrVq1yEHztZgkOYkV4PA6A8IOqKqcQFlrX\nhLubJuRNE88rLZHj9tuB//s/4HOf89c2GTzbQ1ZwxVqKsvciT6Ziyt6qDvYodmOVjjg56fGxLStz\nhp1iMbvPgt7rwS3lQrO9dOE/LudyWE21XsbnsTym3IKXF9FleaJnGh0V5B2L2cbR17+ef18ieFLu\nQ0O2pe9VVoPLWn+/XVBPJTMquXfzNlUEHyRvAiGGaJo2j2YQDZ3oFnws5hSUQksq8FQzP666H4J/\n5BHz77qBT0JZUFVE4meNQH52+fkKSZMkhCFLOoLnSpCHHaqrnfJSUuJvDP2Ahw0KXQrTGSycpInQ\nV68GTjlF9MkNNzi/z8NSMsHzkgNyiEYOfRCRErXU17vLBHnDzc3iWbgS9TrJS1ZU5D2oEin8ruWo\nQjpB8iYQIsFTZTQ/K8Lr1qmtM6/T1CcCOKkV6vLzrvITfvETQzYtsHXvvaL8rS6FkWc9yPfXxaRN\nQ1dyeEq2svwSPJ/QqklVjMIA9JYx7yPej01Nznh0SUl4O5G5VWlC8PJ6ilxTnoMrtuFhIb977gk8\n8IDoczlTxS1E09FhvyYPjgwlWfnRM5EiaGpyKm55UZOu3dCQL5teJ3nJ1yIZVhksfsuFVFUJo27v\nvcXZAJs2FcabbgiN4Cs3j4Lp6SQvvyw2D2y7rViN5q7NRLfgAaf7WagFzxcu/ShwP/XTTXeMfvOb\novztnDnqnGRuwcuWi07ByQSv+558P5mU/aahcs9IRVbFFvji5MwJhJM6f6aGBuciXTwe7AIoB7e+\nTQwBeQOPmyzzcRkeFvLL+1quScT7XiZIntVC3EZ9KY839TcRd3OzU5bk3cBkMDQ05F+LZ8Ok0/kh\nFdmCp/5QyZEq3OMW/isrExvLXnlFlF+4/nr/vOmF0Ajebyzpgw+E1m9rAy66SORWE6h2CDBx69Hw\niVHoYha3FvhEeeghsb2eFnhk+HHvSWC9rDmKiXZ0iHK4Mnj4QbaKdGMkE/wee6i/J8fg6+qcv/Wb\nEsrbqpo3xR5lx5+fW4TcKuVkX1Gx5QieE7SJXMresltqpUzwsgUrE55bCOPtt+3XMsHLbSBZoO+l\nUs7n5P0O2PJUW5svmzRev/2tMJSOOcb5ubygT4pEJYOqtUI3WS0rc8ped/ckisE3bjYrTRsqd/x7\n79mvOXkGnSccFLiXUVUlBGaffYAf/9j8Gtxa4AR/3HHCfbvySvXvCknzc6v4ODbmFEyVxfzGG/r7\n60IessVuSvDxuJMQ/VrwPL6rIvhiaxxx4tpnH/s1JxpuwVdV5ZcQCOs4Sq441q0D/vEP9+/LcuGW\nAsoNi5ERQbL8ueRTydwMEfKiZs60x4iIXFYMRPg0R6qrnfwhW/DU99XV+bLT0iKOEzztNOciKkEm\neFKSKq9PZTS5yWpFhZPgOzr886YXQjsjqW6zOjetjEaEUlIiCIJ3IBeMtrZwjnQrFrxNp57qJJX/\n/m+za/DBVi2y6nbyFrJA53aod3e3c8KoJjkJfm1tPsHrQh6yJTZvnvp7qmefNs3OXfab9cLniiqd\nr9Ad0mNjok3cAKFUPMD2gsbGnIqlsdE5ZlVVxadq6sCv+5nPCEXS06O/n6xo3DwLLvMU2pg+3VbQ\nMsGbhC7574ngZc+ACJ8vsnIi5sYh4CR4WZnPmOE8H3bvvZ2fyyEbemaVEaMK0dAahmpOVFU5293T\n4583vRCaBV9WVoby8nL0euzYeest4OCD7WPOSOPxjuVakC9cbQm8/rr3GaOAUwg5uftZeOGDrVLg\nsuASCtmJ6bauIcfA3dJUW1vN89TlEI1q0xGgXmDm7fVL8FwEVQRfyEERmQywcKGI4b70kv0+z8p4\n4AFBBKqNW7xPk0ln39xyS3Cb2OQQQTrtblWalD4mcPLnFjhB9gZMQmv19fbiKcm1vAZDJEsWfF2d\n0ziSzwom5VpVlT8WM2c6kxu8vGFSzCrC1hUg0z13U1O+127Km6YIjeABsWDgtViw667AM8/YoQwS\nPp2bLxf8DwsffACceKIIIxx0kHfaoi4L4m9/E+6fSYEsToxtbfbgk0DrlHohGRhumUnykLlZcTvs\nkD8pdAQvW3C6Nqj6ik8evxvJeKhkZETI0GGH2e/JHqMXrrlGWLovvij+5jF+Ttzr1gG77w5IlbTR\n0OC04JNJ5/0//3ngyCPN2+MG1fRzCwfJn7l9l487eXx8TGULnhMd9xDlcw0oDEtyLXvsJEfEgTU1\nTr6QjUCSl97e/HEuK3Na9V5zieaCKr9dfl4vNDU55wr1gwlvmiJUgq+qqvI8QFaOvdNA6bTeliL4\niy4Cfv978dqyvGttuLmfv/2tqD3tFVYjoaFnJ0HffXf33xVC8G5hHdnKcdsMsv32+dfSLVrKZKHa\nVQqoy/jyGLzfmip8PWPjRjEW8sKx6SE6b7whiolx8Nrgshy8+aYwFDiqq52WcU1N/g5WVdneQiDz\nhFd4Ux4jt9AR/4zmMVdw8rU4CXNS5YQci9ljQcaebrGW5kd9vVPm5F2q5PFRiI+3OxZzkr5X8gF5\nEyojRkfwbhY8v7e9uOzNm6YIleArKiow5NP/JZLTZWKYbjEuFjJJFVp+gLB+vVjMcQNNABIU8hp2\n2SX/u14xchmjoyJkQHAL0XC3MR7Pd5G5cG+3XT5p6AheVkS6ySRnQQBOJeJ3QZKT99q1akVruqal\nIl5O1iZhCDnro6wsvFLYsqfvlfctK083Q4DLHT03t+CJoDdtEgYSbwt/Xu4dj46KsYjF7LbIckKe\nA/2uudl5PdkIJKubQpxyoTcuz6p9ABwk6yoZ180pnUy0tDg9AaLKQnhThwlH8F4w3ahTLOSwRBBZ\nDl7lJYi0SYiIdFSnx3Ch9EodffppEQo75RT7PbcDCzjBJ5P5AsonZGVlfoxet2gp92l5udrzUVnw\nnGiIUDs7xZFyf/ub+n4ETt66tpluLFONIZ+ksiegGrvqaieRlpebH8zsF/L08wojyGPtFp5T1dPh\nXlkiIcJtLS2inj/PgOOEzImfjt2rqbFJVLfRieSwsdE5H+QQDVnsdH/+jPG4UyZkWZbnPSkHlRyp\n1oYyGb3BU1vr9GomHcEnk0mMeKQo8Bx3Dt1CkG6hMWjIVrFXWMAkQ8Br3YQGm8iXaoNzq4LIirve\nXrJw0UXAO+8433NbmOVDpnouTnKpVL4br0s7lC2xdFq9qDU8nH8N7iWkUsD994v3rrzSuWdCBd5e\nXWhTXpjTQWXpc1mVP1cp3xkznCRDmS1hQJYNrwPR5bUvt5PCuEdGHi4fp3jcXn94912nZc29NJ7F\nQv1XXW0rAVluaC7Sd2trnQQvPwO1k8aYz+WSEicB6+L9gFAsZFyYhsjdPLOSEqd8kMyb8KYpQiX4\neDyOUY/arnPn6n6rfn9LxeBlgtdZ8Nks8KMfASecYL+nc8lMJzFl3hAx8fAQhai4gHlZfyrCdVNI\nfLKo4vucxEpL859XR/CyNdjbqw8ByAvKnGgSCeDLX7b/9np+3le6hWpVWEgF1RhyBWISamloyO8z\nlWcQBOnLfaNa2B4ZEbsozzwzvx/cYtJchkhGeZgiFnPKEs8M4gr1f/7Hfk1Wdm2t/VomXZIjOVYv\nY+1aYdiQAUJjz7llbMzO4JPbz58LcCpHncxZFvDHPwLXXiue3c9aKcXjTXjTFKETfNYjPUGXjy0P\nmlfsMGjIFq6KEEdHRTbGN77hfL9YgiflQmTBr0ebOLiAeXkGqvbE4yIFVCWofMhUKZKcjFSejY7k\n5JTRwUF9CKCnR7jaZF3J3+NWmletIv6MnZ3qNQsK0YyNAQ8/rI/JqxQBDynJRK3qX1m2LEtNBEEc\nRC7LBl+sHhoSKZ7z5gFf+xpw333i4AoOt0VZbvSQHMg1YXgqLG+L7nBx+k5Fhd3Xcqom3Yv6WjU3\nv/Y1ES7aYw9byZBxyOXwuedEqjZB9lj4M3KuUsn422+L1NnjjxdK68UX3Y2PdeucGWOc4L140xSh\nEnzMYMVJF6KRLXidpS9jeNi57blQcMuytFTtUWzapI7/6lI8TUs8EwGQEHESpYnBCcHruqq2P/+8\nEP5ttsnP8/cqzsWFUkXQOkUmE2tXl57gly8XbZsxA/jNb/IXufkzee0DkMv2qvYmUKrjyScDxx4r\nYvsqqNaAdCUJALWikL0iHQn4leN168SmOh6ykq89bZog0QsuEP22777OxAU5pOS2gC8vFAPO9YhV\nq5wZMLrqmhwkO6WlNtnLMWwy/mhcVfJ9/fXi/5ER+3t0PW6J336783fyGgUneC8L/sYbxfGSBDm+\nL2OvvZy706nvTXjTFKETvOWxAqhahALyLXjd9zjefltYHPPmOS0RyxJlcg87zPygCx4b1pFQc7O/\nuiheoQR6ZhIycmm5AJOlycMgXmEBlZVNk7q7Gzj0UKeS4M+kCtHwRVZV6ErXHtka7O3Vk/Orrwqy\n6OgATj/dOW7ZrFM+vM4L4JNs/Xq10iLr7sEHxf+6hVtVjj4nePlzVTqdvO6gIwE/CQWPPy4s1uuu\nA374Q5toZWU7Z44Ixdx+u9oQkeXZbaMelytSwFyZrlzpJHy+yOqVtUTXqa1VrzFZlj0HvJIM5DUi\n7hHIPCOHaHTVOFUyLo9jdbU7wctKjhO8F2+aIlSCN4FpDN7Egj/rLNuypQEYGgIOOURs0/7rX4G7\n7jJrF7c8dHnmiYS/zQ1e8Th6Zro3ZQNw65Umxp/+ZL/36qvu11URvLz5hysMLvR/+APw8587BZUr\nKlX6qC7OLRNxd7d+p6QcWuDpiWLHn/23V7VRuaCTajzl++k2GsnE1NLifM/rKLj6ejHOnFx1601E\n8A8+CFx1lftmrL//3b5mPG6Pi7xBb84cp8W+007Oz1W1WnTg5CfHxQGRJKBrs9dOfFIGra1qgh8d\nFYQYj3vXEpKNEO6VeOX9c2OAK2aTjaa1tVv+rGQZoRJ8NptFiUdlKB4/5Z0rx9V23NF+rVNuK1fa\nr2nSv/uu2ClLMO1wbrm4bV/2ykrgMF1wIQIiJcWfl7Q+P9TZqwCbysqWrQdO1PKQfeELwp0ka0w+\nfk6eYLpyEioLXrfYK1uXvL2dnXrXWQXZila1r6vLqcR0ZRTkGHx9vbM/vHYsE2FwgtOFtNasAc47\nT9QKX7RIKHKV7Le1CaudcPzx9mvKxCJUVztJV66eKM8Pt01unChpvsinkun2jxBB6sKZ9Pk776jJ\nlJR2WZm3spDDTJxnvKp46s55MFlPq6vzR/Ck6E140xShEnw6nUaZRwI5d3u4NSd3PCdZHVHyMAMp\nDlkjmxbm4lah205RPwTvFUohYaCwBT0PnwQqAvHy5lQkKhMVnwQq2Xr7bdsa5M9x2235Vp5s5XZ0\niDbKY+pmwcvgz93e7nwmr5288vFs3BAgvPee0xPSLerLYbamJud7XhOf+oDn3euUwtKlwJ132n/v\nu69QPHfdZT/TXXc5a5oDwkoHhOXL66wTuDwdd5zzM1mWdN5FJiMOryZUVorUVR5K6+jQx/BJtnVh\nS96PqjRYUtqlpd71qX79a+ffXHa8Ngly/uDf9QoxVVaKf34iLUSVJrxpilAJfnh4GEmPMnk8fMA7\nzW01mzT6j38MHH64cDkzGXvQy8vtySzHeE3rbvMVc64U/vlPsTBCgumH4L1i8NQXRHp9fWKV/+ij\n7e+oFqe8ruu1EJpIOPtXN6mJ+Plv7703f5GXK4DFi4UCuPbafGXb3q5PP5WVOE+r6+pyXstLaZue\nv7nffvZrnQUvY+ZMf4elkzLiY6azQFWhidWrgXPPFXF0QBCq/L3Vq8Xeh+pq9c5vvidi332dn8lj\nr7Ke//AHEdq54AL7vY9/XFRRveMO+72NG/UERzKkq/HEf6cKQZI3mUh4h0vkPSBcXryMA07q3MhU\nETxvM93Dz9kQpHhMeNMUoZULBoBMJuM4iFsFnQsjh0X4ZdraRLrXd78riGD77Z0boGbMEGT/9NPi\nCDEO0w7n4QRODF/4gjiB5R//AJ56yn1HqAyZtNJpMUEXLxY7/eTDBN58E/jYx5y/6ejIj2t2dopr\n6chSZaVwkpaHSDdkzz0nMky8JtTAgLAS29rsfHWVFdbdrb+X3Ffc6rzqKuc4upU+7u21vRU/ZZV1\ndXJk1NY687u9wgU0iQs9npHw8svif5W39ZvfuP+WiKihQbSnpcX2IuSQimqR+IILzM5l+M9/8hcZ\na2rEmHR2iqqZvBKnH9AzlJSY72HgbSB4RUJ0h5ar7snHnkJWfmSO5q8Jb5pi3AleB1mBcfIiYeTC\nR24pIKyWmhohXFdcYdeYB8wteO4283g8TUaqK8Pv6wVuaY2NAQsWAK+9Jv5+5x07LOCW9tfbq45b\nbtgg2vzII8BllzkVlOqZuYBygR8b08fFjztOFNsyWUvgxOr2HVnhJpMipt/RkX/iDWH5cuffOjJe\nscJZNVKXmqeC6bGLVVVmC278+4B+q74XqH+IiIvhAZKzVMqWd9lIUHlzpidgtbfnK4jGRvG8770n\nqmYGAb8bwuTdtm7gBM9fq1KTeaioEIKn7wZJ8OMeotEhmRTkt3ixIGr+vOSuu2lfshyuvdb5vld9\ndssSVifP2uGTceFC5z1MLT3AuRg5MGCTOyA8Eppcbl7G4KA6M+GSS4CPfhT43vfsVD+CSsj4M/EC\nSnPniowjHXp7zeqmHH008J3vuH9nYCDfg6N9EW1t5spYN6aPPuq0rgs5GMUL1dX+6smrFln9nO3A\nvbuXXir8/F/Atki5lS1b3KoMFdPCe21t+ePrJ6RpAssy319C4ATvVT2SyyBXBiqlrLLgKyvNF1qp\nb4IM0YRK8CMjIygvUALr6kSs+8tfFuELLlQdHSIGZ1qugVu8brU1AJFKuXix8z1OaFyptLX5I3je\nXnmS1NU5N4vowi2WpX7uxx6zX/vJZQZsS+7DD8UCpJvlPWuW2RF3zz5rl1vWoa8vv60UbunrMy/w\npvueLHpe+fKFoK7OHg8TBXLvvWJxkodl/HgWXJkcfHBxB4YTufFrymOrstZNCZ4W1zmCOo2NnxtB\nHqUpkcpHDZqCP4sqBs/DVrRjuKTE3FAhbiqGN2WMexaNDjU19gLRW285txN3dBRefe8rXxGTjGPj\nRpGOdt996gnzu9+JRa333nMK0fr1zq3fXuCTRfbAqqttQkqn3d1gL4tRJnRTgjcpmdvaWvwh1YT+\nfv1Gk/7+wurcd3barrIc6grjqEduaHiVTCCcdJKzvEUhMXhAzIFidm2T7PLxlGVLNdZ+ogeyTAWR\n/ZdO22sQo6O2Ba/yfFU8yeXKj/VPfWFZ+emngNOC556KqedI8jppsmjS6XTBmqiiwrnAyONb3d2F\nHbNGOOssZ9zuT38S6Whnngl8+tPq39x3n6jLzsMfa9b4c/t5PDIedyqLujp74ngRqJfV0dIiBPcH\nPxCFj7y8FootmhC8XD/bDV7xzf7+/PWEpibRLwMDZnXuZSxYIJTEBRfkT94wzj0dHrZlsdDwQ6EE\nDxSnbMmC51zitrhN8EPwcnzcrb2m5P+Nb9hrK6Ojdv+rfq8ac/6en7Ot6fpdXern4JzkJwwk/6YY\n3pQRGsFns1kMDw+jqpBZivwT2vkOwSBOs/r3v+3X3C11y1uViW3dOn8bGTIZvh3ZaXHU1dnKwmvR\nyCtEsnatILkrrxS1SbxqtVA7TEMFpm6tlyvf2ZlPIKmUrZBMQkEcY2N23vfttwtvjaOQXYXPPOOs\nNihjYMCe2H4yqjiK2ZXulkHkdS8KHXBilJWiSpmbZKKRspNj1W5FEk0z3HgIdXTUvodK8ajCI9wo\n83PQCnn3KutdBid1UyNwxozieVNGaATfv7nnUoWcCA3hRvGO4fmyQ0P+Ymcq8HhZoXFM0xriHJzQ\nuEA3NdnPq1rAufFGm0B0u/8Iv/2trQT/8x9vS5raYSK4gLmCVREXt57a2vK/U1Fhu9BezynDy5o1\njYUShoeBT30K+K//EpuOVNfnBO9nPSYo+L2nfJgL4HwuOfTpt9Q0gdKT5QVkNyVrOj68LtXoqDtJ\nq3ahc4/WK9TLZZBkVZ73Ks+toUG06wc/MOcJ8ZvieFNGaAQ/vFkyvFaDddbLZz/rzG3nC3+Dg+Zu\nom5HIlcYXgQIqN0/r4O4VeAWDO+aadPsv1Xhpx13tCeHG/G1tDgtyVmzvJUhuYamgqgLj911F3DP\nPfbfKmuNx6kHB9UWPMm2qUVHCpqIKpVSTzqdUXTwwfYB2hwdHfY1r71Wrdh6e+3+8KqJ4wem5bH9\npOkCTjIkouNWuywrqrE2GRcqLSIbK27eiumaC19T4lk0qmur5jaXA68UVa7gSBbkeb/XXuo2XnWV\n8KJVO6dVaGoy501ThG7BV3uMmo583nnHSRY8VjY87B12IOgmAI/pm6xnqKwLr7ojKnDS45ZQS4t7\niKax0RZWtxj47NkihY7fwysfnYbINESjI/iaGmd/qxSRHFKQn6WhwXZvTUMXdLITeSBlZeqFZd04\nDwwA++yT/z4nvuFhNRlQmCkW817MdoNMmqYZP34teK6kSN7c+lk1P73mXmWl3Rey1+PmLZvGquXS\nDDRfVM+huh8PHfISxipwWSd5kOe9Shk3NYnNkH5QX2/Om6YIjeD7NveGl6thmg3Dtebjj+fHV3XQ\nxSh5TJ+I022RR2W1BB2ioUmtIuRp0+z2caKRPZnly53F1aqqvGP6JPAmCouXaZXR0pI/+WTIJQDk\na23YYJd65plTbujpEa7w7ruLv8vL1XFPnSjqFBYnw1hMnXFBchmPF5b1o4NJxk95ufcCugweMjFR\noCpjwis8XF9v97Vq97YOpllO223n/NutcJlXQoCqlAMHby/NI7m+j0rJ1tf7298AiLaa8qYptliI\n5tVXxQaanXd2EolpNoycacBrXrhBF37hFjx9xy22qIodrlvnHfeV769a5AKcIRoVgfLa81xwvEJV\nqZR3iIYsWy8LPhZzf96ZM90tz7Ky/M/lCcgPTPC7NkLfLy1Vx151fdXTk09EGzc6wxljY2pZJWVc\nUlIcwfNnjcfNrNn6ev+bt/hJSibGlSrM5qVUGhvzD60huI2paUhO9qZJJlVkPjamTmUmT1NWCm5p\notRf8rGhqihBY6NZOQeOp56aRCGa7s0sVFtbixdfFJuVHnlEZK/cf79NdKYE77ezCLpFHU7wJoKl\nIra1a71jpbJAc5LhSrqlxSYlVSigrMy+Fid4r0yVxx8H/t//c/8OkYkXwcuHG8toaRHt0S2WNTbm\nhzFMiluZYI897Nc6gtdZrG1t+WT3l7842zI2plaUJJexWHEEz+Wrvt67QiIgxsPvPTk5FbqXxMvS\nbm62jSW5z9yyaNzKcnPolITuABPZ4gf0nv3atcANN9jhPn5Nmnfy2KjWXsrL841Sr7U+YfXbvBkE\nQiP4zs2S39TUhNdfd3522WXAAQeIWu2mQiYTi6kw6CY1D4MQ6bqFaHTEYzIRObiXwMewsZGXC1X/\nliYLtzaD2DhCJOG1aNzS4u7Wk6LUWXjTp+e797Jy9ZsaSeCTuLxcHUbgxMCVazotxoITl2zB6zwh\nmvSxmLlMqsAVdV1dvjepWiNIpfxb8CtW2K8LVaZeIRq+K1seT7eNRaYWvGoc4nH9piyVBa/rtx13\nBC6/XNTXz2Sc0QbqLzluX1Njto7nZZTX1zt5MwgURQ8bN25EWsNGpInq6uqUGvcf/xCdWahlbrpd\nWmfBcyuYJpMbWfpN2VMhkXDeg18zldJPCvn7vEiT3zifCsmkWemHGTPc09zoM114oaEhPw4uX69Q\nq5ITfDyuJiFOaHLOemmpOKgZEGWDv/IV5/eTSXX/0JiMjblPYC8C4L+tqckvbazqU1XIywucnPr6\nvHmYWKMAACAASURBVOPwqvH28hpaWuznlT30sAheZ9VblprgvZTxzTcLeeclpDMZMU9kTzeZNEvx\n9Lont+DrilmxZyiI4NevX49TTjkF06dPx7x58/C6bKLD1kT19fWu5PjJTxbSAvNNK7qB5znYZMGH\nfbyWLMAymdGk1lnSQZUIkFFaapbbXl9vlp6qC1s1NOT/XiZFU8UtQ86N9qqBn8mISqOVleL8VT72\nBx8s2skncizmHp7KZt0JyosA5E1vMhmoiDyV8r87ly9cd3V5h0hV4+FlwfPCefL13fowFjPLsfez\nk7a7W60cTe4jz8NMRqy7yZxSVmYWKvM63rOhwcmbQaAggj/00EPxwAMPoL6+Hsceeyw++tGPYrlU\nw3XdunWor69HWVlZwdavW8zKlAhk4nroITGZhodFOuHxxwMPPCA+K6b8gQnkySjHI1WVBjl0FjZX\n9oUQZDJpZjnX1JhdX1efZ8aM/DGVlWqhISdOkLGYTaicDLgcjo0JN7yjQxxWAfAj08T/PIxnQvBu\nxOeVFCEffyf3kypmXFvrP8zC95Z0d3t7gCoy9ZIB3n55/rnF4Pv7zTwSP3wyMKDuez8HcRDkkA1B\nlc3EjTHqLy8lMHOmkzeDgG86WLZsGT744AN0bJb+hoYGzJgxA0cffTTWrVuXq2Pc3t6O1s27Tfzm\ngxKmTdPvrjTZnAQ4Beq220QmT2OjWEzZbTfnd01qrFRW+i+VcOihwBNP5E8WeaLQmOpCAbrQTVOT\nPVFTKf9hm7o6M4JPpcwml67oVkWFN8EXWoaCXzeRsAkzmbTH1atWCa9QaFnOdNOREXdyAtytQi+C\n59emg7k55JOXANH2QjbbEXp7vRfWVTzjZQg1NOgJ3m2NZe1aITteO6rlsNKppwJf/aoIqz37rPOz\nbFbd94UYyOm0muDjceeh8ICzTHV1tfCWvO7Z0ODkzSDg21665ZZbcOqpp6KhoQENm32fo446Cu3t\n7RhlUtrR0ZFzM371KzEo8j8vuKWKmZYq4JOOQiDyAoufKoN+S0QccIBQLED+c8spZKQAVESi2vVJ\n4AJcyOJ7dbVZTQ7T1FydBd/Rka/UZNL1IlEd+HWyWZu4uaXGyUVlhRIpjY6K7CN+2pBb/5tAN7lJ\nwfH+Vx0tqNtRWshmO4LJoSwqpeWlhFta1M8FuBP8f/7jHcYA8mXm8MPFWQgqGS4tVRuDhRD88LA6\nPVtlOKxaZb8mJUm8reJCyxLPznkzCPgm+KVLl+Koo45yvLdkyRIcdNBBqGCzqaurC42Njejv70cs\nFlP+Axa53sttUcLU0uODSxsSZKJy2w0of+Y35llWZgve+vXAIYfY1rJsNRPBqxai3AiYxwTPPVf8\nP3++eRsTCTOCtyyzdQodwXd25hOr3J9+ct914ZwVK2xFwdvCd0OriIsm4g03AEcc4fzsiSeA0093\nb4+bB6gLCViW6BP+3M3N+cpEtVmtsrI4C37DBm+CV1nwXmtBdXW2nPjxyFatMiNeWWbIC1F5romE\n2pgsJONp6dL8w3T4/TlU58iS8bVo0aI8LqS0SOLNoFB0kt3LL7+MRYsW4Uc/+pHj/fb2dtTX16PL\npR6nKj+VQz4MmsPUgucWA3kEqo0SOrKQFwz9FqxavdpJns88A3z96+K1TAgk3KrysW5lTfn1SZgp\nrmyCeNxs/cHL66J26DyvTCbfmpL70091P74rlh8ADdht5R4Xv7aK4Is9Jc2NzNwMg/p6pzyr4u0q\ngk8m/afpcnR3eyuI8nLgf/8XuPBCWwl5hUe5F+mH4Ht7zer58NPWADsTT0W0o6Nqo8HvPOb44hfF\noeOEV17J/87zz9uv6f669F3ATosk3gwKvgl+p512ws9+9jMsX74ct99+O0488UT89Kc/xYIFCxzf\n6+zsRFNTE9pdJMgr1JRM6q0FU3eZTywiHtmaevNNPXmdeqrzb79lmtetyw87UPlZ2V11IwG39DIe\n76Ot17xvjz9enKeqQyJhRvBjY2Jy65QhlW7QWUfZrDfB+1k0dFuQo/Hk9+PPqAo3FUvwbmPktnhc\nU+OUcxXJqfqltLTwNOOKCjEevGy2CvfeC1x9NfDLXwJ/+IN4z2vu8UV/v4vAXjt4t98+f5xIiaiU\nydBQ/vwbHS1uU9p3vuNUyKpCdU88Yb+msS8t1UcLiNSJN4OCb4K/6aab0NPTg/322w+LFy/GI488\ngnMpLrAZlmVhYGAAqVQKe+45H4CV9++ooyx87GOLXO8VRDkGPrHczjy1LBEvl7HHHsDZZ9t/+81S\nGRrKF3KdNV4owXNhozxnHsu0LDE5dWQQj5t5RNmscL1140JWpiqGTL+XiU5+ZtPDRAD39RCVwubj\noJKBQjIrONxq/rhZvXJ/NjXlt18Vay8t9efxEGpq7EXbv/zF/buPPGK/pu/y+LIKtbW2wvJbbkJH\nvGRodXXlG1lu/Z7N5vddX1/hp3tNmybmFp/DqlO1uDdBY1lWJryzq65aBJkPU6kXHbwZFHwT/I47\n7ogXX3wRr776Kl577TXsJqeiABjcrEqrqqq0lst++3lbAqp4l1/wSU2dTqRy4on2Zy0tomCVjL4+\nUQaXlh1Ms3c45J1v5FZzMnvhBeemChmmri5Z8Jxk6X66BdiysnxPSaVsqL06RUQKRLdHY/Vq70VV\nP4TgphCprVzGuDOpCv0VcogOJwo3L8gtJMAV1YEHipPDZIJXhWIqKgoL0SQSosY9YFvwJmtLZNx4\nnUBVXl74GQu6+UWh0sHB/LZ6FdOT+6inp/DTt+bMEc/G78kLFxJU411Wpp+DqZSTN4NCQTH4WCyG\nPffcc/NCaT4GNq8ekia66CLn54mESGvikzuILfcq8EUkEkx55+puuwFvvKG2PDs6hNVKFq5MAkcf\n7d0G2XIeGhKCyqtRLljgjNvJMN3CT7VGuJtPAqjr4/LyfMtZdToRKQGdlUUEqiOLDz7wznv3k/7r\nRkqkELmMcYJXeWKF1HeiuueAe3oqeQeq+7a2ApdeKnbSPv64+lhEFamWlzuLh5mivz/fmDApbkay\nryI0GYVuytMZfVQue2QkX3l4hRflKPHAQH547+KLzdq34475yRG8vg8ZN9wgIxlMJvWLyJWV+bwZ\nBEKhVdJElFVz661Oq+7QQ/Pdaz+HV/sBd2Ffe038T4JKxLvDDs5qjRwUGqFBlUlARwp8oUyV+rdy\npb+j2viEMdnFzAXYK5UukcifkKoJ70Xw5LbqMm2GhvI/kyenn0W5ykpxNJ/KGiOFyENP/Nq//33+\nbwrJrNh7b/s130Qko7JShFRUslBeLjJ3fvxjW57k76lItbS0sCyadDq/to1JyOLPfxahnTfecP+e\nZRVuwetkp6bGlkk5xOkV1pP7qL0938Leay/1vWWZ2GmnfBnlMqwicM4dOu4W+2ucvBkEQiF40kTc\n1aBITlkZ8O1vi9ec4EzyXwsBvy5ZiyTMFM6gZroRPD9RnUNncX7+8yKP+v771ZPHa3FLBldUn/2s\n81xKFfgRgF7WVCyWTygqJUJk6VaVkaBztWULTbZ6/cTgk0ngvPPU1i09s5f7zlHIvJo3z37tdnjE\n6KhegaiMBJlEdCmAfkI03PucOdMp7yYhi/fec+4NkNtC6OoqfM+Abvyrq21ylI0Cr/UjWSG0twvl\nyPtj223VfSBb+nPmuHsMqrEkWUwm9etGVVVq3iwWoRA8Fczh6T7XXy+OPdu0Cdh/f/EeD1EUU4nP\nDaWlYhLceCNw0kniPVrJpqO03DwiIlYicnnBdOZM/aDtvTdw8slqq9TroAEZXKgGB93TyVIp0Sa3\nNNRDDrFfqwheRXbUBl24iLuquvbJlh0d0lEI3LJeqI1+sjgKmVc77GC/luuEc4yM6D2vQnLN6Tt+\nFBiX83jcaW0Wm3rNyXLDBn+KmkP33JWVttUtzyeaozoFLXuw9H1O3tOmqftAPm5vp53cN+Pp1loa\nGkQYVicDvNDYuKZJmqBns9TVMNb+6EeBb37TSeR8QhTq0nnh6quFFf+lL9n3luPL9L7KMiVhIjKR\nCaO1VaRJHXyw831uwahIRt7a7AVuwQ0NuStEElQe9pKfbWDAtrJVLrXbFnWddcbTNXUxXdna2m47\nZ5jDD9wymug+fqpTFhL65ATPn18GlSRWQbW4a5K2SgaSrjSEDDm0xpUw55SbbjK7nu7aq1eb71OR\noZOFqiqbkGWvhebGmWc636cwqbyZi7xyPt78LGA3bLutu3eiMjquvFKcArXjjvp5W1Oj5s1iESrB\nexWt5xa82w5Jbm2a4vLLxdZnlcDIm5dIq6oGjrQ9EaA88erqxClVS5Y43+eEqrJK/v53bdOV4BNm\neFifikhtApxpn48+6vzO8LCz9opM8CoZI0Wls5S4K6yL08sTf3i48E0nnOBlkiQry099+ULKPPAq\nlu++q/9eJuOexSTDpN1EXKYZITL58DAafy0fu2iSGsyfYe1a/4usZ58tPHx5ExOhosJ+TlmR0hyV\n26k71IMseu5lmJ6pW1mZzxM33CA4AHBugALEWs/VV9vX1s2dujpz3vSDCUPwOoE480zg2GP9t2H6\ndP0AyxYPWS+qLBNZeORYKFmq8iTlBK+arCaZCISzznIqpZER900+RK6XXWa/p4rVklK98057dy1B\n5SVSX5gcoqITZNmiHhws3J3npCQriUIODqmqKi4X3m0xO53WW28qAjXxPChJwJQPZCXIla2qZhPB\n78He7e3+C95df73w8HWpqqmU3S55zUUOoxK48uVYskT0Hfese3vN+jGRyCf4+fPFZsmXXgJ+9jPn\nZ3KbdGFAbsFPeIL3u1jgdiLK3Lnu2/Q5+ERxy1eXJxr9rSIFmmi6EA39VvZAuBD43W0oC+Z3v+tc\nZO3oEGGmAw9UpzOSfNTUCMv87bfFYiTH4KCtVPkBIgRViIWeXWfR8WwFXXxc3k7e3+9e7dEN/Psy\nMfuJTRMKOZWpp8e5n0KHgQE9gahk1YTgSaGaekCy0uaKmntw8vd0RMnBx3vjRv8ET+FaHcFXVtqh\nR9lYIUXF+6GkRJ+48fzz4pm4vPb1qcde2qAPIH+uDw+L9/beW7SBzx15Hrhl0UyaRVZK96k0lLxM\nRk8Ic+d6lzQl8MUPtwkuEwm5T6pJRcKjIzXdWHCPxKugE8dppwFnnOF8z7KciqWjQ5DCM88A3/9+\n/jV4eKSkRLiNslB67YBUuasUntIp495ee81C119yKmFfn60M/RI8fybd9nW/8BuH7+8Hfv5zYRVe\ncon+e4OD+lRgFcGbGAUk46bhEHl3M7+vW50evwTf3p5/3KEM2YCgpAOdYZZK2QTPPX8OPhcbGtw3\nrsn8MDiolnlV2FaWU3ku8XGm8g4EXRioqso/b5ogFILv7e1FMpnM1YbXgTTm0JB6kTWVAo45prDt\n2G6LVLGYc5Ao3KGy4Gny6EhNF2vm9zf1QABRNEu2fkpKnATPFZGKFE2s0FhMeAAEOfapmhx0xJtb\n/j6NlW6iylkJhS7GAU6Cl0mkkA1AgNmGH0CkqV58sQhfNDcDhx3mbnX39OivrQoLmcgMkZRpOIrL\nSmenU3a5fMlVQ+WYvAp8qvf0OK+nCifuuqvzbwq76OSmosImTp3xxmWgvR247jr3NnP09qo9IVUo\nVTZe5PnKvepbbnF+pvPiamrMedMPQiN4k5VgbqmqNOXatWLQvvlNkW982GHmbfCTh0sTT6UUKHyg\nI3jdRiceVzadgH/4A/CJT+RvzCgpcV6jt9cm2bPOAn7xC+f3TTKSslngqafsv485xvm5yuoaHhYL\n17pJWFpqk5VOCfAj4wBB8EQ8fmvB83vI7S30QGnTEM3JJ4vJy4lQlY9PZZo2btQvjLvlTruBrHzT\nfuNezbp1znHktY6Gh51EZGLBc3Ls7nYqblWfynF92kOgS7aorPQ+7clPZEP2pvr61DKvWruSuUCe\nr24JELyfeBuI4IPMoAFCIvh0Oo1yg8Ie118PHHmkWAyUteIZZ9hCtuuuYhHDpCwAwctF5wqAYo5u\nBK97HJqcMolza840HkwVH+WQjlztMZu17xeLCauf/8YkxW5oSExwsuIeesj5uW6yvPWW+sSdDRvE\nDke6nm6iqnKSC02R5QQvewaFHhximnKomviq0AEtjrtZ8LLsm+4EpTE33RHNZfLDD50yzQl+aMjZ\nJt0ZuxzcG+jpcSoo1XPPmeP82ys0Vl7unavvh+Blr2RgwHyBXeYC2YJ3exZ+D64IRD0oM970g1AI\nfmRkxKihyaTY/rx4sdnJL34Wbty+K1v3FBdTETG9xwWeb/OmQZInPH8e0xATkZJM0KWl+e6/nHnC\nJ5GJe08TkJ79/fedn+u8jjvvFOQwfbrTMmttddZlMY2n9/fbLm1rqzjN3hScRPg6zf77m8XSVcTo\nZn1xqBQIf4/khaZBOq2/tuwR+V0/MFWQXAmuXeskTC776bRTQZtUr+WWaU+P8xlUBM9r+W+zjZ0E\noFPMJmmMfkLXxaTrytQmz1c3RcTvwZ+nvNycN/0gtCwavwsFsnurWmTyE8t2s+A4OVKNEEC91Zys\nl6OPFq+vu06UHyCQuyZbbzxEoCP4H/7Q+TdNMpm8y8ry33NTYCYZGKTkdDFBXZt/8xvx/ymnuO9d\nMCWddFrUGn/8ceEZyKEiN+ji93/9q780Qw5TD1lFwtdfb78m8bcs20PUWZiyBe93zcl0kZVvLFy9\n2qnEuVeWyTgtcpM6UXKIhnuUKhnjKcy33WZ7Tm7PMn26XuZGR/0RvOzxpdPmxea8QjRuWY68jSed\nZBs3lEUT5AIrEBLBDw8PI+mzNJ+8eUHlAusEX5Wn6xb35sTAN4nw9LDRUXFSCxHaSSeJSfHVr6qv\nLZMFKaNsVk9E8+Y5F2SI4FVlA2SCd8vMMSE3mij8+X/yEztl0q3+fDwOnHOO0wKWrWFTgs9mxY7K\nT31KXNdP8SyVEj/wQKGwTaxg1c5T06J3qj4+8khBojfdZJ8hMDZm97HOq5EJnl971qz8FFcZhdR9\n+fBD5xhz4ymTcba1uhqQDmzLAw899PY6F7lVSpPLPfe23Ayz0lK9khwc9BeikQ2k0VHzSqYVFc7+\nkQ1PN47m8+Tkk8Ui7o03imsWwpteCIXgM5mM75VgkxKgOuLac8/899zi0JxwuXL43OeEhb5xoyCb\nvfZSu5cqgpfvRwJO71dW5k/kWMwpHDTh+ISNxch9c/7WbRHRhFxpeLglFYvZ7qWKaBMJcQDKU0+J\nPufWltwn/LOf/1wsHqsgKwY/Xppqg9TMme7KiUMVMzfdFaq7x7bbiiPdyFjIZm1rWCeTbrJfXe0d\nYy+E4Feu1CtBVdryl77kfj0/h3DI3+ec5iW7ujDNyEhxx/CNjZmf6BWLOe8ly75blIXLDfV/Q4O4\ndyG86YVQCH5sbAwlRRZ4V7lqukmlWuV3O5RAZ8E3NAgL3auyJdf+NPlk74I8ELJIEol8gh8ddSot\nGnBOXOXlQqDoM105WQ6TricB5c9aWmrHiVUE//rrouTygQeKPuSCLStf3r6KCrFZSwWZvPxsUOIE\nRVvEL7/cvP6MKiRnsl0d8N7bQHH1bNa2UHXyqzovgFBd7U16XgS/YIGoBbVhg6iLUlEh3nNTOLI1\nnEy6W7huvKT6jLeZy4pX2qzbaWHFEDwdfm4KHsKSx0/uJ/58XAaI4+rrxZwNgjdlhHTMBrSHgZhC\nNUl18TnVoobbtnHe4aZZExxcCGny6UJC/KBiWYDk35CS4O8TOVB/UHzPjeB1Xa8iAv78IyO24KrI\niN9T3nwmW4P8GXgc2qtNfkI0XEaWLBFF3/bbz7v+PUFV/dH0OEyv1Fcaa14mWKe85HABf666Ou97\neSmAc84RuzdbW0VdlJ4eEY5zO+9YNV5u6xNuBK/ygvkYmaxXEXTzdXS0sFpChJISf2fy8pCUPH6y\nouHhZu6hUv/zcFWxvCkjFIIvKSlB1mfumxxmUQm1LiyhGvSODr3gcwu5kPNtueVDg+Tlfqtyx2Uy\nIwLg16JJReTq5e4D+tou3NKg3/MUOH7eqkrB8uvKFqzcHv732Jg+PirLsx+C5xNru+3s7CbTHHgq\nEMWhq18kw0u8yU0fGbHJUhd+khUS77tp07yzx7wUgEyapaWi33VyQoery3Aj+E98Qk+Qqs2A/Jm5\nseCV4qsLoWWzhc1lQmmpPwueh5Vkr0MO0fBoAp83xB3Tp4s+L4Q3vRAKwcfjcd8NffRRYOFC+2+V\n8KkE+Re/sA8Q5rAsPVlwoivkoBH+e2qnjlQ4wctdIocIyGrmGp+OViMi5K6/DjrLjFuQ1G4+YSzL\nPZzglvqp2vpNKCnRk4NM8Kbxc0BfwsI0zCOXhADMs2i8iIjc9HTaljFd2FBe7OXk19DgHbbw6jNd\nW3WxfV24wi189elPiz0Sp52W/5kqv5zLD1d8XspZZ6WXlBR2XjIhHtf/XmUIcIKX55tM8Hw3rGoO\nNjaK+VwIb3ohFIIvLS3FqM+dJjNnOmuqZ7P57rssqBs3ik0+utibTvC5QJlsw5bBLSoiMtUkHB11\nErz8PHKIgJ7v+OMFKfz858CvfuW/fbqu56SbTouJzL0fHkpRWZvc+pBDMnJf8wmcSOjHSA45elmr\n3BrUEbzJzmFdiQnTEKhXnJ88lv5+u091ikeO4fK/6+u9lYkXwevaqiP4WEwdb/dSfjvsAPz618AD\nD7h/D3COMx8vrxCN7iCZYkPX5eX6ecNr/tMc5h5pe7tzbssErwtH0f3ImyqEN70QCsGXlZUhXcCp\nu/LmFJk0Oal87GO2ZaTbIajLEuDXlXfUmYBPGJpcXiGlWCzf6pbrXFC7fvtbYe1dcoltVfDDOYB8\nZcHvrxN2uY3ZbH54iywkFdHy5/ay4HkfJxL6HGPZavKKn3M3XBfyMBE9v4uCBCI5L4Ina7e72/la\nBTncxfu2rs47s0jFCXzDmN8SyCUl6t+YLkDLVrjKMubKmc9Tr/RW3SaiYqx3QDyv7t5cqdBY8L4Y\nG3MqYVl+eCSBe+c8SpFOF86bbgiF4MvLyzFcQEFuLlQqF5EPAFcGfi14vqmjkEVWPgGJyFSWB2/v\n2Fj+RJTJjJcfkEMXXqEZPmF05CPLjpz7u3KlLcyqFEJO2l6HZcsCb5oDXkiVSxkmoueWyubm1ZHl\n70VEXFHq1jWo7+VdxKtW2a/33NO7T1RhxhNOEGR0990is0gFXcw5FlOHQkxPkuPzrqRErdx5njzv\nF91RlmTQ6NpcTHZhKiXOnNDxBbfWySiT5ZDLnKxsqAQJ4JxXfC6PjBTOm24IheCTyWRBDeWa8qGH\n3Gt8807XDa7OYqKzSt3q0LtBFaJRrRl0d9tErVLM8mYuty4jwSahcMs+0QmqiuA5PvxQCG4ioSYV\nfl230gmW5SRAtz6WJ4NXuMEkRm4Sx3cjeFWNfQIpKq9QKSd1XRbN3XeLdjz4oPN9OqD+/PNF2QWv\ne/3jHyK9l6OpSYzl2Wfrifmb31S/X1LiNJrIazRdr+JjX1OjDl9yeeXfv/JK9TW96KSQeUx48EGx\nf0FXolmudQ/kG5XcIOPGzA03ALvskv97IP9QoEJ50w0TKkRz6qnAY48Jy+Ooo5yfjY05SYTHUBMJ\nUVVRhm5x6sgjgT/+UZ0HbQJOrkSEKoLnRKPKm3azemXIBC+DX7+3V61QdLnq1JeHHioEVRfyUuXs\nE/j95Da6kals3XmFPkxS4UzCmG4bBnXxecCc4MkA4UWs5PE99VQhPyec4Hz/8suFFf/LX4q/vbyF\nhgaxQY8rCpOMkIsvFhbzF7/ofF8up03zyNSC58+pWyTm8mpSMpr6QJfz7yZjXkYByYJuEZz3BbWb\nG5jHHefsG7eNaTzcxpXCpArRVFRUYMikpKGERAI44gi1Gy5PWjkWd9dd4t9114lFkXjcttRllJQI\nl8w0piiDCxkRkkpIR0bcY4NeefGq7+rIS47TqixwOSREJPXcc+KosUsvFX/rJjKXPfl5+d9yG93I\nVP7My4CRx0w14U2OAHSLS7tteSdL0Yt06frDw3abVbF0FRHHYrbs9vZ6h2iovws5bnDWLFFagSvO\nRMKpaInU+JyTD7jm48BlIZVSjxG3lk0InkIbKv7zymHnawIq5U3v6bLuuMxRO847T3hHS5aIMt+c\nrPnzyoYAvwdXBJlM4bzphglF8G6QY9LyjraSElFq4KtfBf79b6GNedXHIMEHhia6injp7NS5c4XS\nkeu2/9d/ifRQmsxu2+RJgKlb5Zi2CcHLoQtq8667iho01Me6Mzj5pHcL98hC7VbZUf7Mizhla0xF\nmiaE4dYmN4InQvZSRHT9/n733cFeUJVmlkF9X4zxxxVtVZUzAYDazUNX99zjjO3rFkpVh1QD/i14\n6gcVrdBubx04D7z9dv6eGyJ43WI2bx+1e9ddRYhNdUYFl3+3Uhy8zZOK4GtqapBOpzFSzHE9Ergl\nPG+efZCC7rvFbHrwAiczt41OAwOCEN5+WxxVd/rpzs/feEOEolauFGUArrhCf0/50G+Z4OXMBbfj\nBwk6d1cXa+VrD7KVzPtEVna6kMenPpWfN+1F8HKs1e0cXTe4Wbsmnp1XiIZitDwG7yfHn+BWcoNA\nJFJM+JYrzmTSmcJLpMQ9Ysty7j/h8sf7JpVS9xVPCtDRBCfAd94R/6vG1stz+elPhTJ6+WWxkW3x\nYufnNFa6owC5MjI5PpTPDbfj/fhnQ0Ph8GYoBF+/2cfv8lM5ygMlJaLq2jXXiIHS5cNuCfD+dyN4\nmng0kPK4Ud5/SYlYWHML5/BFnd13Bw44wPm5PLlVZCKvA+hISheD58IpW7myu8khL0jNmyfOAViy\nJF+ZmIYjCCqFUGyIJgiCp/4ZHHTuavULHelwFHpoOQd/5tJSZ5YLKXY+Vn19To+Tjxvv/4oKd+sa\n0Pfl8uX263/9S/yvIniv2v9NTYLk58/Pbx9gjxX3lp5/3n7Nn82kDAaXUTl0xI0kHp4bGAiJORZd\nyQAAIABJREFUNwO7EgMdO9Vb6LlpGnzpS8C3vmVetzksqNIFVQTPzzwF8oVTDtm44fzzhcV7//3A\nq6/mh3NOPNFZdE1FJnIWh24xSJc6ysn0S19yehzcAlSVO+a4+Wax0C1PfMvyJkE58yiMEI3bohxd\n24u04nFB7NTHsZjZSVsyTA65ccvkMgVX6qmU8/Q0Mha4AfL++87fcBLkoSIqi+AGVXjzuONEgTTC\n66+L/1UKna8fqDx32cDg1+DlOTjB87nC561JYgb/Pvc4h4achhgn/3Q6HN70UX3BHLWbe7zHT2nA\nSQQ+SGQVk1BTSYJbbsm3TjnxzJjhL4z0mc+IfzrMmwesWCEWmd94A9h77/zvyNaxbuLpYvBybPWH\nPxRrHkuXiq3qBNkikxWyjkBNdmnLFpQqRm1CdG5xdjcDgu9V8EJ9vQixdHUJsi8khGJy+AfJlenR\nfSrwNa1kErj9drEJ8P777RIinODfeAM45BB1O/k4JhLem5Co3fx3PLUQsC14lQXNCb6uLn+tQx5r\nrmjr6uz2UT82NDgVK//+v/+tfgYOHcHLaZjcghfF0oLnzVBDNN1+ztibROAWBwmTvBHDq0a9LgxS\nDJJJkUf8m9+oCUz2IHQTT7dbUEX8zc0izU93vmcikW/B68jRhJhlC151aIeJonDrfzcLnhS5CcGT\nl7VhQ+F52rr1BL55hmTPbXHPC7w/yspEuOfKKwWx0lGM/BnefNPp6fGQDm9HPG4eOuLzihZG6QAW\nUo6qg4C4N6YaOzkrjLdPNU9nz3YSfH+//awmSpob4HxeyN4Hn3/pdDi8GQrB07FTA6aFuScZuIDI\ni2BuO065Yg748HQjyKELXa60zoL/v/8zuw8f9oYGQYb83rqFQ5P89TffdP4th8EAswMw3M7NdMu1\np+cwuQdV6pQPuPYDOQT105+Kmvw//an9HrVFPpjdD3h/mBwtuHq1k7ife85+zZWfWxEv+fulpcJr\naGy0Q0SyJ6rKROJzSaV4ZQVDi8U77gjcd1/+92fMcC6m9vcLD9kUXP65Nygb5rxf4vFweDMUgq/a\nLCFbK8HzyS0PGgmYyoLi1ofXCfFhQB4OXe6w6ti6Bx4wDynxGCxZV/zZdXHMsjKzDTqzZglXedky\nccKUDC9ymzfPPlJPBbeNTuSdmFillFa4Zk3h60ayW3/00eKZufyQ0lFVKjQF97J0csEJSa6jxBdE\ned+UlHiPKclLSYkI/axebXsLvMIsoK42yUNFJvKz//5iI9mKFXZ5aC4z227rDP0NDZkf5Qg4x4HH\n/932h5SXh8OboRB89eYZ0u/39OBJAm4lyC4bkb/KkuDxw0Jq4BQLr7rVBJUF6+e0HBXBc6/z5ZfV\nvysvF4dvA2KCq0oGXH+9OCt3xx3FRFXBzbr+wQ/E/d223ass2L32cv69YIH+9wSy4NesKWwTEpDv\n1hOB8bAE1W/h1n7A6dR5eOMN598rVtivZeWnUhjnnmuHYXjtH0F09t/nn+/8HSfPWEysFVxyif2e\naU337bZzfpfLZ02Nc66m0/42RXIvgz/LwoX22baJhHPdIx4PhzeLIvjVq1fjoIMOwo9//GPH+xWb\npTnopP2JAk4gcl66biMS4JyAul22YcL07EiVS+3n3E+5Fgng9GhoC74K//Vf9mvVgtqXv+w92aiW\ny9FHi4wjjv/5H29rWhUv//SnhQfz0Y+KchqLFrlfA7Bj8G1t6r42iZOffLKzvW5yxcnPL0eYji+N\nZ3e3ID6ybPmYc/kZG1P35y67AC+8ILJjeOVLGTNmiHWAN94Q1+JEvH692FHKjSk/h3Zw8JBMebmz\nL0dH/VnwfA2KE3wsBvz3fwvvc+lS59pARUU4vFkUwXd1dWHp0qX4xje+gT//+c+591ObzYug0yQn\nCvhkIG0t7yJUWS184OX6I1sCfHLsvbe/vGk/Lj9XJGRl77GHvcHkf/9X/1sivXRavTPTJAXyyitF\naOORR5yLaKZKVSaJffYBvvAFEWZ6/nlRTsNkkZW8hLY2NfGYxMlPPVVs8rnwQrFRTuX5kQLQ1Vg3\ngWk1Rrke/6GH5n+HP6tlqZXb8LDow91281a4u+8udo52dzvnnkp+uTXvB9z4qqnJzwrSrUupwPlZ\n5Q3uuKN9kA8hmQyHN4si+L322gv7778/LMvCYrY9rLS0FNXV1VttFg2f3CQIfKHn9NPVseHTThM7\nV595xs5M2JI4/HBh/f7gByJ+7Qd+ziHYdVfRH2eeCXz/++K9WAy47DIx4b/3Pf1vZe+Bk3JrqxkR\nxWLqejqmYTF+jxNOAP75T2FJ+t1IRKGuf/9b7RWZbnzadlvgttvEgiC1gVv/9JobEH5j8KZZN/xQ\nnrY2e+Gdh8s4wY+N6QneL2SPTqVkTz0V+Mtf/F+bW+zNzU6CHxszP8oRcD6b6TmxlhUObxadB3/S\nSSfh2Wefzdte29zcjA0m+6wnIbgA02BeeaVw3T//eTtEIGP33UXtmfHCLbeIU6IK2fHoJ0Sz445i\nf0Ah7rJM4I2N4t977wF/+lNxBzvoxkUGv8fXvla420+L0qtWqTOHhobc8/HdwK1EWsDmJOU3i8Z0\nXe+XvxRktHy5iJ3X1wsy5DLFQzKWpV7vKITg5QVnnSwcfjhw7bX2Rjxx3qn7teV9AFxZ6p5BBz4O\npunQtKYTNG8aT/UVK1YoE/BP2BxreO655xwk39TUhI6ODoyMjCAWiyn/LTIJZE5AcBKiY/mOPlrE\nEk1JZDwgl4HVQWXN+c2rLpQUYzHnxE2lgGefFXnWqrN3/cBkYRRw9lGhi6OAM9NFRWhB7Uh/9lnx\nP7dw/ShkwDzLp6JCHMv3/vu2l1RV5ewnTvBjYyJM9sQTwLvv2u+b1HSRIfOemzfHw3sm+4b22Ud4\n3qWloqQBt+BjMX8xeF7HR6WEFi1alMeFFH8n3gwKRgS/cuVKzJ8/H42NjUgmk0gmkzjppJMA2Cu/\n6XTacZ5gKpVCf38/ysvLUVrMcSsTEPKi0da2lqwKHRSzS9IvuEufSgnycUtddAMnVtOwWFAE72Wl\nFrNhsbJSrAUAtuIqJk2ShxKKHWuuLMbGBEF+8pPO4zFXrvR/XTk07Was8DmpO8hDxr33iu/OmeM0\nBrq6nJa417gSwZus09jXHMbo6GiON4OCkZ01e/ZsLF26FBkmNXuqtoAxtLa2Yvnm5Njp06fjQ/mE\n6UkMVQ1zr4JHkwkq68qPsBaLZNLOyig0hEHgrrZpPJQ/azF1j7zKAxdbNPBPfxLxZio8xzNZ/IZA\nuMU7MFCcPMsErwI/ltAUcgzezW7kfbFunfPgbB14XZo//hG4807g6qvFugOXia4uM4t+9mzv73B0\ndHQ4eDMIGBF8LBbDQnnHgQdaW1vRtnlEVhd6dNIEheoUojDLE29pqDzEQkMuhYB7SH7y71XgE72Q\nhe1i7u9F8MUUBwOE+081gCzLvV6/F3jVyqGh4gieK2UdwZvUuZchlylwW4/h1r6qvIEXEgmRuXTh\nhfmfuRE89xx23139nUWLFinD06+84uTNIPD/2zv36Kiq649/ZybzyjwyycwkEAK0YKCIAoryswYW\nxRSoCDWEyENh1a6qQEURqxhAUGFpbbVUrMrC+uRRaJeAol22ClQIApZIS4t0STBIEPKYTDKvzHuy\nf38M9+TOJCaZzL0zSbyftWYlmTu5954z537PPufss7dgG53kcjkUvBo3mUzweDyIJLra0weIf+i7\nOwTsK3RUnmRyXiYKf4omWQt+yJBoLPDXXuu+Nc4Xpe5a/R3R1VSqkNnZfL7YhdVEOw9+55Ds2kBn\nuze5qY6eTAPFdzqdTdHwyyO0M19nXoz8Tj2RhVkgOsUmtG4mbZdlZWWhoqIC2dnZ0PCeIP6urKxk\nnpJeSLzodJWkoq/R0dxwKkM08wU+WQteJouN29Id+NM6PY0hA3QdO1zAsN/t1oEScWsFYkWwoSE5\nN16+y3B8JzZoUM8NojVrolb7E090vaeB/0wKHdS2s8Qt/LIlGlDQ7RZeN5O24GUyGSZOnIjRo0fH\nvJ9zpXS2niyX93LiF/ySHWr3NjpymUulBc/vTNKxtiGUE0NXQiakZRkvpInOwfPX9ZI1WPgGUPx6\nYUfhJ7qLSgWsWxedcqms7PyzfCtbQKcUAJ27lPLLm2jO5+hirrC6KUosGgAYcCXlUn1PJsB6OfGR\nIBO1lno7PcmaIyR8z5V0rG3EB9LqKR0JPN+NVkiBj2+DiZ6bL+o9SS3Ih99W4tuSEFFUc3O7bhf8\nMggd87CzxXF+PSYajuTyZeF1UzSBt1z5Bhp7kmm4lxM/9PouWPCJWiPJwJ+WScT/WCg6ijHfEzra\nr8Ivj5ACHz9lm6hI863cZL30+CPc+KkjTpiT2bDWHfh1K/QaWWfPO79zmzUrsfPa7cLrpmgCz80f\n9cd4NPxdb3l5375a3lfpaM4ylcsofIFPx/INF50xWTqaGuCXR0jLMt5bJdH5fb4gJjtFwy9j/LnW\nrImGrPjqq+Su0RX8jkXofSqdCfz//V80Fee5c4mPeh0O4XVTNOc3Lr+gO9nxXi/k5pujQ+0bbgCe\nfz49sd3FpKO2lco5eP4cbjry7wolPh1ZjvyRkJCLrPFeKYmKGl+IkxV4vrDFt6WBA9uCzokJv26F\nXmTtbEpWLgfuuadn53W7hddN0QSeWw3ujxZ8bm5bEuD+SEdD9GS8SRKFL4LJukn2BKH8AjoSFr4x\nMGOGMNfpiERFmj8jkOwGLH6gN5cruuCbyo7a748tg5DuqEDiYSC6i8slvG6KNkWj1WqhUCj6pcD3\ndzr6ylK50ak76ePEpCfb6DuiI4EvK4vGvH/99Wh0UbFIRIQikVhBTFbgtdrOk02LTXznJrQbsxhO\nFa2tUQteaN0U7bGVyWSCx1WQSA3xo8NVq7ofL1wI+HO46XCT5AQpmWsTdbyIes01wNGjPT9vd0lE\n4OMFXQhtMRrbRgWpztwZL8BCS5AYThV2e/S5E1o3RbPgAUCtVrcLIyzR+4m3eNatS+31+cLa0yBj\nQtBFuKVO8fnai6xcnrqprkSszHifeSG0he+I0JPQwMkQPyUj9JSKGJvzL19u62iF1E1RBT4zMxPe\n/rbN8ztA/FeWaj9//hRNOhZZObicoT2hI6tVoxEvaFv8eRN57OLvVQhtueLOnfC9CEH8/fcFgf/m\nm7Z6ElI3RRV4nU4naIZwidQQb8GJtaj0bfAXWVM59w/EPrxc0uyekOrNYvH1lMjCYrw3jxCLkvy9\nIqlehosXYCHaL99LSYw26XS2tRkhdVNUgVepVAgKvYQtITrxgUNTLfDptNr5C4LJbO7qyE1xzZqe\nn68r4tdIEpknjl8MFiL2f7wnTSrh7p8TYiEsbn599iQjWle4XG0dq5C6KarAK5XKmBjyEn2DN94A\n9u6NJuaeMEGY7eWJ0FE+1VTBF+ZkfP+HDQOmTIl978EHe36+roi/10SmWeJHbEIIfGfhClJFR/lr\newpfxsRIfmO3t11DSN0UVeAzMjJisjxJ9A1kMqCkBPj882juTTEsls7gb+dP9eiBn5g7Gc8htRo4\ncAB4//1oNqNXX03+3jojftSTiMCLYYMNG9b2e6oznsWvRwjRhvj1KYbAO51ta11C6qaoM5xyuRyt\nqX5CJfo8/DlOmy0aDiJVKBTAkSPAn/8MzJ6d3Llksmiu3pkzhbm3zogPq5xILPL4R1QIAeOPXlK9\nyMoJPDfLIcTCdnwSbqFpbm77HoTUTVFtM7lcDkplMk+JfgG/yQiY3KbbFBUBL76Y2t27ySKTteU8\nvf/+xEYMYljw/Ojh/GxRqSB+xCnE4jZ/nUIMSXO52s4rpG6KKvCyVCbylOiXpFoc+jKffhoNdfzS\nS4ktEMfPBggxJSeTAY8+Gu0kO0p7JyZZWW1W+y23ADt2JH9O/ihEDFlramoTeCF1U9QpGiKSRF4i\nYfjGSz+MNi0aiaaI44g3FvmblJLh2WeBDRtSPxKyWICPP45a7gmmkv5W+Ba8GJL2zTdt5xVSN0UV\n+EgkAlUqwxBK9Av4Xh3pmKL5rsFPnrFmDfDYY8KcN5U7d+MpLhb2fPzQz2LEsq+raxs5Cambogu8\nQuzI/hL9Dr4vujRFIz5Tp0bj44wcmXge0e8KfFdPMeIyORxtrq5C6qaoc/DhcBgZqd6KKNHn4e+s\nFDJmukTHyGTRCJeSuH87/Dl4MSYl1qxp8x4TUjdFVd9gMChN0UgkDF/UpTl4id4Af4pGjGmn9evb\ngtsJqZuiWvChUAjKVMaZlegX8Le2X7qUvvuQkODgJ4HhJ4UXCrkcuOOO6O9C6qaoAh8IBKDuS87E\nEr0CvsB3lLhaQiLV8L1oxE5CI6Ruiirwfr8fmnRGjpLok/C9aGpq0ncfEhIc/Dl4MSx4PkLqpqgC\n7/V6kRm/h1pCogvuuAMYPLjtb2kztES64SctEVvghdRNURdZPR4P9OnIuSbRpxk5Ejh9GnjqqeiO\nTGmvnES6SaXAC6mbogl8MBhEOByWBF6iRxiNwO9+l+67kJCIwp82FHNSQmjdFG2KpvmKr5spmawJ\nEhISEr0A/sK/mIusQuumaALfcGWPuVWowBYSEhISaYK/k1XMSQmhdVM0gXde8SvKTmd6HgkJCQkB\n4FvwYgq80LopmsD7r6xKSH7wEhISfR1uDl6rFXfRX2jdFG2R1XWlyzMYDGJdos/w5JNPdvj7dxGp\nLqJI9dBGb68Lvg88P52kGAitmzISKeXSa6+9hnvvvRcXLlzAkCFDxLhEn4Ef2/m7nuFKqosoUj20\n0dvr4vLltly9118fzVUsFkLrpmhTNPX19QCkRVYJCYm+DT+l4d13i3stoXUzKYE/fvw4xowZgyVL\nlrQ7ZrfbodPpoBV7V4CEhISEiAwZAixZArzwArBsmbjXElo3eyzwzz33HCZOnIja2lps2bIFFy9e\njDne0tICndhReSQkJCRERiYDNm8Gli8Xf1e10LrZo0XWv/3tbygvL8fvf/97FBcX45prrkEwGIz5\njMvlgtFoFOQm+zqRSARyITIZS0j0Q5YuXQqPxwMPf7vod5Rly5bh4MGDgp2vRwJ/ww034K233sKi\nRYuwe/dumM1m5OXlxXymP1jwoVAIHo8HPp8PbrcbLS0t8Hq9aG5uhtPphMvlQnNzM1wuF3w+H3w+\nH4LBIPx+PwKBAILBIEKhECKRCFpbWwEAcrkcSqUSGo0GOp0OarUaKpUKRqMRRqMRmZmZ0Ol0yMrK\nYi+TyQSdTseOZWZm9tlk5l9//TUMBkO7tkFErD5bWlrg8XhYnbe0tMDn88HlcsHhcLBjHo+Hbe2O\nRCKIRCIxv3NwdSWTyaBQKKBUKpGRkYGMjAyoVCpkZmZCq9VCr9dDr9fDaDQiKysLBoMBFosFJpMJ\nJpMJVqsVWq1WkLrPz89HOBxuZxilAiJCU1MT7HY7XC4XXC4XnE4nHA4HGhsb4XA44Ha74fV6WXvm\n2nY4HGZtmVsQlcvlUCgUkMvlrE6VSmXMS6vVwmAwIDs7G2azmbXj7OxsPPDAA9Dr9dDpdP0+h0Qo\nFILNZkNtbS2am5vR2NgIm83GtMXtdqfOgvf7/Vi7di1Gjx7NYiOMHTsWhYWFWLRoEWw2G37xi1/g\n5Zdfbhc7we/3s4hoCxcuhE6ng06nY51BTk4OjEYjNBoNsrKykJubi+zsbMEzQBERE2iuITc2NqKu\nrg5NTU1wu91wOByor69HfX09PB4Pq3S3292ta6hUKmi1Wmi1WqhUKmg0Gmg0GqhUKmRkZEChUECh\nUICIEA6H0dLSAr/fD6/Xi0AggEAgAJfLxXxgu4NOp0NOTg4sFgt7aCwWC/Lz85GTkwOTycR+ajQa\nGAwG6PV6aDQa6PV6wXPlBoNBOBwO2O12OJ1OeL1etLS0oLGxEU6nEx6PBzabjQm0x+OB3W6H3W5n\nnWR3PSgyMzNZeeLrmP87J8REBCJCa2srIpEIQqEQwuEwE1iv1wu/3w+Px4NAINDptbVaLSwWCywW\nC6xWKxN+vV6P7Oxs1o6zs7NhNBpZhxHfti8lmcmEiFjbcbvdsNvtaGhowOXLl5lwc3Xc0NCAxsZG\n1nna7fZOy6lQKKDX65GZmcnaNr89y+XymBFpMBhkRgxXp1wdh0IhBINB+Hw+eDyemI7328jKyoLV\namVt22QywWg0QqfTwWg0wmq1wmw2sw4jJyeH6Yterxe1g4hEIqxT9Hq98Pl8aGlpgdPpZEZHXV0d\n6urq4HQ64XQ6YbPZ2HfQlaYYDAaMGTNGsPvt1E3yxIkTmDdvHuuxAeDBBx/Eww8/DADYuHEjnnvu\nOVy+fLmdVfPUU09h//79AIC6ujq4XC5mjXVGRkYGa1wajQZqtRpKpbLdgyyTydhDyzUsrjFxvWEg\nEOhWo9Lr9cjNzcWAAQOg1+thNpsxYMAA1rA4645rRCaTCVlZWezBFWpTQigUYsLINQ6+Vcs9JJyF\nyz20nEjW19ezWBZdwT2wKpUKarUaGo2GWVvcMa4T4ASSs4w5i45r3Fxdd4XBYEBOTg6zlLOzs5lI\ncqMV7mE1GAwwGAxs1MJZgFlZWaImcg+Hw8wQcLlcsNlscDqdrJ5tNhsTTbvdzixht9uNEN/dogM4\nsdTr9dBqtVAoFOw9pVIJuVzOniOuM+LaNff9B4NBBINBeDwehMPhTq/HWci5ubmwWCysbs1mMwYN\nGgSLxQKDwcBGLCaTCRaLBXq9XpQRIjdKa2pqYu2mqakJzc3N8Hg88Hq9zABraGhgxxwOB+uwvHyn\n9E7KzdUx1765kTKnHdyLryFcJxWJRBAMBhEIBOD3+9nL5/N1WeccZrOZaYTFYkFeXh77DjitMZvN\n7G+9Xg+1Wg25XI7S0lJMmDAB5eXlyVZ5z/3gv/rqKxQWFuLw4cOYOHFit//P6/WisbGRPRSBQADN\nzc3trDvOqgoEAky4ueE3v8ORyWRsaMh1BNzwT61WMwvKaDSyxsxVuNVqhU6n61eJwQOBABwOBxtu\nO51OZp263W72e0tLC5tG4hpyKBRCKBSC3+9HMBhk9cw97JyFrFar2YiFq2uj0QiTycSsV+4h4/7m\nrO3+DGdJc8LEiZLb7Y6ZyuPEmrN2fT4fQqFQu3bN1bdSqWSdHNcpc+2aa+OcUOTn5/fLds0RCoXY\nyM/tdqOpqYnpBmf4cH9zHSI30uFGGvGSx2kIN82kUChYh8CNxjUaDRvJcB0iZ4Tq9XpkZWUxw8Vq\ntXaasOOLL77AgQMH0NjYiKamJjQ0NMDtdrNpyJKSEqxbt06Q+uqRwEciEdx22234+9//jq1bt0Kn\n0yEvLw9FRUWC3JSEhISERPL0SOBdLhdycnJARMjPz0dTUxMikQjcbne/XiCRkJCQ6Ev0yHfPaDSi\nqqoK1dXVuHjxImpqavDll1/GiPuRI0ewfft2+Hw+wW62N3L8+HHceuut+PDDD9sdO336NF599VU4\nHI52x4gIBw8exK5du7o1f91XaGhowEsvvYTTp0+3O3bmzJlO6+Mf//gHdu7c2a/qozOOHTuGrVu3\ndjivHIlE8P777+O9997r1sJkX6a1tRUffPAB3n333Q7nuH0+H7Zv344jR46k4e7Eg3NCiKe6uhqb\nN29moYPjOX78ON5+++1urUeABKalpYVmzZpFAAgAjR8/ns6fPy/0ZXoNr732GgEgk8lE1dXVREQU\nDodpyZIlJJfLCQANGTKETp48yf6nubmZiouLWR1NmjSJLl68mK4iCEZFRQX94Ac/IAA0ZswY9n4k\nEqH777+f1UdBQQFVVlay4w6Hg6ZOncrqo6ioiGpqatJRhKTxer20YcMGWrVqFXtNnTqVSktLWZm8\nXi/NmTOHlXfs2LF07tw5do6LFy/Sddddx47PmTOHGhsb01WkHuN0OmnSpEnkcrmIiMhms1FRUREF\nAgH2mUuXLtH48eNZWWfPnk02m40dP3bsGA0aNIgdX7VqFfl8vpSXJVmcTic98cQTrE1MmzaNzGYz\nabVa2rx5MxERtba20po1a0ipVBIAMpvNdOjQIXYOn89HZWVlrC6uvfZaqqqq6vS6ggv82rVrKS8v\nj6ZNm0Y7duygu+66i0aPHk1NTU1CX6pX4PF4KCcnhwBQeXk5ERFt376d9Ho9zZw5k5599llau3Yt\nWa1W+uqrr4iIaPny5VRQUEDFxcW0d+9eKikpoRtuuIHcbnc6i5IUhw8fJq1WS9OmTaMVK1aQWq1m\nx3bu3Ek6nY5mzZpFzzzzDD355JNkNpuZqD388MM0aNAgKi4upj179lBpaSldf/31TBj6Ert27aKM\njAzKyMigyZMnM3EyGAz0zTffEBHR008/TVarlaZPn05vvfUW/fznP6cRI0YwES8pKaERI0bQpEmT\n6PDhw3TTTTfRzJkzKRgMprNoCfPFF18QANbuX3jhBQJAHo+HfaasrIwKCwtp0qRJ9Mknn1BRURHd\neuutFAwGqbW1lYYOHUoTJkygadOm0aeffkoFBQW0atWqdBWpx3z88cdMmDMyMmjChAk0e/Zsmj17\nNq1cuZKIiPbv308qlYpmzZpF5eXl9OKLL5JOp6N///vfRET061//miwWC/3kJz+hN954g+655x4q\nLCyM6RDjEVTgvV4vDRw4kLZt28beC4fDNH78eNqwYYOQl+pVLFy4kADQ8uXLiYjopptuoscffzzm\nM/PmzaOFCxeSw+Gg7Oxseu+999ixQCBAI0eOpE2bNqX0voXkzTffpHnz5lEkEqGtW7dSdnY2O3bz\nzTe3eyjvvPNOWrBgATmdTsrJyaE9e/awY8FgkEaNGkUbN25M2f2LgdvtpuLiYjIYDHT06FEiin7X\nQ4YMoT/+8Y/sc5FIhH74wx/SmjVr6Pz586RQKNhDTRS1/kwmE7377rspL0My2O32GIE/ffo0abVa\nOnjwIBER1dTUUEZGRsxozuVykdlspnfeeYf27dtHarWa/H4/O37q1CkCQBcuXEhtYZIg/JwYAAAH\nBElEQVRkw4YNJJfL6ciRI9/6mZKSErrnnnti3lu+fDn9+Mc/pmAwSEOHDqUtW7awY5FIhIqKijrt\n8ATdP19RUYHW1lbMmTOHvadQKGC1Wvt1XPiysjIA0fJXV1fj5MmTWLx4ccxnBgwYAIPBgI8++ghZ\nWVmYMWMGO6ZSqdjGjb7K3XffjV27dsHlcmHTpk1YvXo1gOjO1RMnTrQLSMfVx/79+6HX6zFr1ix2\nTKlUMh/tvkogEMD48eNx4MABfP/732ehX48dO4aWlhbMnz+ffVYulyMvLw8GgwHvvPMObr75Zowd\nO5Yd5zb59LX6+N///hfz9+jRo5Gfn88iJu7evRs33ngjxo8fzz7DbWAzGAzYtWsXFi9eHLPPZODA\ngcxVty/BuR1v27YNL7zwAsvcxGG32/HRRx9hWVw0M+45OX78ONxuNxYsWMCO8dvNtyGowFdVVWHo\n0KExkdAqKytx9OhR3H777UJeqlfBJch1Op2orq6GVqtFQUEBO15bW4vt27dj3rx5OHfuHK666qoY\nH+VDhw7hzJkzuO2221J+70Kzbds22O12PPTQQwCA8+fPQ61Wx8S2rq+vx9atW1l9DB8+PKY+Kioq\n8J///AczZ85M+f0LRTgcxujRo3H77bfDZrNh3Lhx2L9/P6qqqlBQUBCz8/vUqVM4ePAgSktLUVVV\nhZEjR8aca8eOHWhtbcXkyZNTXYykqKqqgtVqxUBelgwiQk1NDTseX9Zdu3bB7/djypQpHR7/zW9+\ng+Li4nahUXo7ZWVlUCqV2LJlC1asWIEJEybELBrX1tbC6/XGlNfj8eDll1/G/Pnzce7cOeTn58eI\n+X//+1/s378fpaWl33pdwXdCeL1eEBFkMhkuXryI+fPnY926dfje974n9KV6LdyGIaVSCbfbjbKy\nMpSWlmLy5Mk4evRozOr3uXPnsHDhQvz2t79Fbm5uGu86eaqqqrB+/Xo888wzMYLNrw+Px4OysjL8\n9Kc/xS233ILKysqY+qiursZdd92FZ599FgMGDEhHMZImHA5DpVJhz549AKI7uQcOHIhXXnkFM2bM\ngM/nY89IbW0t5s6di5UrV6KwsBAAYurjn//8Jx544AH86U9/EnUHrxhwhl5NTQ2am5tx6NAh1NTU\n4LHHHmMWOL+slZWVuP/++/H2228zjzz+8R07duD1119HZWVlCkshDGPGjMH58+dx6tQphMNhbNq0\nCdOnT0d1dXVMZ+X1eqHRaBAIBLBgwQJcd911mDt3Lt58882YdlNXV4e5c+fikUceadcJxpDczFIs\nZ8+eJaVSSfPnz6df/epXNHjwYPrDH/5Ara2tQl6m1/HJJ58QABo+fDi5XC6yWCw0depUWr16NQ0b\nNozKy8vZAtm//vUvUigU9LOf/YyWL19OQ4YModdffz3NJRCGBx98kIYOHRrzntvtptzcXCouLqbV\nq1fT8OHDaeXKlcyT4tSpU6RQKGjRokX00EMPtZuf7os8/vjjVFBQQMuWLaO1a9fSlClTaNiwYVRT\nU0Nff/01qdVqKisro0ceeYQGDx5MGzduZM/I+++/T3K5nJYtW0b33XcfDR8+PGa9pi/hcrlIrVaz\nxUXuNW7cOKqtraUPP/yQZDIZ/fKXv6QlS5bQsGHDaO/evez/169fTwaDgcrLy2nu3Ll0/fXX04kT\nJ9JYIuHYvn07AaAvv/ySiKJrM4WFhWwt5uqrr6bFixezBemamhrSaDQ0Z84cevTRR2nw4MH0/PPP\nd6mtgnvRHDlyhIqKiuhHP/oRHT9+XOjT90o4gR85ciQRRb0HSkpKaMyYMfTBBx+0+/yBAwfoxhtv\npOnTp8e4T/ZlPvvsM/YAW61WysvLozlz5hAR0ZkzZ6i0tJSuvfZa2rdvX7v/PXjwIPOU+Pzzz1N9\n64JTXV1NmzZtolGjRpHFYqGlS5cyDxqiqOvf5MmTaeLEifTpp5+2+/+//OUvNGrUKJo/f36M+2Rf\nZOfOnbRx40b661//SpcuXaL8/Hx69dVX2fHdu3fT1VdfTXPnzqWzZ8/G/G84HKbnn3+ehg4dSitW\nrOjTnnitra3sFQqF6M477ySVSkX19fXsM+fPn6dFixbRiBEjaMeOHe3O8dlnn9GUKVOoqKiIKioq\nunVd0XKyfpdwuVw4fPgwhg8fjlGjRqX7dtJCc3Mz1q1bh0GDBuGaa67ByZMnYbVasXTp0nTfmkQv\n4qqrrsK0adPwyiuvpPtWUsq4ceNw6tQpFszwwoULePrpp5kzglhIAi8hIZEy9u7di0mTJsFisaT7\nVlLK5cuXsW/fPlRUVODs2bO47777cO+994p+XUngJSQkJPopUh45CQkJiX6KJPASEhIS/RRJ4CUk\nJCT6KZLAS0hISPRTJIGXkJCQ6KdIAi8hISHRT/l/dx3YMbLYhvEAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x10f6fafd0>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "As you see, the error grows with larger $n$, and we have to find out why.  \n",
      "**Important point** is that it is not a problem of the algorithm: it is a problem of representing  \n",
      "the matrix in the memory. The error occurs in the moment when the matrix elements are evaluated approximately."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Linear systems and inverse matrix\n",
      "What was the problem in the previous example? Why the error grows so quickly?  \n",
      "And here is one of the main concepts of numerical linear algebra: the concept of **condition number** of a matrix.  \n",
      "But before that we have to define the **inverse**. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "The inverse of a matrix $A$ is defined as a matrix $X$ denoted by $A^{-1}$ such that  \n",
      "$$\n",
      "   AX = XA = I, \n",
      "$$\n",
      "where $I$ is the identity matrix (i.e., $I_{ij} = 0$ if $i \\ne j$ and $1$ otherwise).\n",
      "The computation of the inverse is linked to the solution of linear systems.  Indeed, $i$-th column of the product gives  \n",
      "$$\n",
      "A x_i = e_i,\n",
      "$$\n",
      "where $e_i$ is the $i$-th column of the identity matrix. Thus, we can apply Gaussian elimination to solve this system. Moreover, if there are no divisions by zero in this process (and the pivots do not depend on the right-hand side), the it is possible to solve the system.\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Inverse matrix and linear systems\n",
      "If we have computed $A^{-1}$, the solution of linear system  \n",
      "$$Ax = f$$\n",
      "is just $x = A^{-1} f$.  \n",
      "Indeed,  \n",
      "$$   \n",
      "    A(A^{-1} f) = (AA^{-1})f = I f = f.\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Neumann series \n",
      "To study, why there can be such big errors in a solution (see the example above on the Hilbert matrix)  we need an important auxilary result:  **Neumann series**:  \n",
      "If for a matrix $\\Vert F \\Vert < 1$ then the matrix $(I - F)$ is invertible and\n",
      "$$(I - F)^{-1} = I + F + F^2 + F^3 + \\ldots = \\sum_{k=0}^{\\infty} F^k.$$\n",
      "Note that it is a matrix version of the geometric progression. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "##Proof\n",
      "The proof is constructive. First of all, the show prove that the series $\\sum_{k=0}^{\\infty} F^k$ converges.  \n",
      "Like in the scalar case, we have  \n",
      "$$\n",
      "   (I - F) \\sum_{k=0}^N F^k = (I - F^{N+1}) \\rightarrow I\n",
      "$$\n",
      "We can also estimate the **norm of the inverse**:\n",
      "$$\n",
      "  \\Vert \\sum_{k=0}^N F^k \\Vert \\leq \\sum_{k=0}^N \\Vert F \\Vert^k \\Vert I \\Vert = \\frac{\\Vert I \\Vert}{I - \\Vert F \\Vert} \n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "## Small perturbation of the inverse\n",
      "Using this result, we can estimate, how the perturbation of the matrix influences the inverse matrix. We assume\n",
      "that the perturbatin $E$ is small in the sense that $\\Vert A^{-1} E \\Vert < 1$. Then\n",
      "$$(A + E)^{-1} = \\sum_{k=0}^{\\infty} (-A^{-1} E)^k A^{-1}$$\n",
      "and moreover, \n",
      "$$\n",
      "  \\frac{\\Vert (A + E)^{-1} - A^{-1} \\Vert}{\\Vert A^{-1} \\Vert} \\leq \\frac{\\Vert A^{-1} \\Vert \\Vert E \\Vert}{1 - \\Vert A^{-1} E \\Vert}.\n",
      "$$\n",
      "As you see, the norm of the inverse enters the estimate."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Condition number of a linear system\n",
      "Now consider the **perturbed** linear system:\n",
      "$$\n",
      "   (A + \\Delta A) \\widehat{x} = f + \\Delta f.\n",
      "$$\n",
      "Then the algebra begins!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "##Estimates\n",
      "$\\widehat{x} - x = (A + \\Delta A)^{-1} (f + \\Delta f) - A^{-1} f = ((A + \\Delta A)^{-1} - A^{-1})f + (A + \\Delta A)^{-1} \\Delta f = $  \n",
      "$ = \\Big[\\sum_{k=0}^{\\infty} (-A^{-1} \\Delta A)^k\\Big] A^{-1} f + \\Big[\\sum_{k=0}^{\\infty} (A^{-1} \\Delta A)^k \\Big] A^{-1} \\Delta f,$  \n",
      "therefore\n",
      "$\\frac{\\Vert \\widehat{x} - x \\Vert}{\\Vert x \\Vert} \\leq \\frac{\\Vert A \\Vert \\Vert A^{-1} \\Vert}{1 - \\Vert A^{-1}\\Delta A\\Vert}\\Big(\\frac{\\Vert\\Delta A\\Vert}{\\Vert A \\Vert} + \\frac{\\Vert \\Delta f \\Vert}{ \\Vert f \\Vert}\\Big)$.\n",
      "\n",
      "The crucial role is played by the **condition number** $\\mathrm{cond}(A) = \\Vert A \\Vert \\Vert A^{-1} \\Vert$.  \n",
      "The larger the condition number, the less number of digits we can recover.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "The spectral norm of the matrix is equal to **largest singular value**, and the singular values of the inverse matrix are equal to the inverses of the singular values.  Thus,  the condition number is equal to the ratio of the largest singular value and the smallest singular value.\n",
      "$$\n",
      "    \\mathrm{cond(A)} = \\Vert A \\Vert \\Vert A^{-1} \\Vert = \\frac{\\sigma_{\\max}}{\\sigma_{\\min}}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Hilbert matrix (again)\n",
      "We can also try to test how tight is the estimate, both with ones in the right-hand side, and with a random vector in the right-hand side. The results are strickingly different"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "\n",
      "n = 100\n",
      "a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]\n",
      "a = np.array(a)\n",
      "rhs = np.ones(n) #Right-hand side\n",
      "f = np.linalg.solve(a, rhs)\n",
      "\n",
      "#And check if everything is fine\n",
      "er = np.linalg.norm(a.dot(f) - rhs) / np.linalg.norm(rhs)\n",
      "cn = np.linalg.cond(a)\n",
      "print 'Error:', er, 'Condition number:', cn"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error: 2.81450551264e-08 Condition number: 3.33972881047e+19\n"
       ]
      }
     ],
     "prompt_number": 10
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "And with random right-hand side..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import matplotlib.pyplot as plt\n",
      "%matplotlib inline\n",
      "\n",
      "n = 500\n",
      "a = [[1.0/(i + j + 1) for i in xrange(n)] for j in xrange(n)]\n",
      "a = np.array(a)\n",
      "rhs = np.random.randn(n) #Right-hand side\n",
      "f = np.linalg.solve(a, rhs)\n",
      "\n",
      "#And check if everything is fine\n",
      "er = np.linalg.norm(a.dot(f) - rhs) / np.linalg.norm(rhs)\n",
      "cn = np.linalg.cond(a)\n",
      "print 'Error:', er, 'Condition number:', cn"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Error: 292.997678103 Condition number: 3.5650102203e+20\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "** Can you think about an explanation?**  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Todays lecture\n",
      "- Gaussian elimination\n",
      "- Condition number"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##### Questions?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "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()"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "skip"
      }
     },
     "outputs": [
      {
       "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: 160%;\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",
        "    .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"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 1,
       "text": [
        "<IPython.core.display.HTML at 0x104457990>"
       ]
      }
     ],
     "prompt_number": 1
    }
   ],
   "metadata": {}
  }
 ]
}