lectures/lecture-9/lecture-9.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 9: 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 today\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": [
    "## 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": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "DeviceArray([-1.47168069, -0.75157939,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ,\n",
       "              0.        ,  0.        ,  0.        ,  0.        ],            dtype=float64)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGiCAYAAACh/hJSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAASlklEQVR4nO3dXYhU9/nA8Wc1dQzBnSCa3Z1mNZpXCM0WkroIyUVwye5eSCy5MJKLTQgthFIIS5AG6ktoQJJCkbRiLq03Nr2Jl15kqUqpMSRBetWidksUs5tGcGfdkgno+V+U7J+tmmg9u8/M7OcDhzgvOj8Ow3xzZs7M01EURREAsMCWZC8AgMVJgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASNGSAdq3b1/cd999sXz58ujv74+PPvooe0lta/fu3dHR0TFne+SRR7KX1TaOHz8emzdvjlqtFh0dHXH48OE5txdFETt37oyenp648847Y2BgIE6fPp2z2DbxXfv8xRdfvOY5PzQ0lLPYNtdyAXrvvfdidHQ0du3aFZ9++mn09fXF4OBgfPHFF9lLa1uPPvpofP7557Pbn//85+wltY2ZmZno6+uLffv2Xff2t99+O955551499134+TJk3HXXXfF4OBgfPXVVwu80vbxXfs8ImJoaGjOc/7QoUMLuMJFpGgxGzZsKH72s5/NXr5y5UpRq9WKPXv2JK6qfe3atavo6+vLXsaiEBHF+++/P3v56tWrRXd3d/HrX/969rpLly4VlUqlOHToUMIK289/7/OiKIqRkZHi2WefTVnPYtNSR0Bff/11fPLJJzEwMDB73ZIlS2JgYCBOnDiRuLL2dvr06ajVarF+/fp44YUX4rPPPste0qIwPj4eExMTc57v1Wo1+vv7Pd/n2dGjR+Oee+6Jhx9+OF555ZW4ePFi9pLaUksF6Msvv4wrV65EV1fXnOu7urpiYmIiaVXtrb+/Pw4cOBBHjhyJ/fv3x/j4eDz11FMxPT2dvbS2981z2vN9YQ0NDcXBgwdjbGws3nrrrTh27FgMDw/HlStXspfWdu7IXgDNbXh4ePbPjz32WPT398fatWvjj3/8Y7z88suJK4P58fzzz8/++Qc/+EE89thjcf/998fRo0dj06ZNiStrPy11BLRq1apYunRpTE5Ozrl+cnIyuru7k1a1uNx9993x0EMPxZkzZ7KX0va+eU57vudav359rFq1ynN+HrRUgJYtWxaPP/54jI2NzV539erVGBsbi40bNyaubPG4fPlynD17Nnp6erKX0vbWrVsX3d3dc57v9Xo9Tp486fm+gM6fPx8XL170nJ8HLfcW3OjoaIyMjMQTTzwRGzZsiL1798bMzEy89NJL2UtrS6+99lps3rw51q5dGxcuXIhdu3bF0qVLY9u2bdlLawuXL1+e83/W4+PjcerUqVi5cmWsWbMmXn311XjzzTfjwQcfjHXr1sWOHTuiVqvFli1b8hbd4r5tn69cuTLeeOONeO6556K7uzvOnj0b27dvjwceeCAGBwcTV92msk/D+1/89re/LdasWVMsW7as2LBhQ/Hhhx9mL6ltbd26tejp6SmWLVtWfP/73y+2bt1anDlzJntZbeNPf/pTERHXbCMjI0VR/OdU7B07dhRdXV1FpVIpNm3aVPz973/PXXSL+7Z9/u9//7t45plnitWrVxff+973irVr1xY/+clPiomJiexlt6WOoiiKrPgBsHi11GdAALQPAQIghQABkEKAAEghQACkECAAUggQAClaNkCNRiN2794djUYjeymLgv29sOzvhWefL7yW/SJqvV6ParUaU1NT0dnZmb2ctmd/Lyz7e+HZ5wuvZY+AAGhtAgRAiqb7NeyrV6/GhQsXYsWKFdHR0XHD+9Xr9Tn/ZX7Z3wvL/l549nl5iqKI6enpqNVqsWTJjY9zmu4zoPPnz0dvb2/2MgC4TefOnYt77733hrc33RHQihUrspdwXVNTU9lLAGgJ9Xo9ent7v/P1vOkC9G1vu2VyVgzArfmu13MnIQCQYt4CtG/fvrjvvvti+fLl0d/fHx999NF8PRQALWheAvTee+/F6Oho7Nq1Kz799NPo6+uLwcHB+OKLL+bj4QBoQfNyFlx/f3/86Ec/it/97ncR8Z9Tq3t7e+PnP/95/OIXv/jWv/vNt5GbTZOdLAjQtG72VyVKPwL6+uuv45NPPomBgYH/f5AlS2JgYCBOnDhR9sMB0KJKPwvuyy+/jCtXrkRXV9ec67u6uuJvf/vbNfdvNBpzfvzPl8AAFof0s+D27NkT1Wp1dvMlVIDFofQArVq1KpYuXRqTk5Nzrp+cnIzu7u5r7v/666/H1NTU7Hbu3LmylwRAEyo9QMuWLYvHH388xsbGZq+7evVqjI2NxcaNG6+5f6VSic7OzjkbAO1vXn4JYXR0NEZGRuKJJ56IDRs2xN69e2NmZiZeeuml+Xg4AFrQvARo69at8a9//St27twZExMT8cMf/jCOHDlyzYkJACxeTfdr2L4HBNDa0r4HBAA3Q4AASNF04xiaVVljIryVB/AfjoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApGjakdxTU1PR2dl5W/9GWWO0y1Tmmoz3BlqZIyAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSNO1E1DKUOTG0naermqwKZHAEBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiBFW4/kLlNZY6vbebR3hPHewM1zBARACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKE1EXWJkTQ9t5uqrJqtD+HAEBkEKAAEghQACkECAAUggQAClKD9Du3bujo6NjzvbII4+U/TAAtLh5OQ370UcfjQ8++OD/H+QOZ3sDMNe8lOGOO+6I7u7u+finAWgT8/IZ0OnTp6NWq8X69evjhRdeiM8+++yG9200GlGv1+dsALS/0gPU398fBw4ciCNHjsT+/ftjfHw8nnrqqZienr7u/ffs2RPVanV26+3tLXtJADShjmKef/Pk0qVLsXbt2vjNb34TL7/88jW3NxqNaDQas5fr9Xr09vbG1NRUdHZ2zufSWl4z/hRPWfwUD7Suer0e1Wr1O1/H5/3sgLvvvjseeuihOHPmzHVvr1QqUalU5nsZADSZef8e0OXLl+Ps2bPR09Mz3w8FQAspPUCvvfZaHDt2LP75z3/GX/7yl/jxj38cS5cujW3btpX9UAC0sNLfgjt//nxs27YtLl68GKtXr44nn3wyPvzww1i9enXZDwVACys9QH/4wx/K/icBaEN+Cw6AFAIEQAo/0tbCyvquTDN+n6jMNflOETQnR0AApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMBGVUieGtvN0VZNVoVyOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMJKbUpU1trqdR3tHGO8NEY6AAEgiQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASGEiKk2pzImh7Txd1WRVWpkjIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApjOSm7ZU1trqdR3tHGO/NwnMEBEAKAQIghQABkEKAAEhxywE6fvx4bN68OWq1WnR0dMThw4fn3F4URezcuTN6enrizjvvjIGBgTh9+nRZ6wWgTdxygGZmZqKvry/27dt33dvffvvteOedd+Ldd9+NkydPxl133RWDg4Px1Vdf3fZiAWgft3wa9vDwcAwPD1/3tqIoYu/evfHLX/4ynn322YiIOHjwYHR1dcXhw4fj+eefv73VAtA2Sv0MaHx8PCYmJmJgYGD2umq1Gv39/XHixIkyHwqAFlfqF1EnJiYiIqKrq2vO9V1dXbO3/bdGoxGNRmP2cr1eL3NJADSp9LPg9uzZE9VqdXbr7e3NXhIAC6DUAHV3d0dExOTk5JzrJycnZ2/7b6+//npMTU3NbufOnStzSQA0qVIDtG7duuju7o6xsbHZ6+r1epw8eTI2btx43b9TqVSis7NzzgZA+7vlz4AuX74cZ86cmb08Pj4ep06dipUrV8aaNWvi1VdfjTfffDMefPDBWLduXezYsSNqtVps2bKlzHUD0OJuOUAff/xxPP3007OXR0dHIyJiZGQkDhw4ENu3b4+ZmZn46U9/GpcuXYonn3wyjhw5EsuXLy9v1QC0vI6iyX6DvV6vR7VajampKW/H0VSacRxDmZrspYAWdrOv4+lnwQGwOAkQAClMRIWbVOZbVM34dl5Za/JWHjfLERAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFEZyQ4Kyxla382jvCOO9250jIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFKYiAotrMyJoe08XdVk1ebkCAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACiO5gYgob2x1O4/2jjDeu0yOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEhhIipQqjInhrbzdFWTVR0BAZBEgABIIUAApBAgAFIIEAApbjlAx48fj82bN0etVouOjo44fPjwnNtffPHF6OjomLMNDQ2VtV4A2sQtB2hmZib6+vpi3759N7zP0NBQfP7557PboUOHbmuRALSfW/4e0PDwcAwPD3/rfSqVSnR3d//PiwKg/c3LZ0BHjx6Ne+65Jx5++OF45ZVX4uLFize8b6PRiHq9PmcDoP2VHqChoaE4ePBgjI2NxVtvvRXHjh2L4eHhuHLlynXvv2fPnqhWq7Nbb29v2UsCoAl1FLfxexAdHR3x/vvvx5YtW254n3/84x9x//33xwcffBCbNm265vZGoxGNRmP2cr1ej97e3piamorOzs7/dWlAG2jGn+IpSzv/FE+9Xo9qtfqdr+Pzfhr2+vXrY9WqVXHmzJnr3l6pVKKzs3POBkD7m/cAnT9/Pi5evBg9PT3z/VAAtJBbPgvu8uXLc45mxsfH49SpU7Fy5cpYuXJlvPHGG/Hcc89Fd3d3nD17NrZv3x4PPPBADA4OlrpwAFrbLQfo448/jqeffnr28ujoaEREjIyMxP79++Ovf/1r/P73v49Lly5FrVaLZ555Jn71q19FpVIpb9UAtLzbOglhPtzsh1dA+3MSQmtqmpMQAOB6BAiAFEZyA02rrLepmvGtvDLX1Kpv5zkCAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIWJqEDbK3NiaDtPV13oyaqOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMJIb4BaUNba6nUd73yxHQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQwERUgQVmTVSOac7rqzXAEBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQIqmm4j6zZTAer2evBIAbsd3TX1tugBNT09HRERvb2/ySgC4HdPT01GtVm94e0dR5mDyEly9ejUuXLgQK1as+NY55/V6PXp7e+PcuXPR2dm5gCtcnOzvhWV/Lzz7vDxFUcT09HTUarVYsuTGn/Q03RHQkiVL4t57773p+3d2dnqyLCD7e2HZ3wvPPi/Htx35fMNJCACkECAAUrRsgCqVSuzatSsqlUr2UhYF+3th2d8Lzz5feE13EgIAi0PLHgEB0NoECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFL8H4nN3o57QjHaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAGiCAYAAACh/hJSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAASlklEQVR4nO3dXYhU9/nA8Wc1dQzBnSCa3Z1mNZpXCM0WkroIyUVwye5eSCy5MJKLTQgthFIIS5AG6ktoQJJCkbRiLq03Nr2Jl15kqUqpMSRBetWidksUs5tGcGfdkgno+V+U7J+tmmg9u8/M7OcDhzgvOj8Ow3xzZs7M01EURREAsMCWZC8AgMVJgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASNGSAdq3b1/cd999sXz58ujv74+PPvooe0lta/fu3dHR0TFne+SRR7KX1TaOHz8emzdvjlqtFh0dHXH48OE5txdFETt37oyenp648847Y2BgIE6fPp2z2DbxXfv8xRdfvOY5PzQ0lLPYNtdyAXrvvfdidHQ0du3aFZ9++mn09fXF4OBgfPHFF9lLa1uPPvpofP7557Pbn//85+wltY2ZmZno6+uLffv2Xff2t99+O955551499134+TJk3HXXXfF4OBgfPXVVwu80vbxXfs8ImJoaGjOc/7QoUMLuMJFpGgxGzZsKH72s5/NXr5y5UpRq9WKPXv2JK6qfe3atavo6+vLXsaiEBHF+++/P3v56tWrRXd3d/HrX/969rpLly4VlUqlOHToUMIK289/7/OiKIqRkZHi2WefTVnPYtNSR0Bff/11fPLJJzEwMDB73ZIlS2JgYCBOnDiRuLL2dvr06ajVarF+/fp44YUX4rPPPste0qIwPj4eExMTc57v1Wo1+vv7Pd/n2dGjR+Oee+6Jhx9+OF555ZW4ePFi9pLaUksF6Msvv4wrV65EV1fXnOu7urpiYmIiaVXtrb+/Pw4cOBBHjhyJ/fv3x/j4eDz11FMxPT2dvbS2981z2vN9YQ0NDcXBgwdjbGws3nrrrTh27FgMDw/HlStXspfWdu7IXgDNbXh4ePbPjz32WPT398fatWvjj3/8Y7z88suJK4P58fzzz8/++Qc/+EE89thjcf/998fRo0dj06ZNiStrPy11BLRq1apYunRpTE5Ozrl+cnIyuru7k1a1uNx9993x0EMPxZkzZ7KX0va+eU57vudav359rFq1ynN+HrRUgJYtWxaPP/54jI2NzV539erVGBsbi40bNyaubPG4fPlynD17Nnp6erKX0vbWrVsX3d3dc57v9Xo9Tp486fm+gM6fPx8XL170nJ8HLfcW3OjoaIyMjMQTTzwRGzZsiL1798bMzEy89NJL2UtrS6+99lps3rw51q5dGxcuXIhdu3bF0qVLY9u2bdlLawuXL1+e83/W4+PjcerUqVi5cmWsWbMmXn311XjzzTfjwQcfjHXr1sWOHTuiVqvFli1b8hbd4r5tn69cuTLeeOONeO6556K7uzvOnj0b27dvjwceeCAGBwcTV92msk/D+1/89re/LdasWVMsW7as2LBhQ/Hhhx9mL6ltbd26tejp6SmWLVtWfP/73y+2bt1anDlzJntZbeNPf/pTERHXbCMjI0VR/OdU7B07dhRdXV1FpVIpNm3aVPz973/PXXSL+7Z9/u9//7t45plnitWrVxff+973irVr1xY/+clPiomJiexlt6WOoiiKrPgBsHi11GdAALQPAQIghQABkEKAAEghQACkECAAUggQAClaNkCNRiN2794djUYjeymLgv29sOzvhWefL7yW/SJqvV6ParUaU1NT0dnZmb2ctmd/Lyz7e+HZ5wuvZY+AAGhtAgRAiqb7NeyrV6/GhQsXYsWKFdHR0XHD+9Xr9Tn/ZX7Z3wvL/l549nl5iqKI6enpqNVqsWTJjY9zmu4zoPPnz0dvb2/2MgC4TefOnYt77733hrc33RHQihUrspdwXVNTU9lLAGgJ9Xo9ent7v/P1vOkC9G1vu2VyVgzArfmu13MnIQCQYt4CtG/fvrjvvvti+fLl0d/fHx999NF8PRQALWheAvTee+/F6Oho7Nq1Kz799NPo6+uLwcHB+OKLL+bj4QBoQfNyFlx/f3/86Ec/it/97ncR8Z9Tq3t7e+PnP/95/OIXv/jWv/vNt5GbTZOdLAjQtG72VyVKPwL6+uuv45NPPomBgYH/f5AlS2JgYCBOnDhR9sMB0KJKPwvuyy+/jCtXrkRXV9ec67u6uuJvf/vbNfdvNBpzfvzPl8AAFof0s+D27NkT1Wp1dvMlVIDFofQArVq1KpYuXRqTk5Nzrp+cnIzu7u5r7v/666/H1NTU7Hbu3LmylwRAEyo9QMuWLYvHH388xsbGZq+7evVqjI2NxcaNG6+5f6VSic7OzjkbAO1vXn4JYXR0NEZGRuKJJ56IDRs2xN69e2NmZiZeeuml+Xg4AFrQvARo69at8a9//St27twZExMT8cMf/jCOHDlyzYkJACxeTfdr2L4HBNDa0r4HBAA3Q4AASNF04xiaVVljIryVB/AfjoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApGjakdxTU1PR2dl5W/9GWWO0y1Tmmoz3BlqZIyAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSNO1E1DKUOTG0naermqwKZHAEBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiBFW4/kLlNZY6vbebR3hPHewM1zBARACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKE1EXWJkTQ9t5uqrJqtD+HAEBkEKAAEghQACkECAAUggQAClKD9Du3bujo6NjzvbII4+U/TAAtLh5OQ370UcfjQ8++OD/H+QOZ3sDMNe8lOGOO+6I7u7u+finAWgT8/IZ0OnTp6NWq8X69evjhRdeiM8+++yG9200GlGv1+dsALS/0gPU398fBw4ciCNHjsT+/ftjfHw8nnrqqZienr7u/ffs2RPVanV26+3tLXtJADShjmKef/Pk0qVLsXbt2vjNb34TL7/88jW3NxqNaDQas5fr9Xr09vbG1NRUdHZ2zufSWl4z/hRPWfwUD7Suer0e1Wr1O1/H5/3sgLvvvjseeuihOHPmzHVvr1QqUalU5nsZADSZef8e0OXLl+Ps2bPR09Mz3w8FQAspPUCvvfZaHDt2LP75z3/GX/7yl/jxj38cS5cujW3btpX9UAC0sNLfgjt//nxs27YtLl68GKtXr44nn3wyPvzww1i9enXZDwVACys9QH/4wx/K/icBaEN+Cw6AFAIEQAo/0tbCyvquTDN+n6jMNflOETQnR0AApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMBGVUieGtvN0VZNVoVyOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMJKbUpU1trqdR3tHGO8NEY6AAEgiQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASGEiKk2pzImh7Txd1WRVWpkjIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApjOSm7ZU1trqdR3tHGO/NwnMEBEAKAQIghQABkEKAAEhxywE6fvx4bN68OWq1WnR0dMThw4fn3F4URezcuTN6enrizjvvjIGBgTh9+nRZ6wWgTdxygGZmZqKvry/27dt33dvffvvteOedd+Ldd9+NkydPxl133RWDg4Px1Vdf3fZiAWgft3wa9vDwcAwPD1/3tqIoYu/evfHLX/4ynn322YiIOHjwYHR1dcXhw4fj+eefv73VAtA2Sv0MaHx8PCYmJmJgYGD2umq1Gv39/XHixIkyHwqAFlfqF1EnJiYiIqKrq2vO9V1dXbO3/bdGoxGNRmP2cr1eL3NJADSp9LPg9uzZE9VqdXbr7e3NXhIAC6DUAHV3d0dExOTk5JzrJycnZ2/7b6+//npMTU3NbufOnStzSQA0qVIDtG7duuju7o6xsbHZ6+r1epw8eTI2btx43b9TqVSis7NzzgZA+7vlz4AuX74cZ86cmb08Pj4ep06dipUrV8aaNWvi1VdfjTfffDMefPDBWLduXezYsSNqtVps2bKlzHUD0OJuOUAff/xxPP3007OXR0dHIyJiZGQkDhw4ENu3b4+ZmZn46U9/GpcuXYonn3wyjhw5EsuXLy9v1QC0vI6iyX6DvV6vR7VajampKW/H0VSacRxDmZrspYAWdrOv4+lnwQGwOAkQAClMRIWbVOZbVM34dl5Za/JWHjfLERAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFEZyQ4Kyxla382jvCOO9250jIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFKYiAotrMyJoe08XdVk1ebkCAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACiO5gYgob2x1O4/2jjDeu0yOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEhhIipQqjInhrbzdFWTVR0BAZBEgABIIUAApBAgAFIIEAApbjlAx48fj82bN0etVouOjo44fPjwnNtffPHF6OjomLMNDQ2VtV4A2sQtB2hmZib6+vpi3759N7zP0NBQfP7557PboUOHbmuRALSfW/4e0PDwcAwPD3/rfSqVSnR3d//PiwKg/c3LZ0BHjx6Ne+65Jx5++OF45ZVX4uLFize8b6PRiHq9PmcDoP2VHqChoaE4ePBgjI2NxVtvvRXHjh2L4eHhuHLlynXvv2fPnqhWq7Nbb29v2UsCoAl1FLfxexAdHR3x/vvvx5YtW254n3/84x9x//33xwcffBCbNm265vZGoxGNRmP2cr1ej97e3piamorOzs7/dWlAG2jGn+IpSzv/FE+9Xo9qtfqdr+Pzfhr2+vXrY9WqVXHmzJnr3l6pVKKzs3POBkD7m/cAnT9/Pi5evBg9PT3z/VAAtJBbPgvu8uXLc45mxsfH49SpU7Fy5cpYuXJlvPHGG/Hcc89Fd3d3nD17NrZv3x4PPPBADA4OlrpwAFrbLQfo448/jqeffnr28ujoaEREjIyMxP79++Ovf/1r/P73v49Lly5FrVaLZ555Jn71q19FpVIpb9UAtLzbOglhPtzsh1dA+3MSQmtqmpMQAOB6BAiAFEZyA02rrLepmvGtvDLX1Kpv5zkCAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIWJqEDbK3NiaDtPV13oyaqOgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQQIABSCBAAKQQIgBQCBEAKAQIghQABkEKAAEghQACkMJIb4BaUNba6nUd73yxHQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQAoBAiCFAAGQQoAASCFAAKQwERUgQVmTVSOac7rqzXAEBEAKAQIghQABkEKAAEghQACkECAAUggQACkECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFIIEAApBAiAFAIEQIqmm4j6zZTAer2evBIAbsd3TX1tugBNT09HRERvb2/ySgC4HdPT01GtVm94e0dR5mDyEly9ejUuXLgQK1as+NY55/V6PXp7e+PcuXPR2dm5gCtcnOzvhWV/Lzz7vDxFUcT09HTUarVYsuTGn/Q03RHQkiVL4t57773p+3d2dnqyLCD7e2HZ3wvPPi/Htx35fMNJCACkECAAUrRsgCqVSuzatSsqlUr2UhYF+3th2d8Lzz5feE13EgIAi0PLHgEB0NoECIAUAgRACgECIIUAAZBCgABIIUAApBAgAFL8H4nN3o57QjHaAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import jax.numpy as jnp\n",
    "import jax\n",
    "import matplotlib.pyplot as plt\n",
    "from jax.config import config\n",
    "config.update(\"jax_enable_x64\", True)\n",
    "#Generate a random tridiagonal matrix\n",
    "\n",
    "n = 20\n",
    "d = jax.random.normal(jax.random.PRNGKey(0), (n, ))\n",
    "sub_diag = jax.random.normal(jax.random.PRNGKey(1), (n-1,))\n",
    "\n",
    "mat = jnp.diag(d) + jnp.diag(sub_diag, -1) + jnp.diag(sub_diag, 1)\n",
    "mat1 = jnp.abs(mat)\n",
    "mat1 = mat1/jnp.max(mat1.flatten())\n",
    "plt.spy(mat)\n",
    "q, r = jnp.linalg.qr(mat)\n",
    "plt.figure()\n",
    "b = r.dot(q)\n",
    "# b[abs(b) <= 1e-12] = 0\n",
    "b = b.at[abs(b) <= 1e-12].set(0.0)\n",
    "#b = jax.ops.index_update(b, jax.ops.index[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": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA/oAAAMwCAYAAACQlRtzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACoX0lEQVR4nOzdd3xUVf7/8Xd6AUIIhCogRZqioGsBCyKCKyiLgKioiIp9C+JasOzqrrvoWnZB1sJasC0KrAoKIigqCggiIEgTEKSEXhLSJ5n7+yO/O9+gIC0z5957Xs/HIw8TSOCNhwzzns+558Y5juMIAAAAAAAEQrzpAAAAAAAAoOpQ9AEAAAAACBCKPgAAAAAAAULRBwAAAAAgQCj6AAAAAAAECEUfAAAAAIAAoegDAAAAABAgFH0AAAAAAAKEog8AAAAAQIBQ9AEAAAAACJBE0wGOxcKFCzV+/Hh9/PHH2rx5s3bv3q3atWurfv366tChg7p27aru3burfv36pqMCAAAAABATcY7jOKZDHKnt27dr2LBhevPNNw/5uXfccYdGjx4dg1QAAAAAAJjnu4n+hg0bdP7552vdunWRH2vdurXat2+v2rVrq7CwUGvXrtXixYtVWFhoMCkAAAAAALHnq6Kfm5urrl27Rkp+165d9a9//Usnn3zyzz63tLRUM2fO1L59+2IdEwAAAAAAY3y1df+mm27Siy++KEm64oor9OabbyohIcFwKgAAAAAAvMM3RX/x4sXq2LGjJKlx48ZatmyZatSoYTgVAAAAAADe4pvb6z3//POR9++44w5KPgAAAAAAB+CLiX55ebmysrKUl5cnSVq9erVatmxpOBUAAAAAAN7ji4n+d999Fyn5NWvWVIsWLVRWVqZXXnlF3bp1U/369ZWSkqJGjRrp4osv1nPPPaeSkhLDqQEAAAAAiD1fTPRffPFF3XTTTZKkk046SVOnTlX//v01f/78g35NkyZNNHHiRJ1++umxigkAAAAAgHG+uL3exo0b9/v44osv1rJlyyRJbdq00emnn66EhAQtWbJECxculCRt2LBB559/vmbNmqXTTjvtoL92SUnJftN/x3FUWlqqOnXqKC4uLgp/GgAAAAAAjpzjONq3b59q1Kjxi33VF0V/7969kfe/++47SVJ6errGjh2ryy+/fL/P/fTTTzVgwADt3LlThYWFuuKKK7R8+XIlJycf8NceMWKEHnnkkZ/9+MaNG5WRkVF1fwgAAAAAAI5BXl6eGjdurL1796pmzZoH/TxfbN0fMmSIXnrppf1+7J133tFll112wM+fO3euzjnnHIXDYUnSyy+/rOuvv/6An/vTif7mzZvVrl27KkoOAAAAAEDV2rhxo4477riD/rwvJvqpqan7fdypU6eDlnz35/v27auJEydKkt5+++2DFv2UlBSlpKREPnZf9/DDRD8UCmn69Onq0aOHkpKSTMdBlLDOwdGsWTPt3r1bX331ldq2bRv58aKiItWvX1+StG7dOmVlZZmKiCowZswY3X333frNb36j1157bb+fGzFihB577DFde+21Gj16tKGEqAp79uzR8ccfL0nauXNn5PH5+uuv1zvvvKNBgwbp6aef5nE7AIYOHapXXnlFw4cP13333SdJ+uGHH9SxY0elpqZq06ZNrHOA8TzMDn5aZ3eif6jbzfui6FevXn2/j3+p5Ff+HLfoz5kz57B/L/c6h4yMDF8U/fT0dGVkZHj+LySOHuscHGVlZZKkrKys/R5fKq+rHx578MsSEyv+aXW/bytzX7hOTk5mnX2uuLg48n5WVlbk+YP7/eyuMY/b/hcfX3GTqurVq0e+b90n2HFxcaxzwPE8zA5+XOdDnSfni9vr1a5de7+PD2drfeVp2b59+7Rv374qzwUAR6K0tFSSDnpmCIIhFApJYp2Dzl3nxMREDu8NOHet/fLkHwAknxT9Nm3a7PfxTyf8B/LTrQwUfQCmUfTt4K4zpSDYKH/2YK0B+JEviv5JJ52038f5+fmH/JqfFvtfOpEQAKKtvLw8ckAoRT/YKAV2YJ3twVoD8CNfFP1mzZqpWbNmkY+XL19+yK9ZsWJF5P2srCxVq1YtKtkA4HC4U16Joh90bN23A+XPHqw1AD/yRdGXpL59+0bef++99w75+ZU/57zzzotCIgA4fBR9e7B13w6UP3uw1gD8yDdF/7bbbos8wM6ZM0eTJ08+6OfOnz9f77zzTuTjwYMHRzseAPyiykWfJ4vBRimwg3sXDdY5+PieBuBHvin6LVq00O233x75eODAgfuVedfnn3+uSy65ROXl5ZKks846S717945ZTgA4kMondLu3akIwUQrsUPl7GsHGWgPwI189Yj3++ONauHChvvjiCxUUFKhfv35q27atTj/9dCUkJGjJkiX65ptvIp/foEEDjR8/ntveADCOE/ftwTX6duAFHXuw1gD8yFdFPyUlRe+//75uu+02jRs3TlLFoXuVD95znXnmmZowYYIaN24c65gA8DMUfXtwjb4dKH/24DINAH7ku/2jNWvW1H//+199/vnnuvHGG9W6dWtVr15daWlpOv7443XllVfqnXfe0dy5cyn5ADyD8mcPCqAdWGd7sNYA/MhXE/3KzjvvPE7TB+AbTPTtwdZ9O1D+7MFaA/Aj3030AcCPKPr2YPeGHSh/9mCtAfgRRR8AYoCibw9KgR1YZ3uw1gD8iKIPADFA0bcHpcAOrLM9WGsAfkTRB4AYoOjbg7W2g3sSO/dWDz6KPgA/ougDQAxwQJs9KAV2YJ3twVoD8COKPgDEQElJiSSeKNqAUmAH1tkerDUAP6LoA0AMuEU/NTXVcBJEG7s37ED5swdrDcCPKPoAEAPFxcWSKPo24PZ6dqD82YO1BuBHFH0AiAGKvj0oBXZgne3gOA5rDcCXKPoAEANs3bcHpcAOrLMdysvLI+9zhwUAfkLRB4AYYKJvD26vZweKvh3cdZZYawD+QtEHgBig6NuDAmgHd52Z8gZbWVlZ5H2+pwH4CUUfAGLALfopKSmGkyDaKPp2cAsg6xxsTPQB+BVFHwBigIm+Pbi9nh14QccOlYt+QkKCwSQAcGQo+gAQAxR9e3B7PTtQ9O1QeZ3j4uIMpwGAw0fRB4AYoOjbgVtx2YN1tgPrDMCvKPoAEAPcXs8OlW/FRTEINgqgHVhnAH5F0QeAGGCibwd3275EMQg6CqAdWGcAfkXRB4AYoOjboXLR5w4LwUYBtAPrDMCvKPoAEAPcXs8O7iUaEsUg6CiAdmCdAfgVRR8AYoCJvh0qrzMndAdbWVmZJCkxMdFwEkQTRR+AX1H0ASAGKPp2cCf67NwIPgqgHVhnAH5F0QeAGKDo24FLNOxBAbQD6wzAryj6ABADFH07cBtFe1AA7eCuM5doAPAbij4AxAAF0A5s3bcHa20HXtAB4FcUfQCIASb6dmDrvj0o+nZwD12k6APwG4o+AMQABdAO7NywB0XfDkz0AfgVRR8Aoqy8vDzyZPFQBdBxnFhEQpRQ/uzhrnVycrLhJIimQxV9HrMBeBVFHwCizC0E0oGLPvdbDw52btiDF3XscLCiz+M2AK+j6ANAlLnlT2JLd9Cxdd8epaWlkij6QcfWfQB+RdEHgChzi35CQgK3aAo4prz2YK3tQNEH4FcUfQCIMqa89uDuCvag6NuBog/Aryj6ABBllD97UP7s4DgOa20Jij4Av6LoA0CUFRUVSaIQ2IDyZwe3/EmsddC5ZzFwdwUAfkPRB4Aoc4t+enq64SSINnZv2KHynTQo+sHGoYsA/IqiDwBRVlhYKImibwMm+naoXPSZ9Aabu9asMwC/oegDQJS5E/20tDTDSRBtFH07uOuckJCghIQEw2kQTWzdB+BXFH0AiDIm+vZg674d2M5tD168A+BXFH0AiDKKvj0oBXZgne3BRB+AX1H0ASDKKPr2oADagXW2B7s3APgVRR8Aooxr9O3B1n07UPTtwWF8APyKog8AUcZE3x4UQDuwzvZgog/Aryj6ABBlFH17UADtwJTXHqw1AL+i6ANAlLF13x5s3bcDL+jYg8P4APgVRR8AooyJvj0ogHZgO7c9+J4G4FcUfQCIMoq+PSgFdmCd7cFEH4BfUfQBIMoo+vZg674dKPr2YPcGAL+i6ANAlHGNvj0ogHZgne3BYXwA/IqiDwBRxkTfHu6LOkz0g42ibw+27gPwK4o+AEQZRd8ebtFnrYONom8P1hqAX1H0ASDK3KLP1v1gcxwnco0+ax1sbOe2BxN9AH5F0QeAKGPKawe35EsU/aBjymsPDuMD4FcUfQCIMrbu28F9QUei6Acd5c8e7N4A4FcUfQCIMoq+Hdx1TkxMVGJiouE0iCYm+vbgRR0AfkXRB4AochyH2+tZgnW2B0XfDuFwWGVlZZKY6APwH4o+AERRaWmpwuGwJCb6QUfRtwdF3w7uNF+i6APwH4o+AERRQUFB5H2KfrBR9O3BWtvBfUFH4kUdAP5D0QeAKHKLfnJyspKSkgynQTRxdwV7uHdYSE1NNZwE0VR5os/jNwC/oegDQBTl5+dLkqpXr244CaKNKa89WGs7uEU/KSlJ8fE8ZQbgLzxqAUAUUfTtQfmzh7vWTPSDjVvrAfAzij4ARJFb9KtVq2Y4CaKNom8Pd+s+ax1s3FoPgJ9R9AEgipjo26OwsFAS5c8GvKhjByb6APyMog8AUeQexkfRDz7Knz04jM8O7kSfog/Ajyj6ABBFTPTtQdG3B2ttB3eiz9Z9AH5E0QeAKKLo24PyZw8O47MDE30AfkbRB4Ao4jA+e7jlLz093XASRBuH8dmBw/gA+BlFHwCiiGv07cFE3x6stR04jA+An1H0ASCK2LpvD8qfHcrKylReXi6JrftB5xZ91hmAH1H0ASCKjrToO44TzTiIoiO5vR7r7F/uCzoSL+oEHXdXAOBnFH0AiKLDuUY/Li4uVnEQRYcz0Wet/a9y0T/QtduscXD80qGLrDMAr6PoA0AUsXXfHmzdt4M75U1JSVF8PE+jgoyJPgA/418oAIgiDuOzB0XfDqyzPSj6APyMog8AUcRE3x7uizrcSjHYKH/2YK0B+BlFHwCi6HCu0UcwUPTtwETfHhR9AH5G0QeAKKLo24Oib4dfOqANwULRB+BnFH0AiKJ9+/ZJkmrWrGk4CaKNom8Ht/wx0Q8+1hqAn1H0ASBKysvLIxP9jIwMw2kQbRR9OzDRtwcTfQB+RtEHgChxS75E0Q+68vLySCmg6AcbU157UPQB+BlFHwCiJC8vT5KUnJyslJQUw2kQTYWFhZH3KfrBxmF89qDoA/Azij4ARIlb9JnmB5+7bT8uLo4CGHBs3bcHaw3Azyj6ABAlFH17VL4+Py4uznAaRBNTXnuw1gD8jKIPAFFC0bcHt1G0h3uZBmsdfBR9AH5G0QeAKKHo24MT9+3hrnV6errhJIg2ij4AP6PoA0CUUPTtQdG3BxN9e1D0AfgZRR8AooSibw+Kvj2Y6NuDog/Azyj6ABAlFH17UPTtwUTfHm7R504aAPyIog8AUULRtwdF3x5M9O3BRB+An1H0ASBKKPr2oOjbg4m+PSj6APyMog8AUULRtwdF3x5M9O0QDodVUlIiiaIPwJ8o+gAQJRR9e1D07cFE3w5uyZco+gD8iaIPAFGSm5srSapRo4bhJIg2ir49mOjbwd22L1H0AfgTRR8AomTv3r2SpFq1apkNgqjLz8+XRNG3ARN9O7hFPz4+XomJiYbTAMCRo+gDQJTs2bNHkpSZmWk2CKJu3759krhMwwZM9O1Q+SC+uLg4w2kA4MhR9AEgSpjo28Mt+lymEWyO4zDRtwQn7gPwO4o+AERBcXFx5IkiE/3go+jbobS0VOXl5ZKY6Aed+4IO6wzAryj6ABAF7jQ/Li6O7dwWcO+wQNEPNrf8SUz0g46iD8DvKPoAEAVu0a9Zs6bi43moDTqu0beDe31+YmKikpKSDKdBNFH0AfhdYJ59Dhs2THFxcZG3448/3nQkABZzD+Lj+nw7sHXfDlyfb4+ioiJJFH0A/hWIoj9//nyNHDnSdAwAiHAn+lyfH3yO41D0LcGJ+/Zgog/A73xf9EOhkIYMGaJwOGw6CgBEUPTtUVJSolAoJImiH3RM9O1B0Qfgd74v+o8//riWLl0qSRo4cKDhNABQga379nCn+RJFP+iY6NuDog/A73xd9FeuXKlHH31UknT11Vere/fuhhMBQAUm+vZwi356eroSEhIMp0E0uUWfiX7wUfQB+J1vi77jOBoyZIhKSkpUq1YtPf3006YjAUAEE317cH2+PfLz8yVJ1atXN5wE0UbRB+B3vi36zz33nGbPni1JeuKJJ1S3bl3DiQDg/zDRt0deXp4kir4N3KLPWgcfRR+A3/my6G/cuFH33XefJOncc8/VDTfcYDgRAOxv165dkqSsrCzDSRBt7kQ/IyPDcBJEm7vWTPSDj6IPwO98WfRvv/127du3T8nJyXrhhRcUFxdnOhIA7Mct+rVr1zacBNHG1n17MNG3B0UfgN/5rui/9dZb+uCDDyRJ9957r9q2bWs4EQD8HEXfHhR9ezDRtwdFH4Df+aro79q1S7///e8lSa1atdIDDzxgOBEAHBhF3x5co28PXtSxB0UfgN8lmg5wJO68807t2LFDkvT8888rJSXlmH/NkpISlZSURD52n7CFQiGFQqFj/vWjyc3n9Zw4Nqyz/ziOEyn6GRkZh1y7yj/vh8ce7M+9w0L16tUPuXbhcDjyX9bZf9znCOnp6QddP3eNJR63/cy9lWJycvIB1/Gnj9sILp6H2cFP63y4GX1T9KdPn67XX39dknTdddepa9euVfLrjhgxQo888sgBfz+/vIo7Y8YM0xEQA6yzfxQVFUUehBcsWKDvvvvuFz+/8gP2zJkz2RbsM4sXL5ZUsYtj6tSpv/i5a9askSRt3rz5kJ8L7/nhhx8i/z3Y+m3ZsiXyPo/b/pWTkyNJWr58+QHXetu2bZIqXthlne3AOtvBD+vs7jg6FF8U/YKCAt1yyy2SKrbBPvnkk1X2aw8fPlzDhg2LfJyXl6fGjRurR48enj9BORQKacaMGerevbuSkpJMx0GUsM7+s379eklSamqqLrvsskMeGFpaWhp5v2vXrsrOzo5mPFSxd955R5J06qmnqmfPnr/4uYsWLZIkNWrU6JCfC+956qmnJEmdO3c+6Pr997//jbzP47Z/Pfjgg5Iq7u7UrVu3n/28+zgvsc5Bx/MwO/hpnd3dZYfii6L/wAMPRB5Qn3rqKdWpU6fKfu2UlJQDXgKQlJTk+UV2+Skrjh7r7B+5ubmSKl6YTE5OPuTnO44TeZ919h/3H9zatWsfcu3i4+Mj/2Wd/cc9dT8zM/Og6+euscT3s58VFRVJqrj86kBrWPnHWGc7sM528MM6H24+zx/Gt3DhQj3zzDOSKiZd1113neFEAPDLOIjPLu4LOzVr1jScBNHG7fXswWF8APzO8xP9JUuWRA622bBhg84666yDfq57UJ9UcY1c5c996KGH1KtXr+gFBYD/j6Jvl71790qqmPIi2Dh13x4UfQB+5/miX9natWu1du3aw/rc0tJSzZs3L/Jx5RcBACCaKPp2YaJvD3eiz4GZwUfRB+B3nt+6DwB+s3PnTkkUfVsw0bdDOBxm674lKt/mlKIPwK88X/QHDx4sx3EO6+2VV16JfF3Tpk33+7nBgweb+0MAsIq7g6hevXqGkyDaHMdhom8J977qEhP9oGOtAQSB54s+APjN9u3bJUl169Y1nATRVlBQoPLycklM9IPOnebHx8crLS3NcBpEk7vWSUlJh3XnFADwIoo+AFQxir493Gl+YmIiW3wDrvJBfHFxcYbTIJo4iwFAEFD0AaCKUfTt4V6fX7NmTcpfwOXl5Uni+nwbUPQBBAFFHwCqmFv0s7OzDSdBtLkTfbbtBx9nMdiDog8gCCj6AFCFSktLtWfPHklM9G1QeaKPYKPo24OiDyAIAlX0K5/Qv379etNxAFjIvbVefHy8srKyDKdBtLkv6jDRDz536z5FP/go+gCCIFBFHwBMq7xtPz6eh9ig2717tySpdu3ahpMg2pjo24OiDyAIeBYKAFWIg/js4hZ9dm8EH0XfHhR9AEFA0QeAKrRt2zZJFH1bUPTt4Rb9jIwMw0kQbRR9AEFA0QeAKrR582ZJUqNGjQwnQSxQ9O3BRN8eFH0AQUDRB4AqlJOTI4mibwuKvj0o+vag6AMIAoo+AFQht+g3bNjQcBLEAofx2YOibw+KPoAgoOgDQBVyt+5T9O2wa9cuSUz0bUDRtwdFH0AQUPQBoAox0bcLW/ftQdG3B0UfQBBQ9AGgioTDYW3ZskUS1+jbIBwOa8+ePZIo+jag6NuDog8gCCj6AFBFdu7cqVAoJEmqX7++4TSItry8PIXDYUlSrVq1DKdBtOXl5Umi6NuAog8gCCj6AFBF3G37devWVVJSkuE0iDZ32356erpSU1MNp0E0hUIhFRYWSpIyMjIMp0G0UfQBBAFFHwCqCLfWs8vOnTslsW3fBnv37o28n5mZaSwHYmPfvn2SeFEHgL9R9AGginDivl3cop+dnW04CaLN3b1Rs2ZNJSQkGE6DaHPPY6DoA/Azij4AVBFO3LfLjh07JFH0beAeushZDMFXUlKikpISSRR9AP5G0QeAKsLWfbtQ9O1B0beHu21fkmrUqGEwCQAcG4o+AFQRtu7bhaJvD3frPucxBJ+7bb969epcpgHA1yj6AFBF2LpvF7fo16lTx3ASRBsTfXu4t1Fk2z4Av6PoA0AVYeu+XTiMzx4UfXtQ9AEEBUUfAKpAKBTS9u3bJUkNGjQwnAaxwNZ9e7B13x7u1v2aNWsaTgIAx4aiDwBVYMOGDXIcR6mpqapbt67pOIgBir49mOjbg4k+gKCg6ANAFVi3bp0k6fjjj1dcXJzhNIgFrtG3hzvRp+gHH0UfQFBQ9AGgCrhFv1mzZoaTIBZKS0sjhYCJfvC5E3227gcfW/cBBAVFHwCqAEXfLtu2bZMkJSYmMuW1AFv37cFEH0BQUPQBoApQ9O3iFv169eopPp5/SoOOrfv2oOgDCAqenQBAFaDo22Xr1q2SpPr16xtOgmhzHIdbKVqEog8gKCj6AFAFKPp2cYt+vXr1DCdBtO3bt0+hUEiSVLt2bcNpEG1cow8gKCj6AHCMCgoKtH37dkkUfVu4W/eZ6AefO81PT09Xenq64TSINrfoM9EH4HcUfQA4RuvXr5dUMQE61mt4HcepgkSItmPdus86+4db9I/0NoqssT8d6cGLrDMAr6LoA8AxOtZt+3FxcVUZBzFwtBN91tp/jrTos8b+drhFn3UG4HUUfQA4Rlyfbx+u0bfH0U704U979+6VJGVmZhrNAQDHiqIPAMeIom8fTt23B0XfHqFQSPn5+ZK4lSIA/6PoA8Axcot+8+bNDSdBLDiOo5ycHElSgwYNDKdBtFH07eFO8yUm+gD8j6IPAMdo9erVkij6tsjNzVVBQYEk6bjjjjOcBtFG0beHe31+RkaGEhISDKcBgGND0QeAY1BWVqbvv/9ektS2bVvDaRALmzZtklRxT/W0tDTDaRBtFH17cH0+gCCh6APAMVi3bp1CoZDS0tLUpEkT03EQAxs3bpTENN8W27dvl0TRt8GR3loPALyMog8Ax2DFihWSpNatWys+nodUG7gTfYq+HTh40R4UfQBBwrNSADgGK1eulMS2fZu4Rb9x48aGkyAWKPr2cIs+W/cBBAFFHwCOgTvRb9OmjeEkiBUm+vbIz8+PHLzIHRaCz71Gn4k+gCCg6APAMXCLPhN9e3CNvj3caX61atVUvXp1w2kQbWzdBxAkFH0AOEqO40S27jPRtwcTfXuwbd8uFH0AQULRB4CjtHXrVuXm5io+Pl4nnHCC6TiIEYq+PbZs2SKJom8L91aKtWvXNpwEAI4dRR8AjpI7zW/WrJlSU1MNp0Es5OXlad++fZKkRo0aGU6DaHMn+lyfb4ddu3ZJ4laKAIKBog8AR4nr8+3jTvMzMzO5ZtsCbN23i1v0megDCAKKPgAcpW+//VaSdOKJJxpOglhxD+Lj1np2oOjbha37AIKEog8AR2nx4sWSpI4dO5oNgpjh+ny7UPTt4TgOW/cBBApFHwCOQllZmZYsWSKJom8Tir5duEbfHrm5uSovL5fERB9AMFD0AeAofP/99youLla1atXUsmVL03EQIxs2bJBE0bcFp+7bw53mp6enc7gqgECg6APAUVi0aJEk6ZRTTlF8PA+ltvjhhx8kSc2bNzecBNFWXl6u7du3S6Lo24Bt+wCChmenAHAU3KLPtn27rF27VpLUokULw0kQbbt27VJ5ebni4uKUnZ1tOg6ijIP4AAQNRR8AjoJ7EF+HDh2M5kDslJSURK7RZ6IffO71+XXq1FFSUpLhNIg2bq0HIGgo+gBwhBzHYaJvoXXr1slxHFWrVk1169Y1HQdRxkF8dmGiDyBoKPoAcIQ2btyo3bt3KyEhQSeeeKLpOIiRytv24+LiDKdBtHEQn1127NghSbyIByAwKPoAcITmzZsnSWrfvj2nM1vEPYiP6/PtQNG3i3vwIkUfQFBQ9AHgCM2ZM0eS1LlzZ8NJEEscxGeXH3/8UZLUpEkTw0kQC27Rr1evnuEkAFA1KPoAcITmzp0riaJvG4q+XTZs2CBJatq0qeEkiAUm+gCChqIPAEeguLhYCxculCR16tTJcBrEEkXfLkz07ULRBxA0FH0AOALffPONQqGQ6tatq2bNmpmOgxgJh8Nco28Rx3EiRZ+Jvh22bdsmiaIPIDgo+gBwBCpfn8/J6/bYsmWLSkpKlJCQwITXArm5ucrPz5ckNW7c2HAaRFtBQYEKCwslUfQBBAdFHwCOgHt9Ptv27eJu22/atKkSExMNp0G0udP8OnXqKD093XAaRJu7bT81NVXVq1c3nAYAqgZFHwAOk+M4mj17tiQO4rPNqlWrJEknnHCC4SSIBQ7is0vl6/PZqQUgKCj6AHCYvvvuO23fvl1paWk6/fTTTcdBDC1fvlySdOKJJxpOgljgID67cH0+gCCi6APAYfrkk08kSeedd55SUlIMp0EsLVu2TJLUrl07w0kQC0z07ZKTkyNJatiwoeEkAFB1KPoAcJjcot+tWzfDSRBr7kSfom8HJvp22bx5sySpUaNGhpMAQNWh6APAYSgrK9Pnn38uiaJvm9zc3EgRaNu2reE0iAUm+nZhog8giCj6AHAYvv76a+3bt09ZWVnq0KGD6TiIoRUrVkiqmPZlZmaaDYOYYKJvFyb6AIKIog8Ah8Hdtt+1a1fFx/PQaRO27dulpKREW7ZskUTRtwVFH0AQ8WwVAA7D1KlTJUndu3c3nASxxkF8dtm0aZOkinuqZ2dnG06DWGDrPoAgougDwCFs375dX331lSTpkksuMZwGscZE3y5r1qyRJDVv3px7qlugqKhIu3fvlsREH0CwUPQB4BCmTJkix3F06qmnRv2JoOM4Uf31ceQo+nZZvXq1JKlVq1aGkyAW3Gl+WloaZ3AACBSKPgAcwvvvvy9JuvTSS6Py6zM19K7c3NzICewUfTt8//33kqQTTjjBcBLEQuVt+0fyWMzjNgCvo+gDwC8oLi7W9OnTJUWv6MO7Fi9eLKniULasrCyzYRAT7kSfom8HDuIDEFQUfQD4BTNnzlRBQYEaNmyoU0891XQcxNjChQslSaeddprhJIgVd6LP1n07uEWfg/gABA1FHwB+wdtvvy1Juuyyy9iqaSG36PMijx1KS0u1fv16SUz0beFu3WeiDyBoKPoAcBDFxcV69913JUlXXXWV4TQwgaJvl3Xr1ikcDqtatWpq0KCB6TiIAbbuAwgqij4AHMSHH36offv2qXHjxurUqZPpOIixgoICrVy5UhJF3xaVD+JjB48dKh/GBwBBQtEHgIMYN26cJOmKK65QfDwPl7ZZsmSJwuGwGjRooPr165uOgxjgID77MNEHEFQ8cwWAA8jNzdUHH3wgSbryyisNp4EJbNu3j1v0OYjPDo7jcI0+gMCi6APAAbz11lsqKipS27ZtKXqWoujbp/LWfQTf9u3bVVxcrLi4OLbuAwgcij4AHMCLL74oSRoyZAjX6lpq/vz5kri1nk1WrVoliaJvi3Xr1kmqmOanpKQYTgMAVYuiDwA/8e2332rBggVKSkrStddeazoODNi7d6+WLVsmSRzEaIk9e/ZErtc+8cQTDadBLLhFv1mzZoaTAEDVo+gDwE+89NJLkqTevXsrOzvbcBqY8NVXX8lxHLVs2VJ169Y1HQcx8N1330mSmjRpopo1axpOg1ig6AMIMoo+AFSSl5ensWPHSpJuuukms2FgzJw5cyRJZ599tuEkiJWlS5dKkk466STDSRArbtFv3ry54SQAUPUo+gBQydixY7Vv3z61adNG3bt3Nx0HhsyePVuS1LlzZ8NJECtu0W/fvr3hJIgVJvoAgoyiDwD/X3l5uUaNGiVJ+sMf/qD4eB4ibVRWVqZ58+ZJoujbxN26T9G3B0UfQJDxLBYA/r8pU6Zo7dq1qlWrFofwWWzJkiUqKChQzZo11a5dO9NxEAOO4zDRt0xZWZk2bNggiaIPIJgo+gCgiif6f/vb3yRJt9xyi6pVq2Y4EUz58ssvJUlnnXUWuzossWnTJuXm5ioxMVFt2rQxHQcxsGnTJpWVlSk5OVkNGzY0HQcAqhzPYABA0scff6z58+crLS1Nd955p+k4MGjmzJmSpPPPP99sEMSMO81v1aqVkpOTDadBLLjb9o8//nhe0AMQSDyyAYCkRx99VJJ08803czs1i5WVlemzzz6TJHXr1s1sGMQM2/btw/X5AIKOog/Aep988olmzZql5ORk3X333abjwKBvvvlGubm5yszM1Kmnnmo6DmJk8eLFkqSTTz7ZbBDEDEUfQNBR9AFYLRwO695775VUcW1+o0aNDCeCSZ988omkim37CQkJhtMgVr7++mtJ0q9+9SvDSRArFH0AQUfRB2C1CRMm6JtvvlGNGjX00EMPmY4Dw9yif+GFFxpOgljZvXu31q5dK4mib5PVq1dLklq0aGE4CQBEB0UfgLWKioo0fPhwSdLdd9+t7Oxsw4lgUlFRkWbPni2J6/NtsmDBAkkVhS8rK8twGsSC4zhauXKlJKlt27aG0wBAdFD0AVjrscce07p169SoUSNO2oc+++wzlZSUqFGjRmrdurXpOIgRd9v+6aefbjgJYmXLli3Ky8tTQkICE30AgUXRB2ClNWvW6PHHH5ck/etf/1L16tUNJ4JpH3zwgSTpkksuUVxcnOE0iBWKvn3caX7z5s2VkpJiOA0ARAdFH4B1HMfR7373O5WUlKh79+7q16+f6UgwzHGc/Yo+7EHRt8+KFSsksW0fQLBR9AFYZ8KECZo2bZqSk5M1evRoprfQ0qVLtWHDBqWlpXF9vkVycnKUk5Oj+Ph4bqdoEXei36ZNG8NJACB6KPoArLJlyxbddtttkqR7771XrVq1MpwIXuBO87t166a0tDTDaRAr8+fPlyS1a9dO1apVM5wGscJEH4ANKPoArOE4jm688Ubt3r1bHTt21IMPPmg6Ejzi/ffflyRdeumlhpMglr744gtJ0tlnn204CWKJiT4AGySaDnAk1q9frxkzZujzzz+PbLPMz89XjRo1dNxxx6lTp04aOHCgunTpYjoqAA964YUX9OGHHyolJUVvvPGGkpOTTUeCB2zcuFFfffWV4uLiuD7fMm7RP/fccw0nQazk5eVp8+bNkij6AILNF0V/0aJFuvXWWyNb7H5qz5492rNnj5YuXaoxY8bo/PPP16uvvqomTZrEOCkAr1qyZImGDRsmSRoxYoTatWtnOBG8YuLEiZIqyl7Dhg0Np0Gs5Ofna+HChZIo+jZZtWqVJKl+/frKzMw0GwYAosgXRX/VqlU/K/mtWrXSSSedpDp16mjv3r2aM2eONm3aJKniXsidOnXSF198oebNm5uIDMBD9u7dq759+6qoqEgXXXSR/vCHP5iOBA8ZP368JGnAgAGGkyCW5s6dq/LycjVt2pTBgEW4Ph+ALXxR9F0tW7bUkCFDdM0116hRo0b7/Vw4HNbYsWP1u9/9ToWFhcrJydHVV1+tOXPmcKI2YLFwOKxrr71Wa9eu1fHHH68333xT8fEcT4IKP/74Y2TbPrdZtAvb9u20dOlSSWJXF4DA80XRb9CggV555RVde+21SkhIOODnxMfH64YbblCtWrXUt29fSdJXX32l6dOn66KLLoplXAAe8uc//1kffPCBUlNT9b///U+1a9c2HQkeMmHCBElSly5dVL9+fcNpEEsUfTstWrRIktSxY0fDSQAgunwx1urSpYsGDx580JJf2WWXXaYzzjgj8vGUKVOiGQ2Ah7344ot69NFHJUnPP/8898nGfhzH0auvvipJuvLKKw2nQSwVFxdr3rx5kij6NnEcR4sXL5YkdejQwWgWAIg2XxT9I1X5Njnr1683FwSAMR9++KFuvfVWSdJDDz2k6667znAieM2iRYv03XffKSUlRVdccYXpOIih2bNnq6ioSA0aNODkdYts3rxZu3btUmJiok488UTTcQAgqgJZ9Ctfk19eXm4wCQATZs+ercsvv1zl5eW67rrr9Mgjj5iOBA8aO3asJKlPnz6cvm2Z6dOnS5J69OjBOT4Wcaf5bdu2VWpqqtkwABBlgSz67kErktS4cWODSQDE2ldffaWLL75YBQUF6tGjh8aMGcMTefxMSUmJ3nzzTUnS4MGDzYZBzFUu+rCHe30+2/YB2CBwRX/Dhg2aOXNm5OMLL7zQYBoAsfT111/roosu0r59+9S1a1e9++67Sk5ONh0LHjR58mTt3r1bDRo0UPfu3U3HQQxt27YtMtll7e3irjsH8QGwgS9O3T8Sw4YNi2zXb9KkiS699NJf/PySkhKVlJREPs7Ly5MkhUIhhUKh6AWtAm4+r+fEsWGdD89nn32mvn37Kj8/X+eee67eeecdJSUl+eL/W1lZWeR9Pzz2BMGzzz4rqWKaHw6HFQ6HY/L7ur9POBxmnQ358MMPJVWUvczMzCpfh8p/l1hjb3GL/kknnXTMa1P561nnYON5mB38tM6HmzFQRf/VV1/V//73v8jHI0aMUEpKyi9+zYgRIw54/e706dOVnp5e5RmjYcaMGaYjIAZY54ObO3eunnrqKZWVlal9+/a6/fbb9fnnn5uOddgqnyXy6aefKiMjw2Ca4Nu4caM+++wzxcfHq3nz5po6dWrMfu/Vq1dLkjZt2hTT3xf/xz2bIVprn5OTI6nihHcet72joKBAP/zwg6SKXR3HuvY7duyQxDrbhHW2gx/WubCw8LA+LzBFf8GCBZETtiXpqquu0sCBAw/5dcOHD9ewYcMiH+fl5alx48bq0aOH559sh0IhzZgxQ927d1dSUpLpOIgS1vngHMfRqFGj9MQTTygcDqtPnz567bXXfHfIUuWif/7553M/9yi78847JUmXXHJJzO/G4E4UGzVqpJ49e8b090bF7pkbbrhBknTbbbfpvPPOq/Lf46233oq8z+O2d3zyySeSpOOPP75K7rKxcePGyPusc7DxPMwOflpndwf6oQSi6K9bt06XXnqpiouLJUknn3yynn/++cP62pSUlANO/ZOSkjy/yC4/ZcXRY533V1JSottuu02vvPKKJOnmm2/Ws88+q4SEBMPJjlx8/P8dl8I6R1dubq5ef/11SdJvf/vbmP+/dv9+xsfHs84GzJ49W7t371bt2rXVpUsXJSZW/dMgvp+9acGCBZKkTp06VcmaVP41WGc7sM528MM6H24+3xf9LVu2qHv37tq6daukiq1406ZN8/w0HsDR27Ztm/r166fZs2crPj5eTz/9tH7/+99zuj4O6YUXXlBeXp7atm2rbt26mY6DGHvvvfckSZdeemlUSj68a+7cuZIqij4A2MDX/8rt2rVL3bt319q1ayVJDRo00Mcff6wGDRoYTgYgWj755BNdc8012rp1q2rWrKm3335bF110kelY8IGSkhL961//kiTdfffd+01eEXyO42jSpEmSpD59+pgNg5hyHEdfffWVJIo+AHv49llOXl6eLrroIi1btkySVKdOHX388cdq1qyZ4WQAoiEUCumBBx6I7OBp166dvvrqK0o+Dtsbb7yhLVu2qFGjRrr66qtNx0GMLVmyROvXr1daWhq31bPM999/r927dys1NVUnn3yy6TgAEBO+nOgXFBSoZ8+e+uabbyRJNWvW1LRp09SuXTvDyQBEw8qVK3X99ddHJjI33XST/vWvf/nmzhgwLxwO64knnpBUcRhfcnKy4USINfeuPN27d+exwzLutv1f/epXfO8DsIbvJvrFxcXq3bu3Zs+eLUlKT0/XlClTdNpppxlOBqCqlZWV6fHHH1eHDh301VdfKSMjQ2+//bbGjBnDE3UckUmTJmnVqlWqWbOmbr75ZtNxEGOO42jcuHGSVCUnrsNf2LYPwEa+muiHQiH169dPM2fOlFRxYv6kSZN09tlnG04GoKotWrRIt9xyi77++mtJ0kUXXaQxY8aoSZMmhpPBb8LhsB5++GFJ0u23364aNWqYDYSY++abb7RmzRqlpaWpd+/epuMgxubMmSOJog/ALr6Z6JeXl2vgwIGaOnWqJCkxMVHjx4/XhRdeaDgZgKq0Y8cO3XLLLTrttNP09ddfKzMzU6+88oo+/PBDSj6Oyvjx47VkyRJlZGTorrvuMh0HBrjT/N69e6t69eqG0yCWdu7cqaVLl0oSgyEAVvFF0XccRzfeeKMmTpwoqeIeta+//jqvygMBEgqFNGrUKLVq1UpjxoyR4zgaOHCgli1bpsGDB3PrPByVUCikhx56SJL0xz/+UbVr1zacCLEWDof19ttvS5Kuuuoqw2kQa5999pkkqX379qpbt67ZMAAQQ77Yuv/cc8/p1VdfjXzcokULffnll/ryyy8P6+tHjx4drWgAjlE4HNb48eP15z//Wd9//70kqUOHDnrmmWd0zjnnGE4Hv3v11Ve1Zs0aZWdna+jQoabjwICZM2dq8+bNqlmzpn7961+bjoMYcy/37Nq1q+EkABBbvij627dv3+/j1atXa/Xq1Yf99RR9wHscx9H777+vhx56SEuWLJFUcZvMRx99VEOGDFFCQoLhhPC7oqIiPfLII5Kk+++/n2vzLfXSSy9JkgYOHKiUlBTDaRBrn376qSTpggsuMJwEAGLLF0UfQHCUl5frvffe02OPPaYFCxZIkjIyMvTHP/5RQ4cOpYyhyjz55JPatGmTGjdurFtvvdV0HBiwe/duvfvuu5KkG2+80XAaxFpOTo5Wrlyp+Ph4denSxXQcAIgpX1yj//DDD8txnKN+A2BecXGxxowZo7Zt26p///5asGCB0tPTNXz4cK1bt04PPfQQJR9VZsOGDRoxYoQk6R//+IdSU1MNJ4IJb775pkpKSnTKKafo1FNPNR0HMeZO8zt27KjMzEyzYQAgxpjoA4iqTZs2acyYMRozZoy2bdsmSapVq5buuOMO/e53v+NwJETFPffco6KiIp177rncN91SjuPoP//5j6SKaT4HetpnxowZkti2D8BOFH0AVc5xHH3yySd69tlnNXnyZJWXl0uSGjdurGHDhmnIkCHc4gpRM2vWLL399tuKj4/XqFGjKHiW+uKLL7R06VKlpaXpmmuuMR0HMRYOh/Xhhx9Kki6++GLDaQAg9ij6AKrMpk2b9MYbb+iVV16JnKAvSeeff75uv/129enTR0lJSQYTIuhKS0t1++23S5JuuukmdejQwWwgGPPMM89Ikq655hrVqlXLcBrE2oIFC7R9+3bVqFGDO7gAsBJFH8AxKSgo0DvvvKPXXntNn3zySeRcjBo1amjQoEG67bbbdOKJJxpOCVs89thjWrZsmbKzs/W3v/3NdBwYsnHjxsghfL/73e8Mp4EJU6ZMkST16NGDF5gBWImiD+CIFRYWatq0aZo4caLef/995efnR37uvPPO06BBgzRgwAAO10NMLV++XI8++qgkadSoUapdu7bhRDDl2WefVXl5ubp06aL27dubjgMD3KLfq1cvw0kAwAyKPoDDUlBQoKlTp2rixImaMmWKCgoKIj/XokULDRo0SNdee62aNWtmMCVsFQ6HddNNNykUCqlXr14cwGex3NxcPfvss5KkoUOHmg0DI7Zu3apvvvlGEtfnA7AXRR/AQf3444+aOnWqpk6dqk8++URFRUWRn2vatKn69++v/v3768wzz+TAMxj1z3/+U3PmzFH16tX17LPP8vfRYs8++6zy8vLUrl079e7d23QcGDBp0iRJ0umnn6769esbTgMAZlD0AUSEQiHNnj07Uu6XLVu23883b95cl19+ufr376/TTjuNMgVP+Pbbb3X//fdLkp588kk1adLEcCKYUlhYqH/+85+SpPvuu0/x8fGGE8GECRMmSJL69etnOAkAmEPRBywWDof17bff6tNPP9Wnn36qWbNmKS8vL/LzCQkJ6ty5s3r27KmePXuqffv2lHt4SlFRka6++mqVlpaqd+/euvnmm01HgkEvv/yyduzYoaZNm+rKK680HQcG7NixQ59++qkk6fLLLzecBgDMoegDFnEcR8uXL9fMmTP16aef6vPPP9fu3bv3+5zs7GxdfPHF6tmzp3r06MFtqeBpw4cP17Jly1SvXj29+OKLvBBlsVAopCeeeEKSdM8993DSuqXee+89hcNhnXrqqWrevLnpOABgDEUfCLDCwkJ9/fXXmjt3rubMmaO5c+dq586d+31O9erVdd5556lr167q2rWrOnbsyHZX+ML06dM1cuRISRWT3OzsbMOJYNKLL76oDRs2qF69err++utNx4Eh7rZ9pvkAbEfRBwLCcRytW7dO8+fPj5T6xYsXq6ysbL/PS0tL0znnnBMp9qeddhqTL/jOpk2bdM0110iS7rjjDvXs2dNwIphUUFCgv/zlL5Kkhx56SGlpaYYTwYRt27Zp5syZkqT+/fsbTgMAZlH0AR9yHEdr167VN998o4ULF0b+u2fPnp99bsOGDdW5c2d16tRJnTt3VseOHZWSkmIgNVA1SktLNWDAAO3YsUMdOnSIbNeGvUaOHKmtW7eqWbNmuummm0zHgSFvvvmmysvLddZZZ6lly5am4wCAURR9wOOKi4v1ww8/6PXXX9eKFSsipT43N/dnn5ucnKyTTz45Uuo7d+6sxo0bc90yAuWee+7R3LlzVbNmTU2cOJHpreV27dqlxx9/XJL017/+VcnJyYYTwZRXX31VknTdddcZTgIA5lH0AY8oLy/XDz/8oKVLl2rp0qX67rvvtHTpUq1evVrhcPhnn++W+tNOOy3ydtJJJ/EkF4E2fvz4yHX5r732mlq0aGE4EUwbMWKE8vLydPLJJ+uqq64yHQeGLF68WEuWLFFycrIGDBhgOg4AGEfRB2KsqKhIq1ev1qpVqyJvK1eu1LJly1RUVHTAr6levbo6duyok08+WR06dNBpp52mE088kVIPqyxdulQ33nijpIp7pPfu3dtwIpi2cuXKyAs/jz32GAeJWsyd5vfu3VtZWVmG0wCAeRR9IArC4bBycnL2K/Pu248//ijHcQ74dampqWrXrp3at2+vk046Se3bt1ebNm20aNEi9erVi0PzYK3t27fr0ksvVX5+vi644AL99a9/NR0JhjmOoz/84Q8qKytTr169dPHFF5uOBEOKiooiRX/w4MFmwwCAR1D0gaNUUFCgdevW6YcffvjZ27p161RcXHzQr83MzFTr1q33ezvppJPUsmVLJSQk7Pe5oVBIixcvjvKfBvCukpIS9e3bVz/++KNatmypCRMmKDGRf75sN2nSJE2fPl3Jycn617/+ZToODBo/frz27Nmjpk2b6te//rXpOADgCTxTAg6iuLhYmzZt0oYNG/Tjjz/+rNRv27btF78+ISFBzZs3V5s2bX5W6rOzszkgDzgMjuPo1ltv1ezZs1WzZk29//77bMuFioqKdOedd0qS7rrrLk5Yt9xzzz0nSbrlllt+9mI5ANiKog8rlZeXa9u2bdqwYYM2bNigjRs3/uy/27dvP+SvU6tWLTVv3vyAb40bN2arPXCMnnzySY0dO1YJCQkaP3682rRpYzoSPOAvf/mL1q9fr0aNGun+++83HQcGLVy4UPPmzVNSUlLkDA8AAEUfAeM4jvLy8rRlyxbl5OQoJyfnZ+9v3LhRmzZtUllZ2SF/vfT0dDVp0kSNGzf+WZFv1qyZatWqFYM/FWCnd999V/fee68k6Z///Kd69OhhOBG8YNGiRXriiSckSaNHj1b16tUNJ4JJ7mUb/fv3V926dc2GAQAPoejDF8LhsHbv3q3t27dr27ZtByzw7vuFhYWH9WsmJCSoUaNGaty4sZo0aRIp9JX/W6tWLbbYAwbMmjVLV111VWTr/m9/+1vTkeABZWVluvHGG1VeXq7LL79cffr0MR0JBm3atEnjxo2TJA0bNsxwGgDwFoo+jCksLNT27dv3e9u2bdsBf2znzp0qLy8/7F+7Zs2aatiwoRo0aKCGDRtG3ho0aKDjjjtOTZo0Uf369TnQC/CgpUuXqnfv3iopKdFvfvMbPfPMM7zgBknSU089pUWLFqlWrVp65plnTMeBYc8884zKysrUpUsX/epXvzIdBwA8hZaDKlFaWqpdu3Yd8q1yeS8oKDji3ycrK0t169Y9YIGv/H56enoU/pQAom39+vW66KKLlJubq3POOUfjxo3jBTlIkr799lv96U9/kiQ9/fTTqlevnuFEMCkvL08vvPCCJOmPf/yj4TQA4D08e8J+ysvLlZubq7179x5WcXff8vPzj+r3S0lJUb169VS3bt3I208/dt+ys7M53A4IsJ07d+qiiy7Sli1bdOKJJ2ry5MlKS0szHQseUFxcrKuvvlqlpaXq3bu3rrvuOtORYNjIkSOVm5urtm3bqmfPnqbjAIDnUPQDxnEc7du3T3v37tWePXu0d+/eyNvhfLxv376j/r3j4+OVlZWl2rVrH/Ttp2W+evXqbMkFoIKCAvXq1Uvff/+9GjdurGnTpnHYJSLuu+8+LVu2TPXq1dOLL77IvxuWy83N1dNPPy1J+tOf/qT4+HjDiQDAeyj6PrV06VL96U9/0tq1a/XII48oNzdXe/bsUW5ursLh8DH/+unp6Ycs7T99y8zM5B9b4Bg5jmM6QsyFQiFdfvnlmj9/vrKysvTRRx/puOOOMx0LHjF9+nSNHDlSkvTyyy8rOzvbcCKYNnLkSO3du1ft2rXT5ZdfbjSLjY/ZAPyBou9T+fn5eu+99w7688nJyapVq5YyMzMj/6389tMfq/xxzZo1lZycHLs/DGA5m6eT4XBYN954oz788EOlpaVpypQpatu2relY8IgdO3Zo8ODBkqTbb7+dLdrQ3r17I9P8P//5z0pISDCSw+bHbQD+QNH3qZYtW2rUqFFat26dunbtqjp16uxX2lNTU01HBIBDuu+++/T6668rISFBEyZM0FlnnWU6EjyivLxc11xzjbZs2aI2bdroiSeeMB0JHvCvf/1Lubm5OvHEE9W/f3/TcQDAsyj6PpWdna1bb71VU6dO1a9//WsOqQPgO08//XSkvL300kvq1auX4UTwkkcffVTTp09XWlqaJkyYwN1UoJycnMhjxsMPP8zlggDwC3iEBADE3Jtvvqm77rpLkvTYY49xijr289FHH+mRRx6RJL3wwgs66aSTDCeCFzz44IMqLCxU586d1a9fP9NxAMDTKPoAgJiaPn165LrroUOH6p577jEbCJ6yceNGXX311XIcRzfffLOuvfZa05HgAQsXLtTYsWMlVewG4hp5APhlFH0AQMx8/fXX6tu3r8rKynTVVVfpqaee4gk7IoqKitS3b1/t2rVLp556auS0fdjNcRzdddddchxHAwcO1Jlnnmk6EgB4HkUfABAT33//vXr27KmCggJdeOGFGjt2LNfYIsJxHA0ZMkQLFixQVlaWJk6cyMGykCRNmDBBn332mVJTUzVixAjTcQDAF3iGBQCIuq1bt+qiiy7Szp07ddppp+mdd97hNp7Yz+OPP67//ve/SkxM1MSJE9WsWTPTkeABubm5+sMf/iBJuvfee9WkSRPDiQDAHyj6AICoysvL08UXX6z169erRYsWmjp1qmrUqGE6Fjzk/fff1/333y9JGjVqlLp27Wo4Ebzi/vvv19atW9WqVSsNHz7cdBwA8A2KPgAgakpLS9WvXz8tXrxYdevW1UcffaS6deuajgUPWbZsmQYOHCjHcXTrrbfqtttuMx0JHvHVV1/pueeekyQ9//zzSklJMZwIAPyDog8AiArHcXTjjTfq448/VrVq1TRlyhS1aNHCdCx4yPbt23XppZcqPz9fXbp00ahRo0xHgkeUlJTopptukuM4uu6669jlAQBHiKIPAIiK+++/X2+88YYSEhI0ceJE/epXvzIdCR5SUFCgSy65ROvWrVOzZs00ceJEJSUlmY4Fj/jzn/+s7777TtnZ2XriiSdMxwEA36HoAwCq3OjRo/XYY49Jkl588UX9+te/NpwIXuLeXvHrr79WVlaWPvzwQ9WpU8d0LHjEnDlzIuV+zJgxys7ONpwIAPyHog8AqFKTJ0/W73//e0nSo48+qsGDB5sNBE9xHEe///3v9f777yslJUWTJ09W69atTceCRxQUFOi6665TOBzWoEGD1KdPH9ORAMCXKPoAgCqzePHiyMFqt9xyS+QkdcD1xBNP6LnnnlNcXJzefPNNnX322aYjwUN+//vfa82aNTruuOM0cuRI03EAwLco+gCAKrF161b17t1bBQUFuvDCC/XMM88oLi7OdCx4yLhx43TvvfdKkp5++mn169fPcCJ4yRtvvKGXX35Z8fHxeu2115SZmWk6EgD4FkUfAHDMioqK1KdPH23cuFGtWrXS+PHjOVgN+/n8888jl3EMHTpUQ4cONZoH3rJq1SrdeuutkqQ//elPnLIPAMeIog8AOCbubfTmzZunWrVq6YMPPlCtWrVMx4KHLF++XH369FFpaan69u2rJ5980nQkeEhRUZEGDBiggoICde3aVQ8++KDpSADgexR9AMAxGTFihMaNG6fExET973//0wknnGA6Ejxky5Ytuvjii7V371517tw5cstFwDV06FAtWbJE2dnZevPNN/n7AQBVgKIPADhqM2bMiEzfRo8ezXZb7Gffvn3q1auXNmzYoBNOOEGTJk1SWlqa6VjwkBdeeEFjxoxRXFyc3njjDTVo0MB0JAAIBIo+AOCobNy4MXLC/o033qhbbrnFdCR4SCgU0oABA7Ro0SJlZ2frww8/VJ06dUzHgod88cUX+u1vfytJ+tvf/qYePXoYTgQAwUHRBwAcsZKSEvXv3187d+7UqaeeqtGjR5uOBA9xHEe33367pk2bprS0NH3wwQdq0aKF6VjwkA0bNqhfv34qKyvTFVdcofvuu890JAAIFIo+AOCI3XXXXZo/f75q1aqliRMnKjU11XQkeMjf//53vfjii4qPj9dbb72lM844w3QkeEhhYaEuu+wy7dixQx06dNBLL73ErTgBoIpR9AEAR+Tdd9/Vv//978g1tc2aNTMdCR7y5ptvRs5tGDVqlHr37m04EbzEcRwNGTJECxcuVJ06dfTee++pWrVqpmMBQOBQ9AEAhy0nJ0c33XSTJOnuu+9Wz549DSeCl3z22We6/vrrJUl//OMfdccddxhOBK8ZPnx45C4dEydOVNOmTU1HAoBAougDAA5LOBzW4MGDtWvXLnXs2FF//etfTUeCh6xYsUKXXXaZQqGQLr/8cj3++OOmI8FjRo4cGfl7MWbMGHXp0sVwIgAILoo+AOCwjBo1SjNmzFBaWprefPNNJScnm44Ej9i6dasuvvhi7d27V507d9arr76q+HieYuD/vP3227rzzjslVZyw7+78AABEB/8KAwAOaenSpbr33nslSU899ZTatm1rOBG8oqCgQJdccol+/PFHtWzZUpMmTVJaWprpWPCQmTNnatCgQXIcR3fccYeGDx9uOhIABB5FHwDwi8rKynTDDTeotLRUl1xyiW699VbTkeAR5eXluuqqq/TNN9+odu3a+vDDD1WnTh3TseAhixcvVp8+fVRaWqr+/ftr5MiRnLAPADFA0QcA/KJRo0ZpwYIFqlmzpl544QWepCNi+PDhev/995WSkqLJkyerZcuWpiPBQ9atW6eLL75Y+/btU5cuXfT6668rISHBdCwAsAJFHwBwUD/88EPkVmlPPvmkGjZsaDgRvOL111/XE088IUkaO3asOnfubDgRvGTHjh266KKLtHXrVrVv317vvfeeUlNTTccCAGtQ9AEAB+Q4jm655RYVFRWpa9euuvHGG01HgkfMmzcvcpvFBx54QFdeeaXhRPCS/Px89erVS6tXr1aTJk00bdo0ZWZmmo4FAFah6AMADuh///ufPv74Y6WmpmrMmDFs2YckafPmzerTp49KSkr0m9/8Rn/5y19MR4KHhEIhDRgwQF9//bWysrL00UcfsRMIAAyg6AMAfqawsFB33XWXJOnee+/l2mtIkoqKitSnTx9t3bpVJ510kl5//XVuo4cIx3E0ZMgQffjhh0pLS9OUKVPUpk0b07EAwEr86wwA+JknnnhCGzZsUJMmTXTPPfeYjgMPcBxHt99+uxYsWKDatWtr8uTJqlGjhulY8JD7779fr732mhISEjR+/HidddZZpiMBgLUo+gCA/fz444967LHHJFUcwJeenm44Ebzg5Zdf1tixYxUfH6/x48erWbNmpiPBQ0aNGhV53BgzZowuueQSw4kAwG4UfQDAfh544AEVFxfr/PPPV//+/U3HgQcsXrxYd9xxhyTp0Ucf1QUXXGA4Ebzk7bff1tChQyVV/P244YYbzAYCAFD0AQD/Z+nSpfrvf/8rSXr66ac5gA/au3ev+vXrp5KSEvXq1Uv33nuv6UjwkJkzZ2rQoEFyHEd33HGH7r//ftORAACi6AMAKnnooYfkOI4GDBigjh07mo4DwxzH0Q033KAffvhBxx9/vF577TUO30PE4sWL1adPH5WWlqpfv34aOXIkLw4C8B3HcVRQUKBdu3YpFAqZjlNlEk0HAAB4w7x58zRp0iTFx8dzyzRIkl566SW9++67SkpK0oQJE5SVlWU6Ejxi3bp1uvjii7Vv3z516dJFb7zxhhISEkzHAmAht6jv2bNHe/fujbzl5ubu9/Ev/XhZWZkkqWPHjoEZdFD0AcBDHMcx9nu7W24HDx6s1q1bG8thA5PrfLhWr16tP/zhD5Kkv//97/rVr35lOBG8YseOHbrooou0detWtW/fXu+9955SU1NNxwLgc47jKC8vT7t27dLu3buP6K0qJvHx8fHKz8+vgj+JN1D0AcAwL2x1/fLLLzVz5kwlJyfrT3/6k+k4geWFtT4coVBI11xzjQoLC9W1a1cNGzbMdCTf8MsaH63CwkJdeumlWr16tZo0aaJp06YpMzPTdKyYC/o6A8fKcRzt3btXO3bsOODbzp07f1bW9+zZo/Ly8qP+PZOSklSrVi3VrFlTmZmZB3w72M9Vq1ZNn332WaBuC0rRBwDoH//4h6SKaX7Tpk0Np4Fpf/3rXzV//nxlZmbq1Vdf5bp8SJLKy8s1cOBAzZs3T7Vq1dK0adPUsGFD07EAxIDjONqzZ4+2bt160PL+0yLvboc/Uunp6crKyvrFt9q1a+/3ca1atZSenn7UL8KFQqHAvYBH0QcAyy1fvlzvv/++4uLidNddd5mOA8PmzZunv/3tb5Kk559/Xo0bNzacCF7gOI6GDh2qSZMmKSUlRZMmTVLbtm1NxwJwjIqKirRt2zZt3bpVW7Zs0datWyNvP/34aLbH16hRQ9nZ2T97q1Onzs/Keu3atVWrVi0uBaoiFH0AsNyTTz4pSbrsssvUqlUrw2lgUmlpqW688UaFw2ENHDhQV1xxhelI8IinnnpKo0ePliS99tprOvfccw0nAvBLysrKtGXLFm3atGm/t5ycnP1KfG5u7hH9urVq1VLdunX3K+wHKvLuz1HazaHoA4DFNm3apDfeeEOSdM899xhOA9NGjBihZcuWKTs7W6NGjTIdBx7x9ttv6+6775ZU8cLggAEDDCcC7FZSUqKcnJyflfjKb1u3blU4HD6sXy8lJUX169dXgwYNVL9+/cjbTz+uV6+eUlJSovynQ1Wh6AOAxUaNGqVQKKTzzjtPZ555puk4MGjZsmWRLfujRo1S7dq1DSeCF3zxxRcaNGiQJOn3v/89BzMCMVBSUqKNGzdq/fr1kbd169ZF3s/JyTmsXycxMVGNGjXScccdF3lr2LDhfgW+QYMGqlmzZuCuTwdFHwCsVVRUpBdffFGS9Mc//tFwGpgUDod18803KxQK6ZJLLmHLPiRJ69evV9++fVVaWqrLLrtMTz/9NGUAqAKO42jr1q1avXq11qxZEynx7n9zcnIOeRvWlJSU/Qr8gd7q1q3LYaoWo+gDgKXefvtt7dmzR02bNlXPnj1Nx4FBb7zxhubMmaNq1arp2WefpcxB+fn5+s1vfqOdO3fq1FNP1RtvvKGEhATTsQDfcBxH27Zt0+rVqyOFvvL7BQUFv/j16enpOv744w/41rRpU2VnZ/NYjV9E0QcASz333HOSpFtuuYUn8BbLy8uLnM/w0EMPcco+FA6HNWjQIC1ZskT16tXTe++9p/T0dNOxAE8qKyvT2rVrtXz5ci1fvlwrVqzQ8uXLtXr1auXn5x/06+Lj43X88cerZcuWatGixX5FvlmzZqpTpw5FHseEog8AFlq4cKHmz5+vpKQk3XjjjabjwKC//OUv2rZtm0444QQNHTrUdBx4wF/+8he9++67Sk5O1rvvvsuLP4Aq7rO+atWqSKF3S/2qVasOetu5+Ph4NW3aVC1bttQJJ5wQeWvZsqWaNWum5OTkGP8pYBOKPgBY6KWXXpIk9evXT3Xr1jWcBqasWLFCI0eOlCSNHDmS05Sh//3vf3rkkUckSc8//7w6depkOBEQe3v37tW3336rxYsXa+HChfryyy+1adMmlZaWHvDz09PT1bZtW7Vr107t2rVT27Zt1bp1azVr1ozHVRhD0QcAy5SUlGjcuHGSpOuvv95wGpg0fPhwlZWV6ZJLLtHFF19sOg4MW7lypa677jpJ0p133snjA6ywdetWzZ8/XwsXLtTixYu1ePFi/fjjjwf83Bo1aujEE0/cr9C3a9dOTZo04dA7eA5FHwAs88EHH2jPnj1q1KiRunXrZjoODJk7d64mTZqk+Ph4/eMf/zAdB4YVFhbq8ssvV0FBgc4//3z+TiCQ8vPz9c0332j+/PmaN2+e5s+fr40bNx7wc5s2baoOHTropJNOUnl5ua6//nqdcMIJXDcP36DoA4BlXnvtNUnSNddcwyF8lnIcR/fdd58kafDgwWrbtq3hRDDtd7/7nb777jvVq1dP48aNU2IiTxHhb47j6IcfftCsWbP05Zdfav78+Vq+fLnC4fB+nxcXF6cTTzxRp512mjp27KgOHTro5JNPVq1atSRVXJs/depUNWvWjJIPX+FRHAAssnPnTk2dOlWSIlt0YZ9p06Zp1qxZSklJ0cMPP2w6Dgx79dVX9fLLLys+Pl7jxo1T/fr1TUcCjlg4HNby5cs1a9YszZo1S1988YVycnJ+9nmNGzfWGWecoTPOOENnnnmmTj31VNWoUcNAYiC6KPoAYJF3331XZWVl6tixI1NcSzmOo/vvv19SxRSXE9XttnLlSt12222SpIcfflhdu3Y1nAg4PI7jaM2aNZo+fbo+/vhjzZo1S7t3797vc5KSknTGGWfonHPOUadOnXTGGWeoQYMGhhIDsUXRBwCLTJgwQZI0YMAAw0lgypQpU7R48WLVqFEjsn0fdgqFQrr22mtVVFSkCy+8MPICEOBVe/bs0cyZMzV9+nRNnz5d69ev3+/nq1Wrpk6dOum8887TeeedpzPOOENpaWlmwgKGUfQBwBI7d+7UzJkzJUn9+/c3nAYmOI6jv//975Kk2267TbVr1zacCCb9/e9/14IFC1SrVi2NHTuWMzvgOY7jaMWKFZo8ebImT56sefPm7XeNfXJyss4++2z16NFDF1xwgTp27KikpCSDiQHvoOgDgCXee+89lZeXq0OHDmrZsqXpODBg1qxZmjt3rlJSUnTnnXeajgODFixYoL/+9a+SpGeffVaNGjUynAioUFZWpjlz5mjy5MmaNGmS1qxZs9/Pt2vXTj169FCPHj103nnnqVq1aoaSAt5G0QcAS0ycOFGSdPnllxtOAlNGjBghSbrhhhs4cM1iRUVFGjRokMrLy3XFFVfoyiuvNB0JlguHw5o9e7beeustTZgwQTt27Ij8XHJysrp166bevXurV69enCsCHCaKPgBYYN++fZFt+3379jWcBiYsXLhQH330kRISEnT33XebjgODHnroIa1YsUINGjTQv//9b9NxYCnHcbRw4UK99dZbevvtt/e7n31WVpZ69eql3/zmN+rRowen4gNHgaIPABb4+OOPFQqF1KJFC7Vu3dp0HBjw9NNPS5KuvPJKNWvWzHAamLJw4UL985//lCT95z//4ZwGxNz27dv12muv6aWXXtLKlSsjP56RkaHLLrtMV111lbp166bERGoKcCz4DgIAC0yZMkWS1KtXL8XFxRlOg1jbvn175I4LXJtvr/Lyct1yyy0Kh8O68sor1atXL9ORYIny8nJNnz5dL774oiZPnqyysjJJUmpqqi699FJdddVVuvjii5Wammo4KRAcFH0ACDjHcTR16lRJ4om9pV566SWVlpbqjDPO0GmnnWY6Dgx59tlntWDBAtWsWTMy1Qeiae/evfrPf/6j0aNHa8OGDZEfP/PMMzVkyBANGDBAGRkZBhMCwUXRB4CAW7RokbZs2aJq1aqpS5cupuMgxsrLy/X8889Lku644w7DaWBKTk6OHnjgAUkVhzJyGCOiafXq1Ro5cqTGjh2rgoICSRXX3V977bW68cYb1b59e8MJgeCj6ANAwLnT/AsvvFApKSmG0yDWpkyZog0bNigrK0sDBgwwHQeGPPjgg9q3b5/OOOMM3XLLLabjIKDmzZunv/3tb/rggw/kOI4kqX379ho6dKgGDhzI1nwghij6ABBwn3zyiSTp17/+teEkMOHZZ5+VJN144408ybbUwoULNXbsWEnSyJEjFR8fbzYQAmfu3Ll65JFH9NFHH0V+7JJLLtHQoUN1wQUXcDYMYABFHwACrKioSHPmzJEkXXDBBYbTINbWrVunjz76SHFxcbr11ltNx4EBjuPorrvukuM4uuqqq3TWWWeZjoQAmT9/vh588EHNmDFDkpSQkKBBgwbp3nvv5Q4vgGEUfQAIsLlz56q0tFSNGjXSCSecYDoOYuz111+XJHXr1k3Nmzc3nAYmTJ48WZ999plSU1P12GOPmY6DgFi/fr3uv/9+jRs3TpKUmJiowYMHa/jw4TzWAB7hu71bpaWlev3119WzZ081bdpUqampatCggTp37qwnn3xSO3fuNB0RADxj5syZksTWSQs5jhMp+oMGDTKcBiaUlZXp3nvvlSQNGzZMTZo0MZwIfpebm6t7771Xbdq00bhx4xQXF6frrrtOq1ev1n/+8x9KPuAhvpror1y5UldddZUWL168349v3bpVW7du1dy5c/XEE0/olVdeUc+ePc2EBAAPcYt+165dDSdBrH311Vdas2aNqlWrpssuu8x0HBgwbtw4rVq1SllZWZHCDxwNx3H01ltvaejQodq+fbukiheQn3zySXXs2NFwOgAH4puiv2nTJnXr1k05OTmSpLi4OJ133nlq0aKFduzYoY8//lhFRUXavn27+vTpo2nTpnE9KgCr7du3T19//bUkrs+30WuvvSZJ6tevn6pXr244DWKtrKxMf/nLXyRJf/zjH7lXOY7a2rVrddttt0Wuw2/durWeeuop9ezZk51igIf5pugPHDgwUvKbNm2qSZMm6ZRTTon8/M6dO3XllVfqk08+USgU0uWXX661a9cqMzPTUGIAMGvOnDkqKytTs2bN1LRpU9NxEEOhUEgTJkyQJF177bWG08CEN998U2vWrFGdOnX029/+1nQc+FA4HNYzzzyj++67T8XFxUpJSdGDDz6ou+++m1u1Aj7gi2v0p06dqi+++EKSlJycrPfff3+/ki9JderU0aRJkyLXBu3evVv/+Mc/Yp4VALxi7ty5kqRzzjnHcBLE2meffaZdu3YpOztb559/vuk4iLGysjL99a9/lSTdfffdqlGjhuFE8JuNGzeqR48eGjp0qIqLi9WtWzctXbpUDz74ICUf8AlfFP1///vfkfevu+46tW/f/oCfV61atcg2NUl64YUXVFZWFvV8AOBFbtHv1KmT4SSItYkTJ0qSLrvsMiUm+mbzHqrIuHHjtHbtWmVnZ+uOO+4wHQc+M378eLVv316ffPKJ0tPT9eyzz2rGjBncuQXwGc8X/fz8fH3yySeRj6+//vpf/PzK1yLu3r1bs2bNimo+APCicDisefPmSRL3zbZMWVmZ3n33XUlS//79DadBrDmOo6eeekqSdOedd6patWqGE8EvQqGQhg4dqiuuuEK5ubk644wztGjRIt12221ciw/4kOeL/pw5c1RSUiKpYmJ/+umn/+Lnp6am7je9ck+cBgCbrFixQrm5uapWrdpBd0EhmGbNmqUdO3aodu3abNu30GeffaZvv/1W6enpuuWWW0zHgU/k5OSoa9euGjlypCTpvvvu0+zZs9WqVSvDyQAcLc8X/RUrVkTeb9++/WFtQTz11FMP+PUAYIuvvvpKknT66aezddsy77zzjiTpN7/5jZKSkgynQaw9/fTTkqTBgwcrKyvLcBr4wYIFC3Tqqadq9uzZysjI0HvvvacRI0bwbwfgc57/Dl61alXk/cM9NbpJkyaR91euXFnlmQDA67g+306O42jKlCmSKoo+7PL999/rgw8+kCT94Q9/MJwGfjBlyhQNGDBAhYWFat++vd555x21bNnSdCwAVcDzE/1du3ZF3q9Xr95hfU39+vUj7+/evbvKMwGA17lFn+vz7bJy5UqtX79eKSkp6tatm+k4iLF//etfkqRLL72ULdc4pDFjxqh3794qLCxUjx49NHv2bEo+ECCen+jn5+dH3k9LSzusr6n8eZW//kBKSkoiZwBIUl5enqSKA0lCodCRRI05N5/Xc+LYsM7B5zhO5P3S0tJjXuv8/PzIZUunnnoqf3c8pLy8XFLFYYnRWJfJkydLkrp06aLk5GTW3oBwOBx5P5b//3Nzc/Xqq69Kkn73u9+x9lHm/v91HMeX/6//9re/6ZFHHpFUcUerZ599VklJSb78s0Qbz8Ps4Kd1PtyMni/6xcXFkfeTk5MP62sq39+zqKjoFz93xIgRkQe6yqZPn6709PTDTGnWjBkzTEdADLDOdvj888+VmZl5TL/GypUr5TiOsrKy9M0331RNMFSJNWvWSKo4+Grq1KlV/uu/8cYbkiouYYvGr49Dy8nJibwfy8ftDz/8UIWFhTruuONUUFDA+kdZ5R2jfvr32XEcjRs3TuPHj5ckXXHFFerTp4+v/gym8P/IDn5Y58LCwsP6PM8X/dTU1Mj7paWlh/U1lSf0h9oFMHz4cA0bNizycV5enho3bqwePXooIyPjCNPGVigU0owZM9S9e3cOXAow1tkuXbp0UaNGjY7p19iwYYOkim37PXv2rIpYqCLffvutJKlhw4ZVvja5ubmRc2nuuusutWjRokp/fRyet99+O/J+LB+3H374YUnS0KFD1atXr5j8njar/IKOX/59dhxHf/rTnyIl/7HHHtvvOTAOjOdhdvDTOrs70A/F80W/evXqkfcPNZ0/0OdV/voDSUlJ2W8HgCspKcnzi+zyU1YcPdbZDlWxzkuWLJEkdezYkb8zHpOQkCBJio+Pr/K1+eKLL1RWVqZWrVqpTZs2Vfpr4/DFx//f8Uexetz+5ptvtHjxYiUnJ2vw4MF838dA5f/Hfvn3+bHHHtPjjz8uqeLuDHfeeafhRP7il3XGsfHDOh9uPs8fxle7du3I+9u2bTusr9m6dWvkfW4tA8A2ixcvliR16NDBaA7E1qeffiqpYroIu7z44ouSpH79+u33vAlwvfzyyxo+fLgk6cknn6TkAxbwfNFv3bp15P0ff/zxsL7G3bYqiakGAKuEQiEtXbpUUsVEH/aYOXOmJOmCCy4wnASxVFpaGtmKfcMNNxhOAy96//33ddNNN0mS7rnnHt11112GEwGIBc8X/bZt20beX7p0qcrKyg75NQsXLjzg1wNA0K1atUolJSWqUaOGmjVrZjoOYmTbtm1atmyZ4uLi1KVLF9NxEEPTpk3T7t271aBBA3Xt2tV0HHjM0qVLddVVVykcDuv666/XY489ZjoSgBjxfNHv3Llz5Br6goICLViw4Bc/v6SkRF999VXkYyYbAGzibts/5ZRT9rtWGMHmbts/5ZRT2LptmTfffFOSdNVVV0XOgAAkadeuXfrNb36jgoICXXDBBXrhhRcUFxdnOhaAGPH8s8Dq1aurW7dukY/Hjh37i5//zjvvaN++fZIqrs8/77zzohkPADxl0aJFkrg+3zZu0efFbbvk5eVp8uTJkqSBAwcaTgMvKSsr04ABA7Ru3To1a9ZM48eP9/wBYwCqlueLviTdfvvtkffHjh2rZcuWHfDzCgsL9ac//Sny8c0336zERM/fWAAAqox7ff4pp5xiOAliievz7TRp0iQVFxerdevWOvXUU03HgYc88sgjmjlzpqpVq6ZJkyax0wewkC+Kfq9evXTuuedKqtiaf8kll0RuH+XatWuX+vTpozVr1kiqmObfe++9Mc8KACatWLFCktSuXTvDSRAr27Zt05o1axQXF6dzzjnHdBzE0P/+9z9J0oABA9iSjYjPP/9cf/vb3yRV3JGhffv2hhMBMME34+7//ve/OuOMM7RlyxatX79eHTp0UJcuXdSiRQvt2LFDH3/8sQoLCyVJiYmJGj9+vDIzM82GBoAY2rdvnzZt2iSJO47YZO7cuZKkk046STVr1jScBrFSUFCgjz76SJLUt29fw2ngFbt379Y111wjx3E0ePBgXXnllaYjATDEN0X/uOOO08yZM3XVVVdp8eLFchxHn332mT777LP9Pi87O1uvvPLKftf1A4ANVq1aJUmqW7eusrKyDKdBrLhFv1OnToaTIJY++ugjFRcX6/jjj+dSHUTceuut2rRpk0444QQ988wzpuMAMMg3RV+qmFDNmzdPb731lsaNG6dly5Zp27ZtyszMVPPmzdW3b19df/31qlOnjumoABBzK1eulMRtRW1D0bfTu+++K0m67LLL2LYPSRVnNkyYMEEJCQkaN26cqlevbjoSAIN8VfQlKTk5WYMGDdKgQYNMRwEAT3Gvz2fbvj1KS0v19ddfS6Lo2yQUCun999+XVFH0gby8PN1xxx2SpLvvvlunnXaa4UQATPPFYXwAgENjom+fb7/9VsXFxcrKylKrVq1Mx0GMzJ49W7m5uapTp446d+5sOg484IEHHtDmzZvVokWL/e5ABcBeFH0ACAgm+vZxt+2fddZZbN+2iHsI30UXXaSEhATDaWDa/Pnz9e9//1uS9MILLygtLc1wIgBeQNEHgAAIhUKR24sy0bfHV199JYlt+7aZNm2aJOnXv/614SQwzXEcDRs2TI7j6Nprr+UwagARFH0ACIAffvhBoVBI6enpOu6440zHQYwsXLhQkvSrX/3KcBLEypYtW7R48WJJUo8ePcyGgXHvvfeeZs+erbS0NI0YMcJ0HAAeQtEHgABwr89v06aN4uN5aLdBQUGBvv/+e0lSx44dDadBrEyfPl2SdNppp6lu3bqG08CkUCike++9V5J01113qVGjRoYTAfASng0CQACsXr1akjiQzSJLliyR4zhq0KCB6tWrZzoOYsS9Pp9t+xgzZoxWr16tunXr6p577jEdB4DHUPQBIADWrVsnSWrevLnhJIiVRYsWSWKabxPHcfTJJ59IYtu+7QoKCvTII49Ikh5++GHVqFHDcCIAXkPRB4AAcIt+s2bNDCdBrLhFv0OHDmaDIGZWrVql7du3KzU1VWeeeabpODBozJgx2rFjh5o3b64hQ4aYjgPAgyj6ABAAFH37uAeyMdG3x+effy6p4naKKSkphtPAlOLiYj3xxBOSpOHDhyspKclwIgBeRNEHAJ8Lh8MUfcuEQiEtXbpUEkXfJm7R79Kli+EkMOnll1/Wli1b1LhxYw0aNMh0HAAeRdEHAJ/bunWrSkpKFB8fr8aNG5uOgxhYuXKlSkpKlJGRwYs7lnAch6IPlZeX6+mnn5Yk3XPPPUpOTjacCIBXUfQBwOfcaX6TJk3YwmkJd9v+Kaecwu0ULbF27Vrl5OQoKSlJZ511luk4MGTKlClau3atMjMzdf3115uOA8DDeHYAAD7Htn37LF++XJJ00kknGU6CWJk1a5Yk6YwzzlBaWprhNDBl5MiRkqSbb75Z1apVM5wGgJdR9AHA5yj69lmxYoUkqW3btoaTIFbmzJkjSTr33HMNJ4EpS5cu1cyZM5WQkKA77rjDdBwAHkfRBwCfo+jbxy367dq1M5wEsTJ//nxJYtu+xV544QVJ0mWXXaYmTZoYTgPA6yj6AOBzP/zwgySpefPmhpMgFkpLS7V27VpJTPRtkZ+fr2XLlkmq2LoP+xQVFemNN96QVLFtHwAOhaIPAD7HRN8uq1evVnl5uTIyMtSgQQPTcRAD33zzjcLhsBo3bsyaW+qdd95Rbm6umjZtqm7dupmOA8AHKPoA4GOhUEibNm2SRNG3ReXr8+Pi4gynQSy42/aZ5tvrpZdekiTdcMMN3GkDwGHhkQIAfGzTpk0Kh8NKSUlRvXr1TMdBDHAQn30o+nZbv369Pv30U8XFxXFLPQCHjaIPAD62efNmSVKjRo2Y7lqCom8fir7d3n77bUnSBRdcoMaNGxtOA8AvKPoA4GM5OTmSKoo+7EDRt8vWrVu1YcMGxcXF6bTTTjMdBwa89dZbkqQrr7zScBIAfkLRBwAfqzzRR/CFw2GtWrVKEkXfFgsXLpQktWnTRjVq1DCcBrG2cuVKLV68WImJierbt6/pOAB8hKIPAD7mTvQbNmxoOAliYePGjSoqKlJycjKHL1piyZIlkqQOHTqYDQIj3Gl+jx49lJWVZTgNAD+h6AOAjzHRt8sPP/wgqeIOCwkJCYbTIBbcon/yyScbTgITJkyYIEm64oorDCcB4DcUfQDwMbfoM9G3Q+WiDztQ9O21du1aLV++XImJierdu7fpOAB8hqIPAD7GYXx2cYt+8+bNDSdBLJSUlGjlypWSKPo2ev/99yVJ5557rjIzM82GAeA7FH0A8CnHcdi6b5l169ZJoujbYsWKFSovL1etWrX4HreQW/QvvfRSw0kA+BFFHwA8xHGcw/7cvXv3qqioSJLUoEGDaEVCFBzJOlfGRN8/jnaNK6u8bT8uLu6Yfz1UvapY5wPJzc3VrFmzJFH0ARwdij4A+JS7bT8rK0tpaWmG0+BwHGtZo+h7X1UW8m+//VYS2/a9KNovvEybNk1lZWVq06aNWrZsGdXfC0AwUfQBwKc4iM8u+/bt044dOyRxGJ8tOIjPXtOmTZMkXXLJJYaTAPArij4A+BQH8dnFvT6/du3aysjIMJwGsUDRt5PjOPrkk08kSd27dzecBoBfUfQBwKeY6NuFbft22bNnj7Zv3y5Jatu2reE0iKU1a9Zo48aNSk5O1jnnnGM6DgCfougDgE9x4r5dKPp2+f777yVVvJBXo0YNw2kQS+40v1OnTkpPTzecBoBfUfQBwKfYum8Xbq1nl1WrVkmSWrVqZTgJYs0t+t26dTOcBICfUfQBwKfcos/WfTsw0beLO9Fv3bq14SSIpXA4rJkzZ0qi6AM4NhR9APCpbdu2SZLq1atnOAliYf369ZKk448/3mgOxAYTfTstWbJEu3fvVvXq1XX66aebjgPAxyj6AOBTO3fulCRlZ2cbToJY2LRpkyTpuOOOM5wEseBO9Cn6dpk9e7YkqXPnzkpKSjKcBoCfUfQBwIcKCwtVWFgoSapTp47hNIi2/Px85eXlSeJMBhuEw2GtXr1aElv3bTNnzhxJFUUfAI4FRR8AfGjHjh2SpOTkZE7ktoB7h4WMjAzW2wKbNm1SUVGREhMTuVTDMnPnzpVE0Qdw7Cj6AOBDlbftx8XFGU6DaHO37TPNt4O7bb9FixZs37bI1q1btW7dOsXFxenMM880HQeAz1H0AcCH3Ik+2/bt4E70Kfp24CA+O7nT/JNOOkkZGRmG0wDwO4o+APiQW/Q5iM8ObtHnID47cGs9O7nX53fq1MlwEgBBQNEHAB/ixH27MNG3y5o1ayRJLVu2NJwEseRO9Cn6AKoCRR8AfIit+3bhGn27/Pjjj5KkZs2aGU6CWCkrK9M333wjiaIPoGpQ9AHAh5jo24WJvj0cx9H69eslSU2bNjUbBjGzcuVKFRcXq0aNGjrhhBNMxwEQABR9APAhJvp2oejbY/fu3SooKJAkNWnSxHAaxMqiRYskSaeccori43l6DuDY8UgCAD7EYXz2CIVC2rp1qyQO47OBu22/Xr16SktLM5wGsbJ48WJJUocOHYzmABAcFH0A8CG27ttj69atchxHSUlJrLcF2LZvJ3ei37FjR8NJAAQFRR8AfIit+/Zwt+03aNCALb0WcCf6xx9/vNkgiBnHcSITfYo+gKrCMwYA8JmysjLt3r1bEhN9G3B9vl2Y6Ntnw4YN2rNnjxITE9WuXTvTcQAEBEUfAHzGLflxcXHKysoynAbR5t5aj+vz7cBE3z7uNP/EE09USkqK2TAAAoOiDwA+427br1WrlhITEw2nQbRt2bJFktSwYUPDSRALbtFnom8Prs8HEA0UfQDwGU7ct8v27dslSXXr1jWcBLHgbt1nom+PJUuWSOLEfQBVi6IPAD7jnrjPQXx2cF/YoegHX15envbu3SuJib5NVqxYIUlcnw+gSlH0AcBnmOjbxZ3os97B527br127tqpXr244DWIhFAppzZo1kqS2bdsaTgMgSCj6AOAz7kSf4mcHXtixByfu22ft2rUqKytTtWrVuLMGgCpF0QcAn9m1a5ckceK+JSj69nDvsNC4cWPDSRArK1eulCS1adNGcXFxhtMACBKKPgD4zJ49eyRVnLqPYCsqKlJ+fr4krtG3QU5OjiTusGAT9/p8tu0DqGoUfQDwGfewLop+8LnT/KSkJGVkZBhOg2hzb6XYoEEDw0kQK5Un+gBQlSj6AOAzTPTtUXnbPtt6g4+ibx8m+gCihaIPAD5D0bcH1+fbhaJvF8dxmOgDiBqKPgD4jLt1PzMz02gORJ97az2uz7cDRd8uOTk52rdvnxISEtSyZUvTcQAEDEUfAHyGib49mOjbo7y8PPLCDkXfDu40v0WLFkpOTjacBkDQUPQBwEdCoZAKCgokUfRtQNG3x/bt2xUOhxUfH88ODkt8//33kqRWrVoZTgIgiCj6AOAj7jRfkmrWrGkwCWKBrfv2cLft161bVwkJCYbTIBbWrVsnSWrevLnhJACCiKIPAD7iXp+fkZFBGbAAE317cH2+fSj6AKKJog8APsL1+Xah6NuDom8ft+g3a9bMcBIAQUTRBwAfoejbhaJvj5ycHEkUfZtQ9AFEE0UfADzEcZxf/HmKvl24Rt8eTPTtkpubq927d0ui6AOIDoo+AHhAXFzcYX1ebm6uJA7is0FxcbHy8/MlMdG3gVv0GzZsaDgJDsfhPmYfjDvNr1OnjqpXr14VkQBgPxR9APCRvLw8SRR9G7jb9pOSklhvCzDRtwvb9gFEG0UfAHzELfoZGRmGkyDa3G372dnZxzw9hPdR9O1C0QcQbRR9APARtu7bw71+t3bt2oaTINocx9HWrVslUfRtQdEHEG0UfQDwESb69nCLflZWluEkiLbdu3crFApJkurVq2c4DWKBog8g2ij6AOAjFH17cIcFe7jnMWRkZCglJcVwGsQCRR9AtFH0AcBH2LpvD3eiT9EPvl27dkmqOIEdwec4jtavXy+Jog8geij6AOAjTPTt4U702boffG7R5zwGO+zYsUOFhYWKi4tTkyZNTMcBEFAUfQDwEYq+Pdi6b4+dO3dKYqJvi5ycHElS3bp1uVQDQNRQ9AHAR9i6bw+27tuDib5dNm/eLElq2LCh4SQAgoyiDwA+wkTfHmzdtwdF3y7uRJ+iDyCaKPoA4BOlpaUqLi6WRNG3AVv37cHWfbu4Rb9Ro0aGkwAIMoo+APjEvn37Iu/XqFHDYBLEgrt1n4l+8DHRtwtb9wHEAkUfAHzCvT6/WrVqSkxMNJwG0cZE3x4UfbuwdR9ALFD0AcAnuD7fHqFQSPn5+ZIo+jZg675d2LoPIBYo+gDgE27RZ9t+8O3duzfyfmZmprEciA0m+nZh6z6AWKDoA4BPuBNein7wuUW/Ro0aSkhIMBsGUeU4TqToM9EPvlAopO3bt0ui6AOILoo+APiEW/SrV69uOAmizT2PoWbNmoaTINpyc3NVXl4uiYm+DbZs2SJJSkpK4oUdAFFF0QcAn2Cibw/3Mg2KfvC50/z09HSlpqYaToNoc6/Pb9CggeLjeRoOIHp4hAEAn2Cibw8m+vZg275dOHEfQKxQ9AHAJyj69nCLPndYCD73xH227duBg/gAxApFHwB8gqJvDyb69uDEfbtwaz0AsULRBwCfoOjbg6JvD7bu24Wt+wBihaIPAD5B0bcHRd8ebN23C1v3AcQKRR8AfIKibw+Kvj2Y6Ntl27ZtkqT69esbTgIg6Cj6AOATFH17UPTtwTX6dtmxY4ckKTs723ASAEFH0QcAn6Do24Oibw+27tvDcZzIelP0AUQbRR8AfIKib4+8vDxJFH0b7N69W5KUlZVlOAmiLTc3V+Xl5ZK4VANA9FH0AcAnKPr2YKJvD3eta9WqZTgJos3dtl+9enWlpqYaTgMg6BJNBzhc69ev14wZM/T5559r6dKl2rBhg/Lz81WjRg0dd9xx6tSpkwYOHKguXbqYjgoAUUHRtwdF3x579+6VxFrbwN22zzQfQCx4vugvWrRIt956q+bPn3/An9+zZ4/27NmjpUuXasyYMTr//PP16quvqkmTJjFOCgDRRdG3h7t1PyMjw3ASRFM4HI6sdWZmptkwiDqKPoBY8nzRX7Vq1c9KfqtWrXTSSSepTp062rt3r+bMmaNNmzZJkj777DN16tRJX3zxhZo3b24iMgBUOcdxKPqWCIfDKigokCTVqFHDcBpEU35+vsLhsCSKvg04cR9ALHm+6LtatmypIUOG6JprrlGjRo32+7lwOKyxY8fqd7/7nQoLC5WTk6Orr75ac+bMUVxcnKHEAFB1ioqK5DiOJIp+0Lkv6EisddC52/aTk5O5ZtsCTPQBxJLni36DBg30yiuv6Nprr1VCQsIBPyc+Pl433HCDatWqpb59+0qSvvrqK02fPl0XXXRRLOMCQFS45S8uLk5paWmG0yCa3LVOSEig/AWcexYD03w7UPQBxJLnT93v0qWLBg8efNCSX9lll12mM844I/LxlClTohkNAGKmsLBQkpSenq74eM8/dOMY7Nu3T1LFtn12pQUbB/HZha37AGIpcM8Wzz777Mj769evNxcEAKqQW/SZ5gcfZzHYg4m+XZjoA4ilwBX9ytOP8vJyg0kAoOpUnugj2CpP9BFs7kSfom8Ht+gz0QcQC4Er+kuXLo2837hxY4NJAKDqUPTt4RZ9JvrB50702bpvB3frPhN9ALEQqKK/YcMGzZw5M/LxhRdeaDANAFQdir493K37TPSDj4m+Xdi6DyCWPH/q/pEYNmxYZLt+kyZNdOmllx7ya0pKSlRSUhL5OC8vT5IUCoUUCoWiE7SKuPm8nhPHhnW2y8Eee9zHprS0NP4u+Jj7b1Q4HD7oOrrlLz09nbX2oXA4HHn/UOu3e/duSRUv6rDW/lJ5vQ5n7UKh0H5nMrDe/sHzMDv4aZ0PN2Ngiv6rr76q//3vf5GPR4wYoZSUlEN+3YgRI/TII4/87MenT5/um8nZjBkzTEdADLDOdvj888+VlZX1sx+fO3euJKmgoEBTp06NdSxUkTVr1kiSNm/efNB1/PrrryVVbOFnrf1n8+bNkiTHcQ75uO1ebrht2zbW2mf27Nkj6fDWWfq/F3Xi4+M1Z86cw7qbFLyF52F28MM6u7s8DyUQRX/BggW69dZbIx9fddVVGjhw4GF97fDhwzVs2LDIx3l5eWrcuLF69OihjIyMKs9alUKhkGbMmKHu3bsrKSnJdBxECetsh7i4ODmOo/POO09NmjT52c/n5ORIqtit1LNnz1jHQxVxi13Dhg0Puo7ffPONJKl169astQ9NmDAh8v6hHrffeOMNSdLpp5/OWvvMtm3bIu8fzr/P7vd+VlbWYe04hXfwPMwOflpnd5fnoRx10R89erRGjx59tF9+QI8++qj69+9/RF+zbt06XXrppSouLpYknXzyyXr++ecP++tTUlIOOPlPSkry/CK7/JQVR491tsPB1tm9xKh69er8PfAxd4oXHx9/0HV0X6nPyMhgrX0oPv7/jj861OO2ex5DZmYma+0zldfrcP59drftZ2dns9Y+xfMwO/hhnQ8331EX/Z07d2rVqlVH++UH5F6XeLi2bNmi7t27a+vWrZKk5s2ba9q0aZ6fxAPAkeIwPntwGJ89uJWiPThxH0Cs+fbU/V27dql79+5au3atJKlBgwb6+OOP1aBBA8PJAKDqUfTtwe317OGuNQOK4HNP3M/OzjacBIAtjrroP/zww3Icp0rfhgwZcli/d15eni666CItW7ZMUsWrox9//LGaNWt2tH8cAPA0ir49mOjbg4m+PdyiX7t2bcNJANjCdxP9goIC9ezZM3JYUc2aNTVt2jS1a9fOcDIAiB6Kvj2Y6NvDPVCJoh987in9B7qrCgBEg6+KfnFxsXr37q3Zs2dLqnjCO2XKFJ122mmGkwFAdFH07eFO9Cn6wcdE3x5u0a9Vq5bhJABs4ZuiHwqF1K9fP82cOVNSxWn5kyZN0tlnn204GQBEH0XfHgUFBZIo+kFXWlqq0tJSSRR9G7gHTmdmZhrNAcAevij65eXlGjhwoKZOnSpJSkxM1Pjx43XhhRcaTgYAsUHRt4db9FnrYHOn+RJF3wZu0WeiDyBWPF/0HcfRjTfeqIkTJ0qquD/t66+/rt69extOBgCxQ9G3h7vW1apVM5wE0eRen5+WlqbExKO+2zF8wt26z0QfQKx4/l+W5557Tq+++mrk4xYtWujLL7/Ul19+eVhfP3r06GhFA4CYoejbg4m+Hbg+3y5M9AHEmueL/vbt2/f7ePXq1Vq9evVhfz1FH0AQuEU/LS3NcBJEk+M4TPQtQdG3CxN9ALHm+a37AAAm+rYoKiqKvM9aBxtF3x5lZWWR9aboA4gVzxf9hx9+WI7jHPUbAAQBRd8O7jpLrHXQudfoZ2RkGE6CaHPXWqLoA4gdzxd9AABb923hXp+fkpKihIQEw2kQTUz07eFu269WrZqSkpIMpwFgC4o+APhAcXGxJIp+0HF9vj0o+vbgID4AJlD0AcDjysvLFQqFJFH0g44T9+1B0bcHB/EBMIGiDwAeV1JSEnk/NTXVYBJEGxN9e+Tn50uSqlevbjgJoo2JPgATKPoA4HHutn2p4tptBJc70afoBx9rbY/c3FxJHLwIILYo+gDgce4t1xITE5WYmGg4DaKJuyvYg6JvD/fU/Zo1axpOAsAmFH0A8Dh3os+2/eCj/NmDtbYHt1IEYAJFHwA8jqJvDyb69qDo24OiD8AEij4AeBxF3x6UP3twGJ89KPoATKDoA4DHUfTtwUTfHryoYw+KPgATKPoA4HEUfXtQ/uzBWtuDog/ABIo+AHicW/TT0tIMJ0G0ueWPiX7wUfTtQdEHYAJFHwA8jom+Pdyt+5S/4KPo24OiD8AEij4AeBxF3x5u0Wf3RvBR9O1B0QdgAkUfADyuqKhIEkXfBu5as3U/2MrLy1VSUiKJom8Dij4AEyj6AOBxTPTtwYs6dnCn+RJFP+gcx6HoAzCCog8AHkfRtwcHL9rBLfpxcXF8XwdcYWGhwuGwJIo+gNii6AOAx1H07eFO9Cn6weYW/erVqysuLs5wGkSTO82Pj4/nkhwAMUXRBwCPo+jbg7W2Awfx2aPytn1e1AEQSxR9APA4yp89mOjbIT8/XxJF3wZu0a9Ro4bhJABsQ9EHAI/jum17cBifHZjo28N9Uad69eqGkwCwDUUfADyOib49eFHHDhR9e1D0AZhC0QcAj6Po24Ot+3YoLCyURNG3QeWDFwEglij6AOBxFH07OI7DWlvCLfq8oBN8TPQBmELRBwCP47ptO4RCocj9timAwcbODXtQ9AGYQtEHAI9jymsHt/xJrHXQuWvNfdWDjzssADCFog8AHuI4zs9+jKJvB3ed4+LilJKSYjgNoomJfnAc6DG7Mq7RB2AKRR8APCAuLu6gP+cWQMpfsFW+ROOX/j7A/yj6/ne436Ns3QdgCkUfADyutLRUEhP9oGPnhj04jM8eFH0AplD0AcDjSkpKJEnJycmGkyCamPLag7W2B9foAzCFog8AHudO9Cn6wcbdFexB0bcHE30AplD0AcDj3Ik+1+gHm7t1n/IXfBR9e3AYHwBTKPoA4HFM9O1A+bMHt9ezBxN9AKZQ9AHA49yiz0Q/2DiMzx68qGMPij4AUyj6AOBxHMZnB8qfPVhre3AYHwBTKPoA4GGO4ygUCklioh90HMZnD4q+PbhGH4ApFH0A8DB3277ERD/oOIzPHoWFhZJY66BzHIet+wCMoegDgIdR9O3BlNcerLUdiouLFQ6HJVH0AcQeRR8APMy9Pl+i6Acdh/HZg6JvB3fbvsQdFgDEHkUfADzMnegnJCQoISHBcBpEE0XfHhR9O7iXaKSkpPD4DSDmKPoA4GHcWs8e7u4N1jrYHMeJvKjDlDfYeEEHgEkUfQDwMG6tZw/W2g5uyZcogEHnTvR5QQeACRR9APAwJvr2YKJvB3fKK1H0g46iD8Akij4AeBhTXnvwoo4d3KKfmJioxMREw2kQTdxGEYBJFH0A8DC3/FH0g4+Jvh0of/ZwX9Rhog/ABIo+AHgY5c8erLUdOKDNHmzdB2ASRR8APIyJvj0o+nZwD+Oj6AcfRR+ASRR9APAwrtu2B0XfDm7RT01NNZwE0cbuDQAmUfQBwMM4jM8eFH07sM72YKIPwCSKPgB4GFv37cGLOnZgom8Pij4Akyj6AOBhTP/swVrbgXW2B0UfgEkUfQDwMCb69uA8BjtQ9O3BNfoATKLoA4CHUf7sQQG0A1v37cFEH4BJFH0A8DCu27YHRd8OrLM9KPoATKLoA4CHMdG3BwXQDqyzPdi6D8Akij4AeBgTfXtQAO3A1n17MNEHYBJFHwA8jMP47MGLOnbgBR17UPQBmETRBwAPoxTYIRwOq6ysTBJrHXR8T9uDog/AJIo+AHgYE307uOssUQCDjq379uAafQAmUfQBwMM4jM8O7pRXYq2Djom+PZjoAzCJog8AHsZ123aoXPRZ62BzJ/oU/eBziz4TfQAmUPQBwMPYum+Hyi/oxMXFGU6DaHLXmq37wee+qEPRB2ACRR8APMwt+klJSYaTIJrYzm0P1toenMcAwCSKPgB4mHsSO0U/2Ch/9mDrvh3C4bBCoZAkij4AMyj6AOBh7hNFin6wcRaDPdi6b4fK526w1gBMoOgDgIdR9O3A3RXswe4NO7g7NySKPgAzKPoA4GEUfTtQ/uzB1n07uOuckJCgxMREw2kA2IiiDwAe4jjOfh9T9IPpp+tM0bcHW/ftwAs6AEyj6AOABxzslmoU/WA52DpT9IPjULdHZK2D4VDrzIn7AEyj6AOAh1H07UD5sweTXjtQ9AGYRtEHAA+j6NuBU/ftwdZ9O1D0AZhG0QcAD6Po24FT9+3B7g078IIOANMo+gDgYRR9O7DO9mDSawfWGYBpFH0A8DAKoB1YZ3sw0bcDRR+AaRR9APAwCqAd3HXmGv1gKysrUzgclkTRDzoOXQRgGkUfADzMLYCJiYmGkyCa3Gv0eUEn2NzyJzHpDTom+gBMo+gDgIcx0bcD62wH9wUdid0bQUfRB2AaRR8APKysrEwSBTDo2LpvB7fox8fHs0sn4Cj6AEyj6AOAhzHptQNb9+3AOtuDog/ANIo+AHiU4zhM9C3BCzp2cIs+OzeCz727AkUfgCkUfQDwKLfkSxTAoGPrvh1YZ3sw0QdgGkUfADzKLQUSRT/o2NJtByb69qDoAzCNog8AHkXRtwdb9+3ACzr2cIt+SkqK4SQAbEXRBwCPoujbg6JvByb69mCiD8A0ij4AeJRb/uLj4xUfz8N1kFEA7cA1+vag6AMwjWeOAOBRTHntwVrbgRd07EHRB2AaRR8APIryZw/W2g5co28Pbq8HwDSKPgB4FOXPHkx67cA624OJPgDTKPoA4FEUfXuw1nbgGn17cOo+ANMo+gDgUZQ/e7DWdmCibw93rSn6AEyh6AOAR7nlLzEx0XASRBsF0A5co28P9xp9vqcBmELRBwCPKisrk0QpsAETfTvwgo49WGsAplH0AcCjKH/2YK3twDX69mDrPgDTKPoA4FGUP3sw/bMD62wP1hqAaRR9APAoir49WGs7cI2+PSj6AEyj6AOAR1H+7MFa24HyZw8O4wNgGkUfADyK8mcPCqAduEbfHnxPAzAtEEV/2LBhiouLi7wdf/zxpiMBwDGj6NuDtbYD5c8OjuNwGB8A43xf9OfPn6+RI0eajgEAVY7yZw/W2g5co2+H8vJyOY4jiRd1AJjj66IfCoU0ZMgQhcNh01EAoMpR/uzBpNcOrLMd3HWWWGsA5vi66D/++ONaunSpJGngwIGG0wBA1aLo24O1tgPX6NvBPYhPYq0BmOPbor9y5Uo9+uijkqSrr75a3bt3N5wIAKoW5c8OjuOorKxMEmsddEz07VB5op+YmGgwCQCb+bLoO46jIUOGqKSkRLVq1dLTTz9tOhIAVDmKvh3cdZYogEHHNfp2qHwQX1xcnOE0AGzly6L/3HPPafbs2ZKkJ554QnXr1jWcCACqHkXfDpWLPmsdbEz07cA6A/AC3xX9jRs36r777pMknXvuubrhhhsMJwKA6HALIFs/g42ibw+u0bcDRR+AF/iu6N9+++3at2+fkpOT9cILL7AlCkBgcd22HSpfz8taBxsF0A7uYXysMwCTfFX033rrLX3wwQeSpHvvvVdt27Y1nAgAosct+kz0g63yzg1evA42rtG3Ay/oAPAC3xT9Xbt26fe//70kqVWrVnrggQcMJwKA6GKibwfOYrAHBdAOlQ/jAwBTfDMmuvPOO7Vjxw5J0vPPP19lD54lJSX73e80Ly9PUsUTr8rXTXqRm8/rOXFsWGe7VH7sqbylm/UPhvLycklSOByOrGlBQYGkivLHOvtfOByOvP/T9XS/p+Pj41lrn6u8fpUfqyWpsLBQUsWLd6yz//E8zA5+WufDzeiLoj99+nS9/vrrkqTrrrtOXbt2rbJfe8SIEXrkkUcO+Hump6dX2e8TTTNmzDAdATHAOtth1qxZWr58uSRpzZo1kqQff/xRU6dONRkLVWT16tWSpC1btkTWdMOGDZIqbh3LOvvf5s2bI+//9HF7586dkqTFixfHMhKiwB0MSRXrXPmym4ULF0qSioqK+J4OEJ6H2cEP6+y+mHgoR130R48erdGjRx/tlx/Qo48+qv79++/3YwUFBbrlllskSbVr19aTTz5Zpb/n8OHDNWzYsMjHeXl5aty4sXr06KGMjIwq/b2qWigU0owZM9S9e3e2fAYY62wH90nieeedp+OPP15SxQuOktS6dWv17NnTVDRUoe+++06S1KBBg8iauqWvWrVqrHMATJw4MfL+Tx+33bsGnX322erSpUvMs6HquC/aSBXrXPlyDPeyq+zsbL6nA4DnYXbw0zpXfqHxlxx10d+5c6dWrVp1tF9+QHv37v3Zjz3wwANav369JOmpp55SnTp1qvT3TElJOeBlAElJSZ5fZJefsuLosc52SExMjKyzuwU4OTmZtQ+IhIQESRVbt3+6pnyPB0N8/P8df/TTNXW3W6anp7PWPld5/X66zu5jd0pKCuscIDxG28EP63y4+Tx9GN/ChQv1zDPPSJK6du2q6667znAiAIgd93puTt0PNu6tbg/W2g4cxgfAC4762ePDDz+shx9+uAqj/NySJUsir4pu2LBBZ5111kE/1z2oT6q49rHy5z700EPq1atX9IICQBS42z/dKTCCidso2oNT9+3AOgPwAt88q1i7dq3Wrl17WJ9bWlqqefPmRT6u/CIAAPgFBdAOrLM93ALo9W2hODYUfQBe4Omt+wBgM7bu24GdG/Zw15qiH2zubZsp+gBM8nTRHzx4sBzHOay3V155JfJ1TZs23e/nBg8ebO4PAQBHiQJoByb69nCv0Wetg42JPgAv8HTRBwCbMdG3A0XfHqy1HTiMD4AXUPQBwKMoBXZgne3B1n07MNEH4AUUfQDwKLbu24Gib4fy8nI5jiOJtQ46ij4AL6DoA4BHsXXfDhR9O7jrLLHWQcdhfAC8gKIPAB7FRN8OFH07VC76bN0PNib6ALwgMEW/8gn969evNx0HAI4ZBdAO7Nywg3vivsRaBx2H8QHwgsAUfQAIGgqgHXhBxw5s3beHW/TZuQHAJIo+AHgUW/ftQNG3g7vO8fHxio/n6VeQsXUfgBfwLw0AeBQTfTtQ9O3gbt1nnYPPXWsm+gBMougDgEdRAO3AOtuBdbYHaw3ACyj6AOBRbN23A6XADu46M+UNPib6ALyAog8AHsXWfTtQ9O3A1n17UPQBeAFFHwA8igJoB9bZDqyzPSj6ALyAog8AHsXWfTtQAO3A1n178D0NwAso+gDgUWzdtwOlwA5s3bcHE30AXkDRBwCPYqJvB4q+HVhne1D0AXgBRR8APIpiYAfW2Q5s3bcHaw3ACyj6AOBRbN23A0XfDmzdtwdrDcALKPoA4FFs3bcDRd8OrLM92LoPwAso+gDgUUz07UABtAPbue1B0QfgBRR9APAoCqAdWGc7sJ3bHnxPA/ACij4AeBRb9+1AKbADE317MNEH4AUUfQDwKLbu24GibwfW2R4UfQBeQNEHAA9yHCdS9JnoBxsF0A5s3bcHRR+AF1D0AcCD3JIvUQyCjqJvB7bu24PvaQBeQNEHAA+i6NuDUmAH1tkeTPQBeAFFHwA8yC0FElv3g44CaAe27tvBcRyKPgBPoOgDgAcx0bcHRd8ObN23Q+XHbtYagEkUfQDwoMoTfQpgsFH07cA624HHbgBeQdEHAA+q/GQxPp6H6iCjANqBrft2cNdZYqIPwCyePQKAB1W+tV5cXJzhNIgmir4d2LpvB4o+AK+g6AOAB1H+7MFa24F1tgMHqQLwCoo+AHiI4ziS/u/JIk8Ug8ldZ4kCGFSV11hi636QVV7ryuvMbiwAJlH0AcADfvqE0N26TykIlgM98afoB8vByh1b94PlYOvMrfUAeAVFHwA8iIm+PSj6dmCd7UDRB+AVFH0A8CBKgT1Yazuwdd8OfD8D8AqKPgB4EFv37UExsANb9+3ARB+AV1D0AcCD2LpvD4q+HVhnO1D0AXgFRR8APIiJvj0ogHZg674dKPoAvIKiDwAeRPmzB2ttB7bu24HvZwBeQdEHAA9i6749WGs7/L/27j8oqzL///gLUcmfiFCiSaaUYpOuU4OJlULq5u8xtE2pNn/NbspuO7XuuDPOrs5Os/341G7bspVTG2SZZq6blsYWsZpJoWWGormJP1AxDVBECFC4v38w9/lyiwIicK77vp6PmTNzzs057rt574Hrfd7XfR06+nagow/AFBT6AGAgpu7bgw6gHcizHSj0AZiCQh8ADESX1x4UgHZg6r4dKPQBmIJCHwAMRPFnD3JtB6bu24H7GYApKPQBwEBM3bcHhYEdyLMd6OgDMAWFPgAYiKn7dvB4PKqpqZFEARjomLpvBwp9AKag0AcAA9HRt4M3zxK5DnRM3bcDD3QAmIJCHwAMxDRfO3jzLJHrQMc9bQce6AAwBYU+ABiIqft2oNC3B4W+HZi6D8AUFPoAYCDvlG4K/cBGoW8P7mk7UOgDMAWFPgAYiKLADnULfXId2Lin7cB39AGYgkIfAAzkXYmdoiCw1V2Mr107/iQHMgp9O/AdfQCmYFQBAAbyFgUUf4GN4s8e5NoOTN0HYApGkABgIDr6dqD4swe5tgOLLgIwBYU+ABiIjr4deKBjDwpAO3h/d5NnAG5jBAkABqIAtAMPdOxBR98OvBoVgCkYWQCAgSgA7UDxZw9ybQc6+gBMwQgSAAxER98OFH/2INd2oKMPwBQU+gBgIDr6dqD4swe5tgMdfQCmYAQJAAaio28Hij97kGs70NEHYAoKfQAwEB19O/BAxx4U+nagow/AFIwgAcBAFAV24IGOPSgA7UBHH4ApGFkAgIG8nV4KwMDGAx17UADagQc6AEzBCBIADEQBaAfybA9ybQce6AAwBYU+ABiI727bgeLPHuTaDuQZgCko9AHAQHx32w480LEHBaAdvB19pu4DcBsjSAAwEAWgHXigYwfv/SxxTwc6HugAMAUjCwAwEAWgHSgK7ODNs0SuAx2L8QEwBSNIADAQHX07UOjbgULfHizGB8AUFPoAYCA6+nag0LeDt/iT6PQGOjr6AEzBCBIADERH3w4U+nago28POvoATEGhDwAGoqNvBx7o2IFC3x509AGYghEkABiIAtAOPNCxA4W+PejoAzAFIwsAMIjH45FEAWgLpu7bwZvnoKAgBQUFuRwNWpr397ZERx+AORhBAoABLh7809G3A4W+Hchz4LncAxs6+gBMQaEPAAaio28HCkA7kGd70NEHYApGkABgIDr6diDPdqDQtwcdfQCmoNAHAAPR0bcDebaDt/ijyxv4eKgDwBSMLADAQAwW7UCe7UCe7cFDHQCmoNAHAAN5p3TT6Q1sFIB2IM/2INcATMEIEgAMxGDRDuTZDuTZHizGB8AUFPoAYCAWabMDBaAdyLM9WIwPgCko9AHAQCzSZgce6NiBQt8edPQBmIIRJAAYiALQDjzQsQOFvj3o6AMwBSMLADAQBaAdKADtQJfXHuQagCkYQQKAgejo24FC3w50ee1BrgGYgkIfAAxER98OFPp2IM/2oKMPwBSMIAHAQHT07UCe7UChbw86+gBMQaEPAAaio28H8mwHCn170NEHYApGFgBgIDq9dqAAtAN5tgcdfQCmoNAHAAPR6bUDBaAdyLM9yDUAUzCCBAAD0dG3A0WBHcizHWpqauTxeCQxdR+A+yj0AcBAdPTtQAFoB+90boq/wOa9nyXuaQDuYwQJAAaio28H8mwHHujYoW6hz0MdAG6j0AcAA9HRtwN5tgOFvh28Mzckcg3AfYwsAMBAdHrtQAFoB/JsBzr6AExCoQ8ABqLTawcKQDuQZzvQ0QdgEkaQAGAgCgM7kGc7kGc71O3o85AWgNv4LQQABmLqvh3Isx0o9O3g7egHBwcrKCjI5WgA2I5CHwAMxNR9O5BnO3jzzPe2Axt5BmASRhYAYCA6vXag02uHup1eBC7yDMAkFPoAYCA6vXag0LcDebYDHX0AJmEECQAGoqNvBwpAO5BnO9DRB2ASCn0AMBAdfTtQANqBPNuBPAMwid/OLdq5c6fWrFmjjIwMHT9+XMXFxQoPD1dkZKSGDRumhIQEjRs3TpGRkW6HCgBXjI6+HcizHSgA7cDUfQAm8bvfRKdOndITTzyhlStX1vvZiRMndOLECX399ddKTU1VcnKyUlJSXIgSAK4OHX07kGc7UOjbgan7AEziV4V+fn6+4uPjdejQIeezQYMGaciQIQoPD1d5ebny8vK0a9culZeXuxgpAFwdOr12oAC0A3m2Ax19ACbxm99EJSUlSkhIcIr8hIQEvfDCCxo6dGi9c6uqqpSZmanS0tK2DhMAWgSdXjtQANqBAtAOdPQBmMRv/uIsWrRIBw8elCQ98MADWrly5WV/kXbs2FHjx49vy/AAoEXR0bcDhb4dKADtwAMdACbxi1bRrl279Nprr0mSoqKi9Oqrr/LHEkBAo6NvBx7o2IEHOnbggQ4Ak/jFCPKVV15x9pOTk9WtWzcXowGA1kcBaAce6NiBQt8OdPQBmMT4kUV1dbVWrVrlHE+fPt3FaACgbVAA2oEC0A7k2Q509AGYxPgR5J49e3T27FlJUmhoqKKjo3XhwgWlpqZqzJgxioyMVEhIiK6//npNmDBBL7/8siorK12OGgCuDh19O1AA2oE824GOPgCTGP+baMeOHc5+VFSUjh07phkzZmj79u0+5xUUFKigoEDp6el6+umntXbtWsXGxrZ1uABw1Twej1Po09EPbBSAdiDPdqCjD8Akxhf6R48e9TmeMGGCcnNzJUkxMTGKjY1VcHCwcnJytHPnTklSfn6+4uPj9emnn+r2229v8N+vrKz0mQFQUlIiSSouLtb58+db8j+lxZ0/f17l5eUqKipShw4d3A4HrYQ828Hj8UiSTp8+rcLCQufzkpISBo0BpKysTJJUUVGhoqIiVVVVOZ8XFRW5GRpaSEVFhaTa8YX397Y3797P4P9Onz7t7BcWFqpTp07OZzU1NeQ5gDAOs4M/5dn7Cnnv2PGyPIZ77LHHPJJ8ts6dO3vWrFlT79zMzExPRESEc150dLSnsrKywX9/6dKl9f59NjY2NjY2NjY2NjY2NjZTt6NHjzZY5wZ5PI09CnDX/Pnz9c9//tPns3Xr1um+++675Pmff/657rrrLmfa6+uvv645c+Zc9t+/uKNfU1Oj4uJihYeHKygoqAX+C1rP2bNnFRUVpaNHj6p79+5uh4NWQp7tQJ7tQJ7tQJ7tQJ7tQJ7t4E959ng8Ki0tVZ8+fRr8imezp+6npKQoJSWluZdf0pNPPqkZM2b4fHbNNdf4HMfFxV22yPf+PDExUWvXrpUkvfPOOw0W+iEhIQoJCfH5rEePHlcYubu6d+9u/P8hcfXIsx3Isx3Isx3Isx3Isx3Isx38Jc+hoaGNntPsQr+wsFD79+9v7uWXdObMmXqfde3a1ee4oSK/7jneQj8rK6tFYgMAAAAAwB8Yv5xzeHi4z/Ett9zS6DWDBw929ktLS50FCwAAAAAACHTNLvSXLVsmj8fTotv8+fPr/e/ExMT4HF/c4b+Ubt26+RwHaqEfEhKipUuX1vvqAQILebYDebYDebYDebYDebYDebZDIObZ+MX4Dh06pAEDBjjHH3zwgSZNmtTgNV9//bVuu+025/jcuXPq0qVLq8UIAAAAAIApjJ+6379/f/Xv39853rt3b6PX7Nu3z9nv2bMnRT4AAAAAwBrGF/qSlJiY6Oy/9957jZ5f95xRo0a1QkQAAAAAAJjJLwr9BQsWqEOHDpJqV9HfsGHDZc/dvn271q1b5xzPnj27tcMDAAAAAMAYflHoR0dHa+HChc5xUlKSTzHvtWXLFk2ePFnV1dWSpBEjRmjq1KltFicAAAAAAG4zfjE+r8rKSo0bN05bt251Phs8eLBiY2MVHBysnJwcffXVV87PevfurezsbEVFRbkRLgAAAAAArvCbQl+SSkpKtGDBAq1atarB8+644w69++671hf51dXVys3N1Y4dO/Tll19qx44dysnJ0fnz5yVJo0eP1ubNm90NEvVUVVXpnXfe0apVq5Sbm6uTJ08qLCxM/fv3V2JiombPnq2IiAi3w8RV4N60w+HDh/Xxxx9ry5Yt2r17t/Lz83Xu3Dl169ZNffv2VVxcnJKSkjR69Gi3Q0UzFRYWatu2bdq+fbt2796tvLw8FRQU6Ny5c+rQoYPCwsJ06623Kj4+Xj//+c91/fXXux0yWsETTzyhv/71r85xv379dPjwYfcCQrOkpaVpzpw5V3TNvHnz9Nprr7VSRGgrO3fu1Jo1a5SRkaHjx4+ruLhY4eHhioyM1LBhw5SQkKBx48YpMjLS7VCvjMcPbdmyxTNv3jzPoEGDPF27dvV06tTJc+ONN3pmzpzpWbdunaempsbtEF3373//29O5c2ePpMtuo0ePdjtMXGTfvn2eYcOGNZi36667zrNx40a3Q0UzcW8Gvp07d3qGDx/eYI7rbvHx8Z4jR464HTaaYdKkSU3Oc0hIiGfZsmWe6upqt8NGC8rOzva0a9fOJ9f9+vVzOyw0Q2pqapPvZ+82b948t8PGVTh58qTnwQcfbFKuk5OT3Q73irW/6icFLhg1ahSr6TfizJkzKi8vdzsMXIFjx45pzJgxKigokCQFBQVp1KhRio6O1g8//KCMjAz9+OOPOnXqlKZNm6b09HTdc889LkeNK8W9Gfj279+v7du3+3w2cOBA3XrrrYqIiNCZM2eUlZWlY8eOSZI2b96suLg4bd26VQMGDHAjZLSAiIgIDR48WP369VPXrl1VXl6uAwcOaPv27bpw4YIqKyu1bNkyHTx4UG+88Ybb4aIFnD9/XvPnz1dNTY3boaCFxcTEaMyYMY2eN3LkyDaIBq0hPz9f8fHxOnTokPPZoEGDNGTIEIWHh6u8vFx5eXnatWuX347b/LLQR9P16tVLsbGxzvaf//xHf/vb39wOC5eQlJTkFPn9+vXT+vXr9ZOf/MT5eWFhoWbOnKlPPvlE58+f1/3336+8vDz16NHDpYhxNbg3A99NN92k+fPn66GHHqo3ZbumpkZpaWn69a9/rfLychUUFOjBBx9UVlaWgoKCXIoYVyo+Pl5TpkzRmDFjdNNNN13ynJMnT+rxxx93vna4YsUKTZkyRTNmzGjLUNEKnnnmGe3evVtS7d/wt99+2+WI0FLuuOMOpaSkuB0GWklJSYkSEhKcIj8hIUEvvPCChg4dWu/cqqoqZWZmqrS0tK3DvGoU+gFq/PjxOnLkiG644Qafz7Ozs12KCA3ZtGmTs9Bkx44d9f7772vIkCE+50RERGj9+vUaOnSoDh48qOLiYj377LP685//7EbIaCbuzcDXu3dvpaam6uGHH1ZwcPAlz2nXrp3mzp2rsLAwJSYmSpK++OILffTRR7r33nvbMlxchUWLFjV6Tq9evbRy5UqdPHlSmZmZkqTly5dT6Pu5b7/9Vk8++aQk6cEHH9TYsWMp9AE/sWjRIh08eFCS9MADD2jlypWX/XvdsWNHjR8/vi3DazF+8Xo9XLnIyMh6hQTM9Y9//MPZf+SRR+oV+V5dunTRn/70J+d4+fLlunDhQqvHh5bDvRn4Ro8erdmzZ1920FDXfffdp+HDhzvHGzdubM3Q4JKgoCCfRb6+/vprF6PB1fJ4PJo/f74qKysVFhamv/zlL26HBKCJdu3a5SygGBUVpVdffbVJf6/9EYU+4LJz587pk08+cY4bW/F1+vTp6tq1qySpuLhYn376aavGB6B13Xnnnc4+K3UHrmuvvdbZ98cpoPj/Xn75ZW3btk2S9H//93+67rrrXI4IQFO98sorzn5ycrK6devmYjSti0IfcFlWVpYqKysl1XbsY2NjGzz/mmuuUVxcnHPsnQoKwD/V/U5+dXW1i5GgNe3du9fZv/HGG90LBFfl6NGj+v3vfy9JuvvuuzV37lyXIwLQVNXV1T6vaZ8+fbqL0bQ+vqMPuGzfvn3O/pAhQ9S+feO35W233aaPP/643vUA/I93MS+pdhohAk9BQYGee+4555jv5/uvhQsXqrS0VB07dtTy5ctZPDNAnTlzRu+++65yc3NVUlKi7t27q0+fPoqLi9OQIUPIu5/as2ePzp49K0kKDQ1VdHS0Lly4oDfffFNvvfWWcnNzdfr0aUVERGjo0KGaOnWq5s6dq5CQEJcjbx4KfcBl+/fvd/b79evXpGvqfsf722+/bfGYALSN/Px8n1k5Y8eOdTEatKTy8nIdPnxYH374oZ599lmdOnVKkjR48GCnIwz/snr1an3wwQeSpMWLF2vw4MEuR4TWsn79eq1fv/6SP7v55pu1ePFizZ07l4Lfz+zYscPZj4qK0rFjxzRjxox6r8QtKChQQUGB0tPT9fTTT2vt2rWNzrg1EYU+4LKioiJnv1evXk26JjIy0tkvLi5u8ZgAtI0nnnjCma5/ww03aMqUKS5HhOb67LPPdPfddzd4zsSJE7Vy5cqA/k5ooCoqKtJjjz0mSRo4cKCWLFnickRwy3fffaf58+frvffe0+rVq9WlSxe3Q0ITHT161Od4woQJys3NlSTFxMQoNjZWwcHBysnJ0c6dOyXVPpCPj4/Xp59+qttvv73NY74afEcfcNm5c+ec/U6dOjXpmrrn1b0egP9444039K9//cs5fuqpp/x2eiAaFhYWplWrVmnjxo3q0aOH2+GgGR5//HH98MMPkmoX8+JeDUw33HCDfvvb32rTpk06evSoKioqVFZWpv379+ull15STEyMc+4HH3ygpKQk1dTUuBgxrsSZM2ec/T179ig3N1edO3fWmjVrtG/fPq1YsUKpqan66quvlJmZqYiICEm1M7QeeOABVVVVuRR581DoAy6rqKhw9jt27Nika+oOMH788ccWjwlA6/ryyy/16KOPOsezZs1SUlKSixHhavXp00fJyclKTk7WwoUL9fDDD2v48OFq3769Tp8+rVmzZumee+7R//73P7dDxRX66KOP9Oabb0qqfQVuQkKCyxGhNUybNk2HDh3Sc889pwkTJqhv374KCQlR586dNXDgQC1YsEDffPONz9uRNmzYoLffftvFqHElysrK6n321ltv6f7776/3eUJCgjZs2KB27WrL5by8PK1cubLVY2xJTN1vQykpKUpJSWnRf/PJJ59kUR8/d8011zj7TX1S6F2lX2r6LAAAZjh06JCmTJniPOQbOnSoz+t+4J8GDBhwyb/xBQUFWrJkidLS0vTf//5XI0aM0ObNmzV06FAXosSVKisr0y9/+UtJUnh4uM+iiggsTZlt07FjR7322ms6cOCAtm7dKkl65pln9NBDD7VydGgJdcfckhQXF6f77rvvsufHxcUpMTFRa9eulSS98847jb4G2yQU+m2osLDQZ+G1llB3Cgr8U9euXZ39pnbn655X93oAZjtx4oTGjRun77//XlJtcZienq7u3bu7HBlaS58+fZSamqru3bvrxRdf1OnTpzVz5kzt3r1bwcHBboeHRixZskSHDx+WJD3//PPOVF7Yq127dlq6dKmzeOqePXt07Ngx9e3b1+XI0JiLx8wNFfl1z/EW+llZWa0SV2th6j7gsvDwcGf/5MmTTbrGWyRIUs+ePVs8JgAtr6ioSOPGjVNeXp4kqXfv3srIyFDv3r1djgxt4amnnnIe6Ozbt08ffvihyxGhMTt37tTf//53SbXTeB955BGXI4IpRo0apQ4dOjjHvOrYP9Qdc0vSLbfc0ug1dd+uUVpaqtLS0haPq7XQ0W9Dy5Yt07Jly9wOA4YZNGiQs3/kyJEmXZOfn+/s110YBoCZzp49q3vvvddZ3TciIkIZGRnq37+/y5GhrXTu3FkjR45Uenq6JGnbtm2aPHmyy1GhITk5Oc5Ca/n5+RoxYsRlz/Uu1CfVztype+4f/vAHTZo0qfUCRZvr0KGDIiIidOLECUm1s3ZhvovHzE2ZFXvxW1JKS0v95s0pFPqAy+o+Kdy9e7cuXLig9u0bvjW9r/y4+HoA5ikrK9PEiRP11VdfSZJCQ0OVnp7epE4CAktYWJizX/fVqjBfXl6eMxunMVVVVcrOznaO6z4EQOCou7Abr9jzD7feeqvPcVPeXHVxBz80NLRFY2pNTN0HXDZy5EhnFf2ysjJ9+eWXDZ5fWVmpL774wjm+5557WjU+AM1XUVGhqVOnatu2bZJqu7obN270u3fxomV4u38SX7sC/NnBgwd19uxZ57hPnz4uRoOm6t+/v89Mur179zZ6Td2vZfTs2dOvHupQ6AMu69q1q8aMGeMcp6WlNXj+unXrnKeLPXv21KhRo1ozPADNdP78eU2fPl2ZmZmSal+LuX79et15550uRwY3FBUV6fPPP3eOmY1lvtmzZ8vj8TRpS01Nda7r16+fz89mz57t3n8EWsXrr7/u7IeGhmrYsGHuBYMrkpiY6Oy/9957jZ5f9xx/G3NT6AMGWLhwobOflpbmfI/3YuXl5frjH//oHP/iF79odJo/gLZXXV2tpKQkbdq0SZLUvn17rVmzxlmlGf6vuLi4yefW1NToV7/6lfNq1JCQEL6fDxikKVO4vbKysvT88887xzNnzmQs5kcWLFjgLKSYlZWlDRs2XPbc7du3a926dc6xvz20o9AHDDBp0iTdfffdkmqn5k+ePFk5OTk+5xQVFWnatGk6cOCApNpu/uLFi9s8VgAN83g8mjdvnvM6nnbt2unNN9/U1KlTXY4MLWnFihWKjY3VihUrfKbwXiwnJ0cTJ07U6tWrnc9+97vf1Vv9GYB71q5dq+HDh2vFihUqKSm55DkVFRV68cUXNXbsWFVUVEiSevTooaVLl7ZlqLhK0dHRPg22pKQkn2Lea8uWLZo8ebKqq6slSSNGjPC7v+NBHo/H43YQaB0TJ05UQUGBz2fff/+98wq3Ll266Kabbqp33aZNm/iukQuOHTum4cOHO9/hDAoK0ujRoxUdHa0ffvhBGRkZKi8vl1TbHUxPT/eZ8g//wb0Z2F566SUlJyc7xzfffLN++tOfNvn6lJSU1ggLLeyFF17Q448/Lqn2d3JMTIwGDRqksLAwBQUFqaioSDk5Oc7DWa/p06dr9erVdAADTFpamubMmSOpdur+4cOH3Q0IV6Ru/rz3c0xMjMLCwlRdXa3jx4/r888/93mo16lTJ6Wnp/vddG7UNtXGjRunrVu3Op8NHjxYsbGxCg4OVk5OjrOArlT7Otzs7GxFRUW5EW6z8VcmgO3du7fB17WVlZXpm2++qfd5VVVVa4aFy+jbt68yMzM1a9Ys7dq1Sx6PR5s3b9bmzZt9zrv22muVmppKke/HuDcD26lTp3yOv/vuO3333XdNvp5C3z94F1GVpAsXLmjPnj3as2fPZc/v1q2bli1bpt/85jcKDg5uixABNENT7ufhw4crLS2NtTb8VEhIiN5//30tWLBAq1atklS76F7dhfe87rjjDr377rt+V+RLFPqAUWJiYpSdna3Vq1dr1apVys3N1cmTJ9WjRw8NGDBAiYmJmjNnjiIiItwOFQCstmDBAo0ZM0YZGRnKzs5Wbm6u8vPzdebMGUlS9+7d1bt3bw0bNkxjx47V9OnTm/TOZgBtb9asWRo4cKCysrL0xRdfKC8vT4WFhSoqKlJNTY1CQ0PVv39/jRgxQjNmzNBdd93ldsi4SqGhoXr77bf16KOPasWKFfrss890/PhxVVdXq1evXhoxYoR+9rOfadq0aQoKCnI73GZh6j4AAAAAAAGExfgAAAAAAAggFPoAAAAAAAQQCn0AAAAAAAIIhT4AAAAAAAGEQh8AAAAAgABCoQ8AAAAAQACh0AcAAAAAIIBQ6AMAAAAAEEAo9AEAAAAACCAU+gAAAAAABBAKfQAAAAAAAgiFPgAAAAAAAYRCHwAAAACAAEKhDwAAAABAAKHQBwAAAAAggPw/P1N2yMavCOIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1000x800 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import jax.numpy as jnp\n",
    "\n",
    "lm = jnp.array([1, 2, 3, 4])\n",
    "M = len(lm)\n",
    "D = jnp.array(lm)\n",
    "a = jnp.min(lm)\n",
    "b = jnp.max(lm)\n",
    "t = jnp.linspace(-1, 6, 1000)\n",
    "u = 0.5 * jnp.ones(M)\n",
    "rho = 1\n",
    "def fun(lam):\n",
    "    return 1 + rho * jnp.sum(u**2/(D - lam))\n",
    "res = [fun(lam) for lam in t]\n",
    "plt.figure(figsize=(10,8))\n",
    "plt.plot(t, res, 'k')\n",
    "#plt.plot(jnp.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",
    "Lets explain a little bit..."
   ]
  },
  {
   "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 $N$ steps we have the factor\n",
    "\n",
    "$$\\left(1 - \\frac{1}{N}\\right)^{\\frac{N}{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 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  },
  "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": {},
  "rise": {
   "scroll": true,
   "theme": "sky",
   "transition": "zoom"
  },
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}