cerebro.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Written and Copyright by Mohit Agarwal\n",
    "#### Georgia Institute of Technology\n",
    "#### Email: me.agmohit@gmail.com\n",
    "\n",
    "## Please cite the following publication if you are using the codes for your study and publication\n",
    "# Agarwal, Mohit, and Raghupathy Sivakumar. \n",
    "# \"Cerebro: A Wearable Solution to Detect and Track User Preferences using Brainwaves.\" \n",
    "# The 5th ACM Workshop on Wearable Systems and Applications. ACM, 2019.\n",
    "\n",
    "### Code to decode Product Preference data [10 products] and rank it\n",
    "### Uses fix rule based on N200, Min and ERSP features"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "##Importing Libraries\n",
    "\n",
    "import numpy as np\n",
    "import math\n",
    "import scipy.io\n",
    "import os\n",
    "from sklearn.decomposition import FastICA\n",
    "import mne\n",
    "from mne.time_frequency import psd_multitaper\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "import itertools\n",
    "from scipy import stats\n",
    "\n",
    "from sklearn import svm\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.linear_model import LinearRegression\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.tree import DecisionTreeClassifier\n",
    "from sklearn.linear_model import ElasticNet\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.metrics import confusion_matrix\n",
    "from sklearn.metrics import accuracy_score\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "mne.set_log_level('ERROR')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Parameters\n",
    "\n",
    "chan_list = ['Fp1','F3','F7','C3','T7','P3','P7','O1','Pz','Fp2','Fz','F4','F8','Cz','C4','T8','P4','P8','O2']\n",
    "\n",
    "selected_chan = [1,3,8,10,11,13,14]\n",
    "total_chan = len(selected_chan)\n",
    "\n",
    "total_sub = 14  # Anything from 1-14\n",
    "total_prod = 10 \n",
    "freq = 256.0\n",
    "time_len = 768\n",
    "\n",
    "time = [(x-freq)*1.0/freq for x in xrange(1,time_len+1)]\n",
    "time = np.array(time)\n",
    "\n",
    "n200_ind = [idx for idx, t_ind in enumerate(time) if (t_ind>=0.2 and t_ind<=0.3)]\n",
    "n200_ind = np.array(n200_ind)\n",
    "\n",
    "erp_ind = [idx for idx, t_ind in enumerate(time) if (t_ind>=0.2 and t_ind<=0.3)]\n",
    "erp_ind = np.array(erp_ind)\n",
    "erp_len = len(erp_ind)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loading data\n",
    "data_dict = {}\n",
    "data_dir = 'data/'\n",
    "list_of_files = [f for f in os.listdir(data_dir) if os.path.isfile(os.path.join(data_dir, f))]\n",
    "for curr_file in list_of_files:\n",
    "    if '_data.mat' in curr_file:\n",
    "        sub_id = int(curr_file[1:3])\n",
    "        data_dict[sub_id] = scipy.io.loadmat(os.path.join(data_dir,curr_file))\n",
    "    if curr_file == 'WA2.mat':\n",
    "        WA2 = scipy.io.loadmat(os.path.join(data_dir,curr_file))['WA2']\n",
    "        WA2 = np.delete(WA2,11,0)  # to remove 12th subject\n",
    "    if curr_file == 'outChoices.mat':\n",
    "        outChoices = scipy.io.loadmat(os.path.join(data_dir,curr_file))['out']\n",
    "# WA2 contains 14(sub) x 10 (prod) : value represents # of times the product was chosen total\n",
    "# outChoices : column 0: if prod_i was chosen? column 7: sub_id*100+prod_1 , col 8: sub_id*100+prod_2\n",
    "choices = np.zeros([total_sub+1, total_prod+1, total_prod+1])\n",
    "for idx in range(outChoices.shape[0]):\n",
    "    sub_id = int(outChoices[idx, 7]//100)\n",
    "    sub_id_2 = int(outChoices[idx, 8]//100)\n",
    "    if sub_id == 12:\n",
    "        continue\n",
    "    if sub_id > 12:\n",
    "        sub_id = sub_id - 1\n",
    "        sub_id_2 = sub_id_2 - 1\n",
    "    assert sub_id>0 and sub_id <= (total_sub+1) and sub_id == sub_id_2, \"Error 1: error decoding\"+str(sub_id)\n",
    "    prod_1 = int(outChoices[idx, 7]%100)\n",
    "    prod_2 = int(outChoices[idx, 8]%100)\n",
    "    assert prod_1 > 0 and prod_1 <= total_prod and  prod_2 > 0 and prod_2 <= total_prod, \"Error 2: error decoding \"+str(prod_2)\n",
    "    if prod_1 > prod_2 or prod_1==prod_2:\n",
    "        print \"check it baby\", prod_1, prod_2\n",
    "    if outChoices[idx, 0] == 0:\n",
    "        choices[sub_id, prod_1, prod_2] = choices[sub_id, prod_1, prod_2] + 1\n",
    "    elif outChoices[idx, 0] == 1:\n",
    "        choices[sub_id, prod_2, prod_1] = choices[sub_id, prod_2, prod_1] + 1\n",
    "        \n",
    "        \n",
    "pref = np.zeros([total_sub+1, total_prod+1])\n",
    "pref[1:,1:] = WA2\n",
    "#data_dict[sub_id]['sig'], ['lab'] contains original eeg signals (25x768x500) and labels (1x500)\n",
    "# What are channels?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "ncomp = 1\n",
    "sig_data = np.zeros([total_sub+1, total_prod+1, 768, 7])\n",
    "ica_data = np.zeros([total_sub+1, total_prod+1, 768, 1])\n",
    "n200_val = np.zeros([total_sub+1, total_prod+1])\n",
    "ersp_val = np.zeros([total_sub+1, total_prod+1])\n",
    "feat_val = np.zeros([total_sub+1, total_prod+1])\n",
    "for sub_id in range(1,total_sub+1):\n",
    "    sig = data_dict[sub_id]['sig']\n",
    "    lab = data_dict[sub_id]['lab'].flatten()\n",
    "    for prod_id in range(1,total_prod+1):\n",
    "        sig_prod =  sig[selected_chan, :, np.argwhere(lab==prod_id)]     \n",
    "        avg_sig_chan = np.transpose(np.mean(sig_prod, axis=0))\n",
    "        # compute ICA\n",
    "        ica_sig = FastICA(n_components=ncomp)\n",
    "        S_ = ica_sig.fit_transform(avg_sig_chan)  # Get the estimated sources\n",
    "        A_sig = ica_sig.mixing_  # Get estimated mixing matrix\n",
    "        A_sig_norm = np.linalg.norm(A_sig,axis=0)\n",
    "        S_ = S_*A_sig_norm\n",
    "    \n",
    "        if sum(A_sig)<0:\n",
    "            S_ = -1*S_\n",
    "            \n",
    "        \n",
    "        \n",
    "        n200val =  np.mean(S_[n200_ind])\n",
    "        featval = min(S_[erp_ind])\n",
    "        \n",
    "        info = mne.create_info(ch_names=['ica'], sfreq=freq, ch_types=['eeg'])\n",
    "        raw = mne.io.RawArray(np.transpose(S_[256:]), info)\n",
    "        psds, freqs = psd_multitaper(raw, low_bias=True, tmin=0.1, tmax=0.5, fmin=13, fmax=26, proj=True, n_jobs=1)\n",
    "        raw2 = mne.io.RawArray(np.transpose(S_), info)\n",
    "        psds2, freqs2 = psd_multitaper(raw2, low_bias=True, tmin=0.5, tmax=1.0, fmin=13, fmax=26, proj=True, n_jobs=1)\n",
    "        erspval = np.mean(10*np.log10(psds)) - np.mean(10*np.log10(psds2))\n",
    "\n",
    "        n200_val[sub_id, prod_id] = 10*np.log10(1 + abs(n200val)**2)\n",
    "        ersp_val[sub_id, prod_id] = erspval\n",
    "        feat_val[sub_id, prod_id] = 10*np.log10(1 + featval**2)\n",
    "        \n",
    "        sig_data[sub_id, prod_id, :, :] = avg_sig_chan\n",
    "        ica_data[sub_id, prod_id, :, :] = S_\n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pairwise Accuracy [Max - Mean - std]\n",
      "Validation \t  0.2531710586881039 0.21766517997334361 0.016270837736846972\n",
      "Train \t  0.19017083810008625 0.15890894110548684 0.014288618519160958\n",
      "Test \t  0.06923970322558004 -0.02751868787131508 0.03926872634262003\n",
      "Ranking Accuracy Mean [Train - Test - Random]\n",
      "Dist \t  [0.46996599 0.67698413 0.88809524]\n",
      "Tau \t  [0.51007937 0.27103175 0.00595238]\n",
      "NDCG \t  [0.95695327 0.92002218 0.87499508]\n",
      "Ranking Accuracy Max [Train - Test ]\n",
      "Dist \t  [0.36462585 0.33333333]\n",
      "Tau \t  [0.62585034 0.66666667]\n",
      "NDCG \t  [0.97544729 0.97425431]\n",
      "Ranking Accuracy Std [Train - Test]\n",
      "Dist \t  [0.0393621  0.11796065]\n",
      "Tau \t  [0.04224192 0.1390764 ]\n",
      "NDCG \t  [0.00860341 0.03838353]\n"
     ]
    }
   ],
   "source": [
    "# Ranking Evaluation\n",
    "\n",
    "nsim = 1\n",
    "\n",
    "train_split = 7\n",
    "test_split = 10 - train_split\n",
    "\n",
    "train_split_cv = 6\n",
    "test_split_cv = train_split - train_split_cv\n",
    "\n",
    "total_splits = int(scipy.special.comb(10, train_split))\n",
    "\n",
    "total_samples_train = int(scipy.special.comb(train_split, 2))*2\n",
    "total_samples_test = int(scipy.special.comb(test_split, 2))\n",
    "\n",
    "total_samples_train_cv = int(scipy.special.comb(train_split_cv, 2))*2\n",
    "total_samples_test_cv = int(scipy.special.comb(test_split_cv, 2))\n",
    "total_samples_train_rank_cv = int(scipy.special.comb(train_split_cv, 3))\n",
    "\n",
    "total_samples_train_rank = int(scipy.special.comb(train_split, 3))\n",
    "total_samples_test_rank = int(scipy.special.comb(test_split, 3))\n",
    "\n",
    "X_train, y_train = np.zeros([total_samples_train, 3]), np.zeros([total_samples_train,])\n",
    "X_train_rank, y_train_rank, y_train_rel = np.zeros([total_samples_train_rank, 9]), np.zeros([total_samples_train_rank,3]), np.zeros([total_samples_train_rank,3])\n",
    "\n",
    "X_test, y_test = np.zeros([total_samples_test, 3]), np.zeros([total_samples_test,])\n",
    "X_test_rank, y_test_rank, y_test_rel = np.zeros([total_samples_test_rank, 9]), np.zeros([total_samples_test_rank, 3]), np.zeros([total_samples_test_rank, 3])\n",
    "\n",
    "rank_metrics = np.zeros([total_sub + 1, total_splits, 3, 3]) # rows: ML-train, ML-test, random col: avg dist, tau, NDCG\n",
    "rank_count = np.zeros([total_sub + 1, total_splits, 3, 3])\n",
    "\n",
    "# Function to calculate NDCG score\n",
    "def ndcg_score(rel, t_rank, p_rank):\n",
    "    dcg = sum(rel*1.0/np.log2(p_rank + 1))\n",
    "    idcg = sum(rel*1.0/np.log2(t_rank + 1))\n",
    "    return dcg*1.0/idcg\n",
    "\n",
    "def cv_train(sub_id, train_prod):\n",
    "# Step 2: Preparing Training Dataset\n",
    "    cv_train_acc = []\n",
    "    for cv_train_prod_s in itertools.combinations(train_prod, train_split_cv):\n",
    "        cv_test_prod = list(set(train_prod)-set(cv_train_prod_s))\n",
    "        cv_train_prod = list(cv_train_prod_s)\n",
    "\n",
    "        X_train_cv, y_train_cv = np.zeros([total_samples_train_cv, 3]), np.zeros([total_samples_train_cv,])\n",
    "        X_train_rank_cv, y_train_rank_cv, y_train_rel_cv = np.zeros([total_samples_train_rank_cv, 9]), np.zeros([total_samples_train_rank_cv,3]), np.zeros([total_samples_train_rank_cv,3])\n",
    "#         X_test_cv, y_test_cv = np.zeros([total_samples_test_cv, 3]), np.zeros([total_samples_test_cv,])\n",
    "\n",
    "        sample_idx_train_cv = 0\n",
    "        for i in range(train_split_cv):\n",
    "            for j in range(i+1, train_split_cv):\n",
    "                prod_i, prod_j = cv_train_prod[i], cv_train_prod[j]\n",
    "\n",
    "                y_train_cv[sample_idx_train_cv] = np.sign(pref[sub_id, prod_i] - pref[sub_id, prod_j])\n",
    "                X_train_cv[sample_idx_train_cv,:] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]]) -  np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "\n",
    "                sample_idx_train_cv = sample_idx_train_cv + 1\n",
    "                y_train_cv[sample_idx_train_cv] = np.sign(pref[sub_id, prod_j] - pref[sub_id, prod_i])\n",
    "                X_train_cv[sample_idx_train_cv,:] = np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]]) - np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]])\n",
    "\n",
    "                sample_idx_train_cv = sample_idx_train_cv + 1\n",
    "\n",
    "        idx_train_cv = range(total_samples_train_cv)\n",
    "        random.shuffle(idx_train_cv)\n",
    "        X_train_cv = X_train_cv[idx_train_cv, :]\n",
    "        y_train_cv = y_train_cv[idx_train_cv, ]\n",
    "\n",
    "        # Step 2A: Preparing Training dataset for training evaluation\n",
    "        sample_idx_train_rank_cv = 0\n",
    "        for i in range(train_split_cv):\n",
    "            for j in range(i+1, train_split_cv):\n",
    "                for k in range(j+1, train_split_cv):\n",
    "                    prod_i, prod_j, prod_k = cv_train_prod[i], cv_train_prod[j], cv_train_prod[k]\n",
    "                    if not ((pref[sub_id, prod_i] == pref[sub_id, prod_j]) or (pref[sub_id, prod_i] == pref[sub_id, prod_k]) or (pref[sub_id, prod_j] == pref[sub_id, prod_k])):\n",
    "\n",
    "                        ground_pref_cv = [pref[sub_id, prod_i], pref[sub_id, prod_j], pref[sub_id, prod_k]]   # [rank_i, rank_j, rank_k]\n",
    "#                         print ground_pref_cv\n",
    "                        ground_ordr_cv = np.argsort(ground_pref_cv)\n",
    "                        ground_rank_cv = 3 - np.argsort(ground_ordr_cv)  # rank 1 means top, rank 3 means bottom\n",
    "                        y_train_rank_cv[sample_idx_train_rank_cv,:] = ground_rank_cv\n",
    "                        y_train_rel_cv[sample_idx_train_rank_cv,:] = ground_pref_cv\n",
    "\n",
    "                        X_train_rank_cv[sample_idx_train_rank_cv,0:3] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]])\n",
    "                        X_train_rank_cv[sample_idx_train_rank_cv,3:6] = np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "                        X_train_rank_cv[sample_idx_train_rank_cv,6:9] = np.array([n200_val[sub_id, prod_k], feat_val[sub_id, prod_k], ersp_val[sub_id, prod_k]])\n",
    "                        sample_idx_train_rank_cv = sample_idx_train_rank_cv + 1\n",
    "\n",
    "\n",
    "\n",
    "        clf = ElasticNet(l1_ratio=0.05, alpha=0.1, normalize=True)\n",
    "        clf.fit(X_train_cv, y_train_cv) \n",
    "\n",
    "        cv_train_acc.append(clf.score(X_train_cv, y_train_cv))\n",
    "\n",
    "        for eval_idx_cv in range(sample_idx_train_rank_cv):\n",
    "        # Distance Metric\n",
    "\n",
    "            pred_pref_cv = [clf.predict(np.reshape(X_train_rank_cv[eval_idx_cv,0:3],[1,-1]))[0], clf.predict(np.reshape(X_train_rank_cv[eval_idx_cv,3:6],[1,-1]))[0], clf.predict(np.reshape(X_train_rank_cv[eval_idx_cv,6:9],[1,-1]))[0]]\n",
    "            ord_pref_cv = np.argsort(pred_pref_cv)\n",
    "            pred_rank_cv = 3 - np.argsort(ord_pref_cv)\n",
    "            true_rank_cv = y_train_rank_cv[eval_idx_cv,]\n",
    "            dist_ml_train_cv = np.mean(abs(true_rank_cv - pred_rank_cv))  \n",
    "\n",
    "            true_rel_cv = y_train_rel_cv[eval_idx_cv,]\n",
    "\n",
    "    return cv_train_acc\n",
    "\n",
    "\n",
    "\n",
    "acc_metric = np.zeros([total_sub + 1, total_splits, 3])\n",
    "acc_count = np.zeros([total_sub + 1, total_splits, 3])\n",
    "\n",
    "test_comb = []\n",
    "\n",
    "\n",
    "\n",
    "for sub_id in range(1,total_sub+1):\n",
    "    \n",
    "    prod_cnt = pref[sub_id,1:]\n",
    "    prod_ord = np.flip(np.argsort(prod_cnt))\n",
    "    sort_idx = 10 -  np.argsort(prod_ord)\n",
    "    \n",
    "    for sim_idx in range(nsim):\n",
    "\n",
    "        acc_idx = 0\n",
    "\n",
    "        # Step 1: Choosing products for training/testing : All possible permutations/combinations\n",
    "        for subset in itertools.combinations(range(10), train_split):\n",
    "            train_prod = 1 + prod_ord[list(subset)]\n",
    "            \n",
    "            test_prod = list(set(1+prod_ord)-set(train_prod))\n",
    "            \n",
    "            subset_test = 1 + np.sort(list(set(range(10))-set(subset)))\n",
    "   \n",
    "\n",
    "            if sub_id==2:\n",
    "                test_comb.append(subset_test)\n",
    "\n",
    "            cv_train_acc = cv_train(sub_id, train_prod)#0\n",
    "            \n",
    "            \n",
    "            \n",
    "            sample_idx_train = 0\n",
    "            for i in range(train_split):\n",
    "                for j in range(i+1, train_split):\n",
    "                    prod_i, prod_j = train_prod[i], train_prod[j]\n",
    "                    if not pref[sub_id, prod_i] == pref[sub_id, prod_j]:\n",
    "\n",
    "                        y_train[sample_idx_train] = np.sign(pref[sub_id, prod_i] - pref[sub_id, prod_j])\n",
    "                        X_train[sample_idx_train,:] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]]) -  np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "\n",
    "                        sample_idx_train = sample_idx_train + 1\n",
    "                        y_train[sample_idx_train] = np.sign(pref[sub_id, prod_j] - pref[sub_id, prod_i])\n",
    "                        X_train[sample_idx_train,:] = np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]]) - np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]])\n",
    "\n",
    "                        sample_idx_train = sample_idx_train + 1\n",
    "\n",
    "            idx_train = range(total_samples_train)\n",
    "            random.shuffle(idx_train)\n",
    "            X_train = X_train[idx_train, :]\n",
    "            y_train = y_train[idx_train, ]\n",
    "            \n",
    "            # Step 2A: Preparing Training dataset for training evaluation\n",
    "            sample_idx_train_rank = 0\n",
    "            for i in range(train_split):\n",
    "                for j in range(i+1, train_split):\n",
    "                    for k in range(j+1, train_split):\n",
    "                        prod_i, prod_j, prod_k = train_prod[i], train_prod[j], train_prod[k]\n",
    "                        if not ((pref[sub_id, prod_i] == pref[sub_id, prod_j]) or (pref[sub_id, prod_i] == pref[sub_id, prod_k]) or (pref[sub_id, prod_j] == pref[sub_id, prod_k])):\n",
    "                            \n",
    "                            ground_pref = [pref[sub_id, prod_i], pref[sub_id, prod_j], pref[sub_id, prod_k]]   # [rank_i, rank_j, rank_k]\n",
    "                            ground_ordr = np.argsort(ground_pref)\n",
    "                            ground_rank = 3 - np.argsort(ground_ordr)  # rank 1 means top, rank 3 means bottom\n",
    "                            y_train_rank[sample_idx_train_rank,:] = ground_rank\n",
    "                            y_train_rel[sample_idx_train_rank,:] = ground_pref\n",
    "                            \n",
    "                            X_train_rank[sample_idx_train_rank,0:3] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]])\n",
    "                            X_train_rank[sample_idx_train_rank,3:6] = np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "                            X_train_rank[sample_idx_train_rank,6:9] = np.array([n200_val[sub_id, prod_k], feat_val[sub_id, prod_k], ersp_val[sub_id, prod_k]])\n",
    "                            sample_idx_train_rank = sample_idx_train_rank + 1\n",
    "    \n",
    "            # Step 3: Preparing Testing Dataset for pairwise evaluation\n",
    "            sample_idx_test = 0\n",
    "            for i in range(test_split):\n",
    "                for j in range(i+1, test_split):\n",
    "                    prod_i, prod_j = test_prod[i], test_prod[j]\n",
    "                    if not pref[sub_id, prod_i] == pref[sub_id, prod_j]:\n",
    "                        y_test[sample_idx_test] = np.sign(pref[sub_id, prod_i] - pref[sub_id, prod_j])\n",
    "                        X_test[sample_idx_test,:] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]]) - np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "\n",
    "                        sample_idx_test = sample_idx_test + 1\n",
    "                        \n",
    "            # Step 3B: Preparing Testing Dataset for Ranking\n",
    "            sample_idx_test_rank = 0\n",
    "            for i in range(test_split):\n",
    "                for j in range(i+1, test_split):\n",
    "                    for k in range(j+1, test_split):\n",
    "                        prod_i, prod_j, prod_k = test_prod[i], test_prod[j], test_prod[k]\n",
    "                        if not ((pref[sub_id, prod_i] == pref[sub_id, prod_j]) or (pref[sub_id, prod_i] == pref[sub_id, prod_k]) or (pref[sub_id, prod_j] == pref[sub_id, prod_k])):\n",
    "                        \n",
    "                            ground_pref = [pref[sub_id, prod_i], pref[sub_id, prod_j], pref[sub_id, prod_k]]   # [rank_i, rank_j, rank_k]\n",
    "                            ground_ordr = np.argsort(ground_pref)\n",
    "                            ground_rank = 3 - np.argsort(ground_ordr)  # rank 1 means top, rank 3 means bottom\n",
    "                            y_test_rank[sample_idx_test_rank,:] = ground_rank\n",
    "                            y_test_rel[sample_idx_test_rank,:] = ground_pref\n",
    "                            \n",
    "                            X_test_rank[sample_idx_test_rank,0:3] = np.array([n200_val[sub_id, prod_i], feat_val[sub_id, prod_i], ersp_val[sub_id, prod_i]])\n",
    "                            X_test_rank[sample_idx_test_rank,3:6] = np.array([n200_val[sub_id, prod_j], feat_val[sub_id, prod_j], ersp_val[sub_id, prod_j]])\n",
    "                            X_test_rank[sample_idx_test_rank,6:9] = np.array([n200_val[sub_id, prod_k], feat_val[sub_id, prod_k], ersp_val[sub_id, prod_k]])\n",
    "\n",
    "                            sample_idx_test_rank = sample_idx_test_rank + 1\n",
    "\n",
    "\n",
    "            # Step 4: Computing Accuracy for Pairwise \n",
    "            clf = ElasticNet(l1_ratio=0.05, alpha=0.1, normalize=True) \n",
    "\n",
    "            clf.fit(X_train, y_train) \n",
    "\n",
    "            classification_acc = [np.mean(cv_train_acc), clf.score(X_train, y_train), clf.score(X_test, y_test)]\n",
    "            acc_metric[sub_id, acc_idx, :] = acc_metric[sub_id, acc_idx, :] + classification_acc\n",
    "            acc_count[sub_id, acc_idx, :] = acc_count[sub_id, acc_idx, :] + [1,1,1]\n",
    "        \n",
    "            # Step 5: Computing Accuracy for Ranking \n",
    "            \n",
    "            \n",
    "            coef = clf.coef_.ravel() / np.linalg.norm(clf.coef_)  \n",
    "            \n",
    "            # Step 5A: For Training Data\n",
    "            dist_ml_train, tau_ml_train, ndcg_ml_train = 0, 0, 0\n",
    "            for eval_idx in range(total_samples_train_rank):\n",
    "                # Distance Metric\n",
    "                pred_pref = [np.dot(X_train_rank[eval_idx,0:3],coef), np.dot(X_train_rank[eval_idx,3:6],coef), np.dot(X_train_rank[eval_idx,6:9],coef)]\n",
    "                ord_pref = np.argsort(pred_pref)\n",
    "                pred_rank = 3 - np.argsort(ord_pref)\n",
    "                true_rank = y_train_rank[eval_idx,]\n",
    "                dist_ml_train = dist_ml_train + np.mean(abs(true_rank - pred_rank))  \n",
    "                # Tau Metric\n",
    "                tau, _ = stats.kendalltau(pred_rank, true_rank)\n",
    "                tau_ml_train = tau_ml_train + tau\n",
    "                # NDCG Metric\n",
    "                true_rel = y_train_rel[eval_idx,]\n",
    "                ndcg_ml_train = ndcg_ml_train + ndcg_score(true_rel, true_rank, pred_rank)\n",
    "            \n",
    "\n",
    "            rank_metrics[sub_id, acc_idx, 0, 0] = rank_metrics[sub_id, acc_idx, 0, 0] + dist_ml_train*1.0/total_samples_train_rank\n",
    "            rank_metrics[sub_id, acc_idx, 0, 1] = rank_metrics[sub_id, acc_idx, 0, 1] + tau_ml_train*1.0/total_samples_train_rank\n",
    "            rank_metrics[sub_id, acc_idx, 0, 2] = rank_metrics[sub_id, acc_idx, 0, 2] + ndcg_ml_train*1.0/total_samples_train_rank\n",
    "            rank_count[sub_id, acc_idx, 0, 0] = rank_count[sub_id, acc_idx, 0, 0] + 1\n",
    "            rank_count[sub_id, acc_idx, 0, 1] = rank_count[sub_id, acc_idx, 0, 1] + 1\n",
    "            rank_count[sub_id, acc_idx, 0, 2] = rank_count[sub_id, acc_idx, 0, 2] + 1\n",
    "            \n",
    "            dist_ml_test, tau_ml_test, ndcg_ml_test = 0, 0, 0\n",
    "            dist_rand, tau_rand, ndcg_rand = 0, 0, 0\n",
    "            for eval_idx in range(total_samples_test_rank):\n",
    "                pred_pref = [np.dot(X_test_rank[eval_idx,0:3],coef), np.dot(X_test_rank[eval_idx,3:6],coef), np.dot(X_test_rank[eval_idx,6:9],coef)]   \n",
    "\n",
    "                ord_pref = np.argsort(pred_pref)\n",
    "                pred_rank = 3 - np.argsort(ord_pref)\n",
    "                true_rank = y_test_rank[eval_idx,]\n",
    "                dist_ml_test = dist_ml_test + np.mean(abs(true_rank - pred_rank))\n",
    "                tau, _ = stats.kendalltau(pred_rank, true_rank)\n",
    "                tau_ml_test = tau_ml_test + tau\n",
    "                true_rel = y_test_rel[eval_idx,]\n",
    "                ndcg_ml_test = ndcg_ml_test + ndcg_score(true_rel, true_rank, pred_rank)\n",
    "                \n",
    "                pred_rand = np.random.uniform(0,1,[1,3])\n",
    "                ordr_rand = np.argsort(pred_rand)\n",
    "                rank_rand = 3 - np.argsort(ordr_rand)\n",
    "                rank_rand = rank_rand.flatten()\n",
    "                rank_worst = np.array([true_rank[1], true_rank[0], true_rank[2]])\n",
    "                dist_rand = dist_rand + np.mean(abs(true_rank - rank_rand))\n",
    "                tau, _ = stats.kendalltau(rank_rand, true_rank)\n",
    "                tau_rand = tau_rand + tau\n",
    "                ndcg_rand = ndcg_rand + ndcg_score(true_rel, true_rank, rank_rand)\n",
    "    \n",
    "            rank_metrics[sub_id, acc_idx, 1, 0] = rank_metrics[sub_id, acc_idx, 1, 0] + dist_ml_test*1.0/total_samples_test_rank\n",
    "            rank_metrics[sub_id, acc_idx, 1, 1] = rank_metrics[sub_id, acc_idx, 1, 1] + tau_ml_test*1.0/total_samples_test_rank\n",
    "            rank_metrics[sub_id, acc_idx, 1, 2] = rank_metrics[sub_id, acc_idx, 1, 2] + ndcg_ml_test*1.0/total_samples_test_rank\n",
    "            rank_count[sub_id, acc_idx, 1, 0] = rank_count[sub_id, acc_idx, 1, 0] + 1  \n",
    "            rank_count[sub_id, acc_idx, 1, 1] = rank_count[sub_id, acc_idx, 1, 1] + 1\n",
    "            rank_count[sub_id, acc_idx, 1, 2] = rank_count[sub_id, acc_idx, 1, 2] + 1\n",
    "        \n",
    "            rank_metrics[sub_id, acc_idx, 2, 0] = rank_metrics[sub_id, acc_idx, 2, 0] + dist_rand*1.0/total_samples_test_rank\n",
    "            rank_metrics[sub_id, acc_idx, 2, 1] = rank_metrics[sub_id, acc_idx, 2, 1] + tau_rand*1.0/total_samples_test_rank\n",
    "            rank_metrics[sub_id, acc_idx, 2, 2] = rank_metrics[sub_id, acc_idx, 2, 2] + ndcg_rand*1.0/total_samples_test_rank\n",
    "            rank_count[sub_id, acc_idx, 2, 0] = rank_count[sub_id, acc_idx, 2, 0] + 1  \n",
    "            rank_count[sub_id, acc_idx, 2, 1] = rank_count[sub_id, acc_idx, 2, 1] + 1\n",
    "            rank_count[sub_id, acc_idx, 2, 2] = rank_count[sub_id, acc_idx, 2, 2] + 1\n",
    "        \n",
    "            acc_idx = acc_idx + 1\n",
    "\n",
    "\n",
    "acc_metric[1:,:,:] = acc_metric[1:,:,:]/acc_count[1:,:,:]\n",
    "rank_metrics[1:,:,:,:] = rank_metrics[1:,:,:,:]/rank_count[1:,:,:,:]\n",
    "\n",
    "\n",
    "\n",
    "# acc_metric_persub = acc_metric\n",
    "acc_metric_subavg = np.mean(acc_metric[1:,:,:],axis=0)\n",
    "# rank_metrics_persub = rank_metrics\n",
    "rank_metrics_subavg = np.mean(rank_metrics[1:,:,:,:],axis=0)\n",
    "\n",
    "print \"Pairwise Accuracy [Max - Mean - std]\"\n",
    "print \"Validation \\t \", max(acc_metric_subavg[:,0]), np.mean(acc_metric_subavg[:,0]), np.std(acc_metric_subavg[:,0])\n",
    "print \"Train \\t \", max(acc_metric_subavg[:,1]), np.mean(acc_metric_subavg[:,1]), np.std(acc_metric_subavg[:,1])\n",
    "print \"Test \\t \", max(acc_metric_subavg[:,2]), np.mean(acc_metric_subavg[:,2]), np.std(acc_metric_subavg[:,2])\n",
    "\n",
    "print \"Ranking Accuracy Mean [Train - Test - Random]\"\n",
    "print \"Dist \\t \", np.mean(rank_metrics_subavg[:,:,0],axis=0)\n",
    "print \"Tau \\t \", np.mean(rank_metrics_subavg[:,:,1],axis=0)\n",
    "print \"NDCG \\t \", np.mean(rank_metrics_subavg[:,:,2],axis=0)\n",
    "\n",
    "# For random, max and std - does not make sense\n",
    "print \"Ranking Accuracy Max [Train - Test ]\"\n",
    "print \"Dist \\t \", np.min(rank_metrics_subavg[:,0:2,0],axis=0)\n",
    "print \"Tau \\t \", np.max(rank_metrics_subavg[:,0:2,1],axis=0)\n",
    "print \"NDCG \\t \", np.max(rank_metrics_subavg[:,0:2,2],axis=0)\n",
    "\n",
    "print \"Ranking Accuracy Std [Train - Test]\"\n",
    "print \"Dist \\t \", np.std(rank_metrics_subavg[:,0:2,0],axis=0)\n",
    "print \"Tau \\t \", np.std(rank_metrics_subavg[:,0:2,1],axis=0)\n",
    "print \"NDCG \\t \", np.std(rank_metrics_subavg[:,0:2,2],axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Top 5 Subjects: Avg NDCG, std 0.9729758779950863 0.0007864680056827411\n",
      "Top 5 Subjects: Avg MD, std 0.42857142857142855 0.05216405309573009\n",
      "Top 5 Combinations: Avg NDCG, std 0.9609647540217887 0.0044628265143876475\n",
      "Top 5 Combinations: Avg MD, std 0.4777777777777775 0.020786985482077337\n"
     ]
    }
   ],
   "source": [
    "avgrank_persub = np.mean(rank_metrics[1:,:,1,:], axis=0)\n",
    "avg_ndcg_persub_top5 = np.sort(avgrank_persub[:,2])[-5:]\n",
    "print \"Top 5 Subjects: Avg NDCG, std\", np.mean(avg_ndcg_persub_top5), np.std(avg_ndcg_persub_top5)\n",
    "avg_md_persub_top5 = np.sort(avgrank_persub[:,0])[:5]\n",
    "print \"Top 5 Subjects: Avg MD, std\", np.mean(avg_md_persub_top5), np.std(avg_md_persub_top5)\n",
    "\n",
    "avgrank_percomb = np.mean(rank_metrics[1:,:,1,:], axis=1)\n",
    "avg_ndcg_percomb_top5 = np.sort(avgrank_percomb[:,2])[-5:]\n",
    "print \"Top 5 Combinations: Avg NDCG, std\", np.mean(avg_ndcg_percomb_top5), np.std(avg_ndcg_percomb_top5)\n",
    "avg_md_percomb_top5 = np.sort(avgrank_percomb[:,0])[:5]\n",
    "print \"Top 5 Combinations: Avg MD, std\", np.mean(avg_md_percomb_top5), np.std(avg_md_percomb_top5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAE5ZJREFUeJzt3X+s3fd91/Hnq9e10FZCguyWEcexJzmUIcqSnrlcQrq7hhQDU/IHU3XTTWsmUQvaBFSUokRIZXJUBTFNZWNWkdttIkDrRdYaeZDVsZJddVQ3la9pSrFNUs/d8DWDeGkyFIpw477545yrnZzYPsf2uffccz/Ph3R07vfz/Xyv3+er49f9nM/5/khVIUlqw9smXYAkae0Y+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGbJp0AYO2bNlSO3bsmHQZkjRVjh8//kdVtXVYv3UX+jt27GBpaWnSZUjSVEnyB6P0c3pHkhoyUugn2ZPkxSSnkzxyifWfSfJC7/FSktf61m1P8kySU0lOJtkxvvIlSVdj6PROkhlgP3APsAwcS3K4qk6u9KmqT/T1fwi4ve9XPAF8uqqOJnkH8P1xFS9JujqjjPR3A6er6kxVXQAOAvddof/9wBcBkvwIsKmqjgJU1etV9d3rrFmSdI1GCf2bgbN9y8u9trdIciuwE3iu13Qb8FqS30zy9SS/0PvkMLjd3iRLSZbOnz9/da9AkjSycX+ROw8cqqqLveVNwF3Aw8CPAT8MPDC4UVUdqKpOVXW2bh16xJEk6RqNEvrngFv6lrf12i5lnt7UTs8y8EJvaugN4CngjmspdBSLi/D4491nSdJbjXKc/jFgV5KddMN+HvjwYKck7wZuAhYHtr0xydaqOg98AFiVg/AXF+Huu+HCBdi8GZ59FmZnV+NfkqTpNXSk3xuhPwgcAU4BT1bViST7ktzb13UeOFh9N93tTfM8DDyb5JtAgM+N8wWsWFjoBv7Fi93nhYXV+FckabqNdEZuVT0NPD3Q9qmB5Z+/zLZHgfdcY30jm5vrjvBXRvpzc6v9L0rS9Fl3l2G4VrOz3SmdhYVu4Du1I0lvtWFCH7pBb9hL0uV57R1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDVkpNBPsifJi0lOJ3nkEus/k+SF3uOlJK8NrL8hyXKSXxlX4ZKkq7dpWIckM8B+4B5gGTiW5HBVnVzpU1Wf6Ov/EHD7wK95DPjKWCqWJF2zUUb6u4HTVXWmqi4AB4H7rtD/fuCLKwtJ3gu8C3jmegqVJF2/UUL/ZuBs3/Jyr+0tktwK7ASe6y2/DfhF4OHrK1OSNA7j/iJ3HjhUVRd7yx8Dnq6q5SttlGRvkqUkS+fPnx9zSZKkFUPn9IFzwC19y9t6bZcyD3y8b3kWuCvJx4B3AJuTvF5Vb/oyuKoOAAcAOp1OjVi7JOkqjRL6x4BdSXbSDft54MODnZK8G7gJWFxpq6qf7lv/ANAZDHxJ0toZOr1TVW8ADwJHgFPAk1V1Ism+JPf2dZ0HDlaVI3VJWqey3jK60+nU0tLSpMuQpKmS5HhVdYb184xcSWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+kMsLsLjj3efJWnabRqlU5I9wC8BM8Dnq+qfD6z/DPATvcUfAN5ZVTcm+VHgs8ANwEXg01X1G+MqfrUtLsLdd8OFC7B5Mzz7LMzOTroqSbp2Q0M/yQywH7gHWAaOJTlcVSdX+lTVJ/r6PwTc3lv8LvCzVfWtJH8eOJ7kSFW9Ns4XsVoWFrqBf/Fi93lhwdCXNN1Gmd7ZDZyuqjNVdQE4CNx3hf73A18EqKqXqupbvZ//B/AysPX6Sl47c3PdEf7MTPd5bm7SFUnS9Rlleudm4Gzf8jLwvkt1THIrsBN47hLrdgObgd+7+jInY3a2O6WzsNANfEf5kqbdSHP6V2EeOFRVF/sbk/wQ8G+Bj1TV9wc3SrIX2Auwffv2MZd0fWZnDXtJG8co0zvngFv6lrf12i5lnt7UzookNwD/EfinVfX8pTaqqgNV1amqztatUzP7I0lTZ5TQPwbsSrIzyWa6wX54sFOSdwM3AYt9bZuBLwFPVNWh8ZQsSbpWQ0O/qt4AHgSOAKeAJ6vqRJJ9Se7t6zoPHKyq6mv7EPB+4IEkL/QePzrG+iVJVyFvzujJ63Q6tbS0NOkyJGmqJDleVZ1h/TwjV5IaYuhL0jqwVpd8Gfchm5Kkq7SWl3xxpC9JE3apS76sFkNfkiZsLS/54vSOJE3YWl7yxdCXpHVgrS754vSOJDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDRkp9JPsSfJiktNJHrnE+s8keaH3eCnJa33rPpLkW73HR8ZZ/DRbq5sgS1K/oTdRSTID7AfuAZaBY0kOV9XJlT5V9Ym+/g8Bt/d+/rPAPwM6QAHHe9u+OtZXMWXW8ibIktRvlJH+buB0VZ2pqgvAQeC+K/S/H/hi7+e/CRytqu/0gv4osOd6Ct4I1vImyJLUb5TQvxk427e83Gt7iyS3AjuB565225as5U2QJanfuO+ROw8cqqqLV7NRkr3AXoDt27ePuaT1Zy1vgixJ/UYJ/XPALX3L23ptlzIPfHxg27mBbRcGN6qqA8ABgE6nUyPUNPXW6ibIktRvlOmdY8CuJDuTbKYb7IcHOyV5N3AT0H88yhHgg0luSnIT8MFemyRpAoaO9KvqjSQP0g3rGeDXqupEkn3AUlWt/AGYBw5WVfVt+50kj9H9wwGwr6q+M96XIEkaVfoyel3odDq1tLQ06TIkaaokOV5VnWH9PCNXkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CXpKiwuwuOPd5+n0bjvnCVJG9biItx9d/fe1ps3d++AN203Q3KkL0kjWljoBv7Fi93nhYVJV3T1DH1JGtHcXHeEPzPTfZ6bm3RFV8/pHUka0exsd0pnYaEb+NM2tQOGviRdldnZ6Qz7FU7vSFJDDH1JaoihL0kNMfQlqSEjhX6SPUleTHI6ySOX6fOhJCeTnEjyhb72f9FrO5Xkl5NkXMVLkq7O0KN3kswA+4F7gGXgWJLDVXWyr88u4FHgzqp6Nck7e+1/DbgTeE+v638CfhxYGOeLkKRLWVyc7sMrV8Moh2zuBk5X1RmAJAeB+4CTfX0+CuyvqlcBqurlXnsBfwrYDAR4O/C/xlO6JF3eRrhkwmoYZXrnZuBs3/Jyr63fbcBtSb6a5PkkewCqahH4HeAPe48jVXXq+suWpCvbCJdMWA3jOjlrE7ALmAO2AV9J8peBLcBf7LUBHE1yV1X9bv/GSfYCewG2b98+ppLa5MdZqWvlkgkrI/1pvGTCahgl9M8Bt/Qtb+u19VsGvlZV3wO+neQl/uSPwPNV9TpAkt8GZoE3hX5VHQAOAHQ6nbr6lyHw46zUbyNcMmE1jDK9cwzYlWRnks3APHB4oM9TdAOeJFvoTvecAf478ONJNiV5O90vcZ3eWSV+nJXebHYWHn3UwO83NPSr6g3gQeAI3cB+sqpOJNmX5N5etyPAK0lO0p3D/2RVvQIcAn4P+CbwDeAbVfVbq/A6xMa4AqCk1ZWq9TWb0ul0amlpadJlTC3n9KU2JTleVZ1h/bzK5gYz7VcAlLS6vAyDJDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQlzRxi4vw+OPdZ60ub5coaaIWF+Huu+HCBdi8GZ591lt+riZH+pImamGhG/gXL3afFxYmXdHGZuhLmqi5ue4If2am+zw3N+mKNraRQj/JniQvJjmd5JHL9PlQkpNJTiT5Ql/79iTPJDnVW79jPKVL2ghmZ7tTOo895tTOWhg6p59kBtgP3AMsA8eSHK6qk319dgGPAndW1atJ3tn3K54APl1VR5O8A/j+WF+BpKk3O2vYr5VRRvq7gdNVdaaqLgAHgfsG+nwU2F9VrwJU1csASX4E2FRVR3vtr1fVd8dWvSTpqowS+jcDZ/uWl3tt/W4Dbkvy1STPJ9nT1/5akt9M8vUkv9D75CBJmoBxfZG7CdgFzAH3A59LcmOv/S7gYeDHgB8GHhjcOMneJEtJls6fPz+mkiRJg0YJ/XPALX3L23pt/ZaBw1X1var6NvAS3T8Cy8ALvamhN4CngDsG/4GqOlBVnarqbN269Vpeh1aZJ89IG8MoJ2cdA3Yl2Uk37OeBDw/0eYruCP/Xk2yhO61zBngNuDHJ1qo6D3wAWBpX8VobnjwjbRxDR/q9EfqDwBHgFPBkVZ1Isi/Jvb1uR4BXkpwEfgf4ZFW9UlUX6U7tPJvkm0CAz63GC9Hq8eQZaeMY6TIMVfU08PRA26f6fi7gH/ceg9seBd5zfWVqklZOnlkZ6XvyjDS9vPaOhlo5eWZhoRv4Tu1I08vQ10g8eUbaGLz2jiQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9TYyXdpDWnodsaiK8tIM0GY70NRFe2kGaDENfE+F9UaXJcHpHE+GlHaTJMPQ1MV7aQVp7Tu9IUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDRgr9JHuSvJjkdJJHLtPnQ0lOJjmR5AsD625IspzkV8ZRtCTp2gy9DEOSGWA/cA+wDBxLcriqTvb12QU8CtxZVa8meefAr3kM+Mr4ypYkXYtRRvq7gdNVdaaqLgAHgfsG+nwU2F9VrwJU1csrK5K8F3gX8Mx4SpauzJuzSJc3ygXXbgbO9i0vA+8b6HMbQJKvAjPAz1fVl5O8DfhF4GeAv3H95UpX5s1ZpCsb1xe5m4BdwBxwP/C5JDcCHwOerqrlK22cZG+SpSRL58+fH1NJapE3Z1ldfoqafqOM9M8Bt/Qtb+u19VsGvlZV3wO+neQlun8EZoG7knwMeAewOcnrVfWmL4Or6gBwAKDT6dQ1vRKJP7k5y8pI35uzjI+fojaGUUb6x4BdSXYm2QzMA4cH+jxFd5RPki10p3vOVNVPV9X2qtoBPAw8MRj40jit3JzlsccMpXHzU9TGMHSkX1VvJHkQOEJ3vv7XqupEkn3AUlUd7q37YJKTwEXgk1X1ymoWLl2ON2dZHX6K2hhStb5mUzqdTi0tLU26DEmXsLjoLS7XqyTHq6ozrJ+3S5Q0Mj9FTT8vwyBtQB5lo8txpC9tMB5loytxpC9tMB5loysx9KURTcuUycpRNjMzHmWjt3J6RxrBNE2ZrJyr4FE2uhRDXxrBpaZM1nOYepSNLsfpHWkETploo3CkL41gNadMPOFJa8nQl0a0GlMm0/RdgTYGp3ekCfLwSq01Q1+aIL8r0FpzekeaIA+v1Foz9KUJ8/BKrSWndySpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JD1t2N0ZOcB/7gOn7FFuCPxlTORuT+Gc59dGXun+EmsY9uraqtwzqtu9C/XkmWRrkjfKvcP8O5j67M/TPcet5HTu9IUkMMfUlqyEYM/QOTLmCdc/8M5z66MvfPcOt2H224OX1J0uVtxJG+JOkypjb0k8wk+XqS/3CJdQ8kOZ/khd7j702ixklK8vtJvtl7/UuXWJ8kv5zkdJL/kuSOSdQ5KSPsn7kkf9z3HvrUJOqcpCQ3JjmU5L8lOZVkdmB96++hYftnXb6HpvnSyv8IOAXccJn1v1FVD65hPevRT1TV5Y4V/lvArt7jfcBne88tudL+AfjdqvrJNatm/fkl4MtV9VNJNgM/MLC+9ffQsP0D6/A9NJUj/STbgL8DfH7StUyx+4Anqut54MYkPzTporQ+JPkzwPuBXwWoqgtV9dpAt2bfQyPun3VpKkMf+JfAPwG+f4U+f7f3kfNQklvWqK71pIBnkhxPsvcS628GzvYtL/faWjFs/wDMJvlGkt9O8pfWsrh1YCdwHvj13jTq55P84ECflt9Do+wfWIfvoakL/SQ/CbxcVcev0O23gB1V9R7gKPBv1qS49eWvV9UddD+CfzzJ+ydd0DozbP/8Z7qntf8V4F8BT611gRO2CbgD+GxV3Q78H+CRyZa0royyf9ble2jqQh+4E7g3ye8DB4EPJPl3/R2q6pWq+n+9xc8D713bEievqs71nl8GvgTsHuhyDuj/BLSt19aEYfunqv53Vb3e+/lp4O1Jtqx5oZOzDCxX1dd6y4fohly/lt9DQ/fPen0PTV3oV9WjVbWtqnYA88BzVfUz/X0G5hXvpfuFbzOS/GCSP73yM/BB4L8OdDsM/GzvCIy/CvxxVf3hGpc6EaPsnyR/Lkl6P++m+3/llbWudVKq6n8CZ5P8hV7T3cDJgW7NvodG2T/r9T00zUfvvEmSfcBSVR0G/mGSe4E3gO8AD0yytgl4F/Cl3vttE/CFqvpykr8PUFX/Gnga+NvAaeC7wM9NqNZJGGX//BTwD5K8AfxfYL7aO5PxIeDf945MOQP8nO+hNxm2f9ble8gzciWpIVM3vSNJunaGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDfn/7EajSIEsm8sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Plot the MD pattern as per combinations used in test rank\n",
    "\n",
    "chosen_metric =  rank_metrics[:,:,1,0]\n",
    "res = np.zeros([120*14,2])\n",
    "\n",
    "idx = 0\n",
    "for sub_id in range(1, total_sub+1):\n",
    "    for comb_id in range(0,120):\n",
    "        prod_cnt = np.sort(pref[sub_id,1:])\n",
    "#         te_comb = test_comb[comb_id]-1\n",
    "#         tr_comb = list(set(range(10))-set(te_comb))\n",
    "#         myprod_pref = prod_cnt[te_comb]\n",
    "#         res[idx, 0] = np.mean(np.diff(myprod_pref))\n",
    "        te_comb = test_comb[comb_id]\n",
    "        tr_comb = list(set(range(1,11))-set(te_comb))\n",
    "#         print tr_comb, te_comb\n",
    "        res[idx, 0] = np.mean(tr_comb)\n",
    "        res[idx, 1] = chosen_metric[sub_id, comb_id]\n",
    "        idx = idx + 1\n",
    "\n",
    "resolution = 0.15\n",
    "val_min = int((min(res[:,0]) + 0.5)*100)\n",
    "val_max = int((max(res[:,0]) - 0.5)*100 + resolution*100)\n",
    "interval = 1.0\n",
    "\n",
    "for idx in range(val_min,val_max,int(resolution*100)):\n",
    "    range_min = idx/100.0 - 0.5\n",
    "    range_max = idx/100.0 + 0.5\n",
    "    xx = (res[:,0]<range_max) *(res[:,0]>range_min)\n",
    "    yy = np.mean(res[np.argwhere(xx),1])\n",
    "    plt.plot(idx/100.0,yy,'b.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Best and Worst testing/Training Accuracies\n",
      "0.4714285714285715 0.8809523809523808\n",
      "0.7304761904761905 0.19918367346938773\n"
     ]
    }
   ],
   "source": [
    "max_ind = 10\n",
    "\n",
    "best_acc, worst_acc = [], []\n",
    "\n",
    "best_acc_tr, worst_acc_tr = [], []\n",
    "\n",
    "for sub_idx in range(1, total_sub+1):\n",
    "    chosen_metric_test = (rank_metrics[sub_idx, :,1,0])\n",
    "#     chosen_metric_train = (rank_metrics[sub_idx, :,0,2])\n",
    "    chosen_metric_train = (rank_metrics[sub_idx, :, 0, 0])# - (acc_metric[sub_idx, :, 0])) #abs(rank_metrics[sub_idx, :, 0,2])\n",
    "#     chosen_metric_train = (rank_metrics[sub_idx, :, 0, 2]) #acc_metric[sub_idx, :, 0]\n",
    "\n",
    "    sort_ind_train = np.argsort(chosen_metric_train)\n",
    "\n",
    "#     print chosen_metric_train[sort_ind_train]\n",
    "\n",
    "    worst_ind_train = sort_ind_train[:max_ind]\n",
    "    best_ind_train = list(sort_ind_train[-max_ind:])\n",
    "\n",
    "    chosen_metric_test_val =  np.array(chosen_metric_test).flatten()\n",
    "    chosen_metric_train_val =  np.array(chosen_metric_train).flatten()\n",
    "\n",
    "    best_train_idx =  np.array(best_ind_train).flatten()\n",
    "    best_acc.append(np.mean(chosen_metric_test_val[best_train_idx]))\n",
    "    best_acc_tr.append(np.mean(chosen_metric_train_val[best_train_idx]))\n",
    "\n",
    "    worst_train_idx =  np.array(worst_ind_train).flatten()\n",
    "    worst_acc.append(np.mean(chosen_metric_test_val[worst_train_idx]))\n",
    "    worst_acc_tr.append(np.mean(chosen_metric_train_val[worst_train_idx]))\n",
    "    \n",
    "# print best_acc\n",
    "# print worst_acc\n",
    "\n",
    "print \"Best and Worst testing/Training Accuracies\"\n",
    "print np.mean(best_acc), np.mean(worst_acc)\n",
    "print np.mean(best_acc_tr), np.mean(worst_acc_tr)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}