juliaL28demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## PCA (Part 2)\n",
    "\n",
    "CPSC 340: Machine Learning and Data Mining\n",
    "\n",
    "The University of British Columbia\n",
    "\n",
    "2018 Winter Term 1\n",
    "\n",
    "By Mike Gelbart (converted to Julia by Jason Chau)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [],
   "source": [
    "using Printf\n",
    "using PyCall\n",
    "pickle = pyimport(\"pickle\")\n",
    "plt = pyimport(\"matplotlib.pyplot\")\n",
    "using PyPlot\n",
    "using Statistics\n",
    "SkDecomposition = pyimport(\"sklearn.decomposition\")\n",
    "using LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAEhCAYAAADfxcKRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deXRU1Zo28CeEpDIQCkLmgRARZAggMgkyqwwCCspFxFYUri0KXJn6toDdgqJB7KvSgNgXGRWBvgqiKGCQSRbDJYwCMilDgIQxMxCSsL8/+FKdIme/lTokJ9PzWytrSb11qnadOmdnW6nnvB5KKQUiIiIii1Qr6wEQERFR1cLFBxEREVmKiw8iIiKyFBcfREREZCkuPoiIiMhSXHwQERGRpbj4ICIiIktx8UFERESW4uKDiIiILMXFRxn47//+b3h4eCAuLq6sh+K2zZs3w8PDA5s3by7roTicPn0affr0QWBgIDw8PDBmzJiyHhKRIQ8Pj2L9bN682XGuff3112U9bPr/3nrrLdStWxfVq1dHrVq1yno4FVr1sh5AVbRgwQIAwOHDh7Fr1y60a9eujEdUsY0dOxa7du3CggULEBYWhvDw8LIeEpGhHTt2OP373XffxaZNm7Bx40an25s0aYK9e/daOTRyYfXq1XjvvfcwefJk9O7dGzabrayHVKFx8WGxxMREHDhwAH369MEPP/yA+fPn39Pi48aNG/D19S3BEVY8hw4dQtu2bdG/f/+yHgqR6OGHH3b6d3BwMKpVq1bk9qqsvM5phw4dAgD85S9/QUhISBmPpuLjn10sNn/+fADA9OnT0aFDByxfvhzXr18v1rb16tVD3759sXLlSrRs2RI+Pj6YOnUqAGDOnDno3LkzQkJC4O/vj2bNmmHGjBnIzc11eoyuXbsiLi4Ou3fvRqdOneDn54f77rsP06dPx+3bt53ue/ToUfTq1Qt+fn4ICgrCiBEjkJmZaTi2BQsWoEWLFvDx8UFgYCAGDBiA3377zek+L730EmrUqIGjR4+iZ8+e8Pf3R3h4OKZPnw4A2LlzJzp27Ah/f380bNgQixcvFvdHwcfSJ0+exNq1ax0fWZ8+fRo3b97E+PHj8eCDD8JutyMwMBDt27fH6tWrizzO7du3MWvWLDz44IPw9fVFrVq18PDDD+O7775zut+KFSvQvn17+Pv7o0aNGujZsyf27dsnjpHoXuXm5mLy5MmIiIhAzZo18dhjj+HYsWNF7rdhwwY8+uijqFmzJvz8/PDII4/g559/dvn4xT1XWrZsiU6dOhXZPj8/H5GRkXj66acdt926dQvTpk1Do0aNYLPZEBwcjJdffhmXL1922rYk5jSlFN5//33ExMTAx8cHrVu3RkJCArp27YquXbs63TcjIwMTJkxAbGwsvL29ERkZiTFjxiA7O1vcR/Xq1cNbb70FAAgNDYWHhwemTJkC4M680KNHD4SHh8PX1xeNGzfGm2++afiYu3btQr9+/VCnTh34+Pigfv36Rf5MfOLECQwZMgQhISGw2Wxo3Lgx5syZI46vQlJkmevXryu73a7atGmjlFLq888/VwDUokWLirV9TEyMCg8PV/fdd59asGCB2rRpk/rnP/+plFJq7Nixau7cuWrdunVq48aN6uOPP1ZBQUHq5ZdfdnqMLl26qDp16qgGDRqozz77TCUkJKjXX39dAVCLFy923C8lJUWFhISoyMhItXDhQvXjjz+q559/XtWtW1cBUJs2bXLc9/3331cA1HPPPad++OEHtWTJEnXfffcpu92ujh8/7rjf0KFDlbe3t2rcuLGaOXOmSkhIUC+//LICoCZOnKgaNmyo5s+fr9avX6/69u2rAKjExETt/khPT1c7duxQYWFh6pFHHlE7duxQO3bsUDdv3lRpaWnqpZdeUl988YXauHGjWrdunZowYYKqVq2a0+tUSqkXXnhBeXh4qD//+c9q9erVau3ateq9995TM2fOdNznvffeUx4eHmrYsGFqzZo1auXKlap9+/bK399fHT58uFjvH9Hdhg4dqvz9/Q1rmzZtUgBUvXr11PPPP69++OEHtWzZMlW3bl3VoEEDlZeX57jvF198oTw8PFT//v3VypUr1ffff6/69u2rPD091YYNG8QxFPdcmTlzpgLgdE4rpdSPP/6oAKjvvvtOKaVUfn6+6tWrl/L391dTp05VCQkJ6vPPP1eRkZGqSZMm6vr1645tS2JOmzhxogKg/vVf/1WtW7dOzZs3T9WtW1eFh4erLl26OO6XnZ2tHnzwQRUUFKQ++ugjtWHDBjVz5kxlt9tV9+7d1e3bt7X7aO/evWr48OEKgFq3bp3asWOHSkpKUkop9e6776qPP/5Y/fDDD2rz5s3qs88+U7Gxsapbt25Oj7Fu3Trl5eWlmjdvrhYtWqQ2btyoFixYoAYPHuy4z+HDh5XdblfNmjVTS5YsUT/99JMaP368qlatmpoyZYr4PlY0XHxYaMmSJQqA+uyzz5RSSmVmZqoaNWqoTp06FWv7mJgY5enpqY4dOybeLz8/X+Xm5qolS5YoT09Pde3aNUetS5cuCoDatWuX0zZNmjRRPXv2dPz73//935WHh4fav3+/0/0ef/xxp8VHamqq8vX1VU888YTT/c6ePatsNpsaMmSI47ahQ4cqAOqbb75x3Jabm6uCg4MVALV3717H7VevXlWenp5q3LhxLvbKnf3Sp08f8T55eXkqNzdXDR8+XLVs2dJx+9atWxUANXnyZO22Z8+eVdWrV1ejR492uj0zM1OFhYWpQYMGuRwjkZHiLD7uPrf+93//VwFQO3bsUErd+aUaGBio+vXr53S//Px81aJFC9W2bVu3xqQ7V65cuaK8vb3VpEmTnO4/aNAgFRoaqnJzc5VSSi1btqzIea6UUrt371YA1Keffuq47V7ntGvXrimbzaaeffZZp/vv2LFDAXBafMTHx6tq1aqp3bt3O93366+/VgDUjz/+KI7h7bffVgDU5cuXtfe5ffu2ys3NVVu2bFEA1IEDBxy1+vXrq/r166sbN25ot+/Zs6eKiopS6enpTrePGjVK+fj4OM3lFR3/7GKh+fPnw9fXF4MHDwYA1KhRA3/605/wyy+/4MSJE8V6jObNm6Nhw4ZFbt+3bx+efPJJ1KlTB56envDy8sKLL76I/Px8HD9+3Om+YWFhaNu2bZHHPXPmjOPfmzZtQtOmTdGiRQun+w0ZMsTp3zt27MCNGzfw0ksvOd0eHR2N7t27F/nY18PDA0888YTj39WrV8f999+P8PBwtGzZ0nF7YGAgQkJCnMbkrn/84x945JFHUKNGDVSvXh1eXl6YP3++05+D1q5dCwAYOXKk9nHWr1+PvLw8vPjii8jLy3P8+Pj4oEuXLuUq+UOVz5NPPun07+bNmwOA49zYvn07rl27hqFDhzodn7dv30avXr2we/dul39WKM65UqdOHfTr1w+LFy92/Ik2NTUVq1evxosvvojq1e98hXDNmjWoVasW+vXr5zSeBx98EGFhYUXOl3uZ03bu3ImcnBwMGjTIaduHH34Y9erVc7ptzZo1iIuLw4MPPug0rp49e95Tgu+PP/7AkCFDEBYW5hhnly5dAMCx/44fP47ff/8dw4cPh4+Pj+Hj3Lx5Ez///DMGDBgAPz8/pzE+8cQTuHnzJnbu3GlqjOURFx8WOXnyJLZu3Yo+ffpAKYW0tDSkpaVh4MCBAP4vAeOKUZLj7Nmz6NSpE86fP4+ZM2fil19+we7dux1/J7xx44bT/evUqVPkMWw2m9P9rl69irCwsCL3u/u2q1evascVERHhqBfw8/MrcvJ5e3sjMDCwyPbe3t64efNmkduLY+XKlRg0aBAiIyPx5ZdfYseOHdi9ezeGDRvm9JiXL1+Gp6en4WstcPHiRQBAmzZt4OXl5fSzYsUKXLlyxdQYiYrj7vO1IGVRcL4WHJ8DBw4scnx+8MEHUErh2rVr2scv7rkCAMOGDcP58+eRkJAAAFi2bBlycnKc/ufj4sWLSEtLg7e3d5HxpKSkFDlf7mVOK5hfQkNDizzG3bddvHgRBw8eLDKmgIAAKKVMncdZWVno1KkTdu3ahWnTpmHz5s3YvXs3Vq5c6TTOgu+6REVFaR/r6tWryMvLw6xZs4qMseB/2CrTXMO0i0UWLFgApRS+/vprw9z+4sWLMW3aNHh6eoqP4+HhUeS2b7/9FtnZ2Vi5ciViYmIct+/fv9/0eOvUqYOUlJQit999W8HEmJycXOS+Fy5cQFBQkOkx3Isvv/wSsbGxWLFihdM+y8nJcbpfcHAw8vPzkZKSoo3oFryGr7/+2mn/EpUHBcfnrFmztKkZo1/OBYp7rgBAz549ERERgYULF6Jnz55YuHAh2rVrhyZNmjiNp06dOli3bp3h8wUEBDj9+17mtIL5p2ABVlhKSorTpx9BQUHw9fXV/o+emblq48aNuHDhAjZv3uz4tAMA0tLSnO4XHBwMADh37pz2sWrXrg1PT0+88MIL2k9iY2Nj3R5jecXFhwXy8/OxePFi1K9fH59//nmR+po1a/C3v/0Na9euRd++fd1+/IKTt3DuXCmFefPmmR5zt27dMGPGDBw4cMDpTy9fffWV0/3at28PX19ffPnll/jTn/7kuP3cuXPYuHGj45Mdq3l4eMDb29tpYktJSSnyDf7evXsjPj4ec+fOxTvvvGP4WD179kT16tXx+++/45lnninVcRO565FHHkGtWrVw5MgRjBo1yu3ti3uuAHD8cvzkk0/wyy+/IDExEf/zP//jdJ++ffti+fLlyM/PN30ZgeLOae3atYPNZsOKFSuc0jY7d+7EmTNnnBYfffv2xfvvv486deqU2C9xo3ECKLJPGjZsiPr162PBggUYN26c4TVC/Pz80K1bN+zbtw/NmzeHt7d3iYyxvOLiwwJr167FhQsX8MEHHxSJfgFAXFwcZs+ejfnz55tafDz++OPw9vbGc889h7/+9a+4efMm5s6di9TUVNNjHjNmDBYsWIA+ffpg2rRpCA0NxdKlS3H06FGn+9WqVQv/8R//gUmTJuHFF1/Ec889h6tXr2Lq1Knw8fHB22+/bXoM96Igvvf6669j4MCBSEpKwrvvvovw8HCn79d06tQJL7zwAqZNm4aLFy+ib9++sNls2LdvH/z8/DB69GjUq1cP77zzDiZPnow//vgDvXr1Qu3atXHx4kX885//hL+/vyMeSGS1GjVqYNasWRg6dCiuXbuGgQMHIiQkBJcvX8aBAwdw+fJlzJ07V7t9cc+VAsOGDcMHH3yAIUOGwNfXF88++6xTffDgwVi6dCmeeOIJvPHGG2jbti28vLxw7tw5bNq0CU899RQGDBggvqbizmmBgYEYN24c4uPjUbt2bQwYMADnzp3D1KlTER4ejmrV/u+bBWPGjME333yDzp07Y+zYsWjevDlu376Ns2fP4qeffsL48ePdXix16NABtWvXxogRI/D222/Dy8sLS5cuxYEDB4rcd86cOejXrx8efvhhjB07FnXr1sXZs2exfv16LF26FAAwc+ZMdOzYEZ06dcJrr72GevXqITMzEydPnsT3339f5GJ0FVpZftu1qujfv7/y9vZWly5d0t5n8ODBqnr16iolJUV7HynV8f3336sWLVooHx8fFRkZqf7t3/5NrV27tkgstkuXLqpp06ZFth86dKiKiYlxuu3IkSPq8ccfVz4+PiowMFANHz5crV69ushjKnUnNty8eXPl7e2t7Ha7euqpp4pEUHXf7NeNqTgpFul+06dPV/Xq1VM2m001btxYzZs3z/GN9cLy8/PVxx9/rOLi4hzjb9++vfr++++d7vftt9+qbt26qZo1ayqbzaZiYmLUwIEDXUYZiXSKk3b5xz/+4XT7qVOnFAC1cOFCp9u3bNmi+vTpowIDA5WXl5eKjIxUffr0KbK9keKeKwU6dOigAKjnn3/esJ6bm6v+67/+yzEn1ahRQzVq1Ei9+uqr6sSJE477lcScdvv2bTVt2jQVFRWlvL29VfPmzdWaNWtUixYt1IABA5weMysrS7311lvqgQcecJzrzZo1U2PHjhXnXqX0aZft27er9u3bKz8/PxUcHKz+/Oc/q7179xq+Rzt27FC9e/dWdrtd2Ww2Vb9+fTV27Fin+5w6dUoNGzZMRUZGKi8vLxUcHKw6dOigpk2bJo6vovFQSqmyWvgQERGVtFOnTqFRo0Z4++23MWnSpLIeDhng4oOIiCqsAwcOYNmyZejQoQNq1qyJY8eOYcaMGcjIyMChQ4fEL9tS2eF3PoiIqMLy9/dHYmIi5s+fj7S0NNjtdnTt2hXvvfceFx7lGD/5ICIiIkvxImNERERkKS4+iIiIyFJcfBAREZGlSu0Lp59++ik+/PBDJCcno2nTpvjkk0/QqVMnl9vdvn0bFy5cQEBAgOFld4mo9CmlkJmZiYiICKcLNZU2s/MGwLmDqKy5NW+UxsVDli9frry8vNS8efPUkSNH1BtvvKH8/f3VmTNnXG6blJSkAPCHP/wpBz9JSUmlMUUYupd5QynOHfzhT3n5Kc68USppl3bt2uGhhx5yuqRv48aN0b9/f8THx4vbpqeno1atWli5ciX8/f2L1I2uiV9AasqWmZmprRW0gjaSl5enrWVkZGhr6enp2lqtWrW0NQDw9fXV1qS3q0aNGtpabm6utnZ3E6TCpP+DlDrOSttJ7b3tdru2lpSUpK1J7+/dXX0LK+g2aeSPP/7Q1qR9BqBIN9/CpP0mvU9STTpOpeO7oDX63ZRSyMvLc0QXrXAv8wbwf3PHJ598YngOSa9D1+YckPef1H9Dmo+k55POf2n+c/W40ljNvg7p/26lucrsdrrj1ZX8/HxTjyltZ/ZcleY/ALh165a2lpWVZWo76TVK86O7Y7lx4wbeeOONYs0bJf5nl1u3bmHPnj148803nW7v0aMHtm/fXuT+OTk5Tt0TC36J+Pv7Gy4+pJNLOkmknW928SHVpDfNz89PW3NVN7v4kMYj1aRFhDSBSNtJ74X02qX3XjrZpeeTJl3puHDVfVjaN2b3m9W14tRLirvzBqCfO3x9fQ1/gUvHlvQL3+ziQ9pOOpalcd7L4kPalosP97eT5hxpTnX15wijbsIFpH3j5eWlrUmvURqPmf9xAYo3b5T4H3OvXLmC/Pz8Ihd3CQ0NNWzRHh8fD7vd7viJjo4u6SERUTnn7rwBcO4gqshK7Ztkd698lFKGq6GJEyciPT3d8SN9tE5ElVtx5w2AcwdRRVbif3YJCgqCp6dnkf9buXTpkuGlbm02m8uPEomocnN33gA4dxBVZCW++PD29karVq2QkJCAAQMGOG5PSEjAU089VezH8fT0NPxbo/RlHemLM9KX/BYtWqStSX9DlfoGJCcna2uuvvMh/Q0uLi5OW5P+bh0SEqKtSa9D2qfXrl3T1u6//35tTfoypvTFUenviBcuXDC1nfR/y9IXTlNTU7U1QP5bsfR3W+nvyGaZ+XtvKXwPXVRS84ZEmjukL/JJ56P0PkvHnbRoMvt9EED+Um1gYKCp7Vw9p47Z7wuYJR2z0vfzpPdemv+kuUo61qSgAiB/X0Q63qSxSsei2e88Gj2f9Hu2yPMW+55uGDduHF544QW0bt0a7du3x9///necPXsWI0aMKI2nI6JKgPMGUdVRKouPZ599FlevXsU777yD5ORkxMXF4ccff0RMTExpPB0RVQKcN4iqjlK7wunrr7+O119/vbQenogqIc4bRFUDe7sQERGRpbj4ICIiIktx8UFERESWKrXvfNyr9PR0w5jPjh07tNskJiZqa1LUtHPnztqaFKfSXfYZkPu3nDt3TlsD5Muk//zzz9qaFLWqV6+ettagQQNt7dFHH9XWpGivFN2T+t5IsTcpSidFbffs2aOtmb0suxRBA+RomxSXky6RLG0njaciRG1L0rVr1wzjqmajmFJ8UIpFSqSxSMeOUcuJ4taDgoK0tbp162pr0twhnedSnNjs+SEdy9K5LMX7L126ZKomXU5BmuOkeRowH0OWLoMvxWmlmhSzNnoP3Tkf+MkHERERWYqLDyIiIrIUFx9ERERkKS4+iIiIyFJcfBAREZGluPggIiIiS5XbqO2iRYsMY4d9+vTRbvPcc89pa2FhYdqaFE+TYkhnz57V1qSI1po1a7Q1ADh69Ki25u3tra1JES4pviVFf6XolLTfpNjb5cuXTdWOHz+ure3du1dby8nJ0dak1ydF16TYIgBERkZqa9KxWKdOHW3t7nbzhUkdeKVjURclV0qJkcfyLCMjw/B9vX79unYbKVIvHT9SNFo658x2PXYVgZbi71KEX4qqS1FjKaIrnSPSvCrFcKV5RXp/r1y5oq1J860Up5XOR2ksrmL60nssdViWos1m97d0fBvN/9K5cjd+8kFERESW4uKDiIiILMXFBxEREVmKiw8iIiKyFBcfREREZCkuPoiIiMhS5TZqO2LECMMoT8OGDbXbSNFPKb6UmpqqrZ0/f15bO3nypLYmRbukSBQANGvWTFurWbOmtibFomrXrq2tSTHD0NBQbU2K9V27dk1by8jI0NZOnz6trUlx2uzsbG1N2mdNmjTR1qTOxOHh4doaIEfipJr0OiTSazTTYVMpJcYay7PMzEzDqK0UcZQi7FJEVYowSvtPqknRVlfHh3QuS8edNAeeOnVKW5P2m9luqdJ2UjReeu1SZFaqSdF/6byS3kNpngZg2JG5gPQ7TtpOek7p95F0zrh7jt2Nn3wQERGRpbj4ICIiIktx8UFERESW4uKDiIiILMXFBxEREVmKiw8iIiKyVLmN2gYHBxtG3KTOrWY7M0rdaY8dO6atSd1gpdiX1EUQkKN9UtdTqVtq/fr1tbU2bdpoawEBAdqaFIk7ePCgtiZF2/bv329qO+m4kF5DVFSUtiaR3l9AjnZLMWzpcaVOqPdyvFU2fn5+hvFBs92rJVIUU4qvSnFzqSOqq5i+9F5LXZqlx5XGevXqVW0tODhYW5Ni+lIkWDoHpDlAei+k+HJWVpa2JsV+pd83UpdpQN430lwWERFh6jGl2LO0nVHXZuk9uFuJf/IxZcoUeHh4OP1IvxSJiDhvEFUtpfLJR9OmTbFhwwbHv6UVNxERwHmDqCoplcVH9erV+X8tROQWzhtEVUepfOH0xIkTiIiIQGxsLAYPHow//vhDe9+cnBxkZGQ4/RBR1ePOvAFw7iCqyEp88dGuXTssWbIE69evx7x585CSkoIOHTpov5gUHx8Pu93u+ImOji7pIRFROefuvAFw7iCqyEp88dG7d28888wzaNasGR577DH88MMPAIDFixcb3n/ixIlIT093/CQlJZX0kIionHN33gA4dxBVZKUetfX390ezZs1w4sQJw7rNZjOMeWVmZhpGeaR4a2JiorYmReIOHz6srV28eFFbkyJxUsdDqVuqq7rUubBly5baWlBQkLYmdWiVYr9SfE3aN9J7aDZOazZKvWbNGm3tlVde0dakYw2QjykpoifVjM6HAmY7tuoeU4oKW8HVvAHo546QkBDDiKB0/EjHpDR32O12bU3qei1F0aW4bKNGjbQ14M6XdnXS0tK0NSnCabaTbE5OjrYmHZNS7FfqBiw9n/TeS9tJ55UUUZXO1U2bNmlrrsYjRV/btWunrUljlY7hmJgYbc3o94Y7kfVSv8hYTk4OfvvtN5ctyImICnDeIKrcSnzxMWHCBGzZsgWnTp3Crl27MHDgQGRkZGDo0KEl/VREVElw3iCqWkr8zy7nzp3Dc889hytXriA4OBgPP/wwdu7cKX58Q0RVG+cNoqqlxBcfy5cvL+mHJKJKjvMGUdXCxnJERERkKS4+iIiIyFLltqvttWvXDOOTu3bt0m4jXRHx559/1tYGDBigrUmRWalTrBSnctURVYpMSVdxTElJ0dakaJs0HinaK8XXDh06pK1JEUQpoivFP6V4otQJU3q+efPmmXpMQO6GKXXtlPqZmO32aSa+W9ZR23tRq1Ytw2ihdC5L547ULVvqMi11Gd25c6e2Jh0DrvrdSF2apZiqFCmVjjsp4m42wiode1KE1dU5qSN1O5Y60EpzqhR5Xr9+vTgeaZ9KUdvGjRtra9IlI6SYuRR7NupaLL2vd+MnH0RERGQpLj6IiIjIUlx8EBERkaW4+CAiIiJLcfFBREREluLig4iIiCxVbqO2V65cMYwVnTlzRruNFLXt1auXtvbqq69qa1Kc9sMPP9TWpBiW1EESkON70uWma9eura0dP37c1HikaJcU0ZLiu1KXUHeiWoWZjeFK8UPpNbiKPErPKW0rRQmlyKP0fNKxKMUhK2rctlatWobxyfPnz2u3CQ0N1dak80NqfFezZk1trU2bNtqaFGF3dbl5qauol5eXtibFTaXjwGz8WzrvpJrZrs/SfpFeu9RhXHo+6TXExcVpa4D8/ktRaqkm/U6R4sTSMWzUDVc6xu7GTz6IiIjIUlx8EBERkaW4+CAiIiJLcfFBREREluLig4iIiCzFxQcRERFZiosPIiIislS5vc5HbGysYf66Y8eO2m1atmyprUltiufMmaOtSS3XW7Vqpa1J+XcpOw0A2dnZ2lpgYKC2ZpS7LiBdz0HKuUvX8vjtt9+0NekaGdLrk7Lz0rUszNbMtu92lWevV6+etia9fuk6BlL+32zbc52Keo0P4M45UqNGjSK3S/NDgwYNtDVpv2dlZWlr0jUiBg4cqK1Jx4B0vQpAfq+lY1KaA6RjS3qN0lilc0u6PojZ69kEBQWZekxpnNJxIZ0/ffv21dYAeX9L75M0/9tsNm2tefPm2prReVTAx8enyG3S3HY3fvJBREREluLig4iIiCzFxQcRERFZiosPIiIishQXH0RERGQpLj6IiIjIUm5Hbbdu3YoPP/wQe/bsQXJyMlatWoX+/fs76kopTJ06FX//+9+RmpqKdu3aYc6cOWjatKlbz5OTk2MY1YqMjNRuYxT9KY5GjRppa1LsTYqvnjx5UlvLyMgQxyO9Duk5pVqLFi20NSkWJr1+qTV8amqqtmY2+iqRxiJF9yTSOKX9Asjvcf369bU1KaqWkpJiajupvbcu1qeUErdzl1XzBnAncmgUEa0ZbFgAACAASURBVAwICNBuYzb6KZ1z0nZmI5yu5g7pPZOOZ2msUkxTine60169uGORXoM0TukSBb6+vtpaeHi4tiZdTkB6f6V4MiDvUylOLEWbpZr0PkmxX6M5V9onRR672Pf8/7Kzs9GiRQvMnj3bsD5jxgx89NFHmD17Nnbv3o2wsDA8/vjjbg2KiCoXzhtEVJjbn3z07t0bvXv3NqwppfDJJ59g8uTJePrppwEAixcvRmhoKL766iu8+uqr9zZaIqqQOG8QUWEl+p2PU6dOISUlBT169HDcZrPZ0KVLF2zfvt1wm5ycHGRkZDj9EFHVYWbeADh3EFVkJbr4KPibdGhoqNPtoaGh2r9Xx8fHw263O36io6NLckhEVM6ZmTcAzh1EFVmppF3u/lKQUkr7RaGJEyciPT3d8ZOUlFQaQyKics6deQPg3EFUkZVoY7mwsDAAd/5PpvC3hC9dulTk/2oK2Gw28VvKRFS5mZk3AM4dRBVZiS4+YmNjERYWhoSEBEcHyVu3bmHLli344IMP3Hqs3Nxcl3HGu0kxJClmJ3UKzMvL09akb+JLcbGLFy9qa64e96GHHtLW6tatq61J/wcpRfukyJgU67t69aq2JnUKNhu1lSJhErNxSGl/AkBaWpq2tmfPHm1N6r5pNoInnUe641sphfT0dO12Jakk5w3gTlTdKD4pxbElZjsfS8eWVJPeLyky6upxJdI8J3XKlc47s11tpZr0XkiXKJAeU5r/pcWt9Nql/elq7pCeU4rFSjVX3ZDNMHovpNd9N7dHlJWV5XQNi1OnTmH//v0IDAxE3bp1MWbMGLz//vto0KABGjRogPfffx9+fn4YMmSIu09FRJUE5w0iKsztxUdiYiK6devm+Pe4ceMAAEOHDsWiRYvw17/+FTdu3MDrr7/uuFjQTz/9JH7yQESVG+cNIirM7cVH165dxY++PDw8MGXKFEyZMuVexkVElQjnDSIqjL1diIiIyFJcfBAREZGluPggIiIiS5V8/qaE1KxZ0zACJUXNpEiYFG+qU6eOtnbp0iVtTYovSXE5qVMicOfv42aeU4p+SXHaM2fOmHo+KVYldft09bf/kt7ObMTyXjrlStE+aZ9KEV3p9Rt1cS1gptvz7du3LYvaljTd9T+kKKLZ485sFN3s8eHn56etuSKdr9Kl6aVz2WznWmmfSueW2fdCIkVbpY63Urxd2i+uxinN42ajzdJ2rqK/Okavw53LY/CTDyIiIrIUFx9ERERkKS4+iIiIyFJcfBAREZGluPggIiIiS3HxQURERJYqt1Hb++67z7CvgxR9lSJT2dnZ2tqpU6e0NSlmJsXeIiMjtbXg4GBtDZAjU1LUTIpTSZ1ypciYtN+ysrK0taioKG2tcIOxu1kd0S2tzpRmo4RSDE+KvkqdgqXjTRfDNdsdtTyoVq2a4ftqNsIokfaT2ai22fiqq8eVYpBSLFh6HWY7qUqPKb0XZiPKZiO60mNKr116fdJ862pbad4pjePb3c7e7lzWgJ98EBERkaW4+CAiIiJLcfFBREREluLig4iIiCzFxQcRERFZiosPIiIislS5jdr6+voaRlkjIiJMPZ4UC5UiU1KnXClWJHUSleJbAHDt2jVtLSwsTFuTYlFSFDM6Olpbu3DhgrYmdeeVYn3Sa5C6a5Zk98XikN6ne+lMKW1rNhInRSVTUlK0NV3XTrP7rDzIy8szjElLx4/ZjrfSMWK20/K9dG6VzjupJkV0zXZSlTrCSvOjFEU1+xqkyyJIc470Goy6rhcwG8MFzMdpzV5uwOzxZlRzZ97gJx9ERERkKS4+iIiIyFJcfBAREZGluPggIiIiS3HxQURERJbi4oOIiIgs5XbUduvWrfjwww+xZ88eJCcnY9WqVejfv7+j/tJLL2Hx4sVO27Rr1w47d+50b2DVq7uMpBptoyPFkMxEEV09plSTuqUCQO3atbU1Kaa1f/9+ba1hw4baWv369bU1KRJ37tw5be38+fPaWmBgoLZmNqIqRU3NRmZdvU9mSa/RbOTR3e6TBXRdi6X4nRlWzRvAnc7IrjqH3s1sl1WzHUjNdrx11W1YOg6kDtVmOzhLlyIw27lXOl+l81x6PrNzdVpamrZmNi7s6pIBZuO0pXEpAmkeMHov3Jkz3f7kIzs7Gy1atMDs2bO19+nVqxeSk5MdPz/++KO7T0NElQjnDSIqzO1PPnr37o3evXuL97HZbOKFpIioauG8QUSFlcp3PjZv3oyQkBA0bNgQr7zyCi5duqS9b05ODjIyMpx+iKjqcWfeADh3EFVkJb746N27N5YuXYqNGzfib3/7G3bv3o3u3bsjJyfH8P7x8fGw2+2OH+lS30RUObk7bwCcO4gqshLv7fLss886/jsuLg6tW7dGTEwMfvjhBzz99NNF7j9x4kSMGzfO8e+MjAxOIkRVjLvzBsC5g6giK/XGcuHh4YiJicGJEycM6zabDTabrbSHQUQViKt5A+DcQVSRlfri4+rVq0hKShI7oBq5deuWYQdDKWol/Y1Yim8FBQVpa1J8SxqLFKXLzMzU1gAgNTVVW7ty5Yq21qJFC21N6pQrxcliYmK0tbNnz2prUvfJY8eOaWuuooQ60vsrxb+kWJsUM7uXuJzEVcdLM0o6NmsFs/MGcOe8NDo3zUauSyP6aPY9cbWdFH2VatJcJs050j6V9ptUkxaU0twoxfulmLHUKV36s5/UKVzqouvqEhJmu9pKzB5v0nxsNMdL8/7d3F58ZGVl4eTJk45/nzp1Cvv370dgYCACAwMxZcoUPPPMMwgPD8fp06cxadIkBAUFYcCAAe4+FRFVEpw3iKgwtxcfiYmJ6Natm+PfBX9zHTp0KObOnYtff/0VS5YsQVpaGsLDw9GtWzesWLECAQEBJTdqIqpQOG8QUWFuLz66du0qfoSzfv36exoQEVU+nDeIqDD2diEiIiJLcfFBREREluLig4iIiCxV6lFbs65fv24YO5Q6vprtQOsq+qojxbekaKsUXQPuJAN0pO60Fy9e1NbOnDmjrUnxrcIJhbsdPnxYW5NifVJ8rTQ6t0rfNTD7mK7icmajm2aZjZFWRvn5+YbxaikGKHUhleLPpRGZNVu7l/FIc5I0B0jHlhR9ve+++7S1mjVramtS1PbAgQPamhS3l2LGwcHB2pq7HV8LuLo2jdn4v9maFKeV9ptRrVS72hIRERHdCy4+iIiIyFJcfBAREZGluPggIiIiS3HxQURERJbi4oOIiIgsVW6jtjVr1jSMXEmRISluKG0nRZ+kGJIUiZLieVKUFpC780rxzyNHjmhrISEh2poUCz516pS2dvz4cW1N6vgoddGVeHl5aWvSey9FLM123nTV1dZVFNcMKQ5qtgOv7vhWSomR6PJM1xFbOn7Mdqg2e4yYjUW6ijKajZVL77W0byTSeSdFe2/cuKGtSRF+f39/bU26LIJUCw0N1dak/SnVXHXuvpd5x8xzSseU9DqMjgt3jhV+8kFERESW4uKDiIiILMXFBxEREVmKiw8iIiKyFBcfREREZCkuPoiIiMhS5TZqm5eX51aHPECO+ZiNL0kxOymiKsVppRgqAAQGBmprAQEB2toDDzygrR07dkxbkzpMdu/eXVtLT0/X1ubNm6etSfvUbHxNIu0zqda2bVttTequDMjxNakzZ0pKirYmxROl402KJ0pR24oqJyfHMFYrdb02G42UjmVJaURiXW0rRY39/Py0NSluKkXKa9Wqpa1JsVhpHvfx8TFVk+Zc6Vw22+3YbLQVkH8fmf09Jo1HOmaksRo9pjtzND/5ICIiIktx8UFERESW4uKDiIiILMXFBxEREVmKiw8iIiKyFBcfREREZCm3orbx8fFYuXIljh49Cl9fX3To0AEffPCBU8QzJycHEyZMwLJly3Djxg08+uij+PTTTxEVFeXWwK5fv24YZZLiRFIsSoqZSXEiqcOiFO2SIlFSh1lAjn9KUTspMivFwqRInBQJlJ7vX/7lX7S1bdu2aWu///67tiZ1HzbqgFxAeu+l+KUUUQ0LC9PWADlOayb6CgCZmZnamhRPNNO1sqSjtlbOHfn5+YavWTp3pLmjNLr7SpFEqRusq7FI20pzgDTnSHFaaT6W4rTSWMzO49J20uuT5nFpnNK+lsYpnf+AfGxIr9HsY0rnujQfWdrVdsuWLRg5ciR27tyJhIQE5OXloUePHk4ticeMGYNVq1Zh+fLl2LZtG7KystC3b1+XbYSJqPLi3EFEhbn1yce6deuc/r1w4UKEhIRgz5496Ny5M9LT0zF//nx88cUXeOyxxwAAX375JaKjo7Fhwwb07Nmz5EZORBUG5w4iKuyevvNRcIXLgity7tmzB7m5uejRo4fjPhEREYiLi8P27dsNHyMnJwcZGRlOP0RUuXHuIKraTC8+lFIYN24cOnbsiLi4OAB3Lg3t7e1d5HK1oaGh2stGx8fHw263O36io6PNDomIKgDOHURkevExatQoHDx4EMuWLXN5X6WU9otJEydORHp6uuMnKSnJ7JCIqALg3EFEphYfo0ePxnfffYdNmzY5fRM9LCwMt27dKvJN/0uXLmmbE9lsNtSsWdPph4gqJ84dRAS4+YVTpRRGjx6NVatWYfPmzYiNjXWqt2rVCl5eXkhISMCgQYMAAMnJyTh06BBmzJjh1sCUUoYRICn2JUURpaiVFF+Svmkv1aRxSl0iXZHia1KEKzIyUls7ceKEtma2i6L0EXijRo20tcuXL2trUsxQeu+l91eKhl29elVbk147IMeQzR43EjNxWitZOXfk5uYavq/SPpJik1IUU9pO2u/ScSfVpOg/IJ8j0lil81yKuLs6D3Sk89VsV1epZjZqKh0zZiPYrtJb7naSLWC2W7j0fO5G+N3pRO/WTDdy5Eh89dVXWL16NQICAhx/i7Xb7fD19YXdbsfw4cMxfvx41KlTB4GBgZgwYQKaNWvm+AY7EVU9nDuIqDC3Fh9z584FAHTt2tXp9oULF+Kll14CAHz88ceoXr06Bg0a5LhQ0KJFi8RVGRFVbpw7iKgwt//s4oqPjw9mzZqFWbNmmR4UEVUunDuIqDD2diEiIiJLcfFBREREluLig4iIiCxlLtdnATNxOSkS504EqDApnnbx4kVtTepcK0XCALkjohRRu379urZWcDlrIwWXuDZy5swZbe3kyZOmHlOKPUuxPul7A1lZWdqa9N5LNanjbWJiorYGAFeuXNHWpBiy9N5LEUTpdUj7zWxUsjzTzR1SNNLsfpe2MxuZleYcV1Fb6Tmlx72XTro60n6Tjknp+aTIqKt5VcedTqyFSfvM7JwDyPOO2TixFNE1G8M1ev3SPrkbP/kgIiIiS3HxQURERJbi4oOIiIgsxcUHERERWYqLDyIiIrIUFx9ERERkqXIbtc3LyzOM+UhR06CgIG3NbNxQiq9KMTupvberqJUUG5VqtWrV0tYuXLigrUn7RuqyKkXCfv/9d23t3Llz2prdbtfWpP0t1aROuVJ0UYqnSa8dAH777TdtraCpmpGwsDBtTTpOzUZmde99cS6HXl7dvHnTMJIoHSPScS71ljHbEVXqQC2d49K5AwAZGRmmnlMaq3SOSNFPad9Inaal98ksaZzSey+NRbosgPR80jkOAMHBwdqadCxKNWl+kKLG7h4zrjr2FsZPPoiIiMhSXHwQERGRpbj4ICIiIktx8UFERESW4uKDiIiILMXFBxEREVmq3EZtlVKGUS0pTiTFfKTIlBRPO336tLYmRaakeJoUMwPkeJM01q1bt2prUuxN2m9paWnamhTRunbtmraWmpqqrUmdiSVSF11pf0rdfqXtwsPDxfFI778Ue5b2jat4rxm6SGBFjtpmZmYaHtPSsSXFJqXjwGxXW6n7p3SOS520ASA5OVlbk45J6f2W5lwppioxO1dLHVjNHrNSDFWaq6VjxuwlAwD5OJXeC7Nzp3R8Sx2GjS5DIV2a4m785IOIiIgsxcUHERERWYqLDyIiIrIUFx9ERERkKS4+iIiIyFJcfBAREZGl3Fp8xMfHo02bNggICEBISAj69++PY8eOOd2na9eu8PDwcPoZPHhwiQ6aiCoWzh1EVJhbIe0tW7Zg5MiRaNOmDfLy8jB58mT06NEDR44cgb+/v+N+r7zyCt555x3Hv81cn6BatWqG1yAwmyuX8sdSC2uzLc7/+OMPbe3SpUvaGoAik3Jh0jUipJx3ZGSktibt01q1amlr0jVAgoKCtDXpvZCucyFdk+P8+fPamnRtAOm1S9vFxMRoawAQEBCgre3Zs0dbk67FIF03wuw1DnSvsaSv82Hl3JGammp47QI/Pz/tNq6uvaMj7SfpWh7utB4vzNW1HEJCQkp8PNJzSq3jpflR2m/SeSc9X2lcc8SobXwBaZzS8eRqnGaPDen9lcYqXcvD3evRuHOdD7ferXXr1jn9e+HChQgJCcGePXvQuXNnx+1+fn7iL20iqlo4dxBRYff0nY+C/xO9++qSS5cuRVBQEJo2bYoJEyYgMzNT+xg5OTnIyMhw+iGiyo1zB1HVZvry6kopjBs3Dh07dkRcXJzj9ueffx6xsbEICwvDoUOHMHHiRBw4cAAJCQmGjxMfH4+pU6eaHQYRVTCcO4jI9OJj1KhROHjwILZt2+Z0+yuvvOL477i4ODRo0ACtW7fG3r178dBDDxV5nIkTJ2LcuHGOf2dkZCA6OtrssIionOPcQUSmFh+jR4/Gd999h61btyIqKkq870MPPQQvLy+cOHHCcAKx2Wxigx4iqjw4dxAR4ObiQymF0aNHY9WqVdi8eTNiY2NdbnP48GHk5ua67AJKRJUX5w4iKsytxcfIkSPx1VdfYfXq1QgICEBKSgqAO+2DfX198fvvv2Pp0qV44oknEBQUhCNHjmD8+PFo2bIlHnnkEbcGVpDzLzJgIaYkRbtq1qyprUkxJOlLbGfPntXWNm7cqK3Vrl1bWwPkeGvdunW1NSmKWaNGDW2tcNTxbteuXdPWJNJ+k94LqWW0FAnLzs42tZ3ZWN+GDRu0NQBo3769thYaGqqtFZxTRqSxmiWdMyXJyrkjOzvbMK4oRcPNRm2lKKZ0/JiNabr6tEiKOkqxSWms0nik40d6jdJYpH0qzf/SOKU29tJrNxsJluLJrs45aTzSXCa9Rul3gxTvly5DYXQ+SY91N7fSLnPnzkV6ejq6du2K8PBwx8+KFSsA3Hnzf/75Z/Ts2RMPPPAA/vKXv6BHjx7YsGGD+EuFiCo3zh1EVJjbf3aRREdHY8uWLfc0ICKqfDh3EFFh7O1CREREluLig4iIiCzFxQcRERFZiosPIiIispTpK5yWtuvXrxtGjqRuodKX2qT41tWrV7U1KU4lRWIHDhyordntdm0NkONUUldbKd4qxbuOHj2qrSUlJWlrkt9++01ba9iwobYmdR6VopJSlMwsqbukq+eTOtdKcWkpqia9fol0DFfGJElWVpbh+yO9Vun8kN5r6XiVnk+aq6QLp7nqaivF2KXjWTpGpJq0b6SYpjQfS6T3Sdpv0vskdU6WYq9mO/qa7VoLyPtb2qdGXZ4LSJcpkDqJJycnu/U8d+MnH0RERGQpLj6IiIjIUlx8EBERkaW4+CAiIiJLcfFBREREluLig4iIiCxVbqO2GRkZhpEkKWorRcKkaJsUT8vMzNTWpE6xUuRIitICckRNitpJkWEp9idF1KT9Jo1F6txrtvtiZGSktiZ1u3QVbda5ePGituaqe6MUe0tNTdXWpEig9JxSrTLGaSXXr193O84o7SPpmJSOc6nLqnTOmY13uqpL54h0nkv70mxkWDo/pJiqNOdIkVlpO7Mdb81Gl129h2Y7/krPKf1OkTqXSzWj53On8zY/+SAiIiJLcfFBREREluLig4iIiCzFxQcRERFZiosPIiIishQXH0RERGSpchu1DQkJMYzVSpFCKdomde6T4mJmY6FSdO/y5cvaGiBHca9fv66t3XfffdqaFIEKDAzU1qTXeOXKFW0tODhYW5MieNJrlyLR/v7+2prUmfH8+fPamtnYKyC//1JcTjoW3Ymx3etjStuUdx4eHoZxTSn+LkXqJdJ7IkXxpZim2YguIB93Uk2K4UpzgLSd9Dqk81U6P6TIrLRvpJq0X8xevkHazlVHX+k4lY4b6RiWOmJLv1Ok+djouHCnuzg/+SAiIiJLcfFBREREluLig4iIiCzFxQcRERFZiosPIiIishQXH0RERGQpt6K2c+fOxdy5c3H69GkAQNOmTfGf//mf6N27N4A7sZwJEyZg2bJluHHjBh599FF8+umniIqKcntgt2/fNowVSdFXKf4oRYCkmhQXk+JUUnxLiosBcsQzJCREW5O6OkodWqXIrESKdkkxNKlT4oULF7Q1KS4nHRfS65OiZFI8W4oYAvKxIXWtLK1OmWXN6rnDaF9J57L0nkhRbSk2Kc1HUgxX6sLs6n2WOsJKx4/ZTrIS6TVKY5HOAen1S3OAWdI8ZrYTsBRtBe50dNcxO19JjymRXofRWFx1bC/MrRkrKioK06dPR2JiIhITE9G9e3c89dRTOHz4MABgzJgxWLVqFZYvX45t27YhKysLffv2dbu9NRFVLpw7iKgwtz756Nevn9O/33vvPcydOxc7d+5EVFQU5s+fjy+++AKPPfYYAODLL79EdHQ0NmzYgJ49e5bcqImoQuHcQUSFmf6sNj8/H8uXL0d2djbat2+PPXv2IDc3Fz169HDcJyIiAnFxcdi+fbv2cXJycpCRkeH0Q0SVF+cOInJ78fHrr7+iRo0asNlsGDFiBFatWoUmTZogJSUF3t7eqF27ttP9Q0NDkZKSon28+Ph42O12x090dLT7r4KIyj3OHURUwO3FxwMPPID9+/dj586deO211zB06FAcOXJEe3+llPhlpokTJyI9Pd3xk5SU5O6QiKgC4NxBRAXcbizn7e2N+++/HwDQunVr7N69GzNnzsSzzz6LW7duITU11en/YC5duoQOHTpoH89ms7lslkREFR/nDiIqcM9dbZVSyMnJQatWreDl5YWEhAQMGjQIAJCcnIxDhw5hxowZbj/u9evXDWNV0v8JSZFZKS4mdViU4ktGXXeLw1VXQ6kjrBQJlGKqWVlZ2pq036RYmBTtleLCycnJ2prUYffs2bPamhRrlF6D9F6429GxMCkSKL2HZqN9EjMRXSu62pbW3FG9enXDKLTZbspS4sZsB1aznY1dHXdS52fpfHUVHdcx2y3W7PElPab0u0F6D81G383GrKXfKYB8CQPpcaVuuNJxI0WU3Y0Tu5NOc+uImzRpEnr37o3o6GhkZmZi+fLl2Lx5M9atWwe73Y7hw4dj/PjxqFOnDgIDAzFhwgQ0a9bM8Q12IqqaOHcQUWFuLT4uXryIF154AcnJybDb7WjevDnWrVuHxx9/HADw8ccfo3r16hg0aJDjQkGLFi0SV6tEVPlx7iCiwtxafMyfP1+s+/j4YNasWZg1a9Y9DYqIKhfOHURUWPm+JjMRERFVOlx8EBERkaXuOe1S0gq+Xav7RrDZtIv0zWRpO+mbya5SK2YeE5C/JW92rFLiQ/oGtZQCkL4lbvZ9kr4tbbYmpT2k2r0kT8xuW15qBbdbkXopKQVj1R2XUgJJOs6l7cx+J0V6TKkBnnReAXJqpTSauUnncmVIu5hNw5lN37na1uxcbTbN5e5jFuzL4ry/HqqczS7nzp3jlQqJyomkpCRTnWXLAucOovKhOPNGuVt83L59GxcuXEBAQAA8PDyQkZGB6OhoJCUliRn2qob7RY/7xpg7+0UphczMTERERLhs415ecO4oHu4XY9wvesXdN+7MG+Xuzy7VqlUzXDHVrFmTB4QB7hc97htjxd0vdrvdgtGUHM4d7uF+Mcb9olecfVPceaNi/C8NERERVRpcfBAREZGlPKdMmTKlrAfhiqenJ7p27Wq6/0Blxf2ix31jrKrtl6r2eouL+8UY94teSe+bcveFUyIiIqrc+GcXIiIishQXH0RERGQpLj6IiIjIUlx8EBERkaW4+CAiIiJLlevFx6efforY2Fj4+PigVatW+OWXX8p6SJbbunUr+vXrh4iICHh4eODbb791qiulMGXKFERERMDX1xddu3bF4cOHy2i01omPj0ebNm0QEBCAkJAQ9O/fH8eOHXO6T05ODkaPHo2goCD4+/vjySefxLlz58poxNaZO3cumjdv7rgaYfv27bF27VpHvSrsF84dnDt0OHcYs3reKLeLjxUrVmDMmDGYPHky9u3bh06dOqF37944e/ZsWQ/NUtnZ2WjRogVmz55tWJ8xYwY++ugjzJ49G7t370ZYWBgef/xxZGZmWjxSa23ZsgUjR47Ezp07kZCQgLy8PPTo0cOps++YMWOwatUqLF++HNu2bUNWVhb69u0rdnGsDKKiojB9+nQkJiYiMTER3bt3x1NPPeX4xVLZ9wvnjjs4dxjj3GHM8nlDlVNt27ZVI0aMcLqtUaNG6s033yyjEZU9AGrVqlWOf9++fVuFhYWp6dOnO267efOmstvt6rPPPiuLIZaZS5cuKQBqy5YtSiml0tLSlJeXl1q+fLnjPufPn1fVqlVT69atK6thlpnatWurzz//vErsF84dRXHu0OPcoVea80a5/OTj1q1b2LNnD3r06OF0e48ePbB9+/YyGlX5c+rUKaSkpDjtJ5vNhi5dulS5/ZSeng4ACAwMBADs2bMHubm5TvsmIiICcXFxVWrf5OfnY/ny5cjOzkb79u0r/X7h3FE8nDv+D+eOoqyYN8rlNWSvXLmC/Px8hIaGOt0eGhqKlJSUMhpV+VOwL4z205kzZ8piSGVCKYVx48ahY8eOiIuLA3Bn33h7e6N27dpO960qx9Cvv/6K9u3b4+bNm6hRowZWrVqFJk2aYP/+/ZV6v3DuKB7OHXdw7nBm5bxRLhcfBTw8PJz+rZQqSTBt1gAAArdJREFUchtxP40aNQoHDx7Etm3bXN63quybBx54APv370daWhq++eYbDB06FFu2bNHev7Ltl6p+ThRXVd9PnDucWTlvlMs/uwQFBcHT07PIiurSpUtFVupVWVhYGABU6f00evRofPfdd9i0aROioqIct4eFheHWrVtITU11un9V2Tfe3t64//770bp1a8THx6NFixaYOXNmpd8vnDuKh3MH5w4jVs4b5XLx4e3tjVatWiEhIcHp9oSEBHTo0KGMRlX+xMbGIiwszGk/3bp1C1u2bKn0+0kphVGjRmHlypXYuHEjYmNjneqtWrWCl5eX075JTk7GoUOHKv2+MaKUQk5OTqXfL5w7iodzB+eO4ijVecP892BL1/Lly5WXl5eaP3++OnLkiBozZozy9/dXp0+fLuuhWSozM1Pt27dP7du3TwFQH330kdq3b586c+aMUkqp6dOnK7vdrlauXKl+/fVX9dxzz6nw8HCVkZFRxiMvXa+99pqy2+1q8+bNKjk52fFz/fp1x31GjBihoqKi1IYNG9TevXtV9+7dVYsWLVReXl4Zjrz0TZw4UW3dulWdOnVKHTx4UE2aNElVq1ZN/fTTT0qpyr9fOHfcwbnDGOcOY1bPG+V28aGUUnPmzFExMTHK29tbPfTQQ44oVFWyadMmBaDIz9ChQ5VSdyJzb7/9tgoLC1M2m0117txZ/frrr2U7aAsY7RMAauHChY773LhxQ40aNUoFBgYqX19f1bdvX3X27NmyG7RFhg0b5jhvgoOD1aOPPuqYQJSqGvuFcwfnDh3OHcasnjc8lFLK3GcmRERERO4rl9/5ICIiosqLiw8iIiKyFBcfREREZCkuPoiIiMhSXHwQERGRpbj4ICIiIktx8UFERESW4uKDiIiILMXFBxEREVmKiw8iIiKyFBcfREREZKn/B4mpNRTZ9z7+AAAAAElFTkSuQmCC",
      "text/plain": [
       "Figure(PyObject <Figure size 640x480 with 2 Axes>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "We have 32 faces, each 32x32 pixels."
     ]
    }
   ],
   "source": [
    "fo=open(\"img/faces.pkl\",\"r\")\n",
    "faces = pickle.loads(pybytes(read(fo)))\n",
    "\n",
    "# helper functions\n",
    "function display_face(f, cmap=\"Greys_r\", vmin=nothing, vmax=nothing)\n",
    "    plt.imshow(reshape(f, 32,32), cmap=cmap, vmin=vmin, vmax=vmax)\n",
    "    plt.tick_params(axis=\"both\", which=\"both\", bottom=\"off\", left=\"off\", labelbottom=\"off\", labelleft=\"off\")\n",
    "    plt.grid(false)\n",
    "end\n",
    "function reconstruct(model, X)\n",
    "    Z = model.transform(X)\n",
    "    try\n",
    "        return model.inverse_transform(Z)\n",
    "    catch(e) # no inverse_transform function\n",
    "        W = model.components_\n",
    "        return Z*W\n",
    "    end\n",
    "    # Note on the above: you can't just do Z@W for PCA because of the centring of the data, \n",
    "    # which needs to be undone for the reconstruction\n",
    "end\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.subplot(1,2,1)\n",
    "display_face(faces[rand(1:size(faces)[1]),:])\n",
    "plt.title(\"A random face\")\n",
    "plt.subplot(1,2,2)\n",
    "display_face([Statistics.mean(faces[:,i]) for i in 1:size(faces)[2]])\n",
    "plt.title(\"The average face\")\n",
    "\n",
    "N,D = size(faces)\n",
    "Printf.@printf(\"We have %d faces, each 32x32 pixels.\", M)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "## Testing it out with scikit-learn's PCA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "PyObject PCA(copy=True, iterated_power='auto', n_components=None, random_state=None,\n",
       "    svd_solver='auto', tol=0.0, whiten=False)"
      ]
     },
     "execution_count": 93,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pca = SkDecomposition.PCA()\n",
    "pca.fit(faces)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With  2 components, variance explained = 36%\n",
      "With  3 components, variance explained = 71%\n",
      "With  4 components, variance explained = 76%\n",
      "With  5 components, variance explained = 79%\n",
      "With  6 components, variance explained = 80%\n",
      "With  7 components, variance explained = 82%\n",
      "With  8 components, variance explained = 83%\n",
      "With  9 components, variance explained = 84%\n",
      "With 10 components, variance explained = 85%\n",
      "With 11 components, variance explained = 85%\n",
      "With 12 components, variance explained = 86%\n",
      "With 13 components, variance explained = 86%\n",
      "With 14 components, variance explained = 87%\n",
      "With 15 components, variance explained = 87%\n",
      "With 16 components, variance explained = 88%\n",
      "With 17 components, variance explained = 88%\n",
      "With 18 components, variance explained = 88%\n",
      "With 19 components, variance explained = 89%\n",
      "With 20 components, variance explained = 89%\n",
      "With 21 components, variance explained = 89%\n"
     ]
    }
   ],
   "source": [
    "cumulative_var_expl = cumsum(pca.explained_variance_ratio_)\n",
    "for i in 1:20\n",
    "    Printf.@printf(\"With %2d components, variance explained = %.0f%%\\n\", i+1, cumulative_var_expl[i]*100)\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Our own implementation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [],
   "source": [
    "pixel_means = [mean(faces[:,i]) for i in 1:size(faces)[2]]\n",
    "faces_centered = faces .- pixel_means'\n",
    "F = svd(faces_centered, full=false)\n",
    "Z = F.U*Matrix(Diagonal(F.S));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024,)"
      ]
     },
     "execution_count": 109,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(pixel_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2414, 1024)"
      ]
     },
     "execution_count": 110,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(Z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1024, 1024)"
      ]
     },
     "execution_count": 112,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "size(F.Vt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With  1 components, variance explained = 36%\n",
      "With  2 components, variance explained = 71%\n",
      "With  3 components, variance explained = 76%\n",
      "With  4 components, variance explained = 79%\n",
      "With  5 components, variance explained = 80%\n",
      "With  6 components, variance explained = 82%\n",
      "With  7 components, variance explained = 83%\n",
      "With  8 components, variance explained = 84%\n",
      "With  9 components, variance explained = 85%\n",
      "With 10 components, variance explained = 85%\n",
      "With 11 components, variance explained = 86%\n",
      "With 12 components, variance explained = 86%\n",
      "With 13 components, variance explained = 87%\n",
      "With 14 components, variance explained = 87%\n",
      "With 15 components, variance explained = 88%\n",
      "With 16 components, variance explained = 88%\n",
      "With 17 components, variance explained = 88%\n",
      "With 18 components, variance explained = 89%\n",
      "With 19 components, variance explained = 89%\n"
     ]
    }
   ],
   "source": [
    "for i in 1:19\n",
    "    reconstr = Z[:,1:i]*F.Vt[1:i,:]\n",
    "    var_expl = 1-sum((reconstr-faces_centered).^2)/sum(faces_centered.^2)\n",
    "    Printf.@printf(\"With %2d components, variance explained = %.0f%%\\n\", i, var_expl*100)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.7.0",
   "language": "julia",
   "name": "julia-0.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}