CEsprey - RandomForest/Facies_Tree_Ensemble_Classifier.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "from sklearn import preprocessing\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.ensemble import GradientBoostingClassifier\n",
    "from sklearn.ensemble import ExtraTreesClassifier\n",
    "from sklearn.ensemble import VotingClassifier\n",
    "from sklearn import cross_validation\n",
    "from sklearn.cross_validation import KFold\n",
    "from sklearn.cross_validation import train_test_split\n",
    "from sklearn import metrics\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from classification_utilities import display_cm, display_adj_cm\n",
    "import matplotlib.pyplot as plt\n",
    "import pywt\n",
    "from scipy.signal import upfirdn\n",
    "from sklearn.model_selection import LeaveOneGroupOut"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "filename = 'training_data.csv'\n",
    "training_data = pd.read_csv(filename)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Create a difference vector for each feature e.g. x1-x2, x1-x3... x2-x3...\n",
    "\n",
    "feature_vectors = training_data.drop(['Formation','Facies'], axis=1)\n",
    "feature_vectors = feature_vectors[['Depth','GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE']]\n",
    "\n",
    "def difference_vector(feature_vectors):\n",
    "    length = len(feature_vectors['Depth'])\n",
    "    df_temp = np.zeros((25, length))\n",
    "                          \n",
    "    for i in range(0,int(len(feature_vectors['Depth']))):\n",
    "                       \n",
    "        vector_i = feature_vectors.iloc[i,:]\n",
    "        vector_i = vector_i[['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE']]\n",
    "        for j, value_j in enumerate(vector_i):\n",
    "            for k, value_k in enumerate(vector_i): \n",
    "                differ_j_k = value_j - value_k          \n",
    "                df_temp[5*j+k, i] = np.abs(differ_j_k)\n",
    "                \n",
    "    return df_temp\n",
    "\n",
    "def diff_vec2frame(feature_vectors, df_temp):\n",
    "    \n",
    "    heads = feature_vectors.columns[1::]\n",
    "    for i in range(0,5):\n",
    "        string_i = heads[i]\n",
    "        for j in range(0,5):\n",
    "            string_j = heads[j]\n",
    "            col_head = 'diff'+string_i+string_j\n",
    "            \n",
    "            df = pd.Series(df_temp[5*i+j, :])\n",
    "            feature_vectors[col_head] = df\n",
    "    return feature_vectors\n",
    "            \n",
    "df_diff = difference_vector(feature_vectors)    \n",
    "feature_vectors = diff_vec2frame(feature_vectors, df_diff)\n",
    "\n",
    "# Drop duplicated columns and column of zeros\n",
    "feature_vectors = feature_vectors.T.drop_duplicates().T   \n",
    "feature_vectors.drop('diffGRGR', axis = 1, inplace = True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## sort by Formation and put the rank of each formation into a separate column. \n",
    "## Set rank of all rows not in specific formation to zero\n",
    "\n",
    "feature_vectors['Formation'] = training_data['Formation']\n",
    "# feature_vectors['Facies'] = training_data['Facies']\n",
    "\n",
    "def rank_formation(feature_vectors):\n",
    "    formation_labels = np.sort(feature_vectors['Formation'].unique())\n",
    "    feature_vec_origin = feature_vectors[['Formation','GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE']]\n",
    "    \n",
    "    for i, value_i in enumerate(formation_labels):\n",
    "        \n",
    "        formation = feature_vec_origin[feature_vec_origin['Formation']==value_i]\n",
    "        form_rank = formation.rank(pct=True, numeric_only=True)\n",
    "        \n",
    "        for j, value_j in enumerate(form_rank.columns):\n",
    "            \n",
    "            rank_form_vector = form_rank[value_j]\n",
    "            feature_vectors[value_j + '_' + value_i +'_rank'] = rank_form_vector\n",
    "            feature_vectors[value_j + '_' + value_i +'_rank'] = feature_vectors[value_j + '_' + value_i +'_rank'].fillna(0)\n",
    "            \n",
    "    return feature_vectors\n",
    "\n",
    "feature_vectors = rank_formation(feature_vectors)      "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## create DWT coeffs for scaled original feature vectors for first 8 scales. Upsample and hold the coeffs.\n",
    "import pywt\n",
    "from scipy.signal import upfirdn\n",
    "\n",
    "def DWT_facies(feature_vectors):\n",
    "    \n",
    "    feature_vec_origin = feature_vectors[['GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE']]\n",
    "    cols = feature_vec_origin.columns\n",
    "    \n",
    "    scaler = preprocessing.StandardScaler().fit(feature_vec_origin)\n",
    "    scaled_features = scaler.transform(feature_vec_origin)\n",
    "    scaled_features = pd.DataFrame(scaled_features, columns = cols)\n",
    "    \n",
    "    for i, value_i in enumerate(cols):\n",
    "        \n",
    "        sig = scaled_features[value_i]\n",
    "        \n",
    "        for j in np.arange(1,8):\n",
    "\n",
    "            fir_coeffs = np.ones((2**j,), dtype=np.int) \n",
    "            dwt_coeffs = pywt.downcoef('a', sig, 'db1', mode='symmetric', level=j)    \n",
    "            dwt_coeff_up = upfirdn(fir_coeffs, dwt_coeffs, 2**j)\n",
    "                 \n",
    "            feature_vectors[value_i + '_DWT_level' + str(j)] = dwt_coeff_up[0:len(sig)]\n",
    "\n",
    "    return feature_vectors\n",
    "\n",
    "feature_vectors = DWT_facies(feature_vectors)\n",
    "        \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Scale feature vectors\n",
    "\n",
    "feature_vectors_d = feature_vectors.drop(['Formation'], axis=1)\n",
    "scaler = preprocessing.StandardScaler().fit(feature_vectors_d)\n",
    "scaled_features = scaler.transform(feature_vectors_d)\n",
    "scaled_features = pd.DataFrame(scaled_features, columns = feature_vectors_d.columns)\n",
    "\n",
    "# One hot encode categoricals\n",
    "\n",
    "scaled_features['NM_M'] = training_data['NM_M']\n",
    "scaled_features['Formation'] = training_data['Formation']\n",
    "scaled_features['NM_M'] = scaled_features['NM_M'].astype('category')\n",
    "scaled_features = pd.get_dummies(scaled_features)\n",
    "\n",
    "#NM_M and Formation - remove from feature vector, scale feature vector and add NM_M and formation back before one hot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "predictors = scaled_features.columns\n",
    "correct_facies_labels = training_data['Facies'].values\n",
    "# Train-test split\n",
    "X_train, X_test, y_train, y_test = train_test_split(scaled_features,  correct_facies_labels, test_size=0.2, random_state=0)\n",
    "\n",
    "alg1 = RandomForestClassifier(random_state=1, n_estimators=315, min_samples_split=2, min_samples_leaf=1, max_features= None)\n",
    "alg2 = GradientBoostingClassifier(random_state=1, n_estimators=500, max_depth=6)\n",
    "alg3 = ExtraTreesClassifier(n_estimators=450, max_depth=None, min_samples_split=3, min_samples_leaf=1, random_state=1)\n",
    "\n",
    "eclf = VotingClassifier(estimators=[('rf', alg1), ('gbc', alg2), ('etc', alg3)], voting='hard')\n",
    "eclf.fit(X_train, y_train)\n",
    "\n",
    "predicted_vote = eclf.predict(X_test)\n",
    "\n",
    "print(\"prediction from tree ensemble:\")\n",
    "print(metrics.accuracy_score(list(y_test), predicted_vote))\n",
    "\n",
    "print(\"f1 score:\")\n",
    "print(metrics.f1_score(list(y_test), predicted_vote, average = 'weighted'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# alg1 = RandomForestClassifier(random_state=1, n_estimators=315, min_samples_split=2, min_samples_leaf=1, max_features= None)\n",
    "# alg1.fit(X_train, y_train)\n",
    "# predicted_vote = eclf.predict(X_test)\n",
    "# print(\"prediction from tree ensemble:\")\n",
    "# print(metrics.accuracy_score(list(y_test), predicted_vote))\n",
    "\n",
    "# print(\"f1 score:\")\n",
    "# print(metrics.f1_score(list(y_test), predicted_vote, average = 'weighted'))\n",
    "\n",
    "# print(\"Features sorted by their score:\")\n",
    "# print(sorted(zip(map(lambda x: round(x, 4), alg1.feature_importances_), predictors), \n",
    "#              reverse=True))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "f1_ens = []\n",
    "\n",
    "wells = training_data[\"Well Name\"].values\n",
    "logo = LeaveOneGroupOut()\n",
    "\n",
    "X = np.array(scaled_features)\n",
    "y = np.array(correct_facies_labels)\n",
    "for train, test in logo.split(X, y, groups=wells):\n",
    "    well_name = wells[test[0]]\n",
    "    eclf.fit(X[train], y[train])\n",
    "    pred_cv = eclf.predict(X[test])\n",
    "    sc = metrics.f1_score(y[test], pred_cv, labels = np.arange(10), average = 'micro')\n",
    "    print(\"{:>20s}  {:.3f}\".format(well_name, sc))\n",
    "    f1_ens.append(sc)\n",
    "    \n",
    "print(\"-Average leave-one-well-out F1 Score: %6f\" % (sum(f1_ens)/(1.0*(len(f1_ens)))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def facies_confusion(result, y_test):\n",
    "    facies_labels = ['SS', 'CSiS', 'FSiS', 'SiSh', 'MS',\n",
    "                     'WS', 'D','PS', 'BS']\n",
    "    conf = confusion_matrix(y_test, result)\n",
    "    display_cm(conf, facies_labels, hide_zeros=True, display_metrics = True)\n",
    "\n",
    "    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",
    "    print(accuracy(conf))\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",
    "    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))\n",
    "\n",
    "    print(accuracy_adjacent(conf, adjacent_facies))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "## Run on test data\n",
    "# read in Test data\n",
    "\n",
    "filename = 'validation_data_nofacies.csv'\n",
    "test_data = pd.read_csv(filename)\n",
    "\n",
    "# apply feature engineering to test_data\n",
    "feature_vectors_test = test_data.drop(['Formation'], axis=1)\n",
    "feature_vectors_test = feature_vectors_test[['Depth','GR', 'ILD_log10', 'DeltaPHI', 'PHIND', 'PE']]\n",
    "\n",
    "df_diff_test = difference_vector(feature_vectors_test)    \n",
    "feature_vectors_test = diff_vec2frame(feature_vectors_test, df_diff_test)\n",
    "\n",
    "# Drop duplicated columns and column of zeros\n",
    "feature_vectors_test = feature_vectors_test.T.drop_duplicates().T   \n",
    "feature_vectors_test.drop('diffGRGR', axis = 1, inplace = True)\n",
    "\n",
    "feature_vectors_test['Formation'] = test_data['Formation']\n",
    "feature_vectors_test = rank_formation(feature_vectors_test)  \n",
    "feature_vectors_test = DWT_facies(feature_vectors_test)\n",
    "\n",
    "# Scale feature vectors\n",
    "\n",
    "feature_vectors_test_d = feature_vectors_test.drop('Formation', axis=1)\n",
    "scaler = preprocessing.StandardScaler().fit(feature_vectors_test_d)\n",
    "scaled_features_test = scaler.transform(feature_vectors_test_d)\n",
    "scaled_features_test = pd.DataFrame(scaled_features_test, columns = feature_vectors_test_d.columns)\n",
    "\n",
    "# One hot encode categoricals\n",
    "\n",
    "scaled_features_test['NM_M'] = test_data['NM_M']\n",
    "scaled_features_test['Formation'] = test_data['Formation']\n",
    "scaled_features_test['NM_M'] = scaled_features_test['NM_M'].astype('category')\n",
    "scaled_features_test = pd.get_dummies(scaled_features_test)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "eclf.fit(scaled_features, correct_facies_labels)\n",
    "predicted_random_forest = eclf.predict(scaled_features_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "test_data['Facies'] = predicted_random_forest\n",
    "test_data.to_csv('test_data_prediction_ver2_CE.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}