PCA-week4-FINAL.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Principal Component Analysis (PCA)\n",
    "We will implement the PCA algorithm. We will first implement PCA, then apply it (once again) to the MNIST digit dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objective\n",
    "1. Write code that implements PCA.\n",
    "2. Write code that implements PCA for high-dimensional datasets"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first import the packages we need for this week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import numpy as np\n",
    "import scipy\n",
    "import scipy.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from ipywidgets import interact\n",
    "\n",
    "plt.style.use('fivethirtyeight')\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from load_data import load_mnist\n",
    "\n",
    "MNIST = load_mnist('./')\n",
    "images, labels = MNIST['data'], MNIST['target']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's plot a digit from the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD1CAYAAAB9TzjVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAANd0lEQVR4nO3df4wU9RnH8c8jLUal7KFGvIimP4JJTaOHKCGxqdSmjUUTaEyrxHA0MdHEkqAhpGpPy38aozRqKjFV0rNSsQk1YNQq4YjGHzEoUkVpCzWEXrlwot55xEQrPP3jBnIC+51lZ/bHPbxfyWV352F2Hlc+zNx8Z+dr7i4AMZ3U6gYANA4BBwIj4EBgBBwIjIADgX2tUW88PDzM6XmgiSqVih25rNAe3MyuNLN/mtlOM7utyHsBKF/dATezCZJ+L+mnki6QtMDMLiirMQDFFTlEnyVpp7t/IElmtkbSPEnvH/kHOzo6CmwGQDVDQ0PJepFD9HMk/WfM6/5sGYA2USTgR/1CL4kTa0AbKRLwfknnjnk9TdKeYu0AKFORgG+WNN3MvmVmEyVdJ2l9OW0BKEPdJ9nc/UszWyzpBUkTJK1y9/dK6wxAYYUudHH35yQ9V1IvAErGpapAYAQcCIyAA4ERcCAwAg4ERsCBwAg4EBgBBwIj4EBgBBwIjIADgRFwIDACDgRGwIHACDgQGAEHAiPgQGAEHAiMgAOBEXAgMAIOBNaw6YMR08yZM5P1xYsXV611d3cn13388ceT9YceeihZ37JlS7J+ImIPDgRGwIHACDgQGAEHAiPgQGAEHAiMgAOBmbs35I2Hh4cPv3FHR0dDtoHydXV1Jet9fX3J+uTJk8ts5yuGh4eT9TPOOKNh225XQ0NDh59XKhU7sl7oQhcz2yVpRNIBSV+6+yVF3g9Aucq4ku2H7r6vhPcBUDJ+BwcCKxpwl/Simb1lZjeW0RCA8hQ9RL/M3feY2VmSNpjZP9z95TIaA1BcoT24u+/JHgclPS1pVhlNAShH3QE3s9PM7BuHnkv6iaRtZTUGoLgih+hTJT1tZofe58/u/rdSukLDzJqVPshau3Ztsl6pVJL11HUVIyMjyXW/+OKLZD1vnHv27NlVa3nfFc/b9nhVd8Dd/QNJF5XYC4CSMUwGBEbAgcAIOBAYAQcCI+BAYNw2eRw69dRTq9Yuvvji5LpPPPFEst7Z2VlXT7XYsWNHsn7vvfcm62vWrEnWX3311aq1np6e5Lp33313sj5esQcHAiPgQGAEHAiMgAOBEXAgMAIOBEbAgcAYBx+HHnnkkaq1BQsWNLGT45M3Rj9p0qRk/aWXXkrW58yZU7V24YUXJteNij04EBgBBwIj4EBgBBwIjIADgRFwIDACDgTGOHgbmjlzZrJ+1VVXVa1lt7GuW95Y8zPPPJOs33fffVVre/bsSa779ttvJ+uffPJJsn7FFVdUrRX9XMYr9uBAYAQcCIyAA4ERcCAwAg4ERsCBwAg4EJilpnstYnh4+PAbd3R0NGQb41VXV1ey3tfXl6xPnjy57m0///zzyXre98kvv/zyZD31vetHH300ue6HH36YrOc5cOBA1dpnn32WXDfvvytv+uFWGRoaOvy8UqkcNdifuwc3s1VmNmhm28YsO93MNpjZjuxxSmkdAyhNLYfof5R05RHLbpO00d2nS9qYvQbQZnID7u4vS/r4iMXzJPVmz3slzS+5LwAlqPck21R3H5Ck7PGs8loCUBbOogOB1RvwvWbWKUnZ42B5LQEoS70BXy9pUfZ8kaR15bQDoEy53wc3syclzZF0ppn1S/qtpHsk/cXMbpC0W9LPG9nkeHP++ecn68uWLUvWK5VKsr5v376qtYGBgeS6vb29yfr+/fuT9WeffbZQvVVOOeWUZH3p0qXJ+vXXX19mO02TG3B3r3blw49K7gVAyTjJBgRGwIHACDgQGAEHAiPgQGDcNrkOJ598crKeunWwJM2dOzdZHxkZSda7u7ur1t58883kunnDRSeq8847r9UtNAR7cCAwAg4ERsCBwAg4EBgBBwIj4EBgBBwIjHHwOsyYMSNZzxvnzjNv3rxkPW+KX+AQ9uBAYAQcCIyAA4ERcCAwAg4ERsCBwAg4EBjj4HVYsWJFsm521CyuX5E3js04d31OOqn6/urgwYNN7KR9sAcHAiPgQGAEHAiMgAOBEXAgMAIOBEbAgcAYB6/i6quvrlrr6upKruvuyfr69evr6glpqbHuvP8nW7duLbudtpC7BzezVWY2aGbbxixbbmb/NbOt2U+xOxwAaIhaDtH/KOnKYyz/nbt3ZT/PldsWgDLkBtzdX5b0cRN6AVCyIifZFpvZO9kh/JTSOgJQmnoDvlLSdyR1SRqQdH9pHQEoTV0Bd/e97n7A3Q9K+oOkWeW2BaAMdQXczDrHvPyZpG3V/iyA1skdBzezJyXNkXSmmfVL+q2kOWbWJckl7ZJ0UwN7bInUPNoTJ05Mrjs4OJisP/XUU3X1FF3evOvLly+v+737+vqS9dtvv73u925nuQF39wXHWPxYA3oBUDIuVQUCI+BAYAQcCIyAA4ERcCAwvi7aAJ9//nmyPjAw0KRO2kveMFhPT0+yvmzZsmS9v7+/au3++9MXW+7fvz9ZH6/YgwOBEXAgMAIOBEbAgcAIOBAYAQcCI+BAYIyDN8CJfFvk1C2l88axr7322mR93bp1yfo111yTrJ+I2IMDgRFwIDACDgRGwIHACDgQGAEHAiPgQGCMg1dhZnXVJGn+/PnJ+pIlS+rqqR3ceuutyfqdd95ZtVapVJLrrl69Olnv7u5O1nE09uBAYAQcCIyAA4ERcCAwAg4ERsCBwAg4EBjj4FW4e101STr77LOT9QcffDBZX7VqVbL+0UcfVa3Nnj07ue7ChQuT9YsuuihZnzZtWrK+e/fuqrUXXnghue7DDz+crOP45e7BzexcM9tkZtvN7D0zW5ItP93MNpjZjuxxSuPbBXA8ajlE/1LSUnf/rqTZkn5lZhdIuk3SRnefLmlj9hpAG8kNuLsPuPuW7PmIpO2SzpE0T1Jv9sd6JaWvzwTQdMd1ks3MvilphqQ3JE119wFp9B8BSWeV3RyAYmoOuJlNkrRW0i3u/mnjWgJQlpoCbmZf12i4V7v7X7PFe82sM6t3ShpsTIsA6pU7TGaj3418TNJ2d18xprRe0iJJ92SP6XvankAmTJiQrN98883Jet7tfz/9tPoB1PTp05PrFvXaa68l65s2bapau+uuu8puBzlqGQe/TNJCSe+a2dZs2R0aDfZfzOwGSbsl/bwxLQKoV27A3f0VSdXucPCjctsBUCYuVQUCI+BAYAQcCIyAA4ERcCAwvi5axeuvv161tnnz5uS6l156aaFt533ddOrUqXW/d+qrppK0Zs2aZH083/L5RMQeHAiMgAOBEXAgMAIOBEbAgcAIOBAYAQcCs7xbANdreHj48Bt3dHQ0ZBut0tnZmazfdNNNyXpPT0+ynjc9cer/2QMPPJBcd+XKlcn6zp07k3W0l6GhocPPK5XKUX9x2IMDgRFwIDACDgRGwIHACDgQGAEHAiPgQGCMgwPjGOPgwAmMgAOBEXAgMAIOBEbAgcAIOBAYAQcCyw24mZ1rZpvMbLuZvWdmS7Lly83sv2a2NfuZ2/h2ARyPWiY++FLSUnffYmbfkPSWmW3Iar9z9/sa1x6AImqZH3xA0kD2fMTMtks6p9GNASjuuH4HN7NvSpoh6Y1s0WIze8fMVpnZlJJ7A1BQzQE3s0mS1kq6xd0/lbRS0nckdWl0D39/QzoEULeaAm5mX9douFe7+18lyd33uvsBdz8o6Q+SZjWuTQD1qOUsukl6TNJ2d18xZvnYW4v+TNK28tsDUEQtZ9Evk7RQ0rtmtjVbdoekBWbWJckl7ZKUvlcwgKar5Sz6K5KOdaPu58pvB0CZuJINCIyAA4ERcCAwAg4ERsCBwAg4EBgBBwIj4EBgBBwIjIADgRFwIDACDgRGwIHAavm6aGFjZ0AE0DzswYHACDgQmLl7q3sA0CDswYHAmhZwM7vSzP5pZjvN7LZmbbcWZrbLzN7NpmB6s8W9rDKzQTPbNmbZ6Wa2wcx2ZI8tuQd9ld5aPoVVYnqtln9urZ76qymH6GY2QdK/JP1YUr+kzZIWuPv7Dd94Dcxsl6RL3H1fG/TyA0n7JT3u7t/Llt0r6WN3vyf7x3GKu/+6TXpbLml/K6ewyu7w2zl2ei1J8yX9Ui3+3BK9/UJN+NyatQefJWmnu3/g7l9IWiNpXpO2Pa64+8uSPj5i8TxJvdnzXo3+BWm6Kr21nLsPuPuW7PmIpEPTa7X8c0v01hTNCvg5kv4z5nW/2mt+M5f0opm9ZWY3trqZY5iazRF3aK64s1rcz5HaZgqrI6bXaqvPrRVTfzUr4Me67XI7nb6/zN0vlvRTSb/KDkVRm7aZwuoY02u1jVZN/dWsgPdLOnfM62mS9jRp27ncfU/2OCjpabXfNEx7D80kkz0Otrifw9plCqtjTa+lNvncWjn1V7MCvlnSdDP7lplNlHSdpPVN2naSmZ2WnfyQmZ0m6Sdqv2mY1ktalD1fJGldC3v5inaYwqra9Fpqg8+t5VN/uXtTfiTN1eiZ9H9L+k2ztltDX9+W9Pfs571W9ybpSY0esv1Po0c+N0g6Q9JGSTuyx9PbqLc/SXpX0jsaDVRnC/r6vkZ/5XtH0tbsZ247fG6J3pryuXElGxAYV7IBgRFwIDACDgRGwIHACDgQGAEHAiPgQGAEHAjs/3400ze7++IPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,4))\n",
    "plt.imshow(images[0].reshape(28,28), cmap='gray');\n",
    "plt.grid(False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 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 $S$.\n",
    "   Sort by the largest eigenvalues and the corresponding eigenvectors (`eig`).\n",
    "3. Compute the orthogonal projection matrix and use that to project the data onto the subspace spanned by the eigenvectors."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data normalization `normalize`\n",
    "\n",
    "We will first implement the data normalization mentioned above.\n",
    "\n",
    "Before we implement the main steps of PCA, we will need to do some data preprocessing.\n",
    "\n",
    "To preprocess the dataset for PCA, we will make sure that the dataset has zero mean. Given a dataset $\\mathbf{X}$, we will subtract the mean vector from each row of the dataset to obtain a zero-mean dataset $\\overline{\\mathbf{X}}$. In the first part of this notebook, you will implement `normalize` to do that.\n",
    "\n",
    "To work with images, it's also a common practice to convert the pixels from unsigned interger 8 (uint8) encoding to a floating point number representation between 0-1. We will do this conversion for you for the MNIST dataset so that you don't have to worry about it.\n",
    "\n",
    "Data normalization is a common practice. More details can be found in\n",
    "[Data Normalization or Feature Scaling](https://en.wikipedia.org/wiki/Feature_scaling)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-10a35b1c4186ecb4",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "def normalize(X):\n",
    "    \"\"\"Normalize the given dataset X to have zero mean.\n",
    "    Args:\n",
    "        X: ndarray, dataset of shape (N,D)\n",
    "    \n",
    "    Returns:\n",
    "        (Xbar, mean): tuple of ndarray, Xbar is the normalized dataset\n",
    "        with mean 0; mean is the sample mean of the dataset.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    ### Uncomment and modify the code below\n",
    "    N, D = X.shape\n",
    "    mu = np.zeros((D,)) # <-- EDIT THIS, compute the mean of X\n",
    "    mu = np.mean(X, axis=0)\n",
    "    Xbar = (X-mu)            # <-- EDIT THIS, compute the normalized data Xbar\n",
    "    return Xbar, mu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-8ed3ab2f7f38aab4",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "\"\"\"Test data normalization\"\"\"\n",
    "from numpy.testing import assert_allclose\n",
    "\n",
    "\n",
    "X0 = np.array([[0, 0.0], \n",
    "               [1.0, 1.0], \n",
    "               [2.0, 2.0]])\n",
    "X0_normalize, X0_mean = normalize(X0)\n",
    "# Test that normalized data has zero mean\n",
    "assert_allclose(np.mean(X0_normalize, 0), np.zeros((2,)))\n",
    "assert_allclose(X0_mean, np.array([1.0, 1.0]))\n",
    "assert_allclose(normalize(X0_normalize)[0], X0_normalize)\n",
    "\n",
    "\n",
    "X0 = np.array([[0, 0.0], \n",
    "               [1.0, 0.0], \n",
    "               [2.0, 0.0]])\n",
    "X0_normalize, X0_mean = normalize(X0)\n",
    "# Test that normalized data has zero mean and unit variance\n",
    "assert_allclose(np.mean(X0_normalize, 0), np.zeros((2,)))\n",
    "assert_allclose(X0_mean, np.array([1.0, 0.0]))\n",
    "assert_allclose(normalize(X0_normalize)[0], X0_normalize)\n",
    "\n",
    "# Some hidden tests below\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute eigenvalues and eigenvectors `eig`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-6d74b9253f3a0461",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\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",
    "    # YOUR CODE HERE\n",
    "    # Uncomment and modify the code below\n",
    "    # Compute the eigenvalues and eigenvectors\n",
    "    # You can use library routines in `np.linalg.*` \n",
    "    # https://numpy.org/doc/stable/reference/routines.linalg.html\n",
    "    # for this\n",
    "    eigvals, eigvecs = np.linalg.eig(S)\n",
    "    # The eigenvalues and eigenvectors need to be\n",
    "    # sorted in descending order according to the eigenvalues\n",
    "    # We will use `np.argsort` (https://docs.scipy.org/doc/numpy/reference/generated/numpy.argsort.html)\n",
    "    # to find a permutation of the indices\n",
    "    # of eigvals that will sort eigvals in ascending order and\n",
    "    # then find the descending order via [::-1], which reverse\n",
    "    # the indices\n",
    "    sort_indices = np.argsort(eigvals)[::-1]\n",
    "    # Notice that we are sorting the columns (not rows) of\n",
    "    # eigvecs since the columns represent the eigenvectors.\n",
    "    return eigvals[sort_indices], eigvecs[:, sort_indices]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some test cases for implementing `eig`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-a8db750754a119a2",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "def _flip_eigenvectors(B):\n",
    "    \"\"\"Flip the eigenvectors.    \n",
    "    \"\"\"\n",
    "    signs = np.sign(B[np.argmax(np.abs(B), axis=0), range(B.shape[1])])\n",
    "    return B * signs\n",
    "\n",
    "def _normalize_eigenvectors(B):\n",
    "    # Normalize eigenvectors to have unit length\n",
    "    # Also flip the direction of the eigenvector based on\n",
    "    # the first element\n",
    "    B_normalized = B / np.linalg.norm(B, axis=0)\n",
    "    for i in range(B.shape[1]):\n",
    "        if (B_normalized[0, i] < 0):\n",
    "            B_normalized[:, i] = -B_normalized[:, i]\n",
    "    return B_normalized\n",
    "\n",
    "\n",
    "A = np.array([[3, 2], [2, 3]])\n",
    "expected_eigenvalues = np.array([5., 1.])\n",
    "expected_eigenvectors = np.array(\n",
    "    [[ 0.70710678, -0.70710678],\n",
    "     [ 0.70710678,  0.70710678]]\n",
    ")\n",
    "actual_eigenvalues, actual_eigenvectors = eig(A)\n",
    "# Check that the eigenvalues match\n",
    "assert_allclose(actual_eigenvalues, expected_eigenvalues)\n",
    "# Check that the eigenvectors match\n",
    "assert_allclose(\n",
    "    _normalize_eigenvectors(actual_eigenvectors),\n",
    "    _normalize_eigenvectors(expected_eigenvectors),\n",
    ")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute projection matrix\n",
    "\n",
    "Next given a orthonormal basis spanned by the eigenvectors,\n",
    "we will compute the projection matrix. This should be the same\n",
    "as what you have done last week."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-3498a1f49501ed3d",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\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",
    "    # YOUR CODE HERE\n",
    "    # Uncomment and modify the code below\n",
    "    P = np.eye(B.shape[0])\n",
    "    P = B @ np.linalg.inv(B.T @ B) @ B.T\n",
    "    return P  # <-- EDIT THIS to compute the projection matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-65a86398cff0c3ac",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "B = np.array([[1, 0],\n",
    "              [1, 1],\n",
    "              [1, 2]])\n",
    "\n",
    "assert_allclose(\n",
    "    projection_matrix(B), \n",
    "    np.array([[5,  2, -1],\n",
    "              [2,  2,  2],\n",
    "              [-1, 2,  5]]) / 6\n",
    ")\n",
    "\n",
    "# Some hidden tests below"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compute principal component analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-98795bac9f7e1f82",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\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",
    "        the reconstructed data, the sample mean of the X, principal values\n",
    "        and principal components\n",
    "    \"\"\"\n",
    "\n",
    "    # YOUR CODE HERE\n",
    "    # your solution should take advantage of the functions you have implemented above.\n",
    "    ### Uncomment and modify the code below\n",
    "    # first perform normalization on the digits so that they have zero mean and unit variance\n",
    "    X_normalized, mean = normalize(X) # EDIT THIS\n",
    "    # Then compute the data covariance matrix S\n",
    "    S = np.cov(X_normalized, rowvar = False, bias = True) # EDIT THIS\n",
    "\n",
    "    # Next find eigenvalues and corresponding eigenvectors for S\n",
    "    eig_vals, eig_vecs = eig(S)\n",
    "    # Take the top `num_components` of eig_vals and eig_vecs,\n",
    "    # This will be the corresponding principal values and components\n",
    "    principal_vals, principal_components = eig_vals[0:num_components], eig_vecs[:,0:num_components]\n",
    "\n",
    "    # reconstruct the data from the using the basis spanned by the principal components\n",
    "    # Notice that we have subtracted the mean from X so make sure that you add it back\n",
    "    # to the reconstructed data\n",
    "    P = projection_matrix(principal_components)\n",
    "    reconst = np.zeros_like(X_normalized)\n",
    "    reconst = (P @ X_normalized.T).T + mean\n",
    "    return reconst, mean, principal_vals, principal_components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def draw_vector(v0, v1, ax=None, label=None):\n",
    "    \"\"\"Draw a vector from v0 to v1.\"\"\"\n",
    "    ax = ax or plt.gca()\n",
    "    arrowprops=dict(arrowstyle='->',\n",
    "                    linewidth=2,\n",
    "                    shrinkA=0, shrinkB=0, \n",
    "                    color='k')\n",
    "    ax.annotate('', v1, v0, arrowprops=arrowprops, label=label)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some test cases that check the implementation of PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-90d1f11031fdf7d8",
     "locked": true,
     "points": 0,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "D = 2\n",
    "N = 10\n",
    "# Generate a dataset X from a 2D Gaussian distribution\n",
    "mvn = scipy.stats.multivariate_normal(\n",
    "    mean=np.ones(D, dtype=np.float64), \n",
    "    cov=np.array([[1, 0.8], [0.8, 1]], dtype=np.float64)\n",
    ")\n",
    "\n",
    "X = mvn.rvs((N,), random_state=np.random.RandomState(0))\n",
    "reconst, m, pv, pc = PCA(X, 1)\n",
    "# Check the shape returned by the PCA implementation matches the specification.\n",
    "assert reconst.shape == X.shape\n",
    "assert m.shape == (D, )\n",
    "assert pv.shape == (1, )\n",
    "assert pc.shape == (D, 1)\n",
    "\n",
    "# Check that PCA with num_components == D gives identical reconstruction\n",
    "reconst, m, pv, pc = PCA(X, D)\n",
    "assert reconst.shape == X.shape\n",
    "assert m.shape == (D, )\n",
    "assert pv.shape == (2, )\n",
    "assert pc.shape == (D, 2)\n",
    "assert_allclose(reconst, X)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-0c7ba302f4b562c7",
     "locked": true,
     "points": 4.5,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize PCA\n",
    "We will first visualize what PCA does on a 2D toy dataset. You can use the visualization\n",
    "below to get better intuition about what PCA does and use it to debug your code above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAakAAAFzCAYAAACEifYSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de3hU1b038O+aSTKZkCEJITGESwiQWqwConJJoFRFkWLV0pfWS7UtR3vwlB6LgkiNr1pRQAFbr1SUHu89UlFbFFFfLBawKiAqghggBMMQJ+Q6SSaZycx6/5gk5DL37Jm9Z8/38zw+D0wmM2uDzDdrrd/+LSGlBBERkRYZ1B4AERGRPwwpIiLSLIYUERFpFkOKiIg0iyFFRESalaT2AELV0NDAMkQiIh3LyMgQvR/jTIqIiDSLIUVERJrFkOqmrKxM7SFEld6vD9D/Ner9+gD9X6Perw9Q9hoZUkREpFkMKSIi0qy4qe7zR0qJpqYmeDyefr9WamoqGhoaFBiVNqlxfQaDAenp6RCiT9EOEVFQcR9STU1NMJlMSElJ6fdrmUwmpKamKjAqbVLj+pxOJ5qammCxWGL6vkSkD3G/3OfxeBQJKIqOlJQURWa5RJSY4j6kiIhIvxhSCluxYgUeffRRv1/fvHkzvvrqqxiOiIgofjGkYuzNN9/EoUOH1B4GEVFcSLiQqrC7cNP2Wly+pRo3ba9Fhd3V79dcvXo1zj//fFx55ZU4fPgwAODZZ5/FhRdeiJKSElx//fVoaWnBRx99hC1btuCuu+7CtGnTUF5e7vN5RETklVAhVWF34aqtNdh41IEdVU5sPOrAVVtr+hVU+/btw6ZNm/DBBx/g+eefx969ewEAP/rRj/D+++9j586dOPPMM/H8889j8uTJmD17Nu677z7s2LEDhYWFPp9HRARE54fqeBP3JejhWL7XjnK7u8dj5XY3lu+1Y/2MQRG95q5duzBnzhykpaUBAGbPng0AOHDgAO6//340NDSgqakJF198sc/vD/V5RJRYOn+o7v6ZtbvaiddnZaPAkqziyGIroWZSJ1vcPh+v8vN4qHzdqPpf//VfePDBB7Fr1y4sXboUra2tPr831OcRUWIJ9EN1IkmokBqSZvT5eJ6fx0NRXFyMzZs3w+FwwG634+233wbgvck4Ly8PLpcLGzdu7Hp+eno67PbT/5P5ex4RJbZo/VAdbxJqua90ogW7q509fjoptBhROjHybggTJkzA3LlzMX36dAwfPhxTp04FANx55524+OKLMXz4cJx11lloamoCAPzkJz/BLbfcgj//+c947rnn/D6PiBJbNH6ojkdCyvg48NbfybwNDQ3IyMgI+XUq7C4s32tHVYsbeWnegOpc321tbdV1WyS1ri/cv6P+KCsrQ1FRUUzeSw16vz5A/9cY6vX52pMqtBjjYk8q0r9DXyfzJtRMCgAKLMkRF0kQEcVKgSUZr8/K9vtDdaJIuJAiIooX/KGaIUVEFLeEzYqUTRsg6k5BZg2Gc+58yNx8tYelKIYUEVEcEjYrzA8thsFm7XrMeOQAHEtW6yqoEqoEnYhIL1I2begRUABg6JhZ6QlDiogoDom6U74fr6+J8UiiiyGlA2vWrFHsterr6/H000+H/X3BjighImXJrMG+H8/MjvFIooshpSAppSqn0K5du9bn45GMp6GhAc8884wSwyKiKHLOnQ9Pr70nT24+nHPnqzSi6Ei4kBI2K0zrliN1xe9gWrccoteabrgqKiowadIk3Hbbbfj+97+PBx98EBdeeCGKi4vxwAMPdD3v5ZdfRnFxMUpKSvDrX/8aAHD8+HFcccUVKC4uxhVXXIFvvvkGAHDzzTfj9ttvx6WXXorx48fjjTfeAABUVVVh9uzZmDZtGqZOnYpdu3bhnnvugcPhwLRp03DTTTf1GU9lZSWGDh3aNY433ngDN998MwDAZrPhuuuuQ0lJCUpKSvDRRx/h3nvvRXl5OaZNm4a77roLAPDII4/4vCZfR5QQUWzI3Hw4lqyGa+pMtI89F66pM3VXNAEkWHVftKphysrK8Pjjj2POnDl44403sG3bNkgpcc0112Dnzp0YNGgQ1qxZg61btyI7Oxt1dXUAgCVLluDqq6/Gtddei+effx5Lly7FSy+9BAD49ttv8fbbb+Prr7/GNddcgyuvvBIbN27ExRdfjMWLF8PtdqOlpQXFxcVYv349duzYAcAbmp3jCbYMuHTpUpSUlODFF1+E2+1GU1MT7r77bhw8eLDr9bZt24YjR470uaYBAwZ0HVHS3t6OGTNmYMKECRH/GRJR+GRuPtoWlKo9jKhKqJAKVA3Tn7/o4cOH44ILLkBpaSm2bduG6dOnAwCam5tx5MgR7N+/H1deeSWys71rxVlZWQCATz75BC+88AIA4Oqrr8bdd9/d9Zpz5syBwWDAd7/7XVRXVwMAJk6ciIULF8LlcmHOnDkYN25cwPEE88EHH2DdunUAAKPRiIyMDNTX1/d4zrZt23xeU1NTk88jSoiIlJRQy33RqoYZMGAAAO8e0K233oodO3Zgx44d+PTTT3HDDTdASunzOI8+4+j2HJPJ1PXrzv6KJSUleOutt5Cfn48FCxbg5ZdfDjgeX68b7lEg/q6p9+sSEUVDQoVUtKthLr74YrzwwgtdncytViuqq6sxY8YMvPbaa6itrQWAruW+SZMm4dVXXwUAvPLKK5gyZUrA1z9+/DhycnLwi1/8Aj//+c/x2WefAQCSkpLgcvk/sTMnJweHDh2Cx+PB5s2bux6fMWNGV5GE2+1GY2MjLBZLj6NE/F2TvyNKiIiUlFDLfc6582E8cqDHkp+S1TAXXXQRDh06hEsvvRSAd0bz1FNPYezYsbjtttu6lvDGjRuHJ598EqtWrcLChQvxyCOPYPDgwXj88ccDvv6OHTvw6KOPIikpCenp6XjyyScBAL/85S9RUlKC8ePHo7S077LlPffcg5/97GcYMmQIzj777K7AWblyJW655Ra88MILMBgMWLt2LSZNmoQpU6Zg6tSpmDlzJu677z6f1+TviBIiIiUl3FEdXb2u6msgM7N79LriUR3RwaM6lKP36wP0f416vz6AR3X0SyJUwxAR6UXChRQRaVPngaRHT5kwylqbkGcnUV8MKSJSXc9TaI3Y0+jA7mpnXJxCS9GVUNV9RKRNy/faexyTDgDldjeW77X7+Q5KFHEfUgaDAU6nU+1hkB9OpxMGQ9z/b0ZRdrLF7fPxKj+PU+KI++W+9PR0NDU1weFw9Pu1GhsbMXDgQAVGpU1qXJ/BYEB6enpM35Piz5A0o8/H8/w8Tokj7kNKCAGLxaLIa9lsNgwfPlyR19IivV8fxa/SiRbsrnb2WPIrtBhROlGZf9sUv7gOQ0SqK7Ak4/VZ2Zg3yozzMtyYN8rMogkCoIOZFBHpQ4ElGetnDEJZWQ2KikaoPRzSCM6kiIhIsxhSRESkWQwpIiLSLIYUERFpFkOKiIg0SxMhJYRIFUJ8LIT4TAjxpRDiXrXHRERE6tNKCXobgIuklE1CiGQAO4QQW6SU/1Z7YEREpB5NhJT0nrzY1PHb5I7/4uM0RiIiihrNnMwrhDAC2ANgDIDHpZRLu3+9+8m8ZWVlMR4dERFFQ/cTfH2dzKuZkOokhMgE8BqA30op93c+7u/4eCXp/VhnvV8foP9r1Pv1Afq/Rr1fH6Ds8fGaKJzoTkpZD+CfAC5TeShERKQyTYSUECKnYwYFIYQZwEwAX6k7KiIiUpsmCicADAHwbMe+lAHAK1LKzSqPiYiIVKaJkJJSfg7gXLXHQURE2qKJ5T4iIiJfGFJERKRZDCkiItIshhQREWkWQ4qIiDSLIUVERJrFkCIiIs3SxH1SRKQtFXYXlu+142SLG0PSjCidaEGBJVntYVECYkgRUQ8Vdheu2lqDcru767Hd1U68PiubQUUxx5Aioh6W77X3CCgAKLe7sXyvHetnDFJpVJELdVYobFakbNoAUXcKMmswnHPnQ+bmqzBi6o4hRUQ9nGxx+3y8ys/jWhbqrFDYrDA/tBgGm7XrMeORA3AsWc2gUhkLJ4iohyFpRp+P5/l5XMsCzQq7S9m0oUdAAYChY2ZF6mJIEVEPpRMtKLT0DKRCi3eZLN6EOisUdad8Pk/U1yg+JgoPl/uIqIcCSzJen5WN5XvtqGpxIy+Oq/tCnRXKrME+nyczsxUfE4WHIUVEfRRYkuOySKK30okW7K529ljyK7QYce+IJpjWPdFVJOGacTmMRw70WPLz5ObDOXe+GsOmbhhSRKRbvmaF96cfReGKuyDaHF3PMx45gNb5tyN5+2aI+hrIzGxW92kEQ4qIdK37rFDYrDCX9gwowFskkbx9M9oWlKoxRAqAIUVEcS/Ue6FSNm2AoVdAdWKRhDYxpIgoroXTIcNfFR/AIgmtYgk6EcVUhd2Fm7bX4vIt1bhpey0q7K5+vV6o90IBAar4TGYWSWgUZ1JEFDPR6AsYTocM59z5fav4TGa0LlrBIgmNYkgRUcxEoy+gv3uhznGfgmndn/v04nMsWe3t0ccqvrjAkCKimIlGX0Bf90JNM5zCQ++vRPKpk12Pde/FF6iKj8eUaAv3pIgoZvzNer6qb494f6rzXqh5o8yYnpeCeaPMeL3xDZi6BRQQWi++zuXIjUcd2FHlxMajDly1tabf+2YUOc6kiHRIq7OB0okWfFjVisoW2ePx6lYPNh51YHe1Ew9/R6AozNft3SEjdVetz+cFKzPX2zElesCQItIZzR9aKAQA6fNL5XY31h1Pwg/G9e8tIu3Fp6djSvSCy31EOhNOSXYsdC85/9HbNahs9gR8frWz/x9Lzrnz4elVDBFKLz49HVOiF5xJEemMlmYDvmZ1weSkBA6xUERaxeevIW08HlOiFwwpIp2J9Wwg0P6Xr1ldIIUWIxaM8N22KFzBqvh80dMxJXrBkCLSmVjOBoLtf/mb1XU3IAkYm5mEwoHJKJ1ogbNKnWXJTno5pkQvGFJEOhPL2UCwajh/s7oR6UYUpBt9jq2syvd7abVikaKLIUWkQ7GaDQTb//I3qwu30rD3jK3AYcPhrX9DQVoTTINz2DVCxxhSRBSxYPtfSs3qOmdsBQ4b1h5+HpfUfoE02XGDbVnPbhKkLwwpIopYKPtfSszqTra4UVx3AJu/WI2BnrY+X+/sJsFDC/WHIUVEEYvV/tcPGg7i95+tQAr8l6fz0EJ9YkgRUb9Ee/9L2Ky4a9sDMAYIKICHFuoVQ4qIfNJKNV3KS4/B2NYa8DmhdJOg+MSQItKh/gZMf/v/KRVwwmZF0v7dfr8uhQHtE6bAee1CFk3oFEOKSGeUaDDbn27gSry/6CiEMH65B8Ll9PkcaTDCcfsaeMZOCOk1KT6xwSyRzijRYLY//f/6+/7J9dUwr1iE5A/fg6GxzudzpBAMqATBkCLSGSUazPan/19/33/45udgqP024HPaJ0xlQCUIhhSRzijRYLZ0ogWFlp7PD7X/X3/eX9isGFhxKOBzPLn5cF67MOhrkT5wT4pIZ5RoMNuf+58ief/Te1C7IaTvAxGl0Yj2SReG3QJJK1WKFBmGFJEGKPlBqtQNtpHe/xTu+wubFeaHFsNgswZ83fZxk8PuKKH5U4opKIYUkcqi8UGq9nET4bx/ykuPBQ0oT4opoiW+/lQpkjZwT4pIZVo77j1WhM0K0x9/j6RPPwz4vCZjKo4tuD+i+6C0dEoxRYYzKSKVJeIHqbBZYV5xCwy11X6fU5U8EP8v62z838J5mNRUgPURvE+sTykm5WkipIQQwwE8ByAPgAfAU1LKP6k7KqLYSMQP0pSXHgsYUC0iGSUT70WFORcAMCLCwI7lKcUUHZoIKQDtAG6TUu4VQlgA7BFCvCulPKD2wIiiqcLuQpPTg1Qj0Nrtc1jPH6TCZoXxC/+tjgDgnaxzugIKiDywY3lKMUWHJkJKSnkSwMmOX9uFEAcBDAXAkCLd8lUwkWoALhpqworJGbr8IO2q5Gv33eoIAI6Zc3Fb0fVdv+9vYKtdREL9o4mQ6k4IMRLAuQA+UnckRNHlq2Ci1QMMSDbEZUCFUkafsmlDwEo+d9ZgHPnJIkzyFGAEZz4EQEg/N86pQQiRDmA7gPullJu6f62hoaFroGVlZbEeGpHiFnxuwp7GvstY52W4se6cvqfPatkJh8DCL02obD1dMDws1YPHvteGoebTnzFjnl8Ni5+OEm0DB6HshsVwZeaE9b7rjiehus2AHJMHC0a093g/0r6ioqKuX2dkZIjeX9fMTEoIkQzgVQAv9g6o3rpflJLKysqi9tpaoPfrA+LrGkdZa7Gn0dH38ex0FBWN8Pk9Wr2+B7fXorK157VUthrwznEX1lS+ClF3CjJrMDAoG6jo+/3uwXloX7oWI3PzQ77GCrsLi3oslxpxqDVV8zfqavXvUElKXqMmQkoIIQA8A+CglHKt2uMhigU9VZ75KqMvrjuAVf96CMnu0/tPnoHZ8Aw6o0cDWU9uPlqXrA77PijeqJsYNBFSAEoAXA/gCyHEvo7Hfi+lfEvFMRFFVTxVngXbb+pdRl/gsGHL5w/CLF09Hjc01sB11kS4zzwHor4GMjM77F58nRLx/rJEpImQklLuANBnLZJI7+Kh8iyUtk29Z4V/+fJxDOgVUJ2ajx7BidWr+h3GiXh/WSJiWyQiCiiUtk2ds8J5o8z4L1GG6U2H/b5eq1viqq01qLD7DrFQ9ec4EYofDCkiCijUZbXOWeHDX6wLuCzyb8sYRXoTdg/G6XkpmDfKrPmiCQqfJpb7iEgdodzbFOqyWueZUIYam9/3a0ZS1426SuwdxcNyKfUPQ4ooQYV6REgoVYiGg/uQ+vAyGNr6ltR3cgOYPX5pV7sj7h1RKLjcR5SgQj0iJNiymuHgPpgfvC1gQHkAXPPd32BX1lkAuHdEoeNMiigOdC7LHT1lwihrrSKl6uGUcPtbVhM2K1IfXgbh8f1a0mCAZ1AuKq5dDGPTKEzXeKk9aQ9Dikjjei7LGbGn0aHIEej9LeH2LvHdAUNbq9/ntE++CG0LSpELRHQeFBGX+4g0Llon9/anhPv0Ep//gPKYzHDOnd+vMRJxJkWkcdHqrBBpxwvvEt8dfpf4AEAajGhdtCKiThJE3TGkiDQump0Vwi3h7joPKsgMqnXRCnjGTuj3+Ii43EekcVrprCBsVqTe/98Bz4PymFLhWP4MA4oUw5kUkcZ1X5Y7WtOEUdnpMa+OMxzch9S1S2Fw+j/nyrvEt5JLfKQohhRRHOhclisrq/F71lS0GA7ug3nlIgj4P0zQnZKKtltX9msGFY0ye4p/DCki8kvYrEhdszRgQLWIZPxw2io8PuwsFET4PtEqs6f4xz0poghV2F24aXstLt9SjZu21/a7q7cWpWzaAIMr8FH27w46Bzs8g/Gjt2si/rOIVpk9xT/OpIgiEGrfu3gn6k4F/PoxUzZuHeNtGHu8yY3jTd4/j3D/LHiAIfnDmRRRBBLlJ3+ZNdjv16zJGbh4QmlXw9juwv2z4AGG5A9nUkQh6H2kRbm93efz9PaTv3PufBgPfQZDbXWPx6tMWZg+4f/6DKiu54TxZxFKp3VKTAwpoiB8Le0NSPJ9rJ/efvKXuflwLPsTUl56DMYjBwEA7tFj0XjlAkw6no4RLW5UdFvm6y7cP4uxmUlocnngdrsx+QwzVk7JUHzpNJTzs0hbGFJEQfha2mtulxiQBDR3m1DF80/+wmZFyl/WIOnrzyGlB9KShbYFpfCMnQCZm4+23z3Q4/n5ANYXen/tO8SBo43ewpJgQdD3+w34qsH3TLU/EmUfUW+4J0UUhL9N/bGZSTE7uryzknDB5ybFKwmFzQrz/b9F8oE9EO0uGNxuGOtPwbxyEQwH9wX9/u7nTV2Qk4wBSQLN7cCeU+3YeNSBq7bWBBxvrPb3EmUfUW8YUkRB+NvULxzovcH2H7NzsH7GoKgG1FVba7DxqAN7Go0hffCHStisSF11Kwz1NX2/BgnT0ytDep0CSzJKJ1rwrcOD5vae91QFC4JYVfaxgjA+MaSIglC7d160ZgCGg/uQVvofMJ6q8vsc0dIU0mt1BqmvvSkgcBDEqrKPFYTxiSFFFESw49OjLRozgK4TdQMc+Q4AMi09pNfzFaTdBQqCX37HjAG9dsej8UNA6UQLhqX1LHgZlibidh8xUbBwgigE4R5poSSlZwDCZkXqA/8NQ7CAgkDbjXeE9Jr+ghQIHDgVdhcW7mzoUYBiNnjwWEmUlk+FALq3eBK+qzRJOziTItI4JZcbDQf3Ie33v4IxQCcJNwSspkyU//dDITeM9RekI9KNAWedvmZgDo8B//N14ACNxPK9dlQ2e3o8VtnsYeGExjGkiDSu+3LjeRnuiJcbu5b4AvTiazSY8J3JazFi6uO4s2lUyK/tL0j/cVngccaymIGFE/GJy31EcaC/R3XY9uzGkCfuhKHdf0C5IHD5OYu7ukiE8+Ed6VH0sSxmYOFEfGJIEemcbc9uDH90KZJl4NB5c9C52JV1Vtfvw/3wjmTfzlc7pGGpnqgUM7D1UnxiSBHpmLBZMeSJO4MGVKU5G7cVXd/1+2DFDkq1FvI1A7suqzYqRRORzvZIXQwpogDiudebsFlhfmhxwCU+KQTax05Ey09v6erFF+jDOxqthXrPwMrK+t5YrBQ1qzQpMgwpIj/iuddbVyeJADfqOpJM8Kz4C2Rufo9efIEEurGYH/4UDazuI/IjXnu9dc6gAnWScMGAz371B8jc/LBemxVyFGucSRH5Ea8fyCmbNsBgs/r9ut1gwpxzFmOouwjre30t2PImK+Qo1hhSRH7E6wdyoCPfD6fmYtb4Zagw52J6r7ANZXmTFXIUa1zuI/JD7caykfJ35Hu5aXBXQAF9wzaU5U21+xhS4uFMisgPLZQs915+uy5LoCjI9zjnzofxyIEeS34VaWdg5jl3dAWUr7ANdXlT6xVywZYs47liMxExpIgCUPMD2dfy24epJrxZ6Ar4oSpz8+FYshopmzZA1NdAZmajbdYNQUvM43V5s7tgS5bxXLGZqBhSRBrla/mtstUQUrm3zM1H24LSrt+HUmKuh/2mYCXyLKGPPwwpIg3pvhR1qL7d53OiVV2oheXN/gq2ZBmvFZuJjCFFpBG+lqJ8Ocd9CqZ1f4aoOwWZNRjOufPDvt/Jn1CXN7W6rxNsyVIPS5qJhiFFBG186AY73RYApniq8NC7K5FcX931mPHQZ3As+5NiQRWMlvd1gi1Z6mFJM9EwpCjhRfKhG41Q87cUlZMq8N3MZOSlGfGnd9bD1C2gAMBQW42Ulx5D2+8e6Nf7hyrSfZ1Y/CAQbMlSD0uaiYYhRQkv3A/daM0k/C1F/SA/FetnDPK2O/rLIZ/PMR45GPH7hiuSfZ1Yzr6CLVlqvYSeeuLNvJTwwv3QjVZPP383D987ogmmP/4e5mW/hEHKfr2HEiLZ14nXPoikPs6kKOGF+6EbrQoxX0tR945owphHboOhtjrg97pHj+3Xe4cjkn0dVtVRpPoVUkKIG6SUzyk1GCI1hPuhG0qoRbr/0n0pynvcxu1BA0qmmOC8dmHQ11ZKJPs6rKqjSIUUUkKIs3w9DOA/ATCkKK6F+6EbLNTC3X/xFWgjHdXeAwtPVcEjJQxC+ByLx2RG66IVMavs6xTuvg6r6ihSoc6k/g3gb/AGU3cFSg1ECLEBwOUAbFLKs5V6XaJQhPOhGyzUAu2/lE609AikX37HjIU7G/oE2icnn8YAmxVvf9uIX+ytwKLRubjjO2f0eE334Dy0Ll0b84CKBKvqKFKhhtRBAEuklD3OdRZCvKngWP4HwGPgzIziQKBQ87f/csze3meG9dZxB5p7NZbw2E4C+/fgXZsdP/m4HG0eiWRDz58PPYPOiJuA6sSqOopEwOo+IURnt69LANR3e3wKAEgp5yg1ECnlBwBqlXo9IrX423850ezuM8PqHVDFdQfw6Sd3YO83lfjxR0fR5pH4z5HZuHV0DgDAY0yC69xiOJY9HFcBRRSpYCXo+4QQ86WUjVJKtxDCKIS4H8AHsRgcUTzyVUo+bIAB1S2egN9X4LBh8xersc9Wg6s+KkerR+I/Cgbh0XHDIISAJzcfB27+A9p+9wADihKGkAHuuxBCtAAwAdgMYG3Hf+fCu2+Up/hghBgJYLOvPamGhoaugZaVlSn91kSKOuEQWHc8CdVOA3JSPGhpBz6o873/YjZ44PAY8OyBxzHy4Dv44YdH0ez24BfDB2H9ucNhEALOARZ8/atlcGXmxPhKiKKrqOj0CWkZGRl9KoSC7UmNB7Ae3oKGyzse+18Asat39aH7RSmprKwsaq+tBXq/PkA711gE4AfjTv/+8i3VAJx9npdqBDZPbIfx1Q1wHv4X5nQE1HXDsvBUR0ABgBg3GSMvKNbM9UWT3q9R79cHKHuNwULqKID3AEyHt7LPAeCd3gUURBSYr32qAocNrxx7BhO3fYndtc24YtdhNLk9+NnQTGyYOALGjoDymMxwzp2PCrsLdx1KRtPhak11HvdFzYa9WmgWTMoJFlK7AYwDcBLA0wAWAXhaCPEzKeVlSg5ECPEygB8AGCyEqARwt5TyGSXfgyhS/f3g632fUIHDhh377sGQtgbsrW/B7A+PoLHdg5/kZ+DZiQVdASU77oM6Zs7pqAxMRueMLNzed7H68FazS7qWO7RTZEJZ7nsFwM1SyjohxLMAnoW32k9RUsprlH5NIiUo8cHX+z6hP77/DIa0NWBfQwtm7TqCepcbVw3JwAvnjUSSQcCTlg73+CldZ0Ut317brxNlY/nhrebptzx5V3+CVfddL6W8WkpZBwBSynIAMwAsi/rIiDRCqeaonfcJvTnyBMZ9ux9fNDowa9cR1LncuDxvIF46v6Drfij3+CloW1DaVcXX3953sWzwqmafPvYI1J+AMykp5Ys+HpMAHozaiIg0RskPPmGzwvzwMhxsbMWlO4+gxunG7DMG4n/PH4kUg/dnRtmxB9Vdf3vfxfLDW80+fewRqD88qoMoCCU/+FI2bcChU3W4ZNdhVDvbcUmOBRsvGAmT0ftP0WMwwOGjF5+/YzxC7X0Xyw/v/o41Xt+boi324YIAABk5SURBVINHdRAFoWRz1CNHjuKSnYfxbVs7Lhycjk2TC5HaEVDSYEDr7WvhGTuhz/d17mkt3X4CzcYBYfe+i2WDVzX79LFHoP4wpCgiiVTmq9QHX3l5Oeb87zs42daO72cPwOuTC2HumkEZ0Xr7Gp8B1X0c953pQlFR+Df0xvrDW80+fewRqC8MKQpbIpb5RvrBJ2xWpGzagIpj5bjir+/A2tiEkjOy8Pfzh2FAkndZqrPMPFBAKYEf3hSPGFIUNpb5hsZwcB/MDy/DN3UNuGznYXzT4sSU3Ey88uKLMO3ehvb6GsjM7K4ycyLqiyFFYWOZb3DCZkXqw8twor4BM3cexrEWJy7ITMOb5w9H2u5taFtQqvYQieICQ4rCFm9lvrHeP/Me+34rquobccnOIzja4sR5GWZsKR6FjGQj2utD7yrWe+zXZQnou+sbUU8MKQpbPB0FrvT+WbDAEzYrzA8thq3yG1yy6zDKmtswIcOMLcWjkZns/ecmM7MjHvuHqSa8WejS7d4fUW8MKQpbPJX5Krl/FizwOmdQp058g0t3HcGhpjacMzAVb08djUEpHQHl40bdcMZe2WpQfO8vkSo1Kf4wpCgi8VIppuT+WaDA25CyF6lPrUCNow2zdh3BAXsrzrKkYmvxaAw2ef+ZeTqq+EItkojF3l8iVmpSfGHHCdI1JffP/IVGTvnnSF23HHWtbZi16zC+aGzFd9NNeLd4NHJN3g969+A8OJY/E1aZeSz2/mLZ048oEgwp0jUl2+T4OxNq1b9WoN7Zjst2HcFnja0oGmDCuyVjcEaqN6A8ufloXbo27DJzn8fQp3oU3ftjpSZpHZf7SNeU3D/rXTBSXHcAb+1fA2drC2Z/eBR7GxwYPSAF75WMxpDU0zOo1iWrI7oPytfYr8uqVXQZLt4qNSnxMKRI95TaP+seGrBZsWHHGrham/HDD49gd30LRqal4N3iMRhqTgEASIMxohlUoLGXlSl7KHY8VWpSYmJIEYWhMzRM655AW2szLv/3UXxU14IR5mS8VzIaI9I6AgpA66+Xab6TRDxValJiYkgRhaGzF1/r7h244t9Hsau2GcNSk/FeyRiMTDMBADymVLQuWhn1XnxKiZdKTUpMDCnShc57fY6eMmGUtTYqs4HOG3VbrZWY+9FR/KumGfmpyXi3ZAxGDegMKDMcy5/R/AyKKF4wpCju9bzXx4g9jY6o3OuTsmkDnCcrMffjcrx/qgl5piS8WzwaRenegJJh3gdFRMExpCjuRdpVItROC51LfO17duL/fFyO96rtyDUl4Z3iMTjTkgpPWjrc46ewmzlRFDCkKO5Fcq9PqJ0WOpf42qsqcfXHx7DVZsfgFCPeKR6NswamAgDc46ewqzlRlPBmXop7kdzrE0qnhc5efO6qE7jmkwq8+W0jBiUbsbV4DM4eaAbgvVE31F58RBQ+zqQo7kVyr0+w2Zfh4D6kPrwMbkcLrtt9DG9UNSAz2Yi3i0djfIa5xxLfMXMOlm+vZYNWoihgSFHc636vz9GaJozKTg8aFIFmX8JmhbkjoH6xtwKbTjZgYJIBW6aOxsTMNACnl/hi0aC1+95ZujsZq/J4VAclDoYU6ULnvT5lZTUoKhoR9Pn+Zl/3jmhC6qrb4WltwfxPj+OVE/WwJBnw1tTRuCDLG1Ddl/iUOAokUAFH3xBMxqGtNexSTgmDIUUJyVenhXtHNGHME3cA1Sdx46fH8VJlHQYYDdg8ZRSmDBoAoGcvvgq7C/+0tvp8/VAbtAabiSl5HhZRPGJIUcLq3WnBtO4J4NsTWLDvGzz3TR3SjAb8Y8oolGSnA+g4D6qjF19nuFS3Sp+vHWqD1mAhxC7llOgYUpTwrOXHYXtuPcYd3407Pq/EhuO1MBsF3phciO8P9gZU7xt1fYVLp3AatAYLIXYpp0THkKKEZi0/DtODS1DSXIVbvjiBp47VINUg8NqkUbgwxxs07sF5fbqZ+z0AMdUQ1n5RsBBil3JKdLxPihKWtfw4zKtuxYjmKty2/wSeKD+FFIPAq5MLMTPXGwL+Diz0Fy4/yDeFVdAQ7FDGzr2zeaPMmJ6XgstyXCyaoITCkKKEZC0/jrQVizCspRpLD1jxyNFTSBYCGy8YiVm5A2FPToNr6kw4/BxYqNSJv71DaN4oc58QKrAko3SiBXlpRlS3GbB8rx0VdldkF04UZ7jcRwmp4elHUNR6CncePIm1h6uRJIBXLhiJOXkZAIDPhp+HCQFaHSl5DlOwozJ2nnTgp+/VobldIpoNdIm0iCFFCWlUdRnu+aoKD5bZYBTAy+ePxI+GeAOqIu0M5N5wU9DXiMU5TBV2F376Xi2a23s+zjJ0ShQMKUpIaw8cw6qvv4UBwAvnFeDH+ZkAgFaRhLbbH0J+YfAbgmNh+V57n4DqxDJ0SgQMKdKNCrsLdx1KRtPh6oA99NauXYtV+ytgAPDseQWYNzSr62vOsy/QTEAB/qsIAZahU2JgSJFiQj2fKVrv7e3ckAzACaBn5wbDwX0wPb0SD+/5Cn/YdwxCCDw17Xu4ZtDpfwLOrFwYb/gtfN+eqw5/VYQDkgTL0CkhMKRIEbFotBpIoM4Nz+Qeh/nB2/BI2Uncsd8KAHjq3BH42e/vgeuzf0PU10BmZsOlwUMLfd0nNSAJeGVmFosmKCEwpHQuVh201e4xF6hzg+nplXjicBVu6wiodeOH41fDs+D+29NwrPlr1MfWH72rCAe4m7FqxlAGFCUMhpSOxbKDtto95nwtixU4bLjv401Yv/tL3PLFCQDAo+OG4caR2QAA0dIUk7H1V/cqwrKyegYUJRTezKtjoZw+qxS1e8x1v7l25v5Xcf8zl+HAh7fg0PuvY+Fn3wAAHj57KG4uHNz1PTItPSZjI6LIcSalY7Gc3ajdY65zWezvz/0Nb722CnfWNqPZ6cSKr78FADz0vXz8dnRO1/OlwYi2G++IydiIKHIMKR2L5exGyQ4M/RnDdTvX4a7aZiQLgQc6AmrFWUOwaEwuAEAaDPAMykXbjXfAM3ZCzMZGRJFhSOlYrGc3sejAEIiwWfGPb04BAFzSW0h+fqYZFS1OtHskkgwC7ZMvQluAdkdEpC0MKR3zXRmWq9uN95RNG/DUsZoej+2ud+CAvQ0PnJUP+8B8ZHYc+x7P1LwfjSjWGFI6l0iVYa5T3+JL++nj3EenpWDe0CxcPyILn2V/Fzef/RtY3wEyTVV4cloGSoaYVRxtZE44BBapeD8aUawxpEg3kgefgUtzLfBIifvGDsH5mWkQQmD7gDG4eOI93ie5JBpdbly5tRZvzBoUd0G17niSqvejEcWaZkrQhRCXCSEOCSEOCyFYdkVhc86dj00/nIK3i8fggqwBEELgRPoZ+M34hX2e2y6Bm3c0qDDK/qlu8/1Pls1mSa80MZMSQhgBPA7gEgCVAD4RQvxdSnlA3ZHplx73NWRuPg5ftwhn7tnW1eooc+58WN8B4Orbka+hzRP7QfZTjskDoG91JpvNkl5pIqQATAJwWEp5FACEEH8FcCUAhlQUqN1nL5pcmTl9qvcyTVVodPWdaWSYgi8kaC3MF4xox6HWVNXuRyOKNa2E1FAA33T7fSWAySqNRffU7rMXTKBgEDYrUjZtgKg7BZk1GM4QmsI+OS0DV26tRXu3yVSS8D4ebBxaC/OhZqn6/WhEsSSkVP9gAiHEPACzpJQ3dvz+egCTpJS/7XxOQ0ND10DLyspiP0gdWfC5CXsa+y4PnZfhxrpz2lQY0WknHAILvzShsvX0LGdYqgePfa8NRbaDGPPXx2B0nR5ja1YODl+3CK7MHF8v12VPnQH3Hk6Gvd0AS5IHd49x4byswMt9dx1KxtvVfT/8vz+oHWvOcoZ5ZUTkS1FRUdevMzIyRO+va2UmVQlgeLffDwNg9ffk7helpLKysqi9thZ0Xt8oay32NDr6fH1UdjqKitQ98O/B7bWobO05tspWA9457sIP//4EDK6eIZpaV40z92zrWuLz93dYBODqSeGNpelwNTrPpuru4/okpOTlqzJ70fv/o4D+r1Hv1wcoe41aqe77BECREKJQCJEC4GoAf1d5TLrVvRlrJ63sa/jqN1hcdwArN/0Ohra+wQoAor7G5+P95a+tVKsHUWnSS0R9aWImJaVsF0IsBLAV3tKlDVLKL1Uelm5poc+eP72DobjuAN77bAVS4H9pTmZmR2UspRMt+EeFA60+qrujWfKttWINIjVpIqQAQEr5FoC31B5HolC7z54/nf0GPbaT+EP5Rvy4+uPAAWUywxmlVkcFlmRcOMSELZV99+miVfKtxWINIjVpZbmPCIA3GN48z4mdB1biOtsupMl2v8/1mMxwLFoR1SPfV07JiOnSaCzPACOKB5qZSRF1Ktz6HJLt3wZ8jseUCsfyZ6IaUEDsl0bVPuGYSGsYUqQ5ou5UwK9LgxHH/nM57jyYipN7qqO+bxPLpVG1Tzgm0hqGFGmOzBrs+3GjEZ6sHFRcuxizyoeh3H662q9z3ybeqX3CMZHWcE+KNMc5dz48vZbxPLn5aFn5PBxr/oo7m0bpdt+mc3lx3igzpuelYN4oM4smKKFxJkWaI3Pz4Viy2tv+qKNRbPf2R3rft9Fq5SWRGhhSpEkyN9/vMe/ctyFKHFzuo7ij5Y4ZRKQszqQo7gQqCy+rUnt0RKQkhhTFJe7bECUGLvcREZFmMaSIiEizuNwXB9gVm4gSFUNK4/TQFTuSI9+JiACGlOYF6oodD4UDhoP7YH54GUS3AwuNRw7AsWQ1g4qIguKelMbFc3cFYbMitVdAAYChY2alFRV2F27aXovLt1Tjpu21qLC71B4SEXXgTErj4rW7grBZkbrq1pgf+R4uPSynEukZZ1IaF4/dFYTNCvNDi2E85f/O2mgd+R4uHjJIpG2cSWlcrA/dU0LKpg0w2Kx+vx7NI9/DFc/LqUSJgCEVB+Ktu0KgQws9JjNaF63AMXMOlm+vVb2s3t9yakWTGxV2l6Z/GCBKBAwpnfJ1b1Ws+Du00D04D61L1+KYOUcz+0C+DhkEgONNbly1tYZ7U0Qq456UDnUWA2w86sCOKic2HnXgqq01OOEQMXl/f4cWti5dC5mbr6l9oM7l1BHpfWdU3JsiUh9DSof8hcC647GZOHceWuiaOhPtY8+Fa+rMHvdFaW0fqMCS7DOkAO5NEamNy3065C8Eqp3BfyZRqgVTvB1aqMUxERFDSpf8feDmpHgCfl+s7hnytQ+kdlm9FsdERFzu0yV/91YtGNEe8PtitVfUuQ80b5QZ0/NSMG+UWfUCBS2OiYg4k9Ilf/dWOasCh00s94q0WFavxTERJTqGlE75+sANdrQ692WISGsYUirR4hlRgfZleNwGEamBIaUCrTY19bdMONJRDfNDi3u0OuJxG0QUCyycUIGWbmbtrXOZ8B+zc7B+xiAUWJJ99uLT2nEbRKRPDCkVaO1m1mD89eLTynEbRKRfDCkVxFuBgr9efFo5boOI9IshpYJ4OyPKXy8+rRy3QUT6xcIJFWj5jCjDwX0wPb0SotkOOcCCthvvgGfsBDiWrPZW99XXQGZms7qPiGKCIaUSLd44avzwPaSuW46uXumOZpgfvA2O29fAM3aC3158RETRwuU+AuA98j31qRXofZiH8LhhenqlKmMiImJIEQAg5aXHIDy+qwtFS1OMR0NE5MWQIgibFUn7d/v9ukxLj+FoiIhOY0iRtyDC5fT5NQmg7cY7YjsgIqIODCnye7OuBNC6oBSesRNiOyAiog4MKfJ7s277ucVwT50Z49EQEZ3GkCL/N+teu1ClERERefE+KYLMzefNukSkSQwpAuANKt6sS0Raw+U+IiLSLIYUERFpFpf7YkCLR8UTEcUD1UNKCDEPwD0AxgKYJKX03/ogDmn1qHgionigheW+/QDmAvhA7YFEg5aPiici0jrVZ1JSyoMAIETv/tv6oMZR8cJm9ZaT152CzBrMcnIiiluqh5TexfqoeGGzwvzQYhhs1q7HjEcOwLFkdVTej4gomoSUMvpvIsR7APJ8fOlOKeUbHc/5J4DF/vakGhoaugZaVlYWjWFGxQmHwMIvTahsPb2yOizVg8e+14ahZuX/7AtefxqD9n/U5/Hasyej4qobFX8/IqL+KCoq6vp1RkZGnyW1mMykpJSKNoDrflFKKisrU/y1iwC8WeiK2VHxqe1tPh8f6PZ2OY/Wn50valQ1RuPvUEv0fn2A/q9R79cHKHuNXO6LgVgeFe+vWazMzI7J+3diVSMRKUH16j4hxI+FEJUApgJ4UwixVe0xxTO/zWLnzo/pOFjVSERKUH0mJaV8DcBrao8jHvmr4vPbLLYhdnt5alQ1EpH+qB5SFJlAVXxaaBYb66pGItIn1Zf7KDIpmzb0CCgAMHTMrLSgdKIFhZaegVRo8RZPEBGFijOpOOXvyHdRXxPjkfhWYEnG67OyY1bVSET6xJCKU1qp4gskllWNRKRPDKk45Zw7H8YjB3os+alRxRcJdoUnolAxpOJA2FV8Gsb7p4goHAwpjdN6FV+4At0/xaVBIuqN1X0ap/UqvnDx/ikiCgdDSuO0XsUXLt4/RUThYEhpXDxU8YWD908RUTgYUhqnlV58Sum8f2reKDOm56Vg3igziyaIyC8WTmhcvFbxBcL7p4goVAypOBCPVXxERErgch8REWkWQ4qIiDSLIUVERJrFkCIiIs1iSBERkWYxpIiISLMYUkREpFkMKSIi0iyGFBERaRZDioiINIshRUREmsWQIiIizWJIERGRZrELOoWkwu7C8r12nGxxY0ia95BCngFFRNHGkKKgKuwuXLW1BuV2d9dju6udPKyQiKKOy30U1PK99h4BBQDldjeW77WrNCIiShQMKQrqZIvb5+NVfh4nIlIKQ4qCGpJm9Pl4np/HiYiUwpCioEonWlBo6RlIhRZv8QQRUTSxcIKCKrAk4/VZ2Vi+146qFjfyWN1HRDHCkKKQFFiSsX7GILWHQUQJhst9RESkWQwpIiLSLIYUERFpFkOKiIg0iyFFRESaxZAiIiLNYkgREZFmMaSIiEizGFJERKRZQkqp9hhC0tDQEB8DJSKiiGRkZIjej3EmRUREmsWQIiIizYqb5T4iIko8nEkREZFmMaS6EUI8JIT4SgjxuRDiNSFEptpjUpoQYp4Q4kshhEcIcb7a41GKEOIyIcQhIcRhIcQdao9HaUKIDUIImxBiv9pjiQYhxHAhxPtCiIMd/3/eovaYlCaESBVCfCyE+KzjGu9Ve0zRIIQwCiE+FUJsVuL1GFI9vQvgbCnlOABfA1im8niiYT+AuQA+UHsgShFCGAE8DmA2gLMAXCOEOEvdUSnufwBcpvYgoqgdwG1SyrEApgD4jQ7/DtsAXCSlHA9gAoDLhBBTVB5TNNwC4KBSL8aQ6kZK+Y6Usr3jt/8GMEzN8USDlPKglPKQ2uNQ2CQAh6WUR6WUTgB/BXClymNSlJTyAwC1ao8jWqSUJ6WUezt+bYf3Q26ouqNSlvRq6vhtcsd/uioKEEIMAzAHwNNKvSZDyr/5ALaoPQgKyVAA33T7fSV09gGXSIQQIwGcC+AjdUeivI6lsH0AbADelVLq7Rr/COB2AB6lXjDhjo8XQrwHIM/Hl+6UUr7R8Zw74V1+eDGWY1NKKNeoM31uAITOfkJNFEKIdACvAvidlLJR7fEoTUrpBjChY7/7NSHE2VJKXewzCiEuB2CTUu4RQvxAqddNuJCSUs4M9HUhxC8AXA7gYhmn9fnBrlGHKgEM7/b7YQCsKo2FIiSESIY3oF6UUm5SezzRJKWsF0L8E959Rl2EFIASAFcIIX4IIBXAQCHEC1LKn/fnRbnc140Q4jIASwFcIaVsUXs8FLJPABQJIQqFECkArgbwd5XHRGEQQggAzwA4KKVcq/Z4okEIkdNZMSyEMAOYCeArdUelHCnlMinlMCnlSHj/DW7rb0ABDKneHgNgAfCuEGKfEGKd2gNSmhDix0KISgBTAbwphNiq9pj6q6PYZSGArfBuuL8ipfxS3VEpSwjxMoAPAZwphKgUQvyH2mNSWAmA6wFc1PFvb1/HT+R6MgTA+0KIz+H9wepdKaUiZdp6xo4TRESkWZxJERGRZjGkiIhIsxhSRESkWQwpIiLSLIYUERFpFkOKiIg0iyFFpCIhxM+FELLj+IZkIUSREKJFCGEXQoxWe3xEaku4tkhEWiKlfEEI8SMAPwVwB7xdCMwAfi2lPKLq4Ig0gDfzEqlMCDEIwBfwdiQQAN6UUl7e8bX7AXwfwLcAbmC7Lko0XO4jUpmUshbe83c6u7mvAQAhxNkARksppwN4D97jY4gSCkOKSGVCiEIAiwC0djy0tqMj+HScPtNsC4BpKgyPSFUMKSIVCSEMAJ6Ft7HxQgCvwHu0+L0AsgA0dDy1AcAgNcZIpCaGFJG6lsA7Y3pPSvkMgN/Ae2rr7fAu/2V0PC8DOj4+nsgfFk4QaZQQ4hwAy6SU1wohfg3AJKV8VO1xEcUSS9CJNEpK+YUQokII8S94Z1c3qD0moljjTIqIiDSLe1JERKRZDCkiItIshhQREWkWQ4qIiDSLIUVERJrFkCIiIs1iSBERkWYxpIiISLMYUkREpFn/H2CIRjnOxrRMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mvn = scipy.stats.multivariate_normal(\n",
    "    mean=np.ones(2), \n",
    "    cov=np.array([[1, 0.8], [0.8, 1]])\n",
    ")\n",
    "\n",
    "X = mvn.rvs((100,), random_state=np.random.RandomState(0))\n",
    "\n",
    "num_components = 1\n",
    "X_reconst, mean, principal_values, principal_components = PCA(X, num_components)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "# eig_vals, eig_vecs = eig_vals[:num_components], eig_vecs[:, :num_components]\n",
    "ax.scatter(X[:, 0], X[:, 1], label='data')\n",
    "for (princial_variance, principal_component) in (zip(principal_values, principal_components.T)):\n",
    "    draw_vector(\n",
    "        mean, mean + np.sqrt(princial_variance) * principal_component, \n",
    "        ax=ax)\n",
    "ax.scatter(X_reconst[:, 0], X_reconst[:, 1], label='reconstructed')\n",
    "plt.axis('equal');\n",
    "plt.legend();\n",
    "ax.set(xlabel='$\\mathbf{x}_0$', ylabel='$\\mathbf{x}_1$');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also compare our PCA implementation with the implementation in scikit-learn (a popular machine learning library in Python that includes implementation of PCA)\n",
    "to see\n",
    "if we get identical results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "difference in reconstruction for num_components = 1: 2.643993664852573e-28\n",
      "difference in reconstruction for num_components = 2: 4.477357638322725e-28\n",
      "difference in reconstruction for num_components = 3: 7.867635675115741e-29\n"
     ]
    }
   ],
   "source": [
    "random = np.random.RandomState(0)\n",
    "X = random.randn(10, 5)\n",
    "\n",
    "from sklearn.decomposition import PCA as SKPCA\n",
    "\n",
    "for num_component in range(1, 4):\n",
    "    # We can compute a standard solution given by scikit-learn's implementation of PCA\n",
    "    pca = SKPCA(n_components=num_component, svd_solver=\"full\")\n",
    "    sklearn_reconst = pca.inverse_transform(pca.fit_transform(X))\n",
    "    reconst, _, _, _ = PCA(X, num_component)\n",
    "    # The difference in the result should be very small (<10^-20)\n",
    "    print(\n",
    "        \"difference in reconstruction for num_components = {}: {}\".format(\n",
    "            num_component, np.square(reconst - sklearn_reconst).sum()\n",
    "        )\n",
    "    )\n",
    "    np.testing.assert_allclose(reconst, sklearn_reconst)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA for MNIST digits\n",
    "\n",
    "Once you have implemented PCA correctly, it's time to apply to the MNIST dataset.\n",
    "First, we will do some preprocessing of the data to get it into a good shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-55fb94072c74b006",
     "locked": true,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "## Some preprocessing of the data\n",
    "NUM_DATAPOINTS = 1000\n",
    "X = (images.reshape(-1, 28 * 28)[:NUM_DATAPOINTS]) / 255."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "Image data of dtype complex128 cannot be converted to float",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-17-443d0a611d8e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfigsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnum_images_to_show\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m3\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0mnum_images_to_show\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m28\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m28\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"gray\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mimshow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mreconst_images\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"gray\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m;\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/__init__.py\u001b[0m in \u001b[0;36minner\u001b[0;34m(ax, data, *args, **kwargs)\u001b[0m\n\u001b[1;32m   1563\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0minner\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdata\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1564\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mdata\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1565\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mmap\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msanitize_sequence\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1566\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1567\u001b[0m         \u001b[0mbound\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnew_sig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbind\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0max\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/cbook/deprecation.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    356\u001b[0m                 \u001b[0;34mf\"%(removal)s.  If any parameter follows {name!r}, they \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    357\u001b[0m                 f\"should be pass as keyword, not positionally.\")\n\u001b[0;32m--> 358\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    359\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    360\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/cbook/deprecation.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwargs)\u001b[0m\n\u001b[1;32m    356\u001b[0m                 \u001b[0;34mf\"%(removal)s.  If any parameter follows {name!r}, they \"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    357\u001b[0m                 f\"should be pass as keyword, not positionally.\")\n\u001b[0;32m--> 358\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    359\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    360\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mwrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/axes/_axes.py\u001b[0m in \u001b[0;36mimshow\u001b[0;34m(self, X, cmap, norm, aspect, interpolation, alpha, vmin, vmax, origin, extent, shape, filternorm, filterrad, imlim, resample, url, **kwargs)\u001b[0m\n\u001b[1;32m   5613\u001b[0m                               resample=resample, **kwargs)\n\u001b[1;32m   5614\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 5615\u001b[0;31m         \u001b[0mim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_data\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   5616\u001b[0m         \u001b[0mim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_alpha\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   5617\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mim\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_clip_path\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/conda/lib/python3.7/site-packages/matplotlib/image.py\u001b[0m in \u001b[0;36mset_data\u001b[0;34m(self, A)\u001b[0m\n\u001b[1;32m    692\u001b[0m                 not np.can_cast(self._A.dtype, float, \"same_kind\")):\n\u001b[1;32m    693\u001b[0m             raise TypeError(\"Image data of dtype {} cannot be converted to \"\n\u001b[0;32m--> 694\u001b[0;31m                             \"float\".format(self._A.dtype))\n\u001b[0m\u001b[1;32m    695\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    696\u001b[0m         if not (self._A.ndim == 2\n",
      "\u001b[0;31mTypeError\u001b[0m: Image data of dtype complex128 cannot be converted to float"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAwMAAADBCAYAAACT6mwBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOy9WXBc13U2+p2e53kG0JhHEpxAUgwlkposyopcYlJ2FMcPrtR9SJXrPiUPzlte/1Sqkjzcl5Qrrgyu3Pgqjp1YiSXHkiKJEidxJgEQEwF0o9HzPE/nPtBr63SzMRDE0CDPV8Xi1N04Z/c+a6/hW9/ieJ6HCBEiRIgQIUKECBEinj9I9voCRIgQIUKECBEiRIgQsTcQgwERIkSIECFChAgRIp5TiMGACBEiRIgQIUKECBHPKcRgQIQIESJEiBAhQoSI5xRiMCBChAgRIkSIECFCxHMKMRgQIUKECBEiRIgQIeI5xVMFAxzHvclx3AOO4+Y4jvvz7booESJEiBAhQoQIESJE7Dy4rc4Z4DhOCmAGwDcA+AFcA/BdnucnASCVSokDDESIECFChAgRIkSIaCMYjUZO+PenqQycBDDH8/wCz/NlAP8K4J2nuTgRIkSIECFChAgRIkTsHp4mGOgA4BP83f/bfxMhQoQIESJEiBAhQsQ+gOwp3su1+LeW1KDZ2dmn+DEiRIgQIUKECBEiRIjYKgYHB9f8v6cJBvwAugR/7wQQWOvFJ06ceIofJeJZwrVr18T9IKIB4p4QIYS4H0QIIe4HEc0Q98STI5lMrvl/T0MTugZgkOO4Xo7jFAD+EMB/PsXniRAhQoQIESJEiBAhYhex5coAz/NVjuP+bwAfApAC+DHP8/e37cpEiBAhQoQIESJEiBCxo3gamhB4nv9vAP+9TdciQoQIESJEiBAhQoSIXYQ4gViECBEiRIgQIUKEiOcUT1UZeJbAca3EkYCtDmV7ntBq7cR1E7EZCPcO/VncOyJEiNhNrHX+A6I9EvF84LkNBpRKJTQaDeRyOex2O8xmM1QqFWw2G9RqNdLpNGKxGIrFIqLRKBKJBKrVKnK5HKrV6l5f/p5CIpHAYDCw9RodHYVOp2P/XygUMD09jWAwiGKxiEwmg1qttodXLKJdQHtHq9XCZDJhaGgIJpMJAPDuu++iVCphaWkJ4XAYlUoF2WwW1WoVtVrtuX/uRIgQ8fTgOA4KhQIymQwmkwldXV1Qq9UwmUzMFhHS6TRCoRAKhQICgQCCwaAYHIh4JvHcBgNqtRpOpxM6nQ5Hjx7F6OgorFYrDh8+DKvVisXFRdy7dw/JZBI3btzAgwcPkMvlUC6Xn3unRCaTwWazwWaz4ejRo/j+97+Prq5HKrMcxyEUCuEf//EfcfnyZSQSCRSLRTEYEAEAkEqlcDqd8Hg8GBwcxPe+9z309fUhGAzihz/8IZLJJH71q1/h2rVryOVyCAQCyOfzbA+JB7EIESKeBhKJBBqNBmq1GoODgzh//jzsdjsGBgbQ398PqVTKXjs/P48rV64gFovhs88+QzgcFs8yEc8knrtgQCKRgOM4lgkwGo2w2Wyw2+2w2WxwuVyw2+0olUqIRCJQKBSw2WyIRCKQSCSQyZ67JWOQSqWQSqVQqVQwGAwwm82wWCxwOp1wuVzgOI79slqtMJvNqFarUKvV4Hn+mczuchwHiUQCqVQKpVIJuVy+4Xuq1SoqlQp4nkelUnnuDheZTAalUgmtVgu73Q6Px4NgMAiXywWVSgWHwwG73Q6lUolsNgupVAqe51EoFPb60ncMZJdoP3EcB7lcDrlcvi6FoRkULPE8j1KphHK5DJ7nUa/Xd+rSdwy0DhKJBBLJo/Y2mUwGmUzWYGsqlQqKxSKzMWLA+Dg4joNUKoVCoWD2itaXUKvVUKlUUK/XUa1WnxlbTftELpdDoVBALpfDbDZDrVbDbrfD4XDA4XDA5XLB4/FAIpGA53nwPI9cLgen0wmJRAK9Xg+FQsHW5nnYZyqVCgqFosEuEda6f1o7nudRLpdRqVQAgNmg52HdADSsmfCXXC6HVCpt+Ddar3q9znyCer2+q/bsufJsZTIZjEYjVCoVDh06hDfeeAMWiwUdHR1wOp1Qq9XQarUAALPZjLGxMRSLRTgcDhw/fhxzc3P46U9/ikwms8d3svvgOA52ux0dHR3Q6/U4evQouru74fV6odVqGzasWq3GkSNHoNfrsbq6CqfTiVQqBb/fD7/fv4d3sb3gOA4GgwEGgwFWqxXnzp1DX1/fug8vz/MIBAJ48OABMpkM5ufn4fP5nhsD2QwyggCg1WohlUrx4osvoq+vD8lkEjMzM0ilUrh79y6uXr2Kcrm8x1e8/ZBIJDAajVCr1VCr1bBarVCpVBgaGsLIyMimEhB0SJNDVyqVcOXKFdy8eROlUgmZTGbfrR1RNilxo1Qq0d3dja6uLsjlcmi1Wsjlcjx48ABffvklUqkUEokEUqnUXl9620BIibHb7ZiYmIDFYoHFYoHD4WjIgkciEUxPTyOVSmFxcREPHz7c94kKmUwGtVoNhUKBAwcO4NChQ9BoNHC73cxu9/X1QaPRwGw2g+M41Ot15PN5VKtVqFQqHD58GNlsFpFIBNFoFNlsFisrK8+8H6BUKnHq1CkcPXoUKpWKPYPA+g59oVBALpdDLpfDrVu3MD8/j0qlgnw+z+zTfrNFWwHRqBUKBVQqFVQqFdRqNYaHh+FwOKBUKqHX6yGVSlEqlZidvnfvHsLhMFKpFILBIMrlMmq12o4/i89lMKDX6zE2NoY333yTZSBVKhUAsCjNaDTCYDCA53n09vaiVCrhq6++wocffrjHd7F3sFgsGBoagsViwfHjxzE0NAS9Xg+1Wg3gawOhVCoxMjICl8sFv98PmUyGeDyOarWK1dVVAF+v834HZbd7e3vx+7//+zh9+jSARie3GXfu3MFHH32EUCiETCYDv9//TKzF04Icv4mJCfA8j2g0is7OTsTjceTzedy4ceOZPESkUil0Oh2MRiNMJhN6e3thMBjw6quv4o033mAH8HqgYICy5NlsFhzHYWlpCdlsFoVCYd+tndlsZjxuj8cDnU6HkydPYmJiAmq1GhaLBWq1Gr/5zW8QDAYRCoVQqVSQTqfF5+m3oEykSqWCx+PBmTNn0N3djZ6eHgwNDTUEmvPz8/jwww8RDAZRrVaxtLS074MBqVTKKEHj4+N45513YDKZ0N3dDYvFAgCs6iQMqIvFIorFIqMTVyoVzM/PY3FxEdFoFPF4/JkPBhQKBY4cOYJvf/vb0Ov16OzsZL2B61UaU6kUYrEYYrEYeJ5HJpNBoVAAx3GsUkmV8WcZSqUSdrsdGo0Ger0eRqMRZrMZr776KoaHh6HT6eB0OiGXy5HNZpHNZhEOh6HT6TA9PY1AIIBUKsV8CTEY2AZQuV2v16Orqws2mw1ut5tlDISleKETR3+WSCRQKBSs6Vij0TCqB73uWYVMJoNGo4FCoYDFYmH0H2oCValUzJgSOI6DSqViTaIejwcajQaLi4tQKBTsNft93TiOg16vh9PphMPhgEajgUwmWzcQAACdTgeXywW5XI6uri6Ew2GUy2UUCgVWgiajud9BZVKtVguj0QilUgmv14uuri44nU4oFIrH7rOZAlIqlZ6ZsjzHcZDJZJBIJFAqlSwR0d3dDZvNBqPRCK/XC51OB4vFwrK6m/lcQq1Wg0qlgt1uZxUWslVCilq7QEiP0uv1MBgMAIC+vj5G5XS5XNBoNOxwJfoCJXh6enqg0+lYtq0V6vU6MpkMy1CWSiXU6/V9S6NaDyqViglkmEwmaLVadHR0MPut1WoZ5Yqg0WgYJcbpdMJut7Ms734KJKVSKQwGAzQaDXQ6HTweD7RaLbxeL0wmEzu/M5kMisUiUqkUo0cRLaNaraJer8NqtUKn00Emk8FqtaKnpwd6vR7xeBxyuRz5fB6pVOqZ2j8ymQwKhQJarRZqtRpKpZL5SXS+NZ/5QigUCmg0GlQqFXR2dmJoaAj5fB7BYBCFQgGxWIyJQtDz96yA4zhGF7ZarfB6vTAajdDpdIxFYLFY2NrSmlLATvu1WCxCLpcjlUohm80iHo8jmUxu6Fs8DZ75YEAikcBkMsFiscDj8eCdd97B0NAQXC4XLBYLVCpVQyBQLpdRr9fZL+DrUqPBYEBHRwf7YmKxGCvfPEsbWgi9Xo/R0VEYjUaMjIzg8OHDMBgM6O/vh8fjgVQqfYwnTyoNer0eVqsVnZ2dyOfzSKfTmJycBPAoQCuVSntxS9sGqVSKgYEBvP7666yhWsiXXAtutxvnzp1DNptlahbpdBpzc3NIp9OIx+MsO7efQQ6vVCrF6OgoTp8+zTLfHo8HRqMRVqt1zfeXSiWEQiGsrq4yQ7jfQc6ZQqGAx+OB1+ttoN1ptVo4HA5WlhfSONYDrY1wzU+ePAmLxYLV1VX8+te/xvz8PFKpFCKRSFvtLeq3USgUOH78OCYmJgAAf/zHf8yy/waDAQqFAkaj8TFHtr+/H3/wB3+AQqGAeDzOsmnNKBQKDWIQwWAQ+Xwe5XL5mQm+CS6XCyMjI9DpdOju7obdbofL5cLExASjoRE3nmC1WvHiiy8in89DqVSiXq8jkUjg/v37rKK7H6BSqTAxMYHBwUF4PB5MTEzAbDbDarXCbrejUqkgEAggmUxicXERX375JRKJBPL5PPL5PFQqFXp6emAymTA+Ps6C9GPHjqGnpweRSAQ9PT3w+XyYmZnB5cuXkcvl9vq2tw2UrDKZTA3BN4BNnW9qtZoF6d/61rdw5swZxGIxTE5OIh6P4+7du7hy5QqrvlCi4lmATCaDx+OBxWLB4OAg3n77bTidTmi1Wuh0OtavQolDCqpUKhVkMhlUKhVee+01nDp1CouLi+jv70c8HseVK1dw69YtFkTthK165oMB4NFC6/V62Gw2DA4OYnx8nJUOhREuBQCUFahWqyyTR81XlLkql8tIJpOo1+vPRJZ7LSiVSthsNlitVrjdbrhcLuj1epZtIvA8z4IqatDmeR5arRYWiwXFYhF2u50ZlfUyC/sFFGhSyZnumbDWntBqtdBqtSiXy4jFYsjn88yJkUgkKBaLT9Q02q6gpkW5XA6LxYKRkRFYLBb09PTA5XI10PNagaR8M5kMSqXSM/GMUQM+9QZ0dXXBbDZjZGQEAwMDjP4iDLCf9L6pOc3tdrNM+c2bNxGJRFAul9tub5GNVSgUcLlcGB0dBQAcPHgQVquVNZu3shlE6RweHkatVkM6nUYul2u5Zvl8HslkEvF4HAqFglVM9jsVphkSiQQ6nQ5utxtGo5ElbqxWK5xOJ6u8AI17i6hEtVoNDx8+REdHBxQKBRYWFvbiNrYMuVwOl8uFgYEB9PT04OTJkw1Jh2w2i1wuh3A4jIcPH+Lq1auIRCLIZDLIZDLQ6XTIZDLweDwseKBqCVVMIpEIlEol0un0pkQj9hMo6CanlRpem4PHVs8Y2XwK1HU6HXieRygUQq1WQzQaRTAYhFKpRLVafSb8ACGI8mmxWOByuTA0NITOzk5oNBpotdo1ZzIR+4QqMsCjoCqTySASiWBmZgZSqZRVr3YCz1wwQLQEKlUplUocOHAAo6OjcLlcrFGYOuSFqNVqrCRK2VmO4zA0NISuri6YTCa88MIL6OjowMOHDzE5OYlCoYDV1VXE4/E9uuOnAz28tBmpLOj1euFwOGA2mzE8PAyj0QiPx8MUX8iJq1arKBaLLCiicv9m6Q37BcJGPAoqKVvkcrkYReFJIJFI4HA4MDQ0hHA4jNXVVebQPAtGkgyjSqViFDOLxQK9Xv9YZqQVdDodhoaGYLVaWVNjLpdDMpncd8pC9IxZLBYcOXKElZAHBgag1WrZHpLL5dvy3RNVj1S/SDWlXC5DKpW2VTaOEg46nQ6dnZ3o7u4G8KgqSRmz9QIYCjhlMhlbw1aOilqtxsGDB6HRaJBKpeByuZBMJuH3+zE3N8eqA/sp6CQKHlEMqN+ku7sbw8PD0Gg0LBCgtdnMZ7pcLhw+fBixWAzpdBoKhYI50O1URaGzRiqVoru7G/39/TCZTDh58iT6+/thNBoZLYxsbDabxdTUFILBIPx+P6sKCBXeIpEIKpUKenp6kMlkoNfr2R6j85EcZY7jkMlkEI1G2WyicDiMfD6/18uzadCzo1Ao0NfXh6NHj8JisaCrqwtarZapUG0VarWa9f3MzMywIGG/swOaoVQqMTo6ikOHDqGrq4uxT55UFQ4Aq+oZjUYMDg5iZWWFVTR3ohL17HhraOTjGgwGlsU+d+4cXnvtNVb+0ul0zHEVolqtIp1OI5vNYnp6GpcvX2YHTWdnJ5xOJ37v934P5XIZN2/exKeffopIJIJLly7t22BASCsghQmz2Yy33noLJ06cgEqlgtlsZlQppVLZILFaqVSQSCRYxpG+A6vV+kwFA5RtU6vVOHz4ML73ve+xQ9ZqtUIul7NG6s1CJpOhp6cHnZ2d8Pl8WFlZQa1WYxWC/Q7iTVIg2dnZycr1er1+Q+NoNpvxwgsvMDUFCrpnZmb2VTBAz4RcLkdHRwfefPNN9Pb2wuVyobOzk0kektO7WWrQRj9Tp9NBo9GgVqvB6/WiWCyiXC5jbm5uG+5q+0B8bqoeHTp0CPPz83A4HA3SfGuBmkR5nodKpVrTUa3X63A4HHjxxReRSCRw9+5dxGIxfPHFFwgEAqwa3C6O7mZAdttut6OnpwcXLlxAd3c39Ho9zGYzpFIptFots9ub2VsSiQR9fX1wOp2Ix+MsOFheXkYul2Nc73bgyZMTq1arcfbsWbz77rswGo2w2+0wmUzI5XKIRCLI5/O4fPkyPv/8c6TTadYIXK1WUSgUGu6nVCrh4cOHkEql8Hg8iMVijOah0+mg1Wpx8OBB1Go1jI+P46WXXkIul8ONGzcwOTmJcDiMS5cu7atggKopRqMRR48exVtvvQWr1QqPxwOTybTpvbMWKLFTqVQwNzfHhrvlcrl9tU4bQaPR4KWXXsKFCxdYVXajZMZaoN6CYrGIeDzesJfFYGADUPMG0XnMZjPjJZNMHRnFtb4caiYrlUpIp9PgOA7FYhHVahVSqZRtYqvVCpPJhHK5/MQZ4XYAaXfL5XKW4WieHeDxeFjZiqhSzQaBqimlUqlBR9doNO7Rne0MhOVPrVYLp9MJt9vd0ET9pA487VeiQSiVyg2z5fsJQvoH3SfdI+2j9TKx5MgQf95sNqNer7Ns8W4oLGwVwv0ilUqZ7RFyl8mGCL/vJz00aO1avU/4jNN3sB2BxnaB7LBSqWRKSlQ1ArDpZIKQnkhrudaeIgEDqspJpVImN02Z4XZwcjcC3Ss5wyaTidGAPB4Pk8kmeiut5WYrHyqVitES6PzMZDJQq9WM572X60TfOZ1bOp2OzSwRVh7r9TrS6TTS6TSrDBD1IpFItPxs6h0EwHjtxWIRCoWCVdRo75I0JE0opn4cqoiSlGa7Qnhm63Q6dv5TFZfWcaOAfDM/RyqVolarbems3A+g844ka+n53Oq6EaUUQEMT906t3TMVDOh0OkxMTKCrq4vxtfR6PQYGBmAymR5TT2gGcZt1Oh38fj9rokomk1heXoZarWa6sTSwhJQ79hNoCiwFS/39/ayS0t3dDZ1Oh8HBQeaoUImr1SZMJpO4c+cOq4xQVYYoCs8ShBQo4lQKB0M9iwbuaUBVJ5VKxYKArRgziUSCnp4enD9/HrFYjFEAs9ksQqFQW5WahUNluru7MTo6Cq1Wy5rKnE4nxsbGYLVaodFotuT8kwIHKU8RJaidHP2NIJPJGCd5bGwM3/jGN+B0OjE4OLgrz5FarWaNtdFoFD6fD/F4HHNzcwgEAjv+858WxEvW6/U4c+YMDh8+DJvNxhpfyS5t1fEiO2cwGHDkyBF0dXVhZmYGpVIJsVgMS0tL8Pl8exIQ0H5XKBTo7+/Ht771LaZaY7PZUK/XMT09zZz/mzdvIpFIwO/3Y2lpiSm3bQbpdBqzs7PIZrMwGAzQ6/WsGl4sFuF0OtHX18ekS51OJ0KhEBQKBfx+P5aXlzE9Pd1WNopAFDONRgOHw4Fz586hv78fXV1dTAHwSavdayGVSuHBgwdIJpO4f/8+wuEwstlsW67LVkDPGzVPP23wBDwKRDOZDLLZLJaWltgMkJ1qVn+mggGNRoMjR47g8OHD6OrqwpEjR1jj2WYMIh1Q9XqdySBShWBlZYXNHqAhEhaLBZVKZV8GAzabDd3d3XA6nTh16hQcDgc6OzsxMDDAOJjChqG1skmkEEQHqEQigc1mw8jIyK7dz26B9pFCoWBSYYTtePifNZBjTFkNIRXmST/H6/XCarUiFoshGo2iUCggEokgHo+31YFC+0Mul6O3txevvPIKaxAmShBJFW5lv1AwQNXLUqnEqIz7KRiQSqUwm80wm80YHBzESy+9hI6Ojl2rXqhUKnR0dKBer7NggPbTfggGaHCW1WrFyZMn8dprrzFaQnNfwFYqlvRLp9NhbGwM9XodZrMZ0WgUgUAAxWJxzwZICqWre3t7ceHCBRw4cIBRWZLJJObn53Hr1i0sLS3h008/ZZr3TzoFN5vNsl4lqkAUCgUsLS0hk8ng6NGjrKfOYDBgaGiINcv6fD589dVXmJ+fbysbRSChD7PZjI6ODrzwwgvMZyJq8HYhk8lgamoKgUAAc3NziMfj+3L2yVogGvV2JmXK5TIbouj3+zE/P49CobBjFNlnIhigjD9pmZvNZjY6XGgY6/U6a3atVqvsASVuLQDGh6ThPcViEclkEpFIpEFhSPhrv4CCGGoq6+3tZSVgs9nMGoiaI1uhUlAzyMDWajWmqkSZy/20NhuBpNJorTZjKGngSiaTYQ3IlC1Xq9X7ynnbCCT/KDSKJEFLDaxCipBQgUtI81gL1LhHVQbKDLbbHmtumqb9ImyGpfvfCojCQL0lyWSyIdhqd5BdUSqVcDgc6OrqgtvtZrNMtuKA0CwKsjnCCt56jjC9job/cBwHjUbDkiDt3DtAWXtyUGkPtJr58rQgSgwpXVUqlTWVUXYDMpkMDocDbrebDcKSy+Uol8vI5XJIpVKIRqMIhUKIx+OM5rsVUCWEBpHl83kUCgU2IXZpaQn3799nVQOtVotCoQC9Xg+XywWv14vh4WF2TclkcptXY+uQSCSw2+3o7+9n/VwbUalb2e1Wr2n+P5JULpVKzB5KpVJkMpm2pXo+CUgS2uFwQK1Wb4t/SNXfWq3G5I93kp6374MBjuPYIIeuri4MDg5idHSUObZClMtlBINBZLNZlg0CgImJCRw4cICNIa9UKixDlM1mUalU4Pf70dvbi/7+flgslrZzQjYDq9XKrv/NN9/ECy+8wJw2crCEgcBm7nEtR26/BUobQa1W4+jRoxgbG8Pg4GCDrOpaqFaruH37Nm7evAmZTMY0mzs6OjAyMrJtJdh2QFdXF06dOsXUp/r6+lgWnGhV1Ggt3GNCw7befiFHsVKpsP6BXC7XdgGVVqvFwMAAUw2amJhgTgLNNHkaGkypVEIqlUKxWMTU1BSmpqZgMpnw6quvrjlsq51AVQybzYbz58/j9OnTsFqtbLbCVmxGKpWCz+dDpVJhCSCaHEtJnrXAcRw6Ojpw5swZhEIh3Lt3D3fv3mVc73YNCCwWC0ZHR2Gz2Zgoxk5zsc1mM44cOYJkMomFhYU9o0VqtVq8/vrreOWVV+BwOGC328HzPMLhMJaXlxEOh/Hll1/i8uXLKBQKT9WgGolEcPXqVajVaiYrns/nsbS0hGQyiRs3buCDDz6AVqvF8ePHceDAAej1evT09GB8fByjo6M4deoUYrEY3nvvPXz88cdts6fkcjnOnTuHd999F3q9nqn9rJWweBKfoPl1JAZRKBQgl8uRSCSQSCQwOzuLYrG4rfe12yAa60svvQSn04murq5tqaqQDaIgNJPJoFwu79iMmGciGCBZNWGzcKuyea1WQzabRSKRQDAYxOzsLFNOoIw2ZZlo8VOpFDiOQy6Xg0ajaShr7TeHV6PRwOVywW63Y2hoCIcPH14z+m928tcyYJu9//WatvcDSG2hv7+fTQ9eD9TcGg6HMTU1BblcjlwuB71eD7Va/UxkQ4SgKoDT6cTJkydx6NChDfcW/X2jw1HYjCuc2vu0cnc7ARoq43Q64XK5WJP5VtG8NlTdzOfzCIVCWFhYgM1m2zdDj6ipTqPRoL+/H0ePHmXTN7f6XZZKJZYBpgpDrVZbd6CdEFRFpmw7HeTtqixEZ57NZoPD4WDUoJ22r2q1mk2C3owa2E5BLpejr68PJ06cgFqtZgFfPp9HOBxmkqFLS0tP/bPy+TwCgUBDozIFA6lUir2OuOKkWHjw4EF0dnbCarWio6MDkUgEFy9e3JS92y1QX9OpU6fYJPj1zvnN+gStoFar0dHRgVqthtnZWTidTkilUpaQ3c/gOA5msxn9/f2w2WwwGAzb8myQD1GtVlEul1EqlXZ0WOSGwQDHcT8G8DaAMM/zB3/7bxYAPwXQA2ARwB/wPN+6NX+HQNkfhUKB4eFhHD58GE6nEw6Hg2UfgUcLSkNGYrEYrl+/jpWVFTZGXCqVYmlpCQaDAcViET6fD5lMBrdu3WIyTpVKBel0GlarlQUDpA1OpTaLxcIak9rJ0VMoFKz8RwEAyWECjRMFaeMBX+uikyMm3NzC91ATIw0QWStTu5/GjpNToFAo0NnZCa/Xy7JipB3czPkul8sIhULIZrNIJpNMC/jGjRuYn59ng8ZInUq4FsLPaTfndiugvdHKIApLzKVSCYlEgtH1bt26BQAsi65SqWA0GhvoQKTV39nZyZy/bDbLJoHvFSj5YDAY0NnZyRyBzRwK1WoVoVAIqVSKlYR5nodGo2GZOpPJxO6VJglHo9F980wRSN2GfgkrRZsF6ZNTE+LDhw9x48YNFAoFtm90Oh18Ph9MJhNsNhvr1yC7JgTZrXbvu1CpVHC5XNBqtRgdHcXQ0BAsFguMRuOGiamt9Om0ygxTUK7RaGAwGFAoFJ6KhvMkIIqrcFaJcKDV5OQkLl26xPqKtgO1Wg2FQoFViUiet1khqFarYWVlBbdu3YLb7YbJZEI0GoVWq4XJZNK36W4AACAASURBVILRaMT4+Dhee+01JBIJzMzMIJ1Ob8s1bhU8zyMWi2F+fp4lV2gKbvOQP6HdbhU00BBWocQ4x3EwGo0N1TmOezQMcWJiAtFolEnVEh27nSajbwRSh1OpVPB6vQ1Uq7VANM9qtcqU9tolSbqZysA/APh/APyT4N/+HMBHPM//H47j/vy3f//h9l/e2lAqlbDb7dDr9Xj55Zdx4cIFaLVa2Gy2hoZeys76/X74fD784he/wOTkJEwmE9xuN1QqFW7evIlgMIh4PI4bN26wgzkcDjMpLI57JCNGJUdSbDCZTOjp6YHH40E2m0UwGGyrYECj0WBwcBB2ux3Hjx/HN7/5TdZXIXToa7Uaq4gAX9MyiAMuPCTpPfV6nXHZSqVSg1JMu3Nu14NWq2WqHOfPn8fv/u7vMseMGhybH+JCoYA7d+5gcXERs7Oz+Oyzz9hwrEKhwDSDKTgV8pv3Y/9JKwi/8818/7lcDnNzc0gkEnC73Xj//fchkUjQ2dnJpl6PjIw0VGFkMhk6OzshkUhQr9eZlnixWNyz2QPEgSe1sbGxMQwMDDAHdCMUi0VMTk5idnYWhUIBiUQCtVoNXV1d8Hq90Ov1GBkZgUajQTQaxZUrV9hQpHayNZsBSRmTUgkdnJvd+9Q8nU6ncfv2baysrODOnTv46KOPmPQlBQOkiDYxMYHz588zOmRzMCCcML9dSiA7AYPBgImJCXR0dODo0aM4ffr0mvckxFbup9V7aI0oSHc4HMhms0xOc6eh0+ng8XjYFGCTyYREIoHFxUUkEgl88skn+MUvfsGSf9uBSqWCTCYDAGz+C/UONr/u/v37mJ+fh9PpRDabRUdHBw4fPoxz585Bp9Ph9ddfx+DgIB48eIC///u/3/NgoFarYXl5GVevXoVGo2Hy61arFZ2dnWv2ILWy7blcDvPz86wngnrjBgcHH6P/0VRe2jdSqRTRaBTT09PIZrM7d8PbDJqdYLFYMD4+joMHDz4mLNKMer2ObDaLQqEAtVoNo9HYNvOYNrwKnuc/4ziup+mf3wHw8m///I8A/he7FAzQppLL5WwQCI1+puZCYSMrRZ3pdJo1AodCIdTrdRgMBiYdynEcEokEAoEAgsEgSqXSYxncYrGIWq0GnudZdoTGTOt0OtTr9bbJ7BJ/lKJXm83GHCxqZgS+LkVRBoRKn1R1IVpGc5agUqkwSlWpVEK5XH7s3vdrMED6vmQgu7u7N1SMqlarSCaTCIfDWFlZwfLyckOzGClTCR24VqXX/QJhpp7jvh6qJczyrjW6nv5M2aB4PA63241gMNgwW0BIp6L3SCQSxt8Vzr/Yy/Uj2oZWq4XBYIDRaGQZsc003heLRSQSCYTDYRQKBdawSPMnyCHR6/VIp9OIx+OIxWIso91cnWsXG9QM2ifN/UlPAuLRFgoFxONxhMNhhMNhVpWjM0Cv1zPll2g0imw2y3TzKRAXrhVlvalC0I4CCFQhohkVRqNxQwradgY2tGZSqZQFv+TM7QaUSiV7voiWQ+d3NBpFJBJBJBLZ1qSAcJbJeoE3z/Ms8SOTyRAKhSCXy1m1jyr0tVoNiUSiLRrVeZ5HPp9HLBZDqVRija/lcvmxSr7QptOaCJXN8vk8kskkEokEOwuIdk3vI1D1jud5GAwG1rTcbs/bRiBREZJoJyGaVs69MPFKFHSe59matwONmtvMRvxtMPC+gCaU5HneJPj/BM/zZuF7UqkU++DZ2dntul4RIkSIECFChAgRIkQ8AQYHB9mfjUZjQ/Sxa/WJEydOPPVnkCKLVqvF4OAg3njjDTgcDhw4cABjY2MAHjX8lMtlNrAhn8/jiy++wJUrVxCPx3Hr1i2EQqHHMrJCLtxaAdKpU6fw13/913jhhRcAfJ0NuHXrFh4+fIipqSn8y7/8y7Y0Lj0N5HI5hoaG0Nvbi4GBAXz/+9/H2NhYyyxRJpNhHOT/+q//wj//8z+jUqng6NGj6O/vx/DwML797W+jq6uLvYdee+vWLSwuLuLixYtIJBI4cOAAxsfH4Xa78Z3vfAdHjx5l7ykWi/i7v/s7/OhHP8JPfvITnD17ds9oHRuB5iRYLBZcuHAB3/nOdzZU/nnw4AF++MMf4le/+lXLKaZWqxUvvfQSvF4vJiYm8Pbbb8NkMrHvZHl5GT//+c8xMzOD6elppoTRbiC6HPVA9PX1wWw24/jx43jnnXcYh5KmbxKfmLJmpVIJwWAQ0WgUDx8+xM9+9jMsLi7igw8+wNmzZ6HVavHyyy/j0KFD8Hq9eO211+B2u9nPBr7uP7l48SL+8i//Eg8ePEA6nUYikdiTLJvJZMKf/Mmf4MKFC9DpdPB6vSxj21wlqVQqrA/p3r17uH79OqLRKD799FNMTk42qNgQjYZkk1UqFaNlVCoVNlnY6XTiD//wD3Ho0CF4PB4cOHCggacbDAbx/vvvY3Z2Fvfv38fnn3++J+V4mUyGN998E9/61rdgs9nYgMhWuHHjBo4dO9bwb+VyGV988QWuXbuGQCCADz/8EPPz82tODaZny+l0oru7G1qtFl6vFy6XCw6HA8ePH4fdbmeSwZlMBu+99x4+//xzhEIhfPXVV2tOqd0LHDt2DH/1V3+FM2fOAHhydZcnzUI37136vVqt4tq1a7h+/ToCgQB++ctfYnp6+klv54lw7do1/O3f/i3+9E//FDabjQ0H/eyzz/Bnf/ZnuHXrFstStwOo6vTyyy/ju9/9LhwOB8bGxtDT04P79+/jb/7mb9jwLb/fvydceWp8NRgM8Hq9ePfddzE2NgaHwwGv1wulUtnwWtoPoVAIi4uLyGQyuHfvHubn5xGLxXDnzh2k02m8+OKLePvtt2GxWDA2NrbmMx6LxfCTn/wEX3zxBcLhMO7evcsoWZvBtWvXtsWv3CpOnTqFv/iLv8DBgweh1Wqh1+tbVmV5nsfq6ipWV1cRCoXwb//2b7h+/TrGxsbw1ltvweFwoL+/H319fQ3vp+F+8Xgc//AP/4D33nvvqffJetK2Ww0GQhzHuXmeX+U4zg0gvMXP2TQ4jmNKHUQLslqtzPmgxoxoNMqCAOJsUwmZmlWAxw3cVq8J+Fq/f69LXdRULdTwpoYzmUzWcnopjWmPxWJsqiI5JNVqlVGBqHGKqB3hcBiRSIQ1DZHKSTqdhlarZTxummC8n6BQKNhIdnLq1tondADReq5VSqZSoEwm25d0Krp2mUzGVEzIuSLqglKpbGjCrNVqyGQyTA4zEomgWCwiEokgkUhgdXW1YQol/Qy9Xs9Kr632jpCiJPx9t0GUElI0oWmerUq+9Nxls1ksLy+zxtdAIIBkMskkjIk6RFQ8iUTCZpsQPYaoZu3CNd0IZAOIvrPedytEOp1m/Uj5fB6lUglLS0tYXV1ldJCNqBvUbJzJZFCtVhvoVfSZpBokkUhgNpvR2dkJ4JHKEAWzzQ2juwlyLOn5o+erlXPfSqp3PXtDFA/q/SL7L2zuprkZws+lRt5iscjoEbROO2XP6Fmja6GfR7a3nSDsbfH7/SiXy/B6vcyPMZlMsFqtyOfze+YzSCQSGAwGdHR0wOVywWAwMDWuVnuG+gRTqRRWVlaQTqcRDAaZH0DN5Bz3aFhdqwF4zZ8pTNq2SyDXCmTriaaqUChYYprOPQq4aU8CX9NBo9EolpaWEIlEEIvFkM1mWQJLoVAgn8+zWU3CBLVQ3IDoizslxrLV0+Q/AXwfwP/57e//sW1X1AJkDPV6PU6cOIGxsTF4PB6m4a1Wq1EulxGPx/HJJ59gdnaWbVRyUNPpNBsitp2QSqXQaDRsaNdeqlFYrVYMDAzAZDLh7NmzOHHiBIxGI+x2O4DGIRY0AMXv9+M3v/kNm3BXLpcb7qFUKiEajUIikcDn8+Hhw4eIx+P47LPPMDU1xXSc6/U6wuEwm9EwPz8Ph8MBnU7Hfv5+gcViwYsvvgiv14vBwcF1DRoFStlstuWBRHtX2Di510HjVkADhzQaDc6cOYOTJ0+yhj4afkQqS7R/8vk8ZmZmEI1GMTU1hcuXLyObzTLHI5/PY3V1lT2TlKE9cuQIXn/9dSZh2K7Q6XSsD8fpdMJqtbKDQQhySKvVKqanp/FP//RPWFhYYNLFJI1JjpSwilCv15kClUQiYYpfPM839LZQNaYdB4+p1WrY7XbWcHfo0CHWX7Eebt68iZWVFQQCAdy/fx+ZTAaRSATRaBTFYhGxWGxTPz+bzcLv90MqlWJ1dRVKpRL9/f3o6elhQQpVYo4fP46+vj5MTU0hEolArVYjkUggGo3uibNC/SitFF6aQYIOQtELsj3N5xI5FYVCgU2Qj0QimJ6eRi6XYwG/2WzGiy++iL6+PvZeiUQCj8cDhUIBh8OBe/fusX6zUCi0Y4ETzSvR6/UolUosiG63QECI+fl5vPfee7Db7bDZbBgfH4fBYMCxY8fgcDhw5coVLCws7EmwKZfLcfLkSbzxxhtsSrrNZms5/K9SqTBbdfv2bfznf/4nk2iPxWJMcbFer0Or1WJoaAg2m23dnpZKpQKfz4d79+6xc7RdodfrYbPZmCBLR0cH+vr6mHQ09T9VKhWEQiEkk0nWW1GpVHDx4kV8/PHHyGQyLNuv0Whw584dJoTT29vLgnAKAqhvkwbqFotF5HK5Hdkvm5EW/X/xqFnYxnGcH8Bf4FEQ8P9xHPd/AVgG8J1tv7LGa2DNsJ2dnRgdHYXT6URnZyfT/q9UKqyj/fbt24jFYlhZWUG5XGYZlZ3I8HAc1zaVARpRb7PZmNwqHXZAo25tNptFPB5HKBTCzMwMHj58iGQyiVqt1nBw0GwGpVKJQCCA6elp5uwvLy83lGbJ0ZNIJIjH46wx22Kx7Ml6bBVarRbd3d0YGBiAzWZbtwm0Wq2iWCyiVCq1zFJSpE/RvVCxZD9BOOG7v78fJ06cYJNbhdUTYcaCKnV+vx9TU1O4ePHiuvQUUoFxu93o7+9v+zUi58RkMjGpQ2EwJAQ5atFoFJcuXcK9e/c2/Pz1aA9CR49+0c9utwZiarSjYM/lcm1qAvfKygpmZ2cxOzuLzz//HKlUimXd1qIGtQIFn0LI5XIkk0nkcjkm5SqXy9HZ2QmPx4NarQaHw4FkMsmynXsFmsOgUqnWXTNhxZeaoenfW72WspjhcBgLCwvw+/24dOkSkskkU7JyOp04ePBgw3uJKkhBlN1uZ82xkUhke29eAGpcViqVjHZI3127IhaLIRaLwWazMSlgpVIJj8cDuVyOubm5PXteJRIJOjo6cPz4cZbYWGtIH1X+c7kcVldXcffuXcTjcaRSqQblJvqOnE4nLBbLutQ0mqIeCoV25P62E0qlktn5/v5+DA0Nwel0MjUvQr1eRyaTQSwWY4yBSqWC2dlZXL9+ne1bEhwJBoOoVCpIJpOsEkyBGFWRKCAgBsxO0Yc3oyb03TX+67VtvpY1oVQqmWqQw+GAy+WC0Whk6hB+vx+Li4uMy0a0ICrxUnRWq9XauhT1tDAYDBgYGIDD4YDNZntMaSWVSiEQCCCXy+HBgwdYXl5GMBhsWC86ZLPZLGKxGHiex40bN6DT6bC4uIj5+XlkMhlkMpnHylXNdBmiPOw3EH0qkUhAr9evacyIC7iwsIDl5eXHpOIkEgmsVissFgucTieGhoZYNmG/UDwISqUSFosFZrOZGUXis28X6LAJBoNYXFxk1QihsRWCehYkEgkePnzIlDt2CxQYEYeZkg7NjmO1WmUZ7dXV1ccc062A4zhotVrY7XY2hdXhcMBgMLSdVj5J9VK2uTlgEZa+KbAGwOw6VXiFdoXet92gIEuj0cDr9TIbFggEdnVvUQCgUCgwMDAAr9eL/v5+GI3Gx15La5LJZFjFiYJDlUqFrq6ux6qzNKcimUyySdbRaBSJRAL5fB75fJ59F63OTCF1QUhf2o2gqV6vI5fLIZFIIJPJ7IszhuYQ3L17l/kiRDfZq0CT6CszMzOs6rtWJj+TyeCrr75iUr6ZTIZRyoTgeR6BQAAXL16EzWZDb28vXC5Xy89UqVQ4cOAAXnnlFcTjcczNzbXl8ESO42C32zE+Pg6TyYTh4WFGj6Wziaq8mUwGN2/exOLiIrMd1WoVi4uLTE2Q7BZR2IvFIuLxOLLZLOr1OnuWKAHH8zysVitcLhcymQzrvdtutL1HQlkIt9uNzs5ODAwMYGRkhJVQKpUKbt68iffffx/xeBz3799n0qGUQRI+bM9yMOB2u3HmzBm43W44nU4mCUoHbyAQwGeffYZ4PI4vv/wS9+7dQ7lcZmUnykQSjahWq2FxcRG3b98G8OgASaVSLBvevJZ0KJHsKHEM9xuKxSJCoRDL/K7ldFSrVUxOTuK///u/EYvFEAwGG/5fKpXC6/VifHwcHR0dePnllzE4OMgqSfsJGo2GORWdnZ2PDffbDmQyGUilUkxNTcFkMsHhcODYsWNrBgM2mw2nT5/G4OAgPv74YywsLOyaY8BxHFwuF1566SU4HA4mP9uqqbNcLmNubo7NEtiOA48qbv39/ejo6MDg4CD6+/uZY9ZOoAZzcmibOcn1ep05ncSplUgkuHXrFm7fvo1sNotMJsOCqJ3MBNN1WSwWHD9+HN3d3ahWq7h///6uUjnUajUL7l599VWcPXsWZrP5McdKKPO8urrKKJ9EHzMYDPjGN77REAzwPI+VlRXcvHkTkUgE//M//4M7d+4wu06Sj+l0Gnq9vuV9k8NCtDiiIu2GY1ur1RCLxeDz+RAKhbYluN5plEol3L17l2WY+/r6YDQaH9Pg303UajXMz8/jk08+QVdXF9xuNxwOR8vXhkIh/PSnP8XVq1eZNGar5Gq9XsedO3fwox/9CE6nE3/0R38El8vV8h71ej3efPNNjI+P49atW/jxj3/ctsFAf38/Lly4AIvFgo6ODlgsFkbTBB6tz+TkJKLRKN5//33cuHGjIcmRz+eRy+UaEqjpdBrT09NQq9U4fPgwotEoDAYDY5mo1WrI5XJWjThw4AAikQji8fiOzKhor1NjDSgUCuh0Ouh0OqZiItTppox3MplEKpXacRUWahRpF5ADolKpYDKZ2BQ80jGmCD6VSiEWiyESiSAYDCIQCKxZQiZONw0jozkEuVxu3XsXPgDttEZPAuHshbWCGQo2qQFbOEmXQI1UNPHZYrHAZDK1/Lx2Bw2fo56H5gFNFEgKqWiUWSQe80YQasinUik2u6MZtLckEgm0Wi0KhcKeUPSUSiWjv6hUqsfK/bQm5XIZ6XQasViMNbJuBxQKBftO6Pd2BDVMWiwWpqsthHDPkPCDxWJBKpVivV6U2NlO0PPdao/J5XIYjUZUq9V150XsBIh6SgOMHA4HOjo6oNVqW/aEVCoV1ogZjUYRDAZZMEA9c0R1orUmKkM0GmVUFuHPp2QONbO3Qisbvxs2X8jF3i8JJxIPoHWm73gvAnchdZWU3lol94Cvz3NKkPl8vg0/P5/PM+rPer6YVCqFyWRCvV7H8vJy21U0haIfer2eMS7MZjNjDJAfkMvlEI/HmW+1urq6IZWxVqsxH4sEbmQyWcMcIprXQD1DNC9lJ9D2wYBEIkF3dzfOnTvHyuEA2KTObDYLn88Hn8/HqC47iWYp0r0GbVTKOFBpmbI0wWAQH374IVPhmJ+fZw/rego5+XyeHcD0+5MoRWznsJvdBk18JcpVs5NH8o6ZTIZRrUg1RwiJRAK3243Dhw+zYW/7DUTpoCZ5oUMn/I7L5TKbrHj79m1MT08jlUphdnaWDfPbbAZP+Iw1qysQha1YLGJxcRGXLl1CIBDA/Px829EFiO6UTCZx6dIlxrNtx+zXToJs1FqZ0EKhgNXVVeTzeVy9ehVXr17FD37wA8zNzbHgabsdvkKhgMXFRcbhbhYKoGuu1Wq7lr2VSCRsGFt/fz/Onj0Lm82GgwcPssbO5ub0UqmEqakpzM3Nwe/3486dO1hdXYVMJoNCoWDBOzlx5KCQzGw6nX5inj/P84jFYojH41hZWcHKygpWV1eRTqd35RkkJ7JSqSAYDLZl03wzhMo5FBi7XC6YTKZd7xmgXi+dTofTp0/j9OnTMJvNsFqtj9F+w+EwkskklpaWNpVk5TgOJpMJ/f39jFMPtA4SqZq1tLQEv9+/p4pdraDT6Rg1b3h4mMkQC6lBCwsLSKfTuHnzJi5dusQS009ir2q1GmZnZ/HBBx/A7Xbj/PnzLemAO419EQz09PTgzJkzMJvNLBgoFotYXV1FMplkwUDzxOCdQLOTu9dOLzXnGQwGxmETlmxDoRD+9V//FV988UWDbOF6m5XneeRyOeTzefb3rUIob7dfggMKBmg2QzPy+TyWlpaQSCRYMNCKxyeRSOByuTA+Ps4qW/sJHMcxOgBJCRK3lL5PYTCQSCSQTCbx8ccf4/3330ehUGiQnn0SR6H584Va6STLtrS0hEuXLmFpaWnXHJHNgud5hMNhXL9+HZFIBFevXmVUjJ3ge7YzhMFAq3kdhUIBwWAQiUQCX3zxBf7jP/4DP/jBD7CwsLBjFUYKBorFItxu92N7h64ZwK5VBki2UKlUoq+vD9/4xjfgcDjgcDhgtVpbTpcmhar//d//RTgcxv379xGNRtn1qlQq5PN5TE5OAvg6m+/3+7G8vMykW58EPM8jHo9jYWGBBQOhUGhNEYXthkQigclkglwux/Ly8r6QrhYGA9VqFWazmSmx7XYwQMo1FosFv/M7v4NvfvObLYUPqtUqotEoFhYWsLS0tOkkhsFgQF9fH6O5NVO1CRQMzMzMwO/3t51d1Ov1OHDgADo6OjA6Ogq73c7OPuARzefevXsIBoO4fPkyPvnkE9Yc/CQ2i4KBTCaD3t5eHDt2rGE42G6hbYOBZlkl0vCmshqVZoSNLLuZqSejKuRs7mZkS2oRGo2Gde5bLJaGqgDwtQrAVqhTW1lPKq01Z3XbHVKplI1GJ818ms1Av4SNequrq4jFYszZbWUAhM40zRcQrodQk5gakNLpNAqFQltUnYRrQo3DRqOxZSaOqDBE1aP5AaT0tdmfJ5fLmeykUEqR1p5kOGnOBVUjdoJGstY1ajQaRn0h6qLQISHqCenbUyWDGge36zqJAkL7hnouiCLYTtjILgjndZBTsJP0D3rmlEolEokE4vE4qtUqc8b3wn6RYp5Go2HJA51Ox3q/WlWlhU6msPeLUKlUmHIcgQJqkrvdjAPfvA5ChSOFQgGFQsF00ncKRFktl8sNNJv9cMYIvy/ga8np3UqS0d6SyWQwm81MCpnsuVQqZeccPYPFYpHRzqLR6KbtuEajYdRYOita2Tw6SyORSFslc+h7USgUTDKbziKhXaUEGF0/NQhvxW4J9/Ze0d7aNhjQarXo6upizROdnZ1sCifwKCqbnZ1lTsFuLWBzpookolZWVthk0J0GNVXr9Xp4PB68/fbbTPvWarUy2b7d5u2T09tK5nCnBmVsF2jy7fj4OPr7+2GxWBqMNPFsi8UiJicn8bOf/Qx+vx/hcJjpK6/FBRcakOZm9kgkgkgkgvn5eVy5cgVTU1NsANVeQ6vV4tixY+jq6sLIyAgrJ1MJs16vs3uLRCIsC764uMiUfZ6EH6/X62GxWDA8PIxTp04x3jNVqjKZDBKJBD766CPcvXsXq6urCAQCu6oootfrcejQIdjtdpw6dQrj4+MwGo3QarXgOI711lQqFSarGo1GmSzmdj2T9XodKysryGQyCAQC6OvrQ7lchsPhwMDAwGNDovbSYdqMw7Pb15fL5TAzM4Pl5WWoVCpGgTt48CC6u7vZNe3m2ikUCrjdblitVjYpmRwqes6ar4fneTYvRthkTahWq6yCLnyPUCloo/3Y/DOlUimjmRiNRvT09LBZDMLBntuNTCYDn8+HfD7PEhN7Pdtnv4Bma1itVoyOjuKVV16BzWZDT0/PY4mDWCzG1NkuXryI69evI5lMbmquB8dx6OjowJkzZ5jK2looFou4d+8ePvroI9aPuNfgOI4NtnO5XBgbG8PQ0BA6OztZwoeemUQigRs3bmBmZoZJhG7VDxWyNvbKT2rbYEChULBmDavVyvjwlK2lYVgki7lbaP6ySPc1nU6zZpCdBsdxUCqVMBgMsNlsGBsbw5EjR1j1ZD196Z1E85TM/UQNUigUGBwcxKlTpxr0loXXT9KXoVAIt27dwtLS0qY/v9U68DyPTCaDcDiM1dVV+Hw+LC4uPvW9bBcUCgU6OzsxMjKCwcFB9PT0PDYoivYY9e6EQiHE4/EtUfaoScrpdKK7u5utGTW0UzAwMzOD69evs6z4bpaXVSoVOjo60N3djZ6eHrjd7gb6F1V7SqUScrkckskk06rfTtvA8zxT96pUKlheXmb8456enpbvoez8XqC5MtD8f7uNcrmMcDgMjuOwuLiIubk5pgq1VzaLKnF03un1+oa91UybA74+f2gCbHPwTdnXTCazbddJwgg6nQ6VSgU2mw0mk4kNxtspFItFJJNJRo0VVnFErA+5XA673c6GZR08eBBWq5UlMYSgHp5YLIaZmRncuXOH9YRtBI7jmOTyRnMGqtUqgsEg5ubmtuUetwPUtKtSqRpULIW9HUSzzufz8Pl8WFhY2Dbfby8Tpm0bDNDhQU4lHSREB8hkMggGgyw7tlOa0zKZDDabjWVAiPNKjmE6ncbCwgImJyexsrKyK1P0iIs+PDzMqgE0RnyvDCMdENQYZbVaYTabodVqIZVK26YE2AzKPlPJVKfTNdAshDQhwkZ7TaVSwWAwQK/Xw2AwsHJ280FZq9UQj8dZc3e7TWCUSqVs0JiQK9kKVFalfhWJRPJU1aDm7GcoFML9+/cRiUSwsrLC6FTtpiRCNKZsNotQKIRQKIREIrHjCmebWeeNeoV2Es0qY8Lvl/S2qSl/Lw/E3c7Mkd3UarVsWKTX64XX6205yZreI1TUslgsLDvfSshgp7Fba5bP5xEMBlGr1WCzTlZTgwAAIABJREFU2QA8slE0+I+oLe1cgd4rCBWkhJSWVueyQqGA2WwGz/NMwUpYBV4PPM8jGo3izp07sNlscLvdMJvNLV8rl8vR09ODI0eOIJPJ7Jr/tB5kMhm6u7vR3d2N3t5eNjldLpczhbyVlRVEo1FMTk6yZNRu0VR3Em0bDFCWmSgnFAwUi0WW1ZmcnGSNnNt9yFEgolarMT4+jgMHDqCnp4dJQ6ZSKfh8PoTDYXzxxRf49NNPGYd3pyGTyTA2Noa33noLVqsVvb29MJlMe575czgcGBkZYRNku7q6GHWoHYMBymJ0dHTAbrfD4/HAbrczfqAwEHiSB53UB8xmM5MEJN6v8HMqlQoWFxdx+fJlhEKhXdk7TwKZTAar1cp0ldeTwZPJZIzrTBxUYHvoYfV6HTMzM/j3f/93xONxzM7OIhgMrkvN2ikIub+t9kWhUMDCwgKzT1NTUy3pGzuBtSR997L0DGy8B1KpFKamphAIBNj3upto/i6bv+OdBCmOdXd3o6OjA+fPn8fw8DD0ej2bOLreNUulUvT29uKFF16A3+9nqj7bjfWuY6PXbBdisRgmJyfhcDjg8XjQ3d0NtVoNt9uNZDLJZB33u1O2EyDlpUwmA5fLhXQ6DY1G03LejU6ng9frhcFgYOchx3GbqgzwPI/Z2Vn8/Oc/h91ux/nz5zExMdHytVqtFqdPn4bVasXMzAx++ctf7nkwoFQqcfLkSbz++uuwWq0YGhqCyWRigUA6ncYnn3yCa9euIRgMMhXLdktKbQVtGwwAX5dFhZUB4iFTcx7xuLfbAEgkEsjlciiVSpjNZrjdbuYQEU2Jfj7pNdP03Z0G9Qw4HA6YTKaGxmqCsCFvN66HlDDMZjOrCAgNzVac6t2AXC6HRqNhhpHUmNaSjt3M9SsUigb1FGr2Fmb0arUao75Qtma3HduNQJWxzehh0/MinMK7XpWqlcpUq4Y64lJms1mEw2HE4/E9yX6uB+GeEM7joOZhUuXajWto/r3d0MoOkN42NcDuJVqt306tJWm9k1KX1WqFw+GA0+mEXC5/LHnQ3DgMPHpmdDodbDYbstkso9KuNT9hM9dEz6FwgF3zGlCAJ5wKvdP2vVwuI5VKQalUsuBaJpNBp9NBr9cjk8m0LWVI6MvsxTWSoIFEImk5zVuI5kTsk/bO0BT5er3OKqLNlV6qDppMJrjdbkSj0bYYlkiytR6Ph1XEqTJCfmc0GoXP50Mikdh2+mczaK3oHNzJ3su9X/0ngHBRaCjWdh8g9LC6XC709PTAbDbj5MmTOHbsGHQ6HSQSCdLpNHw+H27cuMG67akLvB0O4Wq1iocPH8Ln82FmZqaheWy7QeovSqUSIyMjOHPmDCwWCyvj0lCTXC7XwGltl0hao9HAZrPBbrfDZDLBYDAwqsuTgugxBoMBXq+XDSmhA5oMYiAQwPT0NBKJBK5fv47Z2VkUCoUdp5I8KcrlMoLBIBYWFlCv1+H1etd8rUajgcvlYsZUqVSyhlmhseQ4jimQqNVqtj4AcPDgQaZaBDzKskciEeTzeRQKBdhsNkilUkSj0Z29cRG7CtLz1ul0CAaDu+oskZPTLHiw09BoNCyRMzExgVOnTjGK5VoqOXTu1ev1hiC9q6sLcrkcNpsNiUQCPp8Pq6urTDp0o8CAAg9S8JPL5Th48CDOnTsHq9XaMPWYPmt1dZUNrpyfn4fP59vxRv5kMom5uTnWKwQAVqsVZ8+exeDgIC5evIhAINCWVWiJRMI46EKVnb26lmYHn/qyqtUqfD4fbt++jXg8jgcPHiAejz+RKhywfo9QqVRiNr1dpkdTAExnPyXv6NoLhQIb2Or3+7G4uNhSTnyrkMlk7FwkG0R9GrlcDn6/H3Nzc6w3aCewr4IBAA3BAC3UdjnglKkhTv6JEyfgcDhw6tQpHD16lGUoM5kMlpaWcP36dcTj8bYbiV6tVjE7O4vLly8jEAgww7kToEmdOp0OIyMjOHv2bIOmfqVSYRlSYTDQDkETKQdQMGA2m2EwGLZEtxI6FCaTCb29vXA6nbDZbFAqlQ3BRSAQwCeffIJQKISbN2+ygVntEiARSAdarVZDrVavexhoNBq43e6G4Xccxz32HlpznU4Hs9mM0dFR5vwfOHAABoOBcUwLhQKWl5eZAbTb7ZDL5W3VZC3i6UHBgNlsxuTk5K455M2BwG7KVGo0Gng8HphMJhw/fhxvvPEGVCoVo9m1AmVaK5UKSzDIZDJ4vV54PB643W6USiUEAgHcvHkTkUiEOXlrQVj9k8vl0Ol0UKlUOHjwIN544w32nAJfJ+Oq1Sr8fj9u3bqF1dVVNvBspyvjiUSC6bGnUinwPA+r1YqXX34Z2WwWyWQSn3zyyY79/KcBzUZwu93Mju02mhvQhb9on5TLZSwvL+PixYuIRCJ48OABYrHYEyU6m4ON5meqWCyyrHq7+E0UCAhp6UI7RIkpmsK8uLi45epbM+gZVCqVDX2fpVIJiUSCUdLn5uaQz+d3rNK8r4IBKnVRlnm7ypIUCVLJUaFQwG63My1ehUKBSqWCQqGAQCDAFGWSySTS6XTbbGgCKZqQfNxOZEooaCL1F9LjJboNNZDSQ5ROpxuoMO0QDACNZfG1MhnA13svn8+jVCo9dv2kBKJSqWCxWJj0XavppVQtoc9q1+Yj0qamJrJmJ41oAvV6Hblcjum1FwqFNYMbofNBgRhVkahZiyY8VqtVpNNpZhCJcrMXdCqyEeSw0fyFvWzYVyqVUCqVbOggrR1dE80+IcpSKpVqywZLqVQKpVLZQNMDdnbKOyV+6Ds1Go3Q6/W7lrGl54CoqKSM06oiSd9jsVhEOBxGsVhk05FJnIGcOVpLcu5lMtm6swSI3klUSeKI22w2aLVaRnMEvp7tUywWEYvFEAqFEIvF2DO50xl50mIvFouMFVCv15m9IP18lUrVdlRCqoiS3aB72e35KLTPyfGl75bojYVCAZlMhqmgkUjDk1wj2XWr1fpYTwLP88jn8wiHw0w4gIZV7iVNlgJhEhAROuak2BWNRhGNRrftDCI6FjXBU/KQ1ozO12q1yoYD7uQcgn0VDFQqFfh8PgQCAfj9/m1xwiUSCXtAzGYzxsfHYbPZMDQ0hOPHj0OlUqFcLmNqagp+vx+//vWv4fP5EIlEsLq6ynjf7YZKpYJ8Po9isbjtm0c4jGpwcBDf/OY34Xa7MTIyAoPBAKlUyuQV5+bm8Ktf/QrBYBB3795lVYp2LOWuB9p7NLmz+aDR6XQ4fvw4uru7MTw8jOPHj7NKQ7MTTdWS3Wos3SpUKhWGh4dx4sQJOJ1OdugSKpUKYrEYCoUCbt68iQ8++ADRaBSzs7Osqar5e6Z+F6fTiZ6eHrzyyivo7u4Gz/P/f3t3Ghvndd4L/H+4DJcZksNtONwpUaRkSZQiOZZkO3FcB43TII0ToDGcAm1aFOht0RS3wC3Q3H5o+6VAcHFToBcoWtziBm1ve7MAaWLDjePEslNFkrVRGxdxX4azcPZ9OAvJ934gz9HLESmRImc41Px/gGByOCZfDt953/Oc85znwWuvvaZK4AEPOzzabDbYbDZMT0+rgUg+ydSvmpoadHd348yZM+jv70dra2veZ/j0g9ienh709vbCarXik5/8JJ577jkVEOh7WNhsNty8eRMjIyMIhUIFNUgCHlYvKSkpQWNjo6pNXl5enrP3R2lpqXqtDh8+rNJ0ZGCaa/q8bFnKUKbr6K2srGBxcVH107l16xbcbjfa29tx7NgxVFRUqIGaHCSXlJSoAEfuH9hqVa+8vByHDh1Ce3s72tvb1YbOzs5OtUdODhjlBt5gMIjLly/j6tWrKj98L3tobEXuGSgvL8fMzAyGhoZgMpnQ3t6O6upqnDlzBm+++SZ8Ph8uXbqE8fHxnB3LTpWVlaGjowMnT55ERUUFIpGIagaXjxXh6upqnD59Gj09PTh+/Dja29tRX1+vrunJZBKzs7Pwer0YGhrCyMgI/H7/jqt7CSFw+PBhvPHGG2hqakJ7e/sj+4PGxsbw3nvvIRwOw+l0qsme/SygISt5ybGf7Csg31+zs7N477334PV6d1RS/HEqKytVGfNPfepTeOWVV9DQ0IC2tjYAG4PfWCyGcDisUrly4UAFA8vLywgEAnA6napj5G7JtAWz2YyWlhacPHkSnZ2dOHToEI4dO4bS0lKMjY3B6XRicnISH374YUHVxd2MjCh30xHvceRrVltbi7a2Njz//PPo7u6G2WxGVVUVNO1hUxuPx4PBwUEVQMklrkJLiXmS5eVlBINBde5lD3IrKirQ09ODgYEBHD58GD09Paipqdn0e8k0t1xvPtqtsrIyWK1WHD58GCaT6ZENXjJtLhKJwGaz4c6dO6rB1lYDTjlDVlNTg+bmZhw7dgz9/f24e/cunnvuuQ3PTSaTsNvtmJmZgc1mw8zMzL7NHslZPTlQ6u3thdFozPumN31qi8ViQV9fn9rflF29S5ZfdjqdmJubw/T0dN6KHOyErES1srKius3Lx3MVDMhVr6qqKjQ2NqK3t1elq+WDPqjTz+JnW11dRSQSgcvlgt1ux40bNzA/P4/+/n41q59NpjzIGegnVQFrampCd3c3+vv78eqrr8JqtT6STy5ndOfn5+F2u/HgwQMMDQ3ldcO3PHdlMYGFhQVYLBZ0dnbCaDSiq6sL586dw+LiIkZHRwsqGCgpKUF9fT06OztVo1K5ly4f90J9zxh5rzaZTBtShHw+HxwOB5xOJxwOx1PtNRRCoKWlBWfPnt20z8Dq6ioWFxdx69YtdT8NBoN7+as+FdnkVm7gr6+vV+f/6uoqfD4fhoaG4PP59ux4ZYfjuro69PX14fz582r1Rv5cuSqQTCY37SOylwo6GMhegpf52C0tLbDb7bvqPGg0GlFfX4/KykpVBrOxsRHd3d2wWCwoLy+H3+9HJpPB2NiYWhkopE2e+tdH/3FZWRlaW1sxMDAAm82GkZGRXf0cmULT2NiolmFlfwPZ1VCmTCwtLSGZTGJqagoejwcPHjxQddcLbUZyJ5aXl1UVAdmaXQihqifJmW55Hm02cJYrAX6/H6FQSC1lF1rahiQHAJFIRN3M9OeZ7EMgy8r29PTAZDLB6/UiEAio2VeDwaA2TFZUVKC9vX3DTVyfx7qysoJQKIRYLAa73Q6v1wu/37+ne4OeRiqVUnnJLpcLtbW1Kj0un+V8q6qq0NzcjOrqavT29qK/vx+NjY2oqal5pFKJ3Lyt/7cfr6FMI4tGo5sGyLKEMwC18gIAFosFHo8Hy8vLauZ5r1RVVaGzs1O9d7PvJbLCUSwWQzqdzvvrJvvYyFK1sseG3+9Xg+GJiYlHVuuAtb+77Ei81QBCBo0yNaijo0Pd9/TnkJwEkYHA9PQ0PB5PTsp5b5dsVjU+Po5EIoHDhw+rlDmLxYKSkhIcO3ZM9QFyOp37VrJSvs6VlZXqn74ZYSQSycu5pZ9INZlMiMfjG3LkS0tLVXaELOKgn9SjvSEnw+RE28DAABoaGtDR0YGqqioYDAY1WRgMBjE9PQ2/349gMJjz86SggwE9mSPb1dUFs9kMr9e76YVwu/TR6/nz53H69GlUVVWhoaEBFRUVavYjEAjg3XffxZUrV1QHxEKy2SYdg8GAs2fP4ujRo7h//z6uXbv21JsuZTOpsrIynD59Gp/97GdhNptx9OhRtLa2qlWV8vJyJBIJBINBeL1evP3227hz5w78fj+mp6dVnl2hDnyfJJVKYWJiAteuXVMdZUtLS3HkyBGcPXsWVqtVVbWQy4t6yWQS8/PzCAQCmJiYUEFFPrtn75QMgOx2O1ZXV9WMoWQwGGCxWLC6uopYLAa/349AIIDp6Wk4HA4YjUYcOXJEbageGBhQN0SDwQCDwaD6YwBr53I6ncbExARmZ2cxMTGB4eFhLCwsbNkgJx/k7Gw8HgcA3L59G4FAAAMDA2hvb89rqlBDQwOef/55NDU14ZVXXsGnP/1ptY8hewOs3OAmyzHv9YB6u2Q5PlkCMvsYDAYDGhsbsbKygk9/+tM4dOgQgLUN5bLTZzAY3NO/v9lsxksvvYSenh6cOHHikb+hnCmVZX/z/bqFQiHMzs4iGAzi4sWLuHLlyoZc62AwiKmpqU0DUVn6WqaIbhYMyL1AsoP9iy++qKob6aVSKYyPj8Nut+PBgwf44IMP4PV6EQ6H9y0YSKVSuHPnDhYWFjAwMID+/n4YDAaYTCacPHlSTRycOHECY2NjePvtt/ctGJDHZTab1b+lpSXMzs7CbrfD4XDkZaUumUxicnJSpeKcPHkSAFSXa4PBAKvViurqavh8Phw6dAi1tbVwuVzwer05P75iUVZWhvr6etTU1OD06dN488030draCqvVqjbqJ5NJlbZ18eJF1aV534MBIUQngH8BYAWwCuB/a5r2t0KIBgDfB9ADYA7Am5qm5WS9R1+HvLq6WnVtlBtgNmsYs9XOefm5yWRCc3OzmhXp6elRg7jS0lJ4PB6Ew2H4/X7VXKKQBrL6Wb/s45L7IOrq6rC4uIjq6upHbnbyQp5dxzo7sJDL2AaDAQ0NDejq6kJDQwMOHz6s8gH1NaflTLLT6cTs7CxisRhisdi+1w/fLTm7GQqFkEgkIIRQlZRaW1vVZnP90qj+XyaTQSQSQSAQUANLWe+5UOk3TW82OypTLQCoBjXl5eWIRCJIJBIqjUyeL8eOHduwEgA82thJ0zTEYjE1CxqNRtUgfD/JFIWlpSWEw2EEg8G8LfEDDyt0yLQWi8WilrQ3SwPRV37JdX3qJ5GbE+Xq4MrKyqb9JWTKSmlpKcLhMOrr61FdXY3V1VV1/drt7yF/ZkVFBZqammC1WlFXV/dIx3F9SeR8rQzo03KWl5fVpm+v1wun06mKZ+j77Oz0+8vXuqKiQm0aNpvNqK+vh8lkUvdTGUjKwEPuWfB4PGpldD/Pp3A4jOXlZVitVrUCIjeBypVxucooN0HL3ysf5N9S38dGvzk+Go0iEAggkUjk5XWUEzYGg0GtSOtTiOU+k5WVFZhMJtTV1SGTySAQCKhiIFv9nvK9I8vTFkK/gEIjXyd5PphMJjQ0NKi9OnI8u7KyovZbRiIRuN1uuFyuvJwn2/mrLQP4b5qm3RZC1AAYFEL8HMDvALioadq3hBDfBPBNAH+Wi4PUN6mQnVzb2tpw4cIFdHZ2qo1ycgZsZWVF5bTL3NrGxsYNNfHb29tx9OhRmEwmtLS0qGpBsnTbxMQE7t27h2AwmPcNi0+ysrKiNgVarVaYzWbVNTd7dtBsNuMzn/kMWlpaNvz/Pp9PVfeRgxqLxYKOjo4Nb+bS0lK1lHj69GkcOXJEVZmQuY/z8/OIx+Mq11DOfgeDQaRSqYLfH/C4wFGed7Ibb0dHh8r3LS8vx4ULF3Du3DnU1dWpyD4711DWJ/7oo48wNzenymUehGBApitsthFdf3Eym804duwYlpaW0N7ejkAggIqKCjXbJAMFfbCe3YAqkUggGo3C6XRiampKDYCKncz3ra+vR39/Py5cuKAmMbJnhuVKQCqVgt1ux9DQEJxO575uzgsEArhx4wZmZmaQSCRgMBjUrLQsQQw8XEI3m80Ih8P41Kc+hba2Nvj9fszMzGBpaQkOhwMul+upyvrJVDe5Ebyrq0ulCpWUlGy4LtpsNly5cgUulwuTk5M5f59mp3zK3HhZQUvel57mWio7gjc0NKgiDy0tLWhtbUVtba0qmiFXgFdXV3Hr1i1cu3ZNldGWJbSj0aiqILZfZPrK6uoq5ubm8NOf/hT37t3DwMAAzp07h5KSEjQ1Nak9Ey+++CK6u7s3lITMJTlJVFFRgd7eXpWB0N7eriY6RkdHMTY2BpfLlZd7gEwTSiaT8Pl8iMfjG6oZyTShiooK9PX14Vd+5VcQCARw9epVVQI0u4qO3BRttVpVBSHZN2M3WRvPEhkEtLe3o6OjQ62Wt7S04PDhw6p6F7BWvjQSiWBwcBAOhwNjY2OYm5tDMBhELBbb/2BA0zQXANf6x1EhxAMA7QDeAPDq+tP+GcAvsMfBQPYvr29Y1NXVhVdeeQV+vx/Dw8MYGRlRJ2wmk1GbdUwmE06cOIH+/n4YjUZYrVZVQqq+vh6lpaUqNzMej8NmsyEajWJ0dBS3b99WMzOFtCqwsrKCubk5CCFU9RqZ75ldH7ehoQGf+9zn8MILLwCAqiwxMTGB+fl5JJNJeL1epFIpnDx5Ei+//LIKuICH0X5paSlaW1vR2dmpbhjpdBoejwfXr1+H1+vF5OQkJicnkUgk4HK5VC3oQnrtHmer2siapqkqN93d3WrG22g04syZM/jEJz6h8i+z/z+Px4PJyUnMz8/jvffew8jIiIr+5XMKlSypJmfetso5l3snZAqIfJ6+8Yw8j7JXAbKDAdnQTzZl2+/29IVA3qAPHTqE5557Dq+88orqX7FZMCDLUNpsNtULZT8rnvn9fly5ckUNzmReskxv0pM9Lebn5/Haa6/hhRdegN1ux/3791VQIbu97/TaUlpaiubmZrS3t6O3txeHDx/G4cOHUV1djdLSUqTTabhcLtWs8aOPPoLNZlOz8bmy2TUnk8moCiIyIHiaNCkhBAwGAyoqKtDZ2YnXX38dXV1d6OnpUSmNMliQx5BOp3H16lV8+9vfRiwWUytLcqVpv69Z8rokUyl+/OMfw2g04stf/jJ6e3tRW1sLi8WCyspKVFVVIRgMwu124/Lly3tWhfBxZAO4uro6nDt3Dl/96ldhNpthMBjUbO/9+/dx9+5dtbqfa5lMBj6fD0IIuN3uR/r+lJWVoba2FpqmqfRPuadhYWFBjav074Py8nL09fXh1KlTaGhowPHjx1FfX6/6zRQ7/WpAV1cXXn75ZdTX1+PMmTPo6OhQ2Snl5eWq1LjH48GlS5dw9+5deDweTE9PqzL6OT/eHZaN6gFwCcBJADZN08y6rwU1TauXn4fDYfWNJycn9+JYiYiIiIhoh/r6+tTHdXV1Gyr0bDu5SwhhAvBDAH+iaVpkp8125Mz0drW0tOD06dNobGzEF7/4RXz5y19W0aascuLz+ZBMJuF2u+HxeNSMWCaTUU0cZEfUuro61XijrKxMzUKm02k4HA71Pa5evapmtWVupKzXXyiEEDAajaiurkZDQwMuXLiAjo4OHD9+HJ/5zGdUUxHZcEZG9XImTVa20eefyhm75ubmR1pxSzLvN5lMqqoSCwsLuHr1qqpyEYlENjRU2czNmzd3fD7kQklJCU6dOoXz58/DYrHg13/913H27NlNy+qlUinYbDYEg0GV9yc3A8kqO/pZWplv+8Mf/hDvvvuuWsE6SJux2tra8PWvfx2f/OQn0dbWhuPHj2/YYLjZCsp2JhfkTKP+3FxcXMTbb7+NYDCIwcFBTE9Pq9nRQkqlamlpwWuvvYauri48//zzeP311zfMbicSCSwuLiIej2NhYQEzMzMIhUK4du0aHjx4oJqAbTUbWFlZifr6ehgMBgwMDOD06dMwGo3o7u5WVT56enpQWVn5yH6dlZUVTExM4NKlS/D5fLh9+zaGh4dVfvl+9bSQaXb65XKLxYIvfelLGBgYUCtt+oZft2/fxrFjx1RVn2AwqPpZyH4lY2Njal+J3+/f8jwxm82wWCyora3F5z//eZw/fx5msxl9fX2qL0pZWRkCgQD+6Z/+CT/72c8QDAYxMzODaDS6ab+M3bJarTh69CgaGhrw1ltv4Y033tiQnik7+0ajUQwNDWFiYgI+nw83b96Ex+PZ9s8xmUx49dVXMTAwgLa2Nrz00ksqNaG2thYrKytwu92IRCIIhUKqCtH777+Pd999t2DS9La6Z8gCF6Wlparkb21tLS5cuIAjR46gurpaVXfzeDxYXFzE0tIS5ufnEQwGEQwGYbPZkEgkNqSgPY4cRxgMBtTU1Kj0mIGBARw9ehRVVVWwWCyquEZzczNSqRQ+/PBDDA4Owu12q9X0fJGbq61WK86cOYOvfOUraGpqQk1NjUpT0e+XkY3u5ubm4HA4Hrmny308FotFXa/k6odckSkpKVHf0+Vy4fr162qMdfHiRdVN92n3E+7lOEKuaDQ0NOA3f/M38bnPfU4VlCkvL8fk5CTu3bunSnzKVcTr16+rMvfZaXwWiwXHjx9HXV0dTp48qYrU1NfXw2g0qntgJpPB9evX1arn9evXYbfb1R6hvUxre1wBnG0FA0KIcqwFAv+madq/rz/sFkK0aprmEkK0Atj+FWqX5GDYaDRC0zS0t7erZVSZh20ymdDU1LRlpQ+5GTQej8PhcKid/YODg09deSdf5CbLWCyGeDwOg8GAubk5lJSU4Pnnn0dlZaUqlVdeXr6hfrZ8c1qt1i2//2aDvNXVVTgcDlVd48GDB5idnYXT6VTNjA6q7NSVzRgMBhw5cmTb3wdYu2C6XC4MDw8jHo8XdOWgzejLzck9KXtBDq7S6bTajwAAH3/8MYLBoOpJcRBVVFSgpaUFy8vLaG5uRn9/vyoN53K5VK7zVjdAfS79Jz7xCbz++uuq9rQsSSe7DMvcbflaLi8vw+Vy4caNG3A6nbDZbHC5XAWR1iF/X5vNBqfTiZaWFpw4cQJNTU1YWVmBxWJ5pPuvDDxlkzo5EKmsrFR5z/L7P672t2zu09TUpFIhZaMv/TmdyWQwMzODGzduqImlfOXGZwfSlZWVaGxsRHV1NVKpFEwmE2w2G8bGxnYUDMimYufOnUNTUxN6enpQW1urvp7JZBAOh+HxeODxeDA0NIRgMKgGgIVO7msCgLGxMYyNjcFsNqO0tBTLy8uq90tjY6PqYxKLxXD//n1VU7+0tFSVM15cXHziz5SbbeX+lvb2dpjNZly4cEFV95JllOUEWTQaxdjYGD744AMkk0nVfDNfKioq0N3djSNHjuDw4cMqrVPfsVyef2VlZSgrK0NVVRVOnTqFU6dOAXigt9NVAAAdIElEQVQ4LnjSpM9mE4mxWAyjo6OYnZ3F5OQk/H5/QTXcjMViKkVcdkXOZDKora2FwWBQm3z1abDy91hYWEAqlXokcJYTRs3Nzejr68OxY8c2BPyyOEcikcD09DSuXLmi9mDuR++F7VQTEgD+D4AHmqb9je5L7wD4OoBvrf/37b08MLmBSlaWiEQiqKqqUjdE3fGpijfywi1blMsLvcx3lJ1f0+k0FhYWMDo6qi4AXq8XPp+voPoIbMfKyoraHDg/P4/BwUHU19ejsbERZrNZ3VSy24I/jjzh9bO36XQac3NzGB4eRjQaxdzcHFwuFwKBwIGuFCTPC33VFXlO7fT7yO+xtLSk+gjIUoAHYSN1tlQqhfn5edXQp6amRt1AdnI+AQ+7o8oqKYlEQm06j8ViOHXqlFpdKqSbRDa5eVGuFqZSKVWBTO6PkOeO7K+wvLyMvr4+NaMUj8e3fM/Imcyqqip0d3erLrlypU/2YtA0DX6/H36/H6lUSjX0Gx0dhdvtVk2NCo0MXmS+t6z3L5vvyMorm5Hdq2VJ45MnT6K5uRmRSASHDh3a8jWVmzfr6urU3ioAarVUrnB6PB5V6jHXA+Hsa+zy8vKGPTZlZWUwGo1qr5IsDnH27FnVoXs76urq0Nvbi6amJpUTLmd95SrVzMwM5ubmEAgEMDc3h2g0ilAotO9B5NPKZDJwuVwwmUxIJpNoampCQ0MDamtrUVtbi0wmA6PRqCpxybGGPE+edJ0uKytT5TjlOWU0GtHW1qaui16vF6urq/D7/XA4HKrBlrwG5PtekMlk4Ha7VSWpnp4e1XhMrnJvVnVwr2xWNa6Q6Esw66uvyeMsLS1FZWXlhscbGhrw3HPPobGxUTUH0/9eFotFBYomk2nDeFROYExOTiIcDqusg8fdG3JtOysDLwP4LQBDQoi764/9OdaCgB8IIX4PgA3AV/fywJaWluB0OhEOhzE/Pw+bzYba2lpYrdZHGtfITcXAwxdablgE1gY16XRaReey5fYHH3ygWjzLi+N+Vt14Gul0WgUz4XAYU1NTMJlMOH78OPr6+mCxWPDiiy8+diUgm9ygJV8zu92OWCyGy5cv4+LFixtKbMoqTAeRXBKVdX3leaIf2G2XLHcnu+bKDY/Dw8MqYDpoQVMkEsHFixcxODiIM2fOqI1lVqsVLS0tO3p9QqGQqiZjs9mwuLgIn8+HwcFB+P1+/Ou//iumpqbULHehktcIn8+nNuYKIVBdXa1mmuW1qKKiQtUY/9KXvoTPfvaz6pzbitxwJhu2yZuITN2TG/eXl5cxOjqqNgjfuXMHDodDpdRsdnMqBPob4cWLF3H9+nUMDAxACAGr1aqq/GymtLQUHR0daGpqQiaTwalTp1SJRBlsbzbIkp2G5cbOyspKNQBMJpP4xS9+ge9+97sIBAKw2+1IJpM5L0MpJyHS6fSGngAysJSpJpqmwWKxqIH78ePHd7QZ3GAwoKenBy0tLSoAkRuTQ6EQ/H4/3nvvPdy6dQtLS0vqWiULBhxES0tLGBwcxOjoKNra2jA3N6c2uJ48eRKVlZWqSWImk8ELL7ygAtTtpEXJgF8GbfL9KlNu5HXf4/FgamoKt27dUlWZZPOofL+28Xgcd+7cwejoKBwOB1ZXV9HQ0IDm5mY0NDSgpqYGR44cURXxHmcnG/ezi0QUKpkGJccA2RMCstKW/veorq5GS0uLus5m/476UrdyclrTNPX9Z2Zm8KMf/QgulwsTExOYmZnZ13HCdqoJXQaw1V3/s3t7OA/JXFFgbQknGo2ipKRERdXZlV/0y73yj6KvtS2XZAKBALxeLxwOB6anp/O+XLfX9Muk8uOqqipVy1Z2BZYn9uMGcPo9BbIFdiKRQCgUQjgcxuLiIubn59Ws6EEb3G5GzgjIN6EcqMlqODv5PrJLrZwFkukh6XR637q/7kYmk4HH40EoFEJLSwtCoRAqKys39FLY7DXKHkjJ4FK/FO9wOFQzFZn2cBDSqOT+EX3wmJ0rqq/EBTxcIdiOrV5PSZ8WJM8zj8eD8fFx2Gy23fxqeSNviG63G16vF/X19fD5fKr2vz4vGdiYdiDTMzRNg9lsVq+H/pqf/T7Tl1yWexfkMaRSKdVgUl8KOdfvVXmdlTnD6XR6w/HJVBQAqjtzTU0NDAaDut5v5xjLysrUYC+dTqt9colEQvXLcDqdmJube+J+loNidXVVdfhdXV1FY2MjIpGI6vQq8/yzu5/r/7td+vNNfpzJZFQFo4WFBUxPT6sAbL/umXKyCljLj3e5XEilUurY5f1L/57Ty04lyv5Y/7zs65V+b1Mhk9cSOR6QAYGcXM5Ok5WD/e2Qr4GczMlkMmqPjmxAms++NZsp2O4QMhiQM2DV1dWora3F8ePHVbe2zs7OLW+yoVAIDocDS0tLsNvtcLlciEajmJycVHnJhbI5aq/oN2ROT08jHo/DarWivLwc7e3taGpqUhv15E11ZWVFbST2eDxwOp1Ip9MqzUB2woxEIpiamlIpQwf9hiEFg0FMTEzA6/Wivb0dyWQSzc3NeO6557b9RgfWZtE//PBDTExMwO/3Y25uTpUK2++mT09LXrhWVlZgt9tx+fJlNDY24ty5c6pBn8yplOSMh81mUwOdlZUVuFwu1Yna5XKp7ssy4D8oUqmUqtvf1taGxcVFZDIZtcyea8FgEGNjYwiHw7h9+zbu37+vcpIPGvme8Hg8+Pjjj2E2mzE5OYnOzk709/fjypUrqK6uRk1NjSrBJxsgZpet1Q/Isukni+TNPhQKYXh4eEOHdBnU5eN9urS0BI/Hg1gshrt378JkMqnUFaPRqNJPsnu+yD452z1G/UpVNBpV+7tsNpuaDJudnUUymVTXqWdJIpHA/Pw8PB6PWpWsqqpSzeaMRiOam5tVp2A5gbYdmqYhGo3C5/OpjtUyPXR6ehqhUAiLi4sqyCyUybNgMIihoSEYjUY0NDTAbDajqakJmqaho6MDdXV1aG5u3nXzMBmUxeNxuFwuOJ3ODSXHC1Emk8Ho6KhaPUqlUhs23Musge3un5MBosPhUCVaFxYWEAgEMDs7i6mpKXXO7PdrUrDBgGyGlUgkcO/ePczPz6O2thYOhwM9PT3o7e1Fc3PzljfgQCCAe/fuwe/34+7duxgZGVGDMzn7UShvzr0io1khBKLRKCYmJtDU1IRoNAqr1Yrjx4/jpZdeQk1NjdpHIDeQLS0tYWRkBNeuXUM8HlcpEPJmIau67Gfnyb2maRq8Xi8CgQBqampQXV0Nt9uNo0ePoqura0fBQCgUwn/8x3/g5z//uRpw6GcaDiKZHy+EwOzsLOLxuFr2lHmyFRUVG4KBTCaDkZERXLp0SVWxSaVS8Pv9qsa3DAKelDJTiGT9ftn40G63I5PJqPdUrnm9Xly/fh0ulwu3bt3C7du3VcB10MiBp9PpxEcffaTy45uamvDXf/3X+OCDD2CxWNDW1oaTJ0+qG7I83/T7ejabucwmVwMymQz8fj8GBwcxPz+PBw8eqCZM+SJXV8vLy3Hjxg1EIhHU1tbi2LFjaGhoQE9PD+rr6x8JBmSK7HZTL/SBUDgcxt27d7GwsIDh4WEMDg5uSIt4Vq7rerFYDNPT0xBCYHx8XDWLlIFAa2srTpw4gdraWrS2tsJqte6oUILceyhTkKenp1UKmH71p5BSZWTQIqsiVVdXo7W1FQaDAaFQCJ2dnTCbzbsOBlZWVlSzOpvNBpvNhoWFBbViU4hSqRTu3r0Lm82mqrbJNDuDwaAmJHYSDCwvL2N2dhaXL1+G3+/HzZs3MT8/r7oMy0IQ+31+FGwwADy84CWTSTXz5ff7UVVVBaPRqDbkAI8uT7lcLrjdbgQCAZWukUwmVZ7os0q+ZvKCFI/HEQwGVWk1p9OJ2traDQM1uXHa7Xar5SrZ9S4ajR7IGdztkkt3ssKDz+dDc3MzwuGw2lAlBx/ytZUpVPoBmLzAHrQ9J9shl75lF0S/3w+Xy4Wamho1QJGvjf48kmlTcqVJfiz3oxxE8uIuhFApdDJ1Q5++uBdk+og851ZWVuD3+1VX61gspsoCH2Ry4kfmX8vXTzZ7lNcuk8mkJjxk6k9ZWZkqDZpdAU3Sb9KVDSbl5utAIIBEIpH3wYk+bzwWi6nSqH6/H0IINDc3P3JMOynlq99HIe8HgUAAwWBQXaei0egzNyGWTT/hIH9XGUSm02mUl5era1V2et92eDwe+Hw+tUJwEDZeyxVf/eeRSESlUNfU1CCRSKj32FbFNPTBtRxv6K+Bsimp3F91ENKL5XgzEokgGAzC4/GovWBy8stoNKrSqfK6o9+jpR/cy8/9fj+8Xq+67oRCoYLb11XQwYAk9wmkUincvn0bExMTuHnzJi5fvrxlp7tIJAKPx4NkMqmCgWdxNeBJEokEJiYmsLCwgNnZWdy+fVuleBgMBrWBVkbxHo9nw4bqdDp9YAduO5FOp1XvhGAwiPr6elitVhw5cgT9/f1q70Umk8HCwoKazZPkPpRnlcw3TiQS+OijjzA1NYXy8nJV8QR4eONdWFjA4uKiuunI964s/3vQB6/A2u9qt9vx4YcfwmKxoKamBh0dHY8MaHfz/d1uN+x2O+LxOKamptQ5dvfuXYTDYfj9/oKdYdsJGYyXlJSoHigA8Mtf/hJVVVVobGzE4OAgqqur0dbWhpaWFphMJnR1daGuru6JaQ3hcBhut3vDtdDtduPjjz+G1+tVs3P7YXV1FW63G0tLS6iurobf70dNTQ1WVlZw5swZVQN+J+S5OTo6qlI94/E4vF4v7t+/r/LZn4Vz52msrq4ikUioVZFYLKY2mVdXV+/ovSsnzmQOeKEM7LZLDuC9Xi8+/vhjjI+P49SpU6ioqFD3wJaWlk3/33Q6jYmJCczPzyMej6trFQA1QNb3AZmdnVUFWwr13JPnhhz/vPvuuzCZTKqHRU1NDXp7e9HR0QGj0Yj29nZUVFRgcXERc3NzSCaT8Hg8Ku1HZgfcu3cPd+7cQSKRgNfrzdvepJ04EMGAnNWREZv0pM2w9DDHOVv2a1fsr5es0e5yubC8vIzOzk54PB4YjUYcOXJEbV5PJpNwOBy4evXqhnrfcnbvWaXf0D88PIzh4WEAT97w+izz+/0YGhpCU1MTzp8/v2HmcaelabNpmoZQKIS5uTn4/X5cu3YNs7OzGyreFNrN5GnpqwDpU3VGR0cBrG2cnZiYQGVlJQ4dOoSuri7U19dD0x72S3lcipZsBBcOh3Hv3j2MjIzA7/djbGxs3/ujyL9zKBRCRUUFotEojEYjenp6djUJ4/P5MDQ0hGg0qooZxONx1XRLNlMrRrIIAADVvLRYybHV8vIyxsfH1X6cI0eOIJFIwGQywWKxbHqdl72H5N6b4eHhDfXxV1dX1eplKpVSQVOhk8cri10IIdDd3Q2n0wmz2awmDhobG9V+k2AwiOnpaUSjUczMzKixgQwGZA+MQk6LPRDBwFaehRvhfuFrt7WlpSU4HA4kk0k1MyfTQlKpFKanp7G4uIhAIKD+n2QyWRQrKNmK+TxaXl7G0tISIpEIxsfHcfnyZZSVle0opzSbvOmurq5ifn5edcF1u90Ih8OqqMKzEghsh5wIkqkuZWVliEQiMBgMWFxcRH19Pex2O8rLyzdNnfH5fKo8ss1mU2kdhTYwkeUNhRBwOBy4du0aGhoaNjxnq1QoPU3TMD4+jrm5OcTjcfh8PkQiEVUdTu4RIJL06VSBQADj4+OqM7Xdbt9QuVG+x1KpFIaHh1WBEZm6qP+e8r54kPcayvTXQCCgMgjS6TRqamrg9/tRXV2N2dlZTExMqAIZ+n0RcrWh0IPvAx0MEOWC3+/Hxx9/jPLycly+fFnVepcpLslkUrUgl/QlXqk4JJNJlSf8/e9/H++//z4A7LhHxVaWlpZUapXc4CrLihb6jWUvpVIpBAIBlJSUIBQKYXp6GmVlZbh+/bpKy9J3Ut3s/5d7LuSeMRlgFBJZ/jEajeI///M/MTY29kiDTb3HDa4SiQQikciGMon6vO5iCibpyVZXV7G0tISSkhKMjY1hcXERZWVlqKysVKnY2Zv05cy/nKBYWlp6JMjUn2sHOT00GAxiaWkJpaWlmJqaUs1v5X6CZDKpSoOmUqkNG/Jl4FTo12wGA0RZZJk4oseRwWEqlVK5srT39BseC20Av5fkRn0AcLvdcLvd+3xEVEzkwD0cDh/4/kt7Tb93Up8K9Sx5urVsIiIiIiI68BgMEBEREREVKQYDRERERERFisEAEREREVGRYjBARERERFSkRK7Ki4XDYdYtIyIiIiIqIHV1dRtqFXNlgIiIiIioSDEYICIiIiIqUjlLEyIiIiIiosLGlQEiIiIioiLFYICIiIiIqEjlNBgQQnxeCDEuhJgSQnwzlz+LCpcQYk4IMSSEuCuEuLX+WIMQ4udCiMn1/9bv93FSbgghviOE8AghhnWPbfn3F0L89/VrxrgQ4vX9OWrKlS3Oh78SQjjWrxF3hRBf0H2N58MzTAjRKYT4SAjxQAgxIoT4r+uP8xpRpB5zTvA6kSM52zMghCgFMAHgVwHYAdwE8DVN00Zz8gOpYAkh5gB8UtM0n+6x/wEgoGnat9YDxXpN0/5sv46RckcI8QqAGIB/0TTt5Ppjm/79hRDHAXwXwDkAbQA+ANCvadrKPh0+7bEtzoe/AhDTNO1/Zj2X58MzTgjRCqBV07TbQogaAIMAvgzgd8BrRFF6zDnxJnidyIlcrgycAzCladqMpmlpAN8D8EYOfx4dLG8A+Of1j/8Za290egZpmnYJQCDr4a3+/m8A+J6maSlN02YBTGHtWkLPiC3Oh63wfHjGaZrm0jTt9vrHUQAPALSD14ii9ZhzYis8J3Ypl8FAO4AF3ed2PP6PSc8uDcDPhBCDQojfX3+sRdM0F7D2xgdg2bejo/2w1d+f143i9Q0hxP31NCKZEsLzoYgIIXoAnAFwHbxGEB45JwBeJ3Iil8GA2OQx1jEtTi9rmnYWwK8B+KP1NAGizfC6UZz+HkAvgE8AcAH49vrjPB+KhBDCBOCHAP5E07TI4566yWM8J55Bm5wTvE7kSC6DATuATt3nHQCcOfx5VKA0TXOu/9cD4EdYW75zr+cFyvxAz/4dIe2Drf7+vG4UIU3T3JqmrWiatgrgH/FwiZ/nQxEQQpRjbdD3b5qm/fv6w7xGFLHNzgleJ3Inl8HATQB9QohDQggDgLcAvJPDn0cFSAhhXN8ABCGEEcDnAAxj7Vz4+vrTvg7g7f05QtonW/393wHwlhCiQghxCEAfgBv7cHyUR3LQt+4rWLtGADwfnnlCCAHg/wB4oGna3+i+xGtEkdrqnOB1InfKcvWNNU1bFkJ8A8D7AEoBfEfTtJFc/TwqWC0AfrT23kYZgP+nadpPhRA3AfxACPF7AGwAvrqPx0g5JIT4LoBXATQJIewA/hLAt7DJ31/TtBEhxA8AjAJYBvBHB6UihBDiOwC+CMAjq+RkfV0A+FsAXwCQAPA7cpNcMdnifHhVCPEJrC3tzwH4L8DBPh9o214G8FsAhoQQd9cf+3M8g9cI2ratzomv8TqRGzkrLUpEVEw2K5mZ9fUvAPhjrAUD5wH8raZp5/N7lERERBvlbGWAiKiYaJp2ab3yxVbewFqgoAG4JoQwCyFaZcUUAAiHw5yd2aG6urrNNg8SEdE25bQDMRERKSx/R0REBYfBABFRfrD8HRERFRwGA0RE+XFgy99NTk7u9yEQEVGOMBggIsqPdwD8tlhzAUBYv1+AiIhoP3ADMRHRHtiiZGY5AGia9g8AfoK1SkJTWCst+rv7c6REREQPMRggItoDmqZ97Qlf1wD8UZ4Oh4iIaFuYJkREREREVKQYDBARERERFSkGA0RERERERYrBABERERFRkWIwQERERERUpBgMEBEREREVKQYDRERERERFisEAEREREVGRYjBARERERFSkGAwQERERERUpBgNEREREREWKwQARERERUZFiMEBEREREVKQYDBARERERFSkGA0RERERERYrBABERERFRkWIwQERERERUpBgMEBEREREVKQYDRERERERFisEAEREREVGRYjBARERERFSkGAwQERERERUpBgNEREREREWKwQARERERUZFiMEBEREREVKQYDBARERERFSkGA0RERERERYrBABERERFRkWIwQES0B4QQnxdCjAshpoQQ39zk668KIcJCiLvr//5iP46TiIhIr2y/D4CI6KATQpQC+DsAvwrADuCmEOIdTdNGs576S03Tvpj3AyQiItoCVwaIiHbvHIApTdNmNE1LA/gegDf2+ZiIiIieiCsDRES71w5gQfe5HcD5TZ73ohDiHgAngD/VNG1kq284OTm5t0e4S4V0PH19fft9CEREzwwGA0REuyc2eUzL+vw2gG5N02JCiC8A+DGALUe1hTTgnZycLKjjISKivcM0ISKi3bMD6NR93oG12X9F07SIpmmx9Y9/AqBcCNGUv0MkIiJ6FIMBIqLduwmgTwhxSAhhAPAWgHf0TxBCWIUQYv3jc1i7/vrzfqREREQ6TBMiItolTdOWhRDfAPA+gFIA39E0bUQI8QfrX/8HAL8B4A+FEMsAlgC8pWladioRERFRXgnei4iICkM4HC7IC3Ih7xmoq6vbbL8GERFtE9OEiIiIiIiKFIMBIiIiIqIixWCAiIiIiKhIMRggIiIiIipSDAaIiIiIiIoUgwEiIiIioiLFYICIiIiIqEgxGCAiIiIiKlIMBoiIiIiIihSDASIiIiKiIsVggIiIiIioSDEYICIiIiIqUgwGiIiIiIiKFIMBIiIiIqIixWCAiIiIiKhIMRggIiIiIipSDAaIiIiIiIoUgwEiIiIioiLFYICIiIiIqEgxGCAiIiIiKlIMBoiIiIiIihSDASIiIiKiIsVggIiIiIioSDEYICIiIiIqUgwGiIiIiIiKFIMBIiIiIqIixWCAiIiIiKhIMRggIiIiIipSDAaIiPaAEOLzQohxIcSUEOKbm3xdCCH+1/rX7wshzu7HcRIREekxGCAi2iUhRCmAvwPwawCOA/iaEOJ41tN+DUDf+r/fB/D3eT1IIiKiTTAYICLavXMApjRNm9E0LQ3gewDeyHrOGwD+RVtzDYBZCNGa7wN9Gn19fft9CERElCMMBoiIdq8dwILuc/v6Yzt9DhERUV4xGCAi2j2xyWPaUzyHiIgor8r2+wCIiJ4BdgCdus87ADh3+py6urrNAgYiIqKc4coAEdHu3QTQJ4Q4JIQwAHgLwDtZz3kHwG+vVxW6ACCsaZor3wdKRESkx2CAiGiXNE1bBvANAO8DeADgB5qmjQgh/kAI8QfrT/sJgBkAUwD+EcD/LbRSpNsoj/qqECIshLi7/u8vcnw83xFCeIQQw1t8neVaiYh2SWgaU1aJiPJpvRTpBIBfxVr60E0AX9M0bVT3nC8A+GMAXwBwHsDfapp2fp+P6VUAf6pp2hdzdRxZx/QKgBjWqjCd3OTreX2NiIieRVwZICLKv0IsRbqdY8orTdMuAQg85ikHtlwrEVGhYDBARJR/hViKdLs/70UhxD0hxHtCiBM5PJ7tYLlWIqJdYjUhIqL8K8RSpNv5ebcBdGuaFltP0fkx1joq7xeWayUi2iWuDBAR5d+elCLN9zFpmhbRNC22/vFPAJQLIZpyeExPku/XiIjomcNggIgo/wqxFOkTj0kIYRVCiPWPz2HtHuLP4TE9Ccu1EhHtEtOEiIjyTNO0ZSGELEVaCuA7shTp+tf/AWulSL+AtVKkCQC/WwDH9BsA/lAIsQxgCcBbWg5L0gkhvgvgVQBNQgg7gL8EUK47nry+RkREzyKWFiUiIiIiKlJMEyIiIiIiKlIMBoiIiIiIihSDASIiIiKiIsVggIiIiIioSDEYICIiIiIqUgwGiIiIiIiKFIMBIiIiIqIixWCAiIiIiKhI/X+iM2S6gq2OSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 2160x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "reconst, _, _, _ = PCA(X, num_components=10)\n",
    "num_images_to_show = 10\n",
    "reconst_images = np.reshape(reconst[:num_images_to_show], (-1, 28, 28))\n",
    "fig, ax = plt.subplots(2, 1, figsize=(num_images_to_show * 3, 3))\n",
    "ax[0].imshow(np.concatenate(np.reshape(X[:num_images_to_show], (-1, 28, 28)), -1), cmap=\"gray\")\n",
    "ax[1].imshow(np.concatenate(reconst_images, -1), cmap=\"gray\");"
   ]
  },
  {
   "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 $10.0$ for our dataset?\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have provided a function in the next cell which computes the mean squared error (MSE), which will be useful for answering the question above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mse(predict, actual):\n",
    "    \"\"\"Helper function for computing the mean squared error (MSE)\"\"\"\n",
    "    return np.square(predict - actual).sum(axis=1).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "loss = []\n",
    "reconstructions = []\n",
    "# iterate over different number of principal components, and compute the MSE\n",
    "for num_component in range(1, 100, 5):\n",
    "    reconst, _, _, _ = PCA(X, num_component)\n",
    "    error = mse(reconst, X)\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\n",
    "loss = np.asarray(loss)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "# create a table showing the number of principal components and MSE\n",
    "pd.DataFrame(loss, columns=['no. of components', 'mse']).head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also put these numbers into perspective by plotting them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(loss[:,0], loss[:,1]);\n",
    "ax.axhline(10, 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 dont't tell us everything_! Just what does it mean _qualitatively_ for the loss to decrease from around\n",
    "$45.0$ to less than $10.0$?\n",
    "\n",
    "Let's find out! In the next cell, we draw the the leftmost image is the original dight. Then we show the reconstruction of the image on the right, in descending number of principal components used."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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",
    "    # concatenate the actual and reconstructed images as large image before plotting it\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 through the reconstructions for other digits. Once again, `interact` becomes handy for visualing the reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@interact(i=(0, 10))\n",
    "def show_pca_digits(i=1):\n",
    "    \"\"\"Show the i th digit and its reconstruction\"\"\"\n",
    "    plt.figure(figsize=(4,4))\n",
    "    actual_sample = X[i].reshape(28,28)\n",
    "    reconst_sample = (reconst[i, :]).reshape(28, 28)\n",
    "    plt.imshow(np.hstack([actual_sample, reconst_sample]), cmap='gray')\n",
    "    plt.grid(False)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA for high-dimensional datasets\n",
    "\n",
    "Sometimes, the dimensionality of our dataset may be larger than the number of samples we\n",
    "have. Then it might be inefficient to perform PCA with your implementation above. Instead,\n",
    "as mentioned in the lectures, you can implement PCA in a more efficient manner, which we\n",
    "call \"PCA for high dimensional data\" (PCA_high_dim).\n",
    "\n",
    "Below are the steps for performing PCA for high dimensional dataset\n",
    "1. Normalize the dataset matrix $X$ to obtain $\\overline{X}$ that has zero mean.\n",
    "2. Compute the matrix $\\overline{X}\\overline{X}^T$ (a $N$ by $N$ matrix with $N << D$)\n",
    "3. Compute eigenvalues $\\lambda$s and eigenvectors $V$ for $\\overline{X}\\overline{X}^T$ with shape (N, N). Compare this with computing the eigenspectrum of $\\overline{X}^T\\overline{X}$ which has shape (D, D), when $N << D$, computation of the eigenspectrum of $\\overline{X}\\overline{X}^T$ will be computationally less expensive.\n",
    "4. Compute the eigenvectors for the original covariance matrix as $\\overline{X}^TV$. Choose the eigenvectors associated with the `n` largest eigenvalues to be the basis of the principal subspace $U$.\n",
    "    1. Notice that $\\overline{X}^TV$ would give a matrix of shape (D, N) but the eigenvectors beyond the Dth column will have eigenvalues of 0, so it is safe to drop any columns beyond the D'th dimension. \n",
    "    2. Also note that the columns of $U$ will not be unit-length if we pre-multiply $V$ with $\\overline{X}^T$, so we will have to normalize the columns of $U$ so that they have unit-length to be consistent with the `PCA` implementation above.\n",
    "5. Compute the orthogonal projection of the data onto the subspace spanned by columns of $U$. \n",
    "\n",
    "Functions you wrote for earlier assignments will be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-6b8c115e4c40bd67",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "def PCA_high_dim(X, num_components):\n",
    "    \"\"\"Compute PCA for small sample size but high-dimensional features. \n",
    "    Args:\n",
    "        X: ndarray of size (N, D), where D is the dimension of the sample,\n",
    "           and N is the number of samples\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` pricipal components.\n",
    "    \"\"\"\n",
    "    # YOUR CODE HERE\n",
    "    # Uncomment and modify the code below\n",
    "    N, D = X.shape\n",
    "    # Normalize the dataset\n",
    "    X_normalized, mean = normalize(X)\n",
    "    # Find the covariance matrix\n",
    "    M = np.dot(X_normalized, X_normalized.T) / N\n",
    "    S = np.cov(X_normalized.T)\n",
    "    # Next find eigenvalues and corresponding eigenvectors for S\n",
    "    # Make sure that you only take the first D eigenvalues/vectors\n",
    "    # You can also take a look at the eigenvalues beyond column (D-1) and they should be \n",
    "    # zero (or a very small number due to finite floating point precision)\n",
    "    eig_vals, eig_vecs = eig(S)\n",
    "    eig_vals = eig_vals[0:D]\n",
    "    eig_vecs = eig_vecs[:, 0:D]\n",
    "    # Compute the eigenvalues and eigenvectors for the original system\n",
    "    \n",
    "    # Normalize the eigenvectors to have unit-length\n",
    "    # Take the top `num_components` of the eigenvalues / eigenvectors\n",
    "    # as the principal values and principal components\n",
    "    principal_values = eig_vals[0:num_components]\n",
    "    principal_components = eig_vecs[:, 0:num_components]\n",
    "    # reconstruct the images from the lower dimensional representation\n",
    "    # Remember to add back the sample mean\n",
    "    P = projection_matrix(principal_components)\n",
    "    reconst = (P @ X_normalized.T).T + mean\n",
    "    return reconst, mean, principal_values, principal_components"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "nbgrader": {
     "grade": false,
     "grade_id": "cell-ddfc3a4390a8b957",
     "locked": false,
     "schema_version": 3,
     "solution": true,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# YOUR CODE HERE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-8334eae5a3ead42f",
     "locked": true,
     "points": 2.5,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "# Some hidden tests below\n",
    "### ..."
   ]
  },
  {
   "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`. Given the same dataset, `PCA` and `PCA_high_dim` should give identical results.\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": 23,
   "metadata": {
    "nbgrader": {
     "grade": true,
     "grade_id": "cell-e1f35e6bd20a7b6a",
     "locked": true,
     "points": 1,
     "schema_version": 3,
     "solution": false,
     "task": false
    }
   },
   "outputs": [],
   "source": [
    "random = np.random.RandomState(0)\n",
    "# Generate some random data\n",
    "X = random.randn(5, 4)\n",
    "pca_rec, pca_mean, pca_pvs, pca_pcs = PCA(X, 2)\n",
    "pca_hd_rec, pca_hd_mean, pca_hd_pvs, pca_hd_pcs = PCA_high_dim(X, 2)\n",
    "# Check that the results returned by PCA and PCA_high_dim are identical\n",
    "np.testing.assert_allclose(pca_rec, pca_hd_rec)\n",
    "np.testing.assert_allclose(pca_mean, pca_hd_mean)\n",
    "np.testing.assert_allclose(pca_pvs, pca_pvs)\n",
    "np.testing.assert_allclose(pca_pcs, pca_pcs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Congratulations_! You have now learned how PCA works!"
   ]
  }
 ],
 "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.7.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {
    "height": "calc(100% - 180px)",
    "left": "10px",
    "top": "150px",
    "width": "169px"
   },
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}