pca/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": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/matplotlib/__init__.py:1401: UserWarning:  This call to matplotlib.use() has no effect\n",
      "because the backend has already been chosen;\n",
      "matplotlib.use() must be called *before* pylab, matplotlib.pyplot,\n",
      "or matplotlib.backends is imported for the first time.\n",
      "\n",
      "  warnings.warn(_use_error_msg)\n"
     ]
    }
   ],
   "source": [
    "# PACKAGE: DO NOT EDIT\n",
    "import 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": 28,
   "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": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "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": 30,
   "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 0x7fb64574c780>"
      ]
     },
     "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": 31,
   "metadata": {
    "collapsed": true
   },
   "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",
    "    mean = np.zeros(D)\n",
    "    for n in range(N):\n",
    "        mean += X[n, :] # EDIT THIS\n",
    "    mean = mean / N\n",
    "    return mean\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",
    "    covariance = np.zeros((D, D))\n",
    "    mean = mean_naive(X)\n",
    "    for n in range(N):\n",
    "        diff = np.asmatrix(X[n, :] - mean)\n",
    "        covariance += diff.T @ diff # EDIT THIS\n",
    "    covariance = covariance / (N-1)\n",
    "    return covariance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "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, bias=True) # EDIT THIS\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": 33,
   "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 0x7fb64576ebe0>"
      ]
     },
     "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": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 112 ms, sys: 0 ns, total: 112 ms\n",
      "Wall time: 131 ms\n",
      "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n",
      "Wall time: 4.17 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": 35,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.04 s, sys: 8 ms, total: 3.04 s\n",
      "Wall time: 3.53 s\n",
      "CPU times: user 20 ms, sys: 24 ms, total: 44 ms\n",
      "Wall time: 42.3 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": 36,
   "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": 37,
   "metadata": {
    "collapsed": true
   },
   "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": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAEGCAYAAAAT/1CLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5+PHPM1tCCEsRRAQUinHfQGURQRZRQC2Cdf1W\nAa0tVm1rW1utrVqtdavaavlJq1JRQRQVQUVRi7i0omhdCqIGEJQdFwJZZybz/P64dyYTyEwmyWQm\nmTzv12temTnnnnvPPQz3mXvuveeIqmKMMcbkKk+2K2CMMcY0Jwt0xhhjcpoFOmOMMTnNAp0xxpic\nZoHOGGNMTvNluwLNqaSkxG4pNcaYHNepUydJlm9ndMYYY3KaBTpjjDE5zQKdq7i4ONtVaNGsfRKz\ntknM2iYxa5vE0t02GQt0IjJWRD4VkdUicnUd+SIi97j5H4nIgPrKisjjIvKB+1onIh9kan+MMca0\nDhm5GUVEvMB0YAywAVguIgtV9eO4xcYBRe5rEHAfMChZWVU9J24bdwIlmdgfY4wxrUem7rocCKxW\n1bUAIjIXmADEB7oJwMPqDL65TEQ6i0gPoE99ZUVEgLOBUalURlUpLS0lEonE0vLz8ykpsTiZSDra\nx+PxUFhYiPPPZYwxmZGpQNcT+DLu8wacs7b6lumZYtlhwFZVTaljt7S0lLy8PAKBQCwtLy+P/Pz8\nVIq3Selon2AwSGlpKR06dEhTrYwxpn658hzdecBjyRaIv7iZn59Pt27dqKysrLXM7p9Nbelonx07\ndrBly5Y01KZlsRsLErO2SczaJrH62qaoqCjldWUq0G0Eesd97uWmpbKMP1lZEfEBk4BjklUgvlFK\nSkr2ODuprKxMesbS+Z/OJndM7ZlsMzmrvvZJVceOHendu3f9C7YixcXFDfpP15ZY2yRmbZNYutsm\nU3ddLgeKRKSviASAc4GFuy2zELjQvftyMFCiqptTKHsS8Imqbmj+3UifGTNmMHDgQC655JIGldux\nYwcPPPBAM9XKGGNyT0bO6FQ1LCKXA4sBLzBTVVeKyDQ3fwawCBgPrAbKganJysat/lzq6bZsiR58\n8EGeeeYZevZs2BliSUkJDz74ID/84Q+bqWbGGNO8CiePAKB01tKMbC9jz9Gp6iJVPVBV+6nqzW7a\nDDfIoY7L3PwjVPXdZGXj8qZE19FaXHnllaxbt46zzjqL6dOn89577zFmzBiGDRvGySefHOubXrVq\nFaNGjeKEE07g+OOPZ82aNfzhD3/g888/54QTTuD3v/99rfWuX7+e4447jksvvZRjjjmGSy65hKVL\nl3LKKacwYMAA3nvvPQDKysq47LLLGDVqFMOGDeP555+PlR83bhzDhw9n+PDhvP322wC88cYbTJw4\nkQsvvJDjjjuOSy65BJuZ3hjTWuTKzShNEr3+lq5l67uOd/fdd/PKK6/w7LPPstdee7Fz505eeOEF\nfD4fS5cu5cYbb+SRRx5h5syZTJs2jbPPPptgMEh1dTXXX389q1at4s0336xz3WvXruWhhx7ikEMO\nYeTIkcybN48XX3yRRYsWceeddzJnzhzuvPNOhg8fzvTp09mxYwejR49mxIgRdOvWjfnz55Ofn8+a\nNWu4+OKLWbp0KQArVqxg2bJl9OjRg1NOOYVly5YxZMiQlNvNGGOyxQJdC7Bz504uvfRS1q5di4gQ\nCoUAGDhwIHfeeSebNm3i9NNPp1+/fvWua//99+ewww4D4OCDD+bEE09ERDjssMP44osvAFiyZAkv\nvPAC9957LwBVVVVs2LCBffbZh6uuuooVK1bg8XhYs2ZNbL39+/ePdbMeccQRfPHFFxbojDEJZbp7\nMhkLdDhnYNm86/Lmm29m2LBhzJ49m/Xr13PaaacBcNZZZ3HssceyePFizjrrLO6++2769OmTdF15\neXmx9x6PJ/ZZRKiurgacB+YffvjhPe5quuWWW9h777158803iUQidO/ePZYX/8yh1+slHA43aZ+N\nMSZTbFDnFmDnzp306NEDgDlz5sTS161bR58+fZg2bRrjx49n5cqVdOjQgV27djVpe6NHj+Yf//hH\n7Drbhx9+GKtH9+7d8Xg8zJ07NxYYjTGmNbNA1wL87Gc/48Ybb2TYsGG1gsv8+fMZMmQIJ5xwAqtW\nreK8886jS5cuDB48mCFDhuxxM0qqrrrqKkKhEEOHDmXw4MH86U9/AuCHP/whjz32GEOHDqW4uJj2\n7dunZf+MMSabJJfvnks0w3hJSQmdOnWqlWYPjCeXrgfG62r71s4e/E3M2iaxXGibZNfhGpsHDW+b\n+mYYt2t0KWqrAc4Y07a1pJtKGsu6Lo0xxuQ0C3TGGGNymgU6Y4wxOc0CnTHGtHGFk0fErsXlIgt0\nKcr1L4IxxuQqC3StxMyZM3nssVY3SYMxxmSdPV7QSlx00UXZroIxxjSa7Pga3ztL8S37V8a3bWd0\nWbB+/XoGDhzIT3/6UwYPHszEiROpqKgAYNasWYwcOZKhQ4dywQUXUF5eDjjjUN5777189tlnjBo1\nqta6jj/+eAA++OADxo8fz4knnsikSZPYsmXLHtu+9NJL+cUvfsFJJ53EUUcdxRtvvMFll13GwIED\nufTSS2PLLVmyhDFjxjB8+HAmT55MWVkZALfddhsjR45kyJAh/OxnP4sNI3bqqady/fXXM2rUKI45\n5hj+85//NE/jGWNaFd/ri8i//ZcU/Pws8mbfi3fNx5mvQ8a32AIVTh5BYQOWrU8qD1auWbOGBx54\ngHvuuYcpU6awcOFCzjnnHE4//XQmT54MwB//+EceeeQRfvzjH8fKHXjggQSDwdg4mPPnz2fixImE\nQiF+/etfM2fOHLp27crTTz/NTTfdxPTp0/fY9o4dO3j55ZdZtGgR559/Pi+++GJsWp+PPvqInj17\ncscdd/DMM8/Qvn17/vKXvzBjxgyuvfZafvSjH/Gb3/wGgB/96Ee8+OKLjBs3DoBwOMySJUt46aWX\nuO2221iwYEEKLWqMyTWy89vY+/wHbwdAvT7CRw4hPHgU+ffdlNH6WKDLkv33358jjzwSgKOPPjo2\nhc7HH3/MzTffTElJCaWlpYwePXqPshMnTmT+/PlceeWVPP300/zzn/+kuLiYVatWccYZZwDsMftA\nvLFjx8am7unWrVutaX2++OILNm3axKeffsopp5wCQCgUYsCAAQC8/vrr3HPPPVRUVPDtt99yyCGH\nxALd6aefvsf+GGPakEgE32vPk/fE32NJ4UP6Ex48mvCxw6Gwo5NogS7zSmctrXcsx3QPgxM/nY7X\n6411Xf7kJz9h9uzZHHHEEcyePbvOCVYnTZrE5MmTOf300xER+vXrx8qVKzn44IN5+eWXU962iOwx\nrU91dTVer5eRI0fy4IMPxvIqKyuprKzkV7/6Fa+++iq9evXilltuobKyco/12jQ+xrQ8zT2Ul+fL\nteQ9dBfe1StqpVdefXezbK8h7BpdC1NaWso+++xDKBRi3rx5dS7Tt29fvF4vt99+O5MmTQKgqKiI\nr776infeeQdwzsJWrVrVqDocd9xxvP3226xduxaAsrIy1qxZEwtqe+21F6WlpSxcuLBR6zfG5JbA\n43+n3fWX4F29gkinLlT+5LpsV6mWjAU6ERkrIp+KyGoRubqOfBGRe9z8j0RkQCplReQKEflERFaK\nyO2Z2JfmdO211zJ69GhOOeWUpKN3T5o0iSeeeCLWVRkIBJg1axbXX389Q4cOZdiwYbGg11Bdu3Zl\n+vTpXHzxxRx//PGMGTOG1atX07lzZyZPnsyQIUM488wz6d+/f6PWb4xJv4w/6xusir0NLHoMIhGC\no8+g/JZZhAeNSlIw8zIyTY+IeIHPgDHABmA5cJ6qfhy3zHjgCmA8MAj4q6oOSlZWREYC1wKnqmqV\niOytqtui60znND25MIJ3U9g0PYnlwnQrzcXaJrF0t019x6i0TKlTUYbvw2V4330D30fLkCqnl6d6\nv35UTfklkX6HNnl70Hqn6RkIrFbVtQAiMheYAMTfZzoBeFidyLtMRDqLSA+gT5KylwK3qmoVQHyQ\nS7e2GuCMMQYg/66r8a58DwmH9siruOHv4G25t3xkquuyJ/Bl3OcNbloqyyQreyAwTETeFpHXROS4\ntNbaGGPaqvJSAvPuj330fbgMqsNUH3gkVedfRtmdc2uWbcFBDlr/XZc+oAswGDgOeEJEvqt19McW\nFxfH3ufn59e62zAq/g5Cs6d0tM/OnTvZtq3ZTryzJv77ZWqztkmsoW3T/4+XAPD+7+7fM6+edSbL\nr5VXHabrf19nnzeexV9eGlvmi/EXUHLg0YSjjwh8uyv1dTawrvXlAQ3q2sxUoNsI9I773MtNS2UZ\nf5KyG4Cn3cD2johEgK7A9t0rEN8oJSUle1xvStc1qFyVrvbp2LEjvXv3rn/BVsSuQyVmbZNYU9om\nWbn61pks/+CSLeQ98Xc8WzcAUF10ON5i53GBLudcTJc01ydRXrq/N5nqulwOFIlIXxEJAOcCu9+b\nvhC40L37cjBQoqqb6yn7DDASQEQOBALAV/VVxuPxEAwG07FfpgGCwSAejz3RYkxL1e7e3+PZuoFI\n915UXHETFdfem+0qpUVGzuhUNSwilwOLAS8wU1VXisg0N38GsAjnjsvVQDkwNVlZd9UzgZkisgII\nApPr6rbcXWFhIaWlpbGHtMHpUuvYsWN6djgHpaN9PB4PhYWpDrZmjGl24TCBp2oGhtAOnQieMYXQ\niNPB19qvbNXI2J6o6iKcYBafNiPuvQKXpVrWTQ8CP2hoXUSEDh061Erbtm1bznWppZO1jzG5RbZt\nIv++m/CurRlYouz22VCQez9GcydkG2OMSYnv7SXk/fNOpKKMyF7d8Xy91cnIwSAHNgSYMca0KM09\nwknezDvI/383IhVlhI8ZRvlNDzTbtloKO6Mzxpg2xP/a86jfT9V5lxMe9T2QpIOKNItMD8Bhgc4Y\nYzKscPII+pOZA75nfTH+fz0T+xzZd38qL72OyH79mn3bLYUFOmOMyTWhIL53luJfsgDv6pW1sspv\nmAF57Zq9Ci1p2EQLdMYYk0MCT/wD/+vPI7tKANCC9oROGEvgpaecBTIQ5FoaC3TGGNPKyZYNsfeB\n5+cAzowCodETCQ8ZDXntagJdG2SBzhhjWinZvpnAgofx/XtxLC005CRCJ010pszJwo0mLZEFOmOM\naQbNOYelfLONwMJH8b2+CKkOo3FD61VN+13at5dIS7oOl4w9R2eMMa1Mwa//D/+rCyESIXT8yZTf\n+nC2q9Si2RmdMca0dKr43ngh9lFCIUKDRhI8Ywq67/5ZrFjrYIHOGGNaMNm0nvyH7sL76YextPKb\nHmxTz8E1lQU6Y4xpiYJVBJ6bg/+52Uh1mEjH7+DZ+S1ARoNca7kOl4xdozPGmEZqznEpC37/QwIL\nZiHVYUInnkb5LbOaZTttgZ3RGWNMSxCswvf+f2IfPVu+pHrfPlRN/QWRA4/MYsVaPwt0xhiTLap4\nVq/E/+/F+N5+FSkvjWVVnXkxofHngs+fxQrmBgt0xhiTRHM+D1fwmwvwbK0Z1aS6z4F4130GQOh7\nF6R9e4nkwnW4ZCzQGWNMplRVEHh+buyjZ+sGIp33Inz8GMJDTyHSq2+zzkXXVlmgM8aY5haJ4Pv3\nYgJPPoBnx9ex5Ipf3kb1YceA1w7FzSljd12KyFgR+VREVovI1XXki4jc4+Z/JCID6isrIjeIyEYR\n+cB9jc/U/hhjTCo8n3xAuxumkf/AbXh2fE1134NiedVHDspYkCudtTTnuygTyUgLi4gXmA6MATYA\ny0Vkoap+HLfYOKDIfQ0C7gMGpVD2blX9cyb2wxhjGqrglp8DEPlOV4Jn/YjwkJMonDoqy7VqWzJ1\nvjwQWK2qawFEZC4wAYgPdBOAh1VVgWUi0llEegB9UihrjDGNls4bTuTrrQQWPBL7rIF8gqeeR2jc\nOZCX3+T1m4bLVKDrCXwZ93kDzllbfcv0TKHsFSJyIfAu8EtV/TZdlTbGmFTJjq/xPzcb/6vPIuFQ\nLL38tkfQLt2yWDPT2q+A3gfcBKj7907goroWLC4urndlqSzTlln7JGZtk1hraJv+7t+66ppKXtmM\nW+i2/FU84SCK8O2hx/Gdj5cD8NnXO+DrHWndXp1t+rv7o5l75rVC9X1vioqKUl5XpgLdRqB33Ode\nbloqy/gTlVXVrdFEEbkfeC5RBeprlOLi4gY1XFtj7ZOYtU1ira1tktV19zz5Zlvsffe3nIlPw8cM\nIzhxKv7e3wW3O7Qh60xHXi5I9/cmU4FuOVAkIn1xgtS5wPm7LbMQuNy9BjcIKFHVzSKyPVFZEemh\nqpvd8hOBFc2/K8aY1igt1+FU8X78Hv5/LcD7/r9jyeGjBhOcOJVI3B2VpuXISKBT1bCIXA4sBrzA\nTFVdKSLT3PwZwCJgPLAaKAemJivrrvp2ETkap+tyHfDjTOyPMabt8b/0JP4lC/Bsdm4ZUK83llf5\ni1szUoe2+nhAU2XsGp2qLsIJZvFpM+LeK3BZqmXd9MyNkWOMaXt21lxby5v9N8B5TCA08nuETzyV\n9j87M1s1Mw3Q2m9GMcYYIM1jUqri+8/L5M35WywpfOgAQqPPoLr/8TaSSStj/1rGGBNHtm8m76G7\n8K1YXiu98jd3ZWT71j2ZfjbxqjHGuPwvzqPgt1PxrViOtu9A5Q9/k+0qmTSwQGeMMa68x6YjwUpC\ng0ZSfssswsPGZbtKJg1S6roUkTzgOuA8YC9V7SQiJwMHqurfkpc2xpgWqrSEwDOzYh8jXbpRNflK\nqo8+vnk3O2up86xYs27FRKV6je5unKG4/g94wU1b6aZboDPGZETabjgJh/EvWUDgmYeQsl2x5PI/\nzYJ2BU1bt2lxUg10E4EDVLVMRCIAqrpRRHo2X9WMMSb9vB8uI++x6bHn4cKHDsD38X+dTAtyOSnV\nQBfcfVkR6QZ8XffixhjTMrW7y5nSMtK9J1Xn/oTq/sdTOGVk2rdjd0+2HKkGunnALBG5Epyht4C/\nAHOTljLGmAZK6/Nw4RDe9/+D/80XY0la0J7ghCmETjoDfP6mb8O0eKkGut8CtwH/AwqAYuB+4MZm\nqpcxxjSOKp4vVuN74wX8b72ClO6slV12+2zo0DlLlTPZkFKgU9UgcCVwpdtl+ZU7ZJcxxrQo7X7/\nQ7xfrol9ru7dj/CwseTNme4kpDHIWfdk65DyyCgiUgAcABTizCYAgKr+p3mqZowxKaiqwP/y07GP\n3i/XoIUdCQ0ZQ3jYWCL7HQAiNYHOtDmpPkd3Ic5jBEGgIi5Lgf2aoV7GGJNcKIj/1WfxP/conpJv\nY8kVV9xI9VGDwR/IYuVMS5LqGd3twJmq+nJzVsYYY1Lhe+15AgsexvO1M/dydd+D8X7+ifP+2OFp\n3ZZ1T7Z+DXm8YGkz1sMYY5KLuy0gf+YdAFT36kvwzIup7j+0WR4RMLkh1UD3e+AuEfmDqn7VnBUy\nxuS+hj5C4Fn7CXmP1Vxji3TvSXDiRYQHjQSPDdlrkks10H2G8yjBT6I3oQCCM1+qN2EpY4xpAvlm\nG4F5D+D/z0u10sv/NAt8NsuYSU2q35RHgIeBx6l9M4oxxtSpcPII+tP4a1yBp/+J/4W5SLAK9fkJ\nnfx9AoseczLTHOTsOlxuS/XbshdwnT07Z4xpVuFw7G1ggTOrQOi4EQTP+THarUdNoDOmAVLt3P4n\ncEFTNiQiY0XkUxFZLSJX15EvInKPm/+RiAxoQNlfioiKSNem1NEYkyXhML7XF1Fwdc1hprrvQZRf\new9Vl9+AduuRxcqZ1i7VM7qBwOUici2wNT5DVeu9l1dEvMB0YAywAVguIgtV9eO4xcYBRe5rEHAf\nMKi+siLSGzgZ+CLFfTHGpEmTx6UMh/H95yUCCx/Bs31zrayK6+6zG01MWqQa6O53X401EFitqmsB\nRGQuMAGID3QTgIfd7tFlItLZHTy6Tz1l7wZ+DSxoQv2MMVlQcPWFeLZvAiDSozfBCZPJn/FHJzPN\nQc6uw7VdqY51Oav+pZLqCXwZ93kDzllbfcv0TFZWRCYAG1X1w7i7QY0xLZV7Bhfl2b6JyD69CU64\nkPDgUeDxQjTQGZMmCQOdiFygqo+47y9KtJyqzmyOitXHHXvztzjdlvUqLi5OyzJtmbVPYm21bfq7\nf+va//g8qQ7T5aO36P7vReTtqHkUd92Ei/n2sIHO2duatQ1a5x5+d380s+E7kiVt9XuTivrapqio\nKOV1JTujOw/nsQJIfCOKAqkEuo1A77jPvdy0VJbxJ0jvB/QFomdzvYD/ishAVd2yewXqa5Ti4uIG\nNVxbY+2TWK63TSrX4ZLt/yEbPiHw3KN4vnIu70d69I7N7t110gUkuoMs2Tpzob1z/XvTFOlum4SB\nTlXHx71v6tg6y3FmPOiLE6TOBc7fbZmFODe8zMXpmixR1c0isr2usqq6Etg7WlhE1gHH2sgtxrQA\nwarY2/yH7gQgsu/+BL93IeFBIyicOrpRq7XrbKYxUp294H1V7V9H+ruqemx95VU1LCKXA4sBLzBT\nVVeKyDQ3fwawCBgPrAbKganJyqa0d8aYzKoow79kIf7F82JJ1fv2IXTGhYSPO9G5BmdMhqV61+UB\nuyeI01/43VQ3pKqLcIJZfNqMuPcKXJZq2TqW6ZNqXYwxtTX5MYGdOwi89CT+fz2DlJfWyqq4eaY9\nJmCyKmmgE5GH3beBuPdRfQA7szLG0P6X5yBud2X1QUcRPO182t35GyfTgpzJsvrO6NYkeK/Av4F5\nGGPaHPl6K4FnH635HKwifPQQgqeeT+TAI7JYM2P2lDTQqeofAERkmaouzkyVjDEtlXyzjcCzs/G9\n9jxSXTMuZflNDxLZr18Wa2ZMYqk+MG5BzphWrMnX4IDAI3/Fv/Q5JBxCRQgNHo1/2b8AGhzkktXD\n7qw06WYTOhljEpJvtsfeB16ZD0Bo4EhCEy4k0qtvLNAZ05JZoDPG7EG2biSwaC6+N1+MpYWPHU7w\njClEeqd8s7UxLYIFOmNyRDq6Jz0bPsf/3Gx8y5YgGkHjxpCtvOLGJtbQmOxI9YHxRGNdVuEMsrxM\nVasSLGOMaSUKrp0KgHo8hIaeQvDU82l/zeQs18qYpkn1jO5CYAjOXHQbcMaV7A68i/M8HSIyQVXf\nbYY6GmOaiWzfTODxv8c+q99PaPiphMad0+TJTktnLXXGLGxqJY1polQD3UrgaVW9J5rgDst1MHAC\ncC1wL04wNMY0k3R0TwJQUU7gudn4Fz+BhEKx5PI/z0U779W0dRvTwqQa6M4Hdv/23wd8paqXi8gd\nwFVprZkxJv0iEXxvvkjgyfvxlHwLQGjISfjfegXAgpzJSakGuq3A6dSexftUYJv7Ph8I7V7IGNOy\ntLvhx3jXO/N8Vfc7hKrzLydywGGxQGdMLko10P0UmCciK3Bm++4NHA6c5eYPwum6NMa0MJ5PP4q9\n964vJtKlG8Gzf0x40Kgmj0NpD3eb1iDVkVFeEpHv4kyjsy/OTALPq+rX0XzgpWarpTGmYVTxrlhO\nYOGjeD+rCXRVZ0whNP5cyMvPYuWMyayUn6Nzg9oj9S5ojGmSpt5w4n33DQLPPoJ33WcAaEFhbOqc\n0MQpaaihMa1Lqs/R9QVuBo4GCuPzVHW/ZqiXMaYhwjUDLLe79/cARDp+h9DYswiNmkDhtFMbvWrr\nnjStXapndHNwpun5Jc7s38aYlqCqAv9ri/C/+EQsKdJlb0LjzyV04qkQyMti5YxpGVINdIcBQ1U1\n0pyVMcakaNcOAq/Mx//KfKR0Z62s8jtmg8+fpYoZ0/KkesvV60D/5qyIMW1J4eQRsWtxjdH+F+cQ\neGYWUrqT6n6HUvHTm2oyLcgZU0uqZ3TrgBdFZD6wJT5DVa9LZQUiMhb4K+AFHlDVW3fLFzd/PE73\n6BRV/W+ysiJyEzABiOA80zdFVTeluE/GtB7BKvwvPB77KMEqwkcOcmb0PuhIiBt8uTHsOpzJZakG\nuvbAc4Af5xm6BhERLzAdGIMzVuZyEVmoqh/HLTYOKHJfg3BGXhlUT9k7VPX37jZ+ClwHTGto/Yxp\nsVTxvvs6eXPvw/NVzW9Mm9HbmNSl+hzd1CZuZyCwWlXXAojIXJwzsfhANwF4WFUVWCYinUWkB86g\n0XWWVdX4ixPtAW1iPY1pMTxfriUw+158q94HoLp3P7xfrgEaPqO3MW1ZwkAnIn1UdZ37PuFMi9EA\nVI+eOCOqRG3AOWurb5me9ZUVkZtxZlcoAUamUBdjMqKpz8O1+/0PnTnh2nek6syLCY84lcKLTkpf\nBY1pI5Kd0f0P6OC+X41ztrT7hQDFuW6WNap6LXCtiFwDXA5cX9dyxcXF9a4rlWXaMmufxOpqm/4N\nzJNwiK7vvlrrQdXtx45k84kTqG7XHtZ+3uB1xvzu/mhm0v1oDva9SczaJrH62qaoKPUJoBIGOlXt\nEPe+aQPiwUZqX9vr5aalsszu1wXrKgswG2dosjoDXX2NUlxc3KCGa2usfeqWyllbsnYrKioCVXzv\nLCUw7x94tm+O5VX88QHa9foudXWn1LvOFsK+N4lZ2ySW7rZJeQiwJloOFLkjrGwEzsWZ+ifeQuBy\n9xrcIKBEVTeLyPZEZUWkSFWjYX8C8Enz74ox6eP57CPy5t6Hd80qAKr37YN30zoAIr0SXjEwxjRA\nRoYAU9WwO1HrYpyuzpmqulJEprn5M3DOxsbjdJOWA1OTlXVXfauIHITzeMF67I5L08oU3PxTACKd\nvkNw4kWEh49r9HU4e0TAmLplbAgwVV2EE8zi02bEvVfgslTLuulnNqYuxmSL7Pga/8KasdE1kEdo\n3DkEx50L7QqyWDNjcpcNAWZMPZp69yQAFWUEXngC/4uPI1WVseTy2x5Fu3RrWgWNMUmlGuiiQ4C9\n14x1MSYn+V96isDCh5FdJQCE+w/F9/6/ASzIGZMBGRsCzJg2JRyKvc2bfS8A1QccTtU5PyJy4JGN\nHufSrsMZ03AZGQLMmLZCvt6K/9Vn8b32fCytet8+BM+6hOr+xzd5TEpjTMNlaggwY1q0Js/qveJd\n/P+aj/d13GOsAAAfBUlEQVT9t5DdLmVX/PEB8GbqSR5jzO5S/t8nIp2Ag9jz8YIl6a6UMa1CsCr2\ntt0dvwJAvT5Cx44gdNIZsUcHLMgZk12pPkc3BWcGgVJqP16gUOfADcbkrkgE31uvEHjqwZqkLt0I\njTid8Imnop33atLq7TqcMemV6k/Nm4Hvq+oLzVkZY1o67/+WE3hiBt4v1tRKL//zY3bmZkwLler/\nTB/wUnNWxJjmlI5n4fJv/xW+le8CEOmyN8EzLyL/fnf+YAtyxrRYqQ7WfBvwOxFp6uDOxrQqEjfI\nsm/lu2hBe6rO/jHltz1C+ISxWayZMSZVqf4MvRLYB/i1iHwdn5HKWJfGtDZS8g3+hY/gf/XZWFrw\nlLMIfu8HUNipyeu363DGZE6qge4HzVoLY9IgLUN1le0i8MLj+Bc/iQQr0bjn3oLn1zkUqzGmhUv1\nObrXmrsixrQE7a86HynbBUB4wFCCZ15MwbUXZblWxpimSPXxghsT5dkQYKZVq6rE/9pzsY9Stovw\nwUcTPOsSIgcclsWKGWPSJdWuy92H/doHOBGYn97qGJMhFWX4lyzA/8ITeHbtqEn+1R1UH35sWobq\nsutwxrQMjR4CTETGAuelvUbGNKfSnQRefgr/y0/Huiir+x6M93NncvrqI47LZu2MMc2gKQ//vAQ8\nnq6KGJOKpt5w0v6X5yCVFQBUH3QUwdN/QPXhx1I4ZWSaamiMaWlSvUa3+zBfBcD5wJdpr5Ex6VS6\nk8ALNb/HpLKC8OHHETz9B0QOPqrpq5+1lOLiYoqavCZjTHNJ9YxuNc64ltELF+XA+8Dk5qiUMU1W\nXkpg8TznMYGKsprk6+4j0u+QLFbMGJNpqQY6v6pWN2tNjHE1tXvS/+yjBF54vOYxgcOPw7diOYAF\nOWPaoHqH9BIRL1AqInlN2ZCIjBWRT0VktYhcXUe+iMg9bv5HIjKgvrIicoeIfOIuP19EOjeljqYV\nC4djb/OefAAp20X1wUdR/tt7qLzqjixWzBiTbfWe0alqtYh8BuwFbGrMRtxgOR0YA2wAlovIQlX9\nOG6xcUCR+xoE3AcMqqfsy8A1qhoWkduAa4DfNKaOpvXyfriMvMemxz5XH3AYwUkXUX3oAHtMwBiT\nctflbOA5EfkrTrDRaEaKE68OBFar6loAEZkLTADiA90E4GFVVWCZiHQWkR5An0RlVTV+RoVlwPdT\n3B+TA2TjOvIe+3/4/vdOrfSK3/0tLQHOGJMbUg10l7p/b9gtPdWJV3tS+w7NDThnbfUt0zPFsgAX\nkeRxh+Li4normcoybVk626f/Hy8B4P3f3b9nXpLtRfMq7v0DXd97DdEI4bx2bBl2Gr1emeeUW726\nQetMJb8+9t1JzNomMWubxOprm6Ki1O91TvWB8b4przELRORaIIxz5lmn+hqluLi4QQ3X1jRX+yRb\n5x557s0lAN3efRUVD6GR36Nq0kV07tgZ3EDXoHW6ot2TjdlD++4kZm2TmLVNYulum0zNFrmR2sOI\n9XLTUlnGn6ysiEwBTgNGu92eJsfI9s34X3oK/+vPx9LChw4geP7lRHqn0qFgjGnLMhXolgNFItIX\nJ0idi/PAebyFwOXuNbhBQImqbhaR7YnKusOQ/Ro4UVXLM7MrJlM8a1bhf/EJfMudLsp4lb++067D\nGWNSkpFA594VeTmwGPACM1V1pYhMc/NnAIuA8TgPp5cDU5OVdVf9NyAPeFmcg94yVZ2WiX0y9WvU\n83BxJ+UFNzqXhtXrJTRoDKGxZ1NwnXNtz4KcMSZVmTqjQ1UX4QSz+LQZce8VqHNmy7rKuukHpLma\nJos8a1aR90TsK4EWtCc04nuExkxEu+zdpHXbIwLGtF0ZC3TGJCJbNxJ48gH877xaK73srnnQriBL\ntTLG5AoLdCarAo/ei3/JAqQ6jPoDhE7+PoHn5ziZFuSMMWlggc40SaOuw1XW3DcUePkpVITQCWMJ\nTroI3WvvmkDXQNY9aYypiwU6kznBKvxLFuJ/ruZxx/ARAwme/WMi+/XLYsWMMbnMAp2pV+HkEfSn\nCWdMoSC+1xcRWPgonh1f1cqq/NXtTa6fMcYkY4HONLuCqy/A89VWAKr3LyI46SLa3X1NlmtljGkr\nLNCZ9ItE8C1fGvvo+Wor1T37ODMKDDgBPPXODpWQXYczxjSUBTqTPqp4V7xLYN79eNd/FkuunPZ7\nwoNGgMebvboZY9osC3QmLTxrVhGY9w98q94HINK5a+x6XHjI6GxWzRjTxlmgM2kRG66roJDgaf9H\naMwkCi85pVHrsu5JY0w6WaAzDX8WLlKN96N38C9ZEEvSQB6hMWcSPPU8aN8h/ZU0xphGskBnUiYl\n3+B7fRH+pc/G7qKMKr99NvqdrlmqmTHGJGaBziQXP5vAlWcj1WEAIt32JTTydPKe+LuzmAU5Y0wL\nZYHOJCTbN5P3yF9rEiIRwgOGEho1gerDjgWPJxboGsquwxljMsUCXRvRoOtw4RD+F54gsPBhJFgV\nSy6/cy66V9OmyzHGmEyzQGdq8XzyIfmz7sKzaT0AoSEn4X/rFQALcsaYVskCnYnJu/9W/G++CECk\ney+qJv+c6sOOjQU6Y4xpjSzQ5ZBGTZlTtiv21v/mi6jfT/DU/yN06nkQyGtSfew6nDGmJbBA10bJ\nV1vwL56H/7XnY2nhw46l6sKfo/v0ymLNjDEmvRo/um4DichYEflURFaLyNV15IuI3OPmfyQiA+or\nKyJnichKEYmIyLGZ2pfWzLPuM/L+340UXHU+gZeeQqoqY3mVV91hQc4Yk3MyckYnIl5gOjAG2AAs\nF5GFqvpx3GLjgCL3NQi4DxhUT9kVwCSgcfe4txXxz8Jd/yMnyesldPwYQmPPpuC6S5xMkQav2ron\njTEtXaa6LgcCq1V1LYCIzAUmAPGBbgLwsKoqsExEOotID6BPorKquspNy9BuZF9Dr8N5Pv2IvKcf\njH3W/AJCI04jdPL37S5KY0ybkKlA1xP4Mu7zBpyztvqW6Zli2XoVFxenZZls6+/+rauu8XkFG9fS\nY+kCCj//uNYyH11+C5H8AvimxHnVs85U8o21TTLWNolZ2yRWX9sUFRWlvK42czNKfY1SXFzcoIZL\nRaPugkxRsroe8dxMfB+8BbhncGPPIvDMLAD6HXFUg9dZOmtps7RPrrC2SczaJjFrm8TS3TaZCnQb\ngd5xn3u5aaks40+hbJvnWV/z68f3wVtoIJ/QmEkEx58DhZ1igc4YY9qaTAW65UCRiPTFCVLnAufv\ntsxC4HL3GtwgoERVN4vI9hTKtlmeNR8TWPAwvg+XxdKCY88mdOp5aMfvZLFmxhjTMmQk0KlqWEQu\nBxYDXmCmqq4UkWlu/gxgETAeWA2UA1OTlQUQkYnAvUA34HkR+UBVGzfbZwuSSpen57P/OQFuxXIA\nNJCPBJ1HBYLn/aRB27M7J40xuSxj1+hUdRFOMItPmxH3XoHLUi3rps8H5qe3pq1Dwc1XAKD57QiN\nnkhw7NkUXnFGlmtljDEtT5u5GaW186wvJjDv/thnLWjvzOh98plQ2CmLNTPGmJbNAl0TNfbOylTL\nyZYNBJ5+EP/br9ZKL/vzXGjfIeXtWfekMaatskDXguU9dCe+155HIhHU7yc0eiKBF59wMhsQ5Iwx\npi2zQNfCyI6vY+/9rz6LiofQ8PEEz5iC7rV3TaAzxhiTEgt0LYR8tQX/orn4X6+55yZ87HCqzrwY\n3Xf/LNbMGGNaNwt0WSabvyDw3Bx8b72MVFfXyqu84sYGrcuuwxljzJ4s0GWA7PwWQkEIBZGg8zeq\n4JrJiKrTRTnkJEKnnU/BtRdlsbbGGJNbLNCloMF3VkYi+JY+G/vY/oqJiZf1eAmdMJbgqeei3W0u\nOGOMSTcLdGkmm9aT/88/4/3sf7G0SIfO4A+AP4C6f72ffwJA+Z/noF1Sny7HuieNMaZhLNClSziE\n//nHCCx8BAmHiHT6Dp6SbwEo/9szeywePUtsSJAzxphc0Pmfzrj8O6b2zMj2PBnZShvQ7vofk/f0\nTCQcIjR8POV/stkCjDEN1/mfG2OBIJe3mUl2RtcUpTtjb70b1hLpti9VU39J9WHHNG211j1pTNpl\n+iwiGacuBexo4JRr2diHltRujWWBrjHKdhFYPA//S0/FkoLjzyV4xhTIy89evYxpw5pyQM6Fg3l9\nmmMfk62zJbWpBbqGKNtFYPGT+F96Eqkoq5UVPGdaliplWqOWdBDIhtZygITG1yfT+9HS2i1eRVjZ\nWBZmQ1k1X5ZW118gzSzQNUD7X52LlDsBLnzoAIITp1Bw808btS7rnmxZWtJBoiXVpSka2z3XmjTk\n3ypYrZSGIuwKaSztpS8r2RWKsDOosb9Rv162g3yvkOcV92/NuuZ/Xk57n4f2fqG9Tyj0Syxvc3k1\n5SGlLByhPKxUhJWycM167/xwFyXBiPtSSoKRWN7QZ7aiCgpEFGpyYPiCbbF6ROsVNfXVbwirEo5A\nRJW4zdFvzma+ropfU+ZZoGsAKS9zAtwZU4gcdGS2q5M2Le3A2prq01oO5vW1aUtr80R2r2dpKMLW\n8ghbKmrOEv784S7KQhFKQ0ppWGPvo85++SvyvUK+zwkg+XEH7BveLUm47ds+2InfI/g9EPAIAU9N\nuTs/3MWOYIQdVU4A2REXsA54bDO7QhGq6jiROfuVr/dMdP1jVVnCvKlLv02Yd8jjWxLmAdz0350J\n81Z+G06Y99E3oYR589dVJMz7uiqC3wM923vp5b7mrkm8fHOwQNcA5df8lcjBR2W7Ggm1loNVc2mN\n+18WitDOJ3hE6l+Ypu9jaShCsFqprIZgRKmsVqqqaw7KH38bIs8jBLyQ53UO5n3mbAZg3fk9CKsS\nikAo4vx6H/DUVgAWju3KrmCEnSFlV7D2WctP3viWirBSHo5Q5p5hRPV/cguRuDMIjTsTOGreFgJe\nJ7BE6xJ1zFNb2FoeoTT+1MH1xyQHcoCXNlQlzPvL/0oT5t3y/q6EecmCx1eVztmMV6CDXyj0e9hQ\n5kS9k3rm0THgoYNf6OD30DEg/Mndzq2DOlFVXfNvVBFW7vvYCX4T+uRTFnLO1Mrcs7c1O511dm/n\nocAntPM5Z3sFPufzC19WAnDlEYV0CnjoFPDQOU/oFPAw6SUn4L4xYW88gEdAxLktf+D8bQAsPb1b\nrC6V1VBZrUx+9RsAHjzxO/g8glfA5wGfCN9/2VnnJ+fsw97tPLW+43PXZPYOTwt0DZDOINeaDsrJ\nzloae60lm/uvqpQElU3l1Wwur2ZTWc1P7b98tAuvB7zi/KeN+7HPg5+UOt05ume3zuzislhXUoF7\ngGnvF76pjLCutJrPd4b5fFeYdbuq+XxXza/mno86QaSdV2IHp6hhC7YhOAec6N+oU1/YHus+auf+\nnbO6HICz+7VzuqSqIrEuqqhe7vYSOf6ZbQnzogGvLt978auEedF61eXzXYmv16xPci0nelDP90L3\ndl72KfDy9jZnaL1fHFlIe5+HQr/zb9DB76G9r+bAO/ekLlSGnQN1pRtArnnHOZO7/piOe2zrD+85\nQeyqozoQijiBPhhRQtXKQ585+xYNHp3zPHR2A8gZi53tfXbuPnTwe8j3grj/iNHv/5Mnd91je9FA\nN+3Qwj3yooFu1si99siLrvPTc3vU2WbR/OuPTTxR8xFd/Anzju4aSJh35ncLEubtU+BNmJcpFuha\nmXQGiOgBP2rJxkr8HiHgwf0l7byPWr8rHDvYivuK+mxHiGDEuQ5RFVGCcceop9eWO3kR59dg/BnE\nte+UxH7pl4eUiri8M1/6yjlYuQGkg1+46yPnF/evj+5AWcg9S3C7qKLGLdoOuL9KAY8Ir292fsX3\n6+hlc7lz7aIuN7yX+Jf5L99K3K112Zs7EuYl084rVFRr7EXcycb/knQV/XtLMGHeE0m6hQp8zr9p\nNEhGX9FtHdzZR1W18+9XWa0EIxo7O+vgF3wenO47cd5Hg9HQfQJ09HvoEBDnr1+42z07undoZ9q7\nQTx6djHmeeff6L9ndo8F8ei/1xHznLPED77f3alLBEJuXcYucgLqsol7s087L50CskfwuO6YxAdy\ngLG92+2RFg10Vx655zyP0UB37YA9g2A00CULHnu3y/6BPhta0o/4jAU6ERkL/BXwAg+o6q275Yub\nPx4oB6ao6n+TlRWRLsDjQB9gHXC2qibuvG6kp7sex6x9huN/9RvyvdDOV/NLusP+k4iIUPZuSexg\nHqpWQgp6yKX4NILvPzvI9zkHtei1gagF6ypiZwAF7llAgc+DAGVxB/HoxeWoW9/fyU734vVOt8so\nauzz252ukIDzq7aDe+DxeYRNZdV8WVpz91N8gIh2XyRy1JNbE+ZFuzfqctFrif9Jpq9M3FX0r42J\nu5hu/yBxN9JbWxMHgeiZQKFP6NHeS48CLz0KPDzuBoefHl5ItUI4okQUqhVmfur8ir7ooPa1unQ8\nUvML+9x+7SgPx3claSx4TOzTjr4dvfTp4KNvBx99O3g53D2Yb75wXyLqnFVUVDtlo+289PRugNOt\nF+3SG/WcEyAWju1aq1sr+jccgY4BD50CEuue6hSQWPDYdMG+dbZLNEgsm9g9Yd6XP6i7bCLRQHfB\nge0TLvPdjokPQX067JnXkg6e0LLqk426tKT9TyYjgU5EvMB0YAywAVguIgtV9eO4xcYBRe5rEHAf\nMKieslcD/1LVW0Xkavfzb9Jd/08LevB81wFQ1wXXvmc6f+vq2+9+gruCxBeVo33cDXVrkgP9sm2J\nD/S7K/RJLNiN3DfPDdROwA663TTFJU5XW6/2Nb9MnTuzlE3lTvA9oKMvdiaY515XedM965jQJ9+9\n7iPutZaaAHHTsR3dX/pCe7/za/8st4tp3pi9KHVvJIgGkOi1kN8c3YHCuDLt/cIFS5y2XDSuKxG3\njhEFUCa43UhvT9ybHgVeOgZqDwr0uHvN4Mbj9vxlHg10dx3feY+86H7MGN4l5TbfnUec+rf3Q9e4\nxzCTdRUN75HX6O01VGs5mDWXxu5/c7RbU9aZrGxj85qjLs0hU2d0A4HVqroWQETmAhOA+EA3AXhY\nVRVYJiKdRaQHztlaorITgBFu+VnAUpoh0J21/W0OK9vIjstvosL99VxZDZVhJfL0Q3hV4ftT8Qv4\n3YvmPg/89N9Od9YdgztR6f5ir9ztovLp++e7XXDOwbw87qJyv47eWLdd9KD+jBtsf3VUBzr5xf31\n7lzEjp6RLRrXlV2hmluZozcJhCJKjwIvvQudO5/2K/TRKSB856FNAMw/Zc/rBVDzi37F2fskzHv3\nzMRnAnVdT4ju/xVH7NlVFDWm154P3//yqMTLRx2/z55BIFcO1s110G2ug1lxcXGDt5eNwJLp70dj\n28Y0TqYCXU/gy7jPG3DO2upbpmc9ZburavQK+RZgz6OtK9GXKpVl+lds5YCKrbxf7d4p5HVfAei/\nzhkd5f124+oo6VygHeHbvEdL3+fmXde7YWd0z6xzyp3TIa4bMey+3HV22/UF3aJ5HiDffUVVOa/t\n38D2uHombqNk+YnzlrsntA3Nq78+dUu+zvo0x36kf3utUUvZj5bYri2pLi1NfW1TVJT6Mz05czOK\nqqqI1H2HAfU3SnFxcb3LJMuvM+/NjY3LSyLp81qNXGd95XYUJW6f5nh+LBvPpDVlm6l8d9K5vdak\nMW3TVljbJJbutslUoNsI9I773MtNS2UZf5KyW0Wkh6pudrs5E98R0Ux8I2YD0NB77lpSv7d1lRhj\nclmmpulZDhSJSF8RCQDnAgt3W2YhcKE4BgMlbrdksrILgcnu+8nAguaovG/E7FhAM8YY07pk5IxO\nVcMicjmwGOfq1kxVXSki09z8GcAinEcLVuM8XjA1WVl31bcCT4jIxcB64OzmqH9LulvJGGNMw2Ts\nGp2qLsIJZvFpM+LeK3BZqmXd9K+B0emtqTHGmFxiM4wbY4zJaRbojDHG5DQLdMYYY3KaBTpjjDE5\nzQKdMcaYnGaBzhhjTE6zQGeMMSaniWrC4SFbvZKSktzdOWOMMQB06tRJkuXbGZ0xxpicZoHOGGNM\nTsvprktjjDHGzuiMMcbkNAt0gIiMFZFPRWS1iFyd7fpkgojMFJFtIrIiLq2LiLwsIsXu3+/E5V3j\nts+nInJKXPoxIvI/N+8eEUl6Ubg1EJHeIvKqiHwsIitF5GdueptvHxHJF5F3RORDt23+4Ka3+baJ\nEhGviLwvIs+5n61tABFZ5+7TByLyrpuWmbZR1Tb9wpn6Zw3wXSAAfAgcmu16ZWC/hwMDgBVxabcD\nV7vvrwZuc98f6rZLHtDXbS+vm/cOMBgQ4AVgXLb3LQ1t0wMY4L7vAHzmtkGbbx93Pwrd937gbXf/\n2nzbxLXRL4A5wHPuZ2sbZ5/WAV13S8tI29gZHQwEVqvqWlUNAnOBCVmuU7NT1deBb3ZLngDMct/P\nAs6IS5+rqlWq+jnOnIED3VndO6rqMnW+gQ/HlWm1VHWzqv7Xfb8LWAX0xNoHdZS6H/3uS7G2AUBE\negGnAg/EJVvbJJaRtrFA5xzAvoz7vMFNa4u6qzOrO8AWoLv7PlEb9XTf756eM0SkD9Af58zF2odY\n19wHwDbgZVW1tqnxF+DXQCQuzdrGocArIvKeiPzITctI22Rs4lXTuqiqikibviVXRAqBp4Cfq+rO\n+EsBbbl9VLUaOFpEOgPzReTw3fLbZNuIyGnANlV9T0RG1LVMW20b1wmqulFE9gZeFpFP4jObs23s\njA42Ar3jPvdy09qirW7XAO7fbW56ojba6L7fPb3VExE/TpCbrapPu8nWPnFUdQfwKjAWaxuAocD3\nRGQdziWQUSLyKNY2AKjqRvfvNmA+zmWjjLSNBTpYDhSJSF8RCQDnAguzXKdsWQhMdt9PBhbEpZ8r\nInki0hcoAt5xuxx2ishg986nC+PKtFruvjwIrFLVu+Ky2nz7iEg390wOEWkHjAE+wdoGVb1GVXup\nah+c48gSVf0B1jaISHsR6RB9D5wMrCBTbZPtO3FawgsYj3Nn3Rrg2mzXJ0P7/BiwGQjh9HNfDOwF\n/AsoBl4BusQtf63bPp8Sd5cTcKz7hV0D/A13EILW/AJOwLme8BHwgfsab+2jAEcC77ttswK4zk1v\n822zWzuNoOauyzbfNjh3tX/ovlZGj7OZahsbGcUYY0xOs65LY4wxOc0CnTHGmJxmgc4YY0xOs0Bn\njDEmp1mgM8YYk9Ms0BnTSCLyWxF5oP4l07rNP4rIVyKyJcXlVUQOaO56GdOS2eMFxrQSIrIfzjNF\n+6szukQqZRQoUtXV9SzXB/gc8KtquIlVzfp2jIlnZ3TGtB77AV+nGuSMMQ4LdMbUQ0R+IyIbRWSX\nOwnkaDf9BncsQ0TkbyJSGvcKi8gNbt6+IvKUiGwXkc9F5KdJttVJRB52l10vIr8TEY+InAS8DOzr\nrv+hBOWvEpHNIrJJRC7aLe9UcSYE3SkiX0br53rd/bvDXf8QEeknIktE5Gu3u3R2dPivetrFIyJX\ni8gat+wTItIl0XZS+kcwpgks0BmThIgcBFwOHKeqHYBTcCaQrEVVL1fVQlUtxBlC7FtggYh4gGdx\nhj7qCYwGfi5xMybv5l6gE86QSSfijOU3VVVfAcYBm9ztTKmjrmOBX+GMP1kEnLTbImXu+jrjzJl2\nqYhE5/Ia7v7t7K7/LZyJLW8B9gUOwRlk94YU2uUKnDnCTnTLfgtMT7IdY5qVBTpjkqvGmeX4UBHx\nq+o6VV2TaGER6QY8A1yhqu8DxwHdVPVGVQ2q6lrgfpxBf3cv63XTr1HVXaq6DrgTuCDFup4N/FNV\nV6hqGW5QilLVpar6P1WNqOpHOOOdnphoZaq6WlVfVmfyy+3AXXHLJ2uXaThjGW5Q1Sq3Ht8XEZsW\nzGSFBTpjknBv4vg5zsF6m4jMFZF961pWnKl9ngTmqOpcN3l/nO7GHdEX8FtqJpiM1xVnxu71cWnr\nSX3SzX2pPVll/HoQkUEi8qrbLVqCE5C6JlqZiHR393ejiOwEHo0uX0+77I8zT110f1fhBMa69tmY\nZmeBzph6qOocVT0B5wCuwG0JFr0X2An8Li7tS+BzVe0c9+qgquPrKP8VzmwS+8el7Ufqc5FtpvYc\nXvvtlj8HZ/qT3qraCZiB0z0Jzn7t7k9u+hGq2hH4QdzyydrlS5zR5uP3OV+d+cjsNm+TcRbojElC\nRA4SkVEikgdUAhVApI7lfozTrfd/qhqf/w6wy71xo52IeEXkcBE5bvd1qDNz9xPAzSLSQUT2B36B\ncyaViieAKSJyqIgUANfvlt8B+EZVK0VkIHB+XN52d7++u9vypUCJiPQErorb32TtMsPdh/3dZbuJ\nyIQk2zGmWVmgMya5POBWnLOtLcDewDV1LHcezsF7U9ydl791g9dpwNE4z499BTyAc8NJXa7AuWlk\nLfAmzlnYzFQqqqovAH8BlgCr3b/xfgLcKCK7gOtwAmO0bDlwM/Bvt8txMPAHYABQAjwPPB23rmTt\n8lecM8eX3G0tAwYl2Y4xzcoeGDfGGJPT7IzOGGNMTrNAZ4wxJqdZoDPGGJPTLNAZY4zJaRbojDHG\n5DQLdMYYY3KaBTpjjDE5zQKdMcaYnGaBzhhjTE77/+77LhZUCGTlAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb64570bdd8>"
      ]
     },
     "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": 39,
   "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": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbMAAAEGCAYAAADv6ntBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XecFdX5+PHPc8vuAksRUKQJaCA2FBVRWOnYQEUwNhJF\no0ls8asxlhi7MSYmv+RrDF+JsQQVoiaC3ShEEUVULIggwtKLdKQsW257fn/M3Lt3++zuvXfb8369\nljt35sycM+fe3Yczc+YcUVWMMcaYpszX0AUwxhhj6suCmTHGmCbPgpkxxpgmz4KZMcaYJs+CmTHG\nmCYv0NAFSIU9e/ZYl0xjjGnm2rdvL1Vts5aZMcaYJs+CmTHGmCavxQWz/Pz8hi5Co2V1UzWrm+pZ\n/VTN6qZqqaybFhfMjDHGND8WzIwxxjR5zaI3ozGmZVBVCgoKiMViDV0Uz3JyctizZ09DF6NRKl83\nPp+P3NxcRKrstFglC2bGmCajoKCA7OxssrKyGroonmVnZ5OTk9PQxWiUytdNKBSioKCAtm3b1vpY\ndpnRGNNkxGKxJhXITO1kZWXVudVtwcwY06x1eGoTHZ7a1NDFMGlmwcwYY2ph6tSpDBo0iJ/85Ce1\n2m/37t08/vjjaSpV7bzxxhv8+c9/buhipJTdMzPGmFp44okneOmll+jevXuNaX1rltMaiPX5Pnv2\n7OGJJ57gyiuvTH8hqxGJRBg7dixjx45t0HKkmrXMjDHGoxtvvJG1a9dy/vnnM2XKFD777DNOPfVU\nhg4dymmnnZZ4CHjZsmWMGjWKvAt/xOALfsiqVau49957WbNmDaeccgp33nlnhWP/85//ZMiQIeTl\n5fHTn/4UgHXr1nH22WczZMgQzjnnHDZs2MCePXs4+uijE/eW9u/fz1FHHUU4HGbatGmMHDmSvLw8\nLrnkEgoLCwG4+uqrufHGGxk9ejR33XUX06dP5+abbwbgzTffZPTo0QwdOpTx48ezbds2AB588EGu\nvfZaxo0bx7HHHsvUqVOrLeuOHTu45JJLGDlyJCNHjuSjjz5K06dQOWuZGWOarNrcC/OSdvfl1be2\n/vznPzNnzhxeffVVOnXqxN69e3nzzTcJBALMnTuX++67j2eeeYYnn3ySq666iotOPJZQOEy4Wzfu\nvvtuli1bxgcffFDhuMuWLeOPf/wjb7/9Np06deK7774D4JZbbuHiiy9m0qRJPPPMM9x6663MmDGD\n/v3788EHHzBs2DDeeustRo0aRTAY5Oyzz2by5MkA/OY3v+GZZ57hZz/7GQDffvstb7/9Nn6/n+nT\npyfyHjx4MHPmzEFEePrpp3n44Yd54IEHAGeEjldffZWCggIGDhzIFVdcwcqVKyst62233cY111zD\n4MGD2bBhA+eddx6ffPKJh08mNSyYGWNMHe3du5err76a1atXIyKEw2EABg0axP/7f/+PzaeN5uxR\nIzi039HVHmfevHmce+65dOrUCYADDjgAgIULF/Lss88CcNFFF3H33XcDMHHiRGbNmsWwYcN48cUX\nE5cuv/76ax544AH27NlDQUEBo0ePTuQxfvx4/H5/hbw3bdrE5ZdfztatWwmFQvTq1Sux7bTTTiM7\nO5vs7GwOPPBAtm3bVmVZ586dyzfffJPYd9++fRQUFJCbm+u9QuvBgpkxpsmqqSUFpS0yL2lr64EH\nHmDo0KFMnz6ddevWcdZZZwFw/vnnM3DgQN6e8Qw/+PmN/PmvU+jdu3fK8j3zzDO57777+O677/jy\nyy8ZNmwYANdccw3Tp0+nf//+TJ8+vUwrsE2bNpUe65ZbbuHaa69l7NixvP/++/zud79LbMvOzk4s\n+/1+IpFIlWWKxWLMmTOnwZ6ps3tmxhhTR3v37qVr164AzJgxI7F+7dq19O7dm6snXci4EcNYunQp\nbdu2Zd++fZUeZ9iwYbz00kvs2rULIHHpbtCgQbz44osAvPDCCwwePBiA3Nxcjj/+eG699VZOP/30\nRIuroKCAgw8+mHA4zL/+9S/P59CtWzfAuRdWk6rKOmrUKB577LFEusWLF3vKP1UsmBljTB39z//8\nD/fddx9Dhw4lGo0m1s+aNYvBgweTd+GP+HrVai6++GI6duzIySefzODBgyt0ADniiCO46aabGDdu\nHHl5edx+++0APPTQQ0yfPp0hQ4bw/PPPl2k1TZw4kRdeeIEJEyYk1v36179m9OjRnH766fTt29fT\nOdx2221MnjyZ4cOHJy4dVqeqsv7+97/niy++YMiQIZx00kk89dRTnvJPFVFt+pM012am6fz8fM8f\ncktjdVM1q5vqZap+9uzZQ/v27Wu1TzovM9bEt2Y54HTNNxUVFxdXuCxZ3Wdc3UzTds/MGNOsNUQQ\nM5lnlxmNMcY0eRkLZiJyhogsF5GVInJbJdt/KCKLReQrEflQRI5N2rbWXb9IRD7NVJmNMcY0DRm5\nzCgifmAKcCqwEVgoIq+o6tdJydYAw1X1OxE5E3gMOClp+0hV3ZGJ8hpjjGlaMtUyGwSsVNXVqhoC\nngPGJydQ1Q9V9Tv37UdAjwyVzRhjTBOXqWDWHdiQ9H6ju64qVwBvJr1XYI6IfCYiP01D+YwxzVTu\n5BHkTh7R0MUwadboejOKyEicYHZK0upTVHWTiBwEzBaRb1R1XmX7xwf6rI6XNC2V1U3VrG6ql4n6\nycnJKTMqhRfxwZSKi4tTX6AatAae+NdM/Ad354ILLsh4/vUxadIkHn300Vo/ClFb5T+XvXv3JgY7\nBjw/8pGpYLYJ6Jn0voe7rgwROQZ4HDhTVXfG16vqJvd1m4jMwrlsWWkwq+nE7XmhqlndVM3qpnqZ\nfM6srsMlNdQwS1ecP7FJPWemqqgqM2fOrPexanrOrrLnzNq1a0fPnj0rTV9tXrUvXp0sBPqKSB8R\nyQIuAl5JTiAihwAzgUtUdUXS+jYi0ja+DJwGLMlQuY0xJmHdunUMGjSI66+/npNPPpkJEyZQVFQE\nUOX0K7+d+nceeeQRVqxYwahRo8oca8iQIQAsWrSIsWPHMnz4cCZOnMiWLVsq5L1t2zZ++MMfkpeX\nR15eHh9//DEAf/3rXxk8eDCDBw/m//7v/wC45557+Pvf/57Y98EHH+SRRx6hoKCAc845h2HDhjFk\nyBBef/31RFkGDhzIz372MwYPHszGjRvp378/O3c6bYpJkyYxfPhwTj75ZP7xj38kjtu9e3fuv/9+\n8vLyGDNmTKJFFS/rkAucn3hZn3/+eUaNGsUpp5zCDTfcUGbUlPrKSMtMVSMich3wFuAHnlTVpSJy\nlbt9KnAX0An4PxEBiKjqQKALMMtdFwBmqOp/MlFuY0zjVpt7YV7SFkybW2OaVatW8fjjj/OXv/yF\nyy67jFdeeYULL7yw0ulXrj6tNM9+/foRCoUS4zbOmjWLCRMmEA6HueWWW5gxYwadO3dm5syZ3H//\n/UyZMqVMvrfeeit5eXlMnz6daDRKQUEBixYtYsaMGcyZMwdVZcyYMeTl5TFhwgR+9atfJWbDfuml\nl3jxxRfJycnh2WefpV27duzcuZMxY8YkJulctWoVjz76KCeeeGKFc54yZQoHHHAARUVFjBo1inPO\nOYeOHTuyf/9+Bg4cyJ133sldd93FtGnTuPnmmxNl/edv7iIajbL3wO4sX76cmTNn8tZbbxEMBrnp\nppt48cUXufTSS2uscy8yds9MVd8A3ii3bmrS8pVAhSlYVXU1cGz59cYY0xB69erFMcccA8CAAQNY\nv349UMX0K0nBDGDChAnMmjWLG2+8kZkzZ/LUU0+Rn5/PsmXLOPfccwFn9PkuXbpUyHfevHmJCTL9\nfj/t27dnwYIFjBs3LjEi/llnncWCBQu46qqr2LFjB5s3b2bHjh106NCBHj16EA6Huf/++5k/fz4+\nn4/NmzcnWlM9e/asNJABTJ06lddeew1wpoxZtWoVHTt2JCsrizPOOCNRF++++27Zsn67NlHW559/\nni+//JKRI0cCziXGDh061O1DqESj6wBijDFeeWlJxVtkXtJ6UX5alPhlxuqmX4mbOHEikydP5uyz\nz0ZEOOyww1i6dCmHH344s2fPTkn54saPH8/LL7/Mtm3bEoMRv/DCC+zYsYP33nuPYDBI//79Ex0w\nqpoi5v333+e9995j9uzZtG7dmnHjxiX2CQaDuFfNapwiRlW5+OKLE3OyQWo75dhwVsYYkwJepl/p\n06cPfr+fhx56iIkTJwJOp7UdO3YkZmUOh8MsW7aswr7Dhw/niSeeACAajbJnzx4GDx7M66+/TmFh\nIfv37+f1119PTBMzceJEZs6cycsvv5xo9e3du5fOnTsTDAaZN28eGzZsqJBPeXv37qV9+/a0bt2a\nFStW8OmnNQ/CVFlZhw8fzssvv8z27dsBZ+oYL/l7ZcHMGGNSwOv0K/GpW+IBJisri2nTpnH33XeT\nl5fH0KFDE4Et2e9+9zvef/99hgwZwvDhw1m+fDkDBgxg0qRJjB49mjFjxnDJJZdw7LHOXZkjjjiC\ngoICunbtysEHHwzABRdcwKJFixgyZAjPPfcc/fr1q/G8xowZQzQaZdCgQdxzzz0MHDiwxn3iZT35\n/EkMmzSZ5cuXc/jhh3PHHXcwYcIEhgwZwrnnnlumC3592RQwJsHqpmpWN9VrzFPApPoyY21kegqY\nxjblTF265tsUMMYYU4mGCGIm8+wyozHGmCbPgpkxxpgmz4KZMcaYJs+CmTGmyfD5fIRCoYYuhkmT\nUCiEz1e3sGQdQIwxTUZubi4FBQWJB5Ubu8A3iwGIdDy4WeZXk5rKs3fvXtq1a5d47/P5yM3NrTRt\njXnVaS9jjGkAIkLbtm0buhie5T7njK9YcOb5zTK/mtRUnm3bttVphPzK2GVGY4wxTZ4FM2OMMU2e\nBTNjjDFNngUzY4wxTZ4FM2OMMU2eBTNjjDFNngUzY4wxTZ4FM2OMMU2eBTNjjDFNngUzY4xJh1is\noUvQolgwM8aYNPAv/7Khi9CiWDAzxpg0kK2bSt9Ewg1XkBbCgpkxxqSBb8eW0uVNaxuuIC2EBTNj\njEkDSQ5mG1Y1YElaBgtmxhiTBmVaZustmKWbBTNjjEkDa5lllgUzY4xJtXAI2b0z8da/fiWoNmCB\nmr+MBTMROUNElovIShG5rZLtPxSRxSLylYh8KCLHet3XGGMaE9m1DUkKXlKwF/luu+f9cyePIHfy\niDSUrPnKSDATET8wBTgTOBK4WESOLJdsDTBcVfsD9wOP1WJfY4xpNHzbt1Rcl6L7ZhboKpepltkg\nYKWqrlbVEPAcMD45gap+qKrfuW8/Anp43dcYYxqT5PtlcXbfLL0yFcy6AxuS3m9011XlCuDNOu5r\njDENyldZMLMejWkVaOgClCciI3GC2Sl12T8/Pz8laVoqq5uqWd1Uz+qnVK/VK+hYbl101TLPdXSc\n+1pZ+rpuawheylNTWfv27espr0wFs01Az6T3Pdx1ZYjIMcDjwJmqurM2+8bVdOL5+fmeK6elsbqp\nmtVN9ax+ympVUljmvfr9ZO/aSt9DekB2K8/Hqa5O67qtIVRVnlR+bzJ1mXEh0FdE+ohIFnAR8Epy\nAhE5BJgJXKKqK2qzrzHG1CSTHSfK3zOLde2FqOLbuCYj+TcG/sUfZzS/jAQzVY0A1wFvAcuAF1R1\nqYhcJSJXucnuAjoB/ycii0Tk0+r2zUS5jTGm1sIhfLt3oL7SP6+xQw4DWtZ9s8CHszObX6YyUtU3\ngDfKrZuatHwlcKXXfY0xpjGSndsA0I4HJVposZ6HAbNbTo/G4kICn32Q0SxtBBBjjEkh3/bNAMQ6\nH5xYFzvke4A7EkgLEPh8PhIqzmieFsyMMSaF4q0xPbBrYl3iMuOGVS1iBurAgjkZz9OCmTHGpFD8\nGbPklpm2O4BYh05IcRHittyaK9n7Hf4lC1G/P6P5WjAzxpgUSrTMkoIZxO+bNf+RQAIfv4vEYkT7\nD8povhbMjDEmhSprmUHppUZ/M+/RGO/FGBk8JqP5WjAzxpgUKr1nVr5l5nQC8W1wOoGk/Lm3ktIO\nF/6ln6buuLUgWzbiX70MzWlF5Li8jObtKZiJSLaIPCAiq0Vkj7vuNBG5Lr3FM8aUZ6OmN2KhEny7\nd6J+P9qhU5lN0V5uMEtTj8asmU8mloPvNMy4EsEFbqvshGGQnZPRvL22zP4MHA38EIhP0rMUuDod\nhTLGmKZIdm4FnGfM8Jd9jFe7dEeDWfh2bIX9+1Kar2/lUoJv/Svx3v/FfGTPrpTmUSPVRC/GyJBT\nM5s33oPZBGCSqi4AYgCqugkbvd4YYxKqul8GgD9ArEcfJ92G1anLNFRCzuO/LzsZaDRK4IP/pC4P\nD3yrv8G3dROx9h2JHnlczTukOn+P6UKUGy1ERA4Edlae3BjTEOwSZMOq7BmzZPEejf4U9mjMevlp\nfJvXE+vas8z64HuvZ/SZtkD8EuNJo8CX2W754D2Y/QuYJiJ9AESkK/BXnIkyjTHGUDrDdKUtM0pH\nAqn2vlkknFiUTWurz2/NcoJv/BMVofiKW0vz6Xggvq2b8H+zKLEurf/RiUQIfPSOs9gAlxjBezC7\nHVgDfAV0APKBb4H70lQuY4zJmFT9oa/qGbO4qIdgFnzr34nl1vdcVfWAvZEw2Y//HonFCJ/2A2J9\njy7dNHQsAIG5r9Wq/HXlX/oZvn27iXXtSax3v4zkWZ6ngYZVNQTcCNzoXl7coZp0gdYY02LlTh7B\ncUDBtLkNXZQGV+09MyDW81An3abKp4KRXdvIenla6ftQMTl/e4DwisWUTLoOsrIT27JefRb/xtXE\nDupG6LwryhwnPHwswVeeJvDZ+5Ts2w1tO9TrvKqT/J+A8OBTQSRteVXH83NmItLanTyzLzBYRIaI\nyJD0Fc0YY5qWmlpmtM4l1vlgJByudHPWc48iSc+LFV92ExoMEnz3VVr95jpka+m8xMFXnwWg5Mc3\nV+gGr526EO0/CImECc5/uz6nVCuZflA6mdfnzC4FtgDvAM8n/dg9M2OMAecZsz27nGfMDuhUZbL4\nSCDl+Zd9QfDjd9Gk1ldk5NkU3TGF2IHd8K/Lp/U9P01sk2iU0OhziR5Rec/B8PCzAAjOfQ0ycCEt\n+r2j0IO6pT2fqnhtmT0EnKeqnVW1Z9LPIeksnDHGNBWJVlnHLtX25ouPBFJGJELW0w8DEDrrh2XT\n9+5H4b1/I3L8KUjh/tL1nbsQOv+nVCU6YDCx9gfg27weX/5XtTmVCrzcUww3UMePuNp0zZ+bxnIY\nY0yTlrhfdmAVlxhd8U4gyYKzX8T/7VpiB3UjfOaFFXdq05bi6++n5OJrEqtKLr8ZWrWuOqNAINER\nJPhuGjqCRCP4Py+dgDMyaETq86gFr8HsTuBPItI5nYUxxpimqsb7Za7ylxnlux1kvfQPAEp+dH2Z\nTh5lEwrhMy5IvI0ePbDGMoWHub0aF86tNl1tenPKzq1kzXyK1jddRKuH7yjdkMZOJl546s0IrMDp\nhn+NlPZUEUBVNfNPxxljTCNTU0/GOO18MJrTGikuBCDr+alIcRGR4/KIHntySsukXboTOfJ4Al9/\nnpLj5fzpNvyLP06MNhLr0gPf1o0pOXZ9eQ1mzwBP43T6KEpfcYwxpmny2jLD5yPW81D8+UsACC6Y\ngwazKPlhesZtj4w4u/bBTBXZvB7/8sX4VyxOrA58+RHqDxAeOIzIyLOJHj6A3MtGprjEdeM1mHUC\n7rJny4wxpnJe75mBc98sHswAQuMmVTkEVn1Fjs9D27ZH9u2pMW3wzefxr1iMf8VXSMHeCttLLrqa\ncN7p0K5hLylWxmswewq4BKd1ZowxphzPLTNKh7UCiB3YlfC4i9NWLoJZhPNOJ+s/L5Sui4TxrVmO\n/5sv8S8vHfIq+7lHS8vVoRPRfscQ69ef7Gf/AlB555RGwmswGwRcJyK/BrYmb1DVYSkvlTHGNCWh\nEnx7vkP9gQrzmFUmPuAwUGFkj3QIjzgrEcxyHvol/vwlSKi4YrphY4l+/xii/Y5xWopuH4l4MKtM\nYxn5xWsw+7v7Y4wxppxEq6zTQZ5GjI/16ptYjh6X/oGUtGvpI8EBdxbqWNdDiB4+gOjhx5Lz6P0A\nlFxxS9rLki5ex2acVnMqY4xpmWoaLb+CQNKf3gyPZVh8zd1EDz8Wbd+xdKUbzJqyKoOZiFyiqs+4\nyz+uKp2qPlnVNmOMaQmqmsesYNpc8vPz6VvZTg0kclLtex82lkuJ1amuZXYxTpd8cDp/VEYBC2bG\npFj8AdbG8EekMZWlsfL6jFlj1Rw+2yqDmaqOTVpuHA8SGGNMI1SbnowmPbyOmv9FFes/9ZqRiJwh\nIstFZKWI3FbJ9sNFZIGIlIjIL8ttWysiX4nIotrkaYwxmdDUW2bNgdfejBVGxhRnXKtDvewsIn5g\nCnAqsBFYKCKvqOrXScl2AdcD51ZxmJGqusNjeY0xpnKRMASCKT2ktcwaXrXBTETiD0lnJS3H9QaW\nesxnELBSVVe7x30OGA8kgpmqbgO2icg4j8c0xiRLHqAnFgOf57l3W5Scv9xJ8bX3VJjQss5KivDt\n9f6MmUmPmlpmq6pYVmA+8C+P+XQHNiS93wic5HHfeH5zRCQK/E1VH6sqYX5+fo0H85KmpbK6qVom\n6yY+3WJleVa17aAFb5HrLm+c/y5FB6dmusHqyuJle2MQ2Leb/vHlLz+C+3/OqguvI5bjTKFSn3PI\n2f4tRwAl7TqSv2pVpWlq8zkCcMff4xtrtV9dt6WLlzxrKk/fvt76glYbzFT1XgAR+UhV3/J0xPQ4\nRVU3ichBwGwR+UZV51WWsKYTz8/P91w5LY3VTdUaqm6qyzN5m3/hXFr999+J9721hEiKy1vT+Tfm\n705w9szEcuyAzuRuyOfof02h6JcPlRlnsC7n4C/cCUCgW89K96/pu1PXeqvrMRvT9ziVv1eerkOk\nIJBtAnomve/hrvNEVTe5r9uAWTiXLY0xgG/V1+T87bdl161d0UClaZwCC99LLBf9+hFiB3XDv24F\nrX97PbJzW72OXXq/LD0DBRtvMnVRfSHQV0T6iEgWcBHwipcdRaSNiLSNLwOnAUuq38uYlkG2bybn\nf3+NhEOEh5+VWO9fV7tLSbWZnLGpkd078SVNY6IHdqXo148Q7XEovs3rafXAz+t1/NqMlm/SJyPB\nTFUjwHXAW8Ay4AVVXSoiV4nIVQAicrCIbAR+AdwhIhtFpB3QBfhARL4EPgFeV9X/ZKLcxjRq+/fR\n6k+34dv7HZGjBlJy6Q2JTb71KyEWbcDCNR6Bhe8lJpOM0w6dKLr9YaLfOwrfzq1V7FmqumCf9cZz\nzjGtJ2OD8to1v95U9Q3gjXLrpiYtb8G5/FjeXuDY9JbOmKYn56934/t2HdHuvSm+7p4y4/1JqATZ\nvAHt3rvBytdYJF9iLKNNW4pu+SM5D9+ZGHzXv/RTokcNrFM+9oxZRZkcWcRTMKtmbMYSnJ6JH6lq\nScpKZYypUeDrz4m1P4DiX/wOWudW2O5fu4JICw9m8UuMGgwi4XDFBNmtKL7xt+ReeRoAOX+4hdB5\nPyY8blKtH22wllnD8toyuxQYjDOX2UacFlQX4FOc580QkfGqaqNzGJMhmpVN8Q0PVvlH1Ld2BeSd\nluFSNS6BT+chqkSOHkTgi/mVJwpmJRZFY2T/+3H8q5ZR/JPboE3b6jMoLkwslhmFvhFqDuMvVsfr\nfz2WAjer6iGqOkRVDwFuAr7ACWyPAo+kqYzGNGqZ7Dwh2zcnlot/dgexQw+vMm1tO4E0R4FP5gIQ\nOXG4p/RFNz6Its4l8MV8Wt/zM3zrK39uLM63I+l+mz2k3qC8tswmAeUfbX8U2KGq14nIH4CbU1oy\nY0wFwbdfTCxHBw6tNq1vXX6LHgkk+RJj5Pg8T/tEBwym8N7HyHnkLvzrV9Lq/msoueymsolCJfjW\nLMe/4isCX32ShpJXr7m3sOrKazDbCpwNvJy0bhwQf0AjB6jkgrQxJmX27yM473VPSWMdOuPbvQPZ\n9i16cGX9qpq/5EuMtGrjeT89qBtFd04he9qfCH7wFjmPlT7D1+r+a/GtWY5EI/UunwWl1PL6X7br\ngadFZL6IPCci83HmOos/oHESdpnRmLQKvvc6UlzkKW2stzOqgn9dy314OrBwLuD9EmMZWdmUXHkb\nxZf9Ak0alNi/cinEokR7HkZ41HiKf/brFJXW1Jenlpmqvi0ihwJjgW44XexfV9Wd8e3A22krpTEt\nXSRCcPaLNadzxXr3g0UL8K3Nh5NGpbFgjZPs3olv+WI0ECRy3JA6HkSIjDyHWK9+tL73KgCKbvo9\n0cOOLNsx5G8PpKDEpr48P2fmBq5nakxojEm5wCfv4tu1nVjXQ/BtXl9j+mivfgD4WmjLrPQS44mV\nPrZQG8mdbKLH1GZ8dJNJXp8z6wM8AAwAynwz3J6Nxph0USX4nxcACJ1xATlP/bHGXRKXGdeucKaG\nEUlL0fxub8HGJnGJcdCIBi1HprXk+3BeW2YzcKaAuQkorCGtMSaF/N8swr8un1jbDkSGnAoegpke\ncCCxth3w7duN7NiCHpjiQXAjYbKen0pWUu9KQiWQlZ3afOogJZcYTZPjNZgdBeSpaiydhTHGVBRv\nlYVHn+s9WIgQ690X31cL8a3LJ5rCYCa7tpEz5V78K5ei/kCiZ19w7quET/tByvKpq1ReYvSiJbeG\nGhOvvRnnUTrPmjEmQ+TbdQQWLUCDWU4wq4WYe9/Mn8LpYPxLPqX1XT/Bv3IpsY4HUXT7w4ltwdem\nQ0nZ3paZHo0/d/IIsp9xytTSLjG2dF5bZmuB/4jILGBL8gZVvSvVhTLGOLLecibcjOSdXmYSSS+i\n7n0zX4pGAgm+NI2sl/6RaPUUX/VraFtaJt+e7wjOeYnwuItTkl992CXGlsdrMGsDvAYEKTvJpjEm\njQLznXlxQ6fX/vJdrPf3AXeMxhR0Asme9RQqQujcyYTGXwo+f4U0Wa//k/Coc2r1kHJl4q25ul7C\nix49MCOXGE3j4fU5s8vTXRBjTEUSDhEZMBjt1qvW+2rng9HWufj2fofs3oke0LnWx/Ct+Kr0eLnt\nKL7qDqLOwX2aAAAgAElEQVT9K5/oPdqvP/4VXxF8+0XC4y+tdV5eVRboZMtGgu+/mXgfGTQypXna\nfbHGr8pgJiK9VXWtu3xoVelUdXUaymVMo5E7eQTH0XB/0MJnXFC3HUWI9u5H4OvP8a1bQbQ2wSwW\nI/jm82T9+++JVYX3/R3t1KXKXUrOu4LWD95A1pvPO/f3ctvVrdxelRQTWPgewXlv4F/+ZZlNkRNO\nSW/eptGprmX2FRB/zH0loED56xQKVLzWYIypkZdLadFefYkePsDzMcsfK9arL3z9Ob41K4gO8HgP\nqWAPOY89SODLj8qsri6QAcQOH0DkqBMILP2MrP+8QOgHV3out2eR0jER21w/EXGnYNGsHCKDRhD8\nwJ2EPqd16vM2jVqVwUxV2yYtt8xht41pCLHSJ2DCZ1xQr3tdiR6N6/I9jQTuW7mUnCn34tu1DW3T\nluKf/IpW/3u75/xCE39MYOlnBN/+N6FUdNNXRbZuJLDkU/xLPsW/7IvEJikuJHrYkYSHjSVy0kho\n1aY0mJkWx/NwVsaYzPAvWZhYru+9n9IejSs8tQRb/fZ6JBoletiRFF97d42tsfJi3zuKyIDBBBYt\nIOv1GXUsdGnrq/UvLyo7Z1iSwgeeItajT5l1dm+r5bLhrIxpZIKzZ5a+CdTv/5vapQea0wrfru2e\n0ks0SujMCwn94Cd1zjs08ccEFi0g+N+Xqk1XPrjKrm0E33udwNzSaW58O7aiue2IHDWQ6NEnEj36\nBNrc6NxDLB/ITMtmw1kZ04jIlg0EFn+cugP6fMQO6Yt/xeIqk/i//jyxXPQ/DxD1OJFlVWK9+hIZ\nOIzAp/M8pfd/9QnBd17Bv+hDJFZ2kKHCe/7m3PdroROMGu9sOCtjGpHgnFkpP2a0dzXBrLiQ7Cce\nKk1bz0AWF5pwOf7P3kdUAefBbSksgKL9SGGBs+xq9cdbAFC/n/CgkURGnUOr390IQKzP9ysc2y4l\nmsp4DWbx4aw+S2NZjGnZCgvKPCuVKvFOIJXJfv5v+HZsqXJ7nfPs0YfIyaMJLpgDQOu7flJ12s5d\nCI84m8jQM9EOnVJeFtMy2HBWxtSDbN+csmMFP/gPUlxE5PABBL5ZlLLjxqeDKc+/9DOC77xcZrDg\nVAqdd0UimEV7Hgatc9HETxuy3HuDhX+YUeloIsbUhg1nZUxdqZLz2G9Tc6xYLNHxI3zqeakNZl0P\nQYNZSDhUurKokOwnncuLofGXkj3zyZTlF5c87UzRb56osD0ezCyQmVSw4ayMqaPAh7PxJw33RCxa\n5z/M/sUf49v2LbHOXYgen+IBcv0BYocchn/VssSq7OcfxbdjK9Fe/QiPm5SWYGZMJnnueysi7YHv\nU7Fr/jupLpQxjV7RfrKen1pmlW/lUmL9jqnT4RKtstET0tJSifXqlwhm/qWfEnz3VdQfoOQnt9W7\n+39dVdeRI9OdPKxTSdPnqb+riFwGfAu8CjyR9PO414xE5AwRWS4iK0Xktkq2Hy4iC0SkRER+WZt9\njcm0rFn/wLdnF9HvHZVYF/jiwzodS75dR2DJQjQrh/DwcakqYhnRXqX3zbKf+AMAoQmXEetZ5bCr\nzU7BtLkWtJoxrw9vPAD8QFW7qGqfpB9Pvwki4gemAGcCRwIXi8iR5ZLtAq4H/liHfY3JGN+G1QRn\nv4iKj5JLb0isD3wxv07Hi987iuSdCm3a1pC6bpK7uPt2biXa5/uEx16UlryMaQheg1kAeLse+QwC\nVqrqalUNAc8B45MTqOo2VV0IFYaQq3Ff0/xkeoZiz1TJfuZhJBYjPOoc54Fel2/zBmTz+lofMvCB\nM2dZeMzElBWzvFj33ollDQQpufI28Ntodqb58BrMfg/cISJ1fQy/O7Ah6f1Gd1269zUmpQIf/Rf/\n8i/Rtu0JnXdFxe11uNQooWIiR52Q3uGZAsHEYujcy2woKNPseP2v2Y3AwcAtIrIzeUNjG5sxP7/m\nKeK9pGmpGkvdHOe+NobyxMuyaslijnz2EQDWDz+XXd9uAbYktgOEP/wv+X1PqNVxAdYdNZi9Seda\n3fnXtW7i+33d7wQot2998ktHWZubln7+1ampbvr2rfw5yfK8BrMfeUxXlU2UfT6th7su5fvWdOL5\n+fmeK6elaYx105jKc8RXHxAs2EP0sCPodN5kOpUbL1D9AdpsXEXfrgdBbnvPx40d2I0uZ06gSyW9\nGKs7/7rWTd/vH16nY9aUXzrK2hw0xt+rxiKVdeP1ObP36pnPQqCvO/r+JuAiYFIG9jUmZZxOH0LJ\nJTdUOvBt9PABBJZ+SuDLj4nkneb5uOEx6emOb0xL4nUKmPuq2uZlOCtVjYjIdcBbODNTP6mqS0Xk\nKnf7VBE5GPgUaAfEROQG4EhV3VvZvl7KbUwqSTRKeOQ5lQ5+CxA9bogTzL6YX2Mw860rvbQSHnZm\nSstpTEvk9TJj+SGsDgaGA56H+FbVN4A3yq2bmrS8BecSoqd9jck0zW1HyQ+urHJ75LghZD/7F/xf\nfQLhEASzqjiQkvXco6XvW+dWnq4ZsOe6TKbUeTgrETkDuDjlJTKmkSqZeAXktqtyu3Y+mOghh+Ff\nvwr/N4uI9h9UaTr/4k8IJM0hVlsWIIypqD4z3r0NnJuqghjTGCU/NxY5peb7YNHjnPnA/FV10Y9G\nyrbKjDEp4XU4q0PL/RwN/Iayz38Z0+wE5yeNFZDdqsb0kQHOIMGBL+aDOzFlssC8N/F/u5ZY0ojy\nxpj689oyWwnku68rgY+AocDkNJXLmIYXixKY/1btdundj1iHTvh2bce3fmXZjUWFZM1yRqcPnf/T\nVJXSGIP3YBZUVb+q+tyfXFUdqqo287RptvzLvsC3a3vtdvL5iLqts/KXGrPefA7fnu+IHnYkkUEj\nUlTK9LMBek1TUGMwcwf6LRCR7AyUx5hKL881hMD7/6nTfpHjki41umTXNoJvPg9AycXXgEj9C2iM\nSagxmKlqFFgBdEp/cYwBf3JPv0j5caczpGg/gc/er9Ou0SOPR7Ny8K9dgezaBkDWi08ioRIiJw4n\n1vfoVJbUGIP358ymA6+JyMM4A/0m/utsk3OaVPMv+bR0+fP5RBvgklzgk7lIqITo4cfi/+bL2u2c\nlU306IEEPv+ANjdeAICKOJNhnv+TNJS2/upzGbFg2lxnWKLUFceYWvMazK52X+8pt16BljO7n8kI\n/9LSYBac+1qtgll82pj63uMJxqdlOeWM2gczIHJcHoHPP0i8F1VCo89Fu1Q6LkBG2H0v05x56gBS\nbkLOWk/OaYxXsvc7/ElDPQWWfops35zZMmzdhH/FYjQrh8jA4XU6RnTAyWjSfTFtnUto/KWpKqIx\nppz6PDRtTMr5l1bsIBt87/WMliH+bFnkxOHQqnWdjqHtDiB22FGJ96Hxl1Y7eogxpn4smJlGxb9k\nYYV1gfffhGgkMwWIxQjMd3oxRk453dMuVXVdjxw/JLEcHm2D5RiTThbMTOOhin9J2ZZZrOsh+Hbv\nxP/lRxkpgn/5l/h2bCXWqQvRwwfU61jhkeeUvqlq0OEq2LNdxtSOBTPTaPg2rcW3ewex9gck1oWH\njwOcjiCZEPjAbZXlnVbpnGW10oxHwzemsbFgZhqNeJf86FEDE+vCeaej/gD+xZ8kntnKnTwi0Wsx\npYoLCSx8L5GvMabpsGBmGo14l/zo0SeWrmzXgcgJQxGNEZj3Ztryzp08gtyfjUVKion2PRo9uOG6\n0Btjas+CmWkcwiH83ywCIHrUCWU2RUa4lxrnvQGxaPqLcsoZac/DGJNaXh+aNiat/PlLnBE3ehyK\ndig7clr0iOOJHdgV3/bNZUYHSQcNZqV0EGDrxGFMZljLzDQKiftlRw+suNHnK+0IUo9nzqq61yY7\ntiSWIycMzVjHDeuxaEzqWDAzjUJpMDux0u2RoWeiPh/+pJHoKwiHSo/32fvI7p1Vpy0pJjD/bXJ+\n/wva3HRR6SFGnl27ghtjGgW7zGga3t7d+NetQINBot8/Bqh4eU47dCI6YDCBzysJZqr4P3uf7Oen\nJla1+sudAMQ6H0z0sCOJfa90NI7sJx5yBhIuLnR2DwaRsDM6f6yez5YZYxqGBTPT4AJfOw9KR/sd\nA1lVT5sXHn5WhWDmW/MN2TP+D/+KxWXWR448Hv/qZfh2bMG3Ywt8XDq5Q3DeG05+hx1JeOgZRAaN\nJPcaa5EZ05RZMDMNrqZLjHHRYwYR63hgYvbn7L/9luCHzjiK2rY9JRN+TM7Tfwag+NY/QSyKb9M6\nfCuX4l/1NcH3na79oXGTCJ9yOtqtl6fy2RQnxjR+ds/MNCzVxHiMlXb+SObzExk6NvE2+OHbaCBI\naOzF7H9oOpHR4yukj/U8lMjIsym58tbE6tAFP/UcyIwxTYO1zEyDks3r8X23g1i7A4j1qHlGofDw\nsWS9PM1ZHjTSCUwHdq13OaxXoTFNmwUz06AC8VbZUSd4GgtRO3VJLJdce3faymWMaVosmJkG5fV+\nWSpY68uY5svumZmGEw7hX1b5EFbGGFMbGQtmInKGiCwXkZUiclsl20VE/uJuXywixydtWysiX4nI\nIhFJ73hGJmP8K5cioWKiPfqgB3Ru6OIYY5qwjFxmFBE/MAU4FdgILBSRV1T166RkZwJ93Z+TgEfd\n17iRqrojE+U1mZHJS4zGmOYtUy2zQcBKVV2tqiHgOaBcP2rGA0+r4yOgg4jUv5uaSav6zC2W9dp0\nwC4xGmPqL1PBrDuwIen9Rned1zQKzBGRz0Tkp2krpUm5SoOdKr5VyxJvo98/NrOFMsY0O02lN+Mp\nqrpJRA4CZovIN6o6r7KE+fn5NR7MS5qWqrZ1c1wN+8W3r168iLZrltFu1RLarVpCsGBPaZ7rN1S6\nb13yrKk89WHfm+pZ/VTN6qZqNdVN377ext7JVDDbBPRMet/DXecpjarGX7eJyCycy5aVBrOaTjw/\nP99z5bQ09ambMvupIru24duwOrGq/59/gcRiifexDp3wuaPapyRPV7z7fao/YfveVM/qp2pWN1VL\nZd1kKpgtBPqKSB+cAHURMKlcmleA60TkOZyOH3tUdbOItAF8qrrPXT4NuC9D5TYeBWfPxLdxDb6N\nq/FtXJMYkT5BhMjhA4geM4ho/5OI9TyU3MtGNkxhjTHNTkaCmapGROQ64C3ADzypqktF5Cp3+1Tg\nDWAssBIoBC53d+8CzBKReHlnqOp/MlFuUw3VMnOLZT/7lzKbneGp+hD4+nMA9k95BVq1KZOmrg8x\n28PPxpjyMnbPTFXfwAlYyeumJi0rcG0l+60GrIdAI+JbsZjs5x/Dv3JJYl14+DhiPfoQ63EosR59\n0HYHAJR2/igXyIwxJpWaSgcQ41E8eKSj9eLbuIasfz9OwG2Rxdp2wLdvNwAlP7455fkZY4xXFswM\n4ATB46gkCEYj5P54DAAqPkRjaHYO4TMvJHTGheReNbbCsZLZJUFjTCZYMDNlxWL4NqzCv+wL52d5\n0gzOPiE04lzC4y9F23dsuDIaY0w5FswMqCYW21x3LrJ/b6XJCh+chnbpkalSGWOMZxbMWjJV/Es+\nJevff0+skv17iXXqQvSI4xI/bX5xgZPcApkxppGyYNZC+VYuJftfj+H/5ssy6/f/YYYzc7PzKARQ\n/X0vuydmjGkMLJg1V6plAhKU9nSMDBhCYNGHTrI2bQmNm0T2C39z3h/ULaPFNMaYVLBg1oz4Vi5N\nLOdeNhINBiGYhQayIJiV2BZY9KHTI/G0HxA680Jo0zYRzIwxpimyYNYcxGIE33yOrBefKLNawmEI\nhxH2l1kfOvU8wmf/sEyPxIJpc51x0jJSYGOMSS0LZk2c7N5J9mMPElhadgLugqfegUgYwiEkHIJw\niDa/vBiA0I9+3hBFNcaYtMnUfGamlrxMeun/6hNa3XklgaWfom3bU/SL35Vu9PkgKxvatEU7dHI6\ndRhjTDNlLbMmqHyQixx5PCU/vR09oHO1+1nPQ2NMc2XBrCkoKcK3YTW+dSvxr1+ZWK0+H6GJPyY8\n7mLw+RuwgMYY07AsmDVGRaUdNlrfdgmyZSOSNEpHItntfyHW9+hMlswYYxolC2aNSXEhwTmzyHrj\n+cQq3+YNqN9PtFtvYr2+R+yQ75E9YwqABTJjjHFZMGsMQiUE33mZ4GszElOqxBXe+xix7r3LPCcW\nPv38DBfQGGMaNwtmDSjekSPWoRO+3TsBiB52BKGJV9DqD790tvXuV6tjWicPY0xLZMGsAcj2zQQ+\neCvx3rd7J9Fe/QhNvJzosSdXGIbKGGNM9SyYZUpJMYFP5xH44D8Evv68zKain99P9IRTLIgZY0wd\nWTBLJ1VyLxvpLLZqg7i9FDWYRWTgMIIL5gAQHTi0wq52udAYY7yzYJYOxYUEFswh+N+XE6ukaD/R\nw44gPPRMIoNGQpu2lFx1RwMW0hhjmg8LZikk364j+N+XCH7wFlJcWGbb/t/+A+3eu2EKZowxzZwF\ns/qKhPF/MZ/gf18msOyLxOpov/6ER51LZODQMt3qjTHGpJ4Fs7qIxfDlLyG4YDbBd19NrNbsHCJD\nTiU8ajyxQ77XgAU0xpiWxYJZLfjWryLw0RwCH72Db+fWxPpojz5ERpxNOO80aJ3bgCU0xpiWyYJZ\nDWT7ZgIfv0NgwRz8G9ck1sc6HkRk8Ggig08l1vPQBiyhMcYYC2aV2bebwCdzCS6Ygz9/SWK1tmlH\nZNAIwoPHOOMi+mw6OGOMaQwsmMVFwgQWvue0wJYsRKJRADQrm8hxeUQGjyHa/0QIBBu4oMYYY8rL\nWDATkTOAhwE/8Liq/q7cdnG3jwUKgctU9XMv+6aEKtlP/y9SWID6fESOPZnIyaOJHJ8HOa1Tnp0x\nlSmOKNuLo+wojuET6Jjto2O2j9YBQWyEGGOqlJFgJiJ+YApwKrARWCgir6jq10nJzgT6uj8nAY8C\nJ3nct/6CWYTO/pHTI/HEEdCuQ60PEVNlS2GM9QUR1hdE2VAQ5dvCKEURJRRTSqJKKAahqPM+HIUY\nFecpi/Mh+H0QEAj4hICA330N+IRAuW1Bn5PeJ6AKStIrzj8lMWVvSNkXjrE3FGNfWNkbirE3pBRF\nWxFcsImge+ygTwhK6XLAPXZAypYrOb+45PcBH2T5hCy/kJW0HPQ5dVEUVYoiSmHEeS2KKIVRJVqu\napL/lPultExZPudYAZ8Q9AmqSkkUiqJOnRdFleKIsxyOOaOG+XBeyywjKFqm7PF601gr2n+xhbZB\noW3QR7us0tfWASEcI3EeRZHSPIujil8gyy9k+4VsX+ly0CfsDcXYURxje1GU7cXO51GZbL8T2A5w\ng1sbN7j5xKkXn+AuC0LZkdGS6y2qJMoWL2eRW85wTMkNCu2yfLRzz81ZFtoEfIRjTrpQTCmOQknU\neV8SVfbvz6bV6h2JOkuefi/+GWX7hSw/ZCctq0I4hvP7EINwTAm5n5O65+UX8Lvn6nfPM6Yk8i6K\nlP28o6rk+IVWASefVn5JvA/6hKgqkRhEEq8QjWmZ15hCVJ3vYCTm/G4n/54F4t8393fSJ2XrOvkz\nKCzIot3GXc7vZuKzKt2nzO+paunvK855+5PqwO9+Z/2J7664310S3wVxfx9jSb/7zrKW+VzKf9Pi\n+0Tc83bOXd26KD22zz1B5/em9DySjxPXrbWfO05oV+l3OtUy1TIbBKxU1dUAIvIcMB5IDkjjgadV\nVYGPRKSDiHQFenvYNyVu6HgG3+6PwsIosLPMtoqBQRPriqPKxv1RNu6PEo6lulSZJESjzvk4qg60\nzVN15yvs2x9NewmCPjgwx8e3hc4XqWtrH7tKYpREYXNhjM2F6f2CbS0CqMt5+mF3SYpL01wEYGdR\nQxeiwfxx8T52X9497flkKph1BzYkvd+I0/qqKU13j/sm5Ofn11iYqtLMXpvD2qL6deroGFS6Zsfo\nmqN0y1a6ZCs5fiVLnD9UQYGgz3kfcP+nVpl40Iz/DynxPyWVpOXkV0ksx/8kl/9fouDknxtQ2viV\n3ACJ11y/ku0rPV7iJyaJ5Xi+MZz8Sv8H6+Yl5fJ0lyMKYYWwe6xQzH2vTnla+SHHp+T4IMddzvY5\n//tMro/k5ViifE6dhJPKLEC2D7Ld42QlLQek4n9Mkv8HW9n/rMVNUxgV9kehIOK87o8KBRGhMOrk\nkeNTsv1O3jlJeccUQokWiLh14azL9cMBQeWAoNIxS2nrrzjetCoUx2BPRNgTdl6Lo1Km/OB8LjEt\nt2+571W8bnJ8So7fqZccH4n6LoxCQVQoiJS+7o865xiUeF1CVtL5BUUrbZkgTgGcz1+czz0Wrwvn\nvdPSh4BPnd8N9/ciIKX1Hkt6japTHyKln3G8TPHPWICSGJTExH0tXQ7HnGP7y/0ERBPLPnHuZ5S2\nCtUpC6W/D2V/TyTxOVX2XY1/RuXPJf5ZJX5XpHQ5Xo/RxL5C1M0j/juY/B0u/51OHEdK/8Yk/x2o\nSkDiLUclQGmLP/67mPhdcV/j5ahwTHdFWz+M7Byt9u9yTX+z+/btW+32RNk9pWpCajrx/Pz8KtP8\nuXUxBZGq/3ee/AUpfXUucfVo46dHrp/WgabbwzE/P58jPX5xWpr8/HwGHmkPwlelut+rls7qpmqp\nrJtMBbNNQM+k9z3cdV7SBD3smxIju+ek47DGGGPSLFPNiIVAXxHpIyJZwEXAK+XSvAJcKo6TgT2q\nutnjvsYYY1qwjLTMVDUiItcBb+Fcjn5SVZeKyFXu9qnAGzjd8lfidM2/vLp9M1FuY4wxTUPG7pmp\n6hs4ASt53dSkZQWu9bqvMcYYE9d0eysYY4wxLgtmxhhjmjwLZsYYY5o8Ua36uaqmYs+ePU3/JIwx\nxlSrffv2VT7zbS0zY4wxTZ4FM2OMMU1es7jMaIwxpmWzlpkxxpgmr8UEMxE5Q0SWi8hKEbmtocuT\nKSLypIhsE5ElSes6ishsEcl3Xw9I2vYrt46Wi8jpSetPEJGv3G1/kSY+U6SI9BSRd0XkaxFZKiL/\n465v8XUDICI5IvKJiHzp1s+97nqrH5w5GkXkCxF5zX1v9eISkbXueS0SkU/ddemvH1Vt9j84w2Ct\nAg4FsoAvgSMbulwZOvdhwPHAkqR1DwG3ucu3Ab93l4906yYb6OPWmd/d9glwMs6EAW8CZzb0udWz\nXroCx7vLbYEV7vm3+Lpxz0mAXHc5CHzsnqPVj3NOvwBmAK+5761eSutmLdC53Lq0109LaZklJgdV\n1RAQn+Cz2VPVecCucqvHA9Pc5WnAuUnrn1PVElVdgzNO5iB3ktR2qvqROt+yp5P2aZJUdbOqfu4u\n7wOW4cyd1+LrBpzh5VS1wH0bdH8Uqx9EpAcwDng8aXWLr5capL1+Wkowq2riz5aqizozEgBsAbq4\ny9VNkLqxkvXNgoj0Bo7DaX1Y3bjcS2mLgG3AbFW1+nH8L3ALztyUcVYvpRSYIyKfichP3XVpr59m\nNzmnqR1VVRFpsV1aRSQXeBG4QVX3Jl+Wb+l1o6pRYICIdABmicjR5ba3uPoRkbOAbar6mYiMqCxN\nS6yXck5R1U0ichAwW0S+Sd6YrvppKS0zL5ODtiRb3WY87us2d31V9bTJXS6/vkkTkSBOIJuuqjPd\n1VY35ajqbuBd4AysfvKAc0RkLc7tilEi8ixWLwmqusl93QbMwrnNk/b6aSnBzCb4LOsVYLK7PBl4\nOWn9RSKSLSJ9gL7AJ+7lgb0icrLbo+jSpH2aJPc8ngCWqeqfkja1+LoBEJED3RYZItIKOBX4hhZe\nP6r6K1Xtoaq9cf6OvKOqP6KF10uciLQRkbbxZeA0YAmZqJ+G7vmSqR+ciT9X4PSW+XVDlyeD5/1P\nYDMQxrnufAXQCfgvkA/MATompf+1W0fLSeo9BAx0v5SrgL/iPnDfVH+AU3Cu7S8GFrk/Y61uEud0\nDPCFWz9LgLvc9VY/pec1gtLejFYvzjkditM78UtgafxvbSbqx0YAMcYY0+S1lMuMxhhjmjELZsYY\nY5o8C2bGGGOaPAtmxhhjmjwLZsYYY5o8C2bG1EBEbheRx2tOmdI8fyMiO0Rki8f0KiLfS3e5jGms\nrGu+MY2MiByC88xNL3VGUfCyjwJ9VXVlDel6A2uAoKpG6lnUBs/HmDhrmRnT+BwC7PQayIwxFsyM\nSRCRW0Vkk4jscycKHO2uv8cdfw8R+auIFCT9RETkHndbNxF5UUS2i8gaEbm+mrzai8jTbtp1InKH\niPhEZAwwG+jmHv8fVex/s4hsFpFvReTH5baNE2fiyL0isiFePtc893W3e/zBInKYiLwjIjvdS5vT\n40NZ1VAvPhG5TURWufu+ICIdq8rH04dgTB1ZMDMGEJHvA9cBJ6pqW+B0nEkGy1DV61Q1V1VzcYbE\n+g54WUR8wKs4w/h0B0YDN0jSzLnlPAK0xxn+ZzjO2HOXq+oc4EzgWzefyyop6xnAL3HGS+wLjCmX\nZL97vA44825dLSLxuaCGua8d3OMvwJn88EGgG3AEzsCv93iol5/jzDE13N33O2BKNfkYkzYWzIxx\nRHFmuz1SRIKqulZVV1WVWEQOBF4Cfq6qXwAnAgeq6n2qGlLV1cDfcQajLb+v313/K1Xdp6prgf8H\nXOKxrBcAT6nqElXdjxt44lR1rqp+paoxVV2MMz7n8KoOpqorVXW2OhMkbgf+lJS+unq5CmfsvY2q\nWuKW4wciYlNLmYyzYGYMzh904AacP8jbROQ5EelWWVpxpo75NzBDVZ9zV/fCuTS4O/4D3E7pJITJ\nOuPM3Lwuad06vE/O2I2yExomHwcROUlE3nUvYe7BCTqdqzqYiHRxz3eTiOwFno2nr6FeeuHMcxY/\n32U4wa+yczYmrSyYGeNS1RmqegrOH2kFfl9F0keAvcAdSes2AGtUtUPST1tVHVvJ/jtwZjHolbTu\nELzPZ7WZsnNAHVJu+wycqTV6qmp7YCrOpURwzqu837rr+6tqO+BHSemrq5cNOKOcJ59zjjrzWVk3\naVfmWrwAAAFGSURBVJNRFsyMwbk3JCKjRCQbKAaKgFgl6X6Gcwnuh6qavP0TYJ/bWaKViPhF5GgR\nObH8MdSZwfkF4AERaSsivYBf4LSIvHgBuExEjhSR1sDd5ba3BXaparGIDAImJW3b7p7XoeXSFwB7\nRKQ7cHPS+VZXL1Pdc+jlpj1QRMZXk48xaWPBzBhHNvA7nFbTFuAg4FeVpLsY5w/0t0k9Gm93A9RZ\nwACc56t2AI/jdPKozM9xOmqsBj7AaU096aWgqvom8L/AO8BK9zXZNcB9IrIPuAsn+MX3LQQeAOa7\nlwdPBu4Fjgf2AK8DM5OOVV29PIzTAnzbzesj4KRq8jEmbeyhaWOMMU2etcyMMcY0eRbMjDHGNHkW\nzIwxxjR5FsyMMcY0eRbMjDHGNHkWzIwxxjR5FsyMMcY0eRbMjDHGNHkWzIwxxjR5/x/VtAJkhYcw\nwAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb6463f2a20>"
      ]
     },
     "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": 41,
   "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": 42,
   "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": 43,
   "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": 44,
   "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": 45,
   "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": 46,
   "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)."
   ]
  }
 ],
 "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": {
    "3c215f4376c7471c8388a5f47e675a4d": {
     "views": [
      {
       "cell_index": 9
      }
     ]
    }
   },
   "version": "1.2.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}