tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# CNN tutorial\n",
    "\n",
    "This tutorial will introduce a simple Convolutional Neural Network (CNN) architecture to classify images.\n",
    "\n",
    "As an example we will use the UC Merced dataset ( http://weegee.vision.ucmerced.edu/datasets/landuse.html ).\n",
    "It consists of 21 categories, with 100 images each. The RGB images will be resized to `256 x 256` pixels.\n",
    "\n",
    "## Tasks\n",
    "\n",
    "1. Go through the code and run cell by cell (using the \"Run\" button or \"shift+enter\"). \n",
    "2. Complete the CNN architecture\n",
    "3. Answer the questions\n",
    "\n",
    "**Windows users:** Adjust the relative paths for your system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import itertools\n",
    "\n",
    "from sklearn.metrics import confusion_matrix, classification_report\n",
    "\n",
    "import keras\n",
    "from keras.models import Sequential, Model, load_model\n",
    "from keras.layers import Input, Dense, Conv2D, MaxPooling2D, GlobalAveragePooling2D, AveragePooling2D, ZeroPadding2D, Dropout, Flatten, add, Reshape, Activation\n",
    "from keras.layers.normalization import BatchNormalization\n",
    "from keras.optimizers import Adam, SGD\n",
    "from keras.utils import np_utils, plot_model\n",
    "from keras import backend as K\n",
    "from keras import initializers\n",
    "from keras.preprocessing.image import ImageDataGenerator\n",
    "from keras.callbacks import ModelCheckpoint\n",
    "\n",
    "import os\n",
    "from shutil import copyfile\n",
    "\n",
    "# optional to load the original image files instead of numpy files\n",
    "from PIL import Image\n",
    "\n",
    "\n",
    "\n",
    "# set random seed\n",
    "seed = 4321"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the dataset\n",
    "\n",
    "Set the (relative) `data_path` on your system.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data source: http://weegee.vision.ucmerced.edu/datasets/landuse.html\n",
    "\n",
    "category_names = [\n",
    " 'agricultural',\n",
    " 'airplane',\n",
    " 'baseballdiamond',\n",
    " 'beach',\n",
    " 'buildings',\n",
    " 'chaparral',\n",
    " 'denseresidential',\n",
    " 'forest',\n",
    " 'freeway',\n",
    " 'golfcourse',\n",
    " 'harbor',\n",
    " 'intersection',\n",
    " 'mediumresidential',\n",
    " 'mobilehomepark',\n",
    " 'overpass',\n",
    " 'parkinglot',\n",
    " 'river',\n",
    " 'runway',\n",
    " 'sparseresidential',\n",
    " 'storagetanks',\n",
    " 'tenniscourt']\n",
    "\n",
    "\n",
    "# load all data\n",
    "N_images = 100  # per category\n",
    "nb_classes = 21\n",
    "patch_size = 256  # width and height in pixel\n",
    "channels = 3\n",
    "\n",
    "'''\n",
    "# Load the original dataset\n",
    "data_dir = 'data/UCMerced_LandUse/Images'\n",
    "images_orig = []\n",
    "labels_num = []\n",
    "for i in range(len(category_names)):\n",
    "    \n",
    "    # load all images per class\n",
    "    for j in range(N_images):        \n",
    "        img_path = os.path.join(data_dir, category_names[i], category_names[i]+'{:02d}.tif'.format(j))\n",
    "        \n",
    "        img = Image.open(img_path)\n",
    "        \n",
    "        img = img.resize((patch_size, patch_size), Image.ANTIALIAS)\n",
    "        arr = np.asarray(img)\n",
    "        images_orig.append(arr)\n",
    "        labels_num.append([i])\n",
    "\n",
    "images_orig = np.array(images_orig)\n",
    "labels_num = np.array(labels_num)\n",
    "\n",
    "# Save the numpy datasets\n",
    "np.save('data/UCMerced_images_orig.npy', images_orig)\n",
    "np.save('data/UCMerced_labels_num.npy', labels_num)\n",
    "'''\n",
    "\n",
    "# Load the numpy datasets\n",
    "images_orig = np.load('data/UCMerced_images_orig.npy')\n",
    "labels_num = np.load('data/UCMerced_labels_num.npy')\n",
    "\n",
    "# Convert labels to categorical 'one hot encoded vector' --> required format for softmax cross-entropy loss.\n",
    "# Create an array of zeros with length nb_classes and put a 1 at the index of the true label.\n",
    "labels = np_utils.to_categorical(labels_num, nb_classes)\n",
    "print('example of a label with an integer coding: \\n{}'.format(labels_num[333]))\n",
    "print('the same example as a one-hot encoded: vector: \\n{}'.format(labels[333]))\n",
    "\n",
    "print(images_orig.shape, images_orig.dtype)\n",
    "print(labels.shape, labels.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Organize the raw data in new directories train, val, test\n",
    "\n",
    "For training we want to load the data from individual directories.\n",
    "\n",
    "Instead of running this you can download the splitted data directories from the Google drive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "copy_data = False\n",
    "\n",
    "\n",
    "# Collect all image paths\n",
    "data_dir = 'data/UCMerced_LandUse/Images'\n",
    "image_paths = []\n",
    "for i in range(len(category_names)):\n",
    "\n",
    "    # load all images per class\n",
    "    for j in range(N_images):        \n",
    "        img_path = os.path.join(data_dir, category_names[i], category_names[i]+'{:02d}.tif'.format(j))\n",
    "\n",
    "        image_paths.append(img_path)\n",
    "\n",
    "# Split data into training and test\n",
    "image_paths = np.array(image_paths)\n",
    "\n",
    "nb_images = image_paths.shape[0]\n",
    "shuffled_indices = np.arange(nb_images)\n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(shuffled_indices)\n",
    "\n",
    "training_indices = shuffled_indices[:int(0.7*nb_images)]               # 70% for training\n",
    "val_indices = shuffled_indices[int(0.7*nb_images):int(0.8*nb_images)]  # 10% for val\n",
    "test_indices = shuffled_indices[int(0.8*nb_images):]                   # 20% for test\n",
    "\n",
    "train_paths = image_paths[training_indices]\n",
    "val_paths = image_paths[val_indices]\n",
    "test_paths = image_paths[test_indices]\n",
    "\n",
    "# Copy images to train, val, test directories\n",
    "\n",
    "# make directories:\n",
    "for split in ['train', 'val', 'test']:\n",
    "    if not os.path.exists('data/UCMerced_LandUse/'+split):\n",
    "        os.makedirs('data/UCMerced_LandUse/'+split)\n",
    "\n",
    "def copy_images_to_split_directory(paths, split):\n",
    "    for p in paths:\n",
    "        dst=p.replace('UCMerced_LandUse/Images', 'UCMerced_LandUse/'+split)\n",
    "        if not os.path.exists(os.path.dirname(dst)):\n",
    "            os.makedirs(os.path.dirname(dst))\n",
    "        copyfile(src=p, dst=dst)\n",
    "\n",
    "if copy_data:\n",
    "    copy_images_to_split_directory(train_paths, split='train')\n",
    "    copy_images_to_split_directory(val_paths, split='val')\n",
    "    copy_images_to_split_directory(test_paths, split='test')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize the images\n",
    "\n",
    "**Task**:\n",
    "* Which categories might be easy to classify for a CNN? Note three and explain why.\n",
    "* Which category-pairs might be confused? Note three and explain why."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# VISUALIZE IMAGES\n",
    "def plotImages( images, n_images=5):\n",
    "    fig, axes = plt.subplots(1, n_images, figsize=(10, 10))\n",
    "    axes = axes.flatten()\n",
    "    for img, ax in zip(images, axes):\n",
    "        ax.imshow(img)\n",
    "        ax.set_xticks(())\n",
    "        ax.set_yticks(())\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "# Plot some examples for each category    \n",
    "def plot_examples_per_category():\n",
    "    for i in range(nb_classes):\n",
    "        print(category_names[i])\n",
    "        images_to_print = images_orig[i*N_images:i*N_images+N_images]\n",
    "        np.random.seed(seed)\n",
    "        np.random.shuffle(images_to_print)\n",
    "        print(images_to_print.shape)\n",
    "        plotImages(images_to_print)\n",
    "\n",
    "plot_examples_per_category()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN Architecture\n",
    "\n",
    "**Task**: Implement the following convolutional neural network (CNN). Consider the information in the Figure as well as the additional remarks:\n",
    "\n",
    "* Use `strides = (1, 1)` for all convolutional layers. \n",
    "* Use `padding='same'` to pad the inputs and preserve the spatial dimensions after the convolution.\n",
    "* Use 'relu' as activation for both conv and dense layers, except the last layer that should output 'softmax' activations (pseudo class probabilities).\n",
    "* The spatial extent is reduced by using max pooling with `pool_size = (3, 3)`. Set the parameter `strides`, such that the spatial extent is reduced by factor 2 in both spatial dimensions.\n",
    "* Add a dropout layer with `prob_drop_hidden = 0.5` after the first dense layer.\n",
    "\n",
    "![Figure simple CNN](figures/simple_CNN.png)\n",
    "*Figure 1: Visualization of a simple CNN architecture*\n",
    "\n",
    "\n",
    "**Questions**:\n",
    "1. What is the shape of a single filter kernel in the second convolutional layer?\n",
    "2. What is the shape of the output feature maps after the third convolutional layer (before the MaxPooling)?\n",
    "3. `model.summary()` shows the dimension of each layer output. What is the first dimension (None) of our tensors corresponding to?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# parameters:\n",
    "input_shape = (patch_size, patch_size, channels)\n",
    "pool_size = (3, 3)                  # size of pooling area for max pooling\n",
    "prob_drop_hidden = 0.5              # drop probability for dropout @ fc layer\n",
    "\n",
    "def simple_CNN():\n",
    "\n",
    "    model = Sequential()\n",
    "\n",
    "    # conv1 layer\n",
    "    # NOTE: the input shape is only needed for the first layer\n",
    "    model.add(Conv2D(filters=64, kernel_size=(3, 3), strides=(1,1), padding='same', activation='relu', \n",
    "                   input_shape=input_shape))   \n",
    "    \n",
    "    model.add(MaxPooling2D(pool_size=pool_size, strides=(TODO,TODO), padding='same'))\n",
    "    \n",
    "\n",
    "    # TODO ...\n",
    "    \n",
    "    \n",
    "    # reduce the spatial dimensions to 1x1 by averaging the activations for each depth\n",
    "    model.add(GlobalAveragePooling2D())\n",
    "    \n",
    "    # fc1 layer\n",
    "    model.add(Dense(128, activation='relu'))\n",
    "    model.add(Dropout(rate=prob_drop_hidden))\n",
    "\n",
    "    # The output of the model (softmax)\n",
    "    \n",
    "    # TODO ...\n",
    "\n",
    "    return model\n",
    "\n",
    "\n",
    "# clear the session first, such that layer names start with index 1\n",
    "K.clear_session()\n",
    "\n",
    "model_simple_CNN = simple_CNN()\n",
    "model_simple_CNN.summary()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We could also plot the model graph\n",
    "# Note: You might need further installations\n",
    "\n",
    "# plot_model(model_simple_CNN, show_shapes=True, to_file='simple_CNN.png')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Some helper functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A custom Data generator to iterate over the training data and load the samples as batches\n",
    "# (We will not use it today...)\n",
    "def generator(features, labels, batch_size, patch_size, channels, nb_classes):\n",
    " # Create empty arrays to contain batch of features and labels#\n",
    " batch_features = np.zeros((batch_size, patch_size, patch_size, channels))\n",
    " batch_labels = np.zeros((batch_size, nb_classes))\n",
    "\n",
    " while True:\n",
    "    for i in range(batch_size):\n",
    "         # choose random index in features\n",
    "         index = np.random.choice(a=features.shape[0], size=1)\n",
    "\n",
    "         #Optional: Add some data augmentation here\n",
    "\n",
    "         batch_features[i] = features[index, :, :, :]\n",
    "         batch_labels[i] = labels[index]\n",
    "    yield batch_features, batch_labels\n",
    "    \n",
    "    \n",
    "# Function to re-initialize the weights before starting a new experiment/training\n",
    "def initialize_weights(model, layer_name=None):\n",
    "    session = K.get_session()\n",
    "    if layer_name is None:\n",
    "        for layer in model.layers: \n",
    "            if hasattr(layer, 'kernel_initializer'):\n",
    "                print('initialize weights of layer: {}'.format(layer.name))\n",
    "                layer.kernel.initializer.run(session=session)\n",
    "    else:\n",
    "        layer = model.get_layer(name=layer_name)\n",
    "        if hasattr(layer, 'kernel_initializer'):\n",
    "            print('initialize weights of layer: {}'.format(layer.name))\n",
    "            layer.kernel.initializer.run(session=session)\n",
    "            \n",
    "            \n",
    "# Helper function to plot training curves\n",
    "def plot_train_val_accuracy_and_loss(history, legend_suffix=\"\", figsize=(8,6)):\n",
    "    # list all data in history\n",
    "    print(history.keys())\n",
    "\n",
    "    # summarize history for accuracy\n",
    "    plt.figure(figsize=figsize)\n",
    "    plt.plot(history['acc'])\n",
    "    plt.plot(history['val_acc'])\n",
    "    plt.title('Accuracy')\n",
    "    plt.ylabel('accuracy')\n",
    "    plt.xlabel('epoch')\n",
    "    plt.legend(['train_'+legend_suffix, 'val_'+legend_suffix], loc='lower right')\n",
    "    plt.show()\n",
    "\n",
    "    # summarize history for loss\n",
    "    plt.figure(figsize=figsize)\n",
    "    plt.plot(history['loss'])\n",
    "    plt.plot(history['val_loss'])\n",
    "    plt.title('Loss')\n",
    "    plt.ylabel('loss')\n",
    "    plt.xlabel('epoch')\n",
    "    plt.legend(['train_'+legend_suffix, 'val_'+legend_suffix], loc='upper right')\n",
    "    plt.show()            "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocess the input images\n",
    "We zero-center and normalize each input feature (RGB channel) separately.\n",
    "Thus, we substract the training mean and divide by the training standard deviation for each RGB channel separately. \n",
    "\n",
    "You don't have to run this cell. The training mean and std are already provided."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# PREPROCESSING: NORMALIZE THE INPUT DATA\n",
    "# as calculating the Mean and STD over the entire training data takes some time, the values are provided.\n",
    "\n",
    "run_preprocessing = False\n",
    "\n",
    "\n",
    "def get_mean_and_std_per_channel(data_train):\n",
    "   \n",
    "    if not data_train.dtype == np.float64:\n",
    "        print('CONVERTING: to float64 for preprocessing')\n",
    "        data_train = np.asarray(data_train, dtype=np.float64)\n",
    "\n",
    "    mean = np.nanmean(data_train, axis=(0, 1, 2))\n",
    "    std = np.nanstd(data_train, axis=(0, 1, 2))\n",
    "\n",
    "    print('mean: ', mean)\n",
    "    print('mean.shape', mean.shape, 'mean.dtype: ', mean.dtype)\n",
    "\n",
    "    print('std: ', std)\n",
    "    print('std.shape', std.shape, 'std.dtype: ', std.dtype)\n",
    "        \n",
    "    data_norm = normalize_images_per_channel(data_train, mean, std)\n",
    "    \n",
    "    # Control mean and std after preprocessing\n",
    "    print('Check mean and std after preprocessing. Means should be close to zero. Stds should be close to one.')\n",
    "    mean_control = np.nanmean(data_norm, axis=(0, 1, 2))\n",
    "    std_control = np.nanstd(data_norm, axis=(0, 1, 2))\n",
    "    print('mean_control', mean_control)\n",
    "    print('std_control: ', std_control)\n",
    "\n",
    "    return np.asarray(mean, dtype=np.float32), np.asarray(std, dtype=np.float32)\n",
    "\n",
    "\n",
    "def normalize_images_per_channel(images, mean_train, std_train):\n",
    "    normalized_images = images - mean_train\n",
    "    normalized_images = normalized_images / std_train\n",
    "    print('normalized_images.dtype', normalized_images.dtype)\n",
    "\n",
    "    return normalized_images\n",
    "\n",
    "\n",
    "def denormalize_images_per_channel(images, mean_train, std_train):\n",
    "    denormalized_images = images * std_train\n",
    "    denormalized_images = denormalized_images + mean_train\n",
    "    print('denormalized_images.dtype', denormalized_images.dtype)\n",
    "\n",
    "    return denormalized_images\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "if run_preprocessing:\n",
    "    print('calculating mean and std per channel from the training data...')\n",
    "    MEAN_train, STD_train = get_mean_and_std_per_channel(data_train=X_train)\n",
    "    print(MEAN_train)\n",
    "    print(STD_train)\n",
    "    # Save mean and std to file\n",
    "    # mean_train_file = data_dir + '/UCmerced_mean_train.npy'\n",
    "    # std_train_file = data_dir + '/UCmerced_std_train.npy'\n",
    "    # np.save(mean_train_file, MEAN_train)\n",
    "    # np.save(std_train_file, STD_train)\n",
    "\n",
    "    # Load mean and std from file\n",
    "    # MEAN_train = np.load(mean_train_file)\n",
    "    # STD_train = np.load(std_train_file)\n",
    "\n",
    "    X_train_prepro = normalize_images_per_channel(images=X_train, mean_train=MEAN_train, std_train=STD_train)\n",
    "    X_test_prepro = normalize_images_per_channel(images=X_test, mean_train=MEAN_train, std_train=STD_train)\n",
    "    X_val_prepro = normalize_images_per_channel(images=X_val, mean_train=MEAN_train, std_train=STD_train)\n",
    "\n",
    "    # Save preprocessed data split\n",
    "    # np.savez('pretrained_model_simpleCNN/intermediate_results/dataset_TrainValTest_float32.npz',\n",
    "    #        X_train_prepro=X_train_prepro, y_train=y_train,\n",
    "    #        X_val_prepro=X_val_prepro, y_val=y_val,\n",
    "    #        X_test_prepro=X_test_prepro, y_test=y_test,\n",
    "    #        MEAN_train=MEAN_train, STD_train=STD_train)\n",
    " \n",
    "    \n",
    "print('done')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preprocessing function which is used by the data generator\n",
    "\n",
    "As the kreas data generator expects a preprocessing function with only 1 input argument, we have hardcoded the training mean and std of this particular dataset and split here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MEAN_train = np.array([123.92658,  125.36055,  115.218414], dtype=np.float32)\n",
    "STD_train = np.array([55.47791355, 51.37877802, 49.92980842], dtype=np.float32)\n",
    "\n",
    "def normalize_UCMerced(image):\n",
    "    np.array([image])\n",
    "    \n",
    "    mean_train = np.array([123.92658,  125.36055,  115.218414], dtype=np.float32)\n",
    "    std_train = np.array([55.47791355, 51.37877802, 49.92980842], dtype=np.float32)\n",
    "    \n",
    "    normalized_image = image - mean_train\n",
    "    normalized_image = normalized_image / std_train\n",
    "    return normalized_image"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train the model\n",
    "\n",
    "The model was trained on a GPU (Titan Xp 12GB) for 200 epochs within 20 minutes (6 seconds per epoch). This code loads the weights after the epoch 200. You will train one epoch more and plot the training/validation curves for the entire 201 epochs. \n",
    "\n",
    "If you have not implemented the model correctly the weights of the pre-trained model will not match the model architecture. \n",
    "\n",
    "* Does the model learn anything better than a random classifier?\n",
    "* Is the model overfitting?\n",
    "* After how many epochs should we stop training? Or should we train even longer?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set hyper parameters\n",
    "# Set nb_epoch to 201 if you want to train 1 more epoch on your machine. Set it to 200 for no further training \n",
    "nb_epoch = 201\n",
    "base_learning_rate=0.0001\n",
    "batch_size = 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the number of batches per epoch. \n",
    "# One epoch is defined as an update pass through the full training data. \n",
    "# One iteration is defined as an update through one mini-batch.\n",
    "batches_per_epoch = train_paths.shape[0]//batch_size\n",
    "# the number of batches to see the full validation data:\n",
    "validation_steps = val_paths.shape[0]//batch_size\n",
    "test_steps = test_paths.shape[0]\n",
    "print('number of images per batch: {}'.format(batch_size))\n",
    "print('batches per epoch: {}'.format(batches_per_epoch))\n",
    "print('validation steps: {}'.format(validation_steps))\n",
    "print('test steps: {}'.format(test_steps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Data augmentation\n",
    "\n",
    "Which data augmentation techniques are useful for our dataset?\n",
    "\n",
    "**Task**: Implement some data augmentation operations. Look at the implemented keras tools (https://keras.io/preprocessing/image/) and choose the one that make sense for this problem.\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data Generator provided by keras to load the training data in batches\n",
    "# Usage if all data is in memory: generator=image_gen.flow(X_train_prepro, Y_train, batch_size=batch_size),\n",
    "\n",
    "image_gen = ImageDataGenerator(\n",
    "    preprocessing_function=normalize_UCMerced\n",
    "    # TODO: Add some arguments for data augmentation ...\n",
    "\n",
    ")\n",
    "\n",
    "# At Test/Validation time without data augmentation\n",
    "image_gen_test = ImageDataGenerator(preprocessing_function=normalize_UCMerced)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create iterators that read the data in batches from the train, val, or test directory\n",
    "\n",
    "This needs much less memory than loading the entire training dataset into memory.\n",
    "Especially the images converted to float32 consume much more memory.\n",
    "\n",
    "However, data loading can be the bottleneck when training on a GPU. We would have to make sure that the data reading is fast enough and does not slow down the training process. This is usually regulated with a queue and multiple workers that read the images in parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_it = image_gen.flow_from_directory('data/UCMerced_LandUse/train/', \n",
    "                                       classes=category_names, class_mode='categorical', \n",
    "                                       batch_size=batch_size, target_size=(256, 256))\n",
    "\n",
    "val_it = image_gen_test.flow_from_directory('data/UCMerced_LandUse/val/', \n",
    "                                       classes=category_names, class_mode='categorical', \n",
    "                                       batch_size=batch_size, target_size=(256, 256))\n",
    "\n",
    "# Let's look at the first return from the training iterator\n",
    "batchX, batchy = train_it.next()\n",
    "print('Batch shape=%s, min=%.3f, max=%.3f' % (batchX.shape, batchX.min(), batchX.max()))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimize\n",
    "\n",
    "**Attention**: If your model architecture is not as expected, there will be an error when loading the weights, because the number of parameters is not matching."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "load_pre_trained_model = True  # load the pretrained model weights\n",
    "continue_training = True if nb_epoch > 200 else False  # to load the history (loss curves) for the 200 pretrained epochs\n",
    "\n",
    "# set your architecture\n",
    "model = model_simple_CNN\n",
    "\n",
    "# Define an optimizer\n",
    "opt = Adam(lr=base_learning_rate, beta_1=0.9, beta_2=0.999, epsilon=None, decay=0.0, amsgrad=False)\n",
    "\n",
    "if load_pre_trained_model:\n",
    "    print('Loading weights from pre-trained model...')\n",
    "    weights_path='pretrained_model_simpleCNN_v2/weights_simpleCNN_epoch200.h5'\n",
    "    model.load_weights(filepath=weights_path, by_name=False, skip_mismatch=False)\n",
    "    \n",
    "    # Note: We could also load the full model (architecture + weights)\n",
    "    #  model = load_model('pretrained_model_simpleCNN_v2/model_simpleCNN_epoch3.h5')\n",
    "    initial_epoch = 200\n",
    "else:\n",
    "    # re-initialize weights of the model to train from scratch\n",
    "    print('initializing weights...')\n",
    "    initialize_weights(model)\n",
    "    initial_epoch = 0\n",
    "\n",
    "model.compile(optimizer=opt, loss='categorical_crossentropy', metrics=['accuracy'])\n",
    "\n",
    "# Save model checkpoints/parameters if the performance has improved.\n",
    "checkpoint_best_path = 'weights_best.hdf5'\n",
    "checkpoint_best = ModelCheckpoint(checkpoint_best_path, monitor='val_loss', verbose=1, save_best_only=True, mode='min', save_weights_only=True, period=1)\n",
    "        \n",
    "# Fit the model parameters\n",
    "print('training...')\n",
    "history = model.fit_generator(\n",
    "                generator=train_it,\n",
    "                steps_per_epoch=batches_per_epoch, epochs=nb_epoch, \n",
    "                validation_data=val_it,\n",
    "                validation_steps=validation_steps,\n",
    "                callbacks=[checkpoint_best],\n",
    "                initial_epoch=initial_epoch,\n",
    "                workers = 2,\n",
    "                max_q_size=10)\n",
    "\n",
    "print('done')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the training history"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open('pretrained_model_simpleCNN_v2/trainHistoryDict_epoch200.pkl', 'rb') as f:\n",
    "    history_pretrained = pickle.load(f, encoding='latin1')\n",
    "\n",
    "        \n",
    "if continue_training:\n",
    "    # Load the previous history\n",
    "    # Merge previous history until epoch 200 with new history\n",
    "    for key in history_pretrained.keys():\n",
    "            history.history[key] = history_pretrained[key] + history.history[key]    \n",
    "else:\n",
    "    history.history = history_pretrained\n",
    "    # # Save the history\n",
    "    with open('trainHistoryDict_epoch{}.pkl'.format(nb_epoch), 'wb') as file:\n",
    "        pickle.dump(history.history, file)\n",
    "        \n",
    "# Save the model state of the last epoch    \n",
    "model.save_weights('weights_simpleCNN_epoch{}.h5'.format(nb_epoch))\n",
    "model.save('model_simpleCNN_epoch{}.h5'.format(nb_epoch))\n",
    "print('done')\n",
    "\n",
    "# Plot history\n",
    "plot_train_val_accuracy_and_loss(history.history, legend_suffix='original', figsize=(8,6))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Prediction on the test data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# create another iterator that load all test images and labels with one iteration (needed for later analysis)\n",
    "test_it_load = image_gen_test.flow_from_directory('data/UCMerced_LandUse/test/', \n",
    "                                       classes=category_names, class_mode='categorical', \n",
    "                                       batch_size=test_steps, target_size=(256, 256),\n",
    "                                       shuffle=False)\n",
    "X_test_prepro, y_test = test_it_load.next()\n",
    "print('test steps:', test_steps)\n",
    "print(X_test_prepro.shape)\n",
    "\n",
    "# convert the test labels from one hot encoded to integers\n",
    "y_test_num = np.argmax(np.squeeze(y_test), axis=1)\n",
    "y_test_cat = [category_names[l] for l in y_test_num]\n",
    "print('num. test labels: ', len(y_test_cat))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_prediction = True  # set this to false to load the saved predictions\n",
    "\n",
    "if run_prediction:\n",
    "    # load the weights form the best epoch\n",
    "    model.load_weights(filepath='pretrained_model_simpleCNN_v2/weights_best.hdf5', by_name=False, skip_mismatch=False)\n",
    "    print('predicting...')\n",
    "    predictions = model.predict(x=X_test_prepro, batch_size=batch_size, steps=None)\n",
    "    np.save('predictions.npy', predictions)\n",
    "else:\n",
    "    print('load predictions...')\n",
    "    predictions = np.load('pretrained_model_simpleCNN_v2/predictions.npy')\n",
    "\n",
    "# convert predictions    \n",
    "predictions_num = np.argmax(predictions, axis=1)\n",
    "predictions_cat = [category_names[pred] for pred in predictions_num]\n",
    "print('num. predictions: ', len(predictions_cat))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# To visualize some random samples shuffle the dataset \n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(X_test_prepro)\n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(predictions_cat)\n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(y_test_cat)\n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(predictions_num)\n",
    "np.random.seed(seed)\n",
    "np.random.shuffle(y_test_num)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluation on the test data\n",
    "\n",
    "1. Which categories are mainly confused?\n",
    "2. Which categories are easy for the classifier?\n",
    "3. Compare the results to the assessment you did at the beginning! How was your intuition?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def plot_confusion_matrix(cm, classes,\n",
    "                          normalize=False,\n",
    "                          title='Confusion matrix',\n",
    "                          cmap=plt.cm.gray,\n",
    "                          filepath=None,\n",
    "                          figsize=(12,12)):\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",
    "    fig = plt.figure(figsize=figsize)\n",
    "\n",
    "    plt.imshow(cm, interpolation='nearest', cmap=cmap, vmin=0, vmax=1.0)\n",
    "    plt.title(title)\n",
    "    plt.colorbar()\n",
    "    tick_marks = np.arange(len(classes))\n",
    "    plt.xticks(tick_marks, classes, rotation=90, fontsize=16)\n",
    "    plt.yticks(tick_marks, classes, fontsize=16)\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\", fontsize=9)\n",
    "\n",
    "    plt.ylabel('True label', fontsize=22)\n",
    "    plt.xlabel('Predicted label', fontsize=22)\n",
    "    plt.tight_layout()\n",
    "\n",
    "    if filepath is not None:\n",
    "        fig.savefig(filepath)\n",
    "                \n",
    "\n",
    "# classification report\n",
    "report_string = classification_report(y_test_num, predictions_num, target_names=category_names)\n",
    "print(report_string)\n",
    "                \n",
    "# Overall accuracy\n",
    "print('num. correct samples:  {} \\nnum. all test samples: {}'.format(np.sum(y_test_num == predictions_num) , predictions_num.shape[0]))\n",
    "overall_acc = np.sum(y_test_num == predictions_num) / float(predictions_num.shape[0])\n",
    "print('overall accuracy: {:.2f}'.format(overall_acc))\n",
    "\n",
    "# Confusion matrix\n",
    "confmat = confusion_matrix(y_test_num, predictions_num )\n",
    "plot_confusion_matrix(cm=confmat, classes=category_names, normalize=True, cmap=plt.cm.Blues, filepath=None)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization of some test examples\n",
    "\n",
    "**Task**: Let's look at some specific error cases (highlighted with a red frame). Can we explain why the network might fail in some of the cases?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# VISUALIZE IMAGES\n",
    "def plotImages_categories( images, labels, predictions, n_rows=5, n_cols=4, figsize=(10, 10)):\n",
    "    fig, axes = plt.subplots(n_rows, n_cols, figsize=figsize)\n",
    "    axes = axes.flatten()\n",
    "    for i in range(len(axes)):\n",
    "        \n",
    "        axes[i].imshow(images[i])\n",
    "        axes[i].set_xticks(())\n",
    "        axes[i].set_yticks(())\n",
    "        if predictions is not None:\n",
    "            title = 'GT:   {}\\nPRED: {}'.format(labels[i], predictions[i])\n",
    "            \n",
    "            if predictions[i] != labels[i]:\n",
    "                [j.set_linewidth(5) for j in axes[i].spines.values()]\n",
    "                [j.set_color('red') for j in axes[i].spines.values()]\n",
    "\n",
    "        else:\n",
    "            title = 'GT: {}'.format(labels[i])\n",
    "        axes[i].set_title(title, fontdict={'family':'monospace'}, loc='left')\n",
    "    plt.tight_layout()\n",
    "    plt.show()\n",
    "\n",
    "# prepare the images for visualization: 1. denormalize, 2. convert to unint8\n",
    "X_test_vis = np.asarray(denormalize_images_per_channel(images=X_test_prepro, mean_train=MEAN_train, std_train=STD_train),\n",
    "                        dtype=np.uint8)\n",
    "\n",
    "plotImages_categories(images=X_test_vis, labels=y_test_cat, predictions=predictions_cat, \n",
    "                      n_rows=20, n_cols=4, figsize=(12, 60))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bonus ... or homework\n",
    "\n",
    "1. What is the effect of preprocessing/normalizing the input data? Train and test without preprocessing, but still convert the images from uint8 to float32.\n",
    "2. Similar, test the effect of the data augmentation. Run a baseline without data augmentation.\n",
    "3. Try to improve the performance by modifying/extending the CNN architectrue.\n",
    "\n",
    "Our implemented CNN is a pretty simple, yet powerful enough to learn generalizing features that solve this scene classification task with a decent performance.\n",
    "\n",
    "Below, is the code for the famous ResNet50 architecture [1].\n",
    "It is a very deep network constructed with building blocks. Within each block a skip connection directly forwards the input to the output. These skip connections allow to train much deeper networks and help to propagate the gradient during the optimization procedure.\n",
    "\n",
    "As a homework you may want to look at the code and the model summary (i.e. look at the number of trainable parameters). \n",
    "Furthermore, you could try to train the ResNet50 on the UC Merced dataset from scratch or fine-tune the weights trained on ImageNet.\n",
    "\n",
    "**Note**: Once you start using CNNs for your own projects and datasets, it might be a good idea to take a state-of-the-art architecture such as the ResNet50 and apply it to your new problem. Usually we can also find some pre-trained weights that can be used to initialize the network. This may be crucial, especially when your reference dataset is too small to train such a deep network from scratch.\n",
    "\n",
    "After testing the state-of-the-art you may want to modify the existing architecture to better suite your specific task. Think about what is specific about your data and the task you want to solve.\n",
    "\n",
    "**Keras already provides many CNN architectures including pre-trained weights**\n",
    "\n",
    "For an easy start, check out: https://keras.io/applications/\n",
    "\n",
    "[1] He, Kaiming, et al. \"Deep residual learning for image recognition.\" Proceedings of the IEEE conference on computer vision and pattern recognition. 2016."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def identity_block(input_tensor, kernel_size, filters, stage, block):\n",
    "    \"\"\"\n",
    "    The identity_block is the block that has no conv layer at shortcut\n",
    "    Arguments\n",
    "        input_tensor: input tensor\n",
    "        kernel_size: defualt 3, the kernel size of middle conv layer at main path\n",
    "        filters: list of integers, the nb_filters of 3 conv layer at main path\n",
    "        stage: integer, current stage label, used for generating layer names\n",
    "        block: 'a','b'..., current block label, used for generating layer names\n",
    "    \"\"\"\n",
    "\n",
    "    nb_filter1, nb_filter2, nb_filter3 = filters\n",
    "    conv_name_base = 'res' + str(stage) + block + '_branch'\n",
    "    bn_name_base = 'bn' + str(stage) + block + '_branch'\n",
    "\n",
    "    x = Conv2D(filters=nb_filter1, kernel_size=(1, 1), name=conv_name_base + '2a')(input_tensor)\n",
    "\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2a')(x)\n",
    "    x = Activation('relu')(x)\n",
    "\n",
    "    x = Conv2D(filters=nb_filter2, kernel_size=(kernel_size, kernel_size), \n",
    "               padding='same', name=conv_name_base + '2b')(x)\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2b')(x)\n",
    "    x = Activation('relu')(x)\n",
    "\n",
    "    x = Conv2D(filters=nb_filter3, kernel_size=(1, 1), name=conv_name_base + '2c')(x)\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2c')(x)\n",
    "\n",
    "    x = add([x, input_tensor])\n",
    "    x = Activation('relu')(x)\n",
    "    return x\n",
    "\n",
    "\n",
    "def conv_block(input_tensor, kernel_size, filters, stage, block, strides=(2, 2)):\n",
    "    \"\"\"\n",
    "    conv_block is the block that has a conv layer at shortcut\n",
    "    # Arguments\n",
    "        input_tensor: input tensor\n",
    "        kernel_size: defualt 3, the kernel size of middle conv layer at main path\n",
    "        filters: list of integers, the nb_filters of 3 conv layer at main path\n",
    "        stage: integer, current stage label, used for generating layer names\n",
    "        block: 'a','b'..., current block label, used for generating layer names\n",
    "    Note that from stage 3, the first conv layer at main path is with subsample=(2,2)\n",
    "    And the shortcut should have subsample=(2,2) as well\n",
    "    \"\"\"\n",
    "\n",
    "    nb_filter1, nb_filter2, nb_filter3 = filters\n",
    "    conv_name_base = 'res' + str(stage) + block + '_branch'\n",
    "    bn_name_base = 'bn' + str(stage) + block + '_branch'\n",
    "\n",
    "    x = Conv2D(filters=nb_filter1, kernel_size=(1, 1), strides=strides,\n",
    "                      name=conv_name_base + '2a')(input_tensor)\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2a')(x)\n",
    "    x = Activation('relu')(x)\n",
    "\n",
    "    x = Conv2D(filters=nb_filter2, kernel_size= (kernel_size, kernel_size), padding='same',\n",
    "                      name=conv_name_base + '2b')(x)\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2b')(x)\n",
    "    x = Activation('relu')(x)\n",
    "\n",
    "    x = Conv2D(filters=nb_filter3, kernel_size=(1, 1), name=conv_name_base + '2c')(x)\n",
    "    x = BatchNormalization(axis=bn_axis, name=bn_name_base + '2c')(x)\n",
    "\n",
    "    shortcut = Conv2D(filters=nb_filter3, kernel_size=(1, 1), strides=strides,\n",
    "                             name=conv_name_base + '1')(input_tensor)\n",
    "    shortcut = BatchNormalization(axis=bn_axis, name=bn_name_base + '1')(shortcut)\n",
    "\n",
    "    x = add([x, shortcut])\n",
    "    x = Activation('relu')(x)\n",
    "    return x\n",
    "\n",
    "\n",
    "\n",
    "def ResNet50(img_input):\n",
    "    model = Sequential()\n",
    "\n",
    "    x = Conv2D(64, kernel_size=(7, 7), strides=(2, 2), name='conv1', padding='same')(img_input)\n",
    "    x = BatchNormalization(axis=bn_axis, name='bn_conv1')(x)\n",
    "    x = Activation('relu')(x)\n",
    "    x = MaxPooling2D((3, 3), strides=(2, 2))(x)\n",
    "\n",
    "\n",
    "    x = conv_block(x, 3, [64, 64, 256], stage=2, block='a', strides=(1, 1))\n",
    "    x = identity_block(x, 3, [64, 64, 256], stage=2, block='b')\n",
    "    x = identity_block(x, 3, [64, 64, 256], stage=2, block='c')\n",
    "    \n",
    "    x = conv_block(x, 3, [128, 128, 512], stage=3, block='a', strides=(1, 1))\n",
    "    x = identity_block(x, 3, [128, 128, 512], stage=3, block='b')\n",
    "    x = identity_block(x, 3, [128, 128, 512], stage=3, block='c')\n",
    "    x = identity_block(x, 3, [128, 128, 512], stage=3, block='d')\n",
    "    \n",
    "    x = conv_block(x, 3, [256, 256, 1024], stage=4, block='a', strides=(1, 1))\n",
    "    x = identity_block(x, 3, [256, 256, 1024], stage=4, block='b')\n",
    "    x = identity_block(x, 3, [256, 256, 1024], stage=4, block='c')\n",
    "    x = identity_block(x, 3, [256, 256, 1024], stage=4, block='d')\n",
    "    x = identity_block(x, 3, [256, 256, 1024], stage=4, block='e')\n",
    "    x = identity_block(x, 3, [256, 256, 1024], stage=4, block='f')\n",
    "\n",
    "    x = conv_block(x, 3, [512, 512, 2048], stage=5, block='a', strides=(1, 1))\n",
    "    x = identity_block(x, 3, [512, 512, 2048], stage=5, block='b')\n",
    "    x = identity_block(x, 3, [512, 512, 2048], stage=5, block='c')\n",
    "\n",
    "    # Fully Connected Softmax Layer\n",
    "\n",
    "    x_fc = GlobalAveragePooling2D()(x)\n",
    "    x_fc = Dense(nb_classes, activation='softmax', name='your_output')(x_fc)\n",
    "\n",
    "    model = Model(img_input, x_fc)\n",
    "    return model\n",
    "\n",
    "\n",
    "input_shape = (patch_size, patch_size, channels)\n",
    "img_input = Input(shape=input_shape)\n",
    "bn_axis = 3\n",
    "\n",
    "model_resnet50 = ResNet50(img_input)\n",
    "model_resnet50.summary()"
   ]
  },
  {
   "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}