exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "70877d70",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Introduction to Machine Learning"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c30f301e",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "JWC5tKyP3x1e"
   },
   "source": [
    "Written by Morgan Schwartz and David Van Valen.\n",
    "\n",
    "---\n",
    "\n",
    "In this exercise, we are going to follow the basic workflow that is the foundation of any machine or deep learning project\n",
    "1. Data wrangling\n",
    "2. Model configuration and training\n",
    "3. Model evaluation\n",
    "\n",
    "Along the way, we will implement a linear classifier, test a random forest classifier and explore the role of feature engineering in traditional machine learning.\n",
    "\n",
    "We are going to look at a collection of images of Jurkat cells published in the Broad Bioimage Collection ([BBBC048](https://bbbc.broadinstitute.org/BBBC048)). The cells were fixed and stained with PI (propidium iodide) to quantify DNA content and a MPM2 (mitotic protein monoclonal #2) antibody to identify mitotic cells. The goal is to predict the stage of the cell cycle from images like those shown below.\n",
    "\n",
    "![](https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41467-017-00623-3/MediaObjects/41467_2017_623_Fig2_HTML.jpg?as=webp)\n",
    "\n",
    "[Eulenberg et al. (2017) Reconstructing cell cycle and disease progression using deep learning. Nature communications](https://www.nature.com/articles/s41467-017-00623-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a483d8e",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-danger\">\n",
    "Set your python kernel to <code>01_intro_ml</code>\n",
    "\n",
    "![](images/kernel-change.png)\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63d9b674",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Part A: The Linear Classifier\n",
    "\n",
    "While deep learning might seem intimidating, don't worry. Its conceptual underpinnings are rooted in linear algebra and calculus - if you can perform matrix multiplication and take derivatives you can understand what is happening in a deep learning workflow. In this section, we will implement a simple linear classifier by hand and train it to predict cell cycle stages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "534b1a2d",
   "metadata": {
    "id": "_DwKG_Gi3x1f"
   },
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "import os\n",
    "\n",
    "import imageio as iio\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import skimage\n",
    "import sklearn\n",
    "import sklearn.model_selection\n",
    "import sklearn.ensemble\n",
    "import tqdm.auto"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a767e5b",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "rZ2jTjsd3x1g"
   },
   "source": [
    "## The supervised machine learning workflow\n",
    "Recall from class the conceptual workflow for a supervised machine learning project.\n",
    "- First, we create a <em>training dataset</em>, a paired collection of raw data and labels where the labels contain information about the \"insight\" we wish to extract from the raw data.\n",
    "- Once we have training data, we can then use it to train a <em>model</em>. The model is a mathematical black box - it takes in data and transforms it into an output. The model has some parameters that we can adjust to change how it performs this mapping.\n",
    "- Adjusting these parameters to produce outputs that we want is called training the model. To do this we need two things. First, we need a notion of what we want the output to look like. This notion is captured by a <em>loss function</em>, which compares model outputs and labels and produces a score telling us if the model did a \"good\" job or not on our given task. By convention, low values of the loss function's output (e.g. the loss) correspond to good performance and high values to bad performance. We also need an <em>optimization algorithm</em>, which is a set of rules for how to adjust the model parameters to reduce the loss\n",
    "- Using the training data, loss function, and optimization algorithm, we can then train the model\n",
    "- Once the model is trained, we need to evaluate its performance to see how well it performs and what kinds of mistakes it makes. We can also perform this kind of monitoring during training (this is actually a standard practice).\n",
    "\n",
    "Because this workflow defines the lifecycle of most machine learning projects, this notebook is structured to go over each of these steps while constructing a linear classifier."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba93b37b",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "rbIDLvJ23x1g"
   },
   "source": [
    "## Create training data\n",
    "\n",
    "During the initial setup of this exercise, we downloaded the data and unzipped the relevant files using the script `setup.sh`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca4a2761",
   "metadata": {},
   "outputs": [],
   "source": [
    "data_dir = \"data/CellCycle\"\n",
    "sorted(os.listdir(data_dir))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22b681d0",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "The command above should generate the following output. If you see something different, please check that the `setup.sh` script ran correctly.\n",
    "```\n",
    "['66.lst~',\n",
    " 'Anaphase',\n",
    " 'G1',\n",
    " 'G2',\n",
    " 'Metaphase',\n",
    " 'Prophase',\n",
    " 'S',\n",
    " 'Telophase',\n",
    " 'img.lst',\n",
    " 'img.lst~']\n",
    "```\n",
    "\n",
    "The metadata for each file is stored in `img.lst` so we will first load this information to inform how we load the rest of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "190d354d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the csv file using pandas\n",
    "df = pd.read_csv(os.path.join(data_dir, \"img.lst\"), sep=\"\\t\", header=None)\n",
    "\n",
    "# Rename columns to make the data easier to work with\n",
    "df = df.rename(columns={1: \"class\", 2: \"filepath\"})\n",
    "\n",
    "# Extract the channel information from the filepath column and create a new column with it\n",
    "df[\"channel\"] = (\n",
    "    df[\"filepath\"]\n",
    "    .str.split(\"/\", expand=True)[2]\n",
    "    .str.split(\"_\", expand=True)[1]\n",
    "    .str.slice(2, 3)\n",
    ")\n",
    "\n",
    "# Extract the file ID from the filepath and save in its own column\n",
    "df[\"id\"] = df[\"filepath\"].str.split(\"/\", expand=True)[2].str.split(\"_\", expand=True)[0]\n",
    "\n",
    "# Look at the first few rows in the dataframe\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "215d0434",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check the total number of unique classes in the dataset\n",
    "df[\"class\"].unique()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "068c3a17",
   "metadata": {},
   "outputs": [],
   "source": [
    "class_lut = [\"Ana\", \"Meta\", \"Pro\", \"Telo\", \"G1\", \"G2\", \"S\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6793598",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "For each `id` there are three images. One for each of the channels: phase, PI and MPM2. Let's take a look at a single image"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9d69e796",
   "metadata": {},
   "outputs": [],
   "source": [
    "im_id = \"12432\"\n",
    "\n",
    "# Load channel 3\n",
    "filepath = df[(df[\"id\"] == im_id) & (df[\"channel\"] == \"3\")][\"filepath\"].values[0]\n",
    "im3 = iio.imread(os.path.join(data_dir, filepath))\n",
    "\n",
    "# Load channel 4\n",
    "filepath = df[(df[\"id\"] == im_id) & (df[\"channel\"] == \"4\")][\"filepath\"].values[0]\n",
    "im4 = iio.imread(os.path.join(data_dir, filepath))\n",
    "\n",
    "# Load channel 6\n",
    "filepath = df[(df[\"id\"] == im_id) & (df[\"channel\"] == \"6\")][\"filepath\"].values[0]\n",
    "im6 = iio.imread(os.path.join(data_dir, filepath))\n",
    "\n",
    "# Create a matplotlib subplot with one row and 3 columns\n",
    "# Plot each of the three images\n",
    "fig, ax = plt.subplots(1, 3, figsize=(10, 3))\n",
    "ax[0].imshow(im3, cmap=\"Greys_r\")\n",
    "ax[0].set_title(\"phase\")\n",
    "ax[1].imshow(im4, cmap=\"Greys_r\")\n",
    "ax[1].set_title(\"PI\")\n",
    "ax[2].imshow(im6, cmap=\"Greys_r\")\n",
    "ax[2].set_title(\"MPM2\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d17f255",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Now we can load all of the images into a dataset. We will want to load each of the three channels for each image and create an array with the shape (w, h, ch). Then we will combine all images in the dataset into a single array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd9cc41c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create empty lists to hold images and classes as we load them\n",
    "ims = []\n",
    "ys = []\n",
    "\n",
    "# Iterate over each unique id in the dataset\n",
    "for i, g in df.groupby(\"id\"):\n",
    "    im = []\n",
    "\n",
    "    # Each row in the group corresponds to a different channel for a single image/id\n",
    "    for _, r in g.iterrows():\n",
    "        # Set the complete filepath for this image\n",
    "        path = os.path.join(data_dir, r[\"filepath\"])\n",
    "\n",
    "        # Load in the data using imageio and append to a list\n",
    "        im.append(iio.imread(path))\n",
    "\n",
    "    # Stack a list of three images into a single\n",
    "    im = np.stack(im, axis=-1)\n",
    "    ims.append(im)\n",
    "    ys.append(r[\"class\"])\n",
    "\n",
    "X_data = np.stack(ims)\n",
    "y_data = np.stack(ys)\n",
    "print(\"X shape:\", X_data.shape)\n",
    "print(\"y shape:\", y_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e64e5600",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "WG4mmJ173x1h"
   },
   "source": [
    "In the previous cell, you probably observed that there are 4 dimensions rather than the 3 you might have been expecting. This is because while each image is (66, 66, 3), the full dataset has many images. The different images are stacked along the first dimension. The full size of the training images is (# images, 66, 66, 3)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a3a2a97",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Let's take a look at a sample image from each class and plot each of the three channels separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "47e5c8fb",
   "metadata": {
    "id": "5RIxP_hq3x1i"
   },
   "outputs": [],
   "source": [
    "# Iterate over each class in the dataset\n",
    "for c in np.unique(y_data):\n",
    "    # Select a random index for the class of interest\n",
    "    i = np.random.choice(np.where(y_data == c)[0])\n",
    "\n",
    "    # Create a matplotlib subplot with one row and 3 columsn\n",
    "    fig, ax = plt.subplots(1, 3, figsize=(10, 3))\n",
    "    ax[0].set_ylabel(class_lut[c])\n",
    "\n",
    "    # Plot each of the three channels\n",
    "    for j, ch in enumerate([\"phase contrast\", \"PI\", \"MPM2\"]):\n",
    "        ax[j].imshow(X_data[i, ..., j], cmap=\"Greys_r\")\n",
    "        ax[j].set_title(ch)\n",
    "        ax[j].xaxis.set_tick_params(labelbottom=False)\n",
    "        ax[j].yaxis.set_tick_params(labelleft=False)\n",
    "        ax[j].set_xticks([])\n",
    "        ax[j].set_yticks([])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eeae968",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "2D_YhVir3x1i"
   },
   "source": [
    "For this exercise, we will want to flatten the training data into a vector and select a single channel to work with. We work with the phase channel first."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af92cb98",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Record the original width so that we can use this for reshaping later\n",
    "image_width = X_data.shape[1]\n",
    "image_width"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1100cb72",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "J3MJPnh-3x1j",
    "outputId": "c280504b-94c4-4888-dd02-875c789dee9f"
   },
   "outputs": [],
   "source": [
    "# Flatten the images 1d vectors\n",
    "X_flat = np.reshape(X_data[..., 0], (-1, image_width * image_width))\n",
    "print(X_flat.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "194679c6",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "### Checking Class Balance"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0c55414",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 1.1\n",
    "\n",
    "Let's check the balance of classes in this dataset (stored in `y_data`). There are at least three ways you could do this. Pick one to try.\n",
    "\n",
    "- Count the number of items in each class using `np.unique` ([see docs](https://numpy.org/doc/stable/reference/generated/numpy.unique.html)).\n",
    "- Use the `Counter` object which is imported from `collections` ([see docs](https://docs.python.org/3/library/collections.html#collections.Counter))\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "619591a5",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "# Add your code to check class balances here\n",
    "# You should end up with a count of number of items in each of the 7 classes"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7130d69c",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "This dataset is highly inbalanced so we will want to correct the class balance before training."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bf8b3f6",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "l2yrGjOL3x1j"
   },
   "source": [
    "### Split the training dataset into training and testing datasets\n",
    "How do we know how well our model is doing? A common practice to evaluate models is to evaluate them on splits of the original training dataset. Splitting the data is important, because we want to see how models perform on data that was not used to train them. We split into\n",
    "- The <em>training</em> dataset used to train the model\n",
    "- A held out <em>testing</em> dataset used to evaluate the final trained version of the model\n",
    "\n",
    "While there is no hard and fast rule, 80%/20% splits are a reasonable starting point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d65e087",
   "metadata": {
    "id": "hLQUiSoj3x1j"
   },
   "outputs": [],
   "source": [
    "# Split the dataset into training, validation, and testing splits\n",
    "seed = 10\n",
    "train_size = 0.8\n",
    "X_train, X_test, y_train, y_test = sklearn.model_selection.train_test_split(\n",
    "    X_flat, y_data, train_size=train_size, random_state=seed\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d91a6f6",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "### Correct class imbalance\n",
    "\n",
    "![](https://miro.medium.com/v2/resize:fit:1400/format:webp/1*6DlvkNY2TlFw2Veyvb4qdQ.jpeg)\n",
    "\n",
    "(*Image by [Angelica Lo Duca](https://towardsdatascience.com/how-to-balance-a-dataset-in-python-36dff9d12704)*)\n",
    "\n",
    "There are several ways to correct class imbalance. In this example, we are going to oversample underrepresented classes until we have an equal number of samples for each class.\n",
    "\n",
    "It's important to note that we need to correct class imbalance after generating the train/test split in our dataset. When we are oversampling, we want to prevent samples that are used in our training dataset from also appearing in our testing dataset."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "929cbb5a",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 1.2\n",
    "\n",
    "Complete the `balance_classes` function following the structure outlined in the comments. You can use `sklearn.utils.resample` to generate a new random set of `n_samples`.\n",
    "    \n",
    "Hint: You may want to use boolean indexing to select subsets of arrays. For example we can select all samples in class 3 with the following:\n",
    "```python\n",
    "xx = X_data[y_data == 3]\n",
    "yy = y_data[y_data == 3]\n",
    "```\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4bf212b-79d3-4111-b6f8-488900738efd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "sklearn.utils.resample?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cc603792",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "def balance_classes(X, y):\n",
    "    \"\"\"For a given multiclass dataset, upsample underrepresented classes\n",
    "    to match the number of samples in the majority class\n",
    "    \n",
    "    Args:\n",
    "        X (np.array): Array of raw data\n",
    "        y (np.array): Array of class labels\n",
    "        \n",
    "    Returns:\n",
    "        np.array: X\n",
    "        np.array: y\n",
    "    \"\"\"\n",
    "    classes = np.unique(y)\n",
    "    \n",
    "    # Identify which class has the most samples\n",
    "    # Hint: use your code for counting the number of samples in each class\n",
    "    maj_samples = ...\n",
    "    maj_id = ...\n",
    "    \n",
    "    # Collect new samples in an array\n",
    "    new_X, new_y = [], []\n",
    "    for c in classes:\n",
    "        # Resample the minority classes to match majority number of samples using sklearn.utils.resample\n",
    "        \n",
    "        # Store the new samples in new_X and new_y\n",
    "        ...\n",
    "    # Concatenate the list of arrays to create a single array\n",
    "    new_X = np.concatenate(new_X)\n",
    "    new_y = np.concatenate(new_y)\n",
    "    \n",
    "    # Shuffle arrays to randomize sample order\n",
    "    new_X, new_y = sklearn.utils.shuffle(new_X, new_y)\n",
    "\n",
    "    return new_X, new_y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0825ecc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, y_train = balance_classes(X_train, y_train)\n",
    "X_test, y_test = balance_classes(X_test, y_test)\n",
    "print(f\"Train shape: X {X_train.shape}, y {y_train.shape}\")\n",
    "print(f\"Test shape: X {X_test.shape}, y {y_test.shape}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "801fd1b4-7804-404f-b7e3-6d252fa91e9e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Run this cell to check your class balance code\n",
    "def check_balance(y):\n",
    "    cls, n = np.unique(y, return_counts=True)\n",
    "    max_n = max(n)\n",
    "    if not np.all(n == max_n):\n",
    "        raise ValueError('Sample is not balanced!')\n",
    "        \n",
    "check_balance(y_train)\n",
    "check_balance(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "81836f7d",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "### One-hot encoding\n",
    "\n",
    "Currently, our data have labels that range from 0 to 6. While we know that each of these 7 classes is comparable, this encoding implies that some classes have more weight than others. Alternatively, we want to use a binary encoding so that all classes are seen as equivalent by the model.\n",
    "\n",
    "Instead of representing each label with a number from 0 to 6, we will use an array of length 7 where each position in the array is a binary value encoding the class.\n",
    "\n",
    "For example, `5` is encoded as `[0, 0, 0, 0, 1, 0, 0]`"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a4fb005",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 1.3\n",
    "\n",
    "In order to transform our data from integer to one-hot encoding, we will use `sklearn.preprocessing.LabelBinarizer` ([docs](https://scikit-learn.org/stable/modules/generated/sklearn.preprocessing.LabelBinarizer.html#sklearn.preprocessing.LabelBinarizer)). Take a look at the documentation to learn how to initialize and fit the `LabelBinarizer` and add your code below.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d35f31bc",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "# Initialize and fit the LabelBinarizer\n",
    "lb = ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "19210cb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "lb.classes_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c23fba76",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Run the following cell to check that you set up the `LabelBinarizer` correctly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78756745",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Check that the detected classes are correct\n",
    "print(lb.classes_, \"\\n\")\n",
    "assert np.all(lb.classes_ == np.arange(7))\n",
    "\n",
    "# Test a transformation\n",
    "sample = lb.transform([1, 4])\n",
    "print(sample)\n",
    "\n",
    "target = np.zeros((2, 7))\n",
    "target[0, 1] = 1\n",
    "target[1, 4] = 1\n",
    "assert np.all(sample == target)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aace4a2",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Now we can apply the transformation to our train and test data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b956add5",
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = lb.transform(y_train)\n",
    "y_test = lb.transform(y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "040abd24",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-success\">\n",
    "\n",
    "## Checkpoint 1\n",
    "\n",
    "We have completed the data wrangling phase of this exercise\n",
    "- Reshaping the data to fit our model\n",
    "- Balancing classes\n",
    "- Applying one hot encoding\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14cc0cb8",
   "metadata": {
    "cell_marker": "r\"\"\"",
    "id": "2fhlPKpe3x1j"
   },
   "source": [
    "## The linear classifier\n",
    "The linear classifier produces class scores that are a linear function of the pixel values. Mathematically, this can be written as $\\vec{y} = W \\vec{x}$, where $\\vec{y}$ is the vector of class scores, $W$ is a matrix of weights and $\\vec{x}$ is the image vector. The shape of the weights matrix is determined by the number of classes and the length of the image vector. In this case $W$ is 7 by 4356. Our learning task is to find a set of weights that maximize our performance on our classification task. We will solve this task by doing the following steps\n",
    "- Randomly initializing a set of weights\n",
    "- Defining a loss function that measures our performance on the classification task\n",
    "- Use stochastic gradient descent to find \"optimal\" weights"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b3185c",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "W0LQ99KR3x1j"
   },
   "source": [
    "### Create the matrix of weights\n",
    "Properly initializing weights is essential for getting deep learning methods to work correctly. We are going to start with weights initizalized with zeros, but will invesigate other methods later in the exercise.\n",
    "\n",
    "Lets create the linear classifier using object oriented programming, which will help with organization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72303dc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "class LinearClassifier:\n",
    "    def __init__(self, image_size=image_width * image_width, n_classes=7):\n",
    "        self.image_size = image_size\n",
    "        self.n_classes = n_classes\n",
    "\n",
    "        # Initialize weights\n",
    "        self._initialize_weights()\n",
    "\n",
    "    def _initialize_weights(self):\n",
    "        self.W = np.zeros((self.n_classes, self.image_size))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21dfda05",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "7yeA5geI3x1k"
   },
   "source": [
    "### Apply the softmax transform to complete the model outputs\n",
    "Our `LinearClassifier` class needs a method to perform predictions - which in our case is performing matrix multiplication and then applying the softmax transform.  The softmax functions transforms a vector of arbitrary real numbers and turns it into probabilities that we can use for our final prediction.\n",
    "\n",
    "![](https://miro.medium.com/max/1400/1*gctBX5YHUUpBEK3MWD6r3Q.png)\n",
    "\n",
    "(*Image by [Thomas Kurbiel](https://towardsdatascience.com/derivative-of-the-softmax-function-and-the-categorical-cross-entropy-loss-ffceefc081d1)*)\n",
    "\n",
    "There is an excellent [derivation of the softmax function](https://towardsdatascience.com/derivative-of-the-softmax-function-and-the-categorical-cross-entropy-loss-ffceefc081d1) available on Towards Data Science if you are interested in the details of the math."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33e7f236",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def predict(self, X, epsilon=1e-5):\n",
    "    # Matrix multiplication of X by the weights\n",
    "    y = np.matmul(X, self.W.T)\n",
    "\n",
    "    # Apply softmax - epsilon added for numerical stability\n",
    "    y = np.exp(y) / np.sum(np.exp(y) + epsilon, axis=-1, keepdims=True)\n",
    "    return y\n",
    "\n",
    "\n",
    "# Assign methods to class\n",
    "setattr(LinearClassifier, \"predict\", predict)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8037cd7",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Before we train the model, let's take a brief moment to check what the untrained model predictions look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3408015a",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "4tgJymU33x1l",
    "outputId": "736e325f-56eb-4bf7-9caf-bc6591c2a448"
   },
   "outputs": [],
   "source": [
    "lc = LinearClassifier()\n",
    "\n",
    "fig, axes = plt.subplots(2, 4, figsize=(20, 10))\n",
    "for i, j in enumerate(np.random.randint(X_test.shape[0], size=(8,))):\n",
    "    # Get an example image\n",
    "    X_sample = X_test[j, ...]\n",
    "\n",
    "    # Reshape flattened vector to image\n",
    "    X_reshape = np.reshape(X_sample, (66, 66))\n",
    "\n",
    "    # Predict the label\n",
    "    y_pred = lc.predict(X_sample)\n",
    "\n",
    "    # Display results\n",
    "    axes.flatten()[i].imshow(X_reshape, cmap=\"gray\")\n",
    "    axes.flatten()[i].set_title(\n",
    "        \"Label \" + str(np.argmax(y_test[j])) + \", Prediction \" + str(np.argmax(y_pred))\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b338c9e",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "1iKRfu3J3x1l"
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 2.1\n",
    "What do you notice about the initial results of the model?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3261f90c",
   "metadata": {
    "cell_marker": "r\"\"\"",
    "id": "jv4Rc_xS3x1l"
   },
   "source": [
    "## Stochastic gradient descent\n",
    "To train this model, we will use stochastic gradient descent. In its simplest version, this algorithm consists of the following steps:\n",
    "- Select several images from the training dataset at random\n",
    "- Compute the gradient of the loss function with respect to the weights, given the selected images\n",
    "- Update the weights using the update rule $W_{ij} \\rightarrow W_{ij} - lr\\frac{\\partial loss}{\\partial W_{ij}}$\n",
    "\n",
    "Recall that the origin of this update rule is from multivariable calculus - the gradient tells us the direction in which the loss function increases the most. So to minimize the loss function we move in the opposite direction of the gradient.\n",
    "\n",
    "Also recall from the course notes that for this problem we can compute the gradient analytically. The gradient is given by\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{\\partial loss}{\\partial W_{ij}} = \\left(p_i - 1(i \\mbox{ is correct}) \\right)x_j,\n",
    "\\end{equation}\n",
    "where $1$ is an indicator function that is 1 if the statement inside the parentheses is true and 0 if it is false.\n",
    "\n",
    "A complete derivation of $\\frac{\\partial loss}{\\partial W_{ij}}$ is included in the Towards Data Science [article](https://towardsdatascience.com/derivative-of-the-softmax-function-and-the-categorical-cross-entropy-loss-ffceefc081d1) recommended above if you are interested in the details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8502b29",
   "metadata": {
    "id": "XXHGyfz33x1l"
   },
   "outputs": [],
   "source": [
    "def grad(self, X, y_true, y_pred):\n",
    "    # Compute the gradients for each class and save in list\n",
    "    gradients = []\n",
    "    for i in range(self.n_classes):\n",
    "        # Calculate the difference between the class probability and true score\n",
    "        difference = y_pred[..., i] - y_true[..., i]\n",
    "        difference = np.expand_dims(difference, axis=-1)\n",
    "        grad = difference * X\n",
    "        gradients.append(grad)\n",
    "\n",
    "    gradient = np.stack(gradients, axis=1)\n",
    "\n",
    "    return gradient\n",
    "\n",
    "\n",
    "def loss(self, X, y_true, y_pred):\n",
    "    loss = np.mean(-y_true * np.log(y_pred))\n",
    "\n",
    "    return loss\n",
    "\n",
    "\n",
    "def fit(self, X_train, y_train, n_epochs, batch_size=1, learning_rate=1e-5):\n",
    "    loss_list = []\n",
    "\n",
    "    rng = np.random.default_rng()\n",
    "\n",
    "    # Iterate over epochs\n",
    "    for epoch in range(n_epochs):\n",
    "        print(f\"Epoch {epoch}\")\n",
    "        n_batches = int(np.floor(X_train.shape[0] / batch_size))\n",
    "\n",
    "        # Generate random index\n",
    "        index = np.arange(X_train.shape[0])\n",
    "        np.random.shuffle(index)\n",
    "        X_shfl = X_train[index]\n",
    "        y_shfl = y_train[index]\n",
    "\n",
    "        # Iterate over batches\n",
    "        for batch in tqdm.trange(n_batches):\n",
    "            beg = batch * batch_size\n",
    "            end = (\n",
    "                (batch + 1) * batch_size\n",
    "                if (batch + 1) * batch_size < X_train.shape[0]\n",
    "                else -1\n",
    "            )\n",
    "            X_batch = X_shfl[beg:end]\n",
    "            y_batch = y_shfl[beg:end]\n",
    "\n",
    "            # Skip empty batch if it shows up at the end of the epoch\n",
    "            if X_batch.shape[0] == 0:\n",
    "                continue\n",
    "\n",
    "            # Predict\n",
    "            y_pred = self.predict(X_batch)\n",
    "\n",
    "            # Compute the loss\n",
    "            loss = self.loss(X_batch, y_batch, y_pred)\n",
    "            loss_list.append(loss)\n",
    "\n",
    "            # Compute the gradient\n",
    "            gradient = self.grad(X_batch, y_batch, y_pred)\n",
    "\n",
    "            # Compute the mean gradient over all the example images\n",
    "            gradient = np.mean(gradient, axis=0, keepdims=False)\n",
    "\n",
    "            # Update the weights\n",
    "            self.W -= learning_rate * gradient\n",
    "\n",
    "            if np.count_nonzero(np.isnan(self.W)) != 0:\n",
    "                print(epoch, batch)\n",
    "                break\n",
    "                \n",
    "        print('Final loss', loss)\n",
    "\n",
    "    return loss_list\n",
    "\n",
    "\n",
    "# Assign methods to class\n",
    "setattr(LinearClassifier, \"grad\", grad)\n",
    "setattr(LinearClassifier, \"loss\", loss)\n",
    "setattr(LinearClassifier, \"fit\", fit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d61bcc65",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "We're ready to train our model!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "231c175c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "lc = LinearClassifier()\n",
    "loss_log = lc.fit(X_train, y_train, n_epochs=10, batch_size=16)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c64fb0da",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "Let's plot the loss curve to see how the model trained."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "22e800a7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def smooth(scalars, weight):\n",
    "    \"\"\"Compute the exponential moving average to smooth data\n",
    "\n",
    "    Credit: https://stackoverflow.com/questions/42281844/what-is-the-mathematics-behind-the-smoothing-parameter-in-tensorboards-scalar/49357445#49357445\n",
    "\n",
    "    \"\"\"\n",
    "    last = scalars[0]  # First value in the plot (first timestep)\n",
    "    smoothed = list()\n",
    "    for point in scalars:\n",
    "        smoothed_val = last * weight + (1 - weight) * point  # Calculate smoothed value\n",
    "        smoothed.append(smoothed_val)  # Save it\n",
    "        last = smoothed_val  # Anchor the last smoothed value\n",
    "\n",
    "    return smoothed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58472143",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(smooth(loss_log, 0.9))\n",
    "ax.set_ylabel(\"loss\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98186b13",
   "metadata": {
    "cell_marker": "\"\"\"",
    "id": "Hh-GWNb-3x1m"
   },
   "source": [
    "## Evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "904e65d8",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "T_wsv4zW3x1m",
    "outputId": "1d5928f9-1acb-43ec-8c58-ad756e538171"
   },
   "outputs": [],
   "source": [
    "# Visualize some predictions\n",
    "fig, axes = plt.subplots(2, 4, figsize=(20, 10))\n",
    "for i, j in enumerate(np.random.randint(X_test.shape[0], size=(8,))):\n",
    "    # Get an example image\n",
    "    X_sample = X_test[j]\n",
    "\n",
    "    # Reshape flattened vector to image\n",
    "    X_reshape = np.reshape(X_sample, (66, 66))\n",
    "\n",
    "    # Predict the label\n",
    "    y_pred = lc.predict(X_sample)\n",
    "\n",
    "    # Display results\n",
    "    axes.flatten()[i].imshow(X_reshape, cmap=\"gray\")\n",
    "    axes.flatten()[i].set_title(\n",
    "        \"Label \" + str(np.argmax(y_test[j])) + \", Prediction \" + str(np.argmax(y_pred))\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a243b701",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "In addition to inspecting the results of individual predictions, we can also look at summary statistics that capture model performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a10dfde8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def benchmark_performance(y_true, y_pred):\n",
    "    \"\"\"Calculates recall, precision, f1 and a confusion matrix for sample predictions\n",
    "\n",
    "    Args:\n",
    "        y_true (list): List of integers of true class values\n",
    "        y_pred (list): List of integers of predicted class value\n",
    "\n",
    "    Returns:\n",
    "        dict: Dictionary with keys `recall`, `precision`, `f1`, and `cm`\n",
    "\n",
    "    \"\"\"\n",
    "    _round = lambda x: round(x, 3)\n",
    "\n",
    "    metrics = {\n",
    "        \"recall\": _round(sklearn.metrics.recall_score(y_true, y_pred, average=\"macro\")),\n",
    "        \"precision\": _round(\n",
    "            sklearn.metrics.precision_score(y_true, y_pred, average=\"macro\")\n",
    "        ),\n",
    "        \"f1\": _round(sklearn.metrics.f1_score(y_true, y_pred, average=\"macro\")),\n",
    "        \"cm\": sklearn.metrics.confusion_matrix(y_true, y_pred, normalize=None),\n",
    "        \"cm_norm\": sklearn.metrics.confusion_matrix(y_true, y_pred, normalize=\"true\"),\n",
    "    }\n",
    "\n",
    "    return metrics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c01c42b3",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 2.2\n",
    "\n",
    "For each of the 4 metrics above, describe in your own words what this metric tells you about model performance.\n",
    "\n",
    "- Recall\n",
    "- Precision\n",
    "- F1\n",
    "- Confusion Matrix\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6376bbf-67fe-40ab-b8bb-2e311dc78413",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "*Write your answers here*\n",
    "- Recall\n",
    "- Precision\n",
    "- F1 Score\n",
    "- Confusion Matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "660b8fe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_metrics(metrics, name, ax=None):\n",
    "    \"\"\"Plots a confusion matrix with summary statistics listed above the plot\n",
    "\n",
    "    The annotations on the confusion matrix are the total counts while\n",
    "    the colormap represents those counts normalized to the total true items\n",
    "    in that class.\n",
    "\n",
    "    Args:\n",
    "        metrics (dict): Dictionary output of `benchmark_performance`\n",
    "        name (str): Title for the plot\n",
    "        ax (optional, matplotlib subplot): Subplot axis to plot onto.\n",
    "            If not provided, a new plot is created\n",
    "        classes (optional, list): A list of the classes to label the X and y\n",
    "            axes. Defaults to [0, 1] for a two class problem.\n",
    "    \"\"\"\n",
    "    if ax is None:\n",
    "        fig, ax = plt.subplots(figsize=(5, 5))\n",
    "    cb = ax.imshow(metrics[\"cm_norm\"], cmap=\"Greens\", vmin=0, vmax=1)\n",
    "\n",
    "    classes = np.arange(metrics[\"cm\"].shape[0])\n",
    "    ax.set_xticks(range(len(classes)), class_lut)\n",
    "    ax.set_yticks(range(len(classes)), class_lut)\n",
    "    ax.set_xlabel(\"Predicted Label\")\n",
    "    ax.set_ylabel(\"True Label\")\n",
    "\n",
    "    for i in range(len(classes)):\n",
    "        for j in range(len(classes)):\n",
    "            color = \"green\" if metrics[\"cm_norm\"][i, j] < 0.5 else \"white\"\n",
    "            ax.annotate(\n",
    "                \"{}\".format(metrics[\"cm\"][i, j]),\n",
    "                (j, i),\n",
    "                color=color,\n",
    "                va=\"center\",\n",
    "                ha=\"center\",\n",
    "            )\n",
    "\n",
    "    _ = plt.colorbar(cb, ax=ax)\n",
    "    _ = ax.set_title(\n",
    "        \"{}\\n\"\n",
    "        \"Recall: {}\\n\"\n",
    "        \"Precision: {}\\n\"\n",
    "        \"F1 Score: {}\\n\"\n",
    "        \"\".format(name, metrics[\"recall\"], metrics[\"precision\"], metrics[\"f1\"])\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46ad5e16-2406-4472-a586-29fed8e3b3e3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "def summarize_performance(model, X_train, y_train, X_test, y_test, title=''):\n",
    "    \"\"\"Quick function to generate predictions on the train and test splits, \n",
    "    benchmark and plot the results\n",
    "    \"\"\"\n",
    "    \n",
    "    # Generate predictions and metrics for training data\n",
    "    y_pred = model.predict(X_train)\n",
    "    # Convert from one hot encoding to original class labels\n",
    "    y_pred = np.argmax(y_pred, axis=-1)\n",
    "    y_true = np.argmax(y_train, axis=-1)\n",
    "    train_metrics = benchmark_performance(y_true, y_pred)\n",
    "\n",
    "    # Generate predictions and metrics for test data\n",
    "    y_pred = model.predict(X_test)\n",
    "    # Convert from one hot encoding to original class labels\n",
    "    y_pred = np.argmax(y_pred, axis=-1)\n",
    "    y_true = np.argmax(y_test, axis=-1)\n",
    "    test_metrics = benchmark_performance(y_true, y_pred)\n",
    "\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(10, 4))\n",
    "\n",
    "    plot_metrics(train_metrics, f\"{title} Training\", ax[0])\n",
    "    plot_metrics(test_metrics, f\"{title} Testing\", ax[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "59816c1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "summarize_performance(lc, X_train, y_train, X_test, y_test, \"Linear Classifier\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d07077dc",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 2.3\n",
    "\n",
    "What do you notice about the results after training the model?\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b131f659",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "## Checkpoint 2\n",
    "\n",
    "We have written a simple linear classifier, trained it and considered a few ways to evaluate model performance. Next we'll look at some other machine learning methods for classification.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b398cfc8",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Part B: Random Forest Classifier\n",
    "\n",
    "Decisions trees are a useful tool for generating interpretable classification results. As shown in the example below, trees are constructed such that the data is split at each node according to a feature in the data. At the bottom of the tree, the leafs should correspond to a single class such that we can predict the class of the data depending on which leaf it is associated with.\n",
    "\n",
    "![](https://miro.medium.com/v2/resize:fit:1262/format:webp/1*LMoJmXCsQlciGTEyoSN39g.jpeg)\n",
    "\n",
    "(*Image by [Tony Yiu](https://towardsdatascience.com/understanding-random-forest-58381e0602d2)*)\n",
    "\n",
    "A random forest classifer is much like what it sounds. It takes predictions from many different decision trees and assigns the class with the most votes. Ultimately this ensemble method of considering many different trees performs better than any one decision tree.\n",
    "\n",
    "![](https://miro.medium.com/v2/resize:fit:1052/format:webp/1*VHDtVaDPNepRglIAv72BFg.jpeg)\n",
    "\n",
    "(*Image by [Tony Yiu](https://towardsdatascience.com/understanding-random-forest-58381e0602d2)*)\n",
    "\n",
    "In this exercise, we will use `scikit-learn's` implementation of `RandomForestClassifier` ([docs](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b7b9090",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "rfc = sklearn.ensemble.RandomForestClassifier()\n",
    "rfc.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b6acfc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "summarize_performance(rfc, X_train, y_train, X_test, y_test, \"Random Forest Classifier\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57f8f16d",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "## Parameter Optimization\n",
    "\n",
    "Our initial random forest classifier was trained using the default parameters provided by `sklearn`, but these often won't be the right values for our problem. In many situations, we may have some idea what a reasonable parameter value might be, but most of the time we will need to perform a grid search to select the optimal value. `sklearn` provides a class `RandomizedSearchCV` ([docs](https://scikit-learn.org/stable/modules/generated/sklearn.model_selection.RandomizedSearchCV.html)) to perform a random search over a provided grid of parameters, which we can use to improve the performance of the random forest classifier."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "deb23125",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 3.1\n",
    "\n",
    "For this task, your job is to select a set of parameters to use for optimization. You should start by looking at the [docs](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestClassifier.html) for `RandomForestClassifier` to see what parameters are available. For each parameter that you choose to optimize, think about what range of values you should try.\n",
    "\n",
    "For example:\n",
    "```\n",
    "parameter_grid = {\n",
    "    'criterion': ['gini', 'entropy', 'log_loss'],\n",
    "    'max_depth': [1, 2, 3]\n",
    "}\n",
    "```\n",
    "\n",
    "Fill out the dictionary below with your parameter selection.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "713c5a6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "# Here's a few possible parameter types to consider, but there are more possibilities\n",
    "\n",
    "parameter_grid = {\n",
    "    \"n_estimators\": ...,\n",
    "    \"criterion\": ...,\n",
    "    \"max_depth\": ...,\n",
    "    \"min_samples_leaf\": ...\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "330bf7d6",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "If we had more time, we could use `RandomizedSearchCV` to look for the best set of parameters:\n",
    "\n",
    "```python\n",
    "rf_random = sklearn.model_selection.RandomizedSearchCV(\n",
    "    estimator=sklearn.ensemble.RandomForestClassifier(),\n",
    "    param_distributions=parameter_grid,\n",
    "    n_iter=5,\n",
    "    verbose=3)\n",
    "rf_random.fit(X_train, y_train)\n",
    "print(rf_random.best_params_)\n",
    "```\n",
    "\n",
    "However to save on time we are going to distribute the work among the group.\n",
    "\n",
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "## Checkpoint 3\n",
    "\n",
    "We'll come back together as a group to discuss the parameter space that we want to explore and to sign up for parameter configurations.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e042b4f0",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 4.1\n",
    "\n",
    "Configure the `RandomForestClassifier` below with a set of parameters according to the parameter space we discussed during the previous checkpoint. Then train and evaluate your classifier!\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2cb39563",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "random_rfc = sklearn.ensemble.RandomForestClassifier(\n",
    "    # Add your parameter configuration here\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8add9e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%time\n",
    "random_rfc.fit(X_train, y_train)\n",
    "\n",
    "summarize_performance(rfc, X_train, y_train, X_test, y_test, \"Random Forest Classifier\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e37f97c",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "## Checkpoint 4\n",
    "\n",
    "Report the performance of your model on the test split in the group spreadsheet (check the course chat for a link). We'll discuss what seems to be working best as a group.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "56825297",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "# Part C: Feature Engineering\n",
    "\n",
    "Classical machine learning methods often turn to manual feature engineering to extract elements of the data that the model will use for prediction. So far we have relied on the raw data alone, but in some cases well designed features can produce a better model than raw data alone."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "667bce02",
   "metadata": {
    "cell_marker": "r\"\"\""
   },
   "source": [
    "## Introduction to image filters\n",
    "\n",
    "Image filters operate by taking a small kernel (or matrix) and applying it to each pixel in the image to compute a new value. For example, this 3x3 kernel will sharpen an image\n",
    "\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "    0 & -1 & 0 \\\\\n",
    "    -1 & 5 & -1 \\\\\n",
    "    0 & -1 & 0 \\\\\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "Filters can produce a variety of effects on images depending on how the kernel is configured. This can range from blurring an image to extracting edges. Filtered images can contain data that is more informative to the model when distinguishing between classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f9553ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "im = X_data[np.random.randint(X_data.shape[0]), ..., 0]\n",
    "fig, ax = plt.subplots(1, 3, figsize=(10, 4))\n",
    "\n",
    "ax[0].imshow(im, cmap=\"Greys_r\")\n",
    "ax[0].set_title(\"Original\")\n",
    "\n",
    "ax[1].imshow(skimage.filters.gaussian(im), cmap=\"Greys_r\")\n",
    "ax[1].set_title(\"Gaussian\")\n",
    "\n",
    "ax[2].imshow(skimage.filters.laplace(im), cmap=\"Greys_r\")\n",
    "ax[2].set_title(\"Laplace\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30c3d513",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "A variety of filters are made available through the `skimage.filters` [module](https://scikit-image.org/docs/stable/api/skimage.filters.html). In this part of the exercise, we are going to explore how filters can be applied to images in order to extract features for model prediction. While we are going to work with the filters that are easily available through skimage, there are many other transformations that can be applied to images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c83da262",
   "metadata": {},
   "outputs": [],
   "source": [
    "from inspect import getmembers, isfunction\n",
    "\n",
    "# List all functions in the skimage.filters module to get a list of available filters\n",
    "[m[0] for m in getmembers(skimage.filters, isfunction)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a464e38f",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 5.1\n",
    "\n",
    "Test a variety of the available filters from the `skimage` module. Whenever you are making a modification to an image, you should check the results to make sure that errors are not introduced while generating the transformation. The easiest way to plot an image is just to run `plt.imshow(image)`.\n",
    "\n",
    "Ultimately we are going to use model performance to select the best features for our classification task, but you should be familiar with the output of any filters that you are using. The goal of this next section is to identify a set of candidate filters from which one will be chosen that you think will lead to better classification results on the two classes we are trying to distinguish.\n",
    "\n",
    "Keep the following things you may want to keep in mind as you approach this problem\n",
    "- Look at several randomly selected images from the different classes when you are testing a filter\n",
    "- Explore the effect of parameters available for each filter\n",
    "\n",
    "**Challenge**: While this task could be approached by testing filters one at a time, consider writing a for loop to rapidly test filters in an automated fashion.\n",
    "\n",
    "\n",
    "**Tip**: Within a jupyter notebook, you can run `function?` to look at the documentation for that function. Check out the example below.\n",
    "\n",
    "</div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "20d0fd27",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "skimage.filters.gaussian?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "043fee3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "##########################\n",
    "######## To Do ###########\n",
    "##########################\n",
    "\n",
    "# Put your code for testing filters here\n",
    "\n",
    "# Hint\n",
    "# Pick a test image\n",
    "im = ...\n",
    "\n",
    "# Apply a filter \n",
    "filtered = ...\n",
    "\n",
    "# Plot the before and after\n",
    "fig, ax = plt.subplots(1, 2)\n",
    "ax[0].imshow(im, cmap='Greys_r')\n",
    "ax[1].imshow(filtered, cmap='Greys_r')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397f5666-3b44-4c3d-a31e-d0ee0ac022d7",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "## Checkpoint 6 \n",
    "\n",
    "When you have selected a few filters that you think are interesting please post a screenshot to the group chat. \n",
    "\n",
    "We'll take a look at the results as a group and then move on to training a model on filtered data.\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b6dab14",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "\n",
    "#### Task 6.1\n",
    "\n",
    "Your challenge now is to combine the features that you have explored with one of the two models we've trained so far. It's up to you to select a set of features and model design that you think will produce the best results. We'll compare results as a group at the end.\n",
    "\n",
    "As you approach this task, you are encouraged to copy/paste code from earlier in the exercise and modify it as needed. Here's a rough outline of what you will need to do:\n",
    "- Build a new dataset with your selected filters. The easiest way to combine data from multiple filters is to concatenate images together along the axis corresponding to the image size. For example if our initial training data has the shape `(n_images, image_width * image_width)`, the new data will have the shape `(n_images, n_filters * image_width * image_width)`. Here's a code snippet to get you started:\n",
    "```python\n",
    "# Select single image (10) and channel (0)\n",
    "im = X_data[10, ..., 0]\n",
    "gaus = skimage.filters.gaussian(im)\n",
    "lapl = skimage.filters.laplace(im)\n",
    "\n",
    "# Flatten each image\n",
    "im = np.reshape(im, (image_width * image_width))\n",
    "gaus = np.reshape(gaus, (image_width * image_width))\n",
    "lapl = np.reshape(lapl, (image_width * image_width))\n",
    "\n",
    "# Concatenate together\n",
    "ims = np.concatenate([im, gaus, lapl], axis=-1)\n",
    "```\n",
    "\n",
    "- Finish assembling your dataset by splitting into train/test split (`sklearn.model_selection.train_test_split`), balancing classes (`imblearn.over_sampling.RandomOverSampler`) and one hot encoding (`sklearn.label_processing.LabelBinarizer`)\n",
    "\n",
    "- Configure your model. If you choose the `LinearClassifier` make sure you update the `image_size` parameter to match the dimensions of your new dataset.\n",
    "\n",
    "- Train and evaluate the results\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3dac5844",
   "metadata": {
    "cell_marker": "\"\"\""
   },
   "source": [
    "<div class=\"alert alert-block alert-success\">\n",
    "\n",
    "## Checkpoint 7\n",
    "\n",
    "Share the results of your model on the spreadsheet (see the course chat for the link) and we'll compare results as a group.\n",
    "</div>"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_markers": "\"\"\"",
   "formats": "ipynb,py:percent"
  },
  "kernelspec": {
   "display_name": "Python [conda env:01_intro_ml]",
   "language": "python",
   "name": "conda-env-01_intro_ml-py"
  },
  "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}