Seismic Classifiers-Pub.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we'll import all the libraries we need down the line. We also set the \"random seed\", so results can be reproduced by avid readers. Keras should report using the Tensorflow backend, otherwise reproducibility cannot be guaranteed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import obspy\n",
    "import keras\n",
    "import time\n",
    "from keras_tqdm import TQDMNotebookCallback\n",
    "from tqdm import tnrange, tqdm_notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "import tensorflow as tf\n",
    "from obspy.io.segy.segy import _read_segy\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "np.random.seed(42)\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For experimentation with network models, we keep the keras imports separate, to reduce loading time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.models import Sequential, Model, clone_model\n",
    "from keras.layers import Conv2D, Dense, Activation, Flatten, Dropout, Input\n",
    "from keras.layers.normalization import BatchNormalization\n",
    "from keras.preprocessing.image import ImageDataGenerator\n",
    "from keras.optimizers import SGD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to define some parameters. As we are using Transfer learning, we have to adjust these parameters to fit into the network that we use and test. \n",
    "\n",
    "| Model   |      Channels      |  Patch-Size |\n",
    "|----------|:-------------:|------:|\n",
    "| Waldeland |  1 | 64 |\n",
    "| VGG16 |    3   |  64 |\n",
    "| ResNet50 | 3 | 244 |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "patch_size = 64 # for ResNet50 put 244\n",
    "batch_size = 256\n",
    "num_channels = 1\n",
    "num_classes = 9\n",
    "all_examples = 158812\n",
    "num_examples = 7500\n",
    "epochs = 20\n",
    "steps=450\n",
    "sampler = list(range(all_examples))\n",
    "\n",
    "opt = 'adam'\n",
    "lossfkt = ['categorical_crossentropy']\n",
    "metrica = ['mae', 'acc']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we test, whether we are running on CPU or GPU. We want to run on GPU, if it's not in the device list. It will be slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from tensorflow.python.client import device_lib\n",
    "print(device_lib.list_local_devices())\n",
    "# It should say GPU here. Otherwise your model will run sloooow."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data Loading\n",
    "Now let's load the F3 data and read three slices. The labeled data, as well as, a distal inline and a crossline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "filename = 'data/Dutch Government_F3_entire_8bit seismic.segy'\n",
    "\n",
    "t0=time.time()\n",
    "stream0 = _read_segy(filename, headonly=True)\n",
    "print('--> data read in {:.1f} sec'.format(time.time()-t0)) #Thanks to aadm \n",
    "\n",
    "t0=time.time()\n",
    "\n",
    "labeled_data = np.stack(t.data for t in stream0.traces if t.header.for_3d_poststack_data_this_field_is_for_in_line_number == 339).T\n",
    "inline_data = np.stack(t.data for t in stream0.traces if t.header.for_3d_poststack_data_this_field_is_for_in_line_number == 500).T\n",
    "xline_data = np.stack(t.data for t in stream0.traces if t.header.for_3d_poststack_data_this_field_is_for_cross_line_number == 500).T\n",
    "\n",
    "print('--> created slices in {:.1f} sec'.format(time.time()-t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Helper Functions\n",
    "From these slices, we need to extract patches. While, we could do that before and save them as array or image data, using a generator that utilizes the CPU, while the GPU trains the network is a bit more storage- and memory-friendly. `patch_extractor2D()` automates the patch-extraction and pads sides, where necessary.\n",
    "\n",
    "Then we build `acc_assess()` to format our test accuracy assessment nicely, because we're lazy and retyping it for every model we build is a nuisance.\n",
    "\n",
    "All functions are accompanied with a little sanity check. While this is not automated testing (like TDD), it does help to make sure, our function works as intended."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def patch_extractor2D(img,mid_x,mid_y,patch_size,dimensions=1):\n",
    "    try:\n",
    "        x,y,c = img.shape\n",
    "    except ValueError:\n",
    "        x,y = img.shape\n",
    "        c=1\n",
    "    patch= np.pad(img, patch_size//2, 'constant', constant_values=0)[mid_y:mid_y+patch_size,mid_x:mid_x+patch_size] #because it's padded we don't subtract half patches all the tim\n",
    "    if c != dimensions:\n",
    "        tmp_patch = np.zeros((patch_size,patch_size,dimensions))\n",
    "        for uia in range(dimensions):\n",
    "            tmp_patch[:,:,uia] = patch\n",
    "        return tmp_patch\n",
    "    return patch\n",
    "image=np.random.rand(10,10)//.1\n",
    "print(image)\n",
    "\n",
    "patch_extractor2D(image,10,10,4,1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def acc_assess(data,loss=['categorical_crossentropy'],metrics=['acc']):\n",
    "    if not isinstance(loss, list):\n",
    "        try:\n",
    "            loss = [loss]\n",
    "        except:\n",
    "            raise(\"Loss must be list.\")\n",
    "    if not isinstance(metrics, list):\n",
    "        try:\n",
    "            metrics = [metrics]\n",
    "        except:\n",
    "            raise(\"Metrics must be list.\")\n",
    "    out='The test loss is {:.3f}\\n'.format(data[0])\n",
    "    for i, metric in enumerate(metrics):            \n",
    "        if metric in 'mae':\n",
    "            out += \"The total mean error on the test is {:.3f}\\n\".format(data[i+1])\n",
    "        if metric in 'accuracy':\n",
    "            out += \"The test accuracy is {:.1f}%\\n\".format(data[i+1]*100)\n",
    "    return out\n",
    "print(acc_assess([1,2,3],'bla',[\"acc\", \"mae\"]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exploratory Data Analysis\n",
    "We need to load and check our labels. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labels = pd.read_csv('data/classification.ixz', delimiter=\" \", names=[\"Inline\",\"Xline\",\"Time\",\"Class\"])\n",
    "labels.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "labels[\"Xline\"]-=300-1\n",
    "labels[\"Time\"] = labels[\"Time\"]//4\n",
    "labels.describe()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "labeled_data.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "vml = np.percentile(labeled_data, 99)\n",
    "img1 = plt.imshow(labeled_data, cmap=\"Greys\", vmin=-vml, vmax=vml, aspect='auto')\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('labeled_data.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "vmx = np.percentile(xline_data, 99)\n",
    "plt.imshow(xline_data, cmap=\"Greys\", vmin=-vmx, vmax=vmx, aspect='auto')\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('xline_data.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "vmy = np.percentile(inline_data, 99)\n",
    "plt.imshow(inline_data, cmap=\"Greys\", vmin=-vmy, vmax=vmy, aspect='auto')\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('inline_data.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(labeled_data, cmap=\"Greys\", vmin=-vml, vmax=vml, aspect='auto')\n",
    "img1 = plt.scatter(labels[\"Xline\"],labels[[\"Time\"]],c=labels[[\"Class\"]],cmap='Dark2',alpha=0.03)\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('label.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train the Network\n",
    "Now we perform a test-train split. Then we can validate the results of our experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data, test_data, train_samples, test_samples = train_test_split(\n",
    "    labels, sampler, random_state=42)\n",
    "print(train_data.shape,test_data.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the `keras` data generator that wraps the `patch_extractor2D()`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SeismicSequence(keras.utils.Sequence):\n",
    "    def __init__(self, img, x_set, t_set, y_set, patch_size, batch_size, dimensions):\n",
    "        self.slice = img\n",
    "        self.X,self.t = x_set,t_set\n",
    "        self.batch_size = batch_size\n",
    "        self.patch_size = patch_size\n",
    "        self.dimensions = dimensions\n",
    "        self.label = y_set\n",
    "    \n",
    "    def __len__(self):\n",
    "        return len(self.X) // self.batch_size\n",
    "    \n",
    "    def __getitem__(self,idx):\n",
    "        sampler = np.random.permutation(len(self.X))\n",
    "        samples = sampler[idx*self.batch_size:(idx+1)*self.batch_size]\n",
    "        labels = keras.utils.to_categorical(self.label[samples], num_classes=9)\n",
    "        if self.dimensions == 1:\n",
    "            return np.expand_dims(np.array([patch_extractor2D(self.slice,self.X[x],self.t[x],self.patch_size,self.dimensions) for x in samples]), axis=4), labels\n",
    "        else:\n",
    "            return np.array([patch_extractor2D(self.slice,self.X[x],self.t[x],self.patch_size,self.dimensions) for x in samples]), labels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define several callbacks for keras. The training should be stopped early, if the validation loss or the categorical cross entropy do not improve within the defined patience. Checkpoints are written to `tmp.h5` for every epoch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earlystop1 = keras.callbacks.EarlyStopping(monitor='val_loss',\n",
    "                              min_delta=0,\n",
    "                              patience=3,\n",
    "                              verbose=0, mode='auto')\n",
    "\n",
    "earlystop2 = keras.callbacks.EarlyStopping(monitor='val_acc',\n",
    "                              min_delta=0,\n",
    "                              patience=3,\n",
    "                              verbose=0, mode='auto')\n",
    "\n",
    "checkpoint = keras.callbacks.ModelCheckpoint('tmp.h5', \n",
    "                                     monitor='val_loss', \n",
    "                                     verbose=0, \n",
    "                                     save_best_only=False, \n",
    "                                     save_weights_only=False, \n",
    "                                     mode='auto', \n",
    "                                     period=1)\n",
    "\n",
    "callbacklist = [TQDMNotebookCallback(leave_inner=True, leave_outer=True), earlystop1, earlystop2, checkpoint]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Waldeland CNN\n",
    "The model introduced by Waldeland, reproduced from MalenoV. Compared to today's standards this is a relatively shallow CNN. We train the network from scratch."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf.logging.set_verbosity(tf.logging.ERROR)\n",
    "\n",
    "model_vanilla = Sequential()\n",
    "model_vanilla.add(Conv2D(50, (5, 5), padding='same', input_shape=(patch_size,patch_size,1), strides=(4, 4), data_format=\"channels_last\",name = 'conv_layer1'))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Conv2D(50, (3, 3), strides=(2, 2), padding = 'same',name = 'conv_layer2'))\n",
    "model_vanilla.add(Dropout(0.5))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Conv2D(50, (3, 3), strides=(2, 2), padding= 'same',name = 'conv_layer3'))\n",
    "model_vanilla.add(Dropout(0.4))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Conv2D(50, (3, 3), strides=(2, 2), padding= 'same',name = 'conv_layer4'))\n",
    "model_vanilla.add(Dropout(0.2))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Conv2D(50, (3, 3), strides=(2, 2), padding= 'same',name = 'conv_layer5'))\n",
    "model_vanilla.add(Flatten())\n",
    "model_vanilla.add(Dense(50,name = 'dense_layer1'))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Dense(10,name = 'attribute_layer'))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('relu'))\n",
    "model_vanilla.add(Dense(num_classes, name = 'pre-softmax_layer'))\n",
    "model_vanilla.add(BatchNormalization())\n",
    "model_vanilla.add(Activation('softmax'))\n",
    "\n",
    "model_vanilla.compile(loss=lossfkt,\n",
    "                  optimizer=opt,\n",
    "                  metrics=metrica)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "t0=time.time()\n",
    "\n",
    "hist_vanilla = model_vanilla.fit_generator(\n",
    "    SeismicSequence(\n",
    "        labeled_data,\n",
    "        train_data[\"Xline\"].values,\n",
    "        train_data[\"Time\"].values,\n",
    "        train_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        1),\n",
    "    steps_per_epoch=steps,\n",
    "    validation_data = SeismicSequence(\n",
    "        labeled_data,\n",
    "        test_data[\"Xline\"].values,\n",
    "        test_data[\"Time\"].values,\n",
    "        test_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        1),\n",
    "    validation_steps = len(test_samples)//batch_size,\n",
    "    epochs = epochs,\n",
    "    verbose = 0,\n",
    "    callbacks = callbacklist)\n",
    "\n",
    "print('--> Training for Waldeland CNN took {:.1f} sec'.format(time.time()-t0)) #Thanks to aadm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_vanilla.save(\"vanilla_model.h5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "vanillascore=model_vanilla.evaluate(np.expand_dims(np.array([patch_extractor2D(labeled_data,labels[\"Xline\"][x],labels[\"Time\"][x],64) for x in test_samples]), axis=4),keras.utils.to_categorical(labels[\"Class\"][test_samples], num_classes=9), verbose=0)\n",
    "print(acc_assess(vanillascore,lossfkt,metrica))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the metric on training as well as validation gives a good overview, if we are doing appropriate training or if we are overfitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(hist_vanilla.history.keys())\n",
    "plt.plot(hist_vanilla.history['acc'])\n",
    "plt.plot(hist_vanilla.history['val_acc'])\n",
    "plt.title('model accuracy')\n",
    "plt.ylabel('accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# summarize history for loss\n",
    "plt.plot(hist_vanilla.history['loss'])\n",
    "plt.plot(hist_vanilla.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = xline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "predx = np.full_like(xline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max),leave=False, desc='Time'):\n",
    "        predx[depth,space] = np.argmax(model_vanilla.predict(np.expand_dims(np.expand_dims(patch_extractor2D(xline_data,space,depth,patch_size), axis=0), axis=4)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('vanilla_predx.npy',predx,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(predx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(xline_data, cmap=\"Greys\", vmin=-vmx, vmax=vmx, aspect='auto')\n",
    "img1 = plt.imshow(predx, aspect='auto', cmap=\"Dark2\", alpha=0.5)\n",
    "plt.savefig('pred1_x.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = inline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "predi= np.full_like(inline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max),leave=False, desc='Time'):\n",
    "        predi[depth,space] = np.argmax(model_vanilla.predict(np.expand_dims(np.expand_dims(patch_extractor2D(inline_data,space,depth,patch_size), axis=0), axis=4)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('vanilla_predi.npy',predi,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "predi = np.load('vanilla_predi.npy')\n",
    "plt.imshow(predi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "np.save('vanilla_predi.npy',predi,allow_pickle=False)\n",
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(inline_data, cmap=\"Greys\", vmin=-vmy, vmax=vmy, aspect='auto')\n",
    "img1 = plt.imshow(predi, aspect='auto', cmap=\"Dark2\", alpha=0.5)\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('pred1_i.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## VGG16 Transfer Learning\n",
    "We import the VGG16 model trained on the ImageNet dataset. We freeze all layers and cut off the classification part. We can then retrain the classification neurons, to see if the filters generalize to seismic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.applications.vgg16 import VGG16\n",
    "from keras import backend as K\n",
    "K.set_image_dim_ordering('tf')   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_tensor = Input(shape=(patch_size,patch_size,3))\n",
    "base_model = keras.applications.vgg16.VGG16(include_top=False, weights='imagenet', input_tensor=input_tensor, input_shape=None)\n",
    "\n",
    "for layer in base_model.layers[:8]:\n",
    "    layer.trainable = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = base_model.output\n",
    "x = Flatten()(x)\n",
    "x = Dense(256,name = 'dense_layer1')(x)\n",
    "x = BatchNormalization()(x)\n",
    "x = Activation('relu')(x)\n",
    "x = Dropout(.5)(x)\n",
    "x = Dense(num_classes, name = 'pre-softmax_layer')(x)\n",
    "x = BatchNormalization()(x)\n",
    "x = Activation('softmax')(x)\n",
    "\n",
    "vgg = Model(input=base_model.input, output=x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sgd = SGD(lr=1e-4, decay=1e-6, momentum=0.9, nesterov=True)\n",
    "vgg.compile(loss=lossfkt,\n",
    "                  optimizer=sgd,\n",
    "                  metrics=metrica)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "vgg_hist = vgg.fit_generator(\n",
    "    SeismicSequence(\n",
    "        labeled_data,\n",
    "        train_data[\"Xline\"].values,\n",
    "        train_data[\"Time\"].values,\n",
    "        train_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        3),\n",
    "    steps_per_epoch=steps,\n",
    "    validation_data = SeismicSequence(\n",
    "        labeled_data,\n",
    "        test_data[\"Xline\"].values,\n",
    "        test_data[\"Time\"].values,\n",
    "        test_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        3),\n",
    "    validation_steps = len(test_data)//batch_size,\n",
    "    epochs = epochs,\n",
    "    verbose = 0,\n",
    "    callbacks = callbacklist)\n",
    "\n",
    "print('--> Training for VGG transfer took {:.1f} sec'.format(time.time()-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vgg.save('vgg_model.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "vggscore=vgg.evaluate(np.array([patch_extractor2D(labeled_data,labels[\"Xline\"][x],labels[\"Time\"][x],64,3) for x in test_samples]), keras.utils.to_categorical(labels[\"Class\"][test_samples], num_classes=9))\n",
    "print(acc_assess(vggscore,lossfkt,metrica))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(hist_vanilla.history.keys())\n",
    "plt.plot(vgg_hist.history['acc'])\n",
    "plt.plot(vgg_hist.history['val_acc'])\n",
    "plt.title('model accuracy')\n",
    "plt.ylabel('accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize history for loss\n",
    "plt.plot(vgg_hist.history['loss'])\n",
    "plt.plot(vgg_hist.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = xline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "vgg_predx = np.full_like(xline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max),leave=False, desc='Time'):\n",
    "        vgg_predx[depth,space] = np.argmax(vgg.predict(np.expand_dims(patch_extractor2D(xline_data,space,depth,patch_size,3), axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('vgg_predx.npy',vgg_predx,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(vgg_predx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vgg_predx=np.load('vgg_predx.npy')\n",
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(xline_data, cmap=\"Greys\", vmin=-vmx, vmax=vmx, aspect='auto')\n",
    "img1 = plt.imshow(vgg_predx, aspect='auto', cmap=\"Dark2\", alpha=0.5)\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('vgg1_x.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = inline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "vgg_predi = np.full_like(inline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max),leave=False, desc='Time'):\n",
    "        vgg_predi[depth,space] = np.argmax(vgg.predict(np.expand_dims(patch_extractor2D(inline_data,space,depth,patch_size,3), axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('vgg_predi.npy',vgg_predi,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(vgg_predi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vgg_predi= np.load('vgg_predi.npy')\n",
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(inline_data, cmap=\"Greys\", vmin=-vmy, vmax=vmy, aspect='auto')\n",
    "img1 = plt.imshow(vgg_predi, aspect='auto', cmap=\"Dark2\", alpha=0.5)\n",
    "plt.yticks(np.arange(0, 462, 100), np.arange(0, 462*4, 400))\n",
    "plt.xlabel('Trace Location')\n",
    "plt.ylabel('Time [ms]')\n",
    "plt.savefig('vgg1_i.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ResNet50 Transfer Learning\n",
    "We import the ResNet50 that was trained on the ImageNet data and freeze all layers, like we did for the VGG16. Then we retrain the classifier to see if the learned filters generalize on seismic data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from keras.applications.resnet50 import ResNet50\n",
    "from keras import backend as K\n",
    "K.set_image_dim_ordering('tf')   \n",
    "patch_size=244"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_tensor = Input(shape=(patch_size,patch_size,3))\n",
    "res_base = ResNet50(include_top=False, weights='imagenet', input_tensor=input_tensor, input_shape=None, pooling=None)\n",
    "\n",
    "for layer in res_base.layers[:45]:\n",
    "    layer.trainable = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = res_base.output\n",
    "q = Flatten()(q)\n",
    "q = BatchNormalization()(q)\n",
    "q = Activation('relu')(q)\n",
    "q = Dense(10,name = 'attribute_layer')(q)\n",
    "q = BatchNormalization()(q)\n",
    "q = Activation('relu')(q)\n",
    "q = Dense(num_classes, name = 'pre-softmax_layer')(q)\n",
    "q = BatchNormalization()(q)\n",
    "q = Activation('softmax')(q)\n",
    "resnet = Model(input=res_base.input, output=q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sgd = SGD(lr=1e-3, decay=1e-6, momentum=0.9, nesterov=True)\n",
    "resnet.compile(loss=lossfkt,\n",
    "                  optimizer=sgd,\n",
    "                  metrics=metrica)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = time.time()\n",
    "\n",
    "batch_size=50\n",
    "res_hist = resnet.fit_generator(\n",
    "    SeismicSequence(\n",
    "        labeled_data,\n",
    "        train_data[\"Xline\"].values,\n",
    "        train_data[\"Time\"].values,\n",
    "        train_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        3),\n",
    "    steps_per_epoch=steps,\n",
    "    validation_data = SeismicSequence(\n",
    "        labeled_data,\n",
    "        test_data[\"Xline\"].values,\n",
    "        test_data[\"Time\"].values,\n",
    "        test_data[\"Class\"].values,\n",
    "        patch_size,\n",
    "        batch_size,\n",
    "        3),\n",
    "    validation_steps = len(test_data)//batch_size,\n",
    "    epochs = epochs,\n",
    "    verbose = 0,\n",
    "    callbacks = callbacklist)\n",
    "\n",
    "print('--> Training for ResNet transfer took {:.1f} sec'.format(time.time()-t0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resnet.save('resnet_model.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "resnetscore=resnet.evaluate(np.array([patch_extractor2D(labeled_data,labels[\"Xline\"][x],labels[\"Time\"][x],patch_size,3) for x in test_samples]), keras.utils.to_categorical(labels[\"Class\"][test_samples], num_classes=9))\n",
    "print(acc_assess(resnetscore,lossfkt,metrica))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = xline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "resnet_predx = np.full_like(xline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max),leave=False, desc='Time'):\n",
    "        resnet_predx[depth,space] = np.argmax(resnet.predict(np.expand_dims(patch_extractor2D(xline_data,space,depth,patch_size,3), axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('resnet_predx.npy',resnet_predx,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(resnet_predx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(xline_data, cmap=\"Greys\", vmin=-vmx, vmax=vmx, aspect='auto')\n",
    "img1 = plt.imshow(resnet_predx, aspect='auto', cmap=\"Dark2\", alpha=0.8)\n",
    "plt.savefig('resnet_x.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "t_max, y_max = inline_data.shape\n",
    "\n",
    "half_patch = patch_size//2\n",
    "\n",
    "resnet_predi = np.full_like(inline_data,-1)\n",
    "\n",
    "for space in tqdm_notebook(range(y_max-400,y_max-300),desc='Space'):\n",
    "    for depth in tqdm_notebook(range(t_max-400,t_max-300),leave=False, desc='Time'):\n",
    "        resnet_predi[depth,space] = np.argmax(resnet.predict(np.expand_dims(patch_extractor2D(inline_data,space,depth,patch_size,3), axis=0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.save('resnet_predi.npy',resnet_predi,allow_pickle=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(resnet_predi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig2 = plt.figure(figsize=(15.0, 10.0))\n",
    "img2 = plt.imshow(inline_data, cmap=\"Greys\", vmin=-vmy, vmax=vmy, aspect='auto')\n",
    "img1 = plt.imshow(resnet_predi, aspect='auto', cmap=\"Dark2\", alpha=0.8)\n",
    "plt.savefig('resnet_i.png', bbox_inches='tight')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(res_hist.history.keys())\n",
    "plt.plot(res_hist.history['acc'])\n",
    "plt.plot(res_hist.history['val_acc'])\n",
    "plt.title('model accuracy')\n",
    "plt.ylabel('accuracy')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# summarize history for loss\n",
    "plt.plot(res_hist.history['loss'])\n",
    "plt.plot(res_hist.history['val_loss'])\n",
    "plt.title('model loss')\n",
    "plt.ylabel('loss')\n",
    "plt.xlabel('epoch')\n",
    "plt.legend(['train', 'test'], loc='upper left')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plot_model(resnet, to_file='model_resnet.png')\n",
    "plot_model(resnet, to_file='model_resnet_shapes.png', show_shapes=True)\n",
    "SVG(model_to_dot(resnet).create(prog='dot', format='svg'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Summary\n",
    "We can see the summaries of the layers in the model definitions. Leveraging high-dimensional CNNs that are already trained can be very valuable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "model_vanilla.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vgg.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "resnet.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}