week4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 4 Assessment: Principal Component Analysis (PCA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning Objective\n",
    "\n",
    "In this notebook, we will implement PCA. We will implement the two versions of PCA as described in the lectures, which handles the when the dataset size exceeds the dataset dimensionality, as well as the case when we have the dimensionality greater than the size of the dataset.\n",
    "\n",
    "We will break down the task of implementing PCA into small components and combine them in the end to produce the final algorithm. We will apply PCA to the MNIST dataset and observe how the reconstruction changes as we change the number of principal components used.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import numpy as np\n",
    "import timeit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/matplotlib/__init__.py:1401: UserWarning:  This call to matplotlib.use() has no effect\n",
      "because the backend has already been chosen;\n",
      "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    }
   ],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import matplotlib as mpl\n",
    "mpl.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's plot a digit from the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD1CAYAAAB9TzjVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAEA5JREFUeJzt3X2MHPV9x/HPF5dn4jUQ1Zx4TKzlD4TgsHPUVUJxiRoR\ngzCBAjHIpg0lSAREqG0VUkNOiioQOggPKpYFWDGF2iBBhImh1oERKDxUC4cBAwrniwg5fMUYfIst\nJCjHt3/c3Hb37vY3593Zh/v5/ZJWN/P77sx9Wfy5mZ3ZnTF3F4A47dfqBgA0DgEHIkbAgYgRcCBi\nBByI2F80asXFYpHD80AT5XI5Gz9W1xbczM42sz+Y2TYzu6GedQHIXs0BN7MZkv5d0g8lnSRpsZmd\nlFVjADLg7jU9JP21pE1l8zdKunFsfnh42MceklySFwqF0nS7PeiN3qZjX+U5myyn9eyiHy3pz2Xz\ng8kYgDZhtX5U1cz+XtLZ7v5PyfwSSX/l7tdIlQfZ+vv7M2gVwHj5fL40PdlBNnbR6Y3epnFfjdxF\nL0jKm9m3zOwAST+WtKGO9QHIWM3nwd39KzO7RtImSTMkrXH3tzPrDEDd6vqgi7s/JempjHoBkDE+\nqgpEjIADESPgQMQIOBAxAg5EjIADESPgQMQIOBAxAg5EjIADESPgQMQIOBAxAg5EjIADESPgQMQI\nOBAxAg5EjIADESPgQMQIOBAxAg5ErGG3D0ac5s2bFxy75pprqi67dOnS4LoffPDBYP2ee+4J1vv6\n+oL1fRFbcCBiBByIGAEHIkbAgYgRcCBiBByIGAEHIsZ5cFTo7OwM1nt7eyvmBwYGKsZmzpxZdVl3\nD657yZIlwfp5550XrB955JHB+r6oroCb2fuSdksakfSVu38ni6YAZCOLLfjfuvvODNYDIGO8Bwci\nVm/AXdIzZvaamf00i4YAZMfSDnwEFzY72t0/NLO/lNQr6Vp3f0GSisViacX9/f11Nwpgonw+X5rO\n5XI2vl5XwCtWZNYtaY+790iVAZ81a5YkqVAoqKurK5PflzV6G5V2FH3z5s0V8wMDA5ozZ05pPnQU\nvV7FYjFYH38UvV3/n2bZ1/DwcGl6soDXvItuZoea2TfGpiX9QNLWWtcHIHv1HEWfLem3Zja2nv90\n9//KpCs0zOmnnx6sP/bYY8F6LpcLjoX2CHfv3h1c95dffhmsp53nnj9/ftWxtO+Kp/3u6armgLv7\nHyWdmmEvADLGaTIgYgQciBgBByJGwIGIEXAgYnxddBo65JBDqtbmzp0bXPahhx4K1js6OmrqaSrS\nPtF42223Bevr168P1l988cWK+b6+vtLYypUrg8vecsstwfp0xRYciBgBByJGwIGIEXAgYgQciBgB\nByJGwIGIcR58Glq9enXV2uLFi5vYyd5JO0d/2GGHBevPP/98sL5gwYKqtVNOOSW4bKzYggMRI+BA\nxAg4EDECDkSMgAMRI+BAxAg4EDHOg7ehefPmBcfOOeecqssml7GuWdq55ieffLJi/tJLL9Xy5ctL\n8z09PVWX3b59e3Ddr7/+erC+a9euYP2ss86aMLbffqPbsHpfl+mKLTgQMQIORIyAAxEj4EDECDgQ\nMQIORIyAAxHjPHgLdHZ2Buu9vb0V8wMDAxVjM2fOrLps6Pa9kvT0008H62nfJz/zzDMnjL333nul\n6dD1x++///7guj/++ONg/Y033gjWv/7666pjoc8OSOnfVU+7/XC7St2Cm9kaM9thZlvLxo4ws14z\n609+Ht7YNgHUYiq76L+RdPa4sRskPevueUnPJvMA2kxqwN39BUmfjhteJGltMr1W0vkZ9wUgA5b2\nnk2SzOwESb9z95OT+WF3n5VMm6RdY/NjisViacVp96QCUJt8Pl+azuVyEz5wX/dBNnd3Mwv+lejq\n6pIkFQqF0nS7aWZvaQfZNm/eXDE/MDCgOXPmlOZDB9nSZH2Qrbu7W93d3aX50MUN6z3IlmZkZKRi\nvq+vr3Tw7PPPPw8uO9nBw/HrykqW/9aGh4eD9VpPk31kZh2SlPzcUeN6ADRQrQHfIOnyZPpySU9k\n0w6ALKXuopvZOkkLJH3TzAYl/VLSrZIeNbMrJP1J0sWNbHK6OfHEE4P1FStWBOu5XC44tnPnzqrL\nDg0NBde9du3aYH3Pnj3B+saNGyvmu7u7K8bG19vFwQcfHKwvW7YsWL/sssuybKdpUgPu7tXelH0/\n414AZIyPqgIRI+BAxAg4EDECDkSMgAMR4+uiNTjwwAOD9dClgyVp4cKFwfru3buDY0uXLq267Kuv\nvhpcd9rpon3Vcccd1+oWGoItOBAxAg5EjIADESPgQMQIOBAxAg5EjIADEeM8eA1OO+20YD3tPHea\nRYsWVcz39PRUjKXd4hcYwxYciBgBByJGwIGIEXAgYgQciBgBByJGwIGIcR68BnfccUewPno3p+rS\nzmNPVufcd7r99pu4vRobm+zWwvsCtuBAxAg4EDECDkSMgAMRI+BAxAg4EDECDkSM8+BVnHvuuVVr\nnZ2dwWXdPVjfsGFDTT0hbLJz3WNjaf9PtmzZ0pCeWi11C25ma8xsh5ltLRvrNrMPzWxL8qjvCgcA\nGmIqu+i/kXT2JOO/dvfO5PFUtm0ByEJqwN39BUmfNqEXABmztPcmkmRmJ0j6nbufnMx3S/pHSUVJ\nr0pa5u67ypcpFoulFff392fWMID/l8/nS9O5XG7ClyBqDfhsSTsluaRfSepw95+UL1Me8FmzZkmS\nCoWCurq69v6/ognG9xY6yPboo48G13XAAQcE68uXLw/W77zzzmBv7aSdehsZGamY7+vr09y5cyWl\nH2RbtWpVsH7ttdfW11yZLF+z4eHh0vRkAa/pNJm7f+TuI+7+taT7JJ1ec4cAGqamgJtZR9nsjyRt\nrfZcAK2Teh7czNZJWiDpm2Y2KOmXkhaYWadGd9Hfl3RVA3tsidB9tNN2wXfs2BGsP/LIIzX1FLu0\n+653d3fXvO7NmzcH6zfeeGPN625nqQF398WTDD/QgF4AZIyPqgIRI+BAxAg4EDECDkSMgAMR4+ui\nDfDFF18E60NDQ03qpL2knQZbuXJlsL5ixYpgfXBwsOrY7bffHlx2z549wfp0xRYciBgBByJGwIGI\nEXAgYgQciBgBByJGwIGIcR68AfblyyKHLimddh77kksuCdafeOKJYP3CCy+smC8UCjr++OODy8SO\nLTgQMQIORIyAAxEj4EDECDgQMQIORIyAAxHjPHgVZhNuEjGlmiSdf/75wfp1111XU0/t4Prrrw+O\n3XTTTVWXzeVywXU//PDDwfrSpUtTusN4bMGBiBFwIGIEHIgYAQciRsCBiBFwIGIEHIgY58GrcPea\napJ01FFHBet33313sL5mzZoJY+Xfs/7kk0+qLjt//vzgupcsWRKsn3rqqcH6McccUzHf19ennp6e\n0vwHH3xQddlNmzYF133vvfcG69h7qVtwMzvWzJ4zs3fM7G0zuy4ZP8LMes2sP/l5eOPbBbA3prKL\n/pWkZe5+kqT5kn5mZidJukHSs+6el/RsMg+gjaQG3N2H3L0vmd4t6V1JR0taJGlt8rS1ksKfzwTQ\ndJb2frLiyWYnSHpB0smSPnD3Wcm4Sdo1Ni9JxWKxtOL+/v6M2gVQLp/Pl6ZzudyEL0lM+SCbmR0m\n6TFJP3f3z8q/cOHubmZV/1J0dXVJGr0I3th0uxnf20UXXVT1uevWrQuua2RkJFhfvXp1sD7+INt9\n992nK6+8sjTfbgfZ5s6dW5oPHWR75ZVXguu+6667gvW05cdr139vWfY1PDwcrE/pNJmZ7a/RcD/s\n7o8nwx+ZWUdS75C0o44+ATRA6hY82f1+QNK77n5HWWmDpMsl3Zr8DF/Tdh8yY8aMYP3qq68O1sdf\n/nf79u3auHFjaf6zzz6rumz5LlsjvPTSSxXzBx10UMXYc889V3XZm2++uWF9YXJT2UX/rqQlkt4y\nsy3J2C80GuxHzewKSX+SdHFjWgRQq9SAu/vvJVW7wsH3s20HQJb4qCoQMQIORIyAAxEj4EDECDgQ\nMb4uWsXLL79ctVYoFILL1vsppfFfN92+fXvF2OzZs2ted+hTcJK0fv36YH38JZ8LhYLOOOOMmvtB\nY7EFByJGwIGIEXAgYgQciBgBByJGwIGIEXAgYpwHr2JwcLBq7YILLggue9VVVwXrK1eurKmnqUi7\nKsqqVauC9W3btmXZDlqMLTgQMQIORIyAAxEj4EDECDgQMQIORIyAAxHjPHgNhoaGgvXu7u666uMV\nCoXUa60Dk2ELDkSMgAMRI+BAxAg4EDECDkSMgAMRI+BAxFIDbmbHmtlzZvaOmb1tZtcl491m9qGZ\nbUkeCxvfLoC9MZUPunwlaZm795nZNyS9Zma9Se3X7t7TuPYA1GMq9wcfkjSUTO82s3clHd3oxgDU\nb6/eg5vZCZJOk/TfydC1Zvamma0xs8Mz7g1Anczdp/ZEs8MkPS/p39z9cTObLWmnJJf0K0kd7v6T\nsecXi8XSivv7+zNtGsCofD5fms7lcjbhCe6e+pC0v6RNkv65Sv0ESVvLx4aHh33sodE/Al4oFErT\n7fagN3qbjn2V52yybE7lKLpJekDSu+5+R9l4R9nTfiRpa9q6ADTXVI6if1fSEklvmdmWZOwXkhab\nWadG/5K8Lyl8rWAATTeVo+i/lzRx3156Kvt2AGSJT7IBESPgQMQIOBAxAg5EjIADESPgQMQIOBAx\nAg5EjIADESPgQMQIOBAxAg5EjIADEZvyFV32VvkVXQA03mRXdGELDkSMgAMRa9guOoDWYwsORKxp\nATezs83sD2a2zcxuaNbvnQoze9/M3kpuwfRqi3tZY2Y7zGxr2dgRZtZrZv3Jz5Zcg75Kby2/hVXg\n9lotf91afeuvpuyim9kMSe9J+jtJg5IKkha7+zsN/+VTYGbvS/qOu+9sg17+RtIeSQ+6+8nJ2G2S\nPnX3W5M/joe7+7+0SW/dkva08hZWyRV+O8pvryXpfEn/oBa/boHeLlYTXrdmbcFPl7TN3f/o7l9K\nWi9pUZN+97Ti7i9I+nTc8CJJa5PptRr9B9J0VXprOXcfcve+ZHq3pLHba7X8dQv01hTNCvjRkv5c\nNj+o9rq/mUt6xsxeM7OftrqZScxO7hEnSf8jaXYrm5lE29zCatzttdrqdWvFrb84yDbqe+7eKemH\nkn6W7Iq2JR99T9VOpz5WSfq2pE6N3qTy9lY1ktxe6zFJP3f3z8prrX7dJumtKa9bswL+oaRjy+aP\nScbagrt/mPzcIem3Gn1L0U4+GruTTPJzR4v7KXH3j9x9xN2/lnSfWvTamdn+Gg3Qw+7+eDLcFq/b\nZL0163VrVsALkvJm9i0zO0DSjyVtaNLvDjKzQ5ODHzKzQyX9QO13G6YNki5Ppi+X9EQLe6nQDrew\nqnZ7LbXB69byW39N5eaDWTwkLdTokfQBSf/arN87hb6+LemN5PF2q3uTtE6ju2z/q9FjFVdIOlLS\ns5L6JT0j6Yg26u0/JL0l6U2NBqqjBX19T6O7329K2pI8FrbD6xborSmvG59kAyLGQTYgYgQciBgB\nByJGwIGIEXAgYgQciBgBByJGwIGI/R9qLKwhi6VZXgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f33a8f1d4e0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from ipywidgets import interact\n",
    "\n",
    "from sklearn.datasets import fetch_mldata\n",
    "MNIST = fetch_mldata('MNIST original', data_home='./MNIST')\n",
    "%matplotlib inline\n",
    "\n",
    "plt.figure(figsize=(4,4))\n",
    "plt.imshow(MNIST.data[0].reshape(28,28), cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before we implement PCA, we will need to do some data preprocessing. In this assessment, some of them \n",
    "will be implemented by you, others we will take care of. However, when you are working on real world problems, you will need to do all these steps by yourself!\n",
    "\n",
    "The preprocessing steps we will do are\n",
    "1. Convert unsigned interger 8 (uint8) encoding of pixels to a floating point number between 0-1.\n",
    "2. Subtract from each image the mean $\\mu$.\n",
    "3. Scale each dimension of each image by $\\frac{1}{\\sigma}$ where $\\sigma$ is the standard deviation of this dimension across the whole dataset.\n",
    "\n",
    "The steps above ensure that our images will have zero mean and one variance. These preprocessing\n",
    "steps are also known as [Data Normalization or Feature Scaling](https://en.wikipedia.org/wiki/Feature_scaling)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. PCA\n",
    "\n",
    "Now we will implement PCA. Before we do that, let's pause for a moment and\n",
    "think about the steps for performing PCA. Assume that we are performing PCA on\n",
    "some dataset $\\boldsymbol X$ for $M$ principal components. \n",
    "We then need to perform the following steps, which we break into parts:\n",
    "\n",
    "1. Data normalization (`normalize`).\n",
    "2. Find eigenvalues and corresponding eigenvectors for the covariance matrix $\\boldsymbol S$.\n",
    "   Sort by the largest eigenvalues and the corresponding eigenvectors (`eig`).\n",
    "\n",
    "After these steps, we can then compute the projection and reconstruction of the data onto the spaced spanned by the top $M$ eigenvectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def normalize(X):\n",
    "    \"\"\"Normalize the given dataset X\n",
    "    Args:\n",
    "        X: ndarray, dataset\n",
    "    \n",
    "    Returns:\n",
    "        (Xbar, mean, std): ndarray, Xbar is the normalized dataset\n",
    "        with mean 0 and standard deviation 1; mean and std are the \n",
    "        mean and standard deviation respectively.\n",
    "    \n",
    "    Note:\n",
    "        You will encounter dimensions where the standard deviation is\n",
    "        zero, for those when you do normalization the normalized data\n",
    "        will be NaN. Handle this by setting using `std = 1` for those \n",
    "        dimensions when doing normalization.\n",
    "    \"\"\"\n",
    "    mu = np.mean(X, axis=0) # EDIT THIS\n",
    "    X = X - mu\n",
    "    std = np.std(X, axis=0)\n",
    "    std_filled = std.copy()\n",
    "    std_filled[std==0] = 1.     \n",
    "    Xbar = (X-mu) / std_filled # EDIT THIS\n",
    "    return Xbar, mu, std\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def eig(S):\n",
    "    \"\"\"Compute the eigenvalues and corresponding eigenvectors \n",
    "        for the covariance matrix S.\n",
    "    Args:\n",
    "        S: ndarray, covariance matrix\n",
    "    \n",
    "    Returns:\n",
    "        (eigvals, eigvecs): ndarray, the eigenvalues and eigenvectors\n",
    "\n",
    "    Note:\n",
    "        the eigenvals and eigenvecs SHOULD BE sorted in descending\n",
    "        order of the eigen values\n",
    "        \n",
    "        Hint: take a look at np.argsort for how to sort in numpy.\n",
    "    \"\"\"\n",
    "    eig_val, eig_vec = np.linalg.eigh(S)\n",
    "    sorted_eig  = np.argsort(-eig_val)\n",
    "    eig_val = eig_val[sorted_eig]\n",
    "    eig_vec = eig_vec[:, sorted_eig]\n",
    "    return (eig_val, eig_vec) # EDIT THIS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def projection_matrix(B):\n",
    "    \"\"\"Compute the projection matrix onto the space spanned by `B`\n",
    "    Args:\n",
    "        B: ndarray of dimension (D, M), the basis for the subspace\n",
    "    \n",
    "    Returns:\n",
    "        P: the projection matrix\n",
    "    \"\"\"\n",
    "    P = B @ B.T # EDIT THIS\n",
    "\n",
    "    return P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, with the help of the functions you have implemented above, let's implement PCA! __When you implement `PCA`, do take advantage of the functions that you have implemented above.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn import decomposition\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def PCA(X, num_components):\n",
    "    \"\"\"\n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the data,\n",
    "           and N is the number of datapoints\n",
    "        num_components: the number of principal components to use.\n",
    "    Returns:\n",
    "        X_reconstruct: ndarray of the reconstruction\n",
    "        of X from the first `num_components` principal components.\n",
    "    \"\"\"\n",
    "    # Compute the data covariance matrix S\n",
    "    M = (X - np.mean(X.T, axis=1))\n",
    "    S = np.cov(X.T)\n",
    "\n",
    "    # Next find eigenvalues and corresponding eigenvectors for S by implementing eig().\n",
    "    eig_vals, eig_vecs = eig(S)\n",
    "    \n",
    "    # Select eigen vectors\n",
    "    U = eig_vecs[:, range(num_components)]\n",
    "    \n",
    "    # Reconstruct the images from the lowerdimensional representation\n",
    "    # To do this, we first need to find the projection_matrix (which you implemented earlier)\n",
    "    # which projects our input data onto the vector space spanned by the eigenvectors\n",
    "    P = projection_matrix(U) # projection matrix\n",
    "    \n",
    "    # Then for each data point x_i in the dataset X \n",
    "    #   we can project the original x_i onto the eigenbasis.\n",
    "    X_reconstruct =  np.dot(X, P)\n",
    "    return X_reconstruct\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Some preprocessing of the data\n",
    "NUM_DATAPOINTS = 1000\n",
    "X = (MNIST.data.reshape(-1, 28 * 28)[:NUM_DATAPOINTS]) / 255.\n",
    "Xbar, mu, std = normalize(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The greater number of of principal components we use, the smaller will our reconstruction\n",
    "error be. Now, let's answer the following question: \n",
    "\n",
    "\n",
    "> How many principal components do we need\n",
    "> in order to reach a Mean Squared Error (MSE) of less than $100$ for our dataset?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def mse(predict, actual):\n",
    "    return np.square(predict - actual).sum(axis=1).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "n = 1, reconstruction_error = 445.176488\n",
      "n = 2, reconstruction_error = 403.941979\n",
      "n = 3, reconstruction_error = 377.592431\n",
      "n = 4, reconstruction_error = 353.300124\n",
      "n = 5, reconstruction_error = 335.442895\n",
      "n = 6, reconstruction_error = 320.515895\n",
      "n = 7, reconstruction_error = 307.284207\n",
      "n = 8, reconstruction_error = 294.953574\n",
      "n = 9, reconstruction_error = 283.976615\n",
      "n = 10, reconstruction_error = 274.013351\n",
      "n = 11, reconstruction_error = 264.535178\n",
      "n = 12, reconstruction_error = 255.411519\n",
      "n = 13, reconstruction_error = 246.977773\n",
      "n = 14, reconstruction_error = 238.612489\n",
      "n = 15, reconstruction_error = 230.823945\n",
      "n = 16, reconstruction_error = 223.745087\n",
      "n = 17, reconstruction_error = 217.252885\n",
      "n = 18, reconstruction_error = 211.109805\n",
      "n = 19, reconstruction_error = 205.191266\n",
      "n = 20, reconstruction_error = 199.459376\n",
      "n = 21, reconstruction_error = 193.830646\n",
      "n = 22, reconstruction_error = 188.729930\n",
      "n = 23, reconstruction_error = 183.680207\n",
      "n = 24, reconstruction_error = 178.716447\n",
      "n = 25, reconstruction_error = 173.838899\n",
      "n = 26, reconstruction_error = 169.240087\n",
      "n = 27, reconstruction_error = 164.709929\n",
      "n = 28, reconstruction_error = 160.399040\n",
      "n = 29, reconstruction_error = 156.242176\n",
      "n = 30, reconstruction_error = 152.314700\n",
      "n = 31, reconstruction_error = 148.585512\n",
      "n = 32, reconstruction_error = 144.967513\n",
      "n = 33, reconstruction_error = 141.463610\n",
      "n = 34, reconstruction_error = 138.049052\n",
      "n = 35, reconstruction_error = 134.804194\n",
      "n = 36, reconstruction_error = 131.637698\n",
      "n = 37, reconstruction_error = 128.514868\n",
      "n = 38, reconstruction_error = 125.449114\n",
      "n = 39, reconstruction_error = 122.445685\n",
      "n = 40, reconstruction_error = 119.617802\n",
      "n = 41, reconstruction_error = 116.861161\n",
      "n = 42, reconstruction_error = 114.120986\n",
      "n = 43, reconstruction_error = 111.494982\n",
      "n = 44, reconstruction_error = 109.043607\n",
      "n = 45, reconstruction_error = 106.677698\n",
      "n = 46, reconstruction_error = 104.411285\n",
      "n = 47, reconstruction_error = 102.289186\n",
      "n = 48, reconstruction_error = 100.209621\n",
      "n = 49, reconstruction_error = 98.186397\n",
      "n = 50, reconstruction_error = 96.202592\n",
      "n = 51, reconstruction_error = 94.230480\n",
      "n = 52, reconstruction_error = 92.321572\n",
      "n = 53, reconstruction_error = 90.452334\n",
      "n = 54, reconstruction_error = 88.719880\n",
      "n = 55, reconstruction_error = 87.010172\n",
      "n = 56, reconstruction_error = 85.357146\n",
      "n = 57, reconstruction_error = 83.733070\n",
      "n = 58, reconstruction_error = 82.134603\n",
      "n = 59, reconstruction_error = 80.557153\n",
      "n = 60, reconstruction_error = 79.043511\n",
      "n = 61, reconstruction_error = 77.591119\n",
      "n = 62, reconstruction_error = 76.179282\n",
      "n = 63, reconstruction_error = 74.839740\n",
      "n = 64, reconstruction_error = 73.537097\n",
      "n = 65, reconstruction_error = 72.263037\n",
      "n = 66, reconstruction_error = 71.014394\n",
      "n = 67, reconstruction_error = 69.789344\n",
      "n = 68, reconstruction_error = 68.600844\n",
      "n = 69, reconstruction_error = 67.430337\n",
      "n = 70, reconstruction_error = 66.293432\n",
      "n = 71, reconstruction_error = 65.181586\n",
      "n = 72, reconstruction_error = 64.095697\n",
      "n = 73, reconstruction_error = 63.024775\n",
      "n = 74, reconstruction_error = 61.964528\n",
      "n = 75, reconstruction_error = 60.928472\n",
      "n = 76, reconstruction_error = 59.927952\n",
      "n = 77, reconstruction_error = 58.936304\n",
      "n = 78, reconstruction_error = 57.959615\n",
      "n = 79, reconstruction_error = 56.997041\n",
      "n = 80, reconstruction_error = 56.062799\n",
      "n = 81, reconstruction_error = 55.152727\n",
      "n = 82, reconstruction_error = 54.257452\n",
      "n = 83, reconstruction_error = 53.377407\n",
      "n = 84, reconstruction_error = 52.519199\n",
      "n = 85, reconstruction_error = 51.663204\n",
      "n = 86, reconstruction_error = 50.851429\n",
      "n = 87, reconstruction_error = 50.049822\n",
      "n = 88, reconstruction_error = 49.255963\n",
      "n = 89, reconstruction_error = 48.485538\n",
      "n = 90, reconstruction_error = 47.726620\n",
      "n = 91, reconstruction_error = 46.986703\n",
      "n = 92, reconstruction_error = 46.268002\n",
      "n = 93, reconstruction_error = 45.556610\n",
      "n = 94, reconstruction_error = 44.859183\n",
      "n = 95, reconstruction_error = 44.184838\n",
      "n = 96, reconstruction_error = 43.525914\n",
      "n = 97, reconstruction_error = 42.894227\n",
      "n = 98, reconstruction_error = 42.270628\n",
      "n = 99, reconstruction_error = 41.655122\n"
     ]
    }
   ],
   "source": [
    "loss = []\n",
    "reconstructions = []\n",
    "for num_component in range(1, 100):\n",
    "    reconst = PCA(Xbar, num_component)\n",
    "    error = mse(reconst, Xbar)\n",
    "    reconstructions.append(reconst)\n",
    "    print('n = {:d}, reconstruction_error = {:f}'.format(num_component, error))\n",
    "    loss.append((num_component, error))\n",
    "\n",
    "reconstructions = np.asarray(reconstructions)\n",
    "reconstructions = reconstructions * std + mu # \"unnormalize\" the reconstructed image\n",
    "loss = np.asarray(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   1.        ,  445.17648843],\n",
       "       [   2.        ,  403.94197923],\n",
       "       [   3.        ,  377.59243065],\n",
       "       [   4.        ,  353.30012434],\n",
       "       [   5.        ,  335.44289452],\n",
       "       [   6.        ,  320.51589546],\n",
       "       [   7.        ,  307.28420723],\n",
       "       [   8.        ,  294.95357442],\n",
       "       [   9.        ,  283.97661502],\n",
       "       [  10.        ,  274.01335119],\n",
       "       [  11.        ,  264.53517783],\n",
       "       [  12.        ,  255.4115194 ],\n",
       "       [  13.        ,  246.97777313],\n",
       "       [  14.        ,  238.61248858],\n",
       "       [  15.        ,  230.82394538],\n",
       "       [  16.        ,  223.74508667],\n",
       "       [  17.        ,  217.25288549],\n",
       "       [  18.        ,  211.10980526],\n",
       "       [  19.        ,  205.19126619],\n",
       "       [  20.        ,  199.45937602],\n",
       "       [  21.        ,  193.83064587],\n",
       "       [  22.        ,  188.72992977],\n",
       "       [  23.        ,  183.68020736],\n",
       "       [  24.        ,  178.7164468 ],\n",
       "       [  25.        ,  173.83889867],\n",
       "       [  26.        ,  169.24008665],\n",
       "       [  27.        ,  164.70992918],\n",
       "       [  28.        ,  160.39904031],\n",
       "       [  29.        ,  156.24217637],\n",
       "       [  30.        ,  152.31470003],\n",
       "       [  31.        ,  148.58551212],\n",
       "       [  32.        ,  144.96751302],\n",
       "       [  33.        ,  141.4636097 ],\n",
       "       [  34.        ,  138.04905222],\n",
       "       [  35.        ,  134.80419402],\n",
       "       [  36.        ,  131.63769824],\n",
       "       [  37.        ,  128.51486788],\n",
       "       [  38.        ,  125.44911449],\n",
       "       [  39.        ,  122.44568498],\n",
       "       [  40.        ,  119.61780193],\n",
       "       [  41.        ,  116.861161  ],\n",
       "       [  42.        ,  114.12098626],\n",
       "       [  43.        ,  111.49498167],\n",
       "       [  44.        ,  109.04360723],\n",
       "       [  45.        ,  106.67769811],\n",
       "       [  46.        ,  104.41128462],\n",
       "       [  47.        ,  102.28918562],\n",
       "       [  48.        ,  100.20962071],\n",
       "       [  49.        ,   98.18639674],\n",
       "       [  50.        ,   96.20259239],\n",
       "       [  51.        ,   94.23047999],\n",
       "       [  52.        ,   92.32157216],\n",
       "       [  53.        ,   90.45233407],\n",
       "       [  54.        ,   88.71988007],\n",
       "       [  55.        ,   87.01017231],\n",
       "       [  56.        ,   85.35714647],\n",
       "       [  57.        ,   83.73306986],\n",
       "       [  58.        ,   82.1346028 ],\n",
       "       [  59.        ,   80.55715323],\n",
       "       [  60.        ,   79.04351123],\n",
       "       [  61.        ,   77.59111949],\n",
       "       [  62.        ,   76.17928152],\n",
       "       [  63.        ,   74.83974007],\n",
       "       [  64.        ,   73.53709719],\n",
       "       [  65.        ,   72.26303721],\n",
       "       [  66.        ,   71.01439426],\n",
       "       [  67.        ,   69.78934379],\n",
       "       [  68.        ,   68.60084356],\n",
       "       [  69.        ,   67.43033735],\n",
       "       [  70.        ,   66.29343157],\n",
       "       [  71.        ,   65.18158647],\n",
       "       [  72.        ,   64.09569727],\n",
       "       [  73.        ,   63.02477536],\n",
       "       [  74.        ,   61.96452806],\n",
       "       [  75.        ,   60.92847182],\n",
       "       [  76.        ,   59.92795243],\n",
       "       [  77.        ,   58.93630438],\n",
       "       [  78.        ,   57.95961501],\n",
       "       [  79.        ,   56.9970415 ],\n",
       "       [  80.        ,   56.06279858],\n",
       "       [  81.        ,   55.15272708],\n",
       "       [  82.        ,   54.2574517 ],\n",
       "       [  83.        ,   53.37740707],\n",
       "       [  84.        ,   52.5191994 ],\n",
       "       [  85.        ,   51.66320389],\n",
       "       [  86.        ,   50.85142887],\n",
       "       [  87.        ,   50.04982152],\n",
       "       [  88.        ,   49.25596342],\n",
       "       [  89.        ,   48.48553832],\n",
       "       [  90.        ,   47.72661966],\n",
       "       [  91.        ,   46.98670333],\n",
       "       [  92.        ,   46.26800157],\n",
       "       [  93.        ,   45.55661019],\n",
       "       [  94.        ,   44.85918259],\n",
       "       [  95.        ,   44.18483827],\n",
       "       [  96.        ,   43.52591415],\n",
       "       [  97.        ,   42.89422674],\n",
       "       [  98.        ,   42.27062771],\n",
       "       [  99.        ,   41.65512189]])"
      ]
     },
     "execution_count": 163,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "loss"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also put these numbers into perspective by plotting them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAa0AAAEXCAYAAAAeBxskAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4m9WV+PHvseTdiZfE2ROSgFkSlkAhZSuEfS2BFmg6\n0IG2006naQszpQstvyktDdBOCy1DoRtQBlogZSmUsoVA2CGBLDQLiQJJSBwvceI1tmVLPr8/3tex\nJEuOJFuSl/N5Hj3Wu10dvZJ1dO97da+oKsYYY8xQkJXpAIwxxph4WdIyxhgzZFjSMsYYM2RY0jLG\nGDNkWNIyxhgzZFjSMsYYM2RY0jLDhogsE5E/ZjqOSCLyTRHZISJdInJjCsrfKiI3DHCZKTmXqYjV\njCyWtNJIRP4kIioij0fZNt/dFohYf5GIvC4ie0Rkr4hsFpE/i8hod/t097hot+vS9dxMdCIyCfgV\ncAswGfhFCh7mOOD2AS7zM8B/DXCZZoC5nwc3ZjqOdPJmOoAR6GPgQhEZr6o1Iev/HdgGTOleISKn\nA48DNwFfBfzAQcDFQG5EufOB5RHrmgY29JFJRHJUtSPJw2fifDl8SlWrBjCsfXGp6q6BLBdAVfcM\ndJnGDAhVtVuabsCfgBeBV4HvhayfBnQCPwICIet/Bby7nzKnAwqcnEAcZwFBYErE+s8BrcBod/kH\nwEc4yXIX8DyQ30e5W4GfAL8G9gA1ODUAb8g+y4A/Rhx3A7A1ynn6JrADaAH+CGQDX8NJ7vXA74Gc\niLLvBW4F6nCS9u+BvIjH+ybwAdAO+IAfRsS4FfgpcBewG3inj+d8PvCee45q3WMK3W03uq9N6G16\nH+dukfs8m9z4bway9heXu/6GRF4Hd7+FwPqQ2B+L9TrFc25x3lfL3MdsBF4B5kZ5njfEOp/uPgcC\nj7rltALvAxfGc85T/f5xy7gVqAQ63PP3LxHxK/B14AGg2Y3h+oh9st33xxac9+E64N8TKceNt9f7\nyy37Nnd/P1AFPJzuz7tU3TIewEi6hfwzXYnzYSnu+p8AzwFXE560vgc0RP7jR5Q5ncSTVpb7hv5e\nxPpngL+49z/j/tN+GiepzgGuZf9Jqx74PlABXI6TjL8css8y4ktaTcD9wGFuDO3As8D/uesuANqA\n/4gouwn4Q8hxtcDtIfvciPOhdQkwA+cD8GPgpojn0eTuezAwK8bzPRII4CSEQ4Hz3LIecLcXuedR\ngaOBCYCnj3PX5L4XDgG+AOwFrtlfXERPWvt7HX6M82H+DbesOcAPYr1OcZ7bS9zHOgSYjZMo9gBj\nImKLmbTcc1SD839yMk5N9ULgvHjOeRreP/+D84XhMve8/QDoAs4I2Ufd5/AVnAS80F13RkSM7wNn\n47wPP4fzv/7leMsBynCS3i/c8zYB8OA06+4A5uH87x4HXJvpz7+BumU8gJF0oydp5blv/NPcN9kO\nnA+3qwlPWgXAU+4btQp4Ergm4kNguru9FedDKPR2Qh+x3AqsDVke734YnOMu/yewCchO4PltxWkG\nC133LPBQyPIy4ktatYR/C/4Hzrff3JB1TwKPRpS9lZDEgNOs2g4UuuezFTg34vH/FWiIeB5L43i+\nDwDLI9bNx/kQO8Bdnue+PlP2U9ZW4LWIdTcD2/cXF9GTVszXwT0XbcB1fcQT9jrt79zGKCMLJ3le\nESvWKMfcBFT3UWY85zyV7x8/8PWIx38CeClkWYE7IvbZANzi3p/hxntoxD7/DayOtxx3eTNwY8Q+\nvwZewv1SPNxu1hEjA1S1Heef7ys43/i8wN+j7NeqqhfhvMmvx2mSuB7YKCKHRez+RZxvy6G3VX2E\ncT8wW0SOcZevwPlHf9FdXozTzLDN7UDyBREZFcfTWx2xvBMnISZqg4ZfR6oGNqqqP2LduIjjlqtq\nMGT5DZzrfwfifPvPBx4TkZbuG/A7oFhEykPLiSPG2ThNvaFeAQSYFcfxkd6KWH4DmNLd6SaBuKDv\n12E2zhenFxKMr69zi4jMEJEH3M4BTTi1lmLggAQe4xPAm6q6N8b2eM95Kt4/BwE5MR5/dsS6vs7/\nsW6870a8D3+AUzOOt5xY7gOOADaLyG9F5LMikrOfY4YM64iROb8HVgJTgftUtVNEou6oqltxvj3+\nSUR+iFMD+i5OoupWqaqb431wVd0gIu/i1DJWun8f7P6HVdVKETkUpzZ4OvD/gJ+JyCdVdXsfRUd2\nWFDCe6l24fzDhsqOUk5nlHKirUvki1f3vpfhnMNIoZ0PYn1oZlq8ce3vdUiFp3FqMwuB7W4Mr+N8\n0KdbKt4/iejr/Hf/PRGn5h+5X7zlRKWqq0VkBs41xtNwal43icjxqjrkO2dZTStDVHU9sAI4Caft\nP97j6on+DTEZ9wOfd2tbR+G094c+ll9Vn1PV7+J8cyvA6bnYH7XApIh1x0TbMUnHiYgnZPlEnCad\nD3EudrcDM1V1c5RbMFqBfVgHnBKx7lScD5Z1ScR+fMTyiThfRgb6g2Y9znk4O8HjYp5bERmDU9O5\nVVWfd9/f7ST+Pn0POFFECmNsH+hzHqmv989m9360x1+bwGO85/6dFuU9+GGC8XbgXGIIo6otqvqE\nqn4Lp2Z3mBvnkGc1rcw6B6dnUtTuxe7vL4pw2uO3uvevAg7H+fYUqkxEJkSs26uqzX08/kM4vYzu\nAVaq6r5/PBH5Ms6XmuU4F4jPAEbhfOD1x4vA3SJyGU7z5aXAp9zHGAhjgN+IyK9xLuLfBPyuu7lJ\nRG4GbhYRdWPx4iTko1X1ewk+1v8AK0XkdpwmxunA/wJ/VtWPk4h9jvua/wXng+YanBrugFLVFhH5\nJXCjiLQBS3CaTc9X1Vv6ODTmuXXL2QV8RUQ+dPf9Oc61s0TchfPzjydF5Ec4zWGzgaCqPsvAn/NI\n+3v/3IFTa9kFrMF5/87HqdXERVU3i8i9wB9E5Ls4zcKFOE2j5ar6swTi3QKcJCLTcGpte4Bv45y3\n1e66z+P0Fo7WujDkWNLKIFVtpXfzQKhXcLro3ofTM6gF2Ahcqap/jtj3ySjH/wand1isx98tIv/A\nqT1dG7G5HrgO54MnF6fr+1dVdWkf8cbjfpyk+xucZqM/A3fgNE8OhEdxugd3N0s9gtOLDgBVvUlE\nqnDOyy9xPlQ34TS/JkRV3xeRi3A+2L6Ocw3nUZzzloz/xbn+8y5OU9ad9P5yMlD+H06S+RZOT7x6\nel+riRTz3Kpql/tF5A6cXnHbcK7RJPIBjKpWicjJ7nHP4DQd+3Cu5abinEfq8/2D8/OILpyfo5Tj\n1L6uTOL/4qs4yeWHOMmxCaemeGeC5fwI51LDRpzrlDPcsv4L5/pYFk7njc+q6sYEyx6UurtcG2My\nSES24vTW+2mmY4lGRJYBm1X13zIdS6qMhOc4HNg1LWOMMUOGJS1jjDFDhjUPGmOMGTKspmWMMWbI\nGLK9BxsbG62KaIwxw1xxcXHYYARW0zLGGDNkWNIyxhgzZIz4pOXz+QZFGQNVjsUyuMsYqHIslsFd\nxkCVMxxj6a8Rn7SMMcYMHWlNWiLiEZFVIvK0u3yjiFSKyGr3dn7Ivte7UxxsFJFz0hmnMcaYwSnd\nvQevwRkHK3R+oNtV9RehO4nILGABzkCZk4AXReTgJEbhNsYYM4ykraYlIlNwJjyMZxqO+cDD7tQY\nW3AGpZybyviMMcYMfmkbEUNEHgVuwZne4jpVvdCdhuGLQCPOyNbfVtV6EbkTeFtVH3SPvQd4VlUf\n7S4v9Hdag+UCoTHGmP6rqOiZwDnyd1ppaR4UkQuBWlV9T0TmhWy6G2eKAXX//hL4UqLlhz7BRPl8\nvn4dP1BlWCyDP5bh9nwsltSVYbGkTrqaB08CLnKnX3gYOF1EHlTVGlUNqmoX8Ad6mgArcaah7zbF\nXTdg2gPKmt0dPFPr4Werh/wM1MYYMyKkpaalqtfjTuLm1rSuU9UrRWSiqla5u11Cz5TVTwF/EZHb\ncDpiVODMoDsgAl3KjL9U0RZUnPkNm/nqYUWU5tovAIwxZjDL9NiDPxeROTjNg1txptlGVdeJyGKc\nqd0DwMKB7DnozRJmjPawvj6wb92G+k5OnJA7UA9hjDEmBdKetFR1GbDMvf+FPvZbBCxKVRyzSrPD\nktZ6S1rGGDPojdj2sMNKssOWNzQEYuxpjDFmsBixSWtWaXglc319Z4YiMcYYE68Rm7QOK42oadV3\nYrM4G2PM4DZik9a0Ig+F3p7frDV0KNVtXRmMyBhjzP6M2KSVJcKhJdZEaIwxQ8mITVrQu4nQkpYx\nxgxuIzppzep1Xct6EBpjzGA2wpNWePPghgaraRljzGA2opNW5G+1PqgP0GU9CI0xZtAa0UlrXH4W\nxd6eJNUWVLY22zyTxhgzWI3opCUiHFgQ3s3dOmMYY8zgNaKTFsBBheFJa4MlLWOMGbRGfNI6sCD8\nGpaNQWiMMYOXJS2raRljzJAx4pPWzIhrWr7GAB1B60FojDGDUVqTloh4RGSViDztLpeJyBIR8bl/\nS0P2vV5ENovIRhE5J1UxjfLC5ALPvuWAOonLGGPM4JPumtY1wIaQ5e8DS1W1AljqLiMis4AFwGzg\nXOAuEfGQIpE/Mv7nHmsiNMaYwShtSUtEpgAXAH8MWT0fuN+9fz9wccj6h1XVr6pbgM3A3FTFduSY\n8B8Zv1btT9VDGWOM6Yd01rR+BXwXCL2INF5Vq9z71cB49/5kYHvIfjvcdSlxysS8sOVllX6bW8sY\nYwYhSceHs4hcCJyvql8XkXnAdap6oYg0qGpJyH71qloqIncCb6vqg+76e4BnVfXR7n0bGxv3Be7z\n+foVX0cXnP52Pv6unvm1Fh/TxowCS1zGGJNuFRUV++4XFxdL6DZvr71T4yTgIhE5H8gDRovIg0CN\niExU1SoRmQjUuvtXAlNDjp/irosq9AkmyufzMfuQCk7eVsfSyp5mwY+yJ3B2RVHcZfQnhoEsx2IZ\n3GVYLIM/luH2fAZbLP2VluZBVb1eVaeo6nScDhYvqeqVwFPAVe5uVwFPuvefAhaISK6IzAAqgOWp\njHHepNyw5Zd32nUtY4wZbNJV04rlVmCxiHwZ2AZcDqCq60RkMbAeCAALVTWlI9meNikPaNq3/EaV\nn84uJTtLYh9kjDEmrdKetFR1GbDMvb8bOCPGfouARemKa3apl3H5WdS2Of1EWgLKitoOTpyQu58j\njTHGpMuIHxGjm4hYE6ExxgxylrRCOE2EPV7e2Z6hSIwxxkRjSStEZE1rZV0nDf6uGHsbY4xJN0ta\nISYWeDispOcyX5fCq1XWRGiMMYOFJa0IkbWtZXZdyxhjBg1LWhEir2strWy3IZ2MMWaQsKQV4aQJ\nOeSEnJVtLUGbqsQYYwYJS1oRCrOzODnit1nP77BehMYYMxhY0ori7KnhTYQvbLekZYwxg4ElrSjO\nmRKetN6q6aCxw7q+G2NMplnSimLGaC8VxT1d3wNqvQiNMWYwsKQVw9kRta3nrYnQGGMyzpJWDJFJ\na8mOdrqs67sxxmSUJa0YThifw6jsnmlJdrV3sbquM4MRGWOMsaQVQ45HOG2SdX03xpjBxJJWH3p1\nfbekZYwxGZWWpCUieSKyXETWiMg6Efmxu/5GEakUkdXu7fyQY64Xkc0islFEzklHnJHOmhyetFbV\ndVLTmtIJlI0xxvQhXTUtP3C6qh4FzAHOFZHj3W23q+oc9/YMgIjMAhYAs4FzgbtExJOmWPcZX+Dh\n6LHZYev+trUt3WEYY4xxpSVpqaPFXcx2b311xZsPPKyqflXdAmwG5qY4zKgumJYftnz/pr02gK4x\nxmSIpOsD2K0pvQccBPxGVb8nIjcCXwQagXeBb6tqvYjcCbytqg+6x94DPKuqj3aX19jYuC9wn8+X\nsrhr/cKnV+TRRU9PwvuOaufwUTZChjHGpEJFRcW++8XFxRK6zdtr7xRR1SAwR0RKgCdE5HDgbuAm\nnFrXTcAvgS8lWnboE0yUz+fr8/gK4Jzq3Twb8uPipa1lXHJMadxlDFQs6SrDYkldGRbL4I9luD2f\nwRZLf6W996CqNgAvA+eqao2qBlW1C/gDPU2AlcDUkMOmuOsy4upDCsOWH9/SZmMRGmNMBqSr92C5\nW8NCRPKBs4APRGRiyG6XAGvd+08BC0QkV0Rm4FR4lqcj1mjOnJzL5IKefiCtAeXRj1ozFY4xxoxY\n6appTQReFpH3gRXAElV9Gvi5iPzTXX8a8J8AqroOWAysB54DFrrNixnhyRK+cHBB2Lr7NrZahwxj\njEmztFzTUtX3gaOjrP9CH8csAhalMq5EXFlRwM/XNNPl5qm1ezpZVdfJMeU5mQ3MGGNGEBsRI05T\nirycFTGI7r0b92YoGmOMGZksaSXg6ogmwsUftlJtI2QYY0zaWNJKwFlT8phW1NMho6MLfrOupY8j\njDHGDCRLWgnwZgnXHFEUtu7eD/bSYDOWGGNMWljSStAVBxUyLr/ntO0NKIt3ZvdxhDHGmIFiSStB\neV7hG7PDa1uPVHlp7rQfGxtjTKpZ0krCFw8tpDinZzispoDwpw+sJ6ExxqSaJa0kjMrO4t9nhde2\n7lzXQnvAfmxsjDGpZEkrSV87rJBCb09tq6ati3vsd1vGGJNSlrSSVJbn4YsRA+n+Yk2TDaRrjDEp\nZEmrH649sohR2T21rXq/csc/mzMYkTHGDG+WtPphbJ6Hbx0efm3rrnV7bZQMY4xJEUta/fT12UWU\nZfd0wGgLKj9b3ZTBiIwxZviypNVPhdlZfGVa+JAY/7epFV+jDZNhjDEDzZLWALh4fICZo3rGJAwq\n/Phdq20ZY8xAS9fMxXkislxE1ojIOhH5sbu+TESWiIjP/Vsacsz1IrJZRDaKyDnpiDNZ3iz4708U\nh617+uN2Xtnpz1BExhgzPKWrpuUHTlfVo4A5wLkicjzwfWCpqlYAS91lRGQWsACYDZwL3CUinqgl\nDxLzp+dxbHn4GITXL28g0GU/ODbGmIGSlqSlju45PLLdmwLzgfvd9fcDF7v35wMPq6pfVbcAm4G5\n6Yg1WSLCrZ8sCVu3vj7A/21qzVBExhgz/KTtmpaIeERkNVALLFHVd4Dxqlrl7lINjHfvTwa2hxy+\nw103qB1bnsPnDswPW/fTlU00+O0Hx8YYMxBENb3NVyJSAjwBfBN4XVVLQrbVq2qpiNwJvK2qD7rr\n7wGeVdVHu/dtbGzcF7jP50tb/PtT6xcufS+Ptq6eHx0vmNTJt2dab0JjjIlHRUXFvvvFxcUSus2b\n7mBUtUFEXsa5VlUjIhNVtUpEJuLUwgAqgakhh01x10UV+gQT5fP5+nV8ZBkVwHXBZm5a2dN78K9V\n2XzjuMkcXtb3vFsDHUumyxlusQy352OxpK4MiyV10tV7sNytYSEi+cBZwAfAU8BV7m5XAU+6958C\nFohIrojMwMkFy9MR60BYOLuIaUXhXeC/+Ua9dcowxph+Stc1rYnAyyLyPrAC55rW08CtwFki4gPO\ndJdR1XXAYmA98BywUFWHzNhIeV7hlrnhXeBX1XVy9/qWGEcYY4yJR1qaB1X1feDoKOt3A2fEOGYR\nsCjFoaXMBQfkc/H0fP62tW3fuptXNnPBtHxmjk57q6wxxgwLNiJGCv38+GJKQmY4bgsq33qjnnR3\nfjHGmOHCklYKjcv3cEvEb7der+7gTxvtt1vGGJMMS1optuDAfM6YnBu27ocrGm1AXWOMSYIlrRQT\nEW4/sYRCb08zYWtA+cor9XQErZnQGGMSYUkrDaYVebn1k+G9CVfv7mTRShsJ3hhjEmFJK02urChg\n/vS8sHV3rG2xkeCNMSYBlrTSRET49YmlTC7o+dGxAl97bQ+724fMT9CMMSajLGmlUUluFr89pZTQ\ngbSqWrv42qv1dFk3eGOM2S9LWmn2qYm5XHtEUdi6JZV+fv1PGy3DGGP2x5JWBvzgmNF8clxO2Lqf\nrmxiVaO9HMYY0xf7lMyA7CzhnlNLKcvtOf1BhR9uzKHOrm8ZY0xMlrQyZEqRl99+qjRs3a6OLP7t\nlXqCNhq8McZEZUkrg86emtfr+taynX5uWdWcoYiMMWZws6SVYT88ZjQnjA+/vvWL95t55uO2GEcY\nY8zItd+kJSKXRSwfErF87UAHNZJkZwn3zStjfH74S/G1V+v5qCmQoaiMMWZwiqemdU/E8lsRyz8Z\noFhGrAkFHu6bV4aHnmtZTZ3KlS/tprmzK4ORGWPM4BJP0pIEl3sXIDJVRF4WkfUisk5ErnHX3ygi\nlSKy2r2dH3LM9SKyWUQ2isg5ccQ5pJ04IZdvzQgf+X19fYCvWMcMY4zZJ56kFfmJub/laALAt1V1\nFnA8sFBEZrnbblfVOe7tGQB32wJgNnAucJeIeKIVPJx8flKAz8zID1v33PZ2fvSuDaxrjDEQZ0cM\ncWR1J47I5f1R1SpVXenebwY2AJP7OGQ+8LCq+lV1C7AZmBvPYw1lIvC/J5VwZFl22Po717Vw/8a9\nGYrKGGMGD9nf1O8i0kV4bUpClgVQVY27FiQi04FXgcOB/wK+CDQC7+LUxupF5E7gbVV90D3mHuBZ\nVX20u5zGxsZ9Mfl8vngffkio9QtXr8llV0fPdwqPKP87289xJXaNyxgzvFVUVOy7X1xcHHYJyhvH\n8TMGKhARKQIeA65V1SYRuRu4CScJ3gT8EvhSouWGPsFE+Xy+fh0/UGWEllMB/HVCB+c9U0ebO1Fk\nUIXrN+Wz5IJyDi7J3m8ZAxVLpssYTLEMt+djsaSuDIsldfbbPKiq26LdgKaQ+/slItk4CevPqvq4\nW3aNqgZVtQv4Az1NgJXA1JDDp7jrRow5Y3P43SnhI2Y0diiXv7jbpjIxxoxY8fxO619De++JyLEi\nsh2oc3v2HdLH4d3HCE7X+Q2qelvI+okhu10CrHXvPwUsEJFcEZkBVADL43pGw8hF0/P50SdGh63b\n2hzkiqV78AetR6ExZuSJpyPGdUB1yPLvgReBI92//xNHGScBXwBOj+je/nMR+aeIvA+cBvwngKqu\nAxYD64HngIWqOiKrF9ceUcQVFQVh696u7eAbr9scXMaYkSeea1pTgX+C83sr4AjgTFXdIyLfx+nZ\n1ydVfZ3ov+d6po9jFgGL4ohvWBMRbj+hhG3NAV6v7ti3/q8ftTEu38NPjxuNU5E1xpjhL56aVgDo\nHhzvROADVd3jLrcC+VGPMgMmxyM8cPoYDhod/h3jN+tauGOtTR5pjBk54klarwCLRORI4JvA30O2\nHUp406FJkdLcLBafNYbyvPCX7EfvNvHAJvsNlzFmZIgnaV0DHA28gVOz+lnIti/gXHMyaTBztJe/\nnjWGUdnhzYHXvNnAk1ttVHhjzPAXzzUtD3A1PT8qLhaRYnfbXSmKy8QwZ2wOfz5jDJe+UEeH+zvj\nLoUvLdvDvfPKmNX34cYYM6TFU9PaCmwBPnL/bo24bUlFYCa2Uybm8sdTy8gKqXAF3cT1Yt2wH6LR\nGDOCxZO01gA+4AZgOpAdccuJeaRJmYum53PnSSVhXTKDCjd8kMMTW1ozFpcxxqRSPCNiHA1cCpTh\nXNd6BmcE9hx3NIsR+fupweBfKgq561Ol4YkL4cuv1PMnG2DXGDMMxTXKu6quVdXv4NS0bgMuBKpE\n5JgUxmbi8PmDCvjtKaVhTYVdCte+2cDPVzexvwGRjTFmKIkraYWoAE4FTgBWAfUDHpFJ2OcOLOC3\nnyrFE/Eb45tXNfPdtxttEkljzLARz9iDZSKyUESWA38DWoBTVPU0d64rMwhcfmABD55eRm5WeIL6\nwwd7uerlPbQGbEoTY8zQF09NayfwDZyEtRB4GzhIRE7vvqUyQBO/86blc+dsP8U54VWupz9u59PP\n1lHbZpcfjTFDWzy/06oG8oCvuLdICswcyKBM8uYUd/HMeeVcuqSOqtae2tV7dZ2c8fQu/nrWGA7t\nYz4uY4wZzOLpPThdVWf0cbOENcjMLstmyQXlzCoN/06yvSXI2U/v4sUd7RmKzBhj+ifRjhhmiJhS\n5OW588s5fVJu2PqmTmciybvXtVjPQmPMkGNJaxgbnZPFI2eN4aqDw+fj6lK4fnkj177ZQIdNJmmM\nGULSkrREZKqIvCwi60VknYhc464vE5ElIuJz/5aGHHO9iGx2Z0c+J3bppi/ZWcKvTizhp8eNDvst\nF8D9m1qZ/7x10DDGDB3pqmkFgG+r6izgeGChiMwCvg8sVdUKYKm7jLttATAbOBe4S0RsUL0kiQjf\nOHwUD5/Re4T4t2o6mPdULSt3dcQ42hhjBo+0JC1VrVLVle79ZmADMBmYD9zv7nY/cLF7fz7wsKr6\n3d+CbQbmpiPW4ezsqXksubCc6aPC8//O1i7Oe3YXf/HZ0E/GmMEt7de0RGQ6zvxc7wDjVbXK3VQN\njHfvTwa2hxy2w11n+unQkmxe/vQ45kV00PAH4euvN/Cdt+w6lzFm8JJ09iATkSLcmZBV9XERaVDV\nkpDt9apaKiJ3Am+r6oPu+nuAZ1X10e59Gxsb9wXu8/nS9hyGi4DCnVuz+XNl799sHTU6yK2H+hlr\n4/cbYzKgoqJi3/3i4uKwaxrx/Lh4QIhINvAY8GdVfdxdXSMiE1W1SkQmArXu+kpgasjhU9x1UYU+\nwUT5fL5+HT9QZWQilt8cDKd+2Mq33qinPaQvxpomD1euyufe08dy6qS8tMSS6nIGSxkWy+CPZbg9\nn8EWS3+lq/egAPcAG1T1tpBNTwFXufevAp4MWb9ARHJFZAbOQL3L0xHrSHP5gQU8f0E504rCr3Pt\n7hQufn43/72i0ZoLjTGDRrquaZ0EfAE4XURWu7fzgVuBs0TEB5zpLqOq64DFwHrgOWChzduVOkeN\nyWHZp3v/EFmBO9a2cNY/duFr7MxMcMYYEyItzYOq+jogMTafEeOYRcCilAVlwpTlefjrWWO4ZXUz\nv1zTTGjdas3uTk59ahc3zy3mqoMLcCrOxhiTfjYihtnHkyXccMxo/n7eWMblhE9l0hpQrn2zgSte\n2kNdu1V6jTGZYUnL9HLyhFweOqadi6fn99r2zMftnPS3WpbYoLvGmAywpGWiGu2F++aV8puTSyjy\nhjcH1rThC8xuAAAdfElEQVR1cdmS3Vz3VoNNLmmMSStLWiYmEeGKikJemz+OY8t7/57rjx/s5ZQn\nd9kQUMaYtLGkZfZrxmgvz55fznfnjMIT0Qdjc1OAs/6xi5tXNdHZZV3jjTGpZUnLxCU7S/jB0aN5\n7vxyZkaMXRhU+PnqZs58ehcb6q1rvDEmdSxpmYQcNy6HV+eP4+qIObqgu2t8Lb9Y02y1LmNMSljS\nMgkrys7iVyeV8siZYxifH/4W6uiCn65sYt5TtaxvtreXMWZg2aeKSdo5U/N46+JxXBKla/y6+gBf\nXJPLDcsb2dtpPQyNMQPDkpbpl7I8D/edVsa9p5YyNi/87dSFcOe6Fk78Wy0vVdrvuowx/WdJywyI\nz8ws4J1LxnH5gb1rXdtagnzmhd189ZU91LbZaBrGmORZ0jIDZkyeh9+fUsbiM8cwpdDTa/vij9o4\n9vEafr++haB11DDGJMGSlhlwZ0/N4+1LxrFgUmevUZKbOpTvvtPIaX/fxTs1/ozEZ4wZuixpmZQo\nys7i2zM7efHCcmaX9p5M4P09nZzzTB1ffXUPVa3WZGiMiY8lLZNSnyjPYdlF41g0t7jXGIYAiz9s\n49jHarjt/WbaAtZkaIzpmyUtk3LZWcLC2UWs+Ox4Lp3Zu6PG3oDyk/eaOO7xGh7a3EqXWvIyxkSX\nlqQlIveKSK2IrA1Zd6OIVEbMZNy97XoR2SwiG0XknHTEaFJvYoGHP55axt/PHcusKE2GO/YG+Y/X\n6pn31C6WVrajlryMMRHSVdP6E3BulPW3q+oc9/YMgIjMAhYAs91j7hKR3l3RzJD1qYm5vHrROH5x\nfDElOb2bDN/f08lnX9jNhc/VWWcNY0yYtCQtVX0V2BPn7vOBh1XVr6pbgM3A3JQFZzLCmyX822FF\nrLp0AgtnF5ET5Z34RnUH5zxTx7XrclldZ9OfGGNA0tUEIyLTgadV9XB3+Ubgi0Aj8C7wbVWtF5E7\ngbdV9UF3v3uAZ1X10dDyGhsb9wXu8/nS8RRMClW2C3dtzeaFut7Nht3mlQX4yrRODi6yZkNjhrOK\niop994uLi8OaY2J/QqTe3cBNgLp/fwl8KZmCQp9gonw+X7+OH6gyRnosFcC8I2B1XQc3rWxiaWXv\nZsFle7ws2+Plwml5fGfOKI4ak5OSWFJZhsUy+GMZbs9nsMXSXxnrPaiqNaoaVNUu4A/0NAFWAlND\ndp3irjMjwJyxOTx29liePm8sx4+LnpSe/ridU5/axYIXd9usycaMMBlLWiIyMWTxEqC7Z+FTwAIR\nyRWRGThfwpenOz6TWSdPyOXZ88fy2NljmF0U/cfHz21v5/Snd3HJ83W8Xu233obGjABpaR4UkYeA\necBYEdkB/AiYJyJzcJoHtwL/DqCq60RkMbAeCAALVdWGTBiBRIQzJucx9Sg/W/KncuvqJlbV9Z4Z\n+eWdfl7e6ef4cTl864gizp2aR5b07pVojBn60pK0VPXzUVbf08f+i4BFqYvIDCUiztxdZ0/J5cVK\nPz9b3cS7u3onr7drO3h76R4qir18Y3YRnzuwgLwoo3AYY4YuGxHDDBkiwllT8lhyQTlPnD2GE8dH\nv+blawxwzZsNHPHXan62uom6dquoGzNcWNIyQ46IcNrkPJ45v5xnzx/LmZNzo+63q72LW1Y1M3tx\nNYt8OWyo7107M8YMLZa0zJB2wvhcHj17LK9eVM7lM/PxRGkN9AfhbzVeTvhbLZ95vo4lO9ptfENj\nhihLWmZYOHJMDr8/tYzVl47n67MLGZUd/VrWSzv9XLZkN598opY/bGihubMrzZEaY/rDkpYZVqYW\nebl5bgnrLp/AornFTC2KPmylrzHAd95uZPYj1Vz/TgMfNQXSHKkxJhmWtMywNDoni4Wzi1j12fH8\naV4ZR46K3hmjqVO5e/1ejnmshsteqOOF7dZ0aMxglslhnIxJOW+WcPGMfGYH/DSVHMBv17fwxJY2\nos03uaTSz5JKPwcUefjXgwu5oqKACQU2wYAxg4nVtMyI8YnyHP5wahn/vHwC3zlqFGPzor/9t7UE\nuWllE7MXV3PF0t08v72dQJfVvowZDCxpmRFnYoGHHx4zmrWXTeDuT5VyzNjsqPsFFf7xcTufe3E3\nR/y1mru2ZrPFrn0Zk1GWtMyIlecVPn9QAS99ehxLLyzncwfmkxujNbCqtYv7dmRz9GM1XPjsLh7e\n3EprwHoeGpNulrSMwWk6/N0pZXzwuYnc+sliDiuJfbn39eoOvvZaPYc+XM03X6/nTRus15i0sY4Y\nxoQozc3ia7OK+PfDClmxq4MHNrXy+JY29kbpudHUqTzga+UBXysHFHm4bGYBn52Zz2Gl0ZsbjTH9\nZzUtY6IQEeaOy+V/Ty7lgwUTuOOkkpjd5sHpvPGL95s54W+1nPhEDb9Y08zmRhs2ypiBZjUtY/Zj\nVHYW/3pwISfITrrKp/NnXysPf9hKbVv0a1rrGwKsX9nET1c2MbvUy8XT87loej6HlFgNzJj+sqRl\nTAIOKcnmJ8cV89+fGM3LO/08tLmVf3zchj9GJWxdfYB19c0sWtVMRbGXTx+Qx5GSxYGqNueXMUlI\n1ySQ9wIXArWqeri7rgx4BJiOMwnk5apa7267HvgyEAS+parPpyNOY+LlzXKmSTlrSh6NHV38Y1sb\nj21pY9lOP8EYfTJ8jQFue78FyGP8pmrOnpLHOVPzmDcpl6Jsa6k3Jh7p+k/5E3BuxLrvA0tVtQJY\n6i4jIrOABcBs95i7RMSGJTCDVnFOFv9SUchjZ49l44IJ3HZCCadMzCWrj4pUTVsXD/haufKlPcz8\nSxWXPF/H3eta+LDRfgdmTF/SNXPxqyIyPWL1fGCee/9+YBnwPXf9w6rqB7aIyGZgLvBWOmI1pj/G\n5nn40qGFfOnQQuragzy9rZ2/b2vjlZ3+qENHAXR0wcs7/by808/1yxuZMcrDmVPyOHNyHidPyKHQ\namHG9FDVtNxwmgHXhiw3hNyX7mXgTuDKkG33AJdGltfQ0KDdN4V9N/9VV0VdP1K211588YCUv2LF\nikH5/Ib69ufO+pwW3LFWuX1d1O2/O/5S5fZ1Mbe32+uT8vf/ihUrBu3zG0nb9+2n4Z/9g6Ijhqqq\niOhAlNXY2Mg2nw+AY0fgdgDfAJUfTaaf31DffkyhnyWfbGNVY3K1p0c2t/DXf2zjhFLPoHx+md4O\n/X//d5cxGJ/fSNleUVERZatD0vVLfrd58OmQjhgbgXmqWiUiE4FlqnqI2wkDVb3F3e954EZVDWse\nbGxsHJDAffs5Qekqw2IZ/LGk4vmoKpubAryww88L29t5q8ZPRwKjQx1Zls1pk3KZNymX48fnku9N\nrEfiYD0vQ70Mi2XgFBcXh72pM1nTegq4CrjV/ftkyPq/iMhtwCSgAliekQiNSTERoaI4m4ribBbO\nLqKls4vXqvy8WOnnxR3tbGuJ/YNmgPf3dPL+nk5+vbaFXA/MLc/h1EnOtbCjx+aQ67Fu9WZ4SVeX\n94dwOl2MFZEdwI9wktViEfkysA24HEBV14nIYmA9EAAWqmrf/7nGDBNF2VmcNy2f86blo6p81BRk\naWU7Syvbeb26I+pwUt38QXituoPXqjsAyPPAseU5nDghlxPH53BcuXXqMENfunoPfj7GpjNi7L8I\nWJS6iIwZ/ESEA4u9HFhcxFdnFeEPKu/UdvBSZTvPftTExr19J6D2oDO47+tuEvMKHDUmm+PH5/LJ\ncTkcPz4nHU/DmAE1KDpiGGP2L9cjnDIxl1Mm5nJFcS2lU2fyyk4/y3b6eaXKz8f7aUoMKLxX18l7\ndZ38Zp2zbmpeHidV7eH48bnMHZfDoSVeG6nDDGqWtIwZosbmefjszAI+O7MAgK3NAV7Z6eeNaj9v\nVHdQ2br/VvXt7Vk8/GEbD3/YBsDobOHY8hyOG5fDJ8fl8InyHIpzrEnRDB6WtIwZJqaP8jL9EC9X\nHVKIqrKtJcib1X7equngzRo/HzbtP4k1dSov7fTz0k4/4PyA8tASL8eW53BsuZPEDi3x4u1ruA9j\nUsiSljHDkIg4SWyUl3+pKASgpjXI27UdvF3j553aDtbs7ow5TmI3BTY0BNjQEOABXysABV7hqDHZ\nHD02m6PH5FDSKjYAsEkbS1rGjBDjCzzMn57P/On5AOzt7OLvq7dQmTOOd2r9LK/toKFj/z9/bA0o\nb9V08FZNB7AXyGfUP6s4oiybOWOzObIshyPHZHNwsdXIzMCzpGXMCFWYncUnSrpYUDEKGEWXKpsb\nA7xT28GKXR2sqO3gg4YA8fyKv7lTebOmgzf3JTLI9cBhJdkcXhZyK82mJNeukZnkWdIyxgCQJcLB\nJdkcXJLNFw52mhQbO7pYuauD9+o6eXdXB+/t6mBXe3xDdviDsHp3J6t3h8/gPLnAw+wyL7NKs/fd\nDi62jyITH3unGGNiKs7J4rTJeZw2OQ9whp3avjfI6rpOVtV1sLKuk1W72mkKxN8MWNkapLI1yAs7\n/PvWeQWm5udxdOUeDivxclhpNoeWeJkxyovHmhhNCEtaxpi4iQjTirxMK/JykXttbNMmHzkTZ7Bm\ndyer6zpYu6eTf+7ppLot/kEUAwpbWrPYsqUtbH2uBw4a7SSxQ4q9HFLiJrPRXrItmY1IlrSMMf0i\nwr6eit2dPABq24Ks29PpJLF656+vMUBnAgMC+4Owrj7AuvrwyTG9AgeO9lJR7OWQEi8VxU4To9gc\nmsOeJS1jTEqMy/cwbrJnX9MiQEdQ8TUGWFffyfp9twA79iY2vGhAYWNjgI2NAZ7+OHRLAePXVHFQ\nsZeK0V7nb3E2B432csAoj/VmHAYsaRlj0ibHI8wuy2Z2WXbY+gZ/Fy+u3UJz4QTWN3TyQX0nHzQE\n4u70EaqmrYuatg7ecMdc7OYVmDHay4GjvRzkJrSZ7vKkgizEfmc2JFjSMsZkXEluFkeN7qLC/SF0\nt93tQTY0BNjUEOCDBieRbWpI7HpZt4CCrzGAr7F3G2K+R5gxysOM0V5Kg9l8IriXGaM9zBjlZUqh\nxzqDDCKWtIwxg9aYPA8nT/Bw8oTcsPUN/i58jQE2NnayuTHAJjcZfdTUSVATTzBtQWV9Q4D1DQEg\nmwcrG/Zty86CA4q8zBjlYfpop0fjjFEepo9ymhwLvPa7s3SypGWMGXJKcrM4bpwzsG+oDRt9ZE+Y\nzqbGAJsbA2xu6vlbm0TtDKCzC6ecpgBU+nttH5+f5SSwIg/Tuv8WedF2YXqXWi/HAWZJyxgzbHiz\n4KDibA4qzu61ramjiw9DkthHTQFnuSlAUxzDV8XSfQ3tndrILflkvbeTSQUephW5t1FephY696cW\neZlc6LHZpROU8aQlIluBZiAIBFT1WBEpAx4BpgNbgctVtT5TMRpjhr7ROVkcPTaHo8eG185UlXp/\nFx81B/moKcCKLbU05oxmS1OAj5qC7PYnV0MD6FLYsTfIjr1B3qzpvV1wampT3SQ2tdDDlCIPkwud\nW0enE591EumR8aTlOk1V60KWvw8sVdVbReT77vL3MhOaMWY4ExHK8jyU5Xk4tjyHo7s6qago27e9\nsaOLrc0BtjYH2dIUYGtzgC3NQbY0B6jcG9zvSPl9UaC6rYvqti5W7OqMskcB+e9WMakwi8mFTi/H\nKYVeJrlJbVKhh8kFWZTmjpzej4MlaUWaD8xz798PLMOSljEmA4pzsjhqTA5Hjem9rbNLqdwb3JfU\nPm4JsK05yLaWAFsaOtjd2f9E0hZUPmwK9jkfWp4HJhZ4mFjgJLJ99ws8TCzIor1dmBbUYdEUKar9\n+JowEAGIbAEacZoHf6eqvxeRBlUtcbcLUN+93K2xsXFf4D6fL50hG2NMXNqDUO0XdvqF6vYsdvqF\nqnahyi9U+4W6DkFJXyIp9ipjc5TyHOdv5G1MjjI2W8nzpC2kqCoqKvbdLy4uDjtBg6GmdbKqVorI\nOGCJiHwQulFVVUT6zKyhTzBRPp+vX8cPVBkWy+CPZbg9H4sldWV0l3PEoRUc0cc+HUFlZ2uQj1uC\nbG8JsL0lSOXentvHzZ20dQ1cUmsMCI0B4cPWvvcbnS2ML/BQnpfF+HwP5flZeFvrOXTKOMrzshiX\n72Gs+zffm97aW8aTlqpWun9rReQJYC5QIyITVbVKRCYCvfrlGGPMUJfj6ZlhGnJ7bd+0yce4Aw7c\nl8R2tjqdOna6y1Wtzt+9gYFtMWvqVJoaA/gaw6KFbQ299i3yCmPzsxiX52FsfhYnjM/hm4ePGtB4\nQmU0aYlIIZClqs3u/bOBnwBPAVcBt7p/n8xclMYYkxkizm/SSnKzeg191U1VaepUqlqDVO11pn2p\nbu2iqtVJblWtQSqb/ezuzKIrBVeDWgJKS3OQrc3ONbfcFP8uLdM1rfHAE26vFy/wF1V9TkRWAItF\n5MvANuDyDMZojDGDlohQnCMU52RxaEn0xObz+Zh54EHUtndR0+okspo2J7HVtAapbuuipi1Ibavz\ntz8Vt7H5qR0hJKNJS1U/Ao6Ksn43cEb6IzLGmOHJkyX7ehXO6WO/Lvd3a9WtXexqD1Lb5iS6TVW7\n6cwbza72LmrbuqhrD7KrratXgivPG8ZJyxhjzOCSJcKYPA9j8jxAT83Nl1cd9vs1cJomGzqUXW1B\ndrV3UdfexcHFqU0rlrSMMcYkRUQozRVKc7M4OE2PacMTG2OMGTIsaRljjBkyLGkZY4wZMixpGWOM\nGTIsaRljjBkyMj5gbrJCB8w1xhgzPEUOmGs1LWOMMUOGJS1jjDFDxpBtHjTGGDPyWE3LGGPMkDFi\nk5aI3CsitSKytp/llIjIoyLygYhsEJETkn18EblMRNaJSJeIHJvscxCRb7rxrBORn++njKki8rKI\nrHf3vybRWGKVkUQseSKyXETWuPv/OIlYopaRaCzu/h4RWSUiTycaR1/lJBnLVhH5p4isFpF3k4kn\nWhmJxhLt/Z7MeYn1fxNvLCJyiPs8um9NInJtEuckajlJnJf/dPdbKyIPue/DZM5Lr3KSiOUa9/h1\nIc8l0fPSq4xE40gZVR2RN+AU4BhgbT/LuR/4N/d+DlCS7OMDhwGHAMuAY5Ms4zTgRSDXXR63nzIm\nAse490cBm4BZicTSRxmJxiJAkXs/G3gHOD7BWGKVkVAs7j7/BfwFeDqZ16ePcpKJZSswNmJdou+X\naGUk+hr1er8nc15ilJPweXH38wDVwAHJvkZRyok7FmAysAXId5cXA1cn8frEKieRWA4H1gIFOGPL\nvggclEgsfZSR1Osz0LcRO2Cuqr4qItP7U4aIFOMkjqvdMjuAjmQfX1U3uOXG9fgxnsN/ALeqqt/d\np89Zn1W1Cqhy7zeLyAZgsqouiTeWWGUAX0kwFgVa3MVs96aJnJdYZZDgeRGRKcAFwCKcpJPw6xOr\nnERjiSWZeKKIO5Y+3u8NicQRqxwRSfa8nAF8qKrbQh4jzkOjlyMi/5NgLF4gX0Q6cT7sdyb5+vQq\nh8TeL4cB76hqq/vYrwCfUdWfJxBL1DKAYxOII2VGbPPgAJkB7ALuc5t//ijODMyZdDDwKRF5R0Re\nEZHj4j3QTYBH49ROkhJRRsKxiNOUthqoBZaoasKxxCgj0Vh+BXwX6Er08eMoJ5nXSIEXReQ9Eflq\nkrFEKyORWAbq/R6rnGTfuwuAh5KIo69y4o5FVSuBXwAf43x5a1TVFxJ98D7KSeS8rHX3HSMiBcD5\nwNQEQ4lVRtKfLQPJklb/eHGa5+5W1aOBvcD3MxsSXqAMp0nsOzgzQO/365WIFAGPAdeqalMyDxyl\njIRjUdWgqs4BpgBzReTwROOIUUbcsYjIhUCtqr6X6GPHWU4yr9HJ7nM6D1goIqckEVK0MhKJZaDe\n77HKSfi8iEgOcBHw1yTi6KucRN4vpcB8nGQ8CSgUkSuTiCFWOXHH4tbufga8ADwHrAaCicTRRxlJ\nfbYMNEta/bMD2BFSG3gU558xk3YAj6tjOc43/LF9HSAi2TjJ5s+q+ngyDxqjjIRj6aaqDcDLwLnJ\nxBOljERiOQm4SES2Ag8Dp4vIg0mEEKuchM+L+y28u0nmCWBuosHEKCORWAbq/R6rnGTeL+cBK1W1\nJok4+ionkVjOBLao6i5V7QQeB05MIoZY5SR0XlT1HlX9hKqeAtTjXGNOSIwykv5/HkiWtPpBVauB\n7SJyiLvqDGB9BkMC+BvOBVNE5GCci9x1sXZ2vyndA2xQ1duSecA+ykg0lnIRKXHv5wNnAR8kGEus\nMuKORVWvV9Upqjodp8noJVVN+JtzH+Ukel4KRWRU933gbJwmnLj1UUYi52VA3u99lJPQeXF9noFp\nGowsJ5FYPgaOF5EC93/hDGBDEjHEKifR98s49+80nGtRf0k0kBhlJPP6DDzNQO+PwXDDeYNWAZ04\n3yC+nGQ5c4B3gffdF7U02ccHLnHv+4Ea4PkkysgBHsT5QFoJnL6fMk7GudbxPk4zwGqcNuy4Y+mj\njERjORJY5ZazFvhvd30iscQqI6FYQsqbR0+vv4Renz7KSfS8zATWuLd1wA+TOC+xykg0ll7v92TO\nS4xyEo2lENgNFIesSyaWaOUkGsuPcb4crQUeAHKTjCVaOYnG8hrOl4A1wBnJnJcYZST1PzTQNxsR\nwxhjzJBhzYPGGGOGDEtaxhhjhgxLWsYYY4YMS1rGGGOGDEtaxhhjhgxLWsYYY4YMS1rGmAElzhQo\nZ2Y6DjM8WdIyxhgzZFjSMsOS+23/OhF5X0QaReQRcSbmu1pEXo/YV0XkIPf+n0TkLhF5VkRaROQN\nEZkgIr8SkXpxJsA7Oo7Hnyoij4vILhHZLSJ3uuuzROQGEdkmzgSe/yfOVB2IyHQ3li+KyHb38b4m\nIse5z6Ohuxx3/6vd+O50n+MHInJGyPZJIvKUiOwRkc0i8pWQbTeKyGL38ZvFmdTv2IhjH3Pj3yIi\n34rnWBF5AJgG/N09f991z/uD7nloEJEVIjI+8VfVGEtaZni7HGew3Bk4wztdncBxN+AMBuoH3sIZ\ntmYszuCufY7RKCIe4GlgGzAdZ26xh93NV9Mzsd9MoAi4M6KITwIVwOdwpjb5Ic5gqrOBy0Xk1Ih9\nP3Rj+xHwuIiUudsexhm6ZxJwKXCziJwecuxF7j4lwFPdcYhIFvB3nCF8JuOMgXetiJyzv2NV9Qs4\nY+h9WlWL1JnH6SqgGGd6izHA14C2Pk6hMTFZ0jLD2R2qulNV9+B8CM+J87gnVPU9VW3HGQ29XVX/\nT1WDwCM484X1ZS5OoviOqu5V1XZV7a7dXQHcpqofqWoLcD2wQERCJ2S9yT3mBZxpOx5S1Vp1Rml/\nLeLxa4FfqWqnqj4CbAQuEJGpOKPMf88tazXwR+BfQ459XVWfcZ/XA8BR7vrjgHJV/YmqdqjqR8Af\ncAb93d+x0XTiJKuD1Jk25j1NcvobY0bszMVmRKgOud+Kk0jiETrNRVuU5aL9HD8V2KaqgSjbJuHU\nwLptw/k/DG0uS+TxKzV8ANFt7mNMAvaoanPEtmNDliPPT56bPA8AJolIQ8h2D07C7PPYGM/5AZxz\n8rA7Av+DOAP2dkbZ15g+WU3LjDR7caYxB0BEJqTgMbYD0yJqT9124iSFbtOAAOGJKRGT3WksQsvb\n6d7KuqcjCdlWGUeZ23HmdSoJuY1S1fPjjClsFG63FvhjVZ2FMz/UhYTX+IyJmyUtM9KsAWaLyBwR\nyQNuTMFjLMeZMuZWdx6rPBE5yd32EPCfIjJDnJmebwYeiVFDicc44Fsiki0ilwGHAc+o6nbgTeAW\n9/GPxJm6Jp7JLJcDzSLyPRHJFxGPiBwu8U+vXoNzvQ4AETlNRI5wr/U14TQXdsX/FI3pYUnLjCiq\nugn4CfAi4ANe7/uIpB4jCHwaOAinU8IOnE4VAPfiNJe9CmwB2oFv9uPh3sHptFEHLAIuVdXd7rbP\n43QE2Ylzbe5HqvpinPFfiHMNcItb9h9xOlPE4xbgBren4HXABJwOLE04kxq+gnMOjEmYzadlzBAl\nIlcD/6aqJ2c6FmPSxWpaxhhjhgzrPWhMEkRkGs505NHMUtWP0xmPMSOFNQ8aY4wZMqx50BhjzJBh\nScsYY8yQYUnLGGPMkGFJyxhjzJBhScsYY8yQYUnLGGPMkPH/ASeof3mb5vtVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f33b7646048>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(loss[:,0], loss[:,1]);\n",
    "ax.axhline(100, linestyle='--', color='r', linewidth=2)\n",
    "ax.xaxis.set_ticks(np.arange(1, 100, 5));\n",
    "ax.set(xlabel='num_components', ylabel='MSE', title='MSE vs number of principal components');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But _numbers don't tell us everything_! Just what does it mean _qualitatively_ for the loss to decrease from around\n",
    "$450.0$ to less than $100.0$?\n",
    "\n",
    "Let's find out! In the next cell, we draw the original eight as the leftmost image. Then we show the reconstruction of the image on the right, in descending number of principal components used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQMAAACaCAYAAAAdFl52AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XusndV95vFnlVsAE5MANsYYbMA2mEsI2JCGQEqYISEd\nNZeWKqitmJk2jVonymiaajqRqrSKKkWjmU5RapFp0miolGmD1EShbZoxgaRoaLk4Lthgg83FYIPB\nmIu5x1ze+cPHznmf9dj79fY5Z+/j/f1IUbw26+z97v1711rvfu3zrNI0jQAAAAAAAAAc/H5u0AcA\nAAAAAAAAYGpwMxAAAAAAAAAYEdwMBAAAAAAAAEYENwMBAAAAAACAEcHNQAAAAAAAAGBEcDMQAAAA\nAAAAGBGHHsgPl1I+Iuk6SYdI+kbTNF/Z/d927NjRHOCxAQAAAAAAAOjTzJkziz/W978MLKUcImmF\npKskLZF0TSllSf+HBwAAAAAAAGAyHcivCV8k6aGmaR5pmmanpL+R9LGJOSwAAAAAAAAAE+1Afk14\nrqTN49pbJF2cOm7cuPEAXgbT0dINyyRJqxbdPeAjAQAAAAAAGB0LFy7c538/oMzArpYvX64VK1Zo\n+fLlU/FyGAJ3fn7X/1PzwWPsAVPr7bff3vPn66+/Xr/zO7+jn/s59us62I2vuyRqPgQme/3zmkvU\nHdiN609gMBh72G3lypX7/O8HcsXyhKR549onjz0GAAAAAAAAYAgdyM3AuyUtLKUsKKUcLulTkm6a\nmMMCAAAAAAAAMNH6/jXhpmneLKV8VtL/lXSIpG82TXP/hB0ZAAAAAAAAgAl1QJmBTdN8X9L3J+hY\nAAAAAAAAAEyiKdlABACAUeEbCLChwGigzqOHmgMAgOmKqxgAAAAAAABgRHAzEAAAAAAAABgR3AwE\nAAAAAAAARgSZgUAHpZQJ6dOPpmkmpA8m3kTVPD0PdR9e+1v3UkpVqy7PQX2HRz/ZcBNVP86DA/f2\n22+32ukzfeutt1rtLjXvMo4POeSQ/f45ar7/vMbuzTffrB7rZ15O+smJJW9y4qVzwB/zcS7VdU/j\nz/uk+lH3wfAadzkPJmp+79KHmk88aj5xps+RAgAAAAAAADgg3AwEAAAAAAAARgQ3AwEAAAAAAIAR\nwc1AAAAAAAAAYESwgQhGXgoC9ceOOOKIfbYl6cgjj9xnW5IOP/zwnq+9c+fOVvv111+v+vhjr732\n2j6f59BDD42hyaMSUt4rWFzq77M49NB6Cj3qqKNa7WOPPbbq44+l88nDzl9++eWqz44dO1rtV199\nterz05/+tNXu8lkcDLqECSc+JlONfWzPnDmz1X73u9/dai9ZskTHHHPMPl9Hqmvl9ZWkF154odV+\n5ZVXqj4+hxysNe/yvtKmAb2eJ9XGx+g73/nOqs9xxx3Xaqex72tAOj4f6+k8eOmll1rttAb4nHYw\nzPep5v6+3njjjVb79ddfrz5nHyNSPf5SAPiMGTNa7Tlz5lR9TjzxxH3+jFRvKuJzt49zSXrxxRdb\n7VRzX+enY837Hdd+XdSrLdXnQaq5r+knnHBC1ccf8/leqtcSP0+luqZpfvf3kWo8HevuumwQ4GPW\n21I9ttJ54HP+0UcfXfWZNWtWqz179uyqj9c9XT943VONfQ04WGvez1jvMnenseVzY9rYyevnNfc1\nXqrnh/S8Xr90re41Z37/mVRzfyzV3J+nn5pLdd295um50zzjdU/f67zuU/HdnX8ZCAAAAAAAAIwI\nbgYCAAAAAAAAI4KbgQAAAAAAAMCIIDNwisyfP7967PTTT2+1L7/88qrPO97xjv1+Lf8d9Jtvvrnq\nc/vtt+/3805HXfIAU4aA5zp5poBngknd8qI8hyQdn2cKbN++verz/PPP77MttfOFjjjiiJilcjDk\nT7gumVJdMhgOO+ywqo+Px3e9611VH8+Qmjt3btXHMylS/oTnYaTzYMuWLa32k08+WfV57rnnWu0u\n2TPDrt8ad8mG8zzAVGPPBPP5fd68ea32pZdequOPP771WMoS8hy4VM9Nmza12o8//njVp1fNpYOj\n7t5ONfccmXTueC265EWlNf20005rtdPY9/MrZeE8/fTTrfZjjz1W9fHzYNu2bVUfz545GGqe3kOv\nfKGdO3dW50E6V7zuKQ9wyZIlrfaZZ55Z9fFzI10L+DE+8cQTrfbDDz9c/cxDDz3Uaj/11FNVH88f\nm44171LjlKfVKxsuzfee/XnSSSdVfRYvXtxqL1q0qOpz8sknt9opJ9Lfw9atW6s+Pp+n+f2ZZ55p\ntbvkhQ67NC/759Ul87HLZ5HyXn1NP/vss6s+Z511Vqvt871U1z3lmG3evLnVfvTRR6s+veb3g6Hm\nUrc13Ouevst0yf7skvfqNT7nnHNa7VRzP5/S8fm1utdXkh555JFW29cEqb5GpOY/0yUP0Me51Lvm\nUl33NIf48fg4l+q6e82luu5ec2ni635ANwNLKZskvSTpLUlvNk2zdCIOCgAAAAAAAMDEm4h/GXh5\n0zT1P1kBAAAAAAAAMFTIDAQAAAAAAABGRDmQ3zsupTwqaYd2/Zrw/2qa5i92/7cdO3bseeKNGzce\nyDFiGlq6YZkkadWiuwd8JAAAAAAAAKNj4cKFe/48c+bMKjz3QG8Gzm2a5olSyixJN0v6XNM0t0nt\nm4FXXnmlVqxYoeXLl/f9WsPMNwKRpA9+8IOtdgocTmGXU+Xuu9s36W688cYJff47P3+XJOni6y6a\n0OfdXykw2sPi0yYtHg56wgkntNqzZ8+ufsY3CEgBo0cccUSrnYJTPYT0hRdeqPp4mHAKj98dMv2l\nL31Jf/zHf9zaUGQ3D2Ud9jDaFDLdhQdRp3Bh38Qjhb57SPipp55a9fGNBtJmBL6JQOrj80MKkfWN\nBTZs2FD18YDaZ599ttXu9zOdSl02EfCx5DWX6rr7pj9S781BJGnBggWttm8oND6o+rjjjtOzzz5b\nnV8ebCzV81XaGMj/cu2+++6r+vjmA15zafjr3mWTGK9xCpn2n0ljzTeJOeWUU6o+HiDt8316nrQG\n+LyS1h/fBCKFja9du7bVfvDBB6s+vtFAGhPDpMs5md6DPzb+eW644QZde+211YZQM2fOrJ5n/AWz\nlDcR8E1h0hzidU/ngZ+HXvP169dXP7N69epWO41933xm2Greb409UD5toOBrgI+tLuHxvkGMVG80\nkOrp83uaZ/x40nvwcbxu3bqqj58b6fpvWOq+t+9+vTaDkvJGQM7fZ9r4zceobwgjSeeee26r7dd6\nUrf53TepSOeBX8ulGt9zzz2t9nSq+W79XLelGvtj6Vzxa2r/zibVYz2dBz7WfZ3osqanazvfSCb9\nI6k1a9a02n4OSPU1f/q8dhvEfZd+NvfrUvN07njN0zWZ1zxt+OXrQro28MdSH6972jzINwG79957\nqz5e97SB3L7qnqxcuXLPn9PNwAP6NeGmaZ4Y+/9tkr4rabB3fgAAAAAAAADsVd83A0spR5dSjtn9\nZ0lXSqr/mhIAAAAAAADAUDiQ3YRnS/ru2K81HSrp/zRN84MJOSoAAAAAAAAAE67vm4FN0zwi6T0T\neCzThv8O+q//+q9XfQ4//PBWO2XXdcln87w4/33zlEfhx5csW7as1X7iiSeqPrfffnvP5xk2/jmn\nXEavTcp38OyBk046qdVOeRT+PJ5NKNU1T9l1/nOeiyDVeQUpP2B8xs4xxxwT87S65DIcDLrU3LP+\nUs6n54Z5VpxU50Km3Aj/3D0vKh2j59VI9fmdMmM8K/Lll19utVNu0SD1mynlOSRpbHl2UMp79cyo\nlAvp+VB+zK+++mrrNV999dXqc9+xY0f1vD7vpNxKPwdT/Xzd8Nfe288NUpe6e43T2HI+jrrkhnkm\npFTPuWlN91yzNK49BzKtJT7WU26lz+cpX9LPsWHPlOq3j8+Dns121FFHVZ/zeeedVz3POeec02qn\nc8UzydLa69dTTz31VNXHc4d9LfFxLrXnFUnavn171cfH/qBr3k+Nu2TDpWs7Hzfz5s1rtS+6qE4y\n8lzINB59LfFxLtWf+3PPPdfz+NL87tcd6ZrMXyutJYOu+3j95L+mn0trun+mfh0nSRdccEGr7eNc\nqrPi0vnlx+j5nFKd5ddlfk/Zdf5afh03bDVPNe6SC+lreMpR91qk624f60uXLq36+Dqfct/89f3a\nKdXcz0s/l6Q6g9LzaaX6u1/KCfd5JeVBT6UuNfbHutTcP4uU1ejZzhdeeGHVxzMC05zbq+ZSXfc0\nF3ndU+6o1z3dJ/C6p30DPA/6QB1QZiAAAAAAAACA6YObgQAAAAAAAMCI4GYgAAAAAAAAMCK4GQgA\nAAAAAACMiAPZTfig5MHwKWT6F3/xF1tt35ygX7fcckv12L/8y7+02h4a62HWkvTRj3601f7ABz7Q\n87V904PpykPd0+fjm674pgJSvWGIhxKnjVtcCiH10PcuYbmJB6TPmDGj6jP+sRkzZnTaRCAFd083\nKZDVg4I9eFaqQ5zT5hIeXJxCZD3kPQWJe/B7+tz9nEsbWfjxpMBaD7PfunVrqz1sG0kkHoydxo3X\nIoVMe5jw+eefX/XxzRrS/O4B8l7jJ598cs+f582bp40bN1ZBz2mcn3HGGa12CpnuFZIvSZs2bWq1\nvebS8Nc9haF78LQHU6d50OvpGwZI9VhPYdV+zqXg7s2bN7faaeMIX6O6bGLTJZDca56Ox+emYZM2\nS+iyAYWvh75+z5s3rxr75557bvU8/pmmse8btTz66KNVn0ceeaTVThvJ+AYhPhf5e5CkuXPnttrp\nvPDzIL32IKUa+/qX5neXguB9De+ycYRv5JKuH7zm4+f33fxzT6HvvplE2jjCr03Stcpjjz3Waj/+\n+ONVn2Gru/Map3Hua3ra+M2vi9L87o+l6ySfQ9IGHT6f+jiX6nMljWM/D1ONfd3ymntbGr6ae027\nbPyWNubysZ429/P5M411H39pnvENYPxzfuCBB6qf8XXe52lJuvjii1ttn5ukuuZ+vknSxo0bW+1B\nbyDi0jjudf0+6Jr75iDpWsrrnjaS8bp7zaV6k5P0vc6vDbzmEhuIAAAAAAAAAOgTNwMBAAAAAACA\nEcHNQAAAAAAAAGBEkBloPFfml3/5lyftte69995We+XKlVWfXnk5nqEkSXfccUer3SUz8P3vf3/1\n2K233trz5wYp5bV53suRRx5Z9fEssRNPPLHq4zkynhGYcm88Y8SzQ6Q6RzA9j79Wyq/yvLGUbTT+\nvR955JGxT/oMh1nKWfNcyJQl5LlqZ511VtVnwYIFrbbnDEr1eEx5HZ4fknI//OdSlopnDqVz2Y/5\n+OOPr/p4JuZ0zAf1zz3lOvk48blcqvM6PKNPqsdEyn7yzKh9ZTh96EMf0po1a2J2pPP5wTNs0zGn\nzFPPSTkYai7V2TKeEehZK1KdDZdyunxseSakVNc85co8/PDDrXbKlfGxns4Lz69Kc5o/5jWXpl/d\nu2SBpjXL13Af+2eeeWZ1HqQsLz+/PD9Kqut+//33V338PHjxxRerPp754+t+ykfzx9J87+fOsEk1\n9uuglCHlYz3lqS5durTV9qy4NEb8tVPNPXf3oYceqvr4YynTya8p0vv0Pp4RK9VZkcNe83Sd69L1\nqV/bpWw9z1dN87tfz6fP3a/JvOaS9OCDD7baGzZsqPr4nJ8yRf3101j3unse2bDXPElrul/Tp/fl\nubrvec97qj6eF5fmRs8mTdd2W7ZsabW9fmnse590Xvi8l9Z0z8BLeXJd5pAu420qed37qXnK9/b5\nPV0Le05/qrmP9S7Xdmlce93TWud175Jbnr6Let0PtObT604AAAAAAAAAgL5xMxAAAAAAAAAYEdwM\nBAAAAAAAAEYENwMBAAAAAACAEdFzA5FSyjcl/TtJ25qmOWfssXdL+rak+ZI2SfrVpmnq3RKGnIdP\nStLHP/7x/X6eFIrqIaQpANJDK3ttFjKZ0vENGw/MTEHiHpieNt+YNWvWPttSHebtQaAphNSDp7dv\n3171eeWVV1rtFJrsIbepj4dpp6D48Y8dccQRVRizNPwbiHjN03vwQFbfLESqNwzxjTekuuZeK6kO\nBU8bBPi4TkHivtlMCpr1jQZS4L2HSvt7kOqxPeybCqTNVPw89fNfqmuagnl98410PnUJAPdNYXxz\nCa/5k08+2anmXhvfjEaqQ9TT/OCB5Ol9Dhuve5qbPGjaz38PEZekRYsWtdppI56XXnqp1fb1W2pv\nCpPaUn3upDXANydJmyL566dQfA9IT2N/mELm0/VNGuvOz920CYR/Pl7zRYsWVRuHpdf2NTzV2IPE\n07ni49bHfuKh5WmDDH8Pqb4+hwxbwHz63P140vvy+T1tIuCh8z4mfIMmqb6WSxtH+GM+30t1zdO1\ngb9WOpd9jUobDfj8njaamuiA+f3hYz2tdT7vpesS3yjFx7VUr4fpu0yXjSO8Xmns+1jfunVr1ccf\nSxtE+YYmft5K9QYi/l0mXQMNsuZJl/nd1+N0be6bxHS5fk9j3efhNHf7ZnB+LZfmcq9x2oDM65mu\nVfx7StpYxueDQc/v/lr91Dxt/Ob3Z3xTDakeA6nmvnlXqrmvvV2+s6Vx7XVPmz/5uZw2s/T3ldaA\nQWwg8r8lfcQe+wNJtzRNs1DSLWNtAAAAAAAAAEOs583Apmluk+S3QD8m6YaxP98gaf//OR0AAAAA\nAACAKVW6/NPCUsp8SX8/7teEX2ia5tixPxdJz+9u77Zjx449T7xx48YJPGRMB0s3LJMkrVp094CP\nBAAAAAAAYHSMjx2ZOXNm9bvlPTMDe2mapiml7POO4vLly7VixQotX778QF/ugPjvan/qU5+q+qQ8\nIed5GCk/5Lrrrmu1UzbIs88+u8/n7ZfnUXzhC1/o63l+//d/v+9juPPzu/5/ImvuGVIpY8R/tz7l\nrHn2QMoMdJ4PkPJgPCskZUu88cYbrXbK+/JcoJNPPrnqc9xxx7XaKadh9/n1yU9+Ut/5zndiBoqf\nuyknb6qk7AvPlUkZkF7PlCvjOTwpY8ezJVJem39eKTPGz5VXX3216tPlc/ZMsGXLllV9LrroolY7\nZVTce++9rfbNN9/caj/yyCM9j2Uqeb6PVI/1U089terjn8+FF15Y9fFcmZQfsmHDhlY7fT4+/j1z\n7rXXXtvz569+9av63Oc+V/VJY/bEE09sta+44oqqzyWXXNJqp3N59erVrfaPf/zjqo9n40ylLvlx\nhx5aX574fO65Muedd171Mz5npBw/z4zxtlSP9ZQ7Nb7uUl4DPFcmZVx98IMfbLWvvPLKqo+vC+vW\nrav6+Fj3vLuplGru62HKtvS1zq/jJOncc89ttcePoxNOOEHPPPNMdX6lzB9fI1NeqOcJpefx7MGU\nZeTXKpdffnmr/eEPf7j6GZ8f7rvvvqrPLbfc0mqn+WuQ+XEpT8vXeX+fUj2/p0xYXzN9nfW6SPU6\nnzIDfaw//3wdke5rQurj7zOtUVdddVWrna5nPLP2+9//ftXH54NB1fz666/XZz7zmapPlyxQz9NK\nn4XXPL1PPw/SmPVrgZQb5o+l678u3+suvvjiVvujH/1o1cfnNL/2/Lu/+7vqZ9auXdtqp6zGycoJ\nT6/l62Ga332sL126tOrj+aDpO9tPf/rTVjutzz7W0zrvNfZ6ppr3yoyW6muVX/mVX6n6/NIv/VKr\nnT7Tf/iHf2i1f/CDH1R9dp9zU3HfxY/Ray7Vdfcs0DQPdqm5ryVdskBTzX1dSPOD1z3d9/G6p2uV\nq6++utX2mkv1Z+o1l+q697p/tHLlyn3+935nhadLKXMkaez/6xUWAAAAAAAAwFDp92bgTZKuHfvz\ntZK+NzGHAwAAAAAAAGCy9LwZWEr5a0n/ImlxKWVLKeU3JX1F0r8tpWyU9G/G2gAAAAAAAACGWM/M\nwKZprtnLf6rDjAAAAAAAAAAMrQPeQGRYLViwoHrsmmva9zXf8Y539HyeFMr4ox/9qNVOAZ4uBYpO\nlrSxRi8p1HnYePhyqp+HsadQYg/uToG6HkTqIaQpZNpDR9Nn6uHGKezYg87TRgMubb4x/rmnMjy6\nX+k9+LmcgsVPP/30VjuNfd9kIQXEeii3t6U6gNg3hZDqDUNSCLAH56f37j/nPyPV51jafMNrn15r\nmKSNI3xTnRQk7mG9aWMGr59vFiLVwespMNo/Z2+nMdxlDHpt0rnT5bzwx6bD+Pe6ezC81N4RTZLO\nPPPMVjvV3IP80yYafh6kmntAefrcXTqXfR1Loep+3ZHGtc8z02EN7+WYY46pHjvttNNaba+5VAeS\ne23efPPNqqYbN26snsc31Unzu28Gkmrj1xTpGsPPA9dlvk9h7WnOGHa+8VWq8eLFi1tt31hGquvl\na3jaTMXXhDSOfAOKtCmYz/mpDl7zNC/7PJNq7I8Ne83T+e9zdbpu88f8OkCqx1/auMVrnDYA9HMn\nbUbgdU/fD30NT+Pcz5UuNR72mqfj8c9nxowZVR/f3M/XeKn+HpfmXK9x2hyuyyZgvmGIf+7pXPY1\nPG1E6mO9yzyTPlN/nnQ8E7UZaRd+jOm1ve4TVXMfx/3W3L8PpvHon3O6buuyAa3XPW0m2eUz9eM5\n0JpPzrZCAAAAAAAAAIYONwMBAAAAAACAEcHNQAAAAAAAAGBEHLSZgZdddln1WJeMQM+bSHmAq1ev\n7v/AJljK0Pi1X/u1nj/nmRVf+9rXJuyYJkLKNfPf0U/5E54H6O30PCm7wXMjPJsgZc51yXvw99Xl\n9/xT5og/lp5n/GNvv/12p/yJQUrZF16/U045pepz6qmnttopN+zFF19stR9//PGqz6OPPtpqP/nk\nkz2fJ+U6eS1Sbbrk+Pl8leYvPw/Suez5F12yzqaSv/ejjz666uN175IxkvK+PCfsgQceqPr4uZHy\nQ/yY+xlHaVx7Zto73/nOqs/hhx/eaqea+3tP72GQ0lj38/vkk0+u+syfP7/V9vXP89ykOics5UR6\nzdNn2iVrs0sfX3/SfOWPpXwaX2927NhR9Rn2HEF/XykT1sf6vHnzqj4+x47PCZozZ462bt1ajXWf\n76U6Tyjl9fpjXca+j1mprrGP9TTf+zhOeUPDNtZdel9z585ttT0nUpJmz57daqex7uPY5/uUKeXP\n0yWDNeVX+XmR5ne/ZvWsRKnON/YMQame31OG4SD5/J5yzGfNmtVq+3Vc6pNq43nsnh0n1ddyKTfM\nx02aO7tk8XpOrNdTqq9rU9aY193n9y7jPK2zU8k/C6+nVI91nwuk+rtUyvT1sZ9y3f37vV/PS/VY\n8nHd5btp+t7pGadp3ff3mfIv/TyYynzALlJOstd9utVcquuersm87inX1uuevpf7Madru4muO/8y\nEAAAAAAAABgR3AwEAAAAAAAARgQ3AwEAAAAAAIARwc1AAAAAAAAAYEQcNBuIePhsCpjvYtWqVa32\nMG0Wklx44YXVY2lTEbdu3bpWOwVvDlIKafWQXQ/bl+pA5hTe68GbKUzYQ4k90DOF93roaAr49HDV\nFDLtYeMpfNl/LoVMjw+53rlzZwxFHabw2RSyftJJJ7XaaQMRD6dN4d5+fqfweA+jTQHlHiDda+OW\nvfFg57Rphp/LKWzcz6e0aYYH5g5bwLyfy74RiFRvHJE2EXCbN2+uHuuyiYDXPc1FadyOlzaR8edJ\noeEeOJwCtz2AP21m5I+ljQYGKQWb+wYBKVDeNxXxAPe06Y9vIrBp06aqjwdIp+PzmnfZLCSFaXv9\nUo19I4206YKP6+3bt1d90nwwKGle9GsVH+eSdPrpp7fa6bPwsT6+5hdccIE2btyohx56qNUnjRsP\n804bBHifLhuFpc3OfKx7O62Hfh3i1ylSXreGSZrfPVA+rfP72iRmt15jPY0Hr1/apML7dNlYJq0R\nvob7nCfV17Wpnr5JRgrFHybpO4nP72nDKJ+H0/num8Kk86LL5hu+YUiXTWLSeeBzftpMYs6cOa12\n2ijM6+6bGKZNBQYpfRb+3tPGQAsWLGi105rp53vaANAfS2PCr4PSxjs+/v19pfnB5+p0Pe/nu58D\n6XjShkdpnR8k/3zS+e5193V+2Gue+qT12euermG97ul4vO5TUXP+ZSAAAAAAAAAwIrgZCAAAAAAA\nAIwIbgYCAAAAAAAAI6JnZmAp5ZuS/p2kbU3TnDP22B9J+rSk3QEOX2ya5vuTdZBd/PZv/3arnXLW\nnOevSNJPfvKTCTumqZAyGLpkGT3yyCOTcTgTpkvWUsrh8d/ZT7ktnheSzgN/zHMHUm5El6w4P56U\nf+TvK+UeepZKygwc/z5fe+212CflHQ1KylCbO3fuPttS/Rlu27at6uOZUil/wmvebx6gS/ljPj8d\nf/zxVR9/ryl3x48n5WD555EyKgbJP4tUY8/eSJksnv2Z8gD9sZS74/Nnmotcygj09mGHHdZ6LNXc\ns1Q8Oy5JOUqet5Jyp9J5OVXSvOd1T3krXnd/n6nmvtal3DCfl9O64edFypNzKVfGs+FSTp5npSZe\n95SXOEz5cX7+S3V+zqJFi6o+Pu+l9drn8w0bNlRt/7xSLpCPiTTfe93TeeDvNc3dnovnWXqeYSbV\nuWFprUvXJoPkYz2Na8+FTPO71y+N9YcffrjV9gypVHOvX7om8j7peVy6butyPePzTFrTfawPWyas\nn/9pPvPvLqnmnpWV5jg/D9Jn4WMpjRHPCEzjusv1n5/v6b37Y2md8PPdc6+HKQ92b/x9pvnd8xJT\n7pvX/bHHHqv6+LVcWvv8POiSo96rLXU7388444xWO2Xe+zVsyr1+9tlnq8eGSXrvXvfpVvP0WLqe\n8fncay7VdU97Fvj1TFoDJlqXbwP/W9JHwuP/s2ma88f+N9AbgQAAAAAAAAB663kzsGma2yRN/m1J\nAAAAAAAAAJOqdPlVl1LKfEl/b78m/B8k7ZC0StLvNU3T+t2NHTt27HnijRs3TtgBY3pYumGZJGnV\norsHfCQAAAAAAACjY+HChXv+PHPmzCpLrt+bgbMlbZfUSPqypDlN0/zH8T8z/mbglVdeqRUrVmj5\n8uV9vYn4NkhnAAAeuElEQVQu/uRP/qTV7jcz8Gtf+1qrPey/n/+Zz3ymemx80ffmu9/9bqt9++23\nT9gxSdKdn79LknTxdRf19fMpQ8MzdVLeiufwvOtd76r6eIaAZ3FI0qZNm1ptz+9IeT5d8vf8vOyS\nJeRtqc5XSefp7tyB3/qt39I3vvGNmLfiWQSDzBvyLC1JuvTSS1vt9773vVUfr3HKA1y9enWrff/9\n91d9BpkZOGvWrKrPkiVL9tmW6vyJrVu3Vn3uueeeVnvNmjWt9qDzhvxcPu+886o+F198cavtGVNS\nnR9yxx13VH18nvPMOalbZmCvbLHxuSRf//rX9elPf7p6jpQH+L73va/VvuSSS6o+XvN777236vPj\nH/+41U5/ITfIzMCUp3X++ee32hdccEHVxz8zr1/K/PXH0rrvGYEpD6ZLZqCvAWm+8HnO37ckLV26\ntNX2jB2pzke78847qz6ep9XP/DVRUpaxZwl96EMfqvp47k6qn4+Bu+/+2V9C/uEf/qG+/OUvV/l6\nE5UZmJ7Hz4O0zp977rmtttc8/Yxfl6xatarq4zVPa3qX6/6J4tdy55xzTtXH572UO+XXYPfdd1/V\nxx+brMxAz5eT6mvClIvqa5tf30j1HJfWdB/rDzzwQNXH3+tU1nz8/Plnf/Zn+sY3vlH1ueii9neD\ndA3kmYGeBSrV8+BkZgb6z6Ucbl87zjzzzKqPr+spR9PP97vuuqvVTtewfnxTWfOURe3XrFdccUXV\nx6/lUv0eeuihVnv9+vVVn37y49I49vHv4yh99/PPOc1fH/jAB1rtNA/6+eTX7pK0du3aVjuth7vX\nrcm+7yLVdU/fU7zu063m6XnS2PJ7Eun63eue5pB//dd/bbXTWtfl++p4K1eu3PPndDOwdzJ60DTN\nnhTjUsrXJf19P8/Tr3nz5lWP+UVcl0nwW9/6VvXYsN388y8O73//+1vttIFIlwuXdIExTLpsIJIC\nh/0mSzoPfPJKAZ4+qfgi2+XGX3oP/iX42GOPrfr4za30Pv310yQ4foJ95ZVX+r6BOVXS+/QLxHRz\n1yfBNIY9iLdL8Hqqn3+hTV9w/efSxij+RS9tIuCPpQ1zPEQ6jWsPnU+LzyB53X1TASlvtuH8ffqX\nBKm+eZTGRKqX67Xw+jlw6KGHVjfx0g3NxYsXt9rpfPcvSClk2s/vtAb4l/SpvDmYxrrXPX05dP4+\nu9Q88bHVZROuxG8qzpw5s+rjN7dSqLqvC2mjG/9LrLSRjJ/fXW5sT5Y0rvwvutJf8Pkxp7n7wQcf\nbLX9L4Qef/zx6rPosjFXF2mzGT+/01/o+bWbn4NpHduyZUur7XOBNNgbvol/FunGh98AS9clfr6n\nv+DwNcCfJ93k78I/03Rd2WUjLP+Lev8Lbqm+Dkl/ietzWlrTB/mXPT7WTz755KqPz/fpM/WxnjY+\n9M8inTs+n/e7MVCXTQT8XPY1Xarfe7r54O/Lr+3StYubynMgrS0+76VrO58/05zmf8HRZa1L1zz9\n3Cz1mqfP1K/T0o0+n+/T83iN0z9sGLaNY/w6Mq11fnN02Gue+ni90rX52Wef3Wqn+zP+POk7m9d9\nKmre10xRShk/oj8hqb5tCQAAAAAAAGCo9PyXgaWUv5b0C5KOL6VskfQlSb9QSjlfu35NeJOk+ndV\nAQAAAAAAAAyVnjcDm6a5Jjz8l5NwLAAAAAAAAAAmUV+ZgYOWMnZSXoHzkN30+/jDxkOlL7vssp4/\n479Hn7IRU/bGMEk5PJ4ZmDYZ8d/HT1kcnhmYgm/957rkRvhjKePNMyDTJgLeJ53bniGVcg/H5wy8\n9NJLMVdmKgOGe+mSrZdynvy9p3wt3yjltddeq/p0yVfpstGAZ8OlvDvPPU35cX5upDwMz1F67LHH\nqj6ex+HPM8hsIamue8pR8s80fRaeKeL5UVJ9bqSx1SVzyz8zPw98w4cTTjihqvl73vOe6nk9TysF\nnfvcnXLyPFx4mLJBpTyOZ8+e3WqnTTM8R81rnPK1fBOB9NpdMsE8dyqNfT/mlBvmAdupj9fLN46Q\n6kwdn+Okeh3rcp00WdL87lliXTb8Shk7Hjbuc8EzzzxTZbp1+SxSdqSvAel9+bmcNnXzec7fZ7pG\n82vWdO0ybJmBPt5SuL6PmzR3e93T5+NjPW3i4frJB01j369VzjrrrKqP5wCn61xfw7vM76nmw5QZ\nmLLiPBc1ZYH6tUuaB31NTzXvZ95La4A/T5qvvO4LFiyo+vhclL6L+pzma1/6bjPI+T1t2ul1Txnp\nPu+l88DHRNo0w6X6ddlUx8eSzw9po0PfJCZtGuP5wSnL2Gue5kH/vNIcMpV8Lkxj3d/7VNa8y6Zu\nvWou1XVPNfbHUma0191rLtWfR8oHnei6D/ZbIAAAAAAAAIApw81AAAAAAAAAYERwMxAAAAAAAAAY\nEdwMBAAAAAAAAEbEtNxApMsmGsmWLVta7WELWv7d3/3d6rFTTjllv59nw4YNrbZvnDIdpPBjD8xM\nYbneJwXne4hsCgv15/Zw3PQzHpqcgmZ9U4gUtuqB2ykk3IOVU3j8yy+/3Ppz2nRhmKR6ejhtCin2\nWqRNF1zafMZfPwXD+6YwaZODLjX2ce3h45KqDV+6hEynPuPPg2HktUgbPPhjHhQv1ZvCpE1iumwU\n4XNPCiT34/H6+TmwePFinX322dVjzuevFC68du3aVjvV3AOHB71JjEvj2DeJOfroo6s+HqLu53Y6\nL/yzSGPfz4u0QYDPB2nMet1TyHSXNd03R1i/fn3Vx69nUsj0MNW9y4ZaqY+P47TW+UZKaQMRn6vT\nZ+NrQBr7XTaI8k3u0qYZvh77e0gbR/i6P2wbAyW95kqp/kxTuL4HyKdwfZ8PPMDdX0eqa5yu7XzO\nSMHwvjFQ2jTGXytthuNj3ce5VF8bDBt/n2mjDa+Fz+1SvTlIOi98zk/XZL6WdNl8MK0B/j7SJjE+\n56d1zOenBx98sOqzefPmVtvn92Ha/E/K18s+TlIfn9/TBoD+efk8KNXXFGmDhS7roa8BPr+nmvt8\nn+Y4X7fSmu6bgqXvfsPGa5o2ifE+k1Xz9B2yywZRPh902RgobWTrP5euVbzuXnOprvtUjPXhuVIE\nAAAAAAAAMKm4GQgAAAAAAACMCG4GAgAAAAAAACNiWmYGpgyGLr9TvXr16sk4nE5ShsA111zTas+b\nN6/q47kH/rvk69atq37mpptu6ucQh95E/d685wykjDLPFeySOeL5GCkz0PMn0rnseQpPPvlk1cfz\nclLmwviMkddffz1+fsOUO5KOpUv+n+fTpMyYWbNmtdopD8YfS2PWnyflRXndTzjhhKqPn4OvvPJK\n1WfTpk2tdsqV8bwJz1WS6s+wS4bGVPLj6TfHr0t2iWd6pM/C+3gOnFTX1NunnXZaq71s2bIqNyzl\nm3hW3KpVq6o+ngH7wgsvVH08Ay/NV1PJP+d+snukOlPO58809n1spef155k9e3bVx8e+tyVp/vz5\nrXaaH7w2Ps4lac2aNa22nxdSzkd06b1OFa95yon0x9K47nLu+mt5Ht/OnTurdTXlE/q1QJc1fMGC\nBVUf/7m0jnn+mZ8H27Ztq37Gr0uGaf2WuuUvp3PS5+5Umy7nsl8f+zzT5bVTrqCvJWeccUbV5+ST\nT26107nsNU2ZsD7WX3rppaqPZ0UOcpwnfjwpv83HWpq7/edSzrt/PilH08dJurbzcy718Zyw9J3N\nnyetz17jLvO7zyHDlAcr5fm9y9hPP+d83ktZej6/d8n8TvOMz++e/elrvFRfM6a12cf6xo0bqz5+\nrqRsUP8MB70G+DyXxnqvDP6kS809RzPNuT6vpJr7WE95r571nL5f+FyU5neve5ofvO5p3Ex03Ydr\nNgEAAAAAAAAwabgZCAAAAAAAAIwIbgYCAAAAAAAAI4KbgQAAAAAAAMCI6Jk6W0qZJ+mvJM2W1Ej6\ni6ZpriulvFvStyXNl7RJ0q82TVOn10+CFLjomzckHvj4zDPPTMjxpA0CPvnJT7baKYw2PeY8FPXb\n3/52q502EDkYeHioVIeCp1BuDw9OQaUe7p0CT/1z97DOFGruwaQpiNo3qfDXkeoNQ9IGIn7upnDV\n8Z/FW2+9NfCg2V7SJhrbt29vtdPn5aGxKdTZz50U8OvP02XMpnnHg4vTefrss8+22mkTAR/bmzdv\nrvr4+0jjZtg2DHEeAtxlQ4wUAuybdqQad5lzfZ1IG0X4BhO+BvjYnzt3blWbFCB91113tdr3339/\n1cfHfqq5zzNpjhtkAHkKxvZNkLps7nLmmWfu8zkk6fHHH2+1Uxizh4anTWP8sbQ5iI/9l19+ueqz\nZcuWVvvee++t+nigfL+bhXjdB1nzFOzv83laV338nXXWWVWfCy+8sNX2uXzRokXVe0+bxPgmP95O\nP5euMfx9pTnt4YcfbrV9Q5Eum2d1WdMHve772pvGhH+Gc+bMqfosWbKk1U7rYdpka7wU+u7Xg2ns\n+3mQ1hEff77GS/WcnzaOeO6551rttGnGMI1rqZ6r/ZjTtZ2vUV1qnq6Ter22VNc9zd1z585ttdNG\nEf5zafNBn6vTJgIbNmxotdN3Ub8u8veZ5vtBngfptb3uaU338efXcZK0ePHinq/v10FpQxoft2kN\n8I2AfP1JNfeNI9K1nV/Pp/nB58phu25LvKZprHsfr3nahMtrnt63f17p+5hfM3apebqn45vMpXXM\n657uz3jd/T1Ig5nfu7zCm5J+r2maJZLeJ2l5KWWJpD+QdEvTNAsl3TLWBgAAAAAAADCket4MbJpm\na9M0q8f+/JKk9ZLmSvqYpBvGut0g6eOTdZAAAAAAAAAADlzZn18hKKXMl3SbpHMkPd40zbFjjxdJ\nz+9uS9KOHTv2PHH6J7M4uC3dsEyStGrR3QM+EgAAAAAAgNGxcOHCPX+eOXNm9fv6PTMDdyulzJD0\nt5L+U9M0L47/HfCmaZpSyl7vKi5fvlwrVqzQ8uXLu77cPl1xxRXVYx/5yEd6/pzflPze975X9fHf\n/U+5b+edd16rvXTp0qqPZxOkrAS/Eev5aJJ00003tdrr16+v+gyjOz+/6//7rXnKDvLf40/ZcP5Y\nynbxPKGUMZJyuMZL9fTf80/P65lWnh8l1RmBKVvCcxlS7sDu8+vP//zP9dnPfnbg2UG9eI6EJJ1/\n/vmt9s///M9XfTxvIuWReSZLyntwKRvEz8v0mfpzP/HEE1Ufz75J2UFbt25ttT1DJr1+l4yRYcsc\n8TH63ve+t+pz6aWXttqeFSfVYyBlbT799NOtdjpXfH5I54FnXPlYHz/OTz/9dD388MN67LHHWn3W\nrl1bPa/nCz3/fB3D2yW/1B9LNZ/K88DnyzQvv+9972u1L7vssqqP5zj55/Poo49WP+PnQaq51zNl\nUvp5kdYIX8M9r1CSHnjggVY7zQ+eAZtq5TVO58Ewjf2UxXb55Ze32ldddVXVxzOk0pp5993tv3Qc\nnx131VVX6R//8R+rLC+vp1Rf76Wx7/yaUarXm3TMfq74/NXlmjHlMHqftEb5fDVR50U6Zs9o8ppL\n0oc//OFW2zOcpHrNvOOOO6o+fn3s+Xsp39FrnnIF/Tzosu6nfDvPOUwZpz6vpDU9PTZIXvfx1+pf\n+tKX9KMf/aj6GR/rZ5xxRtXH18x//ud/rvr4OprWTB/rqcb+/cIzwqS67ql+Xve0JnXJDeuVEZjG\nmo/j9B1kstaA9F3Zr9d9nEvtmxSStG3btqqP5+p63qpU1yLN754pl7IjfY7wOTZ9H/MMyFRzn4vS\n9YOP6wMd5xN532VvvO7pO5rX3Wvu1+WStGbNmla735p7dmTKA/Sap9p4/dI/cvPvcf4z6bmnan5f\nuXLlPv97p1mhlHKYdt0I/FbTNN8Ze/jpUsqcsf8+R1I9ggEAAAAAAAAMjZ43A8d+BfgvJa1vmuZP\nx/2nmyRdO/bnayXV/8wOAAAAAAAAwNDo8mvCl0j6DUlrSyn3jD32RUlfkXRjKeU3JT0m6Vcn5xAB\nAAAAAAAATISeNwObpvl/kupwgl3q8D4AAAAAAAAAQ6nzBiLDJAXpd+GhlV/4wheqPh7mPXfu3L5e\ny6XQcg8S//a3v131ScGyoyAF33pAs4d0S3WIegrZ9Y0qjjzyyKqPB4p6cHAKGPVNPdKGMH5+peBU\nD0V97bXXqj7++inQevxjw755iFSH5kt1GHMKdPdzJQXEdgkA99d/4YUXevZJ56BvWJA2COiyoYnX\nOAU/H3bYYa12l1DpYeOfaQpe9xDgtFmCbx6U5m4PkE7h/76+pPPS6/fUU0+12uODz08//XTdeuut\n1ftKY9/nkFQ7D0n2YPH0c8N2DqQ5zQOivVZS/V59c4JFixZVPzNr1qxWu8u4TueFbwaSNqjxPr4J\nkFQH3Ke1xN9n2vjA+wx6k5he0mfqwetpzPrY9w2jpDoIfsmSJXv+/Prrr+sTn/hENWZ9Q5F0jGkz\nAh+3aRx7CH5672+88Ub1WC9pwxDXZa2fyvPC3/v4zV12881l0uYNfv2ewv/POeecVttr02W+T3OT\nnytpEwEf6+n6L72+6zKufZ0ftus7X8d8Yyyp3vgjbSrlNff5Xqo3meuyrqbvY34epHXCx3XaGMhf\nPz2Pj/1UY99EwGveZQxP5ThPn6lf89x///1VH/8+dsopp1R99jW/7+ZjMt0n8M/dz4v0PD6u06Zg\nfl502Rioy3XbdOB1T9fvXnev+amnnlr9jH+PG2TNpbruaaMbr3uXa7thqflwHAUAAAAAAACAScfN\nQAAAAAAAAGBEcDMQAAAAAAAAGBHTMjPw9ttvrx7zrJlly5b19dwTlRG4Zs2aVvu2226r+ozPlUJb\nysbxXKf0O/ueGZCy2DwfIGVTeUaT/+5/yn7x5025Mp4pkPLIPOMgfRbDlhEzEVLmiGc3pEwIz+aZ\nP39+1WfGjBmtdspy8Fwgz4GT6nMuZQZ6jdMx++unrD/PlkgZI25Y8if2h9c9ZbH5Z5jGlmeKpEwp\n/5zT+PPMn5Qf4o9527PG7rnnniqrJI3rLnmA0/E88PkqjYnNmzf37OPjbfHixa32scceW/2Mf84p\nz8dzPf1YpHo+SGPf15s0p/k5mDJrPRNz2PKh+tElJzLNyz7W07Wd50KOn+9ff/11zZgxo9Mc4mM/\nZYJ53dP55GM95QP2WsPTZ+FStvIgpffkdX/00UerPj5GU87a0qVLW+2UDeyZcp6pm/J7vVZdciJT\nHuBzzz3Xaqfz3ce+58J1NWzXf348/t7Td51bb7211U7X6hdffHGr7dmSkjRnzpxW29dQqa57Og98\nfk9j36//0hzi7yPVys/Lg0Gar/xzTt/dfY284IILqj4+rj1HNknXkT5G07Wd191rnuamdO66frLi\n0vw+bOu81z2NLa/7sNXcjzndW/B1IX13cP3mQnrdp6Lmw3VWAQAAAAAAAJg03AwEAAAAAAAARgQ3\nAwEAAAAAAIARwc1AAAAAAAAAYERMyw1EkhtvvLFnn343FXEeRP3DH/6w6vNP//RPE/JaoyqF7npQ\nqQc/S3VwcQr49XD2FPKZNnQYL4X/e0h4Cg3399AlCHrYwqInSwrL9Xru3Lmz6uMBsevWrav6+IYw\nidcmhf/78aQa+7mRziUPue51vknDFxw8UbpsLuGhv15zSbrvvvta7bQxgwd3p3POXz8Fwfc6L/1c\n2rlzZ3UOppqPynmQPncPZO4S4L569epW++ijj65+xj+vNHd7zVM4dJex7/U76qijeh5PFwdDzdM6\n5ps2rV+/vurj50HajM03Cxo/9q+++mpdf/311c+kce1B8H58qU+ar/wcS+d7r80/0uc1Ha8F/JjT\nhisPPvhgq51C3u+6665WO20mccwxx7TaPm5Srfx40tri4fFpAzlfA1J9fX44WK//utTcx1/asO3O\nO+9stU866aSqj28skOZKn89TjX1TmNTH656uRydrrh62zYJcv5sH+ed+zz33VH18M4m0zvvrp++H\nvvlHqrGfq/48aaOUg2F97levzYOkuu7TreZSvm500/k8mL5HDgAAAAAAAGC/cDMQAAAAAAAAGBHc\nDAQAAAAAAABGRM/MwFLKPEl/JWm2pEbSXzRNc10p5Y8kfVrSM2Ndv9g0zfcn60D3V8oQ7JIriOHl\n+QApo8JzNVK+g+fGdMnp6mI6ZrsMO/9MUz39sZRZMVk1PuSQQ6o+6TF01yUrq0uOX5eaT9SY9awQ\nzyX1Nmr9jHXPgkoZsZOFsb9/uozrlMHlNU2ZPymLarerr7465gx2PcZ++uBnuly3+TVZyvbznKmU\nL9lPZlOXLLYu7wE/459Pl2zefmvez7Vdyv8a9rE/7HlkEzW/ez6nJG3YsKHVTjXv5/vhsOeyDnvN\npW6fu9edmu/bIOreZQORNyX9XtM0q0spx0j6SSnl5rH/9j+bpvnvk3d4AAAAAAAAACZKz5uBTdNs\nlbR17M8vlVLWS5o72QcGAAAAAAAAYGKV/fnnkaWU+ZJuk3SOpP8s6T9I2iFplXb968E9/9Zzx44d\ne55448aNE3O0mDaWblgmSVq16O4BHwkAAAAAAMDoWLhw4Z4/z5w5s/qd6843A0spMyT9k6Q/aZrm\nO6WU2ZK2a1eO4JclzWma5j/u7j/+ZuCVV16pFStWaPny5f2+D0wzd37+LknSxdddNOAj6e1gzwxk\n7B38NUZtKjMD94XxB0ycNK73Nda/+tWv6nOf+1yn5x723LBR1iVHiczA4XMg65/XM9V3VDIDR0Wq\nsWfxjlJ+3IGYLtee1HzyrVy5cs+f083ALpmBKqUcJulvJX2raZrvSFLTNE+P++9fl/T3B3qwwCCM\nwkQw6qjx6KHmwMGnn4v5Ljd8MNy61JA6H1y8ntT34JdqTN0PbtR88Hr+NVrZdTv2LyWtb5rmT8c9\nPmdct09Ium/iDw8AAAAAAADAROnyLwMvkfQbktaWUu4Ze+yLkq4ppZyvXb8mvEnSZyblCAEAAAAA\nAABMiC67Cf8/SSmU4fsTfzgAAAAAAAAAJsv+p+0CAAAAAAAAmJa4GQgAAAAAAACMCG4GAgAAAAAA\nACOCm4EAAAAAAADAiChN00zKE+/YsWNynhgAAAAAAABATzNnzqw2BeZfBgIAAAAAAAAjgpuBAAAA\nAAAAwIiYtF8TBgAAAAAAADBc+JeBAAAAAAAAwIjgZiAAAAAAAAAwIib9ZmAp5SOllAdLKQ+VUv5g\nsl8PGHWllE2llLWllHtKKavGHnt3KeXmUsrGsf9/16CPEzgYlFK+WUrZVkq5b9xjex1vpZT/OrYe\nPlhK+fBgjhqY/vYy9v6olPLE2Pp3Tynlo+P+G2MPmCCllHmllB+VUtaVUu4vpXx+7HHWP2AS7WPs\nsf5hv01qZmAp5RBJGyT9W0lbJN0t6ZqmadZN2osCI66UsknS0qZpto977L9Jeq5pmq+M3ZR/V9M0\n/2VQxwgcLEopl0l6WdJfNU1zzthjcbyVUpZI+mtJF0k6SdIPJS1qmuatAR0+MG3tZez9kaSXm6b5\n79aXsQdMoFLKHElzmqZZXUo5RtJPJH1c0r8X6x8wafYx9n5VrH/YT5P9LwMvkvRQ0zSPNE2zU9Lf\nSPrYJL8mgNrHJN0w9ucbtGvRAHCAmqa5TdJz9vDextvHJP1N0zQ/bZrmUUkPadc6CWA/7WXs7Q1j\nD5hATdNsbZpm9difX5K0XtJcsf4Bk2ofY29vGHvYq8m+GThX0uZx7S3a98kK4MA1kn5YSvlJKeW3\nxx6b3TTN1rE/PyVp9mAODRgJextvrInA5PtcKWXN2K8R7/4VRcYeMElKKfMlvVfSnWL9A6aMjT2J\n9Q/7iQ1EgIPPB5qmOV/SVZKWj/0q1R7NrmyAycsHALAH4w2YUtdLOk3S+ZK2Svofgz0c4OBWSpkh\n6W8l/aemaV4c/99Y/4DJE8Ye6x/222TfDHxC0rxx7ZPHHgMwSZqmeWLs/7dJ+q52/VPwp8cyJnZn\nTWwb3BECB729jTfWRGASNU3zdNM0bzVN87akr+tnvwrF2AMmWCnlMO26GfGtpmm+M/Yw6x8wydLY\nY/1DPyb7ZuDdkhaWUhaUUg6X9ClJN03yawIjq5Ry9FiYrEopR0u6UtJ92jXurh3rdq2k7w3mCIGR\nsLfxdpOkT5VSjiilLJC0UNJdAzg+4KC0+ybEmE9o1/onMfaACVVKKZL+UtL6pmn+dNx/Yv0DJtHe\nxh7rH/px6GQ+edM0b5ZSPivp/0o6RNI3m6a5fzJfExhxsyV9d9c6oUMl/Z+maX5QSrlb0o2llN+U\n9Jh27TgF4ACVUv5a0i9IOr6UskXSlyR9RWG8NU1zfynlRknrJL0paTm7uQH92cvY+4VSyvna9auJ\nmyR9RmLsAZPgEkm/IWltKeWesce+KNY/YLLtbexdw/qH/VV2xTkAAAAAAAAAONixgQgAAAAAAAAw\nIrgZCAAAAAAAAIwIbgYCAAAAAAAAI4KbgQAAAAAAAMCI4GYgAAAAAAAAMCK4GQgAAAAAAACMCG4G\nAgAAAAAAACPi/wOuPf5escm8aAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f33a8ff67f0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(image_idx=(0, 1000))\n",
    "def show_num_components_reconst(image_idx):\n",
    "    fig, ax = plt.subplots(figsize=(20., 20.))\n",
    "    actual = X[image_idx]\n",
    "    x = np.concatenate([actual[np.newaxis, :], reconstructions[:, image_idx]])\n",
    "    ax.imshow(np.hstack(x.reshape(-1, 28, 28)[np.arange(10)]),\n",
    "              cmap='gray');\n",
    "    ax.axvline(28, color='orange', linewidth=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also browse throught the reconstructions for other digits. Once again, `interact` becomes handy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAARkAAACaCAYAAACdQzzVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFI5JREFUeJztnXuwVNWVh78VkIBKQBAUAUUNmiCgk0gUEyMqEFDLKJZR\nw0xZJuWQmKlKfIQyTMVHJZMYa2ocK5VUGTWJU4zBUVCQImMA4yuhFDVGCYhP3i+faGK4V3DPH32a\nOWfvze2+fXtzu5vfV0Xds3avPr3uvX0X5/x67bXMOYcQQqTiY90dgBCitVGSEUIkRUlGCJEUJRkh\nRFKUZIQQSVGSEUIkpWdXnmxmU4BbgR7AHc65m8qPbd++XZ+NC7GP0a9fP/PXar6SMbMewM+AqcAo\n4BIzG1V7eEKIVqQrt0ufA15xzr3mnGsH5gBfrk9YQohWoSu3S0OB9Tl7A3BSzPHll1/uwssIIRqZ\nkSNHdvh4lzSZapk2bRrz5s1j2rRpe+Pl6kKzxQuKeW/RbDGnjveFF17o8PGu3C5tBIbn7GHZmhBC\n7KYrSWY5MNLMjjSzXsDFwIL6hCWEaBVqvl1yzu00s38BHqL0EfYvnXN/qVtkQoiWoEuajHNuEbCo\nTrEIIVoQVfwKIZKiJCOESIqSjBAiKUoyQoik7JViPCGajY8++qhg+72wd+3aRc+exT8fs2BvYPA8\n/7wAO3fuLNgf//jHazqPTyye7kBXMkKIpCjJCCGSoiQjhEiKkowQIikSfhuIsWPHBms33XRTwf7Y\nx8L/F6655prdx6NHj2bFihX1D64BqXUwoS+I7tq1K/DxhdXevXsX7F69egUC7d///vfgPOvWrSvY\nQ4cODXzOP//8Dm0IxeH58+cHPkuXLi3Yf/vb34Ln+++fmDhcb8FYVzJCiKQoyQghkqIkI4RIijSZ\nbmTixIkF+7zzzgt8qinC+sY3vlE4XrJkSeCzYEGx1U81xVzNSDV6gq/B+EV1EGog7e3tge3/DP3n\nABx//PEF+5RTTgl8pkyZUrAHDx4c+Bx66KEFu0+fPoGPH+PDDz+8+7hXr14AfPjhhwWfmMZXb7o6\nEmUN8D6wC9jpnDuxHkEJIVqHelzJnO6ce7MO5xFCtCDSZIQQSbFaaw0AzOx1YDul26XbnHO/KD+W\nnyCpkShCtC75kSixCZJdvV36gnNuo5kNBhab2YvOucd8J41ECUXe2NqYMWMqnqejnb5tbW3RHbwA\nl112WcHetm1bxdfaG+zp5xz7Pn2hNSZaViP8+kKvL4ZCWFiXF1WXLFnCxIkTg9c67rjjgvN88Ytf\nLNjHHHNM4HPkkUcW7L59+wY+Bx98cMWYly1bVrBvu+02AG655RauvPJKICwOjNHZC4+UI1Fwzm3M\nvm4D7qc0VVIIIXbTlVnYB5hZ3/IxMBnYN+rZhRBV05XbpUOA+7NLxp7A3c65/61LVEKIlqErc5de\nA46v6NiCHHDAAQX76KOPDnyuvvrqgt2vX7/AZ7/99qv4WuvXry/YMR0itumu2alGF4j5+Gsxjaqt\nra1DG8Lf8eGHH16wP/nJT/KZz3ymsHbqqacG5zniiCMK9qBBgwKf/EZGgLVr1wY+H3zwQfD6Pn7B\nXt4uH69atargEyvq0wZJIURToSQjhEiKkowQIilKMkKIpGgXdgXGjx8frE2dOrVgf/aznw18fPGs\n1srq++67r6JPudCqlahG1I3tJPeF8Zio++677xbsww47LPD51Kc+VbBPO+20gj19+nROPvnkwtrw\n4cOD85R3P5fxBVyAl156qWA/++yzgc/+++9fsGMfJPhFhnlRt3zsd/jbG+hKRgiRFCUZIURSlGSE\nEEmRJuNxxhlnFI7zkwA6g68NpOxEl9d/GmU0aXfhbxz0u8VB2Hlu3Lhxgc+kSZMKdr7wbseOHUya\nNIkePXoUfLZs2RKcx9dOBg4cGPj4esvGjRsDH19LeeeddwIfX3/asWNHcOx3BfQ1I6hdP9wTupIR\nQiRFSUYIkRQlGSFEUpRkhBBJqSj8mtkvgXOAbc650dnaAOAeYASwBviKcy5UohqcvMhbxh8vEhPB\nfDHRL+6CcHdrrNtZpfNC2KEttms2H2O9RbvuIvZ9+LvW/d3LEAq/w4YNC3z8AstLLrkk8PEL6/LC\n65YtW+jfvz+PPPJIwef9998PznPOOecU7AEDBgQ+/vfxxhtvBD6+gBzbwe9/2BB7X7z99tsFnxEj\nRgTn8QsGu/phQjVXMr8Gpnhr1wJLnXMjgaWZLYQQARWTTNaz921v+cvAXdnxXUA4lUwIIahyWoGZ\njQAW5m6X3nXO9c+ODXinbJfRtAIh9g1STyvAOefMrMNM1SjTCvx78e9///uBTzWd//3u7NdeG94t\n+pMIqtnEWO4un8cfLxubelA+dytNK4jhF47Fvof+/Qv/13HuuecGPhdffHHBzv+RlPGL3ebNm7f7\nePLkyfzud79j+fLlBZ9YUZ8/rSCm2zz00EMF29d6AD796U8X7B/96EeBj19oN2vWLABmzpzJzTff\nDMBvf/vbgk9Ms/rrX/8arPnkdZpU0wq2mtmQ7MWGAI3xjhVCNBy1JpkFwKXZ8aXA/PqEI4RoNSom\nGTP7DbAMONbMNpjZ14GbgElm9jIwMbOFECKgoibjnAuLCEqcWedYhBAtSMvuwo4JpN/85jcrPi9f\nENfe3s7q1asDn5///Oedjue1114L1hYvXlywFy5cWPE8jz/+eLB21llnAaXCqtWrV3Psscd2Or5G\nI1Zslt9VDKHIC3D66acX7OnTpwc+Q4YMKdibN28OfObOnVuwZ8+evft48uTJ3H777cHYFH/cLIS/\nY3/EDcCrr75asN96663Axy/Qi4mz/m7u/M7/PXUBqGbUb1cLPLWtQAiRFCUZIURSlGSEEElpWU0m\ndi9eTaf2OXPmAHDRRRcxZ84c7rnnnppef+XKlQXbL9wC2L59e6fPG+u+n98U6G8QbCV27txZsP2J\nAgBTphS32cXGufqjWh977LHA59577y3Yvpayfv36YMpBbMqAr6/E9A1/E2zsfeFvovQ3OkJYrJjX\nW8rH/s8wpsn47I0NkkIIUTNKMkKIpCjJCCGSoiQjhEhKywi/Rx11VMH2C5MgFLDKRWwxLrrooppF\nX4BNmzbV/NzO4o9EaYWxKLEugX63ulGjRgU+Y8aMKdixXc9//OMfC/YDDzwQ+PiCrV/4FysEjI0y\n8QVbv4APQvE1NqbE9/HFYoD33nuvYB944IHBsT/GxS9wTIGuZIQQSVGSEUIkRUlGCJGUWqcV3ABc\nDpR3bc1yzi1KFaTPEUccEaz5Xe7y96NlWqGTf6zzXc+ePQvHzfh9+jH7RWMQdrA78cQTAx9/msMz\nzzwT+PjFd6+88krg84lPfKJg+1pG3759g5h9TQRCfSWmNXVURFfG/73HtCZfJ8rHXD72tcqYJlNv\nTa/WaQUAtzjnTsj+7bUEI4RoLmqdViCEEFVR67SCG4DLgO3A08DV/nA3TSsQYt+g0rSCWpPMIcCb\ngAN+AAxxzn0t/5x8khkzZkxdpxXENJkbbrihYA8ePLjiec4+++w9PtYI0xVixDSZH//4x0CpadWa\nNWuiTasafVpBNZqM35Aq1pjstNNOK9ixTYu/+tWvCrZfNwMdazL33nsvF154YRBzrHalGr3F94lt\nkBw9enTBjjVg8+uI7rzzTgBmzJixexLG/PnFdtyx91NnNZn8tIK6jURxzm3NBXQ7ULmlWx254oor\ngrVqkkorcOqppwZr5aTS1tbWEl3xIP4He/jhhxfs2C5sf5f6888/H/g8/fTTBTsmfg4cOLBg+8mh\nZ8+eQSKMFdH5f7CxJOOPOfYTHITjZGPFgB2N8S0f+wJ2LJ56f3BQ00fY5XEoGecDK+oTjhCi1ajm\nI+zfABOAg81sA3A9MMHMTqB0u7QGmJEwRiFEE1PrtII7E8QihGhBWmaDZDXccccd3R1Cpxk6dGjB\nvvzyyys+Z+vWrcFarKNeI+FrF7ENrv5I1QEDBgQ+vr6xZcuWwMcfQVuN+OnrLzt37gy0i5hg6j8v\npjX5rx8rMpw0aVLBjmmQvtaU/1S3fOyLyoMGDQrOU8331Rm0rUAIkRQlGSFEUpRkhBBJUZIRQiRl\nnxJ+YztXGwlf5AW48cYbC3bfvn0Dn7LY2bt3b7Zv384Pf/jDwKeW8SvdSUyM9Yvdqil+O+iggwIf\nf5RJbIyMfx5fOG9rawu63sVGwfrn9sVrgDPPLI6Vj1Wi+1Xu/mhbgKVLlxbs119/PTj2d6nvDXQl\nI4RIipKMECIpSjJCiKQ0pSYTKw6qpmDoqquuKthLliypW0yViGkMM2fOLNjjx4+veJ5YcVm5K+BP\nf/pTvvvd70a75jc6fgFYrIPc5s2bC3ZsVKtfxDdhwoTAxx85G+uM55/bj6e9vT34nfbr1y84j7+J\n80tf+lLg4xffxbo6vvjiiwX7wQcfDHz83eT+FAsIpyXENkN2R2c8IYSoGSUZIURSlGSEEElRkhFC\nJKWafjLDgf8CDqHUP+YXzrlbzWwAcA8wglJPma/4fX5Tcffddwdrs2bNKtixcaA+N998c7CWF8J+\n8pOfsGzZssBnw4YNBfvCCy8MfHzxLD+2pIzfxS62U3rOnDkF+4knngh88kJvM4q+MWKFbfk2jwBH\nH3104OOPqT3uuOMCn+nTpxfs2NiUtWvXFuxdu3YV7LPPPjvYCR3bGe2LurGYfTH/kUceCXz8Dyme\nfPLJwMcXcWNjamNtTX38925XxeFqrmR2UmoUPgo4GfiWmY0CrgWWOudGAkszWwghClQzEmWzc+7Z\n7Ph9YBUwFPgycFfmdhdwXqoghRDNS1XTCnY7l6YWPAaMBtY55/pn6wa8U7ZBI1GE2FeoNBKl6mI8\nMzsQmAt8xzn3Xv6ezDnnzGyP2WratGnJR4z4IyOuu+66wMfXaWL3leWk29bWFi2gq5Zq7mt9jSFW\nHNiZgsFGHePSEXuKOaYd+L+/cePGBT4XXHBBwY5Nb/A7/b/55puBj1+wl//9HXbYYWzatCkovot1\nmfPfQ8uXLw98Hn300YL91FNPBT7r1q0r2DGNzy9ELGt8DzzwAOedV7rRqEZL6awm47+Pfar6dMnM\n9qOUYP7bOTcvW95anlqQfW2MQT5CiIaiYpLJboXuBFY55/4j99AC4NLs+FJgvv9cIYSo5nbp88A/\nAS+Y2XPZ2izgJuB/zOzrwFrgK2lCFEI0M9WMRHkC2NON3Jl7WBdCCKBJd2HHWLGiOMQyNit46tSp\nBfurX/1qsnj8sRt+fAC33nprwf7ggw+SxdNsxIRNXwyOiai+SOnPxgY46aSTCnbv3r0DH3+udH6c\na3t7OyNGjAgKBtesWROcx4/x8ccfD3xWrVpVsGM70P2OdrEPJfzn5X8W5ePYWNrUaFuBECIpSjJC\niKQoyQghktIymozPW2+9FazNnj27YMe6zJU3Ow4aNIgNGzZEu8v7hVr33Xdf4LNp06aCvXLlyspB\niw7x9YSYhuUXL+Y79pfxNZCxY8cGPv5UiPz7ady4cfzhD39g9erVBZ/YZlp/qkA1RYZ+sSCEG0Z3\n7NgR+Ph6VP7n1RUtRmNqhRANjZKMECIpSjJCiKQoyQghktKywm81dLTred68ecyYMWNvhyQyYjt/\nfdG0R48egY9fWBcb5+oLv7HRrf658+NmFy1axPXXX1/xORDujI75+KJurBgv1inQpyNxtyze+j/X\neo8/iaErGSFEUpRkhBBJUZIRQiSlK9MKbgAuB97IXGc55xalClTsW3TUtbAzxDY/+mvV6BL+hs0+\nffoEEwyqOY//nGqpl3ayNzQYn2qE3/K0gmfNrC/wjJktzh67xTn37+nCE0I0O9X0k9kMbM6O3zez\n8rQCIYSoSFemFVwFXAZsB56mdLWzu4mKphUIsW9QaVpB1Ukmm1bwKPBvzrl5ZnYI8CYlneYHwBDn\n3NfK/vkkM2bMmKbrpN9s8YJihtp0m8528J87dy4XXHBBVZpMd2ggPqnfF/lpBTWPRIlNK3DObc09\nfjuwsKvBCtFVUv1R+8nLOdctXeaakZqnFZTHoWScD4T9JYUQ+zxdmVZwiZmdQOl2aQ2gGnwhREBX\nphWoJkYIURHdVAohkqIkI4RIipKMECIpSjJCiKR0quK3M+SL8YQQ+waxYjxdyQghkqIkI4RISrLb\nJSGEAF3JCCESoyQjhEhK8iRjZlPMbLWZvWJm16Z+vVows1+a2TYzW5FbG2Bmi83s5ezrQd0Zo4+Z\nDTez35vZSjP7i5l9O1tvyLjNrLeZPWVmf87ivTFbb8h485hZDzP7k5ktzOyGjtnM1pjZC2b2nJk9\nna11W8xJk4yZ9QB+BkwFRlHaVDkq5WvWyK+BKd7atcBS59xIYGlmNxLltqijgJOBb2U/20aNuw04\nwzl3PHACMMXMTqZx483zbSA/rKkZYj7dOXeCc+7EzO6+mJ1zyf4B44GHcvb3gO+lfM0uxDoCWJGz\nV1NqxAUwBFjd3TFWiH8+MKkZ4gb2B54FTmr0eIFhlP4ozwAWNsN7g1JXhIO9tW6LOfXt0lBgfc7e\nQPP0Bz7ElfobA2yhNK2hIcnaov4D8CQNHHd22/EcsA1Y7Jxr6Hgz/hOYCeRHODZ6zA5YYmbPmNk/\nZ2vdFvM+Paa2Wpxzzswa8rP+rC3qXOA7zrn38p3hGi1u59wu4AQz6w/cb2ajvccbKl4zOwfY5px7\nxswmxHwaLeaMLzjnNprZYGCxmb2Yf3Bvx5z6SmYjMDxnD8vWmoGt5e5/2ddt3RxPQKwtKk0Qt3Pu\nXeD3lHSwRo7388C5ZrYGmAOcYWazaeyYcc5tzL5uA+4HPkc3xpw6ySwHRprZkWbWC7gYWJD4NevF\nAuDS7PhSSppHw7Cntqg0aNxmNii7gsHM+lDSj16kQeMFcM59zzk3zDk3gtJ792Hn3D/SwDGb2QHZ\nfDTM7ABgMqXWuN0X814Qoc4CXgJeBf61u0WxPcT4G0qzpT6kpBt9HRhISfB7GVgCDOjuOL2Yv0Dp\n3vt54Lns31mNGjcwFvhTFu8K4LpsvSHjjcQ/gf8Xfhs2ZuAo4M/Zv7+U/+a6M2ZtKxBCJEUVv0KI\npCjJCCGSoiQjhEiKkowQIilKMkKIpCjJCCGSoiQjhEjK/wEqjhQgatrWIwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f33a9312630>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(i=(0, 10))\n",
    "def show_pca_digits(i=1):\n",
    "    plt.figure(figsize=(4,4))\n",
    "    actual_sample = X[i].reshape(28,28)\n",
    "    reconst_sample = (reconst[i, :] * std + mu).reshape(28, 28)\n",
    "    plt.imshow(np.hstack([actual_sample, reconst_sample]), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. PCA for high-dimensional datasets\n",
    "\n",
    "Sometimes, the dimensionality of our dataset may be larger than the number of data points we\n",
    "have. Then it might be inefficient to perform PCA with the implementation above. Instead,\n",
    "as mentioned in the lectures, we can implement PCA in a more efficient manner, which we\n",
    "call PCA for high-dimensional data (PCA_high_dim)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the __normalized__ data matrix $\\boldsymbol{\\bar{X}}$ of size $N \\times D$ where $D > N$. To do PCA we perform the following steps:\n",
    "\n",
    "* We solve the following eigenvalue/eigenvector equation for the matrix $\\frac{1}{N} \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T$, i.e. we solve for $\\lambda_i$, $\\boldsymbol c_i$ in\n",
    "$$\\frac{1}{N} \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i = \\lambda_i \\boldsymbol c_i.$$\n",
    "\n",
    "* We want to recover original eigenvectors $\\boldsymbol b_i$ of the data covariance matrix $\\boldsymbol S = \\frac{1}{N} \\boldsymbol{\\bar{X}^T} \\boldsymbol{\\bar{X}}$.\n",
    "\n",
    "* Left-multiply the eigenvectors $\\boldsymbol c_i$ by $\\boldsymbol{\\bar{X}}^T$ yields\n",
    "$$\\frac{1}{N} \\boldsymbol{\\bar{X}}^T \\boldsymbol{\\bar{X}} \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i = \\lambda_i \\boldsymbol{\\bar{X}}^T \\boldsymbol c_i$$ and we recover $\\boldsymbol b_i=\\boldsymbol{\\bar{X}}^T \\boldsymbol c_i$ as eigenvector of $\\boldsymbol S$ with the eigenvalue $\\lambda_i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "def PCA_high_dim(X, num_components):\n",
    "    \"\"\"Compute PCA for small sample size. \n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the data,\n",
    "           and N is the number of data points in the training set. You may assume the input \n",
    "           has been normalized.\n",
    "        num_components: the number of principal components to use.\n",
    "    Returns:\n",
    "        X_reconstruct: (N, D) ndarray. the reconstruction\n",
    "        of X from the first `num_components` principal components.\n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    M = 1 / N * X @ X.T # EDIT THIS, compute the matrix \\frac{1}{N}XX^T.\n",
    "    eig_vals, eig_vecs = eig(M) # EDIT THIS, compute the eigenvalues. \n",
    "    U = eig_vecs[:, range(num_components)] # EDIT THIS. Compute the eigenvectors for the original PCA problem.\n",
    "    # Similar to what you would do in PCA, compute the projection matrix,\n",
    "    # then perform the projection.\n",
    "    P = projection_matrix(U) # projection matrix\n",
    "    print(P.shape)\n",
    "    X_reconstruct = P @ X  # EDIT THIS.\n",
    "    return X_reconstruct"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the same dataset, `PCA_high_dim` and `PCA` should give the same output. \n",
    "Assuming we have implemented `PCA` correctly, we can then use `PCA` to test the correctness\n",
    "of `PCA_high_dim`.\n",
    "\n",
    "We can use this __invariant__\n",
    "to test our implementation of `PCA_high_dim`, assuming that we have correctly implemented `PCA`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1000, 1000)\n",
      "correct\n"
     ]
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(PCA(Xbar, 2), PCA_high_dim(Xbar, 2))\n",
    "# In fact, you can generate random input dataset to verify your implementation.\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's compare the running time between `PCA` and `PCA_high_dim`.\n",
    "\n",
    "__Tips__ for running benchmarks or computationally expensive code:\n",
    "\n",
    "When you have some computation that takes up a non-negligible amount of time. Try separating\n",
    "the code that produces output from the code that analyzes the result (e.g. plot the results, comput statistics of the results). In this way, you don't have to recompute when you want to produce more analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def time(f, repeat=100):\n",
    "    times = []\n",
    "    for _ in range(repeat):\n",
    "        start = timeit.default_timer()\n",
    "        f()\n",
    "        stop = timeit.default_timer()\n",
    "        times.append(stop-start)\n",
    "    return np.mean(times), np.std(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "times_mm0 = []\n",
    "times_mm1 = []\n",
    "\n",
    "for datasetsize in np.arange(4, 784, step=20):\n",
    "    XX = Xbar[:datasetsize]\n",
    "    mu, sigma = time(lambda : XX.T @ XX)\n",
    "    times_mm0.append((datasetsize, mu, sigma))\n",
    "    \n",
    "    mu, sigma = time(lambda : XX @ XX.T)\n",
    "    times_mm1.append((datasetsize, mu, sigma))\n",
    "    \n",
    "times_mm0 = np.asarray(times_mm0)\n",
    "times_mm1 = np.asarray(times_mm1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAEGCAYAAADylEXaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VNX5wPHvmTV7QtgEZJVQAUEWQRAE3BCCSpWiICiC\n1apVf2JttbZaK1gt1aKlVluVrdaNRUVFFhd2o6yyiBJkRwiBkD2Z9fz+mMlkErJMklmyvJ/nycOd\ne8+9950B5s09q9JaI4QQQjQ1hkgHIIQQQkSCJEAhhBBNkiRAIYQQTZIkQCGEEE2SJEAhhBBNkinS\nAYRSTk6OdHEVQohGLjExUdXmPHkCFEII0SRJAhRCCNEkSQKsRnp6eqRDqDGJOTwk5vBoaDE3tHih\nYcYcDJIAhRBCNEmSAIUQQjRJYesFqpQaBbwEGIHXtdbPlTuuvMdTgULgDq31NqVUe2Ah0BrQwH+0\n1i95z3kKuAvI9F7mca318upi0VqTn5+P2+2uNu6oqChycnICe5P1hMQcHAaDgbi4ODz/NIUQjU1Y\nEqBSygi8DFwDHAM2K6WWaa2/8ys2Gkjx/lwKvOL90wn8xpsM44GtSqnVfufO1lo/X5N48vPzsVqt\nWCyWastarVaioqJqcvmIk5iDw263k5+fT3x8fKRDEUKEQLiqQAcC+7XWB7TWduAdYGy5MmOBhdoj\nDUhSSrXRWp/QWm8D0FrnAXuBdnUJxu12B5T8RNNmsVgCqiUQQjRM4aoCbQcc9Xt9DM/TXXVl2gEn\nSnYopToBfYGv/co9oJS6HdiC50nxbEUB+PdyioqKwmq1Bhx8cXFxwGXrC4k5OHJzczl16lSlxxti\n7zmJOfQaWrzQ8GJOSUmp8zUazEwwSqk4YAnwkNY617v7FWAGnrbBGcALwLSKzvf/sHJycgKubisu\nLq53VXPVkZiDJyEhgfbt21d4LD09PSj/CcNJYg69hhYvNMyYgyFcVaDHAf9vkfO9+wIqo5Qy40l+\n/9NaLy0poLXO0Fq7tNZu4DU8Va1CCCEiJGnecZLmlf96r3mZcAhXAtwMpCilOiulLMAEYFm5MsuA\n25XHICBHa33C2zv0DWCv1vrv/icopdr4vbwR2B26txBeb775JtOnT6d3795Mnz6defPmRTokIYRo\nVMKSALXWTuB+YCWeTizvaa33KKXuUUrd4y22HDgA7MfzNHefd/8Q4DbgSqXUDu9PqvfYLKXULqXU\nTuAKYHqo3kMofmOZO3cuDz/8sO/1zJkzufvuuwGYPHkyTz75JMOGDWP27NlMnTq1yvL+ioqKSE1N\nxeVyAZCcnMzQoUMZPHgwU6ZMobCwEICMjAymTZtGnz59GD58OOPHj2f//v2+63z88cckJSWxb98+\n3z673c7o0aNxOp1B/SyEECLcwtYG6B2ft7zcvlf9tjXw6wrO2wBUOBBLa31bkMMMqwkTJjBgwACy\ns7NJS0tj5cqVrFq1ynd8x44d9OnTJ+DyJd5++22uv/56jEYjANHR0WzYsAGAu+66i7lz5/LrX/+a\nyZMnM3HiRObOnQvArl27OHXqFF27dgVgyZIlDB48mMWLF/P4448Dnp6Rw4cPZ+nSpdx8882h+WCE\nECIMZCaYCIqJiWHcuHHMnDmTRx99lIULFxIdHe07Xj4BVle+xNKlS0lNTT1nP8DgwYM5ePAg69at\nw2QyMW1aaZ+hXr16cdlllwGesZJpaWnMmTOHpUuXlrnGmDFjWLRoUZ3euxBCRFqD6QXaWE2ePJmB\nAwfy1ltv0blz5zLHdu/ezT333BNwefBUUR4+fJiOHTuec8zpdLJ69Wquvvpq9u7dWya5lrd8+XKu\nuuoqunbtSnJycplk3KNHD7Zt21abtyuEEPWGJECoUdteIGWzpwY+Tn/WrFm0aNGiwja1N954o0bl\nAc6cOUNiYmKZfUVFRQwdOhTwPAHedtttvmrPyixZssSXfG+66SYWL17sS4BGoxGLxUJeXp7MkiJE\nE1Ly/VeT77j6TBJgBM2ZM4fi4mLmzZvHs88+yw033FDn8tHR0ecMKPdvAyzRvXt3Pvzwwwrvc/bs\nWdatW8eePXtQSuF2u1FKMWPGDN+8mDabrV6O2xNCiEBJAqTq32ZKBmgH+zeftWvX8tZbb7Fq1Sri\n4+N5/PHH2blzJ717965T+aSkJNxud7UDy4cNG8bTTz/N/PnzueOOOwBPlWtubi779u3jlltu4cUX\nX/SVT01NZdOmTQwZMoSsrCyaN2+O2Wyu+wchhBARIp1gIuDo0aM8+OCDzJ8/31eFeM899/DKK68E\npfzw4cNJS0urMgalFG+++SZr1qyhT58+DBo0iKeffprWrVuzePFirrvuujLlb7jhBpYsWQLA+vXr\nGTlyZI3esxBClHCumYRzzaRIhyFPgJHQvn17vv322zL7Jk2axKRJFf+DqGn5qVOn8vrrrzNixAgA\njh+vuN2yTZs2zJ8//5z9H3/88Tn7/DvjLFq0iKeeeqrCawohREMhT4CNUO/evbn88st9A+GDyW63\nM2bMGN9YQSGEaKjkCTBADa3X0223hWaOAIvFwsSJE0NybSGECCd5AhRCCNEkSQIUQgjRJEkCFEII\nETZurVnVrFekwwAkAQohhAiTDSdtXPFRJqkXP8bqZhdFOhzpBCOEECJ4Ssb35U9d49v3Y46TJ7fk\n8MkRzyxVbW1ZFBsiP5GGJEAhhBAhcdbm5q87cnl9bwFODTEmxf/1iuOxl6YR67aRH+H4JAEKIYQI\nKrsy8vKefP62I5dsu0YBk1Ni+EO/BNrEGIl12yIdIiAJUAghRIAqqt70p7VmWfP+/O6CW9n/TQ4A\nw9tYmTEggd7NLeEKM2CSABuYN998k61bt/L5559z1VVX0bt3b6ZOnRrpsIRoMhrbkkDBcjDXyaNf\nZ7Oq18MAdEs0MWNAIiPPt/pWkalvpBdogOKmjCBuyoigXnPu3Lk8/PDDvtczZ87k7rvvrvLY5MmT\nefLJJxk2bBizZ8+uMPkVFRWRmprqmwotOTmZoUOHMnjwYKZMmUJhYSEAGRkZTJs2jT59+jB8+HDG\njx/P/v37fdf5+OOPSUpKYt++fVW+j8OHDzN48OAKj1U3aXZV5wbi2WefZc6cOQHdqzy73c7o0aMr\nXVtRCFG9Yqfmue25DPogg1XHbCQ4C5mdvpCNP2/Fte2j6m3yA0mAETVhwgRWrlxJdnY2K1asYOXK\nlbz00kvVHvNfnb0ib7/9Ntdffz1GoxEoXQ/wq6++wmKxMHfuXLTWTJ48maFDh7Jjxw7Wrl3Lk08+\nyalTp3zXWbJkCYMHD2bx4sW1fo+rVq2q9bmhvpfFYmH48OEsXbo0RBEJ0bitOlrMoA8yeG5HHjYX\n3HJBNN998wgPHF+J2VB/E18JSYARFBMTw7hx45g5cyaPPvooCxcuJDo6utpj1SXApUuXkpqaWuGx\nwYMHc/DgQdatW4fJZGLatGm+Y7169eKyyy4DID8/n7S0NObMmRNQgnC5XDz44IMMGjSIG2+8kaKi\nIgDatSutJpo1axaXXHIJo0aN4s477/Q9ublcLn7zm9+cc25lnn/+efr378+oUaPKPLGW3Ovw4cMM\nGDCAe++9l/79+3PXXXexZs0arr32Wvr168fWrVt954wZM4ZFixZV+/6EEKWO5DuZ9PkZbv7sDIfy\nXHRPMvHJ6Bb8e1gy59lzIh1ewCQBRtjkyZN5/fXXefbZZ+ncuXNAx3bv3k3Pnj0rvJ7dbufw4cN0\n7NjxnGNOp5PVq1fTo0cP9u7dW2USXb58OVdddRVdu3YlOTmZHTt2VPk+fvzxR375y1+SlpZGYmIi\ny5YtK3N827ZtLFu2jA0bNrBo0SK2b99e5typU6dWeq6/HTt2sHTpUtavX897773Htm3bKix34MAB\n7r//fjZv3sy+fftYtGgRK1asYMaMGbzwwgu+cj169Kj0GkI0BEnzjvvaJUPNpkw82+EGLl16ik+O\nFBNnUswckMC6sa0Ycp41LDEEk3SCgSrb9uJqULZE/oI1Ad971qxZtGjRosJ2qMqOvfHGG5Ve78yZ\nMyQmJpbZV1RUxNChQwHPE+Btt93G3Llzq4xryZIlvjUAb7rpJhYvXlxlwuzYsaNvdfo+ffpw5MiR\nMsfT0tJITU0lKiqKqKgoRo0aVebciy66qNJz/W3atIkxY8YQExMDwOjRoyuNp+SXhAsvvJDhw4ej\nlKJnz55lrm80GrFYLOTl5fkWGxZCVOy+btNY0GY4uDTjOkczc2AibWKMkQ6r1iQBRtCcOXMoLi5m\n3rx5PPvss9xwww0BHatKdHQ0xcXF5+zbsGFDmX3du3fnww8/rPAaZ8+eZd26dezZswelFG63G6UU\nM2bMqLRB22ot/e3PaDRWW40ZrHMDuabBYPC9Vkqds06izWYjKiqqzvcUojFzujVLWg4EYMnI5lzV\nruH/n5EESNVPbMXFxURFRfme/GrydFeVtWvX8tZbb7Fq1Sri4+N5/PHH2blzJ717967yWHWSkpJw\nu92+uCszbNgwnn76aebPn88dd9wBeKpWc3Nz2bdvH7fccgsvvviir3xqaiqbNm1iyJAhtXq/gwYN\n4qGHHuLhhx/G6XSycuVK331r4rLLLuO+++7zXWfFihV1GgaSlZVF8+bNMZsjPy2TEPXZt2cc5Jui\n6Vb4E1e1axxDQKQNMAKOHj3Kgw8+yPz5833Vbvfccw+vvPJKlccCNXz4cNLS0qoso5TizTffZM2a\nNfTp04dBgwbx9NNP07p1axYvXsx1111XpvwNN9zAkiVLavhOS/Xr14/Ro0czZMgQxo8fT48ePUhI\nSKjxdfr06cNNN93E0KFDGT9+PP369at1TADr16+v8fAJIZqiDSc9s7cMy/6+9hcpyMO8uv70ulZa\n60jHEDI5OTkVvrmcnJxz2skqE6onwFD65ptveP311/nPf/4T6VDKyM/PJy4ujsLCQlJTU3nxxRd9\n7YrVPbGGyuTJk3nqqafo2rVrhcer+reSnp5OSkpKKMMLOom57qobCB/ueIMxMD+QmG9efZpVx2ws\n/O5lbvjbXyotd853pdYYv9+Bae0nmLasRTkcvrLB+j5NTEys1ZgLqQJthHr37s3ll1+Oy+XyjQWs\nDx566CG+//57bDYbEydOrLJTTTjY7XbGjBlTafITQni43Jq0DDsAw7P3gtsNhqorEFX2GUzrV2Be\nvxxDRmkvVWfPSzDt2RLSeAMlCTBADeHJz99tt90W9GtmZWVV2Bln2bJlJCcnV3v+66+/HrZ7BcJi\nsTBx4sSgXEuIxmxXloNch2bhdy/Tzn4WPe0qiI5Fx8ShY+PRMXEQE+f50ytm+niU2w2AO7klzstT\ncVw+Ct2yTdBn1aotSYAiYMnJyef0Jm0M9xIi2BrbfKEbTtowul38/PRmAJTWUJiPKsyH0ycrPkkp\nnJcMwzFsDK5el4Ch/tRGlZAEKIQQokobTtrpk3+YGLen/S5/7mdQVIjyJkFVkOdLiFFz/wZA4exF\n6MTg1NaESth6gSqlRimlflBK7VdKPVbBcaWU+of3+E6lVD/v/vZKqS+VUt8ppfYopf7P75xkpdRq\npVS6989m4Xo/QgjRFLjcmq8ybJ62vxJGE8QloFu1xd2pG66e/XENGI5z+Bhfkfqe/CBMCVApZQRe\nBkYDPYCJSqke5YqNBlK8P3cDJf3+ncBvtNY9gEHAr/3OfQz4XGudAnzufS2EECJI9px1kGPXXJu/\nt/rCDUy4ngAHAvu11ge01nbgHWBsuTJjgYXaIw1IUkq10Vqf0FpvA9Ba5wF7gXZ+5yzwbi8Afh7q\nNyKEEPXRWZs7JNfdeNKO0e1i8NkfQnL9SApXG2A74Kjf62PApQGUaQecKNmhlOoE9AW+9u5qrbUu\nOX4SaF1ZAOnp6b5tk8mEUgqLJbAVistPLdYQSMx1Z7fbycrKKrNEVHn+/64aCom5rjzz0AYSU6jj\n9gzj9sQz56vDTGxX+7UtK4t15Y8W+uQfIcZRWG1Z8HxBB6NMdYIx1rLBdIJRSsUBS4CHtNa55Y9r\nrbVSqtJR/f4fltaa/Px8bDZbtffNzc2t1YwlkSQxB4fBYKBTp06Vzn9a3wZoB0JiDoINnh6elcXk\n/6Ue6rj3ZTsAzy9omwrieDKlZa2uU9ln7NaabzefYGp22erPQN5XsMqEUrgS4HGgvd/r8737Aiqj\nlDLjSX7/01r7z6OTUVJNqpRqQ8m/hGoopQKe+f/UqVO0b9+++oL1iMQsRNOw9kTpL/Fpp+z8VOCi\nbWzwhhvsPevkrE1zbV7ja/+D8LUBbgZSlFKdlVIWYAJQftG3ZcDt3t6gg4Acb2JTwBvAXq313ys4\nZ4p3ewpQ8fIGQgjRCK39qWwt1rLDdV9Jxd9G7/i/wdmNr/0PwpQAtdZO4H5gJZ5OLO9prfcope5R\nSt3jLbYcOADsB14D7vPuHwLcBlyplNrh/SlZ7vw54BqlVDpwtfe1EEKEnMsd2XmUXW7N+pNlE+CH\nh4KcADNsnvF/jkLcrdoG9dr1QdjaALXWy/EkOf99r/pta+DXFZy3AaiwEUZrfQa4KriRCiFExRx+\nSW/SF1m8PrwZcebILKrz7RnP8IQSViOkZdg5UegKyiK1Wms2nrRzu7f9z3VhHwynfqrzdesTWQ5J\nCCEC9NreAt/2iqPFpC4/zYlCVxVnVC5p3nHflGm14d/+B3BVuyg08FGQngJ/yHFyutjNyLzSBNjY\nSAIUQogAZBS6eG57aQf0LvFGdmY5uPqjTHZnOao4MzTWlGv/+3mnaAA+CFICLGn/u8zb/icJUAgh\nwmzAhpg6PSmVqOsT11Nbc8l1lFY5rr6uJYNaWThe6GL08kw+Oxa+cazFTs3Xp8omwFHto7Aa4asM\nOydr+VTqb2PJ/J/e9j/dvBWmEf/DNOJ/db52fSEJUAghqvF1ho239xdi8fvGbB5l5INrWzCuczR5\nDs0tn51h3vcFlV8kmPGcslPsgouSzb59CRYDV7b1VoPWsTeop/2vdP7Pxvj0B5IAhRCiSi635rdp\nOQA82Kvs+OEok+K14c14pHc8Lg3Tv8rmpYPmii4TVOtOeJ42h7exltn/887BqQb9MddJRpGba8q1\n/2VPbddolngCSYBCCOFTUTXp/H0F7MxycH6skYd7x51zjkEp/tg/gX8OTcKk4M3joU+AJe1/I9pa\nca6ZhHPNJMBTDWoxwKaTdjJqWA3q/943eOf/bMztfyAJUAghKnWm2MWMrZ6OL88MTCTGVPlX5uSU\nWJaMbEGcsbSdMNce/Amqs21utp9xYFIwuHXZ+YwTLQaubFf3atCNJz3j/2L92v8Clb9gDfkL1tT6\n3uEkCVAIISoxY2su2XbNiLZWbugYVW354W2tvHFxaWeY33+TE/SYNp604dYwoJWlwjGIte0NWvIk\nGY72v/qSJCUBCiEavLr28KzI9tN2FuwrxKRg1qWJlU6KXl6XmNInwP+lF/LpkeDOzrLGO/6vfPtf\nidEdvNWgGXZOFdW8N+ihPBc/Fbq5uhGP/yshCVAIIcpxa81v07LRwH094+iWVPt2vf/blM2Z4roP\nSyixztv+N7xtxQkw0WLginZRuHXtqkHXe8f/DWnk7X8gCVAIIc7xv/RCtmQ6OC/awG/7BLZyTEWG\nnGfhVJGb33yVg2e2x7o5UejihxwnsSZF/xaVr2fqqwY96EmANXlCrkv7X0MjCVAIIcr5s7fjy4wB\nicTXYa7Pl4c2I86k+OBQEUsP1r0qtGT1hyHnWbAYK6+SHd0+CrMBNmbYyaxhNejGk/ZGP/6vhCRA\nIYQo53Sxm8taW/hFl+g6XadTvIlnBiYC8Juvsms9b2iJkvk/h1XS/lciyWrgyrZWbzVo4DPUHIpq\nwbECF1fnSgIUQogmyahg1qCkgDu+VOX2bjFc085Ktl3zfxvP1roqVGvt1/5XfY/UsX69Qf3HClZl\nXWJ3T/tfTuNv/wNJgEKIJi6r2MW87wsY82mmb98vL4wtM81YXSil+MfQZiRZFKuO2fhvemHA5/q3\n3e3PdXK80EVzq4GezapfyS61QzRmA2w4aeOUOSGg+61N6t5k2v9AEqAQognKc7h598dCbl59mm7v\nnGT6V9lsPGn3Hf9Dv8ASRqDaxBh5fnASAI9/ncOhPGeNr7HWr/enIYAn0ySrgSu81aAftLgkoHus\nS+reZNr/QBKgECJEQjE2Lxju+DKLbm+f5FfrzrLqmA0NXNXOyr+GJvnKJFiC/9U4rnM0YztFke/U\n/HrD2Rqfv7aa8X8VKakGXdzq0mrLHrUmczC6FVc1kfY/kAQohGjktNZ86DcrygeHiihyaQa3tvD8\noER+mHAeS0a24NaU2JDGoZTi74OTaBllKPO0GQiXW7PuRNXj/yoyxlsNuiapB5nmqodzrEvytP8N\nbSLtfxBgAlRKWZVSzyilDiilcrz7Riql7g9teEII4fGX7bk17kW5NdNO6qenmfJllm/f05cksGt8\naz5Nbckvu8fRIsoY7FAr1TzKyEtDkqovWM7OLAc5dk3HOCOd4qtv/yuRZDUwoo0VtzJUWA3q1pr0\nHAfLTxn5V9trmlT7H0Cgn+RsoB0wCfjUu2+Pd/8/QxCXEEKUmcpr1o48Xvg2j+s6RvHLC+MYep6l\n0l6aR/KdzNiay6IDnie/FlEGThd7JqYuv6RRuKV2iObWrjG8td/TGcbu0lWO6YOy7X81NbZzNKuP\n21jc8lKuLHSx9bSdbZkOtpy2s+20nVy7BqyQmMLDRz4BmsbTHwSeAG8EumqtC5RSbgCt9XGlVONZ\nGEoIUe/8bUeeb3tspyg+PlzMh4c8P92TTNx5YSy3dI3xlcm1u5m9M49/fZePzQVWI9zXI47pvePp\n8L8TQYmpZDhB/tQ1tb7Gs5cm+hLguFWn+e+VzUmyVl4hV938n1UZ0yGah9yn+Ty5Fz979+Q5x9vG\nGOgWbeeazYv51fHPAUmA5dnLl1VKtQTOBD0iIYQADuY6mfdD6QrrC65ozk8FLubvK2DBDwXszXby\nSFqOb9YWgH5LMnxPeuO7RPNE/wQ6xAVeZRguiX6dbNaftDNqeSbvXt2cjpVUb6ZlBDYAviLNrAZu\nOr2Z91oNJt6s6NvCQv8WZvq1tNC/hYW2sUbS09Ppu+Qj3zmSAMtaBCxQSk0HUEq1AV4E3glVYEKI\npm3mtlyc5caMt4018njfBB7pHc/Hh4t47fsCvsoo7VByutjNoFYWnhmYSP+Wlc+VWZ90TzKxN9vJ\nNZ94kmDfCub4LHZBz2YmWkZX0V7pcoKx4q/0eXtf5ZkD75L8ynsYDVVXtzaV9j8IPAE+DvwV2AXE\nAOnAa8DTIYpLCNEElAyTyJ5atjVlx2k7Sw4WYTWCrYJ+Lxaj4qYuMdzUJYY9WQ6GfHgKgAVXJHND\nx6igzOASLp+mtuT2L7NYd8LGmE9P88bwZozu4Bm+4FwziT92Hs9zHX9ecfuf36wysb8ciW7eGnfL\nNuiWbT1/tmqDu0UbrNpJ5+JM8g0K3G5wu8Dl8m2bCkqrmpvK0x8EmAC11nZgOjDdW/V5WgdjanMh\nhKhASbXm3d3jmLM7v8qyPf1mbCkZ99aQJFkNLL6mOf+3KZu39xcy6YssnhuYyN094gD4otlFAIxo\nc+70Z6aNK33byu1GZZ7AkHkC2FbhvWLvuAJVwVd3L79tSYAVUErFAF2BOCCl5DcsrfWm0IQmhGiK\nvjxezJc/2UiwKB7uHV9tAgym8k+itVXTjjIWo+JfQ5PoFG/k2e15/O7rHA7lO3nKFMuW+C6YFFx2\nXtmqUZV9ButbL/te57+2EnX6JIbMkxgyf/IlQ5V5AuOR/Z5zvMlPG01gMIDRCAYjTsBc6PmcnX0G\n1/HdNxwBJUCl1O14hjvYAf81PTTQIQRxCSGaILfWPOV9+pveK55mVfSMbGyUUjzax9Np58GNZ/nX\nngK+7Psn3MrAwJYW4vyXZdIa68IXUX5Vl1is6LYdcbXtSPla47gpIwDIn/eFJ/GVk56eTt+Zd3le\nxEZ2mEg4BfoEOAsYp7VeHcpghBBN2wcHi/j2jIM2MQZ+1SO0M7PUVxO7xtA2xshtX55hL54n0mHl\n2v+Mm9di2roeHRWDKg58cu2Kkl9TFuinYQfWhDAOIUQTZ3dpZmzzPP091ieBGFPT/bIe3tbKqjEt\n6VjsWaFidHu/9r/8HKz/fQkA2y2/ikR4jUag/8KeAP6ulGoRymCEEE3Xwn0FHMxzkZJoYlJKTPUn\nRFqI+wFemGRm2+bfs2XL42WGRlj/908MuWdxXXgxzhHXhzSGxi7QBLgPuAHIUEq5vD9upVTAE/Mp\npUYppX5QSu1XSj1WwXGllPqH9/hOpVQ/v2NzlVKnlFK7y53zlFLquFJqh/cnNdB4hBD1R77DzV+9\ns7480S8BUzVj1eqDqJefAv82uBBIdBXRJ/+w77Vxx1eYN61GW6wUT/utVGnWUaCf3n+BhcDFQDfv\nT4r3z2oppYzAy8BooAcwUSnVo1yx0d5rpgB3A6/4HZsPjKrk8rO11n28P8sDejdCiHrl5T35ZBa7\nuaSlmes7Vr/aeaQYDnzv2zZtXkvME7/EsG9XeG5emI91/gsA2MfdiW59fnju24gFmgCbA09qrXdr\nrX/0/wnw/IHAfq31Ae+YwneAseXKjAUWao80IMk74wxa63VAFkKIBieruPqKojm7PF3wn7oksV4P\nYrd88pZv29WlO4YzGUT/5f8wf7DAM7g8hKzv/hvD2dO4LuiOY+S4kN6rqQg0Ac4DbqvDfdoBR/1e\nH/Puq2mZijzgrTKdq5RqVocYhRAh0OXtkwx6P4Ppm87y3o+FHM0/dzX0fKdm5PlWhp5X87kuw0Wd\nOIJx63rf66I/zME+5lZAY31/HtHPTkedyQjJvY3fbcO85iO00YRt2u/AEL4lnBqzQIdBDATuV0r9\nASjzN6y1Hhb0qAL3CjADz3jEGcALwLSKCqanp9f6JnU5N1Ik5vCQmCvm6R/i6chiVprvs518n+1k\n3g+eLvutrW76JLgp+QpSaO5omU16ekUrpccEEHf1ZTYPJYDrVK79xwuI9ev4kn7wIPS9gvjE1nRc\nNhfzvp1e03gbAAAgAElEQVRYH59K4pjbS8vUMea+3j8N//4LACeGjiGjyAl+55SUCeQ6gbz3hvJv\nOiUlpc7XCDQBvub9qa3jQHu/1+d799W0TBlaa18yVkq9BnxcWdnafljp6elB+aDDSWIOD4m5ct9n\nOwDP/JzHbmvHjtN2vsqw81WGjbRTdjJsBlZmllZA3XxBDKl9KmnT2uD5Gqgy7kDKBMA3YHzBmjL7\nVVYmMbvS0MqA0u6y90pJwTZ4OOr15zB9m0aXJa/6zgtWzNbs07g6XEDC7Q+QYKr4azuQ61RWxj/p\nNbR/03UR6FygC+p4n814pk/rjCepTQBuLVdmGZ6nzHeAS4EcrXWVC3gppdr4lbkR2F1VeSFEeKzz\nLuAKYDUqLm1t5dLWVh4iHrfW7D3r5KsMG4+k5QDwp/6JlV5r89DCiH8pm1cuQrmcOAZegfmbL88t\nkJBE8fRnMa9eivmdVzC4vNW8udmQUPMV4MvTBgO2Ox+FSpKfqJ1KP02l1G1a6/96tyusVgTQWs+t\n7iZaa6dS6n5gJWAE5mqt9yil7vEefxVYDqQC+4FCYKpfLG8DI4AWSqljwJ+01m8As5RSffBUgR4C\nZFSoEPXA+pO2So8ZlKJnspmeyWZfAmwbW4/btPJzMa/xrJXnGDOx4gQIoBSOkeM4GNOMC1/zLJRj\nffdVbHedM+orME6Hb9MxegLuTgF1uq9Q+Sda4VHVrxMT8Qx/gMo7wGig2gQI4B2isLzcvlf9tjXw\n60rOnVjJ/rp0zBFChIBba9afqDwBhkKwJrGuiPnzD1DFRTgvGhBQEipqXdqSY96wAsfwVNzdetf8\nvstLl1u1/3xKjc8X1as0AWqtU/22rwhPOEKIhm5XloNseyNZLc1WjGX1EgAc15VvtQmMdcFsiv78\n2jnVl1WtGKEyjmFZtrB0h6X+9o5tyAIaBqGU2l7J/i3BDUcI0dCF++kvlMzrlqPycnB16V6rdfLc\nLdtiPHYQ82dLAz9Ja6wLZqMcjurLijoJdBxg1/I7lGe0apfghiOECKUBG2J8q7CHSqNJgE4n5hXv\nAmAfMxFqMUDfdtuDAFjen4fKygzoHNNXn2HasxUdm1Dj+4maqbJLkVKq5Bnc4rddohOwJxRBCSEa\nJodbs/GkPdJhBIXp6y8wnM7A3aY9rn5Dfftr0qHEdfEgnP0vx7R1PZa3/4Xt13+q+oT8XCzeRW5t\nE+4h6o1ZtQldBKi6PrU/VrKtgY3AoqBHJISIqJInxNp0LNlx2kG+U9M1wcT+3HNnfGkw3G7M3mnP\n7KkT6zTptG3S/Rh3bcb8zZc4h6Xi6jWg0rLW9/6NIS8b188uxnn5aAggAQarh2dT7ClaZQLUWv8Z\nQCmVprVeGZ6QhBANVcnwh2FtrA06ARp3pmE8fgh3sxY4B19dp2vp5q2x//x2rO/9B+t/X6Lwmblg\ntpxTzvDDTsxrP0EbTRTf8XCtqlxFzQT0a40kPyFEINZ52/8ub3PuF3xDYvnY8/TnGHVzhcmqphzX\njsfdtiOGjGNlhjf4OB1EeVd6cFx3K7ptxzrfU1RPFpMSQgSFzaVJy/AkwPo8qXUgjOm70bHxOIZf\nF5wLmszYbn8IAMtHb6Iyy05yZV7+DoafDuNu3Q77dZOCc09RLUmAQjQSSfOOh7yHZ1W2ZNopdkGP\nZiZaRtfjmV0C5Lj6RogO3sr0ru59cQy+GuWwY/3vS779/mP+bFMeljF/YSQJUIgGINjJTevgD1T3\nVX9G4OkvbsoI30TWwaAtVuxX3xS065WwT7gXHR2L6ds0376SMX+Oy0bi6tk/6PcUlQtoZtUq5gK1\n4Vm3L01r3UgG/wjR+BzMdfLPPfm+189sy+OP/YM7zqwkAQ5rE3gCDOUUZjVWXOjbdAwfE5RJrMvT\nSc2xj7sT65v/8O0rGfNnm3hf0O8nqhbo1OK3A4PxrAV4DM9SRa2BLXjGA6KUGqu1lplhhKhHdpy2\n89KufD48XITb76Hv+Z159G5u5oZO0UG5T6HTzZZMOwYFQ+pp+19lSx15duYS/fdHfS8dqRVOPxyw\nqoYUOK68AdO65RiP7Pfts024JyQJV1Qt0CrQPcBvtdYdtNaXaa07AL8BtuNJhq8Ac0IUoxCiBrTW\nfH68mBtWnGbER5m8f6gIo4Jbu5Ztz7p3/Vn2ng3OdFtfZ9hxuKF3spkkq+drJXtqu/r1hFcJlX2G\n6GcfwvjjXt8+ndwydDc0mrBNme576RvzJ8Iu0CfAW4Hm5fa9ApzWWt+vlPob8NugRiZEE1GXgecV\nuXxZJruzPIkt3qy442ex3NMjjnaxRt7a76nmG98lmkUHipj0+Rm+uL6VL2nVVm2qP+sDlXmC6Fm/\nwXDqJ9xtOmA4cSQs93V37enbrmrMX1McnB5Ogf6rzwCuL7dvDCVLPkMUIDO3ChFmWmu2n7bzu7Rs\n377dWQ5aRRv4U/8Edo0/jxkDEmlXbr29l4Yk0SvZzIE8F3etzcLlrlunmIaYANVPh4l+5gEMp37C\n1bEbhY//o/qTQkDG/EVOoE+ADwKLlFK7gaNAe+AiYLz3+KVIFagQYXMs38miA0W8s7+QH3LKzrjy\njyFJ3NwlhihT5TOJxJgMvHllMld8lMnq4zb+sj2XJ6pYlb0qOXY32884MCkY1LphDIA3HNpH9PO/\n9az00K03RdP/AjFxkQ5LhFlACVBrvUop1QXPiu1t8Sxs+4nW+kzJcWBVyKIUQgDw9v5C3tlfyLoT\nNkqe2VpEGfhFl2he/a4AgNu7xQZ0rY7xJuaNaMZNq87wws58eje3MLYWnWK+yrDh1jCglYU4c/0f\nWWX4YSfRs3+PKirA2ftSiu//M1ijIh2WiIBAnwDxJrv/VltQCBEy964/C4DVCKnto5nQNYYr21kx\nG5QvAdbE8LZRPD0gkT98k8N968+SkhjwV4JP6fRnDaP6M/r536LsNhwDr8D2q8fBZI50SCJCAh0H\n2Bl4BugDlKkn8PYIFUKEyE8FLt/2oFYWJnaNYWyn6Dp3XClxX49Yvj1t5z1vp5hA+HfcWX/Cs/xR\nbQbAVzk0wavvzLuqLVMtZ2k1sbLbcAwfg+2Oh8HQ8GesEbUX6K97b+FZDuk3QGE1ZYUQQfT7b0o7\nuKwYU7fu+ZuHFpKSklJmn1KKl4Y044ccJ9+eqVlftqxiF7uyHFiNMLBVPWz/0xrjtg1YF/3Ht8s+\n6mbsE+6V1RZEwAmwJzBEa+0OZTBCiLJWHS3mw0PFIb9PtEnx5pXJjFiWyRlb4P/N13sXvx3Y0kJ0\nFZ1uIsGQvhvrO69i3L+7zH5JfqJEoHUo64C+oQxECFFWodPNI37DG0KtfZyJ+Vck+17P2ZVX7Zyh\nG+ph+586cYSoOU8SM/N+jPt3445Pwjb5Qb8CkvyER6BPgIeAFUqp94GT/ge01k8GOyghBMzakceR\nfBcXJZt9A9urEoyB9P6J7IktuRwpcPHcwESMhoqTRn0b/2ddMBvTmo9QbjfaYsUx6mbsqRMgOrbM\n/JtCQOAJMBb4GDDjGQMohAih7846+OfufBTw4mVJXP1xZthjsBjgtb0FHMt38frwZsRWMMThhxwn\nsSZFvxYRbP/z6+Bi/uJDtDLgGH4d9hvvQDdrEbm4RL0X6DjAqaEORAjh4daa6ZuycWq488JYLmkZ\nmeTy/rUtmPT5GT49Wsz1K07zztXNaVXBOn+DWluwGCNTrWg4+APWuX/zvXb2uQzbzXej23WKSDyi\nYak0ASqlOmmtD3m3u1RWTmt9IARxCdFk/XdfIV+fstM62sAT/YK7ZFFNDDnPyqoxLRm/+gzbTju4\n5uNMFo9sTkpi2XFzEan+tBVjeX8e5hWLUH5984qn/yX8sYgGq6onwF1AvHd7P6CB8r/maUAG0ggR\nJJlFLv60JQeAvwxMDNpYv9rqlmRm9XUtmfBZaRJ866qy8+KHOwEa92zFOu8FDJk/oZUB+6ibsax4\nLyjXDufk0zLRdeRVmgC11vF+2/V/fiMhGoE/bM4h2665sq2VmzoHZ62+umoVbeSjUS24c+1ZVhwt\n5ucrT/uOJVgUvZPDNJNKfi7Wd17BvP5TAFznd8E27be4L+geUAKUhCPKq/m8R0IIAAZsiIENx4O2\njNHan4p578cioozwwuAklF93/UivqxdrNvC/K5N57OscXvu+dMq1Ia2tlfYQDbaY30/BkHsWbTZj\nHzsFx+gJYJKvMFF7MhWaEPXEw195xvw9cnECnRPq3xe70aCYNSiRDvFGnticC4Rg/J/WqKxTGA6n\nYzyUjuFwuu+QIfcsrm69KZ72CLqNfO2IupOp0ISoQLAXqQ3Ej7kufpZo4sGL6u+yPEopHrgo3pcA\ny68yX1uWd/+N4fA+jIfTUfm5FZYpvuNhnMOvA4O0yIjgCNtUaEqpUcBLeDrNvK61fq7cceU9noon\nyd6htd7mPTYXuA44pbW+yO+cZOBdoBOewfo3a63P1jZGISLt75clRWxIQW3UqZOO3ebbtCx/27et\n4xJwdeyGu2MK7k4pRP3raQCcV9xQ+3sJUYFAE2DJVGhba3MTpZQReBm4BjgGbFZKLdNaf+dXbDSQ\n4v25FHjF+yfAfOCfwMJyl34M+Fxr/ZxS6jHv60drE6MQkTYpJYYhtVhRob6qcqUHt5uofz/je2kf\nOwVXpxTcHVPQya3KTlfmTYBCBFu4pkIbCOwvGTOolHoHGAv4J8CxwELtmXwwTSmVpJRqo7U+obVe\np5TqVMF1xwIjvNsLgDVIAhQNyOfHSye6/vMlkRvzF26WRa9h2rLO99p+k8y1IcIvXFOhtQOO+r0+\nRunTXVVl2gEnqrhua611yfGTQOtaxCaamEi071XE6db88Zsc3+sWUU1jSK3py4+wLH8bbTSiXK7q\nTxAiRBrNVGhaa62UqnTq+vT09MoOVasu50ZKU415wAZPp4zNQ6vqqxUTwP2CVaZyS08a2ZtdWuUZ\nrr+z6u8TnPdesnyMf5n4H3dzwTtzADgyahIdP1lYq+vUpkw41Zc4aqKhxVx+XcvaCLivtVIqEfgZ\n5w6D+CKA049T9snxfO++mpYpL6OkmlQp1QY4VVnB2n5Y6enpQfmgw6lJx7zB80+mymuFs0wlcu1u\nXt+SAZT2KwvH31lAn3MA78s5cwQA+VPXVHvPkusYjvxI9PuvobQb+3WTaD5+GngTYCDvPVhlQq1J\n//9rYALqwqWUugP4CfgIeMPv5/UA77MZSFFKdVZKWYAJwLJyZZYBtyuPQUCOX/VmZZYBU7zbU4AP\nA4xHiFqzuTSfHSttu9sVwFJF5b20K4/MYjcDIzTRdbips6eJmv0YqrgQx6VXYB93Z6RDEiLgJ8Bn\ngF9orT+tzU201k6l1P3ASjzDIOZqrfcope7xHn8VWI5nCMR+PMMgfNWuSqm38XR2aaGUOgb8SWv9\nBvAc8J5S6k7gMHBzbeITojrZNjerjhWz/Egxnx8vJs9RWtt+8+rTfHZdK9rFBtaGdyTfyT/35APw\nzMBErvkk/EsdhVVxIVGzH8eQlYmr60XYfvmYjOUT9UKgCdAErKrLjbTWy/EkOf99r/pta+DXlZw7\nsZL9Z4Cr6hKXEFX593f5LD9SzMaTNpx+Lcw9m5nYc9azDt2JQje3fHaGT1NbEF/Bmnnlzdiai80F\n4zpHM6BV438CjHplJsbD+3C3akvRQzPB0niGeoiGLdAE+Ffgj0qpGXUZDC9EXYSr9+bR/NIFVh/9\n2tNL06jg8vMspHaIZnSHKDrFm3zxXJBgZHeWg6lfZvHO1c0xVTE35tZMO4sOFGE1wpP96+ewh2B/\nvqYdm9Cx8RT95q8QnxTUa5eQia5FbQSaAKcD5wG/U0qd8T8gc4GKxqLYqfnnnnxe+DbPt29spyhS\nO0Qz8vwomlUy68mia1pwzceZfHbcxu/ScnhhcKJvImv/pK215g/eYQ/39YijY3z9m+8zaHKzfZva\nZKbowZno82ozgkqI0An0f+DkkEYhRIStPFrMY19nczCv7Li0BVc0r+SMUl0STLx1VTI3rDzN3B8K\n6Bxv5IFe8eeUW3a4mLRTdlpEGZje+9zj4dB35l1ACJ+YCvKwrHgP86rFvl22O3+H+8KLa33J7X98\nrUn2UBShF+g4wLWhDkQ0bcFeWihQB3Kd/P7rbFYe88xLeWGSiVmDkrhhxelqzizr0tZWXr28GVPX\nnOWJLbl0iDcxtlPpen42l+bJzZ6nv8f7JpBgaWSdQIoLMa9eimX5O6jC/DKHnJddE6GghKhaoMsh\nVToZX4BToQlR78zcmss/dudhd0OCWfFY3wTu6h6LuZbr293YOYZDeS7+vDWXX63Lom1MS9+x/+zN\n53C+iwuTTNzeLTgrKNQLdhvmL5dh/vgtDLmeeeidF/bBPu5OYp55IMLBCVG1QKtAy1fenwcMB94P\nbjhC1N7xgtLqy0uXZmA0gEkpTN4/Pa9Lyz+/09PWN7FrDE/1T6B1TN2nInuoVxwH85ws3FfIxM9L\nm8v/5m1XnDEgscpOMg1NzO8mYTjreVp2demO/Rd34urRv+xk1lWQzisikmo9FZp3eaMKhycIEU47\nTtv51558lh4s8u37IcdZxRkeFzc387dBiQxsFbxu+UopXhicxLF8F1/8VLrcT65dc2VbK1e3a1xD\nAAxnT+NqfwH2cXfi6jM44MQnRH1Ql25oq/CsxSdE2Lm1ZuXRYl7ek8+Gk3YA/B+s0m5shdPtmXDa\npT1/OjU43TB2peeJ5YvrWmIMwdOY2aCYd0Uyoz/J5Ltspy+2GQNKe4eGQpXLDwWTrfQXjeK7H8c5\n+GoZ2C4apEDbALuU2xUD3ErZ1RuECIu53xfwrz357M/1JJd4s+L2brHc3T2WixdnAHBhkrna64Qi\n+ZVItBh495rm9Frkief2lBh6JlcfU0Ng2rLet+0cMjKCkQhRN4E+Ae4HNFDyjVEIbKd0Hk4hQirb\nVjr/wsNfecaYnR9r5J4esdzeLTYivSo3Dy2ssnt++7jS/15PXZIYjpDCwrRhRaRDECIoAk2AZq21\nLNwlwi6r2MW/vivgP9+Vdq3v38LM/RfFcX3H6AbToSSpkkH0DY06k4Fx7/ZIhyFEUFSbAJVSRiBf\nKZWktbZVV16IYMgscvHynnxe31tAvrPsMo+fXdcypG1ponKmjatQutJlN4VoUKpNgFprl1JqH9Ac\nz5JIQoTMyUIXc3bnM/f7Aopcni/aq9tZ+e3F8Vy73NN5JRzJL5wD8iO9Mn3AtMa8YWXAxWWIg6jv\nAq0C/R/wsVLqJeAYnvZAIOAFcYUIyMWLT2LzVraPah/Fby+Op38TWTMvpOw2TNs3+V4ajvyIu8MF\nNbqEYf8eDBnHcCc1x5B9pvoThKjnAk2A93r/fKrcfg2U7yEqRK3ZXHB9xygeuTiei5tL4qsTrTGk\n78a8cRWmb75AFRb4Dlnefpni371Qo3F75vWezi/Oy0ZiWf520MMVItwCHQjfOdSBiKbLf/mhL65r\nST954jtHTcb4qcwTmDatxrxxJYaM4779rk7dMB7aB4Dpu20Yd2/B1WtAYAHYbZi++RIAx9BrJQGK\nRqERr8ciGop/7i7t4dnYkl8k2vdiHymdoMmd1ALnZdfgHDIS9/mdfYkUwLLoPxT17B/QIHbT1g2o\nogJcXbqj23UKQdRChJ8kQBFRZ4pdLNxXGLTrNZgOJUFm+P5b37a2WHH2vxznkGtx9ewHhnPnOHUn\ntcB4OB3TN1/iHHRVtdcvGfvnGHotIB1cROMgCVBE1L/3lvb2DJfGliSNe7cT9fff+14X/GMpRMdW\neY79xjuImvc8liVv4LxkGJhKZ6kpX92qsk5h3LMFbTLjvPTKoMcvRKQ0jtG5okHKd7jLDHAXNWfc\ns5Wovz+GsheX7qwm+QE4Lx+Fu017DKd+wrT2kyrLmjauRmmNs+8QiEuoa8hC1BuSAEXEzP+hgGy7\nZlCrxtXuFy7G3VuImv17lN2GY1hqDU82YRv3SwAsHy4oM8F1GVpj3ujt/emt/hSisZAEKCLC5tK8\nvMfz9PdQ77iAzsme2q7RVV+Cp8rRv3NKIIw7vybqxd+jHHYcI67HNvWRGt/XdckwXF26Y8g5i3nl\n4grLGH78DsOJo7gTmwXeY1SIBkISoIiI934s5EShmx5JJq49PyrS4TQoxm/TiHrpjyiHA8eVY7FN\nmV675YiUwn7z3QBYlr8D+TnnFDFvKB37h1G6DIjGRRKgCDuXW/OP3SVPf/Eyr2cNGLdvIuofT6Cc\nDuxX34jt9ofqtBafq3tfnL0GoIoKsHz0v7IH7TZMX3smepLqT9EYSQIUYffxkWLSc5x0iDNyU+fo\nSIfToETNedKT/EaOwz75waCswG4f73kKNH/2PupMhm+/aftGVGEBrk7dcJ8vEz6JxkcSoAgrrTUv\n7soD4IGL4hrMckYRY7dh3LXZ91K5nNhH3Yz91vuDkvwA3B1TcAy6CuV0YFk6z7ffVDL12dBRQbmP\nEPWNVOqLkEua55mOK3tqO9adsLH9tIMWUQYmp5R2169ucdkmQ2tUxnFMu77BuOsbjHu3o+ylq5DZ\nUyd62u1qmfy2//G1Cj9n+03TMG1eg2njKt8+4+4taKMJx+DqB8oL0RBJAhRhNXuXp+3v3h5xRJvk\n6c+fZeGLmHZ+gyGz7Kpjro4pGA+nA2C/5Vchubdu3Q7HiOuxfP6Bb5/Sbpx9h0Jc41nNXgh/kgCb\nKP+nsnDZftrOmp9sxJsVd15Y/WDtpqYk+ejYBJwXXYKr96W4LroEndQ8oGESdZ2ezDH2dswbVqBs\npYPqHVL9KRoxSYCiTmqSSGfv9LT9Tf1ZLElWaX4GUCeP+rbtY6fgvPhS3J1/VuH8naGmE5NxjLoZ\ny4cLAXAnNMPVa2DY4xAiXORbSITNR4eLsRjgvp6BDXxvCixL5/q27TdNxX1Bj4gkP18Mo2/xbTsv\nuwZM8juyaLzClgCVUqOUUj8opfYrpR6r4LhSSv3De3ynUqpfdecqpZ5SSh1XSu3w/tRwPigRThq4\ntWsM58VE7gu+PjEc2of56y8jHUZZfvOIOkb+IoKBCBF6YUmASikj8DIwGugBTFRK9ShXbDSQ4v25\nG3glwHNna637eH+Wh/adiLowKHiwV3ykw6g3LItfj3QIVdLNW0U6BCFCKlxPgAOB/VrrA1prO/AO\nMLZcmbHAQu2RBiQppdoEeK7wM2BDjK9trjq2KpYiSpp3PODrVCTf4WaOd8wfwNiO0XRJkCo18Kzf\nZ9r1DToqJtKhCNFkhSsBtgOO+r0+5t0XSJnqzn3AW2U6VynVLHghN045drevMwpA/yUZzP+hAHsQ\n1+TLtrmZtSOXXotO8sSWXN/+Pw+QpXQA0BrroteAsm1uQojwaui/jr8CzMDTvDQDeAGYVlHB9PT0\nWt+kLucG24ANnieGzUOrWkXdU8Y/7jN2ePsnM4tPmChwlY6/O1bg4qFN2czamsW09g7GtHJhMlR+\nnaruddYBbx03s8jvHr3jXezM87T52U4cpKor1afPue/MuwDPwPGq1CbmhPRvuWD/bhwx8XzftR8X\nM6/aa/Wtw/3Kq+4awbxXsNSnWALR0OKFhhdzMCbOCFcCPA6093t9vndfIGXMlZ2rtfZNXKiUeg34\nuLIAavthpaen168ZSjZ4PrYqY/IrcyjPyT935/NmegHFLs/hy8+zsP6kHYC5w5vx1x15/JDjZOZ+\nK29mGHnk4ngmXBADG34K+F5zs1oy/4dC3+ruw9tYeeTieIaeZ6HZ/OqvU+8+Z6+6xlx+dXXcbqLn\nP+vZ/PntXHBRL9+xQN59XT+jmnzO9eXvo77+26hMQ4sXGmbMwRCuBLgZSFFKdcaTvCYAt5Yrswy4\nXyn1DnApkKO1PqGUyqzsXKVUG631Ce/5NwK7Q/9WGo6712ax5GARJbWbqR2imN4rngGtLL62vZu6\nxDC2UzTvHyrirzvySM9xcv+GbF74Nq/MtQqdbo7kuziU5+RQXumfJV75rgCAa9tH8Uhvzz0au5Kn\nxJoMQDelfY7x2AHcya1wXHFDiCITQgQiLAlQa+1USt0PrASMwFyt9R6l1D3e468Cy4FUYD9QCEyt\n6lzvpWcppfrgqQI9BIRmnqgGwOXWbMm0s+pY6Swe7x0owqhgwgXRPNQ7nguTzBWeazQoftElhhs7\nRbP4YBF/3Z7LAb/kduE7JzhZ5K7y/j/vFM3DvePo3bzxJ75a85ts2n7jHWCxBnxqXWd5EUKcK2xt\ngN4hCsvL7XvVb1sDvw70XO/+24IcZoOSVezis+M2Vh8r5rPjxZy1le3Iclf3WB64KI4OcYH9NRsN\nilsuiGFc52je+7GQ+zZkA3CyyI3ZAB3ijHSKN9ExzkSneCMd401M+TILgPlXJAf3zTVCprXLMWT+\nhLtNB5xDRkY6HCGavIbeCabJuvaTTDZn2nH75bzO8UZGnh/Fv/d6qiP/NiipVtc2GRS3psT6EuCu\n8a1pG2PEKEsX1Z6tCMuHCzyb4+6U1dWFqAfkf2EDUezUvPJdvu/116fsmA0wrI2VkedHMfJ8K10T\nPVWcJQkwWNoH+ARZmXBOuH1Op5N6wrx6KYacLFydf4brkmGRDkcIgSTAeqf85NJaa5YcLOLPW3M5\nml/aLvfmlckMb2sl3ly7oZzhTEr1TU2SpCntc5w9+kNC7Z6mS1g+eRsA+y/uCtpCtkKIupEEWI99\nnWHjD5tz2JLpAKBHMxPfnXUCcF3H6EiGFhHheLpTJ45gfe/fvtdRr8wAwNWhK66e/XH1vARXt15g\njarZdQvzcXbvi6tn/6DGGwr17elZiFCRBFhPTfnyDB8e8vTobB1t4A/9EpjUNYbmC36q5kxRGyon\nC8sHCzCt+QjlLu3x6uzZH+O+nRiP7Md4ZD98+i7abMaV0stXxnBon2dD+3dCOndmHfsvfilPf0LU\nI5IAw6i6tfMyi0qrOD88VEy0UXH/RXH8X6844mpQ1bl5aGGTHNRaK7ZizCvew7L8bVRxEVoZcAy/\nDpYQkQUAABGmSURBVPNaz5wKxb97Aew2jOm7MO7ZinH3VgxH0jF9t813iZg/3V3tbZz9huLu2jNk\nb0MIUXOSAOuBb8/YefW7ApYcKJ3ebMIF0TzRP5F2sZFbOqghtxMaDnyPjo1HxyV4lvgxnPsLhGnd\ncixL5mLIPg2As89gbDf/Ct2uky8BAmCxeqo+e14CNwN52Rj37iD65acAcHX0/2VDldk0ep8ObRPu\nDfI7FELUlSTACHG6NZ8cKebV7/L5KsMzLZl/5dirwxrGuLr6kCRVxnFMOzZh3L7Jty/mz/f4trVS\nEBOPjotHx5ZOyB31xiwAXJ26YZ9wL67ufQlIfBKugSM8i3QBRU9XPl9oSbulbh35z0kIUZYkwCCp\nrnrT30u78nhtbwHHCjxVnglmxaSUGO7uHkffJRlVnlsfEk7EuV0YDnyPadtGjDs2YTx+6Jwirk7d\nUAV5qIJcVGEBFOSiCnLxn4LW3aI19l/chfPSKyt8QhRCNG6SAMMkPcfh2/6Td4mgCxKM/Kp7HBNT\nYmo9nKGxqHZeTb8OJjH/9wsMuWdLD0XH4ux9Ka6+lxH16kwAiv78n9JzXU4ozPckxPxcYmZ4Jhwq\nfHZhjaYjE0I0LpIAA1CTp7vytp+2M3tnHh8dLp2j88q2Vu7pEcfV51sxSK/AKqmcLEwbVpZpkzPk\nnsXd4jycfYfg6nsZrp/1BpN3nlNvAizDaIL4JHR8Utm+mZL8hGjSJAGGgNaadSdszN6Vz5qfbABY\nDGD39q5fem2LSs+VKk7A7ca4ZyvmNR9h3L4R5XKVOVw4cy7u8zvLkAIhRJ1IAgyyZYeKeHFXHttO\ne6o840yKaRfGcm/POLq/ezLC0dV/5mX/xbz2EwynPZ+VVgacfYfgGDGG6NmPA+Bu3yXkcchgcCEa\nP0mAtaS1Js+hybK5OVNcOnD6du/qCC2iDNzbI447L4wlydq02/eq5bD7Nq1L3gA8HVQcw8bgvHw0\nOrllpCITQjRikgBr4PpPMzljc5NV7OaMzY2jgiXy2scZefCiOCanxBJtkiq66hiOHcT679J2O+cl\nw3CMuM4zZZghcmMghRCNnyTAGlh/0l7mdaxJkRxloLnVwI4znirPbeNaY5Zlg6qnNebP3sfy7qso\nvyfA4geertNlA6m6DFb15vY/viYz7gjRgEkCrIEPrm1OstVA8ygjyVZDmSe8kp6iVSU/6eDiobLP\nYH39r5h2fQOAY1gq5nXnrHfcKEhbohD1lyTAGhjRtmYrAAgP/1UcjNs2EDX3b6i8HHRsAsXTHsF1\nybBGmwCFEPWXJMAgca6ZBED+1DWRDaQes8593jeez9nzEmx3PYZuVvmQkPLkaUoIEUySAJuoQNbW\nC0YZdfa0b9u89mO02Yz95l/huPqmMtOPSXuaECLcJAEGg91Wuq110x6grTWGI/sxbt/kmaD64A++\nQ67zu2C7549hGccnhBDVkQQYgOo6r1jn/s23HTN9PK7el+K8eBCuHv0hOsZ3LBhPVIFco8p5NYsK\nMBz90ffS8t5/0FHREBWDjopBR3u3raUrzqusTLTJDEajZ8oxk8kzRMEv0Ru/TcO0fRPGb7/CkJXp\n268tVpT3F4SiP70i048JIeoNSYB1ZPrqc8xffeZ7bTh7GsPaTzCv/QRtNOH6WW9cFw/C2fvSsMem\nMo5hOPIjxqMHMBzdj+HIj74ZVv6/vXsPjqss4zj+/eXaS0rCTbTQhooRreiUOiBFRpQCUlCqgzpU\ny00YRKUKjmABB5HxAo4woqIdQVSUwpSiglfKTRkZAREKlBZMS1tKKbRUSS9p0ib7+Mf7bjgsu5tt\n2GzOdp/PzJk9+573nH3ONs2T8+7Z98lq+tP8QY8z9vxP5m237PybwOir5w6sZ9r2pH/K4fQdPI3+\nd06l5ezj4ot58nPOpYcnwDdAG1+i+carX9PWffl14Wro8YeoW7GUhqWP0rD0UZpv/slAn1FXnA/N\no7CmUfGxGZpHh8eocdFCUB1ImOrC52WJK67GO29F3VtClYPEwtYtA33GXjj7dTFbYyOZ8ZOoXx0L\ntX7iLLStG3q6Uc821LsNtnWjnm7qVywFQkJTfx/09YXKCn07UCaD+l6tcNHf3jGQ9DLtb/fyQs65\n1PMEOFSZDM0/+y7q3krflMNpWByKsWbaO8i0d7DjxFNgSxcNTz5C/RMP0vDkw2hzFwANyx4b9PDN\nN/24+Pb51w4eYtueZCYcEJaJ8fHNE6ChYWAodcdHX58ks7J9uq+5Lc/BM9DfR8tZxwLFi8L63ZvO\nuTTyBDhEjX9dQMPTi8nstju9Z15Aw5yPv75TSyt906bTN206vZl+Ws6YDsC2C6+C3h60vSc+9kLv\nNrS9l6bbbwRg+zEngWXADGX6IWNgmYHvy20/5iQY04IlFsaGxzGXhs8A8yaucqmrg7qm4Tu+c84N\nM0+AQ1D33HKa4qTNvWdeiO22++BXOYl5Lfvf9d6C3QYS4Ow5ebcPJMAC20vlV2XOuVrnCXBnbe+l\ned63UN8OdnzoRPqnTBvpiF6nXN+pq+S8ms45V2meAHdS063XUb92FZk3T6B31ud3at9yJBRPOM45\nVx5+q95OqF/yCE2LFmL19fSccwkkvivnnHOuungCLNWWTTRfdwUA2z92OplJ7xjhgJxzzr0RFUuA\nko6T9Iyk5ZLm5tkuST+M25+QNHWwfSXtIekuSZ3xcfdhCd6MUb+8irpXXqb/bQex44RZw/Iyzjnn\nKqciCVBSPXAtMAOYDMySNDmn2wygIy5nAz8tYd+5wD1m1gHcE5+XXcMDi2j419+xUaPp+dzFUO8f\nnTrnXLWr1BXgocByM3vWzLYDtwAzc/rMBG604EGgTdJbBtl3JvCruP4r4GPlDrxhyyaaf30NAL2z\nv4S9aXy5X8I559wIqNSlzL7AmsTz54HcyTHz9dl3kH33MbN1cf1FYJ9CAXR2du581ABjx7H62JMZ\nt/oZntvnABjqcSpsyOc7gjzmyvCYh1+1xQvVF3M5vuq1y4zlmZlJskLbh/pmdXZ2stdJp4ZjDC20\niuvs7Ky62noec2V4zMOv2uKF6oy5HCo1BLoWmJB4vl9sK6VPsX1fisOkxMf1ZYzZOefcLqxSCfBf\nQIekSZKagJOBO3L63AGcGu8GPQzoisObxfa9Azgtrp8G3D7cJ+Kcc27XUJEhUDPrk3QucCdQD9xg\nZk9JOidunwf8GTgeWA50A2cU2zce+gpggaQzgdXApypxPs4556pfxT4DNLM/E5Jcsm1eYt2AL5a6\nb2zfCEwvb6TOOedqgc8E45xzriZ5AnTOOVeTPAE655yrSQofve2aurq6dt2Tc845B0Bra6uGsp9f\nATrnnKtJngCdc87VpF16CNQ555wrxK8AnXPO1SRPgAUMVsB3pEi6QdJ6SUsSbQULA0u6KJ7DM5I+\nPEIxT5B0n6Slkp6S9OW0xy1plKSHJT0eY/5m2mOOMdRLekzSH6sh3hjHKklPSlos6ZFqiFtSm6SF\nkp6WtEzStDTHLOnA+P5ml02Szkt5zOfH/3tLJN0c/0+WN14z8yVnIUy5tgJ4K9AEPA5MHum4Ymwf\nAKYCSxJt3wPmxvW5wJVxfXKMvRmYFM+pfgRifgswNa6PA/4TY0tt3ICAlrjeCDwEHJbmmGMcXwHm\nA3+shp+NGMsqYK+ctlTHTag/elZcbwLa0h5zIvZ6Qvm49rTGTCiDtxIYHZ8vAE4vd7x+BZhfKQV8\nR4SZ3Q/8N6e5UGHgmcAtZtZrZisJ86weWpFAE8xsnZk9Gtc3A8sIP+CpjduCLfFpY1yMFMcsaT/g\nBOD6RHNq4x1EauOW1Er4Q/TnAGa23cxeSXPMOaYDK8xsNemOuQEYLakBGAO8UO54PQHmV6g4b1oV\nKgycuvOQtD9wMOGKKtVxx+HExYQyW3eZWdpj/gFwIZBJtKU53iwD7pb0b0lnx7Y0xz0J2AD8Ig43\nXy9pLOmOOelk4Oa4nsqYzWwt8H3gOWAdoTrQIsocryfAXYyF8YBU3torqQW4DTjPzDYlt6UxbjPr\nN7MphBqUh0o6KGd7amKW9BFgvZn9u1CfNMWb44j4Ps8AvijpA8mNKYy7gfAxxE/N7GBgK2E4bkAK\nYwZAoaTcicCtudvSFHP8bG8m4Y+N8cBYSbOTfcoRryfA/Eop4JsmhQoDp+Y8JDUSkt9NZvbb2Jz6\nuAHi8NZ9wHGkN+b3AydKWkUYsj9K0m9Ib7wD4l/7mNl64HeEoas0x/088HwcEQBYSEiIaY45awbw\nqJm9FJ+nNeajgZVmtsHMdgC/BQ4vd7yeAPMrpYBvmhQqDHwHcLKkZkmTgA7g4UoHJ0mEz0uWmdnV\niU2pjVvS3pLa4vpo4Bjg6bTGbGYXmdl+ZrY/4ef1XjObndZ4sySNlTQuuw4cCywhxXGb2YvAGkkH\nxqbpwFJSHHPCLF4d/oT0xvwccJikMfH3x3TCvQPljbdSd/VU20Iozvsfwt1El4x0PIm4biaMie8g\n/CV6JrAncA/QCdwN7JHof0k8h2eAGSMU8xGEoYongMVxOT7NcQPvAR6LMS8BLo3tqY05EccHefUu\n0FTHS7jT+vG4PJX9v1YFcU8BHok/H78Hdq+CmMcCG4HWRFtqYwa+Sfijcwnwa8IdnmWN12eCcc45\nV5N8CNQ551xN8gTonHOuJnkCdM45V5M8ATrnnKtJngCdc87VJE+Azg2RpIslXT94z7K+5rckvSzp\nxRL7m6S3DXdczlUj/xqEc1VC0kTCd5zaLcyaUso+BnSY2fJB+u1PmH2/0cz63mCoI/46zpXCrwCd\nqx4TgY2lJj/nXHGeAJ0bhKSvSVoraXMstjk9tl8W59tE0o8lbUksfZIui9vGS7pN0gZJKyV9qchr\ntUq6MfZdLenrkuokHQ3cBYyPx/9lgf0vkLRO0guSPpuz7YRYvWCTpDXZ+KL74+Mr8fjTJB0g6V5J\nG+Ow603Z6eEGeV/qJM2VtCLuu0DSHoVep6R/BOeGgSdA54qI8z2eCxxiZuOADxMKuL6GmZ1rZi1m\n1kKY+u1/wO2S6oA/EKb62pcwp+F5Klyx+kdAK2GKsCOBU4EzzOxuwkTGL8TXOT1PrMcBXyXMW9pB\nmFA4aWs8XhuhbuDnJWXrqWUrMLTF4/+TUBT4u4TZ+N9JmGz4shLelzmEOm1Hxn3/B1xb5HWcGxGe\nAJ0rrp8wB+FkSY1mtsrMVhTqLGlvwtyQc8zsMeAQYG8zu9xC4dRngesIE1bn7lsf2y8ys81mtgq4\nCjilxFg/BfzCzJaY2VZissoys7+Z2ZNmljGzJwjzyh5Z6GBmttzM7rJQZHQDcHWif7H35RzCnJ7P\nm1lvjOMTCoVNnUsNT4DOFRFvHjmP8Et8vaRbJI3P11eh5NNCYL6Z3RKb2wnDlq9kF+BiXi3kmbQX\nofL86kTbakovRDqe1xYFTR4HSe+TdF8cXu0iJKq9Ch1M0j7xfNdK2gT8Jtt/kPelHfhd4nyXERJm\nvnN2bsR4AnRuEGY238yOIPxiN+DKAl1/BGwCvp5oW0Ooa9aWWMaZ2fF59n+ZUOWjPdE2kdLrsK3j\ntTXRJuZsn08oGzPBzFqBeYRhTshfWPQ7sf3dZrYbMDvRv9j7soYwG3/ynEdZqPvnt5271PAE6FwR\nkg6UdJSkZqAH2AZk8vT7HGF48DNmltz+MLA53jAyWlK9pIMkHZJ7DDPrBxYA35Y0TlI78BXClVcp\nFgCnS5osaQzwjZzt44D/mlmPpEOBTye2bYjn9dac/luALkn7AhckzrfY+zIvnkN77Lu3pJlFXse5\nEeEJ0LnimoErCFdnLwJvAi7K028W4Zf6C4k7QS+OSe0jhPpxK+Nxrifc6JLPHMLNKs8C/yBctd1Q\nSqBm9hfgB8C9wPL4mPQF4HJJm4FLCQkzu2838G3ggTh0eRihHttUoAv4E6Eqd1ax9+UawpXmovha\nDwLvK/I6zo0I/yK8c865muRXgM4552qSJ0DnnHM1yROgc865muQJ0DnnXE3yBOicc64meQJ0zjlX\nkzwBOuecq0meAJ1zztUkT4DOOedq0v8Bd449yevb7nQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f33a957deb8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='size of dataset', ylabel='running time')\n",
    "bar = ax.errorbar(times_mm0[:, 0], times_mm0[:, 1], times_mm0[:, 2], label=\"$X^T X$ (PCA)\", linewidth=2)\n",
    "ax.errorbar(times_mm1[:, 0], times_mm1[:, 1], times_mm1[:, 2], label=\"$X X^T$ (PCA_high_dim)\", linewidth=2)\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first benchmark the time taken to compute $\\boldsymbol X^T\\boldsymbol X$ and $\\boldsymbol X\\boldsymbol X^T$. Jupyter's magic command `%time` is quite handy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we benchmark `PCA`, `PCA_high_dim`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "invalid syntax (<ipython-input-172-6d03ca2c8880>, line 7)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-172-6d03ca2c8880>\"\u001b[0;36m, line \u001b[0;32m7\u001b[0m\n\u001b[0;31m    mu, sigma = time(???)\u001b[0m\n\u001b[0m                     ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m invalid syntax\n"
     ]
    }
   ],
   "source": [
    "times0 = []\n",
    "times1 = []\n",
    "\n",
    "for datasetsize in np.arange(4, 784, step=100):\n",
    "    XX = Xbar[:datasetsize]\n",
    "    npc = 2\n",
    "    mu, sigma = time(???)\n",
    "    times0.append((datasetsize, mu, sigma))\n",
    "    \n",
    "    mu, sigma = time(???)\n",
    "    times1.append((datasetsize, mu, sigma))\n",
    "    \n",
    "times0 = np.asarray(times0)\n",
    "times1 = np.asarray(times1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, use the `time` magic command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%time Xbar.T @ Xbar\n",
    "%time Xbar @ Xbar.T\n",
    "pass # Put this here, so that our output does not show the result of computing `Xbar @ Xbar.T`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare the running time for __PCA__ and __PCA\\_high\\_dim__ directly. Spend some time and think about what this plot means. We mentioned in lectures that `PCA_high_dim` are advantageous when\n",
    "we have dataset size $N$ < data dimension $D$. Although our plot for the two running times does not intersect exactly at $N = D$, it does show the trend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.set(xlabel='number of datapoints', ylabel='run time')\n",
    "ax.errorbar(times0[:, 0], times0[:, 1], times0[:, 2], label=\"PCA\", linewidth=2)\n",
    "ax.errorbar(times1[:, 0], times1[:, 1], times1[:, 2], label=\"PCA_high_dim\", linewidth=2)\n",
    "ax.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, with the magic command `time`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%time PCA(Xbar, 2)\n",
    "%time PCA_high_dim(Xbar, 2)\n",
    "pass"
   ]
  }
 ],
 "metadata": {
  "coursera": {
   "course_slug": "mathematics-machine-learning-pca",
   "graded_item_id": "CXC11",
   "launcher_item_id": "ub5A7"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}