analysis-n-of-2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "'''\n",
    "t_res: number of raw values to include in each power spectrum\n",
    "bns: number of bins for each power spectrum\n",
    "cmprs: number of consecutive power spectra to average\n",
    "'''\n",
    "t_res = 250\n",
    "bns   = t_res//4\n",
    "cmprs = 1\n",
    "\n",
    "'''\n",
    "Channel 0 (N1P) - Right Helix\n",
    "Channel 1 (N2P) - Right ear canal, front facing\n",
    "Channel 2 (N3P) - Right ear canal, back facing\n",
    "Channel 3 (N4P) - Left Helix\n",
    "Channel 4 (N5P) - Left ear canal, front facing\n",
    "Channel 5 (N6P) - Left ear canal, back facing\n",
    "Channel 6 (N7P) - Approximately FP1\n",
    "Channel 7 (N8P) - Right mastoid\n",
    "'''\n",
    "#right_ear_nomastoid = [0,1,2]\n",
    "right_ear = [0,1,2,7]\n",
    "left_ear = [3,4,5]\n",
    "fp1 = [6]\n",
    "#config for which chans to use\n",
    "channels = fp1\n",
    "outFile = 'fp1'\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "files = os.listdir('data')\n",
    "\n",
    "def subject_id (filename):\n",
    "    return filename.split('.')[0]\n",
    "\n",
    "def task_name (filename):\n",
    "    return filename.split('.')[1]\n",
    "\n",
    "subject_ids = list(set(map(subject_id, files)))\n",
    "task_names = list(set(map(task_name, files)))\n",
    "\n",
    "\n",
    "def find_with_id (id):\n",
    "    return lambda filename: subject_id(filename) == id\n",
    "\n",
    "def find_with_task (task):\n",
    "    return lambda filename: task_name(filename) == task\n",
    "\n",
    "# both args are optional\n",
    "def find_recordings (id=None, task=None):\n",
    "    lst = files\n",
    "    if id:\n",
    "        lst = list(filter(find_with_id(id), lst))\n",
    "    if task:\n",
    "        lst = list(filter(find_with_task(task), lst))\n",
    "    return lst\n",
    "\n",
    "%pylab inline\n",
    "import pandas as pd\n",
    "from os.path import join\n",
    "\n",
    "def load_recording (filename):\n",
    "    return pd.read_csv(join('data', filename), header=None)\n",
    "\n",
    "#load_recording(files[3])\n",
    "\n",
    "\n",
    "# find_recordings(task='breathe_o')\n",
    "# find_recordings('002', 'breathe_o')\n",
    "'''utilities'''\n",
    "def withold_final (percent, lst):\n",
    "    n = int(len(lst)*percent)\n",
    "    return lst[:-1*n, :]\n",
    "#withold_final(0.4,\n",
    "#              np.array([[1,0],[2,0],[3,0],[4,0],[5,0]]))\n",
    "def get_final (percent, lst):\n",
    "    n = int(len(lst)*percent)\n",
    "    return lst[-1*n:, :]\n",
    "#get_final(0.4,\n",
    "#          np.array([[1,0],[2,0],[3,0],[4,0],[5,0]]))\n",
    "def difference (a, b):\n",
    "    return list(set(b) - set(a))\n",
    "#difference([1,2,3], [1,2,3,4,5])\n",
    "\n",
    "import itertools\n",
    "from scipy.interpolate import interp1d\n",
    "from functools import partial\n",
    "\n",
    "def grouper (n, iterable):\n",
    "    \"grouper(3, 'ABCDEFG', 'x') --> ABC DEF Gxx\"\n",
    "    args = [iter(iterable)] * n\n",
    "    return list(itertools.zip_longest(*args))[:-1]\n",
    "\n",
    "def spectrum (vector):\n",
    "    '''get the power spectrum of a vector of raw EEG data'''\n",
    "    A = np.fft.fft(vector)\n",
    "    ps = np.abs(A)*2.0\n",
    "    ps = ps[:len(ps)//2]\n",
    "    return ps\n",
    "\n",
    "def binned (pspectra, n):\n",
    "    '''compress an array of power spectra into vectors of length n'''\n",
    "    l = len(pspectra)\n",
    "    array = np.zeros([l,n])\n",
    "    for i,ps in enumerate(pspectra):\n",
    "        x = np.arange(1,len(ps)+1)\n",
    "        f = interp1d(x,ps)#/np.sum(ps))\n",
    "        array[i] = f(np.arange(1, n+1))\n",
    "    index = np.argwhere(array[:,0]==-1)\n",
    "    array = np.delete(array,index,0)\n",
    "    return array\n",
    "\n",
    "def feature_vector (bins, readings): # A function we apply to each group of power spectra\n",
    "  '''\n",
    "  Create 100, log10-spaced bins for each power spectrum.\n",
    "  For more on how this particular implementation works, see:\n",
    "      http://people.ischool.berkeley.edu/~chuang/pubs/MMJC15.pdf\n",
    "  '''\n",
    "  bins = binned(list(map(spectrum, readings)),\n",
    "                bins)\n",
    "  return np.log10(np.mean(bins, 0))\n",
    "\n",
    "\n",
    "from functools import reduce\n",
    "\n",
    "def features (raws, bins=50, time_resolution=250, spectra_compression=2 ):\n",
    "    featureVectorF = partial(feature_vector, bins)\n",
    "    spectra = map(spectrum, grouper(time_resolution, raws))\n",
    "    groupedSpectra = grouper(spectra_compression, spectra)\n",
    "    vs = list(map(featureVectorF, groupedSpectra))\n",
    "    return np.array(vs)\n",
    "\n",
    "#breathe = load_recording(\n",
    "#    find_recordings(task='breathe')[0])\n",
    "\n",
    "#features(breathe[1], bins=25, spectra_compression=1).shape\n",
    "\n",
    "def featureify (acc, df, bins=bns, time_res=t_res, spectra_compress=cmprs):\n",
    "    f = lambda df, n: features(df[n],\n",
    "                               bins=bins,\n",
    "                               time_resolution=time_res,\n",
    "                               spectra_compression=spectra_compress)\n",
    "\n",
    "    fs = [f(df, chan) for chan in channels]\n",
    "    return acc + np.concatenate(fs, 1).tolist()\n",
    "\n",
    "def vectors (filenames):\n",
    "    return np.array(reduce(featureify,\n",
    "                  # a list of dataframes of recordings\n",
    "                  map(load_recording, filenames),\n",
    "                  []))\n",
    "\n",
    "def load_feature_vectors(id=None, task=None):\n",
    "    return vectors(find_recordings(id, task))\n",
    "\n",
    "breathe_features = load_feature_vectors(id='001', task='breathe')\n",
    "\n",
    "# np.array(breathe_features).shape\n",
    "\n",
    "# here's a realy simple test of authentication power\n",
    "from sklearn import preprocessing\n",
    "\n",
    "def vectors_labels (list1, list2):\n",
    "    def label (l):\n",
    "        return lambda x: l\n",
    "    X = np.concatenate([list1, list2])\n",
    "    y = np.array(list(map(label(0), list1)) + list(map(label(1), list2)))\n",
    "    return X, y\n",
    "\n",
    "# scaler  = preprocessing.RobustScaler()\n",
    "# max_breathe = load_feature_vectors(id='001', task='breathe')\n",
    "# john_breathe = load_feature_vectors(id='002', task='breathe')\n",
    "# X, y = vectors_labels(max_breathe, john_breathe)\n",
    "# X = scaler.fit_transform(X)\n",
    "# cv = cross_val_svm(X,y,7)\n",
    "# (cv.mean(),cv.std())\n",
    "\n",
    "# here's a fancier version of above using xgboost\n",
    "import xgboost as xgb\n",
    "from xgboost.sklearn import XGBClassifier\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn import metrics\n",
    "#from sklearn.grid_search impo\n",
    "\n",
    "# eval_metrics:\n",
    "# http://xgboost.readthedocs.io/en/latest//parameter.html\n",
    "metrics = ['auc', 'map']\n",
    "\n",
    "def cv (alg,X,y,nfold=7):\n",
    "    xgtrain = xgb.DMatrix(X,y)\n",
    "    param = alg.get_xgb_params()\n",
    "    cvresults = xgb.cv(param,\n",
    "                      xgtrain,\n",
    "                      num_boost_round=alg.get_params()['n_estimators'],\n",
    "                      nfold=nfold,\n",
    "                      metrics=metrics,\n",
    "                      early_stopping_rounds=50)\n",
    "    alg.set_params(n_estimators=cvresults.shape[0])\n",
    "    alg.fit(X,y,eval_metric=metrics)\n",
    "    return alg, cvresults\n",
    "\n",
    "def inspect (alg):\n",
    "    return alg.booster().get_score(importance_type='weight')\n",
    "    #Predict training set:\n",
    "    #dtrain_predictions = alg.predict(X)\n",
    "    #dtrain_predprob = alg.predict_proba(X)[:,1]\n",
    "    ##Print model report:\n",
    "    #print(\"Accuracy : %.4g\" % metrics.accuracy_score(y, dtrain_predictions))\n",
    "    #print(\"AUC Score (Train): %f\" % metrics.roc_auc_score(y, dtrain_predprob))\n",
    "    #feat_imp = pd.Series(features).sort_values(ascending=False)\n",
    "    # return alg.booster().get_fscore()\n",
    "    #feat_imp.plot(kind='bar', title='Feature Importances')\n",
    "    #plt.ylabel('Feature Importance Score')\n",
    "\n",
    "\n",
    "def fresh_xgb ():\n",
    "    return XGBClassifier(\n",
    "         learning_rate =0.1,\n",
    "         n_estimators=1000,\n",
    "         max_depth=5,\n",
    "         min_child_weight=1,\n",
    "         gamma=0,\n",
    "         subsample=0.8,\n",
    "         colsample_bytree=0.8,\n",
    "         objective= 'binary:logistic',\n",
    "         nthread=4,\n",
    "         scale_pos_weight=1,\n",
    "         seed=27) \n",
    "#\n",
    "#X, y = vectors_labels(\n",
    "#    max_breathe,\n",
    "#    john_breathe)\n",
    "#\n",
    "#alg, cvresults = cv(fresh_xgb(), X,y)\n",
    "## see result of last cv round\n",
    "##cvresults[-1:]\n",
    "#inspect(alg)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'FAR': 0.0,\n",
       " 'FRR': 0.0,\n",
       " 'f-scores': {'f1': 1},\n",
       " 'right_person_wrong_task_far': 0.0,\n",
       " 'right_person_wrong_task_frr': 0.0}"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# returns FAR, FRR\n",
    "def far_frr (alg, X, y):\n",
    "    false_accept = 0\n",
    "    false_reject = 0\n",
    "    predicted_labels = alg.predict(X)\n",
    "    for idx, predicted_label in enumerate(predicted_labels):\n",
    "        correct_label = y[idx]\n",
    "        if (predicted_label==correct_label):\n",
    "           continue\n",
    "        elif (predicted_label==0):\n",
    "           false_accept+=1\n",
    "           continue\n",
    "        false_reject+=1\n",
    "    def rate (num):\n",
    "        return num/len(X)\n",
    "    return { 'FAR': rate(false_accept),\n",
    "             'FRR': rate(false_reject),\n",
    "             'f-scores': inspect(alg), }\n",
    "\n",
    "def test_inherence (alg, stats, right_person_wrong_task_vectors):\n",
    "    labels = [0 for i in(right_person_wrong_task_vectors)]\n",
    "    right_person_wrong_task_stats = far_frr(alg, right_person_wrong_task_vectors, labels)\n",
    "    return right_person_wrong_task_stats\n",
    "\n",
    "def test_auth (task, percent_to_withold=0, against_all_tasks=True):\n",
    "\n",
    "    def apply_if_witholding (f, vectors):\n",
    "        if (percent_to_withold):\n",
    "            return f(percent_to_withold, vectors)\n",
    "        return vectors\n",
    "\n",
    "    def trainset (vectors):\n",
    "        return apply_if_witholding(withold_final, vectors)\n",
    "\n",
    "    def testset (vectors):\n",
    "        return apply_if_witholding(get_final, vectors)\n",
    "\n",
    "    # get `task` by this `subject`\n",
    "    # HACK just using subj 1 for now\n",
    "    this_persons_features =  load_feature_vectors('001', task)\n",
    "\n",
    "    # we can test against only this task from others\n",
    "    if not against_all_tasks:\n",
    "        # get `task` NOT from this `subject`\n",
    "        not_this_persons_features =  load_feature_vectors('002', task=task)  \n",
    "\n",
    "    # or against all tasks from others\n",
    "    elif against_all_tasks:\n",
    "        this_persons_recordings = find_recordings('001')\n",
    "        not_this_persons_recordings = difference(this_persons_recordings, files)\n",
    "        not_this_persons_features = vectors(not_this_persons_recordings)\n",
    "\n",
    "\n",
    "    # turn these pos/neg features into vectors, labels X, y\n",
    "    trainX, trainy = vectors_labels(trainset(this_persons_features),\n",
    "                                    trainset(not_this_persons_features))\n",
    "\n",
    "\n",
    "    # fit an alg over 7 cv rounds\n",
    "    alg, cvresults = cv(fresh_xgb(),\n",
    "                        trainX, trainy,\n",
    "                        nfold=7)\n",
    "    #inspect(alg)\n",
    "\n",
    "    # turn these pos/neg features into vectors, labels X, y\n",
    "    testX, testy= vectors_labels(testset(this_persons_features),\n",
    "                                 testset(not_this_persons_features))\n",
    "\n",
    "\n",
    "    # get far, frr\n",
    "    stats = far_frr(alg, testX, testy)\n",
    "\n",
    "    # test right person wrong task\n",
    "    right_person_right_task_recordings = find_recordings('001', task)\n",
    "    right_person_wrong_task_vectors = vectors(difference(right_person_right_task_recordings,\n",
    "                                                         find_recordings('001')))\n",
    "    inherence_stats = test_inherence(alg, stats, right_person_wrong_task_vectors)\n",
    "\n",
    "    # HACK\n",
    "    stats['right_person_wrong_task_far'] = inherence_stats['FAR']\n",
    "    stats['right_person_wrong_task_frr'] = inherence_stats['FRR']\n",
    "    return stats\n",
    "\n",
    "test_auth(task_names[2],\n",
    "          against_all_tasks=True)\n",
    "\n",
    "#alg, cvresults = cv(fresh_xgb(), X,y)\n",
    "##far_frr(alg, X, y)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "          task  FAR  FRR   f-scores  right_person_wrong_task_far  \\\n",
       "0       speech  0.0  0.0  {'f1': 1}                          0.0   \n",
       "1         face  0.0  0.0  {'f1': 1}                          0.0   \n",
       "2      breathe  0.0  0.0  {'f1': 1}                          0.0   \n",
       "3  listennoise  0.0  0.0  {'f1': 1}                          0.0   \n",
       "4       song_o  0.0  0.0  {'f1': 1}                          0.0   \n",
       "5         song  0.0  0.0  {'f1': 1}                          0.0   \n",
       "6     sequence  0.0  0.0  {'f1': 1}                          0.0   \n",
       "7   listentone  0.0  0.0  {'f1': 1}                          0.0   \n",
       "8    breathe_o  0.0  0.0  {'f1': 1}                          0.0   \n",
       "9        sport  0.0  0.0  {'f1': 1}                          0.0   \n",
       "\n",
       "   right_person_wrong_task_frr  \n",
       "0                          0.0  \n",
       "1                          0.0  \n",
       "2                          0.0  \n",
       "3                          0.0  \n",
       "4                          0.0  \n",
       "5                          0.0  \n",
       "6                          0.0  \n",
       "7                          0.0  \n",
       "8                          0.0  \n",
       "9                          0.0  "
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running sport\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running breathe_o\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running listentone\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running sequence\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running song\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running song_o\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running listennoise\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running breathe\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running face\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "running speech\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "when I do \n",
    "\n",
    "    for task in task_names:\n",
    "        print(task, test_auth(task))\n",
    "\n",
    "the results are perfect...........too perfect\n",
    "I think we are overfitting.\n",
    "\n",
    "its confusing, but in the cv round, we are doing 7 rounds of cv to generate some parameters\n",
    "then, we fit those parameters to the algorithm, and return it\n",
    "so, when we test on that same data later, we are effectively testing on the train set,\n",
    "since the paramaters were generated from there.\n",
    "\n",
    "\n",
    "as a simple (dumb) solution, lets withold the last n samples of each group (positive, negative) for testing\n",
    "'''\n",
    "\n",
    "stats = ['FAR','FRR', 'f-scores', 'right_person_wrong_task_far', 'right_person_wrong_task_frr']\n",
    "df = pd.DataFrame(columns=['task'] + stats)\n",
    "\n",
    "for idx, task in enumerate(task_names):\n",
    "    print('running', task)\n",
    "    results = test_auth(task,\n",
    "                        percent_to_withold=0.7)\n",
    "    results_stats = [results[stat] for stat in stats]\n",
    "    row = [task]+results_stats\n",
    "    df.loc[idx] = row\n",
    "    #df.append(row)\n",
    "\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/ffff/.local/lib/python3.5/site-packages/scipy/stats/stats.py:3658: RuntimeWarning: invalid value encountered in absolute\n",
      "  prob = distributions.t.sf(np.abs(t), df) * 2  # use np.abs to get upper tail\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "far diff? nan frr diff? nan\n"
     ]
    }
   ],
   "source": [
    "from scipy.stats import ttest_rel\n",
    "\n",
    "wrong_task_frr = df['right_person_wrong_task_frr'].tolist()\n",
    "right_task_frr = df['FRR'].tolist()\n",
    "pvalue_frr = ttest_rel(right_task_frr, wrong_task_frr)[1]\n",
    "#wrong_task_far = df['right_person_wrong_task_far'].tolist()\n",
    "#right_task_far = df['FAR'].tolist()\n",
    "#pvalue_far = ttest_rel(right_task_far, wrong_task_far)[1]\n",
    "print('far diff?', pvalue_far, 'frr diff?', pvalue_frr)\n",
    "#chisquare(\n",
    "#    df['right_person_wrong_task_frr'].tolist(),\n",
    "#    f_exp=df['FRR'].tolist())\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "df.to_csv(outFile+'-'+str(pvalue_frr)+'.csv')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.2"
  },
  "name": "analysis-n-of-2.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 1
}