machine_learning_approach.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Projet CentraleSupelec - Dreem \n",
    "# Machine Learning course 3A OBT\n",
    "## Alexis Tuil et Adil Bousfiha\n",
    "### Approche Machine Learning Classique\n",
    "### Preprocessing + Features Extraction + Random Forest"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load useful libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import math\n",
    "import numpy as np\n",
    "from scipy.stats import moment, skew, kurtosis\n",
    "from scipy.signal import butter, lfilter, freqz, freqs\n",
    "import h5py\n",
    "from pyeeg import dfa, hurst\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define utils fonctions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "CHUNK_SIZE = 1000\n",
    "\n",
    "# frequence pour la decomposition du signal eeg\n",
    "\"\"\"In order to extract Energy features from our signals, \n",
    "    we filter each eeg between the following frequencies \n",
    "    to extract the associated wave type : \n",
    "    alpha [8–13 Hz], beta [12–24 Hz], theta [4–8 Hz], delta [0.5–2 Hz],\n",
    "    \"\"\"\n",
    "frequencies = [ [8,13],\n",
    "                [12,24],\n",
    "                [4,8],\n",
    "                [0.5,2]]\n",
    "\n",
    "order = 4 # ordre du filtre\n",
    "fs = 50.0 # frequence du filtre"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chunks(l, n):\n",
    "    \"\"\"Yield successive n-sized chunks from l.\"\"\"\n",
    "    for i in range(0, len(l), n):\n",
    "        yield l[i:i + n]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_features_accel(df):\n",
    "    res = [[l.var() , skew(l), kurtosis(l), min_max_d(l)] for l in df]\n",
    "    \n",
    "    return res\n",
    "\n",
    "def get_stats(l, with_pyeeg=True):\n",
    "    \"\"\"Compute statistis for the array argument\n",
    "       return an array \n",
    "       l : numpy 1D array\n",
    "    \"\"\"\n",
    "    if with_pyeeg:\n",
    "        return [l.var() , skew(l), kurtosis(l), min_max_d(l), energy(l), dfa(l), hurst(l)]\n",
    "    else:\n",
    "        return [l.var() , skew(l), kurtosis(l), min_max_d(l), energy(l)]\n",
    "\n",
    "def min_max_d(l):\n",
    "    \"\"\"Compute the min-max distance of the array argument\n",
    "        return a scalar\n",
    "        l : numpy 1D array\n",
    "    \"\"\"\n",
    "    a2 = math.pow(max(l) - min(l), 2)\n",
    "    b2 = math.pow(np.argmax(l) - np.argmin(l), 2) \n",
    "    \n",
    "    return math.sqrt( a2 + b2 )\n",
    "\n",
    "def energy(l):\n",
    "    \"\"\" Compute the energy of the array argument\n",
    "        return a scalar\n",
    "        l : numpy 1D array\n",
    "    \"\"\"\n",
    "    return sum(l*l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### On définit des fonctions utiles qui filtrent le signal échantilloné à la fréquence fs au travers d'un filtre passe-bande d'ordre order, et de bande passante [low, high]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def butter_lowpass(low, high, fs=fs, order=order):\n",
    "    nyq = 0.5 * fs\n",
    "    normal_low = low / nyq\n",
    "    normal_high = high / nyq\n",
    "    #DEBUG #print('butter_lowpass: %s %s' % (normal_low, normal_high))\n",
    "    b, a = butter(order, (normal_low, normal_high), btype='band', analog=False)\n",
    "    return b, a\n",
    "\n",
    "def butter_lowpass_filter(data, low, high, fs=fs, order=order):\n",
    "   #DEBUG #print('butter_lowpass_filter: %s %s' % (low, high))\n",
    "    b, a = butter_lowpass(low, high)\n",
    "    y = lfilter(b, a, data)\n",
    "    return y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def wave_extraction(l, freq, order=order, fs=fs):\n",
    "    \"\"\" Extract the wave with frequence range freq from l\n",
    "        return an array\n",
    "        l : numpy 1D array\n",
    "        freq : array of size 2 \"\"\"\n",
    "    #DEBUG #print('wave_extraction: %s' % (freq))\n",
    "    y = butter_lowpass_filter(l, freq[0], freq[1])\n",
    "    return(y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Feature Extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = h5py.File('all/train.h5', 'r')\n",
    "names = ['eeg_1', 'eeg_2', 'eeg_3', 'eeg_4', 'eeg_5', 'eeg_6', 'eeg_7']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Features sur les raw EEG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Extraction des features pour les Raw eeg:\n",
    "- min-max distance\n",
    "- Moment d'ordre 2\n",
    "- Energy\n",
    "- Variance\n",
    "- Kurtosis\n",
    "- Skewness"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for eeg_name in names:\n",
    "    res = []\n",
    "    for chunk in tqdm(chunks(data[eeg_name], CHUNK_SIZE)):\n",
    "        for i in tqdm(range(len(chunk))):\n",
    "            temp = [stats for freq in frequencies for stats in get_stats(chunk[i,:], with_pyeeg=False)]\n",
    "            res.append(temp)\n",
    "            np.savez(\"features/backup_stats_on_raw_\" + eeg_name + \"_train\", res)\n",
    "\n",
    "    np.savez(\"features/stats_on_raw_\" + eeg_name + \"_train\", res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### On fait la moyenne sur variance, skew et kurtosis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# load all matrices\n",
    "mat1 = np.load(\"features/stats_on_raw_\" + eeg_1 + \"_train.npz\")['arr_0'][,:3]\n",
    "mat2 = np.load(\"features/stats_on_raw_\" + eeg_2 + \"_train.npz\")['arr_0'][,:3]\n",
    "mat3 = np.load(\"features/stats_on_raw_\" + eeg_3 + \"_train.npz\")['arr_0'][,:3]  \n",
    "mat4 = np.load(\"features/stats_on_raw_\" + eeg_4 + \"_train.npz\")['arr_0'][,:3]  \n",
    "mat5 = np.load(\"features/stats_on_raw_\" + eeg_5 + \"_train.npz\")['arr_0'][,:3] \n",
    "mat6 = np.load(\"features/stats_on_raw_\" + eeg_6 + \"_train.npz\")['arr_0'][,:3]  \n",
    "mat7 = np.load(\"features/stats_on_raw_\" + eeg_7 + \"_train.npz\")['arr_0'][,:3]\n",
    "\n",
    "X = np.hstack((mat1, mat2, mat3, mat4, mat5, mat6, mat7))\n",
    "\n",
    "# compute the mean\n",
    "X_mean = np.vstack((np.mean(X[:,[0+0,0+3,0+6,0+9,0+12,0+15,0+18]], axis=1),\n",
    "                    np.mean(X[:,[1+0,1+3,1+6,1+9,1+12,1+15,1+18]], axis=1),\n",
    "                    np.mean(X[:,[2+0,2+3,2+6,2+9,2+12,2+15,2+18]], axis=1))).transpose()\n",
    "\n",
    "np.savez('features/stats_mean_train', X_mean)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Features sur les signaux EEG filtrés"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### On décompose chaque 30 sec epoch eeg en 4 sous signaux, puis on calcul les features statistiques pour chacun des sous signaux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for eeg_name in names:\n",
    "    res = []\n",
    "    for chunk in tqdm(chunks(data[eeg_name], CHUNK_SIZE)):\n",
    "        for i in tqdm(range(len(chunk))):\n",
    "            temp = [stats for freq in frequencies for stats in get_stats(wave_extraction(chunk[i,:], freq))]\n",
    "            res.append(temp)\n",
    "            np.savez(\"features/backup_stats_on_filter_\" + eeg_name + \"_train\", res)\n",
    "\n",
    "    np.savez(\"features/stats_on_filter_\" + eeg_name + \"_train\", res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Features sur les signaux de l'accelerometre 3D (x, y z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### On extrati 4 statistiques:\n",
    "- Variance\n",
    "- Skewness\n",
    "- Kurtosis\n",
    "- Min-max distance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = [i for chunk in chunks(data['accelerometer_x']) for i in compute_features_accel(chunk)]\n",
    "y = [i for chunk in chunks(data['accelerometer_y']) for i in compute_features_accel(chunk)]\n",
    "z = [i for chunk in chunks(data['accelerometer_z']) for i in compute_features_accel(chunk)]\n",
    "\n",
    "np.savez('features/stats_accel_train', np.hstack((x, y, z)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combine all features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mat2 = np.load('features/MMD_train.npz')['arr_0'] # (, 7)\n",
    "mat3 = np.load('features/M2_train.npz')['arr_0'] #(, 7)\n",
    "mat4 = np.load('features/stats_mean_train.npz')['arr_0'] #(, 3)\n",
    "mat5 = np.load('features/energy_train.npz')['arr_0'] #(, 42)\n",
    "mat6 = np.load('features/stats_accel_train.npz')['arr_0'] #(, 12)\n",
    "mat7 = np.load('features/stats_on_filter_eeg_1_train.npz')['arr_0'] #(, 28)\n",
    "mat8 = np.load('features/stats_on_filter_eeg_2_train.npz')['arr_0'] #(, 28)\n",
    "mat9 = np.load('features/stats_on_filter_eeg_3_train.npz')['arr_0'] #(, 28)\n",
    "mat10 = np.load('features/stats_on_filter_eeg_4_train.npz')['arr_0'] #(, 28)\n",
    "mat11 = np.load('features/stats_on_filter_eeg_5_train.npz')['arr_0'] #(, 28)\n",
    "mat12 = np.load('features/stats_on_filter_eeg_6_train.npz')['arr_0'] #(, 28)\n",
    "mat13 = np.load('features/stats_on_filter_eeg_7_train.npz')['arr_0'] #(, 28)\n",
    "X = np.hstack((mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, mat11, mat12, mat13))\n",
    "# X size : (, 267)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = pd.read_csv('all/train_y.csv').sleep_stage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=.2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Choice of the classifier\n",
    "voir en fin de notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "p = CompareClassifiers(X,y)\n",
    "p.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters Tuning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf = RandomForestClassifier(random_state = 42)\n",
    "# Number of trees in random forest\n",
    "n_estimators = [20]\n",
    "# Number of features to consider at every split\n",
    "max_features = ['auto', 'sqrt']\n",
    "# Maximum number of levels in tree\n",
    "max_depth = [int(x) for x in np.linspace(10, 110, num = 11)]\n",
    "max_depth.append(None)\n",
    "# Minimum number of samples required to split a node\n",
    "min_samples_split = [2, 5, 10]\n",
    "# Minimum number of samples required at each leaf node\n",
    "min_samples_leaf = [1, 2, 4]\n",
    "# Method of selecting samples for training each tree\n",
    "bootstrap = [True, False]\n",
    "# Create the random grid\n",
    "random_grid = {'n_estimators': n_estimators,\n",
    "               'max_features': max_features,\n",
    "               'max_depth': max_depth,\n",
    "               'min_samples_split': min_samples_split,\n",
    "               'min_samples_leaf': min_samples_leaf,\n",
    "               'bootstrap': bootstrap,\n",
    "               'class_weight':['balanced', 'balanced_subsample', None]}\n",
    "\n",
    "rf_random = RandomizedSearchCV(estimator = rf, param_distributions = random_grid, n_iter = 100, cv = 5, \n",
    "                               scoring=\"f1_macro\", verbose=1, random_state=42, n_jobs = -1)\n",
    "rf_random.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(rf_random.best_params_)\n",
    "print(rf_random.best_score_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the parameter grid based on the results of random search \n",
    "param_grid = {\n",
    "    'bootstrap': [False],\n",
    "    'max_depth': [20],\n",
    "    'max_features': ['auto'],\n",
    "    'min_samples_leaf': [4, 10, 15],\n",
    "    'min_samples_split': [10, 20, 30],\n",
    "    'n_estimators': [20],\n",
    "    'class_weight':['balanced_subsample']\n",
    "}\n",
    "\n",
    "# Instantiate the grid search model\n",
    "grid_search = GridSearchCV(estimator = rf, param_grid = param_grid, \n",
    "                          cv = 5, n_jobs = -1, verbose = 1, scoring='f1_macro')\n",
    "\n",
    "grid_search.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(grid_search.best_params_)\n",
    "print(grid_search.best_score_)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = RandomForestClassifier(n_estimators = 200, \n",
    "                             min_samples_split = 10,\n",
    "                             min_samples_leaf = 4,\n",
    "                             max_features = 'auto',\n",
    "                             max_depth = 20,\n",
    "                             class_weight  = 'balanced_subsample',\n",
    "                             bootstrap = False)\n",
    "clf.fit(X_train, y_train)\n",
    "print(f1_score(y_pred=clf.predict(X_test), y_true=y_test, average=\"macro\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Train on all data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf.fit(X, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Predict on the real test set"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"features has been extracted from the test set the same way than the train set\n",
    "   we don't show it again here\n",
    "\"\"\"\n",
    "mat2 = np.load('features/MMD_test.npz')['arr_0']\n",
    "mat3 = np.load('features/M2_test.npz')['arr_0']\n",
    "mat4 = np.load('features/stats_test_mean.npz')['arr_0']\n",
    "mat5 = np.load('features/energy_test.npz')['arr_0']\n",
    "mat6 = np.load('features/stats_accel_test.npz')['arr_0']\n",
    "np.nan_to_num(mat6, copy=False)\n",
    "mat7 = np.load('features/stats_on_filter_eeg_1_test.npz')['arr_0']\n",
    "mat8 = np.load('features/stats_on_filter_eeg_2_test.npz')['arr_0']\n",
    "mat9 = np.load('features/stats_on_filter_eeg_3_test.npz')['arr_0']\n",
    "mat10 = np.load('features/stats_on_filter_eeg_4_test.npz')['arr_0']\n",
    "mat11 = np.load('features/stats_on_filter_eeg_5_test.npz')['arr_0']\n",
    "mat12 = np.load('features/stats_on_filter_eeg_6_test.npz')['arr_0']\n",
    "mat13 = np.load('features/stats_on_filter_eeg_7_test.npz')['arr_0']\n",
    "\n",
    "X_true_test = np.hstack((mat2, mat3, mat4, mat5, mat6, mat7, mat8, mat9, mat10, mat11, mat12, mat13))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from utils.submit import Submit\n",
    "res = Submit(clf.predict(X_true_test))\n",
    "res.submit('submissions/to_submit.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "results on submit : 0.63203"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Annexe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class pour comparer les algorithmes de classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare Algorithms\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn import model_selection\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.naive_bayes import GaussianNB\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier, BaggingClassifier\n",
    "\n",
    "\n",
    "class CompareClassifiers:\n",
    "    def __init__(self, X, y):\n",
    "        self.X = X\n",
    "        self.y = y\n",
    "\n",
    "    def run(self):\n",
    "        # prepare configuration for cross validation test harness\n",
    "        seed = 7\n",
    "        # prepare models\n",
    "        models = []\n",
    "        models.append(('LR', LogisticRegression(random_state=0, solver='lbfgs')))\n",
    "        models.append(('BAG', BaggingClassifier(n_estimators=20, n_jobs=-1)))\n",
    "        models.append(('GB', GradientBoostingClassifier(n_estimators=20,)))\n",
    "        models.append(('SVC', SVC(gamma='auto')))\n",
    "        models.append(('RF', RandomForestClassifier(n_estimators=100, random_state=0, n_jobs = -1, bootstrap=True, max_features='sqrt')))\n",
    "        \n",
    "        # evaluate each model in turn\n",
    "        results = []\n",
    "        names = []\n",
    "        scoring = 'f1_macro'\n",
    "        for name, model in models:\n",
    "            print(\"running for %s\" % (name))\n",
    "            kfold = model_selection.KFold(n_splits=10, random_state=seed)\n",
    "            cv_results = model_selection.cross_val_score(model, self.X, self.y, cv=kfold, scoring=scoring, verbose=2)\n",
    "            results.append(cv_results)\n",
    "            names.append(name)\n",
    "            msg = \"%s: %f (%f)\" % (name, cv_results.mean(), cv_results.std())\n",
    "            print(msg)\n",
    "        # boxplot algorithm comparison\n",
    "        fig = plt.figure()\n",
    "        fig.suptitle('Algorithm Comparison')\n",
    "        ax = fig.add_subplot(111)\n",
    "        ax.set(ylabel='F1 score')\n",
    "        plt.boxplot(results)\n",
    "        ax.set_xticklabels(names)\n",
    "        plt.show()\n",
    "        fig.savefig('compare_clf.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Class pour créer le csv de submit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "class Submit:\n",
    "    def __init__(self, y_pred, y_true=None):\n",
    "        self.y_pred = y_pred\n",
    "        self.y_true = y_true\n",
    "        self.results = pd.DataFrame()\n",
    "\n",
    "    def score(self, average='weighted'):\n",
    "        if self.y_true == None:\n",
    "            raise ValueError(\"true labels not found\")\n",
    "        return f1_score(self.y_true, self.y_pred, average=average)\n",
    "\n",
    "    def submit(self, filename):\n",
    "        res = pd.DataFrame()\n",
    "        res['id'] = range(len(self.y_pred))\n",
    "        res['sleep_stage'] = self.y_pred\n",
    "        res.set_index('id').to_csv(filename)\n",
    "        return \"csv created at %s\" % filename"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}