LA_Team/Facies_classification_LA_TEAM_05.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Facies classification using Machine Learning #\n",
    "## LA Team Submission 5 ## \n",
    "### _[Lukas Mosser](https://at.linkedin.com/in/lukas-mosser-9948b32b/en), [Alfredo De la Fuente](https://pe.linkedin.com/in/alfredodelafuenteb)_ ####"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this approach for solving the facies classfication problem ( https://github.com/seg/2016-ml-contest. ) we will explore the following statregies:\n",
    "- Features Exploration: based on [Paolo Bestagini's work](https://github.com/seg/2016-ml-contest/blob/master/ispl/facies_classification_try02.ipynb), we will consider imputation, normalization and augmentation routines for the initial features.\n",
    "- Model tuning: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Libraries\n",
    "\n",
    "We will need to install the following libraries and packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied (use --upgrade to upgrade): pandas in /home/alfredo/anaconda2/lib/python2.7/site-packages\n",
      "Requirement already satisfied (use --upgrade to upgrade): python-dateutil in /home/alfredo/anaconda2/lib/python2.7/site-packages (from pandas)\n",
      "Requirement already satisfied (use --upgrade to upgrade): pytz>=2011k in /home/alfredo/anaconda2/lib/python2.7/site-packages (from pandas)\n",
      "Requirement already satisfied (use --upgrade to upgrade): numpy>=1.7.0 in /home/alfredo/anaconda2/lib/python2.7/site-packages (from pandas)\n",
      "Requirement already satisfied (use --upgrade to upgrade): six>=1.5 in /home/alfredo/anaconda2/lib/python2.7/site-packages (from python-dateutil->pandas)\n",
      "Requirement already satisfied (use --upgrade to upgrade): scikit-learn in /home/alfredo/anaconda2/lib/python2.7/site-packages\n",
      "Requirement already satisfied (use --upgrade to upgrade): tpot in /home/alfredo/anaconda2/lib/python2.7/site-packages\n",
      "Requirement already satisfied (use --upgrade to upgrade): deap in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): numpy in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): update-checker in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): scipy in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): tqdm in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): scikit-learn in /home/alfredo/anaconda2/lib/python2.7/site-packages (from tpot)\n",
      "Requirement already satisfied (use --upgrade to upgrade): requests>=2.3.0 in /home/alfredo/anaconda2/lib/python2.7/site-packages (from update-checker->tpot)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "You are using pip version 8.1.2, however version 9.0.1 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\n",
      "You are using pip version 8.1.2, however version 9.0.1 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\n",
      "You are using pip version 8.1.2, however version 9.0.1 is available.\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\n"
     ]
    }
   ],
   "source": [
    "%%sh\n",
    "pip install pandas\n",
    "pip install scikit-learn\n",
    "pip install tpot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "import numpy as np\n",
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import KFold , StratifiedKFold\n",
    "from classification_utilities import display_cm, display_adj_cm\n",
    "from sklearn.metrics import confusion_matrix, f1_score\n",
    "from sklearn import preprocessing\n",
    "from sklearn.model_selection import LeavePGroupsOut\n",
    "from sklearn.multiclass import OneVsOneClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from scipy.signal import medfilt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Preprocessing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#Load Data\n",
    "data = pd.read_csv('../facies_vectors.csv')\n",
    "\n",
    "# Parameters\n",
    "feature_names = ['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE', 'NM_M', 'RELPOS']\n",
    "facies_names = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS', 'WS', 'D', 'PS', 'BS']\n",
    "facies_colors = ['#F4D03F', '#F5B041','#DC7633','#6E2C00', '#1B4F72','#2E86C1', '#AED6F1', '#A569BD', '#196F3D']\n",
    "\n",
    "# Store features and labels\n",
    "X = data[feature_names].values \n",
    "y = data['Facies'].values \n",
    "\n",
    "# Store well labels and depths\n",
    "well = data['Well Name'].values\n",
    "depth = data['Depth'].values\n",
    "\n",
    "# Fill 'PE' missing values with mean\n",
    "imp = preprocessing.Imputer(missing_values='NaN', strategy='mean', axis=0)\n",
    "imp.fit(X)\n",
    "X = imp.transform(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We procceed to run [Paolo Bestagini's routine](https://github.com/seg/2016-ml-contest/blob/master/ispl/facies_classification_try02.ipynb) to include a small window of values to acount for the spatial component in the log analysis, as well as the gradient information with respect to depth. This will be our prepared training dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Feature windows concatenation function\n",
    "def augment_features_window(X, N_neig):\n",
    "    \n",
    "    # Parameters\n",
    "    N_row = X.shape[0]\n",
    "    N_feat = X.shape[1]\n",
    "\n",
    "    # Zero padding\n",
    "    X = np.vstack((np.zeros((N_neig, N_feat)), X, (np.zeros((N_neig, N_feat)))))\n",
    "\n",
    "    # Loop over windows\n",
    "    X_aug = np.zeros((N_row, N_feat*(2*N_neig+1)))\n",
    "    for r in np.arange(N_row)+N_neig:\n",
    "        this_row = []\n",
    "        for c in np.arange(-N_neig,N_neig+1):\n",
    "            this_row = np.hstack((this_row, X[r+c]))\n",
    "        X_aug[r-N_neig] = this_row\n",
    "\n",
    "    return X_aug\n",
    "\n",
    "\n",
    "# Feature gradient computation function\n",
    "def augment_features_gradient(X, depth):\n",
    "    \n",
    "    # Compute features gradient\n",
    "    d_diff = np.diff(depth).reshape((-1, 1))\n",
    "    d_diff[d_diff==0] = 0.001\n",
    "    X_diff = np.diff(X, axis=0)\n",
    "    X_grad = X_diff / d_diff\n",
    "        \n",
    "    # Compensate for last missing value\n",
    "    X_grad = np.concatenate((X_grad, np.zeros((1, X_grad.shape[1]))))\n",
    "    \n",
    "    return X_grad\n",
    "\n",
    "\n",
    "# Feature augmentation function\n",
    "def augment_features(X, well, depth, N_neig=1):\n",
    "    \n",
    "    # Augment features\n",
    "    X_aug = np.zeros((X.shape[0], X.shape[1]*(N_neig*2+2)))\n",
    "    for w in np.unique(well):\n",
    "        w_idx = np.where(well == w)[0]\n",
    "        X_aug_win = augment_features_window(X[w_idx, :], N_neig)\n",
    "        X_aug_grad = augment_features_gradient(X[w_idx, :], depth[w_idx])\n",
    "        X_aug[w_idx, :] = np.concatenate((X_aug_win, X_aug_grad), axis=1)\n",
    "    \n",
    "    # Find padded rows\n",
    "    padded_rows = np.unique(np.where(X_aug[:, 0:7] == np.zeros((1, 7)))[0])\n",
    "    \n",
    "    return X_aug, padded_rows\n",
    "\n",
    "X_aug, padded_rows = augment_features(X, well, depth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Initialize model selection methods\n",
    "lpgo = LeavePGroupsOut(2)\n",
    "\n",
    "# Generate splits\n",
    "split_list = []\n",
    "for train, val in lpgo.split(X, y, groups=data['Well Name']):\n",
    "    hist_tr = np.histogram(y[train], bins=np.arange(len(facies_names)+1)+.5)\n",
    "    hist_val = np.histogram(y[val], bins=np.arange(len(facies_names)+1)+.5)\n",
    "    if np.all(hist_tr[0] != 0) & np.all(hist_val[0] != 0):\n",
    "        split_list.append({'train':train, 'val':val})\n",
    "    \n",
    "        \n",
    "def preprocess():\n",
    "    \n",
    "    # Preprocess data to use in model\n",
    "    X_train_aux = []\n",
    "    X_test_aux = []\n",
    "    y_train_aux = []\n",
    "    y_test_aux = []\n",
    "    \n",
    "    # For each data split\n",
    "    split = split_list[5]\n",
    "        \n",
    "    # Remove padded rows\n",
    "    split_train_no_pad = np.setdiff1d(split['train'], padded_rows)\n",
    "\n",
    "    # Select training and validation data from current split\n",
    "    X_tr = X_aug[split_train_no_pad, :]\n",
    "    X_v = X_aug[split['val'], :]\n",
    "    y_tr = y[split_train_no_pad]\n",
    "    y_v = y[split['val']]\n",
    "\n",
    "    # Select well labels for validation data\n",
    "    well_v = well[split['val']]\n",
    "\n",
    "    # Feature normalization\n",
    "    scaler = preprocessing.RobustScaler(quantile_range=(25.0, 75.0)).fit(X_tr)\n",
    "    X_tr = scaler.transform(X_tr)\n",
    "    X_v = scaler.transform(X_v)\n",
    "        \n",
    "    X_train_aux.append( X_tr )\n",
    "    X_test_aux.append( X_v )\n",
    "    y_train_aux.append( y_tr )\n",
    "    y_test_aux.append (  y_v )\n",
    "    \n",
    "    X_train = np.concatenate( X_train_aux )\n",
    "    X_test = np.concatenate ( X_test_aux )\n",
    "    y_train = np.concatenate ( y_train_aux )\n",
    "    y_test = np.concatenate ( y_test_aux )\n",
    "    \n",
    "    return X_train , X_test , y_train , y_test "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data Analysis\n",
    "\n",
    "In this section we will run a Cross Validation routine "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alfredo/anaconda2/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n",
      "Optimization Progress: 100%|██████████| 20/20 [04:24<00:00, 34.90s/pipeline]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation 1 - Current best internal CV score: 0.585100621865\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimization Progress:  95%|█████████▌| 38/40 [07:55<00:24, 12.43s/pipeline]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation 2 - Current best internal CV score: 0.585100621865\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimization Progress:  93%|█████████▎| 56/60 [12:25<00:45, 11.38s/pipeline]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation 3 - Current best internal CV score: 0.585100621865\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimization Progress: 100%|██████████| 80/80 [18:27<00:00, 18.42s/pipeline]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation 4 - Current best internal CV score: 0.58540811006\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimization Progress: 100%|██████████| 100/100 [24:51<00:00, 20.19s/pipeline]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generation 5 - Current best internal CV score: 0.588497131428\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Optimization Progress:  96%|█████████▌| 115/120 [28:34<01:18, 15.78s/pipeline]"
     ]
    }
   ],
   "source": [
    "from tpot import TPOTClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = preprocess()\n",
    "\n",
    "tpot = TPOTClassifier(generations=5, population_size=20, \n",
    "                      verbosity=2,max_eval_time_mins=20,\n",
    "                      max_time_mins=100,scoring='f1_micro',\n",
    "                      random_state = 17)\n",
    "tpot.fit(X_train, y_train)\n",
    "print(tpot.score(X_test, y_test))\n",
    "tpot.export('FinalPipeline.py')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alfredo/anaconda2/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.\n",
      "  \"This module will be removed in 0.20.\", DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "from sklearn.ensemble import  RandomForestClassifier, VotingClassifier\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.naive_bayes import BernoulliNB\n",
    "from sklearn.pipeline import make_pipeline, make_union\n",
    "from sklearn.preprocessing import FunctionTransformer\n",
    "import xgboost as xgb\n",
    "from xgboost.sklearn import  XGBClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Train and test a classifier\n",
    "def train_and_test(X_tr, y_tr, X_v, well_v):\n",
    "    \n",
    "    # Feature normalization\n",
    "    scaler = preprocessing.RobustScaler(quantile_range=(25.0, 75.0)).fit(X_tr)\n",
    "    X_tr = scaler.transform(X_tr)\n",
    "    X_v = scaler.transform(X_v)\n",
    "    \n",
    "    # Train classifier\n",
    "    #clf = make_pipeline(make_union(VotingClassifier([(\"est\", ExtraTreesClassifier(criterion=\"gini\", max_features=1.0, n_estimators=500))]), FunctionTransformer(lambda X: X)), XGBClassifier(learning_rate=0.73, max_depth=10, min_child_weight=10, n_estimators=500, subsample=0.27))\n",
    "    #clf =  make_pipeline( KNeighborsClassifier(n_neighbors=5, weights=\"distance\") ) \n",
    "    #clf = make_pipeline(MaxAbsScaler(),make_union(VotingClassifier([(\"est\", RandomForestClassifier(n_estimators=500))]), FunctionTransformer(lambda X: X)),ExtraTreesClassifier(criterion=\"entropy\", max_features=0.0001, n_estimators=500))\n",
    "    # * clf = make_pipeline( make_union(VotingClassifier([(\"est\", BernoulliNB(alpha=60.0, binarize=0.26, fit_prior=True))]), FunctionTransformer(lambda X: X)),RandomForestClassifier(n_estimators=500))\n",
    "    clf = make_pipeline ( XGBClassifier(learning_rate=0.12, max_depth=3, min_child_weight=10, n_estimators=150, seed = 17, colsample_bytree = 0.9) )\n",
    "    clf.fit(X_tr, y_tr)\n",
    "    \n",
    "    # Test classifier\n",
    "    y_v_hat = clf.predict(X_v)\n",
    "    \n",
    "    # Clean isolated facies for each well\n",
    "    for w in np.unique(well_v):\n",
    "        y_v_hat[well_v==w] = medfilt(y_v_hat[well_v==w], kernel_size=5)\n",
    "    \n",
    "    return y_v_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Prediction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Load testing data\n",
    "test_data = pd.read_csv('../validation_data_nofacies.csv')\n",
    "\n",
    "# Prepare training data\n",
    "X_tr = X\n",
    "y_tr = y\n",
    "\n",
    "# Augment features\n",
    "X_tr, padded_rows = augment_features(X_tr, well, depth)\n",
    "\n",
    "# Removed padded rows\n",
    "X_tr = np.delete(X_tr, padded_rows, axis=0)\n",
    "y_tr = np.delete(y_tr, padded_rows, axis=0) \n",
    "\n",
    "# Prepare test data\n",
    "well_ts = test_data['Well Name'].values\n",
    "depth_ts = test_data['Depth'].values\n",
    "X_ts = test_data[feature_names].values\n",
    "\n",
    "# Augment features\n",
    "X_ts, padded_rows = augment_features(X_ts, well_ts, depth_ts)\n",
    "\n",
    "# Predict test labels\n",
    "y_ts_hat = train_and_test(X_tr, y_tr, X_ts, well_ts)\n",
    "\n",
    "# Save predicted labels\n",
    "test_data['Facies'] = y_ts_hat\n",
    "test_data.to_csv('Prediction_XX_Final.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [Root]",
   "language": "python",
   "name": "Python [Root]"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}