COMP5318_Assignment_1_SVD_PCA.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\jouse\\Anaconda2\\envs\\py36\\lib\\site-packages\\h5py\\__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from `float` to `np.floating` is deprecated. In future, it will be treated as `np.float64 == np.dtype(float).type`.\n",
      "  from ._conv import register_converters as _register_converters\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import h5py\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import scipy\n",
    "from collections import Counter\n",
    "from ipywidgets import interact, interactive, fixed\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "from numpy.random import shuffle\n",
    "from scipy.spatial.distance import pdist, squareform\n",
    "from scipy.spatial import KDTree\n",
    "\n",
    "mpl.style.use('seaborn')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_accuracy(y_true, y_pred): \n",
    "    \"\"\"\n",
    "    Calculate the accuracy score.\n",
    "    \"\"\"\n",
    "    y_true = np.array(y_true)\n",
    "    y_pred = np.array(y_pred)\n",
    "    return np.sum(y_true == y_pred)/len(y_true)\n",
    "\n",
    "def flat(data):\n",
    "    \"\"\"\n",
    "    Flatten images in the data array. \n",
    "    \"\"\"\n",
    "    return data.reshape(data.shape[0], data.shape[1]**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load training data and labels\n",
    "with h5py.File('images_training.h5','r') as H:\n",
    "    data = np.copy(H['data'])\n",
    "    data = np.array([x/np.amax(x) for x in data]) # normalize images\n",
    "with h5py.File('labels_training.h5','r') as H:\n",
    "    label = np.copy(H['label'])\n",
    "    \n",
    "# load testing data and labels\n",
    "with h5py.File('images_testing.h5','r') as H:\n",
    "    data_test = np.copy(H['data'])\n",
    "    data_test = np.array([x/np.amax(x) for x in data_test]) # normalize images\n",
    "with h5py.File('labels_testing_2000.h5','r') as H:\n",
    "    label_test = np.copy(H['label'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_svd = 10 # number of singular values to save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Singular Value Decomposition on training images\n",
    "\n",
    "u, s, vh = np.linalg.svd(data)\n",
    "data_svd = np.array([u[i][:,:k_svd] @ np.diag(s[i][:k_svd]) @ vh[i][:k_svd:,] for i in range(s.shape[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Singular Value Decomposition on testing images\n",
    "\n",
    "u, s, vh = np.linalg.svd(data_test)\n",
    "data_test_svd = np.array([u[i][:,:k_svd] @ np.diag(s[i][:k_svd]) @ vh[i][:k_svd:,] for i in range(s.shape[0])])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Principal Component Analysis\n",
    "\n",
    "Note: PCA is applied to the entire dataset, while SVD is applied to each image. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PCA(n_components, data):\n",
    "    \"\"\"\n",
    "    Principal Component Analysis\n",
    "    \n",
    "    Parameters: \n",
    "    \n",
    "    n_components: int\n",
    "        Number of principal components to use\n",
    "        \n",
    "    data : array-like, shape = (n_samples, n_features)\n",
    "        Data to perform PCA on\n",
    "    \n",
    "    \"\"\"\n",
    "    n_components = n_components\n",
    "    data = data\n",
    "    X = data - np.mean(data, axis=0) # data matrix normalized by mean\n",
    "    S = np.cov(X.T) # compute covariance matrix of X\n",
    "    L, V = np.linalg.eig(S) # L - array of eigenvalues; V - matrix of eigenvectors\n",
    "    V = V[:,np.argsort(-L)][:,:n_components] # sort eigenvectors by descending order on eigenvalues\n",
    "    L = -np.sort(-L)[:n_components] # sort eigenvalues by descending order\n",
    "    X_PCA = V.T @ X.T\n",
    "    X_PCA = X_PCA.T\n",
    "    return X_PCA, V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(30000, 20)\n",
      "(5000, 20)\n"
     ]
    }
   ],
   "source": [
    "n_components = 20 # number of principal components to use\n",
    "X_PCA, V = PCA(n_components=n_components, data=np.append(flat(data), flat(data_test), axis=0))\n",
    "data_pca = X_PCA[:data.shape[0]]\n",
    "data_test_pca = X_PCA[-data_test.shape[0]:]\n",
    "print(data_pca.shape)\n",
    "print(data_test_pca.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-fold Cross Validation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kfold(N, k=10):\n",
    "    \"\"\"\n",
    "    Generate lists of indices for data divided into k equal parts for cross validation. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    N: int\n",
    "        number of samples in the data. \n",
    "    \n",
    "    k: int\n",
    "        number of parts the data will be split into. \n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    indices: array, shape(k,) \n",
    "        An array of arrays, each containing the indices for one part of data. \n",
    "    \"\"\"\n",
    "    \n",
    "    arr = np.arange(N)\n",
    "    np.random.shuffle(arr)\n",
    "    indices = np.array([arr[i::k] for i in range(k)])\n",
    "    \n",
    "    return indices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}