notebooks/ML4H_ML_intro.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction to ML using ML4H"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prerequisites\n",
    "- Basic comfort with python, some linear algebra, some data science\n",
    "- Follow the instructions in the main [README](https://github.com/broadinstitute/ml4h) for installing ML4H\n",
    "- Now we are ready to teach the machines!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "import random\n",
    "from typing import List, Dict, Callable\n",
    "from collections import defaultdict, Counter\n",
    "\n",
    "import csv\n",
    "import gzip\n",
    "import h5py\n",
    "import shutil\n",
    "import zipfile\n",
    "import pydicom\n",
    "import numpy as np\n",
    "\n",
    "from keras import metrics\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Dense, Conv2D, Flatten\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import gridspec"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Linear Regression\n",
    "We explore machine learning on Bio medical data using Cloud computing, Python, Tensorflow, and the ml4h codebase.\n",
    "\n",
    "We will start with linear regression.  Our model is a vector, one weight for each input feature, and a single bias weight.\n",
    "\n",
    "\\begin{equation}\n",
    "y = xw + b\n",
    "\\end{equation}\n",
    "\n",
    "For notational convenience absorb the bias term into the weight vector by adding a 1 to the input data matrix $X$\n",
    "\n",
    "\\begin{equation}\n",
    "y = [\\textbf{1}, X][b, \\textbf{w}]^T\n",
    "\\end{equation}\n",
    "\n",
    "#### The Dense Layer is Matrix (Tensor) Multiplication\n",
    "![Matrix Multiplication](https://www.mathwarehouse.com/algebra/matrix/images/matrix-multiplication/how-to-multiply-2-matrices-demo.gif)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def linear_regression():\n",
    "    samples = 40\n",
    "    real_weight = 2.0\n",
    "    real_bias = 0.5\n",
    "    x = np.linspace(-1, 1, samples)\n",
    "    y = real_weight*x + real_bias + (np.random.randn(*x.shape) * 0.1)\n",
    "\n",
    "    linear_model = Sequential()\n",
    "    linear_model.add(Dense(1, input_dim=1))\n",
    "    linear_model.compile(loss='mean_squared_error', optimizer='sgd')\n",
    "    linear_model.summary()\n",
    "    linear_model.fit(x, y, batch_size=1, epochs=6)\n",
    "\n",
    "    learned_slope = linear_model.get_weights()[0][0][0]\n",
    "    learned_bias = linear_model.get_weights()[1][0]\n",
    "    print('Learned slope:',  learned_slope, 'real slope:', real_weight, 'learned bias:', learned_bias, 'real bias:', real_bias)\n",
    "\n",
    "    plt.plot(x, y)\n",
    "    plt.plot([-1,1], [-learned_slope+learned_bias, learned_slope+learned_bias], 'r')\n",
    "    plt.show()\n",
    "    print('Linear Regression complete!')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "linear_regression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now Logistic Regression:\n",
    "We take the real-valued predictions from linear regression and squish them with a sigmoid.\n",
    "\n",
    "\\begin{equation}\n",
    "\\textbf{y} = \\sigma(X\\textbf{w} + b)\n",
    "\\end{equation}\n",
    "\n",
    "where \n",
    "\\begin{equation}\n",
    "\\sigma(x) = \\frac{e^x}{1+e^x} = \\frac{1}{1+e^{-x}}\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid(x):\n",
    "    a = []\n",
    "    for item in x:\n",
    "        a.append(np.exp(item)/(1+np.exp(item)))\n",
    "    return a\n",
    "\n",
    "x = np.arange(-10., 10., 0.2)\n",
    "sig = sigmoid(x)\n",
    "plt.plot(x,sig)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def logistic_regression(epochs = 600, num_labels = 10):\n",
    "    train, test, valid = load_data('mnist.pkl.gz')  \n",
    "    \n",
    "    train_y = make_one_hot(train[1], num_labels)\n",
    "    valid_y = make_one_hot(valid[1], num_labels)\n",
    "    test_y = make_one_hot(test[1], num_labels)\n",
    "\n",
    "    logistic_model = Sequential()\n",
    "    logistic_model.add(Dense(num_labels, activation='softmax', input_dim=784, name='mnist_templates'))\n",
    "    logistic_model.compile(loss='categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])\n",
    "    logistic_model.summary()\n",
    "    \n",
    "    templates = logistic_model.layers[0].get_weights()[0]\n",
    "    plot_templates(templates, 0)\n",
    "    print('weights shape:', templates.shape)\n",
    "\n",
    "    for e in range(epochs):\n",
    "        trainidx = random.sample(range(0, train[0].shape[0]), 8192)\n",
    "        x_batch = train[0][trainidx,:]\n",
    "        y_batch = train_y[trainidx]\n",
    "        logistic_model.train_on_batch(x_batch, y_batch)\n",
    "        if e % 100 == 0:\n",
    "            plot_templates(logistic_model.layers[0].get_weights()[0], e)\n",
    "            print('Logistic Model test set loss and accuracy:', logistic_model.evaluate(test[0], test_y), 'at epoch', e)\n",
    "\n",
    "\n",
    "def plot_templates(templates, epoch):\n",
    "    n = 10\n",
    "    templates = templates.reshape((28,28,n))\n",
    "    plt.figure(figsize=(16, 8))\n",
    "    for i in range(n):\n",
    "        ax = plt.subplot(2, 5, i+1)\t\t\n",
    "        plt.imshow(templates[:, :, i])\n",
    "        plt.gray()\n",
    "        ax.get_xaxis().set_visible(False)\n",
    "        ax.get_yaxis().set_visible(False)\n",
    "\n",
    "    plot_name = \"./regression_example/mnist_templates_\"+str(epoch)+\".png\"\n",
    "    if not os.path.exists(os.path.dirname(plot_name)):\n",
    "        os.makedirs(os.path.dirname(plot_name))\t\t\n",
    "    plt.savefig(plot_name)\n",
    "\n",
    "\n",
    "def make_one_hot(y, num_labels):\n",
    "    ohy = np.zeros((len(y), num_labels))\n",
    "    for i in range(0, len(y)):\n",
    "        ohy[i][y[i]] = 1.0\n",
    "    return ohy\n",
    "\n",
    "\n",
    "def load_data(dataset):\n",
    "    ''' Loads the dataset\n",
    "    :param dataset: the path to the dataset (here MNIST)'''\n",
    "    data_dir, data_file = os.path.split(dataset)\n",
    "    if data_dir == \"\" and not os.path.isfile(dataset):\n",
    "        # Check if dataset is in the data directory.\n",
    "        new_path = os.path.join(\"data\", dataset)\n",
    "        if os.path.isfile(new_path) or data_file == 'mnist.pkl.gz':\n",
    "            dataset = new_path\n",
    "\n",
    "    if (not os.path.isfile(dataset)) and data_file == 'mnist.pkl.gz':\n",
    "        from urllib.request import urlretrieve\n",
    "        origin = ('http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz')\n",
    "        print('Downloading data from %s' % origin)\n",
    "        if not os.path.exists(os.path.dirname(dataset)):\n",
    "            os.makedirs(os.path.dirname(dataset))\t\n",
    "        urlretrieve(origin, dataset)\n",
    "\n",
    "    print('loading data...')\n",
    "    f = gzip.open(dataset, 'rb')\n",
    "    if sys.version_info[0] == 3:\n",
    "        u = pickle._Unpickler(f)\n",
    "        u.encoding = 'latin1'\n",
    "        train_set, valid_set, test_set = u.load()\n",
    "    else:\n",
    "        train_set, valid_set, test_set = pickle.load(f)\n",
    "    f.close()\n",
    "\n",
    "    return train_set, valid_set, test_set\n",
    "\n",
    "def plot_mnist(sides):\n",
    "    train, _, _ = load_data('mnist.pkl.gz')\n",
    "    print(train[0].shape)\n",
    "    mnist_images = train[0].reshape((-1, 28, 28, 1))\n",
    "    sides = int(np.ceil(np.sqrt(min(sides, mnist_images.shape[0]))))\n",
    "    _, axes = plt.subplots(sides, sides, figsize=(16, 16))\n",
    "    for i in range(sides*sides):\n",
    "        axes[i // sides, i % sides].imshow(mnist_images[i, ..., 0], cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Look B4 U Learn!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_mnist(25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cross Entropy Loss:\n",
    "Our favorite loss function for categorical data.\n",
    "\\begin{equation}\n",
    "L(true, model) = -\\sum_{x\\in\\mathcal{X}} true(x)\\, \\log model(x)\n",
    "\\end{equation}\n",
    "\n",
    "Binary cross entropy with $N$ data points $x$ each with a binary label: \n",
    "\\begin{equation}\n",
    "true(x) \\in \\{0, 1\\} \\\\\n",
    "L(true, model) = -\\frac{1}{N}\\sum^N_{i=1} true(x_i)\\log(model(x_i)) + (1-true(x_i))log(1-model(x_i))\n",
    "\\end{equation}\n",
    "\n",
    "This is the Kullback Leibler divergence between the true distribution and the predicted. \n",
    "This function emerges in many fields as diverse as probability, information theory, and physics.\n",
    "What is the information difference between the truth and our model?  How much data do I lose by replacing the truth with the model's predictions. What is the temperature difference between my predictions and the truth?!\n",
    "\n",
    "Categorical cross entropy with $K$ different classes or labels: \n",
    "\\begin{equation}\n",
    "true(x) \\in \\{0, 1, 2, ..., K\\} \\\\\n",
    "L(true, model) = -\\frac{1}{N}\\sum^N_{i=1}\\sum^K_{j=1} y_{ik}\\log(q_k(x_i)))\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "logistic_regression()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Deep Models: \"Hidden\" Layers and The MultiLayerPerceptron"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def multilayer_perceptron():\n",
    "    train, test, valid = load_data('mnist.pkl.gz')\n",
    "\n",
    "    num_labels = 10\n",
    "    train_y = make_one_hot(train[1], num_labels)\n",
    "    valid_y = make_one_hot(valid[1], num_labels)\n",
    "    test_y = make_one_hot(test[1], num_labels)\n",
    "\n",
    "    mlp_model = Sequential()\n",
    "    mlp_model.add(Dense(500, activation='relu', input_dim=784))\n",
    "    mlp_model.add(Dense(num_labels, activation='softmax'))\n",
    "    mlp_model.compile(loss='categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])\n",
    "    mlp_model.summary()\n",
    "    mlp_model.fit(train[0], train_y, validation_data=(valid[0],valid_y), batch_size=32, epochs=3)\n",
    "    print('Multilayer Perceptron trained. Test set loss and accuracy:', mlp_model.evaluate(test[0], test_y))\n",
    "\n",
    "multilayer_perceptron()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convolutions Flip, Slide, Multiply, Add\n",
    "Convolutions look for their kernel in a larger signal.\n",
    "\n",
    "In convolution, you always and only find what you're looking with.\n",
    "\n",
    "Convolution and cross correlation are deeply related:\n",
    "\n",
    "\\begin{equation}\n",
    "f(t) \\circledast g(t) \\triangleq\\ \\int_{-\\infty}^\\infty f(\\tau) g(t - \\tau) \\, d\\tau. = \\int_{-\\infty}^\\infty f(t-\\tau) g(\\tau)\\, d\\tau.\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "\n",
    "![title](https://upload.wikimedia.org/wikipedia/commons/2/21/Comparison_convolution_correlation.svg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def convolutional_neural_network(filters=32, kernel_size=(3,3), padding='valid', num_labels = 10):\n",
    "    train, test, valid = load_data('mnist.pkl.gz')\n",
    "\n",
    "    train_y = make_one_hot(train[1], num_labels)\n",
    "    valid_y = make_one_hot(valid[1], num_labels)\n",
    "    test_y = make_one_hot(test[1], num_labels)\n",
    "    \n",
    "    print(train[0].shape)\n",
    "    mnist_images = train[0].reshape((-1, 28, 28, 1))\n",
    "    mnist_valid = valid[0].reshape((-1, 28, 28, 1))\n",
    "    mnist_test = test[0].reshape((-1, 28, 28, 1))\n",
    "\n",
    "    cnn_model = Sequential()\n",
    "    cnn_model.add(Conv2D(input_shape=(28, 28, 1), filters=filters, kernel_size=kernel_size, padding=padding, activation='relu'))\n",
    "    cnn_model.add(Conv2D(filters=filters, kernel_size=kernel_size, padding=padding, activation='relu'))\n",
    "    cnn_model.add(Conv2D(filters=filters, kernel_size=kernel_size, padding=padding, activation='relu'))\n",
    "    cnn_model.add(Flatten())\n",
    "    cnn_model.add(Dense(16, activation='relu'))\n",
    "    cnn_model.add(Dense(num_labels, activation='softmax'))\n",
    "    cnn_model.compile(loss='categorical_crossentropy', optimizer='sgd', metrics=['accuracy'])\n",
    "    cnn_model.summary()\n",
    "    cnn_model.fit(mnist_images, train_y, validation_data=(mnist_valid, valid_y), batch_size=32, epochs=3)\n",
    "    \n",
    "    print('Convolutional Neural Network trained. Test set loss and accuracy:', cnn_model.evaluate(mnist_test, test_y))\n",
    "\n",
    "convolutional_neural_network()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Why (and When!) is Convolution Helpful?\n",
    "- Decouples input size from model size\n",
    "- Translationally Equivariant (Not Invariant), so we can find features wherever they might occur in the signal\n",
    "- Local structure is often informative\n",
    "- But not always! (eg Tabular data)"
   ]
  }
 ],
 "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}