JLOWE/JLowe_NN.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Facies classification using Machine Learning\n",
    "## Joshua Lowe\n",
    "### https://uk.linkedin.com/in/jlowegeo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains my submission to the SEG Machine Learning contest 2016/17.\n",
    "I have implemented code to train a Neural Network and predict facies in a well from a variety or wireline logs.\n",
    "\n",
    "I have used bits of code from the original tutorial by Brendon Hall and from PA_Team, where I have used the 'blind well test' implemented by using leaveonegroupout. \n",
    "\n",
    "Thanks for all the different teams submissions as I have been able to learn a lot of skills around implementing machine learning algorithms in Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "np.random.seed(1000)\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import time as tm\n",
    "import pandas as pd\n",
    "from scipy.signal import medfilt\n",
    "\n",
    "from keras.models import Sequential\n",
    "from keras.constraints import maxnorm\n",
    "from keras.layers import Dense, Dropout\n",
    "from keras.utils import np_utils\n",
    "\n",
    "from sklearn.metrics import f1_score, confusion_matrix\n",
    "from sklearn.model_selection import LeaveOneGroupOut\n",
    "from sklearn import preprocessing\n",
    "\n",
    "#Cross Val of final model\n",
    "from sklearn.model_selection import cross_val_score, StratifiedKFold\n",
    "from keras.wrappers.scikit_learn import KerasClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "training_data = pd.read_csv('../training_data.csv')\n",
    "blind_data = pd.read_csv('../nofacies_data.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def accuracy(conf):\n",
    "    total_correct = 0.\n",
    "    nb_classes = conf.shape[0]\n",
    "    for i in np.arange(0,nb_classes):\n",
    "        total_correct += conf[i][i]\n",
    "    acc = total_correct/sum(sum(conf))\n",
    "    return acc\n",
    "\n",
    "adjacent_facies = np.array([[1], [0, 2], [1], [4], [3, 5], [4, 6, 7], [5, 7], [5, 6, 8], [6, 7]])\n",
    "\n",
    "\n",
    "def accuracy_adjacent(conf, adjacent_facies):\n",
    "    nb_classes = conf.shape[0]\n",
    "    total_correct = 0.\n",
    "    for i in np.arange(0,nb_classes):\n",
    "        total_correct += conf[i][i]\n",
    "        for j in adjacent_facies[i]:\n",
    "            total_correct += conf[i][j]\n",
    "    return total_correct / sum(sum(conf))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# 1=sandstone  2=c_siltstone   3=f_siltstone \n",
    "# 4=marine_silt_shale 5=mudstone 6=wackestone 7=dolomite\n",
    "# 8=packstone 9=bafflestone\n",
    "facies_colors = ['#F4D03F', '#F5B041','#DC7633','#6E2C00',\n",
    "       '#1B4F72','#2E86C1', '#AED6F1', '#A569BD', '#196F3D']\n",
    "\n",
    "facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS',\n",
    "                 'WS', 'D','PS', 'BS']\n",
    "#facies_color_map is a dictionary that maps facies labels\n",
    "#to their respective colors\n",
    "facies_color_map = {}\n",
    "for ind, label in enumerate(facies_labels):\n",
    "    facies_color_map[label] = facies_colors[ind]\n",
    "\n",
    "def label_facies(row, labels):\n",
    "    return labels[ row['Facies'] -1]\n",
    "    \n",
    "training_data.loc[:,'FaciesLabels'] = training_data.apply(lambda row: label_facies(row, facies_labels), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sorting the data and dropping unwanted columns from the training and test data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Leave the depth in as a predictor - can the NN recognise depth trends? - Other teams gone much further and have taken into account a predictors relationship/change with depth."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "X = training_data.drop(['Formation', 'Well Name', 'Facies', 'FaciesLabels'], axis=1).values\n",
    "y = training_data['Facies'].values - 1\n",
    "X_blind = blind_data.drop(['Formation', 'Well Name'], axis=1).values\n",
    "wells = training_data[\"Well Name\"].values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Scaling predictors in the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "scaler = preprocessing.RobustScaler().fit(X)\n",
    "X_scaled = scaler.transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Defining the neural network model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def DNN():\n",
    "    # Model\n",
    "    model = Sequential()\n",
    "    model.add(Dense(205, input_dim=8, activation='relu',W_constraint=maxnorm(5)))\n",
    "    model.add(Dropout(0.1))\n",
    "    model.add(Dense(69, activation='relu',W_constraint=maxnorm(5)))\n",
    "    model.add(Dropout(0.1))\n",
    "    model.add(Dense(69, activation='relu'))\n",
    "    model.add(Dense(9, activation='softmax'))\n",
    "    # Compilation\n",
    "    model.compile(loss='categorical_crossentropy', optimizer='adam', metrics=['accuracy'])\n",
    "    return model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross Validation using a 'Blind Well Test'. Code adapted from PA_Team submission "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     CHURCHMAN BIBLE f1_weigthted:0.443 | acc:0.495 | acc_adj:0.752\n",
      "      CROSS H CATTLE f1_weigthted:0.270 | acc:0.347 | acc_adj:0.856\n",
      "            LUKE G U f1_weigthted:0.450 | acc:0.423 | acc_adj:0.889\n",
      "               NEWBY f1_weigthted:0.464 | acc:0.447 | acc_adj:0.834\n",
      "               NOLAN f1_weigthted:0.485 | acc:0.453 | acc_adj:0.867\n",
      "          Recruit F9 f1_weigthted:0.929 | acc:0.868 | acc_adj:1.000\n",
      "             SHANKLE f1_weigthted:0.423 | acc:0.461 | acc_adj:0.973\n",
      "           SHRIMPLIN f1_weigthted:0.574 | acc:0.616 | acc_adj:0.941\n",
      "Avg F1 50.4949479748 Avg Acc 51.37284348 Avg Adj 88.914034357\n",
      "Blind Well Test Run Time: 163.758193 seconds\n"
     ]
    }
   ],
   "source": [
    "logo = LeaveOneGroupOut()\n",
    "t0 = tm.time()\n",
    "f1s_ls = []\n",
    "acc_ls = []\n",
    "adj_ls = []\n",
    "\n",
    "for train, test in logo.split(X_scaled, y, groups=wells):\n",
    "    well_name = wells[test[0]]\n",
    "    X_tr = X_scaled[train]\n",
    "    X_te = X_scaled[test]\n",
    "   \n",
    "    #convert y array into categories matrix\n",
    "    classes = 9\n",
    "    y_tr = np_utils.to_categorical(y[train], classes)\n",
    "    \n",
    "    # Method initialization\n",
    "    NN = DNN()\n",
    "    \n",
    "    # Training\n",
    "    NN.fit(X_tr, y_tr, nb_epoch=15, batch_size=5, verbose=0) \n",
    "    \n",
    "    # Predict\n",
    "    y_hat = NN.predict_classes(X_te, verbose=0)\n",
    "    y_hat = medfilt(y_hat, kernel_size=7)\n",
    "    \n",
    "    try:\n",
    "        f1s = f1_score(y[test], y_hat, average=\"weighted\", labels=[0, 1, 2, 3, 4, 5, 6, 7, 8])\n",
    "    except:\n",
    "        f1s = 0\n",
    "\n",
    "    try:\n",
    "        conf = confusion_matrix(y[test], y_hat, labels=[0, 1, 2, 3, 4, 5, 6, 7, 8])\n",
    "        acc = accuracy(conf) # similar to f1 micro\n",
    "    except:\n",
    "        acc = 0\n",
    "\n",
    "    try:\n",
    "        acc_adj = accuracy_adjacent(conf, adjacent_facies)\n",
    "    except:\n",
    "        acc_adj = 0\n",
    "\n",
    "    f1s_ls += [f1s]\n",
    "    acc_ls += [acc]\n",
    "    adj_ls += [acc_adj]\n",
    "    print(\"{:>20s} f1_weigthted:{:.3f} | acc:{:.3f} | acc_adj:{:.3f}\".format(well_name, f1s, acc, acc_adj))\n",
    "\n",
    "t1 = tm.time()\n",
    "print(\"Avg F1\", np.average(f1s_ls)*100, \"Avg Acc\", np.average(acc_ls)*100, \"Avg Adj\", np.average(adj_ls)*100)\n",
    "print(\"Blind Well Test Run Time:\",'{:f}'.format((t1-t0)), \"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Cross Validation using stratified K-fold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.49613602  0.42967543  0.33436533  0.53715171  0.45510836]\n",
      "Cross Validation Run Time: 95.910309 seconds\n"
     ]
    }
   ],
   "source": [
    "#Another robustness test of the model using statified K fold\n",
    "X_train = X_scaled\n",
    "Y_train = np_utils.to_categorical(y, classes)\n",
    "t2 = tm.time()\n",
    "estimator = KerasClassifier(build_fn=DNN, nb_epoch=15, batch_size=5, verbose=0)\n",
    "skf = StratifiedKFold(n_splits=5, shuffle=True)\n",
    "results_dnn = cross_val_score(estimator, X_train, Y_train, cv= skf.get_n_splits(X_train, Y_train))\n",
    "print (results_dnn)\n",
    "t3 = tm.time()\n",
    "print(\"Cross Validation Run Time:\",'{:f}'.format((t3-t2)), \"seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Final Model which uses all the training data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By using all the training data I may be potentially increasing the variance of the model but I believe it’s best to use all the data in the model as the data available is limited."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "f1 training error:  0.750747\n",
      "f1 test error:  50.494948\n"
     ]
    }
   ],
   "source": [
    "NN = DNN()\n",
    "NN.fit(X_train, Y_train, nb_epoch=15, batch_size=5, verbose=0)\n",
    "\n",
    "y_predicted = NN.predict_classes(X_train, verbose=0)\n",
    "y_predicted = medfilt(y_predicted, kernel_size=7)\n",
    "\n",
    "f1s = f1_score(y, y_predicted, average=\"weighted\")\n",
    "Avgf1s = np.average(f1s_ls)*100\n",
    "print (\"f1 training error: \", '{:f}'.format(f1s))\n",
    "print (\"f1 test error: \", '{:f}'.format(Avgf1s))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "My variance is high and my bias is too low.\n",
    "\n",
    "I haven’t found the optimum bias-variance trade off. --> Back to the drawing board."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Predicting the lithologies in the unknown test wells "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_blind = scaler.transform(X_blind)\n",
    "y_blind = NN.predict_classes(x_blind, verbose=0)\n",
    "y_blind = medfilt(y_blind, kernel_size=7)\n",
    "blind_data[\"Facies\"] = y_blind + 1  # return the original value (1-9)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "blind_data.to_csv(\"J_Lowe_Submission.csv\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda root]",
   "language": "python",
   "name": "conda-root-py"
  },
  "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}