COMP5318_Assignment_1_GPC_Submission.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "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 choice\n",
    "from numpy import linalg as LA\n",
    "from scipy import optimize\n",
    "from scipy.special import expit\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 from 2D to 1D 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",
    "    label = np.array(label, dtype=np.int32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 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'])\n",
    "    label_test = np.array(label_test, dtype=np.int32)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Singular Value Decomposition "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "k_svd = 10 # number of singular values to save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Singular Value Decomposition on training images\n",
    "\n",
    "u, s, vh = LA.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": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Singular Value Decomposition on testing images\n",
    "\n",
    "u, s, vh = LA.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": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "def PCA(n_components, data):\n",
    "    \"\"\"\n",
    "    Principal Component Analysis\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    n_components: int\n",
    "        Number of principal components to use\n",
    "        \n",
    "    data : array, shape = (n_samples, n_features)\n",
    "        Data on which to perform PCA\n",
    "        \n",
    "    Returns \n",
    "    -------\n",
    "    \n",
    "    X_PCA: array, shape = (n_samples, n_components)\n",
    "        PCA transformed data\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 = LA.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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(30000, 20)\n",
      "(5000, 20)\n"
     ]
    }
   ],
   "source": [
    "n_components = 20\n",
    "X_PCA = 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": [
    "## Gaussian Process Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_sum2(x0, x1):\n",
    "    \"\"\"\n",
    "    Generator for computing individual sum of squares entries. \n",
    "    \"\"\"\n",
    "    for i in range(x0.shape[0]):\n",
    "        for j in range(x1.shape[0]):\n",
    "            yield i, j, np.sum(np.square(x0[i] - x1[j]))\n",
    "\n",
    "def kernel(x0, x1, theta, deriv=False):\n",
    "    \"\"\"\n",
    "    Calculate the kernel matrix. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    x0, x1: arrays, shape = (n_samples, n_features)\n",
    "        input arrays\n",
    "    \n",
    "    theta: array, shape = (3,)\n",
    "        hyperparameters for kernel function\n",
    "    \n",
    "    deriv: boolean\n",
    "        whether to return derivatives of kernel matrix wrt theta\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    \n",
    "    k: array, shape = (x0.shape[0], x1.shape[0])\n",
    "        kernel matrix, optional\n",
    "        returned if deriv == False\n",
    "    \n",
    "    K: array, shape = (len(theta)+1, x0.shape[0], x1.shape[0])\n",
    "        kernel matrix and derivatives wrt theta, optional\n",
    "        returned if deriv == True\n",
    "    \"\"\"\n",
    "    # Take the exponent of theta to ensure positive values\n",
    "    theta = np.exp(theta)\n",
    "    \n",
    "    # Sum of squares for input values\n",
    "    sum2 = np.zeros((x0.shape[0], x1.shape[0]))\n",
    "    for i, j, s in gen_sum2(x0, x1):\n",
    "        sum2[i][j] = s*theta[1]\n",
    "\n",
    "    k = theta[0] * np.exp(-0.5*sum2)\n",
    "    \n",
    "    if deriv:\n",
    "        K = np.zeros((len(theta)+1, x0.shape[0], x1.shape[0]))\n",
    "        # K[:,:,0] is the original covariance matrix\n",
    "        K[0] = k + theta[2]*np.eye(x0.shape[0], x1.shape[0])\n",
    "        K[1] = k\n",
    "        K[2] = -0.5*sum2*k\n",
    "        K[3] = theta[2]*np.eye(x0.shape[0], x1.shape[0])\n",
    "        return K\n",
    "    else:\n",
    "        return k + theta[2]*np.eye(x0.shape[0], x1.shape[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logPosterior(theta, *args):\n",
    "    \"\"\"\n",
    "    Calculate the negative posterior log likelihood. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    \n",
    "    theta: array, shape = (3,)\n",
    "        hyperparameters for kernel function\n",
    "    \n",
    "    data: array, shape = (n_samples, n_features)\n",
    "        input array\n",
    "    \n",
    "    t: array, shape = (1, n_samples)\n",
    "        label column\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    \n",
    "    -logp: float\n",
    "        negative of posterior log likelihood\n",
    "    \n",
    "    \n",
    "    \"\"\"\n",
    "    data, t = args\n",
    "    k = kernel(data, data, theta, deriv=False)\n",
    "    L = LA.cholesky(k)\n",
    "    alpha = LA.solve(L.T, LA.solve(L, t))\n",
    "    logp = -0.5*t.T @ alpha - np.sum(np.log(np.diag(L))) - 0.5*data.shape[0]*np.log(2.0*np.pi)\n",
    "    return -logp\n",
    "\n",
    "def gradLogPosterior(theta, *args):\n",
    "    \"\"\"\n",
    "    Calculate the gradient of the posterior log likelihood. \n",
    "    \n",
    "    theta: array, shape = (3,)\n",
    "        hyperparameters for kernel function\n",
    "    \n",
    "    data: array, shape = (n_samples, n_features)\n",
    "        input array\n",
    "    \n",
    "    t: array, shape = (1, n_samples)\n",
    "        label column\n",
    "    \n",
    "    Returns:\n",
    "    --------\n",
    "    \n",
    "    -dlogpdtheta: array, shape = (3,)\n",
    "        negative of posterior log likelihood gradient\n",
    "    \n",
    "    \"\"\"\n",
    "    data, t = args\n",
    "    K = kernel(data, data, theta, deriv=True)\n",
    "\n",
    "    L = LA.cholesky(K[0])\n",
    "    k_inv = LA.solve(L.T, LA.solve(L, np.eye(data.shape[0])))\n",
    "    dlogpdtheta = np.zeros(len(theta))\n",
    "    for d in range(1, len(theta)+1):\n",
    "        dlogpdtheta[d-1] = 0.5*t.T@k_inv@K[d]@k_inv@t - 0.5*np.trace(k_inv@K[d])\n",
    "\n",
    "    return -dlogpdtheta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training Class 0\n",
      "Counter({0: 9049, 1: 951})\n",
      "Training Class 1\n",
      "Counter({0: 8986, 1: 1014})\n",
      "Training Class 2\n",
      "Counter({0: 8994, 1: 1006})\n",
      "Training Class 3\n",
      "Counter({0: 9019, 1: 981})\n",
      "Training Class 4\n",
      "Counter({0: 8928, 1: 1072})\n",
      "Training Class 5\n",
      "Counter({0: 8998, 1: 1002})\n",
      "Training Class 6\n",
      "Counter({0: 8998, 1: 1002})\n",
      "Training Class 7\n",
      "Counter({0: 9047, 1: 953})\n",
      "Training Class 8\n",
      "Counter({0: 9023, 1: 977})\n",
      "Training Class 9\n",
      "Counter({0: 8958, 1: 1042})\n",
      "0.85\n",
      "Wall time: 3h 44min 20s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "sub = choice(data.shape[0], 10000)\n",
    "train = flat(data)[sub]\n",
    "train_label = label[sub]\n",
    "test = flat(data_test)[:2000]\n",
    "y_true = label_test\n",
    "\n",
    "theta = np.log([1, 0.005, 1e-10])\n",
    "values = np.array(list(set(train_label)))\n",
    "P = []\n",
    "for i in values:\n",
    "    # create binary class labels for class i\n",
    "    t = np.zeros(train_label.shape[0], dtype=int)\n",
    "    t[train_label == i] = 1\n",
    "    t = t.T\n",
    "    print('Training Class {}'.format(i))\n",
    "    print(Counter(t))\n",
    "    \n",
    "    K = kernel(train, train, theta, deriv=False)\n",
    "    K_s = kernel(test, train, theta, deriv=False)\n",
    "    L = LA.cholesky(K)\n",
    "    alpha = LA.solve(L.T, LA.solve(L, t))\n",
    "    mu = K_s @ alpha\n",
    "    logit = expit(mu)\n",
    "    prob = (logit-np.amin(logit))/(np.ptp(expit(mu)))\n",
    "    P.append(prob)\n",
    "P = np.array(P)\n",
    "y_pred = np.argmax(P, axis=0)\n",
    "y_pred = np.array([values[i] for i in y_pred])\n",
    "print(get_accuracy(y_pred, y_true))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('predicted_labels.h5','w') as H:\n",
    "    H.create_dataset('label',data=y_pred)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "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
}