Classifying SSVEP extended covariance matrices with MDM algorithm

covarianceSSVEPclassification.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "from numpy import nan_to_num, array, empty_like, empty, vstack, concatenate, linspace, tile\n",
    "from scipy.signal import filtfilt, butter\n",
    "import matplotlib.pyplot as plt\n",
    "import pickle\n",
    "import gzip"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classifying extended SSVEP covariance matrices for EEG-based BCI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first step is to load the covariance matrices obtained from the first notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with gzip.open('data_sample_covariance.pz', 'rb') as f:\n",
    "    o = pickle.load(f)\n",
    "cov_train = o['cov_train']\n",
    "cov_test = o['cov_test']\n",
    "y_train = o['y_train']\n",
    "y_test = o['y_test']\n",
    "classes = o['classes']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Covariance matrices are in separated in training and testing set, with `y` the associated labels."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing mean covariance for each class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate the mean covariance, this notebook relies on the implementation of [pyRiemann toolbox](https://github.com/alexandrebarachant/pyRiemann). The function mean_riemann takes a numpy ndarray of shape (n_sample, n_dim, n_dim), where n_sample is the number of covariance matrices used to compute the mean and n_dim is the size of the covariance matrices. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyriemann.utils.mean import mean_riemann"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean of each class is computed with:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cov_centers = empty((len(classes), 24, 24))\n",
    "for i, l in enumerate(classes):\n",
    "    cov_centers[i, :, :] = mean_riemann(cov_train[y_train == l, :, :])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting mean covariance matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaUAAAGxCAYAAADClakjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+UXWV97/HP98xkMjPJ5AeQBEg0IEGCIEZ6zVVD4vij\nEi3ecL0tRe0tVi+3axVar9x6+aGWrIgUSkvV1bJ6l0IaivijWH61FZFiyKRLrrkCAhIkegkISQaR\nAAkJSWbme/84e5LJcJ5nn5nzYz8z5/1aaxZn9nOevR9O9ne+Z+/9/DB3FwAAKSgV3QAAAIaRlAAA\nySApAQCSQVICACSDpAQASAZJCQCQDJLSJGNmrzOzl83Mim4LAIwVSamJzGyrme3JksY2M1trZt01\n7vNJM3vP8O/u/kt3n+EMQANkZheY2SYze9XMbhix/eRs+wtm9mszu9vMTh5RvtbM1oza10IzGzIz\n/m42EB9uc7mk33L3GZKWSHqrpEuLbRIwqT0r6QuSrq+w/Rx3P0LSUZLulPTNKvbHl70GIyk1n0mS\nuz8n6XsqJyeZWYeZ/aWZPWVm283sOjObmpUdaWZ3mtnO7Fvdfdn2GyW9XtKd2dXXn47+NmdmPzCz\nNWa2MXvPXWZ2xMHGmP1+dgX3KzP73OgrL2Aic/fb3P0OSS+M2v6yuz+Z/domaUjSCdXu18yOMbNd\nWUy9bGavmNlg/VreukhKBTGzBZI+IGlLtulqSYsknZb9d76kP8vK/qekX0o6UtJcSZdJkrv/vqSn\nJZ2V3bL7y+z9o7/NfUTSeZLmSJoq6U+zNrxJ0t9m5cdIminp2Hr+fwIpM7OdkvZI+rKkL+a9ffiF\nu293954s7mZIulXSNxrX0tZBUmq+28zsZZWTSb+k1dn28yV92t1fcvdXJF2lcrKQpAMqJ43j3X3Q\n3f991D7zOjWsdfdfuPs+Sd9WdnUm6b9IusPdf+juAzqUBIGW4O6zVf4ydqGkn4wq/kz2zOkFM3uh\nQrkkycwulnSSpE82tLEtgqTUfKuyb1a9khZLOsrM5kjqlvTjEQHwXZWvjCTpGkm/kHS3mf08C4Kx\n2DHi9R5J07PXx6p8BSZJcve9kn49xn0DE1p23v9vSTea2VEjiq5x9yOGf1S+i3EYM/uApD9WOa73\nNafFkxtJqfmGnyltkLRO0l9Kel7lZHHKiCCY5e4zs/fudvc/dfcTJP0nSReZ2buz/dXy4HW7pAUH\nG2bWpUOJEGglbSp/MZxfbQUzO0nSWkm/4+7bGtWwVkNSKtaXJP2mpDdL+qqkL2VXTTKz+Wb2/uz1\nb5nZ8EPYXZIGJA0/VO2X9IZR+612jNItkj5kZm83syk6dCsRmBTMrM3MOlVOOu1mNjXb9j4zW2Jm\nJTObIelalTtDbM7bZbbfHkm3Sfqsu/+wkf8PrYak1FyHXdW4+/OSbpT0eUkXS/q5pPvN7EVJd0t6\nY/bWEyXdY2a7JP27pL/NrrQk6c8lfT677XdRheMEr6Tc/TGVbz18S9I2SS9Lek4StyEwWXxO5bsQ\nF0v6WPb6s5Jmqdwx4UWVOxsdL2mlu+/P6oXiZnj76SrH519nve92Zc+KUSNjjCWGmdk0lYN0kbs/\nVXR7ALQerpRanJmdZWZdWUL6K0kPk5AAFIWkhFUq37p7RuXBg+cW2xwArYzbdwCAZLTXUtnMVqrc\ng6wk6Xp3v7rCe8h6mPDcvSmzrhNTaCWV4mrcV0rZ3GpPSHqvyrd/Nkk6190fH/U+3/8P/3rw9zX/\ndJP+7MO/Vy47ckr0GEP/b0+wrHTCtHjdLbvDdU/pOez3NevW6c/OO0+S5NtzOp4NhD+vKb+7Il53\nDFZfcaVWf+6yuu2vHurVprxzbiyrboylTeM910vdM5qSlMYSUx/QHx38fYt+pBO1VJJ0YdfR0WN8\nf2+4g9gVnzglWveH3wkPxVlx+ZsO+/0Ld31Tn1956E6w7xkI1i2dMiN63Ckr3xktr9Zkjql6alab\nrKunYlzV8kxpqaQt7v6Uux9QeYbdVTXsD2h1xBRaXi1Jab5GTFGj8oPyqkdDA3gNYgotr6ZnStVa\n8083HXw9qzt+260I73rLW4puwmv0rlhedBNeo1XatH5Dn9Zv6Kv7futpi3508HW7phbYkspWLDq1\n6Ca8Rqucv7VqVJuqjatanim9XdJqd1+Z/X6JJB/9YHb0M6XDyhJ5pjRSKs+UJrN6PlOq53FDmvhM\nqeqYGvlMaaRUnimNlsIzJaQl9EypliulTZIWmdlClSf2PFeHllo4/OCB5OM7D0QPUDotfLL68/uD\nZZKkzrZw2d6hYJEd1xXf7yAdnyaqRiW7Oqo6pkLJZ+nK+HJYr398ZrBsyulHBMskaemB8Llvx3VH\n68Y++tK8zmhdtJZxJyV3HzSzC1Weo224+2reZIYAAogpoMZnSu5+l8qLWwGoA2IKrY5phgAAySAp\nAQCSQVICACSDpAQASAZJCQCQjKbM6BAaBBsbhyRJ333ffcGyD9x5RrTu3lt/FSyb9j/CA/22ffon\n0f0+vy08oPetj7wtWhfFmkzLtIQGwcbGIUnShsf7g2Unbp8Xrbvj4fC5f/xvhgfHSpLvCpd7V2RM\nIVoOV0oAgGSQlAAAySApAQCSQVICACSDpAQASAZJCQCQjKZ0CQ+tfZS3/ESs2/fAv26P1h14Ndz9\nd+jxF4Jlcz90VHS/c54Jd4tF2ibA0hVVC619lLf8RKzbd9sZ8XP/uBOmB8tKi8Jl5Z1H/tRMy6mL\nlsKVEgAgGSQlAEAySEoAgGSQlAAAySApAQCSQVICACSDpAQASEZzlq7YsrtyQWd8yvrY8hOxcUiS\ndP0tPw+Wfeo9c4JlW254JrrfTY8/Fyz7g6uiVYG6+eF3tlXcvvRAPC5iy0/ExiFJ0v4N4XjsnDs1\nWtcH9wXLSvOGonXRWrhSAgAkg6QEAEgGSQkAkAySEgAgGSQlAEAySEoAgGSYe7wLac0HMPMD//Zv\nlQv3xruC2jE9wbLY8hOS5C8PhAvbwksYlE6svMzGQfvDbW5f8R/jdSFJyjvnUltiwrp65O7JNMrM\nfN+1t1YuO647Xnl3OC7ylp/wF8JLzdj8WfHjtkdGn3R2Rau2LZgf3zcmpFBc1TROycy2SnpJ0pCk\nA+6+tJb9Aa2OmEKrq3Xw7JCkXnffWY/GACCm0NpqfaZkddgHgEOIKbS0Wk9+l/R9M9tkZufXo0FA\niyOm0NJqvX23zN23m9kclQNps7tvHP2mNevWHXz9rre8Re9asqTGwwKNs35Dn9Zv6Cvq8FXF1Bfu\n+ubB1ysWnap3LTq1mW0ExqzauKpb7zszu1zSLne/dtR2et/hMPS+q/K4kZii9x0mulBcjfv2nZl1\nm9n07PU0Se+X9Oj4mwi0NmIKqO323TxJt5qZZ/v5urvfXemNvr3ytPV2XPwb0rZP/yRYNvdDR0Xr\nxpagWHzdm4NlX3zP96L7PUrh5TYu2NNaV0qxK55GXe008iqr0WP2qlB9TO2pfMWT97/vu2J3EOJ/\nDmLLT1jsSkiS9ofrqiO+7AVay7iTkrs/KYmHQ0CdEFMAXU8BAAkhKQEAkkFSAgAkg6QEAEgGSQkA\nkIxaZ3SozkCgq+1gvAvu89v2BMvmPBMuk6RNjz8XLFscGQAb6/ItSf/iLwbLLojWBOqndMqMytvn\ndUbreVfk/J4WHzxbmhcZ7J4zADbW7ds6421Ga+FKCQCQDJISACAZJCUAQDJISgCAZJCUAADJICkB\nAJJBUgIAJKNui/wFD2DmvndXQ4+Riu/N+VK0fNlvLwiWbdkYHv8kSW/6yDHBMj8QHj9SmhMfA7Lj\n1vB4ru458WFsR938u9HyyaKoRf5CWimmJMkHIstt5CliwchS/Lt+aotYFqXui/wBAFBvJCUAQDJI\nSgCAZJCUAADJICkBAJJBUgIAJKM5S1e0iFiXb0lq7wx/Bzj5w/OidW1eeGkA6wjv13+9L7rfuUu6\ng2WlxTOjdYHC5Q1pSbD7dWwYDt3FuVICACSEpAQASAZJCQCQDJISACAZJCUAQDJISgCAZJCUAADJ\nyB2nZGbXSzpLUr+7n5Ztmy3pW5IWStoq6Rx3f6mB7ZwQ8pafiI1FuuOazdG6q75wWrDMZkwJlpUW\nhMchSdITN20Lls18dG+07oJPhstqGYvR6OVUikZMjVHsfMkb15PguB/GIsVVc6W0VtKZo7ZdIuke\ndz9J0r2SLq13w4BJjJgCAnKTkrtvlLRz1OZVktZlr9dJOrvO7QImLWIKCBvvM6W57t4vSe6+Q9Lc\n+jUJaEnEFKD6zX0XfQiw+oorD77uXbFcvSuW1+mwQP2t39Cn9Rv6im4GMYVJpdq4smoeKpvZQkl3\njngou1lSr7v3m9nRkn7g7icH6rrv3TWmxk9UD775a9HyIjo62OyO6H6fuOrnwbKZc6ZG6y747keD\nZZOpo0Ope4bcva5Pp4mp6vng4PgrF9GpIOeYdHQos66einFV7e07y36G3SHp49nr8yTdXlPrgNZD\nTAEV5F4pmdnNknolHSmpX9Llkm6T9I+SXifpKZW7r1bsD91K3+r2Xfkv0fLY8hP+4v543a62cNmR\nkSuanK8dg4+Fex2XFs2I1u342LviO58kQt/oxr0/YmpMGnalNDQ0vnp5SvGg40qpLBRXuc+U3D10\nj+Z9NbcKaEHEFBDGjA4AgGSQlAAAySApAQCSQVICACSDpAQASAZJCQCQjHpNMwRJfiAy7kGSdYS/\nA8RmZcgr3/tv/cGy7t9ZEN2vBiIzL7QzngKTWGy8UN6sIjWMNapltpNWwJUSACAZJCUAQDJISgCA\nZJCUAADJICkBAJJBUgIAJIMu4XVUmtMZLfdf7wvXXdAd33mkG2ms2/fAv26P7ra0OLI8RSffWdCi\nGtg1m27fcfzVAQAkg6QEAEgGSQkAkAySEgAgGSQlAEAySEoAgGTQJbyOdtz6XLR87pJwt+8nbtoW\nrXvCh44KF0Zm+o52+Zb09U//32DZitPnx9u06oxoeRFqmYE5VhcFKjXou3OjumY3cIbxVsCVEgAg\nGSQlAEAySEoAgGSQlAAAySApAQCSQVICACSDpAQASEbuOCUzu17SWZL63f20bNvlks6XNDww5zJ3\nv6thrZwguufEP87S4pnBspmP7o3XXRQeb2TtkXEPOctPxMYildom3niKWpYFaNaSAsTU2MT+XfLG\nlo3337Sm/dYwHo5lLaq7Ulor6cwK269199OzH4IHqB4xBQTkJiV33yhpZ4UiUjowDsQUEFbLM6UL\nzewhM/uamYXvSwGoFjGFljfeue+uk7TG3d3MrpB0raRPht68+oorD77uXbFcvSuWj/OwQOOt39Cn\n9Rv6mn1YYgqTWrVxZdVMQmlmCyXdOfxQttqyrNx97678Fk8Cz3/0W9HyGe+dEyx77pYd0bpzfy/c\nIaGWjg5PX/NksCyvo8Px9/1etHyysK4euXtdb60RU/WRZEeHGvbdSh0dQnFV7e0704j73WZ29Iiy\nD0t6tLbmAS2HmAIqyL1SMrObJfVKOlJSv6TLJb1b0hJJQ5K2SvpDd+8P1J9Q3+pqWb6glb7lSNK+\nvw53ECudHO7C7r+Md3+3Y7vCZT1t8UZ1RL5n7R+K1+2uvO8py99d1yulVospVM8HBnLeEPn71Kgl\nPvKOO86/e6XpsyrGVe4zJXf/aIXNa8fVCgDEFBDBjA4AgGSQlAAAySApAQCSQVICACSDpAQASAZJ\nCQCQjPFOMzRpFTXWqJEjyBslNhZJL4XHW5RO7onveF9kPFFXzjilFw8Ei+x1s+N1Ozvj5UCj5Y2T\nHP8wyvx9J4IrJQBAMkhKAIBkkJQAAMkgKQEAkkFSAgAkg6QEAEhGU7qE17IcxERSS7ftFLt854kt\nQRHr9r3vG09H99vx7nnhwpzlJ0qvD3dTH9r8fLSuze2IlgMN18jlJybI3xiulAAAySApAQCSQVIC\nACSDpAQASAZJCQCQDJISACAZJCUAQDKaMk5pIo7BabaJuHSFHdsVLowsPxEdhyTJX4osPzEjfsr6\nK+GxU3bklGhd62YlFyQu9nciwb8R48GVEgAgGSQlAEAySEoAgGSQlAAAySApAQCSQVICACQjtw+s\nmS2QdKOkeZKGJH3V3b9iZrMlfUvSQklbJZ3j7i81sK1NUcsyGyl2224k62kLF3ZFynKWn4h2+941\nEK/7upnBMrc90bpqj3cZr5dWiylgLKq5UhqQdJG7nyLpHZIuMLPFki6RdI+7nyTpXkmXNq6ZwKRC\nTAEBuUnJ3Xe4+0PZ692SNktaIGmVpHXZ29ZJOrtRjQQmE2IKCBvTMyUzO07SEkn3S5rn7v1SOcgk\nza1344DJjpgCDlf1vCpmNl3SLZI+5e67zWz0w5fgw5jVV1x58HXviuXqXbF8rO0Emmb9jx/QfT9+\noOHHIabQStZv6NP6vo2577NqHuybWbukf5b0XXf/crZts6Red+83s6Ml/cDdT65Q133vrrG2vzBF\ndXSYiHPfDWz4P+HCSEcHf35/fMdtkf/XnI4OpTcdGT7u7nhHB+ucWnF7+9J3yt3r+g/QSjGF6vng\nYM4bJs/cd6XpsyrGVbW3726Q9Nhw8GTukPTx7PV5km6vqYVAayGmgAqq6RK+TNLHJD1iZg+qfEvh\nMklXS/q2mX1C0lOSzmlkQ4HJgpgCwqq6fVfTAbjVUJWibt/Fjpt3zIH7N4ULI7fZ7NieeJtiy0/0\nTIvWHfrFzmBZaclx0boqVb5x0H7CG+t++64WxNTk5QPx29NRjbx914DbhrXevgMAoOFISgCAZJCU\nAADJICkBAJJBUgIAJIOkBABIRtXTDNUi1O04r8txLd2VUZ2aPsfIEhT2utnBsqHNz0d3a0eGl5DI\nW34i1u174JZHo3XbVy6MlgMNl+rftVi7huJL0YSGWgTfPqZ3AwDQQCQlAEAySEoAgGSQlAAAySAp\nAQCSQVICACSDpAQASEZTxilNJEWtPFuUmsaCdYdXl1VnZ3i/czuiu7XuyGnZHh7DJCk6JiJvHJK/\nzHIQwJiNcRxS7u7qujcAAGpAUgIAJKPpSWn9hr5mHzIXbapOkm368QNFN6FwSf670Kaq0KbXIimJ\nNlUrxTbdR1JK8t+FNlUnyTb1bSz0+Ny+AwAkg6QEAEiG1dIFuqoDmDX2AEATuHsy/f2JKUwWleKq\n4UkJAIBqcfsOAJAMkhIAIBkkJQBAMpqWlMxspZk9bmZPmNnFzTpujJltNbOfmNmDZvajgtpwvZn1\nm9nDI7bNNrO7zexnZvY9M5uZQJsuN7NnzOyB7Gdlk9u0wMzuNbOfmtkjZvYn2fZCP6siEVPBNhBT\n1bUpyZhqSlIys5Kkv5F0pqRTJH3EzBY349g5hiT1uvtb3X1pQW1Yq/LnMtIlku5x95Mk3Svp0gTa\nJEnXuvvp2c9dTW7TgKSL3P0USe+QdEF2DhX9WRWCmIoipqqTZEw160ppqaQt7v6Uux+Q9E1Jq5p0\n7BhTwbcw3X2jpJ2jNq+StC57vU7S2Qm0SSp/XoVw9x3u/lD2erekzZIWqODPqkDEVAAxVZ1UY6pZ\nJ898Sb8c8fsz2baiuaTvm9kmMzu/6MaMMNfd+6XyiSNpbsHtGXahmT1kZl8r8jaZmR0naYmk+yXN\nS/SzajRiamyIqYiUYqrVOzosc/fTJX1Q5UvXM4puUEAKg8muk/QGd18iaYeka4tohJlNl3SLpE9l\n3+5GfzYpfFatjJiqHjFVQbOS0rOSXj/i9wXZtkK5+/bsv7+SdKvKt0RS0G9m8yTJzI6W9FzB7ZG7\n/8oPjbT+qqS3NbsNZtaucvD8g7vfnm1O7rNqEmJqbJI7T4ipypqVlDZJWmRmC82sQ9K5ku5o0rEr\nMrPu7BuCzGyapPdLerSo5ujwe8t3SPp49vo8SbePrtAEh7UpOzmHfVjFfFY3SHrM3b88YlsKn1UR\niKmc5oiYqkZ6MeXuTfmRtFLSzyRtkXRJs44bac/xkh6S9KCkR4pqk6SbJW2TtE/S05L+QNJsSfdk\nn9fdkmYl0KYbJT2cfWa3qXzfuZltWiZpcMS/2QPZOXVEkZ9VkT/EVLAdxFR1bUoyppj7DgCQjFbv\n6AAASAhJCQCQDJISACAZJCUAQDJISgCAZJCUAADJICkBAJJBUgIAJIOklAgz68hmCt5qZi+NXPTL\nzKaY2T+a2ZNmNmRmK0bVXWtma0ZtW5i9l39jtCwzuyCbsfxVM7thxPaPmtkuM3s5+3kli5e3ZuXE\nVEH4cNPRrvL0I8vdfaakz0v6tpkNT7rZJ+ljkraPYZ9M14FW96ykL0i6fuRGd7/Z3XvcfYa7z5D0\nR5J+4e4P5uyPmGqw9qIbgDJ33yNpzYjf/8XMnpT0G+5+q6SvSJKZDY1132Z2jKQndCig2iR1untb\nzQ0HEubut0mSmb1N8fWmzlN5LrqqEFONQ1JKVDZ1/ImSfjreXQy/8PJyAj0j9n1Tba0DJg8zWyhp\nucqTpEbfOvyCmGocklKCsjVObpL09+7+RJXVPmNmF474veI3NjO7WNJJklJdfA1ott+X1OfuT43a\nTkwVgGdKiTEzUzkh7ZP0x2Ooeo27HzH8I+m0Cvv+QLbPVe6+ry4NBia+/yrp7ytsJ6YKwJVSeq6X\ndJSkD7r7YL12amYnSVor6T+7+7Z67ReYyMxsmaRjJH1nHHWJqQYgKSXEzP5O0mJJ73P3/aPKOnTo\nynaqmU2t4puZZXV7VF5E7LPu/sM6NxtIlpm1SZqi8q23djObKmlgxBe+8yR9x91fqXaX2X6JqQbh\n9l0isq7f/13SEkn9I8ZQfCR7y88kvSLpWEl3Sdozort4yHDPoNMlvVHSX2f73GVmL9f//wJIzuck\n7ZF0scpDKvZI+qwkZQnqt1X51l2o6zcx1WCsPAsASAZXSgCAZJCUAADJICkBAJJRU++7bMLQL6mc\n3K5396srvIeHVpjw3N3y31U7YgqtpFJcjbujQzZT7hOS3itpm6RNks5198dHvc+Hfv3cwd9XX/0X\nWn3x/yo3aCg+DMc6poYLDxyIN7A9MgXVwOHHXX31NVp98WfKv3RMie/Xwn+brL1+PexXX3GlVn/u\nsrrtrx5atU3W1dOUpDSmmHrx+YO/r/7zq7X60ourO0gs3NtybpzE6o76dFZfeZVWX3bJiLpjqDy6\nNC8mq9Sq5+9YNatNobiq5fbdUklb3P0pdz8g6ZuSVtWwP6DVEVNoebUkpfmSfjni92cUn4UXQBwx\nhZbXlBkdVl/9Fwdfz5oxsxmHHJPeZe8sugmv0btiedFNeI1WadP6DX1av6Gv7vutp9V/fuhR06yZ\nMwpsSWW9y9Obm7RVzt9aNapN1cZVLc+U3i5ptbsPr456iSQf/WB29DOlkVJ5pnSYRJ4pIR1NfKZU\nfUyNeKY0Jk16pvTausU/U0JaQnFVy1/RTZIWZWuRbJd0rqSPVHpjMPnsz0ksHR3BIs9JSlYKn+i+\nLzxlnE3J+UgiSQmoUdUxFZT3HbNRp2/ud1viBtUZd1Jy98FsrZG7daj76ua6tQxoMcQUUOMzJXe/\nS+XFrQDUATGFVseMDgCAZJCUAADJICkBAJJBUgIAJIOkBABIRlNGewYHwUbGIUmShsKDH6y7K143\nMljPpk+L1x3nfoGmCZ2GecOBhobCZaUGfkf1yHGBEbhSAgAkg6QEAEgGSQkAkAySEgAgGSQlAEAy\nSEoAgGQ0ZwGgwDITuctPRLp9+67d8bqd4bWY/NXI0hXTuqP71RTWdkEC8tY+Col1+65ltENeV/Q2\n1hpDdbhSAgAkg6QEAEgGSQkAkAySEgAgGSQlAEAySEoAgGSQlAAAyWjO4IH2toqbrZQzuCG2/ERk\nHJIk+XiXvShVbms1bQKapojTMLb8hOV8v40tmWF5g5zQSrhSAgAkg6QEAEgGSQkAkAySEgAgGSQl\nAEAySEoAgGQ0p0v4wGDFzb4vvISEJNn0acGy2PIT0viXvbCZM6L7lTEFPxIQ6kXdyOUnYt2+84ZK\nxJbMAEao6S+smW2V9JKkIUkH3H1pPRoFtCpiCq2u1q/9Q5J63X1nPRoDgJhCa6v1mtrqsA8AhxBT\naGm1nvwu6ftmtsnMzq9Hg4AWR0yhpdV6+26Zu283szkqB9Jmd984+k2rr77m4OveZe9U7xnLajws\n0DjrN/Rp/Ya+og5fXUxdedXB173Lz1Dv8jOa2UZgzKqNK/M6TTBqZpdL2uXu147a7kPP91esU1Pv\nu0gPOqmBve/aw3nc2nImc8WEZF09cvemzxoajamXX6hcqZG972Ia2PvOIjGHiSsUV+M+U8ys28ym\nZ6+nSXq/pEfH30SgtRFTQG237+ZJutXMPNvP19397orv7JhScbNNGf/hbVp3/A2RJShiV0O+f3/8\nuLFvfFwpoTbVx1TwyiRvOZjIEhJtOfEYW34i70ooMFZRUv4VGldKLWXc/9ru/qSkJXVsC9DSiCmA\nrqcAgISQlAAAySApAQCSQVICACSDpAQASEZz+lpaoM9naPuw2IC8KZW7mVdVN7L8RLTLtyQNRrrF\n5jQJqJ+mj+XNj9do3UhZfcbvY5LgSgkAkAySEgAgGSQlAEAySEoAgGSQlAAAySApAQCSQVICACSj\nKeOUJtQiXXnLT0TGIg3tfDFn35HvALFlASTZlI5gmQ9FlgWoYTxL3qKFsYUU0VgWWA4mWTX8Dfhg\n18XBsjMsvihn7Fv3WYuPCZa9+upAdL/HntATLJt70fHRujZnerCs7dRTo3VbAVdKAIBkkJQAAMkg\nKQEAkkFSAgAkg6QEAEgGSQkAkIwJ1Fd7Aoh1+Zbi3b7b4v8U7rEu45Fu39F6knVMDReWClgeARgl\n1u17SWd3tG5HR3hYw5w3hOv6YHw9jZ754SEamp7zZzVv2EmL40oJAJAMkhIAIBkkJQBAMkhKAIBk\nkJQAAMkgKQEAkkFSAgAkI3eckpldL+ksSf3uflq2bbakb0laKGmrpHPc/aUGtnNiyFl+IjoW6cCB\naFWbGh5PFBvDZHnjn/btCxfmjLuyjshYDQQRU2MTOwtj45AkyT083qitIzwOb+DVeJvapsbGBsbH\nOCGumiuwl1XPAAAL1klEQVSltZLOHLXtEkn3uPtJku6VdGm9GwZMYsQUEJCblNx9o6SdozavkrQu\ne71O0tl1bhcwaRFTQNh4nynNdfd+SXL3HZLm1q9JQEsipgDVb+676E3U1VdcefB174rl6l2xvE6H\nBepv/YY+rd/QV3QziClMKtXGlcUeBB58k9lCSXeOeCi7WVKvu/eb2dGSfuDuJwfquu/dNabGT1RD\nL74Yf0OjOjpEOlhYKX4x7AMD4cKcjg6l6dOj5ZOFdfXI3es6Oy0xVb2rur8YLPsPM3qidWN/3976\n3qODZQOvxv8uzjo+3Mmn/XdfF61rPeG4aTu54j/5pBSKq2pv35kOn4r6Dkkfz16fJ+n2mloHtB5i\nCqigmi7hN0vqlXSkmT0t6XJJV0n6RzP7hKSnJJ3TyEZOFDYl3kU62nU7ciUkST40GK4buRqKXUVJ\nknV1RctRf8TU2Jy1+JhgWWz5CSne7Xv6SZG6Ocu2lI7qDJbZ9HibFFsuBvlJyd0/Gih6X53bArQE\nYgoIY0YHAEAySEoAgGSQlAAAySApAQCSQVICACSDpAQASEa9phmC4mOJysJjH2JjmKScsUgDkePm\nzMqQu9wGULBXXw3POuKD8ZkXoktQxMYi7c+Ji1hYDbF0RS24UgIAJIOkBABIBkkJAJAMkhIAIBkk\nJQBAMkhKAIBk0CW8rnLWgYstXRFbAFA5S1DEun3HFvGTpNiSGVbXde2AcTn2hPBCfj3z48vFtE0N\nn8Ox5Sdyv67PmBIuy1mGRlMidcGVEgAgHSQlAEAySEoAgGSQlAAAySApAQCSQVICACSDLuFNZB3h\nrqK+b1+8bldXuDDWXTyne6q/+FK4sDNe13rCXXWBepl70fHhwuk5f8I8PGO3Te8O18ub6TsSVzZz\nVrxuW1u8vMVxpQQASAZJCQCQDJISACAZJCUAQDJISgCAZJCUAADJICkBAJKRO07JzK6XdJakfnc/\nLdt2uaTzJT2Xve0yd7+rYa2cICxv/EEpshREbPmJWuQtP5EzFgn1R0yNjc2ZHi6sZcxPZNxgrtjy\nE3ltio0rRFVXSmslnVlh+7Xufnr2Q/AA1SOmgIDcpOTuGyXtrFDECnDAOBBTQFgt94wuNLOHzOxr\nZjazbi0CWhcxhZY33rnvrpO0xt3dzK6QdK2kT4bevPqKKw++7l2xXL0rlo/zsEDjrd/Qp/Ub+pp9\nWGIKk1q1cWUembDw4JvMFkq6c/ihbLVlWbn73l35LZ4EfM/e+Bvaww9Aff/+aNXYZK7RB6eRY0qS\n781pc0SpRSZkta4euXtdb60RU9UbfPTRcGGCHR1s2rR43Ui8luYcNd4WTTihuKr29p1pxP1uMzt6\nRNmHJUXOGgAVEFNABdV0Cb9ZUq+kI83saUmXS3q3mS2RNCRpq6Q/bGAbJwzrjiwvkVe3o6OOLRnD\ncWu42jmn67Jg2ZrTTgqWffbhx6P7/cZfvSdY9uy3+6N1X/+VNwfLdn/liWjdzkU533DrhJgam7ZT\nTy26CU3jr8aXsIl2hcm76RW7a5LXTT02tKSWuhXkJiV3/2iFzWvHdBQABxFTQBgzOgAAkkFSAgAk\ng6QEAEgGSQkAkAySEgAgGSQlAEAyxjvNEBAdi3T8qvDI9DWDb4zut3RieKmCY8/OGRMxtTNY1HVG\nfLS8zY4sRwA0Qy3zhtQ050gNlXPHIY1t31wpAQCSQVICACSDpAQASAZJCQCQDJISACAZJCUAQDLo\nEo5xiy1BEev2/ZmfxpeuuHXz/GDZM7c8F6173LvmBMtevjted/qJ3dFyoOHylp+ILl2Rv2DruA9s\nkeuX3KUrxtYSrpQAAMkgKQEAkkFSAgAkg6QEAEgGSQkAkAySEgAgGSQlAEAyzGvq217FAczc9+5q\n6DFQjP3XfT9YFlt+YnBz/HxoP+PIYJm/eCBa146bG667NT5OSe2VB1RMOfN9cveaFgaoJ2Jq8vKB\ngcbtPLbERF4eyF2eYuxK02dVjCuulAAAySApAQCSQVICACSDpAQASAZJCQCQDJISACAZuUtXmNkC\nSTdKmidpSNJX3f0rZjZb0rckLZS0VdI57v5SA9uKxDz77f5g2bFnh6ezz1t+YuGpM4Jlgw/HT7H2\nBUcEy4ae2RutazOnRMvrhZhCUN4yEPG1K+JVSzUsP1FL3TGuXVHNldKApIvc/RRJ75B0gZktlnSJ\npHvc/SRJ90q6dExHBloXMQUE5CYld9/h7g9lr3dL2ixpgaRVktZlb1sn6exGNRKYTIgpIGxMz5TM\n7DhJSyTdL2meu/dL5SCTFB5KD6AiYgo4XNXLoZvZdEm3SPqUu+82s9E3MIM3NFdfceXB170rlqt3\nxfKxthNomvt++rDue+zhhh+HmEIrWd+3Ues3bsx9X1Vz35lZu6R/lvRdd/9ytm2zpF537zezoyX9\nwN1PrlCXebomqSd7bwqWHXv2nGBZbkeHNScGy3I7Onzw+HDd+5+N1g11dOg49wN1n/uOmEIlvn9/\nzjsmT0eH0swjapr77gZJjw0HT+YOSR/PXp8n6fYq9wWAmAIqqqZL+DJJH5P0iJk9qHI6vkzS1ZK+\nbWafkPSUpHMa2VBgsiCmgDCWrsC4DT78k3Dh1M5w2Z5X4jueOStclje1f1d3uOylF+N1S5VvM7T/\nxttZugJN4QfiS7NE5S0vUcvSFQ1Q6pnN0hUAgLSRlAAAySApAQCSQVICACSDpAQASAZJCQCQjKqn\nGQJG2/2VJ4JlXWccFSx7+e74jA6z/tvCYFne8hNtp4dnkhi4d1u0bmlupBs70Ax5syPEunXn1W2P\n/LmvZUaHgcGcuvVfugIAgKYgKQEAkkFSAgAkg6QEAEgGSQkAkAySEgAgGSQlAEAyGKeEcetcNC1Y\nZrMrr+IqSdNPjCwvIUnt4XENodVhD4qMicgdh5TM4hRoWXnLT8RO0prO3xoqj3EcUu7u6ro3AABq\nQFICACSj6Ulp/Ya+Zh8yF22qToptuu+nDxfdhMKl+O9Cm6qTZJv6NhZ6fJKSaFO1UmzTfY+RlFL8\nd6FN1UmyTRtbLCkBABBCUgIAJMPcvbEHMGvsAYAmcPdkOowTU5gsKsVVw5MSAADV4vYdACAZJCUA\nQDJISgCAZDQtKZnZSjN73MyeMLOLm3XcGDPbamY/MbMHzexHBbXhejPrN7OHR2ybbWZ3m9nPzOx7\nZjYzgTZdbmbPmNkD2c/KJrdpgZnda2Y/NbNHzOxPsu2FflZFIqaCbSCmqmtTkjHVlKRkZiVJfyPp\nTEmnSPqImS1uxrFzDEnqdfe3uvvSgtqwVuXPZaRLJN3j7idJulfSpQm0SZKudffTs5+7mtymAUkX\nufspkt4h6YLsHCr6syoEMRVFTFUnyZhq1pXSUklb3P0pdz8g6ZuSVjXp2DGmgm9huvtGSTtHbV4l\naV32ep2ksxNok1TgPNruvsPdH8pe75a0WdICFfxZFYiYCiCmqpNqTDXr5Jkv6Zcjfn8m21Y0l/R9\nM9tkZucX3ZgR5rp7v1Q+cSTNLbg9wy40s4fM7GtF3iYzs+MkLZF0v6R5iX5WjUZMjQ0xFZFSTLV6\nR4dl7n66pA+qfOl6RtENCkhhMNl1kt7g7ksk7ZB0bRGNMLPpkm6R9Kns293ozyaFz6qVEVPVI6Yq\naFZSelbS60f8viDbVih3357991eSblX5lkgK+s1sniSZ2dGSniu4PXL3X/mhkdZflfS2ZrfBzNpV\nDp5/cPfbs83JfVZNQkyNTXLnCTFVWbOS0iZJi8xsoZl1SDpX0h1NOnZFZtadfUOQmU2T9H5JjxbV\nHB1+b/kOSR/PXp8n6fbRFZrgsDZlJ+ewD6uYz+oGSY+5+5dHbEvhsyoCMZXTHBFT1Ugvpty9KT+S\nVkr6maQtki5p1nEj7Tle0kOSHpT0SFFtknSzpG2S9kl6WtIfSJot6Z7s87pb0qwE2nSjpIezz+w2\nle87N7NNyyQNjvg3eyA7p44o8rMq8oeYCraDmKquTUnGFHPfAQCS0eodHQAACSEpAQCSQVICACSD\npAQASAZJCQCQDJISACAZJCUAQDL+P7Uc+rcP/sMLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10d6cb690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 7))\n",
    "for i, l in enumerate(classes):\n",
    "    plt.subplot(2, 2, i+1)\n",
    "    plt.imshow(cov_centers[i, :, :], cmap=plt.get_cmap('RdPu'), interpolation='nearest')\n",
    "    _ = plt.title(l)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As seen in the previous notebook, the mean covariance matrices for each class concentrate the highest values in the block corresponding to the filtered signal in the associated bandwith. The resting signal shows high correlations split across all frequencies."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Minimum distance to mean"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The classification scheme is simple, each test sample $\\hat{\\Sigma}$ is associated to the class with the smallest distance to mean $\\Sigma^{c}_{\\mu}$ :\n",
    "\\begin{equation}\n",
    "c^{*} = \\mathrm{argmin}_{c} \\delta(\\hat{\\Sigma}, \\Sigma^{c}_{\\mu})\n",
    "\\end{equation}\n",
    "The Riemannian distance used here is the AIRM, implemented in pyRiemann"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from pyriemann.utils.distance import distance_riemann"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Evaluation accuracy on test set is 91.67%\n"
     ]
    }
   ],
   "source": [
    "accuracy = list()\n",
    "for sample, true_label in zip(cov_test, y_test):\n",
    "    dist = [distance_riemann(sample, cov_centers[m]) for m in range(len(classes))]\n",
    "    if classes[array(dist).argmin()] == true_label:\n",
    "        accuracy.append(1)\n",
    "    else: accuracy.append(0)\n",
    "test_accuracy = 100.*array(accuracy).sum()/len(y_test)\n",
    "            \n",
    "print ('Evaluation accuracy on test set is %.2f%%' % test_accuracy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The evaluation score is high, but it is computed on only one session, so there is no source of variation from the electronic setup. Also, the subject selected for this notebook is the one with the highest result and has very constant performance across session. [Other subjects](https://github.com/sylvchev/dataset-ssvep-exoskeleton/) have lower score and offer more challenge for machine learning techniques."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "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.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}