week1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Week 1: Mean/Covariance of a data set and effect of linear transformation\n",
    "\n",
    "In this week, we are going to investigate how the mean and (co)variance of a dataset changes\n",
    "when we apply affine transformation to the dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Learning objectives\n",
    "1. Get Farmiliar with basic programming using Python and Numpy/Scipy.\n",
    "2. Learn to appreciate implementing\n",
    "   functions to compute statistics of dataset in vectorized way.\n",
    "3. Understand the effects of affine transformations on a dataset.\n",
    "4. Understand the importance of testing in programming for machine learning."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are a few links for your reference. You may want to refer back to them throughout the whole course.\n",
    "\n",
    "- If you are less comfortable with programming in Python, have a look at this Coursera course https://www.coursera.org/learn/python.\n",
    "- To learn more about using Scipy/Numpy, have a look at the [Getting Started Guide](https://scipy.org/getting-started.html). You should also refer to the numpy [documentation](https://docs.scipy.org/doc/) for references of available functions.\n",
    "\n",
    "- If you want to learn more about creating plots in Python, checkout the tutorials found on matplotlib's website \n",
    "https://matplotlib.org/tutorials/index.html. Once you are more familiar with plotting, check out this excellent blog post http://pbpython.com/effective-matplotlib.html.\n",
    "\n",
    "- There are more advanced libraries for interactive data visualization. For example, [bqplot](https://github.com/bloomberg/bqplot) or [d3.js](https://d3js.org/). You may want to check out other libraries if you feel adventurous.\n",
    "\n",
    "- Although we use Jupyter notebook for these exercises, you may also want to check out [Jupyter Lab](https://github.com/jupyterlab/jupyterlab) when you want to work on your own projects."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's import the packages that we will use for the week. Run the cell below to import the packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "matplotlib.use('Agg')\n",
    "import matplotlib.pyplot as plt\n",
    "plt.style.use('fivethirtyeight')\n",
    "from sklearn.datasets import fetch_lfw_people, fetch_mldata, fetch_olivetti_faces\n",
    "import time\n",
    "import timeit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "from ipywidgets import interact"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we are going to retrieve Olivetti faces dataset.\n",
    "\n",
    "When working with some datasets, before digging into further analysis, it is almost always\n",
    "useful to do a few things to understand your dataset. First of all, answer the following\n",
    "set of questions:\n",
    "\n",
    "1. What is the size of your dataset?\n",
    "2. What is the dimensionality of your data?\n",
    "\n",
    "The dataset we have are usually stored as 2D matrices, then it would be really important\n",
    "to know which dimension represents the dimension of the dataset, and which represents\n",
    "the data points in the dataset. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "downloading Olivetti faces from http://cs.nyu.edu/~roweis/data/olivettifaces.mat to /home/jovyan/scikit_learn_data\n",
      "Shape of the faces dataset: (400, 4096)\n",
      "400 data points\n"
     ]
    }
   ],
   "source": [
    "image_shape = (64, 64)\n",
    "# Load faces data\n",
    "dataset = fetch_olivetti_faces()\n",
    "faces = dataset.data\n",
    "\n",
    "print('Shape of the faces dataset: {}'.format(faces.shape))\n",
    "print('{} data points'.format(faces.shape[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When your dataset are images, it's a really good idea to see what they look like.\n",
    "\n",
    "One very\n",
    "convenient tool in Jupyter is the `interact` widget, which we use to visualize the images (faces). For more information on how to use interact, have a look at the documentation [here](http://ipywidgets.readthedocs.io/en/stable/examples/Using%20Interact.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD3CAYAAAAwh5neAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnWuMXdWV5/+rXq7yowqXbbCNCwPBEMCiwUKJCUlEQ54k\nachklGdHnhHSfJlpJUpPBqdHGqmlRIPyoTWtTGs0rSSNJw3uPDoZII9u0SSdSaMOAQcnsTFQYGz8\nLttVLhvb9bq150Pdc/jvVXcv7ypX3YKT9ZMsn3vPOfvss8/Zdf9rr7XXlhACHMepJi0LXQHHceYP\n7+COU2G8gztOhfEO7jgVxju441QY7+COU2EuqoOLyAdE5AUReUlEts5VpRzHmRtktn5wEWkF8CKA\n9wI4COBpAJ8KITwHAMPDw+5gd5wm0tPTI/q7i/kFfxuAl0IIe0MIYwD+DsA9F1Ge4zhzzMV08MsB\nHKDPB+vfOY7zBqGtGRfp7+9vxmUc5/eODRs2mPsvpoMfAtBHn9fVv2vI17/+9ejz4sWLy+2lS5dG\n+9auXVtuX3bZZeV2T09PdFx7e3u5LRKbHzy20NraCgDo6OjAyMhIdFyxDwBaWmJBMzk52bC8tra4\n2fS1GS6zKE9EEEKIyud74TrpMjo7O6N9HR0dDY/TYyv8ubjuyZMnsWLFiujafNzExERUBtd3bGzs\nguU3qsf4+Dg0k5OT09o+F932XE6tViu3dZvyPl0PLlPXN3Wf+r3i80ZHR6N9Z8+eLbcPHHhdBA8M\nDAAAPvGJT+Db3/42jh07Vu577bXXojKK623fvn3afTAXI9GfBrBBRK4SkQ4AnwTw6EWU5zjOHDPr\nX/AQwoSI/CcA/wigFcA3Qwi756xmjuNcNBdlg4cQfgzgx1kXUrKWZeGyZcuifYsWLWq4rWUWSylL\n4uVKaH0cX28uptUW5U1OTqK1tTW6tmUqpOqkj+X6s3zUFObRyZMnsXjx4uhYS9ayZOfnAsQSldtq\nPqYj55pVfJxuj0bvTvFcuA30e8uw1J7JM2NWrFhRbh85cqTcrtVqZpnafErhkWyOU2G8gztOhWmK\nmwyYLpFY4mkZlJLlluy0pKAlkViqWVI+dY6+tt7Hddaj3HzffJ6+Fz5O7+Pzurq6ym09SswSmus0\nOTkZlWnJQkuupswDXQ+WltoLknqG+nur/FRbWffFZdRqteh90eVbzzqFPo7f7xMnTpTb7E1aunQp\nTp06VX7Wbc/mgYX/gjtOhfEO7jgVxju441SYptngVvQXuwoAoLu7u+F5lg1u2cXFNkd9NTrOKpO3\nLZtel9coYqogZZ/r8q3rafsxVQ+OgNO2aSrqim1kXUd9Xykb13oueiwgZdPmRqvp8nPR7kvrPWNS\nEYCAPaaSGn/i8nTb63eXIzwt/BfccSqMd3DHqTBNk+jnzp2LPi9fvrzcZvcOEMsTdqtomW9F86TO\n0/LLitxiLHeaJRktV0pKTmp5xuj6W9dmuA14osj4+HhD06VR/SwTgGV+rsQtzIGOjg6Mjo4m79uS\nv5qU2TMTU8G6ltXGuXAZ7BrjySW6HmfOnIk+64lXKfwX3HEqjHdwx6kwTZPoel4tj5zzqDkQyyJL\nhltyieVZKmLsQuWlIrf0XGhrNLWRTAwhTKsHf86V2vq8lOxsVK8crLbS+7gN+Fnr554a9a/Vaklp\nn+vp0Mdapo5FrrTnd3Mm7ZsyIzhHgp4EpM3HoaGhrGv5L7jjVBjv4I5TYbyDO06FWTAbXOcWY9i2\n4W1tm1puLbZVi5lUXV1d03JncVSRjsZqlNdN10nXy3Lp6LGAVA41y9Wm93G7WjOm+Dg9Q4/bhMvX\n5fF9a3uZI+Bmm+QhlcvNclVZ74BlI7N9XjyXWq2Gtra2qK30mEduYg4m19XG7mLtOtb14JlmFv4L\n7jgVxju441SYpkl0LWtZglgpeC3JmEoZDMSScXBwEMBU9Jw2FVJRXEB+dBZfS0vGRnJVRDAxMZFM\nm6zvpVGqYav84hop+Lrj4+PT2j9Vj9wortzkG1yPmbjkuHxd95SLyzKrCvkrIhgZGYneCW3SsVTm\na1lyPddk0eWx+ajroVMxJ8vMOspxnDcl3sEdp8J4B3ecCtM0G1wvvWIlU0zZetrOsWys06dPl9ts\nv+jjUskIdR0tV5g1461RGGtvby+Ghoam5RZPlW/Z/1wvvpZ1L1zf0dHRZHvr8YTUmIGuh/V9bq5v\nyyVnuQa5zlbShUZhrN3d3Thz5ozpEuXny7a6fpY8NnD+/PloX+p56vvnMvVYA7/fFv4L7jgV5oId\nXES+KSIDIrKLvusVkcdFpL/+/3KrDMdxFoYcif4ggP8J4P/Qd1sBPBFCeEBEttY/3z+TC+eu4Ggl\nWmDZYq2maUkdPk7Xg5NUaFmbQkv0JUuWlNuFpOvt7cXp06cjiWe1B8tJXf9U++i24jJ13vnU7CaN\nlS88FVGmnwubRMPDwwCmVpMdHByMorOsvN/W0kJcD2437Q7l9imSLnR3d+Ps2bOmW4v3cTSmrge3\nqZXTPCX5Ozo6ometZ49pt1mKC/6ChxD+H4BB9fU9ALbVt7cBuDfrao7jNJXZ2uCXhRCKldKOAshL\n8eg4TlORnCgbEbkSwA9DCBvrn0+FEC6h/UMhhMgOHx4eLgvu7++fq/o6jkNs2LCh3O7p6ZkWDjhb\nN9kxEVkTQjgiImsADFzohK985SvR54997GPlthUuyjaKNctKuzMarW+2du1avPTSS8lraffFyZMn\ny222u7WNzPXQs+S4HsUyyW9961vx/PPPR/biJZdc0vAcwF4TjI/la2sbnMsotk+fPo3u7u7k+ILl\ngtJwm6SWEgZet7uB19fluv7667Fnz57IzmQ3kE7Yyej2TmXv4bEQoPHy1cVz4Ta1Ek9aa8ZxGfr9\n5gSKjdyv1157LV588UXs2lWOa+OJJ56Iyjhw4AAAYM+ePbCYrUR/FMCW+vYWAI/MshzHceaRHDfZ\ndgD/CuA6ETkoIvcBeADAe0WkH8B76p8dx3mDcUGJHkL4VGLXXTO5UK4MB2KpzPLRmjlkuXesifRc\npq4H5263kj/muix0pFJuwgcrkWBqKV3rOL2sU2oml75nK6KO93Gb6mdmJV3k9uB84frd4cgw/TxZ\nsqfy6wOxfGc5vWjRItMk4jpa98loky41i45NoLa2tujedBtkz1DLOspxnDcl3sEdp8I0bbKJNXFB\nS7/URAMtS3hkVEt03sdyb926dclr6Qi11Mi5HrlNjSADsQxlOdbT0zNNGqbuhaWmlqRcpo6EStVD\n5y1PmTpWogkr5zhfW0tjPo+X3+np6UmugKrhe9H32SjXWqN6MPyu9PT0RG2s5Tq3CY/068lU/O5Y\nE1as1UX5s15N9ODBg4m7ifFfcMepMN7BHafCeAd3nArTNBuc7WDAjgJKubx0GWyjXHrppdG+IjqJ\nrzU4OBit/wQgmRNc15Gjj/RMJ54FpfNVc5mFjdnX14djx46ht7cXjbCSRmjblz/zvek2TCW9yJ05\nBaTHE/Rntov12Au7uNiebW9vj+xWPk6PjeS6FNmW1svt8phHynUHTB9vSbWxbm+O2LMSN1rtz++H\njsRLvTsa/wV3nArjHdxxKkzTJDovFwzEkkbLOJadLP20vGYpqJcgZgnGkxi0e4Gl1PHjx6N9R48e\nLbdfeeWVcpsnhgDAsWPHym2W8gCwcuXKafW99dZb8cwzz0Rtwu67VatWRWWwpLYkuiU1rUi2VD61\n3FzwQPwMrXx7fF4x2WT16tU4ceJE1P7sgtImEbfPvn37on36GRZo04nNveLd+fznP4+HHnoIN954\nY7lv9erV0Xlr166ddh4wXcpz/bVrk0ktS1Wr1aJ3R0t069kw/gvuOBXGO7jjVJimSXSeuAHEElrL\nTi3FC7RM4eP0SOWrr75abhdzZjdv3oynnnoqOo5H31mGA7H841Hda665JjqOZb6O/tq4cWO5zSbF\nxo0bsXfv3obX1pL02muvRYrU3Gst4bQs5+3c1Tr5PCv9tFUGt1VhLm3cuBEHDx6M5Ds/2/Xr10dl\ncJIDDV/7iiuuKLe3b98eHcdmFick6e/vjyT6jh07ovN+97vfldu33npruc1yWtdfv5spL5E2c1L5\nAgA7Mi+6VtZRjuO8KfEO7jgVxju441SYptngHFkGxPaGnhGUWhLGsjv2798ffWZbqXCdbN68GZdf\nfnl03E033VRuv+Utb4n2PfLI65mo2GX25JNPRsexK0WPH/DYA1/72muvjZIEHDp0qNzW98l2q7b1\n2F5nd4yV162wdVtaWjA6OppcElc/F54xpW18PZsqBee507PpePyCZ0/p8RtuY+1SZHuaxzj0c/ng\nBz9Ybu/evbvcvu2223D77beXn19++eXovF/84hfl9q9+9atye9OmTdFxHDlnLXudimRraWmJztM2\nuL7vFP4L7jgVxju441SYpkl0HezPslDLyZRE18dxhJqOVOLIIpbQWhJdd9115baezHLVVVeV21/9\n6lfLbS33uF7WpJfBwcHymMHBwUh6pvLQAa9HfAHTIwJZ5nNklZWjrnDbLF68GCMjI1FbsSvSkvkD\nA3Gm7NRkFu025Agvds/pHHVslujJNxwteMMNNyTrwa7Svr6+6Dh2oXEU4fvf/35zUg2/E1wPNrGA\n+LnrdkwtpaVztfHnNWvWRGXoyLkU/gvuOBXGO7jjVBjv4I5TYZpmg6dmLAHT7cVUzmo9K4dtU21b\nszuJr10s+VLw2GOPldt33HFHtI/tQA5P1Tmw2R7SriW2CXkJntOnT0chuuwGsXJss5sJiO1Trode\n7ie1BFStVouul8rVruuoZ5PxGIiVHIPHCfjZrl69OmoPLkOPSfCMQJ34gJ872//63WH3V+FyWrVq\nFV577TUcOXKk3KevnVpiykrYqcd2UmsE8LXa2tqi+ms3mU4emsJ/wR2nwuQsXdQnIj8TkedEZLeI\nfK7+fa+IPC4i/fX/l1+oLMdxmkuORJ8A8KchhF+LyDIAO0TkcQD/DsATIYQHRGQrgK0A7k9eSLk6\nWIanls7R+3hGFxDLIO2KYEnG0lXPzGIZ/pOf/CRZD3YfabnEslObIizRuR6dnZ1Jea2jwljm6jZg\nd4wlr62ED6lkDTp5RSqHt64/o6Ur11+7L7mNWZ5qmcwSt3A9Nroel6ffMTZL+JzR0dHIDaoj9vi+\n+Z3TpoiVU57rkloNtb293cydd/311yfLZy74Cx5COBJC+HV9+wyAPQAuB3APgG31w7YBuDfrio7j\nNI0Z2eAiciWAWwA8BeCyEEIxGnEUwGWJ0xzHWSAkd5VCEVkK4OcAvhJC+L6InAohXEL7h0IIpR0+\nPDxcFswTABzHmTs4+UVPT880WzfLTSYi7QD+HsBDIYTv178+JiJrQghHRGQNgIF0CfHMHiAOudRu\nhFS+aW1/ctijziLCNmFhU4nItDBTtp20zZlaX0rbZWxv6TpyvQoba8OGDejv709mSNEZQPh6eiyD\n245DX7UrhtuxGCfo7u7G6dOno7biGW88PqGvrevIyQ7ZlcfPCIjDi4u2ufPOO/HTn/40GtvgOulx\nDf5srePGNrKVa764bldXF86fPx89M+2W5GdtrZun16FjuF6NsvAsW7YMZ86cid5N/c4dPnwYgJ3d\nBsgbRRcA3wCwJ4TwF7TrUQBb6ttbADyiz3UcZ2HJ+QW/HcBnAfxORHbWv/szAA8A+I6I3AdgP4CP\nz08VHceZLRfs4CGEfwGQ8mPdlX0hw02mYTnCUlBHeFluG5aohbvk3Llz01xcbCpoOcaSnSPDtAzn\n+uqIKZZ0emYcSzC+T91Wqdl1QCxluQ20dNUzlXibJa+W3gzv08+vUXsDdgQjuwNbWlqS0lvfS8q1\npOuRWk4JSM/2WrRoUXSeNh+tKEPGunbq3df3zKaCvq5+D1J4JJvjVBjv4I5TYZo22USPclsrRFp5\ntVPH6VFM/swyWUeJpY4DYinFx2mJlYoE0+exzBobG4uOTUXeAfmjwVaeO743a1kjHjnX92LlT0+Z\nS/o4lrx8XytWrIi8FtZEEW5TXX4q0YV+to0SLSxbtgxnz541k2WkTAc9yp2bXz4V0akluo7m0+9B\nsr5ZRzmO86bEO7jjVBjv4I5TYZpmg88EtkXYDrFcPxrex3nAdXSWNesnNZvHsq+s5BV8LW0jW3ar\nZYOn7H9tZ1shydYsNIbdZJabhvfp9eTYzuZ76erqisYeUlGEQNx2ur1Tz0a3B7s9i2dURJDxs0gl\nZ9DXtsZlLPs8NRY1MTFhupJzQ8z9F9xxKox3cMepMAsm0a0IJ2YmcjhFIfe6urqmyT0rB3YqskrL\nZC2tGK6/lqS8j2Whll98nuUeSS0RrLEipqx89Yy1jJQVYcjSW5sXLIe5DOte9PvB17PceimzRERM\nF26qXlb5ltlm3Vsq+rDR5xT+C+44FcY7uONUGO/gjlNh3hB50S37gkMKtT1nuclS9qh2i1k2LV/P\nWjvMsn1TLi6deCJ1LcB2x+SOQ+gwyEbX1eh7sVxGqbEMbcfrfO18Tuo8fc/cPtY7wOfp8F9+hjoJ\nouU2TLlwNdY+LoPbTb9H1juXG87tv+COU2G8gztOhVkwN5nlKkjJLi33rLzRqcnyugzep6Ugy/nc\nPNeWuaFdJ6nzLCmsTYpUvWYr5Vky5spfXS9+LtoEYNOEZ/ZNTk5GZeRKXC1d+XOjPHSNysh9LkB6\nJpiurxXJlhs5mOtKtvBfcMepMN7BHafCLNgouiUZU5FbWp7mTsxPTfgAYjmsl9nhMjjpgCXXrVFp\nvpY12USXwTJO1zEl47QsbJRLrKOjAyMjI9H1rOV+rNH8lNS0pDynE+7p6cmWpJaU59Fyfk7WZBC+\nL/1crMQWVqQZY72nKZnf2tqaHSFp4b/gjlNhvIM7ToXxDu44FaZpNniuKwyI7SA+L9eW0efxtuWC\nsmZPcX11GbzPSkpR2G8hhGkRU2w76nqcOnWq3NbL5aZcOtpGY5cU29yjo6PJ2V+5s+usa2vbkZNG\n9Pb2lttLlixpuIwPMH38wEqAmZrtZR2XG9mny8ltq9wEjJpUlJveZ+G/4I5TYXLWJusUkV+JyG9E\nZLeI/Hn9+14ReVxE+uv/L79QWY7jNJcciT4K4M4Qwmv1VUb/RUR+AuDfAHgihPCAiGwFsBXA/alC\nLOmqJQxLYD7OikayyI1MsuRZ7mQZK5KtkIXj4+Noa2uL7pMlupa1J06cKLf1CqgpeaknV/BnnT/9\nQrK0wHoWqQQQuaZTCCFy0Wl3YM619PW4bazEEIyV1x5Iu8k0VrRao3yBmsnJyeg4XV9riamoHhc6\nIExRGHDt9X8BwD0AttW/3wbg3qwrOo7TNCQnO6OItALYAeAaAH8VQrhfRE6FEC6p7xcAQ8VnABge\nHi4L7u/vn/OKO44Trw/e09MzTdJmjaKHEGoAbhaRSwD8QEQ2qv1BRMy/FPv27Ys+8+qOesVP/mzJ\nwtz54Jw22VriyJobbkl0a7SW9xWj4+Pj42hvb8+W6K+88kq5ffr06WhfrkRvlD+sp6cHw8PD0QQQ\nji6zJt/o9khF9+njBgYGyu1iFL2rqwvnz59PSnQ9ip6ayw3EbcDH5Uj09vb2afeh35fUvP3cyD59\nXiOpXaRv5tVudb2K87iDN2JGbrIQwikR+RmADwA4JiJrQghHRGQNgAHr3NzZTED+pHpGN0CjhIa1\nWi3bvZNbP31ta1ZRsT0+Po7W1tbkjCPdifklsHKEp1w/us7aLcbrhaVsWL3PeoGtDtjd3d3wOJ0H\n3Mp9bi1jnBpTyV3DrKWlxRwfSqGfy2zscy5jfHzcbGMrXDq61oUOEJFV9V9uiEgXgPcCeB7AowC2\n1A/bAuCRrCs6jtM0cn7B1wDYVrfDWwB8J4TwQxH5VwDfEZH7AOwH8PF5rKfjOLPggh08hPBbALc0\n+P4kgLtyL6Qli5VYgKWKZZvmunf4mlYkmxUtZNlz1ow3K7EAl892pnadcL10LjeOerMiq1Kzp7Qk\nZbTtmzt7KpVzTF9bf8/X49l7udfSdbaWdc4tz3LvWjKcn4WW0/wu8fty9uzZcnt0dDQp3/V5Fh7J\n5jgVxju441SYpk02mcmILMs6a4Sa0dK70Qjk2NjYtDKs/F6p4yxy3SOTk5NJCWaZItos4RFwKyED\nf9YmQCoRhZbo1gh7asVM7QZKTcyZnJyMrpeacKTP08+Ty7A8N3xecV9FXriUxwFIj15b75XelzsC\nzm1ltYGF/4I7ToXxDu44FcY7uONUmKbZ4JYdYtm+1sR5to+0OyYVsjgT11oqjNVaukiTmrRvhd1q\n+5ntbHYfAdNDUlOwO61wx9RqNSxZsiQKEeVrWQkedf3Z9rWSdHDb6bEWPpbfD/3MUkkXgNiNmEog\nodGzAXOXTUqVobHCnFNRbSEEMwQ6F/8Fd5wK4x3ccSrMgrnJrJUTZxPsbyWNsFa7tBILpOSTlRRA\ny7FUgoNarRbVy5LoPEFj+fI4cc758+fLbY6E0vXgWWKFzD927BhWrlwZJZFgaazlKLedXiWUy+Dj\ntAnBbcplhBAilxo/P/3MrHcnlXzDivwqjqvVamhtbc06FrCXJ7LMFO0ebHSOluj6PnNNM/8Fd5wK\n4x3ccSqMd3DHqTBNs8Gt0ElrCVi2NSxXgQ5VZTsnlSv7QvtS7q/c7C5AOjTTCj3UCR9WrFhRbusx\ng9/+9rflNtvB69evj47j9lm9enX0Pdu41iw/zjDC20CcGJLzuL/rXe+KjuNnODw8XG6Pj48n22q2\nSyFbYzmNwl3HxsawaNGiqA20PZ5a+yw3CaWuVypjUa1WM13EnBHJwn/BHafCeAd3nArTNIn+85//\nPPp82223ldtWtBDLFB21ZMnrlAlguTNyky7q46wcZEwhC4vccI1mNAGxdAWAJ598stxmea2vvXLl\nymQ9WPbzdY8fPx5JdGsmGJfBSyHpenBixf3790fHcdRcsQzT6tWrMTg4GLkA+VlbM7qs9rZcrNY+\nK8+9tUxVCivpBd8LRwOOjY2ZkXg6P34K/wV3nArjHdxxKkzTJPrDDz8cfX7HO96RPJYlDEs/HdFk\nSanUCKSVeMKCZaKesG9FTDWSgoVET5kE69atiz5zTnkdQcaTQ3hbl82SuqhvZ2cnTp06FeWhZ1nI\nkXFAWuYDsfS+8sory+2jR49Gx/FkGX6ebW1tUT2s/HK5S0VZ0WQpb0GRsz5FzkIhunxrH3tFdCQb\n38szzzwTlbF7924AwN13323Ww3/BHafCeAd3nArjHdxxKkzTbPDDhw9Hnw8dOlRuX3fdddG+2SZo\nYFJL2Mw24YMV4cVYUVd6pha7odim1ckOeTaZzovOEWp6WWAmtaxRe3t79JltdW2Dp5IzArFLzXI9\nsg3OUWGrVq2K7jM1G1CXb+2zItka2eojIyNob2833aVWkhHGcqelxgl0Hnd2lz733HNRGXv37k1e\nm/FfcMepMNkdXERaReRZEflh/XOviDwuIv31/5dfqAzHcZrLTCT65wDsAVDoqK0AngghPCAiW+uf\n70+dzAkHAODHP/5xuX3jjTdG+1j6zGRVUiYl8XITSOh6sNlg1cnK16Yj2dhFYuU7Y6wEFamINCC+\n70Iaj42NTcvxxuex6wuwl+1JTe7Rrs1UIodFixYl3Y36Wmxm6WfB7ZNrSunVVtmUys3/rl2n1sQi\nrpc2xxiOCDxy5Ei0z5rMEtUj5yARWQfgQwC+Tl/fA2BbfXsbgHuzrug4TtOQHMe9iHwPwH8HsAzA\nfw4hfFhEToUQimWFBcBQ8RkAhoeHy4L7+/vnvOKO4wAbNmwot3t6eqZFbV1QoovIhwEMhBB2iMgd\njY4JIQQRMf9SvOc974k+cyTbF7/4xWgfSzceWdVy0op2aiTRX331VfT19VnVjNDB/wVWDrmcecBF\ntBTLSa7/0NBQdB7nXdMmhpa5jbZ1+SzROzo6IonKkxh01Jwl0flYvi9dDx45L57n4sWLce7cuSiS\njc0DS2pbEt2Kgmw0H3xwcBC9vb2mRE/lSdOmE+/TP6Lcxvysi/u8+uqrsXfvXjz//PPlvu3bt0dl\nFFGFenRdk2OD3w7gj0TkbgCdALpF5G8BHBORNSGEIyKyBsCAVYhOFrhr165yW3cKfoktezQ3ZNFa\nU8t6EIxecjenTkA6iaG+59msV6XP4z9qlk1YvIwiMu3FtNw77ArTzyX1nHQZqbERK2TYStKhyV1W\n13Kh8b3MNt85l2ktH5z6fmJiInKFaZflnK1NFkL4UghhXQjhSgCfBPDTEMIfA3gUwJb6YVsAPJJ1\nRcdxmsbF+MEfAPBeEekH8J76Z8dx3kDMKJIthPDPAP65vn0SwF3ZF1LSiXN6HTx4MNp39dVXl9tW\nsgZGy6yUO8xye+QuMWPN1LLkWLHd0dGB0dHRpPyzcrxbObYZa2mkwh4NIaClpcWcdZUq00qwoWeJ\nMak88SGE5BLKM3GVppJo6DZlc0NHN2p7neGxGGt5X762doWllnni9+i1117Dyy+/XH7W75WOJEzh\nkWyOU2G8gztOhWnaZBM9csguL3YHAGmJrmWKFbnFWMsCsTyzZJYVaWZNfphNQgktEfmzlVPOkoyN\nogM7Oztx/vz5SO5Zo+jcVlYEFqf01c+MP+tIMB7RbxR51+g8i9SkDiA9aaSIMizQ7Z07wp4yRRqV\nWVCYqps2bcLBgwejSDbLXWfhv+COU2G8gztOhfEO7jgVpmk2uHaXsN334osvRvve9773NSzDsj81\nqaSLuTaxxlpqyXLHNAqJnJycRGdnZ7ZLh9G2bypSTkc+cZ2LpYWuuuoqHDlyJJnIwXLF6HECfr5W\nsgauL5/T1tYW1dGy8S2XJT93q/5W21v5yFM2uJU/XUcLchncHi+99FK0zUku9X3mvsf+C+44FcY7\nuONUmKZJdCtRgY5kY/cA5yDLdTfoY3Mjjqz8XrkTHKzcX4VkPH78OFasWJGMirIityzJaCUZ4Gtx\nFNfZs2cjOc/SWyd8aDQjLfW5UXlAWsp3dnYml4fKzTEOpNtgJktWsalgRQRaCR9SphOQXqaKl3na\nv3+/mfhG7tbaAAAQNElEQVTEahPGf8Edp8J4B3ecCuMd3HEqTNNscG2bsl1SLCNbwGtx8TpXmtQa\nZjPBCilkuP5WthTLftb52VP1t5JS5CaN5HXKLLTtbGWFsZIkpMYo9PfcBuw21PnIrTGP3LBktoP1\nfaZs9RBCdgYZy1bnMQ89tsPnHThwoNzmsaeBgYGojfX7bc14i+qbdZTjOG9KvIM7ToVpmkS38miz\nqwCIZ5dt3ry53NaRSZzojxP2AelINo21JFHqPJ0gwIrAYpnI5509eza6H3ZdcQQTELsYOQGjPs9C\nn5eCy7PyxOnr8r1xDnwtTxvJ8NHRUbS3t0dSmdsmNy8fkHZdaTctz2bUiSZSS1Zp+H3Rx1mRj/yZ\nkzrwc9fvwGzxX3DHqTDewR2nwjRNolujy3pEkCefHD16tNzWkVUs97ScTOVM1yOwjSaDFKRGr/W1\nOBKsmMhRcPz4cWj6+vqwa9euaGSXzRRdhiU1UznTOZ+8LqNog2uuuQaHDx+etnRPwUwmOHBb8Wiw\nHs1nU6qQ8q2trTh37lzUHmwCWCu2WpGD1vuRyvE+OTlpvi+piSIabg9tznDuNW4rnXjCmrQzp0sX\nOY7z5sQ7uONUGO/gjlNhmmaDazcC207a/fXCCy+U27xwoY5qs6LQUm4yK2GitjlTNpx2YfA4gZ4Z\nx/sKW6yvrw87d+5Mumq0nZ1abwuIXYVc3xMnTkTHNbId3/3ud0dJBoB4nGMms7gYbm/tUmQ7e9Wq\nVQCmXKFPP/00VqxYUe6zcquzG06PNaTOs2Yi6vcolWxTf7aWvbKSb/B4C9vnenmsuUjm6b/gjlNh\nsn7BRWQfgDMAagAmQgi3ikgvgG8DuBLAPgAfDyEMpcpwHKf5zESi/2EIgXXfVgBPhBAeEJGt9c/3\np07OzeMMxLLuscceK7ff/va3R8exVLOC73OljrXPco+wHOvt7Y32sQxn90hfX18kz1juadcS19+S\n7yyNrYg6PqezszO5TLK+TytqjuvPMly3Kbdj4Q7cvHkz+vv7cejQoYZl8HMGgNWrVze8LgCsXLmy\n3NbmAZNKBjE5OZmd09xKyMBLBPNzB9LJN/S7aMlwa0IMczES/R4A2+rb2wDcexFlOY4zD0jOL6uI\nvAJgGFMS/X+HEP5aRE6FEC6p7xcAQ8VnABgeHi4L5oEyx3Hmjg0bNpTbPT09037ycyX6O0MIh0Tk\nUgCPi0i01lAIIYiI+Zfiox/9aHxhY74vjwxfccUV5faXv/zl6Li+vr5yW0e5NRrJPXDgANavX5+8\nrhUVxRFjPDIOxNFqWsbyeYVUu/nmm7Fz586kRNfMl0T/9Kc/jYcffjg61orYm2uJXjyzz372s/jW\nt74V5d/LlejFSHxBSqLr9uX3r5C7IyMj6OzszM61xu+pluj8TuhoxqGh14eqvvvd75bbhw8fBgD8\n6Ec/woc+9CEz3XdRZ73slyarg4cQDtX/HxCRHwB4G4BjIrImhHBERNYAGLDKmMkSsOw2e+6558rt\nX/7yl9Fxlr3FL4iVAM8Ke0y5grRbj9fistw2/JLddNNNUUgqb1s5wa0XjtEvNP8h4O1LL700ahPu\nZLqt2K7UbZPKd67h87htent7kx1cj0lYdUwlWshdFrlWq0VlWu3IfwiscFTL1ZaL9QfK4oK9TkSW\niMiyYhvA+wDsAvAogC31w7YAeCS/uo7jNIOcPwOXAfhB/a9hG4CHQwj/ICJPA/iOiNwHYD+Aj89f\nNR3HmQ0X7OAhhL0A/qDB9ycB3DXbC+cO87P99eCDD0b77rjjjnKb3VFALKNZnlkT83NdaFqG87W0\nvG40O2toaAgrVqyIIrdSch2I7XhLkrKMs9qXj1u7dm1kS1rJMbiNtWnAZXJ7WFKSZfj69eujcRTe\np12gvC83d9tMcurnJvBgdNISfmZWHa2IulSe+JngkWyOU2G8gztOhfEO7jgV5g2R0UXbxWxvsO2r\nZ0i9+uqr5Ta7ToB0ZhJrbS8r2wvbgdpOZZtQzxyysn6wrcr+XO0W4jJ1W7H9zMfp66ZcM0uXLo3a\nKtc3r8c8uE34WWiXIl+Lj+vr64vaOHfpXysc1fKDc3sUdraIYGxsLHKFaRs8lQvdSpKo3xf+zO+j\nDhnm9tH9pxmhqo7jvMHxDu44FeYNIdE1LH0sCcaRVVq+s8RjOWlFsuVGuVnJK7SpwBJal59aMljf\nM7uPtCuP740j6nRkFdeZ9y1btiwqP7WkMZDO8a7rlTKxgLT7q6urK9pnJeJIXUvX0TK/mEJ2t7e3\nY2JiInJxWe5ANom0lOf3zwqBTi0H1draauZWn7NINsdx3rx4B3ecCtM0ia4lRSovNRDLKZauenTZ\nilDjFUs5YkzPxrJGWlOy3Mr/Zq2myfepV9O0Vu7MHem3ZnGlcrwvWbIkOWtpJvVgKcvb2jTjOmqz\nJDW6bC0VpSU0S2UrJxufV5gs7e3tGBkZMctn+c4Rh/q5873pfaloQS3RGcvTZOG/4I5TYbyDO06F\n8Q7uOBWmaTa4xrJRUraptkN4Ur22sdhWKuxuEZlmg1vL1LIbx1ofjOurZz6l7K2WlpboPCsJgFVH\ny63F6CVyi3InJiaiMq31sFLJK/SxfM/aBm8UoVar1dDW1hYda83y42et3VOpbCxWBGBRxrJlyzAy\nMhKVqZdd5neOy9cZhZjc91vfM7exHofwSDbHcbyDO06VWbBINpYfViC95S6x8qnxeUXEW3d3t+km\ns8pP1V0fp02FVNKBEEIyQYMuw5K8LJs5Qs1yp1n5t61c34yVHCOVbONCpGSnluF8n1p65y6ry9Kb\n20pHsmmJztGTnIyE3X9AbMZZ8jrl6tX3pdvRly5yHMc7uONUGe/gjlNhmmaDWzOpcrHsOcvOYXtI\nT8y33HUp95E1u8lC275cvtUes8mjbWEtl2slbrRCKVPhxXosoJE7rVarobW1NZm/PqeMglQZOhFH\naq250dHRqA20Dc7X45l8+t20lsfmY1PjJiGE6DwdMmvNsoyulXWU4zhvSryDO06FWTCJbkljlki5\nOaW1DGJZbkkdnhGk852nlr6xXH6W+0i7/1LRWlb+MMvNZ2GZANwm1tpkXA99nyl3T67ZMz4+nnRx\n5Ubo6c9WbrVUm46MjETRarq92TVmzVazXIV87dS7o9vXShph4b/gjlNhsjq4iFwiIt8TkedFZI+I\n3CYivSLyuIj01/9fPt+VdRxnZuRK9L8E8A8hhH8rIh0AFgP4MwBPhBAeEJGtALYCuD9VgBWtZu2z\nJmRYub/0pI8CnTTi5MmTDa8LxJKdJZiOWrJGjXmUlGWWFXGl91kTSlJpjrW0TE280JNNrGQH3KaW\nKWKlqWbpWlyrtbXVnLxirWSqz0vts9If68g1noiic+yllk2yvArafEy1sY6qtHLP5UYI5qwu2gPg\n3QC+Ua/4WAjhFIB7AGyrH7YNwL1ZV3Qcp2mIlW0SAETkZgB/DeA5TC1CuAPA5wAcCiFcUj9GAAwV\nnwFgeHi4LLi/v3/ua+44DjZs2FBu9/T0TBt5y5HobQA2AfiTEMJTIvKXmJLjJSGEICLmX4rPfOYz\n0WdrBJxlJ0sTHXTwhS98odzetGlTtK/RH64lS5ZMk2os0bUcY4nO0syS6HpUm6VyIeMGBwfR29s7\nrX4FM5HoeqJEo+sCjSV6d3c3Tp8+nQykmQuJro/jZ13UvbW11Vz11VrhM1eis+zW+4r36pZbbsGz\nzz4bTSix3gkOdLEkul4FZseOHeX21772tWnnPPbYY/jIRz5ijqoX79yuXbtgkdPBDwI4GEJ4qv75\ne5jq4MdEZE0I4YiIrAEwYBViTfzXD5A7ECdP1B0rNzmBNeuMH5iOckslTNR/PLheuoPzHyudMCHl\nAtTlc/voTpdKiGHl8+ZzxsfHkwkTrRl6VuSW5eZslPO+6OA5s6x0Ha1kityp9R88blPtJuNnqMds\nUu4vy11pJaxIJcMEbHfjnLnJQghHARwQkevqX92FKbn+KIAt9e+2AHgk64qO4zSN3FH0PwHwUH0E\nfS+Af4+pPw7fEZH7AOwH8PH5qaLjOLMlq4OHEHYCuLXBrrtme2Er4QPvGxoaKrcvv/zy6DhrBUqm\nkMZjY2PmkjsajmjiOullgdietqSUnkzAEiyVuw2ITQArcisVeafL13I9tUyQ5bqybN/UqqxALDsL\n+7ajowPj4+PZbkSur34WbGZZK6WmovJEJMqvpm1wbh/e1u+iNjtT1049P/1crFz8Fh7J5jgVxju4\n41QY7+COU2EWLC96KhwVAI4fP15us42iXRbabcakZn/p/NUpnzsQ21GcMEDby1zHlStXRvtSrjwR\nia5nuYX4PJ7NBKRDbfXYQsq2tu7FGhfQdWQbkcc59DNKza7Tn61kGOzWGh4ejvaxDW4FcXEd+Rkt\nXrw4es/0O8H3k1rTDbBDSXOSbepEHBrPi+44zoVDVWcLh6o6jjP/NApV9V9wx6kw3sEdp8LMm0R3\nHGfh8V9wx6kw897BReQDIvKCiLxUz/zSFETkmyIyICK76LumppkSkT4R+ZmIPCciu0XkcwtUj04R\n+ZWI/KZejz9fiHpQfVpF5FkR+eFC1UNE9onI70Rkp4g8s4D1mNd0aPPawUWkFcBfAfgggBsAfEpE\nbpjPaxIPAviA+m4rptJMbQDwBNS89nlgAsCfhhBuALAZwH+s33+z6zEK4M4Qwh8AuBnAB0Rk8wLU\no+BzAPbQ54Wqxx+GEG4OIRTzLBaiHkU6tLdiKqHKnjmtRwhh3v4BuA3AP9LnLwH40nxeU13/SgC7\n6PMLANbUt9cAeKFZdalf8xEA713IemAqn96vAbx9IeoBYF39pb0TwA8X6rkA2AdgpfquqfUA0APg\nFdTHwuajHvMt0S8HcIA+H6x/t1BcFkI4Ut8+CuCyZl1YRK4EcAuApxaiHnVZvBNTiTkeD1MJPBai\nPf4HgP8CgEOxFqIeAcA/icgOEfkPC1SPqwAcB/A3dZPl6yKyZC7r8Xs7yBam/jw2xYUgIksB/D2A\nz4cQorQxzapHCKEWQrgZU7+gbxORjc2uh4h8GMBACGFH6pgmPpd31tvjg5gynd69APUo0qH9rxDC\nLQDOokE6tIupx3x38EMA+ujzuvp3C8Wxenop5KSZmgtEpB1TnfuhEML3F6oeBWEqI+7PMDU+0ex6\n3A7gj0RkH4C/A3CniPztAtQDIYRD9f8HAPwAwNsWoB6N0qFtmst6zHcHfxrABhG5qp4N5pOYSvW0\nUDQ1zVQ92+w3AOwJIfzFAtZjlYgUGXC7MDUO8Hyz6xFC+FIIYV0I4UpMvQs/DSH8cbPrISJLRGRZ\nsQ3gfQB2NbseoRnp0JowmHE3gBcBvAzgv8739ei62wEcATCOqb+U9wFYgakBnn4A/wSgd57r8E5M\nyavfAthZ/3f3AtTjJgDP1uuxC8B/q3/f1HqoOt2B1wfZmt0eVwP4Tf3f7uK9XIj2wJRX45n6s/m/\nAJbPZT08ks1xKszv7SCb4/w+4B3ccSqMd3DHqTDewR2nwngHd5wK4x3ccSqMd3DHqTDewR2nwvx/\nRbO63Ql6l9cAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e919766a0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "@interact(n=(0, len(faces)-1))\n",
    "def display_faces(n=0):\n",
    "    plt.figure()\n",
    "    plt.imshow(faces[n].reshape((64, 64)), cmap='gray')\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Mean and Covariance of a Dataset\n",
    "\n",
    "You will now need to implement functions to which compute the mean and covariance of a dataset.\n",
    "\n",
    "There are two ways to compute the mean and covariance. The naive way would be to iterate over the dataset\n",
    "to compute them. This would be implemented as a `for` loop in Python. However, computing them for large\n",
    "dataset would be slow. Alternatively, you can use the functions provided by numpy to compute them, these are much\n",
    "faster as numpy uses machine code to compute them. You will implment function which computes mean and covariane both\n",
    "in the naive way and in the fast way. Later we will compare the performance between these two approaches. If you need to find out which numpy routine to call, have a look at the documentation https://docs.scipy.org/doc/numpy/reference/.\n",
    "It is a good exercise to refer to the official documentation whenever you are not sure about something."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "__When you implement the functions for your assignment, make sure you read\n",
    "the docstring which dimension of your inputs corresponds to the number of data points and which \n",
    "corresponds to the dimension of the dataset.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def mean_naive(X):\n",
    "    \"\"\"Compute the mean for a dataset by iterating over the dataset\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean: (D, ) ndarray which is the mean of the dataset.\n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    calc = np.zeros(D)\n",
    "    for n in range(N):\n",
    "        calc += X[n] # EDIT THIS\n",
    "    mean = calc / N\n",
    "    return mean\n",
    "\n",
    "\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def cov_naive(X):\n",
    "    \"\"\"Compute the covariance for a dataset\n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    covariance: (D, D) ndarray which is the covariance matrix of the dataset.\n",
    "    \n",
    "    \"\"\"\n",
    "    N, D = X.shape\n",
    "    calc = np.zeros((D, D))\n",
    "    mean = mean_naive(X)\n",
    "    for n in range(N):\n",
    "        diff = np.asmatrix(X[n] - mean)\n",
    "        calc += diff.T @ diff # EDIT THIS\n",
    "    covariance = calc / N\n",
    "    return covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def mean(X):\n",
    "    \"\"\"Compute the mean for a dataset\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean: (D, ) ndarray which is the mean of the dataset.\n",
    "    \"\"\"\n",
    "    mean = np.mean(X, axis=0) # EDIT THIS\n",
    "    return mean\n",
    " \n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def cov(X):\n",
    "    \"\"\"Compute the covariance for a dataset\n",
    "    Arguments\n",
    "    ---------\n",
    "    X: (N, D) ndarray representing the dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    covariance_matrix: (D, D) ndarray which is the covariance matrix of the dataset.\n",
    "    \n",
    "    \"\"\"\n",
    "    # It is possible to vectorize our code for computing the covariance, i.e. we do not need to explicitly\n",
    "    # iterate over the entire dataset as looping in Python tends to be slow\n",
    "    N, D = X.shape\n",
    "    covariance_matrix = np.cov(X, rowvar=False)\n",
    "    return covariance_matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the `mean` function implemented, let's take a look at the _mean_ face of our dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAAD3CAYAAAAwh5neAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnW2MXdWVpt/l+qCq7KLsMpg42MIMFG0hGJMWygAhre7Q\ntOhMBCiKECg9cUYoEyUziJEYDU5PNFL/GA1/0ppW0hpNCxIcNQ2DaDJYqKc7bjo9E0QnBBLAJjgU\nw4ewY1wGV5VdZdcne37UvZd1Vp311r7XVffCyXoky+fcfe4++3zsuu/aa+21JaWEIAiqybpONyAI\ngrUjOngQVJjo4EFQYaKDB0GFiQ4eBBUmOngQVJiz6uAicpOI/EpEXhORPavVqCAIVgdp1Q8uIl0A\nXgVwI4DDAH4G4I6U0i8BYHJyMhzsQdBGhoaGxH52Nr/gnwTwWkrp9ZTSHIBHANxyFvUFQbDKnE0H\nvxDA22r/cO2zIAg+JHS34ySjo6PtOE0Q/MYxMjJCy8+mgx8BsF3tb6t9Vsr9999f2D///PNLtwFg\neHi4dHvjxo2F4zZs2NDY7u/vL5T19vY2tnt6egAAv/71r7F9+/bCcV1dXaXbALBu3brSMv253Rcp\nmkF6v7794osvYteuXcuOzYGNmbz//vtNfe/gwYO44oorstuRe5w+lz2vbmN9+9ChQ9i5c2dpmd1e\nqWxxcbF0e2FhoXCc3q8fd+rUKQwODmJ+ft79ni6bnZ1tbM/MzBSOO336dGP75MmThbLx8fHG9vHj\nxxvbJ06cAAB8+ctfxoMPPtjYt9/Rx/7whz8E42wk+s8AjIjIxSLSC+B2APvOor4gCFaZln/BU0oL\nIvLvAPwdgC4A300pvbxqLQuC4Kw5Kxs8pfQ3AP4m59i6TG6cuLu7dNseq8uYNM4ts8cx6e2Vna1E\nr3+nFWnMJLpuhz3O+15XV5fbjmZMCK+NVkJ798OaR+zcun52vzXsmek2isiqvFeeeQcU32ltSto+\noftBM/en0I6so4Ig+EgSHTwIKkxb3GTAcomh961E9yQvkyW5Uo3Ryihx2f7Z1p8r81m7ctvEpJ89\nF5PG3vlyryX3ulYq8+6BHlEHirLcfscrs/vZMpmYfkyiM/MxJHoQBNHBg6DKRAcPggrTNhvc2hDa\n3mARZK3axV60k7XFvO/YdrTiCmPfm5ubyx5rsO41r4yR6wprZpxD49m+LApNR4ktLCycVRRa2b7e\n1hFotg7dptnZWbeNK53bI9eFxmxwO05l991zZx0VBMFHkujgQVBh2ibRWSQbkzAMLQWtlPLK7KSA\nXPmbazbkSvSpqamWou2Yu5G1wytbXFx0zYNcFxTgm0FMQtvnomW0t2332WSTXClfr6OnpwfT09OF\nMnud3rnZ+8fcfN47tm7dusJ+SPQgCJYRHTwIKkzHJDqL0vFGYa3M0nIpd0RzcnLSLWsl+m2lOjwJ\ndvz4cVd6W/lVNre9bJ/JfI0+1/z8vHvdTKLnSm8mr/V86vHx8cL+3Nxc6XG2jEWolcnwsmups2XL\nFoyPj9N59V6dzGSx8t2jmUg2zzSzxC94EFSY6OBBUGGigwdBhWmbDW7tSmZDeLaetcWYO6PMPhoc\nHMS7777rHmfxbCwWnWUps7OHh4dx+PDhgv2st88555xCHX19faXb9liWIMCz92dnZ10bnLmgWGQY\ns5+1m7Ket6yvrw/Hjh3DmTNnSo+zrk2938wsMU3Z2MiWLVtw7Nix7Bl1Xn0rnVvfO11mx1BWw5Uc\nv+BBUGGigwdBhemYRGf5w7Ts0nKPuRtyJh0MDg7i2LFjheO8SQe2jEVn6e8xd0b9HgwPD+Ptt98u\nyGsmw9evX9/YHhgYKJTpdNH6eyxyULfp9OnTrtxjEy30cwF8Sa1lNwBMT08Xzg0AO3bswDvvvFMo\n09+zdeS6yby8a0DR9ajvx4kTJ7ITUWhyJTPgmzpW1rNItnCTBUEQHTwIqkx08CCoMB1L+KCx9pG2\nsaytp2GuCM+tderUqcJxbNaSTUjgHcfcZGW5v3ft2oWjR48WbHBtS1s7W9vgehsoLt+ky6wdr21y\nbX/OzMy4dqa9Tu3ysu4vbT/re2zt56mpqWXbO3bswNGjRwvL/XguM4C/ExqWVES3Xx83NTVFZ+9p\nchNbsjBZb1xjbm6ucFzuLEJL/IIHQYVZsYOLyHdFZExEDqrPhkVkv4iM1v7ftLbNDIKgFXIk+oMA\nvgPg++qzPQCeSindJyJ7avv3skpYpI91x3gzcXIlkT2f/p6Vd2wlSW+GlK2DudD097RMPnXqlOta\nYpFb9tyeK8/KQi3Z9f2YnZ11XZZWojP3l5bl3jZQlPJark9MTLgSnd1vizcrT5slgG9WLSws0CQg\nOXUA/D56rjEbAZh7nYwVf8FTSv8XwAnz8S0A9ta29wK4NetsQRC0lVZt8AtSSkdr2+8AuGCV2hME\nwSoiObJXRHYAeDKldEVtfyKltFGVj6eUCnb45ORko+LR0dHVam8QBIqRkZHG9tDQ0DJ3SKtusmMi\nsjWldFREtgIYW+kL+/fvL+xrt5ANw/Psymbyp5fZ4Ndeey2eeeaZwnG5M9J0m1joZI699aUvfQnf\n//73C9et7XN9b4Ci+2twcLBQdu655za2h4eHG9vafQYUXW91e3RmZgZ9fX2Fduhrtravvm5tLwNL\nNnTZ9smTJwvHaZu8XscXvvAFPPbYYwXXlbb37dgIs5G9ZamtDV523M0334x9+/bRpa01rbw7QPHe\nleXvr78f7J2rZyb63ve+57YPaF2i7wOwu7a9G8ATLdYTBMEakuMmexjAPwH4LRE5LCJ3ArgPwI0i\nMgrg92v7QRB8yFhRoqeU7nCKbmjmRGzGGHMt5boscuW7jfDy2mT3tZRiLj9bh81Bruuzcr6OvRYm\nNb3kCjYaLjcRArvHbDkh3Q4tJ7VbDPBdfnNzc9n3m0Wa6fvDZtd5ySr7+voKdbJn7S2PtVKZxnvX\nWYIKIPKiB0GA6OBBUGk6NtmESbDcOlgecC86KzfwH/Dlrx1B1iPDbGKEnfSir1tLS+a6ZAkldB12\nJF5PbGHykZ1bH2uv0xtht/eKSXR9j71oL4BLdH1/9LmYhNZyfWFhgSbw8OR27nG2zVpqe/nZcvY9\n4hc8CCpMdPAgqDDRwYOgwnTMBs9d3lZv587eAYr2tLbtbF50/T1rL46Pjze29Zpm9jhvTS3Aj6w6\nfvx44ThtP9totaGhIXjo+6rt7mZmYDH3jIa5A701x1hyRpskQtvxbIlgOxtO47lVra2u79WmTR9E\nWR85cqQQOWhz1OtzlyXzWKkdZcfWYVGbrRK/4EFQYaKDB0GFaZtEb8b1kxu9xpIweHnB3nnnncJx\n+ntakgPAa6+9Vvo9Fq1m+fjHP97Y/tjHPtbYTim5UtNeM5uI4uWNY+aMdS96UX+5yx9ZvMlCtszm\neNOmj5cow5a9/vrrhTI9uUVfi5Xyuv0XXXQRAOArX/kKnnvuucY+sDwHnm6LNrnsc9H7ue5dXXdf\nX1/BnAk3WRAEy4gOHgQVpm0S3cpYPWpsy7xIIjZya+cdv/fee6VlR44cKRxnJZhG5wzT8p3lO7OT\nQfR87QsvvLCwrdulzQg7kcBb4qjs2DpsEkn9O7Ozs+ju7nblto3sY6uXetKeLaFkV1fVbdT3UZs5\n9loOHz5cKBsb+yA1gX4/bOSdvj/aSzExMYFt27Y19u1SV/rZ6zn31tOh74Gd+OPdb92m7u5u1xxg\ndVjiFzwIKkx08CCoMNHBg6DCdMwGZ0ujem4ymyCBLS3sLX9kI5O0XWztGh31pt0sNomBdols3ry5\nUHbZZZc1ti+55JLG9qWXXlqw73SknHWBaDt+48aNhTJ9PXqbzbKy0XXa1mNRbXqcgCWl0DZnbqIC\nG72nxxq0TQwUx03smIq2tfU4jB0b0XXo53LJJZfg4osvduvXNr4eo7HvlR4vss/Tm9mn6e3tLdw7\nO4ZgxzY84hc8CCpMdPAgqDBtk+hWmrEoIC2ntBS00Ugam1ZWSyQtcS+99NLCcVqiWxmnTQI9IcFK\ndC2XdOpiAIWoKC1dBwYGCnUyV5uW5Vai63N7ct3Wabe9ZY2sfNR1sjTEZSmaVyrbsmVLwe3kucyA\n4r26/PLL3TbqCT3WvNP3Xkv0yy+/HBdc8ME6Hkzaa9PPusJ0O6yc1mVefrxzzjmn8O5b92juCqvx\nCx4EFSY6eBBUmOjgQVBh2maDMztE28hA0SbXrgJrdzCb07qr6uzYsaOwr20q28adO3c2ts8777zG\ntnVZaHvRtkPblTYEV59bX7Mdk9AuJLskkXduW4dn4/f29hb2dRttaDDLz67rZ7PatJ2p69i8eXPB\njmU2pi7bsmWL2w49o9C6UfW59LjJtm3b6Huln4Vuh3WFsaQR+j3z7ndvby8N+bX33yN+wYOgwuQs\nXbRdRH4kIr8UkZdF5O7a58Misl9ERmv/b1qpriAI2kuORF8AcE9K6eciMgjgeRHZD+DLAJ5KKd0n\nInsA7AFwr1eJlRRaQloXg9735AxQlF3a7QEU3SL6uPPPP79wHIvc0jOEtMyy8jE3z7idceWZB/Z+\naNlppbdddqeOnSWn77+NZMvNY8ZWc/UkupXGXu7vgYGBwvm8fPIAjxLT5p41ZzSeubF+/fpCmXVP\n2fvvtYNFFXq57bTp19vbW7h3tv+smkRPKR1NKf28tn0KwCsALgRwC4C9tcP2Arg164xBELSNpmxw\nEdkB4BMAfgrggpTS0VrROwAucL4WBEGHEJYrrXCgyAYA/wfAf0kpPS4iEymljap8PKXU0MmTk5ON\nikdHR1exyUEQ1BkZGWlsDw0NLUsOmOUmE5EeAH8N4KGU0uO1j4+JyNaU0lER2QpgzK8BOHDgQGFf\n27fWftb2o7aBcpe2Bcpzc/f19RVcJ0DRRrThjHpf24T2OG/ZW6Boq9Ztu8suuwyvvvqqGz5q3SrM\nntP16+Os/ant+nrZG2+8gYsvvrhg/+s/+HasQWeg0TO1gGLGGz3Lio1J1O/byMgIRkdHXdvd3u9c\nN5zGjhmUZUvZtWsXXnzxxYINznKa6+PY7Dp7bi/xZN0Gv/LKK3HgwIGCTT4xMVGoo/4e33uvO+y1\n1H5autQ4AfAAgFdSSn+qivYB2F3b3g3giZXqCoKgveT8gn8KwL8CcEBEXqh99scA7gPwqIjcCeAt\nALetTRODIGiVFTt4SulpAF7i7xtyT8TyOjN548lYWweb+VSXeLOzs8ui5tgyO1oa6nOzZYDsdXpJ\nBrdu3eoua8Su00pQ3RbtwrF1eDnNu7q63PtopaW3VDFQvLZcF46uf3Bw0J1ZZV1t+rkwc8lzUQL+\n/RgYGKDJDr3nyUwn+8y8KD07m4zlhvcSZ1giki0IKkx08CCoMB3LycbKPPnUjEzRx9Zl20oSnS2z\nw6LV2HI/Xj7y8847r3SE3X4H4CPbnrS37fDyrtl7z8q8EWRbvzcqb+u0SzJ5I88sks3Kd+bR0HjL\nB23atIlep9dGlqecmXS6TJsbXV1drglnz82IX/AgqDDRwYOgwkQHD4IK0zYbnEWhMTst153Gltz1\n1pMCipFy1lZidrfXXmbT1tt4/PjxZWMBzKbSdjfLM86W981dbjb3Wmz92v73ZvKxOrq6ugrPwhsz\nAPzZWID/zHLfsU2bNtHlq70ZdWx9PdtGva+frb0fue1gxC94EFSY6OBBUGE6JtFZFJoXTdXM8kce\ndgI/c6t4kVUWJq+979kJJRobnaXrYGaKLst1fzXjvtSwiD32LLzoMmZ6WPdorjxlEr3suJmZGaxf\nv57W790TFmFoyzxTyrruWCReuMmCIIgOHgRVJjp4EFSYttng1r7NTd6Qa0fZ+jwbxdqH2gZirrDc\nzDf2OG+tr66uLtfdY11LbIZXrktxtWFhrLk2rB3jyHVLskQLbBaXpqysboPnjr2wddy8NjHYeXPf\nP0v8ggdBhYkOHgQVpm0SneXHZlFA2mXEIoKsRPIkKpNLTNbmSjV7HHP5ebm5mpHoHrnHWVo1RbyZ\ncdbl55lcKSXXTLFmlV0CymsHy+OusW7IXKnMZDkry2lHSin7vWLEL3gQVJjo4EFQYdom0XOTKQD+\nZAIrXXOlIIua07Qy2mn3WVSezSV25syZxr6W6Ox87Ny57bDb3veYWcXqz4kaK/tcn4+t3Om1CeCR\nYR7NmDPePbDvdyumTu4SWABffVUTv+BBUGGigwdBhYkOHgQVpmORbNpmYZP22eR+vZ/r4rJ2pYZF\niTE7mF2Ltru1nT09PV2wwdmSuKyN+rpz69BYW05fG1syiF0nm0mV+2x1fc0kU/BylbPZjDapJbPJ\nvXPnvsOsPuY2ZEtWM+IXPAgqTM7aZH0i8qyIvCgiL4vIn9Q+HxaR/SIyWvt/00p1BUHQXnIk+iyA\nz6SUpmqrjD4tIv8bwOcBPJVSuk9E9gDYA8Bd6tDKHub+8qRgM/mxtTxjUjA3wonV4bl37L7enp6e\nLuyzpA4sN5dXB4vs0+2dmZnJdv14q60CRfOD5Zr3coQvLCy40Y0sAQYzAZhE95JjLCwsZL8TzK3F\nVkfV9847bn5+Pjv/O2PFX/C0RH092J7avwTgFgB7a5/vBXBr1hmDIGgbkuOQF5EuAM8DuBTAn6eU\n7hWRiZTSxlq5ABiv7wPA5ORko+LR0dFVb3gQBEtrq9cZGhpaJj2yRtFTSosArhKRjQB+ICJXmPIk\nIvQvxQsvvFDYHxwc1A0rlOnUxnpbL2IPFPOa2bxdZfLsrbfewkUXXVQ4rhMSfX5+Hj09Pasi0b2l\nndi89/r22NgYtmzZ0pJE1x4AYMnkqNOsRN+yZQvGxsayr4Xla/PmiudI9KmpKWzYsKElid7MCqin\nT59ubE9NTTW2T548CQDYuXMnDh06VCg7ceJEoY7x8XEAwDe/+U23rUCTbrKU0oSI/AjATQCOicjW\nlNJREdkKYIx9l4XeMRuc5dguW3+s7Hz64bI83bkhqM08TGZj6X3m9mDJFD2b0yZ19BIhnjx5Mjvc\nlXVw/dLOzMy4dej69XOZnZ1d1lnr5Ca5sPutzKhjblQgP4yajVfofe/dse/HmtngInJ+7ZcbItIP\n4EYAhwDsA7C7dthuAE9knTEIgraR8wu+FcDemh2+DsCjKaUnReSfADwqIncCeAvAbWvYziAIWmDF\nDp5SegnAJ0o+fw/ADbknsjJFSw4m0Zl9y6LXvLJcaQPkRy2xaymTWd3d3VhYWHC/x1x5zFWoj9NL\nCdXPaY/r7+9v2HJ1mES37jWNlujaBmdLLWs5PT09XTAr2HNiyUI0+twsaYSug5l6dj/X1bsaEp3N\nlmREJFsQVJjo4EFQYT4Uk01yJTqTv2wFSi0Fc6WfbbOXuMHus1F0K9F1W7xJKbbNLKqL5Sorc5Nt\n374d7777LjxyV8UEipJdt5/lbtOS3I7K5y7b04qb036PeVnYxCLmBdH3xz5PL7rRbnsRb7YdjPgF\nD4IKEx08CCpMdPAgqDAds8HZbBvPpmXhl4xmbCyNbjOzh9i1eO61xcXF7FBYlpzRO7e93/re1e/H\n9u3bcfz4cfceMJdRjjsQ4MsCa2ZnZ90y+7munyW2YGG3mvp7ZMOHy/DcmfZ+eLa13dfPM2zwIAia\nIjp4EFSYji1dxCK3tDTxcmwB+fmy9feYBGNuISbDWY53zx2zbt06163FJpS0mmhBo+/Ne++911R0\nX1kdgH8PrLz2nouN3GKrvrKowtyVQW2SB2DpXWP56W1bWLQac3t6z8xK9NxluxjxCx4EFSY6eBBU\nmOjgQVBhOrY2GZvM7tkvbKnY3OR71h7KnT2VmwCP5S23IZbaPmVhpl57bVvYEsSeG256etrNY27b\nwZYx9kJL2TOzeDOk2LPNtcFZMo86GzZswMzMTHYbWRh1rpuMbee6khnxCx4EFSY6eBBUmLZJdDZh\nnbnQmETPjWKqy+Suri66VA9L5MBcM6xNbDlbL5EDSzJo8eRkrrkxPz/vJn9kUtWW6WsrSy5Rtm+3\nW5klxp4Fy6/myfdWly5iuc+ZvPai1RYWFtqTky0Igo8u0cGDoMJ0LJItN+GDhuULt5SNSg8MDNBV\nGlk72Iislqe5E2DY96wk19F8NrWwvh79PTZRhE1KyR3NX40kDPY7uSmPc0fzve/Y42wKb2am5Ep0\nJq9zEpqERA+CYEWigwdBhYkOHgQVpm02OMuLbu0XL/qLHcfcX9q2a8YGz849TWxVTa6dmrv+mD1W\nl7FINn1dPT09hbbkjiGwiL2y5BJlxzG3IYswZBGMXhmLZGPLaFm89yUn2eZKZTbRCXO1rdTOOvEL\nHgQVJruDi0iXiPxCRJ6s7Q+LyH4RGa39v2ntmhkEQSs0I9HvBvAKgHNr+3sAPJVSuk9E9tT27/W+\nbCWFjlCzstCTePa43Nxcuj4mw3Nlj4XJyTJ3Rk9PT7Zss/Xb69RLFOnvWSmv69D3gyVksDD3lJXb\nXnuZhNb3iuWhZ26s3IQPHmzyim0XS1qSW8a2WUKJVZXoIrINwL8EcL/6+BYAe2vbewHcmnXGIAja\nhuT8lRORxwD8VwCDAP5DSulzIjKRUqovKywAxuv7ADA5OdmoeHR0dNUbHgQBMDIy0tgeGhpaNnK7\nokQXkc8BGEspPS8iv1t2TEopiQj9S/Hoo48W9s8999zG9sDAQKGsv7+/sa0lqF3Uvq+vr/Q4oFyG\nDg8P48SJE4XjVkOiaznMUjvXt/v7+5ct1eNNsAGKaZOnp6cLZXpVT2+FT6B88sPNN9+Mffv2uRF7\nllyJziLvyp7Zddddh2eeeabwfPVxzIRbTYk+MjKy7McoV6Kz+eDWxNDPc2pqqrFdf7Zf/OIX8dBD\nD5WW2fofeOABckV5NvinANwsIp8F0AfgXBH5SwDHRGRrSumoiGwFMMYqYUnpmD3qJUUA4Npsdt/b\ntvs5SQHKYA/TvuDAUgfXnRHg637pfea6Y244lvxRd6DcmXK54ahsVpt1gXp12j+abE06DzZGk7t2\nnS1nMwzPdu09trx02b7HijZ4SukbKaVtKaUdAG4H8A8ppT8CsA/A7tphuwE8kXXGIAjaxtn4we8D\ncKOIjAL4/dp+EAQfIpqKZEsp/SOAf6xtvwfghtzv6uVlgaKdbWWs3s+VQWy2UK5Ez51xxEwFFllV\nZ/PmzZicnHRdRixPPJtVlBt51yq5cpjlTPNmss3OzrrjECwffm6Zleie6bewsNBSwgf2brbqJstd\nlpoRkWxBUGGigwdBhWnbZBMr0bUEsy4uLUeYy4XJIG/CympIdJZ/i0loXcfExES25G0lKip32R6b\naEHTzJJGnolh76m3WuzMzIw7Sp/rkrN16jJ7jd4SSouLiy2tWNrMKLr3vtgljXLzujHiFzwIKkx0\n8CCoMNHBg6DCdMwG15FcNgTVs7tZsgNmA+Xa4LlL4lr7hy0Z5Nlip0+fzrZxc5fLzV0y2OLZlexe\nWbyZfWxcQ9u+p0+fzp6Rpt8DFn3IEll4EXUsuYQtb9UG96IWm3GTsWWwNfELHgQVJjp4EFSYjkl0\nvW9nPmnJnpu7zZZ5Lhcrl3Lle+5KjyxZg5VgnhzOlXStou/N/Py8a4o0Y8547kBmGtjr8iR12YSd\nOmxSjecqBfhkpNxIttxlr9jkIe+4ubk5d4kjYHl/8ohf8CCoMNHBg6DCRAcPggrTNhvc2tnahrD2\nhM7mwRIElmVLqePZUdZWYjZ47mwyBrNNVyOBn97X94Ctn6Z5//333TY2E9bbyjpubNwhN9kGexas\nDq+N9rlYvPEWlpAhd/lgu+2Vle17xC94EFSY6OBBUGE+lBJd7zOJzqKdPKyUYlFXuckUcpckshLX\ni4TKnZFm62zm3Ppc7UyOoWEzvHIlOrsfufXlXjPgPyeW2ILJa90vrCT3ysr2PeIXPAgqTHTwIKgw\nbZPodiQxV6Ln5hxnS+SwdjBp70k8trJm7iqW3d3drjS0bdJlud4C+zkbKWcRal6drI0MNnrNUjtr\n9H4ryxPZ79nnkJtTLjf1N0v4oGW47RMxih4EASU6eBBUmOjgQVBh2maDW5tB29bMBtczy1jCh1xb\n2rajlXWumL1vy7TNrM/V39/vjiEwm5PNhmN496e3t5fanBpW5uV4Z+2zOcz1/dCJOO3Yi95nedE1\nbGwk130J+Ha3fa9yXVy6Ppt0sS1LFwVB8NEl6xdcRN4EcArAIoCFlNLVIjIM4H8C2AHgTQC3pZTG\n16aZQRC0QjMS/fdSSu+q/T0Ankop3Scie2r793pftpJCy6JWJTqbbKJhEp25oDxYkgErJ3Wb9bms\nRNffy5V0tk4W8Za7aiiLQstdoVNj74f3nb6+PjfJQzPP3Uv0wSLeNIuLi9nuL+bGYpFsXi70jqwu\nSrgFwN7a9l4At55FXUEQrAGSEyggIm8AmMSSRP8fKaW/EJGJlNLGWrkAGK/vA8Dk5GSjYruoehAE\nq8PIyEhje2hoaNmIZq5Evz6ldEREtgDYLyKHdGFKKYkI/Uvx1a9+tbA/ODjY2B4YGCiUbdzY+DuB\noaGh0u8AxRVK9RxyoDgKW5d41113HZ599tnCcWxE1oNJdFZWl4VXXnklDhw44Ervdkn0z3/+83j8\n8cfd49iEklxyJpvcfvvteOSRR9ZUojPvRp1Pf/rT+PGPf7wqEv3MmTON7ampqULZ9PR0Y3tiYqKx\nferUKQDAt7/9bdx1110FU1WnGQc+MGNfeumlZdehyergKaUjtf/HROQHAD4J4JiIbE0pHRWRrQDG\nWB0sQQBLKKc7rn3QubO49ItpZ7UxG9x7WVi4qC3Tf0B0O84555yWZsNZvHzkzbi7vHOz2V657jmG\nrqO3t9ft1Oy5MJclm03m/SGzSSjZsr25Nrh957w/2M3Y4Lk59Vd8q0RkvYgM1rcB/AGAgwD2Adhd\nO2w3gCeyzhgEQdvI+QW/AMAPan8ZuwH8VUrpb0XkZwAeFZE7AbwF4La1a2YQBK2wYgdPKb0OYFfJ\n5+8BuCH3RFZieMvI2mO9mTdAvpssN684cwvpNjJZyNxkGrtkMnNPManp5eZmYwFW4ubOBGOyvJWZ\nYPqZ9fe0ACfnAAAJLUlEQVT3u8k9mPmSG33I3GQ2icPZ5lMDiu+qldNeXnSWG87WkfvMIpItCCpM\ndPAgqDDRwYOgwrRtNhlLXsdsFL3djD/Us7+YDW7r8Oxn5iZjNqf9nK255cFs69wsH8ytl5vEMXcc\nwtbvhefaUFWWrSfXtmbjGtqWtkkoczOp5NrgzNXmJW5cXFx07fOyfY/4BQ+CChMdPAgqzIdColsX\nmjdZXs8ss3WwEE4m0W0iRK+MJX9ksOgvzw2Xm4OdlVlp7EnBnp4eOutKw5JeePKaRfbZ6/eSb7AE\nj7lLF7Eln6xbLDeRQ647jUl0b/kmu7x07lJUlvgFD4IKEx08CCpM2yQ6g01Eyc2PxUbRNVZyWbnK\n2lWHSfnc0U4muVo1AXLltSb3vllYbngWhebJ9+7ubtfEaGYWmye9WSSljU5jM/u8CDU22YTJd9Ze\nT8oDMYoeBAGigwdBpYkOHgQV5kNpg3vuDTajxtpYnnuKudNy7Rxrq7diK9m1uDQ52UfKztdKPvKe\nnh436SKz83KTXrBkDfY4L4EHawdzf7WSJWclGzw34QOzz72xAbZGWkSyBUGwjOjgQVBh2ibRmczM\nnUzA5JiVxt4kjNykhfbcGpZDjrm4bLSaF42Um6zB1rnaEj0nYWLZfm7ONLvtyfJmoiA995R97p68\nnp2dzZb2XuIGez7bfu+dsxKd5RwMiR4EQXTwIKgy0cGDoMK0zQZnM6nYsbm2GEuYyGad6TIWtqnb\nZNdS88I0LbqNc3Nz7vesfZWbUKKZENc61j2V6ybLfZ7sOHsu71mzxQdykyKy3OTWbZW79G/u2mEs\nh7nnprU2uCV72eiso4Ig+EgSHTwIKkzbJLqVj8yVoqUrkyIs0sdz8VjJpc/FTABNq3m6dcKKubm5\nlnKyMXLvKVuuqVWJnisZvediI7dyXWFMQmupbSW6J/Pn5uaylyRibcxd3jc3uQfLnceIX/AgqDBZ\nHVxENorIYyJySEReEZFrRWRYRPaLyGjt/01r3dggCJojV6L/GYC/TSl9QUR6AQwA+GMAT6WU7hOR\nPQD2ALjXPRGR6BZPorPRZRbJpsmd2FLW5hxYVJ6u//Tp0wXJruV6bvSXJTfxRG6EVCsJJAB/pJzV\nb3OQ5U7WsNLbk+VsVVl7HIt89NrYzOqfnrxuZrXZ3KWuc1YXHQLwOwAeAICU0lxKaQLALQD21g7b\nC+DWrDMGQdA2ZKWYVhG5CsBfAPgllhYhfB7A3QCOpJQ21o4RAOP1fQCYnJxsVDw6Orr6LQ+CACMj\nI43toaGhZfIuR4N2A/htAHellH4qIn+GJTneIKWURIT+pfja175W2GdypL+/v7Hd19fX2LZpk3WZ\nXa2zTMJ8/etfx3e+853CZ0yueqtd2nbo9g4MDLhl9e9dc801+MlPfrLqEp0Fh5RJy6uvvhrPPfec\nOye5mVxommYl+vXXX4+nn3561SX6mTNn3OPKJPo999yDb33rW3RVEq+NLODGogOlpqamGtv163r4\n4Ydxxx130FwF9WMPHjzongfI6+CHARxOKf20tv8Yljr4MRHZmlI6KiJbAYzRE2XOsgL4sjWaXNcM\ngy3Vo1+sVnOVe+c6c+ZM4UXVf6Ds/WDuEs+mzV0qitm+lly3Gcvn7bVxbm7OtWmZq4rN7MvtqCwf\nOftDyZYWYnj55XNn4TVzvhVt8JTSOwDeFpHfqn10A5bk+j4Au2uf7QbwRNYZgyBoG7nDxHcBeKg2\ngv46gH+NpT8Oj4rInQDeAnDb2jQxCIJWyergKaUXAFxdUnRD9olWIZKtmaQRHkzqsIQPzSQ/8LBu\nMm13awnZjET3YEkGrP3JJu1oWrXJNd65ZmZmXJs2d0JJvZ6yslyJPj8/n2X72jpyx0YAvlyW/ryZ\nlXA9IpItCCpMdPAgqDDRwYOgwnQs6WKuS6BVl5lno7AEBCzv+mq4yawNrtuo3WS5OccZbHad3p6e\nns72fbdig9vveC406zZkSR1YssPcpIjezLXFxUWarMG7j7m56y02Aabezs37z4hf8CCoMCuGqraK\nDlUNgmDtKQtVjV/wIKgw0cGDoMKsmUQPgqDzxC94EFSYNe/gInKTiPxKRF6rZX5pCyLyXREZE5GD\n6rO2ppkSke0i8iMR+aWIvCwid3eoHX0i8qyIvFhrx590oh2qPV0i8gsRebJT7RCRN0XkgIi8ICLP\ndbAda5oObU07uIh0AfhzAH8I4HIAd4jI5Wt5TsWDAG4yn+3BUpqpEQBPwcxrXwMWANyTUrocwDUA\n/m3t+tvdjlkAn0kp7QJwFYCbROSaDrSjzt0AXlH7nWrH76WUrkop1edZdKId9XRoO7GUUOWVVW1H\nfUWJtfgH4FoAf6f2vwHgG2t5TnP+HQAOqv1fAdha294K4FftakvtnE8AuLGT7cBSPr2fA/gXnWgH\ngG21l/YzAJ7s1HMB8CaA88xnbW0HgCEAb6A2FrYW7VhriX4hgLfV/uHaZ53igpTS0dr2OwAuaNeJ\nRWQHgE8A+Gkn2lGTxS9gKTHH/rSUwKMT9+O/AfiPAHR4VyfakQD8vYg8LyL/pkPtuBjAcQDfq5ks\n94vI+tVsx2/sIFta+vPYFheCiGwA8NcA/n1K6WQn2pFSWkwpXYWlX9BPisgV7W6HiHwOwFhK6XnS\nznY9l+tr9+MPsWQ6/U4H2lFPh/bfU0qfADCNknRoZ9OOte7gRwBsV/vbap91imO19FLISTO1GohI\nD5Y690Mppcc71Y46aSkj7o+wND7R7nZ8CsDNIvImgEcAfEZE/rID7UBK6Ujt/zEAPwDwyQ60oywd\n2m+vZjvWuoP/DMCIiFxcywZzO5ZSPXWKtqaZqmWbfQDAKymlP+1gO84XkXoG3H4sjQMcanc7Ukrf\nSCltSyntwNK78A8ppT9qdztEZL2IDNa3AfwBgIPtbkdqRzq0NgxmfBbAqwD+H4D/tNbnU+d9GMBR\nAPNY+kt5J4DNWBrgGQXw9wCG17gN12NJXr0E4IXav892oB3/HMAvau04COA/1z5vaztMm34XHwyy\ntft+/DMAL9b+vVx/LztxP7Dk1Xiu9mz+F4BNq9mOiGQLggrzGzvIFgS/CUQHD4IKEx08CCpMdPAg\nqDDRwYOgwkQHD4IKEx08CCpMdPAgqDD/HwK/qMTK8KRLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e91217d68>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def mean_face(faces):\n",
    "    \"\"\"Compute the mean of the `faces`\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    faces: (N, 64 * 64) ndarray representing the faces dataset.\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    mean_face: (64, 64) ndarray which is the mean of the faces.\n",
    "    \"\"\"\n",
    "    mean_face = mean(faces)\n",
    "    return mean_face\n",
    "\n",
    "plt.imshow(mean_face(faces).reshape((64, 64)), cmap='gray');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To put things into perspective, we can benchmark the two different implementation with the `%time` function\n",
    "in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 128 ms, sys: 0 ns, total: 128 ms\n",
      "Wall time: 149 ms\n",
      "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n",
      "Wall time: 5.03 ms\n"
     ]
    }
   ],
   "source": [
    "# We have some huge data matrix, and we want to compute its mean\n",
    "X = np.random.randn(100000, 20)\n",
    "# Benchmarking time for computing mean\n",
    "%time mean_naive(X)\n",
    "%time mean(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.65 s, sys: 0 ns, total: 3.65 s\n",
      "Wall time: 4.6 s\n",
      "CPU times: user 32 ms, sys: 0 ns, total: 32 ms\n",
      "Wall time: 29.2 ms\n"
     ]
    }
   ],
   "source": [
    "# Benchmarking time for computing covariance\n",
    "%time cov_naive(X)\n",
    "%time cov(X)\n",
    "pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively, we can also see how running time increases as we increase the size of our dataset.\n",
    "In the following cell, we run `mean`, `mean_naive` and `cov`, `cov_naive` for many times on different sizes of\n",
    "the dataset and collect their running time. If you are less familiar with Python, you may want to spend\n",
    "some time understanding what the code does. __Understanding how your code scales with the size of your dataset (or dimensionality of the dataset) is crucial__ when you want to apply your algorithm to larger dataset. This is really important when we propose alternative methods a more efficient algorithms to solve the same problem. We will use these techniques again later in this course to analyze the running time of our code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def time(f, repeat=100):\n",
    "    \"\"\"A helper function to time the execution of a function.\n",
    "    \n",
    "    Arguments\n",
    "    ---------\n",
    "    f: a function which we want to time it.\n",
    "    repeat: the number of times we want to execute `f`\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    the mean and standard deviation of the execution.\n",
    "    \"\"\"\n",
    "    times = []\n",
    "    for _ in range(repeat):\n",
    "        start = timeit.default_timer()\n",
    "        f()\n",
    "        stop = timeit.default_timer()\n",
    "        times.append(stop-start)\n",
    "    return np.mean(times), np.std(times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fast_time = []\n",
    "slow_time = []\n",
    "\n",
    "for size in np.arange(100, 5000, step=100):\n",
    "    X = np.random.randn(size, 20)\n",
    "    f = lambda : mean(X)\n",
    "    mu, sigma = time(f)\n",
    "    fast_time.append((size, mu, sigma))\n",
    "    \n",
    "    f = lambda : mean_naive(X)\n",
    "    mu, sigma = time(f)\n",
    "    slow_time.append((size, mu, sigma))\n",
    "\n",
    "fast_time = np.array(fast_time)\n",
    "slow_time = np.array(slow_time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "list indices must be integers or slices, not tuple",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-33-506e0e05078a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      1\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;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merrorbar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfast_time\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfast_time\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[0mfast_time\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'fast mean'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlinewidth\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0merrorbar\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mslow_time\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mslow_time\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[0mslow_time\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlabel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'naive mean'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlinewidth\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0max\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_xlabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'size of dataset'\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[0mset_ylabel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'running time'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mTypeError\u001b[0m: list indices must be integers or slices, not tuple"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAD5CAYAAACOAorsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAADTZJREFUeJzt3H+o3Xd9x/Hna4kFxZIWu4okKWZy1QVmofbXRtlay9ak\n/wShsLZiWVEkzMr+bLc/dNA/NpGBiK3ZKKH4T8OYReOIdsLQCjUuCP2VlvTcpdCkLRRbiVD/KJe8\n98c9heNdkvO9zcl935z7fEDgfs/53Ju3b9L7vOferzdVhSRJHf6gewBJ0sZlhCRJbYyQJKmNEZIk\ntTFCkqQ2RkiS1GZqhJLsT/JGkufP8nySfCvJYpJnk1wz+zElSfNoyCuhR4Fd53h+N7Aw/vMl4Dvn\nP5YkaSOYGqGqehJ46xxH9gDfrWWHgcuSfGRWA0qS5tfmGXyMrcCJieuT48denzx06tQpfzWDJM25\nLVu2ZDXnvTFBktRmFhF6Fdg+cb1t/JgkSec0iwgdBO4Z3yV3I3Cqql6f9k46t9Fo1D3CRcNdDeeu\nhnNXa2Pqz4SSPAbcDFyR5CTwNeB9AFW1DzgE3A4sAr8D7r1Qw0qS5svUCFXVXVOeL+DLM5tIkrRh\neGOCJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqY4QkSW2MkCSpjRGSJLUxQpKkNkZIktTG\nCEmS2hghSVIbIyRJamOEJEltjJAkqY0RkiS1MUKSpDZGSJLUxghJktoYIUlSGyMkSWpjhCRJbYyQ\nJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqY4QkSW2MkCSpjRGSJLUxQpKkNoMilGRXkmNJ\nFpM8cIbntyT5YZJnkhxNcu/sR5UkzZupEUqyCXgI2A3sBO5KsnPFsS8DL1TV1cDNwL8kuWTGs0qS\n5syQV0LXA4tVdbyq3gEOAHtWnCng0iQBPgi8BSzNdFJJ0tzZPODMVuDExPVJ4IYVZ74NHAReAy4F\n/rqqTp/tA45Go1WOuTG5p+Hc1XDuajh3NczCwsJ7ft8hERriNuBp4DPAx4CfJPl5Vf32TIfPZ+CN\nYjQauaeB3NVw7mo4d7U2hnw77lVg+8T1tvFjk+4FHq9li8DLwCdnM6IkaV4NidARYCHJjvHNBney\n/K23Sa8AtwIk+TDwCeD4LAeVJM2fqd+Oq6qlJPcBTwCbgP1VdTTJ3vHz+4AHgUeTPAcEuL+qfn0B\n55YkzYFBPxOqqkPAoRWP7Zt4+zXgr2Y7miRp3vkbEyRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLa\nGCFJUhsjJElqY4QkSW2MkCSpjRGSJLUxQpKkNkZIktTGCEmS2hghSVIbIyRJamOEJEltjJAkqY0R\nkiS1MUKSpDZGSJLUxghJktoYIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJ\nUhsjJElqY4QkSW2MkCSpzaAIJdmV5FiSxSQPnOXMzUmeTnI0yc9mO6YkaR5tnnYgySbgIeAvgZPA\nkSQHq+qFiTOXAQ8Du6rqlSRXXqiBJUnzY8groeuBxao6XlXvAAeAPSvO3A08XlWvAFTVG7MdU5I0\nj4ZEaCtwYuL65PixSR8HLk/y0yS/SnLPrAaUJM2vqd+OW8XH+TRwK/B+4BdJDlfVS2c6PBqNZvTX\nzjf3NJy7Gs5dDeeuhllYWHjP7zskQq8C2yeut40fm3QSeLOq3gbeTvIkcDVwxgidz8AbxWg0ck8D\nuavh3NVw7mptDPl23BFgIcmOJJcAdwIHV5z5AXBTks1JPgDcALw421ElSfNm6iuhqlpKch/wBLAJ\n2F9VR5PsHT+/r6peTPJj4FngNPBIVT1/IQeXJF38Bv1MqKoOAYdWPLZvxfU3gG/MbjRJ0rzzNyZI\nktoYIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqY4QkSW2MkCSp\njRGSJLUxQpKkNkZIktTGCEmS2hghSVIbIyRJamOEJEltjJAkqY0RkiS1MUKSpDZGSJLUxghJktoY\nIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqMyhCSXYlOZZkMckD\n5zh3XZKlJHfMbkRJ0ryaGqEkm4CHgN3ATuCuJDvPcu7rwH/NekhJ0nwa8kroemCxqo5X1TvAAWDP\nGc59Bfge8MYM55MkzbHNA85sBU5MXJ8Ebpg8kGQr8FngFuC6aR9wNBqtYsSNyz0N566Gc1fDuath\nFhYW3vP7DonQEN8E7q+q00mmHj6fgTeK0WjkngZyV8O5q+Hc1doYEqFXge0T19vGj026FjgwDtAV\nwO1Jlqrq+zOZUpI0l4ZE6AiwkGQHy/G5E7h78kBV7Xj37SSPAv9pgCRJ00yNUFUtJbkPeALYBOyv\nqqNJ9o6f33eBZ5QkzalBPxOqqkPAoRWPnTE+VfU35z+WJGkj8DcmSJLaGCFJUhsjJElqY4QkSW2M\nkCSpjRGSJLUxQpKkNkZIktTGCEmS2hghSVIbIyRJamOEJEltjJAkqY0RkiS1MUKSpDZGSJLUxghJ\nktoYIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqY4QkSW2MkCSp\njRGSJLUxQpKkNkZIktTGCEmS2hghSVKbQRFKsivJsSSLSR44w/OfS/JskueSPJXk6tmPKkmaN1Mj\nlGQT8BCwG9gJ3JVk54pjLwN/UVV/AjwI/NusB5UkzZ8hr4SuBxar6nhVvQMcAPZMHqiqp6rqN+PL\nw8C22Y4pSZpHQyK0FTgxcX1y/NjZfAH40fkMJUnaGDbP8oMluYXlCN10rnOj0WiWf+3cck/Duavh\n3NVw7mqYhYWF9/y+QyL0KrB94nrb+LHfk+RTwCPA7qp681wf8HwG3ihGo5F7GshdDeeuhnNXa2PI\nt+OOAAtJdiS5BLgTODh5IMlVwOPA56vqpdmPKUmaR1NfCVXVUpL7gCeATcD+qjqaZO/4+X3AV4EP\nAQ8nAViqqmsv3NiSpHkw6GdCVXUIOLTisX0Tb38R+OJsR5MkzTt/Y4IkqY0RkiS1MUKSpDZGSJLU\nxghJktoYIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2RkiS1MYISZLaGCFJUhsjJElqY4QkSW2M\nkCSpjRGSJLUxQpKkNkZIktTGCEmS2hghSVIbIyRJamOEJEltjJAkqY0RkiS1MUKSpDZGSJLUxghJ\nktoYIUlSGyMkSWpjhCRJbYyQJKmNEZIktTFCkqQ2gyKUZFeSY0kWkzxwhueT5Fvj559Ncs3sR5Uk\nzZupEUqyCXgI2A3sBO5KsnPFsd3AwvjPl4DvzHhOSdIcSlWd+0Dyp8A/VtVt4+u/B6iqf5o486/A\nT6vqsfH1MeDmqnr93TOnTp06918kSbrobdmyJas5P+TbcVuBExPXJ8ePrfaMJEm/xxsTJEltNg84\n8yqwfeJ62/ixVZ1Z7Us0SdL8G/JK6AiwkGRHkkuAO4GDK84cBO4Z3yV3I3Bq8udBkiSdydQIVdUS\ncB/wBPAi8O9VdTTJ3iR7x8cOAceBReAx4Epv555uwK3vnxvv6LkkTyW5umPO9WLavibOXZdkKckd\naznfejJkV0luTvJ0kqNJfrbWM64XA/473JLkh0meGe/q3o4514Mk+5O8keT5szy/+s/vVTWzP8Am\n4H+BPwIuAZ4Bdq44czvwIyDAjcAvZznDxfJn4K7+DLh8/PbujbqrofuaOPffLH9hdEf33Ot1V8Bl\nwAvAVePrK7vnXse7+gfg6+O3/xB4C7ike/amff05cA3w/FmeX/Xn91nfmHA9sFhVx6vqHeAAsGfF\nmT3Ad2vZYeCyJB+Z8RwXg6m7qqqnquo348vDLP+sbaMa8m8L4CvA94A31nK4dWbIru4GHq+qVwCq\naqPua8iuCrg0SYAPshyhpbUdc32oqidZ/t9/Nqv+/D7rCHk793Cr3cMXWP4KY6Oauq8kW4HP4v9Z\nesi/rY8Dlyf5aZJfJblnzaZbX4bs6tvAHwOvAc8Bf1dVp9dmvIvOqj+/D7k7Ts2S3MJyhG7qnmWd\n+yZwf1WdXv6iVeewGfg0cCvwfuAXSQ5X1Uu9Y61LtwFPA58BPgb8JMnPq+q3vWPNh1lHaCa3c28Q\ng/aQ5FPAI8DuqnpzjWZbj4bs61rgwDhAVwC3J1mqqu+vzYjrxpBdnQTerKq3gbeTPAlcDWy0CA3Z\n1b3AP9fyDz0Wk7wMfBL4n7UZ8aKy6s/vs/52nLdzDzd1V0muAh4HPu9XqNP3VVU7quqjVfVR4D+A\nv92AAYJh/x3+ALgpyeYkHwBuYPnu141myK5eYfkVI0k+DHyC5buB9f+t+vP7TF8JVdVSkndv594E\n7K/x7dzj5/exfNfS7Szfzv07lr/K2HAG7uqrwIeAh8df3S9V1bVdM3cauC8xbFdV9WKSHwPPAqeB\nR6rqjLfdzrOB/64eBB5N8hzLd33dX1W/bhu6UZLHgJuBK5KcBL4GvA/e++f3qb/AVJKkC8XfHSdJ\namOEJEltjJAkqY0RkiS1MUKSpDZGSJLUxghJktoYIUlSm/8DbKzkleaRRjMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb36521a2e8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(fast_time[:,0], fast_time[:,1], fast_time[:,2], label='fast mean', linewidth=2)\n",
    "ax.errorbar(slow_time[:,0], slow_time[:,1], slow_time[:,2], label='naive mean', linewidth=2)\n",
    "ax.set_xlabel('size of dataset')\n",
    "ax.set_ylabel('running time')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "## === FILL IN THIS, follow the approach we have above ===\n",
    "fast_time_cov = []\n",
    "slow_time_cov = []\n",
    "for size in np.arange(100, 5000, step=100):\n",
    "    X = np.random.randn(size, 20)\n",
    "    f = lambda : cov(X)               # EDIT THIS\n",
    "    mu, sigma = time(f) # EDIT THIS\n",
    "    fast_time_cov.append((size, mu, sigma))\n",
    "    \n",
    "    f = lambda : cov_naive(x)         # EDIT THIS\n",
    "    mu, sigma = time(f) # EDIT THIS\n",
    "    slow_time_cov.append((size, mu, sigma))\n",
    "\n",
    "fast_time_cov = np.array(fast_time_cov)\n",
    "slow_time_cov = np.array(slow_time_cov)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEGCAYAAAAT/1CLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4VGXawOHfMzPphVCko4AGFQHLIlW6Kwoqwtqwsbpr\nd62r4uraXcvu6mdn14KoYFtFbKyKgl1RV1R6aCKdAOllMjPv98eZTCZlkpNkZjKZPPd1zZXJKXPe\nczI5z3m7GGNQSiml4pWjpROglFJKRZIGOqWUUnFNA51SSqm4poFOKaVUXNNAp5RSKq65WjoBkZSf\nn69NSpVSKs61a9dO6luvOTqllFJxTQOdUkqpuKaBzi8nJ6elkxDT9PqEptcmNL02oem1CS3c10YD\nnVJKqbimgU4ppVRci+tWl0qptsEYQ1FRET6fr6WTYltycjL5+fktnYyYVPPaOBwO0tPTEam3cWVI\nGuiUUq1eUVERSUlJJCYmtnRSbEtKSiI5ObmlkxGTal4bt9tNUVERGRkZTfo8LbpUSrV6Pp+vVQU5\n1TiJiYnNyq1roFNKtUlZs7eSNXtrSydDRYEGOqWUCoNZs2YxZMgQLrzwwkbtl5eXx9NPPx2hVDXO\ne++9x0MPPdTSyQg7raNTSqkweOaZZ3jzzTfp0aNHo/bLz8/nmWee4Y9//GOEUmaPx+Nh0qRJTJo0\nyfY+P+S6ATiyU2wXG2uOTimlmumaa65h06ZNnHbaaTz++ON8//33/Pa3v2XUqFEcd9xxgQ7Qq1at\nYvz48RxzzDGMGzeO9evXc8cdd7Bx40aOOeYY/vrXv9b67JdeeokRI0YwcuRILrroIgB++eUXTjrp\nJEaMGMHJJ5/Mr7/+Sn5+PgMGDAjUZRUXF3PYYYdRUVHBnDlzGDduHCNHjuTcc8+lpKQEgEsvvZRr\nrrmGCRMmcOuttzJ37lyuv/56ABYuXMiECRMYNWoUU6ZMYdeuXQDce++9XH755UyePJmzxg3m9TlP\n1ZvW3Nxczj33XMaNG8e4ceP4+uuvI/RXCE1zdEqpJqus48o7v3G5mEhrTN2bnW0bOr+HHnqIRYsW\n8fbbb9OxY0cKCgpYuHAhLpeLJUuWcOedd/LCCy/w7LPPcskll3D66adTUFBAQkICt912G6tWreLz\nzz+v9bmrVq3iH//4Bx988AEdO3Zk3759ANxwww1Mnz6ds846ixdeeIEbb7yRefPmMXDgQD7//HNG\njx7N+++/z/jx40lISOCkk05ixowZANx999288MILXHzxxQBs27aNDz74AKfTydy5cwPHHj58OIsW\nLUJEeP7553n44Ye55557AGvkkrfffpuvftnLeceN4LYrL2LdunV1pnXmzJlcdtllDB8+nF9//ZXf\n/e53LF261MZfJnw00CmlVJgVFBRw6aWXsmHDBkSEiooKAIYMGcI///lPtm3bxnHHHUf//v3r/ZxP\nP/2UU045hY4dOwLQvn17AL799ltefPFFAM4880xuu+02AKZNm8b8+fMZPXo0r7/+eqA4dOXKldxz\nzz3k5+dTVFTEhAkTAseYMmUKTqez1rG3bt3K+eefz86dO3G73RxwwAGBdccddxxJSUm069CRrI6d\n2LVrV8i0LlmyhNWrVwf2LSwspKioiPT09EZc0ebRQKeaJVaf6FXbZuf7GMnv7j333MOoUaOYO3cu\nv/zyCyeeeCIAp512GoMHD+b999/n7LPP5v/+7//o3bt3WI75Q66bHkMmsOjOO9m3bx8//vgjo0eP\nBuCyyy5j7ty5DBw4kLlz51bLPaalpdX5eTfccAOXX345kyZN4rPPPuO+++4LrEtKSgq8dziceDye\nkOny+XwsWrSoRfsMah2dUkqFWUFBAd26dQNg3rx5geWbNm2id+/eXHLJJUycOJEVK1aQkZFBYWFh\nnZ8zevRo3nzzTfbu3QsQKA4cMmQIr7/+OgCvvvoqw4cPByAlLZ2jjjqKG2+8kYkTJwZyakVFRXTt\n2pWKigpee+012+fQvXt3wKp7a0iotI4fP55///vfge1++uknW8cPJw10SikVZldddRV33nkno0aN\nwuv1BpbPnz+f4cOHc8wxx7BmzRqmT59Ohw4dGDZsGMOHD6/VGOXQQw/luuuuY/LkyYwcOZK//OUv\nADzwwAPMnTuXESNG8Morr1TLbU2bNo1XX32VqVOnBpbdfPPNTJgwgYkTJ5KdnW3rHGbOnMmMGTMY\nM2ZMoDiyPqHSev/99/PDDz8wYsQIhg4dyuzZs20dP5zEmPidhLsxM4zn5OTY/gK0RaGujxZdtu3v\nTkN//2hdm/z8fNq1a9eofVr6u1tWVhbW4ryWaOofqWPWdW3q+xs3NMO41tEppdqk4ADXWvqDqabR\nokullFJxTQOdUkqpuKaBTimlVFzTQKeUUiquaaBTSrVJ6TPGkj5jbEsnQ0WBBjqllGoBc+bMsdUR\nO9acdtpp5OXltXQyGkW7F6i40dL9opRqjBkzZjS6H11LdoMwxmCMsT2ySizRHJ1SSjXTL7/8wpAh\nQ7jyyisZNmwYU6dOpbS0FCDkFDl///vfefTRR1m7di3jx4+v9lkjRowAYNmyZUyaNIkxY8Ywbdo0\n9uzaWevYu3bt4uyzz+YPJ47lDyeO5ZtvvgHgscceY/jw4QwfPpwnnngCgNtvv52nnqqaVufee+/l\n0UcfpaioiJNPPpnRo0czYsQI3n333UBaBg8ezMUXX8zw4cPZsmULAwcOZM+ePQDccsl5XDTlWIYN\nG8Zzzz0X+NwePXpw1113MXLkSI499tjAFD+VaR05ciQjR44MpPWVV14JTF909dVXVxtNJhw0R6dU\nK9Kacq1WWlPJa4FBYxpT95Y+YyyjGtimaM6SBj9n/fr1PP300zzyyCP8/ve/56233uKMM86od4oc\ngH79+uF2uwPjYM6fP5+pU6dSUVHBDTfcwLx58+jUqRNvvPEGTz/4N2687+Fqx73xxhsZOXIkf354\nNl6vl+xkN8uWLWPevHksWrQIYwzHHnssI0eOZOrUqdx0002BWdDffPNNXn/9dZKTk3nxxRfJzMxk\nz549HHvssYEJWNevX8+TTz7J0UcfXeucb7jvYTKz2nNImpfx48dz8skn06FDB4qLixk8eDB//etf\nufXWW5kzZw7XX399IK1z587F6/VSVFTEmjVreOONN3j//fdJSEjguuuu4/XXX+e8885r8JrbpYFO\nKaXC4IADDmDQoEEAHHHEEWzevBmof4qcSlOnTmX+/Plcc801vPHGG8yePZucnBxWrVrFKaecAliz\nAKS236/Wvp9++imzZs1iZSE4nU7atWvHV199xeTJkwMzE5x44ol89dVXXHLJJeTm5rJ9+3Zyc3PJ\nysqiZ8+eVFRUcNddd/HFF1/gcDjYvn17IBfWq1evOoMcwBtznuKzD98jxSls3bqV9evX06FDBxIT\nEzn++OMD12Lx4sXV0gpVaX3llVf48ccfGTduHGAN/5WVldW0P0IIGuiUUnHHTg6sMtdXNGdJvXVf\nP+S6IdfdYL1Y8NQ1TqczUHRZ3xQ5laZNm8aMGTM46aSTEBEOPPBAVqxYwSGHHMKHH35YPS3NNGXK\nFBYsWMCuXbsCAz+/+uqr5Obm8sknn5CQkMDAgQMpKysDQk/j89lnn/H9l5/y+GvvMbxXFpMnTw7s\nk5CQgIgErkV90/gYY5g+fXpgTj0g8DnhonV0SikVQXamyOnTpw9Op5MHHniAadOmAZCdnU1ubm5g\nNu6Kigo2rl1da98xY8bwzDPPAOD1esnPz2f48OG8++67lJSUUFxczLvvvhuYymfatGm88cYbLFiw\nIJBbLCgooFOnTiQkJPDpp5/y66+/NnheBQUFZGRmkZySytq1a/nuu+8a3KeutI4ZM4YFCxawe/du\nwJrex87xG0MDnYo5WbO3BuqilGrt7E6RUzm9TmXwSUxMZM6cOdx2222MHDmSUaNGseJ/39ba7777\n7uOzzz7jgkljuPiUY1mzZg1HHHEEZ511FhMmTODYY4/l3HPP5fDDDwes6XSKioro1q0bXbt2BeD0\n009n2bJljBgxgpdffpl+/fo1eF7HHnssXq+HGRNHcvvttzN48OAG96lM64gRIxgzZgxr1qzhkEMO\n4ZZbbmHq1KmMGDGCU045JVBsGi46TY9fW55qxY6mTNPT1IYT0d6vuaL53Ym1xiiR+Ps3RVOm6WlU\n0WWIdc3RlGl6WiKdTU1Pc+g0PUopFQZ26vFUfNCiS6VijBbdKhVeGuiUUiqG/JDrDkvrSlVFA51S\nSqm4poFOtQlaHBjfHA4HbrfmguKV2+3G4Wh6uNLGKEqpVi89PZ2ioqJAJ+3G+nazNf5k34TURq1r\njoKCAjIzM8OWlkilsz7RujYOh4P09PQmf54GOtWqxFrzehUbRISMjIwm73/Dj0UAXHRU7ebr9a1r\njl27dtGrV6+wpSVS6axPtK9NU2nRpVJKqbgWtUAnIseLyBoRWSciM+tYLyLyiH/9TyJyVEP7ishd\n/m2XicgHItI9WuejVHNonWHbljV7K0d/Hr0ixrYuKoFORJzA48AJQH9guoj0r7HZCUC2/3UR8KSN\nff9ujBlkjDkCeAe4NdLnopSKTfrwoEKJVo5uCLDOGLPBGOMGXgam1NhmCvC8sXwNZIlIt/r2NcYU\nBO2fBsTveGZKKdWKxNKDR7Qao/QAgoej3gIMtbFNj4b2FZF7gPOAfGBcqATk5OQ0mEg724RLZbHF\nt8eURO2YzVX39Ult4rr6NPUzI5EWe2LjHFtCLKW1OceL9nm0nu94/feqyKWnof0aM75sq291aYy5\nGbhZRG4CrgBuq2u7hi5K1Ad1/tx60mktA0mHvD71nUdTz7GpnxmJtNjQpO9ODJ5HKPW2dI2ltDbn\neNE+j9b0HW+B9IT7fhytosutQHBb0Z7+ZXa2sbMvwFzgd81OqVJKqbgSrUD3LZAtIn1EJBE4E3ir\nxjZvAef5W18OA/KNMdvr21dEgkP+FKD2rIRKKaXatKgUXRpjPCJyBfA+4ASeNcasEJFL/OtnAe8B\nk4B1QAlwfn37+j/6PhE5GPABvwCXRON8lFLNox3/VTRFrY7OGPMeVjALXjYr6L0BLre7r3+5FlUq\npZSql46MopRqE2KpubuKrlbf6lKplmTdOFPJax2NZ5VqkzRHp5RSKq5poFMKLdZSKp5poFMqQjR4\nKhUbNNAppZSKaxrolFJKxTUNdEo1QIsglWrdtHtBKxOJESUa+kxtQq+Uas00R6eUUiquaaBTSrUa\nWoysmkIDnVJKqbimgU4ppTklFVXR/r5poFNKKRXXNNAppZSKaxrolFKqDWjLxdMa6JRSSsU1DXRK\nKaXimgY6pZRScU0DnVJKqbimgU4ppVRc00DXRrTlFldKqbbNVqATkSQRuUdENohIvn/ZcSJyRWST\np5RSrZM+XMYOuzm6h4ABwNmA8S9bAVwaiUQppZRS4WJ3PrqpwEHGmGIR8QEYY7aKSPgmRVMBTZ1z\nLhJz1SmlVGtnN0fnpkZQFJH9gD1hT5FSSrUS8VI8GS/nEYrdQPcaMEdE+gCISDfgMeDlSCVMKdU4\n8XKzipfzULHDbqD7C7AR+BnIAnKAbcCdEUqXUkqpGGaMaXijGGEr0Blj3MaYa4wx6UAXIMP/e3lk\nk6eUUioW/bS3IvB+Q4GnBVPSMNv96EQkVUQGAdnAcBEZISIjIpc0pZRSseqldSWB9y/mFLdgShpm\nq9WliJyHVSfnBkqDVhlg/wikSymllE155b7A+02FHnpn2G1Q3zRur+G19VWhYF5OCX85MhOXQyJ6\n3Kaym6N7APidMaaTMaZX0EuDnFKtgDbwiD8+Y/hkWxl//GQvB7+yPbD8kZ+LIn7sD7eUsScouO4o\n9fHhlrKIH7ep7IZ9N7AkgulQbYz2+Wv9thZ7A+9fWFvM9INSY/aJPt7cv6yAuTklbC6y/gbBV/3T\n7ZFvOjEvqNiy0gs5JZywf0rEj90UdnN0fwUeFJFOkUyMUqr1eHJFVc7hT1/kMfLNXbzzS2mrao3X\nWt37QyGbi7z0Sncy84gMfjytS2DdugJPtYeQcMst8/L+r2UEP9O4BN7/tYwdJZE7bnPYDXRrgZOB\nnSLi9b98IhKbZ6WUiqh8t485a6saIByQ7mRNvodzPt7LxHdz+WKHNsgOt3x3VVHhKb1TmH9cR348\ntQszj8xk//TqhXORzNX9Z0MpHgPH9kgKLDu+VzJeAy/XkdOLBXYD3QvA88DhQD//K9v/UynVxjy3\nppjCiqqc27fTuvDA0Hbsl+xg6W43kxfmtmDqwsNnDPuC6qFa2qyVVTno58Z1YFyPZBxSd1FxJANd\nZWvLsw5KCyw7t5/1/oWc4pjM0duto+sI3Gpi8QyUUlHl9ppqN12ARKdwUf90pmen8sSKIh79uYgi\nj3W7yCv3kZUUuzOClXiqgtmJC3ezt8xHbrmPPWU+vDFyx8sr9/H4CvuNTD7bXo4xBgkRCJtq+d4K\nftxTQbtE4fheyYHlE3ok0T3VwfoCL1/udDOya1I9nxJ9dr99s4FzI5kQpVTr8J8NJWwv8XFoVu3n\n5IwEBzcekcmyoDqjmkEx1jyzuqoI9vMdblbmedhVagW5zMSqQBEcEKPtiZVFFLjtRd2OSQ62FHvZ\nUBD+mqXK3NypfVNJdlVdG5dDAjm859fGXp86u4FuCPC0iKwRkU+DX3YPJCLH+/dfJyIz61gvIvKI\nf/1PInJUQ/uKyN9FZLV/+/kikmU3PUqpxjPG8NhyK3BdMSA95Hadkp2B90+sLKrWzyuW+IxhdlCg\nWzCxI59P6cyqM7qy67zubD67e2DdVzvdLZFE9pX7mNWI3NzoblZuKtzFlxU+w2sbrEA3/aDUWuvP\n6Wcte2tTWcz9ve0GuqeAC4G/Ac/UeDVIRJzA48AJQH9guoj0r7HZCVj1ftnARcCTNvb9EBhgjBmE\n1WDmJpvno5Rqgo+2lrMyz0O3VAen9a19s6tLgds0qtgtmj7dXs6Gwqqcz5juyQzokEC3VCeJzurF\nfou3tkwDm8eXF1FQYRjX3V5xYGWg+yTMge6jrWXsKvXRr52L33RKqLW+d4aL0d2SKPUaXt8YW41S\nbNXRGWPmNPM4Q4B1xpgNACLyMjAFWBm0zRTgeX894NcikuWfJaF3qH2NMR8E7f81cGoz06nigM+Y\nkJX0qnke8efmLumfXisQ1GfWyiIuOyyd9jFWV/fsavvFbB9vKwPaRS4xddhb5g0U/d50ZAaLtzUc\nvCoD3Wfby8P6v1BZbDn9oNSQdX/n9Uvl0+3lvLA2tgJdyG+diJwb9P6CUC+bx+kB/Br0+xb/Mjvb\n2NkX4AJgoc30qDjj9VXVX8z8Jj8mW341R5nH8PXOqptc8PlGy7JcN59uLycjQfj9wWkN7+A3vnsS\nhRWGx5fHVq5ue4mXdzeXYTder9zniXo/sUeXW416ju2RxJDO9nJ0fTOd9Exzsqfcx8p94RtseeFm\nq+/cGQeGzsmfuH8KWYnCsj0VIbdpCfXl6KZjdSuA0A1RDPBsWFPUBCJyM+AB5obaJicnp8HPsbNN\n+KTWc8xYWhdr6al73evbXUAiAP9eVYyvKI+LD6hocL+WSKvddU9+tYGfCxz8VOhgdZEDj6m6Iz/9\nzUbGd6p5043s3/hvX20HXJy0n5tdv6xnl839zumUz8fbknliRQETU3aSldCctIbvHJ/Z7MJrEhnX\n0cPiPa4GPxfglR9+YVJnb611TU3rVyty6JRY97qlK3OYtSIFEM7plE9Ozj5bn7lu3ToOT0tkS7GL\n13/aQlIPT7X1TU2r2wfDsrwUb9tATj37TeyYwCvbq4o2m/a3avh+nJ2dXe/6YCEDnTFmUtD7cbY/\nsW5bgV5Bv/f0L7OzTUJ9+4rI74ETgQn1dX9o6KLk5OQ06sI12+fWKdR5zFhaF+Vj+oyBz7cBUNz+\nAI6oeReoY7/cMi+zlu7Eeu4Ch8DTvybQt1snLjssPebOMdS6Cp9h/sZSYB8AN62ueoIXoH+Wi5V5\n1k3r+4osLs7uEJ602Nz3oz0uXAJ/GdmTnpUdlG3sN+03BzJvTy6LtpbzXklnbhvcrsW/x16f4Z0f\ndgJerhrchcXv72nwcwFWeLO4Kvi6NyE9+W4ffG6NTXn2j+n8c3g7pvZJrbXfO8WdKfUVMbFnEqcc\n1aNRxztRSnh31z5WejLJzu7Y5LQGrwP44+GdyO5bO63B+13RsYJXFuwK/N6Uv1W478e2CsxF5IcQ\ny7+zeZxvgWwR6SMiicCZwFs1tnkLOM/f+nIYkG+M2V7fviJyPHADcLIxJrYKhVWTLA+a4+qCJXsp\nrGi49dZt3xWQF9T0+pGRVuPbvyzNZ16MTx8CUOD28djyQo78z04u+nRfYPn47knMPCKDN47ryKaz\nu/Hl1Kom+x/8WkZFlIsvvQam9U2pCnKNcNORmYCV284ta/kBlT7YUsaWYi99M5yBOi07lmwrb3ax\nuPUwY9lb7uP8Jfs4f/Fe9tS4Lk/76w8rr11jVJ7TlzvK8YTpe5KZIEy2MZblwA4JHFlHY5WWZLdm\n+KCaC8SqjexrZ2djjAe4AngfWAW8aoxZISKXiMgl/s3eAzYA67BaeV5W377+fR4DMoAPRWSZiMyy\neT4qRi0JqmzfUOjl2i/z6r2xfL2znLk5JSQGfZPPyU7jniFWo4E/fZHHO7+Uhti75d32bT4DXtvB\nLd8WsKXYS792VUHkjYmdmHlkJuN7JNMusfq/ap7btEhz9z8NyGjSfr/ZL5GJPZMo9piojK7fkNlr\nrCBy/sFpthtrdEt1sKvUx4pm1nsFz9324PAs0lzC/E2lDJu/i3eDvqslHsMJvZJrl2rY0CPNyUGZ\nLgorDD/khqe+bGqfFFJc9q7Vudn263Cjod5AJyLPi8jzQGLl+6BlnwAr6ts/mDHmPWNMP2PMgcaY\ne/zLZhljZvnfG2PM5f71A40x39W3r3/5Qf7pgo7wvy6pfWTVmgS3KktzCa9tKGVuiPHzPD7DdV/l\nAXDlwOo34MsPS+fPh2fgNVbOMJYEB+6Hl1sdgUd0SeSlCR34empn25/z3uboBvDx3ZMY2KHpT+oz\n/TmTpxvR0jESfin08OGWcpKccFa2vS4SAGO7WyOBLN7W9OloVu2r4LvdQaUWh6TxxSmdGdElkd1l\nPs7+uPp3deaRTXuwABjTvfn96YJzg2fV0XculN/1rcr5fR4D4542lKNb738Fv1+Pleuai9XMX6mw\nKPMYvgpqWfj3YVau7Iav81mTV/up9N+rilmxz8MB6U6uG1T7hnDzkRlceEga7tjqu8pTq6pu9Kf0\nTuGjE/fjvUn7ccL+KY1qCv7e5rKIty4tDxoD68qBoTuI23Fkp0RO6JVMiadlW8Q+v7YYA0zpnULH\noI7tDansx9ac/nRzc2o/tPXOcPHOCZ24d0g7gpNz4v7JHN6x8bm5SuHoT/fahqqHqSGd7acluATi\nis/3UWyjCiKS6g10xpg7jDF3YNWB3RH0utMY8y9jTGw9KqtW7etd5QRXU5yVncYZB6ZQ4jGcv2Qv\npUE3yO0lXu79oQCA+4e1q7NIRUS4f1g7Tgt6ugyuA2yu4LomuwFnWa6bW77ND/z+3LgO/Ga/xt/M\n9kt2sLnI2+xitIY8sKwg8H5MI+qyQmlODsWuhuoun/f38bqgEV0kAMb6A92XO8spa0KwrvAZXl5f\nd+mEQ4RLD0vn8ylVOfqbj2p83VywY7pa36tvdjU9vff9UPX3b+q4mZsKvdz+fUHDG0aQrTo6Y8z7\nkU6Iih8r9lZw+3f5DW9YQ11Pyv8cnsVBmS5W7vNw89Kqz7xlaT6FFVYdxvG9QleQO0R4YlT7wO8n\nLtzNt7vCU7d1d9A/773LChvcvsDt44Ile8OSw6wcUDeSxZdf7CjnwZ+q6tPCMUDw4R0Tmbx/csMb\nNlLwjXzSe7v5pTD0A8DuMh/9s1wMbUQOBaBzipOBHRIo81oPZY31/q9l5JbVPUZopYPaVRUNH9q+\neQ06OiZb6S33wtLdjf/Ov7i2hF+Kmt9wyCVWKcZnUZgQNpTYGqZAtVq/Fnl46KdCRry5k5ELdvF/\nQQ0OarYmC6WuUR/SExw8O7Y9iQ54dk1Vkd/rG0tJcQr3DW14pIqEoBki89yGU97PZUkz6lnAyhk+\nH1QM9cCywnoHLzbGcM2XeWwo9DKgGfVclSbtXxnomnceoeSV+7j4031EopBxZlArwu1h6IDt9Rn+\n+ElV4dK3uysY9dYuFmwK/RBw/iFpTQrczSm+fNH/fWlMvWBzBca9tDGiSrAyj+HvP4YnF/bnw61c\n/BWf76OohYowNdCpJgsurhv42k7u+L6Alfs8tE8Szj+46p/5zXpuOJVyy7z8tLeCpDqqTAZ1TOTu\no2sHtD8fnsEBGY1r6n7GgSkUewynf7iHt5vYGtMYw81L86lZQjbzm3xeDVE09UJOCa9vLCXNJTw3\ntn2d2zTG2O7JpLqsESjCPZu0MVYjny3F3jrHNGyu4AYtD/7YcE64PsYYrv0qj3eCAv6k/ZMpcBtm\nLN7LtV/mVSvyBkh1CafXM7pHfcb7Jxv9uJGBY2eJlw+3lOFqYGSRcGvqAM/PrClmW4mvWY2PKl13\neAYDOyTwS5GX279rmSJMDXSqyV4Kag2Z7ISpvVOYN6EDa87oxkMjqm7mr61vOKB84r9xDO9Sdz3Q\nhYemcWJQkddBma56R88P5clR7bnwUKuByozFe5vUz27hr2V8sr2crKApXO4cbOVSLvtsHx/8Wj2X\ntXJfBTd8bbUOfWhEVrXiqaZKcUkgd7EwzMWXr24oDQTlp8Z0aHiHZnhubTGbi5pez3jP/wqZs7ak\nWiOOueM7cP/QdoFSgAlv72J1UGOm0/qm1OquYdewzkkkO+HnvRXsKrX/gPHK+hK8Bib2SqZziv0G\nMM01omsiToHvc+0XXRZV+HjoJ+sB5JZm1hOCVaLyxKj2uMRqcftJIx8SwsFuh/FQY12eLSJjRCS2\nZtlTEVfmMfzth6qn8bVndmP2uA5M2j+l1mC/X+9y11tnAlXFlqFGaBcRHjumKnj+Y3g7khoxqHAl\nhwgPDG1iJt/jAAAgAElEQVTH9Ydn4DNw2ed5jZovze013OKvKwwugrtyYAZXDUjHY6wAGjwu5fmL\n91LmhXOyU5uck6hLJIovNxV6+LO/y8Z9Q9vRN7PxncMbo8IHf7dRv1mXWSuL+MdPhTgFZo+tCsgi\nwsX90/nwxP2s+t08D+Pe2h1Yf34jG6EES3YJI/wPY3Zv2MaYQLHlOVEstgRrfsDfdEps1ASy/1pZ\nTG6Zj6P3S+C4nuG5tQ/skMD1R1hFmH/6Yl8DW4ef3cea87Cmzbkd+KP/55NYnbrnARtEZHAE0qdi\n1NOri9gSVGSW2cAT8n82hM51GGMCHcXH1jMVSfAs1ZV9mppCRLj5qEzuPtoKVDO/sd9w5qnVxWwo\n9JLdzsUfDql+w7x9cCbnZqdS6jWcsWhPYPmafA+HZLm430Z9YmMc3ysZh8BnO8qtYaXC4OJP91FY\nYTj5gOSo3JSdAvPWlbChoHG5utfWlwT+bo+MzOKEOkbsOLxjIktO3o/pB1l/k0pN6YAdrPJhzG7x\n5be73azN99AlxcFve4a/IU5DGjPyS165j0eWV+bm2oV1hvJrB2UwqEMCm8PQwKWx7Aa6FcD1xpj9\njTEjjDH7A9cBP2CNPfkk8GiE0qhiTL7bxz9/atxT+GsbSkI2wV9X4GFLsZdOyY6w1AnYdcWADB4d\nmUVQW5Vq/cZq2lPm5X5/c/u7j25XrZELWAH0oRFZTN4/mfygIclSnMLssR1ISwhvTUHHZCdDOydS\n4YOPtoQnV/fNLjfdUx08PLJ9WG9yoUw/KBWvgfuWNa7u5tLPrFzBnYMzObueUTjSExw8Oao9s0Y1\nv1600rgeVrCy26CpMjd35oGpuByRv6Y1jbY5jx3AY8uLyHcbRndLCnQ4D5fKIsww/xvYYveQZ2EN\ntxXsSeBs/0DKf8eaFFW1AY/+XMS+csPwLvaejDskOVid52F5iD5flS3YxnRLivo8cuf2S6tW7HXS\nwtyQLQHv+6GQArdhfPekkEU6LofwzJgOgT5MYPXza25T8VACxZe/hifQCfDkqA5Rmzfu+sMzSHBY\n9bir6xgUIBSPgT8NSK81Ik4oZzZiVI+GHNbeRecUB9tL7OWi3/CXZpwd5WLLSkP2S8Ruv/gn/cX4\ntxwVmf6OAzokcP3hVZ8drem07H6bdwIn1Vg2GagcojoZiK0JiFRE7Cjx8oT/n+GOwfYqqqf2sYqV\nXgvRInGxjWLLSJrSu6rYa+luN+Pe2sXSOvpJPbumGIfA3UPqL9JJdgnzJnQM/H5uBG9wlYPsftCM\nHF1wseefBqSH/Um+PgdkuDivXxoG60GiPv8L6gs2/aBU29+/cBORRn1XizyGoZ0T6ZfVMgMdJ7uE\noTbnsiv2GCb2tD/3XVNcEzSKUTRKDcB+oLsSeF5EvhCRl0XkC6y56v7kXz8ULbpsE/7+YyElHsOk\n/ZNt/zNUjkzy+oZSaxqeIBU+ExgLL1RDlGga0SWRHaU+Ji/M5fm11Vtkeo3VkKG/jdxZcJ1lJP+Z\n+2a6OCTLRYG7aU/G20u8TF6YG/i9uaNxNMV1gzJIclrdUH4OMXLND7luTvmgKp2PjMxq0Vnkxzey\njrilcnOVguvpTv3A+m6H6t8a6e9AzSL/aLA7MsoHWDMVzMKql/sX0Ne/HGPMB/6hwlQEBY8XF9yq\nL1o2FHiY48/V3Pob+/8MQzsnsn+6k60lXr7YUb2Z83e73RRWGLLbuZo0/Uu4LTi+ExcemkaFD678\nIo9rv8wLrMtMFG6KwhBWjTWpiSONrNpXwW/f2V1tWLSmtGRtru5pzsBwXH/7X+26umW5bk55P7da\nMG+Jm2WwxuToUl0SKNVoKcGBdtHWcq78Io9+L+9gyn9zeXZ19bFXBzVjfM1YZbsg3hizxxjzgjHm\nfmPM88aYPQ3vpeoS3P9mX7n91nLPBH0hb/22IGrl25Xu/l8BHmMVGx3SiGIYEQnk6l7bUL34sqWL\nLWtKcAh/H5bFY8dk1RqN5frDM+jUiEGAo2WSjTnCavp8RznHv7ebLcVehjRhrM1wu2ZQBqkuYWGN\nusbKIJfvNtX6Uba0rqlO+re392B2Su8UMlqiBUaQrqlV39tHRmZxbI8kBGvA52u/qnqY+0sMPsiF\ng91+dH1EZJ6IrBSRzcGvSCcwHm0K6lP2fzZbLxZX+HhkeVV/r6W73bz1S2SGfwrljY2lJDnhpiMa\n/89wmr//2IJNpdVaNi7ZGjvFlsHOyU5j4aT96J5a9S9y0aHNG70/Uo7qlEDXFPs30jc2lDAtKHgs\nOL5TBFNnT+cUJxcfWr31ZGWQy3MbJu+fzLNjI9t5vbHGBRVffrmjnOfWFDPzmzymvZ/LYa/sCKyL\ndt+5hpzXL43/HNeJddO78cQxWYFxU4EWq0eMNLv/HfMAH1aXgnNrvFQjbSuuysX9a1UR22wM4fTs\nGqsTZ7Dbv8vH3ZieoGFw4SHpTSpiPCQrgYEdEsh3m2oNJ77PdeMUOKZrbAU6sCYLXXJy1WjyLVGs\nZ4dDpNrNqj6P/lzIBZ/sw+2Diw5NY864DrYn04y0Pw1IJyOhKi2VQW7S/snMHtuh1kAELa1yODCA\nSQtzufrLPGatLObjbeVsDWq5a7d1crRlJTk4KzuNl4/t2PDGrZzdQHcYcJ4xZqEx5pPgVyQTF6+2\nBf0TlHkJ9M0KpcTjqzUrc792LjYWeqsVZ0ZaZoJw7aCm52pO71u79aXXwNH7JTbY4bylRHO4puao\nr/iyxFP1gPRX/1iDdx2dyf1D2+Fs4bquYB2SnVx2WNX3K89tOL5XMs/FYJADAiOkgJWrPvPAFG7/\nTSbzJnTg+2ldAuui1bJQhWb30fxT4Ejg+wimpc0IzsE5xepQesWAdLJDjIH47Opidpf5OLJTAj/k\nWg0H7hicyfSP9vLAjwVMPyi12qghjRFczzds/k56pjmtV7rL/7PqRn/VwAw6NKOO6nd9U7n1uwLe\nr9EUPlbq51qz4FZ15y/ey+4yLztLfews8VJQUfU3TnRY433+rm9sFadVuuywdO73Dwk2sVcyc8bF\nZpADquWEPz7J/szwKvrsBrpNwH9FZD6wI3iFMebWcCcq3gV3SD4nO5U5a0u4+38FzBlXuwihxFNV\nN3fjERmcuciajuT4XsmM7JrIFzvcPPRTIXfUMbp/Q9xew/VfV1VEr87zsDov9FBMl/Rv+hiBYLWu\nO6ZrIp/VaHkZXASkmiY56KY7v8ZsEYkOAnPgzZ/YiZExWExcKXiw5efHdYjZ4mLVutgNdGnAO0AC\n0CtyyWkbgosubzwik1fWl7BgUxn/2+3mqBot4J5bU8KuUh9HdExgYtA4eSLC3Ue3Y9zbu5m1qog/\nHNq4ILS3zMu5i/dWa+7/+ZTObCn2sKXIy5Zi/6vIy9f+iUrDMYTVaQemVgt0mYnCUc0ce1BV96/R\n7emS4qBLipOuqU6yEoX2z20DiOkgV5MGORUutgKdMeb8SCekLQkuuuye5uTiQ9N5eHkRd3xfUK0F\nXKnH8PDPVjHOjUdk1CrrP7JTIqf2TeE/G0qrzXbdkNV5FZy5aA+bCr10TXGwo9R63B/QIaHOSUGz\nZm9t1PnV5+QDUvjzV3mBHMaorkktMv5fPIvmfGdKtQYhH9FFpHfQ+76hXtFIZDwxxtQaS/HqQRm0\nSxQ+2V5ebaDY59YUs7PUx+EdE0K2qvvrUZkkOqw5xOz4cEsZx72zm02FXo7omBD1uoWsJAcTg84l\n1roVKKXiT31lUT8HvV8H5Ph/Br9yIpe0+LSn3EfNGVXaJzm42j84bfAMvPXl5iodkOHi4v72W0Ke\nsWgPBRWGU3qn8N6kTnRPi36rwlODGkKMa8Z0O0opZUfIQGeMyQh67zDGOP0/g1+to+11DAnVZ+7i\n/ml0TXGwbE/VcEw7Sn0M6pDACQ30kbpuUEa12a7BasTy1c5yHl9RxB+W7A0s9xkrcD47tj2prpZp\n0h9c19g3U79CSqnIavnBBduYbSGmgEl1ObjxiEyuCRqOB+CGenJzlbKSHFx/RCY3+2e+PmbBLlbt\nq6hzVuFnx7RnWgs3LQ9uIah9jJRSkaZDgEXZ9uLQY1ue0y+VvhlVOZwBHRKYbHN8vz8GzXZdOUjv\ngA4JzOiXysMjsgLrWjrIKaVUtNnN0c0D1mMNAVb3pGLKlq0hcnRgDSh8y1GZXPCJNXtyfXVzNQU3\nxf7vpE4M6phQrWjyqi/z6tpNKaXint1Adxgw0hhjf6h9VadQs1dXOqVPSiDQ2c3N1TSsi7ZkVEqp\nSnZbI1QOAaaaqaEBnIMnk2zJiSWVUipe6BBgUdZQjk4ppVR46RBgUWZnSh6llFLho0OARVFhhY+C\nCkOy05qeRymlVOTZ7kcnIu2Ag4Fqw3AYYz4Od6Li1XZ/bq5bqpONhRrplFIqGmwFOhH5PfA4UET1\n7gUG0PEubaqsn+uepoFOKaWixW6O7h7gVGPMwkgmJt5tK7F6Z3RP1WGvlFIqWux2L3ABH0QyIW3B\ntqCiS6WUUtFhN9DdD9wiIi0zCnCcCC66VEopFR12iy6vAboCN4jInuAVxpj9w56qOLVVc3RKKRV1\ndnNo5wDHApOAc2u8bBGR40VkjYisE5GZdawXEXnEv/4nETmqoX1F5DQRWSEiPhEZbDctLaUyR9dD\nc3RKKRU1dvvRfdKcg4iIE6vV5m+BLcC3IvKWMWZl0GYnANn+11DgSWBoA/suB6YB/2pO+qKlcooe\nzdEppVT02O1ecGeodTaHABsCrDPGbPB/3svAFCA40E0BnjfGGOBrEckSkW5A71D7GmNW+ZfZOY0W\n5fYadpf6cAh0SdGqTqWUiha7dXQ1h/3qCowB5tvcvwfwa9DvW7BybQ1t08PmvjFvR6kXA3RNceBy\nxH5gVkqpeNHkIcBE5HhgethTFCE5OTlh2aapfixwAMl0cFT4j5NazzFjaV2spUfPMfzrYi09eo7h\nXxdr6WkorQ3fj7Ozs+tdH8z2EGB1+AB4xea2W6meK+zpX2Znm5oDSde1b4Mauig5OTmNunCNtXxj\nCbCPPh3SyM7eHz7fGjpdsbQu1tKj5xj+dbGWHj3H8K+LtfQ0kNZw34/t1tHVHOYrFTiL6kWK9fkW\nyBaRPlhB6kz//sHeAq7w18ENBfKNMdtFZLeNfWNe5ago3bTFpVJKRZXdHN06rHEtKyuXSoAfgBl2\ndjbGeETkCuB9wAk8a4xZISKX+NfPAt7D6r6wzv/559e3L4CITAUeBfYD3hWRZcaYiTbPKaoqB3Tu\noS0uVQM8S84GoOj8JS2bEKXihN1Al2CMadYoxMaY97CCWfCyWUHvDXC53X39y+djv0FMiwp0LdAc\nnVKqifLO79HSSWiVGmzn7u/HViQiSVFIT9wKDP+lOTqllIqqBnN0xhiviKwFOgLbIp+k+FQ5/JcG\nOqVaF81FtX52iy7nAu+IyMNY/dhM5QqdeLVhPmPYESi61M7iqm3IO79HRLvstGZ6baLLbqC71P/z\n9hrLdeJVG/aU+XD7ICtRSHVpoIsnesNS8aCpudbWktu122G8T6QTEs+2tfLpebQVYOvXWm5Iqun0\nbxyaZi+iYFsc1895lpwdCIRKKVUplu4NzRkZRdm03d9ZvLXm6JRSbVc85BQ10EXBNp1wVam4FA9B\noDnSZ4wFoGjOkhZNR0M00EXBthjoQ6f1bKqp2vrNvKli7brFUnqinRYNdFHQ2hujKBVNsXRDVvFB\nA10UbNeiS1VDvNzMY6lZupZaqFA00EVB5fBfPTRHpyKoOcGjtQReDWaqKTTQRVhhhY+CCkOy0+ow\nrhqmT/tKqXDSQBdh24P60IlooIukthDMWkvOq7XR6xrfNNBFWKjpedrCTbk+emOJLa2lmXi8SJ8x\nliPR6x0tGugibJtOuNoqRDvwamBRKno00EXYNv+oKNrisvVqau47XnKt8RKU4+XvoRpPA12EbQ/z\nzOL6z6piSbwEwVjS1qs1IkEHdY4wnXBVRUv6jLGBwKOUqqI5ugjb3oRRUfSJTsWDeGlwUV+uVXO0\nrYPm6CIsnqfoUSpetIXccFs4x1A00EWQ22vYXebDIdA5RS91W9OWbywqMvQ71TR6942gHaVWbq5r\nigOXo+11Fo/ExIuxNJmjUqp10Dq6CNJ56OJfvNRDqdYvUvWF8VAPqTm6CGpKQ5S2ItaKYGItPUqp\n8NFAZ0PW7K1kzd7a6P22ao5OxYhYC+Sxlp5IaAvn2FpooLOhvnqh+r7M2/2jooSzxaX+86h4F6nv\nuP7vxI5o/y20ji6CHl9RBESv6DJeRk2JhzoBpVTs0EAXQSPy1/Blu4OjVnTZ1ABRX4CU7Zurfikv\nhaSUJqRMhYs+BCjVeFp0GUFbkzoArXdmcdmxhZT7rgn8nnbNGSS+MRsK8lowVVHm9bR0CpSKuHgv\n1tVAFyE+Y9iW2B5onY1RZNc2Uu6/BkfenqplxQUkLphD2rWnkzTnIWTnlhZMYeQ51v5Eym0XB353\nLX4LfL7IHdBTQcJ7Lwd+TZj/HPi8kTueUm2EBroIyS3zUeFw0aGikBRXIzqLl5ZELlGABAWukNvk\n7rCC3N7dePsNDCwvufkRPEeMQCrcJHy8gNQbz6vaKY5yPpK/l6Sn7iX1nitx/ro+sDz5uQdJufsK\nHJvXhf2YzpX/I/WvfyTplVmBZUlvPkfyP25A8veG/XgRF/x98LTi70Z5aUunQIWB1tE1g8dnWJ/c\nmQ0pnVm9upiNhR42FXrYVOhlU6H1z92jfF/9H+Lz4ti4Fufyb3Et/w7H+hWBVckP/YXy6ZdhuvZs\nchodv24g8T9PB35PvfEc3CedQ8Vxp0JiElC93kf27iLlvmtx5O7Ee2B/Sq+9n/RLJllJ7TeIsn6D\nkK2bSPzvq7i++AC8Vg4n9apT8Q4ehWfIOLwHDwJn6/xqJSyaT+LrTyMlxRhXAhWTp5O44HkAfFkd\nca5fScptF1Fx3Km4p/6+2ceTvbtJfPlJEr752DpGl544/DllX0YWrhXf47j1Qsou/Su+Q45o9vEi\nzhicP35F4iv/DixKvX46FcedSsWYyZCaHp7jlBRBQiK4EkCaOepQUFBO/M/TSO4OHLu2IbnbceRX\n/f8mPX4HnuET8A4cYh1bBThyliMF+5CCPOtnYVX1hmPjGnx9Dm7B1Gmgs2WPK41yRwL7HriDAg8U\nVkC+BwoqwNN+ACXOJApeegaHM4VOzhSSXMl0d6bgwcEfdizG+bMT3OWI2w0V5Yi7PPDZaVecghQX\nBn43UpXJdi37EufPS6mYeCruk8+FlDTbaZadW0mcPxvX1x8hxlQtLysl6bWnSFjyDuVnXIJ38Oiq\ndftyrSC3exvePgdT+ucHICW11mebHr0p/8MNuKddQNrVpwLgKMzDsfhtEha/jS8jq9rnUlwITic4\nnFU/myLohuT64gMoL0XKSpGyEiirevJOeupeK9A6nJgax0t49yVMWgYmNR1S0zFp6db7yn1feBgA\nz8AhlJ97JaZLz0CgK7l3DolvPGsFw/++imvpkqadR5DUm85DykoxiUm4Tz6XiuNPJ/2PxwFQetfT\nJD95F841P5Jy37W4f3dBs48XSY71q0h6ZRbONT9WX753N0kvP0nim3OoGHui9ZDVGMbg2LwO19cf\nBxalX3pi1eqExGqBJ+mZB/AMGor3sN+EDqw+H461P5Pwzcc4v/0ksDjx7RerH9qVgHgqAEhYupiE\npYsxqel4Bo/GM2xC484DoLgQ5+plgV9dHy3A128gvh69wRH+AjbZl4tj4xqcG1dXLSwqgPTMRn+W\n87vPcOzcEnjJzqq+xal3XxFyv9TbL8Z7YH8qJpyCZ8jYFnlI0EBnw/cZfTlu3890X7G4aR/wj6Uh\nV0lxIb79uuMdMBjPgMF4Dz2S9MtOAqBi1AkkfLaQxPdexvX5+7hPu9DW4ZKeexDXp+8iXi/GlYB7\n3Mkkfvg6AKU3/IPEuY/h3LqJlMduw3vI4YH9Uu6/FsfOLXj3P4jS6//R4NO3ad8p8L7krqdxLV2C\na+liHDu34lj8VmBd5fnUJfneqzHd9sfX/QB8/p+Ba7N1E86Na3BsWmv93JxTtd+//xbyMxM+fz/k\nuqRX/1XvOfk6dqH8rCvw/uaY2jmF1HTc51yJZ+RxJD33IM5Na6tWXXsGJikFklMwydbPQFofnImU\nlkBZiRWUg4qnpawUz1HHUH7W5Zj9ulU7nGnfidIb/0niG7NJfGcuSUE5c9m727rZbN+MY8cWHDt+\nDaxLu/D4wANFzUCfOvNcTOUDh/9hAGfV+uQHZ4LThXG6wOWqljN3ffEBvi49MJ27YzKyql2fpMfv\nIGGp9f9h0jJxTzmXpHmPA1B6zb0kLHwF1+plJP73VRI++E9gP8eWjfjatYe0zFo3eseWjbi++RjX\nN4sDudzAtUlJgwo34qlAKtxQ4Q6sS/j0PRI+fQ/jdOLLHmAFvYFDA+sT5z6Ga+kSHHm51OSeMgNf\n5274OnXDdO6GyepE+vnjASg/4xJcXy3CuXld4BiBz3z13/g6dsF07Izp2AVfx86Bdc4fv8G56n84\nV/6AY3NOtQfP5Ocfss4nNQ3vQQOqVRXIzq1IUQFSXGD9LCqoOt5/nsYkJUNSMiYxuVpr6IQFz1v/\nLxvX1HmOaVefimfIOCrGn4zvwP5154gL8nAt+xLX958FFqU8+tfa2/l5DzwUk9Eek5mFyWyPadee\npLmP+c8tHef6lTjXr8T30hN4xp4Y8nMiRUzQRY83+fn5tk8uJyeH7OzsOtd9ddUVfJ/Zh84H9qFz\nktA5WdgvSeiUJGQ8/yAA7innQWmxlcMoLYbSElzLvwXAc9hgSEyynjwTEyEhiQR/ICh+YC6mS/Xm\n/cFFiY6Nq0l68TGc65ZX26bs/D8jxYX+fwLrp+u7TwPrjTjwjDwO99TfYzp1rd4s3evBteRdkt54\npto/D4C3Z19KZz4IGVl1pqemWusqn76XLiHxnbnWotQ08HqthhVeLxKGBh0VwyZAcmogqJjkVJJe\nftK6Nn+4MehYXvB6SHrpCQDcx5+OlBQhJUVQUoQUFyElhTh2b7fO498La3WhqPP8fV4SPlpA0ouP\nNOs8Sq+9D+/hwxo8nvPHb0j+9z21/l4txSSn4OvcHefmqjpMk5BoFelOng5pGbXOw7FxDQn/fRXX\n0sW1vgNGHJjMdpiM9ji3bKh1PF9mezxHjyHxozerfSY+H3gqoMIdeKAqP/VCXD99jWPdinq/a75O\nXfEMHYdn6HhSb72w+ucGqXkesu0XEr7+CNfXH+HY2bgRk4zThe+gwwK53orhx+Jc+zOOPTsb9Tm2\nj5eShrfPwfh6H0ziey/VWu/t2RfPuJMCJRnl0y/H9f1nVlGkqX7tPIcPsx50uvTE16UHvi49Sbv+\nLKCB6/bv/+L6ahEJH82v9n0JtR/Ufz+uS7t27eotv9YcnQ2/zVvOb/OWU3TXktorKwPdtNrFSpV/\n6LIb/lFrXWWgqxnkavL1OYTSWx7F9fXHJL46C8fe3QAkz679mZU8g0dTPu0CTI/edW/gdOGZMAXP\nsPEkvjmHRP8Ttrd7b8pu/Ge1INdoIvgOyMZ9QHYg0BU/+W71bYwh/ffjACi97n4c2zfj2LYZx/bN\nyLZfcPjL930du+DrczDe3v3w9T4Yb59+pF8+BYDyS2s/XVYGOs/oE2qvqwx00y+rM9mBptV2+wk6\nnFT8dlog0BX/4yWkvBTKrOJUykoDT8ClV//NCsYpaZjkVEhOIe2q3wHUCnKheA8fSsldT5N2zekA\nmIx2+Lr0wtetF76uPfF13T9wvKJ/Lww8WIjXCvhp15xmpfPeOdYy/4MAXg94vaTee1VVWr0exOux\nGpF4vSQ/cz8AFUPG4di1FceurUhJcbWbVsWoE3BPPR8TlJOpydfnYMov/Svu0y4k7bozrWXdeln1\nOsWFSP4+CKoTM2mZVhHh0HFWyYPTFQh0VX8Hh1XX7K9vBqg46WwqTjrbKiZc8T2un77B+fPSQAti\n98TT8Awdh6/voU2q3zPdD8A97QLcU88PfI/Lp12AY88uZO9O62fuTsRdBoD3wP54Dz0Sb/8j8R40\nAJKSA9+38ktuAUD27MKZ8zOOnOUkLppvXZtOXTHpmZi0TOtnembg/MunXWB938rLkPIyKC8L5Kjd\nx52Kr+8hePscjOncI5BTrgx0xQ/MJeGTd3B9uhDnlg04/UEOIOklKxdunC48/Qfj+c0xJD9n3ePK\nrr2v0dfK+tBkPGNPxDNmMo6c5SR89CYJX3/UtM9qIg10rYEInuET8Bw1gvSLrJt4xehJVl1TWiYm\n3fqZ8vjtAJT96U57n5uWgfvsKwKBruymhzCZ7SNxBtUF3Vy8g4biHTS02urKm0DJg69EPi1hYvbr\nRqjiA++RI8JzjA5VQaT4sQWhNwwK1jXTZLofEDKdECKt/kBXfvlt/g8xUFyAY+c2Uu+81Fr3xxvr\nS3r1NHTqGnhfct8L1htPhRXwCvMCuaviR96wik+bKi0D75CxeIeMrfZw5T7r8qZ/ZrCg73HFlPOq\nrwt+mLv1iQY/ynTsjKfjBBg2IRDoSv75cq3tKgNdreNBVaA7O3R9GVgP1+7TL8Y97QJc33+Ga/Hb\nuFb9YH3u0HF4fzMKz6ChVW0C/IGu2UTw9RtIeb+BGuhUPYJuYOV/uKH2+seb9/FRCXKq9ROB9Hb4\n0tuF7zNdCZgO+2E67Be0LIy3p+a2zIz14zWFKwHP0PF4ho6vymFedlvLpilCotaPTkSOF5E1IrJO\nRGbWsV5E5BH/+p9E5KiG9hWRDiLyoYjk+H/qnVoppVQ1UQl0IuLEym+cAPQHpotI/xqbnQBk+18X\nAU/a2Hcm8JExJhv4yP+7UkopFRCtHN0QYJ0xZoMxxg28DEypsc0U4Hlj+RrIEpFuDew7BZjjfz8H\nOCXSJ6KUUqp1iVYdXQ/g16DftwBDbWzTo4F9uxhjtvvf7wC6hEpATk5OqFUNbnNkPevrW8ctT1Wu\nDHeaKEIAAAoxSURBVNtnRmJdfemMtbRG5PxjLD16/uFfF2vp0fOvP60NrQMa1f0gbhqjGGOMiIRs\nUNbQRamv34ZrrNVMPi+7dleAyn4g9i+5/XRFe11DYimtkTr/WEqPnr+eY7jXxVp6Qq1rbD+6hkSr\n6HIr0Cvo957+ZXa2qW/fnf7iTfw/d4UxzUoppeJAtHJ03wLZItIHK0idCZxVY5u3gCtE5GWsosl8\nY8x2Edldz75vATOA+/w/6+lc1HStZebueJmMM17OIxL02ijVeFEJdMYYj4hcAbwPOIFnjTErROQS\n//pZwHvAJGAdUAKcX9++/o++D3hVRP4A/AKcHo3zibS2cDOLtXOsLz2xltamaIlzKJqzxCqCCrGu\nvv0ilR5Vt3i/NlGrozPGvIcVzIKXzQp6b4A6hyyoa1//8j1AE4YQb73i/QvZkFg7//pu5qr1iETg\njbXvalNF4jyifW3ipjGKij3x8o+uVGvQ0P9bW/5/1ECnVIS05RuLUrFEA51SSrUS+vDUNBrolFKq\niTTwtA4a6JSKE3rTVapuUZu9QCmllGoJGuiUUkrFNS26bCGtqZhJ+4qFX2v7+6vwampnetU0GuhU\n3Ij30U0iqa13ilbxTQNdG6E3JBUP9HusmkIDnVIqLmgQVKFoYxSllFJxTXN0qln0KVopFes00Cml\nlGqS1tIATIsulVJKxTXN0SkVY2LpSVipeKA5OqWUUnFNA51SSqm4poFOKaVUXNNAp5RSKq5poFNK\nKRXXNNAppZSKaxrolFJKxTUNdEoppeKaGGNaOg0Rk5+fH78np5RSCoB27dpJfes1R6eUUiquaaBT\nSikV1+K66FIppZTSHJ1SSqm4poEOEJHjRWSNiKwTkZktnZ5oEJFnRWSXiCwPWtZBRD4UkRz/z/ZB\n627yX581IjIxaPlvRORn/7pHRKTeSuHWQER6ichiEVkpIitE5Cr/8jZ/fUQkWUSWisiP/mtzh395\nm782lUTEKSI/iMg7/t/12gAissl/TstE5Dv/suhcG2NMm34BTmA90BdIBH4E+rd0uqJw3qOBo4Dl\nQcseAGb6388E7ve/7++/LklAH//1cvrXLQWGAQIsBE5o6XMLw7XpBhzlf58BrPVfgzZ/ffznke5/\nnwB84z+/Nn9tgq7RtcA84B3/73ptrHPaBHSqsSwq10ZzdDAEWGeM2WCMcQMvA1NaOE0RZ4z5FNhb\nY/EUYI7//RzglKDlLxtjyo0xG4F1wBAR6QZkGmO+NtY38PmgfVotY8x2Y8z//O8LgVVAD/T6YCxF\n/l8T/C+DXhsARKQnMBl4OmixXpvQonJtNNBZN7Bfg37f4l/WFnUxxmz3v98BdPG/D3WNevjf11we\nN0SkN3AkVs5Frw+BorllwC7gQ2OMXpsq/wfcAPiClum1sRhgkYh8LyIX+ZdF5droDOOqTsYYIyJt\nukmuiKQDrwNXG2MKgqsC2vL1McZ4gSNEJAuYLyIDaqxvk9dGRE4EdhljvheRsXVt01avjd8xxpit\nItIZ+FBEVgevjOS10RwdbAV6Bf3e07+sLdrpLxrA/3OXf3moa7TV/77m8lZPRBKwgtxcY8wb/sV6\nfYIYY/KAxcDx6LUBGAmcLCKbsKpAxovIi+i1AcAYs9X/cxcwH6vaKCrXRgMdfAtki0gfEUkEzgTe\nauE0tZS3gBn+9zOABUHLzxSRJBHpA2QDS/1FDgUiMszf8um8oH1aLf+5PAOsMsY8GLSqzV8fEdnP\nn5NDRFKA3wKr0WuDMeYmY0xPY0xvrPvIx8aYc9Brg4ikiUhG5XvgOGA50bo2Ld0SJxZewCSslnXr\ngZtbOj1ROueXgO1ABVY59x+AjsBHQA6wCOgQtP3N/uuzhqBWTsBg/xd2PfAY/kEIWvMLOAarPuEn\nYJn/NUmvjwEYBPzgvzbLgVv9y9v8talxncZS1eqyzV8brFbtP/pfKyrvs9G6NjoyilJKqbimRZdK\nKaXimgY6pZRScU0DnVJKqbimgU4ppVRc00CnlFIqrmmgU6qJROQvIvJ0w1uG9Zh3i0iuiOywub0R\nkYMinS6lYpl2L1CqlRCR/bH6FB1grNEl7OxjgGxjzLoGtusNbAQSjDGeZia1xY+jVDDN0SnVeuwP\n7LEb5JRSFg10SjVARG4Uka0iUuifBHKCf/nt/rEMEZHHRKQo6OURkdv967qLyOsisltENorIlfUc\nq52IPO/f9hcRuUVEHCJyLPAh0N3/+c+F2P96EdkuIttE5IIa6yaLNSFogYj8Wpk+v0/9P/P8nz9c\nRA4UkY9FZI+/uHRu5fBfDVwXh4jMFJH1/n1fFZEOoY5j64+gVDNooFOqHiJyMHAFcLQxJgOYiDWB\nZDXGmCuMMenGmHSsIcT2AQtExAG8jTX0UQ9gAnC1BM2YXMOjQDusIZPGYI3ld74xZhFwArDNf5zf\n15HW44E/Y40/mQ0cW2OTYv/nZWHNmXapiFTO5TXa/zPL//lfYU1seS/QHTgUa5Dd221clz9hzRE2\nxr/vPuDxeo6jVERpoFOqfl6sWY77i0iCMWaTMWZ9qI1FZD/gTeBPxpgfgKOB/Ywxdxpj3MaYDcBT\nWIP+1tzX6V9+kzGm0BizCfgncK7NtJ4OzDbGLDfGFOMPSpWMMUuMMT8bY3zGmJ+wxjsdE+rDjDHr\njDEfGmvyy93Ag0Hb13ddLsEay3CLMabcn45TRUSnBVMtQgOdUvXwN+K4GutmvUtEXhaR7nVtK9bU\nPv8B5hljXvYvPgCruDGv8gX8hf9v7/5dowjCMI5/HyyiwqlFVFDMgY0gCiJoLKwsRbD2R2EnCDaC\nhSKKQkQ7xSaFWIUUKexEsNBGQfwHbEQJwR9gEC4RTeVr8e7h5NhsrjmL5fnAcTAzO3szxb7MMHvv\nvwSTpXEyY/d8UTbP8Ek3d7E6WWXZD5ImJb2qtkV7ZEAaX6szSTur8X6WtATM9NuvMy9dMk9df7zv\nycBYN2azkXOgM1tHRMxGxHHyAR7A/TWaPgKWgBtF2QLwKSK2FZ9ORJysuX6RzCbRLcomGD4X2VdW\n5/CaGKifJdOf7ImIrcA0uT0JOa5Bd6vygxGxBThftG+alwXy3+bLMW+MzEfmY9723znQmTWQtE/S\nCUljwArwG/hT0+4iua13LiLK+nfAcnVwY5OkDZIOSDoy2Edk5u45YEpSR1IXuEKupIYxB1yQtF/S\nZuDWQH0H+BERK5KOAmeLuu/VuPYOtP8J9CTtBq4W422al+lqDN2q7XZJpxvuYzZSDnRmzcaAe+Rq\n6xuwA7hW0+4M+fD+Upy8vF4Fr1PAIfL9sUXgMXngpM5l8tDIR+A1uQp7MswPjYjnwAPgJfCh+i5d\nAu5IWgZukoGxf+0vYAp4U205HgNuA4eBHvAMeFr01TQvD8mV44vqXm+ByYb7mI2UXxg3M7NW84rO\nzMxazYHOzMxazYHOzMxazYHOzMxazYHOzMxazYHOzMxazYHOzMxazYHOzMxazYHOzMxa7S/sdesx\nKnTMyQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f5e9111a278>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.errorbar(fast_time_cov[:,0], fast_time_cov[:,1], fast_time_cov[:,2], label='fast covariance', linewidth=2)\n",
    "ax.errorbar(slow_time_cov[:,0], slow_time_cov[:,1], slow_time_cov[:,2], label='naive covariance', linewidth=2)\n",
    "ax.set_xlabel('size of dataset')\n",
    "ax.set_ylabel('running time')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Affine Transformation of Dataset\n",
    "In this week we are also going to verify a few properties about the mean and\n",
    "covariance of affine transformation of random variables.\n",
    "\n",
    "Consider a data matrix $\\boldsymbol{X}$ of size (N, D). We would like to know\n",
    "what is the covariance when we apply an affine transformation $\\boldsymbol{A}\\boldsymbol{x}_i + \\boldsymbol{b}$ with a matrix $\\boldsymbol A$ and a vector $\\boldsymbol b$ to each datapoint $\\boldsymbol{x}_i$ in $\\boldsymbol{X}$, i.e.\n",
    "we would like to know what happens to the mean and covariance for the new dataset if we apply affine transformation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def affine_mean(mean, A, b):\n",
    "    \"\"\"Compute the mean after affine transformation\n",
    "    Args:\n",
    "        mean: ndarray, the mean vector\n",
    "        A, b: affine transformation applied to x\n",
    "    Returns:\n",
    "        mean vector after affine transformation\n",
    "    \"\"\"\n",
    "    affine_m = A @ mean + b # EDIT THIS\n",
    "    return affine_m\n",
    "\n",
    "# ===YOU SHOULD EDIT THIS FUNCTION===\n",
    "def affine_covariance(S, A, b):\n",
    "    \"\"\"Compute the covariance matrix after affine transformation\n",
    "    Args:\n",
    "        S: ndarray, the covariance matrix\n",
    "        A, b: affine transformation applied to each element in X        \n",
    "    Returns:\n",
    "        covariance matrix after the transformation\n",
    "    \"\"\"\n",
    "    affine_cov = A @ S @ A.T # EDIT THIS\n",
    "    return affine_cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the two functions above are implemented, we can verify the correctness our implementation. Assuming that we have some matrix $\\boldsymbol A$ and vector $\\boldsymbol b$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "random = np.random.RandomState(42)\n",
    "A = random.randn(4,4)\n",
    "b = random.randn(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we can generate some random dataset $\\boldsymbol{X}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = random.randn(100, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming that for some dataset $\\boldsymbol X$, the mean and covariance are $\\boldsymbol m$, $\\boldsymbol S$, and for the new dataset after affine transformation $ \\boldsymbol X'$, the mean and covariance are $\\boldsymbol m'$ and $\\boldsymbol S'$, then we would have the following identity:\n",
    "\n",
    "$$\\boldsymbol m' = \\text{affine_mean}(\\boldsymbol m, \\boldsymbol A, \\boldsymbol b)$$\n",
    "\n",
    "$$\\boldsymbol S' = \\text{affine_covariance}(\\boldsymbol S, \\boldsymbol A, \\boldsymbol b)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X1 = ((A @ (X.T)).T + b)  # applying affine transformation once\n",
    "X2 = ((A @ (X1.T)).T + b) # and again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One very useful way to compare whether arrays are equal/similar is use the helper functions\n",
    "in `numpy.testing`. the functions in `numpy.testing` will throw an `AssertionError` when the output does not satisfy the assertion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correct\n"
     ]
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(mean(X1), affine_mean(mean(X), A, b))\n",
    "np.testing.assert_almost_equal(cov(X1),  affine_covariance(cov(X), A, b))\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fill in the `???` below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "correct\n"
     ]
    }
   ],
   "source": [
    "np.testing.assert_almost_equal(mean(X2), affine_mean(mean(X1), A, b))\n",
    "np.testing.assert_almost_equal(cov(X2),  affine_covariance(cov(X1), A, b))\n",
    "print('correct')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check out the numpy [documentation](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.testing.html)\n",
    "for details.\n",
    "\n",
    "If you are interested in learning more about floating point arithmetic, here is a good [paper](http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.22.6768)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "coursera": {
   "course_slug": "mathematics-machine-learning-pca",
   "graded_item_id": "YoDq1",
   "launcher_item_id": "vCPZ0"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.0"
  },
  "widgets": {
   "state": {
    "9b7b6bf11f38461ab695892799fb14b8": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}