{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tutorial 3 - Dimensionality reduction and clustering\n",
    "\n",
    "In this tutorial, we'll look at techniques of reducing the dimension of data and techniques for clustering data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.decomposition import PCA\n",
    "from sklearn.manifold import TSNE\n",
    "import sklearn\n",
    "from torchvision import datasets, transforms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dimensionality reduction\n",
    "\n",
    "We are focusing on Unsupervised learning, which is a category where the model is trained on a sample of the input data only and generally tries to find structure or patterns in the data. Dimensionality reduction, specifically PCA, falls into this category.\n",
    "\n",
    "### Principal Component Analysis (PCA)\n",
    "\n",
    "First, let's load the Iris dataset, a very simple dataset with 150 samples, 3 classes and 4 features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "iris = load_iris()\n",
    "X = iris.data  # Input features (Sepal length/width, Petal length/width)\n",
    "y = iris.target # Target labels (Setosa, Versicolor, Virginica)\n",
    "\n",
    "print(\"Input shape X:\", X.shape)\n",
    "print(\"Output shape y:\", y.shape) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll create a PCA model using scikit-learn. This model follows the standard fit-transform interface, which we'll use on the whole dataset here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_sklearn = PCA(n_components=2) \n",
    "pca_sklearn.fit(X) \n",
    "\n",
    "X_reduced_sklearn = pca_sklearn.transform(X)\n",
    "\n",
    "print(\"Original data shape:\", X.shape)\n",
    "print(\"Reduced data shape:\", X_reduced_sklearn.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the resulting 2D represenation of the dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_pca(X, y, title, axis = \"Principal Component\"):\n",
    "    plt.figure(figsize=(8, 6))\n",
    "    plt.scatter(X[:, 0], X[:, 1], c = y, cmap = \"viridis\") \n",
    "    plt.title(title)\n",
    "    plt.xlabel(axis + \" 1\")\n",
    "    plt.ylabel(axis + \" 2\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pca(X_reduced_sklearn, y, \"Iris Data projected by PCA\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 1: Manual implementation of PCA\n",
    "\n",
    "The input data $X$ must be handled as a matrix $X \\in \\mathbb{R}^{n \\times d}$, where $n$ is the number of samples and $d$ is the number of features.\n",
    "\n",
    "**Goal:** To find a projection matrix $W$ that transforms the high-dimensional data $X$ into a lower-dimensional space $X_{reduced} \\in \\mathbb{R}^{n \\times k}$ ($k < d$) while retaining the maximum variance.\n",
    "\n",
    "#### 1. Centering the Data\n",
    "\n",
    "PCA is sensitive to the mean of the data. The first step involves centering the data matrix $X$ by subtracting the mean of each feature (column) from all values in that feature.\n",
    "\n",
    "$$\n",
    "\\text{Centering: } X_{\\text{centered}} = X - \\mu_X\n",
    "$$\n",
    "\n",
    "(Where $\\mu_X$ is the vector of means for each feature.) This ensures that the first principal component is not simply representing the mean of the data.\n",
    "\n",
    "#### 2. Calculating the Covariance Matrix ($\\Sigma$)\n",
    "\n",
    "The covariance matrix $\\Sigma$ captures the variance of each feature and the covariance (relationship) between all pairs of features. This step is critical because PCA aims to find directions (components) that maximize variance.\n",
    "\n",
    "For centered data $X_{\\text{centered}}$, the covariance matrix is calculated using matrix multiplication and scaling:\n",
    "\n",
    "$$\n",
    "\\Sigma = \\frac{1}{n-1} X_{\\text{centered}}^T X_{\\text{centered}}\n",
    "$$\n",
    "\n",
    "#### 3. Eigendecomposition\n",
    "\n",
    "This step finds the fundamental directions of variance in the data. Eigendecomposition is applied to the covariance matrix $\\Sigma$.\n",
    "\n",
    "The process finds the **eigenvectors** and **eigenvalues** of $\\Sigma$:\n",
    "\n",
    "$$\n",
    "\\Sigma \\vec{v} = \\lambda \\vec{v}\n",
    "$$\n",
    "\n",
    "*   The **eigenvectors** ($\\vec{v}$) of $\\Sigma$ are the **Principal Components (PCs)**. They represent the new, orthogonal axes (directions of variance).\n",
    "*   The **eigenvalues** ($\\lambda$) are scalar values that represent the magnitude of variance captured along their corresponding eigenvector (PC).\n",
    "\n",
    "#### 4. Selecting Components and Projection\n",
    "\n",
    "1.  **Sort:** The eigenvectors (PCs) are sorted in descending order based on their corresponding eigenvalues. The components with the largest eigenvalues capture the most variance and are thus the most important.\n",
    "2.  **Select $k$:** A user selects the top $k$ principal components (e.g., $k=2$ for visualization). These $k$ eigenvectors form the **projection matrix** $W$.\n",
    "3.  **Project:** The original centered data is multiplied by the projection matrix $W$ to reduce its dimensionality:\n",
    "\n",
    "$$\n",
    "X_{\\text{reduced}} = X_{\\text{centered}} W\n",
    "$$\n",
    "\n",
    "The result, $X_{\\text{reduced}}$, is the data projected onto the lower-dimensional subspace, ready for data visualization or further analysis. This projection step is a matrix multiplication."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Original data shape:\", X.shape)\n",
    "print(\"Projected data shape (Manual):\", X_reduced.shape)\n",
    "plot_pca(X_reduced, y, \"Manual PCA on Iris\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### t-SNE: An alternative to PCA\n",
    "\n",
    "t-SNE is a much newer and a much more complicated alternative to PCA.\n",
    "\n",
    "The t-SNE algorithm comprises two main stages. First, t-SNE constructs a probability distribution over pairs of high-dimensional objects in such a way that similar objects are assigned a higher probability while dissimilar points are assigned a lower probability. Second, t-SNE defines a similar probability distribution over the points in the low-dimensional map, and it minimizes the Kullback–Leibler divergence (KL divergence) between the two distributions with respect to the locations of the points in the map. (Wikipedia)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_reduced_tsne = TSNE().fit_transform(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Original data shape:\", X.shape)\n",
    "print(\"Projected data shape:\", X_reduced_tsne.shape)\n",
    "plot_pca(X_reduced_tsne, y, \"t-SNE on Iris\", \"t-SNE\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### More complex datasets: MNIST\n",
    "\n",
    "Let's also compare PCA and t-SNE on a much more complex dataset: MNIST."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mnist = datasets.MNIST(\n",
    "    root = \"data\",\n",
    "    download = True,\n",
    "    transform = transforms.Compose([transforms.ToTensor(), transforms.Normalize(0, 1)])\n",
    ")\n",
    "mnist_flat = np.reshape(mnist.data, (mnist.data.shape[0], mnist.data.shape[1] * mnist.data.shape[2]))[:10000]\n",
    "\n",
    "plt.figure(figsize=(10, 7))\n",
    "for i in range(24):\n",
    "    plt.subplot(4, 6, i + 1)\n",
    "    plt.xticks([])\n",
    "    plt.yticks([])\n",
    "    plt.grid(False)\n",
    "    plt.imshow(mnist[i][0].reshape(28, 28), cmap = plt.cm.bone)\n",
    "    plt.title(mnist[i][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pca_mnist = PCA(n_components=2) \n",
    "mnist_reduced_pca = pca_mnist.fit_transform(mnist_flat) \n",
    "print(\"Original data shape:\", mnist_flat.shape)\n",
    "print(\"Reduced data shape:\", mnist_reduced_pca.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pca(mnist_reduced_pca, mnist.targets[:10000], \"PCA on MNIST\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given that MNIST is an image dataset, we can also visualize the individual principal components:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_components = 10\n",
    "pca_mnist = PCA(n_components).fit(mnist_flat)\n",
    "fig = plt.figure(figsize=(n_components*np.power(n_components, 1/5),3))\n",
    "axes = []\n",
    "for i in range(n_components):\n",
    "    ax_ = fig.add_subplot(1, n_components, i+1)\n",
    "    ax_.imshow(pca_mnist.components_[i].reshape(28,28), interpolation='nearest', clim=(-.15, .15));\n",
    "    ax_.set_title(f'PC {i+1}')\n",
    "    ax_.axis('off')\n",
    "    axes.append(ax_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mnist_reduced_tsne = TSNE().fit_transform(mnist_flat)\n",
    "print(\"Original data shape:\", mnist_flat.shape)\n",
    "print(\"Reduced data shape:\", mnist_reduced_tsne.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pca(mnist_reduced_tsne, mnist.targets[:10000], \"t-SNE on MNIST\", \"t-SNE\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## K-means clustering\n",
    "\n",
    "K-means Clustering is a partitioning method used in unsupervised learning to divide $n$ data points into $K$ distinct, non-overlapping subsets (clusters). The fundamental goal is to find structure or patterns in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_kmeans = KMeans(n_clusters = 3, random_state = 42) \n",
    "model_kmeans.fit(X) \n",
    "cluster_labels = model_kmeans.predict(X)\n",
    "centroids = model_kmeans.cluster_centers_\n",
    "print(\"Centroids shape:\", centroids.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot, side-by-side, the true labels and the assigned clusters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kmeans(X, y, cluster_labels, centroids):\n",
    "    pca = PCA(n_components=2).fit(X)\n",
    "    X_reduced = pca.transform(X)\n",
    "    centroids_reduced = pca.transform(centroids)\n",
    "    \n",
    "    fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n",
    "\n",
    "    axes[0].scatter(X_reduced[:, 0], X_reduced[:, 1], c = y, cmap = \"viridis\")\n",
    "    axes[0].set_title(\"True labels\")\n",
    "    axes[0].set_xlabel(\"Principal Component 1\")\n",
    "    axes[0].set_ylabel(\"Principal Component 2\")\n",
    "\n",
    "    axes[1].scatter(X_reduced[:, 0], X_reduced[:, 1], c = cluster_labels, cmap = \"viridis\")\n",
    "    axes[1].set_title(\"K-means clusters and centroids\")\n",
    "    axes[1].set_xlabel(\"Principal Component 1\")\n",
    "    axes[1].set_ylabel(\"Principal Component 2\")\n",
    "\n",
    "    axes[1].scatter(centroids_reduced[:, 0], centroids_reduced[:, 1], marker='X', s=200, color='red', edgecolor='black')\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kmeans(X, y, cluster_labels, centroids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Task 2: Manual implementation of K-means\n",
    "\n",
    "The Core Concept: The algorithm operates iteratively to minimize the inertia, which is the sum of squared distances between each point and the centroid of its assigned cluster.\n",
    "\n",
    "The K-means Iterative Process:\n",
    "\n",
    "1. **Initialization:** $K$ initial centroids ($\\mu_1, \\mu_2, \\dots, \\mu_k$) are chosen randomly from the data (or via a smart selection method like k-means++).\n",
    "2. **Assignment Step:** Each data point is assigned to the nearest centroid.\n",
    "3. **Update Step:** The centroids are recalculated as the mean of all data points assigned to that cluster. This new mean becomes the new centroid.\n",
    "4. **Convergence:** The process repeats (steps 2 and 3) until the cluster assignments no longer change, or a maximum number of iterations is reached.\n",
    "\n",
    "\n",
    "Reimplementing K-means manually requires an iterative loop, distance calculations, and recalculation of means, all leveraging NumPy. This follows the spirit of the Manual optimisation of linear regression, which showed the capability of NumPy to handle complex, iterative mathematical tasks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "np.random.seed(42)\n",
    "max_iterations = 100\n",
    "\n",
    "n_samples, n_features = X.shape\n",
    "\n",
    "random_indices = np.random.choice(n_samples, k, replace=False)\n",
    "centroids = X[random_indices]\n",
    "cluster_labels = np.zeros(n_samples, dtype=int)\n",
    "\n",
    "for i in range(max_iterations):\n",
    "    old_centroids = centroids.copy()\n",
    "        \n",
    "    distances = np.sum((X[:, np.newaxis] - centroids)**2, axis=2)\n",
    "    # TODO\n",
    "    # cluster_labels = ...\n",
    "    # centroids = ...\n",
    "\n",
    "    if np.allclose(old_centroids, centroids):\n",
    "       break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kmeans(X, y, cluster_labels, centroids)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nsn",
   "language": "python",
   "name": "nsn"
  },
  "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.11.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}