Skip to content

Instantly share code, notes, and snippets.

@hrit-ikkumar

hrit-ikkumar/utf-8''week1.ipynb

Created March 25, 2019 10:12
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save hrit-ikkumar/26daa175dfefafe1ef6a4530436ff7f9 to your computer and use it in GitHub Desktop.
Save hrit-ikkumar/26daa175dfefafe1ef6a4530436ff7f9 to your computer and use it in GitHub Desktop.
Download ZIP
Display the source blob
Display the rendered blob
Raw
utf-8''week1.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Mean/Covariance of a data set and effect of a linear transformation\n",
"\n",
"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": [
"First, let's import the packages that we will use for the week"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# PACKAGE: DO NOT EDIT THIS CELL\n",
"import numpy as np\n",
"import matplotlib\n",
"matplotlib.use('Agg')\n",
"import matplotlib.pyplot as plt\n",
"matplotlib.style.use('fivethirtyeight')\n",
"from sklearn.datasets import fetch_lfw_people, fetch_olivetti_faces\n",
"import time\n",
"import timeit"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"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. \n",
"\n",
"__When you implement the functions for your assignment, make sure you read\n",
"the docstring for what each dimension of your inputs represents the data points, and which \n",
"represents the dimensions of the dataset!__. For this assignment, our data is organized as\n",
"__(D,N)__, where D is the dimensionality of the samples and N is the number of samples."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Shape of the faces dataset: (4096, 400)\n",
"400 data points\n"
]
}
],
"source": [
"image_shape = (64, 64)\n",
"# Load faces data\n",
"dataset = fetch_olivetti_faces('./')\n",
"faces = dataset.data.T\n",
"\n",
"print('Shape of the faces dataset: {}'.format(faces.shape))\n",
"print('{} data points'.format(faces.shape[1]))"
]
},
{
"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).\n",
"\n",
"We have created two function which help you visuzlie the faces dataset. You do not need to modify them."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def show_face(face):\n",
" plt.figure()\n",
" plt.imshow(face.reshape((64, 64)), cmap='gray')\n",
" plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "05ffd545d42e47edaa1c162d4cfb6123",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"interactive(children=(IntSlider(value=0, description='n', max=399), Output()), _dom_classes=('widget-interact'…"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"@interact(n=(0, faces.shape[1]-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"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this week, you will need to implement functions in the cell below which compute the mean and covariance of a dataset.\n",
"\n",
"You will implement both mean and covariance in two different ways. First, we will implement them using Python's for loops to iterate over the entire dataset. Later, you will learn to take advantage of Numpy and use its library routines. In the end, we will compare the speed differences between the different approaches."
]
},
{
"cell_type": "code",
"execution_count": 42,
"metadata": {},
"outputs": [],
"source": [
"# ===YOU SHOULD EDIT THIS FUNCTION===\n",
"def mean_naive(X):\n",
" \"\"\"Compute the mean for a dataset by iterating over the dataset\n",
" \n",
" Arguments\n",
" ---------\n",
" X: (N, D) ndarray representing the dataset.\n",
" \n",
" Returns\n",
" -------\n",
" mean: (D, ) ndarray which is the mean of the dataset.\n",
" \"\"\"\n",
" N, D = X.shape\n",
" \n",
" mean = np.zeros(D)\n",
" for n in range(N):\n",
" for i in range(D):\n",
" mean[i] = mean[i] + X[n, i]\n",
" \n",
" mean = mean/N\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",
" mat = np.zeros((N, D))\n",
" mean = mean_naive(X)\n",
" for i in range(N):\n",
" mat[i] = X[i,:] - mean\n",
" for i in range(D):\n",
" for j in range(D):\n",
" covariance[i, j] = covariance[i, j] + mat[:,i]@mat[:,j]\n",
" return covariance/N\n",
"# 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": [
"Now, let's see whether our implementations are consistent"
]
},
{
"cell_type": "code",
"execution_count": 43,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"X:\n",
" [[0 1 2]\n",
" [3 4 5]]\n",
"Expected mean:\n",
" [[ 1.]\n",
" [ 4.]]\n",
"Expected covariance:\n",
" [[ 0.66666667 0.66666667]\n",
" [ 0.66666667 0.66666667]]\n"
]
},
{
"ename": "AssertionError",
"evalue": "\nArrays are not almost equal to 7 decimals\n\n(shapes (3,), (2, 1) mismatch)\n x: array([ 1.5, 2.5, 3.5])\n y: array([[ 1.],\n [ 4.]])",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-43-b07e0fb49c61>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Expected covariance:\\n'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpected_test_cov\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 10\u001b[0;31m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtesting\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massert_almost_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmean\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_test\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpected_test_mean\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 12\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtesting\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massert_almost_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmean_naive\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_test\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpected_test_mean\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/opt/conda/lib/python3.6/site-packages/numpy/testing/utils.py\u001b[0m in \u001b[0;36massert_almost_equal\u001b[0;34m(actual, desired, decimal, err_msg, verbose)\u001b[0m\n\u001b[1;32m 561\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mactual\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m\\\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 562\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdesired\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mndarray\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlist\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 563\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0massert_array_almost_equal\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mactual\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdesired\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdecimal\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0merr_msg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 564\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 565\u001b[0m \u001b[0;31m# If one of desired/actual is not finite, handle it specially here:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/opt/conda/lib/python3.6/site-packages/numpy/testing/utils.py\u001b[0m in \u001b[0;36massert_array_almost_equal\u001b[0;34m(x, y, decimal, err_msg, verbose)\u001b[0m\n\u001b[1;32m 960\u001b[0m assert_array_compare(compare, x, y, err_msg=err_msg, verbose=verbose,\n\u001b[1;32m 961\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Arrays are not almost equal to %d decimals'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mdecimal\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 962\u001b[0;31m precision=decimal)\n\u001b[0m\u001b[1;32m 963\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 964\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/opt/conda/lib/python3.6/site-packages/numpy/testing/utils.py\u001b[0m in \u001b[0;36massert_array_compare\u001b[0;34m(comparison, x, y, err_msg, verbose, header, precision, equal_nan, equal_inf)\u001b[0m\n\u001b[1;32m 713\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 714\u001b[0m names=('x', 'y'), precision=precision)\n\u001b[0;32m--> 715\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mAssertionError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmsg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 716\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 717\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misnumber\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0misnumber\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mAssertionError\u001b[0m: \nArrays are not almost equal to 7 decimals\n\n(shapes (3,), (2, 1) mismatch)\n x: array([ 1.5, 2.5, 3.5])\n y: array([[ 1.],\n [ 4.]])"
]
}
],
"source": [
"# Let's first test the functions on some hand-crafted dataset.\n",
"\n",
"X_test = np.arange(6).reshape(2,3)\n",
"expected_test_mean = np.array([1., 4.]).reshape(-1, 1)\n",
"expected_test_cov = np.array([[2/3., 2/3.], [2/3.,2/3.]])\n",
"print('X:\\n', X_test)\n",
"print('Expected mean:\\n', expected_test_mean)\n",
"print('Expected covariance:\\n', expected_test_cov)\n",
"\n",
"np.testing.assert_almost_equal(mean(X_test), expected_test_mean)\n",
" \n",
"np.testing.assert_almost_equal(mean_naive(X_test), expected_test_mean)\n",
"\n",
"np.testing.assert_almost_equal(cov(X_test), expected_test_cov)\n",
"\n",
"np.testing.assert_almost_equal(cov_naive(X_test), expected_test_cov)\n",
"print(mean_naive(X_test))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now test that both implementation should give identical results running on the faces dataset."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [],
"source": [
"np.testing.assert_almost_equal(mean(faces), mean_naive(faces), decimal=6)\n",
"np.testing.assert_almost_equal(cov(faces), cov_naive(faces))"
]
},
{
"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": 45,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAPcAAAD3CAYAAADBjMJTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+MHdWV57/H/cPdbbfbbkMTB1vYg5txEKxJhLL8SEYzMBllsyNAERmB0MasUDZKZhGrYbU4s9FKI+1qGa0yySgZZXYECY7EQFiGLBaazMRiyGwQCwQSjE0gNAskGBvb4O62u9v9k7t/dNXLqeM65533+vVrKM5Hslz1btWt+6rq9vuee849l1JKCIKgeqxa6QYEQbA8ROcOgooSnTsIKkp07iCoKNG5g6CiROcOgoqypM5NRJ8mol8S0StEtLtVjQqCYOlQs35uIuoA8DKATwE4BOCnAG5MKf0CAMbHx8OBHgRtZGBggPj+Un65Pw7glZTSqymlWQD3A7h2KY0LgqB1LKVznwvgDbZ/KPssCIL3AJ1LOJdKPiuV4iMjI0u4TBAEGsPDw2rZUjr3IQBb2P5mAIe1g++6667C/tlnn126DQCDg4Ol2+vXry8ct3bt2tp2b29voay7u7u23dXVBQA4fPgwtmzZUjiuo6OjdBsAVq1aVVrGP5f7RMW/eXw/396/fz927tx5xrEerDGSd999t6HzDh48iIsuusjdDu9x/FryuryN+fZLL72EHTt2lJbJ7XplCwsLpdvz8/OF4/h+ftypU6fQ39+Pubk59TxeNjMzU9uenp4uHDc1NVXbPnnyZKFsdHS0tn38+PHa9okTJwAAN998M+65557avjyHH/ujH/0IGkuR5T8FMExE24ioG8ANAPYuob4gCFpI07/cKaV5Ivr3AP4RQAeA76SUXmhZy4IgWBJLkeVIKf09gL/3HJtL49qFOztLt+WxvMySw94yeZwlt7Wypcry/Jxm5LAly3k75HHaeR0dHWo7GjEbtDZK2azdD2kSWdfm9Vv3m2M9M95GImrJe6WZdEDxnebmo+wTvB80cn9qbah7RBAE70uicwdBRVmSLG8EKSv4vpTlmsy1pIhXnlk0Mxpctr/U+r3S3mqXt02W3JPXsuSwdj3vd/F+r3pl2j3gI+dAUYrLc7Qyue99npa5Z8lyy2QMWR4EH2CicwdBRYnOHQQVpW02t7QZuH1hRYY1awdrUUzS9tLOke1oxt1lnTc7O+seW5AuNK3MwuvuamRcg6PZulZ0GY/+mp+fX1J0Wdk+3+aRZbIO3qaZmRm1jfWureF1k1k2txyXkvul13W1LgiC9x3RuYOgorRNllsRapZsseDyT8onrUwG+Hslr9dU8MryiYmJpqLoLJei1Q6tbGFhQTUJvG4mQDd9LNksnwuXztq23Lcmjnjle15HV1cXJicnC2Xye2rXtt4/y5WnvWOrVq0q7IcsD4KgRnTuIKgoKybLregbbbRVSisukbwjl+Pj42pZM1Ft9erQZNfx48dVuS0lV9nc9LJ9S9pz+LXm5ubU723Jcq/ctiQ1nw89Ojpa2J+dnS09TpZZkWdl0rvsu+QMDQ1hdHTUnBev1WmZKVKyazQSoaaZY5z45Q6CihKdOwgqSnTuIKgobbO5pR1p2QyabSdtL8tlUWYP9ff34+2331aPk2g2lRV1JSmzqwcHB3Ho0KGCvcy3V69eXaijp6endFsea03u1+z7mZkZ1ea23ExWxJdlL3NXZJ5nrKenB0ePHsXp06dLj5PuS77fyGwvTtlYyNDQEI4ePeqeGafVV+/a/N7xMjlmslR3cfxyB0FFic4dBBVlxWS5le+LSy0u8SyXgmcCQX9/P44ePVo4TptAIMusqCt+nuWyyO/B4OAg3njjjYKktqT3mjVratt9fX2FMp7SmZ9nRQTyNk1NTakSz5o0wZ8LoMtoLrUBYHJysnBtANi6dSveeuutQhk/T9bhdYVpedKAonuR348TJ064k0hwvFGVgG7eSClvRaiFKywIPsBE5w6CihKdOwgqyoola+BIe4jbVNK241juBs11derUqcJx1uwjmUxAO85yhZXl5t65cyeOHDlSsLm57Sztam5z822guKQSL5N2O7fBub05PT2t2pXye3K3lnRxcXuZ32NpL09MTJyxvXXrVhw5cqSwBI/mFgPsd4JjJQTh7efHTUxMmLPwON4klFboqzaOMTs7WzjOOxuQE7/cQVBR6nZuIvoOER0jooPss0Ei2kdEI9n/G5a3mUEQNIpHlt8D4FsAvsc+2w3g0ZTSnUS0O9u/w6rEiuCRLhdtRo1XBsnr8fOkpLNWdNRmOsk6LDcZP49L41OnTqnuIysiS15bc9dJKchlOr8fMzMzqltSynLLxcWluLYNFOU7l+hjY2OqLLfut0SbXcdNEUA3pebn580EHp46APs+au4vGdnn/Z4adX+5U0r/B8AJ8fG1APZk23sAXFf3SkEQtJVmbe5zUkpHACD7f6h1TQqCoBWQR+oS0VYAj6SULsr2x1JK61n5aEqpYHePj4/XKh4ZGWlVe4MgYAwPD9e2BwYGCrZvs66wo0S0KaV0hIg2AThW74R9+/YV9rnrR4bWaXZkI/nNy2zuyy+/HE888UThOO/MMt4mKxzSY199/vOfx/e+973C9+b2OL83QNHF1d/fXyhbt25dbXtwcLC2zV1kQNG9ltuf09PT6OnpKbSDf2dp6/Lvze1jYNFmLts+efJk4Thug+d1XH/99XjwwQcL7ilu38uxEMsm1paGljZ32XHXXHMN9u7day4vzWnm3QGK964sv37+fljvXJ5R6Lvf/a7avmZl+V4Au7LtXQAebrKeIAiWCY8r7D4A/xfAbxPRISK6BcCdAD5FRCMAPpXtB0HwHqKuLE8p3agUXd3IhayZX5b7yOuW8Ep2GbmltUnuc/lkufVkHTJHOK9PSvgc+V0seaklRpBRbt4kBtY9tpb44e3gEpK7vgDdrTc7O+u+31YEGb8/1iw5LbFkT09PoU7rWWtLVtUr42jvupVcAoi85UHwgSY6dxBUlBWbOGLJLm8dVp5uLerKG8QP6JJXjhTzEWBrkoOcwMK/N5eTlnvSSgbB65Aj7nySiiUZrWvzY+X31EbS5b2yZDm/x1oUF2DLcn5/+LUs2cwl+vz8vJl8Q5PY3uNkm7m81vKpefZL21T3iCAI3pdE5w6CihKdOwgqyorZ3N4lZvm2dxYOULSfuS0n85bz86R9ODo6Wtvma4zJ47Q1rgA9Yur48eOF47i9LKPQBgYGoMHvK7ezG5lJZblgOJbLT1sDzEqkKBM8cLvdWqZXzmrjaK5TaZvze7Vhw28ip998881CRKDMIc+vXZaIo147yo7NsaIxmyF+uYOgokTnDoKK0jZZ3oh7xxuVZiVQ0PJ4vfXWW4Xj+HlchgPAK6+8UnqeFYUm+fCHP1zb/tCHPlTbTimp8lJ+Z2tSiZbnzTJhpAtRi+bzLkkk0Sb+yDKZk42bO1qSC1n26quvFsr4RBX+XaR85+0/77zzAABf+MIX8Mwzz9T2gTNz1vG2cDNLPhe+73Xh8rp7enoKJky4woIgqBGdOwgqSttkuZSufHRYlmkRQtYIrZw3/M4775SWvfnmm4XjpOzi8BxfXLJb+cnkxA4+3/rcc88tbPN2cdNBTgrQlh0qOzbHmhCSnzMzM4POzk5VYsuIPWsVUU3OW8sayVVOeRv5feSmjfwuhw4dKpQdO/ab1AL8/ZARdfz+cG/E2NgYNm/eXNuXy0/xZ8/nzEuPBr8HchKPdr95mzo7O1UTwKqDE7/cQVBRonMHQUWJzh0EFWXFbG5reVLNFSaTG1jL+2pLEsmII24HSzuGR7NxV4pMQMDdHhs3biyUXXDBBbXt888/v7a9ffv2gj3HI+Ckm4Pb7evXry+U8e/Dt63ZUjJqjtt2VrQaHxewEkpwG9ObZEBG5fGxBW4DA8VxEjmGwm1rPu4ix0J4Hfy5nH/++di2bZtaP7fp+ZiMfK/4+JB8ntoMPU53d3fh3skxAzmWUUb8cgdBRYnOHQQVpW2yXMoxK7qHSygu/2SUEUemfuWyiMva7du3F47jslxKN24G8MkFUpZzicTTCwMoRDtxudrX11eo03KncSkuZTm/tibRZZ1yW1tqSEpGXqeVKrgsjXK9sqGhoYJrSXOLAcV7deGFF6pt5JNzpEnH7z2X5RdeeCHOOecc9dpcznNzT7q7eDukhOZlWj671atXF9596QL1rHQav9xBUFGicwdBRYnOHQQVpW02t2V3cJsYKNrg3B0g7QzLxpQuqZytW7cW9rkNJdu4Y8eO2vZZZ51V25ZuCW4fynZwO1KG1fJr8+8sxyC4m0guE6RdW9ah2fTd3d2Ffd5GGe5r5U/n9Vuz07hdyevYuHFjwW61bEpeNjRUXIOSt4PPDJSuUn4tPk6yefNm873iz4K3Q7q7rIQP/D3T7nd3d7cZxivvfxnxyx0EFcWznNAWInqMiF4koheI6Lbs80Ei2kdEI9n/G+rVFQRB+/DI8nkAt6eUfkZE/QCeJaJ9AG4G8GhK6U4i2g1gN4A7tEqkjOCyUboR+L4mYYCi1OKuDaDo+uDHnX322YXjrIgsPtOHSyspGb15wOXMKc0kkPeDS00pt+VSODlythu//zJCzZt3zFpVVZPlUg5rubn7+voK19PyvQN29Bc38aQJw9FMjDVr1hTKpAtK3n+tHVa0oJaLjpt73d3dhXsn+09LZHlK6UhK6WfZ9ikALwI4F8C1APZkh+0BcF3dqwVB0DYasrmJaCuAjwJ4CsA5KaUjwOIfAABD+plBELQbsnKbFQ4kWgvgnwH8t5TSQ0Q0llJaz8pHU0o1bTw+Pl6reGRkpIVNDoIgZ3h4uLY9MDBQsJVcrjAi6gLwdwDuTSk9lH18lIg2pZSOENEmAMf0GoADBw4U9rk9K+1lbi9ym8e7vCxQnju7p6en4B4BijahDFHk+9wGlMdpS88CRds0t+UuuOACvPzyy2pIqHSdWPYbr58fJ+1NbsfnZa+99hq2bdtWsPf5H3s5tsAzx/AZV0AxUw2fLWWNQeT3bXh4GCMjI6qtLu+319XGkWMEZVlOdu7cif379xdsbivnOD/OmiUnr60licxt7osvvhgHDhwo2OBjY2OFOvL3+I471GEu12g5AbgbwIsppb9gRXsB7Mq2dwF4uF5dQRC0D88v95UA/g2AA0T0XPbZnwK4E8ADRHQLgF8D+NzyNDEIgmao27lTSo8D0BJzX+29kJV32ZI0mnSVdVgzmHJZNzMzc0Y0nLX0DZeD/NrW0jzye2oJATdt2qQuNWR9Tyk7eVu4m0bWoeUc7+joUO+jlJPacsFA8bt53DSy/v7+fnWGlHSn8edimUiaGxLQ70dfX5+ZmFB7npa5JJ+ZFn0nZ4VZudu1pBeciFALgooSnTsIKsqK5VCzyjTJ1Ig04cfmUq2eLLeWvrGi0KwleLR84WeddVbpSLo8B7BHsDU5L9uh5UmT994q00aKZf3a6LusUy6TpI0wWxFqUrJbnguOtqTPhg0bzO+ptdHKI26ZcbyMmxgdHR2q2SavrRG/3EFQUaJzB0FFic4dBBWlbTa3FV1m2WVel5m17K22vhNQjICTtpFlZ2vttWzYvI3Hjx8/w/a3bChuZ1t5wK0ldj1Lvsr6re8i6+f2vjYjz6qjo6Oj8Cy0MQJAn1UF6M/M+45t2LDBXEJamxlnrXcn28j3+bOV98PbDo345Q6CihKdOwgqyorJciu6TIuSamRJIg05+d5ynWgRUxJLUmvnyckhHBl1xeuwTBNe5nVxNeKi5FiReNaz0KLGLHNDukC9JoYly8uOm56expo1a8z6tXtiRQ7KMs18ku45K8IuXGFB8AEmOncQVJTo3EFQUdpmc0t71pt4wWs3yfo0m0Tag9zmsdxd3ow18jht7a2Ojg7VpSPdR9ZMLa/bsNVYoalem1WOaXhdj1aSBGs2FqesLLe5vWMt1rpqWpssrOt63z9O/HIHQUWJzh0EFaVtstzKX21F93C3kBXpI2WRJkstiWRJWa88k8dZbj0tl1YjslzDe5ykWfNDm+Em3XqamZVSUk0TaUrJZZm0dlh51jnS1eiVx5YUt8o87Ugpud8rjfjlDoKKEp07CCpK22S5NxECoE8MkHLVK/+saDhOM6Oact+KtpO5v06fPl3b57Lcup51bW875LZ2nmVKWfV7osHKPufXs1bQ1NoE2BFfGo2YMNo9kO93M+aNd1kqwF4FNSd+uYOgokTnDoKKEp07CCrKikWocRvFmnBvTczn+143lrQjOVb0l2X3Wt+F29ncrp6cnCzY3NaytFYb+ff21sGRthv/btYyPtb3tGZEeZ8tr6+RRAhaLnFrVqJMQGnZ4Nq1ve+wVZ/lGrSWjdaIX+4gqCietcJ6iOhpItpPRC8Q0Z9ln28joqeIaISIvk9EvmUmgiBoCx5ZPgPgqpTSRLba5+NE9EMAfwLg6yml+4norwHcAuDbWiVS6lguLk3+NZK/mksyS/55I5esOjQXjtzn25OTk4V9KyGDlUtLq8OK2OPtnZ6edrt3tFVPgaLJYeWC13J4z8/Pq1GLVvIKS/ZbslxLbDE/P+9+JyzXlbVKKb932nFzc3Pu/OwadX+50yL5mqxd2b8E4CoAD2af7wFwXd2rBUHQNsjjbCeiDgDPAtgO4K8A/A8AT6aUtmflWwD8MKV0UX7O+Ph4reKRkZEWNzsIAmBxbfOcgYGBguRwjZanlBYAXEJE6wH8AMBHyg6z6njuuecK+/39/bxRhTKefphv8wXkgWIeMplnq0yS/epXv8J5551XOG4lZPnc3By6urpaIsu15Zaseev59rFjxzA0NNSULOcj/cCimZHTqCwfGhrCsWPH3N/Fyq+mzfX2yPKJiQmsXbu2KVneyEqkU1NTte2JiYna9smTJwEAO3bswEsvvVQoO3HiRKGO0dFRAMBXv/pVta0NucJSSmNE9GMAlwFYT0SdKaV5AJsBHLbOtcLpLJvbyoFdth5Y2fX4g7XyaHvDSht5kJZNxfct14aV+FCzMWUCRi1p4cmTJ90hrFbn5i/s9PS0Wgevnz+XmZmZMzpqjjdBhdxvZmac5SoF/KHR1vgE39feHfl+LIvNTURnZ7/YIKJeAL8P4EUAjwG4PjtsF4CH614tCIK24fnl3gRgT2Z3rwLwQErpESL6BYD7iei/Avg5gLuXsZ1BEDRI3c6dUnoewEdLPn8VwMe9F5LShMsMS5Zb9qwVlaaVeeRMjjcayfouZdKqs7MT8/Pz6nmWu85yB/Lj+PI++TXlcb29vTXbLceS5dKFxuGynNvc1nLHXEJPTk4WTAnrOVmJPjj82lbCB16HZd7Jfa87txWy3Jr1qBERakFQUaJzB0FFeU9MHPHKckvyWitBcvnnlXuyzVrSBblvjZZLWc7bok0wkW22orWs3GJlrrAtW7bg7bffhoZ3dUqgKNN5+61ca1yGy9F371I6zbgy5XmWN8WaJGR5O/j9kc9Ti1qU21okm2yHRvxyB0FFic4dBBUlOncQVJQVs7mtWTOaDWuFVFo0YlNxeJst+8f6LpoLbWFhwR3eaiVS1K4t7ze/d/n92LJlC44fP67eA8st5HH5AfbSvJyZmRm1TH7O67eSUlihtJz8PZIhwWVoLkt5PzRbWu7z5xk2dxAELqJzB0FFWbHlhKyILC5HtJxYgD+fNT/Pkl2W68eS3lYOds3lsmrVKtV1ZU0OaTZJAoffm3feeaehqL2yOgD9HkhJrT0XGZFlrb5qRQt6V+iUCRqAxXfNyh8v22JFoVmuTe2ZSVnuXUpLI365g6CiROcOgooSnTsIKsqKrRVmTUTX7BVruVZvojxp/3hnQXmT1Vl5xWXYJLdHrdBRrb2yLdYywJqrbXJyUs0zLtthLSWshYtaz0yizXSynq3X5rYSceSsXbsW09PT7jZaodFeV5i17XUXa8QvdxBUlOjcQVBR2ibLrcnmlpvMkuXe6KRcGnd0dJjL51hJGCz3i9Uma0lZLQmDlRBQoklIr4kxNzenJmq05Kks49+tLDFE2b7cbma2l/UsrHxommRvdjkhKze5Jam1KLT5+fnlz6EWBMH7k+jcQVBRVixCzZusgWPl85aUjT739fWZqyVa7bBGXrkk9U5msc6TMpxH6cn0v/z78POsSR/WBBPvqH0rEijIc7xpib2j9to58jiZZtsyTbyy3JLUnmQkIcuDIFCJzh0EFSU6dxBUlLbZ3FbecmmvaFFd1nGWi4vbco3Y3J7c0LJ+C69d6l0PTB7Ly6wINf69urq6Cm3xjhlYkXhliSHKjrNcg1bkoBWZqJVZEWrW0lYS7X3xJMasVyaTlFjutHrtBOKXOwgqi7tzE1EHEf2ciB7J9rcR0VNENEJE3yei7np1BEHQPhqR5bdhcQHAddn+nwP4ekrpfiL6awC3APi2drKUETzyTEpBTdbJ47y5tHh9lvT2SJ0yLAlZ5rLo6upySzVZv/yefNkgfp6U77wOfj+sZAoSywUlJbbWXks283tl5Ym3XFXeZA0a1kQU2S4r4Yi3zNq2kkG0TJYT0WYA/xrAXdk+AbgKwIPZIXsAXOepKwiC9kCev25E9CCA/w6gH8B/BHAzgCdTStuz8i0AfphSuig/Z3x8vFbxyMhIa1sdBAEAYHh4uLY9MDBQkDF1ZTkR/SGAYymlZ4nod/OPSw41/0o88MADhf1169bVtvv6+gplvb29tW0uO+WC8j09PaXHAeXSc3BwECdOnCgc1wpZziWwlX453+7t7T1j+RxtsgxQTG08OTlZKOOra2orbQLlExmuueYa7N27V43Ek3hluRVRV/bMrrjiCjzxxBOF58uPs8y2Vsry4eHhM36IvLLcms8tzQr+PCcmJmrb+bO96aabcO+995aWyfrvvltfOdtjc18J4Boi+gyAHiza3N8AsJ6IOlNK8wA2AzhsVWIlkLPsTy2hAQDVRpP72rbc90zoL8N6kPLlBhY7N++IgL0OF9+33HOWq81K1Mg7j3fGmzfE1JqdJt2cWp3yD6a1RpyGNSbjXUtOllszBZe6Fp61xHPZfhl1be6U0ldSSptTSlsB3ADgn1JKNwF4DMD12WG7ADxc92pBELSNpfi57wDwJ0T0CoCNAHR9EARB22koQi2l9GMAP862XwXwce+5fIlXoGhXS+nK973Sx5r145Xl3plDlnlgRUzlbNy4EePj46pbyMrjbs0O8kbUNYtXAls5zrQZaTMzM+q4g5Wv3lsmZblm7s3PzzeVrMF6N5t1hXmXhtaICLUgqCjRuYOgorRt4oiU5Vx2STcWlyCWW8WSPtrkk1bIcitfliWbeR1jY2NumdtMtJN3KR2ZJIHTyDJDmlkh76m2auv09LQ6Gu91u8k6eZn8jtqyRgsLC02tHNrIaLn2vshlhrx52DTilzsIKkp07iCoKNG5g6CirJjNzSO0ZFipZmdbiQosm8drc3uXpZX2jrWMj2Z7TU1NuW1a75K13mV7JZodad0riTZDzxrH4Lbu1NSUe2YZfw+sqEIrCYUWKWclhpDlzdrcWjRiI64waynqnPjlDoKKEp07CCrKislyvi9nMHGZ7s21Jss0t4qUSF7J7l1x0Uq0IGWXJoG9Mq5Z+L2Zm5tTzY9GTBjN5WeZA/J7aTK6bPJNjjVBRnOHAvbEIm+EmncpKmsikHbc7OysuuwQcGZ/KiN+uYOgokTnDoKKEp07CCpK22xuaVdzm0HaDzwLh5XMryzLSY5mN0nbyLK5vbPCLCxbtBXJ9vg+vwfWemacd999V21jI6G6zayrZo0zeBNlWM/CqkNro3wuEm18xUqm4F3CV25rZWX7ZcQvdxBUlOjcQVBR3pOynO9bstyKYtKQ8smKpvImQvAuEyRlrRbh5J1ZJuts5Nr8Wu1MbMGxZmp5Zbl1P7z1eb8zoD8nKymFJal5v5AyXCsr2y8jfrmDoKJE5w6CitI2WS5HDL2y3JsT3Fq2xmqHJec1WWetcOldTbKzs1OVg7JNvMzrFZCfWyPiVuSZVqfVRgtrlNpKv8zh+80sGSTPk8/BmwPOm57bStbApbfsEzFaHgRBKdG5g6CiROcOgorSNptb2gjclrZsbj5DzErW4LWdZTuaWXfKsu9lGbeR+bV6e3vVMQPLxrRmtVlo96e7u9u0MTlWmZaD3WqfzDHO7wdPminHWvi+lbecY42FeF2UgG5ny/fK68bi9ckEicu+nFAQBO9PXL/cRPQ6gFMAFgDMp5QuJaJBAN8HsBXA6wD+KKU0ujzNDIKgURqR5b+XUnqb7e8G8GhK6U4i2p3t36GdLGUEl0LNynJr4gjHkuWWm0nDShAgJSRvM7+WlOX8PK+Mk3VakWze1Tut6DLvSpkceT+0c3p6etQEDY08dy1JhxXJxllYWHC7uCxXlRWhpuUqb/sqnwbXAtiTbe8BcN0S6gqCoMWQJwiAiF4DMAogAfifKaW/IaKxlNJ6dsxoSmlDvj8+Pl6rWC5oHgRBaxgeHq5tDwwMFCSlV5ZfmVI6TERDAPYR0UuNNuKLX/xiYb+/v7+23dfXVyhbv772NwMDAwOl5wDFlUL5HHCgONqay7orrrgCTz/9dOE4a+RVw5LlVlkuBS+++GIcOHBAldvtkuWf/exn8dBDD6nHWZNDvHgmjtxwww24//77l1WWW16MnE9+8pP4yU9+0hJZfvr06dr2xMREoWxycrK2PTY2Vts+deoUAOCb3/wmbr311oJ5ylOBA78xXZ9//vkzvkeOq3OnlA5n/x8joh9gceneo0S0KaV0hIg2AThm1WFN7reSv/FOKx+ydzYWfynl7DTL5tZeFCsEVJbxPx68HatXr25qVptEyxfeiEtLu7Y1a8vrgrPgdXR3d6sd2noullvSmhWm/RGTCSOtpXO9Nrd857Q/1o3Y3J6c93XfKCJaQ0T9+TaAPwBwEMBeALuyw3YBeLju1YIgaBueX+5zAPwg+4vYCeBvU0r/QEQ/BfAAEd0C4NcAPrd8zQyCoFHqdu6U0qsAdpZ8/g6Aq70XkrJCW8pVHqvNoAH8rjBv3m/L9cPbaElByxXGkctz4DQjAAAJgklEQVQWWy4oS15qubMt21/KWu+MLkuKNzOjiz+z3t5eNTGHZbJ4owotV5hMwLDU/GdA8V2VElrLW27lcpN1eJ5ZRKgFQUWJzh0EFSU6dxBUlLbNCrMSzVk2Cd9uxN+p2VuWzS3r0OxlyxVm2Zjyc2sNLA3LlvauI2a57rwJF73jDrJ+LeRWhp9aWXa8trQ1jsFtZ5kw0psBxWtzW+40LcniwsKCao+X7ZcRv9xBUFGicwdBRXlPyHLpJtMmuvMZYrIOKyzTkuUyaaFWZiVqtLCiujRXmzdHulUm5bAm/7q6uszZUxwrYYUmqa2IPfn9tcQZVjJG73JC1jJM0vXlTcLgdZlZslxbUkku8exdHqpwTN0jgiB4XxKdOwgqSttkuYU1qcSbz8oaLedImSUlqtWuHEu+e0c1LZnVrOz3SmqO975JrNztVnSZJtk7OztVs6KR2Wia3LYiJGXUmTVDT4s8syaOWJLdaq8m34EYLQ+CDzTRuYOgokTnDoKK8p60uTUXhjUzRtpUmgvKcpl57RppmzdjG8m1sTierCFl12smX3hXV5eaINGy67wJK6xEC/I4LfmG1Q7LxdVMdpt6Nrc3WYNlj2tjAdaaZRGhFgRBjejcQVBR2ibLLWnpnRhgSTAph7UJFd4Eg/LaHCvnm+XGklFoWpSRN9GCrLPVstyT3LBs35vjTG5rUryR6EbNBSWfuyapZ2Zm3HJeS7ogryfbr71zUpZbOQJDlgfBB5jo3EFQUaJzB0FFaZvNbc2Iso712l5WckNr9hgvs0IxeZvk2mZa6KWEt3F2dlY9T9pT3mQQjYSt5kgXlNcV5n2e1nHyWtqzthYG8CYwtHKHS9eUd/ld71peVo5xzRUrbW6JJ298/HIHQUWJzh0EFaVtslxKRstdwuWqJT+sCB7NjSNlFr+WJfs5zebR5skmZmdnm8qhZuG9p9YSSs3Kcu/yQtpzkRFZXneXJZu5vJayXJP2s7Oz7mWCrDZ6ltgF/Ik5rFx3GvHLHQQVxdW5iWg9ET1IRC8R0YtEdDkRDRLRPiIayf7fUL+mIAjahVeW/yWAf0gpXU9E3QD6APwpgEdTSncS0W4AuwHcoV7IkOUSTZZbo8hWhBrHO0mlrM0erGg7Xv/U1FRBpnOJ7o3qkniTRngjn5pJ/gDoI+JW/TJnmHfihZTbmhS3VneVx1kRjVobG1mFU5PUjaz66llu2rPK5zoAvwPgbgBIKc2mlMYAXAtgT3bYHgDX1b1aEARtg+rFqBLRJQD+BsAvsLgg4LMAbgPwZkppPTtuNKVUk+bj4+O1ikdGRlrc7CAIAGB4eLi2PTAwUJBhHt3ZCeBjAG5NKT1FRH+JRQneEF/60pcK+5YE6e3trW339PTUtmVqY14mV80sky1f/vKX8a1vfavwmSVRtVUnZTt4e/v6+tSy/LzLLrsMTz75ZMtluRX4USYnL730UjzzzDPqnOJGcpdxGpXln/jEJ/D444+3XJafPn1aPa5Mlt9+++342te+Zq4WorXRCqaR8CCoiYmJ2nb+ve677z7ceOONZq6B/NiDBw+q1/F07kMADqWUnsr2H8Ri5z5KRJtSSkeIaBOAY1Yl3tlSgL2UDMfrfrGwls/hL1WzucS1a50+fbrwkvI/TvJ+WC4RzYb1Lt9k2boSr2vMyrettXF2dla1YS13lDVDz9tJrXzh1h9Ja7kfCy3/u3c2nfd6dW3ulNJbAN4got/OProaixJ9L4Bd2We7ADxc92pBELQN73DwrQDuzUbKXwXwb7H4h+EBIroFwK8BfG55mhgEQTO4OndK6TkAl5YUXe2+UAsi1BpJ+KBhyRsrWUMjiQs0pCuM29lcNjYiyzWsBAHS3rQm4HCatcE52rWmp6dVG9Y7OSSvp6zMK8vn5uZctq6swzsWAthLWPHPG1mRtoyIUAuCihKdOwgqSnTuIKgoK5Yg0Tvs36xbTLNJrOQBVl70VrjCpM3N28hdYd6c4BbWLDm+PTk56fZtN2Nzy3M0N5l0DVoJGazEhN4EhtoMtIWFBTPRgnYfvbnlJTJZJd/25uVX6657RBAE70vqhp82Cw8/DYJg+ZHhp/HLHQQVJTp3EFSUZZPlQRCsLPHLHQQVZdk7NxF9moh+SUSvZBlb2gIRfYeIjhHRQfZZW1NDEdEWInosS031AhHdtkLt6CGip4lof9aOP8s+30ZET2Xt+H42d2DZIaIOIvo5ET2yUu0goteJ6AARPUdEz2SftT112HKmMFvWzk1EHQD+CsC/AnAhgBuJ6MLlvCbjHgCfFp/txmJqqGEAj6KJeekNMg/g9pTSRwBcBuCPs+/f7nbMALgqpbQTwCUAPk1ElwH4cwBfz9oxCuCWZW5Hzm0AXmT7K9WO30spXZJSyudNtPu5AL9JYbYDi8lQXmxZO/KVHpbjH4DLAfwj2/8KgK8s5zXF9bcCOMj2fwlgU7a9CcAv29WW7JoPA/jUSrYDi/nvfgbgXwJ4G0Bn2bNaxutvzl7YqwA8AoBWqB2vAzhLfNbW5wJgHYDXkI19tbodyy3LzwXwBts/lH22UpyTUjoCANn/Q+26MBFtBfBRAE+tRDsyKfwcFpNq7APw/wCMpZTy0Kd2PZtvAPhPAPLwq40r1I4E4EdE9CwR/bvss3Y/l98CcBzAdzMz5S4iWtOqdix35y6Ly/zADc8T0VoAfwfgP6SUTq5EG1JKCymlS7D4y/lxAB8pO2w520BEfwjgWErpWf5xu9uRcWVK6WNYNBn/mIh+pw3XlOQpzL6dUvoogEm00BRY7s59CMAWtr8ZwOFlvqbF0SwlFDypoVoBEXVhsWPfm1J6aKXakZMWM9f+GItjAOuJKJ9f0I5ncyWAa4jodQD3Y1Gaf2MF2oGU0uHs/2MAfoDFP3jtfi5lKcw+1qp2LHfn/imA4Ww0tBvADVhMz7RStDU1FC3O+LgbwIsppb9YwXacTUTrs+1eAL+PxYGbxwBc3652pJS+klLanFLaisV34Z9SSje1ux1EtIaI+vNtAH8A4CDa/FzScqcwa8PAxWcAvIxFG+8/L/f12HXvA3AEwBwW/0LegkX77lEAI9n/g8vchk9gUWI+D+C57N9nVqAd/wLAz7N2HATwX7LPfwvA0wBeAfC/AKxu4/P5XQCPrEQ7suvtz/69kL+X7X4u2TUvAfBM9mz+N4ANrWpHRKgFQUWJCLUgqCjRuYOgokTnDoKKEp07CCpKdO4gqCjRuYOgokTnDoKKEp07CCrK/we4YHuBGfBGvgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"def mean_face(faces):\n",
" return faces.mean(axis=1).reshape((64, 64))\n",
"\n",
"plt.imshow(mean_face(faces), cmap='gray');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Loops in Python are slow, and most of the time you want to utilise the fast native code provided by Numpy without explicitly using\n",
"for loops. 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": 46,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 3.61 ms, sys: 4.09 ms, total: 7.71 ms\n",
"Wall time: 67.5 ms\n",
"CPU times: user 1.02 ms, sys: 41 µs, total: 1.06 ms\n",
"Wall time: 779 µs\n"
]
}
],
"source": [
"# We have some HUUUGE data matrix which we want to compute its mean\n",
"X = np.random.randn(20, 1000)\n",
"# Benchmarking time for computing mean\n",
"%time mean_naive(X)\n",
"%time mean(X)\n",
"pass"
]
},
{
"cell_type": "code",
"execution_count": 47,
"metadata": {
"scrolled": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 2.39 s, sys: 0 ns, total: 2.39 s\n",
"Wall time: 4.81 s\n",
"CPU times: user 34.9 ms, sys: 39.3 ms, total: 74.2 ms\n",
"Wall time: 184 ms\n"
]
}
],
"source": [
"# Benchmarking time for computing covariance\n",
"%time cov_naive(X)\n",
"%time cov(X)\n",
"pass"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As you can see, using Numpy's functions makes the code much faster! Therefore, whenever you can use something that's implemented in Numpy, be sure that you take advantage of that."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## 2. Affine Transformation of Datasets\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 $(D, N)$. We would like to know\n",
"what is the covariance when we apply affine transformation $\\boldsymbol A\\boldsymbol x_i + \\boldsymbol b$ for 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.\n",
"\n",
"For this assignment, you will need to implement the `affine_mean` and `affine_covariance` in the cell below."
]
},
{
"cell_type": "code",
"execution_count": 48,
"metadata": {},
"outputs": [],
"source": [
"# GRADED FUNCTION: DO NOT EDIT THIS LINE\n",
"def affine_mean(mean, A, b):\n",
" \"\"\"Compute the mean after affine transformation\n",
" Args:\n",
" x: ndarray, the mean vector\n",
" A, b: affine transformation applied to x\n",
" Returns:\n",
" mean vector after affine transformation\n",
" \"\"\"\n",
" ### Edit the code below to compute the mean vector after affine transformation\n",
" affine_m = np.zeros(mean.shape) # affine_m has shape (D, 1)\n",
" ### Update affine_m\n",
" \n",
" ###\n",
" return affine_m\n",
"\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",
" ### EDIT the code below to compute the covariance matrix after affine transformation\n",
" affine_cov = np.zeros(S.shape) # affine_cov has shape (D, D)\n",
" ### Update affine_cov\n",
" \n",
" ###\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 $\\boldsymbol A$ and $\\boldsymbol b$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"random = np.random.RandomState(42)\n",
"A = random.randn(4,4)\n",
"b = random.randn(4,1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Next we can generate some random matrix $\\boldsymbol X$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"X = random.randn(4,100) # D = 4, N = 100"
]
},
{
"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": null,
"metadata": {},
"outputs": [],
"source": [
"X1 = (A @ X) + b # applying affine transformation to each sample in X\n",
"X2 = (A @ X1) + b # twice"
]
},
{
"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`.\n",
"\n",
"Check the Numpy [documentation](https://docs.scipy.org/doc/numpy-1.13.0/reference/routines.testing.html)\n",
"for details. The mostly used function is `np.testing.assert_almost_equal`, which raises AssertionError if the two arrays are not almost equal."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"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))"
]
}
],
"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.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment