lectures/Lecture_22.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Deep Learning: II\n",
    "\n",
    "<font color='red'>UPDATE TO FIX ISSUES FOR APPLE M1</font>\n",
    "\n",
    "*S. R. Taylor (2022)*\n",
    "\n",
    "This lecture and notebook are based on the \"NeuralNetworks\", and \"TensorFloeCodeOnly\" lectures of of G. Richards' \"Astrostatistics\" class at Drexel University (PHYS 440/540, https://github.com/gtrichards/PHYS_440_540), which in turn are based on materials from Andy Connolly, and Ivezic et al. Chapter 9, Andy Connolly's [blog](http://connolly.github.io/introAstroML/blog/regression.html), and Aurelien Geron's [book](https://github.com/ageron/handson-ml2). \n",
    "\n",
    "##### Reading:\n",
    "\n",
    "- [Textbook](http://press.princeton.edu/titles/10159.html) Chapter 9.\n",
    "- Many blogs and videos.\n",
    "- Free online book! http://neuralnetworksanddeeplearning.com/index.html\n",
    "\n",
    "---\n",
    "\n",
    "## Contents\n",
    "* [Photo-z With Deep Learning](#one)\n",
    "* [Convolutional Neural Networks (CNNs)](#two)\n",
    "* [Autoencoders](#three)\n",
    "* [Generative Adversarial Networks (GANs)](#four)\n",
    "    \n",
    "---\n",
    "\n",
    "***Exercises required for class participation are in <font color='red'>red</font>.***\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Before starting today, let's install the two big Deep Learning packages out there. [Tensorflow](https://www.tensorflow.org/) and [PyTorch](https://pytorch.org/)-- the rivalry essentially boils down to Google versus Facebook. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install --upgrade pip"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tensorflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install tensorflow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### PyTorch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip3 install torch torchvision"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Keras\n",
    "\n",
    "**[Keras](https://keras.io/)** is a deep learning API. Essentially it is Scikit-Learn for deep neural networks.\n",
    "\n",
    "Keras needs a computational backend to handle the heavy computation.  Three popular (open sources) deep learning libraries are [TensorFlow](https://www.tensorflow.org/), Microsoft Cognitive Toolkit, and [Theano](http://www.deeplearning.net/software/theano/). TensorFlow now comes bundled with a version of Keras and that's what we'll use here (actually TensorFlow 2). PyTorch is another option. If you apply for a data-science job in industry knowing one of these tools might be the most useful thing for you to have learned.\n",
    "\n",
    "In short, for neural networks:\n",
    "\n",
    "> numpy -> tensorflow\n",
    "\n",
    "> sklearn -> keras\n",
    "\n",
    "In the same way that you can build a linear regression algorithm in numpy without using sklearn, you can build a neural network algorithm (not to mention linear regression) in tensorflow without using keras.  But just as sklearn makes our life easier, so too does keras."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Photo-z with deep learning <a class=\"anchor\" id=\"one\"></a>\n",
    "\n",
    "An important problem in galaxy imaging is being able to use observed galaxy color data to estimate its redshift. This is sometimes called \"photo-z\". Now, on your humble computer, you will build a deep neural network that is trained on a sample of SDSS galaxy color data. You'll see that `pytorch` makes this quite easy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'\n",
    "\n",
    "import torch\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.utils.data as torchdata\n",
    "\n",
    "from astroML.datasets import fetch_sdss_specgals\n",
    "from astroML.utils.decorators import pickle_results\n",
    "\n",
    "from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read in and prepare the SDSS data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fetch and prepare the data\n",
    "data = fetch_sdss_specgals()\n",
    "\n",
    "# put magnitudes into array\n",
    "# normalize to zero mean and unit variance for easier training\n",
    "datanormed = np.zeros((len(data), 6), \n",
    "                      dtype=np.float32)\n",
    "for i, band in enumerate(['u', 'g', 'r', 'i', 'z']):\n",
    "    band = 'modelMag_' + band\n",
    "    datanormed[:, i] = ((data[band] - data[band].mean()) / \n",
    "                        data[band].std())\n",
    "\n",
    "# put redshifts into array\n",
    "datanormed[:, 5] = data['z']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's define a little helper class that will make our network. It is initialized with a given number of neurons in the single hidden layer. These neurons in the hidden layer connect to $5$ input variables (colors) and spits out 1 variable (redshift). You can also see that it uses ReLU activation for the neurons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define structure of neural net\n",
    "class Net(nn.Module):\n",
    "    def __init__(self, nhidden):\n",
    "        super(Net, self).__init__()\n",
    "        self.fc_h = nn.Linear(5, nhidden)\n",
    "        self.fc_o = nn.Linear(nhidden, 1)\n",
    "\n",
    "    def forward(self, x):\n",
    "        h = F.relu(self.fc_h(x))\n",
    "        z = self.fc_o(h)\n",
    "        return z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Do a train-test split on the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# split data into 9:1 train:test\n",
    "dataset = torchdata.TensorDataset(torch.tensor(datanormed[:, 0:5]),\n",
    "                                  torch.tensor(datanormed[:, 5]).view(-1, 1))\n",
    "trainnum = datanormed.shape[0] // 10 * 9\n",
    "traindata, testdata = torchdata.random_split(dataset, \n",
    "                                             [trainnum, datanormed.shape[0] - trainnum])\n",
    "traindataloader = torchdata.DataLoader(traindata, \n",
    "                                       batch_size=128, \n",
    "                                       shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define a function to train the neural network. This looks like a lot of code, but the first block is the key. The network will have $4$ neurons in the hidden layer, and will use stochastic gradient descent to train the weights and biases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "@pickle_results('NNphotoz.pkl')\n",
    "def train_NN():\n",
    "    # 4 neurons in hidden layer\n",
    "    model = Net(4)\n",
    "    # MSE loss\n",
    "    criterion = torch.nn.MSELoss(reduction='sum') \n",
    "    # stochastic gradient descent\n",
    "    optimizer = torch.optim.SGD(model.parameters(), lr=0.001) \n",
    "    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, \n",
    "                                                           verbose=True, \n",
    "                                                           patience=5, \n",
    "                                                           threshold=1e-3)\n",
    "\n",
    "    min_valid_loss = float('inf')\n",
    "    badepochs = 0\n",
    "    for t in range(100):\n",
    "        train_loss = 0\n",
    "        for i, databatch in enumerate(traindataloader, 0):\n",
    "            photometry, redshifts = databatch\n",
    "            optimizer.zero_grad()\n",
    "            z_pred = model(photometry)\n",
    "            loss = criterion(z_pred, redshifts)\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            train_loss += loss.item()\n",
    "\n",
    "        with torch.no_grad():\n",
    "            photometry = testdata[:][0]\n",
    "            redshifts = testdata[:][1]\n",
    "            z_pred = model(photometry)\n",
    "            valid_loss = criterion(z_pred, redshifts)\n",
    "            if t % 10 == 0:\n",
    "                print('Epoch %3i: train loss %0.3e validation loss %0.3e' % (t, \\\n",
    "                        train_loss / len(traindata), valid_loss / len(testdata)))\n",
    "            # stop training if validation loss has not fallen in 10 epochs\n",
    "            if valid_loss > min_valid_loss*(1-1e-3):\n",
    "                badepochs += 1\n",
    "            else:\n",
    "                min_valid_loss = valid_loss\n",
    "                badepochs = 0\n",
    "            if badepochs == 10:\n",
    "                print('Finished training')\n",
    "                break\n",
    "        scheduler.step(valid_loss)\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's train this thing. It should take about five minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = train_NN()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot the results\n",
    "with torch.no_grad():\n",
    "    photometry = testdata[:][0]\n",
    "    redshifts = testdata[:][1]\n",
    "    z_pred = model(photometry)\n",
    "    fig = plt.figure(figsize=(5, 5))\n",
    "    fig.subplots_adjust(wspace=0.25,\n",
    "                        left=0.1, right=0.95,\n",
    "                        bottom=0.15, top=0.9)\n",
    "\n",
    "    ax = plt.axes()\n",
    "    H, zs_bins, zp_bins = np.histogram2d(redshifts.numpy().flatten(), \n",
    "                                         z_pred.numpy().flatten(), 151)\n",
    "    ax.imshow(H.T, origin='lower', interpolation='nearest', aspect='auto',\n",
    "               extent=[zs_bins[0], zs_bins[-1], zp_bins[0], zp_bins[-1]],\n",
    "               cmap=plt.cm.binary)\n",
    "    ax.plot([-0.1, 0.4], [-0.1, 0.4], ':k')\n",
    "    rms = np.sqrt(np.mean((z_pred-redshifts).numpy()**2))\n",
    "    \n",
    "    ax.text(0, 0.35, 'RMS error = %0.3f' % np.sqrt(np.mean((z_pred-\n",
    "                                                            redshifts).numpy()**2)))\n",
    "    ax.set_xlim(-0.02, 0.4001)\n",
    "    ax.set_ylim(-0.02, 0.4001)\n",
    "    ax.xaxis.set_major_locator(plt.MultipleLocator(0.1))\n",
    "    ax.yaxis.set_major_locator(plt.MultipleLocator(0.1))\n",
    "    \n",
    "    ax.set_xlabel(r'$z_{\\rm true}$')\n",
    "    ax.set_ylabel(r'$z_{\\rm fit}$')\n",
    "    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not bad at all! Still, this was a relatively simple fully-connected network. Other problems will require more complexity..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Convolutional Neural Networks (CNNs) <a class=\"anchor\" id=\"two\"></a>\n",
    "\n",
    "Recent interest in neural networks surged in 2012 when a team using a deep **[convolutional neural network (CNN)](https://en.wikipedia.org/wiki/Convolutional_neural_network)** acheived record results classifying objects in the [ImageNet](http://image-net.org/) data set.\n",
    "\n",
    "The idea behind CNNs is inspired by human visual perception. ***Each neuron in your visual cortex doesn't \"see\" all of what your eye can see at once*** and some neurons are more sensitive to one pattern over another (e.g., horizontal lines vs. vertical lines).  \n",
    "\n",
    "Moreover, the simplest deeply connected neural networks would choke on \"real\" data which has far more than 28x28 pixels, and would require following tens of millions of connections. So we use a combination of \n",
    "- \"convolution\"\n",
    "- \"pooling\" \n",
    "\n",
    "to reduce the dimensionality of the data first.\n",
    "\n",
    "![Ivezic Figure 9.19](https://www.astroml.org/_images/fig_cnn_1.png)\n",
    "\n",
    "### Convolutional Layers\n",
    "\n",
    "In a convolutional layer, each neuron is not connected to each neuron in the previous layer, but only those that are within its \"field of view\" as defined by a kernel (filter).  We slide the kernel over the input layer and the value in the next layer depends only on those pixels.\n",
    "\n",
    "![](https://miro.medium.com/max/4800/1*GcI7G-JLAQiEoCON7xFbhg.gif)\n",
    "\n",
    "Here's another perspective that helps to visualize going from one layer to the next.\n",
    "\n",
    "![Convolution example](https://developer.apple.com/library/content/documentation/Performance/Conceptual/vImage/Art/kernel_convolution.jpg)\n",
    "\n",
    "Choosing a filter (kernel) with a certain pattern can help recognize certain types of features (like horizontal aor vertical lines).\n",
    "\n",
    "![vertical filter](https://miro.medium.com/max/1338/1*7IEbBib7Yc7qST5IGQoYXA.jpeg)\n",
    "\n",
    "![horizontal filter](https://miro.medium.com/max/1238/1*PSSAaH2pZbl5bK3Ef_zk4A.jpeg)\n",
    "\n",
    "The purpose of the convolutional layers is really to capture high-level features in the images, like edges. \n",
    "\n",
    "### Pooling Layers\n",
    "\n",
    "It is common to follow the convolutional layer by a so-called \"pooling layer\", essentially to reduce the amount of data that needs to be processed.  The full architechture of a CNN might look something like this:\n",
    "\n",
    "![CNN Example](https://www.researchgate.net/profile/Xian_Wei2/publication/331986652/figure/fig1/AS:740547106988032@1553571597647/The-classic-structure-of-CNN-It-consists-of-two-modules-Feature-extraction-module-and.ppm)\n",
    "\n",
    "Where the pooling layers are reducing the number of pixels by averaging, summing, taking the max, etc. \"**Max Pool**\" is the most popular, and literally corresponds to taking the maximum of the convolution.\n",
    "\n",
    "![Pooling example](https://miro.medium.com/max/1000/1*ydNsGDxMldAiq7b96GDQwg.jpeg)\n",
    "\n",
    "or another way to see this is...\n",
    "\n",
    "![](https://miro.medium.com/max/792/1*uoWYsCV5vBU8SHFPAPao-w.gif)\n",
    "\n",
    "When we are done with convolution and pooling, the output gets fed into a regular, fully connected neural network that outputs the predictions.\n",
    "\n",
    "- This is clearly much more sophisticated than our basic fully-connected multi-layer perceptron. \"Deep\" networks consist of tens of layers with thousands of neurons. **These large networks have become usable thanks to two breakthroughs: the use of sparse layers, and the power of graphics processing units (GPUs).**\n",
    "\n",
    "\n",
    "- Sparse layers or convolutional layers in a deep network contain a large number of hidden nodes but very few synapses. The sparseness arises from the relatively small size of a typical convolution kernel (15x15 is a large kernel), so a hidden node representing one output of the convolution is connected to only a few input nodes. Compare this the our previous perceptron, in which every hidden node was connected to every input node.\n",
    "\n",
    "\n",
    "- Even though the total number of connections is greatly reduced in the sparse layers, the total number of nodes and connections in a modern deep network is still enormous. Luckily, training these networks turns out to be a great task for GPU acceleration.\n",
    "\n",
    "\n",
    "For further study, there are lots of resources for CNNs online.  For example, see\n",
    "https://medium.com/analytics-vidhya/convolutional-neural-networks-cnn-explained-step-by-step-69137a54e5e7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following cells are from Geron, Chapter 10, see \n",
    "https://github.com/ageron/handson-ml2/blob/master/10_neural_nets_with_keras.ipynb\n",
    "\n",
    "We'll start by introducing the **Fashion MNIST data set**. The next cells load the data, define test, trainging, and validation sets; normalize the data; display an example image; list the possible target values ($y$), and show a 4x10 grid of examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fashion_mnist = keras.datasets.fashion_mnist\n",
    "(X_train_full, y_train_full), (X_test, y_test) = fashion_mnist.load_data()\n",
    "\n",
    "print(X_train_full.shape)\n",
    "print(X_test.shape)\n",
    "\n",
    "X_valid, X_train = X_train_full[:5000] / 255., X_train_full[5000:] / 255.\n",
    "y_valid, y_train = y_train_full[:5000], y_train_full[5000:]\n",
    "X_test = X_test / 255.\n",
    "\n",
    "print(X_train.shape)\n",
    "print(X_valid.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This data set is a huge collection of images of items of clothing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's make a list of class names that we can refer to.\n",
    "class_names = [\"T-shirt/top\", \"Trouser\", \"Pullover\", \"Dress\", \"Coat\",\n",
    "               \"Sandal\", \"Shirt\", \"Sneaker\", \"Bag\", \"Ankle boot\"]\n",
    "\n",
    "# Plot 4x10 array of images from the Fashion MNIST database\n",
    "n_rows = 4\n",
    "n_cols = 10\n",
    "plt.figure(figsize=(n_cols * 1.2, n_rows * 1.2))\n",
    "for row in range(n_rows):\n",
    "    for col in range(n_cols):\n",
    "        index = n_cols * row + col\n",
    "        plt.subplot(n_rows, n_cols, index + 1)\n",
    "        plt.imshow(X_train[index], cmap=\"binary\", interpolation=\"nearest\", origin='upper')\n",
    "        plt.axis('off')\n",
    "        plt.title(class_names[y_train[index]], fontsize=12)\n",
    "plt.subplots_adjust(wspace=0.2, hspace=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll use `tensorflow` on this data set, through it's handy API `keras`. Since these are images and we want to preserve contextual information like edges, we'll build a CNN."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow import keras\n",
    "\n",
    "keras.backend.clear_session() # Make sure that we are starting a new model and not adding to an earlier one\n",
    "np.random.seed(42) # Set the numpy and tensorflow random seeds so that we all get the same answer\n",
    "tf.random.set_seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Normalize and prepare the data to pass into a network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(X_train_full, y_train_full), (X_test, y_test) = keras.datasets.fashion_mnist.load_data()\n",
    "X_train, X_valid = X_train_full[:-5000], X_train_full[-5000:]\n",
    "y_train, y_valid = y_train_full[:-5000], y_train_full[-5000:]\n",
    "\n",
    "X_mean = X_train.mean(axis=0, keepdims=True)\n",
    "X_std = X_train.std(axis=0, keepdims=True) + 1e-7\n",
    "X_train = (X_train - X_mean) / X_std\n",
    "X_valid = (X_valid - X_mean) / X_std\n",
    "X_test = (X_test - X_mean) / X_std\n",
    "\n",
    "\n",
    "# Need to reshape for CNN\n",
    "X_train = X_train[:, :, :, np.newaxis]\n",
    "X_valid = X_valid[:,  :, :, np.newaxis]\n",
    "X_test = X_test[:, :, :, np.newaxis]\n",
    "\n",
    "print(len(X_train))\n",
    "print(X_train.shape)\n",
    "print(X_valid.shape)\n",
    "print(X_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next cell looks like a lot, but really it is just a sequential declaration of different types of layer in the network. It is a combination of **2D Convolution layers** (feature finding), **Max Pooling layers** (data reduction), folloed by some **Dense layers** (fully-connected layers at the end of the network for classification) and some **Dropout** incorporated to strengthen the overall network and make sure all the neurons are actually pulling their weight."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = keras.models.Sequential([\n",
    "    keras.layers.InputLayer(input_shape=[28, 28, 1]),\n",
    "    keras.layers.Conv2D(filters=64, kernel_size=7, activation='relu', padding='same'),\n",
    "    keras.layers.MaxPooling2D(pool_size=2),\n",
    "    keras.layers.Conv2D(128, 3, activation='relu', padding='same'),\n",
    "    keras.layers.Conv2D(128, 3, activation='relu', padding='same'),  \n",
    "    keras.layers.MaxPooling2D(pool_size=2),\n",
    "    keras.layers.Conv2D(256, 3, activation='relu', padding='same'),\n",
    "    keras.layers.Conv2D(256, 3, activation='relu', padding='same'),\n",
    "    keras.layers.MaxPooling2D(pool_size=2),\n",
    "    keras.layers.Flatten(),\n",
    "    keras.layers.Dense(units=128, activation='relu'),\n",
    "    keras.layers.Dropout(0.5),\n",
    "    keras.layers.Dense(units=64, activation='relu'),\n",
    "    keras.layers.Dropout(0.5),\n",
    "    keras.layers.Dense(units=10, activation='softmax'),\n",
    "])   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model.summary())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**That's right...your CNN has ~1.5 million parameters that you are going to try to optimize on your machine...**\n",
    "\n",
    "We'll only go through one epoch of optimization, but hopefully that will be enough to get us above $50\\%$ accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### This may take a while! Hopefully no more than 5 minutes for 1 epoch.\n",
    "model.compile(loss=\"sparse_categorical_crossentropy\", optimizer=\"nadam\", metrics=[\"accuracy\"])\n",
    "clf = model.fit(X_train, y_train, epochs=1, validation_data=(X_valid, y_valid))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Evaluate the test set\n",
    "model.evaluate(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Predict values for first 3 test objects"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_new = X_test[:3]\n",
    "y_proba = model.predict(X_new).round(2)\n",
    "print(y_proba)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, you can see that the output is a probability that the object belongs to each class (which has to sum to 1 across all the classes).  If we just want an \"answer\", we assign it to the class with the highest probability (done here with `predict_classes`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_pred = np.argmax(model.predict(X_new), axis=-1) #New way\n",
    "print(y_pred)\n",
    "print(np.array(class_names)[y_pred])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot those and see if the predictions make sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(7.2, 2.4))\n",
    "for index, image in enumerate(X_new):\n",
    "    plt.subplot(1, 3, index + 1)\n",
    "    plt.imshow(image, cmap=\"binary\", interpolation=\"nearest\", origin='upper')\n",
    "    plt.axis('off')\n",
    "    plt.title(class_names[y_test[index]], fontsize=12)\n",
    "plt.subplots_adjust(wspace=0.2, hspace=0.5)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that after you fit your model, you can save it and reload it at some later time (which is good because some models might take hours to train!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save model\n",
    "model.save(\"my_keras_model.h5\")\n",
    "\n",
    "# Reload model\n",
    "model = keras.models.load_model(\"my_keras_model.h5\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's predict the values for 10 random objects. Display them with their actual labels first, then predict and display the predictions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx10 = np.random.choice(np.arange(len(y_test)), size=10, replace=False)\n",
    "X_new = X_test[idx10]\n",
    "y_pred = np.argmax(model.predict(X_new), axis=-1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Correct answers (y_test)\n",
    "plt.figure(figsize=(15, 5))\n",
    "for index, image in enumerate(X_new):\n",
    "    plt.subplot(1, 10, index + 1)\n",
    "    plt.imshow(image, cmap=\"binary\", interpolation=\"nearest\", origin='upper')\n",
    "    plt.axis('off')\n",
    "    plt.title(class_names[y_test[idx10[index]]], fontsize=12)\n",
    "plt.subplots_adjust(wspace=0.2, hspace=0.5)\n",
    "plt.suptitle('ACTUAL LABELS',y=0.7,fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Predicted answers (y_pred)\n",
    "plt.figure(figsize=(15, 5))\n",
    "for index, image in enumerate(X_new):\n",
    "    plt.subplot(1, 10, index + 1)\n",
    "    plt.imshow(image, cmap=\"binary\", interpolation=\"nearest\", origin='upper')\n",
    "    plt.axis('off')\n",
    "    plt.title(class_names[y_pred[index]], fontsize=12)\n",
    "plt.subplots_adjust(wspace=0.2, hspace=0.5)\n",
    "plt.suptitle('PREDICTED LABELS',y=0.7,fontsize=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not too shabby!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Autoencoders <a class=\"anchor\" id=\"three\"></a>\n",
    "\n",
    "**Autoencoders are neural networks that copy their input to their output, but after passing the data through a compression bottleneck**. For example if there are 28x28 = 784 inputs, there will also be 784 outputs, but there will be one or more (odd, but symmetric) hidden layers with fewer neuron than that.\n",
    "\n",
    "***Therefore the encoder layers learn a compressed representation of the data, with the decoder layers transforming the encodings back to the full data dimensionality.*** \n",
    "\n",
    "For example see \n",
    "\n",
    "![autoencoder structure](https://miro.medium.com/max/1400/1*44eDEuZBEsmG_TCAKRI3Kw@2x.png)\n",
    "\n",
    "from\n",
    "https://towardsdatascience.com/applied-deep-learning-part-3-autoencoders-1c083af4d798\n",
    "\n",
    "You can think of this as doing **PCA with a neural network** -- breaking our data down into the only the most important features that we actually *need* (finding the intrinsic dimensionality). In fact, if the network uses only linear (or no) activation functions and $L2$ cost function, then we have exactly PCA.  \n",
    "\n",
    "How might this be useful?  Well, for example we can use it to reconstruct MNIST digits that have had noise added to them:\n",
    "\n",
    "![autoencoder example](https://miro.medium.com/max/1400/1*SxwRp9i23OM0Up4sEze1QQ@2x.png)\n",
    "\n",
    "More exciting possibilities are: \n",
    "\n",
    "- Using autoencoders for \"**unsupervised pretraining**\". For example you have data that is only partially labeled (at least not enough to do traditional supervised classification).  We can train an autoencoder on the full data set, then used the encoder part as the base of a regular neural network that is trained on the labeled data that we do have. This is found to be a much more efficient way of initializing weights and biases than starting from random, because you are \"cheating\" a bit by already figuring out some of the data structure using the autoencoder. See Geron Figure 17.6.\n",
    "- **[Anomaly detection](https://scikit-learn.org/stable/modules/outlier_detection.html)**, e.g. https://towardsdatascience.com/anomaly-detection-with-autoencoder-b4cdce4866a6 and https://towardsdatascience.com/anomaly-detection-using-autoencoders-5b032178a1ea. See also\n",
    "https://www.pyimagesearch.com/2020/03/02/anomaly-detection-with-keras-tensorflow-and-deep-learning/ and \n",
    "https://towardsdatascience.com/autoencoder-neural-network-for-anomaly-detection-with-unlabeled-dataset-af9051a048.\n",
    "\n",
    "### Variational Autoencoders\n",
    "\n",
    "More common in astronomy are **variational autoencoders**, partly because the \"latent space\" (i.e. the compressed representation of the data) that results from a standard autoencoder doesn't necessarily map continuously to the data (e.g., if your training data don't represent the possibilities or span of the full data space). \n",
    "\n",
    "We won't go into detail, just realize that these are something that you might try if you were otherwise going to try an autoencoder to tackle your problem. **Instead of the encoder layers compressing the data down to a single point in the \"latent space\", it is mapped to a continuous distribution** (imagine if PCA gave you a probabilistic compression). The prior on this continuous distribution is a Gaussian with mean and variance. ***Therefore the intuition for variational autoencoders is that the mean vector centers the average encoding of the data point, while the variance samples around this average encoding.*** This allows the encoding and decoding to interpolate much more smoothly across the training data. \n",
    "\n",
    "![Ivezic Figure 9.21](https://www.astroml.org/_images/fig_vae_1.png)\n",
    "\n",
    "Let's build a variational autoencoder and train it on a sample of SDSS spectra. We'll switch back to `pytorch`. Grab the data, normalize and prepare for the network, then do a 3:1 train-test split."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astroML.datasets import sdss_corrected_spectra\n",
    "from astroML.plotting import setup_text_plots\n",
    "setup_text_plots(fontsize=8, usetex=False)\n",
    "\n",
    "# Fetch and prepare the data\n",
    "data = sdss_corrected_spectra.fetch_sdss_corrected_spectra()\n",
    "spectra = sdss_corrected_spectra.reconstruct_spectra(data)\n",
    "wavelengths = sdss_corrected_spectra.compute_wavelengths(data)\n",
    "\n",
    "# normalize spectra by integrated flux and subtract out mean, for easier training\n",
    "spectranorms = np.mean(spectra, axis=1)\n",
    "normedspectra = spectra / spectranorms[:, None]\n",
    "meanspectrum = np.mean(normedspectra, axis=0)\n",
    "normedspectra -= meanspectrum[None, :]\n",
    "\n",
    "# split data into 3:1 train:test\n",
    "torch.manual_seed(802)  # seed used for book figure\n",
    "dataset = torchdata.TensorDataset(torch.tensor(normedspectra))\n",
    "trainnum = normedspectra.shape[0] // 4 * 3\n",
    "traindata, testdata = torchdata.random_split(dataset, \n",
    "                                             [trainnum, normedspectra.shape[0] - trainnum])\n",
    "traindataloader = torchdata.DataLoader(traindata, batch_size=128, shuffle=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define the structure of the variational autoencoder. \n",
    "- In the top block you can see that the first layer reduces the $1000$-D input down to $250$.\n",
    "- We then pass to a $2$D mean and variance layer from which the compressed \"latent\" parameters can be probabilistically drawn. This means that we compressing down from a 1000 flux channels to a 2D compressed representation.\n",
    "- The decoding layers decompress the sampled latent parameters and scale back up to the full 1000 flux channel dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define structure of variation autoencoder\n",
    "class VAE(nn.Module):\n",
    "    def __init__(self, nhidden=250):\n",
    "        super(VAE, self).__init__()\n",
    "\n",
    "        # encoding layers\n",
    "        self.encode_fc = nn.Linear(1000, nhidden)\n",
    "        self.mu        = nn.Linear(nhidden, 2)\n",
    "        self.logvar    = nn.Linear(nhidden, 2)\n",
    "\n",
    "        # decoding layers\n",
    "        self.decode_fc = nn.Linear(2, nhidden)\n",
    "        self.output    = nn.Linear(nhidden, 1000)\n",
    "\n",
    "    def encode(self, x):\n",
    "        h = F.relu(self.encode_fc(x))\n",
    "        return self.mu(h), self.logvar(h)\n",
    "\n",
    "    def reparameterize(self, mu, logvar):\n",
    "        std = torch.exp(0.5*logvar)\n",
    "        eps = torch.randn_like(std)\n",
    "        return eps.mul(std).add_(mu)\n",
    "\n",
    "    def decode(self, z):\n",
    "        h = F.relu(self.decode_fc(z))\n",
    "        return self.output(h)\n",
    "\n",
    "    def forward(self, x):\n",
    "        mu, logvar = self.encode(x)\n",
    "        z = self.reparameterize(mu, logvar)\n",
    "        return self.decode(z), mu, logvar"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The loss is a modified version of the usual MSE, where we use the **Kullback-Leibler divergence** to asses how close the latent parameter distribution is to a Gaussian. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add KL divergence to loss function\n",
    "def VAEloss(criterion, recon_x, x, mu, logvar):\n",
    "    return criterion(recon_x, x) - 0.5 * torch.sum(1 + logvar - \n",
    "                                                   mu.pow(2) - \n",
    "                                                   logvar.exp())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another large block that is essentially just training the network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train_model():\n",
    "    model = VAE()\n",
    "    criterion = torch.nn.MSELoss(reduction='sum')\n",
    "    optimizer = torch.optim.Adam(model.parameters(), lr=1e-4)\n",
    "    scheduler = torch.optim.lr_scheduler.ReduceLROnPlateau(optimizer, \n",
    "                                                           verbose=True, \n",
    "                                                           patience=5, \n",
    "                                                           threshold=1e-3)\n",
    "\n",
    "    min_valid_loss = float('inf')\n",
    "    badepochs = 0\n",
    "    for t in range(1000):\n",
    "        train_loss = 0\n",
    "        for i, databatch in enumerate(traindataloader, 0):\n",
    "            specbatch = databatch[0]\n",
    "            optimizer.zero_grad()\n",
    "            recon, mu, logvar = model(specbatch)\n",
    "            loss = VAEloss(criterion, recon, \n",
    "                           specbatch, mu, logvar)\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "            train_loss += loss.item()\n",
    "\n",
    "        with torch.no_grad():\n",
    "            testspec = testdata[:][0]\n",
    "            recon, mu, logvar = model(testspec)\n",
    "            valid_loss = VAEloss(criterion, recon, \n",
    "                                 testspec, mu, logvar)\n",
    "            if t % 10 == 0:\n",
    "                print('Epoch %3i: train loss %6.1f validation loss %6.1f' % \\\n",
    "                        (t, train_loss / len(traindata), valid_loss / len(testdata)))\n",
    "            # stop training if validation loss has not fallen in 10 epochs\n",
    "            if valid_loss > min_valid_loss*(1-1e-3):\n",
    "                badepochs += 1\n",
    "            else:\n",
    "                min_valid_loss = valid_loss\n",
    "                badepochs = 0\n",
    "            if badepochs == 10:\n",
    "                print('Finished training')\n",
    "                break\n",
    "        scheduler.step(valid_loss)\n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's train this thing! (Takes less than a minute)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = train_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot results\n",
    "with torch.no_grad():\n",
    "    # sort latent parameters from most constrained to least constrained\n",
    "    testspec = dataset[:][0]\n",
    "    recon, mu, logvar = model(testspec)\n",
    "    zorder = np.argsort(np.mean(logvar.numpy(), axis=0))\n",
    "\n",
    "    fig = plt.figure(figsize=(10, 6))\n",
    "    fig.subplots_adjust(hspace=0, wspace=0)\n",
    "    parvalues = [-2.,0.,2.]\n",
    "    for i, z1 in enumerate(parvalues):\n",
    "        for j, z2 in enumerate(parvalues):\n",
    "            # get z1 to vary left to right, z2 bottom to top\n",
    "            ax = fig.add_subplot(3, 3, (2-j)*len(parvalues)+i+1)\n",
    "\n",
    "            z = np.zeros((1,2), dtype=np.float32)\n",
    "            z[0, zorder] = z1, z2 # set z1 is more constrained of the two latent parameters\n",
    "            spectrum = model.decode(torch.tensor(z))\n",
    "            ax.plot(wavelengths, meanspectrum+spectrum.numpy()[0,:])\n",
    "            ax.text(6750, 3, '(%i, %i)' % (z1,z2))\n",
    "\n",
    "            ax.set_xlim(3000, 8000)\n",
    "            ax.set_ylim(-1, 4)\n",
    "\n",
    "            if i == 0 and j == 1:\n",
    "                ax.set_ylabel('flux')\n",
    "            else:\n",
    "                ax.yaxis.set_major_formatter(plt.NullFormatter())\n",
    "            if j == 0 and i == 1:\n",
    "                ax.set_xlabel(r'${\\rm wavelength\\ (\\AA)}$')\n",
    "            else:\n",
    "                ax.xaxis.set_major_formatter(plt.NullFormatter())\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The spectra above are generated from our variational autoencoder. As mentioned, the encoder was limited to two components (2D compressed latent space). As neurons are progressively activated in the latent space, we generate a smooth transition from spectra consistent with quiscent to star-forming galaxies. The numbers in each panel indicate the activation value of the two neurons in the latent space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Generative Adversarial Networks (GANs) <a class=\"anchor\" id=\"four\"></a>\n",
    "\n",
    "Finally here's a brief discussion on GANs, since I teased them so much last time. GANs are pure evil, see\n",
    "https://thispersondoesnotexist.com. The picture you see is not a real person-- it was generated by a neural network.\n",
    "\n",
    "But they are also brilliant, incredibly useful, and relatively new (2014). The idea builds logically on autoencoders.\n",
    "1. We have a generator (like the decoder part of an autoencoder) that can produce fake data (e.g., an image). \n",
    "2. Then we have a discriminator (a standard binary classifier) that tries to distinguish fake data from real.  \n",
    "3. Then the generator learns to produce more and more accurate images to trick the discriminator -- without ever seeing any real images -- it just has the feedback from the discriminator.\n",
    "4. We are pitting artifical neural networks against themselves to make them stronger. What could possibly go wrong, right?\n",
    "\n",
    "![](https://learning.oreilly.com/library/view/java-deep-learning/9781788997454/assets/2cf8b4f1-7163-4af1-aa4b-6066329d554a.png)\n",
    "\n",
    "For more, see\n",
    "https://towardsdatascience.com/understanding-generative-adversarial-networks-gans-cd6e4651a29"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Astronomy has seen some interesting uses of autoencoders and GANs in recent years.  For example:\n",
    "    \n",
    "- https://arxiv.org/abs/1702.00403    \n",
    "- https://www.aanda.org/articles/aa/full_html/2017/07/aa30240-16/aa30240-16.html\n",
    "\n",
    "and more generally applying deep networks to GW detection (my main interest):\n",
    "\n",
    "- https://arxiv.org/pdf/1711.03121.pdf\n",
    "- https://arxiv.org/pdf/1812.02183.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'm afraid that is all the time we can spend on deep learning. I'm not saying that you are ready for a data science job after these two lectures, but you now know more than the majority of astrophysicists that talk about neural networks, and you have used the two big codes. There are tons of useful videos, tutorials, and online courses that can take you further if you are interested in deep learning.\n",
    "\n",
    "Another tool you may be interested in: https://astronn.readthedocs.io/en/latest/index.html"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.7.3"
  },
  "nbpresent": {
   "slides": {
    "a6340146-092b-47d0-8584-e84bb64c0952": {
     "id": "a6340146-092b-47d0-8584-e84bb64c0952",
     "prev": null,
     "regions": {
      "9f3c3dc8-1276-4b96-9b97-e4b352cc7fb5": {
       "attrs": {
        "height": 1,
        "width": 1,
        "x": 0.008349570712902259,
        "y": -0.008482103581361025
       },
       "id": "9f3c3dc8-1276-4b96-9b97-e4b352cc7fb5"
      }
     }
    }
   },
   "themes": {}
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}