Homework 1 - PCA.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5d5e7fb7-bd74-4b73-899b-8a58c7dfaf94",
   "metadata": {},
   "source": [
    "<h1>Homework 1: Principal component analysis.</h1>\n",
    "\n",
    "This notebook introduces the useful concept of principal component analysis, which is also a good application of linear algebra techniques. This notebook introduces the useful concept of principal component analysis, which is also a good application of linear algebra techniques. \n",
    "\n",
    "The flow is guided by [this Wikipedia description of how PCA works](https://en.wikipedia.org/wiki/Principal_component_analysis#Computation_using_the_covariance_method), but using SVD instead of eigenvectors of the covariance matrix, [as described here](https://stats.stackexchange.com/questions/134282/relationship-between-svd-and-pca-how-to-use-svd-to-perform-pca).\n",
    "\n",
    "The data is from the dome imaging project, dated July 18, 2023. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "695ea4a9-6446-4a3e-8f6b-82ad439190cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "using FileIO, Images, LinearAlgebra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30779ef9-8d86-4073-b827-34e4af8060ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The first step is to import the files we want to work with. \n",
    "\n",
    "imgDir = raw\"C:\\Users\\anonc\\Downloads\\PCA_imgs\"    # Change this to the location of the image directory in your system.\n",
    "cd(imgDir)    # Change directories into the one with the images.\n",
    "fnames = readdir()    # This gets all filenames in the working directory\n",
    "imgs = load.(fnames);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58923180-7006-446a-91e6-816e97252bba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look through some of the images to get a feel for what they are.\n",
    "imgs[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c3b0f347-e8b7-4539-90e0-ac32e9745ee0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We next want to convert all the images into numeric vectors. We will break this into two steps.  \n",
    "# First we convert to numeric arrays.  \n",
    "matrices = [convert.(Float64,img) for img in imgs]\n",
    "# Check the type of matrices[1] and the type of imgs[1]. Do you understand the difference?\n",
    "# Next we flatten the matrices to vectors. \n",
    "vs = [m[:] for m in matrices]\n",
    "size(vs[1]), size(imgs[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8faf774-d8ce-4aee-9ddd-6662b870351e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Next we want to combine all the vectors vs into one data matrix of size (length(vs),length(vs[1])), i.e. (32,441100).\n",
    "dat = # Your code here\n",
    "size(dat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d2c7bbb-5b63-45bb-bec7-c3b495810170",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the mean along each column of dat and subtract it from dat.\n",
    "# WARNING: This cell mutates dat, so you should run the previous cell before running this one multiple times. \n",
    "datMeans = # Your code here\n",
    "dat .-= datMeans\n",
    "size(datMeans)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ddb4a921-3e5d-4787-97e8-161b46156d24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute the SVD:\n",
    "U,S,Vt = # Your code here\n",
    "S'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3e09273-d3eb-4cb2-9a25-5cfcd8e71ed4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The columns of Vt are the principal components.  Reshape them to have the same shape as the \n",
    "# original images and normalize them to have maximum value 1. \n",
    "pcs = # Your code here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4ce283f8-d69c-4b37-af93-f9a8aa748ad2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Can you interpret the first three principal components? \n",
    "Gray.(pcs[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd5958c6-1099-4855-808d-bdd90bf5f707",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gray.(pcs[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cb48599b-fe5f-43e9-b54d-a0344546be95",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gray.(pcs[3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7c716c2-2ca9-45bc-bc5a-ffc7a2eedca2",
   "metadata": {},
   "outputs": [],
   "source": [
    "Gray.(pcs[4])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0979f32b-ab36-493d-ad82-5d60d08382fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Heatmaps are better in some ways for viewing the principal components, as the values are both positive and negative. \n",
    "using Plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "669fa4c6-25ed-41a2-9e03-bd517981f29e",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(pcs[1],yflip = true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1892645a-a1e7-49c2-b76d-741e94141416",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(pcs[2],yflip = true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5912654e-ba47-4ca4-96d6-67abf6fc0e83",
   "metadata": {},
   "outputs": [],
   "source": [
    "heatmap(pcs[3],yflip = true)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "01e4a918-69db-4709-80ba-c2a3183e1858",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.3",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}