notebooks/08-pca-mnist-dataset.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dimension reduction using PCA\n",
    "\n",
    "We need to calculate projections $\\tilde{x_n}$ of original data points $x_n$ into a lower dimensional space. Our dataset $\\chi = \\{x_1, x_2, ..., x_N\\}, x_n \\in \\rm{R}^D$\n",
    "\n",
    "$\\bar{x} = \\frac{1}{N} \\sum_{n=1}^N x_n$ \n",
    "\n",
    "$cov = \\frac{1}{N} \\sum_{n=1}^N (x_n - \\bar{x})(x_n - \\bar{x})^T$ \n",
    "\n",
    "We need to find a low-dimensional representation: \n",
    "\n",
    "$z_n = B^T x_n \\in \\rm{R}^M$ where projection matrix $B = [b_1, b_2, ..., b_M] \\in \\rm{R}^{D\\times M}$\n",
    "\n",
    "We assume columns of B are orthonormal: $b_i^Tb_j = 0$ if $i\\neq j$ and $b_i^Tbi = 1$ \n",
    "\n",
    "We seek an M-dimensional subspace $U \\subseteq \\rm{R}^D, dim(U) = M < D$\n",
    "\n",
    "We have to find a $B$ for which the $|x_n-\\tilde{x}_n|^2$ loss is minimum. \n",
    "\n",
    "$z = B^T x$ and $\\tilde{x} = B z$\n",
    "\n",
    "Dataset is available at <http://yann.lecun.com/exdb/mnist/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load necessary packages  \n",
    "import io, gzip \n",
    "from urllib.request import urlopen\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt \n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load images\n",
    "def mnist_images(url):  \n",
    "    with gzip.open(io.BytesIO(urlopen(url).read()), 'rb') as f:\n",
    "        data = np.frombuffer(f.read(), np.uint8, offset=16)\n",
    "        # consult http://yann.lecun.com/exdb/mnist/ for details \n",
    "    data = data.reshape(-1,784)\n",
    "    return data/256  # we are using int8, squared distance could be large\n",
    "\n",
    "# load image labels\n",
    "def mnist_labels(url):        \n",
    "    with gzip.open(io.BytesIO(urlopen(url).read()), 'rb') as f:\n",
    "        data = np.frombuffer(f.read(), np.uint8, offset=8)\n",
    "    return data\n",
    "\n",
    "# plot/view a specific image\n",
    "def mnist_view(im):\n",
    "    plt.axis('off')\n",
    "    plt.imshow(im.reshape((28,28)), cmap=plt.cm.gray, origin='upper')\n",
    "    plt.show()\n",
    "    return\n",
    "\n",
    "# calculate squared distance\n",
    "def dist(x1, x2):\n",
    "    return np.sum(np.square(x1-x2))\n",
    "\n",
    "# find nearest neighbor \n",
    "def find_nn(x):\n",
    "    dist_vec = [dist(x, train_data[i,]) for i in range(len(train_labels))]\n",
    "    print(\"Index: \", np.argmin(dist_vec))\n",
    "    return train_labels[np.argmin(dist_vec)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data = mnist_images('https://pranabdas.github.io/drive/datasets/mnist-handwritten-digits/train-images-idx3-ubyte.gz')\n",
    "test_data    = mnist_images('https://pranabdas.github.io/drive/datasets/mnist-handwritten-digits/t10k-images-idx3-ubyte.gz')\n",
    "train_labels = mnist_labels('https://pranabdas.github.io/drive/datasets/mnist-handwritten-digits/train-labels-idx1-ubyte.gz')\n",
    "test_labels  = mnist_labels('https://pranabdas.github.io/drive/datasets/mnist-handwritten-digits/t10k-labels-idx1-ubyte.gz')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute covariance matrix\n",
    "Sigma = np.cov(train_data, rowvar=0, bias=1)\n",
    "\n",
    "# Compute coordinate-wise variances, in increasing order\n",
    "coordinate_variances = np.sort(Sigma.diagonal())\n",
    "\n",
    "# Compute variances in eigenvector directions, in increasing order\n",
    "eigenvector_variances = np.sort(np.linalg.eigvalsh(Sigma))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEUCAYAAADeJcogAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABI4ElEQVR4nO3dd3hUZfbA8e8hhE5QihRBUFRAREU6CATYFVSaiGVtKAg2bKjYUBHLooKiKIiIvf0sWFbXLqiLbRUEBBQbqwLSlN5Ccn5/vHcmk8lMcifMZGaS83me+8yd97aTyzBn7r1vEVXFGGOMiaRCsgMwxhiTuixJGGOMicqShDHGmKgsSRhjjInKkoQxxpioLEkYY4yJKmlJQkTaiMhSERnsY90LReQzEZkrInNEpFsphGiMMeVexWQcVET+AZwEPO9j3VOAXkAPVc0RkTrAayJypqquSGykxhhTvpX6lYSI9AbaA6cCOT42uRwYrao5AKq6AbgFuCxRMRpjjHFKPUmo6oeqeqWq5gKVilpXRGoCoqprwxbNAY5JVIzGGGOcZD+4lmKWNwN+DC9U1T1AhogUt70xxpi9kOwkUZwsYGOUZduBqqUXijHGlD9JeXAdo2hXCxHL64rongpNyKmyIeJGdevWpV69evGKzRhjyoSvv/56vaoW+nJM9SSxCagVZVk1YEd4YTOgXt4Mxr5xHL16JTAyY4wpQ0Tkf5HKU/1208/AQeGFIpIJqEbp5zyTHK67DqwXdGOM2TspnSRUdTuwU0Tqhy3KBj6Ptl3drBy++AJeey2R0RljTNmXcklCREaIyPCQoruBB0Skkre8DnArcF+0fZw6xDW/uOEGyM1NYLDGGFPGpVySAOoDDQJvVPVt4F/AHBGZC7wKjFPV76PtoE/33TRrBkuXwtNPJzZYY4wpy6SsDV/aXkS/euQRnqo0grPPhgMOgOXLoXLlZEdmjDGpS0S+VtX24eWpXrupZHJyOP0cuOsu+PZbeOghuMw68TA+bN68mbVr15KT46fHGGPSQ2ZmJvvttx9ZWVkxb1tmk0RGBtx+OwwaBLfdBsOHQ82ayQ7MpLLNmzezZs0a9t9/f6pWrYo16DdlgaqyY8cOVq5cCRBzokjFZxJ7z/sVOGAAdOkC69fDPfckOSaT8tauXcv+++9PtWrVLEGYMkNEqFatGvvvvz9r14Z3g1e8spkkdu8GQAQmTnRFkybBunVJjMmkvJycHKpWtZ5eTNlUtWrVEt1GLZtJIuRE9OgBxx0HW7e6207GFMWuIExZVdLPdplPEgD//Ke7qpg+HX7+OUkxGWNMGioXSeLII+HMM13xuHFJismYMq5fv37B+UGDBiUxkuKdc845/PHHHwBMmzaNd955p9SOnernJly5SBIAt94KlSrBc8/B/PlJiMmYcuTpOLZinTx5ctz2Fcm5557L3/72t4Tse968eXz+ecEehOJ5bkpD2UwS3oPrUE2bwujRbv6aa0o5HmPKmZpxqm+em5vLokWL4rKvaKpWrUpGRkZC9j0/wi/SeJ2b0lI2k0SUJ/jXXw+1asH778N775VyTMYkiKoyceJEOnXqRM+ePenVq1ewquOaNWs49dRT6dSpE+3bt+fCCy9k586dwW337NnD2LFj6dixI506dWLAgAGsWLEiuHz8+PF88cUXDBgwgC5duvDcc88B8MEHH9ClSxe6detGr169mDdvXoGY+vfvH5w/99xz+fTTT+nXrx89evSgS5cuzJ49O7h8w4YNnH/++fTt25fs7GzOOOMMdu7cydatWznttNP417/+RXZ2NieffHJwmwULFnDsscfSo0cP+vTpw9y5c6Oenz///JNTTz2VDh060K1bNyZMmFBg+aRJk4Lbz507l1mzZjF69Gi6devGlVdeCcD69es5/fTTg3/v/fffX2Afy5Yto3///nTp0oWePXsyadIknnvuOe69914uuugisrOzefPNNwudG4C7776bjh070qVLF/r06cM333wTXPb444/z9NNPc+aZZ9K9e3e6du3KuJB75hs2bGDo0KH06tWL7Oxs+vbty8/xfvCqqmVqageqI0ZoNBMnqoJq27aqublRVzPl0NKlS5MdQok88MADOmHCBM3Ly1NV1dyQD/agQYP0008/Db5//PHH9fLLLw++v/3223XKlCnB9/Pnz9c+ffoE93XzzTdr//79ddWqVcF1/vzzTz3mmGN0/fr1qqq6a9cuPe+887Rv377BdU444YTg/LBhw3TIkCG6ZcsWVVXdtm2b9ujRQzdv3qyqqrt379bVq1cH13/ooYd0+vTpwfeh+w1sP2jQIN2wYYOqqm7dulVPOOEEXbduXcTzc9ZZZ+k777wTfP/iiy9q7dq1g8e8++67dc6cOaqqOmfOHD3yyCP1k08+KbCPYcOG6eLFi4Pvx4wZE9zn9u3btVevXvrTTz8Fl+/Zsyd4/j777LMC+wo9N88884yOGTMmeL5//vln7datm27fvl1VVR977DFt2bKl/vzzz8Ftzj//fJ07d66qqo4dO1bffvvt4LItW7bo7t27I54H1aI/48BXGuE7tVxdSQBceinsvz8sWADPP1+KMZm0JJKcKRazZ89m3LhxwSqOFSq4/9bLli2jUaNGdOnSJbjusGHDWLJkCbt27QLgvffe47KQPmvatm1L165d+fjjj4Nlffv2pWHDhsH3zz33HJdddhl16tQBoFKlSlx88cVFxnjxxRdTo0YNAKpVq0bPnj1Zvnw54LqMaNAg2Kcn/fv3L/IW03PPPcdFF11E7dq1AahevTpnnXUW70W4PbB582Y2bdrEscceGywbOnQoTZo0ibr/gw8+mGOOOSb4/scff6RRo0YcfvjhwbKxY8fy8ssvB+MZPnw4Bx2UP/SN39tXTz31FLfffnvw3+7AAw9k+PDhvPDCC8F1Tj75ZA488MDg+1NOOYWvvvoKgPr16xdo+1CjRg0yMzN9HduvMtstRzRVq8KECTBihOtK/KSTrPM/k762bNlCvXr1ItaBX758OUcddVSh8tatW7NixQoaNGhQ4Ms/oF27dnz33Xf07NkToMCXI7gvzT59+hQoO/jgg4uMs3nz5gXeZ2VlsWXLFgC2b9/OlClTmDt3Lrt37yYnJ4dDDjkk6r6WLl3KI488wh133BEs27lzJ8OHDy+07ooVKzjssMMKlRcVb/jfu2zZMl555RU+/fTTAuWBpLBkyRLOO++8qPsriohQpUqVAmXt2rXj+ZBfsEWdu8suu4xbbrmFOXPmcOaZZ9K2bdsSxVEU30lCRNoBpwJNVPUfXlkboIqq/jfuke2NCA+uQ519Nkye7LoSt87/TFFSvZPk3NzcqI2kdu/ejUb4A1QVESl2eUD4l1hOTg4VKxb86ghcvURTuYhfYhdffDEDBgzg2muvpUKFCqxZs4Zriqhdkpuby6xZsyJ++YeLFGtx8Yb/vbm5uVxyySVcdNFFUeMp7u+PZneE76rw81/UucvIyGDChAls2bKFRx99lPvuu48HH3yQ6tWrlyieSHz9ZSJyCvBvoApwSsiiOsCDcYsmXopJEhUr5nfXceutsGlTKcRkTALss88+rFy5MuKXfevWrfn6668LlS9ZsoRmzZpRt25dVq1aVWjb+fPn06JFi6jHbNasGcuWLStQ9uOPP5bwL4DVq1czZMiQ4BftDz/8UGB5eHwtWrTgv//197u0WbNmfPfdd4XKY4m3uOO1bNkyePsnnIiQl5cXddsKFSqwffv2AmXFnf9IatasyWWXXcbAgQN59tlnY9q2OH7T383AAFW9FAj92TIfaBnXiOIhpPZGNP37Q/fusGED3H13KcRkTIKccMIJ3H777cH3qkpubi6HHXYYv//+e4HnC08++SRt2rShUqVKiAi9e/dm0qRJweXffPMNn332GT169Ih6vH/84x9MmjSJjRs3Au7X8LRp04r8xVuUXbt2sXnzZgD++usvZs2aVWD5pk2bCtx3P+OMM5g+fTo//fRTsGzbtm3s2LGj0L7r1KlDXl4en3zySbBs9uzZwYZ0frRq1Yr169fz7rvvBsvy8vL4888/g/FMmzatQK2wQLwNGzZk8eLFUfd9xhlncM011wQTyYoVK3jsscc45ZRTom4T6pdffinwfsmSJXG9igD/SeJAIPAkKTStV45hH6XHeyhXFBG48043f889sGpVgmMyJkGuuuoqtm/fTrt27YLVIANfYE899RTTp0+nc+fOdOzYkc8//7zAvfxx48axceNGOnbsSOfOnbnpppuYOXNmkf38NGrUiJtuuokhQ4bQs2dPjj/+eIYNG0azZs1KFP/48eM5/vjj6dGjByNHjmTChAkFEs7IkSPp2rUr/fv3Z8uWLWRlZfHAAw8watQoOnfuTO/evTnvvPPYtm1bxP3PmDGDhx56iOzsbLp3786yZcu45JJLYopx1qxZPPbYY3To0IFevXoxcODA4NVIzZo1uf/++xkxYgRdunShd+/ePPigu8Fyyimn8MILL9CnTx8eeeSRQvsdNmwYzZs3p3PnznTp0oXzzjuPqVOn+u5o8sknn+SYY44hOzubY445hm3btnHaaafF9LcVx9fIdCLyH+ABVX1eRHJVNcMr/yfQVlX7Fb2H0tNeRL/q2BG++MLX+iedBLNnw6hRMGNGgoMzKW3ZsmW0atUq2WEYkzBFfcajjUzn9yrgYuBOEZnt9iUzROQb4HTg8pKFm0A+riQC7rgDMjJg1iyIcOvSGGPKNV9JQlUXAi2AN4B7gO24B9YtVTX1vlp9PJMIaNECzjsPcnPhuusSGJMxxqQhX1VgRaSqqu4AHg0rrygilVXV/0/30hDDlQTAzTfDU0/Bq6/Cp59C166JCcsYY9KN39tN80SkQ4TyTsAnEcqTK4YrCYCGDWHMGDc/dmzq1403xpjS4jdJtASWRCifD7SOXzhxEuOVBMDVV0PdujBvHvzrXwmIyRhj0pDfJLEFKNx+35XF9rO9NMR4JQGQlQU33ujmr70W9uyJc0zGGJOG/CaJF4FHRSTYcE5EWgGzgJcSEdhe2bmzRPeMLrgADjoIli2Dxx+Pf1jGGJNu/CaJq3CN6RaKyHYR2Q4sBH4AxiQquBJTLdGlQKVKEGi4etNNEKVtjjHGlBt+q8DuVNVLgH2BzkAXoI6qjlLV1PoqDXS0VYJbTgCnnAIdOsDq1a4ltjHGlGcxdamhqttVdZGqLlTVLYkKaq8EuhMowcNrcDkm0JfTXXfBmjVxissYY9KQ315gG4vI8yLyo4isDZ8SHWRM9vJKAqBnTxgwALZuhVtuiVNcxiTIOeecQ9euXcnOziY7O5uhQ4cGB/QJyMvLY+rUqfTq1Ys+ffrQq1cvTjnlFFauXFlgvQcffJCZM2eWKI5+/fJ757n22msL9RSbSIMGDSq1Y5U3fseTmInr/fVyXE2n1LWXVxIBEyfCm2/Cww+78SZi7LnXmFI1e/bs4Ohu3377Leeccw4ffvhhcGyEs846iy5duvD+++8HR01bsWJFcHS3gM8++4w9e/YwcuTIvYpn3Lhxce+NNOCZZ56hT58+BUaze/rppxNyLOP/dlNX4AxVfUNVPwqfEhlgzOJwJQFw2GH53XVce20c4jKmlBx++OH0798/OAbC888/T+vWrRk9enSBYTWbNWtWoLfR5cuX06RJExo1asSvv/66VzHUqFGjyJ5k98b8+fMLldWsWTMhxzKxtZOIW5fgInKSiMwTkbki8rGIFHmtKCJXeOt/JCKfiEj0fn4DH8y9TBIA48dD9equu47//Gevd2dMqWncuDFrvAdq//d//8ell15a7DZPPfUUw4YN4/TTTy924JoVK1Zwwgkn0KlTJ4455hhmhHWhPHr06OD4Co8//jhvvvkmp59+Ol27dmXy5MmAGwth0KBBwa6uQ8d1Bvj000/p06cP3bp1o0ePHjz99NNMmTKFZ599liFDhpCdnR0cVKl///7B7fbs2cPYsWPp2LEjnTp1YsCAAQXGehg/fjxvvfUWgwcPpnv37nTu3Jlp06YFl//0008cf/zx9O7dO9gteGDsjHJJVYudgKuBh4DKftYvZl9dgDlATe99Ndyod52irH++d+xM730l3O2vcyOt3656dVVQnTdP4+Hmm93uOndWzcuLyy5Nilq6dGmyQyiRYcOG6erVqwuUXXPNNfrNN9+oqurxxx9f7D7y8vL0tNNOC74fOnRokev369dP58+fH3w/adIkbdiwYfD9xRdfrL/88ouqqj722GParVs3XbZsWYHjnXjiifrbb7+pqmpOTo6efvrpumTJElVVXb16tfbp00fXrFkT3GbPnj1R/94TTjghOH/77bfrlClTgu/nz5+vffr00TzvP/DNN9+s3bp107Vr1waPPWDAAP3pp59UVfXUU08tEOuGDRuKPBfppKjPOPCVRvhO9Xt10BY4G1gnIt+KyJehU4x5aQwwRr3aUaq6HbgMiDao7UBgoqrmeOvvBu7wyguL45UEwFVXQf368Pnn8PLLcdmlSSciyZlKKC8vj1deeYVffvmFI488MlhWnE8++YSePXsG37dt25aFCxdGXHfp0qU0b96ctm3bBssuu+yyIsd57ty5My1b5g9i+eGHH9KvXz8aN24MQMWKFbn00kt57bXXAJg2bRrjx49nv/32C24TequsKO+99x6XhQxc37ZtW7p27VpghL7hw4dTr1694LEHDx7MggULANhvv/3YE9LOKvy5TXnj98H1O94UDy1VdUFogar+ICLNRES8jBZqLdAEWBFSVh9YF3HvgQ/qXj64DqhRw9VwuuAC92xi4EDX6M6YVDJkyBAyMzOpVKkSvXv35oknnggu85Mknn32WRYtWsTzzz8PwM6dO9m4cWMw0YT68ccfadOmTYGyihUrFjky3eGHH17g/dKlS3n00UcL3Nbas2cPvXv3Di6/rgR992/atImGDQv3INSuXTu+++67YCJs3rx5geVZWVls2eLq5EycOJFx48ZRtWpVzjnnHA455JCY4yhLfCUJVX2i+LWKJyI1gY1RFq/C9QUVPpDonbguQUaq6hIRaQbcAER+LhGnB9ehRoyAKVPcoEQzZkCMIx+adJYmXQKH1m4KV6FCBXbs2BF1SMydO3eyefNmPv300wLlAwYMIC8vr9AVQk5ODhUrFv7qKOpKIlDLKiA3N5fbb7+d448/PuL6ubm5JXrwvXv3bgr/znS31UP3V9R43NWqVeOee+5h/fr1TJ8+nXXr1nHvvff6vpIpa3w/jBaRRiLSW0SOD59iOF4W0ZPERmCf8EJ1gxoNBZ4RkY+AZ3HPI1ZE2smu3bsBuP7KK2nfvn2h6eGHH44hXKdixfzxsG+5BTZtinkXxiTNwIEDmTVrVtTlr732Gn379i1U3rVrV+bOnVuovFmzZoXaQOTm5vLLL7/4jqlFixbB2lfRlgceSocTkahXR3Xr1mXVqlWFEsX8+fNpEWM99rp163LjjTfSrFkzPvzww5i2LUv8NqYbCfwCPAn8C5gGvAo8B4yI8ZjRfh5ELBeR/YDpuIfXZwEfAM+LyP6R1q/s/Vq646ab+OqrrwpNo0aNijFcZ8AA6N4dNmzITxjGpIORI0fy6quv8uyzzxb48gzcd3/hhRcYPHhwoe1OPvlknnnmmULlRx99NPPnz2fp0qXBsqlTp8b0S/vYY4/l/fffL5AIdu/ezebNmwG48MILueGGG1i/fn1weU5ODgANGzZk8eLFEfcrIvTu3ZtJkyYFy7755hs+++wzevTo4Su2n3/+OTifl5fHsmXLEtbmIx34vZK4Dhisqo0BUdVmQG1craTZMRxvE1AryrJa3vJwTwLjVPUhVf1VVW/04nku4l7i1Jgu0m4Dn7t774Xffovr7o1JmIoVK/Lqq6+yYMECOnfuTHZ2Nn/7298466yzWLt2LapKrVqF/1sefPDB/PTTT+wMu3UrIjz66KNcf/31ZGdn06NHD2rUqBFTq+eMjAwee+wxxo8fT8eOHenduzennnoqa9e6DhyaNm3KDTfcwODBg+nWrRu9e/fmZa/myMiRI5kwYQJ///vfef311wvte9y4cWzcuJGOHTvSuXNnbrrpJmbOnOn79tVdd91F9+7d6dWrFz169KBVq1Z0LcfDVUqk+3eFVhLZBuyrqru9+TqqulNEmgCfeEnD3wFFFqhq2wjlC4G2qpoXUrYP8LqqFvoJ4N16Gqyqf4WWt69fX79au9Y9RAip4RAvp54KL7wA55wDjz0W992bJFq2bBmtWrVKdhjGJExRn3ER+VpV24eX+72S+Bk40Jv/DTjCm98M1IsxzsUicnRYcIcAv4cmCM8OIDPKfqoCWwuVxrkKbLg77oDMTHjiCYhSQ9AYY8oMv0niEaCZN/9/wFMicg3wPLGPcT0JuFtEsgBEpBowFbjbez9CRIYDqOouYJWIDAtsLCIVROQqYH6g7UTBvyi+VWDDNW8OF13kKr1cE61lhzHGlBF+x5O4T1UD7SRuwT0nOAb4Hvcw2TdVXYRLFG+IyFzgfeAJVZ3rrVIfCK3Ldx7QzuuOYw4wF3cVEfleUoKvJADGjXPDnb7zDrz3XsIOY4wxSee3MV2Qd0vo9r05qKq+BbwVZdkdYe//AorveCYgcCWxY0eJ4ytO3bpw/fWucd3VV8P8+fmHNcaYsqTIJCEig1X1VRE5u6j1VPXJ+Ia1FwLf1tu3J/Qwl14KDz7onks8/TScXeQZMsaY9FTclcTZuPYQRbUxVtztp9RQClcSAFWrwm23wbBh7vbTySe7MpPewlvmGlNW+KnJGkmRN0lUdYj32qGIqWOJjpwopXQlAXDGGXDkka7NxP33J/xwJsEyMzPZkeAfF8Yky44dO8jMjFZZNDq/La7niEh6/E4uxSSRkZE/HvYdd0BI41CThvbbbz9WrlzJ9u3bS/yry5hUo6ps376dlStXFuhV1y+/D65bAzVw7RZSWykmCYC//x2OPRbefdfdfpoypVQOaxIgKysLgFWrVgW7gDCmLMjMzKR+/frBz3gs/La4vhw4GhgRsW1CCmnfqpV+9d130KEDfBnrUBcls3AhtG3rOgJctsy1pTDGmHSyty2uawPdgP+JyBsi8kLoFNdI91YpX0mAey5x9tmQk+OqxhpjTFnhN0nk4mowzQD+CywJm1JHKdVuCnfrrVCliuvX6YsvSvXQxhiTMH4HHbol0YHETRKuJACaNIHLL4eJE10Du48+2qtRKI0xJiXEMuhQRRE5SEQOC58SGWDMkpQkwLXArlMHPvkEIvRgbIwxacdvFdiBwHrgR2Ax8K03zQeeSlh0JZHEJFGrFtx0k5sfO9Y9ozDGmHTm90riTuBqoDKuRlQFoBVulLjU6gtVxDVg2LMnKd/SF1wAhxwCy5e78bCNMSad+U0SBwBPetVfd4lIJVX9HrgCeCBh0ZVUtWruNQmtZytVgrvucvPjx8PGjaUegjHGxI3fJPE7+d13rwJahJQ3jXdQey2QJJJwywlg0CDo0cONh33HHcWvb4wxqcpvkngJCAw5+i/gIRE5FTdY0IJEBLZXAj3tJSlJiMDkyW7+vvvgl1+SEoYxxuw1v4MO3aCqr3pvbwCWAuOAOsCwaNslTZKvJADat4czz4Tdu12tJ2OMSUd+azcNFZHKAKq6TVVHqmobVR2sqj8kNsQSSIEkAe5WU6CB3WefJTUUY4wpEb+3m+4C1ojIoyLSR1K9w/0USRJNmsCYMW5+zBg3LrYxxqQTv7ebDgL6An8BjwK/i8hkETk6kcGVWOCZRAqMDXDttbDffvD55/Dii8mOxhhjYuO7xbWqfqGqV6pqU2CoV/y8iCxNTGh7IUWuJABq1nT9OoFLGLt2JTceY4yJhe8kEaYikOFtnxe/cOIkhZIEwPDh0Lq1q+U0dWqyozHGGP9i6bups4jcKyIrgWeAXcBJqnp4wqIrqerV3evWrcmNw1OxIkya5OZvu81GsDPGpA+/tZv+B/wbNzrdGUBTVb1GVRcmMrgSq1nTvW7Zktw4QvTr50aw27QJJkxIdjTGGOOP3yuJy4EGXtXXuZrqAwCnYJIAdzVRoQJMnw7ff5/saIwxpnh+aze9oqq7Ex1M3KRokmjTBkaMcH0Pjh2b7GiMMaZ4JX1wndoCg31v3pzcOCKYMAFq1HDjTcydm+xojDGmaGUzSaTolQRAgwZwjde5+pVXQl7q1Q0zxpggSxJJMGYM7L8/zJ8PTz+d7GiMMSa6qGNci8iLgK8H1Kp6StwiiofA7aYUTRLVqrl+nYYNg+uvh6FD85t2GGNMKomaJHDDk6anFL+SANdD7H33uauJe+6BceOSHZExxhQWNUmo6i2JOqiInASMAXJwt7wmq+prxWxzLvAPQHBxX6SqyyKuHEgSKfjgOqBCBTfmRK9eMHEinHeee15hjDGppKgriQJEpCJuGNMq4ctU1Xf/TSLSBRgN9FPVLSJSDXhJRP5Q1S+ibHMPbhS8E7whVIuWBlcSANnZbhS7116DG2+EmTOTHZExxhTkt8X1QGA98COwGHcr6ltgPvBUjMccA4xR1S0AqroduAy4JsqxuwGbVPUeXwkCUv6ZRKi77nLddjz6KCxenOxojDGmIL+1m+4ErgYqA6KqFYBWwAdE+XIvQktVLTDkqTdwUbMo41QMB6bEdITKlSEzE3JyUr7b1UMPhQsvdFVhr7oq2dEYY0xBfpPEAcCT3i/5XSJSSVW/B64AHvB7MBGpCWyMsngV0DBCeT0gQ0QeEJEPReRdETmn2IOlwXOJgJtuglq14N134e23kx2NMcbk85skfgcCj1VXAS1CypvGcLwsoieJjcA+oQXekKl1cSPj3a+qvYEhQB8RuaTII6XJcwmAunXzazdddZXrtsMYY1KB3yTxEtDWm/8X8JCInApMBRZE3SqyaEOfRiqvDbQErlLV5QCquhU435sKWbduHe3bt+eHtWsB+Ef//rRv3z44PfzwwzGGWzouuQQOPBCWLLEH2MaY1CGxdugqItVxzwg6Az8BV3vPFPxsWwN4S1W7R1j2BnC+qq4MKWsAPK6q/SKs/yFuPIu/Qsvbt2+vX331FRxzDMybBx9/DN0LHS4lzZ4NJ50EderADz/AvvsmOyJjTHkhIl+ravvw8pi75VDVbV6X4W1UdbDfBOFtuxU3JkUkTYDVYWXrgFpR1s/EDXwUWeAb9s8//YaXdCeeCD17woYNNuaEMSY1+GonISLHF7VcVf8dwzEXi8jRqjo/ZP+HAL+raoHu7lQ1V0Q2iEgbVV0csn51oKpXfTayOnXc64YNMYSWXCIwZQocfTQ88ABccAG0aFHsZsYYkzB+G9M9GPa+Gu6B8i7gU9yodX5NAu4VkRNVdbPXmG4qMBFAREYAqqqPeuvfCjwuIgNUdZXXqG8yML3Io6RhkgA46ijX+nrmTNdL7BtvJDsiY0x55itJqOqB4WUiUhf35f5hLAdU1UUiMgl4Q0TygErAVFWd661SP2z9L0TkOuBFrxlFVeBZVZ1V5IFq13avaZYkAG69FZ5/Ht58E955B/r2TXZExpjyyne3HOFUdb2IXIm7koipPo6qvgW8FWXZHRHK3gXejSnANL2SAKhf33XTMXYsXHEFLFzo2gYaY0xp29vxJLbhGrulnkCSSKMH16EuvRSaN4dly+Chh5IdjTGmvPLbd9NhEaaOuNbWXyc2xBJK4ysJcD2LTJ7s5m++OW1znTEmzfm9kviWgh37fQu8j6u2ekFiQttLaZ4kAAYOhN694a+/YPz4ZEdjjCmPfCUJVa2gqhnea2DKUtXsWLoJL1VlIEkEqsRWqADTpsHS1DzTxpgyrGyOcQ0Fk0SMrcpTSZs2MGoU5Oa6sbHT+E8xxqShqN1yiMijERdEoKrD4xbRXgp2ywFu4OgdO1wnfzWiNfROfevWwSGHwKZNrt3ECSckOyJjTFlTkm45NoRM24HTgebAbmAr0Ag4C/hf3KONlzJwywmgXj3XnTi4q4kcf0MvGWPMXouaJFT16sAE1AQuUdWeqnqBql7qdbp3LnBEaQUbszRuUBdu9Gg3QNHy5fBgePt3Y4xJEL/PJAYDz0cofwVI3fbAad5WIlSlSvlVYm+5BdavT248xpjywW+S2A0cEqH8UGBn/MKJszJyuynghBPg73+HjRvzbz8ZY0wi+U0SjwCviMgFItJJRDp6HfG9BPh+wF3qyliSEIF774WMDJgxAxYvLn4bY4zZG36TxPXAbcAI4G1cP0oX4Xp0vTYxocVBGUsSAK1buy7E8/Jcv05WJdYYk0h+G9Opqs5U1Q6quq+q7qOq7VR1evgYECllv/3c65o1yY0jzm65BfbZBz74AF5/PdnRGGPKsiKThIikd9+jjRq519XhA96ltzp18rvpGDMGdqbuUyFjTJor7kpiCYCI5IlIboQpT0RySyHOkmnY0L2uWpXcOBLgoovgsMPg559h0qRkR2OMKauKSxJ9vNc2uPYQ4VOgPDUFriTKYJLIzISpU938HXfAr78mNx5jTNlUZJJQ1d+81yVFTaUTagkEriRWry6TT3h794aTT3Y9j1x1VbKjMcaURX7HkxgrIh1C3l8uIhtFZIGItEpceHupalXYd1/Xj0UZquEUatIk92e++CJ8GNNAssYYUzy/VWAvAHaAG4AIGA+cCLwMTElEYHFThp9LABxwAFx/vZu/5BLr18kYE19+k0RDYLk3fwkwTVXn4BJExwTEFT9l+LlEwFVXwUEHufEmrF8nY0w8+U0SfwBNRaQRcBpu2FKADG9KXeUgSVSp4gYnAjfUaRlrFmKMSSK/SWI6MBc3nvUMVQ1843Yk/wojNZWDJAHQvz8cdxxs3gzXpm4beGNMmvHb4vou3HgS/1DV0K+gtcD5iQgsbkJrOJVhInDffa632Mcfh88/T3ZExpiywPfwpar6EfBfEWkZUrZQVb9OSGTxUk6uJMCNXjdmjJu/5BLXv5MxxuwNv1Vg64jI68AW4NuQ8tNF5L5EBRcXgSSxcmVy4yglN9wA++8PX30Fj6Zu/7zGmDTh90piKvA7UC9sm8+Ak+IdVFwdcIB7LSdNkmvUgLvvdvPXXQd//ZXceIwx6c1vkjgWuEZVNwChTZdXA3XjHlU8NWwIFSu6Kj87diQ7mlJx2mnQo4cbve7GG5MdjTEmnflNEnuAfSKUHw2kdoXLjIxydzUhAg884P70adPcrSdjjCmJWEame1xEjgAQkQYiMhh4EpiWoNjip1kz97piRTKjKFVt2sDll7suq84/H3JTt69eY0wK85skbgReB94ABFiJa1B3n6remaDY4qccJglwY040aQLz57srCmOMiZXfJHEwcL+qHgBkAbVVtbGqTk1caHHUtKl7/d//khtHKatRI7878RtuKBe1gI0xceY3SSwEqgCo6lZV3bQ3BxWRk0RknojMFZGPRWRQDNteLCKLYzpgOb2SABg0CAYOhC1b3JjYxhgTC79J4jvg8HgcUES6AKOBfqqaDfQDzheRTj62PQTXFUhsDfjKcZIAuP9+qFYNXngB3n472dEYY9KJ3yRxDjBZRAaKSF0RqRY6xXjMMcAYVd0CoKrbgcuAa4raSEQygEnA1TEer9wniaZN88fEvvjiclMT2BgTB36TxBtAJ+BVXJXXLd601XuNRUtVXRBaoKo/AM1ERIrYbizwpKqujfF4rtV1Robrv2nnzpg3LwsuvxwOP9yNiX377cmOxhiTLvwmie7AocCBwEEhU+C9LyJSE9gYZfEq3LgVkbY7AjhUVV/2e6wCKlbMv5r4+ecS7SLdZWbCjBlu/q67YNmy5MZjjEkPFf2spKrxqhaURfQksRHXYK9AHRwRqQTcCZzh5wDr1q2jffv2hcqnrFvHMcC799/PsQ895DvgsqRrVxg5EmbOhAsvhDlzXMM7Y4yJxleSiLNoX0vRyscDU1X1Tz87r1evHl9FamJ85ZVwzz0cG6gOW05NnAivvAIffQRPPQVnn53siIwxqcx3V+FxsgmoFWVZLW95kIh0xrXJ+PdeH7ml18P5d9/t9a7SWe3aMHmym7/ySvjTV+o1xpRXUZOEiIwUkRrxPJiqbgWi7bMJrsPAUBcAR3rtKYIT0N+bf9L3wS1JBJ11FvTs6ToAtFHsjDFFKepK4l68W0AisiiOx1wsIkeHFnjtH35X1QLD5KjqOaraRVWzQyfgW2/e/82S0CShWvS6ZZwITJ/uHmbPnAn/+U+yIzLGpKqiksSfuF5eIU4N6TyTgLtFJAvAa2cxFbjbez9CRIbH8XhO3bqw775uEOg//oj77tNNq1ZwjdcyZeRI2LUrufEYY1JTUUliAvCuiPwJqIisjTbFckBVXYRLFG94t47eB55Q1bneKvWBBjH/JcURsVtOYW64AVq0cKfD2k4YYyKJWrtJVR8RkVdx7SNeB06O10FV9S3grSjL7vCxfXaJDtyqFXz2GSxZAr16lWgXZUmVKu52U48e8M9/wsknuy7GjTEmoMgqsKq6HlgvIter6kelFFPiHHGEe10Uz0cs6a17d7jgAnjoIXfbad481zjdGGPAZxVYVX0YQER6ichob+qT2NAS4Mgj3evChcmNI8VMnAj77w9ffAEPPpjsaIwxqcRXkhCR/UXkv8AzwABvekpEvhaRxokMMK4CVxKLF9tQbSFq1coflOj668vdsBvGmCL4bUw3FfgSOEBV+6pqX6Ap8Lm3LD3Urg2NG7tuUH/8MdnRpJSBA90ziW3bXJcd5byWsDHG4zdJ/A24WVX3BApUNQe4BUiv2052yymq+++HffaBt96C555LdjTGmFTgN0nkApkRyjOBnPiFUwoCScIeXhfSoEF+lx2XXeZaZBtjyje/SeJN4D4RCfa7JCL7AvcBsxMRWMIcdZR7/Tq2we3Ki3PPhd69XYIYMybZ0Rhjks1vkrgEqAasFpFlIvIdrkvvysDlCYotMTp5o6R+8YXdeI9AxI07UaWK6yX2nXeSHZExJpn8VoH9S1X7A0fhhhm9AThCVQeo6rYExhd/TZpAw4bw11/www/JjiYlHXwwTJjg5keNcj2ZGGPKp5i6ClfV5ar6uqq+7A05mn5EoHNnN//558mNJYVdcQW0awe//uq6FDfGlE+lPZ5EarAkUayKFeHJJ6FyZXjkEXjzzWRHZIxJBksSJqrDDsvv+O+882DDhuTGY4wpfeUzSbRr5zooWrTItR4zUV1+uevf6Y8/YPToZEdjjCltvpOEOI1F5LDwKZEBJkT16q6LjtxciDQetgnKyIDHHnOn7Pnn4YUXkh2RMaY0+e27qSewEvgV+DZsWpyw6BKpa1f3+vHHyY0jDTRvDpMmufmLLrIxm4wpT/xeSUwBZgFZqlohbErPjqV793avH36Y3DjSxPnnw7HHuucSw4dDXl7x2xhj0p/fJNEC+Keqbk1kMKUqO9tVh/30U9fhnymSCDz6qOsj8a23XD9Pxpiyz2+S+BnYP5GBlLrataFtW9i92yUKU6z994dZs9z8NdfAggXJjccYk3h+k8SVwAwR6SYidUWkWuiUyAATKnDL6YMPkhtHGhk82HUlvns3nHYabC0715bGmAj8JonHgGzgE2ANsMWbtnqv6cmeS5TI5MnQujUsX+56izXGlF1+k0QX4EBvOihkCrxPT927u6bF//2v68vJ+FK1qqsOW6WKe07xf/+X7IiMMYnit4O//xU1JTrIhKlRA3r0cFV1rN+JmBx+ONxzj5sfOdL6SjSmrIqlMd0RIvKEiHzlTU+JyJGJDK5UDB7sXl95JalhpKMLLoChQ2HLFjjpJNi+PdkRGWPizW9jukHAx8A6YLI3rQHmesvSVyBJvP22VYWNkYir7XToobB4sWtoZ0N0GFO2+L2SuA0YpqpXqepz3nQVcK63LH01aQLt27ufwe+9l+xo0k5WFrz0kntO8cQTrsdYY0zZ4TdJHAK8H6H8PeDg+IWTJHbLaa+0aeNGswO45BKYPz+58Rhj4sdvkliJG5UuXDtcf07pbcgQ9/rKK7BzZ3JjSVNnneW67ti1yz2f+PPPZEdkjIkHv0niHuB5ETlfRDqJSGcRuRh4DrgjceGVklat4OijYdMm+Ne/kh1N2poyxd25W7HCNbTbsyfZERlj9pbfKrAPAlcAw4C3cLeZ/gFcpqpPJC68UnTWWe71qaeSG0caq1IFXn4Z6tVzj3euvTbZERlj9pZoGauO0r59e/2qJGNErFnjOicSgd9+gwYN4h9cOfHxx9Cnj7uSePppOOOMZEdkjCmOiHytqu3Dy4u8khCRzAQFc5KIzBORuSLycVHVaEXkABF51FvvIxF5X0Q6xj2o+vVh4ED3zTZzZtx3X5706AH33efmzzsPvv46ufEYY0quuNtNSwBEJE9EciNMeSKSG8sBRaQLMBrop6rZQD/gfBHpFGWT84CpqtpDVXt62z4hInVjOa4vF1/sXh96CHJy4r778uTCC12C2LkTTjwR1q5NdkTGmJIoLkn08V7bAEdEmALlsRgDjFHVLQCquh24DLgm0sqqepOqLgh5/x3wPNArxuMWr3dvaNkSVq2C116L++7LExF44AHo3NndvTv5ZNdzrDEmvRSZJFT1N2+2uqouCZ+A5UBWjMdsGfql7x3nB6CZiIjPfawFasZ43OKJuGbDAA8+GPfdlzeVK7sH2Q0buucU559vLbKNSTd+q8B+FqU8D1fbyRcRqQlsjLJ4FdDQ5646A9/4PW5Mzj4bqleHuXPhm8Qcojxp1MhdlFWtCo8/DrfemuyIjDGxKO7BdTURqe5mpWrYYEM1gKFALB0eZRE9SWwE9iluByLSHGiiqolp11urluvWFOwbLU46dHBdi1eoADffDE8+meyIjDF+FVkFVkQ+wv1qrxhllR3AJar6mK+DiewPzFDV/hGWPQPcrqpLi9i+CvA2cFG09Zo2bar16tWLGsOoUaMYNWpU0YGuXg0HHeSeui5cCEfE+tjFRPLAA67bjooVXX+KffoUv40xpnREqwIb7csfAFXtKSL74J49dAhbnAv8oaqxtKvdBNSKsqyWtzwiEakAPAlMKSqR1KtXjxK1kwjVsCGMGgX33w+33OJurJu9Nno0/PKLG4di8GA3IGCH8E+VMSalFPtMQlU3AtkRBhv6HUBEKvs9mKpuBWpEWdwEWF3E5g8Ac1X1Vb/H2yvXXOOaEM+eDf/5T6kcsjy4+244/XQ3Nna/frA0aro3xqQCvw+unxaRSL/5OuHGvY7FYhE5OrRARA4BflfVvEgbiMgEb/m0GI9Vco0awdVXu/nLL3ej15m9VqGCe4B9wgmuE8C//91dXRhjUpPfJNESr2FdmPlA6xiPOQm4W0SywD0cB6YCd3vvR4jI8MDKInIJUFVVS78jwWuuccni66/dYAkmLjIz4cUXXcvsVatcolhd1DWkMSZp/CaJLUSuntoQiKlvbVVdhEsUb4jIXNw4FU+o6lxvlfpAAwheYUwGOnhdeIROvqvellj16nDnnW7+qqtc/04mLqpWdR3uHn00/PQT9O1r3Ysbk4p8dfAnIg/gWlef77V4RkRaAdOA5ap6fkKjjEGJO/iLRtXdPH/3XTeg84svxm/fhnXr3BXFd99B27au99g6dZIdlTHlT4k6+AtxFbAIWCgi20VkO7AQ+AHXzUbZJQIPPww1arhxOl96KdkRlSn16sH778Mhh8CCBa5nlHXrkh2VMSbA73gSO1X1EmBfXLuJLkAdVR2lqtsSGWBKaNo0/7bT+efDr+k/GF8q2X9/18C9RQtYtMglCusQ0JjU4PdKAnCd8anqIlVdGOigr9y44AI4/nh34/zUU623ujhr1AjmzHGDBH77LfTqBX/8keyojDG+koSItBCRf4vIryKyNmwqH09zK1Rw/Uk0aQKffw5jxyY7ojKnYUOXKFq3du0nsrNd7SdjTPL4vZKYCfwBjADqAqcAN+J6Y+2bmNBSUJ068MILrl+J++5zzypMXNWv7xJFmzbw/ffQs6e1ozAmmfwmiaNxfTS9h+uO42NVnQFcBExJUGypqXNnmDHDzV90kav1ZOKqXj3XZUfbtvDjj9C1q+tCyxhT+vwmib+AKt78WqCxN/8lrtV1+TJ8OFx3HeTmwkknwbx5yY6ozKlb1z3MDjyb6NHDvTfGlC6/SWIO+aPU/Qe4QUTqA2cCKxMRWMq77TY488z8Toisf6e4y8qCt95yzVM2b3YN7qwGsjGlK5Z2EoHxPMfhhg5dDUwELk1AXKkv0AnRGWfkJwr7qRt3lSu7sSguvthVKDv5ZLj9dhvhzpjS4jdJtAH2gBtqVFUPBeoB9VT134kKLuVlZLg+nc44A7Ztcz91n3022VGVORkZMHWqa6oiAuPGwVlnueE+jDGJ5TdJ/AuoFFqgqhvUT58eZV0gUVx6qfupe8YZ7laUnZq4EnG1jl95xXWp9cwzrorsb78Vu6kxZi/4TRL/BbonMpC0lpHhqsTee6/7NrvxRhgwANavT3ZkZc6gQa6ewAEHwBdfuBpQb76Z7KiMKbv8JonhwG0icqmItBeRw0KnRAaYVi6/HF57Dfbd131zHXWUq8tp4urII+Grr9xjoA0boH9/10mvNYI3Jv78JonlQHtcm4gvgW9DpsUJiSxdDRgA33zjKvevXOkGcr7gAlc9x8RNvXouD0+c6C7kJk+G7t3h55+THZkxZUvUJCEiR3jjSqOqGapaIcqUUXrhpokDDoCPPoIJE9wIOzNmuL4mXn3VnlXEUYUKblyojz92vaV8+SUccYQ73XaajYmPoq4kPsR1wYGI3F864ZQhFSu6ZxPz50OHDvD773Diia512Pz5yY6uTOna1V28nXKKq2R2wQVw3HHuQs4Ys3eKShICVPTmLy6FWMqmww+HTz+F++93fT999BG0b+/qcH7/fbKjKzNq14b/+z/XpqJ2bXjnHXfqH3rINYw3xpRM1JHpRGQ6MAD3DGIwELWtq6qekojgSiLuI9PF08aNriXY/fe7p6wirtvx6693PdqZuFi9GkaNgjfecO/btYMHH4RO5a8DGWN8i3lkOlW9ENeB31fALmBJEZPxY5994O673VidI0e6W1LPP+9upB97rHtmsWdPsqNMew0bwuuvu5FmGzeGr792/TKOHGm1ko2Jld8xrh9W1VGlEM9eS+kriXC//eaSxsyZ+c2HGzd232Znnw3NmiU1vLJg61Z38TZ5MuTkuNrJN97oOvCtXDnZ0RmTOqJdSfhKEukkrZJEwJ9/ulbb06fDDz/kl3fv7joRHDrU3Wg3Jfb99zB6tBtPG1wFtAkT3OnNsPp5xliSSAt5ea7x3axZrlHejh2uPDMT/vY3GDjQTY0aJTfONKXq2lZcd50bIhWgZUvX3ccZZ0ClSkVvb0xZZkki3WzZ4joqevpp+OADl0ACOnRw/VMcd5xrfmw/hWOSm+v6frrpJvjf/1xZ48ZwxRUwYgTUqpXc+IxJBksS6WztWldV57XX4L338q8wwN1k79ULevd2rbtbtHC1pkyxcnJcvYGJE92Y2uA6DzzzTNc1uVU4M+WJJYmyYvt2d2P99dfda+CncECdOq4qT9eu0KWLu+qoUSM5saaJvDyXg6dMceNrB/To4eoQDBkC1aolLTxjSoUlibLq55/dc4wPP3SDHq1eXXB5RoarYtuhg+sy9aij3E/k6tWTEW3KW7oUpk1z9Qi2bnVlWVmuOcu557r8axdqpiyyJFEeqMKvv7oW3p995l6/+aZwk2MROPRQlzCOOsr1K9WqFRx4oD3f8Gze7MaPeuwx1ydUQIsWcNpprguQw6z/Y1OGWJIor7Ztc/1qL1jgEsY338CSJZEb7VWq5JJHq1au2k/gtUWLcn2/ZckSN1LtU0/BmjX55a1bu+FUhwxxXYDYFYZJZ+UmSTRt2lT/F36f3hS0axcsW5afNJYuda3AixrmrWFDOOigyFODBq5L1jJuzx53V++FF2D2bPjrr/xljRq50Wv79nW1levUKX5/Dz/8MKNGpUUb1bRh57Tkyk2SqF69um7bti3ZYaSnLVtcq7Nly1zSWLYMli1jz3ffBXt6jKhKFdc6vEkTNzVuXHg+K6uU/ojSkZPjaia/8AL8+98FrzBEXH9Rxxzjpm7dXB4N1759e+yqN77snJacJQlTYh3btePL2bPdQ/JIk58OkbKyXMJo3NhdlTRoAPXru9fQ+X33Tbv7NqqwaBG8/bbrffY//3FJJFTz5q7CWdu2rmnLkUdC3772hRZvliRKLqWShIicBIwBcnCdDE5W1deKWP9C4GxcR4MKjFPVeZHWtSQRf8X+x9u82VXF/e03N27Gb7/lT4H3oW07ipKZWTh51K3r7t/UqZM/H3jdd9+Ue9i+dSt8/rlLFvPmuToEkT6SmZl/0LdvA1q3hkMOcY+DDj0U9tsv7fJkyrAkUXLRkkSRdxESFEgXYDTQT1W3iEg14CUR+UNVv4iw/ilAL6CHquaISB3gNRE5U1VXlGrwJrKsLFetNlrrM1XXP1UgaaxZA3/8Ufj1jz9cwvn9dzf5IeISRXgS2WcfN9Wqlf8aOh94zcyMyykIVaOGey7xt7+593v2uCuNL7+EhQvdY6BFi2D79ga88UZ+l+YBWVkuWRxyCDRtmn/nLjDVrm1JxJSeUr+SEJEXgTtUdUFI2SHAnao6JML6nwKDVXVtSNnfgeNV9Yrw9e1KIv5K9dfZjh2Fk8eGDfnT+vUFX0OfHpdE1aqFE0eNGoWn6tWLL6te3dUQ8/ENnpsLRx55Ejfd9DLffw/Ll+dPGzcWvW21avmPeho3dlce9esXfN1vPzcOeAJyYEqzK4mSS5krCaBlaIIAUNUfRKSZiIiGZC0RqYlLZGvD9jEHuKMUYjWlrWpV9xDcbzfpe/a4RBGaODZscN+0mzblv4bOh5bt2OGm8EaIJVWxonuQX7Wqew2dDynLqFqVO/74iIFzL3JlTarAoVXRylXYllOJtX9l8seGTNZtymTdxkzW/FWJNRsyWb0+k03bM9m9vBI7lmeylEwWksluKpFDZnDaTSX2UJEqVStQPSuDGrUyqFmrAjVqufmsWkKtfYRatVxuq1q1YJjhU+XKLuEEpooVC85nZNjVTVlVqknC+9LfGGXxKqCh9xrQDPgxfEVV3SMiGeFJxZRDFSu6n8z16sW+rap7WBCeRLZuddO2bfnz0cpC32/Z4pJW4H0xBoLrHj6EADW86aDY/6LCdnjTmsKLcqlALhnkea/h84XfV2C3N69IoQnEZYrg5P4iFcl/DZ0nbJ1Cy9y8WyVsP4HjhmWmf26qxJd1jss/mVEUlc/C91lo3RLsN3yf8VE6Wbm0rySyiJ4kNgL7UDBJFLX+dqCq95pfuH37LhEpalTjdYCNTxabuiJi5yz+6pLUz2KeN8WRelPy1OVP+/9dQk0jFSbjdlO09BeXclWtEnNExhhjIirtZrKbgGi99dfylvtdvxruQtoYY0yClGqSUNWtuNutkTQBwp8e/kyEW7Mikul2Z88jjDEmkZLR4c5iETk6tMCrAvu7qha4Qaqq24GdIlI/bB/ZwOdh+zhJROaJyFwR+VhEBiUg9jJBRGqLyBQR+UREPvJe+4Wt08Nb9pF3XkdG2M/hIvJeyDo3iFgdFwAR6S8im8PKij1fItJERF71/k3micg9IlKuB1YVkUEi8paIvOP9/84OW26f1URS1VKdgCOAD4As73014G0g23s/Ahgesn4/4EWgkve+Di5BtAhZpwuuWmzNkH3+G+hU2n9fOkzAFUCvkPeNgG9w1ZMBDgS+Bup77zOAx4GTQrbJAhYCh3rvBbgVuDLZf1+yJ+8z+irwVCznyzvPXwKdQ8pGAlOT/Tcl8VxeAdwJVIuy3D6rif43SNI//HHAx8Bc4FPgHyHLrgeuD1v/bGCet/4nwN/Clr8ItA0rOwSYnewTnC4TMAy42pu/GzgxbHktYF7I+4uBK8LWyfD+w2Yk++9J8rl8AmgFPB7L+QJOAO6LsL85QN1k/11JOI8HAjOKWcc+qwmektK/s6q+pao9VDVbVbuq6nMhy+5Q1TvC1n9SVbt563dX1ffDdhmxgR7QzC4pfVsL1PTmuwPvhi5U1U3ANhEJNEjIBt4KWycX9x+vwO3E8kRE/gEsUtVlYYuyKf58FVrH8x6ua5ryZgRwbzHr2Gc1wdJ+EACfDfRM8TrjbjkBVFHVSH2bLMddoYFr6PhTMeuUKyLSCDiNyF9szSj+fB0E/FDMOuXJEcD/ROR2EfnAm8aISOj3ln1WEyztkwT+GuiZIojIvrhbgIGu5qLVGttI/vmsqKo5xaxT3twHjNGwChgeP+cr2mc5dJ3ypDZwD/CKqvbBfUbr424xBdhnNcGS0ZguEWJtiGc83u24WcB1qrq7uNX97HLvo0o/IjIKmKOqkX6xFrlpMe+jlZUHbYFhgXOqqrtF5FpgvojUUFelPhr7rMZJWbiSiLWBniloMu7L7QMf64aezxyvvUpR65QLInIgrhbe9CJW83O+on2Wy9059fwQnnTVPXVeQPG3ieyzGidpfyWhqltFJJYGesYjImOBHFWdGrZoq4hUj3CvtwX593Z/BA4Gwh/QtgD+E/dgU9sI3GdtTlg9icNFZC6wG3/nazlwKIXvn7eIUFYe7BCRChFu31UCdnrz9llNsLJwJQExNNAzjoiMAJqr6jURFs8B+oatXwvYV1X/8Irex90jDl0nA/cA/Ov4R5y6VHWcqnbwat8FJ+Bbb/5Y/J2vQut4+gIfJe4vSFnzcdWCg7zbo62AX7wi+6wmWFlJEpOAu0UkC0DcaHdTKfiAy3hEZDDQG7gwyirTgbEi0tBbPwP3UPb+kHWeA04VkRYhZeOBF6M8JCzv/JyvObirj26BFUTkPOBHVS2PV8R3AxNFpFVI2Q3Av1U1cCVhn9UES8oY14kgIscB1+H6Pq6Ea6X6XNFblT8iUh3XPfV83Bjj4Qap6iYR6YBr6VoBqAI8p6r3he3rYNx/xupAZdwvtpu9OujlnojM9a4oAu+LPV9eFzRTgQZAJu7++1Xquqgpd0SkPXAX7v90NVybkRtUdU/IOvZZTaAykySMMcbEX1m53WSMMSYBLEkYY4yJypKEMcaYqCxJGGOMicqShDHGmKgsSRhjjInKkoQxxpioLEkYk2QiUkNEVESaJTsWY8JZkjBpR0Q6i8gXcdrX12Hv47bvkhKRoSLycjJjKIqI3Ccilyc7DlM6rMW1KbdE5DDgY1Wtm+Q4agBbgANVdUUyYzEmnF1JmIQRkbkicqaIvCEiO0RkZegvUBFpJiK/iciRIrJARHaKSAsRqSAit3jr7xSRL0Wke8h2/UVkRcj7CiJym4isEZE/ReR5EakbsryOiDwtIn+JyDZvGMxuwBdAHe9Wj4pIgyj7LiqWbBH5TETOFpGfRWS7iHwqIkcUcV6OEpEvRGSXiHyHG4sidPlor4vx0PN4joi8453Hn0RksIgcJyJLvWP+x+v5OLBNFRGZ7v3Na7z56iHL80Skh4j81/u7lovI0JDlh4rIuyKyWUQ2ishHInKQt+wlERkf9u/4unduN4rIk+J6Yg0sf1xErhCRR739rRORhyXyGA8mxViSMIl2D/AgUA8YCdwqIn8PWV7fW34hUAc3xvMYYCium+j9gBnAG+L19BnBJUA3oBNwIG5YysdClj8OVMWNmdwIGKuq84ABwAZVFW/6g8L8xHI0MATohRtTfTnwSKRAvS/GV4GPcZ34DQIuj/J3hboT1ytqXW/fTwI3Aid6f9N6CvZ6/E/c8J8tgTbeOhNDQwEewJ332rjBp54MSa4PA//1YmyO61n11wh/TwXgNWAFbizpNkBNCp5/gNuBz4DGQE9c190jffzdJtlU1SabEjIBc4GHwsomAS95881wYxSfGLbOauD4sLKXgVu9+f7AipBlvwFHhbyv4+23Jm4Qn51A7QjxZQPrw8rC911cLNnesQ4IWd4G18NuRoRj9vf2WTGkrLe3j2be+9HA3GjnEdeTqQInhZSdEIgb11vq9tC/GWgX+rd6258dFtsGoIc3/xFwW5R/15eA8d58X2BV2N9T1zvnB3vvHwc+DNvHZODRZH9GbSp+sisJk2jLw95/jRslLNS8wIx3m6IBrovs8O1ah+9cRGrifp0uCNw2wv2qBjdaXEvgV1X9M9bAY4hlu6qG/sreghv1sWqE3R4MLNKQrq4pPGJaJD8GZjR/FLbQ0eq2AYERGpt7x94Qck6+wt1aqxayzXdhx9gCZHnzFwNDvduAI0SkcpS4WgKLQ/8eVV2Pu+o4zOexTAqzJGESLfy+c6TPXOhYCYH1/Q5Sn+G9Ntb820aBaal3vJLWzvAby44Y9ikUjme3j+12RijbE6EM3DnZDVSIcE5Cz3XUMSpU9Vvcl/wtwHBgobixLsJl4u/fqlyOh1EWWJIwidYy7P3RuOcO0WzA3b5oF1bejgi/uFV1I+72Taco+1sONBWR2hGW5VL0/4GYYvHpJ6C1dy8/4NAS7iuaH3Bf3G33ZieqmqeqrwLHAGuA0yOstghoE/oQWkTqAE0p+TkyKcSShEm0k0RkoLgGY/1wDytnRFtZ3Q3rKcBdItJWRLLEjcf9d2BalM0mApNEpJuIVBWRA0XkDG9/S3GD3T8kIo1EpKaIZHvbrQD2EZGuIlJJRCrGIZbivO293ioi+4hIc9xD5rjxbkc9AMwUV3Osioi0FJEhfvfh1dY62LvN1AJ3S29FhFU/BNYCU0RkPxFpDMwE3lLVon4MmDRhScIk2j242kfrgVnAjar6fjHbTMLVjnkNWAecDwxQ1ZVR1p+K+2J6DtiEe8bROWT5mbhf1t8DfwDXAqjqb8DNwBteeeM4xFIkVd2NqwnVF/fl+j7wKPnPUeLlWtyY2e8Am4F3gVZFblFQO+ATb9u3cQ+ZXwlfyXsW0Q9XS+1H4FtgK3Du3gRvUoc1pjMJ49X1f1xVH4/zfvsDD6hqs3ju1xhTmF1JmHQkRH9oa4yJo4rFr2JMahCRKrjPbBfglySHY0y5YEnCpJMjcQ3LVgJnJTcUY8oHeyZhjDEmKnsmYYwxJipLEsYYY6KyJGGMMSYqSxLGGGOisiRhjDEmKksSxhhjovp/yNtMrRLCVxIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Compute fraction of overall variance lost when projecting to k coordinate directions\n",
    "total_coordinate_variance = np.cumsum(coordinate_variances)\n",
    "total_coordinate_variance = total_coordinate_variance/total_coordinate_variance[783]\n",
    "\n",
    "# Compute fraction of overall variance lost when projecting to k eigenvector directions\n",
    "total_eigenvector_variance = np.cumsum(eigenvector_variances)\n",
    "total_eigenvector_variance = total_eigenvector_variance/total_eigenvector_variance[783]\n",
    "\n",
    "# Plot these results\n",
    "plt.plot(np.arange(1,784), total_coordinate_variance[784:0:-1], 'b-', lw=2)\n",
    "plt.plot(np.arange(1,784), total_eigenvector_variance[784:0:-1], 'r-', lw=2)\n",
    "plt.xlabel('projection dimension', fontsize=14)\n",
    "plt.ylabel('fraction of residual variance', fontsize=14)\n",
    "plt.xlim(0,784)\n",
    "plt.ylim(0.0,1.0)\n",
    "plt.legend(['coordinate directions', 'PCA directions'], fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvalues, eigenvectors = np.linalg.eigh(Sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function that returns the project-and-reconstruct operations as a single matrix\n",
    "def projection_and_reconstruction(k):\n",
    "    U = eigenvectors[:,(784-k):784]\n",
    "    P = np.dot(U,U.T)\n",
    "    return P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_digit(x):\n",
    "    # Make sure all entries of x are in the range [0,255]\n",
    "    for i in range(784):\n",
    "        x[i] = max(0.0, x[i])\n",
    "        x[i] = min(255.0, x[i])\n",
    "    # Now display\n",
    "    plt.axis('off')\n",
    "    plt.imshow(x.reshape((28,28)), cmap=plt.cm.gray)\n",
    "    plt.show()\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_effect_of_PCA(x, k_list):\n",
    "    print(\"Original:\")\n",
    "    show_digit(x)\n",
    "    for k in k_list:\n",
    "        if (k > 0) and (k < 784):\n",
    "            print(\"Projection to dimension \", k)\n",
    "            P = projection_and_reconstruction(k)\n",
    "            show_digit(P.dot(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original:\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAGmklEQVR4nO3dTYjNbQPH8TnmHsWsDE1q7IwpyUgzLCzZzMqGSEp5WUqapCzMBmElZTbKhmS8JMrKysIslJSyMEI2CnlZmBKZzrN6nnrq/K+TM4zfmfvzWfp13eef+t7/cnVmavV6vQPIs+BvPwDQmDghlDghlDghlDgh1D+lsVar+adc+MPq9Xqt0Z97c0IocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUIocUKo4o/G5M84cuRI5Xbs2LHi2SVLlszqs0+fPl3cz507V7l9/PhxVp/Nr/HmhFDihFDihFDihFDihFDihFDihFC1er36t/z5FYCtGR4eLu6Tk5OVW1dXV/Fss7vGnp6e4r5gQfn/x7dv367cdu7cWTw7MzNT3GnMrwCENiNOCCVOCCVOCCVOCCVOCCVOCOWeswXN7hKnpqaK+9KlSyu3M2fOFM+Wvm/Z0dHRsW/fvuJ+9OjR4l76vujY2Fjx7MmTJ4s7jbnnhDYjTgglTgglTgglTgglTgglTgjlnrMF27ZtK+43b94s7m/evKnchoaGime/fPlS3JvZvXt3cb9y5Url9vnz5+LZkZGR4v748ePi/m/lnhPajDghlDghlDghlDghlDghlDghlN/P2YLBwcFZnb948WLlNtt7zGYmJiaK+/r16yu30dHR4tm7d+8W976+vuLO//PmhFDihFDihFDihFDihFDihFC+MtaC58+fF/eBgYGW95cvX7b0TL9Lb29v5Xb//v3i2TVr1hT3/fv3F/fLly8X9/nKV8agzYgTQokTQokTQokTQokTQokTQrnnbMFs7zm7u7srt2/fvrX0THNh+/btxf3GjRvFfWZmprivWLGicnv//n3xbDtzzwltRpwQSpwQSpwQSpwQSpwQSpwQyo/GbKCrq6u4d3Z2ztGTZHn48GFxf/fuXXFfvnx5ca/VGl73/Wt5c0IocUIocUIocUIocUIocUIocUIo95wNbNiwobivXLlyjp4kS7N7zBMnThT38fHx4r5ly5bK7erVq8Wz85E3J4QSJ4QSJ4QSJ4QSJ4QSJ4TyozEb2LRpU3Fv9tWpZtr1R2M2s3jx4uI+PT1d3J89e1a5DQ4OtvRM7cCPxoQ2I04IJU4IJU4IJU4IJU4IJU4I5Stjf8HWrVsrt+vXr8/hk2Tp7++v3NatW1c8+/Tp09/9OH+dNyeEEieEEieEEieEEieEEieEEieEcs/ZwNTUVHF/+/Ztce/r6yvuPT09v/xM7WBmZqa4z+bvbdGiRS09Uzvz5oRQ4oRQ4oRQ4oRQ4oRQ4oRQ4oRQ7jkb+PTpU3G/dOlScR8bGyvuw8PDv/xM7eD79+/FfXJysrjv2LHjdz5O2/PmhFDihFDihFDihFDihFDihFDihFDuOVtw69at4t7snrN0n3fq1Kni2devXxd35g9vTgglTgglTgglTgglTgglTgjlKqUFzb5S1kx3d3fltnfv3uLZ48ePz+qzaR/enBBKnBBKnBBKnBBKnBBKnBBKnBDKPWcLvn79WtxfvHhR3AcGBiq3w4cPF8/eu3evuD969Ki4/0kLFy4s7mvXrp2jJ5kfvDkhlDghlDghlDghlDghlDghlDghlHvOFkxPTxf3CxcuFPfz589XbqXvenZ0NL/nHBkZKe5Pnjwp7vV6vbiXNLvnXL169R/77PnImxNCiRNCiRNCiRNCiRNCiRNCiRNC1Up3S7VazcXTH7Br167KbWhoqHh2dHR0Vp99586d4j4+Pl65vXr1qnh2YmKiuG/cuLG4P3jwoHLbvHlz8Ww7q9frtUZ/7s0JocQJocQJocQJocQJocQJocQJodxzhuns7Czue/bsKe5nz54t7suWLfvlZ/qvDx8+FPfe3t6W/9sdHR0dq1atqtya3bG2M/ec0GbECaHECaHECaHECaHECaFcpcwz/f39xf3QoUPF/eDBg5VbrdbwX/z/5+fPn8X9wIEDxf3atWuV248fP4pn25mrFGgz4oRQ4oRQ4oRQ4oRQ4oRQ4oRQ7jnhL3PPCW1GnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBBKnBCqVq/X//YzAA14c0IocUIocUIocUIocUIocUKo/wA83zYgdWwCEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection to dimension  100\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJOElEQVR4nO3dTW+NXRvG8VVtb7vVVipShAoS7UiDeilq6i0dEAMDifgC5hIfoZ/ARGIgQULSIIIO0DRhgBAVFUEkJK0qTVXVRu/ZM3IdJ73sex/b8/8NHVntVd2HK3FmrVU1OzubAPiZV+4HAPBzlBMwRTkBU5QTMEU5AVM1KqyqquK/coESm52drfrZn/PmBExRTsAU5QRMUU7AFOUETFFOwBTlBEzJOSf8VFdXyzzaZfTjx48/+TgoId6cgCnKCZiinIApygmYopyAKcoJmKKcgCnmnGVQKBQys4aGBrm2vr5e5t+/f5f56OiozIvFoszx3+HNCZiinIApygmYopyAKcoJmKKcgKmKHaXMm6f/Xamq+ulpgymleNyQV3t7u8x7enoys+7ubrl26dKlMp+YmJD5rVu3ZH7x4sXMbHh4WK7Fn8WbEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzBlO+eM5pi1tbVzzqempuTa6HjJSFdXl8z37duXmUVzzGhGu3btWpmvX79e5suWLcvMTp48KdcODQ3JHL+HNydginICpignYIpyAqYoJ2CKcgKmKCdgynbOGV1VF80i1Zw07xxz06ZNMt+1a5fM1Qz2/Pnzcm20H7OtrU3mhw8flvnWrVszs2g/58jIiMzHxsZkrq43jK4+jOa/pd7DWwq8OQFTlBMwRTkBU5QTMEU5AVOUEzBFOQFTtnPOSN45aB47duyQeUdHh8yfPn2amfX398u1g4ODMr9586bMp6enZX7s2LHMTJ23m1JKL168kPnVq1dlrmaR0f7ef/75R+Zfv36d8/cuF96cgCnKCZiinIApygmYopyAKcoJmKKcgKmKnXNGc68oV9Tdniml1NraKvNoZnbv3r3M7OXLl3JtXleuXJF5Z2dnZnbw4EG59ujRozKP9oOqOWmxWJRrozln3v2g5cCbEzBFOQFTlBMwRTkBU5QTMEU5AVMVO0qJrgDMs2WssbFR5kuWLJH5+/fvZX737t3M7O3bt3JtZP78+TIfHx+X+alTpzKz1atXy7Xd3d0yj0Yxvb29Mle+ffsm85oa/VGPtpSVA29OwBTlBExRTsAU5QRMUU7AFOUETFFOwFTFzjmjbV15tgC1tLTIfPny5XP+2imlNDk5mWu9Es1gX79+LfNHjx5lZmfOnJFr29vbZb5z506Z9/X1ZWbPnj2Ta2dmZmTe1NQk82iLYXQUaynw5gRMUU7AFOUETFFOwBTlBExRTsAU5QRMVeycM9qvmWd/XjQTW7Bggcyja/ai9XlE875ov6dar470TCmloaEhmUfzYzUnjeackejz4rjfkzcnYIpyAqYoJ2CKcgKmKCdginICpignYKpi55zR/ro8+zmjvX3RdXMjIyMyj861zSP63vX19XP+2hMTEzK/f/++zPfs2SPztra2336mXxV9HqJzkJlzAvgfygmYopyAKcoJmKKcgCnKCZj6a0cpatzx5csXubahoUHm0dGY0batYrEo81LKM2IaGxuT+cDAgMx3794t840bN2Zm0ZGf0QgpuiKwUCjIvBx4cwKmKCdginICpignYIpyAqYoJ2CKcgKmKnbOGc2tGhsbM7Nozvnp0yeZR8dLrlmzRubNzc0yV6JjNaMZa0Rtl4tmy9HRmNG1jevWrcvMWltb5dpozhnNluvq6mReDrw5AVOUEzBFOQFTlBMwRTkBU5QTMEU5AVMVO+fMc9RhNG97/vy5zG/fvi3zrq4umauZ3Z07d+TaaL6bN1+0aFFmNj4+LtdGR4ZGvzM1/80zG04pvgIwmuGWA29OwBTlBExRTsAU5QRMUU7AFOUETFFOwFTFzjkjas9mNI/7+PGjzC9duiRztS8xpZQ6Ojoysxs3bsi1pb6KTn39aB9rtGfy1atXMm9pacnMampK+1HNc55vqfDmBExRTsAU5QRMUU7AFOUETFFOwBTlBEz9tXPOycnJzCyac0aePHki82gWuX379szs3Llzcu3jx49lrvaxphTPIqenpzOzaBYYnZkbzTlXrVqVmUVnDefFfk4Av4xyAqYoJ2CKcgKmKCdginICpko6SlFHUEZHFZZSNE6IRFcERvmWLVsys23btsm10Sgl+tmiKwTVVXhjY2NybWThwoUyr6+vz/X182DLGIBfRjkBU5QTMEU5AVOUEzBFOQFTlBMwVbY557x5+t+F6Jq+aL3atpV3xhodnfnw4UOZb9iwITM7cuSIXBttVxscHJT51NSUzPNYvHixzDs7O2WurkYsFApzeqZKxpsTMEU5AVOUEzBFOQFTlBMwRTkBU5QTMFXSOaeaJ0ZzzOrqaplH19Gp4y+jWV80Bx0dHZV5X1+fzNVVd9Es8MSJEzI/ffq0zK9duybzaIarrFixQubq504ppXfv3s0p+xXRXNxR5T0x8H+CcgKmKCdginICpignYIpyAqYoJ2CqbHPOaJYYXckWnTNayrnW58+fZT4wMCDzDx8+ZGb79++Xa/fu3Svz48ePy3z9+vUyv379emamrlVMKaUDBw7IPPqd9ff3Z2bDw8NybST6PJXzHOUsvDkBU5QTMEU5AVOUEzBFOQFTlBMwRTkBU1XBnsuyDX/ynlurfq5o5pVX9Ozq2ZqamuTaQ4cO5cpXrlwpc/V309jYKNdG+YMHD2Te29ubmV2+fFmujUSfl1J/JpTZ2dmffmB4cwKmKCdginICpignYIpyAqYoJ2DKdpRSyWpra2VeLBZL9r27u7tl3tPTI/PNmzdnZtHRlm/evJH52bNnZX7hwoXMLNquFskz3io1RilAhaGcgCnKCZiinIApygmYopyAKcoJmGLOWQLR9YXqesLp6ek//Ti/pbm5OTMrFApybXRk6MTExJye6W/HnBOoMJQTMEU5AVOUEzBFOQFTlBMwRTkBU8w5S0DNMVNKqa6uLjOL9oLOzMzIPO++R/z3mHMCFYZyAqYoJ2CKcgKmKCdginICpignYKqm3A/wN4rOpVWz5Zoa/SuJzl8t55m5+LN4cwKmKCdginICpignYIpyAqYoJ2CKLWNAmbFlDKgwlBMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwRTkBU3I/J4Dy4c0JmKKcgCnKCZiinIApygmYopyAqX8BLTFtSmjMAkQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection to dimension  50\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJJklEQVR4nO3dy0uVXR/G8XtnntJSI/GARhGJaEV0QLImVkSDyFHQqL+tmdGgwFFUEJ2UICQCK8hS1IxKPGWZmvrOXnh5va+fj/ez29eO72fYxUrddLXAH2ut3MbGRgLAz45CfwMANkc5AVOUEzBFOQFTlBMwtVOFuVyOX+UCebaxsZHb7M/ZOQFTlBMwRTkBU5QTMEU5AVOUEzBFOQFTcs4JP7ncpiOx/+KU0d+DnRMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwxZyzAMrLy1Ozqqqqba9NkiRZXV2V+ezsrMzX1tZkjj+HnRMwRTkBU5QTMEU5AVOUEzBFOQFTjFK2YedO/bG1t7fLvLu7OzU7deqUXNvU1CTz6elpmT9//lzmDx8+TM3GxsbkWvy72DkBU5QTMEU5AVOUEzBFOQFTlBMwRTkBU7ZzTucrIBsaGmR++vRpmZ89ezY1O3z4sFzb3Nws8+rqapkfOXJE5mqO2tfXJ9eOjIzIHP8MOydginICpignYIpyAqYoJ2CKcgKmKCdgqmBzzh079P8LZWVlMl9fX5f5ysrKP/6eturo0aMy7+jokPni4mJq1t/fL9cuLS3JPDpLevLkSZlfuHAhNVtYWJBrb926JfMvX77IHP+LnRMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwVbA5Z3T3a2VlpczzOedsaWmReVdXl8w7OztlPjw8nJrdv39frh0aGpJ5ND++ceOGzK9fv56a9fT0yLVTU1Myv3v3rsx///4tc8X5/O92sXMCpignYIpyAqYoJ2CKcgKmKCdginICpvI651Szp6znOZeXl7f1PW1FbW2tzFtbW2UezWjn5uZSs5mZGbk2Es1/b968KXN1b+7Vq1fl2itXrsh8dHRU5i9fvpS5Es3No89lbW1t2187X9g5AVOUEzBFOQFTlBMwRTkBU5QTMJXXUUqWYzrR2uiKyCwqKipkXlpaKvPJyUmZq6fyCv0r/WfPnqVmly5dkmu7u7tl/vr1a5mrz0WNn5Ik/vdSUlIi80J/7pth5wRMUU7AFOUETFFOwBTlBExRTsAU5QRMFexqzGiuFF2TmM8jY7t375Z5eXm5zCcmJmSursaMZqR79+6VeTTvm52dlfng4GBq9uTJE7m2t7dX5sePH5f5gQMHUrNXr17JtdG/l+hImSN2TsAU5QRMUU7AFOUETFFOwBTlBExRTsBUwYY/0VWFq6urf+g7+X/RLLGhoUHmCwsLMl9cXEzNojlldHVmR0eHzKPPXZ2b7Ovrk2v3798v85qaGpkfPHgwNYvmnJEszwsWCjsnYIpyAqYoJ2CKcgKmKCdginICpignYMp2zlnIe0R37dol85aWFplHZwfVPG9sbEyujUxPT8t8z549Mp+fn0/Noif6BgYGZH758mWZHzp0KDWLZqTq+y5W7JyAKcoJmKKcgCnKCZiinIApygmYKtgoJToaFY1aousps1ydqY50bUU0amlra0vNHj16lOlrf/36Veb19fWZ/n7l06dPMo+uHFWfS1NTk1zLKAXAH0M5AVOUEzBFOQFTlBMwRTkBU5QTMGX7Llo0By0rK5N5ljnn+Pi4zEdHR2Xe2Ngoc/W9V1RUyLW/fv2SeSSXy2Var/z48UPmpaWlMlezzOgzfffuncyLETsnYIpyAqYoJ2CKcgKmKCdginICpignYKpgc85o3hbNOXfs0P+vqPOe0Qz08+fPMn///v22v3aS6Ofoop87UlJSIvPa2lqZf//+PTVbWlqSa799+ybz6KypmoNGZ0H/RuycgCnKCZiinIApygmYopyAKcoJmKKcgCnb85yRaE6qzkxGc04160uS+Jm96upqmatzi9GZxyznVJMkPger5qDRnHNiYkLmHz9+lHlra2tqFj3LGP1cKysrMnfEzgmYopyAKcoJmKKcgCnKCZiinIApygmYyuucU80iozll1vtVo3ONSnQ37NzcnMyjtyR7enpSs3v37sm1AwMDMl9bW5P5zMyMzLPcixvNQaN7bdVZ1uiMbJQz5wTwr6GcgCnKCZiinIApygmYopyAKdsjY+vr6zJX10smSbZRzOrqqsxHRkZk/vPnT5mfO3cuNbt27ZpcOzw8LPP5+XmZT01NyTyL6KhcXV2dzKuqqlKz6GnELKMzV+ycgCnKCZiinIApygmYopyAKcoJmKKcgKm8zjmzPGcXHX2KZpFZj5wpg4ODMn/w4IHMm5ubU7Pz58/LtUNDQzK/c+eOzKNjW1ns27dP5urnThI9J42efIzm4sWInRMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwVbDznNEMNMqjOag63xfNQKOvHT0BePv2bZk3NjamZm1tbXKtulYzSeKzpC9evJD55OSkzJVjx47JvL6+XuaLi4up2ezsrFyb5UpPV+ycgCnKCZiinIApygmYopyAKcoJmKKcgCnbe2uzUnNOdT9qksRP2UUz1jdv3si8v78/Nbt48aJcG52J7O3tlfmJEydkruacZWVlcu2ZM2dkHs2X3759m5p9+PBBri3GJ/4i7JyAKcoJmKKcgCnKCZiinIApygmYopyAqaKdc0b3lKpZZHQHamVlpcyjt0HVucQkSZKnT5+mZnNzc3JtV1eXzDs7O2UenblUs8zonGt0pjKa/6rPJZpz/o3YOQFTlBMwRTkBU5QTMEU5AVOUEzCVU78ez+Vy23/Dr8DUuEQdJ9uK6MhYPp+jq6urk3k0aomu1uzo6EjNohHUxMSEzKOnEx8/fpyajY2NybXFbGNjY9OzdOycgCnKCZiinIApygmYopyAKcoJmKKcgKm/ds6J7Wlvb0/NWlpaMv3d0fOC4+PjqVn0tGExY84JFBnKCZiinIApygmYopyAKcoJmKKcgCnmnNiympoamZeXl8s8ulJUzTKXl5fl2ujaTmfMOYEiQzkBU5QTMEU5AVOUEzBFOQFTlBMwxZwTW1ZaWirz6F7bLM825vMu4EJjzgkUGcoJmKKcgCnKCZiinIApygmYYpQCFBijFKDIUE7AFOUETFFOwBTlBExRTsAU5QRMUU7AFOUETFFOwBTlBExRTsAU5QRMUU7AFOUETMnznAAKh50TMEU5AVOUEzBFOQFTlBMwRTkBU/8BnuVf0J62MooAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection to dimension  25\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAIVUlEQVR4nO3dS0+TXRvF8bsUqmBaFZBTIUJUjJgYE40MjImJEz6AX9cJJko8DlCi4llQpIrIQRR6eGZP8iTc68pL39pV+v8NWdlSleVOvLL3ztRqtQSAn45mfwAA+6OcgCnKCZiinIApygmY6lRhJpPhv3KBBqvVapn9vs7OCZiinIApygmYopyAKcoJmKKcgCnKCZiinIApygmYopyAKcoJmKKcgCnKCZiinIApygmYkuc58fflcjmZHzlyROaZzL5HA/+1s7Mj8729PZnj72HnBExRTsAU5QRMUU7AFOUETFFOwFTLjlKikUFHR/q/O9G4Ip/Py7y/v1/mxWJR5qdPnz7w2kKhIPOtrS2Zv3jxQubPnz9Pzd68eSPXbm9vyxz/G3ZOwBTlBExRTsAU5QRMUU7AFOUETFFOwJTtnLOeOWaSJElnZ/pvraenR64dHh6W+dTUVF35uXPnUrORkRG5tq+vT+bZbFbmS0tLMr9//35qNjs7K9fOz8/LfHV1VebValXm7YadEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzDVtDlnvXPMWq0m83K5nJpF1z9WKhWZR6L1Hz58SM0+fvwo1548eVLmZ8+elfno6KjMb968mZp1dXXJtdHf6dOnT2VeKpVk3m7YOQFTlBMwRTkBU5QTMEU5AVOUEzBFOQFTtuc5ozlmPWf/ortdd3d3ZR6dmVQz1iRJknfv3qVmCwsLcu3a2prMz58/L/OZmRmZT05OHvjX/v79u8y/fv0qc+ac/8XOCZiinIApygmYopyAKcoJmKKcgCnKCZhq2TlnPaIZab3vc0brV1ZWUrPoDczl5WWZP3v2TOYbGxsyv337dmqm3hVNEn0fb5IkyevXr2W+uLiYmv369UuuPYzYOQFTlBMwRTkBU5QTMEU5AVOUEzDVlqOU6AnAgYEBmQ8ODsr8xIkTMu/t7U3Noif+olFK5M6dOzJXV2cWCgW5tlgsyvzSpUsyV8flXr58KdceRuycgCnKCZiinIApygmYopyAKcoJmKKcgCnbOWcjHT9+XObRHDOa93V26j9WNcucmpqSa6NjV9HRqui43L1791KzsbExuXZ6elrmFy5ckLmagzLnBGCDcgKmKCdginICpignYIpyAqYoJ2CqaXPORp7XjHR3d8v82LFjMs9kMjLf3t6WuXpicGhoSK6N5qCPHj2SeUTNEx88eCDXRldjRnPSy5cvp2Zzc3Ny7adPn2Teitg5AVOUEzBFOQFTlBMwRTkBU5QTMEU5AVNteZ4zOtPY0aH/zYpmtOvr6zJXd89G5zGjs6KNFM0ab9y4IfNoRnvmzJnUbGJiQq5lzgngr6GcgCnKCZiinIApygmYopyAqbYcpezs7Mi8XC7L/OjRo3Xl6kjZysqKXLu0tCTzRopGRKVSSebZbFbmp06dSs2iZxkPI3ZOwBTlBExRTsAU5QRMUU7AFOUETFFOwFRbzjm3trZkXqlUZN7f319Xrr5/NEuMjpQ10+/fv2UezTnVlaT5fP5An6mVsXMCpignYIpyAqYoJ2CKcgKmKCdginICppo254ye0YvU84Rg9ETf5uamzAuFgszHx8dlrj776uqqXNvT0yPztbU1mTdTND9WV5Lmcjm5NpqhRt/bETsnYIpyAqYoJ2CKcgKmKCdginICpignYKot55yRaFYY3Ws7PDwsc3WeM7qfNZpzNlN0H3A0P46eXlSipxGZcwL4v6GcgCnKCZiinIApygmYopyAKcoJmGronFPNMuudczZS9Ebmly9fZF6tVmWuZpnRnbfd3d0yb6Zozrm7uytz9TMRzZaj85ytiJ0TMEU5AVOUEzBFOQFTlBMwRTkBU7ZPADZz1PL582eZLywsyHx5eVnmfX19qdnQ0JBc29vbK/NmqufIV5LoK0uj60yj8VUrYucETFFOwBTlBExRTsAU5QRMUU7AFOUETDVtzhldbdnIqy8jpVJJ5k+ePJH5lStXZD49PZ2aRddqXrx4UeaPHz+W+cbGhszrEc1ou7q6ZL6+vp6aRdeV/vnzR+atiJ0TMEU5AVOUEzBFOQFTlBMwRTkBU5QTMNXQOWc9s8pmzjkj8/PzMn/48KHMx8bGUrNisSjX3rp1S+bRuce7d+/KXD1PGM1vozx6vvDHjx+pWTTndP55OSh2TsAU5QRMUU7AFOUETFFOwBTlBExRTsCU7XlOZ9ETgLOzszIfHR1NzaI55tWrV2U+MjIi85mZGZmrZ/rGx8fl2ij/9u3bgfPNzU259jBi5wRMUU7AFOUETFFOwBTlBExRTsAU5QRM2b7P2crm5uZkrt6xjOZ5169fl/nExITMJycnZa7OXOZyObk2uu93cXFR5q9evUrNfv78KdceRuycgCnKCZiinIApygmYopyAKcoJmMqoo1uZTKZ1z3W1qHw+L/Nr167JPDpSFj0hODg4mJqp42RJkiRv376VeTRiUleKvn//Xq4tl8syd1ar1TL7fZ2dEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzDFnLPNqOcHkyRJBgYGUrNsNivXRs8PRkfK1DN/rTzHjDDnBFoM5QRMUU7AFOUETFFOwBTlBExRTsAUc078h7oas7NT36QaPeu4t7d34LxSqci1rYw5J9BiKCdginICpignYIpyAqYoJ2CKcgKmeAKwzWQy+47U/lXPrDGac1ar1QOvjz539L1bETsnYIpyAqYoJ2CKcgKmKCdginICpjgyBjQZR8aAFkM5AVOUEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzBFOQFTlBMwRTkBU5QTMEU5AVOUEzAlz3MCaB52TsAU5QRMUU7AFOUETFFOwBTlBEz9A8+NCmUi8dNmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Projection to dimension  10\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAOcAAADnCAYAAADl9EEgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAI+ElEQVR4nO3dS09UaxuE4dVyRkVARESEKAokGiEmJsaRU/+Bv9SRI6cGMSZOtqCoiMhBjspJ6D3byZewqgj9Ear1voa78mqD1F4JT553VarVagEgz7mz/gAAjkY5gVCUEwhFOYFQlBMI1ajCSqXCr3KBU1atVitH/XeenEAoygmEopxAKMoJhKKcQCjKCYSinEAoygmEopxAKMoJhKKcQCjKCYSinEAoygmEopxAKMoJhKKcQCjKCYSinEAoygmEopxAKMoJhJJXYyarVI68TfA/jY3lX1pra6s8e+nSJZl3d3fXlKs/v7m5WZ7d3d2V+fLysswXFhZkvrS0VJptbm7Ks/j/4skJhKKcQCjKCYSinEAoygmEopxAKMoJhIqdc547p/+/4eaBnZ2dpVl/f788OzQ0VFPe19cnczUHbWtrk2d///4t85WVFZl//vxZ5jMzM6XZ+/fv5dm5uTmZ//r1S+b4Xzw5gVCUEwhFOYFQlBMIRTmBUJQTCEU5gVCVarVaHlYq5WGtf7HZx2xqapJ5R0eHzAcGBkqzsbExeXZ4eFjmvb29Mnf7omqGW+uuaa1z0k+fPpVmU1NT8qzL1Qy1KIpidXW1NDs8PJRn61m1Wj2yDDw5gVCUEwhFOYFQlBMIRTmBUJQTCEU5gVCx+5wNDQ0yb29vl7na5+zq6pJnL168KHM1Gy6KolhcXJT59+/fS7OdnR151n320dFRmY+MjMh8YmKiNHOz5ZaWFpm7HV01B1Uz0KLw89t6xJMTCEU5gVCUEwhFOYFQlBMIRTmBUJQTCBU756zl/ZvHyZX9/X2ZHxwcyNzdDTs5OVmaTU9Py7Pu+/Lw4UOZP336VOYPHjwoza5fvy7P3r9/X+bb29syV993N8dcX1+XeT3ug/LkBEJRTiAU5QRCUU4gFOUEQlFOIFTsKMVxvxpXq1duLWtvb0/mbvVpa2tL5mplbGNjQ551Xrx4IXM3ilH5rVu35Fm1plcU/spRtRa2ubkpz+7u7sq8Hl8/yJMTCEU5gVCUEwhFOYFQlBMIRTmBUJQTCHVmc053vaRby3JrXWqW6VaX3Mysubm5pvzChQsyP02vXr2SeU9PT2nm5r/qbFH4K0fVStrCwoI8667OdP/m7ufxLPDkBEJRTiAU5QRCUU4gFOUEQlFOIBTlBELF7nO6Oaebuak5p5tjuhmqe9Wde1XejRs3SrPZ2Vl51n12Z21tTeYfP34szc6fPy/Pus/m5pzq++ZmqPPz8zJ3V2e6n6ezwJMTCEU5gVCUEwhFOYFQlBMIRTmBUJQTCFW3c043U1NzK/c6ObeP6eaY7e3tMm9oaCjN3Lzt5cuXMneamppkru4DXllZkWfd/NdRd+a6Gav7nruvmzkngGOjnEAoygmEopxAKMoJhKKcQKjYUYq7qtBddahyN0pxV1eqla+i8KOY/v7+0sytVblX2U1PT8u8t7dX5teuXSvNGhv1j4v7bO71hurViu7VhW6M40YpiXhyAqEoJxCKcgKhKCcQinICoSgnEIpyAqFi55yOu75ya2vrxGfdetLQ0FBN59X1lG6+62asbtbY2dkpczcHVdS6WVH4z6a4fzNHzVCPk7uv7TTw5ARCUU4gFOUEQlFOIBTlBEJRTiAU5QRC1e2c01G7g+41eG7f080Kr1y5InM1M3NzTvVqw6Ioiq9fv8rczRoHBwdLs8uXL8uz7vtWS+6uQnXfNzfHdPuiZ4EnJxCKcgKhKCcQinICoSgnEIpyAqEoJxDqj51zqlcILi4uyrMud7uF7u5Zdcfq+vq6PPvt2zeZu/MuV8bHx2Xe3d0t81rujnX7lC53c0w3B3WvpDwNPDmBUJQTCEU5gVCUEwhFOYFQlBMIRTmBUH/snFOZn5+X+T///CPzL1++yHx4eFjm7v2fyubm5onPHsfs7Gxp5nYqJyYmZO7u3FXzXzcjbWhoqClnnxPAsVFOIBTlBEJRTiAU5QRCUU4g1F85SnFXNL57907mr1+/lvno6KjM1dWafX198mxPT4/MZ2ZmZF4Lt67mrgR1X1tHR0dp1tzcLM+6UYtbCUtUf58Y+EtQTiAU5QRCUU4gFOUEQlFOIBTlBEL9lXNO58OHDzKfnJyUuVudevToUWl29+5defbJkycyf/Pmjczd2lct3PzYzSLb2tpKs9bWVnnWzUHrEU9OIBTlBEJRTiAU5QRCUU4gFOUEQlFOIBRzziNUq1WZu6sz3759K/OxsbHSbGRkRJ599uyZzN3e4vPnz2U+NzdXmnV1dcmz9+7dk/nAwIDM1azSXQla69WZ7vumrs50Py8nxZMTCEU5gVCUEwhFOYFQlBMIRTmBUJQTCMWc8wQWFxdl7nYq1b227l5atyva3d0t88ePH8tcvQLw8PBQnlX3zhaF37n88eNHaebmkO6zuTlnY6Ouwv7+fml2cHAgz54UT04gFOUEQlFOIBTlBEJRTiAU5QRCUU4gFHPOE1DzuKIoiqmpKZmr+1k3Njbk2fHxcZm7Oeft27dlfvXq1dLs58+f8uzW1pbMl5eXZb60tHSirCj8vqeaUxbF6e1k1oInJxCKcgKhKCcQinICoSgnEIpyAqEq6lfIlUol7/fLfwA17rh586Y8e+fOHZkPDg7K3F1vqVandnZ25Nm1tTWZu1U7dS2nG6Wsrq7K3I1a3NemRjG1jmGq1eqR927y5ARCUU4gFOUEQlFOIBTlBEJRTiAU5QRCMeesM+6KyM7OTpm7OadaZ3PXT7pZoVs5297eLs329vbkWXc9pcvd13aamHMCdYZyAqEoJxCKcgKhKCcQinICoSgnEIo55x+mUjlyZPaflpYWmTc1NZ3476511qhydzbxasvjYs4J1BnKCYSinEAoygmEopxAKMoJhKKcQCjmnH8YN+d0+6DqvPuz3U6km0WqvJ7nmA5zTqDOUE4gFOUEQlFOIBTlBEJRTiAUoxTgjDFKAeoM5QRCUU4gFOUEQlFOIBTlBEJRTiAU5QRCUU4gFOUEQlFOIBTlBEJRTiAU5QRCUU4glNznBHB2eHICoSgnEIpyAqEoJxCKcgKhKCcQ6l9mHXvM3iocXwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "index = 1000\n",
    "show_effect_of_PCA(train_data[index,], [100, 50, 25, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}