3_eeg_T120D101_ML_2_draft_for_1_patient.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def warn(*args, **kwargs):\n",
    "    pass\n",
    "import warnings\n",
    "warnings.warn = warn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load necessary libraries\n",
    "import mne\n",
    "from mne.decoding import Vectorizer\n",
    "import os\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.model_selection import cross_val_score, train_test_split, GridSearchCV, StratifiedKFold\n",
    "from sklearn.metrics import classification_report, accuracy_score, precision_recall_fscore_support\n",
    "\n",
    "# Models\n",
    "from sklearn import svm\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.ensemble import RandomForestClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read the raw data\n",
    "raw_no_med = mne.io.read_raw_fif(\"raw_no_med.fif\", preload=True)\n",
    "raw_with_med = mne.io.read_raw_fif(\"raw_with_med.fif\", preload=True)\n",
    "\n",
    "#Check information\n",
    "display(raw_no_med.info)\n",
    "display(raw_with_med.info)\n",
    "\n",
    "#Plot the data\n",
    "raw_no_med.plot();\n",
    "raw_with_med.plot();\n",
    "\n",
    "# Define the sampling frequency (same for both datasets)\n",
    "fs = raw_no_med.info['sfreq']\n",
    "\n",
    "# Create fixed length epochs\n",
    "epochs_no_med = mne.make_fixed_length_epochs(raw_no_med, duration=5.0, overlap=1.0, preload=True)\n",
    "epochs_with_med = mne.make_fixed_length_epochs(raw_with_med, duration=5.0, overlap=1.0, preload=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"\\n# Define the length of an epoch in seconds\\nepoch_length = 5.0  # for example\\n\\n# Calculate the number of samples per epoch\\nsamples_per_epoch = int(epoch_length * raw_without_medication.info['sfreq'])\\n\\n# Create an array of event times\\nevent_times = np.arange(0, len(raw_without_medication.times), samples_per_epoch)\\n\\n# Create an array of event ids\\nevent_ids = np.ones(len(event_times))\\n\\n# Combine the event times and ids into a single events array\\nevents = np.column_stack([event_times, np.zeros(len(event_times)), event_ids]).astype(int)\\n\\n# Now you can create the epochs\\nepochs_without_medication = mne.Epochs(raw_without_medication, events=events, tmin=0, tmax=epoch_length - 1/raw_without_medication.info['sfreq'], baseline=None)\\n\\nepochs_with_medication = mne.Epochs(raw_with_medication, events=events, tmin=0, tmax=epoch_length - 1/raw_without_medication.info['sfreq'], baseline=None)\\n\""
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "# Define the length of an epoch in seconds\n",
    "epoch_length = 5.0  # for example\n",
    "\n",
    "# Calculate the number of samples per epoch\n",
    "samples_per_epoch = int(epoch_length * raw_no_med.info['sfreq'])\n",
    "\n",
    "# Create an array of event times\n",
    "event_times = np.arange(0, len(raw_no_med.times), samples_per_epoch)\n",
    "\n",
    "# Create an array of event ids\n",
    "event_ids = np.ones(len(event_times))\n",
    "\n",
    "# Combine the event times and ids into a single events array\n",
    "events = np.column_stack([event_times, np.zeros(len(event_times)), event_ids]).astype(int)\n",
    "\n",
    "# Now you can create the epochs\n",
    "epochs_without_medication = mne.Epochs(raw_no_med, events=events, tmin=0, tmax=epoch_length - 1/rraw_no_med.info['sfreq'], baseline=None)\n",
    "\n",
    "epochs_with_medication = mne.Epochs(raw_with_med, events=events, tmin=0, tmax=epoch_length - 1/raw_with_med.info['sfreq'], baseline=None)\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_without_medication = mne.make_fixed_length_epochs(raw_without_medication,duration=5,overlap=0)\n",
    "\n",
    "epochs_with_medication = mne.make_fixed_length_epochs(raw_with_medication,duration=5,overlap=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Now you can create the epochs with the first sample as the baseline\n",
    "#epochs_without_medication = mne.Epochs(raw_without_medication, events=events, tmin=0, tmax=epoch_length - 1/raw_without_medication.info['sfreq'], baseline=(0, 0))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Split the data into epochs\n",
    "#epochs_without_medication = mne.Epochs(raw_without_medication, events=None, tmin=-0.1, tmax=0.8)\n",
    "#epochs_with_medication = mne.Epochs(raw_with_medication, events=None, tmin=-0.1, tmax=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(48, 63, 5120)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs_without_medication.get_data().shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "epochs_without_medication.events[:,-1][epochs_without_medication.events[:,-1] == 1] = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(48,)"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "epochs_without_medication.events[:,-1].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "train_data_UN, test_data_UN, labels_train_UN, labels_test_UN = train_test_split(data_UN, labels_UN, test_size=0.3, random_state=42)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load necessary libraries\n",
    "import mne\n",
    "from mne.decoding import Vectorizer\n",
    "\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.model_selection import cross_val_score, train_test_split, GridSearchCV, StratifiedKFold\n",
    "from sklearn.metrics import classification_report, accuracy_score, precision_recall_fscore_support\n",
    "\n",
    "# Models\n",
    "from sklearn import svm\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "\n",
    "# Logistic Regression\n",
    "clf_lr_pip = make_pipeline(Vectorizer(), StandardScaler(), LogisticRegression(random_state=42))\n",
    "parameters = {'logisticregression__penalty':['l1', 'l2']}\n",
    "gs_cv_lr = GridSearchCV(clf_lr_pip, parameters, scoring='accuracy')\n",
    "gs_cv_lr.fit(train_data_UN, labels_train_UN)\n",
    "\n",
    "print('Best Parameters: {}'.format(gs_cv_lr.best_params_))\n",
    "print('Best Score: {}'.format(gs_cv_lr.best_score_))\n",
    "\n",
    "#Predictions\n",
    "predictions_lr = gs_cv_lr.predict(test_data_UN)\n",
    "\n",
    "#Evaluation\n",
    "report_lr = classification_report(labels_test_UN, predictions_lr, target_names=['Unpleasant', 'Neutral'])\n",
    "print('LR Clasification Report:\\n {}'.format(report_lr))\n",
    "\n",
    "acc_lr = accuracy_score(labels_test_UN, predictions_lr)\n",
    "print(\"Accuracy of LR model: {}\".format(acc_lr))\n",
    "\n",
    "precision_lr,recall_lr,fscore_lr,support_lr=precision_recall_fscore_support(labels_test_UN,predictions_lr,average='macro')\n",
    "print('Precision: {0}, Recall: {1}, f1-score:{2}'.format(precision_lr,recall_lr,fscore_lr))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Features extraction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Power Spectral Density (PSD)\n",
    "psds_without_medication, freqs = mne.time_frequency.psd_array_welch(epochs_without_medication.get_data(), sfreq=raw_without_medication.info['sfreq'], fmin=1., fmax=40.)\n",
    "psds_with_medication, freqs = mne.time_frequency.psd_array_welch(epochs_with_medication.get_data(), sfreq=raw_with_medication.info['sfreq'], fmin=1., fmax=40.)\n",
    "\n",
    "psds_without_medication.shape\n",
    "\n",
    "# Flatten the PSDs\n",
    "psds_without_medication = psds_without_medication.reshape(len(psds_without_medication), -1)\n",
    "psds_with_medication = psds_with_medication.reshape(len(psds_with_medication), -1)\n",
    "\n",
    "# Combine the features and labels\n",
    "features = np.concatenate((psds_without_medication, psds_with_medication))\n",
    "labels = np.concatenate(([0] * len(psds_without_medication), [1] * len(psds_with_medication)))\n",
    "\n",
    "\n",
    "\n",
    "# Extract features\n",
    "#features_without_medication = epochs_without_medication.get_data().mean(axis=2)\n",
    "#features_with_medication = epochs_with_medication.get_data().mean(axis=2)\n",
    "\n",
    "# Create labels\n",
    "#labels_without_medication = [0] * len(features_without_medication)\n",
    "#labels_with_medication = [1] * len(features_with_medication)\n",
    "\n",
    "# Combine the features and labels\n",
    "#features = np.concatenate((features_without_medication, features_with_medication))\n",
    "#labels = np.concatenate((labels_without_medication, labels_with_medication))\n",
    "\n",
    "\n",
    "\n",
    "# Split the data into training and test sets\n",
    "features_train, features_test, labels_train, labels_test = train_test_split(features, labels, test_size=0.33, stratify=labels, shuffle=True)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(48, 630)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psds_without_medication.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import welch\n",
    "import pywt\n",
    "import mne\n",
    "#from mne.time_frequency import tfr_morlet, psd_multitaper, psd_welch\n",
    "\n",
    "# Helper function to calculate Hjorth parameters\n",
    "def hjorth_params(data):\n",
    "    first_deriv = np.diff(data, axis=-1)\n",
    "    second_deriv = np.diff(first_deriv, axis=-1)\n",
    "    \n",
    "    var_zero = np.mean(data ** 2, axis=-1)\n",
    "    var_d1 = np.mean(first_deriv ** 2, axis=-1)\n",
    "    var_d2 = np.mean(second_deriv ** 2, axis=-1)\n",
    "    \n",
    "    activity = var_zero\n",
    "    mobility = np.sqrt(var_d1 / var_zero)\n",
    "    complexity = np.sqrt(var_d2 / var_d1) / mobility\n",
    "    return np.vstack([activity, mobility, complexity]).T\n",
    "\n",
    "# Expanded feature extraction function\n",
    "def extract_features(epochs):\n",
    "    features = []\n",
    "    \n",
    "    # Extract data from MNE epochs object\n",
    "    data = epochs.get_data()\n",
    "    \n",
    "    # Power Spectral Density (PSD)\n",
    "    #psds, freqs = psd_multitaper(epochs, fmin=1, fmax=40)\n",
    "    #psds = psds.reshape(len(psds), -1)  # Flatten the PSDs\n",
    "    #features.append(psds)\n",
    "    \n",
    "    # Wavelet Transform (WT)\n",
    "    coeffs = pywt.wavedec(data, wavelet='db4', level=5, axis=-1)\n",
    "    wt_features = np.array([c.ravel() for c in coeffs]).reshape(len(data), -1)\n",
    "    features.append(wt_features)\n",
    "    \n",
    "    # Time-Frequency Distributions (TFD)\n",
    "    freqs = np.arange(1, 40, 2)\n",
    "    n_cycles = freqs / 2.\n",
    "    power = tfr_array_morlet(data, sfreq=epochs.info['sfreq'],\n",
    "                             freqs=freqs, n_cycles=n_cycles,\n",
    "                             output='power')\n",
    "    tfd_features = power.reshape(len(data), -1)\n",
    "    features.append(tfd_features)\n",
    "    \n",
    "    # Auto Regressive Methods (ARM)\n",
    "    # Implement ARM feature extraction here (e.g., using Yule-Walker equations)\n",
    "    \n",
    "    # Hjorth Parameters\n",
    "    hjorth_features = hjorth_params(data)\n",
    "    features.append(hjorth_features)\n",
    "    \n",
    "    # Combine all features\n",
    "    combined_features = np.concatenate(features, axis=1)\n",
    "    \n",
    "    return combined_features\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (6,) + inhomogeneous part.",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_25008\\3643545301.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0mextract_features\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mepochs_with_medication\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;32m~\\AppData\\Local\\Temp\\ipykernel_25008\\2842907425.py\u001b[0m in \u001b[0;36mextract_features\u001b[1;34m(epochs)\u001b[0m\n\u001b[0;32m     33\u001b[0m     \u001b[1;31m# Wavelet Transform (WT)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     34\u001b[0m     \u001b[0mcoeffs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mpywt\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mwavedec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwavelet\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;34m'db4'\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mlevel\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m5\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 35\u001b[1;33m     \u001b[0mwt_features\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mc\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mc\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mcoeffs\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mreshape\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdata\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m-\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     36\u001b[0m     \u001b[0mfeatures\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwt_features\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     37\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (6,) + inhomogeneous part."
     ]
    }
   ],
   "source": [
    "extract_features(epochs_with_medication)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.95036023e-12 1.18768722e-12 9.58905346e-13 ... 3.40992924e-15\n",
      "  4.99635282e-15 3.06046996e-15]\n",
      " [2.10727097e-12 3.12391312e-12 1.37626678e-12 ... 7.92290425e-15\n",
      "  2.94275846e-15 6.18057075e-16]\n",
      " [1.32228334e-12 1.17702013e-12 1.15200473e-12 ... 8.78299922e-15\n",
      "  3.52778775e-15 3.23877256e-16]\n",
      " ...\n",
      " [2.07861209e-12 2.00458775e-12 1.69308339e-12 ... 3.30442721e-14\n",
      "  1.75636294e-14 2.35121609e-15]\n",
      " [1.62115564e-12 3.75501098e-12 1.69829903e-12 ... 3.21813523e-14\n",
      "  1.32846982e-14 1.99658105e-15]\n",
      " [2.51829578e-12 4.06714988e-12 1.34456959e-12 ... 3.21751190e-14\n",
      "  1.42002684e-14 2.96830228e-15]]\n"
     ]
    }
   ],
   "source": [
    "print(features)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "96\n"
     ]
    }
   ],
   "source": [
    "print(len(labels))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Machine learning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf_svm_0 = make_pipeline(Vectorizer(), StandardScaler(), svm.SVC(kernel='rbf', C=1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Accuracy of 1th fold is 0.9\n",
      "\n",
      "Accuracy of 2th fold is 1.0\n",
      "\n",
      "Accuracy of 3th fold is 0.9473684210526315\n",
      "\n",
      "Accuracy of 4th fold is 0.9473684210526315\n",
      "\n",
      "Accuracy of 5th fold is 1.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "clf_svm_0 = make_pipeline(Vectorizer(), StandardScaler(), svm.SVC(kernel='rbf', C=1))\n",
    "scores = cross_val_score(clf_svm_0, features, labels, cv=5)\n",
    "for i in range(len(scores)):   \n",
    "    print('Accuracy of ' + str(i+1) + 'th fold is ' + str(scores[i]) + '\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#svm\n",
    "clf_svm_pip = make_pipeline(Vectorizer(), StandardScaler(), svm.SVC(random_state=42))\n",
    "parameters = {'svc__kernel':['linear', 'rbf', 'sigmoid'], 'svc__C':[0.1, 1, 10]}\n",
    "gs_cv_svm = GridSearchCV(clf_svm_pip, parameters, scoring='accuracy', cv=StratifiedKFold(n_splits=5), return_train_score=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best Parameters: {'svc__C': 0.1, 'svc__kernel': 'linear'}\n",
      "Best Score: 1.0\n"
     ]
    }
   ],
   "source": [
    "gs_cv_svm.fit(features_train, labels_train)\n",
    "print('Best Parameters: {}'.format(gs_cv_svm.best_params_))\n",
    "print('Best Score: {}'.format(gs_cv_svm.best_score_))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SVM Clasification Report:\n",
      "               precision    recall  f1-score   support\n",
      "\n",
      "  Unpleasant       1.00      1.00      1.00        16\n",
      "     Neutral       1.00      1.00      1.00        16\n",
      "\n",
      "    accuracy                           1.00        32\n",
      "   macro avg       1.00      1.00      1.00        32\n",
      "weighted avg       1.00      1.00      1.00        32\n",
      "\n",
      "Accuracy of SVM model: 1.0\n",
      "Precision: 1.0, Recall: 1.0, f1-score:1.0\n"
     ]
    }
   ],
   "source": [
    "#Prediction\n",
    "predictions_svm = gs_cv_svm.predict(features_test)\n",
    "\n",
    "#Evaluate\n",
    "report_svm = classification_report(labels_test, predictions_svm, target_names=['Unpleasant', 'Neutral'])\n",
    "print('SVM Clasification Report:\\n {}'.format(report_svm))\n",
    "\n",
    "acc_svm = accuracy_score(labels_test, predictions_svm)\n",
    "print(\"Accuracy of SVM model: {}\".format(acc_svm))\n",
    "\n",
    "precision_svm,recall_svm,fscore_svm,support_svm=precision_recall_fscore_support(labels_test,predictions_svm,average='macro')\n",
    "print('Precision: {0}, Recall: {1}, f1-score:{2}'.format(precision_svm,recall_svm,fscore_svm))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best Parameters: {'logisticregression__penalty': 'l2'}\n",
      "Best Score: 1.0\n",
      "LR Clasification Report:\n",
      "               precision    recall  f1-score   support\n",
      "\n",
      "  Unpleasant       1.00      1.00      1.00        16\n",
      "     Neutral       1.00      1.00      1.00        16\n",
      "\n",
      "    accuracy                           1.00        32\n",
      "   macro avg       1.00      1.00      1.00        32\n",
      "weighted avg       1.00      1.00      1.00        32\n",
      "\n",
      "Accuracy of LR model: 1.0\n",
      "Precision: 1.0, Recall: 1.0, f1-score:1.0\n"
     ]
    }
   ],
   "source": [
    "# Logistic Regression\n",
    "clf_lr_pip = make_pipeline(Vectorizer(), StandardScaler(), LogisticRegression(random_state=42))\n",
    "parameters = {'logisticregression__penalty':['l1', 'l2']}\n",
    "gs_cv_lr = GridSearchCV(clf_lr_pip, parameters, scoring='accuracy')\n",
    "gs_cv_lr.fit(features_train, labels_train)\n",
    "\n",
    "print('Best Parameters: {}'.format(gs_cv_lr.best_params_))\n",
    "print('Best Score: {}'.format(gs_cv_lr.best_score_))\n",
    "\n",
    "#Predictions\n",
    "predictions_lr = gs_cv_lr.predict(features_test)\n",
    "\n",
    "#Evaluation\n",
    "report_lr = classification_report(labels_test, predictions_lr, target_names=['Unpleasant', 'Neutral'])\n",
    "print('LR Clasification Report:\\n {}'.format(report_lr))\n",
    "\n",
    "acc_lr = accuracy_score(labels_test, predictions_lr)\n",
    "print(\"Accuracy of LR model: {}\".format(acc_lr))\n",
    "\n",
    "precision_lr,recall_lr,fscore_lr,support_lr=precision_recall_fscore_support(labels_test,predictions_lr,average='macro')\n",
    "print('Precision: {0}, Recall: {1}, f1-score:{2}'.format(precision_lr,recall_lr,fscore_lr))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LDA Clasification Report:\n",
      "               precision    recall  f1-score   support\n",
      "\n",
      "  Unpleasant       1.00      1.00      1.00        16\n",
      "     Neutral       1.00      1.00      1.00        16\n",
      "\n",
      "    accuracy                           1.00        32\n",
      "   macro avg       1.00      1.00      1.00        32\n",
      "weighted avg       1.00      1.00      1.00        32\n",
      "\n",
      "Accuracy of LDA model: 1.0\n",
      "Precision: 1.0, Recall: 1.0, f1-score:1.0\n"
     ]
    }
   ],
   "source": [
    "# Linear Discriminant Analysis\n",
    "clf_lda_pip = make_pipeline(Vectorizer(), StandardScaler(), LinearDiscriminantAnalysis(solver='svd'))\n",
    "clf_lda_pip.fit(features_train,labels_train)\n",
    "\n",
    "#Predictions\n",
    "predictions_lda = clf_lda_pip.predict(features_test)\n",
    "\n",
    "#Evaluation\n",
    "report_lda = classification_report(labels_test, predictions_lda, target_names=['Unpleasant', 'Neutral'])\n",
    "print('LDA Clasification Report:\\n {}'.format(report_lda))\n",
    "\n",
    "acc_lda = accuracy_score(labels_test, predictions_lda)\n",
    "print(\"Accuracy of LDA model: {}\".format(acc_lda))\n",
    "\n",
    "precision_lda,recall_lda,fscore_lda,support_lda=precision_recall_fscore_support(labels_test,predictions_lda,average='macro')\n",
    "print('Precision: {0}, Recall: {1}, f1-score:{2}'.format(precision_lda,recall_lda,fscore_lda))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "accuracies, f1_scores = [], []\n",
    "accuracies.append([acc_svm, acc_lr, acc_lda])\n",
    "f1_scores.append([fscore_svm, fscore_lr, fscore_lda])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotEvalMetrics(tasks, labels, evalMetric, metricName):\n",
    "    width = 0.2  # the width of the bars\n",
    "\n",
    "    # Set position of bar on X axis\n",
    "    rects1 = np.arange(len(evalMetric[:][0]))\n",
    "    rects2 = [x + width for x in rects1]\n",
    "    rects3 = [x + width for x in rects2]\n",
    "\n",
    "    plt.bar(rects1, list(zip(*evalMetric))[0], color='#87CEFA', width=width, edgecolor='white', label=labels[0])\n",
    "    plt.bar(rects2, list(zip(*evalMetric))[1], color='#FFE4E1', width=width, edgecolor='white', label=labels[1])\n",
    "    plt.bar(rects3, list(zip(*evalMetric))[2], color='#CD5C5C', width=width, edgecolor='white', label=labels[2])\n",
    "\n",
    "    plt.xlabel('Classification Tasks')\n",
    "    plt.xticks([r + width for r in range(len(evalMetric[:][0]))], tasks)\n",
    "    plt.ylabel(metricName)\n",
    "\n",
    "    plt.legend(bbox_to_anchor=(1.01, 1), loc='upper left', )\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcQAAAEGCAYAAAD7f+WgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAAsTAAALEwEAmpwYAAAZgklEQVR4nO3de5RV5Znn8e9TlAhqwkVKRUAxUTFENApjYmfoTkx6gh2DtpCODEquokkbpzXR2CuZjE23uXZMVlpsRZOoiCGoa9IYMWS8jRrHtCAqoLElBqMIioKiAuKxnvnjnEqfFHU5ddlVVPH9rMXy7H3e8+7n1F7Wr9533yIzkSRpd1fX2wVIkrQrMBAlScJAlCQJMBAlSQIMREmSAKjv7QI6asSIETl27NjeLkOS+pTly5e/mJkNvV3HrqzPBeLYsWNZtmxZb5chSX1KRDzd2zXs6pwylSQJA1GSJMBAlCQJ6IPHECVJ3W/58uX71dfXXw0cSf8cLDUCq0ql0ucmTpz4QksNDERJEvX19VcfcMAB72poaNhcV1fX725y3djYGBs3bhy/YcOGq4GpLbXpj38FSJI67siGhoYt/TEMAerq6rKhoeEVyiPgltsUtfGI+HFEvBARq1p5PyLihxGxJiIejYhji6pFktSuuv4ahk0q36/V3CtyhHgNMKWN908EDqv8mw38a4G1SJLUpsKOIWbmPRExto0mJwPXZfmBjA9ExNCIGJmZ64uqSZJUm1JjHl1fF92WEaXGLNXXxSPttfvKV75ywM0337xvXV1d1tXV8dGPfnTz9u3b6+bOnbuuqc39998/+PTTT3/HU089tXrUqFETDjjggB3Lly9/oun9I444Yvxbb70VTz755OqO1NibJ9WMAp6pWn62sm6nQIyI2ZRHkRx00EGd3mCpMamvi05/vj3Z2EjUFTPobtyxg7qBAwvpuyf6L4r7tPf6L4r7tPf6r1ZfF/XfWvFmt/V30TF7tJs3t99++95Lly4dunLlyscGDx6c69evr3/44YcHnXnmmYdUB+L1118//NRTT93UtPz6668PWLNmzR6HHnromw899NCgztbYJ84yzcx5wDyASZMmdXqOu74u6M4d3NxFx+wB6/5QSN91ow7ioTPOKKRvgGPnzy+s7yK5T1vnPm2Z+3TXtW7duj2GDx9eGjx4cAKMHDmyNHLkyNeGDBlSuvPOO/c+4YQTXgdYvHjx8Ntuu+0/mj53yimnbLruuuuGz5kz5/nrrrtu+LRp0zYtWrRo345uvzfPMl0HjKlaHl1ZJ0naDZ1yyilbnnvuuYFjx4498vTTTz/o1ltv3Qdg2rRpmxYsWDAc4I477th76NChpQkTJrzR9LkZM2ZsvuWWW4YBLF26dOipp576cme235uBuBiYVTnb9H3AKx4/lKTd15AhQxpXrVr12GWXXfZ0Q0ND6ZOf/OQ7f/jDH+47a9asTbfeeuuwt956iwULFgyfNm3apurP7bfffm8NGTKkNG/evGGHHnrotn322aexM9svbMo0In4KfAAYERHPAv8L2AMgM68AlgB/BawBtgKfLqoWSVLfUF9fz0knnfTqSSed9OpRRx21bf78+fuee+65L40ePfqNJUuWvG3JkiXDfv3rXz/e/HPTp0/ffOGFFx58+eWX/77T2+5a6a3LzBntvJ/A3xa1fUlS3/LII4/sWVdXR9N06IoVKwaPHj16B8DHP/7xTRdccMGYMWPGvPHOd75zp4PMM2fO3Lx+/fo9Tj311C1PP/30Hp3Zfp84qUaS1LNKjVmq5czQjvTX3tnDW7ZsGXDuuecetGXLlgEDBgzIsWPHvnHttdc+DTBr1qzNX/3qV8d84xvfeKalzw4bNqzxkksu2dCVGg1ESdJOarlmsIP9tdtm8uTJW1esWPHblt4bOXJkqVQqPdR8/bp161Y2Xzdu3LgdHb0GEbyXqSRJgIEoSRJgIEqSBBiIkiQBBqIkSYCBKEkS4GUXkqSWZB5NdN/jn8gsEW1fyrHXXnsds3Xr1hXV684///wDr7/++hHDhw8vvfnmm3HhhReuP+ussza11kdXGIiSpJ1F1HfrU0FGHdTpvDn77LOfnzNnzvMrV67c8/jjjx//qU99avOee+7Z6ScftcYpU0lSnzBhwoQ3Bg0a1Pjiiy8OKKJ/A1GS1Cfcd999ex188MHbR40aVSqif6dMJUm7tCuuuGL/G264YcTatWv3XLhw4ZqituMIUZK0Szv77LOfX7Nmzeprr732d1/4whfGbt26tf0bo3aCgShJ6hNmzpz5yoQJE16fO3fuvkX075SpJGlnmaWunBnaYn/R9sBu+/btdfvvv/9RTcuf//znn2/e5uKLL15/xhlnvOP8889/ccCA7j23xkCUJO2snWsGO9Ffu00aGxuXt9dm8uTJW9euXbuqW2pqxilTSZIwECVJAgxESZIAA1GSJMBAlCQJMBAlSQK87EKS1ILGHTuOrhs4sNsyonHHjlLdwIFdevzTtm3b6saNG7ftm9/85rqJEydub2pz//33D37/+98//sYbb3xy+vTpWzpbo4EoSdpJ3cCB9Q+dcUa39Xfs/PldfvwTwFVXXTXsIx/5yLhHH3109YEHHlgCmD9//vBjjz32tRtuuGF4VwLRKVNJUp9x5plnbp48efIrP/rRj4YDNDY2cssttwy/7rrr1t53331v78p9Tg1ESVKfcswxx2z97W9/Owjg9ttv33vMmDFvvPvd737jve9976uLFi0a0tl+DURJUp+SmX98ff311+87ffr0TQCnnXbapoULFw7vbL8eQ5Qk9SkPP/zwXhMnTtxaKpW47bbbhv7qV78aeumll47MTF5++eX6zZs31w0bNqyxo/06QpQk9RnXXHPN0HvvvXfIZz7zmU2LFy9++7hx47Zt2LDh0XXr1q187rnnVk6ZMmXzggULhnWmb0eIkqSdNO7YUerKmaEt9Vc3cGCbbVp7/NMVV1yx/6JFi/bdtm1b3eGHH75t6dKlTxx44IGlG264YfjUqVNfru5j2rRpm6+88sr9zjnnnJc6WqOBKEnaSXvXDHaiv3bbtPb4p0svvfS5ltbfdNNNa5uvmzlz5iszZ858paP1QcFTphExJSKeiIg1EXFRC+8fFBF3RcSKiHg0Iv6qyHokSWpNYYEYEQOAucCJwHhgRkSMb9bsa8CizDwGOA24vKh6JElqS5EjxOOANZn5VGbuABYCJzdrk8DbK6+HAC0OiyVJhWtsbGzs9EXtfUHl+7V69mmRgTgKeKZq+dnKumoXA6dHxLPAEuCLLXUUEbMjYllELNu4cWMRtUrS7m7Vxo0bh/TXUGxsbIyNGzcOAVa11qa3T6qZAVyTmd+LiOOB+RFxZGb+SYJn5jxgHsCkSZOyhX4kSV1QKpU+t2HDhqs3bNhwJP3zkrxGYFWpVPpcaw2KDMR1wJiq5dGVddU+C0wByMz/FxGDgBHACwXWJUlqZuLEiS8AU3u7jt5U5F8BDwKHRcQhETGQ8kkzi5u1+QPwIYCIeBcwCHBOVJLU4woLxMwsAecAS4HHKZ9Nujoi5kRE018hXwLOjIhHgJ8Cn8rqm9RJktRDCj2GmJlLKJ8sU73u61WvHwPeX2QNkiTVoj8eOJUkqcMMREmSMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCDERJkgADUZIkwECUJAkwECVJAgxESZIAA1GSJMBAlCQJMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCDERJkgADUZIkwECUJAkwECVJAgxESZIAA1GSJMBAlCQJMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCagjEiPhYRBickqR+rZag+wTwZER8JyKO6EjnETElIp6IiDURcVErbf4mIh6LiNURcUNH+pckqbvUt9cgM0+PiLcDM4BrIiKBnwA/zcxXW/tcRAwA5gJ/CTwLPBgRizPzsao2hwF/D7w/MzdHxH5d+zqSJHVOTVOhmbkFuAlYCIwE/hp4KCK+2MbHjgPWZOZTmbmj8tmTm7U5E5ibmZsr23mhg/VLktQtajmGODUi/jdwN7AHcFxmnggcDXypjY+OAp6pWn62sq7a4cDhEfHriHggIqZ0pHhJkrpLu1OmwDTg+5l5T/XKzNwaEZ/thu0fBnwAGA3cExETMvPl6kYRMRuYDXDQQQd1cZOSJO2slinTi4F/b1qIiMERMRYgM+9o43PrgDFVy6Mr66o9CyzOzDcz8/fAf1AOyD+RmfMyc1JmTmpoaKihZEmSOqaWQLwRaKxafquyrj0PAodFxCERMRA4DVjcrM3PKY8OiYgRlKdQn6qhb0mSulUtgVhfOSkGgMrrge19KDNLwDnAUuBxYFFmro6IORExtdJsKfBSRDwG3AVckJkvdfRLSJLUVbUcQ9wYEVMzczFARJwMvFhL55m5BFjSbN3Xq14ncH7lnyRJvaaWQDwbWBARlwFB+czRWYVWJUlSD6vlwvzfAe+LiH0qy68VXpUkST2slhEiEfFR4N3AoIgAIDPnFFiXJEk9qpYL86+gfD/TL1KeMv04cHDBdUmS1KNqOcv0zzJzFrA5M/8BOJ7y5RGSJPUbtQTi9sp/t0bEgcCblO9nKklSv1HLMcRbImIo8F3gISCBq4osSpKkntZmIFYeDHxH5d6iN0fEL4BBmflKTxQnSVJPaXPKNDMbKT/TsGn5DcNQktQf1XIM8Y6ImBZN11tIktQP1RKIZ1G+mfcbEbElIl6NiC0F1yVJUo+q5U41b+uJQiRJ6k3tBmJE/HlL65s/MFiSpL6slssuLqh6PQg4DlgOnFBIRZIk9YJapkw/Vr0cEWOAHxRVkCRJvaGWk2qaexZ4V3cXIklSb6rlGOK/UL47DZQD9D2U71gjSVK/UcsxxGVVr0vATzPz1wXVI0lSr6glEG8CtmfmWwARMSAi9srMrcWWJklSz6npTjXA4KrlwcDtxZQjSVLvqCUQB2Xma00Lldd7FVeSJEk9r5ZAfD0ijm1aiIiJwLbiSpIkqefVcgzx74AbI+I5IIADgE8UWZQkST2tlgvzH4yII4BxlVVPZOabxZYlSVLPanfKNCL+Ftg7M1dl5ipgn4j4QvGlSZLUc2o5hnhmZr7ctJCZm4EzC6tIkqReUEsgDqh+OHBEDAAGFleSJEk9r5aTan4J/CwirqwsnwXcVlxJkiT1vFoC8SvAbODsyvKjlM80lSSp32h3yjQzG4HfAGspPwvxBODxYsuSJKlntTpCjIjDgRmVfy8CPwPIzA/2TGmSJPWctqZMfwvcC5yUmWsAIuK8HqlKkqQe1taU6anAeuCuiLgqIj5E+U41kiT1O60GYmb+PDNPA44A7qJ8C7f9IuJfI+K/1dJ5REyJiCciYk1EXNRGu2kRkRExqYP1S5LULWo5qeb1zLwhMz8GjAZWUD7ztE2V6xXnAicC44EZETG+hXZvA/4H5RN3JEnqFbVcmP9Hmbk5M+dl5odqaH4csCYzn8rMHcBC4OQW2v0j8G1ge0dqkSSpO3UoEDtoFPBM1fKzlXV/VHms1JjMvLWtjiJidkQsi4hlGzdu7P5KJUm7vSIDsU0RUQdcCnypvbaVUemkzJzU0NBQfHGSpN1OkYG4DhhTtTy6sq7J24AjgbsjYi3wPmCxJ9ZIknpDkYH4IHBYRBwSEQOB04DFTW9m5iuZOSIzx2bmWOABYGpmLiuwJkmSWlRYIGZmCTgHWEr5Vm+LMnN1RMyJiKlFbVeSpM6o5ebenZaZS4AlzdZ9vZW2HyiyFkmS2tJrJ9VIkrQrMRAlScJAlCQJMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCDERJkgADUZIkwECUJAkwECVJAgxESZIAA1GSJMBAlCQJMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCDERJkgADUZIkwECUJAkwECVJAgxESZIAA1GSJMBAlCQJMBAlSQIMREmSAANRkiTAQJQkCTAQJUkCCg7EiJgSEU9ExJqIuKiF98+PiMci4tGIuCMiDi6yHkmSWlNYIEbEAGAucCIwHpgREeObNVsBTMrMo4CbgO8UVY8kSW0pcoR4HLAmM5/KzB3AQuDk6gaZeVdmbq0sPgCMLrAeSZJaVWQgjgKeqVp+trKuNZ8FbmvpjYiYHRHLImLZxo0bu7FESZLKdomTaiLidGAS8N2W3s/MeZk5KTMnNTQ09GxxkqTdQn2Bfa8DxlQtj66s+xMR8WHgq8BfZOYbBdYjSVKrihwhPggcFhGHRMRA4DRgcXWDiDgGuBKYmpkvFFiLJEltKiwQM7MEnAMsBR4HFmXm6oiYExFTK82+C+wD3BgRD0fE4la6kySpUEVOmZKZS4AlzdZ9ver1h4vcviRJtdolTqqRJKm3GYiSJGEgSpIEGIiSJAEGoiRJgIEoSRJgIEqSBBiIkiQBBqIkSYCBKEkSYCBKkgQYiJIkAQaiJEmAgShJEmAgSpIEGIiSJAEGoiRJgIEoSRJgIEqSBBiIkiQBBqIkSYCBKEkSYCBKkgQYiJIkAQaiJEmAgShJEmAgSpIEGIiSJAEGoiRJgIEoSRJgIEqSBBiIkiQBBqIkSUDBgRgRUyLiiYhYExEXtfD+nhHxs8r7v4mIsUXWI0lSawoLxIgYAMwFTgTGAzMiYnyzZp8FNmfmocD3gW8XVY8kSW0pcoR4HLAmM5/KzB3AQuDkZm1OBq6tvL4J+FBERIE1SZLUosjMYjqOmA5MyczPVZbPAN6bmedUtVlVafNsZfl3lTYvNutrNjC7sjgOeKKQonvWCODFdlupL3Gf9k/9Zb8enJkNvV3Erqy+twuoRWbOA+b1dh3dKSKWZeak3q5D3cd92j+5X3cfRU6ZrgPGVC2PrqxrsU1E1ANDgJcKrEmSpBYVGYgPAodFxCERMRA4DVjcrM1i4JOV19OBO7OoOVxJktpQ2JRpZpYi4hxgKTAA+HFmro6IOcCyzFwM/AiYHxFrgE2UQ3N30a+mgAW4T/sr9+tuorCTaiRJ6ku8U40kSRiIkiQBBmLhImJs5XrL6nUXR8SXI+KaiFgXEXtW1o+IiLW9Uqg6pIb9+vuIeDgiHoqI43urTtUuIjIivle1/OWIuLjy+uLK/6sPR8SqiJjaa4WqMAZi73sL+ExvF6Fud0Fmvge4CLiyl2tRbd4ATo2IEa28//3KPv048OOI8PdnP+MO7X0/AM6rXIep/uce4NDeLkI1KVE+o/S8thpl5uOVtq0Fp/ooA7H3/QG4DzijtwtRIT4GrOztIlSzucDMiBjSWoOIeC/QCGzssarUIxyVFK+161qq138T+Dfg1uLLUTdpb79+NyK+RvmX5md7piR1VWZuiYjrgHOBbc3ePi8iTgdeBT7hTUT6HwOxeC8Bw5qtGw78vmkhM5+MiIeBv+nButQ17e3XCzLzpp4tSd3kB8BDwE+arf9+Zv5zz5ejnuKUacEy8zVgfUScABARw4EplKdJq10CfLmHy1MndWC/qo/JzE3AIhzZ73YMxJ4xC/iflVHgncA/ZObvqhtk5mrKf5Wq72h3v6rP+h6eNLPb8dZtkiThCFGSJMBAlCQJMBAlSQIMREmSAANRkiTAQNQuJiIOiIiFEfG7iFgeEUsi4vCWni7Rxe3MiYgPV15PjojVlScZjIqITl1QHxGfiogDq5avjojxXazz05W6Ho6IHRGxsvL6Wx3o4+6ImNSVOqTdgZddaJcREQHcD1ybmVdU1h0NvB14BvhFZh5ZwHavAO7LzOu72M/dwJczc1m3FLZz/2uBSZn54q5Ul9RfOELUruSDwJtNYQiQmY9k5r3VjSqjxXsrzxp8KCL+rLJ+ZETcU/XMuskRMaDyfMJVldHVeZW210TE9Ij4HOVb5v1jRCyoHolWPvvPlc8+GhFfrKz/ekQ8WFk/L8qmA5OABZXtD64emUXEjMr2V0XEt6u+y2sRcUlEPBIRD0TE/rX8oCLi55UR9OqImF1V707fteozdZX3/6m9ttLuyHuZaldyJLC8hnYvAH+Zmdsj4jDgp5TD6L8DSzPzkogYAOwFvAcY1TSyjIih1R1l5tUR8V8pjz5vioixVW/PBsYC78nMUuX2bACXZeacSn/zgZMqnz2HqpFYecALlWnUbwMTgc3AryLilMz8ObA38EBmfjUivgOcCfxTDT+Dz2TmpogYDDwYETdXam3tu9YDC4BVlZ/PxLZ+LtLuyBGi+qI9gKsiYiVwI9B0nO5B4NNRfsr5hMx8FXgKeEdE/EtETAG2dGA7HwauzMwS/PEelwAfjIjfVLZ/AvDudvr5L8Ddmbmx0tcC4M8r7+0AflF5vZxyqNXi3Ih4BHgAGAMcRtvf9UoqYVhZ7srPReqXDETtSlZTHkW15zzgeeBoyiPDgQCZeQ/loFkHXBMRszJzc6Xd3cDZwNVdKTAiBgGXA9MzcwJwFTCoC12+WfUYobeoYdYmIj5AOayPz8yjgRXAoHa+6/2Ug3wQQHf/XKT+wEDUruROYM+mY2IAEXFURExu1m4IsD4zGyk/WHlApe3BwPOZeRXlX/DHRsQIoC4zbwa+BhzbgXr+D3BWRNRX+h/Of4bfixGxDzC9qv2rwNta6Offgb+IiBGVqdwZwP/tQB3NDQE2Z+bWiDgCeF+lvra+64+AJcCiiKjv4s9F6pc8hqhdRmZmRPw18IOI+AqwHVgL/F2zppcDN0fELOCXwOuV9R8ALoiIN4HXKD+NYhTwk4ho+uPv7ztQ0tXA4cCjlT6vyszLIuIqYBWwgfI0bZNrgCsiYhtwfNX3Wh8RFwF3AQHcmpn/1oE6mvslcHZEPA48QXnaFNr5rpl5aZSfBD8f+FZbbaXdkZddSJKEU6aSJAEGoiRJgIEoSRJgIEqSBBiIkiQBBqIkSYCBKEkSAP8fMlGgCS9FGjYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[1.0, 1.0, 1.0]]\n"
     ]
    }
   ],
   "source": [
    "#Plot Accuracies\n",
    "tasks = ['UN', 'UP', 'NP']\n",
    "labels = ['SVM', 'LR', 'LDA']\n",
    "plotEvalMetrics(tasks, labels, accuracies, 'Accuracy')\n",
    "print(accuracies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "              precision    recall  f1-score   support\n",
      "\n",
      "           0       1.00      0.19      0.32        16\n",
      "           1       0.55      1.00      0.71        16\n",
      "\n",
      "    accuracy                           0.59        32\n",
      "   macro avg       0.78      0.59      0.51        32\n",
      "weighted avg       0.78      0.59      0.51        32\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Train a classifier\n",
    "clf = RandomForestClassifier(random_state=42)\n",
    "clf.fit(features_train, labels_train)\n",
    "\n",
    "# Evaluate the classifier\n",
    "labels_pred = clf.predict(features_test)\n",
    "print(classification_report(labels_test, labels_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering raw data in 1 contiguous segment\n",
      "Setting up band-pass filter from 1 - 50 Hz\n",
      "\n",
      "FIR filter parameters\n",
      "---------------------\n",
      "Designing a one-pass, zero-phase, non-causal bandpass filter:\n",
      "- Windowed time-domain design (firwin) method\n",
      "- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation\n",
      "- Lower passband edge: 1.00\n",
      "- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)\n",
      "- Upper passband edge: 50.00 Hz\n",
      "- Upper transition bandwidth: 12.50 Hz (-6 dB cutoff frequency: 56.25 Hz)\n",
      "- Filter length: 3381 samples (3.302 sec)\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   2 out of   2 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   4 out of   4 | elapsed:    0.0s remaining:    0.0s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Filtering raw data in 1 contiguous segment\n",
      "Setting up band-pass filter from 1 - 50 Hz\n",
      "\n",
      "FIR filter parameters\n",
      "---------------------\n",
      "Designing a one-pass, zero-phase, non-causal bandpass filter:\n",
      "- Windowed time-domain design (firwin) method\n",
      "- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation\n",
      "- Lower passband edge: 1.00\n",
      "- Lower transition bandwidth: 1.00 Hz (-6 dB cutoff frequency: 0.50 Hz)\n",
      "- Upper passband edge: 50.00 Hz\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done  63 out of  63 | elapsed:    3.0s finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- Upper transition bandwidth: 12.50 Hz (-6 dB cutoff frequency: 56.25 Hz)\n",
      "- Filter length: 3381 samples (3.302 sec)\n",
      "\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n",
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   2 out of   2 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   3 out of   3 | elapsed:    0.0s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   4 out of   4 | elapsed:    0.1s remaining:    0.0s\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Not setting metadata\n",
      "241 matching events found\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done  63 out of  63 | elapsed:    4.2s finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "Not setting metadata\n",
      "241 matching events found\n",
      "No baseline correction applied\n",
      "0 projection items activated\n",
      "Using data from preloaded Raw for 241 events and 1024 original time points ...\n",
      "1 bad epochs dropped\n",
      "Effective window size : 0.250 (s)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using data from preloaded Raw for 241 events and 1024 original time points ...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:   25.5s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:   25.5s finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 bad epochs dropped\n",
      "Effective window size : 0.250 (s)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using data from preloaded Raw for 240 events and 1024 original time points ...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:   15.1s remaining:    0.0s\n",
      "[Parallel(n_jobs=1)]: Done   1 out of   1 | elapsed:   15.1s finished\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using data from preloaded Raw for 240 events and 1024 original time points ...\n",
      "Using data from preloaded Raw for 240 events and 1024 original time points ...\n",
      "Using data from preloaded Raw for 240 events and 1024 original time points ...\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "all the input arrays must have same number of dimensions, but the array at index 0 has 3 dimension(s) and the array at index 2 has 2 dimension(s)",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-27-25fe3b67e116>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m     60\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     61\u001b[0m \u001b[1;31m# Combine the features and labels\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 62\u001b[1;33m \u001b[0mfeatures\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpsds_without_medication\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpsds_with_medication\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmean_without_medication\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mmean_with_medication\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstd_without_medication\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mstd_with_medication\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m     63\u001b[0m \u001b[0mlabels\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mconcatenate\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpsds_without_medication\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpsds_with_medication\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m     64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mconcatenate\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[1;31mValueError\u001b[0m: all the input arrays must have same number of dimensions, but the array at index 0 has 3 dimension(s) and the array at index 2 has 2 dimension(s)"
     ]
    }
   ],
   "source": [
    "import mne\n",
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import classification_report\n",
    "#import pyeeg\n",
    "#import nolds\n",
    "\n",
    "\n",
    "# Preprocess the data\n",
    "raw_without_medication.filter(1., 50., fir_design='firwin')\n",
    "raw_with_medication.filter(1., 50., fir_design='firwin')\n",
    "\n",
    "# Define the length of an epoch in seconds\n",
    "epoch_length = 1.0  # for example\n",
    "\n",
    "# Calculate the number of samples per epoch\n",
    "samples_per_epoch = int(epoch_length * raw_without_medication.info['sfreq'])\n",
    "\n",
    "# Create an array of event times\n",
    "event_times = np.arange(0, len(raw_without_medication.times), samples_per_epoch)\n",
    "\n",
    "# Create an array of event ids\n",
    "event_ids = np.ones(len(event_times))\n",
    "\n",
    "# Combine the event times and ids into a single events array\n",
    "events = np.column_stack([event_times, np.zeros(len(event_times)), event_ids]).astype(int)\n",
    "\n",
    "# Now you can create the epochs\n",
    "epochs_without_medication = mne.Epochs(raw_without_medication, events=events, tmin=0, tmax=epoch_length - 1/raw_without_medication.info['sfreq'], baseline=None)\n",
    "epochs_with_medication = mne.Epochs(raw_with_medication, events=events, tmin=0, tmax=epoch_length - 1/raw_with_medication.info['sfreq'], baseline=None)\n",
    "\n",
    "# Feature extraction\n",
    "\n",
    "# Power Spectral Density (PSD)\n",
    "psds_without_medication, freqs = mne.time_frequency.psd_array_welch(epochs_without_medication.get_data(), sfreq=raw_without_medication.info['sfreq'], fmin=1., fmax=50.)\n",
    "psds_with_medication, freqs = mne.time_frequency.psd_array_welch(epochs_with_medication.get_data(), sfreq=raw_with_medication.info['sfreq'], fmin=1., fmax=50.)\n",
    "\n",
    "\n",
    "# Bandpower\n",
    "delta = (1, 4)\n",
    "theta = (4, 8)\n",
    "alpha = (8, 13)\n",
    "beta = (13, 30)\n",
    "gamma = (30, 50)\n",
    "\n",
    "#bandpower_without_medication = np.array([pyeeg.band_power(epoch, delta, theta, alpha, beta, gamma) for epoch in epochs_without_medication.get_data()])\n",
    "#bandpower_with_medication = np.array([pyeeg.band_power(epoch, delta, theta, alpha, beta, gamma) for epoch in epochs_with_medication.get_data()])\n",
    "\n",
    "# Statistical features\n",
    "mean_without_medication = np.mean(epochs_without_medication.get_data(), axis=2)\n",
    "mean_with_medication = np.mean(epochs_with_medication.get_data(), axis=2)\n",
    "\n",
    "std_without_medication = np.std(epochs_without_medication.get_data(), axis=2)\n",
    "std_with_medication = np.std(epochs_with_medication.get_data(), axis=2)\n",
    "\n",
    "# Non-linear features\n",
    "#hurst_without_medication = np.array([nolds.hurst_rs(epoch) for epoch in epochs_without_medication.get_data()])\n",
    "#hurst_with_medication = np.array([nolds.hurst_rs(epoch) for epoch in epochs_with_medication.get_data()])\n",
    "\n",
    "# Combine the features and labels\n",
    "features = np.concatenate((psds_without_medication, psds_with_medication, mean_without_medication, mean_with_medication, std_without_medication, std_with_medication))\n",
    "labels = np.concatenate(([0] * len(psds_without_medication), [1] * len(psds_with_medication)))\n",
    "\n",
    "# Split the data into training and test sets\n",
    "features_train, features_test, labels_train, labels_test = train_test_split(features, labels, test_size=0.2, random_state=42)\n",
    "\n",
    "# Train a classifier\n",
    "clf = RandomForestClassifier(random_state=42)\n",
    "clf.fit(features_train, labels_train)\n",
    "\n",
    "# Evaluate the classifier\n",
    "labels_pred = clf.predict(features_test)\n",
    "print(classification_report(labels_test, labels_pred))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}