lecture-6.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b008d89785c15a8b0fa027155aca67a27b2bb863fb9f9fb9b443842ceb10eb9a"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Lecture 6: Matrix decompositions"
     ]
    },
    {
     "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  \n",
      "**Week 4:** Matrix decompositions: LU, QR, SVD + test"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##Recap of the previous lecture\n",
      "- Eigenvectors and their applications (Schrodinger + PageRank)\n",
      "- Gershgorin circles\n",
      "- Computing eigenvectors using power method, $x := Ax$\n",
      "- Schur theorem, $A = U T U^*$\n",
      "- Normal matrices $AA^* = A^* A$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Today lecture\n",
      "Today we will talk about the matrix factorizations\n",
      "- LU decomposition and Gaussian elimination\n",
      "- QR decomposition and Gram-Schmidt orthogonalization"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "##What is the LU factorization\n",
      "\n",
      "From the title, it is clear, that LU factorization is the representation of a given matrix $A$ as a product  \n",
      "of a **lower triangular matrix** $L$ and an **upper triangular** matrix $U$.  \n",
      "The matrix $L$ has to be non-singular with ones on the diagonal, and the matrix $U$ has to be non-singular.  \n",
      "\n",
      "**Main goal** of the LU decomposition is to solve linear systems, because\n",
      "$$\n",
      "    A^{-1} f = (L U)^{-1} f = U^{-1} L^{-1} f, \n",
      "$$\n",
      "and this reduces to the solution of a linear system \n",
      "$$\n",
      "     L y = f, \n",
      "$$\n",
      "and \n",
      "$$\n",
      "   U x = y.\n",
      "$$\n",
      "Solving linear systems with triangular matrices is **easy**: you just solve the equation with one unknown, put into the second equation, and so on. The computational cost of such algorithm in $\\mathcal{O}(n^2)$ arithmetic operations, once $LU$ factorization is known. \n",
      "\n",
      "Note that the inverse of triangular matrix is a triangular matrix.\n",
      "\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## LU decomposition: the algorithm\n",
      "We can do some **symbolic computations** using the Python sympy package for our Hilbert matrix, to illustrate how LU decomposition looks like."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import sympy\n",
      "from sympy.matrices import Matrix\n",
      "import IPython\n",
      "sympy.init_printing(use_latex=True)\n",
      "n = 5\n",
      "w = Matrix(n, n, lambda i, j: 1/(i + j +  sympy.Integer(1)/2))\n",
      "L, U, tmp = w.LUdecomposition()\n",
      "#Generate the final expression\n",
      "expr = (L + U)\n",
      "fn = '%s \\\\times %s = %s' % (sympy.latex(L), sympy.latex(U), sympy.latex(w))\n",
      "IPython.display.Math(fn)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "latex": [
        "$$\\left[\\begin{matrix}1 & 0 & 0 & 0 & 0\\\\\\frac{5}{7} & 1 & 0 & 0 & 0\\\\\\frac{5}{9} & \\frac{14}{11} & 1 & 0 & 0\\\\\\frac{5}{11} & \\frac{189}{143} & \\frac{9}{5} & 1 & 0\\\\\\frac{5}{13} & \\frac{84}{65} & \\frac{198}{85} & \\frac{44}{19} & 1\\end{matrix}\\right] \\times \\left[\\begin{matrix}\\frac{2}{5} & \\frac{2}{7} & \\frac{2}{9} & \\frac{2}{11} & \\frac{2}{13}\\\\0 & \\frac{8}{441} & \\frac{16}{693} & \\frac{24}{1001} & \\frac{32}{1365}\\\\0 & 0 & \\frac{128}{127413} & \\frac{128}{70785} & \\frac{256}{109395}\\\\0 & 0 & 0 & \\frac{512}{8690825} & \\frac{2048}{15011425}\\\\0 & 0 & 0 & 0 & \\frac{32768}{9256590525}\\end{matrix}\\right] = \\left[\\begin{matrix}\\frac{2}{5} & \\frac{2}{7} & \\frac{2}{9} & \\frac{2}{11} & \\frac{2}{13}\\\\\\frac{2}{7} & \\frac{2}{9} & \\frac{2}{11} & \\frac{2}{13} & \\frac{2}{15}\\\\\\frac{2}{9} & \\frac{2}{11} & \\frac{2}{13} & \\frac{2}{15} & \\frac{2}{17}\\\\\\frac{2}{11} & \\frac{2}{13} & \\frac{2}{15} & \\frac{2}{17} & \\frac{2}{19}\\\\\\frac{2}{13} & \\frac{2}{15} & \\frac{2}{17} & \\frac{2}{19} & \\frac{2}{21}\\end{matrix}\\right]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 39,
       "text": [
        "<IPython.core.display.Math at 0x10c6e90d0>"
       ]
      }
     ],
     "prompt_number": 39
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "A simple implementation of the LU factorization can be done in 3 cycles:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "L = Matrix(n, n, lambda i, j: 0)\n",
      "U = Matrix(n, n, lambda i, j: 0)\n",
      "a = w.copy() #Our matrix\n",
      "for k in xrange(n): #Eliminate one row\n",
      "    \n",
      "    L[k, k] = 1\n",
      "    for i in xrange(k+1, n):\n",
      "        L[i, k] = a[i, k] / a[k, k]\n",
      "        for j in xrange(k+1, n):\n",
      "            a[i, j] = a[i, j] - L[i, k] * a[k, j]\n",
      "    for j in xrange(k, n):\n",
      "        U[k, j] = a[k, j] "
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 40
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\left[\\begin{matrix}\\frac{2}{5} & \\frac{2}{7} & \\frac{2}{9} & \\frac{2}{11} & \\frac{2}{13}\\\\0 & \\frac{8}{441} & \\frac{16}{693} & \\frac{24}{1001} & \\frac{32}{1365}\\\\0 & 0 & \\frac{128}{127413} & \\frac{128}{70785} & \\frac{256}{109395}\\\\0 & 0 & 0 & \\frac{512}{8690825} & \\frac{2048}{15011425}\\\\0 & 0 & 0 & 0 & \\frac{32768}{9256590525}\\end{matrix}\\right]$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAVsAAAB/CAMAAAB2W5FLAAAAq1BMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAC+7Q8WAAAAOHRSTlMAMquZ\ndlQQQO0wRIm7Is3d72b5r/un5YGlw9/b8aHVwam3fIWxxb3J2dHp84uR13CTm/Xn4Wxg09Br0F0A\nABFpSURBVHgB7V1/g9u2DZUt2Wks+3xO2mzJmmVts7Zpl3a/uvn7f7IBJB9EQrRFiLprbrH+OIkU\n+Ai+kyUBeqKaJl76dt9u44qL2+WWWYjK5hpzYTgN3yyA3x+aZn03Qs5UlFtmGjdNZXONuTCchp/t\n7ursljUDHonb1Zn+TC7lllmoyuYac2E4DW/gxTU9ekqZy7ajZcO197umOZzpz+RSbpmFqmyuMReG\n0/AGXlzTExO6PjO3qwRrcy474TZNuWXSAQqVzQGD9cJwgJW1Ff+U4fbYCdzERrllFqiyucZcGE7D\nN1b8DLd9O0K9UFFumQWobK4xF4bT8HSnMKq6XjHm9lQMUW6Z9aGyucZcGE7DN3b8Ebe7vml26fl3\n1I2vKLfMAlQ215gLw2n4Zga+5nZ1t9ls9iXXsnLLkZ9cUdlcYy4Mp+Fnuau5vXc3ZiPoTEW5ZaYx\n3ewVd5RtrisXhtPws9zV3I5AbxWzGbhxO5u6yYY3bicpmm2guPVhMP2dBCy3zEJVNo8xBQob8c4F\ntgEr61JMxW1ps5tdAQOa29V6dTidCho2/WZLS4ll3mbT9nQrvcyy2nscrJdBHVACbt/1a8OINbc7\nOvJdsnEAvrC1dr+RMKgLNleqO+qlKw4BrwDRrXK7PrIB1leNZ+wEbkvHXGsY8IjbzarwP7Pe0VIU\nZWSHw3nMQ3HCLQsxVG4ct5SXC+thzzJbHrel3xnzVbqMuC1J3TpwTvduiq21PweX1yxKFOummTI4\nxTpjUlU14LZFz2R8ZyNuKU1e+LysabbzT5hbz63Lx1cN2zXG2LGuR0wRBHd7X5Rq8a01t1t6nLMp\nPadYTuypsxTzEq2rc9FVUzcdlzF2rMcWdTUBd9vdFae2qUPNLfuwLXpeRtSU/g9yA+N/YPcEzwmG\nMWtu93TcbguHXHcPteq6soeeuf+MqsPxirXaXV0ccFfn8tOY5pafl5Zevk0/ED0+vhvZGQ4C3T4p\nY+xYJzsXKDhcd649GE5jmtuexlz48OL04mXxVU+Pb7Xf05lrv+r7boFTbr//0nuyuV/qv5U43P/x\nBV9ajm27f/114S0qAWhum77n27iC5bQ67rbzhsI34zvq6HC3aw73BX1dN2m/Op5fvaLfW0uygD9c\nt52z96vXx/ObPzXN27brX/25HGHEbXFTPjM3hru9BNj/eDk221quvAmGFAbdx7ZfKhoRcNpAzDD0\nE++9vD2f2/a4ag5l4fG4e88t34ctsAy6j277ENyyixwzDP2UOT2f2+3xXHj2yHjiuT1vDIFKBiWq\ncrqM1eahuJWYwaL/UNxu275rC0OP092ZfyWzFsft4Uy5mt0yKQCny+jp9rH8UlPueBQzWPQfits7\n+pVuyxjrT83qeD9zKI7brbtXLLyZvs6Eu7U5HR6IW+o75L8Kb6G8sym3J3fR7kuuUFtnxP+LOUs4\nJ3Cqx3DDeLEnp8vgq+LDHLfUr48ZbPqPlNu9uzYVnVMO7vbLix8vjvniDs/tHd/bLnDcel3GpqXl\nzEnWZZchZjDqP1JufVp8VxTX3fG5dm62ZuN+IByXbUp+JNfJinQfpSHldUC1l8+RJzr7Rf0oi3wx\n4XbLFxe6vJxL7jn5ssdq7RkL3+SvuY9T288O7YZ+B93HaX9ezzxLDXCjrRU/y6GRDv2MTLIVCbeH\nswvJVn6Vtb9VljNw47acK6tlwu3WH7Bl5wRrT5+fPXP77PwsDNyfbzfXrmWigNAbZdzpVlIuax5b\nSdOLG7G1ffsirOyYwPw20eQf3Y3V6UGCmwlH/g93J+eEpnO3RtPPMpFa9Gvjs52Qs4X2wyioSP4H\noeeOMLge67qnTVEPWvPhXS/Vx6TcNkeOeSefZSK48AmSoLyIXLq6GXK20H5YBRURNjQZJyJ2Q1EP\n1qiPTGdtAkdcDK6X6mMUt9ueXouaytUgS4q1D7IK3UfOlsMxvtNHcrSweWrme6ZkJ4V32wZrYnqZ\nBFDAgYtwfVeoj1Hcpq5fKCFLirVpJCFnm2g/pk9CFzxxPfvH0veUX+RIxuGbPLqAzdUDjsvfhrBk\nV6h4mcEtsqRYRx5ccRO7Qs421n5Mn4TQWK/d2Ffu0ns8YU1GAye6ha0sOM5FpJt3hfqYGdwiS4q1\naSSSsxXtR5QctQ0cHAZO6Yk8JzxdglU4MSOmDQKOd1FcL9XH2LlFlhRr21EiOdtY+2ERByaDd2P3\nKof7oHa45zzFwtwSIrkorrMLJfoYM7fIkmJtHYlLKXLONtJ+WAQVPDBZHIfh9LLDmvYuzy27CNdL\n9TFmbpElxdo6EuRs+fdLOcYhOSqEGTY8h0eXYqf7hLBemFtxEa5zzrEkmWnmlkcO4LD2ydhCSnb7\n1Z5ztvt37f4VxdrHb/vOIqhIuuGeV3uWSG9+bPc/vKX1HT3w25o8ShDTAuHQfB8/P2+a7758dvj4\ni8Pv3/1apo9hbred7akXsqR+LcnY1K1LpdXrN19zzvYvlLr9hjj+7m3Xv//rJeur9dzz9z+QZpzk\nKh9o/Ouf2v41/buevwrp4autp3cy/o9/a5rvf+r615REWHO6+QPJVV6U6WMOdOrX8ydM91plEX7I\n9MPiq1CVhIT9WBhODY1//+wm4gXED8osW5x1TsgiFVd6MiCkqJKQcJ8Lw6lhwE3ECxa5yu/GraOF\n7kirJCQDt4vBKW6Bi3gB8UPGbFSluDVoPwymaa/+QOM6d59fIyFhkIXhGFIt7GaIFyR+UCbZouLW\noP0wmKY9CxlOSFElIWHgheFSX7kkeg+KF5L4YWya1qTcGrQfBtO0RyHDCzbcI/S5EhJGBrcLwSln\n+Vk0VyFeQPwwMstUpNwatB8GU9VtIMMLKeokJIy8MJxyFvN9IF5A/KDNcuWUW4P2w2Cq+vU39hBS\nVElIGHlhOOUs3ISe3iJXSbg1aD8Mpqm3CDUgpKiQkDDwwnCpr1SCm6KnN8hVEm4N2g+D6cjfz6Xi\nxu3D/acTbg3aD4Ppwzn/iSMztybtB8ZTIBMJpiKVKN5AH5fXxVA0POtSjg3LCz3M1n7cZCIXGB2q\nk3NCsfaD2pfKRFxX0FD4OaHjGUMS9YhBBiIYQe8h5RPnbyl3VTx/ycDFsKU0JVqmEspTfaTcFmo/\nnBNlMhE2hYYCcytHM4Yk6hHRWAyDvLgFDOg9pExaBX6jsFSfkesA/gJby1RQnupDcVuk/fD+GExD\n7IRc6DBjSKoegcYiN15dBwzoPVBGOhj5Vt2urOxjPWBrmQrKU30obsu6Nlt5X5ELZQmFnzFkrB4p\nlYEEDNF7ABPpYORbza66Bs5fYGuZipSn+nhMbtntICXzM4aM1SMWGQhjRHoPPwsJ0sHIt1ZwC2wt\nU5HyVB+PzW14982/g6LVIzYZCGNEeo/wXktIB5fqM/Lcu+NWsLVMBeWpPh6Z25AL9dfhnHqkXAbi\nMKD7IIo8ZpwOLtFnXOXW6cs6N8cMpihh/2LZyrU+HpdbvPvmZgzJqkfKZSAeY5iuyZeRDka+Nc/d\nVG043zp9JIketEwllKf6eFRu5d03N2OIVo+IxmJq5H6/n3Vk0Hu4sqSDkW8tw9JW/toLbFYYxDKV\nUHbffIBUQyNw+ZG4dXp05EIpc4eXFxP1CPPB78iVLR7j729Wh3+4g/XnX7jdHc3N0b7fNsi3lmFp\nK58T7j60+/17eubwTd/9m45h9u8X8m//su9YtjLVB3Nr1X5oT6bKOsdKXgYFa6oewTtyU3h+v8fg\n+/c3fFFr//kvrv/4Vd/9yHNzFM9f4tHiv/C3b9v+zQfC/KFvfyNdOvn3/kwak49v+/Y5B5QTffwO\n2o94GNXbw/17+B0HrUY1MAEg1kEMQUehm1cE8clUH49zTpjyYv7+4f7dc4v4ZD7i0DJgIYagHT7W\nQXwyWOa3njy3ouH23PIo8apLfsS2WsJCDEEnBT+vCOKTKSTFrUHQYTDNOVHZHJDD/fvAbYhPYFK1\nJiyJIRrEOoVyFcWtQdBhMM2NrrJ5DOl/tMKtaDVim5nbjCXxCWKdOD65hptyaxB0GExz/Vc2F8jh\n/h3cIj4Rk4oNrycJ8QlineIZT1JuDYIOg2lucJXNBXKIEQK3Ep+IyfyNgBViCMQ6Ep9MAafcGgQd\nBtOcD5XNBXK4f081IGJQsYFYx+nS/UxoLtYplask3BoEHQbT3Ogqm0eQ4f4d9/ui1YhM5m4KFsnE\n/fQkPtYplask3BoEHQbT3NAqm+cgP726G7cP9z9JuDUIOgymOecrm+cgP7065vbZF5SAcEu5oKMx\nmOZGPa85xBb2dc6HtM6OOWqRAja/fcE3x5wE5sUg6DCYeuz0b2XzFOwTLSXnBIugw6T9GA9+bvMg\n7eD8npvxGbIMmfIDBuEb21P6jMQz/1hImgimr6eD0CtVuI3bmugj5fZhtB+J/yiUS0fQgtenIO1o\n1OwhkGmIAXQmU/qMCBySDzQBJuqx5iYuHTTVh+LWIOgwmEYDkM15zZE6xSt0kGFIihUGyL0O+V3p\n+fKGD+3QRDAh+5d1yONO9aG4vdzvp7EHqdPwCh1kGEOKVQzoyQZ/e3zI7xYMIHDrH4oMmMIpUhYh\nj4tc8aU+nhi3SJ2q2UOGFCsMmEn+4U7pMxLGA7c+JTxgjrgVzcr1Pp4Yt41PncordEGGMaRYg4Gj\nLHrnLqHwYsFzG1LCESaOV6yRxyWga308MW5D6lReoQsyDEmxNlFuVfK4/td9kdFhB7ijB2Pnw4Cp\nj1vkcanh1T6eFreSOuVrmJvx2cswZMoPMdDv3A0EXtny3IaUsGDSDz9MieXXyOMSUPpen0Z+WtxK\n6hSv0IkMw3G9bcRAv3Onx50te+6QEob0Q3OLPC4JQkjtR983RwON+bS4bZA6xSt0/ssx9ACWp/7g\nFCsMkHsd8rt65JmyTwH7Jn33jnQPmzXHB/+h6T/465i0DkHL4fxr37384L9vfqmPJ8atpE7DK3Tu\nyzHMElKsMBhyr2VzdBAEUsBO0sEa9g+v+vbjoTm8eXFe/9cnD9adC1ooj/v9pvna1Tm6s1+6eGLc\nMo2PskDDjvgAunQELYhRrjlz4/YyO26ewBCDQJcegpYGMcrl1nSloxzjkAe7ZvnZ7RMNe9CSOK07\n5v2IY5RLzChuDYoMg2mu88rmGnJhOBYgyqfvgpaEdeoStMQxivYFZcWtQZFhMEVn8bqyeQzF2wvD\nOfigYQ/xgUswStASxSjaFSmn3BoUGQZT6SzaqGweIbnNheE8fPqdF/91zBC0RDGKdmUop9waFBkG\n06G3Yauy+QDktxaGi6YyFC2JP0WEoGWIUbQnUTnl1qDIMJhG3clmZXPBCRsLw3ntbfqdF69TR9CC\nGEU7EpcTbg2KDINp3B22K5sDBuuF4QgWGnaJQaB1D0ELXzunvniTcGtQZBhMQUC8rmweQ/H2wnAa\nfmb5xu1M4gqaJdwaFBkG05wXlc015MJwGn5mOeHWIuiYJ94QLyubC07YqIAbKTjqK4JTKbcGRYbB\nVBPB5crmGnJhOA0/r5xya1BkGExznlU215BVcNCPjKQc6hMvsMO3ZUQkEhrqyUZSbj957YfmVMrz\npCSuOabxGEk51CdeYCffloFIBA31ZCOKW4Miw2AqDEQblc0jJLdZASf6kfDaH7QfSNWiDDv5tgx2\nIMcrkzsH5xS32uXPoAz9iHySG1IOpGpDWexEFCKGIcerhSY3bmUaDz6OYrkIUrVBPiJ2IgqJdSUu\nxxtXENiNW/d1o1X4AnMk5ZBULV4PDDqT4V0+7OB/isvxxhU3bpmWaBqPWMohqVq24fMA7CJRiN9B\n+6Vh/DGY23FLeZkufIsnlXJAXxK0IGInohDZETQgMrkz/zdouXHbQD+ipRzQl0DaIXYQhWAHGmoR\niOfWhXn+7TRP+ef0F/oRLRdBqhbSDtiJ0AQ70DCaCOToA2d+7OYWlv18jgv0I5KqxXQemKI5lGE3\nCE3CDmk4TARy8pQ2/wOGO9tYAnlL7gAAAABJRU5ErkJggg==\n",
       "prompt_number": 41,
       "text": [
        "\u23a12/5   2/7    2/9     2/11       2/13   \u23a4\n",
        "\u23a2                                       \u23a5\n",
        "\u23a2              16      24         32    \u23a5\n",
        "\u23a2 0   8/441   \u2500\u2500\u2500     \u2500\u2500\u2500\u2500       \u2500\u2500\u2500\u2500   \u23a5\n",
        "\u23a2             693     1001       1365   \u23a5\n",
        "\u23a2                                       \u23a5\n",
        "\u23a2             128      128       256    \u23a5\n",
        "\u23a2 0     0    \u2500\u2500\u2500\u2500\u2500\u2500   \u2500\u2500\u2500\u2500\u2500     \u2500\u2500\u2500\u2500\u2500\u2500  \u23a5\n",
        "\u23a2            127413   70785     109395  \u23a5\n",
        "\u23a2                                       \u23a5\n",
        "\u23a2                      512       2048   \u23a5\n",
        "\u23a2 0     0      0     \u2500\u2500\u2500\u2500\u2500\u2500\u2500   \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500 \u23a5\n",
        "\u23a2                    8690825   15011425 \u23a5\n",
        "\u23a2                                       \u23a5\n",
        "\u23a2                               32768   \u23a5\n",
        "\u23a2 0     0      0        0     \u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u2500\u23a5\n",
        "\u23a3                             9256590525\u23a6"
       ]
      }
     ],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Strictly regular matrices and LU decomposition \n",
      "A matrix $A$ is strictly regular, if all of its **principal minors** (i.e, submatrices in the first $k$ rows and $k$ columns) are non-singular. In this case, there always exists a LU-decomposition.  The reverse is also true. "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "**Proof**. If there is a LU-decomposition, then the matrix is strictly regular -- this follows from the fact that to get a minor you multiply a corresponding submatrix of $L$ by a corresponding submatrix of $U$, and they are non-singular (since non-singularity of triangular matrices is equivalent to the fact that their diagonal elements are not equal to zero).\n",
      "\n",
      "The other way can be proven by **induction**. Suppose that we know that for all $(n-1) \\times (n-1)$ matrices will non-singular minors LU-decomposition exists."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "Then, consider the block form\n",
      "$$\n",
      "    A = \\begin{bmatrix}\n",
      "          a & c^{\\top} \\\\\n",
      "          b & D\n",
      "    \\end{bmatrix},\n",
      "$$\n",
      "where $D$ is $(n-1) \\times (n-1)$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Then we do \"block elimination\": \n",
      "$$\n",
      "     \\begin{bmatrix}\n",
      "     1 & 0 \\\\\n",
      "     -z & I\n",
      "     \\end{bmatrix}\n",
      "     \\begin{bmatrix}\n",
      "     a & c^{\\top} \\\\\n",
      "     b & D \n",
      "     \\end{bmatrix}=\n",
      "     \\begin{bmatrix}\n",
      "     a & c^{\\top} \\\\\n",
      "     0 & A_1\n",
      "     \\end{bmatrix},\n",
      "$$\n",
      "where $z = \\frac{b}{a}, \\quad A_1 = D - \\frac{1}{a} b c^{\\top}$.\n",
      "We can show that $A_1$ is also strictly regular, thus it has (by induction) the LU-decomposition:  \n",
      "$$\n",
      "   A_1 = L_1 U_1, \n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "And then\n",
      "$$\n",
      "   L = \\begin{bmatrix}\n",
      "    1 & 0 \\\\\n",
      "    z & L_1\n",
      "     \\end{bmatrix}, \\quad \n",
      "     U = \\begin{bmatrix}\n",
      "     a & c^{\\top} \\\\\n",
      "     0 & U_1\n",
      "     \\end{bmatrix}\n",
      "$$\n",
      "we get $LA = U$ and $L$ and $U$ have the necessary triangular structure."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## When LU fails\n",
      "What happens, if the matrix is not strictly regular (or the **pivots** in the Gaussian elimination are really small?). There is classical $2 \\times 2$ example of a bad LU decomposition.  The matrix we look at is  \n",
      "$$\n",
      "    A = \\begin{pmatrix}\n",
      "    \\varepsilon & 1 \\\\\n",
      "    1 & 1 \n",
      "    \\end{pmatrix}\n",
      "$$\n",
      "If $\\varepsilon$ is sufficiently small, we **might** fail."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "eps = 1.1 * 1e-16\n",
      "a = [[eps, 1],[1,  1]]\n",
      "a = np.array(a)\n",
      "a0 = a.copy()\n",
      "n = a.shape[0]\n",
      "L = np.zeros((n, n))\n",
      "U = np.zeros((n, n))\n",
      "for k in xrange(n): #Eliminate one row   \n",
      "    L[k, k] = 1\n",
      "    for i in xrange(k+1, n):\n",
      "        L[i, k] = a[i, k] / a[k, k]\n",
      "        for j in xrange(k+1, n):\n",
      "            a[i, j] = a[i, j] - L[i, k] * a[k, j]\n",
      "    for j in xrange(k, n):\n",
      "        U[k, j] = a[k, j]\n",
      "\n",
      "print 'L * U - A:' \n",
      "print np.dot(L, U) - a0\n",
      "print L"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "L * U - A:\n",
        "[[ 0.  0.]\n",
        " [ 0.  1.]]\n",
        "[[  1.00000000e+00   0.00000000e+00]\n",
        " [  9.09090909e+15   1.00000000e+00]]\n"
       ]
      }
     ],
     "prompt_number": 47
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## The concept of pivoting\n",
      "We can do pivoting, i.e. permute rows and columns to maximize $A_{kk}$ that we divide over.  \n",
      "The simplest but effective strategy is the **row pivoting**: at each step, select the index that is maximal in modulus, and put it onto the diagonal. \n",
      "\n",
      "It gives us the decomposition \n",
      "\n",
      "$$A = P L U,$$\n",
      "where $P$ is a **permutation matrix**.\n",
      "\n",
      "\n",
      "What makes pivoting good?  \n",
      "It is made good by the fact that  \n",
      "\n",
      "$$\n",
      "   | L_{ij}|<1,\n",
      "$$\n",
      "but the elements of $U$ can grow, up to $2^n$! (in practice, this is very rarely encountered).\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Cholesky factorization\n",
      "A very important class of matrices is the class of **symmetric positive definite matrices**. \n",
      "By positive definiteness we mean that the associated **quadratic form**  is always positive\n",
      "$$\n",
      "   (Ax, x) > 0, \\quad A > 0\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "In this case, the LU-factorization always exists and moreover, it can be written in the form of **Cholesky factorization**:  \n",
      "$$\n",
      "   A = R R^*\n",
      "$$\n",
      "And we always can select pivots on the diagonal!\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Summary on the LU for dense matrices\n",
      "- Principal submatrices non-singular - there exists LU\n",
      "- Pivoting is needed to avoid error growth\n",
      "- Complexity is $\\mathcal{O}(n^3)$ for the factorization step (can be speeded by using Strassen ideas) and $\\mathcal{O}(n^2)$ for a solve"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## LU-factorization for sparse matrices\n",
      "The $\\mathcal{O}(n^3)$ cost can be reduced, if the matrix $A$ is sparse. \n",
      "For example, for a three-diagonal matrix, the factors $L$ and $U$ will be bidiagonal, and could be computed in $\\mathcal{O}(n)$ cost.\n",
      "Some demo..."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "import scipy as sp\n",
      "from scipy.sparse import csc_matrix, linalg as sla\n",
      "%matplotlib inline \n",
      "import matplotlib.pyplot as plt\n",
      "plt.xkcd()\n",
      "#A = csc_matrix([[1,2,0,4],[1,0,0,1],[1,0,2,1],[2,2,1,0.]])\n",
      "\n",
      "n = 10\n",
      "ex = np.ones(n);\n",
      "lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csc'); \n",
      "e = sp.sparse.eye(n)\n",
      "A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)\n",
      "A = csc_matrix(A)\n",
      "\n",
      "fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n",
      "ax[0].spy(A)\n",
      "ax[0].set_title('Sparsity pattern of a matrix')\n",
      "lu = sla.splu(A)\n",
      "ax[1].spy(lu.L)\n",
      "ax[1].set_title('Sparsity of an L factor')"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 26,
       "text": [
        "<matplotlib.text.Text at 0x113033490>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAk8AAAEoCAYAAABB3hURAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXd8U9X//19Jmo4kHXRDW8peIih7lb1X+TK1gIAMWRWh\ngCgq8wOooCBLmYIoynAhKCLInioCAgUKQlm2dDdNm9G8f3/0dw9NkzZJm6RpOM/HIw/oveeece99\nv+/7vM/7nCMiIgKHw+FwOBwOxyLE5V0BDofD4XA4nIoEN544HA6Hw+FwrIAbTxwOh8PhcDhWwI0n\nDofD4XA4HCvgxhOHw+FwOByOFXDjicMpgYcPH2Lnzp2Ij48v76pwOBwbcufOHezcuRP//vtveVfF\nYahUKuzatQsnTpwo76pUeCqk8bRv3z6MHj0aO3futGm+ycnJePfddzF69Gj8/vvvNs3bEezatQuj\nR4/GDz/8UC7lX758Gb/99pvBLzMz06ZlJCYm4o8//rBpnkVRqVT47rvvEBUVhfDwcMTExGDatGlW\n5aHRaDBnzhyMHz8e9+/ft0m9du7ciZMnTxod//nnnzF69Gjs2LHDJuVwKi7379/H+PHj8eabb0Kj\n0dgsX51Oh61bt2L06NF4//33bZavo8nKysI333yD5s2bo2bNmoiJicHcuXOtyiMnJwfTp0/HxIkT\nkZKSYvF1Wq0WR44cQU5OjrXVLhNEhGPHjmHEiBHw9vbGsGHD0K9fP2i1WoM0H3zwAcaMGYOLFy+a\nzfPRo0dISEiwZ7WdH3IQ9+7do+nTp1NoaCiFhoZSly5d6NKlS0RElJycTH///TctW7aMXn/9dRo2\nbBjNnDmT9uzZY5TPnj17SCKREAACQH/++afBea1WS0uXLqUXX3yRlVWnTh369NNPSaVSFVu/7du3\nk1wuJ4VCQevXr6fc3FwiIvr7779pzpw5lJycbFV7jx07RqdPn7bqmqI8fvyYLl68SIsXL6YpU6ZQ\nTEwMzZw5k3766SejtFu3biWRSEQASCKR0M2bN8tUtiXo9Xr69NNPqWbNmhQSEsKeSeHfjh07iIhI\nrVbT4sWLqVGjRuy51K9fnzZt2kR5eXkG+SYmJtKjR49MltmhQwfy8PCgBw8eWFXXtLQ0Gj16NCt7\n5MiRdPbsWaN0Z8+epdDQUFb/bt260fvvv09//fWXxWWp1Wrq378/y6N169aUn59v8fXXrl2jH374\nweDYnj17CAD16tXL4Pj3339PUqmUlWWqTZaSlpZGc+bMoZo1a7L71Lx5c9q7d69R/ZOSkmjWrFn0\nzz//lLo8TgEZGRn09ddfU9WqVZm+2rZtGxER5eXl0aVLl+irr76iuLg4GjhwIMXFxdGiRYtIrVYb\n5JOYmEg1a9Zk78Lbb79tVNaPP/5IAwYMYM83NDSUpkyZUqI8JSQkUJ06dQgAjR07lqXVaDQ0Z84c\nOnbsmFXtTU5Ops2bN5Ner7fquuJ49OgRDRs2jEJDQ6ly5co0btw4unz5slG6Q4cOkb+/P7s//fv3\npw8++ICuXbtmcVlKpZI6duzI8ujdu7dBO5RKJV29etVk23bs2EEAaMmSJWbLefz4MdWrV499I0tL\nUlISNWvWjNW3Xr16tGTJEvrxxx9ZGr1eT9OmTWNpqlSpwr5/Rbl48SJ17tyZPDw8SCKR0Jtvvlmm\n+hERHT58mFq0aEGhoaFUr149evfddyknJ8co3ZkzZ2jkyJEG7+6wYcOsen62xCHG06FDh8jNzY0A\nUP369WnGjBlUv359Wr9+PRERhYeHm/zwAqCXX36ZvYjXrl0jiURCAwYMoD179pCPjw/VrVuXdDod\nERUYTkOHDiUAJBaL2Q12d3cnABQREUGnTp0yqt+CBQtIJBJRly5d6N69ewbnBg4cSABo7dq19Omn\nn9Lzzz9v9nfjxg2KjY0luVxOf/zxh0F++fn5lJGRYfae5efnk4+Pj8l7IhaLacKECSztuXPnSCQS\nUUxMDO3cuZO8vLyoWbNm7IN36tQpi+q9d+9e2rBhg0Vpr1+/Tvv372d1evXVV2nu3Ll04cIFunTp\nEl26dIlu3LhBRAXGRHR0NDPshOcifPSrVatGFy5cYO1p27Yt+fj4mLxPQUFB5OfnV6xwm0Kv11Pf\nvn0JALVq1YomTJhAMpmMFAoFnTx50iCdn58fM1JKK5RxcXHk6elJ69ato6lTpxIAWr58ucXXx8bG\nkkgkYh/Qx48fU1hYGLvXWVlZRER069Ytkkql1K9fP9q7dy/5+flRjRo1SKvVElGBvFjyLDdu3Ehp\naWlMyXp4eFBoaCiFhIQwg/zFF1+ku3fvsjouXLiQAND48eOJiOjKlSt05swZg3bodDpavnw5NW7c\nmF5//fVS3UtXJyMjg2rXrk0ASKFQ0KRJk6h79+7Up08fIiKaPn16sbqxTp06dOfOHSIq0BdNmjSh\nWrVq0Z49e6hNmzYkFosN9M+OHTtILBYTAAoICKDQ0FCmY7y8vGjx4sVGH/2zZ8+Sr68vhYSE0C+/\n/GJwfufOnQSAOnToQAkJCdSoUSOz79ratWvp+PHjBIAWLVpkdD/S09OtMqr0ej21adOGAFCnTp1o\n9OjR5OHhQQEBAfT333+zdFlZWeTp6UkAaNiwYXTr1i2LyyjM6NGjydfXl7Zs2UIjR44kALR161Z2\n/vPPPycAtGHDBqNr3377bQJAX3zxhdlyzp49SwAM9FNp6NOnDwGgGjVq0MmTJ0mj0Ril+eyzzwgA\nLViwgD744AMCQK+99ppRuj///JMqVapEAKhjx440depUUigUzNAZNWqU2effqVMngzz//fdf9n2e\nMGECde7cmT3LwgbU8ePHSS6XEwDy9fWl0NBQZghLJBKaMGEC03uOwu7Gk16vp6ZNm5KHhwfNmzfP\nZA9cJpMRABo0aBD98ssv9Ntvv9HOnTvJ19eXANDRo0eJiOidd94hkUhE169fJyKijz/+mADQ559/\nTkRER48eJQDUrFkzgw9fYmIizZo1i8LCwkihUBh8qLdv304ikYgmTpxosv7CR2LNmjV08eJFatGi\nBXl5edGQIUPIw8ODKTKJREISiYSaNGlCDx48oMTERJJIJDRt2jSD/N566y164YUXzN43rVbL8h45\nciQdPHiQfvvtN9q2bRu7XxcvXiQioqlTp5K7uzsz/ObNm0cA6LvvviMioidPnlB0dDS5ublR7969\nDTwrQr0jIiLo+PHjdOnSJWrZsiV5enrSkCFDmMIpnPbFF1+kxMREevDgAQGguLi4Etvy888/EwBq\n06YNM6iICgRnxowZFBoaSr6+vqyXJSjD9u3bGwl7UFAQhYWFmb1/hREUWtWqVen27dtERHT37l2q\nVq0aKRQKOn/+PBEVfIAEY84a46woAQEB1LlzZ5ZnvXr1KCgoiBk95hAMLplMRjqdjmbNmkVVq1Zl\nz1Xw5i1YsIAAsF72mjVrDBR3dnY2DR8+nNzc3KhLly5UrVo1o2cZHBxMP/30E61atYoA0ODBg+nh\nw4esLpcvX6YJEyaQQqGgatWqMY/gjz/+SABo3LhxRETk7u5Onp6elJCQQEQFhtPo0aNZeR07diz1\n/XRlli1bxjwYpryto0aNIgDUoEED2r59O/3222908OBBateuHQGgKVOmEBHRpUuXCAAtXryYiIhu\n375NUqmUGWFEBbIjl8tp165dzEDJzc2lNWvWUNu2bQ2uJyrwOEVERFBERAQplUqjul2/fp3JqUql\notGjR5NUKqWOHTsaeMCEdy0wMJC+/fZb0uv1FBUVRdWqVTPwnp0/f54kEolVhs2KFSuYR+Xx48dE\nVNBpCAkJoYCAAPYdSE5OJgDUuHHjUn9k9Xo9iUQiGjJkCBEVeAWrVKlCNWrUYN7zrVu3EgCSy+VG\nnQnBeNq/f7/ZsmxlPLVu3ZoAMLksLk3VqlXZt7lv374kFovp33//ZWmuXLnCDCfB6UFU4K0W+PLL\nLyksLIzCw8Opd+/e7PmLRCL2DgwfPpyl1+v17L179dVX2Tu5YcMGZkAJz6pdu3YkFotp1apV7J3J\nz8+nr776ihmII0eOLNO9sha7G08fffQRs2qLQyaTUa1atYx6HO+99x4BoF27dhFRwQ0MCAhg5+/d\nu0eenp7Us2dPIiK6evUqiUQiZjQU5euvvyYA1L17dyIiZuCEhoYW+2ETlNeaNWuMzk2aNIkUCgX9\n999/Jq8dMmQIyeVydv7Ro0fk6elZ4r0QEIynpk2bGp0TXKwHDx4kIqKGDRtSjRo12Plr166Rm5sb\nDRs2zGTegjFx5MgRs/WYOnUqyWQyppiKAoDeeOONEvP466+/CAD98ssvJs8LCic6OpqInhpPpnpp\nQUFB1KJFC7P1FsjLy6Pq1atTVFSU0TDH3bt3KSIigrp27UpET40nS4zb4rh8+TIBoLlz57JjggFe\ndCiuOAq70Js3b05ubm60efNmevLkCQFgHqlOnTqRt7c3u+7Ro0fk5eXFDLeiCJ7Cb775xuicIBuC\ncVmUt956iwDQrFmziOhpx0Uwnjp16kQAaMCAAZSRkcE+FACobdu2Bl4rTgE3btwgT09PatiwYbHe\nFkH/CB0lgStXrhAA9iEXjN/CMt2hQweSSqWUmppKRAV6IiYmxmQ5Go2GxGIxeXl5MVkXPmw///yz\nyWt++uknZjwV5ciRIwYd26Ls2rXL4F3Oz8+nbt26UcOGDU2mN0VWVhYFBgZSv3792OiDwLVr1ygw\nMJDpQMF46tKli8X5F0XonH/88cfsmCCrgpEj6DJT90WQiaLP0hS2Np6ys7NNnk9JSSG5XE5Dhw5l\nxwSP4ocffkhEBSMHUVFRBBSMwFiCXq+nmjVrUocOHYo1VoX3Z8uWLUbnPv30U4Nv/6BBg6hJkyYm\n5USv11OtWrUIAJ07d86i+tkCuxtPkyZNIgB06tQpys7OptWrVxvFKclkMpo9e7bBMbVaTXXq1KGQ\nkBBSqVSkUqnIz8+PGUoCb731FkmlUvb3888/z17axMRE2r9/P73++usUGxvLXIKC9ZuVlcVuevv2\n7U2OLwvj26ZemjfeeIPCw8OLbfv169fJ3d2d9eaGDRtGnp6elJiYWNItI6KnxlPR8fGcnBwKCwtj\nwzOZmZnk7u5uZChNnjyZgoKCTOYtCEdhD0NxxMXFUeXKlYs9LygKHx8f8vPzIz8/P/L39zcaHq1V\nqxb16NGDiAqMlh9//JFiY2MpNjaW2rdvzz7EOp2O6tatSzKZjLy9vSk4OJgpc6VSSXK53OhdKYnN\nmzeTp6dnsb3ZjRs3kpeXFz18+JAeP35s4O2LjY2l6dOnU1JSksXlCcrz+++/Z8d0Oh0FBgZaPHR1\n+PBhAkDBwcEEFAwh6vV6ys7OJi8vL9q2bRvl5uZSYGCgkaEkeKdM8dtvvxEAg+EMAaVSSV5eXhQb\nG0tERPHx8fTVV1+xZ9SoUSMCQB988AERFbz7wNNhu1u3bpGXlxcBMPDIjhkzxqp4r2cJwcB44403\nSK/X0/79+2nHjh0GH5tRo0ZRjRo1jD4as2fPNuhAjRw5kry8vAw6ckKnRXgXFy5cyDzUSqWSTp8+\nTXPnzqXY2FiaMGECAaDAwEDKzMwkoqce0LCwMIOhKYG1a9cW61U8deoUATDyvgjk5eVR69atqXHj\nxkREtHfvXgJA69ats/j+ffDBB+Tv71+sHlu2bBlVqlSJMjIy6J9//mFesLFjx1JsbCzNmjXLohAK\ngeXLlxMAOnHihEE73NzcWIdYGPYSvPufffYZSzt+/HiqVKmSRWXZwnjS6XQsVq1bt25MlguPvAjP\nSTCUBFq0aMG+tdu3b2cdeWuGVBs0aECvvvpqsedbtWpFffv2NTJ8iQoMohYtWrDOwTfffEMA6Pjx\n46TRaOj8+fO0dOlSio2NpSlTppBMJiOxWFzq4djSYHfjSfByhISEsDFKkUhE0dHRFB8fT0QFxlPD\nhg2Z0KamprIYmV9//ZWICgKwAVBUVBTdunWLtm3bRi+++CJ5e3uTWCxmbuWFCxeSVCqlzZs3k0Kh\nYEo8IiKCatasScuWLWPlEBW8/MeOHaOaNWuSl5cX+zgICMaTqZ6zOeOJiGjixImkUCioV69eJJVK\nLVYOgvHUqlUr1mtITk5mBqAQGPzDDz+wD+ytW7fos88+o+eff54ZH0W9LUS2M56EoYKoqCiKjY2l\nGjVqULt27eiTTz4x+mC+9dZbJJPJaNOmTewjKwyl1a5dm1asWEHZ2dl0+/ZtAgriCC5dukQeHh7U\nv39/UqvVdPDgQQJgVZBigwYNjIZOCyMoj507dzIFXvS3c+dOi8t76aWXmLF99epVmjRpEhvCGDNm\njEV5nDt3jgDQ77//Tl27djWYrNC9e3fatm0bnT59mr0fN2/epC+++IKaNGnCYlhMeVJLMp6ICnp3\nVatWpZUrVxpMPoiMjKRGjRrR1q1b2fskGE+FjcSVK1ca3LdBgwZxw6kEEhMTKSAggGQyGUVGRrL7\nVq9ePdq7dy8RFRhPYrGYdUa0Wi0bnhUMVyKi4OBg8vDwoNOnT9Px48epb9++VKVKFQPv7bVr1wgo\nCCRv0aIFK8/f35+qVatGw4cPNxhWJyoYuhPiSHv37m3gwRCMp6IfXiLzxhPR005Ct27dKCgoiDp3\n7mzxx1mv11NISAgtXLiw2DTC0PKvv/7KhoKK/oTviyV069aNANCXX35Jly5dotGjR1P16tUJAM2c\nOZOIiDp37kxt2rQhtVpN3bp1I39/fyZvVatWJX9/f4vKsoXx9OjRIybHhX+TJ09maYTRnZkzZ9Kt\nW7fo3Xffpbp165JEIqG2bdsSEdHcuXMJAO3bt8+q8ksynk6ePEkikchkYL/A6NGjKSAggNRqNalU\nKnJ3d6eRI0fSoEGDWFu8vb2pWrVq1Lt37xLfNXvgkIDx33//nSZOnEhTp06l3377jb7//nuSy+VU\npUoVInoa8xQQEEBt2rRhsU6FDRnBeCr8c3d3Zy/H7t27iaigx1w0XUnWr4BGo6GZM2eSWCxmrkIi\nopYtWxIAkz0US4ynnJwcioiIYAacqYA9UxSOeQoJCaHWrVszY7Bwb0Ywngr/PDw82H05dOiQUd62\nMp6++OILAmDRjKuLFy8a1VOI1yiMYDwJQ69Cb69Hjx40btw4AkDz5s0zWx5RwdCGQqEosZ3r1q0j\nDw8PunfvHun1emYsjxgxgo4dO0Znz561ygAQjKfCPyFuLDAw0KI8BOOpqOIUPG+Fjaei8iAEBH/1\n1VdG+ZoznoShlMK/lStXmkz72muvEfA0HpGI6MGDB2xiCACaOnWqRe19lklMTKQ5c+bQa6+9Rtu3\nb6fjx4+z4bITJ06wYTuxWEwtW7ZkhnjHjh0N4pAEL6XwE4vFbEJG4SEZwYMo/IKCgkzObCrKrl27\nKCwsjHkCiJ7K5saNG43SW2I8ERV0LoW6FPbomOPYsWPk7+9foudo0aJF5OvrS6mpqaTVaqlx48YE\ngCZNmkTHjx+nCxcuWOVJEYynoroWANWtW5eICoyn3r17E1HBsKzgNFixYoXBN88cthq2EwyfVq1a\n0YEDB+jkyZMGBrBgPBX+Fe7cpqWlsTwKe6wsoSTjadq0aQbvZVHUajW1bNmSBgwYwI4NGTLEqK6W\nfMPshUPWeerYsSPWr1+P1atXo0uXLoiOjkZkZCQ8PT0N0mVlZSEwMBBdunTBunXrMGvWLKO8IiMj\nMWzYMOzatQuZmZnYtGkTAODevXsAgLp166JRo0YAgOjoaABAUlISMjIySqyjVCrFhx9+iPHjx2Pi\nxInIzMxEUlISzp07V6a2y2QyjBs3DgAwcOBASKVSq/PIzMxEcHAwunfvjm3btmHChAlGaWrVqoWY\nmBh8++23yM7OxooVKwA8vS/2QFhHJjg42GzaF154AbVr14ZYLEafPn0AFKwVkpWVVeJ1cXFxWLJk\nCQ4ePMie9eDBgy2q3549e1C1alVUqVLF5Pndu3cjNjYW7733HqpWrQqRSITKlSsDAObOnYv27duj\nZcuWEIutF5NOnTph9uzZuHDhAlQqFZo3b47s7GykpaVZnZdAbm6u0RoxERERGDp0KL755htkZGRg\n27ZtAEr33AcMGAA3NzeEhoaiRYsWAAoWEiy8HoxA0bXEnjx5gp49e0Kn07FjW7ZswYULF6yux7NE\nREQEli5dik8//RQjR45EVFQUunbtCrFYDA8PD5ZOr9dDJBKhcePGeO211/DTTz9BLpcb5OXu7o7o\n6Gh89NFHuH37NhITEyGTyXDv3j0QEQBg6NChAIA+ffpALBYjPz8ft2/fNlvPIUOGYNu2bdi9eze+\n++47AMbvQGl4/fXXAQBhYWFo27atxdft3r0bzz33HHx9fU2e37hxI9577z0sW7YM/v7+cHNzQ0hI\nCABgwYIFiIqKQrNmzSASiayuc/fu3fHuu+/i0qVLyMnJQa1atZCSkmIkm3Xq1MHRo0chEokQFxeH\nnJwci3WXrRB039ixY9GrVy+0bdsWCoXCKF3z5s0xadIktg5VTEwMgAI90qlTJwDAypUrbVInIsLu\n3bvRunVrk+fVajWGDh2KhIQE/O9//2PHhXe3U6dOrA03btywSZ1KRXlYbMK4sOAGlMlk5OfnR/fv\n3y/2GsHzVLR38vDhQwJAffv2ZceEAF2tVksxMTEEgJo0acICJ0ti9erVBBQEXu7evZtZuKaGQSzx\nPBE9De4UXPGWIHieKleuXGLMjeB5Kur+FFz0I0aMMLrGVp4nYS2jSZMm0dq1a2n8+PEGv6JB5m+9\n9RZ5eXmRRqNhS0C0atXKoPdY1PMk3AshiLxjx44W9xbHjx9vcs2SkydPUteuXQkAzZ8/3+Dcyy+/\nTADYjE5reemll0wGAL/zzjsEWDbTpjjPU+GAccHzdPjwYZNphPiywpjzPBERDR48mBo3bkwZGRls\naGfw4MEGHtP4+Hjm4Tp+/Djl5+ezXn316tXpxo0bNHbsWAIKphULsxk55hGmg7/yyitE9DRgfOvW\nrSW+98HBwSZj6gS5EWRAmCG3e/du+vrrr0ksFpNMJjPwIBaHEBM4adIkSk9PZyMEmzdvNkprqedJ\nr9dTvXr1ip3cUhwDBw6k+vXrGwUj//rrr8xzt3r1aoNz3bt3JwBWr9kn0K1bN5PB8UJc79mzZw08\nTwKbNm0iAOTm5kZXr161qCzB87RixQpauHChgV7t2rWrgQewJIShVVPeQaICz5NMJjO6j0Kc08cf\nf0x6vZ5NCPn2228N0sXHx1NYWJjJb3dxnidhYs6cOXMMjufk5NAnn3xCVapUoYCAAKMYZLVaTVKp\nlBYvXkynTp0ihUJBIpGItm/fbtG9sDVu9jHJCkhISEDfvn3h4eGB3r17QywW49SpUzh27Bj69++P\n7t27s7Tdu3dHeHh4sXmlp6fDx8cH9evXNzheuXJlNGnSxOCYWq2Gn58f3NzcWE/8q6++Qrdu3XDo\n0CH4+/uDiJCQkABfX19IpVIkJSXh8uXLmDlzJurXr4+WLVvik08+AQC0b98e3t7eJutkCYcPHwaA\nUvVy+vbtW6JnJz09HcHBwahevbrB8bp166JWrVrFXmMpJaXNy8sDAKxfv54dCwgIQGBgIKKiohAY\nGGiQXnguUqkUX3/9NV566SV8++236NmzJ3755Zdie5Fubm745ZdfsGHDBowaNcri+xgTE4ONGzei\nbt26ePnll3H27Fncvn0bDx8+xMSJE7FgwQK0adPGorwsJT09Ha1atTKqY79+/bB48WKL8qhevbqR\nV6Ewfn5+SE9Ph1wuR8OGDQ3OBQQEoGXLlsXWzRzCM/L19cWvv/6K7t27Y8+ePQAKZEgqleLff/+F\nXq9HQEAAWrVqhdOnT+PSpUuoXr06fv/9d0RGRmLNmjW4du0azpw5g379+uHx48elev9dmejoaFy6\ndAm9evWCv78/Hj16hJ07d8Lf399g1Ws/Pz8MGDCg2PuXl5cHlUqFVq1aGZ3r168fTp8+zf5Wq9Us\nz65duwIokJPevXvjwIED6NChAwDgv//+Q05ODiIjI5GQkICcnByMGjUK7u7uiImJQWpqKts9oG/f\nvkblWqpjHjx4gBs3buCFF16wKL1ATEwMBg8ejOeeew4DBw7EsWPHkJiYiMePH2PatGn45JNPjL4L\nZYGIkJ6eji5duhid69evn4EOLMrYsWNRqVIl+Pn5oUGDBlaVGxcXBwBGIzU9evSwKp/iSE9PR7Nm\nzeDmZmgKFH6mIpEIu3btQpcuXTB48GBER0djyJAhePToEd5991288MILCAsLM7hep9MVO6ogFosR\nExODZcuW4fLly6hXrx6+//57ptOmT5+OMWPGICAgwOA6tVoNIoKfnx/atGmDgwcPokePHhg1ahQA\nYOTIkba4JZZjT8ssJyeHVq9eTe3bt2fj7/Xr16d58+YZBDLLZLISA/+ICrwI9erVM3lu5cqVBuOn\nPXv2NIjVKOyB6tGjB+n1ekpLSyORSESenp7k6+tLMpmMqlSpQkuXLmWekKVLlxp5tQoTERFB3bp1\nM3sfFi1aZNJqLwnB8/TJJ5+UmG7o0KHUqlUrk+fmz59vsJimwJAhQ9gsRnNUq1at2Knvhw8fpqVL\nl9L58+fpwoULdOHChWKnxBIVBN9/+umn7O/CHihhbNuU56ksnDhxgrp160ahoaH02muv0bZt20pc\n/HL58uVUt25dg/VLLEV4ZqtWrTJ5rkGDBvT7779blNfy5cuNPIN6vZ6uXLlCer2epkyZQjVr1jR5\n7dq1a+n//u//jI5PnjyZvL296cmTJ8WWGxkZyWZwEZGBB0qY5Sis2yV4XfV6PX322WdGkyr+/vtv\nNkmkuKUunmUOHjxIw4cPZ7GMfn5+NG7cOIOgbWHhwZL49ddfCQBduXLF6Nzdu3dJLpczz+H69euN\nArMLe6CEtaaE+J7KlSuTRCIhf39/eumll5hHNiEhgXnlTa2JNmPGDDaLtSSEGNWXXnqpxHRF0ev1\ndPDgQYqKiqLw8HCaPHkyff311yXurDB37lxq3Lix0Y4GlpCUlETA0zXWCqNSqSg8PJwuX75s0vNU\nGvLy8ujFfxPvAAAgAElEQVSdd96h2bNnWx2bVZijR49SlSpVivUA1qhRw2TsKRFRv379DOJrs7Ky\n6MMPP6Tw8HDy9/cnf39/+r//+z+Tz1h4rqYmExA9tQ2qV69OLVq0oBkzZtCxY8dK1E2HDh2iBg0a\nGEz4KuyBKsvuCqXBYcN29+/fLzay3tfXl5YtW1bi9ceOHaN33nnH5DmdTmd2SE6r1dK6desoPDyc\nfeD/+OMP2r9/Px04cMCkIbFixQqSSqXFfshv3Lhh0fo1y5YtIzc3N6u2+NDpdCQWiw2MDVP8/PPP\nxS73r1arTQZUJiUllTh0U5ibN28aLJZmazQaDa1cuZKqVq1KeXl5lJiYaFPjydG8++67xRpIT548\nsdmWFCdPnqS33nrL5Lni5CE9Pb1UQ2iZmZn05ptvUuvWrYno6SypwgsqFsedO3csGhJ6llGr1XTk\nyBGTExPGjRtncq23wqSlpdG4ceOKHd63ZJjq7NmzVK9ePbYW2YMHD2j//v20f/9+k/J/7949Agpm\nkJp6pzMzMy36mAmdpffee89s2vJEp9NRXFyc0Y4RAsI97tGjh02MJ0exbt26YmcTp6enWzzByRSn\nTp0yubiqrYmPj6eoqCirdnGwBSKi/x9J6CCIyMj9fO7cOTRo0MDk0Jgjyi8pLRFBLBaDiKBSqUwO\np+j1euh0Ori7u5vMJysrC/Hx8SwI11JOnTqFJk2awMvLy6rrygIRQalUOuRZFC1XeC7r1q3DgAED\nig30Lsz58+dx7tw5TJ061ei5qtVqSCQSI5d0RUSj0YCIDIKIHU3hZ6TX60sVSF+RKQ/ZePToETIy\nMqwe7ikKEUGj0cDd3d1IToTnao1uBGz3Dhw5cgQdO3a02fukUqng4eEBiURik/ys4Z9//sE///yD\nl156yaL0KSkp+PzzzzFs2DBEREQYnCMiaLXaYr8rFQl7y461766tCrU7er2edu/eTQMHDqTQ0FC2\nqmxFQaVSUZ8+fUgsFhts9qvVamnz5s0UEBBAlStXNhjucDZ++ukntg8WUcGQStElBs6dO8cWrHzl\nlVdK5d52JEIQtrBivEB+fj69//77FBgYSEFBQU7vxbp27RoNGTKEoqKiaP369QY9+fz8fProo48o\nODiYAgICTK4OzrE/58+fZ8tYDB8+3OllozBKpZItCVB0+YiDBw9S+/btaevWrZSdnW0zz6it+fff\nfw02RE9KSjKShYyMDJo6dSp5eHhQ3bp1TQ5jOhPnzp0jHx8fg621BI4dO0YtW7YkkUhEkyZNcvi+\nbdaQl5dHb7/9NrVr147Gjx9v5Ok8deoUC+IfO3ZsmbxZzoRDjKf333+fjY+PGTOG/P39LdpZ2hnY\nt28fSaVSatKkCUmlUpoxYwY7V3iNksmTJ5NcLjc5Ju4MtGjRgu2ynpubS0FBQQaxOdevX2frbTVv\n3pzatm1LzZo1c9qPhDCDcdCgQUauYWH1ZaBgg0tvb2+D/ZiciU2bNpFcLqfAwEBq0qQJATAYwi68\nDsv48ePJ19e32LWXOPbh5s2bbFPSpk2bUvv27alJkyZOKxuFSUpKooYNG5JMJqPnn3/eYNeBo0eP\nkkQiIaBgocrnnnuO+vfv75QLm65evZoiIiLY31OnTqVmzZoZpBGMWy8vL5owYQLJ5XKTG8E7Cy+8\n8AJVqlTJaGHSP/74g22W2759e2ratKnBPm/OREJCAtWtW5cAULt27cjT05MaNWrEVg2/fPkyW+eu\ndevW1KpVK2rTpo1LGFB2N550Oh3VqlWL1q1bR/fu3SOdTkcXLlwghUJh0d5q5Y2wrYZWq6XevXsz\n40mYqnvgwAEW97Rv3z7y8vIyEgZnIDo6mhlPKpWKgoKCDLY5mTt3LnXu3JkyMjLov//+I6VSSW3a\ntDG5u7Yz8OabbxIAoxgElUpFISEh9M0339Ddu3dJr9fTkSNHyMvLy6qYM0dw4MABAkAtW7akW7du\nkV6vp+bNm5OnpyfdvHmT1Go1hYeH0+eff87acvLkSZLJZA5fTfdZZt68edShQwfKyMigx48fU05O\nDkVFRdHYsWPLu2pmmT59OtWpU4cePXpEO3fuNDCeRo4cSa+88gqlpaVRSkoKJScnU506deh///tf\nOdbYNHv37jVYFmby5MkGiw5fvXqV3Nzc6ObNm2z7qxUrVlBQUFCJk1jKC2EPTFO7H0yZMoX+7//+\njzIyMig5OZnS09OpcePGVm1L5QhUKhVFRERQQEAAWy7gk08+IeDptmIzZ86kXr16se9KVlYWNW/e\n3OKtqpwZuxtPP//8M/n5+RkFLjds2NDkfknOxvXr1yknJ4eePHlCXl5ebGhu7dq1VK9ePSMLWiaT\n0bFjx8qjqiXSpUsXZjwRFcwAFDbE1Wg0FBgYaOQG/9///mc0JOYM3Lp1iyQSCbVp08bo3I4dOygi\nIsJoBlDlypXpxx9/dFQVLUJYMVfYKFSr1bLtd/7880/as2cPhYSEGHnWqlevbrAKPsd+aLVaCg0N\npS+//NLg+Pvvv1/sLFRnIi0tje7cuUNarZaGDBnC1n3LyMgwCkMgKlg5Xtjs2Zn4/PPPDYynu3fv\nEgA2O/D111+nPn36GFxz//59AlDsxu3lSVRUFInFYqPZgSqViry8vIw2UZ89e7bFazs5CmH7m8Kj\nMcKaVtOnTye1Wk2+vr4GWzgRFXjT+/Xr5+jq2hy7R3weOXIEXbt2NVjD586dO0hMTETnzp3tXXyZ\nqVevHmQyGUaMGIGhQ4eytamOHDmC/v37G6wYfv78eUgkEqsDwx3B0aNHDdbNkEgkyM/PB1AQ5JiS\nkoKBAwey83q9HidOnHDKZ5SSkoL8/HxMnjzZ6Nzhw4fRp08fgzVRrly5goyMDHTs2NGBtTRPbGws\n5HI5OnfujDp16iAyMhJHjhzB7Nmz0aRJExw5cgQ9e/Y0mKQQHx+PpKQktuovx75cv34d//33n4Fs\nEBGOHz/ulLJRlEqVKqF69er45ptvcOjQIaxduxZAgT4ICQkxWOU5NzcX58+fd8p2FdVfwiQQQYcd\nPnwYgwYNMrjm2LFjeO655yzaAcHRPH78GN26dUPt2rUNjp89exZisdhgHSetVotTp06ZXGOqPOnQ\noQNatWqFVatWoXbt2qhduzYmTJiA2rVrY+HChfjrr7+QnZ2N/v37s2t0Oh1OnjzplO+YtTh8uoxa\nrcbEiRMxYsQIVK1a1dHFl4qzZ8/ixIkTGD9+fLFp0tPTMW7cOMybN89o2xlnoVmzZgAKZqMkJCTg\n+PHjeO+990ym3bhxI65du4ZJkyY5sooWceDAAQDAmDFj4O3tjVq1auGzzz6DXq83SqtUKjFu3DjM\nnDnT4TMIzREVFYXY2Fikp6fj1q1bSElJAQD89ddfJtuiUqkwYcIETJkyxWgBUo7j2LJlCy5duoSp\nU6eWd1UsIicnB2vWrEFMTAx8fHyKTbd48WJIJBIjI8RZaN68OYACYyI+Ph4A8PLLL+Px48dGae/c\nuYMZM2Zg2bJlTrc4682bN5GQkICDBw/C29sbvr6+iI2NxX///Wcy/UcffYS0tDTHLwJpBh8fH2ze\nvBn5+flISEjA/fv3odfrcffuXTx58oSlK3z/161bh3v37rEtyyo09nZtrVu3joKDg+nEiRN09uxZ\n6tq1K/Xt29eizSidAb1eT02bNjUab545cybVrVuXLl68SEeOHKFmzZrR+PHjWaCcs9GuXTsKCwuj\ndu3aGey0PXr0aHr8+DF5eHjQRx99RPHx8fTee+9RREQE/fnnn+VdbSNUKhWFhoYSAPr5559ZXQHQ\nmjVraNGiRRQZGUnnz5+n48ePU5s2bSgmJsZgUVZnQa1WU9WqVVnsVlZWFtvWZPny5fThhx9SWFgY\nnTlzhk6ePEkdOnSgwYMHm1yUkGMfkpKSyMvLiz788EOKj4+nBQsWUHh4eIXacubjjz+myMhIAxm4\ndOkSicVi+vLLL+nq1as0adIkqlevHt2+fbsca1o8n332GUmlUmrfvj0FBQUx/VW3bl26efMm9evX\njzp16kTXrl2j77//nmrVquWUsVtERLGxsWy4a+PGjWyrqJYtW1JCQgJJJBLauHEjXb9+neLi4qhG\njRolLuxbngiTc6ZPn046nY72799PAKhx48aUmJhIUqmU1qxZQ/Hx8fTWW29RZGSk0bYrFRW7G08q\nlYqGDRtGAEgmkxW74qizsmTJEgJAbdq0oZkzZ9KsWbPo9u3blJKSwvb78ff3N4qJcDZ2795NNWvW\npI4dO9LixYupU6dOVKNGDWbsbdiwgXx8fAgA9enTx2AVV2dCCAgvbPcL+7x9/PHHlJmZSb179yYA\n5OPjQxs2bCjH2paMEBAOgO7fv08qlYpat25NlStXphs3blB2djZFR0cTAFIoFLRmzZryrvIzyebN\nm9k+bj179qT09PTyrpLFXL9+neRyOfn5+VFcXBzNnDmTdu/eTUQFe4B6eHiwTpQzdjAE0tLSqH79\n+tS2bVuaMmUKzZw5k+2tSFSwYHHDhg0JAEVERDhl3KnAlClTCAALstbr9VSzZk1q0qQJERUEugsz\nn4cOHerUjoa4uDiDAPEVK1YQAFq7di0REa1Zs4atoD9gwACnDN4vLQ5bYfzvv/92ysA9c7z55psk\nk8moV69e1L9/f1IoFEz56PV6On36tMlNg52dJUuWUHBwsIGnLC0tjS5cuFCOtbKMt99+m3lnFixY\nQLVr16YGDRoYfNTOnz9fqi1WHM3Ro0fJx8eHpFIpeXp6UuXKlSk+Pt4gzZ9//kkpKSnlVEMOUcHs\n2nPnzpV3Nazmzz//JIVCQY0aNaKhQ4dSZGQkDRo0iJ1/8OCB0XpvFYEHDx4QADpw4AA7lp+fT8eO\nHXNqI5Do6ebM0dHRtGnTJurZsydJpVL6+eefWZqkpCQ2kcSZycjIYJtPe3t7sxGAwqSkpDjlKEZZ\ncfgK4xzn4Mcff8SBAwewfv16p4sJMEdWVhYGDhyIw4cPw9vbGyNGjMC8efMQEhJS3lUrFenp6fjr\nr78AFExQKLrJJofDMUSpVGLAgAHYsmVLhYmdLcyHH36Id955BxqNBr169cK0adNsttmvoyEinD59\nGnl5efDz80PTpk3Lu0oOwWbG0507d3D27FkMGDAAMpnMFllyOMVCRLh58yaCgoLg7+9f3tXhVHBy\nc3Pxww8/oEmTJqhTp055V4fzDPDff/8hNzcX1atXL++qcEpBmY0nvV6PN998EytXroROp0OtWrWw\nb98+1KtXz1Z15HA4HLtx6NAhjBo1Co8fP4ZYLMaGDRswZsyYZ27vPg6HYzll3jF19+7dWLNmDaKi\notCuXTsolUp07twZ58+fR3h4uC3q6BRotVoolUrk5uYiOzsbOTk5UKlUSE9PR2ZmJrKyspCeno6s\nrCzk5uYiNzcXGo0GeXl5UKvV0Gg00Gq1yM/PZ9PQxWIxpFIpPD09IZfL4eHhAXd3d/j4+MDHxwcy\nmQxyuRy+vr7s5+fnB7lczs7JZLIKN+xmDiJi9zMnJwdKpZLd85ycHOTm5iIrKwsZGRnsnFKphEaj\ngU6nQ35+PvLz8w3+LyDcK5FIBIlEAqlUCjc3N7i5ucHd3R0ymQxeXl5QKBRQKBTw8fGBr68vvL29\nERgYCD8/P/j5+SEoKAheXl4V/t4TEdLS0pCamoqsrCxkZWUhMzMTGRkZSElJQUZGBrKzs6FSqdj7\nLLzbOp2OvctCH0wsFkMikUAsFsPNzQ379u1zaiNEp9MhOjoaDRs2RKNGjTB58mTExMTAzc0No0aN\nKu/q2RSuwxwH12GOw5467Jdffin2/pTZ8yQYTUuXLmUVGDhwIHx8fLBt2zbk5+dDLBajZ8+e8PDw\ngFwuh4+PD+RyOQICAhASEgJ/f3/4+PjA09MTvr6+CA4ORqVKlWy+mzQRMcUh3OCUlBT8999/SEtL\nQ3Z2NjIyMpCUlISkpCQolUqkpKTgyZMnyM7OtqgMd3d3eHl5wcvLC+7u7vD09ISnpyfc3d3h5uYG\niUQCiUQCIoJer4dWq0VeXh5UKhXUajXUajWysrKQl5dncbvkcjn8/f0RGBiISpUqISAgAIGBgahS\npQr8/f3h5+fH/vX09IS3tzcUCgU8PT2hUChsvvu4RqNBRkYGUlNTkZmZCZVKhZycHKSkpCAzMxNK\npRJPnjxhikOpVCI1NRWpqalMeVv6WspkMtaeove48P8LC4Bw7/Pz86HVaqHT6aDT6aDRaKBSqZCX\nlwelUgm1Wl1i2V5eXggMDERgYCCCgoKYQlIoFKhUqRJ7jytVqgQfHx+myGz9bhMRe3eys7ORmpqK\n5ORkPHr0iCkU4R4nJycjJSWFKfXU1NQS2ymRSKBQKCCTydi7Xfh9FovFBsZRfn4+3NzcoFarodVq\n8ccff9isnfZgw4YNWLRoERITE9k78t1332HIkCFQq9VsMVmuw7gO4zrs2dBhEokEarUaOp0O586d\nKzbfMnme7t27h3PnzmHHjh3smEgkQmRkJBOc/v37Izk5GSdPnkTHjh1x/fp1ZGVlMQu8JNzc3Fij\nPT094eHhAalUavSCiUQi9jLp9XrodDpotVpoNBrWu1Kr1VAqlQYWvCkUCgWCg4MRGhoKPz8/1KxZ\nk/3fx8eHWfRyuRxyuRx+fn7w9fVlL5SHh0dZbilDq9Uygc3MzDToGQo9F6VSyXo1wsuUmpqKK1eu\nICkpCenp6RaVJbxI7u7u8PDwgKenJ6RSKetRuru7M+VEBTM0WW9IsOJzc3MN7rU5vL294e/vz3pH\nQUFBaNCgAbufvr6+8Pf3h1wuh7e3N7y9vVlP1cvLiy0uZ2ulWRidTsc+UFlZWXjy5AkyMzPZfX7y\n5AkT5tTUVPz7779ITU1FdnY2tFptiXkLQqxQKODl5QWJRMKOSaVSiMVipiwFJSm818Lz12g00Gg0\nUCqV0Ol0JZYnk8mYMgwMDET16tXZxz8sLAyBgYHw9vZmvVQ/Pz8EBgZCoVCUqmfap08fJCcnQ6fT\nsdWgnZGdO3diypQpBm2MjIw0WOiW67DSwXUY12EVVYdZor/KpNVu3rwJuVyOyMhIdiw5ORk7d+7E\nV199BaDg4YnFYnh4eODMmTMG16tUKqSkpLCHpVarkZ6ebmTRC5a00Jst7NoUXG4ikYhZkW5ubkxB\nCW5hDw8PZjX7+PiwmxwYGIiQkBAEBQVBLpdbpOg3bdqE5ORkqFQqzJ07tyy3sFikUikTwCpVqpQq\nD7VajYyMDOa+zMzMZD2S7Oxs9v+cnBzmller1cjLy4NWq2U9So1GY+CmBwAPDw+4ubnBw8OD9VKF\ne+3j4wM/Pz/WY5HJZKwnU6lSJdbDsjXBwcHsvcjMzCxzfm5ubvD39y9VQLrQe0pLS2M9UeG+Fx4a\nEZSI0GvMzc2FVqs1WGFcuM/Cey0oX+FjIbzXwjvu7++P4OBgVKlSxar3ujQsWLAA8+fPZ39LJBLk\n5OQwuXdmwwko2H4lLi7O4NjKlSsxePBg9lHjOqx0cB1mHeHh4ey9SEtLs0meXIeVzMKFC7FgwQKD\nuubl5Vmkv8pcG2EMXCKRQKVS4eWXX0aPHj3YPjx6vZ69rMHBwZBKpczlHRISgsDAQPj7+6N///7o\n1q1bWavjEApv02IvxWMLPDw8EBISgpCQENStW9dkmqNHjzLF3LhxYwfX0LYU3hKgvBF6mtWqVbP6\n2sIy4+wUHZoQhrgqUhsKexnWrl2LX375BdeuXWPHuA4rP8zpsH379jGPjjPuKWoNDx8+LO8qGPAs\n6DDB01sYS/VXmWKeMjIyUL16dXTo0AEtWrTAjh070LVrV3z44YfM9du5c2doNBqcPHkSbm5uyM/P\nt3gs2Bo6dpyPo0fn2zxfUxQde67I8LY4H67QjsJy78yMGjUKR48excSJE/H333/j3r172LRpExo2\nbMjSOEqHORJXeMcA12kHwNviTFiiv8rkefLz88Phw4fx5ptvYuvWrViyZAmGDBlikEav1zP3d+Gx\n+tGj5+Pzz+ez/9+9a768atXArrEEe+XL4XBKprDcOzNr167Fe++9h48++givvPIKtmzZArlcbpDG\nlA7juoXDcV0s0V9lHrZr0qQJDh06ZLWbrrBCsZdyuXsXDvNGcTiciodCocBHH32E5cuXl1p/cTic\nZw+bDUoWp3gq+hoSHA7Heiqa3JdkOFW0tnA4nLJhicw7JKKr6CwHoMDt7UgcXR6H86xTNBCzImNK\nh3E4HNfFnP6y+zxiIXIdAObNm8f+b0m8gC2xZXnz5s2zXWblDG+L8+EK7Sgs9xWd4nRYRcYV3jHA\nddoB8LY4E5boL5ttDFwcPXr0QEZGhtFKnaZmx5UmiFzIo/C1psoo/H8e7Mnh2Jfi5L4iYqotRXWV\nvSe/cDgcx2GJ/rK758nd3R0ajcaitGUJIrcmPVdaHI59sUbunR1TbSlOV3HdwuFUfCzRX3YfwPfy\n8rJqjyMOh1PxcSW5d6W2cDgc81gi83Y3njw9PZGbm2tRWh5EzuG4BtbIvbNjqi2FdQfXIxyOa2GJ\n/rK78WSN+74iB5FzOJynuPqwXWHdwfUIh+NaWKK/7B7zJJPJoFKpSn19WYLIrdmShwd6cji2o6xy\n70yUpi08iJzDqbhYIvN2N568vb2hVCpBRKVabI4HkXM4FY+yyr0zYa4tpnYx4EHkHE7FxRL9Zfdh\nOz8/P+Tn50OpVNq7KA6H4yS4kty7Uls4HI55LJF5uxtPgYGBAICUlBSD45YOqdkqGNPR5XE4zzLF\nyX1FxFxbuM7gcFwLS/SXQzxPAJCRkWFw3FJXtq2CMR1dHofzLFOc3FdEzLWF6wwOx7WwRH/ZPebJ\n19cXAJCZmWmzPB0VgMkDPTmc0mEPuS8vStMWUzse8CByDqdiYInM29148vHxAQBkZWXZLM+7d00H\nadoarrQ4nNJhD7kvL0rTFlO6gweRczgVA0tk3u7DdkIlsrOz7V0Uh8NxElxJ7l2pLRwOxzyWyLzd\njaeAgAAAlgWOWrMuky3gQeQcjn2wRu6dHVNtMac7uM7gcCoulugvuxtP/v7+8PDwwMOHD82mdbQr\nmweRczj2wRq5d3ZMtcWc7uA6g8OpuFiiv+we8yQSiRAQEIDU1FQ8fJiG4cM/KTZtWYIjHRWMyQM9\nORzzFJb7io6ptpgKCDcHDyLncCoGlugvuxtPACCXy5GTk4OwMH+7BXo7KhiTKy0OxzIEuXcFiral\nNHqAB5FzOBUHc/rL7sN2AFCpUiWkpaU5oigOh+MkuJLcu1JbOByOeczJvEOMp8DAwFIHjvKVyDmc\niklZ5N7ZKKktXGdwOK6HOf3lEOMpODgYycnJpbqWr0TO4VRMyiL3zkZJbeE6g8NxPczpL4fEPAUE\nBCA9Pd3gWFmCI63p6fEgcg6nfDAl9xUVe7WFB5FzOM6JOZl3iPGkUCigUqmQn58PiUQCoGyrhFuq\nJHgQOYdTfpiS+4qKvdpiyjDiQeQcTvljTuYdMmwnrNapVCodURyHw3ECXEnuXaktHA7HPOZk3iHG\nk7Ba55MnTxxRHIMHkXM45Ud5yb09sLYtZdktgesRDqf8MSfzDjGeqlSpAgB49OiRI4pj8CByDqf8\nKC+5twfWtqUsw21cj3A45Y85mXdIzFNwcDAAmJ15Y++gblOBmTyInMOxD5bKfUWgvNvCg8g5HMdi\nTuYdNtsOgNlF5soSRG4Oc4GZtoYrLc6zjqVyXxEo77YU1l9ct3A49seczDtk2M7b2xsAD7bkcJ4l\nXEnuXaktHA7HPOZk3iGeJw8PDwCAWq22+BpLXdOCp6o0AZqlcYVbWx53sXOeVUoj986KLdpSliBy\nU3DdwuHYD3My7zDjyc3NDVlZWRZfY62wl3WjTnuVZ8+hSA7HmSmN3DsrtmiLrQ0YbhBxOPbDnMw7\nZNhOJBJBoVAYuL9s3QvjcDjOhSm5r6jYqy1cD3I4zok5mXeI8QQA7u7u0Gg07G9Le032XPPE0eup\n8PVbOM8aReW+ImOPtnDvEYfjvJQk8+VmPFmKPdc8cfR6Knz9Fs6zBjeeOBxORaUkmXdIzBMAeHl5\nITc3t0x5uEoQuTV582BPTkXGFnLvLDiqLaXd+JyvA8Xh2JaSZN5hxpNcLkdOTk6Z8nCVIPLS5M3h\nVERsIffOgqPaUpodDfhmwhyO7SlJ5h02bCeVSqHVaktMY8/gSR6YyeE4HkvkvqJQnm3h+ovDcTwl\nybzDjCc3NzfodLoS09gziNwZAtSdoTwOx5FYIvcVhfJsC/cmcTiOpySZd5jxJJFIkJ+fb5O8eBA5\nh1MxsKXclzeu1BYOh2OekmTeqpins2fPYtWqVXj55ZfRv39/dlyr1WLbtm04deoUpk6diqZNmxpd\nKxaLQURWVt08jgiO5EHkHE7psJfclwalUonZs2dDLpfjgw8+gEgkYucuXLiA9evXo2XLlhg/fjzE\nYuN+pTO1pSjFbXzOg8g5nNJTksxbbDytXLkS77zzDnJzc3Hu3DlmPKWlpaFr1664ePEiPD09sXfv\nXqxfvx7Dhw+3Te3N4AhB50HkHE7F5vLlyxgyZAgSEhKg1+sRHR2Ndu3aAQDeffddLF68GACwe/du\nHD16FFu2bIGXl1d5VtkqzG18znULh2NbLBq227p1K+Li4jBnzhysXr0aKSkpBueUSiWWLl2Kffv2\n4YsvvsDkyZNx5MgRgzyctcfG4XDshzPI/ePHj9G5c2dUq1YN586dg5eXF9NhycnJWLJkCebOnYsV\nK1YgISEBly9fRlxcnFE+ztAWDofjOEqSeYuMp4iICKxatQrvvPMOcnJyUKVKFXZu9erViI2NxZw5\nc9C1a1dER0fj5ZdfxpYtWwzy0Ov1Jl3h5rBmzRNnytsZyuNwypvSyr0tkclkGDhwIH788UeEhIRA\nq9UiLCwMALBmzRo0bdoUixcvxowZMxASEoKNGzdi/fr10Ov1Bvk4Q1s4HI7jKEnmLdIEXbt2xdSp\nUwwKFs0AACAASURBVHHz5k2sWLECH3/8MQDgjz/+QHJyMsaMGWOQXqvVokaNGkbHpFIp8vLyIBKJ\njH7z5883WXZp1jyxFHvm7QzlcTj2Zv78+SblWRjyEuS+PPH19cWGDRvg7u6O2NhYDBkyBM2bNwdQ\nMEw3depUg/RarRbVqlUziIkSjpdGh3E4HOekLPrLqoDxjRs3IjIyEr169QIApKamwsvLCwqFgqW5\ncOEC9u7dizNnzhhcm5ubi8qVKyMvL8+qxpUGZw0iL015PNCTUxHx9PQE8FTunYF//vkHP/zwA27c\nuMGOpaamIjAwkP2t1Woxbdo0vPbaa0bGkyN1mCWY82TzIHIOp3RYor8sNp6OHz+O9evX4/Dhw+yY\nVCqFSqXCkydPEBQUhKtXr6J///5YtGgR6tevb3C9RqNx2N5Qzh5EXtpyOJyKgru7O4Cncl/eEBGi\no6Mxc+ZM1K5dmx2XSqW4d+8eACAvLw/Dhw+HQqHAjBkzjPJwpA6zBHO6wZRhVNyq5BwO5ymW6C+L\njKe8vDwsXLgQ/v7+CAgIwO3btxEYGIgOHTqgbt26aNKkCWrUqIErV67giy++QJ8+fYzyyM/Ph5ub\nG4KDg3ngJYfjIsyfP7/E4SpB7subjRs34t9//0Xz5s1x584deHp6IiwsDG+88QamTZuGXbt24fr1\n6xg+fDi2bdtmUmFyHcbhuBZl0V8WxTydOXMGhw8fxv3791G7dm3UqlULPXr0gEQiwYkTJxAbG4vn\nnnsOly5dMmk4AQXuL1tN/eVB5BxOxcCWcl8WFi1aBAAYNmwYatWqhfDwcNy4cQMzZ87EDz/8gNDQ\nUGzfvh3Lly83CEMojLO0pSwI689xOBzzlCTzFnUJ27Vrhy+//BJVq1ZFgwYN8Mcff8DPzw8A4O3t\njdmzZ5vNQ6lUQi6XW1Ht4uFB5BxOxcCWcl8Wdu7ciUePHqFz585ITk7GlStXUKtWLYhEIvTq1YvF\ncZaEs7SFw+E4hpJk3iLjSSqVIiYmhv3dvXt3qytR3r02HkTO4Tie8pZ7AWFBTAAIDAxEgwYNrM7D\nWdpiCms82DyInMOxjDJ7nsqKRqOBRqOBj4+PI4ozyd279ndZ8yByDucpziD3tsLZ22KpTjCVjq9E\nzuEYY07mHbLiW3p6OgCgUqVKjiiOw+E4Aa4k967UFg6HYx5zMu8Q4ykrKwtAQXyUvbBmA1575c2D\nyDmcpzhC7h2Fq7TFlM7geoTDMcaczDvEeFIqlQBQ7CwWW+AMw2Q8iJzDeYoj5N5RuEpbzK39xOFw\nCjAn8w6JeRIsOGeLF+BB5ByO/XBWuS8NFbEtZfHG8yByzrOOOZl3iPHkrPECPIicw7Efzir3paEi\ntqUsOoEHkXOedZwi5kmohL+/vyOK43A4ToAryX1FbEvHjvPLuwocToXFnMw7xHhKTk4GAAQEBDii\nOAYPIudwyo/yknt7UBHbYk/9x+G4OuZk3mGeJ6lU6vCZKs4wTMaDyDnPKuUl9/agIraFD7dxOKXH\nnMw7JOZJqVSWu9KxxnNjz6BuU6v78iByjiviDHJvK1ypLZbAVyLnPOuYk3mHGE9Pnjwpd3e3NZ4b\newWRm1I2zuAd43DsgTPIva1wpbZYgjldxXULx9UxJ/MOGbbLyMioUIGWHA6n7LiS3LtSWzgcjnnM\nybzDPE+VK1d2RFE2xVLXtOCpKk2AZlnc384wFMnhFEdFlXtTuEpbTOmMsm4qzHULxxUxJ/MOMZ5S\nUlLQuHFjRxRlU6wV9tIoh7IMEzrDUCSHUxwVVe5N4SptMbcxcEk4OuyAwylPzMm83YftiAgpKSlO\nGS/Ap/JyOPbBmeXeWlypLabgy5twOIZYIvN2N55yc3ORm5uLoKAgexdlNZb2mrhy4XCsw5nl3lpc\nqS2m4N4jDscQS2Te7sbTkydPAFSsxeWKwtdO4nCswxXkXsCV2mIKvhI5h2OIJTJv95injIwMABVr\nT6iScFQQubXrqVhbHg/05NgTV5J7V2qLKWy1gXDhY1y3cCoylsi8w4wnX19fexflEBwRRF6W9VR4\nEDnHGXAluXeltpiiLAaMKSOJ6xZORccSmbf7sF1KSgoAVJh4AWfYD4/DqehUNLkvCVdqi6VwXcV5\nlrFE5u1uPGVnZwMAfHx87F2UTbBnELmzBqjzgHiOralocl8SrtQWS+HDaZxnGUtk3u7GU2ZmJgDX\nc3nbM4jc0QHqPCCeY2tcSe5dqS0cDsc8lsi83WOesrKyALh2r81Rq4TbO4jcmrx5sCenJFxJ7l2p\nLWXBUZNfuG7hlDeWyLzdjSelUgkPDw+4uTlkMfNyoSyCbum1jggiL03eHI4pXEnuXaktZcHRk184\nnPLCEpm3+7CdUqmEQqGwdzEcDseJcCW5d6W2lAW+HhTnWcESmXfICuOenp72LsauWDO0Zi94EDmn\nIuEKci/gSm0pC6UdtuNwKhqWyLzdjSeVSgWZTGbvYuxKWTbOtBU8iJxTkXAFuRdwpbaUhdJufM7h\nVDQskXm7D+JrNBq4u7vbuxinwxHBkTyInOOsuJLcu1JbbI0pPWJOt/Agco6zY4nM29140ul0z2Sg\npSMEnQeRc5wVV5J7V2qLrTGlG8zpCx5EznF2LJF5uw/bccXD4Tx7uJLcu1JbOByOeZzCeNJqtZBK\npfYuxuHYM4jc0QHqzhAQz3EtXEnuXaktHA7HPJbIvN2Np/z8fEgkEnsX43DsGUTu6AB1ZwiI57gW\nriT3rtQWDodjHktk3u6+aCKCWGx3G61C4KxB5KUpjwd6ckrCleTeldriDBTWUaaOcd3CKW8skXmH\nDOSLRCJHFOP03L37dLabvXBUMCZXWhxzuJLcu1JbyhtThhEPIuc4G+Zk3u7dKSKydxEcDsfJcCW5\nd6W2cDgc81gi89wXbUNKswKvrfPmQeQcDsdZ4HqE46o4xHh6VnpuzjBMxoPIOc6CK8m9K7XFkXA9\nwqmomJN5u8c8SSQSaLVaexdTYXFUELkjAjB5oCdHwJXk3pXaUpHgQeSc8sISmXeI8ZSfn2/vYios\njggid0QZAFdanKe4kty7UlucFVP6iQeRc8oLS2Te7sN2bm5u0Ol09i6Gw+E4Ea4k967UFg6HYx5L\nZN6s8XTmzBkolUr29+XLl5GUlGSQJj4+HsOGDUOVKlWwatUqg3PPeq/NGYLIHV0eD/7kOIvcP3jw\nANevX2d/Z2Rk4MKFCwZpVCoVlixZgvDwcAwaNAhpaWkG552lLa4M1xkcZ8ImnqcRI0Zg69atAIDs\n7Gx07twZR44cYeevXbuGVq1aYdeuXahUqRKWL1+O1157zapKuDLOEETu6PJ48CfHWeR+8+bNGDVq\nFPt71qxZmDdvHvubiNCvXz/MnTsXqampSE1NRfPmzfHo0SOWxlna4spwncFxJiyRebMxTyEhISyT\nzMxMZGZmIj09nZ3//PPP0aJFCyxduhT+/v4AgPbt22PVqlWYNm0ad3kXwhrPDV+JnFORcRa5DwkJ\nMahHYmIisrKy2N9Xr17FmTNn8Pvvv0Mmk6Fp06YYPXo0+vfvjwsXLkAkEjlNW1wFa/QgX4mcUx5Y\nIvNmjafCGYSHh2Ps2LE4efIkJk+eDJ1Oh/Xr12Pr1q1o2rQpSzdo0CCcP38eAODu7g6NRlPaNrgU\n1nhu+ErknIqMs8h9UQU4b948tG3blu2avnr1anTr1g0dO3ZkaaZPn46mTZtCr9dDIpE4TVtchbJ4\nsHkQOccRWCLzJQ7bZWVl4cKFC5DJZOyYTCZj6x8cO3YMIpEIgwYNYuf1ej3i4+PRsmVLAIBCoWAx\nU/Xr14dIJDL5mz9/fqkayeFwyo/58+cbyfLNmzcN5L48OXTokJH+Ap6u4bJv3z68+uqrBtdcu3YN\nzZo1Y3tbcR3G4bgmZdFfJRpPwnBdhw4d2LHU1FT8+uuv2LNnD/Lz8+Hm5mawB8y2bdsQHx+PESNG\nAAC8vb2RnZ0NoEAJcewLX4mcU95kZ2cbyH15kp+fb+BVSk1NBQC8/vrryMvLYzpM4OHDh5gxYwbe\nffddpte4DuNwnh0s1V8lGk/e3t6IjIzE2LFjMWLECISGhmL79u1IS0vDt99+i+DgYGRlZeH48ePQ\naDRYs2YN4uLisGPHDhb/pFAooFKpkJ+fzxWPA+ArkXPKG6VSaSD35UmjRo2wadMmjB07Fo0aNUKf\nPn0AAD/99BMeP36M4OBgHDhwADqdDhcuXECnTp0watQo9OvXj+Xx/9o78yA3qnPtP5JG+zKSRuPB\ndmwWsxQXCJRviHGKZQykMAGDqRsMcVzGONeGkOCwFSSpJCNMqEoFu1wGsxhyyY0z7DcGYrOUgRgo\nQzFsdgqyOIRgE0yYfUbbaJmRvj/m647k0ahb6u7Tp4/eX9XUTLWkPv22+nnn7XOePodyGEE0D2rz\nV03PU0tLC66++mr88pe/xKFDh7BkyRK8/vrrsNvt6O7uht1ux5o1a9DZ2YlSqYQFCxZg//79aG9v\nl/cRDocBTJrNd+/erVN44sDCRK7UxnTGTDKRE0rE4/Gqw1Xvv/8+gEndSzdSZvDtb38bd999N154\n4QWcfvrpOPLII7Fz505s374dRx99NDZs2IBly5bhvvvuQ2trKx577DFceOGFFfugHMYPZCIn9ERL\n/lI0jHd1deHWW2+F1+sFANx0001yAQUA999/P9asWYODBw9i6dKlFUN4ABAKhQBMdoWZmUR5hYWJ\nXKkNJWOm3lDSEh9edH/SSSehv79fzl/79u3Dzp074ff7AQAXXHABPvnkEzz77LNYunQpIpHIlH3w\nEgtBJnKCDWo0r2p5FinxAMDs2bNx7LHHVrw+f/58zJ8/v+Znx8bG1DRFEIQA8KT78vwViURwzDHH\nyMUTAESjUVx99dWKn+chFoIgjEeN5ute2+7mm2+u6/1Skkqn0/U2RUwDr/NAST1jPM1nRZgDr7o/\n8sgj8fHHH9f1GV5jEQktqyWwth0Q4qNG84YvDCwdRCaTMbqppoGF2LV0hfM0nxVhDiLpXqRYeEVL\nTqv2WSqICC2o0bzhCwMHg0EA4OKxZavAes06gtAbkXQvUixWgvIgYRZqNG948SQ92kuJRz1q75qM\nnDuJ9bxMNA+UWIike5FisRI85EGiOVGjecOLJ/ILGIeRcyexnpeJ5oESC5F0L1IsVqKzM67qfZQ7\nCL3hwvMkPfo7NDRkdFNNjVGm7sPbMNpEXs++yezJLyLpXqRYrIReJnKaB4qoFzWaN7x4CoVCsNls\nGB0dNbqppsYoU/d0nzGyPUpa1kck3YsUi5XQkgfKCySaB4qoFzWaN3zYzmazIRQKUeLRiJHmSTJm\nEnojku5FisXqUK4iWKBG84YXT8BkF9jIyAiLpoTFSPMkr8ZMMoJaG5F0L1IsVqaRnqNqeYRyC6GE\nkuaZFE9tbW3o7+9n0VTTQyZyghdE0r1IsVgZtSbycqrlEcothBJKmjfc8yQdxPDwMIumiDJEmYm8\nnn2T0ZMfRNK9SLFYGSOG7Si3ENVQ0jyT4ikQCODQoUMsmiLKEGUm8kb2TZiPSLoXKRYrY0QeoNxC\nVENJ80yG7QKBAFKpFIumCILgBJF0L1IsBEEoo6R5JsVTOBymOVIMoJ4FeHnaNw/tEcYjku5FikU0\n6Ak8wgiUNM+keGpvb0cymUShUGDRXNNQzwK8PO2bh/YI4xFJ9yLFIho07EYYgZLmmXmegMl1YqLR\nKIsmiWng1UTeSHtk9OQbkXQvUizNAKsVFCi3iIuS5pkUT62trQCA4eFhSjwmw7uJvNF2CP4QSfci\nxdIMsF5BgRAPJc0zGbaLxWLyQRAE0RyIpHuRYiEIQhklzTMpnqQViulpFeMhEznBCyLpXqRYmo1q\neYRM5oQSSppnUjz5fD4AQCaTYdFcU0MmcoIXRNK9SLE0G9XyCA3LEUooaZ6554ngDzKRE0Ygku5F\nioWoTnmOqraNcktzoaR5Zsuz1DoIwlwOHPj3kilGQSby5kMk3YsUC1GdaoURmcibFyXNM5skEwCt\nSk4QTYRIuhcpFoIglFHSPJPiyeVywe12I5FIsGiOOAwjzZFkIiemQyTdixQLQRDKKGmeSfEETJqv\nyGxpDjwMk5GJvDkRSfcixUIQhDK1NM/E8wRMPvaXTqdZNUdohEzkhB6IpHuRYiHUQSby5qaW5pkV\nT16vF2NjY6yaIzRCJnJCD0TSvUixEOqoVhiVb6PcIja1NM9s2I4SD0E0HyLpXqRYCIJQhoviyePx\nIJfLsWqOUIBM5AQLRNK9SLEQjWN0jzzBD7U0z6x4cjgcGB8fZ9UcoQAPw2RkIhcfkXQvUiwEQShT\nS/PMPE8OhwMTExOsmiPqoJ6eGyNN3dWMmWQitzYi6V6kWIip1JMHyUTeHNTSPLPiyWazsWqKqJN6\nem6M6rJWmt1XbyhpsUEk3YsUCzEVtTlBaa08yi3iUEvzzIbtbDYbSqUSq+YIguAAkXQvUiwEQShT\nS/PMep4IsVDbNS31VDViUG+kK7ze9qiLnSCIelAativfRrlFXJgVTxMTE3C5XKyaIwymXrE3khy0\ndIXzMBRJiKV7kWIhGketxYAKIutTS/PMhu3y+TwlHgth5FQGRPMgku5FioVQB+XB5qaW5pkVT9ls\nFh6Ph1VzhEbU3jUZOXcSzctkfUTSvUixEOqg3qPmppbmmRVPhUIBTqeTVXMEI4ycO4nmZbI+Iule\npFgIdXR2xs0+BMJEammemeeJEo/4sDKRN2rCJBM5e0TSvUixEOrQa+6n8m2UW6wDF8UTdXmLDwsT\nuRaDN5nI2SOS7kWKhVCHltULyERufXQbttu3bx+++93v4u9///uU1zKZDP70pz9N+9lcLge3211P\ncwRn8LAeHmEteNL98PAw1q9fj+7u7qqv//nPf0Y6nZ728zzFQhCE8dTSvOri6ZlnnsFXvvIV/O53\nv0N7e3vFa7/5zW9w7LHH4uSTT8ZFF11UNQHRkyrWx0gTOa93ZGRa1wYvuu/r68MZZ5yBrq6uKd3w\nf/3rX9HZ2YmTTjoJ8+bNwzvvvFN1H7zEQlgDyh3WR5en7X7wgx9gzpw5+Pjjj9Ha2ipv37VrF1at\nWoVcLocrr7wSxWIRZ5xxBjKZzJSDoLu25kAko7dIsZgBL7q/++678be//Q3PP/88rrjiCnl7sVjE\nueeei7179+KYY47B9773PSxcuBB79uyZsg9eYiGsAeUO61NL86o8Ty+99BI+/fRTrF+/HsFgsOK1\nrVu34pprrsH999+PbDYLm82G8847Dz/60Y+wefNmAJMTTWWzWfj9fo2hEFaDhUGSxUzk9eybzJ6T\n8KL7bDaL7u5uHHfccbjwwgsrXtuxYwdyuRx6e3tRKpXg9XrR2tqKSy65BAMDA7DbJ+8veYmFMB8t\nFoPDcxUtJswvSppXVTzdeOON8Hq9WLt2bcX2gYEB/P73v8c777wDm80Gr9cLADj99NMxPDwsvy+V\nSgEAAoFAQ0EQ1oWF0FnMRN7IvpsdXnS/efNmHDx4EPfee++U1371q19h9erVFabQM888E6OjoxXv\n4yUWwny05IHpchXlFv5Q0rzisF2xWMTQ0BA6Ozths9kwODgov/bee+8hFArhtNNOk7clk0m8+OKL\nuOSSS+Rt2WwWAODxePDUU0/BZrNN+YnH4w0FSBCEecTj8ap6XrZsWYXuzUTKWRdeeCH6+vpQKBTk\n13p6erBo0aKK9z/11FO4+OKL5V4ngHIYQYiIlvylWDz19PTgX//6F1544QV0dHQgFothyZIl2L9/\nP4DJ4kqiWCzi5ptvxsyZM3HZZZfJ26UKLhgMyn8TYlDPPCg87VsLZARVJhQKVejeLIrFIp555hkA\nwDHHHIOOjg6ccMIJeOyxxyreI7Fr1y5s2bIFGzdurNgP5TBCD8pzB+URflGTvxSH7bZs2QIA2LRp\nE0455RQ8/vjjePjhh3Ho0CE888wzSKVS+OlPf4oFCxbgwQcfRCqVwpNPPllx15ZMJgFMdn999tln\nmgMj+EHLPChm7lsLZARVxu/3V+jeLHbt2oWPPvoIZ555Jh588EE8/fTT2LJlC5YvX45TTjkFp556\nKn7xi1/Abrdj7969+PWvf43HH38cxx57bMV+KIcRelCeOyiP8Iua/KVYPM2aNQsAMHPmTJx33nk4\n++yz8fzzzyMajWLu3Lno7u7GLbfcgo0bN+KWW25BPB6vKJyAyi7vdevWYd26dZoCI6wNq1nCyURu\nPPF4fNrhqrfeeguAucN2M2fOBDB5J3niiSfixBNPxNDQEDZv3oxQKISHHnoI3/nOd3DRRRfh61//\nOnbv3o05c+ZM2Q/lMKIWes1TRyZytmjJX4rF07XXXouNGzdiy5Yt2Lt3L55++mmMjY1hw4YNAIAr\nrrgCl156KZLJ5JT5nyRGRkYAoGKKA6J50ctwqfZ9ZCI3Bx50f+qpp+LMM8/EG2+8gdtvvx1vvvkm\nXnnlFaxfvx5z584FALzyyiv45z//WbVokuAhFoJf9MoDZCLnByXNK3qe5s2bh+7ubvzlL3/B5s2b\nsXjxYrz++usVJnGPxzNt4QQAQ0NDAIBYLFbXwRMEYV140X13dzeOP/54xONxTExM4PHHH8ePf/zj\nivfUKpwAfmIhCIINSppXNUnm8uXL0dvbi0Qigc2bN+Pkk0+u6yCkCi4cDtf1OcKa8GD01mvfPMRi\nVXjR/ZFHHomenh5ks1m8/PLL+OY3v1n3PniJheCTavpvdOFzgg+UNK96YWCHwwGHw9HQQUgVXCQS\naejzhLXgweit1755iMWq8KR7m82maXZwnmIh+EPtwsCN7IcwByXNqy6etPD5558jEonQulDEtPA6\nE3kj7ZHRcxKRdC9SLAQbOjvj8kMoWiATuTkoaZ5J8TQwMICOjg4WTREW5cAB6JJoasHKjElJaxKR\ndC9SLAS/VMuBZCI3ByXNq14YWAuDg4PU3U0QTYZIuhcpFoIglFHSPJPiaXh4GG1tbSyaIjhFr3lQ\ntOybTORsEUn3IsVCmEd57qiWR6rljGbPI2ahpHkmxdPAwADdtTU5PAyTkYmcLSLpXqRYCPNQGoKr\nljOaPY+YhZLmmXiehoaGaH4Uom7IRG5tRNK9SLE0A4fr+nBNmaXHchO5FkM5mciNR0nzhhdPpVIJ\n6XTa1PWtCGtCJnLrIpLuRYqlWVBaQ87qeiQTubGo0bzhw3aZTAbA5EJ7BEE0ByLpXqRYCIJQRo3m\nDe95SqfTAMxdWZ3gE15M5HrcubFuj3dE0r1IsTQL5T3W1Xqv9R7qqqZ/vfLbdMffDHnELNRo3vDi\nSargvF6v0U0RFoOHYTIykRuDSLoXKRZikgMHgNdei6t4p5r3VNe/kfmtWfKIWajRPLOeJ+ryJmpR\nT88NCxM5CwOmyEZPkXQvUiwiMp0hvJZh3CyUpiqohtrjJxO5fqjRvOHFk7S4Hj3mS9Sinp4bo03k\nLNoAxE5aIulepFhEROnx/mqvq9e32vepoxGjt9peJjKR64cazRtuGB8dHQUAhEIho5siCIITRNK9\nSLEQBKGMGs0b3vMkHURra6vRTRFNhtquaekusxEDZyNd4fW2J2IXu0i6FykWYhJeNafFZM7TfFZW\nR43mqXgiLEu9Ym8kOWjpCudpKJI1IulepFiISfQ2jOuFlgJGxPmszEKN5g0ftiOzJdEoRk5lQBiL\nSLoXKRaCIJRRo3lmk2T6fD6jmyIEQ+1dk5ELZ9KinI0hku5FioUgCGXUaN7w4imRSMDj8cDpdBrd\nFNGkGDnnCc2n0hgi6V6kWAiCUEaN5g33PCUSCXpKhWAGmcj5QCTdixQLIQZa5oiyem5hgRrNG148\n5fN5uN1uo5shCABkIucFkXQvUizEJGbN86QX1fJItZjIRN4YajRv+LBdLpejxEPoBgsTORnVtSOS\n7kWKRUSU1pWzup6rHT95MY1FjeaZLM9CRktCL+oxkTd6h6X0OdbLPfC0vIRaRNK9SLGIiNK6ctVe\nt9LQldp18xpZxoWojhrNG148ZbNZeDweo5shiApEMpFb0bQuku5FioWYhNd5nqrR2RmfMiRXbZva\nPGHFfMIaNZo3vHgqFAr0lAphKqwWEzbaRF7Pvs2+YxZJ9yLFwgKrXKMELSY8HWo0b3jxVCwWYbcb\nbq0iiGlhIXQWJvJG9m0WIulepFhYYKVenWakvMeKFhOujhrNM8kINpuNRTMEQXCESLoXKRaCIJRR\n0rzhxZPdbsfExITRzRBNTD1znvC0by3w/rSNSLoXKRYWvPpqHKUSFH9YTcXBu1aUUHqasF7Kz4fV\nz41RqNG84cWTw+GgxEMYSj1zJ/G0by3wbvoUSfcixdKM8K4VJdQ+baeW8vNh9XNjFGo0b7jnyel0\nYnx83OhmCKIuyERuLCLpXqRY9IZmsOYbVisoiPZdqtG84cWTy+VCPp83uhmCqAsykRuLSLoXKRa9\nqfaPlQzj/MB6BQVRUKN5w4ft3G43stms0c0QBMERIulepFiI5qKzM272IVgSNZo3vOfJ4/FQ4iGY\no2XhTDP3zUN7eiCS7kWKRW+qmb5Zrxs33dBhLS1Makq5fSst7aK3sbyZUaN5GrYjhEQkEzmvpvVa\niKR7kWIRkemGDmvBy02GnuhtLG9m1Gje8OLJ6/VibGzM6GYIQjO8msgbac9so6dIuhcpFi2o7eEx\n+9pTQ2dnXJUv65xzpi6Dwitql3FRQul7tsL3qxU1mqfiiSD+P7ybyBttxwxE0r1IsWhBbQ8PGcat\nTSMLLYuGGs0bbhgPhULI5/PI5XJGN0UQBCeIpHuRYiEIQhk1mje85ykSiQAAhoeHccQRRxjdHEFM\nC5nI2SGS7kWKRQtqzeF6GcYbMYJXY3ojtfJ+rGS4VmsYL99mpfhYokbzhhdPoVAIAJBIJJo68RDm\nQyZydoike5FisRKNGMGr0SxGarVxNtsQXCOo0bzq4imVSuHAgQOYM2cOWltbp7z+9ttvY8+eGEOU\nXgAAFJdJREFUPVixYgVmzJghb5feOzo6WtfBEwQPkIm8MXjTfalUwv79++HxeHBUldvtgYEBdHd3\n46tf/Sq+9rWvVbzGWyx6U48R3IqG8emM1M1qGC/fVu31Rr573r5zrajRvKri6dlnn8XatWvR19eH\ncDiMRx55BN/4xjcAAIVCAatXr0Z3dzcA4K677sL27duxcOFCAP/u/hoZGWk8EoIwiQMHjF/AVEQT\nOU+6/8c//oFrr70WL730Emw2G6655hrce++9sNsnLZ9PPPEE1qxZg2QyCQC48847ceutt6KlZTI9\n8hSLEdRjBFe7jQzj1qWR794KBVE9qNG8omH8zjvvxNKlSxGNRvHoo49i7ty5uPzyy3Ho0CEAwCOP\nPIKdO3fiyiuvxJYtW3DDDTfg4osvxv79+wEAPp8PAJBOpzUHRBCENeBF93v37sXpp5+O3bt3Y+PG\njbjuuuvwwAMP4N577wUAZLNZrFixAmeffTZWrlyJt956C5s2bcI999wj74OXWAiCYIMazdfseXr/\n/ffxk5/8BJdeeikefvhhRKNRtLS0YNmyZdi0aRM2bNiAe+65BzfeeCN+9rOfyZ/74IMPcMcdd6C7\nuxt+v1/xIAjCLIw0TDaziZwX3a9YsQIA0NPTg/nz5yObzaK7uxt33HEHrr/+ejzwwAM46qijsHPn\nTvkz3d3dWLx4Ma6//nq0tLRwE4uEaEMkWmnEFN3MhvHynnQtDwFU28bTKgdaUKP5msVTX18fAKBY\nLOL9998HAGzduhV2ux1nnXUWPvroI3z44Yd47rnnKj4XjUblv4PBIIBJz9Rjjz2G5cuXV22rq6sL\n8XhcISSC0BcehsmsbCKPx+O4/fbbK7YtX74cd911F4BJ3ZtJX18fSqUSPv/8cwwPD+Ptt9/G6Ogo\nli1bBmByyO773/9+xWfa2toQiUTkYT3echgNi1XSiClahH/wh8NDnDytcqAGLfmrZvG0ePFi3H77\n7ejq6sKOHTvk7QsWLMCll16Kl156CT6fr8KN/umnn+K3v/2tXFB5vV4AwNjYGDo6OuoMjSD4g0zk\ntclkMhW6N5OXX34ZZ511FpYsWVKx/ec//zkA4JNPPsG8efMqXlu/fj2WL18uF0+Uw/imEVN0MxvG\nD88tej4scPjrVuwhVZu/FA3jp512GgBg1qxZWLZsGXbs2IGenh50d3ejo6NDnkjK7XZjeHgYl112\nGVauXCk/sRIIBABMPvJX7Sk9grAaZCKvTSqVqtC9mRx//PHw+XzIZrNYtmwZRkZG8Nxzz+HGG2/E\n9u3bAUA2igOTd6J//OMf5QdgAMphhFiUFzTViptGpjyQKM+LPBVE9aA2fykWT1u3bgUAbNiwAd/6\n1rewevVqfPnLX8YDDzyAXbt2IRQK4fzzz8f8+fOxfft2XHvttfjhD38of97pdCIYDGJkZASLFi1C\nqVTSGhtBEJwQj8enHaqSdG8mzz77LHp7e7F48WK5IDrttNPw4osv4sCBA7j88suxbt069PT0YN++\nffB4PHjttdfkeV4AymEEISpa8pdi8XTFFVfg+eefx5NPPolSqYRNmzbhuOOOwxNPPAGfz4c9e/ag\nq6sL7777LrZt24ZFixZN2Ud7ezu++OIL9RERhImQiVwfeND9+eefj2g0in379mHbtm3Yu3cvPvzw\nQzz66KM4/vjjsXHjRnzpS1/Ctm3bcNVVV+GGG26Ay+Wash8eYpHQawZv3mlkxuxq28gwXh96DdtZ\nHSXNKxZPK1euRC6Xw6ZNm9DV1YUZM2Zg9+7dmD17NgBg3rx5FV3c1YjFYhgcHKzz0AnCHHgYJrOy\niVyCB93HYjH84Q9/wM0334w77rgDNpsNjz76qGwYd7lcuO2223Dbbbcp7sfsWJoNLcNHZBiffpsS\nWub4EgklzauaJHPNmjVYs2ZNwwcRCARMf+qGIBqhnp4bI03dLGdx1mvfvOj+1FNPxcsvv6xpH2bH\nYuQDBHr18FRua6xdwtqINBO5kuZVFU9a6ejoQE9PD4umCEJX6um5McpErtbUqRd67Vsk3Zsdi5Gz\nOevVw1NrG9EciLSWnpLmFWcY1+sgpDmjCIJoDkTSvUixGE1nZxw2GxR/OjvjZh8qQUyLkuaZ9DyF\nw2GkUilMTEzA4XCwaJIgTENt17TUU6XV1GlUe1q72EXSvVGxNNP8OYT4WHXx6GooaZ5J8VQ+Qy/N\nk0KITr1ibyQ5aOkKZzUUKZLujYqlESMuzTBO8IpIQ7tKmmcybCct19Lf38+iOYIwFDK6qkMk3YsU\nC9HcKE3jQEyipHkmxZO0fEtvby+L5gjCUNTeNa1ape59jWDkvvVqTyTd1xtLtfPF+jszi1dfjaNU\nguKPVZY9EQ2lXmutCwOLgpLmmRRPsVgMADAwMMCiOYLgAiPnQWE9x0oj7Ymk+3pjacZ5cQhCpGtc\nSfNMPE/SeKHZ61wRhFmIYiJXu++LLpqPpUtPBCCG7imHqccKZmBievSeYVzN3E88zgOlpHkmxZO0\nTlT5ApwE0UyIYiKv571Sd7cIuqccph4qiKyN3qbveuZ+4unaUdI8k2E7ybVOd22EqPCwHh5viKR7\ns2Ox6jVAEFZFSfNMiiev1wuHwyFEEiWIahhpIufBoN5IeyLp3uxYeLojJ4hmQEnzTIonm81m+tpQ\nBMEDzWQiF0n3ZsdCs3ETBFuUNM/E8wQAbrcbuVyOVXMEwT0szJFa2tBiIp87N4Zt274vlO6NiEXt\nYrtTX483tG+1aDEN82T6JerHTMP4dK+bdT3V0jyz4snn8yGTybBqjiC4h8U/Di2zhNdzfIfPej13\n7uRvkXRvRCyNGHFZXTdattEM6NZFy3dfjUYWC+bFRF5L80yG7QDA7/cjnU6zao4gCBORpC6S7o2I\npdpwnNIQHS28SxBsqKV5Zj1PLpcL+XyeVXMEYTr1LMBb790V66ev6m1PkrpIuq8nFrXDcY1sYzFs\nR4iL0pCZFkS75mppnlnx5HQ6USgUWDVHEKZTzwK8Ru1bL+ptT5K6SLqvJxYtc+XwPIxBWJ/yfKP3\nQyaiXZu1NM+seGppacH4+Dir5gjCkrAykRttxpSkLpLuRYqFIIygEbN5o59hYSKvpXlmxZPdbkex\nWGTVHEFYEhZ3bix6MSSpi6R7kWIhCCOoVtBoNYyr/YwR1NI8M8O43W5HqVRi1RxBECYiSV0k3YsU\nC0EQytTSPLOeJ5vNxqopguAWHkzkeplEa7Un5RuRdM9LLDSPEkGwoZbmmRVPpVKJm+RDEGbBg4lc\nL5NorfYkqYuke15ioXmUCIINtTTPrHiamJiAy+Vi1RxBCAOLWcK1mjEPf3ze4bAD+JlQuhcpFiW0\nT6sQb6gNkTDCKK3HNiW0fPfVYHUejOiRran5EiPOOOOM0vnnn8+qOYKwDOec01VzW7XX9W7PKETS\nPS+xnHNOV2lyYLT2D8vvmZiKkq713o/afYuSW1hQS/PMDOPj4+NoaWHW0UUQBAeIpHuRYiEIQpla\nmmdWPOXz+abp8iaIeuBh+GLVqrgh+xVJ9yLFQphPNc3prcNq+5s6W73yZ1gcK4/U0jyz4qlQKMDp\ndLJqjiAsAw9PROk907CESLoXKRbCfLQstqulDaX5kvReGNjK1NI8sz7oXC4Ht9vNqjmCEBqrzEQu\nku5ZxaKX8ZcgqtHZGcerr8an/K0H0/VkmWWc15ona2meWfGUzWbh8XhYNUcQQnPgAHRNetXQY0Zf\nkXTPKpZmvcsnqmvaSJ2rbU/tNr1mC+dlHchammc2bJfJZODz+Vg1RxAEB4ike5FiIQgj6OyMm30I\nulJL88yKp1QqhUAgwKo5grA0RprI65mfRSsi6d7MWMrv8o86CjjnnLjiDw8PIhD1wdqYrcUcrmRG\nF4FammcybJfP5zE+Pi5MEiUIozHS38JqJnKRdM9TLOR9EhfWQ7ZazOFqFwG2KkqaZ1I8DQ8PAwDC\n4TCL5ghCKOrpKTLKRD6dMXO69s4992Rcc81ZAMTQvdk5zOiZlAn90TIrt5H7NrJ3qFFTN8uZ1tXo\nR03+YlI89fX1AQDa29tZNEcQQlFPT5FR5tJG7jI/+OADAGLo3uwcVn7+aW07a6DW4KyEFmN2NYws\nqNXmCVaGcbXtVUMpfzHxPI2OjgIAIpEIi+YIguAAkXQvUiyEmKg1a4tm6jYKJc0z6XnKZrMAIMx8\nLwTBO2qHdqS72ka68mt1ha9evQizZhUAiKF7ymGE3ui1sG6t/Rk5dKgW3obtyrdNlyPV5C8mxVMi\nkQAABINBvPjii1i8eDGLZg0jHo9X/duKUCz8oUcc9XbNN9KVr9QVvn37dgCTurc6lMP4xMpxHK6Z\neDyOzs7J3+WxaDFr6zV0WC/lx8/bsF2tbeUo5i8WKxM/9NBDJQClgwcPllauXMmiSUMBIP9YHYqF\nP9TEcdVVXVO2GblCeiP7Lte91WGZw5S+23PO6SoBJcWfWt9ZM2nFKkwXS7XvUa9tRu1bhO9FKX8x\n8Tz19vYCmDReDQwMsGiSIIRG7d3jqlVxdHYq/zQyl0z5Zw5v5733/lGhe6tjVA6rdt7pCTlCJJTm\niNJ7rqlGtjWSv5gM2w0ODsLv98Pr9coHRBCE8Rj5ZJbSE2Dlurc6RuWwasMveq83RhBmojTsqPdc\nU41sayR/Mel5SqfT8Pv98t8EQfCPFmNpIFCpe6tDOYzgnUb0qsVkLtps4uWoyV/MDOOhUEj+2+qs\nXbsWfr9fmH8MotDf34+WlhY4HA6zD4VrvF4gEgF8PmDGDGDmTGDhwv+c8r7y4aPu7uunvP7KKz+T\n/5Z6SopFYGICaGkBNm/ejDfeeEP34zcDymGE0Xz22Wdwu91oaan8t1yuMwm1cz9t3Lhyyrb/+q8F\nNT/Dk/Fcidmzo1O2XX31oinb/ud/vjtl265dP5H/Lj+fxSJgswFbtmzBq6++Om3bTIqn8gruvffe\nY9Gk7hQKBaRSKYyNjeGmm25COp1GJpPBjh07MDo6ikQigeHhYSQSCYyNjWFsbAz5fB7ZbBa5XA75\nfB6FQgETExMoFosAALvdDqfTCY/HA7/fD7fbDZfLhVAohFAoBJ/PB7/fj9bWVvknHA7D7/fLr/l8\nPthstobjOvfcc5FKpZDL5fQ6VZoplUry+Uyn00ilUkgmk0in00in0xgbG0MikcDIyIj8WiqVkqfT\nn5iYwMTERMXfEtK5stlscDgccDqdaGlpQUtLC1wuF3w+H7xeLwKBAAKBAEKhEFpbWxEMBhGLxRAO\nhxEOh9He3g6v16vp3E+Hw+FAsVhEqVSq63OXXXbClG133vktFApAMAjEYkB7+2TRJFEqlTA0NITB\nwRPw7rvvIpFIYHR0FCMjIxgYGMDIyAiSySQymYx8PUvX9vj4uHwtS8dqt9vhcDhgt9vR0tICp9PZ\n+IngCKNymMMxtfO/2raFC2Py35dcchzOPff/4PW6EAh40NrqQ1tbEMGgBz6fG36/G21tQbS2Vv7T\nFTGHLVmyRM4RPMFTDtuzpzKHnX32LBw8eHDaHPbJJ304+ugZFdv+939fxapVnRXbrrvuV7jvvv+u\n2PYf/+EBUJnDqn23bW3/foJt3ryOKa+XF3gSW7eunbLtkUfWleWwQSQSCcyZM4Gnn35alxzm8Xim\ntClhK9WboRvgggsuQDKZxJtvvokVK1bIdzxtbW3o6OhANBpFKBSCx+NBa2srZsyYgUgkApfLpetx\nlEoljI2NIZlMyv8kBgYG8MUXX2BoaAjJZBIjIyPo7e1Fb28vUqkUBgYG0N/fj2QyqaoNl8sFr9cL\nr9cLl8sFj8cDj8cDl8sl94o4HA6USiUUi0UUCgVks1lkMhnkcjnkcjkkEgl5Xhk1+P1+RKNRxGIx\nRCIRtLW1IRaLYdasWYhGowiHw/Jvj8eDYDCIQCAAj8eDQCCge09NPp/HyMgIBgcHMTo6ikwmg3Q6\njYGBAYyOjiKVSqG/v19OHKlUCoODgxgcHJSTt9rL0ufzyfEcfo7L/y4XsHTuJyYmUCgUMD4+jvHx\nceTzeWQyGWSzWVUFpdfrRSwWQywWQ3t7u1xUBQIBRCIR+TqORCIIhUJyMab3tV0qleRrJ5lMYnBw\nEH19ffj888/lhCKd476+PgwMDMhJfXBwsGacDocDgUAAPp9PvrbLr2e73Q67/d//7KVEXywWUSwW\n4XQ68eabb+oWq1lQDqMcRjmsuXJYqVSC3W6fNn8xKZ4mJibQ2dkJAPjiiy+QSCTkCrwWLS0tctAe\njwdutxtOp3PKBWaz2eSLqVgsYnx8HIVCAfl8Xr67yuVySKVSFRV8NQKBAGbMmIEjjjgCgUAAbW1t\nOOKIIxAOhxEKheSKXkqe4XAYra2t8gWl1yR6hUJBFuzo6GjFnaF055JKpeS7GuliksTb29srr8el\nhHQhuVwuuN1ueDweOJ1O+Y7S5XLJyalUKqFUKskXmFTFj42NVZxrJYLBIKLRqHx3FIlEZPFKd6jR\naBR+vx/BYBDBYFC+U/V6vQgGg2htbTV0iG58fFz+B5VIJNDf34/R0VH5PPf398tiHhwclO9+kskk\nCoVCzX1LIg4EAvB6vXA4HPI2p9MJu90uJ0spSUrXtfT95/N55PN5pFIpjI+P12zP5/PJyTAWi8nn\ntq2tDbNnz0YsFkMwGJR72sLhMGKxGAKBgKbetYmJCSGGUSmH1Q/lMMphVs9htfIXk+JpOjKZDAYG\nBuQvK5fLYXh4eEpFL1XSuVxOTihSd6bU5QZMdmOWDxe4XC65W9jtdstVcygUkk9yLBZDR0cH2tvb\n4ff7p4w3W5lcLoeRkRG5+3J0dFS+I0kmk/Lf6XRa7pbP5XLIZrMoFAryHWU+n5fPs3QRSndFbrdb\nvkuVznUoFEI4HJbvWHw+n5xcIpGIfIclMtLd09DQkHwnKp338qERKYlId41jY2MoFApTrmvpfDud\nTjn5Sv8spOtausaj0ShmzJiBWbNmCXld8wTlMGOhHGYelMNqY2rxRBAEQRAEYTWYTFVAEARBEAQh\nClQ8EQRBEARB1AEVTwRBEARBEHVAxRNBEARBEEQd/D/ZeXoYkWWo5wAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x104336ad0>"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "<5x5 sparse matrix of type '<type 'numpy.float64'>'\n",
        "\twith 13 stored elements in Compressed Sparse Column format>"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## QR-decomposition\n",
      "The next decomposition: **QR** decomposition. Again from the name it is clear that a matrix is represented as a product  \n",
      "$$\n",
      "    A = Q R, \n",
      "$$\n",
      "where $Q$ is an **orthogonal (unitary)** matrix and $R$ is **upper triangular**.  \n",
      "The matrix sizes: $Q$ is $n \\times m$, $R$ is $m \\times m$."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "This algorithm plays a crucial role in many problems:\n",
      "- Computing orthogonal bases ior a linear space\n",
      "- Used in the preprocessing step for the SVD\n",
      "- QR-algorithm for the computation of eigenvectors and eigenvalues (one of the 10 most important algorithms of the 20th century) is    based on the QR-decomposition"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Theorem \n",
      "Every rectangular $n \\times m$ matrix, $n \\geq m$ matrix has a QR-decomposition. \n",
      "There are several ways to prove it and compute it.\n",
      "\n",
      "- (mathematical) Using the Gram matrices and Cholesky factorization\n",
      "- (geometrical) Using the Gram-Schmidt orthogonalization\n",
      "- (practical) Using Householder/Givens transformations \n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Mathematical way\n",
      "$$A = QR,$$\n",
      "then $A^* A = R^* R$, the matrix $A^* A$ is called **Gram matrix**, and its elements are scalar products of the columns of $A$.  "
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "It is simple to show that $A^* A$ is positive definite:\n",
      "$$\n",
      "   (A^* A y, y) = (Ay, Ay) = \\Vert Ay \\Vert^2  \\geq 0.\n",
      "$$\n",
      "Therefore, $A^* A = R^* R$ always exists.\n",
      "Then the matrix $A R^{-1}$ is orthogonal:  \n",
      "$$\n",
      "   (A R^{-1})^* (AR^{-1})= R^{-*} A^* A R^{-1} = R^{-*} R^* R R^{-1} = I.\n",
      "$$\n",
      "\n",
      "**Gram matrices are not good for numerical stability**!"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Let us show that on random and Hilbert matrices :)\n",
      "The Gram-matrix method may lead to the loss of orthogonality."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "n = 10\n",
      "r = 5\n",
      "#a = np.random.randn(n, r)\n",
      "a = [[1.0/(i+j+0.5) for i in xrange(r)] for j in xrange(n)]\n",
      "a = np.array(a)\n",
      "q, Rmat = np.linalg.qr(a)\n",
      "e = np.eye(r)\n",
      "print 'Built-in QR orth', np.linalg.norm(np.dot(q.T, q) - e)\n",
      "gram_matrix = a.T.dot(a)\n",
      "Rmat1 = np.linalg.cholesky(gram_matrix)\n",
      "q1 = np.dot(a, np.linalg.inv(Rmat1.T))\n",
      "print 'Via Gram matrix:', np.linalg.norm(np.dot(q1.T, q1) - e)"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Built-in QR orth 8.10893700925e-16\n",
        "Via Gram matrix: 1.78151575967e-08\n"
       ]
      }
     ],
     "prompt_number": 50
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Second way: Gram-Schmidt orthogonalization\n",
      "QR-decomposition is a mathematical way of writing down the Gram-Schmidt orthogonalization process.  \n",
      "Given a sequence of vectors $a_1, \\ldots, a_m$ we want to find orthogonal basis $q_1, \\ldots, q_m$ such that every $a_i$ is a linear combination of such vectors.  \n",
      "\n",
      "**Gram-Schmidt:**\n",
      "1. $q_1 := a_1/\\Vert a_1 \\Vert$\n",
      "2. $q_2 := a_2 - (a_2, q_1) q_1, \\quad q_2 := q_2/\\Vert q_2 \\Vert$\n",
      "3. $q_3 := a_3 - (a_3, q_1) q_1 - (a_3, q_2) q_2), \\quad q_2 := q_3/\\Vert q_3 \\Vert$\n",
      "4. And go on\n",
      "Note that the transformation from $Q$ to $A$ has triangular structure, since from the $k$-th vector we subtract only the previous ones."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## Gram-Schmidt is unstable\n",
      "Gram-Schmidt can be very unstable (i.e., the produced vectors will be not orthogonal, especially if $q_k$ has small norm).  \n",
      "This is called **loss of orthogonality**.  \n",
      "There is a remedy, called **modified Gram-Schmidt** method. \n",
      "Instead of doing \n",
      "$$q_k := a_k - (a_k, q_1) q_1 - \\ldots - (a_k, q_{k-1}) q_{k-1}$$\n",
      "we do it step-by-step. First we set $q_k := a_k$ and orthogonalize sequentially:\n",
      "$$\n",
      "   q_k := q_k - (q_k, q_1)q_1, \\quad q_k := q_{k} - (q_k,q_2)q_2, \\ldots\n",
      "$$\n",
      "In exact arithmetic, it is the same. In floating point it is absolutely different!\n",
      "\n",
      "Note that the complexity in $\\mathcal{O}(nm^2)$ operations"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## QR-decomposition: the (almost) practical way\n",
      "\n",
      "If $A = QR$, then  \n",
      "$$\n",
      "R = Q^* A,\n",
      "$$\n",
      "and we need to find a certain orthogonal matrix $Q$ that brings a matrix into orthogonal form.  \n",
      "For simplicity, we will look for an $n \\times n$ matrix such that\n",
      "$$\n",
      "   Q^* A = \\begin{bmatrix}\n",
      "   * & * & *  \\\\\n",
      "   0 & * & * \\\\\n",
      "   0 & 0 & * \\\\\n",
      "   &0_{(n-m) \\times m}\n",
      "   \\end{bmatrix}\n",
      "$$\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "We will do it column-by-column.  \n",
      "First, we find a Householder matrix $H_1 = (I - 2 uu^{\\top})$ such that (we illustrate on a $4 \\times 3$ matrix)\n",
      "\n",
      "$$\n",
      "   H_1 A = \\begin{bmatrix}\n",
      "    * & * & * \\\\\n",
      "    0 & * & * \\\\\n",
      "    0 & * & * \\\\\n",
      "    0 & * & *\n",
      "   \\end{bmatrix}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Then, \n",
      "$$\n",
      "   H_2 H_1 A = \\begin{bmatrix}\n",
      "    * & * & * \\\\\n",
      "    0 & * & * \\\\\n",
      "    0 & 0 & * \\\\\n",
      "    0 & 0 & *\n",
      "   \\end{bmatrix},\n",
      "$$\n",
      "where\n",
      "$$\n",
      "  H_2 = \\begin{bmatrix}\n",
      "  1 & 0 \\\\\n",
      "  0 & H'_2, \n",
      "  \\end{bmatrix}\n",
      "$$\n",
      "and $H'_2$ is a $3 \\times 3$ Householder matrix."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "source": [
      "Finally, \n",
      "$$\n",
      "   H_3 H_2 H_1 A = \\begin{bmatrix}\n",
      "    * & * & * \\\\\n",
      "    0 & * & * \\\\\n",
      "    0 & 0 & * \\\\\n",
      "    0 & 0 & 0\n",
      "   \\end{bmatrix},\n",
      "$$\n",
      "\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "A simple implementation of the Householder method in Python looks as follows."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#Initialization\n",
      "import numpy as np\n",
      "n = 4\n",
      "r = 3\n",
      "a = np.random.randn(n, r)\n",
      "k = 0"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def householder(x):\n",
      "    u = x.copy()\n",
      "    u = u / np.linalg.norm(u)\n",
      "    u[0] = u[0] - 1 #Or plus\n",
      "    u = u / np.linalg.norm(u)\n",
      "    u = np.reshape(u, (-1, 1))\n",
      "    return u\n",
      "\n",
      "\n",
      "x = a[:, k:k+1]\n",
      "u0 = householder(a[k:, k:k+1])\n",
      "u = np.zeros((n, 1))\n",
      "u[k:] = u0\n",
      "H = np.eye(n) - 2 * u.dot(u.T)\n",
      "a = H.T.dot(a)\n",
      "print a\n",
      "k = k + 1\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "fragment"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  2.55804239e+00  -3.02487125e-01   1.22104370e-01]\n",
        " [  2.18724250e-16   1.12507545e+00  -3.09474748e-01]\n",
        " [ -2.64480936e-17  -3.45664708e-16   4.80083115e-01]\n",
        " [  2.76406866e-17   1.07140777e-16  -5.55111512e-17]]\n"
       ]
      }
     ],
     "prompt_number": 31
    },
    {
     "cell_type": "markdown",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "source": [
      "## The QR algorithm\n",
      "The QR algorithm (Kublanovskaya and Francis independently proposed it in 1961). Do not **mix** QR algorithm and QR decomposition. \n",
      "QR-decomposition is the representation of a matrix, whereas QR algorithm uses QR decomposition computes the eigenvalues!\n",
      "\n",
      "Computation of the eigenvectors and eigenvalues of a matrix.  \n",
      "\n",
      "\n",
      "- Initialization: $$A_0 := A$$\n",
      "- Iteration step: $$Q_k R_k = A_k, \\quad A_{k+1} = R_k Q_k,$$\n",
      "so we compute the QR-factorization of a matrix, and then multiply them in the opposite order.  \n",
      "Good news: the spectrum of $A_k$ is the same as of $A_0$, and $A_k$ converges to a triangular matrix!  \n",
      "It is too tempting not to implement it :)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "n = 4\n",
      "a = [[1.0/(i + j + 0.5) for i in xrange(n)] for j in xrange(n)]\n",
      "niters = 10\n",
      "for k in xrange(niters):\n",
      "    q, rmat = np.linalg.qr(a)\n",
      "    a = rmat.dot(q)\n",
      "    print 'current a:'\n",
      "    print a\n"
     ],
     "language": "python",
     "metadata": {
      "slideshow": {
       "slide_type": "slide"
      }
     },
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "current a:\n",
        "[[  2.40047183e+00   1.43485636e-01   4.99605047e-03  -5.56291523e-05]\n",
        " [  1.43485636e-01   3.59286592e-01   1.60534687e-02  -1.94225770e-04]\n",
        " [  4.99605047e-03   1.60534687e-02   1.60692682e-02  -2.80885670e-04]\n",
        " [ -5.56291523e-05  -1.94225770e-04  -2.80885670e-04   2.40684296e-04]]\n",
        "current a:\n",
        "[[  2.41031150e+00   2.09437993e-02   3.18722176e-05   5.45279213e-09]\n",
        " [  2.09437993e-02   3.50196113e-01   6.87806955e-04   1.28079871e-07]\n",
        " [  3.18722176e-05   6.87806955e-04   1.53250849e-02   4.17732654e-06]\n",
        " [  5.45279205e-09   1.28079871e-07   4.17732654e-06   2.35678648e-04]]\n",
        "current a:\n",
        "[[  2.41051991e+00   3.04114601e-03   2.02621961e-07  -5.33227597e-13]\n",
        " [  3.04114601e-03   3.49989111e-01   3.00995528e-05  -8.62074652e-11]\n",
        " [  2.02621961e-07   3.00995528e-05   1.53236760e-02  -6.42429491e-08]\n",
        " [ -5.33147641e-13  -8.62074827e-11  -6.42429491e-08   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052431e+00   4.41545673e-04   1.28806658e-09   1.32059800e-16]\n",
        " [  4.41545673e-04   3.49984720e-01   1.31786000e-06   5.80333420e-14]\n",
        " [  1.28806664e-09   1.31786000e-06   1.53236733e-02   9.88053919e-10]\n",
        " [  5.21260158e-17   5.80510077e-14   9.88053906e-10   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   6.41081268e-05   8.18816607e-12  -7.99356474e-17]\n",
        " [  6.41081268e-05   3.49984627e-01   5.77009674e-08  -2.14109052e-17]\n",
        " [  8.18822358e-12   5.77009674e-08   1.53236733e-02  -1.51962427e-11]\n",
        " [ -5.09637195e-21  -3.90911829e-17  -1.51962302e-11   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   9.30787458e-06   5.19949264e-14   7.99300814e-17]\n",
        " [  9.30787458e-06   3.49984626e-01   2.52637036e-09  -1.76560777e-17]\n",
        " [  5.20524342e-14   2.52637031e-09   1.53236733e-02   2.33729939e-13]\n",
        " [  4.98273388e-25   2.63237618e-20   2.33717423e-13   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   1.35141258e-06   2.73389069e-16  -7.99300126e-17]\n",
        " [  1.35141258e-06   3.49984625e-01   1.10614263e-10   1.76826923e-17]\n",
        " [  3.30896669e-16   1.10614211e-10   1.53236733e-02  -3.60708065e-15]\n",
        " [ -4.87162969e-29  -1.77262591e-23  -3.59456477e-15   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   1.96211922e-07  -5.54040645e-17   7.99300027e-17]\n",
        " [  1.96211922e-07   3.49984625e-01   4.84316731e-12  -1.76827548e-17]\n",
        " [  2.10350596e-18   4.84311567e-12   1.53236733e-02   6.78001457e-17]\n",
        " [  4.76300288e-33   1.19367537e-26   5.52842647e-17   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   2.84880569e-08  -5.74941943e-17  -7.99300013e-17]\n",
        " [  2.84880568e-08   3.49984625e-01   2.12101875e-13   1.76827614e-17]\n",
        " [  1.33719609e-20   2.12050235e-13   1.53236733e-02  -1.33661508e-17]\n",
        " [ -4.65679822e-37  -8.03813647e-30  -8.50269816e-19   2.35677492e-04]]\n",
        "current a:\n",
        "[[  2.41052440e+00   4.13618795e-09  -5.75074807e-17   7.99300011e-17]\n",
        " [  4.13618792e-09   3.49984625e-01   9.33601458e-15  -1.76827623e-17]\n",
        " [  8.50053870e-23   9.28437503e-15   1.53236733e-02   1.25289581e-17]\n",
        " [  4.55296169e-41   5.41283161e-33   1.30771163e-20   2.35677492e-04]]\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Take home message\n",
      "- LU decomposition, Cholesky, their complexity is $\\mathcal{O}(n^3)$, used to solve dense linear systems\n",
      "- QR decomposition, find orthogonal basis, the basic component for the computation of eigenvalues and eigenvectors"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Next time\n",
      "- How do we compute LU factorizations for large-scale problems? \n",
      "- And: the SVD (with some examples)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "n = 100\n",
      "a = np.random.randn(n, n)\n",
      "%timeit np.linalg.qr(a)\n",
      "%timeit scipy.linalg.lu(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 777 \u00b5s per loop\n",
        "1000 loops, best of 3: 218 \u00b5s per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit l, u = np.linalg.cholesky?"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 35
    },
    {
     "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: 110%;\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": 137,
       "text": [
        "<IPython.core.display.HTML at 0x1224ad390>"
       ]
      }
     ],
     "prompt_number": 137
    }
   ],
   "metadata": {}
  }
 ]
}