lectures/lecture8/lecture-8.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 8: Symmetric eigenvalue problem and SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the first part\n",
    "\n",
    "- QR decomposition and Gram-Schmidt algorithm\n",
    "- Schur decomposition and QR-algorithm (basic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Plan for the second part\n",
    "\n",
    "Today we will talk about:\n",
    "\n",
    "- Algorithms for the symmetric eigenvalue problems\n",
    "    - QR algorithm (in more details)\n",
    "    - Divide-and-Conquer\n",
    "    - bisection\n",
    "- Algorithms for SVD computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Matrices with complex spectrum\n",
    "\n",
    "- What matrices have complex spectrum?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### QR algorithm for matrices with partially complex spectrum\n",
    "\n",
    "- QR algorithm computes Schur form of the given matrix\n",
    "\n",
    "$$ A = UTU^* $$\n",
    "\n",
    "- Can we get complex eigenvalues with QR algorithm?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  0 -1]\n",
      " [ 1  2  0]\n",
      " [ 3 -1  5]]\n",
      "[4.20556943+0.j         1.89721528+0.66545695j 1.89721528-0.66545695j]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "A = np.array([[1, 0, -1],\n",
    "             [1, 2, 0],\n",
    "             [3, -1, 5]])\n",
    "print(A)\n",
    "print(np.linalg.eigvals(A))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4.20556943  2.31899016 -3.07679678]\n",
      " [ 0.          2.48538322 -0.87692801]\n",
      " [ 0.          0.89947459  1.30904735]]\n"
     ]
    }
   ],
   "source": [
    "niters = 20\n",
    "for k in range(niters):\n",
    "    Q, R = np.linalg.qr(A)\n",
    "    A = R @ Q\n",
    "\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.89721528+0.66545695j 1.89721528-0.66545695j]\n"
     ]
    }
   ],
   "source": [
    "print(np.linalg.eigvals(A[1:, 1:]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Schur form computation\n",
    "\n",
    "- Recall that we are trying to avoid $\\mathcal{O}(n^4)$ complexity. \n",
    "- The idea is to make a matrix have a simpler structure so that each step of QR algorithm becomes cheaper.\n",
    "\n",
    "- In case of a general matrix we can use the **Hessenberg form**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Hessenberg form\n",
    "\n",
    "The matrix $A$ is said to be in the Hessenberg form, if\n",
    "\n",
    "$$a_{ij} = 0, \\quad \\mbox{if } i \\geq j+2.$$\n",
    "\n",
    "$$H = \\begin{bmatrix} * & * & * & * & * \\\\ * & * & * & * & * \\\\ 0 & * & * & * & *\\\\ 0 & 0 & * & * & *\\\\ 0 & 0 & 0 & * & * \\\\ \\end{bmatrix}.$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Reduction any matrix to Hessenberg form\n",
    "\n",
    "- By applying Householder reflections we can reduce any matrix to the Hessenberg form\n",
    "\n",
    "$$U^* A U = H$$\n",
    "\n",
    "- The only difference with Schur decomposition is that we have to map the first column to the vector with two non-zeros, and the first element is not changed.\n",
    "\n",
    "- The computational cost of such reduction is $\\mathcal{O}(n^3)$ operations.\n",
    "\n",
    "- In a Hessenberg form, computation of one iteration of the QR algorithm costs $\\mathcal{O}(n^2)$ operations (e.g. using Givens rotations, how?), and the Hessenberg form is preserved by the QR iteration (check why)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Symmetric (Hermitian) case\n",
    "\n",
    "- In the symmetric case, we have $A = A^*$, then $H = H^*$ and the upper Hessenberg form becomes tridiagonal matrix.\n",
    "\n",
    "- From now on we will talk about the case of symmetric tridiagonal form.\n",
    "\n",
    "- Any symmetric (Hermitian) matrix can be reduced to the tridiagonal form by Householder reflections.\n",
    "\n",
    "- **Key point** is that tridiagonal form is preserved by the QR algorithm, and the cost of one step can be reduced to $\\mathcal{O}(n)$!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR algorithm: iterations\n",
    "\n",
    "- The iterations of the QR algorithm have the following form:\n",
    "\n",
    "$$A_k = Q_k R_k, \\quad A_{k+1} = R_k Q_k.$$\n",
    "\n",
    "- If $A_0 = A$ is <font color='red'> tridiagonal symmetric matrix </font>, this form is preserved by the QR algorithm.\n",
    "\n",
    "Let us see.."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.46901023, -0.82910604,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
       "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ])"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD7CAYAAAChbJLhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACelJREFUeJzt3E1IlP0ax/Hf+BJKJgW5iEyZghYqERW10p3MLioMs5jAGdoLRWH0MjWhbmoj0UJa9UJFKxfSIkGECheBijPULgMX4sJFI2FPzTyL5xzrdDrHmfH+37cz1/cDQvNQM9dNffmP81zeoVwulxMAMyqCHgCAv4geMIboAWOIHjCG6AFjiB4wxvfos9msbty4oe7ubkWjUc3Pz/s9glMnTpxQNBpVNBpVf39/0ONs2MzMjKLRqCRpfn5ePT09Onv2rG7evKlsNhvwdMX79bpSqZTa29vX/t7GxsYCns6tKr9f8PXr1/r27ZueP3+u6elpDQ0N6cGDB36P4cTq6qok6dGjRwFP4o2RkRGNjo6qtrZWkjQ4OKi+vj4dO3ZMN27c0Pj4uDo7OwOesnC/X1c6nVZvb69isVjAk/nD95P+/fv3am9vlyQdPHhQc3Nzfo/gzIcPH/T161fFYjGdP39e09PTQY+0IU1NTRoeHl57nEqldPToUUlSR0eH3r59G9RoG/L7dc3NzWliYkLnzp3T1atXlclkApzOPd+jz2QyqqurW3tcWVmp79+/+z2GEzU1NYrH43r48KFu3bqlS5culfS1RSIRVVX9fDOYy+UUCoUkSVu3btWXL1+CGm1Dfr+uAwcO6PLly3ry5In27Nmj+/fvBzide75HX1dXp5WVlbXH2Wz2P/4CSlk4HNbx48cVCoUUDoe1fft2LS0tBT2WZyoqfv5zWVlZUX19fYDTeKezs1NtbW1rv06n0wFP5Jbv0R86dEiTk5OSpOnpae3fv9/vEZx5+fKlhoaGJEmLi4vKZDJqaGgIeCrvtLS0aGpqSpI0OTmpI0eOBDyRN+LxuGZnZyVJ7969U2tra8ATueX7EdvZ2ak3b97ozJkzyuVyGhgY8HsEZ7q6utTf36+enh6FQiENDAyUzbsYSbpy5YquX7+ue/fuae/evYpEIkGP5IlEIqFkMqnq6mrt3LlTyWQy6JGcCvFTdoAtLOcAxhA9YAzRA8YQPWAM0QPGBBr98+fPg3x5Z7iu0lPO1/Y7oneA6yo95Xxtv+PtPWCML8s5+/bt++MPnlRWVurHjx9FP+/hw4c3MpYzy8vL2rFjR9BjeK5cr0sqj2tbWFhYW5P+f3zZEf3+/bs+f/7s+fOW2w04gI04depUXr+vqOiz2awSiYQ+fvyoLVu26M6dO2pubi7mqQD4rKjv6X+9+83FixfXfrIMwOZXVPTlfPcboNwVFX053/0GKHdFRV/Od78Byl1R0Zfz3W+AclfU8VzOd78Byl1R0VdUVOj27dtezwLAByX9jfi/b8ecD+4KBvyD3XvAGKIHjCF6wBiiB4whesAYogeMIXrAGKIHjCF6wBiiB4zxJfrDhw8rl8vl9eVKKBTK+wsoZ5z0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8Ys+nuhlvIKq6rldl8n5c77KIUcdIDxhA9YAzRA8YQPWAM0QPGED1gDNEDxhA9YAzRA8YQPWDMplvDLUTQK7uFPCcru9gsOOkBY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWNKeiOvEGzvAf/gpAeMIXrAmKLf3p84cULbtm2TJDU2NmpwcNCzoQC4U1T0q6urkqRHjx55OgwA94p6e//hwwd9/fpVsVhM58+f1/T0tNdzAXCkqJO+pqZG8Xhcp0+f1qdPn3ThwgW9evVKVVVm/mcAULKKqjQcDqu5uVmhUEjhcFjbt2/X0tKSdu3a5fV8ADxW1Nv7ly9famhoSJK0uLioTCajhoYGTwcD4EZRJ31XV5f6+/vV09OjUCikgYEB3toDJaKoUrds2aK7d+96PQsAH3A8/0G+a7Au1nULfV5WdlEoNvIAY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWOIHjCGNdwNCPoOu4U+Lyu7kDjpAXOIHjCG6AFjiB4whugBY4geMIboAWOIHjCG6AFjiB4whjVcn7Cyi82Ckx4whugBY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWNYw92EWNmFS5z0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8YwxpuiWNlF4XipAeMIXrAmLyin5mZUTQalSTNz8+rp6dHZ8+e1c2bN5XNZp0OCMBb60Y/MjKia9euaXV1VZI0ODiovr4+PX36VLlcTuPj486HBOCddaNvamrS8PDw2uNUKqWjR49Kkjo6OvT27Vt30wHw3LrRRyIRVVX9/JA/l8utfVq7detWffnyxd10ADxX8Ad5FRU//8jKyorq6+s9HQiAWwVH39LSoqmpKUnS5OSkjhw54vlQANwpOPorV65oeHhY3d3d+uuvvxSJRFzMBcCRvDbyGhsb9eLFC0lSOBzW48ePnQ4FwB3WcA1hZRcSG3mAOUQPGEP0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxjDGi7+iJXd8sVJDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8YQ/SAMazhYsNKaWWXdV1OesAcogeMIXrAGKIHjCF6wBiiB4whesAYogeMIXrAGKIHjGENF74KemWXO+xy0gPmED1gDNEDxhA9YAzRA8YQPWAM0QPGED1gDNEDxrCRh02L7T03OOkBY4geMCav6GdmZhSNRiVJqVRK7e3tikajikajGhsbczogAG+t+z39yMiIRkdHVVtbK0lKp9Pq7e1VLBZzPhwA76170jc1NWl4eHjt8dzcnCYmJnTu3DldvXpVmUzG6YAAvLVu9JFIRFVVP98QHDhwQJcvX9aTJ0+0Z88e3b9/3+mAALxV8Ad5nZ2damtrW/t1Op32fCgA7hQcfTwe1+zsrCTp3bt3am1t9XwoAO4UvJyTSCSUTCZVXV2tnTt3KplMupgLgCN5Rd/Y2KgXL15IklpbW/Xs2TOnQwFwhzVclIV812BdrOsW+rxBr+yykQcYQ/SAMUQPGEP0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDGs4cKUoO+wW+jzuljZ5aQHjCF6wBiiB4whesAYogeMIXrAGKIHjCF6wBiiB4whesAYogf+h1wul/eXK6FQKO+vfBE9YAzRA8YQPWAM0QPGED1gDNEDxhA9YAzRA8YQPWAM0QPGcDdcwAOb4S67+eKkB4whesAYogeMIXrAGKIHjCF6wBiiB4whesAYogeMIXrAmFDO5a08/+XYsWPavXv3f/335eVl7dixw/XL+47rKj3lcG0LCwuamppa/zfmAnTy5MkgX94Zrqv0lPO1/Y6394AxRA8YU5lIJBJBDtDW1hbkyzvDdZWecr62X/nyQR6AzYO394AxRA8YQ/SAMUQPGEP0gDF/AzEElrIFHqU3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD7CAYAAAChbJLhAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAACelJREFUeJzt3E1IlP0ax/Hf+BJKJgW5iEyZghYqERW10p3MLioMs5jAGdoLRWH0MjWhbmoj0UJa9UJFKxfSIkGECheBijPULgMX4sJFI2FPzTyL5xzrdDrHmfH+37cz1/cDQvNQM9dNffmP81zeoVwulxMAMyqCHgCAv4geMIboAWOIHjCG6AFjiB4wxvfos9msbty4oe7ubkWjUc3Pz/s9glMnTpxQNBpVNBpVf39/0ONs2MzMjKLRqCRpfn5ePT09Onv2rG7evKlsNhvwdMX79bpSqZTa29vX/t7GxsYCns6tKr9f8PXr1/r27ZueP3+u6elpDQ0N6cGDB36P4cTq6qok6dGjRwFP4o2RkRGNjo6qtrZWkjQ4OKi+vj4dO3ZMN27c0Pj4uDo7OwOesnC/X1c6nVZvb69isVjAk/nD95P+/fv3am9vlyQdPHhQc3Nzfo/gzIcPH/T161fFYjGdP39e09PTQY+0IU1NTRoeHl57nEqldPToUUlSR0eH3r59G9RoG/L7dc3NzWliYkLnzp3T1atXlclkApzOPd+jz2QyqqurW3tcWVmp79+/+z2GEzU1NYrH43r48KFu3bqlS5culfS1RSIRVVX9fDOYy+UUCoUkSVu3btWXL1+CGm1Dfr+uAwcO6PLly3ry5In27Nmj+/fvBzide75HX1dXp5WVlbXH2Wz2P/4CSlk4HNbx48cVCoUUDoe1fft2LS0tBT2WZyoqfv5zWVlZUX19fYDTeKezs1NtbW1rv06n0wFP5Jbv0R86dEiTk5OSpOnpae3fv9/vEZx5+fKlhoaGJEmLi4vKZDJqaGgIeCrvtLS0aGpqSpI0OTmpI0eOBDyRN+LxuGZnZyVJ7969U2tra8ATueX7EdvZ2ak3b97ozJkzyuVyGhgY8HsEZ7q6utTf36+enh6FQiENDAyUzbsYSbpy5YquX7+ue/fuae/evYpEIkGP5IlEIqFkMqnq6mrt3LlTyWQy6JGcCvFTdoAtLOcAxhA9YAzRA8YQPWAM0QPGBBr98+fPg3x5Z7iu0lPO1/Y7oneA6yo95Xxtv+PtPWCML8s5+/bt++MPnlRWVurHjx9FP+/hw4c3MpYzy8vL2rFjR9BjeK5cr0sqj2tbWFhYW5P+f3zZEf3+/bs+f/7s+fOW2w04gI04depUXr+vqOiz2awSiYQ+fvyoLVu26M6dO2pubi7mqQD4rKjv6X+9+83FixfXfrIMwOZXVPTlfPcboNwVFX053/0GKHdFRV/Od78Byl1R0Zfz3W+AclfU8VzOd78Byl1R0VdUVOj27dtezwLAByX9jfi/b8ecD+4KBvyD3XvAGKIHjCF6wBiiB4whesAYogeMIXrAGKIHjCF6wBiiB4zxJfrDhw8rl8vl9eVKKBTK+wsoZ5z0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8Ys+nuhlvIKq6rldl8n5c77KIUcdIDxhA9YAzRA8YQPWAM0QPGED1gDNEDxhA9YAzRA8YQPWDMplvDLUTQK7uFPCcru9gsOOkBY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWNKeiOvEGzvAf/gpAeMIXrAmKLf3p84cULbtm2TJDU2NmpwcNCzoQC4U1T0q6urkqRHjx55OgwA94p6e//hwwd9/fpVsVhM58+f1/T0tNdzAXCkqJO+pqZG8Xhcp0+f1qdPn3ThwgW9evVKVVVm/mcAULKKqjQcDqu5uVmhUEjhcFjbt2/X0tKSdu3a5fV8ADxW1Nv7ly9famhoSJK0uLioTCajhoYGTwcD4EZRJ31XV5f6+/vV09OjUCikgYEB3toDJaKoUrds2aK7d+96PQsAH3A8/0G+a7Au1nULfV5WdlEoNvIAY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWOIHjCGNdwNCPoOu4U+Lyu7kDjpAXOIHjCG6AFjiB4whugBY4geMIboAWOIHjCG6AFjiB4whjVcn7Cyi82Ckx4whugBY4geMIboAWOIHjCG6AFjiB4whugBY4geMIboAWNYw92EWNmFS5z0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8YwxpuiWNlF4XipAeMIXrAmLyin5mZUTQalSTNz8+rp6dHZ8+e1c2bN5XNZp0OCMBb60Y/MjKia9euaXV1VZI0ODiovr4+PX36VLlcTuPj486HBOCddaNvamrS8PDw2uNUKqWjR49Kkjo6OvT27Vt30wHw3LrRRyIRVVX9/JA/l8utfVq7detWffnyxd10ADxX8Ad5FRU//8jKyorq6+s9HQiAWwVH39LSoqmpKUnS5OSkjhw54vlQANwpOPorV65oeHhY3d3d+uuvvxSJRFzMBcCRvDbyGhsb9eLFC0lSOBzW48ePnQ4FwB3WcA1hZRcSG3mAOUQPGEP0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxjDGi7+iJXd8sVJDxhD9IAxRA8YQ/SAMUQPGEP0gDFEDxhD9IAxRA8YQ/SAMazhYsNKaWWXdV1OesAcogeMIXrAGKIHjCF6wBiiB4whesAYogeMIXrAGKIHjGENF74KemWXO+xy0gPmED1gDNEDxhA9YAzRA8YQPWAM0QPGED1gDNEDxrCRh02L7T03OOkBY4geMCav6GdmZhSNRiVJqVRK7e3tikajikajGhsbczogAG+t+z39yMiIRkdHVVtbK0lKp9Pq7e1VLBZzPhwA76170jc1NWl4eHjt8dzcnCYmJnTu3DldvXpVmUzG6YAAvLVu9JFIRFVVP98QHDhwQJcvX9aTJ0+0Z88e3b9/3+mAALxV8Ad5nZ2damtrW/t1Op32fCgA7hQcfTwe1+zsrCTp3bt3am1t9XwoAO4UvJyTSCSUTCZVXV2tnTt3KplMupgLgCN5Rd/Y2KgXL15IklpbW/Xs2TOnQwFwhzVclIV812BdrOsW+rxBr+yykQcYQ/SAMUQPGEP0gDFEDxhD9IAxRA8YQ/SAMUQPGEP0gDGs4cKUoO+wW+jzuljZ5aQHjCF6wBiiB4whesAYogeMIXrAGKIHjCF6wBiiB4whesAYogf+h1wul/eXK6FQKO+vfBE9YAzRA8YQPWAM0QPGED1gDNEDxhA9YAzRA8YQPWAM0QPGcDdcwAOb4S67+eKkB4whesAYogeMIXrAGKIHjCF6wBiiB4whesAYogeMIXrAmFDO5a08/+XYsWPavXv3f/335eVl7dixw/XL+47rKj3lcG0LCwuamppa/zfmAnTy5MkgX94Zrqv0lPO1/Y6394AxRA8YU5lIJBJBDtDW1hbkyzvDdZWecr62X/nyQR6AzYO394AxRA8YQ/SAMUQPGEP0gDF/AzEElrIFHqU3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "sns.set_style('white')\n",
    "#Generate a random tridiagonal matrix\n",
    "\n",
    "n = 20\n",
    "d = np.random.randn(n)\n",
    "sub_diag = np.random.randn(n-1)\n",
    "\n",
    "mat = np.diag(d) + np.diag(sub_diag, -1) + np.diag(sub_diag, 1)\n",
    "mat1 = np.abs(mat)\n",
    "mat1 = mat1/np.max(mat1.flatten())\n",
    "plt.spy(mat)\n",
    "q, r = np.linalg.qr(mat)\n",
    "plt.figure()\n",
    "b = r.dot(q)\n",
    "b[abs(b) <= 1e-12] = 0\n",
    "plt.spy(b)\n",
    "#plt.figure()\n",
    "#plt.imshow(np.abs(r.dot(q)))\n",
    "b[0, :]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tridiagonal form\n",
    "\n",
    "- In the tridiagonal form, you do not have to compute the $Q$ matrix: you only have to compute the **triadiagonal part** that appears after the iterations \n",
    "\n",
    "$$A_k = Q_k R_k, \\quad A_{k+1}  = R_k Q_k,$$\n",
    "\n",
    "in the case when $A_k = A^*_k$ and is also tridiagonal.\n",
    "\n",
    "- Such matrix is defined by $\\mathcal{O}(n)$ parameters; computation of the QR is more complicated, but it is possible to compute $A_{k+1}$ directly without computing $Q_k$.\n",
    "\n",
    "- This is called **implicit QR-step**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Theorem on implicit QR iteration\n",
    "\n",
    "All the implicit QR algorithms are based on the following **theorem**: \n",
    "\n",
    "Let \n",
    "\n",
    "$$Q^* A Q = H$$ \n",
    "\n",
    "be an irreducible upper Hessenberg matrix. Then, the first column of the matrix $Q$ defines all of its other columns. \n",
    "It can be found from the equation\n",
    "\n",
    "$$A Q = Q H. $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convergence of the QR-algorithm\n",
    "\n",
    "- The convergence of the QR-algorithm is a very **delicate issue** (see E. E. Tyrtyshnikov, \"Brief introduction to numerical analysis\" for details). \n",
    "\n",
    "\n",
    "**Summary.** If we have a decomposition of the form\n",
    "\n",
    "$$A = X \\Lambda X^{-1}, \\quad A = \\begin{bmatrix}A_{11} & A_{12} \\\\ A_{21} & A_{22}\\end{bmatrix}$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "\\Lambda = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ \n",
    "0 & \\Lambda_2 \\end{bmatrix}, \\quad \\lambda(\\Lambda_1)=\\{\\lambda_1,\\dots,\\lambda_m\\}, \\ \\lambda(\\Lambda_2)=\\{\\lambda_{m+1},\\dots,\\lambda_r\\},\n",
    "$$\n",
    "\n",
    "and there is a **gap** between the eigenvalues of $\\Lambda_1$ and $\\Lambda_2$ ($|\\lambda_1|\\geq \\dots \\geq |\\lambda_m| > |\\lambda_{m+1}| \\geq\\dots \\geq |\\lambda_r| >0$), then the $A^{(k)}_{21}$ block of $A_k$\n",
    "in the QR-iteration goes to zero with \n",
    "\n",
    "$$\\Vert A^{(k)}_{21} \\Vert \\leq  C q^k, \\quad q = \\left| \\frac{\\lambda_{m+1}}{\\lambda_{m}}  \\right |,$$\n",
    "\n",
    "where $m$ is the size of $\\Lambda_1$.\n",
    "\n",
    "So we need to increase the gap!  It can be done by the **QR algorithm with shifts**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## QR-algorithm with shifts\n",
    "\n",
    "$$A_{k} - s_k I = Q_k R_k, \\quad A_{k+1} = R_k Q_k + s_k I$$\n",
    "\n",
    "\n",
    "The convergence rate for a shifted version is then\n",
    "\n",
    "$$\\left| \\frac{\\lambda_{m+1} - s_k}{\\lambda_{m} - s_k}  \\right |,$$\n",
    "\n",
    "where $\\lambda_m$ is the $m$-th largest eigenvalue of the matrix in modulus. \n",
    "- If the shift is close to the eigenvalue, then the convergence speed is better.\n",
    "- There are different strategies to choose shifts. \n",
    "- Introducing shifts is a general strategy to improve convergence of iterative methods of finding eigenvalues. \n",
    "- In next slides we will illustrate how to choose shift on a simpler algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Shifts and power method\n",
    "\n",
    "- Remember the power method for the computation of the eigenvalues.\n",
    "\n",
    "$$x_{k+1} := A x_k, \\quad x_{k+1} := \\frac{x_{k+1}}{\\Vert x_{k+1} \\Vert}.$$\n",
    "\n",
    "- It converges to the eigenvector corresponding to the largest eigenvalue in modulus. \n",
    "\n",
    "- The convergence can be very slow.\n",
    "\n",
    "- Let us try to use shifting strategy. If we shift the matrix as \n",
    "\n",
    "$$  A := A - \\lambda_k I$$\n",
    "\n",
    "and the corresponding eigenvalue becomes small (but we need large). \n",
    "- That is not what we wanted!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Inverse iteration and Rayleigh quotient iteration\n",
    "\n",
    "- To make a small eigenvalue large, we need to **invert the matrix**, and that gives us **inverse iteration**\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda I)^{-1} x_k,$$\n",
    "\n",
    "where $\\lambda$ is the shift which is approximation to the eigenvalue that we want. \n",
    "\n",
    "- As it was for the power method, the convergence is linear.\n",
    "\n",
    "- To accelerate convergence one can use the **Rayleigh quotient iteration** (inverse iteration with adaptive shifts) which is given by the selection of the **adaptive shift**:\n",
    "\n",
    "$$x_{k+1} = (A - \\lambda_k I)^{-1} x_k,$$\n",
    "\n",
    "$$\\lambda_k = \\frac{(Ax_k, x_k)}{(x_k, x_k)}$$\n",
    "\n",
    "- In the symmetric case $A = A^*$ the convergence is **locally cubic** and **locally quadratic** otherwise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular values and eigenvalues (1)\n",
    "\n",
    "- Now let us talk about singular values and eigenvalues. \n",
    "\n",
    "- SVD \n",
    "\n",
    "$$A = U \\Sigma V^*$$\n",
    "\n",
    "exists for any matrix.\n",
    "\n",
    "- It can be also viewed as a reduction of a given matrix to the diagonal form by means of two-sided unitary transformations:\n",
    "\n",
    "$$\\Sigma = U^* A V.$$\n",
    "\n",
    "- By two-sided Householder transformation we can reduce any matrix to the **bidiagonal form** $B$ (how?)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Singular values and eigenvalues (2)\n",
    "\n",
    "- **Implicit QR algorithm** (with shifts) gives the way of computing the eigenvalues (and Schur form).\n",
    "\n",
    "- But we cannot apply QR algorithm directly to the bidiagonal matrix, as it is not diagonalizable in general case.\n",
    "\n",
    "- However, the problem of the computation of the SVD can be reduced to the **symmetric eigenvalue problem** in two ways:\n",
    "\n",
    "\n",
    "1. Work with the tridiagonal matrix \n",
    "\n",
    "$$T = B^* B$$ \n",
    "\n",
    "2. Work with the extended matrix \n",
    "\n",
    "$$T = \\begin{bmatrix} 0 & B \\\\ B^* & 0 \\end{bmatrix}$$\n",
    "\n",
    "- The case 1 is OK if you **do not form T directly**!\n",
    "\n",
    "- Thus, the problem of computing singular values can be reduced to the problem of the computation of the eigenvalues of symmetric tridiagonal matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Algorithms for the SEV  (symmetric eigenvalue problem)\n",
    "\n",
    "Done:\n",
    "- QR algorithm: the \"gold standard\" of the eigenvalue computations\n",
    "- RQI-iteration: Rayleigh quotient iteration is implicitly performed at each step of the QR algorithm\n",
    "\n",
    "Next slides:\n",
    "- Divide-and-conquer\n",
    "- Bisection method\n",
    "- Jacobi method"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Divide-and-conquer\n",
    "\n",
    "- Suppose we have a tridiagonal matrix, and we split it into two blocks:\n",
    "\n",
    "$$T = \\begin{bmatrix} T'_1 & B \\\\ B^{\\top} & T'_2 \\end{bmatrix}$$\n",
    "\n",
    "- We can write the matrix $T$ as\n",
    "\n",
    "$$T = \\begin{bmatrix} T_1 & 0 \\\\ 0 & T_2 \\end{bmatrix} + b_m v v^*$$\n",
    "\n",
    "where $vv^*$ is rank 1 matrix, $v = (0,\\dots,0,1,1,0,\\dots,0)^T$.\n",
    "\n",
    "- Suppose we have decomposed $T_1$ and $T_2$ already:\n",
    "\n",
    "$$T_1 = Q_1 \\Lambda_1 Q^*_1, \\quad T_2 = Q_2 \\Lambda_2 Q^*_2$$\n",
    "\n",
    "- Then (check),\n",
    "\n",
    "$$\\begin{bmatrix} Q^*_1 & 0 \\\\ 0 & Q^*_2 \\end{bmatrix} T\\begin{bmatrix} Q_1 & 0 \\\\ 0 & Q_2 \\end{bmatrix} = D + \\rho u u^{*}, \\quad D = \\begin{bmatrix} \\Lambda_1 & 0 \\\\ 0 & \\Lambda_2\\end{bmatrix}$$\n",
    "\n",
    "- I.e. we have reduced the problem to the problem of the computation of the eigenvalues of <font color='red'> diagonal plus low-rank matrix </font>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Diagonal-plus-low-rank matrix\n",
    "\n",
    "It is tricky to compute the eigenvalues of the matrix\n",
    "\n",
    "$$D + \\rho u u^* $$\n",
    "\n",
    "The characteristic polynomial has the form\n",
    "\n",
    "$$\\det(D + \\rho uu^* - \\lambda I) = \\det(D - \\lambda I)\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 0.$$\n",
    "\n",
    "Then (prove!!) \n",
    "\n",
    "$$\\det(I + \\rho (D - \\lambda I)^{-1} uu^*) = 1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "Hint: find $\\det(I + w u^*)$ using the fact that $\\text{det}(C) = \\prod_{i=1}^n\\lambda_i(C)$ and $\\text{trace}(C) = \\sum_{i=1}^n \\lambda_i$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Characteristic equation\n",
    "\n",
    "$$1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{d_i - \\lambda} = 0$$\n",
    "\n",
    "How to find the roots?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtgAAAJMCAYAAADAJNFQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xd4VFXiPvA3jSSQhARS6J1cCC10VBQERQRZpImAX6W5ICjs2oW1oCiKKO4C6hpQQEARBTS0qFSB0AyhJlcgIaSQnpCemWTu74+Q/GSlk3PPzL3v53nmSZ95uUMm75w59xwnTdNARERERETVw1l2ACIiIiIiI2HBJiIiIiKqRizYRERERETViAWbiIiIiKgasWATEREREVUjFmwiIiIiomrkqseNKIriDGAygPEA2gFwAxALIAzA56qqcq1AIiIiIjIEJ9HrYCuK4gFgA4CBAGyoKNZeAJpc/pa1AMawZBMRERGREegxReQDVJTrRACdVVVtp6pqUwBDABQCGA1gnA45iIiIiIiEE1qwFUVpAWA6gDIAD6uqerzya6qqbgLw0eUPJ4rMQURERESkF9FzsMcAcAGwXFXVU1f5+lcASgFcEJyDiIiIiEgXogt2/8tvf7zaF1VVPQ/gPcEZiIiIiIh0I7pgt7/8NkZRlNoAJgC4DxUnOZ4G8IWqqqcFZyAiIiIi0o2wVUQurx5SfPnDBwEsB9Dwf76tDMA0VVXDhIQgIiIiItKZyJMcvf/0/jeoKNsPA/AE0BjAQlSMoH+uKEo/gTmIiIiIiHQjcoqIx5/erwWgh6qq8Zc/TgLwvKIogahYou89AL2ud2Xdu3fXAgIChAS9Hk3T4OTkpPvtkv4c+b5OSUlBfn4+FEWp+lxiYiJKS0vRqlUricns163c31lZWUhNTUXbtm3h7FwxLpGZmYm0tDSEhIQ47P8be5GTk4OUlBQEBwfDzc0N6enpyMjIQLt27arl+h35d7s65Ofn48KFC2jRogU8PT2Rl5eHxMREtGzZEh4eHje+Agdi9vvaTGTd1+fOnctUVfXGhVTTNCGX4OBg/+DgYO3y5b/X+J6QP31P4PWub9iwYZoMp0+flnK7pD9Hvq+ffPJJrVmzZld8bsSIEVq7du0kJbJ/t3J/v//++xoArbCwsOpz8+bN0wBoxcXFIuKZyueff64B0JKTkzVN07Q33nhDq/jzVD0c+Xe7OmzcuFEDoEVFRWmapmk//PCDBkA7duyY5GTVz+z3tZnIuq+Dg4OPaDfRg0VOEckDUDnB+/g1vucPANbL7zcTmIXI0CwWC2rUqCE7hmFZrRUPU25ubpKTGBOPr1iVx9fVVfS6BkRUSVjBVlXVAiD+Bt+m4f+XcOv1vpGIro0FWywWFLFYsMUqKysDwONLpCfRW6Ufuvy22zW+3hRADQA2AOcFZyEyLBZssSwWC9zc3Di3UxAWbLH4BJFIf6IL9trLb0cqivK/S/QBFduoA8BuVVVzBGchMiwWbLGsVivLn0As2GJxBJtIf6IL9k8AIlGxscwmRVFaVn5BUZTRAJ69/OG7gnMQGRoLtlgs2GKxYItVWbA5gk2kH6G/baqq2hRFGQVgO4BQAKqiKKdRUbibX/6211VV3S4yB5HRWSwW1KxZU3YMw7JarXwCI1BZWRlcXFw4BUcQPoEh0p/oEWyoqpoMoAuANwDEAmiFik1oIgAMVFV1rugMREbHEWyxKudgkxhWq5WjqwJxBJtIf7r8tqmqWgTgncsXIqpmLNhicYqIWDy+YvEkRyL9CR/BJiLxWLDFYgEUi8dXLJ7kSKQ/FmwiA2DBFovHVywWbLE4RYRIfyzYRAbAOcJisQCKxeMrFqeIEOmPBZvIADjCKhYLoFg8vmJxlRYi/bFgExkAC7ZYXKZPLBZssXh8ifTHgk1kAKWlpfDw8JAdw7A4BUessrIyTl8QiMeXSH8s2EQOTtM0lJSUsGALxBFAsXh8xeLxJdIfCzaRg7NardA0jQVbIBYUsXh8xeIINpH+WLCJHFxJSQkAsGALxDnuYrFgi8WCTaQ/FmwiB1dZsN3d3SUnMS4WQLF4fMXi8SXSHws2kYPjCLZ4LChi8fiKxRFsIv2xYBM5uNLSUgAs2CKxAIrF4ysWjy+R/liwiRwcR7DF4xxssVgAxeIINpH+WLCJHBwLtngsgGKxAIpVVlbG/79EOmPBJnJwLNjicZ1xsfgERiyr1conMEQ6Y8EmcnAs2OKVlpZylRaBWLDF4isERPpjwSZycJUnObIAimGz2WC1Wnl8BWLBFovHl0h/LNhEDo4j2GJZLBYAfAIjEgugWBzBJtIfCzaRg2PBFovHVzxOwRGLJzkS6Y8Fm8jBsQCKxSk44rFgi8WTHIn0x4JN5OBYsMViwRaPBVssjmAT6Y8Fm8jBcSdHsfgERiybzYaysjJu5CMQR7CJ9MeCTeTgWADF4gi2WDyJVDyOYBPpjwWbyMGVlJTA2dmZI1SCsGCLxeMrHkewifTHgk3k4LjLoFgsgGLx+IrHZRCJ9MeCTeTgiouLUbNmTdkxDItTcMRiwRbPYrHw+BLpjAWbyMEVFRXB09NTdgzDYgEUi8dXvNLSUp5ESqQzFmwiB1dUVMQRbIFYAMWqPMmRBVAcjmAT6Y8Fm8jBsWCLxSkiYvEJjHgcwSbSHws2kYMrLi7mFBGBWADF4vEVq7y8HOXl5SzYRDpjwSZycBzBFosFUCweX7G4zjiRHCzYRA7ueic5apqmcxrjqZwiwoIiBgu2WNeb487HByJxWLCJHNy1lulzcnKSkMZ4uBW9WJXHl1MYxLhaweZjA5F4LNhEDo5TRMTiCKtYnMIgFv//EsnBgk3k4LgOtlilpaXcil4gFkCxuAwikRws2EQOjiPYYnGVFrFYsMXi8SWSgwWbyIFpmsat0gVjwRaLBVAsjmATycGCTeTAKle4YAEUhwVbLJ7kKBYLNpEcLNhEDqyoqAgAOIItEAu2WDzJUSy+QkAkBws2kQMrLi4GwIItEgu2WCyAYnEEm0gOFmwiB1Y5gs0CKA5XaRGrtLQUTk5OXKVFEE7BIZKDBZvIgRUWFgIAatWqJTmJcXEEW6zK48vNT8TgFBwiOViwiRxYQUEBAMDLy0tyEuNiwRarpKSEu2QKxBFsIjlYsIkcGAu2eFwGUSw+gRGLI9hEcrBgEzkwThERjwVQLB5fsXiSI5EcLNhEDowj2OKxAIrFKSJicYoIkRws2EQOjAVbPBZssXh8xeIUESI5WLCJHFjlFBEWbHFYAMXiCLZYHMEmkoMFm8iBFRQUwMnJiQVFkLKyMlitVhZsgfgERizOwSaSgwWbyIEVFBTAy8uLawgLUrlTJgugOCzYYpWUlMDZ2Zkb+RDpjAWbyIEVFhZyeohALNjicYqIWJXHl0/CifTFgk3kwCpHsEmMyoLNdbDF4Qi2WCUlJTy+RBKwYBM5sIKCAq6BLVBRUREAjmCLxBFssXh8ieRgwSZyYJwiIhY38hGPI9hisWATycGCTeTA8vPzWbAFYsEWr7i4mAVQIBZsIjlYsIkcWF5eHnx8fGTHMCwWbLHKyspQXl7OEWyBWLCJ5GDBJnJgLNhisWCLxVVaxOMrBERysGATOTAWbLG4U6ZYJSUlAMACKBBHsInkYMEmclDl5eUoLCxkwRaII9hicQRbPBZsIjlYsIkcVH5+PgCwYAvEgi1W5TKIXGdcHBZsIjlYsIkcVF5eHgAWbJEqCzZHWMVgwRaPBZtIDhZsIgfFgi1eYWEhatasCWdnPlSKwFcIxONOjkRy8K8GkYNiwRavsLCQ5U8gjmCLxxFsIjlYsIkcFAu2eCzYYnEEWzwWbCI5WLCJHBQLtngs2GJxBFssTdNYsIkkYcEmclAs2OKxYItVWbB5jMUoKyuDzWZjwSaSgAWbyEHl5OQAAPz8/CQnMa6bKdiapumUxngqp4hwBFuMynXGWbCJ9MeCTeSgsrOz4ebmxtE/gfLz86/5CoGTk5POaYznWiPYPLbVgztlEsnDgk3koHJycuDn58cyIlBeXh68vb1lxzCswsJCuLq6ws3NTXYUQ+IINpE8LNhEDio7Oxt16tSRHcPQ8vPzWbAFKioq4iswAlUWbB5jIv25yrhRRVHuBvAbgERVVZvJyEDk6CpHsEkcFmyxKjfyITG4SguRPLqPYCuK4g5gmYzbJjISjmCLZbFYUFpayoItEEewxWLBJpJHRsmdA6CNhNslMpTs7OwbjmBzhYvbl5+fDwAs2AJxBFusGxVsPj4QiaNrwVYUpQuAFwAU63m7REaUk5Nz3RFsnvx4Z1iwxSsqKmLBFqiyYHt6el7xeT42EImnW8FWFMUNwFcANABz9bpdIiMqLy/HpUuXOEVEoMqCzY18xCkoKOATGIE4RYRIHj1HsGcB6AjgAwDHdbxdIsPJzc0FwE1mROIItnj5+fnw8vKSHcOwWLCJ5NGlYCuK0h4VBTsWHL0mumPZ2dkAwBFsgSq3omfBFocj2GKxYBPJI7xgK4riAuBLAG4AJquqWir6NomMrrJgcwRbHI5gi8cRbLFYsInk0WME+wUA3QF8qqrqPh1uj8jwcnJyAHAEWySOYIvHdcbFqizY3MmRSH9CN5pRFKU1KpblSwTw2p1cl8ViQUxMTLXkuhUlJSVSbpf050j39YkTJwBUjGRfK3NeXh5KS0sd5t+ktxvd36qqAgDS09OrdsT7s/T0dABAbGwsC8xtsFgssFgsKC4u/sv9kJGRAQDV9n/XkX63q1NiYiI8PT0RGxv7l88DQFxcHNzd3WVEE8as97UZ2ft9LaxgK4rihIoNZTwAPKOqav6dXF+NGjXQtm3basl2K2JiYqTcLunPke7rHTt2AAC6d++OwMDAq36Pj48P3N3dHebfpLcb3d81atSAs7MzunXrBmfnv77YV3nc27Rp85dl0OjGKqc5tWjR4i/3Q0BAAABU2/9dR/rdrk41a9ZErVq1/vJvr3zyeLVj7+jMel+bkb3f1yKniEwHcC+Ab1RV3SzwdohMh3OwxcvJyYGvr+9VyzXduYKCAgDgHGyBuM44kTwip4iMvPx2jKIoY67xPU0VRancSqq5qqrnBeYhMozs7Gx4eXnBzc1NdhTDqizYJAZPIhWPBZtIHpEF+8R1rt8PQAiAUgBHLn+uRGAWIkPJzs7mCY6C5ebm8hUCgSoLNkewxSkqKuL0JSJJhBVsVVWfu9bXFEV5BEA4gFRVVXuLykBkVGlpadece03VIycnhwVboMopIhzBFqewsJBPYIgk4eRCIgeUlpaGoKAg2TEMjQVbLI5gi1dQUMDjSyQJCzaRA2LBFi83N5dzsAXiHGzxWLCJ5GHBJnIwNpsN6enpLNgCaZrGEWzBcnNzAXAlHJFYsInkEbrRzLWoqroJgJOM2yZydNnZ2SgvL0e9evVkRzGs4uJiWCwWlj+BKgt27dq1JScxLm5FTyQPR7CJHExaWhoAcARbII6uipebmwsvLy+4ukoZ5zE8TdM4gk0kEQs2kYNhwRYvJycHADgHWyDOcRfLYrGgrKyMBZtIEhZsIgeTmpoKgAVbpMqCzRFscViwxeJOmURysWATORiOYIvHgi0eC7ZYLNhEcrFgEzmYtLQ0uLm5sfwJVDkHmwVQHBZssViwieRiwSZyMJW7ODo5cSEeUTiCLR4Ltlgs2ERysWATOZi0tDQu0ScYT3IUjwVbLBZsIrlYsIkcDHdxFC83Nxfe3t5cQk4Qm82GS5cusWALxIJNJBcLNpGDSUlJQf369WXHMLTs7GxODxGooKAANpuNBVugvLw8ANyKnkgWFmwiB1JaWorU1FQ0btxYdhRD41b0YvEkUvEuXboEgDtlEsnCgk3kQJKTkwEATZo0kZzE2NLT0xEYGCg7hmGxYIvHgk0kFws2kQO5cOECABZs0ViwxWLBFi83NxceHh5wd3eXHYXIlFiwiRxIYmIiAHCKiECaprFgC8aCLd6lS5c4ek0kEQs2kQOpHMFmwRYnLy8PFouFBVsgFmzxuEoLkVws2EQOJDExEf7+/vD09JQdxbDS09MBgAVbIBZs8XJzczmCTSQRCzaRA7lw4QLnXwvGgi1eZcFmARSHU0SI5GLBJnIgiYmJnB4iGAu2eLm5ufDy8uJGPgJxigiRXCzYRA6EI9jisWCLl5WVxY18BOMUESK5WLCJHMSlS5eQl5fHEWzBKgu2v7+/5CTGlZGRwScwgnGKCJFcLNhEDiIuLg4A0KJFC8lJjC09PR1+fn6oUaOG7CiGlZGRgYCAANkxDMtisaC4uJhTRIgkYsEmchBnzpwBALRq1UpyEmPjGtjisWCLxV0cieRjwSZyEGfPngXAgi0aC7Z4LNhisWATyceCTeQgzpw5gwYNGqBWrVqyoxhaWloaC7ZAhYWFKCoq4jEWqLJgc4oIkTws2EQO4syZMxy91kFSUhIaNmwoO4ZhZWRkAABHsAXiOuNE8rFgEzmIs2fPonXr1rf0M5qmCUpjTHl5ecjPz+dKLQKxYIt3syPYfHwgEocFm8gB5OXlIS0t7ZZGsJ2cnAQmMqakpCQAQKNGjSQnMS4WbPFuNAebjw1E4rFgEzmAc+fOAcAtj2DTrWHBFo8FWzxOESGSjwWbyAH88ccfAFiwRWPBFq+yYPMkR3EqR7B9fHwkJyEyLxZsIgdw8uRJuLi4QFEU2VEMrbJgN2jQQHIS40pPT4e7uzu8vLxkRzGs3NxceHt7w8XFRXYUItNiwSZyACdOnEBwcDDc3d1lRzG0pKQkBAUFcRdHgSrXwOY8YHEyMzPh7+8vOwaRqbFgEzmAkydPon379rJjGF5SUhKnhwjGTWbEY8Emko8Fm8jOFRYWIi4ujgVbB4mJiSzYgrFgi8djTCQfCzaRnTt9+jQ0TUOHDh1kRzG8pKQkroEt2MWLF1G/fn3ZMQwtMzOTBZtIMhZsIjt34sQJAOAItmAFBQXIzc3lLo4ClZeXIyUlhcdYsIyMDE4RIZKMBZvIzh0/fhyenp5o0aKF7CiGFh8fDwBo3ry55CTGlZ6ejvLychZsgQoLC1FcXMwRbCLJWLCJ7NyRI0cQGhrKJbcEq9zMp2XLlpKTGFdycjIAsGALlJmZCYAb+RDJxoJNZMfKyspw9OhRdO/eXXYUw6ss2HylQJzKgs0TScWp3MiHU0SI5GLBJrJjMTExKCoqYsHWQVxcHHx9fVGnTh3ZUQyLI9jicQSbyD6wYBPZscOHDwMAC7YOzp07x+khgiUlJcHV1ZXbpAvEEWwi+8CCTWTHjhw5Ah8fH7Ru3Vp2FMNjwRYvOTkZ9evXh7Mz//SIUlmwOYJNJBcf5Yjs2MGDB9GtWzcWEsHKyspw/vx5zr8WLDk5mdNDBMvMzISrqytq164tOwqRqfGvNpGdysvLQ3R0NHr37i07iuElJiairKyMI9iCJScn8wRHwSrXwHZycpIdhcjUWLCJ7NT+/fths9lw3333yY5ieFyiTx8cwRaPuzgS2QcWbCI7tWfPHri6uqJXr16yoxieqqoAAEVRJCcxrkuXLiE/P58FWzDu4khkH1iwiezUnj170K1bN9SqVUt2FMM7ffo0ateujfr168uOYljcKVMfFy9eRL169WTHIDI9FmwiO1RYWIhDhw7h3nvvlR3FFE6fPo2QkBDOWxWosmDzRFJxNE1DSkoKXyUgsgMs2ER2aNeuXbBarXjwwQdlRzGFyoJN4sTFxQFgwRYpJycHJSUlLNhEdoAFm8gObdu2DZ6enhzB1kFmZibS09NZsAWLi4uDn58ffH19ZUcxrJSUFABAgwYNJCchIhZsIjsUERGB+++/Hx4eHrKjGF5MTAwAsGALFhcXx/nXgnEreiL7wYJNZGfi4uJw5swZDBw4UHYUUzh9+jQAFmzR4uLiOD1EMI5gE9kPFmwiO7Nx40YAwKBBgyQnMYcTJ07A29sbjRs3lh3FsGw2G3fK1EHlCDYLNpF8LNhEdmb9+vXo1KkTNz3RSVRUFEJDQ7mCiEApKSmwWCws2IKlpKSgbt26cHd3lx2FyPRYsInsyMWLF7F//34MHz5cdhRTKC8vx7Fjx9ClSxfZUQytcgURzsEWiztlEtkPFmwiO7JhwwZomoYRI0bIjmIKf/zxB4qKiu6oYGuaVo2JjCk2NhYAEBwcfEs/x2N7a1JSUjg9hMhOsGAT2ZGVK1eiXbt2POFOJ1FRUQBwWwWbU0puXmxsLGrWrIkmTZrc1Pfz2N4ejmAT2Q8WbCI7ERMTg4MHD2LChAksGDqJioqCh4cH2rRpIzuKocXExEBRFDg780+OKGVlZUhLS+MINpGd4KMdkZ1YsWIFXFxc8MQTT8iOYhq///47OnbsCFdXV9lRDC0mJgZt27aVHcPQkpKSYLPZ0LRpU9lRiAgs2ER2oby8HF9//TUGDRqEoKAg2XFMwWq14vDhw+jVq5fsKIZWWFiIhIQEvkogWHx8PACeSEpkL1iwiezAtm3bkJKSgvHjx8uOYhqqqqKoqAh333237CiGpqoqAHAEWzCu1EJkX1iwiezAJ598ggYNGuCRRx6RHcU0jh49CgC45557JCcxtsqt6FmwxYqPj4eLiws3TCKyEyzYRJIdO3YMv/76K2bMmIEaNWrIjmMa0dHRaNy4MRo1aiQ7iqGdPn0aLi4uaNWqlewohhYfH4/GjRvzfAIiO8GCTSTZwoULUbNmTfz973+XHcVUjh49yukhOoiOjkbbtm25u6Bg8fHxnB5CZEdYsIkkSk5Oxpo1azBx4kT4+fnJjmMacXFxSE1NRe/evWVHMbzo6Gh07txZdgzDY8Emsi8s2EQSzZ07FwDw/PPPS05iLtu3bwcA9O/fX3ISY0tPT0dKSgoLtmDFxcVITU1lwSayIyzYRJLExcVh6dKlePrpp/mHUWfbt29HQEAAl44TLDo6GgAQGhoqOYmxnT9/HgBXECGyJyzYRJK89dZbcHNzw7/+9S9ht6FpmrDrdlQ2mw07duxAr169uGOmYJUrtbBgi3Xu3DkAQIsWLW7p5/j4QCQOCzaRBEeOHMGqVavw7LPPon79+kJug+Xx6k6cOIGMjAxuMKODo0ePomnTpjy/QLDY2FgAuOlXZPjYQCQeCzaRzsrLy/HMM88gKCgIs2fPlh3HdLZs2QKA61/r4fDhw+jatavsGIYXExODoKAgPpEhsiMs2EQ6W7p0KY4cOYIFCxagdu3asuOYzqZNm9C1a1cEBgbKjmJo6enpiIuLw1133SU7iuHFxsbyfAIiO8OCTaSjpKQkvPrqq+jTpw/Gjh0rO47pZGZmIjIyEkOGDJEdxfAiIyMBgAVbME3TEBMTw4JNZGdYsIl0YrPZMGHCBFgsFixdupTzICXYsmULNE3jlvQ6iIyMhJubG7p06SI7iqFlZGQgJyeHW9ET2Rld9lRVFKUNgJcB9ANQH0AxgGMAlqqq+rUeGYhkW7JkCX799Vd89tln3DZaknXr1qFJkybo3LkzVFWVHcfQIiMjERoaCk9PT9lRDC0mJgbAzZ/gSET6ED6CrSjKEABHAUwAEAQgFkAJgPsArFQUZbWiKBzKI0M7dOgQXnzxRQwaNAhTpkyRHceUsrOzERERgdGjR8PZmS/eiWSxWHD48GFOD9HBra4gQkT6EPpXRlGUIACrAXgACANQV1XVTqqq1gMwDEA+gLEAnhOZg0imzMxMjBw5EvXr18fKlSs5NUSS9evXw2q14vHHH5cdxfAOHTqE4uJi9O3bV3YUwzt58iS8vLzQuHFj2VGI6E9ED+NMBuANIArAVFVViyq/oKrqRgCvXf7wn4JzEElRXFyMYcOGIT09HT/88APq1q0rO5Jpffvtt2jdujW37dbBzp074eTkhD59+siOYnjR0dHo1KkTX5UhsjOifyP7Xn67XlVV21W+vuny22aKonABTzKU8vJyjBs3Dvv27cPKlSu5HrBEqamp2LlzJx5//HG+gqCDHTt2IDQ0FHXq1JEdxdBsNhuio6O5UyaRHRJdsF8HMB7Axmt8vdaf3tflhEsiPWiahhkzZmDDhg1YuHAhHnvsMdmRTG3VqlWw2WycHqKDkpISREZGol+/frKjGF5cXBwKCgr4qgyRHRJaalVVPQDgwHW+ZejltxkAMkVmIdKLpmmYOXMmPv30U7z44ouYOXOm7EimpmkavvjiC9xzzz0ICQmRHcfw9u7di9LSUhZsHURHRwMAR7CJ7JC0SVuKotRDxdJ9ALBGVVVNVhai6mKz2TB16lQsWrQIzz//PObPny87kunt3LkTZ86c4eotOtm8eTPc3d15gqMOoqOj4eLignbt2smOQkT/Q8q0DEVRaqFi2ogvKkau593oZywWS9V6n3oqKSmRcrukvzu9r0tLSzF79mxs2bIFTz/9NCZNmlS1hJYMeXl50n5v7MmHH34IHx8ftG/f/opjcaf3d3p6OoCKZdK41vP/t3HjRnTv3h0JCQm3fR0ZGRkAKtZ4ro4580Z9HP/tt9/QokULxMfH39LPJSYmAgDi4+MN93/XqPc1/ZW939e6F2xFUbwAhAPoCaAcwBOqqqbd6Odq1KghZaeqmJgY7pBlEndyX2dkZGDo0KGIjIzEBx98gJdeekn6yXQ+Pj7Sfm/sRVpaGrZv347p06f/ZZ7qnf5uBwYGAqhYf7hmzZp3lNMozp49i/Pnz+OFF164o2MbEBAAAGjbtm21/B4Z8XFc0zScOnUKQ4cOveV/25kzZwAAzZs3N9xxMeJ9TVdn7/e1rgVbUZQAVKwc0gOADcAEVVUj9MxAVN0OHTqExx57DGlpafj+++8xYsQI2ZHosv/85z8oKyvDtGnTZEcxhfDwcADAoEGDJCcxvjNnziArK4ub+RDZKd3mYCuK0gJAJCrKdRkqRq65TTo5LE3T8Mknn6B3794AgD179rBc25H8/Hx8+umnGD58OFq3bi07jimsW7cOoaGhaNGihewohnfgQMX6Ab169ZKchIiuRpeCrShKRwD7ALQEUARgqKqq3+hx20QipKam4tFHH8U///lPDBo0CEePHkX37t1lx6I/CQsLQ25uLl566SXZUUzhwoULiIyM5JKUOomMjISPjw9XxiGyU8ILtqIorQH8AqAegBwAD6iqukX07RKJoGkuLUQEAAAgAElEQVQaVqxYgZCQEEREROCTTz7Bhg0b4OfHfZLsSWlpKRYuXIg+ffqgZ8+esuOYwvfffw8ALNg6OXDgAHr06MEdHInslNDfTEVRaqLihMZAVKwW0ldV1UiRt0kkSmxsLAYOHIjx48ejXbt2OHbsGGbOnCn9ZEb6qy+++AJJSUmYPXu27CimsXbtWnTt2hUtW7aUHcXwLl26hOPHj+Puu++WHYWIrkH0U9/ZABRUnNA4SlXV44Jvj6jaZWVlYcaMGejQoQMOHDiARYsWYffu3VAURXY0uorCwkLMnTsXffv2xQMPPCA7jimcP3++6mRfEm/Pnj2w2Wy4//77ZUchomsQtoqIoijuAKZf/rAIwNwbFJKRqqqmispDdKsKCwvx2Wef4d1330VeXh7+/ve/Y86cOVXLs5F9WrRoEdLT07Fhwwa+uqCTFStWwMnJCaNHj5YdxRR27NgBDw8PnuBIZMdELtPXAUDty+97AbjnBt/vITAL0U2rXH1iwYIFyMzMxEMPPYQFCxagffv2sqPRDaSlpeH999/H4MGD+fK5TsrLy7Fs2TIMGDAATZs2lR3HFHbu3Im7774bHh78s0lkr4QVbFVVjwDg8BE5jPT0dHzzzTdYsmQJsrOzMXDgQLz++ussag7ktddeQ1FRET766CPZUUwjIiICiYmJWLhwoewoppCVlYVjx45h7ty5sqMQ0XVI2SqdyF5omobIyEgsWrQI33//PcrLyzFkyBDMnj0bPXr0kB2PbsHBgwfx1Vdf4aWXXuL8eB2FhYUhICAAQ4YMkR3FFH7++WcAQP/+/SUnIaLrYcEmU8rIyMCaNWuwfPlyREdHo3bt2hg3bhzeeOMNbpLhgGw2G5577jnUr18fr7/+uuw4ppGcnIzw8HA8//zzqFGjhuw4prB582b4+/tz3X0iO8eCTaZRUlKCrVu3YuXKldi0aRPKysrQtWtXfPbZZ3jiiSeQmJjIcu2gPv/8cxw+fBirVq2Ct7e37DimsWjRImiahqlTp8qOYgrl5eXYtm0bBg0aBBcXF9lxiOg6WLDJ0AoLC7Ft2zZ8//332LRpEwoKChAUFIR//OMfeOqpp3jiogHEx8fj5ZdfxoABAzB27FjZcUwjPz8f//3vfzFixAg+MdXJoUOHkJWVhcGDB8uOQkQ3wIJNhpOcnIxt27Zhy5Yt2Lp1K4qLi+Hv748xY8ZgxIgR6N+/P1xd+V/fCGw2GyZOnAhnZ2eEhYVxWT4dffnll8jNzcULL7wgO4ppbNiwAa6urhgwYIDsKER0A2wZ5PAsFgv279+PrVu3YuvWrThx4gQAoGHDhpgwYQJGjhyJe++9l6XagD7//HPs2rULYWFhaNKkiew4pmG1WrFw4UL07t2bW9HrRNM0fPfdd3jwwQfh5+cnOw4R3QAbBzmc4uJiHDx4ELt378bu3bsRGRmJkpISuLq64t5778X8+fPx8MMPo127dhzRNLCTJ0/ihRdewEMPPYRJkybJjmMqX331FRISErBkyRLZUUzj4MGDSEhIwJw5c2RHIaKbwIJNdi89PR2HDx/GgQMHsHv3bhw8eBAWiwVOTk7o1KkTpkyZgr59+6J///48wc0kCgsL8dhjj6F27dpVuwiSPkpLSzF37lz07NkTgwYNkh3HNNauXYsaNWrg0UcflR2FiG4CCzbZlYKCAkRFReHQoUNVl4SEBACAi4sLunbtihkzZqBPnz7o3bs3fH19JScmGZ599lnExsbil19+QVBQkOw4phIWFobExEQsW7aMT2x0YrPZsG7dOgwcOBC1a9e+8Q8QkXQs2CSFpmlITEzE8ePHr7ioqgqbzQYAaNasGXr27InnnnsOPXr0QJcuXVCrVi3JyUm2lStXYvny5XjjjTe42YbO8vLyMHfuXNx777144IEHZMcxjb179yI5ORnz58+XHYWIbhILNgmlaRrS09OhqipiYmJw4sSJqjJ96dKlqu9r1qwZOnbsiFGjRqFnz57o3r07AgICJCYnexQVFYWpU6eiT58+eOONN2THMZ13330XaWlpCA8P5+i1jsLCwlC7dm1ODyFyICzYVC2Ki4tx9uxZqKr6l8ufi7S3tzc6duyIsWPHomPHjujYsSPat28PHx8fienJEaSmpmLo0KEICAjAd999x402dHbmzBksXLgQ48eP5y6COsrOzsa6deswefJk1KxZU3YcIrpJLNh0UzRNQ2pqKuLj4/9yiYuLw4ULF6BpWtX3N2rUCIqiYNy4cVAUperStGlTjnzRLSspKcGwYcOQnZ2Nffv2ITAwUHYk03nhhRfg7u6OefPmyY5iKl9//TVKS0vx9NNPy45CRLeABZsAVKxrm5KSgqSkpKrLnwv0+fPnUVJScsXP1KtXD82bN8c999yDCRMmVJXo4OBgzpWmaqNpGqZMmYIDBw7ghx9+QGhoqOxIprNlyxaEh4fjgw8+QL169WTHMQ1N0/DFF1+gZ8+e6NSpk+w4RHQLWLBNoLCwEBcvXkRSUhKSk5OvKNGVl7S0tCtGoAHA19cXzZs3R0hICAYPHozmzZtXXZo1awZPT09J/yIyk/nz52PlypWYM2cOhg8fLjuO6eTl5WHKlCkICQnBzJkzZccxlR07duD06dNYtmyZ7ChEdItYsB1UaWkp0tLSkJqaitTU1Kr3r/a5goKCv/y8r68vGjVqhEaNGiE0NLTq/YYNG1a9zyXwSLavv/4ar776Kh5//HG8/vrrsuOY0iuvvIKUlBR8//33cHd3lx3HVObPn4969eph7NixsqMQ0S1iwbYDmqYhLy8PmZmZV1yysrL+8rn09HSkpaUhNzf3qtfl5+eHevXqoV69eujWrRvq1auHoKAg1K9f/4oS7eXlpfO/kmT431clHElERAQmTpyIfv36Yfny5Zy7L8GuXbvw+eef4/nnn+eW6DqLjo7Gzz//jHnz5sHDw0PIbTjy4wORvWPBrkY2mw2XLl1Cbm4ucnJy/vI2JyfnqqU5KysLZWVlV71OFxcX+Pv7V106dOiABx98EEFBQVVFuvL9wMBAjjBRFUcupEeOHMGIESPQrl07bNiwgf+vJbh06RImTpyIli1b4p133pEdx3Q+/PBDeHl5YerUqdV+3Y782EDkKFiwr8Nms2H//v2Iioq6amH+37d5eXnXHRFwcXFB3bp1q8pycHAw7r777isK9J+/7u/vDx8fHz4YkqmcPXsWgwYNQkBAALZu3colHCXQNA1Tp07FhQsX8Ntvv3F5OJ2pqoq1a9fiH//4B6fqETkoFuzr2LRpEyZPnnzF52rWrAk/Pz/4+vrCz88PjRo1QocOHa743LXeenl5sSwTXUdaWhoGDhwIm82GiIgI1K9fX3YkU1qxYgW+/fZbvPvuu7jrrrtkxzGdN998Ex4eHnj55ZdlRyGi28SCfR1DhgzBxo0b0bZtW/j6+sLX1xc1atSQHYvIkPLz8zF48GBcvHgRO3bsQHBwsOxIpqSqKp599lncf//9eOWVV2THMZ3o6GisXbsWs2fP5nrvRA6MBfs6nJycEBwczD/0RIJZLBaMHDkS0dHR+PHHH3lCnSRFRUUYPXo0PDw88PXXX3O3TAlmz54NX19fvPjii7KjENEdYMEmIqlsNhsmTZqEn3/+GcuWLcPgwYNlRzIlTdMwefJkHD9+HJs3b0bDhg1lRzKdzZs3Y8uWLZg/fz7nXhM5OBZsIpLqlVdewapVqzB37lxMnDhRdhzT+vjjj/HNN9/gvffew8MPPyw7jumUlJRg5syZUBSFG/oQGQALNhFJ8/HHH2PBggWYPn06Zs2aJTuOaf366694+eWXMXLkSLz66quy45jSxx9/jHPnziEiIoLn+hAZgLPsAERkTqtXr8YLL7yAkSNH4t///jdX2JEkPj4eo0ePRkhICL766iveDxJcuHAB7777LoYNG4YBAwbIjkNE1YAFm4h09/PPP2P8+PHo27cvT6aTKCcnB4MHD4bNZsOGDRu4w6sElXPfnZyc8PHHH8uOQ0TVhFNEiEhXv//+O0aMGIGQkBBs3LhR2DbQdH0WiwXDhw/H2bNn8csvv6BVq1ayI5lSWFgYfvnlF3z66ado1qyZ7DhEVE1YsIlINxcuXMDgwYNRt25dbN26FbVr15YdyZQqR0137dqFVatWoU+fPrIjmVJCQgJeeOEF9OvXD1OmTJEdh4iqEQs2EekiPz8fQ4YMQUlJCXbu3IkGDRrIjmRab775Jr7++mvMnTsX48aNkx3HlMrLyzF+/HgAwLJly+DszBmbREbCgk1EwpWXl2PcuHE4deoUtmzZgrZt28qOZFpffvkl3nnnHUyePJkrt0j07rvvYteuXfjyyy85NYTIgFiwiUi4V199FeHh4Vi8eDFXSZDol19+wZQpUzBgwAB8+umnXDFEkl27dmHOnDl44oknqkaxichY+JoUEQm1fPnyqrWup0+fLjuOaZ04caLq5NJ169bBzc1NdiRTSktLw7hx49CqVSs+ySEyMI5gE5Ewx44dwzPPPIN+/frhk08+kR3HtJKTkzFo0CD4+Phg8+bN8PHxkR3JlCwWC0aOHImcnBxs3rwZ3t7esiMRkSAs2EQkxKVLlzBy5EjUqVMH33zzDVxd+XAjQ35+PgYPHozc3Fzs3bsXjRo1kh3JtGbOnIm9e/dizZo1CA0NlR2HiATiXzwiqnaapmHixImIj4/H7t27ERgYKDuSKVmtVowaNQonT57E5s2b0alTJ9mRTOuLL77A559/jpdffhljxoyRHYeIBGPBJqJq95///Afr16/HRx99hHvuuUd2HFPSNA3Tpk1DREQEli5dioceekh2JNPatWsXnn32WQwcOBDvvfee7DhEpAOe5EhE1erUqVN45ZVXMGTIEPzzn/+UHce05s2bh6VLl2L27NmYNGmS7Dimdfz4cQwdOhStW7fGmjVr4OLiIjsSEemABZuIqo3FYsH//d//wcfHB0uXLuUKCZKsWbMGs2fPxrhx4/DOO+/IjmNaCQkJePjhh+Ht7Y1t27bBz89PdiQi0gmniBBRtZk7dy6OHj2KDRs2cN61JLt378aECRPQt29fLFu2jE9yJMnKysLAgQNRWFiIvXv3onHjxrIjEZGOWLCJqFpER0fjvffew1NPPYVHH31UdhzdaJomO0KVc+fOYdiwYWjZsiXWr18Pd3d32ZHuiKZpDvkEoaioCEOGDEF8fDx+/vlntG/fXnYkItIZCzYR3TGbzYZnnnkGdevWxcKFC2XH0YW9Fb/8/HwMHToUTk5O2LRpk0NPR7C3Y3srysrKMGbMGBw4cADr1q3DfffdJzsSEUnAgk1Ed2zp0qU4cOAAVq5c6dDFzlHZbDY8+eSTiI2NRUREBFq0aCE7kilpmobp06fjp59+wuLFizFixAjZkYhIEhZsIroj6enpePXVV9GnTx888cQTsuOY0ttvv42NGzfik08+Qf/+/WXHMa23334bX3zxBWbNmoXp06fLjkNEEnEVESK6I//617+Qn5+PTz/91KFf2ndUGzduxJw5czB+/HjMmDFDdhzTCgsLw1tvvYXx48dj7ty5suMQkWQs2ER022JiYrBs2TJMmzYNISEhsuOYTnx8PMaPH49u3brhs88+4xMcSX766SdMnToVDz/8ML744gveD0TEgk1Et2/WrFmoVasW/vWvf8mOYjoWiwWPP/44NE3D2rVr4eHhITuSKe3fvx+jR49G165dsW7dOri5ucmORER2gHOwiei2REZGYuPGjZg7dy4CAgJkxzGd2bNn49ChQ/juu+94UqMkMTExGDJkCBo3bozNmzejVq1asiMROazy8nLk5eUhNzf3L5ecnJwrPg4JCcHQoUNlR74uFmwiui1vv/02AgIC8I9//EN2FNPZunUrFixYgGeeeQajRo2SHceULl68iIcffhhubm7Ytm0bn2QSXWa1WpGTk4OsrCxkZWUhOzu76v3//Tg7O7uqNOfl5V33ep2cnFC7dm34+vo6xCtFLNhEdMt+//13bNu2DfPmzeOonc6ys7MxceJEdOjQAR9//LHsOKaUn5+PwYMHIzMzE7t37+YrCGRYmqahsLAQ6enpSEtLQ3p6+hWXjIyMv5Tm6xVlV1dX1K1bF3Xq1EHdunXRrFkz1KlTB76+vje8eHt7w9n5/89sjomJ0eMQ3DYWbCK6Ze+99x58fX0xbdo02VFMZ8aMGcjMzMTWrVs571oCq9WKUaNG4fjx4wgPD0fXrl1lRyK6Zfn5+UhJSUFycjIuXrz4l+L85zJdXFx81euoXbs2AgICULduXQQEBKBNmzZXlOerve/t7W2ak4BZsInolsTGxmL9+vV4/fXX4ePjIzuOqfz4449YvXo13nzzTYSGhsqOYzqapmHq1KmIiIhAWFgYHn74YdmRiK5QWlqKlJSUKy7Jycl/eb+goOAvP+vm5obAwEAEBgYiKCgIbdu2rfq48nOV7wcEBMDd3V3Cv9BxsGAT0S1ZvHgx3N3d8dxzz8mOYiq5ubmYMmUKQkNDMWvWLNlxTOmdd97Bl19+iddffx2TJ0+WHYdMKDc3FwkJCde8pKen/+Vn3N3d0aBBAzRo0AChoaEYNGhQ1ccNGzZE/fr1ERQUhNq1a5tmdFkPLNhEdNPy8vKwYsUKjB49mid16eyNN95ARkYGtmzZgho1asiOYzpr167Fm2++iSeffBJz5syRHYcMqqysDAkJCThz5gzOnj2LM2fOIC4urqpA/+/8Zg8PDzRp0gRNmzZFp06d0KRJEzRq1OiKAu3n58fiLAELNhHdtK+//hoFBQXcBlpnx44dw5IlSzB16lR06dJFdhzTiYqKwoQJE9C7d2+EhYWxrNAdsdlsuHDhAmJiYqpKdOXb8+fPo6ysrOp7a9WqhZYtW6J58+bo27dvVZmuvAQGBvL/o51iwSaim6JpGpYsWYJu3bqhR48esuOYhqZpePbZZ1GnTh1uwS1BWloahg4dCn9/f/zwww989YBums1mQ0JCAk6dOoXTp09XvY2JiUFhYWHV93l7e6NVq1bo0qULHnvsMbRq1QqtW7dGq1atEBQUxALtoFiwieimHDhwoGprdNLP6tWrsXfvXixduhR+fn6y45iKxWLBiBEjkJWVhX379iEwMFB2JLJTRUVFOHHiBKKionD06FEcPXoUp06dumIFjgYNGqBdu3aYPHkyQkJCEBISguDgYAQEBLBEGxALNpGBaZpWbde1atUqeHp6YuTIkdV2nXR9paWlmD17Nrp06YIJEybIjmM6s2bNwr59+/Dtt9+ic+fOsuNUu+p8fDCTkpISREVF4eDBg1WFOiYmBjabDQDg5+eHzp07Y8qUKWjXrl1Vmfb19ZWcnPTEgk1kUNU5ImKxWLB27VoMHTqUS/Pp6L///S8uXLiApUuXXrHBAom3adMmfPTRR5g+fTpGjx4tO0614mjprUlMTERkZGTVJSoqClarFUDFqHSXLl0wfPhwdOnSBZ07d0aTJk14jIkFm4huLCIiAllZWXjiiSdkRzGN/Px8zJ07F/369cMDDzwgO46pJCYm4qmnnkJoaCgWLFggOw7p7MKFC9ixYwd27NiBnTt3IikpCQDg6emJbt264Z///Cfuuusu9OrVC/Xq1ZOcluwVCzYR3dCqVavg7++PAQMGyI5iGosXL0ZGRgbee+89jobpqKysDGPGjIHFYsF3333H3TJNIDs7Gz///HNVqT537hwAICAgAPfffz969+6Nu+66C506dYKbm5vktOQoWLCJ6LqKiooQHh6OCRMm8I+LToqLi/HJJ59gwIAB6Nmzp+w4pjJ//nzs27cPq1evRuvWrWXHIQE0TUNsbCw2bdqETZs2Yd++fSgvL4ePjw/69u2LGTNmoF+/fmjXrh2f3NJtY8Emouvavn07iouL8eijj8qOYhpffvkl0tPT8dprr8mOYionT57EW2+9hcceewxjx46VHYeqkaZpOHLkCL799lts3LgRcXFxAIBOnTrhtddew+DBg9GtWze4urIWUfXg/yQiuq7w8HB4e3ujT58+sqOYgtVqxYcffohevXrxmOvIarVi/Pjx8PX1xeLFi2XHoWpy8uRJfPPNN/j2228RFxcHNzc3PPDAA3jxxRfxyCOPoHHjxrIjkkGxYBPRNdlsNmzatAkPPfQQN9jQyYYNG5CQkIB///vffHlaRwsWLMDvv/+OdevWISAgQHYcugO5ublYvXo1wsLCcOzYMTg7O6N///6YPXs2hg0bxvXkSRcs2ER0TVFRUbh48SKGDBkiO4ppfPbZZ2jWrBkeeeQR2VFMIyEhAe+88w6GDx/Odd4dlKZpOHDgAObPn4+IiAgUFxejc+fOWLRoEUaNGoWgoCDZEclkWLCJ6JrCw8Ph7OyMQYMGyY5iCqdPn8auXbvw/vvvw8XFRXYc03jxxRcBAAsXLpSchG5VWVkZNmzYgI8++ggHDx5ErVq18OSTT+Lpp59G165dZccjE2PBJqJr2rJlC3r16gV/f3/ZUUzhs88+Q40aNTBx4kTZUUxj+/bt+P777/H222+jSZMmsuPQTSotLcWyZcvw4Ycf4vz582jVqhWWLFmCHj16oFu3brLjEYFbgxHRVeXm5iIqKgoPPvig7CimUFhYiBUrVuCxxx7jHGCdWK1WzJgxA82bN8dLL70kOw7dhLKyMixbtgzBwcGYPn066tevjw0bNiA2NhbTpk1DrVq1ZEckAsARbCK6hj179sBms6Ffv36yo5jCTz/9hPz8fEyaNEl2FNNYsWIFTp8+jfXr13NDGTunaRo2bNiAV155BWfPnkWPHj2wdOlSPPDAAzwZmOyS8IKtKIofgDcBDANQD0AGgAgAb6uqmiD69ono9uzYsQMeHh7c6EQnq1evRqNGjXDffffJjmIKJSUlmDNnDnr27Mk13u3cqVOnMHPmTGzfvh3t27fHTz/9hEceeYTFmuya0Ckil8v1PgAzAfgBOAGgJoCJAKIVReko8vaJ6Pbt3LkTvXv3hru7u+wohpeRkYGIiAiMHTsWzs6cuaeH//73v0hKSsK7777LomanSkpK8Oqrr6JTp06IiorC4sWLcfToUQwZMoT3Gdk90Y/kYQDaAtgCoKGqqt0ANACwHIAvgG8VReGp8kR2JicnB8ePH0ffvn1lRzGF7777DmVlZRg3bpzsKKZQXFyM9957D/369UP//v1lx6GrOHToELp06YIPPvgATz31FP744w9Mnz6dOy2SwxBWsBVFaQNgOIACAP+nqmo+AKiqWgJgMoAYVJTvYaIyENHtOXjwIADg7rvvlpzEHNasWYP27dujY0e+qKeH5cuXIz09HW+88YbsKPQ/ysvL8dZbb+Guu+5Cfn4+tm7dimXLlnElI3I4IkewnwDgBCBcVdXsP39BVdVyAF9d/nC0wAxEdBsiIyPh7OyM7t27y45ieKmpqdi/fz9Gj+ZDoR7Ky8uxYMEC9OzZk/Pd7UxmZiYGDRqEOXPmYNy4cTh58iQGDhwoOxbRbRH5WkvlmVH7r/H1A5ff3iswAxHdhsjISHTo0AFeXl6yoxje5s2bAYC7Zerkhx9+QFxcHD788EPO47Ujhw4dwogRI5CRkYGwsDBMmjSJ9w85NJEj2K0uv42/xtcrVxAJUhSFf8WJ7ITNZsPBgwfRq1cv2VFMITw8HI0bN+b0EJ189NFHaN26NYYOHSo7Cl32448/om/fvnBzc8P+/fsxefJklmtyeE6apgm5YkVRCgDUAnCXqqoHrvJ1LwD5lz9srqrq+etd3/Dhw7X169dXe87r+eH3JCzfo6JmzZq63i7JUVRUZKj7OiYmBvn5+ejRo8ct/VxRUSEOHz6CNm3aICgoSFA6+e70/k5MTERcXBzuvbc3nJ1v71xtm82Gffv2oV69emjduvVtZzGahIQEnD9/Hn369KmW66u8rwsKCvD777+jVatWaNiwYbVctyPKysrCyZMn0aVLF3h7e0vNkpKSgjNnzsDb2xsdOnSAm5vbHV2f0R7H6eoe69YYITXz0bZtW91vW1GU3y8v2nFdIqeIeF5+W3yNr//5857X+J4qFosFMTExdxzqVqRczIfNZkNRUZGut0tyGO2+Lisvg6bd+r8pK6vilAlXV1dDHY//daf3t8VqAVDxB/12C/alS5dgs9ng5eVl6GN9q6xWK4CKJ3sVp/Lcmcr7+sKFC3B2doaPj7epj3dpaSmAimXwXFzkLeSVmpqKxMRE+Pr6omXLlrBarVX3/e0y2uM4XV3KxRS0aOimey+8JZqmCbkEBwdbgoODteDg4E7X+Lrb5a9rwcHBbW50fcOGDdNkOH36tJTbJf0Z7b4eM2aM1rp161v+ueeff17z8PDQrFargFT2407v7/nz52sAtIKCgtu+jmeeeUarVauWVlxcfEdZjGbOnDkaAK28vLxaru/06dNaXl6e5uXlpY0fP75artORhYeHawC0w4cPS8uwYMECDYA2evToan2sMdrjOF2brPs6ODj4iHYTPVjkHOzCy2+vtf/sn3evuNYoNxHpLDo6Gh06dOB6szqIiIhA//79uU23Dr755hsUFBRgypQpsqOY3pIlS/Diiy9i1KhRWLVqFR9ryJBEFuysy2/rXOPrdf/0fobAHER0kzRNw9GjRxEaGio7iuElJCQgLi4O/fr1kx3FFJYvX4727dujZ8+eN/5mEiY8PBwzZszA3/72N6xevZrlmgxLZMGOvfy22TW+3vTy24uqqnLCFJEdSExMRE5ODgu2Dnbu3AkALNg6SE5ORmRkJMaNG8fVKSSKiorC448/ji5dumDNmjV3fEIjkT0TWbCPXH57rbW+Kj9/UGAGIroF0dHRAMCCrYOdO3fC398f7dq1kx3F8LZt2wYA3MxHorS0NAwZMgT+/v4IDw9HrVq1ZEciEkpkwa5cU2+YoihXTBNRFMUFwPjLH64SmIGIbsHx48cBAB06dJCcxNg0TZhEd9QAACAASURBVMOOHTtw//33w9lZ5MMwAcCWLVvQs2dPNG/eXHYUUyovL8cTTzyB7Oxs/PTTT6hXr57sSETCCXtkV1X1OIAtALwBfK8oSl0AUBTFA8BSAG0BqAA2iMpARLcmJiYGTZo0kb42rtElJiYiKSmJW3XrQFVVxMTE4PHHH5cdxbTmzZuHX3/9FYsWLUKnTp1kxyHSheizC6YA2AvgfgAXFEWJAdACgB+ASwCGq6pqE5yBiG5SbGyslIX7zSYyMhIAcNddd0lOYnzr1q2Dk5MTRo0aJTuKKUVFReGtt97C2LFjMWnSJNlxiHQj9LVJVVWTAHQF8B9UrBTSEUAZgG8AdFdV9bTI2yeim2ez2RAbG4s2bdrIjmJ4kZGR8PT05PboOti8eTPat29v6p0bZbFYLJgwYQICAwOxePFinmBKpiJ8fRxVVbMAzLx8ISI7lZSUhKKiIo5g6+DAgQPo1q0bV1EQLDMzEwcPHsS0adNkRzGl999/H8ePH8ePP/4IPz8/2XGIdMWza4gIQMX0EAAcwRaspKQEUVFRnB6ig4iICGiaxrnuEiQkJGDevHkYPXo0/va3v8mOQ6Q7FmwiAlBxgiPAgi1adHQ0rFYrevW61gqmVF22bt2KgIAALoUowcsvvwwnJyd8+OGHsqMQScGCTUQAKkaw/fz8EBgYKDuKoUVFRQEAunbtKjmJsZWXl2Pbtm0YOHAgl0LU2W+//YbvvvsOr7zyCho3biw7DpEUfNQhIgDAH3/8geDgYJ6IJFh0dDTq1KnD4iHY4cOHkZWVhUGDBsmOYiqapmHWrFlo2LAhXnrpJdlxiKRhwSYiAEB8fDxatGghO4bhHT16FKGhoXwiI9j27dsBAA888IDkJOayc+dO7N27F7NmzULNmjVlxyGShgWbiFBWVoYLFy5wpzvBrFYrTpw4gc6dO8uOYnh79uxB+/bt4e/vLzuKaWiahrfeegsNGjTAxIkTZcchkooFm4iQlJSE8vJyFmzBYmNjUVpayoItmNVqxb59+9CnTx/ZUUxlz549+O233/Dqq6/Cw8NDdhwiqViwiQjx8fEAwIItWHR0NAAgNDRUchJjO3r0KAoLC7k8n84WLlwIf39/TJ48WXYUIulYsImIBVsnJ0+ehJubGxRFkR3F0Hbv3g0ALNg6On/+PMLDw/H3v/8dnp6esuMQSceCTUSIi4uDs7MzV7YQLCYmBsHBwXB1Fb6JrqlFRkaiZcuWqFevnuwopvH5558DAKZOnSo5CZF9YMEmMjBN027q++Lj49G4cWNu3S1YTEwMt6LXwaFDh9CzZ0/ZMezezT4+3EhxcTHCwsLw6KOP8kk60WUs2EQGdSvLwMXHx3N6iGAlJSWIi4tjwRYsOTkZycnJ6NGjh+wodqu6l4jcuHEjsrOzMW3atGq9XiJHxoJNRDh//jwLtmBnzpyBzWZjwRbs8OHDAMCCraPVq1ejUaNGuP/++2VHIbIbLNhEJme1WpGamsqXdgWLiYkBABZswQ4dOgRXV1eu1KKTzMxMREREYMyYMdySnuhP+NtAZHKpqanQNA0NGzaUHcXQYmJi4OTkxBVEBDt06BA6duzIlSx0sm7dOpSVlWHcuHGyoxDZFRZsIpNLSUkBABZswWJjY9G0aVMWP4E0TUN0dDS6dOkiO4pprFmzBiEhIejYsaPsKER2hQWbyOSSk5MBAA0aNJCcxNjOnTuH1q1by45haBcvXkRWVhY6deokO4opJCUlYe/evRg7dmy1nzhJ5OhYsIlMjiPY+oiPj0eLFi1kxzC048ePAwBHU3WyadMmAMDw4cMlJyGyPyzYRCaXnJwMNzc3+Pv7y45iWHl5ecjMzORKLYJVFuwOHTpITmIO4eHhaNmyJdq0aSM7CpHdYcEmMrmUlBTUr1+fKwAIVLkVPUewxTpx4gQaN24MPz8/2VEMr7CwENu3b8eQIUM4PYToKvgXlcjkkpOTOf9aMBZsfRw/fpzTQ3SyY8cOlJaW4pFHHpEdhcgusWATmVxKSgrnXwsWFxcHAJwiIpDVakVMTAynh+jk119/haenJ3r37i07CpFdYsEmMjmOYIsXFxeH2rVrc+qCQHFxcbBarQgJCZEdxRS2b9+Oe++9F+7u7rKjENklFmwiEysoKEBeXh5HsAWrXEGEc1XF+eOPPwAAwcHBkpMYX2pqKk6dOoX+/fvLjkJkt1iwiUzs4sWLAID69etLTmJscXFxnH8tmKqqwP9r786j4zrLPI//KpJs2bJjO47thCxWNt44KwmBkLBMJt1kwmGb0GEyadKsyQycwDA9nKYbaLqbBgLpGYaZPmcYhgAJSydNM8OWQE/INCTDEBOakIUO8mMIXvCSWHIsySpJpZJU88d7r12WS7KW+1Sp6n4/5/iUpapSvXpf3arffe573ysCdj384Ac/kCRdffXVDW4JsHgRsIEc27dvnyRpw4YNDW5J66pUKtq5c6dOP/30RjelpW3dulXr1q1jGk4dPPTQQ1q1apUuueSSRjcFWLQI2ECO9fX1SZLWrVvX4Ja0rv7+fg0PD+vUU09tdFNa2tatW6le18nDDz+sK664Qm1tbY1uCrBoEbCBHOvt7ZUkLjKzAJVKZcb700vRM8997o7Vt9XMTCEEx9ZAkgYGBvTUU0/pyiuvbHRTgEWNgA3kWBqwqWDP3WxPWCRgz91cTwYdHBzUM888QwW7Dh555BFVKhVdccUVjW4KsKgRsIEc6+vrU1dXl5YtW9boprQsAra/X/3qV5I4wbEeNm/erEKhoBe/+MWNbgqwqBGwgRzr7e1leoizNGCz1rifp59+WpJ01llnNbglrW/z5s268MILdfzxxze6KcCiRsAGcqyvr4/pIc527dqldevWcUEORzt27JAkdXd3N7YhLa5SqejRRx/VZZdd1uimAIseARvIMSrY/nbv3s30EGc7duzQ6tWrqao627Nnj/r6+lieD5gFAjaQY729vVSwnRGw/W3fvp3qdR089thjkkTABmaBgA3kGFNE/BGw/e3YsUMbN25sdDNa3uOPP65CoaCLLrqo0U0BFj0CNpBTIyMjKhaLTBFxVCqV1Nvby0VmHFUqFW3fvp2AXQePPfaYzj77bK1cubLRTQEWPQI2kFOsge1v7969klhBxNOBAwc0NDREwK6Dxx9/XC94wQsa3QygKRCwgZxKL5NOBdvPvn37JEkbNmxocEtaFyuI1MfBgwf1m9/8RhdffHGjmwI0BQI2kFNUsP2lfbx+/foGt6R1bd++XZKoYDszM0nSeeed1+CWAM2BgA3kVFrBJmD7SSvY9LGftIJNwPbV09MjSdq0aVODWwI0BwI2kFPPPfecJGnNmjUNbknr4iiBv127dqmzs1Nr165tdFNa2pYtW9Te3s7VMoFZImADOTUwMCBJWr16dYNb0rp6e3vV2dmprq6uRjelZe3du1cnn3yyCoVCo5vS0rZs2aKzzjpLHR0djW4K0BQI2EBO9ff3a9myZXxgOtq3b5/Wr19P+HOUBmz46unpYXoIMAcEbCCnBgYGqF4740qZ/gjY/srlsn7961/r3HPPbXRTgKZBwAZyamBgQKtWrWp0M1oaAdsfAdvftm3bVC6XCdjAHBCwgZzq7++ngu2MgO1reHhYAwMDBGxnW7dulSQ9//nPb3BLgOZBwAZyigq2v3QONnykV8okYPvatm2bJOnMM89scEuA5kHABnKKgO2rWCxqZGSECrYjAnZ9bNu2TcuXL2dnEZgDAjaQU0wR8cUa2P4I2PWxbds2dXd3sxoOMAcEbCCnqGD7ImD7I2DXx7Zt23TGGWc0uhlAUyFgAzlUKpU0OjpKwHaUXiadw+p+9u7dq/b2dp144omNbkrLqlQqBGxgHgjYQAurVCo1v89VHP319fVJEuHP0d69e7VhwwYddxwfZfMx3ftDtQMHDmhwcFDd3d3+DQJaCO9KQIuaab5kGrCpYPvp7++XJJ1wwgkNbknreuaZZ5geMg9zmUudriBCBRuYGwI2kENp+KOC7Sft4+OPP77BLWldrDPuj4ANzA8BG8ghKtj++vv7dfzxx6utra3RTWlZ+/fvZwqOs+3bt0siYANzRcAGcoiA7a+/v5/+dbZ//36tXbu20c1oabt27dKKFSs42gXMEQEbyCGmiPhjnXFfpVJJQ0NDBGxne/bs0fOe97xGNwNoOgRsIIeoYPsjYPvav3+/JFZp8bZnzx6dcsopjW4G0HQI2EAODQwMqFAoaOXKlY1uSssiYPtKAzYVbF+7d++mgg3MAwEbyKH0BDzWD/ZDwPZFwPZXqVSYIgLME5+uQA4NDg6yfJwzArYvAra/5557TmNjYwRsYB4I2EAODQ0NacWKFY1uRsuanJzUwMAAAdsRV8r0t3v3bkliDjYwDwRsIIeGh4fV1dXV6Ga0rIMHD6pSqRCwHVHB9rdnzx5JooINzAMBG8ihYrFIwHbEMoj+9u/fr+XLl6uzs7PRTWlZBGxg/gjYQA4RsH0RsP1xFUd/BGxg/gjYQA4Vi0UtX7680c1oWQRsf319fUwPcbZ7926tXbtWS5cubXRTgKZDwAZyiAq2LwK2Py6T7m/Pnj06+eSTG90MoCkRsIEcImD7ImD7Y4qIv3379umkk05qdDOApkTABnKIgO2LgO2PCra/vr4+dmKAeSJgAzkzMTGh0dFRArajNGBzMR8flUpFBw4c0Jo1axrdlJZGwAbmr937BUII50p6v6SrJZ0saUTSE5I+b2Zf8X59AEcaGRmRJAK2o/7+fq1cuVLt7e5vsblULBZVqVS0atWqRjelZZXLZfX39xOwgXlyffcPIbxW0t9J6pQ0KmmLpA2SXiHpFSGEayXdZGYVz3YAOKxYLEoiYHviUvS+BgYGJHGEwFN6IZ9169Y1uCVAc3KbIhJC2CDpbxTD9R2S1prZxWZ2kqTrJB2U9PuS3uPVBgBHI2D741L0vgYHByURsD1xKXpgYTznYN8saaWkn0t6p5kNp3eY2bckfSD58g8d2wBgijRgsw62n2KxSMB2RMD2R8AGFsYzYF+V3H7DzCZr3H9fctsdQuBMFaBOqGD7Gxoaon8dpQGbOdh+0oDNFBFgfjznYH9Y0lcl/Wya+6s/fTgTCKgTAra/YrGo9evXN7oZLYs52P56e3slUcEG5sst2JrZTyT9ZIaHvD657ZXU59UOAEciYPtjnXFfTBHxl1awWWscmJ+GrIMdQjhJcek+SbqbVUSA+iFg+2OKiC8Ctr++vj6tWrVKS5YsaXRTgKZU94AdQuiS9C1JqxUr15+odxuAPBsejucbEwD9cJKjL6aI+OMiM8DCzGqKSAjhLklvmeXPvN/Mrp3m56yQdK+kyyVNKK6B/exsfujY2Jh6enpm2YTsjI6ONuR1UX+tNtYDAwM1t5unn35akrRr1y4NDQ01ommLwkLH+9ln41vXli1bjtpZOXjwYMv9PdVTOv+3p6dHbW1tR92/bds2LVu2TFu3bp3Vz2MsjrRz505JsR9XrlxZ8zHbt29XV1dX0/UbY50fi32sZzsHuyypNMvHjtX6ZghhneLKIS+WNCnpbWZ2/yx/ppYsWaJNmzbN9uGZ6enpacjrov5abazTw7tTf6f0A/XSSy/V0qVLG9G0RWGh471hwwZJ0rnnnntEtbpcLqtcLmvjxo0t9fdUT+nKFZs2baoZsDs6OrR69epZ92+rbdsLtW3bNknSGWecMW2/jIyM6LTTTmu6fmOs82Oxj/WsAraZ3SLplvm+SAjhTEnfl3SWpHFJbzaze+b78wDMX7FYVFtbG3MrnTDH3d/AwABL9Dnr6+vTxRdf3OhmAE3LfXm8EMJFku6XdJKkYUlvNLPveb8ugNqKxaKWL1+uQqHQ6Ka0pDRgMwfbD5ei97d//36dcMIJjW4G0LRcA3YI4RxJD0haL+mApFeb2WbP1wQwM5aQ80UF2x8B21e5XNbw8LBWr17d6KYATcttFZEQwnLFExrXK64WchXhGmg8Arav9MRR+tjPwMAAAdtRukoL03CA+fNcpu9DkoLiCY1vNLMnHV8LwCwRsH0xRcTf4OAg4c8RARtYOJcpIiGEpZJuTb4clvSxEMJMT7nezJ7xaAuAIw0PDxOwHVHB9scUEV8EbGDhvOZgXygp3TJXSHrpMR7f6dQOAFNQwfZFBdvX5OSkDh48SMB2RMAGFs4lYJvZzySxRAGwCBWLxUPrDCN7VLB9DQ0NqVKpEP4cEbCBhav7pdIBNNbw8LCWL1/e6Ga0LCrYvgYHByVxmXRPacBmFRFg/gjYQM6MjIxo2bJljW5Gy2KZPl8EbH9UsIGFI2ADOTM6OqrOTk578DI0NKRCocBOjJODBw9K4giBpzRgsxMDzB8BG2hhlUrlqO8RsH2lJ5FypUwfw8PDkjhCkIVa7w9SDNjLly9XR0dHnVsEtA4CNtCipgt4BGxfQ0NDhD9HacDmPIL5O9bO38DAANNDgAUiYAM5Mj4+rvHxcQK2o2KxyPQFRwRsfwRsYOEI2ECOlEolSSJgO2KdcV8EbH8EbGDhCNhAjoyOjkoiYHsaGhqigu2IgO2PgA0sHAEbyBECtj8q2L4I2P4I2MDCEbCBHCFg++MkR19pwGYZRD8EbGDhCNhAjhCw/Y2MjFBddTQ8PKyOjg6WkHNEwAYWjoAN5AgB2x9XyvTFDoyvcrms4eFhAjawQARsIEdGRkYkEbA9sc64r+HhYQK2Iy6TDmSDgA3kSFrBpsLqhwq2LwK2rzRgr169usEtAZobARvIEaaI+KpUKlSwnRGwfR08eFCStHLlyga3BGhuBGwgRwjY2atUKof+Xy6XNTk5SQU7I9V9myJg+yoWi5LESjjAAhGwgRwhYGenUCgc9T3muGejVt+mCNi+0mUQCdjAwhCwgRwhYPtijrs/ArYvKthANgjYQI4QsH1RwfZHwPaVBmz6GFgYAjaQI2nAXrp0aYNb0pqoYPsjYPtiigiQDQI2kCNUsH1RwfZHwPZFBRvIBgEbyJGxsTFJUnt7e4Nb0prSgE0F2w8B2xdzsIFsELCBHBkbG9PSpUtnXKUB88cRAl+VSoWA7Wx4eFjt7e1asmRJo5sCNDUCNpAjpVKJD05HVLB9lUolVSoVArajYrFI/wIZIGADOTI2NkbAdkQF21d6Ah4B0M/w8DDTQ4AMELCBHCmVSqwg4ogKti8Ctj8q2EA2CNhAjlDB9sUyfb4I2P6KxSIVbCADBGwgR9KTHOGDZfp8EbD9MUUEyAYBG8gRTnL0xRxsX2nA5giBH6aIANkgYAM5QgXbV6lUksSVMr1QwfbHFBEgGwRsIEeoYPsqlUoqFApcyMcJAdsf64wD2SBgAznCSY6+0lVauJCPDwK2PyrYQDYI2ECOMEXEF0cIfBGw/RGwgWwQsIEcIQD6YgfGFyc5+mOKCJANAjaQIwRAX1zIxxfrjPsql8sql8tUsIEMELCBHKGC7YuA7YtVWnwxBQfIDgEbyBEq2L4I2L5GR0dVKBTU0dHR6Ka0pGKxKElUsIEMELCBHKGC7YuA7YtVWnwRsIHsELCBHGGZPl8EbF/0ry+miADZIWADOcIUEV8cIfBVKpW4DL0jKthAdgjYQI4QAH2xA+NrdHSU/nWUBmwq2MDCEbCBnJiYmNDExAQBxRFTGHzRv77SKSJUsIGFI2ADOTE2NiZJVLAdEQB9MUXEF1NEgOwQsIGcIGD7I2D7YoqIL05yBLJDwAZaWKVSOfR/ArY/ArYv+jdb1e8PkjQyMiKJgA1kgYANtKipawUTsP0RAH0xRSQb060jnl6Knr9hYOEI2EBOlMtlSeIqeI4I2L7oX19pwGYnBlg4AjaQE+Pj45Kk9vb2BrekdbEMoi/mYPsaHR1VW1sb7xFABgjYQE5QwfZVqVRYB9sZU0R80b9AdgjYQE4QsH2Nj4+rUqkQsB0xRcTX6OgoARvICAEbyAkCtq9SqSSJE8Q8MUXEF/0LZIeADeQEAdsXq7T4YwqDLyrYQHYI2EBOELB9EbD9MUXEFwEbyA4BG8gJArYv+tdXpVIhYDvjCAGQHQI2kBMEQF/0r69yuaxKpUIAdEQFG8gOARvICQKgL6aI+KJ//RGwgewQsIGc4EIzvtiB8ZUGbKaI+GEVESA7BGwgJwiAvuhfX2nApn/9UMEGskPABnKCAOgr7V+mMPhgiog/AjaQHQI2kBMEbF9UWH2xA+OPVUSA7BCwgZwgYPuif32xA+OPCjaQHQI2kBMEQF/0ry+miPjjJEcgOwRsICcIgL4IgL6YIuKPCjaQHQI2kBMEbF/0ry+miPianJxUuVwmYAMZIWADOUEA9EX/+uIIga9SqSSJ/gWyQsAGcoIA6IspDL4I2L64kA+QLQI2kBPplRzb2toa3JLWxBQGX+zA+CJgA9kiYAM5US6X1d7erkKh0OimtCSOEPhiB8YXRwiAbBGwgZwol8uEE0cEbF8EQF/0L5AtAjaQEwRsXwQUX/SvL05yBLLVXu8XDCFcKelHkn5rZt31fn0grwjYvqhg+6J/fbEDA2SrrhXsEMJSSV+o9+sCIGB7IwD6IgD6on+BbNU76H5E0rl1fk0AImB7SwN2e3vdDwzmAgHQF6uIANmqW8AOIVwq6X2SRur1mgAOI2D7GhsbU0dHB6u0OCFg+6J/gWzVJWCHEDok3SmpIulj9XhNAEciYPuif30xBccXJzkC2apXBfuDki6SdLukJ+v0mgCqjI+PM33BEQHbF+tg+6KCDWTLPWCHEC5QDNhbRPUaaBgCoK+xsTHCiaOxsTG1t7fruOM4R94DARvIlus7VQihTdIXJXVIutnMSp6vB2B6VLB9sQPjK53jDh+c5Ahky7sU8D5JL5L0GTP7sfNrAZgBAdsXAdtXuVymuuqICjaQrVl92oYQ7pL0lln+zPvN7NoQwjmKy/L9VtIH5te8w8bGxtTT07PQHzNno6OjDXld1F+rjfXAwIDK5fKh32lwcPCIr/NuoeP97LPPSpLMTF1dXdq/f78k0b8Z6O3tlRT7Mt0pfPbZZ9XW1jav/m21bXuhdu7cKUnavn27Vq1adej/krRjxw6Njo42qmkLxljnx2If69mWs8qSZju9YyyEUFC8oEynpHeZ2cH5NK7akiVLtGnTpoX+mDnr6elpyOui/lptrFetWqWOjo5Dv9PSpUvV2dnZUr/jQix0vDds2CBJCiFo5cqV6uzsVFdXF/2bgXXr1kmSNm3adChgd3V1zfvvt9W27YXasWOHJKm7u/tQv6xdu1aSdP755x/q/2bEWOfHYh/rWQVsM7tF0i2z/aEhhHdLermke8zsu/NsG4AMTUxMqLOzs9HNaFnj4+NMEXHESaS+mCICZMtrQub1ye2NIYQbp3nMxhBCJfn/GWa23aktAMQcbG/0ry/muPviJEcgW16fBr+Y4WevkXSe4pSTnyXfa94JX0CTGB8fV1tbW6Ob0bII2L44ydEX64wD2XL5NDCz90x3XwjhNZLulfSMmb3M4/UBHG1iYoIA6IiA7Yv+9TU2Nqa2tjZ2woGMsGI/0MIqlcqh/1PB9kUA9FUul+nfjFW/P5RKJY4QABkiYAMtqlAoHPE1AdAX/euLk0izM/W9QeIkUiBrBGwgJ5gi4ouA7Yv+9TU2NsYJjkCG6v5uZWb3STp69xmAK6aI+CIA+mKKiC8q2EC2qGADOUEF2xcB2xdTRHyNjY3Rv0CGCNhAThAAfXGEwBd/v75YBhHIFgEbyAkCoI90JQYCYPaqV7lgiogv+hfIFgEbyAmmiGSLVVr81FrlgikivrhSJpAtAjaQE1SwfRGwfdG/vgjYQLYI2EBOEFB8cYTAF1MYfPH+AGSLgA3kBAHQFwHFF1NEfFHBBrJFwAZygikivgjYvuhfXwRsIFsEbCAnqGD7IgD6YoqIL44QANkiYAM5MDk5qUqlQkBxRMD2RQD0xQ4MkC0CNpAD4+PjksQUEUcEbF/0ry+miADZImADOTAxMSFJBBRHBEBfVFh9cYQAyBYBG8gBKti+KpUKAdsZAdAXOzBAtgjYQA6kAZsPUB+Tk5OS6F9P7MD4YooIkC0CNpADTBHxxQ6Mr4mJCU7SdUbABrJFwAZygCkivgjYvtL+JQD64QgBkC0CNpADVLB9EbB90b/+qGAD2SJgAzlAQPFF//oql8uS6F9PBGwgWwRsIAeYIuKLgO2LKSL+WKUFyBYBG8gBpoj4ImD7on/9sUwfkC0CNpADBBRf9K8vpoj4mpyc1OTkJBVsIEMEbCAHmCLii4DtiykivuhfIHsEbCAHmCLii4Dti/71xRECIHsEbCAHqGD7IgD6IgD6SvuXCjaQHQI2kANUsH0RsH0xhcEXARvIHgEbyAECoC/61xf964v+BbJHwAZygCkivggovpgi4osKNpA9AjaQA0wR8UXA9sUUEV8EbCB7BGwgBwiAvuhfX/SvL3ZggOwRsIEcYIqILwKgL6aI+KJ/gewRsIEcYIqILwK2L/rXF1NEgOwRsIEcoILti/71le4g0r8+mCICZI+ADeTA5OSkJAKKFwKgL/rXF1NEgOwRsIEcSAPKccexyXsgAPqif30xRQTIHp+2QA5QwfZFAPTFOQS+mOMOZI+ADeQAFWxfBGxfzHH3xQ4MkD0+bYEcoILti/71xQ6ML3ZggOwRsIEWVqlUJFHB9kb/+iJg+5j6/kAFG8gOnwZAiyoUCof+T4XVFwHQF/2brer3BokKNuCBgA3kABVWXwRAX/SvL/oXyB6ftkAOUMH2RUDxRf/6YhURIHsEbCAHqGD7YgfGFwHbF/0LZI9PWyAHCIC+2IHxRQD0xUmOQPb4NABygADoiwDoiwDoi5McgezxaQvkABVsXwRsXwRAX+zAANkjYAM5QAXbFwHbF/3rix0YIHt82gI5QAXbFwHQF/3riwo2kD0CNpADVLB9sQPji4Dtiwo2kD0+bYEcIAD6YgfGFwHbFxVsvo035wAAEdlJREFUIHt8GgA5QAD0RQD0NTExoUKhcNQlvpENKthA9vi0BXIgrWATsH2wA+NrYmKC8OeICjaQPT4NgByYmJgg/DlK+5cKqw8Cti8q2ED2+MQFcmBycpIPT0f0r6+JiQmqq444AgNkj60JyAEq2L7oX1/j4+PswDhK+5cjMEB2+EQAcoAKqy+mMPiif33Rv0D2CNhADlBh9UVA8UX/+mIKDpA9PnGBHKCC7YsA6Iv+9cUUHCB7BGwgB6hg+2IHxhcB2xcVbCB7fOICOUAA9MUOjC8Cti8q2ED2+EQAcoAA6IsA6Iv+9UUFG8gen7hADlDB9kUA9EX/+qKCDWSPgA3kABVsXwRAX1RYfdG/QPb4xAVygAq2LwK2LyqsvuhfIHsEbCAHqGD7YgfGFzswvqhgA9njExfIAQKgL3ZgfBGwfVHBBrLHJwKQAwRAXwRAX/SvL/oXyB6fuEAOUMH2RUDxRf/6Gh8fZ4oIkDECNpADVLB9EQB90b++6F8ge3ziAjlABdsX/euLAOiLkxyB7BGwgRyggu2L/vVFwPbFSY5A9tx3WUMIx0m6WdJbJZ0vqUPSFkl3SPqsmVW82wDkHRVWXwRAXxMTE+rs7Gx0M1oWFWwge64llxBCp6TvSvofki6XtEtSr6RLJH1G0j0hhIJnGwBQYfVGwPZFhdUX/Qtkz/sT93ZJ10r6raRLzOx8M9so6bWSipJukPQm5zYAuUcF2xcB2xf964sKNpA9t4AdQjhT0q2SxiW9ysyeTO8zs/skfSr58u1ebQAQUcH2xQ6MLwK2LyrYQPY8d1lvlNQm6S4ze6rG/XdKKkna6dgGACIAeqMC6IuA7Yu/XyB7nlvU7yS33651p5ltl3Sb4+sDSFDB9kUA9EX/+qKCDWTPM2BfkNz2hBBWSXqbpFdIWiHpl5I+Z2a/dHx9AInJyUkqVI4IgL7oX19UsIHsuZS0ktVD1iVfnibpKUmflnSdpFdKeq+kJ0IIt3i8PoAjUcH2RQD0Rf/6ooINZM/rE3dl1f/vkTQi6VWSlikG7k8rVs8/G0K42qkNQO5VKnGZeeZg+6lUKvSvk/Tvl4Dtg/4F/BTSDWwmIYS7JL1llj/zfkm36PDJiyOSzjezbVN+5lcVl+h7xMxeMos29EraMcs2AAAAAFnbaGbrjvWg2U66Kiuu+DEbY4qhOvWVqeE6cZtiwL48hLDezPbN9ENn88sAAAAAjTargG1mtyhWpWclhLBEUkVSQdKT0zxsq2Jw75DULWnGgA0AAAA0A5c52GY2JqlW1bpaJfknxaANAAAAND3PZQV+mtxeNs39GyUtkTQpabtjOwAAAIC68QzYX0turw8hnFLj/luT24fM7IBjOwAAAIC68QzY35G0WfHCMveFEM5K7wgh3CDp3cmXH3dsAwAAAFBXs1qmb76SyvU/SAqSJhSv4LhC0hnJQz5sZh9za8A8hRDWSPpzxQvjnCSpV3H5wb80M5YKXKRCCOdKer+kqyWdrLiazROSPm9mX5nmOW9W3Nm7UHGlnMcl/Wcz+84Mr3O+4t/HP1dc832HpL+VdLuZDWf2C2HWQghXSvqRpN+aWXeN+9skvUfSOySdI2lI0iOSPmlmPzrGz/1TSVcqTmn7taQ7Jf21mU1k/GughhDCcZJulvRWSecrnhi/RdIdkj5rZkd9iLFdN6cQwqmK29urFN/DByU9rDgGP57heYx3E0i25c2SzjKzE2d4XF3GM4SwTNKfSLpRcdpyv6QHJX3MzH4x199vKteALUkhhOWS3ifpBklnSipKelTSp83sftcXn4ckXP9Y0iZJBxVXOzlT0hrFzv9nZjbdyihokBDCayX9naROSaOK47Yh+SdJd0u6qfrDOITwSUl/rHiy7VPJc89O7v4zM/tojdd5oaSHJHVJekbSHkkXKIavJyS93MwOZv37YXohhKWKb8DnStoxNWCHEAqS/kbxTXRC0i8krVW86NWkpHeY2V01fu6rJX1bUpvim/UBSRcpHvn7vqRXm9m4yy8FSYeuCvxNSdcqjtUWxSLN6clDvibpRrbr5pcEpAclnajD7+GnKG6rk5JuMbMv1nge490kQgi3SfqApP3TBex6jWey2t0Dkl6hGOL/SdKpiplhVNLrzOyBhfy+7tdONrNhM/uomV1gZsvNbJ2ZXbsYw3XiDsVw/T1Jp5jZZZKeJ+kuSasl/W1SDcMiEULYoBigOhXHb62ZXWxmJykehTgo6fcVK5jpc16juBE/J+lyM7vQzM6R9C8VN7aPhBBeOuV1OhUDV5ekjyj+fbxQcQfsUUkXS/qU5++Kmj6iGK6nc6tiuP6NpE1mdoliteJdikuJfjaEcEb1E0IIJymGt+MUP9i7k+ddoLhC0jWKR0vg63bFcP1bSZeY2flmtlHSaxWLNTcoXk9BEtt1k/uyYrj+oaTTzexixbDzccXt8L/X2E4Z7yYQQiiEEP5cMVzP9Lh6judtiuH6HyV1J1nvFMW/t05Jd4cQjp/nryypDgG7mSRTDN6gePj4D9I9HjMbVTxE2aMYvq9rWCNRy82Kh4R+Lumd1YeDzOxbOrxR/2HVcz6Y3P6Jmf1j1eO/rbiRFnT0m8GbFTfAh83sL8xsMnnObkm/p7jc5NtCCCdn9YthZiGESxWPkI1Mc3+7pD9Kvnynmf1KksysYmafVdwhW5r8jGr/TvEN+x4z+3z6TTPr0eFA976keg4HIYQzFXeOxiW9qvrIoZndp8Mfmm+vehrbdRMKIZwn6VLFquVNZtYrSWY2YWZ/qhiClijuKFdjvBe5pFjxTUl/MYuH12U8QwgnSHqn4pGRN5vZM8lz0r+37yvu7P2buf22RyJgH+kmxQG818yeq74jmW95Z/LlDfVuGGZ0VXL7jXTjmuK+5LY7hLAmhHC2pCsUN7y7azz+C8ntNSGE1VXff3Ny+6WpT0jm5j+gePGmN8yt+ZiPEEKH4jZZkTTduRxXKU4n2DvN4b40PF+fzA9M/UFyW2usNyueT3KCpN+de8sxSzcqTs/5qpk9VeP+OyV9SNIXJYntuqmlK43tN7M9Ne7/eXKbTg1ivJtACOEaxak+r1ecujFtBbvO4/kGxQLKI2a2pcZrpZ8LC8p6BOwjXZ7cPjzN/T9Jbl9eh7Zg9j6seALUt6a5v6vq/+06PM5Pmllx6oPNbJ/idIIOSS+RDp0k98LkIfx9LA4fVJwTfbumv2Lssbbpn0saUzwU/Xzp0MnZpyb3b57meYy1v99Jbr9d604z225mt5nZV5NvsV03r13J7YnTLOt7fnK7s+p7jPfid57iORNfUZxe95MZHlvP8Zxt1rs0hNA1zWOOaVaXSs+RdCL9dFehTFcQ2RBCWGFmQ3VoE47BzH6imTfc1ye3vZL6dOxxluJYn1n12FMV52VJ018YKf37OHua+5GREMIFigF7i2L1+pXTPHTGsTaziRDCbsWVjc5Ofl76nL4ZTnRirP1dkNz2hBBWSXqb4pzJFYpHED5nZr+sejzbdZMys54QwsOKq/V8KYTwr82sLzlB+Y8kvUxx6uaXq57GeC9+P5V0qZk9LkkhhJkeW8/xPNZr7Vacmtau+NnwTzO0aVpUsI+0LrndP8391dNGpl1iBotHMv8rPRnt7mS1gWONs3R4rNNxTp9TmmHHaupz4CCpWnxRsZJxs5mVZnj4QsZ6Ls9BhpITl9JxOE1xNYFPK57/8kpJ75X0RAjhlqqnsV03t+sk/R/FIxc7QwiPS9qreISqR9I1Zrar6vGM9yJnZg+n4XoW6jmeM75WMtV0oMbz5oSAfaRlyW3NE6amfH/ZNI/BIpEc2vmW4uovfZI+kdx1rHGuvm/ZlNu5PAc+3ifpRZI+YzOsjZtgrJvTyqr/36PY369S7O/TFMN2u+IKMFcnj2Osm1tJ8UjkiGJfp6uISHHZtak70ox3a6nneNbl74CAfaRjXTiiur98FxDHgoQQViie3Hi54rjeZGbPJnfXOhFyqnSs03Gez3OQsRDCOYpnk/9Wx1jyKcFYN6fOqv93KVYv/7eZjZrZLjP7D4pLcx6nuNyWxFg3reSktQcVLzLzI0mXKK7uc6akv1asav/fEMIVVU9jvFtLPcezLn8HzME+UlGx2tk5zf3VS3LNtOeDBgohrFMM1y9W3JDeNmXd9fRQ0nTjLB0e63Sc5/McZCiZj/kFxTF41ywvBMFYN6fqfv2KmdWaK3mb4pKJl4cQ1ouxbmbvl/QCxYtAvcbMysn3t0l6bwihpDgX+78pLucnMd6tpp7jWZe/AyrYR0rn45wwzf1rq/7f69wWzEOydu5mxXA9rli5nnqZ9GONs3R4rNNxTp/TmVyddDbPQbZuVTwT/B4z++4sn7OQsZ7Lc5CtQR2uHE23QsxWxSW9JKlbbNfN7PeS2/9YFa6rfULxSOQlyXJuEuPdauo5njO+VnKez6oaz5sTAvaR0vUQu6e5f2Nyu9emubY9GieEcJHiZe7PkjQs6fVmdk+Nhx5rnKXDY/1rSUpOrhmact+Mz0Hmrk9ubwwhVKr/Sbo3uW9j1fe7dYyxTi5Eky4Llo5b+pz1M7xpM9aOzGxMM68mIMUAnobwstium1nax7XWJJaZHdDhoDP1sd2z+LmM9+JXz/E81mudojjDY1JxacB5IWAf6WfJ7UumuT/9/iN1aAvmIJmb+4CkkyQdkPS7Zva9aR6ejvPFyWoFU3/WesWleSYVryCWejS55e+jMX6huANV61+6XFup6nujOjzWl6u2SxVXI9mvw2/afTq8tNN0z2Os/f00ub1smvs3Kl7db1JxiS626+Y1mNzWvHrilFVl0scy3q2lnuM526z3ZHIl73khYB/pG8ntdcmlNA9JDhm8Nfnyq8KikVQZ75W0XnG1kKuSq+3VZGbbJT2mOMfqphoPuTm5/V5SOUmlfx/vqNGGjYpX9RuT9PU5/gqYBTN7j5m9rNY/SX+cPOyZqu8/o3jCVJ9iZbvWVRfTsb5nylVAZxrrKyVtUqyofT+L3w01fS25vX6ai4/cmtw+ZGYH2K6b2oPJ7VFjkLhJ8aqeA5KekHgfbzV1Hs97FaeQvizUXpw7/VkLynoE7Cpm9qSk7ykuEfU/QwhrpUN7z59X/FA1Sd9sWCNRy4ckBcU92zcm43gs6ZJ9nwohXJV+M4TwOkl/pnjo+fYpz7lT8XKvLw0hfCqZXqAQwvMk/S/FSuiXp7nULxrAzCYk/VXy5Z3JNCJJUgjh3yq+aZcUl32r9l8UT255Uwjh3ycnWCqEcK4Ov+n+VzPjRCg/31E8n2KFpPtCCGeld4QQbpD07uTLj1c9h+26Od2uGHheF0L4q+qr54UQrpf0qeTLTybTh1KMd2upy3ia2X5Jn1PMwF9PgrhCCG0hhI9KukbxqOYdC/llCpUKK9FUCyGcKun/KR5+HFZc4P5MSWsU956vnHL1MDRQCGGppGcVT0gYUlLdmMH1SWVTIYQvSHp78v1fKm6I5yRff8jMbpv65BDCKxX3fpdK2qd4id/zk68fk/TyWpd5ha8QwmsUx2WHmXVPua8jue9fKO6E/UJxez49ecibzOzuGj/zLYpv3gXFce6TdKFiJe3vJb02CfBwklSu/0FxB3pCcTtdoXioWJI+bGYfm/IctusmlGxvdyiO15DiSayn6fDUkC8prghVmfI8xrtJJKH5h5L2m1nNC7jUazyTq8P+UHFJyLLi58IpimuvlxWXBn1wnr+qJCrYR0kmzb9Qce3NXkkXKe5Z3yPpRYTrRedCHT7bd4Wklx7jX/XcrpsVN+SfKp7scKpixezGWhuxJJnZA4pzQr+uGLwuVLwIwn9SnJrCm/Iik6xK8BrFC9Q8pRjW1igGt2tqhevkeV9SvDT33yuuxXyepF8prtV7HeHan5ntVpwn/2c6fBn7lZLul3Tt1HCdYLtuQsn29iLFI0QDimPQrnh1x39lZm+dGq4TjHdrqct4mtmA4qpUH1U8ofoCxXM6vqNYSH1wob8IFWwAAAAgQ1SwAQAAgAwRsAEAAIAMEbABAACADBGwAQAAgAwRsAEAAIAMEbABAACADBGwAQAAgAwRsAEAAIAMEbABAACADBGwAQAAgAwRsAEAAIAM/X93hg8xB4t+ugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "lm = [1, 2, 3, 4]\n",
    "M = len(lm)\n",
    "D = np.array(lm)\n",
    "a = np.min(lm)\n",
    "b = np.max(lm)\n",
    "t = np.linspace(-1, 6, 1000)\n",
    "u = 0.5 * np.ones(M)\n",
    "rho = 1\n",
    "def fun(lam):\n",
    "    return 1 + rho * np.sum(u**2/(D - lam))\n",
    "res = [fun(lam) for lam in t]\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(res, 'k')\n",
    "plt.plot(np.zeros_like(t))\n",
    "plt.ylim([-6, 6])\n",
    "plt.tight_layout()\n",
    "plt.yticks(fontsize=24)\n",
    "plt.xticks(fontsize=24)\n",
    "plt.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "- The function has only one root at $[d_i, d_{i+1}]$\n",
    "\n",
    "- We have proved, by the way, the **Cauchy interlacing theorem** (what happens to the eigenvalues under rank-$1$ perturbation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## How to find the root\n",
    "\n",
    "- A Newton method will fail (draw a picture with a tangent line).\n",
    "\n",
    "- Note that Newton method is just approximation of a function $f(\\lambda)$ by a linear function.\n",
    "\n",
    "- Much better approximation is the **hyperbola**:\n",
    "\n",
    "$$f(\\lambda) \\approx c_0 + \\frac{c_1}{d_i - \\lambda} + \\frac{c_2}{d_{i+1} - \\lambda}.$$\n",
    "\n",
    "- To fit the coefficients, we have to evaluate $f(\\lambda)$ and $f'(\\lambda)$ in the particular point.\n",
    "\n",
    "- After that, the approximation can be recovered from solving **quadratic equation**\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Important issues\n",
    "\n",
    "- First, stability: this method was abandoned for a long time due to instability of the computation of the eigenvectors.\n",
    "\n",
    "- In the recursion, we need to compute the eigenvectors of the $D + \\rho uu^*$ matrix.\n",
    "\n",
    "- The exact expression for the eigenvectors is just (let us check!)\n",
    "\n",
    "$$(D - \\alpha_i I)^{-1}u,$$ where $\\alpha_i$ is the computed root.\n",
    "\n",
    "- The reason of instability:\n",
    "    - if $\\alpha_i$ and $\\alpha_{i+1}$ are close, then the corresponding eigenvectors are collinear, but they have to be orthogonal\n",
    "    - if $\\alpha_i$ and $\\alpha_{i+1}$ are close, then they close to $d_i$, therefore matrices $D - \\alpha_i I$ and $D - \\alpha_{i+1} I$ are close to singular"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Loewner theorem\n",
    "\n",
    "- The solution came is to use a strange **Loewner theorem:**\n",
    "\n",
    "If $\\alpha_i$ and $d_i$ satisfy the **interlacing theorem** \n",
    "\n",
    "$$d_n < \\alpha_n < \\ldots < d_{i+1} < \\alpha_{i+1} \\ldots$$\n",
    "\n",
    "Then there exists a vector $\\widehat{u}$ such that $\\alpha_i$ are exact eigenvalues of the matrix\n",
    "\n",
    "$$\\widehat{D} = D + \\widehat{u} \\widehat{u}^*.$$\n",
    "\n",
    "- So, you first compute the eigenvalues, then compute $\\widehat{u}$ and only then the eigenvectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Divide and conquer and the Fast Multipole Method\n",
    "\n",
    "In the computations of divide and conquer we have to evaluate the sums of the form\n",
    "\n",
    "$$f(\\lambda) = 1 + \\rho \\sum_{i=1}^n \\frac{|u_i|^2}{(d_i - \\lambda)},$$\n",
    "\n",
    "and have to do it at least for $n$ points. \n",
    "\n",
    "- The complexity is then $\\mathcal{O}(n^2)$, as for the QR algorithm.\n",
    "\n",
    "- Can we make it $\\mathcal{O}(n \\log n)$? \n",
    "\n",
    "- The answer is yes, but we have to replace the computations by the approximate ones\n",
    "by the help of **Fast Multipole Method**.\n",
    "\n",
    "I will do a short explanation on the whiteboard."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Few more algorithms\n",
    "\n",
    "- Absolutely different approach is based on the **bisection**.\n",
    "\n",
    "- Given a matrix $A$ its inertia is defined as a triple \n",
    "\n",
    "$$(\\nu, \\zeta, \\pi),$$\n",
    "\n",
    "where $\\nu$ is the number of negative, $\\zeta$ - zero and $\\pi$ - positive eigenvalues. \n",
    "\n",
    "- If $X$ is non-singular, then \n",
    "\n",
    "$$Inertia(A) = Inertia(X^* A X)$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Bisection via Gaussian elimination\n",
    "\n",
    "- Given $z$ we can do the Gaussian elimination:\n",
    "\n",
    "$$A - zI = L D L^*,$$\n",
    "\n",
    "and inertia of the diagonal matrix is trivial to compute.\n",
    "\n",
    "- Thus, if we want to find all the eigenvalues in the interval $a$, $b$\n",
    "\n",
    "- Using inertia, we can easily count the number of eigenvalues in an interval. \n",
    "- Illustration: if $Inertia(A)=(5,0,2)$ and after shift $Inertia(A-zI)=(4,0,3)$, $z\\in[a,b]$ then we know that $\\lambda(A)\\in[a,z]$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method\n",
    "\n",
    "- Recall what a Jacobi (Givens rotations) are\n",
    "\n",
    "- In a plane they correspong to a $2 \\times 2$ orthogonal matrix of the form\n",
    "\n",
    "$$\\begin{pmatrix} \\cos \\phi & \\sin \\phi \\\\ -\\sin \\phi & \\cos \\phi \\end{pmatrix},$$\n",
    "\n",
    "and in the $n$-dimensional case we select two variables $i$ and $j$ and rotate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method (cont.)\n",
    "\n",
    "- The idea of the Jacobi method is to minimize sum of squares of off-diagonal elements:\n",
    "\n",
    "$$\\Gamma(A) = \\mathrm{off}( U^* A U), \\quad \\mathrm{off}^2(X) = \\sum_{i \\ne j} \\left|X_{ij}\\right|^2 = \\|X \\|^2_F - \\sum\\limits_{i=1}^n x^2_{ii}.$$\n",
    "\n",
    "by applying succesive Jacobi rotations $U$ to zero off-diagonal elements. \n",
    "\n",
    "- When the \"pivot\" is chosen, it is easy to eliminate it. \n",
    "\n",
    "- The main question is then what is the order of **sweeps** we have to make (i.e. in which order to eliminate).\n",
    "\n",
    "- If we always eliminate the largest off-diagonal elements the method has quadratic convergence.\n",
    "\n",
    "- In practice, a cyclic order (i.e., $(1, 2), (1, 3), \\ldots, (2, 3), \\ldots$) is used."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi method: convergence\n",
    "\n",
    "- To show convergence, we firstly show that \n",
    "\n",
    "$$ \\text{off}(B) < \\text{off}(A), $$\n",
    "\n",
    "where $B = U^*AU$.\n",
    "\n",
    "- In this case we use the unitary invariance of Frobenius norm and denote by $p$ and $q$ the indices that is changed after rotation:\n",
    "\n",
    "$ \\Gamma^2(A) = \\text{off}^2(B) =  \\|B\\|^2_F - \\sum\\limits_{i=1}^n b^2_{ii} = \\| A \\|^2_F - \\sum\\limits_{i \\neq p, q} b^2_{ii} - (b^2_{pp} + b^2_{qq}) = \\| A \\|^2_F - \\sum\\limits_{i \\neq p, q} a^2_{ii} - (a^2_{pp} + 2a^2_{pq} + a^2_{qq}) = \\| A \\|^2_F - \\sum\\limits_{i =1}^n a^2_{ii} - 2a^2_{pq} = \\text{off}^2(A) - 2a^2_{pq} < \\text{off}^2(A)$\n",
    "\n",
    "- We show that the ''size'' of off-diagonal elements decreases after Jacobi rotation.\n",
    "\n",
    "- If we always select the largest off-diagonal element $a_{pq} = \\gamma$ to eliminate (pivot), then we have\n",
    "\n",
    "$$ |a_{ij}| \\leq \\gamma, $$\n",
    "\n",
    "thus  \n",
    "\n",
    "$$ \\text{off}(A)^2 \\leq 2 N \\gamma^2, $$\n",
    "\n",
    "where $2N = n(n-1)$ is the number of off-diagonal elements.\n",
    "\n",
    "- Or rewrite this inequality in the form\n",
    "\n",
    "$$2\\gamma^2 \\geq \\frac{\\text{off}^2(A)}{N}.$$\n",
    "\n",
    "Now we use relations $\\Gamma^2(A) = \\text{off}^2(A) - 2\\gamma^2 \\leq \\text{off}^2(A) - \\dfrac{\\text{off}^2(A)}{N}$ and get\n",
    "\n",
    "$$ \\Gamma(A) \\leq \\sqrt{\\left(1 - \\frac{1}{N}\\right)} \\text{off}(A). $$\n",
    "\n",
    "- Aften $K$ steps we have the factor\n",
    "\n",
    "$$\\left(1 - \\frac{1}{N}\\right)^{\\frac{K}{2}} \\approx e^{-\\frac{1}{2}},$$\n",
    "\n",
    "i.e. linear convergence. However, the convergence is locally quadratic (given without proof here)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Jacobi: summary\n",
    "\n",
    "Jacobi method was the first numerical method for the eigenvalues, proposed in 1846.\n",
    "\n",
    "- Large constant\n",
    "- Very accurate (high relative error for small eigenvalues)\n",
    "- Good parallel capabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary for this part\n",
    "- Many algorithms for the computation of the SEV solution:\n",
    "    - QR\n",
    "    - Divide-and-conquer\n",
    "    - Bisection\n",
    "    - Jacobi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next lecture\n",
    "- We start **sparse and/or structured** NLA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Questions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\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",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}