Chapter 4: Basics of Deep Learning.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Chapter 4 - Basics of Deep Learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import sklearn\n",
    "import torch\n",
    "import torchvision\n",
    "import tqdm\n",
    "\n",
    "# Convenience imports\n",
    "from random import choice \n",
    "from numpy import array, dot, random\n",
    "import sklearn.datasets\n",
    "import sklearn.linear_model\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "import copy\n",
    "import itertools\n",
    "import pandas as pd\n",
    "\n",
    "np.random.seed(1)\n",
    "%matplotlib inline  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. MLP on Various Datasets\n",
    "This example explores a two-layer neural network implemented in python and numpy.\n",
    "Uncomment the dataset to be explored:\n",
    "- Moons Dataset\n",
    "- Circles Dataset\n",
    "- XOR\n",
    "- AND\n",
    "\n",
    "Modified example from: \n",
    "https://github.com/dennybritz/nn-from-scratch"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Generate a dataset and plot it\n",
    "\n",
    "## Moons Dataset ##\n",
    "# X, y = sklearn.datasets.make_moons(200, noise=0.20)\n",
    "\n",
    "## Circles dataset ##\n",
    "# X, y = sklearn.datasets.make_circles(noise=0.2, factor=0.5, random_state=1)\n",
    "\n",
    "## XOR with noise ##\n",
    "xx, yy = np.meshgrid(np.linspace(-3, 3, 50),\n",
    "                     np.linspace(-3, 3, 50))\n",
    "rng = np.random.RandomState(0)\n",
    "X = rng.randn(400, 2)\n",
    "y = np.logical_xor(X[:, 0] > 0, X[:, 1] > 0).astype(int) \n",
    "\n",
    "## AND with noise ##\n",
    "# xx, yy = np.meshgrid(np.linspace(-3, 3, 50),\n",
    "#                      np.linspace(-3, 3, 50))\n",
    "# rng = np.random.RandomState(0)\n",
    "# X = rng.randn(400, 2)\n",
    "# y = np.logical_and(X[:, 0] > 0, X[:, 1] > 0).astype(int) \n",
    "\n",
    "\n",
    "# Plot the dataset\n",
    "plt.scatter(X[:,0], X[:,1], s=40, c=y, cmap=plt.cm.coolwarm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "num_examples = len(X) # training set size\n",
    "nn_input_dim = 2 # input layer dimensionality\n",
    "nn_output_dim = 2 # output layer dimensionality\n",
    "\n",
    "# Gradient descent parameters (I picked these by hand)\n",
    "epsilon = 0.01 # learning rate for gradient descent\n",
    "reg_lambda = 0.01 # regularization strength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function to plot a decision boundary.\n",
    "# If you don't fully understand this function don't worry, it just generates the contour plot below.\n",
    "def plot_decision_boundary(pred_func):\n",
    "    # Set min and max values and give it some padding\n",
    "    x_min, x_max = X[:, 0].min() - .5, X[:, 0].max() + .5\n",
    "    y_min, y_max = X[:, 1].min() - .5, X[:, 1].max() + .5\n",
    "    h = 0.01\n",
    "    # Generate a grid of points with distance h between them\n",
    "    xx, yy = np.meshgrid(np.arange(x_min, x_max, h), np.arange(y_min, y_max, h))\n",
    "    # Predict the function value for the whole gid\n",
    "    Z = pred_func(np.c_[xx.ravel(), yy.ravel()])\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    # Plot the contour and training examples\n",
    "    plt.contourf(xx, yy, Z, cmap=plt.cm.Spectral)\n",
    "    plt.scatter(X[:, 0], X[:, 1], c=y, cmap=plt.cm.Spectral)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function to evaluate the total loss on the dataset\n",
    "def calculate_loss(model):\n",
    "    W1, b1, W2, b2 = model['W1'], model['b1'], model['W2'], model['b2']\n",
    "    # Forward propagation to calculate our predictions\n",
    "    z1 = X.dot(W1) + b1\n",
    "    a1 = np.tanh(z1)\n",
    "    z2 = a1.dot(W2) + b2\n",
    "    exp_scores = np.exp(z2)\n",
    "    probs = exp_scores / np.sum(exp_scores, axis=1, keepdims=True)\n",
    "    # Calculating the loss\n",
    "    corect_logprobs = -np.log(probs[range(num_examples), y])\n",
    "    data_loss = np.sum(corect_logprobs)\n",
    "    # Add regulatization term to loss (optional)\n",
    "    data_loss += reg_lambda/2 * (np.sum(np.square(W1)) + np.sum(np.square(W2)))\n",
    "    return 1./num_examples * data_loss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function to predict an output (0 or 1)\n",
    "def predict(model, x):\n",
    "    W1, b1, W2, b2 = model['W1'], model['b1'], model['W2'], model['b2']\n",
    "    # Forward propagation\n",
    "    z1 = x.dot(W1) + b1\n",
    "    a1 = np.tanh(z1)\n",
    "    z2 = a1.dot(W2) + b2\n",
    "    exp_scores = np.exp(z2)\n",
    "    probs = exp_scores / np.sum(exp_scores, axis=1, keepdims=True)\n",
    "    return np.argmax(probs, axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This function learns parameters for the neural network and returns the model.\n",
    "# - nn_hdim: Number of nodes in the hidden layer\n",
    "# - num_passes: Number of passes through the training data for gradient descent\n",
    "# - print_loss: If True, print the loss every 1000 iterations\n",
    "def build_model(nn_hdim, num_passes=20000, print_loss=False):\n",
    "    \n",
    "    # Initialize the parameters to random values. We need to learn these.\n",
    "    np.random.seed(0)\n",
    "    W1 = np.random.randn(nn_input_dim, nn_hdim) / np.sqrt(nn_input_dim)\n",
    "    b1 = np.zeros((1, nn_hdim))\n",
    "    W2 = np.random.randn(nn_hdim, nn_output_dim) / np.sqrt(nn_hdim)\n",
    "    b2 = np.zeros((1, nn_output_dim))\n",
    "\n",
    "    # This is what we return at the end\n",
    "    model = {}\n",
    "    \n",
    "    # Gradient descent. For each batch...\n",
    "    for i in range(0, num_passes):\n",
    "\n",
    "        # Forward propagation\n",
    "        z1 = X.dot(W1) + b1\n",
    "        a1 = np.tanh(z1)\n",
    "        z2 = a1.dot(W2) + b2\n",
    "        exp_scores = np.exp(z2)\n",
    "        probs = exp_scores / np.sum(exp_scores, axis=1, keepdims=True)\n",
    "\n",
    "        # Backpropagation\n",
    "        delta3 = probs\n",
    "        delta3[range(num_examples), y] -= 1\n",
    "        dW2 = (a1.T).dot(delta3)\n",
    "        db2 = np.sum(delta3, axis=0, keepdims=True)\n",
    "        delta2 = delta3.dot(W2.T) * (1 - np.power(a1, 2))\n",
    "        dW1 = np.dot(X.T, delta2)\n",
    "        db1 = np.sum(delta2, axis=0)\n",
    "\n",
    "        # Add regularization terms (b1 and b2 don't have regularization terms)\n",
    "        dW2 += reg_lambda * W2\n",
    "        dW1 += reg_lambda * W1\n",
    "\n",
    "        # Gradient descent parameter update\n",
    "        W1 += -epsilon * dW1\n",
    "        b1 += -epsilon * db1\n",
    "        W2 += -epsilon * dW2\n",
    "        b2 += -epsilon * db2\n",
    "        \n",
    "        # Assign new parameters to the model\n",
    "        model = { 'W1': W1, 'b1': b1, 'W2': W2, 'b2': b2}\n",
    "        \n",
    "        # Optionally print the loss.\n",
    "        # This is expensive because it uses the whole dataset, so we don't want to do it too often.\n",
    "        if print_loss and i % 1000 == 0:\n",
    "          print(\"Loss after iteration %i: %f\" %(i, calculate_loss(model)))\n",
    "    \n",
    "    return model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Build a model with a 3-dimensional hidden layer\n",
    "model = build_model(3, print_loss=True)\n",
    "\n",
    "# Plot the decision boundary\n",
    "plot_decision_boundary(lambda x: predict(model, x))\n",
    "plt.title(\"Decision Boundary for hidden layer size 3\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Experiment with the hidden layer sizef\n",
    "plt.figure(figsize=(16, 32))\n",
    "hidden_layer_dimensions = [1, 2, 3, 4, 5, 20, 50]\n",
    "for i, nn_hdim in enumerate(hidden_layer_dimensions):\n",
    "    plt.subplot(5, 2, i+1)\n",
    "    plt.title('Hidden Layer size %d' % nn_hdim)\n",
    "    model = build_model(nn_hdim)\n",
    "    plot_decision_boundary(lambda x: predict(model, x))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Spoken Digit Dataset\n",
    "Spoken digits classification and generation.\n",
    "Original dataset can be found here: https://github.com/Jakobovski/free-spoken-digit-dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import torchvision\n",
    "from torch.utils import data\n",
    "import torchaudio\n",
    "from torch.utils.data.sampler import SubsetRandomSampler\n",
    "\n",
    "import torch.utils.data\n",
    "import torch.nn as nn\n",
    "import torch.nn.functional as F\n",
    "import torch.optim as optim\n",
    "from torch.autograd import Variable\n",
    "from torchvision import datasets, transforms\n",
    "from torchvision.utils import make_grid , save_image\n",
    "import os\n",
    "import re\n",
    "import concurrent\n",
    "\n",
    "import numpy as np\n",
    "from scipy.fftpack import fft\n",
    "from scipy import signal\n",
    "from scipy.io import wavfile\n",
    "import matplotlib.pyplot as plt\n",
    "import copy\n",
    "import itertools\n",
    "import pandas as pd\n",
    "from tqdm import tnrange, tqdm_notebook\n",
    "import librosa\n",
    "import librosa.display\n",
    "\n",
    "import IPython.display as ipd\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "use_cuda = torch.cuda.is_available()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download FSDD from Github"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!git clone https://github.com/Jakobovski/free-spoken-digit-dataset.git"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspect the data downloaded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotRawWave(plotTitle, sampleRate, samples, figWidth=14, figHeight=4):\n",
    "    plt.figure(figsize=(figWidth, figHeight))\n",
    "    plt.plot(samples)\n",
    "    plt.title(\"Raw sound wave of \" + plotTitle)\n",
    "    plt.ylabel(\"Amplitude\")\n",
    "    plt.xlabel(\"Time [sec]\")\n",
    "    plt.show()  # force display while in for loop\n",
    "    return None\n",
    "\n",
    "def computeLogSpectrogram(audio, sampleRate, windowSize=20, stepSize=10, epsilon=1e-10):\n",
    "    nperseg  = int(round(windowSize * sampleRate / 1000))\n",
    "    noverlap = int(round(stepSize   * sampleRate / 1000))\n",
    "    freqs, times, spec = signal.spectrogram(audio,\n",
    "                                            fs=sampleRate,\n",
    "                                            window='ham',\n",
    "                                            nperseg=nperseg,\n",
    "                                            noverlap=noverlap,\n",
    "                                            detrend=False)\n",
    "    return freqs, times, np.log(spec.T.astype(np.float32) + epsilon)\n",
    "\n",
    "def plotSpectrogram(plotTitle, freqs, times, spectrogram, figWidth=14, figHeight=4):\n",
    "    fig = plt.figure(figsize=(figWidth, figHeight))\n",
    "    plt.imshow(spectrogram.T, aspect='auto', origin='lower', \n",
    "               cmap=\"jet\",   #  default was \"viridis\"  (perceptually uniform)\n",
    "               extent=[times.min(), times.max(), freqs.min(), freqs.max()])\n",
    "    plt.colorbar(pad=0.01)\n",
    "    plt.title('Spectrogram of ' + plotTitle)\n",
    "    plt.ylabel(\"Frequency [Hz]\")\n",
    "    plt.xlabel(\"Time [sec]\")\n",
    "    fig.tight_layout()\n",
    "    plt.show()  # force display while in for loop\n",
    "    return None\n",
    "\n",
    "def computeMelSpectrogram(samples, sampleRate, nMels=128):\n",
    "    melSpectrum = librosa.feature.melspectrogram(samples.astype(np.float16), sr=sampleRate, n_mels=nMels)\n",
    "    \n",
    "    return melSpectrum\n",
    "\n",
    "def computeLogMelSpectrogram(samples, sampleRate, nMels=128):\n",
    "    melSpectrum = librosa.feature.melspectrogram(samples.astype(np.float16), sr=sampleRate, n_mels=128)\n",
    "    \n",
    "    # Convert to dB, which is a log scale.  Use peak power as reference.\n",
    "    logMelSpectrogram = librosa.power_to_db(melSpectrum, ref=np.max)\n",
    "    \n",
    "    return logMelSpectrogram\n",
    "\n",
    "def plotMelSpectrogram(plotTitle, sampleRate, melSpectrum, figWidth=14, figHeight=4):\n",
    "    fig = plt.figure(figsize=(figWidth, figHeight))\n",
    "    librosa.display.specshow(melSpectrum, sr=sampleRate, x_axis='time', y_axis='mel', cmap=\"jet\")\n",
    "    plt.title('Mel frequency spectrogram: ' + plotTitle)\n",
    "    plt.colorbar(pad=0.01, format='%+02.0f dB')\n",
    "    plt.tight_layout()  \n",
    "    plt.show()  # force display while in for loop\n",
    "    return None\n",
    "\n",
    "def plotLogMelSpectrogram(plotTitle, sampleRate, logMelSpectrum, figWidth=14, figHeight=4):\n",
    "    fig = plt.figure(figsize=(figWidth, figHeight))\n",
    "    librosa.display.specshow(logMelSpectrum, sr=sampleRate, x_axis='time', y_axis='mel', cmap=\"jet\")\n",
    "    plt.title('Mel log-frequency power spectrogram: ' + plotTitle)\n",
    "    plt.colorbar(pad=0.01, format='%+02.0f dB')\n",
    "    plt.tight_layout()  \n",
    "    plt.show()  # force display while in for loop\n",
    "    return None\n",
    "\n",
    "def computeMFCC(samples, sampleRate, nFFT=512, hopLength=256, nMFCC=40, norm=None):\n",
    "    mfcc = librosa.feature.mfcc(y=samples.astype(np.float16), sr=sampleRate, \n",
    "                                n_fft=nFFT, hop_length=hopLength, n_mfcc=nMFCC,\n",
    "                                norm=norm)\n",
    "    \n",
    "    # Let's add on the first and second deltas \n",
    "    mfcc = librosa.feature.delta(mfcc, order=2)\n",
    "    return mfcc\n",
    "\n",
    "def plotMFCC(plotTitle, sampleRate, mfcc, figWidth=14, figHeight=4):\n",
    "    mfcc = sklearn.preprocessing.scale(mfcc, axis=1)\n",
    "    fig = plt.figure(figsize=(figWidth, figHeight))\n",
    "    librosa.display.specshow(mfcc, sr=sampleRate, x_axis='time',cmap=\"jet\")\n",
    "    plt.colorbar(pad=0.01)\n",
    "    plt.title(\"Mel-frequency cepstral coefficients (MFCC): \" + plotTitle)\n",
    "    plt.tight_layout()\n",
    "    plt.show()  # force display while in for loop\n",
    "    return None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def showWavefile(filename):\n",
    "    sampleRate, samples = wavfile.read(filename)  \n",
    "    plotRawWave(filename, sampleRate, samples)\n",
    "    \n",
    "    melSpectrogram = computeMelSpectrogram(samples, sampleRate)\n",
    "    \n",
    "    logMelSpectrogram = computeLogMelSpectrogram(samples, sampleRate)\n",
    "    plotLogMelSpectrogram(filename, sampleRate, logMelSpectrogram)\n",
    "    \n",
    "    return sampleRate, samples, logMelSpectrogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample a file from the dataset\n",
    "filename = 'free-spoken-digit-dataset/recordings/9_theo_16.wav'\n",
    "\n",
    "ipd.display( ipd.Audio(filename) )\n",
    "sampleRate, samples, logMelSpectrogram = showWavefile(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Increase the amount of data: Data Augmentation\n",
    "The original dataset contains 1,500 examples (it has since grown to 2,000 at the time of this writing). We can create additional examples by augmenting the original files. Here we apply time stretching and pitch shifting to increase the number of examples. \n",
    "\n",
    "This has already been computed for convenience, so it is left commented out here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Wav File Paths\n",
    "original_wavfiles=[]\n",
    "\n",
    "root='free-spoken-digit-dataset/recordings/'\n",
    "for path, dirs, files in os.walk(root, topdown=True):\n",
    "    for file in files:\n",
    "        if file.endswith('.wav'):\n",
    "            original_wavfiles.append(os.path.join(path, file))\n",
    "print(\"Number of wavfiles: \", len(original_wavfiles))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "augmented_dir = './data/wav/'\n",
    "os.makedirs(augmented_dir, exist_ok=True)\n",
    "default_sample_rate = 8000 # Reducing everything to 8kHz phone quality for efficiency\n",
    "\n",
    "# Create Augmented examples\n",
    "for file_path in tqdm(original_wavfiles):\n",
    "    y, Fs = librosa.load(file_path) \n",
    "    for ts in [0.75,1,1.25]:\n",
    "        for ps in [-1,0,+1]:\n",
    "            y_new = librosa.effects.time_stretch(y, ts)\n",
    "            y_new = librosa.effects.pitch_shift(y_new, Fs, n_steps=ps)\n",
    "            basename = '_'.join([os.path.splitext(os.path.basename(file_path))[0],str(ts),str(ps)])+'.wav'\n",
    "            output_file = os.path.join(augmented_dir, basename)\n",
    "            librosa.output.write_wav(output_file, y_new, default_sample_rate)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "augmented_wavfiles=[]\n",
    "\n",
    "root='./data/wav/'\n",
    "for path, dirs, files in os.walk(root, topdown=True):\n",
    "    for file in files:\n",
    "        if file.endswith('.wav'):\n",
    "            augmented_wavfiles.append(os.path.join(path, file))\n",
    "print(\"Number of wavfiles: \", len(augmented_wavfiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note: Loading audio and performing the FFT to convert it to the frequence spectrum requires repetative computation. We can save time on this conversion by computing it once and saving the features in numpy arrays, which will be much faster during loading. Furthermore, even doing this process once can take a significant amount of time, so we parallelize the process here. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simplify the network architecture, we fix the size of all examples in the dataset to 1.5s and save them if we need to inspect them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trim_files(file_path, output_dir, max_length):\n",
    "    samples, Fs = librosa.load(file_path) \n",
    "    short_samples = librosa.util.fix_length(samples, 12000)# int(max_length * Fs))\n",
    "\n",
    "    basename = os.path.basename(file_path)\n",
    "    output_file = os.path.join(output_dir, basename)\n",
    "    librosa.output.write_wav(output_file, short_samples, Fs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "executor = concurrent.futures.ProcessPoolExecutor(10)\n",
    "\n",
    "trimmed_dir = './data/trimmed_wav/'\n",
    "max_length = 1.5  # Maximum length in seconds\n",
    "os.makedirs(trimmed_dir, exist_ok=True)\n",
    "futures = [executor.submit(trim_files, file_path, trimmed_dir, max_length) for file_path in augmented_wavfiles]\n",
    "\n",
    "concurrent.futures.wait(futures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trimmed_wavfiles=[]\n",
    "\n",
    "root='./data/trimmed_wav/'\n",
    "for path, dirs, files in os.walk(root, topdown=True):\n",
    "    for file in files:\n",
    "        if file.endswith('.wav'):\n",
    "            trimmed_wavfiles.append(os.path.join(path, file))\n",
    "print(\"Number of wavfiles: \", len(trimmed_wavfiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Save features as numpy arrays for faster loading."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cache_spectrogram_features(file_path, output_dir):\n",
    "    samples, Fs = librosa.load(file_path) \n",
    "\n",
    "    melSpectrum = librosa.feature.melspectrogram(samples.astype(np.float16), sr=Fs, n_mels=40)\n",
    "    \n",
    "    # Convert to dB, which is a log scale.  Use peak power as reference.\n",
    "    logMelSpectrogram = librosa.power_to_db(melSpectrum, ref=np.max)\n",
    "    basename = os.path.splitext(os.path.basename(file_path))[0]+'.npy'\n",
    "    output_file = os.path.join(output_dir, basename)\n",
    "    np.save(output_file, logMelSpectrogram)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "executor = concurrent.futures.ProcessPoolExecutor(10)\n",
    "\n",
    "feature_dir = './data/npy/'\n",
    "os.makedirs(feature_dir, exist_ok=True)\n",
    "futures = [executor.submit(cache_spectrogram_features, file_path, feature_dir) for file_path in trimmed_wavfiles]\n",
    "concurrent.futures.wait(futures)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "npyfiles=[]\n",
    "\n",
    "root='data/npy/'\n",
    "for path, dirs, files in os.walk(root, topdown=True):\n",
    "    for file in files:\n",
    "        if file.endswith('.npy'):\n",
    "            npyfiles.append(os.path.join(path, file))\n",
    "print(\"Number of npy files: \", len(npyfiles))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creat audio dataset splits and data loaders"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class AudioDataset(data.Dataset):\n",
    "    'Characterizes a dataset for PyTorch'\n",
    "    def __init__(self, filenames, labels, sr=8000):\n",
    "        'Initialization'\n",
    "        self.labels = labels\n",
    "        self.filenames = filenames\n",
    "        self.sr = sr\n",
    "\n",
    "    def __len__(self):\n",
    "        'Denotes the total number of samples'\n",
    "        return len(self.filenames)\n",
    "\n",
    "    def __getitem__(self, index):\n",
    "        'Generates one sample of data'\n",
    "        # Select sample\n",
    "        filename = self.filenames[index]\n",
    "        # Load data and get label\n",
    "        X = (torch.Tensor(np.load(filename))+80)/80  # 0 centered [-1,1]\n",
    "        y = self.labels[int(os.path.basename(filename)[0])]  # First character in filename is class\n",
    "        return X, y\n",
    "    \n",
    "    def pad_data(self, s):\n",
    "        padded = np.zeros((self.maxlen,), dtype=np.int64)\n",
    "        if len(s) > self.maxlen: padded[:] = s[:self.maxlen]\n",
    "        else: padded[:len(s)] = s\n",
    "        return padded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "audio_dataset = AudioDataset(npyfiles, list(range(10)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "validation_split = .1\n",
    "testing_split = 0.1\n",
    "shuffle_dataset = True\n",
    "random_seed= 1\n",
    "\n",
    "# Creating data indices for training and validation splits:\n",
    "dataset_size = len(audio_dataset)\n",
    "indices = list(range(dataset_size))\n",
    "\n",
    "if shuffle_dataset :\n",
    "    np.random.seed(random_seed)\n",
    "    np.random.shuffle(indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_indices = indices[:int(-dataset_size*(validation_split+testing_split))]\n",
    "val_indices = indices[int(-dataset_size*(validation_split+testing_split)):int(-dataset_size*testing_split)]\n",
    "test_indices = indices[int(-dataset_size*testing_split):]\n",
    "\n",
    "# # Creating PT data samplers and loaders:\n",
    "train_sampler = SubsetRandomSampler(train_indices)\n",
    "val_sampler = SubsetRandomSampler(val_indices)\n",
    "test_sampler = SubsetRandomSampler(test_indices)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Training examples: ' + str(len(train_sampler)))\n",
    "print('Validation examples: ' + str(len(val_sampler)))\n",
    "print('Testing examples: ' + str(len(test_sampler)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exploratory Data Analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create DataLoaders\n",
    "dataloader_args = dict(batch_size=32, num_workers=8, pin_memory=False)\n",
    "train_loader = torch.utils.data.DataLoader(audio_dataset, **dataloader_args, sampler=train_sampler)\n",
    "val_loader = torch.utils.data.DataLoader(audio_dataset, **dataloader_args, sampler=val_sampler)\n",
    "test_loader = torch.utils.data.DataLoader(audio_dataset, **dataloader_args, sampler=test_sampler)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 1. Let's look at a sample of our data to make sure that it looks reasonable. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def imshow(img, title=''):\n",
    "    img = img # unnormalize\n",
    "    npimg = img.numpy()\n",
    "    plt.title(title) if title is not None else plt.title(str(c).zfill(3))\n",
    "    plt.imshow(np.transpose(npimg, (1,2,0)),interpolation='nearest', aspect='auto',origin='lower')\n",
    "\n",
    "dataiter = iter(train_loader)\n",
    "X, y = dataiter.next()\n",
    "X = (X+40)/80\n",
    "# X = X / torch.Tensor([40]).expand_as(X)\n",
    "grid_img = torchvision.utils.make_grid(X.unsqueeze(1), normalize=True)\n",
    "imshow(grid_img, title='Sample from Spoken Digits')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 2. Let's examine the balance between the classes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_labels(loader):\n",
    "    return np.asarray([label for batch in loader for label in batch[1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_labels = get_labels(train_loader)\n",
    "val_labels = get_labels(val_loader)\n",
    "test_labels = get_labels(test_loader)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show_hist(data, title):\n",
    "    d = np.diff(np.unique(data)).min()\n",
    "    left_of_first_bin = data.min() - float(d)/2\n",
    "    right_of_last_bin = data.max() + float(d)/2\n",
    "    plt.hist(data, np.arange(left_of_first_bin, right_of_last_bin + d, d), rwidth=0.5)\n",
    "    plt.grid(True)\n",
    "    plt.xticks(range(0, 10))\n",
    "    plt.title(title)\n",
    "    plt.ylabel(\"Number of examples\")\n",
    "    plt.xlabel(\"Example label\")\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show_hist(train_labels, \"Histogram for Training Data\")\n",
    "show_hist(val_labels, \"Histogram for Validation Data\")\n",
    "show_hist(test_labels, \"Histogram for Testing Data\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_features = 960 # 40 logmel features * 24 time steps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Training a Supervised Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define model architecture in pytorch\n",
    "class Model(nn.Module):\n",
    "    def __init__(self, input_size):\n",
    "        super(Model, self).__init__()\n",
    "        self.input_size = input_size\n",
    "        \n",
    "        self.fc1 = nn.Linear(self.input_size, 128)\n",
    "        self.bc1 = nn.BatchNorm1d(128)\n",
    "        \n",
    "        self.fc2 = nn.Linear(128, 128)\n",
    "        self.bc2 = nn.BatchNorm1d(128)\n",
    "        \n",
    "        self.fc3 = nn.Linear(128, 10)\n",
    "        \n",
    "        \n",
    "    def forward(self, x):\n",
    "        x = x.view((-1, self.input_size))\n",
    "        h = self.fc1(x)\n",
    "        h = self.bc1(h)\n",
    "        h = torch.relu(h)\n",
    "        h = F.dropout(h, p=0.5, training=self.training)\n",
    "        \n",
    "        h = self.fc2(h)\n",
    "        h = self.bc2(h)\n",
    "        h = torch.relu(h)\n",
    "        h = F.dropout(h, p=0.5, training=self.training)\n",
    "        \n",
    "        h = self.fc3(h)\n",
    "        out = torch.log_softmax(h,dim=1)\n",
    "        return out\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model = Model(n_features)\n",
    "# Track metrics when training\n",
    "optimizer = optim.Adam(model.parameters(), lr=0.01)\n",
    "n_epoch = 40\n",
    "train_losses = []\n",
    "val_losses = []\n",
    "best_val = np.inf\n",
    "best_model = copy.deepcopy(model)\n",
    "accuracies = []\n",
    "\n",
    "device = torch.device(\"cuda\" if use_cuda else \"cpu\")\n",
    "model.to(device)\n",
    "\n",
    "# Train the model\n",
    "for epoch in range(n_epoch):\n",
    "    print('Epoch: ' + str(epoch + 1)) #Index at 0\n",
    "    # Training loop\n",
    "    model.train()\n",
    "    correct = 0\n",
    "    train_loss = 0\n",
    "    with tqdm_notebook(enumerate(train_loader), total=len(train_loader)) as pbar:\n",
    "        for batch_idx, (data, target) in pbar:\n",
    "            # Get Samples\n",
    "            data, target = data.to(device), target.to(device)\n",
    "\n",
    "            # Clear gradients\n",
    "            optimizer.zero_grad()\n",
    "\n",
    "            # Forward Propagation \n",
    "            y_pred = model(data) \n",
    "\n",
    "            # Error Computation\n",
    "            loss = F.cross_entropy(y_pred, target)\n",
    "\n",
    "            # Backpropagation\n",
    "            loss.backward()\n",
    "            optimizer.step()\n",
    "\n",
    "            # Track training loss\n",
    "            train_loss += loss.data.sum()\n",
    "            \n",
    "            #Added to look at training accuracy\n",
    "            pred = y_pred.argmax(dim=1, keepdim=True) # get the index of the max log-probability\n",
    "            correct += pred.eq(target.view_as(pred)).sum().item()\n",
    "            \n",
    "            pbar.set_description(\"Current loss %.4f\" % (train_loss/(len(target)*(1+batch_idx))))\n",
    "    train_losses.append(train_loss/len(train_loader.sampler.indices))\n",
    "    print('Training set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\\n'.format(\n",
    "        train_losses[-1], correct, len(train_loader.sampler.indices),\n",
    "        100. * correct / len(train_loader.sampler.indices)))\n",
    "    \n",
    "    # Validation Loop\n",
    "    model.eval()\n",
    "    val_loss = 0\n",
    "    correct = 0\n",
    "    with torch.no_grad():\n",
    "        for batch_idx, (data, target) in enumerate(val_loader):\n",
    "            data, target = data.to(device), target.to(device)\n",
    "            output = model(data)\n",
    "            val_loss += F.nll_loss(output, target, reduction='sum').item() # sum up batch loss\n",
    "            pred = output.argmax(dim=1, keepdim=True) # get the index of the max log-probability\n",
    "            correct += pred.eq(target.view_as(pred)).sum().item()\n",
    "\n",
    "    val_loss /= len(val_loader.sampler.indices)\n",
    "    val_losses.append(val_loss)\n",
    "    if val_loss < best_val:\n",
    "        best_val = val_loss\n",
    "        best_model = copy.deepcopy(model)\n",
    "    accuracy = 100. * correct / len(val_loader.sampler.indices)\n",
    "    accuracies.append(accuracy)\n",
    "    print('\\nValidation set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\\n'.format(\n",
    "        val_loss, correct, len(val_indices), accuracy))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_learning_curves(train_losses, val_losses, accuracies):\n",
    "    plt.figure()\n",
    "    fig, ax1 = plt.subplots()\n",
    "    ax2 = ax1.twinx()\n",
    "    \n",
    "    ax1.set_ylabel(\"Loss\")\n",
    "    ax1.plot(range(len(train_losses)), train_losses, label='Training Loss')\n",
    "    ax1.plot(range(len(train_losses)), val_losses, label='Validation Loss')\n",
    "    ax2.plot(range(len(accuracies)), accuracies, label='Validation Accuracy', color='g')\n",
    "    ax2.set_ylabel('Accuracy')\n",
    "    \n",
    "    plt.title(\"Training/Validation Curves\")\n",
    "    plt.legend(loc=\"best\")\n",
    "    plt.xlabel(\"Epoch\")\n",
    "    \n",
    "    handles, labels = ax1.get_legend_handles_labels()\n",
    "    handles2, labels2 = ax2.get_legend_handles_labels()\n",
    "    ax2.legend(handles+handles2, labels+labels2, bbox_to_anchor=(1.1, 1), loc=2, borderaxespad=0.)\n",
    "    plt.grid()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_learning_curves(train_losses, val_losses, accuracies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Test Set\n",
    "model.eval()\n",
    "val_loss = 0\n",
    "correct = 0\n",
    "with torch.no_grad():\n",
    "    for batch_idx, (data, target) in enumerate(test_loader):\n",
    "        data, target = data.to(device), target.to(device)\n",
    "        output = model(data)\n",
    "        val_loss += F.nll_loss(output, target, reduction='sum').item() # sum up batch loss\n",
    "        pred = output.argmax(dim=1, keepdim=True) # get the index of the max log-probability\n",
    "        correct += pred.eq(target.view_as(pred)).sum().item()\n",
    "\n",
    "val_loss /= len(val_loader.sampler.indices)\n",
    "# val_losses.append(val_loss)\n",
    "if val_loss < best_val:\n",
    "    best_val = val_loss\n",
    "    best_model = copy.deepcopy(model)\n",
    "accuracy = 100. * correct / len(val_loader.sampler.indices)\n",
    "# accuracies.append(accuracy)\n",
    "print('\\nTest set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\\n'.format(\n",
    "    val_loss, correct, len(val_indices), accuracy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Restricted Boltzmann Machine\n",
    "Train an RBM to generate examples from the MNIST dataset.\n",
    "\n",
    "Modificaiton of the example shown: https://github.com/odie2630463/Restricted-Boltzmann-Machines-in-pytorch/blob/master/RBM.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def show(file_name,img):\n",
    "    npimg = np.transpose(img.numpy(),(2,1,0))\n",
    "    f = \"./%s.png\" % file_name\n",
    "    plt.imshow(npimg)\n",
    "#     plt.imsave(f,npimg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class RBM(nn.Module):\n",
    "    def __init__(self,\n",
    "                 n_vis=256,\n",
    "                 n_hin=500,\n",
    "                 k=5):\n",
    "        super(RBM, self).__init__()\n",
    "        self.W = nn.Parameter(torch.randn(n_hin,n_vis)*1e-2).to(device)\n",
    "        self.v_bias = nn.Parameter(torch.zeros(n_vis))\n",
    "        self.h_bias = nn.Parameter(torch.zeros(n_hin))\n",
    "        self.k = k\n",
    "    \n",
    "    def sample_from_p(self,p):\n",
    "        return F.relu(torch.sign(p - torch.rand(p.size()).to(device)))\n",
    "    \n",
    "    def v_to_h(self,v):\n",
    "        p_h = torch.sigmoid(F.linear(v.to(device),self.W.to(device),self.h_bias.to(device)))\n",
    "        sample_h = self.sample_from_p(p_h)\n",
    "        return p_h,sample_h\n",
    "    \n",
    "    def h_to_v(self,h):\n",
    "        p_v = torch.sigmoid(F.linear(h.to(device),self.W.t().to(device),self.v_bias.to(device)))\n",
    "        sample_v = self.sample_from_p(p_v)\n",
    "        return p_v,sample_v\n",
    "        \n",
    "    def forward(self,v):\n",
    "        pre_h1,h1 = self.v_to_h(v)\n",
    "        \n",
    "        h_ = h1\n",
    "        for _ in range(self.k):\n",
    "            pre_v_,v_ = self.h_to_v(h_)\n",
    "            pre_h_,h_ = self.v_to_h(v_)\n",
    "        \n",
    "        return v,v_\n",
    "    \n",
    "    def free_energy(self,v):\n",
    "        vbias_term = v.mv(self.v_bias)\n",
    "        wx_b = F.linear(v,self.W,self.h_bias)\n",
    "        hidden_term = wx_b.exp().add(1).log().sum(1)\n",
    "        return (-hidden_term - vbias_term).mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbm = RBM(n_vis=n_features, n_hin=128, k=1)\n",
    "device = torch.device(\"cuda\" if use_cuda else \"cpu\")\n",
    "rbm.to(device)\n",
    "\n",
    "train_op = optim.Adam(rbm.parameters(), 0.001)\n",
    "\n",
    "num_epochs = 20\n",
    "\n",
    "for epoch in range(num_epochs):\n",
    "    loss_ = []\n",
    "    total_loss = 0\n",
    "    for _, (data,target) in enumerate(train_loader):\n",
    "        data = Variable(data.view(-1,n_features)).to(device)\n",
    "        sample_data = torch.bernoulli(data)\n",
    "        \n",
    "        v,v1 = rbm(sample_data)\n",
    "        loss = rbm.free_energy(v) - rbm.free_energy(v1)\n",
    "#         total_loss += loss\n",
    "        train_op.zero_grad()\n",
    "        loss.backward()\n",
    "        train_op.step()\n",
    "    print('epoch [{}/{}], loss: {:.4f}'.format(epoch + 1, num_epochs, torch.mean(loss)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Show real and generated features that have been learned by the RBM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(\"real\",make_grid(v.view(32, 1, 40, 24).cpu()))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "show(\"generate\",make_grid(v1.view(32, 1, 40, 24).cpu().data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes,\n",
    "                          normalize=False,\n",
    "                          title='Confusion matrix',\n",
    "                          cmap=plt.cm.Blues):\n",
    "    \"\"\"\n",
    "    This function prints and plots the confusion matrix.\n",
    "    Normalization can be applied by setting `normalize=True`.\n",
    "    \"\"\"\n",
    "    if normalize:\n",
    "        cm = cm.astype('float') / cm.sum(axis=1)[:, np.newaxis]\n",
    "        print(\"Normalized confusion matrix\")\n",
    "    else:\n",
    "        print('Confusion matrix, without normalization')\n",
    "\n",
    "    print(cm)\n",
    "\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=45)\n",
    "    plt.yticks(tick_marks, classes)\n",
    "\n",
    "    fmt = '.2f' if normalize else 'd'\n",
    "    thresh = cm.max() / 2.\n",
    "    for i, j in itertools.product(range(cm.shape[0]), range(cm.shape[1])):\n",
    "        plt.text(j, i, format(cm[i, j], fmt),\n",
    "                 horizontalalignment=\"center\",\n",
    "                 color=\"white\" if cm[i, j] > thresh else \"black\")\n",
    "\n",
    "    plt.tight_layout()\n",
    "    plt.ylabel('True label')\n",
    "    plt.xlabel('Predicted label')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Extract features from our training set\n",
    "batch_size = train_loader.batch_size\n",
    "output_size = 128\n",
    "input_size = n_features\n",
    "\n",
    "train_features = np.zeros((len(train_loader)*batch_size, output_size))\n",
    "train_labels = np.zeros(len(train_loader)*batch_size)\n",
    "test_features = np.zeros((len(test_loader)*batch_size, output_size))\n",
    "test_labels = np.zeros(len(test_loader)*batch_size)\n",
    "\n",
    "for i, (batch,labels) in enumerate(train_loader):\n",
    "    batch = Variable(batch.view(-1,input_size))  # flatten input data\n",
    "\n",
    "#     train_features[i*batch_size:i*batch_size+len(batch)] = rbm_2.v_to_h(rbm_1.v_to_h(batch)[0])[0].cpu().data.numpy()\n",
    "    train_features[i*batch_size:i*batch_size+len(batch)] = rbm.v_to_h(batch)[0].cpu().data.numpy()\n",
    "#     train_features[i*batch_size:i*batch_size+len(batch)] = model(batch)[0].cpu().data.numpy()\n",
    "    train_labels[i*batch_size:(i+1)*batch_size] = labels.numpy()\n",
    "\n",
    "for i, (batch, labels) in enumerate(test_loader):\n",
    "    batch = Variable(batch.view(-1,input_size))  # flatten input data\n",
    "\n",
    "#     test_features[i*batch_size:i*batch_size+len(batch)] = rbm_2.v_to_h(rbm_1.v_to_h(batch)[0])[0].cpu().data.numpy()\n",
    "    test_features[i*batch_size:i*batch_size+len(batch)] = rbm.v_to_h(batch)[0].cpu().data.numpy()\n",
    "#     test_features[i*batch_size:i*batch_size+len(batch)] = model(batch)[0].cpu().data.numpy()\n",
    "    test_labels[i*batch_size:i*batch_size+len(batch)] = labels.numpy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Train a Logistic Regression classifier on the extracted features\n",
    "clf = LogisticRegression()\n",
    "clf.fit(train_features, train_labels)\n",
    "predictions = clf.predict(test_features)\n",
    "\n",
    "print('Result: %d/%d' % (sum(predictions == test_labels), test_labels.shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compute confusion matrix\n",
    "cnf_matrix = sklearn.metrics.confusion_matrix(predictions, test_labels)\n",
    "np.set_printoptions(precision=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display confusion matrix\n",
    "plt.figure()\n",
    "class_names = [0,1,2,3,4,5,6,7,8,9]\n",
    "plot_confusion_matrix(cnf_matrix, classes=class_names,\n",
    "                      title='Confusion matrix, without normalization')\n",
    "\n",
    "# Plot normalized confusion matrix\n",
    "plt.figure()\n",
    "plot_confusion_matrix(cnf_matrix, classes=class_names, normalize=True,\n",
    "                      title='Normalized confusion matrix')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "#  Autoencoder\n",
    "Train an autoencoder to learn a compressed representation of the spoken digits.\n",
    "\n",
    "Modification of example shown: https://github.com/SherlockLiao/pytorch-beginner/tree/master/08-AutoEncoder"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PyTorch model definition \n",
    "class autoencoder(nn.Module):\n",
    "    def __init__(self, input_size):\n",
    "        super(autoencoder, self).__init__()\n",
    "        self.input_size = input_size\n",
    "        \n",
    "        self.e_fc1 = nn.Linear(self.input_size, 512)        \n",
    "        self.e_fc2 = nn.Linear(512, 128)\n",
    "        self.e_fc3 = nn.Linear(128, 64)\n",
    "        self.e_fc4 = nn.Linear(64,64)\n",
    "        \n",
    "        self.d_fc1 = nn.Linear(64, 64)\n",
    "        self.d_fc2 = nn.Linear(64, 128)\n",
    "        self.d_fc3 = nn.Linear(128, 512)\n",
    "        self.d_fc4 = nn.Linear(512, self.input_size)\n",
    "    \n",
    "    def forward(self, x):\n",
    "        # Encoder\n",
    "        h = F.relu(self.e_fc1(x))\n",
    "        h = F.relu(self.e_fc2(h))\n",
    "        h = F.relu(self.e_fc3(h))\n",
    "        h = self.e_fc4(h)\n",
    "        \n",
    "        # Decoder\n",
    "        h = F.relu(self.d_fc1(h))\n",
    "        h = F.relu(self.d_fc2(h))\n",
    "        h = F.relu(self.d_fc3(h))\n",
    "        h = self.d_fc4(h)\n",
    "        out = torch.tanh(h)\n",
    "\n",
    "        return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert tensor to image\n",
    "def to_img(x):\n",
    "    x = 0.5 * (x + 1)\n",
    "    x = x.clamp(0, 1)\n",
    "    x = x.view(x.size(0), 1, 40, 24)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "num_epochs = 20\n",
    "batch_size = 128\n",
    "learning_rate = 0.001\n",
    "\n",
    "model = autoencoder(n_features)\n",
    "device = torch.device(\"cuda\" if use_cuda else \"cpu\")\n",
    "model.to(device)\n",
    "optimizer = torch.optim.Adam(\n",
    "    model.parameters(), lr=learning_rate, weight_decay=1e-5)\n",
    "save_dir = './autoencoder_imgs/'\n",
    "os.makedirs(save_dir, exist_ok=True)\n",
    "\n",
    "\n",
    "for epoch in range(num_epochs):\n",
    "    for data, target in train_loader:\n",
    "        # Get samples\n",
    "        input = data.view(-1,n_features).to(device)  # We will reuse the formatted input as our target\n",
    "        \n",
    "        # Forward Propagation\n",
    "        output = model(input)\n",
    "\n",
    "        # Error Computation\n",
    "        loss = F.mse_loss(output, input)\n",
    "\n",
    "        # Clear gradients\n",
    "        optimizer.zero_grad()\n",
    "\n",
    "        # Backpropagation\n",
    "        loss.backward()\n",
    "\n",
    "        # Parameter Update\n",
    "        optimizer.step()\n",
    "    print('epoch [{}/{}], loss:{:.4f}'\n",
    "          .format(epoch + 1, num_epochs, loss))\n",
    "    if epoch % 10 == 0:\n",
    "        pic = to_img(output.cpu().data)\n",
    "        grid_img = torchvision.utils.make_grid(pic, normalize=True)\n",
    "        imshow(grid_img, title='Autoencoder Reconstructed Output from Spoken Digits')\n",
    "        save_image(pic, os.path.join(save_dir, 'image_{}.png'.format(epoch)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}