{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Homework 8\n", "\n", "## Distributional Approximation and Gaussian Processes\n", "
\n", "This notebook is arranged in cells. Texts are usually written in the markdown cells, and here you can use html tags (make it bold, italic, colored, etc). You can double click on this cell to see the formatting.
\n", "
\n", "The ellipsis (...) are provided where you are expected to write your solution but feel free to change the template (not over much) in case this style is not to your taste.
\n", "
\n", "Hit \"Shift-Enter\" on a code cell to evaluate it. Double click a Markdown cell to edit.
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***\n", "### Link Okpy" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from client.api.notebook import Notebook\n", "ok = Notebook('hw8.ok')\n", "_ = ok.auth(inline = True)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "from scipy.integrate import quad\n", "#For plotting\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Problem - Back to Quasar\n", "\n", "In HW5, we performed Principal Component Analysis (PCA) on the quasar (QSO) spectra from the Sloan Digital Sky Survey (SDSS); we filtered for high $S/N$ to apply the standard PCA and selected 18 high-$S/N$ spectra of QSOs with redshift 2.0 < z < 2.1, trimmed to $1340 < \\lambda < 1620\\ \\mathring{A}$. Then, using the first three principal eigenvectors from the covariance matrix, we reconstructed each of the 18 QSO spectra.\n", "

\n", "In this assignment, we do Expectation Maximization PCA with and without per-observation weights. We use a simple noise fit of PCA components to individual spectra. Finally, using a Gaussian process, we compute the posterior distribution of the QSO's true emission spectrum and sample from it. \n", "

\n", "The following analysis is based on https://arxiv.org/pdf/1208.4122.pdf, and https://arxiv.org/pdf/1605.04460.pdf\n", "

" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Load data\n", "wavelength = np.loadtxt(\"HW5_Problem2_wavelength.txt\")\n", "X = np.loadtxt(\"HW5_Problem2_QSOspectra.txt\")\n", "ivar = np.loadtxt(\"HW5_Problem2_ivar_flux.txt\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Data dimension\n", "print( np.shape(wavelength) )\n", "print( np.shape(X) )\n", "print( np.shape(ivar) )" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "In the above cell, we load the following data: wavelength in Angstroms (\"wavelength\"), a 2D array of spectra x fluxes (\"$X$\"), and another 2D array of inverse variances ($1/\\sigma^2$) of the flux array (\"ivar\").\n", "

\n", "We have 824 wavelength bins, so \"$X$\" is a 18 $\\times$ 824 matrix, each row containing fluxes of different QSO spectra and each column containing fluxes in different wavelength bins. (e.g. X[i,j] is the measured flux of QSO $i$ in wavelength bin $j$.) Similarly, \"ivar\" is a 18 $\\times$ 824 matrix. (e.g. ivar[i,j] is the inverse variance of the flux of QSO $i$ in wavelength bin $j$.)\n", "

\n", "Remember that in HW5, we computed the eigenvectors of the covariance of the quasars, sorted by their descending eigenvalues; we call them the principal components (henceforth denoted by $\\phi$). Suppose that we have $k$ eigenvectors, each of length 824. Construct the matrix of eigenvectors $\\phi = [\\phi_1\\ \\phi_2\\ ...\\ \\phi_k]$, with $\\phi_i$ the $i$th principal eigenvector.

\n", "We can reconstruct the data as:

\n", "$$ \\hat{X} = \\mu + \\sum_k c_k \\phi_k $$\n", "
\n", "where $\\mu$ is the mean of the initial dataset and $c_k$ is the reconstruction coefficient for eigenvector $\\phi_k$.\n", "

\n", "More specifically, we define $\\mu$ as:\n", "

\n", "$ \\mu$ $ =\n", " \\begin{bmatrix}\n", " \\overline{x}_1 & \\overline{x}_2 & \\dots & \\overline{x}_{824} \\\\\n", " \\end{bmatrix}.$\n", "

\n", "The mean-centered data matrix $X_c$ can be defined as:\n", "

\n", "$X_c = X - \\mu$ $ =\n", " \\begin{bmatrix}\n", " x_{(1,1)} - \\overline{x}_1 & x_{(1,2)} - \\overline{x}_2 & \\dots & x_{(1,824)} - \\overline{x}_{824} \\\\\n", " x_{(2,1)} - \\overline{x}_1 & x_{(2,2)} - \\overline{x}_2 & \\dots & x_{(2,824)} - \\overline{x}_{824} \\\\\n", " \\vdots & \\vdots & \\vdots & \\vdots \\\\\n", " x_{(18,1)} - \\overline{x}_1 & x_{(18,2)} - \\overline{x}_2 & \\dots & x_{(18,824)} - \\overline{x}_{824}\n", " \\end{bmatrix}$\n", "

\n", "where $x_{m,n}$ denote the flux of $m$th QSO in $n$th wavelength bin, and $\\overline{x}_k$ is the mean flux in $k$th wavelength bin.\n", "

\n", "$ \\mu$ $ =\n", " \\begin{bmatrix}\n", " \\overline{x}_1 & \\overline{x}_2 & \\dots & \\overline{x}_{824} \\\\\n", " \\end{bmatrix}$\n", "

\n", " 1. Plot $\\mu$ as a function of wavelength $\\lambda$.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "\"Expectation Maximization (EM) is an iterative technique for solving parameters to maximize a likelihood function for models with unknown hidden (or latent) variables. Each iteration involves two steps: finding the expectation value of the hidden variables given the current model (E-step), and then modifying the model parameters to maximize the fit likelihood given the estimates of the hidden variables (M-step).\" (https://arxiv.org/pdf/1208.4122.pdf)\n", "

\n", "Now, do Expectation Maximization PCA. In this case, we wish to solve for the eigenvectors, and the latent variables are the coefficients $c$. The likelihood is \"the ability of the eigenvectors to describe the data.\"\n", "

\n", "First, find the eigenvector $\\phi_1$ with the highest eigenvalue (the first principal eigenvector):\n", "

\n", "1. Initialize: Let $\\phi$ is a random vector of length 824.\n", "

\n", "2. E-step: For each QSO $j$, $$c_j = X_{row\\ j} \\cdot \\phi$$
Here, \"$\\cdot$\" represents a dot product, so $X_{row\\ j}$ and $\\phi$ are vectors of length 824, so $c_j$ is a number. $c = [c_1\\ c_2\\ ...\\ c_{18}]$ is a vector of length 18 (because we have 18 QSOs in this problem). So for each QSO $j$, we solve the coefficient $c_j$ which best fits that QSO using $\\phi$.\n", "

\n", "3. M-step: $$\\phi = \\frac{\\sum_j c_j\\ X_{row\\ j}}{\\sum_j c_j^2} $$\n", "
\n", "Using the coefficients $c_j$, we update $\\phi$ to find the vector which best fits the data given $c_j$.\n", "

\n", "4. Normalize: \n", "$$ \\phi = \\frac{\\phi}{|\\phi|} $$\n", "

\n", "5. Iterate until converged. Once converged, $c_1 = c$, and $\\phi_1 = \\phi$\n", "

\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After you get $\\phi_1$, subtract the projection of $\\phi$ from $X$ ($X - c_1 \\otimes \\phi_1$, where \"$\\otimes$\" is the outer product (https://en.wikipedia.org/wiki/Outer_product). $c_1$ is a vector of length 18, and $\\phi_1$ is a vector of length 824, so $c_1 \\otimes \\phi_1$ is a $18 \\times 824$ matrix.) and repeat the EM algorithm.\n", "

\n", "(So to find $\\phi_2$, you should use a data matrix $X - c \\otimes \\phi_1$. To find $\\phi_2$, use $X - c_1 \\otimes \\phi_1 - c_2 \\otimes \\phi_2$), and so on.\n", "

\n", " 2. Using EM PCA, find the first three principal eigenvectors $\\phi_1, \\phi_2, \\phi_3$ and plot them as a function of wavelength.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, reconstruct the data using the first principal eigenvectors:\n", "

\n", "$$ \\hat{X} = \\sum_{k = 1}^3 c_k \\otimes \\phi_k $$\n", "
\n", " 3. For any one QSO spectra, plot the original and reconstructed spectra, using the above equation.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Alternatively, you can also reconstruct the data using \"PC scores.\" (Call the PC score matrix $Z$) \n", "

\n", "$$ Z = X_c \\phi $$\n", "

\n", "Then, we can reconstruct the data by mapping it back to 824 dimensions with $\\bf \\phi^T$:\n", "

\n", "$$ \\hat{X} = \\mu + Z \\phi^T $$\n", "
\n", " 4. For any one QSO spectra, plot the original and reconstructed spectra, using PC scores.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, include noisier QSO spectra." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Load data\n", "wavelength = np.loadtxt(\"HW5_Problem2_wavelength_300.txt\")\n", "X = np.loadtxt(\"HW5_Problem2_QSOspectra_300.txt\")\n", "ivar = np.loadtxt(\"HW5_Problem2_ivar_flux_300.txt\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "ivar[ivar==0] = 1.e-4" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "(824,)\n", "(2562, 824)\n", "(2562, 824)\n" ] } ], "source": [ "# Data dimension\n", "print( np.shape(wavelength) )\n", "print( np.shape(X) )\n", "print( np.shape(ivar) )" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "We now have 2562 quasars (including 18 high $S/N$ quasars we had before). The below cell plots the spectra of two quasars; you can see how noisy they are." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3sAAAEKCAYAAACrLEyuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXm4HEW5/7/VPXOWLCSEhCCBEHZBEJGACyrI7sqV6wZe\ndy+i/vR6RRGvAnrFi8iiAqIgIKIoyiKggYQ9CVnJSvZ9IevJyUly9jkz0/X7o6e6q6urunvmzJyZ\nOef9PE+enOm1pqe66n3r3RjnHARBEARBEARBEMTgwqp2AwiCIAiCIAiCIIjyQ8oeQRAEQRAEQRDE\nIISUPYIgCIIgCIIgiEEIKXsEQRAEQRAEQRCDEFL2CIIgCIIgCIIgBiGk7BEEQRAEQRAEQQxCSNkj\nCIIgCIIgCIIYhJCyRxAEQRAEQRAEMQghZY8gCIIgCIIgCGIQkqp2A4pl7NixfNKkSdVuBkEQBDEA\nLFy4sJVzPq7a7agXaI4kCIIYGiSdH+tO2Zs0aRIWLFhQ7WYQBEEQAwBjbEu121BP0BxJEAQxNEg6\nP5IbJ0EQBEEQBEEQxCCElD2CIAiCIAiCIIhBCCl7BEEQBEEQBEEQgxBS9giCIAiCIAiCIAYhpOwR\nBEEQBEEQBEEMQiqm7DHGHmCMtTDGlsccdyZjLMcY+3il2kIQBEEQBEEQBDHUqKRl70EAl0QdwBiz\nAdwM4LkKtoMgCIIgCIIgCGLIUTFlj3M+A0BbzGHfBPA4gJZKtYNwmb+pDet2d1S7GQRBEARBEESN\n8NSS7TjQk612M4gKUrWYPcbYBAAfA/DbarVhKPHJe+bgwl/OqHYzCIIgCIIgiBpg3e4O/NcjS3DN\nY0ur3RSigqSqeO9fAfg+59xhjEUeyBi7EsCVADBx4sQBaBpBEARBEARBDD52t/fC4RxdfXkAwM4D\nvVVuEVFJqqnsTQbwSEHRGwvgg4yxHOf8SfVAzvm9AO4FgMmTJ/MBbSVBEARBEARBDBLe8X8vAgCe\n+Pq7AQBxRheivqmassc5P1r8zRh7EMC/dIoeQRAEQRAEQRDlxXFc+4lFut6gpmLKHmPsrwDOBTCW\nMbYNwA0A0gDAOf9dpe5LEARBEARBEEQ0mZwDALDIsjeoqZiyxzm/vIhjv1CpdhCl89jCbbjklMMw\norGa3r4EQRAEQRBEuenNujF7pOoNbqqWjZOoDL3ZPM695WX87bWt/brO4q378N1Hl+JH/1hWppYR\nBEEQBEEQtUJPQdkjy97ghpS9Qcabr5uKzXu78f3H+6ekdfTmAACtnX3laBZBEARBEARRQ/RmXTdO\n0vUGN6TsFVi7uwN/mrO52s2oGfKFoF2bonYJgiAIgiAGHWTZGxpQMFaBi381A5wDn33XpGo3pSbI\nFZS9FCl7BEEQBEEQg46MUPbI9DOooZ+3AKfqfQHyjmvaJ8seQRAEQRDE4KOnjyx7QwFS9oYwi7bu\n8zIxqeRdXY+UPYIgCIIgiEHCQ3M2e3/35lwZUITuEIMTUvaGKDv29+Cyu2fjB0/oE7nkFMve/E1t\n2HmgZ8DaRxAEQRAEQZSP3e29uP6pFd7n7ftcuU7U2yMGJ6TsDVE6M262zeXbD2j355WYvU/eMwcX\n3DZ9YBpHEARBEARBlBXVgrd8RzsAoI+UvUHNkFb2PnXPHPxj8bZqN8Nj/qY2PLVke7WbAcBP0GJL\nUbtdfXqXzyieXroDF94+HZyCIgmCIAiCIKqGGpq3ZW8XAKCrL1eF1hADxZDOxjlvUxvmbWrDx04/\nwtvGOQerUqDqJ++ZAwC49G0TBuyeJhWsvScLoP/ZOL/9yGI43F1NStkU/0cQBEEQBFENVMteNs/R\nmLKwcU8XNuzpxLHjRlSpZUQlGbKWPccQjEoGKJcbp6wCANj9VNC48n+98o/F27C3M1PtZhAEQRAE\nQZSELhHLOSeMAwC8uq4VD8/bgmyeXDoHG0NW2cs6+s4sXoOf/HMFrnxowcA1qEaxGSuLC6ZTx1r0\nzgM9+O+/LcVVf15Y7aYMGp5ZttOYCZYgCIIgiPKTzYdlsVHNaQDADU+vwA//sRzH//DZgW4WUWGG\nrrKn6fAAPMXmD7M247mVuweySQOCTnFr783ipmdXaQN0bYuhPxl5xe3qWNdDNuc2fld7b5VbMjiY\nvb4VX394EW6ZtqbaTSEIgiia3e29uHnqaqOHEEHUKjmNoUMoe8TgZegqe4bMQ4Np6P70vXNC29bu\n7gxtu23aGtwzfSOeLCSHkRXClMXKUn9lMNRwqWeFtZbYU3CH3dNBbrEEQdQf3/n7Evz2lQ1YuHVf\ntZtCEEWRi7DsEYOXoavsFVY31AQkg0mgn7uxLbTt4l/NAADI37o36z4LoZDJipldJmWvnt04Rb6e\nOv4KNYWo59OQGrLDD0EQdYzwgqE5gag3chp5blhjyqupLCCr9eBiyEpbwo1TzRDJS7TtrW/pwLrd\nHf1uVzURTyLg4sqAfFli9vp9CWKQIJS9RlL2CIKoY6ikEFEP5PIObnpmFXYe6EFOk3wlZTGMaAwm\n519UsFpvau3CrPWtA9JOonIMWWlLuHGmrOAjMI3dnZmc1/l1XHD7DFz4yxlla99Aoiq4fdJgwDmQ\nN8Q3FnWPOp4Uq1SJY9CSKSRmaUzZVW4JQRAEQQxuXljVgntmbMSt09Zq81VYGmXv479zw4Def+sr\n+Mx98waknUTlGLrKXkGhUU3XJv7fXxbhsrtno703m/geeYdj7sa93mfdiko1uOL3c/Ha5rDiKpQa\nUWQTcL+DLqC3WAaDZa+eFdZagtw4CYIYDFSrJi9BFMPagtfZ2BEN2rAcnWWPGFwMWWmrT1L2ZCHe\nJM8v23YAAJDJJld8LvrldHz63rmYvaEVy7cfwHE/fBYvr2kpvdFlYvaGvfiffywLbV+4ZR8mXTsF\nH71rlrfN4bxMbpzha7ywcjfauvr6fe1KQxN6eemrkhvnW388Ddc8tnRA70kQxOCFFgCJemDL3m4A\nwLiRjdqyYzZjaG4gT5vBzJBV9oQp27ZYwKxtitkTAn/SmL7Vu9qxYY9rIWtpz2BxwQX0+Ros5yDm\nq78v2Bba5zg8tBK0fX9P0ZOcqux19GbxlYcW4IsPvlZcY6sITevlwYvZSw/M8NPRm8VX/vga2ntz\n2j5OEARRDAy0AEjUD8IjLe9wbTZO22KhZIXE4GLIKnvCpTJtMc+lEzBb9orNyLivy3f3ZMx3WdPV\nsqsGuve6wQ53hzwPKntb9nbh7J+/hN+8vL6o+6nPTSwubdoTLgVRa4hHRYu45SGTc2P2dP3NRE9f\nHtv395R0vyeX7MALq6pvUScIGcbYA4yxFsbY8ohjzmWMLWGMrWCMTR/I9hFmSk3kRhADRWtnBhsL\n8lVXJgfADV/K6yx7Fksc0kTUJ0NW2fPcOG1F2TMcbxWp7PUWBFrAtQoKZS9bhbg9nRVOJKZZ39KJ\nRxe61g5dDJXDg6UYhNtlsQXnQ37idVTOgHv/10Fj6wBh2Vu67QDeaOtOdM5n75+Hs3/+Ukn3oymM\nqFEeBHCJaSdjbDSAuwF8lHP+FgCfGKB2EQRR55zzi5dx3m3T8asX1mJft2t86MtzbYIWUvYGP0NW\n2fNKL1iWkn3S4MZZEBmT1ouTY/sYgAbb9YdWLXvrWzoqXs9E12RL88trlT2HB+qyDC8E8Xb05opq\ng/rcPEtpUVepDqJP1INiWg+Id+CfS3fgvb94OdE5C7aUXryYQi6JWoRzPgNAuBiqzxUAnuCcby0c\nT+bpGoHcOIlap6vPNTj86oV1WLWzHYBrbNAl3CNlb/AzdJW9nJ+gJRizp6dY5aQ3mw98FopURlL2\nXt+2HxfcPgP3v7op4VWDPDxvC069YVps/Jxur62RgLVunA4PKKPiVh1FZCWVzxP4rpG1r0HVQRPr\niqQLJpXixn+txM4DpbmEEsQAcgKAgxljrzDGFjLGPmc6kDF2JWNsAWNswZ49ewawiUMT8vIg6pG+\nnGOM2dMpe7WSQZ7oP0NX2cuLOnvMU/wA12J11Z8Who63RIKWhIJqR8a3fDHmF2+X3TjXt7j+1Mu2\nHyiy9S7XPbkcHZmcNpWujE64tjQv9q72Xs25CFj2xLXae/pn2ePe9qIuUxVE0+ugqfXBAD9IdRX+\nvlc34Tt/o6ycRM2TAnAGgA8BuBjAdYyxE3QHcs7v5ZxP5pxPHjdu3EC2kSCIOsG17GmUPca0BoBM\njeSYIPrPkFT2/r7gDcwp1L+zWDBm7y/ztmDqil3Gc5MaJWTLF4Nf3kF+eXoK1r+mfmYljGuSrs1J\nDfaOkqBF/N1X5IqPuEQml8eJP3oWTy/Z4batDlSoemhjPTHQlr0Fm8OecplcXnMkQdQU2wBM45x3\ncc5bAcwAcFqV20SA3DiJ+iSbd7TWOtvWW/ZqJaEg0X+GpLJ3zWOv46E5WwC4sWuy4qKuenDO8eE7\nZ3qZAJPKqXLMHgf3sk/2KVZEwC15sGZXR9HfQ5SDiBKef/7sau1+XZCujrxSekG3KpQE0YaW9gwy\nOQc3P7saQH24SHqWvTpoaz2QsOuVhR37e/DE4u1VbQNBlMhTAN7DGEsxxoYBeAeAVVVuEyFBwwhR\nT/TlDAlamF7ZI8ve4GFIKnsqcudXk6Vkcg6Wb2/392skfp1rZ0BByvuFyWVlTy4o/u+/nV10u5Mk\nDvnd9A3a7UktGw4PJmjRpe3Vsacjg5nr/NgR0VZxLeFGWg8KVHxE5+DjtufW4KklYSWpHCTtQ+Wg\nK6N3N650UiSCiIMx9lcAcwCcyBjbxhj7MmPsKsbYVQDAOV8FYCqA1wHMB3Af59xYpoEgCCIKt/SC\nLju7SdkjD5jBQqraDag2DEE3TtXCrQqLOhFRZ+3KSgJtNu+g0SnU2ZOtiJKS2R9zeZzCpFNQ42Td\n9x4/Ftv39cDhPHB+Uovgp+6Zg42tXd5n8bXFQCMGlnpwkaxkNs7NrV1Ys7sDF7/lsPJfvEQ457jz\nJbeO4qVvm1D26+sCxAeauDhXgqg0nPPLExxzC4BbBqA5RAmQMydRa6za2Y5hDbZ2X1/OCcimAsug\n7MlyKefc8yYj6o8hb9ljDIEELarVoSsTXNnQCYm62nn5vGwN8y17jibZCQA0asoexJHEjRMoTUn5\nzoUnIG1byDs8IJwnFZJlRU9uo0j7a9elZa/8vP+2V/BVTUKgaiJbnCuB2of+75mB90yrdkZQgiBq\nn7auPvxi6mrjvEejCFENOOfo7tN7rXzg1zNxzi2vaPf15f1snFe8Y6K33WTZ++hds7y/aX20vhny\nyh4QtLblFSGwU7XsJYx/k619WYd7L4p8pHytxrR+JSYJcYJrKYJtyrJgWSxUVL2/MXtioLETKqq1\ngGfZq8i1K3DREpm6fBeO/+EzWLWz+PjRYlD70L0zNiY+t5hSHe29WfRm9RZzsuwRBBHHDU+vwN2v\nbMCMtVTOgqgdHnntDZx8/TRs3dtd1HkvrW7B7c+vRcpiOGbscG+7ZemzcfZIJcTqQVYjzJCyB2iz\nTQq6lNUTnYyos+zJhStzecez6AlhdfHWfV5GUKA0y15Um2RKeUVti8FiriVSVoBLrbsiLpELuXHW\nPn6ClnpobencMm01snnuFWBt6EefjEJX1DUpxfwEb/3xc/jIXa9q96mLOgRBECo9hfnftDhETm1E\nNXh2uZsxfkNrZ0nnj2hK4UtnH+19thnzyoOZoAXS+mbIK3sMQWVJjSdSY/aSunHK18nl/YyWYuvH\n7p6Ntbv9F1Uuv1CsUhFbVL0E2TpVSMWb5zzg2tpfy16+jt045aZuau3Cw/O2VKM5FUP8tGKBr1IF\nVXUxe9NW7MKka6dgX4wLabm6CyVoIYihTVtXH56LKLOUBBpFiHpkeEMKlsVw+KgmAGJxP1rZqwdZ\njTAz5JS9kGLEWLRlT4nZ05YxyMW5cTqeJcH0wjSmfDfOpC+VeDXjLXvFv6W2xcCY68ZZSsyeijjN\nc+PU+IfXKrrSCx+961X88B/LB5W1T+3bldKHdH1IZI1dvyd6pbJcz5ssewQxtPnyH1/DlX9aiAM9\nWeMxxmGifqYvYhBTajcc2eTmZhR5H5LIY+TGWd9UTNljjD3AGGthjGlTRTPGPsMYe50xtowxNpsx\nNiDFYnWWKVmADMfsBScCXYfXCY45JeumH/elf2Fkl7liX6n4mL0iL4hCwG7BjVO+fv8te+7/daTr\naX8zUSNxMBmIBmow12avLbwvcd1CnPn4wm244anSs9APYPUHgiBqkE2FJGJJFjBLnfcIohKUsugp\ny1zDG4Wy535OWSxW7iRlr76ppGXvQQCXROzfBOAczvmpAH4K4N4KtsVDN7DLm9T9fYrLmbxbCKi6\numHy5JDLO17pAdP7It936bb9mL2hVX+ghvhsnKVZ9myLlVxnT0W4zWWVmL16IOrxDZQf+81TV+OZ\nZTsreo+BUoB0z0ykeI5L7Sz6+tWPLsUf55TuRkvxBwQxtBHjepKZSK43dqA7i417uiKOJoiBoZhS\nCA9+8SxccNKhAIARirJnWSxWTqQF0vqmYsoe53wGgLaI/bM55/sKH+cCOKJSbZFRV+iWvrEfbV0Z\n77MqBKpxS0LYfH3bfhz/w2fxypqWUG0+9zyOCaObAbgKTpwbp3j5AOCyu2fjit/PS/aFIq4JAMMa\n7BItexZYwcVVfiZJ6+ypiEuIMhe2VT8exI732+kWCgZGafjtKxvw9YcXVfQeXLG+Vgq9ZS+ZxTfJ\n416+/YCxmLqA3DgJYmhTzCJoT5+v7H34rplo7cxEHE0QA0MxS+Zp20Jzgytnespe4Qo2Y7HKXM5x\n0NLRW0oziRqgViTuLwN4diBulNcoK7+Yusb7WxVEVcFXWKgWbHb11FfW7NEKxzmHY1RzGs1pO5CN\nU+XQkY0AgGPGDdfuT0KUwuGa50u07DHXshcV0wi4k+bRP5iC259fG9tGYQ21a6XnJcCL2dPsq2cL\n0aKt+zDp2il4fdt+AH6a5Uq7LOmsw7Jlb86GvSUnUFm0dR8+fOerXgygCUrQQhBDGzECRM2fYk+v\nlIL+jbYefz8NI0QVKKXfpW3mZX0Xyp5YXOWIz+1w18vrcdbPXsSO/T2RxxG1SdVFbsbY++Eqe9+P\nOOZKxtgCxtiCPXv6V+8mG+NyGbLsqcpe4aNIU6tavvzzHKRsN51tNs+NE0pj2sK4kY39iv2KKg3A\neWkDQ8pz41RcUjUN3bavB5wDd7y4DgDw1JLtoWPE9xc1DXU1XWoezXOsp1gOzjlue24NthcG6xdX\n7QYATF+zBzPW7sG+bjc+tTqWPbdfvLS6BZf/fi4emrNZe25cXxbfKe53IcseQQxxCkNA1Fgg5tQe\nQ73OwZSgixjcpGzLyw0xQknQ4nAeO7eKWpObW8mFuR6pqrLHGHsrgPsAXMo532s6jnN+L+d8Mud8\n8rhx4/p1T50gKw/2cW6cfoIR9yXJKXXoHIdjX1cf8g73FCYuWcdUpc9ibj27nQfiV0s2tXZh14Gw\nGV1cUyff5hO8xDpsm4Ex9/vKVhBdOn5Rl+2UCQfhp/9aif96ZEnoGNEGz10vwl/vn0t34OsPLyy+\n0RUiyrJXTxaiFTvacedL6/HtRxYD8F04OBAoGlwuZW/59gN4afXu0Pao8iViwWBLm75YLAfHgs1G\n73Cvf8lu0Trq2SJLEET/ESNA1Pyos+zp9hNErZOyGNIFuctL0FLYx3l8Xx5WcAHtiAmRIGqTqil7\njLGJAJ4A8FnOudn/r8zEZuNU9qsxauLYVOGlcZxgHbpbnluD03/6PPZ0ZJCyLFjMrVXnx30F7+0q\newyvrIm3WL7/1lfwzpte9D4L41iUZS/vmK2KUdjMlKDFHG/VmLJx/6ubtNcTbcglsOx986+L8cyy\n/tU/imPG2j144NVN6EwwcEW5N+Q5R1tXH87++Uue0luriCQD4jeU+0+g1qTUn/ujzH74zlfxpQcX\nePe+7snlaOvq09bZE26cAlNBd4cD1zz+uvdZ7fNCaYwPNicxjSCI6IWfTMGiZ1T2aBghqoCQSYpx\nkGpIWZ4sN7Kg7J0yYRQAYHijHSsnDmtwy4PFxcMTtUklSy/8FcAcACcyxrYxxr7MGLuKMXZV4ZDr\nARwC4G7G2BLG2IJKtUVGF7MXFZMWitnjwILNbbj2iWUACpY9SU59eskOAMD+7ixStqvIORzaJC6A\n+7LGFbOMI8qy5/BSIvaCMXtxpRdURU7fRvd/P2avum6cn3tgPv73Xytxyg3TsLs9OuhYtF2boMXh\nmL62Bdv398TGiUXRX3ege6ZviC0Q3FeoB5kuBEzKv4DpN+4rU2H1p5fswJ/mbsEvpq6OXDAQNBiC\nOjnnmHSIH98aXpwRyl50e8iNkyCGNtybNyMW8wrji0nZo3T0RDVhSoqWKDkiZTFvUVVY9m7+97fi\nsavehTeNag6Y9q54x8TQ+eIcUXaKqC+ifZ36Aef88pj9XwHwlUrd34QuZi+q9EIoYQvnmC65vKkJ\nTDKFl8nhHLblumjKtepU1ctirKjVGR2eMqJR61Q3zKSkCkXV8058UXWxpT1iEPBj9tz/q63syew6\n0IvxBzUZ9/s1EsPIbrJJv1F3Xw4nXz8Nd15+unSP4lbpVG56djUAYPPPP2Q8RihCqtWMgwet29Lv\nnck6aErbpTdMXFNyY9YtGKhKpSmlNAcwccww73PO4UhJzRN9Na7LUxppghjaiCEiaiwQXg495MZJ\n1BAmnS5q3kvblicDNKVdGaC5wcbkSWMK57on3/7J0/CuYw/BX+ZtDZwvztnfHaw9TdQHVU/QMtDE\nxerklJFfV3rhUEkxUBO09BVc5XqzeaRty3OFdBy9EGqV0bKnGwAcXpqriVtnT7ipyjF7ZjfY/d19\n5jY6Qetff79zOYlrSlRsR96RlL2E30nUaPrNy+tD96gkYqAXlj3xxcNunP6HPZ39T7WcyeW961uM\nBdyev3j2JO05UW5TGcnlM6u8n0JpjLNn15plz5ToiSCIyiCGgEjLXmFXb8Gd84t/mB/YT5Y9ohqY\nup0qv8qkbObNj7owCXFJizFtmI0IrenMkLJXjww5ZU+nrMjEZuN0OPKSgJnnwQQtQhDtzTkFy55r\nHTPV2RMJWnT0ZvP40B0z8Z6bX9Lu9xJsGNxRRFxhKYItk2L2gglswoOJ2H2gxzwIODz4fxIGKq4q\nTvGMenylWIiEz/vIJt+wPhBCg6/sud/XC85W7i+/A+tbXMX0T3M24+O/nZ34Xvu6fMW/szfnKV8v\nrNrtZf0EgGaD1dAYS8mDBY7V9zmb0LJXa4rVmT97AZNvfL7azSCIIYMYk6Lmx7xi2XtZja2vrWGE\nGOJEzWspy/JCOXRhEuJUxqIT6NXY1EkkpGJunLVK1MoHEB741ePVUgT5fDBBi1g56cs5SNsMloWA\nZU/1qWaFBC06bpyyEit2xCf98GPKgttTNkPO4ZGxdDqmfvu9XtvyiptqVuvGGS9gO4pCmiSSMM85\nrKLKhpZGlK7X3ZfDvTPcWDytm2wJSlpXn6vIDJcyRkYrlOUZXYULrRezV/jed7y4LrDSJ/fvja2d\nAIDrnlqR+D5b9nbhnFte8T53ZXzLXmtn0PprchHtzuSwYU8ntu7txvvffKi3nYMHLXvK+5kV+2J+\nlxoyLAMA2rrMVnGCIMpPVGIzgVhM6s3mDSEMJPkS1UPuf62dGUy+8QXjsQ2SG2daZ9mT3oNUpLJH\nfb4eGXqWvTg3TsVSoA7wNz2zKpBMwlWG/P3ye2BblpTkxD9exmJmwVO4+5nwlSy9ZS9tuT9vsbXg\nRAIMmzE4DjfG7PkKbPw11Qk1yTkDZX2JsuzdMm0Npq0Ilw8Q5B1/uE2qP3Rm3FXigLIXk/GzVO6d\nsQHzN7mlCoQi1OAlaPFbLGfDlBcHujN6d8ooZq0PVlHpyGSNApUcfyfT3pvD+bdNxxcffC2w3eF+\nhjy3rcHrisWWWnPT7A+ZXB6fumcOlr6xv9pNIYhBgxghotZCxZw6c10r3vrjaeH9FPtLVAHdAvuW\nvdHyYspmOOOogwEARx4cnnc9OYaxaMsemfbqkiGn7BWrQKjC5MbWroAlLyrWJu25cfqukGo2UCvC\nsmeKW+rN5nHKDdM8pdOPPQgeJwq/x7muqoj2iKLqsuCcUxRd3X11eG6cRSiItaDsyZmn9DGRUnKT\nhNqecOMc0ZDMslfMc7jyoWBS2/97ZjU+ec8cAJqYvQT3K6VovOrT39mbM04Qh49u1m5/aXWL97es\nKHLOI904hdIa1+4aM+xFsnpnB+ZtasOPnlxe7aYQxOAhQcyePI509YXnYxJ7iWoi990GOzqRWspm\n+Mb7j8NLV5+D4w4dEdov5lmLRZfGIl2vPhlyyl6c4qOO+zqh8e8Ltnl/Ryl7tuWukHAuJShRjo1K\n0NKb1S8bHujJBmKavBdeaYbtWfaKW34UzRFF1WXlNqcoukCysgG++6bLnI17zQeL6w+QdSYqMai8\nS9eaoGUvmQrR2VucG2ecsic//+dW+lbI7r5g3JvvwlGI2TM0N1hXsfil604lK2t3Nq+dIE6dMMpb\nkJBpVFxMZLdNrnxW3TiFsqcrsWJi+/4ebG6NXhWtJuL7mmoPEgQRTS7v4NcvrAvMm17MXsT4GmfF\n6G/JHIIoBZ0Lctz8kLbchIHHjAsrevI1rUK+BhPkxlmfDDnpQU6KoUN1p1Oz/QHA1rZu7281gYlM\nyrZgsWBhc/VYFlF6oTent+yp/tSmrGJegpYil2I8y17BBVV+BLIioCpwUTiGNgIIlLIInFNrS0ia\n5pRiffRj9vyVOLXfrd3dgUnXTsHGPZ2xSq+pDSIOTPQvoTTo6uwFrhdI1hLcd9dL6yLbAgAdSnKV\nfJ6HfveJY4bhn998jzY24PyTDg18DghoPKjsqYs3mXxCy5700p3985dw7q2vRB5fTbySGTEWWYIg\n9Dy1ZAd++cJa3P7cWm+bLzCbz4sbR2ptiiKGFvJaZ1z8aJRrJuDLZgzRpbFI2atPhpz0cMqEUfjc\nu45KfHycMB9Vxy5liaLqvsKku57JspcxWPZU4d+kdHlunIV7fvyMI7TXUxGt8bJxytY8TfxeIste\nhPvm5x8slFwYAAAgAElEQVSYH96I0lwIS0F3m7zDsW53R2CbbjB1OPcefNKkH+JZyAqH2oYnFm0H\nAExdsStW6VWf05/mbsFrm9uwr8t1pxzVnAbgZ6psUBK0hK4nZe1ULXu3SsKSCdWyJ9ciFDzwhTMB\nuBnCVIY1BBdk5KBzDo6M5N4cKr3guXEOnmCaPrLsEUS/EAtsfXl/7PBi9iLmr/gFRxJ8iephCrEp\nBVFreGRTOtqNc/BMrUOKISk9RK1aJHHjVPebjknZvrLnWfYcHopB0si7AMyFXFWFMV/IuBlK0FIQ\n6sUgcP6bgxYTE74bZ6GouiGGS7z0SRZ6/FqAyQck00TLOcftz63BygSZShPdR9OmW59bgwt/OQOb\npaBnXdNzDveUwKRxYGr8onvt4MWFwpWyWOyCg1qQ/Lonl+MTv5uDvV0ZAL6y16dY9kyXFfdrsC2t\nIBT1G25q7cKjC7cFtjlO0LLXlLa8mIG0xo1zRKPZ+s65+32HN7hWUfXdE3Uu455ZPcXskRsnQfQP\nsdAlLy6ZShbJxM3/ZOQgqoHodvJc3F+L2w8+cBJu/cRpOPu4Q2JKL1Cnr0eGpPQQtWqhEle2wImx\n7NkFYV0WPuW/OaJi9vTKnmrx+8ETy3DcD581unEKK0fSr828BC3BshFu24M1BsV3iMMrD5GsCQDM\nE21XXx53vLQe/3H/vCKuZkY3eIkMlntjUuI7gaLqye6X1wgZagvEd7ctvcIVOFazojd2RAPaRWxg\nwVImVrXFOG56vmJ7yra0SlNUPcXFW/eFtuWlbLRAMBZV14JhDeZAc87d/j+soBDm8g4enLUJawtW\nWKH46p7J4aOajNd1r117k5jjcLT3us+blD2CKI2cUmMU8MeeKOtd3KIRuXES1UTuf/31hGpusPHx\nM44IeBzpGEyZrocSQ1J6iLTsKZ/7Z9mzYBUyWgaKVgcsezC+XHI6fBnV4rd6V4d3LfX+gCz4FmfP\nEDF78vfLapTWJCs9vmUv+f1NE61wE4yqBVMMJotdknuUErMnhItAP1B+6rx0/zi3Cd2CxPiDmiR3\nUXGc+L3EPQxuwoUD0gZl72dTVmHNro7QdgA4ZESj9nqmPtKtyXAX9UxFnT1h/cvmOX78z5X48J2v\nAojOxnnW0WMCn/d2ZgKfTZb0anL782txzWOvAwAaKWaPIEpCXsASqFmsH1u4LVTeJE6wpTp7RFUo\ndLuvP7wIiwoLrAOV44B0vfpkSEoPcYGqMnF+0Gt2dWDB5jbtPjdmDwE3TkC17HFjNkjTAssHfj1T\nu12NX/Ite+79itWNWKFsRKDt0vMoRoGLc5nRFZU2HdtRsHSMiEm2kxTdbYQiZMtuP5pzZWGgL+ck\nUv68uohKP5DxhRMWG3+mK3QvX1PszuaDv4HRsufF9jHtMY8u3IYP36nvgzrFM1CeQuGUww8KbetU\nErzIcO66ajYXirGLMgx9OQecc68uoK7d8qucczjOuPGFQBbOjl7zfavFY5JLLFn2CKI0VBd2GRFa\n8d1Hl+LS38wK7tPM//99wQne30nk63kb94YWlgiiXHztzwsBDGSOA9L26pEhKT1EuXGGYqdiBO2e\nbB5PLN6u3ZeSiqqbapdxbnbjlFP5t7T3RrYDCFsCRYKWvGfhKdKyZ4WLqgfT8pfPjfPtP30ej8zf\nGthmGryEe+LIiNiuYtANXkJnkWUDncLiOP53enLJDnz94YWx9xPXlr9eKFa0cJBb2D76eiZXYz+m\nUih3TuC+poWMvLQKblJes4ZzxYKDUMbc65mFopRt4VvnHx/YFqnsFdrXXHD1lF2d93VnvUygOqul\nbuVz274e7+/2CPfUaiEv0MTVR1SZdO0U3PFifPZUghjs+Bltw3Mg5xztPfoxR2fZa0yH4/6i+NS9\nc71apwRRDuTFYRGHOlCWvYGqf0yUlyGp7EVZ9tRu3J+OLRK0uDF7/nZZyOY8mXPle37xcuwxaqKO\ntFJnr1jLnldUXXoGckHrUursRa0K3fXy+uA5hmcvLHvDK6rsiQQp8ZY9+fRpK3ZrjtLfT43dDN7f\n3WJbLNaVSFceZMWOdlz96FL32opl78VVu7FwS5tZkXMcMOb232L7v7hHsxR3pyZoUVEFsE9OPtJ4\nrONwZB3uKZOy6+UmyUr3zLJdoXN1CzdycqT2Kln2ot4feYGmIWXhT3M24/mV8X1McPvz8dlTCcJE\ne28W//abWfjba1vjD65hdG6cgu6+PE773+cAhL1pcg4PufKrdUCTsGFP7dbxJOob22LoyuQS1S4u\nB2TYq0+GpLI3ccww806lI/dL2bMYLMu1asgCnSycJzWJm+L3ZNTELV7phXwwdkvlW+cdp93OmKvM\nyMqGrph7kq/gJeWIOFZ15Yyz7EVlbSwG3W1yCa2WjpSNMym6WEe1H4j7m+LmZExWNvXawgK4rqUT\n//7bOUar9eKt+8G5a1Ustv+LazZJAtE1j7+OP87ebDxHtVidfdxYzP/h+cbj8w5HU2F1vafP/w7b\n9nWbTimcF94mK/PVck+Jcr+R39mGlIXrnlqB/3xowQC0iiCARVv2Yckb+/HH2Vuq3ZR+kZWyGzsO\nR5c0j23Y0+n9zZXFzbzDQ8pdY0payIqL6SPJmKgwKZvh2ieW4c6X1scf3E8mjhlGbpx1ypBU9t59\n7CHGfaoVpT9+0Laos+fwUGyXTLleHdWy5yVo8bJx6rU9y2K457NnhLbbhbbLz2DxVj+A/ZxbXgGQ\nbEITSVWiBgp1BdWkaJQ/Zk9n2dPE1Wlj+4r/9cQzeGrJDuO1Zcte3OAaF1cq9qoWwGXbg6UrVMuv\nnaDsg0o25x7fpGTU3NdtdpHUuSc22vqMnMIluqlg2esu1M+yGNCjSfYio3PtfH6lbwEcKDcYlajf\nT3bxpqLqxEBTTBKuWka8Y2nbwi9fWIu33DDN26fO+TdPXQ3AnRdkl3GBrPzpHguXYpTJ440oJ305\nB1/788JAgrSUxbChpTPirPJRircPURsMSelheINZSVAntf507LRteW548nVkoZvz8q3+hSx7SoIW\nkxcnA8OE0c2h7V5R9Qhh1HXRi2+bUNCivqoq9OsEjMcXbvPirMol/Gote3lh2Yv+crqC4fH3c0+Q\nSxiEE7RIJS4MD1icn40J6hP9S03komaeU7PUxil7y7cfCG0TbZFj9uJIa9yihjfaGH9QOLOnsGIK\nZU/E7FmMeS6dutp9Y0c0QNeNfz9zk/d3teYwdZFGRv5NiknQQhYFohyUGu9da2S90gsWnlwSjLFX\nF3leXt3ibi9sli15gBqzF77X1x9ehKN/8EzhGvQeEuVj6bb9eHb5rkDIgW1ZgYXvg8q0CK6jwbYw\nfe0ePDhrU/zBRE0xJJU9ebBWUeXm/ih7nmWPBwd9WbiLUiaKdQ8UddTk+wO+4mJKBGMx/T5WaHuU\ndbM3l0/Uyg7Psmc+JqUI6ep993RkcPWjS/HbVza47S5b6YUoy1742EcXvOF9dt04i0PXp8IJWvwV\ndd3xy7YdwGk/eQ5TXt8Zb9nj4prRSqGq7KUsfTZOwYfvfDWUZU60pRhlT5c0IWVbmPc/F2DsiIbA\ndjUBjFDwGPNdhXX3fvX75xlLTQiqpSBF/S6lytgkYxLlQMxbZRpqq4bs6v5GW09gn/r6iTFFjLtN\nirwgLzKu3Bn0jgCAZ5dL3gL0IhJlZL/GQ2bVzvZAJukPnvqmst7zexef6P2dshkyOQc//ufKst6D\nqDxDU9mLWCEvr2WvUHrB4QGlQbi6AQXLXsl3CKK6hwoLh2fZM0zYViG2UEVkLY3KSNqVyScSkjsy\nBctexLdNKY2QV1xnr28NZetUvw7nHFv2Fh8Ir/uJdYXPAWD62j34XqHumXdcgu//2fvn4aZnV7nn\naB5nSNlzopW93YXsrHe+tE6boEXGj9mLbqeapTaJC+k6xX1EtKWpGMtehIVWtd7613e3izp9TLLs\n6e7NWO0WSI6KuZR/kWLcTEnIrH0YYw8wxloYY8tjjjuTMZZjjH18INrV2pnxsj+LocW0UFgviDns\n0YVvhPapi0DCHVyMFw2KZU/2RLj/1WgLB72GRDnZ3x0uUQW4Cp+g3CV63nv8WO/vqEz2RG0zJJW9\nKJcUVSCMy4QYhW1Zviuk0bJnnhCKvXVGLb1QUJ7yXjZOVtge/v66l1gcFqVM9PTlE7VTrDxFu3Ga\nY/auuG8eblMyC6oCyGMLt+GcW17B3CKzUkVZ9tRdB5T0/G6mVf2Xuv35tV7q+5nrWnHP9I3G+6lK\nsLhm3tH3QbvwrHa398Yqe17MXoyyoHPjjFMQ5f3Ltx/A7wrfsWzKXkqv7DWH3DiBTDaPprQV+h5A\nshIW1SqQbPr9dh3oVVx9k0NhFXXBgwAuiTqAMWYDuBnAcwPRIAC4+u9LvSRApWZyrhZ5JfmKQLxj\nOstIRnn/xKKR+O7qvNRYRPiAPDfMXLcn8XkEoUOVP3SUO7ZblrNk2Vkue0TUPkNS2YtCFZLihN0o\n0jbzC5MbY/bKJ2KGlL3CJJWV3Dhf/u65mPODYKZDizGtAizcJHWZQD9VSI8/d+PeREJyu6fsmY8V\nQrpoykNztiQ6XrBo6z4AwPoig5X1MXuiJp0SS6f0B4ebYxbveHGdNvW9TnkL9TvHv7/OoiPakcvz\n0O8evrY4tjg3zmLLPlx292y0Ftw61aQGURSn7LntaRRunH3BmL2mtK2NTU3yXaqlIOlcZd9o68Y7\nb3oRrZ3+Sq7aF1/ftt+oKJbTsrdo6z5MunYKVu4Iu6wRpcM5nwGgLeawbwJ4HEBL5VvkIr8roh/V\nS8zedU8tx1tumBZetI0oE6TGugtlTywOqYujxVhO5Pfws/fPT3weQehIouzpYuD7g/zqyzJCe2/t\n1aUlzJCyp6BOBtv39xiOjMe2GGzG3HTOhtILPMK0V6y4llFWWkSdPXE/iwFHjx2OcSODiS8sFhb0\nRfsBoE9RcB760ln44Ftdv/BrHn89kWWvtcNVAqIOFUK/WEmasmwnpq91V0OTTLCiHcW6HDmF7Gny\nSpWQFeLcenP56BpyOvQxe/r7OAbLoVDcso4Tu8Lmx+yFr3Pjv52Cww5qAlB8ghYAeGzRNmzd65Y8\nkPt4MbWoGlLm38tk2WtSYva6+/J4aM4WNKdtrWDKWHT8IVA910edwqazTsvNW779AD561yxjuu1y\nfpVnl+0EoLdMLNq6L5C6nigfjLEJAD4G4LcJjr2SMbaAMbZgz57+WZDc2rDu374bZ78uOWD87TXX\nTdM0bum2hj1iguELanhB1OKUClnYiXLSkaAWrGzZO+qQiDJjCZHlAnkcaE+geBK1w5BX9q6+8ITA\n5yRum6dMOCjRtVOWBctyr+lw/0WRV/iyjlM2y16vskJpK3X2TAlNGNNP5p4bpzIZvu+EcRgmF82O\n+QINKQu72nuRyzvRpRcK7ZWbIrvpqagKkr8KHd0eFYdz3Dx1Dd583VRP8BYTvSmWTj5XFix0LrKq\nZU73DF5e3RL4Pp2Zwuoy1/dJ4Qqcy/PYkgNeUXWNH+OJh43EsEb3t1SVZNuyYhWkKa/vxGW/nQUg\nWPdQlxHTBDPmiQ27pOS8bJyFOntKn29K29r4UyA+5q1aCVp0yp7OOi23TljZtrXpawv2x/08dK1C\n83QLQpfdPRvn3za9bPciAvwKwPc557FFVjnn93LOJ3POJ48bN65fN7WY/y6Id6ZeYvY8C55hVtW9\nFn9VYsHFmCPeITVxGGXFJapFkgXJkWXKxvm+E8bhM++YaHTjPNATr3gStcOQV/auPOeYwOckSRAu\nOGl8omunpDp7jsNDKZwBV3gt13ygWnjSlnDjDMbsqViMafeJbTphVBbC49w4J44ZhrzDsbsjE/ld\nxQqq3BQxuOgUAvVS4qezGPCtvy7Gb15OVmSUc+CfS92ad5ta3QQvpvpSajC/GlOnE4h7c8HfRRc7\ndt1TKzBtxW7v3iLgOs+DyX32FQrPe26cDvesW+bv57t8qqRtyxjLKQoQx9Ha2YeXV7coLibJhcMT\nDxsZyngnUAUroeQ2Fd6lXkXRbUxZRuUx1rIXK1JXBp0VQvRDGVlw3FVIoDF+VJP2muW0Uoo+r/Zt\nqrdUcSYDeIQxthnAxwHczRj7t0rfVLboi3emXJmPBwpT90/yXogxRzyDlB1t2ZPfS/nvpGWJCCIp\nSfrv6GF+Buv+TAMPfeks/Oxjpwa8dOTcDmTZqy+GvLKnKmBJBuekq5wpm/kJWhyuFWijMvEVuyqo\nCv1ikhIrl6b5mjGmncxtSVlUBb24WkMybyoIpC3tvZHPV9xCXj0Sz1r3yNWBT3xkjOHppTtwy7Q1\n0Q0T54HjxMNGAvCzWvnKXvBYXQIfWSHK5BycfP1Uz7UR8DNGyufo2NrmCvgtHb3Sd+KBbKhf/fNC\nAMEMqapFV0U0z6S0i+eu9gErgesj4CZL+eKDrwW2FWMIOHx0M6Z9+33afWomPPEdGlJuu9U+39xg\nG/t5nOJqmkhX7mjHW66f6mVALTe6ZukUeLl5Qtk7ZHhD6DgAiLcFJUf0V3UM2NOR0R1OlAnO+dGc\n80mc80kAHgPwdc75k5W+ryXF7In/60zX61fSs5Cyp859ygKU/P7K1885xbv4E4RKd1/Om3uSdKeD\nh6VD2849sXRrv1zKSJ7XMzlK0FJPVK76Yp2SZHBOOvGlLNdqki9k43TjjIKrITnHKVsWQFOCFpEc\nJarOnj4bp5+gJVVY7RUT35sPOwjjD2pE3uGxSqmYHN1naz5WTJpySxjC2wTqbUU7ipVLHAcYN8KN\nY9xbSIhhitkLuXE6PFRKobsvj3++vsP7rLpZmpSOroLrppwQJ68IDG8U3PbkRYKbp67Wf7ECon/p\nFLeGFOu3ZW9Ygx1STor9DUx9U63BJ5S9lMWQsq1QnaumlD5mD4h3bTTt/cOsTejqy+OVNS341JkT\nI69RCvqYTE1cp9TC7kx0dtvyWvbc/9XfSMQz15siUCswxv4K4FwAYxlj2wDcACANAJzz31WrXTbz\n3/t6c+MU9GdOFZY7k7KnehvkHe4thMhjjDp2E0SxXHj7dK+80eaff6hoy544rz80SsqePC/lHI7p\na/fgxPEjcZjBw4SoHUjZU0jimpQ0M1nKFm6criKiS1oR5cZZ7DShunGOaAj+vDoXQ8AVzNW4BEDK\nxpnnSNsWMjknYP254KTxmLp8V+xqk3DPXLe7M9JVzlPWpKaI+CvdM1d/qlInVodzz1KpKi3qJXVZ\n3uKUCNmy15nJ4cXV+sR6+7r7sLczE1L2ZGUy5cVhJjfd7G7PoLUzoz2nwfYHctWyZ1ssssaiwFTX\nrhhMx5uycaZsC2mLhTLFNjfYxmvFvduccxzoyeJzD8zHle89Bh8qJCES59mmYMB+olss0T33gAWh\n8L8uDtM9tnxCpmMQeoVL8RiDdZGIhnN+eRHHfqGCTQkgZ+PM16myp77q4nVI4i2Tti1sbu3C6l0d\nAMJjkBpHLI8r8nuX5+UL0SCGJmod2yTZ4WXLXjkMCbJlT7bm5fIcn39gPo4c04yZ15zX7/sQlWXI\nunEuvf4iLLn+wtD2JINz0onPK6ruuXGGheJs3inbhKDWF5qoZGKKKqoeVXsvm3c8JUM+rCFloS/n\nxLq+inOvfWJZpOlfXEYXEKy37ClunIX/40oRqDjcV4RDVjg1Zk+TFCbO+iWnKP7ZlJXG4x6aswUX\n/2pGoP2cAx3S+X6G1eI6zf/7yyLtOWnJsqcrqp5k8UNXZkEXnxqF2bKnz8YpLHsqTWnLaFUU3+UT\nZxyh3e9w13K69I39+MZfFnnbc56yF/kVsLczU5Jri+4ZqxlwAcVFTCq9oaOcsUKeK58yRojnUi9p\n+YlkiAVKwO+btf4T/+21rZh07RTvs0mpi3stzjjqYHDOce6tr+Cqgsv8WUePCRyjpraXF2bk2+ZL\nyNRMEFHovHM+rsxnB5cpZk8gJ1uT7y/kkjfaSs9YTwwcQ1bZGzUsHTJ3A8my2CXNvNxg21JR9aA5\nXJBzzGsvxb6oalreow4ZHvhsEqgZ0wvOQrZz3TiDZREAV9nL5BzM2dga2S45oD2qsLeYGGULpOfG\nqc3GqZ7v/p+kFk3wOr7CpsbXxdXZS2LZ23XAj/XSFfWVae3sCyh7ec6xWYr/82onFplNpK2rT2st\narAtz3qqWn5TCWrTAeFFBgD4wrsnGS3JOpJa9l5e7aaVT9n6BYoGKeGMiinhgsDhXBvXmNSyd8aN\nL+DKhxZGHqO/b3ibzgorC7BCqTRZeIuJ931k/lZMvvF54zni3VAXAzxFIPGdiHrAtvzf1o/Zq+1f\n+Y4Xg8m4TNNM3Gsxccyw0Jh35MHD8Pl3HQXAnRPVcUderJH/zjnxC6EEUQw6Ze97F58Y+Cyya5cL\neTFPvv2+gizTkLKowHodMGSVPRPJYvYSWvZSoqi6K3zpAmeTxLwlRS1yOUmx7EW6cWr2iZX8bN7x\nVndkga/RttCXdzB3Y3Rd4MB9I76qroitZ9nTunEqlr3C52KVPYf7g1hPNhfaJ6MOtm5ZjfCXkif9\n2RvcmmknjB+RqH/JrokO59iy18/MKBTnJO4cMmnb0mfjTPnKkdo/LIshn+A+Ow+EE5ccOWYY1v/s\nA4nbZ3qnjjw42IdfWOVmLLU1LpyA22bTtS592+EAgDHDw+8h4P7WOuunKXZHh6gLKejpy2Nvp5/I\npK2rD5OunYLnV+6W7ps0Zk+6bmFyNSXQKUbIvPaJZWjt7DMKwr5lU7Xsuc+/xvUAokj0bpylXeu6\nJ5d7mY4rSagcgjE0IvrFkK2a8rXF/CPi8GXkd1B140wS80wQSVGzgQPhMIqUZeHDhRCE/iRm0SH3\n533drht/X87Bm6+bGpjTiNqDlD2FJHpXUreltG3Btnw3TpsxXPb2Cf1soZl2pe7JqOagUGuasC3G\ntIVihWKXc7jvxildRGep1CG7AUS5BYqJUra8RNXOczgwf1Mbfv3COgD+b7e/MAglhcO3zqlunKFa\nfpoELboJXbYQ/X2BW+i3uSGVSAjvk851HB7IepjylL3goB9X+2njni7s7Qo/lwZbFmRKs+ypvHT1\nOQCKc+8zHXrVuceGVi4B993SZay0GTNe67sXnYjlP7kYBzXplT3OudYNM0lGUtOCzUfvehVn3PiC\n93lNIQ7ovpkbvW06ZU9nuZXvISzQZmWv+N/NdI4pG6dv2SNtbzAhygUB/U/Q8qe5W/DNvy4uW9tM\nqHU9TX057lW2mJ94SNAgLYhZVnQJEvn6lKCFKCeOw7WLgM2KHGYx4K4r3o5Z156HGz7ylrK2QZYH\nVC+ll1aTslfLkLIH4Hf/cQY+cMphiY9PusrZYFtuZrOCsmdZTJv1slys2e0KktdcciK+fcHxIWE7\nMhun5kvJbn9pz43T36/GU5lISe5vUcpDNu/gZ1NWBtwCv/iH17Bo6z59zB44vv7wQvzyhbX4zcvr\nvYk1zlVSxeGIcOMMt1Em7wA6Tzo1yQrg3iOJFTd4blCg92onKg0brombkzHV4muw/Ri3UOmFhDF7\nKseMG1H0Oaa+mbYtbYydmzwm3DbbYl6//+Tk4HmWxTCiMWW8V2cmh8/ePz+0Xaym6lw8BaYYSjXA\n3svaZxAQ/etFJ2gRixKmNpWm7Bm2G3bUSzwXURxWYc4C6icuU12sNPX/uPFXNw82pvzyNK5lL7hf\nHofk6+fy4Tp7ZOkjkqL2lTzn2jlPXegQ7+qE0c3aRfx+tUnq3/uURfVaHyOGOqTsAbjklMPwuXdN\nSnx80lgk4f4hVvhsxrTnqrF2/eV9x4/Dty84AYCbFl9gjtnTX2fXAdmiFHb1i7MmqecC4QFswuhm\n7+8Ne7rw+5mbQgrXTc+s0radcz8u8ZZpazB3o+suub9IN862zgweec21vqlKkSo0qMlfTG6cuiQx\nbmbN4t045XMsi2F9SyeyyvV1yX+S4Lo9un+rCxGpEpW9OB676l3hdkRMFLp3Ji0tIAT6uOXbmUyx\neaZbtXYEJ69nl+3E1r3d3gQb5TqbJGsp4H+XgOtX0tILsrIn3DgNbTLJtHe8uA6Lt+7T7jNa9sRi\nhZqsqE4zNRLRyImZ6sUypc5FaqvF57ivoxNYG1KWtxDGWPiYNbvaMenaKZi1vjUUv6c+v74isigT\nQxvVuyPv6JW9gVSygm6cQTmLZoHahpS9AsW8L8W5cTJw7md81BUvV91G+otJIYtK0KJjl1REWljn\n5GOTKnvy6pI6WCVRnF/bvE/rgsg5x1Fj/JguUadOWD1Ml1ZXd//3X36GzJBlz4merB2DAqeb1B0e\nXunV0Zf326Aqe/M3teGC26djxQ6lvlyJyh4AyUUpnI0zl+dlD76ePGlMaFtUL9ApbXK/kbOPWcwv\n12Ex4DsXnoDHv/buwLnGBC5Kv/jaw4vwgV/P8J5/pGUvFzy3U5O0BvBdZU3p2r3raZQ4+TjZjbMz\nkwu5n5qE9NufX4uP3T1bu093SktHL9YWPAbiak4SgwPXsuf+LfppueLKK0Xa0lv2pq/dg+37e/zS\nCzHX0a0Pua7uYn947HhhlVtK5zP3zQssFrqlF0jZI0pDnQNyDtfG7A0k8hzZ0h6M1e/ozVGilhqm\nYsoeY+wBxlgLY2y5YT9jjN3BGFvPGHudMfb2SrUlCcWsSiR14xSlF/KFuC7G9Of++cvvKOLu8QQS\nnEjb5flQbodJ+P3G+4/1r6krvaDMjMeOC2b/1LVHFRiLydio4nAEvqCYSOMGHFMWTyCs7KnHqklB\nWjszWsFeWN6OGes/E8dgBVSJsuwJtu3vDnzW1XBU+fgZR+DUCaMAAFedc6xXbNVUVL0xZSOTd7Cj\njIsRBzXpS3tGWYdUNxV12+Gj/YKuNmNeDJnNGL51/vE446iDlXvp76PrN119eU+picokK6/CPr9y\nN065YRoWaSxoOjdObVF1bZ09/7jePj8b5yk3TMOld81Sjg23Ub3Pvq4+/FRa6BDX783mvUQz77n5\nZatIAl0AACAASURBVC8brNqkSlh9ieqjy8ZZ6790OqW81IUGf/6B+bjklzOMyuqDXzwTK35ysfdZ\nNw7JMXu6MAw5plosjADCshc8VvXIIAgTalx+Pq+P2RtI5DlATcz29NIdoXmIqB0qadl7EMAlEfs/\nAOD4wr8rAfy2gm2JpRhTeFK3JeH+4ZZe4K5lTzl34phheM/xY/HWI0YV1d7I9mmyWQLBdgdq2Rmu\nc8TBw3DCeDf+Sljn5MlupCK4v3j1udrryEqEOlj1Q9czKkK9BSuHw4Hv/H0JWqVsiEB03GCHktE0\nzo3zkdfewD8Wbw9dRyiechIbN/Oq8dYea3f7cV5qzJ63XXmOadvCq99/f+R1v/DuSV7fkJ+7btX6\n7OMO8eoobmkLKpalcv/nJ2PBj8K1LQGARYxEKU3JA9nad9go3xWYMd8tVWdFj9re3ae3xgnrbpSg\nJvfrWevdUiRLtu43Hh8Xs6eb1OV+0K1k4xQFoP1r+seu3NGOvBMsK3GgO4vTf/o87n91U+icKa/v\nxOcfmI9nlu0Mxo+a3DjJP2RQYcnZOPOlW/YGwhrYl3OQzTuhhUf5nerI5Ly+rY7nY0c0YnhjCs//\n9/vw5DfOjlD23L91coKcbVf2sMhp6uyRZY9IitpXDvRkQ4vRgpnXvB/Tv3cuXiwkRys3S6+/CIuv\nuzB2sXrN7o7I/UT1qNg0zTmfASAqJ/+lAB7iLnMBjGaMvalS7YkjSn8TyS8aUxYmjG7GhSePT3TN\ndCFNs5yNU+a4Q0fgyW+cXXKbZT7zjone36YkMCZlL0pYE0qeUNhkQfncEw9N1DZZMFcFxv5Y9jjc\nCXXM8GC9xJ4+f5B8YtF23Pbc2sD+KItEe4+q7AX39yUsmi2sRE1pObNoshgYWQDn3HXdUBXrnMND\nlrgjlDIFKmnbF1qChevd/+Xf4uGvvNOzFv7nHxfEtjkJTWnb6PobtYCiK3kwUXLfTVvBfi0EsmLd\nluV+I+NZ9gqT7wsrd+MFJc20zsVT90v7iS/c4/88d4tXvFlGJxTKCn6ca+mcQrkPAPjgHTPxu+kb\nAgsVf5q7WdM293+h9H794UWB/arwTtk4Bye2lI2zP5a9gbD8nnjdszj/tumhRBRqiQXRFHX4FZl5\njx8/Em87crR2zGi0bd/VXdPVA99T+tPhPGQNr7Zlhqgf1L7yvltexrLtB7THHjlmGI46ZDiOLSE5\nWhJGDUvj4OEN5M1Rx1RzTXYCgDekz9sK26pClLgiShgcP34EZl17HsaOaEx0TZF9U1h0VIvCieNH\neopKfxdB5TILuqLkbnv8v6983zH+9ihBWyh7dliALiVBizpY9Ce5Ay8UwVaVPZOFBnCF4z/P3WLc\n35HJBeL0QiuzCd1wRPygnObfZIk0IRIl5ByOc04YF3CTzTscB0m/eZLH6LoVh4UWk4uSUPZMltDT\njhyNk990UKLvIt9Hvy/iPGXnvP85H7bFMLpQt1Lu77alV16T3EutsagilL6vPLQAX3logXafe333\nBjrLhuceV/j/nhkb9PfSKXsx9RxlfvRk0Ht+9a4OTzG0mN5aLNrbZxBI1VNy0nehLIODB7nOnvhd\ni5mfMrk8bnp2VdH1TkuBc2BrWzfSqbBlT37/xDiutkldRNPF7DWmrcgFJDm2Snb1zlHpBaIfRMWI\nVwvqz/VLXTjgMMauZIwtYIwt2LNnT/wJJd3DvM8uKCulrGC7bpyuQKQKmeV8cUZKSkXQjVNqi/Th\nuxefiLOOHlM4xvy9Ggrf3XPjLMESJwfPq8Jpfyx7zyzbhedW7kZT2gpkZDS5OgDAvTM24sYpq4z7\nOQ9mRw3F7CUcgEUcnyxMJHXjFDSlLOSFVdgKZnLN5p1A/FuSp5iWyizo3HtVpUooeymL4YiDm6Fy\nyuEH4ehCTGISN+SUJvZOEPduiSKxDbaF8Qe5MXqzrz0PK35yceC5XPyWw7zvY+rWJqXT1G/E4VGK\nvqycRY0lovsLRUmnLM9ct0fr2qmv55isQzWnLa/9KcvSKomcuzUqTUKGeo4Yv7bv78F/PrQAX/jD\nfCx5w+y6StQHFnOTinHOA5a97z26FK+saYk9/5H5b+Ce6RsD9SUrTaNq2VMW1kzj7ghF2dO5eDcE\nPCLC15DfF3kcyDuOJj6chGUiGbWp7FW7BUSpVFPZ2w7gSOnzEYVtITjn93LOJ3POJ48bN65CzYlI\nEOFloiz+qmJyyBUEdvkS8rivup2Y+Oo5x2i3ywkrdG5vbluY8tn9P+priXiplKbOXlICpRfK6Mbp\nX8PCsAZ94g8VOZjeRLsUt6e2tzebbAAW15Ctb06RK70p28Kr61qxu723oOz5r2sm5wSsuUlIS1nl\ndG6cugQtgKtQ6Or1WIx5v21TKj4baJRlL+7duuKsgpuydNywhhSGN6a8PvTVc47B5EljvEPM7sz6\ne5iUPaFkrd3dgUcXvKE9Rla6vMx/8vutxAz5Fr7gdTjnXq2/oPVVb43TJZURJUhkmtO2p+yZahQu\n2LIPb/vf5/HMsp2hfaJt2/Z1e9luZTejF1e34JU1e3D135dozyXqB/GeikVKwP3tH124DV/4w2ux\n5yf1fignITdOrpQ3MYy76nm6McrSxNrLyEJ5NqDshe9LwjKRlCQLeZe8JXl96HJAbpz1SzWVvacB\nfK6QlfOdAA5wzvVSxgAQJWymPMte8QiBM5t3QhNGSat8hlNMiVjkVqvCrxCSo9043X2pBMear2G2\n7JWjRlfaYoHYuCiSDFayq09Y2UsWs3egUINGduM01eQzYTG3KHdv1kHKYgFlLJNz0Ji2cc4JhcWP\nBM9RVrqDCVoMlj3pmeoWEGyLeYsATTFF3cXxJmL7QcTChLiu2keLj9nT/7ZCaHx2+S5877HXDcc4\nob/l31pNZW+qYyZ/FILoaUeMwiEjGrVunLpMsJ++d25oW1OD7QmiKYtprYSi/t7r2/RxIXmH4z03\nv4yTrp+KhVvatGnAqbBu/SOG68/cN1dS9pKfn3ThMgl3vbQOSxNYi9Wxi/PgWJ903C12gQhQ3Dil\nv3OOo1H2SFgmkhFn2Zt17Xn47X8MbBJ7ctevXypZeuGvAOYAOJExto0x9mXG2FWMsasKhzwDYCOA\n9QB+D+DrlWpLf/EsKpqJYOkNF0WeKyahXN4JCbvya9PfOUD2kDMJ1WrGQ2FnTGJc08XsJUVWFMrp\nxilfozmizpzc5CTFr7skAVod20wKgYpQGGU3zmLHSfl427JCz6o5beMTk48AkNyNUyALR6ai6nKG\nO12tO8Z8i3JzAmXbZHGW22DC76ualXcRc+i9p4XtRfYtU6xnkqQKutX9QCZLRblTlT6BLAyK63C4\nz0fNwBrVZpVG2/Ks0ratt+zF9U/5nM/dP19f4DdRa4haRrw3cze2eZZnk5Jy5UML8NLqYLKicsmD\neYfj1ufW4tLfxKdz14VI5DWLLcVeRyAWMXSLGbJrvzy/6Iqq13q9QqJ2iFP23PqPAzvi0mJF/ZJI\n2WOMnazZdm7UOZzzyznnb+KcpznnR3DO7+ec/45z/rvCfs45/wbn/FjO+amc8/Kk/CuRaFdGs2VP\nTfmsIgTRrMNDgqo88Pdb2ZOTVGhc9OS2qPuixgtfkI6Og4oiUtkrh2XPthIXFU8iuEcVi+5JaNkT\n11Bj9qKufcNHgq+ZvIqmWvYAN9NnqojfRe6run6hKpOyZa9BE29nS26cJmV76Q0XaV1HVeIWEfxr\nmPellZ0mwc0kcJlcdJMsEMgr+iLrZa+UudVL/e6Iaxose9Lf4ud3OAeDXkETiYDiuOOl9fjIXa8C\ncPuStrZfjHAhu7l29eW1lkbxO7a092Lexr1DpshuKXNkrSKPyaIP6/qL43A8t3I3vvRgcOoulzxY\nTMxSaG5FcHEkqZJlWiCKGp6CCz1qgpbgsSQrE0mJc+PU1Z+tNFFlqwRk/atNklr2/s4Y+37B5bKZ\nMXYngJsq2bCBJmqFxHPj1ByiSzox5VvvwYNfPBOAFLOnceMs58AfsNRIv6p8R5MbZ9R3txVFt6QE\nLRFunKZbD9e4BV52+gSMbAzH5qVslnjgi1K2BFFCRlJlD3C/g/y82rr6jC5yQNh6Jg+saoIWwC1l\nUIylNejGqbHsGWL2dG0Tx3tunAZlb1RzGsML8ZRRfSehF6e2r4rkPN5PW/jftJAgHuunJh8Z2G76\nbZO4/uY0rlxyX1OLVIvP6uKDbuWUc/c30u1T3Th/8ITezVTGNih7+2OyJ3YpVkSdpVE88hdXt+BT\n987Fvu6+2PYMEgbNHCm/p6IUjdz3Zq1vxdNLdxiTVZXLjVOtaRq6j9Qm9VV3eHBh7bXN+xLd0zSe\nRsbsSe2U32dHk6WW5GAiKXGLHeXwiioWdd2zUZORvRgZiRg4kip774CbTGU2gNcA7ABQngJxNULU\na+MlaNHs07mmveXwUV4NOvFCOlykhfePlyfQpHOA6ThZsNUVoQbCE2JcbJO7T5yrt/4kIVB6IWGC\nFt32Yw8dgXEjw2UvUlIA/TCNkthaSMrym5fXawugq0QJ98VYKoY3poqyXKrH5lXLnqLQNqftouJJ\n5b4qHy9+W7Uvy6U1dP2cMeYdE2VZFUJZtLIX/Q0sy/w9RaIakRRHCJumVXrxWNXXJC5mL4qsNAuK\nWoyBunhKghZhLVTfB90CkMPd56O37AUVsL/O1yeQkUlZlvZa+7r0itmbDxsJIPx81ELuMuL61RBI\nqsSgmSPl+UC4o8v95TP3zcO3/roYGeEWrImXKwdxwq6sbOrcJb//ePTCh65vqts+eKqbAMOfB3Xt\nDMbp+X+HLXsO58jlHfxsyspAMXaCUInr/7qkaZVGna/U0iUAKXu1StLekgXQA6AZQBOATZzz2ssL\n2w8iSy9EWMBihVTZghIRs5cUzrl2kjJZ9kxtcT/71zThxSoon4shVULphS6N4G1bTCvtpwrF6wFg\ntCZD5XMrd2P7/h7cMm1NovZGuU8kzcYJuGm9i7K8RQhNts2wubU7sF+27CX5XRhj2tg38Zep9AKg\nt2Dblq/kJRHq+yP4+/0vvE8oe0IwFc/N9Ox9y0Nwv8lSYVL+c3n9ir7oP4GYIeG2qcTqqRO6/Jsf\n4tXgDJ4baEMJpgLXshf+rnsNyp7oW2q2UpHQRXeseDZp02A0+Bg0c6T8noqFMl3fy+Td/lAOV3wd\nscqewaIGuO/RS6ujy0TcefnpoW3yVzl94mj8+tPuMVHjuPyey3PH/E1taFPeKYdzzNvUht/P3ITv\nP74ssn3E0CYTI2tExcBXCnUulEt+mY4haoOkM/FrcCeyMwG8F8DljLFHK9aqKhBV56s/2TiD6dOD\nV5DfiWICt3WTq7zNNPmqwrYQzEQ7fvyRk/Hdi04IHOOdIlz9SngI6RKKqusGDNmCJ2PbvhKoG3wA\ns9VCR5IYrSSMaExplZPJRx2MpTdchHcdc0hgu6psBdw4GcP2/T2B/Y1py1Oki/1ZtDF7SmNlN05T\n6QURq5d3OCaMbsZ1Hw6FLnn0Ryj040vD1xjV7CpF+7uFZa9wP8Po5ln2lEvlHY6mtIUXrz4nsN2k\nUAUSM2gStPz9Nd/KpiZkEf+raeodznHUIcNwzNjh+OZ5x7nfh7sLOOWaRFOG0guqlVAgnpNq2dMp\nx6t2tuO/HlmM1oLVIqq24iBj0MyR8jjUWhg3df1F9F1Vn08as9ObzUeWwsnm/Ovc9MwqrNsdtCTL\n747qcZGkCR889U2hbfIYNaIx5Y17UYtp8nspK6j3v7oJ3/jLosCxnPvWkM17u+IbSQxZ5BJQOmrB\na2KEJqymlAVIovIkVfa+zDm/nnOe5Zzv5JxfCrd0wqAhsvRCP5KTmBKnAMkUPPWenOvbEbiP7K4X\nEOqD56iWvS+cfTT+33nHK8cIRbd0N075nKSWvQtPHh/aZjGmTbqRtphX6mCCpvg3ADQnKA0gSJLE\nRaUhZYWS9QxvSGm/31GHDMeo5jTOP+nQwHZ1pU5141RpStklD/hBi7PYFjwm6MZpUvbc7dm8g1nX\nnocvv+fo0HGe8tUvwV9YJMN7Dh6W9toA+P3ZuBov7T/7uKDCnbbCv6MpcYm88pqVfisRayRbp/1s\nnPD+n3TtFGzdG7TWOpwjk3UwedLBaCwo0hxugpZyKXu2xbS10EzuN+IxqjF7Jgv4U0t2eDFS1XA1\nqhKDZo6U5ynRT3R9T+xTx4akvfQLf5iPM39mLrzel/f74z0zNoZq/MkxfarV+fLfh8uPJMEcs5fs\n/LgkR5z7Y8AOZfGuvxzozuKmZ1eVvRh3d1+OXE6rwIGYGOpaKHMjwmo+JC2cxL0DRHVIOhO3MMYm\nyv8ATK9kw2oJ27OeFP9yBerfKTNGkjS2WktWjBunfE+5zergIBfPNSHqxJksKxdplDIVueB5Esve\neW8+FHddEXaxsQ2WvZRtoaWjFwBw/PgR2jYU405Z0mTJw4kJ4tw41d8xqlSArVG2mhvsyFiSKHR1\n9hhj+OCph+EPX3CTC8kZNnUJcCzGPIUk6pmJbl4py97kSWPw1fcdg/+77NRQ+3SILsgY8OAXzwqU\nT3GT/QSftUkJyuSiLXsyasyeoEOxpnG4GRCb0r4i7xQWeJJkQkuCbTGtO7IpHlWMIaai8zpEYo9q\nuBpViUEzR+p+Mq2yl9fH7CVNzz53YxsAoKM3i10HesPXzwWvoy6OBix7yuKF6j6ZFHkMlm/nu8tH\nnx+XQdHh3HOhLuZ9SsLPp67GPdM3Ysrr5S1X/NG7ZuGMG81KOVEZ4pS9anDkmOBi+qEFZe+Ig5vx\n60+/DQBZ9mqVpMreFAD/Kvz/Itz6eM9WqlHVIIllz6Tr3f7J0/CbK/TFLaMte/HtUideDr2yV4og\n7St75oZcc8mJ+Nq5x3rFu9Vb32X43jLDG32lIc85zpx0sPdZt/B//+cnB1wIBWlDXZmUxTz3zePG\n6ZW9YurD9OacopUnDh76PU1unAL1d4z6DXXucE2S5U0I5EtvuChgkTMRWBwQ/zPg7s+cgfe/2bU4\nytbQY8YND1+D+eUcklhD++N24luY9df9wQdPwoTR7kTEvXP013Iky17atgKZX1O2FXrWpjhNWdjU\n1dmTMdXVU+GOq3Q1pW1v3HG4W7alXJNoV18OG1s7Q9ujakhaLHlNP8B3qa0FV6MBYtDMkbpFJ91C\nQzEJWta3dGDStVPw6rrW0L5LfjUT77zpxdB29T1Sx37ZjbgvVx7FydRdE1v2YkIAHM5jFcJSEUpk\npkzPQrC+JTxWEJWnFpW9x656t7cYDPjx8g0pP6SEYvZqk0TKXqEO3lsL/x8P4Cy4BdMHDVFWO7X8\ngMplbz8Cbz1ilHafPEmoQmQgG6fh/dBZJ4TiFWyj/vwoZcOShEkTh4xoxPcveTMOG9Xkfh4ezIaZ\nRLGQ/brzDg88a50waHJPSNl6N86UzfCnL5+FX33qbd7go/LQ7M2x7RRc9+TyojPKnT7x4NBzHN6o\nd7P0skUaSmEIZKuI7jpNaTvkMjWqOa31o1cJWJwNWVnlzKbHjhuB1354AY6VlL6cw5Eu/P6m5CaA\n/33Lk6Al/hqeJdFwP/W3lb93ymKJk4rIQpUswOmUQyEDxi065DlHb9ZBU8ryyl04DgcDtElVSuGN\nth5s3BOOF4rKomYxVpQl4kBPFmmb1YSr0UAwmOZI3aJTlGVPV+MucD2LeW69/1y6I3QdNRZZEOdh\ncd/Mjd7fcWUakmIaM9R+vPSGi/DkN8LJVuMWvba2dZcUJpAE4f1hUibbe7PY3EpxgvXCgZ7ki2sD\nxfiDmrzFYAAYVQih6M3mvXen3G7ERHkoKaCCc74IbqrpQUOkTJLAVc40SViKUC1fQ5b7TLWJdMre\nrZ84DZefdWTscXH4MXvxx04+6mD8+CMnh1zlkjBcUj4ch0cWeo+iwba0x6csC0cdMhz/dvoEo/L5\nxzlbkje4SJ76xtm47/OTvV/w8IJiPKIxrf9+hQPVGjUhS58VVEJUmhts/7eTdquC2a8//TZ863w3\nFlMXn2cyXMvujGnbwriRjTh8tO/G0Zd3vELm5XDjPHbccPzoQydp97EE76B3P+8ckxtnUNmWD0vZ\nDOlUsj75H/fP8763LGzqLGQmN04VYS1sTNueG20m5wAMqPQcalqQ5eCwLBZp+VPpkSb/oUg9z5G6\n3y06Zi94/B0vrgt8tpj/7hfjYaEuIKlrMH9fsC3Ulv4ijxnynKzGro9qTgfSzotHEGe1+69HluCB\nWZsij8k7HD9+egXeaOuOPE5FuNqbYqYuv3cuzr31laKuSVSPAz1ZHDwsjY+cdni1m2Lk0JGurNPS\nkfH6H1n2apN4EwAAxth3pI8WgLfDrSM0aIgSIj2rQoT1z6S0RAnspbhxAq5FZ+KYoFudbTH838dO\nxVNLgnXkosStJG6c3nUYwxfOPjr2OB3DpZi9nKLsFTP5p2ym/Z3k56pz/6w0px05GoD/ewqLzDAp\npk5GfGM1aYzaP1IWgwiLty2Gmde8H+/9xcve/saU7Qkk8plnHT0Gz6/c7X2+9G0TvL9Ff1QXIYDo\nmEExkN/x6dPxmfvmYeXOdvTlHO+7Rq1WJ03Q8uLV5xr3+SUjIi/h3i9BXT/5WqywCMO5m6DFVKdS\nZXd7Bptbu3D8+JGBeDedu6OYAKct3x3aJyPObUxZnmU1k3MwvDFVNstesXzunZNw/dPLi44xGkJl\nFwbVHKl144xQ9qIScAHBRc5i4k5VBS5qYbBcyp5pQUo3lAQ8AmwLfTknUSbnuJIQy7YfwIOzN2PJ\nG/u11kMTYtwyuXuv2NEOwF2YG0KJk+qW3mweJ4wfif86/3itRbyavPf4sejpy3sxe3s6MpJlj5S9\nWiTpGz9S+tcINy7h0ko1qhokUeQiXSKNvv6SUB2qs1e8G6c4TlWSLIvhindMxN+++q7A9ig3KrGr\n0gsxsvuqiD8SFGMtMMXsyUpEErfSSiN+5rTBEimUkWalELnaP9SFgiPHDAvsb0pLMXvSqXd8+nSj\nhcx32fS3qbUUdQjh4ODhDfjUma5V2RUaCgN8gk5UjgQtSSzBfp09/X7x7mjLT1hMm4zGhPiNMgFl\nL9ynHc6xelc7pq7YFXk94UrZlLa9xYDebL7gxhkT71euatYSX37P0fjkmUfCZqyomD1gSJVdAAbR\nHKnzItD1vYxG2dNZ+BmT3xMH983ciKdjhNd9XX1Ytu1AYFukslcms7dJB9LdW96SxMNBx52KFRTw\n3+Ni3+dUQstKew3GghFhHMetqVyLSa7+9OV34LGvvRsnHjYSAHD2cWO1MXu92Txa2sPJl4iBJ5Fl\nj3P+k0o3pNpEWvYSuJCZlCpLEdhlpTKRkqVcViiISWsZRVGMZa8c9wHClj1VMP7exScar5M2xOzJ\n1jzVNbIaiN/cba9G2Sv8ryp78qA+silYtkG32t6ctr2BVe5XzQ023nzYQdq2efGncsZWycJlIqW4\ndAKuYCMStGQTrKyXI3wrySVMMZHefk8ZlJ5B4f+UYUHBhHiecjZAXfxQ3uGYvnZP7PV+P8ONQ2pK\n255FvCebT5SgpRJZ0MRvbTFW9IptaghZDwbTHDlMU6YmKmZPHrd0SheDPw5OWbYTU5bFZ4v899/N\nDsWVym+lqgjFFaBOiund120Oun9bAPJFvyO3Pb8W3yy42DuO6y4t2lDs2yzGItM4MLIxhY5MDvt7\nsjhkRKP2GKLycM6x+I39OP3I0dH1G7mr7MkeQE9942xc+ptZA9HMRIw/qAmLr7sQo5rTmL/Zza4r\nW7f/86EFmLmuFZt//qFqNZEoEDkbM8b+yRh72vRvoBo5EESJd8ncOOO3hyx7coIW43X1F1bH81KU\nP9GcciiOSe4DFOoESs9RTQox6ZBw1keByVImW7hqQtljQtmzQnEmMk2KUCV/t+f++32BcgvaOntS\nghb1sZju69WMDNw3+L+OtLRTWE8zOceL4ZMzrKp89X3HuNfoh/DvW9eLV8RUvO6usewVY9WTz+vN\n5iP7Xt7h6OyNt4w9ucS1eDSlfTdOUVszzkWsEokfxPOIcvE1nluDK9LlZjDOkXIsmkCnQAhrttw3\nTIs+xfYfXQIhU7w7UEbLnlHZC3v3yPNYXLxcHPfN3IiTrp+K1ph6drPXt+Iv87Zq94n3zTQOHFRI\nXiYy5RLV4fFF23HZ3bPx7HK9l8fu9l78c+kOV/lnQWXvtCNH422FsJFa4eDhDbAkC6Tc/2Zqsu9G\nkXd4bDF5ojTiLHu3DkgraoAkVrtoN079TnnyUAX2JDqWKa21ao0rJShWruNVSVQBXf6ounFG6QMp\ny9IqJLKFrBoxeyopz7JncuN0/29S2ip+j0OGN+BNo5oD/UWX4bEpbeOYccMxYXQzrr7ohMC+uBhS\nXcxeVP9OS4qMUPayeY5JY4fjle+eG3Ixlbn6ohNx9UVmi20SikrQorHcyagJWuTrFusyI9673mwe\nIxpTyOT09b0cziOzXao0pexQTGc+RplLEi9ULL5lr/hzh4hlb9DNkSMawxmNdQuCv3rBdUGUx1+t\nZY/1z4VbIL+vamsyRbxbkfeQumywzl74WPkrefFKJU6mN05ZBcCNfRKX1TncXHHfPPf/d0wM7bO9\nmD39OHBQcxrb9/fgQE9pNQiJ8rCupQOAm5lVx+cfmI/VuzpwzNjhOGREI4YpHkBPfO3dFW9jKYjx\n/nMPzMevP/22QK6AXN6N79/U2oXFW/fhsrcfob3G9U8tx8PztmLdzz5AcaVlJk7Z28Q51y8jDTr8\nkfudx4zBSW86CH+Ytdndk2CeMgmWsqKjTniiCCVg9s83CVnq8SZXzCRKbCXdOCeMbg5tk5+JcONs\nSLkB7lFxGQ0pfSr3JlnZS1d/gBBtTMW5cSrCvPhNJxzsPjNZ0V+zuyN0nea0jZFNacy69rzQvlhl\nT3pMor26c846egzmb2rDUYf4ylxDYRVb1LaaNNZsjS0XJSl7CRO0uH8XfrMik4rkPGXPTaKyI/tq\nLQAAIABJREFU11DMOe/wohKcNEqWPdG+uOQWlbHsuc+jlMyatRhrUgEG3Rw5QmPZ0/U9UTJBVvay\nufBxDNGLeDoabCucjVN6+dU5q1yWPdO46ZVBke4bcOMUilY/25GXwhyEO/qKHQfw+MLtxjhsv43R\nbpwjCvVuOxJ4GBCVQyycmBZAdhdi3Lr6crAYCy2aleJlMRDI4/0DszYHlL1sniNlAx+6Yya6+/JG\nZe8fi90Eg5kcJREqN3FP80nxB2Ps8Qq3parI792x40bgOxf6lhLhrhGZ7MTwJANF1aW/f/Shk3BU\nhMui365kbpymOSZJUplKJHYAgJv//VSDIuL/LbIYCoEh6hm7lrLw9mA2zuoNEN++wI29EM9TLhVx\n9Njh+Py7jgrsV2P2WjpcFx6hIMuWvi9pMqE2RSi2cQmDZBckT5HSHH/Hp0/H9O+d66VYBoKWvYHC\nz8YZP9GJJAkNBpfMi04eDwD40Kl+SmvxvNSkInG3E8JfTzavjXWSj+stxrKXtgOWX8bC1ns1GZGI\nmUjKRSePj31fhHtaKfXyhkiClpLnSMbYA4yxFsbYcsP+zzDGXmeMLWOMzWaMndbfxiZBV6czMh5U\n+plNSleSKUa2Hg5rDL9LUVmcyzUWmcYX4SIpG83kd8J34+xfO3JKHVoA+N6jr+OBWZuwYU90cXMx\nD5q8fMR1KxyiT8QgXhGT0ia8kzp7c0UvklQTWb5tVmQTkS1XLHiaQoe8uFOq1Vd24rqS3BuPqWRD\nqo38Rd1U0dHC8NgRDYHzzZYU+W+z8GOO2dNvD7lxljCCf+v843HhyePxb6dPiD+4jMhfSbRaCMpR\ni1auG2f4AHkFKElB8XKRslhAwP/2BScEApFTNvMsaE1pG2dMGgPAnKDl3BPH4ZwTxuF/Pvj/2Tvv\neDmquv9/vttuLyn3pvfeCyEJJKQQhECQoiBNpT48IMgPUQTBgoIY9UF9rBFQsIDgY6UpCCJFikCA\nhBIghAChpPdyy+75/bFzZs/MnDM7s3f33i3f9+t1X3f3zJmZMzNn55zv+bYJ9vEB4PcXHILRrfWe\n81fFzcKFaSCRx3SsUMv/mnvbv6nasygh73e+wp0HwU8gdSMHlIZqrzkaAIzp14D1y5dh4sBMEBtb\ns+caXd3P6Jgp/XHL2Qfb34/8/qN4+LVNaZ89n+eRTIWLPFsdizqeoU7Yc7ftsTeyB4Bx49a+uY9p\na/ZyEfYqI/VCV8bIWwEs9dn+FoCFQogpAK4FcGPI4+eE7h3q5yagTtxM74Qg41O2Oo4ceAUSWNT3\nonoO+V5Q74P64LtixqnePzW9ijx/a2M6mMorH+zy7NuZTHlSzRgjgpI8B0t7PYl8xqa1MLmIt7e9\ntHKVqv7u7nHEvQjUYTA1luNRd84tKoVso7EwfC473MKdU/jLlEsevGwhHrl8kaeO33HTkbbSn4MO\nVl0N0OIXVKa1oRo3fXqWcVJcKNRr+t9Tp+PCRaMwpl+DZ5ubtBmnt1wV9nL1EzpiQj985+NTQ+3z\n63Nn45VveOdqdr49l3Bqf5I+ewlnWxuq4/jVObNt/7eZw9JBT1SBUhU03C9UFdN9jGhWf8OYSALA\nnBF9cML0gfjmiZOD7ZAH1Jx42djTljZTCiX4S82e64fs1tYNaq7BlEFNjrJv3P0K2jpSqM4SoCWU\nz55rZTRC5CvsDWiqxsZd/sEddLgnE401znvWFZ+9sMFuSpScx0ghxKMAjOpYIcQTQojt1tenAOht\nn/JM0KTqElXrp0+94O27OpIpgbtefB8X/OY5u9+pbgBqswrlerDRECbeTmtgMOO0IxTnMElVj9mR\nFPZ3WdxiRc50t00IgdFX/w1fv/sVAJl7YtIuyubmsjDM5A95/6OGuYpqbRHEkqVYUAPKuV1UPMKe\noY/KY+giWjNdI9vMeBoR7SKi3QCmWp93EdFuIvIuM5Uwbs2e+iOTHVCt01ybcGg8cgnQEgRz+Piu\nB2jpKdRLGta7DlcsHa9sNO9n0uzlw1zspk8fhE9Y+eOCYnphSV+LREzvsydJZBFMrzthMu6+eD4G\nNGUmPAdZAmCE/CfTpq4m+6NTsxfcRBJIrzz+4NQZgcyQ802QJsp8cHUhhL2Mz176//j+DThl1hBv\n4vtoxGM+2ZlK4UBn0uE76iZ0gBbrWL86ZzYe+vxCbZ49VRBNC3vh8hkJeBdHBrp8bLsSjbNCArR0\n1xh5LoC/5fF4vvzz8wsD101mE/YQbHz608r3cMnvnsffX/7Qju4rrRyAtKm0XPH/cGf43F1naIKa\nuDl19lDMHdnbUx635gDquOtMqm6Zn2k0Fg3VMfzwtBnGc6rv4mRKeO6V/Hr9fWsc5XJR67dPvW3t\nC6sN3nu9ZU8bdlr59bK5bWzadQAf7NzvW4fJHfmcTNYSatyBUtLsqfPbatstJ/3dvQhiWhSxNXts\nxpl3fEdjIURUCNEohGgQQsSsz/K7PpFXiaIGoADckbZkmfmHZ9qki3qoRXn/Hj6+1Xhc+aKWA8QJ\n0weioTqGw8b29W3XYWP023sCnYmsvC71cr967ETHfvGYPgean9ATdGEsF5+kbKu47uihMqR5H8sE\nmIhw9TFmp/vqeBRTBju1SLIPVcej/nnxDCZ0UVuz591WzMOKnJ8EaaMUwnUh5E3IiZC8v3+/dAG+\nfdJU1Madx4gSeYT0ZDLtj+fnQ5lMhfPZk6u7C8e2YFRLPYgIe11moKrZ6IDmmtDCHuCdTAzu5XwP\nqnn2nOXZn0QlBGjpjjGSiBYjLexd4VPnfCJ6loie3bw5vDmvm5EtGbPxbK/GzpTA21v3oiOZ0ptf\naUyQdVz159X25wMdKYxurXcsaKzbvBcnr3gCAHD4DY9kPZ6bIItZjdVxfPbwMZ5yXcJypxmnjITp\nvc5DRvbBSJ8gVqp82JkS9tguj2Qyy9xmBYPqTAnMuu5B7NjXbrXRW3/WdQ9izYe7re3GpgAAZl//\nEA751j817SydBeVixg7QYhgu1PElH1Fsu4uYxoxTtt+r2dN3wiibcRaMilh6DQIRYdG4FgDSZy+z\nLabR7Lkxm82pxyH7R2CaLD30+YX47kkZc0L3ceXrVr53Jw9qwuprjnIE0FCRe3/j+O4zucucW3+N\nulJduPxz5o9w1IlH9EnV/XyDwppB/PxTB2Wt02TlKzINffJa0nn2MudfOLYFyz82BV86OiPgLbT6\nXFDky9PPhBPwmlHY+9vCnvdlWqxRvoDMvQ7zPHPx33zh3R2O75MGOufr0Qh5ooRt2t2Gd7ftR1+f\nRMVJIUL57Ln9/3SXrTrBD8zBjFMIr0A2sMn5HjGZcWbTSqv7MrlDRFMB3AzgeCHEVlM9IcSNQohZ\nQohZLS3h3inZiGfxvdyypw0Lv/svfO2ul40r8mEtT3Yf6EAsQp7FiBc37Ax1HBX1UEN6e6NES3Rv\nmLhG2FMr+uW4q01EfQVm1X/pxXd34OQVTwLILICaoviqkX+37GnDMzKpdZYgMbmawBYyanclIc04\nI0RYvWEnHn3duTijpo8q5jHZjToPk9YvEYPwZnpPxKIs7BWK7otmUQLYCafJKagEyUNmFPYc5qCE\nzyweBQGB01wmJepr1OHnZzinFDhkotRsFCripu85DeJQ2kw2LbDKS9XlPnNjylvnZ8YZjfj7i8wc\n2ozzDsvEVThqUn/UJaIeLYrKdSdMxqsf7MJig6AmzxZ3pV4gIpw62/ncw77K1YAvfpi0TBlhTzVH\nku0L2ZhuJEj/cOMXHdPEgU7nc//aRycBANZ8uBuvfLALUc0EtDMl0JlK4sQZg3CbIeFxqos+e7qr\nrnGYcZonr364r8UtoJnMOKvi/r8R3bGZcBDRUAB/AvApIcTrPdWOeJTg96i3W0LH429swaKx3nci\nIbyf2DPrt6O+KlYw7fCDl2U3U1XHr5htxpnZrs4R5Bik01iYxi3JvrbMzf3Zv970bDdZBGzb40zz\nImW8bEFichf2ctqNcZG0NXuEj/74cQBwBHZzmHFa/eaRyxcVvR+b+lv9xysbsWBMi49mz9uZ3t22\nD29v3aetz3QdXnpVyISldwpZGWEiBx8pV+qF2kQMlx813pj8231uj2bP+o1csHAUrj1hMj5uyFdi\nH68IZ/BEqgCd/m+b6fk0N523zlk2bUgzZg7t5Sh7/IrFOG5aOqx+tsnCqJZ6HDNlgKMsm0DRpy6B\nLy4db/RJkoK1KVWESibHnX89d9uy5RM0af7k/uq7Vl6HXzCfniZI/3CTS993ryg21cbxvVOmo58V\nEc9PgGmuTRi3hc2z59ac6fqk+owHNus1+9lw/z7cE8F4TG/GGSTFSYUEaMkZIvodgCcBjCOiDUR0\nLhFdQEQXWFW+CqAPgJ8S0QtE9GxPtDOe5VlL08XOZEo7IQ0aoMXNnrZOrWZj6Q8eDX0swCkEmcZf\nANphXhegRW2aHaDFkFTe71W0tz2T9063/752b168VEpg+752TxmQ3b3g1ifWY9WGHb51dLBmLz+o\nwp4O9d0v+/+wPnUYawWxK1aiyvt+w/b9OPvWZzIRYt0+e5p+/vn/e9H+zJq9/MPCnkJMWcV2atey\na/ZME0u3Zs/EvNF9AKQ1dTqfNonMP5SIRfCpucMCr54X02taHfxk6+Uqqt/VpAdUZ40/XHCIR8s1\nuFetHVY/m827brMsU++t6mCf7Z7Le21Kqq47l8kM101QM06T5i+TiynzMpX93qSJLQbCzDOWf2wK\njp7cP6fzmASyIGakfj57uw902j422Vi/fJnnfaI7bXU8u2bvlrMO1panEZ6+7J6TSxM+9+9Imumc\nM28Eprr8SiUVknohZ4QQpwkhBggh4kKIwUKIXwghVgghVljbzxNC9BJCTLf+ZvVEO4Oa476/8wA+\n+7vntdtyDSCmW6yTvmc63CmRVIJGypw0sAmxCOGixaPtsrgu9YIm+JrOZ08I//fG3raMMKfbf3+H\nt91JTd5OKYxl04qs27wXZ/7yP751dPS0rJdKCVx7zytYlyXfYLEjFwxMc0bVIqWULOF15t4yiJC7\nT2YT5kzbhRDsO5ojJdSVCo8adVP9GdrmnTkcUx2r/ISEr310Eh65fBH61lcZNXsTBzTi0iVjNXub\n+cKR40AEDMzRzKsQEGWE6YwZZ2abCd2kw6S5k5PTXGze5T4ylP7IvnWYN6pP5pzZNBaKz1425ZIM\nDCJzKQVuWxZhz6R5iUa9kxK5ktjVhMCFJJVlgFQ5dfZQ/OyT2X0vwyAnOn6Dr98zWdvFCYpO6+pO\nvaBjhE9gCMArkLnNvTNJ1d3tSdOrNo7TXabJ9rFZs1cWBPHPzEauwl4Yv8/5o/viF2eaFzeCmoY1\n1cSx9vpjcNiYjEmqX35SIGMdoRMoBfwtN4wLTNapdGacyZTwaFHlLQ6iFfHzLzbR05q9NzfvwS8e\nfwsX/nZlj7ajqySTMrWG/n6qY3MpBWjxm9+6NXlvbdnruX51b5PJ6q1PrMfIq+4LvHDKZGBhTyHj\ns+cM0BKJ6Cc8QXCbcZqIRyN2KHuTNvCixaONgTdMLJs6AG99a1no/fKBX4AW208M0oxT+mTpj3Xd\nCZOte+F6QRgeiizuSroLe/JOTv+obNo62UL12Zle7MOtZ36JJgKciWiEfLVIgPm+2KkXlAFFTvg7\nS8BOvtBuYGcdOlxbLic6UR9tlSp8Nboigb6qSYgcBnWiKiffap/sXafXaGSboLnfSe7a0oRPLh6M\n79+A28+bY69OJ2IR23/Yja+pHFMydNUclyh3Yc+d5sSPSMTfkmL3Aa85ZFCk0OmIz6KcSo4zJn85\nv0WqbL5YOmGvUyfsWecOIuy1NJSesCevN99BS25+bB3mXP9gXo/ph3x3qvdz14EOvPReOvhQKaXR\nUvGba7n75Gd/9zzueOZdc33DXOROa59cUq9UOizsKUQVoU5nohHkN/hRy1dMEtSMU0UdF/xMOkuJ\n7ygRRg90JD2msRmfLP1FuqMiZsPtE2hCJ5DKEjVPW7Uycc1mnqamkcjmB9e7LoH1y5fhiIn9fOup\nRCi7GacJXeoFO7hAEQ8yuQRoyYVrjpvku91v3qtq9tyT1H+9tjlQKgjTgCkjeV7z0Ym2KbfaB9z+\no9KcrbXR3zzYrX1zm8hI05wmyx9xQFM1Dh3d116dTsQixiBRuQTIYYoPP+1aXYBnTHCaXob5CYfR\n7EXI/9jvbNsX/MQutKkXXGM7kVnQ8ntvtRkCsEizep0wmEwJz35SiPjP+m12GgYTubxGNQGcuxUp\nAIRZAAjCdfe+io272nDLv9/Cvas+yOuxdch3rDoGn33LMzj2R48jlRIOzV7xjshe/ITwdo3V0Mq3\ntzu+q33SbMYpzxW+fZUO3zIFOdFy91k5Qc4W0fKVbxyFH5wy3VHmEPYCvmEdERwD7VGcqD5gahCV\nB1/d5AiGk66bxnSLwkb2M5mfudH77KULG2vSk/P2zpTjRRbUZ4+IMr6IeRRSIkSe0Pxh9gWcgQbk\nvQrq09IT5BKgpRDIZ//9U6Y5yuNRZ5ROnenb3JF9PGVA2vyswVpYMEUtlEEc6qvj9mTBbTaqCoqn\nzx6K9cuX+aafEMLblz0+e7H09mZLoLN9lxTNXmO1SbPHw0s54CdwmYJUqRAR/vz8e/b36hAa3zB9\nyJ0yyc2mHPJQSnQLfA7rHyL0qk3YZvmeuj7HPpBlYqub+OrMONX5ySV3vIDhV96LFY+8iZ/+a62m\nPeFfpLlo9tZ8uCunYDA65H1IFMg8/Ot3v4KLbi+8iahOs/ecJfjsbe90LCj0tJ9kWFZ88iCt3+yW\n3W1Y/rc1jjK/BRCTsJfspkXfcoRHYwVbs+d6EerC1euoTcQ8kyd1jMhFs6fuX8rd2+TzIz+ksphx\nhhX2pDCUizOvPFVDVXoS686PllXYk4KJoyx/b+1ohEJNmFSOn57WPB+naKBjPgmBiwU54ezlE/Gy\nkGRWFNNPdcIAp6bZ/Tx0k2DV71Pls4ePxu/On+s4jxvZB+sSUVsgcwt7YZLI2+3MFo3Tvu8uYc8a\nixNRsxlnEEGAKX78onEGMfF0v/uymaCrBNHi2Dm9yDwJPHpyf9xytl+wIn9016mWRAi+fkR+k9O3\nNu/1Pbdu4tuZ8kY+VV/f72xNH/Pmx9bhO39/zbM/EYyCqYlchL2lP3gMx/3436H309FRIM1edyPn\nkbo5we4DnS7NXvGOyTqWTu5vuyOpfOOeV7DiEWdaEV9NYKde251tnsiYKe1fTZ6Rk3h3J7T9nHL4\n3UVdph5BUAcGZ5628OfvSXS5CgHg9vPmKNE4nffWpAELqyWQE+FsA5TudLKtcvK81xX6Ottz/NFp\nM3DkxH4Y3KumIOkMokSoSeT20x3ZUo/1y5dhdGu9XRb3yRFVLEwY0ICvHzcJ33dpzrsLOehK4cgt\nJLk1rbouMmVws8eXD0ibiGULmCAndtXxKDoteyr3byIX8x+Pz55b2HOZccp3Y0rV7Fka8EQ0gmtP\nmIzzF6TzVnLqhfLAT5MSJOKqexEpjAl6kOAw8nhuX3uVq46ZgIOG9Q58Xjc67aYuYrfkhpOn4eef\nOghHTGjFZR8Z6zt2f/9BfQpFgfTvUee/lEp5zTvVxegGS9tu0pA89sYWTPv6A3j+ne3a7Tp6ei1Q\nXksY096e5IOd+7Fhu9d0OKkx45Tsaet0RMouNc0eEDxOgt9vwuSzZy+6ltpkuAgojV9NNyEnMO5J\nVNQ1wQmDurqdi7DnFH6Ks4N/56SpuO6Eyb511EufO7KPJ+jNFUvHoV9jFcYZcslkXvDB7oF8hrk4\nO8u21lsT8wOu0NfZXmbThjTjxk/PQiwaKcjK3Kmzh2DJhOA+ftmIlUA0TiLCmYcONwYiKTRy/JW/\nTfdg40mCrhmM+tQl8ODnF+LCRaMc5dFIJOt1qavaVy+bCAAeE82wASgEzGaccpItzTgb7d9CesVV\nBvNJRCP2xHLOyN741Nxh9r6ceqE88DfjzP4+di8ihTFBz5bjD8gIe36ava4G9dBdJxm/AEdO6oej\nJvXHzWcejIHNNTmdX4iMqebSSc5UMmnNnstnzyHsxax6/u/0l94PHjgqn9YpudChvHMkH+48gF0H\nwmkou4tDvvVPzP/2w55y+Zx088ndBzod43Dxjshm3JrXOSP0iyzun4S6MK5Lug6Ei8rNOOHRWEE1\nTVLJ+OyFP6a6uh3YjFM9t/KlWPv3J2YNwSfnDvOt49BWaqKmHTqqL56+6ghHUBQVOeGQj+b/LRmD\nRy9fbDxfRrOXtfkebJ89gy9SWJPSfHP1solYPK41a71/X3k4nvrSkqz1MpHkilez19NIoV0+e3cf\ncL8zdD2kV20CrQ3VmDGk2VEeJcpqmiQnbTK/5vrlyzwT0N+eO8ee+AZ5VxG80UXdC17ydycn6NLv\n0BYKYxFEI4R/fG4BVljpLmQ/koIiU9rohD05rgVZxXdP3MJYaQTR7MlAQBEy21F0tSfqcoi5LXBW\nfHKm/d3jzpFDAwQyGo5Zw3vhEMXnV+ez16m8v2sTlrCXZQEvTLPcY+mHOw9gxSNvOoTAzmQKC77z\nMO5bnf9AJ/J61QWAud96CEtueCTv5yok8h2rE/bSmr3S9dkDvAsuJvcCP+2caZFe3rOeXngoRVjY\nU5Dqc/fEy2TeGQT1WEEDtLgdv+3y0GcvHtyXLm9l0N+snHDI+9GrNo6hfWqN9asDa/a8d1W+SHQm\nd0A4Ya8QZpxBGdRcg/6G/Gsq8RLQ7BWaIya04rz5I7LWs98Frg7tWUHXPHY56Ln7j/p98iB91FnZ\njxM+lgLzx/TFZ1xaQxPLpg7AdSdO9kzW5VVUxZ3C3kkzB2PFJ2fibCs1hUyWK7V6Y/o12As1sh/p\nJshM6fBfh43AyL51Wp9M+d4O4pfpfgeH8bkKYgpcFVeEPeV3uX75MvtzV82+dGO/c5wGxipWKe7z\n6caBH58+w/ec6zbvxdRrHgCQvmeqcNCZEmhzWZyoQrVMy5Itt2CY2+IOHnXhbc9h+d/W4E3F53DX\ngU68s20frvrz6uAHDoi8vipXn9u8uy3v5yokcqzQxRP44UNvFLXvfBBk6xPKnO2ixd5xyfMbUb4a\nhT2rO5f4LeoReDRWsDV7BmEvFxeUqmjGZCV4gJaM74HZpLO08A62Ulsa7FdrO+EHFLzlBCCbsKe7\npVv2pB3tR7Q4HY0HNdcEOmapEbMTrVeuZu/mMw/Gl4+daNyeSaqu1+y5+4Sud0YM+8r7/8o3jsKf\nLpzn20713dSVCewNJ0/DgKYa4ztJ5siTk+2aRBRLJw/w/O50ydyluRUnVS9trl42Ef/8wiJtQBXZ\n23PJY6pq9nT97+6L59ufg4x5NfHM2GCqni9jjD4Gc+sIEfz863XnD5OsPhF1CnuplPCYcarRlKsL\nEMTELZzsa7NMupVxQ74HC+FTZUfjDHht/3ptE+5+8f28t6OryMUw3drqc29vd40lpTfXkHO6Jiuo\nV2dK4NOHDPfU8+sirNnLP+HDt4WAiJYC+F8AUQA3CyGWu7Y3AfgtgKFWW/5HCHFLIdvkh63ZM5hx\n5mK+p76YwuUMIiSFcObcC3324sH98rc1ewH3T7g0e9kmAXKCkjVAi6ZMroZOGeQ0t/v1ubPxi8ff\nwkBL6AuC9MU6dFTfwPt0N9K3ymQnzyjCHunN19ymKn6THZNmT5pe+ZEvYS9u+9U5j/GFI8chlRI4\nedZg3P3iB1mDafTT5PGT/YijcZYHbq2U2u10Y9rIljqs84kwqfbhlvoqfOhKiTCgObs1gorUYsWj\n5qTq+VgovfXsgx3aO3c+XPV37dVaeM8fZj6R1uxlvuuSqqtavEIIW+6hVC78qBYhcqGnEK4OYaNx\nnnXLMwC8uY9VwggNBzqSuOuF93HizEGh5nIdyZSjvrwO07k7S9yMU43gvHl3G9o6k57I0UAWM07D\nhWdMYPPQ0AqjYKMxEUUB/ATA0QAmAjiNiNxL5xcBeEUIMQ3AIgA3EFHPRGBAZqXFvSItJ3hdFfbC\nRCGTZ3IOKKFPXzS4b538oQcNeiM1DPLRZHsUUjNhemkEYURfp2ZvVEs9rj9xSqh+0L+pGo9cvghX\nHTM+53YUmkRMDtqVq9nLhvTZM2mWpVB/y9kH48vLJjh+q9MGN2Hm0MzCgUezF2bS5zDjDLybB9MC\nVu+6BJZ/fCoOGtYb1xw3KeskWTfxygRvKeEXFmPj7gJRInuiqtPe/vBUf/NEtQ+3Nnqj0IY1/5Vp\nT+KRiHFBNB+yx6JxrY6FPqfPnvM+uV02dOcP4xZSFYs6xkqdz15bZ8o2uQ3qfy0F+QMdSezc3+ER\nQNTv7rFaPka1XLYpqMtKGKRmr60j5UmHlCthhIb7X/4QX/zjKnz3fm8qCz/cgbP8ArQAwKsfZILm\nlKKwJ/ugjHnQ1pHSzn09AVqU76aUWbI4l2CJ+WJPWye+e/+aoo5erqOQS6+zAawVQqwTQrQDuAPA\n8a46AkADpWcU9QC2AQgXUi6PpAwmCHJAy0XYU/epCpFfSBf1r5SFPTlplEKb1IQE/c3GXJq9bCuX\nUrOX7fi6w/z90sNw06dn5W11clifuqLWcrBmLztuzZ57MtPf0nAtHteK8w4b6dCG/PxTs/Cnz2TM\nM91RKsOswgfV7PlFgVWD9uRihgcAn1k0CoeP1wcJumDRKIzr14AjJ/bXbmdKC69JYqZAJ5iZNC/y\nfapu1wVrCWv+K48X89HsFULTpR4x4tLs+d0zSRiBKBGLOCbAyZTAvjbvVEnmzgvqf33PqrSZ42k3\nPYVpX3/AYz6nfvUIe1b71X2kQFaIGGZSc3nns+9i/rf/mZdjhnFdkBqrh17dGOocu5VooZ3JlB09\nNIisUGp59oCMsCf7S1tnSp+n0qf/m/wWU1kE5e7gew+8jp88/Cb+tHJDj7UhFwppxjkIwLvK9w0A\n5rjq/BjAXQDeB9AA4BQhRI+Jy/KlZZrkuyPXhSVMFDL5O3AMICVsyJkJepJ+EfzyrIM06ol3AAAg\nAElEQVRxz6oPtD4/ftialYCavVwY378R4/vrA2WUI8OtJKgLx7b0cEuKH12wpk/NHYZLjxhj3Mfd\nV8PmjFTxC9ACBFsQUoP25PpO++JSs6Z6VEs97v/cgpyOyxQf7klZJALAGqV1fdBk4lYdi2Bve9Kx\njy49R1hhT40ca/bZK4Cw59Ls+fnX607vvndE5sVJnRnnHo2wJ9sSVOvwxJtbsa+9E8+/s8M+rooq\nDLnn3+6cmwBsP8JoAbT6as7ArT4J7LOxbvMefOWvL+GmT88yzql27u9AIhqxTYSBzD013XcTu/Zn\n6l90+0ps3JUOKBNEYCllzZ5cPG7vTGkFO3eJ+iyy++zloaE50p5M93FTDstipafVDUcBeAHAQADT\nAfyYiDyzbCI6n4ieJaJnN2/eXLDGdBqEPVsI7OL7K4zfgNY3rdRkPaW9Uvg6YcYgAMCQ3rW4cNGo\n0L4UQVcMw0R88+O8+SNwlhWBsFwZ2qcWz335CJx32IiebkrRIscWnfnjtSdMRh9XUnS3P4+KOoEI\nSzbNXo3l9xfUZDzmo41gGMA77ESJMgFaNIOiKYJmtRI1U6Ld37AAcZzB9yphC3vmAC1UgJmOeioi\ns1ZRbldZMLbFU79fg3nh0x2gJZkS2GPIqxmNUM5WGntdgoyq+HKbeErNnnquQppxus1Wc+X6+17F\nv9duxcSv3o9HXt+krTPt6w/giO894iiT1xk2n6kUDgDg/pczWkGTqaJKCcp6GN43HSW9b33atcEd\nSEjiF+jOLOyl//eksBfRaLRLgUJq9t4DMET5PtgqUzkbwHKRfousJaK3AIwH8B+1khDiRgA3AsCs\nWbMKdofly9T9osoIgd0nG8s5WKSEZT31TdVUG8czVx/R5aTYmR+af73aeBRj+9XjzEOH4+o/v2Ss\nl01b6hehsZxwCyuME/vdYEfm9e836uTOvXjkFsTCzIscwp5m5eNTc4dhf3snzjtspGfbhYtGYdJA\n51qaugJfiAkaU/roTBLlZEunxTNFmZQaOIcZqKauu1//43MLUFcVw/X3vao9bsaMM2JcPJTnbGmo\nyluofneKJL+FSHXb01ctQVNNHKs27HTUGdBc7QlWI6mKRxwT3M5kCnvaTZo9CmWeqB73zFscUy+H\nz7vb/10+J1WLKLUde9qSuP3pdwK3wc3wK+/Ffy8ciS8dPcEuMwkNYVH731U+c4P3dux3fJfXuS+k\nv2B7p37K+qLr+VfHIzjgSqdRipq9E6YPghBp95WHX9tsFNL9UkzIvnbRbSsxc1gvnGulRfLLUdhd\n2HPQEns2hZRengEwhohGWEFXTkXaZFPlHQBLAICI+gEYB2BdAdvkiy3UWROgCxeNwi/PmmWvwOTq\n35ILctAql9QLQHqg7aofnM50xFTvgc8txLFTzZG4ANZmMMGQ3U3+HrP1Y6c/j3ObOzKZLlKZCfUd\npBPOErEILj58jOOYLQ1pQf6KpeM9vwfH8brx/caUDu4FsYjDDDOEGafVJ6Mh+9yYfg0Y2FxjNnGM\nymicEaPAJYsfvGwh/n3l4VnPGQS3GafftajjeL/GalTHo562DmwyR3l2a/Z2Hegw3o9ohELlTJ30\ntfvtzy+9t8uxLakc54MdB/DEm1sy57GuSTVnkwLZlj1tOefakxrEnz/inArmS7OnPoswQcnUutmi\neKpaO5NJ7YMu3z9dNOZS9NkjInxs5mDbnNMs7Jnvvex3967+ANfe84pdXgw+e/J3XmrpHwom7Akh\nOgFcDOB+AK8C+L0Q4mUiuoCILrCqXQvgUCJaDeAhAFcIIbboj1h4Ura5ZvphXrF0PA4f388WAnNJ\nqp4rmTx7Slm3nb14kfcj6I89zKScYUz4mXHqcKz6u+qqect+919z0TeEVtUxwQz49n7si4vx6jeW\nareFnXgzlYdfzjidYBc3mNAntJq99Of6KudEt39jNa5w+YXu79BrVGQT4pHsAVqaauJ2vtR8YI/T\nEfJdjNW1y/1ekIsyOqpcSdX9tJMRorxFClS1eWff+gxOv+lp+7t8XziEvY6un9ckHOTj2IDzvZnN\nFG/T7oymtV0RfPdm0e6paTCCPgvV4kOXhLzUqLXcFdoMv1v3goT6EzGnXnD+7wnkzzaoGed7O/bj\nK395qcejnRfULlEIcZ8QYqwQYpQQ4ptW2QohxArr8/tCiCOFEFOEEJOFEL8tZHuy0Vhj+bu4fGoy\nJlzd1xadzx5roTKCeBB7d8A5MZkwoHKCrjB5xmXGmU0ucpt4qahat0NG9cm5SUHNLqvjUaOfYDZN\nIcN4fPZUzV4In71MDjb12Om6R0xwRnZ96qoluHCRc8KrRjVUkUJTNEq2CalbeCxU15aHpSyaPW2A\nFgou7LkDtJjMPYFwAVqy4TehlfddFWza83BeU1qFfJlxqnMqP1NCADhDEW7Ve7p5dxu+/JfVOGAQ\nZHIR9tRFwBF969MfSkt55KDO0lSa+sQdz7yLB1/RRzYt5qTqGeuyYPUvu/MF/Oapt/Hc29sL2Krs\n9HSAlqLi+hOn4KvHTsSsYb0c5UmXxq870PvsldhkrADNlS/qoFYq6kRb94IoddNYpntwa/ay9Rv1\nt+qeA4ZJyOtHPiwNVD/k7rRcYEoHTzROH39UwBxgpd6KxKzuLydvMrCQn8CkRjWUqBEVhQB61SXw\nizNn4Xf/NdfY5nyiulv4/XyCROPs4+PPXhWLOhY4X/twj7bevNF90maceVJ96CxoXnov7WsmZfp8\na/b2Wr6I7qjF+TLjVOdx2bQz72zbZ39WNTPL//YqfvvUO7j7xfe1+6n3pD3gZEUuyEUjmdGjhGU9\n+3pkYBu5sPidk6badc779bPafZMp4ZivddqJ6NPfe1azF8yVSCL7rcniobtgYU+huTaBc+aP8Axu\nyS4GaGmsjoU2HdH57DHh7aXV21diJtZMESH7jvvX2GpYjffT7OWLfPgQs88ekw1391X7iU6w0y0a\n9K5L2MGB1O1ybJUmbH59WuaQU6lRfN/kmLBkQj9MGdzkbFOBura6KOv3O9eacbrKGixhWEfCZcb5\n77VbECFgyiDndf709IMQjVDOYeHd2kWd0Hjsjx7HE29usfuBIxpnHjR7MgCK25fZLezlqt1R+0I2\nobhWsYhQhTbZRtMzV+9/R8BnURu3FjwoE1W21PzCVKTZ9rKpAwBkrABGt9ZjVEudXe+JtVvQ1pl0\naHSTKeEQxGW6i2II0KLLL+mH1OyaFsG6i0JG4ywbpKA2UumgYXj+q0eG3icziJSwGWcefo//+NwC\nx0s+rL20bhVZZVSOz5SpLKSjvPobvO28ORjTWq+tP3FAox1tLx/C3vj+DVjz4W5HmXSA7wpywlZf\nFcP3PjGty8djyg+pZ0hEI2hPpnxTJ5hM5Q8e3sueJEUd7+T0fzmp9hP2dmnMONXAKH5DQsE0eyAA\nAgT/1AvapOqua62tMgdqcptx7u9Ion9jNRqqnVO46kTEisbZtdQL0hzW5C7x9tZ9is9eZpJu8s8K\nQ0bYi+DJN7di7sjeICLPsdWmJVMi8GJVNgsG1exSDZqiavbknMSUE9Kp2Qtoxqlo9mSqqjDBu4qR\n57/yEdRZJtXxSAQHkEI8EnHk1zz95qcxvE8t1m/NaFGTKeG4b/J+F4OwFzRIoET2p7D5Q/MNa/YC\nsHRyf9xx/lycMWdoTvtHIxR61VwODuqLqdRkvXwwpl8DJiurlxkVerD9HWacrm1/uOAQfHLusK42\nkakAMpq9TH+aN7ovWhv1ubGuOW6S/TkfWoW7Lp6PV75xlKOsubZraUyAzOT6LxfNw6JxrVlqM5WI\nfIVW2SkOFM2eyyT5tNlDYELn7ypsM04rB5/Pj+WWsw7G7BG9XY3LjIt+k6+CLZQqmj3/aJzeMrc1\n94whzcZjJKIRzwJnPOYVMBPRiBWNMzcNmxS05L00LaqqwlWHRuPVFfZZAufGXW047aan8MeV6Yxd\nbs3eLx7PROv002S++O4Ox/dsgr8aCEiXVF09n1Gzl1OAlnSHaKiO4ahJ/XDJ4aPxpWMmZNmruOlV\nl1BSo5D93y34qIIekNa46syDiyPPXvr/Dx58Azv36f2IVXLNeZlvWNgLABFh7sg+3erfpQ0EUWrS\nXgHaG3ZVRb1/7n1mDe/NPntMIOQgE7S7qCuy+TCPTMQintDcjdVdN8yQaWa6M60MU1rIniEnbWoe\nPXe/Mb1PhXBGrpTId3J1ADPOOSP74KdnzHSURSgTBdNvAbBQ73l5VMrqs5fdjLO5NoE3rz9Gu388\nSh6Tvngk4nkfyeTuXZ1gyudi0hAKIez2q4KNOyl7Lvzr9c2O78+8tQ2AV9i7/r419ud1W/Q+jABw\n/E/+jY1WQJtkSuAPz23wPf8BRWBVFzPUeyoFEVN/dWj2lM9+3VAef1RLPWLRCC47clxerDeKhZh1\nL+NRsj+bSAmnZu9Ap1urq++Xf1q5AX9a6f98u4r6u/XrdxIp7H/rb2t6NrBMj52Z8cXO56WacZac\ntJd/7NQLAVV7jkFW2eWxLy7OY6uYSiHMnPHa4yehLhHVTvRq4lGM7Ns1E+JsA2agYwRMJcFULm7N\nXkIJNODug6ZuJABfM041Mbofbr8XQmas7NGJFPmnXtARJlATEeEMywqlv2VNEIvqTUcjka5H48yE\nuDdr9uS5H3tjM55ZnxbI9uRB2LvxUWd+vQ070lofv2icT6zd6ntMmarinlX6gCoqqnZSvbsdDrPC\ndB3Te1MVTFUhcVy/BgDO35BkaO9aAChba6N4RC4sRrIuLpo0exLT9O+y37+Iy37/oqNs94EO/OTh\ntYFdf7LhTj2zdtNuPLHWnDFO9ptHX9+Mt10azO6Ehb0iRY5pDjNOno9lUi/k8Lsd1Cvte/n9U6Zh\niPViZZgg5DKR/NQhw/GyIb/dS18/Cg9etrCrzeoyMugUC3uMCSnEVMUzycsl7kmbyaxNCCVFgkaz\nJ9/r2SaBbvMvgfC5V/OJPGMuP5+wPjyfWTQK664/xg6iYkoiH+2Cz57K21v3GqNrJkUmF9pT67bh\n5BVPAui6Gaf6npV9YeuedgD+kT7vfPZd3+NKITRItFDVjFPOGQCXGaedRiS7Zs9t/nns1AGeSKNA\n2s/75a8fZQc0KTekFQlRdqOvnfvacdFtK+3ve9s7Hfc0zG/92ntewXfvfw0Pr9kUqr0m1Dl5VSyC\nI773KE6/+WljfVXY3+iTMqXQcICWIsX22VN+FaU2HVswpiXvx8ykXgg/mDVUx7B++bJ8N4mpIPKl\nXS8W4SoWCTbJZioX24wzqtPsuYU901EymiAir7A3qFcNjp8+EOfMG+HbFp2AFDbvVV4RzjZk49z5\nmeuLhYzOR1aUxjorkEvM8s9zE4lQ4AiQfiz87r8wb7Q+D2gqJZDUmIp2VbOnRlyVYzxR2oR1X7v+\n2A1VMazdtAcX3bYSN3ximjaoyQ7LtyrIc3p9YyYQlpozslNjxnnur57FT8+YiWOmOAU0k89eezKF\nRMzrfynbVldVvlNyqZXvSIqsiosXreBmklNvfMrxXbf4aup77+3YD0CvTc2FsMGe1Oe/ydIw9wSs\n2StSorawp2r2SmdCtn75MvRv0gev6AphUy+olHAUY6aHEfakrmfb4WbFJ2fiF2fOynn/CQMaMWlg\nIxrLyDeEyS92GiDr3atqJdymiH6LIbKqukvKmgfFI4T/PXUGpg1p9m2L24xTCKDOCqJR2wORC3VR\nek2sX74MXzl2ov3dlHw+G/V2dEOn6aiMzBklQkeq68IeADy7Xp8IOimEdsE1V5+9ZErg4dc22RNz\nwJnuZu2mPdhuCIYxwoqofe/qD/DQq3rtzd0vvo+te9oCaXf+3x0v2J8dfnqaAC0A8JnbVnquW93+\n+Not2L633TpeColoRKt5zVP61aLli0vHIR4lDGyu7vJiZyqVjs4r7ysAfKD0HRUp6Nf5RLsFgGfX\nb7NzSAYliAZdFfZ6UrNX5t2rdJEDKzmEvZ5qTXDOnje8oMcPm3oBgO3Uz8IekytSA1FsfrNLJw/A\nkgn9ct5/9ojeuPeSw0o+xDdTONw5v/zMON1jlEznIYTTD10mEE8p2psguDUzAgJnzB2Gy48ah/9a\nMDLQMQpBLqkdTP6JvzpnNhaMNVvFSO2PquWcPby3bRaejsYZfrDTaT5MAULcedAkuQp7tz6xHmff\n8gzufMZrjhmJwJN2RmV4n4zv8wFD6od7V3+As255Bveu/iBUuzpNGjqX5nTS1+53fFe3P//ODnz6\nl/+xjiEQ10RWBco/p/LSyQPwxjePQW0iFlqr7SYlBGZ84x+Yce0/7DKTZk8KhNmmjCeteBLH/ujx\nrOdWn12Qeai6YLB9X0Y43bGvPed8mLnAwl6RYjuzK0+oFF4FX/vopIKaSh5tmUscN31g4H3kfRP5\nSPzHVCSLx6fTEvSqYw0YU1m4F8kc0Tg9AVqco5RtiYGMQEdEeOJLh2PNtUvtY3dlpT8ejeCixaN7\nZMHC1vjn4rNn2Gnh2BYM72P2KbfzlkUjdiLqj0zsh35W4BbK0WevWSPYNdfq33cdyZRe2MvRZ+/d\nbenAFTrNSoTId1Kt5hp0R21UWe2jtTHlh1TvoypAuyODumlPOtvx8vvpc3d0pozCXrGY9ncHfv6q\nC30WOiQCXkHL3ec3bN+HNzfvwfs709q0XBZAdKRCCntqHVXwO3nFk7jkd8/npU1BKF8D4RJHjpnR\nEtPsFZpRLfWhhcnMynQBGsRUBFcsHY+zDh2O1ob8myYzTDGTtEwC5fsz4TDjdPnsuZaPpRZQNbsn\ngp00WpoCFmqee/6CkZ7ojvkkE6AlB81egIs+f8FInHTQYEdZvSLsbbW0FqoGLldzwObauMenyORD\n9oMH3/CUCSFs4TMMb2zcjbteTEfJ3Lgrff6Gqhh2W5oagr+5nKppPhAgAIuOrxw7Af98dRNufvwt\nR7kuKAuytAfwav6kINdm+ezpiFbQBM+v7wf5XejceNRndaAjifnffhj9GqvssnwFcFIfvWrKLITI\naqGg9otte9tx8Iiu58oNCmv2ipSoxoyzNHR7xUhmdZlhciEaIQxsrslekWHKDLkaLS0j1Mm1Wxth\n0uyZGGZpsBqqC6Mxv+qYCQW1NBEhzVBV/O6NnEMOaq7BWCtcv6TOyrcZoUwCclW75Sc0nDhjED42\nc5B2m85kM4yZWXsyFTrlw0W3rcRHvv8otllCq/TZa1YsKATM5pmAc/FBl54hyKNJRCP42MzBnnJV\nGxRUWHh32z78+OG1rjakg8ykffb0DQoa5KccMJlxnjB9YKD7oJO11Wcl+4tcPACcAvo/12zEE2+a\n0yX4oQp4qpavPUDfl3VSKYHt+9rRu5aFvYrH9m/g1AtdhjV7DMMUG0T0SyLaREQvGbYTEf2QiNYS\n0SoimqmrV2jkSrsMq+/Q7Lnz3rkGqZhixql7AV93wmTc/OlZGNe/wbMtED38Tu9K6oUgAqKuigw0\n0daZss0WaxIZE1bdZPm4aWm3BwKMCdebarwTzzCaup37O3yFvV89sR5/ef49R5nOhy4aIdRXZYS9\nVRt24st/0f5EADg1QTrN3pjWeoxs8c9pmj6nU4uZiEXQEdJkDwC+fvfLeHebM1hI1DJFFcKcX7GS\nNHtRjcD790sPww9OnRFIs6cTvDtTqn+ld3tS2X7Orc/i9JvM6RJ8z62a9rpMNIdfeS8uvn2lbjcA\nmcWTXQc6kBJArzoW9ioeXU6iynkV5JfMfWNpj2GYouFWAPpEjGmOBjDG+jsfwM+6oU0eBlga7cFW\nzrGPKnnATKkXjpnSH589fLRD0yfN8qRmCgBqEzEcMTH3AEPF8kbPd3ANP/9yKZS0dabsvHG1CX/N\nXszOcUbYtlcf/l3nn7ffR6Pm5st/fslXE/i1u17GpXe+gFRKYOU7+iifQFrDGMYU1WnG6W1vPBrB\nkF7+eXWjEUKtEq1x1rBeOGx0X3QmU/bkPpUKJtQ3a7Q10QjZAkjcZMZZUZo9TR/VzHlN3PDA654y\nU+RUSd589kyaPavv37PKHARI1pHm131Y2GPkGKq+t0sp9UIxIe8ba/YYhikWhBCPAtjmU+V4AL8W\naZ4C0ExE3Z5xeZAl7PWpq8L65ctw6Oi+9jZTgJafnnEQPn/kOHulTYhMjqlWxY+m1OlKgJYgx9Ud\ntlYR9qRwU6MEp9EJntm0XwA8mq103eDC3ltb9hq1hsOvvNf+/IfnNuBjP33Co+WT1FVFQwnPqt9o\nm6a9sWhEm8RcJULkWIQY2rsWsSjh5fd3YcJX/471W/YiKUSgXG196jMT+FpL4xqhjACiBjhSNY6V\nbsYZtcpUYe+VbxyFw8b09dR9a8tex/e1m3Zjl5KjUbfoIDWzYSK561DNODs1wp7KAy9/6PguNd8y\nQihr9hi7wzsCtPRUY0qcTDROhmGYkmEQADUW/QarrFs5fHwrTps9BFcePd6zLe5OveDarr57N1k5\npvIZ5KhYFvAKthCrOW69NOPsSOKANcF0mnH6H1L6tZ0xZ6ijPBohXH7UOHxsRqaLhdHsDWyuCeS3\n9MU/rgIAXHrnC9rtVbFoqPuZLUBLPEJZI7XGooTqeMQW2uPRiC2QtHWm8Nzb25ESwmO27ObpdVvx\nTyXXnzxvNEK2MKBq9r570jR7MaWSzDiDavaqY1HcYKVvMdHWmcQR33sUl/9hlV2mE7ykYKamPwCA\nB1/ZiO//I6MpFELg+vtexRsb9ek+HNE4hVnYe2/Hfpz/m+ccZX976UNs3t1m5/7TRcAtFCzsFSly\nZcuZVL2nWlPajB+Q9geRfgsMwzDlBBGdT0TPEtGzmzdvzuuxE7EIvvWxqRjS22sK59bsmSbpQghM\nG5xOmJ7Nf6oU6erY3GCIeqlDaqDakyl7gqkKe0+tyyiLde2Spp9uEzIhgIsWj8aMYb3ssjDRLTdb\nmludhjAMVbFIKE2p2gd1AVqiEcqq2YtSOkG9vLeJWMRhovzaxt1IpYTRBBMAlv9tDU658Sm8sWmP\nXSaD3nSmhK3VSUTJFmiqYhH7cxdTz5UUSyf395S5FRzxKCESIbQ2VOPqYyYYj/X+Dm+icp3vqNTo\nbd3jFPbO+/Wz+N+HMtFld+zrwI2PrsMnfv4kvnnvK7jiD6vw95cyppmqYlAV/Pa2O/P8mZK8X3T7\nSntRpCrefQ+9grpXaZER9jJlxZbQuVQY3KsW65cvwwkzun1RnGEYJlfeAzBE+T7YKvMghLhRCDFL\nCDGrpSV7nqp84fbZcwsXqvB3xdHj8cjli+x8cOVEV3327v/cAsd3P4WlTIfQpghitQbN1eiWegCZ\nVBcCwg7q4k6rIP0Ec7UmlGa6dVVdy3eYFvaCN0Jtry4lQiIWyTqpliaUag5D1dTwtQ93IymEJ9WI\nyopH3vSU3XzmLABA77qELYDEoxHUWcJ5QhFsK0mzt3h8K5ZYuWslUrMn3ymquavf89uwfZ+nTJcH\n0aTZcyN98rbv68BNj72FO599Fxf8diXe2Zo+T9IQoGWnYkYKmKNzbt/bbveFriaXDwMLe0VKZrWn\ncl4ADMMwjM1dAD5tReWcC2CnEMLs/d8DuM3a/Mwq49EIhvXJr1bPL5BJd9LVYTpMWhc1QItE1eyp\n/PSMmbjxUwc5ji+FRI+wZ/sfBruYg4b1cmgHt+yxhL1EVzV70VD3U60qJ+JqHrZYhGxh14QUNGSQ\nlniMHILd6xt3I5kyR9J0UxWL4FfnzMaolnqcevAQHOhIOoQ9+QxjEdJGXq8E3FYA9pzXKlf9I/2E\nInfkU0Bvxim1cI+94W/5YPLp67CieZoCtOxyCXumgDDRCNnbgkQezRcs7BUp8nfAZpwMwzDlBxH9\nDsCTAMYR0QYiOpeILiCiC6wq9wFYB2AtgJsAfKaHmmrEO0F1TnBsn70CyGSHjuqDuSP75P/AOZBv\nnz2/AC22Zk8xWTSZKTbVxnHkpP4Z4UlkhEQ1Nx+QmeQGnX+eeehwHO7SzgBAfXXXhD2icPdTXRDf\nsa8DW/a0OUztopHsAVqkMGibcUadZpwf7DyAXfs7AgVoAdLC+8KxaQ17dTyK/e1JtHemG5WIRexn\n2JEUdvsrbWHfnT5BCnS1itZT4qdRdWvUAEM0zpTA8+9sx08e9mpg3fV0yOYmDT57uw84zThNaUji\n0QgeWrMRgNcyopB07VfJFIyoxlmVYRiGKQ+EEKdl2S4AXNRNzcmJRMw5PrnnSXaO0wJo4L5z0tTA\nmpZCk+/UC/KydMeVWiF1rmwSjnS+/1JIVDVwTTVx/NdhI32PBQBHTOiHNzbtxttb96GhOmb7y0Uo\n8+y7qtnrSKZCzXvUmk+u24pZ1z2IN755tF0Wj/pr9ka11KGfFSFWRjWNRyNIRJ37bN7dhhF9g2mm\n97ZlJv41iSgOdKQc0TiXf3wqrr3nFQzrU2ubb1aSGSfgFfZk7j2ZY1G9H36/c51QpY/GmdL693nr\n6d9VbZ1JvL9jv8MnU9XeuYVOk7C3+r2dWP3eTgDBNcX5gIW9IoUDtDAMwzDFjMkUUFJIP/Pu9HfJ\nRr7XZL9w5DgIAXxsptfPvDoewYWLRmHZlAHYtb8DT7/lzN7xxwsPxcd/9oTVLplfL7NdavbUZ/fg\nZQvR0lDl2EdHNAJ0WPs3VsdsjUtVLGpH7szFZ2/yoEbMHdEHNz/+FjqSwVIc2Gjau/QHj9qfY1F/\nn71JA5tsAVdqWuLRiK1h6t9YjQ93HUB7MoV4LNiDlnnUgLQA2Z5M2Qnq49EIDhrWC3+5aB4AfU7l\nSsAtU0mTRqkZblcEKT8NmC4ojykaZ5BUIibNXltnCocu/6ejTD3ergNun73sC1zdacbJwl6RwgFa\nGIZhmGKmyRU63L1aLzUm04c05/3cRSTr5V2z11ybwDdPnKLdRkS4YmkmDYaa9xBwmrxFXYvGAhmB\nXAplA5uqbUEvXdfcrmiE7ElsQ3Xc1kxUxSO2sBdKUAPwiVmD8Z2TpuHvL31oCSvR6OMAAB/+SURB\nVHspRCj41FTX3Dc3Z/KwxaPka0asCqfyeuJRsv0ge9Ul8KGVNiSoJkbVDklt4W5LGHCbJMpDun87\n5Y5wa/asjtdgm7hmBDa/1Hj72zWavaRXqEumhMPP1fST/Z8HXtOW6wTIK/+02v4sNXsysExngDQk\n7mjGhYSFvSJF9gF1tYc1ewzDMEyx0FjtL+yN6deA+y9dgNGt9Xk/dzFp9oppbHbMGWxz0Mz2W84+\nGH94boNRk+QnuG7e3WZPwhuqY7bwo+4R1n9Rnk/6EHYkU9i21z9iYtD2AukgQvtdYfFN+0tBVU2J\n0FSTmSYncpicV1tC4y7Lp8stDMvzdzXZd6nhMeMkp2avWtHGpnzujS4XZBDNnkmrdu8qfQwsVVCM\nRsjzvN62onV2pFJ4c/Meoxmnip8vYr5hYa9IkS8aYjNOhmEYpoi46+J5ePT1zbapm0SnnBjXv6Eg\nbSgms7dcNXuzh/fO+7iu3hdbs2eVCSGweFwrFo9rxZoPd2n392vPxl0ZYa++KmYLLupjD3sv5BxH\nDVrSmQqezD3b6aJRwp428/HU+yX9IZMpYfdtdUEjiNayqSaO286bY3+Xmj0ZrdGtHZT3q8JkPaRc\nspDso/IZqL6ffoKwzjRTlx8ymRLY3ZER+sP20zblPPGoV9hb80H69yQEsOSGR3DtCZOzHrM7F6xY\n2CtS2IyTYRiGKUamDm7G1MFe08zuTIVQjMLeT8+Y6RGA/fj9BYfkvS2qxiLjs+e9V6ZtuhxlkkQs\nglvPno3bn347LexFvVqpsE9Fyj5ykt/embI1Hg1VMexuM2vlghCPEPYF1OxdftQ4bNnThtkj+uD1\njbsBAFVKDsMgPlYLxrZg8qAm+7sU9l7csAOAV2CU/bjSzDiThuuVZrW1inmtqS4A2xfSUaYRADuT\nwhFEJay/nPq7SEQjHoHy/Z3O4C9tAfwDOfUCo829wpo9hmEYplhxr9YXku6cKGVDNuWYKQOwaJw3\nHUH3tkXR7EW8i8aZevr9d7iSTn952QT7c3tnCrNH9MYPTp0BIrK1VOpcfEjv4DkD1fZKIbkzlbK1\nXNmSobvPrSMWjeCzS8aguTajofvZGTPtz/VKoJqBzTX4zblzMHFgY0YISwk7dYPOZ2+ayx910sBG\nx/eaRHqfv77wvvYY0Qo143T77Ek6LJ9QVbMnzTiPntzfU18n2OkEwGQqhW1723JqK+A0DQ3ia6dr\ng5vuTLfBwl6RojXj7KnGMAzDMIyGv140DzOHpie83TldzXdQlK6Q7zx7XUE1DZMaMjVAi8R0/7bv\nc0YVVIUQd+TDjLCXqXPC9EG48/y52mPLue1ps4fakUYjLjPOCQMa7XNmS4YOZNeIxaMRDGquwW/P\nzZhWzhnZB2uuXYoLF43ChYtGafeTiwnJVCY6qM6M030Xz7dSWEiq485rcPtpnTN/OAAUxK+1mDE9\nthlDm3HC9IH47knT7LJWK9CT7h7phL29Lm1wLEJICmEH2gHMUTfdnD5nKACnZi+IFnafpl0HD+8V\n6JyFgIW9IkW+h0lTxjAMwzDFwLQhzRjRNz0J605TtGLU7BUDqhsQudxB1MTPJmFPaiQarUAZqgmd\nO/BF3BJ+1Hnz4F61mGNIdi+jt1bHI/bzk81oqonj/y44BD86bYYtPGZLhg5kX2CQwlXUYd6aFsKu\nWDrekz5EIrU3nSlhC526AC3qbexVG/doa2pcwp77vi+dPADrly9D3/oqVBKmd0VVLIofnDoDQ/vU\n2mWLx7XiN+fOxn8v9ArmOp8994JFNELoTAls2pXR7CVTwjdi5uwRvQFknvk/12zKtN0lKLqjEgN6\nzV7YSLX5hIW9IsX9IkxTRCMKwzAMwyAjTJhMswpyziKSsIpJy6gL+jDXEr4mDMgEyzH5PF6yZAw+\nOXcozpk/AkBaWLn/0gUAvP580mdP9dWs8fFZlJPiqlhUMTHNtOPg4b3RUB23Bcwgk+MGg7Amkedx\n+DIG6DtSrkumUrbQqcv3prZf1/tVzd6yKQPQ2lBZQp2JMFarRITDxrRoF3j2tyc9frLbXabIsQih\nM+nV7B35/UdhQi4SyD744KsbjW1XTYTVdrkJoqkuFCzsFSnyBdEeIC8IwzAMw/QUlRpRUFJEWSC0\nbRnWpw5PX7UElywZk3X/3nUJXHfCFFy4aBSuWDoeZ8wZhl7WZNat2RvcK619kZrdbKHk1bD6uiB0\nEttnL4uwd/6CkThu2kDfOtLUNOLQ7GWfTB08vDdGt9bjso+M8/XZU4+kW+vo11gNALjmoxPxkzNm\ndmtutWImFysA3WPb35HEkF61jrKte9sxcUDGdzIaIWzd04Z9LgFs3Za9MCEXTXTaXHfbdZo9nRln\nT1ojcK8rUmot51S1c7KsxzAMwxQbUsDoDjPOliLUjBRTpGxTOPd+jdWhNAtVsSguXDQKiVjEzhXn\n9nOaN7ovbjtvDi45fLTvuSXV1vmr43rNnkRqiLP5Qp4+e6hHS/fNEyfjDxccYrepzmq7OtGOBhD2\nGqrjePCyhZgyuEnJvxfF0N61OGxMJpG9eiidZrt3XQJrv3k0zpo3Ius5K4nRLekFguF9arPUzKDr\nKwc6kmiojuGBzy3ARy3Bf9veNoe2LRaNYPV7O7XHPPmgwdryjL+rVwvuDqajN+P0RoB9+LVNnrLu\nglMvFClSLa2GDC4mJ3CGYRiGAYBz54/Aw2s24yMT+xX8XHdfPB9vbt5T8POEoZiG5iCCTFjcfmcq\n80b3xT/XpE3csmkuZHTNqpii2dPsIyfTOn8sFfck/JLDR+P02UNBRJg4sBEpAZw+Z5inbthblNHs\nER794mIAwPAr700fSxH0TWsdrM3z8s0Tp+DEmYMwf3TfrBFVJTphb397ErEoYWy/Bswb1Qd3v/g+\ntu/twOSBGQEsGiG8uVmvxVs6uT/+77kNnnK5cCGQfv6q4iWQsNeRRJ+6BP5+6QIc/M0HAaQjv7r9\nCbuLgvZAIlpKRK8R0VoiutJQZxERvUBELxPRI4VsTymREfZYs8cwDMMUL6NbG/DUVUvQ2lBd8HP1\nb6rGvNF9s1fsRorJZy+axZQyF3Tmi45zRrymkjqqtJo9bz05l9ZFWlRxn++yI8fZi+K1iRi+cNQ4\nWyvnSDYf0pxOHkN7fapmL9RRK5uaRBSHjWkBEQX2v9VV29+R9Jjq7mnrdAhgchFiQFM1HvviYpw2\ne6j9XZrZupH+mUJ4F1DcGm7d739fe7pdqiXCnf+d/7yaQSmYsEdEUQA/AXA0gIkATiOiia46zQB+\nCuA4IcQkACcXqj2lhnRyVp08i2g8YRiGYRgGxSXsFcov6MvLJuDui+f7njPbuaUfVd/6Kl8zThnt\nMFuusjBazGhInz0VKaTqzuf02WNxr5Co1m2XHzUOQHphQPrVqc9n5rBe+M9VS/Do5Yvt5z1pYBOG\n9K7F+P7pQEVTFRNdN/KYKSFs4U6X5w8Adh/wauukxlFlbL8GfHHpuOwXWgAKqdmbDWCtEGKdEKId\nwB0AjnfVOR3An4QQ7wCAEKLnDFqLDJ1mj2EYhmGY4qKIAoOGFmSCVj/vsJGYMrhJu00Kedk0Zpd9\nZCwe/sIiLBnfardT554ifT+za/ayNltpY6Zy2Odla/Y0+x2rBIhhUa/7GN6nzv4sNXtq/zv5oMFo\nbazG0D619vMb1ZLeR1ZrrkloA7AAimYPmWizkwel+//IvnWOup+cOwzzR/fFoaMyKUf2dyS1x/7M\notGBrzGfFFLYGwTgXeX7BqtMZSyAXkT0LyJ6jog+XcD2lBQ1cV2AliIaURiGYRiGKSp/+qCaPZnX\nLR+Tz5gml52OSIQwom8dIhGyJ9y6XWTqhaw+e2E0e0rdsM/LDtbhytUHAEvGt2LlVz4CwOyzx+Sf\n6UObM76U1n/5WJdNGeB4xlLoklHu5by6sSZm+5ECmdySQCZdyZwRve3n2tpQhUcuX4RfnTPb0ZZ+\njdX47Xlz0F8xCd1zoLNH8+q56ekALTEABwFYAqAGwJNE9JQQ4nW1EhGdD+B8ABg6dGi3N7InkJGE\n1NDDRTSeMAzDMAyDItPsBWxMTSKK9cuX5eWcUmsWxheObGFPZ8aZ/v+JWUNw29PvGI8R5nxd8WWM\nacwEY9EI2jtTiEbIntQL1u11G801cUwe1ITn3t6OuOKrB6SFOBX52KSwt8syu2yojtsRYuX2XQfS\nx5g3ui9WXXMkGqvjju3D+ji1ekBmsUMVMLfubcfYfpm8lh+b4dZ1dS+FFDvfAzBE+T7YKlPZAOB+\nIcReIcQWAI8CmOY+kBDiRiHELCHErJaWloI1uJiYM6I3rj5mAq47YXJPN4VhGIZhGAPFlOC9J4gG\nNONUkZZKuj2kGeeFi0bhzeuPMR4jzH03mesFQWpL1fPFFZ9DeeyDh/fO+RxMOKIRsvM/SjPOXfst\nYa/aGR1Tauak8iRTz6nZq1c0e9EI2ceRIny1EpVWXSjJRO50Cvu96xJ23e+dMj3kFeaXQmr2ngEw\nhohGIC3knYq0j57KXwH8mIhiABIA5gD4fgHbVDIQEf5rwUhXWQ81hmEYhmEYLRUu6wU249Sh04VJ\nM854NOJ7zDBmnF0xqdP5hKW1fUlEKH3sez47H8P7erU+TGGIRyNoqkkLUw2WkCY1do2aVAhARliT\niwm966sciwANVRmRyGEO7RIW3ci6KVeUzl51+nb0BAUT9oQQnUR0MYD7AUQB/FII8TIRXWBtXyGE\neJWI/g5gFYAUgJuFEC8Vqk2lTjFF/GIYhmEYprh89nqSMJFA5S3T+bnJMnm8X541C0SEs295xlEv\nF+EyF3QBWqQfnxRMZfAOprD0rktg2952RCjzPPo3pX3lpg9pBpC2jFORfU0Ka5cfNQ696xI4ZnJ/\nRw5Et2ZPIjV21YZ8kzG7L7jaWpvQ1r/j/LkFi5proqA+e0KI+wDc5ypb4fr+XQDfLWQ7ygV3GFeG\nYRiGYXqWSl+I7bRmuWHug6zp5+cmzeMOH99Pu7277rvU/qjnO3f+SHz772vQUFU82ptK4IkrD8eG\n7ftBRDjQmXbulOaSR03qj+e/8hH0qnMKWXLxQAprzbUJfP5IbwqEuoSq2csIgZn9TZo9K02DS7PX\nbBD2ZPCX7qSnA7QwIYiFiTPMMAzDMEzBqXQzTmkWF2pBmjJJq01kO153afakX1d7MmWXXbhoFC5c\nNKpbzs9kqI5HMbq1HkAmD2ONonFzC3oqJjNMiUmzl9nfX7PXmUo5ytWE6j0NSw8lBGv2GIZhGKa4\nqHTN3qDmGgDAqQd7o6Xfe4k+EXtGs2cmayqHbrrtVZZmr70zlaUm0518ZvEo9K5LYE5ATVmVQTMn\nMfnsZQK0+PvsJV2avZEtxePDyZq9EiLOmj2GYRiGKQqiEUIyJSo+eFqf+ipjGgeTwGbfMx/VXjxL\nBM3u8pWssjRHbSzsFRUzh/aycxwGweRz595O5Iy8KjXXJs1e1CDsjSiigD0sPZQQ3WWywDAMwzCM\nPzEl/D6jx3RvZOoFP81emCmPSeuSD6T5X1sHC3ulTDYzTjnHHt1Sr91u6mNyUcIdoKU2UTz6tOJp\nCZOVOJtxMgzDMExRkIhG0NaZKjrN3pVHj8eYVv2EtbuYPKgRL723KydBuF9jFTbuagusuVtz7dLQ\n5wiDjMbZnkwW9DxMYcmWfkMKe/NG93WUuwO8mPZLpop3MYCFvRKCwzszDMMwTHEQj0WAtuLT7F2w\nsOcDh9xx/iHYvLsNwmCm6Zd64a8XzceaD3cFPlc287yuwpq98oDg/zvdZwV8kT6okuOnD8RfX3jf\nkZNPxeSzV0ywsMcwDMMwDBMSaW3DHhZe6qtiqK+K4e2te7Xb/VIv9G+qtnOnmZhm5VQLw/2XLsD7\nO/aH3k9qhNhnr7zZsS+dlL2p1plO439OnoavHDvR4cenIhUxqrD3PydPK1Arc4OFPYZhGIZhmJDI\ndEhsdWNGaj3d82Q/zV4Qzpk3PPQ+4/o3YFz/htD7ycAcbZ1sxlmK/OSMGbjp0bfslA1uzp43HBt3\nHcDO/e0AgOYap7AXj0bQtz57GoVOS9j744WH4KBhvbPU7l44QAvDMAzDMExIsvkAMZmohu4Ac1JA\nztXwze1XVUgWjm3BvNF9cMXS8d12TiZ/jG5twLdPmmoMcvi1j07CT884yNbsmZKhuzlsjLMPyqTq\nxWbWDbBmj2EYhmEYJjTSV6fTHYaPsZHza9MEOBfNninNQ6Goq4rhtvPmdus5me6nt5WQPWgy9JvP\nnGUndgeApNWZizFyPgt7DMMwDMMwIZEh1zuS7MtlImqbcbo1e+n/Op89E3dfPB+b9xzIW9sYRuU7\nJ03FcWu3BM6PVxWLOnLvyUWfYhT22AaBYRiGYboZIlpKRK8R0VoiulKzvYmI7iaiF4noZSI6uyfa\nyZiRAVpY2DMjzTU9Zpywpb3ATBnchMPH98tX0xjGQXNtAsdOHZjz/jJAi/TlLSaKr0UMwzAMU8YQ\nURTATwAcDWAigNOIaKKr2kUAXhFCTAOwCMANRBTMmYTpFi5cNBoAAmsCKpFoxD9iKRvAMuVCxoyz\nhxuioQibxDAMwzBlzWwAa4UQ64QQ7QDuAHC8q44A0EBp1Ug9gG0AOru3mYwfSyf3x/rlywIHdKhk\n3Jq9iB2Nk8U9pjyYPjidDqSxOp6lZvfDPnsMwzAM070MAvCu8n0DgDmuOj8GcBeA9wE0ADhFCMH2\ngkxJkTREKOxq6gWGKTau/9gUnD1vBFob/XNE9gSs2WMYhmGY4uMoAC8AGAhgOoAfE1GjriIRnU9E\nzxLRs5s3b+7ONjKML3VV6QAWJ8wY5CiXPnss6zHlQnU8iimDm3q6GVpYs8cwDMMw3ct7AIYo3wdb\nZSpnA1gu0nZua4noLQDjAfzHfTAhxI0AbgSAWbNm8fyZKRpqEzGsuuZI1CWc003W7DFM98GaPYZh\nGIbpXp4BMIaIRlhBV05F2mRT5R0ASwCAiPoBGAdgXbe2kmHyQGN13BiOPkzqBYZhcoM1ewzDMAzT\njQghOonoYgD3A4gC+KUQ4mUiusDavgLAtQBuJaLVAAjAFUKILT3WaIbJI2RIss4wTP5hYY9hGIZh\nuhkhxH0A7nOVrVA+vw/gyO5uF8N0BwkrR2GiGOPUM0yZwcJeCfCzM2aiKs4vRIZhGIZhSp9PHDwE\nG7bvx2eXjOnppjBM2cPCXglw9JQBPd0EhmEYhmGYvFAVi+JLx0zo6WYwTEXA6iKGYRiGYRiGYZgy\nhIU9hmEYhmEYhmGYMoSFPYZhGIZhGIZhmDKEhT2GYRiGYRiGYZgyhIU9hmEYhmEYhmGYMoSFPYZh\nGIZhGIZhmDKEhT2GYRiGYRiGYZgyhIU9hmEYhmEYhmGYMoSEED3dhlAQ0WYAb3fzafsC2NLN5ywG\n+LorC77uyqJUrnuYEKKlpxtRKvAY2a3wdVcWlXjdlXjNQOlcd6DxseSEvZ6AiJ4VQszq6XZ0N3zd\nlQVfd2VRqdfN5J9K7Ut83ZVFJV53JV4zUH7XzWacDMMwDMMwDMMwZQgLewzDMAzDMAzDMGUIC3vB\nuLGnG9BD8HVXFnzdlUWlXjeTfyq1L/F1VxaVeN2VeM1AmV03++wxDMMwDMMwDMOUIazZYxiGYRiG\nYRiGKUMqVtgjol8S0SYiekkpu5aIVhHRC0T0ABENtMqHE9F+q/wFIlqh7HMQEa0morVE9EMiop64\nnqDorlvZ9nkiEkTUVyn7knVtrxHRUUp52V53uTxvQx+/hojeU67tGGVb2T5r03WXy7MGzH2ciD5L\nRGuI6GUi+o5SXhbPm8k/PD7y+KhsK8vxEeAxksdIu7z8x0ghREX+AVgAYCaAl5SyRuXzJQBWWJ+H\nq/Vcx/kPgLkACMDfABzd09cW9rqt8iEA7kc6P1Nfq2wigBcBVAEYAeBNANEKuO6yeN6GPn4NgC9o\n6pb1s/a57rJ41j7XvRjAgwCqrO+t5fa8+a/b+hKPjzw+lv3z5jHSU7fcn3dFjJEVq9kTQjwKYJur\nbJfytQ6Ar0MjEQ1AegB8SqR7wK8BnJDvtuYT3XVbfB/AF+G85uMB3CGEaBNCvAVgLYDZFXDdWkrt\nun2uWUclPOvAlNF1XwhguRCizaqzySovm+fN5B8eHz3w+JiFMrtuHZXwvANTRtddEWNkxQp7Jojo\nm0T0LoAzAHxV2TTCUmE/QkSHWWWDAGxQ6mywykoKIjoewHtCiBddmwYBeFf5Lq+v3K8bKOPnDeCz\nlDbH+iUR9bLKyvpZW+iuGyjvZz0WwGFE9LR1fQdb5ZXwvJk8w+Ojg7L+DVXw+AjwGMljZJk9bxb2\nXAghrhZCDAFwG4CLreIPAAwVQkwHcBmA24mosafamE+IqBbAVXAO3GVPlusu2+cN4GcARgKYjvR1\n3tCzzek2TNddzs8aAGIAeiNtcnI5gN+XhH8BU5Tw+FgZVPD4CPAYyWNkGY6RLOyZuQ3AxwHAUuNu\ntT4/h7Tt7lgA7wEYrOwz2CorJUYhbY/8IhGtR/oaVhJRf6SvZYhSV15fWV93OT9vIcRGIURSCJEC\ncBOA2damcn7Wxusu52dtsQHAn0Sa/wBIAeiLMn/eTMHh8bG8f0MVOT4CPEbyGFmeYyQLewpENEb5\nejyANVZ5CxFFrc8jAYwBsE4I8QGAXUQ011oJ+DSAv3Zzs7uEEGK1EKJVCDFcCDEc6Y4/UwjxIYC7\nAJxKRFVENALp6/5PuV93OT9vy95cciIAGZWqbJ81YL7ucn7WFn9B2gEdRDQWQALAFpT582byD4+P\nPD6W+/gI8BhpwWNkuT1vUQRRYnriD8DvkFZPdyD9IjsXwB+R7uCrANwNYJBV9+MAXgbwAoCVAD6q\nHGeWtc+bAH4MK1F9sf7prtu1fT2sqFvW96uta3sNSsShcr7ucnnehj7+GwCrrT5+F4ABlfCsTddd\nLs/a57oTAH5rXcdKAIeX2/Pmv27rSzw+8vhY9s/bNFaU+/M2XXcFPO+KGCPJajjDMAzDMAzDMAxT\nRrAZJ8MwDMMwDMMwTBnCwh7DMAzDMAzDMEwZwsIewzAMwzAMwzBMGcLCHsMwDMMwDMMwTBnCwh7D\nMAzDMAzDMEwZwsIeU7IQ0feJ6FLl+/1EdLPy/QYiuizP59yTz+NZx5xORMco368hoi8E2C9JRC8Q\n0UDXsQQRLS1AO5uJ6DP5Pq5y/FOIaC0R3VOoczAMw1QCPD7y+MgwEhb2mFLm3wAOBQAiigDoC2CS\nsv1QAE/0QLvCMh3AMVlredkvhJguhHhfKTsNwOPW/3zTDEA7mBFRrKsHF0LcCeC8rh6HYRiG4fGR\nx0eGScPCHlPKPAHgEOvzJKSTXO4mol5EVAVgAoCVRFRPRA8R0UoiWk1ExwMAES0noovkwdQVQyK6\nnIieIaJVRPR13cl1dYhoOBG9SkQ3EdHLRPQAEdVY2w626r5ARN8lopeIKAHgGwBOscpPsQ4/kYj+\nRUTriOiSIDeDiAjAyQDOAvARIqrOpU1W+SQi+o9VvoqIxgBYDmCUUncRET1GRHcBeMXa7zLrul6S\nq8rW+dcQ0a1E9DoR3UZERxDRv4noDSKaHeT6GIZhmMDw+OhsD4+PTOXS01nd+Y//uvIH4C0AQwH8\nN4ALAFyL9CrgPACPWXViABqtz30BrAVAAGYAeEQ51isAhgA4EsCNVp0IgHsALLDq7LH+a+sAGA6g\nE8B0q97vAXzS+vwSgEOsz8sBvGR9PgvAj5V2XIP0QF1ltXcrgLjm2ve4vs8D8JD1+XYAH7c+59Km\nHwE4w/qcAFBjHecl5XyLAOwFMML6fhCA1QDqANQDeNm6x/L8U6x79RyAX1r37ngAf3Ed856e7lf8\nx3/8x3+l/sfjo+M7j4/8V7F/rNljSp0nkDZHORTAk9af/P5vqw4BuJ6IVgF4EMAgAP2EEM8DaCWi\ngUQ0DcB2IcS7SA9URwJ4HsBKAOMBjHGd16/OW0KIF6zPzwEYTkTNABqEEE9a5bdnua57hRBtQogt\nADYB6BfgXpwG4A7r8x1wmqqEbdOTAK4ioisADBNC7Dec8z9CiLesz/MB/FkIsVcIsQfAnwAcppx/\ntRAihfQg95AQQiA9+A0PcG0MwzBMOHh8zMDjI1OxdNmOmGF6GOmXMAXpVbh3AXwewC4At1h1zgDQ\nAuAgIUQHEa0HUG1t+z8AJwHoD+BOq4wAfEsI8XOf82rrENFwAG1KURLpVb+wuI/h+1sloiiAjwM4\nnoiuttrXh4gaDMfzbZMQ4nYiehrAMgD3EdF/A1inqbrX9yoyqOdPKd9T4PcQwzBMIeDxETw+Mgxr\n9phS5wkAxwLYJoRICiG2Ie0ofQgyzudNADZZA9liAMOU/e8EcCrSA9r/WWX3AziHiOoBgIgGEVGr\n67xB6tgIIXYg7S8xxyo6Vdm8G0CDd69QLAGwSggxRAgxXAgxDMAfAZyYS5uIaCSAdUKIHwL4K4Cp\nAdr5GIATiKiWiOqscz/WlYtiGIZhcobHxzQ8PjIVDQt7TKmzGmm7/adcZTv/f3t3rFpVEIQB+J8+\nL+A7pMwDpNdasLHxAdJY+AB29jZaiJBOsbdJFSQWIphAXkCsRBT7SXE2IJckJrnk3nDyfeWyZ5nq\nDLM7e85o8UiS3SRbVfUtyeMkx6cTu/so0wv6e3f/GGMfM7VsfBrPvMvCS/wyc87wJMmrqvqaqW//\n9xjfy3Th/N8L6Ff1KMmHhbH3+f9Xx86L6WGSwzG+meRtd/9Msj8ul79YXKi7vyR5k+RzkoMkr0cr\nEACrJz9O5EfutJragoGbVlUbo1c/VfUsyb3u3llivb/dvXGbYlpWVW0nedrdD9YVAwCrJT9eKp7t\nyI9cg5M9WJ37Y3fyMNPF7OdLrvenFn4aewtiuraxa/syya91xQDAWsiPF5AfWYaTPQAAgBlysgcA\nADBDij0AAIAZUuwBAADMkGIPAABghhR7AAAAM6TYAwAAmKETHcVrP0eI2iwAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(1,2,figsize=(15,4))\n", "ax = axes[0]; i = 50\n", "ax.plot(wavelength, X[i,:])\n", "ax.set_xlabel('Wavelength [Angstrom]'); ax.set_ylabel('Flux')\n", "\n", "ax = axes[1]; i = 500\n", "plt.plot(wavelength, X[i,:])\n", "ax.set_xlabel('Wavelength [Angstrom]'); ax.set_ylabel('Flux')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "Now, perform EM PCA on 2562 quasars. \n", "

\n", " 5. Using EM PCA, find the first 10 principal eigenvectors $\\phi_1, \\phi_2, ..., \\phi_{10}$ and reconstruct the data using them. ($ \\hat{X} = \\sum_{k = 1}^{10} c_k \\otimes \\phi_k $) For any two spectra, plot the original and reconstructed spectra.
\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "So far we treated all data equally when solving for the eigenvectors. However, we find that some data have considerably larger measurement noise, and they can unduly influence the solution. Now, we perform EM PCA with per-observation weights (called weighted EMPCA) so that the high $S/N$ data receive greater weight. (See https://arxiv.org/pdf/1208.4122.pdf for more detailed explanation. The following description is paraphrased from this paper.)\n", "

\n", "Basically, we add weights $w$ to the measured data in M-step: $\\phi = \\sum_j w_j\\ c_j\\ X_{row\\ j}$\n", "

\n", "In this case, the situation is more complicated since the measured flux in each wavelength bin for each quasar has a different weight. So we cannot do a simple dot product to derive $c$; instead, we must solve a set of linear equations for $c$. Similarly, M-step must solve a set of linear equations to update $\\phi$ instead of just performing a simple sum. Hence, the weighted EMPCA starts with a set of random orthonormal vectors and iterates over.\n", "

\n", "1. Initialize: Let $\\phi$ is a set of random orthonormal vectors." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Create an aray of random orthonormal vectors\n", "# Reference: https://github.com/sbailey/empca\n", "def _random_orthonormal(nvec, nvar, seed=1):\n", " \"\"\"\n", " Return array of random orthonormal vectors A[nvec, nvar] \n", " Doesn't protect against rare duplicate vectors leading to 0s\n", " \"\"\"\n", "\n", " if seed is not None:\n", " np.random.seed(seed)\n", " \n", " A = np.random.normal(size=(nvec, nvar))\n", " for i in range(nvec):\n", " A[i] /= np.linalg.norm(A[i])\n", "\n", " for i in range(1, nvec):\n", " for j in range(0, i):\n", " A[i] -= np.dot(A[j], A[i]) * A[j]\n", " A[i] /= np.linalg.norm(A[i])\n", "\n", " return A\n", "\n", "# Number of quasars\n", "nQSO = len(X)\n", "# Number of wavelength bins\n", "nLambda = len(wavelength)\n", "# Number of eigenvectors we want\n", "nEigvec = 10\n", "\n", "# A set of random orthonormal vectors\n", "phi = _random_orthonormal(nLambda, nEigvec, seed=1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "2. E-step: $X_{row\\ j} = \\phi\\ c_{col\\ j}$. ($X_{row\\ j}$ refers to $j$th row of $X$, and $c_{col\\ j}$ is $j$th column of $c$. Note that $X$ is a matrix of dimension \"nQSO\" x \"nLambda\", $\\phi$ is a matrix of dimension \"nLambda\" x \"nEigvec\", and $c$ is a matrix of dimension \"nEigvec\" x \"nQSO\".) Solve for $c$ assuming weights $w$.\n", "

\n", "We define weight $w$ as the inverse variance (\"ivar\"). (So $w$ is a matrix of dimension \"nQSO\" x \"nLambda\") This makes sense. \"We weight the measured data by the estimated measurement variance so that noisy observations do not unduly affect the solution, while allowing PCA to describe the remaining signal variance.\"\n", "

\n", "Now, solve $X_{row\\ j} = \\phi\\ c_{col\\ j}$ for $c_{col\\ j}$ with weights $w_{row\\ j}$. More generally, let $A = \\phi, x = c_{col\\ j}, b = X_{row\\ j}, w = w_{row\\ j}$:\n", "

\n", "$$ b = Ax $$
\n", "$$ wb = wAx$$
\n", "$$ A^T wb = (A^T w A)x$$
\n", "$$ (A^T w A)^{-1}A^T wb = x$$
\n", "

\n", "Hence, we get:

\n", "$$ c_{col\\ j} = (\\phi^T w_{row\\ j}\\ \\phi)^{-1}\\ \\phi^T w_{row\\ j}\\ X_{row\\ j} $$\n", "

\n", "In the below cell, we define the function \"_solve.\"
_solve(A, b, w) solves $Ax = b$ with weights $w$. This function solves $Ax = b$ with weights $w$ using $x = (A^T w A)^{-1}A^T wb$" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Solve Ax = b with weights w using the above set of equations\n", "# Reference: https://github.com/sbailey/empca\n", "import scipy\n", "def _solve(A, b, w):\n", " \"\"\"\n", " Solve Ax = b with weights w; return x\n", " \n", " A : 2D array\n", " b : 1D array length A.shape[0]\n", " w : 1D array same length as b\n", " \"\"\"\n", " \n", " #- Apply weights\n", " # nvar = len(w)\n", " # W = dia_matrix((w, 0), shape=(nvar, nvar))\n", " # bx = A.T.dot( W.dot(b) )\n", " # Ax = A.T.dot( W.dot(A) )\n", " \n", " b = A.T.dot( w*b )\n", " A = A.T.dot( (A.T * w).T )\n", "\n", " if isinstance(A, scipy.sparse.spmatrix):\n", " x = scipy.sparse.linalg.spsolve(A, b)\n", " else:\n", " x = np.linalg.lstsq(A, b)[0]\n", " \n", " return x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, in the E-step, for each QSO $j$, we can solve $X_{row\\ j} = \\phi c_{col\\ j}$ for $c_{col\\ j}$ with weights $w_{row\\ j}$ using the function \"_solve\".\n", "

\n", "Similarly in the M-step, for each wavelength bin $j$, we can solve $X_{col\\ j} = c^T \\phi_{row\\ j}$ for $\\phi_{row\\ j}$ with weights $w_{col\\ j}$ using the function \"_solve\".\n", "

\n", " 6. The below cell uses the weighted EMPCA to find $\\phi$. Fill in the blank and run the weighted EMPCA.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "C = np.zeros( (nEigvec, nQSO) )\n", "\n", "W = ivar\n", "\n", "# Number of iteration for EMPCA\n", "niteration = 20\n", "\n", "for jj in range(niteration):\n", " print(\"iteration\", jj+1, \"/20\")\n", " \n", " # E-step\n", " for i in range(...):\n", " b = ...\n", " A = ... \n", " w = ... \n", " ... = _solve(A, b, w) \n", "\n", " # M-step\n", " for j in range(nLambda):\n", " b = ...\n", " A = ... \n", " w = ...\n", " ... = _solve(A, b, w)\n", " \n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "Reconstruct the data using $\\phi$:

\n", "$$ \\hat{X} = (\\phi c)^T $$\n", "

\n", "$\\phi$ is a matrix of dimension \"nLambda\" x \"nEigvec\", and $c$ is a matrix of dimension \"nEigvec\" x \"nQSO\". So $\\hat{X}$ is a matrix of dimension \"nQSO\" x \"nLambda\" as expected.\n", "

\n", " 7. Reconstruct the data using the above equation. Remember that you chose two spectra in Part 5. For the same two spectra, plot the original and reconstructed spectra. Part 5 uses EMPCA without weights. Compared to Part 5, does your reconstructed spectra become less noisy?
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n", "The following analysis is based on https://arxiv.org/pdf/1605.04460.pdf.\n", "

\n", "In Part 6, we reconstruct the QSO spectra from the noisy data. This reconstructed spectra is closer to the true spectra of QSO. Note that in reality, the true spectra can never be directly observed, both due to measurement error and due to absorption by intervening matter along the line of sight. So we wish to perform inference about the true spectra of QSO using a non-parametric technique called Gaussian processes (GP). We henceforth call the measured spectra as $y(\\lambda)$ and the true spectra as $f(\\lambda)$ (where $\\lambda$ refers to wavelength).\n", "

\n", "A gaussian process is fully specified by its first two central moments: a mean function $\\mu(\\lambda)$ and a covariance function $K(\\lambda, \\lambda')$:

\n", "$$ \\mu(\\lambda) = \\mathbb{E}[f(\\lambda)\\ \\vert\\ \\lambda] $$
\n", "$$ K(\\lambda, \\lambda') = \\mathrm{cov}[f(\\lambda), f(\\lambda')\\ \\vert\\ \\lambda, \\lambda'] $$.\n", "
\n", "In this problem, we can derive the posterior distribution of $f$ conditioned on the observed values of $y$:

\n", "$$ p(f^*\\ \\vert\\ \\lambda^*, \\lambda, y, \\sigma(\\lambda)^2) = \\mathcal{N}(f^*\\ \\vert\\ \\mu_{f|y}(\\lambda^*), K_{f|y}(\\lambda^*, \\lambda^{*,})) $$\n", "
\n", "where $\\mathcal{N}(f\\ \\vert\\ \\mu, K)$ is a multivariate Gaussian given by:
\n", "$$ \\mathcal{N}(f\\ \\vert\\ \\mu, K) = \\frac{1}{\\sqrt{(2\\pi)^d \\mathrm{det}K}} \\mathrm{exp}\\big(-\\frac{1}{2}(f-\\mu)^TK^{-1}(f-\\mu)\\big) $$
\n", "where $d$ is the dimension of $f$.\n", "



\n", "In other words, for the QSO $i$, the measured spectra $y$ is $X_{row\\ i}$. Then, we can compute the posterior distribution of $f$ given $X_{row\\ i}$ as:\n", "

\n", "$$ \\mathcal{N}\\big(f\\ \\big\\vert\\ \\mu_{f|X_{row\\ i}}, K_{f|X_{row\\ i}}\\big) $$\n", "
\n", "where $\\mu$ is given by:\n", "

\n", "$ \\mu$ $ =\n", " \\begin{bmatrix}\n", " \\overline{x}_1 & \\overline{x}_2 & \\dots & \\overline{x}_{824} \\\\\n", " \\end{bmatrix}.$\n", "

\n", "The mean function $\\mu_{f|X_{row\\ i}}$ and the covariance function $K_{f|X_{row\\ i}}$ are defined as:

\n", "$$ \\mu_{f|X_{row\\ i}} = \\mu\\ +\\ K(K + V)^{-1}(X_{row\\ i} - \\mu) $$
\n", "$$ K_{f|X_{row\\ i}} = K - K(K + V)^{-1}K$$\n", "
\n", "where $K = \\phi \\phi^T$ (We can use $\\phi$ from Part 7. $\\phi$ is a matrix of eigenvectors, its dimension is \"nLambda\" x \"nEigvec\").

$V$ is a diagonal matrix whose entries are $\\sigma(\\lambda)^2$ e.g. for the QSO $i$, $V$ = np.diag(1/ivar[i,:]).\n", "

\n", "Finally, we can plot $f(\\lambda)$ by sampling from $\\mathcal{N}\\big(f\\ \\big\\vert\\ \\mu_{f|X_{row\\ i}}, K_{f|X_{row\\ i}}\\big)$.\n", "

\n", " 8. For any two spectra, plot $f(\\lambda)$ using Gaussian processes. You can use np.random.multivariate_normal (https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.random.multivariate_normal.html) to sample from a multivariate Gaussian.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "***" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## To Submit\n", "Execute the following cell to submit.\n", "If you make changes, execute the cell again to resubmit the final copy of the notebook, they do not get updated automatically.
\n", "__We recommend that all the above cells should be executed (their output visible) in the notebook at the time of submission.__
\n", "Only the final submission before the deadline will be graded. \n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "_ = ok.submit()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 2", "language": "python", "name": "python2" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.1" } }, "nbformat": 4, "nbformat_minor": 2 }